NUMERICAL ANALYSIS OF MODULAR GRAD-DIV STABILITY METHODS FOR THE TIME-DEPENDENT NAVIER-STOKES/DARCY MODEL

In this paper, we construct a modular grad-div stabilization method for the Navier-Stokes/Darcy model, which is based on the first order Backward Euler scheme. This method does not enlarge the accuracy of numerical solution, but also can improve mass conservation and relax the influence of parameters. Herein, we give stability analysis and error estimations. Finally, by some numerical experiment, the scheme our proposed is shown to be valid.


1.
Introduction. Numerical methods of Navier-Stokes/Darcy have attracted a lot of attention. So far, a great deal of numerical methods are proposed to solve this model by virtue of different ways, such as finite element methods [12], discontinuous Galerkin finite element methods [5], two-grid methods [1,15,16], modified two-grid methods [6], partitioned time stepping method [7], characteristic stabilized finite element methods [8], mortar finite element methods [2], grad-div stabilized projection finite element method [14], modular grad-div method [11] and so on. The grad-div stabilized method is first introduced in [4], which can penalize mass conservation and improve the solution quality efficiently. Recently, the effectiveness of the graddiv stabilized method has been proved in finite element simulation of Stokes and Navier-Stokes equation [10]. However, this method leads to a singular matrix stemming from grad-div term, and the larger stabilized parameter will cause solver breakdown. As a alternative method for grad-div stabilization, the modular graddiv stabilization method is introduced in [3]. The modular grad-div stabilization method for the Stokes/Darcy model is proposed in [13]. The modular grad-div stabilization method is not only easy to implement, but also avoids the influence of large parameters on the solution as well as preserving the advantages of the grad-div stabilization method.
In this paper, we extended the modular grad-div stabilization methods from Stokes/Darcy model to Navier-Stokes/Darcy model. Compare to Navier-Stokes model, the convection term causes some difficulties in theoretical analysis and numerical simulation. To deal with convection term, the assumption u · n f > 0 is used [6,9]. Our proposed method not only relax the influence of parameters but also improves the mass conservation.
The rest of this paper is organized as follows: In section 2, some notations and the time-dependent Navier-Stokes/Darcy model are introduced; In section 3, we give the modular grad-div stabilization method, and stability analysis is provided ; In section 4, under some regularity assumptions imposed on the true solution, error estimates are also given; In section 5, Some numerical experiments are given to verify the theoretical result, we compared with the standard scheme, standard grad-div scheme and modular grad-div scheme in the numerical experiment; Finally some conclusions are obtained in section 6.
2. Functional setting of the time-dependent Navier-Stokes/Darcy model. The model we considered is confined in a bounded domain Ω ∈ R d (d = 2 or 3), which is decomposed into a fluid flow region Ω f and porous media flow region Ω p , see figure 1.
In the sake of simplicity, we assume that ∂Ω f and ∂Ω p are smooth enough throughout this paper. The time-dependent Navier-Stokes equation govern the fluid flow in Ω f is expressed as is the fluid velocity filed, p = p(x, t) is the pressure, function f 1 is the external force, and coefficient ν > 0 is the kinetic viscosity. The Darcy equation govern the porous media flow in Ω p is expressed as : here the first equation is the saturated flow model and the second equation is the Darcy's law. S 0 is the specific mass storativity coefficient, u p = u p (x, t) the velocity, φ = φ(x, t) the hydraulic head, K is the hydraulic conductivity tensor. We assume that K is a positive symmetric tensor and the function f 2 is a source term.
We know that is the piezometric head, where p p denotes the dynamic pressure, ρ is the height from a reference level. Substituting the second formula of (2) into the first formula of (2), the following Darcy equation can be obtained: The interface conditions of the conservation of mass, balance of forces, and the Beavers-Joseph-Saffman condition are imposed on the interface Γ by: where, n f and n p are the unit outward normal vectors on ∂Ω f and ∂Ω p , respectively, and τ i , i = 1, · · ·, d − 1, are the orthonormal tangential unit vectors on the interface Γ, g is the gravitational acceleration, α is a positive parameter depending on the properties of the porous medium and must be determined experimentally.
For simplicity of analysis, we impose the following boundary conditions on Γ f ,Γ p : In this paper, we assume that: u · n f > 0 on Γ.
The assumption is not hold for general case of Navier-Stokes/Darcy Model. But for the gentle river, the water infiltration satisfies the assumption u · n f > 0 on the interface. It is a special case of Navier-Stokes/Darcy Model.
Next, Hillbert spaces will be introduced: where (·, ·) D denotes the L 2 inner produce in the domain D, with corresponding norm || · || D . In the rest of this paper, we neglect the subscript. The spaces H f and H p are equipped with the following norms: Some discrete norms are defined as, Due to ∇ · u = 0, we define a trilinear form a f,c (·; ·, ·) as follows under the condition (8), we have Then, we have the following estimates for a f,c : In addition, we recall the Poincaré inequality and trace inequality. There exist positive constants c p and c t which only depend on the domain Ω f and existc p and c t which only depend on the domain Ω p .
Thus, the weak formulation of the time-dependent Naiver-Stokes/Darcy model is to find u : Bilinear forms a f and a p are continuous and coercive: The interface coupling term a Γ satisfies the following estimates: Lemma 2.1. We assume that and K is uniformly bounded and positive definite in Ω p , there exist two constants k min > 0, k max > 0 such that (7) is also the solution to the equation (12). The converse of the statement is also true.
Lemma 2.2. (Discrete Gronwall Lemma). Let ∆t, H, a n , b n , c n , d n be nonnegative mumbers for n ≥ 0 such that for N ≥ 1. If 3. The modular grad-div stabilization algorithms. We construct τ h is a quasiuniform triangulation of the domain Ω f Ω p , depending on a positive parameter This leads to t m = m∆t and T = N ∆t as a uniform distribution of discrete time levels.
). The modular grad-div scheme our proposed can be written as: Step2 Lemma 3.1. For Algorithm 1, we can obtain the following result, Refer to Lemma 6 of [3] for proof details.
Where the constant C = exp(∆t (16), using the property (10), and we can obtain Then seeing [13] Theorem 3.1 for a detail proof.
4. Error analysis. In this section, we will give some error estimates of our proposed method. Denote u n , p n and φ n be the true solution at time t n = n∆t. Assuming that the true solution have the following regularities, φ ∈ L ∞ (0, T ; Hp ∩ H k+1 (Ωp)), φt ∈ L ∞ (0, T ; H k+1 (Ωp)), φtt ∈ L 2 (0, T ; L 2 (Ωp)) Then define a projection operator [12]: when some regularity conditions on (u(t), p(t), φ(t)) are statisfied, (P h u(t), P h p(t), P h φ(t)) is an approximation of (u(t), p(t), φ(t)), with the following properties: Next, we define the following error equations Lemma 4.1. For Algorithm 1, The following inequality holds Proof. For the detailed proof process, please refer to Lemma 10 in [3].

Remark 1.
Here we can propose a second-order backward differentiation formula(BDF2) methods for Stokes/Darcy model and Navier-Stokes/Darcy model, respectively. We will analyze it in the future.
Here, the parameters n, ρ, g, ν, K, S 0 and α are set to 1. The initial conditions, boundary conditions, and the forcing terms follows the exact solution, and we set h = ∆t in experiments. The famous Taylor-Hood element(P2-P1) for the Navier-Stokes problem and the piecewise quadratic polynomials(P2) for the Darcy flow are used, respectively. The numerical results are presented in Table 1-Table 4. Table 1 , Table 2 and Table 3 shows the error and convergence order of velocity u, pressure p and hydraulic head φ by using the standard scheme, standard graddiv scheme and modular grad-div scheme, respectively. where γ = 1 and β = 0.2 for standard grad-div scheme and modular grad-div scheme. By observation and comparison, it can be found that the divergence velocity errors of the standard graddiv and modular grad-div scheme are smaller than that of the standard scheme, and the modular grad-div scheme is more accurate. Moreover both numerical results is consist with the theoretical analysis.
Next, we set up a numerical experiment by using the different parameter K to show the superiority of modular grad-div scheme. Let's fix ∆t = h = 1 40 . Table 4 shows the ||∇ · e u || f for the standard scheme without stabilization, standard grad-div scheme and modular grad-div scheme with hydraulic conductity K = I,0.1I,0.01I,0.001I. We observe that the divergence of velocity error for all three schemes increase with the decrease of K. Yet their growth is relatively small, in particular there's almost no change in modular grad-div scheme. Such results are also consistent with our theoretical analysis.   Table 4 The ||∇ · e u || f for the standard without grad-div scheme, standard scheme and grad-div scheme with vaying hydraulic conductivity tensor K. K Non-stabilized Standard grad-div modular grad-div I 0.0169173 0.0108085 0.077557 1e − 1I 0.0202974 0.0130437 0.0775562 1e − 2I 0.0464625 0.0301879 0.077554 1e − 3I 0.124238 0.0824988 0.0775521 6. Conclusion. In this paper, we extend the grad-div stabilized method from Stokes/Darcy model to Navier-Stokes/Darcy model. Stability and error estimates are provided. Numerical experiments confirm the theoretical analysis, and show that the modular grad-div scheme is more efficient than that of the standard graddiv scheme.