The Crank-Nicolson Extrapolation Stabilized Finite Element Method for Natural Convection Problem

This paper studies a fully discrete Crank-Nicolson linear extrapolation stabilized finite element method for the natural convection problem, which is unconditionally stable and has second order temporal accuracy ofO(Δt +hΔt+h). A simple artificial viscosity stabilized of the linear system for the approximation of the new time level connected to antidiffusion of its effects at the old time level is used. An unconditionally stability and an a priori error estimate are derived for the fully discrete scheme. A series of numerical results are presented that validate our theoretical findings.


Introduction
Natural convection flow has many thermal engineering applications such as in double-glazed windows, solar collectors, cooling devices for electronic instruments, gas-filled cavities around nuclear reactor cores, and building insulation.Typically, fluid flow and heat transfer are governed by the partial differential equation system of momentum, mass, and energy conservation, but in the case of natural convection, the so-called Boussinesq approximation is generally used.The natural convection problem which we consider is for bounded, polyhedral domains Ω  ⊂ Ω in R  ( = 2, 3) with dist(Ω  , Ω) > 0, the simulation time  * , and the force field  : Ω × (0,  * ] → R; find the velocity  : Ω × (0,  * ] → R  , the pressure  : Ω × (0,  * ] → R, and the temperature  : Ω × (0,  * ] → R satisfying [ where  is a unit vector in the direction of gravity,  is the outward unit normal to Ω, and Γ  = Ω \ Γ  where Γ  is a regular open subset of Ω, Pr is Prandtl number, Ra is Rayleigh number, and  > 0 is thermal conductivity parameter.Moreover,  =   in Ω  and  =   in Ω  , where   and   are positive constants.A global-in-time existence result for a more general natural convection problem (Navier-Stokes/Fourier model) is given in [2].Many authors have worked hard to study for a great variety of efficient numerical schemes for the natural convection problem [3][4][5][6][7][8][9][10][11][12][13][14][15][16][17] and relevant research [18,19].We mention only a few papers here.[3,4] are the early papers by using mixed finite element (FE) method.C ¸ıbık and Kaya [5] have formulated a projection-based stabilization FE technique for solving the steady-state natural convection problems.The global stabilizations are added for both velocity and temperature variables and these effects are subtracted from the large scales.Galvin et al. [7] consider the problem of poor mass conservation in mixed FE algorithms for flow problems with large rotation-free forcing in the momentum equation.Zhang et al. [8] have presented a subgrid stabilized defect-correction method for steady-state natural convection problem.Shi and Ren [11] have proposed a least squares Galerkin-Petrov nonconforming mixed FE method for stationary conduction-convection problems.Luo et al. [12] have given an optimizing reduced Petrov-Galerkin least squares mixed FE for the nonstationary conductionconvection problem.Boland and Layton [1] have derived stability properties and error estimates for the mixed FE spatial discretization case when used to approximate heat flow in a fluid enclosed by a solid medium.Benítez and Bermúdez have presented a second order Lagrange-Galerkin method for natural convection problems in [17].In [20,21], a stability analysis of thermal natural convection in superposed fluid and porous layers is carried out.
Our goal in this paper is to solve time-dependent natural convection problem efficiently and accurately.Usually fully implicit schemes are (almost) unconditionally stable, but one has to solve a system of nonlinear equations at each time step.Although an explicit scheme is much easier in computation, it suffers a restricted time step size from the stability requirement.A popular approach is based on an implicit scheme for the linear term and a semi-implicit scheme or an explicit scheme for the nonlinear term.There are numerous works on the Crank-Nicolson and relevant high order scheme for the Navier-Stokes (NS) equations [22][23][24][25][26][27][28][29][30].The Crank-Nicolson linear extrapolation (CNLE) scheme for NS equations was first studied by Baker in [23].The second and third order CNLE methods are introduced and analysed in [24].A stabilized extrapolated trapezoidal FE method is given in [25] for the NS equations.A variational multiscale method based on the CNLE scheme for the NS equations is proposed in [26].He et al. [27][28][29] have studied the NS equations based on the Crank-Nicolson extrapolation (Crank-Nicolson/Adams-Bashforth, or two level methods) schemes.In [31], we have studied fully implicit Crank-Nicolson scheme for natural convection problem.
We consider herein a simple, second order accurate, and unconditionally stable fully discrete Crank-Nicolson linear extrapolation stabilized (CNSLE) FE method for natural convection problem which requires the solution of one linear system per time step.Suppressing the spatial discretization, the method is where the time step Δ > 0, the constant  =  (1), and  +1/2 = (3/2)  − (1/2) −1 is the linear extrapolation of the velocity to  +1/2 from previous time levels.It is a three time levels scheme.Artificial viscosity stabilizations are introduced into the linear systems for  +1 and  +1 by adding −ℎΔ +1 and −ℎΔ +1 to the left-hand sides (LHS) and correcting them by −ℎΔ  and −ℎΔ  on the righthand sides (RHS), respectively.To the best of the authors' knowledge, there are no papers dealing with the error analysis of the fully discrete CNLE FE method for natural convection problem.
The paper is organized as follows.Section 2 collects some preliminaries for the analysis that follows.In Section 3, we give the fully discrete CNSLE FE method and prove it is unconditionally stable in Theorem 5.In Section 4, error estimates for velocity and temperature are derived in Theorem 6. Numerical experiments are described in Section 5. Conclusions follow.Problem.Let (⋅, ⋅) and ‖ ⋅ ‖ denote the  2 (Ω) inner product and norm, respectively.Define the velocity space , the pressure space , the temperature space , and the divergence-free space  as follows:

Finite Element Approximation.
Assume Ω ℎ = {} to be a quasiuniform mesh of Ω with mesh size 0 < ℎ < 1; let ( ℎ ,  ℎ ,  ℎ ) ⊂ (, , ) be a pair of conforming velocitypressure-temperature finite element spaces which contain piecewise continuous polynomials of degree ,  − 1, and , respectively, and satisfy the usual inf-sup condition and the following approximation properties [33]: The subspace  ℎ of  ℎ is given by

The Modified Stokes Projection.
The following modified Stokes projection operators are similar than the one in [5].
For the reader's convenience, we only present the definition and the error results of the projection operators.We can easily derive the result.

Numerical Scheme and Its Stability
where   ℎ is a known approximation to (,   ).The fully discrete CNSLE FE method of (1) is presented as follows.
Algorithm 4. Consider the following steps.
We find that CNLES FE method requires the solution of only one linear problem at each time step; thus it needs less time than fully implicit Crank-Nicolson scheme.Denote  min = min(  ,   ), and  max = max(  ,   ).

Stability of the Method
Theorem 5. Suppose 0 <  min ≤  max < ∞, and  ∈  2 (0, ;  −1 (Ω)).Algorithm 4 is unconditionally stable in the following sense, for any ℎ, Δ > 0,  ⩾ 0, and 0 ⩽  ⩽  − Proof.We first derive the stability of temperature  and then of velocity .Choosing  ℎ = ( 1 ℎ +  0 ℎ )/2 ∈  ℎ in (17) yields Applying the Cauchy-Schwarz and Young's inequalities gives where (ũ +1 , p+1 , T+1 ) are the modified Stokes projections of (( +1 ), ( +1 ), ( +1 )) into spaces ( ℎ ,  ℎ ,  ℎ ), respectively.For  ⩾ 2, we assume the exact solution of problem (1) satisfies following regularity assumptions: Theorem 6. Assume that the solution (, , ) of problem (1) satisfies regularity assumptions (32).Let the finite element spaces ( ℎ ,  ℎ ,  ℎ ) include continuous piecewise polynomials of degree ,  − 1, and  ( ⩾ 2), respectively.If Δ satisfies the condition then there is a constant Ĉ = Ĉ(, , , Ω, , , , ) such that, for any 0 Remark 7. We see that our proof is conditionally convergent.We think that the condition ( 33) is not necessary theoretical, and it might be removed.But we cannot make it herein and will study it in the future.Now we give the outline of the proof.The proof is proved in Part 1 (establishing the error estimate for Step 1 in Algorithm 4) and Part 2 (deriving the error estimate for Step 2 in Algorithm 4).In each parts, we first derive the error equations of the momentum equation and the energy equation; then give the error estimates of the momentum error equation and the energy error equation, respectively, and finally derive Theorem 6 by using the Gronwall lemma and the triangle inequality.
Proof.Part 1. Derive the error estimate of Step 1 in Algorithm 4.
Step 1.1.Establishe the error equations of the momentum equation and the energy equation.Subtracting equations ( 16) and ( 17) from ( 4) Adding and subtracting to ( 35) and (36), respectively, we have where Using the error decomposition (31) and choosing Step 1.2.Derive the error inequalities of the momentum equation (40).From the definition of the modified Stokes projection ( 11)-( 13), we have For the linear terms on the RHS of (40), applying the Cauchy-Schwarz, Young inequalities.and the approximation properties (9) gives where we have used the assumption (33) in the last equation.

Numerical Experiments
We first verify the convergence rates and the effectiveness of our methods by an analytical solution.Then we test the stability of the method in case Ra is large in a squared cavity with the left wall heating.The code was implemented using the software package FreeFEM++ [37].The pairs of continuous ( 2 ,  1 ,  2 ) elements are chosen for the FE spaces ( ℎ ,  ℎ ,  ℎ ).All computations are carried out in the domain Ω = [0, 1] 2 .The uniform mesh is obtained by dividing Ω into squares and then drawing a diagonal in each square in the same direction.We set the parameter  = (1).
Example 1 (analytical solution).As in [38][39][40][41], to obtain an analytical solution for the considered problem, a right-hand side function is added to the momentum equation of (1).The analytical solution is The initial velocity field and temperature field are equal to the analytical solution at time  = 0. We choose the parameters that are  = 1.0,Pr = 1.0, and Ra = 100, respectively.The mesh and time step sizes scalings are set that for a refinement, each of ℎ and Δ get cut in half, where the final time is chosen as  * = 0.1 and Δ = (1/10) ℎ.We list the errors and the convergence rates (CR) for velocity  and temperature  in  ∞ (0,  * ;  2 (Ω)) (denoting as  ∞,0 ) and  2 (0,  * ;  1 (Ω)) (denoting as  0,1 ) norm in Tables 1-4, respectively.Tables 1  3 and 4 are the results of Crank-Nicolson stabilized (CNS) FE method, which the nonlinear system is solved by Newton iterative at each time step and the iterative tolerance is 10 −8 .We find that the errors listed in Tables 1 and 3 and in Tables 2 and 4 are similar, respectively, which means that the CNSLE method is comparable to the CNS method.From Tables 1-4, we find that a quadratic convergence rate for the computed velocity and temperature in  1 seminorm is obtained, which test and verify the theoretical results.However in Tables 1-4, it is easy to see that a cubic convergence rate both for the computed velocity and temperature in  2 norm is obtained, which indicate that our error estimate in  2 norm is suboptimal.As expected, since CNSLE is the linearized version of the CNS method, which does not include any iterations in computing, CPU cost by CNSLE method is relatively less than by CNS method.The numerical results agree well with the theoretical predictions.
Example 2 (natural convection in a squared cavity with the left wall heating).Now we present natural convection in a squared cavity with the left wall heating.The boundary conditions are given in Figure 1.We choose the initial conditions  0 = 1,  0 = 1, the parameters  = 1.0,Pr = 0.71, 10 4 ⩽ Ra ⩽ 10 7 , the right hand functions  = 0, the time step Δ = 0.1, and the uniform mesh size ℎ = 1/64.We performed the following study: calculating the problem from the rest   where  + 1,  denote  +1 ,   , respectively.We describe the final results of the problem at its steady state in Tables 5-7 and Figures 2-4, which are according to the results of [16,17,[34][35][36].The numbers in parenthesis of Tables 5-7 correspond to the used mesh.Tables 5 and 6 present the maximum vertical velocity at midheight and horizontal velocity at midwidth and compare the quantitative results with the available benchmark solutions [16,[34][35][36].Figure 2 presents the vertical velocity distribution at the midheight and the horizontal velocity distribution at the midwidth for Ra = 10 4 , 10 5 , 10 6 , 10 7 , respectively.We find that the differences in the profiles are getting larger along with the increase of Rayleigh numbers.It is seen that the agreement is excellent even at higher Rayleigh numbers.
The average Nusselt numbers at the hot wall of the cavity are given in Table 7, which are according to the results of [16,34,35].Figure 4 presents streamlines and isotherms for Ra = 10 4 , 10 5 , 10 6 , 10 7 , respectively.For the streamline patterns, it is easy to see that elliptical vortex at the cavity center break up into two vortices tending to approach to the corners differentially heated sides of the cavity with Rayleigh number increasing.With the increase of Rayleigh numbers Ra, the temperature convection becomes increasingly prominent, the isotherms gradually transform into horizontal except for the immediate neighborhood of the hot and cold walls which remain parallel to the isothermal vertical walls.These results are keeping with the results of [16,17,[34][35][36].

Conclusions
In this paper, a fully discrete stabilized finite element method based on the Crank-Nicolson linear extrapolation scheme in time for natural convection problem is given.The method is unconditionally stable.Optimal error estimates in H1 seminorm and suboptimal error estimates in  2 norm are derived for velocity and temperature with a constrain condition.The derived theoretical results are supported by two numerical examples.

Figure 1 :
Figure 1: Natural convection in cavity: the physical domain with its boundary conditions.

∫ 1 0
The heat transfer coefficient in terms of the local Nusselt number is defined byNu local (, ) = −   .(82)Variation of local Nusselt number at hot wall and cold wall of cavity for different Rayleigh numbers is plotted in Figure 3.We calculate the average Nusselt number on the vertical boundary of the cavity at  = 0 by Nu = Nu local (, ) .

Figure 2 :
Figure 2: Natural convection in cavity: comparison of vertical velocity at the midheight (a) and horizontal velocity at midwidth (b) for different Rayleigh numbers.

Figure 3 :
Figure 3: Natural convection in cavity: comparison of local Nusselt numbers along the hot wall ( = 0) (a) and the cold wall ( = 1) (b) for different Rayleigh numbers.

Table 7 :
Comparison of average Nusselt number on the vertical boundary of the cavity at  = 0 with mesh size used in computation for Example 2.

Table 5 :
Comparison of maximum vertical velocity at  = 0.5 with mesh size used in computation for Example 2.

Table 6 :
Comparison of maximum horizontal velocity at  = 0.5 with mesh size used in computation for Example 2.