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.


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 Ra (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 Ra 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 Ra 1, 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 Ra 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 M : X → Y, where X and Y are typical linear spaces, and assume that our goal is to find a good approximation of x 1 , where M(x 1 ) = 0. Assume that we actually solve a related problem M 1 (x) = 0, whose solution we denote by x 0 . Then, we find the residual or defect d 0 = M(x 0 ) and use it to solve the error or correct e 0 , either by nonlinearly updating M 1 (x 2 ) − M 1 (x 0 ) = d 0 or by linearly updating M 1 (x 0 )e 0 = d 0 to obtain the corrected solution x 2 , setting x 2 = x 0 + e 0 ; obviously, this process can be repeated. If the mapping M is nonlinear, then the linearization of M 1 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 P 1 -P 0 -P 1 are used to discretize NC equations. The finite element approximation of the pseudo-stress tensor sigma = Pr∇u − pI + κdiag(∇T) is the discontinuous slice constant of sigma h = Pr∇u h − p h I + κdiag(∇T h ). We use the superconvergent slice to recover the amount of continuous space G(σ h ). The approximate solution computed by the local recovered error estimator ||σ h − G(σ h )|| K 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][36][37]). 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.

Preliminaries
Suppose that Ω = [0, 1] × [0, 1] is a bounded, convex open area whose boundary ∂Ω is Lipschitz continuous, so Γ T = ∂Ω \ Γ B , where Γ B is a regular open subset of ∂Ω. We consider the following NC equations [7,9]: in Ω, u = 0, in ∂Ω, γ is the external force function, Pr and Ra > 0 are the Prandtl number and Rayleigh number, respectively, κ > 0 is the thermal conductivity parameter, and j = (0, 1) T is the two-dimensional unit vector. Next, we introduce some Hilbert spaces: We denote the inner product and norm of L 2 (Ω) by (·, ·) and || · || 0 . We define the inner product in space X and W. H m,2 (Ω) and its norm || · || m are written as H m (Ω) and || · || m,2 . In addition, the dual space for For simplicity, we define u and v such that they belong to the same finite element space X. We define two continuous bilinear forms a(·, ·) and d(·, ·) and a trilinear form b(u; v, w): In addition, we need to define two bilinear forms a(·, ·) and a trilinear form b(·; ·, ·) in W × W and X × W × W, respectively.
The variational form of the NC equations (3) is as follows: solving (u, p, T) ∈ X × M × W, ∀(v, q, s) ∈ X × M × W: Then, there exists a unique result for the following solution [2,3].

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 τ h , whose minimum angle θ min is bounded away from zero, N(K) and N(e) to represent the union of all elements K and all edges e of the elements in τ h , respectively. In addition, h e and h K denote, as usual, the diameters of an edge e and element K, respectively. Then, velocity-pressure-temperature finite element spaces can be constructed on τ h (Ω). First, we define the finite element space where P 1 (K) represents piecewise linear polynomials on K, and P 0 (K) represents piecewise constant polynomials on K. According to the above definition, the discrete form of the original equation is given: For the finite element approximation (3), we present some results in [2,3] as follows.
Under the assumption of Theorem 1, the finite element For convenience, we define Since we choose an unstable finite element pair P 1 -P 0 -P 1 , we consider the following penalty jump stabilization method: where Γ h is the set of inner edges of τ h , [q] e is the jump of q on edge e, and β 0 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, R h : X → X h satisfies the following elementwise error estimate (see [24,27]). For ϕ ∈ X, there is

Application to Natural Convection Equations
We define A h and D h analogously: In this way, formula (5) can be written as It is well known that Ra 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 A h (α) is poor. To overcome this problem, let α 0 ≥ α, A(α 0 ) be a more stabilized or regularized approximation of A h (α). The DCM is given as follows: First, solveũ 1 h to satisfy Next, the above form is corrected by correction steps, which makes the obtainedũ j+1 h more closely approximateũ for j = 1, 2, ..., N: With A h , α, α 0 (K) := max(α, h K ) and D h (·), the defect-correction discretization is applied to the incompressible NC equations (see [6]), where the trilinear term b(·; ·, ·) and b(·; ·, ·) are discretized by the Oseen scheme. In the calculation process, in order to produce better results, D h (·) 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:

Error Analysis
Similarly, the term ||R h || L(X, is a consistency error, which is a higher-order term, so we only need to estimate ||( Theorem 3. Suppose that (u, p, T) and (u j+1 , p j+1 , T j+1 ) are the solutions of (2) and (5), respectively: where RF(ũ, α) = A(α) + RD(ũ), RD(·) is Lipschitz continuous in some region aroundũ, and diag(u) is a diagonal matrix with diagonal elements from the vector u.
Proof of Theorem 1. First, we consider ||( We integrate by parts over each element K ∈ τ h (Ω) and denote the collection of edges of τ h (Ω) in the interior of Ω by E h . Using the Cauchy-Schwarz inequality on each element K and edge e: We assumeφ = (ϕ, q, s) This completes the proof for Theorem 1.
Note that this amounts to replacing ||RF(ũ, The residual term is r j+1 , and C = C(C 1 , C 2 , θ min (τ h (Ω))) is a computable constant. In addition, let S h be the set of all interior sides in Ω, and for any piecewise constant q, let [q] e = q| K + e − q| K − e denote their jumps on the side e ∈ S h , where K + e and K −

Recovery Error Estimator
In this section, we construct a recovery error estimator and analyze its properties. Assuming that (u h , p h , T h ) is the numerical result, we consider the pseudo-stress tensor σ := Pr∇u − pI + κdiag(∇T) and its finite element approximation σ j+1 h The main idea of the Zienkienwicz-Zhu estimator is to restore the discontinuous finite element gradient to a continuous recovery term (see [37]).
In τ h , N, N h , φ i and ω v are respectively defined as the set of all vertices in τ h , the set of all vertices within Ω in τ h and the base function of i (∀i ∈ N) that shares a collection of all units of a vertex (ω i := supp φ i ).
Combined with the superconvergent piece recovery technique in the literature [37], the piecewise constant tensor σ j+1 h is restored to continuous, and a recovery pseudo-stress tensor is obtained.
For any vertex i ∈ N and its patch ω i , the following is defined: where |K| 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 η K 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 ||[σ [σ · n e ] := σ| K + e · n e + σ| K − e · n e , [σ · ι e ] := σ| K + e · ι e + σ| K − e · ι e denote the jumps of the normal and tangential component of σ and G(σ j+1 h ), there are two grid-independent constants C 1 and C 2 .  h ] e || S h error estimation can be obtained.
Proof of Lemma 2. According to the Brezis-Gallonet inequality in [37] and the method in [37], the conclusion can be easily obtained.

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.  (2) and (5), respectively, Proof of Theorem 2. According to (12), with the Cauchy-Schwarz inequality, we have According to the definition of B h , the approximate solution (u At the same time F h (ũ j+1 h ),φ X×M×W = 0. With the above lemma, we obtain Theorem 4 is proved.

Effectiveness Analysis
From Lemma 1, the following can be obtained:  (2) and (5), respectively. There is a normal number C that is not related to the grid: Proof of Lemma 3. From the definition of σ h , the following can be obtained: Specific details can be found in [37].

||[∇T
Proof of Theorem 3. Combined with the above estimates, Lemma 1 and Lemma 3 are applied.
The proof of validity is complete.

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: • DOF := the number of degrees of freedom of triangulation τ j h ; denotes the relative value of L 2 norm; denotes the relative value of global recovery-based estimator; • e r rate := 2 log(e j+1 r /e j r ) log(DOF j /DOF j+1 ) denotes the convergence rate of the error; • η r rate := 2 log(η j+1 r /η j r ) log(DOF j /DOF j+1 ) denotes the convergence rate of the error; • I e f f := η r e r denotes effectivity index for the global recovery-based estimator η r , where the validity index I e f f denotes the ratio of the estimator error to the true error η r /e r , and if I e f f tends to 1, it means that our estimator error is asymptotically equivalent to the true error, thus verifying our validity and reliability.

Smooth True Solution
The purpose of the first example is to solve a smooth true solution problem in the Ω = [0, 1] 2 region and verify our method, which is effective for the smooth true solution. We define the true solution u = (u 1 , u 2 ), pressure p and temperature T as follows:  Tables 1 and 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 η r for similar degrees of freedom with different Ra. From the table, we can see that the error results without the defect-correction method when Ra = 10 4 are very large, and when Ra is greater than 10 4 , the results are not counted. Our method, on the other hand, is still relatively stable at Ra = 10 5 . It can be seen that the posteriori error estimator based on the defect-correction method is applicable. I e f f tends to 1, which indicates that our estimators and real errors are effective and progressive, indicating that our method is effective.

L-Shape Domain Problem
The second example is a flow problem in the L-shape domain Ω = (−1, 1) 2 − [0, 1] 2 , with Pr = 0.71, Ra = 100 and 1000, κ = 1. The NC equations (3) is determined by exact velocities u 1 and u 2 , 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 (0.1, 0.1) and along the line y = −1.05, 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). Tables 4 and 5 show the numerical results for uniform/adaptive refinements with the recovery-based estimator η r .  According to in Figure 2c, we can see that near the singular point (0.1, 0.1) of u and the line of y = −1.05 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 Tables 4-6, our adaptive method can obtain better accuracy with less vertex information, which is better than the uniform grid. For example, when Ra = 100, 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 Ra = 1000, 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 Ra becomes larger, our method calculates results that are superior to FEM.

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 Ω = [0, 1] 2 , with Pr = 0.71 and κ = 1, and the boundary definition is as follows: the left boundary and the bottom are T = 0, and the bottom edge is ∂T ∂n = 0. The rest is T = 4y(1 − y), and the speed is the 0 Dirichlet condition. Figure 4 shows the initial mesh and the adaptive encryption result grid. Figures 5 and 6 are the velocity streamline and isothermal chart for different Rayleigh numbers: 10 3 , 10 4 and 10 5 . Figures 7 and 8 present the results of the streamline of the velocity and the isotherm diagram calculated without the defect-correction method at different Ra numbers, and it can be seen that as Ra increases, the calculated results become more unstable [38]. In Figure 9, we also show the vertical velocity distribution at x 2 = 0.5 and the horizontal velocity distribution at x 1 = 0.5, 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].
Tables 7 and 8 present the peak values of vertical velocity at x 2 = 0.5 and horizontal velocity at x 1 = 0.5 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].

Bernard Convection Problem
In this experiment, we consider the domain Ω = [0, 5] × [0, 1] with the forcing f = 0 and γ = 0. Figure 10 displays the initial and boundary conditions for velocity u and temperature T: u = 0 on ∂Ω, ∂T/∂n = 0 on the lateral boundaries, with a fixed temperature difference between the top and bottom boundaries.
The results of numerical streamlines and numerical temperature are shown in Figures 11 and 12 with Pr = 0.71, κ = 1 and Ra = 10 4 . The results coincide with the phenomenon described in the literature [38], indicating that our adaptive DCM is suitable.

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 P 1 -P 0 -P 1 . 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.