Open access peer-reviewed chapter - ONLINE FIRST

# Porous Flow with Diffuse Interfaces

By Paul Papatzacos

Submitted: September 7th 2020Reviewed: December 12th 2020Published: January 19th 2021

DOI: 10.5772/intechopen.95474

## Abstract

This chapter presents a model developed by the author, in publications dated from 2002 to 2016, on flow in porous media assuming diffuse interfaces. It contains five sections. Section 1 is an Introduction, tracing the origin of the diffuse interface formalism. Section 1 also presents the traditional compositional model, pointing out its emphasis on phases and questioning the concept of relative permeabilities. Section 2 presents the mass, momentum, and energy balance equations, for a multicomponent continuous fluid, in their most general form, at the pore level. The existence of constitutive equations with phase-inducing terms is mentioned, but the equations are not introduced at this level, and phases are not an explicit concern. Section 3 is about the averaging of the pore level equations inside a region containing many pores. There is no explicit mention of phases and therefore not of relative permeabilities. Section 4 is the technical basis from which the constitutive equations of the model arise, and it is shown that many models can exist. Section 5 introduces constitutive equations and presents a minimal model for multicomponent, multiphase, and thermal flow in neutrally wetting porous media, i.e., a model with a minimal amount of phenomenological parameters.

### Keywords

• flow in porous media
• Marle averaging
• diffuse interface
• multiphase flow
• phase segregation
• relative permeabilities

## 1. Introduction

The model presented in this chapter was developed in Refs. [1-3]. Ref. [1] covers the one-component two-phase case, Ref. [2] is a generalisation to an arbitrary number of components and three phases (two liquids, one gas), Ref. [3] generalises to variable temperature. All three show applications to neutrally wetting media. The way wetting can be accounted for is discussed in general terms in Ref. [1]. A partial practical implementation, valid for incomplete wetting, is suggested in Ref. [4] but a fully satisfactory solution that accounts for total wetting (or capillary condensation in adsorption terms) is still a matter for further research.

The diffuse interface theory was initiated in 1893 in a paper by van der Waals (see the translation by Rowlinson [5]) where he proposes to replace the old assumption of a surface between phases by the assumption of a continuous transition inside a thin interphase region, where certain quantities, notably the density, vary continuously. The core of his theory consists of a Helmholtz function modified by the addition of a term proportional to the squared gradient of the density, thus accounting for the energy stored in the region. In the model presented in the following pages, the van der Waals theory is introduced at the upscaled level. See Section 4. Generalising the van der Waals expression to νchemical components [2], one obtains

FF=FbF+α=1ν12ΛαRα2,E1

where Rαis the upscaled density of fluid component α, Fis the Helmholtz function per unit volume, the superscript F referring to the fluid, while bF means “bulk fluid” and refers to fluid regions that are far away from interphase layers. The Λαare constants (Assumption A7, Appendix 1). For later reference FbFand its differential are

FbF=PbF+α=1νCbFαRα,dFbF=SbFdT+α=1νCbFαdRα,E2

where PbFis the pressure in the bulk fluid, SbFthe entropy per unit volume, and CbFαthe chemical potential of component α, divided by its molar mass (the Cα/Mαof Refs. [2, 3]).

The van der Waals paper was followed, in 1901, by a paper by Korteweg [6] about the equations of motion of a fluid with large but continuous density changes, where Korteweg showed that, for such a fluid, the usual scalar pressure must be replaced by a symmetric second order tensor.

The van der Waals-Korteweg papers were apparently forgotten, then rediscovered in the nineteen seventies, the diffuse-interface method being introduced as a novel way to solve fluid mechanics problems in two-phase flow. For a review, see [7]. While van der Waals and Korteweg assumed that the gradients in Eq. (1) are small, modern advances have shown that this limitation can be lifted [7, 8].

Ref. [8] is, to the author’s knowledge, the first formulation of flow in porous media with diffuse interfaces.1 The purpose with such a novel formulation is to avoid some of the known weaknesses of the traditional compositional models treating multiple phase problems. The mathematical core of the compositional models used in reservoir engineering, consists of equations expressing mass balance for each chemical component. The distribution of the phases in the reservoir is essential to the formulation, and is determined at the beginning of each time step. A component, present in a phase, is transported with the Darcy velocity of the phase. Each phase-dependent Darcy velocity is written with a permeability that is modified by a multiplicative factor, the relative permeability. The above description emphasises two of the weaknesses of the traditional models. The first is the assumption that well-defined phases exist at all times. The second is the use of relative permeabilities, a concept that is “seriously questioned”, as expressed by Adler and Brenner in a 1988 paper [9]. Concerning the understanding of relative permeabilities in the framework of the model presented here, see the comments following Eq. (100) below.

The mathematical core of the model (actually a family of models) to be presented consists of mass balance equations, one per chemical component, of a momentum balance equation, and of an entropy balance equation. The thermodynamical description of the fluid mixture involved is part of the core. The central purpose is to calculate the component densities, and other characteristic quantities such as fluid velocity and fluid temperature, as functions of space and time. If the approximation of constant temperature is valid, only the mass and momentum balance equations are necessary.

Phases (and thereby relative permeabilities) do not take part in the formulation. They result from the solutions of the model equations, and are detected by rapid variations of densities, and by regions of approximately uniform densities. They can be shown to exist in static equilibrium or steady state dynamical situations [1, 2, 3].

A minimal model is presented in Section 5, consisting of a minimal amount of parameters.

A note on the appendices: Some concepts are grouped in appendices for easy reference. Appendix 1, for example, lists all the assumptions the model is built on.

A note on wetting: The problem of accounting for the wetting properties of the pore surface remains to be solved. There are two approaches to the problem: through boundary conditions to the Navier–Stokes equations at pore level, or through the theory of adsorption at the upscaled Darcy level. The first approach has been used in the publications considered here, and it is explained in Ref. [10] that the diffuse interface theory presented in Ref. [1] is consistent for neutral wetting, i.e., for pore level wetting angles around 90: see Appendix 1, Assumption A0. Refs. [2, 3] assume this limitation.

A note on notation: Right-handed Cartesian coordinates x1x2x3are assumed, the plane x1x2being horizontal, and axis x3pointing upwards. Any vector Ahas components Akwhere kis 1, 2, or 3. In addition, the notation t=/tand k=/xkis used. The summation convention applies to latin indexes i, j, k, l: an index that is repeated in a term (as in AjBj, or Ckk) indicates summation (j=13AjBj, or k=13Ckk). Note that an index that is repeated, but not in the same term (as in Ai=Bi), means that the expression is valid for all values of the index in the set 1,2,3. Symbols are otherwise defined when introduced.

## 2. Pore level equations

The fluid is a mixture of νchemical components, it is continuous, and if phases exist, there are interphase regions where quantities vary continuously, possibly rapidly. With assumptions A1 to A3 (Appendix 1) the balance equations for mass, momentum, and energy, are written below in a most general manner:

tρα+kραvk+ikα=0,α=1ν,E3
tρvi+kρvivktki=ρfi,E4
tε+kεvk+jktkivi=ρfkvk.E5

For the theoretical basis of these equations see chapter 11 of the book by Hirschfelder et al [11]. Greek superscripts indicate the species so that ραis the mass density of component α, while ρis the total density:

ρ=α=1νρα.E6

Further, vis the local velocity, defined as the total momentum divided by the total mass, both per volume; ιαis the non-convective mass current of component α, with the property

α=1νια=0,E7

Jis the non-convective energy current, and tijis the stress tensor; ε=u+12ρv2is the total energy, uis the internal energy, both per unit volume, and fis an external force per unit mass of fluid, thus being species-independent. In the case of gravity

f=g,g=00g,g=9.81m/s2.E8

It is convenient, for later use, to introduce the gravitational potential

W=gx,wherebyfi=iW.E9

The symmetry tij=tjiis assumed below. It is reminded that it expresses angular momentum balance for a fluid where no other torques exist than the one due to the external force f, and the one due the surface stress ntacting on every surface element where the normal vector is n.

Within the van der Waals theory, one expects that ιkα, Jk, and tki, contain terms whose magnitudes are important in the interphase regions, but are otherwise negligible. No constitutive equations are introduced at the pore level, but it is mentioned for later reference, that for a simple one-phase fluid,

tij=pδij+θij,E10

where pis the pressure and θijis the viscous stress tensor (symmetric and linear in the gradients of the vi, with coefficients of shear and bulk viscosity).

The upscaling, i.e., the averaging over many pores, is done in the next section by the method due to Marle [12]. This method assumes that the physical quantities that appear in the balance equations above are treated as distributions [12, 13], the underlying reason being that such quantities are discontinuous, and that one needs to average their partial derivatives. Taking ρas an example, it is (i) a fluid density in a pore, (ii) a rock density in the rock matrix, (iii) undefined on the pore surface, and one needs to average tρand kρ.

The physical quantities appearing in Eqs. (3) to (5) are listed in the first line of Table 1. Their values in the pores, or in the rock, are denoted with a superscript F, or S, as shown in the second and third lines of the table. In the third line, a missing entry indicates non-existence, the first two 0-values indicate no material transport in the solid, the third zero value follows from Assumption A4 (Appendix 1): as shown by Marle [12], this assumption implies that the momentum balance equations in the solid and on the pore surface have the form 0=0. It is important to keep in mind, especially when averaging, that any quantity with an F superscript is equal to 0on the pore surface (denoted by Σ) and in the solid; likewise, any quantity with superscript S is equal to 0on Σand in the fluid. The quantities in the first line of the table are not defined on Σ, but limit values are, as in

Genericραρvkιkαtkifkujk
In fluidρFαρFvkFιkFαtkifkuFjkF
In solidρS000uSjkS

### Table 1.

First line, left of the vertical: quantities appearing in Eqs. (3) and (4); right of the vertical: quantities only appearing in Eq. (5). Second and third lines: notation when specialising to the fluid and the solid. Concerning the missing entries and the three 0-values, see beginning of paragraph containing Eq. (11).

vFΣ=0,E11

expressing that the fluid velocity vanishes on Σ. Expressions of the type XΣfor some X, used in this section and in the next, define the limit of the quantity Xat Σ, along the line carrying the normal to Σpointing towards the fluid. The existence of the normal implies some idealisation of the pore surface.

Note finally that, according to Eq. (6),

ρF=α=1νρFα.E12

As a preliminary to upscaling, Eqs. (3)(5) are now explicitly written in terms of distributions, using the notation of Appendix 2 where Ais replaced by Fαor F, as appropriate, while Bis replaced by S. As an example, ραxt, is a distribution depending on space xand time t. It is continuous in time but not in space, being equal to ρFαxtif xis inside a pore and to ρSxtif xis in the rock. It is not defined when xis on Σbut we assume that ρFαΣand ρSΣexist.

The generalised mass balance equation is obtained from Eq. (3), using Eqs. (112), (113), and (11):

tρα+kραvk+ikα+ιαnΣδΣ=0,α=1ν.E13

From this equation one now gets, using Table 1, the equations that are separately valid inside the pores, on Σ, and in the solid:

3tρFα+kρFαvkF+ikFα=0,α=1ν,inFE14
ιFαnΣδΣ=0,      α=1ν,onΣE15
tρS=0.inSE16

Turning to the generalised momentum balance equation, one must account for Assumption A4 (Appendix 1), implying that the momentum balance equations in the solid and on the pore surface have the form 0=0. Inside the pores, one simply re-writes Eq. (4), with the superscript F on the density and velocity:

tρFviF+kρFvkFviFtki=ρFfi,inF.E17

The generalised energy balance equation is obtained from Eq. (5), using Eqs. (112), (113), and (11):

tu+12ρv2+kuvk+12ρv2vk+jktkivi+JFnΣJSnΣδΣ=fkρvk.E18

Using Table 1, one then obtains:

tuF+12ρFvF2+kuFvkF+12ρFvF2vkF+JkFtkiviF=fkρFvkF,inFE19
JFnΣJSnΣδΣ=0,onΣE20
tuS+kJkS=0,inS.E21

## 3. Averaging

The Marle averaging process [12] is followed in all essentials, except in the assumption that well-defined phases exist, separated by interphase surfaces. The averaging volume is a sphere of radius r, large when compared to a pore radius, small when compared to a linear dimension of the reservoir. A Cfunction mxis introduced, somewhat flat around x=0and equal to zero for xr, normalised so that its integral over all space is equal to 1. (See Eq. (17) in Ref. [12] for an example of such a function.) Given any function of space and time, fxt, its average Fxtis obtained by the convolution

Fxtfmxt=R3fytmxydy,E22

where the integration is over all of space. The convolution ensures that Fis C[12].

The averaged balance equations are differential equations in the averaged quantities (averaged densities, velocities, ...). These equations are established by a three-step process. Step 1: the generalised equations for mass, momentum, and energy balance are each in turn convoluted with m. Step 2: the following rules [12] are applied, allowing to take the differential operators out of the averaging convolutions:

tfZm=tfZm,E23
ifZm=ifZmεZfZΣniδΣm,E24

where fZfZxtand Zis F, Fα, or S; εZis 1 if Zis F or Fα, −1 if Zis S. Eq. (11) is also applied at step 2. Step 3: the remaining convolutions are used to define averaged quantities, where the following constraints should be obeyed: (i) except for the definition of porosity, the way to define the averages is suggested by the equations obtained after completion of step 2; (ii) the averaged equations have essentially the same forms as Eqs. (3) to (5).

Differential equations in the averaged quantities result from the three steps. It is shown below that the mass and momentum balance equations should be treated together, and that the energy equation can be treated as an addition.

The averaging of the mass and momentum balance equations follows. Steps 1 and 2 are applied to Eqs. (14) to (17), and lead to:

tρFαm+kρFαvkF+ιkFαm=ιFαnΣδΣm,E25
ιFαnΣδΣm=0,E26
tρSm=0,E27
tρFviFm+kρFvkFviFtkim=nktkiΣδΣm+fiρFm.E28

Note that the last term on the right-hand side of the last equation originally is ρFfim, but fican be taken out of the convolution integral because of Assumption A5 (Appendix 2). The definitions of averaged (or upscaled) quantities now follow (step 3). Porosity Φis defined first, then follow definitions suggested by the convolution operations in the equations above.

Porosity Let a function χxbe 1 when xis in F, and let it be 0 otherwise.

Then

χm=Φ.E29

Note that, according to definition (22), Φcan depend on x.

Species adsorption The left-hand side of Eq. (26) defines the amount, KΣα, of component αadsorbed at the pore surface:

ιFαnΣδΣm=KΣα.E30

Density of solidEq. (27) suggests defining the solid density RSby:

ρSm=1ΦRS.E31

Density of componentαThe first term on the left-hand side of Eq. (25) suggests defining the density Rαof fluid component αby

ρFαm=Rα.E32

Note that Eq. (12) implies that the averaged total fluid density, R, is

R=α=1νRα=ρFm.E33

Fluid velocity The first term on the lef-hand side of Eq. (28) has the convolution of a product of two term, where the average of one of them is known from Eqs. (32) and (33). The averaged fluid velocity, denoted Vi(without the F superscript since there is no velocity in S or on Σ) is then defined by

ρFviFm=ΦRVi.E34

Diffusive mass current of componentαThe second convolution on the left-hand side of Eq. (25) is used to define the upscaled diffusive current in the fluid, denoted Ikα(without the F superscript since there are no diffusive currents in S or on Σ). It is not defined as the average of ιkFαbecause both pore level effects of convection and diffusion contribute to it [12]. Keeping in mind the constraint that the averaged equations should have the same form as the original ones, Ikαis defined by

ρFαvkF+ιkFαm=ΦRαVk+ΦIkα.E35

Note that, summing this equation over αfrom 1 to ν, and using previous results, one gets

α=1νIkα=0.E36

Stress tensor The second convolution on the left-hand side of Eq. (28) suggests defining Tki, the upscaled version of tki, by

ρFvkFviFtkim=ΦRViVkTki.E37

Note that tki=tikimplies Tki=Tik.

Frictional force per unit volume The first convolution on the right-hand side of Eq. (28) suggests defining the frictional force per unit volume FiFby

nktkiΣδΣm=FiF.E38

The upscaled mass and momentum balance equations now follow, in the following order: mass balance for the solid, mass balance at the pore surface, mass balance for the fluid in the pores, and momentum balance for the fluid in the pores:

t1ΦRS=0,E39
KΣα=0,α=1ν,E40
tΦRα+kΦRαVk+ΦIkα=0,α=1ν,E41
tΦRVi+kΦRVkViTki=FiF+ΦRfi.E42

The first equation states that porosity and solid density do not vary with time (consistently with Assumption 4 (Appendix 2) but can vary in space. The second equation states that adsorption is negligibly small, for any component, consistent with Assumption A0 (Appendix 2). The remaining two equations determine the νcomponent densities and the three velocity components. This means that the components of the diffusive mass current, of the stress tensor, and of the frictional force must be provided. (Constitutive equations).

The averaging of the energy balance equations now follows. Steps 1 and 2 are applied to Eqs. (19) to (21), giving:

tuF+12ρFvF2m+kuFvkF+12ρFvF2vkF+JkFtkiviFm=JFnΣδΣm+fkρFvkFm,E43
JFnΣδΣmJSnΣδΣm=0,E44
tuSm+kJkSm=JSnΣδΣm.E45

As in the case of Eq. (28), Assumption A4 (Appendix 2) has been used to take fkout of the convolution on the right-hand side of Eq: (43). Six definitions are introduced below, built on the definitions that were introduced in connection with the averaging of the mass and momentum balance equations. An underscore indicates the defined quantity.

Internal energy per unit volume of solid

uSm=1ΦUS¯,E46

Internal energy current of solid

jkSm=JkS¯,E47

Solid to fluid energy transfer

jSnΣδΣm=QSF¯,E48

Fluid to solid energy transfer

jFnΣδΣm=QFS¯,E49

Internal energy per unit volume of fluid

uF+12ρFvF2m=ΦUF¯+12ΦRV2,E50

Internal energy current of fluid

uFvkF+12ρFvF2vkF+jkFtkiviFm=ΦUFVk+12ΦRV2Vk+JkF¯TkiVi.E51

Using these definitions in Eqs. (43) to (45) one obtains:

tΦUF+12ΦRV2+kΦUFVk+12ΦRV2Vk+JkFTkiVi=QFS+ΦRVkfk.inFE52
QFS+QSF=0.onΣE53
t1ΦUS+kJkS=QSF.inS.E54

Eq. (52) contains a redundancy in the form of a balance equation for kinetic energy. This equation can be obtained directly by multiplying both sides of Eq. (42) with Vi, summing over i, and using Eq. (41). Subtracting the equation thus obtained from Eq. (52), one gets

tΦUF+kΦUFVk+JkF=QFSFkFVk+TjijVi.E55

## 4. Basis for constitutive equations

In practical application, the upscaled balance equations will be used to calculate, primarily, the densities Rα, the velocity components Viand the temperature T. Assuming that the approximation of constant Tis valid, one needs only focus on the mass and momentum balance equations, (39) to (42). One sees in this case that one needs expressions for the Ikα, the Tij, and the Fiin terms of Tand the Rα. If the approximation of constant temperature is not valid, one must first distinguish between the temperatures of the solid and the fluid, and one needs an equation for the energy transfer in case of a temperature difference. One introduces the simplifying Assumption A6 (see Ref. [12] and Appendix 1), from which it follows that just one additional equation is needed for calculating Txt. Such an equation usually describes the evolution of either total energy or total entropy. In either case, expressions are needed for the currents JkSand JkF. (Expressions for the energy transfers QSFand QFSare unnecessary since they cancel when taking the sum of the solid and fluid energies). Most of what is needed is obtained in Section 5 by applying the theory of irreversible processes, starting from the evolution equation for entropy, although it seems that some preliminary work is unavoidable to directly obtain an expression for the pressure tensor Pijthat replaces the usual scalar pressure.

The derivation of the pressure tensor is given below, followed by the derivation of the evolution equation for the total entropy. It is essential to use expression (1) in both derivations. 2

### 4.1 The pressure tensor

One considers the upscaled fluid, consisting of a mixture of νcomponents in a container with surface ∂Ωand volume Ω, inside a large bath at uniform and constant temperature T. One looks for conditions of equilibrium in the presence of gravity. The fluid has the Helmholtz free energy density given by Eq. (1), and it is assumed that the bounding surface is neutrally wetting so that there is no energy stored on ∂Ω. The total energy stored in the fluid is

F=ΩFF+Wα=1νRα,E56

where Wis the gravitational potential (see Eqs. (8) and (9)).

One now looks for the conditions the Rαsatisfy when Fis at its minimum, given that the total mass, ΩR, is constant. This is minimisation with constraint, a standard problem in variational calculus. It is easily found that

FbFRα+WκΛα2Rα=0,inΩ,E57
nRα=0,on∂Ω,E58

where κis a Lagrange multiplier, possibly a function of Tbut not of the Rα. See Ref. [14] for the technique and the theorems involved: essentially, the expression on the left-hand side of Eq. (57) must be continuous. Note that Eq. (58) expresses the non-wetting property of the outer boundary: the density neither increases nor decreases along the normal.

The next step consists in multiplying Eq. (57) with iRαand summing over α. Each term of the resulting equation can then be re-written as a gradient (first term), or as a sum of a gradient and the component of a force (sum of second and third terms), or as a sum of a gradient and a divergence (fourth term). The result is the following expression:

kPikfi=0.E59

where

Pik=PbFα=1νΛαRα2Rαα=1ν12ΛαRα2δik+α=1νΛαiRαkRα.E60

The first equation obviously generalises the classical P=f, and it is natural to generalise Eq. (10) by setting

Tij=ΦPij+Θij.E61

Tijand Pijare symmetric, so that Θij=Θji. Eqs. (60) and (61) are used in the long calculation that lead to the evolution equation for the entropy. They are also used to get expressions (99) and (100) for the modified Darcy equation, where the presence of the Λ‘s is exclusively due to their presence in Eq. (60).

### 4.2 The entropy evolution equation

Eqs. (1) and (2), together with Eqs. (41), (54), (55), (60), and (61), are needed to get the entropy equation for the fluid.

Taking the differential of Eq. (1) and using Eq. (2), one obtains, keeping in mind the assumption that the Λαare constants,

dFF=SbFdT+α=1νCbFαdRα+α=1νΛαiRαdiRα.E62

This equation leads immediately to two conclusions. There is no additional entropy, and no additional chemical potentials due to large density gradients, since: (i) SF, defined as FF/T, is equal to SbF; (ii) CFα, defined as FF/Rα, is equal to CbFα. The differential of FFcan now be re-written with a simpler notation:

dFF=SFdT+α=1νCFαdRα+α=1νΛαiRαdiRα.E63

To obtain dSFin terms of dUFand dFFone differentiates UF=TSF+FF, obtaining

TdSF=dUFα=1νCFαdRαα=1νΛαkRαdkRα.E64

This expression is now used to construct an evolution equation for the total entropy (including fluid and solid). Such an equation must be of the form

tΦSF+1ΦSS+kΦSFVk+Pk=Q,E65

where SSis the entropy of the rock, Pkis a diffusive current, and Qis a source term.

Eq. (64) directly gives the two equations that follow, by replacing dwith t(first equation), then with Vii(second equation):

TtSF=tUFα=1νCFαtRαα=1νΛαkRαtkRα.E66
TiSFVi=TSFiVi+ViiUFα=1νCFαViiRαViα=1νΛαkRαikRα.E67

Since it is ΦSFthat is needed to get Eq. (65), one needs to multiply the two equations above with Φand then commute Φwith tand i. The same commutations are necessary on the right-hand sides of the above equations, so as to get ΦUFand thus allow the use of Eq. (55). Assumption A8 (Appendix 2) has been added to Assumption A4 so as to avoid the proliferation of iΦ-terms in the minimal model to be presented. With this simplification in place, and remembering to replace Tijin Eq. (55) with the expression obtained from Eqs. (60) and (61), one obtains, after elementary but somewhat long calculations:

TtΦSF+iΦSFVi=kJkF+QSFFkFVk+ΘjijVi+Φα=1νCFαiIiα+Φα=1νΛαkRαkiIiα+Φα=1νΛαRαiVikkRα+kRαkRαiVi.E68

Keeping in mind that one is looking for an entropy equation of the form of Eq. (65), and in anticipation of using the methods of irreversible processes, one now considers each term, or group of terms, on the right-hand side above and writes it either as a divergence, or as the sum of a divergence and a scalar product. Most terms on the first line already have the required form. The general term in the first sum on the second line is easily transformed as required, the general term in the second sum also, although with somewhat more work. As to the third line, it is easily seen to be a sum of divergences. One then gets:

TtΦSF+iΦSFVi= kJkF+Φα=1νAαIkα+Φα=1νΛαkRαRαiVi+iIiα+QSFFkFVk+ΘjijViΦα=1νIiαiAα,E69

where Aα=CFαΛα2Rα. Note that Aαoccurs twice on the right-hand side of Eq. (69), once as Aα, and once as iAα, both times multiplying Ikαin a sum over α. Eq. (36) then implies that one can modify the above expression of Aαby the addition of any expression that does not depend on α. To conform with the notation of Ref. [2], and especially Ref. [3], one then sets

Aα=CFα+WΛα2Rα.E70

as the expression to substitute on the right-hand side of Eq. (69).

The entropy equation for the solid is much easier to obtain because of Assumption A4 (Appendix 2). Indeed, Eq. (64) is replaced by

TdSS=dUS,E71

which implies

Tt1ΦSS=t1ΦUS,E72

and, using (54):

Tt1ΦSS=kJkSQSF.E73

Taking the sum of Eqs. (69) and (73) one sees that QSFcancels, and that the right-hand side consists of a sum of scalar products and of a divergence, say kKk. A last step remains because the left-hand side has a multiplicative factor T(see Eqs. (69) and (73)). Dividing both sides by Tdoes not affect the scalar products but it modifies the divergence, at least if one assumes that Tis variable, to kKk/T. This is easily remedied by writing

1TkKk=kKkT+KkT2kT.E74

One finally obtains

tΦSF+1ΦSS+kΦSFVk+KkT=Q,E75

where

Kk=JkF+JkSΦα=1νΛαkRαRαiVi+iIiαΦα=1νAαIkα.E76
TQ=Kk/TkT+ΘjijViFkFVkΦα=1νIkαkAα.E77

QFSand QSFdo not appear in the entropy equation, and are thus not determinable inside the model. Similarly, the currents JkFand JkSonly appear as summed, so that they are not determined individually inside the model. (See the next section.)

## 5. Constitutive equations and the minimal model

The source term in the entropy equation plays a central role in what follows. It has been written as a sum of scalar products. Each term of this sum is the scalar product of a force (explicit or generalised) and a current. FkFis an explicit force, while the gradient of some quantity (temperature, velocity component, ...) is a generalised force. The theory of irreversible processes states that linear relations exist (at least for processes not far from equilibrium), between forces and currents. The coefficients, called phenomenological coefficients, are parameters whose signs must be such that the source term cannot be negative, to ensure against a decrease of the entropy when the system is isolated. The linear relations just mentioned are constitutive equations, and it is implied that the phenomenological coefficients must be provided as input. For various reasons, some coefficients can be put equal to zero, one often recurring reason being the belief that they are negligible in the physical situation considered. Thus there is a family of models, each member being characterised by its set of non-zero phenomenological coefficients.

The previously mentioned minimal model (see the text after Eq. (67)), is the model that contains the least possible number of non-zero phenomenological coefficients. “Least possible” means that one must obey the constraints that exist for these coefficients (Onsager symmetry, isotropy, sign) and arbitrarily setting some of them equal to zero is not always possible. A certain amount of trial and error is also required to avoid unduly reducing the model’s predictive power.

The usual vector notation is used for tensors of order 1, i.e., a bold faced letter is used when the subscript can be suppressed. Vector forces are here denoted X(with components Xi), and vector currents are denoted Y(with components Yi), with a superscript to discriminate between the different currents and forces in an alphabetic order aa shown below. (Yhas been chosen instead of the traditional Jso as to avoid confusion with the currents related to the energy equations.) The usual tensor conventions are assumed for latin subscripts (see “A note on notation” in the Introduction). Concerning the second order tensors Θjiand jVi, on the right-hand side of Eq. (77), one sets ΔjijViand, referring to Appendix 3, especially to Eq. (115), one writes

Δij=XAδij+XijB+εijkX˜k,E78
Θij=YAδij+YijB,E79

where XAand XijBare found by replacing Zby Δin Eqs. (115). Note that YijB, being symmetric, has no antisymmetric part.

Referring now to the vectors on the right-hand side of Eq. (77), one introduces the following notation:

XkC=kTXkD=FkFXkα=ΦkAα,YkC=Kk/TYkD=VkYkα=Ikα.E80

Using Eq. (116), one easily finds that

TQ=3XAYA+XijBYijB+XCYC+XDYD+α=1νXαYα.E81

Note that Qis a proper scalar since the source of entropy is independent of the coordinate system, and does not change sign when the coordinate system changes handedness. The right-hand side of Eq. (81) is then a sum of scalar products of proper tensors, the only pseudo vector, X˜k, having dropped out.

Referring to the first sentence of this section, one writes the currents as linear combinations of the forces:

YiC=LiCAXA+LiklCBXklB+LikCCXkC+LikCDXkD+β=1νLikXkβYiD=LiDAXA+LiklDBXklB+LikDCXkC+LikDDXkD+β=1νLikXkβYiα=LiαAXA+LiklαBXklB+LikαCXkC+LikαDXkD+β=1νLikαβXkβ,α=1ν.

In these expressions, the Lare the phenomenological coefficients. They are independent of the generalised forces, but they can depend on the temperature and the component densities. They are tensors, their orders being equal to the number of subscripts. They obey the Onsager relations [11, 15]: any coefficient with nsubscripts, denoted say, LnMN, obeys

LnMN=LnNM,MN=ABCD1ν,E83

where the nsubscripts are the same but not necessarily in the same order. Details concerning subscripts are not needed in what follows.

### 5.1 The viscosity tensor Θij

With hindsight, one knows that the upscaled viscosity tensor is not required, so that one is justified in setting equal to zero all the L‘s in the first two lines of the system of Eqs. (82), and also setting to zero all the L‘s related to the zeroed ones by the Onsager symmetry:

LAZ=LBZ=LZA=LZB=0,foranyZ.E84

Then YA=YijB=0, so that

Θij=0,E85

and the source term of the entropy equation reduces to

TQ=XCYC+XDYD+α=1νXαYα.E86

### 5.2 Vector currents and forces

System (82) now reduces to linear relations between vectors, the coefficients being second order tensors:

YiC=LikCCXkC+LikCDXkD+β=1νLikXkβE87
YiD=LikCDXkC+LikDDXkD+β=1νLikXkβE88
Yiα=LikXkC+LikXkD+β=1νLikαβXkβ,α=1ν,E89

where Onsager symmetry is accounted for. Remembering that Yiα=Ikα, Eq. (36) leads to three constraints on the L-coefficients of Eq. (89):

α=1νLik=α=1νLik=α=1νLikαβ=0.E90

The linear combinations above show the possibilities of constructing models where interactions between thermal conduction, fluid flow, and mass diffusion are quantified by choosing the Lcoefficients. Also, thermal conduction and permeability can be modelled by second order tensors, through tensors LikCCand LikDD.

However, if one limits oneself to the minimal model where: (i) the diffusive entropy current, Ki/T, is only due to the kTforce, (ii) the FkFfrictional force is only due to fluid velocity, Vi, (iii) each non-convective mass current, Iiα, is only due to the gradients of the Aβ, then one requires

LikCD=0,andLik=Lik=0forallα,E91

not violating the constraints in Eqs. (90).

Note that it is not possible to simplify the model to the extent that all cross-couplings are eliminated, since that would imply that Likαβare zero except when the superscripts are equal: such a matrix would not obey the third constraint in (90). One could of course set Likαβ=0for all αand β, but that would eliminate all the non-convective mass currents from the model, and probably make it useless.

In the minimal model one can add a fourth requirement to the three above: the upscaled fluid and the upscaled medium are isotropic. Then further simplifications result since the remaining second order tensors, LikCC, LikDD, and Likαβ, that are properties of the solid and the fluid, must be invariant under rotations of the coordinate axes. Such tensors are called isotropic and it can be shown that an isotropic second order tensor is proportional to the Kronecker delta (see Ref. [15]). Isotropy thus introduces the following restrictions:

LikCC=CCδik,LikDD=DDδik,Likαβ=αβδik,E92

and one gets

YC=CCXC,XD=1/DDYD,Yα=β=1ναβXβ.E93

The source term of the entropy equation is now

TQ=CCXC2+1/DDYD2+α=1νβ=1ναβXαXβ.E94

Keeping in mind that the αβ-matrix is symmetric, and that the sum of all elements on the same line or the same column is zero, it can be shown that

α=1νβ=1ναβXαXβ=α,β,α<βαβXαXβ2.E95

According to the expression for the source of entropy above, the remaining phenomenological coefficients must then satisfy

CC>0,DD>0,αβ<0,α<β.E96

These coefficients are determined below.

Using the notation of expressions (80) in the first of expressions (93) one gets Kiproportional to iT. Now Ki/Tis an entropy transport by diffusion, so that Kiis heat transport by diffusion. To recover Fourier’s law one sets

CC=k/T,implyingKi=kiT,E97

where k>0is the thermal conductivity of the averaged solid–fluid system. See Ref. [3].

The second of expressions (93) gives FiF=Vi/DD, and it is shown below that the choice

DD=KΦ2η,impliesVi=KΦηiPbFRfi,E98

away from interphase regions. In this expression, ηis the viscosity of the pore fluid, and Kis the absolute permeability. As already mentioned, the phenomenological coefficients can be functions of the component densities and of the temperature, and it is known that ηis such a function (see Section 3.4 in Ref. [3]). The possibility of letting Kbe such a function is not used in the minimal model. To prove the implication (98), one uses Eq. (42), with the upscaled stress tensor given by Eqs. (60), (61), and (85). Using Eq. (41) one obtains:

Vi=KΦηiPbFRfi+RtVi+VkkViα=1νΛαRαi2Rα.E99

The sum over αvanishes in the bulk fluid and one is left with the Darcy formula with an additional term, proportional to the material derivative of the velocity. One can carry out order of magnitude estimates of the three first terms in the square brackets above, in the manner of Section 4.1 of Ref. [3]. Using the numerical values given in the Appendix of the same reference, one easily finds that, if the gradient of pressure is of order 1, then the gravity term Rfiis of order 103, and the material derivative term is of order 109. Neglecting the material derivative, one obtains the Darcy formula, modified by terms that only become significant inside the interphase regions. Specialised to a one component fluid, this modified Darcy formula is:

Vi=KΦηiPbFRfi+ΛKRΦηi2R.E100

It is shown in Ref. [4] that the added non-Darcy term can, in some well-defined flow types, produce a relative permeability when its numerical contribution is taken away as an added term, then put back as a multiplicative factor to the Darcy term. However, it is concluded in Ref. [4] that relative permeabilities cannot capture the full complexity of two-phase flow.

Given below is another version of Vi, that follows by application of the Gibbs-Duhem equation, obtained by differentiating FbFand using the expression of dFbF(Eqs. (2)):

Vi=KΦηα=1νRαiAβ+SFiT,E101

where Aαis given by Eq. (70). Using now the notation of expressions (79) in the third of expressions (93) one gets the diffuse mass current of component α:

Ikα=Φβ=1ναβkAβ,E102

where it is reminded that the αβ-matrix is symmetric, its off diagonal elements are negative, and the elements of any line (or column) sum to zero. In addition, the elements can be functions of the densities. This last turns out to be extremely useful because the kCFβbehave as 1/Rβwhen Rβgoes to zero, which is not acceptable when the differential equations of the model are solved numerically. In the minimal model one then sets

αβ=aRαRβ2,αβ,E103

using the same positive number afor all elements. Each diagonal element ααis then the negative sum of the off diagonal elements of line α. For an example of such a matrix in the case ν=3, see Refs. [2, 3].

### 5.3 Closing details

It is easy to see that the first two restrictions displayed in Eqs. (92) do not implicate any other assumptions done in the minimal model so that the restrictions can be lifted, either singly or together, thus allowing thermal conductivity and/or permeability to be represented by a second order tensor when experiments indicate that such upgrading is required.

A non-thermal version of the model consists of νmass balance equations, see Eqs. (41), where the Darcy velocity is given by Eq. (101) and the mass diffusion velocities by Eqs. (102) and (103), the auxiliary variables being defined by Eqs. (70). Fot the thermal version of the model one needs to re-write the entropy equation, Eq. (75), in terms of temperature: see Section 2.3 of Ref. [3].

Start and boundary conditions must be supplied for the numerical solutions of the differential equations of the model. Special attention must be taken with the boundary conditions since the equations are of the fourth degree in the space variables. See Ref. [3] for a detailed presentation.

The transport coefficients, kand η, of Eqs. (97) and (98), are needed as functions of Tand of the Rα. See Ref. [3] and references given there.

The central thermodynamical function of the model is the Helmholtz function, especially in the bulk, introduced by Eq. (1). It is calculated from the equation of state of the mixture considered, which must be van der Waals or related (Redlich-Kwong, ...) so that, for temperatures less than the critical, regions of unstable fluid insure the existence of interphase regions; association terms must be included for the polar molecules of the mixture (see Ref. [16] and references given there).

It is shown in Ref. [16] that, independently of the equation of state that is chosen, the Helmholtz function contains the sum of the Helmholtz functions of the components, each considered as a gas where molecular interactions are neglected (ideal). Each ideal gas Helmholtz function contains a function of Tthat drops out under differentiation when Tis assumed constant, but is important for the thermal model in accounting for the energy stored in the internal degrees of freedom of the molecules. The Helmholtz function of one mole of an ideal gas of component αis (see Refs. [2, 16] and references given there):

fid.gas,α=RGTlnVαTevα,E104

where vαis the molar volume, RGis the gas constant, e=exp1, and VαTis a function of temperature with the dimension of a volume per mole. Vαdepends on the atomic masses, the principal moments of inertia of the molecule, .... It can be obtained from statistical physics, as shown in Ref. [16].

It is briefly shown below that Vαcan also be obtained from an expression of fid.gas,αin terms of the heat capacity of the component. Assuming Pvα=RGT, one easily finds that, with Tand vαas independent variables, the molar internal energy uα(the “id.gas” superscript is dropped for simplicity) is independent of vα, and that one can write

duα=cVαTdT,E105

where cVαis the molar heat capacity of component α. The molar entropy, sα, obeying Tdsα=duα+Pdvα, one obtains

dsα=cVαTdT/T+RGdvα/vα.E106

Integrating the last two displayed equations, one obtains fid.gas,αas uαTsα:

fid.gas,α=u0αTs0αRGTlnvαv0α+T0T1TTcVαTdT,E107

where u0αand s0αare the internal energy and the entropy at a reference state where volume, temperature, and pressure are v0α, T0, and P0=RGT0/v0α. Equating the right-hand sides of Eqs. (104) and (107), one obtains an expression for TlnVαin terms of the experimentally measurable function cVα. As stated in Ref. [16], it is in fact the derivative with respect to Tof TlnVαthat is needed in the differential equations of the thermal model. Assuming s0α=0, one obtains

ddTRGTlnVαT=RGlneRGT0P0T0TcVαTTdT.E108

Concerning examples of numerical solutions of the equations of the minimal model, see Ref. [2] for phase segregation, and for coning at uniform temperature; see Ref. [3] for an injection-production situation at variable temperature.

## A.1 Appendix 1: Assumptions

A0: None of the chemical species completely wets the rock

A1: There are no sources or sinks

A2: There is no loss of energy by radiation

A3: There are no chemical reactions between the chemical components

A4: The solid is perfectly rigid

A5: The external force per unit volume fiis approximately constant in the averaging volume

A6: At each point, the difference between the solid and fluid temperatures is negligible

A7: The Λαare constant numbers inside the porous medium

A8: The porosity Φis uniform in space

## A.2 Appendix 2: Derivatives of distibutions – Formulas

The following is a set of formulas for the space derivatives of distributions having discontinuities across a given surface Σ. For proofs, see [12, 13].

Consider a function fxtsymbolising a physical quantity, continuous in time but dicontinuous in space across a surface Σthat divides space in two regions, called A-side and B-side; f=fBon the B-side, f=fAon the A-side. Σhas a normal vector nat each point, pointing towards A:

fxt=fAxtxontheAsideofΣfBxtxontheBsideofΣ,E109

One now defines the following regular distributions:

kfxt=kfAxtxontheAsideofΣkfBxtxontheBsideofΣ,E110
tfxt=tfAxtxontheAsideofΣtfBxtxontheBsideofΣ.E111

Then [12, 13].

tf=tf,E112
kf=kf+fAΣfBΣnkδΣ,E113

where δΣis a surface-Dirac distribution, its action δΣφon a so-called test function φ[13] being defined as:

δΣφ=Σφ.E114

## A.3 Appendix 3: Formulas concerning tensors

Let Zijbe an arbitrary second order tensor. It can be shown (see Ref. [15]) that the following expression has general validity:

Zij=Zδij+Ẑij+εijkZ˜k,whereZ=13Zkk,Ẑij=12Zij+Zji13Zkkδij,Z˜k=12εklmZlm.E115

If Zijis a tensor, then Zis a scalar, Ẑis a symmetric traceless tensor, and Z˜kis a pseudo-vector. A most useful property is as follows: if Uijand Vijare two tensors then

UijVij=3UV+ÛijV̂ij+2U˜iV˜i.E116

## Notes

• The author first became aware of this paper in 2015. There is one important difference with what is presented in the pages that follow. See the footnote in Section 4.
• This is where the present paper and Ref. [8] differ most.

chapter PDF

## More

© 2021 The Author(s). Licensee IntechOpen. This chapter is distributed under the terms of the Creative Commons Attribution 3.0 License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

## How to cite and reference

### Cite this chapter Copy to clipboard

Paul Papatzacos (January 19th 2021). Porous Flow with Diffuse Interfaces [Online First], IntechOpen, DOI: 10.5772/intechopen.95474. Available from: