Fractional derivative modeling for axisymmetric consolidation of multilayered cross-anisotropic viscoelastic porous media

Computers and Mathematics with Applications xxx (xxxx) xxx

Contents lists available at ScienceDirect

Computers and Mathematics with Applications journal homepage: www.elsevier.com/locate/camwa

Fractional derivative modeling for axisymmetric consolidation of multilayered cross-anisotropic viscoelastic porous media ∗

Zhi Yong Ai , Yong Zhi Zhao, Wen Jie Liu Department of Geotechnical Engineering, Key Laboratory of Geotechnical and Underground Engineering of Ministry of Education, College of Civil Engineering, Tongji University, Shanghai, China

article

info

Article history: Received 3 November 2018 Received in revised form 27 August 2019 Accepted 31 August 2019 Available online xxxx Keywords: Cross-anisotropic Multilayered porous media Viscoelasticity The fractional Merchant model Axisymmetric consolidation

a b s t r a c t This paper presents a fractional derivative modeling for axisymmetric consolidation of multilayered cross-anisotropic viscoelastic porous media. The elastic solution for multilayered porous media is acquired by the extended precise integration method. The flexibility coefficient of the fractional Merchant viscoelastic model is deduced in the Laplace transformed domain, and the viscoelastic solution of the problem is further obtained by the elastic–viscoelastic correspondence principle. The numerical results are compared with those of published literatures to verify the proposed method, and the influences of fractional derivative order, viscoelasticity of solid skeleton and layering of materials on the axisymmetric consolidation of multilayered viscoelastic porous media are investigated by several examples. © 2019 Elsevier Ltd. All rights reserved.

1. Introduction Accurate prediction of the time-dependent behavior of porous media is one of the important subjects in geotechnical and geology engineering, which is helpful for optimizing design and guaranteeing engineering safety. Therefore, it is meaningful to investigate the consolidation problems of porous media. Biot [1,2] presented a 3D consolidation theory to describe the relationship between the deformation of the solid skeleton and the dissipation of pore fluid pressure. Based on the Biot’s theory, many scholars conducted research on the consolidation of porous materials. Rajapakse & Senjuntichai [3] studied the consolidation of an elastic porous half-space with the compressibility of fluid. Singh et al. [4,5] considered anisotropic permeability and solved the plane strain and axisymmetric Biot’s consolidation problem. The above studies are all aimed at a single medium layer. However, geomaterials are layered naturally. Hence, it is more reasonable to adopt a layered model. There is also much research about the consolidation of layered poroelastic media [6–11]. What is more, geomaterials often reveal different properties in different directions due to its compression in the deposition process, which is named as anisotropy. Booker & Randolph [12] considered the anisotropy of soil when solving the consolidation problem. Chen et al. [13] presented an exact solution of a semi-infinite transversely isotropic saturated soil subjected to a uniform circular loading at the ground surface. Later, Cai and Geng [14] solved a similar problem with time-varying loadings. Ai & Cheng [15] presented the extended precise integration solution for the consolidation of transversely isotropic poroelastic layered media with high computational efficiency and accuracy. The studies mentioned above show the necessity to consider the anisotropy and the layered property of porous media. On the other hand, the intrinsic rheological properties of geomaterials have important effects on the consolidation deformation. Taylor & Merchant [16] used the Kelvin viscoelastic model to produce an early solution to the consolidation ∗ Corresponding author. E-mail address: [email protected] (Z.Y. Ai). https://doi.org/10.1016/j.camwa.2019.08.033 0898-1221/© 2019 Elsevier Ltd. All rights reserved.

Please cite this article as: Z.Y. Ai, Y.Z. Zhao and W.J. Liu, Fractional derivative modeling for axisymmetric consolidation of multilayered cross-anisotropic viscoelastic porous media, Computers and Mathematics with Applications (2019), https://doi.org/10.1016/j.camwa.2019.08.033.

2

Z.Y. Ai, Y.Z. Zhao and W.J. Liu / Computers and Mathematics with Applications xxx (xxxx) xxx

with rheological properties. Thereafter, many scholars [17–20] studied the problem of 1D rheological consolidation and the experimental investigations for 1D rheological consolidation also had many results [21–24]. Compared with the research on 1D rheological consolidation, those on multi-dimensional rheological consolidation were fewer due to its complexity. Tan [25] first solved the 3D problem for the consolidation of the clay-layer with the Maxwell model. Fu et al. [26] took advantages of the finite element method to deal with engineering problems of viscoelastic soils. But overall, the study of multi-dimensional viscoelastic consolidation is still limited because of its complexity. The above classical viscoelastic models have a widespread application and could reflect the rheological properties of soft clay on the whole. However, the creep compliance and relaxation moduli of these models cannot agree well with the experimental data in some cases, especially during the initial stages of creep or relaxation [27]. A constitutive model with fractional derivatives was proposed by Gemant [28], and some recent studies [29,30] have shown that it can agree with the tests well and simulate the viscoelasticity of materials to good extent. In conclusion, the previous studies for Biot’s consolidation with rheology or creep mainly focus on the classical viscoelastic model and 1D consolidation, and it is difficult to describe properties of a porous medium and predict timedependent behaviors during the engineering construction. It is obvious that the model for multi-dimensional Biot’s consolidation of cross-anisotropic multilayered viscoelastic porous media is more in line with the real geomaterials. What is more, the fractional derivative modeling can describe creep properties of porous media to a good degree and make it consistent with the observed experimental data. Hence, it is meaningful to study the multi-dimensional consolidation of multilayered cross-anisotropic viscoelastic porous materials with a fractional derivative model. Based on the extended precise integration solution of cross-anisotropic multilayered porous media [15] and the elastic– viscoelastic correspondence principle, this paper investigates the axisymmetric Biot’s consolidation of cross-anisotropic multilayered viscoelastic porous media with a fractional Merchant model. The comparisons between the results from the presented method and existing literatures are carried out to demonstrate the efficiency and accuracy of the proposed method, and several numerical examples are designed to illustrate the effects of fractional derivative order, viscoelasticity of a solid skeleton and stratification of materials on the time-dependent behaviors of viscoelastic porous media. 2. Precise integration solution for cross-anisotropic multilayered porous media The relevant variables in the problem of Biot’s consolidation for cross-anisotropic multilayered porous media subjected to axisymmetric vertical loads are symmetrical about the z axis and are independent of θ . Hence, the static equilibrium differential equations of this problem can be written as:

∂τrz σr − σθ ∂σr + + =0 ∂r ∂z r

(1a)

∂σz ∂τrz τrz + + =0 (1b) ∂z ∂r r in which, τrz denotes shear stress components in the r-z plane; σr , σθ , σz are normal stress components in the r, θ , z directions, respectively. The constitutive equations of a cross-anisotropic elastic medium can be written in matrix form as:

⎡

c11 ⎢c12 ′ σs = ⎣ c13 0

c12 c11 c13 0

⎤

c13 c13 c33 0

0 0⎥ ε 0⎦ s c44

(2)

2 2 2 where /c11 [ = λζ ( 1 − ζ νvh , c122 )]= λζ νh + ζ νvh , c13 = λζ νvh (1 + νh ), c33 = λ 1 − νh , c44 = Gv ; in which, λ = Ev (1 + νh ) 1 − νh − 2ζ νvh , ζ = Eh /Ev , and Gv , Ev and Eh are the vertical shear modulus, vertical Young’s modulus and horizontal Young’s modulus, respectively, νh is the Poisson’s ratios of horizontal strain for the effect of vertical stress, and νv h denotes Poisson’s ratios of horizontal strain for the effect of horizontal stress. According to the principle of effective stress, we obtain

(

)

(

)

(

σ s = σ ′s − σ

)

(3)

]T

where σ s = [σr , σθ , σz , τrz ]T is the total stress vector, σ ′s = σr′ , σθ′ , σz′ , τrz′ denotes the effective stress vector and σ = [σ , σ , σ , 0]T is the excess pore fluid pressure vector, in which σ represents the excess pore fluid pressure (pressure is taken as the positive). The seepage continuity equation can be expressed as:

[

[ ( 2 ) ] ∂e ∂σ 1 ∂ σ 1 ∂σ ∂ 2σ + nβ = kh + + k v ∂t ∂t γw ∂r2 r ∂r ∂ z2 where e =

∂ ur ∂r

+

ur r

+

∂ uz ∂z

(4)

is the bulk strain; kh and kv are the coefficients of permeability in the horizontal and vertical

directions, respectively; γw , β and n denote the unit weight of fluid, compressibility coefficient and void ratio, respectively. Please cite this article as: Z.Y. Ai, Y.Z. Zhao and W.J. Liu, Fractional derivative modeling for axisymmetric consolidation of multilayered cross-anisotropic viscoelastic porous media, Computers and Mathematics with Applications (2019), https://doi.org/10.1016/j.camwa.2019.08.033.

Z.Y. Ai, Y.Z. Zhao and W.J. Liu / Computers and Mathematics with Applications xxx (xxxx) xxx

3

The expression of the total flow from time 0 to time t can be written as: t

∫ Qz = 0

kv ∂σ dt γw ∂ z

(5)

The integral transform technique is used in this paper. The Laplace transform with respect to time t and its inverse transformation are defined as [31]:

˜ , z , s) = g(r

∞

∫

g(r , z , t)e−st dt

(6a)

0

g(r , z , t) =

a+i∞

∫

1 2π i

˜ , z , s)est ds g(r

(6b)

a−i∞

˜ , z , s) is the corresponding function of g(r , z , t) in the Laplace transform domain, and the superscript ‘‘∼’’ in which g(r means that the variate √ has been in the Laplace transform domain; s represents the Laplace transform parameter with respect to t, and i = −1. The mth order Hankel transform with respect to radius r and its corresponding inverse transformation can be expressed as [31]: ∞

∫

⌢m

g (ξ , z , t) =

g(r , z , t)Jm (ξ r)rdr

(7a)

0

g(r , z , t) =

∞

∫

⌢m

g (ξ , z , t)Jm (ξ r)ξ dξ

(7b)

0

where ⌢ g m (ξ , z , t) is the corresponding function of g(r , z , t) in the Hankel transform domain, and the superscript ‘‘⌢ ’’ means that the variate is in the Hankel transform domain; ξ represents the Hankel transform parameter with respect to r, and Jm (ξ r) denotes the mth order Bessel function. By combining the Laplace transform and Hankel transform, the Laplace–Hankel transform is obtained: g (ξ , z , s) = m

∫

∞

1 2π i

g(r , z , t)Jm (ξ r)re−st drdt

(8a)

0

0

g(r , z , t) =

∞

∫

∫

c +i∞ c −i∞

∞

∫

g (ξ , z , s)Jm (ξ r)ξ est dξ ds m

(8b)

0

where g (ξ , z , s) is the corresponding function of g(r , z , t) in the Laplace–Hankel transform domain, and the superscript ‘‘-’’ means that the variate is in the Laplace–Hankel transform domain. In the Laplace transform domain, combining Eq. (2) and the principle of effective stress, we obtain: m

∂ u˜ r 1 ∂ u˜ z = τ˜rz − ∂z c44 ∂r ( ) u˜ r 1 c13 ∂ u˜ r ∂ u˜ z = + (σ˜ z + σ˜ ) − ∂z c33 c33 ∂r r

(9) (10)

Substituting Eq. (1) into Eqs. (2)–(3) and employing the Laplace transform to them, we get:

[ ( 2 ) ( ) ] 2 c11 c33 − c13 ∂ τ˜rz c13 ∂ σ˜ ∂ u˜ r 1 ∂ u˜ r u˜ r c13 ∂ σ˜ z =− + − 2 + + −1 ∂z c33 ∂r2 r ∂r r c33 ∂ r c33 ∂r

(11)

Similarly, employing the Laplace transform to Eqs. (1b) and (5) results in the following equations:

∂ σ˜ z ∂ τ˜rz τ˜rz =− − ∂z ∂r r ∂ σ˜ sγw ˜ = Qz ∂z kv

(12) (13)

Substituting Eqs. (10) and (13) into the transformed Eq. (4) leads to:

( )( ) [ ( 2 )] ∂ u˜ r 1 ∂ ∂ Q˜ z c13 u˜ r 1 kh 1 ∂ = 1− + + σ˜ z + + nβ − + σ˜ ∂z c33 ∂r r c33 c33 sγw ∂ r 2 r ∂r

(14)

Employing the Hankel transform to Eqs. (9)–(14) with respect to r and then writing the results in matrix form, we can obtain: d W (z ) = T · W (z ) (15) dz Please cite this article as: Z.Y. Ai, Y.Z. Zhao and W.J. Liu, Fractional derivative modeling for axisymmetric consolidation of multilayered cross-anisotropic viscoelastic porous media, Computers and Mathematics with Applications (2019), https://doi.org/10.1016/j.camwa.2019.08.033.

4

Z.Y. Ai, Y.Z. Zhao and W.J. Liu / Computers and Mathematics with Applications xxx (xxxx) xxx

Fig. 1. Three blocks of the multilayered porous media.

Eq. (15) is an ordinary differential matrix equation, and W (z ) = [V , U ]T , V = τ 1rz , σ 0z , σ 0

[

]T

[

0 T

and U = ur , uz , Q z 1

0

]

represent the generalized stress and displacement vectors in the Laplace–Hankel transform domain, respectively. And the [ ] [ ] 0 b2 ξ b3 ξ A D 0 0 , superscripts ‘‘0’’ and ‘‘1’’indicate the order of Hankel transformation. Besides, T = , A = −ξ B −AT 0 0 0 ⎤ ⎡ ] [ b4 0 0 b5 ξ 2 0 0 c c −c T A is the transpose of matrix of A, D = ⎣ 0 0 0 ⎦, B = 0 b1 b1 , in which, b1 = c1 , b2 = c13 , b3 = 13c 33 , 33 33 33 0 b1 ϑ 0 0 Fz b4 =

1 , c44

b5 =

2 c11 c33 −c13

c33

, Fz =

sγw kv

,ϑ =

1 c33

+ nβ +

kh sγw

ξ 2.

The problem of Biot’s consolidation, which satisfies the relationship similar to Eq. (15), can be solved by the extended precise integration method [15]. For the convenience of analysis, the multilayered porous media subjected to an external load can be divided into three layer blocks according to the depths of the calculated point HC and applied load HF as shown in Fig. 1. When HF < HC , the generalized stress and the generalized displacement of the calculated point can be deduced as:

(

V c = − I + Gˆ 12 Qˆ 3

)−1

Rˆ 1 P , U c = Qˆ 3 V c

(16)

When HF > HC , the solutions can be written as:

(

U b = − I + Qˆ 23 Gˆ 1

)−1

Rˆ 4 P , V b = −Gˆ 1 U b

(17)

where P denotes the external load vector, I is the unit matrix, the detailed expressions of the matrices Gˆ 12 , Qˆ 3 , Rˆ 1 , Qˆ 23 , Gˆ 1 and Rˆ 4 have been written in Ref. [15]. The boundary conditions BC can be expressed as follows in the Laplace–Hankel transform domain: BC = [V a ,U d , P ]

(18) 0

in which V a = [τ 1rz , σ 0z , σ 0 ] is the generalized stress vector at the surface, and U d = [ur , uz , Q z ] is the generalized displacement vector at the base. In this paper, U d = [0, 0, 0], and V a = [0, −Fz , 0], P = [0, 0, 0] when the external load acts at the surface, or V a = [0, 0, 0], P = [0, Fz , 0]. The uniformly vertical circular load is employed, and the expression in the Laplace–Hankel transform domain of a uniformly vertical circular load with radius r and intensity q can be written qbJ (ξ b) qbJ (ξ b) as 1ξ s , hence Fz = 1ξ s . 1

0

3. Solution for viscoelastic porous media with fractional derivative In order to describe the viscoelasticity of the solid skeleton, relevant viscoelastic models should be introduced. A viscoelastic model usually consists of a spring H and a dashpot N, in which elastic modulus E and viscosity coefficient η are used to describe the characteristics of spring and dashpot, respectively. The structure of the classical Merchant model Please cite this article as: Z.Y. Ai, Y.Z. Zhao and W.J. Liu, Fractional derivative modeling for axisymmetric consolidation of multilayered cross-anisotropic viscoelastic porous media, Computers and Mathematics with Applications (2019), https://doi.org/10.1016/j.camwa.2019.08.033.

Z.Y. Ai, Y.Z. Zhao and W.J. Liu / Computers and Mathematics with Applications xxx (xxxx) xxx

5

Fig. 2. Structures of viscoelastic models.

is shown in Fig. 2(a), and its constitutive equation can be expressed as:

σs′ +

η1 dσs′ E0 E1 η 1 E 0 dε = ε+ E0 + E1 dt E0 + E1 E0 + E1 dt

(19)

The spring that represents the elastic element satisfies Hooke’s law, i.e. σs′ (t) ∼ ε (t), which could be considered that stress is proportional to the zero-order derivative of the strain versus time. And the dashpot that represents the viscous fluid satisfies the Newtonian friction law, i.e. σs′ (t) ∼ d1 ε (t)/dt1 , which means that stress is proportional to the first order derivative of the strain versus time. Based on this view, Blair [32] and Gerasimov [33] indicated a constitutive model that satisfies σs′ (t) ∼ dα ε (t)/dt α (0 < α ≤ 1), and this is a constitutive model with fractional derivative. It means that stress is proportional to the fractional derivative of the strain versus time and can reflect the developments of strain over time better. Many scholars defined different fractional derivatives [34], and the definition of Riemann–Liouville [34] is employed in this study. During the interval (0, T ), the fractional differential of function g(t) is defined as follows: Dα [g(t)] =

=[

1

1

·

d

Γ (1 − α ) dt

Γ (1 − α ) ·

tα

∫ 0

t

g(τ )

g(0) · t −α 1 dτ = + α (t − τ ) Γ (1 − α ) Γ (1 − α )

]dg = Iα (t) · dg

t

∫

g ′ (t − τ )dτ · τ −α 0

(20)

(0 ≤ t ≤ T )

where α (0 < α ≤ 1) denotes the order, Dα is the derivative differential operator; Iα (t) represents the Abel kernel; Γ is the Gamma function, which is defined as:

Γ (x) =

∫

∞

t x−1 · e−t dt

Re(x) > 0;

Γ (1 + x) = x · Γ (x)

(21)

0

The structure of the fractional Merchant model is shown in Fig. 2(b), and it is obtained by replacing the ideal viscous dashpot shown in Fig. 2(c) in the classical Merchant model with an Abel viscous dashpot shown in Fig. 2(d). The relationship between stress and strain in an ideal viscous dashpot can be expressed as:

σs′ (t) = η1

dε (t)

dt The relationship between stress and strain in an Abel viscous dashpot can be expressed as:

σs′ (t) = η1

dα ε (t)

, (0 < α ≤ 1) dt α Hence, the constitutive equations of a fractional Merchant model can be deduced as follows: σs′ +

η1 dα σs′ E0 E1 η1 E0 dα ε = ε+ α E0 + E1 dt E0 + E1 E0 + E1 dt α

(22)

(23)

(24)

The fractional derivative viscoelastic model owns the advantages of the classical viscoelastic model, i.e., parameters are clear and easy to understand. What is more, it can describe the results of creep tests better with the combination of only a few components [35,36]. Please cite this article as: Z.Y. Ai, Y.Z. Zhao and W.J. Liu, Fractional derivative modeling for axisymmetric consolidation of multilayered cross-anisotropic viscoelastic porous media, Computers and Mathematics with Applications (2019), https://doi.org/10.1016/j.camwa.2019.08.033.

6

Z.Y. Ai, Y.Z. Zhao and W.J. Liu / Computers and Mathematics with Applications xxx (xxxx) xxx

Lee [37] presented the elastic–viscoelastic correspondence principle, i.e. the basic governing equations of viscoelastic bodies in the Laplace transform domain are identical to those of elastic bodies in form. Thus, the solution of the corresponding viscoelastic problem in Laplace transform domain can be obtained by displacing elastic parameters of the elastic solutions with corresponding viscoelastic parameters. The relationship between stress and strain of isotropic elastic media is only decided by elastic modulus under the assumption that Poisson’s ratio keeps constant. Therefore, the corresponding conversion between an elastic problem and a viscoelastic one can be carried out by replacing the elastic modulus of the elastic solution with the reciprocal of flexibility coefficient V (s) = ε˜ (s) /σ˜ s′ (s) of the viscoelastic medium in the Laplace transform domain. By combining the fractional differential definition of Riemann–Liouville, and taking the Laplace transform to Eqs. (19), (24) with respect to time t, the flexibility coefficients of a classical Merchant model and a fractional Merchant model can be expressed as follows, respectively: V (s) = V (s) =

1 E0 1 E0

+ +

1

(25)

E 1 + sη 1 1

(26)

E1 + sα η1

It is necessary to state that the elastic–viscoelastic correspondence principle is only suitable for an isotropic elastic medium. The equilibrium differential equations and geometric equations of transversely isotropic media are consistent with homogeneous isotropic materials, but different in constitutive equations. In order to extend the elastic–viscoelastic correspondence principle to a cross-anisotropic medium, the assumption is employed that the ratio of horizontal and vertical elastic modulus is constant for a cross-anisotropic medium, which is also applied to the viscosity coefficient. Although the assumption is a special case of transversely isotropic media, research based on this hypothesis also has good engineering application value [38,39]. Hence, the flexibility coefficients in Eqs. (25) and (26) can be written as:

ζ · Vh (s) = ς · Vg (s) = Vv (s) = ζ · Vh (s) = ς · Vg (s) = Vv (s) =

1 E0v 1 E0v

+ +

1

(27)

E1v + sη1v 1

(28)

E 1 v + sα η 1 v

E

where ς = GE v , ζ = Eh ; E0v , E1v and η1v denote the vertical parameters of E0 , E1 and η1 , respectively. v v The solutions of viscoelastic porous media in the Laplace transform domain can be acquired by displacing elastic modulus of the precise integration solution for Biot’s consolidation of elastic porous media obtained above with reciprocals of the flexibility coefficient obtained in Eqs. (27) and (28). Then the real solutions of viscoelastic porous materials can be obtained by the numerical inverse transformation. 4. Numerical analysis 4.1. Verification Two examples are put forward to verify the feasibility of the proposed theory in this study. First, the proposed theory for viscoelastic cross-anisotropic porous media is degenerated into elastic cross-anisotropic porous media, and the solution of an elastic cross-anisotropic porous medium subjected to a uniformly vertical circular load was obtained by Booker and Randolph [12]. As depicted in Fig. 3, the parameters are as follows: Eh = 0.5Ev , Gv = 0.3Ev , νv h = νh = 0.25, kv = kh . c k Besides, the following dimensionless parameters are introduced: dimensionless time factor τ = γ11b2v t, dimensionless displacement wr∗ =

Ev uz (r) . bq

w

Fig. 3 shows the surface settlement along the radius with time. It can be seen that the results from this study fit well with those by Booker and Randolph [12], which proves the credibility of the proposed method when degenerated into an elastic state. Then the verification in a viscoelastic state is carried out here. Since the research on the consolidation of a viscoelastic porous medium mainly focuses on one-dimensional consolidation and there are few multi-dimensionalsolutions, the theory of this paper is degraded into one-dimensional viscoelastic consolidation for verification. A single porous medium layer is investigated in this example. The top of medium layer is permeable and the bottom is impermeable. The materials parameters are assumed as follows: H = 19 m, γw = 0.01 MN/m2 , kv = 1 × 10−7 m/s, η = 230 MPa · d, Ev = 12 MPa, q = 0.1 MPa. Besides, the Merchant model in this paper is degenerated to the Kelvin–Voigt constitutive model with E0 → ∞. Li et al. [18] obtained the analytical solution of one-dimensional consolidation of a viscoelastic porous medium layer. Results in Li et al. [18] are recomputed and compared with the results using the proposed method in Fig. 4. As shown in Fig. 4, the results of this study are in good agreement with the results of Li et al. [18], which validates the accuracy of the presented theory further. Please cite this article as: Z.Y. Ai, Y.Z. Zhao and W.J. Liu, Fractional derivative modeling for axisymmetric consolidation of multilayered cross-anisotropic viscoelastic porous media, Computers and Mathematics with Applications (2019), https://doi.org/10.1016/j.camwa.2019.08.033.

Z.Y. Ai, Y.Z. Zhao and W.J. Liu / Computers and Mathematics with Applications xxx (xxxx) xxx

7

Fig. 3. Surface settlements of an elastic cross-anisotropic porous medium.

Fig. 4. Pore fluid pressure curves of viscoelastic porous medium layer.

4.2. Influence of the fractional derivative order Compared with the classical viscoelastic model, the fractional order α is the key parameter in fractional derivative viscoelastic model. This example will analyze the influence of the fractional derivative order on the consolidation and time-dependent behavior of viscoelastic porous media. Four cases are adopted in the analysis, i.e., α = 0.4, 0.6, 0.8, 1.0. As shown in Fig. 5, a two-layered model is employed, and a vertical uniformly distributed circular load acts at the surface of the porous media. The material parameters of porous media are listed in Table 1. The calculated point for displacement is origin (0, 0). In the following examples, the dimensionless parameters are adopted as follows: dimensionless time factor uz (o) τ = Eγ0vbk2v t, dimensionless displacement w0∗ = E0vbq . w As shown in Fig. 5, the fractional derivative order has no influence on the initial settlement and the final settlement, while it affects the rate of consolidation heavily, and the rate of consolidation increases with α increasing. When the value of α gets smaller, a little change of α will bring great effects on the consolidation rate, which means that the creep properties of porous media are standing out. This deduction is also consistent with previous research [30], which also illustrates the correctness of this theory. Obviously, according to test results of different porous media, we can change the Please cite this article as: Z.Y. Ai, Y.Z. Zhao and W.J. Liu, Fractional derivative modeling for axisymmetric consolidation of multilayered cross-anisotropic viscoelastic porous media, Computers and Mathematics with Applications (2019), https://doi.org/10.1016/j.camwa.2019.08.033.

8

Z.Y. Ai, Y.Z. Zhao and W.J. Liu / Computers and Mathematics with Applications xxx (xxxx) xxx

Fig. 5. Central settlements with different fractional orders.

Table 1 Parameters of viscoelastic porous media with fractional order derivative. Layer

1 2

Parameters E0v /MPa

E1v /MPa

η1v /MPa · s

kz (= kh ) /cm/s

ζ

ς

νvh (= νh )

h/m

25 10

5 4

5E7 2E8

1E−5 1E−6

2.0 1.5

0.40 0.37

0.25 0.35

5.0 5.0

Table 2 The parameters of viscoelastic porous media in 8 cases. Case

Parameters E0v /MPa

E1v /MPa

η1v /MPa · s

1 2 3 4 5 6 7 8

10.0 15.0 20.0 20.0 20.0 20.0 20.0 20.0

5.0 5.0 5.0 10.0 15.0 5.0 5.0

1.0E6 1.0E6 1.0E6 1.0E6 1.0E6 1.0E7 1.0E8 0

∞

Table 3 The parameters of viscoelastic porous media with fractional order derivative. Layer

1

Parameters kz (= kh ) /cm/s

ζ

ς

νvh (= νh )

H/m

α

5E−6

1.5

0.37

0.35

10

0.7

fractional derivative order to simulate the consolidation behavior of various porous media with different characteristics. The fractional Merchant model has few parameters, but it can simulate different porous media well, which leads to efficient analysis for porous media with complex characteristics. 4.3. Influence of the viscoelasticity of solid skeleton The influences of the viscoelastic parameters E0v , E1v and η1v of the Merchant model on the consolidation and timedependent behavior of viscoelastic porous media are analyzed in this section. A single-layered model shown in Fig. 6 is employed. Eight cases are designed as shown in Table 2, and other parameters adopted are listed in Table 3. Please cite this article as: Z.Y. Ai, Y.Z. Zhao and W.J. Liu, Fractional derivative modeling for axisymmetric consolidation of multilayered cross-anisotropic viscoelastic porous media, Computers and Mathematics with Applications (2019), https://doi.org/10.1016/j.camwa.2019.08.033.

Z.Y. Ai, Y.Z. Zhao and W.J. Liu / Computers and Mathematics with Applications xxx (xxxx) xxx

9

Fig. 6. Central settlements with different E0v .

Fig. 7. Excess pore fluid pressures at point d with different E0v .

4.3.1. Influence of parameter E0v The central displacements and the excess pore fluid pressures at point d (b, b) for different cases are plotted in Figs. 6 E k and 7, respectively. The dimensionless parameters are adopted as follows: dimensionless time factor τ = γ1vb2v t, and w

dimensionless displacement w0∗ = 1vbqz . For case 8, E1v is adopted as 5 MPa when calculating the dimensionless parameters, especially. It can be seen from Fig. 6 that the settlement differences of cases 1∼3 are mainly in initial settlements. The consolidation rate and later settlements are almost the same, so that the differences between the final settlements are equal to those of the initial settlements. Case 8 simulates the consolidation of porous media. The initial settlements of case 1 and case 8 with the same E0v are in good agreement. It is because E0v controls the primary consolidation and has few effects on creep and the consolidation rate. As a result, the initial settlements decrease with E0v increasing. As shown in Fig. 7, excess pore fluid pressure increases slightly at the beginning of consolidation, which is called the Mandel–Cryer effect. When τ ≈ 1, the excess pore fluid pressure is basically dissipated. And the consolidation is completed for poroelastic media as shown in Fig. 6 (case 8). But the creep deformation is very impressive for cases 1∼3, which shows that the Merchant model in this paper can simulate the creep of material better. E

u (o)

Please cite this article as: Z.Y. Ai, Y.Z. Zhao and W.J. Liu, Fractional derivative modeling for axisymmetric consolidation of multilayered cross-anisotropic viscoelastic porous media, Computers and Mathematics with Applications (2019), https://doi.org/10.1016/j.camwa.2019.08.033.

10

Z.Y. Ai, Y.Z. Zhao and W.J. Liu / Computers and Mathematics with Applications xxx (xxxx) xxx

Fig. 8. Central settlements with time with different E1v .

Fig. 9. Excess pore fluid pressures at d with time with different E1v .

4.3.2. Influence of parameter E1v The central displacements and the excess pore fluid pressures at point d (b, b) are calculated, respectively. The E k dimensionless parameters are adopted as follows: dimensionless time factor τ = γ0vb2v t, dimensionless displacement

w0∗ =

E0v uz (o) . bq

w

Fig. 8 illustrates that the final settlements decrease with E1v increasing, and E1v rarely affects the initial settlements and the rate of consolidation. Fig. 9 shows that the excess pore fluid pressure dissipates with a same pattern. It means that the time required to complete the primary consolidation are the same, which is in accordance with Fig. 8. When the excess fluid pressure’s dissipation is completed, settlements between cases (3∼5, 8) have significant differences, which indicates E1v controls the creep stage. In conclusion, different E1v leads to different final settlements and a same primary consolidation stage, respectively. 4.3.3. Influence of parameter η1v In this section, the calculated points and dimensionless parameters are the same with those in Section 4.3.2. As shown in Fig. 10, η1v has little impact on initial settlements and final settlements, while it affects the rate of consolidation obviously, and the rate of consolidation decreases with η1v increasing. This is because η1v reflects the thickness of the solid particle combined with the fluid film to some extent. The thickness increases with the increasing of η1v , which means it is more difficult for the solid particles to overcome the cohesive force and slip, and the rate of consolidation gets Please cite this article as: Z.Y. Ai, Y.Z. Zhao and W.J. Liu, Fractional derivative modeling for axisymmetric consolidation of multilayered cross-anisotropic viscoelastic porous media, Computers and Mathematics with Applications (2019), https://doi.org/10.1016/j.camwa.2019.08.033.

Z.Y. Ai, Y.Z. Zhao and W.J. Liu / Computers and Mathematics with Applications xxx (xxxx) xxx

11

Fig. 10. Central settlements with time with different η1v .

Fig. 11. Excess pore fluid pressures at point d with time with different η1v . Table 4 The parameters of viscoelastic porous media for three cases. Case

Eo1v /MPa

Eo2v /MPa

Eo3v /MPa

1 2 3

25 10 30

20 20 20

10 25 10

slower. Fig. 11 illustrates that η1v has no influence on the dissipation of excess pore fluid pressures. It is also accordance with Fig. 10, where different cases have the same primary consolidation stage. 4.4. Influence of the layering of materials Natural geomaterials are usually layered and this property will affect the engineering construction greatly. A threelayered model plotted in Fig. 12 is adopted to analyze the influence of layering on consolidation. The surface of porous media suffers a uniformly vertical circular load. We employ proper values of E0v to simulate layered materials. Three cases are employed and the parameters are listed in Table 4, in which E0i v denotes E0v of the ith layer. E1v = 5 MPa, η1v = 1.0 × 107 MPa s. Other parameters adopted are listed in Table 3. For the sake of analysis, the dimensionless time Please cite this article as: Z.Y. Ai, Y.Z. Zhao and W.J. Liu, Fractional derivative modeling for axisymmetric consolidation of multilayered cross-anisotropic viscoelastic porous media, Computers and Mathematics with Applications (2019), https://doi.org/10.1016/j.camwa.2019.08.033.

12

Z.Y. Ai, Y.Z. Zhao and W.J. Liu / Computers and Mathematics with Applications xxx (xxxx) xxx

Fig. 12. Central settlements with time.

Fig. 13. Surface displacements along radius.

factor and displacement are defined as: τ =

E1v kv γw b 2

t, w0∗ =

E1v uz (o) , bq

respectively. The central settlements with time, the

surface displacements along radius and central displacements along with depth in three cases are plotted in Figs. 12–14, respectively. It can be seen from Fig. 12 that although the average moduli of case (1) and case (2) are the same, the central settlement of case (1) is smaller than that of case (2) for the upper layer’s modulus is smaller. The average moduli of Layer 2 and Layer 3 in case (2) and case (3) are the same, but the modulus of the upper layer in case (3) is larger. Therefore, the central settlements of case (3) are smaller. Hence, the reinforcement of the upper layer of materials can reduce the surface settlements effectively. Figs. 13 and 14 illustrate the surface displacements along radius and central displacements along depth, respectively. It is noted that the displacements have a further development as time goes on, and the displacements near the load are drastically different. So, it is important to consider the settlement differences near the load in engineering design. 5. Conclusions Based on the extended precise integration method and the elastic–viscoelastic correspondence principle, we present a solution for Biot’s consolidation of cross-anisotropic multilayered viscoelastic porous media subjected to a uniformly Please cite this article as: Z.Y. Ai, Y.Z. Zhao and W.J. Liu, Fractional derivative modeling for axisymmetric consolidation of multilayered cross-anisotropic viscoelastic porous media, Computers and Mathematics with Applications (2019), https://doi.org/10.1016/j.camwa.2019.08.033.

Z.Y. Ai, Y.Z. Zhao and W.J. Liu / Computers and Mathematics with Applications xxx (xxxx) xxx

13

Fig. 14. Central displacements along depth.

vertical circular load by introducing the Merchant model with fractional derivatives. Firstly, the ordinary differential matrix equation in the transform domain for cross-anisotropic porous media is derived based on the fundamental governing equations of an axisymmetric poroelastic body and the Laplace–Hankel transform. Then the extended precise integration solution for Biot’s consolidation of a cross-anisotropic multilayered porous media is obtained. By introducing the fractional Merchant model and the elastic–viscoelastic correspondence principle, we acquire the solution for viscoelastic porous media. The proposed method is verified by the comparison with those in published literatures. What is more, the influences of fractional derivative order, viscoelasticity of solid skeleton and layering of materials on the time-dependent behavior of multilayered viscoelastic porous media are discussed by several examples. The following results can be concluded from the numerical examples: (1) The fractional derivative order can express creep well with its change, and it has no influence on initial settlements, final settlements and dissipation of excess pore fluid pressures. Besides, the rate of consolidation increases with the fractional derivative order α increasing. (2) The initial settlements are controlled by E0 , and the creep deformation is controlled by E1 . In addition, η1v has a great influence on the rate of creep. Compared with the elastic medium model, the viscoelastic medium model can reflect the creep of porous media better. (3) Layering of materials cannot be ignored. The reinforcement of the upper layer of porous media can reduce the surface settlements effectively. It is important to consider the differences of settlements near the load. Acknowledgment This research is supported by the National Natural Science Foundation of China (Grant No. 41672275). References [1] M.A. Biot, General theory of three-dimensional consolidation, J. Appl. Phys. 12 (2) (1941) 155–164. [2] M.A. Biot, General solution of the equations of elasticity and consolidation for a porous material, J. Appl. Mech. 23 (3) (1956) 91–95. [3] R.K.N.D. Rajapakse, T. Senjuntichai, Fundamental solutions for a poroelastic half-space with compressible constituents, J. Appl. Mech. 60 (4) (1993) 847–856. [4] S.J. Singh, R. Kumar, S. Rani, Quasi-static deformation of a poroelastic half-space with anisotropic permeability by two-dimensional surface loads, Geophys. J. Int. 170 (3) (2007) 1311–1327. [5] S.J. Singh, R. Kumar, S. Rani, Consolidation of a poroelastic half-space with anisotropic permeability and compressible constituents by axisymmetric surface loading, J. Earth Syst. Sci. 118 (5) (2009) 563–574. [6] J.R. Booker, J.C. Small, Finite layer analysis of consolidation i, Int. J. Numer. Anal. Methods Geomech. 6 (2) (1982) 151–171. [7] J.R. Booker, J.C. Small, Finite layer analysis of consolidation II, Int. J. Numer. Anal. Methods Geomech. 6 (2) (1982) 173–194. [8] G.J. Chen, Consolidation of multilayered half space with anisotropic permeability and compressible constituents, Int. J. Solids Struct. 41 (41) (2004) 4567–4586. [9] T. Senjuntichai, R.K.N.D. Rajapakse, Exact stiffness method for quasi-statics of a multi-layered poroelastic medium, Int. J. Solids Struct. 32 (11) (1995) 1535–1553. [10] Z.Y. Ai, Y.C. Cheng, W.Z. Zeng, Analytical layer-element solution to axisymmetric consolidation of multilayered soils, Comput. Geotech. 38 (2) (2011) 227–232.

Please cite this article as: Z.Y. Ai, Y.Z. Zhao and W.J. Liu, Fractional derivative modeling for axisymmetric consolidation of multilayered cross-anisotropic viscoelastic porous media, Computers and Mathematics with Applications (2019), https://doi.org/10.1016/j.camwa.2019.08.033.

14

Z.Y. Ai, Y.Z. Zhao and W.J. Liu / Computers and Mathematics with Applications xxx (xxxx) xxx

[11] Y.J. Chiou, S.Y. Chi, Boundary element analysis of biot consolidation in layered elastic soils, Int. J. Numer. Anal. Methods Geomech. 18 (6) (1994) 377–396. [12] J.R. Booker, M.F. Randolph, Consolidation of a cross-anisotropic soil medium, Quart. J. Mech. Appl. Math. 37 (3) (1984) 479–495. [13] S.L. Chen, L.Z. Chen, L.M. Zhang, The axisymmetric consolidation of a semi-infinite transversely isotropic saturated soil, Int. J. Numer. Anal. Methods Geomech. 29 (13) (2005) 1249–1270. [14] Y.Q. Cai, X.Y. Geng, Consolidation analysis of a semi-infinite transversely isotropic saturated soil under general time-varying loadings, Comput. Geotech. 36 (3) (2009) 484–492. [15] Z.Y. Ai, Y.C. Cheng, Extended precise integration method for consolidation of transversely isotropic poroelastic layered media, Comput. Math. Appl. 68 (12) (2014) 1806–1818. [16] D.W. Taylor, W.A. Merchant, A theory of clay consolidation accounting for secondary compressions, Stud. Appl. Math. 19 (3) (1940) 167–185. [17] L.H. Lan, K.H. Xie, Semi-analytical solution of one-dimensional visco-elastic consolidation for layered soft clay, China Civ. Eng. J. 36 (4) (2003) 105–110. [18] X.B. Li, X.L. Jia, K.H. Xie, Analytical solution of 1-D viscoelastic consolidation of soft soils under time-dependent loadings, Rock Soil Mech. 27 (Sl) (2006) 140–146. [19] K.H. Xie, X.Y. Xie, X.B. Li, Analytical theory for one-dimensional consolidation of clayey soils exhibiting rheological characteristics under time-dependent loading, Int. J. Numer. Anal. Methods Geomech. 32 (14) (2008) 1833–1855. [20] J.C. Liu, G.H. Lei, X.D. Wang, One-dimensional consolidation of visco-elastic marine clay under depth-varying and time-dependent load, Mar. Georesour. Geotechnol. 33 (4) (2015) 342–352. [21] H.E. Wahls, Analysis of primary and secondary consolidation, J. Soil Mech. Found. Div. 88 (6) (1965) 207–234. [22] G. Mesri, Y.K. Choi, Settlement analysis of embankments on soft clays, J. Geotech. Eng. 111 (4) (1985) 441–464. [23] J.H. Yin, F. Tong, Constitutive modeling of time-dependent stress–strain behaviour of saturated soils exhibiting both creep and swelling, Can. Geotech. J. 48 (12) (2011) 1870–1885. [24] Z.Y. Yin, J.H. Wang, A one-dimensional strain-rate based model for soft structured clays, Sci. China-Technol. Sci. 55 (1) (2012) 90–100. [25] T.K. Tan, Three-dimensional theory on the consolidation and flow of the clay-layers, Sci. China Math. 6 (1) (1957) 203–215. [26] Y.L. Fu, Z.J. Luo, X. Liao, J.M. Zhang, A three-dimensional full coupling model to simulate and predict land subsidence caused by high-rise building, J. Jilin Univ.-Earth Sci. Ed. 46 (6) (2016) 1781–1789. [27] W.M. Zhang, A new rheological model theory with fractional order derivatives, Natur. Sci. J. Xiangtan Univ. 23 (1) (2001) 30–36. [28] A. Gemant, A method of analyzing experimental results obtained from elasto-viscous bodies, J. Appl. Phys. 7 (1) (1936) 311–317. [29] D.S. Yin, H. Wu, C. Cheng, Fractional order constitutive model of geomaterials under the condition of triaxial test, Int. J. Numer. Anal. Methods Geomech. 37 (8) (2013) 961–972. [30] Z.Y. Liu, Q. Yang, One-dimensional rheological consolidation analysis of saturated clay using fractional order kelvin’s model, Rock Soil Mech. 38 (12) (2017) 3680–3687, 3697. [31] I.N. Sneddon, The Use of Integral Transform, McGraw-Hill, New York, 1972. [32] G.W.S. Blair, The role of psychophysics in rheology, J,. Colloid Sci. 2 (1) (1947) 21–32. [33] A.N. Gerasimov, A generalization of linear laws of deformation and its application to inner friction problems, Akad.Nauk Sssr.Prikl.Mat.Meh 12 (3) (1948) 251–260. [34] W. Chen, H.G. Sun, X.C. Li, Fractional Derivative Modeling in Mechanics and Engineering, Science Press, Beijing, 2010. [35] H. Schiessel, R. Metzler, T.F. A. Blumen, Generalized viscoelastic models: their fractional equations with solutions, J. Phys. A 28 (23) (1995) 6567–6584. [36] N. Shimizu, W. Zhang, Fractional calculus approach to dynamic problem of viscoelastic material, JSME Int. J. 42 (4) (1999) 825–837. [37] E.H. Lee, Stress analysis in visco-elastic bodies, Quart. J. Mech. Appl. Math. 21 (2) (1955) 671–674. [38] X.Z. Liu, H.H. Zhu, Back analysis on staged construction in transversely isotropic viscoelastic soil and its application to geotechnical engineering, Chin. J. Geotech. Eng. 24 (1) (2002) 89–92. [39] H.Z. Kuang, X.Z. Liu, Back analysis on transversely isotropic viscoelasticity of rock and soil medium, Chin. J. Undergr. Space Eng. 2 (1) (2006) 112–114.

Please cite this article as: Z.Y. Ai, Y.Z. Zhao and W.J. Liu, Fractional derivative modeling for axisymmetric consolidation of multilayered cross-anisotropic viscoelastic porous media, Computers and Mathematics with Applications (2019), https://doi.org/10.1016/j.camwa.2019.08.033.