Study on Accuracy of the High-Resolution Schemes

Abstract

The high-resolution (HR) schemes have been widely used as they can achieve the numerical solution without oscillation and artificial diffusion, especially for convection-dominated problems. However, there still have arguments about the order of accuracy of HR schemes, especially about the extreme value of the solution. In this paper, it is proved that any HR scheme designed in the NVD diagram has second-order accuracy when its combined segments totally locate in the BAIR region. In other words, it has been verified in our study that the segments, which have low-order accuracy when independently employed, have at least second-order accuracy when locate in BAIR region by analysis of two implementation methods of HR scheme and also a number of numerical examples. Meanwhile Richardson extrapolation has been used to estimate the order of accuracy of HR schemes which achieve the same conclusion.

1. Introduction

How to discretize the convective term in the governing equation, in order to acquire bounded and accurate numerical solution in the convection-dominated problems, has become one of the research keynotes in the field of computational fluid dynamics and numerical heat transfer. Some certain lower-order schemes, such as the first-order upwind (FUD) scheme, can present unconditionally stability in the numerical simulation procedure. However, when such lower-order schemes are used to discretize the convective term, numerical false diffusion would be caused, thus it can degrade the accuracy of the results, especially in the computational domain with sharp gradients, which means that these lower-order schemes cannot illustrate the physical problem perfectly due to the false diffusion. Therefore, some schemes, such as QUICK (Quadratic Upstream Interpolation for Convective Kinematics) scheme and CD (second-order central difference) scheme, are popular due to their higher-order accuracy. When these higher-order schemes are adopted, they can guarantee the order of accuracy of results, that means we can obtain more accurate numerical solution; but it would cause unphysical oscillation or overshoot/undershoot (unbounded solution) when grid Peclet number exceeds certain limit or a sharp change of profile exists in the computational domain, which lead to the low-accuracy of result, especially for the CD scheme which would result in obvious oscillation, overshoot/undershoot and instability.

Therefore, researchers in this field make attempts to develop the discrete scheme which can guarantee the boundedness and accuracy of the numerical solution. The high-resolution (HR) schemes proposed in recent years can overcome the above drawbacks which means HR schemes cannot only eliminate nonphysical oscillation and overshoot/undershoot, but also weaken the false diffusion and guarantee the stability of the solution procedure. Early, HR schemes are developed on the basis of total variation diminishing (TVD) proposed by Harten [1]; after that Gaskell and Lau proposed normalized variable formulation (NVF) and the convective boundedness criterion (CBC) [2], and then the definitions of HR schemes can be illustrated in normalized variable diagram (NVD) [3] after normalized. Some HR schemes developed in recent years, such as HOAB [4], HLPA [5], MINMOD [1], MUSCL [6], SMART [2] and STOIC [7], are found to obtain the results which are accurate with the absence of unphysical oscillation and artificial diffusion.

The research on the boundedness and accuracy of HR schemes has become a key point. It has been generally considered that the HR scheme has at least second-order accuracy on the whole when its characteristic line passes through node Q [8] in NVD (see Figure 2), that is when HR scheme is used in the physical problem without extreme value of solution, it has at least second-order accuracy. However, when there is the extreme point in the physical problem, researchers hold different views on the accuracy of HR schemes: Leonard indicated that the characteristic line which is near the extreme value of solution is of first-order accuracy, and the adoption of HR scheme at the extreme point would not degrade the order of accuracy of the overall algorithm [3]; Gaskell and Lau pointed out that the adoption of HR scheme at the extreme point of solution is a second-order approximation only through the adoption of the form containing QUICK scheme and correction term [2]; Sweby held the view that at the extreme point of the solution the second-order accuracy must be lost [9]; Wei et al. concluded that the scheme possesses at least second-order accuracy in any condition when its characteristic line locates within BAIR (Boundedness accuracy and interpolative reasonableness) region [4]. Thus it can be seen that the discussions about the order of accuracy of HR schemes, especially about that at the extreme point of solution, have not reached an agreement yet. In this paper the accuracy of HR schemes are discussed to investigate whether the HR schemes locating in the BAIR region, particularly at the extreme point, have at least second-order accuracy or not. The layout is as follows. In Section 2, the order of accuracy of HR schemes has been analyzed. And a number of numerical examples are shown in Section 3 to verify the analysis in Section 2. Finally, related conclusions are given in Section 4.

2. Analysis of the Accuracy of HR Schemes

2.1. Normalized Variable Formulation

The value of cell face is involved when a finite volume method is used to discretize the governing equation. Generally the value of cell face is acquired by the interpolation of the values of its nearby nodes. In a one-dimensional uniform grid system, the symbol U, C, and D are referred as the upstream, central and downstream node, respectively, and the symbol frepresents the cell face (see Figure 1). Then the interpolation of cell face values can be written as: ϕ_f = F(ϕ_U,ϕ_C,ϕ_D). Introducing the following normalized variable:

{\tilde{ϕ}}_{} = \frac{ϕ - ϕ_{U}}{ϕ_{D} - ϕ_{U}} .

(1)

Note that ${\tilde{ϕ}}_{U} = 0$ and ${\tilde{ϕ}}_{D} = 1$ , and the function of cell face value is only associated with ${\tilde{ϕ}}_{C}$ after normalized, that is: ${\tilde{ϕ}}_{f} = F ({\tilde{ϕ}}_{C})$ . The difference scheme can be defined by normalized variable, for example SUD scheme and CD scheme can be represented respectively as: $(1) {\tilde{ϕ}}_{f} = 1.5 {\tilde{ϕ}}_{C}$ for SUD scheme; $(2) {\tilde{ϕ}}_{f} = 0.5 + 0.5 {\tilde{ϕ}}_{C}$ for CD scheme. The normalized expressions of several common used HR schemes (MINMOD [1], MUSCL [6], SMART [2], STOIC [7], HOAB [4] and HLPA [5]) are illustrated in Table 1.

Table 1:

Normalized expressions of HR schemes.

HR Schemes	${\tilde{ϕ}}_{f}$	${\tilde{ϕ}}_{C}$

HLPA	${\tilde{ϕ}}_{C} (2 - {\tilde{ϕ}}_{C})$	${\tilde{ϕ}}_{C} = 0.5$
	${\tilde{ϕ}}_{C} (2 - {\tilde{ϕ}}_{C})$	$0 \leq {\tilde{ϕ}}_{C} \leq 1$
	${\tilde{ϕ}}_{C}$	Else

MUSCL	$2 {\tilde{ϕ}}_{C}$	$0 \leq {\tilde{ϕ}}_{C} \leq \frac{1}{4}$
	$\frac{1}{4} + {\tilde{ϕ}}_{C}$	$\frac{1}{4} \leq {\tilde{ϕ}}_{C} \leq \frac{3}{4}$
	1	$\frac{3}{4} \leq {\tilde{ϕ}}_{C} \leq 1$
	${\tilde{ϕ}}_{C}$	Else

HOAB	$3.5 {\tilde{ϕ}}_{C}$	$0 \leq {\tilde{ϕ}}_{C} \leq \frac{1}{6}$
	$\frac{1}{2} + \frac{1}{2} {\tilde{ϕ}}_{C}$	$\frac{1}{6} \leq {\tilde{ϕ}}_{C} \leq \frac{1}{2}$
	$\frac{1}{4} + {\tilde{ϕ}}_{C}$	$\frac{1}{2} \leq {\tilde{ϕ}}_{C} \leq \frac{3}{4}$
	1	$\frac{3}{4} \leq {\tilde{ϕ}}_{C} \leq 1$
	${\tilde{ϕ}}_{C}$	Else

MINMOD	$\frac{3}{2} {\tilde{ϕ}}_{C}$	$0 \leq {\tilde{ϕ}}_{C} \leq \frac{1}{2}$
	$\frac{1}{2} + \frac{1}{2} {\tilde{ϕ}}_{C}$	$\frac{1}{2} \leq {\tilde{ϕ}}_{C} \leq 1$
	${\tilde{ϕ}}_{C}$	Else

SMART	$3 {\tilde{ϕ}}_{C}$	$0 \leq {\tilde{ϕ}}_{C} \leq \frac{1}{6}$
	$\frac{3}{8} + \frac{3}{4} {\tilde{ϕ}}_{C}$	$\frac{1}{6} \leq {\tilde{ϕ}}_{C} \leq \frac{5}{6}$
	1	$\frac{5}{6} \leq {\tilde{ϕ}}_{C} \leq 1$
	${\tilde{ϕ}}_{C}$	Else

STOIC	$3 {\tilde{ϕ}}_{C}$	$0 \leq {\tilde{ϕ}}_{C} \leq \frac{1}{5}$
	$\frac{1}{2} + \frac{1}{2} {\tilde{ϕ}}_{C}$	$\frac{1}{5} \leq {\tilde{ϕ}}_{C} \leq \frac{1}{2}$
	$\frac{3}{8} + \frac{3}{4} {\tilde{ϕ}}_{C}$	$\frac{1}{2} \leq {\tilde{ϕ}}_{C} \leq \frac{5}{6}$
	1	$\frac{5}{6} \leq {\tilde{ϕ}}_{C} \leq 1$
	${\tilde{ϕ}}_{C}$	Else

Figure 1:

Uniform grid system.

Figure 2:

NVD of SUD scheme and CD scheme.

Taking ${\tilde{ϕ}}_{C}$ as an independent variable and ${\tilde{ϕ}}_{f}$ as a dependent variable, the normalized scheme can be expressed in the form of characteristic line, and the diagram that can illustrate the relationship between ${\tilde{ϕ}}_{C}$ and ${\tilde{ϕ}}_{f}$ is called normalized variable diagram (NVD). The definitions of SUD scheme and CD scheme in the NVD are shown in Figure 2. Here for brevity, the region where ${\tilde{ϕ}}_{C} \in [0,1]$ in NVD is named as inner domain, and ${\tilde{ϕ}}_{C} \notin [0,1]$ as outer domain, and according to CBC, the HR schemes are all presented in the form of FUD scheme in the outer domain of the NVD, that is ${\tilde{ϕ}}_{f} = {\tilde{ϕ}}_{C}$ .

2.2. BAIR Condition

As SUD and CD schemes, the HR schemes can be represented by normalized variable as well, which means that the definition of HR schemes can be illustrated in NVD. The BAIR condition proposed by Wei et al. [4] refers to the region in NVD that satisfies the following conditions:

CBC proposed by Gaskell and Lau [2].

The line of a scheme passing through O(0, 0), Q(0.5, 0.75), and P(1, 1) in the NVD.

Interpolative reasonableness.

Thus the region that satisfies the BAIR condition is the shaded area shown in Figure 3. From Figure 3, it can be seen that the HR composite schemes satisfying BAIR condition own the following properties: (1) the line of a scheme is located between the lines of SUD scheme and CD scheme in the inner domain; (2) the line happens to coincide with that of FUD scheme in the outer domain; (3) any HR scheme would pass through node O, Q, and P mentioned above. The characteristic lines of several common-used HR schemes, such as MINMOD [1], HLPA, and SMART, all fall into the BAIR region.

Figure 3:

The region of BAIR.

2.3. Interpolation Method

As mentioned above, cell face value can be obtained by the interpolation of adjacent node values via various difference schemes. From the first two properties of HR schemes satisfying BAIR condition, it is known that any difference scheme falling into this region can be written as a combination of SUD scheme and CD scheme (as SCSD scheme proposed in reference [10], here we name it the rewritten scheme for short), then the cell face value is:

ϕ_{f} = α ϕ_{SUD} + (1 - α) ϕ_{CD} .

(2)

Considering that ϕ_SUD = 1.5ϕ_C − 0.5ϕ_U, ϕ_CD = 0.5ϕ_C + 0.5ϕ_D, (2) can be normalized as:

{\tilde{ϕ}}_{f} = α {\tilde{ϕ}}_{SUD} + (1 - α) {\tilde{ϕ}}_{CD} .

(3)

The first and the second term on the right side of (3) represent the normalized SUD scheme and CD scheme respectively, in which ${\tilde{ϕ}}_{SUD} = 1.5 {\tilde{ϕ}}_{C}$ and ${\tilde{ϕ}}_{CD} = 0.5 + 0.5 {\tilde{ϕ}}_{C}$ . The value of coefficient α represents the proportion of SUD scheme which depends on the type of scheme. When α = 1.0 it means that SUD scheme is employed for the cell-face value interpolation, and α = 0.0 means CD scheme is employed.

For schemes satisfying the BAIR condition, the characteristic lines lie between the lines of SUD scheme and CD scheme, which means the value of ${\tilde{ϕ}}_{f}$ is between the values of the above two schemes with the same ${\tilde{ϕ}}_{C}$ . Thus it can obtain α ∈ [0, 1], which can guarantee at least second-order accuracy of the rewritten HR schemes. In addition, the truncation error term of the SUD scheme according to the Taylor series expansion is − (3/4)Δx²(∂²ϕ/∂x²), and that of CD scheme is (1/4)Δx²(∂²ϕ/∂x²), thus that of the rewritten scheme is:

\begin{matrix} o (Δ x^{2}) = [- α \frac{3}{4} + (1 - α) \frac{1}{4}] Δ x^{2} \frac{\partial^{2} ϕ}{\partial x^{2}} \\ = (\frac{1}{4} - α) Δ x^{2} \frac{\partial^{2} ϕ}{\partial x^{2}} . \end{matrix}

(4)

On the basis of value range of α and the above truncation error term shown in (4), the value range of truncation error of the rewritten scheme is [− 0.75Δx²,0.25Δx²], indicating that the rewritten scheme is at least second-order accurate, while it is third-order accurate particularly when α = 0.25.

2.3.1. Derivation of α

For brevity, QUICK scheme is taken as an example here to present interpolation using (2), in order to explain the process of getting α value. The original form of QUICK scheme is:

ϕ_{f} = \frac{1}{8} (3 ϕ_{D} + 6 ϕ_{C} - ϕ_{U}) .

(5)

Substitute (5) into (2), we have

α = \frac{1}{4} .

(6)

Then, the rewritten form of QUICK scheme is obtained,

ϕ_{f} = \frac{1}{4} ϕ_{SUD} + \frac{3}{4} ϕ_{CD} .

(7)

Equation (7) is also the expression of cell face value when the QUICK scheme is used to discretize governing equation.

HR schemes usually consist of several difference schemes, and any difference scheme satisfying BAIR condition can be rewritten in the form like (2), that is to say they can be written in a form combining with α, SUD scheme and CD scheme together. Here in order to explain the process of derivation of α for HR schemes, we take SMART as example, which is:

{\tilde{ϕ}}_{f} = {\begin{cases} 3 {\tilde{ϕ}}_{C} & 0 \leq {\tilde{ϕ}}_{C} \leq \frac{1}{6} \\ \frac{3}{8} + \frac{3}{4} {\tilde{ϕ}}_{C} & \frac{1}{6} \leq {\tilde{ϕ}}_{C} \leq \frac{5}{6} \\ 1 & \frac{5}{6} \leq {\tilde{ϕ}}_{C} \leq 1 \\ {\tilde{ϕ}}_{C} & else . \end{cases}

(8)

It can be seen that SMART is composed of four difference schemes, among which every scheme can be expressed in the form shown in (3), and different intervals of ${\tilde{ϕ}}_{C}$ have a corresponding α, respectively. Here one of the SMART segments is taken as an example to derive the α, when $0 < {\tilde{ϕ}}_{C} \leq 1 / 6$ we have:

3 {\tilde{ϕ}}_{C} = α {\tilde{ϕ}}_{SUD} + (1 - α) {\tilde{ϕ}}_{CD},

(9)

where ${\tilde{ϕ}}_{SUD} = 1.5 {\tilde{ϕ}}_{C}$ and ${\tilde{ϕ}}_{CD} = 0.5 + 0.5 {\tilde{ϕ}}_{C}$ . After simplifying (9) we get:

α = \frac{5 {\tilde{ϕ}}_{C} - 1}{2 {\tilde{ϕ}}_{C} - 1} .

(10)

Then the other segments of SMART can be obtained in the same way, from which the SMRAT scheme can be rewritten as:

{\tilde{ϕ}}_{f} = {\begin{cases} \frac{5 {\tilde{ϕ}}_{C} - 1}{2 {\tilde{ϕ}}_{C} - 1} {\tilde{ϕ}}_{SUD} + (1 - \frac{5 {\tilde{ϕ}}_{C} - 1}{2 {\tilde{ϕ}}_{C} - 1}) {\tilde{ϕ}}_{CD} & 0 \leq {\tilde{ϕ}}_{C} \leq \frac{1}{6} \\ \frac{1}{4} {\tilde{ϕ}}_{SUD} + \frac{3}{4} {\tilde{ϕ}}_{CD} & \frac{1}{6} \leq {\tilde{ϕ}}_{C} \leq \frac{5}{6} \\ \frac{1 - {\tilde{ϕ}}_{C}}{2 {\tilde{ϕ}}_{C} - 1} {\tilde{ϕ}}_{SUD} + (1 - \frac{1 - {\tilde{ϕ}}_{C}}{2 {\tilde{ϕ}}_{C} - 1}) {\tilde{ϕ}}_{CD} & \frac{5}{6} \leq {\tilde{ϕ}}_{C} \leq 1 \\ \frac{{\tilde{ϕ}}_{C} - 1}{2 {\tilde{ϕ}}_{C} - 1} {\tilde{ϕ}}_{SUD} + (1 - \frac{{\tilde{ϕ}}_{C} - 1}{2 {\tilde{ϕ}}_{C} - 1}) {\tilde{ϕ}}_{CD} & else . \end{cases}

(11)

2.3.2. Values of α for Different HR Schemes

Through above analysis, it is found that HR schemes lying in the BAIR region can be described by expressions containing coefficient α. The values of α in different intervals of several common used HR schemes (MINMOD [1], MUSCL [6], SMART [2], STOIC [7], HOAB [4] and HLPA [5]) are illustrated in Table 2.

Table 2:

Value of α.

HR Schemes	α	${\tilde{ϕ}}_{C}$

HLPA	0.5	${\tilde{ϕ}}_{C} = 0.5$
	$\frac{{\tilde{ϕ}}_{C} (3 - 2 {\tilde{ϕ}}_{C}) - 1}{2 {\tilde{ϕ}}_{C} - 1}$	$0 \leq {\tilde{ϕ}}_{C} \leq 1$
	$\frac{{\tilde{ϕ}}_{C} - 1}{2 {\tilde{ϕ}}_{C} - 1}$	Else

MUSCL	$\frac{3 {\tilde{ϕ}}_{C} - 1}{2 {\tilde{ϕ}}_{C} - 1}$	$0 \leq {\tilde{ϕ}}_{C} \leq \frac{1}{4}$
	$\frac{1}{2}$	$\frac{1}{4} \leq {\tilde{ϕ}}_{C} \leq \frac{3}{4}$
	$\frac{1 - {\tilde{ϕ}}_{C}}{2 {\tilde{ϕ}}_{C} - 1}$	$\frac{3}{4} \leq {\tilde{ϕ}}_{C} \leq 1$
	$\frac{{\tilde{ϕ}}_{C} - 1}{2 {\tilde{ϕ}}_{C} - 1}$	Else

HOAB	$\frac{6 {\tilde{ϕ}}_{C} - 1}{2 {\tilde{ϕ}}_{C} - 1}$	$0 \leq {\tilde{ϕ}}_{C} \leq \frac{1}{6}$
	0.0	$\frac{1}{6} \leq {\tilde{ϕ}}_{C} \leq \frac{1}{2}$
	$\frac{1}{2}$	$\frac{1}{2} \leq {\tilde{ϕ}}_{C} \leq \frac{3}{4}$
	$\frac{1 - {\tilde{ϕ}}_{C}}{2 {\tilde{ϕ}}_{C} - 1}$	$\frac{3}{4} \leq {\tilde{ϕ}}_{C} \leq 1$
	$\frac{{\tilde{ϕ}}_{C} - 1}{2 {\tilde{ϕ}}_{C} - 1}$	Else

MINMOD	1.0	$0 \leq {\tilde{ϕ}}_{C} \leq \frac{1}{2}$
	0.0	$\frac{1}{2} \leq {\tilde{ϕ}}_{C} \leq 1$
	$\frac{{\tilde{ϕ}}_{C} - 1}{2 {\tilde{ϕ}}_{C} - 1}$	Else

SMART	$\frac{5 {\tilde{ϕ}}_{C} - 1}{2 {\tilde{ϕ}}_{C} - 1}$	$0 \leq {\tilde{ϕ}}_{C} \leq \frac{1}{6}$
	$\frac{1}{4}$	$\frac{1}{6} \leq {\tilde{ϕ}}_{C} \leq \frac{5}{6}$
	$\frac{1 - {\tilde{ϕ}}_{C}}{2 {\tilde{ϕ}}_{C} - 1}$	$\frac{5}{6} \leq {\tilde{ϕ}}_{C} \leq 1$
	$\frac{{\tilde{ϕ}}_{C} - 1}{2 {\tilde{ϕ}}_{C} - 1}$	Else

STOIC	$\frac{5 {\tilde{ϕ}}_{C} - 1}{2 {\tilde{ϕ}}_{C} - 1}$	$0 \leq {\tilde{ϕ}}_{C} \leq \frac{1}{5}$
	0.0	$\frac{1}{5} \leq {\tilde{ϕ}}_{C} \leq \frac{1}{2}$
	$\frac{1}{4}$	$\frac{1}{2} \leq {\tilde{ϕ}}_{C} \leq \frac{5}{6}$
	$\frac{1 - {\tilde{ϕ}}_{C}}{2 {\tilde{ϕ}}_{C} - 1}$	$\frac{5}{6} \leq {\tilde{ϕ}}_{C} \leq 1$
	$\frac{{\tilde{ϕ}}_{C} - 1}{2 {\tilde{ϕ}}_{C} - 1}$	Else

In order to express the value of α clearly, the relationship between α and ${\tilde{ϕ}}_{C}$ in the inner domain of the HR schemes shown in Table 2 is presented in Figure 4. α of different schemes in the outer domain can be presented in an identical way, and from $α = ({\tilde{ϕ}}_{C} - 1) / (2 {\tilde{ϕ}}_{C} - 1)$ we know that when ${\tilde{ϕ}}_{C} \to + \infty$ or ${\tilde{ϕ}}_{C} \to - \infty$ we have α → 0.5, which means when ${\tilde{ϕ}}_{C} \in (1, + \infty)$ we have α ∈ [0, 0.5) and when ${\tilde{ϕ}}_{C} \in (- \infty, 0)$ we get α ∈ (0.5, 1] (see Figure 5). Therefore it can be verified again that α ∈ [0, 1], and HR schemes have at least a second-order accuracy after rewritten.

Figure 4:

The value of α for different HR schemes.

Figure 5:

The value of α in the outer domain.

2.4. The Implementation of HR Schemes

The implementation of HR schemes is achieved through discretizing the governing equation using HR schemes. The difference between the implementations of HR schemes lies in the way to interpolate cell face value. There are two kinds of implementation methods of rewritten HR schemes (here the interval $0 \leq {\tilde{ϕ}}_{C} \leq 1 / 6$ of SMART is taken as an example):

(1) Direct Implementation Method. In direct implementation method, HR schemes are employed in the interpolation directly. As α of rewritten SMART in the above interval is $α = (5 {\tilde{ϕ}}_{C} - 1) / (2 {\tilde{ϕ}}_{C} - 1)$ , the interpolation of cell face value can be presented as:

\begin{array}{l} ϕ_{f} = α ϕ_{SUD} + (1 - α) ϕ_{CD} \\ = \frac{5 {\tilde{ϕ}}_{C} - 1}{2 {\tilde{ϕ}}_{C} - 1} (\frac{3}{2} ϕ_{C} - \frac{1}{2} ϕ_{U}) - \frac{3 {\tilde{ϕ}}_{C}}{2 {\tilde{ϕ}}_{C} - 1} (\frac{1}{2} ϕ_{C} + \frac{1}{2} ϕ_{D}) \\ = A_{C} ϕ_{C} + A_{D} ϕ_{D} + A_{U} ϕ_{U}, \end{array}

(12)

where $A_{C} = 3 (4 {\tilde{ϕ}}_{C} - 1) / [2 (2 {\tilde{ϕ}}_{C} - 1)]$ , $A_{D} = - 3 {\tilde{ϕ}}_{C} / [2 (2 {\tilde{ϕ}}_{C} - 1)]$ and $A_{U} = (1 - 5 {\tilde{ϕ}}_{C}) / [2 (2 {\tilde{ϕ}}_{C} - 1)]$ . And the solution mentioned in reference [2] is employed, that is the SUR (successive under-relaxation) iterative method, in which A_C, A_D and A_U containing α are regarded as the coefficients of discretized equation. Thus, it is equivalent to use the scheme that possesses at least second-order accuracy in the discrezition of convective term.

(2) Deferred Correction (DC) Method. In DC method, HR schemes have been introduced into the source term of discretized equation in the form of correction term, that is to say, DC method is adopted to discretize the convective term. Then the interpolation of cell face value can be presented as:

ϕ_{f} = ϕ_{C}^{new} + {[α ϕ_{SUD} + (1 - α) ϕ_{CD} - ϕ_{C}]}^{old},

(13)

where superscripts new and old represent numerical results of the present and last iteration layer, respectively, and the second term on the right side of (13) is the correction term which can be put into the source term of discretized equation. Since the first term on the right side of (13) is independent of the second term, the form of (13) is equivalent to use FUD scheme to discretize convective term. Substitute $α = (5 {\tilde{ϕ}}_{C} - 1) / (2 {\tilde{ϕ}}_{C} - 1)$ into (13), we get:

ϕ_{f} = ϕ_{C}^{new} + {[(3 ϕ_{C} - 2 ϕ_{U}) - ϕ_{C}]}^{old} .

(14)

If DC method is introduced into the original HR scheme, the expression of interpolation is:

ϕ_{f} = ϕ_{C}^{new} + {[F ({\tilde{ϕ}}_{C}) \cdot (ϕ_{D} - ϕ_{U}) + ϕ_{U} - ϕ_{C}]}^{old} .

(15)

As is known, in this interval the original SMART is ${\tilde{ϕ}}_{f} = 3 {\tilde{ϕ}}_{C}$ , and when introducing it to (15) the form of interpolation same with (13) can be obtained. It indicates that when DC method has been introduced into the discretization of convective term, the same expression of convective term can be gained in the discretized equation by the rewritten and the original HR schemes respectively, which means the rewritten and original HR schemes are of equivalence. However, considering the superiority of the latter in computation efficiency, the expression in (15) is employed when DC method is introduced. In addition, the introduction of DC method can improve the stability of solution only, and has no influence on the accuracy. Thus when the above two implementation methods of HR schemes are employed respectively and the solution reaches final convergence, the same results can be obtained, that is they are equivalent. Using the first implementation method, every segment of composite schemes, especially the FUD scheme segment in the outer domain, is replaced by schemes that possess second-order accuracy as α ∈ [0, 1]. On the other hand, employing the second implementation method, the interpolation of cell face value is ϕ_f = ϕ_C^new in the outer domain, and it has been verified that the two kinds of implementation methods are equivalent and the segment of FUD scheme is replaced by second-order accurate schemes, that means the segment of FUD scheme is of second-order accuracy, so is the inner domain. Therefore, the employment of (15) has second-order accuracy, in other words, every segment of HR schemes that satisfies the BAIR condition is of at least second-order accuracy, though it may be indicated by Taylor series expansion that some segments of HR scheme, such as ${\tilde{ϕ}}_{f} = 3 {\tilde{ϕ}}_{C}$ , have low-order accuracy when independently employed. The above statements are proved in Section 3 by a number of numerical examples.

3. Test Examples

In this section, the accuracy of HR scheme is explained in detail. First, the above two implementation methods of HR schemes are employed in the following test examples to verify the accuracy of HR schemes: (1) a pure convection problem of a stepwise profile in an oblique uniform velocity field (it is called 2D pure convection problem for short); (2) lid-driven flow in a square cavity with the given temperature boundary condition (it is called lid-driven flow in a square cavity for short). Second, by applying the Richardson extrapolation [11] in the second problem, the order of accuracy of HR schemes has been given. The following HR schemes have been adopted in both two problems: HOAB, HLPA, MINMOD, MUSCL, SMART and STOIC. The sum of the residual of discretized equation is:

RES = \sum | a_{P} ϕ_{P} - (\sum_{NB = E, W, N, S} a_{NB} ϕ_{NB} + b) |,

(16)

which can be referred as the convergence criterion, thus the convergence will be reached when the sum of the residual error is smaller than the set value (depending on specific problem).

3.1. Explanations of Test Examples

3.1.1. Problem 1: 2D Pure Convection Problem

This problem is shown in Figure 6. In an oblique uniform velocity field, the configuration of transportation scalar is a stepwise profile; the flow angle is 30°; the boundary condition of transportation scalar is also presented in Figure 6. A uniform grid system with 22 × 22 grids is used, and the convergence criterion is: RES < 10⁻⁵. In Figure 7, the scalar profiles at the mid-horizontal plane (y = 0.5) predicted by HLPA, MINMOD and SMART (representative and typical of the six HR schemes) are presented for comparison. Here, Solution I and Solution II represent the numerical solution of the first (see (12)) and the second (see (15)) implementation method, respectively. Figure 7 indicates that the results of two implementation methods agree well with each other. To analyze the results of the above two implementation methods quantitatively, the definitions of average deviation and the maximum deviation are given at first, which are:

\begin{matrix} \bar{err} = \frac{\sum_{all nodes} | ϕ_{DC} - ϕ_{C - S} |}{N^{2}}, \\ er r_{\max} = \max (| ϕ_{DC} - ϕ_{C - S} |), \end{matrix}

(17)

where ϕ_DC represents the results obtained by DC method shown in (15); ϕ_{C − S} represents the result obtained by the employment of rewritten scheme shown in (12), that is the original HR schemes have been replaced by the combination of two schemes of second-order accuracy, then the discretized equation contains α which represents the proportion of SUD scheme, and the coefficient matrix is a nine-diagonal matrix, and SUR iterative method has been used in solution; $\bar{err}$ refers to average deviation and err_max refers to maximum deviation; N stands for the number of grid nodes. The derivation of the above two implementation methods of different HR schemes (still taking HLPA, MINMOD and SMRAT as examples) is presented in Table 3. When the computation reached the final convergence, the lines that are carried out by the two implementation methods of the same HR scheme are totally same in the same location of the computational domain, and every segment of composite schemes has been carried out by the two implementation methods during the solution. Moreover, the deviation in Table 3 shows that the difference between two implementation methods are less than 5 × 10⁻⁵, which can be convincingly explained that such differences are caused by round-off error. Then the equivalence of two implementation methods can be verified by the same results they acquired.

Table 3:

Deviation of two implementation methods for Problem 1.

	Deviation	HLPA	MINMOD	SMART

Θ	$\bar{err}$	8.27 × 10⁻⁷	8.71 × 10⁻⁷	1.42 × 10⁻⁶
Θ	err_max	7.17 × 10⁻⁵	3.34 × 10⁻⁵	4.30 × 10⁻⁵

Figure 6:

A 2D pure convection problem.

Figure 7:

Comparison for scalar profiles at the mid-horizontal plane (y = 0.5) of two different implementation methods of HR schemes for problem 1.

3.1.2. Problem 2: Lid-Driven Flow in a Square Cavity

This problem is illustrated in Figure 8, from which we know the velocity of top moving wall is u_top; Re = 1000; T_h and T_l represent the given high-temperature condition and the low-temperature condition respectively. A uniform grid with 66 × 66 cells is adopted, and the diffusion term is discretized by CD scheme. The convergence criterion is RES < 10⁻¹⁰ for velocity field and RES < 10⁻¹³ for temperature field. Equations (18) and (19) are the dimensionless governing equation of this physical problem:

\begin{matrix} \frac{\partial U}{\partial X} + \frac{\partial V}{\partial Y} = 0, \\ \frac{\partial (U U)}{\partial X} + \frac{\partial (U V)}{\partial Y} = - \frac{\partial P}{\partial X} + \frac{1}{Re} (\frac{\partial^{2} U}{\partial X^{2}} + \frac{\partial^{2} U}{\partial Y^{2}}), \end{matrix}

(18)

\begin{matrix} \frac{\partial (V U)}{\partial X} + \frac{\partial (V V)}{\partial Y} = - \frac{\partial P}{\partial Y} + \frac{1}{Re} (\frac{\partial^{2} V}{\partial X^{2}} + \frac{\partial^{2} V}{\partial Y^{2}}), \\ \frac{\partial (U Θ)}{\partial X} + \frac{\partial (V Θ)}{\partial Y} = \frac{1}{Re \cdot Pr} (\frac{\partial^{2} Θ}{\partial X^{2}} + \frac{\partial^{2} Θ}{\partial Y^{2}}) . \end{matrix}

(19)

In this paper, coupling solution is used to calculate the variables U and V in (18), then distribution of temperature field is obtained based on the convergent velocity field, which is equivalent to computing a convection-diffusion problem with governing equation (19). The comparison of the velocity and temperature distributions for the two implementation methods is shown in Figure 9 (here we take HLPA, SMART and MINMOD as examples); Table 4 shows deviation of the above two implementation methods in order to compare the results quantitatively.

Table 4:

Deviation of two implementation methods for Problem 2.

	Deviation	HLPA	MINMOD	SMART

U	$\bar{err}$	1.12 × 10⁻⁷	3.75 × 10⁻⁸	1.30 × 10⁻⁷
U	err_max	4.02 × 10⁻⁷	1.39 × 10⁻⁷	4.63 × 10⁻⁷

V	$\bar{err}$	1.06 × 10⁻⁷	3.61 × 10⁻⁸	1.23 × 10⁻⁷
V	err_max	3.94 × 10⁻⁷	1.38 × 10⁻⁷	4.57 × 10⁻⁷

Θ	$\bar{err}$	2.41 × 10⁻⁹	8.15 × 10⁻¹⁰	2.78 × 10⁻⁹
Θ	err_max	5.75 × 10⁻⁹	1.94 × 10⁻⁹	6.63 × 10⁻⁹

Figure 8:

Lid-driven flow in a square cavity.

Figure 9:

Comparison of results for two implementation methods.

Although the convergence processes of two implementation methods are different, the calculation shows that the profiles obtained by the two implementation methods of HR schemes are the same in any location of the computational domain when the calculation reaches convergence, and every segment of composite schemes has been carried out by both implementation methods during the solution process, for example the percentage of the adoption of the segment ${\tilde{ϕ}}_{f} = {\tilde{ϕ}}_{C}$ in the outer domain is at least 1.6%. It can be seen from comparison and deviation that the results of the two implementation methods agree well. Similarly, it can be indicated that the deviation of the two implementation methods in this problem is caused by round-off error. In addition, it can be seen from (18) that velocity field is calculated by non-linear equations and that of temperature field is solved by linear equation, and the results of velocity field and temperature field by two implementation methods agree well respectively, so it indicates that two implementation methods of HR schemes can get same results no matter whether it is a non-linear problem or a linear problem. Therefore, the application of the two implementation methods has the same effect.

From the above two test examples it can be concluded that the adoption of DC method and the combination of SUD scheme and CD scheme to interpolate the cell face value respectively can obtain the same results as expected, which means the results obtained by the two implementation methods of HR schemes are of the same order of accuracy, and thus the two implementation methods of HR schemes are proved to be equivalent. The expression used in the first implementation method is (12), which is expressed by the two schemes that have second-order accuracy, thus the numerical solution is certain to have second-order accuracy when (12) is used. Moreover, every segment of composite scheme is presented in the form of (12), and it can be known from the equivalence of two implementation methods that the segments of composite schemes that don't pass through node Q, such as ${\tilde{ϕ}}_{f} = {\tilde{ϕ}}_{C}$ and ${\tilde{ϕ}}_{f} = 3 {\tilde{ϕ}}_{C}$ , and the line of which is curved rather than straight, such as ${\tilde{ϕ}}_{f} = {\tilde{ϕ}}_{C} (2 - {\tilde{ϕ}}_{C})$ , have second-order accuracy as well.

3.2. Richardson Extrapolation

Through the two different implementation methods of HR schemes employed to obtain cell face value, it has been verified in the last section that the two implementation methods are equivalent, and then verified that HR schemes satisfying the BAIR condition have second-order accuracy on the whole. In this section, order of accuracy of HR schemes is explained by Richardson extrapolation [11] combined with numerical experiments. The Richardson extrapolation can estimate the order of truncation error of the discretized equation, and its idea is:

Numerical computation would be carried out on three sets of grids with different densities, and the ratio between two adjacent grid systems is 3.

The value of ϕ in the same location of three sets of grids should be recorded.

The order of accuracy of truncation error is.

n = \frac{\ln [(φ_{3 h} - φ_{9 h}) / (φ_{h} - φ_{3 h})]}{\ln (3)},

(20)

where ϕ_h, ϕ_3h, and ϕ_9h represent the numerical solution of the grids with step size of h, 3h and 9h respectively; the value of n is the order of accuracy, for example n = 2 means that the order of truncation error is second. Three sets of grids are shown in Figure 10, node P is the communal node of three sets of grids, and symbols e_9h, e_3h and e_h represent the east face of control volume of node P in coarse grid, middle grid and refined grid respectively; s_9h, w_9h and n_9h represent the south face, west face and north face of control volume of node P in the coarse grid respectively. The numerical computation is carried out in three sets of grids respectively, if the schemes used to interpolate the value of east face of control volume of communal node P in coarse grid, middle grid and refined grid respectively are the same, (20) can be used to solve n. Meanwhile variables ϕ_h,ϕ_3h, and ϕ_9h in (20) refer to the value of ϕ on the communal node, such as node P, calculated in three sets of grids respectively.

Figure 10:

Three sets of uniform grid systems.

The Richardson extrapolation will be employed in the test example 2 to estimate the order of accuracy of HR schemes, and DC method will be used for the computation. The control volume number of three sets of grids are 98 × 98 (the coarse grid), 290 × 290 (the middle grid) and 866 × 866 (the refined grid) respectively. And convergence criteria for both the coarse grid and the middle grid are RES < 10⁻¹⁰, and that for the refined grid is RES < 10⁻¹³. The results of the estimation for the overall order of accuracy of schemes (the value of n) are shown in Table 5, and those of the segments of composite schemes are presented in Table 6, where the data in parentheses represent the adoption percentage of every segment of composite scheme used to interpolate the value of e_9h of control volume when computing ϕ in coarse grids. And what shown in Table 5 are the averaged results extrapolated by all nodes, which means it is neglected whether the schemes used to interpolative east face of communal node in different girds are the same or not. In Table 6, in order to satisfy the requirement (the value of ϕ in the same location) of Richardson extrapolation, for those segments whose available nodes are not enough, such as ${\tilde{ϕ}}_{f} = 2 {\tilde{ϕ}}_{C}, {\tilde{ϕ}}_{f} = 1$ and ${\tilde{ϕ}}_{f} = {\tilde{ϕ}}_{C}$ of MUSCL, aggregation measures have been taken, that is the Richardson extrapolation has been adopted to estimate the overall order of accuracy of these segments. Moreover, in Table 6, in the process of temperature field calculation using HLPA in three sets of grids respectively, the communal nodes used in the expression ${\tilde{ϕ}}_{f} = {\tilde{ϕ}}_{C}$ to interpolate the cell face value of e_9h, e_3h and e_h simultaneously are not available, thus the Richardson extrapolation can't be employed for this segment; the estimation result of order of accuracy for the QUICK segment of composite scheme is second-order, that is because only the scheme used to interpolate the east face of control volume has been concerned. However, the value of ϕ used in extrapolation also depends on the way of interpolation for south face, west face and north face, for example, if the other three faces values have been interpolated by schemes which have second-order accuracy, then it can greatly affect the accuracy of numerical solution, which can lead to the estimation result of order of accuracy for QUICK isn't equivalent to third-order. In addition, the diffusion term is discretized by CD scheme, which is second-order accurate and can also affect the value of n.

Table 5:

Order of accuracy of HR schemes.

	HLPA	HOAB	MINMOD	MUSCL	SMART	STOIC

U	2.16	2.01	1.97	2.09	1.97	1.94
V	2.10	2.03	1.98	2.07	1.98	1.97
Θ	2.14	2.06	1.98	2.09	1.97	1.97

Table 6:

Order of accuracy and percentage of segments of HR schemes.

HR Scheme	${\tilde{ϕ}}_{f}$	n (%)
HR Scheme	${\tilde{ϕ}}_{f}$	U	V	Θ

MINMOD	${\tilde{ϕ}}_{f} = {\tilde{ϕ}}_{C}$	2.31 (2.41%)	1.88 (2.04%)	1.97 (1.57%)
	${\tilde{ϕ}}_{f} = 1.5 {\tilde{ϕ}}_{C}$	1.95 (49.05%)	1.97 (51.92%)	1.99 (52.41%)
	${\tilde{ϕ}}_{f} = 0.5 + 0.5 {\tilde{ϕ}}_{C}$	1.97 (48.54%)	1.98 (46.04%)	1.98 (46.02%)

MUSCL	${\begin{matrix} {\tilde{ϕ}}_{f} = 2 {\tilde{ϕ}}_{C} \\ {\tilde{ϕ}}_{f} = 1 \\ {\tilde{ϕ}}_{f} = {\tilde{ϕ}}_{C} \end{matrix}$	2.05 (5.81%)	2.14 (4.23%)	2.08 (3.28%)
MUSCL	${\tilde{ϕ}}_{f} = 0.25 + {\tilde{ϕ}}_{C}$	2.10 (94.82%)	2.06 (95.77%)	2.09 (96.72%)

SMART	${\begin{matrix} {\tilde{ϕ}}_{f} = 3 {\tilde{ϕ}}_{C} \\ {\tilde{ϕ}}_{f} = 1 \\ {\tilde{ϕ}}_{f} = {\tilde{ϕ}}_{C} \end{matrix}$	2.13 (4.22%)	1.99 (3.42%)	1.97 (2.35%)
SMART	${\tilde{ϕ}}_{f} = 0.375 + 0.75 {\tilde{ϕ}}_{C}$	1.97 (95.78%)	1.97 (96.58%)	1.97 (97.65%)

STOIC	${\begin{matrix} {\tilde{ϕ}}_{f} = 3 {\tilde{ϕ}}_{C} \\ {\tilde{ϕ}}_{f} = 1 \\ {\tilde{ϕ}}_{f} = {\tilde{ϕ}}_{C} \end{matrix}$	2.06 (4.57%)	1.94 (3.78%)	1.98 (2.5%)
	${\tilde{ϕ}}_{f} = 0.5 + 0.5 {\tilde{ϕ}}_{C}$	1.93 (47.74%)	1.99 (43.72%)	1.98 (51.75%)
	${\tilde{ϕ}}_{f} = 0.375 + 0.75 {\tilde{ϕ}}_{C}$	1.95 (47.69%)	1.96 (52.50%)	1.95 (45.75%)

HOAB	${\begin{matrix} {\tilde{ϕ}}_{f} = 3.5 {\tilde{ϕ}}_{C} \\ {\tilde{ϕ}}_{f} = 1 \\ {\tilde{ϕ}}_{f} = {\tilde{ϕ}}_{C} \end{matrix}$	2.16 (4.70%)	2.05 (3.71%)	2.06 (2.86%)
	${\tilde{ϕ}}_{f} = 0.5 + 0.5 {\tilde{ϕ}}_{C}$	1.98 (48.18%)	2.06 (43.78%)	2.06 (51.84%)
	${\tilde{ϕ}}_{f} = 0.25 + {\tilde{ϕ}}_{C}$	2.04 (47.12%)	2.01 (52.51%)	2.06 (45.30%)

HLPA	${\tilde{ϕ}}_{f} = {\tilde{ϕ}}_{C} (2 - {\tilde{ϕ}}_{C})$	2.16 (97.59%)	2.09 (97.95%)	2.14 (98.41%)
HLPA	${\tilde{ϕ}}_{f} = {\tilde{ϕ}}_{C}$	2.38 (2.41%)	1.95 (2.05%)	— (1.59%)

From Table 5, it is noted that the values of n representing order of accuracy of different HR schemes all equal 2 approximately; Table 6 shows that the values of n for each segment of HR schemes are also 2 approximately. Then it can be concluded that the above kinds of HR schemes possess not only second-order accuracy on the whole, but also have at least second-order accuracy of each segment, especially in the situation that ${\tilde{ϕ}}_{f} = {\tilde{ϕ}}_{C}$ or ${\tilde{ϕ}}_{f} = 3 {\tilde{ϕ}}_{C}$ and the type of characteristic line is curved.

4. Conclusions

The numerical solutions gained by DC method and the combination of two schemes possessing second-order accuracy respectively have been compared when HR schemes is used to discretize the governing equation. And it has been verified that HR schemes satisfying the BAIR condition have at least second-order accuracy (no matter whether the lines of which are straight or curved). In other words, any scheme of which the line in the NVD lies between the lines of two second-order schemes (SUD scheme and CD scheme) in the inner domain and happens to coincide with FUD scheme in the outer domain dose possesses at least second-order accuracy. The following points should be valued particularly:

For ${\tilde{ϕ}}_{f} = {\tilde{ϕ}}_{C}$ segment, only when its line locates in BAIR region which means the line lies between SUD scheme and CD scheme, the adoption of the expression ${\tilde{ϕ}}_{f} = {\tilde{ϕ}}_{C}$ has second-order accuracy, and in other region it still has first-order accuracy.

Those segments locating in BAIR region without passing through node Q, such as ${\tilde{ϕ}}_{f} = 3 {\tilde{ϕ}}_{C}$ , can be verified by Taylor expansion that they have lower-order accuracy when employed independently. However, when they are parts of composite schemes, these segments are second-order accurate as they lie in BAIR region.

The scheme whose line is curved and passes through node Q in BAIR region is also second-order accurate as the scheme whose line is straight.

Meanwhile, the application of Richardson extrapolation has explained and also proved the above statements. In addition, DC method is recommended to implement HR schemes because of its simplicity and superiority in computation efficiency.

Footnotes

Nomenclature

Conflict of Interests

The authors declare that there is no conflict of interests regarding the publication of this paper.

Acknowledgments

The study is supported by the National Science Foundation of China (nos. 51325603 and 51134006).

References

Harten

, “High resolution schemes for hyperbolic conservation laws,” Journal of Computational Physics, vol. 49, no. 3, pp. 357–393, 1983.

Gaskell

P. H.

and Lau

A. K. C.

, “Curvature-compensated convective transport: SMART a new boundedness-preserving transport algorithm,” International Journal for Numerical Methods in Fluids, vol. 8, no. 6, pp. 617–641, 1988.

Leonard

B. P.

, “Simple high-accuracy resolution program for convective modeling of discontinuities,” International Journal for Numerical Methods in Fluids, vol. 8, no. 10, pp. 1291–1318, 1988.

Wei

J.-J.

Tao

W.-Q.

Kawaguchi

, and Wang

H.-S.

, “A new high-order-accurate and bounded scheme for incompressible flow,” Numerical Heat Transfer B, vol. 43, no. 1, pp. 19–41, 2003.

Zhu

, “A low-diffusive and oscillation-free convection scheme,” Communications in Applied Numerical Methods, vol. 7, no. 3, pp. 225–232, 1991.

van Leer

, “Towards the ultimate conservative difference scheme. II. Monotonicity and conservation combined in a second-order scheme,” Journal of Computational Physics, vol. 14, no. 4, pp. 361–370, 1974.

Darwish

M. S.

, “A new high-resolution scheme based on the normalized variable formulation,” Numerical Heat Transfer B, vol. 24, no. 3, pp. 353–373, 1993.

Leonard

B. P.

, “The ULTIMATE conservative difference scheme applied to unsteady one-dimensional advection,” Computer Methods in Applied Mechanics and Engineering, vol. 88, no. 1, pp. 17–74, 1991.

Sweby

P. K.

, “High resolution schemes using flux limiters for hyperbolic conservation laws,” SIAM Journal on Numerical Analysis, vol. 21, no. 5, pp. 995–1011, 1984.

10.

M. J.

Tao

W. Q.

, and Wang

S. J.

, “Stability-controllable second-order upwind difference scheme for convection term,” Journal of Thermal Science, vol. 7, no. 2, pp. 119–130, 1998.

11.

Richardson

L. F.

, “The approximate arithmetical solution by finite differences of physical problems involving differential equations, with an application to the stresses in a Masonry dam,” Philosophical Transactions of the Royal Society A, vol. 210, pp. 307–357, 1911.