# Coupled Finite Volume Methods and Extended Finite Element Methods for the Dynamic Crack Propagation Modelling with the Pressurized Crack Surfaces

Authors: Shouyan Jiang; Chengbin Du
Journal: Shock and Vibration
 (2017)

---

## Abstract

We model the fluid flow within the crack as one-dimensional flow and assume that the flow is laminar; the fluid is incompressible and accounts for the time-dependent rate of crack opening. Here, we discretise the flow equation by finite volume methods. The extended finite element methods are used for solving solid medium with crack under dynamic loads. Having constructed the approximation of dynamic extended finite element methods, the derivation of governing equation for dynamic extended finite element methods is presented. The implicit time algorithm is elaborated for the time descritisation of dominant equation. In addition, the interaction integral method is given for evaluating stress intensity factors. Then, the coupling model for modelling hydraulic fracture can be established by the extended finite element methods and the finite volume methods. We compare our present numerical results with our experimental results for verifying the proposed model. Finally, we investigate the water pressure distribution along crack surface and the effect of water pressure distribution on the fracture property.


---

## Body

## 1. Introduction

For hydraulic concrete structures, the external dynamic loads, such as strong earthquake, may cause cracking of these structures. The cracking of the structure causes the fluid injecting into the solid medium. The injecting fluid produces fluid pressures along crack surface and affects the deformation of the solid medium and the fracture propagation again. Many researchers [1–4] contributed to the study of hydraulic fracturing problem and these efforts led to a progressive recognition of the multiscale nature of the hydraulic fracturing problem. In simulation of hydraulic fracturing, these important aspects need to be specially concerned, namely, the flow of viscous fracturing fluid, the creation of fracture surfaces in the solid, the formation of a lag between the crack edge and the fluid front, the elastic deformation of the solid, and the leak-off of fluid from the fracture.Some researches (2005) [5] showed that, in hydraulic structures, the water could penetrate an initiated crack, and the crack opening velocity, the magnitude of the opening, and crack mouth pressure had an important effect on the water pressure distribution. Boone and Ingraffea (1990) [6] proposed a finite difference approximation for modelling fluid flow along the fracture. Wu and Wong (2014) [7] incorporated the cubic law into the numerical manifold method for modelling fluid flow through fractures. Lisjak et al. (2017) [8] assumed that the fluid flow in discontinuous, porous rock masses was a viscous, compressible fluid, and the flow was explicitly solved based on a cubic law approximation. The fluid flow along a propagating crack surface satisfies some natural flow law. Some hypotheses [9–12], such as linear distribution of the water pressure along crack case, full reservoir pressure case, for evaluating water pressure along a propagating crack cannot reflect accurately the variation of the water pressure along new developing cracks in structures.In recent years, many numerical methods have been developed for hydraulic fracturing modeling, such as finite element methods [13, 14], generalized finite element methods [15], finite-discrete element methods [8], numerical manifold methods [7], boundary element methods [16], discontinuous deformation analysis methods [17], and extended finite element methods (XFEM). The XFEM shows huge advantage for dealing with discontinuous problems [18–20] and also hydraulic fracture problems. The XFEM mesh does not need to align with a discontinuity. For moving discontinuities, such as crack propagation problem, it does not need to carry on remeshing. Mesh refinement is also unnecessary around a discontinuous feature. The first simulation of hydraulic fracture in XFEM was due to Réthoré et al. [21] and they developed a two-scale numerical model for fluid flow in fractured, deforming porous media. In 2009, Lecampion [22] adopted the XFEM for investigating the solution of hydraulic fracture problems. Gordeliy and Peirce (2013) [23] proposed coupled algorithms that used the XFEM to solve the elastic crack component of the elastohydrodynamic equations that governed the propagation of hydraulic fractures in an elastic medium. Subsequently, they (2013) [24] proposed two novel XFEM schemes for modeling fluid driven fractures both of which exploited an implicit level set algorithm for locating the singular free boundary that occurred when the fluid and fracture fronts coalesce. Their excellent works provide mathematical proof for using XFEM to solve hydraulic fracturing problems. Khoei et al. [25] simulated the crack growth in saturated porous media using XFEM. Taleghani [26] also developed an XFEM code to simulate fracture propagation, initiation, and intersection, and the presented coupled fluid flow-fracture mechanics simulations extended available modeling efforts and provided a unified framework for evaluating fracture design parameters and their consequences. Salimzadeh and Khalili (2015) [27] proposed a three-phase hydromechanical model for hydraulic fracturing and they handled discontinuity by using XFEM while cohesive crack model was used as fracturing criterion. Wang et al. (2015) [28] proposed a hybrid approach combining the XFEM and the finite volume method to simulate hydraulic fracturing in concrete dams. Our current study concerns developing a model for cracking modeling of structure under water pressure along a propagating crack surface and dynamic loads. Additionally, we will compare our present numerical results with our experimental results for verifying the proposed model.This paper is organized as follows. Section2 introduces some governing equations for elastic dynamic responses of the solid medium and fluid flow pressure within the crack. Section 3 discusses numerical approximation of extended finite element methods and finite volume methods of the flow along a crack. Section 4 gives a numerical example for investigating the water pressure distribution along crack surface and the effect of water pressure distribution on the fracture property. We also compare our present numerical results with our experimental results. Section 5 summarises the major conclusions that can be drawn from this study.


## 2. Governing Equations

### 2.1. Elastic Dynamic Responses of the Solid Medium

The boundary of a bounded domain,Ω∈R2, is partitioned into three parts: the displacement boundary Γu, the traction boundary Γt, and the crack boundary Γc that is traction-free. The elastodynamic basic equation is expressed as(1)∇·σ+b=ρu¨inΩ,ε=∇suinΩ,σ=D:εinΩwith the following boundary and initial conditions:(2)ux,t=u-x,tonΓt,σ·n=t-onΓu,σ·n=0onΓc,ux,t=0=u-0,u˙x,t=0=u˙-0,where σ is the Cauchy stress tensor, b is the body force vector, ε is the strain tensor, ρ is the material density, u¨ is the acceleration field vector, ∇s is the symmetric part of the gradient operator, u is the displacement field vector, D is the constitutive matrix, n is the unit outward normal vector to the crack surface, u- is the prescribed displacement, t- is the external traction vector, u-0 is the initial displacement vector, and u˙-0 is the initial velocity vector.


### 2.2. Fluid Flow Pressure within the Crack

In this paper, we model the fluid flow within the crack as one-dimensional flow. Assume that the flow is laminar and the fluid is incompressible. But here we account for the time-dependent rate of crack opening.The conservation of the incompressible fluid in the fracture can be expressed as [6](3)∇·q+w˙+g=0,where ∇ is the divergence operator defined inx direction; q is the fluid flux; w˙ is the time-dependent rate of crack opening, and w˙=∂w/∂t; and g is the fluid loss into the solid media, and here we ignore the fluid loss; that is, g=0.Additionally, Poiseuille’s law [29] gives the following expression:(4)q=-w312μ∇p,where w is the crack opening; μ is the fluid viscosity; and p is the fluid pressure.The pressure boundary conditions at the fluid injection point in the crack are(5)px=0,t=p0,where p0 is the pressure of the fluid injection point.In the fluid lag zone [30], (6)px,t<0lt≤x≤l.According to the fluid pressure continuity, lag condition (6) provides the net-pressure boundary condition at the fluid front x=lt for the fluid flow equations (3).


## 2.1. Elastic Dynamic Responses of the Solid Medium

The boundary of a bounded domain,Ω∈R2, is partitioned into three parts: the displacement boundary Γu, the traction boundary Γt, and the crack boundary Γc that is traction-free. The elastodynamic basic equation is expressed as(1)∇·σ+b=ρu¨inΩ,ε=∇suinΩ,σ=D:εinΩwith the following boundary and initial conditions:(2)ux,t=u-x,tonΓt,σ·n=t-onΓu,σ·n=0onΓc,ux,t=0=u-0,u˙x,t=0=u˙-0,where σ is the Cauchy stress tensor, b is the body force vector, ε is the strain tensor, ρ is the material density, u¨ is the acceleration field vector, ∇s is the symmetric part of the gradient operator, u is the displacement field vector, D is the constitutive matrix, n is the unit outward normal vector to the crack surface, u- is the prescribed displacement, t- is the external traction vector, u-0 is the initial displacement vector, and u˙-0 is the initial velocity vector.


## 2.2. Fluid Flow Pressure within the Crack

In this paper, we model the fluid flow within the crack as one-dimensional flow. Assume that the flow is laminar and the fluid is incompressible. But here we account for the time-dependent rate of crack opening.The conservation of the incompressible fluid in the fracture can be expressed as [6](3)∇·q+w˙+g=0,where ∇ is the divergence operator defined inx direction; q is the fluid flux; w˙ is the time-dependent rate of crack opening, and w˙=∂w/∂t; and g is the fluid loss into the solid media, and here we ignore the fluid loss; that is, g=0.Additionally, Poiseuille’s law [29] gives the following expression:(4)q=-w312μ∇p,where w is the crack opening; μ is the fluid viscosity; and p is the fluid pressure.The pressure boundary conditions at the fluid injection point in the crack are(5)px=0,t=p0,where p0 is the pressure of the fluid injection point.In the fluid lag zone [30], (6)px,t<0lt≤x≤l.According to the fluid pressure continuity, lag condition (6) provides the net-pressure boundary condition at the fluid front x=lt for the fluid flow equations (3).


## 3. Numerical Approximation

### 3.1. XFEM for Dynamic Problems

#### 3.1.1. XFEM Approximation

The XFEM approximation for 2D cracked domains can be written as(7)uhx=∑i∈INixui+∑i∈Iabs∗Ni∗xHx-Hxiai+∑i∈Ibr∗Ni∗x∑j=14Fjx-Fjxibij,where Ni(x) is the standard finite element shape function of node i; ui is the unknown of the standard finite element part at node i; I is the set of all nodes in the domain; and Ni∗(x) is the partition of unity functions, and the function can hold the same form with the standard finite element shape function but is not necessary; ai and bij are the nodal enriched degree of freedom; Iabs∗ and Ibr∗ are the set of enrichment nodes shown in Figure 1, and Iabs∗,Ibr∗⊂I.Figure 1
A strategy for enriched elements and nodes.For these elements that are cut completely by a crack, the nodes of these elements that are the nodal subsetIabs∗ are enriched by Heaviside function H(x). The definition of Heaviside function H(x) follows:(8)Hx=+1,x-x∗·n>0-1,x-x∗·n<0,where x∗ is the projection of a point x on the crack surface; n is the unit outward normal to the crack surface.For these elements that are cut partially by a crack, the nodes of these elements that are the nodal subsetIbr∗ are enriched by the crack tip enrichment function Fj(x). The definition of the crack tip enrichment function Fj(x) follows:(9)Fj=1,2,3,4r,θ=rsin⁡θ2,rcos⁡θ2,rsin⁡θsin⁡θ2,rsin⁡θcos⁡θ2,where r and θ are the local crack tip coordinate system.


#### 3.1.2. Discrete Equations

By the principle of virtual work, the following discrete equations can be obtained:(10)Mu¨h+Kuh=f,where K(M) is the global stiffness (mass) matrix assembled by the element stiffness (mass) matrix; f is the global external force vector; uh and u¨h denote the vector of nodal parameters (which include the classic degrees of freedom, u, and the enrichment degrees of freedom, a, b) and its second time derivative, respectively; and(11)uh=u,a,bT,u¨h=u¨,a¨,b¨T.The element stiffness matrix is expressed by(12)ke=kuukuakubkaukaakabkbukbakbb,where(13)krs=∫ΩeBrTDBsdΩ=∫ΩeB1rB2rB3rB4rTDB1sB2sB3sB4sdΩr,s=u,a,b,Biu=∂Ni∂x0∂Ni∂y0∂Ni∂y∂Ni∂xT,i=1,2,3,4,Bia=∂Ni∗H⌢∂x0∂Ni∗H⌢∂y0∂Ni∗H⌢∂y∂Ni∗H⌢∂xT,H⌢=Hx-Hxi,i=1,2,3,4,Bib=Bib1Bib2Bib3Bib4,Bibj=∂Ni∗F⌢j∂x0∂Ni∗F⌢j∂y0∂Ni∗F⌢j∂y∂Ni∗F⌢j∂xT,F⌢j=Fjx-Fjxi,i,j=1,2,3,4.The element mass matrix is expressed by(14)me=muumuamubmaumaamabmbumbambb,where(15)muu=∫ΩeρNTNdΩ,mua=mauT=∫ΩeρNTNH⌢dΩ,mjub=mjbuT=∫ΩeρNTNF⌢jdΩ,maa=∫ΩeρNH⌢TNH⌢dΩ,mjab=mjbaT=∫ΩeρNH⌢TNF⌢jdΩ,mjkbb=∫ΩeρNF⌢jTNF⌢kdΩ,j,k=1,2,3,4,N=N1N2N3N4,Ni=Ni00Ni,i=1,2,3,4.The element external force vector is(16)fe=fufafbT,where(17)fu=∫ΩeNTbdΩ+∫ΓteNTt-dΓ,fa=∫ΩeN∗H⌢TbdΩ+∫ΓteN∗H⌢Tt-dΓ+∫ΓceN∗H⌢Tp·ndΓ,fbj=∫ΩeN∗F⌢jTbdΩ+∫ΓteN∗F⌢jTt-dΓ+∫ΓceN∗F⌢jTp·ndΓ,j=1,2,3,4.


#### 3.1.3. Time Integration Schemes

The following time-integration scheme is used in dynamic analysis. Equation (10) for a specific time t+Δt is expressed as(18)Mu¨t+Δt+Kut+Δt=ft+Δt,u˙t+Δt=u˙t+1-γu¨tΔt+γu¨t+ΔtΔt,ut+Δt=ut+u˙tΔt+12-βu¨tΔt2+βu¨t+ΔtΔt2,where ut, u˙t, and u¨t are the displacement, velocity, and acceleration vectors at time t, respectively; Δt is the time step; γ and β are parameters that can be determined to obtain integration accuracy and stability, with(19)β=141-α2,γ=12-α,-13≤α≤0.Here, referring to the software ABAQUS (ABAQUS Theory Manual, Version 6.9), we set parameterα=-0.05 to remove the slight high frequency noise in the solution without having any significant effect on the meaningful, lower frequency response.The following steps describe the prescribe integration method procedure, while neglecting the damping effects.(I) Initial Calculations
(i)
Form stiffness matrixK, and mass matrix M.
(ii)
Give the initial displacement vectoru0 and the initial velocity vector u˙0. Then, calculate the initial acceleration vector u¨0 by the equilibrium equation:(20)Mu¨0+Ku0=f0.
(iii)
Select a time stepΔt and the parameters β and γ. Here, β=0.275625 and γ=0.55 are used. Calculate integration constants:(21)c0=1βΔt2,c1=γβΔt,c2=1βΔt,c3=12β-1,c4=γβ-1,c5=Δtγ2β-1,c6=Δt1-γ,c7=γΔt.
(iv)
Form the effective stiffness matrixK~: (22)K~=K+c0M.(II) For Each Time Step
(i)
Calculate effective loadsf~t+Δt at time t+Δt: (23)f~t+Δt=ft+Δt+Mc0ut+c2u˙t+c3u¨t.
(ii)
Solve for the displacement vectorut+Δt at time t+Δt: (24)K~ut+Δt=f~t+Δt.
(iii)
Calculate the acceleration vectoru¨t+Δt and the velocity vector u˙t+Δt at time t+Δt:(25)u¨t+Δt=c0ut+Δt-ut-c2u˙t-c3u¨t,u˙t+Δt=u˙t+c6u¨t+c7u¨t+Δt.


#### 3.1.4. Integration Schemes at the Discontinuities

For these elements partitioned by a crack, the ordinary Gauss quadrature rules cannot accurately calculate the integration of enrichment function. An alternative method that is dividing the enrichment element into a set of subpolygons usually needs to be used [31]. Additionally, some simplified numerical integration methods also had been proposed in literatures. Ventura [32] conducted an important first attempt to simplify numerical integration. His work is based on replacing nonpolynomial functions by “equivalent” polynomials. However, the proposed method is exact for triangular and tetrahedral elements, but for quadrilateral elements, when the opposite sides are not parallel, additional approximation is introduced. Another method is strain smoothing [33]. In strain smoothing, the surface integration is transformed into equivalent boundary integration by use of the Green-Ostrogradsky theorem. Natarajan et al. [34] used the new numerical integration proposed for arbitrary polygons [35] to integrate the discontinuous and singular integrands appearing in the XFEM stiffness matrix. In this paper, the method subdividing the element into subquads is used. For these elements partitioned completely or partially by a crack, the method subdividing these elements into subquads is shown in Figure 2.Figure 2
Element partitioning method for these elements containing a discontinuous interface.
(a)
Elements partitioned completely by a straight crack
(b)
Elements partitioned completely by a broken crack
(c)
Elements partitioned partially by a straight crackTo solve the element stiffness or mass matrix of these enrichment elements, each subquad element is, respectively, transferred into the standard element(-1,1)×(-1,1) by the method of the coordinate transformation. The Gauss integration points are distributed into each subquad. To improve the accuracy of crack tip integration, 15 Gauss integration points are distributed into each subquad for elements containing crack tip. The numerical integration is firstly performed in each subquad element domain, and then the element stiffness or mass matrix of the enrichment element can be obtained by assembling the numerical integration results of each subquad element. It is worthwhile pointing out that these subquads are only necessary for integration purposes. They do not provide additional degrees of freedom for the global stiffness and mass matrix.


#### 3.1.5. Interaction Integral for Computing Stress Intensity Factors

Take field 1,σij(1),εij(1),ui(1), for the actual field, and field 2, σij(2),εij(2),ui(2), for the auxiliary field. The actual field is obtained from numerical solutions computed by using XFEM, and the auxiliary field refers to the asymptotic results of linear elastic fracture dynamics [36]. The interaction integral equation which is used to evaluate the stress intensity factors is as follows:(26)I1,2=∫Aσij1ui,12+σij2ui,11-12σik1εik2δ1jq,jdΩ+∫Γd+∪Γd-12σik1εik2δ1j-σij1ui,12-σij2ui,11qnjdΓ.The second term in (26) denotes the contribution of traction along crack interface. As shown in Figure 3, A denotes the circle domain with centre at the crack tip and the radius R; Γd+∪Γd- consists of a interface starting from the external integration radius to crack tip in a two-way manner. R is defined as(27)R=rkhe,where he is the crack-tip element size; rk is a user-specified scalar multiple; q is the weight function; q=1 if the node lies in A; and q=0 if the node lies outside of A or lies on the boundary of A. The weight function q in the interior of an element is obtained by the interpolation of the nodal value:(28)q=∑i=14Niqi.Figure 3
Elements selection for the interaction integral near the crack tip.Additionally, the interaction integral relates to the stress intensity factors through the relation:(29)M1,2=2E∗KIKIaux+KIIKIIaux,where KIaux and KIIaux are the local auxiliary stress intensity factors for the auxiliary fields, respectively; and the definition of E∗ is(30)E∗=Eplane  stressE1-ν2plane  strain.By settingKIaux=1 and KIIaux=0 as well as M1,2=M1(1,2), we obtain the expression of KI as follows:(31)KI=E∗M11,22.Similarly, we obtain the equality(32)KII=E∗M21,22by setting KIaux=0 and KIIaux=1 and M1,2=M2(1,2).


#### 3.1.6. Crack Propagation Criteria

The maximum circumferential stress criterion [37] is used to determine the crack growth direction. Once KI and KII are calculated, the criterion gives the following crack growth direction:(33)θc=2tan-1⁡14KIKII±KIKII2+8,where θc is the crack growth angle in the local crack-tip coordinate system. If KII=0, then θc=0. It should also be noted that if KII>0, the crack growth angle θc<0, and if KII<0, then θc>0. Sukumar and Prévost (2003) [38] proposed a computationally more amenable expression for θc:(34)θc=2tan-1⁡-2KII/KI1+1+8KII/KI2.The equivalent stress intensity factor then follows:(35)Ke=cos⁡θc2KIcos2⁡θc2-1.5KIIsin⁡θc.The double-K criterion [39] is used for determining crack propagation. The crack propagation processes can be expressed as follows:(i)
IfKe<KIcini, the crack does not propagate.(ii)
WhenKe=KIcini, the performed crack begins to crack initially.(iii)
WhenKIcini<Ke<KIcun, the propagating crack develops steadily.(iv)
WhenKe≥KIcun, the crack propagates unsteadily.K
I
c
i
n
i is the initiation toughness; KIcun is the unstable fracture toughness.


### 3.2. Finite Volume Methods of the Flow along a Crack

Considering (3) and (4) and ignoring the fluid loss into the solid media, the following expression can be obtained:(36)∂∂xk∂p∂x+∂w∂t=0,where (37)k=-w312μ.


#### 3.2.1. Discretised Equation

Equation (36) is one-dimensional steady state diffusion equation. Here, we discretise the equation by finite volume methods. The first step in the finite volume method is to divide the domain into discrete control volumes. As shown in Figure 4, we place a number of nodes in the fluid space between P and Q. The boundaries of control volumes are positioned midway between adjacent nodes. Thus, each node is surrounded by a control volume. The west side face of the control volume is referred to by “w” and the east side face of the control volume is referred to by “e.”Figure 4
Dividing fluid domain into discrete control volumes.Integrating the governing equation in (36) over a control volume, the following discretised equation can be obtained.(38)∫ΔV∂∂xk∂p∂xdV+∫ΔV∂w∂tdV=∮Ank∂p∂xdA+∫ΔV∂w∂tdV=kA∂p∂xe-kA∂p∂xw+∫ΔV∂w∂tdV=kA∂p∂xe-kA∂p∂xw+Δw-Δt·ΔV=0,where A is the cross-sectional area of the control volume face; ΔV is the volume; w- is the average crack opening width over a control volume.According to (36) and observing a typical finite volume element with node 2, we can obtain(39)kA∂p∂xe-kA∂p∂xw=-Δw-2Δt·ΔV2.To calculate gradients at the control volume faces, a linear approximate distribution of pressure is considered here. So (39) yields(40)keAep3-p2L2-kwAwp2-p1L1=-Δw-2Δt·ΔV2.kw, Aw, ke, and Ae can be obtained by the way of linearly interpolated values, and(41)kw=δwk1+1-δwk2,Aw=δwA1+1-δwA2,δw=Lw2L1,ke=δek2+1-δek3,Ae=δeA2+1-δeA3,δe=Le3L2.Equation (40) can be rewritten as(42)a21p1+a22p2+a23p3=b2,where(43)a21=kwAwL1,a22=-kwAwL1-keAeL2,a23=keAeL2,b2=-Δw-2Δt·ΔV2.


#### 3.2.2. Boundary Conditions

(I) Boundary Conditions for Fluid Injection Point. As shown in Figure 5, integration of (36) over the control volume surrounding node 1 gives(44)keAep2-p1L1-kwAwp1-pwL0=-Δw-1Δt·ΔV1.Equation (44) can be rewritten as(45)a11p1+a12p2=b1,where(46)a11=-keAeL1-kwAwL0,a12=keAeL1,b1=-Δw-1Δt·ΔV1-kwAwL0pw,ke=δek1+1-δek2,Ae=δeA1+1-δeA2.kw and Aw can be calculated directly according to crack opening displacement, and pw is the pressure of the fluid injection point.Figure 5
Boundary conditions for fluid injection point.(II) Boundary Conditions for Fluid Front. Similarly, as shown in Figure 6, integration of (36) over the control volume surrounding noden gives(47)keAept-pnLn-kwAwpn-pn-1Ln-1=-Δw-nΔt·ΔVn.Equation (47) can be rewritten as(48)ann-1pn-1+annpn=bn,where(49)ann-1=kwAwLn-1,ann=-kwAwLn-1-keAeLn,bn=-Δw-nΔt·ΔVn-keAeptLn,kw=δwkn-1+1-δwkn,Aw=δwAn-1+1-δwAn.ke and Ae can be calculated directly according to crack width at the fluid front, and pt is the pressure of the fluid front and pt=0.Figure 6
Boundary conditions for fluid front.By setting up discretised equations of the forms (42), (45), and (48), at each node, the fluid pressure distribution along a crack can be obtained.


### 3.3. Coupling Scheme

To correctly solve the system of equations, the elasticity, fluid flow, and fracture growth should be coupled together; see Figure7. The crack opening causes the fluid injecting into the solid medium. The injecting fluid produces fluid pressures along crack surface and affects the deformation of the solid medium and crack propagation. Due to the deformation of the solid medium, the crack opening changes again.Figure 7
The coupling of the elasticity, fluid flow, and fracture growth.According to (36), the relations between the fluid pressure and crack width are nonlinear. It is necessary to use an algorithm that can solve iteratively the fluid pressure and crack width. The staggered Newton algorithm [40] is used here. The iteration processes for solving fluid pressure at time t are stated as follows:(i)
Assume that the fluid pressurepk has been computed at thekth step iteration. Compute dynamic responses of solid medium under dynamic loads and fluid pressure by using XFEM. According to dynamic responses, we can obtain the crack width wk at the kth step iteration.(ii)
According to the known crack widthwk at thekth step iteration, compute the fluid pressure pk+1/2 by using FVM.(iii)
Compute the fluid pressurepk+1 at the (k+1)th step iteration by the flowing expression:(50)pk+1=1-αpk+αpk+1/2.(iv)
Compute dynamic responses of solid medium under dynamic loads and fluid pressurepk+1 by using XFEM. According to dynamic responses, we can obtain the crack width wk+1/2.(v)
Compute the crack widthwk+1 at the (k+1)th step iteration by the flowing expression:(51)wk+1=1-αwk+αwk+1/2.(vi)
Given toleranceε=0.01, compute the following expression:(52)δp=∑i=1npk+1-pk∑i=1npk+1.If δp≤ε, finish iteration, or else go to (ii).


## 3.1. XFEM for Dynamic Problems

### 3.1.1. XFEM Approximation

The XFEM approximation for 2D cracked domains can be written as(7)uhx=∑i∈INixui+∑i∈Iabs∗Ni∗xHx-Hxiai+∑i∈Ibr∗Ni∗x∑j=14Fjx-Fjxibij,where Ni(x) is the standard finite element shape function of node i; ui is the unknown of the standard finite element part at node i; I is the set of all nodes in the domain; and Ni∗(x) is the partition of unity functions, and the function can hold the same form with the standard finite element shape function but is not necessary; ai and bij are the nodal enriched degree of freedom; Iabs∗ and Ibr∗ are the set of enrichment nodes shown in Figure 1, and Iabs∗,Ibr∗⊂I.Figure 1
A strategy for enriched elements and nodes.For these elements that are cut completely by a crack, the nodes of these elements that are the nodal subsetIabs∗ are enriched by Heaviside function H(x). The definition of Heaviside function H(x) follows:(8)Hx=+1,x-x∗·n>0-1,x-x∗·n<0,where x∗ is the projection of a point x on the crack surface; n is the unit outward normal to the crack surface.For these elements that are cut partially by a crack, the nodes of these elements that are the nodal subsetIbr∗ are enriched by the crack tip enrichment function Fj(x). The definition of the crack tip enrichment function Fj(x) follows:(9)Fj=1,2,3,4r,θ=rsin⁡θ2,rcos⁡θ2,rsin⁡θsin⁡θ2,rsin⁡θcos⁡θ2,where r and θ are the local crack tip coordinate system.


### 3.1.2. Discrete Equations

By the principle of virtual work, the following discrete equations can be obtained:(10)Mu¨h+Kuh=f,where K(M) is the global stiffness (mass) matrix assembled by the element stiffness (mass) matrix; f is the global external force vector; uh and u¨h denote the vector of nodal parameters (which include the classic degrees of freedom, u, and the enrichment degrees of freedom, a, b) and its second time derivative, respectively; and(11)uh=u,a,bT,u¨h=u¨,a¨,b¨T.The element stiffness matrix is expressed by(12)ke=kuukuakubkaukaakabkbukbakbb,where(13)krs=∫ΩeBrTDBsdΩ=∫ΩeB1rB2rB3rB4rTDB1sB2sB3sB4sdΩr,s=u,a,b,Biu=∂Ni∂x0∂Ni∂y0∂Ni∂y∂Ni∂xT,i=1,2,3,4,Bia=∂Ni∗H⌢∂x0∂Ni∗H⌢∂y0∂Ni∗H⌢∂y∂Ni∗H⌢∂xT,H⌢=Hx-Hxi,i=1,2,3,4,Bib=Bib1Bib2Bib3Bib4,Bibj=∂Ni∗F⌢j∂x0∂Ni∗F⌢j∂y0∂Ni∗F⌢j∂y∂Ni∗F⌢j∂xT,F⌢j=Fjx-Fjxi,i,j=1,2,3,4.The element mass matrix is expressed by(14)me=muumuamubmaumaamabmbumbambb,where(15)muu=∫ΩeρNTNdΩ,mua=mauT=∫ΩeρNTNH⌢dΩ,mjub=mjbuT=∫ΩeρNTNF⌢jdΩ,maa=∫ΩeρNH⌢TNH⌢dΩ,mjab=mjbaT=∫ΩeρNH⌢TNF⌢jdΩ,mjkbb=∫ΩeρNF⌢jTNF⌢kdΩ,j,k=1,2,3,4,N=N1N2N3N4,Ni=Ni00Ni,i=1,2,3,4.The element external force vector is(16)fe=fufafbT,where(17)fu=∫ΩeNTbdΩ+∫ΓteNTt-dΓ,fa=∫ΩeN∗H⌢TbdΩ+∫ΓteN∗H⌢Tt-dΓ+∫ΓceN∗H⌢Tp·ndΓ,fbj=∫ΩeN∗F⌢jTbdΩ+∫ΓteN∗F⌢jTt-dΓ+∫ΓceN∗F⌢jTp·ndΓ,j=1,2,3,4.


### 3.1.3. Time Integration Schemes

The following time-integration scheme is used in dynamic analysis. Equation (10) for a specific time t+Δt is expressed as(18)Mu¨t+Δt+Kut+Δt=ft+Δt,u˙t+Δt=u˙t+1-γu¨tΔt+γu¨t+ΔtΔt,ut+Δt=ut+u˙tΔt+12-βu¨tΔt2+βu¨t+ΔtΔt2,where ut, u˙t, and u¨t are the displacement, velocity, and acceleration vectors at time t, respectively; Δt is the time step; γ and β are parameters that can be determined to obtain integration accuracy and stability, with(19)β=141-α2,γ=12-α,-13≤α≤0.Here, referring to the software ABAQUS (ABAQUS Theory Manual, Version 6.9), we set parameterα=-0.05 to remove the slight high frequency noise in the solution without having any significant effect on the meaningful, lower frequency response.The following steps describe the prescribe integration method procedure, while neglecting the damping effects.(I) Initial Calculations
(i)
Form stiffness matrixK, and mass matrix M.
(ii)
Give the initial displacement vectoru0 and the initial velocity vector u˙0. Then, calculate the initial acceleration vector u¨0 by the equilibrium equation:(20)Mu¨0+Ku0=f0.
(iii)
Select a time stepΔt and the parameters β and γ. Here, β=0.275625 and γ=0.55 are used. Calculate integration constants:(21)c0=1βΔt2,c1=γβΔt,c2=1βΔt,c3=12β-1,c4=γβ-1,c5=Δtγ2β-1,c6=Δt1-γ,c7=γΔt.
(iv)
Form the effective stiffness matrixK~: (22)K~=K+c0M.(II) For Each Time Step
(i)
Calculate effective loadsf~t+Δt at time t+Δt: (23)f~t+Δt=ft+Δt+Mc0ut+c2u˙t+c3u¨t.
(ii)
Solve for the displacement vectorut+Δt at time t+Δt: (24)K~ut+Δt=f~t+Δt.
(iii)
Calculate the acceleration vectoru¨t+Δt and the velocity vector u˙t+Δt at time t+Δt:(25)u¨t+Δt=c0ut+Δt-ut-c2u˙t-c3u¨t,u˙t+Δt=u˙t+c6u¨t+c7u¨t+Δt.


### 3.1.4. Integration Schemes at the Discontinuities

For these elements partitioned by a crack, the ordinary Gauss quadrature rules cannot accurately calculate the integration of enrichment function. An alternative method that is dividing the enrichment element into a set of subpolygons usually needs to be used [31]. Additionally, some simplified numerical integration methods also had been proposed in literatures. Ventura [32] conducted an important first attempt to simplify numerical integration. His work is based on replacing nonpolynomial functions by “equivalent” polynomials. However, the proposed method is exact for triangular and tetrahedral elements, but for quadrilateral elements, when the opposite sides are not parallel, additional approximation is introduced. Another method is strain smoothing [33]. In strain smoothing, the surface integration is transformed into equivalent boundary integration by use of the Green-Ostrogradsky theorem. Natarajan et al. [34] used the new numerical integration proposed for arbitrary polygons [35] to integrate the discontinuous and singular integrands appearing in the XFEM stiffness matrix. In this paper, the method subdividing the element into subquads is used. For these elements partitioned completely or partially by a crack, the method subdividing these elements into subquads is shown in Figure 2.Figure 2
Element partitioning method for these elements containing a discontinuous interface.
(a)
Elements partitioned completely by a straight crack
(b)
Elements partitioned completely by a broken crack
(c)
Elements partitioned partially by a straight crackTo solve the element stiffness or mass matrix of these enrichment elements, each subquad element is, respectively, transferred into the standard element(-1,1)×(-1,1) by the method of the coordinate transformation. The Gauss integration points are distributed into each subquad. To improve the accuracy of crack tip integration, 15 Gauss integration points are distributed into each subquad for elements containing crack tip. The numerical integration is firstly performed in each subquad element domain, and then the element stiffness or mass matrix of the enrichment element can be obtained by assembling the numerical integration results of each subquad element. It is worthwhile pointing out that these subquads are only necessary for integration purposes. They do not provide additional degrees of freedom for the global stiffness and mass matrix.


### 3.1.5. Interaction Integral for Computing Stress Intensity Factors

Take field 1,σij(1),εij(1),ui(1), for the actual field, and field 2, σij(2),εij(2),ui(2), for the auxiliary field. The actual field is obtained from numerical solutions computed by using XFEM, and the auxiliary field refers to the asymptotic results of linear elastic fracture dynamics [36]. The interaction integral equation which is used to evaluate the stress intensity factors is as follows:(26)I1,2=∫Aσij1ui,12+σij2ui,11-12σik1εik2δ1jq,jdΩ+∫Γd+∪Γd-12σik1εik2δ1j-σij1ui,12-σij2ui,11qnjdΓ.The second term in (26) denotes the contribution of traction along crack interface. As shown in Figure 3, A denotes the circle domain with centre at the crack tip and the radius R; Γd+∪Γd- consists of a interface starting from the external integration radius to crack tip in a two-way manner. R is defined as(27)R=rkhe,where he is the crack-tip element size; rk is a user-specified scalar multiple; q is the weight function; q=1 if the node lies in A; and q=0 if the node lies outside of A or lies on the boundary of A. The weight function q in the interior of an element is obtained by the interpolation of the nodal value:(28)q=∑i=14Niqi.Figure 3
Elements selection for the interaction integral near the crack tip.Additionally, the interaction integral relates to the stress intensity factors through the relation:(29)M1,2=2E∗KIKIaux+KIIKIIaux,where KIaux and KIIaux are the local auxiliary stress intensity factors for the auxiliary fields, respectively; and the definition of E∗ is(30)E∗=Eplane  stressE1-ν2plane  strain.By settingKIaux=1 and KIIaux=0 as well as M1,2=M1(1,2), we obtain the expression of KI as follows:(31)KI=E∗M11,22.Similarly, we obtain the equality(32)KII=E∗M21,22by setting KIaux=0 and KIIaux=1 and M1,2=M2(1,2).


### 3.1.6. Crack Propagation Criteria

The maximum circumferential stress criterion [37] is used to determine the crack growth direction. Once KI and KII are calculated, the criterion gives the following crack growth direction:(33)θc=2tan-1⁡14KIKII±KIKII2+8,where θc is the crack growth angle in the local crack-tip coordinate system. If KII=0, then θc=0. It should also be noted that if KII>0, the crack growth angle θc<0, and if KII<0, then θc>0. Sukumar and Prévost (2003) [38] proposed a computationally more amenable expression for θc:(34)θc=2tan-1⁡-2KII/KI1+1+8KII/KI2.The equivalent stress intensity factor then follows:(35)Ke=cos⁡θc2KIcos2⁡θc2-1.5KIIsin⁡θc.The double-K criterion [39] is used for determining crack propagation. The crack propagation processes can be expressed as follows:(i)
IfKe<KIcini, the crack does not propagate.(ii)
WhenKe=KIcini, the performed crack begins to crack initially.(iii)
WhenKIcini<Ke<KIcun, the propagating crack develops steadily.(iv)
WhenKe≥KIcun, the crack propagates unsteadily.K
I
c
i
n
i is the initiation toughness; KIcun is the unstable fracture toughness.


## 3.1.1. XFEM Approximation

The XFEM approximation for 2D cracked domains can be written as(7)uhx=∑i∈INixui+∑i∈Iabs∗Ni∗xHx-Hxiai+∑i∈Ibr∗Ni∗x∑j=14Fjx-Fjxibij,where Ni(x) is the standard finite element shape function of node i; ui is the unknown of the standard finite element part at node i; I is the set of all nodes in the domain; and Ni∗(x) is the partition of unity functions, and the function can hold the same form with the standard finite element shape function but is not necessary; ai and bij are the nodal enriched degree of freedom; Iabs∗ and Ibr∗ are the set of enrichment nodes shown in Figure 1, and Iabs∗,Ibr∗⊂I.Figure 1
A strategy for enriched elements and nodes.For these elements that are cut completely by a crack, the nodes of these elements that are the nodal subsetIabs∗ are enriched by Heaviside function H(x). The definition of Heaviside function H(x) follows:(8)Hx=+1,x-x∗·n>0-1,x-x∗·n<0,where x∗ is the projection of a point x on the crack surface; n is the unit outward normal to the crack surface.For these elements that are cut partially by a crack, the nodes of these elements that are the nodal subsetIbr∗ are enriched by the crack tip enrichment function Fj(x). The definition of the crack tip enrichment function Fj(x) follows:(9)Fj=1,2,3,4r,θ=rsin⁡θ2,rcos⁡θ2,rsin⁡θsin⁡θ2,rsin⁡θcos⁡θ2,where r and θ are the local crack tip coordinate system.


## 3.1.2. Discrete Equations

By the principle of virtual work, the following discrete equations can be obtained:(10)Mu¨h+Kuh=f,where K(M) is the global stiffness (mass) matrix assembled by the element stiffness (mass) matrix; f is the global external force vector; uh and u¨h denote the vector of nodal parameters (which include the classic degrees of freedom, u, and the enrichment degrees of freedom, a, b) and its second time derivative, respectively; and(11)uh=u,a,bT,u¨h=u¨,a¨,b¨T.The element stiffness matrix is expressed by(12)ke=kuukuakubkaukaakabkbukbakbb,where(13)krs=∫ΩeBrTDBsdΩ=∫ΩeB1rB2rB3rB4rTDB1sB2sB3sB4sdΩr,s=u,a,b,Biu=∂Ni∂x0∂Ni∂y0∂Ni∂y∂Ni∂xT,i=1,2,3,4,Bia=∂Ni∗H⌢∂x0∂Ni∗H⌢∂y0∂Ni∗H⌢∂y∂Ni∗H⌢∂xT,H⌢=Hx-Hxi,i=1,2,3,4,Bib=Bib1Bib2Bib3Bib4,Bibj=∂Ni∗F⌢j∂x0∂Ni∗F⌢j∂y0∂Ni∗F⌢j∂y∂Ni∗F⌢j∂xT,F⌢j=Fjx-Fjxi,i,j=1,2,3,4.The element mass matrix is expressed by(14)me=muumuamubmaumaamabmbumbambb,where(15)muu=∫ΩeρNTNdΩ,mua=mauT=∫ΩeρNTNH⌢dΩ,mjub=mjbuT=∫ΩeρNTNF⌢jdΩ,maa=∫ΩeρNH⌢TNH⌢dΩ,mjab=mjbaT=∫ΩeρNH⌢TNF⌢jdΩ,mjkbb=∫ΩeρNF⌢jTNF⌢kdΩ,j,k=1,2,3,4,N=N1N2N3N4,Ni=Ni00Ni,i=1,2,3,4.The element external force vector is(16)fe=fufafbT,where(17)fu=∫ΩeNTbdΩ+∫ΓteNTt-dΓ,fa=∫ΩeN∗H⌢TbdΩ+∫ΓteN∗H⌢Tt-dΓ+∫ΓceN∗H⌢Tp·ndΓ,fbj=∫ΩeN∗F⌢jTbdΩ+∫ΓteN∗F⌢jTt-dΓ+∫ΓceN∗F⌢jTp·ndΓ,j=1,2,3,4.


## 3.1.3. Time Integration Schemes

The following time-integration scheme is used in dynamic analysis. Equation (10) for a specific time t+Δt is expressed as(18)Mu¨t+Δt+Kut+Δt=ft+Δt,u˙t+Δt=u˙t+1-γu¨tΔt+γu¨t+ΔtΔt,ut+Δt=ut+u˙tΔt+12-βu¨tΔt2+βu¨t+ΔtΔt2,where ut, u˙t, and u¨t are the displacement, velocity, and acceleration vectors at time t, respectively; Δt is the time step; γ and β are parameters that can be determined to obtain integration accuracy and stability, with(19)β=141-α2,γ=12-α,-13≤α≤0.Here, referring to the software ABAQUS (ABAQUS Theory Manual, Version 6.9), we set parameterα=-0.05 to remove the slight high frequency noise in the solution without having any significant effect on the meaningful, lower frequency response.The following steps describe the prescribe integration method procedure, while neglecting the damping effects.(I) Initial Calculations
(i)
Form stiffness matrixK, and mass matrix M.
(ii)
Give the initial displacement vectoru0 and the initial velocity vector u˙0. Then, calculate the initial acceleration vector u¨0 by the equilibrium equation:(20)Mu¨0+Ku0=f0.
(iii)
Select a time stepΔt and the parameters β and γ. Here, β=0.275625 and γ=0.55 are used. Calculate integration constants:(21)c0=1βΔt2,c1=γβΔt,c2=1βΔt,c3=12β-1,c4=γβ-1,c5=Δtγ2β-1,c6=Δt1-γ,c7=γΔt.
(iv)
Form the effective stiffness matrixK~: (22)K~=K+c0M.(II) For Each Time Step
(i)
Calculate effective loadsf~t+Δt at time t+Δt: (23)f~t+Δt=ft+Δt+Mc0ut+c2u˙t+c3u¨t.
(ii)
Solve for the displacement vectorut+Δt at time t+Δt: (24)K~ut+Δt=f~t+Δt.
(iii)
Calculate the acceleration vectoru¨t+Δt and the velocity vector u˙t+Δt at time t+Δt:(25)u¨t+Δt=c0ut+Δt-ut-c2u˙t-c3u¨t,u˙t+Δt=u˙t+c6u¨t+c7u¨t+Δt.


## 3.1.4. Integration Schemes at the Discontinuities

For these elements partitioned by a crack, the ordinary Gauss quadrature rules cannot accurately calculate the integration of enrichment function. An alternative method that is dividing the enrichment element into a set of subpolygons usually needs to be used [31]. Additionally, some simplified numerical integration methods also had been proposed in literatures. Ventura [32] conducted an important first attempt to simplify numerical integration. His work is based on replacing nonpolynomial functions by “equivalent” polynomials. However, the proposed method is exact for triangular and tetrahedral elements, but for quadrilateral elements, when the opposite sides are not parallel, additional approximation is introduced. Another method is strain smoothing [33]. In strain smoothing, the surface integration is transformed into equivalent boundary integration by use of the Green-Ostrogradsky theorem. Natarajan et al. [34] used the new numerical integration proposed for arbitrary polygons [35] to integrate the discontinuous and singular integrands appearing in the XFEM stiffness matrix. In this paper, the method subdividing the element into subquads is used. For these elements partitioned completely or partially by a crack, the method subdividing these elements into subquads is shown in Figure 2.Figure 2
Element partitioning method for these elements containing a discontinuous interface.
(a)
Elements partitioned completely by a straight crack
(b)
Elements partitioned completely by a broken crack
(c)
Elements partitioned partially by a straight crackTo solve the element stiffness or mass matrix of these enrichment elements, each subquad element is, respectively, transferred into the standard element(-1,1)×(-1,1) by the method of the coordinate transformation. The Gauss integration points are distributed into each subquad. To improve the accuracy of crack tip integration, 15 Gauss integration points are distributed into each subquad for elements containing crack tip. The numerical integration is firstly performed in each subquad element domain, and then the element stiffness or mass matrix of the enrichment element can be obtained by assembling the numerical integration results of each subquad element. It is worthwhile pointing out that these subquads are only necessary for integration purposes. They do not provide additional degrees of freedom for the global stiffness and mass matrix.


## 3.1.5. Interaction Integral for Computing Stress Intensity Factors

Take field 1,σij(1),εij(1),ui(1), for the actual field, and field 2, σij(2),εij(2),ui(2), for the auxiliary field. The actual field is obtained from numerical solutions computed by using XFEM, and the auxiliary field refers to the asymptotic results of linear elastic fracture dynamics [36]. The interaction integral equation which is used to evaluate the stress intensity factors is as follows:(26)I1,2=∫Aσij1ui,12+σij2ui,11-12σik1εik2δ1jq,jdΩ+∫Γd+∪Γd-12σik1εik2δ1j-σij1ui,12-σij2ui,11qnjdΓ.The second term in (26) denotes the contribution of traction along crack interface. As shown in Figure 3, A denotes the circle domain with centre at the crack tip and the radius R; Γd+∪Γd- consists of a interface starting from the external integration radius to crack tip in a two-way manner. R is defined as(27)R=rkhe,where he is the crack-tip element size; rk is a user-specified scalar multiple; q is the weight function; q=1 if the node lies in A; and q=0 if the node lies outside of A or lies on the boundary of A. The weight function q in the interior of an element is obtained by the interpolation of the nodal value:(28)q=∑i=14Niqi.Figure 3
Elements selection for the interaction integral near the crack tip.Additionally, the interaction integral relates to the stress intensity factors through the relation:(29)M1,2=2E∗KIKIaux+KIIKIIaux,where KIaux and KIIaux are the local auxiliary stress intensity factors for the auxiliary fields, respectively; and the definition of E∗ is(30)E∗=Eplane  stressE1-ν2plane  strain.By settingKIaux=1 and KIIaux=0 as well as M1,2=M1(1,2), we obtain the expression of KI as follows:(31)KI=E∗M11,22.Similarly, we obtain the equality(32)KII=E∗M21,22by setting KIaux=0 and KIIaux=1 and M1,2=M2(1,2).


## 3.1.6. Crack Propagation Criteria

The maximum circumferential stress criterion [37] is used to determine the crack growth direction. Once KI and KII are calculated, the criterion gives the following crack growth direction:(33)θc=2tan-1⁡14KIKII±KIKII2+8,where θc is the crack growth angle in the local crack-tip coordinate system. If KII=0, then θc=0. It should also be noted that if KII>0, the crack growth angle θc<0, and if KII<0, then θc>0. Sukumar and Prévost (2003) [38] proposed a computationally more amenable expression for θc:(34)θc=2tan-1⁡-2KII/KI1+1+8KII/KI2.The equivalent stress intensity factor then follows:(35)Ke=cos⁡θc2KIcos2⁡θc2-1.5KIIsin⁡θc.The double-K criterion [39] is used for determining crack propagation. The crack propagation processes can be expressed as follows:(i)
IfKe<KIcini, the crack does not propagate.(ii)
WhenKe=KIcini, the performed crack begins to crack initially.(iii)
WhenKIcini<Ke<KIcun, the propagating crack develops steadily.(iv)
WhenKe≥KIcun, the crack propagates unsteadily.K
I
c
i
n
i is the initiation toughness; KIcun is the unstable fracture toughness.


## 3.2. Finite Volume Methods of the Flow along a Crack

Considering (3) and (4) and ignoring the fluid loss into the solid media, the following expression can be obtained:(36)∂∂xk∂p∂x+∂w∂t=0,where (37)k=-w312μ.


### 3.2.1. Discretised Equation

Equation (36) is one-dimensional steady state diffusion equation. Here, we discretise the equation by finite volume methods. The first step in the finite volume method is to divide the domain into discrete control volumes. As shown in Figure 4, we place a number of nodes in the fluid space between P and Q. The boundaries of control volumes are positioned midway between adjacent nodes. Thus, each node is surrounded by a control volume. The west side face of the control volume is referred to by “w” and the east side face of the control volume is referred to by “e.”Figure 4
Dividing fluid domain into discrete control volumes.Integrating the governing equation in (36) over a control volume, the following discretised equation can be obtained.(38)∫ΔV∂∂xk∂p∂xdV+∫ΔV∂w∂tdV=∮Ank∂p∂xdA+∫ΔV∂w∂tdV=kA∂p∂xe-kA∂p∂xw+∫ΔV∂w∂tdV=kA∂p∂xe-kA∂p∂xw+Δw-Δt·ΔV=0,where A is the cross-sectional area of the control volume face; ΔV is the volume; w- is the average crack opening width over a control volume.According to (36) and observing a typical finite volume element with node 2, we can obtain(39)kA∂p∂xe-kA∂p∂xw=-Δw-2Δt·ΔV2.To calculate gradients at the control volume faces, a linear approximate distribution of pressure is considered here. So (39) yields(40)keAep3-p2L2-kwAwp2-p1L1=-Δw-2Δt·ΔV2.kw, Aw, ke, and Ae can be obtained by the way of linearly interpolated values, and(41)kw=δwk1+1-δwk2,Aw=δwA1+1-δwA2,δw=Lw2L1,ke=δek2+1-δek3,Ae=δeA2+1-δeA3,δe=Le3L2.Equation (40) can be rewritten as(42)a21p1+a22p2+a23p3=b2,where(43)a21=kwAwL1,a22=-kwAwL1-keAeL2,a23=keAeL2,b2=-Δw-2Δt·ΔV2.


### 3.2.2. Boundary Conditions

(I) Boundary Conditions for Fluid Injection Point. As shown in Figure 5, integration of (36) over the control volume surrounding node 1 gives(44)keAep2-p1L1-kwAwp1-pwL0=-Δw-1Δt·ΔV1.Equation (44) can be rewritten as(45)a11p1+a12p2=b1,where(46)a11=-keAeL1-kwAwL0,a12=keAeL1,b1=-Δw-1Δt·ΔV1-kwAwL0pw,ke=δek1+1-δek2,Ae=δeA1+1-δeA2.kw and Aw can be calculated directly according to crack opening displacement, and pw is the pressure of the fluid injection point.Figure 5
Boundary conditions for fluid injection point.(II) Boundary Conditions for Fluid Front. Similarly, as shown in Figure 6, integration of (36) over the control volume surrounding noden gives(47)keAept-pnLn-kwAwpn-pn-1Ln-1=-Δw-nΔt·ΔVn.Equation (47) can be rewritten as(48)ann-1pn-1+annpn=bn,where(49)ann-1=kwAwLn-1,ann=-kwAwLn-1-keAeLn,bn=-Δw-nΔt·ΔVn-keAeptLn,kw=δwkn-1+1-δwkn,Aw=δwAn-1+1-δwAn.ke and Ae can be calculated directly according to crack width at the fluid front, and pt is the pressure of the fluid front and pt=0.Figure 6
Boundary conditions for fluid front.By setting up discretised equations of the forms (42), (45), and (48), at each node, the fluid pressure distribution along a crack can be obtained.


## 3.2.1. Discretised Equation

Equation (36) is one-dimensional steady state diffusion equation. Here, we discretise the equation by finite volume methods. The first step in the finite volume method is to divide the domain into discrete control volumes. As shown in Figure 4, we place a number of nodes in the fluid space between P and Q. The boundaries of control volumes are positioned midway between adjacent nodes. Thus, each node is surrounded by a control volume. The west side face of the control volume is referred to by “w” and the east side face of the control volume is referred to by “e.”Figure 4
Dividing fluid domain into discrete control volumes.Integrating the governing equation in (36) over a control volume, the following discretised equation can be obtained.(38)∫ΔV∂∂xk∂p∂xdV+∫ΔV∂w∂tdV=∮Ank∂p∂xdA+∫ΔV∂w∂tdV=kA∂p∂xe-kA∂p∂xw+∫ΔV∂w∂tdV=kA∂p∂xe-kA∂p∂xw+Δw-Δt·ΔV=0,where A is the cross-sectional area of the control volume face; ΔV is the volume; w- is the average crack opening width over a control volume.According to (36) and observing a typical finite volume element with node 2, we can obtain(39)kA∂p∂xe-kA∂p∂xw=-Δw-2Δt·ΔV2.To calculate gradients at the control volume faces, a linear approximate distribution of pressure is considered here. So (39) yields(40)keAep3-p2L2-kwAwp2-p1L1=-Δw-2Δt·ΔV2.kw, Aw, ke, and Ae can be obtained by the way of linearly interpolated values, and(41)kw=δwk1+1-δwk2,Aw=δwA1+1-δwA2,δw=Lw2L1,ke=δek2+1-δek3,Ae=δeA2+1-δeA3,δe=Le3L2.Equation (40) can be rewritten as(42)a21p1+a22p2+a23p3=b2,where(43)a21=kwAwL1,a22=-kwAwL1-keAeL2,a23=keAeL2,b2=-Δw-2Δt·ΔV2.


## 3.2.2. Boundary Conditions

(I) Boundary Conditions for Fluid Injection Point. As shown in Figure 5, integration of (36) over the control volume surrounding node 1 gives(44)keAep2-p1L1-kwAwp1-pwL0=-Δw-1Δt·ΔV1.Equation (44) can be rewritten as(45)a11p1+a12p2=b1,where(46)a11=-keAeL1-kwAwL0,a12=keAeL1,b1=-Δw-1Δt·ΔV1-kwAwL0pw,ke=δek1+1-δek2,Ae=δeA1+1-δeA2.kw and Aw can be calculated directly according to crack opening displacement, and pw is the pressure of the fluid injection point.Figure 5
Boundary conditions for fluid injection point.(II) Boundary Conditions for Fluid Front. Similarly, as shown in Figure 6, integration of (36) over the control volume surrounding noden gives(47)keAept-pnLn-kwAwpn-pn-1Ln-1=-Δw-nΔt·ΔVn.Equation (47) can be rewritten as(48)ann-1pn-1+annpn=bn,where(49)ann-1=kwAwLn-1,ann=-kwAwLn-1-keAeLn,bn=-Δw-nΔt·ΔVn-keAeptLn,kw=δwkn-1+1-δwkn,Aw=δwAn-1+1-δwAn.ke and Ae can be calculated directly according to crack width at the fluid front, and pt is the pressure of the fluid front and pt=0.Figure 6
Boundary conditions for fluid front.By setting up discretised equations of the forms (42), (45), and (48), at each node, the fluid pressure distribution along a crack can be obtained.


## 3.3. Coupling Scheme

To correctly solve the system of equations, the elasticity, fluid flow, and fracture growth should be coupled together; see Figure7. The crack opening causes the fluid injecting into the solid medium. The injecting fluid produces fluid pressures along crack surface and affects the deformation of the solid medium and crack propagation. Due to the deformation of the solid medium, the crack opening changes again.Figure 7
The coupling of the elasticity, fluid flow, and fracture growth.According to (36), the relations between the fluid pressure and crack width are nonlinear. It is necessary to use an algorithm that can solve iteratively the fluid pressure and crack width. The staggered Newton algorithm [40] is used here. The iteration processes for solving fluid pressure at time t are stated as follows:(i)
Assume that the fluid pressurepk has been computed at thekth step iteration. Compute dynamic responses of solid medium under dynamic loads and fluid pressure by using XFEM. According to dynamic responses, we can obtain the crack width wk at the kth step iteration.(ii)
According to the known crack widthwk at thekth step iteration, compute the fluid pressure pk+1/2 by using FVM.(iii)
Compute the fluid pressurepk+1 at the (k+1)th step iteration by the flowing expression:(50)pk+1=1-αpk+αpk+1/2.(iv)
Compute dynamic responses of solid medium under dynamic loads and fluid pressurepk+1 by using XFEM. According to dynamic responses, we can obtain the crack width wk+1/2.(v)
Compute the crack widthwk+1 at the (k+1)th step iteration by the flowing expression:(51)wk+1=1-αwk+αwk+1/2.(vi)
Given toleranceε=0.01, compute the following expression:(52)δp=∑i=1npk+1-pk∑i=1npk+1.If δp≤ε, finish iteration, or else go to (ii).


## 4. Numerical Example

As shown in Figure8, a notched cubic concrete specimen with the dimension of 200 mm × 200 mm × 200 mm subjected to a splitting force was considered in this example. Here, the splitting force was applied on the iron. The iron was fastened on the concrete specimen. The material properties of concrete were Young modulus E = 36 GPa, Poisson ratio v = 0.167, and mass density ρ = 2400 kg/m3. The material properties of iron were Young modulus E = 200 GPa, Poisson ratio v = 0.3, and density ρ = 7800 kg/m3. In numerical model, the specimen was discretised into 10154 elements and 10370 nodes; see Figure 9. The concrete material’s initiation toughness KIcini equalled 0.888 MPa·m1/2 for slow loading case and 0.909 MPa·m1/2 for fast loading case, and its unstable fracture toughness KIcun equalled 2.914 MPa·m1/2 for slow loading case and 3.028 MPa·m1/2 for fast loading case. The parameters were tested by experiment. For the slow loading case, we set the time step Δt = 10 s, while Δt = 0.1 s for the fast loading case.Figure 8
Schematic of a notched cubic concrete specimen subjected to a splitting force (unit: mm).Figure 9
XFEM mesh.In this section, we mainly investigated the water pressure distribution along crack surface and the effect of water pressure distribution on the fracture property.


### 4.1. Splitting Force versus CMOD Curve: Comparison with Our Experiments [41]

As shown in Figures10 and 11, we plotted splitting force versus CMOD curve for slow loading case and fast loading case, respectively, with three different water pressure cases: (i) without water pressure; (ii) water pressure p=0.2 MPa at the crack mouth; (iii) water pressure p=0.4 MPa at the crack mouth. We also compared our present numerical results with our experimental results.Figure 10
Splitting force versus CMOD curve for slow loading case with different water pressures.
(a)
Without water pressure
(b)
Water pressurepn = 0.2 MPa at the crack mouth
(c)
Water pressurepn = 0.4 MPa at the crack mouthFigure 11
Splitting force vesus CMOD curve for fast loading case with different water pressures.
(a)
Without water pressure
(b)
Water pressurepn = 0.2 MPa at the crack mouth
(c)
Water pressurepn = 0.4 MPa at the crack mouthBoth numerical and experimental results show that, (i) with the increase of the applied water pressure at the crack mouth, the peak value of the splitting force decreased dramatically; (ii) under the same case, compared with the slow loading case, the maximum splitting force from the fast loading case had an obvious increase. Additionally, a fairly satisfactory agreement can still be observed for the curve in ascending section. This verification indicated that the proposed model was quite effective for simulating hydraulic fracture problems.We also listed the maximum splitting force with three different water pressure cases in Table1. The maximum error reached to 16.98%.Table 1
The maximum splitting forceFmax/kN.
Slow loading case
Fast loading case
p
n = 0.0
p
n = 0.2 MPa
p
n = 0.4 MPa
p
n = 0.0 MPa
p
n = 0.2 MPa
p
n = 0.4 MPa
Numerical results
15.60
13.31
10.75
19.51
17.02
14.00
Experimental results
18.79
14.28
10.52
19.82
16.16
12.95
Error
16.98%
6.79%
2.19%
1.56%
5.32%
8.11%


### 4.2. Effect of Water Pressure at the Crack Mouth on the Fracture Property

For experiment, the maximum water pressurepn = 0.4 MPa at the crack mouth was applied. Due to the limitation of experimental conditions, it was quite difficult in the experiment to increase the water pressure at the crack mouth for investigating its effect on the fracture property. The numerical simulation showed its huge advantages in extending test conditions. Here, we observed numerically six cases, pn = 0.0, 0.2 MPa, 0.4 MPa, 0.6 MPa, 0.8 MPa, and 1.0 MPa.Figure12(a) showed the effect of water pressure distribution on the fracture property for a slow loading case with six different water pressure cases and Figure 12(b) for fast loading case. We also plotted maximum splitting force versus applied water pressure at the crack mouth curve; see Figure 13. As depicted in Section 4.1, the results shown in these figures were still straightforward. With the increase of the applied water pressure at the crack mouth, the peak value of the splitting force decreased dramatically. Under the same case, compared with the slow loading case, the maximum splitting force from the fast loading case had an obvious increase. It was obvious that if the applied water pressure increased, the mechanical splitting force was designed to decrease to maintain the same CMOD. In Figure 13, it also could be found that if pn ≥ 0.6 MPa, the decrease velocity of the maximum splitting force slowed down.Figure 12
Effect of water pressure at the crack mouth on the fracture property.
(a)
Slow loading case
(b)
Fast loading caseFigure 13
Maximum splitting force versus applied water pressure at the crack mouth curve.The F-CMOD curve contained two phases: the prepeak/peak response and the postpeak response. The CMOD value responding to peak response was always less than 0.2 mm. With the increase of the applied water pressure at the crack mouth, the CMOD value responding to peak response decreased significantly.


### 4.3. Water Pressure Distribution along a Propagating Crack Surface

Figure14 showed the water pressure distribution along a propagating crack surface for slow and fast loading cases. To limit the length of this paper, only the case pn = 0.4 MPa was shown here. It could be seen that the water pressure distribution along crack surface followed parabolic distribution. But, at the beginning, the water pressure distribution closely approximated linear distribution. At last, the water pressure for most crack segments approximated the applied water pressure at the crack mouth. Adjacent to tip, the water pressure dropped rapidly to zero.Figure 14
Water pressure distribution along a propagating crack surface.
(a)
Slow loading case
(b)
Fast loading case


### 4.4. Failure Patterns

As shown in Figures15 and 16, we showed the failure mode for slow and fast loading cases with pn = 0.2 MPa. In numerical model, we assumed that the concrete specimen was ideal homogeneous material. Therefore, the numerical results showed that the crack was pure mode-I crack and propagated in a straight manner. The experimental results also showed that the crack propagated approximately in a straight manner, but a slight deflection could still be observed. This phenomenon could be interpreted as the heterogeneity of the actual concrete specimen.Figure 15
Failure mode for slow loading case withpn = 0.2 MPa.
(a)
x-directional displacement contours/mm
(b)
y-directional displacement contours/mm
(c)
Deformed mesh
(d)
Failure mode of experimentFigure 16
Failure mode for fast loading case withpn = 0.2 MPa.
(a)
x-directional displacement contours/mm
(b)
y-directional displacement contours/mm
(c)
Deformed mesh
(d)
Failure mode of experiment


## 4.1. Splitting Force versus CMOD Curve: Comparison with Our Experiments [41]

As shown in Figures10 and 11, we plotted splitting force versus CMOD curve for slow loading case and fast loading case, respectively, with three different water pressure cases: (i) without water pressure; (ii) water pressure p=0.2 MPa at the crack mouth; (iii) water pressure p=0.4 MPa at the crack mouth. We also compared our present numerical results with our experimental results.Figure 10
Splitting force versus CMOD curve for slow loading case with different water pressures.
(a)
Without water pressure
(b)
Water pressurepn = 0.2 MPa at the crack mouth
(c)
Water pressurepn = 0.4 MPa at the crack mouthFigure 11
Splitting force vesus CMOD curve for fast loading case with different water pressures.
(a)
Without water pressure
(b)
Water pressurepn = 0.2 MPa at the crack mouth
(c)
Water pressurepn = 0.4 MPa at the crack mouthBoth numerical and experimental results show that, (i) with the increase of the applied water pressure at the crack mouth, the peak value of the splitting force decreased dramatically; (ii) under the same case, compared with the slow loading case, the maximum splitting force from the fast loading case had an obvious increase. Additionally, a fairly satisfactory agreement can still be observed for the curve in ascending section. This verification indicated that the proposed model was quite effective for simulating hydraulic fracture problems.We also listed the maximum splitting force with three different water pressure cases in Table1. The maximum error reached to 16.98%.Table 1
The maximum splitting forceFmax/kN.
Slow loading case
Fast loading case
p
n = 0.0
p
n = 0.2 MPa
p
n = 0.4 MPa
p
n = 0.0 MPa
p
n = 0.2 MPa
p
n = 0.4 MPa
Numerical results
15.60
13.31
10.75
19.51
17.02
14.00
Experimental results
18.79
14.28
10.52
19.82
16.16
12.95
Error
16.98%
6.79%
2.19%
1.56%
5.32%
8.11%


## 4.2. Effect of Water Pressure at the Crack Mouth on the Fracture Property

For experiment, the maximum water pressurepn = 0.4 MPa at the crack mouth was applied. Due to the limitation of experimental conditions, it was quite difficult in the experiment to increase the water pressure at the crack mouth for investigating its effect on the fracture property. The numerical simulation showed its huge advantages in extending test conditions. Here, we observed numerically six cases, pn = 0.0, 0.2 MPa, 0.4 MPa, 0.6 MPa, 0.8 MPa, and 1.0 MPa.Figure12(a) showed the effect of water pressure distribution on the fracture property for a slow loading case with six different water pressure cases and Figure 12(b) for fast loading case. We also plotted maximum splitting force versus applied water pressure at the crack mouth curve; see Figure 13. As depicted in Section 4.1, the results shown in these figures were still straightforward. With the increase of the applied water pressure at the crack mouth, the peak value of the splitting force decreased dramatically. Under the same case, compared with the slow loading case, the maximum splitting force from the fast loading case had an obvious increase. It was obvious that if the applied water pressure increased, the mechanical splitting force was designed to decrease to maintain the same CMOD. In Figure 13, it also could be found that if pn ≥ 0.6 MPa, the decrease velocity of the maximum splitting force slowed down.Figure 12
Effect of water pressure at the crack mouth on the fracture property.
(a)
Slow loading case
(b)
Fast loading caseFigure 13
Maximum splitting force versus applied water pressure at the crack mouth curve.The F-CMOD curve contained two phases: the prepeak/peak response and the postpeak response. The CMOD value responding to peak response was always less than 0.2 mm. With the increase of the applied water pressure at the crack mouth, the CMOD value responding to peak response decreased significantly.


## 4.3. Water Pressure Distribution along a Propagating Crack Surface

Figure14 showed the water pressure distribution along a propagating crack surface for slow and fast loading cases. To limit the length of this paper, only the case pn = 0.4 MPa was shown here. It could be seen that the water pressure distribution along crack surface followed parabolic distribution. But, at the beginning, the water pressure distribution closely approximated linear distribution. At last, the water pressure for most crack segments approximated the applied water pressure at the crack mouth. Adjacent to tip, the water pressure dropped rapidly to zero.Figure 14
Water pressure distribution along a propagating crack surface.
(a)
Slow loading case
(b)
Fast loading case


## 4.4. Failure Patterns

As shown in Figures15 and 16, we showed the failure mode for slow and fast loading cases with pn = 0.2 MPa. In numerical model, we assumed that the concrete specimen was ideal homogeneous material. Therefore, the numerical results showed that the crack was pure mode-I crack and propagated in a straight manner. The experimental results also showed that the crack propagated approximately in a straight manner, but a slight deflection could still be observed. This phenomenon could be interpreted as the heterogeneity of the actual concrete specimen.Figure 15
Failure mode for slow loading case withpn = 0.2 MPa.
(a)
x-directional displacement contours/mm
(b)
y-directional displacement contours/mm
(c)
Deformed mesh
(d)
Failure mode of experimentFigure 16
Failure mode for fast loading case withpn = 0.2 MPa.
(a)
x-directional displacement contours/mm
(b)
y-directional displacement contours/mm
(c)
Deformed mesh
(d)
Failure mode of experiment


## 5. Summary and Conclusions

In this paper, we model the fluid flow within the crack as one-dimensional flow and assume that the flow is laminar; the fluid is incompressible and accounts for the time-dependent rate of crack opening. Here, we discretise the flow equation by finite volume methods. The extended finite element methods are used for solving solid medium with crack under dynamic loads. Having constructed the approximation of dynamic extended finite element methods, the derivation of governing equation for dynamic extended finite element methods is presented. The implicit time algorithm is elaborated for the time descritisation of dominant equation. In addition, the interaction integral method is given for evaluating stress intensity factors. Then, the coupling model for modelling hydraulic fracture can be established by the extended finite element methods and the finite volume methods. We compare our present numerical results with our experimental results for verifying the proposed model. Finally, we investigate the water pressure distribution along crack surface and the effect of water pressure distribution on the fracture property. Some valuable conclusions can be drawn from this study.(i)
Some conclusions between numerical results and experimental results were quite identical. A fairly satisfactory agreement could also be observed for F-CMOD curve in ascending section. Therefore, the proposed model was quite effective for simulating hydraulic fracture problems.(ii)
The F-CMOD curve contained two phases: the prepeak/peak response and the postpeak response. The CMOD value responding to peak response was always less than 0.2 mm. With the increase of the applied water pressure at the crack mouth, the CMOD value responding to peak response decreased significantly.(iii)
With the increase of the applied water pressure at the crack mouth, the peak value of the splitting force decreased dramatically. Under the same case, compared with the slow loading case, the maximum splitting force from the fast loading case had an obvious increase.(iv)
The water pressure distribution along crack surface followed parabolic distribution. But, at the beginning, the water pressure distribution closely approximated linear distribution. At last, the water pressure for most crack segments approximated the applied water pressure at the crack mouth. Adjacent to tip, the water pressure dropped rapidly to zero.(v)
In numerical model, we assumed that the concrete specimen was ideal homogeneous material. Therefore, the numerical results showed that the crack was pure mode-I crack and propagated in a straight manner. The experimental results also showed that the crack propagated approximately in a straight manner, but a slight deflection could still be observed. This phenomenon could be interpreted as the heterogeneity of the actual concrete specimen.Our future work will further investigate the water pressure distribution along a propagating crack surface and the effect of water pressure distribution on the fracture property with considering crack opening-closing.



---

*Source: 3751340-2017-11-28.xml*