Numerical solution of the convection–diffusion equation

The following convection diffusion equation is considered here[2]

${\displaystyle c\rho \left[{\frac {\partial T(x,t)}{\partial t}}+\epsilon u{\frac {\partial T(x,t)}{\partial x}}\right]=\lambda {\frac {\partial ^{2}T(x,t)}{\partial x^{2}}}+Q(x,t)}$

In the above equation, four terms represents transience, convection, diffusion and a source term respectively, where

• T is the temperature in particular case of heat transfer otherwise it is the variable of interest
• t is time
• c is the specific heat
• u is velocity
• ε is porosity that is the ratio of liquid volume to the total volume
• ρ is mass density
• λ is thermal conductivity
• Q(x,t) is source term representing the capacity of internal sources

The equation above can be written in the form

${\displaystyle {\frac {\partial T}{\partial t}}=a{\frac {\partial ^{2}T}{\partial x^{2}}}-\epsilon u{\frac {\partial T}{\partial x}}+{\frac {Q}{c\rho }}}$

where a = λ/cρ is the diffusion coefficient.

A solution of the transient convection–diffusion equation can be approximated through a finite difference approach, known as the finite difference method (FDM).

An explicit scheme of FDM has been considered and stability criteria are formulated. In this scheme, temperature is totally dependent on the old temperature (the initial conditions) and θ, a weighting parameter between 0 and 1. Substitution of θ = 0 gives the explicit discretization of the unsteady conductive heat transfer equation.

${\displaystyle {\frac {T_{i}^{f}-T_{i}^{f-1}}{\Delta t}}=a{\frac {T_{i-1}^{f-1}-2T_{i}^{f-1}+T_{i+1}^{f-1}}{h^{2}}}-\epsilon u{\frac {T_{i+1}^{f-1}-T_{i-1}^{f-1}}{2h}}+{\frac {Q_{i}^{f-1}}{c\rho }}}$

where

• Δt = t^f − t^{f − 1}
• h is the uniform grid spacing (mesh step)

${\displaystyle T_{i}^{f}=\left(1-{\frac {2a\Delta t}{h^{2}}}\right)T_{i}^{f-1}+\left({\frac {a\Delta t}{h^{2}}}+{\frac {\epsilon u\Delta t}{2h}}\right)T_{i-1}^{f-1}+\left({\frac {a\Delta t}{h^{2}}}-{\frac {\epsilon u\Delta t}{2h}}\right)T_{i+1}^{f-1}+{\frac {Q_{i}^{f-1}}{c\rho }}\Delta t}$

${\displaystyle {\begin{aligned}h&<{\frac {2a}{\epsilon u}}\\\Delta t&<{\frac {h^{2}}{2a}}\end{aligned}}}$

These inequalities set a stringent maximum limit to the time step size and represents a serious limitation for the explicit scheme. This method is not recommended for general transient problems because the maximum possible time step has to be reduced as the square of h.

In implicit scheme, the temperature is dependent at the new time level t + Δt. After using implicit scheme, it was found that all coefficients are positive. It makes the implicit scheme unconditionally stable for any size of time step. This scheme is preferred for general purpose transient calculations because of its robustness and unconditional stability.[3] The disadvantage of this method is that more procedures are involved and due to larger Δt, truncation error is also larger.

In the Crank–Nicolson method, the temperature is equally dependent on t and t + Δt. It is a second-order method in time and this method is generally used in diffusion problems.

${\displaystyle \Delta t<{\frac {h^{2}}{a}}}$

This time step limitation is less restricted than the explicit method. The Crank–Nicolson method is based on the central differencing and hence it is second-order accurate in time.[4]