Monotonic unbounded schemes transformer (MUST): an approach to remove undershoots and overshoots in family of unbounded schemes using finite-volume method

Abstract

Purpose

This paper presents a Monotonic Unbounded Schemes Transformer (MUST) approach to bound/monotonize (remove undershoots and overshoots) unbounded spatial differencing schemes automatically, and naturally. Automatically means the approach (1) captures the critical cell Peclet number when an unbounded scheme starts to produce physically unrealistic solution automatically, and (2) removes the undershoots and overshoots as part of the formulation without requiring human interventions. Naturally implies, all the terms in the discretization equation of the unbounded spatial differencing scheme are retained.

Design/methodology/approach

The authors do not formulate new higher-order scheme. MUST transforms an unbounded higher-order scheme into a bounded higher-order scheme.

Findings

The solutions obtained with MUST are identical to those without MUST when the cell Peclet number is smaller than the critical cell Peclet number. For cell Peclet numbers larger than the critical cell Peclet numbers, MUST sets the nodal values to the limiter value which can be derived for the problem at-hand. The authors propose a way to derive the limiter value. The authors tested MUST on the central differencing scheme, the second-order upwind differencing scheme and the QUICK differencing scheme. In all cases tested, MUST is able to (1) capture the critical cell Peclet numbers; the exact locations when overshoots and undershoots occur, and (2) limit the nodal value to the value of the limiter values. These are achieved by retaining all the discretization terms of the respective differencing schemes naturally and accomplished automatically as part of the discretization process. The authors demonstrated MUST using one-dimensional problems. Results for a two-dimensional convection–diffusion problem are shown in Appendix to show generality of MUST.

Originality/value

The authors present an original approach to convert any unbounded scheme to bounded scheme while retaining all the terms in the original discretization equation.

Citation

Yap, Y.F. and Chai, J.C. (2024), "Monotonic unbounded schemes transformer (MUST): an approach to remove undershoots and overshoots in family of unbounded schemes using finite-volume method", International Journal of Numerical Methods for Heat & Fluid Flow, Vol. 34 No. 11, pp. 4049-4084. https://doi.org/10.1108/HFF-04-2024-0293

Nomenclature

a_E

= coefficient of the east neighbor;

a_P

= coefficient of the nodal point;

a_W

= coefficient of the west neighbor;

C, D

= positive constants;

k_u

= undershoot indicator;

k_s

= sign(ϕ_E – ϕ_limit);

P

= Peclet number;

S

= source term;

u

= fluid velocity; and

x

= spatial location.

Symbols

Δx

= width of a control volume;

Γ

= diffusion coefficient;

ϕ

= generic dependent variable;

ϕ_desired

= desired value of ϕ;

ϕ_P,ns

= ϕ_P without source terms;

ϕ_P,s

= ϕ_P with source terms;

ϕP*

= the most current value of ϕ_P;

ρ

= fluid density;

ξ

= combined variable;

ξ′

= undershoots variable; and

ξ″

= overshoots variable.

1. Introduction

In the modeling of the convection–diffusion equation using the finite-volume method, profile assumptions are needed to evaluate the dependent variable at the two boundaries of a control volume. When the central difference (CD) scheme is used, the value of the dependent variable at an interface is interpolated linearly from the two neighboring control volumes as shown in Figure 1a. This approach is stable and accurate when the cell Peclet numbers are less than 2. Once the cell Peclet number exceeds 2, physically unrealistic undershoots and overshoots are predicted. To overcome this, one-sided bias (toward the upstream neighbor) to the interpolation is introduced. The simplest one is called the upwind scheme [Figure 1(b)] where the value of the dependent variable at an interface is set to the value of the dependent variable of the upstream neighbor. In multidimensional flows where the velocity is not aligned with the grid, this approach though stable for all cell Peclet numbers, suffers from false diffusion. To reduce false diffusion, higher-order schemes were formulated. These include, but are not limited to the second-order upwind difference (SOUD) scheme [Figure 1(c)] and the QUICK scheme [Figure 1(d)] (Leonard, 1979). Unlike the CD scheme and the upwind scheme, the SOUD scheme involves one more upstream neighbor, namely, the WW neighbor as shown in Figure 1(c). The QUICK scheme includes two more neighbors one upstream (WW) neighbor and one downstream (p) neighbor as depicted in Figure 1(d). These two higher-order upwind family difference schemes reduce false diffusion, but admit physically unrealistic solutions which appear as undershoots and overshoots. In the CD scheme, the point where physically unrealistic solution starts is known. It happens when the cell Peclet number is larger than 2 regardless of the values of the neighboring dependent values. Therefore, it is possible to devise an approach which switch to another stable scheme when the cell Peclet number exceed 2. Spalding (Spalding, 1972) presented such a scheme and it is called the hybrid scheme.

Multiples studies comparing various spatial difference schemes for have been reported in the literature (de Vahl Davis and Mallinson, 1976; Beier et al., 1983; Leschziner, 1980; Leschziner and Rodi, 1981; Pollard and Siu, 1982; Runchal, 1972; Han et al., 1981; Shyy et al., 1992; Vanka, 1987). De Vahl Davis and Mallinson (de Vahl Davis and Mallinson, 1976) evaluated the upwind and CD schemes using a driven cavity problem. Their main concern was to evaluate the effects of the two schemes on false diffusion. Beier et al. (Beier et al., 1983) examined five spatial difference scheme including the upwind, CD and QUICK schemes on their abilities to model recirculating flows. They concluded that QUICK and CD schemes are more accurate than the upwind scheme. Leschziner (Leschziner, 1980) evaluated the hybrid scheme, a hybrid central and skew-upwind scheme and the QUICK scheme on four problems with and without recirculating flows. He concluded that the hybrid central and skew-upwind scheme and the QUICK scheme performed better than the hybrid scheme. However, these two more accurate schemes suffer from boundness or undershoots and overshoots problems. Leschziner and Rodi (Leschziner and Rodi, 1981) extended Leschziner’s work (Leschziner, 1980) to study unconfined turbulent recirculating flows using the same three difference schemes. The same conclusions were drawn. Pollard and Siu (Pollard and Siu, 1982) tested the upwind scheme, the hybrid scheme and the QUICK scheme on three flow geometries, namely, flow between parallel plates, axisymmetric flow with sudden expansion and driven cavity flow. The QUICK scheme was found to be equal or more accurate than the upwind and the hybrid schemes. Runchal (Runchal, 1972) carried out a systematic performance tests comparing the upwind difference scheme, the CD scheme and the hybrid difference scheme. As expected, the upwind difference scheme and the hybrid difference scheme were unconditionally stable but less accurate than the CD scheme when the CD scheme produced physically realistic solutions. Han et al. (Han et al., 1981) studied the hybrid scheme and the QUICK scheme for both laminar and turbulent recirculating flows. They concluded that QUICK scheme, though not unconditionally stable is more accurate than the hybrid scheme. Shyy et al. (Shyy et al., 1992) carried out numerical studies on a driven cavity problem using different implementations of the second-order upwind scheme. They did not notice oscillations in the solutions and mentioned that a critical cell Reynolds (Peclet) number similar to the CD scheme above which wiggles are noticed may not be valid for the SOUD scheme. Vanka (Vanka, 1987) implemented three forms of the SOUD scheme and compared their solutions with the hybrid scheme. They found that while the hybrid scheme is stable while the SOUD scheme is unstable above a certain Reynolds number. He concluded that if the reasons for the overshoots can be identified and removed, SOUD scheme’s second-order accuracy will be beneficial.

In an effort to ensure boundness of higher-order schemes, the flux limiters (FL) and normalized variables (NV) approaches have been proposed. Piperno and Depeyre (Piperno and Depeyre, 1998) proposed criteria for the design of limiters. Waterson and Deconinck (Waterson and Deconinck, 2007) proposed a unified design principles to formulated bounded higher-order spatial differencing schemes. The NV schemes examined by the authors were converted to the FL approach to facilitate direct comparisons. Some higher-order schemes include, but are not limited to, MINMOD (Harten, 1983), CLAM (van Leer, 1974), MUSCL (van Leer, 1979), SMART (Gaskell and Lau, 1988), WENO (Liu et al., 1994). Shu (Shu, 2009) gave an exhaustive review on WENO of various order. WENO was used to solve parabolic problems with degenerate diffusion (Arbogast et al., 2019). Recent developments include but are not limited to, a sixth-order WENO scheme (Abedian, 2023), WENO for nonlinear degenerate parabolic problems (Jiang, 2021), a fourth-order scheme for natural convection in enclosures (Balam and Gupta, 2020), comparisons of different spatial differencing schemes for one-dimensional multiphase flows in porous media (Moshiri and Manzari, 2019).

The common traits in the conclusions of the above studies are (1) upwind scheme and hybrid scheme are stable but less accurate, (2) CD scheme, SOUD scheme and QUICK scheme are more accurate but suffer from boundness problems, not unconditionally stable, (3) that a critical cell Reynolds (Peclet) number similar to the CD scheme above which wiggles are noticed may not be valid for the SOUD scheme and (4) if the reasons for the overshoots can be identified and removed, SOUD scheme’s second-order accuracy will be beneficial.

This article will address the last two items. We propose an approach to remove undershoots and overshoots automatically and naturally. Automatically means our approach capture the critical Peclet number when a higher-order upwind scheme starts to produce physically unrealistic solution automatically as part of our formulation and does not require human interventions. Naturally implies we retain all the terms in the discretization equation of the higher-order upwind scheme.

In the next three subsections, we describe the monotonic behavior of the convection-–diffusion equation. This is followed by a discussion on the problem by examining the discretization equation of the CD scheme which does not ensure monotonic variations. The objectives of the article are then stated. The outline of the article is presented to conclude this section.

1.1 Steady, convection–diffusion equation without source [1]

In this section, we will re-examine the way the discretization schemes are implemented to discretize the steady, convection–diffusion equation without source which, for one-dimensional situations, can be written as:

(1) d(ρuϕ)dx=ddx(Γdϕdx)

where ρ is the fluid density, u is the fluid velocity, ϕ is a generic dependent variable, x is the spatial location and Γ is the diffusion coefficient. Equation (1) is subjected to the following boundary conditions:

(2a)  x=0→ϕ=ϕ0

(2b)  x=L→ϕ=ϕL

Using equation (1) with equation (2) as the boundary conditions, the exact solution can be written as:

(3)  ϕ−ϕ0ϕL−ϕ0=exp(Px/L)−1exp(P)−1

where the Peclet number P is defined as:

(4)  P ≡ ρuLΓ

From the exact solution given by equation (3), the dependent variable varies exponentially between the two end points. As such, the exact solution does not contain undershoots (ϕ < ϕ_min) or overshoots (ϕ > ϕ_max) where ϕ_min = min(ϕ₀, ϕ_L), and ϕ_max = max(ϕ₀, ϕ_L).

Exponential Scheme: The exponential scheme is based on the exact solution [equation (3)], and for constant fluid properties the familiar algebraic discretization equation can be written as:

(5)  aPϕP=aEϕE+aWϕW

where

(6a) aE=Pexp⁡(P)−1

(6b) aW=Pexp⁡(P)exp⁡(P)−1

(6c) aP=aE+aW

The subscripts in equations (5) and (6) refer to the control volume of interest p, its west neighbor W and its east neighbor E as shown in Figure 2. As the exponential scheme [equation (6)] is obtained from the exact solution, it does not admit undershoots and overshoots in its solutions. In another words, the nodal value of the dependent variable is between the values of the upstream neighbor and the downstream neighbor or min(ϕ_E, ϕ_W) ≤ ϕ_P ≤ max(ϕ_E, ϕ_W).

The problem: Unbounded differencing schemes including, but are not limited to, the CD scheme, the SOUD scheme and the QUICK scheme, admit undershoots and overshoots from the convection–diffusion parts of the governing equation (equation (1)). We will use the CD scheme to demonstrate this behavior. Using the CD scheme, and for constant properties, a discretization equation for equation (1) is given by equation (5) and the coefficients of the CD scheme can be written as:

(7a) aE=1−P2

(7b) aW=1+P2

(7c) aP=aE+aW

The values of ϕ_P for P = 3, with different ϕ_E, and ϕ_W obtained using equations (5) and (7) are shown in Table 1. For the CD scheme, it is obvious when the cell Peclet number is great than 2, the nodal value ϕ_P overshoots or undershoots the value of its upstream neighbor ϕ_W. It is possible to use the hybrid scheme (Spalding, 1972) to switch to the upwind scheme when P > 2 to avoid undershoots or overshoots.

In this simple example, it can be seen that a two-step process is needed to eliminate undershoots and overshoots. These are (1) detecting the critical cell P where undershoots or overshoots begin, and (2) switching to another scheme which does not produce undershoots or overshoots.

There are two issues associated with this approach of detecting the critical cell Peclet number a priori. First, for the CD scheme, the demarcation of when to switch is known and it is when |P| ≥ 2. However, this demarcation is not so obvious in other unbounded difference schemes. There is a second issue which will be discussed in the next section.

1.2 Steady, convection–diffusion equation with source

In this section, we consider the steady, convection–diffusion equation with source which, for one-dimensional situations can be written as:

(8) d(ρuϕ)dx=ddx(Γdϕdx)+S

where S is the source term. The discretization equation for equation (8) using the exponential scheme is:

(9) aPϕP=aEϕE+aWϕW+b

where

(10a) aE=Pexp⁡(P)−1

(10b) aW=Pexp⁡(P)exp⁡(P)−1

(10c) aP=aE+aW

(10d) b=S(Δx)2Γ

where Δx is the width of a control volume shown in Figure 2.

Let us examine the effects of the source term. For this discussion, we will consider a situation where u > 0, and ϕ_W < ϕ_E. The exponential scheme without source [equation (10a–c)] produces nodal ϕ_P between ϕ_W and ϕ_E. For ease of discussions, we will call this nodal value without source as ϕ_P,ns. Now consider a situation where there is a positive source, the exponential scheme will predict nodal ϕ_P with positive source of ϕ_W ≤ ϕ_P,ns < ϕ_P,s.

The corresponding discretization equation for the CD scheme is the same as given in equation (9). The coefficients for the CD scheme are:

(11a) aE=1−P2 

(11b) aW=1+P2

(11c) aP=aE+aW

(11d) b=S(Δx)2Γ

It is commonly accepted that in the CD scheme, the critical cell Peclet number, where the nodal ϕ_P values fall below or exceed the upstream neighbor values, is 2. However, this critical cell Peclet number is known only when there is no physical source term. When there are physical sources, this demarcation is not known a priori even for the CD scheme! [2] We will again consider a situation where u > 0, and ϕ_W < ϕ_E. To better understand this statement, let us consider a situation where the convection–diffusion part of the discretization equation of the CD scheme is operating in the undershoot regime or ϕ_P,ns ≤ ϕ_W. For this to happen, the cell Peclet number must be greater than 2. Now, we will introduce a positive source. Using equation (11), there are two scenarios of interest. First, the positive source is smaller than ϕ_P,ns from the convection–diffusion part of the discretization equation. This leads to a nodal ϕ which is ϕ_P,ns < ϕ_P,s < ϕ_W. Second, the positive source is larger than ϕ_P,ns from the convection–diffusion part of the discretization equation. This leads to a nodal ϕ where ϕ_P,ns < ϕ_W < ϕ_P,s. In this example, with the same undershoot from the discretization of the convection–diffusion part, depending on the strength of the source, the nodal ϕ_P value with positive source can be smaller than or larger than the upstream neighbor value. In both cases, due to the undershoot created by the discretization of the convection–diffusion part of the governing equation, ϕ_P,CD < ϕ_P,EXP. For the CD scheme, the most conservative approach is to switch to a bounded scheme when |P| > 2 as done in the hybrid scheme (Spalding, 1972) for all situation.

When there is no source term, for the CD scheme, there is a critical cell Peclet number for the convection–diffusion part of the discretization equation where, ϕ_P falls below or exceed the value of the upstream neighbor. This, however, is not the case in other unbounded upwind scheme; even for the convection–diffusion part of the discretization equation. Source terms makes the problem more difficult to analyze.

In summary, when an unbounded spatial differencing scheme is used in the presence of physical sources, and there are no overshoots (for u > 0, ϕ_P > ϕ_W) or undershoots (for u > 0, ϕ_P < ϕ_W) in the solutions, it is not possible to conclude that there are no overshoots or undershoots in the convection–diffusion part of the governing equation.

As a result, there is a need to device an approach that (1) automatically detects when undershoots and overshoots, (2) eliminates the detected undershoots and overshoots and (3) continues to use all the terms of the same higher-order scheme.

1.3 Objectives of the article

This article presents an approach to (1) define and determine the undershoots and overshoots values, (2) automatically captures the exact locations where undershoots and overshoots occur, (3) limits the nodal value to a value determined in (1) and (4) continues to use the same higher-order scheme. The final discretization equation has the same order as the original discretization equation as all terms in the original discretization equation are retained. These are achieved automatically (without ad hoc intervention), and naturally (retaining all discretization terms associated with the governing equation).

1.4 Outline of the article

The remainder of this article is divided into four (4) sections. The always-positive variable treatment of Patankar (Patankar, 1980) is discussed in the next section. This approach ensures an always-positive variable remains positive throughout the solutions process. In another word, it limits the value of a variable ϕ to zero and above or ϕ = max(ϕ, 0). The extension of this approach to limit a variable ϕ to any arbitrary minimum value ϕ_min or ϕ = max(ϕ, ϕ_min) is then described. In the same section, we review Patankar’s approach to set the value of an internal control volume to any desired value. The incorporation of our approach using extensions of Patankar’s always-positive variable and internal value specification to the CD scheme, the SOUD scheme and the QUICK scheme is presented in the following section. This approach transforms an unbounded spatial difference scheme to a bounded spatial difference scheme while maintaining all the discretization terms of the unbounded scheme and is called Monotonic Upwind Scheme Transformer (MUST). The next section demonstrated the capabilities of MUST using the CD, SOUD and QUICK schemes. This article then concludes with some concluding remarks. We presented initial extension of MUST to a simple two-dimensional problem in Appendix to show the generality of MUST.

2. Treatment of always-positive variables and its generalization

In this section, we will describe Patankar’s (Patankar, 1980) always-positive variable treatment. This is followed by an extension of the approach to eliminate undershoots. We will conclude this section by extending the undershoots elimination approach to eliminate overshoots.

2.1 Patankar’s always-positive variable treatment

As discussed in Patankar (Patankar, 1980), this approach ensures an always-positive variable remains always positive in the presence of a negative source term. Consider a discretization equation with positive coefficients, and a negative source (C is positive) given in equation (12):

(12) aPϕP=aEϕE+aWϕW−C

where

(13) aP=aE+aW

We can re-write equation (12) as:

(14) (aP+CϕP*)ϕP=aEϕE+aWϕW

where ϕP* is the most current value of ϕ_P. Note that once converged, ϕP*=ϕP, equation (14) reverts back to the original discretization equation, namely, equation (12). As there are no negative source terms on the right-side of equation (14), and all coefficients are positive, ϕ_P will remains always positive. By this simple re-arrangement of the discretization equation, we have ensured ϕ_P to be always positive. However, this is not exactly what we need in mitigating undershoots or overshoots in our spatial differencing schemes. In the spatial differencing schemes, we want an approach that will eliminate (1) undershoots by ensuring ϕ_P ≥ ϕ_desired, and (2) overshoots by guaranteeing ϕ_P ≤ ϕ_desired. We will now extend the always-positive variable treatment to eliminate undershoots and overshoots. A combined formulation will then be presented.

2.2 Elimination of undershoots

In subsection 2.1, Patankar presented a way to ensure an always-positive variable to remains always positive. Another way to view this is Patankar described an approach to eliminates undershoots of a variable ϕ to values below zero. In our application, we want to keep a variable from falling below a desired value ϕ_desired. We will follow a three-step approach by (1) define a new variable ξ′ ≡ ϕ – ϕ_desired, (2) re-write the discretization equation in-terms of ξ′ and (3) apply Patankar’s approach to ensure ξ′ ≥ 0, which in-turns means ϕ ≥ ϕ_desired. We define a new variable ξ′ as:

(15) ξ′≡ ϕ−ϕdesired

where ϕ_desired is the desire value we want to set the value of the dependent variable. We will add the same term to both sides of equation (12) as:

(16) aP(ϕP−ϕdesired)=aEϕE+aWϕW−C−aPϕdesired

Using equation (13), equation (16) can be written as:

(17) aP(ϕP−ϕdesired)=aE(ϕE−ϕdesired)+aW(ϕW−ϕdesired)−C

Using equation (15), equation (17) becomes:

(18) aPξP′=aEξE′+aWξW′−C

Now, we can ensure that ξ′ is always positive by re-arranging equation (18) as:

(19) (aP+CξP′*)ξP′=aEξE′+aWξW′

As ξP′≥0, from equation (19), we have developed a way to ensure ϕ_p ≥ ϕ_desired! This has been done by keeping all the terms of the original discretization equation, and without any approximation or alteration to the final converged value.

2.3 Elimination of overshoots

Overshoots are due to the presence of positive sources (D > 0) as written in equation (20):

(20) aPϕP=aEϕE+aWϕW+D

Using equation (15), a final discretization equation can be written as:

(21) (aP−DξP′*)ξP′=aEξE′+aWξW′

Two observations can be made. These are (1) when a_P = a_E + a_W, (aE+aW)/(aP−D/ξP*) maybe greater than 1, and (2) the a_P maybe negative. Both of these violate some of Patankar’s four basic rules (Patankar, 1980). We will not accept this type of drawbacks in our formulation. To overcome these short-comings, we will define another variable as:

(22) ξ″≡ ϕdesired−ϕ

Using equation (22) and equation (20), and after some algebraic rearrangements, we get:

(23) (aP+DξP″*)ξP″=aEξE″+aWξW″

Equation (23) ensures that ξP″≥0 or ϕ_desired ≥ ϕ_P, and eliminates overshoots. Note that the problems encountered by using equation (15) [with equation (21) as the final discretization equation] to eliminate overshoots are not found in equation (23). Now that we have formulated approaches to eliminate undershoots, and overshoots, we will present a combined formulation applicable to both situations.

2.4 Combined formulation

We will consider a situation where undershoots and overshoots (where C > 0 and D > 0) may occur as described in equation (24):

(24) aPϕP=aEϕE+aWϕW+D−C

Undershoots may occur when C > D. Using the approach described in the previous section, we can eliminate undershoots using:

(25) (aP+CξP′*)ξP′=aEξE′+aWξW′+D

As C > 0, the coefficient for ξP′ is always positive; undershoots are eliminated. We can eliminate overshoots using:

(26) (aP+DξP″*)ξP″=aEξE″+aWξW″+C

Again, all coefficients are positive, thus, eliminating overshoots. We will now define a combined variable:

(27) ξ ≡ ku(ϕ−ϕdesired)+(1−ku)(ϕdesired−ϕ)

where k_u is an undershoot indicator given by:

(28) ku ≡ sign(ϕ−ϕdesired,0)

Using equation (27), equation (25) and equation (26) can be written as:

(29) (aP−[min⁡(−ksC,0)+min⁡(ksD,0)]ξP*)ξP=aEξE+aWξW+max⁡(−ksC,0)+max⁡(ksD,0)

where

(30) ks=sign(ϕ−ϕdesired)

2.5 Patankar’s internal value specification

We now describe Patankar’s approach to set the value of an internal control volume to any desired value (ϕ_P,desired). The discretization equation with a source can be written as:

(31) (aP−SPΔx)ϕP=aEϕE+aWϕW+SCΔx

We set S_C = 10³⁰ ϕ_P,desired, and S_P = –10³⁰. Substituting these into equation (31), it reduces to

(32) (aP+1030Δx)ϕP=aEϕE+aWϕW+1030ϕP,desiredΔx

Since all other terms are finite, equation (32) reduces to:

(33) 1030ϕP=1030ϕP,desired

or

(34) ϕP=ϕP,desired

Using this simple approach, we can set the value of an internal control volume to any desired value.

2.6 Remarks

We presented an approach to ensure the dependent variable at a control volume to be larger than or equal to any arbitrary minimum value and smaller than or equal to any arbitrary maximum value. This approach eliminates undershoots and overshoots due to sources. We also describe the approach to set an internal control volume to any desired value of Patankar (Patankar, 1980).

In spatial differencing schemes, undershoots and overshoots can be the result of sources and/or negative coefficients. We present an extension of the approach presented in this section to eliminate undershoots and overshoots as the result of sources and/or negative coefficients.

3. Undershoots, overshoots and limiter value

In this section, we will define what we will call as undershoots and overshoots. If needed, other definitions can be derived. To begin, we use the solution to the steady, convection–diffusion equation without source given by equation (3). We will first describe the undershoots, and overshoots definitions. We will then formulate our undershoots, and overshoots definitions for the general situations when there are nonzero sources.

3.1 Undershoots and overshoots

For the one-dimensional control volume shown in Figure 2, when there is no flow or the Peclet number p = 0, the nodal ϕ_P value is (ϕ_W + ϕ_E)/2 as shown in Figures 3 and 4. For the situation where u > 0, as the velocity (or Peclet number) increases, the nodal ϕ_P value will approach the value of the dependent variable at the immediate upstream neighbor or ϕ_W again, as depicted in Figures 3 and 4. The exact solutions for these two situations are given by equation (3). Figure 3 shows a situation where ϕ_E > ϕ_W. In this situation, the nodal ϕ_P value decreases, and approaches ϕ_W. This ϕ_W, is indeed the correct limiting value when there are no physical source terms or S = 0. Unbounded differencing scheme can lead to undershoots (ϕ_P < ϕ_W) as shown in Figure 3. The other two curves in Figure 3 show two other possible solutions. They are approximate solutions, but there are no undershoots in these solutions. Figure 4 shows a situation where overshoots (ϕ_P > ϕ_W) is encountered. These two physically impossible solutions are called undershoots (Figure 3), and overshoots (Figure 4), respectively. From Figures 3 and 4, we want to eliminate:

(35) UndershootsϕE>ϕWOvershootsϕE<ϕW

3.2 Limiter value

Once undershoots or overshoots are detected, MUST sets the nodal value of the dependent variable to a yet-to-be-derived limiter value. This section presents one such possible limiter. Other forms of limiter can be derived.

When there are no sources, the limiter value is the value of the immediate upstream nodal. So, when u > 0, ϕ_P → ϕ_W or the limit value ϕ_limit,P = ϕ_W.

When there is a positive source term, the limiting ϕ_P value should be higher than ϕ_W. Similarly, when there is a negative source term, the limiting ϕ_P value should be lower than ϕ_W. There is no unique way to formulate the limiting value. We present a possible way to formulate this value. Other limiting values can be formulated.

This limiting ϕ_P value is obtained by performing a balance over the control volume. Performing a balance over control volume p shown in Figure 2, and for constant properties, we can write:

(36) ρuϕw+SΔx=ρuϕe

In equation (36), we assumed unit flow area. Using upwind approximation, equation (36) can be written as:

(37) ϕe=ϕlimit,P=ϕW+SΔxF

where F = ρu. A quick look at what we have achieved with equation (37) is beneficial. When the source term is zero, the limiting ϕ value reduces correctly to ϕ_W or the immediate upstream neighbor. When there is a positive source, the nodal value is limited to the upstream value plus the contribution from the source term. When the source term is negative, the nodal value is limited to the upstream value minus the contribution from the source term. Also note that diffusion is neglected in the balance equation [equation (36)]. Extending the idea in equation (35), we can define undershoots and overshoots in the presence of sources as:

(38) UndershootsϕE>ϕW+SΔxFOvershootsϕE<ϕW+SΔxF

Now that we have described our undershoots definition, overshoots definition and formulated the limit value, we will proceed to incorporate this in MUST.

Remarks: As discussed, the proposed limiter is one possible way to formulate the limiting value. In this approach, diffusion is neglected. For steady, one-dimensional problems, the nodal value approaches the value of its immediate upstream neighbor when (1) there is no source term in the control volume of interest, or (2) when the cell Peclet number is “large” when there is a finite source term in the control volume of interest. Other limiters can be derived by including diffusion, and using higher-order extrapolations of the interface value in equation (36). For two-dimensional convection-diffusion problem without source, one possible limiter is the mixing-cup values given by equation (A1). To maintain focus, this article aims at presenting the concept of MUST. The effects of different limiters on the final solutions are considered outside the scope of this article The study on the effects of the limiter is left to the explorations of interested readers.

4. Incorporation into spatial differencing schemes

In this section, we will demonstrate how to incorporate MUST into three existing higher-order differencing schemes which suffer from undershoots and overshoots. These are (1) the CD scheme, (2) the SOUD scheme and (3) the QUICK scheme. Physically unrealistic solutions in the CD scheme are due to the presence of negative coefficients. In this article, we consider a way to implement the SOUD scheme where physically unrealistic solutions are the results of the source terms from the higher-order terms of the convection-diffusion parts of the governing equation. Physically unrealistic solutions in the QUICK scheme can be due to the negative coefficients or the additional source terms due to the higher-order terms of the convection–diffusion parts of the governing equation. We elect to demonstrate MUST utilities using these three schemes as they cover all possible combinations of the reasons for physically unrealistic results of a spatial difference scheme. We will show the detailed formulations for u > 0. Formulation for u < 0 does not involve new concepts and is left to the explorations of interested readers.

4.1 The central difference scheme

The discretization equation using the CD scheme (when u > 0) is given by equation (9) and equation (11). We notice that a_E will be negative when the cell Peclet number is larger than 2. This negative a_E leads to undershoots, and overshoots. We want to develop an approach to eliminate undershoots, and overshoots when a_E < 0. We will start by removing the negative coefficient from the neighbor. To achieve this, we will add zeroes of different forms to both sides of equation (9) to write:

(39) [aE+aW+(|aE|+aE)]ϕP=aEϕE+aWϕW+b+(|aE|+aE)ϕE

At this point, it is good to know that we have not changed equation (9). Equation (39) is just a more cumbersome form of equation (9) which we re-arranged as:

(40) (|aE|+aW)ϕP=|aE|ϕE+aWϕW+b+2aE(ϕE−ϕP)

Comparing equation (9) and equation (40), it is clear that the additional terms arise when P > 2 or a_E < 0. A combined discretization equation for all Peclet numbers or all values of a_E can be written as:

(41) (|aE|+aW)ϕP=|aE|ϕE+aWϕW+b+2min(aE,0)(ϕE−ϕP)

With equation (41), we have eliminated negative coefficient a_E for the neighbor node. However, in this process, we have created an additional source term due to removal of the negative coefficient. This additional source can potentially be negative which, may leads to undershoots or positive which, may results in overshoots. Since a_E < 0, the additional source term is negative when ϕ_P < ϕ_E, and positive when ϕ_P > ϕ_E.

Undershoots: To eliminate undershoots, we must ensure ϕ_P ≥ ϕ_limit,P given by equation (37). Borrowing the idea of equation (15), we will define:

(42) ξ′ ≡ ϕ−ϕlimit,P=ϕ−ϕW−SΔxF

The main objective to define a new variable ξ′ is to ensure our dependent variable ϕ is larger than a lower limit ϕ_limit,P. Subtracting the same terms from both sides of equation (41) gives:

(43) (|aE|+aW)(ϕP−ϕW−SΔxF)=|aE|(ϕE−ϕW−SΔxF)+aW(ϕW−ϕW−SΔxF)                                           +2min(aE,0)(ϕE−ϕP)+b

Equation (43) can be written more compactly as:

(44) (|aE|+aW)ξ′P=|aE|ξ′E+aWξ′W+2min(aE,0)(ϕE−ϕP)+b

To eliminate undershoots, equation (44) can be written as:

(45)  [|aE|+aW−(min[2min(aE,0)ϕE,0]+min[−2min(aE,0)ϕP,0]+min[b,0])ξ′P*]ξ′P= |aE|ξE′+aWξW′+max[2min(aE,0)ϕE,0]+max[−2min(aE,0)ϕP,0]+max[b,0]

Using equation (41), all coefficients are positive during the iteration process. When the third-term on the right-side of equation (41) becomes zero, it appears as a positive source on the third-term of equation (41). The same is true for the fourth-term. This positive source will not lead to undershoots. It must be stressed that once a converged solution is obtained, (1) for situations where a_E > 0, equation (41) is identical to equation (9), and (2) for situations where a_E < 0, equation (41) sets the nodal value of the dependent variable ξ′P to the limiting value of the dependent variable given by equation (42); ensuring there is no undershoots. We also ensure that negative b does not lead to undershoots.

Overshoots: To eliminate overshoots, we will define:

(46) ξ″ ≡ ϕlimit,P−ϕ=ϕW+SΔxF−ϕ

Using equation (46), and after some algebraic rearrangement, the discretization equation for all Peclet numbers can be written as:

(47) (|aE|+aW)ξP″=|aE|ξE″+aWξW″−2min(aE,0)(ϕE−ϕP)−b

To avoid overshoots, we will write equation (47) as:

(48)  [|aE|+aW+(max[2min(aE,0)ϕE,0]+max[−2min(aE,0)ϕP,0]+max[b,0])ξ″P*]ξP″= |aE|ξE″+aWξW″−min[2min(aE,0)(ϕE−ϕP),0]−min[b,0]

To facilitate combined formulation, we will write equation (48) as:

(49) [|aE|+aW−(min[−2min(aE,0)ϕE,0]+min[2min(aE,0)ϕP,0]+min[−b,0])ξ″P*]ξP″=|aE|ξE″+aWξW″+max[−2min(aE,0)ϕE,0]+max[2min(aE,0)ϕP,0]+max[−b,0]

Now that we have the discretization equations to prevent undershoots and overshoots, and we will formulate a discretization equation applicable to both situations.

Combined Formulation: For this combined formulation, we will define a combined variable ξ:

(50) ξ ≡ ku ξ′+(1−ku) ξ″

where k_u is an undershoot indicator given by:

(51) ku ≡ max[sign(ϕE−ϕlimit),0]

Using equation (45), equation (49) and equation (50), a discretization equation for the combined formulation can be written as:

(52)  [|aE|+aW−(min[ks2min(aE,0)ϕE,0]+min[−ks2min(aE,0)ϕP,0]+min[ksb,0])ξP*]ξP=|aE|ξE+aWξW+max[ks2min(aE,0)ϕE,0]+max[−ks2min(aE,0)ϕP,0]+max(ksb,0)

where

(53) ks=sign(ϕE−ϕlimit)

Equation (52) ensures that undershoots, and overshoots are eliminated during the iteration process. The nodal ϕ_P are limited according to equation (50).

MUST with Pe < 2: In the absence of source terms, and when the cell Peclet number is less than 2, the coefficient a_E in equation (11a) is > 0. As a result, min(a_E, 0) in equation (52) becomes zero. Equation (52) can be simplified to:

(54) (|aE|+aW)ϕE=|aE|ϕE+aWϕW

It can be seen that the discretization equation with MUST reduces to the discretization equation of the original CD scheme.

MUST with Pe > 2: For this discussion, we will neglect the source term, and let ϕ_E > ϕ_W. For this situation, ϕ_limit = ϕ_W, k_s = 1 and k_u = 1. As the cell Peclet number is greater than 2, the coefficient a_E in equation (11a) is < 0. As a result, min(a_E, 0) in equation (52) is a_E. The terms min[k_s2min(a_E,0)ϕ_E, 0] = 2a_Eϕ_E, min[–k_s2min(a_E, 0)ϕ_P, 0] = 0, max[k_s2min(a_E, 0)ϕ_E,0] = 0 and max[–k_s2min(a_E, 0)ϕ_P, 0] = –2a_Eϕ_P. Equation (48) can be simplified as:

(55) (|aE|+aW−2aEϕEξP*)ξP=|aE|ξE+aWξW−2aEϕP

Upon convergence, ξP*=ξP, equation (55) becomes:

(56) (|aE|+aW)ϕP−2aEϕE=|aE|ϕE+aWϕW−2aEϕP

Since a_E < 0, equation (56) can be written as:

(57) (aE+aW)ϕP=aEϕE+aWϕW

The CD scheme with MUST reverts to the CD scheme without MUST [equations (5) and (7)]. This is all good. However, since a_E < 0, equation (57) will produce undershoots and overshoots just like the CD scheme without MUST.

This is avoided as MUST embeds Patankar’s internal value specification treatment described in subsection 2.5 to set the nodal value to the value of the limiter. During the iteration process when the cell Peclet number Pe is greater than 2, the nodal value ϕ_P approaches the limiter ϕ_limit. As a result, ξP*=ϕP*−ϕlimit approaches zero. The last term of the coefficient for ξ_P in equation (55):

(58) 2aEϕEξP*→−∞

Equation (55) becomes:

(59) (|aE|+aW+∞)ξP=|aE|ξE+aWξW−2aEϕP

which reduces to

(60) (∞)ξP=finite

(61) ξP→0

(62) ξP ≡ ϕP−ϕlimit→0

(63) ϕP→ϕlimit

It can be seen that MUST sets the nodal value of the dependent variable ϕ_P to the value of the limiter. The manifestation of this approach will be discussed when the solutions using the CD scheme are presented later in this article.

Remarks: Although the CD scheme is used as demonstration, the above treatment can be carried out for any spatial differencing scheme where undershoots and overshoots are the results of negative coefficient(s).

We presented an approach to eliminate undershoots, and overshoots as a result of physical sources in subsection 2.4. In subsection 3.1, we extended the approach to eliminate undershoots and overshoots due to negative coefficient(s). We will now apply MUST to the SOUD scheme.

4.2 The second-order upwind difference scheme

One possible implementation of the SOUD scheme for the convection–diffusion terms when u > 0 is:

(64) ρu(ϕP−ϕW)+ρu2(ϕP−2ϕW+ϕWW)=Γ[(ϕE−ϕPΔx)−(ϕP−ϕWΔx)]+SΔx

where WW is the west neighbor of control volume W as shown in Figure 2. Equation (64) can be written as:

(65) aPϕP=aEϕE+aWϕW+b

where

(66a) aE=1

(66b) aW=1+P

(66c) aP=1+(1+P)=aE+aW

(66d) b=P2(2ϕW−ϕP−ϕWW)+S(Δx)2Γ

This implementation of the SOUD scheme is chosen as it is a case where the undershoots and overshoots are the results of the sources. Equation (64) can be expanded as:

(67) (aE+aW)ϕP=aEϕE+aWϕW+P2(2ϕW−ϕP−ϕWW)+S(Δx)2Γ

Unlike the CD scheme, since the coefficients a_E and a_W are positive, equation (67) is applicable to all Peclet numbers as far as the coefficients are concern. We seek to eliminate undershoots and overshoots as the result of the source terms.

Undershoots: To eliminate undershoots, and using equation (42), equation (67) with MUST can be written as:

(68) [aE+aW−(P[min(2ϕW,0)+min(−ϕP,0)+min(−ϕWW,0)]/2+min[S(Δx)2/Γ,0])ξ′P*]ξ′P=aEξ′E+aWξ′W+(P[max(2ϕW,0)+max(−ϕP,0)+max(−ϕWW,0)]/2+max[S(Δx)2/Γ,0])

Similar to the CD scheme, the discretization equation to eliminate overshoots can now be formulated.

Overshoots: Using equation (46) to eliminate overshoots, equation (67) can be written as:

(69) [aE+aW−(P[min(−2ϕW,0)+min(ϕP,0)+min(ϕWW,0)]/2+min[S(Δx)2/Γ,0])ξ″P*]ξP″=aEξE″+aWξE″+(P[max(−2ϕW,0)+max(ϕP,0)+max(ϕWW,0)]/2+max[S(Δx)2/Γ,0])

Combined Formulation: A combined formulation that eliminates undershoots and overshoots can be written by using equation (50) and equation (53) to combine equation (68) and equation (69) as:

(70) [aE+aW−(P[min(2ksϕW,0)+min(−ksϕP,0)+min(−ksϕWW,0)]/2+min[S(Δx)2/Γ,0])ξP*]ξP=aEξE+aWξW+(P[max(2ksϕW,0)+max(−ksϕP,0)+max(−ksϕWW,0)]/2+max[S(Δx)2/Γ,0])

Remarks: We presented an example where MUST is used to eliminate undershoots and overshoots as the results of the source terms. We will next demonstrate MUST on a scheme where physically unrealistic solutions can be due to negative coefficients and the source terms.

4.3 The QUICK scheme

One possible implementation of the QUICK scheme when u > 0 is:

(71a) ϕe=68ϕP+38ϕE−18ϕW

(71b) ϕw=68ϕW+38ϕP−18ϕWW

The discretization equation can then be written as:

(72) ρu(68ϕP+38ϕE−18ϕW)−ρu(68ϕW+38ϕP−18ϕWW)=Γ[(ϕE−ϕPΔx)−(ϕP−ϕWΔx)]+SΔx

One possible way to write the discretization equation is:

(73) (2ΓΔx+6F8−3F8)ϕP=(ΓΔx−3F8)ϕE+(ΓΔx+6F8)ϕW+F8ϕW−F8ϕWW+SΔx

Or

(74) aPϕP=aEϕE+aWϕW+b

where

(75a) aE=1−3P8

(75b) aW=1+6P8

(75c) aP=aE+aW

(75d) b=P8ϕW−P8ϕWW+S(Δx)2Γ

In this situation, we demonstrate MUST on a scheme where physically unrealistic undershoots and overshoots can be the results of negative coefficients and/or source terms. Following our approach for the CD scheme, we will add different forms of zeroes to both sides of equation (74) as:

(76) [aE+aW+(|aE|+aE)]ϕP=aEϕE+aWϕW+(|aE|+aE)ϕE+P8ϕW−P8ϕWW+S(Δx)2Γ

Equation (76) can be rewritten as:

(77) (|aE|+aW)ϕP=|aE|ϕE+aWϕW+2min(aE,0)(ϕE−ϕP)+P(ϕW−ϕWW)8+S(Δx)2Γ

Undershoots: Similar to the other schemes, using equation (42), equation (77) with MUST can be written as:

(78) [aE+aW−(min[2min(aE,0)ϕE,0]+min[−2min(aE,0)ϕP,0])ξ′P*−(min[PϕW/8,0]+min[−PϕWW/8,0]+min[S(Δx)2/Γ,0])ξ′P*]ξ′P=aEξ′E+aWξ′W+max[2min(aE,0)ϕE,0]+max[−2min(aE,0)ϕP,0]+max(PϕW/8,0)+max(−PϕWW/8,0)+max[S(Δx)2/Γ,0]

We will take a pause here to examine the resulting discretization equations for the CD scheme and the QUICK scheme. Physically unrealistic solutions in the CD scheme are the result of negative coefficient. The QUICK scheme suffers from undershoots and overshoots because of negative coefficient and negative sources.

Comparing equation (78) and equation (45), it can be seen that the additional terms due to negative coefficient is the same for the two spatial differencing schemes. The exact forms of a_E are different.

Overshoots: Using equation (46), equation (77) can be written as:

(79) [aE+aW−(min[−2min(aE,0)ϕE,0]+min[2min(aE,0)ϕP,0])ξ″P*−(min[−PϕW/8,0]+min[PϕWW/8,0]+min[−S(Δx)2/Γ,0])ξ″P*]ξP″=aEξE″+aWξW″+max[−2min(aE,0)ϕE,0]+max[2min(aE,0)ϕP,0]+max[−PϕW/8,0]+max[PϕWW/8,0]+max[−S(Δx)2/Γ,0]

Combined Formulation: A combined formulation for equation (78) and equation (79) can be written using equation (50) and equation (53) as:

(80) [aE+aW−(min[2ksmin(aE,0)ϕE,0]+min[−2ksmin(aE,0)ϕP,0])ξ′P*−(min[PksϕW/8,0]+min[−PksϕWW/8,0]+min[Sks(Δx)2/Γ,0])ξ′P*]ξP′=aEξE′+aWξW′+max[2ksmin(aE,0)ϕE,0]+max[−2ksmin(aE,0)ϕP,0]+max(PksϕW/8,0)+max(−PksϕWW/8,0)+max[Sks(Δx)2/Γ,0]

Remarks: We showed how MUST is used to eliminate undershoots and overshoots due to negative coefficient and source in the QUICK scheme.

4.4 Monotonic unbounded schemes transformer procedure

The MUST procedure for u > 0, can be summarized as:

• Starting from the usual discretization equation [equations 7, 66 and 75), check the sign of a_E.

• If a_E can be negative, add different zeroes to both sides of the discretization equation as in equations (39) and (76).

• Subtracting the same terms from both sides as in equations (43), (47), (68), (69), (78) and (79).

• Rearranging the source terms to eliminate undershoots [equations 45, 68 and 78] and overshoots [equations 49, 69 and 79].

• Combine into combined formulation for undershoots and overshoots [equations 52, 70 and 80].

The same procedure can be used to derive the discretization equation for u < 0. The discretization equations for u > 0 and u < 0 can then be combined to write a discretization equation for all velocities.

4.5 Remarks

We incorporated MUST into three higher-order unbounded schemes to eliminate undershoots and overshoots. We will test the implementations using different test problems in the next section.

5. Sample applications of monotonic unbounded schemes transformer

In this section, we will demonstrate MUST’s capabilities. For simplicity of writing, Peclet number is used to refer to the cell Peclet number. For ease of analyses, and without loss of generality, we set P > 0. Therefore, in all discussions, the upstream neighbors are ϕ_WW and ϕ_W. The downstream neighbor is ϕ_E. Extension of MUST to when P < 0 is straightforward and does not require any new concept, so it is left to the explorations of interested readers. We divide our discussions into two parts, namely, without sources S = 0 and with sources S ≠ 0.

5.1 Zero source

The limit value: When there are no sources, as the Peclet number increases, the nodal ϕ_P value approaches the value of the upstream neighbor or ϕ_W in this case.

The CD scheme: Figure 5 shows the nodal values of the dependent variable ϕ_P, predicted using the exponential scheme [equation (10)], the CD scheme without MUST [equation (7)] and the CD scheme with MUST [equation (52)]. On the left-side of Figure 5, the upstream neighbor ϕ_W is set to 1 and downstream neighbor ϕ_E is set to 2[3]. Using the exponential scheme, the nodal ϕ_P value approaches the upstream neighbor ϕ_W exponentially as the Peclet number increases. The CD scheme without MUST does the same and the nodal value ϕ_P reaches the upstream neighbor’s value when the critical Peclet number of 2 is reached. After which, the predicted nodal value falls below the upstream neighbor’s value, which is physically not possible for this pure convection-diffusion problem. When MUST is applied to the CD scheme and when the Peclet number is larger than 2, MUST automatically limits the nodal value to the upstream neighbor’s value without any ad hoc interventions. The right-side of Figure 5 shows the nodal ϕ_P value when ϕ_W = 2 and ϕ_E = 1. Overshoots are encountered with the CD scheme without MUST when the Peclet number is larger than 2. MUST limits the nodal ϕ_P values to the upstream neighbor’s values. When the Peclet number is less than 2, the solutions from CD scheme with MUST and without MUST are identical. Table 2 shows the same information on the left-side of Figure 5 in tabular format. It can be seen that MUST sets the nodal value to the upstream neighbor value when the cell Peclet number is greater than 2. Table 3 shows the convergence history with different ϕ_P initial guesses for ϕ_W = 1 and ϕ_E = 2 with the cell Peclet number of 3. Although not shown, when the initial guess is set to ϕ_W which is the correct answer, the solution converges in one iteration. Three different initial guesses, namely, ϕP*=−ϕE, ϕ_E and 1,000ϕ_E are used to study the convergence history. We examine the last term of the coefficient of ξ_P in equation (52) defined as below:

(81) α ≡ −(min[ks2min(aE,0)ϕE,0]+min[−ks2min(aE,0)ϕP,0])ξP*

We do not normally drive the solution of the linear equation until converged solution is obtained during the iteration process. However, this exercise shows the robustness of MUST. From Table 3, it can be seen that the solutions converge well for all three choices of initial guesses. Even when a really bad and unlikely initial guess of 1,000ϕ_E, the solution converges really well. As α increases, ξP*→0. As a result, ϕ_P → ϕ_limit. Since the changes between iterations are not expected to be large, the three initial guesses can be considered unreasonable guesses. Even then, the solutions converge well.

Remarks: All the above are done automatically without any ad hoc interventions, and naturally by retaining all terms in the discretization equation of the CD scheme. For the CD scheme where the critical Peclet number is a constant known a-priori, MUST functions like the hybrid scheme of Spalding (Spalding, 1972). We will demonstrate MUST’s utilities and capabilities when the critical Peclet numbers are not constant and unknown a-priori.

SOUD scheme: Equation (66) is used to calculate the nodal ϕ_P values without MUST. The nodal ϕ_P values for SOUD scheme with MUST are calculated using equation (70). Unlike the CD scheme where the critical Peclet number is known and equals to 2, the critical Peclet number for the SOUD scheme is not constant and not known a-priori.

Figure 6(a) shows how MUST predicts the critical Peclet numbers and eliminates undershoots. In this undershoots study, the upstream neighbor and downstream neighbor are set to ϕ_W = 1 and ϕ_E = 2, respectively. The upstream neighbor of ϕ_W, namely, ϕ_WW are set to 3, 2, 1.5 and 1. The ϕ_WW values are chosen so that the critical Peclet numbers occur around 1, 2, 4 and never. From the figure, we can see that MUST sets the nodal ϕ_P value to the value of the upstream neighbor at the Peclet numbers where undershoots begin with the SOUD scheme without MUST. Table 2 shows these in tabulated format. It can be seen that depending on the values of ϕ_WW, undershoots begin at different Peclet numbers. Without MUST, the nodal values continue to decrease leading to larger undershoots (red). With MUST, the nodal value is limited to the upstream value of 1 (blue). From Table 2 and Figure 6(a), for Peclet numbers where there are no undershoots, as expected, the predictions for SOUD scheme with or without MUST are identical.

Figure 6(b) shows how MUST removes overshoots. For this demonstration, the upstream neighbor is set to ϕ_W = 2, and the downstream neighbor is set to ϕ_E = 1. The values of ϕ_WW are set to 0, 1, 1.5 and 2. Again, these ϕ_WW values are chosen so that overshoots occur at around Peclet numbers of 1, 2, 4 and never. Note that these ϕ_WW are different from the undershoots study. This shows that MUST can capture the critical Peclet numbers for different combinations of ϕ_E, ϕ_W and ϕ_WW. Similar to undershoots, MUST limits the nodal values to the upstream neighbor value when SOUD scheme without MUST starts to predict values higher than the upstream neighbor. Again, when there are no overshoots, the solutions by SOUD scheme with MUST collapse to the solutions of SOUD scheme without MUST.

The QUICK scheme: Equation (75) is used to calculate the nodal ϕ_P values for QUICK without MUST. The nodal ϕ_P values for QUICK with MUST are calculated using equation (80). Similar to the SOUD scheme, the critical Peclet numbers at which undershoots and overshoots begin are not known a-priori and change depending on the combinations of the values of the dependent variables at the neighboring nodes.

Figure 7(a) shows how MUST eliminates undershoots. In this undershoots study, the upstream neighbor is set to ϕ_W = 1, and the downstream neighbor is set to ϕ_E = 2. The values of ϕ_WW are set to 3, 2, 1.5 and 1. The same ϕ_WW values as those used for the CD scheme are chosen to demonstrate the ability of MUST to capture undershoots at different Peclet numbers with the same three ϕ_W, ϕ_W and ϕ_WW nodal values. Unlike the SOUD scheme, undershoots are observed for all ϕ_WW values. For these values of ϕ_WW, undershoots occur in between 1 < P < 3. We can see again, that MUST set the nodal value ϕ_P to the value of the upstream neighbor at the Peclet numbers where undershoots begin with the QUICK without MUST.

Figure 7(b) shows how MUST removes overshoots. For this demonstration, the upstream neighbor is set to ϕ_W = 2, and the downstream neighbor is set to ϕ_E = 1. The values of ϕ_WW are set to 0, 1, 1.5 and 2. Similar to undershoots, MUST limits the nodal value to the upstream neighbor value when QUICK without MUST starts to predict values higher than the upstream neighbor. Again, for Peclet number before overshoots, the solutions by QUICK with MUST collapse to the solutions of QUICK without MUST.

Remarks: We tested MUST on three different spatial difference schemes. MUST was able to predict the critical Peclet numbers where undershoots or overshoots begin and limit the nodal value of the dependent variable to the value of its upstream neighbor. The critical Peclet number for the CD scheme is known and fixed at 2. The critical Peclet numbers for the SOUD scheme and the QUICK scheme are functions of the relative values of the neighbors. MUST predicts these critical Peclet numbers automatically and naturally.

Another point to be noted is the undershoots and overshoots in the CD scheme are due to the negative neighbor coefficient. In its current implementation[4], the overshoots and undershoots in the SOUD scheme are due to the source term (equation 66d). The overshoots and undershoots in the QUICK scheme are due to the combined effects of the negative neighbor coefficient (equation 75a), and the source term (equation 75d). MUST removes the undershoots and overshoots as a result of these two factors very well.

In this section, we demonstrated MUST’s ability (1) to predict the critical Peclet numbers based on the relative values of the neighbors, and (2) to limit the nodal values to the upstream neighbor’s values. Both undershoots and overshoots were avoided. When the Peclet number is smaller than the critical Peclet number, MUST’s solutions are identical to the original discretization equation without MUST.

5.2 Finite source

When the source is zero, the Peclet number is the independent variable. However, when the source is finite, the Peclet number alone is not sufficient. For ease of discussions, and without loss of generality, for the remainder of this subsection, we will set Δx = 1 and Γ = 1. With these, the Peclet number becomes the independent variable. In this subsection, we will study four cases, namely, (1) ϕ_WW = 3, ϕ_W = 1, ϕ_E = 2 and S = 1, (2) ϕ_WW = 3, ϕ_W = 1, ϕ_E = 2 and S = –1, (3) ϕ_WW = 1, ϕ_W = 2, ϕ_E = 1 and S = 1 and (4) ϕ_WW = 1, ϕ_W = 2, ϕ_E = 1 and S = –1.

The limit values and the exponential scheme: When the source term is finite, and constant in a convection–diffusion problem, there is no exact solution even for one-dimensional problems. We will examine the solutions obtained using the exponential scheme [equation (10)] and the limit values [equation (37)]. Figure 8 shows the values of ϕ_P and the limit values for the four cases. For Case 1, when the Peclet number is zero, and as a result of the positive source, the nodal value of 2 is larger than the average of ϕ_W and ϕ_E, which is 1.5. As the Peclet number increases, the nodal value ϕ_P approaches the upstream value of ϕ_W. Note that as the source is positive, the nodal value ϕ_P is approaching but always higher than ϕ_W. The limit value which neglects diffusion, on the other hand approaches infinity when the Peclet number is zero. It decreases inversely with the Peclet number. Both the exponential scheme value and the limit value are identical at higher Peclet numbers where diffusion does not play an important role. In Case 2, the same thing is observed when the source is negative. The two (limit and exponential scheme) values are identical at large Peclet numbers. As a result of the negative source, the nodal value ϕ_P is smaller than the upstream value ϕ_W. As the Peclet number increases, the nodal value recovers and approaches ϕ_W Unlike the exponential scheme, as diffusion is neglected, the limit values approach infinity as the Peclet number approaches zero. Cases 3 and 4 are similar to Cases 1 and 2.

Case 1 (ϕ_WW = 3, ϕ_W = 1, ϕ_E = 2 and S = 1): We will examine the relation between the nodal value ϕ_P and the limit value ϕ_limit. Recalling that Δx = 1 and Γ = 1, the limit value (equation (37)) becomes:

(82) ϕlimit,P=ϕW+SP=1+1P

With ϕ_E = 2, the following can be concluded:

(83) P<1ϕlimit,P>2ϕlimit,P>ϕEP>1ϕlimit,P<2ϕlimit,P<ϕE

Since ϕ_W < ϕ_E, we must prevent:

(84) P<1overshootϕP≤ϕlimit,P=1+1PP>1undershootϕP≥ϕlimit,P=1+1P

The above applies to all differencing schemes.

Figure 9 shows the predictions of the CD scheme, the SOUD scheme and the QUICK scheme with and without MUST. Figure 9(a) shows the CD scheme predictions match the exponential scheme well up to P = 2. For this case, when P > 2, the CD scheme admits undershoot values. The predictions with MUST are identical to those without MUST when P ≤ 2. When P > 2, MUST becomes effective. When the nodal value ϕ_P is less than the limit value, equation (84) limits the nodal value to the limit values. Figure 9(a) shows the CD scheme’s predictions with MUST match those of the exponential scheme very well.

Figure 9(a) shows the SOUD scheme without MUST predicts the nodal ϕ_P values lower than the limit values for all Peclet numbers. As undershoots are encountered when P < 1 and according to equation (84), we want to prevent overshoots. Therefore, nodal values lower than the limit values were admitted as solutions. As a result, the SOUD scheme with MUST follows the SOUD scheme without MUST until P = 1. After this point, MUST ensures that there is no undershoot and predicts nodal ϕ_P values according to equation (84).

Similar to the SOUD scheme, the QUICK scheme without MUST predicts nodal ϕ_P values lower than the limit values for all Peclet numbers. Once MUST is activated, the nodal ϕ_P values are set to the limit values.

Case 2 (ϕ_WW = 3, ϕ_W = 1, ϕ_E = 2 and S = –1): We will examine the relation between the nodal value ϕ_P, and the limit value ϕ_limit. Recalling that Δx = 1 and Γ = 1, the limit value (equation (37)) becomes:

(85) ϕlimit,P=ϕW+SP=1−1P

With ϕ_E = 2, the following can be concluded:

(86) P<1ϕlimit,P<0ϕlimit,P<ϕEP>10<ϕlimit,P<1ϕlimit,P<ϕE

Since ϕ_W < ϕ_E, we must prevent:

(87) P<1undershootϕP≥ϕlimit,P=1−1PP>1undershootϕP≥ϕlimit,P=1−1P

The above applies to all differencing schemes.

Figure 9(b) shows the results for the three schemes tested in this article. At small Peclet numbers, all three schemes predicted values lower than the exponential scheme but higher than the limit value. Once the values of the different schemes intersect the limit values, MUST ensures that the solutions of the various schemes are limited to the limit values; removing undershoots in the process. As Peclet number increases, all solutions (including the exponential scheme) match the limit values. This is because the contribution from the diffusion term diminishes with increasing Peclet number.

Case 3 (ϕ_WW = 1, ϕ_W = 2, ϕ_E = 1 and S = 1): With Δx = 1 and Γ = 1, the limit value [equation (37)] becomes:

(88) ϕlimit,P=ϕW+SP=2+1P

With ϕ_E = 1, the following can be concluded:

(89) P<1ϕlimit,P>2ϕlimit,P>ϕEP>1ϕlimit,P>2ϕlimit,P>ϕE

Since ϕ_W > ϕ_E, we must prevent:

(90) P<1overshootϕP≤ϕlimit,P=2+1PP>1overshootϕP≤ϕlimit,P=2+1P

The above applies to all differencing schemes.

From Figure 10(a), it is clear that at small Peclet numbers, the CD scheme and the QUICK scheme produce solutions that are larger than the exponential scheme. At small Peclet numbers, all three schemes produce solutions which are smaller than the limit values. Therefore, some solutions larger than the exponential scheme (but smaller than the limit values) are admitted. Once the solutions (by the different schemes) intersect the limit values, MUST limits the solutions to the limit values eliminates uncontrolled overshoots.

Case 4 (ϕ_WW = 1, ϕ_W = 2, ϕ_E = 1 and S = –1): With Δx = 1 and Γ = 1, the limit value [equation (37)] becomes:

(91) ϕlimit,P=ϕW+SP=2−1P

With ϕ_E = 1, the following can be concluded:

(92) P<1ϕlimit,P<1ϕlimit,P<ϕEP>1ϕlimit,P>1ϕlimit,P>ϕE

Since ϕ_W > ϕ_E, we must prevent:

(93) P<1undershootϕP≤ϕlimit,P=2−1PP>1overshootϕP≤ϕlimit,P=2−1P

Similar to the other three cases, it can be seen from Figure 10(b) that at small Peclet numbers the solutions with MUST are identical to those without MUST. As Peclet numbers increase, MUST prevents overshoots and limits the nodal values to the limit values.

6. Concluding remarks

We presented an approach to monotonize unbounded higher-order/upwind schemes which we call MUST. We demonstrated MUST using the CD scheme, the SOUD scheme and the QUICK scheme.

For problems with zero sources, the limit value of the nodal point is the value of the upstream neighbor. In our tests, MUST is able to (1) capture the critical Peclet numbers, (2) eliminates undershoots and overshoots by setting the nodal value to the value of its upstream neighbor. For Peclet numbers below the critical Peclet numbers, the solutions obtained using MUST are identical to the solutions without MUST. This is accomplished for known constant critical Peclet number (CD scheme) and unknown Peclet numbers which depends on the relative values of ϕ_WW, ϕ_W and ϕ_E (SOUD scheme and QUICK scheme).

For problem with finite sources, a limit value is proposed to demonstrate the capability of MUST. Other limit values can of course be formulated. Again, below the critical Peclet numbers, solutions with MUST are identical to those without MUST. Beyond the critical Peclet numbers, MUST sets the nodal values to the limit values and eliminates uncontrolled overshoots and undershoots.

The above is achieved automatically and without ad hoc intervention. We also accomplished this naturally by keeping all the terms in the discretization equation of the unbounded schemes.

Figures

Figure 1.

Four spatial difference schemes (a) central; (b) upwind; (c) second-order upwind and (d) QUICK

Figure 2.

One-dimensional control volume

Figure 3.

Steady convection–diffusion with zero source – undershoots

Figure 4.

Steady convection–diffusion with zero source – overshoots

Figure 5.

Values of ϕ_P using CD scheme with and without MUST

Figure 6.

Values of ϕ_P using second-order upwind scheme with and without MUST

Figure 7.

Values of ϕ_P using QUICK scheme with and without MUST

Figure 8.

Values of ϕ_P with exponential scheme and the limiter

Figure 9.

Values of ϕ_P with ϕ_WW = 3, ϕ_W = 1 and ϕ_E = 2

Figure 10.

Values of ϕ_P with ϕ_WW = 1, ϕ_W = 2 and ϕ_E = 1

Figure A1.

Two-dimensional control volume

Table 1.

Predicted ϕ_P using the central difference scheme

ϕW	ϕE	ϕP	Remarks
0.0	1.0	−0.5	Undershoot
0.5	1.0	0.375
1.0	0.0	1.25	Overshoot
1.0	0.5	1.125

Source: Table by authors

Table 2.

Nodal ϕ_P values for the SOUD scheme with ϕ_W = 1, ϕ_E = 2 and different ϕ_WW values

Peclet	Expo	ϕWW = 3		ϕWW = 2		ϕWW = 1.5		ϕWW = 1
		Normal	MUST	Normal	MUST	Normal	MUST	Normal	MUST
0.10	1.48	1.42	1.42	1.44	1.44	1.45	1.45	1.47	1.47
0.20	1.45	1.35	1.35	1.39	1.39	1.41	1.41	1.44	1.44
0.30	1.43	1.29	1.29	1.35	1.35	1.38	1.38	1.41	1.41
0.40	1.40	1.23	1.23	1.31	1.31	1.35	1.35	1.39	1.39
0.50	1.38	1.18	1.18	1.27	1.27	1.32	1.32	1.36	1.36
0.80	1.31	1.06	1.06	1.19	1.19	1.25	1.25	1.31	1.31
0.90	1.29	1.03	1.03	1.16	1.16	1.23	1.23	1.30	1.30
1.00	1.27	1.00	1.00	1.14	1.14	1.21	1.21	1.29	1.29
1.10	1.25	0.97	1.00	1.12	1.12	1.20	1.20	1.27	1.27
1.20	1.23	0.95	1.00	1.11	1.11	1.18	1.18	1.26	1.26
1.80	1.14	0.83	1.00	1.02	1.02	1.12	1.12	1.21	1.21
1.90	1.13	0.81	1.00	1.01	1.01	1.11	1.11	1.21	1.21
2.00	1.12	0.80	1.00	1.00	1.00	1.10	1.10	1.20	1.20
2.10	1.11	0.79	1.00	0.99	1.00	1.09	1.09	1.19	1.19
2.20	1.10	0.77	1.00	0.98	1.00	1.09	1.09	1.19	1.19
4.00	1.02	0.63	1.00	0.88	1.00	1.00	1.00	1.13	1.13
4.10	1.02	0.62	1.00	0.87	1.00	1.00	1.00	1.12	1.12
4.20	1.02	0.61	1.00	0.87	1.00	0.99	1.00	1.12	1.12
4.30	1.01	0.61	1.00	0.86	1.00	0.99	1.00	1.12	1.12
4.40	1.01	0.60	1.00	0.86	1.00	0.99	1.00	1.12	1.12
7.40	1.00	0.51	1.00	0.79	1.00	0.94	1.00	1.08	1.08
7.50	1.00	0.51	1.00	0.79	1.00	0.93	1.00	1.08	1.08
7.60	1.00	0.51	1.00	0.79	1.00	0.93	1.00	1.08	1.08
7.70	1.00	0.51	1.00	0.79	1.00	0.93	1.00	1.07	1.07
7.80	1.00	0.50	1.00	0.79	1.00	0.93	1.00	1.07	1.07

Source: Table by authors

Table 3.

Convergence history for the CD scheme for cell Peclet number of 3 with different ϕ_P initial guesses for a problem with zero source, ϕ_W = 1 and ϕ_E = 2

Initial guess	ϕP*=−ϕE			ϕP*=ϕE			ϕP*=1000ϕE
Iteration	α (Eq. 81)	ξP* (Eq. 52)	ϕP* (Eq. 50)	α (Eq. 81)	ξP* (Eq. 52)	ϕP* (Eq. 50)	α (Eq. 81)	ξP* (Eq. 52)	ϕP* (Eq. 50)
1	1.33E + 01	3.00E−01	1.300	4.00E + 00	5.00E−01	1.500	3.00E−03	6.67E + 02	667.611
5	3.38E + 01	5.91E−02	1.059	2.46E + 01	8.12E−02	1.081	2.32E−01	8.63E + 00	9.633
10	1.63E + 02	1.22E−02	1.012	1.24E + 02	1.61E−02	1.016	9.78E + 00	2.04E−01	1.204
15	7.10E + 02	2.82E−03	1.003	5.46E + 02	3.67E−03	1.004	6.10E + 01	3.28E−02	1.033
20	3.01E + 03	6.64E−04	1.001	2.32E + 03	8.62E−04	1.001	2.78E + 02	7.19E−03	1.007
25	1.27E + 04	1.57E−04	1.000	9.80E + 03	2.04E−04	1.000	1.19E + 03	1.68E−03	1.002
30	5.36E + 04	3.73E−05	1.000	4.13E + 04	4.84E−05	1.000	5.05E + 03	3.96E−04	1.000
35	2.26E + 05	8.85E−06	1.000	1.74E + 05	1.15E−05	1.000	2.13E + 04	9.39E−05	1.000

Source: Table by authors

Table A1.

Nodal ϕ_P values using different schemes with and without MUST with ϕ_W = 1, ϕ_E = 2, ϕ_S = 1.5, ϕ_N = 2, ϕ_WW = 1, ϕ_SS = 2.5, and the same Peclet number in both coordinate directions

Peclet	Expo	Upwind	Hybrid	Central		Second-order upwind		QUICK
				Normal	MUST	Normal	MUST	Normal	MUST
0.10	1.61	1.61	1.61	1.61	1.61	1.59	1.59	1.60	1.60
0.20	1.59	1.59	1.59	1.59	1.59	1.55	1.55	1.58	1.58
1.00	1.45	1.50	1.44	1.44	1.44	1.39	1.39	1.42	1.42
1.10	1.44	1.49	1.42	1.42	1.42	1.38	1.38	1.40	1.40
1.20	1.42	1.48	1.40	1.40	1.40	1.37	1.37	1.39	1.39
2.00	1.34	1.44	1.25	1.25	1.25	1.30	1.30	1.27	1.27
2.10	1.33	1.43	1.25	1.23	1.25	1.29	1.29	1.26	1.26
2.20	1.33	1.43	1.25	1.21	1.25	1.29	1.29	1.25	1.25
2.30	1.32	1.42	1.25	1.19	1.25	1.28	1.28	1.24	1.25
2.40	1.31	1.42	1.25	1.18	1.25	1.28	1.28	1.22	1.25
2.50	1.31	1.42	1.25	1.16	1.25	1.27	1.27	1.21	1.25
3.00	1.29	1.40	1.25	1.06	1.25	1.25	1.25	1.16	1.25
3.10	1.28	1.40	1.25	1.04	1.25	1.25	1.25	1.15	1.25
3.20	1.28	1.39	1.25	1.03	1.25	1.24	1.25	1.14	1.25
3.30	1.28	1.39	1.25	1.01	1.25	1.24	1.25	1.13	1.25
3.40	1.27	1.39	1.25	0.99	1.25	1.24	1.25	1.12	1.25
5.00	1.26	1.36	1.25	0.69	1.25	1.20	1.25	1.00	1.25
5.10	1.26	1.36	1.25	0.67	1.25	1.20	1.25	0.99	1.25
5.20	1.25	1.35	1.25	0.65	1.25	1.19	1.25	0.99	1.25
5.30	1.25	1.35	1.25	0.63	1.25	1.19	1.25	0.98	1.25
5.40	1.25	1.35	1.25	0.61	1.25	1.19	1.25	0.98	1.25

Source: Table by authors

Notes

1. In this article, sink is understood to be negative source. Therefore, source implies positive and negative sources.

2. There is also question if the upstream value is the right limiting value for ϕ_P.

3. Neighbor values of 1 and 2 (instead of 0 and 1) are used to show MUST’s capabilities to set the nodal value to any generic value rather than zero.

4. Although not shown, but we tried other implementations [see (10, 11)], the same conclusion can be made.

Appendix

In this Appendix, we present results for a two-dimensional problem to show the generality of MUST to eliminate undershoots. The inclusion of the complete formulations for two-dimensional problems with proper explanations is not attempted in this article as it will dilute the focus of this article. Furthermore, even in today’s situations, one-dimensional formulations are used in many practical modeling. As a result, the complete two-dimensional MUST will be presented in a follow-up article.

There are no new concepts in the formulation of two-dimensional discretization equation using MUST. Therefore, it is not discussed here. The key is in the formulation of the limiter value. For a two-dimensional Cartesian-coordinates control volume shown in Figure A1 with u > 0, v > 0 and without source, the limiter value becomes:

(A1)  ϕlimit=(ρuΔy)wϕW+(ρvΔx)sϕS(ρuΔy)w+(ρvΔx)s

The limiter value in Eq. (A1) is the mixing-cup value.

Table A1 shows the results for the simplest steady, two-dimensional convection–diffusion problem without source. We examined the exponential scheme, the upwind scheme, the hybrid scheme, the CD scheme with and without MUST, a SOUD scheme with and without MUST and the QUICK scheme with and without MUST. For this test problem, we set u > 0, v > 0 and the cell Peclet number (Pe) in both coordinate directions to be equal. The values of the immediate upstream neighbors are set to ϕ_W = 1 and ϕ_S = 1.5. According to equation (1), the limiter value is 1.25. The downstream neighbors are set to ϕ_E = 2 and ϕ_N = 2. The two other upstream neighbors are set to ϕ_WW = 1 and ϕ_SS = 2.5. The exponential scheme predicts the correct large cell Peclet number limiter value when the cell Peclet number reaches 5.2. The upwind scheme approaches the limiter value slower than the exponential scheme. As expected, the hybrid scheme set diffusion to zero and predicts the limiter value when the cell Peclet number reaches 2. The CD scheme predicts values below (red) the limiter value of 1.25 when the cell Peclet number exceeds 2. The CD scheme with MUST eliminates the undershoot (ϕ_P < 1.25) starting from cell Peclet number of 2 (blue). Unlike the CD scheme where the starting point of the undershoot is known, undershoots (red) in the SOUD scheme begin when the cell Peclet number reaches 3.2. The onset of undershoots is not fixed and known a priori as in the CD scheme. The onset depends on the relative values of the four upstream neighbors, and the two downstream neighbors. SOUD scheme with MUST predicts the onset of undershoot and eliminates undershoots (blue). A similar trend is observed for the QUICK scheme. Here the undershoots begin at cell Peclet number of 2.3.

For this simple test problem, MUST (1) predicts the same solutions (as those without MUST) before the onset of undershoots, (2) determine the onset of undershoots automatically and naturally, (3) set the nodal value to the limiter value once undershoot starts and (4) predicts the correct large cell Peclet number limit.