Effects of permeability conditions on time-dependent fracture of poroelastic media

Mechanics of Materials 138 (2019) 103156

Contents lists available at ScienceDirect

Mechanics of Materials journal homepage: www.elsevier.com/locate/mechmat

Research paper

Effects of permeability conditions on time-dependent fracture of poroelastic media

T

Yu-Yun Lin , Chen-Hsueh Yang ⁎

Department of Civil Engineering, National Cheng Kung University, Taiwan

ARTICLE INFO

ABSTRACT

Keywords: Poroelastic materials Instantaneous energy release rate Cohesive zone elements

The time-dependent fracture of fluid-infiltrated porous materials is known to be affected by fluid diffusion, which is strongly influenced by the initial pore pressure field attributed to loading and the permeability conditions on the boundaries. In this research, the concurrent solid deformation and fluid migration for a tractionfree crack of finite length 2a in a poroelastic medium are analyzed. The remote boundaries of the medium are subjected to constant stress or constant stain in the direction perpendicular to crack faces. Without confining the range of the fluid diffusion to the neighborhood of the crack, four types of permeability conditions are considered on the boundaries: (I) the remote boundaries and the crack faces are both impermeable, (II) the remote boundaries are permeable to constant pore pressure, but the crack faces are impermeable, (III) the remote boundaries are impermeable, but the crack faces are permeable to constant pore pressure, and (IV) the remote boundaries and the crack faces are both permeable to constant pore pressure. The singularities of the stress and the pore pressure gradient at small flow times are obtained from finite element simulations and compared to asymptotic predictions. The instantaneous fracture energy J is evaluated using a method that uses the J-integral around the cohesive zone embedded ahead of the crack tip in the simulations. The numerical results of J are validated with the asymptotes predicted from stress intensities at small flow times. When there is a difference between the imposed pore pressure and the initial pore pressure caused by applying loading on the boundaries, the variation in J(t) is affected by short and long ranges of fluid flow, and J(t) reaches its drained limit at a time on the order of W2/c, where W is the dimensional size and c is the consolidation coefficient of porous materials. The drained limit may not be the maximum value of J, especially when the difference on the remote boundaries is positive under strain control. For stress-control cases, whether the maximum J occurs before the drained state depends not only on the difference in pore pressure but also on the location of the permeable boundaries.

1. Introduction The mechanical properties of polymeric gels are important for many applications, such as drug delivery, tissue engineering, biosensors, and micro-devices (Duncan, 2003; Beebe et al., 2000; Luo and Shoichet, 2004; Nowak et al., 2002). Fracture properties are particularly important for the development of strong and tough polymeric gels. Several studies have noted that the rate dependence of the fracture toughness of polymeric gels implies kinetic processes associated with fracture. Experiments have shown that the fracture of some types of polymeric gel can be delayed (Bonn et al., 1998; Skrzeszeswka et al., 2010). Pizzocolo et al. (2013) demonstrated that a crack in a polymeric gel does not propagate smoothly but has a step-wise pattern; the duration of the pauses in propagation is inversely related to permeability, suggesting a flow transport process around the crack tip.

⁎

The effects of fluid transport on fracture have been studied by many researchers within the framework of poroelasticity developed by Biot (1941). There are numerous studies on crack propagation caused by hydraulic pressure in poroelastic media (Atkinson and Craster, 1991; Boone and Ingraffea, 1990; Cheng et al., 1994) because hydraulic fracturing is commonly used to initiate and propagate cracks in rocks for oil extraction. For a given set of stresses exerted on a crack in a poroelastic material, Rice and Cleary (1976) stated that the amount of energy flow to the crack tip per unit crack advance depends on whether the surrounding material responds in an undrained or drained manner. For a crack in a saturated porous medium subjected to loads, the highly concentrated hydrostatic tension near the crack tip lowers the local fluid pore pressure and causes a converging fluid flow to the crack tip, which acts like a sink. Consequently, the fluid flow effectively increases the stress intensity factor and fracture energy with time. Atkinson and

Corresponding author. E-mail address: [email protected] (Y.-Y. Lin).

https://doi.org/10.1016/j.mechmat.2019.103156 Received 27 April 2018; Received in revised form 2 August 2019; Accepted 27 August 2019 Available online 30 August 2019 0167-6636/ © 2019 Elsevier Ltd. All rights reserved.

Mechanics of Materials 138 (2019) 103156

Y.-Y. Lin and C.-H. Yang

Craster (1991) found the stress intensity factor for a semi-infinite crack subjected to an exponentially decayed normal stress distribution on the crack faces in the Laplace- and Fourier-transformed domain using the Wiener–Hopf technique. Later, Craster and Atkinson (1994) expanded on the previous results at small times and obtained the asymptotic behaviors of the stress singularity for semi-infinite cracks. They also developed an alternative asymptotic approach using the integral equation of dislocation density and pore pressure gradient discontinuities. Atkinson and Craster (1992) derived an invariant integral in the Laplace-transformed domain of poroelasticity and applied it to connect the inner singular fields and the outer fields. Some important results of their studies were summarized in Craster and Atkinson (1996a). Wang and Hong (2012) numerically studied the delayed fracture of largely deformed hydrogels with the assumption that viscoelastic relaxation occurs much faster than the characteristic time of flow-dependent behavior. Bouklas et al. (2015) computed the transient energy release rate of hydrogels using a J-integral formula (the same as that used by Wilson and Yu (1979) for thermoelastic cracks) considering the effects of large deformation and solvent diffusion but ignoring those of viscoelasticity. For steady crack propagation in a polymer gel, Noselli et al. (2016) performed an asymptotic analysis of near-tip fields. Their results are consistent with the findings of Atkinson and Craster (1991). They expressed the energy release rate required in the far field as the sum of the fracture energy needed for chain breaking at the crack tip and the energy dissipation computed using an area integral of the singular pore pressure gradient over the domain encircled by the far field. This expression is valid for both stationary and propagating cracks as long as the small-scale process zone holds. The calculation of the area integral relies on the pore pressure and dilation fields in the process zone. However, these fields for stationary and propagating cracks are different, as indicated by Atkinson and Craster (1991). Moreover, the result of the area integral for a steadily propagating crack is semi-analytical because the size of the process zone cannot be obtained from the asymptotic analysis. They also employed a poroelastic cohesive zone model in finite element method (FEM) simulations to investigate the effect of failure mechanisms on the toughness ratio (i.e., the ratio of the stress intensity factor at infinity to its counterpart at the crack tip). Yang and Lin (2018) developed a numerical method for calculating the instantaneous fracture energy of a mode-I crack of length 2a in an infinite poroviscoelastic medium. This method computes the instantaneous fracture energy using the J-integral around the cohesive zone ahead of the crack tip, and thus the driving energy for crack extension, along with additional energy dissipation through fluid diffusion and viscoelasticity, can be evaluated directly. Unlike the integrals used by Bouklas et al. (2015) and Noselli et al. (2016), the J-integral used by Yang and Lin (2018) needs no derivation for the energy dissipation in the area enclosed by the contour. To obtain an accurate instantaneous energy release rate J, the cohesive zone length must be less than the size of the region dominated by the inner stress intensity factor at all times. This region is particularly small for permeable cracks at the undrained state. The numerical results were verified with the undrained and drained limits for impermeable and permeable cracks in poroelastic media. The limiting values of J caused by drainage or relaxation in poroviscoelastic media were compared with the numerical results for permeable and impermeable cracks. In the study of Yang and Lin (2018), a specific biaxial stress state was applied to the impermeable remote boundaries to maintain zero pore pressure at locations far from the crack tip, so that fluids diffuse mainly to eliminate the initial singularity near the crack tip. As a result, the diffusion is confined to the neighborhood of the crack, and the drainage time is on the order of a2/c, where a is the half crack length and c is the consolidation coefficient of porous materials. When the viscoelastic relaxation time differs from the drainage time in poroviscoelastic media, the J curve exhibits two-stage behavior. Otherwise,

the viscoelastic and the diffusion phenomena are coupled and difficult to distinguish in the J curve. However, when stresses and permeability at infinity vary, fluid diffusion may not be confined to the neighborhood of the crack. The aim of this research is to understand how permeability conditions on the boundaries affect the time-dependent behavior of fracture energy for cracks in poroelastic media subjected to mode-I loadings. In general, fluid diffusion is initially driven by the inhomogeneous stress and pore pressure field near the crack tip as well as the permeability on the crack faces, and later influenced by the condition of permeability applied to the remote boundaries. When the characteristic distance of fluid diffusion varies, drainage time changes correspondingly, affecting the time-dependent behavior of fracture energy for cracks. Mixed permeability conditions on the boundaries of poroelastic media and on the crack faces are considered in the following section. Normal stress and pore pressure fields are analyzed using FEM simulations to find singular behaviors in the near-tip region. The singularities of the stress and pore pressure gradient at small flow times are compared with the predictions from semi-infinite cracks. Furthermore, the time-dependent fracture energy is computed using FEM simulations with embedded cohesive elements. The numerical results are verified with the predictions calculated from the stress intensities at small flow times and at the drained state. 2. Constitutive model and field equations of poroelastic materials In a poroelastic material, fluids infiltrate the isotropic linear elastic solid frame. According to Biot (1941, 1962), the constitutive equations can be expressed as ij

= 2µ

p=

ij

+

c e ij

C

(1a)

ij

(1b)

Ce + M

where σij denotes the total stress, p denotes the increment of fluid pressure, ɛij is the solid-frame strain, ξ is the fluid volume increase per unit volume of poroelastic material, and the dilatation e = kk is the increase in the volume ratio of the poroelastic material. The material properties μ, λc, C, and M are constants. The strain components are related to the displacements ui by ij

=

1 ( u i, j + u j , i ) 2

(2)

Alternatively, the constitutive equations can be written in terms of ɛij and p as ij

= 2µ

=

ij

+ e

p

ij

(3a)

ij

1 p+ e M

(3b) C2

where = M , and = c . Let the stresses be separated into a M deviatoric response and a volumetric response, respectively expressed as C

ij

= Sij

P

(4a)

ij

where

Sij = 2µ

1 e 3

ij

ij

(4b)

and

P=

kk /3

= 2 µ 3

K c e + C or P = 2 µ 3

Ke + p

(4c)

+ . Here, μ represents the shear + c and K = where K c = modulus; Kc refers to the bulk modulus for a closed system that has conserved fluid volume, expressed as = 0 ; and K represents the bulk modulus at zero fluid pressure ( p = 0 ). In (1b) and (4c), M and C denote the contribution of an increase in fluid volume ξ to pore pressure p and total pressure P, respectively, under the condition of e = 0 . In addition, 2

Mechanics of Materials 138 (2019) 103156

Y.-Y. Lin and C.-H. Yang

Fig. 1. Plane-strain poroelastic medium of dimensions 2 W containing a traction-free crack of length 2a is subjected to uniform tension σo or strain ɛo in the ydirection on y = ± W . (a) Remote boundaries and the crack faces are impermeable. (b) Remote boundaries are permeable to pore pressure po, but the crack faces are impermeable. (c) Remote boundaries are impermeable, but the crack faces are permeable to pore pressure po. (d) Remote boundaries and the crack faces are permeable to pore pressure po. The impermeable boundaries and the permeable boundaries are shown in red and black, respectively. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

α is the Biot's coefficient, which indicates the ratio of an increase in fluid volume ξ to an increase in volume ratio e at zero pore pressure ( p = 0 ). If both the solid grains and the fluid are incompressible, and the poroelastic material is fully saturated with fluids, the parameters are M → ∞, C → ∞, and α → 1. Hence, = e . In porous materials, the rate of fluid outflow per unit surface area, Vi, is assumed to obey Darcy's law. Hence,

Vi =

kp, i

t

=

Vi, i = kp, ii

(6)

With negligible inertia forces and in the absence of the body force, the stresses satisfy the equilibrium equations ij, j

(7)

=0 Substituting (1a) and (2) into (7) and using e =

(5)

where k represents the coefficient of the porous material's permeability. If the fluid is incompressible, the rate of fluid volume increase must be equal to the volume of fluid flowing in through the surface of a unit reference volume per unit time, and hence

µui, jj + (µ +

c ) e, i

C

,i

= uk, k gives (8)

=0

Substituting (1b) into (6) and using (2µ + gives 3

kk

c ) e , ii

=C

, ii

from (8)

Mechanics of Materials 138 (2019) 103156

Y.-Y. Lin and C.-H. Yang

t

=c

field remains uniform and independent of time as well as stress and strain fields. Under the condition of plane strain ( 33 = 0 ), because remote stress 11 = 0 , the normal stress and strain, given in (3a), reduce to

(9)

, ii 2µ +

where c = 2µ + Mk . Eqs. (8) and (9) are the governing equations in c terms of ui and ξ. Alternatively, substituting (3a) and (2) into (7) gives

µui, jj + (µ + ) e, i

11

p + t

22

e = kp, ii t

p=

(12a)

p, i = 0

e = kp, ii = ce , ii t

22

ij

22

(13b)

ij

˜ = e˜ + 1 p˜ M

(13c)

where the notation f˜ denotes the operation

f (t ) e

st dt

0

for function f(t).

In addition, taking the Laplace transform of (7), (10), and (11) gives (14a)

˜ij, j = 0 µu˜ i, jj + (µ + ) e˜, i

p˜, i = 0

k 1 p˜ = p˜ + e˜ s , ii M

(14b) (14c)

For saturated incompressible solid and fluid constituents, the transformed governing equations reduce to

µu˜ i, jj + (µ + ) e˜, i

p˜, i = 0

k p˜ = e˜ s , ii

2

o

for t > 0

1 4µ

o

for t > 0

22

=

2µ

o

= 4µ

o

o

(16a)

p

(16b) 22 .

When

22

=

o

is applied on

(17a) (17b)

is applied on the remote boundaries, we obtain

for t > 0

(18a)

for t > 0

(18b)

Hence, for fracture problem I, the stress-control condition is interchangeable with the strain-control condition if σo is replaced with 4μɛo. In the small region near the crack tip, the pore pressure as well as stress and strain fields are non-uniform and vary with time. Eventually, the pore pressure in the near-tip region will reach the equilibrium value of 1 p = 2 o or p = 2µ o , which is the constant value of p at remote locations. In fracture problem II, the remote boundaries are permeable to 1 p = po . If po = 2 o under stress-control or po = 2µ o under straincontrol, problem II is simply reduced to fracture problem I so that the fluid drainage is limited in the near-tip region. If the value po specified 1 on the remote boundaries is different from or 2µ o , the pore 2 o pressure field in the remote region will not be uniform and time-independent. The fluids will drain in the remote region until the pore pressure reaches the equilibrium value of po. Simultaneously, the stress and strain fields change with the flow of the fluids. In fracture problem III, the permeable condition is given only on the crack faces. At the remote locations, the pressure will vary from the initial value at the undrained state to the final value of p = po at the drained state. Given constant stress σo on the remote boundaries, p= o /2 initially, so that the strain behaves as

(13a)

p˜

=

p=

(2µ + ) k . where c Taking the Laplace transform of (2), (3a), and (3b) gives

˜ij = 2µ ˜ij + e˜

1 1

22

1 2

while

(12b)

1 ˜ij = (u˜ i, j + u˜ j, i ) 2

2µ 1

1 2 p 2µ (1 )

+

Based on incompressibility, 11 = the remote boundaries, we obtain

(11)

Eqs. (10) and (11) are the governing equations in terms of ui and p. For saturated incompressible solid and fluid constituents, the governing equations reduce to

µui, jj + (µ + ) e, i

=

22

1

(10)

p, i = 0

Combining (3b) and (6) gives

1 M

=

(15a) (15b)

22

=

22

=

3. Mode-I crack problems in poroelastic medium Consider a plane-strain saturated porous medium of dimensions 2W that contains a traction-free center crack of length 2a and W ≫ a. On the boundaries remote from the crack, either constant tensile stress σo or constant tensile stain ɛo is applied in the direction perpendicular to the crack, whereas the stress components σ11 and σ12 vanish. Without limiting the drainage of pore fluid in the near-tip region, the permeability conditions of the porous medium are considered to include the following four types: (I) the remote boundaries and the crack faces are both impermeable, as shown in Fig. 1(a); (II) the crack faces remain impermeable but the remote boundaries are permeable to constant pore pressure po, as shown in Fig. 1(b); (III) the remote boundaries are impermeable, but the crack faces are permeable to constant pore pressure po, as shown in Fig. 1(c); (IV) both the remote boundaries and the crack faces are permeable to constant pore pressure po, as shown in Fig. 1(d). The permeable and impermeable boundaries are shown as black and red lines, respectively, in Fig. 1. In fracture problem I, the remote boundaries and the crack faces are both impermeable. Thus, most of the porous medium behaves like an incompressible solid, and the fluids drain only within a small region near the crack tip. In the region far from the crack, the pore pressure

1 4µ

o

1 2µ

for t

o

+

0 1

(19a)

2 p for t 2µ o

Given constant strain ɛo on the remote p = 2µ o initially, and thus the stress behaves as 22

= 4µ

22

=

o

2µ 1

for t o

boundaries, (20a)

0

1 1

(19b)

2

po for t

(20b)

If the specified value po is the same as the initial value, the pore pressure in the remote region will be independent of time, so that the fluids only flow within the small region near the crack tip, and the behavior in problem III will be similar to that in problem I; however, the behaviors still differ due to variations in the permeability condition. Detailed comparisons are given in the next section. In fracture problem IV, permeable conditions are given on both the crack faces and the remote boundaries. When the specified value po is 1 the same as 2 o under stress-control or as 2µ o under strain-control, the porous material at remote locations behaves as an incompressible solid at all times and the pore pressure remains constant. Thus, the fluid drainage is limited to a small region near the crack tip, and the result is 4

Mechanics of Materials 138 (2019) 103156

Y.-Y. Lin and C.-H. Yang

Fig. 2. Permeability conditions at zero remote stress for (a) sub-problem IIb, (b) sub-problem IIIb, and (c) sub-problem IVb, where the pore pressure imposed on permeable boundaries is given by qo = po + o/2 .

coordinates, parameter. If X Y coordinates are transformed into r the crack-tip region is located at r → 0, where either pore pressure is negligible or the flow is completed even at a short time and the material behaves in a drained manner. The behaviors in the crack-tip region depend on the permeability of the crack faces. For permeable crack faces, p˜ is negligible as r → 0, and 2p ˜ in (14c) dominates. Therefore, as r → 0, the fields at the first order are expected to be (Atkinson and Craster, 1991)

1

the same as that for problem III, namely po = 2 o or po = 2µ o . Other values of po will cause the drainage of fluids in the remote region. As fluid drainage in fracture problems II, III, and IV is affected by the specified pore pressures on permeable boundaries, we decompose the problems into two sub-problems: (a) the first is the fracture caused by remote stress σo with imposed pore pressure o /2 or by remote 2µ o; (b) the second is the strain ɛo with imposed pore pressure fracture caused by imposed pore pressure qo = po + o/2 at zero remote stress, as illustrated in Fig. 2, or by imposed pore pressure qo = po + 2µ o at zero remote strain, as illustrated in Fig. 3. The solutions of fracture problems II, III, and IV can be obtained by superimposing the solutions for sub-problems a and b. Regardless of stresscontrol or strain-control, the behaviors in sub-problems IIa, IIIa, and IVa depend on whether crack faces are permeable or impermeable. The behaviors in sub-problems IIb, IIIb, and IVb are expected to depend on the specified location of the permeability condition, the size of the porous medium, and the condition of zero stress/strain.

1/2 (c / s ) 1/4 f ( ij

˜ij

K1(2 r )

u˜ i

K1(r /2 )1/2 (c/ s )1/4gi ( )

p˜

K2(r /2 )1/2 (c /s )1/4cos

)

(21a) (21b)

2

(21c)

where K1 is the transformed stress intensity factor and K2 is the intensity factor of the transformed pore pressure gradient in the crack-tip region. For impermeable crack faces, as r → 0, the flow will soon be completed and 2 p˜ and e˜ are of the same order in (14c) in the crack-tip region. Therefore, as r → 0, the fields at the first order are expected to be (Craster and Atkinson, 1996a)

4. Near-crack-tip fields For a center crack that lies on y = 0 and |x| ≤ a, let the dimensions of the near-tip region be rescaled by the parameter c /s , i.e., 2µ + X = (x a)* s /c , Y = y * s /c , where c = 2µ + Mk and s is the Laplace

˜ij

c

5

K1(2 r )

1/2 (c /s ) 1/4f ( ij

)

(22a)

Mechanics of Materials 138 (2019) 103156

Y.-Y. Lin and C.-H. Yang

Fig. 3. Permeability conditions at zero remote strain for (a) sub-problems IIb, (b) sub-problem IIIb, and (c) sub-problem IVb, where the pore pressure imposed on permeable boundaries is given by qo = po + 2µ o .

u˜ i

(22b)

K1(r /2 )1/2 (c/ s )1/4gi ( )

F1 =

P1j, j dV = V

p˜

1 K3r (c /s ) 2 cos

3 K1r 2

3 (c / s ) 4

3 cos + cos 2 2

1˜ k 1 2 ˜ ˜ + p˜, i p˜, i + tij ˜ij + pe p˜ 2 2s 2M

(26) This means that the values of the integral I along any two counterclockwise contours S1 and S3 should be the same. If the crack is subjected to remote loadings, at a small flow time s → ∞, the near-tip region (also known as the inner region) is enclosed by the boundary exhibiting singular behavior in an undrained manner at r ≫ 1. The singular fields behave as (Atkinson and Craster, 1991)

(23)

where t˜ij = ˜ij + p˜ ij , so that (14a) and (14c) can be generated by the Euler-Lagrange equations. An energy momentum tensor is defined by

Plj =

L u˜ i, l + u˜ i, j

L p˜ p˜, j , l

L

lj

=

˜ij u˜ i, l +

k p˜ p˜ s ,j ,l

L

lj

(25)

For a volume V containing a crack, S must be a closed surface that does not enclose any singularities. Let S = S1 + S2 + S3 + S4 , where S1 and S3 are the circular contours around the crack, and S2 and S4 are two paths along the upper and lower crack faces, respectively. If the integrands on the paths S2 and S4 are both zero, then (25) gives the pathindependent integral

(22c)

where K1 is the transformed stress intensity factor in the crack-tip region, but the transformed pore pressure gradient field has no singularity and K3 is the coefficient of the r term. fij(θ) and gi(θ) are the usual elastic angular dependence functions. According to Craster and Atkinson (1992), the Lagrangian in the transformed domain can be written as

L=

P1j nj dS = 0 S

(24)

so that Plj, j= 0 when L does not depend on xl explicitly. An invariant integral of Plj is defined by 6

1/2 (c / s ) 1/4 f ( ij

˜ij

K e (2 r )

u˜ i

K e (r /2 )1/2 (c/ s )1/4gi ( )

)

(27a) (27b)

Mechanics of Materials 138 (2019) 103156

Y.-Y. Lin and C.-H. Yang

p˜

1/2 (c / s ) 1/4cos

K e (2 r )

(27c)

2

where K e is the transformed intensity factor for both stress and pore a for stress-control cases and pressure fields; K˜ e (s ) = o s 4 µ a o K˜ e (s ) = for strain-control cases. For the specific permeability

co =

s

conditions stated in the first sub-problems, the fluids flow mostly in the near-tip region to eliminate the singularity in the pore pressure field. If the contour S1 is set on the circular boundary at r ≫ 1 and the contour S3 is set on the circular boundary of the crack-tip region at r → 0, the integral I can be applied to find the singularities of the stresses and pore pressure gradients at s → ∞. Craster and Atkinson (1992, 1996a) showed that for impermeable cracks, as s → ∞,

1 K˜ e 2

vu 2µ

1 = K˜12

1 i

21

vu 2µ

21

= K1

(28)

v 2µ

+ K2

2

k 8s

1

+

1/2

1

cot 1(

1 (

1/2 +)

ln

( (1

1/2 +) 1/2 ) +

±

= ¯

= ¯

¯=

¯) + ( ¯2 + 4 ¯)1/2 2(2 ¯ 1)

(2

1

1/2

(31b)

2) ± ( ¯2 + 4 ¯)1/2 2(1 2 ¯)

(¯

(31c)

1 v 2(vu v )

(31d)

For a permeable crack of length 2a subjected to 22 = o and p = 0 on crack faces, the singularities of the stress and pore pressure gradient in the crack-tip region are given by

K˜1 (s ) =

ico 2 ¯ (c / s )1/4 a N+ (0)

(32a)

) i2 o a (s / c )1/2 iN+ (0) + 2 ) s (N+ (0)/ ¯ d ) 2¯

(32b)

a s ¯ (N+ (0)/ ¯ o

2( u K˜2 (s ) = (1

d)

1

where the first term in K˜1 (s ) is the limiting value as s → ∞, and together with K˜2 (s ) satisfy the invariant integral in (29). It is noted that 1

1

d=

0

=

i

1 0

1 y1/2 (1

1 y )1/2 N ( iy )

y1/2 1 (1 y )1/2 N ( iy )

N+ (0) dy ¯

N+ (0) dy ¯

1 y )1/2 N ( iy )

p (1 p2 )1/2 1 dp p2 ¯ p+i

N+ (0) dy ¯

(33c)

(33d)

2 qo c1/2ico (s/ c )1/4 2 ) ) s 3/2 (d N+ (0)/ ¯) N+ (0)

N+ (0)/ ¯) N+ (0)

d

N+ (0) ¯

(34a) 2

+ co

iN+ (0) + 2¯

= K1

v 2µ

+ K2

2

k 8s

(35)

Finite element models of the first and second sub-problems of a finite crack of length 2a in a saturated poroelastic plane of dimensions 2W were built using 1/4 of the domain and symmetry conditions in the software ABAQUS. Four-node pore pressure plane-strain elements (CPE4P) are the coupled displacement-pressure elements under the plane strain assumption used in the models to simulate the constitutive model of poroelastic materials, in which the elastic properties of the solid network are described by drained Young's modulus E and Poisson's ratio ν, respectively. Let elastic properties = 0.229 and µ = 0.488 MPa, and permeability k = 8.59 × 10 4 m4/(N · s) in the models. The half crack is located on y = 0 and 0 ≤ x ≤ a. Let the domain size W be much larger than the half-length, i.e., W/a =50, 100, and 200. For first sub-problems I, IIa, IIIa, and IVa, uniform stress 22 = o or uniform displacement u2 = o W is applied on the edges of y = W, depending on whether it is stress-controlled or strain-controlled. For the first subproblems, the impermeable and permeable conditions with p = 0.5 o or p = 2µ o on the boundaries are also specified. For second sub-problems IIb, IIIb, and IVb, 22 = 0 or u2 = 0 is applied on the edges of y = W, depending on whether it is stress-controlled or strain-controlled. The permeable and impermeable conditions on the boundaries are also specified for the second sub-problems. We first ran the simulations without embedding cohesive zone elements between crack faces. These FEM simulations provide a full solution that can be used to extract the behavior of the stress and pore pressure fields near the crack tip, especially the singularities in the crack-tip region. For this task, the mesh needs to be very fine and well arranged in the neighborhood of the crack tip. We embedded cohesive zone elements between the crack faces around the tip in the FEM models, as illustrated in Fig. 4, to compute the instantaneous energy release rate right in front of the crack tip. The instantaneous energy release rate is the direct driving force of solid network cleavage and is computed as

(31a) 1

1

5. Finite element models

(30)

1)1/2 +

1

21

1

cs 2 (s /c ) 2

vu

1

tan 0

y 3/2 (1

s (d

qo 2k

where the first term is the same as that obtained by (28), and N+ (0) = 1 v ,

1 2

1

(34b) ˜ ˜ The relationship between K1 and K2 can be verified by the following invariant integral:

o

n¯1 =

1

2 2 qo ¯ (s/ c )1/4

K˜2 (s ) =

where vu is the undrained Poisson's ratio at r ≫ 1 and equivalent to 1/2 for saturated poroelastic materials with incompressible constituents. With the uniform stress and pore pressure field subtracted, the fracture behaviors of the first sub-problems are equivalent to those of a finite crack subjected to 22 = o on the crack faces at small times, as found by Craster and Atkinson (1996b). For an impermeable crack of length 2a, as time increases (s decreases), the crack-tip singularity K˜1 (s ) increases as

1

exp

1

(1 2( u

K˜1 (s ) =

(29)

n¯ (c/ s )1/2 a K˜1 (s ) = 1+ 1 sN+ (0) 2a

2

0

and for permeable cracks (Craster and Atkinson, 1994), as s → ∞,

Ke

2¯

The fracture behaviors of the second sub-problems are more complicated due to the influence of the long-distance flow of fluids. However, at very short time s → ∞, the crack is considered to be semiinfinite and only the permeability condition specified on the crack faces matters. The singularities in the crack-tip region of a permeable crack can be predicted for a semi-finite crack subjected to 22 = 0 and p = qo on the crack faces. As shown by Craster and Atkinson (1994),

v 2µ

1

N ( )=

(33a)

Lc

G=J=

m 0

(33b)

x

dx

(36)

where σm is the traction in the cohesive zone, δ denotes the opening 7

Mechanics of Materials 138 (2019) 103156

Y.-Y. Lin and C.-H. Yang

Fig. 4. Illustration of the cohesive zone model.

displacement between two surfaces of the cohesive zone, and Lc is the length of the cohesive zone. The traction-separation law in cohesive zone elements is set to m

(37)

= kc 8

where kc is assigned to be a constant of 7.27 × 10 Pa/m. The tractionseparation law is set to be independent of pore pressure because the purpose of using the cohesive zone is to compute the local instantaneous energy release rate using (36), which is not affected by the pore pressure. The cohesive zone model can be described by linear or nonlinear springs; the simulation results remain the same as long as we use the same length Lc. In our simulations, Lc is small in comparison with the size of the crack-tip region. When the size of the cohesive zone is directly controlled, the details of the traction-separation law, including the parameter kc, do not matter. More detailed discussions on the length of the cohesive zone can be found in the study of Yang and Lin (2018).

Fig. 5. (a) Normalized stress and (b) normalized pore pressure near the crack tip for an impermeable crack subjected to the remote stress σo or the remote strain ɛo (problems I and sub-problem IIa) at a small flow time. Note t * = ct / a2 and K e = o a . The remote stress σo can be replaced by 4μɛo for the straincontrol case.

6. Numerical results

p

6.1. Short-time near-tip fields for impermeable cracks

Note that the term σo can be replaced by 4μɛo when the crack is under remote strain control.

For impermeable cracks in problems I and II, the near-tip behaviors at small flow times are affected by two conditions applied at the remote boundaries, namely the loading imposed by stress σo or strain ɛo, and the permeability condition.

6.1.2. Sub-problem IIb The effects of the imposed remote pore pressure qo on the near-tip fields of impermeable cracks appear in the FEM results for problem IIb. The stress σ22 and the pore pressure p at a small flow time are normalized and plotted against r* in Fig. 6. The data in Fig. 6 indicate an empirical formula of stress singularity for fracture problem IIb at t much smaller than W2/c for r* ≪ 1, i.e.,

6.1.1. Problem I and sub-problem IIa The near-tip fields obtained from the FEM simulations of problems I and IIa show the effect of remote loading. The stress σ22 and the pore pressure p at small normalized time t * = ct / a2 are normalized and plotted against the normalized distance from the crack tip, r * = (x a)/ ct , in Fig. 5. At a small flow time, the near-tip fields exhibit two distinct behaviors for r* ≫ 1 and r* ≪ 1. The numerical data agree well with the following predictions by Craster and Atkinson (1996b). For r* ≫ 1, 22

p

K e (ct )

1 4/

2 r * , where K e =

1 4/

2 r*

K e (ct )

o

22

p

22

1 4/

1 4/

2 r* , a 1 qo (ct ) 2 W

(40a) (40b)

0

The stress intensity factor K1 increases with time by which has a much slower rate than the time dependence of K1 caused by remote loading shown in (39a). ct ( W2 )1/2 ,

(38b)

For r* ≪ 1

K1 (t )(ct )

K1 (t )(ct )

where K1 (t )

(38a)

a

(39b)

const

6.2. Short-time near-tip fields for permeable cracks

2 r* ,

where K1 (t ) = K e

1 H (t ) + N+ (0)

n¯1 ct N+ (0) a2

For permeable cracks in problems III and IV, the near-tip fields at small flow times depend on the imposed pore pressure on crack faces together with remote conditions.

1/2

(39a) 8

Mechanics of Materials 138 (2019) 103156

Y.-Y. Lin and C.-H. Yang

Fig. 6. (a) Normalized stress and (b) normalized pore pressure near the crack tip for an impermeable crack subjected to the remote pore pressure qo (subproblem IIb) at a small flow time. Note t * = ct / a2 .

Fig. 7. (a) Normalized stress and (b) normalized pore pressure near the crack tip for a permeable crack subjected to the remote stress σo or the remote strain ɛo (sub-problems IIIa and IVa) at a small flow time. Note t * = ct / a2 , and K e = o a . The remote stress σo can be replaced by 4μɛo for the strain-control case.

6.2.1. Sub-problems IIIa and IVa The near-tip fields obtained from the FEM simulations of problems IIIa and IVa show the effect of remote loading when the crack faces are 2µ o under strain permeable to p = o /2 under stress control or p = control. The normalized stress σ22 and pore pressure p against r* at short normalized time t * = ct / a2 are plotted in Fig. 7. As t* → 0, the numerical results follow the predictions by Craster and Atkinson (1996a) for permeable cracks at a small flow time, i.e., For r* ≫ 1, 22

p

K e (ct )

1 4/

K e (ct )

(41a)

2 r*

1 4/

6.2.2. Sub-problems IIIb and IVb The near-tip fields obtained from the FEM simulations of problems IIIb and IVb show the effect of the imposed pore pressure qo on crack faces when the remote stress or the remote strain is zero. The normalized stress σ22 and pore pressure p against r*at short normalized time t * = ct / a2 are plotted in Fig. 8. The numerical results are the same as the predictions for a semi-infinite crack by Craster and Atkinson (1994), as explained earlier, i.e., for r* ≪ 1, 22

(41b)

2 r*

p

For r* ≪ 1, 22

p

1 4/

K1 (t )(ct )

1

K1 (t ) =

(42b)

K2 (t )(ct ) 4 r */2

qo

where 1 K1 (t ) = K e ¯ (N+ (0)/ ¯

2( u K2 (t ) = K e (1

d)

) 2 )

H (t )

(5/4)

2 (N+ (0)/ ¯

ico 2 N+ (0)(N+ (0)/ ¯

d)

K2 (t ) =

ct 1/4 d) a2

N+ (0) +i 2¯

(ct )

(42c) 1 2

1 4/

1 K2 (t )(ct ) 4

(43a)

2 r*

(43b)

r */2

where

(42a)

2 r*

K1 (t )(ct )

(5/4)(d

2 2¯ (3/4)(d N+ (0)/ ¯) N+ (0) qo (ct )

(42d)

2 ico q (ct )1/4 N+ (0)/ ¯) N+ (0) o

d

N+ (0) ¯

(43c) 2

+ co

iN+ (0) + 2¯

1/4

(43d)

6.3. Instantaneous energy release rates J(t) for problems I and II

Again, the near-tip fields exhibit two distinct behaviors, given respectively by Ke for r* ≫ 1, and K1(0) and K2(t) for r* ≪ 1. Note that K e = o a under remote stress control and K e = 4µ o a under remote strain control.

The instantaneous energy release rates J(t) for the same fracture problems were separately computed using FEM simulations with embedded cohesive elements. For impermeable cracks in problem I and 9

Mechanics of Materials 138 (2019) 103156

Y.-Y. Lin and C.-H. Yang

Fig. 9. Normalized instantaneous energy release rates for problem I and sub1 problem IIa. Note J = 2µ K e2 , where K e = o a under remote stress control and K e = 4µ o given by (45).

a under remote strain control. The prediction of initial J(t) is

Fig. 8. (a) Normalized stress and (b) normalized pore pressure near the crack tip for a permeable crack subjected to the pore pressure qo on the crack faces (sub-problems IIIb and IVb) at a small flow time. Note t * = ct / a2 .

sub-problems IIa and IIb, the numerical results were carried out for dimensional ratios W/a = 50, 100, and 200. Problem I and sub-problem IIa are identical. Their numerical results are plotted in Fig. 9 using the normalization by J∞, which is the value of J at the drained state when t* ≫ 1. At the drained state, J should reach

J =

1 2µ

K e2

(44)

Note that K e = o a under remote stress control, whereas K e = 4µ o a under remote strain control. Based on (39a), the instantaneous energy release rate for problems I and IIa at small flow times (t* ≪ 1) is estimated by

J (t ) =

1 2µ

Ke2

1 H (t ) + N+ (0)

n¯1 ct N+ (0) a2

1/2 2

(45) 1

1/2 K e2 2µ

The limiting value at the undrained state occurs at J (0) = for saturated poroelastic media. The numerical data in Fig. 9 agree well with these predictions. The data for the three ratios of W/a merge into one curve, indicating that the fluid flow range for impermeable cracks in problems I and IIa is not affected by the dimension of the remote boundaries. The numerical results of J for sub-problem IIb are normalized by J * = (1 ) qo2 a/2µ and plotted against the normalized time t ** = ct /W 2 in Figs. 10(a) and (b) for zero remote stress and zero remote strain, respectively. The merging of the results for various dimensional ratios indicates that the change in J is caused by fluid flow in the range of W. At the undrained and drained limits, J should be zero under zero stress. However, when the remote boundaries are subjected to zero ɛ22,

Fig. 10. Instantaneous energy release rates for sub-problem IIb normalized by J * = (1 ) qo2 a/2µ at (a) zero remote stress and (b) zero remote strain. Note qo is the imposed remote pore pressure. The predictions of drained limit and initial J(t) are given by (46) and (47), respectively.

pore pressure reaches qo everywhere at the drained state, so that the 1 2 stress on the remote boundaries should be 22 = 1 qo and the drained limit reaches

10

Mechanics of Materials 138 (2019) 103156

Y.-Y. Lin and C.-H. Yang

less than 0.5 o , the maximum of J(t) occurs before the drained state. For problem II under strain control, the undrained limit remains unchanged, but the drained limit changes with the difference qo between po and 2µ o . When qo ≠ 0, the time required to reach the drained state is also postponed because the range of fluid flow extends from the neartip region to the whole domain. The drained limit at t → ∞ is calculated using

J (t

posed remote pore pressure.

)=

1

1 1

2µ

2

qo

a

1 2µ

qo2a

ct W2

1 1

2µ

2

(po + 2µ o) + 4µ

2 o

a

(48)

2

(46)

6.4. Instantaneous energy release rates J(t) for problems III and IV

At a small flow time, similar to the empirical formula (40a) for K1(t), the instantaneous energy release rate J(t) behaves as

J (t )

1

When the difference qo is negative, the corresponding tensile stress on the remote boundaries increases, and hence J(t → ∞) becomes )8µ o2 a. The J(t) curve is also a higher than the value of J = (1 two-stage increasing curve whose maximum is at the drained state. If the difference qo is positive, the corresponding tensile stress on the remote boundaries decreases, and hence J(t → ∞) becomes lower than the value of J∞, and could be even lower than the undrained limit when the difference is large. The maximum of J(t) occurs at around t* = 1 10, at which point the pore pressure field in the near-tip region is equilibrated in the first stage.

Fig. 11. Normalized instantaneous energy release rates for problem II subjected 1 to the remote stress σo when W/a = 50. Note J = 2µ o2 a , and po is the im-

J (t

)=

Since sub-problems IIIa and IVa are identical, the numerical results of J(t) for permeable cracks in sub-problems IIIa and IVa for dimensional ratios W/a = 50, 100, and 200 are plotted together in Fig. 13 using normalization by J∞. Based on (42c), the instantaneous energy release rate for sub-problems IIIa and IVa at small flow times (t* ≪ 1) is estimated by

(47)

The numerical results agree well with (47) up to t**∼0.1, as shown in Fig. 10. To obtain the values of J for problems II, we superimposed the stress intensity factors K1(t) in the crack-tip region for sub-problems IIa and IIb, which were obtained by back-calculating the numerical data of J(t) 1 with J (t ) = 2µ K1 (t ) 2 . K1(t) for problem II was then substituted into the formula to compute J(t). The normalized J for problem II with W/ a = 50 is plotted in Figs. 11 and 12 for the stress- and strain-control cases, respectively. For problem II under stress control, although the undrained and drained limits remain unchanged, the variations of J(t) curves are influenced by the imposed pore pressure po on the remote boundaries. When the imposed pore pressure po is different from 0.5 o , the time required to reach the drained state is postponed because the range of fluid flow extends from the near-tip region to the whole domain. If po is

J (t ) =

1 2µ

Ke2

1 ¯ (N+ (0)/ ¯

d)

H (t )

(5/4)

ico 2 N+ (0)(N+ (0)/ ¯

2 ct 1/4 d ) a2

The limiting value at the undrained state is J (0) = 1 ¯ (N+ (0) / ¯

1 2µ

(49)

m ( , 1/2) 2K e2 ,

. Again, σo is interchangeable with 4μɛo where m ( , 1/2) = d) for the calculation of Ke under stress and strain control. The numerical results of J(t) for sub-problem IIIb under zero remote ) qo2 a/2µ and plotted in Figs. 14(a) stress are normalized by J * = (1 and (b) against normalized times t * = ct / a2 and t ** = ct /W 2 , respectively. The limits of J should be zero at the undrained and drained states with no applied stress. At a small flow time, based on (43c), we expect that the instantaneous energy release rate (t) can be estimated by

Fig. 13. Normalized instantaneous energy release rate for permeable cracks for 1 sub-problems IIIa and IVa. Note J = 2µ K e2 , where K e = o a under remote

Fig. 12. Normalized instantaneous energy release rates for problem II subjected )8µ o2 a , and po is the to the remote strain ɛo when W/a = 50. Note J = (1 imposed remote pore pressure.

stress control and K e = 4µ o initial J(t) is given by (49).

11

a under remote strain control. The prediction of

Mechanics of Materials 138 (2019) 103156

Y.-Y. Lin and C.-H. Yang

Fig. 14. Normalized instantaneous energy release rates for sub-problem IIIb at zero remote stress against normalized times (a) t * = ct / a2 and (b) t ** = ct /W 2 . ) qo2 a/2µ , where qo is the imposed pore pressure on crack Note J * = (1 faces. The prediction of initial J(t) is given by (50).

J (t ) =

1 2µ

qo2a

(5/4)(d

2 ico N+ (0)/ ¯) N+ (0)

2

ct a2

Fig. 15. Normalized instantaneous energy release rates for sub-problem IIIb at zero remote strain against normalized times (a) t * = ct / a2 and (b) t ** = ct /W 2 . ) qo2 a/2µ , where qo is the imposed pore pressure on crack Note J * = (1 faces. The predictions of drained limit and initial J(t) are given by (46) and (50), respectively.

1/2

(50)

The numerical data in Fig. 14(a) agree well with (50) for t* < 0.1, and all data for the various ratios W/a merge during this period. After t* > 1, the fluids flow from the remote region to the near-tip region. Hence, we switched the normalized time to t** in Fig. 14(b) to show the merging of all data for the various ratios W/a at a long time (t** > 1). The maximum J/J* occurs at t** ∼ 0.05 for sub-problem IIIb. Similarly, the results of J/J* against normalized times t * = ct / a2 and t ** = ct /W 2 for sub-problem IIIb under zero strain are plotted respectively in Figs. 15(a) and (b). The undrained limit J (0) = 0 but the drained limit J∞ follows (46). At a small time, J still follows the prediction by (50) given for a permeable crack. Similar to the results for stress control, all data for the various ratios W/a merge to a function of normalized time t* for t* < 0.1, and (50) agrees well with the data. After t* > 1, the results are influenced by the remote condition and flow range W, and thus merge at normalized time t** for t** > 1 and eventually reach the drained limit J∞. Following the same calculation procedure of J(t) for problem II, the results of J(t) for problem III with W/a = 50 are obtained and plotted using normalization by J∞ in Figs. 16 and 17 for stress- and straincontrol cases, respectively. For the stress-control case, the variations of J(t) are influenced by the difference qo between the imposed pore pressure po on the crack faces and 0.5 o , and the time required to reach the drained state is postponed because the range of fluid flow extends from near-tip region to the whole domain. When po < 0.5 o ,

Fig. 16. Normalized instantaneous energy release rates for problem III sub1 jected to the remote stress σo when W/a = 50. Note J = 2µ o2 a , and po is the imposed pore pressure on crack faces.

the maximum of J(t) occurs before the drained state. For the straincontrol case, the drained limit changed with the difference qo between po on the crack faces and 2µ o . The drained limit at t → ∞ is still given by (48). When qo ≠ 0, the time required to reach the drained state is 12

Mechanics of Materials 138 (2019) 103156

Y.-Y. Lin and C.-H. Yang

Fig. 17. Normalized instantaneous energy release rate for problem III subjected )8µ o2 a , and po is the to the remote strain ɛo when W/a = 50. Note J = (1 imposed pore pressure on crack faces.

postponed because the range of fluid flow extends from the near-tip region to the whole domain. The numerical results of J/J* for sub-problem IVb at zero remote stress are plotted in Fig. 18 against normalized time t* and t**. The maximum of J/J* for sub-problem IVb occurs slightly earlier, at t** ∼ 0.01, compared to that for sub-problem IIIb. For both cases, the maximum value increases with dimensional ratio W/a; however, the value in sub-problem IVb is lower than that in sub-problem IIIb at a given ratio of W/a. The permeable condition on the remote boundaries causes a lower maximum value of J/J* for sub-problem IVb at an earlier time than that for sub-problem IIIb. Similarly, the numerical results of J/J* for sub-problem IVb at zero remote strain are plotted in Fig. 19 against normalized times t* and t**. The calculated results of J(t) for problem IV with W/a = 50 are normalized by J∞ and plotted in Figs. 20 and 21 for stress- and straincontrol cases, respectively. For the stress-control cases, the maximum of J(t) is dominated by the diffusion process driven by the permeable crack faces. When po = 0, a two-stage increasing J(t) curve appears and reaches its maximum at the drained state. For the strain-control cases, the undrained limit remains unchanged, and the drained limit at t → ∞ can be calculated by (48). When the difference qo between the imposed pore pressure po on the permeable boundaries and 2µ o is negative, the J(t) curve is also a two-stage increasing curve whose maximum is at the drained state. If the difference is positive, the corresponding tensile stress on the remote boundaries decreases, and hence J(t → ∞) becomes lower than the value of J∞, and could be even lower than the undrained limit when the difference is large.

Fig. 18. Normalized instantaneous energy release rates for sub-problem IVb at zero remote stress against normalized times (a) t * = ct / a2 and (b) t ** = ct /W 2 . ) qo2 a/2µ , where qo is the imposed pore pressure at all Note J * = (1 boundaries. The prediction of initial J(t) is given by (50).

the singular pore pressure field initiated by remote stress or remote strain. The stress intensity factor in the crack-tip region and the instantaneous energy release rate J(t) monotonically increase until the field equilibrates at ∼ a2/c. (b) In sub-problems IIb, IIIb, and IVb, the fluid flow is driven by the difference between the pore pressure imposed on the permeable boundaries and the initial zero field. Under zero remote stress, the escape of fluids from the permeable crack faces increases the crack opening displacement, so that J increases with time until t ∼ a2/c; however, when remote boundaries are impermeable, the fluids outside the neighborhood of the crack could flow into the near-tip region after t ∼ a2/c to reduce the deformation of the crack faces, so that J decreases with time from t ∼ a2/c until equilibrium is reached at t ∼ W2/c. When remote boundaries are permeable, the escape of fluids from the remote boundaries reduces the deformation of the poroelastic domain, gradually affecting the crack opening displacement, so that J starts decreasing earlier than t ∼ W2/c. When the remote boundaries maintain zero strain, a tensile stress compensates for the reduced deformation caused by fluids flowing out, and hence J increases with time whether the remote boundaries are permeable or not. In contrast, with the flow-in of fluids, zero strain on the remote boundaries is accompanied by a compressive stress, and thus J decreases with time. After superimposing the solutions of two sub-problems, the effects of permeability conditions coupled with stress and strain on J(t) were observed. In the strain-controlled cases, the drained limit occurring at a long time varies with the difference in pore pressure between the imposed value on the boundaries and the initial value caused by the applied strain. A positive difference lowers the drained limit of J, and

7. Conclusions This research used finite element models with cohesive zone elements to directly compute instantaneous energy release rates of center cracks in poroelastic media subjected to mode-I loading and various permeability conditions. The length of the cohesive zone needs to be smaller than the size of the crack-tip region, where the drained singular stress field dominates, particularly at the undrained state. In addition, the stress and pore pressure fields in the neighborhood of cracks were computed and verified with asymptotic analysis at short times. The numerical results of J(t) at short times agree well with those predicted from the analytical stress intensity factors in the crack-tip region. Despite the complexity of the boundary conditions, the time-dependent fracture behavior of a crack in a poroelastic medium is basically attributed to the fluid flow, and thus we divided the problems into two categories. (a) In problem I and sub-problems IIa, IIIa, and IVa, the fluid flow is limited to the neighborhood of the crack and driven to eliminate 13

Mechanics of Materials 138 (2019) 103156

Y.-Y. Lin and C.-H. Yang

Fig. 21. Normalized instantaneous energy release rate for problem IV subjected )8µ o2 a , and po is the to the remote strain ɛo when W/a = 50. Note J = (1 imposed pore pressure on all boundaries.

permeable boundaries and the initial value at remote locations caused by the applied stress. However, the condition changes with the location of the permeable boundaries. The maximum values of J depends on the dimensional ratio W/a. It is worth noting that although the time-dependent J curve may exhibit two-stage increasing behavior similar to the fracture behavior of cracks in poroviscoelastic media, unlike the competition between the flow time and intrinsic viscoelastic time in the poroviscoelastic crack problems, the two-stage behavior in J(t) of a poroelastic crack problem is attributed to short and long ranges of fluid flow. In conclusion, the time required to reach the first stage is related to the crack length, and the time required to reach the second stage depends on the domain size. Acknowledgments

Fig. 19. Normalized instantaneous energy release rates for sub-problem IVb at zero strain against normalized times (a) t * = ct / a2 and (b) t ** = ct /W 2 . Note J * = (1 ) qo2 a/2µ , where qo is the imposed pore pressure on all boundaries. The predictions of drained limit and initial J(t) are given by (46) and (50), respectively.

The authors would like to thank the Ministry of Science and Technology of Taiwan, for financially supporting this research under Contract no. MOST 107-2221-E-006-027. References Atkinson, C, Craster, R.V., 1991. Proc. R. Soc. Lond. A 434, 605–633. Atkinson, C., Craster, R.V., 1992. Philos. Trans. R. Soc. Lond. Ser. A 339, 231–263. Beebe, D.J., Moore, J.S., Bauer, J.M., Yu, Q., Liu, R.H., Devadoss, C., Jo, B.H., 2000. Nature 404, 588. Biot, M.A., 1941. J. Appl. Phys. 12, 155–164. Biot, M.A., 1962. J. Appl. Phys. 33, 1482–1498. Bonn, D., Kellay, H., Prochnow, M., Ben-Djemiaa, K., Meunier, J., 1998. Science 280, 265–267. Boone, T.J., Ingraffea, A.R., 1990. Int. J. Numer. Anal. Method Geomech. 14, 27–47. Bouklas, N., Landis, C.M., Huang, R., 2015. J. Appl. Mech. 82, 081007. Cheng, A.H.D., Sidauruk, P., Abousleiman, Y., 1994. Math. J. 4, 76–82. Craster, R.V., Atkinson, C., 1992. J. Mech. Phys. Solids 40, 887–924. Craster, R.V., Atkinson, C., 1994. Philos. Trans. R. Soc. Lond. A 346, 387–428. Craster, R.V., Atkinson, C., 1996a. Mechanics of Poroelastic Media 23 Kluwer Academic Publishers. Craster, R.V., Atkinson, C., 1996b. Q. J. Mech. Appl. Math. 49, 311–335. Duncan, R., 2003. Nat. Rev. Drug Discov. 2, 347. Luo, Y., Shoichet, M.S., 2004. Nat. Mater. 3, 249. Noselli, G., Lucantonio, A., McMeeking, R.M., DeSimone, A., 2016. J. Mech. Phys. Solids 94, 33–46. Nowak, A.P., Breedveld, V., Pakstis, L., Ozbas, B., Pine, D.J., Pochan, D., Deming, T.J., 2002. Nature 417, 424. Pizzocolo, F., Huyghe, J., Ito, K., 2013. Eng. Fract. Mech. 97, 72–79. Rice, J.R., Cleary, M.P., 1976. Rev. Geophys. Space Phys. 14, 227. Skrzeszeswka, P.J., Sprakel, J., de Wolf, F.A., Fokkink, R., Stuart, M.A.C., van der Gucht, J, 2010. Macromolecules 43, 3542. Wang, X., Hong, W., 2012. Soft Matter 8, 8171. Wilson, W.K., Yu, I.-W., 1979. Int. J. Fract. 15, 377. Yang, C.H., Lin, Y.Y., 2018. Eur. J. Mech. A 69, 78–87.

Fig. 20. Normalized instantaneous energy release rate for problem IV subjected 1 to the remote stress σo when W/a = 50. Note J = 2µ o2 a , and po is the imposed pore pressure on all boundaries.

hence the maximum of J occurs at approximately t ∼ a2/c. In the stresscontrolled cases, the drained limit occurring at a long time does not change. Whether the drained limit is the maximum or not depends on the difference in pore pressure between the imposed value on the 14