# Multidimensional Convection and Diffusion Problems

## 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$ (1)

which can be rewritten as $\frac{\partial (\rho \varphi )}{\partial t}+\frac{\partial J_{x}}{\partial x}+\frac{\partial J_{y}}{\partial y}=S$ (2)

where $J_{x}=\rho u\varphi -\Gamma \frac{\partial \varphi }{\partial x}$ (3) $J_{y}=\rho v\varphi -\Gamma \frac{\partial \varphi }{\partial y}$ (4)

Integrating eq. (2) with respect to t in the interval of (t, t+Δt) and over the control volume P, we have \begin{align} & \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}}} \\ \end{align}

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$ (5)

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}}$

where $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 \begin{align} & \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} \\ \end{align} (6)

where $\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$ (7)

where $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\}$ (8) $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\}$ (9) $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\}$ (10) $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\}$ (11) $b=a_{P}^{0}\varphi _{P}^{0}+S_{C}\Delta x\Delta y$ (12) \begin{align} & 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}) \\ \end{align} (13) $a_{P}^{0}=\frac{\rho _{P}^{0}\Delta x\Delta y}{\Delta t}$ (14)

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$ (15)

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.

## 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 (Patankar, 1980) $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$ (16)

where $a_{E}=D_{e}\left\{ A(\left| \text{Pe}_{\Delta e} \right|)+\left[\!\left[ -\text{Pe}_{\Delta e},0 \right]\!\right] \right\}$ (17) $a_{W}=D_{w}\left\{ A(\left| \text{Pe}_{\Delta w} \right|)+\left[\!\left[ \text{Pe}_{\Delta w},0 \right]\!\right] \right\}$ (18) $a_{N}=D_{n}\left\{ A(\left| \text{Pe}_{\Delta n} \right|)+\left[\!\left[ -\text{Pe}_{\Delta n},0 \right]\!\right] \right\}$ (19) $a_{S}=D_{s}\left\{ A(\left| \text{Pe}_{\Delta s} \right|)+\left[\!\left[ \text{Pe}_{\Delta s},0 \right]\!\right] \right\}$ (20) $a_{T}=D_{t}\left\{ A(\left| \text{Pe}_{\Delta t} \right|)+\left[\!\left[ -\text{Pe}_{\Delta t},0 \right]\!\right] \right\}$ (21) $a_{B}=D_{b}\left\{ A(\left| \text{Pe}_{\Delta b} \right|)+\left[\!\left[ \text{Pe}_{\Delta b},0 \right]\!\right] \right\}$ (22) $b=a_{P}^{0}\varphi _{P}^{0}+S_{C}\Delta x\Delta y\Delta z$ (23) $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$ (24) $a_{P}^{0}=\frac{\rho \Delta x\Delta y\Delta z}{\Delta t}$ (25)

The expressions for conductance at the faces of the control volume are \begin{align} & 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}} \\ \end{align} (26)

and the flow rates are: \begin{align} & 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 \\ \end{align} (27)

The Different discretization schemes for the above three-dimensional problem can be obtained by using different expressions for A(|PeΔ|) from Table 4.3. In addition to the six first order discretization schemes described above, some researchers have used higher order schemes such as second order upwind (Leonard et al., 1981) and QUICK (Quadratic Upwind Interpolation of Convective Kinetics; Leonard, 1979) 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 (Patankar, 1980). 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 (Tao, 2001). 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.