Multidimensional Convection and Diffusion Problems

From Thermal-FluidsPedia

Revision as of 06:35, 14 July 2010 by Founding Editor (Talk | contribs)
Jump to: navigation, search
 Related Topics Catalog
Computational methodologies for forced convection
  1. One-Dimensional Steady-State Convection and Diffusion
    1. Central Difference Scheme
    2. Upwind Scheme
    3. Hybrid Scheme
    4. Exponential and Power Law Schemes
    5. A Generalized Expression of Discretization Schemes
  2. Multidimensional Convection and Diffusion Problems
  3. Numerical Solution of Flow Field
    1. Special Difficulties
    2. Staggered grid
    3. Pressure Correction Equation
    4. The SIMPLE Algorithm
  4. Numerical Simulation of Interfaces and Free Surfaces
  5. Application of Computational Methods

Two-dimensional problem

The heat transfer problems discussed in the preceding subsection are steady-state convection-diffusion problems with the general variable \varphi varying in one dimension only. We now turn our attention to the unsteady state two-dimensional convection-diffusion problem which includes a source term S. The problem is described by

\frac{\partial (\rho \varphi )}{\partial t}+\frac{\partial (\rho u\varphi )}{\partial x}+\frac{\partial (\rho v\varphi )}{\partial y}=\frac{\partial }{\partial x}\left( \Gamma \frac{\partial \varphi }{\partial x} \right)+\frac{\partial }{\partial y}\left( \Gamma \frac{\partial \varphi }{\partial y} \right)+S


which can be rewritten as

\frac{\partial (\rho \varphi )}{\partial t}+\frac{\partial J_{x}}{\partial x}+\frac{\partial J_{y}}{\partial y}=S



J_{x}=\rho u\varphi -\Gamma \frac{\partial \varphi }{\partial x}


J_{y}=\rho v\varphi -\Gamma \frac{\partial \varphi }{\partial y}


Integrating eq. (2) with respect to t in the interval of (t, t+Δt) and over the control volume P, we have

  & \int_{s}^{n}{\int_{e}^{w}{\int_{t}^{t+\Delta t}{\frac{\partial }{\partial t}\left( \rho \varphi  \right)dt}dxdy}}+\int_{t}^{t+\Delta t}{\int_{s}^{n}{\int_{w}^{e}{\frac{\partial J_{x}}{\partial x}dxdydt}}}+\int_{t}^{t+\Delta t}{\int_{e}^{w}{\int_{s}^{n}{\frac{\partial J_{y}}{\partial y}dydxdt}}} \\ 
 & =\int_{t}^{t+\Delta t}{\int_{s}^{n}{\int_{w}^{e}{(S_{C}+S_{P}\varphi )dxdydt}}} \\ 

where the source term is treated as a linear function of \varphi . Assuming the total fluxes are uniform on all faces of the control volume and employing fully-implicit scheme, the above equation becomes

(\rho _{P}\varphi _{P}-\rho _{P}^{0}\varphi _{P}^{0})\Delta x\Delta y+(J_{x}^{e}-J_{x}^{w})\Delta y\Delta t+(J_{y}^{n}-J_{y}^{s})\Delta x\Delta t=(S_{C}+S_{P}\varphi _{P})\Delta x\Delta y\Delta t

where the superscript 0 represents the values at the previous time step. Introducing the integrated total fluxes J_{e}=J_{x}^{e}\Delta y, J_{w}=J_{x}^{w}\Delta y, J_{n}=J_{y}^{n}\Delta x and J_{s}=J_{y}^{s}\Delta x, and dividing the above equation by Δt yields

\frac{(\rho _{P}\varphi _{P}-\rho _{P}^{0}\varphi _{P}^{0})}{\Delta t}\Delta x\Delta y+(J_{e}-J_{w})+(J_{n}-J_{s})=(S_{C}+S_{P}\varphi _{P})\Delta x\Delta y


Defining the following integrated total flux

J_{e}^{*}=\frac{J_{e}}{D_{e}},\text{ }J_{w}^{*}=\frac{J_{w}}{D_{w}},\text{ }J_{n}^{*}=\frac{J_{n}}{D_{n}},\text{ }J_{s}^{*}=\frac{J_{s}}{D_{s}}


D_{e}=\frac{\Gamma _{e}}{(\delta x)_{e}}\Delta y,\text{ }D_{w}=\frac{\Gamma _{w}}{(\delta x)_{w}}\Delta y,\text{ }D_{n}=\frac{\Gamma _{n}}{(\delta y)_{n}}\Delta x,\text{ }D_{s}=\frac{\Gamma _{s}}{(\delta y)_{s}}\Delta x

the integrated fluxes at the east and west faces of the control volume can be evaluated:

J_{e}=D_{e}J_{e}^{*}=D_{e}[B(\text{Pe}_{\Delta e})\varphi _{P}-A(\text{Pe}_{\Delta e})\varphi _{E}]

J_{w}=D_{w}J_{w}^{*}=D_{w}[B(\text{Pe}_{\Delta w})\varphi _{W}-A(\text{Pe}_{\Delta w})\varphi _{P}]

Substituting B(PeΔ) − A(PeΔ) = PeΔ into the two equations above yields

J_{e}=D_{e}J_{e}^{*}=D_{e}[A(\text{Pe}_{\Delta e})\varphi _{P}+\text{Pe}_{\Delta e}\varphi _{P}-A(\text{Pe}_{\Delta e})\varphi _{E}]

J_{w}=D_{w}J_{w}^{*}=D_{w}[B(\text{Pe}_{\Delta w})\varphi _{W}-B(\text{Pe}_{\Delta w})\varphi _{P}+\text{Pe}_{\Delta w}\varphi _{P}]

Similarly, the integrated total flux at the north and south faces of the control volume can be expressed as

J_{n}=D_{n}J_{n}^{*}=D_{n}[A(\text{Pe}_{\Delta n})\varphi _{P}+\text{Pe}_{\Delta n}\varphi _{P}-A(\text{Pe}_{\Delta n})\varphi _{N}]

J_{s}=D_{s}J_{s}^{*}=D_{s}[B(\text{Pe}_{\Delta s})\varphi _{S}-B(\text{Pe}_{\Delta s})\varphi _{P}+\text{Pe}_{\Delta s}\varphi _{P}]

Substituting the above four integrated total fluxes into eq. (5), we have

  & \left\{ \rho _{P}\Delta x\Delta y/\Delta t-S_{P}\Delta x\Delta y+D_{e}A(\text{Pe}_{\Delta e})+D_{w}B(\text{Pe}_{\Delta w}) \right. \\ 
 & \left. +D_{n}A(\text{Pe}_{\Delta n})+D_{s}B(\text{Pe}_{\Delta s})+(F_{e}-F_{w})+(F_{n}-F_{s}) \right\}\varphi _{P} \\ 
 & =D_{e}A(\text{Pe}_{\Delta e})\varphi _{E}+D_{w}B(\text{Pe}_{\Delta w})\varphi _{W}+D_{n}A(\text{Pe}_{\Delta n})\varphi _{N} \\ 
 & +D_{s}B(\text{Pe}_{\Delta s})\varphi _{S}+S_{C}\Delta x\Delta y+(\rho _{P}^{0}\Delta x\Delta y/\Delta t)\varphi _{P}^{0} \\ 



\begin{matrix}{}\\\end{matrix}F_{e}=D_{e}\text{Pe}_{\Delta e}=(\rho u)_{e}\Delta y,\text{ }F_{w}=D_{w}\text{Pe}_{\Delta w}=(\rho u)_{w}\Delta y

\begin{matrix}{}\\\end{matrix}F_{n}=D_{n}\text{Pe}_{\Delta n}=(\rho v)_{n}\Delta x,\text{ }F_{s}=D_{s}\text{Pe}_{\Delta s}=(\rho v)_{s}\Delta x

are the mass flow rates at the four faces of the control volume. Equation (6) can be rearranged to obtain the final form of the following discretized equation:

a_{P}\varphi _{P}=a_{E}\varphi _{E}+a_{W}\varphi _{W}+a_{N}\varphi _{N}+a_{S}\varphi _{S}+b



a_{E}=D_{e}A(\text{Pe}_{\Delta e})=D_{e}\left\{ A(\left| \text{Pe}_{\Delta e} \right|)+\left[\!\left[ -\text{Pe}_{\Delta e},0 \right]\!\right] \right\}


a_{W}=D_{w}B(\text{Pe}_{\Delta w})=D_{w}\left\{ A(\left| \text{Pe}_{\Delta w} \right|)+\left[\!\left[ \text{Pe}_{\Delta w},0 \right]\!\right] \right\}


a_{N}=D_{n}A(\text{Pe}_{\Delta n})=D_{n}\left\{ A(\left| \text{Pe}_{\Delta n} \right|)+\left[\!\left[ -\text{Pe}_{\Delta n},0 \right]\!\right] \right\}


a_{S}=D_{s}B(\text{Pe}_{\Delta s})=D_{s}\left\{ A(\left| \text{Pe}_{\Delta s} \right|)+\left[\!\left[ \text{Pe}_{\Delta s},0 \right]\!\right] \right\}


b=a_{P}^{0}\varphi _{P}^{0}+S_{C}\Delta x\Delta y


  & a_{P}=a_{E}+a_{W}+a_{N}+a_{S}+(\rho _{P}/\rho _{P}^{0})a_{P}^{0} \\ 
 & \text{       }-S_{P}\Delta x\Delta y+(F_{e}-F_{w})+(F_{n}-F_{s}) \\ 


a_{P}^{0}=\frac{\rho _{P}^{0}\Delta x\Delta y}{\Delta t}


If the continuity equation is satisfied, eq. (13) can be simplified as

a_{P}=a_{E}+a_{W}+a_{N}+a_{S}+a_{P}^{0}-S_{P}\Delta x\Delta y


Similar to the case of one-dimensional convection-diffusion, different discretization schemes for the discretized equations (7) – (14) can be obtained by using different expressions for A(|PeΔ|) from the following table.

Table Summary of A(|PeΔ|) for different schemes

Scheme A(|PeΔ|)
Central difference 1-0.5\left| \text{Pe}_{\Delta } \right|
Upwind 1
Hybrid \left[\!\left[ 0,1-0.5\left| \text{Pe}_{\Delta } \right| \right]\!\right]
Exponential \left| \text{Pe}_{\Delta } \right|/[\exp (\left| \text{Pe}_{\Delta } \right|)-1]
Power Law \left[\!\left[ 0,(1-0.1\left| \text{Pe}_{\Delta } \right|)^{5} \right]\!\right]

Three-dimensional problem

The discretized equation for a transient three-dimensional convection-diffusion problem can be obtained by integrating the conservation equation with respect to t in the interval of (t, t+Δt) and over the three-dimensional control volume P (formed by considering two additional neighbors at top, T, and bottom, B). The final form of the governing equation is [1]

a_{P}\varphi _{P}=a_{E}\varphi _{E}+a_{W}\varphi _{W}+a_{N}\varphi _{N}+a_{S}\varphi _{S}+a_{T}\varphi _{T}+a_{B}\varphi _{B}+b



a_{E}=D_{e}\left\{ A(\left| \text{Pe}_{\Delta e} \right|)+\left[\!\left[ -\text{Pe}_{\Delta e},0 \right]\!\right] \right\}


a_{W}=D_{w}\left\{ A(\left| \text{Pe}_{\Delta w} \right|)+\left[\!\left[ \text{Pe}_{\Delta w},0 \right]\!\right] \right\}


a_{N}=D_{n}\left\{ A(\left| \text{Pe}_{\Delta n} \right|)+\left[\!\left[ -\text{Pe}_{\Delta n},0 \right]\!\right] \right\}


a_{S}=D_{s}\left\{ A(\left| \text{Pe}_{\Delta s} \right|)+\left[\!\left[ \text{Pe}_{\Delta s},0 \right]\!\right] \right\}


a_{T}=D_{t}\left\{ A(\left| \text{Pe}_{\Delta t} \right|)+\left[\!\left[ -\text{Pe}_{\Delta t},0 \right]\!\right] \right\}


a_{B}=D_{b}\left\{ A(\left| \text{Pe}_{\Delta b} \right|)+\left[\!\left[ \text{Pe}_{\Delta b},0 \right]\!\right] \right\}


b=a_{P}^{0}\varphi _{P}^{0}+S_{C}\Delta x\Delta y\Delta z


a_{P}=a_{E}+a_{W}+a_{N}+a_{S}+a_{T}+a_{B}+a_{P}^{0}-S_{P}\Delta x\Delta y\Delta z


a_{P}^{0}=\frac{\rho \Delta x\Delta y\Delta z}{\Delta t}


The expressions for conductance at the faces of the control volume are

  & D_{e}=\frac{\Gamma _{e}\Delta y\Delta z}{(\delta x)_{e}},\text{  }D_{w}=\frac{\Gamma _{w}\Delta y\Delta z}{(\delta x)_{w}},\text{   }D_{n}=\frac{\Gamma _{n}\Delta x\Delta z}{(\delta y)_{n}} \\ 
 & D_{s}=\frac{\Gamma _{s}\Delta x\Delta z}{(\delta y)_{s}},\text{   }D_{t}=\frac{\Gamma _{t}\Delta x\Delta y}{(\delta z)_{t}},\text{  }D_{b}=\frac{\Gamma _{b}\Delta x\Delta y}{(\delta z)_{b}} \\ 


and the flow rates are:

  & F_{e}=D_{e}\text{Pe}_{\Delta e}=(\rho u)_{e}\Delta y\Delta z,\text{  }F_{w}=D_{w}\text{Pe}_{\Delta w}=(\rho u)_{w}\Delta y\Delta z \\ 
 & F_{n}=D_{n}\text{Pe}_{\Delta n}=(\rho v)_{n}\Delta x\Delta z,\text{  }F_{s}=D_{s}\text{Pe}_{\Delta s}=(\rho v)_{s}\Delta x\Delta z \\ 
 & F_{t}=D_{t}\text{Pe}_{\Delta t}=(\rho w)_{t}\Delta x\Delta y,\text{  }F_{b}=D_{b}\text{Pe}_{\Delta b}=(\rho w)_{b}\Delta x\Delta y \\ 


The Different discretization schemes for the above three-dimensional problem can be obtained by using different expressions for A(|PeΔ|) from the above table. In addition to the six first order discretization schemes described above, some researchers have used higher order schemes such as second order upwind [2] and QUICK (Quadratic Upwind Interpolation of Convective Kinetics)[3] schemes to overcome the false diffusion problem, which is referred to as error caused by using the discretization scheme with accuracy less than the second order [1]. The error due to false diffusion could potentially be severe for (1) transient problems, (2) multidimensional steady-state problems, or (3) problems with non-constant source terms [4]. While the accuracies of these higher order schemes are better than the first order schemes, their computational time is much greater than that of the first order schemes.


  1. 1.0 1.1 Patankar, S.V., 1980, Numerical Heat Transfer and Fluid Flow, Hemisphere, Washington, DC.
  2. Leonard 1981, B.P., “A Survey of Finite Differences with Upwinding for Numerical Modeling of the Incompressible Convection Diffusion Equation,” Computational Techniques in Transient and Turbulent Flows, Taylor, C., and Morgan, K., Eds., Pineridge Press, Swansea, pp. 1-35.
  3. Leonard, B.P., 1979, “A Stable and Accurate Convective Modeling Procedure based on Quadratic Upstream Interpolation,” Computer Methods in Applied Mechanics and Engineering, Vol. 29, pp. 59-98.
  4. Tao 2001, W.Q., Numerical Heat Transfer, 2nd Ed., Xi’an Jiaotong University Press, Xi’an, China (in Chinese).