Convergence analysis of Euler and BDF2 grad-div stabilization methods for the time-dependent penetrative convection model

1.   Introduction

The grad-div stabilization method was first introduced for the finite element approximation of the incompressible Navier-Stokes equations by adding the grad-div stabilization term to the momentum equation in ^[12]. For the most of Lagrange finite element pairs, such as the Taylor-Hood element and MINI element, the grad-div stabilization term is nonzero, and there is no point-wise mass conservation. Thus, this stabilization term is used to improve the mass conservation and to reduce the velocity error caused by the pressure error. The influence of the grad-div stabilization term on the accuracy of the numerical solution was studied in ^[40] for the steady Stokes problem.

It is well known that the standard Galerkin finite element methods are not suitable for the numerical simulation of the incompressible flow problems in the case of the high Reynolds number or the small viscosity. Different stabilized methods have been developed to overcome the instabilities of numerical simulation, such as Galerkin least square methods in ^[13,25], the residual-free bubbles methods in ^[14,15], the large eddy simulation methods in ^[45], the sub-grid scale methods in ^[23,31], the variational multiscale methods in ^{[26,27,34,49]} and the streamline upwind Petrov-Galerkin (SUPG) method ^[5]. We note that the grad-div stabilization method is also a powerful tool for the high Reynolds flow problems. There have an amount of works for this issue. For example, the SUPG and grad-div stabilization methods were studied for the generalized Oseen problem in ^[37], where it was concluded that the SUPG method is less important for the inf-sup stable pair of velocity and pressure due to the choice of stabilization parameter. Moreover, if one considers only the grad-div stabilization method, it was acknowledged that the grad-div method can lead to stabilized results. The grad-div method for the rotation form of the Navier-Stokes equations can be found in ^[32], where numerical results showed that the difference between the skew-symmetric and the rotation form of the nonlinearity is from the increased error in the Bernoulli pressure, which results in the increasing of velocity error. The use of the grad-div stabilization term ameliorates the effect and reduces the velocity error. For the time-dependent Oseen problem, the semi-discrete and fully discrete schemes based on the grad-div stabilization method were analyzed in ^[17]. Under the assumption of sufficiently smooth solutions and based on the specific Stokes projection, the authors proved the optimal error estimates for velocity and pressure with constants that are independent of the viscosity. Subsequently, this observation was extended to the time-dependent Navier-Stokes equations in ^[18], where the constants in the optimal error bounds do not depend on the inverse power of viscosity. Thus, the results in ^[17] and ^[18] gave a theoretical confirmation about stable simulations of the high Reynolds number flows by using the grad-div stabilization method. Recently, the analysis technique in ^[18] was extended to the second-order backward differentiation formula (BDF2) scheme with variable time-step size of the Navier-Stokes equations in ^[21]. A review study about error analysis of the grad-div stabilization methods can be found in ^[20].

Since there are some advantages of the grad-div stabilization method in numerical simulations of the incompressible flows, it has been applied to some related and coupled systems, such as the magnetohydrodynamic (MHD) equations ^[11] and the time-dependent dual-porosity-Stokes equations ^[33]. In this paper, we will consider the following time-dependent penetrative convection equations in the nondimensional form:

where $\Omega$ is a bounded and convex polyhedral domain in $\mathbb R^3$ with boundary $\partial\Omega$, and $(0, T]$ is a finite time interval. In the above system, the unknown $({{\textbf{u}}}, p, \theta)$ denotes the velocity of the fluid, the pressure and temperature, respectively. $\nu > 0$ and $\kappa > 0$ represent the viscosity coefficient and thermal conductivity parameter. ${{\textbf{f}}}$ and $g$ are two given functions. ${\textbf{i}}_{3}$ is the unit basis vector given by ${\textbf{i}}_{3} = (0, 0, 1)^{\mathrm{T}}$, and $(\gamma_{1}\theta+\gamma_{2}\theta^{2}){\textbf{i}}_{3}$ in (1.1) denotes the buoyancy term. The constraint $\nabla\cdot{{\textbf{u}}} = 0$ represents the incompressibility of the fluid.

The penetrative convection equations (1.1)–(1.3) used to describe the motion in which convection in a thermally unstable region extends into an adjacent stable region. It has a wide range of applications in the fields of atmospheric dynamics, atmospheric fronts, katabatic winds, dense gas dispersion, natural ventilation, cooling of electronic equipment ^{[29,36,39,42]}. The stabilities of penetrative convection equations have been studied in ^{[7,8,16,41,46]}, and numerical simulations have been reported in ^[3,38]. Convergence and error estimates of finite element fully discrete schemes were studied in ^[44] and ^[6], where the BDF2 and Crank-Nicolson schemes were used, respectively. However, the high Reynolds number was not considered in ^[6,44]. By considering a linear Boussinesq approximation, i.e., $\gamma_2 = 0$ in (1.1), a grad-div stabilization finite element method was studied in ^[19]. In addition, a rigorous error analysis was given and error estimates were derived in ^[19], where the constants in error bounds are independent of inverse powers of the viscosity and thermal conductivity coefficients.

In this paper, based on the grad-div stabilization method, we consider and study two finite element schemes by using the first-order backward Euler and BDF2 formulas to discrete the time derivatives, respectively. The proposed schemes are both linearized semi-implicit schemes where we use the implicit-explicit method and the standard extrapolation formula to deal with nonlinear terms. Moreover, these schemes are unconditionally stable without any condition of the time step size and mesh size. Following the analysis techniques developed in ^[18,19], we derive error bounds of the velocity and temperature in which the constants are independent of negative powers of the viscosity and thermal conductivity coefficients.

The rest of this paper is organized as follows. In next section, we state some notations and recall some known results in finite element theory. In sections three and four, the first-order Euler and second-order BDF2 schemes with the grad-div stabilization are proposed, respectively. Unconditional stability of numerical schemes and error analysis are also presented in these sections. In section five, we give numerical results to support the theory analysis and numerical stability of the grad-div stabilization method for the high Reynolds flows.

2.   Materials and methods

For $k\in\mathbb N^+$ and $1\leq p\leq +\infty$, let $L^{p}(\Omega)$ and $W^{k, p}(\Omega)$ denote the standard Lebesgue space and Sobolev space, respectively. When $p = 2$, $W^{k, 2}(\Omega)$ is the Hilber space $H^k(\Omega)$. The norms in $L^{p}(\Omega)$, $W^{k, p}(\Omega)$ and $H^k(\Omega)$ are defined as the classical senses (cf.^[1]) and are denoted by $\|\cdot\|_{L^{p}}$, $\|\cdot\|_{W^{k, p}}$ and $\|\cdot\|_{H^{k}}$. We define $H_0^{k}(\Omega)$ to be the subspace of $H^{k}(\Omega)$ of functions with zero trace on $\partial\Omega$. The dual space of $H_0^1(\Omega)$ is denoted by $H^{-1}(\Omega)$. The boldface notations ${{\textbf{L}}}^p(\Omega)$, ${{\textbf{W}}}^{k, p}(\Omega)$ and ${{\textbf{H}}}^k(\Omega)$ are used to denote the vector-value spaces $L^p(\Omega)^3$, $W^{k, p}(\Omega)^3$ and $H^k(\Omega)^3$, respectively. The corresponding norms are denoted by $\|\cdot\|_{{{\textbf{L}}}^{p}}$, $\|\cdot\|_{{{\textbf{W}}}^{k, p}}$ and $\|\cdot\|_{{{\textbf{H}}}^{k}}$.

For simplicity, we denote ${{\textbf{V}}} = {{\textbf{H}}}_{0}^{1}(\Omega)$, $Y = H_{0}^{1}(\Omega)$,

Introduce the trilinear forms $c_1(\cdot, \cdot, \cdot)$ and $c_2(\cdot, \cdot, \cdot)$, which are given by

for ${{\textbf{u}}}_{i}\in{{\textbf{V}}}$ with $i = 1, 2, 3$ and $\theta_{j}\in Y$ with $j = 1, 2$. By the integration by parts, it is easy to check that the trilinear forms $c_1(\cdot, \cdot, \cdot)$ and $c_2(\cdot, \cdot, \cdot)$ satisfy the skew-symmetric properties, i.e.,

which give

For given ${{\textbf{f}}}\in {{\textbf{L}}}^2(\Omega), g\in L^2(\Omega)$, ${{\textbf{u}}}_0\in{{\textbf{L}}}^2(\Omega)$ and $\theta_0\in L^2(\Omega)$, in terms of ther skew-symmetric properties (2.1), there {holds} the following energy inequalities:

and

for $0\leq t\leq T$.

Let $0 = t_0 < t_1 < \cdots < t_N = T$ be a uniform partition of the time interval $[0, T]$ with time step $\tau = T/N$ and $t_n = n\tau$ with $0\leq n\leq N$. Denote ${{\textbf{u}}}^n = {{\textbf{u}}}(t_n)$, $p^n = p(t_n)$, $\theta^n = \theta(t_n)$, ${{\textbf{f}}}^{n} = {{\textbf{f}}}(t_{n})$ and $g^{n} = g(t_{n})$. For any sequence of functions $\{f^n\}_{n = 0}^N$, denote the backward Euler discretization and BDF2 discretization by

and

where $\widehat {f}^{n}$ is the standard extrapolation formula. For the discrete time derivative $D_{2, \tau}$, there holds the telescope formula ^[35]:

We define the finite element spaces as follows. Let $\mathcal {T}_h = \{K_j\}_{j = 1}^L$ be a quasi-uniform tetrahedron partition of $\Omega$ with the mesh size $h = \max \, \{diam \, K_j, \ j = 1, \cdots, L\}$. In finite element discretization of (1.1)–(1.3), we use the Taylor-Hood ($P_2-P_1$) finite element for the velocity ${{\textbf{u}}}$ and pressure $p$, and the piece-wise quadratic finite element for the temperature $\theta$. The corresponding finite element subspaces of ${{\textbf{V}}}, Q$ and $Y$ are denoted by ${{\textbf{V}}}_h,$ $Q_h$ and $Y_h$, respectively, i.e.,

It is {well known} that the discrete inf-sup condition holds for the inf-sup stable element such that

where $\beta_{0} > 0$ is independent of the mesh size $h$.

Denote the space of discrete divergence-free functions by

In the sequel, let $\pi_h$ and $R_h$ be the $L^2$ orthogonal and Ritz orthogonal projections onto $Q_h$ and $Y_h$, respectively. Then the following error bounds and stability hold ^[4,43]:

Let $I_{h}$ be the Lagrange interpolation operator onto ${{\textbf{V}}}_h$. The following bound can be found in ^[4]:

where $n > 3/p$ when $1 < p\leq\infty$ and $n\geq 3$ when $p = 1$.

Next, we recall some known results about finite element approximations of Stokes problem in ^[17,18]. Consider the Stokes problem:

Let us denote $({{\textbf{u}}}_h, p_h)\in{{\textbf{V}}}_h\times Q_h$, the mixed finite element approximation solution to the Stokes problem (2.10). Then one has error estimates ^[22,30]

where $C > 0$ is independent of $h, \nu$ and $\kappa$. It is clear that the error estimates of the velocity depend on the negative power of $\nu$. To prove error bounds with constants being independent of the negative power of $\nu$, we recall the specific projection of $({{\textbf{u}}}, p)$ introduced in ^[17] that we will denote by ${{\textbf{s}}}_h$ defined by

Let $({{\textbf{u}}}, p, \theta)$ be the solution to the system (1.1)–(1.5) with ${{\textbf{u}}}\in {{\textbf{L}}}^\infty((0, T];{{\textbf{V}}}\cap{{\textbf{H}}}^3(\Omega))$, $p\in L^\infty((0, T];Q\cap H^2(\Omega))$, $\theta\in L^\infty((0, T];Y\cap H^3(\Omega))$ and ${{\textbf{u}}}_t\in {{\textbf{L}}}^\infty((0, T];{{\textbf{H}}}^1(\Omega))$. We note that $({{\textbf{u}}}, 0)$ is the solution to the Stokes problem (2.10) with the right-hand side ${{\textbf{F}}}$ given by

Let $({{\textbf{s}}}_h, l_h)\in {{\textbf{V}}}_h\times Q_h$ be the corresponding mixed finite element approximation of $({{\textbf{u}}}, 0)$. Then, from (2.11)–(2.13), we have

where $C > 0$ is independent of $h$ and $\nu$. In terms of the inverse inequality in the theory of the finite element method and $\|I_{h}{{\textbf{u}}}\|_{{{\textbf{L}}}^\infty}\leq C\|{{\textbf{u}}}\|_{{{\textbf{L}}}^\infty}$, we have

where $C > 0$ is independent of $h, \nu$ and $\kappa$. Furthermore, there holds

where $C > 0$ is independent of $h, \nu$ and $\kappa$.

Finally, we recall the discrete Gronwall inequality established in ^[24].

Lemma 2.1. Let $a_k, b_k, c_k$ and $\gamma_k$, for integers $k\geq 0$, be the nonnegative numbers such that

Suppose that $\tau\gamma_k < 1$ for all $k$ and set $\sigma_k = (1-\tau \gamma_k)^{-1}$, then

Remark 2.1. If the first sum on the right in (2.19) extends only up to $n-1$, then the estimate (2.20) holds for all $\tau > 0$ with $\sigma_k = 1$.

3.   The Euler grad-div stabilization finite element approximation

In this section, we consider the first-order Euler fully discrete scheme for the problems (1.1)–(1.5). This Euler scheme is a semi-implicit scheme, i.e., one only solves two linear systems at each time discrete level. Based on the grad-div stabilization method, we propose the following first-order Euler finite element scheme for $0\leq n\leq N-1$.

Step Ⅰ: Find $\theta_h^{n+1}\in Y_h$ by

with the initial iteration values $\theta_{h}^{0} = R_h\theta_{0}$, ${{\textbf{u}}}_{h}^{0} = I_h{{\textbf{u}}}_{0}$ and for any $\psi_h\in Y_h$.

Step Ⅱ: Find $({{\textbf{u}}}_{h}^{n+1}, p_{h}^{n+1})\in {{\textbf{V}}}_{h}\times Q_{h}$ by

for any $({{\textbf{v}}}_{h}, q_{h})\in {{\textbf{V}}}_{h}\times Q_{h}$, where $\beta > 0$ is the stabilization parameter. The detailed discussion about the choice of $\beta$ can be found in ^[28].

Next lemma gives the stability of the scheme (3.1)–(3.2), where the discrete energy inequalities can be viewed as the discrete version of the energy estimates (2.2)–(2.3). Moreover, since (3.1)–(3.2) are both linearized problems, the discrete energy inequalities imply the existence and uniqueness of numerical solution $({{\textbf{u}}}_h^{n+1}, p_h^{n+1}, \theta_h^{n+1})$ to the Euler scheme (3.1)–(3.2).

Lemma 3.1. For $0\leq n\leq N-1$ and all $\tau > 0$, $h > 0$, the finite element scheme (3.1) and (3.2) has a unique solution $\theta_{h}^{n+1}\in Y_{h}$ and $({{\textbf{u}}}_{h}^{n+1}, p_{h}^{n+1})\in {{\textbf{V}}}_{h}\times Q_{h}$. Moreover, the discrete energy inequalities hold:

and

for all $0\leq m\leq N-1$, where $C > 0$ is independent of $\tau$ and $h$.

Proof. Setting $\psi_{h} = 2\tau\theta_{h}^{n+1}$ in (3.1), using the skew-symmetric property (2.1), the Hölder inequality and Young inequality, we have

Summing up the above inequality from $n = 0$ to $n = m$ gives (3.3).

Similarly, taking $({{\textbf{v}}}_{h}, q_h) = 2\tau({{\textbf{u}}}_{h}^{n+1}, p_h^{n+1})$ in (3.2), we have

Summing up the above inequality from $n = 0$ to $n = m$ and noticing (3.3), we get (3.4).□

Denote

The ${{\textbf{s}}}_h(t)$ is the finite element solution to the Stokes problem (2.10) with the right-hand side (2.14) and ${{\textbf{s}}}_{h}^{n} = {{\textbf{s}}}_{h}(t_n)$.

Thanks to the discrete inf-sup condition (2.5), the finite element scheme (3.2) is equivalent to: Find ${{\textbf{u}}}_{h}^{n+1}\in {{\textbf{V}}}_{0h}$ by

For simplicity, we take ${{\textbf{s}}}_h^{0} = {{\textbf{u}}}_{h}^{0}$. For $0\leq n\leq N-1$, ${{\textbf{s}}}_h^{n+1}$ and $\theta^{n+1}$ satisfy

and

In (3.6), we use $\nabla\cdot{{\textbf{u}}}^{n+1} = 0$ in $\Omega$ and $(\nabla\cdot{{\textbf{v}}}_h, \pi_h p^{n+1}) = 0$ due to ${{\textbf{v}}}_h\in {{\textbf{V}}}_{0h}$.

Next, we give error equations. Subtracting (3.5) from (3.6) leads to

where

Subtracting (3.1) from (3.7), we have

where

The main result in this section is presented in the following theorem.

Theorem 3.1. Let $({{\textbf{u}}}, p, \theta)$ be the unique solution to the time-dependent penetrative convection problems (1.1)–(1.5) and satisfy the following regularity assumptions:

Then, for the first-order Euler grad-div scheme (3.1)–(3.2) under the time step condition $\tau\leq Ch^2$, when $h$ and $\tau$ are sufficiently small, we have the following error estimate:

with $0\leq m\leq N-1$, where $C_0 > 0$ is independent of $h, \tau$, $\nu$ and $\kappa$.

Proof. Taking ${{\textbf{v}}}_h = 2\tau e_h^{n+1}$ in (3.8) and summing up the resulting equation from $n = 0$ to $n = m$, we have

We estimate the right-hand side of (3.12) as follows. For $J_1^{n+1}$, we split it as

According to the Taylor's formula and the regularity assumption (3.10), one has

Notice that

From the Hölder inequality, we have

Thus,

where $C > 0$ is independent of $h, \tau, \nu$ and $\kappa$.

For $J_2^{n+1}$, it follows from the Hölder inequality and (2.6) that

then

where $C > 0$ is independent of $h, \tau, \nu, \kappa$ and $\epsilon_1 > 0$ is some small constant determined later.

For $J_3^{n+1}$, in terms of the skew-symmetric property (2.1), (2.18), the Hölder inequality and the regularity assumption (3.10), there holds

where $C > 0$ is independent of $h, \tau, \nu$ and $\kappa$. Thus, the third term can be estimated by

where $C > 0$ is independent of $h, \tau, \nu$ and $\kappa$.

A similar method to bound $J_4^{n+1}$ leads to

and

where $C > 0$ is independent of $h, \tau, \nu$ and $\kappa$.

For $J_5^{n+1}$, using the Hölder inequality, Young inequality and the regularity assumption (3.10), we have

then

where $C > 0$ is independent of $h, \tau, \nu$ and $\kappa$.

By a similar method, we estimate the last term by

where we use the Sobolev imbedding inequality and $C > 0$ is independent of $h, \tau, \nu$ and $\kappa$. Taking the sum gives

where $C > 0$ is independent of $h, \tau, \nu$ and $\kappa$.

Substituting the estimates (3.13)–(3.18) into (3.12), we obtain

where $C > 0$ is independent of $h, \tau, \nu$ and $\kappa$.

Next, we estimate $e_\theta^{n+1}$. Taking $\psi_h = 2\tau e_\theta^{n+1}$ in (3.9), summing up the resulting equation from $n = 0$ to $n = m$ and noticing the skew-symmetric property (2.1), we have

where we noted $e_{\theta}^0 = 0$. By the Taylor's formula and the Hölder inequality, it is easy to prove that

by noticing

where $C > 0$ is independent of $h, \tau$, $\nu$ and $\kappa$.

For last term $\lambda_3^{n+1}$, in terms of the following splitting:

we estimate it by

Furthermore, we get

where $C > 0$ is independent of $h, \tau$, $\nu$ and $\kappa$.

Substituting the estimates (3.21)–(3.23) into (3.20), we get

where $C > 0$ is independent of $h, \tau$, $\nu$ and $\kappa$.

Taking the sum of (3.19) and (3.24), we have

where $C > 0$ is independent of $h, \tau$, $\nu$ and $\kappa$.

To uniformly bound the norm $\|\theta_h^n\|_{L^\infty}$, we use the method of mathematical induction. Taking $m = 0$ in (3.25) and noting $e_h^0 = 0$, $e_\theta^0 = 0$ and $\|\theta_h^0\|_{L^\infty} = \|R_h\theta_0\|_{L^\infty}\leq \|\theta_0\|_{H^3}$, we get

For sufficiently small $\tau$ with $C\tau < 1/2$, we select small parameter $\epsilon_1\leq 2/3$ in (3.26). We then conclude that there exists some $C_1 > 0$, which is independent of $h, \tau$, $\nu$ and $\kappa$ such that

Thus, the error estimate (3.11) holds for $m = 0$ by taking $C_0 > C_1$. Now, we assume that the error estimate (3.11) holds for $m\leq k-1$ with $1\leq k\leq N-1$. Under the time step condition $\tau\leq Ch^2$, we have

By the inverse inequality, we get

for sufficiently small $h$ with $C_0h^{1/2} \leq 1$. Thus,

As a result, resetting $m$ by $k$ in (3.25), we have

We select small parameter $\epsilon_1\leq 2/3$ in (3.29), then

Applying the discrete Gronwall inequality to the above inequality, we conclude that there exists some $C_2 > 0$, which is independent of $h, \tau$, $\nu$ and $\kappa$ such that

for sufficiently small $\tau$. Thus, we prove that the error estimate (3.11) also holds for $m = k$ by taking $C_0\geq C_2$, and we finish the mathematical induction. □

Based on (3.11) in Theorem 3.1, we get the following error estimate for the first-order Euler grad-div finite element scheme (3.1)–(3.2).

Theorem 3.2 Under the assumptions in Theorem 3.1, we have

with $0\leq m\leq N-1$, where $C > 0$ is independent of $h, \tau$, $\nu$ and $\kappa$.

4.   The BDF2 grad-div stabilization finite element approximation

In this section, we consider the second-order BDF2 fully discrete scheme. Based on the extrapolation method and grad-div stabilization method, we propose the following linearized second-order finite element scheme for $1\leq n\leq N-1$:

Step Ⅰ: Find $\theta_h^{n+1}\in Y_h$ by

Step Ⅱ: Find $({{\textbf{u}}}_{h}^{n+1}, \theta_{h}^{n+1})\in {{\textbf{V}}}_{h}\times Q_{h}$ by

Remark 4.1. Numerical solution $({{\textbf{u}}}_h^1, \theta_h^1)\in {{\textbf{V}}}_{0h}\times Y_h$ can be solved from the Euler finite element scheme (3.1)–(3.2) in the above section. Usually, since the one-step error is first higher order, then one has

where $C_3 > 0$ is independent of $h, \tau$, $\nu$ and $\kappa$.

By a similar proof for Lemma 3.1, we can get the following unconditional stabilities of the numerical scheme (4.1)–(4.2), which imply the existence and uniqueness of numerical solution $({{\textbf{u}}}_h^{n+1}, p_h^{n+1}, \theta_h^{n+1})$ to the BDF2 scheme (4.1)–(4.2).

Lemma 4.1. For $1\leq n\leq N-1$ and all $\tau > 0$, $h > 0$, the BDF2 finite element scheme (4.1) and (4.2) has a unique solution $\theta_{h}^{n+1}\in Y_{h}$ and $({{\textbf{u}}}_{h}^{n+1}, p_{h}^{n+1})\in {{\textbf{V}}}_{h}\times Q_{h}$. Moreover, the discrete energy inequalities hold:

and

for all $1\leq m\leq N-1$, where $C > 0$ is independent of $h$ and $\tau$.

In terms of the discrete inf-sup condition (2.5), we rewrite the discrete variational problem (4.2) as: find ${{\textbf{u}}}_{h}^{n+1}\in {{\textbf{V}}}_{0h}$ by

On the other hand, let ${{\textbf{s}}}_{h}^{0} = {{\textbf{u}}}_{h}^{0}$. For $1\leq n\leq N-1$, ${{\textbf{s}}}_{h}^{n+1}$ and $\theta^{n+1}$ satisfy

and

In (4.7), we use $\nabla\cdot{{\textbf{u}}}^{n+1} = 0$ in $\Omega$ and $(\nabla\cdot{{\textbf{v}}}_h, \pi_hp^{n+1}) = 0$.

Next, we give error equations corresponding to the BDF2 scheme. Subtracting (4.6) from (4.7) leads to

where

Subtracting (4.1) from (4.8) leads to

where

The main result in this section is presented in the following theorem.

Theorem 4.1. Under the assumptions in (3.10), we further assume that

Then, for the second-order BDF2 grad-div scheme (4.1)–(4.2), under the time step condition $\tau\leq Ch$, when $h$ and $\tau$ are sufficiently small, we have the following error estimate:

with $0\leq m\leq N-1$, where $\widehat C_0 > 0$ is independent of $h, \tau$, $\nu$ and $\kappa$.

Proof. Taking ${{\textbf{v}}}_h = 4\tau e_h^{n+1}$ in (4.9) and summing up the resulting equation from $n = 1$ to $n = m$, we have

where $C > 0$ is independent of $h, \tau$, $\nu, \kappa$ and we use

We estimate the right-hand side of (4.13) as follows. Since

by using the Taylor's formula and the regularity assumption (4.11), one has

which results in

where $C > 0$ is independent of $h, \tau$, $\nu$ and $\kappa$. From the Hölder inequality, we have

Thus, we can obtain

where $C > 0$ is independent of $h, \tau$, $\nu$ and $\kappa$.

From the Hölder inequality and (2.6), we estimate $(I_{2}^{n+1}, e_{h}^{n+1})$ by

Then, by the Young inequality, it is easy to see that

where $C > 0$ is independent of $h, \tau, \nu, \kappa$ and $\epsilon_2 > 0$ is a small constant determined later.

For the term of $I_3^{n+1}$, according to the Hölder inequality, (2.6), (2.18), the regularity assumption (3.10) and the Taylor formula with integral remainder

where $(t-t_{n})_+ = \max\{t-t_n, 0\}$, there holds

Thus, we have

where $C > 0$ is independent of $h, \tau$, $\nu$ and $\kappa$.

A similar method leads to

and

where $C > 0$ is independent of $h, \tau$, $\nu$ and $\kappa$.

For $I_5^{n+1}$, using (4.16), the Hölder inequality and the regularity assumption (3.10), we have

{then}

where $C > 0$ is independent of $h, \tau$, $\nu$ and $\kappa$.

We estimate the last term $I_6^{n+1}$ by

where $C > 0$ is independent of $h, \tau$, $\nu$ and $\kappa$. Taking the sum of the above inequality gives

where $C > 0$ is independent of $h, \tau$, $\nu$ and $\kappa$.

Substituting the estimates (4.14)–(4.20) into (4.13), we obtain

where $C > 0$ is independent of $h, \tau$, $\nu, \kappa$ and the estimate (4.3) is used.

Next, we estimate $e_\theta^{n+1}$. Taking $\psi_h = 4\tau e_\theta^{n+1}$ in (4.10), summing up the resulting equation from $n = 1$ to $n = m$ and noticing the skew-symmetric property (2.1), we have

where we noted $e_{\theta}^0 = 0$.

By the regularity assumptions (3.10) and (4.11), it is easy to prove that

then we have

where $C > 0$ is independent of $h, \tau$, $\nu$ and $\kappa$.

For last term $X_3^{n+1}$, we estimate it by

where the skew-symmetric property (2.1) is used. We have

where the $C > 0$ is independent of $h, \tau$, $\nu$ and $\kappa$.

Substituting the estimates (4.23)–(4.25) into (4.22) and using (4.3), we get

where $C > 0$ is independent of $h, \tau$, $\nu$ and $\kappa$.

Taking the sum of (4.21) and (4.26), we have

where $C > 0$ is independent of $h, \tau$, $\nu$ and $\kappa$.

Next, we bound uniformly the norm $\|\theta_h^n\|_{L^\infty}$ by the method of mathematical induction as in the proof of Theorem 3.1.

According to (4.3), we can see that (4.12) holds for $m = 0$ if we take $\widehat C_0\geq C_3$. Now, we assume that (4.12) also holds for $m\leq k-1$ with $1\leq k\leq N-1$. Under the condition $\tau\leq Ch$ and the inverse inequality, one has

for sufficiently small $h$ with $\widehat C_0 h^{1/2}\leq 1$, which further gives

To finish the mathematical induction, we need to prove that (4.12) also holds for $m\leq k$. As a result of (4.28), we have

by setting $m = k$ in (4.27). Select small $\epsilon_2\leq 1$, then we have

where $C > 0$ is independent of $h, \tau$, $\nu$ and $\kappa$.

Applying the discrete Grönwall inequality to the above inequality, we conclude that there exists some $C_4 > 0$, which is independent of $h, \tau$, $\nu$ and $\kappa$, such that

for sufficiently small $\tau$. Thus, we prove that the error estimate (4.12) also holds for $m\leq k$ by taking $\widehat C_0\geq C_4$, and we finish the mathematical induction.□

Based on (4.12) in Theorem 4.1, we get the following error estimate for the second-order BDF2 grad-div finite element scheme (4.1)–(4.2).

Theorem 4.2 Under the assumptions in Theorem 4.1, we have

with $0\leq m\leq N-1$, where $C > 0$ is independent of $h, \tau$, $\nu$ and $\kappa$.

5.   Numerical results

In this section, we present numerical results of the second-order BDF2 grad-div stabilization method for the penetrative convection problems (1.1)–(1.5) with small viscosity coefficient $\nu$. For simplicity, we only consider the 2D problem with the domain $\Omega = [0, 1]^2$. To check the convergence rates, we take the appropriate ${{\textbf{f}}}$ and $g$ such that the exact solution $({{\textbf{u}}}, p, \theta)$ is given by

In numerical computation, we take the physical parameters $\kappa = 10^{-1}, \gamma_{1} = 10^{-1}, \gamma_{2} = 10^{-1}$, the stabilized parameter $\beta = 0.1$ and the final time $T = 1$. In addition, we denote

In terms of error estimates in Theorems 4.1 and 4.2, we have the second-order convergence rates of the velocity and temperature if we select $\tau = \mathcal O(h)$. We use a uniform mesh on $\Omega$ with the spatial mesh size $h$ in each direction. By taking gradually decreasing mesh sizes $h = 1/4, 1/8, \cdots, 1/128$, we choose different iteration numbers $N = 1/h$ such that the time step size $\tau = h$. We present numerical results with small viscosity coefficients $\nu = 10^{-3}$ and $\nu = 10^{-4}$ in Tables 1 and 2, respectively, from which we can see that the second-order convergence rates are reached. Thus, these numerical results are in good agreement with the theoretical analysis.

On the other hand, the second-order convergence rates in the $L^2$ norm is sub-optimal since we use the $P_2$ element to approximate ${{\textbf{u}}}$ and $\theta$. Theoretically, the optimal error estimate should be

for $1\leq n\leq N$. To check (5.1), we take $\tau = h^{3/2}$ such that

By taking different mesh sizes $h = 1/2^2, 1/3^2, \cdots, 1/7^2$, we give numerical results in Tables 3 and 4 with $\nu = 10^{-3}$ and $\nu = 10^{-4}$, respectively. We can see that the almost third-order convergence rates in $L^2$ norm are reached. Thus, the $L^2$ error estimates in (4.29) is not optimal. How to prove the optimal error estimates in $L^2$ norm will be reported in future work.

6.   Conclusions

In this paper, we studied the first-order Euler and second-order BDF2 finite element schemes for the approximation of the the time-dependent penetrative convection equations. In designing of numerical scheme, we used the grad-div stabilization method to overcome instabilities from the high Reynolds number. Main advantages of the proposed schemes are two folds. One is that they both are unconditionally stable without any condition of the time step and mesh size. Another is that one only needs to solve linearized systems at each time step. In terms of the analysis technique in ^[18,19], uniform error estimates in $L^2$ norm were derived in which the constants are independent of inverse powers of the viscosity coefficient and thermal conductivity coefficient.

In some related models, the nonlinear term $2\nu\mathcal D({{\textbf{u}}}):\mathcal D({{\textbf{u}}})$, which is used to describe the viscous dissipation, appears in the temperature equation (1.3) ^[2,9,47]. Here, $\mathcal D({{\textbf{u}}})$ is the strain tensor having components $\mathcal D_{ij}({{\textbf{u}}}) = \frac 12 (\frac{\partial u_i}{\partial x_j}+\frac{\partial u_j}{\partial x_i})$. However, this strongly nonlinear term will result in much more difficulties in well-posedness analysis and numerical analysis. Thus, we consider the Eq (1.3) without the viscous dissipation in this paper. On the other hand, different modified Boussinesq approximations for some practical problems have been considered, such as the heat and mass transfer models in ^[10,48]. In future works, we can use the grad-div stabilization method for these modified models.

On the other hand, for the reason of simplicity, we use the zero Dirichlet boundary condition for the temperature. We remark that uniform error estimates derived in this paper can be extended to the problem with the mixed boundary condition for the temperature in ^[19], where the mixed boundary condition is composed by zero Dirichlet and Neumann boundary conditions.

Use of AI tools declaration

The authors declare they have not used Artificial Intelligence (AI) tools in the creation of this article.

Acknowledgments

This work was supported by National Natural Science Foundation of China (No. 11771337).

Conflict of interest

On behalf of all authors, the corresponding author states that there is no conflict of interest.