Coupled formulation and algorithms for the simulation of non‐planar three‐dimensional hydraulic fractures using the generalized finite element method

This paper presents a coupled hydro‐mechanical formulation for the simulation of non‐planar three‐dimensional hydraulic fractures. Deformation in the rock is modeled using linear elasticity, and the lubrication theory is adopted for the fluid flow in the fracture. The governing equations of the fluid flow and elasticity and the subsequent discretization are fully coupled. A Generalized/eXtended Finite Element Method (G/XFEM) is adopted for the discretization of the coupled system of equations. A Newton–Raphson method is used to solve the resulting system of nonlinear equations. A discretization strategy for the fluid flow problem on non‐planar three‐dimensional surfaces and a computationally efficient strategy for handling time integration combined with mesh adaptivity are also presented. Several three‐dimensional numerical verification examples are solved. The examples illustrate the generality and accuracy of the proposed coupled formulation and discretization strategies. Copyright © 2015 John Wiley & Sons, Ltd.     

An embedded fracture modeling framework for simulation of hydraulic fracturing and shear stimulation

A numerical modeling framework is described that is able to calculate the coupled processes of fluid flow, geomechanics, and rock failure for application to general engineering problems related to reservoir stimulation, including hydraulic fracturing and shear stimulation. The numerical formulation employs the use of an embedded fracture modeling approach, which provides several advantages over more traditional methods in terms of computational complexity and efficiency. Specifically, the embedded fracture modeling strategy avoids the usual requirement that the discretization of the fracture domain conforms to the discretization of the rock volume surrounding the fractures. As fluid is exchanged between the two domains, conservation of mass is guaranteed through a coupling term that appears as a simple source term in the governing mass balance equations. In this manner, as new tensile fractures nucleate and propagate subject to mechanical effects, numerical complexities associated with the introduction of new fracture control volumes are largely negated. In addition, the ability to discretize the fractures and surrounding rock volume independently provides the freedom to choose an acceptable level of discretization for each domain separately. Three numerical examples were performed to demonstrate the utility of the embedded fracture model for application to problems involving fluid flow, mechanical deformation, and rock failure. The results of the numerical examples confirm that the embedded fracture model was able to capture accurately the complex and nonlinear evolution of reservoir permeability as new fractures propagate through the reservoir and as fractures fail in shear.     

A finite element solution of a fully coupled implicit formulation for reservoir simulation

An iterative method is presented for solving a fully coupled and implicit formulation of fluid flow in a porous medium. The mathematical model describes a set of fully coupled three-phase flow of compressible and immiscible fluids in a saturated oil reservoir. The finite element method is applied to obtain the simultaneous solution (SS) for the resulting highly non-linear partial differential equations where fluid pressures are the primary unknowns. The final discretized equations are solved iteratively by using a fully implicit numerical scheme. Several examples, illustrating the use of the present model, are described. The increased stability achieved with this scheme has permitted the use of larger time steps with smaller material balance errors.     

Coupled hydromechanical‐fracture simulations of nonplanar three‐dimensional hydraulic fracture propagation

This paper presents an algorithm and a fully coupled hydromechanical‐fracture formulation for the simulation of three‐dimensional nonplanar hydraulic fracture propagation. The propagation algorithm automatically estimates the magnitude of time steps such that a regularized form of Irwin's criterion is satisfied along the predicted 3‐D fracture front at every fracture propagation step. A generalized finite element method is used for the discretization of elasticity equations governing the deformation of the rock, and a finite element method is adopted for the solution of the fluid flow equation on the basis of Poiseuille's cubic law. Adaptive mesh refinement is used for discretization error control, leading to significantly fewer degrees of freedom than available nonadaptive methods. An efficient computational scheme to handle nonlinear time‐dependent problems with adaptive mesh refinement is presented. Explicit fracture surface representations are used to avoid mapping of 3‐D solutions between generalized finite element method meshes. Examples demonstrating the accuracy, robustness, and computational efficiency of the proposed formulation, regularized Irwin's criterion, and propagation algorithm are presented.     

Numerical analysis of thermal effects during carbon dioxide injection with enhanced gas recovery: a theoretical case study for the Altmark gas field

Prediction about reservoir temperature change during carbon dioxide injection requires consideration of all, often subtle, thermal effects. In particular, Joule?CThomson cooling (JTC) and the viscous heat dissipation (VHD) effect are factors that cause flowing fluid temperature to differ from the static formation temperature. In this work, warm-back behavior (thermal recovery after injection completed), as well as JTC and VHD effects, at a multi-layered depleted gas reservoir are demonstrated numerically. OpenGeoSys (OGS) is able to solve coupled partial differential equations for pressure, temperature and mole-fraction of each component of the mixture with a combination of monolithic and staggered approaches. The Galerkin finite element approach is adapted for space discretization of governing equations, whereas for temporal discretization, a generalized implicit single-step scheme is used. For numerical modeling of warm-back behavior, we chose a simplified test case of carbon dioxide injection. This test case is numerically solved by using OGS and FeFlow simulators independently. OGS differs from FeFlow in the capability of representing multi-componential effects on warm-back behavior. We verify both code results by showing the close comparison of shut-in temperature profiles along the injection well. As the JTC cooling rate is inversely proportional to the volumetric heat capacity of the solid matrix, the injection layers are cooled faster as compared to the non-injection layers. The shut-in temperature profiles are showing a significant change in reservoir temperature; hence it is important to account for thermal effects in injection monitoring.     

Hydro‐mechanical modeling of cohesive crack propagation in multiphase porous media using the extended finite element method

In this paper, a numerical model is developed for the fully coupled hydro‐mechanical analysis of deformable, progressively fracturing porous media interacting with the flow of two immiscible, compressible wetting and non‐wetting pore fluids, in which the coupling between various processes is taken into account. The governing equations involving the coupled solid skeleton deformation and two‐phase fluid flow in partially saturated porous media including cohesive cracks are derived within the framework of the generalized Biot theory. The fluid flow within the crack is simulated using the Darcy law in which the permeability variation with porosity because of the cracking of the solid skeleton is accounted. The cohesive crack model is integrated into the numerical modeling by means of which the nonlinear fracture processes occurring along the fracture process zone are simulated. The solid phase displacement, the wetting phase pressure and the capillary pressure are taken as the primary variables of the three‐phase formulation. The other variables are incorporated into the model via the experimentally determined functions, which specify the relationship between the hydraulic properties of the fracturing porous medium, that is saturation, permeability and capillary pressure. The spatial discretization is implemented by employing the extended finite element method, and the time domain discretization is performed using the generalized Newmark scheme to derive the final system of fully coupled nonlinear equations of the hydro‐mechanical problem. It is illustrated that by allowing for the interaction between various processes, that is the solid skeleton deformation, the wetting and the non‐wetting pore fluid flow and the cohesive crack propagation, the effect of the presence of the geomechanical discontinuity can be completely captured. Copyright © 2012 John Wiley & Sons, Ltd.     

Transient drainage volume characterization and flow simulation in reservoir models using the fast marching method

The Fast Marching Method (FMM) has been applied to characterize the transient drainage volume and to simulate flow as a function of time in porous media using the concept of the “diffusive time of flight” (DTOF). The DTOF (τ) provides a spatial coordinate which reduces the three dimensional pressure diffusivity equation to an equivalent one dimensional formulation. It is obtained from the solution to the Eikonal equation via the FMM. Previous applications of this approach have solved the flow equations numerically or by using an analytic asymptotic approximation. Both solution approaches rely upon three characteristics. (1) Accurate solution for the DTOF irrespective of the degree of heterogeneity within a reservoir model. (2) The approximation of the three dimensional pressure solution in terms of the τ coordinate. (3) Accurate representation and discretization of the drainage volume, on which the asymptotic and numerical flow simulations are based. The second and third of these characteristics are specific to reservoir engineering applications, and provide the focus of this study. Analysis of the drainage volume shows that the near well region requires special treatment, leading to a composite discretization for the drainage volume. This discretization has a direct impact upon the calculation of the well test pressure derivative when the asymptotic approximation is used for pressure transient interpretation. For flow simulation, the discretization directly impacts the calculation of the well index and the intercell transmissibility computed in the τ coordinate, and places additional constraints on the discretization of the drainage volume. These new results are validated by comparison with a commercial finite difference flow simulator, and are shown to be more accurate than earlier computational approaches.

     

Heterogeneity preserving upscaling for heat transport in fractured geothermal reservoirs

In simulation of fluid injection in fractured geothermal reservoirs, the characteristics of the physical processes are severely affected by the local occurence of connected fractures. To resolve these structurally dominated processes, there is a need to develop discretization strategies that also limit computational effort. In this paper, we present an upscaling methodology for geothermal heat transport with fractures represented explicitly in the computational grid. The heat transport is modeled by an advection-conduction equation for the temperature, and solved on a highly irregular coarse grid that preserves the fracture heterogeneity. The upscaling is based on different strategies for the advective term and the conductive term. The coarse scale advective term is constructed from sums of fine scale fluxes, whereas the coarse scale conductive term is constructed based on numerically computed basis functions. The method naturally incorporates the coupling between solution variables in the matrix and in the fractures, respectively, via the discretization. In this way, explicit transfer terms that couple fracture and matrix solution variables are avoided. Numerical results show that the upscaling methodology performs well, in particular for large upscaling ratios, and that it is applicable also to highly complex fracture networks.     

Discrete fracture model for coupled flow and geomechanics

We present a fully implicit formulation of coupled flow and geomechanics for fractured three-dimensional subsurface formations. The Reservoir Characterization Model (RCM) consists of a computational grid, in which the fractures are represented explicitly. The Discrete Fracture Model (DFM) has been widely used to model the flow and transport in natural geological porous formations. Here, we extend the DFM approach to model deformation. The flow equations are discretized using a finite-volume method, and the poroelasticity equations are discretized using a Galerkin finite-element approximation. The two discretizations—flow and mechanics—share the same three-dimensional unstructured grid. The mechanical behavior of the fractures is modeled as a contact problem between two computational planes. The set of fully coupled nonlinear equations is solved implicitly. The implementation is validated for two problems with analytical solutions. The methodology is then applied to a shale-gas production scenario where a synthetic reservoir with 100 natural fractures is produced using a hydraulically fractured horizontal well.     

Embedded discontinuity approach for coupled hydromechanical analysis of fractured porous media

In this paper, a new continuum approach for the coupled hydromechanical analysis of fractured porous media is proposed. The methodology for describing the hydraulic characteristics invokes an enriched form of Darcy's law formulated in the presence of an embedded discontinuity. The constitutive relations governing the hydromechanical response are derived by averaging the fluid pressure gradient and the discontinuous displacement fields over a selected referential volume of the material, subject to some physical constraints. The framework incorporates an internal length scale which is explicitly embedded in the definition of gradient operators. The respective field equations are derived following the general form of balance equations in interacting continua. The conventional finite element method is then employed for the spatial discretization, and the generalized Newmark scheme is used for the temporal discretization. The proposed methodology is verified by some numerical examples dealing with a steady-state flow through fractured media as well as a time-dependent consolidation in the presence of a discontinuity.     

Thermo-hydro-mechanical modeling of fracturing porous media with two-phase fluid flow using X-FEM technique

In this paper, a fully coupled thermo-hydro-mechanical model is presented for two-phase fluid flow and heat transfer in fractured/fracturing porous media using the extended finite element method. In the fractured porous medium, the traction, heat, and mass transfer between the fracture space and the surrounding media are coupled. The wetting and nonwetting fluid phases are water and gas, which are assumed to be immiscible, and no phase-change is considered. The system of coupled equations consists of the linear momentum balance of solid phase, wetting and nonwetting fluid continuities, and thermal energy conservation. The main variables used to solve the system of equations are solid phase displacement, wetting fluid pressure, capillary pressure, and temperature. The fracture is assumed to impose the strong discontinuity in the displacement field and weak discontinuities in the fluid pressure, capillary pressure, and temperature fields. The mode I fracture propagation is employed using a cohesive fracture model. Finally, several numerical examples are solved to illustrate the capability of the proposed computational algorithm. It is shown that the effect of thermal expansion on the effective stress can influence the rate of fracture propagation and the injection pressure in hydraulic fracturing process. Moreover, the effect of thermal loading is investigated properly on fracture opening and fluids flow in unsaturated porous media, and the convective heat transfer within the fracture is captured successfully. It is shown how the proposed computational model is capable of modeling the fully coupled thermal fracture propagation in unsaturated porous media.     

储层流固耦合的数学模型和非线性有限元方程

根据饱和多孔介质固体骨架的平衡方程和多孔介质中流体的连续性方程,建立了储层流固耦合数学模型。模型中引入了Jaumann应力速率公式描述多孔介质固体骨架的大变形效应,并考虑了地应力、初始孔隙压力、初始流体密度和初始孔隙度对耦合模型的影响。基于与微分方程等价的加权余量公式,在空间域采用有限元离散,对时间域进行隐式差分格式离散,导出了以单元节点位移和单元节点孔隙压力为未知量的储层流固耦合的非线性有限元增量方程。该模型在石油工程中有广泛的应用,为储层流固耦合的数值模拟奠定了理论基础。     

Implementation of a Locally Conservative Numerical Subgrid Upscaling Scheme for Two-Phase Darcy Flow

We present a locally mass conservative scheme for the approximation of two-phase flow in a porous medium that allows us to obtain detailed fine scale solutions on relatively coarse meshes. The permeability is assumed to be resolvable on a fine numerical grid, but limits on computational power require that computations be performed on a coarse grid. We define a two-scale mixed finite element space and resulting method, and describe in detail the solution algorithm. It involves a coarse scale operator coupled to a subgrid scale operator localized in space to each coarse grid element. An influence function (numerical Greens function) technique allows us to solve these subgrid scale problems independently of the coarse grid approximation. The coarse grid problem is modified to take into account the subgrid scale solution and solved as a large linear system of equations posed over a coarse grid. Finally, the coarse scale solution is corrected on the subgrid scale, providing a fine grid representation of the solution. Numerical examples are presented, which show that near-well behavior and even extremely heterogeneous permeability barriers and streaks are upscaled well by the technique.     

Large‐scale poroelastic fractured reservoirs modeling using the fast multipole displacement discontinuity method

An effective approach to modeling the geomechanical behavior of the network and its permeability variation is to use a poroelastic displacement discontinuity method (DDM). However, the approach becomes rather computationally intensive for an extensive system of cracks, particularly when considering coupled diffusion/deformation processes. This is because of additional unknowns and the need for time‐marching schemes for the numerical integration. The Fast Multipole Method (FMM) is a technique that can accelerate the solution of large fracture problems with linear complexity with the number of unknowns both in memory and CPU time. Previous works combining DDM and FMM for large‐scale problems have accounted only for elastic rocks, neglecting the fluid leak‐off from the fractures into the matrix and its influence on pore pressure and stress field. In this work we develop an efficient geomechanical model for large‐scale natural fracture networks in poroelastic reservoirs with fracture flow in response to injection and production operations. Accuracy and computational performance of the proposed method with those of conventional poroelastic DDM are compared through several case studies involving up to several tens of thousands of boundary elements. The results show the effectiveness of the FMM approach to successfully evaluate field‐scale problems for the design of exploitation strategies in unconventional geothermal and petroleum reservoirs. An example considering faults reveals the impact of reservoir compartmentalization because of sealing faults for both geomechanical and flow variables under elastic and poroelastic rocks. Copyright © 2015 John Wiley & Sons, Ltd.     

Phase-state identification bypass method for three-phase thermal compositional simulation

In this paper, we propose a strategy to bypass the phase identification of fluid mixtures that can form three, or more, phases. The strategy is used for reservoir simulation of multicomponent, three-phase, thermal compositional displacement processes. Since the solution path in compositional space is determined by a limited number of “key” tie-simplexes, the proposed “bypass” method uses information from the parameterized tie-simplexes and their extensions. The tie-simplex parameterization is performed in the discrete phase-fraction space. Once the phase-fraction space is discretized, a conventional three-phase flash is used adaptively to compute the phase states at the discretization nodes. If all discretization vertices of a given discrete cell, in phase-fraction space, have the same phase state, then this state is assigned to the entire cell and expensive flash calculations are bypassed. We demonstrate the robustness and efficiency of our phase identification bypassing strategy for several cases with three-phase flow, including a six-component ES-SAGD (enhanced solvent SAGD) model.     

A general method for simulating reactive dissolution in carbonate rocks with arbitrary geometry

The two-scale continuum model is widely used in simulating the reactive dissolution process and predicting the optimum injection rate for carbonate reservoir acidizing treatment. The numerical methods of this model are currently based on structured grids, which are not applicable for complicated geometries. In this study, a general numerical scheme for simulating a reactive flow problem on both structured and unstructured grids is presented based on the finite volume method (FVM). The convection and diffusion terms involved in the reactive flow model are discretized by using the upwind scheme and two-point flux approximation (TPFA), respectively. The location of the centroid node inside each control volume is moved by using an optimization algorithm to make the connections with the surrounding elements as orthogonal as possible, which systematically improves the accuracy of the TPFA scheme. Additionally, in order to avoid the computational complexity resulting from the discretization of the non-linear term, the mass balance equation is only discretized in the spatial domain to get a set of ordinary differential equations (ODEs). These ODEs are coupled with the reaction equations and then solved using the numerical algorithm on ODEs. The accuracy and efficiency of the proposed method are studied by comparing the results obtained from the proposed numerical method with previous experimental and numerical results. This comparison indicates that, compared with the previous methods, the proposed method predicts the wormhole structure more accurately. Finally, the presented method is used to check the effect of the domain geometry, and it is found that the geometry of the flow domain has no effect on the optimum injection velocity, but the radial domain requires a larger breakthrough volume than the linear domain when other parameters are fixed.     

Numerical fixed domain mapping solution of free-surface flows coupled with an evolving interior field

The evolution of a gravity-driven free-surface flow of varying horizontal extent which couples with a field evolving within the flow is solved using a finite difference discretization of a mapping of the problem onto the unit square. Since the size of the solution domain may show several orders of magnitude of variation, while the normalized geometry of the domain and the internal field may not vary significantly, this procedure avoids excessively fine or coarse discretizations, as well as interpolations at the boundary. The parabolic and hyperbolic evolution equations for the internal field are considered. The evolution of the coupled system is solved by an implicit marching scheme. The discretizations in space and in time are accurate to second order. Multipoint upwinding is used to avoid an instability arising advective terms are large. The evolution equations are nonlinear, and are solved using a nested Newton–Raphson procedure. The nesting is achieved by using successively better approximations to the ture evolution equations. The matrix equation that arises is solved by a conjugate-gradient-like (ORTHOMIN) iteration procedure with an incomplete Cholesky factorization preconditioning. The method has a wide variety of potential applications in the earth sciences, with the ability to describe glacier flow, lava flow, avalanching and landslides. Some calculations of the thermomechanical evolution of ice-sheets are given as illustrations, and the possible existence of thermally induced instabilities is considered.     

A scalable multistage linear solver for reservoir models with multisegment wells

Wellbore flow and interactions between wells and the reservoir can be complex. Accurate modeling of these behaviors is especially important for multilateral and other advanced wells. This paper describes a new scalable linear solver for flow simulation of detailed reservoir models with advanced wells and well groups. A general purpose research simulator serves as the computational platform, in which a multisegment well (MsWell) model is used to describe wellbore flow. In the MsWell model, the wellbore is discretized into a number of segments. Hence, the MsWell model adds a large number of equations and unknowns, which are fully coupled to the reservoir model. Operating constraints on groups of wells add one more level of complexity to the system. The new linear solver is a generalized two-stage constrained pressure residual preconditioner. A global pressure system is obtained algebraically in the first stage. The system represents the pressure coupling between the reservoir and wells accurately. The well groups are disaggregated into individual multisegment wells, which then are further reduced to a standard well-like form. The two-stage scheme serves as the inner loop of a generalized minimum residual solver. Algebraic multigrid is used to compute the first-stage pressure solution; a special block-based incomplete lower–upper preconditioner is used for the second stage. We demonstrate the superior performance of this new solver compared with state-of-the-art methods using a variety of highly detailed reservoir models with complex wells and well groups.     

A novel finite element double porosity model for multiphase flow through deformable fractured porous media

Based on the theory of double-porosity, a novel mathematical model for multiphase fluid flow in a deforming fractured reservoir is developed. The present formulation, consisting of both the equilibrium and continuity equations, accounts for the significant influence of coupling between fluid flow and solid deformation, usually ignored in the reservoir simulation literature. A Galerkin-based finite element method is applied to discretize the governing equations both in the space and time domain. Throughout the derived set of equations the solid displacements as well as the fluid pressure values are considered as the primary unknowns and may be used to determine other reservoir parameters such as stresses, saturations, etc. The final set of equations represents a highly non-linear system as the elements of the coefficient matrices are updated during each iteration in terms of the independent variables. The model is employed to solve a field scale example where the results are compared to those of ten other uncoupled models. The results illustrate a significantly different behaviour for the case of a reservoir where the impact of coupling is also considered. © 1997 by John Wiley & Sons, Ltd.     

A semi-implicit method for incompressible three-phase flow in porous media

In this paper, we present a semi-implicit method for the incompressible three-phase flow equations in two dimensions. In particular, a high-order discontinuous Galerkin spatial discretization is coupled with a backward Euler discretization in time. We consider a pressure-saturation formulation, decouple the pressure and saturation equations, and solve them sequentially while still keeping each equation implicit in its respective unknown. We present several numerical examples on both homogeneous and heterogeneous media, with varying permeability and porosity. Our results demonstrate the robustness of the scheme. In particular, no slope limiters are required and a relatively large time step may be taken.