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.
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., , 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 symbol denoting generic positive constants, independent of any mesh parameters, may take different values at different occurrences. denotes generic small positive constant.
2. The Model Problem and Corresponding Weak Formulation
In the domain , we denote the velocity by and pressure by . The Brinkman model is described by the parametric equationswhere is the media permeability and and are given data. For simplicity, we consider the homogenous boundary conditions: for , we let and in the limit case of ; here 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 be the space of square integrable function in with inner product . We also use the Hilbert spacethen is defined as the subspace of with zero trace on and is the subspace of with zero normal trace on . is the subspace of consisting of functions with zero mean value on .
If , we have and ; in the singular case that the perturbation parameter , the vector space is replaced by .
Define two bilinear forms by The variational formulation of the model is represented by the following: find , such thatThe associated norm on is given by For , we find that . The L-B-B condition holds, as
Theorem 1. There exists a unique solution to problem (5).
Proof. The coercivity and continuity of can be obtained by simple calculation, which meanswhere is a closed subspace of via 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 and are initially proposed by Han and Wu to treat with Stokes problem [15]. We find it is stable uniformly in the parameter . In brief, we firstly recall that is quasiuniformly rectangular partitioning of the domain ; and are staggered rectangular meshes based on . Then, the discrete spaces are defined asLet Using the subspaces and instead of and in the variational problem (5), we obtain the discrete problem: find , such that
3. Convergence Analysis and Error Estimate
In this section, we give the convergence analysis and error estimate. For , the definition of the interpolating by is expressed by the formulation [15]where is the set that contains all vertical edges of the partition , while contains all horizontal edges. For , the interpolation operator can be redefined as [18]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 projection operator.
Referring to [18], we recall some properties of the interpolating and projection .
Lemma 2. (i) There exist two constants and independent of h, such that(ii) For , there exists a constant independent of h, such that(iii) For any , we have
Theorem 3. The discrete inf-sup condition is valid; namely, there is a constant , such that
Proof. Based on the fact that , we find which means
Theorem 4. If satisfy (5) and is the solution of (11), then the following error estimate holds:
Proof. For the error, we subtract (11) from (5) by and obtain thatTaking , we havewhich is derived by combining with the equality (17) and (23). For convenience, we introduce a norm for the norm on , which is commonly used in the following: Then, from (24) and inequality, we further have where are arbitrary positive constants and combine with the quality of the interpolation , which is based on the fact that as on each . Now, we can obtain thatFor , there exist a and a constant independent of [17], satisfying Then, we havewhich meansBased on (16), we get that
In general, the solution depends on the parameter ; the conclusion above will deteriorate as approaches zero. As the introduction in [2], we recall some notation in the following: withwhere denote the vertices of the domain .
Theorem 5. When the solution 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. In the case , system (5) reduces to elliptical about . Referring to [2], we derive thatwhile the solution in case of has the new regularityApplying with (37), we conclude that Combining with (27), (36), and (37), we have Finally, from (29) and (39), we have
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: . 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 to zero. Lastly, we introduce several permeability cases to test the performance of the proposed method.
We need to compute the following useful errors: , , , and corresponding convergence rates.
Example 1. The analytical solution of the problem is as follows: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.

(a) t = 0

(b) t = 0.25

(a) True solution

(b) Numerical solution when t = 0

(c) Numerical solution when t = 0.25
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: 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 is which is coincident with our theoretical analysis of Theorem 5.

(a)

(b)
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 and .
Example 3. We solve the Brinkman equations on a fine mesh with a high-contrast permeability. In Figure 4(a), the red color denotes regions of media with high permeability , and the blue region is low permeability . 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].

(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 Figures 5(a) and 5(b) on a fine mesh. The computed velocity component 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].

(a)

(b)

(a) Open form, Darcy

(b) Open form, Brinkman

(c) Fibrous, Darcy

(d) Fibrous, Brinkman
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.
Conflicts of Interest
The authors declare that they have no conflicts of interest regarding the publication of this paper.
Acknowledgments
This work is supported by the National Natural Science Foundation of China Grants 61503227 and 11801293 and the Shandong Provincial Natural Science Foundation (ZR2018PA003). The authors would like to thank the editor and referees for their valuable comments and suggestions that led to a considerably improved paper.