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

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.


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 saddlepoint 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 algo-rithms. 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 divergencefree condition.

Preliminaries
Let Ω ⊂ R 2 be a bounded domain, which has a Lipschitz continuous boundary ∂Ω with a regular open subset Γ. Consider the following stationary natural convection problem.
Seek the velocity u = (u 1 (x), u 2 (x)) , the pressure p = p(x), and the temperature T(x), such that where γ is the forcing function, n is the outward unit vector, and j = (0, 1) . In addition, the positive parameter κ 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: Here, the space L 2 (Ω) is endowed with L 2 -scalar product (·, ·) and L 2 -norm · . In addition, the space H 1 (Ω) is used to represent the standard definitions for Sobolev spaces W m,p (Ω), m, p > 0.
Moreover, we recall the Poincaré inequality [21] as follows: where C p is the Poincaré constant. Next, we denote two trilinear forms by which satisfy the following properties [7,22,23] for all u, v, w, ∈ M and T, s ∈ Z. Here, N and N are two fixed positive constants.
The following existence and uniqueness of the solution to (6) are classical results.
Next, we consider a family of quasi-uniform and regular triangulations K h = {K : ∪ K⊂Ω K = Ω} with mesh size h, which is a partition of the domain Ω. Then, we assume that the finite element subspace 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 M h × W h satisfy the following discrete inf-sup condition [21]; for each where the constant β ∈ (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 (u h , p h , The following theorem is established for the stability of the finite element discretization.

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 (∇ · u n+1 h , q) = 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 (∇ · 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 2. Given a relaxation parameter ρ > 0, find In order to make the velocity of Uzawa algorithm have a weakly divergence-free property, let g be a gauge variable [26] and d be a variable, such that u = d + ∇g. If g and p satisfy an elliptic equation Pr∆g = p, then (1)-(4) can be rewritten as for n = 0, 1, . . .
whereh n+1 := g n+1 − g n . So one obtains p n+1 = P r ∆g n+1 = P r ∆h n+1 + Pr∆g n = P r ∆h n+1 + p n , and Now, we are ready to write the Uzawa-type finite element iterative algorithm as follows Algorithm 2.
Step 1. Obtain the initial guess (u 0 Step Step From (21) and Step 4 of Algorithm 2, we obtain 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 Firstly, we recall the convergence results of the initial guess. Note thatû 0 Secondly, we show that the solution sequence generated by Algorithm 2 is bounded.
Thirdly, we are going to develop the convergence analysis for Algorithm 2.

Theorem 4.
Under the assumptions of Theorem 3, the following estimates hold where D ∈ (0, 1 2 ) and H ∈ ( 3 4 , 1) are two constants independent of n and h.

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.
Here, we set the parameters Ra = Pr = κ = 1 and use the stopping rule   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 = 1 32 . Figure 2 expresses different iterative steps of the log errors with different values ρ. From Figure 2, we observe that u n h , p n h , and T n h converged faster when ρ was larger. However, we have an interesting observation that it became slow when ρ was too large (e.g., ρ = 1.7 or 1.9). It is not surprising since from Theorem 3 and 4 the relaxation parameter ρ had a limited interval, and the value ρ = 1.7 or 1.9 may have been out of its interval.   Hence, we should reveal the convergence on the relaxation parameter ρ by showing the values with respect to n and ρ under the mesh size h = 1 32 . From Table 1, we find that Algorithms 1 and 2 converged faster when we chose larger ρ. However, if the ρ 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 ρ = 1.6 with the least iterative steps n = 44 and n = 42, respectively. The mark "/" means that the iterative step was larger than 600.
Based on the previous section, Algorithm 2 produced the divergence-free velocity approximation. Hence, in Table 2 we list the value of ∇ · u n h . 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. The mark "/" means that the iterative step was larger than 600.
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., ∂T ∂n = 1 on inner wall, T = 0 on the other wall, and zero Dirichlet condition on velocity were imposed. Set Pr = 0.7, κ = 1, γ = 0, and h = 1 80 . Figures 3 and 4 express the numerical streamlines, isobars, and isotherms for different radii of inner circle r in based on Ra = 100 and Ra = 250 with ρ = 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 Tables 3 and 4, 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 ρ = 1.6 and Wang's algorithm [29] for r in = 0.2 and r in = 0.8, respectively. From Tables 3 and 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.   Figure 4. Numerical streamlines (the first column), isotherms (the second column), and isobars (the third column) for Ra = 100 (the first line) and Ra = 250 (the second line) with r in = 0.8. The mark "/" means that the iterative step was larger than 600. The mark "/" means that the iterative step was larger than 600.

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].