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.
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. 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.
2. Governing Equations
2.1. Elastic Dynamic Responses of the Solid Medium
The boundary of a bounded domain, , is partitioned into three parts: the displacement boundary , the traction boundary , and the crack boundary that is traction-free. The elastodynamic basic equation is expressed aswith the following boundary and initial conditions:where is the Cauchy stress tensor, is the body force vector, is the strain tensor, is the material density, is the acceleration field vector, is the symmetric part of the gradient operator, is the displacement field vector, is the constitutive matrix, is the unit outward normal vector to the crack surface, is the prescribed displacement, is the external traction vector, is the initial displacement vector, and 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]where is the divergence operator defined in x direction; 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, .
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 arewhere is the pressure of the fluid injection point.
In the fluid lag zone [30], According to the fluid pressure continuity, lag condition (6) provides the net-pressure boundary condition at the fluid front 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 aswhere is the standard finite element shape function of node ; is the unknown of the standard finite element part at node ; is the set of all nodes in the domain; and 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; and are the nodal enriched degree of freedom; and are the set of enrichment nodes shown in Figure 1, and .
For these elements that are cut completely by a crack, the nodes of these elements that are the nodal subset are enriched by Heaviside function . The definition of Heaviside function follows:where is the projection of a point on the crack surface; 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 are enriched by the crack tip enrichment function . The definition of the crack tip enrichment function follows:where 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:where is the global stiffness (mass) matrix assembled by the element stiffness (mass) matrix; is the global external force vector; and denote the vector of nodal parameters (which include the classic degrees of freedom, , and the enrichment degrees of freedom, , ) and its second time derivative, respectively; and
The element stiffness matrix is expressed bywhere
The element mass matrix is expressed bywhere
The element external force vector iswhere
3.1.3. Time Integration Schemes
The following time-integration scheme is used in dynamic analysis. Equation (10) for a specific time is expressed aswhere , , 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 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 , and mass matrix .(ii)Give the initial displacement vector and the initial velocity vector . Then, calculate the initial acceleration vector by the equilibrium equation:(iii)Select a time step and the parameters and . Here, and are used. Calculate integration constants:(iv)Form the effective stiffness matrix :
(II) For Each Time Step(i)Calculate effective loads at time : (ii)Solve for the displacement vector at time : (iii)Calculate the acceleration vector and the velocity vector 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.
(a) Elements partitioned completely by a straight crack
(b) Elements partitioned completely by a broken crack
(c) Elements partitioned partially by a straight crack
To solve the element stiffness or mass matrix of these enrichment elements, each subquad element is, respectively, transferred into the standard element 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, , for the actual field, and field 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:
The second term in (26) denotes the contribution of traction along crack interface. As shown in Figure 3, denotes the circle domain with centre at the crack tip and the radius ; consists of a interface starting from the external integration radius to crack tip in a two-way manner. is defined aswhere is the crack-tip element size; is a user-specified scalar multiple; is the weight function; if the node lies in A; and if the node lies outside of or lies on the boundary of . 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 and are the local auxiliary stress intensity factors for the auxiliary fields, respectively; and the definition of is
By setting and as well as , we obtain the expression of as follows:
Similarly, we obtain the equalityby setting and and .
3.1.6. Crack Propagation Criteria
The maximum circumferential stress criterion [37] is used to determine the crack growth direction. Once and are calculated, the criterion gives the following crack growth direction:where is the crack growth angle in the local crack-tip coordinate system. If , then . It should also be noted that if , the crack growth angle , and if , then . Sukumar and Prévost (2003) [38] proposed a computationally more amenable expression for :
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:(i)If , the crack does not propagate.(ii)When , the performed crack begins to crack initially.(iii)When , the propagating crack develops steadily.(iv)When , the crack propagates unsteadily.
is the initiation toughness; 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:where
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 and . 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.”
Integrating the governing equation in (36) over a control volume, the following discretised equation can be obtained.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, , , and can be obtained by the way of linearly interpolated values, and
Equation (40) can be rewritten aswhere
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 givesEquation (44) can be rewritten aswhere and can be calculated directly according to crack opening displacement, and 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 givesEquation (47) can be rewritten aswhere and can be calculated directly according to crack width at the fluid front, and is the pressure of the fluid front and .
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 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 by using FVM.(iii)Compute the fluid pressure at the th step iteration by the flowing expression:(iv)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 .(v)Compute the crack width at the th step iteration by the flowing expression:(vi)Given tolerance , compute the following expression:If , finish iteration, or else go to (ii).
4. 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 = 0.167, and mass density = 2400 kg/m3. The material properties of iron were Young modulus = 200 GPa, Poisson ratio = 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 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 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 = 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.
4.1. 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 MPa at the crack mouth; (iii) water pressure MPa at the crack mouth. We also compared our present numerical results with our experimental results.
(a) Without water pressure
(b) Water pressure = 0.2 MPa at the crack mouth
(c) Water pressure = 0.4 MPa at the crack mouth
(a) Without water pressure
(b) Water pressure = 0.2 MPa at the crack mouth
(c) Water pressure = 0.4 MPa at the crack mouth
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%.
4.2. 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.
(a) Slow loading case
(b) Fast loading case
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
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.
(a) Slow loading case
(b) Fast loading case
4.4. Failure Patterns
As shown in Figures 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.
(a) x-directional displacement contours/mm
(b) y-directional displacement contours/mm
(c) Deformed mesh
(d) Failure mode of experiment
(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.
Conflicts of Interest
The authors declare that there are no conflicts of interest regarding the publication of this paper.
Acknowledgments
The authors gratefully acknowledge support for this research from the National Natural Science Foundation of China (Grants nos. 51309088, 11372098, 51579084, and 51308188), the Fundamental Research Funds for the Central Universities (Grant no. 2015B01714), the National Sci-Tech Support Plan of China (Grant no. 2015BAB07B10), and China Postdoctoral Science Foundation funded project (Grant no. 2014T70466). The authors also would like to thank Professor Charles Augarde for the discussion during their visit to Durham University, who stayed in Hohai University supported by the Royal Society.