Finite volume method for one-dimensional steady state diffusion

The Finite volume method in computational fluid dynamics is a discretization technique for partial differential equations that arise from physical conservation laws. These equations can be different in nature, e.g. elliptic, parabolic, or hyperbolic. The first well-documented use of this method was by Evans and Harlow (1957) at Los Alamos. The general equation for steady diffusion can easily be derived from the general transport equation for property Φ by deleting transient and convective terms.

General Transport equation can be defined as

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

where

$$\rho$$ is density and $$ \phi $$ is the conserved quantity,

$$\Gamma$$ is the Diffusion coefficient and $$S$$ is the Source term.

$$\operatorname{div}(\rho \phi \upsilon)$$ is the Net rate of flow of $$ \phi $$ out of fluid element (convection),

$$\operatorname{div}(\Gamma \operatorname{grad} \phi) $$ is Rate of increase of $$ \phi $$ due to diffusion,

$$ S_\phi$$ is Rate of increase of $$\phi$$ due to sources.

$$\frac{\partial \rho \phi }{\partial t} $$ is Rate of increase of $$ \phi $$ of fluid element(transient),

Conditions under which the transient and convective terms goes to zero:
 * Steady State
 * Low Reynolds Number

For one-dimensional, steady-state diffusion, General Transport equation reduces to:


 * $$\operatorname{div}(\Gamma\operatorname{grad}\phi)+ S_\phi=0$$,

or,


 * $$\frac {d}{dx} (\Gamma\operatorname{grad}\phi) +S_\phi=0$$.

The following steps comprise the finite volume method for one-dimensional steady state diffusion -

STEP 1

Grid Generation


 * Divide the domain into equal parts of small domain.
 * Place nodal points at the center of each small domain.
 * Create control volumes using these nodal points.
 * Create control volumes near the edges in such a way that the physical boundaries coincide with control volume boundaries (Figure 1).
 * Assume a general nodal point 'P' for a general control volume. Adjacent nodal points to the East and West are identified by E and W respectively. The West-side face of the control volume is referred to by 'w' and the East-side control volume face by 'e' (Figure 2).
 * The distance between WP, wP, Pe and PE are identified by $$\delta x_{WP}$$,$$\delta x_{wP}$$,$$\delta x_{Pe}$$ and $$\delta x_{PE}$$ respectively (Figure 4).

STEP 2 Discretization


 * The crux of Finite volume method is to integrate the governing equation over each control volume.
 * Nodal points are used to discretize equations.
 * At nodal point P, the control volume integral is given by (Figure 3)

$$\int_{\Delta V} \frac{d}{dx}\left (\Gamma \frac{d \phi}{dx}\right ) dV + \int_{\Delta V} S dV = \left (\Gamma A \frac{d \phi}{dx}\right )_e - \left (\Gamma A \frac{d \phi}{dx}\right )_w + \overrightarrow{S} \Delta V= 0$$ ,

where

$$A$$ is Cross-sectional Area Cross section (geometry) of control volume face, $$ \Delta V$$ is Volume, $$ \overrightarrow{S} $$is average value of source S over the control volume.


 * It states that the difference between the diffusive flux Fick's laws of diffusion of $$ \phi $$ through the east and west faces of some volume corresponds to the change in the quantity $$ \phi $$ in that volume.
 * The diffusive coefficient of $$ \phi $$ and $$ \frac{d \phi}{dx} $$ are required in order to reach a useful conclusion.
 * Central differencing technique  is used to derive the diffusive coefficient of $$ \phi $$:

$$\Gamma_w = \frac{\Gamma_W + \Gamma_P}{2}$$ ,

$$\Gamma_e = \frac{\Gamma_P + \Gamma_E}{2}$$.


 * $$\frac {d \phi}{dx}$$ is calculated using the nodal point values (Figure 4):

$$\left ( \frac{d\phi }{dx} \right )_e=\frac{\phi_E - \phi_P}{\delta x_{PE}}$$ , $$\left ( \frac{d\phi }{dx} \right )_w=\frac{\phi_P - \phi_W}{\delta x_{WP}}$$,


 * In some practical situations, the source term can be linearized:

$$\overrightarrow{S} \Delta V= S_u + S_p\phi_p$$.

$$\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 ) + (S_u + S_p\phi_p)$$.
 * Merging the above equations leads to


 * Re-arranging gives

$$\left ( \frac{\Gamma_e}{\delta x_{PE}} A_e + \frac{\Gamma_w}{\delta x_{WP}} A_w - S_p \right )\phi_P = \left ( \frac{\Gamma_w}{\delta x_{WP}} A_w \right ) \phi_W + \left ( \frac{\Gamma_e}{\delta x_{WP}} A_e \right ) \phi_E + S_u$$.


 * Compare and identify the above equation with

$$ a_P\phi_P = a_W\phi_W + a_E\phi_E + S_u $$

where

STEP 3:

Solution of equations
 * Discretized equation must be set up at each of the nodal points in order to solve the problem.
 * The resulting system of linear algebraic equations Linear equation can then be solved to obtain $$\phi$$ at the nodal points.
 * The matrix of higher order can be solved in MATLAB.

This method can also be applied to a 2D situation. See Finite volume method for two dimensional diffusion problem.