# Solvability of the Brinkman-Forchheimer-Darcy Equation

Authors: Piotr Skrzypacz; Dongming Wei
Journal: Journal of Applied Mathematics
 (2017)

---

## Abstract

The nonlinear Brinkman-Forchheimer-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 variable 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.


---

## Body

## 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. These problems are of importance for example in oil reservoir optimization; see [1]. Most 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 modelling some porous medium flows; that is, in closely packed media, saturated fluid flows at slow velocity but with relatively large Reynolds numbers. The flows in such closely packed medium behave nonlinearly and cannot be modelled accurately by the Darcy equation which is linear. The deficiency can be circumvented with the Brinkman–Forchheimer-Darcy law for flows in closely packed media, which leads to the following model: let Ω⊂Rn, n=2,3, 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(1)-divεν∇u-εu⊗u+εϱ∇p+σu=finΩ,divεu=0inΩ,where u:Ω→Rn and p:Ω→R denote the unknown velocity and pressure, respectively. The positive quantity ε=ε(x) stands for porosity which describes the proportion of the nonsolid volume to the total volume of material and varies spatially in general. The expression σ(u) represents the friction forces caused by the packing of spherical particles with one constant diameter. The right hand side f represents an outer force (e.g., gravitation), ϱ the constant fluid density, and ν the constant kinematic viscosity of the fluid, respectively. The expression u⊗u symbolizes the dyadic product of u with itself.The formula given by Ergun [2] will be used to model the influence of the packing on the flow inertia effects(2)σu=150ν1-ε2ε2dp2u+1.751-εεdpuu.Thereby dp 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 is the Forchheimer law. The model (1)-(2) with the pressure term of the form ∇(εp) has been proposed in [3, Section 2.2]. For the derivation of the equations, their limitations, modelling, and homogenization questions in porous media we refer to [4–9] and [10, Chapter 1].To close the system (1) we prescribe Dirichlet boundary condition(3)uΓ=g,where(4)∫Γiεg·nds=0has to be fulfilled on each connected component Γi of the boundary Γ. The distribution of porosity ε is assumed to satisfy the boundsA10<ε0≤εx≤ε1≤1∀x∈Ω,with some constants 0<ε0,ε1≤1.A comprehensive account of fluid flows through porous media beyond Darcy law’s valid regimes and classified by the Reynolds number can be found in, for example, [11]. Also, see [12] for simulating pumped water levels in abstraction boreholes using such nonlinear Darcy-Forchheimer law and [13–15] for recent references 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 channeling effect into account. This kind of porosity distribution obeys assumptionA1.Let us introduce dimensionless quantities(5)u∗=uU0,p∗=pϱU02,x∗=xdp,g∗=gU0,where U0 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:(6)-divεRe∇u-εu⊗u+ε∇p+αReu+βuu=finΩ,divεu=0inΩ,u=gonΓ,where(7)αx=150κ2x,βx=1.75κx,with(8)κx=1-εxεx,and the Reynolds number is defined by (9)Re=U0dpν.The existence and uniqueness of the solution of a system for flow, temperature, and solute transport containing the nonlinear flow model (6) but with constant porosity and without the convective term have been established in [16]. We will extend this result for the flow component to the case (6) when the variable porosity depends on the location and the convective term is included. The recent existence results for the linear Brinkman problem with spatially varying α can be found in [17].Remark 1.
Equation (6) becomes a Navier-Stokes problem if ε≡1.Notation. Throughout the work we use the following notations for function spaces. For m∈N0, p≥1 and bounded subdomain G⊂Ω; let Wm,p(G) be the usual Sobolev space equipped with norm ∥·∥m,p,G. If p=2, we denote the Sobolev space by Hm(G) and use the standard abbreviations ∥·∥m,G and |·|m,G for the norm and seminorm, respectively. We denote by DG the space of C∞(G) functions with compact support contained in G. Furthermore, H0mG stands for the closure of D(G) with respect to the norm ∥·∥m,G. The counterparts spaces consisting of vector valued functions will be denoted by bold faced symbols like HmG≔HmGn or D(G)≔DGn. The L2 inner product over G⊂Ω and ∂G⊂∂Ω will be denoted by ·,·G and ·,·∂G, respectively. In the case G=Ω the domain index will be omitted. In the following we denote by C the generic constant which is usually independent of the model parameters; otherwise dependence will be indicated.


## 2. Existence and Uniqueness Results

In the following the porosityε is assumed to belong to W1,3(Ω)∩L∞(Ω). We start with the weak formulation of problem (6) and look for its solution in suitable Sobolev spaces that are adjusted to the modified momentum and mass balances in (6) and to the smoothness of the weighting function ε.


### 2.1. Variational Formulation

Let(10)L02Ω≔v∈L2Ω:v,1=0be the space consisting of L2 functions with zero mean value. We define the spaces (11)X≔H1Ω,X0≔H01Ω,Q≔L2Ω,M≔L02Ω,V≔X0×M.Let us introduce the following bilinear forms:(12)a:X×X⟶R,au,v=1Reε∇u,∇v,b:X×Q⟶R,bu,q=divεu,q,c:X×X⟶R,cu,v=1Reαu,v.Furthermore, we define the semilinear form (13)d:X×X×X⟶R,dw;u,v=βwu,v,and trilinear form (14)n:X×X×X⟶R,nw,u,v=εw·∇u,v.We set (15)Aw;u,v≔au,v+cu,v+nw,u,v+dw;u,v.Multiplying momentum and mass balances in (6) by test functions v∈X0 and q∈M, respectively, and integrating by parts implies the weak formulation(16)Findu,p∈X×MwithuΓ=gsuch thatAu;u,v-bv,p=f,v,bu,q=0∀v∈X0,∀q∈M.First, we recall the following result from [18].Theorem 2.
The mappingu↦εu is an isomorphism from H1(Ω) onto itself and from H01(Ω) onto itself. It holds for all u∈H1(Ω) that(17)εu1≤Cε1+ε1,3u1,uε1≤Cε0-1+ε0-2ε1,3u1.In the following the closed subspace ofH01(Ω) defined by (18)W=w∈H01Ω:bw,q=0∀q∈L02Ωwill be employed. Next, we establish and prove some properties of trilinear form n(·,·,·) and nonlinear form d·;·,·.Lemma 3.
Letu,v∈H1(Ω) and w∈H1(Ω) with div(εw)=0 and w·n|Γ=0. Then we have(19)nw,u,v=-nw,v,u.Furthermore, the trilinear form n(·,·,·) and the nonlinear form d(·;·,·) are continuous; that is,(20)nu,v,w≤Cεu1v1w1∀u,v,w∈H1Ω,(21)du,v,w≤Cεu1v1w1∀u,v,w∈H1Ω,and for u∈W and for a sequence uk∈W with limk→∞uk-u0=0, we have also(22)limk→∞nuk,uk,v=nu,u,v∀v∈W.Proof.
We follow the proof of [19, Lemma 2.1, §2, Chapter IV] and adapt it to the trilinear form (23)nw,u,v=εw·∇u,v=∑i,j=1nεwj∂jui,vi,which has the weighting factor ε. Hereby, symbols with subscripts denote components of bold faced vectors, for example, u=(ui)i=1,…,n. Let u∈H1, v∈D(Ω), and w∈W. Integrating by parts and employing density argument, we obtain immediately (19). (24)∑i,j=1nεwj∂jui,vi=-∑i,j=1n∂jεwjvi,ui+∑i,j=1nεwjnjui,vi=-∑i,j=1nεwj∂jvi,ui-divεwu,v+εw·nu,v=-nw,v,u.From Sobolev embedding H1(Ω)↪L4(Ω) (see [20]) and Hölder inequality follows (25)εwj∂jui,vi≤ε0,∞wj0,4∂jui0vi0,4≤Cε0,∞wj1ui1vi1,and consequently the proof of (20) is completed. Since limk→∞uikujk-uiuj0,1=0 and ε∂jvi∈L∞(Ω), the continuity estimate (20) implies (26)limk→∞nuk,uk,v=-limk→∞nuk,v,uk=-limk→∞∑i,j=1nεujk∂jvik,uik=-∑i,j=1nεuj∂jvi,ui=-nu,v,u=nu,u,v.The continuity of d(·;·,·) follows from Hölder inequality and Sobolev embedding H1(Ω)↪L4(Ω) (see [20]). (27)du;v,w≤β∞u0,4v0,4w0≤Cεu1v1w1.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. Sinceg∈H1/2(Γ), we can extend g inside of Ω in the form of (28)g=ε-1curlh,with some h∈H2(Ω). The operator curl is defined then as (29)curlh=∂2h3-∂3h2,∂3h1-∂1h3,∂1h2-∂2h1.We note that in the two-dimensional case the vector potential h∈H2(Ω) can be replaced by a scalar function h∈H2(Ω) and the operator curl is then redefined as curlh=(∂2h,-∂1h). Our aim is to adapt the extension of Hopf (see [21]) to our model. We recall that for any parameter μ>0 there exists a scalar function φμ∈C2(Ω¯) such thatExiφμ=1insomeneighborhoodofΓdependingonμ,iiφμx=0ifdΓx≥2exp⁡-1μ,wheredΓx≔infy∈Γx-ydenotesthedistanceofxtoΓ,iii∂jφμx≤μdΓxifdΓx<2exp⁡-1μ,j=1,…,n.For the construction of φμ see also [19, Lemma 2.4, §2, Chapter IV].Let us define(30)gμ≔ε-1curlφμh.In the following lemma we establish bounds which are crucial for proving existence of velocity.Lemma 4.
The functiongμ satisfies the following conditions:(31)divεgμ=0,gμΓ=g,∀μ>0,and for any δ>0 there exists sufficiently small μ>0 such that(32)du+gμ;gμ,u≤δβ0,∞u1u1+gμ0∀u∈X0,(33)nu,gμ,u≤δu12∀u∈W.Proof.
The relations in (31) are obvious. We follow [16] in order to show (32). Since h∈H2(Ω) Sobolev’s embedding theorem implies h∈L∞(Ω), so we get according to the properties of φμ in Ex the following bound: (34)gμ≤Cε0-1∇h+μdΓxh≤CμdΓx+∇h.Defining (35)Ωμ≔x∈Ω:dΓx<2exp⁡-1μ,we obtain from Cauchy-Schwarz and triangle inequalities(36)βu+gμ,gμ·u≤β0,∞u0u·gμ0,Ωμ+β0,∞gμ0u·gμ0,Ωμ,(37)u·gμ0,Ωμ2≤∫Ωμu2gμ2dx≤C∫Ωμu2μdΓx2+2μdΓx∇h+∇h2dx≤Cμ2udΓ0,Ωμ2+2μudΓ0,Ωμu0,4,Ωμ∇h0,4,Ωμ+u0,4,Ωμ2∇h0,4,Ωμ2≤CμudΓ0,Ωμ+u0,4∇h0,4,Ωμ2,and consequently(38)u·gμ0,Ωμ≤CμudΓ0,Ωμ+u0,4∇h0,4,Ωμ.Applying Hardy inequality (see [20]) (39)vdΓ0≤Cv1∀v∈H01Ω,and using Sobolev embedding H1(Ω)↪L4(Ω), estimate (38) becomes(40)u·gμ0,Ωμ≤Cλμu1,where (41)λμ≔max⁡μ,∇h0,4,Ωμ.From (36) and (40), Poincaré inequality, and the fact that limμ→0λμ=0 we conclude that for any δ>0 we can choose sufficiently small μ>0 such that (42)βu+gμgμ,u≤δβ0,∞u1u1+gμ0,holds. Therefore the proof of estimate (32) is completed. Now, we take a look at the trilinear convective term (43)nu,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 [19, Lemma 2.3, §2, Chapter IV], and it satisfies(44)u·∇curlφμh,uΩμ=u·∇εgμ,uΩμ≤δu12,as long as μ>0 is chosen sufficiently small. Using Hölder inequality, Sobolev embedding H1(Ω)↪L6(Ω) yields (45)u·∇εgμ,uΩμ≤Cε1,3gμ·u0u1,which together with (40) implies for sufficiently small μ>0 the bound(46)u·∇εgμ,uΩμ≤δu12.From (44) and (46) follows the desired estimate (33).While the general framework for linear and nonsymmetric saddle point problems can be found in [18], our problem requires more attention due to its nonlinear character. Setting w≔u-gμ, the weak formulation (16) is equivalent to the following problem:(47)Findw,p∈Vsuch  thatAw+gμ;w+gμ,v-bv,p+bw+gμ,q=f,v∀v,q∈V.Let us define the nonlinear mapping G:W→W with(48)Gw,v≔aw+gμ,v+cw+gμ,v-f,v+nw+gμ,w+gμ,v+dw+gμ;w+gμ,v,where [·,·] defines the inner product in W via [u,v]≔(∇u,∇v). Then, the variational problem (47) reads in the space W as follows.Findw∈W such that(49)Gw,v=0∀v∈W.


### 2.2. Solvability of Nonlinear Saddle Point Problem

We start our study of the nonlinear operator problem (49) with the following lemma.Lemma 5.
The mappingG defined in (48) is continuous and there exists r>0 such that(50)Gu,u>0∀u∈Wwithu1=r.Proof.
Let(uk)k∈N be a sequence in W with limk→∞∥uk-u∥1=0. Then, applying Cauchy-Schwarz inequality and (33), we obtain for any v∈W(51)Guk-Gu,v≤1Reε∇uk-u,∇v+1Reαuk-u,v+βuk+gμuk-u,v+βuk+gμ-u+gμu+gμ,v+nuk,uk,v-nu,u,v+nuk-u,gμ,v+ngμ,uk-u,v≤ε1Reuk-u1v1+1Reα0,∞uk-u0v0+β0,∞uk+gμ0,4uk-u0v0,4+β0,∞u+gμ0,4uk-u0v0,4+nuk,uk,v-nu,u,v+Cuk-u1gμ1v1.The boundedness of uk in W, (22), the Poincaré inequality, and the above inequality imply that (52)Guk-Gu,v⟶0ask⟶∞∀v∈W.Thus, employing (53)Guk-Gu1=supv∈Wv≠0⁡Guk-Gu,vv1,we state that G is continuous. Now, we note that for any u∈W we have(54)Gu,u=1Reε∇u+gμ,∇u+1Reαu+gμ,u+βu+gμu+gμ,u+nu+gμ,u+gμ,u-f,u≥ε0Reu12-ε1Re∇gμ,∇u+1Reαu,u-1Reαgμ,u+βu+gμ,u2-βu+gμgμ,u+nu,gμ,u+ngμ,gμ,u-f0u0≥ε0Reu12-ε1Regμ1u1-1Reα0,∞gμ0u0-βu+gμgμ,u-nu,gμ,u-Cgμ12u1-f0u0.From the Poincaré inequality, we infer the estimate (55)v1≤Cv1∀v∈H01Ω,which together with (32), (33), and (54) results in (56)Gu,u≥ε0Re-δ1+β0,∞u12-ε1Regμ1+C11Reα0,∞gμ0+δβ0,∞gμ0+C2gμ12+C3f0u1.Choosing δ such that (57)0<δ<δ0≔ε0Re1+β0,∞-1,and r>r0 with(58)r0≔ε1/Regμ1+1/ReC1α0,∞gμ0+δβ0,∞gμ0+C2gμ12+C3f0ε0/Re-δ1+β0,∞,leads to the desired assertion (50).The following lemma plays a key role in the existence proof.Lemma 6.
LetY be finite dimensional Hilbert space with inner product [·,·] inducing a norm ∥·∥ and T:Y→Y be a continuous mapping such that (59)Tx,x>0forx=r0>0.Then there exists x∈Y, with x≤r0, such that (60)Tx=0.Proof.
See [22].Now we are able to prove the main result concerning existence of velocity.Theorem 7.
The problem (49) has at least one solution u∈W.Proof.
We construct the approximate sequence of Galerkin solutions. Since the spaceW is separable, there exists a sequence of linearly independent elements wii∈N⊂W. Let Xm be the finite dimensional subspace of W with (61)Xm≔spanwi,i=1,…,m,and endowed with the scalar product of W. Let um=∑j=1m‍ajwj,aj∈R, be a Galerkin solution of (49) defined by(62)Gum,wj=0,∀j=1,…,m.From Lemmas 5 and 6 we conclude that(63)Gum,w=0∀w∈Xmhas a solution um∈Xm. The unknown coefficients aj can be obtained from the algebraic system (62). On the other hand, multiplying (62) by aj and adding the equations for j=1,…,m we have (64)0=Gum,um≥1Re-δ1+β0,∞um12-1Regμ1+C11Reα0,∞gμ0+δβ0,∞gμ0+C2gμ12+C3f0um1.This gives together with (58) the uniform boundedness in W(65)um1≤r0;therefore there exists u∈W and a subsequence mk→∞ (we write for the convenience m instead of mk) such that (66)um⇀uinW.Furthermore, the compactness of embedding H1(Ω)↪L4(Ω) implies (67)um⟶uinL4Ω.Taking the limit in (63) with m→∞ we get(68)Gu,w=0∀w∈Xm.Finally, we apply the continuity argument and state that (68) is preserved for any w∈W; therefore u is the solution of (49).For the reconstruction of the pressure we need inf-sup-theorem.Theorem 8.
Assume that the bilinear formb(·,·) satisfies the inf-sup condition(69)infq∈Msupv∈X0bv,qv1q0≥γ>0.Then, for each solution u of the nonlinear problem (49) there exists a unique pressure p∈M such that the pair (u,p)∈V is a solution of the homogeneous problem (47).Proof.
See [19, Theorem 1.4, §1, Chapter IV].We end up this subsection by proving the existence of the pressure.Theorem 9.
Letw be solution of problem (49). Then, there exists unique pressure p∈M.Proof.
We verify the inf-sup condition (69) of Theorem 8 by employing the isomorphism of Theorem 2. From [19, Corollary 2.4, Section 2, Chapter I] follows that for any q in L02(Ω) there exists v in H01(Ω) such that (70)divv,q≥γ∗v1q0,with a positive constant γ∗. Setting u=v/ε and applying the isomorphism in Theorem 2, we obtain the estimate (71)bu,q=divv,q≥γ∗v1q0≥γεu1q0,where γε=γ∗/Cε0-1+ε0-2|ε|1,3. From the above estimate we conclude the inf-sup condition (69).


### 2.3. Uniqueness of Weak Solution

We exploit a priori estimates in order to prove uniqueness of weak velocity and pressure.Theorem 10.
Ifgμ1, f-1≔sup0≠v∈H1Ωf,v/v1 are sufficiently small; then the solution of (49) is unique.Proof.
Assume thatu1,p1 and (u2,p2) are two different solutions of (47). From (19) in Lemma 3 we obtain nw,u,u=0∀w,u∈W. Then, we obtain(72)0=Gu1-Gu2,u1-u2=au1-u2,u1-u2+cu1-u2,u1-u2-f,u1-u2+nu1+gμ,u1+gμ,u1-u2-nu2+gμ,u2+gμ,u1-u2+βu1+gμu1+gμ,u1-u2-βu2+gμu2+gμ,u1-u2≥ε0Reu1-u212-f-1u1-u21+nu1-u2,u2+gμ,u1-u2+βu1+gμu1-u2,u1-u2+βu1+gμ-u2+gμu2+gμ,u1-u2≥ε0Reu1-u212-f-1u1-u21-nu1-u2,u2,u1-u2-nu1-u2,gμ,u1-u2-β0,∞u1+gμ·u1-u2,u1-u2-β0,∞u1+gμ-u2+gμ·u2+gμ,u1-u2.From Cauchy-Schwarz inequality and Sobolev embedding H1(Ω)↪L4(Ω) we deduce(73)u1+gμ·u1-u2,u1-u2≤Cu10+gμ0u1-u212,u1+gμ-u2+gμ·u2+gμ,u1-u2≤Cu20+gμ0u1-u212,and according to (20) we have(74)nu1-u2,u2,u1-u2≤Cu21u1-u212,and by (31) we can find μ such that(75)nu1-u2,gμ,u1-u2≤ε04Reu1-u212.Now, we find upper bounds for u1 and u2. Testing (47) with u results in (76)ε0Reu12≤f-1u1+ε0Regμ1u1+Cgμ0u0+Cgμ12u1+Cβ0,∞gμ0u12+Cβ0,∞gμ0,42u1.From Sobolev embedding H1(Ω)↪L4(Ω) we deduce for sufficiently small ∥gμ∥1(77)u1≤f-1+C1gμ1+C2gμ12ε0/Re-C3β0,∞gμ1≕Cgμ1,f-1.Putting (73), (74), (75), and (77) into (72) and using the inequality (78)f-1u1-u21≤ε04Reu1-u212+2Reε0f-12,we obtain(79)0≥ε02Reu1-u212-2Reε0f-12-Cgμ1,f-1β0,∞u1-u212-ε04Reu1-u212-Cgμ1,f-1u1-u212.For sufficiently small gμ1, f-1 the constant C(gμ1,f-1) in (77) gets small and consequently the right hand side of (79) is nonnegative. This implies u1=u2 and according to Theorem 9  p1-p2=0.


## 2.1. Variational Formulation

Let(10)L02Ω≔v∈L2Ω:v,1=0be the space consisting of L2 functions with zero mean value. We define the spaces (11)X≔H1Ω,X0≔H01Ω,Q≔L2Ω,M≔L02Ω,V≔X0×M.Let us introduce the following bilinear forms:(12)a:X×X⟶R,au,v=1Reε∇u,∇v,b:X×Q⟶R,bu,q=divεu,q,c:X×X⟶R,cu,v=1Reαu,v.Furthermore, we define the semilinear form (13)d:X×X×X⟶R,dw;u,v=βwu,v,and trilinear form (14)n:X×X×X⟶R,nw,u,v=εw·∇u,v.We set (15)Aw;u,v≔au,v+cu,v+nw,u,v+dw;u,v.Multiplying momentum and mass balances in (6) by test functions v∈X0 and q∈M, respectively, and integrating by parts implies the weak formulation(16)Findu,p∈X×MwithuΓ=gsuch thatAu;u,v-bv,p=f,v,bu,q=0∀v∈X0,∀q∈M.First, we recall the following result from [18].Theorem 2.
The mappingu↦εu is an isomorphism from H1(Ω) onto itself and from H01(Ω) onto itself. It holds for all u∈H1(Ω) that(17)εu1≤Cε1+ε1,3u1,uε1≤Cε0-1+ε0-2ε1,3u1.In the following the closed subspace ofH01(Ω) defined by (18)W=w∈H01Ω:bw,q=0∀q∈L02Ωwill be employed. Next, we establish and prove some properties of trilinear form n(·,·,·) and nonlinear form d·;·,·.Lemma 3.
Letu,v∈H1(Ω) and w∈H1(Ω) with div(εw)=0 and w·n|Γ=0. Then we have(19)nw,u,v=-nw,v,u.Furthermore, the trilinear form n(·,·,·) and the nonlinear form d(·;·,·) are continuous; that is,(20)nu,v,w≤Cεu1v1w1∀u,v,w∈H1Ω,(21)du,v,w≤Cεu1v1w1∀u,v,w∈H1Ω,and for u∈W and for a sequence uk∈W with limk→∞uk-u0=0, we have also(22)limk→∞nuk,uk,v=nu,u,v∀v∈W.Proof.
We follow the proof of [19, Lemma 2.1, §2, Chapter IV] and adapt it to the trilinear form (23)nw,u,v=εw·∇u,v=∑i,j=1nεwj∂jui,vi,which has the weighting factor ε. Hereby, symbols with subscripts denote components of bold faced vectors, for example, u=(ui)i=1,…,n. Let u∈H1, v∈D(Ω), and w∈W. Integrating by parts and employing density argument, we obtain immediately (19). (24)∑i,j=1nεwj∂jui,vi=-∑i,j=1n∂jεwjvi,ui+∑i,j=1nεwjnjui,vi=-∑i,j=1nεwj∂jvi,ui-divεwu,v+εw·nu,v=-nw,v,u.From Sobolev embedding H1(Ω)↪L4(Ω) (see [20]) and Hölder inequality follows (25)εwj∂jui,vi≤ε0,∞wj0,4∂jui0vi0,4≤Cε0,∞wj1ui1vi1,and consequently the proof of (20) is completed. Since limk→∞uikujk-uiuj0,1=0 and ε∂jvi∈L∞(Ω), the continuity estimate (20) implies (26)limk→∞nuk,uk,v=-limk→∞nuk,v,uk=-limk→∞∑i,j=1nεujk∂jvik,uik=-∑i,j=1nεuj∂jvi,ui=-nu,v,u=nu,u,v.The continuity of d(·;·,·) follows from Hölder inequality and Sobolev embedding H1(Ω)↪L4(Ω) (see [20]). (27)du;v,w≤β∞u0,4v0,4w0≤Cεu1v1w1.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. Sinceg∈H1/2(Γ), we can extend g inside of Ω in the form of (28)g=ε-1curlh,with some h∈H2(Ω). The operator curl is defined then as (29)curlh=∂2h3-∂3h2,∂3h1-∂1h3,∂1h2-∂2h1.We note that in the two-dimensional case the vector potential h∈H2(Ω) can be replaced by a scalar function h∈H2(Ω) and the operator curl is then redefined as curlh=(∂2h,-∂1h). Our aim is to adapt the extension of Hopf (see [21]) to our model. We recall that for any parameter μ>0 there exists a scalar function φμ∈C2(Ω¯) such thatExiφμ=1insomeneighborhoodofΓdependingonμ,iiφμx=0ifdΓx≥2exp⁡-1μ,wheredΓx≔infy∈Γx-ydenotesthedistanceofxtoΓ,iii∂jφμx≤μdΓxifdΓx<2exp⁡-1μ,j=1,…,n.For the construction of φμ see also [19, Lemma 2.4, §2, Chapter IV].Let us define(30)gμ≔ε-1curlφμh.In the following lemma we establish bounds which are crucial for proving existence of velocity.Lemma 4.
The functiongμ satisfies the following conditions:(31)divεgμ=0,gμΓ=g,∀μ>0,and for any δ>0 there exists sufficiently small μ>0 such that(32)du+gμ;gμ,u≤δβ0,∞u1u1+gμ0∀u∈X0,(33)nu,gμ,u≤δu12∀u∈W.Proof.
The relations in (31) are obvious. We follow [16] in order to show (32). Since h∈H2(Ω) Sobolev’s embedding theorem implies h∈L∞(Ω), so we get according to the properties of φμ in Ex the following bound: (34)gμ≤Cε0-1∇h+μdΓxh≤CμdΓx+∇h.Defining (35)Ωμ≔x∈Ω:dΓx<2exp⁡-1μ,we obtain from Cauchy-Schwarz and triangle inequalities(36)βu+gμ,gμ·u≤β0,∞u0u·gμ0,Ωμ+β0,∞gμ0u·gμ0,Ωμ,(37)u·gμ0,Ωμ2≤∫Ωμu2gμ2dx≤C∫Ωμu2μdΓx2+2μdΓx∇h+∇h2dx≤Cμ2udΓ0,Ωμ2+2μudΓ0,Ωμu0,4,Ωμ∇h0,4,Ωμ+u0,4,Ωμ2∇h0,4,Ωμ2≤CμudΓ0,Ωμ+u0,4∇h0,4,Ωμ2,and consequently(38)u·gμ0,Ωμ≤CμudΓ0,Ωμ+u0,4∇h0,4,Ωμ.Applying Hardy inequality (see [20]) (39)vdΓ0≤Cv1∀v∈H01Ω,and using Sobolev embedding H1(Ω)↪L4(Ω), estimate (38) becomes(40)u·gμ0,Ωμ≤Cλμu1,where (41)λμ≔max⁡μ,∇h0,4,Ωμ.From (36) and (40), Poincaré inequality, and the fact that limμ→0λμ=0 we conclude that for any δ>0 we can choose sufficiently small μ>0 such that (42)βu+gμgμ,u≤δβ0,∞u1u1+gμ0,holds. Therefore the proof of estimate (32) is completed. Now, we take a look at the trilinear convective term (43)nu,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 [19, Lemma 2.3, §2, Chapter IV], and it satisfies(44)u·∇curlφμh,uΩμ=u·∇εgμ,uΩμ≤δu12,as long as μ>0 is chosen sufficiently small. Using Hölder inequality, Sobolev embedding H1(Ω)↪L6(Ω) yields (45)u·∇εgμ,uΩμ≤Cε1,3gμ·u0u1,which together with (40) implies for sufficiently small μ>0 the bound(46)u·∇εgμ,uΩμ≤δu12.From (44) and (46) follows the desired estimate (33).While the general framework for linear and nonsymmetric saddle point problems can be found in [18], our problem requires more attention due to its nonlinear character. Setting w≔u-gμ, the weak formulation (16) is equivalent to the following problem:(47)Findw,p∈Vsuch  thatAw+gμ;w+gμ,v-bv,p+bw+gμ,q=f,v∀v,q∈V.Let us define the nonlinear mapping G:W→W with(48)Gw,v≔aw+gμ,v+cw+gμ,v-f,v+nw+gμ,w+gμ,v+dw+gμ;w+gμ,v,where [·,·] defines the inner product in W via [u,v]≔(∇u,∇v). Then, the variational problem (47) reads in the space W as follows.Findw∈W such that(49)Gw,v=0∀v∈W.


## 2.2. Solvability of Nonlinear Saddle Point Problem

We start our study of the nonlinear operator problem (49) with the following lemma.Lemma 5.
The mappingG defined in (48) is continuous and there exists r>0 such that(50)Gu,u>0∀u∈Wwithu1=r.Proof.
Let(uk)k∈N be a sequence in W with limk→∞∥uk-u∥1=0. Then, applying Cauchy-Schwarz inequality and (33), we obtain for any v∈W(51)Guk-Gu,v≤1Reε∇uk-u,∇v+1Reαuk-u,v+βuk+gμuk-u,v+βuk+gμ-u+gμu+gμ,v+nuk,uk,v-nu,u,v+nuk-u,gμ,v+ngμ,uk-u,v≤ε1Reuk-u1v1+1Reα0,∞uk-u0v0+β0,∞uk+gμ0,4uk-u0v0,4+β0,∞u+gμ0,4uk-u0v0,4+nuk,uk,v-nu,u,v+Cuk-u1gμ1v1.The boundedness of uk in W, (22), the Poincaré inequality, and the above inequality imply that (52)Guk-Gu,v⟶0ask⟶∞∀v∈W.Thus, employing (53)Guk-Gu1=supv∈Wv≠0⁡Guk-Gu,vv1,we state that G is continuous. Now, we note that for any u∈W we have(54)Gu,u=1Reε∇u+gμ,∇u+1Reαu+gμ,u+βu+gμu+gμ,u+nu+gμ,u+gμ,u-f,u≥ε0Reu12-ε1Re∇gμ,∇u+1Reαu,u-1Reαgμ,u+βu+gμ,u2-βu+gμgμ,u+nu,gμ,u+ngμ,gμ,u-f0u0≥ε0Reu12-ε1Regμ1u1-1Reα0,∞gμ0u0-βu+gμgμ,u-nu,gμ,u-Cgμ12u1-f0u0.From the Poincaré inequality, we infer the estimate (55)v1≤Cv1∀v∈H01Ω,which together with (32), (33), and (54) results in (56)Gu,u≥ε0Re-δ1+β0,∞u12-ε1Regμ1+C11Reα0,∞gμ0+δβ0,∞gμ0+C2gμ12+C3f0u1.Choosing δ such that (57)0<δ<δ0≔ε0Re1+β0,∞-1,and r>r0 with(58)r0≔ε1/Regμ1+1/ReC1α0,∞gμ0+δβ0,∞gμ0+C2gμ12+C3f0ε0/Re-δ1+β0,∞,leads to the desired assertion (50).The following lemma plays a key role in the existence proof.Lemma 6.
LetY be finite dimensional Hilbert space with inner product [·,·] inducing a norm ∥·∥ and T:Y→Y be a continuous mapping such that (59)Tx,x>0forx=r0>0.Then there exists x∈Y, with x≤r0, such that (60)Tx=0.Proof.
See [22].Now we are able to prove the main result concerning existence of velocity.Theorem 7.
The problem (49) has at least one solution u∈W.Proof.
We construct the approximate sequence of Galerkin solutions. Since the spaceW is separable, there exists a sequence of linearly independent elements wii∈N⊂W. Let Xm be the finite dimensional subspace of W with (61)Xm≔spanwi,i=1,…,m,and endowed with the scalar product of W. Let um=∑j=1m‍ajwj,aj∈R, be a Galerkin solution of (49) defined by(62)Gum,wj=0,∀j=1,…,m.From Lemmas 5 and 6 we conclude that(63)Gum,w=0∀w∈Xmhas a solution um∈Xm. The unknown coefficients aj can be obtained from the algebraic system (62). On the other hand, multiplying (62) by aj and adding the equations for j=1,…,m we have (64)0=Gum,um≥1Re-δ1+β0,∞um12-1Regμ1+C11Reα0,∞gμ0+δβ0,∞gμ0+C2gμ12+C3f0um1.This gives together with (58) the uniform boundedness in W(65)um1≤r0;therefore there exists u∈W and a subsequence mk→∞ (we write for the convenience m instead of mk) such that (66)um⇀uinW.Furthermore, the compactness of embedding H1(Ω)↪L4(Ω) implies (67)um⟶uinL4Ω.Taking the limit in (63) with m→∞ we get(68)Gu,w=0∀w∈Xm.Finally, we apply the continuity argument and state that (68) is preserved for any w∈W; therefore u is the solution of (49).For the reconstruction of the pressure we need inf-sup-theorem.Theorem 8.
Assume that the bilinear formb(·,·) satisfies the inf-sup condition(69)infq∈Msupv∈X0bv,qv1q0≥γ>0.Then, for each solution u of the nonlinear problem (49) there exists a unique pressure p∈M such that the pair (u,p)∈V is a solution of the homogeneous problem (47).Proof.
See [19, Theorem 1.4, §1, Chapter IV].We end up this subsection by proving the existence of the pressure.Theorem 9.
Letw be solution of problem (49). Then, there exists unique pressure p∈M.Proof.
We verify the inf-sup condition (69) of Theorem 8 by employing the isomorphism of Theorem 2. From [19, Corollary 2.4, Section 2, Chapter I] follows that for any q in L02(Ω) there exists v in H01(Ω) such that (70)divv,q≥γ∗v1q0,with a positive constant γ∗. Setting u=v/ε and applying the isomorphism in Theorem 2, we obtain the estimate (71)bu,q=divv,q≥γ∗v1q0≥γεu1q0,where γε=γ∗/Cε0-1+ε0-2|ε|1,3. From the above estimate we conclude the inf-sup condition (69).


## 2.3. Uniqueness of Weak Solution

We exploit a priori estimates in order to prove uniqueness of weak velocity and pressure.Theorem 10.
Ifgμ1, f-1≔sup0≠v∈H1Ωf,v/v1 are sufficiently small; then the solution of (49) is unique.Proof.
Assume thatu1,p1 and (u2,p2) are two different solutions of (47). From (19) in Lemma 3 we obtain nw,u,u=0∀w,u∈W. Then, we obtain(72)0=Gu1-Gu2,u1-u2=au1-u2,u1-u2+cu1-u2,u1-u2-f,u1-u2+nu1+gμ,u1+gμ,u1-u2-nu2+gμ,u2+gμ,u1-u2+βu1+gμu1+gμ,u1-u2-βu2+gμu2+gμ,u1-u2≥ε0Reu1-u212-f-1u1-u21+nu1-u2,u2+gμ,u1-u2+βu1+gμu1-u2,u1-u2+βu1+gμ-u2+gμu2+gμ,u1-u2≥ε0Reu1-u212-f-1u1-u21-nu1-u2,u2,u1-u2-nu1-u2,gμ,u1-u2-β0,∞u1+gμ·u1-u2,u1-u2-β0,∞u1+gμ-u2+gμ·u2+gμ,u1-u2.From Cauchy-Schwarz inequality and Sobolev embedding H1(Ω)↪L4(Ω) we deduce(73)u1+gμ·u1-u2,u1-u2≤Cu10+gμ0u1-u212,u1+gμ-u2+gμ·u2+gμ,u1-u2≤Cu20+gμ0u1-u212,and according to (20) we have(74)nu1-u2,u2,u1-u2≤Cu21u1-u212,and by (31) we can find μ such that(75)nu1-u2,gμ,u1-u2≤ε04Reu1-u212.Now, we find upper bounds for u1 and u2. Testing (47) with u results in (76)ε0Reu12≤f-1u1+ε0Regμ1u1+Cgμ0u0+Cgμ12u1+Cβ0,∞gμ0u12+Cβ0,∞gμ0,42u1.From Sobolev embedding H1(Ω)↪L4(Ω) we deduce for sufficiently small ∥gμ∥1(77)u1≤f-1+C1gμ1+C2gμ12ε0/Re-C3β0,∞gμ1≕Cgμ1,f-1.Putting (73), (74), (75), and (77) into (72) and using the inequality (78)f-1u1-u21≤ε04Reu1-u212+2Reε0f-12,we obtain(79)0≥ε02Reu1-u212-2Reε0f-12-Cgμ1,f-1β0,∞u1-u212-ε04Reu1-u212-Cgμ1,f-1u1-u212.For sufficiently small gμ1, f-1 the constant C(gμ1,f-1) in (77) gets small and consequently the right hand side of (79) is nonnegative. This implies u1=u2 and according to Theorem 9  p1-p2=0.


## 3. A Channel Flow Problem in Packed Bed Reactors

In this section, we provide an example of the flow problem in packed bed reactors with numerical solutions at small and relatively large Reynolds numbers to show the nonlinear behavior of the velocity solutions. Our numerical tests were conducted using the noncommercial object-oriented finite element package MoonMD [23] that was originally developed by the research group in Magdeburg and used for several benchmarks. The numerical results generated by MoonMD have been also verified by commercial software packages FLUENT® and FEMLAB®; see [24] and [3, Chapter 2.6].Let the reactor channel with Newtonian-fluid be represented by the plain domainΩ=(0,L)×(-R,R) where R=5 and L=60. In all computations we use the porosity distribution from [3, Section 2.2] which is determined experimentally and takes into account the effect of wall channeling in packed bed reactors(80)εx,y=εy=ε∞1+1-ε∞ε∞e-6R-y,where ε∞=0.45. The distribution of the porosity is presented in Figure 1. We distinguish between the inlet, outlet, and membrane parts of domain boundary Γ and denote them by Γin, Γout, and Γw, respectively. Let (81)Γin=x,y∈Γ:x=0,Γout=x,y∈Γ:x=L,Γw=x,y∈Γ:y=-R,y=R.At the inlet Γin and at the membrane wall Γw we prescribe Dirichlet boundary conditions, namely, the plug flow conditions (82)uΓin=uin=uin,0T,uΓw=uw=0,uwTfory=-R,0,-uwTfory=R,where uin>0 and uw≥0. We consider the case of fixed bed reactor, that is, uw=0. At the outlet Γout we set the following outflow boundary condition: (83)-1Re∂u∂n+pn=0,where n denotes the outer normal. This boundary condition results from the integration by parts when deriving the weak formulation, and it is called thedo-nothing boundary condition. In order to avoid discontinuity between the inflow and wall conditions we replace constant profile by trapezoidal one with zero value at the corners. Our computations are carried out on the Cartesian mesh which consists of 100 stretched rectangular cells on the coarse level (see Figure 2) and will be three times uniformly refined. In order to approximate the weak solution of the dimensionless system from (6) subject to the Dirichlet anddo-nothing boundary conditions, we apply biquadratic conforming and discontinuous piecewise linear finite elements for the velocity and pressure, respectively.Figure 1
Varying porosity.Figure 2
Initial mesh for reactor flow problem.The detailed finite element analysis of the Brinkman-Forchheimer-Darcy equation will be presented in the forthcoming work. The plots of velocity magnitude in fixed bed reactor (uw=0) are presented along the vertical axis x=50. In the investigated reactor the inlet velocity is assumed to be normalized (uin=1). Due to the variation of porosity we might expect higher velocity at the reactor walls Γw. This tunneling effect can be well observed in Figure 3 which shows the velocity profiles for different Reynolds numbers. This effect is not present when applying standard Brinkman equation with constant porosity. We remark that the maximum of velocity magnitude decreases with increasing Reynolds numbers in the considered cases Re=1,10,100. Our preliminary results are comparable with those ones obtained in [3, Chapter 4.3]. The weighted areas under each velocity profile u=u(y) are the same due to conservation of mass. In fact, we have according to the divergence theorem (84)∫0,50×-R,Rdivεudx=0⟹∫-RRεyuydy=∫-RRεyvindy.However, we expect increasing maximum of the velocity in the case of high Reynolds numbers. This case requires stabilized finite element methods and will be considered in the forthcoming works.Figure 3
Flow profiles in fixed bed reactor atx=50.


## 4. Conclusion

In this work, we have extended the existence and uniqueness of solution result in literature for the porous medium flow problem based on the nonlinear Brinkman-Forchheimer-extended Darcy law. The existing result is valid only for constant porosity and without the considered convection effects, and our result holds for variable porosity and it includes convective effects. We also provided a numerical solution to demonstrate the nonlinear velocity solutions at moderately large Reynolds numbers for which case the Brinkman-Forchheimer-Darcy law applies.



---

*Source: 7305230-2017-07-20.xml*