Comparison of Three FE-FV Numerical Schemes for Single- and Two-Phase Flow Simulation of Fractured Porous Media

We benchmark a family of hybrid finite element–node-centered finite volume discretization methods (FEFV) for single- and two-phase flow/transport through porous media with discrete fracture representations. Special emphasis is placed on a new method we call DFEFVM in which the mesh is split along fracture–matrix interfaces so that discontinuities in concentration or saturation can evolve rather than being suppressed by nodal averaging of these variables. The main objective is to illustrate differences among three discretization schemes suitable for discrete fracture modeling: (a) FEFVM with volumetric finite elements for both fractures and porous rock matrix, (b) FEFVM with lower dimensional finite elements for fractures and volumetric finite elements for the matrix, and (c) DFEFVM with a mesh that is split along material discontinuities. Fracture discontinuities strongly influence single- and multi-phase fluid flow. Continuum methods, when used to model transport across such interfaces, smear out concentration/saturation. We show that the new DFEFVM addresses this problem producing significantly more accurate results. Sealed and open single fractures as well as a realistic fracture geometry are used to conduct tracer and water-flooding numerical experiments. The benchmarking results also reveal the limitations/mesh refinement requirements of FE node-centered FV hybrid methods. We show that the DFEFVM method produces more accurate results even for much coarser meshes.     

Fracture-Flow-Enhanced Matrix Diffusion in Solute Transport Through Fractured Porous Media

Over the past few decades, significant progress of assessing chemical transport in fractured rocks has been made in laboratory and field investigations as well as in mathematic modeling. In most of these studies, however, matrix diffusion on fracture–matrix surfaces is considered as a process of molecular diffusion only. Mathematical modeling based on this traditional concept often had problems in explaining or predicting tracer transport in fractured rock. In this article, we propose a new conceptual model of fracture-flow-enhanced matrix diffusion, which correlates with fracture-flow velocity. The proposed model incorporates an additional matrix-diffusion process, induced by rapid fluid flow along fractures. According to the boundary-layer theory, fracture-flow-enhanced matrix diffusion may dominate mass-transfer processes at fracture–matrix interfaces, where rapid flow occurs through fractures. The new conceptual model can be easily integrated with analytical solutions, as demonstrated in this article, and numerical models, as we foresee. The new conceptual model is preliminarily validated using laboratory experimental results from a series of tracer breakthrough tests with different velocities in a simple fracture system. Validating of the new model with field experiments in complicated fracture systems and numerical modeling will be explored in future research.     

A dual-porosity model for gas reservoir flow incorporating adsorption behaviour—Part II. Numerical algorithm and example applications

In the present study, an algorithm is presented for the dual-porosity model formulated in Part I of this series. The resultant flow equation with the dual-porosity formulation is of an integro-(partial) differential equation involving differential terms for the Darcy flow in large fractures and integrals in time for diffusion within matrix blocks. The algorithm developed here to solve this equation involves a step-by-step finite difference procedure combined with a quadrature scheme. The quadrature scheme, used for the integral terms, is based on the trapezoidal method which is of second-order precision. This order of accuracy is consistent with the first- and second-order finite difference approximations used here to solve the differential terms in the derived flow equation. In an approach consistent with many petroleum reservoir and groundwater numerical flow models, the example formulation presented uses a first-order implicit algorithm. A two-dimensional example is also demonstrated, with the proposed model and numerical scheme being directly incorporated into the commercial gas reservoir simulator SIMED II that is based on a fully implicit finite difference approach. The solution procedure is applied to several problems to demonstrate its performance. Results from the derived dual-porosity formulation are also compared to the classic Warren–Root model. Whilst some of this work confirmed previous findings regarding Warren–Root inaccuracies at early times, it was also found that inaccuracy can re-enter the Warren–Root results whenever there are changes in boundary conditions leading to transient variation within the domain.     

Nonlinear Multigrid Methods for Numerical Solution of the Variably Saturated Flow Equation in Two Space Dimensions

The need of accurate and efficient numerical schemes to solve Richards’ equation is well recognized. This study is carried out to examine the numerical performances of the nonlinear multigrid method for numerical solving of the two-dimensional Richards’ equation modeling water flow in variably saturated porous media. The numerical approach is based on an implicit, second-order accurate time discretization combined with a vertex centered finite volume method for spatial discretization. The test problems simulate infiltration of water in 2D saturated–unsaturated soils with hydraulic properties described by van Genuchten–Mualem models. The numerical results obtained are compared with those provided by the modified Picard–preconditioned conjugated gradient (Krylov subspace) approach.     

Nonlinear Multigrid Methods for Numerical Solution of the Unsaturated Flow Equation in Two Space Dimensions

Picard and Newton iterations are widely used to solve numerically the nonlinear Richards’ equation (RE) governing water flow in unsaturated porous media. When solving RE in two space dimensions, direct methods applied to the linearized problem in the Newton/Picard iterations are inefficient. The numerical solving of RE in 2D with a nonlinear multigrid (MG) method that avoids Picard/Newton iterations is the focus of this work. The numerical approach is based on an implicit, second-order accurate time discretization combined with a second-order accurate finite difference spatial discretization. The test problems simulate infiltration of water in 2D unsaturated soils with hydraulic properties described by Broadbridge–White and van Genuchten–Mualem models. The numerical results show that nonlinear MG deserves to be taken into consideration for numerical solving of RE.     

A parallel implicit domain decomposition algorithm for the large eddy simulation of incompressible turbulent flows on 3D unstructured meshes

We present a parallel fully implicit algorithm for the large eddy simulation (LES) of incompressible turbulent flows on unstructured meshes in three dimensions. The LES governing equations are discretized by a stabilized Galerkin finite element method in space and an implicit second-order backward differentiation scheme in time. To efficiently solve the resulting large nonlinear systems, we present a highly parallel Newton-Krylov-Schwarz algorithm based on domain decomposition techniques. Analytic Jacobian is applied in order to obtain the best achievable performance. Two benchmark problems of lid-driven cavity and flow passing a square cylinder are employed to validate the proposed algorithm. We then apply the algorithm to the LES of turbulent flows passing a full-size high-speed train with realistic geometry and operating conditions. The numerical results show that the algorithm is both accurate and efficient and exhibits a good scalability and parallel efficiency with tens of millions of degrees of freedom on a computer with up to 4096 processors. To understand the numerical behavior of the proposed fully implicit scheme, we study several important issues, including the choices of linear solvers, the overlapping size of the subdomains, and, especially, the accuracy of the Jacobian matrix. The results show that an exact Jacobian is necessary for the efficiency and the robustness of the proposed LES solver.     

Upscaled Unstructured Computational Grids for Efficient Simulation of Flow in Fractured Porous Media

Discrete fracture modeling (DFM) is currently the most promising approach for modeling of naturally fractured reservoirs and simulation of multiphase fluid flow therein. In contrast with the classical double-porosity/double permeability models, in the DFM approach all the interactions and fluid flow in and between the fractures and within the matrix are modeled in a unified manner, using the same computational grid. There is no need for computing the shape factors, which are crucial to the accuracy of the double-porosity models. We have exploited this concept in order to develop a new method for the generation of unstructured computational grids. In the new approach the geological model (GM) of the reservoir is first generated, using square or cubic grid blocks. The GM is then upscaled using a method based on the multiresolution wavelet transformations that we recently developed. The upscaled grid contains a distribution of the square or cubic blocks of various sizes. A map of the blocks’ centers is then used with an optimized Delauney triangulation method and the advancing-front technique, in order to generate the final unstructured triangulated grid suitable for use in any general reservoir simulator with any number of fluid phases. The new method also includes an algorithm for generating fractures that, contrary to the previous methods, does not require modifying their paths due to the complexities that may arise in spatial distribution of the grid blocks. It also includes an effective partitioning of the simulation domain that results in large savings in the computation times. The speed-up in the computations with the new upscaled unstructured grid is about three orders of magnitude over that for the initial GM. Simulation of waterflooding indicates that the agreement between the results obtained with the GM and the upscaled unstructured grid is excellent. The method is equally applicable to the simulations of multiphase flow in unfractured, but highly heterogeneous, reservoirs.     

Explicit Rate–Time Solutions for Modeling Matrix–Fracture Flow of Single Phase Gas in Dual-Porosity Media

One of the widely used methods for modeling matrix–fracture fluid exchange in naturally fractured reservoirs is dual porosity approach. In this type of modeling, matrix blocks are regarded as sources/sinks in the fracture network medium. The rate of fluid transfer from matrix blocks into fracture medium may be modeled using shape factor concept (Warren and Root, SPEJ 3:245–255, 1963); or the rate–time solution is directly derived for the specific matrix geometry (de Swaan, SPEJ 16:117–122, 1976). Numerous works have been conducted to study matrix–fracture fluid exchange for slightly compressible fluids (e.g. oil). However, little attention has been taken to systems containing gas (compressible fluid). The objective of this work is to develop explicit rate–time solutions for matrix–fracture fluid transfer in systems containing single phase gas. For this purpose, the governing equation describing flow of gas from matrix block into fracture system is linearized using pseudopressure and pseudotime functions. Then, the governing equation is solved under specific boundary conditions to obtain an implicit relation between rate and time. Since rate calculations using such an implicit relation need iterations, which may be computationally inconvenient, an explicit rate–time relation is developed with the aid of material balance equation and several specific assumptions. Also, expressions are derived for average pseudopressure in matrix block. Furthermore, simplified solutions (originated from the complex general solutions) are introduced applicable in infinite and finite acting flow periods in matrix. Based on the derived solutions, expressions are developed for shape factor. An important observation is that the shape factor for gas systems is the same as that of oil bearing matrix blocks. Subsequently, a multiplier is introduced which relates rate to matrix pressure instead of matrix pseudopressure. Finally, the introduced equations are verified using a numerical simulator.     

Time-Dependent Shape Factors for Uniform and Non-Uniform Pressure Boundary Conditions

Matrix–fracture transfer functions are the backbone of any dual-porosity or dual-permeability formulation. The chief feature within them is the accurate definition of shape factors. To date, there is no completely accepted formulation of a matrix–fracture transfer function. Many formulations of shape factors for instantly-filled fractures with uniform pressure distribution have been presented and used; however, they differ by up to five times in magnitude. Based on a recently presented transfer function, time-dependent shape factors for water imbibing from fracture to matrix under pressure driven flow are proposed. Also new matrix–fracture transfer pressure-based shape factors for instantly-filled fractures with non-uniform pressure distribution are presented in this article. These are the boundary conditions for a case for porous media with clusters of parallel and disconnected fractures, for instance. These new pressure-based shape factors were obtained by solving the pressure diffusivity equation for a single phase using non-uniform boundary conditions. This leads to time-dependent shape factors because of the transient part of the solution for pressure. However, approximating the solution with an exponential function, one obtains constant shape factors that can be easily implemented in current dual-porosity reservoir simulators. The approximate shape factors provide good results for systems where the transient behavior of pressure is short (a case commonly encountered in fractured reservoirs).     

An adaptive local grid refinement and peak/valley capture algorithm to solve nonlinear transport problems with moving sharp-fronts

Highly nonlinear advection–dispersion-reaction equations govern numerous transport phenomena. Robust, accurate, and efficient algorithms to solve these equations hold the key to the success of applying numerical models to field problems. This paper presents the development and verification of a computational algorithm to approximate the highly nonlinear transport equations of reactive chemical transport and multiphase flow. The algorithm was developed based on the Lagrangian-Eulerian decoupling method with an adaptive ZOOMing and Peak/valley Capture (LEZOOMPC) scheme. It consists of both backward and forward node tracking, rough element determination, peak/valley capturing, and adaptive local grid refinement. A second-order tracking was implemented to accurately and efficiently track all fictitious particles. Shanks’ method was introduced to deal with slowly converging case. The accuracy and efficiency of the algorithm were verified with the Burger equation for a variety of cases. The robustness of the algorithm to achieve convergent solutions was demonstrated by highly nonlinear reactive contaminant transport and multiphase flow problems.     

Modeling of Solute Transport in a 3D Rough-Walled Fracture–Matrix System

Fluid flow and solute transport in a 3D rough-walled fracture–matrix system were simulated by directly solving the Navier–Stokes equations for fracture flow and solving the transport equation for the whole domain of fracture and matrix with considering matrix diffusion. The rough-walled fracture–matrix model was built from laser-scanned surface tomography of a real rock sample, by considering realistic features of surfaces roughness and asperity contacts. The numerical modeling results were compared with both analytical solutions based on simplified fracture surface geometry and numerical results by particle tracking based on the Reynolds equation. The aim is to investigate impacts of surface roughness on solute transport in natural fracture–matrix systems and to quantify the uncertainties in application of simplified models. The results show that fracture surface roughness significantly increases heterogeneity of velocity field in the rough-walled fractures, which consequently cause complex transport behavior, especially the dispersive distributions of solute concentration in the fracture and complex concentration profiles in the matrix. Such complex transport behaviors caused by surface roughness are important sources of uncertainty that needs to be considered for modeling of solute transport processes in fractured rocks. The presented direct numerical simulations of fluid flow and solute transport serve as efficient numerical experiments that provide reliable results for the analysis of effective transmissivity as well as effective dispersion coefficient in rough-walled fracture–matrix systems. Such analysis is helpful in model verifications, uncertainty quantifications and design of laboratorial experiments.     

A 3D Computational Study of Effective Medium Methods Applied to Fractured Media

This work evaluates and improves upon existing effective medium methods for permeability upscaling in fractured media. Specifically, we are concerned with the asymmetric self-consistent, symmetric self-consistent, and differential methods. In effective medium theory, inhomogeneity is modeled as ellipsoidal inclusions embedded in the rock matrix. Fractured media correspond to the limiting case of flat ellipsoids, for which we derive a novel set of simplified formulas. The new formulas have improved numerical stability properties, and require a smaller number of input parameters. To assess their accuracy, we compare the analytical permeability predictions with three-dimensional finite-element simulations. We also compare the results with a semi-analytical method based on percolation theory and curve-fitting, which represents an alternative upscaling approach. A large number of cases is considered, with varying fracture aperture, density, matrix/fracture permeability contrast, orientation, shape, and number of fracture sets. The differential method is seen to be the best choice for sealed fractures and thin open fractures. For high-permeable, connected fractures, the semi-analytical method provides the best fit to the numerical data, whereas the differential method breaks down. The two self-consistent methods can be used for both unconnected and connected fractures, although the asymmetric method is somewhat unreliable for sealed fractures. For open fractures, the symmetric method is generally the more accurate for moderate fracture densities, but only the asymmetric method is seen to have correct asymptotic behavior. The asymmetric method is also surprisingly accurate at predicting percolation thresholds.     

Implicit schemes with large time step for non‐linear equations: application to river flow hydraulics

In this work, first‐order upwind implicit schemes are considered. The traditional tridiagonal scheme is rewritten as a sum of two bidiagonal schemes in order to produce a simpler method better suited for unsteady transcritical flows. On the other hand, the origin of the instabilities associated to the use of upwind implicit methods for shock propagations is identified and a new stability condition for non‐linear problems is proposed. This modification produces a robust, simple and accurate upwind semi‐explicit scheme suitable for discontinuous flows with high Courant–Friedrichs–Lewy (CFL) numbers. The discretization at the boundaries is based on the condition of global mass conservation thus enabling a fully conservative solution for all kind of boundary conditions. The performance of the proposed technique will be shown in the solution of the inviscid Burgers' equation, in an ideal dambreak test case, in some steady open channel flow test cases with analytical solution and in a realistic flood routing problem, where stable and accurate solutions will be presented using CFL values up to 100. Copyright © 2004 John Wiley & Sons, Ltd.     

A New Coal-Permeability Model: Internal Swelling Stress and Fracture–Matrix Interaction

We have developed a new coal-permeability model for uniaxial strain and constant confining-stress conditions. The model is unique in that it explicitly considers fracture–matrix interaction during coal-deformation processes and is based on a newly proposed internal swelling stress concept. This concept is used to account for the impact of matrix swelling (or shrinkage) on fracture-aperture changes resulting from partial separation of matrix blocks by fractures that do not completely cut through the whole matrix. The proposed permeability model is evaluated using data from three Valencia Canyon coalbed wells in the San Juan Basin, where increased permeability has been observed during CH₄ gas production, as well as using published data from laboratory tests. Model results are generally in good agreement with observed permeability changes. The importance of fracture–matrix interaction in determining coal permeability, demonstrated in this study using relatively simple stress conditions, underscores the need for a dual-continuum (fracture and matrix) mechanical approach to rigorously capture coal-deformation processes under complex stress conditions, as well as the coupled flow and transport processes in coal seams.     

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

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

Experimental study and mathematical simulation of the characteristics of a turbulent flow in a straight circular pipe rotating about its longitudinal axis

The vorticity formed in the cross section of a turbulent flow in a straight circular pipe rotating about its longitudinal axis decreases the values of the turbulent stresses, turbulence energy, and dissipation rate along the pipe. The results of laboratory experiments and calculations by the second-order closure model of turbulent transfer are presented. On the whole, the model using a system of transport equations yields better agreement with experimental data than the models with algebraic relations for second-order moments. Translated from Prikladnaya Mekhanika i Tekhnicheskaya Fizika, Vol. 39, No. 2, pp. 103–116, March–April, 1998.     

An efficient implicit unstructured finite volume solver for generalised Newtonian fluids

An implicit finite volume solver is developed for the steady-state solution of generalised Newtonian fluids on unstructured meshes in 2D. The pseudo-compressibility technique is employed to couple the continuity and momentum equations by transforming the governing equations into a hyperbolic system. A second-order accurate spatial discretisation is provided by performing a least-squares gradient reconstruction within each control volume of unstructured meshes. A central flux function is used for the convective terms and a solution jump term is added to the averaged component for the viscous terms. Global implicit time-stepping using successive evolution–relaxation is utilised to accelerate the convergence to steady-state solutions. The performance of our flow solver is examined for power-law and Carreau–Yasuda non-Newtonian fluids in different geometries. The effects of model parameters and Reynolds number are studied on the convergence rate and flow features. Our results verify second-order accuracy of the discretisation and also fast and efficient convergence to the steady-state solution for a wide range of flow variables.     

Modeling of Two-Phase Flow in Fractured Porous Media on Unstructured Non-Uniformly Coarsened Grids

In this article, we investigate two strategies for coarsening fractured geological models. The first approach, which generates grids that resolve the fractures, is referred to as explicit fracture-matrix separation (EFMS). The second approach is based on a non-uniform coarsening strategy introduced in Aarnes et al. (Adv Water Resour 30(11):2177–2193, 2007a). A series of two-phase flow simulations where the saturation is modeled on the respective coarse grids are performed. The accuracy of the resulting solutions is examined, and the robustness of the two strategies is assessed with respect to number of fractures, degree of coarsening, well locations, phase viscosities, and fracture permeability. The numerical results show that saturation solutions obtained on the non-uniform coarse grids are consistently more accurate than the corresponding saturation solutions obtained on the EFMS grids. The numerical results also reveal that it is much easier to tune the upscaling factor with the non-uniform coarsening approach.     

High‐order implicit Runge–Kutta time integrators for fluid‐structure interactions

This paper presents an approach to develop high‐order, temporally accurate, finite element approximations of fluid‐structure interaction (FSI) problems. The proposed numerical method uses an implicit monolithic formulation in which the same implicit Runge–Kutta (IRK) temporal integrator is used for the incompressible flow, the structural equations undergoing large displacements, and the coupling terms at the fluid‐solid interface. In this context of stiff interaction problems, the fully implicit one‐step approach presented is an original alternative to traditional multistep or explicit one‐step finite element approaches. The numerical scheme takes advantage of an arbitrary Lagrangian–Eulerian formulation of the equations designed to satisfy the geometric conservation law and to guarantee that the high‐order temporal accuracy of the IRK time integrators observed on fixed meshes is preserved on arbitrary Lagrangian–Eulerian deforming meshes. A thorough review of the literature reveals that in most previous works, high‐order time accuracy (higher than second order) is seldom achieved for FSI problems. We present thorough time‐step refinement studies for a rigid oscillating‐airfoil on deforming meshes to confirm the time accuracy on the extracted aerodynamics reactions of IRK time integrators up to fifth order. Efficiency of the proposed approach is then tested on a stiff FSI problem of flow‐induced vibrations of a flexible strip. The time‐step refinement studies indicate the following: stability of the proposed approach is always observed even with large time step and spurious oscillations on the structure are avoided without added damping. While higher order IRK schemes require more memory than classical schemes (implicit Euler), they are faster for a given level of temporal accuracy in two dimensions. Copyright © 2015 John Wiley & Sons, Ltd.