A meshless method for axisymmetric problems in continuously nonhomogeneous saturated porous media

Computers and Geotechnics 62 (2014) 100–109

Contents lists available at ScienceDirect

Computers and Geotechnics journal homepage: www.elsevier.com/locate/compgeo

A meshless method for axisymmetric problems in continuously nonhomogeneous saturated porous media J. Sladek a,⇑, V. Sladek a, M. Schanz b a b

Institution of Construction and Architecture, Slovak Academy of Sciences, 84503 Bratislava, Slovakia Institute of Applied Mechanics, Graz University of Technology, Graz, Austria

a r t i c l e

i n f o

Article history: Received 15 October 2013 Received in revised form 6 June 2014 Accepted 8 July 2014

Keywords: Solid and fluid phases Coupled problem Meshless local Petrov–Galerkin method (MLPG) Moving least-squares approximation 3D axisymmetric problem Borehole

a b s t r a c t A meshless method based on the local Petrov–Galerkin approach is proposed to analyze 3-d axisymmetric problems in porous functionally graded materials. Constitutive equations for porous materials possess a coupling between mechanical displacements for solid and fluid phases. The work is based on the u–u formulation and the incognita fields of the coupled analysis in focus are the solid skeleton displacements and the fluid displacements. Independent spatial discretization is considered for each phase of the model, rendering a more flexible and efficient methodology. Both displacements are approximated by the moving least-squares (MLS) scheme. The paper presents in the first time a general meshless method for the numerical analysis of axisymmetric problems in continuously nonhomogeneous saturated porous media. Numerical results are given for boreholes in continuously nonhomogeneous porous medium with prescribed misfit and exponential variation of material parameters in the excavation zone. Ó 2014 Elsevier Ltd. All rights reserved.

1. Introduction Poroelasticity belongs to the continuum mechanics with two or more phases. The one-dimensional theory of the consolidation of a water saturated elastic porous geomaterial was first developed by Terzaghi [1]. Later Biot [2] formulated a theory for multidimensional problems of porous materials saturated by a viscous fluid. The generalized three-dimensional theory of poroelasticity in anisotropic porous materials has been developed by Biot [3] too. Theory of poroelasticity has been successfully applied in the study of variety of problems in geomechanics, biomechanics, materials engineering, environmental geomechanics and energy resource recovery from geological formations [4,5]. The extension to a nearly saturated poroelastic material has been presented by Aifantis [6] and Wilson and Aifantis [7] for the quasi-static case. The dynamic extension of Biot’s theory to three phases has been published by Vardoulakis and Beskos [8]. A state of the art overview on the theory of dynamic poroelasticity, its numerical approximation, and applications may be found in Schanz’ review paper [5]. Since the coupled differential equations are generally difficult to solve exactly, it appears that numerical approaches have to be adopted to attain solutions. Analytical methods are restricted to

⇑ Corresponding author. E-mail address: [email protected] (J. Sladek). http://dx.doi.org/10.1016/j.compgeo.2014.07.006 0266-352X/Ó 2014 Elsevier Ltd. All rights reserved.

simple boundary value problems. A nice review is given by Selvadurai [9]. Despite the universality and great success of the finite and boundary element methods in their applicability even to multifield problems, there are some restrictions leading to exclusion of the finite elements with equal order interpolation for pressure and displacements in poroelastic problems [10–13]. The dynamic Green’s function of homogeneous poroelastic half-plane has been derived by Senjuntichai and Rajapakse [14] and later applied to a vertical vibration of an embedded rigid foundation in a poroelastic soil [15]. The boundary element method (BEM) has been developed for transient and time harmonic analysis of dynamic poroelasticity problems [16]. Later a simple BEM formulation for poroelasticity via particular integrals has been developed by Banerjee [17]. Time domain BEM has been applied for axisymmetric quasi-static poroelasticity [18]. Dynamic Green’s functions for poroelastic and layered poroelastic half-spaces have been derived in [19] and [20]. In general, material coefficients in poroelasticity are anisotropic [21] and spatially variable. Axisymmetric problems have received considerable attention in the past due to their close relevance to geotechnical and rock testing methods such as uni-axial and tri-axial compression tests, double-punch tests and point load strength tests. In addition, stress analysis of cylinders is also relevant to applications involving biomedical and mechanical engineering. A cylindrical borehole drilled in a soil/rock medium is commonly found in the petroleum industry. Stability of borehole is important because it is the one of major

101

J. Sladek et al. / Computers and Geotechnics 62 (2014) 100–109

problems in oil and gas industries. In the past, the classical theory of elasticity has been used extensively to analyze various elastostatic and elastodynamic problems involving cylinders and boreholes [22]. However, geological materials are normally two-phase materials consisting of a solid skeleton with voids filled with water. Such materials are commonly known as poroelastic materials and widely considered as much more realistic representation for natural soils and rocks than ideal elastic materials [23]. The governing equations for a poroelastic material undergoing axisymmetric deformations are given by Rice and Cleary [24]. A misfit of the radial displacement with respect to the borehole diameter is considered here. A cylindrical borehole in a poroelastic medium with consideration of excavation disturbed zone is considered. Shear modulus and permeability coefficient are assumed to be continuously non-homogenous in radial direction for the disturbed zone. Vrettos [25] derived Green’s functions for vertical point load on a non-homogeneous medium with an exponential variation of the shear modulus decreasing with depth. In spite of the great success of domain and boundary discretization methods for the solution of general boundary value problems, there is still a growing interest in the development of new advanced computational methods. The finite element method (FEM) can be successfully applied to problems with an arbitrary variation of material properties by using special graded elements. In commercial computer codes, however, material properties are considered to be uniform within each element. In recent years, meshless formulations are becoming popular due to their high adaptability and easier preparation of input and output data in numerical analyses. The moving least squares (MLS) approximation is generally considered as one of many schemes to interpolate discrete data with a reasonable accuracy. The order of continuity of the MLS approximation is given by the minimum between the orders of continuity of the basis functions and that of the weight function. Thus, continuity can be tuned to a desired value. In conventional discretization methods, however, the interpolation functions usually result in a discontinuity in the secondary fields (gradients of primary fields) on the interfaces of elements. For modeling coupled fields the approach based on piecewise continuous elements can bring some inaccuracies. Therefore, a model which is based on C1-continuity, such as the meshless method, is expected to be more accurate than conventional discretization techniques. A drawback of meshless methods is higher CPU time compared to regular FEM. However, this drawback can be overcome. Recently, the authors [48–51] have developed a modified MLPG formulation, where Taylor series expansions and analytical integrations over the local sub-domains in two-dimensional elastodynamics are applied. A variety of meshless methods can be derived from a weak-form formulation either on the global domain or on a set of local subdomains. In the global formulation, background cells are required for the integration of the weak-form. In methods based on local weakform formulation, on the other hand, no background cells are required. The meshless local Petrov–Galerkin (MLPG) method is a fundamental base for the derivation of many meshless formulations, since the trial and test functions can be chosen from different functional spaces [26–29]. The MLPG method with a Heaviside step function as the test function [29] has been successfully applied to solve various 3-d axisymmetric problems [30–32]. The MLPG has been successfully applied to porous problems [33,52,53]. Many meshless formulations in poroelastic media have been applied to analyze consolidation problems [54–57]. The MLPG has been applied also to dynamic poroelastic problems, however up to day only as two-dimensional analyses [33,52]. In all early published papers based on the meshless formulations homogeneous material properties are considered. A meshfree algorithm based on the Galerkin approach is proposed for the fully coupled analysis of flow

and deformation in unsaturated poroelastic media by Khoshghalb and Khalili [34]. Temporal discretization is achieved there using a three-point approximation technique with second order accuracy. Sheu [35] analyzed the prediction of probabilistic settlements with the uncertainty in the spatial variability of Young’s modulus to illustrate the preliminary development of a spectral stochastic meshless local Petrov–Galerkin (SSMLPG) method. Generalized polynomial chaos expansions of Young’s moduli and a twodimensional meshfree weak–strong formulation in elasticity are combined to derive the SSMLPG formulation. In the present paper, the MLPG is developed for an axisymmetric 3D boundary value problem in a porous material with continuously varying material properties. It is the first meshless application to such a problem. Because of the axial symmetry, the analyzed domain is the cross-section of the considered body with the plane involving the axis of symmetry. Both governing equations for the balance of momentum in solid and fluid phases are satisfied in a weak form on small fictitious subdomains in the present paper. Nodal points are introduced and spread on the analyzed domain and each node is surrounded by a small circle for simplicity; but in general, it can be of an arbitrary shape. The spatial variations of the displacements in solid and fluid phases are approximated by the moving least-squares scheme [36,37]. After performing the spatial integrations, one obtains a system of ordinary differential equations (ODE) for temporal variations of certain nodal unknowns. The backward difference method is applied for the approximation of ‘‘velocities’’ and the Houbolt method [38] is applied for the accelerations in the ODE. The influence of the material gradation on displacements and stresses in porous medium around the borehole is investigated. 2. Local boundary integral equations In Biot’s theory a fully saturated material is considered, i.e., an elastic skeleton with a statistical distribution of interconnected pores is modeled. The porosity is denoted by

/¼

Vf V

ð1Þ

where V f is the volume of the interconnected pores contained in a sample of bulk volume V. The sealed pores are considered as part of the solid. Full saturation is assumed leading to V = V f + V s with Vs the volume of the solid. There are several possibilities to write governing equations for porous materials: (i) To use the solid displacement usi and the fluid displacement ufi with six (four) unknowns in 3-d (2-d) problems [39]. (ii) Alternatively the solid displacement usi and the seepage displacement wi (wi ¼ /ðufi  usi Þ) with also six (four) unknowns in 3-d (2-d) can be used [40]. Beside the seepage displacement also sometimes the seepage velocity, i.e., the time derivative of wi is applied. (iii) A combination of the pore pressure p and the solid displacement usi with four (three) unknowns in 3-d (2-d) can be established. As shown by Bonnet [41], this choice is sufficient. In the present paper we use the first approach based on solid and fluid displacements. If the constitutive equations are formulated for the elastic solid and the interstitial fluid, a partial stress formulation is obtained [2,3] s ij

s ij

r ¼ 2Ge

2 Q2 þ K  Gþ 3 R

rf ¼ /p ¼ Q eskk þ Refkk ;

!

eskk dij þ Q efkk dij ;

ð2Þ

ð3Þ

102

J. Sladek et al. / Computers and Geotechnics 62 (2014) 100–109

where strain tensors are given by

esij

 1 s ¼ ui;j þ usj;i and 2

3=z

efkk ¼ ufk;k :

The respective stress tensors are denoted by rsij and rfdij. The elastic skeleton is assumed to be isotropic and homogeneous with the two elastic material constants bulk modulus K and shear modulus G. The coupling between the solid and the fluid is characterized by the two parameters Q and R. Considerations of constitutive relations at micro mechanical level as given by Detournay and Cheng [42] shows that these parameters can be calculated from the bulk moduli of the constituents by

R¼

K f ðK s  KÞ þ /K s ðK s  K f Þ f

/½ð1  /ÞK  KK K

2

s

ð5Þ

;

where Ks denotes the bulk modulus of the solid grains and Kf the bulk modulus of the fluid. An alternative representation of the constitutive Eqs. (2) and (3) is used in Biot’s earlier work [2]. There, the total stress rij ¼ rsij þ rf dij is introduced and with Biot’s effective stress coefficient

  K Q ; s ¼ / 1þ R K

ð6Þ

the constitutive equation can be written as



2 3

ð7Þ

The governing equations for porous materials can be written as the balance of momentum for total stresses and for fluid [39]:

rij;j ðx; sÞ þ ð1  /Þfis þ /fif ¼ ð1  /Þqs u€si ðx; sÞ þ /qf u€fi ðx; sÞ; h

rf;i ðx; sÞ þ /fif ¼ /qf u€fi ðx; sÞ þ qA u€fi ðx; sÞ  u€si ðx; sÞ

fis

/2 h

j

1=r Fig. 1. Geometry and generating section of a circular cylinder of radius a and height h.

rsrr 1 2 c11 B rs C 6 c B uu C 6 12 B s C¼6 @ rzz A 4 c13 rsrz 0 0



rij ¼ 2Gesij þ K  G eskk dij  adij p ¼ rijeffectiv e  adij p:

þ

Ω

ð4Þ

;

K f ðK s  KÞ þ /K s ðK s  K f Þ

a¼1

Γ

h

/2 K f K s K s

s

Q¼

a

u_ fi ðx;

sÞ 

u_ si ðx;

ð8Þ

c13

0

30

usr;r

1

c11

c13

usr =r usz;z

C 6Q 7 C 6 7 f C þ 6 7ekk ; A 4Q 5

usr;z þ usz;r

0

c13

c11

B 0 7 7B 7B 0 5@

0

0

c44

0

usr;r

0

1

B

ufr;r

2

Q

3 ð10Þ

1 C

C f C rf ¼ ½ Q Q Q B @ usr =r A þ ½ R R R B @ ur =r A;

ð11Þ

rab ¼ rsab þ dab rf ;

ð12Þ

usz;z

ufz;z

where

i

i sÞ ;

c12

c11 ¼ ð9Þ

fif

where and are denote the volume densities of body forces in solid and fluid medium, respectively. Further, the respective mass densities are denoted by qs and qf. The apparent mass density qA is introduced [39] to describe the dynamic interaction between solid and fluid phases. Usually it has the following form, qA = C/ qf, where C is a factor depending on the geometry of the pores and the frequency of excitation. The term /2/j is the damping factor valid for circular pores, when j denotes the permeability. For an axisymmetric problem in a porous piezoelectric material, it is convenient to use cylindrical coordinates x  ðr; u; zÞ. The axis of symmetry is identical with the z-axis. The angular component of the displacements vanishes and all physical fields as well as material coefficients are independent of the angular coordinate u. A cylinder of height h and radius a is generated by the rotation of the planar domain X bounded by the boundary C around the axis of symmetry as shown in Fig. 1. Thus, in the cylindrical coordinates the nonzero strains are given as

4 Q2 GþK þ ; 3 R

2 Q2 c12 ¼ c13 ¼ K  G þ ; 3 R

c44 ¼ G:

The following essential and natural boundary conditions are assumed for the mechanical field in the solid phase

~ si ðx; sÞ; usi ðx; sÞ ¼ u

on Cu ;

ti ðx; sÞ  rij nj ¼ ~t i ðx; sÞ;

on Ct ;

C ¼ Cu [ Ct ;

whilst for the fluid phase

~ fi ðx; sÞ; ufi ðx; sÞ ¼ u

on Cfu ;

rf ðx; sÞ ¼ r~ f ðx; sÞ; on Cft ; C ¼ Cfu [ Cft : where Cu is the part of the global boundary C with prescribed displacements, while on Ct the traction vector is given. The initial conditions for mechanical displacements are assumed as

  usi ðx; sÞs¼0 ¼ usi ðx; 0Þ and u_ si ðx; sÞs¼0 ¼ u_ si ðx; 0Þ   ufi ðx; sÞ

s¼0

  ¼ ufi ðx; 0Þ and u_ fi ðx; sÞ

s¼0

¼ u_ fi ðx; 0Þ in X

esrr ¼ usr;r ; esuu ¼ usr =r; esrz ¼ ðusr;z þ usz;r Þ=2; eszz ¼ usz;z ; efkk ¼ ufr;r þ ufz;z þ ufr =r:

Assuming for vanishing body forces, the governing Eqs. (6) and (7) in the cylindrical coordinate system take the form

Generally, material properties are varying with Cartesian coordinates. In such a case matrices C, Q and R in the following constitutive equations are functions of (r, z)

rrb;b ðr; z; sÞ þ ½rrr ðr; z; sÞ  ruu ðr; z; sÞ  ð1  /Þqs ðxÞu€sr ðr; z; sÞ

1 r € fr ðr; z; sÞ ¼ 0; þ qf ðxÞ/u

ð13Þ

J. Sladek et al. / Computers and Geotechnics 62 (2014) 100–109

1 r € fz ðr; z; sÞ ¼ 0; þ qf ðxÞ/u

rzb;b ðr; z; sÞ þ rrz ðr; z; sÞ  ð1  /Þqs ðxÞu€sz ðr; z; sÞ ð14Þ

rf;b ðr; z; sÞ  /qf ðxÞu€fb ðr; z; sÞ  qA ðxÞ½u€fb ðr; z; sÞ  u€sb ðr; z; sÞ /2



½u_ fb ðr; z; sÞ  u_ sb ðr; z; sÞ ¼ 0;

j

ð15Þ

103

where pT(x) = [p1(x), p2(x), . . ., pm(x)] is a complete monomial basis of order m; and a(x, s) is a vector containing the coefficients aj(x, s), j = 1, 2, . . ., m and x  (r, z). For the considered axis-symmetric case, the monomial basis on the cross section is expressed by the radial and axial coordinates (r, z)

pT ðxÞ ¼ ½1; r; z; for linear basis m ¼ 3;

ð23Þ

pT ðxÞ ¼ ½1; r; z; r2 ; rz; z2 ; for quadratic basis m ¼ 6:

ð24Þ

with the summation convention being assumed for repeated subscript b e {r, z}. The local weak form of the governing Eqs. (13)–(15) can be written as

The coefficient vector a(x, s) is determined by minimizing a weighted discrete L2-norm defined as

Z

JðxÞ ¼

Z

rrb;b U 1 dX þ

Xs

þ

Xs

Z Xs

Z

Xs

Z Xs

Z Xs

Z

þ

Z

Z Xs

h

Xs

Z



€s s ur ðr; z;

ð1  /Þq

Xs

Xs

sÞU 1 dX ð16Þ

1 rrz ðr; z; sÞU 2 dX  r

Z

€ sz ðr; z; sÞU 2 dX ð1  /Þqs u

Xs

qf /u€ fz ðr; z; sÞU 2 dX ¼ 0

rf;b U 3 dX  

Z

qf /u€ fr ðr; z; sÞU 1 dX ¼ 0;

rzb;b U 2 dX þ

Xs

1 ðrrr  ruu ÞU 1 dX  r

ð17Þ

@ Xs

AðxÞ ¼

i /2 h f u_ b ðr; z; sÞ  u_ sb ðr; z; sÞ U 3 dX ¼ 0;

j

þ

Xs

Z

þ

Z Xs



Z Xs

Z Xs

Xs

Z

1 rrz ðr; z; sÞdX  r

Xs

/2 h

j

€ sr ðr; z; sÞdX ð1  /Þqs u

Z Xs

€ sz ðr; z; sÞdX ð1  /Þqs u ð20Þ

€ fb ðr; z; sÞdX  /qf u i

Z Xs

h

i

qA u€fb ðr; z; sÞ  u€ sb ðr; z; sÞ dX

u_ fb ðr; z; sÞ  u_ sb ðr; z; sÞ dX ¼ 0:

ð21Þ

In the present paper the trial functions are approximated by the MLS method on a number of nodes spread over the influence domain. According to the MLS [36,37] method, the approximation of physical fields u(x, s) (displacements for both phases) over a number of randomly located nodes {xa}, a = 1, 2, . . ., n, is given by

uðx; sÞ ¼ pT ðxÞaðx; sÞ;

ð26Þ

n X wa ðxÞpðxa ÞpT ðxa Þ; a¼1

ð18Þ

ð19Þ

qf /u€ fz ðr; z; sÞdX ¼ 0;

rf nb dC 

@ Xs

Z

1 ðrrr  ruu ÞdX  r

qf /u€ fr ðr; z; sÞdX ¼ 0;

rzb ðr; z; sÞnb dC þ

@ Xs

Z

Z Xs

Z

where n is the number of nodes used for the approximation. It is determined by the weight function wa(x) associated with the node a. A 4th order spline-type weight function is applied in the present ^ a ðsÞ stands for the fictitious nodal values, but work. The symbol u not the nodal values of the unknown trial functions in general. The stationarity of J in Eq. (25) with respect to a(x, s) leads to the follow^ ðsÞ ¼ ½u ^ 1 ðsÞ; . . . ; u ^ n ðsÞT ing linear relation between a(x, s) and u

where

i qA u€ fb ðr; z; sÞ  u€ sb ðr; z; sÞ U 3 dX

rrb ðr; z; sÞnb dC þ

ð25Þ

a¼1

^ ðsÞ ¼ 0; AðxÞaðx; sÞ  BðxÞu € fb ðr; z; sÞU 3 dX /qf u

where Uk(x) are test functions and oXs is the boundary of the local subdomain which consists of three parts oXs = Ls [ Cst [ Csu [37]. Here, Ls is the local boundary which is totally inside the global domain, Cst is the part of the local boundary which coincides with the global traction boundary, i.e., Cst = oXs \ Ct. Similarly Csu is the part of the local boundary that coincides with the global displacement boundary, i.e., Csu = oXs \ Cu. Applying the Gauss divergence theorem to the first domain integrals in Eqs. (16)–(18) and assuming the test functions to be given by the Heaviside unit step functions within the subdomain, one can recast the above equations into the following local integral equations

Z

n X  2 ^ a ðsÞ ; wa ðxÞ pT ðxa Þaðx; sÞ  u

ð22Þ

  BðxÞ ¼ w1 ðxÞpðx1 Þ; w2 ðxÞpðx2 Þ; . . . ; wn ðxÞpðxn Þ :

ð27Þ

The solution of Eq. (26) for a(x, s) and a subsequent substitution into Eq. (22) gives the approximation formulas for the displacements in solid and fluid phases [37]

us ðx; sÞ ¼

n X ^ sa ðsÞ; N a ðxÞu a¼1

uf ðx; sÞ ¼

n X ^ fa ðsÞ; Na ðxÞu

ð28Þ

a¼1 T ^ sa ðsÞ ¼ ðu ^ fa ðsÞ ¼ ^ sa ^ sa where the nodal values u and u r ðsÞ; uz ðsÞÞ  T ^ fa ^ fa u are fictitious parameters for the displacements in r ðsÞ; u z ðsÞ

solid and fluid phases, respectively, and Na(x) is the shape function associated with node a. The number of nodes n used for the approximation is determined by the weight function wa(x). A 4th order spline-type weight function is applied in the present work as below

8       < 1  6 da 2 þ 8 da 3  3 da 4 ; 0 6 da 6 r a a a a r r ra w ðxÞ ¼ ; : a 0; d P ra

ð29Þ

where da = kx  xak and ra is the size of the support domain. It is seen that the C1 – continuity is ensured over the entire domain, and therefore the continuity of gradients of approximated fields of displacements and pore pressure is satisfied. In the MLS approximation the rate of the convergence of the solution may depend upon the nodal distance as well as the size of the support domain [43]. It should be noted that a smaller size of the subdomains may induce larger oscillations in the nodal shape functions [37]. A necessary condition for a regular MLS approximation is that at least m weight functions are non-zero (i.e. n P m) for each sample point x e X. This condition determines the size of the support domain. Substituting the approximations given by Eq. (28) into the local integral Eqs. (19)–(21), we obtain a system of ordinary differential

104

J. Sladek et al. / Computers and Geotechnics 62 (2014) 100–109

equations (ODE) n o ^ sa ^ sa ^ fa ^ fa : u r ; uz ; ur ; uz

Z n X ^ sa u r ð sÞ a¼1

@ Xs

for

the

fictitious

ðc11 ðxÞ þ QðxÞÞnr ðxÞN a;r ðxÞ þ

þc44 ðxÞnz ðxÞN a;z ðxÞ

io

dC þ

Z Xs

unknown

parameters

c12 ðxÞ þ QðxÞ nr ðxÞN a ðxÞ r

 1 1 ðc11 ðxÞ  c12 ðxÞÞðN a;r ðxÞ  N a ðxÞÞ dX r r

Z n X €^ sa ðsÞ ð1  /Þq ðxÞN a ðxÞdX  u s r a¼1

Xs

a¼1

@ Xs

Z n X ^ sa þ u z ð sÞ Z n X ^ fa þ u r ð sÞ þ þ

@ Xs

a¼1 n X

^ fa u z ð sÞ

a¼1 n X

€^ fa ðsÞ u r

Z

@ Xs

Z

Xs

a¼1

Z n X ^ sa ð s Þ u z a¼1

h

3. Numerical examples

i ðc13 ðxÞ þ Q ðxÞÞnr ðxÞN a;z ðxÞ þ c44 ðxÞnz ðxÞNa;r ðxÞ dC

1 ðQ ðxÞ þ RðxÞÞnr ðxÞðN a;r ðxÞ þ N a ðxÞÞd

C

r

ðQ ðxÞ þ RðxÞÞnr ðxÞN a;z ðxÞdC

qf ðxÞ/ðxÞNa ðxÞdX ¼ 0; h

@ Xs

ð30Þ

i ðc11 ðxÞ þ QðxÞÞnz ðxÞNa;z ðxÞ þ c44 ðxÞnr ðxÞNa;r ðxÞ dC

X Z n c44 ðxÞ a €^ sa ðsÞ ð1  /Þq ðxÞNa ðxÞdX u N;r ðxÞdX  þ s z r Xs Xs a¼1

Z n X ^sa þ ðc44 ðxÞnr ðxÞNa;z ðxÞ þ ðc13 ðxÞ þ Q ðxÞÞnz ðxÞðNa;r ðxÞ u r ðsÞ Z

@ Xs

a¼1 Z 1 a c44 ðxÞ a þ N ðxÞÞÞdC þ N ;z ðxÞdX r r Xs Z n X 1 ^fa þ ðQ ðxÞ þ RðxÞÞnz ðxÞðNa;r ðxÞ þ Na ðxÞÞdC u r ðsÞ r @ Xs a¼1 Z n X ^fa u þ ðQ ðxÞ þ RðxÞÞnz ðxÞN a;z ðxÞdC z ðsÞ

þ

a¼1 n X

€^fa ðsÞ u z

Z

a¼1

qf ðxÞ/ðxÞNa ðxÞdX ¼ 0;

@ Xs

ð31Þ

a¼1

@ Xs

  1 nb ðxÞRðxÞ Na;r ðxÞ þ Na ðxÞ dC r @ Xs a¼1 Z n X ^ fa þ nb ðxÞRðxÞNa;z ðxÞdC u z ðsÞ Z n X ^ fa u þ r ðsÞ

  

n h X a¼1 n h X a¼1 n X

€^ fa ðsÞ  u €^ sa ðsÞ u b b ^_ fa ^_ sa u b ðsÞ  ub ðsÞ

€^ fa ðsÞ u b

a¼1

Z Xs

iZ Xs

iZ

Xs

qA ðxÞNa ðxÞdX /2

j

ðxÞNa ðxÞdX

/qf ðxÞNa ðxÞdX ¼ 0;

b 2 fr; zg:

ð32Þ

Eqs. (30)–(32) are considered in the subdomains Xs around each interior node xs and at the boundary nodes with prescribed natural boundary conditions (Ct and Cft ). On the parts of the global boundary Cu and Cfu with prescribed mechanical displacements the collocation Eq. (28) are applied. n X ^ sa ðsÞ ¼ u ~ s ðf; sÞ for f 2 Cu ; Na ðfÞu

ð33Þ

a¼1 n X ^ fa ðsÞ ¼ u ~ f ðf; sÞ for f 2 Cfu ; Na ðfÞu a¼1

R ¼ 4:7  108 N m2 ;

a ¼ 0:8; / ¼ 0:19; j ¼ 1:9  1010 m4 =Ns:

  1 nb ðxÞQ ðxÞ Na;r ðxÞ þ Na ðxÞ dC r @ Xs Z n X ^ az ðsÞ þ nb ðxÞQðxÞNa;z ðxÞdCþ u a¼1

G ¼ 6:0  109 N m2 ;

qs ¼ 2800 kg m3 ; qf ¼ 1000 kg m3 ;

Z n X ^ sa u r ðsÞ a¼1

In the first example, a finite cylinder (Fig. 1) with porous properties is analyzed. The cylinder is under a uniform pressure load on the top surface, r33 = 104 Pa with Heaviside time variation. The bottom surface of the cylinder is fixed in the axial direction, us3 ¼ uf3 ¼ 0. On the lateral sides of the cylinder, i.e., at r = a, the axial component of the traction vector and the radial displacement are vanishing. In such a case the axial symmetric problem corresponds to 1-d problem, where analytical solution is available [44]. Analytical solution can be derived in the Laplace transform domain. A numerical inverse transformation (CQM) is applied to the analytical solution. The radius of the cylinder a = 1 m and its thickness h = 1 m are considered. The mechanical displacements for both solid and fluid phases on the cross section X of the cylinder are approximated by using 121 (11  11) equidistantly distributed nodes. The local subdomains are selected to be circular with a radius rloc = 0.08 m. The material coefficients correspond to Berea sandstone [5]:

K ¼ 8:0  109 N m2 ;

@ Xs

Xs

By collecting the discretized local integral equations together with the discretized boundary conditions for the displacements in solid and fluid phases, we arrive at a complete system of ordinary differential equations. The backward difference method is applied for the approximation of ‘‘velocities’’ in the ODE (30)– (32). The Houbolt method [38] is applied for the accelerations in the same ODE.

Coefficient Q is computed from Eq. (6). Firstly, homogeneous material properties are considered. We have tested computer code for various values of the coupling parameter Q in Berea sandstone material. Other material parameters correspond to the normal Berea sandstone given above. Fig. 2 presents a temporal variation of the axial displacement on the top of cylinder for Q = 0. One can observe a good agreement of the present MLPG and analytical results. Time step Dt = 0, 1  105 s has been selected in numerical analyses. It satisfies the Courant criterion for the time step Dt 6 Dx=c1 , where c1 is the velocity of propagation of dilatational pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi waves in the solid ðc1 ¼ ð3K þ 4GÞ=3qs ¼ 2:39  103 m s1 Þ and Dx ¼ a=11  0:1 m is the distance between two nodes in the spatial discretization. In the next example we have considered the Berea sandstone with Q = 1.675  109 Nm2, which corresponds a = 0.867. Numerical results are presented in Fig. 3. Finally, we have analyzed also normal Berea sandstone cylinder. MLPG and analytical results are presented in Fig. 4. In all material model cases one can observe a good agreement of results. Fig. 5 presents a comparison of axial displacements for various bulk modulus Q. The largest axial displacement is observed for incompressible material model (Q = 0). With increasing value of the bulk modulus Q the axial displacement decreases and frequency of oscillations increases due to higher elastic wave velocity. In the next numerical example we analyze the water saturated soil cylinder, where material properties are following [45]:

K 0 ¼ 2:1  108 N m2 ; R0 ¼ 1:206  109 N m2 ;

ð34Þ

G0 ¼ 9:8  107 N m2 ; Q 0 ¼ 1:259  109 N m2 ;

qs ¼ 2700 kg m3 ; qf ¼ 1000 kg m3 ;

J. Sladek et al. / Computers and Geotechnics 62 (2014) 100–109

Fig. 2. Temporal variation of axial displacement on the top cylinder for material with Q = 0.

105

Fig. 4. Temporal variation of axial displacement on the top cylinder for normal Berea sandstone with Q = 1.511  109 N m2.

Fig. 3. Temporal variation of axial displacement on the top cylinder for material with Q = 1.675  109 N m2.

/0 ¼ 0:48;

Fig. 5. Influence of the bulk modulus Q on the axial displacement on the top surface.

j0 ¼ 3:55  109 m4 =Ns:

In a real soil the material properties are inhomogeneous. Essentially, the shear modulus is dependent on the stress level. However, for small strains (stress) this effect can be neglected in x–y plane. Therefore, for simplicity it can be considered only a dependence of the shear modulus on the depth. We consider the gradation with exponential variation of the shear modulus along axial direction like Vrettos [25]:

GðzÞ ¼ G0 ½1 þ expðm1 ðh  zÞÞ0:5:

ð35Þ

According to the material inhomogeneity model, the shear modulus on the top surface (z = h) is equal to the value corresponding to homogeneous case G0. For exponent m1 = 0.5 m1 we get G(0) = 0, 803G0 on the bottom surface. The geometry and boundary conditions are the same as in the previous examples. The problem is analyzed as 1-d problem with vanishing normal displacements on vertical surfaces of the cylinder. Time variations of the radial displacement on top surface of the cylinder for homogeneous and inhomogeneous cylinder are presented in Fig. 6. The radii of the support domain in the MLS and of the local integration subdo3 1 3 1 main are set to cdi and adi , respectively, where di and di are the distances to the third and first nearest points from node i, respectively. In all the applications that follow, c = 4 and a = 0.6, have

Fig. 6. Influence of the material gradation in soil on the axial displacement on the top surface.

106

J. Sladek et al. / Computers and Geotechnics 62 (2014) 100–109

Fig. 7. Temporal variation of the axial displacement on the top surface for free and fixed cylinder on the vertical surface.

Fig. 8. Temporal variation of the pore pressure on the bottom surface for free and fixed cylinder on the vertical surface.

been selected. One can observe large amplitudes of displacements for the inhomogeneous cylinder due to smaller stiffness. Peak values are shifted to larger instants since wave velocities are smaller. It is observed a good agreement of the present u-u formulation results with analytical ones for a porous cylinder under impact load. Numerical results illustrate that the present formulation is suitable for dynamic problems. The consolidations for unidirectional (fixed normal displacement) and axial symmetric problems (traction free surfaces) in water saturated soil are mutually compared here. Boundary conditions and geometry for unidirectional problems are defined in the previous example. For axisymmetric problem the top and lateral boundaries have vanishing pore pressure. The cylinder is under a uniform static pressure load on the top surface, r33 = 104 Pa. A quasi-static character of the balance of momentum is considered here. It means that the velocity of elastic waves is infinite and mass density has to be vanishing. The axial displacement at the center of the cylinder on the top surface (r = 0, z = h) and pore of pressure on the bottom surface (r = 0, z = 0) are computed. Numerical results are presented in Figs. 7 and 8. Higher axial displacements are observed for the cylinder with free traction on outer vertical surface. Both displacements are approaching to static quantities, uz = 0.393  104 m for axisymmetric case and uz = 0.2935  104 m for 1-d case. However, the pore pressure in this case is lower than for fixed cylinder with impermeable boundary conditions. Pore pressure is decaying for both boundary condition cases. One can observe a strong influence of displacement and pore pressure conditions on outer surface of the cylinder on the axial displacement and pore pressure. Now we consider the problem of an infinite cylindrical borehole drilled in a porous rock layer subjected to axisymmetric radial load as shown in Fig. 9. The excavation process affects the soil/rock properties around a borehole. The shear modulus of the excavation disturbed zone (EDZ) is normally reduced from its original value before excavation. On the contrary, the inflow in the EDZ is larger than that of the undisturbed zone [46,47]. Kaewjuea [23] assumed that the shear modulus and permeability coefficient in the EDZ linearly vary with the radial distance r. However, an exponential variation seems to be more realistic. Therefore, in this paper an exponential variation for both parameters is assumed:

Fig. 9. Borehole with rigid cylinder with radial misfit Dr.

J. Sladek et al. / Computers and Geotechnics 62 (2014) 100–109

lðrÞ ¼ l0 expðm1 ðr  DÞÞ;

ð36Þ

jðrÞ ¼ j0 expðm2 ðr  DÞÞ;

ð37Þ

where l0 and j0 denote the original values of the shear modulus and the permeability coefficient respectively before excavation. Parameters m1 and m2 are non-negative constants representing the gradation of disturbance due to the drilling process in shear modulus and permeability coefficient respectively. It is assumed that both material parameters have their original values l0 and j0 at sufficiently large distance r = D from the center of the borehole. On segment of the borehole surface with length 2h = 2 m is considered a rigid cylinder radial misfit Dr = 105 m. The infinite space is replaced by a finite cylinder with truncated size, where radius D = 6 m, finite high 2H = 10 m and borehole radius a = 1 m. Due to symmetry with respect to x-y plane, it is sufficient to analyze only a symmetric part, see Fig. 9. Impermeable boundary conditions are assumed on outer boundaries and borehole surface is drained. On truncated surfaces normal displacements are vanishing. A quasi-static character of the balance of momentum is considered here. Time variation of the pore pressure and angular stresses at two different radial coordinates on plane z = 0 are presented in Figs. 10

107

and 11, respectively. Both graphs are valid for the homogeneous soil. The pore pressure is growing with time and approaching the static equivalent which is uniform with radius. Angular stresses are relaxing in time and approaching static equivalents, which are different at considered points. The angular stresses at the point far from the borehole surface are negative. Temporal variation of the radial displacement at node with radius r = 2 m and plane z = 0 is presented in Fig. 12. One can observe that radial displacement is approaching a smaller value for large instants in inhomogeneous soil, where m1 = 0.1 m1 and m2 = 0.2 m1 are considered. It is mainly due to increasing value of the shear modulus in inhomogeneous soil with respect to homogeneous one. The time variation of the pore pressure at node with radius r = 2 m and plane z = 0 for both homogeneous and inhomogeneous soils is presented in Fig. 13. It seems that inhomogeneity for the permeability coefficient has a small influence on the pore pressure since both curves are almost parallel (see Fig. 14). One can observe the influence of inhomogeneity on the angular stresses on the borehole surface (r = 1 m, z = 0). The angular stresses are decreasing with time and approaching static equivalents. A smaller value is observed for the inhomogeneous soil.

Fig. 10. Temporal variation of the pore pressure at z = 0 in a homogeneous soil. Fig. 12. Temporal variation of the radial displacement at z = 0 and r = 2 m.

Fig. 11. Temporal variation of angular stresses at z = 0 in a homogeneous soil.

Fig. 13. Influence of inhomogeneity on the pore pressure at z = 0 and r = 2 m.

108

J. Sladek et al. / Computers and Geotechnics 62 (2014) 100–109

Fig. 14. Influence of inhomogeneity on the angular stresses at z = 0 and r = 1 m.

4. Conclusions A meshless local Petrov–Galerkin method (MLPG) is presented for the 3D axisymmetric problem in continuously nonhomogeneous saturated porous media. Stationary and transient dynamic conditions are considered. The present u-u formulation seems to be a powerful tool to analyze a general boundary value 3D axisymmetric problem in dynamic poroelasticity. The moving least-squares (MLS) scheme is adopted for the approximation of the physical field quantities. The order of continuity of the MLS approximation is given by the minimum between the orders of continuity of the basis functions and that of the weight function. The present model has the C1-continuity and one can expect a higher accuracy than in conventional discretization techniques for a medium with continuously varying material properties. The proposed method is a truly meshless method, which requires neither domain elements nor background cells in either interpolation or integration. The proposed method is applied to a porous cylinder and a borehole problem with prescribed misfit radial displacement on a part of hole. For the cylinder an exponential variation of the shear modulus is considered along the axial coordinate. Large amplitudes of displacements for inhomogeneous cylinder due to smaller stiffness are observed and peak values are shifted to larger instants since wave velocities are smaller. The excavation disturbed zone is considered for the borehole problem, where the shear modulus and permeability coefficient are exponentially varying along the radial coordinate r. The consolidation problem with static boundary conditions and vanishing mass densities is analyzed. The influence of inhomogeneity on the radial displacement and pore pressure is less than 10% for considered material inhomogeneity. A large influence is observed only on the angular stresses on the borehole surface. Acknowledgement The first two authors acknowledge the support by the Slovak Science and Technology Assistance Agency registered under number APVV-0014-10 and APVV-0032-10. References [1] Terzaghi K. Erdbaumechanik auf bodenphysikalischer Grundlage. Leipzig: Franz Deuticke; 1925. [2] Biot MA. General theory of three-dimensional consolidation. J Appl Phys 1941;12:155–64. [3] Biot MA. Theory of elasticity and consolidation for a porous anisotropic solid. J Appl Phys 1955;26:182–5.

[4] Selvadurai APS. Mechanics of poroelastic media. Dordrecht: Kluwer Academic; 1996. [5] Schanz M. Poroelastodynamics: linear models, analytical solutions, and numerical methods. Appl Mech Rev 2009;62:030803–30815. [6] Aifantis EC. On the problem of diffusion in solids. Acta Mech 1980;37:265–96. [7] Wilson RK, Aifantis EC. On the theory of consolidation with double porosity. Int J Eng Sci 1982;20:1009–35. [8] Vardoulakis I, Beskos DE. Dynamic behavior of nearly saturated porous media. Mech Matls 1986;5:87–108. [9] Selvadurai APS. The analytical method in geomechanics. Appl Mech Rev 2007;60:87–106. [10] Lewis RW, Schrefler BA. The finite element method in the static and dynamic deformation and consolidation of porous media. 2nd ed. Chichester: John Wiley & Sons; 1998. [11] Zienkiewicz OC, Chan AHC, Pastor M, Schrefler BA, Shiomi T. Computational geomechanics with special reference to earthquake engineering. Chichester: John Wiley & Sons; 1999. [12] Li X, Han X, Pastor M. An iterative stabilized fractional step algorithm for finite element analysis in saturated soil dynamics. Comput Method Appl Mech Eng 2003;192:3845–59. [13] Soares Jr D. A time-domain FEM approach based on implicit Green’s functions for the dynamic analysis of porous media. Comput Method Appl Mech Eng 2008;197:4645–52. [14] Senjuntichai T, Rajapakse RKND. Dynamic Green’s functions of homogeneous poroelastic half-plane. J Eng Mech ASCE 1994;120:2381–404. [15] Senjuntichai T, Mani S, Rajapakse RKND. Vertical vibration of an embedded rigid foundation in a poroelastic soil. Soil Dyn Earthquake Eng J 2006;26:626–36. [16] Chen J, Dargush GF. Boundary element method for dynamic poroelastic and thermoelastic analyses. Int J Solids Struct 1995;32:2257–78. [17] Park KH, Banerjee PK. A simple BEM formulation for poroelasticity via particular integrals. Int J Solids Struct 2006;43:3613–25. [18] Dargush GF, Banerjee PK. A boundary element method for axisymmetric soil consolidation. Int J Solids Struct 1991;28:897–915. [19] Norris AN. Dynamic Green’s functions in anisotropic piezoelectric, thermoelastic and poroelastic solids. Proc R Soc Lond A 1994;447:175–88. [20] Pan E. Green’s function in layered poroelastic half-spaces. Int J Numer Anal Methods Geomech 1999;23:1631–53. [21] Cheng AHD. Material coefficients of anisotropic poroelasticity. Int J Rock Mech Min Sci 1997;34:199–205. [22] Rajapakse RKND, Gross D. Traction and contact problems for an anisotropic medium with a cylindrical borehole. Int J Solids Struct 1996;33:2193–211. [23] Kaewjuea W. Poroelastic solutions for borehole and cylinder. PhD Thesis, Chulalongkorn University, 2010. [24] Rice JR, Cleary MP. Some basic stress-diffusion solutions for fluid saturated elastic porous media with compressible constituents. Rev Geophys Space Phys 1976;14:227–41. [25] Vrettos Ch. Green’s functions for vertical point load on an elastic half-space with depth-degrading stiffness. Eng Anal Boundary Elements 2008;32:1037–45. [26] Zhu T, Zhang JD, Atluri SN. A local boundary integral equation (LBIE) method in computational mechanics, and a meshless discretization approach. Comput Mech 1998;21:223–35. [27] Atluri SN, Sladek J, Sladek V, Zhu T. The local boundary integral equation (LBIE) and its meshless implementation for linear elasticity. Comput Mech 2000;25:180–98. [28] Sladek J, Sladek V, Atluri SN. Local boundary integral equation (LBIE) method for solving problems of elasticity with nonhomogeneous material properties. Comput Mech 2000;24:456–62. [29] Sladek J, Stanak P, Han ZD, Sladek V, Atluri SN. Applications of the MLPG method in engineering & sciences: a review. CMES Comput Model Eng Sci 2013;92:423–75. [30] Sladek J, Sladek V, Krivacek J, Zhang Ch. Meshless local Petrov–Galerkin method for stress and crack analysis in 3-D axisymmteric FGM bodies. CMES Comput Model Eng Sci 2005;8:259–70. [31] Sladek J, Sladek V, Hellmich Ch, Eberhardsteiner J. Heat conduction analysis of 3D axisymmetric and anisotropic FGM bodies by meshless local Petrov– Galerkin method. Comput Mech 2007;39:323–33. [32] Sladek J, Sladek V, Solek P, Saez A. Dynamic 3-D axisymmetric problems in continuously nonhomogeneous piezoelectric solids. Int J Solids Struct 2008;45:4523–42. [33] Soares Jr D, Sladek V, Sladek J, Zmindak M, Medvecky S. Porous media analysis by modified MLPG formulations. CMC – Comput Mater Contin 2012;27:101–26. [34] Khoshghalb A, Khalili N. A mesfree method for fully coupled analysis of flow and deformation in unsaturated porous media. Int J Num Anal Meth Geomech 2013;37:716–43. [35] Sheu GY. Prediction of probabilistic settlements via spectral stochastic meshless local Petrov–Galerkin methods. Comput Geotech 2011;38:407–15. [36] Belytschko T, Krogauz Y, Organ D, Fleming M, Krysl P. Meshless methods; an overview and recent developments. Comput Method Appl Mech Eng 1996;139:3–47. [37] Atluri SN. The meshless method, (MLPG) for domain & BIE discretizations. Forsyth: Tech Science Press; 2004. [38] Houbolt JC. A recurrence matrix solution for the dynamic response of elastic aircraft. J Aeronaut Sci 1950;17:371–6.

J. Sladek et al. / Computers and Geotechnics 62 (2014) 100–109 [39] Biot MA. Theory of propagation of elastic waves in a fluid-saturated porous solid, I. Low frequency range. J Acoust Soc Am 1956;28:168–78. [40] Zienkiewicz OC, Shiomi T. Dynamic behaviour of saturated porous media; the generalized Biot formulation and its numerical solution. Int J Numer Anal Method Geomech 1984;8:71–96. [41] Bonnet G. Basic singular solutions for a poroelastic medium in the dynamic range. J Acoust Soc Am 1987;82:1758–62. [42] Detournay E, Cheng AHD. Fundamentals of poroelasticity. Comprehensive rock engineering: principles, practice & projects, Vol. II. Pergamon Press; 1993. p. 113–71 (Chapter 5). [43] Wen PH, Aliabadi MH. An improved meshless collocation method for elastostatic and elastodynamic problems. Commun Num Method Eng 2008;24:635–51. [44] Schanz M, Cheng AHD. Transient wave propagation in a one-dimensional poroelastic column. Acta Mech 2000;145:1–18. [45] Kim YK, Kingsbury HB. Dynamic characterization of poroelastic materials. Exp Mech 1979;19:252–8. [46] Lai X, Cai M, Ren F, Xie M, Esaki T. Assessment of rock mass characteristics and the excavation disturbed zone in the Lingxin Coal Mine beneath and Xitian river, China. Int J Rock Mech Min 2006;43:572–81. [47] Kwon S, Lee CS, Cho SJ, Jeon SW, Cho WJ. An investigation of the excavation damaged zone at the KAERI underground research tunnel. Tunn Undergr Space Technol 2009;24:1–13. [48] Sladek V, Sladek J, Zhang C. Analytical integrations in meshless implementations of local integral equations. In: Proceedings of the 8th world congress on computational mechanics – WCCM8, Int. Center for

[49]

[50]

[51] [52]

[53] [54] [55]

[56]

[57]

109

Numerical Methods in Engn. (CIMNE), Barcelona, 2008. CD-ROM, ISBN 97884-96736-55-9. Sladek V, Sladek J. Local integral equations implemented by MLSapproximation and analytical integrations. Eng Anal Bound Elem 2010;34:904–13. Sladek V, Sladek J, Zhang C. On increasing computational efficiency of local integral equation method combined with meshless implementations. CMES – Comput Model Eng Sci 2010;63:243–63. Soares D, Sladek V, Sladek J. Modified meshless local Petrov–Galerkin formulations for elastodynamics. Int J Num Method Eng 2012;90:1508–28. Soares Jr D. Dynamic analysis of porous media considering unequal phase discretization by meshless local Petrov–Galerkin formulations. CMES Comput Model Eng Sci 2010;61:177–200. Bergamaschi L. An efficient parallel MLPG method for poroelastic models. CMES Comput Model Eng Sci 2009;29:191–215. Wang JG, Liu GR, Lin P. Numerical analysis of Biot’s consolidation process by radial point interpolation method. Int J Solids Struct 2002;39:1557–73. Wang WD, Wang JG, Wang ZL, Nogami T. An unequal-order radial interpolation meshless method for Biot’s consolidation theory. Comput Geotech 2007;34:61–70. Wang ZL, Li YC. Analysis of factors influencing the solution of the consolidation problem by using an element-free Galerkin method. Comput Geosci 2006;32:624–31. Wang JG, Xie H, Leung CF. A local boundary integral-based meshless method for Biot’s consolidation problem. Eng Anal Bound Elem 2009;33:35–42.