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

Wemodel the fluid flowwithin the crack as one-dimensional flowand assume that the flow is laminar; the fluid is incompressible and accounts for the time-dependent rate of crack opening.Here,we discretise the flowequation byfinite volumemethods.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.


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][2][3][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][10][11][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][19][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.Section 2 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., is partitioned into three parts: the displacement boundary (Γ u ), the traction boundary (Γ t ), and the crack boundary (Γ c ) that is tractionfree.The elastodynamic basic equation is expressed as

Governing Equations
with the following boundary and initial conditions: where  is the Cauchy stress tensor, b is the body force vector,  is the strain tensor,  is the material density, ü 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.

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] ∇ ⋅ q + ẇ +  = 0, where ∇ is the divergence operator defined in x direction; q is the fluid flux; ẇ is the time-dependent rate of crack opening, and ẇ = /; and  is the fluid loss into the solid media, and here we ignore the fluid loss; that is,  = 0. Additionally, Poiseuille's law [29] gives the following expression: where  is the crack opening;  is the fluid viscosity; and  is the fluid pressure.
The pressure boundary conditions at the fluid injection point in the crack are where  0 is the pressure of the fluid injection point.
According to the fluid pressure continuity, lag condition (6) provides the net-pressure boundary condition at the fluid front  =  t for the fluid flow equations (3).

Numerical Approximation
where   (x) is the standard finite element shape function of node ; u  is the unknown of the standard finite element part at node ;  is the set of all nodes in the domain; and  *  (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; a  and b   are the nodal enriched degree of freedom;  * abs and  * br are the set of enrichment nodes shown in Figure 1, and  * abs ,  * br ⊂ .For these elements that are cut completely by a crack, the nodes of these elements that are the nodal subset  * abs are enriched by Heaviside function (x).The definition of Heaviside function (x) follows: 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 subset where  and  are the local crack tip coordinate system.

Discrete Equations.
By the principle of virtual work, the following discrete equations can be obtained: where K (M) is the global stiffness (mass) matrix assembled by the element stiffness (mass) matrix; f is the global external force vector; u h and ü 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 The element stiffness matrix is expressed by where The element mass matrix is expressed by where The element external force vector is where (17)

Time Integration Schemes.
The following time-integration scheme is used in dynamic analysis.Equation ( 10) for a specific time  + Δ is expressed as where u  , u  , and ü  are the displacement, velocity, and acceleration vectors at time , respectively; Δ is the time step;  and  are parameters that can be determined to obtain integration accuracy and stability, with 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 matrix K, and mass matrix M.
(ii) Give the initial displacement vector u 0 and the initial velocity vector u 0 .Then, calculate the initial acceleration vector ü 0 by the equilibrium equation: (iii) Select a time step Δ and the parameters  and .
Here,  = 0.275625 and  = 0.55 are used.Calculate integration constants: (iv) Form the effective stiffness matrix K: (II) For Each Time Step (i) Calculate effective loads f+Δ at time  + Δ: (ii) Solve for the displacement vector u +Δ at time  + Δ: (iii) Calculate the acceleration vector ü +Δ and the velocity vector u +Δ at time  + Δ: 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.
To 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.

Interaction Integral for Computing Stress Intensity Factors.
Take field 1, ( (1)   ,  (1)   ,  (1)   ), for the actual field, and field 2, ( (2)   ,  (2)   ,  (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: (1)    (2)  ,1 +  (2)    (1)  ,1 − 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 ; Γ d + ∪ Γ d − consists of a interface starting from the external integration radius to crack tip in a two-way manner. is defined as where ℎ e is the crack-tip element size;  k is a user-specified scalar multiple;  is the weight function;  = 1 if the node lies in A; and  = 0 if the node lies outside of A or lies on the boundary of A. The weight function  in the interior of an element is obtained by the interpolation of the nodal value: Additionally, the interaction integral relates to the stress intensity factors through the relation: where  aux I and  aux II are the local auxiliary stress intensity factors for the auxiliary fields, respectively; and the definition of  * is By setting  aux I = 1 and  aux II = 0 as well as  (1,2) =  (1,2)   1 , we obtain the expression of  I as follows: Similarly, we obtain the equality by setting  aux I = 0 and  aux II = 1 and  (1,2) =  (1,2) 2 .
3.1.6.Crack Propagation Criteria.The maximum circumferential stress criterion [37] is used to determine the crack growth direction.Once  I and  II are calculated, the criterion gives the following crack growth direction: where  c is the crack growth angle in the local crack-tip coordinate system.If  II = 0, then  c = 0.It should also be noted that if  II > 0, the crack growth angle  c < 0, and if  II < 0, then  c > 0. Sukumar and Prévost (2003) [38] proposed a computationally more amenable expression for  c : The equivalent stress intensity factor then follows: The double-K criterion [39] is used for determining crack propagation.The crack propagation processes can be expressed as follows:

𝐾 ini
Ic is the initiation toughness;  un Ic is the unstable fracture toughness.3) and ( 4) and ignoring the fluid loss into the solid media, the following expression can be obtained:

Finite Volume Methods of the Flow along a Crack. Considering (
where where  is the cross-sectional area of the control volume face; Δ is the volume;  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 To calculate gradients at the control volume faces, a linear approximate distribution of pressure is considered here.So (39) yields w ,  w ,  e , and  e can be obtained by the way of linearly interpolated values, and Equation ( 40) can be rewritten as where (43)

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 Equation ( 44) can be rewritten as where w and  w can be calculated directly according to crack opening displacement, and  w is the pressure of the fluid injection point.
(II) Boundary Conditions for Fluid Front.Similarly, as shown in Figure 6, integration of (36) over the control volume surrounding node n gives Equation ( 47) can be rewritten as where e and  e can be calculated directly according to crack width at the fluid front, and  t is the pressure of the fluid front and  t = 0.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.

Coupling Scheme.
To correctly solve the system of equations, the elasticity, fluid flow, and fracture growth should be coupled together; see Figure 7.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.
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  are stated as follows: (i) Assume that the fluid pressure   has been computed at the kth 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   at the th step iteration.(ii) According to the known crack width   at the kth step iteration, compute the fluid pressure  +1/2 by using FVM.(iii) Compute the fluid pressure  +1 at the ( + 1)th step iteration by the flowing expression: (iv) Compute dynamic responses of solid medium under dynamic loads and fluid pressure  +1 by using XFEM.According to dynamic responses, we can obtain the crack width  +1/2 .
(v) Compute the crack width  +1 at the ( + 1)th step iteration by the flowing expression: (vi) Given tolerance  = 0.01, compute the following expression: If  ≤ , finish iteration, or else go to (ii).

Numerical Example
As shown in Figure 8, 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  = 36 GPa, Poisson ratio V = 0.167, and mass density  = 2400 kg/m 3 .The material properties of iron were Young modulus  = 200 GPa, Poisson ratio V = 0.3, and density  = 7800 kg/m 3 .In numerical model, the specimen was discretised into 10154 elements and 10370 nodes; see Figure 9.The concrete material's initiation toughness  ini Ic equalled 0.888 MPa⋅m 1/2 for slow loading case and 0.909 MPa⋅m 1/2 for fast loading case, and its unstable fracture toughness  un Ic equalled 2.914 MPa⋅m 1/2 for slow loading case and 3.028 MPa⋅m 1/2 for fast loading case.The parameters were tested by experiment.For the slow loading case, we set the time step Δ = 10 s, while Δ = 0.1 s for the fast loading case.In this section, we mainly investigated the water pressure distribution along crack surface and the effect of water pressure distribution on the fracture property.

Splitting Force versus CMOD Curve: Comparison with
Our Experiments [41].As shown in Figures 10 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  = 0.2 MPa at the crack mouth; (iii) water pressure  = 0.4 MPa at the crack mouth.We also compared our present numerical results with our experimental results.
Both 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 Table 1.The maximum error reached to 16.98%.

Effect of Water Pressure at the Crack Mouth on the Fracture
Property.For experiment, the maximum water pressure   = 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,   = 0.0, 0.2 MPa, 0.4 MPa, 0.6 MPa, 0.8 MPa, and 1.0 MPa. Figure 12(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   ≥ 0.6 MPa, the decrease velocity of the maximum splitting force slowed down.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.

Water Pressure Distribution along a Propagating Crack
Surface. Figure 14 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   = 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.15 and 16, we showed the failure mode for slow and fast loading cases with   = 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.

Summary and Conclusions
In this paper, we model the fluid flow within the crack as onedimensional 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.

Figure 1 :
Figure 1: A strategy for enriched elements and nodes.

Figure 2 :
Figure 2: Element partitioning method for these elements containing a discontinuous interface.

Figure 3 :
Figure 3: Elements selection for the interaction integral near the crack tip.

Figure 7 :
Figure 7: The coupling of the elasticity, fluid flow, and fracture growth.

Figure 8 :
Figure 8: Schematic of a notched cubic concrete specimen subjected to a splitting force (unit: mm).

Figure 10 :
Figure 10: Splitting force versus CMOD curve for slow loading case with different water pressures.

Figure 11 :
Figure 11: Splitting force vesus CMOD curve for fast loading case with different water pressures.

2 Figure 12 :Figure 13 :
Figure 12: Effect of water pressure at the crack mouth on the fracture property.

Table 1 :
The maximum splitting force  max /kN.