Finite Element Iterative Methods for the Stationary Double-Diffusive Natural Convection Model

In this paper, we consider the stationary double-diffusive natural convection model, which can model heat and mass transfer phenomena. Based on the fixed point theorem, the existence and uniqueness of the considered model are proved. Moreover, we design three finite element iterative methods for the considered problem. Under the uniqueness condition of a weak solution, iterative method I is stable. Compared with iterative method I, iterative method II is stable with a stronger condition. Moreover, iterative method III is stable with the strongest condition. From the perspective of viscosity, iterative method I displays well in the case of a low viscosity number, iterative method II runs well with slightly low viscosity, and iterative method III can deal with high viscosity. Finally, some numerical experiments are presented for testing the correctness of the theoretic analysis.


Introduction
The double-diffusive natural convection model, which does not only incorporate the velocity vector field as well as the pressure field, but also contains the temperature field and the concentration field, has been widely used in scientific, engineering and industrial applications such as nuclear design, cooling of electronic equipment, aircraft cabins, insulation with double pane windows, and so on. For greater understanding of the physical background, authors can refer to [1][2][3]. In recent years, the impact of nanofluid on free convection heat transfer was investigated by researchers in [4]. The free convective flow of a Nano-Encapsulated Phase Change Material (NEPCM) suspension in an eccentric annulus was investigated numerically in [5]. The authors obtained that the volume fraction of the NEPCM particles and Stefan number effect the thermal and hydrodynamic characteristics of the suspension. The effect of the arrangement of the tubes in a multi-tube heat exchanger during the solidification process was considered in [6], which focused on the natural convection effect in phase change material in this research.
Let Ω ⊂ R 2 be a open bounded domain with a Lipschitz continuous boundary ∂Ω and ∂Γ is a subset of ∂Ω, u = (u 1 , u 2 ) denotes the velocity field, p is the fluid pressure, T is the temperature, C is the concentration, g = (0, 1) is the gravitational acceleration vector, f i is the forcing function, i = 1, 2. Moreover, n represents the outer normal vector, ν > 0 is the viscosity, D a is the Darcy number, γ > 0 is the heat diffusivity, D c is the mass diffusivity, β T and β C are the thermal and solutal expansion coefficients.
The governing equations of this double-diffusive natural convection model are presented as follows [7].
in Ω, −ν ∂ 2 u 2 ∂x 2 + ∂ 2 u 2 ∂y 2 + u 1 ∂u 2 ∂x + u 2 ∂u 2 ∂y = β T T + β C C − D −1 a u 2 − ∂p ∂y , in Ω, in Ω, Many numerical studies were made concerning the double-diffusive natural convection model. A projection-based stabilized finite element method for steady-state natural convection problem was considered in [8]. A stabilized finite element error analysis for the Darcy-Brinkman model of double-diffusive convection in a porous medium was discussed in [9]. An efficient two-step algorithm for the steady-state natural convection problem was presented in [10]. The melting process of a nano-enhanced phase change material in a square cavity was investigated in [11]. In numerical test, the author used the Galerkin finite element method to solve the dimensionless partial differential equations. Based on the idea of curvature stabilization, Çıbık et al. [12] discussed a family of second order time stepping methods for the Darcy-Brinkman equations. A decoupled finite element method called the modified characteristics method was considered in [13]. Rajabi et al. performed the detailed uncertainty propagation analysis and variance-based global sensitivity analysis on the widely adopted double-diffuse convection benchmark problem of a square porous cavity with horizontal temperature and concentration gradients in [14]. In [15], the mixed convection heat transfer of AL2O3 nanofuid in a horizontal channel subjected with two heat sources was considered. In [16], the curvature-based stabilization method was considered for double-diffusive natural convection flows in the presence of a magnetic field and unconditionally stable and optimally accurate second order approximations were obtained. There are several works devoted to the efficient numerical methods for the treatment of nonlinear problems. For example, several iterative methods for the 2D steady penalty Navier-Stokes equations were presented and discussed in [17]. He et al. [18] discussed a combination of two-level methods and iterative methods for solving the 2D/3D steady Navier-Stokes equations. Some iterative finite element methods for steady Navier-Stokes equations with different viscosities were discussed in [19]. Furthermore, the authors refer to the Oseen method [20], the Newton method [21] and the Euler implicit-explicit methods [22]. Recently, Huang et al. [23] have considered and analyzed the Oseen, Newton and Stokes iterative methods for the 2D steady Navier-Stokes equations. He et al. [24] considered and analyzed three iterative methods for the 3D steady MHD equations.
The main work in this paper is to design, analyze, and compare three iteration methods to solve nonlinear equations based on the finite element discretization. Then, we will show the performance of these numerical methods in both theoretical analysis and numerical experiments. By setting , we obtain the conclusion that the three iterative methods are stable and convergent as σ ∈ (0, 1 4 ). Iterative method I and II are valid as σ ∈ [ 1 4 , 1 3 ) and only iterative method I runs well as σ ∈ [ 1 3 , 1). In this paper, by developing some techniques and using some ideas in [7], we prove the existence and uniqueness with a different method, then we obtain a different uniqueness condition. Furthermore, we propose and analyze iterative methods I and III. In addition to this, we use iterative method II to computer a smaller viscosity than them in numerical experiments. Compared with He et al. [24], although the iterative methods are the same, the considered problems are different.
The paper is organized as follows. In Section 2, we describe the considered problem and some mathematical preliminaries. In the next section, we prove the existence and uniqueness of the weak solution to the considered equations. Then, we analyze stability and iterative error estimates of three iterative methods in Section 4. In Section 5, we show some numerical experiments to verify the correctness of theoretical results. In the last section, conclusions are presented.
The variational form of the model (1) is presented as follows:

Existence and Uniqueness
This section gives the existence and uniqueness of (5), which is crucial to consider the discrete scheme. Theorem 1. There exists at least a solution pair (u, p, T, C) ∈ X × M × W × Q which satisfies (5) and Proof. First, for u ∈ X, it is easy to see that a 1 (·, ·) + c 1 (u, ·, ·) and a 2 (·, ·) + c 2 (u, ·, ·) are continuous, elliptic bilinear forms of W × W and Q × Q, respectively. Hence, according to the Lax-Milgram theorem, there exists a unique solution T ∈ W to the second equation of (5), and a unique solution C ∈ Q to the third equation of (5). The theorem will be proved if we can show that there is at least a solution u ∈ X in the first equation of (5). Secondly, a 0 (·, ·) is a continuous and elliptic bilinear form on X × X. Using (2) and (4) we obtain where m = |g| max{|β T |, |β C |}. Then, we define a mapping A : Clearly, u is a solution of the first equation of (5) with v ∈ V, if it is a solution of A(u) = u. Using the Leray-Schauder Principle [26], A(u) = u has at least one solution u ∈ X, if (a) A is completely continuous; (b) there exists M 1 > 0 such that for every λ ∈ [0, 1] and v ∈ X with λAv = v, v satisfies the bound ∇v 0 ≤ M 1 .
Assume u 1 , u 2 ∈ X and subtract the equations obtained from (7) with u = u 1 and u = u 2 . Then, set w = A(u 2 ) − A(u 1 ) and choose v = w to yield Moreover, in order to estimate ∇(T 2 − T 1 ) 0 and ∇(C 2 − C 1 ) 0 , we substitute T 1 and T 2 in the second equation of (5) and subtract the ensuing equations to obtain Taking ψ = T 2 − T 1 and using (3) we obtain Analogously, we have Further, combining (9) and (10), we obtain the bound of (8) as follows and v ∈ X satisfies λAv = v. Then, from (7), we have Using (2) and (4), we arrive at Thirdly, setting ψ = T in the second equation of (5), we have Similarly, taking s = C in the third equation of (5), we obtain Moreover, choosing v = u in the first equation of (5) and using (4), we arrive at which combines with the above two equations to give The proof is completed.

Theorem 2.
Assume that (u, p, T, C) ∈ X × M × W × Q is a solution pair of (5). If ν, D c , γ, C and T satisfy the following uniqueness condition is unique solution pair of (5).

Several Iterative Methods Based on the Finite Element Discretization
In this section, we propose three iterative methods for the double-diffusive natural convection model. Then the stability and convergence of these iterative methods are considered. First, let 0 < h < 1 denote the mesh size which is a real positive parameter and K h = {K : K⊂ΩK =Ω} be a uniform partition ofΩ into non-overlapping triangles. Next, given a K h , we consider the finite element spaces where P i (K) represents the space of the order polynomial on the set K h , i = 1, 2. Please note that the Taylor-Hood element X h × M h satisfies the following discret inf-sup condition where the constant β > 0 is independent of h.
With the above notations, the finite element scheme for the natural convection problem is defined as follows: The following stability and convergence results of the numerical solutions to (12) are showed.
Moreover, the following error estimate holds where c is a positive constant depending on h.
In the following part of this section, we propose and analyse three iterative methods. for For the above three iterative methods, the initial guess Now, we will establish the stability and iterative error estimates of the presented iterative methods for the double-diffusive natural convection model. For the sake of simplicity, let (e n , η n , ξ n , for all n ≥ 0. Furthermore, the following iterative error bounds hold for all n ≥ 0.
Proof. First, the induction method is used to consider the stability of iterative method I.
Hence, we finish the induction method. Moreover, we consider the iterative error estimates of iterative method I. Making use of (12) and (13) yields the error equations a 0 (e n , v) + c 0 (u n−1 h , e n , v) + c 0 (e n−1 , u h , v) + D −1 a (e n , v) − d(η n , v) + d(q, e n ) = (β T ξ n g + β C δ n g, v), Setting ψ = ξ n , s = δ n in the second and the third equation of (20) and using (3), (17), and Theorem 3, we obtain Then, taking (v, q) = (e n , η n ) in the first equation of (20) and using (2), (4), (17), (21) and Theorem 3, we arrive at Hence, using uniqueness condition, we have Furthermore, subtracting (16) from (12), we obtain Applying (4), the Theorem 2 and the Theorem 3, we obtain which combines with (21) and (22), we arrive at for all n ≥ 0. Finally, applying the discrete inf-sup condition, from the first equation of (20) with q = 0, the error estimate of the pressure can be stated as follows.
for all n ≥ 0.

Theorem 5.
Under the assumptions of Theorem 3, suppose that the following condition (the strong uniqueness condition) holds. Then (u n h , p n h , T n h , C n h ) generated by iterative method II satisfies for all n ≥ 0. Furthermore, the following iterative error bounds hold for all n ≥ 0.
Proof. Combining with (19) and (23), it is found that (25) and (26) hold for n = 0. Supposing that (25) and (26) hold for n = k, we shall prove that they are valid for n = k + 1. Subtracting (14) from (12), we obtain the error equations Setting (v, q, ψ, s) = (e n h , η n h , ξ n h , δ n h ) in (27) with n = k + 1 and applying (2), (3), (4) and the assumptions on n = k, we have and ν ∇e k+1 Moreover, imply the strong uniqueness condition (24) on (29), we obtain ∇e k+1 Hence, making use of (30), we rewrite (28) as Combining the first equation of (27) with n = k + 1 and q = 0 and the discrete inf-sup condition, we have Furthermore, subtracting (16) from (14) with n = 1 that Then, taking ψ = T 1 h − T 0 h in the second equation of (33), we observe that Moreover, setting v = u 1 h − u 0 h in the first equation of (33), we obtain Combining (14) with n = 1 and using (34), we obtain In view of the strong uniqueness condition (24), we arrive at Next, taking (v, q, ψ, s) = (u n h , p n h , T n h , C n h ) in (14) with n ≥ 2, and using (2), (3) and (26), we obtain Similarly, we obtain Finally, it has ∇u n h 0 The proof is completed.

Theorem 6.
Under the assumptions of Theorem 3, suppose that the following condition (the stronger uniqueness condition), holds. Then (u n h , p n h , T n h , C n h ) defined by the iterative method III satisfies for all n ≥ 0. Furthermore, the following iterative error bounds hold for all n ≥ 0.

Numerical Experiments
In this section, several numerical experiments are presented to compare these iterative methods for the considered equations. We use the public finite element software FreeFem++ [28].

An Analytical Solution Problem
For numerical implementations, the iterative tolerance is 1.0 × 10 −5 . The first issue to be considered here is to compare these iterative methods for the stationary double-diffusive natural convection in the case of Ω = [0, 1] × [0, 1] ∈ R 2 , to reveal the relationship between the iterative methods and the viscosity. We consider the following exact solutions.
We compare the numbers of iteration and the computational time in Table 1. This table shows that all iterative methods run well in the case of ν = 1. When the viscosity number increases to ν = 1.0 × 10 −2 , iterative method III is divergent. Finally, iterative methods II and III can not export the data with ν = 1.0 × 10 −4 , iterative method I is still convergent. From a computational point of view, the calculation time of iterative method I and iterative method II is similar. However, iterative method II saves about 30% of calculation time compared iterative method III when ν = 1. Iterative method II saves about 35% of calculation time compared with iterative method I when ν = 1.0 × 10 −2 . We can conclude that iterative method III is the simplest method for a high viscosity number. The iterative method II is a fast and high accuracy method for a slightly lower viscosity number. Iterative method I is stable under uniqueness condition of weak solutions in the case of the lowest viscosity number. For three iterative methods, the relative error estimates are presented in Tables 2-4.  (7) -- Table 2. Comparison of three iterative methods using P 2 − P 1 − P 2 − P 2 (h = 1 64 and ν = 1).  Table 3. Comparison of three iterative methods using P 2 − P 1 − P 2 − P 2 (h = 1 64 and ν = 1.0 × 10 −2 ).  Table 4. Comparison of three iterative methods using P 2 − P 1 − P 2 − P 2 (h = 1 64 and ν = 1.0 × 10 −4 ).
Then, we show numerical velocity streamlines, isobars of pressure, isotherms, and isoconcentration lines obtained by three iterative methods with different viscosity numbers. We plot these results in . From these graphs, we obtain that the values of viscosity not only heavily impact on the velocity streamlines and the isobars, but also affect the isotherms and the isoconcentration lines. In fact, three iterations run well with ν = 1.0. However, iterative method III cannot run with ν = 1.0 × 10 −3 while iterative method II cannot export the data with ν = 1.0 × 10 −4 .
To consider the independency of mesh in a square cavity, we use iterative method I to calculate the model (1) under different mesh sizes. The results are presented in Figure 6. We can see that there is no difference in the calculation results under different mesh sizes, so we can verify the independence of the mesh size.

Conclusions
In conclusion, for solving stationary double-diffusive natural convection equations, three iterative methods have their own advantages under different viscosity numbers. In the case of 0 < σ < 1 4 , all methods can export data. Moreover, in the case of 1 4 ≤ σ < 1 3 , iterative method I and II can run well. Finally, in the case of 1 3 ≤ σ < 1, only iterative method I can export data.
From the perspective of physical applications, these finite element iterative methods can be used to simulate different double-diffusive natural convection models, such as the aluminum oxide nanofluid natural convection heat transfer, the natural convection flow of a suspension containing nano-encapsulated. Furthermore, some different boundary conditions of these models with some different calculation areas should be considered, such as the T-geometry enclosure porous cavity, L-geometry cavity, and porous cavity.