A mixed finite element approximation for Darcy-Forchheimer flows of slightly compressible fluids

# A mixed finite element approximation for Darcy-Forchheimer flows of slightly compressible fluids

Thinh Kieu 222Department of Mathematics, University of North Georgia, Gainesville Campus, 3820 Mundy Mill Rd., Oakwood, GA 30566, U.S.A. (thinh.kieu@ung.edu).
today
###### Abstract

In this paper, we consider the generalized Forchheimer flows for slightly compressible fluids in porous media. Using Muskat’s and Ward’s general form of Forchheimer equations, we describe the flow of a single-phase fluid in by a nonlinear degenerate system of density and momentum. A mixed finite element method is proposed for the approximation of the solution of the above system. The stability of the approximations are proved; the error estimates are derived for the numerical approximations for both continuous and discrete time procedures. The continuous dependence of numerical solutions on physical parameters are demonstrated. Experimental studies are presented regarding convergence rates and showing the dependence of the solution on the physical parameters.

Key words. Porous media, error analysis, slightly compressible fluid, dependence on parameters, numerical analysis.

AMS subject classifications. 65M12, 65M15, 65M60, 35Q35, 76S05.

## 1 Introduction

The fluid flow through porous materials, e.g. soil, sand, aquifers, oil reservoir, plants, wood, bones, etc., is a great interest in the research community such as engineering, oil recovery, environmental and groundwater hydrology and medicine. Darcy law, which is the linear relation between the velocity vector and the pressure gradient, is used to describe fluid flow under low velocity and low porosity conditions, see in [3]. It has been observed from many experiments that when the fluid’s velocity is high and porosity is nonuniform, the Darcy’s law becomes inadequate. Consequently, the attention has been attracted to the nonlinear equations for describing of this kind of flow. Dupuit and Forchheimer proposed a modified equation, known as Darcy- Forchheimer equation or generalized Forchheimer equation, by adding the nonlinear terms of velocity to Darcy law, see [16]. Since then, there have been a growing number of articles studying Darcy- Forchheimer equation in theoretical studies (e.g.[40, 41, 8, 7, 37, 32]) and numerical studies (e.g.[4, 10, 36, 18, 30, 38]).

It is well known that the mixed finite element method is among the popular numerical methods for the modeling flow in porous media because it produces the accurate results for both scalar (density or pressure) and vector (velocity or momentum) functions, see [36]. An analysis of mixed finite element method to a Darcy-Forchheimer steady state model was well studied in [35, 38]. The mixed methods for a nondegenerate system modeling flows in porous media was studied in [10, 36, 18, 30]. The authors in [1, 43, 14, 15] analyzed the mixed finite element approximations of the nonlinear degenerate system modeling water-gas flow in porous media. In their analysis, the Kirchhoff transformation is used to move the nonlinearity from coefficients to the gradient.

The objective of this paper is to analyze mixed finite element approximations to the solutions of the system of equations modeling the flows of a single-phase compressible fluid in porous media subject to the generalized Forchheimer law. This is a nonlinear degenerate system with coefficients depending on the density gradient and degenerating to zero as it approaches to infinity. The Kirchhoff transformation is not applicable for this system. For our degenerate equations, we combine the techniques developed in our previous works in [19, 20, 21, 22, 23, 24, 29, 28] and utilize the special structures of the equations to obtain the stability of the approximated solution and the continuous dependence of the solution on parameters. The error estimates also derives for the density and momentum.

The paper is structured as follows. In section §LABEL:Intrsec, we introduce the notations and the relevant results. In section §LABEL:MFEM, we defined a numerical approximation using mixed finite element approximations and the implicit backward difference time discretization to the initial boundary value problem (IVBP) LABEL:rho:eq. In section §LABEL:Bsec, we establish many estimates of the energy type norms for the approximate solution to the IVBP problem (LABEL:semidiscreteform) in Lebesgue norms in terms of the boundary data and the initial data. In section §LABEL:dependence-sec, we focus on proving the continuous dependence of the solution on the coefficients of Forchheimer polynomial . In order to obtain this, we first establish the perturbed monotonicity for our degenerate partial differential equations, see in (LABEL:quasimonotone). It is then proved in Theorem LABEL:DepCoeff that the difference between the two solutions, which corresponds to two different coefficient vectors and is estimated in terms of , see in (LABEL:ssc2). In section §LABEL:err-sec, we study in Theorem LABEL:err-thm the convergence and in Theorem LABEL:par-thm the dependence on coefficients of Forchheimer polynomial of the approximated solution to the problem (LABEL:semidiscreteform). Furthermore, we can specify the convergent rate. In section §LABEL:err-ful-sec, we study the fully discrete version of problem (LABEL:semidiscreteform). In Lemma LABEL:stab-appr, the stability of the approximated solution is proved. Theorems LABEL:Err-ful and LABEL:Dep-ful are for studying the error estimates and the continuous dependence on parameters of the numerical solution. In section §LABEL:Num-result, the numerical experiments in the two-dimensions using the standard finite elements are presented regarding the convergence rates and the dependence of the solution on the physical parameters.

## 2 Preliminaries and auxiliaries

We consider a fluid in porous medium occupying a bounded domain with boundary . Let , and be the spatial and time variables respectively. The fluid flow has velocity , pressure and density .

The Darcy–Forchheimer equation is studied in [2, 19, 20] of the form

 −∇p=N∑i=0ai|v|αiv. \hb@xt@.01(2.1)

where , are real (not necessarily integral) numbers, and the coefficients satisfy and .

In order to take into account the presence of density in the generalized Forchheimer equation, we modify (LABEL:gF) using the dimensional analysis by Muskat [34] and Ward [42]. They proposed the following equation for both laminar and turbulent flows in porous media:

 −∇p=F(vακα−32ρα−1μ2−α), where F is a function of one variable. \hb@xt@.01(2.2)

In particular, when , Ward [42] established from experimental data that

 −∇p=μκv+cFρ√κ|v|v,% where cF>0. \hb@xt@.01(2.3)

Combining (LABEL:gF) with the suggestive form (LABEL:W) for the dependence on and , we propose the following equation

 −∇p=N∑i=0aiραi|v|αiv. \hb@xt@.01(2.4)

Here, the viscosity and permeability are considered constant, and we do not specify the dependence of ’s on them.

Multiplying equation (LABEL:FM) to , we obtain

 g(|ρv|)ρv=−ρ∇p, \hb@xt@.01(2.5)

where the function is a generalized polynomial with non-negative coefficients. More precisely, the function is of the form

 g(s)=a0sα0+a1sα1+⋯+aNsαN,s≥0, \hb@xt@.01(2.6)

where , are real (not necessarily integral) numbers. The coefficients satisfy and . The number is the degree of and is denoted by . Denote the vectors of powers and coefficients by and . The class of functions as in (LABEL:eq2) is denoted by FP(), which is the abbreviation of “Forchheimer polynomials.” When the function in (LABEL:eq2) belongs to FP(), it is referred to as the Forchheimer polynomial.

For slightly compressible fluids, the state equation is

 dρdp=ρκ,κ=const.,κ≫1 \hb@xt@.01(2.7)

which yields

 ρ∇p=κ∇ρ. \hb@xt@.01(2.8)

It follows form (LABEL:eq2) and (LABEL:eq3) that

 g(|ρv|)ρv=−κ∇ρ. \hb@xt@.01(2.9)

Solving for from (LABEL:Gs) gives

 ρv=−κK(|κ∇ρ|)∇ρ, \hb@xt@.01(2.10)

where the function is defined for by

 K(ξ)=1g(s(ξ)), with s=s(ξ) being the unique non-% negative solution of sg(s)=ξ. \hb@xt@.01(2.11)

The continuity equation is

 ϕρt+div(ρv)=f, \hb@xt@.01(2.12)

where is the porosity, and is the external mass flow rate .

By rescaling the coefficients of the conductivity function , we can assume . We introduce the momentum variables

Combining (LABEL:ru) and (LABEL:con-law), we obtain equations in a density-momentum formulation

 {m+K(|∇ρ|)∇ρ=0,ϕρt+∇⋅m=f. \hb@xt@.01(2.13)

We will study the initial boundary value problem (IVBP) associated with the coupled system (LABEL:utporo). We will derive estimates for the solution, and establish their continuous dependence on the coefficients of the function in (LABEL:eq2). As seen below, the generalized permeability tensor is degenerate to zero when . This was not considered in existing papers, e.g. in [13, 30, 35]. Therefore, it creates an additional challenge and requires extra care in the proof and analysis.

Let in FP(). The following numbers are frequently used in our calculations:

 χ(→a)=max{a0,a1,…,aN,1a0,1aN}∈[1,∞), \hb@xt@.01(2.14) a=αNαN+1=deg gdeg g+1∈(0,1),β=2−a∈(1,2),λ=2−a1−a=ββ−1∈(2,∞). \hb@xt@.01(2.15)
###### Lemma 2.1 (cf. [2, 19], Lemma 2.1)

Let be in class FP(). For any , One has

1. and it decreases in .

2. For any , the function increases and .

3. Type of degeneracy

 c0(1+ξ)−a≤K(ξ,→a)≤c1(1+ξ)−a. \hb@xt@.01(2.16)
4. For all

 c2(ξm−a−δm−a)≤K(ξ,→a)ξm≤c3ξm−a, \hb@xt@.01(2.17)

where depend on and only.

In particular, when , , one has

 c2(ξ2−a−1)≤K(ξ,→a)ξ2≤c3ξ2−a. \hb@xt@.01(2.18)
5. Relation with its derivative

 −aK(ξ)≤K′(ξ)ξ≤0. \hb@xt@.01(2.19)

We define

 H(ξ,→a)=∫ξ20K(√s,→a)dsfor ξ≥0. \hb@xt@.01(2.20)

When vector is fixed, we denote and by and , respectively. The function can be compared with and by

 K(ξ)ξ2≤H(ξ)≤2K(ξ)ξ2,c4(ξ2−a−1)≤H(ξ)≤c5ξ2−a, \hb@xt@.01(2.21)

where depend on .

For convenience, we use the following notations: let and be two arbitrary vectors of the same length, including possible length . We denote by and their maximum and minimum vectors, respectively, with components and .

###### Lemma 2.2

Let be the set of admissible . For any coefficient vectors , , and any , one has

(i)

 |K(|y|,→a)y−K(|y′|,→a′)y′|≤(1+a)|y−y′|∫10K(|γ(t)|,→b(t))dt+d0(|y|∨|y′|)|→a−→a′|∫10K(|γ(t)|,→b(t))dt. \hb@xt@.01(2.22)

(ii)

 (K(|y|,→a)y−K(|y′|,→a′)y′)⋅(y−y′) ≥(1−a)|y−y′|2∫10K(|γ(t)|,→b(t))dt \hb@xt@.01(2.23) −d0(|y|∨|y′|)|→a−→a′||y−y′|∫10K(|γ(t)|,→b(t))dt,

 γ(t)=ty+(1−t)y′,→b(t)=→a+(1−t)→a′,d0=N(min{a0,a′0,(1+αN)aN,(1+αN)a′N})−1.

In particular, if then (LABEL:Umono) and (LABEL:quasimonotone) become

(iii)

 ∣∣K(|y′|)y′−K(|y|)y∣∣≤(1+a)|y′−y|∫10K(|γ(t)|,→a)dt. \hb@xt@.01(2.24)

(iv)

 (K(|y′|)y′−K(|y|)y)⋅(y′−y)≥(1−a)|y′−y|2∫10K(|γ(t)|,→a)dt. \hb@xt@.01(2.25)

Proof. (i). Let , and .

Case 1: The origin does not belong to the segment connecting and . Define

 z(t)=K(|γ(t)|,→b(t))γ(t)⋅→k.

By the Mean Value Theorem, we have

 Idef=\joinrel=[K(|y|,→a)y−K(|y′|,→a′)y′]⋅→k=z(1)−z(0)=∫10z′(t)dt.

Elementary calculations give

 I=∫10f1(t)dt+∫10f2(t)dtdef=\joinrel=I1+I2, \hb@xt@.01(2.26)

where

 f1(t) =K(|γ(t)|,→b(t))(y−y′)⋅→k+Kξ(|γ(t)|,→b(t))γ(t)⋅(y−y′)|γ(t)|γ(t)⋅→k, f2(t) =K→a(|γ(t)|,→b(t))(→a−→a′)γ(t)⋅→k.

Estimation of . Since

 |f1(t)|≤∣∣∣K(|γ(t)|,→b(t))(y−y′)+Kξ(|γ(t)|,→b(t))γ(t)⋅(y−y′)|γ(t)|γ(t)∣∣∣|→k|≤(1+a)K(|γ(t)|,→b(t))|y−y′||→k|,

we find that

 I1≤∫10|f1(t)|dt≤(1+a)|y−y′||→k|∫10K(|γ(t)|,→b(t))dt. \hb@xt@.01(2.27)

Estimation of . We find the partial derivative of in . For , taking the partial derivative in of the identity , we find that

 Kai(ξ,→a)=−gai+gs⋅saig2=−K(ξ,→a)gai+gs⋅saig.

From , we have for , which implies

Hence, we obtain

 Kai(ξ,→a)=−K(ξ,→a)gai+gs⋅−s⋅gaig+s⋅gsg=−K(ξ,→a)gaig+s⋅gs=−K(ξ,→a)sαig+s⋅gs. \hb@xt@.01(2.28)

This shows that

 N∑i=0|Kai(ξ,→a)|≤K(ξ,→a)1+sα1+⋯+sαNa0+(1+α1)a1sα1+⋯+(1+αN)aNsαN.

Using inequality , we have Hence,

 N∑i=0|Kai(ξ,→a)|≤K(ξ,→a)N(1+sαN)a0+(1+αN)aNsαN≤d(→a)K(ξ,→a),

where Thus,

 |K→a(ξ,→a)|≤d(→a)K(ξ,→a). \hb@xt@.01(2.29)

Using the estimate (LABEL:Ka), we bound

 |f2(t)|≤|K→a(|γ(t)|,→b(t))| |→a−→a′| |γ(t)| |→k|≤d(→b(t))K(|γ(t)|,→b(t)) |γ(t)| |→a−→a′| |→k|. \hb@xt@.01(2.30)

Since is positive for , the number , for all , can be bounded by Using the fact , (LABEL:ih2) yields

 |f2(t)|≤d0K(|γ(t)|,→b(t))(|y|∨|y′|)|→a−→a′| |→k|, \hb@xt@.01(2.31)

and consequently,

 I2≤∫10|f2(t)|dt≤d0(|y|∨|y′|)|→a−→a′||→k|∫10K(|γ(t)|,→b(t))|dt. \hb@xt@.01(2.32)

Thus, we obtain (LABEL:Umono) by combining (LABEL:Ifm), (LABEL:UI1) and (LABEL:UI2).

Case 2: The origin belongs to the segment connect . We replace by so that and as . Apply the above inequality for and , then let .

(ii) If then

 Kξ(|γ(t)|,→b(t))≥−aK(|γ(t)|,→b(t))|γ(t)|.

This implies

 f1(t)≥(1−a)K(|γ(t)|,→b(t))|y−y′|2. \hb@xt@.01(2.33)

It follows (LABEL:ks) that

 ∫10f1(t)dt≥(1−a)|y′−y|2∫10K(|γ(t)|,→b(t))dt, \hb@xt@.01(2.34)

and from (LABEL:Uh2), we see that

 I2≥−∫10|f2(t)|dt≥−d0(|y|∨|y′|)|→a−→a′||y′−y|∫10K(|γ(t)|,→b(t))|dt. \hb@xt@.01(2.35)

Thus, we obtain (LABEL:quasimonotone) by combining (LABEL:Ifm), (LABEL:int-h1) and (LABEL:I2geq).

Now we derive the trace estimates suitable for our nonlinear problem.

###### Lemma 2.3

Assume is a function defined on .

(i) If then there is a positive constant depending on such that for all ,

 ∫Γ|v|dσ≤C1∥v∥+ε∥∇v∥β0,β+C1ε−1β−1. \hb@xt@.01(2.36)

(ii) If and then there exists such that for all

 |⟨u,v⟩|≤ε(∥v∥2+∥∇v∥β0,β)+C2(ε−1∥u∥2L∞(Γ)+ε−1β−1∥u∥λL∞(Γ)). \hb@xt@.01(2.37)

Consequently,

 |⟨u,v⟩|≤14(∥v∥2+∥∇v∥β0,β)+C3(1+∥u∥λL∞(Γ)) \hb@xt@.01(2.38)

for a constant .

Proof. We recall the trace theorem

 ∫Γ|v|dx≤C∗∫Ω|v|dx+C∗∫Ω|∇v|dx,

for all where is a positive constant depending on . Using Young’s inequality, we obtain (LABEL:sec2).

(ii) We have

 |⟨u,v⟩|≤∥u∥L∞(Γ)∫Γ|v|dσ. \hb@xt@.01(2.39)

Using (LABEL:sec2) and Young’s inequality give

 |⟨u,v⟩|≤∥u∥L∞(Γ)(δ∥v∥2+C2∗|Ω|4δ−1+δ∥∇v∥β0,β+C1δ−1β−1).

If then (LABEL:bnd-eps) clearly holds true.

Otherwise, selecting and the fact that , we obtain (LABEL:bnd-eps).

Estimate (LABEL:bnd-est) follows by choosing in (LABEL:bnd-eps) and using Young’s inequality.

We recall a discrete version of Gronwall Lemma in backward difference form, which is useful later. It can be proven without much difficulty by following the ideas of the proof in Gronwall Lemma.

###### Lemma 2.4

Assume and the nonnegative sequences , satisfying

 an−an−1Δt−ℓan≤gn,n=1,2,3…

then

 an≤(1−ℓΔt)−n(a0+Δtn∑i=1(1−ℓΔt)i−1gi). \hb@xt@.01(2.40)

Proof. Let . Simple calculation shows that

 ¯an−¯an−1Δt=(1−ℓΔt)n−1(an−an−1Δt−ℓan)≤(1−ℓΔt)n−1gn.

 ¯an−¯a0Δt≤n∑i=1(1−ℓΔt)i−1gi,

and hence (LABEL:Gineq) holds true.

Notations: Let be the set of square integrable functions on and the space of -dimensional vectors with all the components in . We denote the inner product in either or . The notation means scalar norm or vector norm and represents the standard Lebesgue norm. Notation means the mixed Lebesgue norm.

For and any nonnegative integer, let denote a Sobolev space endowed with the norm Define with the norm .

Our estimates make use of coefficient-weighted norms. For some strictly positive, bounded function, we denote the weighted -norm by by

 ∥f∥2ω=∫Ωω|f|2dx,

and if thoughout , we have the equivalent

 √ω∗∥f∥≤∥f∥ω≤√ω∗∥f∥.

We will also use weighted versions of Cauchy-Schwarz. With such a weight function , we can bound a standard inner product as

 (f,g)≤∥f∥ω∥g∥ω−1.

Throughout this paper, we use short hand notations, and The letters represent positive generic constants. Their values depend on exponents, coefficients of polynomial , the spatial dimension and domain , independent of the initial data and boundary data, size of mesh and time step. These constants may be different from place to place.

## 3 A mixed finite element approximation

In this section, we will present the mixed weak formulation of the Forchheimer equation. We consider the initial boundary value problem (IVBP) associated with (LABEL:utporo):

 ⎧⎪ ⎪ ⎪⎨⎪ ⎪ ⎪⎩m+K(|∇ρ|)∇ρ=0 in Ω×(0,T),ϕρt+∇⋅m=f in Ω×(0,T),m⋅ν=ψ(x,t) in Γ×(0,T),ρ(x,0)=ρ0(x) in Γ×(0,T). \hb@xt@.01(3.1)

where is a outer normal vector of boundary , and are smooth functions .

Assume that and for all then the system (LABEL:rho:eq) reduces to the equation form

 ρt−∇⋅(ϕ−1K(|∇ρ|)∇ρ)+∇ϕ−1K(|∇ρ|)∇ρ−f(x,t)=0.

This equation is a nonlinear degenerate parabolic equation as the density gradient approaches to infinity. The existence and theory of regularity for degenerate parabolic of this type is studied in [25, 31, 11, 26, 12].

Define and the space

 M=H(div,Ω)={m∈(L2(Ω))d,∇⋅m∈L2(Ω)}

with the norm defined by

The variational formulation of (LABEL:rho:eq) is defined as follows: Find such that

 (m,z)+(K(|∇ρ|)∇ρ,z)=0,z∈M,(ϕρt,r)−(m,∇r)=(f,r)−⟨ψ,r⟩,r∈R \hb@xt@.01(3.2)

with

Let be a family of quasiuniform triangulations of with being the maximum diameter of the mesh elements. Let be the space of discontinuous piecewise polynomials of degree over . Let be the mixed element spaces approximating the space .

For density, we use the standard -projection operator, see in [9], , satisfying

 (πρ−ρ,r)=0,ρ∈R,∀r∈Rh.

This projection has well-known approximation properties, e.g. [6, 27, 5].

• For all , there is a positive constant such that

 ∥πρ∥s≤C0∥ρ∥s. \hb@xt@.01(3.3)
• There exists a positive constant such that for all ,

 ∥πρ−ρ∥0,q≤C1hs∥ρ∥s,q,0≤s≤k+1,1≤q≤∞. \hb@xt@.01(3.4)

When , in short hand we write (LABEL:prjpi) as

 ∥πρ−ρ∥≤C1hs∥ρ∥s.

The semidiscrete formulation of (LABEL:weakform) can read as follows: Find a pair such that

 (mh,z)+(K(|∇ρh|)∇ρh,z)=0,∀z∈Mh,(ϕρh,t,r)−(mh,∇r)=(f,r)−⟨ψ,r⟩,∀r∈Rh \hb@xt@.01(3.5)

with initial data .

Let be the uniform partition of with , for time step . We define . The discrete time mixed finite element approximation to (LABEL:weakform) is defined as follows: For given and . Find a pair in , such that

 (mih,z)+(K(|∇ρih|)∇ρih,z)=0,∀z∈Mh,(ϕρih−ρi−1htΔt,r)−(mih,∇r)=(fi,r)−⟨ψi,r⟩,∀r∈Rh. \hb@xt@.01(3.6)