Approximation of 3D convection diffusion equation using DQM based on modified cubic trigonometric B-splines

Abstract

This paper presents an approximate solution of 3D convection diffusion equation (CDE) using DQM based on modified cubic trigonometric B-spline (CTB) basis functions. The DQM based on CTB basis functions is used to integrate the derivatives of space variables which transformed the CDE into the system of first order ODEs. The resultant system of ODEs is solved using SSPRK (5,4) method. The solutions are approximated numerically and also presented graphically. The accuracy and efficiency of the method are validated by comparing the solutions with existing numerical solutions. The stability analysis of the method is also carried out.

Keywords

3D convection diffusion equation CTB basis functions DQM SSPRK (5 4)stability analysis

1. Introduction

Since the partial differential equations (PDEs) play a key role in mathematical modeling of various real world problems which usually represent a physical, chemical and/or biological phenomenon. Therefore it becomes necessary to solve these equations as accurately as possible in order to predict the absolute behavior of these models. In general, the exact solution for every PDE is not available or is not an easy task to find the exact solutions, and so, the study of the numeric solutions of the PDEs is a very relevant area of research.

The 3D CDE has found many applications in various fields such as astrophysics, oceanography, meteorology, fluid motion, semiconductors, heat transfer, hydraulics, sediment and pollutant transport, and chemical engineering. Especially in fluid dynamics and heat transfer, the CDE plays an important role. The CDE expresses the convection-diffusion of quantities in fluid dynamics and computational hydraulics. We consider 3D CDE as follows:

$\displaystyle\frac{\partial u}{\partial t}+{\beta}^{\prime}_{x}\frac{\partial u% }{\partial x}+{\beta}^{\prime}_{y}\frac{\partial u}{\partial y}+{\beta}^{% \prime}_{z}\frac{\partial u}{\partial z}={\alpha}^{\prime}_{x}\frac{\partial^{% 2}u}{\partial x^{2}}+{\alpha}^{\prime}_{y}\frac{\partial^{2}u}{\partial y^{2}}% +{\alpha}^{\prime}_{z}\frac{\partial{}^{2}u}{\partial z^{2}},({x,y,z})\in% \Omega,t\geqslant 0$ (1)

subject to initial condition

$\displaystyle u(x,y,z,0)=g(x,y,z),({x,y,z})\in\Omega,$

and boundary conditions

$\displaystyle\left.{\begin{array}[]{l}u(0,y,z,t)=f_{1}(y,z,t),u(1,y,z,t)=f_{2}% (y,z,t),\\ u(x,0,z,t)=f_{3}(x,z,t),u(x,1,z,t)=f_{4}(x,z,t),\\ u(x,y,0,t)=f_{5}(x,y,t),u(x,y,1,t)=f_{6}(x,y,t),\\ \end{array}}\right\}t\geqslant 0$

where $\Omega=\{{({x,y,z}):0\leqslant x\leqslant 1,0\leqslant y\leqslant 1,0\leqslant z% \leqslant 1}\}$ is is the computational domain and $u(x,y,\linebreak z,t)$ is a function which represents vorticity, heat, etc., the constant coefficients ${\alpha}^{\prime}_{x},{\alpha}^{\prime}_{y},{\alpha}^{\prime}_{z}>0$ and ${\beta}^{\prime}_{x},{\beta}^{\prime}_{y},{\beta}^{\prime}_{z}>0$ denote the diffusion coefficients and convective velocities respectively.

In the few decades, various numerical schemes have been proposed for solving the convection-diffusion equations. For example, Dawson et al. [1] developed a Lattice Boltzmann method for one dimensional convection-diffusion problem while Ding and Zhang [2] used a semi-discrete and a padé approximation method for proposing a new difference scheme. Tian and Yua [3] developed a HOE scheme for aforesaid problem. He [4] used Taylor series method and Pade approximation to approximate 1D CDE for E reaction arising in rotating disk electrodes. Dhawan et al. [5] used method based on linear and quadratic B-spline functions while authors of [6, 7] presented a fourth-order compact FDM and a new compact unconditionally stable FDM to approximate the one-dimensional advection diffusion equation. Singh et al. [8] used FDM to approximate this equation having nonlinear source/sink term.

Tian and Ge [9] applied an exponential high-order compact ADI method while Dehghan and Mohebbi [10] proposed new class of high-order accurate unconditionally stable methods of order $O({h^{4},t^{4}})$ for approximating two-dimensional CDE. Karaa and Zhang [11] proposed a high order ADI method for solving aforesaid problem while Kalita et al. [12] developed a class of HOC schemes to solve this problem with variable convection coefficients. Tian [13] developed an unconditionally stable rational HOC ADI method of order $O({h^{4},t^{2}})$ , to approximate 2D convection-diffusion problems. A family of unconditionally stable finite difference methods of order $O({h^{3},t^{2}})$ was presented by Cecchi and Pirozzi [14]. Dehghan [15] developed two-level fully explicit and implicit FDM for approximating 3D advection-diffusion equation. Mohanty and Singh [16] employed a fourthorder CFD algorithm to solve the 3D nonlinear singularly perturbed elliptic PDEs. Shukla and Tamsir [17] presented an Expo-MCB-DQM to approximate 2D and 3D convection-diffusion equations. Prieto et al. [18] demonstrated the utility of generalized explicit FDM in order to solve the 3D advection-diffusion equation. Ravnik and Tibuat [19] developed a boundary-domain integral formulation of the CDE.

The DQM was first proposed by Bellman et al. [20]. The various test functions such as spline functions, Lagrange interpolation polynomials and sinc functions [21, 22, 23, 24, 25] were used to calculate weighting coefficients in DQM. Tamsir et al. [27] proposed a DQM based on trigonometric B-spline basis functions to approximate Fisher’s reaction-diffusion equations. Mittal et al. [28] used cubic B-spline quasi-interpolation and differential quadrature methods for approximating the fourth-order parabolic PDEs.

2. Description of the method

Let the computational domain $\Omega$ is uniformly partitioned by taking grid points, $P_{x}$ , $P_{y}$ and $P_{z}$ along with the step sizes $h_{x}=1/{(P_{x}-1)}$ , $h_{y}=1/{(P_{y}-1)}$ , and $h_{z}=1/{(P_{z}-1)}$ . It is assumed that the solution $u({x,y,z,t})$ of the Eq. (1) at the grid $({x_{i},y_{j},z_{k},t})$ is $u_{\textit{ijk}}=u({x_{i},y_{j},z_{k},t})$ , where $i=1,2,\ldots,P_{x}$ , $j=1,2,\ldots,P_{y}$ and $k=1,2,\ldots,P_{z}$ . Partial derivatives of $u({x,y,z,t})$ w. r .t $x$ , $y$ and $z$ are

$\displaystyle\left\{{\begin{array}[]{l}\frac{\partial^{r}u_{\textit{ijk}}}{% \partial x^{r}}=\sum\limits_{p=1}^{P_{x}}{l_{{ip}}^{(r)}u_{\textit{pjk}}},i=1,% 2,\ldots,P_{x}\\ \frac{\partial^{r}u_{\textit{ijk}}}{\partial y^{r}}=\sum\limits_{p=1}^{P_{y}}{% m_{{jp}}^{(r)}u_{\textit{ipk}}},j=1,2,\ldots,P_{y,}\\ \frac{\partial^{r}u_{\textit{ijk}}}{\partial z^{r}}=\sum\limits_{p=1}^{P_{z}}{% n_{{kp}}^{(r)}u_{\textit{ijp}}},k=1,2,\ldots,P_{z,}\\ \end{array}}\right.$ (2)

where $l_{{ip}}^{(r)}$ , $m_{{jp}}^{(r)}$ and $n_{{kp}}^{(r)}$ , $r=1,2$ are the weighting coefficients.

The CTB basis functions [29, 30, 31, 32, 33] at the knots are given as:

$\displaystyle S_{i}(x)=\frac{1}{\zeta}\left\{{{\begin{array}[]{ll}{\xi^{3}({x_% {i-2}}),}&{x\in[{x_{i-2},x_{i-1}})}\\ {\xi({x_{i-2}})({\xi({x_{i-2}})\psi({x_{i}})+\xi({x_{i-1}})\psi({x_{i+1}})})+% \xi^{2}({x_{i-1}})\psi({x_{i+2}}),}&{x\in[{x_{i-1},x_{i}})}\\ {\xi({x_{i-2}})\psi^{2}({x_{i+1}})+\psi({x_{i+2}})({\xi({x_{i-1}})\psi({x_{i+1% }})+\xi({x_{i}})\psi({x_{i+2}})}),}&{x\in[{x_{i},x_{i+1}})}\\ {\psi^{3}({x_{i+2}}),}&{x\in[{x_{i+1},x_{i+2}})}\\ {0,}&{\text{otherwise}}\\ \end{array}}}\right.$ (3)

where $\xi=\sin\left({\frac{x-x_{i}}{2}}\right),\psi=\sin\left({\frac{x_{i}-x}{2}}% \right),\zeta=\sin({0.5h})\sin(h)\sin({1.5h})$ .

The values of $S_{i}(x),S^{\prime}_{i}(x)$ and ${S}^{\prime\prime}_{i}(x)$ at the knots are given by

$\displaystyle S_{i}({x_{j}})=\left\{{{\begin{array}[]{ll}a&{\text{if }i=j,}\\ b&{\text{if }i=j\pm 1,}\\ 0&{\text{else,}}\\ \end{array}}}\right.;{S}^{\prime}_{i}({x_{j}})=\left\{{{\begin{array}[]{ll}c&{% \text{if }i=j+1,}\\ {-c}&{\text{if }i=j-1,}\\ 0&{\text{else,}}\\ \end{array}}}\right.;{S}^{\prime\prime}_{i}({x_{j}})=\left\{{{\begin{array}[]{% ll}e&{\text{if }i=j,}\\ d&{\text{if }i=j\pm 1,}\\ 0&{\text{else,}}\\ \end{array}}}\right.$

where $a=\sin^{2}(0.5h)\text{csc}(h)\text{csc}(1.5h)$ , $b=2\cdot({1+2\cos(h)})^{-1}$ , $c=0.75\text{csc}(1.5h)$ , $d=3(3\cos(h)\linebreak+1)\csc^{2}({0.5h})\cdot({16+({2\cos({0.5h})+\cos({1.5h}% )})})^{-1}$ , $e=-3\text{cot}^{2}(1.5h)\cdot(2+4\text{cos}(h))^{-1}$ .

The CTB basis functions are modified as follows [34]:

$\displaystyle\left.{\begin{array}[]{l}\hat{S}_{1}(x)=S_{1}(x)+2S_{0}(x),\\ \hat{S}_{2}(x)=S_{2}(x)-S_{0}(x),\\ \hat{S}_{m}(x)=S_{m}(x),\text{for }m=3,\ldots,P_{x}-2,\\ \hat{S}_{N-1}(x)=S_{P_{x}-1}(x)-S_{P_{x}+1}(x),\\ \hat{S}_{P_{x}}(x)=S_{P_{x}}(x)+2S_{P_{x}+1}(x),\\ \end{array}}\right\}$ (4)

where $\{\hat{S}_{1},\hat{S}_{2},\ldots,\hat{S}_{P_{x}}\}$ forms a basis over the computational domain.

2.1 Computation of the weighting coefficients

Taking $r=1$ and putting the values of $\hat{S}_{k}(x)({k=1,2,\ldots,P_{x}})$ in Eq. (2), we get

$\displaystyle{\hat{S}}^{\prime}_{k}({x_{i}})=\sum\limits_{j=1}^{P_{x}}{l_{ij}^% {(1)}}\hat{S}_{k}({x_{j}}),\forall i,k=1,2,\ldots,P_{x}.$ (5)

Using Eqs (3) and (4) in above equation, we get a system of linear equations

$\displaystyle A\vec{l}^{(1)}[i]=\vec{R}[i],\text{for }i=1,2,\ldots,P_{x},$ (6)

where tri-diagonal matrix $A$ takes the form of

$\displaystyle A(1,1)=2a+b,A(1,2)=b,$ $\displaystyle A(2,1)=0,A(2,2)=b,A(2,3)=a,$ $\displaystyle A(i-1,i)=A(i+1,i)=a,A(i,i)=b,\text{ for }i=3,4,\ldots,P_{x}-2,$ $\displaystyle A(P_{x}-1,P_{x})=0,A(P_{x}-1,P_{x}-1)=b,A(P_{x}-1,P_{x}-2)=a,$ $\displaystyle A(P_{x},P_{x}-1)=a;A(P_{x},P_{x})=2a+b.$

The $\vec{l}^{(1)}[i]=[{l_{{i1}}^{(1)},l_{i2}^{(1)},\ldots,l_{{iP_{x}}}^{(1)}}]^{T}$ are the first order weighting coefficients at the knot points $x_{i}$ . The vectors $\vec{R}[i]=[{{\hat{S}}^{\prime}_{1,i},{\hat{S}}^{\prime}_{2,i},\ldots,{\hat{S}% }^{\prime}_{P_{x},i}}]^{tr}$ at $x_{i}(i=1,2,\ldots,P_{x})$ are given by

$\displaystyle\vec{R}[1]=\left[{{\begin{array}[]{*{20}c}{-2c}&{2c}&0&0&{0\cdots% }&0\\ \end{array}}}\right]^{tr},$ $\displaystyle\vec{R}[2]=\left[{{\begin{array}[]{*{20}c}{-c}&0&c&0&{0\cdots}&0% \\ \end{array}}}\right]^{tr},$ $\displaystyle\vdots$ $\displaystyle\vec{R}[{P_{x}-1}]=\left[{{\begin{array}[]{*{20}c}0&0&\cdots&0&{-% c}&0&c\\ \end{array}}}\right]^{tr},$ $\displaystyle\vec{R}[{P_{x}}]=\left[{{\begin{array}[]{*{20}c}0&0&\cdots&0&0&{-% 2c}&{2c}\\ \end{array}}}\right]^{tr}.$

Thomas algorithm is used to solve above system, given by Eq. (6); consequently, we have $l_{ip}^{(1)}$ . The weighting coefficients $l_{ij}^{(2)},1\leqslant i,j\leqslant P_{x}$ are computed using the formula [25] as:

$\displaystyle\left\{{\begin{array}[]{l}l_{ij}^{(r)}=r\left({l_{ij}^{(1)}l_{ii}% ^{({r-1})}-\frac{l_{ij}^{({r-1})}}{x_{i}-x_{j}}}\right),\text{ for }i\neq j% \text{ and }r=2,3,\ldots,P_{x}-1;i=1,2,\ldots,P_{x}\\ l_{ii}^{(r)}=-\sum\limits_{j=1,j\neq i}^{P_{x}}{l_{ij}^{(r)},\text{ for }i=j,}% \\ \end{array}}\right.$ (7)

The weighting coefficients $m_{ij}^{(1)}$ and $n_{ij}^{(1)}$ can be obtained in similar way.

3. Discretization of the problem

Using the approximated values of spatial derivatives in Eq. (1) we get

$\displaystyle\left\{{\begin{array}[]{l}\frac{du_{\textit{ijk}}}{dt}={\alpha}^{% \prime}_{x}\sum\limits_{p=1}^{P_{x}}{l_{{ip}}^{(2)}}u_{\textit{pjk}}+{\alpha}^% {\prime}_{y}\sum\limits_{p=1}^{P_{y}}{m_{{jp}}^{(2)}}u_{\textit{ipk}}+{\alpha}% ^{\prime}_{z}\sum\limits_{p=1}^{P_{z}}{n_{{kp}}^{(2)}}u_{\textit{ijp}}-{\beta}% ^{\prime}_{x}\sum\limits_{p=1}^{P_{x}}{l_{{ip}}^{(1)}}u_{\textit{pjk}}\\ \quad-{\beta}^{\prime}_{y}\sum\limits_{p=1}^{P_{y}}{m_{{jp}}^{(1)}}u_{\textit{% ipk}}-{\beta}^{\prime}_{z}\sum\limits_{p=1}^{P_{z}}{n_{{kp}}^{(1)}}u_{\textit{% ijp}},\\ u_{\textit{ijk}}(t=0)=g(x_{i},y_{j},z_{k}),i=1,2,\ldots,P_{x},j=1,2,\ldots,P_{% y},k=1,2,\ldots,P_{z}.\\ \end{array}}\right.$ (8)

At the boundary, the solution can be rewritten as

$\displaystyle\left\{{\begin{array}[]{l}u_{1jk}=f_{1}({y_{j},z_{k},t}),u_{P_{x}% jk}=f_{2}({y_{j},z_{k},t}),\\ u_{i1k}=f_{3}({x_{i},z_{k},t}),u_{iP_{y}k}=f_{4}({x_{i},z_{k},t}),\\ u_{ij1}=f_{5}({x_{i},y_{j},t}),u_{ijP_{z}}=f_{6}({x_{i},y_{j},t}),\\ \end{array}}\right.({x_{i},y_{j},z_{k},t})\in\partial\Omega\text{ and }t>0.$

Using boundary values, Eq. (8) can be written as

$\displaystyle\left\{{\begin{array}[]{l}\frac{du_{\textit{ijk}}}{dt}={\alpha}^{% \prime}_{x}\sum\limits_{p=2}^{P_{x}-1}{l_{{ip}}^{(2)}}u_{\textit{pjk}}+{\alpha% }^{\prime}_{y}\sum\limits_{p=2}^{P_{y}-1}{m_{{jp}}^{(2)}}u_{\textit{ipk}}+{% \alpha}^{\prime}_{z}\sum\limits_{p=2}^{P_{z}-1}{n_{{kp}}^{(2)}}u_{\textit{ijp}% }-{\beta}^{\prime}_{x}\sum\limits_{p=2}^{P_{x}-1}{l_{{ip}}^{(1)}}u_{\textit{% pjk}}\\ \quad-{\beta}^{\prime}_{y}\sum\limits_{p=2}^{P_{y}-1}{m_{{jp}}^{(1)}}u_{% \textit{ipk}}-{\beta}^{\prime}_{z}\sum\limits_{p=2}^{P_{z}-1}{n_{{kp}}^{(1)}}u% _{\textit{ijp}}+F_{\textit{ijk}},\\ u_{\textit{ijk}}(t=0)=g(x_{i},y_{j},z_{k}),i=2,3,\ldots,P_{x}-1,j=2,3,\ldots,P% _{y}-1,k=2,3,\ldots,P_{z}-1.\\ \end{array}}\right.$ (9)

where

$\displaystyle F_{\textit{ijk}}={\alpha}^{\prime}_{x}({l_{i1}^{(2)}u_{1jk}+l_{% iP_{x}}^{(2)}u_{P_{x}jk}})+{\alpha}^{\prime}_{y}({m_{j1}^{(2)}u_{i1k}+m_{jP_{y% }}^{(2)}u_{iP_{y}k}})+{\alpha}^{\prime}_{z}({n_{k1}^{(2)}u_{ij1}+n_{kP_{z}}^{(% 2)}u_{\textit{ijP}_{z}}}){}-{\beta}^{\prime}_{x}({l_{i1}^{(1)}u_{1jk}+l_{iP_{x% }}^{(1)}u_{P_{x}jk}})-{\beta}^{\prime}_{y}({m_{j1}^{(1)}u_{i1k}+m_{jP_{y}}^{(1% )}u_{iP_{y}k}})-{\beta}^{\prime}_{z}({n_{k1}^{(1)}u_{ij1}+n_{kP_{z}}^{(1)}u_{% \textit{ijP}_{z}}}).$

Finally, SSP-RK54 scheme [26] is employed to solve Eq. (9).

4. Stability analysis

The Eq. (9) can be expressed in terms of matrix form as:

$\displaystyle\frac{dU}{dt}=BU+\vec{F},$ (10)

where $U=({u_{222},u_{223},\ldots,u_{22({P_{z}-1})},u_{322},u_{323},\ldots,u_{32({P_{% z}-1})}\ldots,u_{({P_{x}-1})({P_{y}-1})({P_{z}-1})}})^{tr}$ is a solution vector, $\vec{F}=[{F_{\textit{ijk}}}]$ contains the boundary values.

$B={\alpha}^{\prime}_{x}A_{x}^{2}+{\alpha}^{\prime}_{y}B_{y}^{2}+{\alpha}^{% \prime}_{z}C_{z}^{2}-{\beta}^{\prime}_{x}A_{x}^{1}-{\beta}^{\prime}_{y}B_{y}^{% 1}-{\beta}^{\prime}_{z}C_{z}^{1}$ is a square matrix of order $({P_{x}-2})\times({P_{y}-2})\times({P_{z}-2})$ where $A_{x}^{r},B_{y}^{r}$ and $C_{z}^{r}({r=1,2})$ are as given below

$\displaystyle A_{x}^{r}=\left[{{\begin{array}[]{*{20}c}{l_{{22}}^{(r)}I_{x}}&{% l_{{23}}^{(r)}I_{x}}&\ldots&{l_{{2({P_{x}-1})}}^{(r)}I_{x}}\\ {l_{{32}}^{(r)}I_{x}}&{l_{{33}}^{(r)}I_{x}}&\ldots&{l_{{3({P_{x}-1})}}^{(r)}I_% {x}}\\ \vdots&\vdots&\ddots&\vdots\\ {l_{{({P_{x}-1})2}}^{(r)}I_{x}}&{l_{{({P_{x}-1})3}}^{(r)}I_{x}}&\ldots&{l_{{({% P_{x}-1})({P_{x}-1})}}^{(r)}I_{x}}\\ \end{array}}}\right];$ (11) $\displaystyle B_{y}^{r}=\left[{{\begin{array}[]{*{20}c}{M_{y}}&{O_{y}}&\ldots&% {O_{y}}\\ {l_{{2,1}}^{(2)}}&{M_{y}}&\ldots&{O_{y}}\\ \vdots&\vdots&\ddots&\vdots\\ {O_{y}}&{O_{y}}&\ldots&{M_{y}}\\ \end{array}}}\right],\mbox{ where }M_{y}=\left[{{\begin{array}[]{*{20}c}{m_{{2% 2}}^{(r)}I_{z}}&{m_{{23}}^{(r)}I_{z}}&\ldots&{m_{{2({P_{y}-1})}}^{(r)}I_{z}}\\ {m_{{32}}^{(r)}I_{z}}&{m_{{33}}^{(r)}I_{z}}&\ldots&{m_{{3({P_{y}-1})}}^{(r)}I_% {z}}\\ \vdots&\vdots&\ddots&\vdots\\ {m_{{({P_{y}-1})2}}^{(r)}I_{z}}&{m_{{({P_{y}-1})3}}^{(r)}I_{z}}&\ldots&{m_{{({% P_{y}-1})({P_{y}-1})}}^{(r)}I_{z}}\\ \end{array}}}\right];$ (12) $\displaystyle C_{z}^{r}=\left[{{\begin{array}[]{*{20}c}{M_{z}}&{O_{z}}&\ldots&% {O_{z}}\\ {l_{{2,1}}^{(2)}}&{M_{z}}&\ldots&{O_{z}}\\ \vdots&\vdots&\ddots&\vdots\\ {O_{z}}&{O_{z}}&\ldots&{M_{z}}\\ \end{array}}}\right],\mbox{ where }M_{z}=\left[{{\begin{array}[]{*{20}c}{n_{{2% 2}}^{(r)}}&{n_{{23}}^{(r)}}&\ldots&{n_{{2({P_{z}-1})}}^{(r)}}\\ {n_{{32}}^{(r)}}&{n_{{33}}^{(r)}}&\ldots&{n_{{3({P_{z}-1})}}^{(r)}}\\ \vdots&\vdots&\ddots&\vdots\\ {n_{{({P_{z}-1})2}}^{(r)}}&{n_{{({P_{z}-1})3}}^{(r)}}&\ldots&{n_{{({P_{z}-1})(% {P_{z}-1})}}^{(r)}}\\ \end{array}}}\right],$ (13)

where $O_{y}$ , $O_{z}$ are null matrices of order $({P_{y}-2})\times 1$ and $({P_{z}-2})\times 1$ , respectively. $I_{x}$ and $I_{y}$ are identity matrices of order $({P_{x}-2})\times 1$ and $({P_{y}-2})\times 1$ , respectively.

The stability of the method depends on the Eigen values $\lambda_{i}$ of the coefficient matrix $B$ . If the $Re({\lambda_{i}})\leqslant 0$ , the method will be stable. Moreover, for the complex Eigen values, some tolerance may exists i.e. the real part of Eigen values may be significantly small positive. Figure 1 evident that the real parts of the Eigen values of matrix $B$ , for different grid sizes, are all negative and hence the method is unconditional stability.

Figure 1.

Eigen values of the matrix $B$ corresponding to different grid sizes.

5. Numerical results and discussion

In this section, we consider two numerical examples in order to check the accuracy of the method.

Example 1: In this example, the numerical solution of 3D convection-diffusion equation with $\alpha^{\prime}_{x}=\alpha^{\prime}_{y}=\alpha^{\prime}_{z}=0.01$ and $\beta^{\prime}_{x}=\beta^{\prime}_{y}=\beta^{\prime}_{z}=0.8$ over the region [0, 1] $\times$ [0, 1] $\times$ [0, 1] with the exact solution $u({x,y,z})=\exp[-\{({x-0.8t-0.5})^{2}+({y-0.8t-0.5})^{2}+({x-0.8t-0.5})^{2}\}/% (0.01({1+4t})\}]\frac{1}{({1+4t})^{\frac{3}{2}}}$ .

This example is solved for parametric values: $\Delta t=$ 0.01, $h=$ 0.1, 0.05 and 0.025 (where $h_{x}=h_{y}=h_{z}=h$ ) at different time levels and comparative results are presented by Tables 1 and 2 for absolute and root mean error norms, respectively. The obtained findings are compared with those available in literature such as FTBSCS [14], Lax-Wendroff [15], expo-MCB-DQM [17] and also presented in Table 1 which shows that the present method produces highly accurate results than the methods presented in [15] and comparatively better than that those presented in [17]. Table 2 demonstrates the comparison of numerical results with those given by another method called expo-MCB-DQM [17] which shows that present results are comparatively better. Simultaneously Table 2 also depicts that the error norms are decreasing on increasing grid sizes. A graphical representation of absolute errors for $h=$ 0.05, $\Delta t=$ 0.001 at $t=$ 1 is given by Fig. 2. Figure 3 shows the RMS error norms with $\Delta t=$ 0.01 for different grid sizes. As we can see from this figure that the errors are decreasing on increasing the grid sizes.

Table 1
Comparative study for absolute error norms

$(x,y,z)$ such that $x=y=z$	Present method	FTBSCS method [15]	Lax-wendroff method [15]	Expo-MCB-DQM [17] $h=$ 0.1 $t=$ 1
				$\Delta t=$ 0.01	$\Delta t=$ 0.001
0.1	5.28378 $\times$ 10 ${}^{-7}$	4.6 $\times$ 10 ${}^{-3}$	5.4 $\times$ 10 ${}^{-4}$	5.1 $\times$ 10 ${}^{-7}$	7.9 $\times$ 10 ${}^{-7}$
0.2	2.19713 $\times$ 10 ${}^{-7}$	4.6 $\times$ 10 ${}^{-3}$	5.5 $\times$ 10 ${}^{-4}$	2.2 $\times$ 10 ${}^{-7}$	2.2 $\times$ 10 ${}^{-7}$
0.3	2.65127 $\times$ 10 ${}^{-6}$	4.4 $\times$ 10 ${}^{-3}$	5.3 $\times$ 10 ${}^{-4}$	2.6 $\times$ 10 ${}^{-7}$	3.0 $\times$ 10 ${}^{-7}$
0.4	1.85918 $\times$ 10 ${}^{-6}$	4.3 $\times$ 10 ${}^{-3}$	5.2 $\times$ 10 ${}^{-4}$	1.9 $\times$ 10 ${}^{-7}$	2.6 $\times$ 10 ${}^{-6}$
0.5	5.70496 $\times$ 10 ${}^{-5}$	4.2 $\times$ 10 ${}^{-3}$	5.1 $\times$ 10 ${}^{-4}$	5.7 $\times$ 10 ${}^{-7}$	6.9 $\times$ 10 ${}^{-5}$
0.6	1.4507 $\times$ 10 ${}^{-5}$	4.0 $\times$ 10 ${}^{-3}$	4.9 $\times$ 10 ${}^{-4}$	1.5 $\times$ 10 ${}^{-7}$	1.7 $\times$ 10 ${}^{-5}$
0.7	5.51423 $\times$ 10 ${}^{-5}$	4.1 $\times$ 10 ${}^{-3}$	5.0 $\times$ 10 ${}^{-4}$	5.5 $\times$ 10 ${}^{-7}$	6.6 $\times$ 10 ${}^{-5}$
0.8	7.59912 $\times$ 10 ${}^{-6}$	4.4 $\times$ 10 ${}^{-3}$	4.8 $\times$ 10 ${}^{-4}$	7.7 $\times$ 10 ${}^{-7}$	8.8 $\times$ 10 ${}^{-6}$
0.9	2.11202 $\times$ 10 ${}^{-5}$	4.5 $\times$ 10 ${}^{-3}$	5.1 $\times$ 10 ${}^{-4}$	2.1 $\times$ 10 ${}^{-7}$	2.2 $\times$ 10 ${}^{-5}$

Table 2

The root mean squares errors with $\Delta t=$ 001 at different time level

$t$	$h=$ 0.1		$h=$ 0.05		$h=$ 0.025
	Present	Expo-MCB-DQM [17]	Present	Expo-MCB-DQM [17]	Present	Expo-MCB-DQM [17]
0.1	1.6599 $\times$ 10 ${}^{-4}$	1.66211 $\times$ 10 ${}^{-4}$	2.7121 $\times$ 10 ${}^{-6}$	2.54808 $\times$ 10 ${}^{-6}$	5.4896 $\times$ 10 ${}^{-8}$	2.25582 $\times$ 10 ${}^{-7}$
0.2	1.2141 $\times$ 10 ${}^{-4}$	1.21649 $\times$ 10 ${}^{-4}$	2.2794 $\times$ 10 ${}^{-6}$	2.08096 $\times$ 10 ${}^{-6}$	6.3006 $\times$ 10 ${}^{-8}$	2.85001 $\times$ 10 ${}^{-7}$
0.3	9.1353 $\times$ 10 ${}^{-5}$	9.15324 $\times$ 10 ${}^{-5}$	2.7554 $\times$ 10 ${}^{-6}$	2.62743 $\times$ 10 ${}^{-6}$	5.8357 $\times$ 10 ${}^{-7}$	6.60018 $\times$ 10 ${}^{-7}$
0.4	7.4353 $\times$ 10 ${}^{-5}$	7.44655 $\times$ 10 ${}^{-5}$	6.5909 $\times$ 10 ${}^{-6}$	6.53605 $\times$ 10 ${}^{-6}$	1.7348 $\times$ 10 ${}^{-6}$	1.76211 $\times$ 10 ${}^{-6}$
0.5	5.0076 $\times$ 10 ${}^{-5}$	5.01330 $\times$ 10 ${}^{-5}$	5.5877 $\times$ 10 ${}^{-6}$	5.56751 $\times$ 10 ${}^{-6}$	1.4878 $\times$ 10 ${}^{-6}$	1.50829 $\times$ 10 ${}^{-6}$
0.6	3.3564 $\times$ 10 ${}^{-5}$	3.35858 $\times$ 10 ${}^{-5}$	5.5121 $\times$ 10 ${}^{-6}$	5.49092 $\times$ 10 ${}^{-6}$	1.4736 $\times$ 10 ${}^{-6}$	1.48677 $\times$ 10 ${}^{-6}$
0.7	1.9093 $\times$ 10 ${}^{-5}$	1.91012 $\times$ 10 ${}^{-5}$	4.2303 $\times$ 10 ${}^{-6}$	4.22667 $\times$ 10 ${}^{-6}$	1.1309 $\times$ 10 ${}^{-6}$	1.13366 $\times$ 10 ${}^{-6}$
0.8	6.0154 $\times$ 10 ${}^{-6}$	6.01720 $\times$ 10 ${}^{-6}$	1.3013 $\times$ 10 ${}^{-6}$	1.30143 $\times$ 10 ${}^{-6}$	3.3699 $\times$ 10 ${}^{-7}$	3.37305 $\times$ 10 ${}^{-7}$
0.9	1.2081 $\times$ 10 ${}^{-6}$	1.208670 $\times$ 10 ${}^{-6}$	2.0858 $\times$ 10 ${}^{-7}$	2.08647 $\times$ 10 ${}^{-7}$	5.1918 $\times$ 10 ${}^{-8}$	5.19397 $\times$ 10 ${}^{-8}$
1.0	2.4988 $\times$ 10 ${}^{-6}$	2.497980 $\times$ 10 ${}^{-7}$	2.0523 $\times$ 10 ${}^{-8}$	2.05296 $\times$ 10 ${}^{-8}$	4.9152 $\times$ 10 ${}^{-9}$	4.91607 $\times$ 10 ${}^{-9}$

Figure 2.

Absolute error norms for $h=$ 0.05 at $t=$ 1.

Figure 3.

The RMS error norm for Ex. 1.

Figure 4.

The global % error norms with $\Delta t=$ 0001 at $t=$ 1.

Figure 5.

The global % error norm for $\Delta t=$ 0.005, 0.001, 0.002, 0.005 at $t=$ 1.

Example 2: Here, $\alpha^{\prime}_{x}=\alpha^{\prime}_{y}=\alpha^{\prime}_{z}=\alpha$ and $\beta^{\prime}_{x}=\beta^{\prime}_{y}=\beta^{\prime}_{z}=1$ over the similar domain as in Example 1, with the initial condition: $u({x,y,z,0})=\exp\{{x+y+z}\}$ and boundary conditions are extracted from the exact solution $u({x,y,z,t})=\exp\{{-3({1-\alpha})x+y+z}\}$ . This example is solved for different values of $\Delta t$ , $h$ and $\alpha$ . The global % error norm for $\alpha=$ 0.005, 0.01, 0.5 at $t=$ 1 are shown in Fig. 3 with $\Delta t=$ 0.001 and for different values of $h$ . The global % for $\alpha=$ 0.5 at $t=$ 1 with $h=$ 0.1 for $\Delta t=$ 0.0005, 0.001, 0.002 and 0.005 is depicted in Fig. 4. It is evident from Fig. 4 that the global % error decreases with increasing $\alpha$ and decreasing step size $h$ . Further, from Fig. 5 it is observed that the global % error decreases with decreasing time step size. It is also observed that present results are in good agreement with the results presented in [18].

6. Conclusions

The approximate solution of 3D CDE using DQM based on the modified CTB basis functions has been obtained. Two examples have been solved in order to validate the accuracy and efficiency of the method. It has been found that present method gives highly accurate results in comparison of results presented in [15, 17, 18]. The stability analysis shows that the present method is unconditionally stable. Moreover, the implementation of the method is straightforward and reasonable in terms of data complexity and requires low memory storage.

References

Dawson

S.P.

Chen

and Doolen

G.D.

, Lattice boltzmann computations for reaction-diffusion equations, J. Chem. Phys. 98(2) (1993), 1514–1523.

Ding

and Zhang

, A new difference scheme with high accuracy and absolute stability for solving convection-diffusion equations, J. Comput. Appl. Math. 230 (2009), 600–606.

Tian

Z.F.

and Yua

P.X.

, A high-order exponential scheme for solving 1D unsteady convection-diffusion equations, J. Comput. Appl. Math. 235 (2011), 2477–2491.

J.-H.

, A simple approach to one-dimensional convection-diffusion equation and its fractional modification for E reaction arising in rotating disk electrodes, Journal of Electroanalytical Chemistry 854 (2019), 113565.

Dhawan

Kapoor

and Kumar

, Numerical method for advection diffusion equation using FEM and B-splines, J. Comput. Sci. 3 (2012), 429–437.

Mohebbi

and Dehghan

, High-order compact solution of the one-dimensional heat and advection diffusion equations, Appl. Math. Model. 34 (2010), 3071–3084.

Noye

B.J.

, A compact unconditionally stable finite-difference method for transient one-dimensional advection-diffusion, Commun. Appl. Numer. Methods 7 (1991), 501–512.

Singh

Das

Ong

S.H.

and Jafari

, Numerical solution of nonlinear reaction-advection-diffusion equation, Journal of Computational and Nonlinear Dynamics, 2019. doi: 10.1115/1.4042687.

Tian

Z.F.

and Ge

Y.B.

, A fourth-order compact ADI method for solving two-dimensional unsteady convection-diffusion problems, J. Comput. Appl. Math. 198 (2007), 268–286.

10.

Dehghan

and Mohebbi

, High-order compact boundary value method for the solution of unsteady convection-diffusion problems, Math. Comput. Simulat. 79 (2008), 683–699.

11.

Karaa

and Zhang

, High order ADI method for solving unsteady convection-diffusion problems, J. Comput. Phys. 198 (2004), 1–9.

12.

Kalita

J.C.

Dalal

D.C.

and Dass

A.K.

, A class of higher order compact schemes for the unsteady two dimensional convection-diffusion equation with variable convection coefficients, Int. J. Numer. Methods Fluids 38 (2002), 1111–1131.

13.

Tian

Z.F.

, A rational high-order compact ADI method for unsteady convection-diffusion equations, Computer Phys. Comm. 182 (2011), 649662.

14.

Cecchi

M.M.

and Pirozzi

M.A.

, High order finite difference numerical methods for time-dependent convection-dominated problems, Appl. Numer. Math. 55 (2005), 334–356.

15.

Dehghan

, Numerical solution of the three-dimensional advection-diffusion equation, Appl. Math. Comput. 150 (2004), 5–19.

16.

Mohanty

R.K.

and Singh

, A new highly accurate discretization for three dimensional singularly perturbed non-linear elliptic partial differential equation, Numerical Methods for Partial Differential Equations 22(6) (2006), 1379–1395.

17.

Shukla

H.S.

and Tamsir

, An exponential cubic B-spline algorithm for multi-dimensional convection-diffusion equations, Alexandria Engineering Journal 57 (2018), 1999–2006. doi: 10.1016/j.aej.2017.04.011.

18.

Prieto

F.U.

Benito Munozb

J.J.

and Corvinos

L.G.

, Application of the generalized finite difference method to solve the advection diffusion equation, J. Comput. Appl. Math. 235 (2011), 1849–1855.

19.

Ravnik

and Tibuat

, Fast boundary-domain integral method for unsteady convection-diffusion equation with variable diffusivity using the modified helmholtz fundamental solution, Numerical Algorithms 82 (2019), 1441–1466.

20.

Bellman

Kashef

B.G.

and Casti

, Differential quadrature: A technique for the rapid solution of nonlinear differential equations, J. Comput. Phy. 10 (1972), 40–52.

21.

Korkmaz

, Shock wave simulations using sinc differential quadrature method, Int. J. Comput. Aided Engg Software 28(6) (2011), 654–674.

22.

Quan

J.R.

and Chang

C.T.

, New insights in solving distributed system equations by the quadrature methods-I, Comput. Chem. Eng. 13 (1989), 779–788.

23.

Quan

J.R.

and Chang

C.T.

, New insights in solving distributed system equations by the quadrature methods-II, Comput. Chem. Eng. 13 (1989), 1017–1024.

24.

Shu

and Richards

B.E.

, Application of generalized differential quadrature to solve two dimensional incompressible navier-stokes equations, Int. J. Numer. Meth. Fluids 15 (1992), 791–798.

25.

Shu

, Differential Quadrature and its Application in Engineering, Athenaeum Press Ltd., Great Britain, 2000.

26.

Spiteri

J.R.

and Ruuth

S.J.

, A new class of optimal high-order strong stability-preserving time-stepping schemes, SIAM Journal Numer Anal 40(2) (2002), 469–491.

27.

Tamsir

Dhiman

and Srivastava

V.K.

, Cubic trigonometric B-spline differential quadrature method for numerical treatment of fisher’s reaction-diffusion equations, Alexandria Engineering Journal 57 (2018), 2019–2026. doi: 10.1016/j.aej.2017.05.007.

28.

Mittal

R.C.

Kumar

and Jiwari

, A Comparative Study of Cubic B-spline-Based Quasi-interpolation and Differential Quadrature Methods for Solving Fourth-Order Parabolic PDEs, Proceedings of the National Academy of Sciences, India Section A: Physical Sciences, 2020. doi: 10.1007/s40010-020-00684-y.

29.

Abbas

Majid

A.A.

Ismail

A.I.M.

and Rashid

, The application of cubic trigonometric B-spline to the numerical solution of the hyperbolic problems, Appl. Math. Comput. 239 (2014), 74–88.

30.

Alshomrani

A.S.

Pandit

Alzahrani

A.K.

Alghamdi

M.S.

and Jiwari

, A numerical algorithm based on modified cubic trigonometric B-spline functions for computational modelling of hyperbolic-type wave equations, Engineering Computations 34(4) (2017), 1257–1276.

31.

Jiwari

and Alshomrani

A.S.

, A new algorithm based on modified trigonometric cubic B-splines functions for nonlinear burgers’-type equations, International Journal of Numerical Methods for Heat and Fluid Flow 27(8) (2017), 1638–1661.

32.

Jiwari

Pandit

and Koksal

M.E.

, A class of numerical algorithms based on cubic trigonometric B-spline functions for numerical simulation of nonlinear parabolic problems, Computational and Applied Mathematics 38(3) (2019), 140.

33.

Dhiman

and Tamsir

, A collocation technique based on modified form of trigonometric cubic B-spline basis functions for fisher’s reaction-diffusion equation, Multidiscipline Modeling in Materials and Structures 14(5) (2018), 923–939.

34.

Shukla

H.S.

Tamsir

Jiwari

and Srivastava

V.K.

, A numerical algorithm for computation modeling of 3D nonlinear wave equations based on exponential modified cubic B-spline differential quadrature method, International Journal of Computer Mathematics 95(4) (2017), 752–766.

Approximation of 3D convection diffusion equation using DQM based on modified cubic trigonometric B-splines

Abstract

Keywords

1. Introduction

Table 1 Comparative study for absolute error norms

References

Table 1
Comparative study for absolute error norms