Mathematics of fluid flow

The purpose of this page is to review the mathematics of fluid flow. We limit our review to essential aspects of partial differential equations, vector analysis, numerical methods, matrices, and linear algebra. These topics are discussed in the context of two fluid flow applications: analysis of the convection/dispersion equation and diagonalization of the permeability tensor. For more details about the mathematics presented here, consult the literature.

Partial differential equations (PDEs) are frequently encountered in petroleum engineering. We review basic concepts of PDEs by considering the relevant mathematical properties of the continuity equation. (See also this refresher on differential calculus.)

Continuity equation

Fluid flow through a volume can be described mathematically by the continuity equation. The continuity equation has many uses, and its derivation is provided to illustrate the construction of a partial differential equation from physical reasoning. We begin by considering the flow illustrated in Fig. 1. The block in Fig. 1 has length (Δx), width (Δy), and depth (Δz). Fluid flux (J) is the rate of flow of mass per unit cross-sectional area normal to the direction of flow. The notation (J_x)_x denotes fluid flux in the x direction at location x. The cross-sectional area perpendicular to the flux direction is Δy Δz. Fluid flows into the block at x with fluid flux J_x and out of the block at x + Δx with fluid flux J_x + Δx. Applying the principle of conservation of mass, we have the mass balance, which is written as

....................(1)

The mass entering the block in time interval, Δt, for flux across the face of the block at x is

....................(2)

The mass leaving the block in time interval, Δt, across the face of the block at x + Δx is

....................(3)

where we have added a source term q. Flow out of the block through q is represented by q > 0, and flow into the block is represented by q < 0. The source term, q, can represent a variety of important physical systems, including wells, aquifer support, fluid flow into a fracture system from matrix blocks in a naturally fractured reservoir, and gas flow into a cleat system from the coal in a coalbed.

• 

  Fig. 1 – Coordinate system for continuity equation.

Mass accumulation in the block is the change in concentration C of the mass in the block during the time interval Δt, where concentration, C, is defined as the total mass in the block divided by the block volume. The mass accumulation term is

....................(4)

where concentration is evaluated at times t and t + Δt.

Substituting Eqs. 2 through 4 into Eq. 1, dividing by Δx Δy Δz Δt, and rearranging gives

....................(5)

In the limit as Δx → 0, Δt → 0, the differences in Eq. 5 are replaced by partial derivatives. We assume the fluxes and concentrations are sufficiently smooth and continuous to allow the replacement of differences by partial derivatives. Eq. 5 becomes the continuity equation in one space dimension.

....................(6)

Eq. 6 is an example of a partial differential equation.

Partial differential equations

PDEs are an extension of the concept of ordinary differential equations (ODEs). Unlike an ODE, which depends on only one independent variable, a PDE depends on two or more independent variables. In the previous example, Eq. 6 depends on two independent variables: one space dimension (x) and time (t). The order of Eq. 6 depends on the form of concentration and flux. The order of a PDE is the order of the highest derivative that appears in the equation.

....................(7)

for a function (ψ) of two or more independent variables {x,y, ...}. A PDE is linear if it is first order in the unknown function and its partial derivatives, and the coefficients of the partial derivatives, are either constant or depend on the independent variables {x,y, ...}. We illustrate these concepts by considering the continuity equation for flow of a fluid with density (ρ), velocity (v_x), and no source or sink terms. The concentration, C, and flux, J_x, for this example are

and

....................(8)

Substituting Eq. 8 into Eq. 6 gives

....................(9)

Eq. 9 is a linear, first-order PDE if density is the unknown function and velocity is constant. The situation is not so simple in more physically realistic systems.

Consider, for example, a slightly compressible fluid in which density is given by

....................(10)

where P is pressure, c_f is fluid compressibility, and the subscript, 0, refers to a reference value of pressure. Assume, as well, that velocity is proportional to pressure gradient so that

....................(11)

where α is the proportionality constant. Substituting Eqs. 10 and 11 into Eq. 9 gives

....................(12)

Eq. 12 is a nonlinear, second-order PDE. It is second order because of the second-order partial derivative of pressure with respect to x, and it is nonlinear because of the square of the pressure gradient term.

Solutions of PDEs depend on the form of the PDEs and their associated boundary conditions. An important class of second-order PDEs has the form

....................(13)

where the functions {A, B, C, G} are known functions of two independent variables {x,y} and the first-order partial derivatives ∂ ψ(x,y)/∂ x and ∂ ψ(x,y)/∂ y in a region (R) bounded by a surface (S). The mathematical properties of the second-order PDEs depend on the relationship between the functions {A, B, C, G}. A classification scheme for second-order PDEs is given in Table 1.

• 

  Table 1

Boundary conditions for second-order PDEs may be written as

....................(14)

where ψ(x,y) is the unknown function of two independent variables {x,y}, and ∂ ψ(x,y)/∂ n is the derivative normal to a boundary. The functions {α, β, γ} are known functions of {x,y}. All of the functions and applicable derivatives are defined in a domain (R) bounded by a surface (S). A classification scheme for the boundary conditions of a second-order PDE is given in Table 2. The boundary conditions associated with the examples in Table 1 are given in Table 3. The significance of PDE classification is considered further in the discussion of the convection/dispersion equation presented next.

• 

  Table 2

• 

  Table 3

One-dimensional (1D) convection/dispersion equation

The continuity equation is used to describe the mixing of one substance with another by writing flux in the form

....................(15)

The concentration, C, is the concentration of the solute in the solvent. The term v is the velocity of the solute, and D is the dispersion of the solute into the solvent. Substituting Eq. 15 into Eq. 6, the 1D convection/dispersion equation without sources or sinks, gives

....................(16)

If we assume that v and D are constant, Eq. 16 simplifies to the form

....................(17)

Eq. 17 is the 1D convection/dispersion (C/D) equation. The dispersion term is D∂²C/∂x², and the convection term is –v∂C/∂x. If the dispersion term is much larger than the convection term, the solution of Eq. 17 can be approximated by the solution of the equation

....................(18)

Eq. 18 is a parabolic PDE and behaves mathematically like a heat conduction equation. If the convection term is much larger than the dispersion term, the solution of Eq. 17 can be approximated by the solution of the equation.

....................(19)

Eq. 19 is a first-order hyperbolic PDE.

A solution of the 1D C/D equation, given in Eq. 17, is obtained as follows. We assume that solute is moving in the x-direction with constant speed (v_x). The concentration C(x,t) is a function of space and time. We must specify two boundary conditions and an initial condition for the concentration C(x,t) to solve the C/D equation. We impose the boundary conditions, C(0,t) = 1 and C(∞,t) = 0, for all time, t > 0, and the initial condition, C (x,0) = 0, for all values of x > 0. The boundary condition, C (0,t) = 1, says that we are injecting 100% solute at x = 0, and the boundary condition, C(∞,t) = 0, says that the solute never reaches the end of the flow path at x = ∞. The initial condition, C(x,0) = 0, says that there is no solute in the solvent at the initial time, t = 0. The solution of the C/D equation is

....................(20)

where the complementary error function erfc(y) is defined 1 as

....................(21)

The integral in Eq. 21 can be solved using the series expansion on the right side of Eq. 21 or a numerical algorithm. Eq. 20 is illustrated in Fig. 2 for physical parameters v = 1 ft/day and D = 0.01ft²/day. This solution is used to evaluate a numerical solution of the C/D equation. (Look at the page: Numerical methods analysis of fluid flow)

• 

  Fig. 2 – Analytical solution of the 1D C/D equation.

Nomenclature

A, B, C, G	=	functions, Eq. 13
cf	=	fluid compressibility, Eq. 10
C	=	concentration, Eq. 4
D	=	dispersion of the solute into solvent, Eq. 15
Jx, Jy, Jz	=	fluid flux in x-, y-, z-directions
(Jx)x	=	fluid flux in x-direction at location x
(Jy)y	=	fluid flux in y-direction at location y
(Jz)z	=	fluid flux in z-direction at location z
P	=	pressure
q	=	source term, Eq. 3
R	=	region
S	=	surface
t	=	time
tn	=	present time
tn+1	=	future time
v	=	velocity of solute, Eq. 15
vx	=	speed in x-direction, Eq. 9
x,y,z	=	space dimensions
α	=	proportionality constant, Eq. 11
{ α, β, γ }	=	functions, Eq. 14
Δt	=	time interval
Δx	=	length
Δy	=	width
Δz	=	thickness
ρ	=	density, Eq. 8
ψ	=	function, Eq. 7

Noteworthy papers in OnePetro

Matthews, C. S., & Russell. D. G. 1967. Pressure Buildup and Flow Tests in Wells, SPE Monograph Series Vol. 1, Society of Petroleum Engineers, Richardson, TX, 163 pp. SPE Bookstore or Worldcat

External links

Use this section to provide links to relevant material on websites other than PetroWiki and OnePetro

See also

Numerical methods analysis of fluid flow

Vector analysis of fluid flow

Diagonalizing the permeability tensor

PEH:Mathematics_of_Fluid_Flow