Numerical analysis of modular grad-div stability methods for the time-dependent Navier-Stokes/Darcy model

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 grad-div 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 grad-div 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 ${\bf u}\cdot{\bf 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 $\Omega$ $\in$ $R^d(d = 2\ {\rm or} \ 3)$, which is decomposed into a fluid flow region $\Omega_f$ and porous media flow region $\Omega_p$, see figure 1. Here, $\Omega_f\cap\Omega_p = \varnothing$, $\overline{\Omega}_f\cup\overline{\Omega}_p = \overline{\Omega}$ and $\overline{\Omega}_f\cap\overline{\Omega}_p = \Gamma$, $\Gamma_f = \partial {\Omega}_f\cap\partial{\Omega}$, $\Gamma_p = \partial{\Omega}_p\cap\partial{\Omega}$. In the sake of simplicity, we assume that $\partial{\Omega}_f$ and $\partial{\Omega}_p$ are smooth enough throughout this paper.

The time-dependent Navier-Stokes equation govern the fluid flow in $\Omega_f$ is expressed as

here ${\bf u} = {\bf u}(x, t)$ is the fluid velocity filed, $p = p(x, t)$ is the pressure, function ${\bf f}_1$ is the external force, and coefficient $\nu>0$ is the kinetic viscosity.

The Darcy equation govern the porous media flow in $\Omega_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, ${\bf u}_p = {\bf u}_p(x, t)$ the velocity, $\phi = \phi(x, t)$ the hydraulic head, $\mathbf{K}$ is the hydraulic conductivity tensor. We assume that $\mathbf{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, $\rho$ 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 $\Gamma$ by:

where, ${\bf n}_f$ and ${\bf n}_p$ are the unit outward normal vectors on $\partial \Omega_f$ and $\partial\Omega_p$, respectively, and $\tau_i, i = 1, \cdot\cdot\cdot, d-1$, are the orthonormal tangential unit vectors on the interface $\Gamma$, $g$ is the gravitational acceleration, $\alpha$ 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 $\Gamma_f$, $\Gamma_p$:

In this paper, we assume that:

The assumption is not hold for general case of Navier-Stokes/Darcy Model. But for the gentle river, the water infiltration satisfies the assumption ${\bf u}\cdot {\bf n}_f > 0$ on the interface. It is a special case of Navier-Stokes/Darcy Model.

Next, Hillbert spaces will be introduced:

where $(\cdot, \cdot)_D$ denotes the $L^2$ inner produce in the domain D, with corresponding norm $||\cdot ||_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 $\nabla\cdot {\bf u} = 0$, we define a trilinear form ${a_{f, c}}(\cdot; \cdot, \cdot)$ as follows

under the condition (8), we have

Then, we have the following estimates for $a_{f, c}$:

In addition, we recall the Poincar$\acute{e}$ inequality and trace inequality. There exist positive constants $c_p$ and $c_t$ which only depend on the domain $\Omega_f$ and exist $\tilde{c}_p$ and $\tilde{c}_t$ which only depend on the domain $\Omega_p$.

Thus, the weak formulation of the time-dependent Naiver-Stokes/Darcy model is to find ${\bf u}:[0, T]\to H_f, p:[0, T]\to Q, \phi :[0, T]\to H_p,$ such that

Where

Bilinear forms $a_f$ and $a_p$ are continuous and coercive:

The interface coupling term $a_{\Gamma}$ satisfies the following estimates:

Lemma 2.1. We assume that

and $\mathbf{K}$ is uniformly bounded and positive definite in $\Omega_p$, there exist two constants $k_{min}>0, \ k_{max}>0$ such that

Furthermore, ${\bf u}\in L^2(\Omega_p)^d, \phi_0\in L^2(\Omega_p)$, therefore any solution $({\bf u}, p, \phi)\in {\Big(}L^2(0, T;H_f)\cap H^1(0, T;L^2(\Omega_f)^d){\Big)}\times L^2(0, T;Q)\times L^2(0, T;H_p)$ of (1)-(7) is also the solution to the equation (12). The converse of the statement is also true.

Lemma 2.2. (Discrete Gronwall Lemma). Let $\Delta t$, H, $a_n$, $b_n$, $c_n$, $d_n$ be nonnegative mumbers for $n\geq0$ such that for $N\geq1$. If

then for all $\Delta t>0$,

3.   The modular grad-div stabilization algorithms

We construct $\tau_h$ is a quasiuniform triangulation of the domain $\Omega_f\bigcup\Omega_p$, depending on a positive parameter $h>0$, which is made up of triangles if $d = 2$ or tetrahedra if $d = 3$. Then we define the finite element subspace of $H_f$, $Q$, $H_p$ as $H_{fh}$, $Q_h$, $H_{ph}$. We assume that the space pairs $(H_{fh}, Q_h)$ satisfies the discrete LBB condition: there exists a positive constant $\beta$ independent of h, such that $\forall q_h\in Q_h$, $\exists {\bf v}_h\in H_{fh}$, ${\bf v}_{fh}\ne 0$

Next, we divide time interval [0, T] into: $0 = t_0 < t_1 < \cdots < t_N = T$ with $\Delta t = t_i-t_{i-1}$. This leads to $t_m = m\Delta t$ and $T = N\Delta t$ as a uniform distribution of discrete time levels. Let $({\bf u}_{h}^{n+1}, p_{h}^{n+1}, \phi_{h}^{n+1})$ denote the discrete approximation of $({{\bf u}_{h}(t_{n+1})}, p_h(t_{n+1}), \phi_h(t_{n+1}))$. The modular grad-div scheme our proposed can be written as:

Algorithm 1 (The modular grad-div scheme).

For $\forall {\bf v}_{h}\in H_{fh}$, $q_{h}\in Q_h$, and $\psi_h\in H_{ph}$.

Step1. Given ${\bf u}_{h}^{n}\in H_{fh}$, $p_{h}^{n}\in Q_h$, and $\phi_{h}^{n}\in H_{ph}$, find ${\hat{\bf u}}_{h}^{n+1}, p_{h}^{n+1}, \phi_{h}^{n+1}$ satisfying:

Step2. Given $\hat{\bf u}_{h}^{n+1}\in H_{fh}$, find ${\bf u}_{h}^{n+1}\in H_{fh}$ satisfying:

Lemma 3.1. For Algorithm 1, we can obtain the following result,

Proof. Refer to Lemma 6 of [3] for proof details.

Theorem 3.2. (Unconditional Stability) For any $N>1$, the solution of the Algorithm 1, satisfy

Where the constant $C = exp(\Delta t\sum\limits_{n = 0}^{N-1} max{\frac {2C_5g}{h\lambda_{min}}, \frac{2C_6g^2}{h\nu}})$.

Proof. Taking ${\bf v}_h = 2\Delta\hat{\bf u}_h^{n+1}$, $q_h = 2\Delta p_h^{n+1}$ and $\psi_h = 2\Delta \phi_h^{n+1}$ in (16), using the property (10), and we can obtain $a_{f, c}({\bf u}_h^n, \hat{\bf u}_h^{n+1}, \hat{\bf u}_h^{n+1})\geq 0$. 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 ${\bf u}^n, p^n$ and $\phi^n$ be the true solution at time $t^n = n\Delta t$. Assuming that the true solution have the following regularities,

Then define a projection operator [12]:

when some regularity conditions on $({\bf u}(t), p(t), \phi (t))$ are statisfied, $(P_h{\bf u}(t), P_h p(t), P_h \phi (t))$ is an approximation of $({\bf u}(t), p(t), \phi(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].

Theorem 4.2. Under the regularity assumption (18). Suppose $\beta>0$, then there exists a constant $C>0$, such that

Proof. The true solution satisfies the following relations:

Subtracting (16) from (21), we arrive that

Setting ${\bf v}_h = 2\Delta t\theta_{\hat{\bf u}}^{n+1}, q_h = 2\Delta t\theta_p^{n+1}$, and $\psi_h = 2\Delta t\theta_{\phi}^{n+1}$ in (22), the resulting equations are added, this yields:

Next, we bound each term on the right hand side of (23), by virtue of Cauchy-Schwarz-Young inequality,

For the interface terms, we treat them as follows:

So, we have the following estimates:

For the trilinear form, we have

then

Combining Lemma 4.1 with the properties (13), we obtain

Inserting the above results into (29) and let $\varepsilon_1 = \varepsilon_2 = \cdots\cdots = \varepsilon_{10} = \nu$, we can get the following estimates:

Denote $C^{\ast} = exp(\Delta t\sum\limits_{n = 0}^{N-1}max(\frac{6C_5g}{h\lambda_{min}}, \frac{10C_2^2}{\nu}||{\bf u}||_{\infty, 2}^2, \beta, \frac{10C_6g^2}{h\nu}))$. Sum (30) from $n = 0$ to $n = N-1$, together with Theorem 3.2 and Lemma 2.2, we get the following inequality:

Finally, we have

where $C>0$ is a constant.

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.

Algorithm 2 (The BDF2 modular grad-div scheme for Stokes/Darcy model).

For $\forall {\bf v}_{h}\in H_{fh}$, $q_{h}\in Q_h$, and $\psi_h\in H_{ph}$.

Step1. Given ${\bf u}_h^{n-1}, {\bf u}_{h}^{n}\in H_{fh}$, $p_{h}^{n}\in Q_h$, and $\phi_{h}^{n-1}, \phi^n_h\in H_{ph}$, find $\hat{\bf u}_{h}^{n+1}, p_{h}^{n+1}, \phi_{h}^{n+1}$ satisfying:

Step2. Given $\hat{\bf u}_{h}^{n+1}\in H_{fh}$, find ${\bf u}_{h}^{n+1}\in H_{fh}$ satisfying:

Algorithm 3 (The BDF2 modular grad-div scheme for Navier-Stokes/Darcy model).

For $\forall {\bf v}_{h}\in H_{fh}$, $q_{h}\in Q_h$, and $\psi_h\in H_{ph}$.

Step1. Given ${\bf u}_{h}^{n-1}, {\bf u}_{h}^{n}\in H_{fh}$, $p_{h}^{n}\in Q_h$, and $\phi_{h}^{n-1}, \phi^n_h\in H_{ph}$, find ${\hat{\bf u}}_{h}^{n+1}, p_{h}^{n+1}, \phi_{h}^{n+1}$ satisfying:

Step2. Given $\hat{\bf u}_{h}^{n+1}\in H_{fh}$, find ${\bf u}_{h}^{n+1}\in H_{fh}$ satisfying:

5.   The numerical results

In this section, We will compare two grad-div schemes with standard scheme respectively to justify the results of the theoretical analysis. We implement numerical experiments using software Freefem++.

The domain $\Omega$ be decomposed into $\Omega_f = (0, 1)\times(1, 2)$ and $\Omega_p = (0, 1)\times(0, 1)$ with the interface $\Gamma = (0, 1)\times\{1\}$. The exact solution is taken as follows:

Here, the parameters $n, \rho, g, \nu, \mathbf K, S_0$ and $\alpha$ are set to 1. The initial conditions, boundary conditions, and the forcing terms follows the exact solution, and we set $h = \Delta 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 ${\bf u}$, pressure $p$ and hydraulic head $\phi$ by using the standard scheme, standard grad-div scheme and modular grad-div scheme, respectively. where $\gamma = 1$ and $\beta = 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 grad-div 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 ${\bf K}$ to show the superiority of modular grad-div scheme. Let's fix $\Delta t = h = \frac{1}{40}$. Table 4 shows the $||\nabla \cdot e_{\bf u}||_f$ for the standard scheme without stabilization, standard grad-div scheme and modular grad-div scheme with hydraulic conductity ${\bf K} = {\bf I}$, $0.1{\bf I}$, $0.01{\bf I}$, $0.001{\bf I}$. We observe that the divergence of velocity error for all three schemes increase with the decrease of ${\bf 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.

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 grad-div scheme.