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

: Based on the grad-div stabilization method, the first-order backward Euler and second-order BDF2 finite element schemes were studied for the approximations of the time-dependent penetrative convection equations. The proposed schemes are both unconditionally stable. We proved the error bounds of the velocity and temperature in which the constants are independent of inverse powers of the viscosity and thermal conductivity coe ffi cients when the Taylor-Hood element and P 2 element are used in finite element discretizations. Finally, numerical experiments with high Reynolds numbers were shown to confirm the theoretical results.


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 graddiv 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: ∇ • u = 0, for (x, t) ∈ Ω × (0, T ], (1.2) u = 0, θ = 0, for (x, t) ∈ ∂Ω × (0, T ], (1.4) u(x, 0) = u 0 (x), θ(x, 0) = θ 0 (x), for x ∈ Ω, (1.5) where Ω is a bounded and convex polyhedral domain in R 3 with boundary ∂Ω, and (0, T ] is a finite time interval.In the above system, the unknown (u, p, θ) denotes the velocity of the fluid, the pressure and temperature, respectively.ν > 0 and κ > 0 represent the viscosity coefficient and thermal conductivity parameter.f and g are two given functions.i 3 is the unit basis vector given by i 3 = (0, 0, 1) T , and (γ 1 θ + γ 2 θ 2 )i 3 in (1.1) denotes the buoyancy term.The constraint ∇ • 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., γ 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 secondorder 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.

Materials and methods
For k ∈ N + and 1 ≤ p ≤ +∞, let L p (Ω) and W k,p (Ω) denote the standard Lebesgue space and Sobolev space, respectively.When p = 2, W k,2 (Ω) is the Hilber space H k (Ω).The norms in L p (Ω), W k,p (Ω) and H k (Ω) are defined as the classical senses (cf.[1]) and are denoted by ∥ • ∥ L p , ∥ • ∥ W k,p and ∥ • ∥ H k .We define H k 0 (Ω) to be the subspace of H k (Ω) of functions with zero trace on ∂Ω.The dual space of H 1 0 (Ω) is denoted by H −1 (Ω).The boldface notations L p (Ω), W k,p (Ω) and H k (Ω) are used to denote the vector-value spaces L p (Ω) 3 , W k,p (Ω) 3 and H k (Ω) 3 , respectively.The corresponding norms are denoted by For simplicity, we denote Introduce the trilinear forms c 1 (•, •, •) and c 2 (•, •, •), which are given by for u i ∈ V with i = 1, 2, 3 and θ j ∈ Y with j = 1, 2. By the integration by parts, it is easy to check that the trilinear forms c 1 (•, •, •) and c 2 (•, •, •) satisfy the skew-symmetric properties, i.e., For given f ∈ L 2 (Ω), g ∈ L 2 (Ω), u 0 ∈ L 2 (Ω) and θ 0 ∈ L 2 (Ω), in terms of ther skew-symmetric properties (2.1), there holds the following energy inequalities: T be a uniform partition of the time interval [0, T ] with time step τ = T/N and t n = nτ with 0 ≤ n ≤ N. Denote u n = u(t n ), p n = p(t n ), θ n = θ(t n ), f n = f(t n ) and g n = g(t n ).For any sequence of functions { f n } N n=0 , denote the backward Euler discretization and BDF2 discretization by and where f n is the standard extrapolation formula.For the discrete time derivative D 2,τ , there holds the telescope formula [35]: (2.4) We define the finite element spaces as follows.Let T h = {K j } L j=1 be a quasi-uniform tetrahedron partition of Ω with the mesh size h = max {diam K j , j = 1, • • • , L}.In finite element discretization of (1.1)-(1.3),we use the Taylor-Hood (P 2 −P 1 ) finite element for the velocity u and pressure p, and the piece-wise quadratic finite element for the temperature θ.The corresponding finite element subspaces of V, Q and Y are denoted by 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 β 0 > 0 is independent of the mesh size h.Denote the space of discrete divergence-free functions by In the sequel, let π 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 V h .The following bound can be found in [4]: where n > 3/p when 1 < p ≤ ∞ and n ≥ 3 when p = 1.Next, we recall some known results about finite element approximations of Stokes problem in [17,18].Consider the Stokes problem: (2.10) Let us denote (u h , p h ) ∈ V h × 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, ν and κ.It is clear that the error estimates of the velocity depend on the negative power of ν.To prove error bounds with constants being independent of the negative power of ν, we recall the specific projection of (u, p) introduced in [17] that we will denote by s h defined by Let (u, p, θ) be the solution to the system (1.1)-(1.5)with u ∈ L ∞ ((0, T ]; ) and u t ∈ L ∞ ((0, T ]; H 1 (Ω)).We note that (u, 0) is the solution to the Stokes problem (2.10) with the right-hand side F given by (2.14) Let (s h , l h ) ∈ V h × Q h be the corresponding mixed finite element approximation of (u, 0).Then, from (2.11)-(2.13),we have where C > 0 is independent of h and ν.In terms of the inverse inequality in the theory of the finite element method and ∥I h u∥ L ∞ ≤ C∥u∥ L ∞ , we have where C > 0 is independent of h, ν and κ.Furthermore, there holds where C > 0 is independent of h, ν and κ.
Finally, we recall the discrete Gronwall inequality established in [24].Lemma 2.1.Let a k , b k , c k and γ k , for integers k ≥ 0, be the nonnegative numbers such that

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 ≤ n ≤ N − 1.
Step I: with the initial iteration values θ 0 h = R h θ 0 , u 0 h = I h u 0 and for any for any (v h , q h ) ∈ V h × Q h , where β > 0 is the stabilization parameter.The detailed discussion about the choice of β can be found in [28].

Summing up the above inequality from
2), we have Summing up the above inequality from n = 0 to n = m and noticing (3.3), we get (3.4).
The s h (t) is the finite element solution to the Stokes problem (2.10) with the right-hand side (2.14) and Thanks to the discrete inf-sup condition (2.5), the finite element scheme (3.2) is equivalent to: Find For simplicity, we take Next, we give error equations.Subtracting (3.5) from (3.6) leads to where Subtracting (3.1) from (3.7), we have

AIMS Mathematics
Volume 9, Issue 1, 453-480. where The main result in this section is presented in the following theorem.Theorem 3.1.Let (u, p, θ) be the unique solution to the time-dependent penetrative convection problems (1.1)-(1.5)and satisfy the following regularity assumptions: (3.10) Then, for the first-order Euler grad-div scheme (3.1)-(3.2) under the time step condition τ ≤ Ch 2 , when h and τ are sufficiently small, we have the following error estimate: ) and summing up the resulting equation from n = 0 to n = m, we have (J n+1 j , e n+1 h ). (3.12) We estimate the right-hand side of (3.12) as follows.For J n+1 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 where C > 0 is independent of h, τ, ν and κ.
For J n+1 2 , it follows from the Hölder inequality and (2.6) that where C > 0 is independent of h, τ, ν, κ and ϵ 1 > 0 is some small constant determined later.For J n+1 3 , 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, τ, ν and κ.Thus, the third term can be estimated by 2τ m n=0 where C > 0 is independent of h, τ, ν and κ.
A similar method to bound J n+1 4 leads to and 2τ m n=0 where C > 0 is independent of h, τ, ν and κ.
For J n+1 5 , using the Hölder inequality, Young inequality and the regularity assumption (3.10), we have where C > 0 is independent of h, τ, ν and κ.By a similar method, we estimate the last term by where we use the Sobolev imbedding inequality and C > 0 is independent of h, τ, ν and κ.Taking the sum gives where C > 0 is independent of h, τ, ν and κ.Substituting the estimates (3.13)-(3.18)into (3.12),we obtain where C > 0 is independent of h, τ, ν and κ.
Next, we estimate e n+1 θ .Taking ψ h = 2τe 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 0 θ = 0.By the Taylor's formula and the Hölder inequality, it is easy to prove that 2τ by noticing where C > 0 is independent of h, τ, ν and κ.
Taking the sum of (3.19) and (3.24), we have where C > 0 is independent of h, τ, ν and κ.
To uniformly bound the norm ∥θ n h ∥ L ∞ , we use the method of mathematical induction.Taking m = 0 in (3.25) and noting e 0 h = 0, e 0 θ = 0 and (3.26) For sufficiently small τ with Cτ < 1/2, we select small parameter ϵ 1 ≤ 2/3 in (3.26).We then conclude that there exists some C 1 > 0, which is independent of h, τ, ν and κ 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 ≤ k − 1 with 1 ≤ k ≤ N − 1.Under the time step condition τ ≤ Ch 2 , we have By the inverse inequality, we get for sufficiently small h with C 0 h 1/2 ≤ 1.Thus, As a result, resetting m by k in (3.25), we have (3.29) We select small parameter ϵ 1 ≤ 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, τ, ν and κ such that for sufficiently small τ.Thus, we prove that the error estimate (3.11) also holds for m = k by taking C 0 ≥ 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 graddiv finite element scheme (3.1)-(3.2).
Theorem 3.2 Under the assumptions in Theorem 3.1, we have with 0 ≤ m ≤ N − 1, where C > 0 is independent of h, τ, ν and κ.

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 secondorder finite element scheme for 1 ≤ n ≤ N − 1: Step II: Remark 4.1.Numerical solution (u 1 h , θ 1 h ) ∈ V 0h ×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, τ, ν and κ.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 (u n+1 h , p n+1 h , θ n+1 h ) to the BDF2 scheme (4.1)-(4.2).Lemma 4.1.For 1 ≤ n ≤ N − 1 and all τ > 0, h > 0, the BDF2 finite element scheme (4.1) and (4.2) has a unique solution θ n+1 h ∈ Y h and (u n+1 h , p n+1 h ) ∈ V h × Q h .Moreover, the discrete energy inequalities hold: and for all 1 ≤ m ≤ N − 1, where C > 0 is independent of h and τ.
In terms of the discrete inf-sup condition (2.5), we rewrite the discrete variational problem (4.2) as: On the other hand, let In (4.7), we use 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 τ ≤ Ch, when h and τ are sufficiently small, we have the following error estimate: with 0 ≤ m ≤ N − 1, where C 0 > 0 is independent of h, τ, ν and κ.
Proof.Taking v h = 4τe n+1 h in (4.9) and summing up the resulting equation from n = 1 to n = m, we have where C > 0 is independent of h, τ, ν, κ 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, τ, ν and κ.From the Hölder inequality, we have . Thus, we can obtain where C > 0 is independent of h, τ, ν and κ.
From the Hölder inequality and (2.6), we estimate (I n+1 2 , e n+1 h ) by Then, by the Young inequality, it is easy to see that where C > 0 is independent of h, τ, ν, κ and ϵ 2 > 0 is a small constant determined later.
For the term of I n+1 3 , 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, τ, ν and κ.
A similar method leads to and 4τ m n=1 where C > 0 is independent of h, τ, ν and κ.For I n+1 5 , using (4.16), the Hölder inequality and the regularity assumption (3.10), we have where C > 0 is independent of h, τ, ν and κ.
We estimate the last term I n+1 6 by where C > 0 is independent of h, τ, ν and κ.Taking the sum of the above inequality gives 4τ m n=0 where C > 0 is independent of h, τ, ν and κ.
For last term X n+1 3 , we estimate it by where the skew-symmetric property (2.1) is used.We have where the C > 0 is independent of h, τ, ν and κ.Substituting the estimates (4.23)-(4.25)into (4.22) and using (4.3), we get where C > 0 is independent of h, τ, ν and κ.
Taking the sum of (4.21) and (4.26), we have where C > 0 is independent of h, τ, ν and κ.
Next, we bound uniformly the norm ∥θ n h ∥ L ∞ 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 C 0 ≥ C 3 .Now, we assume that (4.12) also holds for m ≤ k − 1 with 1 ≤ k ≤ N − 1.Under the condition τ ≤ Ch and the inverse inequality, one has for sufficiently small h with C 0 h 1/2 ≤ 1, which further gives To finish the mathematical induction, we need to prove that (4.12) also holds for m ≤ k.As a result of (4.28), we have by setting m = k in (4.27).Select small ϵ 2 ≤ 1, then we have where C > 0 is independent of h, τ, ν and κ.
Applying the discrete Grönwall inequality to the above inequality, we conclude that there exists some C 4 > 0, which is independent of h, τ, ν and κ, such that for sufficiently small τ.Thus, we prove that the error estimate (4.12) also holds for m ≤ k by taking C 0 ≥ 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 ≤ m ≤ N − 1, where C > 0 is independent of h, τ, ν and κ.

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 ν.For simplicity, we only consider the 2D problem with the domain Ω = [0, 1] 2 .To check the convergence rates, we take the appropriate f and g such that the exact solution (u, p, θ) is given by In numerical computation, we take the physical parameters κ = 10 −1 , γ 1 = 10 −1 , γ 2 = 10 −1 , the stabilized parameter β = 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 τ = O(h).We use a uniform mesh on Ω with the spatial mesh size h in each direction.By taking gradually decreasing mesh sizes h = 1/4, 1/8, • • • , 1/128, we choose different iteration numbers N = 1/h such that the time step size τ = h.We present numerical results with small viscosity coefficients ν = 10 −3 and ν = 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.
Table 1.Numerical errors and convergence rates with ν = 10 −3 and τ = h.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 u and θ.Theoretically, the optimal error estimate should be (5.1) for 1 ≤ n ≤ N. To check (5.1), we take τ = h 3/2 such that By taking different mesh sizes h = 1/2 2 , 1/3 2 , • • • , 1/7 2 , we give numerical results in Tables 3 and 4 with ν = 10 −3 and ν = 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.

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νD(u) : D(u), which is used to describe the viscous dissipation, appears in the temperature equation (1.3) [2,9,47].Here, D(u) is the strain tensor having components D i j (u) = 1 2 ( ∂u i ∂x j + ∂u j ∂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.