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.     

Variable penalty method for finite element analysis of incompressible flow

A new scheme is applied for increasing the accuracy of the penalty finite element method for incompressible flow by systematically varying from element to element the sign and magnitude of the penalty parameter λ, which enters through ?.v + p/λ = 0, an approximation to the incompressibility constraint. Not only is the error in this approximation reduced beyond that achievable with a constant λ, but also digital truncation error is lowered when it is aggravated by large variations in element size, a critical problem when the discretization must resolve thin boundary layers. The magnitude of the penalty parameter can be chosen smaller than when λ is constant, which also reduces digital truncation error; hence a shorter word-length computer is more likely to succeed. Error estimates of the method are reviewed. Boundary conditions which circumvent the hazards of aphysical pressure modes are catalogued for the finite element basis set chosen here. In order to compare performance, the variable penalty method is pitted against the conventional penalty method with constant λ in several Stokes flow case studies.     

One-Dimensional Matrix-Fracture Transfer in Dual Porosity Systems with Variable Block Size Distribution

Most of the developed models for fractured reservoirs assume ideal matrix block size distribution. This assumption may not be valid in reality for naturally fractured reservoirs and possibly lead to errors in prediction of production from the naturally fractured reservoirs especially during a transient period or early time production from the matrix blocks. In this study, we investigate the effect of variable block size distribution on one- dimensional flow of compressible fluids in fractured reservoirs. The effect of different matrix block size distributions on the single phase matrix-fracture transfer is studied using a recently developed semi-analytical approach. The proposed model is able to simulate fluid exchange between matrix and fracture for continuous or discrete block size distributions using probability density functions or structural information of a fractured formation. The presented semi-analytical model demonstrates a good accuracy compared to the numerical results. There have been recent attempts to consider the effect of variable block size distribution in naturally fractured reservoir modeling for slightly compressible fluids with a constant viscosity and compressibility. The main objective of this study is to consider the effect of variable block size distribution on a one-dimensional matrix-fracture transfer function for single-phase flow of a compressible fluid in fractured porous media. In the proposed semi-analytical model, the pressure variability of viscosity and isothermal compressibility is considered by solving the nonlinear partial differential equation of compressible fluid flow in the fractured media. The closed form solution provided can be applied to flow of compressible fluids with variable matrix block size distribution in naturally fractured gas reservoirs.     

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).     

Up-Scaling of Double Porosity Fractured Media Using Continuous-Time Random Walks Methods

We present a new application of continuous time random walks (CTRW) methods to model fluid flows in fractured rocks. The proposed method allows large scale equivalent permeability tensors and matrix/fractures exchange function to be computed from high resolution maps of fractured porous media. Knowing these parameters allows us to carry out large scale simulations of flows governed by the dual porosity equations of Warren and Root. A direct connection between the exchange function and the time correlation function of the presence in the fractures of a particle undergoing a suitable Brownian motion over the whole medium is derived. This connection allows us to develop an efficient numerical method to compute the transient exchange term within the complete range of time scales of interest. It also gives an alternative probabilistic interpretation of the Warren and Root model. For the sake of simplicity, in the present paper, the method will only be developed to Cartesian structured grids, although it can be adapted for unstructured grids highly suited to describing complex fracture networks.     

Dual Porosity and Dual Permeability Modeling of Horizontal Well in Naturally Fractured Reservoir

This article is the first investigation on the dual permeability flow issue for horizontal well-production in a naturally fractured dual-porosity reservoir. Based on the inter-porosity flow from matrix system to fracture system and treating the media directly connected with horizontal wellbore as matrix and fracture systems, we established a model of horizontal well-production and then solved the model using some modern mathematical methods, such as Laplace integral transformation, separation of variables, eigenvalue, and eigenfunction. Later in the article, we obtained the standard log–log type curves using numerical simulation and analyzed the transient flow behavior thoroughly, which showed it is dual porosity and dual permeability flow behavior. The numerical simulation results showed that there are obvious differences between dual permeability and single permeability models. The dual permeability flow behavior accelerates energy supplement during production and reduces the classical matrix-fracture (V-shaped) response. We also showed that type curves characteristics are affected by external boundary conditions, the parameter κ, ω _f and λ _mf, etc. The research results show that our model would be a good semi-analytical model supplied to users. Because the single permeability modeling ignores the direct fluid supply from matrix to wellbore, we recommend using the dual permeability modeling to make well testing and rate decline interpretation in real case studies.     

Unified semi‐analytical wall boundary conditions for inviscid,laminar or turbulent flows in the meshless SPH method

Wall boundary conditions in smoothed particle hydrodynamics (SPH) is a key issue to perform accurate simulations. We propose here a new approach based on a renormalising factor for writing all boundary terms. This factor depends on the local shape of a wall and on the position of a particle relative to the wall, which is described by segments (in two‐dimensions), instead of the cumbersome fictitious or ghost particles used in most existing SPH models. By solving a dynamic equation for the renormalising factor, we significantly improve traditional wall treatment in SPH, for pressure forces, wall friction and turbulent conditions. The new model is demonstrated for cases including hydrostatic conditions for still water in a tank of complex geometry and a dam break over triangular bed profile with sharp angle where significant improved behaviour is obtained in comparison with the conventional boundary techniques. The latter case is also compared with a finite volume and volume‐of‐fluid scheme. The performance of the model for a two‐dimensional laminar flow in a channel is demonstrated where the profiles of velocity are in agreement with the theoretical ones, demonstrating that the derived wall shear stress balances the pressure gradient. Finally, the performance of the model is demonstrated for flow in a schematic fish pass where both the velocity field and turbulent viscosity fields are satisfactorily reproduced compared with mesh‐based codes. Copyright © 2012 John Wiley & Sons, Ltd.     

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.     

On the application of the Helmholtz–Hodge decomposition in projection methods for incompressible flows with general boundary conditions

This paper is concerned with the analysis of the Helmholtz–Hodge decomposition theorem since it plays a fundamental role in the projection methods that are adopted in the numerical solution of the Navier–Stokes equations for incompressible flows. The paper highlights the role of the orthogonal decomposition of a vector field in a bounded domain when general boundary conditions are in effect. In fact, even if Fractional Time‐Step Methods are standard procedures for de‐coupling the pressure gradient and the velocity field, many problems are encountered in performing the decoupling with higher accuracy. Since the problem of determining a unique and orthogonal decomposition requires only one boundary condition to be well posed, thus either the normal or the tangential ones, result exactly imposed at the end of the projection. Numerical errors are introduced in terms of both the pressure and the velocity but the orthogonality of decomposition guarantees that the former does not contribute to affect the accuracy of the latter. Moreover, it is shown that depending on the meaning of the vector to be decomposed, i.e. acceleration or velocity, the true orthogonal projector can be defined only when suitable boundary conditions are verified. Conversely, it is shown that when the decomposition results non‐orthogonal, the velocity accuracy suffers of other errors. The issue on the resulting accuracy order of the procedure is clearly addressed by means of several accuracy studies and a strategy for improving it is proposed. This paper follows and integrates the issues reported in Iannelli and Denaro (Int. J. Numer. Meth. Fluids 2003; 42 : 399–437). Copyright © 2003 John Wiley & Sons, Ltd.     

Development of an Analytical Time-Dependent Matrix/Fracture Shape Factor for Countercurrent Imbibition in Simulation of Fractured Reservoirs

In dual porosity modeling of naturally fractured reservoirs, fluids exchange between the high porous matrix blocks and high permeable fracture systems is governed by transfer function. Therefore, transfer function, and specially shape factor as the main part of it, control fluids flow behavior, which certainly have significant effects on development and management plan of naturally fractured reservoirs. Also several formulations have been proposed for shape factor by a number of researchers, nearly all of them derived for expansion mechanism. But, shape factor is a phase sensitive parameter that can greatly affect results of simulation. Moreover, several shortcomings are inherent in the derived expressions of shape factor for imbibition process. The main aim of this work is to develop a new time-dependent matrix–fracture shape factor specific to countercurrent imbibition. In this study, fluid saturation distribution within a matrix block is analytically derived by solving capillary–diffusion equation under different imposed boundary conditions for the process where countercurrent imbibition is the dominant oil drive mechanism. The validity of the solutions is checked against literature experimental data (Bourbiaux and Kalaydjian, SPERE 5, 361–368, SPE 18283, 1990) and also by performing single porosity fine grid simulations. Then, the concept of analogy between the transport phenomena is employed to propose a new expression for matrix–fracture transfer function that is used to derive transient shape factor. It is illustrated in this article that time variation of imbibtion rate and shape factor can be used to diagnose different states of imbibition process. Although, the displacement process and employed approaches are completely different in this and other studies (Chang, Technical report, 1993; Kazemi and Gilman (eds.) Flow and contaminant transport in fractured rock. Academic Press, San Dieg, 1993; Zimmerman et al., Water Resour Res, 29, 2127–2137, 1993; Lim and Aziz, J Pet Sci Eng 13, 169–178, 1995), but we arrived at the consistent values of shape factor under limiting condition of pseudo steady state flow. This means that after establishment of pseudo steady state, shape factor is only controlled by matrix geometry regardless of the displacement process, i.e., expansion or imbibition mechanism, However, shape factor is completely phase sensitive and process dependent during unsteady and late-transient states. Finally, boundary condition dependency of shape factor is investigated.     

High-precision analytic solution of the problem on bending vibrations of a clamped square plate

An original iteration algorithm is used to construct new analytic expressions for computing approximate natural frequencies and shape modes of bending vibrations of a square homogeneous plate clamped along its contour. The errors are estimated by comparing with the results of well-known numerical high-precision computations. The results of analytic computations are also compared with experimental data obtained by the author by the resonance method. The proposed research technique and the obtained high-precision expressions for the natural shape modes can be used in the case of rectangular plates and for other types of boundary conditions. A numerical-analytical method is used to show that the small isoperimetric theorem holds.     

Simulation of Pressure Transient Behavior for Asymmetrically Finite-Conductivity Fractured Wells in Coal Reservoirs

Based on Fick’s law in matrix and Darcy flow in cleats and hydraulic fractures, a new semi-analytical model considering the effects of boundary conditions was presented to investigate pressure transient behavior for asymmetrically fractured wells in coal reservoirs. The new model is more accurate than previous model proposed by Anbarci and Ertekin, SPE annual technical conference and exhibition, New Orleans, 27–30 Sept 1998 because new model is expressed in the form of integral expressions and is validated well through numerical simulation. (1) In this paper, the effects of parameters including fracture conductivity, coal reservoir porosity and permeability, fracture asymmetry factor, sorption time constant, fracture half-length, and coalbed methane (CBM) viscosity on bottomhole pressure behavior were discussed in detail. (2) Type curves were established to analyze both transient pressure behavior and flow characteristics in CBM reservoir. According to the characteristics of dimensionless pseudo pressure derivative curves, the process of the flow for fractured CBM wells was divided into six sub-stages. (3) This paper showed the comparison of transient steady state and pseudo steady state models. (4) The effects of parameters including transfer coefficient, wellbore storage coefficient, storage coefficient of cleat, fracture conductivity, fracture asymmetry factor, and rate coefficient on the shape of type curves were also discussed in detail, indicating that it is necessary to keep a bigger fracture conductivity and fracture symmetry for enhancing well production and reducing pressure depletion during the hydraulic fracturing design.     

Large eddy simulation in a turbulent channel flow using exit boundary conditions

The influence of the exit boundary conditions on the vanishing first derivative of the velocity components and constant pressure on the large eddy simulation of the fully developed turbulent channel flow has been investigated for equidistant and stretched grids at the channel exit. Results show that the chosen exit boundary conditions introduce some small disturbances that are mostly damped by the grid stretching. The difference of rms values between the fully developed turbulent channel flow with periodicity conditions and the fully developed channel flow using inlet and the exit boundary conditions is less than 10% for the equidistant grids and less than 5% for the stretched grids. The chosen boundary conditions are of interest because they may be used in complex problems with back flow at the exit. Copyright © 1999 John Wiley & Sons, Ltd.     

Analyzing the 2D fracture problems via the enriched boundary element-free method

In this paper, the enriched boundary element-free method for two-dimensional fracture problems is presented. An improved moving least-squares (IMLS) approximation, in which the orthogonal function system with a weight function is used as the basis function, is used to obtain the shape functions. The IMLS approximation has greater computational efficiency and precision than the existing moving least-squares (MLS) approximation, and does not lead to an ill-conditioned system of equations. Combining the boundary integral equation (BIE) method and the IMLS approximation, a boundary element-free method (BEFM), for two-dimensional fracture problems is obtained. For two-dimensional fracture problems, the enriched basis function is used at the tip of the crack, and then the enriched BEFM is presented. In comparison with other existing meshless boundary integral equation methods, the BEFM is a direct numerical method in which the basic unknown quantity is the real solution of the nodal variables, and the boundary conditions can be implemented easily, which leads to a greater computational precision. When the enriched BEFM is used, the singularity of the stresses at the tip of the crack can be shown better than that in the BEFM. For the purposes of demonstration, some selected numerical examples are solved using the enriched BEFM.     

Efficient computation of the flow around single fluid particles using an artificial boundary condition

Two‐phase flows around fluid particles are often considered to be in infinite domains, to avoid influence of the domain walls. Numerical simulations, however, must be modeled with a bounded domain, thus introducing artificial boundaries. Modeling of fluid flow in a domain with such artificial boundaries requires a careful choice of suitable boundary conditions. Slip boundary conditions for example can have a large impact on the computational results if the domain is chosen to be too small, because they model impermeable walls. This paper introduces an artificial boundary condition for simulations of the flow around single rising or settling fluid particles based on the approximated decay behavior of the velocity and the pressure field in the surrounding liquid. This is applied to the simulation of rising gas bubbles in systems with a Reynolds number of up to 50, and the outcome is compared with experimental results and simulations with slip boundary condition. It is found that domain size can be reduced by a factor of about two compared with slip boundary conditions without loss of accuracy. Copyright © 2014 John Wiley & Sons, Ltd.     

A family of multi‐point flux approximation schemes for general element types in two and three dimensions with convergence performance

A family of flux‐continuous, locally conservative, control‐volume‐distributed multi‐point flux approximation (CVD‐MPFA) schemes has been developed for solving the general geometry‐permeability tensor pressure equation on structured and unstructured grids. These schemes are applicable to the full‐tensor pressure equation with generally discontinuous coefficients and remove the O(1) errors introduced by standard reservoir simulation schemes when applied to full‐tensor flow approximation. The family of flux‐continuous schemes is characterized by a quadrature parameterization. Improved numerical convergence for the family of CVD‐MPFA schemes using the quadrature parameterization has been observed for structured and unstructured grids in two dimensions. The CVD‐MPFA family cell‐vertex formulation is extended to classical general element types in 3‐D including prisms, pyramids, hexahedra and tetrahedra. A numerical convergence study of the CVD‐MPFA schemes on general unstructured grids comprising of triangular elements in 2‐D and prismatic, pyramidal, hexahedral and tetrahedral shape elements in 3‐D is presented. Copyright © 2011 John Wiley & Sons, Ltd.     

Derivation of the Kazemi–Gilman–Elsharkawy Generalized Dual Porosity Shape Factor

Kazemi et al. (SPE Reserv Eng 7(2):219–227, 1992) suggested an empirical matrix-fracture transfer function, verified based on experimental data of Mattax and Kyte (Trans AIME 225(15):177–184, 1962), to model fluid flow in naturally fractured dual porosity petroleum reservoirs using a dual-porosity numerical simulator. Their generalized shape factor should be valid for all possible irregular matrix blocks. The factor is calculated based on the volume of the matrix block, the surface open to flow in all directions and the distances of these surfaces to the centre of the matrix block. The summation is done over all open surfaces of a matrix block. Kazemi et al. (1992) showed that for rectangles and cylinders the formula reduces to the well-known forms of the shape factor. By the time, many authors indicated the validity of the formula, but no theoretical proof was offered for that so far. This study derives the Kazemi et al. (1992) shape factor using control volume finite difference discretization on the fracture-matrix dual continuum. The matrix blocks are handled as Voronoi polyhedra. The derivation is given for both isotropic and tensorial matrix permeability. Based on this derivation the authors conclude that the Kazemi et al. (SPE Reserv Eng 7(2):219–227, 1992) formula is exact under pseudo-steady-state conditions within the dual continuum mathematical concept of natural fractured dual porosity systems.     

A perturbation method for mixed boundary-value problems in pressure transient testing

Data recorded during a well test is interpreted for formation parameters, such as permeability, by comparing the measured pressure transient with that predicted by a mathematical model of the system. In single-phase homogeneous situations, this model is based on the linear diffusion equation. Despite the simplicity of this equation, it is often difficult or laborious to construct exact solutions, in some cases because the relevant problem is of mixed boundary-value character. Here we describe an approximate calculational method which, with little effort, gives good results in a variety of well test problems of this type.We first prove a rather general perturbation theorem, and then apply it in three illustrative cases. The first example provides a test of our basic result for the case of a fully penetrating vertical fracture, for which an exact (though nonelementary) solution is known. The second example is that of partially penetrating or horizontal wells; it is shown that our general result can provide a mathematical basis for the so-called pressure-averaging technique in which the pressure-flux relationship is approximated by assuming that the flux is uniform along the well and then computing the spatial average of the wellbore pressure. Our third example is a new result for the steady-state pressure drop due to flow into a small circular hole on the surface of an impermeable cylinder. This example has relevance for the Schlumberger RFT tool, which withdraws fluid from the formation via a circular probe which penetrates the mudcake surrounding the borehole wall. Numerical results for the shape factor which represents the effect of the borehole curvature are provided.     

On the critical problem of F. D. pressure treatment for laminar flows confined by permeable walls

In the present work the viscous (low Reynolds) flow in plane ducts confined by permeable walls has been studied. A simple model of the filtrating walls has been used, with the normal velocity component proportional to the pressure jump across the wall, resulting in a non-standard boundary value Navier-Stokes problem. A critical analysis of the appropriate boundary condition and pressure problem has led to the conclusions of employing a simple explicit finite volume approach, and of avoiding the use of higher order finite difference schemes. In this paper a special emphasis on the structure of the involved computational matrices has been given to illustrate the chosen algorithm. The latter yields a steady state solution that is second order accurate in space, and it has an accuracy in time of order ≤ Δt (the time step), due to the explicit treatment of the velocity boundary conditions along the membrane. The model has been tested to study the effects of the inlet/outlet conditions, Reynolds number and filtrating wall constant.     

Buoyancy modelling with incompressible SPH for laminar and turbulent flows

This work aims to model buoyant, laminar or turbulent flows, using a two‐dimensional incompressible smoothed particle hydrodynamics model with accurate wall boundary conditions. The buoyancy effects are modelled through the Boussinesq approximation coupled to a heat equation, which makes it possible to apply an incompressible algorithm to compute the pressure field from a Poisson equation. Based on our previous work [1], we extend the unified semi‐analytical wall boundary conditions to the present model. The latter is also combined to a Reynolds‐averaged Navier–Stokes approach to treat turbulent flows. The k ? ? turbulence model is used, where buoyancy is modelled through an additional term in the k ? ? equations like in mesh‐based methods. We propose a unified framework to prescribe isothermal (Dirichlet) or to impose heat flux (Neumann) wall boundary conditions in incompressible smoothed particle hydrodynamics. To illustrate this, a theoretical case is presented (laminar heated Poiseuille flow), where excellent agreement with the theoretical solution is obtained. Several benchmark cases are then proposed: a lock‐exchange flow, two laminar and one turbulent flow in differentially heated cavities, and finally a turbulent heated Poiseuille flow. Comparisons are provided with a finite volume approach using an open‐source industrial code. Copyright © 2015 John Wiley & Sons, Ltd.