Transient Pressure Behavior of Reservoirs with Discrete Conductive Faults and Fractures

Fractures and faults are common features of many well-known reservoirs. They create traps, serve as conduits to oil and gas migration, and can behave as barriers or baffles to fluid flow. Naturally fractured reservoirs consist of fractures in igneous, metamorphic, sedimentary rocks (matrix), and formations. In most sedimentary formations both fractures and matrix contribute to flow and storage, but in igneous and metamorphic rocks only fractures contribute to flow and storage, and the matrix has almost zero permeability and porosity. In this study, we present a mesh-free semianalytical solution for pressure transient behavior in a 2D infinite reservoir containing a network of discrete and/or connected finite- and infinite-conductivity fractures. The proposed solution methodology is based on an analytical-element method and thus can be easily extended to incorporate other reservoir features such as sealing or leaky faults, domains with altered petrophysical properties (for example, fluid permeability or reservoir porosity), and complicated reservoir boundaries. It is shown that the pressure behavior of discretely fractured reservoirs is considerably different from the well-known Warren and Root dual-porosity reservoir model behavior. The pressure behavior of discretely fractured reservoirs shows many different flow regimes depending on fracture distribution, its intensity and conductivity. In some cases, they also exhibit a dual-porosity reservoir model behavior.     

A numerical approach for pressure transient analysis of a vertical well with complex fractures

A new well test model for a vertical fractured well is developed based on a discrete-fracture model in which the fractures are discretized as one dimensional (1-D) entities. The model overcomes the weakness of complex meshing, a large number of grids, and instability in conventional stripe-fracture models. Then, the discrete-fracture model is implemented using a hybrid element finite-element method. Triangular elements are used for matrix and line elements for the fractures. The finite element formulation is validated by comparing with the semi-analytical solution of a single ver-tical fractured well. The accuracy of the approach is shown through several examples with different fracture apertures, fracture conductivity, and fracture amount. Results from the discrete-fracture model agree reasonably well with the stripe-fracture model and the analytic solutions. The advantages of the discrete-fracture model are presented in mesh gen-eration, computational improvement, and abilities to handle complex fractures like wedge-shaped fractures and fractures with branches. Analytical results show that the number of grids in the discrete-fracture model is 10%less than stripe-fracture model, and computational efficiency increases by about 50%. The more fractures there are, the more the com-putational efficiency increases.     

Upscaling in Vertically Fractured Oil Reservoirs Using Homogenization

Flow modeling in fractured reservoirs is largely confined to the so-called sugar cube model. Here, however, we consider vertically fractured reservoirs, i.e., the situation that the reservoir geometry can be approximated by fractures enclosed columns running from the base rock to the cap rock (aggregated columns). This article deals with the application of the homogenization method to derive an upscaled equation for fractured reservoirs with aggregated columns. It turns out that vertical flow in the columns plays an important role, whereas it can be usually disregarded in the sugar cube model. The vertical flow is caused by coupling of the matrix and fracture pressure along the vertical faces of the columns. We formulate a fully implicit three-dimensional upscaled numerical model. Furthermore, we develop a computationally efficient numerical approach. As found previously for the sugar cube model, the Peclet number, i.e., the ratio between the capillary diffusion time in the matrix and the residence time of the fluids in the fracture, plays an important role. The gravity number plays a secondary role. For low Peclet numbers, the results are sensitive to gravity, but relatively insensitive to the water injection rate, lateral matrix column size, and reservoir geometry, i.e., sugar cube versus aggregated column. At a low Peclet number and sufficiently low gravity number, the effective permeability model gives good results, which agree with the solution of the aggregated column model. However, ECLIPSE simulations (Barenblatt or Warren and Root (BWR) approach) show deviations at low Peclet numbers, but show good agreement at intermediate Peclet numbers. At high Peclet numbers, the results are relatively insensitive to gravity, but sensitive to the other conditions mentioned above. The ECLIPSE simulations and the effective permeability model show large deviations from the aggregated column model at high Peclet numbers. We conclude that at low Peclet numbers, it is advantageous to increase the water injection rate to improve the net present value. However, at high Peclet numbers, increasing the flow rate may lead to uneconomical water cuts.     

多尺度嵌入式离散裂缝模型模拟方法

天然裂缝性油藏和人工压裂油藏内裂缝形态多样,分布复杂,传统的离散裂缝模型将裂缝作为基岩网格的边界,采用非结构化网格进行网格划分,其划分过程复杂,计算量大。嵌入式离散裂缝模型划分网格时不需要考虑油藏内的裂缝形态,只需对基岩系统进行简单的网格剖分,可以大大降低网格划分的复杂度,从而提高计算效率。然而,在油藏级别的数值模拟和人工压裂裂缝下的产能分析中,仍然存在计算量巨大、模拟时间过长的问题。本文提出嵌入式离散裂缝模型的多尺度数值计算格式,使用多尺度模拟有限差分法研究嵌入式离散裂缝模型渗流问题。通过在粗网格上求解局部流动问题计算多尺度基函数,多尺度基函数可以捕捉裂缝与基岩间的相互关系,反映单元内的非均质性,因此该方法既有传统尺度升级法的计算效率,又可以保证计算精度,数值结果表明这是一种有效的裂缝性油藏数值模拟方法。     

Modeling of Heat and Mass Transfer in a Fractured Porous Medium

While fractured formations are possibly the most important contributors to the production of oil worldwide, modeling fractured formations with rigorous treatments has eluded reservoir engineers in the past. To date, one of the most commonly used fractured reservoir models remains the one that was suggested by Warren and Root nearly four decades ago. In this paper, a new model for fractures embedded in a porous medium is proposed. The model considers the Navier-Stokes equation in the fracture (channel flow) while using the Brinkman equation for the porous medium. Unlike the previous approach, the proposed model does not require the assumption of orthogonality of the fractures (sugar cube assumption) nor does it impose incorrect boundary conditions for the interface between the fracture and the porous medium. Also, the transfer coefficient between the fracture and matrix interface does not need to be specified, unlike the cases for which Darcy's law is used. In order to demonstrate the usefulness of the approach, a two-dimensional model of a fractured formation is developed and numerical simulation runs conducted.

The proposed model is derived through a series of finite element modeling runs for various cases using the Navier-Stokes equation in the channel while maintaining the Brinkman equation in the porous medium. Various cases studied include different fracture orientations, fracture frequencies, and thermal and solutal constraints. The usefulness of the proposed model in modeling complex formations is discussed. Finally, a series of numerical runs also provided validity of the proposed model for the cases in which thermal and solutal effects are important. Such a study of double diffusive phenomena, coupled with forced convection, in the context of fractured formations has not been reported before.     

An efficient semi‐implicit subgrid method for free‐surface flows on hierarchical grids

We present a new modelling strategy for improving the efficiency of computationally intensive flow problems in environmental free‐surface flows. The approach combines a recently developed semi‐implicit subgrid method with a hierarchical grid solution strategy. The method allows the incorporation of high‐resolution data on subgrid scale to obtain a more accurate and efficient hydrodynamic model. The subgrid method improves the efficiency of the hierarchical grid method by providing better solutions on coarse grids. The method is applicable to both steady and unsteady flows, but we particularly focus on river flows with steady boundary conditions. There, the combined hierarchical grid–subgrid method reduces the computational effort to obtain a steady state with factors up to 43. For unsteady models, the method can be used for efficiently generating accurate initial conditions on high‐resolution grids. Additionally, the method provides automatic insight in grid convergence. We demonstrate the efficiency and applicability of the method using a schematic test for the vortex shedding around a circular cylinder and a real‐world river case study. Copyright © 2015 John Wiley & Sons, Ltd.     

An unstructured-grid numerical model for interfacial multiphase fluids based on multi-moment finite volume formulation and THINC method

We present a practical numerical framework for incompressible interfacial multiphase flows on unstructured grids with arbitrary and hybrid elements. The numerical framework is constructed by combining VPM (volume-average/point-value multi-moment) and UMTHINC (unstructured multi-dimensional tangent of hyperbola interface capturing) schemes. To facilitate accurate and reliable simulations for interfacial multiphase flows on arbitrary and hybrid unstructured grids, we have made the following major new efforts in this work. (1) UMTHINC scheme on prismatic and pyramidal elements to facilitate computations on hybrid arbitrary unstructured grids; (2) Consistent numerical formulation for mass and momentum transports to simulate multiphase flows of large density ratio; (3) Combined FVM-FEM for accurate solution to diffusion equation; (4) Pressure-projection formulation in consistent with the balanced-force model. Integrating all these numerical techniques effectively enhances the accuracy and robustness in interface capturing and numerical solution of multiphase fluid dynamics, which results in a numerical framework of great significance for practical applications. Numerical verifications have been carried out through benchmark tests ranging from surface tension dominant flows of small scale to large scale flows with violently-changing interfaces. Numerical results demonstrate that the present framework is robust with adequate accuracy for simulating multiphase flows in complex geometries.     

Effects of Fracture Boundary Conditions on Matrix-fracture Transfer Shape Factor

The matrix-fracture transfer shape factor is one of the important parameters in modeling naturally fractured reservoirs. Four decades after Warren and Root (1963, SPEJ, 245–255.) introduced the double porosity concept and suggested a relation for it, this parameter is still not completely understood. Even for a single-phase flow problem, investigators report different shape factors. This study shows that for a single-phase flow in a particular matrix block, the shape factor that Warren and Root defined is not unique and depends on the pressure in the fracture and how it changes with time. We use the Laplace domain analytical solutions of the diffusivity equation for different geometries and different boundary conditions to show that the shape factor depends on the fracture pressure change with time. In particular, by imposing a constant fracture pressure as it is typically done, one obtains the shape factor that Lim and Aziz (1995, J. Petrolean Sci. Eng. 13, 169.) calculated. However, other shape factors, similar to those reported in other studies are obtained, when other boundary conditions are chosen. Although, the time variability of the boundary conditions can be accounted for by the Duhamel’s theorem, in practice using large time-steps in numerical simulations can potentially introduce large errors in simulation results. However, numerical simulation models make use of a stepwise approximation of this theorem. It is shown in this paper that this approximation could lead to large errors in matrix-fracture transfer rate if large time-steps are chosen.     

Up-scaling Flow in Fractured Media: Equivalence Between the Large Scale Averaging Theory and the Continuous Time Random Walk Method

In a recent paper, (McNabb, 1978), we set up a method allowing to compute both the transient and steady-state exchange terms between the matrix and fractured regions of a naturally fractured porous medium using the continuous time random walk method (CTRW). In particular, the exchange coefficient parametrizing the large-scale exchange term was computed on physical grounds using a time integration of the so-called time correlation function corresponding to the particle presence in the fractures. On the other hand, the large scale averaging theory (LSAT) as developed by Quintard and Whitaker (Quintard and Whitaker, 1996) gives another definition for this exchange coefficient . It also provides a so-called closure problem allowing to compute from the solution of a well-defined steady state boundary value problem, to be solved over a representative volume of the high resolution fractured map. The goal of the present paper is to show analytically that both definitions coincide, yielding a unique and well defined value of the coefficient. This provides an unification of two approaches whose respective backgrounds are very different at the first glance.     

Modeling Free-Surface Seepage Flow in Complicated Fractured Rock Mass Using a Coupled RPIM-FEM Method

The prediction of the free-surface seepage flow behavior in fractured rock mass is of significance in geotechnical engineering. There are two major issues in solving the seepage flow in complicated fractured rock mass based on the fractured porous medium (FPM) flow model, in which groundwater is assumed to flow simultaneously in both rock matrix and embedded fractures: One is the mesh generation of rock mass in the presence of the fracture network, especially when there exist a large number of stochastic fractures; the other is that a robust iteration algorithm is required since the free surface is unknown at the beginning of solution. Aiming at these two issues, this paper proposes a novel numerical method by coupling radial point interpolation method (RPIM) and finite element method (FEM), in which RPIM is utilized to model the rock matrix and FEM is utilized to model the fractures. On the basis of the variational inequality (VI) theory for free-surface seepage analysis, the computation formulations of the numerical method are derived and the corresponding computation program is developed. Three examples are solved with the present method. It is found that the VI theory can be extended to solve the free-surface seepage problem based on the FPM flow model. A crucial advantage of the present method is that the mesh generation can be greatly simplified. The present method has been verified to be a robust, efficient and reliable method for modeling the groundwater flow in complicated fractured rock mass.     

Assessment of ghost-cell based cut-cell method for large-eddy simulations of compressible flows at high Reynolds number

Large-eddy simulations (LES) of high Reynolds number flows are performed using a non-body conformal method in conjunction with a wall model. We use a simple wall function to model the wall-shear stress and the truncation error of the numerical discretization to model the sub-grid scale turbulence (implicit LES), although these can be easily replaced if necessary. The validation cases are: turbulent flow through an inclined channel, turbulent flow over a wavy surface, and supersonic flow over a circular cylinder. Since the near-wall grids are naturally coarse, the key is to use a method that is capable of capturing the flow dynamics accurately in the vicinity of the interface. Towards the purpose, we develop a Cartesian cut-cell method, referred to as the ghost-cell based cut-cell method (GC-CCM), in the context of fully compressible solutions of Navier–Stokes equations. This method employs ghost-cells inside the solid interface such that the local spatial reconstruction remains consistent everywhere including in the vicinity of the boundary. In order to capture the near-wall flow behavior more accurately with coarse grids, this method decomposes cell faces of merged cells and computes fluxes through each decomposed segment separately. The objective of this work is to qualify whether the proposed method can accurately represent the high Reynolds number flows in the vicinity of immersed interfaces. To analyze the performance of the proposed method, we compare the results to the corresponding numerical results from the two other non-body conformal methods, namely the ghost-cell based immersed boundary method (GCIBM) and standard cut-cell method (S-CCM), that are implemented in the same numerical solver. The comparison demonstrates that the proposed method is capable of capturing near-wall flows relatively accurately with coarse grids.     

A finite-volume ELLAM for non-linear flux convection-diffusion problems

In this paper, a modified finite-volume Eulerian-Lagrangian localized adjoint method (FVELLAM) extended for convection-diffusion problems with a non-linear flux function is introduced. Tracking schemes are discussed using viscous Burgers’ equation. It is shown that in order to have smooth results, only the new time level values should be used in tracking process. Then, the proposed method is employed to study immiscible incompressible two-phase flows in porous media. Various one- and two-dimensional test cases involving internal sources and sinks are solved and accuracy of solution and performance of the method are investigated by comparing the results obtained using FVELLAM with those of fine grid solutions. Finally, it is concluded that although proposed FVELLAM produces satisfactory results even on coarse grids and allows fairly large time-steps, its major advantage is in solving convection-dominated problems in heterogeneous media.     

Matrix-fracture exchange in a fractured porous medium: stochastic upscalingCoefficient d'échange matrice-fracture en milieu poreux fracturé : changement d'échelle par approche stochastique

We present in this Note a stochastic approach to the matrix-fracture exchange in a heterogeneous fractured porous medium. We introduce an intermediate scale, called the unit-scale, between the local-scale (fracture-scale) and the large-scale characteristic of the reservoir mesh (reservoir block). This paper focuses on the problem of upscaling fluid exchange phenomena from the unit scale to the reservoir mesh or block scale. Simplifying the Darcian flow terms enables us to obtain a probabilistic solution of the dual continuum problem, in continuous time, in the case of a purely random exchange coefficient. This is then used to develop several upscaling approaches to the fluid exchange problem, and to analyze the so-called ‘effective’ exchange coefficient. The results are a first contribution to the more general problem of upscaling multidimensional flow-exchange processes in space and time, in randomly heterogeneous dual continua. To cite this article: M. Kfoury et al., C. R. Mecanique 332 (2004).     

A methodology for high performance computation of fully inhomogeneous turbulent flows

A large‐eddy simulation methodology for high performance parallel computation of statistically fully inhomogeneous turbulent flows on structured grids is presented. Strategies and algorithms to improve the memory efficiency as well as the parallel performance of the subgrid‐scale model, the factored scheme, and the Poisson solver on shared‐memory parallel platforms are proposed and evaluated. A novel combination of one‐dimensional red–black/line Gauss–Seidel and two‐dimensional red–black/line Gauss–Seidel methods is shown to provide high efficiency and performance for multigrid relaxation of the Poisson equation. Parallel speedups are measured on various shared‐distributed memory systems. Validations of the code are performed in large‐eddy simulations of turbulent flows through a straight channel and a square duct. Results obtained from the present solver employing a Lagrangian dynamic subgrid‐scale model show good agreements with other available data. The capability of the code for more complex flows is assessed by performing a large‐eddy simulation of the tip‐leakage flow in a linear cascade. Copyright © 2006 John Wiley & Sons, Ltd.     

Multiphase computational method for thermal interactions between compressible fluid and arbitrarily shaped solids

A numerical method was investigated for multiphase fields consisting of compressible gas and arbitrarily shaped solids. Since the proposed model is based on a one‐fluid model in which variables are averaged according to the phase fractions in the computational cells; it enables us to estimate gas‐solid momentum and thermal interactions without setting up adapting grids even if the solids have extremely complicated shapes. The governing equations are derived with the characteristics of an ideal gas assuming the specific heat to be uniform in the multiphase field. The derived equations in conservative form are discretized with a finite volume method. In addition, the pressure is calculated implicitly in a similar way to incompressible flow solvers. Because of these improvements, the proposed method allows us to calculate low Mach number compressible flows free from the Courant‐Friedrichs‐Lewy condition based on the speed of sound and to conserve the mass more accurately. To confirm the validity of the proposed method, it was applied to natural convection around an isothermal cylinder and a heat‐conducting pipe. In comparison with previous studies, it was confirmed that the gas flows and temperature distributions are predicted reasonably. In addition, a numerical experiment was conducted under more complicated conditions, namely, gas leaking from a container including heat sources. As a result, it was demonstrated that the proposed method enables us to predict unsteady variations of pressure and temperature distributions in the container due to the leakage while still conserving mass accurately.     

An efficient method for discretizing 3D fractured media for subsurface flow and transport simulations

We introduce a new method to discretize inclined non‐planar two‐dimensional (2D) fractures in three‐dimensional (3D) fractured media for subsurface flow and transport simulations. The 2D fractures are represented by ellipsoids. We first discretize the fractures and generate a 2D finite element mesh for each fracture. Then, the mesh of fractures is analyzed by searching and treating critical geometric configurations. Based on that search, the method generates a quality mesh and allows for including finer grids. A solute transport problem in fractured porous media is solved to test the method. The results show that the method (i) adequately represents the fractured domain by maintaining the geometric integrity of input surfaces and geologic data, (ii) provides accurate results for both simple and complex fractured domains, (iii) is insensitive to spatial discretization, and (iv) is computationally very efficient. For inclined and vertical fractures, analytical and numerical solutions are shown to be in good agreement. The method is therefore suitable to discretize fracture networks for flow and transport simulations in fractured porous media. Copyright © 2010 John Wiley & Sons, Ltd.     

A Direct Determination of the Transient Exchange Term of Fractured Media Using a Continuous Time Random Walk Method

In two recent papers, Ntinger and Estébenet, 2000; Ntinger et al., submitted, we set-up a method allowing to compute both the transient and steady-state exchange terms between the matrix and fractured regions of a naturally fractured porous medium using continuous time random walk methods (CTRW). The goal of the present paper is to show that a new version of the CTRW algorithm provides a direct determination of the so called transient exchange function f(t) (or its Laplace transform f(s)) widely used in well test interpretation. It is shown that this function is directly linked with the probability density of the first escape time in the fractured region of a Brownian particle launched initially in the matrix region. This new interpretation allows relating directly the exchange coefficient with the mean escape time of brownian particles in the matrix. From a practical point of view, these new results allow to derive a simpler version of the CTRW method. In addition, we obtain a considerable speed up of the CTRW method for up-scaling fractured reservoirs.     

Discretizing two‐dimensional complex fractured fields for incompressible two‐phase flow

A method is introduced to discretize irregular and complex two‐dimensional fractured media. The geometry of the fractured media is first analysed by searching and treating the complex configurations. Based on that, the method generated a good mesh quality and allows for including finer grids. An incompressible two‐phase flow problem is solved to compare the developed method and a public method based on the approximation of a 1D fracture by the edges of a 2D finite element grid of the porous media. The comparison showed that the developed method (i) represents better the fractured domain by maintaining the geometric integrity of input surfaces and geologic data, (ii) provides, for sample and complex fractured domains, excellent and more accurate results, and (iii) is much less sensitive to the grid sizes. Furthermore, the method has to be more efficient than the other methods for transport problems and has to provide better predictable results; this is mainly based on point (ii) and because the method produces optimal triangular grids. Copyright © 2009 John Wiley & Sons, Ltd.     

裂隙岩体水力特征数值模拟试验的初步结果

本文利用文献中的数据进行数值模拟试验,结果表明：流量等效和渗透作用力变形效果等效互不包容;不同级别结构面体系难以用连续介质方法处理,同级别结构面体系在工程作用力变化范围内,可由其静态水力特征推演动态特征;工程岩体结构面间充分连通的假设过于粗糙,常规试验手段获取的量值倾向于局部特征的反映。     

Numerical implementation of Aristov–Pukhnachev's formulation for axisymmetric viscous incompressible flows

In the present work a finite‐difference technique is developed for the implementation of a new method proposed by Aristov and Pukhnachev (Doklady Phys. 2004; 49 (2):112–115) for modeling of the axisymmetric viscous incompressible fluid flows. A new function is introduced that is related to the pressure and a system similar to the vorticity/stream function formulation is derived for the cross‐flow. This system is coupled to an equation for the azimuthal velocity component. The scheme and the algorithm treat the equations for the cross‐flow as an inextricably coupled system, which allows one to satisfy two conditions for the stream function with no condition on the auxiliary function. The issue of singularity of the matrix is tackled by adding a small parameter in the boundary conditions. The scheme is thoroughly validated on grids with different resolutions. The new numerical tool is applied to the Taylor flow between concentric rotating cylinders when the upper and lower lids are allowed to rotate independently from the inner cylinder, while the outer cylinder is held at rest. The phenomenology of this flow is adequately represented by the numerical model, including the hysteresis that takes place near certain specific values of the Reynolds number. Thus, the present results can be construed to demonstrate the viability of the new model. The success can be attributed to the adequate physical nature of the auxiliary function. The proposed technique can be used in the future for in‐depth investigations of the bifurcation phenomena in rotating flows. Copyright © 2009 John Wiley & Sons, Ltd.