# Fundamentals in Peaceman Model for Well-Block radius For Non-Linear Flows
Near Well

A. Ibraguimov1, E. Zakirov.2, I. Indrupskiy3, D. Anikeev4

###### Abstract.

We consider sewing machinery between finite difference and analytical
solutions defined at different scale: far away and near source of the
perturbation of the flow. One of the essences of the approach is that coarse
problem and boundary value problem in the proxy of the source model two
different flows. We are proposing method to glue solution via total fluxes,
which is predefined on coarse grid. It is important to mention that the coarse
solution "does not see" boundary.

From industrial point of view our report provide mathematical tool for
analytical interpretation of simulated data for fluid flow around a well in a
porous medium. It can be considered as a mathematical "shirt" on famous
Peaceman well-block radius formula for linear (Darcy) radial flow but can be
applied in much more general scenario.

As an important case, we consider nonlinear Forchheimer flow. In the article
we rigorously obtained well-block radius, explicitly depending on
$\beta-$Forchheimer factor and total rate of the flow on the well, and provide
generalization of the Dake Formula and evaluation of the $D-$factor.

1 Department of Mathematics and Statistics, Texas Tech University, Oil and Gas
Research Institute of Russian Academy of Sciences

Lubbock, Texas, USA, Moscow,Russia e-mail<EMAIL_ADDRESS>
2 Oil and Gas Research Institute of Russian Academy of Sciences,
Moscow, Russia, e-mail<EMAIL_ADDRESS>
3 Oil and Gas Research Institute of Russian Academy of Sciences,
Moscow, Russia, e-mail<EMAIL_ADDRESS>
4 Oil and Gas Research Institute of Russian Academy of Sciences,
Moscow, Russia, e-mail<EMAIL_ADDRESS>

###### Contents

  1. 1 Introduction
  2. 2 Nomenclature
  3. 3 About well representation in a 3D flow model
  4. 4 Brief review of the connection transmissibility computation problem
  5. 5 Peaceman Well Block Radius and Fundamentals in Finite Difference Solution
    1. 5.1 Origin of the Peaceman Formula Revisited
    2. 5.2 Green function and related interpretation of the well-block radius
    3. 5.3 Radial Darcy flow in two embedded annulus zones sewed by linear balance equation on the finite difference coarse grid, and its link to capacity
    4. 5.4 Two Terms Radial Case
    5. 5.5 Criteria for radial non-linear Forchheimer flow in two annular zones with linear balance equation on the finite difference coarse grid
  6. 6 Impact of the Correction in the Dake Formula
    1. 6.1 Dake Formula and D-factor
  7. 7 Conclusive Remarks and Discussion

##  1\. Introduction

Many industrial simulators of the processes of fluid flows are based on the
numerical solution of the partial differential equations (PDE) models (see for
example [41],[42],[44],[43]). Due to large difference in scales, there is a
need for analytical approximation to replace numerical solution near the
source (boundary) controlling the processes.

Aim of the article is to revisit Peaceman well block radius, which is
routinely used by reservoir engineers to link value of the pressure for the
grid block of numerical solution to the actual value of the pressure on the
well. In the paper we investigate this issue from mathematical point for
linear Darcy flow and non-linear Forchheimer type of the flow.

We reformulate Peaceman method to make analytical and mathematical aspects of
the problem transparent. In this paragraph, we will highlight main topics of
our interest in this project. Consider baseline case of the 2-D steady regime
of the flow initiated in the fully penetrated reservoir of the thickness $h=1$
by sole vertical well of radius $r_{w}$, with given production rate $q$.
Reservoir pressure on the external boundary will be considered fixed. Let
generated grid in the domain to be such that source (final) block $Q_{0}$ of
size $\Delta$ contains the well. Consequently, RHS of the corresponding system
of the algebraic equations is homogeneous (null) in all cells but $Q_{0}$.

Let numerical solution of the PDE to be $\tilde{p}$ defined at each block on
the grid takes value $p_{0}$ at $Q_{0}$. By construction $\tilde{p}$ forms a
finite dimensional matrix and is defined at discreet "points". It does not
"see" the well of the small radius $r_{w}$. In the same time well is
generating flows in the domain, and it is vital to use $p_{0}$ for finding
$p_{w}$, which one can assign as pressure value on the well for further
processing. To do so Peaceman considered material balance (or finite-
difference) equation in the five-spot system of the grid blocks:
$\\{1\\},\\{2\\},\\{3\\},\\{4\\}$, and $\\{0\\}$ (see Fig. 1) [14]:

$\displaystyle\frac{k}{\mu}\cdot\frac{(p_{1}-2p_{0}+p_{3})}{(\Delta x)^{2}})+$
(1.1) $\displaystyle+\frac{k}{\mu}\cdot\frac{(p_{2}-2p_{0}+p_{4})}{(\Delta
y)^{2}}=\frac{q}{\Delta x\cdot\Delta y},$

where $p_{1}=p(-\Delta x,0),\ p_{3}=p(\Delta x,0)\,\ p_{2}=p(0,\Delta y),\
p_{4}=p(0,-\Delta y),\ \text{and }p_{0}=p(0,0).$

Assuming symmetry constraint $p_{1}=p_{2}=p_{3}=p_{4}$, and
$\Delta_{x}=\Delta_{y}=\Delta$, Peaceman then reduced (1) to the equation

(1.2) $p_{1}-p_{0}=q\frac{\alpha}{4},\ \alpha=\frac{\mu}{k}.$

To relate value $p_{0}$ to the value $p_{w}$ on the well Peaceman uses Darcy-
Dupuit formula, which relates pressure drop $p_{\theta}-p_{r}$ between value
of the pressure $p_{\theta}=p(\theta),$ on the imaginary well of radius
$\theta$ and current pressure $p_{r}=p(r)$ at radius $\theta<r\leq\Delta$, to
the total rate $q:$

(1.3) ${-}\alpha^{-1}\frac{\partial p(r)}{\partial r}=v_{r}\leftrightarrow
p_{r}=p_{\theta}+q\alpha\cdot\frac{1}{2\pi}\ln\frac{r}{\theta}.$

Main question is as follows: For given $\Delta$, does auxiliary $R_{0}$ exists
such that $p_{R_{0}}-p_{\Delta}=p_{0}-p_{1}$ for any $q$?

In the above $p_{R_{0}},\ p_{\theta},p_{\Delta}$ are obtained by Dupuit-Darcy
equation of flow for $r=R_{0}$, $r=\theta$ and $r=\Delta$ respectively.

Peaceman obtained the remarkable answer [14]: Such $R_{0}$ exists and is
defined via equation $R_{0}=\Delta\cdot e^{-\pi/2}.$ And this $R_{0}$ is
called the Peaceman well block radius.

Note that Peaceman’s $R_{0}$ does not depends on $q$ and $\theta.$

Figure 1. Five-spot grid of size $\Delta$ with the well at center, and
auxiliary $U(0,R_{w},R_{1})$ and $U(0,R_{w},R_{0})$

Then pressure on the well can be reasonably approximated from any numerical
solution (value $p_{0}$) as shown in Sec. 5.1.

This issue is discussed in detail in Sec. 5.3, in which we also introduce
Peaceman well-posedness in order to make arguments more rigour, which can be
used in different scenarios.

Mathematically material balance equation (1.2) and Dupuit-Darcy equation (1.3)
of flow are split. In the article we further consider non-linear flows in
porous media under hypothesis that block-to-block (global) material balance
can be assumed linear, and Forchheimer non-linearity is included only in the
equation of flow (1.3) and has a form: $-\frac{\partial p(r)}{\partial
r}=\alpha v_{r}+\beta v_{r}^{2}$ . Considering the same arguments as for the
linear Darcy flow, this results in a dependency of $R_{0}$ on $q$. We
implemented this algorithm in Sections 5.5 and 6.1, and obtained formulae for
$R_{0}$ which is very close to traditional one used in simulators if $\Delta$
is large, but quite different when $\Delta$ is small.

In our up-coming research projects we will extend this approach to different
classes of pre-Darcy and post-Darcy equations, based on analytical
formulae(see [9]) for steady state and PSS solutions for generalised non-Darcy
equation. We will also investigate non-linear Forchheimer-type balance
equation of the form $4(p_{1}-p_{0})=\alpha q+f(q),$ for non-linear flow.

Note that domain of the flow $U=\Omega\setminus B(0,r_{w})$ does not contain
singular point, and method of the domain discretization based on the material
balance is such that $\Delta>r_{w}.$ This does not allow unbounded grid
refinement, consequently theoretical comparison between analytical solution
and numerical one is limited. An interesting mathematical insight for Peaceman
formula came from comparison of the numerical solution with Green function
which we did in Sec. 5.2. Namely, finite difference scheme applied for Laplace
equation defined in whole domain with isolated source point generates
numerical solution which is an approximation to the Green function in the
"big" domain $\Omega$ with singularity at $0$ . Then analogous to the Peaceman
well posedness, the question for the analytical Green function $G(x,y)$ is
stated as follows:  Does auxiliary $R_{0}$ exists such that
$\left(G(x,y)\Big{|}_{|X|=\Delta}-G(x,y)\Big{|}_{|X|=R_{0}}\right)\asymp
1+O(\Delta)?$ In Theorem 3 we proved that the answer is $R_{0}=\Delta\cdot
e^{-\pi/2}.$ Which once more conforms Peaceman formula for $\Delta\asymp 0,$
as well.

To investigate non-linear flows we propose splitting method in which balance
equation is the same as for Darcy flow, but non-linearity impact on the
relation between $p_{0}$ and $p_{1}$ and value of the pressure on the well is
considered through Forchheimer law. This results in the dependence of the
auxiliary $R_{0}$ on $q$ (see Section 5.5). Obtained relation justifies
formula for $D-$factor , which is actively used by reservoir engineers (see
Section 6.1) We also formulate general algorithm for well pressure
interpretation. It is worth mentioning that some arguments for numerical data
interpretation are routinely used by reservoir engineers without mathematical
justification. We provide an insight behind picture, which makes clear future
engineering application.

Motivation of the research and brief review are presented in the two sections:
3 \- Well representation in numerical flow simulations, 4 \- Review on the
well transmissibility calculations. Section 7 gives the Concluding remarks.
For convenience we listed all notation in the Section 2.

##  2\. Nomenclature

The following notation is used

$\displaystyle B(0,R)=\\{|X|\leq R\\}$ $\displaystyle X=(x,y)$ $\displaystyle
v_{r}\ -\ \text{radial flow (Darcy) velocity}$ $\displaystyle p(x,y)\ -\
\text{pressure function}$ $\displaystyle\delta\ -\ \text{delta function}$
$\displaystyle\Delta\ -\ \text{grid size}$ $\displaystyle k\ -\
\text{permeability of the porous medium (reservoir)}$ $\displaystyle h\ -\
\text{reservoir thickness}$ $\displaystyle\mu\ -\ \text{viscosity of the
fluid}$ $\displaystyle R_{w}\ -\ \text{well radius}$ $\displaystyle\theta\ -\
\text{imaginary well radius }$ $\displaystyle\tilde{p}\ -\ \text{numerical
solution by discreet scheme of interest defined at discreet points }(i,j).$
$\displaystyle f(s)\asymp g(s),\ \text{if}\ \lim_{s\to
0}\frac{f(s)}{g(s)}=constant>0$

##  3\. About well representation in a 3D flow model

The problem of well representation in numerical flow modelling has long
history and peculiar difficulties. Unfortunately, even today it is not
possible to confirm that all the perplexities have been already overcome. Well
modeling is a classical problem of representing objects of different scale
within a single computational module with the requirement to achieve necessary
accuracy of simulations. For example, in petroleum reservoir simulation grid
block sizes in the horizontal plane are typically 2-3-4 orders of magnitude
larger than the well diameter. A typical well radius is of the order of 10 cm,
and common elementary numerical grid block dimensions are $10^{2}$x$10^{2}$x1
m.

Progress in the availability of computing facilities has led to an increase in
the size of the grid and, to a lesser extent, to a decrease in the size of
individual grid cells. Over the years, the typical size of a grid cell has
decreased insignificantly, while the model sizes (number of grid blocks) have
increased by 2-3 orders of magnitude

At the same time, pressure gradients are maximal in the near-well zones,
requiring their most accurate representation in these areas. The farther the
grid cell is from the well, the less reservoir pressure change in it is
associated with a particular change in well operation (sink/source intensity).

Wells are almost never explicitly simulated. This means that pressure in a
grid block penetrated by a well differs significantly from the bottom-hole
pressure in the same well. This pressure difference causes influx/outflux
to/from the grid block. So the problem of well model arises. This model
relates well flow rate to the pressure difference that forces this rate. To
the best of our knowledge, the problem of transition from the well grid block
pressure to the bottom-hole pressure of an oil well was originally
comprehended in the former USSR. Firstly, it was investigated in relation to
2D flow problems solution on electrical integrators with RC (resistance-
capacitance) grids [10]. Later similar studies were performed for 2D flow
problems solution with finite-difference methods [11, 12]. Same approach to
gas wells simulation was considered in [13]. In much more recent foreign
publications, well connection transmissibility factor was introduced into the
formula for well flow rate in a grid block penetrated by the well. The latter
was multiplied by the difference between the grid block pressure and the well
bottom hole pressure. Well connection transmissibility plays an important role
in the flow modeling. It is usually designated as WI - abbreviation for “well
index”. Or alternatively as CF - abbreviation for "connection transmissibility
factor". The WI (CF) $T_{w}$ is defined as the ratio between the well
volumetric flow rate (at reservoir conditions) $q_{w}$ and the difference
between the grid block pressure $p_{0}$ and the well bottom hole pressure
$p_{w}$ for single-phase flow of a unit-viscosity fluid:

(3.1) $T_{w}=\frac{q_{w}\mu}{p_{0}-p_{w}}$

or

(3.2) $q_{w}=\frac{T_{w}}{\mu}(p_{0}-p_{w})$

Formulas (3.1)-(3.2) express the essence of well modeling in numerical flow
simulations. A well is usually not simulated explicitly with approximation of
its trajectory by a grid mesh. Instead, it only acts as a point-wise
source/sink term in flow (mass balance) equations with the flow rate
calculated according to the formula (3.2). The form of equation (3.2) is very
close to the general expression for the well productivity index PI:

(3.3) $q_{w}=\frac{PI}{\mu}(<p>-p_{w})$

where $Q_{w}$ is the total well flow rate, and $<p>$ is the average reservoir
pressure. Relative affinity of the two equations (3.2) and (3.3) explains why
some methods developed for calculation of PI could also be used for
calculation of $T_{w}$. In different studies, different parameters can be
understood by PI, but in any case this parameter relates the inflow at the
well and the reservoir-well pressure drop.

In general, the inflow to the well is approximated in different ways depending
on the calculation grid. In this paper, we consider only the case of single-
phase incompressible flow for a well penetrating a single grid block of a
uniform Cartesian grid with finite-difference approximation. In this case,
$q_{w}$ is defined as the total flow rate through the cylindrical borehole
surface area with radius $r_{w}$. For a more general case, equation (3.2)
becomes

(3.4) $q_{wi}=\sum_{i}{T_{wi}}M_{i}(H_{i}-H_{wi})$

Here sum is over grid blocks $i$ penetrated by the well. $T_{wi}$ is the CF,
or transmissibility, of the grid block $i$ including contributions of the
block and well geometries, and of the grid block flow properties
(permeability). In many cases $T_{wi}$ is constant, however, it may vary over
time in advanced well models. $M_{i}$ is the mobility in the grid block $i$
including contributions of the fluid properties and fluid-rock interaction. In
the incompressible single-phase case:

(3.5) $M_{i}=\frac{1}{\mu},$

where $\mu$ is the fluid viscosity. In the more complex case, $M_{i}$ also
incorporates characteristics of the two-phase or three-phase flow and depends
on the grid block pressure, fluid saturations (volume fractions) and component
concentrations. $H_{i}-H_{wi}$ is the difference of total fluid potentials
(including pressure and gravity) between the grid block $i$ and the well
connection in it. In the general case, $H_{wi}$ is related to the well
bottomhole pressure $p_{w}$ through a hydrostatic correction and friction
losses, and it depends on the trajectory of the well.

##  4\. Brief review of the connection transmissibility computation problem

In the simplest case, for a vertical well in a homogeneous reservoir, an
analytical formula for calculating $T_{w}$ was presented in the D.Peaceman’s
paper [14]. The Peaceman’s approach is described in detail below. It deals
with calculation of so-called equivalent radius $R_{0}$, so that the well-
block pressure $p_{0}$ can be interpreted as the pressure at $R_{0}$ in the
analytical solution for axisymmetric flow around the well. In the pre-Peaceman
period, the approach of paper [15] was widely used, with areally averaged
pressure used instead of pressure at the equivalent radius. Inaccuracy of this
approach was clearly shown in [14]. Later, flow simulation software (flow
simulators) included only Peaceman-type well models.

For the sake of justice, it was noted above that the problem of individual
well consideration in flow simulations was first solved for vertical wells in
the USSR, long before the appearance of similar papers in the United States.
However, only the papers of D.Peaceman [14, 16] are generally referred to in
this connection.

A number of authors have addressed the problem of well CF calculation. Paper
[16] generalizes results of [14] for the case of a single fully penetrating
vertical well located in a Cartesian grid with orthotropic-anisotropic medium.
Modifications taking into account alternative well orientations were presented
in papers [17, 18]. Although those approaches based on spatial variables
scaling are still used in existing commercial simulators, in some cases they
can lead to significant errors [19]. However, classical Peaceman model does
not cope with complex cases, for example, off-centered wells located not in
the center of a grid block [20, 21], or generalized-form grid geometries, for
example, deformed grid blocks or local grid refinement near the well [20, 22].

Approaches combining semi-analytical solutions and finite difference
calculations for regions exceeding a single grid block are much more accurate.
In [23], productivity index of a horizontal well parallel to a coordinate
axis, as previously derived in [24], was used as a control solution. The
authors of [23] combined this solution with analytical expression for pressure
distribution in a cross-cut problem on a Cartesian grid with large aspect
ratio (width to height) of grid cells. This made possible calculation of well
productivity index by cross-sections of the three-dimensional problem.

There is some gap in representation of the solutions obtained within the
classical Peaceman model when further applied to more complex cases.

Based on the comparison of analytical and numerical results, accurate values
of well connection factor in homogeneous reservoirs were determined in [25].
For inhomogeneous case, reference solution was obtained by small-scale
modeling on Voronoi grids. In [27], a single inclined well was considered in a
laterally infinite layered system. Semianalytical solution derived with the
reflection method and analytical solution of the thin-layer theory was
combined to the finite-difference pressure solution, with lateral boundary
conditions for pressure determined from the analytical solution. In [19],
Laplace equation was solved in elliptic coordinates using geometric
transformation taking into account well inclination in an infinite reservoir.
This result was used together with numerical solution to determine appropriate
well connection transmissibility. Single-layer potential was applied in [27]
to calculate stationary pressure distribution in the vicinity of a well.
Additionally, grid block transmissibilities were also matched to take into
account specifics of radial flow. Later in paper [28] it was proposed to
subdivide well inflow problem into two auxiliary subproblems, so that
superposition of both subsolutions gives a solution of the original problem.
The first subproblem represents a singular flow caused by well operation. And
the second one considers regular flow not depending on the well. More
precisely, the first problem consists in pressure determination in the case of
given inflow rate in an infinite reservoir. The second problem corresponds to
well operation with zero flow rate in a limited reservoir, and boundary
conditions at the outer boundary are adjusted to the first problem solution at
this boundary. The main distinctive feature of this method lies in
transformation of Cartesian coordinates into logarithmic-polar (radial)
coordinates. Then it becomes possible to avoid singularity of the solution in
the vicinity of a well. However, the second auxiliary problem becomes somewhat
worse than on the original grid, since the approximation error increases.
Approximation in the new coordinate system is carried out with standard
methods of two-point or multi-point approximation, adjusted for logarithmic
polarity of the coordinate system. To model a well, a Peaceman-type model was
employed in [28]. The value of the equivalent radius for the well block could
be arbitrary. The only requirement was that the corresponding circle lied
inside the grid block and was larger than the circle of the wellbore.

Exact solutions with grids taking into account well trajectories and near-well
heterogeneity have also been discussed in the literature (confer finite
element-based grids [29] and locally elliptical hybrid grids [30]). Regarding
wells of non-conventional types, i.e. with arbitrary trajectory or
multilateral wells, it should be noted that modern technologies are successful
in their construction. Therefore, they become more and more common. This type
of wells can penetrate grid blocks in any direction. In addition, this type of
well models often use irregular grids to numerically represent various
geological features. These grids cover the range from structured curved
meshes, multi level grids to completely unstructured grids. Combined
application of wells with complex trajectory and progressive grids leads to
significant complications in achieving acceptable accuracy for predicting
future well behavior.

Effective approach to modeling productivity of unconventional well types
during reservoir depletion consists in the use of semi-analytical approaches.
Early studies considered single horizontal wells (of infinite conductivity),
parallel to one of the sides of a rectangular reservoir. Solution methods
employed sequential integral transformation [31, 32] and point Green functions
[33, 34, 35], leading to solutions in the form of series. More complex well
geometries were considered later in [36, 37, 38] using numerical integration
of differential equations. Well hydraulics (i.e., limited well conductivity)
was combined in [39] with flow in the reservoir. All the mentioned semi-
analytical methods have advantages of limited input datasets and high
computational efficiency compared to finite difference modeling. This makes
these methods well suited for preliminary assessment of primary production
during reservoir depletion.

Early semi-analytical methods, however, were limited to homogeneous systems or
almost strictly layered systems [26, 39]. This was a significant limitation,
since productivity of unconventional well types can be significantly affected
by small-scale heterogeneity in the near-well area. Small-scale heterogeneity
could be included in a detailed computational model. But the resulting model
could become very difficult to be built and require significant computational
time to run.

In this paper we recall the original Peaceman problem to provide its accurate
mathematical interpretation and treatment. The corresponding "sewing
machinery" is general and can be extended to more complex cases. As an
example, we consider the practically-important case of nonlinear Forchheimer
flow around the well. Accurate results are obtained for the equivalent radius
and CF (or inflow formula) and compared to the Dake’s formula generally used
in reservoir simulators.

To some extent, this study is based on similar principles with the
mathematical treatment of PI calculations presented in [6] for linear Darcy
flow case in the detail for all three regimes: pseudo-steady, boundary
dominated, and steady state. In case of non-linear Forchheimer flow, the
corresponding framework makes sense only for two regimes: pseudo-steady state
and steady state, and is presented in [8]. For gas, corresponding steady was
presented in [9]. Main issues here were that in spite of essentially time-
dependent nature of the solution, the productivity index, as it was rigorously
proven in cited paper, is time-independent regardless of initial data, or
tends to a steady state value exponentially. Therefore, PI can very naturally
be used in numerical simulations of the flow to tune a numerical model to real
life (observed) data of the reservoir.

##  5\. Peaceman Well Block Radius and Fundamentals in Finite Difference
Solution

###  5.1. Origin of the Peaceman Formula Revisited

Let’s in our own words reformulate main idea of the Peaceman paper [14]. The
aim is to prepare ground for proposed mathematical framework. All other
reports from the Peaceman Well Block radius series of articles [14], [16], and
articles which follow, although having high importance, are actually
extensions and developments of the Peaceman paradigm which we will present in
this section. Peaceman considered the problem of numerical simulation of flow
towards a well, which is subject to Darcy equation. Total rate of production
is fixed. In [14] flow is assumed symmetrical w.r.t. the well located in the
center block $0$ of a 5-point stencil of a numerical finite difference grid as
shown in fig.2.

The finite difference approximation is as in (1). Due to symmetry it is
reduced to equation (1.2), or for arbitrary reservoir thickness

(5.1) $p_{1}-p_{0}=\frac{\mu}{4kh}q,\ \text{ here}\ \alpha=\frac{\mu}{kh},$

Consider the vertical well of radius $R_{w}$ situated in the center of the
block $0$ and causing the flow in radial direction only. Assume that $q$ is
the fixed total production rate. Then flow in the domain
$U(0,R_{1},R_{w})=B(0,R_{1})\setminus B(0,R_{w})$ is uniquely characterised by
the pressure distribution

(5.2) $p(r)=\frac{\mu q}{2kh\pi}\ln\frac{r}{R_{1}}+p_{1},\ \text{ here }p_{1}\
\text{is the value of pressure at}\ r=R_{1}.$

Now let $R_{1}=\Delta$, so that $p_{1}$ is the same as in (5.1), and then

(5.3) $p(R_{w})=\frac{\mu q}{2kh\pi}\ln\frac{R_{w}}{\Delta}+p_{1}.$

Peaceman paradigm in short can be stated as

######  Main Problem 1.

Find $R_{0}$ s.t. the value of $p_{0}$ in (5.1) provides

(5.4) $p(R_{w})=\frac{\alpha q}{2\pi}\ln\frac{R_{w}}{R_{0}}+p_{0}.$

######  Theorem 1.

Assume that total rate of the production $q$ and step size of the grid
$\Delta$ are given. Assume that the single fully penetrating well is located
at a point chosen as the origin and so being the center of the numerical block
$[-\frac{\Delta}{2},\frac{\Delta}{2}]^{2}$. Let $R_{w}$ is such that
$\ln\frac{\Delta}{R_{w}}>\frac{\pi}{2}$. Then necessary and sufficient
condition that guarantees existence of the solution of the system (5.1)-(5.4)
for $R_{0}$ is $\ln\frac{\Delta}{R_{0}}=\frac{\pi}{2}.$

###### Proof.

Proof follows from direct equivalent algebraic operations on the equations
(5.1),(5.4), and (5.3). ∎

######  Remark 1.

In this section interpretation of the Peaceman paper [14] is made directly
without significant modification. But it is already clear that the main idea
is to sew numerical and analytical solutions formulated on different scales to
provide required information for tuning procedure between the numerical
solution and observed well pressure data. In what will follow in this paper,
we will bring more mathematical insight in the Peaceman sewing machinery, and
provide a general framework suitable for non-Darcy flows.

###  5.2. Green function and related interpretation of the well-block radius

Mathematically speaking, Peaceman numerical scheme as presented in [14] is a
finite difference representation for the Green function approximated on the
square block with the source at center. In this subsection we state the
corresponding problem in PDE format and provide its relation to the numerical
solution. Then we will state a problem, provided that $G(x,0)$ is a Green
function with source at the origin, of finding the value of $R_{0}$ as a
function of $\Delta$ such that oscillation of the $G(x,0)$ in the annular
domain $U(0,R_{0},\Delta)=U\cap B(0,R_{0},\Delta)$ almost constant for
$\Delta\asymp 0$ and $\Delta\text{ such that}\ B(0,\Delta)\in U.$

Consider BVP in the bounded domain $U$ containing source point $0$ (see fig.
2) :

(5.5) $\displaystyle\frac{kh}{\mu}\cdot\Delta p=q^{0}\cdot\delta(x)\ \text{in
the domain \ $U$}$ (5.6) $\displaystyle p(x)=0\ \text{on the boundary}\
\partial U$

Figure 2. Five-point stencil for the finite difference scheme, and domain of
the flow with source at $0$

We divide the whole area of flow into $M\times M$ blocks. For all blocks
$0\leq i\leq M,$ $0\leq j\leq M$ (see fig. 2), the finite difference equation
for (5.5) will take form:

$\displaystyle\frac{kh}{\mu}\cdot\frac{(p_{i+1,j}-2p_{i,j}+p_{i-1,j})}{(\Delta
x)^{2}})+$ (5.7)
$\displaystyle+\frac{kh}{\mu}\cdot\frac{(p_{i,j+1}-2p_{i,j}+p_{i,j-1})}{(\Delta
y)^{2}}=\frac{q_{i,j}}{\Delta x\cdot\Delta y}$

or in the form of Peaceman:

$\displaystyle\frac{kh\Delta y}{\mu\Delta
x}\cdot(p_{i+1,j}-2p_{i,j}+p_{i-1,j})+$ (5.8) $\displaystyle+\frac{kh\Delta
x}{\mu\Delta y}\cdot(p_{i,j+1}-2p_{i,j}+p_{i,j-1})=q_{i,j}$

In the above equation, $q_{i,j}=0$ if ${i\neq 0}$ or $j\neq 0$
($q_{i,j}=q\cdot\delta_{i,j}$ \- Kronecker symbol). Evidently, sizes of the
block in $x$ and $y$ direction are correspondingly $\Delta x$ and $\Delta y$
and converge to $0$ as $M\to\infty$. Let us denote $2M\times 2M$ matrix
$P_{M}$

(5.9) $P_{M}=\Big{(}p_{i,j}\Big{)}_{\left((-M\leq i\leq M);(-M\leq j\leq
M)\right)}$

Elements of the matrix represent values of the solution of the discreet
Poisson equation (5.5) with sink/source localized at center $(0,0)$.

Johansen Taousch from SMU proposed that using classical machinery(see [3]) of
splitting Green function $G(x,y)$ into the sum of fundamental solution and
"corrector" one can prove the following.

######  Theorem 2.

Let $(x,y)\neq(0,0)$ is a fixed point of the domain $U$. This point belongs to
one element (block) of the grid $U_{M}$ which approximates the domain $U.$ Let
$p_{M}(x,y)$ is the solution of the system $2M\times 2M,$ extended to
$C^{2}(U)\cap C^{0}(\bar{U})$. Then as $M\to\infty$, function $p_{M}(x,y)\to
G(x),$ where

(5.10) $G(x)=\frac{1}{2\pi}\ln\frac{1}{r}+g(x,y),\ r=\sqrt{x^{2}+y^{2}}.$

is the Green function. Function $G$ has the property

(5.11) $\int_{\gamma}q^{0}\frac{\partial G}{\partial\nu}ds=q^{0},\ \text{where
}\gamma\text{ is a closed contour around the source at}\ 0.$

We do not prove this Theorem in the article but will address this issue in
upcoming publications.

Green function has an interesting property, which directly relate to Peacamn
well-block radius.

######  Theorem 3.

In the axisymmetric case, Green function has the property: There exists
$R_{0}$ s.t. if

(5.12) $4\cdot\left(G(x)|_{|x|=\Delta}-G(y)|_{|y|=R_{0}}\right)=1+O(\Delta)$

then

(5.13) $\ln\frac{\Delta}{R_{0}}=\frac{\pi}{2}+O(\Delta)$

###### Proof.

First observe that

(5.14)
$\left(G(x)|_{|x|=\Delta}-G(y)|_{|y|=R_{0}}\right)=\frac{1}{2\pi}\ln\frac{1}{\Delta}-\frac{1}{2\pi}\ln\frac{1}{R_{0}}+\left(g(x)\big{|}_{|x|=\Delta}-g(x)\big{|}_{|x|=R_{0}}\right).$

After regrouping from above follows

(5.15)
$\left(G(x)|_{|x|=\Delta}-G(y)|_{|y|=R_{0}}\right)=\frac{1}{2\pi}\ln\frac{R_{0}}{\Delta}+\left(g(x)\big{|}_{|x|=\Delta}-g(x)\big{|}_{|x|=R_{0}}\right).$

Using standard appriory estimate for corrector $g(x)$ as a solution of the
Dirichlet problem one has

(5.16)
$\Big{|}g(x)\big{|}_{|x|=\Delta}-g(x)\big{|}_{|x|=R_{0}}\Big{|}\leq\max_{x,y\in{U(0,\Delta)}}\big{|}g(x)-g(y)\big{|}\leq
C\big{|}x-y\big{|}\leq 2C\cdot\Delta.$

Constants $C$ is constant depending on the distance between $0$ and $\partial
U$. ∎

###  5.3. Radial Darcy flow in two embedded annulus zones sewed by linear
balance equation on the finite difference coarse grid, and its link to
capacity

In this section we will highlight sewing features behind the well block
radius, and indicate its link to capacity from the analytical setup which
comes out of the Peaceman method.

Consider auxiliary annular domain $U(0,\theta,\Delta)$, with $\theta<\Delta,$
to be selected by technical arguments. To use Peaceman well block radius as a
sewing machinery, assume that flow is radial around the virtual well-disc of
radius $\theta$ and is generated by given total flux $q$ on
$\Gamma_{\theta}=\\{x:|x|=\theta\\}$. Consequently between $\Gamma_{\theta}$
and exterior boundary $\Gamma_{e}$ pressure drop will occur specified by given
$p_{w}\text{ and}\ p_{e},$ or one of those pressures and the total rate $q$.
In the case of fully radial flow also consider the domain $U(0,R_{w},R_{e})$:
$\Gamma_{e}=\\{r=R_{e}\>>\Delta\\}$ and $\Gamma_{w}=\\{r=R_{w}<\Delta\\}$,
with actual well radius $R_{w}$.

Direct integration of the Darcy equation for radial flow (5.17)

(5.17) $-\frac{\partial p}{\partial r}=\alpha v_{r}=\alpha\frac{q}{2\pi r}$

over the domain $U(0,r,R_{e})$ provides

(5.18) $p(r)-p(R_{e})=\frac{\alpha q}{2\pi}\ln\frac{r}{R_{e}}.$

for any $r>\theta,\ r=\sqrt{x^{2}+y^{2}}$ .

Now consider finite difference scheme with axial symmetry. Let

(5.19) $\Omega_{0}=\\{(x,y):-L\leq x\leq L\ ;-L\leq y\leq
L\\},\Delta_{x}=\Delta_{y}=\frac{L}{M}=\Delta\ ,\ M\in\mathcal{N}.$

We will impose constraint that for given radius $R_{w}$ of the "actual" well
$\Delta<\frac{L-R_{w}}{L}.$ Reasoning behind this constraint will become clear
later in this section.

In the rectangular domain $\Omega_{0}$ consider finite difference BVP

(5.20)
$\displaystyle\frac{kh}{\mu}\cdot(p_{i+1,j}-2p_{i,j}+p_{i-1,j})+\frac{kh}{\mu}\cdot(p_{i,j+1}-2p_{i,j}+p_{i,j-1})=q_{i,j}=q\delta_{i,j}$
(5.21) $\displaystyle p_{-M,j}=p_{M,j}=p_{i,-M}=p_{i,M}=0\text{ - boundary
conditions.}$ In above: (5.22) $\displaystyle\\{-M\leq i\leq M\ ,\ -M\leq
j\leq M\\}=I\times J\,$ (5.23) $\displaystyle\delta_{i,j}-\text{index function
of the block $B(0,0)$}:\delta_{0,0}=1,\delta_{i,j}\text{ equals}\ 0\
\text{otherwise}.$

Denote $2M\times 2M$ matrix

(5.24) $P_{M}=\Big{(}p_{i,j}\Big{)}_{\left((-M\leq i\leq M);(-M\leq j\leq
M)\right)}$

Evidently matrix $P_{M}$ exists and is unique. Assume $\theta$ be fixed
parameter, and

(5.25) $\Delta=\frac{L}{M}>\theta.$

Let $B(0,\theta)$ be the disc model of the well of radius $\theta$. Following
tradition in engineering, we will call the five-point difference equation
(5.20) for $i=0,j=0$ \- the material balance equation.

Matrix $P_{M}$ is $2D$ approximation of the analytical problem

$\displaystyle\Delta p=0\ \text{in}\ \Omega_{0}\setminus B(0,\theta);$ (5.26)
$\displaystyle h\frac{k}{\mu}\int_{\Gamma_{\theta}}\frac{\partial
p}{\partial\nu}ds=q;$ $\displaystyle p=0\ \text{on}\ \partial{\Omega_{0}}.$

Let $\theta<<L,$ then it is reasonable to assume that
$p(x,y)\big{|}_{|x|=\theta}\asymp constant.$

######  Assumption 1.

Assume symmetry condition in the four blocks surrounding the center block
$U(0)=[-\frac{\Delta}{2};\frac{\Delta}{2}]\times[-\frac{\Delta}{2};\frac{\Delta}{2}]$(see
Fig. 2 ):

(5.27) $p_{-1,1}=p_{1,-1}=p_{-1,-1}=p_{1,1}=p_{1}.$

######  Definition 1.

Denote $p_{0,0}=p_{0}$, then finite difference scheme in the five spot grid
containing disc(well) $B(0,\theta)$ in the central block $B_{0,0}$ will take a
form of the symmetric material balance in the form :

(5.28) $p_{1}-p_{0}=\frac{\alpha}{4}q.$

As it was stated, $p_{i,j}$ is uniquely defined by input parameters
$M,q,L,k,\mu,h.$ Note that $p_{1}$ as well as $p_{0}$ depend on the
characteristics of the grid size - $\Delta$, but due to material balance
equation their difference is $\Delta$-independent.

Trace of the solution of the analytical problem (5.26) on the $|x|=\Delta$
depends on $\Delta$. We need to find $R_{0}$ as a function of $\Delta$ such
that the difference between traces of the analytical solution on $|x|=\Delta$
and $|x|=R_{0}$ is $\Delta$-independent as well. To do so let us introduce two
auxiliary BVPs:

$\displaystyle\Delta u_{1}=0\ \text{in}\ U(\Delta,\theta)=B(0,\Delta)\setminus
B(0,\theta);$ (5.29) $\displaystyle
h\frac{k}{\mu}\int_{\Gamma_{\theta}}\frac{\partial u}{\partial\nu}ds=q;$
$\displaystyle u_{1}=p_{1}\ \text{on}\ \partial{B(0,\Delta)}.$

$\displaystyle\Delta u_{0}=0\ \text{in}\ U(R_{0},\theta)=B(0,R_{0})\setminus
B(0,\theta);$ (5.30) $\displaystyle
h\frac{k}{\mu}\int_{\Gamma_{\theta}}\frac{\partial u_{0}}{\partial\nu}ds=q;$
$\displaystyle u_{0}=p_{0}\ \text{on}\ \partial{B(0,R_{0})}.$

######  Definition 2.

Let imaginary well of the radius $\theta$ ($\theta<R_{0}<\Delta$) and total
rate of the production $q$ are fixed. Denote boundary of the well
$\Gamma_{\theta}=\\{|x|=\theta\\}.$ We will say that Peaceman problem in the
annular domains $U(R_{0},\theta)$ and $U(\Delta,\theta)$ is strictly well
posed if $R_{0}$ exists s.t. simultaneously two constrains hold:

  1. (i)

solutions of the problems (5.29) and (5.30) on $\Gamma_{\theta}$ are constant
and

$u_{1}\big{|}_{r=\theta}=u_{0}\big{|}_{r=\theta}=p_{\theta},$

  2. (ii)

$p_{1}$ and $p_{0}$ are solutions of the balance equation (5.28).

Another interpretation of the Definition of well-posedness contains sewing
property of the material balance equation, which allows to sew the analytical
solutions of the problems 5.29 and 5.30 via values on the well boundary.

######  Theorem 4.

Assume balance equation is as given by (5.28). Let imaginary well of the
radius $\theta$ to be at the center of a five spot grid. Let number of the
blocks $M$, size of the square domain $L$, and imaginary well of radius
$\theta$ be such that

(5.31) $\Delta=\frac{L}{M}>e^{\frac{\pi}{2}}\cdot\theta.$

Then for given rate $q$, if

(5.32) $R_{0}=\Delta\cdot e^{-\frac{\pi}{2}},$

then Peaceman problem is strictly well posed.

###### Proof.

Due to radial flow one has (5.17), and consequently (5.18), from which we get
the system

(5.33) $u_{1}(\Delta)-u_{1}(\theta)=p_{1}-u_{1}(\theta)=\frac{\alpha
q}{2\pi}\ln\frac{\Delta}{\theta}.$ (5.34)
$u_{0}(R_{0})-u_{0}(\theta)=p_{0}-u_{0}(\theta)=\frac{\alpha
q}{2\pi}\ln\frac{R_{0}}{\theta}$

Due to balance equation, in order $R_{0}$ to exist it is sufficient that

(5.35) $\frac{\alpha}{4}q-\frac{\alpha
q}{2\pi}\left(\ln\frac{\Delta}{\theta}-\ln\frac{R_{0}}{\theta}\right)=u_{1}(\theta)-u_{0}(\theta).$

or

(5.36)
$\frac{\alpha}{2}q\left(\frac{1}{2}-\frac{1}{\pi}\ln\frac{\Delta}{R_{0}}\right)=u_{1}(\theta)-u_{0}(\theta).$

In order the Peaceman problem to be strictly well-posed, it suffices
expression in the parentheses to be equal zero:

(5.37) $\left(\frac{1}{2}-\frac{1}{\pi}\ln\frac{\Delta}{R_{0}}\right)=0\
\Longleftrightarrow\frac{\pi}{2}=\ln\frac{\Delta}{R_{0}}.$

QED ∎

######  Remark 2.

Due to material balance equation there is no difference between $p_{1}$ and
$p_{0}$ for $\Delta>e^{\frac{\pi}{2}}\cdot r_{w}$. But from point of view of
the engineering it is important to obtain value of the simulated well pressure
from pressure values of the simulated numerical solution, more precisely, from
value of $\bar{p}_{0}$ at nearest block containing well as a function of given
$\Delta$. Algorithm will then consist of two possible scenarios:

  1. (i)

$\Delta>e^{\frac{\pi}{2}}r_{w}$. We will find value of the pressure on the
well by

(5.38) $p_{w}=p_{0}-\frac{\alpha
q}{2\pi}\ln\frac{e^{-\frac{\pi}{2}}\cdot\Delta}{r_{w}}$

  2. (ii)

$r_{w}\leq\Delta\leq e^{\frac{\pi}{2}}r_{w}$. Then we will let
$p_{w}=p_{1}-\frac{\alpha q}{2\pi}\ln\frac{\Delta}{r_{w}}$, or due to balance
equation (5.28) :

(5.39) $p_{w}=p_{0}+q\frac{\alpha}{4}-\frac{\alpha
q}{2\pi}\ln\frac{\Delta}{r_{w}}.$

######  Remark 3.

It is also worth to mention that if $\Delta\to 0$ then $\lim_{\Delta\to
0}p_{0}(\Delta)$ exist and uniquely determined by $q$. It follows from the
following arguments for simple annular domain. Namely let the radius of
imaginary well $\theta=r_{w}$ .Then from Theorem 4 it follows that

(5.40) $p_{w}(\Delta)=p_{0}(\Delta)-\frac{\alpha
q}{2\pi}\ln\frac{R_{0}(\Delta)}{r_{w}}.$

###  5.4. Two Terms Radial Case

Now consider Peaceman problem for nonlinear flow by Forchheimer two-terms law:

(5.41) $-\frac{\partial p}{\partial
r}=\alpha_{1}v_{r}+\beta_{1}v_{r}|v_{r}|,\text{
}\alpha_{1}=\frac{\mu}{k}=\alpha h$

In (5.41) if $\beta_{1}=0$ one can get classical Darcy equation. From (5.41)

(5.42) $\frac{\partial p}{\partial r}=\alpha\frac{q}{2\pi
r}+\beta\frac{q}{4\pi^{2}r^{2}},\text{ }\beta=\frac{\beta_{1}}{h^{2}}$

Consider flow from $\partial B(0,R_{2})$ to $\partial B(0,R_{1})$ in the
annular domain $U$(see Fig. 3):

(5.43) $\displaystyle\begin{cases}&U=B(0,R_{2})\setminus B(0,R_{1}),\
\text{with fixed total rate }\ q=\int_{S}v(r)ds,,\\\ &\text{Given pressure on
one of the boundaries}\ \partial B(0,R_{i}):\\\ &p(r)\big{|}_{r=R_{i}}=f_{i}\
\text{for}\ i=1\ \text{or}\ 2.\end{cases}$

Explicit formula for generic solution of (5.43) for two terms Forchheimer law
follows from basic integration:

(5.44) $p(r)=\frac{\alpha q}{2\pi}\ln
r-\beta\frac{q}{4\pi^{2}}\frac{1}{r}+constant$

Then using boundary conditions in (5.43) one can get a generic formula for
pressure drop between two contours ($\partial B(0,R_{i})$) of the boundary of
annular domain $U(0,R_{1},R_{2})=B(0,R_{2})\setminus B(0,R_{1})$.

(5.45) $p\big{|}_{r=R_{2}}-p\big{|}_{r=R_{1}}=f_{2}-f_{1}=\frac{\alpha
q}{2\pi}\ln\frac{R_{2}}{R_{1}}+\beta\frac{q}{4\pi^{2}}\left(\frac{1}{R_{1}}-\frac{1}{R_{2}}\right)$

Figure 3. General Annual Domain $U$

###  5.5. Criteria for radial non-linear Forchheimer flow in two annular zones
with linear balance equation on the finite difference coarse grid

In this section we will derive Peaceman-type well pressure formula taking into
account non-linearity of the flow near well using Forchheimer correction of
Darcy law by two term equation (5.41). Although we could use same approach as
for the Darcy case, our treatment for non-linear Forchheimer will be
different. Namely we will not introduce imaginary radius $\theta$ and
corresponding Peaceman well-posedness definition. We do it purposely to
highlight main assumption and keep $\theta=R_{w}$.

Since two-term non-linearity of flow is essential only in near-well high
velocity zones, we will adopt conventional engineering hypothesis that on
numerical grid material balance can be still considered linear, while near-
well correction is due to Forchheimer type of non-linearity. From (5.1) and
(5.45) with $\beta\neq 0$ the system of 3 equations follows:

(5.46) $p_{1}-p_{0}=\frac{\alpha}{4}q$ (5.47)
$p_{1}-p_{w}=\frac{\alpha}{2\pi}q\ln\frac{\Delta}{R_{w}}+\beta\frac{q^{2}}{4\pi^{2}}\left(\frac{1}{R_{w}}-\frac{1}{\Delta}\right)$
(5.48)
$p_{0}-p_{w}=\frac{\alpha}{2\pi}q\ln\frac{R_{0}}{R_{w}}+\beta\frac{q^{2}}{4\pi^{2}}\left(\frac{1}{R_{w}}-\frac{1}{R_{0}}\right)$

######  Theorem 5.

Assume that $q$ and $p_{w}$ solve quadratic equation

(5.49)
$p_{1}-p_{w}=\frac{\alpha}{2\pi}q\ln\frac{\Delta}{R_{w}}+\beta\frac{q^{2}}{4\pi^{2}}\left(\frac{1}{R_{w}}-\frac{1}{\Delta}\right)$

Then if $R_{0}$ satisfies equation

(5.50) $R_{0}=\Delta\cdot e^{-\delta\frac{\pi}{2}}$

system (5.46)-(5.48) has a solution for any mutually related $p_{1},\
p_{0},\text{and}\ p_{w}$, if $\delta$ satisfies equation

(5.51)
$\delta+\beta\frac{q}{\alpha\pi^{2}}\left(\frac{e^{\delta\frac{\pi}{2}}}{\Delta}-\frac{1}{\Delta}\right)=1$

###### Proof.

Indeed, let first substitute (5.50) into RHS of the (5.48) to get

(5.52) $\displaystyle
p_{0}-p_{w}=\left(\frac{\alpha}{2\pi}q\ln\frac{\Delta}{R_{w}}+\beta\frac{q^{2}}{4\pi^{2}}\left(\frac{1}{R_{w}}-\frac{1}{\Delta}\right)\right)$
(5.53)
$\displaystyle+\left(\frac{\alpha}{2\pi}q(-\delta\frac{\pi}{2})+\beta\frac{q^{2}}{4\pi^{2}}\left(-\frac{1}{\Delta\cdot
e^{-\delta\frac{\pi}{2}}}+\frac{1}{\Delta}\right)\right)$

First term in parentheses in RHS of (5.52) is equal to $p_{1}-p_{w}$, and the
second one we will rewrite as

(5.54)
$p_{0}-p_{w}=p_{1}-p_{w}+q\left(\frac{\alpha}{4}(-\delta)+\beta\frac{q}{4\pi^{2}}\left(-\frac{1}{\Delta\cdot
e^{-\delta\frac{\pi}{2}}}+\frac{1}{\Delta}\right)\right)$

or

(5.55)
$p_{0}-p_{1}=-q\frac{\alpha}{4}\left(\delta+\beta\frac{q}{\alpha\pi^{2}}\left(\frac{1}{\Delta\cdot
e^{-\delta\frac{\pi}{2}}}-\frac{1}{\Delta}\right)\right).$

Finally we have

(5.56)
$p_{1}-p_{0}=q\frac{\alpha}{4}\left(\delta+\beta\frac{q}{\alpha\pi^{2}}\left(\frac{1}{\Delta\cdot
e^{-\delta\frac{\pi}{2}}}-\frac{1}{\Delta}\right)\right).$

To get (5.46) and consequently to finish the proof it will be sufficient to
assume that $\delta$ satisfies equation(5.51). ∎

######  Remark 4.

  1. (i)

As $\beta\to 0$ factor $\delta\to 1$, and we consequently get Peaceman well-
block radius as in the Theorem 4.

  2. (ii)

Evidently non-linear functional equation (5.51) has a solution for any $q>0.$

##  6\. Impact of the Correction in the Dake Formula

###  6.1. Dake Formula and D-factor

To compute fluid inflow to well in the case of Forchheimer non-linear flow,
reservoir simulators use the so-called D-factor, or dynamic skin-factor. It is
assumed that Peaceman radius is still applicable.

However, a nonlinear (rate-dependent) correction is introduced in the
denominator of the inflow formula. The derivation can be found in [4] and is
briefly given below. As before, we neglect compressibility of the fluid for
simplicity. However, the derivation can be directly extended to compressible
case by introduction of pseudo-pressure [4].

Assume flow is radial and $R_{e}$ is radius of a contour with specified
pressure $p_{e}$. Equation (5.41) is integrated from $R_{w}$ to $R_{e}$ to
relate pressure drop $p_{e}-p_{w}$ to flow rate $q=2\pi rhv_{r}$

(6.1)
$p_{e}-p_{w}=\frac{\alpha}{2\pi}q\ln\frac{R_{e}}{R_{w}}+\beta\frac{q^{2}}{4\pi^{2}}\left(\frac{1}{R_{w}}-\frac{1}{R_{e}}\right)$

Now if $R_{e}$ is large enough to neglect $\frac{1}{R_{e}}$, the dimensionless
pressure drop becomes

(6.2) $\frac{{2\pi}(p_{e}-p_{w})}{\alpha
q}=\ln\frac{R_{e}}{R_{w}}+\frac{\beta}{\alpha}\frac{q}{2\pi}\frac{1}{R_{w}}$

Here the first term in RHS is equivalent to the Darcy case, and the second
term is a function (linear) of $q$, with

(6.3) $D=\frac{\beta}{\alpha}\frac{1}{2{\pi}R_{w}}$

being the D-factor.

In reservoir simulators this result is used for Forchheimer flow as follows.
As in the linear case, inflow to the well is assumed radial within the grid
block containing the well. Formula (6.2) with $R_{e}=R_{0}$ and $p_{e}=p_{0}$
is used to relate well flow rate $q$ to the difference of grid block pressure
$p_{0}$ and well pressure $p_{w}$. The key problem is that $R_{0}$ is computed
by (5.32), just as in the linear case:

(6.4)
$p_{0}-p_{w}=\frac{\alpha}{2\pi}q\left(\ln\frac{R_{0_{Darcy}}}{R_{w}}+Dq\right)$

In fact, the use of $R_{0_{Darcy}}$ in (6.4) is incorrect, since the Peaceman
formula was derived for Darcy flow only. The correct formula for Forchheimer
flow can be obtained by combining (5.48), (5.50) and (5.51).
Recall (5.52)-(5.53) and substitute (5.51) multiplied by $-\frac{\pi}{2}$:

(6.5)
$p_{0}-p_{w}=\frac{\alpha}{2\pi}q\left(\ln\frac{\Delta}{R_{w}}-\frac{\pi}{2}-\beta\frac{q}{2\alpha\pi}\left(\frac{1}{\Delta\cdot
e^{-\delta\frac{\pi}{2}}}-\frac{1}{\Delta}\right)\right)+\beta\frac{q^{2}}{4\pi^{2}}\left(\frac{1}{R_{w}}-\frac{1}{\Delta\cdot
e^{-\delta\frac{\pi}{2}}}\right),$

or rearranging:

(6.6) $p_{0}-p_{w}=\frac{\alpha}{2\pi}q\ln\frac{\Delta\cdot
e^{-\frac{\pi}{2}}}{R_{w}}+\beta\frac{q^{2}}{4\pi^{2}}\left(\frac{1}{R_{w}}-\frac{1}{\Delta\cdot
e^{-\delta\frac{\pi}{2}}}-\frac{1}{\Delta}+\frac{1}{\Delta\cdot
e^{-\delta\frac{\pi}{2}}}\right),$

which finally gives:

(6.7)
$p_{0}-p_{w}=\frac{\alpha}{2\pi}q\ln\frac{R_{0_{Darcy}}}{R_{w}}+\beta\frac{q^{2}}{4\pi^{2}}\left(\frac{1}{R_{w}}-\frac{1}{\Delta}\right).$

This result is very similar to (6.4) except for the last term with
$\frac{1}{\Delta}$. Thus it turns out that for coarse grids with large enough
$\Delta$ the correct formula (6.7) is almost equivalent to the one used by
reservoir simulators. But as $\Delta$ goes smaller, the correction gets
significant, and the D-factor becomes grid-dependent (dependent on $\Delta$).
And note that even for quite large $\Delta$ the correction can become
significant with increasing $q$.

##  7\. Conclusive Remarks and Discussion

We revisited Peaceman well-block radius formula and proved its analog for
Green function in the annular domain. Main algorithm is based on ’sewing’ of
material balance equation (1.2) with the analytical equation of flow. For the
linear case, we introduced the concept of Peaceman well-posedness and
rigorously obtained Peaceman classical well-block radius for Darcy flow near
well and linear material balance equation on coarse numerical grid. We also
considered the following locally nonlinear problem: we kept material balance
equation to be linear but allowed non-linearity in the flow equation (1.3). In
this case, for linear balance equation on the coarse grid and Forchheimer two
terms law near the well, we proved that well-block radius $R_{0}$ should
depend on $q$. We implemented this algorithm in Section 5.5 and Section 6.1,
and obtained formulae for $R_{0}$ which is very close to traditional one used
in flow simulators for big $\Delta$, but quite different for small $\Delta$ or
large $q$. This approach is planned to be extended for different classes of
non-Darcy flows. Wide class of Non-Darcy flows was presented in the articles
[9, 8].

## References

  * [1] Landis, E.M.  _Second Order Equations of Elliptic and Parabolic Type_ , Moscow, Nauka, 1971 (Russian); English transl.: Transl. of Mathematical Monograph, 171, AMS, Providence, RI, 1998.
  * [2] Landkof, N.S.  _Foundations of Modern Potential Theory_ , Moscow, FML, 1966 (Russian); English transl.: Grundlehren der math. Wissenschaften, 180, Springer-Verlag, New York-Heidelberg, 1972.
  * [3] , Evans, L.S.  _Partial Differential Equations_ , American Mathematical Society, 2010.
  * [4] Dake L.P.  _Fundamentals of Reservoir Engineering_ , Volume 8, 1st Edition,1985 ISBN: 9780444418302, eBook ISBN: 9780080568980, Imprint: Elsevier Science.
  * [5] Aulisa, E, Bloshanskaya,L, Hoang, L,Ibragimov, A,  _Analysis of generalized Forchheimer flows of compressible fluids in porous media_ , J. Math. Phys. 50, 103102 (2009).
  * [6] Ibragimov, A., Khalmanova, D., Valkó, P. P., and Walton, J. R.,  _On a mathematical model of the productivity index of a well from reservoir engineering_ ,SIAM J. Appl. Math. 65, 1952 (2005).
  * [7] Ibragimov, A., Soliman, M., Bloshanskaya, L., Siddiqui, R.  _Productivity Index for Darcy and pre-/post-Darcy Flow (Analytical Approach)_ , Transport in Porous Media, (2017).
  * [8] E Aulisa, A Ibragimov, P Valko, J Walton  _Mathematical framework of the well productivity index for fast Forchheimer (non-Darcy) flows in porous media_ Mathematical Models and Methods in Applied Sciences 19 (08), 1241-1275.(2009).
  * [9] Aulisa, E., Bloshanskaya,L., Ibragimov, A.  _Well Productivity Index for Compressible Fluids and Gases_ , Evolution Equations and Control Theory 5(1), (2015).
  * [10] Tolstov, Y.G.  _Application of the method of electrical modeling of physical phenomena to solving some problems of underground hydraulics_ , Tech Magazine. Physics (Russian), volume XII, issue 10, 1942.
  * [11] Vakhitov G.G.  _Solving problems of underground hydrodynamics by the finite difference method_ , Moscow, Proceedings of VNIIneft, issue 10, Gostoptehizda (Russian) pp. 53-88., 1957.
  * [12] Vakhitov G.G.  _Effective ways to solve the problems of developing heterogeneous oil-bearing formations_ , Moscow, Gostoptehizdat (Russian), 1963.
  * [13] Zakirov, S.N.  _On modeling gas wells in violation of the linear filtration law_ , Moscow, Gas business (Russian), No. 6, 1970
  * [14] Peaceman, D.W.  _Interpretation of Well-Block Pressures in Numerical Reservoir Simulation_ , SPEJ, 183-94, June 1978; Trans., AIME, 253. Paper SPE 6893
  * [15] van Poolen, H.K., Breitenbach, E.A., Thurnau, D.H  _Treatment of Individual Wells and Grids in Reservoir Modeling_ , SPEJ (Dec. 1968) 341-346. Trans. AIME 243.
  * [16] Peaceman, D.W.  _Interpretation of Well-Block Pressures in Numerical Reservoir Simulation with Nonsquare Grid Blocks and Anisotropic Permeability_ , SPEJ, June, pp. 531-543, 1983; Trans., AIME, 275. Paper SPE 10528 presented at the 1982 SPE Symposium on Reservoir Simulation, New Orleans, Jan 31-Feb 3.
  * [17] Alvestad, J., Holing, K., Christoffersen, K., Langeland, O., Stave, O.  _Interactive modeling of multiphase inflow performance of horizontal and highly deviated wells_ , Paper SPE 27577 presented at the European Petroleum Computer Conference, Aberdeen, UK, 15–17 March 1994.
  * [18] Mochizuki, S.  _Well productivity for arbitrarily inclined well_ , Paper SPE 29133 presented at the SPE Reservoir Simulation Symposium, San Antonio, TX, 12–13 February 1995.
  * [19] Klausen, R.A., Aavatsmark, I.  _Connection transmissibility factors in reservoir simulation for slanted wells in 3D grids_ , Paper presented at the 7th European Conference on the Mathematics of Oil Recovery, Baveno, Italy, 5–8 September 2000.
  * [20] Ding, Y., Renard, G., Weill, L.  _Representation of wells in numerical reservoir simulation_ , SPE Reservoir Evaluation Engrg. (1998) 18–23.
  * [21] Su, H.J.  _Modelling of off-center wells in reservoir simulation_ , SPE Reservoir Engrg. (1995) 47–51.
  * [22] Peaceman, D.W.  _Representation of a Horizontal Well in Numerical Reservoir Simulation_ , Paper SPE 21217, 1991. Advanced Technology Series, vol 1, No.1 (1993) 7–16.
  * [23] Babu, D.K., Odeh, A.S.  _Numerical simulation of horizontal well_ , Paper SPE 20161 presented at the SPE Middle East Oil Show, Bahrain, 16–19 November 1991.
  * [24] Babu, D.K., Odeh, A.S.  _Productivity of a Horizontal Well_ , SPERE, 417-21, Nov. 1989. Paper SPE 18298 presented at 1988 SPE Annual Technical Conference and Exhibition, Houston, Oct 2-5.
  * [25] Palagi, C.L., Aziz, K.  _The modeling of vertical and horizontal wells with Voronoi grid_ , Paper SPE 24072 presented at the SPE Western Regional Meeting, Bakersfield, CA, 30 March–1 April 1992.
  * [26] Lee, S.H., Milliken, W.J.  _The Productivity Index of an Inclined Well in Finite-Difference Reservoir Simulation_ , Paper SPE 25247 presented at the SPE Symposium on Reservoir Simulation, New Orleans, Feb. 28-Mar. 3, 1993.
  * [27] Ding, Y.  _A generalized 3D well model for reservoir simulation_ , SPE Journal, 1 (1996) 437–450.
  * [28] Ding, Y., Jeannin, L.  _A multi-point flux approximation scheme for well modelling in reservoir simulations_ , Paper presented at the 7th European Conference on the Mathematics of Oil Recovery, Baveno, Italy, 5–8 September 2000.
  * [29] Morita, N., Singh, S.P., Chen, H.S., Whitfill, D.L.  _Three-dimensional well model pre-processors for reservoir simulation with horizontal and curved inclined wells_ , Paper SPE paper 20718 presented at the SPE Annual Technical Conference and Exhibition, New Orleans, LA, 23–26 September 1990.
  * [30] Fung, L.S.-K., Buchanan, W.L., Sharma, R.  _Hybrid-CVFE method for flexible grid reservoir simulation_ , Paper SPE 25266 presented at the SPE Reservoir Simulation Symposium, New Orleans, LA, 28 February–3 March 1993.
  * [31] Goode, P.A., Thambynayagam, R.K.M.  _Pressure Drawdown and Buildup Analysis of Horizontal Wells in Anisotropic Media_ , SPEFE, Dec. 1987, pp. 683-97.
  * [32] Kuchuk, F.J., Goode, P.A., Brice, B.W., Sherrared, D.W., Thambynayagam, R.K.M..  _Pressure Transient Analysis and Inflow Performance for Horizontal Wells_ , Paper SPE 18300 presented at the SPE Annual Technical Conference and Exhibition, Houston, Oct. 2-5, 1988.
  * [33] Clonts, M.D., Ramey, H.J., Jr.  _Pressure Transient Analysis for Wells with Horizontal Drainholes_ , Paper SPE 15116 presented at the SPE California Regional Meeting, Oakland, Apr. 2-4, 1986.
  * [34] Daviau, F., Mouronval, G., Bourdarot, G., Curutchet, P.  _Pressure Analysis for Horizontal Wells_ , Paper SPE 14251 presented at the SPE Annual Technical Conference and Exhibition, Las Vegas, Sept. 22-25, 1985.
  * [35] Ozkan, E., Raghavan, R., Joshi, S.D.  _Horizontal Well Pressure Analysis_ , SPEFE, 567-75, Dec. 1989.
  * [36] Economides, M.J., Brand, C.W., Frick, T.P  _Well Configurations in Anisotropic Reservoirs_ , SPEFE 11, pp. 257-262, Dec. 1996.
  * [37] Maizeret, P.D.  _Well Indices for Nonconventional Wells_ , Master’s report, Stanford University, 1996.
  * [38] Ouyang, L.B., Aziz, K.  _A General Single-Phase Wellbore/Reservoir Coupling Model for Multilateral Wells_ , SPEREE 4(4), pp. 327-335, Aug. 2001.
  * [39] Ouyang, L.B.  _Single Phase and Multiphase Fluid Flow in Horizontal Wells_ , PhD thesis, Stanford University, 1998.
  * [40] Basquet, R., Alabert, F.G., Caltagirone, J.P., Batsale, J.C.  _A Semi-Analytical Approach for Productivity Evaluation of Wells with Complex Geometry in Multilayered Reservoirs_ , Paper SPE 49232 presented at the SPE Annual Technical Conference and Exhibition, New Orleans, Sept. 27-30, 1998.
  * [41] _Eclipse Reference Manual, 96A Release_ , GeoQuest, Schlumberger, 1996.
  * [42] _tNavigator User Manual, 21.1 Release_ , RFD,2021.
  * [43] Naz R., Alsaedi A., Hayat T.  _Flow of fourth grade fluid in a porous medium_ , Applied and Computational Mathematics, an International Journal, Volume 14,2 125-140.
  * [44] _Computer Modeling Group. STARS user’s guide, version 2015_ , Calgary, Alberta, Canada, 2015.