A Uzawa-Type Iterative Algorithm for the Stationary Natural Convection Model

Abstract

In this study, a Uzawa-type iterative algorithm is introduced and analyzed for solving the stationary natural convection model, where physical variables are discretized by utilizing a mixed finite element method. Compared with the common Uzawa iterative algorithm, the main finding is that the proposed algorithm produces weakly divergence-free velocity approximation. In addition, the convergence results of the proposed algorithm are provided, and numerical tests supporting the theory are presented.

1. Introduction

Arising both in nature and in engineering applications, the natural convection model is a coupled system of fluid flow governed by the incompressible Navier-Stokes equations and heat transfer governed by the energy equation. The natural convection problem has been a hot topic in heat transmission science for a long time, because it has been widely used in many fields of production and life, such as room ventilation, general heating, nuclear reaction systems, fire control, katabatic winds, atmospheric fronts, cooling of electronic equipment, natural ventilation, solar collectors, and so on [1,2,3]. In particular with nanofluids, the literature survey in [4] evidences the parameters governing the flow and heat behavior of fluids under natural convection and reveals that there are very few generalized correlations between heat transfer and wall heating conditions in enclosures.

Due to its practical significance, a considerable amount of researchers have put forward many efficient numerical methods to obtain the solution to this problem in different geometries [5,6,7,8,9,10]. For example, Boland and Layton [6,7] have proposed a Galerkin finite element method for the natural convection problem. Several iterative schemes based on the finite element method for the natural convection equations with different Rayleigh numbers have been studied in [9]. The coupled Navier-Stokes/temperature (or Boussinesq) equations [5] were solved by applying a divergence-free low order stabilized finite element method. A unified analysis approach of a local projection stabilization finite element method for solving natural convection problems was given by [8]. However, there still remain some important but challenging problems, especially solving the model effectively with the strong coupling between the velocity, pressure, and temperature fields and the saddle-point problem arising from finite element discretization.

As is known, the Uzawa method [11] is an efficient iterative algorithm for the saddle-point system. Since it is simple, efficient, and has minimal computer memory requirements, it has been widely used in computational science and engineering [12,13,14,15,16]. In particular, some Uzawa iterative methods were designed for the steady incompressible Navier-Stokes equations [17]. Further, the steady magnetohydrodynamic equations [18] and the steady natural convection equations [19] were solved by applying some Uzawa iterative algorithms. However, in these works, the weakly divergence-free constraint on the velocity was not enforced.

Recently, a Uzawa-type iterative algorithm [20] was designed for the coupled Stokes equations, where no saddle point system was required to be solved at each iteration step, and the weakly divergence-free velocity approximation was shown. Inspired by [20], in this article we propose and analyze a Uzawa-type iterative algorithm for the natural convection problem and obtain a numerical velocity, which satisfies the weakly divergence-free condition.

2. Preliminaries

Let $\Omega \subset {\mathbb{R}}^2$ be a bounded domain, which has a Lipschitz continuous boundary $\partial \Omega$ with a regular open subset $\Gamma$. Consider the following stationary natural convection problem. Seek the velocity $\mathbf{u}={(u_1(x),u_2(x))}^{\top }$, the pressure $p=p(x)$, and the temperature $T(x)$, such that

$$\begin{array}{cccc}\hfill \nabla p+(\mathbf{u}·\nabla )\mathbf{u}-Pr\Delta \mathbf{u}& =PrRa\mathbf{j}T,\nabla ·\mathbf{u}=0\hfill & \hfill \mathrm{in}& \Omega ,\hfill \end{array}$$

(1)

$$\begin{array}{cccc}\hfill \mathbf{u}& =0,\hfill & \hfill \mathrm{on}& \partial \Omega ,\hfill \end{array}$$

(2)

$$\begin{array}{cccc}\hfill -\kappa \Delta T& =\gamma -\mathbf{u}·\nabla T,\hfill & \hfill \mathrm{in}& \Omega ,\hfill \end{array}$$

(3)

$$\begin{array}{cccc}\hfill T& =0,\mathrm{on}\partial \Gamma ,\frac{\partial T}{\partial \mathbf{n}}=0,\hfill & \hfill \mathrm{on}& \partial \Omega \backslash \Gamma ,\hfill \end{array}$$

(4)

where $\gamma$ is the forcing function, $\mathbf{n}$ is the outward unit vector, and $\mathbf{j}={(0,1)}^{\top }$. In addition, the positive parameter $\kappa$ presents the thermal conductivity, $Pr$ is the Prandtl number, and $Ra$ is the Rayleigh number.

Next, in order to write the variational form of (1)–(4), we introduce the following necessary function spaces:

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill \mathbf{M}& =H_0^1{(\Omega )}^2=\{\mathbf{v}\in H^1{(\Omega )}^2:\mathbf{v}=0\mathrm{on}\partial \Omega \},\hfill \\ \hfill W& =L_0^2(\Omega )=\{q\in L^2(\Omega ):(q,1)=0\},Z=\{s\in H^1(\Omega ):s=0\mathrm{on}\Gamma \}.\hfill \end{array}\end{array}$$

Here, the space $L^2(\Omega )$ is endowed with $L^2$-scalar product $(·,·)$ and $L^2$-norm $\parallel ·\parallel$. In addition, the space $H^1(\Omega )$ is used to represent the standard definitions for Sobolev spaces $W^{m,p}(\Omega )$, $m,p>0$.

Moreover, we recall the Poincaré inequality [21] as follows:

$$\begin{array}{c}\hfill \begin{array}{c}\hfill \parallel \mathbf{v}\parallel \le C_p\parallel \nabla \mathbf{v}\parallel ,\forall \mathbf{v}\in \mathbf{M},\end{array}\end{array}$$

(5)

where $C_p$ is the Poincaré constant. Next, we denote two trilinear forms by

$$\begin{array}{c}\hfill \begin{array}{c}\hfill b_1(\mathbf{u};\mathbf{v},\mathbf{w})=((\mathbf{u}·\nabla )\mathbf{v},\mathbf{w})+\frac{1}{2}((\nabla ·\mathbf{u})\mathbf{v},\mathbf{w}),b_2(\mathbf{u};T,s)=(\mathbf{u}·\nabla T,s)+\frac{1}{2}((\nabla ·\mathbf{u})T,s),\end{array}\end{array}$$

which satisfy the following properties [7,22,23]

$$\begin{array}{c}\hfill \begin{array}{c}\hfill |b_1(\mathbf{u};\mathbf{v},\mathbf{w})|\le N\parallel \nabla \mathbf{u}\parallel \parallel \nabla \mathbf{v}\parallel \parallel \nabla \mathbf{w}\parallel ,|b_2(\mathbf{u};T,s)|\le \overline{N}\parallel \nabla \mathbf{u}\parallel \parallel \nabla T\parallel \parallel \nabla s\parallel ,\end{array}\end{array}$$

(6)

for all $\mathbf{u},\mathbf{v},\mathbf{w},\in \mathbf{M}$ and $T,s\in Z$. Here, N and $\overline{N}$ are two fixed positive constants.

With the above notations, the weak form of (1)–(4) reads as: find $(\mathbf{u},p,T)\in \mathbf{M}\times W\times Z$ such that

$$\begin{array}{cc}\hfill & Pr(\nabla \mathbf{u},\nabla \mathbf{v})+b_1(\mathbf{u};\mathbf{u},\mathbf{v})-(p,\nabla ·\mathbf{v})=PrRa(\mathbf{j}T,\mathbf{v}),\forall \mathbf{v}\in \mathbf{M},\hfill \end{array}$$

(7)

$$\begin{array}{cc}\hfill & (\nabla ·\mathbf{u},q)=0,\forall q\in W,\hfill \end{array}$$

(8)

$$\begin{array}{cc}\hfill & \kappa (\nabla T,\nabla s)+b_2(\mathbf{u};T,s)=(\gamma ,s),\forall s\in Z.\hfill \end{array}$$

(9)

The following existence and uniqueness of the solution to (6) are classical results.

Theorem 1

([7,19]). There exists at least a solution $(\mathbf{u},p,T)\in \mathbf{M}\times W\times Z$, which satisfies (7)–(9) and

$$\begin{array}{cc}\hfill \parallel \nabla T\parallel & \le {\kappa }^{-1}{\parallel \gamma \parallel }_{-1},\parallel \nabla \mathbf{u}\parallel \le C_p^{2}Ra{\kappa }^{-1}{\parallel \gamma \parallel }_{-1},\hfill \end{array}$$

where ${\parallel \gamma \parallel }_{-1}=\underset{s\in Z}{sup}\frac{|(\gamma ,s)|}{\parallel \nabla s\parallel }.$ Further, if $Pr$, $Ra$, κ, and γ satisfy the uniqueness condition

$$0<P{r}^{-1}\Lambda +\overline{\Lambda }<1,$$

where $\Lambda =C_p^{2}RaN{\kappa }^{-1}{\parallel \gamma \parallel }_{-1}$ and $\overline{\Lambda }=C_p^{2}Ra\overline{N}{\kappa }^{-2}{\parallel \gamma \parallel }_{-1}$, then the solution $(\mathbf{u},p,T)$ of (7)–(9) is unique.

Next, we consider a family of quasi-uniform and regular triangulations $K_h=\{K:{\cup }_{K\subset \Omega }\overline{K}=\overline{\Omega }\}$ with mesh size h, which is a partition of the domain $\Omega$. Then, we assume that the finite element subspace ${\mathbf{M}}_h\times W_h\times Z_h\subset \mathbf{M}\times W\times Z$

$${\mathbf{M}}_h=\{\mathbf{v}\in \mathbf{M}\cap C^0{(\overline{\Omega })}^2:\mathbf{v}{\mid }_K\in P_2{(K)}^2,\forall K\in K_h\},$$

$$W_h=\{q\in W\cap C^0(\overline{\Omega }):q{\mid }_K\in P_1(K),\forall K\in K_h\},$$

$$Z_h=\{s\in Z\cap C^0(\overline{\Omega }):s{\mid }_K\in P_2(K),\forall K\in K_h\},$$

where $P_i(K)$, $i=1,2$ is the set of all polynomials on K of a degree no more than i. As is known, the finite element subspaces ${\mathbf{M}}_h\times W_h$ satisfy the following discrete inf-sup condition [21]; for each $q\in W_h$, there exists $\mathbf{v}\in {\mathbf{M}}_h,\mathbf{v}\ne 0$ such that ${inf}_{q\in W_h}{sup}_{\mathbf{v}\in {\mathbf{M}}_h}\frac{|(\nabla ·\mathbf{v},q)|}{\parallel \nabla \mathbf{v}\parallel \parallel q\parallel }\ge \beta$, where the constant $\beta \in (0,1]$ is proven in [24].

Moreover, according to the above definition of the finite element subspaces, the finite element approximation for (7)–(9) is to seek $({\mathbf{u}}_h,p_h,T_h)\in {\mathbf{M}}_h\times W_h\times Z_h$ such that

$$\begin{array}{cc}\hfill & Pr(\nabla {\mathbf{u}}_h,\nabla \mathbf{v})+b_1({\mathbf{u}}_h;{\mathbf{u}}_h,\mathbf{v})-(p_h,\nabla ·\mathbf{v})=PrRa(\mathbf{j}{T}_h,\mathbf{v}),\forall \mathbf{v}\in {\mathbf{M}}_h,\hfill \end{array}$$

(10)

$$\begin{array}{cc}\hfill & (\nabla ·{\mathbf{u}}_h,q)=0,\forall q\in W_h,\hfill \end{array}$$

(11)

$$\begin{array}{cc}\hfill & \kappa (\nabla T_h,\nabla s)+b_2({\mathbf{u}}_h;T_h,s)=(\gamma ,s),\forall s\in Z_h.\hfill \end{array}$$

(12)

The following theorem is established for the stability of the finite element discretization.

Theorem 2

([6,9,25]). Under the assumptions of Theorem 1, the finite element discretization (10)–(12) has at least a solution $({\mathbf{u}}_h,p_h,T_h)\in {\mathbf{M}}_h\times W_h\times Z_h$, such that

$$\begin{array}{c}\hfill \parallel \nabla {\mathbf{u}}_h\parallel \le C_p^{2}Ra{\kappa }^{-1}{\parallel \gamma \parallel }_{-1},\parallel \nabla T_h\parallel \le {\kappa }^{-1}{\parallel \gamma \parallel }_{-1}.\end{array}$$

3. A Uzawa-Type Iterative Algorithm

In this section, we present a Uzawa-type iterative algorithm for solving the considered problem. Before showing the algorithm, we recall the common Uzawa iterative algorithm based on the mixed finite element method as follows Algorithm 1.

According to the above algorithm, we find that $(\nabla ·{\mathbf{u}}_h^{n+1},q)\ne 0$, which means that the divergence-free constraint on the velocity is not weakly enforced. In fact, from the finite element approximation (10)–(12), we have $(\nabla ·{\mathbf{u}}_h,q)=0$. Although it will result in a saddle problem, it produces weakly divergence-free velocity approximation. Hence, it is interesting to design a Uzawa-type iterative algorithm, which does not only retain the benefits of the common Uzawa iterative algorithm but also retains the velocity in a weakly divergence-free condition.

Algorithm 1: Uzawa iterative algorithm [19].

Step 1. Find initial guess $({\mathbf{u}}_h^0,p_h^0,T_h^0)\in {\mathbf{M}}_h\times W_h\times Z_h$ by $$\begin{array}{c}\hfill \left\{\begin{array}{cc}& Pr(\nabla {\mathbf{u}}_h^0,\nabla \mathbf{v})-(p_h^0,\nabla ·\mathbf{v})=PrRa(\mathbf{j}{T}_h^0,\mathbf{v}),\forall \mathbf{v}\in {\mathbf{M}}_h,\hfill \\ & (\nabla ·{\mathbf{u}}_h^0,q)=0,\forall q\in W_h,\hfill \\ & \kappa (\nabla T_h^0,\nabla s)=(\gamma ,s),\forall s\in Z_h.\hfill \end{array}\right.\end{array}$$ Step 2. Given a relaxation parameter $\rho >0$, find $({\mathbf{u}}_h^{n+1},p_h^{n+1},T_h^{n+1})\in {\mathbf{M}}_h\times W_h\times Z_h$ as solution of $$\begin{array}{c}\hfill \left\{\begin{array}{cc}& Pr(\nabla {\mathbf{u}}_h^{n+1},\nabla \mathbf{v})+b_1({\mathbf{u}}_h^n;{\mathbf{u}}_h^{n+1},\mathbf{v})-(p_h^n,\nabla ·\mathbf{v})=PrRa(\mathbf{j}{T}_h^{n+1},\mathbf{v}),\forall \mathbf{v}\in {\mathbf{M}}_h,\hfill \\ & (p_h^{n+1},q)=(p_h^n,q)-\rho (\nabla ·{\mathbf{u}}_h^{n+1},q),\forall q\in W_h,\hfill \\ & \kappa (\nabla T_h^{n+1},\nabla s)+b_2({\mathbf{u}}_h^n;T_h^{n+1},s)=(\gamma ,s),\forall s\in Z_h.\hfill \end{array}\right.\end{array}$$

In order to make the velocity of Uzawa algorithm have a weakly divergence-free property, let g be a gauge variable [26] and $\mathbf{d}$ be a variable, such that $\mathbf{u}=\mathbf{d}+\nabla g$. If g and p satisfy an elliptic equation $Pr\Delta g=p$, then (1)–(4) can be rewritten as

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill -Pr\Delta \mathbf{d}+((\mathbf{d}+\nabla g)·\nabla )(\mathbf{d}+\nabla g)& =PrRa\mathbf{j}T,\hfill \\ \hfill \nabla ·\mathbf{d}& =-\Delta g,\hfill \\ \hfill -\kappa \Delta T+(\mathbf{d}+\nabla g)·\nabla T& =\gamma .\hfill \end{array}\end{array}$$

Furthermore, begin with $g^0=g^{-1}=0$ and ${\mathbf{d}}^0={\mathbf{u}}_h^0$. Repeat

$$\begin{array}{cc}\hfill -Pr\Delta {\mathbf{d}}^{n+1}+(({\mathbf{d}}^n+\nabla g^{n-1})·\nabla )({\mathbf{d}}^{n+1}+\nabla g^n)& =PrRa\mathbf{j}{T}^{n+1},\hfill \end{array}$$

(13)

$$\begin{array}{cc}\hfill \nabla ·{\mathbf{d}}^{n+1}& =-\Delta g^{n+1},\hfill \end{array}$$

(14)

$$\begin{array}{cc}\hfill -\kappa \Delta T^{n+1}+({\mathbf{d}}^n+\nabla g^{n-1})·\nabla T^{n+1}& =\gamma ,\hfill \end{array}$$

(15)

for $n=0,1,\dots$

Moreover, setting ${\widehat{\mathbf{u}}}^{n+1}={\mathbf{d}}^{n+1}+\nabla g^n$ in (13)–(15), we have

$$\begin{array}{cc}\hfill -Pr\Delta {\widehat{\mathbf{u}}}^{n+1}+({\widehat{\mathbf{u}}}^n·\nabla ){\widehat{\mathbf{u}}}^{n+1}+\nabla p^n& =PrRa\mathbf{j}{T}^{n+1},\hfill \end{array}$$

(16)

$$\begin{array}{cc}\hfill \nabla ·{\widehat{\mathbf{u}}}^{n+1}& =-\Delta {\hslash }^{n+1},\hfill \end{array}$$

(17)

$$\begin{array}{cc}\hfill -\nabla ·(\kappa \nabla T^{n+1})+({\widehat{\mathbf{u}}}^n·\nabla )T^{n+1}& =\gamma ,\hfill \end{array}$$

(18)

where ${\hslash }^{n+1}:=g^{n+1}-g^n$. So one obtains

$$p^{n+1}=P_r\Delta g^{n+1}=P_r\Delta {\hslash }^{n+1}+Pr\Delta g^n=P_r\Delta {\hslash }^{n+1}+p^n,$$

and

$${\mathbf{u}}^{n+1}={\mathbf{d}}^{n+1}+\nabla g^{n+1}={\widehat{\mathbf{u}}}^{n+1}-\nabla g^n+\nabla g^{n+1}={\widehat{\mathbf{u}}}^{n+1}+\nabla {\hslash }^{n+1}.$$

Now, we are ready to write the Uzawa-type finite element iterative algorithm as follows Algorithm 2.

Algorithm 2: Uzawa-type iterative algorithm.

Step 1. Obtain the initial guess $({\mathbf{u}}_h^0,p_h^0,T_h^0)\in {\mathbf{M}}_h\times W_h\times Z_h$ from step 1 of Algorithm 1. Step 2. Find $({\widehat{\mathbf{u}}}_h^{n+1},T_h^{n+1})\in {\mathbf{M}}_h\times Z_h$ as the solution of $$\begin{array}{cc}\hfill & \kappa (\nabla T_h^{n+1},\nabla s)+b_2({\widehat{\mathbf{u}}}_h^n;T_h^{n+1},s)=(\gamma ,s),\forall s\in Z_h,\hfill \end{array}$$ (19) $$\begin{array}{cc}\hfill & Pr(\nabla {\widehat{\mathbf{u}}}_h^{n+1},\nabla \mathbf{v})+b_1({\widehat{\mathbf{u}}}_h^n;{\widehat{\mathbf{u}}}_h^{n+1},\mathbf{v})-(p_h^n,\nabla ·\mathbf{v})=PrRa(\mathbf{j}{T}_h^{n+1},\mathbf{v}),\forall \mathbf{v}\in {\mathbf{M}}_h.\hfill \end{array}$$ (20) Step 3. Find ${\hslash }_h^{n+1}\in W_h$ as the solution of $$\begin{array}{c}\hfill \begin{array}{c}\hfill (\nabla {\hslash }_h^{n+1},\nabla q)=(\nabla ·{\widehat{\mathbf{u}}}_h^{n+1},q),q\in W_h.\end{array}\end{array}$$ (21) Step 4. Compute ${\mathbf{u}}_h^{n+1}$ with ${\mathbf{u}}_h^{n+1}={\widehat{\mathbf{u}}}_h^{n+1}+\nabla {\hslash }_h^{n+1}$. Step 5. Given a relaxation parameter $\rho >0$, find $p_h^{n+1}\in W_h$ from the Richardson update $$\begin{array}{c}\hfill \begin{array}{c}\hfill (p_h^{n+1},q)=(p_h^n,q)-Pr\rho (\nabla {\hslash }_h^{n+1},\nabla q),\forall q\in W_h.\end{array}\end{array}$$ (22) From (21) and Step 4 of Algorithm 2, we obtain $(\nabla ·{\mathbf{u}}_h^{n+1},q)=(\nabla ·{\widehat{\mathbf{u}}}_h^{n+1},q)-(\nabla {\hslash }_h^{n+1},\nabla q)=0$. So the velocity obtained by Algorithm 2 satisfies the weakly divergence-free condition. Moreover, we expect to show the iterative errors between the finite element solutions to (10)–(12) and the Uzawa-type iterative solutions to Algorithm 2. For convenience, assume that ${\mathbf{E}}_h^n={\mathbf{u}}_h-{\mathbf{u}}_h^n$, ${\widehat{\mathbf{E}}}_h^n={\mathbf{u}}_h-{\widehat{\mathbf{u}}}_h^n$, ${\eta }_h^n=p_h-p_h^n$ and ${\theta }_h^n=T_h-T_h^n$. Then, we have ${\widehat{\mathbf{E}}}_h^n={\mathbf{E}}_h^n+\nabla {\hslash }_h^n$.

Firstly, we recall the convergence results of the initial guess. Note that ${\widehat{\mathbf{u}}}_h^0={\mathbf{d}}^0+\nabla g^{-1}={\mathbf{u}}_h^0$, which implies ${\mathbf{E}}_h^0={\widehat{\mathbf{E}}}_h^0$.

Lemma 1

([19]). Let $({\mathbf{u}}_h^0,p_h^0,T_h^0)\in {\mathbf{M}}^h\times W^h\times Z^h$ be the solution of Step 1 of Algorithm 1. Then, under the assumptions of Theorem 2, we have the following results

$$\begin{array}{cc}\hfill \parallel \nabla {\theta }_h^0\parallel & \le k^{-1}\overline{\Lambda }{\parallel \gamma \parallel }_{-1},\parallel {\eta }_h^0\parallel \le 2{\beta }^{-1}Pr\Lambda N^{-1}(P{r}^{-1}\Lambda +\overline{\Lambda }),\parallel \nabla {\mathbf{E}}_h^0\parallel \le \Lambda N^{-1}(P{r}^{-1}\Lambda +\overline{\Lambda }).\hfill \end{array}$$

Secondly, we show that the solution sequence generated by Algorithm 2 is bounded.

Theorem 3.

Let $\{{\mathbf{u}}_h^n,p_h^n,T_h^n\}$ be the solution sequence of Algorithm 2. Then, under the assumptions of Theorem 2, if the relaxation parameter satisfies $\rho \in (0,2(1-\overline{\Lambda }-P{r}^{-1}\Lambda ))$, the sequences $\{\parallel \nabla {\mathbf{u}}_h^n\parallel \}$, $\{\parallel \nabla {\widehat{\mathbf{u}}}_h^n\parallel \}$, $\{\parallel p_h^n\parallel \}$ and $\{\parallel \nabla T_h^n\parallel \}$ are uniformly bounded with respect to $h$.

Proof. 

Subtracting (19) from (12), we have

$$\begin{array}{c}\hfill \begin{array}{c}\hfill b_2({\widehat{\mathbf{E}}}_h^n;T_h,s)-b_2({\widehat{\mathbf{u}}}_h^n;{\theta }_h^{n+1},s)+\kappa (\nabla {\theta }_h^{n+1},\nabla s)=0.\end{array}\end{array}$$

Setting $s={\theta }_h^{n+1}$ obtains

$$\kappa \parallel \nabla {\theta }_h^{n+1}{\parallel }^2=-b_2({\widehat{\mathbf{E}}}_h^n;T_h,{\theta }_h^{n+1}).$$

According to (6) and Theorem 2, we arrive at

$$\begin{array}{c}\hfill \begin{array}{c}\hfill \parallel \nabla {\theta }_h^{n+1}\parallel \le \overline{N}{\kappa }^{-2}{\parallel \gamma \parallel }_{-1}\parallel \nabla {\widehat{\mathbf{E}}}_h^n\parallel .\end{array}\end{array}$$

(23)

Then, subtracting (20) from (10), we have

$$\begin{array}{c}\hfill \begin{array}{c}\hfill Pr(\nabla {\widehat{\mathbf{E}}}_h^{n+1},\nabla \mathbf{v})-({\eta }_h^n,\nabla ·\mathbf{v})=-b_1({\widehat{\mathbf{E}}}_h^n;{\mathbf{u}}_h,\mathbf{v})-b_1({\widehat{\mathbf{u}}}_h^n;{\widehat{\mathbf{E}}}_h^{n+1},\mathbf{v})+PrRa(\mathbf{j}{\theta }_h^{n+1},\mathbf{v}).\end{array}\end{array}$$

(24)

Choosing $\mathbf{v}={\widehat{\mathbf{E}}}_h^{n+1}$ in (24) and combining the ensuing equation with (21) lead to

$$\begin{array}{c}\hfill Pr\parallel \nabla {\widehat{\mathbf{E}}}_h^{n+1}{\parallel }^2=-(\nabla {\eta }_h^n,\nabla {\hslash }_h^{n+1})-b_1({\widehat{\mathbf{E}}}_h^n;{\mathbf{u}}_h,{\widehat{\mathbf{E}}}_h^{n+1})+PrRa(\mathbf{j}{\theta }_h^{n+1},{\widehat{\mathbf{E}}}_h^{n+1}).\end{array}$$

Next, according to (22), we have

$$\begin{array}{c}\hfill \begin{array}{c}\hfill Pr\parallel \nabla {\widehat{\mathbf{E}}}_h^{n+1}{\parallel }^2={(Pr\rho )}^{-1}(p_h^{n+1}-p_h^n,{\eta }_h^n)-b_1({\widehat{\mathbf{E}}}_h^n;{\mathbf{u}}_h,{\widehat{\mathbf{E}}}_h^{n+1})+PrRa(\mathbf{j}{\theta }_h^{n+1},{\widehat{\mathbf{E}}}_h^{n+1}),\end{array}\end{array}$$

which, by using (5), (6), (23), Theorem 2, and the Proposition identity $(u,v)=\frac{1}{2}{(\parallel u+v\parallel }^2-{\parallel u\parallel }^2-{\parallel v\parallel }^2)$, we have

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill 2P{r}^2\rho \parallel \nabla {\widehat{\mathbf{E}}}_h^{n+1}{\parallel }^2+{\parallel {\eta }_h^{n+1}\parallel }^2\le & \parallel {\eta }_h^n{\parallel }^2+{\parallel {\eta }_h^{n+1}-{\eta }_h^n\parallel }^2\hfill \\ & +2Pr\rho (\Lambda +Pr\overline{\Lambda })\parallel \nabla {\widehat{\mathbf{E}}}_h^n\parallel \parallel \nabla {\widehat{\mathbf{E}}}_h^{n+1}\parallel .\hfill \end{array}\end{array}$$

(25)

Then, using (21) and (22), we obtain

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill \parallel {\eta }_h^{n+1}-{\eta }_h^n{\parallel }^2& =(p_h^{n+1}-p_h^n,p_h^{n+1}-p_h^n)=-Pr\rho (\nabla {\hslash }_h^{n+1},\nabla ({\eta }_h^{n+1}-{\eta }_h^n))\hfill \\ & =Pr\rho (\nabla ·{\widehat{\mathbf{E}}}_h^{n+1},{\eta }_h^{n+1}-{\eta }_h^n),\hfill \end{array}\end{array}$$

which leads to

$$\begin{array}{c}\hfill \begin{array}{c}\hfill \parallel {\eta }_h^{n+1}-{\eta }_h^n{\parallel }^2\le {(Pr\rho )}^2\parallel \nabla ·{\widehat{\mathbf{E}}}_h^{n+1}{\parallel }^2\le {(Pr\rho )}^2{\parallel \nabla {\widehat{\mathbf{E}}}_h^{n+1}\parallel }^2,\end{array}\end{array}$$

(26)

where we have applied the fact that $\parallel \nabla ·\mathbf{v}\parallel \le \parallel \nabla \mathbf{v}\parallel$ in [24].

Moreover, substituting (26) into (25) and using the Young inequality, we obtain

$$\begin{array}{c}\hfill \begin{array}{cc}& \parallel \nabla {\widehat{\mathbf{E}}}_h^{n+1}{\parallel }^2(2P{r}^2\rho -P{r}^2{\rho }^2-\varsigma (Pr\rho \Lambda +P{r}^2\rho \overline{\Lambda }))+{\parallel {\eta }_h^{n+1}\parallel }^2\hfill \\ & \le \parallel {\eta }_h^n{\parallel }^2+{\varsigma }^{-1}(Pr\rho \Lambda +P{r}^2\rho \overline{\Lambda }){\parallel \nabla {\widehat{\mathbf{E}}}_h^n\parallel }^2,\hfill \end{array}\end{array}$$

(27)

where $\varsigma >0$ is a parameter to be determined later on.

Furthermore, we solve a quadratic algebraic equation

$${\varsigma }^2(\Lambda +Pr\overline{\Lambda })-\varsigma (2Pr-Pr\rho )+(\Lambda +Pr\overline{\Lambda })=0,$$

to obtain a positive root $\varsigma ={\varsigma }^{*}$, which makes $(2Pr-Pr\rho -\varsigma (\Lambda +Pr\overline{\Lambda }))={\varsigma }^{-1}(\Lambda +Pr\overline{\Lambda })$ hold. In fact, we have

$$\varsigma ={\varsigma }^{*}=\frac{(2Pr-Pr\rho )-\sqrt{\Delta }}{2(\Lambda +Pr\overline{\Lambda })},$$

where $\Delta :=(2Pr-Pr\rho +2(\Lambda +Pr\overline{\Lambda }))(2Pr-Pr\rho -2(\Lambda +Pr\overline{\Lambda }))$.

Next, we set

$$\begin{array}{c}\hfill D_1=Pr\rho (2Pr-Pr\rho -{\varsigma }^{*}(\Lambda +Pr\overline{\Lambda }))=Pr\rho (\Lambda +Pr\overline{\Lambda })/{\varsigma }^{*}=\frac{P{r}^2\rho (2-\rho )+\sqrt{\Delta }}{2}.\end{array}$$

Thus, the inequality (27) is rewritten as

$$D_1\parallel \nabla {\widehat{\mathbf{E}}}_h^{n+1}{\parallel }^2+\parallel {\eta }_h^{n+1}{\parallel }^2\le \parallel {\eta }_h^n{\parallel }^2+D_1{\parallel \nabla {\widehat{\mathbf{E}}}_h^n\parallel }^2,$$

which, along with (23), implies that

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill D_1\parallel \nabla {\widehat{\mathbf{E}}}_h^{n+1}{\parallel }^2+{\parallel {\eta }_h^{n+1}\parallel }^2& \le \parallel {\eta }_h^0{\parallel }^2+D_1{\parallel \nabla {\widehat{\mathbf{E}}}_h^0\parallel }^2,\hfill \\ \hfill \parallel \nabla {\theta }_h^{n+1}\parallel & \le {\overline{N}}^2{\kappa }^{-4}{\parallel \gamma \parallel }_{-1}^2(\parallel {\eta }_h^0{\parallel }^2+D_1\parallel \nabla {\widehat{\mathbf{E}}}_h^0{\parallel }^2).\hfill \end{array}\end{array}$$

(28)

Finally, applying (26) into (22), we obtain

$$\begin{array}{c}\hfill \begin{array}{c}\hfill \parallel \nabla {\hslash }_h^{n+1}\parallel \le C_p^2{(Pr\rho )}^{-1}\parallel p_h^n-p_h^{n+1}\parallel \le C_p^2{(Pr\rho )}^{-1}\parallel {\eta }_h^{n+1}-{\eta }_h^n\parallel \le C_p^2\parallel \nabla {\widehat{\mathbf{E}}}_h^{n+1}\parallel ,\end{array}\end{array}$$

which combines with ${\widehat{\mathbf{E}}}_h^{n+1}={\mathbf{E}}_h^{n+1}+\nabla {\hslash }_h^{n+1}$; then, we have

$$\begin{array}{c}\hfill \begin{array}{c}\hfill \parallel {\mathbf{E}}_h^{n+1}{\parallel }^2\le 2(\parallel {\widehat{\mathbf{E}}}_h^{n+1}{\parallel }^2+\parallel \nabla {\hslash }_h^{n+1}{\parallel }^2)\le 4C_p^4{\parallel \nabla {\widehat{\mathbf{E}}}_h^{n+1}\parallel }^2,\end{array}\end{array}$$

(29)

Finally, combining (29) with (28), we obtain

$$\begin{array}{c}\hfill \begin{array}{c}\hfill D_1\parallel {\mathbf{E}}_h^{n+1}{\parallel }^2\le 4C_p^4(\parallel {\eta }_h^0{\parallel }^2+D_1\parallel \nabla {\widehat{\mathbf{E}}}_h^0{\parallel }^2).\end{array}\end{array}$$

(30)

Hence, using (28), (30), and Lemma 1, we finish the proof of the theorem. □

Thirdly, we are going to develop the convergence analysis for Algorithm 2.

Theorem 4.

Under the assumptions of Theorem 3, the following estimates hold

$$\begin{array}{cc}\hfill P{r}^{2}D{\parallel {\mathbf{E}}_h^{n+1}\parallel }^2& \le 4C_p^4{H}^{n+1}(P{r}^{2}D\parallel \nabla {\widehat{\mathbf{E}}}_h^0{\parallel }^2+\parallel {\eta }_h^0{\parallel }^2),\parallel {\eta }_h^{n+1}{\parallel }_0^2\le H^{n+1}(P{r}^{2}D\parallel \nabla {\widehat{\mathbf{E}}}_h^0{\parallel }^2+\parallel {\eta }_h^0{\parallel }^2),\hfill \\ \hfill P{r}^{2}D{\parallel \nabla {\theta }_h^{n+1}\parallel }^2& \le {\overline{N}}^2{\kappa }^{-4}{\parallel \gamma \parallel }_{-1}^2{H}^n(P{r}^{2}D\parallel \nabla {\widehat{\mathbf{E}}}_h^0{\parallel }^2+\parallel {\eta }_h^0{\parallel }^2),\hfill \end{array}$$

where $D\in (0,\frac{1}{2})$ and $H\in (\frac{3}{4},1)$ are two constants independent of n and h.

Proof. 

By Theorem 3, there exists a positive constant $D_2$, independent of n and h, such that

$$\begin{array}{c}\hfill \parallel \nabla {\widehat{\mathbf{u}}}_h^n\parallel \le D_2.\end{array}$$

(31)

Then, rewrite (24) to obtain

$$\begin{array}{c}\hfill ({\eta }_h^n,\nabla ·\mathbf{v})=Pr(\nabla {\widehat{\mathbf{E}}}_h^{n+1},\nabla \mathbf{v})+b_1({\widehat{\mathbf{E}}}_h^n;{\mathbf{u}}_h,\mathbf{v})+b_1({\widehat{\mathbf{u}}}_h^n;{\widehat{\mathbf{E}}}_h^{n+1},\mathbf{v})-PrRa(\mathbf{j}{\theta }_h^{n+1},\mathbf{v}).\end{array}$$

Applying the inf-sup condition, (5), (6), (23), and Theorem 2 to the above equation, we obtain

$$\begin{array}{cc}\hfill \beta \parallel {\eta }_h^n\parallel \le & Pr\parallel \nabla {\widehat{\mathbf{E}}}_h^{n+1}\parallel +PrRa{C}_p^2\overline{N}{\kappa }^{-2}{\parallel \gamma \parallel }_{-1}\parallel \nabla {\widehat{\mathbf{E}}}_h^n\parallel +PrRa{C}_p^{2}N{\kappa }^{-1}{\parallel \gamma \parallel }_{-1}\parallel \nabla {\widehat{\mathbf{E}}}_h^n\parallel \hfill \\ \hfill & +N\parallel \nabla {\widehat{\mathbf{u}}}_h^n\parallel \parallel \nabla {\widehat{\mathbf{E}}}_h^{n+}\parallel ,\hfill \end{array}$$

which combines with (31) to obtain

$$\begin{array}{c}\hfill \beta \parallel {\eta }_h^n\parallel \le (Pr+N{D}_2)\parallel \nabla {\widehat{\mathbf{E}}}_h^{n+1}\parallel +(\Lambda +Pr\overline{\Lambda })\parallel \nabla {\widehat{\mathbf{E}}}_h^n\parallel .\end{array}$$

Next, using the inequality ${(a+b)}^2\le 2a^2+2b^2$, we have

$$\begin{array}{c}\hfill {\beta }^2\parallel {\eta }_h^n{\parallel }^2\le 2{(Pr+N{D}_2)}^2\parallel \nabla {\widehat{\mathbf{E}}}_h^{n+1}{\parallel }^2+2{(\Lambda +Pr\overline{\Lambda })}^2{\parallel \nabla {\widehat{\mathbf{E}}}_h^n\parallel }^2.\end{array}$$

Hence, one obtains

$$\begin{array}{c}\hfill \begin{array}{c}\hfill \parallel \nabla {\widehat{\mathbf{E}}}_h^{n+1}{\parallel }^2\ge D_3\parallel {\eta }_h^n{\parallel }^2-D_4{\parallel \nabla {\widehat{\mathbf{E}}}_h^n\parallel }^2,\end{array}\end{array}$$

(32)

where $D_3:=\frac{{\beta }^2}{2{(Pr+N{D}_2)}^2}$ and $D_4:=\frac{{(\Lambda +Pr\overline{\Lambda })}^2}{{(Pr+N{D}_2)}^2}$. Obviously, if we let $C_{\rho ,\varsigma }:=Pr\rho (2Pr-Pr\rho -\varsigma (\Lambda +Pr\overline{\Lambda }))$, then (27) becomes

$$\begin{array}{c}\hfill \begin{array}{c}\hfill \delta \parallel \nabla {\widehat{\mathbf{E}}}_h^{n+1}{\parallel }^2+(C_{\rho ,\varsigma }-\delta )\parallel \nabla {\widehat{\mathbf{E}}}_h^{n+1}{\parallel }^2+\parallel {\eta }_h^{n+1}{\parallel }^2\le \parallel {\eta }_h^n{\parallel }^2+{\varsigma }^{-1}(Pr\rho \Lambda +P{r}^2\rho \overline{\Lambda }){\parallel \nabla {\widehat{\mathbf{E}}}_h^n\parallel }^2.\end{array}\end{array}$$

(33)

where $\delta \in (0,C_{\rho ,\varsigma })$ is a parameter to be determined. From (32) and (33), we obtain

$$\begin{array}{c}\hfill \begin{array}{cc}& (C_{\rho ,\varsigma }-\delta )\parallel \nabla {\widehat{\mathbf{E}}}_h^{n+1}{\parallel }^2+{\parallel {\eta }_h^{n+1}\parallel }^2\hfill \\ & \le (1-D_3\delta )\parallel {\eta }_h^n{\parallel }^2+({\varsigma }^{-1}(Pr\rho \Lambda +P{r}^2\rho \overline{\Lambda })+D_4\delta ){\parallel \nabla {\widehat{\mathbf{E}}}_h^n\parallel }^2.\hfill \end{array}\end{array}$$

(34)

Then, we will choose parameters $\varsigma$ and $\delta$ such that

$$\begin{array}{c}\hfill \begin{array}{c}\hfill \frac{C_{\rho ,\varsigma }-\delta }{1}=\frac{{\varsigma }^{-1}Pr\rho (\Lambda +Pr\overline{\Lambda })+D_4\delta }{1-\delta D_3},\end{array}\end{array}$$

(35)

and $1-\delta D_3>0$, which leads to

$$\begin{array}{c}\hfill \begin{array}{c}\hfill D_3{\delta }^2-(1+C_{\rho ,\varsigma }{D}_3+D_4)\delta +C_{\rho ,\varsigma }-{\varsigma }^{-1}Pr\rho (\Lambda +Pr\overline{\Lambda })=0.\end{array}\end{array}$$

(36)

In fact, one finds that

$$C_{\rho ,\varsigma }-{\varsigma }^{-1}Pr\rho (\Lambda +Pr\overline{\Lambda })=(1+C_{\rho ,\varsigma }{D}_3+D_4)\delta -D_3{\delta }^2>C_{\rho ,\varsigma }{D}_3\delta -D_3{\delta }^2>0,$$

which, along with the definition of $C_{\rho ,\varsigma }$, yields

$$(\Lambda +Pr\overline{\Lambda }){\varsigma }^2-(2Pr-Pr\rho )\varsigma +(\Lambda +Pr\overline{\Lambda })<0,$$

and

$$\frac{(2Pr-Pr\rho )-\sqrt{\Delta }}{2(\Lambda +Pr\overline{\Lambda })}<\varsigma <\frac{(2Pr-Pr\rho )+\sqrt{\Delta }}{2(\Lambda +Pr\overline{\Lambda })},$$

where the notation $\Delta$ is defined in the proof of Theorem 3. Note that we have used condition $0<\rho <2(1-\overline{\Lambda }-P{r}^{-1}\Lambda )$. Here, we select

$$\varsigma ={\varsigma }^{+}=\frac{2Pr-Pr\rho }{2(\Lambda +Pr\overline{\Lambda })}.$$

Substituting this parameter into (36), we arrive at $a{\delta }^2-b\delta +c=0$, where $a=D_3$, $b=1+D_4+s_{1}a$, $c=s_1-\frac{P{r}^2{\rho }^2{(\Lambda +Pr\overline{\Lambda })}^2}{s_1}$, and $s_1=P{r}^2\rho (1-\frac{1}{2}\rho )$. Obviously, $b>1+s_{1}a$, $c<s_1$; so, we deduce that

$$b^2-4ac>{(1+s_{1}a)}^2-4a{s}_1\ge 0.$$

Then, the Equation (36) has a real root ${\delta }^{*}=\frac{b-\sqrt{b^2-4ac}}{2a}$.

With the parameter $\epsilon$ and $\delta$ given by ${\epsilon }^{+}$ and ${\delta }^{*}$, it follows from (34) that

$$\begin{array}{c}\hfill \begin{array}{c}\hfill \overline{D}\parallel \nabla {\widehat{\mathbf{E}}}_h^{n+1}{\parallel }^2+\parallel {\eta }_h^{n+1}{\parallel }_0^2\le H(\overline{D}\parallel \nabla {\widehat{\mathbf{E}}}_h^n{\parallel }^2+\parallel {\eta }_h^n{\parallel }^2),\end{array}\end{array}$$

(37)

where $\overline{D}=s_1-{\delta }^{*}$ and $H=1-{\delta }^{*}{D}_3$.

Note that $\overline{D}>0$ and $H>0$. Now, we will prove them. Consider the quadratic function $f(\delta )=a{\delta }^2-b\delta +c$. Because $a>0$, $s_1>0$, $b>1+s_{1}a$ and $c<s_1$, we obtain ${lim}_{\delta \to -\infty }f(\delta )=\infty$ and

$$f(s)=a{s}_1^2-b{s}_1+c<a{s}_1^2-(1+a{s}_1)s_1+s_1=0.$$

Thus, the smallest root ${\delta }^{*}$ of $f(\delta )$ must belong to $(-\infty ,s_1)$. So, the inequality $\overline{D}>0$ holds. Noticing that $C_{\rho ,{\varsigma }^{+}}-{\delta }^{*}=s_1-{\delta }^{*}>0$, it follows readily from (35) that $H>0$.

Finally, note that $0<\overline{D}<s_1=P{r}^2\rho -\frac{1}{2}P{r}^2{\rho }^2\le \frac{P{r}^2}{2}$. If, we choose the $\overline{D}=P{r}^{2}D$ and $0<D<\frac{1}{2}$, the inequality (37) is rewritten as

$$\begin{array}{c}\hfill P{r}^{2}D\parallel \nabla {\widehat{\mathbf{E}}}_h^{n+1}{\parallel }^2+\parallel {\eta }_h^{n+1}{\parallel }_0^2\le H_1(P{r}^{2}D\parallel \nabla {\widehat{\mathbf{E}}}_h^n{\parallel }^2+\parallel {\eta }_h^n{\parallel }^2).\end{array}$$

(38)

According to the definition of $D_3$ and $\beta \le 1$, we arrive at $D_3\le \frac{1}{2P{r}^2}$. Noticing that ${\delta }^{*}<s_1<\frac{P{r}^2}{2}$, we easily find that $1>H=1-{\delta }^{*}{D}_3>\frac{3}{4}$.

Next, using (38) and (29), we obtain

$$\begin{array}{c}\hfill P{r}^{2}D\parallel {\mathbf{E}}_h^{n+1}{\parallel }^2\le 4C_p^4{H}^{n+1}(P{r}^{2}D\parallel \nabla {\widehat{\mathbf{E}}}_h^0{\parallel }^2+\parallel {\eta }_h^0{\parallel }^2),\parallel {\eta }_h^{n+1}{\parallel }_0^2\le H^{n+1}(P{r}^{2}D\parallel \nabla {\widehat{\mathbf{E}}}_h^0{\parallel }^2+\parallel {\eta }_h^0{\parallel }^2).\end{array}$$

Finally, using the above estimates with (23), we finish the proof. □

4. Numerical Study

We will represent some numerical tests to claim the accuracy and performance of the proposed algorithm for the steady natural convection problem in this section. We used the public finite element software FreeFem++ [27] and applied $P_2-P_1-P_2$ element to approximate the velocity, temperature, and pressure, respectively.

In the first numerical test, let the domain $\Omega =[0,1]\times [0,1]$, and the right-hand side of (1)–(4) is selected such that the exact solutions are given by

$$\begin{array}{cc}\hfill p(x,y)& =cos(\pi x)cos(\pi y),T(x,y)=u_1(x,y)+u_2(x,y)\hfill \\ \hfill u_1(x,y)& =2\pi {sin}^2(\pi x)sin(\pi y)cos(\pi y),u_2(x,y)=-2\pi sin(\pi x){sin}^2(\pi y)cos(\pi x).\hfill \end{array}$$

Here, we set the parameters $Ra=Pr=\kappa =1$ and use the stopping rule

$$max\left\{\frac{\parallel {\mathbf{u}}_h^{n+1}-{\mathbf{u}}_h^n\parallel }{\parallel {\mathbf{u}}_h^n\parallel },\frac{\parallel p_h^{n+1}-p_h^n\parallel }{\parallel p_h^n\parallel },\frac{\parallel T_h^{n+1}-T_h^n\parallel }{\parallel T_h^n\parallel }\right\}<1.0\times {10}^{-6}.$$

Figure 1 displays the iteration errors of the velocity, temperature in $H^1$-seminorm, and the pressure in $L^2$-norm for different iterative steps n solved by Algorithm 2. Here, we set the relaxation parameter $\rho =1.6$ and choose five different mesh sizes h. From Figure 1, we observe that the proposed algorithm worked well and kept the convergence when iteration step n became large.

In the above test, we fixed the relaxation parameter and varied the mesh size. Now, we consider different relaxation parameters with the mesh size $h=\frac{1}{32}$. Figure 2 expresses different iterative steps of the log errors with different values $\rho$. From Figure 2, we observe that ${\mathbf{u}}_h^n$, $p_h^n$, and $T_h^n$ converged faster when $\rho$ was larger. However, we have an interesting observation that it became slow when $\rho$ was too large (e.g., $\rho =1.7$ or 1.9). It is not surprising since from Theorem 3 and 4 the relaxation parameter $\rho$ had a limited interval, and the value $\rho =1.7$ or 1.9 may have been out of its interval.

Hence, we should reveal the convergence on the relaxation parameter $\rho$ by showing the values with respect to n and $\rho$ under the mesh size $h=\frac{1}{32}$. From

Table 1. The iterative step n with the relaxation parameter $\rho$.

$\mathit{\rho }$	0.1	0.2	0.3	0.4	0.5	0.6	0.7	0.8	0.9	1.0	1.1	1.2	1.3	1.4	1.5	1.6	1.7	1.8	1.9	2.0
Algorithm 1	509	280	197	153	126	107	93	83	74	67	62	57	53	50	47	44	49	76	159	/
Algorithm 2	531	289	202	156	127	108	94	83	74	67	61	56	52	44	48	42	50	77	154	/

The mark “/” means that the iterative step was larger than 600.

, we find that Algorithms 1 and 2 converged faster when we chose larger $\rho$. However, if the $\rho$ chosen was very large, then these algorithms either need more iterative steps or diverge. In addition, Algorithms 1 and 2 achieved the tolerance error when $\rho =1.6$ with the least iterative steps $n=44$ and $n=42$, respectively.

Based on the previous section, Algorithm 2 produced the divergence-free velocity approximation. Hence, in

Table 2. The value of $\parallel \nabla ·{\mathbf{u}}_h^n\parallel$ with different Rayleigh numbers $Ra$.

$\mathit{R}\mathit{a}$	10	100	150	180
Algorithm 2	1.82 × 10${}^{-8}$	2.65 × 10${}^{-10}$	2.02 × 10${}^{-11}$	4.96 × 10${}^{-12}$
Algorithm 1	3.50 × 10${}^{-18}$	/	/	/

The mark “/” means that the iterative step was larger than 600.

we list the value of $\parallel \nabla ·{\mathbf{u}}_h^n\parallel$. From this table, Algorithms 1 and 2 obtain good numerical results when $Ra=10$. However, when the value of $Ra$ increased, then Algorithm 1 could not achieve the tolerance error and converge. Meanwhile, Algorithm 2 still ran well.

In the second numerical test, we considered the hot cylinder problem solving the proposed algorithm with different Rayleigh numbers. The boundary conditions are given in [28,29], i.e., $\frac{\partial T}{\partial n}=1$ on inner wall, $T=0$ on the other wall, and zero Dirichlet condition on velocity were imposed. Set $Pr=0.7,\kappa =1$, $\gamma =0$, and $h=\frac{1}{80}$. Figure 3 and Figure 4 express the numerical streamlines, isobars, and isotherms for different radii of inner circle $r_{in}$ based on $Ra=100$ and $Ra=250$ with $\rho =1.6$. We observe that it shapes two vortices when $r_{in}=0.2$ and four vortices when $r_{in}=0.8$, which were found to be in good agreement with those reported in [28,29]. Therefore, the given method captured this classical model well.

In

Table 3. Comparisons of numerical results from different algorithms with $h=\frac{1}{80},r_{in}=0.2$.

		Ra = 100			Ra = 250
	x = 0.5	y = 0.5	CPU Time	x = 0.5	y = 0.5	CPU Time
Algorithm 2	0.281	0.284	14.135	0.755	0.760	22.135
Algorithm 1 [19]	0.263	0.465	33.772	/	/	/
Wang’s algorithm [29]	0.274	0.279	51.890	0.714	0.722	56.571

The mark “/” means that the iterative step was larger than 600.

and

Table 4. Comparisons of numerical results from different algorithms with $h=\frac{1}{80},r_{in}=0.8$.

		Ra = 100			Ra = 250
	x = 0.5	y = 0.5	CPU Time	x = 0.5	y = 0.5	CPU Time
Algorithm 2	0.039	0.085	1.811	0.098	0.213	2.191
Algorithm 1 [19]	0.039	0.085	2.077	/	/	/
Wang’s algorithm [29]	0.039	0.086	8.851	0.098	0.214	9.169

The mark “/” means that the iterative step was larger than 600.

, we show the CPU time and the maximum value of velocity at $x=0.5$ and $y=0.5$ by Algorithms 1 and 2 with $\rho =1.6$ and Wang’s algorithm [29] for $r_{in}=0.2$ and $r_{in}=0.8$, respectively. From Table 3 and Table 4, we find that the proposed algorithm took the least computational time among these algorithms to obtain almost the same maximum value of velocity. In particular, Algorithm 1 did not work when $Ra=250$. Therefore, the proposed algorithm solved this model well.

5. Conclusions

In conclusion, we designed a Uzawa-type iterative algorithm based on the mixed finite element method to solve the stationary natural convection model. Compared with the common Uzawa iterative algorithm, a central feature of the proposed algorithm is that it produced weakly divergence-free velocity approximation. This algorithm can be extended to the double-diffusive natural convection [30] and the magnetohydrodynamics flows [31].