ON THE SOLVABILITY OF THE BRINKMAN-FORCHHEIMER-EXTENDED DARCY EQUATION

# On the Solvability of the Brinkman-Forchheimer-Extended Darcy Equation

Piotr Skrzypacz111 Dr. Piotr Skrzypacz, School of Science and Technology, Nazarbayev University, 53 Kabanbay Batyr Ave., Astana 010000 Kazakhstan, Email: piotr.skrzypacz@nu.edu.kz  and  Dongming Wei222 Dr. Dongming Wei, School of Science and Technology, Nazarbayev University, 53 Kabanbay Batyr Ave., Astana 010000 Kazakhstan, Email: dongming.wei@nu.edu.kz
###### Abstract

The nonlinear Brinkman-Forchheimer-extended Darcy equation is used to model some porous medium flow in chemical reactors of packed bed type. The results concerning the existence and uniqueness of a weak solution are presented for nonlinear convective flows in medium with nonconstant porosity and for small data. Furthermore, the finite element approximations to the flow profiles in the fixed bed reactor are presented for several Reynolds numbers at the non-Darcy’s range.

2010 Mathematics Subject Classification (MSC):  76D03, 35Q35

Keywords:  Brinkman-Forchheimer Equation, Packed Bed Reactors, Existence and Uniqueness of Solution

## 1 Introduction

In this section we introduce the mathematical model describing incompressible isothermal flow in porous medium without reaction. The considered equations for the velocity and pressure fields are for flows in fluid saturated porous media. Most of research results for flows in porous media are based on the Darcy equation which is considered to be a suitable model at a small range of Reynolds numbers. However, there are restrictions of Darcy equation for modeling some porous medium flows, e.g. in closely packed medium, saturated fluid flows at slow velocity but with relatively large Reynolds numbers. The flows in such closely packed medium behave nonlinearly and can not be modelled accurately by the Darcy equation which is linear. The deficiency can be circumvented with the Brinkman–Forchheimer-extended Darcy law for flows in closely packed media, which leads to the following model: Let , , represent the reactor channel. We denote its boundary by . The conservation of volume-averaged values of momentum and mass in the packed reactor reads as follows

 −div(εν∇u−εu⊗u)+εϱ∇p+σ(u)=fin% Ω,div(εu)=0inΩ, (1)

where denote the unknown velocity and pressure, respectively. The positive quantity stands for porosity which describes the proportion of the non-solid volume to the total volume of material and varies spatially in general. The expression represents the friction forces caused by the packing and will be specified later on. The right-hand side represents an outer force (e.g. gravitation), the constant fluid density and the constant kinematic viscosity of the fluid, respectively. The expression symbolizes the dyadic product of with itself.

The formula given by Ergun [3] will be used to model the influence of the packing on the flow inertia effects

 σ(u)=150ν(1−ε)2ε2d2pu+1.751−εεdpu|u|. (2)

Thereby stands for the diameter of pellets and denotes the Euclidean vector norm. The linear term in (2) accounts for the head loss according to Darcy and the quadratic term according to Forchheimer law, respectively. For the derivation of the equations, modelling and homogenization questions in porous media we refer to e.g. [2, 4]. To close the system (1) we prescribe Dirichlet boundary condition

 u⏐Γ=g, (3)

whereby

 ∫Γiεg⋅nds=0 (4)

has to be fulfilled on each connected component of the boundary . We remark that in the case of polygonally bounded domain the outer normal vector has jumps and thus the above integral should be replaced by a sum of integrals over each side of . The distribution of porosity is assumed to satisfy the following bounds

 0<ε0≤ε(x)≤ε1≤1∀x∈Ω, (A1)

with some constants .

A comprehemsive account of fluid flows through porous media beyond the Darcy law’s valid regimes and classified by the Reynolds number, can be found in, e.g., [10]. Also, see [11] for simulating pumped water levels in abstraction boreholes using such nonlinear Darcy-Forchheimer law, and [12], [13], and [14] for recent referenes on this model.

In the next section we use the porosity distribution which is estimated for packed beds consisting of spherical particles and takes the near wall channelling effect into account. This kind of porosity distribution obeys assumption (A1).

Let us introduce dimensionless quantities

 u∗=uU0,p∗=pϱU20,x∗=xdp,g∗=gU0,

whereby denotes the magnitude of some reference velocity. For simplicity of notation we omit the asterisks. Then, the reactor flow problem reads in dimensionless form as follows

 ⎧⎪ ⎪ ⎪ ⎪⎨⎪ ⎪ ⎪ ⎪⎩−div(εRe∇u−εu⊗u)+ε∇p+αReu+βu|u|=finΩ,div(εu)=0inΩ,u=gonΓ, (5)

where

 α(x)=150κ2(x),β(x)=1.75κ(x) (6)

with

 κ(x)=1−ε(x)ε(x), (7)

and the Reynolds number is defined by

 Re=U0dpν.

The existence and uniqueness of solution of the nonlinear model (5) with constant porosity and without the convective term has been established in [5]. We will extend this result to the case when the porosity depends on the location and with the convective term in this work.

###### Remark 1

(5) becomes a Navier-Stokes problem if .

Notation  Throughout the work we use the following notations for function spaces. For , and bounded subdomain let be the usual Sobolev space equipped with norm . If , we denote the Sobolev space by and use the standard abbreviations and for the norm and seminorm, respectively. We denote by the space of functions with compact support contained in . Furthermore, stands for the closure of with respect to the norm . The counterparts spaces consisting of vector valued functions will be denoted by bold faced symbols like or . The inner product over and will be denoted by and , respectively. In the case the domain index will be omitted. In the following we denote by the generic constant which is usually independent of the model parameters, otherwise dependences will be indicated.

## 2 Existence and uniqueness results

In the following the porosity is assumed to belong to . We start with the weak formulation of problem (5) and look for its solution in suitable Sobolev spaces.

### 2.1 Variational formulation

Let

 L20(Ω):={v∈L2(Ω):(v,1)=0}

be the space consisting of functions with zero mean value. We define the spaces

 X:=H1(Ω),X0:=H10(Ω),Q:=L2(Ω),M:=L20(Ω),

and

 V:=X0×M.

Let us introduce the following bilinear forms

 a:X×X →R, a(u,v) =1Re(ε∇u,∇v), b:X×Q →R, b(u,q) =(div(εu),q), c:X×X →R, c(u,v) =1Re(αu,v).

Furthermore, we define the semilinear form

 d:X×X×X→R,d(w;u,v)=(β|w|u,v),

and trilinear form

 n:X×X×X→R,n(w,u,v)=((εw⋅∇)u,v).

We set

 A(w;u,v):=a(u,v)+c(u,v)+n(w,u,v)+d(w;u,v).

Multiplying momentum and mass balances in (5) by test functions and , respectively, and integrating by parts implies the weak formulation:
Find with  such that

 A(u;u,v)−b(v,p)+b(u,q)=(f,v)∀(v,q)∈V. (8)

First, we recall the following result from [6]:

###### Theorem 2

The mapping is an isomorphism from onto itself and from onto itself. It holds for all

 ∥εu∥1≤C{ε1+|ε|1,3}∥u∥1and∥∥uε∥∥1≤C{ε−10+ε−20|ε|1,3}∥u∥1.

In the following the closed subspace of defined by

 W={w∈H10(Ω):b(w,q)=0∀q∈L20(Ω)}.

will be employed. Next, we establish and prove some properties of trilinear form and nonlinear form .

###### Lemma 3

Let and with and . Then we have

 n(w,u,v)=−n(w,v,u). (9)

Furthermore, the trilinear form and the nonlinear form are continuous, i.e.

 |n(u,v,w)|≤Cε∥u∥1∥v∥1∥w∥1∀u,v,w∈H1(Ω), (10)
 |d(u,v,w)|≤Cε∥u∥1∥v∥1∥w∥1∀u,v,w∈H1(Ω), (11)

and for and for a sequence with , we have also

 limk→∞n(uk,uk,v)=n(u,u,v)∀v∈W. (12)

Proof. We follow the proof of [7, Lemma 2.1, §2, Chapter IV] and adapt it to the trilinear form

 n(w,u,v)=((εw⋅∇)u,v)=n∑i,j=1(εwj∂jui,vi),

which has the weighting factor . Hereby, symbols with subscripts denote components of bold faced vectors, e.g. . Let , and . Integrating by parts and employing density argument, we obtain immediately (9)

 n∑i,j=1(εwj∂jui,vi)=−n∑i,j=1(∂j(εwjvi),ui)+n∑i,j=1⟨εwjnjui,vi⟩=−n∑i,j=1(εwj∂jvi,ui)−(div(εw)u,v)+⟨(εw⋅n)u,v⟩=−n(w,v,u).

From Sobolev embedding (see [1]) and Hölder inequality follows

 ∣∣(εwj∂jui,vi)∣∣≤|ε|0,∞∥wj∥0,4∥∂jui∥0∥vi∥0,4≤C|ε|0,∞∥wj∥1|ui|1∥vi∥1,

and consequently the proof of (10) is completed. Since and , the continuity estimate (10) implies

 limk→∞n(uk,uk,v)=−limk→∞n(uk,v,uk)=−limk→∞n∑i,j=1(εukj∂jvki,uki)=−n∑i,j=1(εuj∂jvi,ui)=−n(u,v,u)=n(u,u,v).

The continuity of follows from Hölder inequality and Sobolev embedding (see [1])

 |d(u;v,w)|≤|β|∞∥u∥0,4∥v∥0,4∥w∥0≤Cε∥u∥1∥v∥1∥w∥1.

In the next stage we consider the difficulties caused by prescribing the inhomogeneous Dirichlet boundary condition. Analogous difficulties are already encountered in the analysis of Navier–Stokes problem. We will carry out the study of three dimensional case. The extension in two dimensions can be constructed analogously. Since , we can extend inside of in the form of

 g=ε−1curlh

with some . The operator curl is defined then as

 curlh=(∂2h3−∂3h2,∂3h1−∂1h3,∂1h2−∂2h1).

We note that in the two dimensional case the vector potential can be replaced by a scalar function and the operator curl is then redefined as . Our aim is to adapt the extension of Hopf (see [8]) to our model. We recall that for any parameter there exists a scalar function such that

 ∙φμ=1 in some neighborhood of Γ (depending on μ),∙φμ(x)=0 if dΓ(x)≥2exp(−1/μ)\,, where dΓ(x):=infy∈Γ|x−y|denotes the distance of x% to Γ,∙|∂jφμ(x)|≤μ/dΓ(x)  if~{}~{}dΓ(x)<2exp(−1/μ)\,, j=1,…,n.⎫⎪ ⎪ ⎪ ⎪ ⎪⎬⎪ ⎪ ⎪ ⎪ ⎪⎭ (Ex)

For the construction of see also [7, Lemma 2.4, §2, Chapter IV].
Let us define

 gμ:=ε−1curl(φμh). (13)

In the following lemma we establish bounds which are crucial for proving existence of velocity.

###### Lemma 4

The function satisfies the following conditions

 div(εgμ)=0,gμ⏐Γ=g∀μ>0, (14)

and for any there exists sufficiently small such that

 |d(u+gμ;gμ,u)| ≤δ∥β∥0,∞|u|1(|u|1+∥gμ∥0)∀u∈X0, (15) |n(u,gμ,u)| ≤δ|u|21∀u∈W. (16)

Proof. The relations in (14) are obvious. We follow [5] in order to show (15). Since Sobolev’s embedding theorem implies , so we get according to the properties of in (Ex) the following bound

 |gμ|≤Cε−10{|∇h|+μdΓ(x)|h|}≤C{μdΓ(x)+|∇h|}.

Defining

 Ωμ:={x∈Ω:dΓ(x)<2exp(−1/μ)}

we obtain from Cauchy-Schwarz and triangle inequalities

 |(β|u+gμ|,gμ⋅u)|≤∥β∥0,∞∥u∥0∥u⋅gμ∥0,Ωμ+∥β∥0,∞∥gμ∥0∥u⋅gμ∥0,Ωμ, (17)
 ∥u⋅gμ∥20,Ωμ≤∫Ωμ|u|2|gμ|2dx≤C∫Ωμ|u|2{(μ/dΓ(x))2+2μ/dΓ(x)|∇h|+|∇h|2}dx≤C{μ2∥u/dΓ∥20,Ωμ+2μ∥u/dΓ∥0,Ωμ∥u∥0,4,Ωμ∥∥|∇h|∥∥0,4,Ωμ+∥u∥20,4,Ωμ∥∥|∇h|∥∥20,4,Ωμ}≤C{μ∥u/dΓ∥0,Ωμ+∥u∥0,4∥∥|∇h|∥∥0,4,Ωμ}2,

and consequently

 ∥u⋅gμ∥0,Ωμ≤C{μ∥u/dΓ∥0,Ωμ+∥u∥0,4∥∥|∇h|∥∥0,4,Ωμ}. (18)

Applying Hardy inequality (see [1])

 ∥v/dΓ∥0≤C|v|1∀v∈H10(Ω)

and using Sobolev embedding , estimate (18) becomes

 ∥u⋅gμ∥0,Ωμ≤Cλ(μ)∥u∥1, (19)

where

 λ(μ):=max{μ,∥∥|∇h|∥∥0,4,Ωμ}.

From (17), (19), Poincaré inequality and from the fact that we conclude that for any we can choose sufficiently small such that

 |(β|u+gμ|gμ,u)|≤δ∥β∥0,∞|u|1(|u|1+∥gμ∥0)

holds. Therefore the proof of estimate (15) is completed. Now, we take a look at the trilinear convective term

 n(u,gμ,u)=((εu⋅∇)gμ,u)Ωμ=((εu⋅∇){ε−1curl(φμh)},u)Ωμ=((u⋅∇){curl(φμh)},u)Ωμ−((u⋅∇ε)gμ,u)Ωμ.

The first term of above difference becomes small due to [7, Lemma 2.3, §2, Chapter IV], and it satisfies

 ∣∣((u⋅∇){curl(φμh)},u)Ωμ∣∣=∣∣((u⋅∇)(εgμ),u)Ωμ∣∣≤δ|u|21 (20)

as long as is chosen sufficiently small. Using Hölder inequality, Sobolev embedding yields

 ∣∣((u⋅∇ε)gμ,u)Ωμ∣∣≤C∥ε∥1,3∥gμ⋅u∥0∥u∥1,

which together with (19) implies for sufficiently small the bound

 ∣∣((u⋅∇ε)gμ,u)Ωμ∣∣≤δ|u|21. (21)

From (20) and (21) follows the desired estimate (16).
While the general framework for linear and non-symmetric saddle point problems can be found in [6], our problem requires more attention due to its nonlinear character. Setting , the weak formulation (8) is equivalent to the following problem
Find such that

 A(w+gμ;w+gμ,v)−b(v,p)+b(w+gμ,q)=(f,v)∀(v,q)∈V. (22)

Let us define the nonlinear mapping with

 [G(w),v]:=a(w+gμ,v)+c(w+gμ,v)−(f,v)+n(w+gμ,w+gμ,v)+d(w+gμ;w+gμ,v), (23)

whereby defines the inner product in via . Then, the variational problem (22) reads in the space of -weighted divergence free functions as follows
Find such that

 [G(w),v]=0∀v∈W. (24)

### 2.2 Solvability of nonlinear saddle point problem

We start our study of the nonlinear operator problem (24) with the following lemma.

###### Lemma 5

The mapping defined in (23) is continuous and there exists such that

 [G(u),u]>0∀u∈Wwith|u|1=r. (25)

Proof. Let be a sequence in with . Then, applying Cauchy–Schwarz inequality and (16), we obtain for any

 ∣∣[G(uk)−G(u),v]∣∣≤1Re∣∣(ε∇(uk−u),∇v)∣∣+1Re∣∣(α(uk−u),v)∣∣+∣∣(β|uk+gμ|(uk−u),v)∣∣+∣∣(β(|uk+gμ|−|u+gμ|)(u+gμ),v)∣∣+∣∣n(uk,uk,v)−n(u,u,v)∣∣+∣∣n(uk−u,gμ,v)∣∣+∣∣n(gμ,uk−u,,v)∣∣≤ε1Re|uk−u|1|v|1+1Re∥α∥0,∞∥uk−u∥0∥v∥0+∥β∥0,∞∥uk+gμ∥0,4∥uk−u∥0∥v∥0,4+∥β∥0,∞∥u+gμ∥0,4∥uk−u∥0∥v∥0,4+∣∣n(uk,uk,v)−n(u,u,v)∣∣+C∥uk−u∥1∥gμ∥1∥v∥1.

The boundedness of in , (12), the Poincaré inequality, and the above inequality imply that

Thus, employing

we state that is continuous. Now, we note that for any we have

 [G(u),u]=1Re(ε∇(u+gμ),∇u)+1Re(α(u+gμ),u)+(β|u+gμ|(u+gμ),u)+n(u+gμ,u+gμ,u)−(f,u)≥ε0Re|u|21−ε1Re|(∇gμ,∇u)|+1Re(αu,u)−1Re|(αgμ,u)|+(β|u+gμ|,|u|2)−∣∣(β|u+gμ|gμ,u)∣∣+n(u,gμ,u)+n(gμ,gμ,u)−∥f∥0∥u∥0≥ε0Re|u|21−ε1Re|gμ|1|u|1−1Re∥α∥0,∞∥gμ∥0∥u∥0−∣∣(β|u+gμ|gμ,u)∣∣−∣∣n(u,gμ,u)∣∣−C∥gμ∥21∥u∥1−∥f∥0∥u∥0. (26)

From the Poincaré inequality, we infer the estimate

 ∥v∥1≤C|v|1∀v∈H10(Ω),

which together with (15), (16) and (26) results in

 [G(u),u]≥{ε0Re−δ(1+∥β∥0,∞)}|u|21−{ε1Re|gμ|1+C11Re∥α∥0,∞∥gμ∥0+δ∥β∥0,∞∥gμ∥0+C2∥gμ∥21+C3∥f∥0}|u|1.

Choosing such that

 0<δ<δ0:=ε0Re(1+∥β∥0,∞)−1,

and with

 (27)

leads to the desired assertion (25).
The following lemma plays a key role in the existence proof.

###### Lemma 6

Let be finite-dimensional Hilbert space with inner product inducing a norm , and be a continuous mapping such that

 [T(x),x]>0for∥x∥=r0>0.

Then there exists , with , such that

 T(x)=0.

Proof. See [9].
Now we are able to prove the main result concerning existence of velocity.

###### Theorem 7

The problem (24) has at least one solution .

Proof. We construct the approximate sequence of Galerkin solutions. Since the space is separable, there exists a sequence of linearly independent elements . Let be the finite dimensional subspace of with

 Xm:=span{wi, i=1,…,m}

and endowed with the scalar product of . Let  , be a Galerkin solution of (24) defined by

 [G(um),wj]=0,∀j=1,…,m. (28)

From Lemma 5 and Lemma 6 we conclude that

 [G(um),w]=0∀w∈Xm (29)

has a solution . The unknown coefficients can be obtained from the algebraic system (28). On the other hand, multiplying (28) by , and adding the equations for we have

 0=[G(um),um]≥{1Re−δ(1+∥β∥0,∞)}|um|21−{1Re|gμ|1+C11Re∥α∥0,∞∥gμ∥0+δ∥β∥0,∞∥gμ∥0+C2∥gμ∥21+C3∥f∥0}|um|1.

This gives together with (27) the uniform boundedness in

 |um|1≤r0,

therefore there exists and a subsequence ( we write for the convenience instead of ) such that

 um⇀uinW.

Furthermore, the compactness of embedding implies

 um→uinL4(Ω).

Taking the limit in (29) with we get

 [G(u),w]=0∀w∈Xm. (30)

Finally, we apply the continuity argument and state that (30) is preserved for any , therefore is the solution of (24).
For the reconstruction of the pressure we need inf-sup-theorem

###### Theorem 8

Assume that the bilinear form satisfies the inf-sup condition

 infq∈Msupv∈X0b(v,q)|v|1∥q∥0≥γ>0. (31)

Then, for each solution of the nonlinear problem (24) there exists a unique pressure such that the pair is a solution of the homogeneous problem (22).

Proof. See [7, Theorem 1.4, §1, Chapter IV].
We end up this subsection by proving the existence of the pressure.

###### Theorem 9

Let be solution of problem (24). Then, there exists unique pressure .

Proof. We verify the inf-sup condition (31) of Theorem 8 by employing the isomorphism of Theorem 2. From [7, Corollary 2.4, §2, Chapter I] follows that for any in there exists in such that

 (divv,q)≥γ∗∥v∥1∥q∥0

with a positive constant . Setting and applying the isomorphism in Theorem 2, we obtain the estimate

 b(u,q)=(divv,q)≥γ∗∥v∥1∥q∥0≥γε∥u∥1∥q∥0

where . From the above estimate we conclude the inf-sup condition (31).

### 2.3 Uniqueness of weak solution

We exploit a priori estimates in order to prove uniqueness of weak velocity and pressure.

###### Theorem 10

If , are sufficiently small, then the solution of (24) is unique.

Proof. Assume that and are two different solutions of (22). From (9) in Lemma 3 we obtain . Then, we obtain

 0=[G(u1)−G(u2),u1−u2]=a(u1−u2,u1−u2)+c(u1−u2,u1−u2)−(f,u1−u2)+n(u1+gμ,u1+gμ,u1−u2)−n(u2+gμ,u2+gμ,u1−u2)+(β|u1+gμ|(u1+gμ),u1−u2)−(β|u2+gμ|(u2+gμ),u1−u2)≥ε0Re|u1−u2|21−∥f∥−1∥u1−u2∥1+n(u1−u2,u2+gμ,u1−u2)+(β|u1+gμ|(u1−u2),u1−u2)+(β(|u1+gμ|−|u2+gμ|)(u2+gμ),u1−u2)≥ε0Re|u1−u2|21−∥f∥−1∥u1−u2∥1