The Effects of Fractional Time Derivatives in Porothermoelastic Materials Using Finite Element Method

Abstract

In this work, a new model for porothermoelastic waves under a fractional time derivative and two time delays is utilized to study temperature increments, stress and the displacement components of the solid and fluid phases in porothermoelastic media. The governing equations are presented under Lord–Shulman theory with thermal relaxation time. The finite element method has been adopted to solve these equations due to the complex formulations of this problem. The effects of fractional parameter and porosity in porothermoelastic media have been studied. The numerical outcomes for the temperatures, the stresses and the displacement of the fluid and the solid are presented graphically. These results will allow future studies to gain a detailed insight into non-simple porothermoelasticity with various phases.

1. Introduction

Due to various applications in relevant subjects, including the field of geophysics, growing attention is being paid to the interaction between thermoelastic solids and fluids such as porothermoelastic fields and water. The porothermoelastic fields find various applications, mainly in the studying of the impacts of the use of wastes on the disintegration of concrete asphalt. The problems of heat transfer in porous materials especially at small time scales and short heat period are pursued by several researchers. Biot [1] has improved the poroelastic theory. This theory was found to deviate from physical reality in that it predicts infinite speeds of propagations for thermal waves. Lord and Shulman [2] were the first to develop a thermoelastic model that ensures finite wave speed. Their theory is called the of thermoelastic theory with one relaxation time. They obtained their model by modifying the Fourier’s law of heating conduction. Green and Lindsay [3] derived the governing equations of the thermoelastic theory with two thermal relaxation times. They used a generalization of a known thermodynamic inequality. A good review of the subject of thermoelasticity with finite wave speeds can be found in Ignaczak and Ostoja-Starzewski [4]. Sherief [5] applied the state space method to the thermoelastic model with two delay times. El-Karamany and Ezzat [6] studied the effect of modifying Fourier’s law with relaxation time and kernel function on the equations of thermoelasticity. Sherief, Allam, and El-Hagary [7] derived the model of thermo-visco-elasticity and solved a half-space problem for these equations. Sherief and Hussein [8] studied poro-thermoelasticity and obtained equations for a model predicting the finite speed of wave propagation. Sherief and Hussein [9] obtained fundamental solutions for thermoelasticity with two relaxation times under different conditions. Fan et al. [10] discussed in depth the effects of heating and cooling on wellbore pressure and stresses. Sherief et al. [11] studied a 2D axisymmetric thermoelastic problem for an infinite medium containing cylindrical heat sources of different materials with two time delays. Abbas [12] demonstrated the natural-frequencies of poroelastic hollow cylinder material. Cheng and Schanz [13] discussed the propagation of transient waves in a one-dimensional poroelastic column. A theoretical framework for isothermal wave propagation in fluid-saturated elastic porous media was confirmed by [14,15] for cases in both the low frequency and high frequency ranges. According to [15], several authors, e.g., those of papers [16,17,18,19], have participated towards researching several problems in porothermoelastic media. El-Naggar et al. [20] have discussed the impacts of rotations, voids and initial stress on plane waves under magnetothermoelastic theory. Sur [21] studied the propagation of waves on porous asphalt on account of memory response. Hobiny [22] studied the impacts of relaxation time and porosity in porothermoelastic materials under a hybrid finite element method. Abbas and Hobiny [23] studied the thermomechanical responses of a poroelastic material with time delays. Alzahrani and Abbas [24] studied generalized thermoelastic interactions in a poroelastic medium without energy dissipation. Saeed et al. [25] applied the finite element scheme to discuss thermoelastic interaction in a poroelastic medium under the Green and Lindsay model.

Recently, fractional calculus has been successfully applied to solve real-world problems, such as heat transfer, fluid flow and viscoelasticity. Several published pieces of research have proven that mathematical models with fractional derivatives can be applied to describe the process of heat transfer in different media. Various approximations and definitions of fractional derivatives have become the main subject of many pieces of research. The fractional order of strong, normal and weak thermal conductivities under the generalized thermoelastic model was modified by Youssef [26,27]. Youssef [28] used this model to study the thermal shock problem in two-dimensions using Fourier-Laplace transform. Ezzat and Karamany [29,30,31] presented a new theory of fractional heating equations depending on a Taylor expansion of the order of time-fraction. Additionally, Sherief et al. [32] depicted a new model by using the form of the law of thermal conduction. Otherwise, several researchers [33,34,35,36,37,38,39,40,41,42,43,44,45,46] have solved many problems under generalized thermoelasticity theories.

The present work is an attempt to get a new picture of the porothermoelastic model using the fractional calculus with thermal relaxation times. The numerical results are graphically presented to show the impacts of the porosity and fractional parameter in porothermoelastic medium for all considered variables.

2. Basic Equations

Following Abbas and Hobiny [23] with Ezzat and Karamany [29,30,31], the basic equations of porothermoelastic media in the absence of physical forces and the thermal source are expressed as:

Motion equations

$$\mu u_{i,jj}+\left(\mu +\lambda \right)u_{j,ij}+Q{U}_{i,jj}-R_{12}{\mathsf{\Theta }}_{,i}^f-R_{11}{\mathsf{\Theta }}_{,i}^s={\rho }_{11}{\ddot{u}}_i+{\rho }_{12}{\ddot{U}}_i,$$

(1)

$$R{U}_{j,ij}+Q{u}_{j,ij}-R_{21}{\mathsf{\Theta }}_{,i}^s-R_{22}{\mathsf{\Theta }}_{,i}^f={\rho }_{12}{\ddot{u}}_i+{\rho }_{22}{\ddot{U}}_i$$

(2)

Heat equations

$$K^s{\mathsf{\Theta }}_{,ii}^s=\left(\frac{\partial }{\partial t}+\frac{{\tau }_{os}^{\alpha }}{\mathsf{\Gamma }\left(\alpha +1\right)}\frac{{\partial }^{1+\alpha }}{\partial t^{1+\alpha }}\right)\left(F_{12}{\mathsf{\Theta }}^f+F_{11}{\mathsf{\Theta }}^s+T_o{R}_{21}ϵ+T_o{R}_{11}{e}_{kk}\right),$$

(3)

$$K^f{\mathsf{\Theta }}_{,ii}^f=\left(\frac{\partial }{\partial t}+\frac{{\tau }_{of}^{\alpha }}{\mathsf{\Gamma }\left(\alpha +1\right)}\frac{{\partial }^{1+\alpha }}{\partial t^{1+\alpha }}\right)\left(F_{22}{\mathsf{\Theta }}^f+F_{21}{\mathsf{\Theta }}^s+T_o{R}_{22}ϵ+T_o{R}_{12}{e}_{kk}\right)$$

(4)

The constitutive equations

$$e_{ij}=\frac{1}{2}\left(u_{i,j}+u_{j,i}\right), ϵ=U_{i,i},$$

(5)

$${\sigma }_{ij}=\left(\lambda e_{kk}+Qϵ-R_{12}{\mathsf{\Theta }}^f-R_{11}{\mathsf{\Theta }}^s\right){\delta }_{ij}+2\mu e_{ij},$$

(6)

$$\sigma =Q{e}_{kk}+Rϵ-R_{22}{\mathsf{\Theta }}^f-R_{21}{\mathsf{\Theta }}^s,$$

(7)

where ${\tau }_{os},{\tau }_{of}$ are the solid and fluid time delays, $Q,R,\mu , \lambda$ are the poroelastic material coefficients, ${\rho }^f=\beta {\rho }^{f\ast }$ is the fluid density per unit volume of bulk,${\rho }_{12}$ is the dynamics coupling coefficient, ${\rho }_{22}={\rho }^f-{\rho }_{12}$ is the fluid mass coefficient, $ϵ$ are the of the fluid strain components,${\rho }^s=\left(1-\beta \right){\rho }^{s\ast }$ is the solid density per unit volume of bulk, ${\rho }^{s\ast }, {\rho }^{f\ast }$ are the solid and the fluid densities,${\rho }_{11}={\rho }^s-{\rho }_{12}$ is the solid phase mass coefficient, $\beta$ is the material porosity, $\sigma$ is the normal stress of fluid, ${\sigma }_{ij}$ is the solid surface stress, $R_{11}, R_{12}, R_{21}, R_{22}$ are the thermal and mixed coefficients, $K$ is the interface coefficient of the interphase thermal conduction, $u_i, U_i$ are the displacements of the solid and fluid phases,${\mathsf{\Theta }}^f={\mathrm{T}}^f-T_o$ is the increment of fluid temperature,$C_E^f$ are the specific thermal couplings between the fluid and solid phases, ${\mathsf{\Theta }}^s={\mathrm{T}}^s-T_o$ is the increment of solid temperature, $T_o$ is the reference temperature,$e_{ij}$ are the solid strain components, $K^s=\left(1-\beta \right)K^{s\ast }$ is the thermal conductive of the solid, $K^f=\beta K^{f\ast }$ is the fluid thermal conductivity, $K^{f\ast }, K^{s\ast }$ are the thermal conductivities of fluid and solid, $C_E^f,C_E^s$ are the fluid and solid specific heating, ${\alpha }^s, {\alpha }^f$ are the phase coefficients of the thermal expansions, ${\alpha }^{sf}, {\alpha }^{fs}$ are the thermoelastic couplings between the phases with $R_{11}={\alpha }^{fs}Q+{\alpha }^{s}P$, $R_{12}={\alpha }^{sf}P+{\alpha }^{f}Q$, $R_{21}={\alpha }^{sf}R+3{\alpha }^{s}Q$, $R_{22}={\alpha }^{f}R+3{\alpha }^{sf}Q$,$F_{11}={\rho }^s{C}_E^s$, $F_{22}={\rho }^f{C}_E^f$, $F_{12}=-\left(3{\alpha }^s{R}_{12}+{\alpha }^{fs}{R}_{22}\right)T_o$, $P=3\lambda +2\mu$ and $F_{21}=-\left({\alpha }^f{R}_{21}+3{\alpha }^{sf}{R}_{11}\right)T_o$ as in [47] and taking into consideration that:

$$\frac{{\partial }^{\alpha }g\left(\mathit{r},t\right)}{\partial t^{\alpha }}=\{\begin{array}{cc}g\left(\mathit{r},t\right)-g\left(\mathit{r},0\right),& \alpha \to 0,\\ I^{\alpha -1}\frac{\partial g\left(\mathit{r},t\right)}{\partial t},& 0<\alpha <1,\\ \frac{\partial g\left(\mathit{r},t\right)}{\partial t},& \alpha =1,\end{array},$$

(8)

where $I^{\alpha }$ is the fraction of Riemann–Liouville integral introduced as a natural generalization of the m-times repeated well-known integral in $I^{m}g\left(\mathit{r},t\right)$ written in the form of convolution type:

$$I^{m}g\left(\mathit{r},t\right)={{\displaystyle \int }}_0^t\frac{{\left(t-s\right)}^{m-1}}{\mathsf{\Gamma }\left(m\right)}g\left(\mathit{r},s\right)ds, m>0,$$

(9)

where $g\left(t\right)$ is Lebesgue’s integrable function and $\mathsf{\Gamma }\left(m\right)$ is the Gamma function. In the case $g\left(\mathit{r},t\right)$ is continuous absolutely, then,

$$\underset{m\to 1}{\lim }\frac{{\partial }^{m}g\left(\mathit{r},t\right)}{\partial t^m}=\frac{\partial g\left(\mathit{r},t\right)}{\partial t},$$

(10)

The entire spectrum of local thermal conduction described by standard thermal conduction up to ballistic thermal conduction is given in Equation (8). The different parameter values with wide ranges $0<\alpha \le 1$ cover both conductivities, $0<\alpha <1$ for low conductivity and $\alpha =1$ for normal conductivity. The one-dimensional problem is studied, so our calculations are presented in $x$-directions. This means that any function will only depend on $t$ and $x$. The displacement components of solid and fluid can be written as:

$$u_x=u\left(x,t\right), u_y=0, u_z=0, U_x=U\left(x,t\right), U_y=0, U_z=0.$$

(11)

Thus, the governing formulations can be given as

$$\left(2\mu +\lambda \right)\frac{{\partial }^{2}u}{\partial x^2}+Q\frac{{\partial }^{2}U}{\partial x^2}-R_{11}\frac{\partial {\mathsf{\Theta }}^s}{\partial x}-R_{12}\frac{\partial {\mathsf{\Theta }}^f}{\partial x}={\rho }_{11}\frac{{\partial }^{2}u}{\partial t^2}+{\rho }_{12}\frac{{\partial }^{2}U}{\partial t^2},$$

(12)

$$Q\frac{{\partial }^{2}u}{\partial x^2}+R\frac{{\partial }^{2}U}{\partial x^2}-R_{21}\frac{\partial {\mathsf{\Theta }}^s}{\partial x}-R_{22}\frac{\partial {\mathsf{\Theta }}^f}{\partial x}={\rho }_{21}\frac{{\partial }^{2}u}{\partial t^2}+{\rho }_{22}\frac{{\partial }^{2}U}{\partial t^2},$$

(13)

$$K^s \frac{{\partial }^2{\mathsf{\Theta }}^s}{\partial x^2}=\left(\frac{\partial }{\partial t}+\frac{{\tau }_{os}^{\alpha }}{\mathsf{\Gamma }\left(\alpha +1\right)}\frac{{\partial }^{1+\alpha }}{\partial t^{1+\alpha }}\right)\left(F_{11}{\mathsf{\Theta }}^s+F_{12}{\mathsf{\Theta }}^f+T_o{R}_{11}\frac{\partial u}{\partial x}+T_o{R}_{21}\frac{\partial U}{\partial x}\right),$$

(14)

$$K^f \frac{{\partial }^2{\mathsf{\Theta }}^f}{\partial x^2}=\left(\frac{\partial }{\partial t}+\frac{{\tau }_{of}^{\alpha }}{\mathsf{\Gamma }\left(\alpha +1\right)}\frac{{\partial }^{1+\alpha }}{\partial t^{1+\alpha }}\right)\left(F_{21}{\mathsf{\Theta }}^s+F_{22}{\mathsf{\Theta }}^f+T_o{R}_{12}\frac{\partial u}{\partial x}+T_o{R}_{22}\frac{\partial U}{\partial x}\right),$$

(15)

$${\sigma }_{xx}=\left(\lambda +2\mu \right)\frac{\partial u}{\partial x}+Q\frac{\partial U}{\partial x}-R_{11}{\mathsf{\Theta }}^s-R_{12}{\mathsf{\Theta }}^f,$$

(16)

$$\sigma =Q\frac{\partial u}{\partial x}+R\frac{\partial U}{\partial x}-R_{21}{\mathsf{\Theta }}^s-R_{22}{\mathsf{\Theta }}^f.$$

(17)

3. Applications

The initial conditions are defined as

$$u=0, \frac{\partial u}{\partial t}=0,U=0, \frac{\partial U}{\partial t}=0, {\mathsf{\Theta }}^f=0,{\mathsf{\Theta }}^s=0,\frac{\partial {\mathsf{\Theta }}^s}{\partial t}=0, \frac{\partial {\mathsf{\Theta }}^f}{\partial t}=0, t=0,$$

(18)

while the boundary conditions are given by

$${\mathsf{\Theta }}^f\left(0,t\right)=\beta {\mathsf{\Theta }}_{o}H\left(t\right),{\mathsf{\Theta }}^s\left(0,t\right)=\left(1-\beta \right){\mathsf{\Theta }}_{o}H\left(t\right),$$

(19)

$$U\left(0,t\right)=0, u\left(0,t\right)=0,$$

(20)

where ${\mathsf{\Theta }}_o$ is a constant and $H\left(t\right)$ is the Heaviside function. For the sake of relevance, the non-dimensional parameters can be taken by

$$\left({\mathsf{\Theta }}^{s’},{\mathsf{\Theta }}^{f’}\right)=\frac{R_{11}\left({\mathsf{\Theta }}^s,{\mathsf{\Theta }}^f\right)}{\left(\lambda +2\mu \right)}, \left(x^{\prime },u^{\prime },U^{\prime }\right)=\omega c\left(x,u,U\right),\left( t^{\prime },{\tau }_{os}^{’}, {\tau }_{of}^{’}\right)=\omega c^2\left(t,{\tau }_{os}, {\tau }_{of}\right),\left({\sigma }_{xx}^{’},{\sigma }^{\prime }\right)=\frac{\left( {\sigma }_{xx}, \sigma \right)}{\left(\lambda +2\mu \right)},$$

(21)

where $\omega =\frac{F_{11}}{K^s}$ and $c=\sqrt{\frac{\lambda +2\mu }{{\mathsf{\rho }}_{11}}}$. By using these non-dimensional parameters (21), the above formulations can be rewritten as (the dashes were neglected for convenience)

$$\frac{{\partial }^{2}u}{\partial x^2}+y_1\frac{{\partial }^{2}U}{\partial x^2}-\frac{\partial {\mathsf{\Theta }}^s}{\partial x}-y_2\frac{\partial {\mathsf{\Theta }}^f}{\partial x}=\frac{{\partial }^{2}u}{\partial t^2}+y_3\frac{{\partial }^{2}U}{\partial t^2},$$

(22)

$$y_4\frac{{\partial }^{2}U}{\partial x^2}+y_1\frac{{\partial }^{2}u}{\partial x^2}-y_5\frac{\partial {\mathsf{\Theta }}^s}{\partial x}-y_6\frac{\partial {\mathsf{\Theta }}^f}{\partial x}=y_7\frac{{\partial }^{2}u}{\partial t^2}+y_8\frac{{\partial }^{2}U}{\partial t^2},$$

(23)

$$\frac{{\partial }^2{\mathsf{\Theta }}^s}{\partial x^2}=\left(\frac{\partial }{\partial t}+\frac{{\tau }_{os}^{\alpha }}{\mathsf{\Gamma }\left(\alpha +1\right)}\frac{{\partial }^{1+\alpha }}{\partial t^{1+\alpha }}\right)\left(y_9{\mathsf{\Theta }}^s+y_{10}{\mathsf{\Theta }}^f+y_{11}\frac{\partial u}{\partial x}+y_{12}\frac{\partial U}{\partial x}\right),$$

(24)

$$\frac{{\partial }^2{\mathsf{\Theta }}^f}{\partial x^2}=\left(\frac{\partial }{\partial t}+\frac{{\tau }_{of}^{\alpha }}{\mathsf{\Gamma }\left(\alpha +1\right)}\frac{{\partial }^{1+\alpha }}{\partial t^{1+\alpha }}\right)\left(y_{13}{\mathsf{\Theta }}^s+y_{s14}{\mathsf{\Theta }}^f+y_{15}\frac{\partial u}{\partial x}+y_{16}\frac{\partial U}{\partial x}\right),$$

(25)

$${\sigma }_{xx}=\frac{\partial u}{\partial x}+y_1\frac{\partial U}{\partial x}-{\mathsf{\Theta }}^s-y_2{\mathsf{\Theta }}^f,$$

(26)

$$\sigma =y_1\frac{\partial u}{\partial x}+y_4\frac{\partial U}{\partial x}-y_5{\mathsf{\Theta }}^s-y_6{\mathsf{\Theta }}^f,$$

(27)

$${\mathsf{\Theta }}^f=0,{\mathsf{\Theta }}^s=0,U=0,u=0,\frac{\partial {\mathsf{\Theta }}^s}{\partial t}=0, \frac{\partial {\mathsf{\Theta }}^f}{\partial t}=0,\frac{\partial u}{\partial t}=0,\frac{\partial U}{\partial t}=0, t=0,$$

(28)

$${\mathsf{\Theta }}^f\left(0,t\right)=\beta {\mathsf{\Theta }}_{o}H\left(t\right),{\mathsf{\Theta }}^s\left(0,t\right)=\left(1-\beta \right){\mathsf{\Theta }}_{o}H\left(t\right), U\left(0,t\right)=0, u\left(0,t\right)=0,$$

(29)

where $y_1=\frac{Q}{\lambda +2\mu }$, $y_2=\frac{R_{12}}{R_{11}}$, $y_3=\frac{{\rho }_{12}}{{\rho }_{11}}$, $y_4=\frac{R}{\lambda +2\mu }$, $y_5=\frac{R_{21}}{R_{11}}$, $y_6=\frac{R_{22}}{R_{11}}$,$y_7=\frac{{\rho }_{21}}{{\rho }_{11}}$, $y_8=\frac{{\rho }_{22}}{{\rho }_{11}}$, $y_9=1$, $y_{10}=\frac{F_{12}}{F_{11}}$,$y_{11}=\frac{T_o{R}_{11}^2}{F_{11}\left(\lambda +2\mu \right)}$, $y_{12}=\frac{T_o{R}_{11}{R}_{21}}{F_{11}\left(\lambda +2\mu \right)}$, $y_{13}=\frac{F_{21}}{\omega K^f}$, $y_{14}=\frac{F_{22}}{\omega K^f}$, $y_{15}=\frac{T_o{R}_{11}{R}_{12}}{\omega K^f\left(\lambda +2\mu \right)}$, $y_{16}=\frac{T_o{R}_{11}{R}_{22}}{\omega K^f\left(\lambda +2\mu \right)}$.

4. Numerical Scheme

In this section, the complex formulations of waves propagations in porothermoelastic media using the finite element method (FEM) are summarized. The finite element method can be considered in two ways, the first of which is spatial coordinate discretization using the standard weak formulations procedure as in [48,49]. The weak formulations of the basic formulations are derived. The independent weighting functions, to refer to the displacement of fluid $\delta U$, the displacement of solid $\delta u$, the solid temperature $\delta {\mathsf{\Theta }}^s$ and the fluid temperature $\delta {\mathsf{\Theta }}^f$ are confined. The governing formulations are multiplied by an independent test function and after that integrated on the locative domains using the boundary conditions of the problem. Subsequently, the corresponding nodal values for the solid and the fluid displacements and the fluid and the solid temperatures are expressed as

$${\mathsf{\Theta }}^s={\displaystyle \sum }_{j=1}^n{N}_j{\mathsf{\Theta }}_j^s\left(t\right), {\mathsf{\Theta }}^f={\displaystyle \sum }_{j=1}^n{N}_j{\mathsf{\Theta }}_j^f\left(t\right), u={\displaystyle \sum }_{j=1}^n{N}_j{u}_j\left(t\right), U={\displaystyle \sum }_{j=1}^n{N}_j{U}_j\left(t\right),$$

(30)

where $N$ points to the shape function and $n$ refers to the number of nodes per element where both the shape and the test functions are identical as the parts of the standard Galerkin procedures. Therefore,

$$\delta {\mathsf{\Theta }}^s={\displaystyle \sum }_{j=1}^n{N}_j\delta {\mathsf{\Theta }}_j^s, \delta {\mathsf{\Theta }}^f={\displaystyle \sum }_{j=1}^n{N}_j\delta {\mathsf{\Theta }}_j^f, \delta u={\displaystyle \sum }_{j=1}^n{N}_j\delta u_j, \delta U={\displaystyle \sum }_{j=1}^n{N}_j\delta U_j$$

(31)

The second technique states that the fractional derivatives of time for unknown variables should be determined in the next steps using the implicit schemes. Now, the finite element weak formulation of Equations (22)–(25) can be expressed as follows

$${{\displaystyle \int }}_0^L\frac{\partial \delta u}{\partial x}\left(\frac{\partial u}{\partial x}+y_1\frac{\partial U}{\partial x}-{\mathsf{\Theta }}^s-y_2{\mathsf{\Theta }}^f\right)dx+{{\displaystyle \int }}_0^L\delta u\left(\frac{{\partial }^{2}u}{\partial t^2}+y_3\frac{{\partial }^{2}U}{\partial t^2}\right)dx=\delta u{\left(\frac{\partial u}{\partial x}+y_1\frac{\partial U}{\partial x}-{\mathsf{\Theta }}^s-y_2{\mathsf{\Theta }}^f\right)}_0^L,$$

(32)

$${{\displaystyle \int }}_0^L\frac{\partial \delta U}{\partial x}\left(y_1\frac{\partial u}{\partial x}+y_4\frac{\partial U}{\partial x}-y_5{\mathsf{\Theta }}^s-y_6{\mathsf{\Theta }}^f\right)dx+{{\displaystyle \int }}_0^L\delta U\left(y_7\frac{{\partial }^{2}u}{\partial t^2}+y_8\frac{{\partial }^{2}U}{\partial t^2}\right)dx=\delta U{\left(y_1\frac{\partial u}{\partial x}+y_4\frac{\partial U}{\partial x}-y_5{\mathsf{\Theta }}^s-y_6{\mathsf{\Theta }}^f\right)}_0^L,$$

(33)

$${{\displaystyle \int }}_0^L\frac{\partial \delta {\mathsf{\Theta }}^s}{\partial x}\frac{\partial {\mathsf{\Theta }}^s}{\partial x}dx+{{\displaystyle \int }}_0^L\delta {\mathsf{\Theta }}^s\left(\left(\frac{\partial }{\partial t}+\frac{{\tau }_{os}^{\alpha }}{\mathsf{\Gamma }\left(\alpha +1\right)}\frac{{\partial }^{1+\alpha }}{\partial t^{1+\alpha }}\right)\left(y_9{\mathsf{\Theta }}^s+y_{10}{\mathsf{\Theta }}^f+y_{11}\frac{\partial u}{\partial x}+y_{12}\frac{\partial U}{\partial x}\right)\right)dx=\delta {\mathsf{\Theta }}^s{\left(\frac{\partial {\mathsf{\Theta }}^s}{\partial x}\right)}_0^L,$$

(34)

$${{\displaystyle \int }}_0^L\frac{\partial \delta {\mathsf{\Theta }}^f}{\partial x}\frac{\partial {\mathsf{\Theta }}^f}{\partial x}dx+{{\displaystyle \int }}_0^L\delta {\mathsf{\Theta }}^f\left(\left(\frac{\partial }{\partial t}+\frac{{\tau }_{of}^{\alpha }}{\mathsf{\Gamma }\left(\alpha +1\right)}\frac{{\partial }^{1+\alpha }}{\partial t^{1+\alpha }}\right)\left(y_{13}{\mathsf{\Theta }}^s+y_{14}{\mathsf{\Theta }}^f+y_{15}\frac{\partial u}{\partial x}+y_{16}\frac{\partial U}{\partial x}\right)\right)dx=\delta {\mathsf{\Theta }}^f{\left(\frac{\partial {\mathsf{\Theta }}^f}{\partial x}\right)}_0^L,$$

(35)

5. Numerical Result and Discussion

Performing numerical calculations for sandstone saturated with waterlike materials with physical constant values of relevant elastic and thermal parameters for a porothermoelastic medium as in Singh [17,50]:

$${\rho }^{s\ast }=2600\left(kg\right)\left(m^{-3}\right), {\rho }^{f\ast }=820\left(kg\right)\left(m^{-3}\right), {\rho }_{11}=2.137\left(kg\right)\left(m^{-3}\right),$$

$$C_E^s=8790\left(J\right)\left(k{g}^{-1}\right)\left(K^{-1}\right),Q=0.99663\times {10}^{10}\left(Pa\right), \lambda =0.44363\times {10}^{10}\left(Pa\right),$$

$$K^f=125.6 \left(W\right)\left(m^{-1}\right)\left(K^{-1}\right),R=0.07435\times {10}^{10}\left(Pa\right), K^s=167.47\left(W\right)\left(m^{-1}\right)\left(°K^{-1}\right),$$

$$\mu =0.2765\times {10}^{10}\left(Pa\right), C_E^f=7953\left(J\right)\left(k{g}^{-1}\right)\left(K^{-1}\right),$$

$${\alpha }^s={\alpha }^f={\alpha }^{sf}={\alpha }^{fs}=0.1\left(K^{-1}\right).$$

The numerical techniques, exposed above, were used for the displacements, the temperatures and the stress distributions of the fluid and solid phases versus the distances $x$ within the framework of the coupled porothermoelastic model under thermal relaxation times and fractional time derivatives. Figure 1, Figure 2, Figure 3, Figure 4, Figure 5 and Figure 6 present the three curves predicted by different values of the fractional time parameter $\left(\alpha =0.1,0.5,1\right)$ considering the porosity $\left(\beta =0.25\right)$, the time $\left(t=0.4\right)$ and ${\tau }_{os}={\tau }_{of}=0.05$. In these figures, the solid line (^________) and the dashed line (-----) represent the porothermoelastic model under a fractional time derivative while the dotted line (^……) represents the Lord–Schulman model. Figure 1 displays the variations of the fluid temperature as functions of the distance $x$. It is clear that it has maximum value of $x=0$, $\left({\mathsf{\Theta }}^f=0.25\right)$ and reduces with the rising of the distance $x$ to reach zeros on $x=0.35$, which depends on the values of the fractional parameter. Figure 2 shows the variations of the solid temperature with respect to the distances $x$. It is observed that it starts with its maximum values at $x=0$, $\left({\mathsf{\Theta }}^s=0.75\right)$, which satisfies the boundary condition of the problem and then gradually decreases with the increase in the distance $x$ until zero goes beyond the wave front for the porothermoelastic model. Figure 3 displays the variation of fluid displacement with respect to the distance $x$. It is clear that it begins from the zero values which satisfy the problem boundary conditions, after which it rises up to maximum values at a specific position within easy reach of the surface. Figure 4 shows the variation of solid displacement versus the distances $x$. It is clear that it starts from the zero values which satisfy the problem boundary conditions, after which it rises until it reaches maximum value at a particular position near from the surface; after that it reduces again to close to zero. Figure 5 and Figure 6 depict the variations of fluid and solid stresses along the distances $x$. It is observed that the magnitudes of the stresses increase until they reach a maximum value at a particular location near to the surface, and after that decrease with the increase in distance to reach values of zero. As expected, the fractional time derivative has a great effect on the distributions on the variations of displacements, the temperatures and the stress distributions of the liquid and solid phases. Figure 7, Figure 8, Figure 9, Figure 10, Figure 11 and Figure 12 show the effects of porosity $\left(\beta =10,20,30,40\%\right)$ under the fractional porothermoelastic model ($\alpha =0.1)$ with ${\tau }_{os}={\tau }_{of}=0.05$ on the studied fields. As expected, $\mathsf{\beta }$ has great effects on the variation of the displacements, the temperatures and the stress distributions of the fluid and solid phases along the distances $x$.