Recovery-Based Error Estimator for Natural Convection Equations Based on Defect-Correction Methods.

In this paper, we propose an adaptive defect-correction method for natural convection (NC) equations. A defect-correction method (DCM) is proposed for solving NC equations to overcome the convection dominance problem caused by a high Rayleigh number. To solve the large amount of computation and the discontinuity of the gradient of the numerical solution, we combine a new recovery-type posteriori estimator in view of the gradient recovery and superconvergent theory. The presented reliability and efficiency analysis shows that the true error can be effectively bounded by the recovery-based error estimator. Finally, the stability, accuracy and efficiency of the proposed method are confirmed by several numerical investigations.

1. Introduction

Natural convection (NC) equations for buoyancy-driven fluid often appear in practical problems. The stationary NC equation is a coupling equation for the incompressible flow and heat transfer process of viscous fluid, in which the incompressible fluid can be characterized by Boussinesq’s approximation. In atmospheric dynamics, it is an important forced dissipative nonlinear system. It contains the velocity field, pressure and temperature, which we can analyze from the thermodynamic point of view. The motion viscosity of the fluid produces a certain amount of heat, so the motion of the fluid must be accompanied by the conversion of temperature, velocity and pressure. Therefore, the study of this nonlinear system is of great significance. Christie and Mitchell [1], Boland and Layton [2,3] and others have extensively researched the numerical analysis and numerical results of NC equations (see [4,5,6,7,8,9]). In recent years, there has been continuous research on NC equations, such as in [8,10,11,12,13], which studied the natural convection of cavity-filled nanofluids. Ahmad [14] studied the effect of viscosity and thermal conductivity on natural convection in exothermic catalytic chemical reactions on curved surfaces, and in [10,15,16,17,18,19], Singh et al. considered the effect of factors such as Lorentz force on natural convection.

Combined with previous research results, the dimensionless parameter (Rayleigh number) plays an important role in NC equations. It is well known that the buoyancy force is the driving force of NC equations, where the buoyancy force is the density difference caused by the temperature difference. The Rayleigh number represents the ratio of the buoyancy force to the viscous force, which characterizes the relative strength of the buoyancy-driven inertial force to that of the viscous force [20,21]. When , the buoyancy force is much larger than the viscous force, and the convection caused by the buoyancy force is significant, which leads to the convection dominance problem [22]. Therefore, there are a number of effective numerical techniques dedicated to solving this difficulty: the variational multiscale method in [23], the defect-correction method in [24,25,26,27], etc. Among them, defect correction is one of the commonly used methods to address large number problems.

The defect-correction method (DCM) is an iterative improvement technique for improving the accuracy of computational solutions without introducing mesh refinement [28,29,30,31]. In general, the basic idea can be briefly described as follows. Assume that M is a mapping , where X and Y are typical linear spaces, and assume that our goal is to find a good approximation of , where . Assume that we actually solve a related problem , whose solution we denote by . Then, we find the residual or defect and use it to solve the error or correct , either by nonlinearly updating or by linearly updating to obtain the corrected solution , setting ; obviously, this process can be repeated. If the mapping M is nonlinear, then the linearization of is usually chosen in the correction phase to simplify the computation. In the case of solving linear systems, the well-known iterative refinement method is an example of a defect-correction technique.

For computational practicality, the minimum-order coordinated velocity, pressure and temperature finite elements with are used to discretize NC equations. The finite element approximation of the pseudo-stress tensor is the discontinuous slice constant of We use the superconvergent slice to recover the amount of continuous space The approximate solution computed by the local recovered error estimator is closer to the true solution than the approximate solution of the general finite element method. The a posteriori error has been extensively studied (see [32,33,34,35,3637)] Therefore, we use the adaptive defect-correction method to solve NC equations.

This paper is divided into four parts. First, it introduces the main properties of NC equations and the classical results. The second part proposes the defect-correction method and presents the analysis of its error estimation. The third part is combined with the restorative error estimator to solve NC equations. Finally, the validity and reliability of our method are verified by numerical experiments.

2. Preliminaries

Suppose that is a bounded, convex open area whose boundary is Lipschitz continuous, so , where is a regular open subset of . We consider the following NC equations [7,9]: where is fluid velocity, is fluid pressure, is temperature, is the external force function, and are the Prandtl number and Rayleigh number, respectively, is the thermal conductivity parameter, and is the two-dimensional unit vector. Next, we introduce some Hilbert spaces:

We denote the inner product and norm of by and . We define the inner product in space X and W. and its norm are written as and . In addition, the dual space for is . The dual space of X is [2].

For simplicity, we define u and v such that they belong to the same finite element space X. We define two continuous bilinear forms and and a trilinear form :

In addition, we need to define two bilinear forms and a trilinear form in and , respectively.

According to the bilinear form defined above and the trilinear form, the following two known conclusions are obtained [2,3]: where are two fixed constants that depend only on .

The variational form of the NC Equation (3) is as follows: solving , :

Then, there exists a unique result for the following solution [2,3].

There exists at least a solution pair

where

then the solution pair

3. Defect-Correction Method

3.1. The Finite Element Discrete Form of NC Equations

The DCM and the corresponding error estimates for steady-state NC equations are given in this section. For convenience, to represent different constants for different conditions, the constant C is independent of the grid parameter h, but it is always far from the infinity bound.

We apply an edge-to-edge triangulation of into , whose minimum angle is bounded away from zero, and to represent the union of all elements K and all edges e of the elements in , respectively. In addition, and denote, as usual, the diameters of an edge e and element K, respectively. Then, velocity–pressure–temperature finite element spaces can be constructed on . First, we define the finite element space : where represents piecewise linear polynomials on K, and represents piecewise constant polynomials on K.

According to the above definition, the discrete form of the original equation is given: the solution of satisfies the requirement where .

For the finite element approximation (3), we present some results in [2,3] as follows.

Let

For convenience, we define where , .

Since we choose an unstable finite element pair ––, we consider the following penalty jump stabilization method: where is the set of inner edges of , is the jump of q on edge e, and is a given stable parameter.

Therefore, formula (4) can be written as

We assume that there exists an interpolation operator of the Clément type in the finite element space. Specifically, satisfies the following elementwise error estimate (see [24,27]). For , there is , , such that

3.2. Application to Natural Convection Equations

We define and analogously: where , and .

In this way, formula (5) can be written as where , and .

It is well known that is a dimensionless result of the ratio of buoyancy to viscous force in NC equations. When buoyancy is far greater than viscous force, the convection phenomenon driven by buoyancy is more significant than the diffusion phenomenon caused by viscous force; then, convection dominates, and the regularity of is poor. To overcome this problem, let , be a more stabilized or regularized approximation of . The DCM is given as follows:

First, solve to satisfy

Next, the above form is corrected by correction steps, which makes the obtained more closely approximate for :

With , , and , the defect-correction discretization is applied to the incompressible NC equations (see [6]), where the trilinear term and are discretized by the Oseen scheme. In the calculation process, in order to produce better results, is usually disturbed by the local average or flux limiters or by merging the appropriate subgrid model.

The defect-correction method for problem (3) can be described as follows:

Step 1. Find a stabilized finite element iterative solutionsuch thatfor all, .

Step 2. For j = 0,1,2,. . ., find a corrected solutionsuch thatwith.

3.3. Error Analysis

Similarly, the term is a consistency error, which is a higher-order term, so we only need to estimate and .

Suppose that

where

First, we consider . We integrate by parts over each element and denote the collection of edges of in the interior of by . Using the Cauchy–Schwarz inequality on each element K and edge e:

We assume and , where , .

According to the usual Galerkin formulation, is identically zero. For this last term involving , let satisfy . Then,

This completes the proof for Theorem 1. □

Note that this amounts to replacing by , where . The residual term is , and is a computable constant. In addition, let be the set of all interior sides in , and for any piecewise constant q, let denote their jumps on the side , where and are two triangulations sharing the common side e. Let be the length of edge , and define the edge norm as

4. Recovery Error Estimator for NC Equations

4.1. Recovery Error Estimator

In this section, we construct a recovery error estimator and analyze its properties. Assuming that is the numerical result, we consider the pseudo-stress tensor and its finite element approximation .

I is -identity matrix. The main idea of the Zienkienwicz–Zhu estimator is to restore the discontinuous finite element gradient to a continuous recovery term (see [37]).

In , N, , and are respectively defined as the set of all vertices in , the set of all vertices within in and the base function of i () that shares a collection of all units of a vertex ().

Combined with the superconvergent piece recovery technique in the literature [37], the piecewise constant tensor is restored to continuous, and a recovery pseudo-stress tensor is obtained.

For any vertex and its patch , the following is defined: where is the area of triangulation K. Details can be found in [37].

A recovery error estimator is constructed. The local estimator and global estimator are defined as

The recovery error estimator has the following important properties. The proof is similar to that in the literature ([37]) and is not presented in detail here.

Before moving to the next subsection, we illustrate the connection between and , where is the unit outward normal vector on edge e, which plays an important role in the a posterior error of this recovery-based error estimator. In two dimensions, denotes the unit tangential vector along e. denote the jumps of the normal and tangential component of on side e, respectively. Since , we have

Therefore,

From the above relation, we have

Based on the definitions of

According to the definitions of

According to the Brezis–Gallonet inequality in [37] and the method in [37], the conclusion can be easily obtained. □

4.2. The Reliability Analysis

A basic framework of posteriori error estimation for nonlinear stability problems is proposed, and the basic equivalence of the error to residual error is established. We introduce the main results and apply them to estimate the restorative error estimator.

Suppose that

According to (12), with the Cauchy–Schwarz inequality, we have

According to the definition of , the approximate solution satisfies

At the same time .

With the above lemma, we obtain

Theorem 4 is proved. □

4.3. Effectiveness Analysis

From Lemma 1, the following can be obtained:

In order to prove the validity of the recovery error estimator , we first need to estimate and

Suppose that

From the definition of , the following can be obtained:

Specific details can be found in [37]. □

Next, we define the operator : .

Suppose that

Combined with the above estimates, Lemma 1 and Lemma 3 are applied.

The proof of validity is complete. □

5. Numerical Experiment

In this section, four numerical examples are given to verify the effectiveness of the proposed method combined with the adaptive method.

For the sake of convenience, a few definitions are given below: where the validity index denotes the ratio of the estimator error to the true error , and if tends to 1, it means that our estimator error is asymptotically equivalent to the true error, thus verifying our validity and reliability.

the number of degrees of freedom of triangulation ;

denotes the relative value of norm;

denotes the relative value of global recovery-based estimator;

denotes the convergence rate of the error;

denotes the convergence rate of the error;

denotes effectivity index for the global recovery-based estimator ,

5.1. Smooth True Solution

The purpose of the first example is to solve a smooth true solution problem in the region and verify our method, which is effective for the smooth true solution. We define the true solution , pressure p and temperature T as follows:

Figure 1 shows the initial mesh of the adaptive method and the adaptive grid. From Table 1 and Table 2, if only the defect-correction method is used to solve the NC equations, more degrees of freedom are needed to achieve better accuracy, while the adaptive method can yield better accuracy with less vertex information. For example, the accuracy of 5000 degrees of freedom in Table 2 is similar to that of 2773 degrees of freedom in Table 1. This shows that the adaptive method is more economical. Table 3 presents the comparison of the results of our method (DCM) and the adaptive algorithm without the defect-correction method (NDCM). We selected the values of the error estimator for similar degrees of freedom with different . From the table, we can see that the error results without the defect-correction method when are very large, and when is greater than , the results are not counted. Our method, on the other hand, is still relatively stable at . It can be seen that the posteriori error estimator based on the defect-correction method is applicable. tends to 1, which indicates that our estimators and real errors are effective and progressive, indicating that our method is effective.

Figure 1

The initial mesh and adaptive final mesh.

Table 1

Adaptive results of smooth solution, where .

Level	DOF	e_r	e_rrate	η_r	η_rrate	I_eff
0	370	0.1317	-	0.1114	-	1.1824
1	691	0.0920	1.1482	0.0813	1.0056	1.1309
2	1395	0.0614	1.1523	0.0576	0.9820	1.0652
3	2773	0.0427	1.0528	0.0408	1.0018	1.0467
4	5494	0.0305	0.9848	0.0291	0.9898	1.0485

Table 2

Uniform mesh results of smooth solution, where .

Level	DOF	e_r	e_rrate	η_r	η_rrate	I_eff
0	450	0.1211	-	0.1133	-	0.9356
1	800	0.0928	0.9252	0.0886	0.8548	0.9547
2	1250	0.0758	0.9068	0.0732	0.8557	0.9644
3	5000	0.0415	0.8691	0.0409	0.8397	0.9855
4	6050	0.0383	0.8419	0.0379	0.7993	0.9819

Table 3

Comparison of DCM and non-DCM, with ; .

	DOF	Ra = 10^3	DOF	Ra = 10^4	DOF	Ra = 10^5
DCM	2773	0.0408	2670	0.0653	2272	0.0923
NDCM	2392	0.1114	2677	13.91	NAN	NAN

5.2. L-Shape Domain Problem

The second example is a flow problem in the L-shape domain , with , and 1000, . The NC equations (1) is determined by exact velocities and , pressure p and temperature T:

We note that both velocity u and pressure p are continuous in the domain. However, it is clear that u and p are singular at the point and along the line , respectively.

To show the efficiency of the error estimator, we compare numerical results for adaptive refinements with those for uniform refinement. Figure 2 presents the initial mesh (left), the final uniform mesh (middle) and the final adaptive mesh (right). Table 4 and Table 5 show the numerical results for uniform/adaptive refinements with the recovery-based estimator .

Figure 2

The initial mesh, uniform final mesh and adaptive final mesh.

Table 4

Uniform mesh and adaptive refinement results of DCM, where .

Level	Uniform Mesh			Adapt Mesh
	DOF	η_r	η_rrate	DOF	η_r	η_rrate
0	676	0.3523	-	830	0.2547	-
1	1516	0.2627	0.7267	1337	0.1631	1.8698
2	2696	0.2043	0.8735	2344	0.1156	1.5798
3	4282	0.1665	0.8844	4558	0.0833	0.9855
4	6170	0.1419	0.8753	6315	0.0704	1.0135

Table 5

Uniform mesh and adaptive refine results of DCM, where .

Level	Uniform Mesh			Adapt Mesh
	DOF	η_r	η_rrate	DOF	η_r	η_rrate
0	676	0.7098	-	676	0.7098	-
1	1516	0.4955	0.8900	1423	0.4420	1.2728
2	2696	0.3805	0.9174	1952	0.3770	1.0065
3	4282	0.3059	0.9434	3622	0.2778	0.9879
4	6170	0.2602	0.9500	4970	0.2358	1.0362

According to in Figure 2c, we can see that near the singular point of u and the line of with the singularity of p, mesh refinement is carried out effectively, and the desired results are obtained, which shows that the adaptive mesh refinement is effective.

As can be seen from Table 4, Table 5 and Table 6, our adaptive method can obtain better accuracy with less vertex information, which is better than the uniform grid. For example, when , to achieve accuracies of 0.1631 and 0.1665, the adaptive method only needs 1337 degrees of freedom, while the uniform grid requires 4282 degrees of freedom. When , the error is 0.2778 and 0.2602; the adaptive method only needs 3622 degrees of freedom, while the uniform grid needs 6172 degrees of freedom. This shows that our method is more efficient. Furthermore, as becomes larger, our method calculates results that are superior to FEM.

Table 6

FEM with ; .

DOF	Ra=100		Ra=1000
	η_r	η_rrate	η_r	η_rrate
676	0.3610	-	0.4953	-
1516	0.2688	0.8900	0.3637	1.2728
2696	0.2086	0.9174	0.3061	1.0065
4282	0.1695	0.9434	0.2814	0.9879

5.3. Thermally Driven Flow

In the second numerical example, we consider the problem of square cavity flow without a true solution. As shown in Figure 3, the definition of the calculation area is , with and , and the boundary definition is as follows: the left boundary and the bottom are , and the bottom edge is . The rest is , and the speed is the 0 Dirichlet condition.

Figure 3

The model of thermally Driven Flow.

Figure 4 shows the initial mesh and the adaptive encryption result grid. Figure 5 and Figure 6 are the velocity streamline and isothermal chart for different Rayleigh numbers: and . Figure 7 and Figure 8 present the results of the streamline of the velocity and the isotherm diagram calculated without the defect-correction method at different numbers, and it can be seen that as increases, the calculated results become more unstable [38]. In Figure 9, we also show the vertical velocity distribution at and the horizontal velocity distribution at , which are very popular graphical illustrations in experimental studies of thermally driven cavities. We see that the profiles increasingly vary as the Rayleigh number increases, which is consistent with previous studies [5,7,39].

Figure 4

The initial mesh and adaptive final grid.

Figure 5

The streamline of the velocity with different .

Figure 6

The isotherm diagram with different .

Figure 7

The streamline of the velocity without the defect-correction method with different .

Figure 8

The isotherm diagram without the defect-correction method with different .

Figure 9

The values of vertical and horizontal velocity. (a) The vertical velocity at mid-height; (b) The horizontal velocity at mid-width.

Table 7 and Table 8 present the peak values of vertical velocity at and horizontal velocity at for different Rayleigh numbers, respectively, where the number ’m’ in parentheses corresponds to the degree of freedom of the used mesh. It is easy to see that our results are concordant with benchmark data [5,7,39].

Table 7

Comparison of maximum vertical velocity at with mesh size used in the computation.

Ra	Our Method	Ref. [5]	Ref. [39]	Ref. [7]
10^4	19.63(2445)	19.51(1681)	19.79(10,201)	19.74(1024)
10^5	70.22(2537)	68.22(6561)	70.63(10,201)	69.34(1024)

Table 8

Comparison of maximum horizontal velocity at with mesh size used in computation.

Ra	Our Method	Ref. [5]	Ref. [39]	Ref. [7]
10^4	16.37(2445)	16.18(1681)	16.10(10,201)	16.31(1024)
10^5	34.60(2537)	34.81(6561)	34(10,201)	35.83(1024)

5.4. Bernard Convection Problem

In this experiment, we consider the domain with the forcing and . Figure 10 displays the initial and boundary conditions for velocity u and temperature T: on , on the lateral boundaries, with a fixed temperature difference between the top and bottom boundaries.

Figure 10

Physics model of Bernard convection problem.

The results of numerical streamlines and numerical temperature are shown in Figure 11 and Figure 12 with , and . The results coincide with the phenomenon described in the literature [38], indicating that our adaptive DCM is suitable.

Figure 11

Initial mesh and final adaptive meshes.

Figure 12

Calculation results of Bernard convection problem with .

6. Conclusions

In this work, a recovery-type error estimator for NC equations for the DCM finite element method is constructed based on the superconvergent patch recovery technique. Because the construction process is simple and does not involve the numerical solution of discrete NC equations, it is applicable to any stabilization method that can stabilize the lowest-order finite element to ––. The results of numerical experiments fully verify the correctness of the theoretical analysis. The efficiency and stabilization of the new error estimator for the considered problems are proved. Furthermore, our method can be extended to more general fluid dynamics equations.