High-Effectiveness and -Accuracy Difference Scheme Based on Nonuniform Grids for Solving Convection–Diffusion Equations with Boundary Layers

Abstract

In this paper, some rational high-accuracy compact finite difference schemes on nonuniform grids (NRHOC) are introduced for solving convection–diffusion equations. The derived NRHOC schemes not only can suppress the oscillatory property of numerical solutions but can also obtain a high-accuracy approximate solution, and they can effectively solve the convection–diffusion problem with boundary layers by flexibly adjusting the discrete grid, which can be obtained with the singularity in the computational region. Three numerical experiments with boundary layers are conducted to verify the accuracy of the proposed NRHOC schemes. We compare the computed results with the analytical solutions, the results of the rational high-accuracy compact finite difference schemes on uniform grids (RHOC), and the other schemes in the literature. For all test problems, good computed results are obtained with the presented NRHOC schemes. It is shown that the presented NRHOC schemes have a better resolution for the solution of convection-dominated problems.

1. Introduction

The Navier–Stokes equations are very difficult to solve due to their highly nonlinear nature, and in particular, approximate solutions are required to have high accuracy. The convection–diffusion equation is a linearized model of the Navier–Stokes equations, and the study of its high-precision values can help develop algorithms for solving Navier–Stokes equations with high accuracy. At the same time, the convection–diffusion equation itself has a wide range of applications in certain fields of physics, chemistry, biology, environment, finance, biology, etc. [1,2,3]. Therefore, the study of high-precision numerical solutions to convection–diffusion equations is very important and has some academic research value. Researchers have developed various higher-accuracy, compact finite difference methods to solve kinds of convection–diffusion problems [4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23,24].

For the 1D (one-dimensional) convection–diffusion equation

$$\begin{array}{c}\hfill {\displaystyle -\epsilon u_{xx}+c\left(x\right)u_x=f\left(x\right)}\end{array}$$

(1)

the fourth-order high-accuracy compact finite difference scheme (referred to as FOC scheme) obtained by Spotz [7] using the original equation and truncation error correction method is

$$\begin{array}{c}\hfill -\epsilon \left(1+\frac{c^2{h}^2}{12{\epsilon }^2}\right){\delta }_x^2{u}_i+c{\delta }_x{u}_i=\left(1-\frac{c{h}^2}{12\epsilon }{\delta }_x+\frac{h^2}{12}{\delta }_x^2\right)f_i\end{array}$$

(2)

This scheme gives theoretical fourth-order accuracy when solving problems with smooth solutions in the computational region, but when the computational region contains boundary layers or singularities, the computational error of this scheme on coarse grids is large and does not achieve theoretical fourth-order accuracy. If the boundary-layer problem is to be solved with high accuracy, the mesh must be divided finely enough. Too fine grids result in huge computational cost, are memory-consuming, and greatly reduce computational efficiency. From the derivation process, it is not difficult to find that the FOC scheme is a strongly differential scheme. As the mesh Reynolds number, $Pe$$(Pe=ch/\epsilon )$, increases, the dissipation error increases rapidly, thus greatly reducing the resolution of the difference scheme. As we know, in numerical simulation, the resolution of the difference scheme is a key factor affecting the numerical accuracy and reliability of the numerical results. To reduce dissipation errors and improve computational accuracy, we have modified the FOC scheme to obtain a rational fourth-order compact difference scheme on a uniform grid (see Appendix A).

$$\begin{array}{c}\hfill -\alpha {\delta }_x^2{u}_i+c{\delta }_x{u}_i=f_i+c_1{f_x}_i+c_2{f_{xx}}_i\end{array}$$

(3)

where the coefficients are

$$\begin{array}{c}\hfill {\displaystyle \alpha =\epsilon \frac{1+\frac{c^2{h}^2}{6{\epsilon }^2}}{1+\frac{c^2{h}^2}{12{\epsilon }^2}},\begin{array}{cc}& \end{array}{C}_1=-\frac{h}{12}\frac{\frac{ch}{\epsilon }}{1+\frac{c^2{h}^2}{12{\epsilon }^2}},\begin{array}{cc}& \end{array}{C}_2=\frac{h^2}{12}\frac{1+\frac{c^2{h}^2}{6{\epsilon }^2}}{1+\frac{c^2{h}^2}{12{\epsilon }^2}}}\end{array}$$

(4)

For both the FOC scheme and the modified rational scheme in (3), Fourier error analysis reveals that the dissipation error is significantly reduced for the modified scheme for larger grid Reynolds numbers, but the error is still larger at medium and high wavelengths (see Appendix B). By monitoring the calculation process of specific examples, it is easy to see that if the computational grid is not sufficiently subdivided, there are no computational nodes or not enough computational nodes are distributed near and inside the boundary layer, which results in large computational errors in places with large gradient changes. Therefore, a more feasible and reasonable approach to improving the calculation accuracy of boundary-layer attachments is to distribute more calculation nodes in the area with large gradients or boundary layers, while fewer calculation nodes are distributed in the area with a small gradient or relatively gentle changes in physical quantities. This not only can consider the stability of the algorithm and the accuracy of the calculation results but can also reduce the calculation amount, which requires nonuniform mesh generation in the calculation area and that the calculation be carried out on a nonuniform grid. The main methods currently available to solve this type of problem are local grid refinement methods [25,26], mesh adaptive algorithms [27], coordinate transformation techniques [28,29,30], transform-free, nonuniform, grid-based direct discretization methods [31,32,33,34], and other new methods [35,36,37,38].

A fourth-order compact difference scheme for solving the two-dimensional convection–diffusion equation was first proposed by Gupta [4] using the original equation and the truncation error correction method. However, the difference scheme is a strongly dissipative scheme that performs poorly in solving large-gradient or boundary-layer problems. Kalita et al. [31] extended this fourth-order compact difference scheme to a nonuniform grid and developed a polynomial compact difference scheme on nonuniform grids, which improved the solution accuracy of some boundary-layer problems but did not completely improve the strong dissipation of the difference scheme. In this paper, we obtain fourth-order rational compact difference schemes based on nonuniform grids for 1D and 2D convection–diffusion equations using the ideas of the construction of modified differential Equation (3).

The remainder of this study is arranged into four sections as follows.

Section 2 develops NRHOC difference schemes for 1D convection–diffusion equations.

Section 3 constructs NRHOC difference schemes for 2D convection–diffusion equations.

Section 4 carries out numerical experiments to verify the feasibility of the presented NRHOC finite difference schemes.

Section 5 concludes the work of the paper.

2. RHOC Scheme on Nonuniform Grids (NRHOC Scheme) for 1D Convection–Diffusion Equations

Let us consider the following 1D convection–diffusion equation:

$$\begin{array}{c}\hfill \left\{\begin{array}{c}-\epsilon u_{xx}+c\left(x\right)u_x=f,\begin{array}{cc}& \end{array}a<x<b\hfill \\ u\left(a\right)=\varphi \left(x\right),u\left(b\right)=\psi \left(x\right)\hfill \end{array}\right.\end{array}$$

(5)

where $\epsilon$ is constant, with $\epsilon >0$, and is the diffusion coefficient; $c\left(x\right)$ is the convection coefficient; $f\left(x\right)$ is called the source term; $u\left(x\right)$ is the unknown quantity to be computed. Suppose that $c\left(x\right)$, $f\left(x\right)$, and $u\left(x\right)$ are all sufficiently smooth functions of x.

2.1. NRHOC Scheme for 1D Constant-Coefficient Convection–Diffusion Equations

If $c\left(x\right)$ is a constant, Equation (5) is

$$\begin{array}{c}\hfill -\epsilon u_{xx}+c{u}_x=f\left(x\right)\end{array}$$

(6)

First, the solution interval $[a,b]$ is subdivided arbitrarily into N sub-intervals, and the inserted computational nodes are $a=x_0<x_1<x_2<\dots <x_N=b$ in order. We define the step lengths of the left and right sides of node $x_i$ as $h_L=x_i-x_{i-1}$ and $h_R=x_{i+1}-x_i,i=0,1,2,\dots ,N$, respectively.

Assuming that solution u is smooth enough in the computational region, the approximate formula for calculating the first- and second-order derivatives at node $x_i$ is derived using Taylor’s formula as

$$\begin{array}{c}\hfill {\left(\frac{\partial u}{\partial x}\right)}_i={\delta }_x{u}_i-{\chi }_1{\left(\frac{{\partial }^{3}u}{\partial x^3}\right)}_i-{\chi }_2{\left(\frac{{\partial }^{4}u}{\partial x^4}\right)}_i+O\left(\eta \right)\end{array}$$

(7)

$$\begin{array}{c}\hfill {\left(\frac{{\partial }^{2}u}{\partial x^2}\right)}_i={\delta }_x^2{u}_i-{\chi }_3{\left(\frac{{\partial }^{3}u}{\partial x^3}\right)}_i-{\chi }_4{\left(\frac{{\partial }^{4}u}{\partial x^4}\right)}_i+O\left(\omega \right)\end{array}$$

(8)

where

$$\begin{array}{c}\hfill {\chi }_1=\frac{h_{\mathrm{L}}{h}_R}{6},{\chi }_2=\frac{h_L{h}_R\left(h_R-h_L\right)}{24},{\chi }_3=\frac{h_{\mathrm{R}}-h_L}{3},{\chi }_4=\frac{h_R^2-h_L{h}_{\mathrm{R}}+h_L^2}{12}\end{array}$$

(9)

$$\begin{array}{c}\hfill \eta =h_{\mathrm{L}}{h}_R\left(h_R^2-h_{\mathrm{L}}{h}_R+h_L^2\right),\omega =\frac{\left(h_{\mathrm{R}}-h_L\right)\left(h_R^2+h_L^2\right)}{60}\end{array}$$

(10)

$$\begin{array}{c}\hfill {\delta }_x{u}_i=\frac{h_L^2{u}_{i+1}+\left(h_R^2-h_L^2\right)u_i-h_R^2{u}_{i-1}}{h_{\mathrm{L}}{h}_R\left(h_{\mathrm{L}}+h_R\right)}\end{array}$$

(11)

$$\begin{array}{c}\hfill {\delta }_x^2{u}_i=\frac{2\left[h_L{u}_{i+1}-\left(h_R+h_L\right)u_i+h_R{u}_{i-1}\right]}{h_{\mathrm{L}}{h}_R\left(h_{\mathrm{L}}+h_R\right)}\end{array}$$

(12)

where ${\delta }_x{u}_i$ and ${\delta }_x^2{u}_i$ are the second-order central difference operators on nonuniform grids, respectively.

By substituting (7) and (8) in (6), we yield the following equation:

$$\begin{array}{c}\hfill -\epsilon {{\delta }^2}_x{u}_i+c_i{\delta }_x{u}_i+\left(\epsilon {\chi }_3-c{\chi }_1\right){\left(\frac{{\partial }^{3}u}{\partial x^3}\right)}_i+\left(\epsilon {\chi }_4-c{\chi }_2\right){\left(\frac{{\partial }^{4}u}{\partial x^4}\right)}_i+o\left(\eta \right)+o\left(\omega \right)=f_i\end{array}$$

(13)

The derivative function is calculated directly from Equation (6), and using the original equation gives

$$u_{xxx}=\frac{c}{\epsilon }{u}_{xx}-\frac{1}{\epsilon }{f}_x$$

$$u_{xxxx}=\frac{c^3}{{\epsilon }^3}{u}_x-\frac{c^2}{{\epsilon }^3}f-\frac{c}{{\epsilon }^2}{f}_x-\frac{1}{\epsilon }{f}_{xx}$$

And by substituting this in (13), we obtain

$$\begin{array}{c}\hfill \begin{array}{c}-\epsilon {{\delta }^2}_x{u}_i+c_i{\delta }_x{u}_i+\left(-\frac{c^2}{\epsilon }{\chi }_1+c{\chi }_3\right)u_{xxi}+\left(\frac{c^3}{{\epsilon }^2}{\chi }_4-\frac{c^4}{{\epsilon }^3}{\chi }_2\right)u_{xi}+o\left(\eta \right)+o\left(\omega \right)\hfill \\ =\left(1-\frac{c^3}{{\epsilon }^3}{\chi }_2+\frac{c^2}{{\epsilon }^2}{\chi }_4\right)f_i+\left({\chi }_3-\frac{c}{\epsilon }{\chi }_1+\frac{c}{\epsilon }{\chi }_4-\frac{c^2}{{\epsilon }^2}{\chi }_2\right)f_{xi}+\left({\chi }_4-\frac{c}{\epsilon }{\chi }_2\right)f_{xxi}\hfill \end{array}\end{array}$$

(14)

By substituting (7) and (8) in (14), we obtain

$$\begin{array}{c}\hfill \begin{array}{c}-\epsilon \left(1+\frac{c^2}{{\epsilon }^2}{\chi }_1-\frac{c}{\epsilon }{\chi }_3\right){{\delta }^2}_x{u}_i+c_i\left(1-\frac{c^3}{{\epsilon }^3}{\chi }_2+\frac{c^2}{{\epsilon }^2}{\chi }_4\right){\delta }_x{u}_i+{\kappa }_1{u}_{xxxi}+o\left(\eta \right)+o\left(\omega \right)\hfill \\ =\left(1-\frac{c^3}{{\epsilon }^3}{\chi }_2+\frac{c^2}{{\epsilon }^2}{\chi }_4\right)f_i+\left({\chi }_3-\frac{c}{\epsilon }{\chi }_1+\frac{c}{\epsilon }{\chi }_4-\frac{c^2}{{\epsilon }^2}{{\chi }_2}_5\right)f_{xi}+\left({\chi }_4-\frac{c}{\epsilon }{\chi }_2\right)f_{xxi}\hfill \end{array}\end{array}$$

(15)

where

$$\begin{array}{c}\hfill \begin{array}{c}{\eta }_1=\left(\frac{c^3}{{\epsilon }^2}{\chi }_4-\frac{c^4}{{\epsilon }^3}{\chi }_2\right)\eta \begin{array}{cc},& {\omega }_1=\left(-\frac{c^2}{\epsilon }{\chi }_1+c{\chi }_3\right)\omega \end{array}\hfill \\ {\kappa }_1=-{\chi }_3\left(-\frac{c^2}{\epsilon }{\chi }_1+c{\chi }_3\right)-{\chi }_1\left(\frac{c^3}{{\epsilon }^2}{\chi }_4-\frac{c^4}{{\epsilon }^3}{\chi }_2\right)\hfill \end{array}\end{array}$$

(16)

By replacing $u_{xxx}$ in (15) with $u_{xxx}=\frac{c}{\epsilon }{u}_{xx}-\frac{1}{\epsilon }{f}_x$ and omitting the high-order terms, we obtain

$$\begin{array}{c}-\epsilon \left(1+\frac{c^2}{{\epsilon }^2}{\chi }_1-\frac{c}{\epsilon }{\chi }_3\right){\delta }_x^2{u}_i+c\left(1-\frac{c^3}{{\epsilon }^3}{\chi }_2+\frac{c^2}{{\epsilon }^2}{\chi }_4\right){\delta }_x^{}{u}_i+\left({\kappa }_1+\kappa \epsilon \right)u_{xxxxi}+o\left({\eta }_1+{\omega }_1\right)\hfill \\ =\left(1-\frac{c^3}{{\epsilon }^3}{\chi }_2+\frac{c^2}{{\epsilon }^2}{\chi }_4\right)f_i+\left({\chi }_3-\frac{c}{\epsilon }{\chi }_1+\frac{c}{\epsilon }{\chi }_4-\frac{c^2}{{\epsilon }^2}{\chi }_2\right)f_{xi}+\left({\chi }_4-\frac{c}{\epsilon }{\chi }_2-\kappa \right)f_{xxi}\hfill \end{array}$$

Then, the NRHOC difference scheme of constant-coefficient convection–diffusion Equation (6) is yielded.

$$\begin{array}{c}\hfill -\alpha {\delta }_x^2{u}_i+c_i{\delta }_x{u}_i=f_i+c_1{f}_{xi}+c_2{f}_{xxi}\end{array}$$

(17)

where

$$\begin{array}{c}\hfill \begin{array}{c}{\displaystyle \alpha =\epsilon \frac{1+\frac{c^2}{{\epsilon }^2}{\chi }_1-\frac{c}{\epsilon }{\chi }_3}{1+\frac{c^3}{{\epsilon }^3}{\chi }_2+\frac{c^2}{{\epsilon }^2}{\chi }_4},}\hfill \\ \\ {\displaystyle c_1=\epsilon \frac{-\frac{c}{\epsilon }{\chi }_1-\frac{c^2}{{\epsilon }^2}{\chi }_2+{\chi }_3+\frac{c}{\epsilon }{\chi }_4}{1-\frac{c^3}{{\epsilon }^3}{\chi }_2+\frac{c^2}{{\epsilon }^2}{\chi }_4},\begin{array}{cc}& \end{array}{c}_2=\frac{{\chi }_4-\frac{c}{\epsilon }{\chi }_2-\frac{{\kappa }_1}{c}}{1-\frac{c^3}{{\epsilon }^3}{\chi }_2+\frac{c^2}{{\epsilon }^2}{\chi }_4}}\hfill \end{array}\end{array}$$

(18)

$$\begin{array}{c}{\displaystyle \kappa =-{\chi }_3\left(-\frac{c}{\epsilon }{\chi }_1+{\chi }_3\right)-{\chi }_1\left(\frac{c^2}{{\epsilon }^2}{\chi }_4-\frac{c^3}{{\epsilon }^3}{\chi }_2\right)}\hfill \\ {\displaystyle {\chi }_1=\frac{h_L{h}_R}{6},{\chi }_2=\frac{h_L{h}_R\left(h_R-h_L\right)}{24},{\chi }_3=\frac{h_R-h_L}{3},{\chi }_4=\frac{h_R^2-h_L{h}_R+h_L^2}{12}}\hfill \end{array}$$

2.2. NRHOC Scheme for 1D Variable-Coefficient Convection–Diffusion Equations

In this subsection, the RHOC scheme for model Equation (5) will be derived based on difference scheme (17). Suppose that model Equation (5) has the following difference scheme:

$$\begin{array}{c}\hfill -{\alpha }_i{\delta }_x^2{u}_i+c_i{\delta }_x{u}_i=f_i+{c_1}_i{f_x}_i+{c_2}_i{f_{xx}}_i\end{array}$$

(19)

Using Taylor’s formula and the original equation (Equation (5)), the modified differential equation of difference scheme (19) can be derived as

$$\begin{array}{c}\hfill -\epsilon u_{xx}+c{u}_x=f+2c_2{c}_x{u}_{xx}+(c_1{c}_x+c_2{c}_{xx})u_x+O\left(\eta +\omega \right)\end{array}$$

(20)

where $\eta =h_L{h}_R\left(h_R^2-h_L{h}_R+h_L^2\right),\omega =\frac{\left(h_R-h_L\right)\left(h_R^2+h_L^2\right)}{60}$.

To obtain the higher-order difference scheme of the original equation (Equation (5)), its modified differential equation is obtained using residual correction:

$$\begin{array}{c}\hfill -B{u}_{xx}+D\left(x\right)u_x=f\left(x\right)\end{array}$$

(21)

where

$$\begin{array}{c}\hfill B=\epsilon -2c_2{c}_x,D\left(x\right)=c+c_1{c}_x+c_2{c}_{xx}\end{array}$$

(22)

By applying (17) to (21), we obtain the NRHOC scheme for one-dimensional variable-coefficient convection–diffusion Equation (5) as follows:

$$\begin{array}{c}\hfill -{\mathrm{\Lambda }}_i{\delta }_x^2{u}_i+D_i{\delta }_x{u}_i=f_i+C_{1i}{f_x}_i+C_{2i}{f_{xx}}_i\end{array}$$

(23)

where the coefficients are

$$\begin{array}{c}{\displaystyle {\mathrm{\Lambda }}_i=B_i\frac{1+\frac{{D_i}^2}{{B^2}_i}{\chi }_1-\frac{D_i}{B_i}{\chi }_3}{1-\frac{{D_i}^3}{{B_i}^3}{\chi }_2+\frac{{D_i}^2}{{B_i}^2}{\chi }_4},}\hfill \\ \\ {\displaystyle D_{1i}=\frac{-\frac{D_i}{B_i}{\chi }_1-\frac{{D_i}^2}{{B_i}^2}{\chi }_2+{\chi }_3+\frac{D_i}{B_i}{\chi }_4}{1-\frac{{D_i}^3}{{B_i}^3}{\chi }_2+\frac{{D_i}^2}{{B_i}^2}{\chi }_4},D_{2i}=\frac{{\chi }_4-\frac{D_i}{B_i}{\chi }_2-\frac{{\mathrm{K}}_{1i}}{D_i}}{1-\frac{{D_i}^3}{{B_i}^3}{\chi }_2+\frac{{D_i}^2}{{B_i}^2}{\chi }_4},}\hfill \\ \\ {\mathrm{K}}_{1i}=-{\chi }_3\left(-\frac{{D_i}^2}{B_i}{\chi }_1+c{\chi }_3\right)-{\chi }_1\left(\frac{{D_i}^3}{{B_i}^2}{\chi }_4-\frac{{D_i}^4}{{B_i}^3}{\chi }_2\right),\hfill \\ \\ {\displaystyle B_i=\epsilon -2c_{2i}{c}_{xi},\begin{array}{cc}& \end{array}{D}_i=c_i+c_{1i}{c}_{xi}+c_{2i}{c}_{xxi},}\hfill \end{array}$$

$$\begin{array}{c}\hfill \begin{array}{c}{\displaystyle c_{1i}=\frac{-\frac{c_i}{\epsilon }{\chi }_1-\frac{c_i^2}{{\epsilon }^2}{\chi }_2+{\chi }_3+\frac{c_i}{\epsilon }{\chi }_4}{1-\frac{c_i^3}{{\epsilon }^3}{\chi }_2+\frac{c_i^2}{{\epsilon }^2}{\chi }_4},\begin{array}{cc}& \end{array}{c}_{2i}=\frac{{\chi }_4-\frac{c_i}{\epsilon }{\chi }_2-\frac{K_{1i}}{{\epsilon }_1}}{1+\frac{c_i^3}{{\epsilon }^3}{\chi }_2+\frac{c_i^2}{{\epsilon }^2}{\chi }_4},}\hfill \\ \\ {\chi }_1=\frac{h_{\mathrm{L}}{h}_R}{6},{\chi }_2=\frac{h_L{h}_R\left(h_R-h_L\right)}{24},{\chi }_3=\frac{h_{\mathrm{R}}-h_L}{3},{\chi }_4=\frac{h_R^2-h_L{h}_{\mathrm{R}}+h_L^2}{12}.\hfill \end{array}\end{array}$$

(24)

Difference Equation (23) is the rational compact difference scheme on nonuniform grids of model Equation (5). In particular, if $h_L=h_{\mathrm{R}}$, the difference equation is degenerated into difference Equation (A13) on the grid.

3. NRHOC Scheme for 2D Convection–Diffusion Equations

Let us consider the following 2D convection–diffusion equation:

$$\begin{array}{c}\hfill \left\{\begin{array}{c}-{\epsilon }_1{u}_{xx}-{\epsilon }_2{u}_{yy}+c\left(x,y\right)u_x+d\left(x,y\right)u_y=f\left(x,y\right),\begin{array}{cc}& \end{array}a<x<b,c<y<d\hfill \\ u\left(a,y\right)={\varphi }_1\left(y\right),u\left(b,y\right)={\varphi }_2\left(y\right)\hfill \\ u\left(x,c\right)={\psi }_1\left(x\right),u\left(x,d\right)={\psi }_2\left(x\right)\hfill \end{array}\right.\end{array}$$

(25)

where ${\epsilon }_1,{\epsilon }_2$ are constant and ${\epsilon }_1>0,{\epsilon }_2>0$, $c,d$, and f are smooth enough functions of variables x and y. It is a linearized model equation of a two-dimensional, steady, incompressible Navier–Stokes equation.

First, the constant-coefficient problem is considered; then, the NRHOC scheme of Equation (25) with variable coefficients can be obtained by modifying the difference equation.

3.1. NRHOC Scheme for 2D Constant-Coefficient Convection–Diffusion Equations

Let us consider the following two-dimensional convection–diffusion equation with constant coefficients:

$$\begin{array}{c}\hfill -{\epsilon }_1{u}_{xx}-{\epsilon }_2{u}_{yy}+c{u}_x+d{u}_y=f(x,y)\end{array}$$

(26)

We rewrite Equation (26) into the following equivalent equations:

$$\begin{array}{c}\hfill \left\{\begin{array}{c}-{\epsilon }_1{u}_{xx}+c{u}_x=f_1(x,y)\hfill \\ -{\epsilon }_2{u}_{yy}+d{u}_y=f_2(x,y)\hfill \end{array}\right.\end{array}$$

(27)

where

$$\begin{array}{c}\hfill \left\{\begin{array}{c}{f}_1(x,y)=f(x,y)+{\epsilon }_2{u}_{yy}-d{u}_y\hfill \\ f_2(x,y)=f(x,y)+{\epsilon }_2{u}_{xx}-c{u}_x\hfill \end{array}\right.\end{array}$$

(28)

Applying difference scheme (17) to (27) and (28) yields

$$\begin{array}{c}\hfill \left\{\begin{array}{c}-{\alpha }_1{\delta }_x^2{u}_{ij}+c{\delta }_x{u}_{ij}=F_{1ij}\hfill \\ -{\alpha }_2{\delta }_y^2{u}_{ij}+d{\delta }_y{u}_{ij}=F_{2ij}\hfill \end{array}\right.\end{array}$$

(29)

where

$$\begin{array}{c}\hfill \left\{\begin{array}{c}{F}_{1ij}=f_{1ij}+c_1{f}_{1xij}+c_2{f}_{1xxij}\hfill \\ F_{2ij}=f_{2ij}+d_1{f}_{2yij}+d_2{f}_{2yyij}\hfill \end{array}\right.\end{array}$$

(30)

We use Equation (28) to calculate the derivative and substitute it into (29); then, the two obtained equations are added together to obtain the NRHOC scheme of two-dimensional constant-coefficient model Equation (26) as follows:

$$\begin{array}{c}\hfill \left\{-{\alpha }_x{\delta }_x^2-{\alpha }_y{\delta }_y^2+c{\delta }_x+d{\delta }_y+E{\delta }_x^2{\delta }_y^2+G{\delta }_x^2{\delta }_y+H{\delta }_y^2{\delta }_x+R{\delta }_x{\delta }_y\right\}{u}_{ij}=F_{ij}\end{array}$$

(31)

where

$$\begin{array}{c}{\displaystyle {\alpha }_x={\epsilon }_1\frac{1+\frac{c^2}{{\epsilon }_1^2}{\chi }_1-\frac{c}{{\epsilon }_1}{\chi }_3}{1-\frac{c^3}{{\epsilon }_1^3}{\chi }_2+\frac{c^2}{{\epsilon }_1^2}{\chi }_4},{\alpha }_y={\epsilon }_2\frac{1+\frac{d^2}{{\epsilon }_2^2}{\chi }_1-\frac{d}{{\epsilon }_2}{\chi }_3}{1-\frac{d^3}{{\epsilon }_2^3}{\chi }_2+\frac{d^2}{{\epsilon }_2^2}{\chi }_4}}\hfill \\ \\ {\displaystyle c_1=\frac{-\frac{c}{{\epsilon }_1}{\chi }_1-\frac{c^2}{{\epsilon }_1^2}{\chi }_2+{\chi }_3+\frac{c}{{\epsilon }_1}{\chi }_4}{1-\frac{c^3}{{\epsilon }_1^3}{\chi }_2+\frac{c^2}{{{\epsilon }_1}^2}{\chi }_4},c_2=\frac{{\chi }_4-\frac{c}{{\epsilon }_1}{\chi }_2-\frac{{\kappa }_1}{c}}{1-\frac{c^3}{{\epsilon }_1^3}{\chi }_2+\frac{c^2}{{\epsilon }_1^2}{\chi }_4}}\hfill \\ \\ {\displaystyle d_1=\frac{-\frac{d}{{\epsilon }_2}{\chi }_1-\frac{d^2}{{\epsilon }_2^2}{\chi }_2+{\chi }_3+\frac{d}{{\epsilon }_2}{\chi }_4}{1-\frac{d^3}{{\epsilon }_2^3}{\chi }_2+\frac{d^2}{{{\epsilon }_2}^2}{\chi }_4},d_2=\frac{{\chi }_4-\frac{d}{{\epsilon }_2}{\chi }_2-\frac{{\kappa }_2}{d}}{1-\frac{d^3}{{\epsilon }_2^3}{\chi }_2+\frac{d^2}{{\epsilon }_2^2}{\chi }_4}}\hfill \end{array}$$

$$\begin{array}{c}\hfill \begin{array}{c}{\displaystyle E=-c_2{\epsilon }_2-d_2{\epsilon }_1,G=c_{2}d-d_1{\epsilon }_1,H=d_{2}c-c_1{\epsilon }_2,R=c_{1}d+d_{1}c}\hfill \\ \\ {\displaystyle {\kappa }_1=-{\chi }_3\left(-\frac{c^2}{{\epsilon }_1}{\chi }_1+c{\chi }_3\right)-{\chi }_1\left(\frac{c^3}{{\epsilon }_1^2}{\chi }_4-\frac{c^4}{{\epsilon }_1^3}{\chi }_2\right)}\hfill \\ \\ {\displaystyle {\kappa }_2=-{\chi }_3\left(-\frac{d^2}{{\epsilon }_2}{\chi }_1+d{\chi }_3\right)-{\chi }_1\left(\frac{d^3}{{\epsilon }_2^2}{\chi }_4-\frac{d^4}{{\epsilon }_2^3}{\chi }_2\right)}\hfill \\ \\ {\displaystyle F_{i,j}=\left(1+c_1{\delta }_x+c_2{\delta }_x^2+d_1{\delta }_y+d_2{\delta }_y^2\right)f_{i,j}}\hfill \end{array}\end{array}$$

(32)

3.2. NRHOC Scheme for 2D Variable-Coefficient Convection–Diffusion Equations

Assume that model Equation (25) has a difference scheme at point $\left(x_i,y_j\right)$ like (31):

$$\begin{array}{c}\hfill \left\{-{\alpha }_{xij}{\delta }_x^2-{\alpha }_{yij}{\delta }_y^2+c{\delta }_x+d{\delta }_y+E{\delta }_x^2{\delta }_y^2+G{\delta }_x^2{\delta }_y+H{\delta }_y^2{\delta }_x+R{\delta }_x{\delta }_y\right\}{u}_{ij}=F_{ij}\end{array}$$

(33)

where the coefficients are

$$\begin{array}{c}{\displaystyle {\alpha }_{xij}={\epsilon }_1\frac{1+\frac{{c_{ij}}^2}{{{\epsilon }_1}^2}{\chi }_1-\frac{c_{ij}}{{\epsilon }_1}{\chi }_3}{1-\frac{c^3}{{{\epsilon }_1}^3}{\chi }_2+\frac{c^2}{{{\epsilon }_1}^2}{\chi }_4},{\alpha }_{yij}={\epsilon }_2\frac{1+\frac{{d_{ij}}^2}{{{\epsilon }_2}^2}{\chi }_1-\frac{d_{ij}}{{\epsilon }_2}{\chi }_3}{1-\frac{d^3}{{{\epsilon }_2}^3}{\chi }_2+\frac{d^2}{{{\epsilon }_2}^2}{\chi }_4},}\hfill \\ \\ {\displaystyle c_{1ij}=\frac{-\frac{c_{ij}}{{\epsilon }_1}{\chi }_1-\frac{{c_{ij}}^2}{{{\epsilon }_1}^2}{\chi }_2+{\chi }_3+\frac{c_{ij}}{{\epsilon }_1}{\chi }_4}{1-\frac{{c_{ij}}^3}{{{\epsilon }_1}^3}{\chi }_2+\frac{{c_{ij}}^2}{{{\epsilon }_1}^2}{\chi }_4},c_{2ij}=\frac{{\chi }_4-\frac{c_{ij}}{{\epsilon }_1}{\chi }_2-\frac{{\kappa }_{1ij}}{c_{ij}}}{1-\frac{{c_{ij}}^3}{{{\epsilon }_1}^3}{\chi }_2+\frac{{c_{ij}}^2}{{{\epsilon }_1}^2}{\chi }_4}}\hfill \\ \\ {\displaystyle d_{1ij}=\frac{-\frac{d_{ij}}{{\epsilon }_2}{\chi }_1-\frac{{d_{ij}}^2}{{{\epsilon }_2}^2}{\chi }_2+{\chi }_3+\frac{d_{ij}}{{\epsilon }_2}{\chi }_4}{1-\frac{{d_{ij}}^3}{{{\epsilon }_2}^3}{\chi }_2+\frac{{d_{ij}}^2}{{{\epsilon }_2}^2}{\chi }_4},d_2=\frac{{\chi }_4-\frac{d_{ij}}{{\epsilon }_2}{\chi }_2-\frac{{\kappa }_{2ij}}{d_{ij}}}{1-\frac{{d_{ij}}^3}{{{\epsilon }_2}^3}{\chi }_2+\frac{{d_{ij}}^2}{{{\epsilon }_2}^2}{\chi }_4}}\hfill \end{array}$$

$$\begin{array}{c}\hfill \begin{array}{c}{\displaystyle E_{ij}=-c_{2ij}{\epsilon }_2-d_{2ij}{\epsilon }_1,G_{ij}=c_{2ij}{d}_{ij}-d_{1ij}{\epsilon }_1}\hfill \\ \\ {\displaystyle H_{ij}=d_{2ij}{c}_{ij}-c_{1ij}{\epsilon }_2,R_{ij}=c_{1ij}{d}_{ij}+d_{1ij}{c}_{ij}}\hfill \\ \\ {\displaystyle {\kappa }_{1ij}=-{\chi }_3\left(-\frac{{c_{ij}}^2}{{\epsilon }_1}{\chi }_1+c_{ij}{\chi }_3\right)-{\chi }_1\left(\frac{{c_{ij}}^3}{{{\epsilon }_1}^2}{\chi }_4-\frac{{c_{ij}}^4}{{{\epsilon }_1}^3}{\chi }_2\right)}\hfill \\ \\ {\displaystyle {\kappa }_{2ij}=-{\chi }_3\left(-\frac{{d_{ij}}^2}{{\epsilon }_2}{\chi }_1+d_{ij}{\chi }_3\right)-{\chi }_1\left(\frac{{d_{ij}}^3}{{{\epsilon }_2}^2}{\chi }_4-\frac{{d_{ij}}^4}{{{\epsilon }_2}^3}{\chi }_2\right)}\hfill \\ \\ {\displaystyle F_{ij}=\left(1+c_{1ij}{\delta }_x+c_{2ij}{\delta }_x^2+d_{1ij}{\delta }_y+d_{2ij}{\delta }_y^2\right)f_{ij}}\hfill \end{array}\end{array}$$

(34)

Difference Equation (33) is expanded using Taylor’s formula to obtain the corresponding modified differential equation as follows:

$$\begin{array}{c}\hfill \begin{array}{c}{\displaystyle -a{u}_{xx}-b{u}_{yy}+c(x,y)u_x+d(x,y)u_y-(2c_2{d}_x+2d_2{c}_y)u_{xy}-2c_2{c}_x{u}_{xx}-2d_2{d}_y{u}_{yy}}\hfill \\ {\displaystyle \begin{array}{cc}& \end{array}-(c_1{c}_x+c_2{c}_{xx}+d_1{c}_y+d_2{c}_{yy})u_x-(c_1{d}_x+c_2{d}_{xx}+d_1{d}_y+d_2{d}_{yy})u_y=f(x,y)+o(\eta +\omega )}\hfill \end{array}\end{array}$$

(35)

where $\eta =h_L{h}_R\left(h_R^2-h_L{h}_R+h_L^2\right),\omega =\frac{\left(h_R-h_L\right)\left(h_R^2+h_L^2\right)}{60}$.

To obtain the difference scheme with higher-order accuracy, adding the term

$$\begin{array}{c}\hfill [(2c_2{d}_x+2d_2{c}_y)u_{xy}+2c_2{c}_x{u}_{xx}+2d_2{d}_y{u}_{yy}+(c_1{c}_x+c_2{c}_{xx}+d_1{c}_y+\\ \hfill d_2{c}_{yy})u_x+(c_1{d}_x+c_2{d}_{xx}+d_1{d}_y+d_2{d}_{yy})u_y]\end{array}$$

to Equation (25) yields its modified differential equation,

$$\begin{array}{c}\hfill -A{u}_{xx}-B{u}_{yy}+C(x,y)u_x+D(x,y)u_y=F(x,y)\end{array}$$

(36)

where

$$\begin{array}{c}\hfill \begin{array}{c}A={\epsilon }_1-2c_2{c}_x,C=c+c_1{c}_x+c_2{c}_{xx}+d_1{c}_y+d_2{c}_{yy}\hfill \\ B={\epsilon }_2-2d_2{d}_y,D=d+c_1{d}_x+c_2{d}_{xx}+d_1{d}_y+d_2{d}_{yy}\hfill \\ F=f-\left(2c_2{d}_x+2d_2{c}_y\right)u_{xy}\hfill \end{array}\end{array}$$

(37)

Then, we apply (31) to (36), and the NRHOC scheme of two dimensional variable-coefficient convection–diffusion model Equation (25) is obtained as follows:

$$\begin{array}{c}\hfill \begin{array}{c}\left\{-{\mathrm{\Lambda }}_{xij}{\delta }_x^2-{\mathrm{\Lambda }}_{yij}{\delta }_y^2+C_{ij}{\delta }_x+D_{ij}{\delta }_y+E_{ij}{\delta }_x^2{\delta }_y^2+G_{ij}{\delta }_x^2{\delta }_y+H_{ij}{\delta }_y^2{\delta }_x+R_{ij}{\delta }_x{\delta }_y\right\}{u}_{ij}\hfill \\ \begin{array}{cccc}& & & \end{array}=\left\{1+C_{1ij}{\delta }_x+C_{2ij}{\delta }_x^2+D_{1ij}{\delta }_y+D_{2ij}{\delta }_y^2\right\}{F}_{ij}\hfill \end{array}\end{array}$$

(38)

where the coefficients are

$$\begin{array}{c}{\displaystyle {\mathrm{\Lambda }}_{xij}=A_{ij}\frac{\left(1+\frac{{C_{ij}}^2}{A_{ij}^2}{\chi }_1-\frac{C_{ij}}{A_{ij}}{\chi }_3\right)}{\left(1-\frac{C_{ij}^3}{A_{ij}^3}{\chi }_2+\frac{C_{ij}^2}{A_{ij}^2}{\chi }_4\right)},{\mathrm{\Lambda }}_{yij}=B_{ij}\frac{\left(1+\frac{{D_{ij}}^2}{B_{ij}^2}{\chi }_1-\frac{D_{ij}}{B_{ij}}{\chi }_3\right)}{\left(1-\frac{{D_{ij}}^3}{B_{ij}^3}{\chi }_2+\frac{{D_{ij}}^2}{B_{ij}^2}{\chi }_4\right)},}\hfill \\ \\ {\displaystyle C_{1ij}=\frac{-\frac{C_{ij}}{A_{ij}}{\chi }_1-\frac{{C_{ij}}^2}{A_{ij}^2}{\chi }_2+{\chi }_3+\frac{C_{ij}}{A_{ij}}{\chi }_4}{1-\frac{{C_{ij}}^3}{A_{ij}^3}{\chi }_2+\frac{{C_{ij}}^2}{A^2}{\chi }_4},C_{2ij}=\frac{{\chi }_4-\frac{C_{ij}}{A_{ij}}{\chi }_2-\frac{{\kappa }_{1ij}}{C_{ij}}}{1-\frac{{C_{ij}}^3}{A_{ij}^3}{\chi }_2+\frac{{C_{ij}}^2}{A_{ij}^2}{\chi }_4}}\hfill \\ \\ {\displaystyle D_{1ij}=\frac{-\frac{D_{ij}}{B_{ij}}{\chi }_1-\frac{{D_{ij}}^2}{B_{ij}^2}{\chi }_2+{\chi }_3+\frac{D_{ij}}{B_{ij}}{\chi }_4}{1-\frac{{D_{ij}}^3}{B_{ij}^3}{\chi }_2+\frac{{D_{ij}}^2}{B_{ij}^2}{\chi }_4},D_{2i}=\frac{{\chi }_4-\frac{D_{ij}}{B_{ij}}{\chi }_2-\frac{{\kappa }_{2ij}}{D_{ij}}}{1-\frac{{D_{ij}}^3}{B_{ij}^3}{\chi }_2+\frac{{D_{ij}}^2}{B_{ij}^2}{\chi }_4}}\hfill \end{array}$$

$$\begin{array}{c}\hfill \begin{array}{c}{\displaystyle E_{ij}=-C_{2ij}{\mathrm{\Lambda }}_{yij}-D_{2ij}{\mathrm{\Lambda }}_{xij},G_{ij}=C_{2ij}{D}_{ij}-D_{1ij}{\mathrm{\Lambda }}_{xij},}\hfill \\ \\ {\displaystyle H_{ij}=D_{2ij}{C}_{ij}-C_{1ij}{\mathrm{\Lambda }}_{yij},R_{ij}=C_{1ij}{D}_{ij}+D_{1ij}{C}_{ij}+2\left(c_{2ij}{d}_{xij}+d_{2ij}{c}_{yij}\right)}\hfill \\ \\ {\displaystyle {\kappa }_{1ij}=-{\chi }_3\left(-\frac{{C_{ij}}^2}{{\mathrm{\Lambda }}_{xij}}{\chi }_1+C_{ij}{\chi }_3\right)-{\chi }_1\left(\frac{{C_{ij}}^3}{{{\mathrm{\Lambda }}_{xij}}^2}{\chi }_4-\frac{{C_{ij}}^4}{{{\mathrm{\Lambda }}_{xij}}^3}{\chi }_2\right)}\hfill \\ \\ {\displaystyle {\kappa }_{2ij}=-{\chi }_3\left(-\frac{{D_{ij}}^2}{{\mathrm{\Lambda }}_{yij}}{\chi }_1+D_{ij}{\chi }_3\right)-{\chi }_1\left(\frac{{D_{ij}}^3}{{{\mathrm{\Lambda }}_{yij}}^2}{\chi }_4-\frac{{D_{ij}}^4}{{{\mathrm{\Lambda }}_{yij}}^3}{\chi }_2\right)}\hfill \end{array}\end{array}$$

(39)

In particular, difference Equation (38) is the RHOC difference scheme on uniform grids of model Equation (25) when and only when $h_L=h_{\mathrm{R}}$.

4. Numerical Examples

In this section, we select several cases to further verify the high accuracy and validity of the NRHOC schemes and compare the results with those in the literature. For the one-dimensional linear problem, the Thomas algorithm is used to solve a system of tridiagonal linear equations directly. For one-dimensional nonlinear problems and two-dimensional problems, we use the BiCGSTab(2) iterative method [39] to solve the discrete algebraic equations. The iterative process takes zero initial value and terminates when the condition

$${∥u^{n+2}-u^n∥}_2\le {10}^{-14}$$

is satisfied.

Above, n is the iteration number. The convergence order is calculated as

$$Rate=\frac{\log \left(\mathrm{err}\left(N_1\right)/\phantom{\mathrm{err}\left(N_1\right)\mathrm{err}\left(N_2\right)}\mathrm{err}\left(N_2\right)\right)}{\log (N_2/\phantom{N_2{N}_1)}{N}_1)}$$

where E($N_1$) and E($N_2$) are the maximum absolute errors calculated for two different grids $N_1$ and $N_2$.

For the one-dimensional problem, it is assumed that the defined interval to be solved is $\left[X_{\min },X_{\max }\right]$, while the two-dimensional problem is the square region $\left[X_{\min },X_{\max }\right]\times \left[Y_{\min },Y_{\max }\right]$. We use the following grid generation function to obtain the computing grid nodes:

$$x_i=X_{\min }+i\frac{X_{\max }-X_{\min }}{N_x}+\frac{{\lambda }_x}{\pi }\sin \left(\pi ·i\frac{X_{\max }-X_{\min }}{N_x}\right)$$

$$y_j=Y_{\min }+j\frac{Y_{\max }-Y_{\min }}{N_y}+\frac{{\lambda }_y}{\pi }\sin \left(\pi ·j\frac{Y_{\max }-Y_{\min }}{N_y}\right)$$

where $N_x$ and $N_y$ are the numbers of sub-intervals along the x-direction and y-direction. ${\lambda }_x$ and ${\lambda }_y$ are the scaling parameters that control the density of grid points in the x-direction and y-direction, respectively. When $-1\le {\lambda }_x,{\lambda }_y<0$, the computational nodes on the obtained nonuniform grid are dense near the left or lower boundary points of the interval; when $0<{\lambda }_x,{\lambda }_y\le 1$, the computational nodes on the obtained nonuniform grid are dense near the right or upper boundary. The larger the absolute values of ${\lambda }_x$ or ${\lambda }_y$, the more grid points are concentrated in the area that needs to be dense, and vice versa.

When the grid-stretching parameter is equal to 0, the corresponding computational grid is a uniform grid. To represent uniformity, in the following numerical experiments, we use RHOC to represent the rational higher-order compact difference scheme on a unform grid.

Example 1. 

$-\epsilon u_{xx}+\frac{1}{1+x}{u}_x=f\left(x\right),0<\epsilon \le 1,0<x<1$.

For this equation, the boundary conditions are $u\left(0\right)=1+2^{-1/\phantom{1\epsilon }\epsilon },u\left(1\right)=e+2$. The exact solution is $u\left(x\right)=e^x+2^{-1/\phantom{1\epsilon }\epsilon }{(1+x)}^{1+1/\phantom{1\epsilon }\epsilon }$. $f\left(x\right)$ is determined by the exact solution, When $\epsilon$ is small, the solution to the equation has a boundary layer at $x=1$.

Table 1. Comparison of maximum absolute error and convergence order, Example 1.

$\mathit{\epsilon }$	Nodes	FOC [7]		NRHOC
		Error	Rate	Error	Rate	$\mathbf{\lambda }$
1	11	2.18 (−7)	–	4.18 (−8)	–	0.00
	21	1.38 (−8)	3.98	2.71 (−9)	3.94
	41	8.63 (−10)	4.00	1.71 (−10)	3.99
${10}^{-2}$	21	3.35 (−2)	–	6.50 (−3)	–	0.75
	41	2.83 (−3)	3.57	4.08 (−4)	3.99
	81	1.65 (−4)	4.10	2.61 (−5)	3.97
${10}^{-3}$	21	3.60 (−1)	–	1.76 (−1)	–	0.95
	41	7.72 (−1)	−1.10	1.35 (−2)	3.71
	81	3.05 (−1)	1.34	7.93 (−4)	4.08
${10}^{-5}$	201	4.52 (−1)	–	2.95 (−2)	–	1.00
	401	2.18 (−1)	1.05	1.88 (−3)	3.98
	801	1.03 (−1)	1.08	1.15 (−4)	4.03

shows that both the FOC scheme and the NRHOC scheme are stable and effective when $\epsilon =1,{10}^{-2}$, but when $\epsilon ={10}^{-3}$ and $\epsilon ={10}^{-5}$, the error of the FOC scheme is large and no longer decreases as the number of grids increases, while the NRHOC scheme can obtain very exact solutions with appropriate values of grid parameters and can reach the theoretical fourth order. From Figure 1 and Figure 2, it is seen that the NRHOC scheme has obvious advantages over the FOC scheme in terms of computational accuracy and resolution.

Example 2. 

$-\epsilon u_{xx}+u{u}_x=f\left(x\right),0<\epsilon \le 1,0<x<1.$

The exact solution is $u\left(x\right)=e^{-x}+e^{\frac{(x-1)(1+\epsilon )}{\epsilon }}$, and the boundary conditions and source item $f\left(x\right)$ are given by the exact solution. When $\epsilon$ is small, the solution has a boundary layer at $x=1$. We calculated the solution for parameter $\epsilon ={10}^{-1},{10}^{-2},{10}^{-3},{10}^{-4},{10}^{-5}$ using three difference schemes, FOC [7], RHOC, and NRHOC. As shown in

Table 2. Comparison of maximum absolute error and convergence order, Example 2.

Nodes	FOC [7]		RHOC		NRHOC
	Error	Rate	Error	Rate	Error	Rate
$\epsilon ={10}^{-1}$					$\lambda =0.20$
41	9.72 (−6)	–	1.16 (−5)	–	9.41 (−6)	–
81	6.28 (−7)	3.95	7.68 (−7)	3.92	6.00 (−7)	3.97
161	3.99 (−8)	3.98	4.89 (−8)	3.97	3.76 (−8)	3.99
321	2.50 (−9)	4.00	3.07 (−9)	4.00	2.35 (−9)	4.00
$\epsilon ={10}^{-2}$					$\lambda =0.85$
41	1.06 (−2)	–	2.03 (−2)	–	8.34 (−4)	–
81	1.30 (−3)	3.03	7.50 (−4)	4.76	4.92 (−5)	4.08
161	3.32 (−4)	3.03	3.74 (−4)	1.01	2.99 (−6)	4.04
321	2.44 (−5)	3.76	3.03 (−5)	3.62	1.86 (−7)	4.01
$\epsilon ={10}^{-3}$					$\lambda =0.95$
41	6.17 (+0)	–	4.26 (+0)	–	2.65 (−2)	–
81	1.93 (+0)	1.68	3.37 (+0)	0.33	1.69 (−3)	3.97
161	2.59 (−1)	2.89	6.68 (−1)	2.34	1.10 (−4)	3.93
321	1.52 (−2)	4.09	3.87 (−2)	4.11	6.82 (−6)	4.02
$\epsilon ={10}^{-4}$					$\lambda =1.00$
81	3.70 (+1)	–	1.93 (+1)	–	1.11 (−1)	–
161	1.83 (+1)	1.02	9.65 (+0)	1.00	6.09 (−3)	4.19
321	8.43 (+0)	1.12	5.02 (+0)	0.94	3.74 (−4)	4.02
641	3.03 (+0)	1.48	3.26 (+0)	0.62	2.38 (−5)	3.97
$\epsilon ={10}^{-5}$					$\lambda =1.00$
160	1.86 (+2)	–	9.65 (+1)	–	1.62 (−1)	–
320	9.32 (+1)	1.00	4.80 (+1)	1.01	8.24 (−3)	4.30
640	4.66 (+1)	1.00	2.39 (+1)	1.01	5.03 (−4)	4.03
1280	2.31 (+1)	1.01	1.20 (+1)	1.00	3.15 (−6)	4.00

, when $\epsilon ={10}^{-1},{10}^{-2}$, the computational accuracy of all schemes can reach the theoretical fourth-order accuracy, but when $\epsilon ={10}^{-3},{10}^{-4},{10}^{-5}$, the computational errors of the FOC scheme and the RHOC scheme are large, while the NRHOC scheme can still obtain the fourth-order accuracy solution with appropriate values of the grid-scaling parameters.

Figure 3 and Figure 4 give a comparison of the computed and exact solutions, and the error distribution over the computational region using the three difference schemes when $\epsilon ={10}^{-3}$ and $\epsilon ={10}^{-5}$, respectively. As shown in Figure 3 and Figure 4, the computational solutions obtained with the FOC scheme and the RHOC scheme all produce nonphysical oscillations near $x=1$, while the NRHOC scheme matches well with the exact solution because there are enough computational nodes distributed in the boundary layer.

Example 3. 

$-(u_{xx}+u_{yy})+\mathrm{Re}{u}_x=0,0\le x,y\le 1$.

The exact solution is

$$\begin{array}{c}{\displaystyle u(x,y)=\frac{e^{\frac{\mathrm{Re}x}{2}}\sin \left(\pi y\right)[2e^{-\frac{\mathrm{Re}x}{2}}\sinh \left(\sigma x\right)+\sinh \sigma (1-x)]}{\sinh \left(\sigma x\right)}}\hfill \\ \sigma =\sqrt{{\pi }^2+{\mathrm{Re}}^2/4}\hfill \end{array}$$

The boundary conditions are determined by the exact solution. When the Reynolds number ($Re$) is large, The solution to this equation has a boundary layer at $x=0$. We calculated the solutions for Reynolds numbers $Re=1,10,100,1000,5000,7500$ using the NRHOC scheme and the HOC scheme [31].

Table 3. Comparison of maximum absolute error and convergence order, Example 3.

Nodes	HOC [31]		NRHOC		Nodes	HOC [31]		NRHOC
	Error	Rate	Error	Rate		Error	Rate	Error	Rate
$\mathit{Re}=1$, $\lambda =0.00$					$\mathit{Re}=10$, $\lambda =0.25$
$11\times 11$	2.15 (−6)	–	2.15 (−6)	–	$11\times 11$	2.96 (−3)	–	2.22 (−3)	–
$21\times 21$	1.29 (−7)	4.07	1.29 (−7)	4.07	$21\times 21$	1.85 (−4)	4.00	1.37 (−4)	4.01
$41\times 41$	7.94 (−9)	4.02	7.94 (−9)	4.02	$41\times 41$	1.15 (−5)	4.00	8.49 (−6)	4.00
$\mathit{Re}=100$, $\lambda =0.85$					$\mathit{Re}=1000$, $\lambda =0.90$
$21\times 21$	2.23 (−1)	–	1.18 (−2)	–	$32\times 16$	7.82 (−1)	–	1.41 (−1)	–
$41\times 41$	5.07 (−2)	2.13	7.14 (−4)	4.04	$64\times 32$	6.06 (−1)	0.37	1.45 (−2)	3.29
$81\times 81$	5.26 (−3)	3.27	4.38 (−5)	4.03	$128\times 64$	3.65 (−1)	0.73	9.00 (−4)	4.01
$\mathit{Re}=5000$, $\lambda =0.945$					$\mathit{Re}=7500$, $\lambda =0.95$
$128\times 64$	8.16 (−1)	–	3.35 (−2)	–	$128\times 64$	8.74 (−1)	–	7.66 (−2)	–
$256\times 128$	6.66 (−1)	0.29	2.86 (−3)	3.55	$256\times 128$	7.62 (−1)	0.20	8.98 (−3)	3.09
$512\times 256$	4.43 (−1)	0.59	1.79 (−4)	4.00	$512\times 256$	5.81 (−1)	0.39	5.77 (−4)	3.96

lists the computational error and convergence order when different values of $Re$ are taken. The computational data in the table show that compared with the HOC scheme, the NRHOC scheme has higher computational accuracy.

5. Conclusions

In this paper, we develop rational fourth-order compact difference schemes on uniform and nonuniform grids for convection–diffusion equations.

(1) The fourth-order rational compact difference scheme on uniform grids (the RHOC scheme) for one-dimensional convection–diffusion equations was developed by modifying the strongly dissipative polynomial fourth-order compact difference scheme of the one-dimensional convection–diffusion equation in reference [7]. The dispersion and dissipation errors of the RHOC scheme are significantly smaller than those of the FOC scheme; especially, the larger Pe is, the more obvious the advantage of the RHOC scheme in terms of resolution.

(2) To efficiently solve the large-gradient-variation and boundary-layer problems, we extended the obtained difference scheme on uniform grids to nonuniform grids and constructed a rational higher-order compact difference scheme for one-dimensional convection–diffusion equations on nonuniform grids (the NRHOC scheme).

(3) A fourth-order rational compact difference scheme on nonuniform grids for two- dimensional convection–diffusion equations was developed by using the dimensionality reduction method and a partial differential equation correction technique.

(4) Several typical numerical cases were selected for solution to demonstrate the advantages of the NRHOC scheme proposed in this paper in terms of computational accuracy, validity, and stability. Numerical experiments show that the NRHOC scheme has better scale resolution and is more fit to solve convection-dominated problems precisely.

In follow-up research, we aim to focus on the following three aspects: (i) based on the present NRHOC difference scheme of convection–diffusion equations, conducting high-accuracy numerical simulations of flow and heat transfer problems [40]; (ii) combining the present scheme and adaptive methods to improve the computational accuracy of convection–diffusion with convective dominance [41]; (iii) theoretically analyzing the stability and convergence of the proposed differential scheme [42,43]. We look forward to reporting the results of the corresponding studies in the future.