Abstract
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 of . 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.
1. 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 () with , the simulation time , and the force field ; find the velocity , the pressure , and the temperature satisfying [1] where is a unit vector in the direction of gravity, is the outward unit normal to , and where is a regular open subset of , is Prandtl number, is Rayleigh number, and 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–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. Çı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 conduction-convection 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–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–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 , the constant , and is the linear extrapolation of the velocity to from previous time levels. It is a three time levels scheme. Artificial viscosity stabilizations are introduced into the linear systems for and by adding and to the left-hand sides (LHS) and correcting them by and on the right-hand 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.
2. Preliminaries
2.1. The Variational Formulation of Natural Convection Problem
Let and denote the inner product and norm, respectively. Define the velocity space , the pressure space , the temperature space , and the divergence-free space as follows: where denotes the standard Sobolev space [32] with norm . All other norms will be clearly labeled with subscripts.
The weak formulation of problem (1) reads as follows: find for all , for all , such that Here, the skew-symmetric trilinear forms [1, 3, 5] which satisfy the following lemma.
Lemma 1 (see [25, 33]). Let (), for all ; . if, in addition, , and if, in addition, ,
We will use the Poincare inequality: for all , or , .
2.2. Finite Element Approximation
Assume to be a quasiuniform mesh of with mesh size ; let be a pair of conforming velocity-pressure-temperature finite element spaces which contain piecewise continuous polynomials of degree , , and , respectively, and satisfy the usual inf-sup condition and the following approximation properties [33]: The subspace of is given by
2.3. 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.
Definition 2 (see [5]). The modified Stokes projection operator for velocity and pressure is , , such that for all . In spaces (), Definition 2 reduces to the following: given , find such that for all . The modified Stokes projection operator for temperature is , , such that
Lemma 3 (error estimates of the modified Stokes projection). Suppose the inf-sup condition holds; then exists uniquely and satisfies where is a constant independent of and .
3. Numerical Scheme and Its Stability
3.1. Numerical Scheme
Let , , and . Define the linear extrapolation of the convecting velocity to by where is a known approximation to . The fully discrete CNSLE FE method of (1) is presented as follows.
Algorithm 4. Consider the following steps.
Step 1. Let be the modified Stokes projections of into spaces ; then at the first time level, find , for all , such that
Step 2. For , given , find , for all , such that
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 , and .
3.2. Stability of the Method
Theorem 5. Suppose , and . Algorithm 4 is unconditionally stable in the following sense, for any , , , and :
Proof. We first derive the stability of temperature and then of velocity . Choosing in (17) yields
Applying the Cauchy-Schwarz and Young’s inequalities gives
Thus, we get the bound of temperature on the first time level:
For , taking in (19) and using the Cauchy-Schwarz and Young inequalities lead to
Summing the above equation over from 1 to yields
Substituting the bound of (24) into the above relation, we obtain the result of (20).
Now we derive the stability of velocity . Taking in (16) and applying the Cauchy-Schwarz and Young inequalities, we have
Then
For , taking in (18) and making use of the Cauchy-Schwarz inequality and the Young’s inequality lead to
Multiplying the above relation by and summing over from 1 to give
Inserting the bound of (28) into the above relation and using the result (20), we get the desired result (21).
4. Error Estimate
Denoting and splitting the error terms , , and into the approximation errors , , and and the model errors , , and , respectively, where are the modified Stokes projections of into spaces (), respectively. For , 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 , , and (), respectively. If satisfies the condition then there is a constant such that, for any ,
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) at for gives
Adding and subtracting
to (35) and (36), respectively, we have
where
Using the error decomposition (31) and choosing , give
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
We next estimate the term of (40). For the linear terms of , using the Cauchy-Schwarz and Young’s inequalities together with Taylor expansion on , we get
where we have used the assumption (33) in the last equation. For the the nonlinear terms of in (40), by adding and subtracting some terms and using Lemma 1, we have
Using the inequalities of Lemma 1 and the regularity assumptions (32) on , then the second and the third terms of (45) are bounded by
Using Taylor expansion on , we write the first and the fourth terms of (45) as
with . Applying Lemma 1, Taylor expansion on , and the regularity assumptions (32), the terms of (47) are bounded by
Combining these inequalities with (40), we obtain
Step 1.3. Derive the error inequality of the energy error equation (41). Similarly, for the linear terms on the RHS of (41), we have
For the linear terms of , we arrive at
Using Lemma 1, we write the nonlinear terms of in (41) as
From Lemma 1, Taylor expansion on , and the regularity assumptions (32), it follows that
for any . Combining equations (52)-(53) with (41), we get
Step 1.4. Combining (49) with (54), multiplying it by , and using the approximation properties (9) yield
where we have used , . This ended the proof of Part 1.
Part 2. We derive the error estimate of Step 2 in Algorithm 4.
Step 2.1. Establis the error relations of Step 2. For and for all , taking the variational formulation (1) at gives
Then, subtract (18)-(19) from (56) to get
By adding and subtracting some terms to (57), we can obtain following error equations (recall that , ):
with
Using (31) and taking , give
Step 2.2. Derive the error estimate of the momentum equation (60). Basing on Definition 2, we obtain
Similarly, for the linear terms on the RHS of (60), we have
For the linear terms of in (60), we obtain
For the nonlinear terms of in (60), it follows that
where we have used Lemma 1 and rearranged some terms. Using definition (15) of the operator and the regularity assumptions (32) yields
We now estimate each nonlinear terms of (65). For the terms , by making use of Lemma 1, Taylor expansion on and Young’s inequality together with the regularity assumptions (32), we infer that
Adding and subtracting to , using Lemma 1, the inverse, and Young’s inequalities, we have (for detail, see [25])
Combining these inequalities with (60), we arrive at
Step 2.3. Establish the error estimate of the energy error equation (61). For the linear terms on the RHS of (61), we yield that
For the linear terms of , we obtain
By adding and subtracting to the nonlinear terms of in (61), we can obtain
Same as Step 2.2, basing on Lemma 1, Taylor expansion on , Young’s inequality, the inverse inequality, the regularity assumptions (32), and the estimate (66) of the operator , we can derive that
Substituting the above inequalities of Step 2.3 into (61), we obtain
Step 2.4. Combining (69) with (74) and rearranging some terms yield
Multiplying (75) by and summing it over from 1 to , using the approximation property (9), and substituting the error estimate (55) of the first time level into it, we can obtain
Applying condition (33) and the regularity assumptions (32) to (76) yields
Making use of the Gronwall lemma [22], for any , there exists a constant such that
We complete the proof of Theorem 6 by using the triangle inequality.
Corollary 8. Under the conditions of Theorem 6 and letting be continuous finite element spaces, then there exists a constant such that, for any ,
Remark 9. Corollary 8 shows that the optimal convergence rate of Algorithm 4 in space by using elements is for temperature and velocity in seminorm, with the choosing and . However, the error estimates in space are suboptimal for temperature and velocity in norm. We can find these facts from Tables 1–4 in the next section.
5. 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 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 elements are chosen for the FE spaces . All computations are carried out in the domain . 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 .
Example 1 (analytical solution). As in [38–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 . We choose the parameters that are , , and , 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 and . We list the errors and the convergence rates (CR) for velocity and temperature in (denoting as ) and (denoting as ) norm in Tables 1–4, respectively. Tables 1 and 2 are the results of Crank-Nicolson linearized extrapolation stabilized (CNSLE) FE method. Tables 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 . 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 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 norm is obtained, which indicate that our error estimate in 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 , , the parameters , , , the right hand functions , the time step , and the uniform mesh size . We performed the following study: calculating the problem from the rest until the approximation solution reaches a steady state. The criterion to stop this process is
where , denote , , 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–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–36]. Figure 2 presents the vertical velocity distribution at the midheight and the horizontal velocity distribution at the midwidth for , , , , 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 heat transfer coefficient in terms of the local Nusselt number is defined by
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 by
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 , 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–36].
(a)
(b)
(a)
(b)
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)
6. 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 norm are derived for velocity and temperature with a constrain condition. The derived theoretical results are supported by two numerical examples.
Conflict of Interests
The authors declare that there is no conflict of interests regarding the publication of this paper.
Acknowledgments
The authors thank the editor and reviewers for their criticism, valuable comments, and suggestions which helped to improve the results of the paper. This work is in part supported by the NSFC (no. 11171269), the Ph.D. Programs Foundation of Ministry of Education of China (no. 20110201110027), the Fundamental Research Funds for the Central Universities, the China Postdoctoral Science Foundation (no. 2013M531311), the Educational Commission of Henan Province of China (nos. 14B110020, 14B110021, and 14B110025), the Henan Scientific and Technological Research Project (no. 132102310309), the Doctoral Foundation of Henan University of Science and Technology (no. 09001625), the Science Foundation for Cultivating Innovation Ability of Henan University of Science and Technology (no. 2014ZCX009) and the Youth Scientific Foundation of Henan University of Science and Technology (no. 2012QN029).