A new and consistent well model for one-phase flow in anisotropic porous media using a distributed source model
AA new and consistent well model for one-phase flow inanisotropic porous media using a distributed sourcemodel
Timo Koch a, ∗ , Rainer Helmig a , Martin Schneider a a Department of Hydromechanics and Modelling of Hydrosystems, University of Stuttgart,Pfaffenwaldring 61, 70569 Stuttgart, Germany
Abstract
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 byspatially distributed source terms over a small neighborhood region. To thisend, we first present a compact derivation of the exact analytical solution foran arbitrarily oriented, infinite well cylinder in an infinite porous medium withanisotropic permeability tensor in R , for constant well pressure and a giveninjection rate, using a conformal map. The analytical solution motivates thechoice of a kernel function to distribute the sources. The presented model isindependent from the discretization method and the choice of computationalgrids. In numerical experiments, the new well model is shown to be consistentand robust with respect to rotation of the well axis, rotation of the permeabilitytensor, and different anisotropy ratios. Finally, a comparison with a Peaceman-type well model suggests that the new scheme leads to an increased accuracyfor injection (and production) rates for arbitrarily-oriented pressure-controlledwells. Keywords: well model, 1d-3d, mixed-dimension, anisotropic, analytic solution,Peaceman
1. Introduction
Well modeling is essential for various engineering applications, as for exam-ple reservoir simulation, geothermal energy production or energy storage, whereinjection or extraction processes strongly influence the flow behavior. Usuallythe well geometry is not explicitly resolved in the mesh but instead modeled asa line source with given extraction or injection rate. However, this simplifiedapproach introduces singularities, meaning that the logarithmic solution profilesare undefined at the center-line of the well. This leads to a significant deviation ∗ Corresponding author
Email address: [email protected] (Timo Koch)
Preprint July 31, 2019 a r X i v : . [ c s . C E ] J u l etween numerical and analytical solution in the near-well region. For a bet-ter approximation, locally refined meshes around the wells are needed, whichhowever deteriorate efficiency and is therefore often not suitable for field-scalesimulations, especially when multiple wells are present. Similar issues are en-countered for the modeling of vascularized biological tissue perfusion [1, 2] orthe modeling of root water update [3, 4].A common approach is the use of well-index-based well models. Such wellmodels aim to find a relation between well rate, bottom hole pressure and nu-merically calculated pressure (well-block pressure) for each cell (grid-block) thatcontains the well. In reservoir engineering such a relation is denoted as wellindex. The first theoretical derivation of a well index for two-dimensional struc-tured uniform grids with isotropic permeability has been presented by Peace-man [5]. He has shown that the well-block pressure differs from the areal aver-aged analytical pressure and introduced a new relation by using an equivalentwell radius. The equivalent well radius is defined as the distance (relative tothe well location) at which the analytical and numerical pressures are equal.A generalization for structured two-dimensional non-square grids (∆ x (cid:54) = ∆ y )and anisotropic but diagonal permeability tensors has been presented by thesame author in [6]. However, these well models are restricted to uniform gridswhere the well is oriented with one of the grid axes. Furthermore, since theeffective well radius relates numerical and analytical pressure values, the wellindex has to be calculated depending on the used discretization scheme. Thus,Peaceman’s well model is only valid for a cell-centered finite difference schemewith 5-point stencil. Discussion of other discretization schemes can be foundin [7], again with the restriction of two-dimensional grids. Enhancements in-clude, among others, three-dimensional slanted wells [8, 9], Green’s functionsfor the computation of well indices [10, 11, 12], or the singularity subtractionmethod to obtain smooth solutions in the near-well region [13, 14].In this work, a new approach for obtaining more accurate source term fora given well bottom hole pressure is presented. The new model is, in contrastto most of the existence methods, independent of the discretization scheme andcan be used for general unstructured grids. Additionally, the presented methodis not restricted to diagonal tensors and thus works for general anisotropic per-meabilities. In Section 2, we derive a well model, initially for isotropic porousmedia, for which the fluid mass injected by a well is distributed over a smallneighborhood around the well, using kernel functions. The derivation followsthe idea recently presented in [15], however we herein discuss the case withoutmembrane or casing. The model yields a pressure solution without singularity,from which the source term can be reconstructed using a relation found withthe analytical solution for the case of an infinite well in an infinite medium. Themodel generalizes to more complex problems due to the superposition principalvalid for the Laplace operator, a linear operator [15]. In Section 3, the modelis generalized to porous media with general anisotropic permeabilities, basedon an analytical solution constructed in Section 3.1 using a series of coordinatetransformations. We shown that the general model reduces to the model de-rived in Section 2 for isotropic permeabilities. Finally, the new well model is2nalyzed with several numerical experiments in Section 5. The results indicatethat the model is consistent for different anisotropy ratios, robust with respectto rotations of the well relative to the computational grid, and to rotationsof an anisotropic permeability tensor. A comparison with a Peaceman-typewell model in a setup with a K-orthogonal grid and an embedded slanted wellsuggests that the new model more accurately approximates the fluid exchangebetween 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 withisotropic permeability tensor, not including a well casing. The derivation fol-lows [15], where a casing is included in form of a membrane in a different butrelated application, that is modeling fluid exchange between the vascular sys-tem and the embedding biological tissue. Stationary single-phase flow arounda well with radius r ω , in an isotropic porous medium with permeability k , canbe described by the following flow equation − ∇ · (cid:18) ρµ k ∇ p (cid:19) = q Φ Λ in Ω , (2.1)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 (kg s − m − )over a small tubular support region, S (Φ Λ i ), with radius (cid:37) ( s ), around a wellsegment i , such that Φ Λ i = 0 outside the support region. We choose kernelfunctions Φ Λ i ( s ) with the property L i (cid:90) π (cid:90) (cid:37) ( s ) (cid:90) Φ Λ i r d r d θ d s = L i , (2.2)where r , θ , s are the radial, angular, and axial coordinate in a segment-localcylinder coordinate system, and L i is the length of segment i .Assume a radially symmetric zone (distance δ > r ω from center-line) aroundthe 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 < δ isdescribed by the analytical solution p ( r ) = − µkρ q π ln r + C. (2.3)The constant C is determined by fixing a well pressure, p ω , p ω = p ( r ω ) = − µkρ q π ln r ω + C ⇒ C = p ω + µkρ q π ln r ω , so that p ( r ) = − µkρ q π ln (cid:18) rr ω (cid:19) + p ω . (2.4)3 igure 1: An infiltration scenario for an isotropic porous medium. Schematic representationof the introduced symbols. An infinite well with radius r ω , center-line Λ with local cylindricalcoordinate system ( r, θ, s ) is embedded in the porous domain Ω. The kernel function withradius (cid:37) regularizes the pressure solution which can then be evaluated at r = 0: p ( r = 0) = p . Consequently, the source term can be expressed in terms of p ω and p δ as q = 2 πr ω ρkµ ( p ω − p δ ) r ω ln (cid:16) δr ω (cid:17) (2.5)We choose a simple kernel function which regularizes the pressure solution for r ≤ (cid:37) , Φ const ( r ) = (cid:40) π(cid:37) r ≤ (cid:37), r > (cid:37), (2.6)where (cid:37) ≤ δ . The pressure for r < ρ can be obtained by integration fromEq. (2.1), yielding p ( r ) = − µkρ q π (cid:104) r (cid:37) + ln (cid:16) (cid:37)r ω (cid:17) − (cid:105) + p ω r ≤ (cid:37), − µkρ q π ln (cid:16) rr ω (cid:17) + p ω r > (cid:37). (2.7)Figure 1 graphically explains the most important symbols introduced in thissection. As the regularized pressure can be evaluated at the well center-line wecan reformulate Eq. (2.5), q = 2 π ρkµ ( p ω − p )Ξ , with Ξ = ( p ω − p δ )( p ω − p ) 1ln (cid:16) δr ω (cid:17) (2.8)where Ξ is the so called flux scaling factor. The flux scaling factor can be alsoexpressed independent of the pressure. To this end, Eq. (2.7) is evaluated at r = 0, so that p ω is expressed in terms of p , p = − ( p ω − p )Ξ (cid:20) ln (cid:18) (cid:37)r ω (cid:19) − (cid:21) + p ω , (2.9)4here q was replaced by inserting Eq. (2.8). It directly follows from Eq. (2.9)that Ξ = (cid:20) ln (cid:18) (cid:37)r ω (cid:19) − (cid:21) − . (2.10)
3. A well model with distributed source for anisotropic media
In the following section, the developed well model is extended for porousmedia with anisotropic permeability. In Section 3.1, we derive an analyticalsolution for one-phase flow around an infinitely long cylindrical well embeddedin an infinite porous domain in R . This derivation motivates the choice of asuitable kernel function for anisotropic problems, presented in Section 3.3. In the following section, we derive an analytical solution for one-phase flowaround an infinite cylindrical well Γ with radius r ω in an infinite porous domainˆΩ = R \ Γ with anisotropic, homogeneous permeability. We assume, withoutloss of generality, that the well axis passes through the origin of the Cartesiancoordinate system, and denote by ψ a unit vector parallel to the well signifyingthe well orientation. We seek an analytical expression for the hydraulic pressure p such that − ∇ · (cid:18) ρµ K ∇ p (cid:19) = 0 in ˆΩ , (3.1)for a constant well pressure p ω in Pa and some specific pumping rate q inkg s − m − given on ∂ Γ. The total mass flow over the boundary of a wellsegment 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 = QDQ T , (3.2)where D = diag( λ , λ , λ ) is a diagonal matrix composed of the eigenvalues λ i of K , Q = [ ν K, | ν K, | ν K, ] is a rotation matrix with the corresponding eigen-vectors as columns, and A T denotes the transpose of a matrix A . Further usefulproperties derived from the decomposition are det( K ) = λ λ λ , where det( A )denotes the determinant of A , and K n = Q Λ n Q T , where D r = diag( λ r , λ r , λ r ), r ∈ R .It is well known that the anisotropic one-phase flow problem can be trans-formed to an isotropic problem using a coordinate transformation [9, 16, 17, 6,18] U : R → R , x (cid:55)→ u = ˜ S x , (3.3)with the stretching matrix ˜ S = k / I K − / , where k I is an arbitrary scalarconstant, that we choose as k I = det( K ) − / (cf. [9]), rendering the transfor-mation isochoric. The transformation u = ˜ S x deforms the well cylinder such5hat a cross-section orthogonal to the transformed well direction is elliptical.The solution to the isotropic problem − ∇ u · (cid:18) ρµ k I ∇ u p (cid:19) = 0 in ˆΩ u = U ( ˆΩ) , (3.4)is identical on two parallel planes perpendicular to the transformed (normalized)well direction, ψ (cid:48) = ˜ S ψ || ˜ S ψ || − . This motivates the rotation of the coordinatesystem such that the first and second axis are aligned with the major and minoraxis of the well-bore ellipse and third axis is aligned with ψ (cid:48) . To determine thecorresponding rotation matrix ˜ R , we need to characterize this well-bore ellipse.The well cylinder in x -coordinates is given by x T Ψ x = r ω , Ψ = I − ψψ T . (3.5)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, projectinginto the plane perpendicular to the well direction yields the well-bore ellipseequation ˆ v T E ˆ v = ˆ v T P T R ˜ S − Ψ ˜ S − R T P ˆ v = r ω , ˜ S − R T P ˆ v = x (3.6)in ˆ v -coordinates, where R = 2 ( e + ψ (cid:48) )( e + ψ (cid:48) ) T ( e + ψ (cid:48) ) T ( e + ψ (cid:48) ) − I, e = , P = . (3.7)The rotation matrix R can be derived using Rodrigues’ rotation formula asshown in Appendix A. The length of the major and minor ellipse axis are foundas a = r w γ − / and b = r w γ − / , where γ i are the eigenvalues of E , and theaxis orientations are given by ν = P ˆ ν E, , ν = P ˆ ν E, , where ˆ ν E,i denote thecorresponding eigenvectors of E . We assume that the eigenvalues and eigenvec-tors are sorted such that a ≥ b , and oriented such that ψ (cid:48) = ν × ν . Finallythe desired rotation is given by V : R → R , u (cid:55)→ v = ˜ R u = ˆ R T R T u , (3.8)where ˆ R = (cid:20) ν (cid:12)(cid:12)(cid:12)(cid:12) ν (cid:12)(cid:12)(cid:12)(cid:12) ψ (cid:48) (cid:21) (3.9)is rotating about the well direction axis such that the coordinate system isaligned with the principal ellipse axes.Following the derivations from above, we now have to solve a two-dimensionalisotropic Laplace problem with boundary conditions prescribed on an ellipse.To this end, we note that the transformation of a harmonic function f (a func-tion satisfying Laplace’s equation ∆ f = ∇ ·∇ f = 0) with a conformal (angle-preserving) mapping yields another harmonic function [19] (see Appendix B).6sing a Joukowsky transformation, a conformal mapping well-known from aero-dynamics [20], the isotropic problem with a well with elliptic cross-sections, canbe transformed to an isotropic problem with circular cross-sections [16]. Trans-forming into the complex plane (parametrizing the well-bore ellipse plane) Z : R → C , v (cid:55)→ z = ˜ Z v = [1 , i, v = v + iv , (3.10)the (inverse) Joukowsky transformation T : C → C , z (cid:55)→ w = z + (cid:112) z − f (cid:112) z + f , f = (cid:112) a − b (3.11)transforms elliptic isobars into circular isobars, where a and b , a ≥ b , are themajor and minor axis of the well-bore ellipse, as derived above. In particular, thewell-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 symmetricanalytical solution to problem Eq. (3.1) p ( w ) = p ω − µρk I ˆ q π ln (cid:18) | w | r ◦ (cid:19) ζ, ˆ q = qζ = q abr ω , (3.12)where the source scaling factor ζ is necessary to recover the original source q on ∂ Γ. This can be derived from simple geometric considerations as shownin Appendix C. The w corresponding to some x ∈ ˆΩ in original coordinatesis obtained by using all above-mentioned transformations after each other asfollows w = T ( Z ( V ( U ( x )))) = T ( ˜ Z ˜ R ˜ S x ) . (3.13)Such a solution for a slanted well (30 ° with respect to vertical axis) and anisotropicpermeability tensor K A = · − m (3.14)is visualized in Fig. 2. To construct a suitable kernel function for anisotropic problems, we first havea closer look at the properties of the employed Joukowsky transformation. Theeffect of the mapping T , Eq. (3.11), is shown in Fig. 3. Points on the exteriorof a line on the real axis between f and − f are mapped onto the exterior of acircle with radius f . The ellipse with major axis a and minor axis b is mappedonto a circle with radius r ◦ = a + b . Going further away from the well, thedeformation due to the mapping is less and less pronounced. This matches theexpectations for the physical flow problem, since isobars at large distance froman elliptical well-bore become increasingly circular in isotropic media.The inverse transformation is given by T − : C → C , w (cid:55)→ z = 12 (cid:18) w + f w (cid:19) , | w | > f, (3.15)7 igure 2: The analytical pressure solution for a slanted well, r ω = 0 . p ω = 5 . · Pa, total mass injection rate Q ω = 115 .
47 kg s − , and anisotropic permeabilitytensor K A . The top view is oriented in well direction and shows pressure contour surfaceshighlighting their the elliptical shape. Figure 3: Visualization of the (inverse) Joukowsky transformation, T , exemplarily for a = 1 . b = 0 .
9, and thus f = 1 .
2. The points on the left are shown on the complex z -plane, the pointson the right are shown in the complex w -plane, and w = T ( z ). The inner circle on the rightimage has radius f . The other circles have radii of r ◦ = a + b and 3 r ◦ . | w | is necessary to obtain a one-to-one mapping. Thetransformation T − can equally be interpreted as an R → R mapping. TheJacobian of the transformation T − for z = x + iy and w = u + iv has the form J T − = (cid:20) ∂x∂u ∂x∂u∂y∂v ∂y∂v (cid:21) = (cid:20) η (cid:15) − (cid:15) η (cid:21) , (3.16)which follows from the the Cauchy–Riemann equations [21]. Since the trans-formation can be viewed as the composition of a scaling and a rotation, it isangle-preserving. The transformation is associated with a spatially dependentvolume deformation characterized by the determinant of J T − . Furthermore, itcan be shown that the Laplace operator behaves as follows under the transfor-mation z = T − ( w ),∆ w p = ∂ p∂u + ∂ p∂v = (cid:12)(cid:12)(cid:12)(cid:12) ∂T − ∂w (cid:12)(cid:12)(cid:12)(cid:12) ∆ z p = | det( J T − ) | ∆ z p (3.17)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 com-pleteness in Appendix B. From Eq. (3.17) follows that | det( J T − ) | − ∆ w p = ∆ z p, (3.18)for the transformation w = T ( z ). The determinant can be explicitly computed,using Eq. (3.17) and complex differentiation (shown in Appendix D) as | det( J T − )( w ) | = (cid:12)(cid:12)(cid:12)(cid:12) ∂T − ∂w (cid:12)(cid:12)(cid:12)(cid:12) = (cid:12)(cid:12)(cid:12)(cid:12) ∂z∂u (cid:12)(cid:12)(cid:12)(cid:12) = 14 (cid:18) f − f (cid:60) ( w ) | w | (cid:19) := Φ − J , (3.19)where (cid:60) ( w ) is real part of w , and | w | the absolute value of w . We note thatΦ J quickly converges to the value 4 with increasing | w | , that is with increasingdistance from the well. The function Φ − J is plotted in Fig. 4 exemplarily for f = 1 . Instead of excluding the well domain Γ from Ω = R and modeling infiltra-tion or extraction by a flux boundary condition, we will now model the actionof the well on the flow field by a spatially distributed source term, as presentedfor the isotropic problem, − ∇ · (cid:18) ρµ K ∇ p (cid:19) = qζ Φ Λ in Ω . (3.20)From the above derivations, we know that solving Eq. (3.20) in w -coordinates isstraight-forward. Hence, we choose kernel functions in w -coordinates and thentransform to x -coordinates so that the pressure solution satisfies Eq. (3.20).Motivated by the properties of the Joukowsky transform (see Fig. 3), we choose9 J = det( J T ) for f = 1.2 Figure 4: The determinant of the Jacobian of the transformation z = T − ( w ) for | w | > f where f = a − b and a > b . Larger volume deformations only occur very locally in vicinityof the well radius r ◦ = a + b and quickly converge to 0 .
25 with larger distance to the well. a local kernel that is constant on the annulus with inner radius f < (cid:37) i ≤ r ◦ andouter radius (cid:37) o > r ◦ , Φ A ( w ) = (cid:40) π ( (cid:37) o − (cid:37) i ) (cid:37) i ≤ | w | ≤ (cid:37) o , . (3.21)In w -coordinates, we can find a solution to the problem − ∆ w p = ˆ q µρk I Φ A in Ω w = T ( Z ( V ( U (Ω)))) , (3.22)for a given constant well pressure p ω , ˆ q = qζ and constant density and viscosity.By means of integration (cf. [15]), we get p ( w ) = p ω − µKρ ˆ q π (cid:104) ( | w | − (cid:37) o )2 ξ − (cid:37) i ξ ln (cid:16) | w | (cid:37) o (cid:17) + ln (cid:16) (cid:37) o r ω (cid:17)(cid:105) (cid:37) i ≤ | w | ≤ (cid:37) o ,p ω − µKρ ˆ q π (cid:104) − − (cid:37) i ξ ln (cid:16) (cid:37) i (cid:37) o (cid:17) + ln (cid:16) (cid:37) o r ω (cid:17)(cid:105) | w | < (cid:37) i ,p ω − µKρ ˆ q π ln (cid:16) | w | r ω (cid:17) | w | > (cid:37) o , (3.23)where ξ = (cid:37) o − (cid:37) i . This shows that outside the kernel support region ( | w | > (cid:37) o ),we obtain the exact analytical solution derived in Section 3.1. Moreover, thesource term can be reformulated, cf. [15],ˆ q = 2 π ρk I µ ( p ω − p )Ξ , Ξ = (cid:20) ln (cid:18) (cid:37) o r ω (cid:19) − − (cid:37) i ξ ln (cid:18) (cid:37) i (cid:37) o (cid:19)(cid:21) − , (3.24)10 igure 5: Visualization of the coordinate transformation v = V ( U ( x )) = ˜ R ˜ S x . The ellipse E v is orthogonal to the well direction ψ (cid:48) which is equal to e = [0 , , T in v -coordinates. where p := p ( | w | = 0) = p ( | w | = (cid:37) i ) is the fluid pressure evaluated on thewell center-line. Note that for f = 0 and (cid:37) i = f = 0, the isotropic solutionwith for a circular constant kernel (Eqs. (2.7) and (2.10)) is obtained. From thetransformation of the Laplace operator, Eq. (3.18), we see that the problem − ∆ z p = ˆ q µρk I Φ A Φ J in Ω z = Z ( V ( U (Ω))) , (3.25)with altered kernel function Φ Λ = Φ A Φ J is equivalent to Eq. (3.22).The transformation T − changes shape of the kernel support S (Φ A ) froman annulus to an ellipse E Φ ,v . Inverting Z extrudes the solution along the wellcenter-line, and inverting the rotation and stretch described by V and U resultsin a kernel support region in the shape of an elliptic cylinder. Moreover, eachellipse E Φ ,v with normal vector e is transformed to an ellipse E Φ ,x ( s ), thatis the intersection the elliptic cylinder with a plane with the normal vector n E x = ˜ S ˜ R T e , centered at s on the well center-line. The transformation andthe normal vector n E x are visualized in Fig. 5. We note that if none of theprincipal 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 theright-hand side of Eq. (3.20) for a well segment Λ i with length L i is equal tothe integral over the kernel support S (Φ Λ ,i ) which has the shape of the ellipticcylinder given by E := (cid:91) ≤ s ≤ L i E Φ ,x ( s ) . (3.26)11sing ˆ q = qζ , and exploiting that k I was chosen such that det( ˜ R ˜ S ) = 1, it canbe shown that (cid:90) E ˆ q Φ Λ ,i d x = (cid:90) V ( U ( E )) ˆ q Φ Λ ,i d v = ˆ L i (cid:90) (cid:90) E Φ ,v (ˆ s ) ˆ q Φ Λ d ˆ A dˆ s = ˆ q ˆ L i = qL i , (3.27)where 0 ≤ ˆ s ≤ ˆ L i is a local coordinate along the transformed well direction,and the last equality is proven in Appendix C. This is the desired property ofthe kernel function for the anisotropic case corresponding to Eq. (2.2) for theisotropic case.
4. Numerical method
We discretize Eq. (3.20) using a cell-centered finite volume method withmulti-point flux approximation (MPFA) [22]. The domain Ω is decomposed intocontrol volumes K Ω ∈ Ω h such that the computational mesh Ω h is a discreterepresentation of Ω. Furthermore, each control volume boundary, ∂K Ω , can besplit into a finite number of faces σ ⊂ ∂K Ω , such that σ = K Ω ∩ L Ω , with L Ω denoting a neighboring control volume. Integrating Eq. (3.20) over a controlvolume K Ω and applying the Gauss divergence theorem on the left hand sideyields − (cid:90) ∂K Ω (cid:20) ρµ K ∇ p (cid:21) · n K Ω ,σ d A = (cid:90) K Ω ˆ q Φ Λ d x, (4.1)where n K Ω ,σ is the unit outward-pointing normal on face σ ⊂ ∂K Ω . The exactfluxes are approximated by numerical fluxes F K Ω ,σ ≈ − (cid:90) σ (cid:20) ρµ K ∇ p (cid:21) · n K Ω ,σ , (4.2)which are computed using the MPFA-O method described in [22]. The discretesource term is computed as Q K Ω ≈ (cid:90) K Ω ˆ q Φ Λ d x, Q K Ω = Q I |I| (cid:90) K Ω ∩S (Φ Λ , I ) Φ Λ d x, (4.3)where Q I is a numerical approximation of the source term integral over theintersection I = K Ω ∩ Λ, Q I = |I| π ρk I µ ( p ω − p )Ξ , (4.4)and S (Φ Λ , I ) is the kernel support associated with I as depicted in Fig. 6. Insummary, the discrete form of Eq. (4.1) is (cid:88) σ ⊂ ∂K Ω F K Ω ,σ = Q K Ω , K Ω ∈ Ω h . (4.5)12 x K Ω Λ K Ω ∩ Λ S (Φ Λ ,K Ω ∩ Λ ) ψn E x Figure 6: Visualization of the discretization process. The domain Ω is represented by a setof control volumes K Ω ∈ Ω h . The well center-line Λ with direction ψ intersects with a K Ω shown in green. The gray parallelogram is a 2D-projection of the elliptic cylinder that is thepart of the kernel support S (Φ Λ ) associated with K Ω ∩ Λ. We note that due to the dependency of Q I on p the proposed method is non-local in the sense that non-neighbor cells M Ω ∈ Ω h (where M Ω ∩ K Ω is the emptyset or a single point) may have an associated degree of freedom that dependson the degree of freedom of K Ω . The kernel integral in Eq. (4.3) I Φ ,K Ω := (cid:90) K Ω ∩S (Φ Λ , I ) Φ Λ d x, (4.6)is not easily approximated with a quadrature rule, since the intersection K Ω ∩S (Φ Λ , I ), that is the intersection of an elliptic cylinder with for example a hex-ahedron is difficult to compute. However, we use the same idea as in [15], andremark that the integral over the entire support S (Φ Λ , I ) is known exactly; seeEq. (3.27). Hence, the integration problem can be reformulated as the distri-bution of the known integral over all intersected control volumes K Ω weightedwith the respective support volume fractions. Following [15], we create n I inte-gration points x i ∈ S (Φ Λ , I ) with known volume elements V i of similar size andshape, so that I Φ ,K Ω ≈ n I (cid:88) i =1 , x i ∈ K Ω V i Φ Λ ( x i ) , n I (cid:88) i =1 V i ≈ |S (Φ Λ , I ) | . (4.7)Computing the weights for each cell K Ω is a pre-processing step that only hasto be done once for each computational mesh and well geometry.13 . Numerical experiments and discussion We present numerical experiments using the presented method in differ-ent setups. All experiments are conducted with constant fluid density ρ =1000 kg m − and viscosity µ = 1 · − Pa s. The well pressure is constant, p ω = 1 · Pa, and the well radius is r ω = 0 . K ( γ , γ ) = R ( γ ) R ( γ ) K α R T ( γ ) R T ( γ ) , K α = α · − m , (5.1)where R ( γ ) = γ − sin γ γ cos γ , R ( γ ) = cos γ γ − sin γ γ (5.2)are rotation matrices rotating vectors about e , e by the rotation angle γ , γ ,respectively, and α is a given dimensionless K-anisotropy ratio α = K K = K K .The domain Ω = [ − , × [ − , × [ − , is split in two regions,Ω = [ − , × [ − , × [0 , and Ω D = Ω \ Ω. The well center-lineΛ is given by the line through the origin and ψ = R ( β ) R ( β ) e , where R , R are given in Eq. (5.2) and β , β , are rotation angles. The analytical solutionfor all cases is given in Eq. (3.23), q = 1 kg s − m − , and L = | Λ ∩ Ω | (in m).For all setups the inner kernel radius is chosen as (cid:37) i = f . In all of Ω D andon the boundary ∂ Ω the analytical solution is enforced by Dirichlet constraints,modeling the infinite well. The computational mesh Ω h is a structured gridcomposed of regular hexahedra K Ω . Furthermore, we define two error measures. E p = 1 p ω (cid:34) | Ω h | (cid:88) K Ω ∈ Ω h | K Ω | (cid:16) p e, x K Ω − p K Ω (cid:17) (cid:35) (5.3)is the relative discrete L -norm of the pressure, where p e, x K Ω is the exact pres-sure evaluated at the cell centroid x K Ω and p K Ω the discrete numerical cellpressure, and E q = 1 q | Λ ∩ Ω h | (cid:88) K Ω ∈ Ω h K Ω ∩ Λ (cid:54) = ∅ |I| (cid:18) q − Q I |I| ζ (cid:19) (5.4)is the relative discrete L -norm of the source term, where |I| = | K Ω ∩ Λ | is thelength of the intersection of cell K Ω and the well center-line Λ, Q I is the discretesource term given in Eq. (4.4). All setups are implemented in DuMu x [23], anopen-source porous media simulator based on Dune [24, 25].14 h max in m E p = 1= 10= 50= 100( h ) 10 346 h max in m E q = 1= 10= 50= 100( h ) Figure 7: Grid convergence for the relative discrete L -norm pressure E p and source E q fordifferent anisotropy ratios α . None of the principal axis of the permeability tensor is alignedwith the slanted well axis ψ or any of the grid axes. h max α .
32 m 8 .
66 m 4 .
33 m 2 .
17 m1 2.0545 2.0724 1.9454 -10 1.7715 2.0184 2.0763 -50 1.5904 1.9747 2.0925 -100 1.5970 1.9666 2.1218 -
Table 1: Convergence rates for E q for different anisotropy ratios α . In the first numerical experiment grid convergence is investigated for differentanisotropy ratios α . To this end, h max := max K Ω ∈ Ω h h K Ω , where h K Ω is definedas the maximum distance between two vertices of the cell K Ω . Starting at agrid resolution for Ω h of 20 × ×
10 cells ( h max = 10 √ E p and E q for different grid resolutionsand values of α , for β = β = 20 ◦ and γ = γ = − ◦ , so that K is a fulltensor 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 thegiven norm, as expected for the MPFA-O method [26] (super convergence atcell centers). The source term q is a linear function of the pressure p and alsoexhibits second order convergence. We note that the errors for different α are notdirectly 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 thelines in Fig. 7) are presented in Table 1. It can be seen that rates for large grid15
00 300 400 600 o / r E q for different o ( r = 0.1m) Simplified J Exact J ( h ) 2 3 4 6 o / r E q for different o ( r = 10.0m) Simplified J Exact J ( h ) Figure 8: The source error E q for different kernel supports and the same 20 × ×
10 compu-tational grid ( h max ≈ .
32 m). On the left, the case where (cid:37) o is only slightly larger than r ω .On the right, the case (cid:37) o (cid:29) r ω . cells and large α are slightly smaller. This is because the kernel support is stillunder-resolved by the computational grid. For example, for α = 100, the kernelellipse in x -coordinates has major and minor axis of a x ≈ . b x ≈ . h max ≈ .
32 m for the lowest grid resolution. (cid:37) o In [15], it is suggested that increasing the kernel support region (increasing (cid:37) o ), has a similar effect on E q as refining the grid. However, the pressure solutionis then regularized in a larger region, so that there is a trade-off between theaccuracy of the source term and the accuracy of the pressure field with respectto the unmodified problem ( (cid:37) o → (cid:37) i , Φ Λ → δ Λ ). However, every discrete cell K Ω can be also interpreted as a kernel support region, such that the choiceof Φ Λ enables us to better control the discretization error as soon as S (Φ Λ , I )becomes larger than K Ω .As shown in [6] and Fig. 4, isobars become circular, in the transformeddomain U (Ω), with increasing distance to the well. Therefore, a reasonablesimplification is Φ J ≈ (cid:37) o (cid:29) r ω . This is completely analogous to the as-sumption of circular isobars in [6], where an estimate of the error introduced bythe assumptions is given for the two-dimensional case.In the following numerical experiment, we step-wise increase the kernel ra-dius (cid:37) o , for the same 20 × ×
10 grid. This is done once for the case, where (cid:37) o (cid:29) r ω and for the case for which (cid:37) o is only slightly larger than r ω . Further-more, β = β = 20 ◦ and γ = γ = − ◦ . The results are shown in Fig. 8.First, it can be seen that doubling (cid:37) o leads to a 4-times smaller error E q . Thiscan be explained by the fact that the larger the kernel, the more grid cells resolvethe kernel support, and the better is the approximation of p . Furthermore, theresult is consistent with the results in [15]. Moreover, Fig. 8 suggests that for (cid:37) o (cid:29) r ω the simplification of the kernel function (Φ J ≈
4) is not visible in E q ,16 in degree E q for different ( K angle) = 60 = 40 = 20 =0 =20 =40 =60 30 20 10 0 10 20 30 in degree E q for different (well angle) = 30 = 20 = 10 =0 =10 =20 =30 Figure 9: The source error E q for rotations of the permeability tensor (left) and different wellorientations (right). while for kernel radii slightly smaller than the well radius, the simplification in-creases E q by an order of magnitude in comparison to the case using the exactkernel function as derived in Section 3.3. The results show that the presentedmethod is also applicable in cases where the grid resolution is very close to thewell radius. An adaption of the presented method for other applications, suchas the simulation of flow in vascularized tissue, where such ratios of vessel radiusto cell size are typical, cf. [15], is therefore well-conceivable. In the following numerical experiment, we use a single computational meshwith a given resolution for Ω h : 20 × ×
10. First, the well direction is fixed, andthe permeability tensor is rotated by varying γ and γ . Then the permeabilitytensor is fixed and the well is rotated by varying β and β . The results areshown in Fig. 9. It can be seen that the presented well model is rather robustwith respect to rotations. Possible effects influencing the approximation error E q , include the different quality of the kernel integral for different angles withrespect to the grid axes, and differences in the flux approximation quality ofthe MPFA-O method depending on the face co-normal d K Ω ,σ = K n K Ω ,σ . Ad-ditionally, for different well angles the number and size of intersections K Ω ∩ Λcan have an influence on the discrete error.
In particular for petroleum engineering applications, commercial codes typ-ically use Peaceman-type well models [27, 28]. In [6], Peaceman extended hiswell known well-index-based well model for anisotropic diagonal permeabilitytensors and non-cubic but structured rectangular grids. The discrete source17erm in a computational cell K Ω is approximated by Q K Ω = 2 π ρµ ( p ω − p ) L K Ω √ K K ln (cid:16) r r ω (cid:17) , (5.5) r = e − γ (cid:20)(cid:16) K K (cid:17) ∆ x + (cid:16) K K (cid:17) ∆ y (cid:21) (cid:16) K K (cid:17) + (cid:16) K K (cid:17) , (5.6)where ∆ x and ∆ y are the horizontal dimensions of the cell containing the well, L K Ω = |I| the length of the well segment contained in K Ω , and γ the Euler–Mascheroni constant. The Peaceman model has several known limitations. Itsderivation only applies to K-orthogonal structured grids, where the well is ori-ented along one of the grid axes, and perfectly horizontally centered within avertical column of computational cells K Ω . Furthermore, the derivation is spe-cific 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 Peacemanmodel has been generalized for slanted wells with arbitrary orientation, for ex-ample in [8]. The Alvestad model [8] has been adapted for finite volumes, forexample in [9] (formula given in Appendix E, subsequently referred to as pm wellmodel). Such extensions usually constitute a reasonable directional weighting ofthe original Peaceman model but are not directly derived from the mathematicalanalysis of the underlying problem [9].The herein presented model has none of the above-mentioned limitations.In particular, the presented model is valid for arbitrary positive definite andsymmetric permeability tensors, unstructured grids, and is independent of thediscretization scheme. Moreover, the presented model is consistent and we showgrid convergence in the numerical experiments in Section 5.1. However admit-tedly, the Peaceman-type models are cell-local, thus computationally cheaperand easier to implement.Several limitations of the Peaceman well models make it difficult to fairlycompare it with our new model. For cases for which all assumptions of Peacemanare valid, our numerical studies (not shown here) suggest that the Peacemanwell model is generally superior to the presented model with distributed sources.This is because it takes the analytical solution as well as the spatial discretiza-tion method into account. For cases where some assumptions are violated, forexample off-center wells or slanted wells, it is difficult to construct cases wherethe analytical solution is readily constructed but does not feature a singularityon the boundary. Our preliminary numerical studies for such cases (for examplethe slanted well case in Section 5.1 without rotation of the permeability tensor)show large deviations ( >
10 % error in total source term) from the analyticalsolution for the pm well model. However, these errors may be distorted by errorsmade in the discrete approximation of the singular boundary condition, wherethe well intersects the boundary. Finally, for the general case of unstructured18 igure 10: The computational domain for the comparison with a Peaceman-type well model.The well is visualized with a 10-fold increased radius. A selection of pressure iso-surfaces ofthe reference solution are shown with reduced opacity. The domain extent is given in units ofm. grids, simplex grids, and full permeability tenors it is unclear how to applythe original Peaceman model. However, we know that the presented method isconsistent (at least for a single straight well), and thus, the numerical solutionconverges to the exact solution with grid refinement. Therefore, we expect thatthe numerical solution on a very fine grid using the distributed source model isa reasonable reference solution.We compare our model to the pm well model in a numerical experiment.The computational domain Ω = [ − , × [ − , × [0 , containsa slanted straight well Λ with end points at x Λ , = [ − , − , T m, x Λ , =[20 , , T m. The permeability tensor is a diagonal tensor K ( γ , γ ), with γ = 0 ◦ , γ = 90 ◦ , α = 0 .
1. The structured cube grid Ω h is successively, uni-formly refined starting with 10 × ×
10 cells ( h max ≈ .
32 m). The well radiusis r ω = 0 . x/r ω = 200 for the coarsest grid). The kernel support region(chosen as (cid:37) o /r ω = 100) only extends over few cells in the coarsest grid, sothat the regularization effect is minimized. On the boundary ∂ Ω, we specifyNeumann no-flow boundary conditions, that is ( K ∇ p ) · n = 0, except for theplanes perpendicular to the x -axis, where the Dirichlet boundary condition p D ( x = − · Pa, p D ( x = 100) = 3 · Pa are enforced. The ref-erence solution is computed with 160 × ×
160 cells ( h max ≈ .
08 m). Thecomputational domain with pressure iso-surfaces of the reference solution areshown in Fig. 10. In Fig. 11, the relative integral source error E Q = | Q − Q ref || Q ref | , Q = (cid:88) K Ω ∈ Ω h Q K Ω , (5.7)19 h max in m | QQ r e f |/| Q r e f | o )PMDS (const o / x ) Figure 11: Comparison of the relative integral source error between a Peaceman-type model( pm ) and the new model ( ds ) for various grid refinements. The error is computed with respectto a reference solution Q ref . Both axes are logarithmic. with respect to the reference solution Q ref is shown for grids with different re-finement. In a variant of the distributed source model ( ds ), the extent of thekernel support is adapted to the grid size. This is to keep the regularizationeffect of the kernel function minimal in order to get, in addition to a good ap-proximation of the source term, a better approximation of the pressure solutionclose to the well. For (cid:37) o /r ω = 100, the extent of the kernel ellipse E Φ ,x is givenby its major and minor axes, 16 .
12 m and 12 .
54 m. For the reference solutionthis extent is kept constant with grid refinement. While this ensures a very goodapproximation of the source term, the pressure solution is regularized in a largerneighborhood of the well. In the variant, the kernel support is adapted propor-tional to h max , so that for the finest grid shown in Fig. 11 (80 × ×
80 cells, h max ≈ .
17 m), the E Φ ,x major and minor axes measure 2 .
01 m and 1 .
57 m. Itis evident that the numerical solution for the distributed source model convergesto the reference solution. More importantly, the relative error is small ( < . > > r ω / ∆ x ratio. In the variant of the ds model, theerror in the source term with respect to the reference solution also grows withlarger r ω / ∆ x ratio. However, the error is consistently smaller (by a factor > pm well model. Figure 12 shows the numerical pressure solutions20 x -coordinate p i n P a DS (const o )PMDS (const o / x ) 100 50 0 50 100 x -coordinate o )PMDS (const o / x ) Figure 12: Numerical pressure solutions plotted along the x and the x -axis for a Peaceman-type model ( pm ) and the new model ( ds ) for a grid resolution of 160 × ×
160 cells. along the x and the x -axis, for the reference grid resolution. It can be clearlyseen that for the reference solution ( ds ) the pressure solution is regularized. Forthe variant of ds , the regularization is minimized, however in the far field thesolution matches the reference solution better than the pm well model, whichis due to the better approximation of the source term (see Fig. 11). We alsonote that the regularized solution leads to an altered solution in the near-fieldof the well but to a better approximation of the source term and thus the farfield pressure (outside the kernel support), whereas the poor approximation ofthe source term in the pm method leads to a globally poor pressure solution.
6. Summary
A new well model was presented for which the mass exchange between a welland an embedding porous medium is modeled with a source term spatially dis-tributed by a local kernel function. In the spirit of well-index-based well modelsthe source term for a well with given bottom hole pressure is computed basedon the numerical pressure in cells intersecting the well. However, the presentedderivation of the new model is independent of the discretization method andthe type of computational grid. The new model was shown to be consistent in anumerical experiment and exhibited grid convergence with the expected rates.In the same experiment it is shown that the absolute error with respect to ananalytical solution is relatively small, even for coarse computational grids andsmall kernel support. It was shown, that the error in the source term can bedecreased by increasing the region over which the source term is distributed.However, coincidentally, the pressure profile close to the well (inside the kernelsupport) becomes increasingly regularized. A comparison with a Peaceman-type well model generalized for arbitrarily oriented wells, suggested that even21f the region is chosen to be very small (only covering the neighboring cells ofcells with well intersection), thus minimizing the regularization effect, the sourceterm can be approximated with good accuracy ( < > 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 Pro-cesses in Porous Media (SFB 1313, Project Number 327154368).
Appendix A Rodrigues’ rotation formula
We want to rotate a given basis B = { e , e , e } such that e is aligned witha vector ψ . The Rodrigues’ rotation formula [29] describes a rotated vector x rot obtained by rotating a vector x by the angle θ about an axis given by the unitnormal vector kx rot = x cos θ + ( k × x ) sin θ + k ( k · x )(1 − cos θ ) . (A.1)The desired rotation can be described by a rotation by an angle θ = π aboutthe axis k = ( e + ψ )2 | e + ψ | , x rot ,π = 2 k ( k · x ) − x . (A.2)22quation (A.2) can be expressed in matrix notation as x rot ,π = R x with R = (cid:16) kk T − I (cid:17) , R = R T = R − , det( R ) = 1 . (A.3) Appendix B Transformation of Laplace operator
The Joukowsky transformation is a complex function z = T − ( w ) that canbe decomposed in its real and imaginary parts, x = (cid:60) ( z ) and y = (cid:61) ( z ). Fur-thermore, let u = (cid:60) ( w ) and v = (cid:61) ( w ). As a conformal mapping, z satisfies theCauchy Riemann equations [19] ∂x∂u = ∂y∂v and ∂x∂v = − ∂y∂u . (B.1)The Jacobian of the transformation J T − is given by Eq. (3.16). We investigationthe effect of the transformation on the Laplace operator∆ w p = ∂ p∂u + ∂ p∂v , (B.2)where p is analytic in Ω w . Applying the chain rule yields ∂p∂u = ∂p∂x ∂x∂u + ∂p∂y ∂y∂u , (B.3) ∂ p∂u = ∂p∂x ∂ x∂u + ∂p∂y ∂ y∂u + ∂∂x (cid:18) ∂p∂u (cid:19) ∂x∂u + ∂∂y (cid:18) ∂p∂u (cid:19) ∂y∂u (B.4)= ∂p∂x ∂ x∂u + ∂p∂y ∂ y∂u + ∂ p∂x (cid:18) ∂x∂u (cid:19) + 2 ∂ p∂x∂y ∂x∂u ∂y∂u + ∂ p∂y (cid:18) ∂y∂u (cid:19) . Analogously, we arrive at a similar expression for ∂ p/∂v . Using Eq. (B.1) and ∂ y∂v = ∂∂v (cid:18) ∂y∂v (cid:19) (B.1) = ∂∂v (cid:18) ∂x∂u (cid:19) = ∂∂u (cid:18) ∂x∂v (cid:19) (B.1) = − ∂y ∂u , ∂ x∂v = − ∂ x∂u , (B.5)we find that ∂ p∂u + ∂ p∂v = (cid:34)(cid:18) ∂x∂u (cid:19) + (cid:18) ∂y∂u (cid:19) (cid:35) (cid:20) ∂ p∂x + ∂ p∂y (cid:21) . (B.6)With the complex derivative of T − [21], ∂T − ∂w = ∂T − ∂u = ∂x∂u + i ∂y∂u and (cid:12)(cid:12)(cid:12)(cid:12) ∂T − ∂w (cid:12)(cid:12)(cid:12)(cid:12) = (cid:115)(cid:18) ∂x∂u (cid:19) + (cid:18) ∂y∂u (cid:19) . (B.7)From the determinant of the Jacobian of the transformation, we finddet( J T − ) = ∂x∂u ∂y∂v − ∂x∂u ∂y∂v = (cid:18) ∂x∂u (cid:19) + (cid:18) ∂y∂u (cid:19) = (cid:12)(cid:12)(cid:12)(cid:12) ∂T − ∂w (cid:12)(cid:12)(cid:12)(cid:12) , (B.8)23sing Eq. (B.1). Hence,∆ w p = (cid:12)(cid:12)(cid:12)(cid:12) ∂T − ∂w (cid:12)(cid:12)(cid:12)(cid:12) ∆ z p = | det( J T − ) | ∆ z p, (B.9)which also proofs that any harmonic function (∆ f = 0) yields another harmonicfunction after a coordinate transformation with a conformal mapping. (cid:3) Appendix C Source scaling factor in w -coordinates We want to construct a pressure solution in w -coordinates such that the totalmass flux over the well boundary matches the specified boundary condition in x -coordinates. Hence, the total mass flux over the boundary of a well segmentwith length ˆ L in w -coordinates needs to match the total mass flux over theboundary of a well segment with length L in x -coordinates. A ˆ q has to bechosen such that qL = ˆ q ˆ L . A relation between L and ˆ L can be derived bylooking at two related volume integrals. The Joukowsky transformation onlyaffects the two-dimensional well-bore plane such that length ˆ L of a well segmentis not affected. The volume of that well segment in v -coordinates (an ellipticcylinder) is given by V v = πab ˆ L. (C.1)The volume of the same well segment in x -coordinates (a circular cylinder withslanted parallel planar elliptic caps) is given by V x = | ψ T ˜ S ˜ R T e || E ω,x | L = πr ω L, (C.2)where | E ω,x | is the area of the well-bore ellipse described by Eq. (3.6) trans-formed to x -coordinates (as shown in Fig. 5) and the last equality uses thefact that the integral can be transformed to an integral over a regular cylinderwith 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 = L r ω ab , (C.3)and if the source term is chosen asˆ q = q abr ω := qζ, (C.4)then qL = ˆ q ˆ L . (cid:3) Appendix D Determinant of the Joukowsky transformation
In Appendix B, we show that | det( J T − ) | = (cid:12)(cid:12) ∂z∂w (cid:12)(cid:12) . Using complex differen-tiation, ∂z∂w = ∂z∂u = ∂∂u (cid:20) (cid:18) w + f w (cid:19)(cid:21) = 12 (cid:18) − f w (cid:19) = 12 (cid:18) − f w | w | (cid:19) , (D.1)24here we used the identities w − = w | w | − , w denoting the complex conjugateof w , and w = w . Furthermore, (cid:60) (cid:18) ∂z∂w (cid:19) = 14 (cid:32) − f (cid:60) (cid:0) w (cid:1) | w | + f (cid:60) (cid:0) w (cid:1) | w | (cid:33) , (D.2) (cid:61) (cid:18) ∂z∂w (cid:19) = 14 (cid:32) f (cid:61) (cid:0) w (cid:1) | w | (cid:33) , (cid:60) (cid:0) w (cid:1) + (cid:61) (cid:0) w (cid:1) = | w | , (D.3)such that (cid:12)(cid:12)(cid:12)(cid:12) ∂z∂w (cid:12)(cid:12)(cid:12)(cid:12) = (cid:60) (cid:18) ∂z∂w (cid:19) + (cid:61) (cid:18) ∂z∂w (cid:19) = 14 (cid:18) f − f (cid:60) ( w ) | w | (cid:19) . (D.4) (cid:3) Appendix E Extension of Peaceman well model for slanted wells
For a given well direction ψ = [ ψ , ψ , ψ ] T and a cell K Ω with dimensions∆ x × ∆ y × ∆ z , to obtain the generalized well model due to [8], reformulated forcell-centered finite volume schemes in [9], replace k = √ K K in Eq. (5.5) by k = ( ψ K K + ψ K K + ψ K K ) (E.1)and the expression for r by r = e − γ (cid:0) ∆ L + ∆ L (cid:1) (cid:0) √ A + √ A (cid:1) , with (E.2)∆ L = (cid:18) K K (cid:19) ∆ z ψ + (cid:18) K K (cid:19) ∆ x ψ + (cid:18) K K (cid:19) ∆ y ψ , (E.3)∆ L = (cid:18) K K (cid:19) ∆ y ψ + (cid:18) K K (cid:19) ∆ z ψ + (cid:18) K K (cid:19) ∆ x ψ , (E.4) A = (cid:18) K K (cid:19) ψ + (cid:18) K K (cid:19) ψ + (cid:18) K K (cid:19) ψ , (E.5) A = (cid:18) K K (cid:19) ψ + (cid:18) K K (cid:19) ψ + (cid:18) K K (cid:19) ψ . (E.6)We note that the formula reduces to Eq. (5.5), if ψ is aligned with a coordinateaxis. References [1] T. Koch, B. Flemisch, R. Helmig, R. Wiest, D. Obrist, A multi-scale sub-voxel perfusion model to estimate diffusive capillary wall conductivity inmultiple sclerosis lesions from perfusion mri data, bioRxiv (2018). doi:10.1101/507103 . 252] L. Cattaneo, P. Zunino, Computational models for fluid exchange betweenmicrocirculation and tissue interstitium, Networks & Heterogeneous Media9 (1) (2014). doi:10.3934/nhm.2014.9.135 .[3] T. Koch, K. Heck, N. Schr¨oder, H. Class, R. Helmig, A new simulationframework for soil-root interaction, evaporation, root growth, and solutetransport, Vadose Zone Journal 17, 1 (2018). doi:10.2136/vzj2017.12.0210 .[4] C. Doussan, L. Pages, G. Vercambre, Modelling of the Hydraulic Ar-chitecture of Root Systems: An Integrated Approach to Water Absorp-tion—Model Description, Annals of Botany 81 (2) (1998) 213–223 (1998). doi:10.1006/anbo.1997.0540 .[5] D. W. Peaceman, Interpretation of well-block pressures in numerical reser-voir simulation, Society of Petroleum Engineers Journal 18 (03) (1978)183–194 (1978).[6] D. W. Peaceman, Interpretation of well-block pressures in numerical reser-voir simulation with nonsquare grid blocks and anisotropic permeability,Society of Petroleum Engineers Journal 23 (03) (1983) 531–543 (1983). doi:10.2118/10528-PA .[7] Z. Chen, Y. Zhang, Well flow models for various numerical methods., In-ternational Journal of Numerical Analysis & Modeling 6 (3) (2009).[8] J. Alvestad, K. Holing, K. Christoffersen, O. Stava, SPE-27577-MS, Societyof Petroleum Engineers, Aberdeen, United Kingdom, 1994, Ch. InteractiveModelling of Multiphase Inflow Performance of Horizontal and Highly De-viated Wells, p. 16 (1994). doi:10.2118/27577-MS .[9] I. Aavatsmark, R. A. Klausen, Well index in reservoir simulation for slantedand slightly curved wells in 3d grids, SPE Journal 8 (01) (2003) 41–48(2003). doi:10.2118/75275-PA .[10] D. K. Babu, A. S. Odeh, Productivity of a horizontal well, SPE ReservoirEngineering 4 (04) (1989) 417–421 (1989). doi:10.2118/18298-PA .[11] D. K. Babu, A. S. Odeh, A. J. Al-Khalifa, R. C. McCann, The relationbetween wellblock and wellbore pressures in numerical simulation of hor-izontal wells, SPE Reservoir Engineering 6 (03) (1991) 324–328 (1991). doi:10.2118/20161-PA .[12] C. Wolfsteiner, L. J. Durlofsky, K. Aziz, Calculation of well index for non-conventional wells on arbitrary grids, Computational Geosciences 7 (1)(2003) 61–82 (Mar 2003). doi:10.1023/A:1022431729275 .[13] H. B. Hales, SPE-39065-MS, Society of Petroleum Engineers, Rio deJaneiro, Brazil, 1997, Ch. An Improved Method for Simulating ReservoirPressures Through the Incorporation of Analytical Well Functions, p. 6(1997). doi:10.2118/39065-MS . 2614] I. G. Gjerde, K. Kumar, J. M. Nordbotten, B. Wohlmuth, Splittingmethod for elliptic equations with line sources, arXiv e-prints (2018)arXiv:1810.12979 (Oct 2018). arXiv:1810.12979 .[15] T. Koch, M. Schneider, R. Helmig, P. Jenny, Modeling tissue perfu-sion in terms of 1d-3d embedded mixed-dimension coupled problems withdistributed sources, arXiv e-prints (2019) arXiv:1905.03346 (May 2019). arXiv:1905.03346 .[16] C. R. Fitts, Exact solution for two-dimensional flow to a well in ananisotropic domain, Groundwater 44 (1) (2006) 99–101 (2006). doi:10.1111/j.1745-6584.2005.00082.x .[17] I. Aavatsmark, Interpretation of well-cell pressures on hexagonal grids innumerical reservoir simulation, Computational Geosciences 20 (5) (2016)1029–1042 (Oct 2016). doi:10.1007/s10596-016-9575-2 .[18] J. Bear, G. Dagan, The relationship between solutions of flow problems inisotropic and anisotropic soils, Journal of Hydrology 3 (2) (1965) 88 – 96(1965). doi:10.1016/0022-1694(65)90002-8 .[19] Z. Nehari, Conformal mapping, Dover Publications, Inc., New York, 1975,reprinting of the 1952 edition (1975).[20] N. Joukowsky, ¨Uber die Konturen der Tragfl¨achen der Drachenflieger,Zeitschrift f¨ur Flugtechnik und Motorluftschiffahrt 1 (1910) 281–284 (1910).[21] W. Rudin, Real and Complex Analysis, 3rd Ed., McGraw-Hill, Inc., NewYork, NY, USA, 1987 (1987).[22] I. Aavatsmark, An introduction to multipoint flux approximations forquadrilateral grids, Computational Geosciences 6 (3) (2002) 405–432 (Sep2002). doi:10.1023/A:1021291114475 .[23] B. Flemisch, M. Darcis, K. Erbertseder, B. Faigle, A. Lauser, K. Mosthaf,S. M¨uthing, P. Nuske, A. Tatomir, M. Wolff, R. Helmig, DuMu x : DUNE formulti- { phase, component, scale, physics,. . . } flow and transport in porousmedia, Advances in Water Resources 34 (9) (2011) 1102–1112 (2011). doi:10.1016/j.advwatres.2011.03.007 .[24] P. Bastian, M. Blatt, A. Dedner, C. Engwer, R. Kl¨ofkorn, M. Ohlberger,O. Sander, A Generic Grid Interface for Parallel and Adaptive ScientificComputing. Part I: Abstract Framework, Computing 82 (2–3) (2008) 103–119 (2008). doi:10.1007/s00607-008-0003-x .[25] P. Bastian, M. Blatt, A. Dedner, C. Engwer, R. Kl¨ofkorn, R. Korn-huber, M. Ohlberger, O. Sander, A Generic Grid Interface for Paralleland Adaptive Scientific Computing. Part II: Implementation and Testsin DUNE, Computing 82 (2–3) (2008) 121–138 (2008). doi:10.1007/s00607-008-0004-9 . 2726] M. Schneider, D. Gl¨aser, B. Flemisch, R. Helmig, Comparison of finite-volume schemes for diffusion problems, Oil Gas Sci. Technol. - Rev. IFPEnergies nouvelles 73 (2018) 82 (2018). doi:10.2516/ogst/2018064 .[27] Schlumberger N.V., ECLIPSE Technical Description, 2014.[28] Computer Modelling Group Ltd., IMEX User’s Guide, 2014.[29] J. S. Dai, Euler–rodrigues formula variations, quaternion conjugation andintrinsic connections, Mechanism and Machine Theory 92 (2015) 144 – 152(2015). doi:10.1016/j.mechmachtheory.2015.03.004doi:10.1016/j.mechmachtheory.2015.03.004