# Numerical Computation of Brinkman Flow with Stable Mixed Element Method

Authors: Xindong Li; Wenwen Xu
Journal: Mathematical Problems in Engineering
 (2019)

---

## Abstract

Herein, we discuss a mixed finite element method applied to the Brinkman equation of fluid motion in porous medium, which covers a field of situation from the Darcy equation to the Stokes problem associated with various perturbation parameters. The finite element based on staggered meshes is shown to be stable and effective with each case as the corresponding error estimates for both velocity and pressure are established. Finally, we present numerical examples confirming the theoretical analysis and the stability of the finite spaces approximation.


---

## Body

## 1. Introduction

In this paper, we consider an efficient mixed finite method to simulate the Brinkman model that describes the viscous fluid creeping flow in highly heterogeneous porous medium [1]. Firstly, we comment on its mathematical and physical backgrounds. The mathematical framework is that the parameter-dependent model viewed as a generalized Stokes problem can be reduced to be a form of Darcy’s law. The derivation of the equation could be explained in some aspect [2, 3]. It is widely used in calculating the hydrodynamic interaction and hindered transport of solid spherical particles to order or disorder fibrous medium [4] and applied to fuel cell dynamics or the boundary layer close to a solid wall [5].The Brinkman model has been studied by some researchers. The crucial key is to find uniform stable simulate as it covers the fields of both Stokes and Darcy problems. As seen in [2], it does not work as the parameter approaches zero by numerical tests for the common Stokes elements, e.g., P2-P0, Crouzeix-Raviart and Mini element, while it fails as the parameter is far away from zero for the classical Raviart-Thomas element approximating. Reference [2] proposes a new nonconforming robust finite element spaces for the unknown variables and detailed argument for the case of perturbation solution. Using this fact, Xu-Zhang [6] give a divergence-free element based on a special triangulation which could satisfy compatible condition. With different norm, [7, 8] provide a priori and a posteriori error estimates to ensure the Mini element to be uniformly stable. Reference [9] introduces unified stabilized element utilizing orthogonal subscale stabilization and algebraic subgrid scale method with two different weak formulations, both containing a term of the pressure jumps internal boundaries. Reference [10] gives a stable least-square formulation by introducing an auxiliary variable to substitute for velocity flux. Other more novel methods, such as weak Galerkin finite element method [11, 12], multigrid method [13], and virtual method [14], are applied to Brinkman model without considering parameter dependence of solution.In this paper, we will consider a stable mixed finite method for Brinkman equation associated with various perturbation parameters, which is first constructed to deal with Stokes equation [15]; furthermore, we modified it for Darcy equation and made it more efficient [16]. The remainder of this paper is organized as follows. In Section 2, we introduce the model and corresponding weak formulation. In Section 3, we discuss error estimates and prove convergence of the discrete scheme. Finally, in Section 4, some numerical experiments confirming theoretical result close the paper.Throughout this paper, the symbolC,Ci denoting generic positive constants, independent of any mesh parameters, may take different values at different occurrences. εi denotes generic small positive constant.


## 2. The Model Problem and Corresponding Weak Formulation

In the domainΩ⊂R2, we denote the velocity by u and pressure by p. The Brinkman model is described by the parametric equations(1)-t2Δu+K-1u+∇p=f,(2)∇·u=q,where K is the media permeability and f and q are given data. For simplicity, we consider the homogenous boundary conditions: for t>0, we let u=0 and u·n=0 in the limit case of t=0; here n is defined as the unit outward normal vectors of the boundary.We begin by recalling some useful preliminaries and notations. We use the classical definitions for the Sobolev spacesHk(Ω) with the usual norm and seminorm ·k and ·k. Let L2(Ω) be the space of square integrable function in Ω with inner product (·,·). We also use the Hilbert space(3)Hdiv,Ω=u∈L2Ω2;div⁡u∈L2Ω,then H01(Ω) is defined as the subspace of H1(Ω) with zero trace on ∂Ω and H0(div,Ω) is the subspace of H(div,Ω) with zero normal trace on ∂Ω. L02(Ω) is the subspace of L2(Ω) consisting of functions with zero mean value on Ω.Ift>0, we have V=(H01(Ω))2 and S=L02(Ω); in the singular case that the perturbation parameter t→0, the vector space is replaced by V=H0(div,Ω).Define two bilinear forms by(4)au,v=t2∇u,∇v+K-1u,v,bv,w=w,∇·v.The variational formulation of the model is represented by the following: find (u,p)∈V×S, such that(5)au,v-bv,p=f,v∀v∈V,bu,w=q,w∀q∈S.The associated norm on V is given by (6)uV2=t2∇u2+u2+∇·u2.For 0<t<C, we find that uV≤C1u1. The L-B-B condition holds, as(7)supv∈V⁡bv,qvV≥C1supv∈V⁡bv,qv1≥C1βq∀q∈S.Theorem 1.
There exists a unique solution(u,p)∈V×S to problem (5).Proof.
The coercivity and continuity ofa(u,v) can be obtained by simple calculation, which means(8)au,u≥C2uV2,∀u∈V0,au,v≤C3uVvV,∀u,v∈V,where V0 is a closed subspace of V via V0=v∈V∣b(v,q)=0,∀q∈S. Combining with the inf-sup condition (7), the Brezzi condition is satisfied. It is clear that problem (5) is well-posed and the solution is unique [7, 17].Next, we discuss the finite element discretization of the variational formulation (5). The finite element spaces Vh and Sh are initially proposed by Han and Wu to treat with Stokes problem [15]. We find it is stable uniformly in the parameter t∈(0,K]. In brief, we firstly recall that τh is quasiuniformly rectangular partitioning of the domain Ω; τh1 and τh2 are staggered rectangular meshes based on τh. Then, the discrete spaces are defined as(9)Vh1=vh∈C0Ω¯:vh∣T1∈Q1,1T1,∀T1∈τh1,Vh2=vh∈C0Ω¯:vh∣T2∈Q1,1T2,∀T2∈τh2,Sh≔wh:wh∣T=constant,∀T∈τh;∫Ωwhdx=0.Let (10)Vh=V⋂Vh1×Vh2.Using the subspaces Vh and Sh instead of V and S in the variational problem (5), we obtain the discrete problem: find (uh,ph)∈Vh×Sh, such that(11)auh,vh-bvh,ph=f,v∀vh∈Vh,buh,wh=q,wh∀wh∈Sh.


## 3. Convergence Analysis and Error Estimate

In this section, we give the convergence analysis and error estimate. Fort>0, the definition of the interpolating Π:V→Vh by Πu∈Vh is expressed by the formulation [15](12)∫lΠh1u1ds=∫lu1ds∀l∈LV,∫lΠh2u2ds=∫lu2ds∀l∈LH,where LV is the set that contains all vertical edges of the partition τh, while LH contains all horizontal edges. For t=0, the interpolation operator can be redefined as [18](13)∫lΠh1u1ds=∫lu1ds∀l∈LV′,∫lΠh2u2ds=∫lu2ds∀l∈LH′,where LV′ is a set of edges of elements gotten by bisecting the most bottom element edges and most top element edges of LV and LH′ are gotten by bisecting the most left element edges and most right element edges of LH. Furthermore, we denote Ps:S→Sh as the L2 projection operator.Referring to [18], we recall some properties of the interpolating Π and projection Ps.Lemma 2.
(i) There exist two constantsC3 and C4 independent of h, such that(14)u-Πu1,Ω≤C3hj-1uj,Ω,1≤j≤2.∀u∈V.(15)p-Psp≤C4hp1,Ω,∀p∈S.
(ii) Foru∈(H1(Ω))2, there exists a constant C5 independent of h, such that(16)Πu1,Ω≤C5u1,Ω∀u∈V.
(iii) For anyu∈V, we have(17)∫Ωqhdiv⁡u-Πudx=0,∀qh∈Sh.Theorem 3.
The discrete inf-sup condition is valid; namely, there is a constantβ≥0, such that(18)supvh∈Vh⁡bvh,qhvhv≥βqh,∀qh∈Sh.Proof.
Based on the fact thatvhV≤Cvh1, we find (19)supvh∈Vh⁡bvh,qhvhv≥1Cbvh,qhvh1≥1CbΠv,qhΠv1≥1CC5bv,qhv1which means(20)supvh∈Vh⁡bvh,qhvhv≥csupv∈V⁡bv,qhv1≥cβqh,∀qh∈Sh.Theorem 4.
If(u,p) satisfy (5) and (uh,ph) is the solution of (11), then the following error estimate holds: (21)u-uha+p-ph≤Chtu2+hu2+p1.Proof.
For the error, we subtract (11) from (5) by v=vh∈Vh and obtain that(22)au-uh,vh-p-ph,∇·vh=0(23)∇·u-uh,wh=0Taking vh=uh-Πu, we have(24)au-uh,uh-Πu=p-ph,∇·uh-Πu=p-qh,∇·uh-Πu,∀qh∈Sh.which is derived by combining with the equality (17) and (23). For convenience, we introduce a norm ·a for the norm on Vh, which is commonly used in the following: (25)va2=v2+t2∇v2=av,v.Then, from (24) and ε- inequality, we further have (26)uh-Πua2=au-Πu,uh-Πu-p-qh,∇·uh-Πu=au-Πu,uh-Πu-p-Psp,∇·uh-Πu-∇·Π′uh-Πu≤ε1uh-Πua2+1ε1u-Πua2+Chp-Psp0∇·uh-Πu1≤ε1uh-Πua2+1ε1u-Πua2+Chp-Psp0∇·uh-Πu0≤ε1uh-Πua2+1ε1u-Πua2+ε2p-Psp2+1ε2uh-Πu)2where ε1,ε2 are arbitrary positive constants and combine with the quality of the interpolation Π′, which is based on the fact that ∇·uh-Πu1≤C as ∂xi(∇·(uh-Πu))=c on each T∈τhi. Now, we can obtain that(27)uh-Πua≤Cu-Πua+p-Psp)≤Chtu2+hu2+p1.For Psp-ph, there exist a v∈(H01(Ω))2 and a constant c2>0 independent of Psp-ph [17], satisfying (28)∇·v=Psp-ph,v1,Ω≤c2Psp-ph.Then, we have(29)Psp-ph2=Psp-ph,∇·v=Psp-ph,∇·Πv=Psp-p,∇·Πv+au-uh,Πv≤Chp1Πv1+Cu-uhaΠv1which means(30)Psp-ph≤Chtu2+hu2+p1.Based on (16), we get that (31)p-ph≤Chtu2+hu2+p1.In general, the solution(u,p) depends on the parameter t; the conclusion above will deteriorate as t approaches zero. As the introduction in [2], we recall some notation in the following: (32)H+1=g∈H1∩L02∣∫Ωgx2x-xj2dx<∞,j=1,2,…,Nwith(33)g1,+2=g12+∑j=1∫Ωgx2x-xj2dx,where x1,…,xn denote the vertices of the domain Ω. (34)rotv=∂v2∂x1-∂v1∂x2.Theorem 5.
When the solution(u,p) of system (1)-(2) depends on t, then there is a constant C independent of t and h, such that (35)u-uha+p-ph≤Ch1/2frot+q1,+,Proof.
Now, we give some new regularity about the solution. In the caset=0, system (5) reduces to elliptical about (u0,p0). Referring to [2], we derive that(36)u01+p01≤cfrot+q0,while the solution (u,p) in case of t>0 has the new regularity(37)t2u2+tu1+u-u01+p-p01≤Ct1/2frot+q1,+.Applying with (37), we conclude that (38)u-Πu+tu-Πu1≤Ch1/2frot+q1,+.Combining with (27), (36), and (37), we have (39)u-uha≤Cu-Πua+p-Psp≤Cu-Πu+tu-Πu1+hp-p01+p01≤Ch1/2frot+q1,++hfrot+q0+t1/2frot+q1,+≤Ch1/2frot+q1,+.Finally, from (29) and (39), we have(40)p-ph≤p-Psp+Psp-ph≤Chp1+u-uha≤Ch1/2frot+q1,+.


## 4. Numerical Examples

In this section, we present some numerical results for the model problem (1)-(2).For simplicity, we assume that the domain is a unit square:Ω=(0,1)×(0,1). Firstly, we take into account that the analytical solution is independent of t; then, we observe the convergence of the numerical error as the parameter t decreases 1 to zero. Lastly, we introduce several permeability cases to test the performance of the proposed method.We need to compute the following useful errors:uL2, ua, pL2, PspL2 and corresponding convergence rates.Example 1.
The analytical solution of the problem is as follows:(41)u=curlsin2⁡πxsin2⁡πyp=cos⁡πxcos⁡πyf=u-t2Δu-∇p.
Then, the error estimates are listed in Tables1–4, while the convergence rates are plotted with respect to h in Figure 1. Figure 2 gives the true and approximate pressure. We can find that all the error orders are optimal and consistent with the theoretical results of Theorem 4.Table 1
Numerical error ofuL2 for Example 1.
t h 2 - 2 2 - 3 2 - 4 2 - 5 2 - 6 error rate error rate error rate error rate error rate 2 - 2 3.29e-1 — 7.89e-2 2.0587 1.95e-2 2.0133 4.88e-3 2.0031 1.22e-3 2.0008 2 - 4 2.40e-1 — 4.88e-2 2.3006 1.16e-2 2.0752 2.86e-3 2.0184 7.12e-4 2.0050 2 - 6 2.33e-1 — 4.56e-2 2.3569 1.05e-2 2.1107 2.59e-3 2.0284 6.43e-4 2.0065 2 - 8 2.33e-1 — 4.56e-2 2.3572 1.05e-2 2.1144 2.56e-3 2.0377 6.34e-4 2.0134 2 - 10 2.33e-1 — 4.56e-2 2.3572 1.05e-2 2.1143 2.56e-3 2.0378 6.34e-4 2.0145 2 - 12 2.33e-1 — 4.56e-2 2.3571 1.05e-2 2.1143 2.56e-3 2.0377 6.34e-4 2.0144 0 2.37e-1 —– 4.58e-2 2.3757 1.05e-2 2.1206 2.56e-3 2.0395 6.34e-4 2.0145Table 2
Numerical error ofpL2 for Example 1.
t h 2 - 2 2 - 3 2 - 4 2 - 5 2 - 6 error rate error rate error rate error rate error rate 2 - 2 1.64e-1 — 8.17e-2 1.0092 4.03e-2 1.0193 2.01e-2 1.0061 1.00e-2 1.0016 2 - 4 1.58e-1 — 8.01e-2 0.9846 4.01e-2 0.9983 2.00e-3 0.9997 1.00e-4 0.9999 2 - 6 1.58e-1 — 8.00e-2 0.9841 4.01e-2 0.9976 2.00e-3 0.9995 1.00e-4 0.9999 2 - 8 1.58e-1 — 8.00e-2 0.9842 4.01e-2 0.9976 2.00e-3 0.9996 1.00e-4 0.9999 2 - 10 1.58e-1 — 8.00e-2 0.9842 4.01e-2 0.9976 2.00e-3 0.9995 1.00e-4 0.9999 2 - 12 1.58e-1 — 8.00e-2 0.9842 4.01e-2 0.9976 2.00e-3 0.9995 1.00e-4 0.9999 0 1.59e-1 —– 8.00e-2 0.9904 4.01e-2 0.9984 2.00e-3 0.9996 1.00e-4 0.9999Table 3
Numerical error ofua for Example 1.
t h 2 - 2 2 - 3 2 - 4 2 - 5 2 - 6 error rate error rate error rate error rate error rate 2 - 2 3.51e+0 — 1.76e+0 0.9935 8.84e-1 0.9963 4.42e-1 0.9989 2.21e-1 0.9997 2 - 4 3.14e+0 — 1.59e+0 0.9843 7.96e-1 0.9953 3.99e-1 0.9987 1.99e-1 0.9996 2 - 6 3.16e+0 — 1.59e+0 0.9926 7.92e-1 1.0027 3.96e-1 1.0007 1.98e-1 1.0002 2 - 8 3.16e+0 — 1.59e+0 0.9932 7.92e-1 1.0040 3.96e-1 1.0014 1.98e-1 1.0004 2 - 10 3.16e+0 — 1.59e+0 0.9933 7.92e-1 1.0041 3.96e-1 1.0014 1.98e-1 1.0004 2 - 12 3.16e+0 — 1.59e+0 0.9933 7.92e-1 1.0041 3.96e-1 1.0015 1.98e-1 1.0004 0 3.21e+0 — 1.59e+0 1.0092 7.93e-1 1.0076 3.96e-1 1.0020 1.98e-1 1.0005Table 4
Numerical error ofPspL2 for Example 1.
t h 2 - 2 2 - 3 2 - 4 2 - 5 2 - 6 error rate error rate error rate error rate error rate 2 - 2 4.92e-2 — 1.79e-2 1.4621 4.76e-3 1.9068 1.21e-3 1.9813 3.02e-4 1.9962 2 - 4 2.16e-2 — 7.14e-3 1.5973 1.92e-3 1.8963 4.89e-4 1.9727 1.23e-4 1.9931 2 - 6 2.05e-2 — 6.56e-3 1.6422 1.77e-3 1.8935 4.54e-4 1.9587 1.15e-4 1.9856 2 - 8 2.04e-2 — 6.49e-3 1.6510 1.74e-3 1.9033 4.45e-4 1.9633 1.13e-4 1.9829 2 - 10 2.04e-2 — 6.49e-3 1.6516 1.73e-3 1.9044 4.44e-4 1.9652 1.12e-4 1.9853 2 - 12 2.04e-2 — 6.49e-3 1.6517 1.73e-3 1.9045 4.44e-4 1.9653 1.12e-4 1.9985 0 2.57e-2 — 7.04e-3 1.8699 1.80e-3 1.9695 4.51e-4 1.9926 1.13e-4 1.9982Figure 1
Convergence of Example1 when t = 0 and t = 0.25.
(a)
t = 0 (b)
t = 0.25Figure 2
The true and computed pressure solutions of Example1.
(a)
True solution (b)
Numerical solution when t = 0 (c)
Numerical solution when t = 0.25Example 2.
The solution is dependent on perturbation parametert [taking a case of [2] as example], and the analytical solution of the problem is as follows: (42)u=tcurle-xy/t,p=te-x/t.
Then, the error estimates are listed in Tables5–8, while the convergence rates are plotted with respect to h in Figure 3. It is shown that the convergence rate for u is O(h1/2) which is coincident with our theoretical analysis of Theorem 5.Table 5
Numerical error ofuL2 for Example 2.
t h 2 - 2 2 - 3 2 - 4 2 - 5 2 - 6 error rate error rate error rate error rate error rate 2 - 2 1.63e-1 — 5.13e-3 1.6777 1.42e-3 1.8477 2.07e-4 1.9282 5.47e-5 1.9656 2 - 4 4.22e-2 — 1.81e-2 1.2193 6.77e-3 1.4218 2.09e-3 1.6923 5.82e-4 1.8462 2 - 6 1.17e-1 — 4.74e-2 1.3052 1.59e-2 1.5723 6.37e-3 1.3227 2.43e-3 1.3928 2 - 8 1.68e-1 — 1.05e-2 0.6784 5.41e-2 0.9530 2.13e-3 1.3458 6.46e-3 1.7204 2 - 10 1.71e-1 — 1.17e-1 0.5410 8.00e-2 0.5547 5.05e-2 0.6650 2.62e-2 0.9479 2 - 12 1.71e-1 — 1.17e-1 0.5378 8.17e-2 0.5288 5.69e-2 0.5212 3.92e-2 0.5373Table 6
Numerical error ofpL2 for Example 2.
t h 2 - 2 2 - 3 2 - 4 2 - 5 2 - 6 error rate error rate error rate error rate error rate 2 - 2 2.46e-2 — 1.28e-2 0.9439 6.39e-3 0.9972 3.19e-3 1.0021 1.59e-3 1.0010 2 - 4 1.18e-2 — 6.79e-3 0.7994 3.41e-3 0.9935 1.65e-3 1.0529 8.06e-4 1.0302 2 - 6 1.57e-2 — 6.65e-3 1.2433 2.74e-2 1.2787 1.24e-3 1.1456 5.23e-4 1.2441 2 - 8 1.76e-2 — 8.96e-3 0.9766 4.08e-3 1.1342 1.61e-3 1.3455 6.39e-3 1.3293 2 - 10 1.78e-2 — 9.29e-3 0.9373 4.60e-3 1.0130 2.22e-3 1.0534 1.00e-3 1.1452 2 - 12 1.78e-2 — 9.31e-1 0.9352 4.65e-3 1.0022 2.30e-3 1.0130 1.14e-3 1.0196Table 7
Numerical error ofua for Example 2.
t h 2 - 2 2 - 3 2 - 4 2 - 5 2 - 6 error rate error rate error rate error rate error rate 2 - 2 5.97e-2 — 3.24e-2 0.8801 1.69e-2 0.9368 8.68e-3 0.9655 4.39e-3 0.9816 2 - 4 6.70e-2 — 3.89e-2 0.7846 2.30e-2 0.7567 1.29e-2 0.8299 6.93e-3 0.9015 2 - 6 1.21e-1 — 5.81e-2 1.0551 3.00e-2 0.9542 1.79e-2 0.7418 1.10e-2 0.7019 2 - 8 1.68e-1 — 1.05e-1 0.6733 5.61e-2 0.9077 2.72e-2 1.0464 1.43e-2 0.9267 2 - 10 1.71e-1 — 1.18e-3 0.5405 8.01e-2 0.5528 5.08e-3 0.6591 2.72e-2 0.8991 2 - 12 1.71e-2 — 1.18e-1 0.5378 8.17e-2 0.5286 5.69e-2 0.5206 3.93e-2 0.5353Table 8
Numerical error ofPspL2 for Example 2.
t h 2 - 2 2 - 3 2 - 4 2 - 5 2 - 6 error rate error rate error rate error rate error rate 2 - 2 5.78e-3 — 2.38e-3 1.2802 7.44e-4 1.6785 2.07e-4 1.8445 5.47e-5 1.9222 2 - 4 1.02e-2 — 4.60e-3 1.1446 1.66e-3 1.4703 5.00e-4 1.7296 1.37e-4 1.8671 2 - 6 1.57e-2 — 6.63e-3 1.2474 2.64e-3 1.3316 1.07e-3 1.3008 3.67e-4 1.5438 2 - 8 1.76e-2 — 8.96e-3 0.9765 4.08e-3 1.1342 1.61e-2 1.3466 6.32e-4 1.3442 2 - 10 1.78e-2 — 9.29e-3 0.9373 4.60e-3 1.0130 2.21e-3 1.0534 1.00e-3 1.1452 2 - 12 1.78e-2 — 9.31e-3 0.9352 4.65e-3 1.0022 2.30e-3 1.0130 1.14e-3 1.0195Figure 3
Convergence of Example2 when t=2-2 and t=2-12.
(a)
t = 2 - 2 (b)
t = 2 - 12In the following, we give two tests for the reliability of the proposed method on several high contrast of permeability. The corresponding boundary conditions can be set asu=(1,0) and f=0.Example 3.
We solve the Brinkman equations on a100×100 fine mesh with a high-contrast permeability. In Figure 4(a), the red color denotes regions of media with high permeability K-1=1, and the blue region is low permeability K-1=105. The computed first and second components of the velocity are shown in Figures 4(b) and 4(c). The results look rather similar to [11].Figure 4
(a) Profile of permeability; (b) velocity componentu1; and (c) velocity component u2.
(a) (b) (c)Example 4.
Using Darcy and Brinkman models, we test the effects of numerical modeling on fluid flow in porous media with different sets of permeability distribution. A description of the open foam and fibrous can be shown in Figures5(a) and 5(b) on a 100×100 fine mesh. The computed velocity component u1 of Darcy and Brinkman models is presented in Figure 6. We can find that the diffusion effect of Brinkman model is apparent, while the difference in the upper horizontal boundary is visible as no vertical flow is imposed. These results demonstrate that the numerical method is feasible and the simulation is reasonable [13].Figure 5
Two profiles of permeability: (a) open foam; (b) fibrous.
(a) (b)Figure 6
Velocity componentu1 for two sets of permeability distribution.
(a)
Open form, Darcy (b)
Open form, Brinkman (c)
Fibrous, Darcy (d)
Fibrous, Brinkman



---

*Source: 7625201-2019-06-23.xml*