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 Ω is a bounded and convex polyhedral domain in R3 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. i3 is the unit basis vector given by i3=(0,0,1)T, and (γ1θ+γ2θ2)i3 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 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∈N+ and 1≤p≤+∞, let Lp(Ω) and Wk,p(Ω) denote the standard Lebesgue space and Sobolev space, respectively. When p=2, Wk,2(Ω) is the Hilber space Hk(Ω). The norms in Lp(Ω), Wk,p(Ω) and Hk(Ω) are defined as the classical senses (cf.[1]) and are denoted by ‖⋅‖Lp, ‖⋅‖Wk,p and ‖⋅‖Hk. We define Hk0(Ω) to be the subspace of Hk(Ω) of functions with zero trace on ∂Ω. The dual space of H10(Ω) is denoted by H−1(Ω). The boldface notations Lp(Ω), Wk,p(Ω) and Hk(Ω) are used to denote the vector-value spaces Lp(Ω)3, Wk,p(Ω)3 and Hk(Ω)3, respectively. The corresponding norms are denoted by ‖⋅‖Lp, ‖⋅‖Wk,p and ‖⋅‖Hk.
For simplicity, we denote V=H10(Ω), Y=H10(Ω),
Introduce the trilinear forms c1(⋅,⋅,⋅) and c2(⋅,⋅,⋅), which are given by
for ui∈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 c1(⋅,⋅,⋅) and c2(⋅,⋅,⋅) satisfy the skew-symmetric properties, i.e.,
which give
For given f∈L2(Ω),g∈L2(Ω), u0∈L2(Ω) and θ0∈L2(Ω), in terms of ther skew-symmetric properties (2.1), there {holds} the following energy inequalities:
and
for 0≤t≤T.
Let 0=t0<t1<⋯<tN=T be a uniform partition of the time interval [0,T] with time step τ=T/N and tn=nτ with 0≤n≤N. Denote un=u(tn), pn=p(tn), θn=θ(tn),
fn=f(tn) and gn=g(tn). For any sequence of functions {fn}Nn=0, denote the backward Euler discretization and BDF2 discretization by
and
where ˆfn is the standard extrapolation formula. For the discrete time derivative D2,τ, there holds the telescope formula [35]:
We define the finite element spaces as follows. Let Th={Kj}Lj=1 be a quasi-uniform tetrahedron partition of Ω with the mesh size h=max{diamKj, j=1,⋯,L}. In finite element discretization of (1.1)–(1.3), we use the Taylor-Hood (P2−P1) 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 Vh, Qh and Yh, 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 Rh be the L2 orthogonal and Ritz orthogonal projections onto Qh and Yh, respectively. Then the following error bounds and stability hold [4,43]:
Let Ih be the Lagrange interpolation operator onto Vh. 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:
Let us denote (uh,ph)∈Vh×Qh, 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 sh defined by
Let (u,p,θ) be the solution to the system (1.1)–(1.5) with u∈L∞((0,T];V∩H3(Ω)), p∈L∞((0,T];Q∩H2(Ω)),
θ∈L∞((0,T];Y∩H3(Ω)) and ut∈L∞((0,T];H1(Ω)). We note that (u,0) is the solution to the Stokes problem (2.10) with the right-hand side F given by
Let (sh,lh)∈Vh×Qh 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 ‖Ihu‖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 ak,bk,ck and γk, for integers k≥0, be the nonnegative numbers such that
Suppose that τγk<1 for all k and set σk=(1−τγ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 τ>0 with σ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≤n≤N−1.
Step Ⅰ: Find θn+1h∈Yh by
with the initial iteration values θ0h=Rhθ0, u0h=Ihu0 and for any ψh∈Yh.
Step Ⅱ: Find (un+1h,pn+1h)∈Vh×Qh by
for any (vh,qh)∈Vh×Qh, where β>0 is the stabilization parameter. The detailed discussion about the choice of β 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 (un+1h,pn+1h,θn+1h) to the Euler scheme (3.1)–(3.2).
Lemma 3.1. For 0≤n≤N−1 and all τ>0, h>0, the finite element scheme (3.1) and (3.2) has a unique solution θn+1h∈Yh and (un+1h,pn+1h)∈Vh×Qh. Moreover, the discrete energy inequalities hold:
and
for all 0≤m≤N−1, where C>0 is independent of τ and h.
Proof.
Setting ψh=2τθn+1h 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 (vh,qh)=2τ(un+1h,pn+1h) 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 sh(t) is the finite element solution to the Stokes problem (2.10) with the right-hand side (2.14) and snh=sh(tn).
Thanks to the discrete inf-sup condition (2.5), the finite element scheme (3.2) is equivalent to: Find un+1h∈V0h by
For simplicity, we take s0h=u0h. For 0≤n≤N−1, sn+1h and θn+1 satisfy
and
In (3.6), we use ∇⋅un+1=0 in Ω and (∇⋅vh,πhpn+1)=0 due to vh∈V0h.
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 (u,p,θ) 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 τ≤Ch2, when h and τ are sufficiently small, we have the following error estimate:
with 0≤m≤N−1, where C0>0 is independent of h,τ, ν and κ.
Proof. Taking vh=2τen+1h 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 Jn+11, 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,τ,ν and κ.
For Jn+12, it follows from the Hölder inequality and (2.6) that
then
where C>0 is independent of h,τ,ν,κ and ϵ1>0 is some small constant determined later.
For Jn+13, 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
where C>0 is independent of h,τ,ν and κ.
A similar method to bound Jn+14 leads to
and
where C>0 is independent of h,τ,ν and κ.
For Jn+15, using the Hölder inequality, Young inequality and the regularity assumption (3.10), we have
then
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 en+1θ. Taking ψh=2τen+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 e0θ=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,τ, ν and κ.
For last term λn+13, in terms of the following splitting:
we estimate it by
Furthermore, we get
where C>0 is independent of h,τ, ν and κ.
Substituting the estimates (3.21)–(3.23) into (3.20), we get
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 ‖θnh‖L∞, we use the method of mathematical induction. Taking m=0 in (3.25) and noting e0h=0, e0θ=0 and ‖θ0h‖L∞=‖Rhθ0‖L∞≤‖θ0‖H3, we get
For sufficiently small τ with Cτ<1/2, we select small parameter ϵ1≤2/3 in (3.26).
We then conclude that there exists some C1>0, which is independent of h,τ, ν and κ such that
Thus, the error estimate (3.11) holds for m=0 by taking C0>C1. Now, we assume that the error estimate (3.11) holds for m≤k−1 with 1≤k≤N−1. Under the time step condition τ≤Ch2, we have
By the inverse inequality, we get
for sufficiently small h with C0h1/2≤1. Thus,
As a result, resetting m by k in (3.25), we have
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 C2>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 C0≥C2, 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≤m≤N−1, where C>0 is independent of h,τ, ν and κ.
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≤n≤N−1:
Step Ⅰ: Find θn+1h∈Yh by
Step Ⅱ: Find (un+1h,θn+1h)∈Vh×Qh by
Remark 4.1. Numerical solution (u1h,θ1h)∈V0h×Yh 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 C3>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 (un+1h,pn+1h,θn+1h) 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+1h∈Yh and (un+1h,pn+1h)∈Vh×Qh. 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: find un+1h∈V0h by
On the other hand, let s0h=u0h. For 1≤n≤N−1, sn+1h and θn+1 satisfy
and
In (4.7), we use ∇⋅un+1=0 in Ω and (∇⋅vh,πhpn+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 τ≤Ch, when h and τ are sufficiently small, we have the following error estimate:
with 0≤m≤N−1, where ˆC0>0 is independent of h,τ, ν and κ.
Proof. Taking vh=4τen+1h 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 (In+12,en+1h) 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 In+13, 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.