Well modeling is essential for various engineering applications, as for example reservoir simulation, geothermal energy production or energy storage, where injection or extraction processes strongly influence the flow behavior. Usually the well geometry is not explicitly resolved in the mesh but instead modeled as a line source with given extraction or injection rate. However, this simplified approach introduces singularities, meaning that the logarithmic solution profiles are undefined at the center-line of the well. This leads to a significant deviation between numerical and analytical solution in the near-well region. For a better approximation, locally refined meshes around the wells are needed, which however deteriorate efficiency and is therefore often not suitable for field-scale simulations, especially when multiple wells are present. Similar issues are encountered for the modeling of vascularized biological tissue perfusion (Koch2018b, ; cattaneo2014computational, ) or the modeling of root water update (Koch2018a, ; Doussan1998, ).
A common approach is the use of well-index-based well models. Such well models aim to find a relation between well rate, bottom hole pressure and numerically calculated pressure (well-block pressure) for each cell (grid-block) that contains the well. In reservoir engineering such a relation is denoted as well index. The first theoretical derivation of a well index for two-dimensional structured uniform grids with isotropic permeability has been presented by Peaceman (Peaceman1978, ). He has shown that the well-block pressure differs from the areal averaged analytical pressure and introduced a new relation by using an equivalent well radius. The equivalent well radius is defined as the distance (relative to the well location) at which the analytical and numerical pressures are equal. A generalization for structured two-dimensional non-square grids () and anisotropic but diagonal permeability tensors has been presented by the same author in (Peaceman1983, ). However, these well models are restricted to uniform grids where the well is oriented with one of the grid axes. Furthermore, since the effective well radius relates numerical and analytical pressure values, the well index has to be calculated depending on the used discretization scheme. Thus, Peaceman’s well model is only valid for a cell-centered finite difference scheme with 5-point stencil. Discussion of other discretization schemes can be found in (chen2009well, ), again with the restriction of two-dimensional grids. Enhancements include, among others, three-dimensional slanted wells (Alvestad1994, ; Aavatsmark2003, ), Green’s functions for the computation of well indices (Babu1989, ; Babu1991, ; Wolfsteiner2003, ), or the singularity subtraction method to obtain smooth solutions in the near-well region (hales1997improved, ; Gjerde2018, ).
In this work, a new approach for obtaining more accurate source term for a given well bottom hole pressure is presented. The new model is, in contrast to most of the existence methods, independent of the discretization scheme and can be used for general unstructured grids. Additionally, the presented method is not restricted to diagonal tensors and thus works for general anisotropic permeabilities. In Section 2, we derive a well model, initially for isotropic porous media, for which the fluid mass injected by a well is distributed over a small neighborhood around the well, using kernel functions. The derivation follows the idea recently presented in (Koch2019a, ), however we herein discuss the case without membrane or casing. The model yields a pressure solution without singularity, from which the source term can be reconstructed using a relation found with the analytical solution for the case of an infinite well in an infinite medium. The model generalizes to more complex problems due to the superposition principal valid for the Laplace operator, a linear operator (Koch2019a, ). In Section 3, the model is generalized to porous media with general anisotropic permeabilities, based on an analytical solution constructed in Section 3.1 using a series of coordinate transformations. We shown that the general model reduces to the model derived in Section 2 for isotropic permeabilities. Finally, the new well model is analyzed with several numerical experiments in Section 5. The results indicate that the model is consistent for different anisotropy ratios, robust with respect to rotations of the well relative to the computational grid, and to rotations of an anisotropic permeability tensor. A comparison with a Peaceman-type well model in a setup with a K-orthogonal grid and an embedded slanted well suggests that the new model more accurately approximates the fluid exchange between well and rock matrix.
2 A well model with distributed source for isotropic media
First, we derive a well model with distributed source for porous media with isotropic permeability tensor, not including a well casing. The derivation follows (Koch2019a, ), where a casing is included in form of a membrane in a different but related application, that is modeling fluid exchange between the vascular system and the embedding biological tissue. Stationary single-phase flow around a well with radius , in an isotropic porous medium with permeability , can be described by the following flow equation
where is the fluid pressure, the fluid density and its dynamic viscosity. Denoted by is a set of kernel functions that distribute () over a small tubular support region, , with radius , around a well segment , such that outside the support region. We choose kernel functions with the property
where , , are the radial, angular, and axial coordinate in a segment-local cylinder coordinate system, and is the length of segment .
Assume a radially symmetric zone (distance from center-line) around the well, with denoting the pressure at distance from the well center-line, and constant fluid density and viscosity. Then, the pressure for is described by the analytical solution
The constant is determined by fixing a well pressure, ,
Consequently, the source term can be expressed in terms of and as
We choose a simple kernel function which regularizes the pressure solution for ,
where . The pressure for can be obtained by integration from Eq. 2.1, yielding
where is the so called flux scaling factor. The flux scaling factor can be also expressed independent of the pressure. To this end, Eq. 2.7 is evaluated at , so that is expressed in terms of ,
3 A well model with distributed source for anisotropic media
In the following section, the developed well model is extended for porous media with anisotropic permeability. In Section 3.1, we derive an analytical solution for one-phase flow around an infinitely long cylindrical well embedded in an infinite porous domain in . This derivation motivates the choice of a suitable kernel function for anisotropic problems, presented in Section 3.3.
3.1 Analytical solution for anisotropic permeability and slanted well
In the following section, we derive an analytical solution for one-phase flow around an infinite cylindrical well with radius in an infinite porous domain with anisotropic, homogeneous permeability. We assume, without loss of generality, that the well axis passes through the origin of the Cartesian coordinate system, and denote by
a unit vector parallel to the well signifying the well orientation. We seek an analytical expression for the hydraulic pressuresuch that
for a constant well pressure in and some specific pumping rate in given on . The total mass flow over the boundary of a well segment of length is thus given by .
From thermodynamic constraints, is a positive definite and symmetric, second-order tensor field. Hence, can be decomposed such that
is a diagonal matrix composed of the eigenvaluesof ,
is a rotation matrix with the corresponding eigenvectors as columns, anddenotes the transpose of a matrix . Further useful properties derived from the decomposition are , where denotes the determinant of , and , where , .
It is well known that the anisotropic one-phase flow problem can be transformed to an isotropic problem using a coordinate transformation (Aavatsmark2003, ; Fitts2006, ; Aavatsmark2016, ; Peaceman1983, ; Bear1965, )
with the stretching matrix , where is an arbitrary scalar constant, that we choose as (cf. (Aavatsmark2003, )), rendering the transformation isochoric. The transformation deforms the well cylinder such that a cross-section orthogonal to the transformed well direction is elliptical. The solution to the isotropic problem
is identical on two parallel planes perpendicular to the transformed (normalized) well direction, . This motivates the rotation of the coordinate system such that the first and second axis are aligned with the major and minor axis of the well-bore ellipse and third axis is aligned with . To determine the corresponding rotation matrix , we need to characterize this well-bore ellipse. The well cylinder in -coordinates is given by
After stretching, the coordinate system can be rotated with the rotation matrix so that the third axis is aligned with the well direction. Then, projecting into the plane perpendicular to the well direction yields the well-bore ellipse equation
in -coordinates, where
The rotation matrix can be derived using Rodrigues’ rotation formula as shown in Appendix A. The length of the major and minor ellipse axis are found as and , where are the eigenvalues of , and the axis orientations are given by , , where denote the corresponding eigenvectors of . We assume that the eigenvalues and eigenvectors are sorted such that , and oriented such that . Finally the desired rotation is given by
is rotating about the well direction axis such that the coordinate system is aligned with the principal ellipse axes.
Following the derivations from above, we now have to solve a two-dimensional isotropic Laplace problem with boundary conditions prescribed on an ellipse. To this end, we note that the transformation of a harmonic function (a function satisfying Laplace’s equation ) with a conformal (angle-preserving) mapping yields another harmonic function (Nehari1975, ) (see Appendix B). Using a Joukowsky transformation, a conformal mapping well-known from aerodynamics (joukowsky1910, ), the isotropic problem with a well with elliptic cross-sections, can be transformed to an isotropic problem with circular cross-sections (Fitts2006, ). Transforming into the complex plane (parametrizing the well-bore ellipse plane)
the (inverse) Joukowsky transformation
transforms elliptic isobars into circular isobars, where and , , are the major and minor axis of the well-bore ellipse, as derived above. In particular, the well-bore ellipse (where ) is mapped onto a circle with radius . Finally, in the new coordinate system we find the (now) radially symmetric analytical solution to problem Eq. 3.1
where the source scaling factor is necessary to recover the original source on . This can be derived from simple geometric considerations as shown in Appendix C. The corresponding to some in original coordinates is obtained by using all above-mentioned transformations after each other as follows
Such a solution for a slanted well ( with respect to vertical axis) and anisotropic permeability tensor
is visualized in Fig. 2.
3.2 Properties of the conformal mapping
To construct a suitable kernel function for anisotropic problems, we first have a closer look at the properties of the employed Joukowsky transformation. The effect of the mapping , Eq. 3.11, is shown in Fig. 3. Points on the exterior of a line on the real axis between and are mapped onto the exterior of a circle with radius . The ellipse with major axis and minor axis is mapped onto a circle with radius . Going further away from the well, the deformation due to the mapping is less and less pronounced. This matches the expectations for the physical flow problem, since isobars at large distance from an elliptical well-bore become increasingly circular in isotropic media.
The inverse transformation is given by
where the restriction on is necessary to obtain a one-to-one mapping. The transformation can equally be interpreted as an mapping. The Jacobian of the transformation for and has the form
which follows from the the Cauchy–Riemann equations (Rudin1987, ). Since the transformation can be viewed as the composition of a scaling and a rotation, it is angle-preserving. The transformation is associated with a spatially dependent volume deformation characterized by the determinant of . Furthermore, it can be shown that the Laplace operator behaves as follows under the transformation ,
by computing the derivative of the real and the imaginary part of
separately, applying the chain rule and the Cauchy–Riemann equations, as shown for completeness inAppendix B. From Eq. 3.17 follows that
where is real part of , and the absolute value of . We note that quickly converges to the value with increasing , that is with increasing distance from the well. The function is plotted in Fig. 4 exemplarily for .
3.3 A kernel function for anisotropic media
Instead of excluding the well domain from and modeling infiltration or extraction by a flux boundary condition, we will now model the action of the well on the flow field by a spatially distributed source term, as presented for the isotropic problem,
From the above derivations, we know that solving Eq. 3.20 in -coordinates is straight-forward. Hence, we choose kernel functions in -coordinates and then transform to -coordinates so that the pressure solution satisfies Eq. 3.20. Motivated by the properties of the Joukowsky transform (see Fig. 3), we choose a local kernel that is constant on the annulus with inner radius and outer radius ,
In -coordinates, we can find a solution to the problem
for a given constant well pressure , and constant density and viscosity.
By means of integration (cf. (Koch2019a, )), we get
where is the fluid pressure evaluated on the well center-line. Note that for and , the isotropic solution with for a circular constant kernel (Eqs. 2.10 and 2.7) is obtained. From the transformation of the Laplace operator, Eq. 3.18, we see that the problem
with altered kernel function is equivalent to Eq. 3.22.
The transformation changes shape of the kernel support from an annulus to an ellipse . Inverting extrudes the solution along the well center-line, and inverting the rotation and stretch described by and results in a kernel support region in the shape of an elliptic cylinder. Moreover, each ellipse with normal vector is transformed to an ellipse , that is the intersection the elliptic cylinder with a plane with the normal vector , centered at on the well center-line. The transformation and the normal vector are visualized in Fig. 5. We note that if none of the principal axes of the permeability tensor are aligned with the well direction, is not parallel to the well direction in -coordinates. The integral of the right-hand side of Eq. 3.20 for a well segment with length is equal to the integral over the kernel support which has the shape of the elliptic cylinder given by
Using , and exploiting that was chosen such that , it can be shown that
where is a local coordinate along the transformed well direction, and the last equality is proven in Appendix C. This is the desired property of the kernel function for the anisotropic case corresponding to Eq. 2.2 for the isotropic case.
4 Numerical method
We discretize Eq. 3.20 using a cell-centered finite volume method with multi-point flux approximation (MPFA) (Aavatsmark2002, ). The domain is decomposed into control volumes such that the computational mesh is a discrete representation of . Furthermore, each control volume boundary, , can be split into a finite number of faces , such that , with denoting a neighboring control volume. Integrating Eq. 3.20 over a control volume and applying the Gauss divergence theorem on the left hand side yields
where is the unit outward-pointing normal on face . The exact fluxes are approximated by numerical fluxes
which are computed using the MPFA-O method described in (Aavatsmark2002, ). The discrete source term is computed as
where is a numerical approximation of the source term integral over the intersection ,
We note that due to the dependency of on the proposed method is non-local in the sense that non-neighbor cells (where
is the empty set or a single point) may have an associated degree of freedom that depends on the degree of freedom of.
4.1 Kernel integration
The kernel integral in Eq. 4.3
is not easily approximated with a quadrature rule, since the intersection , that is the intersection of an elliptic cylinder with for example a hexahedron is difficult to compute. However, we use the same idea as in (Koch2019a, ), and remark that the integral over the entire support is known exactly; see Eq. 3.27. Hence, the integration problem can be reformulated as the distribution of the known integral over all intersected control volumes weighted with the respective support volume fractions. Following (Koch2019a, ), we create integration points with known volume elements of similar size and shape, so that
Computing the weights for each cell is a pre-processing step that only has to be done once for each computational mesh and well geometry.
5 Numerical experiments and discussion
We present numerical experiments using the presented method in different setups. All experiments are conducted with constant fluid density and viscosity . The well pressure is constant, , and the well radius is if not specified otherwise. The permeability tensor is given as
are rotation matrices rotating vectors about , by the rotation angle , , respectively, and is a given dimensionless K-anisotropy ratio . The domain is split in two regions, and . The well center-line is given by the line through the origin and , where , are given in Eq. 5.2 and , , are rotation angles. The analytical solution for all cases is given in Eq. 3.23, , and (in ). For all setups the inner kernel radius is chosen as . In all of and on the boundary the analytical solution is enforced by Dirichlet constraints, modeling the infinite well. The computational mesh is a structured grid composed of regular hexahedra . Furthermore, we define two error measures.
is the relative discrete -norm of the pressure, where is the exact pressure evaluated at the cell centroid and the discrete numerical cell pressure, and
is the relative discrete -norm of the source term, where is the length of the intersection of cell and the well center-line , is the discrete source term given in Eq. 4.4. All setups are implemented in DuMux (flemisch2011dumu, ), an open-source porous media simulator based on Dune (dunegridpaperI:08, ; dunegridpaperII:08, ).
5.1 Grid convergence for different anisotropy ratios
In the first numerical experiment grid convergence is investigated for different anisotropy ratios . To this end, , where is defined as the maximum distance between two vertices of the cell . Starting at a grid resolution for of cells (), the grid is refined uniformly. Figure 7 shows the errors and for different grid resolutions and values of , for and , so that is a full tensor and none of the principal axis of is aligned with the well direction. For all , the method shows second order convergence for the pressure in the given norm, as expected for the MPFA-O method (Schneider2018, ) (super convergence at cell centers). The source term is a linear function of the pressure and also exhibits second order convergence. We note that the errors for different are not directly comparable since the analytical solution for changes with , although is constant.
However, the convergence rates are shown to be independent of with increasing grid resolution. The convergence rates for (slope of the lines in Fig. 7) are presented in Table 1. It can be seen that rates for large grid cells and large are slightly smaller. This is because the kernel support is still under-resolved by the computational grid. For example, for , the kernel ellipse in -coordinates has major and minor axis of , , respectively, while for the lowest grid resolution.
5.2 Influence of the outer kernel radius
In (Koch2019a, ), it is suggested that increasing the kernel support region (increasing ), has a similar effect on as refining the grid. However, the pressure solution is then regularized in a larger region, so that there is a trade-off between the accuracy of the source term and the accuracy of the pressure field with respect to the unmodified problem (). However, every discrete cell can be also interpreted as a kernel support region, such that the choice of enables us to better control the discretization error as soon as becomes larger than .
As shown in (Peaceman1983, ) and Fig. 4, isobars become circular, in the transformed domain , with increasing distance to the well. Therefore, a reasonable simplification is if . This is completely analogous to the assumption of circular isobars in (Peaceman1983, )
, where an estimate of the error introduced by the assumptions is given for the two-dimensional case.
In the following numerical experiment, we step-wise increase the kernel radius , for the same grid. This is done once for the case, where and for the case for which is only slightly larger than . Furthermore, and . The results are shown in Fig. 8. First, it can be seen that doubling leads to a -times smaller error . This can be explained by the fact that the larger the kernel, the more grid cells resolve the kernel support, and the better is the approximation of . Furthermore, the result is consistent with the results in (Koch2019a, ). Moreover, Fig. 8 suggests that for the simplification of the kernel function () is not visible in , while for kernel radii slightly smaller than the well radius, the simplification increases by an order of magnitude in comparison to the case using the exact kernel function as derived in Section 3.3. The results show that the presented method is also applicable in cases where the grid resolution is very close to the well radius. An adaption of the presented method for other applications, such as the simulation of flow in vascularized tissue, where such ratios of vessel radius to cell size are typical, cf. (Koch2019a, ), is therefore well-conceivable.
5.3 Robustness with respect to rotation
In the following numerical experiment, we use a single computational mesh with a given resolution for : . First, the well direction is fixed, and the permeability tensor is rotated by varying and . Then the permeability tensor is fixed and the well is rotated by varying and . The results are shown in Fig. 9.
It can be seen that the presented well model is rather robust with respect to rotations. Possible effects influencing the approximation error , include the different quality of the kernel integral for different angles with respect to the grid axes, and differences in the flux approximation quality of the MPFA-O method depending on the face co-normal . Additionally, for different well angles the number and size of intersections can have an influence on the discrete error.
5.4 Comparison with a Peaceman-type well model
In particular for petroleum engineering applications, commercial codes typically use Peaceman-type well models (eclipse2014, ; imex2014, ). In (Peaceman1983, ), Peaceman extended his well known well-index-based well model for anisotropic diagonal permeability tensors and non-cubic but structured rectangular grids. The discrete source term in a computational cell is approximated by
where and are the horizontal dimensions of the cell containing the well, the length of the well segment contained in , and the Euler–Mascheroni constant. The Peaceman model has several known limitations. Its derivation only applies to K-orthogonal structured grids, where the well is oriented along one of the grid axes, and perfectly horizontally centered within a vertical column of computational cells . Furthermore, the derivation is specific to cell-centered finite difference schemes with 5-point stencil. Moreover, computational cells may have to be significantly larger than the well radius (depending on the degree of anisotropy) for optimal accuracy. The Peaceman model has been generalized for slanted wells with arbitrary orientation, for example in (Alvestad1994, ). The Alvestad model (Alvestad1994, ) has been adapted for finite volumes, for example in (Aavatsmark2003, ) (formula given in Appendix E, subsequently referred to as pm well model). Such extensions usually constitute a reasonable directional weighting of the original Peaceman model but are not directly derived from the mathematical analysis of the underlying problem (Aavatsmark2003, ).
The herein presented model has none of the above-mentioned limitations. In particular, the presented model is valid for arbitrary positive definite and symmetric permeability tensors, unstructured grids, and is independent of the discretization scheme. Moreover, the presented model is consistent and we show grid convergence in the numerical experiments in Section 5.1. However admittedly, the Peaceman-type models are cell-local, thus computationally cheaper and easier to implement.
Several limitations of the Peaceman well models make it difficult to fairly compare it with our new model. For cases for which all assumptions of Peaceman are valid, our numerical studies (not shown here) suggest that the Peaceman well model is generally superior to the presented model with distributed sources. This is because it takes the analytical solution as well as the spatial discretization method into account. For cases where some assumptions are violated, for example off-center wells or slanted wells, it is difficult to construct cases where the analytical solution is readily constructed but does not feature a singularity on the boundary. Our preliminary numerical studies for such cases (for example the slanted well case in Section 5.1 without rotation of the permeability tensor) show large deviations ( error in total source term) from the analytical solution for the pm well model. However, these errors may be distorted by errors made in the discrete approximation of the singular boundary condition, where the well intersects the boundary. Finally, for the general case of unstructured grids, simplex grids, and full permeability tenors it is unclear how to apply the original Peaceman model. However, we know that the presented method is consistent (at least for a single straight well), and thus, the numerical solution converges to the exact solution with grid refinement. Therefore, we expect that the numerical solution on a very fine grid using the distributed source model is a reasonable reference solution.
We compare our model to the pm well model in a numerical experiment. The computational domain contains a slanted straight well with end points at , . The permeability tensor is a diagonal tensor , with . The structured cube grid is successively, uniformly refined starting with cells (). The well radius is ( for the coarsest grid). The kernel support region (chosen as ) only extends over few cells in the coarsest grid, so that the regularization effect is minimized. On the boundary , we specify Neumann no-flow boundary conditions, that is , except for the planes perpendicular to the -axis, where the Dirichlet boundary condition , are enforced. The reference solution is computed with cells (). The computational domain with pressure iso-surfaces of the reference solution are shown in Fig. 10.
In Fig. 11, the relative integral source error
with respect to the reference solution is shown for grids with different refinement. In a variant of the distributed source model (ds), the extent of the kernel support is adapted to the grid size. This is to keep the regularization effect of the kernel function minimal in order to get, in addition to a good approximation of the source term, a better approximation of the pressure solution close to the well. For , the extent of the kernel ellipse is given by its major and minor axes, and . For the reference solution this extent is kept constant with grid refinement. While this ensures a very good approximation of the source term, the pressure solution is regularized in a larger neighborhood of the well. In the variant, the kernel support is adapted proportional to , so that for the finest grid shown in Fig. 11 ( cells, ), the major and minor axes measure and .
It is evident that the numerical solution for the distributed source model converges to the reference solution. More importantly, the relative error is small () even for the coarsest grid. In comparison, the difference to the Peaceman-type model is large (). In particular, the error grows with grid refinement (to ), signifying that the generalization of Peaceman’s model for arbitrarily-oriented wells is not consistent. The result is comparable with the observations in (Aavatsmark2003, , Table 2), where Alvestad’s well indices are compared to a new numerically computed well index, and it is shown that the difference between those two well indices grow, the higher the ratio. In the variant of the ds model, the error in the source term with respect to the reference solution also grows with larger ratio. However, the error is consistently smaller (by a factor ) than for the pm well model.
Figure 12 shows the numerical pressure solutions along the and the -axis, for the reference grid resolution. It can be clearly seen that for the reference solution (ds) the pressure solution is regularized. For the variant of ds, the regularization is minimized, however in the far field the solution matches the reference solution better than the pm well model, which is due to the better approximation of the source term (see Fig. 11). We also note that the regularized solution leads to an altered solution in the near-field of the well but to a better approximation of the source term and thus the far field pressure (outside the kernel support), whereas the poor approximation of the source term in the pm method leads to a globally poor pressure solution.
A new well model was presented for which the mass exchange between a well and an embedding porous medium is modeled with a source term spatially distributed by a local kernel function. In the spirit of well-index-based well models the source term for a well with given bottom hole pressure is computed based on the numerical pressure in cells intersecting the well. However, the presented derivation of the new model is independent of the discretization method and the type of computational grid. The new model was shown to be consistent in a numerical experiment and exhibited grid convergence with the expected rates. In the same experiment it is shown that the absolute error with respect to an analytical solution is relatively small, even for coarse computational grids and small kernel support. It was shown, that the error in the source term can be decreased by increasing the region over which the source term is distributed. However, coincidentally, the pressure profile close to the well (inside the kernel support) becomes increasingly regularized. A comparison with a Peaceman-type well model generalized for arbitrarily oriented wells, suggested that even if the region is chosen to be very small (only covering the neighboring cells of cells with well intersection), thus minimizing the regularization effect, the source term can be approximated with good accuracy ( error with respect to a reference solution), whereas the Peaceman-type well model for the same case showed larger differences (