Finite volume method for three-dimensional diffusion problem

Finite volume method (FVM) is a numerical method. FVM in computational fluid dynamics is used to solve the partial differential equation which arises from the physical conservation law by using discretisation. Convection is always followed by diffusion and hence where convection is considered we have to consider combine effect of convection and diffusion. But in places where fluid flow plays a non-considerable role we can neglect the convective effect of the flow. In this case we have to consider more simplistic case of only diffusion. The general equation for steady convection-diffusion can be easily derived from the general transport equation for property $\phi$ by deleting transient.

General transport equation is defined as:

${\frac {\partial \rho \phi }{\partial t}}+\operatorname {div} (\rho \phi \upsilon )=\operatorname {div} (\Gamma \operatorname {grad} \phi )+S_{\phi }$ …………………………………………….1

Where,

$\phi$ is a conservative form of all fluid flow,

$\rho$ is density,

$\operatorname {div} (\rho \phi \upsilon )$ is a net rate of flow of $\phi$ out of fluid element represents convective term,

${\frac {\partial \rho \phi }{\partial t}}$ is a transient term,

$\operatorname {div} (\Gamma \operatorname {grad} \phi )$ is a rate of change $\phi$ of due to diffusion,

$S_{\phi }$ is a rate of increase of $\phi$ due to source.

Due to steady state condition transient term becomes zero and due to absence of convection convective term becomes zero, therefore steady state three- dimensional convection and diffusion equation becomes:

$\operatorname {div} (\Gamma \operatorname {grad} \phi )+S_{\phi }=0$ ………………………………………………………….2

Therefore,

${\frac {\partial }{\partial x}}\left(\Gamma {\frac {\partial \phi }{\partial x}}\right)+{\frac {\partial }{\partial y}}\left(\Gamma {\frac {\partial \phi }{\partial y}}\right)+{\frac {\partial }{\partial z}}\left(\Gamma {\frac {\partial \phi }{\partial z}}\right)+S_{\phi }=0$ …………………………………………………………………….3

Flow should also satisfy continuity equation therefore,

$\nabla \cdot (\rho \mathbf {u} )=0$ ………………………………………………………………………………………………………4

To solve the problem we will follow following general steps

Grid formation:

1. Divide the domain into discrete control volume.

2. Place the nodal point between end points defining the physical boundaries. Boundaries/ faces of the control volume are created midway between adjacent nodes.

3. Set up the control volume near the edge of domain such that physical as well as control volume boundaries will coincide with each other.

4. Considering a general nodal point P accompanied by six neighboring nodal point ‘E’ (which represent east), ‘W’ (which represent west), ‘N’ (which represent north), ‘S’ (which represent south), ‘T’ (which represent Top), ‘B’ (which represent bottom). In the considered control volume east side face is referred by ‘e’, west side face is referred by ‘w’, north side face is referred by ‘n’, south side face is referred by ‘s’, top side face is referred by ‘t’, bottom side face is referred by ‘b’.

5. Now the distance between nodes W and P, between nodes P and E, between nodes P and N, between nodes S and P, between nodes P and T, between nodes B and P are denoted as ${\delta x_{WP}},{\delta x_{PE}},{\delta x_{PN}},{\delta x_{SP}},{\delta x_{PT}},{\delta x_{BP}}$ respectively.

Discretisation:

On integration of equation 3 in one dimension over the general control volume gives:

[ $\left(\Gamma A{\frac {d\phi }{dx}}\right)_{e}-\left(\Gamma A{\frac {d\phi }{dx}}\right)_{w}]+[\left(\Gamma A{\frac {d\phi }{dx}}\right)_{n}-\left(\Gamma A{\frac {d\phi }{dx}}\right)_{s}]+[\left(\Gamma A{\frac {d\phi }{dx}}\right)_{t}-\left(\Gamma A{\frac {d\phi }{dx}}\right)_{b}]+{\overrightarrow {S}}\Delta V=0$

Now using central differencing method we can rewrite above equation as

[ $\Gamma _{e}A_{e}\left({\frac {\phi _{E}-\phi _{P}}{\delta x_{PE}}}\right)-\Gamma _{w}A_{w}\left({\frac {\phi _{P}-\phi _{W}}{\delta x_{PE}}}\right)]+[\Gamma _{n}A_{n}\left({\frac {\phi _{E}-\phi _{P}}{\delta x_{PE}}}\right)-\Gamma _{s}A_{s}\left({\frac {\phi _{P}-\phi _{W}}{\delta x_{PE}}}\right)]+[\Gamma _{t}A_{t}\left({\frac {\phi _{E}-\phi _{P}}{\delta x_{PE}}}\right)-\Gamma _{b}A_{b}\left({\frac {\phi _{P}-\phi _{W}}{\delta x_{PE}}}\right)]+(S_{u}+S_{p}\phi _{p})=0$

This can be rearranged to give the discretised equation for interior nodes:

$a_{P}\phi _{P}=a_{W}\phi _{W}+a_{E}\phi _{E}+a_{S}\phi _{S}+a_{N}\phi _{N}+a_{B}\phi _{B}+a_{T}\phi _{T}+S_{u}$

Where

$a_{W}$	$a_{E}$	$a_{S}$	$a_{N}$	$a_{B}$	$a_{T}$	$a_{P}$
${\frac {\Gamma _{w}}{\delta x_{WP}}}A_{w}$	${\frac {\Gamma _{e}}{\delta x_{PE}}}A_{e}$	${\frac {\Gamma _{s}}{\delta x_{SP}}}A_{s}$	${\frac {\Gamma _{n}}{\delta x_{PN}}}A_{n}$	${\frac {\Gamma _{b}}{\delta x_{BP}}}A_{w}$	${\frac {\Gamma _{t}}{\delta x_{PT}}}A_{t}$	$a_{W}+a_{E}+a_{S}+a_{N}+a_{B}+a_{T}-S_{P}$