A new and consistent well model for one-phase flow in anisotropic porous media using a distributed source model
Timo Koch Rainer Helmig Martin Schneider
cs.CE 76S05



1]Timo Kochcor1
timo.koch@iws.uni-stuttgart.de
1]Rainer Helmig
1]Martin Schneider

[cor1]Corresponding author
[1]Department of Hydromechanics and Modelling of Hydrosystems, University of Stuttgart, Pfaffenwaldring 61, 70569 Stuttgart, Germany


A new well model for one-phase flow in anisotropic porous media is introduced,
where the mass exchange between well and a porous medium is
modeled by spatially distributed source terms over a small neighborhood region. To this end,
we first present a compact derivation of the exact analytical solution for an arbitrarily
oriented, infinite well cylinder in an infinite porous medium with anisotropic permeability tensor in ℝ^3, for constant
well pressure and a given injection rate, using a conformal map. The analytical solution motivates the choice of a
kernel function to distribute the sources. The presented model is independent
from the discretization method and the choice of computational grids.
In numerical experiments, the new well model is shown to be consistent and robust
with respect to rotation of the well axis, rotation of the permeability tensor, and different anisotropy ratios. Finally, a comparison with
a Peaceman-type well model suggests that the new scheme leads to an increased accuracy
for injection (and production) rates for arbitrarily-oriented pressure-controlled wells.


 well model 1d-3d mixed-dimension anisotropic analytic solution Peaceman 

A new and consistent well model for one-phase flow in anisotropic porous media using a distributed source model
    [
    Received; accepted; 
===============================================================================================================





§ INTRODUCTION

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 <cit.> or the modeling of root water update <cit.>.

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 <cit.>. 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 (Δ x ≠Δ y) and anisotropic but diagonal permeability tensors has been presented by the same author in <cit.>. 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 <cit.>, again with the restriction of two-dimensional grids. Enhancements include, among others, three-dimensional slanted wells <cit.>, Green's functions for the computation of well indices <cit.>, or the singularity subtraction method to obtain smooth solutions in the near-well region <cit.>.

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 <ref>, 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 <cit.>, 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 <cit.>.
In <ref>, the model is generalized to porous media with general anisotropic permeabilities, based on an analytical solution constructed
in <ref> using a series of coordinate transformations. We shown that the general model reduces to the
model derived in <ref> for isotropic permeabilities.
Finally, the new well model is analyzed with several numerical experiments in <ref>. 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.



§ 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 <cit.>, 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 r_ω, in an isotropic porous medium with permeability k,
can be described by the following flow equation

    - ÷( ρ/μ k  p ) = q Φ_Λ  in Ω,

where p is the fluid pressure, ρ the fluid density and μ its dynamic viscosity.
Denoted by Φ_Λ is a set of kernel functions Φ_Λ_i
that distribute q () over a small tubular support region, 𝒮(Φ_Λ_i),
with radius ϱ(s), around a well segment i, such that Φ_Λ_i = 0 outside the support region.
We choose kernel functions Φ_Λ_i(s) with the property

    ∫_0^L_i∫_0^2π∫_0^ϱ(s)Φ_Λ_i r drdθds = L_i,

where r, θ, s are the radial, angular, and axial coordinate in a segment-local cylinder coordinate system, and L_i is the length of segment i.

Assume a radially symmetric zone (distance δ > r_ω from center-line) around the well, with p_δ denoting
the pressure at distance δ from the well center-line, and constant fluid density and viscosity.
Then, the pressure for r_ω < r < δ is described by the analytical solution

    p(r) = - μ/kρq/2πln r + C.

The constant C is determined by fixing a well pressure, p_ω,

    p_ω = p(r_ω) = -μ/kρq/2πln r_ω + C   ⇒ C = p_ω + μ/kρq/2πln r_ω,

so that

    p(r) = - μ/kρq/2πln( r/r_ω) + p_ω.

Consequently, the source term can be expressed in terms of p_ω and p_δ as

    q = 2π r_ωρ k/μ(p_ω - p_δ)/r_ωln( δ/r_ω)

We choose a simple kernel function which regularizes the pressure solution for r ≤ϱ,

    Φ^const(r) = 1/πϱ^2    r ≤ϱ,
    
        0     r > ϱ,

where ϱ≤δ.
The pressure for r < ρ can be obtained by integration from <ref>, yielding

    p(r) = 
      -μ/kρq/2π[ r^2/2ϱ^2 + ln(ϱ/r_ω) - 1/2] + p_ω    r ≤ϱ,
    
      -μ/kρq/2πln(r/r_ω) + p_ω    r > ϱ.


<Ref> graphically explains the most important symbols introduced in this section.
As the regularized pressure can be evaluated at the well center-line we can reformulate <ref>,

    q = 2πρ k/μ(p_ω - p_0)Ξ,   with Ξ = (p_ω - p_δ)/(p_ω - p_0)1/ln( δ/r_ω)

where Ξ is the so called flux scaling factor.
The flux scaling factor can be also expressed independent of the pressure. To this end,
<ref> is evaluated at r=0, so that p_ω is expressed in terms of p_0,

    p_0 = - (p_ω - p_0)Ξ[ ln(ϱ/r_ω) - 1/2] + p_ω,

where q was replaced by inserting <ref>.
It directly follows from <ref> that

    Ξ = [ln(ϱ/r_ω)-1/2]^-1.




§ 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 <ref>, we derive an analytical solution for one-phase flow around an infinitely long cylindrical well embedded
in an infinite porous domain in ℝ^3. This derivation motivates the choice of a suitable kernel function for
anisotropic problems, presented in <ref>.



 §.§ 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 r_ω in an infinite porous domain Ω̂ = ℝ^3∖Γ 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 pressure p such that

    - ÷( ρ/μ K  p ) = 0   in Ω̂,

for a constant well pressure p_ω in 
and some specific pumping rate q in  given on ∂Γ. The total mass flow over the
boundary of a well segment of length L is thus given by Q = qL.

From thermodynamic constraints, K is a positive definite and symmetric, second-order tensor field.
Hence, K can be decomposed such that

    K = Q D Q^T,

where D = diag(λ_1, λ_2, λ_3) is a diagonal matrix composed of the eigenvalues λ_i of K,
Q = [ ν⃗_K,1|ν⃗_K,2|ν⃗_K,3 ] is a rotation matrix with the corresponding eigenvectors as columns,
and A^T denotes the transpose of a matrix A.
Further useful properties derived from the decomposition are det(K) = λ_1λ_2λ_3,
where det(A) denotes the determinant of A, and K^n = Q Λ^n Q^T,
where D^r = diag(λ_1^r, λ_2^r, λ_3^r), r ∈ℝ.

It is well known that the anisotropic one-phase flow problem can be transformed to an isotropic problem
using a coordinate transformation <cit.>

    U: ℝ^3→ℝ^3, x⃗↦u⃗ = S̃x⃗,

with the stretching matrix S̃ = k_I^1/2 K^-1/2, where k_I is an arbitrary scalar constant, that we
choose as k_I = det(K)^-1/3 (cf. <cit.>), rendering the transformation isochoric.
The transformation u⃗ = S̃x⃗ deforms the well cylinder such that a cross-section
orthogonal to the transformed well direction is elliptical.
The solution to the isotropic problem

    - u( ρ/μ k_I _u p ) = 0   in Ω̂_u = U(Ω̂),

is identical on two parallel planes perpendicular to the transformed (normalized) well direction,
ψ⃗' = S̃ψ⃗||S̃ψ⃗||^-1.
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 R̃, we need to characterize this well-bore ellipse.
The well cylinder in x-coordinates is given by

    x⃗^T Ψx⃗ = r_ω^2,   Ψ = I - ψ⃗ψ⃗^T.

After stretching, the coordinate system can be rotated with the rotation matrix R 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

    v̂⃗̂^T E v̂⃗̂ = v̂⃗̂^T P^T R S̃^-1ΨS̃^-1 R^T P v̂⃗̂ = r_ω^2,   S̃^-1R^T Pv̂⃗̂ = x⃗

in v̂-coordinates, where

    R = 2(e⃗_3 + ψ⃗')(e⃗_3 + ψ⃗')^T/(e⃗_3 + ψ⃗')^T(e⃗_3 + ψ⃗') - I,   e⃗_3 = [[ 0; 0; 1 ]],    P = [[ 1 0; 0 1; 0 0 ]].

The rotation matrix R can be derived using Rodrigues' rotation formula as shown in <ref>.
The length of the major and minor ellipse axis are found as a = r_w γ_1^-1/2 and b = r_w γ_2^-1/2,
where γ_i are the eigenvalues of E, and the axis orientations are given by ν⃗_1 = Pν̂⃗̂_E,1, ν⃗_2 = Pν̂⃗̂_E,2,
where ν̂⃗̂_E,i denote the corresponding eigenvectors of E. We assume that the eigenvalues and eigenvectors are sorted such
that a≥ b, and oriented such that ψ⃗'⃗ = ν⃗_1 ×ν⃗_2. Finally the desired rotation is given by

    V: ℝ^3→ℝ^3, u⃗↦v⃗ = R̃u⃗ = R̂^T R^Tu⃗,

where

    R̂ = [ ν⃗_1 |ν⃗_2 |ψ⃗'⃗]

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 f
(a function satisfying Laplace's equation Δ f = ÷ f = 0)
with a conformal (angle-preserving) mapping yields another harmonic function <cit.> (see <ref>).
Using a Joukowsky transformation, a conformal mapping well-known from aerodynamics <cit.>,
the isotropic problem with a well with elliptic cross-sections,
can be transformed to an isotropic problem with circular cross-sections <cit.>.
Transforming into the complex plane (parametrizing the well-bore ellipse plane)

    Z: ℝ^3→ℂ, v⃗↦ z = Z̃v⃗ = [1, i, 0]v⃗ = v_1 + iv_2,

the (inverse) Joukowsky transformation

    T: ℂ→ℂ, z ↦ w = z + √(z-f)√(z+f),    f = √(a^2 - b^2)

transforms elliptic isobars into circular isobars, where a and b,
a≥ b, are the major and minor
axis of the well-bore ellipse, as derived above. In particular, the well-bore ellipse (where p = p_ω)
is mapped onto a circle with radius r_∘ = a + b. Finally, in the new coordinate system we find the (now) radially symmetric
analytical solution to problem <ref>

    p(w) = p_ω - μ/ρ k_Iq̂/2πln( | w |/r_∘)ζ,   q̂ = qζ = qab/r_ω^2,

where the source scaling factor ζ is necessary to recover the original source q on ∂Γ. This can be derived from
simple geometric considerations as shown in <ref>.
The w corresponding to some x ∈Ω̂ in original coordinates is obtained
by using all above-mentioned transformations after each other as follows

    w = T(Z(V(U(x⃗)))) = T(Z̃R̃S̃x⃗).

Such a solution for a slanted well (30 with respect to vertical axis) and anisotropic permeability tensor

    K_A = [[ 1 0 0; 0 5 4; 0 4 5 ]] 1e-10□

is visualized in <ref>.





 §.§ 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 T, <ref>, is shown in <ref>.
Points on the exterior of a line on the real axis between f and -f are mapped onto the exterior of a circle with radius f.
The ellipse with major axis a and minor axis b is mapped onto a circle with radius r_∘ = a + b.
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

    T^-1: ℂ→ℂ, w ↦ z = 1/2( w + f^2/w),   | w | > f,

where the restriction on | w | is necessary to obtain a one-to-one mapping. The transformation T^-1 can equally
be interpreted as an ℝ^2→ℝ^2 mapping.
The Jacobian of the transformation T^-1 for z = x + iy and w = u + iv has the form

    J_T^-1 = [[ ∂ x/∂ u ∂ x/∂ u; ∂ y/∂ v ∂ y/∂ v ]]
               = [[  η  ϵ; -ϵ  η ]],

which follows from the the Cauchy–Riemann equations <cit.>. 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 J_T^-1. Furthermore, it can be
shown that the Laplace operator behaves as follows under the transformation z = T^-1(w),

    Δ_w p = ∂^2 p/∂ u^2 + ∂^2 p/∂ v^2 = | ∂ T^-1/∂ w|^2 Δ_z p = |det(J_T^-1)| Δ_z p

by computing the derivative of the real and the imaginary part of p separately, applying the chain rule and the Cauchy–Riemann equations, as
shown for completeness in <ref>.
From <ref> follows that

    |det(J_T^-1)|^-1Δ_w p = Δ_z p,

for the transformation w = T(z). The determinant can be explicitly computed,
using <ref> and complex differentiation (shown in <ref>) as

    |det(J_T^-1)(w)| =  | ∂ T^-1/∂ w|^2 = | ∂ z/∂ u|^2 = 1/4( 1 + f^4 - 2f^2(w^2)/| w |^4) := Φ_J^-1,

where (w^2) is real part of w^2, and | w | the absolute value of w. We note that Φ_J quickly converges
to the value 4 with increasing | w |, that is with increasing distance from the well. The function Φ_J^-1 is
plotted in <ref> exemplarily for f = 1.2.




 §.§ A kernel function for anisotropic media


Instead of excluding the well domain Γ from Ω = ℝ^3 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,

    - ÷( ρ/μ K  p ) = qζΦ_Λ  in Ω.

From the above derivations, we know that solving <ref> in w-coordinates is straight-forward. Hence,
we choose kernel functions in w-coordinates and then transform to x-coordinates so that the pressure solution
satisfies <ref>. Motivated by the properties of
the Joukowsky transform (see <ref>),
we choose a local kernel that is constant on the annulus with inner radius f < ϱ_i ≤ r_∘
and outer radius ϱ_o > r_∘,

    Φ_A(w) = 1/π(ϱ_o^2 - ϱ_i^2 )   ϱ_i ≤| w |≤ϱ_o,
    
      0    elsewhere.

In w-coordinates, we can find a solution to the problem

    - Δ_w p = q̂μ/ρ k_IΦ_A   in Ω_w = T(Z(V(U(Ω)))),

for a given constant well pressure p_ω, q̂ = qζ and constant density and viscosity.

By means of integration (cf. <cit.>), we get

    p(w) = 
      p_ω-μ/Kρq̂/2π[ (| w|^2 - ϱ_o^2)/2ξ^2 - ϱ_i^2/ξ^2ln(| w|/ϱ_o) + ln(ϱ_o/r_ω) ]    ϱ_i ≤| w|≤ϱ_o,
    
      p_ω-μ/Kρq̂/2π[ -1/2 - ϱ_i^2/ξ^2ln(ϱ_i/ϱ_o) + ln(ϱ_o/r_ω) ]    | w| < ϱ_i, 
    
      p_ω-μ/Kρq̂/2πln(| w|/r_ω)    | w| > ϱ_o,

where ξ^2 = ϱ_o^2 - ϱ_i^2. This shows that outside the kernel support region (| w| > ϱ_o), we obtain the exact
analytical solution derived in <ref>. Moreover, the source term can be reformulated, cf. <cit.>,

    q̂ = 2πρ k_I/μ(p_ω - p_0)Ξ,   Ξ = [ln(ϱ_o/r_ω)-1/2 - ϱ_i^2/ξ^2ln(ϱ_i/ϱ_o) ]^-1,

where p_0 := p(| w| = 0) = p(| w| = ϱ_i) is the fluid pressure evaluated on the well center-line.
Note that for f = 0 and ϱ_i = f = 0, the isotropic solution with for a circular constant kernel (<ref>) is obtained.
From the transformation of the Laplace operator, <ref>, we see that the problem

    - Δ_z p = q̂μ/ρ k_IΦ_A Φ_J    in Ω_z = Z(V(U(Ω))),

with altered kernel function Φ_Λ = Φ_AΦ_J is equivalent to <ref>.


The transformation T^-1 changes shape of the
kernel support 𝒮(Φ_A) from an annulus to an ellipse E_Φ,v. Inverting Z extrudes the solution along the well center-line,
and inverting the rotation and stretch described by V and U results in a kernel support region in the shape of
an elliptic cylinder. Moreover, each ellipse E_Φ,v with normal vector e⃗_3 is transformed to an ellipse E_Φ,x(s), that is
the intersection the elliptic cylinder with a plane with the normal vector n⃗_E_x = S̃R̃^Te⃗_3,
centered at s on the well center-line. The transformation and the normal vector n⃗_E_x are visualized in <ref>.
We note that if none of the principal axes of the permeability tensor are aligned with the well direction,
n⃗_E_x is not parallel to the well direction ψ⃗ in x-coordinates.
The integral of the right-hand side of <ref>
for a well segment Λ_i with length L_i is equal to the integral over the kernel support 𝒮(Φ_Λ,i)
which has the shape of the elliptic cylinder given by

    𝔼 := ⋃_0 ≤ s ≤ L_i E_Φ,x(s).

Using q̂ = qζ, and exploiting that k_I was chosen such that det(R̃S̃) = 1, it can be shown that

    ∫_𝔼q̂Φ_Λ,idx
    = ∫_V(U(𝔼))q̂Φ_Λ,idv
    = ∫_0^L̂_i∫_E_Φ,v(ŝ)q̂Φ_ΛdÂdŝ = q̂L̂_i = qL_i,

where 0 ≤ŝ≤L̂_i is a local coordinate along the transformed well direction, and the last equality is proven in <ref>.
This is the desired property of the kernel function for the anisotropic case corresponding to <ref> for the isotropic case.



§ NUMERICAL METHOD



We discretize <ref> using a cell-centered finite volume method with multi-point flux approximation (MPFA) <cit.>.
The domain Ω is decomposed into control volumes K_Ω∈Ω_h
such that the computational mesh Ω_h is a discrete representation of Ω.
Furthermore, each control volume boundary, ∂ K_Ω, can be split into a finite number of faces σ⊂∂ K_Ω,
such that σ = K_Ω∩ L_Ω,
with L_Ω denoting a neighboring control volume. Integrating <ref>
over a control volume K_Ω and applying the Gauss divergence theorem on the left hand side yields

    -∫_∂ K_Ω[ ρ/μ K  p ]·n_K_Ω,σ dA = ∫_K_Ωq̂Φ_Λ dx,

where n_K_Ω,σ is the unit outward-pointing normal on face σ⊂∂ K_Ω.
The exact fluxes are approximated by numerical fluxes

    F_K_Ω,σ≈ -∫_σ[ ρ/μ K  p ]·n_K_Ω,σ,

which are computed using the MPFA-O method described in <cit.>.
The discrete source term is computed as

    Q_K_Ω≈∫_K_Ωq̂Φ_Λ dx,   
    Q_K_Ω = Q_ℐ/|ℐ|∫_K_Ω∩𝒮(Φ_Λ,ℐ)Φ_Λ dx,

where Q_ℐ is a numerical approximation of the source term integral over the intersection ℐ = K_Ω∩Λ,

    Q_ℐ = |ℐ| 2πρ k_I/μ(p_ω - p_0)Ξ,

and 𝒮(Φ_Λ,ℐ) is the kernel support associated with ℐ as depicted in <ref>.
In summary, the discrete form of <ref> is

    ∑_σ⊂∂ K_Ω F_K_Ω,σ = Q_K_Ω,    K_Ω∈Ω_h.

We note that due to the dependency of Q_ℐ on p_0 the proposed method is non-local in the sense that
non-neighbor cells M_Ω∈Ω_h (where M_Ω∩ K_Ω is the empty set or a single point)
may have an associated degree of freedom that depends on the degree of freedom of K_Ω.



 §.§ Kernel integration


The kernel integral in <ref>

    𝕀_Φ,K_Ω := ∫_K_Ω∩𝒮(Φ_Λ,ℐ)Φ_Λ dx,

is not easily approximated with a quadrature rule, since the intersection K_Ω∩𝒮(Φ_Λ,ℐ), 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 <cit.>,
and remark that the integral over the entire support 𝒮(Φ_Λ,ℐ) is known exactly; see <ref>.
Hence, the integration problem can be reformulated as the distribution of the known integral over all intersected control volumes K_Ω
weighted with the respective support volume fractions. Following <cit.>,
we create n_ℐ integration points x_i ∈𝒮(Φ_Λ,ℐ)
with known volume elements V_i of similar size and shape, so that

    𝕀_Φ,K_Ω≈∑_i=1, x_i∈ K_Ω^n_ℐ V_iΦ_Λ(x_i),
      ∑_i=1^n_ℐ V_i ≈|𝒮(Φ_Λ,ℐ)|.

Computing the weights for each cell K_Ω is a pre-processing step that only has to be done once for each
computational mesh and well geometry.



§ NUMERICAL EXPERIMENTS AND DISCUSSION



We present numerical experiments using the presented method in different setups.
All experiments are conducted with constant fluid density ρ = 1000 and viscosity μ = 1e-3.
The well pressure is constant, p_ω = 1e6, and the well radius is r_ω = 0.1 if not specified otherwise.
The permeability tensor is given as

    K(γ_1, γ_2) = R_1(γ_1) R_2(γ_2) K_α R^T_2(γ_2) R_1^T(γ_1),
          K_α = [[ 1 0 0; 0 1 0; 0 0 α ]] 1e-12□,

where

    R_1(γ_1) = [[       1       0       0;       0  cosγ_1 -sinγ_1;       0  sinγ_1  cosγ_1 ]],  
      R_2(γ_2) = [[  cosγ_2       0  sinγ_2;       0       1       0; -sinγ_2       0  cosγ_2 ]]

are rotation matrices rotating vectors about e⃗_1, e⃗_2 by the rotation angle γ_1, γ_2, respectively,
and α is a given dimensionless K-anisotropy ratio α = K_33/K_11 = K_33/K_22.
The domain Ω_0 = [-100,100]×[-100,100]×[-50,150] is split in two regions,
Ω = [-100,100]×[-100,100]×[0,100] and Ω_D = Ω_0∖Ω.
The well center-line Λ is given by the line through the origin and ψ⃗ = R_1(β_1) R_2(β_2) e⃗_3,
where R_1, R_2 are given in <ref> and β_1, β_2, are rotation angles.
The analytical solution for all cases is given in <ref>, q = 1, and L = |Λ∩Ω| (in ).
For all setups the inner kernel radius is chosen as ϱ_i = f.
In all of Ω_D and on the boundary ∂Ω the analytical solution is enforced by Dirichlet constraints, modeling the
infinite well.
The computational mesh Ω_h is a structured grid composed of regular hexahedra K_Ω.
Furthermore, we define two error measures.

    E_p = 1/p_ω[1/|Ω_h |∑_K_Ω∈Ω_h| K_Ω|( p_e,x⃗_K_Ω - p_K_Ω)^2 ]^1/2

is the relative discrete L^2-norm of the pressure, where p_e,x⃗_K_Ω
is the exact pressure evaluated at the cell centroid x⃗_K_Ω and p_K_Ω the discrete numerical cell pressure, and

    E_q = 1/q[ 1/|Λ∩Ω_h |∑_K_Ω∈Ω_h 
     K_Ω∩Λ≠∅|ℐ|( q - Q_ℐ/|ℐ|ζ)^2 ]^1/2

is the relative discrete L^2-norm of the source term, where |ℐ| = | K_Ω∩Λ|
is the length of the intersection of cell K_Ω and the well center-line Λ, Q_ℐ is the discrete source term
given in <ref>.
All setups are implemented in  <cit.>, an open-source porous media simulator based on  <cit.>.



 §.§ Grid convergence for different anisotropy ratios




In the first numerical experiment grid convergence is investigated for different anisotropy ratios α.
To this end, h_max := max_K_Ω∈Ω_h h_K_Ω,
where h_K_Ω is defined as the maximum distance between two vertices of the cell K_Ω. Starting at
a grid resolution for Ω_h of 20×20×10 cells (h_max = 10√(3)), the grid is refined uniformly.
<Ref> shows the errors E_p and E_q for different grid resolutions and values of α,
for β_1 = β_2 = 20^∘ and γ_1 = γ_2 = -20^∘, so that K is a full tensor and none of the
principal axis of K 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 <cit.> (super convergence at cell centers).
The source term q is a linear function of the pressure p and also exhibits second order convergence.
We note that the errors for different α are not directly comparable since the analytical solution for p
changes with α, although q is constant.

However, the convergence rates are shown to be independent of α with increasing grid resolution.
The convergence rates for E_q (slope of the lines in <ref>) are presented in <ref>.
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 α=100, the kernel
ellipse in x-coordinates has major and minor axis of a_x ≈55.9m, b_x ≈5.6,
respectively, while h_max≈17.32 for the lowest grid resolution.



 §.§ Influence of the outer kernel radius ϱ_o


In <cit.>, it is suggested that increasing the kernel support region (increasing ϱ_o), has
a similar effect on E_q 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 (ϱ_o →ϱ_i, Φ_Λ→δ_Λ).
However, every discrete cell K_Ω 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 K_Ω.

As shown in <cit.> and <ref>, isobars become circular, in the transformed domain U(Ω),
with increasing distance to the well. Therefore, a reasonable simplification is Φ_J ≈ 4 if ϱ_o ≫ r_ω.
This is completely analogous to the assumption of circular isobars in <cit.>, 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 ϱ_o, for the same
20×20×10 grid. This is done once for the case, where ϱ_o ≫ r_ω and for the case for which
ϱ_o is only slightly larger than r_ω.
Furthermore, β_1 = β_2 = 20^∘ and γ_1 = γ_2 = -20^∘. The results are
shown in <ref>. First, it can be seen that doubling ϱ_o leads to a 4-times smaller
error E_q. 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 p_0.
Furthermore, the result is consistent with the results in <cit.>.
Moreover, <ref> suggests that for ϱ_o ≫ r_ω the simplification of the kernel function
(Φ_J ≈ 4) is not visible in E_q, while for kernel radii slightly smaller than the well radius,
the simplification increases E_q by an order of magnitude in comparison to the case using the exact kernel
function as derived in <ref>.
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. <cit.>,
is therefore well-conceivable.



 §.§ Robustness with respect to rotation

In the following numerical experiment, we use a single computational mesh with a given resolution for Ω_h: 20×20×10.
First, the well direction is fixed, and the permeability tensor is rotated by varying γ_1 and γ_2. Then the permeability
tensor is fixed and the well is rotated by varying β_1 and β_2. The results are shown in <ref>.

It can be seen that the presented well model is rather robust with respect to rotations. Possible effects influencing the
approximation error E_q, 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 d⃗_K_Ω,σ = Kn⃗_K_Ω,σ.
Additionally, for different well angles the number and size of intersections
K_Ω∩Λ can have an influence on the discrete error.



 §.§ Comparison with a Peaceman-type well model



In particular for petroleum engineering applications, commercial codes typically use Peaceman-type well models <cit.>.
In <cit.>, 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 K_Ω is approximated by

    Q_K_Ω   = 2πρ/μ(p_ω - p_0)L_K_Ω√(K_11 K_22)/ln(r_0/r_ω),
    
    r_0    = e^-γ/2[ ( K_22/K_11)^1/2Δ x^2 + ( K_11/K_22)^1/2Δ y^2 ]^1/2/( K_11/K_22)^1/4 + ( K_22/K_11)^1/4,

where Δ x and Δ y are the horizontal dimensions of the cell containing the well,
L_K_Ω = |ℐ| the length of the well segment contained in K_Ω, 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 K_Ω.
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 <cit.>. The
Alvestad model <cit.> has been adapted for finite volumes, for example in <cit.>
(formula given in <ref>, 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 <cit.>.

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 <ref>.
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 <ref> without rotation of the permeability tensor) show large deviations
(>10 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 Ω = [-50,50]×[-100,100]×[0,100]
contains a slanted straight well Λ
with end points at x⃗_Λ,1 = [-20, -50, 25]^T, x⃗_Λ,2 = [20, 50, 75]^T.
The permeability tensor is a diagonal tensor K(γ_1,γ_2), with γ_1 = 0^∘, γ_2 = 90^∘, α = 0.1.
The structured cube grid Ω_h is successively, uniformly refined starting with 10×20×10 cells (h_max≈17.32).
The well radius is r_ω = 0.1 (Δ x / r_ω = 200 for the coarsest grid).
The kernel support region (chosen as ϱ_o / r_ω = 100) 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 (K  p) ·n⃗ = 0, except for the planes perpendicular to the x_1-axis, where the Dirichlet boundary condition p_D(x_2 = -100) = 1e5, p_D(x_2 = 100) = 3e5 are enforced.
The reference solution is computed
with 160×320×160 cells (h_max≈1.08).
The computational domain with pressure iso-surfaces of the reference solution are shown in <ref>.

In <ref>, the relative integral source error

    E_Q = | Q - Q_ref|/| Q_ref|,    Q = ∑_K_Ω∈Ω_h Q_K_Ω,

with respect to the reference
solution Q_ref 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 ϱ_o / r_ω = 100, the extent of the
kernel ellipse E_Φ,x  is given by its major and minor axes,
16.12 and 12.54. 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 h_max,
so that for the finest grid shown in <ref> (80×160×80 cells, h_max≈2.17),
the E_Φ,x major and minor axes measure 2.01 and 1.57.

It is evident that the numerical solution for the distributed source model converges to the reference solution. More importantly,
the relative error is small (<0.5) even for the coarsest grid. In comparison, the difference to the Peaceman-type model is
large (>5). In particular, the error grows with grid refinement (to >8), signifying that the generalization of
Peaceman's model for arbitrarily-oriented wells is not consistent. The result is comparable with the observations in <cit.>,
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 r_ω / Δ x ratio. In the variant of the ds model, the error in the source term
with respect to the reference solution also grows with larger r_ω / Δ x ratio. However, the error is consistently smaller (by a factor >3)
than for the pm well model.

<Ref> shows the numerical pressure solutions along the x_1 and the x_2-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 <ref>).
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.



§ SUMMARY


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 (< 2 error with respect to a reference solution),
whereas the Peaceman-type well model for the same case showed larger differences (> 8)
which also had a negative global effect on the pressure solution. The example showed that the new model gives the choice between a more accurate representation
of the near-well pressure and a more accurate representation of the source term. Additionally, it was shown that if the source term is accurately approximated,
the far-field pressure (outside kernel support) is equally well-approximated while the regularization of the pressure profile only happens locally in the well neighborhood.
On the other hand, a bad approximation of the source term leads to global errors in the pressure profile.
Finally, the new model was shown to be robust with respect to well rotation, as well as robust with respect to rotations of the anisotropic permeability tensor.

In this work, the well model derivation is restricted to one-phase flow and possible modifications for multi-phase flow are yet to be explored.
An extension of the well model for wells with casing is easily conceivable, combining the findings in this work with the derivations presented in <cit.>.
Moreover, the herein presented cases considered the case of a given constant bottom hole pressure. However, the results of <cit.> indicate,
that the presented model may be extended for the coupled 1d-3d case where the well pressure solves an additional one-dimensional
partial differential equation in the well domain. Finally, using the superposition principle as shown in <cit.>, which equally applies for
the case of anisotropic permeabilities, the presented well model is also expected to provide good approximations when multiple wells are present in the domain.



§ ACKNOWLEDGEMENTS


This work was financially supported by the German Research Foundation (DGF), within the Cluster of Excellence in Simulation Technology (EXC 310),
and the Collaborative Research Center on Interface-Driven Multi-Field Processes in Porous Media (SFB 1313, Project Number 327154368).


sectionappsec



§ RODRIGUES' ROTATION FORMULA


We want to rotate a given basis ℬ = {e⃗_1, e⃗_2, e⃗_3 } such that e⃗_3 is aligned with a vector ψ⃗.
The Rodrigues' rotation formula <cit.> describes a rotated vector x⃗_rot obtained by rotating a vector x⃗ by the angle θ about an axis given by the unit normal vector k⃗

    x⃗_rot = x⃗cosθ + (k⃗×x⃗) sinθ + k⃗(k⃗·x⃗)(1 - cosθ).

The desired rotation can be described by a rotation by an angle θ = π about the axis k⃗ = (e⃗_3 + ψ⃗)/2|e⃗_3 + ψ⃗|,

    x⃗_rot,π = 2k⃗(k⃗·x⃗) -x⃗.

<Ref> can be expressed in matrix notation as x⃗_rot,π = Rx⃗ with

    R = (2 k⃗k⃗^T - I),    R = R^T = R^-1,   det(R) = 1.




§ TRANSFORMATION OF LAPLACE OPERATOR


The Joukowsky transformation is a complex function z = T^-1(w) that can be decomposed in
its real and imaginary parts, x = (z) and y = (z). Furthermore, let u = (w) and v = (w).
As a conformal mapping, z satisfies the Cauchy Riemann equations <cit.>

    ∂ x/∂ u = ∂ y/∂ v  and  ∂ x/∂ v = -∂ y/∂ u.

The Jacobian of the transformation J_T^-1 is given by <ref>.
We investigation the effect of the transformation on the Laplace operator

    Δ_w p = ∂^2 p/∂ u^2 + ∂^2 p/∂ v^2,

where p is analytic in Ω_w.
Applying the chain rule yields

    ∂ p/∂ u   = ∂ p/∂ x∂ x/∂ u + ∂ p/∂ y∂ y/∂ u, 
    ∂^2 p/∂ u^2   = ∂ p/∂ x∂^2 x/∂ u^2 + ∂ p/∂ y∂^2 y/∂ u^2
                                         + ∂/∂ x(∂ p/∂ u)∂ x/∂ u
                                         + ∂/∂ y(∂ p/∂ u)∂ y/∂ u
       = ∂ p/∂ x∂^2 x/∂ u^2 + ∂ p/∂ y∂^2 y/∂ u^2
                                         + ∂^2 p/∂ x^2( ∂ x/∂ u)^2
                                         + 2∂^2 p/∂ x ∂ y∂ x/∂ u∂ y/∂ u
                                         + ∂^2 p/∂ y^2( ∂ y/∂ u)^2.

Analogously, we arrive at a similar expression for ∂^2 p / ∂ v^2. Using <ref> and

    ∂^2 y/∂ v^2 = ∂/∂ v(∂ y/∂ v)
                                      (B.1)∂/∂ v(∂ x/∂ u)
                                      = ∂/∂ u(∂ x/∂ v)
                                      (B.1) -∂ y^2/∂ u^2,   ∂^2 x/∂ v^2 = -∂^2 x/∂ u^2,

we find that

    ∂^2 p/∂ u^2 + ∂^2 p/∂ v^2
      = [ ( ∂ x/∂ u)^2 + ( ∂ y/∂ u)^2 ] [ ∂^2 p/∂ x^2 + ∂^2 p/∂ y^2].

With the complex derivative of T^-1 <cit.>,

    ∂ T^-1/∂ w = ∂ T^-1/∂ u = ∂ x/∂ u + i∂ y/∂ u  and  | ∂ T^-1/∂ w| = √(( ∂ x/∂ u)^2 + ( ∂ y/∂ u)^2).

From the determinant of the Jacobian of the transformation, we find

    det(J_T^-1) = ∂ x/∂ u∂ y/∂ v - ∂ x/∂ u∂ y/∂ v = (∂ x/∂ u)^2 + ( ∂ y/∂ u)^2 = | ∂ T^-1/∂ w|^2,

using <ref>. Hence,

    Δ_w p = | ∂ T^-1/∂ w|^2 Δ_z p = |det(J_T^-1)| Δ_z p,

which also proofs that any harmonic function (Δ f = 0) yields another harmonic function after a coordinate transformation with a conformal mapping. 



§ SOURCE SCALING FACTOR IN W-COORDINATES


We want to construct a pressure solution in w-coordinates such that the total mass flux over the well boundary
matches the specified boundary condition in x-coordinates.
Hence, the total mass flux over the boundary of a well segment with length L̂ in w-coordinates needs to match the total mass
flux over the boundary of a well segment with length L in x-coordinates.
A q̂ has to be chosen such that qL = q̂L̂. A relation between L and L̂ can be
derived by looking at two related volume integrals. The Joukowsky transformation only affects the two-dimensional well-bore
plane such that length L̂ of a well segment is not affected.
The volume of that well segment in v-coordinates (an elliptic cylinder) is given by

    V_v = π ab L̂.

The volume of the same well segment in x-coordinates (a circular cylinder with slanted parallel planar elliptic caps)
is given by

    V_x = |ψ⃗^TS̃R̃^Te⃗_3 || E_ω,x| L = π r_ω^2 L,

where
| E_ω,x| is the area of the well-bore ellipse described by <ref> transformed to x-coordinates
(as shown in <ref>) and the last
equality uses the fact that the integral can be transformed to an integral over a regular cylinder with radius r_ω and length L.
From the transformation theorem, we know that V_x = V_v det(S̃). As the parameter k_I
is chosen such that det(S̃) = 1,

    L̂ = Lr_ω^2/ab,

and if the source term is chosen as

    q̂ = qab/r_ω^2 := qζ,

then qL = q̂L̂. 



§ DETERMINANT OF THE JOUKOWSKY TRANSFORMATION


In <ref>, we show that |det(J_T^-1)| = | ∂ z/∂ w|^2.
Using complex differentiation,

    ∂ z/∂ w = ∂ z/∂ u = ∂/∂ u[ 1/2( w + f^2/w)]
    =  1/2( 1 - f^2/w^2) = 1/2( 1 - f^2 w^2/| w|^4),

where we used the identities w^-1 = w| w|^-2, w denoting the complex conjugate of w, and w^2 = w^2.
Furthermore,

    (∂ z/∂ w)^2    = 1/4( 1 - 2f^2 (w^2)/| w|^4 + f^4(w^2)^2/| w|^8), 
    (∂ z/∂ w)^2    = 1/4(f^4 (w^2)^2/| w|^8),   (w^2)^2 + (w^2)^2 =  | w|^4,

such that

    | ∂ z/∂ w|^2 = (∂ z/∂ w)^2 + (∂ z/∂ w)^2
    = 1/4( 1 + f^4 - 2f^2(w^2)/| w |^4).





§ EXTENSION OF PEACEMAN WELL MODEL FOR SLANTED WELLS


For a given well direction ψ⃗ = [ψ_1, ψ_2, ψ_3]^T and a cell K_Ω
with dimensions Δ x ×Δ y ×Δ z,
to obtain the generalized well model due to <cit.>,
reformulated for cell-centered finite volume schemes in <cit.>,
replace k = √(K_11 K_22) in <ref> by

    k = (ψ_1^2 K_22 K_33 + ψ_2^2 K_11 K_33 + ψ_3^2 K_11 K_22)^1/2

and the expression for r_0 by

    r_0    = e^-γ/2( Δ L_1^2  + Δ L_2^2 )^1/2/(√(A_1) + √(A_2)),   with
    Δ L_1^2    = ( K_22/K_33)^1/2Δ z^2 ψ_1^2
                   + ( K_33/K_11)^1/2Δ x^2 ψ_2^2
                   + ( K_11/K_22)^1/2Δ y^2 ψ_3^2, 
    Δ L_2^2    = ( K_33/K_22)^1/2Δ y^2 ψ_1^2
                   + ( K_11/K_33)^1/2Δ z^2 ψ_2^2
                   + ( K_22/K_11)^1/2Δ x^2 ψ_3^2, 
    
    A_1    = ( K_22/K_33)^1/2ψ_1^2
            + ( K_33/K_11)^1/2ψ_2^2
            + ( K_11/K_22)^1/2ψ_3^2, 
    
    A_2    = ( K_33/K_22)^1/2ψ_1^2
            + ( K_11/K_33)^1/2ψ_2^2
            + ( K_22/K_11)^1/2ψ_3^2.

We note that the formula reduces to <ref>,
if ψ⃗ is aligned with a coordinate axis.




elsarticle-num

