1 Introduction

In computing the Navier-Stokes equation, the choices of discretization schemes and grid size are critical for obtaining reasonable solutions. In this paper, a new method of choosing discretization schemes and grid size is presented. This new method mainly addresses two aspects: the practical accuracy of the diffusion term discretization and the behavior of high wavenumber disturbances. An accurate diffusion approximation is essential in studying flows, such as boundary layer flows, turbulent flows, and separated flows. Attention should be paid to diffusion accuracy in order to calculate viscous effects accurately. One example is the convection-diffusion equation

(1)

where δ stands for numerical discretization operator.

The theoretical accuracy of diffusion is determined by . Assume that the leading term in the convection terms' truncation error is diffusive. Then, the numerical diffusion is . The total diffusive truncation error in Eq. (1) is . The practical diffusion accuracy is determined by this total diffusive truncation error, and it is not able to reach its theoretical accuracy if the numerical diffusion overlaps the error of the diffusion terms, that is, if . Therefore, the way in which convection discretization influences the practical diffusion accuracy should be carefully analyzed^[1-3]. The key lies in the combination of the convection scheme, diffusion scheme, and grid size. High wavenumber modes in numerical solutions commonly originate from numerical disturbances, such as non-uniformities of meshes, inaccuracies of boundary conditions, and nonlinear interactions during calculation^[4]. Their scales are much smaller than those we are interested in. Improper choice of schemes and grid size may amplify these high wavenumber modes, making them not only an interference at the length scales of interest^[5-6] but also a source of numerical instabilities^[7-8]. Consequently, attention should be concentrated on how the numerical method influences the high wavenumber modes. Theoretical analysis of the Navier-Stokes equation is very difficult due to its strong nonlinearity, but the basic convection and diffusion effects in the Navier-Stokes equation can be represented by the convection-diffusion equation. Therefore, the convection-diffusion equation is an ideal simplified model of the Navier-Stokes equation for theoretical analysis^[9-10].

Many researchers have studied the influence of numerical diffusion from different angles, and various methods for choosing suitable discretization schemes and grid size have been proposed. Zhang et al.^[11] compared the size of diffusion and numerical diffusion terms of the Navier-Stokes equation qualitatively and put forward a criterion for selecting the scheme and grid size: requiring diffusion terms to be larger than the numerical diffusion. Fu and Ma^[12] adopted the convection-diffusion equation as a model equation and studied the influence of convection discretization on diffusion approximation using Fourier analysis. Wen and Gao^[13] analyzed the truncation error of various convection-diffusion discretization schemes as applied to the Burgers' equation and discussed how numerical diffusion affects diffusion accuracy and numerical stability. Tu et al.^[14] adopted the Burgers' equation as a model equation and discovered numerical oscillation caused by inappropriate second derivative schemes, and a new staggered cell-edge high-order compact scheme was proposed to avoid this problem. Mohammadian^[8] compared the Fourier properties of two schemes for viscous terms in shallow flows and found that a proper treatment of the viscous terms could avoid the oscillations of the shortest waves. Li et al.^[5] proposed a scheme that could compensate for the dissipation deficiency of traditional linear schemes and suppress the spurious energy accumulation that occurs at high wavenumbers.

In this paper, the convection-diffusion equation is adopted as a model of the Navier-Stokes equation. The Fourier analysis^[15] is used in the discussion of both the practical diffusion accuracy and the high wavenumber modes' behavior. A new method for choosing proper schemes for convective flux, diffusive flux, and grid size is proposed. According to this new method, proper scheme combinations and grid sizes satisfy three conditions: (ⅰ) Maintain high resolution near the length scale of interest. (ⅱ) Prevent the diffusion accuracy from being polluted by numerical diffusion near the length scale of interest. (ⅲ) Dissipate high wavenumber numerical disturbances.

2 Fourier analysis of the one-dimensional convection-diffusion equation

For spatial discretization, the error of convective flux discretization includes both numerical diffusion and numerical dispersion, while the error of diffusive flux discretization is diffusive. The Fourier analysis is used to quantify and compare the errors of these different origins. The one-dimensional convection-diffusion equation is

(2)

where c is the field velocity, and ν is the diffusion coefficient. Assume that u is periodic over [0, 2π] and ∆x=2π/N. Decompose the initial value u(x, 0) into Fourier components,

(3)

The exact solution to the convection-diffusion equation is

(4)

Discretize the convection term and the diffusion term by linear difference schemes to obtain the semi-discrete form of Eq. (2),

(5)

The numerical solution for this semi-discrete form is

(6)

where α =k∆x is the scaled wavenumber, and k_e = k_r + ik_i is the modified wavenumber of the convection discretization scheme. In this expression, k_i represents its dispersion property, and k_r represents its dissipation property. k_d is the modified wavenumber of the diffusion discretization scheme. As ∆x → 0, these parameters satisfy

(7)

The numerical solution (6) can be regarded as the product of the numerical influence and the exact solution (see Eq. (4)). Decompose the numerical influence into three parts, i.e., the error of diffusion terms , numerical diffusion , and numerical dispersion . In this way, Eq. (6) can be transformed into Eq. (8) so that the numerical influence of different origins can be compared and analyzed individually,

(8)

3 An approach for selecting discretization schemes and grid size

Each position in a flow contains one characteristic length scale. This length scale reflects the character of the turbulent structure and hence needs to be computed accurately. The convection-diffusion equation is used as a simplified model for theoretical analysis of the Navier-Stokes equation. The grid size and schemes for convective and diffusive flux should satisfy three conditions simultaneously:

(ⅰ) Maintain high resolution near the length scale of interest.

(ⅱ) Prevent the accuracy of the diffusion term from being polluted by numerical diffusion near the length scale of interest.

(ⅲ) Dissipate high wavenumber numerical disturbances.

3.1 Analysis of scheme resolution

Assume that the characteristic length scale is l₀, with the corresponding wavenumber k₀ = 2π/l₀ and characteristic scaled wavenumber α₀ = k₀∆x. The goal of resolution analysis is to find a proper range for α₀. Within this range, the adopted schemes can maintain high resolution quality. Resolution properties are commonly studied by the Fourier analysis. For convection schemes, resolution properties are mainly described by dispersion characters, which are quantified by k_i. For diffusion schemes, resolution properties are described by the modified wavenumber k_d. Plots of k_i and k_d versus α for various schemes are shown in Figs. 1 and 2, respectively.

Define the relative resolution error

(9)

For a given tolerance (ε_{ki max} and ε_{kd max}), there is a maximum well-resolved scaled wavenumber (α_{ki max} and α_{kd max}) for each scheme. Several schemes' maximum well-resolved scaled wavenumbers and their corresponding tolerances are tabulated in Tables 1 and 2. To guarantee high resolution near α₀, the grid size adopted should keep the interested length scales within this range, that is,

(10)

The majority of schemes are capable of maintaining high resolution when α₀∈ [0, 0.8] and α₀∈ [0, 1.25]. Therefore, in the following discussion, [0, 0.8] and [0, 1.25] are adopted as the interested scaled wavenumber ranges.

3.2 Analysis of practical accuracy of diffusion terms

The first criterion concerns practical accuracy of the diffusion terms. As mentioned above, if numerical diffusion originating from the convection term overlaps the discretization error of the diffusion terms, the practical diffusion accuracy is affected by the numerical diffusion. Hence, the practical diffusion accuracy is not able to reach its theoretical accuracy. Under this condition, it is hard to determine how accurate the diffusion term truly is. Therefore, to obtain a reliable solution, the theoretical diffusion accuracy should be maintained within the wave range of interest. The essence is to ensure that numerical diffusion is smaller than the approximation error of diffusion terms.

The Fourier analysis presented in Section 2 shows that the diffusion approximation error and the numerical diffusion can be respectively expressed by and . Their sizes are compared using the ratio

(11)

where

is the cell Reynolds number. The requirement that numerical diffusion should be smaller than the diffusion term's approximation error in the wave range of interest is equivalent to the requirement that

(12)

This inequality is the statement of the first criterion.

Figures 3 and 4 are the plots of against α at Re_∆x=1 for various convection and diffusion scheme combinations. In Fig. 3, the diffusion term scheme is a 2nd-order central scheme, and the convection term scheme changes from a 2nd-order upwind scheme to a 7th-order upwind scheme. In Fig. 4, the diffusion term scheme is a 4th-order central scheme, and the convection schemes are the same as those in Fig. 3. These scheme combinations can be divided into two categories based on . As α → 0, the results of the Fourier analysis are the same as those of the leading truncation error analysis. Therefore, these scheme combinations can be divided into the same two categories based on the truncation error of the convection scheme and the truncation error of the diffusion scheme .

(ⅰ) , which is equivalent to c ≤ d. In this case, numerical diffusion is always larger than the diffusion approximation error. Equation (12) can never be satisfied for this kind of scheme combination regardless of how small ∆x is. The practical diffusion accuracy is always lower than its theoretical accuracy.

(ⅱ) , which is equivalent to c > d. In this case, numerical diffusion is of the same order or higher order than the diffusion approximation error. Equation (12) can be met by choosing ∆x properly, and the theoretical diffusion accuracy can be maintained.

Therefore, for a scheme combination at a fixed Re_∆x, it either will never satisfy Eq. (12) or there exists a largest acceptable ∆x. Similarly, for a scheme combination at a given wave range [0, α₀], it either never satisfies Eq. (12) or there exists a maximum acceptable Re_∆x. Tables 3 and 4 tabulate the largest Re_∆x for several scheme combinations when α₀=0.8 and α₀=1.25. According to the discussion above, to preserve the theoretical diffusion discretization accuracy at the wave range of interest, the convection scheme should be at least one order higher than the diffusion scheme. Adopting a higher-order or compact convection scheme can increase Re_{∆x max} significantly. Increasing Re_{∆x max} is equivalent to increasing the largest acceptable ∆x.

When Eq. (12) is employed as a guideline in the Navier-Stokes equation, the local velocity and viscosity should be used for computing Re_∆x. For most aerodynamic problems, viscosity is essential only in some parts of the flow field where the cell Reynolds number is much smaller than it is in a free stream.

3.3 Influence on high wavenumber numerical disturbances

The second criterion concerns the way in which discretization methods influence high wavenumber modes. High wavenumber modes are close to the smallest length scale that the grid can distinguish and are much higher than our wave range of interest. They are traveling at a wrong speed and direction because of numerical dispersion. Research shows that an improper choice of schemes may amplify these high wavenumber modes and cause non-physical solutions. Therefore, a proper numerical method should prevent high wavenumber modes from increasing.

According to Eq. (6), the amplitude of each mode is determined by . Decompose this term into two parts: the effects of the exact solution e^(−k2νt) and the effects of the schemes . High frequency modes are defined as modes close to the resolution limit of a grid, that is, k→ ± N/2. Then, the requirement that the numerical method dissipates high frequency modes can be expressed as

(13)

Equation (13) is equivalent to

(14)

This expression is determined by both the convection scheme and the diffusion scheme. At high wavenumbers, upwind convection schemes are highly dissipative, while central diffusion schemes suffer from dissipation deficiency, that is,

(15)

Employing Eq. (15) and transforming Eq. (14) into an expression of , we obtain our second criterion, that is,

(16)

This criterion does not conflict with the first criterion, as they are applied to different wavenumber ranges. The first criterion concerns α ∈ [0, α₀], while the second one is for α → π.

Examples of the effects of the schemes ||N_semi|| are plotted in Figs. 5 and 6, where c=1, ν =0.004, ∆x=0.01, and therefore, Re_∆x=2.5. When the diffusion terms are discretized by a 2nd-order central scheme, its combinations with 2nd-order upwind and 3rd-order compact schemes satisfy Eq. (16). When the diffusion terms are discretized by a 4th-order central scheme, its combinations with 2nd-order upwind, 3rd-order compact, and 5th-order compact schemes satisfy Eq. (16). According to Eq. (16), there is a minimum acceptable Re_∆x for each scheme combination. Some values for Re_{∆x min} are listed in Table 5. Therefore, for each scheme combination, Eq. (12) and Eq. (16) introduce an Re_{∆x max} and an Re_{∆x min}. Grid spacings chosen between these two restrictions satisfy these two requirements simultaneously, that is, .

3.4 Influence of time discretization

In obtaining the two criteria above, only space discretization is considered, while other effects of time discretization are neglected. However, when solving an equation, the influence of time discretization is significant. Therefore, a brief discussion of how time discretization affects these two criteria is included. The first criterion is based on comparing space discretization errors of different origins and hence is not affected by time discretization. The second criterion concerns high wavenumber modes, which are affected by both space discretization and time discretization. Thus, this section concentrates on discussing how time discretization affects the second requirement.

The full-discrete form of Eq. (2) is

(17)

Performing time marching with an mth Runge-Kutta scheme, the solution to Eq. (17) can be written as

(18)

Let F denote the discrete Fourier transformation operator and let Z denote the solution's amplitude change from t to (t+∆ t)^[16]. Then, we have

(19)

The expression for Z can be written as

(20)

Notice that c∆t/∆x is the Courant number and denote it by φ. Substituting c∆t/∆x with φ and ν∆t/∆x² with φ/Re_∆x yields

(21)

Based on the differentiation property of the Fourier transformation, we have

(22)

An mth Runge-Kutta scheme satisfies the mth expansion of a Taylor series, which is

(23)

Performing a Fourier transformation on both sides of Eq. (23), we have

(24)

Equation (24) represents the solution's amplitude change from the nth time step to the (n+1)th time step. The exact amplitude change is

and the effects of numerical methods ||N_full||, in accordance with the semi-discrete analysis, can be written as

(25)

Examples of ||N_full|| are plotted in Figs. 7 and 8, where φ=0.1, time marching is performed by a 4th Runge-Kutta scheme, and diffusion terms are discretized by a 4th-order central scheme. As seen in Figs. 6 and 7, the conclusions of the semi-discrete and full-discrete analysis are the same, that is, the combination of 5th-order upwind and 7th-order upwind amplifies high wavenumber modes. By comparing Figs. 8 and 7, it can be seen that high wavenumber modes can be dissipated by increasing Re_∆x. This is also consistent with the results in the previous discussion.

Increasing φ to 0.4 while keeping the other coefficients unchanged, ||N_full|| is plotted in Figs. 9 and 10. At φ=0.4, the values of ||N_full|| are higher than their counterparts at φ=0.1. Comparing Figs. 6 and 9, we can see that the full-discrete analysis draws different conclusions than the semi-discrete analysis. Thus, combinations acceptable in the semi-discrete case may be unacceptable in the full-discrete case.

In conclusion, the effects of time discretization can be summarized as follows: (ⅰ) Re_{∆x min} is higher when time discretization is considered. (ⅱ) Re_{∆x min} is higher when φ is higher, which means that satisfying both criteria becomes more difficult. (ⅲ) At low φ, the semi-discrete conclusion can be used directly, but at high φ, the effects of time discretization must be considered.

4 Numerical experiments and discussion 4.1 Practical diffusion accuracy

Adopt the convection-diffusion equation and analyze the wave packet propagation problem whose initial value is

(26)

The exact solution is

(27)

Let x=0.01, c=1, ν=0.002, k=80, and the corresponding Re_∆x=5 and α=0.8. Time marching is performed by a 4th Runge-Kutta scheme, and let φ=0.01 to avoid the influence of time discretization. Convection terms are discretized by the 2nd and 7th upwind schemes, and diffusion terms are discretized by the 2nd and 4th central schemes. The cell Reynolds number in this case is higher than Re_{∆x max} of the 2nd U-2nd C and 2nd U-4th C scheme combinations, but lower than Re_{∆x max} of 7th U-2nd C and 7th U-4th C scheme combinations, where U and C stand for "upwind" and "central", respectively. Numerical results are displayed in Fig. 11. For 2nd U-2nd C and 2nd U-4th C combinations, the numerical diffusion is larger than the diffusion discretization error, and hence the practical diffusion accuracy is determined by the convection scheme. Therefore, the adoption of high-order diffusion dicretization scheme cannot improve the overall results. Numerical experiments show that the results computed by 2nd U-4th C is worse than that by 2nd U-2nd C. For 7th U-2nd C and 7th U-4th C combinations, the practical diffusion accuracy is determined by the diffusion discretization scheme, and hence a higher-order diffusion scheme leads to better numerical results.

Similar results can be obtained by analyzing the Burgers' equation,

(28)

One exact solution is

(29)

Let x=0.01, ν=0.01, and the maximum Re_∆x=2. Numerical results of various scheme combinations are displayed in Fig. 12. For 2nd U-2nd C and 2nd U-4th C, the numerical diffusion is larger than the diffusion approximation error, and the adoption of 4th central scheme has few effects. For 5th U-2nd C, 7th U-2nd C, and 7th U-4th C, the numerical diffusion is smaller than the diffusion approximation error, and the practical diffusion accuracy is determined by the diffusion scheme. Therefore, the adoption of higher-order diffusion scheme leads to more accurate results (7th U-2nd C and 7th U-4th C), while the adoption of higher-order convection scheme does few effects (5th U-2nd C and 7th U-2nd C). Numerical results of both the convection-diffusion equation and the Burgers' equation are in agreement with the theoretical deductions presented before.

4.2 High wavenumber modes

Analyze the wave packet propagation problem by the convection-diffusion equation and adopt the same initial value as that in the last subsection. Let x=0.01, c=1, ν=0.008, k=80, 300, and the corresponding Re_∆x=1.25, α=0.8, and α=3.0. The time marching method and φ are the same as those in the previous case. The results computed by 7th U-2nd C and 7th U-4th C are displayed in Fig. 13. As Re_∆x=1.25 is lower than Re_{∆x min} for both combinations, both combinations amplify the mode at α=3.0, which is undesirable. Meanwhile, as Re_∆x=1.25 is smaller than Re_{∆x max} of both combinations, the adoption of high-order diffusion scheme at the interested wavenumber (α=0.8) improves numerical accuracy.

Figure 14 shows the numerical results of a combination that meets both criteria, where Re_∆x=10 and 7th U-2nd C are employed. The three wave packages correspond to α=3.0, α=0.8, and α=0.4, respectively. Numerical experiments show that this combination dissipates high wave modes (α=3.0) while preserves high accuracy in the interested wave range (α=0.8 and α=0.4).

5 Conclusions

In this paper, a new approach is put forward for choosing an ideal numerical scheme and grid size. Both diffusion accuracy and high wavenumber performance are considered. This new approach satisfies three requirements: (ⅰ) Maintain high resolution near the length scales of interest. (ⅱ) Maintain theoretical diffusion accuracy at length scales of interest. (ⅲ) Dissipate high wavenumber numerical disturbances. The requirements of this new approach can be summarized as

(30)

Analysis of time discretization shows that high φ leads to high Re_{∆x min}, making it harder to satisfy the last two criteria simultaneously. At low φ, the semi-discrete conclusion can be used directly, but at high φ, the effects of time discretization must be considered. Numerical results are consistent with theoretical conclusions.

The discussion in this paper can also provide insights for the Navier-Stokes equation, as the Reynolds number is the counterpart of c/ν and Re_∆x=Re∆x.