An Explicit Fourth-Order Compact Numerical Scheme for Heat Transfer of Boundary Layer Flow

Abstract

The main contribution of this article is to propose a compact explicit scheme for solving time-dependent partial differential equations (PDEs). The proposed explicit scheme has an advantage over the corresponding implicit compact scheme to find solutions of nonlinear and linear convection–diffusion type equations because the implicit existing compact scheme fails to obtain the solution. In addition, the present scheme provides fourth-order accuracy in space and second-order accuracy in time, and is constructed on three grid points and three time levels. It is a compact multistep scheme and conditionally stable, while the existing multistep scheme developed on three time levels is unconditionally unstable for parabolic and considered a type of equations. The mathematical model of the heat transfer in a mixed convective radiative fluid flow over a flat plate is also given. The convergence conditions of dimensionless forms of these equations are given, and also the proposed scheme solves equations, and results are compared with two existing schemes. It is hoped that the results in the current report are a helpful source for future fluid-flow studies in an industrial environment in an enclosure area.

1. Introduction

The linear reaction convection diffusion equation has a broad range of application in astrophysics, aerospace sciences, biology, and environmental sciences. Energy movements [1,2], fluid clotting [3], and hemodynamics [4] are the areas with relevance to the convection–diffusion equation. The Black–Scholes equation [5,6] is a practical equation for options pricing, transformed into a reaction–convection–diffusion equation on a semi-unbounded domain [7,8].

The mathematical conversion of a linear reaction convection–diffusion equation into a nonstandard finite difference equation was revealed by [9]; however, [10] proposed a study that comprised the family of positivity-preserving finite-difference methods for the classical Fisher–Kolmogorov–Petrovsky–Piscounov equation.

The fourth-order compact finite difference scheme was studied by Sun and Zhang [11] using parameters such as space and time. Additionally, when $\alpha =0\&\alpha \ne 0$, Mohebbi and Dehghan [12] offered some compact methods that suggested a fourth-order compact finite difference in space and time components for Equation (A1). Moreover, Nove [13] studied a three-point third-order implicit finite difference method that is conditionally stable.

Consider the one-dimensional advection–diffusion equation

$$\frac{\partial u}{\partial t}+\alpha \frac{\partial u}{\partial x}=\beta \frac{{\partial }^{2}u}{\partial x^2} \left(x,t\right)\in \left[0,L]\times [0,T\right]$$

(1)

with initial conditions

$$u\left(x,0\right)=\varphi \left(x\right)$$

and boundary conditions

$$u\left(0,t\right)=g_0\left(t\right),u\left(L,t\right)=g_1 \left(t\right)\hspace{1em}t\in \left[0,T\right]$$

The advection–diffusion equation is significant in physical systems such as fluid dynamics [14,15,16]. Equation (1) describes partial differential equations in their one-dimensional version, explaining advection–diffusion of quantities such as mass, heat, energy, vorticity. In [17], the impact of different kinds of parameters on velocity and temperature profiles is given, and it was concluded that, in some cases, a boost of volume fraction of nanoparticle led to a decrease in heat transfer rate, and this was a contradiction to the results of other researchers. Influence of viscous dissipation and heat generation/absorption of non-Newtonian power-law fluid is given in [18]. The equations were solved numerically by implementing the Runge–Kutta–Fehlberg fourth-fifth-order method, and results displayed through graphs. The injection parameter escalation was responsible for the decay of the local Nusselt number for both Newtonian and non-Newtonian nanofluid flow. In existing published research, results were obtained, and those results were often displayed in graphs. Thus, in [19], it was mentioned that the performance of heat transfer of non-Newtonian nanofluid was better than that obtained by Newtonian nanofluid for the case of injection and an impermeable plate.

Using the equation in solving systems consists of difficult initials, and boundary conditions provide little or no knowledge regarding analytical solutions [14,20]. However, [21,22] reveal a significant effort to form stable and accurate numerical techniques for the solution of Equation (1). Effects of Endothermic/Exothermic chemical reaction of MHD free convection flow have been studies in [23]. In the presence of a non-uniform heat source/sink, the Carreau nanoparticle is studied by Khan et al. [24]. In this study, in the presence of suction, heat radiation, and a heat source/sink, Jamalaudin et al. [25] numerically investigate the stagnant convection mixture-point flow nanofluid over a vertical stretching/shrinking sheet.

On the other hand, compact finite difference schemes can be used to find numerical solutions of differential equations on a smaller number of grid points. An existing scheme was updated [26] to three dimensions for solving the three-dimensional convection–diffusion problem. The scheme was called an exponential higher-order compact finite difference scheme. It provided an exact solution of homogeneous convection–diffusion problems with constant coefficients. A numerical scheme [27] consisted of implicit Euler for time discretization, and hybrid scheme for spatial discretization was studied. Since the solution of the type of equations considered was exhibited as a boundary layer, a piecewise uniform Shishkin mesh was considered for discretization in the spatial direction. A new high-order compact ADI method [28] has been considered for solving the time-dependent convection–diffusion equation. In addition, a correction technique was applied for reducing the error of the splitting term. Two numerical schemes, named the upwind-penalty method and SUPG method, were introduced in [29]. The stability analysis was given, and a comparison with the standard SUPG finite element method was given. Most recent research on the proposed numerical schemes can be found in [30,31].

The proposed compact scheme is constructed on three grid points that can be used to find solutions of convection–diffusion type problems. The scheme is explicit and can be used in some situations to obtain the solution without using the additional iterative scheme to solve discretized/difference equations. The existing corresponding implicit discretization is modified, and so, in this manner, an explicit scheme containing a second-order partial derivative in time term is obtained. The existing implicit compact scheme consisted of the mixed partial derivative term in space and time, and this term was responsible for establishing an implicit compact scheme rather than an explicit one. In this contribution, however, this mixed derivative term is turned into a second-order partial derivative in time term, and this term is discretized using a second-order standard–classical central difference formula constructed on three time levels. The first-order partial derivative in the time term is discretized using the standard/classical difference formula. The second-order difference formula is used to obtain the second accurate scheme in time and the fourth order in space. This second-order difference formula, however, can be replaced with the first-order difference formula to get the first-order accuracy in the time scheme. However, the fourth-order space accuracy can still be obtained.

Construction of Compact Scheme

Consider the convection–diffusion equation of the form

$$\frac{\partial U}{\partial t}+\frac{\partial U}{\partial x}=\frac{{\partial }^{2}U}{\partial x^2}$$

(2)

In the first part of the scheme, Equation (1) is discretized as

$$\frac{u_i^{n+1}-u_i^{n-1}}{2\mathsf{\Delta }t}+\frac{u_{i+1}^n-u_{i-1}^n}{2\mathsf{\Delta }\mathrm{x}}=\frac{u_{i+1}^n-2u_i^n+u_{i-1}^n}{{\left(\mathsf{\Delta }x\right)}^2}$$

(3)

Equation (3) is the classical/standard discretization of Equation (2). To use a compact scheme, one remainder term for space discretization of Equation (3) is considered as

$$\frac{u_i^{n+1}-u_i^{n-1}}{2\mathsf{\Delta }t}+\frac{u_{i+1}^n-u_{i-1}^n}{2\mathsf{\Delta }\mathrm{x}}-\frac{{\left(\mathsf{\Delta }x\right)}^2}{6}{u}_{xxx}=\frac{u_{i+1}^n-2u_i^n+u_{i-1}^n}{{\left(\mathsf{\Delta }x\right)}^2}-\frac{{\left(\mathsf{\Delta }x\right)}^2}{12}{u}_{xxxx}$$

(4)

Since the third and fourth partial derivatives of $u$ arenot available, these derivatives are found first. Therefore, finding the third and fourth-order derivative of Equation (2) yields

$$\frac{{\partial }^{3}u}{\partial x^3}=\frac{{\partial }^{2}u}{\partial x\partial t}+\frac{{\partial }^{2}u}{\partial x^2}$$

(5)

$$\frac{{\partial }^{4}u}{\partial x^4}=\frac{{\partial }^{3}u}{\partial x^2\partial t}+\frac{{\partial }^{3}u}{\partial x\partial t}+\frac{{\partial }^{2}u}{\partial x^2}$$

(6)

Substituting the third- and fourth-order derivative into Equation (4) gives a compact scheme that is implicit. Instead of this, the fourth-order partial derivative in Equation (6) can be expressed as

$$\frac{{\partial }^{4}u}{\partial x^4}=\frac{{\partial }^{2}u}{\partial t^2}+\frac{{\partial }^{2}u}{\partial x\partial t}+\frac{{\partial }^{2}u}{\partial t\partial x}+\frac{{\partial }^{2}u}{\partial x^2}$$

(7)

Substituting Equations (5) and (7) into (4) yields

$$\frac{u_i^{n+1}-u_i^{n-1}}{2\mathsf{\Delta }t}+\frac{u_{i+1}^n-u_{i-1}^n}{2\mathsf{\Delta }\mathrm{x}}-\frac{{\left(\mathsf{\Delta }x\right)}^2}{6}\left(\frac{{\partial }^{2}u}{\partial x\partial t}+\frac{{\partial }^{2}u}{\partial x^2}\right)=\frac{u_{i+1}^n-2u_i^n+u_{i-1}^n}{{\left(\mathsf{\Delta }x\right)}^2}-\frac{{\left(\mathsf{\Delta }x\right)}^2}{12}\left(\frac{{\partial }^{2}u}{\partial t^2}+2\frac{{\partial }^{2}u}{\partial t\partial x}+\frac{{\partial }^{2}u}{\partial x^2}\right)$$

(8)

Simplifying Equation (8) gives

$$\frac{u_i^{n+1}-u_i^{n-1}}{2\mathsf{\Delta }t}+\frac{u_{i+1}^n-u_{i-1}^n}{2\mathsf{\Delta }\mathrm{x}}-\frac{{\left(\mathsf{\Delta }x\right)}^2}{6}\frac{{\partial }^{2}u}{\partial x^2}=\frac{u_{i+1}^n-2u_i^n+u_{i-1}^n}{{\left(\mathsf{\Delta }x\right)}^2}-\frac{{\left(\mathsf{\Delta }x\right)}^2}{12}\left(\frac{{\partial }^{2}u}{\partial t^2}+\frac{{\partial }^{2}u}{\partial x^2}\right).$$

(9)

In operator form, Equation (9) is expressed as

$${\delta }_t{u}_i^n+{\delta }_x{u}_i^n-\frac{{\left(\mathsf{\Delta }x\right)}^2}{6}{\delta }_x^2{u}_i^n={\delta }_x^2{u}_i^n-\frac{{\left(\mathsf{\Delta }x\right)}^2}{12}\left({\delta }_t^2{u}_i^n+{\delta }_x^2{u}_i^n\right)$$

(10)

where ${\delta }_t{u}_i^n=\frac{u_i^{n+1}-u_i^{n-1}}{2\mathsf{\Delta }t}$, ${\delta }_x{u}_i^n=\frac{u_{i+1}^n-u_{i-1}^n}{2\mathsf{\Delta }t}$, ${\delta }_x^2{u}_i^n=\frac{u_{i+1}^n-2u_i^n+u_{i-1}^n}{{\left(\mathsf{\Delta }x\right)}^2}$ and ${\delta }_t^2{u}_i^n=\frac{u_i^{n+1}-2u_i^n+u_i^{n-1}}{{\left(\mathsf{\Delta }t\right)}^2}$.

Equation (10) is an explicit compact difference equation that is fourth order in space and second order in time, and it is constructed on three time levels. Thus, it requires an additional scheme evaluated on the first time level.

2. Stability Analysis

The Von Neumann stability criterion is employed to check the stability of the proposed explicit scheme. According to this criterion, some transformations are considered first as

$$u_i^{n\pm 1}=E^{n\pm 1}{e}^{iI\varphi },u_{i\pm 1}^n=E^n{e}^{\left(i\pm 1\right)I\varphi }$$

(11)

Employing transformation (11) into Equation (9) and dividing by $e^{iI\varphi }$ gives

$$\begin{gathered} \frac{E^{n+1}-E^{n-1}}{2\mathsf{\Delta }t}+\frac{\left(e^{I\varphi }-e^{-I\varphi }\right)E^n}{2\mathsf{\Delta }\mathrm{x}}-\frac{{\left(\mathsf{\Delta }x\right)}^2}{6}\left(\frac{e^{I\varphi }-2+e^{-I\varphi }}{{\left(\mathsf{\Delta }x\right)}^2}\right) \\ =\frac{\left(e^{I\varphi }-2+e^{-I\varphi }\right)E^n}{{\left(\mathsf{\Delta }x\right)}^2}-\frac{{\left(\mathsf{\Delta }x\right)}^2}{12}\left(\frac{E^{n+1}-2E^n+E^{n-1}}{{\left(\mathsf{\Delta }t\right)}^2}+\frac{\left(e^{I\varphi }-2+e^{-I\varphi }\right)E^n}{{\left(\mathsf{\Delta }x\right)}^2}\right) \end{gathered}$$

(12)

Equation (12) can be expressed as

$$E^{n+1}=A{E}^n+B{E}^{n-1}$$

(13)

where

$$\begin{gathered} A=\frac{\mathrm{bd}\left(-6cdIsin\theta +\mathsf{\Delta }td\left(cos\theta -1\right)+12d^2\left(cos\theta -1\right)+1\right)}{1+bd} \\ B=\frac{1-\frac{1}{bd}}{1+\frac{1}{bd}}=\frac{bd-1}{bd+1} \end{gathered}$$

Since the proposed scheme is constructed on the three time level, it requires one more equation, which can be constructed as

$$E^n=E^n+O{E}^{n-1}$$

Equations (13) and (14) are expressed as

$$\left[\begin{array}{c}{E}^{n+1}\\ E^n\end{array}\right]=\left[\begin{array}{cc}A& B\\ 1& 0\end{array}\right]\left[\begin{array}{c}{E}^n\\ E^{n-1}\end{array}\right]$$

(14)

Here, the amplification factor is a matrix. The eigenvalues of this matrix are used to construct stability conditions that can be expressed as

$$\begin{gathered} \left|{\lambda }_1\right|\le 1\mathrm{and}\left|{\lambda }_2\right|\le 1 \\ \left|\frac{A}{2}-\frac{1}{2}\sqrt{A^2+4B}\right|\le 1\mathrm{and}\left|\frac{A}{2}+\frac{1}{2}\sqrt{A^2+4B}\right|\le 1 \end{gathered}$$

(15)

Let $A=A_1+I{B}_1$

where

$$\begin{gathered} A_1=\mathsf{\Delta }td\left(cos\varphi -1\right)+12d^2\left(cos\varphi -1\right)+1 \\ {B}_1=-2csin\varphi \end{gathered}$$

Let

$$A=A_1+I{B}_1$$

(16)

$$A^2+4B=A_1^2-B_1^2+2I{A}_1{B}_1+4B=r\left(cos\psi +Isin\psi \right)$$

where $r^2={\left(A_1^2-B_1^2+4B\right)}^2+4A_1^2{B}_1^2$ and $\psi ={\tan }^{-1}\left(\frac{2A_1{B}_1}{A_1^2-B_1^2+4B}\right)$

$${\left(A^2+4B\right)}^{1/2}=r^{\frac{1}{2}}\left(cos\left(\frac{\psi +2k\pi }{2}\right)+Isin\left(\frac{\psi +2k\pi }{2}\right)\right),k=0,1$$

(17)

The condition ${\left|\frac{A}{2}-\frac{1}{2}{\left(A^2+4B\right)}^{\frac{1}{2}}\right|}^2\le 1$ implies

$${\left(\frac{A_1}{2}-r^{\frac{1}{2}}\cos \left(\frac{\psi +2k\pi }{2}\right)\right)}^2+{\left(\frac{B_1}{2}+r^{\frac{1}{2}}\sin \left(\frac{\psi +2k\pi }{2}\right)\right)}^2\le 1$$

(18)

In addition, the condition ${\left|\frac{A}{2}+\frac{1}{2}{\left(A^2+4B\right)}^{\frac{1}{2}}\right|}^2\le 1$ implies

$${\left(\frac{A_1}{2}+r^{\frac{1}{2}}\cos \left(\frac{\psi +2k\pi }{2}\right)\right)}^2+{\left(\frac{B_1}{2}+r^{\frac{1}{2}}\sin \left(\frac{\psi +2k\pi }{2}\right)\right)}^2\le 1$$

(19)

Consider the system of convection–diffusion equations

$$\frac{\partial u}{\partial t}+v\frac{\partial u}{\partial y}=\frac{{\partial }^{2}u}{\partial y^2}+a\theta .$$

(20)

$$\frac{\partial \theta }{\partial t}+v\frac{\partial \theta }{\partial y}=b\frac{{\partial }^2\theta }{\partial y^2}$$

(21)

The system of Equations (20) and (21) can be expressed as

$$\overline{A}\frac{\partial \mathit{U}}{\partial t}+B\frac{\partial \mathit{U}}{\partial y}=C\frac{{\partial }^2\mathit{U}}{\partial y^2}+D\mathit{U}$$

(22)

where $\overline{A}=\left[\begin{array}{cc}1& 0\\ 0& 1\end{array}\right],B=\left[\begin{array}{cc}v& 0\\ 0& v\end{array}\right],C=\left[\begin{array}{cc}1& 0\\ 0& b\end{array}\right],D=\left[\begin{array}{cc}0& a\\ 0& 0\end{array}\right]$, $\mathit{U}={\left[u,\theta \right]}^t$.

The proposed scheme requires a third- and fourth-order derivative, therefore

$$\begin{gathered} \frac{{\partial }^2\mathit{U}}{\partial y^2}=C^{-1}\left(\overline{A}{\mathit{U}}_t+B{\mathit{U}}_y-D\mathit{U}\right) \\ {\mathit{U}}_{yyy}=C^{-1}\left(\overline{A}{\mathit{U}}_{ty}+b{\mathit{U}}_{yy}-D{\mathit{U}}_y\right) \end{gathered}$$

(23)

$$\begin{gathered} {\mathit{U}}_{yyyy}=C^{-1}\left(\overline{A}{\mathit{U}}_{tyy}+b{\mathit{U}}_{yyy}-D{\mathit{U}}_{yy}\right) \\ =A_1{\mathit{U}}_{tt}+A_2{\mathit{U}}_{ty}+A_3{\mathit{U}}_{yy}-C^{-1}\overline{A}{C}^{-1}D{\mathit{U}}_t-B{C}^{-1}D{\mathit{U}}_y-D{\mathit{U}}_{yy} \end{gathered}$$

(24)

where $A_1=C^{-1}\overline{A}{C}^{-1}\overline{A},A_2=C^{-1}\overline{A}{C}^{-1}B+B{C}^{-1}\overline{A} and A_3=B{C}^{-1}B$.

Theorem 1.

The compact scheme for Equation (22) converges subject to the conditions

(i) 

Error at the first time level is bounded.

(ii) 

$\frac{\overline{A}}{2\Delta t}-\frac{\frac{{\left(\Delta y\right)}^2}{6}{C}^{-1}\overline{A}}{\Delta y\Delta t}-\frac{A_2}{\Delta y\Delta t}-\frac{C^{-1}\overline{A}{C}^{-1}D}{2\Delta t}\ge 0$.

Proof. 

The compact scheme for Equation (21) is given as

$$\begin{array}{l}\overline{A}{\delta }_t{\mathit{U}}_i^n+B\left({\delta }_y{\mathit{U}}_i^n-\frac{{\left(\mathsf{\Delta }y\right)}^2}{6}\left(C^{-1}\overline{A}{\mathit{U}}_{ty}+C^{-1}B{\mathit{U}}_{yy}-C^{-1}D{\mathit{U}}_y\right)\right)\\ \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}=C\left({\delta }_y^2{\mathit{U}}_I^n-\frac{{\left(\mathsf{\Delta }y\right)}^2}{12}\left(A_1{\mathit{U}}_{tt}+A_2{\mathit{U}}_{ty}+A_3{\mathit{U}}_{yy}-C^{-1}\overline{A}{C}^{-1}D{\mathit{U}}_t-B{C}^{-1}D{\mathit{U}}_y-D{\mathit{U}}_{yy}\right)\right)\end{array}$$

(25)

Let the difference between exact and numerical approximate solution is expressedas ${\mathit{e}}_i^n=\left|{\mathit{U}}_I^n-{\overline{\mathit{U}}}_i^n\right|$, where ${\overline{\mathit{U}}}_i^n$ shows an exact solution at the $i^{th}$ grid point and at the $n^{th}$ level. The error equation corresponding to Equation (25) is given as

$$\begin{array}{ll}\left(\frac{{\mathit{e}}_i^{n+1}-{\mathit{e}}_i^{n-1}}{2\mathsf{\Delta }t}\right)& +B(\frac{{\mathit{e}}_{i+1}^n-{\mathit{e}}_{i-1}^n}{2\mathsf{\Delta }y}\\ & -\frac{{\left(\mathsf{\Delta }y\right)}^2}{6}(C^{-1}\overline{A}\left\{\frac{{\mathit{e}}_{i+1}^{n+1}-{\mathit{e}}_{i-1}^{n+1}}{2\mathsf{\Delta }y\mathsf{\Delta }t}-\frac{{\mathit{e}}_{i+1}^{n-1}-{\mathit{e}}_{i-1}^{n-1}}{2\mathsf{\Delta }y\mathsf{\Delta }t}\right\}+C^{-1}B\frac{{\mathit{e}}_{i+1}^n-2{\mathit{e}}_i^n+e_{i-1}^n}{{\left(\mathsf{\Delta }y\right)}^2}\\ & -C^{-1}D\frac{{\mathit{e}}_{i+1}^n-{\mathit{e}}_{i-1}^n}{2\mathsf{\Delta }y}\left)\right)\\ & =C(\frac{{\mathit{e}}_{i+1}^n-2{\mathit{e}}_i^n+{\mathit{e}}_{i-1}^n}{{\left(\mathsf{\Delta }y\right)}^2}\\ & -\frac{{\left(\mathsf{\Delta }y\right)}^2}{12}(A_1\frac{{\mathit{e}}_i^{n+1}-2{\mathit{e}}_i^n+{\mathit{e}}_i^{n-1}}{{\left(\mathsf{\Delta }t\right)}^2}+A_2\left(\frac{{\mathit{e}}_{i+1}^{n+1}-{\mathit{e}}_{i-1}^{n+1}}{2\mathsf{\Delta }y\mathsf{\Delta }t}-\frac{{\mathit{e}}_{i+1}^{n-1}-{\mathit{e}}_{i-1}^{n-1}}{2\mathsf{\Delta }y\mathsf{\Delta }t}\right)+A_3\frac{{\mathit{e}}_{i+1}^n-2{\mathit{e}}_i^n+{\mathit{e}}_{i-1}^n}{{\left(\mathsf{\Delta }y\right)}^2}\\ & -C^{-1}\overline{A}{C}^{-1}D\frac{{\mathit{e}}_i^{n+1}-{\mathit{e}}_i^{n-1}}{2\mathsf{\Delta }t}-B{C}^{-1}D\frac{{\mathit{e}}_{i+1}^n-{\mathit{e}}_{i-1}^n}{2\mathsf{\Delta }y}-D\frac{{\mathit{e}}_{i+1}^n-2{\mathit{e}}_i^n+e_{i-1}^n}{{\left(\mathsf{\Delta }y\right)}^2}\left)\right)\end{array}$$

(26)

Applying norm on both sides of the Equation (26) gives

$$\begin{gathered} \frac{\Vert \overline{A}\Vert }{2\mathsf{\Delta }t}{e}^{n+1}\le \frac{\Vert \overline{A}\Vert e^{n-1}}{2\mathsf{\Delta }t}+\frac{\Vert B\Vert e^n}{2\mathsf{\Delta }t}+\frac{{\left(\mathsf{\Delta }y\right)}^2}{6}\left(\begin{array}{c}\frac{\Vert C^{-1}\overline{A}\Vert e^{n+1}}{\mathsf{\Delta }y\mathsf{\Delta }t}+\frac{e^{n-1}}{\mathsf{\Delta }y\mathsf{\Delta }t}\\ +\frac{\Vert C^{-1}B\Vert 4e^n}{{\left(\mathsf{\Delta }y\right)}^2}+\frac{\Vert C^{-1}D\Vert e^n}{\mathsf{\Delta }y}\end{array}\right) \\ +\Vert C\Vert \left(\frac{4e^n}{{\left(\mathsf{\Delta }y\right)}^2}+\frac{{\left(\mathsf{\Delta }y\right)}^2}{12}\left(\begin{array}{c}\Vert A_1\Vert \frac{4e^n}{{\left(\mathsf{\Delta }y\right)}^2}+\\ \Vert A_2\Vert \left(\frac{e^{n+1}}{\mathsf{\Delta }y\mathsf{\Delta }t}+\frac{e^{n-1}}{\mathsf{\Delta }y\mathsf{\Delta }t}\right)+\Vert A_3\Vert \frac{4e^n}{{\left(\mathsf{\Delta }y\right)}^2}+\\ \Vert C^{-1}\overline{A}{C}^{-1}D\Vert \frac{e^{n+1}+e^{n-1}}{2\mathsf{\Delta }t}\\ +\frac{\Vert B{C}^{-1}D\Vert e^n}{\mathsf{\Delta }y}+\frac{\Vert D\Vert 4e^n}{{\left(\mathsf{\Delta }y\right)}^2}\end{array}\right)\right) \end{gathered}$$

(27)

Equation (27) can be expressed as

$$\delta e^{n+1}\le {\delta }_1{e}^{n-1}+{\delta }_2{e}^n+M\left(O\left({\left(\mathsf{\Delta }t\right)}^2,{\left(\mathsf{\Delta }y\right)}^4\right)\right)$$

(28)

where $\delta =\frac{\Vert A\Vert }{2\mathsf{\Delta }t}-\frac{{\left(\mathsf{\Delta }y\right)}^2}{6}\Vert C\Vert \left(\frac{\Vert C^{-1}A\Vert }{\mathsf{\Delta }y\mathsf{\Delta }t}-\frac{\Vert A_2\Vert }{2\mathsf{\Delta }y\mathsf{\Delta }t}-\frac{\Vert C^{-1}A{C}^{-1}D\Vert }{4\mathsf{\Delta }t}\right)$

$$\begin{gathered} {\delta }_1=\Vert A\Vert \frac{1}{2\mathsf{\Delta }t}+\frac{\mathsf{\Delta }y}{6\mathsf{\Delta }t}+\frac{{\left(\mathsf{\Delta }y\right)}^2}{12}\Vert C\Vert \left(\frac{\Vert A_2\Vert }{\mathsf{\Delta }y\mathsf{\Delta }t}+\Vert C^{-1}\overline{A}{C}^{-1}D\Vert \frac{1}{2\mathsf{\Delta }t}\right) \\ {\delta }_2=\frac{\Vert B\Vert }{2\mathsf{\Delta }x}+\frac{{\left(\mathsf{\Delta }y\right)}^2}{6}\left(\Vert C^{-1}B\Vert \frac{4}{{\left(\mathsf{\Delta }y\right)}^2}+\Vert C^{-1}D\Vert \frac{1}{\mathsf{\Delta }y}\right)+\Vert C\Vert \left(\frac{4}{{\left(\mathsf{\Delta }y\right)}^2}+\frac{{\left(\mathsf{\Delta }y\right)}^2}{12}\left(\begin{array}{c}\Vert A_1\Vert \frac{4}{{\left(\mathsf{\Delta }y\right)}^2}+\Vert A_3\Vert \frac{4}{{\left(\mathsf{\Delta }y\right)}^2}+\\ \Vert B{C}^{-1}D\Vert \frac{1}{\mathsf{\Delta }y}+\Vert D\Vert \frac{4}{{\left(\mathsf{\Delta }y\right)}^2}\end{array}\right)\right) \end{gathered}$$

Inequality (28) can be expressed as

$$e^{n+1}\le {\delta }_3{e}^{n-1}+{\delta }_4{e}^n+M_1\left(O\left({\left(\mathsf{\Delta }t\right)}^2,{\left(\mathsf{\Delta }y\right)}^4\right)\right)$$

(29)

Let $e^{n-1}=\max \left(e^{n-1},e^n\right)$, then Inequality (28) can be expressed as:

$$e^{n+1}\le {\delta }_5{e}^{n-1}+M_1\left(O\left({\left(\mathsf{\Delta }t\right)}^2,{\left(\mathsf{\Delta }y\right)}^4\right)\right)$$

(30)

Letting $n=1$ in (30) gives

$$e^2\le {\delta }_5{e}^0+M_1$$

(31)

Since the initial condition is exact, $e^0=0$; therefore, (31) becomes

$$e^2\le M_1$$

(32)

Letting $n=2$ in (30) gives

$$e^3\le {\delta }_5{e}^1+M_1$$

(33)

Since the error at first time level is bounded, $e^1\le B$; therefore, Equation (33) becomes

$$e^3\le {\delta }_{5}B+M_1$$

(34)

Putting $n=3$ in Equation (30) results in

$$e^4\le {\delta }_5{e}^2+M_1\le \left({\delta }_5+1\right)M_1$$

(35)

Putting $n=4$ in Equation (30) gives

$$e^5\le {\delta }_5{e}^3+M_1\le {\delta }_5^{2}B+\left({\delta }_5+1\right)M_1$$

(36)

If it is continued, then for

$$e^{2n+1}\le {\delta }_5^{n}B+\left({\delta }_5^{n-1}+\dots +{\delta }_5+1\right)M_1.$$

(37)

and

$$e^{2n+2}\le \left({\delta }_5^{n-1}+\dots +{\delta }_5+1\right)M_1$$

(38)

This implies

$$e^{2n+1}\le {\delta }_2^{n}B+\left(\frac{{\delta }_5-1}{{\delta }_5-1}\right)M_1$$

(39)

and

$$e^{2n+2}\le \left(\frac{{\delta }_5-1}{{\delta }_5-1}\right)M_1$$

(40)

Since $\dots +{\delta }_5^{n-1}+\dots +{\delta }_5+1$ is an infinite geometric series, it converges; if $\left|{\delta }_5\right|<1$ and it gives convergence conditions, then this is the case when $\max \left(e^n,e^{n-1}\right)=e^{n-1}$. Similarly, for the case when $\max \left(e^n,e^{n-1}\right)=e^n$, the convergence conditions can be found.

The algorithm is constructed that can be used to summarize the use of code for finding the values of unknowns on each grid point and for every time level.

3. Applications

For finding numerical solutions of convection–diffusion problems, two examples are considered. The next part of this contribution comprises the application of the proposed explicit compact scheme for nonlinear PDEs and non-dimensional coupled PDEs for heat transfer in boundary layer flow. An Algorithm 1 is given which summarizes the computer programming for the solution of considered problems.

Algorithm 1. In the beginning, the smallest and largest values of spatial and temporal variables have been defined

Step 1. $\Delta x=\left(x_1-x_0\right)/\left(N_x-1\right)$ and $\Delta t=\left(t_1-t_0\right)/\left(N_t-1\right)$ are the step sizes used in space and time, where $N_x$ and $N_t$ denote grid points and time levels in space and time, respectively, and give values of the parameters contained in the equation(s) considered. Step 2. Use the forward Euler method on the first time level because the proposed scheme is constructed on three time levels. Step 3. Use the proposed fourth order in space and second order in time scheme on the remaining time levels to find one unknown explicitly. Step 4. When values of the dependent variable on each grid point and each time level are computed, a solution can be seen using different plots.

Example 1: Consider the following nonlinear convection–diffusion problem

$$\frac{\partial u}{\partial t}+u\frac{\partial u}{\partial x}=\frac{{\partial }^{2}u}{\partial x^2}$$

(41)

subject to the initial and boundary conditions

$$u\left(x,0\right)=1-\frac{2}{x},x>0$$

(42)

$$u\left(x_0,t\right)=1-\frac{2}{x_0-t}\mathrm{and}u\left(x_1,t\right)=1-\frac{2}{x_1-t}$$

(43)

The exact solution of Equation (41) with initial and boundary conditions (42) and (43), respectively, is given by

$$u\left(x,t\right)=1-\frac{2}{x-t}$$

(44)

The compact discretization for Equation (41) can be expressed as

$$u_i^{n+1}=\left[\begin{array}{c}{u}_i^{n-1}+\frac{2\Delta t}{{\left(\Delta x\right)}^2}\left\{u_{i+1}^n-2u_i^n+u_{i-1}^n\right\}-\\ 2\Delta t\left\{\begin{array}{c}\frac{{\left(\Delta x\right)}^2}{12}\left(\frac{u_i^{n+1}-2u_i^n+u_i^{n-1}}{{\left(\Delta t\right)}^2}\right)+\\ {\left(u_i^n\right)}^2\left(\frac{u_{i+1}^n-2u_i^n+u_{i-1}^n}{{\left(\Delta x\right)}^2}\right)\end{array}\right\}+\\ 2\left(\Delta t\right)u_i^n\left\{\begin{array}{c}\frac{u_{i+1}^n-u_{i-1}^n}{2\Delta x}+\\ \frac{u_i^n}{6}\left(u_{i+1}^n-2u_i^n+u_{i-1}^n\right)\end{array}\right\}\end{array}\right]$$

(45)

The compact numerical scheme (45) can be changed into an explicit scheme by collecting $u_i^{n+1}$ terms on one side. Thus, in this manner, the second order in time and the fourth order in space scheme for linearized Equation (41) are obtained. Note that the standard/classical second order in the time and space central scheme for Equation (41) becomes unconditionally unstable. Scheme (45), however, which is also second order in time and fourth order in space, becomes conditionally stable. Figure 1a,b show solutions obtained by the proposed scheme and comparisons of exact and standard/classical scheme’s solution with one obtained by the proposed scheme. The standard/classical scheme is first order in time and second order in space, and the relative error shown by this scheme is graphed in Figure 1b, which is greater than the relative error obtained by the proposed scheme. $N_x\mathrm{and}{N}_t$ show the number of grid points and time levels, while $x_0\mathrm{and}{x}_1$ show the left and right endpoints of the domain.

Example 2: Boundary Layer Flow Over Stretching Sheet

Consider an incompressible, laminar, Newtonian, and two-dimensional unsteady mixed convective fluid flow over an infinitely long plate. Let $x^{*}$-axis be along the plate, whereas the $y^{*}$-axis is perpendicular to the plate or perpendicular to the direction of the flow. The plate stretches towards the positive $x^{*}$-axis and the fluid is placed in the space $y^{*}>0$. Let the stretching velocity of the plate be $U$. Since the plate is infinitely long, the partial derivative in $x^{*}$ is neglected [23]. Then the governing equations of the flow can be expressed [23] as

$$\frac{\partial v^{*}}{\partial y^{*}}=0$$

(46)

$$\frac{\partial {\mathrm{u}}^{*}}{\partial t^{*}}+v^{*}\frac{\partial u^{*}}{\partial y^{*}}=\nu \frac{{\partial }^2{u}^{*}}{\partial y^{*}{}^2}+g{\beta }_T\left(T-T_{\infty }\right)$$

(47)

$$\frac{\partial \mathrm{T}}{\partial t^{*}}+v^{*}\frac{\partial T}{\partial y^{*}}=\alpha \frac{{\partial }^{2}T}{\partial y^{*}{}^2}+\frac{1}{\rho c_p}\frac{\partial q_r}{\partial y^{*}}$$

(48)

subject to the boundary conditions given as

$$\begin{array}{c}{u}^{*}\left(y^{*},t^{*}\right)=U, T\left(y^{*},t^{*}\right)=T_w when y^{*}\to 0\\ u^{*}\left(y^{*},t^{*}\right)=0, T\left(y^{*},t^{*}\right)=T_{\infty } when y^{*}\to \infty \end{array}\}$$

(49)

where $q_r$ is the heat flux and $u^{*}\mathrm{and}{v}^{*}$ represent the components of velocity in $x^{*}\mathrm{and}{y}^{*}$ directions, respectively; $\nu$ represents the kinematic viscosity, $g$ is used for gravity, ${\beta }_T$ represents the coefficient of volume expansion for temperature, $T$ is the fluid temperature, $\alpha$ represents thermal diffusivity, $T_w$ represents the temperature of the fluid at the wall, and $T_{\infty }$ represents the temperature of the fluid away from the plate.

Equations (46)–(49) are dimensional, so to reduce these equations into dimensionless forms, the following transformations are given as

$$u=\frac{u^{*}}{U},V=\frac{v^{*}}{U},t=\frac{t^{*}{U}^2}{\nu },y=\frac{U{y}^{*}}{\nu },\theta =\frac{T-T_{\infty }}{T_w-T_{\infty }}$$

(50)

Substituting transformations (50) into Equations (46)–(49), the dimensionless equations are given as

$$\frac{\partial v}{\partial y}=0$$

(51)

$$\frac{\partial u}{\partial t}+v\frac{\partial u}{\partial y}=\frac{{\partial }^{2}u}{\partial y^2}+\frac{G{r}_L}{R{e}^3}\theta$$

(52)

$$\frac{\partial \theta }{\partial t}+v\frac{\partial \theta }{\partial y}=\frac{1}{P_r}\left(1+\frac{4}{3}{R}_d\right)\frac{{\partial }^2\theta }{\partial y^2}.$$

(53)

subject to the dimensionless boundary conditions

$$\begin{array}{c}u\left(y,t\right)=1, \theta \left(y,t\right)=1 when y\to 0\\ u\left(y,t\right)=0, \theta \left(y,t\right)=0 when y\to \infty \end{array}\}$$

(54)

where $G{r}_L$ represents Grashoff number of heat transfer, $R_d$ is the radiation parameter, $R_e$ represents the Reynolds number, and $P_r$ represents the Prandtl number; these parameters are defined by

$$G{r}_L=\frac{g{\beta }_T\left(T_w-T_{\infty }\right)L^3}{{\nu }^2},R_d=\frac{4\sigma T_{\infty }^3}{k{k}^{*}},R_e=\frac{\nu }{UL} and P_r=\frac{\nu }{\alpha }$$

(55)

where $\sigma$ and $k\mathrm{and}{k}^{*}$ are, respectively, the Stefan–Boltzmann constant, the thermal conductivity of the fluid, and mean absorption coefficient, and parameters have appeared where radiation parameter $R_d$ was defined in (54). Additionally, in Equation (53), linearized Rosseland radiative flux is considered. The solution of Equation (51) is expressed by

$$v=v_0$$

(56)

where $v_0$ is suction/injection velocity of the plate.

For solving convection–diffusion type Equations (52)–(53) using the proposed scheme, a difference/discretized equation can be developed in the following manner,

$$u_i^{n+1}=\left[\begin{array}{c}{u}_i^{n-1}-\frac{V\Delta t}{\Delta y}\left(u_{i+1}^n-u_{i-1}^n\right)+\\ 2\Delta t\left(\frac{{\left(V\Delta y\right)}^2}{12}+1\right)\left(\frac{u_{i+1}^n-2u_i^n+u_{i-1}^n}{{\left(\Delta y\right)}^2}\right)+\\ \Delta t\left(-\frac{VN{\left(\Delta y\right)}^2}{12}\right)\left(\frac{{\theta }_{i+1}^n-{\theta }_{i-1}^n}{2\Delta y}\right)-\\ \Delta t\left(\frac{{\left(\Delta y\right)}^2}{6}\right)\left(\frac{u_i^{n+1}-2u_i^n+u_i^{n-1}}{{\left(\Delta t\right)}^2}\right)+\\ \Delta t\left(\frac{N{\left(\Delta y\right)}^2}{6}\right)\left(\frac{{\theta }_i^{n+1}-{\theta }_i^{n-1}}{2\Delta t}\right)\\ \Delta t\left(\frac{N}{6}\right)\left({\theta }_{i+1}^n-2{\theta }_i^n+{\theta }_{i-1}^n\right)+\\ 2\Delta tN{\theta }_i^n\end{array}\right]$$

(57)

$${\theta }_i^{n+1}=\left[\begin{array}{c}{\theta }_i^{n-1}-\frac{V\Delta t}{\Delta x}\left({\theta }_{i+1}^n-{\theta }_{i-1}^n\right)+\\ 2\Delta t\left(\frac{{\left(V\Delta x\right)}^2}{12a}+a\right)\left(\frac{{\theta }_{i+1}^n-2{\theta }_i^n+{\theta }_{i-1}^n}{{\left(\Delta x\right)}^2}\right)-\\ \Delta t\left(\frac{{\left(\Delta x\right)}^2}{6a}\right)\left(\frac{u_i^{n+1}-2u_i^n+u_i^{n-1}}{{\left(\Delta t\right)}^2}\right)\end{array}\right]$$

(58)

where $a=\frac{1}{P_r}\left(1+\frac{4R_d}{3}\right)$ and $N=\frac{G{r}_L}{R_e^3}$. Since the term contained in Eq. (58), having the form ${\theta }_i^{n+1}$, which is computed from Equation (57). Thus, two different procedures can be adopted. The one using some iteration scheme and this is an implicit kind of procedure to find the values of unknowns, or another one is an explicit procedure that finds ${\theta }_i^{n+1}$ first and afterward, the term $u_i^{n+1}$ is calculated. In this contribution, an explicit procedure is adopted. Some results are compared with existing schemes. Figure 2 shows the solution obtained by the proposed scheme with the exact solution and two existing schemes. The comparison with the exact solution is given in Figure 2a, while the comparison of absolute errors made by three schemes is drawn in Figure 2b. The proposed scheme produced less error than the first order in time and second order in space scheme, and the absolute error obtained by the DuFort–Frankel method oscillates. Figure 3 shows the comparison of the proposed scheme with the DuFort–Frankel scheme through time. Skin friction coefficient starts with zero due to an initial condition is set to zero, and then for the next time level, it starts with its highest value. Then it decreases and maintains an approximately/exactly constant behavior along the $x$-axis. The difference between skin friction coefficients obtained by both schemes can be seen in the zoomed figure. In addition, it can be observed that the skin friction coefficient oscillates, which is computed from the existing unconditionally stable DuFort–Frankel method. The difference between two local Nusselt numbers obtained by the proposed and existing schemes can be seen in Figure 4. The skin friction coefficient and local Nusselt number decrease over time due to the effects of applied boundary conditions on velocity and temperature profiles, and this effect decreases with time. Figure 5, Figure 6, Figure 7, Figure 8 and Figure 9 show contour plots for velocity profile for only the diffusion parabolic equation with diffusion effect, including the source term and convection–diffusion equations. The mesh plots, Figure 9 and Figure 10, show the surface type of plot that can be used to understand the contour plots. The symbol $x_0$ is used left end point of domain, $x_1$ is used for right end point of domain, $t_f$ is used for final time, $N_y$ is used for number of grid points and $N_t$ is used for number of time levels.

4. Conclusions

A fourth-order explicit compact numerical scheme was proposed for solving convection–diffusion equations. The stability for the scalar case and convergence for the system of partial differential equations were given. A particular case of the existing mathematical model for the heat transfer in mixed convection fluid flow over a stretching sheet with the influence of thermal radiations was given. The proposed compact numerical scheme was used to solve the dimensionless equations in fluid phenomena with heat transfer effects, and some results were compared with existing schemes. It is hoped that the results in the current report are a helpful source for future fluid-flow studies in an industrial environment in an enclosure area.