Abstract
A numerical simulation approach of hydraulic fracture process, considering the couplings of the stress distribution, the fluid flow of the water-air mixture, the compression and dissolution of air, and the element damage evolution, has been developed to investigate the mechanisms of crack initiation and propagation in porous media during hydraulic fracturing. The concept of homogenized pore fluid has been adopted to represent the water air mixture. A large number of numerical analysis on hydraulic fracturing in clay with incipient injection slot have been carried out to study the mechanism of hydraulic fracturing in unsaturated soil with the characteristic of critical model I type of crack loading using stress intensity factor . The results provide a numerical picture depicting the mechanisms of crack initiation and propagation during hydraulic fracturing. The numerical results are in good agreement with the experimental results, which confirms the adequacy and the power of the numerical approach.
1. Introduction
Hydraulic fracturing may be defined as the process of creating a fracture or fracture system in a porous medium by injecting a fluid under pressure through a well bore in order to overcome native stresses. Hydraulic fracturing has been used for more than 50 years to enhance the yield of wells recovering oil at great depth. In the recent past, hydraulic fracturing methods have been developed for creating fractures in soils and making existing fractures larger to enhance the mass transfer of contaminants. The fractures created increase the effective permeability and change paths of fluid flow, thus making in situ remediation more effective and economical. Hydraulic fracturing appears to have useful environmental, agricultural and geotechnical applications.
Hydraulic fracturing techniques are equally applicable to both saturated soils and unsaturated soils to improve the flow of water and air, respectively. Earth’s land is mostly located in the unsaturated vadose zone above the water table, so the hydraulic fracturing of unsaturated soil has important significance. An unsaturated soil is commonly referred to as a 3-phase system composed of solids (soil particles), water, and air. The air-water interface (i.e., the contractile skin) warrants inclusion as an additional phase due to its unique and specific properties. The contractile skin interacts with the soil particles in an independent manner and can significantly change the mechanical behavior of an unsaturated soil. The hydraulic fracturing in unsaturated soil is coupling process of stress distribution, water flow, air flow, water-air interact, and damage evolution, which is arguably one of the most challenging problems.
Especially, the propagation processes of hydraulic fractures enveloped in soil layers are difficult to directly observe, and the details of fracture growth inside soil have generally been inferred from indirect measurements. During the past few years, with the rapid development of computing power, numerical tools have become a good option for gaining some insight into the fracturing process. The main numerical methods coupled with fracture models have been used to simulate fracturing process by many researchers, such as the finite element method (FEM) [1], the mesh-free Galerkin method (EFGM) [2], the extended finite element method (XFEM) [3], the boundary element method (BEM) [4], the discrete element method (DEM) [5], and the displacement discontinuity method (DDM) [6]. The former four are based on the continuum mechanics, the latter two are based on the noncontinuum mechanics.
Hydraulic fracturing is a more complicated process to simulate, as it involves the coupling of four processes [7, 8]: (i) the mechanical deformation induced by the fluid pressure on the fracture surfaces; (ii) the flow of fluid within the fracture; (iii) the fracture propagation; and (iv) the coupling of seepage and stress in the computational domain around the fracture. The simulation of hydraulic fracturing in rock and soil mass, arguably one of the most challenging computational problems in geoengineering, has been the subject of numerous investigations since the pioneering work of Khristianovic and Zheltov [9]. Many numerical methods have been developed for the simulation of hydraulic fracturing [10–14]. Among these methods, FEM method is the most robust one in comparison with BEM, DEM, and DDM which cannot efficiently solve the elasticity-plasticity equation relating the fluid pressure to the fracture opening [15].
Two approaches, namely, the discrete crack model and the smeared crack model were introduced to FEM by Ngo and Scordelis [16] and Rashid [17], respectively. The discrete crack model is aimed at simulating the initiation and propagation of dominant cracks. In contrast, the smeared crack model is based on the idea that many small cracks nucleate which only in a later stage of the loading process link up to form one or more dominant cracks.
The discrete crack approach in its original form has several disadvantages [18]. In the approach, the cracks are forced to propagate along the element boundaries, so that a mesh bias is introduced. Automatic remeshing with sophisticated computer codes allows the mesh bias to be reduced, nevertheless, a computational difficulty, namely, the continuous change in topology, is inherent in the discrete crack approach and is to a certain extent even aggravated by remeshing procedures.
Since each individual crack is not numerically resolved, the smeared crack model captures the deterioration process through a constitutive relation, thus smearing out the cracks over the continuum. It describes a cracked solid by an equivalent anisotropic continuum with degraded material properties in the direction that is normal to the crack orientation and no remeshing is needed. Although the discrete crack approach is sometimes preferred for being based on a mature theory, the smeared crack analysis seems to be more often applicable for usual FEM analyses from engineering practice.
However, many assumptions and simplifications are made in existing numerical studies; for example, fissures need to be placed inside the analyzed domain prior to analysis [19, 20] or the hardening of grout in soil cannot be described [21]. A numerical model for coupled analysis of flow, stress, and damage was proposed to overcome the shortcoming of the placing of preset of fractures [22]. The approach is similar to smeared crack method. The simulated domain is discretized into large numbers of tiny elements to take into account the small local variations of the material properties. At each loading step, the stress and strain in the elements are calculated and examined against the predefined soil strength. Those elements with the stresses above the material’s strength were considered to be isotropically damaged. The material properties of the damaged element are reduced and the width of the facture in the damaged element can be calculated. This approach can effectively simulate hydraulic fracturing [22–24] in rocks. In this paper, the extension of the approach for rocks to unsaturated soil is attempted so that hydraulic fracturing in unsaturated soil can be numerically simulated.
In this paper, the Biot’s consolidation theory is used as a framework that governs the interactions between the soil skeleton and the pore fluid. The concept of homogenized pore fluid is used to describe water-air mixture fluid. The nonlinear constitutive relationship developed by Duncun and Chang is used to describe the soil deformation for its simplicity. The finite element method (FEM) [25] is employed to discretize the Biot equation. The FEM is chosen because it can readily obtain the solutions of boundary value problems on nonconformal domains.
2. The Concept of an Homogenized Pore Fluid
The effective stress equation for saturated soil has been shown to be theoretically valid and practically useful for saturated soils.
For unsaturated soil, Bishop (1959) [26] suggested the following rational equation: in which = effective stress, = total stress, = pore air pressure, = pore water pressure, and = factor whose value ranges from 0 to 1 as degree of saturation varies, and is evaluated experimentally.
However, Sparks [27] has shown that can assume values greater than unity, and that it can also be negative. Coleman [28], Blight [29], and Matyas and Radharkrishna [30] have shown that the value of applicable to volume change behaviour is different from the value of applicable to the strength behaviour of the same soil at the same degree of saturation. This fact indicates that (2) is not fundamentally valid and is not likely to be applicable except under limited ranges of conditions and for limited purposes. Sparks [27] employed a model composed of uniform spheres to develop the following effective stress equation for unsaturated soils: in which = parameters, whose values depend on the degree of saturation, and = surface tension of water.
Both of the preceding effective stress equations involve difficulties for practical application. The concept of a “homogenized pore fluid” presented by [31] offers a simple alternative to these equations which is applicable for a limited but useful range of conditions.
The pores of an unsaturated soil are filled partly with gas (air) and partly with liquid (water), and the air phase might be present in an unsaturated soil either in a continuous or an occluded (bubbles) form. An important and frequently encountered special case is that in which the degree of saturation is sufficiently high so that the air bubbles are occluded. This condition prevails at degrees of saturation greater than about 0.85, which is common in many practical cases.
The concept of a “homogenized pore fluid” is that this mixture of gas and liquid can be considered to be an equivalent homogeneous pore fluid which completely fills the pores of the soil. As the proportions of air and water vary (due to flow of water and air, or due to compression of air and dissolution of air), the mechanical properties of the homogenized pore fluid vary. The compressibility of the pore fluid and changes in the degree of saturation due to changes in pore fluid pressure may be calculated using Boyle’s Law and Henry’s Law [32]. The equation of effective stress may be expressed as in which = effective average pore pressure in the homogenized pore fluid.
Sparks [27] and Barden [33] have found that when the air bubbles are occluded, the pore air pressures and the pore water pressures are very nearly equal. There will inevitably be some small differences between these two pressures as a result of surface tension effects, but these small differences become negligible as the magnitude of the pore pressure increases. Sparks has indicated that for occluded bubbles, the assumption that introduces very small error. Thus, for this condition, the effective stress equation may be written as follows: in which = pore water pressure. Equation (5) has the same form with (1), but the fluid in (5) is water-air mixture fluid which can be compressed, while (1)’s is incompressible water. Equation (5) is used throughout the following developments, and the results are thus applicable to soils containing water and occluded air bubbles.
3. The Compressibility and Permeability of Pore Fluid
The compressibility () of pore fluid containing air and water is defined as in which = volumetric strain in pore fluid due to change () in the fluid pressure. By employing Boyle’s Law and Henry’s Law, the compressibility of an air-water mixture can be expressed in a form suggested by [34]
The first term in the equation accounts for the compressibility of the water; the second term accounts for Boyle’s Law being applied to the free air; and the third term accounts for the air driven into solution in accordance with Henry’s Law. The pore pressure parameter is equal to 1 as the air bubbles become occluded, where = compressibility of the air-water mixture, = volumetric coefficient of air solubility (at 20°C, ), = compressibility of the water, which can be ignored because it is a slight amount, and = degree of saturation.
With the increase in pore pressure pore air dissolved in the water phase, the saturation variation can be calculated as [35] where = initial degree of saturation and = atmospheric pressure.
The permeability of unsaturated clays () depends on the degree of saturation. The relationship between the permeability of saturated soil () and its degree of saturation can often be approximated using the following relationship [36, 37]: where is the threshold degree of saturation at which the water in the pores begins to flow freely. For a degree of saturation less than , the water forms a membrane around the particles and does not flow.
4. Field Equation Governing the Consolidation Phenomenon
Chang and Duncan [31] present an extension of Biot’s theory which makes possible finite element analyses of the consolidation of inelastic, unsaturated soil masses. The basic equations governing consolidation may be stated as follows.
(a) Equilibrium equation. Only changes in stresses and body forces are considered. The geostatic body forces are considered as initial stresses. Consider the following: where is the total stress tensor, is the body force per unit of mass, and is the density of the mass.
(b) Total stresses are the sum of the effective stress and the pore pressure . Consider the following:
(c) Continuity equation. For a compressible pore fluid, the volume decrease of the soil is the sum of the water expelled plus the decrease in volume of the fluid within. The continuity equation can therefore be written as where is the superficial velocity vector, is the displacement vector, and is the compressibility of the pore fluid. “·” indicates partial differentiation with respect to time.
(d) The effective stress-strain relationship may be expressed in a general form as follows: where is the strain tensor, which is related to the displacements through where is the tensor relating stress and strain.
(e) Darcy’s Law. It is assumed that the flow of water is governed by Darcy’s Law as follows: where is the permeability tensor, is the body force per unit of fluid, and is the density of the fluid.
To define the problem, both displacement and flow boundary conditions must be specified. For the displacement boundary conditions, it is assumed that part of the boundary surface is subjected to known applied tractions, , while the remainder of the surface is subjected to specified displacements, . For the flow boundary condition, it is assumed that part of the boundary surface is subjected to a specified flow velocity, , while the remainder of the surface is subjected to known pore pressures, . The boundary conditions may be written in the following form: where is the normal vector to the boundary surface and (−) indicates a prescribed quantity.
5. Constitutive Equation
The nonlinear elastic model proposed by Duncan and Chang is used as the constitutive relationship for solids. The main features of the model are summarized here. The tangential Young’s modulus can be expressed as
The elastic Young’s modulus for unloading and reloading is
The tangential volumetric modulus is described by where is cohesion, is internal friction angle, is atmospheric pressure, , , are initial modulus values, and are parameters, and is a failure ratio that is defined as where is the principal stress difference when the soil is assumed to be in failure and is the asymptotic value of the principal stress difference.
Eight parameters of the constitutive model , , , , , , , and can be determined by conventional triaxial tests.
6. Mechanics of Fracture Initiation and Propagation
The tensile and shear failures are both considered here. An element is considered to be fractured in the tension mode if its minor effective principal stress exceeds the tensile strength of the element as follows: where is the minor effective principal stress and is the tensile failure strength of the element. Compression is assumed to be positive.
An element is considered to be fractured in the shear mode if the shear stress satisfied the Mohr-Coulomb failure criterion as follows: where is the shear stress, the effective normal stress, the effective stress cohesion intercept, and the effective stress angle of friction or shearing resistance.
Once an element is judged to be fractured, that is, mesh element damaged, many smeared cracks occur in this element and the stiffness of the soil is reduced.
The stiffness of the fracture can be described by [38] where is the initial modulus of the soil just prior to the occurrence of cracks, is the total width of smeared cracks in an element (m), and is a scaling parameter indicating the amount of deformation needed to change the stiffness by a factor of . , also follow the formula. The total width of smeared cracks is the difference between the widths of the element after and before cracking. The value of is related to the material types and is 10−4 in this investigation.
The permeability is not a constant and is a function of stresses because the fracture aperture is most likely to change as the stress conditions vary. The effect of stresses on permeability can be expressed as [22] where is the damage factor of permeability reflecting the increase in permeability induced by damage of soil, and are coupling parameters reflecting the influence of stress on the coefficient of permeability.
7. The Simulation of of Unsaturated Soil during Hydraulic Fracture
The critical stress at the onset of propagation depended on the length of the initial slot. The strength was characterized by the stress intensity factor [39, 40] as follows: where is far-field stress, is the half-length of the fracture, and is some function of the geometry of the fracture and enveloping material.
A basic assumption is that propagation under dilation, or mode I, loading occurs when the stress intensity equals a critical value . The magnitude of is commonly regarded as a material property characterizing the resistance of a material to mode I fracturing; here it is called fracture toughness.
Fracture toughness has been applied to prediction of the onset of hydraulic fracturing of soil [41–43], but the validity of that application needs to be examined further.
The stress intensity factor of soil can be defined as follows: Driving pressure is the difference between the fluid pressure within a fracture and the confining stress acting normal to the fracture as follows:
The function in (26) is equal to unity if the fracture is straight, loaded uniformly, and embedded in a homogeneous, isotropic, linearly-elastic material, and it can be taken as a correction factor for practical constraints such as finite sample size or slots intersected by a hole. Singh [37] have given expressions for for a straight slot containing an axial hole and for a straight slot embedded in a rectangular sample of finite size. These conditions represent corrections of <5% for the geometries and dimensions used in this work, so will be set to unity. So, fracture toughness will be determined using the critical driving pressure and the half-length of the initial slot as
The propagation of hydraulic fractures in soil was studied in the laboratory by glycerine injected at a constant rate into rectangular samples of silty clay confined in a triaxial cell [41–43]. In order to examine the validity of the application, a large number of numerical analyses on hydraulic fracturing in clay with incipient injection slot have been carried out to investigate the relation between and . And the results are compared with experimental results [41–43].
A simulation domain 0.1 m 0.39 m with incipient injection slot is the same with the test [41–43]. The soil mass within the simulated area is under the plane strain boundary condition, as shown in Figure 1.
The relationship between initial modulus values and saturation can be obtained from the experimental results [41] as follows:
The input parameters for these simulations are tabulated in Table 1.
Six cases of the initial saturation from 0.85 to 0.95 at interval 0.02, 12 cases of length of the starter slot from 0.008 m to 0.03 m at interval 0.002 m, 2 cases of consolidation durations of 4 days and 14 days are considered. A flow rate of 0.033 cm3/s which was the same with the test was used for all the simulations described below.
The numerical analysis can capture details of the at various saturation in the range >85%. The result is shown in Figures 2 and 3. The laboratory records [42] and the simulations in Figures 2 and 3 show a remarkable similarity. The simulations can accommodate changes caused by the saturation and the length of start slot.
(a) Consolidated for 14 days
(b) Consolidated for 4 days
(a)
(b)
Both saturation and duration of consolidation had a marked effect on . Among samples consolidated for 14 days, for m, is greatest (22 ) at . The fracture toughness decreases abruptly to 2.33 as increases from 0.85 to 0.95. For m, the fracture toughness decreases abruptly from 12.4 to 1.68 . The sharp change in occurs in the range . A slight decrease in occurs with further increase in , although the decrease is small and is almost independent of in the range 0.91–0.95.
Consolidation appears to toughen the clay, according to simulation results of samples consolidated for 4 days and 14 days under a static pressure of 69 kPa. The of the samples consolidated for 14 days is greater by 6.1 than the of 4 days consolidated samples at similar and m. With the increase of saturation, the difference of is reduced, which means that the tough effect of consolidation is weakened.
The forms of the experimental and simulation curves in Figures 2 and 4 are similar; they all show decreasing markedly over a few percent of .
| (a) m |
| (b) m |
| (c) m |
Fracture toughness appears to be independent of the length of the starter slot in the range and to be dependent of the length of the starter slot in the range under the conditions used here. Increase of the length of the starter slot can decrease this dependence.
As we all know, negative pore pressure can affect the strength of the soil. The area of negative pore pressure zone (NPZ) may explain the dependence between fracture toughness and saturation and length of starter slot. Figure 4 shows the NPZ before facture tip. Only half of the simulated domain is displayed because of symmetry. When the saturation is higher (0.95), the area of NPZ and the length of the starter slot are almost independent of each other (Figure 4). When the saturation is lower (0.85), the length of the starter slot has a greater impact on the area of NPZ (Figure 4).
Fracture toughness has been introduced as a parameter that could possibly be used to indicate the onset of propagation. If it is a material property, fracture toughness must be independent of slot length. It can depend on saturation, however, although the dependence is expected to be similar to other properties that also depend on saturation.
From the simulation results, the depends on the length of the initial slot and the saturation of the soil. It is apparent from Figures 2 and 3 that the value of varies over more than a factor of ten, so that some methods of anticipating will be required in order for this to be a useful predictor of the onset of propagation.
8. The Simulation of Hydraulic Fracture in Heterogeneity Unsaturated Soil
Actual soil is often heterogeneous, in order to demonstrate the effectiveness of this approach to simulate hydraulic fracture in heterogeneous unsaturated soil, three simulations are carried out. A simulation domain 0.1 m 0.1 m with initial injection tiny holes in the center is considered here. All boundaries are fixed displacement, undrained boundary. The distribution of heterogeneity remains the same and the saturation varies from one sample to another.
To take the heterogeneity into account, different material properties for different elements are assigned randomly following the Weibull distribution as follows: where is the variable representing a certain material property, is the mean value of the corresponding material property, and is the homogeneity index. is a parameter defined by the shape of the distribution function that describes the degree of material heterogeneity. A larger implies a more homogeneous material and vice versa. Therefore, the parameter is called the homogeneity index. The higher the value of the homogeneity index , the greater the number of elements with properties close to . In the study, , , , , are accorded to the Weibull distribution by .
Figure 5 shows the simulated fracture propagation process. Each column is a simulation result at a saturation. Saturation has a significant impact on fracture propagation and distribution. When saturation is high (), numerous cracks gathered near the initial hole and cross-cuttings occur between the cracks. The lower the saturation, the farther away the cracks spead from the initial hole and the fewer the number of cracks there is.
(a)
(b)
(c)
(d)
Hydraulic fracturing appears to have useful environmental and geotechnical applications, but the details of fracture morphology and propagation in soil are poorly known because the propagation processes of hydraulic fractures enveloped in soil layers are difficult to directly observe. The simulation results may provide some insight.
9. Conclusions
A numerical simulation approach of hydraulic fracture process of unsaturated soil with high saturation has been presented. The approach can consider the coupling of stress distribution, water-air mixture fluid flow, compression and dissolution of air, and element damage evolution. A key to the approach is the concept of homogenized pore fluid which is that the mixture of gas and liquid can be considered to be an equivalent homogeneous pore fluid which completely fills the pores of the soil. As the proportions of air and water vary (due to flow of water and air, or due to compression of air and dissolution of air), the mechanical properties of the homogenized pore fluid vary correspondingly.
A large number of numerical analyses on hydraulic fracturing in unsaturated soil with incipient injection slot have been carried out to study the characteristics of critical stress intensity factor of model I type of crack loading. The numerical results were shown to be in good agreement with the experimental results. The propagation of hydraulic fracturing in unsaturated soil has been simulated also. The simulation results have indicated that saturation is a fundamental constraint on the hydraulic fracturing in clay.
Although this study is a preliminary attempt and needs further development, the simulation results may provide some insight of the crack initiation and propagation mechanisms during hydrofracturing of unsaturated soil with high saturation.
Conflict of Interests
The authors declare that there is no conflict of interests regarding the publication of this paper.
Acknowledgments
The authors gratefully acknowledge the financial support of the National Basic Research Program of China under Grant 2010CB732100 and the National Natural Science Foundation of China under Grant 51278041.