Introduction

UJPR

Universal Journal of Physics Research

Trend Research Publishing

10.31586/ujpr.2022.155

UJPR-155

Research Article

Compact Difference Schemes Combined with Runge-Kutta Methods for Solving Unsteady Convection-Diffusion Problems

Zhu

Zhenwei

1 Chen

Junjie

1 *

1Department of Energy and Power Engineering, School of Mechanical and Power Engineering, Henan Polytechnic University, Jiaozuo, Henan, China

*Corresponding author at: Department of Energy and Power Engineering, School of Mechanical and Power Engineering, Henan Polytechnic University, Jiaozuo, Henan, China

05 11 2021

1 1 05 11 2021 05 11 2021 05 11 2021 05 11 2021

2021

This work is licensed under the Creative Commons Attribution International License (CC BY). http://creativecommons.org/licenses/by/4.0/

The convection-diffusion equation is of primary importance in understanding transport phenomena within a physical system. However, the currently available methods for solving unsteady convection-diffusion problems are generally not able to offer excellent accuracy in both time and space variables. A procedure was given in detail to solve the one-dimensional unsteady convection-diffusion equation through a combination of Runge-Kutta methods and compact difference schemes. The combination method has fourth-order accuracy in both time and space variables. Numerical experiments were conducted and the results were compared with those obtained by the Crank-Nicolson method in order to check the accuracy of the combination method. The analysis results indicated that the combination method is numerically stable at low wave numbers and small CFL numbers. The combination method has higher accuracy than the Crank-Nicolson method.

Runge-Kutta methods; Compact finite difference schemes; Convection-diffusion equations; Stability analysis; High-order accuracy; Crank-Nicolson methods

Introduction

The convection-diffusion equation is a combination of the diffusion convection and diffusion equations, and describes physical phenomena where particles, energy, or other physical quantities are transferred inside a physical system due to convection and diffusion processes [ 1,2]. This equation is of primary importance in understanding transport phenomena within a physical system [ 3,4].

For a transport variable u, the one-dimensional unsteady advection-diffusion equation is usually written as

(1)

\frac{\partial u (x, t)}{\partial t} + c \frac{\partial u (x, t)}{\partial x} = α \frac{\partial^{2} u (x, t)}{\partial x^{2}}, (x, t) \in (0, L) \times (0, T),

with the initial condition

u (x, 0) = ϕ (x),

and the boundary conditions

\frac{\partial u (x, t)}{\partial t} + c \frac{\partial u (x, t)}{\partial x} = α \frac{\partial^{2} u (x, t)}{\partial x^{2}} ， (x, t) \in (0, L) \times (0, T) ， (x, t) \in (0, L) \times (0, T),

in which c is an arbitrary constant indicating the intensity of the convection, the diffusion coefficient α is positive constant, and g₀ and g₁ are assumed in this study to be smooth functions.

A variety of finite difference schemes have been presented to solve unsteady convection-diffusion problems approximately [ 5]. These schemes usually have first-order or second-order accurate in space variables, but have fallen short of expectations for convection-dominated flows if the mesh is not sufficiently refined [ 5]. Discretization with a higher order in space variables is usually associated with large stencils, thereby increasing the band-width of the resulting matrix [ 5,6]. A class of compact finite difference approximations have been recently developed to solve convection-diffusion problems [ 7,8,9]. These compact finite difference schemes have fourth-order accuracy in space variables, but fail when accuracy in time variables is most needed [ 10,11]. Additionally, most of these schemes have to compute the inverse of the block matrix [ 10,11]. It is therefore entirely necessary to develop an effective numerical scheme in order to obtain satisfactory results with a higher order accuracy in time variables and reasonable computational cost.

The primary focus of this study is on developing an effective method to achieve higher accuracy in time variables by a combination of Runge-Kutta methods and high-order compact finite difference discretization schemes. This combination method is then used to solve the one-dimensional unsteady convection-diffusion equation and has fourth-order accuracy in both time and space variables. The Runge-Kutta methods are introduced briefly in Section 2. A fourth-order compact finite difference approximation scheme is presented in Section 3 for the unsteady convection-diffusion equation. A stability analysis is performed on the fully discretized equation in Section 4 to check the stability of the combination method. Some numerical examples are presented in Section 5 and concluding remarks are finally given in Section 6.

Runge-Kutta methods

In numerical analysis, Runge-Kutta methods are a family of implicit and explicit iterative methods [ 12]. When Runge-Kutta methods are used to solve the convection-diffusion equation, time is treated as an independent variable in an ordinary differential equation, for example,

\frac{d Q}{d t} = R (t, Q),

in which Q is an unknown function of time, and R denotes the numerical approximation to the spatial derivatives. The Runge-Kutta method used in this study is given in the following form

Q^{n + 1} = Q^{n} + Δ t \hat{R} (Q^{n}, Δ t),

in which Δt is the time increment. The increment function $\hat{R} (Q^{n}, Δ t)$ is subdivided into N steps on the interval $t^{n} \leq t \leq t^{n + 1}$

Q^{1} = Q^{n} + Δ t (α_{11} R^{n}),

Q^{2} = Q^{n} + Δ t (α_{21} R^{n} + α_{22} R^{1}),

Q^{3} = Q^{n} + Δ t (α_{31} R^{n} + α_{32} R^{1} + α_{33} R^{2}),

Q^{n + 1} = Q^{n} + Δ t (α_{N 1} R^{n} + α_{N 2} R^{1} + \dots + α_{N N} R^{N - 1}),

in which the superscript n, 1, 2, ..., and n + 1 denote the time steps on the time interval $t^{n} \leq t_{1} \leq t_{2} \leq \dots \leq t_{N} \leq t^{n + 1}$ , and α_ij is the weighting factor for the step i and term j.

A four-step algorithm is given by

Q^{1} = Q^{n} + \frac{Δ t}{4} R^{n},

Q^{2} = Q^{n} + \frac{Δ t}{3} R^{1},

Q^{3} = Q^{n} + \frac{Δ t}{2} R^{2},

Q^{n + 1} = Q^{n} + Δ t R^{3} .

This algorithm is convenient to program and no intermediate solution needs to be stored [ 13].

A fourth-order compact finite difference scheme

The differential equation considered here takes the form

α \frac{\partial^{2} u (x)}{\partial x^{2}} - c \frac{\partial u (x)}{\partial x} = f (x), x \in (0, L),

with the boundary conditions

u (0) = g_{0}, u (L) = g_{1} .

A compact difference scheme with fourth-order accuracy can be obtained for Equation (14) as follows [ 14]: the interval 0 ≤ x ≤ 1 is subdivided into n equal subintervals by the grid points $x_{i} = i h$ , in which $h = \frac{1}{n}$ . The mesh function u(ih) is written as u_i at grid point x_i. The second-order central difference schemes for second and first derivatives of u can be written as $δ_{x}^{2} u = \frac{(u_{i + 1} - 2 u_{i} + u_{i - 1})}{h^{2}}$ and $δ_{x} u = \frac{(u_{i + 1} - u_{i - 1})}{2 h^{2}}$ , respectively. The following relation can be derived for Equation (14) at point x_i:

α δ_{x}^{2} u_{i} - c δ_{x} u_{i} - τ_{i} = f_{i},

in which

τ_{i} = \frac{h^{2}}{12} (α \frac{d^{4} u}{d x^{4}} - 2 c \frac{d^{3} u}{d x^{3}}) + O (h^{4}) .

The fourth and third derivatives of u in the above equation should be approximated in order to obtain a compact difference scheme with fourth-order accuracy.

Equation (14) gives:

\frac{d^{3} u}{d x^{3}} |_{i} = \frac{1}{α} (\frac{d f}{d x} + c \frac{d^{2} u}{d x^{2}}) |_{i} = \frac{1}{α} (δ_{x} f_{i} + c δ_{x}^{2} u_{i}) + O (h^{2})

From Equation (14) and Equation (18), it follows that

\frac{d^{4} u}{d x^{4}} |_{i} = \frac{1}{α} (\frac{d^{2} f}{d x^{2}} + c \frac{d^{3} u}{d x^{3}}) |_{i} = \frac{1}{α} δ_{x}^{2} f_{i} + \frac{c}{α^{2}} δ_{x} f_{i} + \frac{c^{2}}{α^{2}} δ_{x}^{2} u_{i}) + O (h^{2})

By substituting Equations (18) and (19) for the corresponding terms used in Equation (17), the solution is given by

τ_{i} = \frac{1}{12} (δ_{x}^{2} f_{i} - \frac{c}{α} δ_{x} f_{i} - \frac{c^{2}}{α} δ_{x}^{2} u_{i}) + O (h^{4})

A fourth-order compact finite difference scheme for Equation (14) can be obtained by substituting Equation (20) in Equation (16). This gives

(α + \frac{c^{2} h^{2}}{12 α}) δ_{x}^{2} u_{i} - c δ_{x} f_{i} = f_{i} + \frac{h^{2}}{12} (δ_{x}^{2} f_{i} - \frac{c}{α} δ_{x} f_{i}) + O (h^{4})

Here, two difference operators are defined as follows:

L_{x} = 1 + \frac{h^{2}}{12} (δ_{x}^{2} - \frac{c}{α} δ_{x}), A_{x} = - (α + \frac{c^{2} h^{2}}{12 α}) δ_{x}^{2} + c δ_{x}) .

Equation (21) can be formulated symbolically as

L_{x}^{- 1} A_{x} u_{i} = f_{i} + O (h^{4})

The above compact finite difference scheme is combined with the Runge-Kutta method described earlier to solve the one-dimensional advection-diffusion equation.

Stability analysis

A von Neumann linear stability analysis is performed to check the stability of the combination method. Let $u_{i}^{n} = b^{n} e^{I k h i}$ to be the value of $u_{i}^{n}$ at x_i , in which $I = \sqrt{- 1}$ , bⁿ is the amplitude at time level n, and k is the wave number. The amplification factor is defined as

g = \frac{b^{n + 1}}{b^{n}}

The amplification factor of Equation (13) discretized by the Runge-Kutta method can be obtained as follows [ 13]:

g = 1 + Z + \frac{1}{2} Z^{2} + \frac{1}{6} Z^{3} + \frac{1}{24} Z^{4}

The variable Z represents the spatial discretization applied to the convection-diffusion terms. The explicit, discretized form is given by

Z = [(24 \frac{F}{P e} + 2 F \cdot P e) (1 - \cos k) - 12 I F \sin k] / (10 + 2 \cos k - I P e \sin k)

in which F is the CFL number and is defined as

F = \frac{c Δ t}{h}

In Equation (26), Pe is the Péclet number and is defined as

P e = \frac{c h}{α}

When the Péclet number is high, the convection term dominates [ 15,16]. When the Péclet number is low the diffusion term dominates [ 15,16].

The amplification factor, as defined above, is rather complex. For stability, the amplification factor has to satisfy the relation $| g | \leq 1$ . When $| g | > 1$ , the combination method for solving the advection-diffusion equation is numerically unstable. Contours of the absolute value of the amplification factor are plotted in the k-F plane in Figures 1-4, respectively, with various Péclet numbers such as 1, 10, 100, and 1000. The results indicate that the combination method for solving the advection-diffusion equation is numerically stable at low wave numbers and small CFL numbers.

Figure 1

Contours of the absolute value of the amplification factor in the k-F plane with a Péclet number of 1.

Figure 2

Contours of the absolute value of the amplification factor in the k-F plane with a Péclet number of 10.

Figure 3

Contours of the absolute value of the amplification factor in the k-F plane with a Péclet number of 100.

Figure 4

Contours of the absolute value of the amplification factor in the k-F plane with a Péclet number of 1000.

Numerical experiments

In numerical analysis, the Crank-Nicolson method is a finite difference method used for numerically solving the heat equation and similar partial differential equations [ 17,18]. Some numerical results are presented in examples 5.1-5.3 to compare the present method, denoted by P.M, with the Crank-Nicolson method, denoted by C-N.M. It is assumed that all exact solutions of the convection-diffusion equations given below are known.

Example 5.1. Given a convection-diffusion equation

\frac{\partial u}{\partial t} + 0.1 \frac{\partial u}{\partial x} = 0.01 \frac{\partial^{2} u}{\partial x^{2}} ， (x, t) \in (0, 1) \times (0, T) .

the exact solution of the above equation is given by

u (x, t) = e^{5 x - (0.25 + 0.01 π^{2}) t} \sin π x

The initial and boundary conditions are defined so that the results are agreed well with the exact solution of the problem. The accuracy of the present method is compared with that of the Crank-Nicolson method. The results are presented inTable 1 for various values of t.Table 1 gives the absolute error along the x direction of the domain, where h is 0.005 and Δt is 0.001.

Table 1

Table 1. The absolute error for various values of t, where h is 0.005 and Δt is 0.001.

t	x = 0.1		x = 0.3		x = 0.5		x = 0.7		x = 0.9
	P.M	CN.M	P.M	CN.M	P.M	CN.M	P.M	CN.M	P.M	CN.M

0.2	6.579E-010	1.897E-005	1.623E-009	3.636E-005	2.108E-009	1.409E-005	2.725E-009	2.066E-004	5.991E-008	1.645E-003
0.4	1.134E-009	3.239E-005	3.028E-009	6.781E-005	3.932E-009	2.630E-005	5.146E-009	3.854E-004	2.455E-007	1.645E-003
0.6	1.474E-009	4.170E-005	4.232E-009	9.472E-005	5.501E-009	3.679E-005	8.038E-009	5.380E-004	5.145E-007	2.105E-003
0.8	1.720E-009	4.820E-005	5.248E-009	1.173E-004	6.838E-009	4.576E-005	1.294E-008	6.654E-004	8.263E-007	2.418E-003
1	1.896E-009	5.268E-005	6.088E-009	1.357E-004	7.951E-009	5.335E-005	2.11E-008	7.685E-004	1.158E-006	2.628E-003

Example 5.2. Given a convection-diffusion equation

\frac{\partial u}{\partial t} + 0.22 \frac{\partial u}{\partial x} = 0.5 \frac{\partial^{2} u}{\partial x^{2}}, (x, t) \in (0, 1) \times (0, T) .

the exact solution is given by

u (x, t) = e^{0.22 x - (0.0242 + 0.5 π^{2}) t} \sin π x

The initial and boundary conditions of the problem are defined so that the results are agreed well with the exact solution of the problem. The accuracy of the present method is compared with that of the Crank-Nicolson method. The results are presented inTable 2 for various values of t.Table 2 gives the absolute error along the x direction of the domain, where h is 0.01 and Δt is 0.5h².

Table 2

Table 2. The absolute error for various values of t, where h is 0.01 and Δt is 0.5h².

t	x = 0.1		x = 0.3		x = 0.5		x = 0.7		x = 0.9
	P.M	CN.M	P.M	CN.M	P.M	CN.M	P.M	CN.M	P.M	CN.M

0.2	1.551E-009	9.770E-006	6.236E-009	2.674E-005	1.655E-008	3.456E-005	3.943E-008	2.921E-005	8.450E-008	1.162E-005
0.4	6.331E-009	7.250E-006	2.185E-008	1.984E-005	4.491E-008	2.563E-005	7.997E-008	2.167E-005	1.306E-007	8.636E-006
0.6	1.052E-008	4.034E-006	3.432E-008	1.104E-005	6.418E-008	1.426E-005	1.025E-007	1.205E-005	1.507E-007	4.807E-006
0.8	1.299E-008	1.995E-006	4.146E-008	5.460E-006	7.459E-008	7.054E-006	1.136E-007	5.963E-006	1.593E-007	2.378E-006
1	1.425E-008	9.253E-007	4.504E-008	2.531E-006	7.965E-008	3.270E-006	1.18E-007	2.765E-006	1.629E-007	1.102E-006

Example 5.3. Given a convection-diffusion equation

\frac{\partial u}{\partial t} + 0.1 \frac{\partial u}{\partial x} = 0.2 \frac{\partial^{2} u}{\partial x^{2}}, (x, t) \in (0, 1) \times (0, T) .

the exact solution is given by

u (x, t) = e^{0.25 x - (0.0125 + 0.02 π^{2}) t} \sin π x

The initial and boundary conditions are defined so that the results are agreed well with the exact solution of the problem. The accuracy of the present method is compared with that of the Crank-Nicolson method. The results are presented inTable 3 for various values of t.Table 3 gives the absolute error along the x direction of the domain, where h is 0.01 and Δt is h².

Table 3

Table 3. The absolute error for various values of t, where h is 0.01 and Δt is h².

t	x = 0.1		x = 0.3		x = 0.5		x = 0.7		x = 0.9
	P.M	CN.M	P.M	CN.M	P.M	CN.M	P.M	CN.M	P.M	CN.M

0.2	3.677E-010	7.143E-006	1.143E-009	1.970E-005	2.687E-009	2.565E-005	9.496E-009	2.183E-005	3.559E-008	8.703E-006
0.4	9.510E-010	9.617E-006	3.896E-009	2.650E-005	1.115E-008	3.447E-005	2.988E-008	2.933E-005	7.205E-008	1.173E-005
0.6	2.363E-009	9.702E-006	9.198E-009	2.672E-005	2.285E-008	3.475E-005	5.017E-008	2.956E-005	9.962E-008	1.184E-005
0.8	4.306E-009	8.698E-006	1.565E-008	2.395E-005	3.480E-008	3.114E-005	6.758E-008	2.649E-005	1.197E-007	1.062E-005
1	6.350E-009	7.309E-006	2.205E-008	2.012E-005	4.558E-008	2.616E-005	8.167E-008	2.225E-005	1.342E-007	8.925E-006

The comparison results presented inTable 1,Table 2 andTable 3 indicate that the present method has higher accuracy than the Crank-Nicolson method.

Conclusions

The primary focus of this study is on developing an effective method to solve the one-dimensional unsteady convection-diffusion equation by a combination of Runge-Kutta methods and high-order compact finite difference discretization schemes. A high-order compact finite difference scheme is used to approximate the spatial derivatives, and Runge-Kutta methods are used to approximate the time derivative. This combination method has fourth-order accuracy in both time and space variables. The numerical results obtained by the combination method have been found to be in good agreement with the exact solutions. Additionally, the combination method has higher accuracy than the Crank-Nicolson method. The combination method has superiority over most of the other high-order compact finite difference schemes in terms of computational cost, and is highly efficient as there is no need to compute the inverse of the block matrix.

Acknowledgement

This work was supported by the National Natural Science Foundation of China (No. 51506048).

References [1]

G. Peng, Z. Gao, and X. Feng. An extremum-preserving finite volume scheme for convection-diffusion equation on general meshes. Applied Mathematics and Computation, Volume 380, 2020, Article Number: 125301.

[2]

L. Wang, X. Yang, H. Wang, Z. Chai, and Z. Wei. A modified regularized lattice Boltzmann model for convection-diffusion equation with a source term. Applied Mathematics Letters, Volume 112, 2021, Article Number: 106766.

[3]

J.N. Reddy and M. Martinez. A dual mesh finite domain method for steady-state convection-diffusion problems. Computers & Fluids, Volume 214, 2021, Article Number: 104760.

[4]

Y. Liu, Y. Feng, and R. Zhang. A high order conservative flux optimization finite element method for steady convec-tion-diffusion equations. Journal of Computational Physics, Volume 425, 2021, Article Number: 109895.

[5]

S. Karaa and J. Zhang. High order ADI method for solving unsteady convection-diffusion problems. Journal of Computational Physics, Volume 198, Issue 1, 2004, Pages 1-9.

[6]

A. Mohebbi and M. Dehghan. High-order compact solution of the one-dimensional heat and advection-diffusion equations. Applied Mathematical Modelling, Volume 34, Issue 10, 2010, Pages 3071-3084.

[7]

Y.-M. Wang and X. Wen. A compact exponential difference method for multi-term time-fractional convec-tion-reaction-diffusion problems with non-smooth solutions. Applied Mathematics and Computation, Volume 381, 2020, Article Number: 125316.

[8]

F. Zhang, X. Gao, and Z. Xie. Difference numerical solutions for time-space fractional advection diffusion equation. Boundary Value Problems, Volume 2019, Issue 1, 2019, Article Number: 14.

[9]

R. Zhang, Z. Wang, J. Liu, and L. Liu. A compact finite difference method for reaction-diffusion problems using compact in-tegration factor methods in high spatial dimensions. Advances in Difference Equations, Volume 2018, Issue 1, 2018, Article Number: 274.

[10]

Y. Zhang and Y. Zhang. A new difference scheme with high accuracy and absolute stability for solving convection-diffusion equations. Journal of Computational and Applied Mathematics, Volume 230, Issue 2, 2009, Pages 600-606.

[11]

M. Dehghan. Weighted finite difference techniques for the one-dimensional advection-diffusion equation. Applied Mathematics and Computation, Volume 147, Issue 2, 2004, Pages 307-319.

[12]

I. Fried. Numerical Solution of Differential Equations. Academic Press, Inc., London, United Kingdom, 1979, ISBN: 978-0-12-267780-9.

[13]

S.-T. Yu, K.-C. Hsieh, and Y.-L.P. Tsai. Simulating waves in flows by Runge-Kutta and compact difference schemes. AIAA Journal, Volume 33, Number 3, 1995, Pages 421-429.

[14]

A.L. Pardhanani, W.F. Spotz, and G.F. Carey. A stable multigrid strategy for convection-diffusion using high order compact discretization. Electronic Transactions on Numerical Analysis, Volume 6, 1997, Pages 211-223.

[15]

T. Han, J.A.C. Humphrey, and B.E. Launder. A comparison of hybrid and quadratic-upstream differencing in high Reynolds number elliptic flows. Computer Methods in Applied Mechanics and Engineering, Volume 29, Issue 1, 1981, Pages 81-95.

[16]

P.P. Chinchapatnam, K. Djidjeli, and P.B. Nair. Unsymmetric and symmetric meshless schemes for the unsteady convec-tion-diffusion equation. Computer Methods in Applied Mechanics and Engineering, Volume 195, Issues 19-22, 2006, Pages 2432-2453.

[17]

H. Fu, Y. Sun, H. Wang, and X. Zheng. Stability and convergence of a Crank-Nicolson finite volume method for space fractional diffusion equations. Applied Numerical Mathematics, Volume 139, 2019, Pages 38-51.

[18]

P. Wang and P. Huang. Convergence of the Crank-Nicolson extrapolation scheme for the Korteweg-de Vries equation. Applied Numerical Mathematics, Volume 143, 2019, Pages 88-96.