- Email: [email protected]
Dynamic stress intensity factor (Mode I) of a penny-shaped crack in an infinite poroelastic solid
International Journal of Engineering Science 40 (2002) 637–646 www.elsevier.com/locate/ijengsci
Dynamic stress intensity factor (Mode I) of a penny-shaped crack in an infinite poroelastic solid Bo Jin *, Zheng Zhong Key Laboratory of Solid Mechanics of MOE, Department of Engineering Mechanics and Technology, Institute of Structure Theory, Tongji University, Shanghai 200092, PR China Received 14 May 2001; received in revised form 20 August 2001; accepted 29 August 2001
Abstract This paper considers the transient stress intensity factor (Mode I) of a penny-shaped crack in an infinite poroelastic solid. The crack surfaces are impermeable. By virtue of the integral transform methods, the poroelastodynamic mixed boundary value problems is formulated as a set of dual integral equations, which, in turn, are reduced to a Fredholm integral equation of the second kind in the Laplace transform domain. Time domain solutions are obtained by inverting Laplace domain solutions using a numerical scheme. A parametric study is presented to illustrate the influence of poroelastic material parameters on the transient stress intensity. The results obtained reveal that the dynamic stress intensity factor of poroelastic medium is smaller than that of elastic medium and the poroelastic medium with a small value of the potential of diffusivity shows higher value of the dynamic stress intensity factor. Ó 2002 Elsevier Science Ltd. All rights reserved. Keywords: Poroelastic solid; Penny-shaped crack; Transient stress intensity factor for Mode I problem
1. Introduction Within the framework of elastodynamic fracture mechanics, a considerable body of literature has been devoted to the study of dynamic stress intensity factor of stationary cracks. Such a problem is encountered when rapidly varying loads are applied to a nominally elastic body that contains stress concentrators such as narrow cuts or slits and, therefore, the computation of the field requires that inertia effects be taken into account. A typical elastodynamic analysis usually aims at determining the crack tip stress intensity factor as a function of time and loading/ge-
*
Corresponding author. Tel.: +86-21-6501-7470. E-mail address: [email protected] (B. Jin).
0020-7225/02/$ - see front matter Ó 2002 Elsevier Science Ltd. All rights reserved. PII: S 0 0 2 0 - 7 2 2 5 ( 0 1 ) 0 0 0 9 1 - X
638
B. Jin, Z. Zhong / International Journal of Engineering Science 40 (2002) 637–646
ometry/material parameters. Of course, such information is essential for a better understanding of the response and fracture behavior of dynamically loaded solids (see [1–4]). In contrast to the intensive investigation into the above elastodynamic problem and the multitude of relative papers, only a few studies exist dealing with the stress and pore pressure fields induced by shear cracks propagating steadily, and quasi-statically, in a poroelastic solid [5–10]. Rice and Simons [5], Simons [6], and Rudnicki [7] have presented the theoretical solutions for the propagation of a steady moving semi-infinite shear crack in a poroelastic solid, which are very useful to study faulting in the earth’s crust. To our knowledge, no previous study considered the transient response of a penny-shaped crack in a poroelastic solid. Notice that in this class of problems both inertial and poroelastic effects should be taken into account according to the poroelastodynamic theory of Biot [11]. The purpose of the present work is therefore to investigate a basic problem lying within the framework of poroelastodynamic fracture mechanics. This is the axisymmetric problem of an infinite poroelastic solid containing a penny-shaped crack, the faces of which are subjected to the transient action of equal and opposite normal tractions. In practical engineering applications, this problem can be considered as an idealization for the case of crack in rock which are known to be subjected to transient loadings. According to the authors’ knowledge, the solutions for the abovementioned problem are not reported in the literature. To obtain transient stress intensity factor, the poroelastodynamics mixed boundary value problem is solved in the Laplace transform domain by using integral transform techniques. Thereafter, Laplace transform solutions are converted to time domain solutions by using a numerical scheme. Selected time domain solutions are presented and compared with the corresponding ideal elastic solutions to demonstrate the influence of poroelastic material parameters.
2. Basic equations At this stage, it is convenient to non-dimensionalize all quantities with respect to length and pffiffiffiffiffiffiffiffi time by selecting the radius of the crack ‘‘a’’ as a unit of length and a q=l as a unit of time. Biot’s equations can be written for axisymmetric case [12] in terms of dimensionless variables as 1 oe of o2 2 ð1Þ r 2 ur þ ðk þ a2 M þ 1Þ aM ¼ 2 ður þ q wr Þ; r or or ot r2 uz þ ðk þ a2 M þ 1Þ
oe of o2 aM ¼ 2 ðuz þ q wz Þ; oz oz ot
ð2Þ
aM
oe of o2 owr M ¼ 2 ðq ur þ m wr Þ þ b ; or or ot ot
ð3Þ
aM
oe of o2 owz M ¼ 2 ðq uz þ m wz Þ þ b oz oz ot ot
ð4Þ
in which ur , uz and wr , wz are the dimensionless solid displacements and the dimensionless fluid displacements relative to matrix; t is the dimensionless time; a, M ¼ M=l, k ¼ k=l, q ¼ qf =q,
B. Jin, Z. Zhong / International Journal of Engineering Science 40 (2002) 637–646
639
pffiffiffiffiffiffi m ¼ m=q and b ¼ ab= ql are the dimensionless material parameters, where k and l are the Lame constants of the bulk material; a and M are Biot’s parameters accounting for compressibility of the two-phased material; q and qf are the mass densities of the bulk material and the pore fluid, respectively; m is a density-like parameter that depends on qf and the geometry of the pores; b is a parameter accounting for the internal friction due to the relative motion between the solid matrix and the pore fluid. The parameter b is equal to the ratio between the fluid viscosity and the intrinsic permeability of the medium. If internal friction is neglected, then b ¼ 0. our ur ouz þ þ ; e¼ or r oz
owr wr owz f¼ þ þ ; or r oz
r2 ¼
o2 1 o o2 þ þ : or2 r or oz2
The constitutive equations of poroelastic medium are given by our ouz srz ¼ þ ; oz or
ð5aÞ
ouz ur our þk rz ¼ ðk þ 2Þ þ ap; oz or r
ð5bÞ
p ¼ ae þ Mf;
ð5cÞ
where rz and srz are the dimensionless total stress components, and p is the dimensionless pore fluid pressure.
3. Laplace transform solutions The Laplace transform of function f ðr; tÞ with respect to time variable t denoted by f~ðr; sÞ is defined by Z 1 ~ f ðr; sÞ ¼ f ðr; tÞest dt; ð6Þ 0
s f~ðr; sÞ ¼ n
Z
1 0
on f ðr; tÞ st e dt otn
ð7Þ
and the inversion by 1 f ðr; tÞ ¼ 2pi
Z
cþi1
f~ðr; sÞest ds;
ð8Þ
ci1
where s is the Laplace transform parameter and the initial state of complete quietude is incorporated in Eq. (7). The parameter c in the inversion formula given pffiffiffiffiffiffiby ffi Eq. (8) is selected such that ~ the line ReðsÞ ¼ c is to the right of all singularities of f and i ¼ 1. Applying Laplace transform
640
B. Jin, Z. Zhong / International Journal of Engineering Science 40 (2002) 637–646
to Eqs. (1)–(4) and after some manipulations Biot’s equations can be reduced to the following equations in the Laplace transform domain:
1 o~ e o~ p ur ða #Þ ¼ 0; r 2 u~r þ ðk þ 1Þ s2 ð1 q #Þ~ r or or 2
r2 u~z þ ðk þ 1Þ r2 p~
o~ e o~ p uz ða #Þ ¼ 0; s2 ð1 q #Þ~ oz oz
ð9Þ ð10Þ
q s 2 q s2 ða #Þ ~ e~ ¼ 0; p # M #
ð11Þ
where # ¼ q s2 =ðm s2 þ b sÞ. In view of Eq. (4), the fluid displacement relative to matrix in the vertical direction can be expressed as 1 ~ z ðr; z; sÞ ¼ w b s þ m s2
! o~ p 2 þ q s u~z : oz
ð12Þ
For axisymmetric problems under consideration, it is natural to introduce Hankel integral transforms with respect to the radial coordinate as 0
f~ ðn; z; sÞ ¼
Z
1
rf~ðr; z; sÞJ0 ðnrÞ dr;
0
1
f~ ðn; z; sÞ ¼
Z
1
rf~ðr; z; sÞJ1 ðnrÞ dr:
0
Making use of the method suggested by the article [13], we obtain 0
p~ ¼ Aec1 z þ Bec2 z ; 0
u~z ¼ c1 a1 Aec1 z þ c2 a2 Bec2 z þ Cec3 z ; 1
n~ ur ¼ Aðv1 þ c21 a1 Þec1 z þ Bðv2 þ c22 a2 Þec2 z þ Cc3 ec3 z ; 0
1 2 c1 z 2 c2 z 2 c3 z c ð1 q s a ÞAe þ c ð1 q s a ÞBe q s Ce ; 1 2 1 2 b s þ m s2 c1 v1 þ ðc21 þ n2 Þa1 c2 v2 þ ðc22 þ n2 Þa2 ðc2 þ n2 Þ c3 z c1 z ¼ Ce ; Ae Bec2 z 3 n n n
~z ¼ w s~rz
1
0
r~z ¼ ðk 2c21 a1 aÞAec1 z þ ðk 2c22 a2 aÞBec2 z 2c3 Cec3 z ;
ð13Þ ð14Þ ð15Þ ð16Þ ð17Þ ð18Þ
where ci ¼
qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi n2 þ L2i ;
i ¼ 1; 2;
ð19Þ
B. Jin, Z. Zhong / International Journal of Engineering Science 40 (2002) 637–646
qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi c3 ¼ n2 þ K2 ;
641
ð20Þ
vi ¼
#M L2i þ q s2 ; q s2 ð# aÞM
ai ¼
k vi þ vi a þ # ; L2i K2
i ¼ 1; 2; i ¼ 1; 2:
ð21Þ ð22Þ
In addition L21 ¼
L22 ¼ b1 ¼
qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi b1 þ b21 4b2 2 qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi b1 b21 4b2 2
;
ð23Þ
;
ð24Þ
ðm s2 þ b sÞðk þ a2 M þ 2Þ þ M s2 2aM q s2 ; ðk þ 2ÞM
ð25Þ
m s4 þ b s3 ðq Þ2 s4 b2 ¼ ; ðk þ 2ÞM
ð26Þ
K ¼ ð1 q #Þ1=2 s;
ð27Þ
where L1 , L2 and K are the dimensionless wave numbers associated with the dilatational wave of the first kind (fast wave), the dilatational wave of the second kind (slow wave) and the rotational wave, respectively. 4. Governing integral equation Let the faces of the crack, described by z ¼ 0, 0 6 r 6 1, be subject to equal and opposite normal tractions r0 ðtÞ which is assumed to be uniformly distributed across the faces of the crack. The faces of the crack (z ¼ 0, 0 6 r 6 1) are impermeable. Because of symmetry, the problem of the infinite poroelastic solid (1 < z < 1, 0 6 r < 1) can be reduced to that of the half space (0 < z < 1, 0 6 r < 1) with the conditions, on the plane z ¼ 0, srz ðr; 0; tÞ ¼ 0;
0 6 r < 1; t > 0;
ð28Þ
wz ðr; 0; tÞ ¼ 0;
0 6 r < 1; t > 0;
ð29Þ
uz ðr; 0; tÞ ¼ 0;
1 < r < 1; t > 0;
ð30Þ
rz ðr; 0; tÞ ¼ r0 ðtÞ;
0 6 r 6 1; t > 0:
ð31Þ
642
B. Jin, Z. Zhong / International Journal of Engineering Science 40 (2002) 637–646
In view of Eqs. (28)–(31), we obtain the following dual integral equations in the Laplace transform domain: Z 1 D1 0 n u~z ðn; 0; sÞJ0 ðnrÞ dn ¼ r~0 ðsÞ; 0 6 r 6 1; ð32Þ D 0 Z 1 0 n~ ð33Þ uz ðn; 0; sÞJ0 ðnrÞ dn ¼ 0; 1 < r < 1; 0
where D1 ¼ ½g4 c1 ð1 q s2 a1 Þ g3 c2 ð1 q s2 a2 Þ ðc23 þ n2 Þ þ c1 c2 ½ðg1 g2 Þ q s2 ða1 g2 a2 g1 Þ ; ð34Þ D ¼ c1 c2 ½ðg2 g1 Þ þ ða2 a1 Þðc23 þ n2 Þ ;
ð35Þ
g1 ¼ v1 þ ðc21 þ n2 Þa1 ;
g2 ¼ v2 þ ðc22 þ n2 Þa2 ;
ð36Þ
g3 ¼ k v1 2c21 a1 a;
g4 ¼ k v2 2c22 a2 a:
ð37Þ
Using now the following identities: Z x rJ0 ðnrÞ sinðnxÞ pffiffiffiffiffiffiffiffiffiffiffiffiffiffi dr ¼ ; n x2 r2 0 Z 1 rJ0 ðnrÞ cosðnxÞ pffiffiffiffiffiffiffiffiffiffiffiffiffiffi dr ¼ ; 2 2 n r x x the following equations can be obtained: Z 1 D1 Lðn; sÞ sinðnxÞ dn ¼ h1 ðx; sÞ; 0 6 x 6 1; nD 0 Z 1 Lðn; sÞ sinðnxÞ dn ¼ 0; 1 < x < 1;
ð38Þ ð39Þ
0
where h1 ðx; sÞ ¼
Z
x
0 0
r~ r0 ðsÞ pffiffiffiffiffiffiffiffiffiffiffiffiffiffi dr ¼ x~ r0 ðsÞ; x2 r2
uz ðn; 0; sÞ: Lðn; sÞ ¼ n~
ð40Þ ð41Þ
It is found that lim
n!1
D1 ¼ l1 ; nD
ð42Þ
B. Jin, Z. Zhong / International Journal of Engineering Science 40 (2002) 637–646
643
where l1 is a constant. Then we define the following integral representation: 2~ r0 ðsÞ Lðn; sÞ ¼ pl1 that is Z 1
Z
1
/~ðy; sÞ sinðnyÞ dy
ð43Þ
0
Lðn; sÞ sinðnyÞ dy ¼
0
r~0 ðsÞ ~ /ðy; sÞH ð1 yÞ; l1
ð44Þ
where H ð1 yÞ is the Heaviside unit step function. Eq. (39) is automatically satisfied from Eq. (44). Substituting Eqs. (43) and (44) to Eq. (38), we have /~ðx; sÞ þ
Z
1
Kðx; y; sÞ/~ðy; sÞ dy ¼ x;
0 6 x 6 1;
ð45Þ
0
where 2 Kðx; y; sÞ ¼ p
Z 0
1
D1 1 sinðnxÞ sinðnyÞ dn: l1 nD
ð46Þ
The stress intensity factor in the Laplace transform domain K~I can be defined by K~I ðsÞ ¼ limþ r!1
pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2~ r0 ðsÞ 2pðr 1Þr~z ðr; 0; sÞ ¼ pffiffiffi /~1 ð1; sÞ: p
ð47Þ
5. Numerical results The Fredholm integral equation (45) was solved numerically for different values of s. For this purpose, Eq. (45) was reduced to systems of algebraic equations in the interval ½0; 1 by means of Simpson’s rule (50 equal subintervals were considered to discretize Eq. (45) in the range ½0; 1 ). The kernel was evaluated by using the extended trapezoidal rule. After obtaining the solutions in the Laplace transform domain, we then use a numerical Laplace inversion scheme to obtain timedomain solutions. There are several numerical Laplace inversion schemes reported in the literature. The present study uses the scheme proposed by Miller and Guy [14]. The latter technique has extensively been utilized in transient crack problems [15]. Variations of the dynamic stress intensity factor KI ðtÞ with the dimensionless time t are shown in Figs. 1 and 2 for the following loadings: (1) r0 ðtÞ ¼ f0 H ðtÞ; (2) r0 ðtÞ ¼ f0 ½H ðtÞ Hðt 2Þ , where H ðtÞ is the Heaviside unit step function in time. The dimensionless parameters of the poroelastic materials are k ¼ 1:2, M ¼ 8, q ¼ 0:53, m ¼ 1:1, and a ¼ 0:95, with b ¼ 0:1 and 10. In the figures ‘‘j’’ and ‘‘*’’ represent b ¼ 0:1 and 10, respectively, and ‘‘’’ represents the solution for an elastic medium. The Lame constants and mass density of the elastic medium equal the Lame constants and mass density of the solid matrix of the poroelastic medium. It is
644
B. Jin, Z. Zhong / International Journal of Engineering Science 40 (2002) 637–646
Fig. 1. Variations of KI ðtÞ with the dimensionless time for r0 ðtÞ ¼ f0 H ðtÞ.
Fig. 2. Variations of KI ðtÞ with the dimensionless time for r0 ðtÞ ¼ f0 ½H ðtÞ H ðt 2Þ .
anticipated that among the different poroelastic properties, b should have the most influence on the dynamic stress intensity factor because it reflects the potential of diffusivity of the material (permeability). The larger the permeability, the smaller the b is. It is noted from Figs. 1 and 2 that some differences exist between the dynamic stress intensity factor for poroelastic medium and for elastic medium. In general, the dynamic stress intensity factor of poroelastic medium is smaller than that of elastic medium. The decrease of the stress intensity factor of poroelastic medium means that coupling between deformation and fluid diffusion reduces the value of the stress intensity factor. It is also found that the poroelastic medium with a larger value of b shows higher value of the dynamic stress intensity factor. This effect of the parameter b is due to the change of the resistance to crack extension by the pore pressure. A high b means a more impermeable medium, so pore pressure is larger in the medium with a higher b than in the medium with a lower b . The resistance to crack extension is increased by a decrease in pore pressure. More specifically, the resistance to crack extension depends on the effective compressive stress, the total stress minus the pore fluid pressure. Thus, an increase in pore pressure increases the stress
B. Jin, Z. Zhong / International Journal of Engineering Science 40 (2002) 637–646
645
intensity factor and a decrease decreases the stress intensity factor. Fig. 1 indicates that KI ðtÞ tends to rise quickly as time increases. The curves reach a peak and then decrease in magnitude. Fig. 2 indicates variations of KI ðtÞ become more oscillatory for the loading r0 ðtÞ ¼ f0 ½H ðtÞ H ðt 2Þ when compared with Fig. 1.
6. Conclusions In this paper, a mathematical formulation is presented for the transient stress intensity factor of a penny-shaped crack in an infinite poroelastic solid. With the aid of Laplace transforms and Hankel transforms, the mixed boundary value problem is formulated as dual integral equations in the Laplace transform domain. By appropriate transforms, it is shown that the dual integral equations can be reduced to a Fredholm integral equation of the second kind. The final Laplace transform inversion is done numerically. Our results give the history of the stress intensity factor and show the effect of the poroelastic material parameters upon the crack-tip stress intensity factor.
Acknowledgements The helpful advice and valuable comments from the reviewer are gratefully acknowledged. The reported work was performed with the financial supported by the Teaching and Research Award Fund for Outstanding Young Teachers in Higher Education Institutions of MOE, PR China.
References [1] [2] [3] [4] [5] [6] [7] [8] [9] [10] [11] [12]
J.D. Achenbach, Wave Propagation in Elastic Solid, North-Holland, Amsterdam, 1973. A.C. Eringen, E.S. Suhubi, Elastodynamics, vol. 2, Linear Theory, Academic Press, New York, 1975. E.P. Chen, G.C. Sih, in: G.C. Sih (Ed.), Elastodynamic Crack Problems, Noordhoff, Leyden, 1977. L.B. Freund, Dynamic Fracture Mechanics, Cambridge University Press, Cambridge, 1990. J.R. Rice, D.A. Simons, The stabilization of spreading shear faults by coupled deformation–diffusion effects in fluid-infiltrated porous materials, J. Geophys. Res. 81 (1976) 5322–5344. D.A. Simons, Boundary-layer analysis of propagating Mode-II cracks in porous elastic solids, J. Mech. and Phys. of Solids 25 (1977) 99–116. J.W. Rudnicki, Boundary layer analysis of plane strain shear cracks propagating steadily on an impermeable plane in an elastic diffusive solid, J. Mech. Phys. Solids 39 (2) (1991) 201–221. J.W. Rudnicki, D.A. Koutsibelas, Steady propagation of plane strain shear cracks on an impermeable plane in an elastic diffusive solid, Internat. J. Solids and Structures 27 (2) (1991) 205–225. R.V. Craster, C. Atkinson, Shear cracks in thermoelastic and poroelastic media, J. Mech. Phys. Solids 40 (4) (1992) 887–924. J.W. Rudnicki, Geomechanics, Internat. J. Solids and Structures 37 (1) (2000) 349–358. M.A. Biot, Mechanics of deformation and acoustic propagation in porous media, J. Appl. Phys. 33 (4) (1962) 1482–1498. B. Jin, The vertical vibration of an elastic circular plate on a fluid-saturated porous half space, Internat. J. Engrg. Sci. 37 (3) (1999) 379–393.
646
B. Jin, Z. Zhong / International Journal of Engineering Science 40 (2002) 637–646
[13] B. Jin, H. Liu, Vertical dynamic response of a disk on a saturated poroelastic half space, Soil Dynamics and Earthquake Engrg. 18 (6) (1999) 437–443. [14] M.K. Miller, W.T. Guy, Numerical inversion of the Laplace transform by use of Jacobi polynomials, SIAM J. Numer. Anal. 3 (1966) 624–635. [15] E.P. Chen, Elastodynamic response of a penny-shaped crack in a cylinder of finite radius, Internat. J. Engrg. Sci. 17 (1979) 379–385.
Dynamic stress intensity factor (Mode I) of a penny-shaped crack in an infinite poroelastic solid Bo Jin *, Zheng Zhong Key Laboratory of Solid Mechanics of MOE, Department of Engineering Mechanics and Technology, Institute of Structure Theory, Tongji University, Shanghai 200092, PR China Received 14 May 2001; received in revised form 20 August 2001; accepted 29 August 2001
Abstract This paper considers the transient stress intensity factor (Mode I) of a penny-shaped crack in an infinite poroelastic solid. The crack surfaces are impermeable. By virtue of the integral transform methods, the poroelastodynamic mixed boundary value problems is formulated as a set of dual integral equations, which, in turn, are reduced to a Fredholm integral equation of the second kind in the Laplace transform domain. Time domain solutions are obtained by inverting Laplace domain solutions using a numerical scheme. A parametric study is presented to illustrate the influence of poroelastic material parameters on the transient stress intensity. The results obtained reveal that the dynamic stress intensity factor of poroelastic medium is smaller than that of elastic medium and the poroelastic medium with a small value of the potential of diffusivity shows higher value of the dynamic stress intensity factor. Ó 2002 Elsevier Science Ltd. All rights reserved. Keywords: Poroelastic solid; Penny-shaped crack; Transient stress intensity factor for Mode I problem
1. Introduction Within the framework of elastodynamic fracture mechanics, a considerable body of literature has been devoted to the study of dynamic stress intensity factor of stationary cracks. Such a problem is encountered when rapidly varying loads are applied to a nominally elastic body that contains stress concentrators such as narrow cuts or slits and, therefore, the computation of the field requires that inertia effects be taken into account. A typical elastodynamic analysis usually aims at determining the crack tip stress intensity factor as a function of time and loading/ge-
*
Corresponding author. Tel.: +86-21-6501-7470. E-mail address: [email protected] (B. Jin).
0020-7225/02/$ - see front matter Ó 2002 Elsevier Science Ltd. All rights reserved. PII: S 0 0 2 0 - 7 2 2 5 ( 0 1 ) 0 0 0 9 1 - X
638
B. Jin, Z. Zhong / International Journal of Engineering Science 40 (2002) 637–646
ometry/material parameters. Of course, such information is essential for a better understanding of the response and fracture behavior of dynamically loaded solids (see [1–4]). In contrast to the intensive investigation into the above elastodynamic problem and the multitude of relative papers, only a few studies exist dealing with the stress and pore pressure fields induced by shear cracks propagating steadily, and quasi-statically, in a poroelastic solid [5–10]. Rice and Simons [5], Simons [6], and Rudnicki [7] have presented the theoretical solutions for the propagation of a steady moving semi-infinite shear crack in a poroelastic solid, which are very useful to study faulting in the earth’s crust. To our knowledge, no previous study considered the transient response of a penny-shaped crack in a poroelastic solid. Notice that in this class of problems both inertial and poroelastic effects should be taken into account according to the poroelastodynamic theory of Biot [11]. The purpose of the present work is therefore to investigate a basic problem lying within the framework of poroelastodynamic fracture mechanics. This is the axisymmetric problem of an infinite poroelastic solid containing a penny-shaped crack, the faces of which are subjected to the transient action of equal and opposite normal tractions. In practical engineering applications, this problem can be considered as an idealization for the case of crack in rock which are known to be subjected to transient loadings. According to the authors’ knowledge, the solutions for the abovementioned problem are not reported in the literature. To obtain transient stress intensity factor, the poroelastodynamics mixed boundary value problem is solved in the Laplace transform domain by using integral transform techniques. Thereafter, Laplace transform solutions are converted to time domain solutions by using a numerical scheme. Selected time domain solutions are presented and compared with the corresponding ideal elastic solutions to demonstrate the influence of poroelastic material parameters.
2. Basic equations At this stage, it is convenient to non-dimensionalize all quantities with respect to length and pffiffiffiffiffiffiffiffi time by selecting the radius of the crack ‘‘a’’ as a unit of length and a q=l as a unit of time. Biot’s equations can be written for axisymmetric case [12] in terms of dimensionless variables as 1 oe of o2 2 ð1Þ r 2 ur þ ðk þ a2 M þ 1Þ aM ¼ 2 ður þ q wr Þ; r or or ot r2 uz þ ðk þ a2 M þ 1Þ
oe of o2 aM ¼ 2 ðuz þ q wz Þ; oz oz ot
ð2Þ
aM
oe of o2 owr M ¼ 2 ðq ur þ m wr Þ þ b ; or or ot ot
ð3Þ
aM
oe of o2 owz M ¼ 2 ðq uz þ m wz Þ þ b oz oz ot ot
ð4Þ
in which ur , uz and wr , wz are the dimensionless solid displacements and the dimensionless fluid displacements relative to matrix; t is the dimensionless time; a, M ¼ M=l, k ¼ k=l, q ¼ qf =q,
B. Jin, Z. Zhong / International Journal of Engineering Science 40 (2002) 637–646
639
pffiffiffiffiffiffi m ¼ m=q and b ¼ ab= ql are the dimensionless material parameters, where k and l are the Lame constants of the bulk material; a and M are Biot’s parameters accounting for compressibility of the two-phased material; q and qf are the mass densities of the bulk material and the pore fluid, respectively; m is a density-like parameter that depends on qf and the geometry of the pores; b is a parameter accounting for the internal friction due to the relative motion between the solid matrix and the pore fluid. The parameter b is equal to the ratio between the fluid viscosity and the intrinsic permeability of the medium. If internal friction is neglected, then b ¼ 0. our ur ouz þ þ ; e¼ or r oz
owr wr owz f¼ þ þ ; or r oz
r2 ¼
o2 1 o o2 þ þ : or2 r or oz2
The constitutive equations of poroelastic medium are given by our ouz srz ¼ þ ; oz or
ð5aÞ
ouz ur our þk rz ¼ ðk þ 2Þ þ ap; oz or r
ð5bÞ
p ¼ ae þ Mf;
ð5cÞ
where rz and srz are the dimensionless total stress components, and p is the dimensionless pore fluid pressure.
3. Laplace transform solutions The Laplace transform of function f ðr; tÞ with respect to time variable t denoted by f~ðr; sÞ is defined by Z 1 ~ f ðr; sÞ ¼ f ðr; tÞest dt; ð6Þ 0
s f~ðr; sÞ ¼ n
Z
1 0
on f ðr; tÞ st e dt otn
ð7Þ
and the inversion by 1 f ðr; tÞ ¼ 2pi
Z
cþi1
f~ðr; sÞest ds;
ð8Þ
ci1
where s is the Laplace transform parameter and the initial state of complete quietude is incorporated in Eq. (7). The parameter c in the inversion formula given pffiffiffiffiffiffiby ffi Eq. (8) is selected such that ~ the line ReðsÞ ¼ c is to the right of all singularities of f and i ¼ 1. Applying Laplace transform
640
B. Jin, Z. Zhong / International Journal of Engineering Science 40 (2002) 637–646
to Eqs. (1)–(4) and after some manipulations Biot’s equations can be reduced to the following equations in the Laplace transform domain:
1 o~ e o~ p ur ða #Þ ¼ 0; r 2 u~r þ ðk þ 1Þ s2 ð1 q #Þ~ r or or 2
r2 u~z þ ðk þ 1Þ r2 p~
o~ e o~ p uz ða #Þ ¼ 0; s2 ð1 q #Þ~ oz oz
ð9Þ ð10Þ
q s 2 q s2 ða #Þ ~ e~ ¼ 0; p # M #
ð11Þ
where # ¼ q s2 =ðm s2 þ b sÞ. In view of Eq. (4), the fluid displacement relative to matrix in the vertical direction can be expressed as 1 ~ z ðr; z; sÞ ¼ w b s þ m s2
! o~ p 2 þ q s u~z : oz
ð12Þ
For axisymmetric problems under consideration, it is natural to introduce Hankel integral transforms with respect to the radial coordinate as 0
f~ ðn; z; sÞ ¼
Z
1
rf~ðr; z; sÞJ0 ðnrÞ dr;
0
1
f~ ðn; z; sÞ ¼
Z
1
rf~ðr; z; sÞJ1 ðnrÞ dr:
0
Making use of the method suggested by the article [13], we obtain 0
p~ ¼ Aec1 z þ Bec2 z ; 0
u~z ¼ c1 a1 Aec1 z þ c2 a2 Bec2 z þ Cec3 z ; 1
n~ ur ¼ Aðv1 þ c21 a1 Þec1 z þ Bðv2 þ c22 a2 Þec2 z þ Cc3 ec3 z ; 0
1 2 c1 z 2 c2 z 2 c3 z c ð1 q s a ÞAe þ c ð1 q s a ÞBe q s Ce ; 1 2 1 2 b s þ m s2 c1 v1 þ ðc21 þ n2 Þa1 c2 v2 þ ðc22 þ n2 Þa2 ðc2 þ n2 Þ c3 z c1 z ¼ Ce ; Ae Bec2 z 3 n n n
~z ¼ w s~rz
1
0
r~z ¼ ðk 2c21 a1 aÞAec1 z þ ðk 2c22 a2 aÞBec2 z 2c3 Cec3 z ;
ð13Þ ð14Þ ð15Þ ð16Þ ð17Þ ð18Þ
where ci ¼
qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi n2 þ L2i ;
i ¼ 1; 2;
ð19Þ
B. Jin, Z. Zhong / International Journal of Engineering Science 40 (2002) 637–646
qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi c3 ¼ n2 þ K2 ;
641
ð20Þ
vi ¼
#M L2i þ q s2 ; q s2 ð# aÞM
ai ¼
k vi þ vi a þ # ; L2i K2
i ¼ 1; 2; i ¼ 1; 2:
ð21Þ ð22Þ
In addition L21 ¼
L22 ¼ b1 ¼
qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi b1 þ b21 4b2 2 qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi b1 b21 4b2 2
;
ð23Þ
;
ð24Þ
ðm s2 þ b sÞðk þ a2 M þ 2Þ þ M s2 2aM q s2 ; ðk þ 2ÞM
ð25Þ
m s4 þ b s3 ðq Þ2 s4 b2 ¼ ; ðk þ 2ÞM
ð26Þ
K ¼ ð1 q #Þ1=2 s;
ð27Þ
where L1 , L2 and K are the dimensionless wave numbers associated with the dilatational wave of the first kind (fast wave), the dilatational wave of the second kind (slow wave) and the rotational wave, respectively. 4. Governing integral equation Let the faces of the crack, described by z ¼ 0, 0 6 r 6 1, be subject to equal and opposite normal tractions r0 ðtÞ which is assumed to be uniformly distributed across the faces of the crack. The faces of the crack (z ¼ 0, 0 6 r 6 1) are impermeable. Because of symmetry, the problem of the infinite poroelastic solid (1 < z < 1, 0 6 r < 1) can be reduced to that of the half space (0 < z < 1, 0 6 r < 1) with the conditions, on the plane z ¼ 0, srz ðr; 0; tÞ ¼ 0;
0 6 r < 1; t > 0;
ð28Þ
wz ðr; 0; tÞ ¼ 0;
0 6 r < 1; t > 0;
ð29Þ
uz ðr; 0; tÞ ¼ 0;
1 < r < 1; t > 0;
ð30Þ
rz ðr; 0; tÞ ¼ r0 ðtÞ;
0 6 r 6 1; t > 0:
ð31Þ
642
B. Jin, Z. Zhong / International Journal of Engineering Science 40 (2002) 637–646
In view of Eqs. (28)–(31), we obtain the following dual integral equations in the Laplace transform domain: Z 1 D1 0 n u~z ðn; 0; sÞJ0 ðnrÞ dn ¼ r~0 ðsÞ; 0 6 r 6 1; ð32Þ D 0 Z 1 0 n~ ð33Þ uz ðn; 0; sÞJ0 ðnrÞ dn ¼ 0; 1 < r < 1; 0
where D1 ¼ ½g4 c1 ð1 q s2 a1 Þ g3 c2 ð1 q s2 a2 Þ ðc23 þ n2 Þ þ c1 c2 ½ðg1 g2 Þ q s2 ða1 g2 a2 g1 Þ ; ð34Þ D ¼ c1 c2 ½ðg2 g1 Þ þ ða2 a1 Þðc23 þ n2 Þ ;
ð35Þ
g1 ¼ v1 þ ðc21 þ n2 Þa1 ;
g2 ¼ v2 þ ðc22 þ n2 Þa2 ;
ð36Þ
g3 ¼ k v1 2c21 a1 a;
g4 ¼ k v2 2c22 a2 a:
ð37Þ
Using now the following identities: Z x rJ0 ðnrÞ sinðnxÞ pffiffiffiffiffiffiffiffiffiffiffiffiffiffi dr ¼ ; n x2 r2 0 Z 1 rJ0 ðnrÞ cosðnxÞ pffiffiffiffiffiffiffiffiffiffiffiffiffiffi dr ¼ ; 2 2 n r x x the following equations can be obtained: Z 1 D1 Lðn; sÞ sinðnxÞ dn ¼ h1 ðx; sÞ; 0 6 x 6 1; nD 0 Z 1 Lðn; sÞ sinðnxÞ dn ¼ 0; 1 < x < 1;
ð38Þ ð39Þ
0
where h1 ðx; sÞ ¼
Z
x
0 0
r~ r0 ðsÞ pffiffiffiffiffiffiffiffiffiffiffiffiffiffi dr ¼ x~ r0 ðsÞ; x2 r2
uz ðn; 0; sÞ: Lðn; sÞ ¼ n~
ð40Þ ð41Þ
It is found that lim
n!1
D1 ¼ l1 ; nD
ð42Þ
B. Jin, Z. Zhong / International Journal of Engineering Science 40 (2002) 637–646
643
where l1 is a constant. Then we define the following integral representation: 2~ r0 ðsÞ Lðn; sÞ ¼ pl1 that is Z 1
Z
1
/~ðy; sÞ sinðnyÞ dy
ð43Þ
0
Lðn; sÞ sinðnyÞ dy ¼
0
r~0 ðsÞ ~ /ðy; sÞH ð1 yÞ; l1
ð44Þ
where H ð1 yÞ is the Heaviside unit step function. Eq. (39) is automatically satisfied from Eq. (44). Substituting Eqs. (43) and (44) to Eq. (38), we have /~ðx; sÞ þ
Z
1
Kðx; y; sÞ/~ðy; sÞ dy ¼ x;
0 6 x 6 1;
ð45Þ
0
where 2 Kðx; y; sÞ ¼ p
Z 0
1
D1 1 sinðnxÞ sinðnyÞ dn: l1 nD
ð46Þ
The stress intensity factor in the Laplace transform domain K~I can be defined by K~I ðsÞ ¼ limþ r!1
pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2~ r0 ðsÞ 2pðr 1Þr~z ðr; 0; sÞ ¼ pffiffiffi /~1 ð1; sÞ: p
ð47Þ
5. Numerical results The Fredholm integral equation (45) was solved numerically for different values of s. For this purpose, Eq. (45) was reduced to systems of algebraic equations in the interval ½0; 1 by means of Simpson’s rule (50 equal subintervals were considered to discretize Eq. (45) in the range ½0; 1 ). The kernel was evaluated by using the extended trapezoidal rule. After obtaining the solutions in the Laplace transform domain, we then use a numerical Laplace inversion scheme to obtain timedomain solutions. There are several numerical Laplace inversion schemes reported in the literature. The present study uses the scheme proposed by Miller and Guy [14]. The latter technique has extensively been utilized in transient crack problems [15]. Variations of the dynamic stress intensity factor KI ðtÞ with the dimensionless time t are shown in Figs. 1 and 2 for the following loadings: (1) r0 ðtÞ ¼ f0 H ðtÞ; (2) r0 ðtÞ ¼ f0 ½H ðtÞ Hðt 2Þ , where H ðtÞ is the Heaviside unit step function in time. The dimensionless parameters of the poroelastic materials are k ¼ 1:2, M ¼ 8, q ¼ 0:53, m ¼ 1:1, and a ¼ 0:95, with b ¼ 0:1 and 10. In the figures ‘‘j’’ and ‘‘*’’ represent b ¼ 0:1 and 10, respectively, and ‘‘’’ represents the solution for an elastic medium. The Lame constants and mass density of the elastic medium equal the Lame constants and mass density of the solid matrix of the poroelastic medium. It is
644
B. Jin, Z. Zhong / International Journal of Engineering Science 40 (2002) 637–646
Fig. 1. Variations of KI ðtÞ with the dimensionless time for r0 ðtÞ ¼ f0 H ðtÞ.
Fig. 2. Variations of KI ðtÞ with the dimensionless time for r0 ðtÞ ¼ f0 ½H ðtÞ H ðt 2Þ .
anticipated that among the different poroelastic properties, b should have the most influence on the dynamic stress intensity factor because it reflects the potential of diffusivity of the material (permeability). The larger the permeability, the smaller the b is. It is noted from Figs. 1 and 2 that some differences exist between the dynamic stress intensity factor for poroelastic medium and for elastic medium. In general, the dynamic stress intensity factor of poroelastic medium is smaller than that of elastic medium. The decrease of the stress intensity factor of poroelastic medium means that coupling between deformation and fluid diffusion reduces the value of the stress intensity factor. It is also found that the poroelastic medium with a larger value of b shows higher value of the dynamic stress intensity factor. This effect of the parameter b is due to the change of the resistance to crack extension by the pore pressure. A high b means a more impermeable medium, so pore pressure is larger in the medium with a higher b than in the medium with a lower b . The resistance to crack extension is increased by a decrease in pore pressure. More specifically, the resistance to crack extension depends on the effective compressive stress, the total stress minus the pore fluid pressure. Thus, an increase in pore pressure increases the stress
B. Jin, Z. Zhong / International Journal of Engineering Science 40 (2002) 637–646
645
intensity factor and a decrease decreases the stress intensity factor. Fig. 1 indicates that KI ðtÞ tends to rise quickly as time increases. The curves reach a peak and then decrease in magnitude. Fig. 2 indicates variations of KI ðtÞ become more oscillatory for the loading r0 ðtÞ ¼ f0 ½H ðtÞ H ðt 2Þ when compared with Fig. 1.
6. Conclusions In this paper, a mathematical formulation is presented for the transient stress intensity factor of a penny-shaped crack in an infinite poroelastic solid. With the aid of Laplace transforms and Hankel transforms, the mixed boundary value problem is formulated as dual integral equations in the Laplace transform domain. By appropriate transforms, it is shown that the dual integral equations can be reduced to a Fredholm integral equation of the second kind. The final Laplace transform inversion is done numerically. Our results give the history of the stress intensity factor and show the effect of the poroelastic material parameters upon the crack-tip stress intensity factor.
Acknowledgements The helpful advice and valuable comments from the reviewer are gratefully acknowledged. The reported work was performed with the financial supported by the Teaching and Research Award Fund for Outstanding Young Teachers in Higher Education Institutions of MOE, PR China.
References [1] [2] [3] [4] [5] [6] [7] [8] [9] [10] [11] [12]
J.D. Achenbach, Wave Propagation in Elastic Solid, North-Holland, Amsterdam, 1973. A.C. Eringen, E.S. Suhubi, Elastodynamics, vol. 2, Linear Theory, Academic Press, New York, 1975. E.P. Chen, G.C. Sih, in: G.C. Sih (Ed.), Elastodynamic Crack Problems, Noordhoff, Leyden, 1977. L.B. Freund, Dynamic Fracture Mechanics, Cambridge University Press, Cambridge, 1990. J.R. Rice, D.A. Simons, The stabilization of spreading shear faults by coupled deformation–diffusion effects in fluid-infiltrated porous materials, J. Geophys. Res. 81 (1976) 5322–5344. D.A. Simons, Boundary-layer analysis of propagating Mode-II cracks in porous elastic solids, J. Mech. and Phys. of Solids 25 (1977) 99–116. J.W. Rudnicki, Boundary layer analysis of plane strain shear cracks propagating steadily on an impermeable plane in an elastic diffusive solid, J. Mech. Phys. Solids 39 (2) (1991) 201–221. J.W. Rudnicki, D.A. Koutsibelas, Steady propagation of plane strain shear cracks on an impermeable plane in an elastic diffusive solid, Internat. J. Solids and Structures 27 (2) (1991) 205–225. R.V. Craster, C. Atkinson, Shear cracks in thermoelastic and poroelastic media, J. Mech. Phys. Solids 40 (4) (1992) 887–924. J.W. Rudnicki, Geomechanics, Internat. J. Solids and Structures 37 (1) (2000) 349–358. M.A. Biot, Mechanics of deformation and acoustic propagation in porous media, J. Appl. Phys. 33 (4) (1962) 1482–1498. B. Jin, The vertical vibration of an elastic circular plate on a fluid-saturated porous half space, Internat. J. Engrg. Sci. 37 (3) (1999) 379–393.
646
B. Jin, Z. Zhong / International Journal of Engineering Science 40 (2002) 637–646
[13] B. Jin, H. Liu, Vertical dynamic response of a disk on a saturated poroelastic half space, Soil Dynamics and Earthquake Engrg. 18 (6) (1999) 437–443. [14] M.K. Miller, W.T. Guy, Numerical inversion of the Laplace transform by use of Jacobi polynomials, SIAM J. Numer. Anal. 3 (1966) 624–635. [15] E.P. Chen, Elastodynamic response of a penny-shaped crack in a cylinder of finite radius, Internat. J. Engrg. Sci. 17 (1979) 379–385.