Numerical Simulation of Fluid-Solid Coupling in Fractured Porous Media with Discrete Fracture Model and Extended Finite Element Method

Fluid-solid coupling is ubiquitous in the process of fluid flow underground and has a significant influence on the development of oil and gas reservoirs. To investigate these phenomena, the coupled mathematical model of solid deformation and fluid flow in fractured porous media is established. In this study, the discrete fracture model (DFM) is applied to capture fluid flow in the fractured porous media, which represents fractures explicitly and avoids calculating shape factor for cross flow. In addition, the extended finite element method (XFEM) is applied to capture solid deformation due to the discontinuity caused by fractures. More importantly, this model captures the change of fractures aperture during the simulation, and then adjusts fluid flow in the fractures. The final linear equation set is derived and solved for a 2D plane strain problem. Results show that the combination of discrete fracture model and extended finite element method is suited for simulating coupled deformation and fluid flow in fractured porous media.


Introduction
The technology of hydraulic fracturing has been widely used for reservoir stimulation, especially for unconventional reservoirs [1].Coupled rock deformation and fluid flow in fractured porous media is important for reservoir simulation because rock deformation exerts an important influence on reservoir production [2].
The general theory of 3D consolidation with elasticity constitutive relationship and Darcy' law has been established by Biot [3], and the effective stress formulation has been put forward by Terzaghi [4].A great number of researches about fluid-solid coupling have been done based on these theories in petroleum engineering, from conventional reservoirs to fractured reservoirs [5].
There exist several methods to simulate fluid flow in fractured porous media [6].Warren and Root introduced the dual continuum concept to characterize naturally fractured reservoirs [7].The dual continuum approaches treat fracture and matrix both as continua distributed within reservoir domain.The fracture-matrix cross flow is based on the analytical solution of pseudosteady-state flow within the matrix system with a simple geometry of matrix blocks.Moreover, shape factors are calculated for different geometries of matrix blocks.The dual continuum approaches consist of dual porosity and dual permeability models [8].Schematics for dual porosity concept and dual permeability concept are shown in Figures 1 and 2 respectively.The difference between dual porosity model and dual permeability model is that dual permeability model takes global matrix-matrix flow into account while dual porosity model does not account for it.An alternative to the dual continuum approaches is the discrete fracture model [9,10].The fractures are represented explicitly within the domain and discretized along with the matrix domain.Lamb [11] presented fracture mapping approach (FM) to simulate fluid flow in fractured porous media.In the FM approach, an element intersected by a fracture is treated as a superposition of two elements, namely a matrix element and a fracture element.The matrix element and fracture element interact via a transfer function.The schematic of fracture mapping approach is shown in Figure 3.The approach adopts the transfer function presented by Barenblatt [12] to account for cross flow between the overlapping matrix and fracture elements, which was a dual continuum model to this extent.Due to stress singularity of fracture tip, it needs mesh refinement around the fracture in the standard finite element framework, which is computational burdensome.An alternative to standard finite element method is the extended finite element method (XFEM).The extended finite element method was introduced by Belytschko [13] to discontinuous problems, which has been widely used in many fields due to flexibility in meshing [14][15][16][17].The fracture is represented by level set method.
Ghafouri and Lewis [18,19] presented a finite element continuum approach to describe the coupled deformation and fluid flow in fractured porous media, which was based on double porosity model.Tran [20] presented high level boundary element method with periodic boundary conditions and flux continuous finite volume element method to simulate coupled fluid flow through discrete fracture network; Al-Khoury [21] used the partition of unity method to describe the fracturing process and double porosity model to describe the resulting fluid flow; Vire et al. [22] presented coupling an adaptive mesh finite element fluid model with a combined finite-discrete element solid model to investigate fluid-solid interactions.The method is flexible in terms of discretization schemes used for each material.Recently, Vire et al. [23] presents an immersed-shell method for modeling fluid-structure interactions.The method consists of immersing the solid structures in an extended fluid domain, and exchanging the coupling forces through a thin shell surrounding the solid structures.Lamb [13] presented FM approach and the extended finite element method for coupled deformation and fluid flow in fractured porous media.The difference between this paper and Lamb's is that here discrete fracture model is used to avoid calculating cross flow and the model captures the change of fractures aperture during the simulation.
The advantages of combination of discrete fracture model and extended finite element method over other methods are that discrete fracture model avoids the computation of shape factor and is more accurate than dual continuum model for simulating fluid flow with large fractures, meanwhile the extended finite element method avoids mesh refinement around the fracture and is well suited for m Ω f Ω discontinuity problems.Furthermore, the model is capable of capturing change of fractures aperture during the simulation.
In this paper, the combination of discrete fracture model and extended finite element method is used to couple deformation and fluid flow in fractured porous media.The governing equations and initial and boundary conditions are presented in Section 2. The numerical solution is presented in Section 3. In the section, the extended finite element method and the discrete fracture model are briefly described to capture deformation and fluid flow respectively, and spatial and temporal discretization are conducted.Finally, numerical simulation and result analysis are performed in Section 4 and the conclusion are drawn in Section 5.

Governing Equations for Rock Deformation
Under the assumptions of small-strain situation, isothermal equilibrium and negligible inertial forces, Biot's theory describes the linear momentum balance equation for a two-phase medium, which is composed by rock and water.T 0 where σ is total stress tensor, g is gravity, ρ is the averaged density of the multiphase system, and ∇ is the symmetric gradient operator matrix.
where φ is porosity, s ρ is the density of solid, and w ρ is the density of water.
The relationship between total stress and effective stress is given by e where e σ is the effective stress, w p is the water pore pressure in the porous matrix, α is the Biot's compressibility coefficient, and The constitutive stress-strain relationship of the solid phase is given by where D is the linear elastic material matrix and ε is the strain of the system.The geometric equation between strain and displacement is given by where u is the displacement of the system.

Governing Equations for Fluid Flow
Fluid flow in fractured porous media is typically simulated using dual-porosity models, but dual-porosity models are not well suited for the modeling of a small number of large-scale fractures, which may dominate the flow.Discrete fracture model, in which the fractures are represented individually, has been broadly applied to simulate flow in fractured porous media.In this paper, discrete fracture model is used for flow simulation.
The equation of motion for fluid flow in porous media is given by Darcy' law as follows ( ) where v is velocity, k is permeability, w μ is viscosity of water.
The continuity equation takes into consideration the grain and fluid volume variation resulting from pressure change (the first term) and total volume strain resulting from solid deformation (the second term), which is given by T where s K is the bulk modulus of the grain material, w K is the bulk modulus of water.

Initial and Boundary Conditions
The initial conditions specify the displacement and water pressure fields at time 0 t = .0 0 , in and on where Ω is the domain of interest and Γ is the boundary.Boundary conditions of solid deformation include displacement condition and force condition, which can be given as follows where l is related to the unit outward normal vector T { , } Boundary conditions of fluid flow include water pressure condition and flux condition, which can be given as follows T ( ) on where w q is the imposed volume flux normal to the boundary.

Application of the Extended Finite Element Method
In the finite element framework, discontinuity modeling needs mesh refinement around the crack, which is computationally burdensome.The extended finite element method was introduced by Belytschko and Black, which had been widely used for discontinuous problem.The advantage of extended finite element method is avoiding mesh refinement around fractures (especially fracture tips) without reducing accuracy.The method is adopted to solve solid deformation in this study.
The extended finite element method exploits the partition of unity property of finite element.The displacement field is decomposed into two parts: the continuous displacement field and the discontinuous part.The displacement approximation can be written as follows where From Equation (13), the opening between the two surfaces of the fracture can be given by 1 ( ) 2 2 where n is the unit outward normal vector.

Application of Discrete Fracture Model
In discrete fracture model, fractures are simplified into 1D line element for 2D problem, and 2D surface area element for 3D problem, as shown in Figure 4.The whole fractured porous media is decomposed into two parts: matrix system and fracture system, which can be given as follows where m represents matrix, f represents fracture, and a is the aperture of fracture.
Assuming that representative element volumes of both matrix and fracture system exist, flow equations are applicable to the whole research region.Then for the discrete fracture model, the integral form of the flow equation can be given as follows.
In the Equations ( 16) and ( 17), the apertures of fractures are calculated from Equation (15).In the implementation of algorithm, fractures are discretized into small segments.Moreover, the aperture of each segment is set equal to average opening of its endpoints, which is easy to calculate from Equation (15).However, the initial aperture of fracture is given in the first time step and updated from the second time step.

Discretization in Space
According to strong form of solid deformation equation, the weak form can be expressed with weighted residual method T 0 ( ) : : where is the trial function space, and cr Γ is the fracture in the domain.
The weighted residual method is applied to the continuity equation for water flow and to its natural boundary condition, which yields where w is the weight function, and Ω is the region of integration expressed by Equation ( 16), which consists of matrix and fracture system.Equation ( 19) turns to where m k is the permeability of matrix, and f k is the permeability of fracture, which can be given by "cubic law" as follows The expression for displacement is given by Equation ( 13), and the expression for water pressure is given as follows w w p p = N p (22) where w p is the vector of the nodal values of water pressure, and p N is the shape function for water pressure.By integrating Equations ( 18) and ( 20), the weak form of the whole system is discretized into the following set of equations where Elements of the aforementioned listed matrices are given in Appendix.

Discretization in Time
Using the fully implicit time discretization scheme, the approximation is given as follows Then, the final discrete equation can be written as follows The scheme is fully implicit and imposes no requirements of the time step size which is usually chosen for both stability and the elimination oscillatory effects in the solution.
In the above equation, the unknown vector X includes standard degrees of freedom of all nodes, enriched degrees of freedom of enriched nodes in two directions and water pressures of all nodes.The matrices are constructed according to Appendix A and the linear equation set is solved for 1 n+ X (value for n + 1 time step) providing that n X (value for n time step) is given.Moreover, the initial value for X is given with the initial conditions by setting nodal displacements and enriched degrees to zero and water pressures to 0 w p .The fracture apertures are calculated according to Equation (15) using the already solved displacements, then the matrix for flow in the fracture is updated.The loop continues until end time is reached.

Numerical Example
The DFM-XFEM model is applied to a 2D plain strain problem shown in Figure 5.The domain is fully saturated and allowed to freely drain at the top, that is, excess pore water pressure is equal to zero.The domain has an area of 10 m × 16 m.The fracture is 8 m long, is inclined at 45° and is centered in the middle of the domain.The right, left and bottom boundaries are assumed to be undrained.The lateral boundaries of the domain are constrained to vertical translation only, and the bottom boundary is constrained to be fixed.Static load is applied at the top boundary and maintained throughout the duration of the simulation.The input parameters are given in Table 1.
The discrete fracture model has been widely used for simulating flow in fractured porous media, and the extended finite element method has been broadly applied to the discontinuity problem.The combination DFM-XFEM model gives full play to their advantages for simulating flow solid coupling in fractured porous media.The solid module shares the same mesh configuration with the fluid flow module.The numerical method is implemented using Matlab.
Lamb presented fracture mapping approach and the extended finite element method to couple deformation and fluid flow in fractured porous media, denoted by FM-XFEM model.Moreover, the model proposed in this paper is denoted by DFM-XFEM model.Because the fracture geometry remains constant throughout the simulation for the FM-XFEM model, that is, the fracture does not close or open during the simulation, the assumption is applied to the DFM-XFEM model for model verification by comparison with FM-XFEM model.The assumption is done by fixing fracture aperture in the fluid flow module during the simulation.The comparison results of the two models are shown in Figures 6 and 7.   discontinuous due to the pre-existing fracture, and the extended finite element method captures discontinuity without mesh refinement around the fracture, which is a great improvement over the standard finite element method.
The excess pore water pressure distributions are shown in Figure 7.The results of two different models are very close, and they are also qualitatively identical.The discrete fracture model represents fracture explicitly, and it does not need to calculate shape factor, which is not easy to determine in the dual porosity model.Because the top boundary is a zero pressure boundary, fluid can freely flow out and the pressure is becoming lower and lower.
The displacement fields along y direction at point A (located at the top shown in Figure 5 and excess pore pressure at point B (located at the bottom shown in Figure 5 for the FM-XFEM model and DFM-XFEM model at varying mesh resolutions during 100 days are shown in Figures 8a and 9a respectively.The plots are zoomed in for 10 days to clarify the discrepancies between the various curves and shown in Figure 8b and Figure 9b.These two figures indicate that the results of DFM-XFEM model show great coincidence with the result of FM-XFEM, which validates the correctness of DFM-XFEM model.Moreover, as the mesh resolutions become higher, the results of DFM-XFEM model show little difference, which validates the convergence of the proposed algorithm.The computational cost of models is shown in Table 2.It can be shown that the computational time of DFM-XFEM is much less than that of FM-XFEM with approximately equal number of elements.The comparison shows the advantage of DFM-XFEM over FM-XFEM.In the above example, the fracture remains constant during the simulation.However, the fracture could close or open during the simulation, so it needs to capture change of the fracture aperture in the fluid flow.The DFM-XFEM model is able to simulate the situation when fracture aperture changes as described in Section 3.2.The model parameters are the same as the above example.The comparison results for a model with changed and fixed fracture aperture are shown in Figures 10 and 11. Figure 10 illustrates the displacement at point A for 100 days and indicates that the displacement for model with fixed fracture aperture is larger than that for model with changed fracture aperture after a short time and the displacements for these two models show little difference after 50 days.Figure 11 illustrates the excess pore pressure at point B for 100 days and indicates that excess pore pressure for a model with fixed fracture aperture is larger than that for model with changed fracture aperture after a short time and excess pore pressures for these two models show little difference after 50 days.These results are seen because the fracture aperture is becoming smaller, as shown in Figure 12. Figure 12 shows that the fracture aperture decreases quickly after a short time and slows down the decreasing rate after a longer time.As the fracture aperture decreases, the fracture conductivity decreases, which results in pore pressure decreasing more slowly.The method is applied to another case, which consists of two intersecting fractures.The domain has an area of 10 m × 10 m.The results of FM-XFEM model and DFM-XFEM are also compared, as shown in Figures 13 and 14.It can be shown from Figure 13 that the displacement fields along y direction are very close for two methods.The y displacements in the area above the fractures are much larger than that below the fractures.It can be shown from Figure 14 that the excess pore pressure distributions of two methods are identical.The pressure around the fracture is smaller than other area because of larger conductivities in the fractures.The water pressures at point (located at the right bottom) are shown in Figure 15a at varying mesh resolutions during 100 days.The plot is zoomed in for 10 days to clarify the discrepancies between the various curves as shown in Figure 15b.Also, the results of FM-XFEM model are shown in the figures.It can be shown that as the mesh resolutions become higher, the results converge to that of FM-XFEM.The differences between them are very small.
With regard to computational cost, the results are shown in Table 3.It can be shown that the computational cost of DFM-XFEM model is much less than that of FM-XFEM model.The results prove that the proposed method is better than FM-XFEM model.

Conclusions
In this paper, the combination of discrete fracture model and extended finite element method is proposed to simulate fluid-solid coupling in fractured porous media.The discrete fracture model is able to capture fluid flow accurately without accessing cross flow between matrix and fracture.The extended finite element method is capable of solving solid deformation without mesh refinement around fractures tips.The model captures change of fracture aperture during the simulation.The results of numerical example show that the proposed method is well suited for simulating fluid-solid coupling in fractured porous media.

Figure 1 .
Figure 1.Schematic illustration of dual porosity model of fractured reservoirs.(a) Actual reservoir; (b) Reservoir model.

Figure 2 .
Figure 2. Schematic illustration of dual permeability model of fractured reservoirs.

Figure 4 .
Figure 4. Schematic of a discrete fracture model.

Figure 5 .
Figure 5. Two-dimensional fractured domain with assigned boundary conditions.

Figure 10 .
Figure 10.Displacement at point A for 100 days.

Figure 11 .
Figure 11.Excess pore pressure at point B for 100 days.

Figure 12 .
Figure 12.Fracture aperture distribution at different time.

Table 1 .
Input parameters for the models.

Table 2 .
Computational cost of models.

Table 3 .
Computational cost of models.