Coupling dynamic response of saturated soil with anisotropic thermal conductivity under fractional order thermoelastic theory

In this paper, a two-dimensional (2D) thermo-hydro-mechanical dynamic (THMD) coupling analysis in the presence of a half-space medium is studied using Ezzat’s fractional order generalized theory of thermoelasticity. Using normal mode analysis (NMA), the influence of the anisotropy of the thermal conduction coefficient, fractional derivatives, and frequency on the THMD response of anisotropy, fully saturated, and poroelastic subgrade is then analyzed with a time-harmonic load including mechanical load and thermal source subjected to the surface. The general relationships among the dimensionless physical variables such as the vertical displacement, excess pore water pressure, vertical stress, and temperature distribution are graphically illustrated. The NMA method does not require the integration and inverse transformation, increases the decoupling speed, and eliminates the limitation of numerical inverse transformation. The obtained results can guide the geotechnical engineering construction according to different values of load frequency, fractional order coefficient, and anisotropy of thermal conduction coefficient. This improves the subgrade stability and enriches the theoretical studies on thermo-hydro-mechanical coupling.


Introduction
The interaction between deformation and temperature distribution can accurately describe the motion of the object while ignoring the pore water.The basic equation of the thermohydro-mechanical coupling process shows that the change of thermal state will be accompanied by the change of the displacement amount and excess pore water pressure, and vice versa.
Biot [1] proposed a unified description method for the linear coupled thermoelastic theory, which allows eliminating the shortcomings of the classical decoupling theory.It is known that the classical theories based on Fourier law have some drawbacks such as the infinite velocity of thermoelastic waves and poor description of the changes of the thermoelastic behavior at low temperature and high flux.Therefore, several non-classical models, referred to as "generalized theories of thermoelastic", have been proposed.For instance, the generalized theory of thermoelasticity was proposed by Lord and Shulman [2] (the L-S generalized theory of thermoelasticity), which is based on the generalized Fourier law.It is also considered as a relaxation time for isotropic media.Green and Lindsay [3] proposed the G-L generalized theory of thermoelasticity which modifies the constitutive relations for stress-train and entropy by introducing two thermal relaxation time parameters.Green and Naghdi [4][5][6] developed three consistent logical theories: the viz.G-N I generalized theory of thermoelasticity, G-N II generalized theory of thermoelasticity, and G-N III generalized theory of thermoelasticity.The G-N-I theory assumes that the heat waves will travel at infinite velocities, while the Type II and III theories solve this problem by assuming that they will travel at finite velocities.However, no energy dissipation is observed in the G-N II theory.The G-N III model, which is mainly a combination of the first two model types, allows the energy dissipation.These are commonly used generalized theories of thermoelastic that can accurately describe the coupling effect among the thermal wave effect, temperature, and elastic fields.They can also represent the finite velocity heat wave in the medium, which is also known as the "second sound" of heat.Other theories, that can describe the coupling effect among the thermal wave effect, temperature, and elastic fields, also exist [7][8][9].These theories have been widely used to solve many kinds of problems [10][11][12][13][14][15][16][17].
In practical engineering analysis and design, the thermo-mechanical response of poroelastic media cannot ignore the influence of liquid in the pore.Biot [18][19][20] proposed the governing equations of fluid-saturated porous elastic solids and the variational principle of thermoelastic coupling problems.In recent years, the thermo-mechanical properties of saturated porous media considering pore water have been widely studied.For instance, Booker and Savvidou [21,22] developed analytical solutions to the basic problem of heat source buried deep in saturated soil, based on the Biot's thermoelastic theory.Bai [23,24] studied lot of elastic and consolidation problems for fully fluid-saturated media that are caused by thermal impact, while considering the Biot's theory.Lu et al. [25] combined these generalized theories of thermoelastic and established the coupled governing equations with porous media for solving different isotropic saturated poroelastic problems [26,27].Guo et al. [28][29][30] studied two-dimensional thermo-hydro-elastic coupling problems for a series of isotropic, homogeneous, fully saturated, porous elastic and porous viscoelastic half-space subgrades.Shakeriaski et al. [31] studied the transient response of a porous medium affected by external traction using the G-L theory of thermoelasticity.Zhu et al. [32] studied the thermo-hydro-mechanical coupling of saturated porous deep-sea sediments under the action of mine cart vibration, based on the G-L generalized theory of thermoelasticity and Darcy's law.
The fractional differential equations have been widely used in many engineering fields.The concept of fractional calculus is a natural extension of classical mathematics, in which several definitions [33][34][35] have been proved to have certain self-similarity.This self-similarity allows it to be used in many application domains.Several definitions of fractional derivatives have been proposed and used in different fields.
With the rapid development of technology, fractional calculus has been introduced and evolved in the field of thermoelasticity.Caputo [36] derived the fractional order generalized theory of thermoelasticity by introducing the Riemann-Liouville fractional integral operator, which naturally introduces the initial conditions in the problem definition.Povstenko [37] proposed an uncoupled thermoelastic theory that can be used to describe the quasi-static state.Youssef et al. [38] proposed a fractional heat conduction equation with a fractional coefficient of 0<α�2 and established a new fractional generalized thermoelastic theory, based on the Riemann-Liouville fractional integral operator.Sherief et al. [39] modified the L-S generalized theory of thermoelasticity and combined it with fractional calculus to develop a novel generalized thermoelasticity model.Youssef and Al-Lehaibi [40] developed a mathematical model which can be used to describe the thermal projectile of semiconductor solid spheres.They then solved a thermoporoelastic problem in the context of G-N theory with fractional-order parameters.Zhang and Ma [41] studied the nonlocal fractional order strain problem which is subject to moving heat source using the fractional order strain theory.Abouelregal et al. [42] proposed a generalized thermoelastic model with two-temperature characteristics: the heat transfer equation with Caputo-Fabrizio fractional differential operator and the phase lags.Abouelregal et al. [43] analyzed the thermoelastic vibrations of a nonlocal isotropic solid medium subjected to a pulsed heat flux based on the Caputo-Fabrizio fractional derivative generalized thermoelasticity.AL-Lehaibi [44] solved a two-dimensions thermoelastic problem for an isotropic and homogeneous nanobeam in terms of the thermoelasticity with one relaxation time and fractional-order strain theory.Sherief and Hussein [45] solved a two-dimensional thermoporoelastic problem for infinitely porous cylinders at certain boundary conditions, using the fractional order generalized thermo-poroelasticity theory.Hu et al. [46] used a fractional dual-phase-lag bio-thermoelastic model to solve the thermoelastic response of a biological tissue during hyperthermia treatment by a moving laser heating.Han et al. [47] studied the thermoelastic transient response of porous microplates subjected to thermal and stress shocks at the left boundary based on the Atangana-Baleanu fractional order generalized thermoelastic theory.This was performed by considering nonlocal effects and fractional order strain combined with the thermoelasticity theory of porous materials and the dual phase-lag heat conduction model.Dutta et al. [48] proposed a generalized thermo-diffusion process in a semiinfinite nonlocal fiber-reinforced double porous thermoelastic diffusive material with voids using the fractional-order Lord-Shulman thermo-elasto-diffusion theory.
The practice has shown that natural soil is formed by the long-term deposition process of rocks, which may lead to the structure anisotropy and the large difference of thermal physical properties in different directions, especially in the horizontal and vertical directions.Yue [49] studied a solution for the thermoelastic problem in vertically inhomogeneous media.Ai et al. [50] derived the 3D thermoelastic problem of layered media around a heat source and obtained its general solution.Ai and Wang [51] studied the three-dimensional thermodynamic behavior of layered materials with anisotropic thermal diffusivity in the Cartesian coordinate system.They then obtained the time-dependent thermal and mechanical responses of the material system.Ai and Wu [52] studied, using the Laplace-Hankel transforms, the impacts of the thermal diffusivity and permeability anisotropy on the thermal consolidation of multilayer porous thermoelastic media.Wang et al. [53] studied the time-dependent response of excess pore water pressure in subsurface saturated soils under the joint action of temperature load and mechanical load, based on the Laplace-Hankel transform and exact integral method.Das et al. [54] studies a three-dimensional thermoelastic problem for an anisotropic rectangular plate.
However, few studies tackle the thermo-hydro-mechanical coupling dynamic problem for anisotropic fully fluid-saturated and poroelastic semi-infinite media with mechanical load or thermal source.Therefore, this study focuses on the thermo-hydro-mechanical coupling dynamics of anisotropic saturated subgrade using the Ezzat's fractional generalized theory of thermoelasticity.Two-dimensional distributions of the on-dimensional excess pore water pressure, vertical displacement, vertical stress, and temperature are analytically derived using normal mode analysis, based on the Caputo fractional derivative of order 0<α�1.The impacts of the load frequency, fractional order parameter, and anisotropy of the thermal conduction coefficient on different physical variables in the subgrade are then detailed.The results obtained in this study have a wide application range in the field of geotechnical engineering.

Basic equations
In this study, a half-space subgrade with anisotropy of thermal conduction coefficient is considered, where the thermal source or mechanical load is subject to the surface.The constitutive equation is given by [55]: Note that in the sequel, a comma followed by a suffix denotes a derivative with respect to the material coordinate, a superposed dot denotes a derivative with respect to time, i,j = x,z and σ ij are the components of the stress tensor, ρ is the density of the medium, and u i denotes the displacement vector components.
The strain-displacement equation is given by: The constitutive equation of the THMD model can be expressed as: Where p is the excess pore water pressure, θ is the increment of the temperature which is equal to T-T 0 , T 0 is the initial reference temperature, and T is the absolute temperature.
Based on Eqs ( 1)-( 3), the constitutive equation of the THMD model can be expressed as: The heat conduction equation for a porous subgrade with fractional derivative, which is modified by the L-S generalized thermoelastic theory, can be expressed as: where m ¼ n 0 r w c w þ ð1 À n 0 Þr s c s denotes the volumetric heat capacity in different directions, n 0 represents the porosity in different directions, ρ w and ρ s are respectively the solid grains and density of pore water, c w and c s are respectively the solid grains and heat capacity of pore water, τ denotes the thermal relaxation time, K 11 and K 33 are the heat conduction coefficients in the horizontal and vertical directions, respectively.The mass conservation equation for anisotropic media is given by:

Problem formulation
The dynamic problem of anisotropy, fully fluid-saturation, and poroelasticity is studied using the Ezzat's fractional order generalized theory of thermoelasticity.A mechanical load (or a thermal source) is considered, and it is assumed that, as z!1, the considered functions are bounded and vanished, as shown in Fig 1 .In this study, the Cartesian coordinate system (x,y, z) and displacement component coordinate system u i = (u,0,w) are used.The constitutive equation is given by: The motion equation is given by:  The following dimensionless variables are introduced to simplify the derivation process and make the derivation result universal [25,56]: ðx; z; ũ; wÞ ¼ VZ 1 ðx; z; u; wÞ ð t; tÞ where According to Eq (9), Eqs ( 5)-( 8) are rewritten in dimensionless forms (the asterisk is dropped for convenience) as follows: where:

Process of normal mode analysis
The normal mode analysis is an inexpensive technique which simulates the low and large amplitude motions.It solves the problem of discrete and truncation errors in the numerical inverse transformation and it can be divided into two parts without integral transformation and inverse transformation, which increases the decoupling speed.In addition, it can be used to characterize the macro-molecular flexibility to predict the directions of compositional changes and interpret the structural experimental data.This method is applied in various fields such as structural engineering, biological physics, molecular spectroscopy, and structural biology.In this study, the NMA weighted residual method is introduced.The solution of the considered variable can be decomposed into the following normal modes [57]: ðu; w; e; s ij ; pÞðx; z; tÞ ¼ ðu 0 ; w 0 ; e 0 ; s 0 ij ; p 0 ÞðzÞexpðot þ iaxÞ ð17Þ where frequency ω is a complex time constant, i ¼ ffi ffi ffi ffi ffi ffiffi À 1 p is the imaginary value, a is the number of waves along the x-direction, and ðu 0 ; w 0 ; e 0 ; s 0 ij ; p 0 ÞðzÞ represents the amplitudes of the field variables.Note that the value of a is related to the wave number, wavelength, and frequency.
(1) Stress condition: where q 0 is the mechanical load 1. Thermal condition: 1. Excess pore water pressure condition:

Thermal source condition
1.
(1) Stress condition: 1. Thermal condition: where Q 0 is the thermal source and ψ(x,t) is the distribution function of the mechanical load or thermal source on the x-axis, which can be written as: 1. Excess pore water pressure condition:

Numerical results and discussion
A numerical example is used to demonstrate the high efficiency of the proposed analytical method.The Riemann-Liouville integral operator is taken into account in the Ezzat's fractional generalized thermoelastic model, and the dimensionless vertical displacement, excess pore water pressure, vertical stress, and temperature changes are then calculated.To study the THMD coupling of fluid-saturated porous roadbed, most of the used parameters are taken from [23,58].They are summarized as: The other constants of the porous problem are given by: In this paper, φ = 1.0 and α are set to 1.That is, the model is degraded to the dynamic response problem of isotropic media under the integer order generalized thermoelastic theory.In addition, the problem is solved using the boundary conditions presented in [30], which results are compared with the obtained dimensionless displacement, excess pore water pressure, vertical stress, and temperature under thermal shock.It can be seen from Figs 2-4 that, for φ = 1.0 and α = 1.0, the results of the two studies are mainly consistent, and the different numerical calculations cause small errors, with maximum and average deviations of 5.07% and 2.5%, respectively.The correctness and rationality of the calculated results are further verified.

Effect of the anisotropy of the thermal conduction coefficient and fractional order parameter with mechanic load or thermal source
This section studies the impact of the anisotropy of heat conductivity and fractional order parameters on the physical variables under mechanical load and thermal sources.Figs 2-5 show the variation rules of dimensionless physical variables including the vertical displacement w 0 = w/q, excess pore water pressure p 0 = p/q, vertical stress σ 0 = σ/q, and temperature θ 0 = θ/q in the subgrade, when its surface is affected by mechanic load.Figs 6-9 show the variations of the dimensionless vertical displacement w@ = w/Q 0 , dimensionless excess pore water pressure p@ = p/Q 0 , dimensionless vertical stress σ@ = σ/Q 0 , and dimensionless temperature θ@ = θ/Q 0 , in the case of a thermal source on the surface.Three different values of the thermal conduction coefficient anisotropy (φ = 0.8, φ = 0.9, and φ = 1.0) and two different values of the fractional order parameter (α = 0.1 and α = 1.0), are specified.These computations are performed for x = 1.0 at time t = 0.5.
Figs 6-9 show the influence of the variation of φ and α on different dimensionless physical variables with mechanical load.The change of φ has a certain influence on all the variables, except the dimensionless excess pore water pressure.In all the cases, φ represents the ratio of   the vertical and horizontal thermal conduction coefficient, while φ = 1.0 denotes the isotropy of the heat conduction coefficient, which is consistent with the results presented in [52].These two sets of figures verify the validity of the proposed model.
In Figs 6 and 8, the dimensionless vertical displacement (w) and vertical stress (σ) are the main differences on the ground surface.When φ increases, the curve of these dimensionless variables quickly decreases.This is mainly due to the near-surface of the vertical displacement.In addition, the vertical stress also increases with the increase of φ.However, at the same depth, all the basic disturbances of the curves disappear.It can be seen from the vertical displacement curve that, when φ increases, the difference between the two adjacent curves gradually increases, while the difference between the vertical stress curves remains the same.In addition, in these two figures, the distribution law and disturbance depth of the curves are similar.This is due to the fact that the dimensionless vertical displacement is mainly caused by dimensionless vertical stress under mechanical load, and the influence of the dimensionless pore water pressure is weak.The dimensionless temperature generated by the mechanical load and the dimensionless temperature are very small.It can be observed from Fig 7 that the change of φ and α has no clear effect on p. Since the boundary conditions assume that the excess pore water pressure is continuous, it is zero at the boundary.When the depth increases, the excess pore water pressure first increases and then decreases under the action of mechanical load.
It can be seen from Fig 9 that φ and α have obvious effects on the dimensionless temperature.When α is small, the curve of φ = 1.0 is larger than the other two curves, the curves of φ = 0.8 and φ = 0.9 are close, while that of φ = 0.8 is slightly larger.When α = 1, the curves of φ = 0.9 and φ = 1.0 are close, and the dimensionless temperature curve gradually increases with the increase of φ.When α increases, the influence of the change of φ on the dimensionless temperature becomes more regular.It can be clearly seen that the influence of the change of φ on all the physical variables in the subgrade only appears at a certain depth.This further illustrates the importance of the fractional order parameter presented in this study.show the impacts of φ and α on four different dimensionless physical variables, after the subgrade surface is subjected to thermal load from the vertical direction.It can be seen that all the dimensionless variables gradually increase with the increase of φ and α.
It can be observed from Fig 10 that the impact of φ on the dimensionless vertical displacement mainly appears at the peak of the curve in the expansion region.The dimensionless vertical displacement generated by the thermal source starts from the compression and moves into the expansion region.
It can be observed from Fig 11 that, when α is small, the impact of φ on the dimensionless excess pore water pressure appears at the main peak of the curve.More precisely, the surface of the subgrade and the surrounding area are completely overlapped.When α increases, the difference between the three curves gradually becomes clear after the curve reaches its peak.The upper surface of the subgrade is still completely overlapped, which is caused by the drainage of the upper surface.
It can be seen from Fig 12 that the variation of φ makes the difference between the dimensionless vertical stress curves more obvious, which is clearer at the peak of the curves.
The dimensionless temperature mainly overlaps on the subgrade surface, as shown in Fig 13 .When the subgrade depth slowly increases, the higher the values of φ and α, the lower the attenuation rate of the curve, and the difference becomes the most obvious within a certain depth interval of the subgrade.The dimensionless temperature gradually increases with the increase of φ, which is consistent with the variation law of presented in [52].This further proves that the developed model is reasonable and the calculated results are correct.In summary, the influence of the anisotropic parameters of the thermal conductivity coefficient on all the physical variables appears after the subgrade reaches a certain depth.This has a non-negligible impact on the safety and stability of the buried structures, pipelines, and structures.and 16 intersect at a point.Before this point, the response curves of the two different physical variables increase with the increase of ω.After this intersection point, when ω increases, the response curves of the two different physical variables decrease.This is because the larger the value of ω, the higher the decay rate of the physical variables.It can be observed from Fig 15 that the increase of ω significantly increases the excess pore water pressure, which is mainly due to the fact that the pore water has not been discharged in time after the increase of ω.It can be seen from Fig 17 that the value of the dimensionless temperature always starts at zero and then decreases along the z-axis.This is due to the fact that the surface is subjected to mechanic load, and the mechanical wave has finite speed propagating in the media.When ω increases, the curve of θ also increases, and the peak value is the most obvious, which basically appears at the same subgrade depth.When ω is constant, θ gradually increases with the decrease of α, and the difference between the two adjacent curves clearly increases.When ω increases, the peak value of the curve with a larger ω is larger than that with a smaller one.In addition, the six curves intersect at one point in the attenuation stage, and the curve with larger ω attenuates faster.It can be seen that the curves with large load frequency increase and decay faster.The influence of α on the curve is mainly at the peak.

Effect of the frequency and fractional order parameter with mechanic load or thermal source
The six curves in Figs 18-21 show the relationship between the four different dimensionless physical variables and the change of the subgrade depth when the upper surface is affected by the thermal source.
It can be seen from Fig 18 that the dimensionless vertical displacement increases (from a negative value to a positive one), reaches a peak value, and then gradually decreases until it becomes null.The six curves successively enter the expansion state within a depth range of 0.5�z�1.When ω is small, the difference between the vertical displacements caused by the change of α is not obvious.When ω increases, the difference between the values of α becomes more obvious, which is mainly reflected in the upper surface and the peak of the curve.When α increases, the dimensionless vertical displacements gradually increase.When ω increases, the difference caused by the change of α becomes more obvious.
It can be seen from Figs 19 and 20 that ω and α significantly affect the dimensionless excess pore water pressure and vertical stress.When ω is constant, the increase of α causes the peak values of the curves of the two physical variables to gradually move deeper into the subgrade.When ω is small, the two curves of α = 0.5 and α = 1.0 are very close.When ω increases, the difference between the two curves with α = 0.5 and α = 1.0 significantly increases.In addition, when α increases, the difference between the two adjacent curves gradually decreases.Moreover, the influence of α on the vertical stress and excess pore water pressure is related to the subgrade depth, to a certain extent.In the case of a constant ω, the difference between the three curves is not obvious when they are close to the upper subgrade surface.When the depth reaches z = 0.5, the difference is very obvious.In addition, the curve attenuates to zero when the depth is z = 4.5.
It can be seen from Fig 21 that α has a certain influence on the dimensionless temperature.More precisely, when w is constant, the temperature curves at the surface with different α values are close to each other.When the subgrade depth increases, the variation of α with the change of the dimensionless temperature becomes more obvious.In addition, when the subgrade depth increases, all the curves monotonically decrease and reach zero at a depth of z = 2.When ω is constant, the attenuation rate of the curve with small α value is slightly faster than that with large α value.The higher the value of ω, the more obvious the influence of the fractional coefficients on the dimensionless temperature.

Conclusion
This paper studies the impact of different variables on the dimensionless vertical displacement, excess pore water pressure, vertical stress, and temperature of the 2D anisotropic fully fluidsaturated semi-infinite subgrade under mechanical load or thermal source, based on the Ezzat's fractional order generalized theory of thermoelasticity, which considers the Riemann-Liouville integral operator and Darcy's law.The impacts of the load frequency, fractional order parameter, and anisotropy of thermal conduction coefficient on different physical variables in the subgrade are detailed.A normal mode analysis method is used to solve this THMD coupling problem.The conclusions drawn from the obtained results can be summarized as follows: 1.The normal mode analysis method is successfully applied to the solution of the eighthorder characteristic equation of anisotropic media.It can accurately analyze all the considered distributions of physical variables.
2. The frequency plays a crucial role for all the dimensionless physical variables regardless of the load type.In particular, it affects the distribution of various physical variables when considering the influence of thermal loads.
3. The fractional order coefficient has only a clear effect on the dimensionless temperature when the mechanical load is considered on the upper surface of the subgrade.However, it has obvious effect on the dimensionless excess pore water pressure, vertical stress, and temperature when the upper surface of the subgrade is subjected to thermal load.The variation of the load frequency will affect the fractional coefficient on physical variables.
4. The variation of the anisotropy of thermal conduction coefficient significantly affects all the considered dimensionless variables.When the anisotropy of thermal conduction coefficient increases, all the physical variables gradually increase.The existence of fractional order parameters makes this effect more obvious, especially in the dimensionless temperature curve under mechanical load, the dimensionless excess pore water pressure curve, and the vertical stress curve under thermal source.

Figs 14 -
show the distributions of the dimensionless vertical displacement, excess pore water pressure, vertical stress, and temperature at frequencies of ω = 1.6 and ω = 2.5, under the THMD coupling model.Different α values (α = 0.1, α = 0.5, and α = 1.0) are also considered.Similar to section 6.1, these are selected at x = 1.0 and observed during the action time t = 0.5.Note that the punctuation is omitted from the Figures.show the distribution of four dimensionless physical variables increasing along the subgrade depth through the change of ω and α.The change of the fractional order parameters in Figs 14-16 has no obvious influence on these three dimensionless physical variables.The vertical displacement and vertical stress curves of the two different frequencies in Figs 14 and 16 intersect at a point.Before this point, the response curves of the two different physical variables increase with the increase of ω.After this intersection point, when ω increases, the response curves of the two different physical variables decrease.This is because the larger the value of ω, the higher the decay rate of the physical variables.It can be observed from Fig15that the increase of ω significantly increases the excess pore water pressure, which is mainly due to the fact that the pore water has not been discharged in time after the increase of ω.