Analysis of a single-phase flow mixture model

# Existence analysis of a single-phase flow mixture model with van der Waals pressure

Ansgar Jüngel A.J.: Institute for Analysis and Scientific Computing, Vienna University of Technology, Wiedner Hauptstraße 8–10, 1040 Wien, Austria Jiří Mikyška J.M.: Department of Mathematics, Faculty of Nuclear Sciences and Physical Engineering, Czech Technical University in Prague, Trojanova 13, 12000 Prague 2, Czech Republic  and  Nicola Zamponi N.Z.: Institute for Analysis and Scientific Computing, Vienna University of Technology, Wiedner Hauptstraße 8–10, 1040 Wien, Austria
July 7, 2019
###### Abstract.

The transport of single-phase fluid mixtures in porous media is described by cross-diffusion equations for the mass densities. The equations are obtained in a thermodynamic consistent way from mass balance, Darcy’s law, and the van der Waals equation of state for mixtures. The model consists of parabolic equations with cross diffusion with a hypocoercive diffusion operator. The global-in-time existence of weak solutions in a bounded domain with equilibrium boundary conditions is proved, extending the boundedness-by-entropy method. Based on the free energy inequality, the large-time convergence of the solution to the constant equilibrium mass density is shown. For the two-species model and specific diffusion matrices, an integral inequality is proved, which reveals a minimum principle for the mass fractions. Without mass diffusion, the two-dimensional pressure is shown to converge exponentially fast to a constant. Numerical examples in one space dimension illustrate this convergence.

###### Key words and phrases:
Cross diffusion, single-phase flow, van der Waals pressure, existence of weak solutions, large-time asymptotics, maximum principle.
###### 2000 Mathematics Subject Classification:
35K51, 76S05
The authors have been partially supported by the bilateral Czech-Austrian project CZ 10/2015. The first and last authors acknowledge partial support from the Austrian Science Fund (FWF), grants P22108, P24304, and W1245. The second author acknowledges support from the Student Grant Agency of the Czech Technical University in Prague, project no. SGS14/206/OHK4/3T/14

## 1. Introduction

The transport of fluid mixtures in porous media has many important industrial applications like oil and gas extraction, dispersion of contaminants in underground water reservoirs, nuclear waste storage, and carbon sequestration. Although there are many papers on the modeling and numerical solution of such compositional models [1, 6, 7, 11, 17, 19], there are no results on their mathematical analysis. In this paper, we provide an existence analysis for a single-phase compositional model with van der Waals pressure in an isothermal setting. From a mathematical viewpoint, the model consists of strongly coupled degenerate parabolic equations for the mass densities. The cross-diffusion coupling and the hypocoercive diffusion operator constitute the main difficulty of the analysis.

Our analysis is a continuation of the program of the first and third author to develop a theory for cross-diffusion equations possessing an entropy (here: free energy) structure [13, 23]. The mathematical novelties are the complex structure of the equations and the observation that the solution of the binary model, for specific diffusion matrices, satisfies an unexpected integral inequality giving rise to a minimum principle, which generally does not hold for strongly coupled diffusion systems.

### Model equations

More specifically, we consider an isothermal fluid mixture of mass densities in a domain (), whose evolution is governed by the transport equations

 (1) ∂tci=div(ci∇p+εn∑j=1Dij(c)∇μj),x∈Ω, t>0, i=1,…,n,

where . The van der Waals pressure and the chemical potentials are given by

 (2) p =ctot1−∑nj=1bjcj−n∑i,j=1aijcicj, (3) μi =logci−log(1−n∑j=1bjcj)+bictot1−∑nj=1bjcj−2n∑j=1aijcj.

These expressions are well defined if a.e., where

 (4) D={(c1,…,cn)∈Rn:ci>0 for% i=1,…,n, n∑j=1bjcj<1}.

Here, is the total mass density and is a (small) parameter. The parameter measures the attraction between the th and th species, and is a measure of the size of the molecules. The diffusion matrix is assumed to be symmetric and positive semidefinite. Moreover, we suppose that the following bound holds:

 (5) D0|Πv|2≤n∑i,j=1Dij(c)vivj≤D1|Πv|2for all v∈Rn, c∈D,

for some , , where is the projection on the subspace of orthogonal to . A property like (5) is known in the literature as hypocoercivity, that is, coercivity on a subspace of the considered vector space. In our case, the matrix in (5) is coercive on the orthogonal complement of the subspace generated by . Bound (5) is justified in the derivation of model (1)-(3), as the diffusion fluxes must sum up to zero (see Section 2).

Equation (2) is the van der Waals equation of state for mixtures, taking into account the finite size of the molecules. Equations (2)-(3) are derived from the Helmholtz free energy of the mixture; see (16) below. For details of the modeling and the underlying assumptions, we refer to Section 2.

We impose the boundary and initial conditions

 (6) μi=0on ∂Ω, t>0,ci(⋅,0)=c0iin Ω, i=1,…,n.

Note that we choose equilibrium boundary conditions. A physically more realistic choice would be to assume that the reservoir boundary is impermeable, leading to no-flux boundary conditions. However, conditions (6) are needed to obtain Sobolev estimates, together with the energy inequality (8) below. Numerical examples for homogeneous Neumann boundary conditions for the pressure in case are presented in Section 7.

Up to our knowledge, there are no analytical results for system (1)-(3) and (6). In the literature, Euler and Navier-Stokes models were considered with van der Waals pressure. For instance, the existence of global classical solutions to the corresponding Euler equations with small initial data was shown in . The existence of traveling waves in one-dimensional Navier-Stokes with capillarity was studied in . Furthermore, in  the existence and stability of shock fronts in the vanishing viscosity limit for Navier-Stokes equations with van der Waals type equations of state was established.

### Main difficulties

A straightforward computation shows that the Gibbs-Duhem relation holds. Therefore, (1) can be written as

 (7) ∂tci=divn∑j=1((cicj+εDij(c))∇μj),i=1,…,n,

which is a cross-diffusion system in the so-called entropy variables . The matrix is of rank one with two eigenvalues, a positive one and the other one equal to zero (with algebraic multiplicity ). Thus, if , system (1) is not parabolic in the sense of Petrovski , and an existence theory for such diffusion systems is highly nontrivial, which is the first difficulty. The property on the eigenvalues is reflected in the energy estimate. Indeed, a formal computation, made rigorous below, shows that

 (8) ddt∫ΩF(c)dx+∫Ω|∇p|2dx+ε∫Ω∇μ:D(c)∇μ dx≤0.

In case we obtain only one gradient estimate for which is not sufficient for the analysis. There exist some results for so-called strongly degenerate parabolic equations (for which the diffusion matrix vanishes in some subset of positive -dimensional measure) . However, the techniques cannot be applied to the present problem. Therefore, we need to assume that . Then the gradient estimates for and together with the boundary conditions (6) yield uniform bounds, which are the basis of the existence proof. The behavior of the solutions for are studied numerically in Section 7.

The second difficulty is the invertibility of the relation between and , i.e. to define for given the mass density vector , where and is defined by (3). A key ingredient for the proof is the positive definiteness of the Hessian of the free energy since . This is only possible under a smallness condition on the eigenvalues of ; see Lemma 6. This condition is not surprising since it just means that phase separation is prohibited. The analysis of multiphase flows requires completely different mathematical techniques; see, e.g.,  for phase transitions in Euler equations with van der Waals pressure.

The third difficulty is the proof of a.e. This property is needed to define and through (2)-(3), but generally a maximum principle cannot be applied to the strongly coupled system (1). The idea is to employ the boundedness-by-entropy method as in [13, 23], i.e. to work with the entropy variables . We show first the existence of weak solutions to a regularized version of (7), define and perform the de-regularization limit to obtain the existence of a weak solution to (1). Since a.e. by definition of , turns out to be bounded. This idea avoids the maximum principle and is the core of the boundedness-by-entropy method. Let us now detail our main results.

### Global existence of solutions

Using the boundedness-by-entropy method and the energy inequality (8), we are able to prove the global existence of bounded weak solutions. We set and .

###### Theorem 1 (Existence and large-time asymptotics).

Let , , be Lebesgue measurable and let such that , where , , is defined by (3). Furthermore, let the matrices and be symmetric and satisfy (5) as well as

 (9) κ:=116mini=1,…,nbimaxi=1,…,nbi−λ∗mini=1,…,nbi>0,K:=1−max1≤i,j≤nb−1iaij>0,

respectively, where is the maximal eigenvalue of . Then:

1. There exists a weak solution to (1)-(6) satisfying the free energy inequality (8) and

 ci−cΓi∈L2(0,∞;H10(Ω))∩H1(0,∞;H−1(Ω)),i=1,…,n, |∇p|∈L2(0,∞;L2(Ω)),logctot∈L∞(0,∞;L2(Ω)).
2. There exists a constant , depending on and such that

 n∑i=1∥ci(t)−cΓi∥2L2(Ω)≤C1+tfor t>0.

The idea of the large-time asymptotics of is to exploit the energy inequality (8). Since it is difficult to relate the free energy and its energy dissipation , we cannot prove an exponential decay rate although numerical experiments in  and Section 7 indicate that this is the case even when . Instead, we show for the relative energy that, for some constant and some nonnegative function ,

from which we deduce that the convergence is of order as . Since the free energy is strictly convex, by Lemma 6 below, we obtain convergence in the norm.

### An integral inequality

If , we obtain only a gradient estimate for . This lack of parabolicity is compensated by the following – surprising – integral identity,

 (10)

for arbitrary functions ; see the Appendix for a formal proof. This means that there exists a family of conserved quantities depending on a function of variables. It is unclear whether this identity is sufficient to perform the limit and to prove the existence of a solution to (1) with .

If , the integral identity (10) does not hold in general. However, for specific diffusion matrices , the following inequality holds in place of (10):

 (11) ∫Ωctot(t)f(c1(t)ctot(t),…,cn−1(t)ctot(t))dx≤∫Ωc0totf(c01c0tot,…,c0n−1c0tot)dx,t>0,

for functions specified in Theorem 3 below. Interestingly, this implies a minimum principle for . A choice of the diffusion matrix ensuring the validity of (11) is, for given , with , in ,

 (12) D(c)=α(c)(F′′)−1+β(c)c⊗c,

where is the Hessian of the free energy . Clearly, is bounded and positive definite (although not strictly) for . In particular, the constraint does not hold, and so the assumptions of Theorem 1 are not satisfied. However, with this choice of , equation (1) becomes

 (13) ∂tci=div((1+εβ(c))ci∇p+εα(c)∇ci) i=1,…,n,

and the existence proof for (13) is simpler than in the case where satisfies (5).

###### Corollary 2 (to Theorem 1).

Let , , be Lebesgue measurable and let , where , , is defined by (3). Furthermore, let the matrices and be symmetric and satisfy (9), (12). Then there exists a weak solution to (1)-(3), (6), satisfying the free energy inequality (8) and, for ,

 ci−cΓi ∈L2(0,∞;H10(Ω))∩H1(0,∞;H−1(Ω)), ∇√ci, ∇p, ∇logci ∈L2(0,∞;L2(Ω)), logci ∈L∞(0,∞;L1(Ω)).

Our second main result reads as follows.

###### Theorem 3 (Integral inequality and minimum principle).

Let for on . Under the assumptions of Corollary 2, the solution to (1)-(3), (6) constructed in Corollary 2 satisfies (11) for all functions such that its Hessian is positive semidefinite in and

 (14) f(cΓ1cΓtot,…,cΓn−1cΓtot)=0,∣∣∣f′(cΓ1cΓtot,…,cΓn−1cΓtot)∣∣∣=0.

Moreover, for any ,

 infΩ×(0,∞)cictot≥min{infΩc0ic0tot, cΓicΓtot}.

### Exponential convergence of the pressure

In the degenerate situation , we are able to show an exponential decay rate for the pressure , at least for sufficiently smooth solutions whose existence is assumed. The key idea of the proof is to analyze the parabolic equation satisfied by ,

 ∂tp=˜DΔp+|∇p|2,where ˜D=n∑i,j=1cicj∂2F∂ci∂cj.

Because of the quadratic gradient term, we need a smallness assumption on at time . Thus, the exponential convergence result holds sufficiently close to equilibrium.

###### Theorem 4 (Exponential decay of the pressure).

Let , , and let be a solution to (1)-(2) with isobaric boundary conditions on , , for some constant . Let . We assume that

 ∇ci∈L4loc(0,∞);L2(Ω)),∇p∈C0([0,∞);L2(Ω))∩L2(0,T;H1(Ω)),

and for any . Then there exists , which depends on and , such that if , then, for some ,

 ∥∇p(c(t))∥L2(Ω)≤∥∇p(c0)∥L2(Ω)e−λt,t>0.

The paper is organized as follows. Details on the modeling of the fluid mixture are presented in Section 2. Auxiliary results on the Hessian of the free energy, the relation between and , and the diffusion matrix (12) are shown in Section 3. In Section 4, we prove Theorem 1 and Corollary 2, while the proofs of Theorems 3 and 4 are presented in Section 5 and 6, respectively. The evolution of the one-dimensional mass densities and the pressure are illustrated numerically in Section 7 for the case . Finally, identity (10) is verified in the Appendix.

## 2. Modeling and energy equation

We consider the isothermal flow of chemical components in a porous domain with porosity . The transport of the partial mass densities is governed by the balance equations for the mass,

 ∂t(φci)+div(civi)=0,i=1,…,n,

where is the partial velocity of the th species. In order to derive equations for the mass densities only, we impose some simplifying assumptions. To shorten the presentation, we set all physical constants equal to one. Moreover, we set to simplify the mathematical analysis. Our results will be also valid for (smooth) space-dependent porosities. Introducing the diffusion fluxes by , where is the barycentric velocity and denotes the total mass density, the balance equations become

 (15) ∂tci+div(civ+Ji)=0,i=1,…,n.

We suppose that the barycentric velocity is given by Darcy’s law , where is the fluid pressure. We refer to  for a justification of this law. The second assumption is that the diffusion fluxes are driven by the gradients of the chemical potentials , i.e.  for ; see, e.g., [14, Section 4.3]. Here, is some number and are diffusion coefficients depending on . According to Onsager’s principle of thermodynamics, the diffusion matrix has to be symmetric and positive semidefinite; moreover, for consistency with the definition , it must hold that for .

The equations are closed by specifying the Helmholtz free energy density

 (16) F(c)=n∑i=1ci(logci−1)−ctotlog(1−n∑j=1bjcj)−n∑i,j=1aijcicj,

where and are positive numbers, and is symmetric. The first term in the free energy is the internal energy and the remaining two terms are the energy contributions of the van der Waals gas [12, Formula (4.3)].

The third assumption is that the fluid is in a single state, i.e., no phase-splitting occurs. Mathematically, this means that the free energy must be convex. This is the case if the maximal eigenvalue of is sufficiently small; see Lemma 6. The single-state assumption is restrictive from a physical viewpoint. It may be overcome by considering the transport equations for each phase separately and imposing suitable boundary conditions at the interface [14, Section 1]. However, this leads to free-boundary cross-diffusion problems which we are not able to treat mathematically. Another approach would be to consider a two-phase (or even multi-phase) compositional model with overlapping of different phases, like in . In such a situation, a new formulation of the thermodynamic equilibrium based upon the minimization of the Helmholtz free energy is employed to describe the splitting of components among different phases.

The chemical potentials are defined in terms of the free energy by

 μi=∂F∂ci=logci−log(1−n∑j=1bjcj)+bictot1−∑nj=1bjcj−2n∑j=1aijcj,

and the pressure is determined by the Gibbs-Duhem equation [4, Formula (64)]

 (17) p=n∑i=1ciμi−F(c)=ctot1−∑nj=1bjcj−n∑i,j=1aijcicj.

This describes the van der Waals equation of state for mixtures, where the parameter is a measure of the attractive force between the molecules of the th and th species, and the parameter is a measure of the size of the molecules. The pressure stays finite if , which means that the mass densities are bounded. In the literature, many modifications of the attractive term have been proposed. Examples are the so-called Peng-Robinson and Soave-Redlich-Kwong equations; see .

Taking the gradient of (17) and observing that , (17) can be written as

 (18) ∇p=n∑i=1ci∇μi.

Therefore, we can formulate (15) as the cross-diffusion equations

 ∂tci=div(n∑j=1(cicj+εDij)∇μj),i=1,…,n.

Multiplying this equation by , summing over , observing again that , and integrating by parts, we arrive at the energy equation

 ddt∫ΩF(c)dx=∫Ωn∑i=1μi∂tcidx=−∫Ω∣∣∣n∑i=1ci∇μi∣∣∣2dx−ε∫Ωn∑i,j=1Dij∇μi⋅∇μjdx.

Since is assumed to be positive definite on , where , and , this gives, thanks to Lemma 5, estimates for and, thanks to the equilibrium boundary condition and Poincaré’s inequality, estimates for .

## 3. Auxiliary results

First we show a result estimating the norms of two vectors from below.

###### Lemma 5.

Let , be such that . Then, for any ,

 |α⋅v|2+|v−(β⋅v)β|2≥14(α⋅β)2|v|2.

The constant is not optimal. For instance, if , we have the theorem of Pythagoras, .

###### Proof.

Let be the projection of on and be the orthogonal part. Then, clearly, . By Young’s inequality with and , we have

 |α⋅v|2+|v−(β⋅v)β|2 =|α⋅(w+w⊥)|2+|w⊥|2 =(α⋅w)2+(α⋅w⊥)2+2(α⋅w)(α⋅w⊥)+|w⊥|2 ≥(1−δ)(α⋅w)2+(1−δ−1)(α⋅w⊥)2+|w⊥|2 =14(α⋅w)2−13(α⋅w⊥)2+|w⊥|2 ≥14(β⋅v)2(α⋅β)2−13|w⊥|2+|w⊥|2.

We deduce from that , and thus,

 |α⋅v|2+|v−(β⋅v)β|2 ≥14(α⋅β)2|w|2+23(α⋅β)2|w⊥|2 ≥14(α⋅β)2(|w|2+|w⊥|2)=14(α⋅β)2|v|2,

finishing the proof. ∎

###### Lemma 6 (Positive definiteness of F′′).

Let , defined in the pressure relation (2), be a symmetric matrix whose maximal eigenvalue satisfies (9). Then the Hessian of the free energy is positive definite, i.e.

 v⋅F′′(c)v≥κn∑i=1v2icifor all c∈D, v∈Rn,

where is given by (9). In particular, is uniformly positive definite.

###### Proof.

A straightforward computation shows that , where

 Bij=(bi+bj)σ+bibjctotσ2+δijci,σ=11−∑nj=1bjcj≥1.

Let . It holds that

 n∑i,j=1Bijvivj =(σ√ctotn∑j=1bjvj+1√ctotn∑i=1vi)2+n∑i=1v2ici−1ctot(n∑i=1vi)2.

Defining , , and for , the quadratic form can be rewritten as

 n∑i,j=1Bijvivj=(ˆα⋅w)2+|w−(β⋅w)β|2.

Since , we may define , which yields

 (19) n∑i,j=1Bijvivj≥(α⋅w)2+|w−(β⋅w)β|2.

The norm of can be estimated from above:

 |ˆα|2≤σ2ctotmaxj=1,…,nbjn∑i=1bici+2σn∑i=1bici+1.

Since , we have or , and . Therefore,

 |ˆα|2≤σ2maxj=1,…,nbjminj=1,…,nbj+2σ2+σ2=(3+maxj=1,…,nbjminj=1,…,nbj)σ2.

We infer that is strictly positive:

 (α⋅β)2=|ˆα|−2(1+σn∑i=1bici)2≥(σ−1+∑ni=1bici)23+maxj=1,…,nbjminj=1,…,nbj=13+maxj=1,…,nbjminj=1,…,nbj.

We apply Lemma 5 to (19) to obtain

 4n∑i,j=1Bijvivj≥(α⋅β)2|w|2≥|w|23+maxj=1,…,nbjminj=1,…,nbj,

which, since , implies that

 (20) 4n∑i,j=1Bijvivj≥minj=1,…,nbj3minj=1,…,nbj+maxj=1,…,nbjn∑i=1v2ici≥minj=1,…,nbj4maxj=1,…,nbjn∑i=1v2ici.

The relation and the definition of allow us to write

 n∑i,j=1aijvivj≤ctotn∑i,j=1aijvi√civj√cj≤λ∗minj=1,…,nbjn∑i=1v2ici.

This, together with (20), yields the desired lower bound for . ∎

###### Lemma 7 (Invertibility of c↦μ).

The mapping , is invertible.

###### Proof.

Since is positive definite in , it follows that is one-to-one and the image is open. We claim that is also closed. Then , and the proof is complete.

Let , , define a sequence in such that as . The claim follows if we prove that there exists such that . Since varies in a bounded subset of , the theorem of Bolzano-Weierstraß implies the existence of a subsequence, which is not relabeled, such that converges to some as , where , . We assume, by contradiction, that . Let us distinguish two cases.

Case 1: There exists such that . If , then (3) implies that , which contradicts the fact that is convergent. Thus it holds that . This means that for some . However, choosing in (3) and exploiting the relation leads to , contradiction.

Case 2: For all , it holds that and . Arguing as in case 1, it follows that for all , which is absurd.

We conclude that , which finishes the proof. ∎

###### Lemma 8.

Let , for , where and satisfies (5). Then for all and ,

 v⋅Bv≥kB(c2tot|v|2+|Πv|2),v⋅ˆBv≥k′B|v|2,

where , only depend on (the constant in (5)) and .

###### Proof.

From (5), , and it follows that

 v⋅Bv ≥(c⋅v)2+εD0|Πv|2 ≥c2tot(1n(ˆc⋅v)2+ε2D0min{b1,…,bn}2|(I−ℓ⊗ℓ)v|2)+ε2D0|(I−ℓ⊗ℓ)v|2,

Applying Lemma 5 with and to the expression in the brackets yields

 v⋅Bv≥14min{1n,ε2D0min{b1,…,bn}2}c2tot(ˆc⋅ℓ)2|v|2+12εD0|(I−ℓ⊗ℓ)v|2.