Numerical Computation of Brinkman Flow with Stable Mixed Element Method

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.


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 parameterdependent 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., 2 − 0 , 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.

Mathematical Problems in Engineering
Throughout this paper, the symbol , denoting generic positive constants, independent of any mesh parameters, may take different values at different occurrences. denotes generic small positive constant.

The Model Problem and Corresponding Weak Formulation
In the domain Ω ⊂ 2 , we denote the velocity by u and pressure by . The Brinkman model is described by the parametric equations where is the media permeability and f and are given data. For simplicity, we consider the homogenous boundary conditions: for > 0, we let u = 0 and u ⋅ n = 0 in the limit case of = 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 spaces (Ω) with the usual norm and seminorm ‖ ⋅ ‖ and | ⋅ | . Let 2 (Ω) be the space of square integrable function in Ω with inner product (⋅, ⋅). We also use the Hilbert space then 1 0 (Ω) is defined as the subspace of 1 (Ω) with zero trace on Ω and 0 (div, Ω) is the subspace of (div, Ω) with zero normal trace on Ω. 2 0 (Ω) is the subspace of 2 (Ω) consisting of functions with zero mean value on Ω.
Define two bilinear forms by (u, k) = 2 (∇u, ∇k) + ( −1 u, k) , The variational formulation of the model is represented by the following: find (u, ) ∈ V × , such that The associated norm on V is given by For 0 < < , we find that ‖u‖ V ≤ 1 ‖u‖ 1 . The L-B-B condition holds, as Theorem 1. There exists a unique solution (u, ) ∈ V × to problem (5).
Proof. The coercivity and continuity of (u, k) can be obtained by simple calculation, which means where V 0 is a closed subspace of V via V 0 = {k ∈ V | (k, ) = 0, ∀ ∈ }. Combining with the inf-sup condition (7), the Brezzi condition is satisfied. It is clear that problem (5) is wellposed and the solution is unique [7,17].
Next, we discuss the finite element discretization of the variational formulation (5). The finite element spaces V h and ℎ are initially proposed by Han and Wu to treat with Stokes problem [15]. We find it is stable uniformly in the parameter ∈ (0, ]. In brief, we firstly recall that ℎ is quasiuniformly rectangular partitioning of the domain Ω; 1 ℎ and 2 ℎ are staggered rectangular meshes based on ℎ . Then, the discrete spaces are defined as Let Using the subspaces ℎ and ℎ instead of and in the variational problem (5), we obtain the discrete problem: find

Convergence Analysis and Error Estimate
In this section, we give the convergence analysis and error estimate. For > 0, the definition of the interpolating Π : where is the set that contains all vertical edges of the partition ℎ , while contains all horizontal edges. For = 0, the interpolation operator can be redefined as [18] ∫ Π 1 where is a set of edges of elements gotten by bisecting the most bottom element edges and most top element edges of and are gotten by bisecting the most left element edges and most right element edges of . Furthermore, we denote : → ℎ as the 2 projection operator.
Referring to [18], we recall some properties of the interpolating Π and projection .
Based on (16), we get that − ℎ ≤ ℎ ( ‖u‖ 2 + ℎ ‖u‖ 2 + 1 ) . In general, the solution (u, ) depends on the parameter ; the conclusion above will deteriorate as approaches zero. As the introduction in [2], we recall some notation in the following: where 1 , . . . , denote the vertices of the domain Ω.
Theorem 5. When the solution (u, ) of system (1)-(2) depends on , then there is a constant independent of and ℎ, such that Proof. Now, we give some new regularity about the solution.

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 ; then, we observe the convergence of the numerical error as the parameter decreases 1 to zero. Lastly, we introduce several permeability cases to test the performance of the proposed method.
Then, the error estimates are listed in Tables 1-4, while the convergence rates are plotted with respect to ℎ 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.

Example 2.
The solution is dependent on perturbation parameter [taking a case of [2] as example], and the analytical solution of the problem is as follows: u = curl − / , Then, the error estimates are listed in Tables 5-8, while the convergence rates are plotted with respect to ℎ in Figure 3. It    is shown that the convergence rate for u is (ℎ 1/2 ) which is coincident with our theoretical analysis of Theorem 5.
In 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 as u = (1, 0) and f = 0.
Example 3. We solve the Brinkman equations on a 100 × 100 fine mesh with a high-contrast permeability. In Figure 4(a), the red color denotes regions of media with high permeability −1 = 1, and the blue region is low permeability −1 = 10 5 . 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 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].

Data Availability
The [numerical results] data used to support the findings of this study are included within the article. These data [numeric results] were obtained by the authors by writing MATLAB codes in order to verify the correctness of the theoretical analysis. There is no data cited.
Mathematical Problems in Engineering

Conflicts of Interest
The authors declare that they have no conflicts of interest regarding the publication of this paper.