Discretization of diffusion term

Consider the steady-state diffusion equation given by:

\[-\nabla\cdot(\Gamma^\phi\nabla\phi) = Q^\phi\]

The equation can be written more generally in term of a diffusion flux $J^{\phi,D}$ as

\[\nabla\cdot\mathbf{J}^{\phi,D}= Q^\phi\]

Integral over a cell (one integration point at each face):

\[\sum_{f}(-\Gamma^\phi\nabla\phi)_f\cdot\mathbf{S}_f = Q_C^{\phi}V_C\]

On non-orthogonal grids the surface vector $\mathbf{S}_f$ should be written as the sum of two vecors $\mathbf{E}_f$ and $\mathbf{T}_f$:

\[\mathbf{S}_f=\mathbf{E}_f+\mathbf{T}_f\]

thus,

\[(\nabla\phi)_f\cdot\mathbf{S}_f=(\nabla\phi)_f\cdot\mathbf{E}_f+(\nabla\phi)_f\cdot\mathbf{T}_f=\mathbf{E}_f\frac{\phi_F-\phi_C}{d_CF}+(\nabla\phi)_f\cdot\mathbf{T}_f\]

The first term on the righthand side represents a contribution similar to the contribution on orthogonal grids, while the second term on the right hand sie is called cross-diffusion or non-orthogonal diffusion.

Minimum Correction Approach

\[\mathbf{E}_f=(\mathbf{e}\cdot\mathbf{S}_f)\mathbf{e}=(S_f\cos\theta)\mathbf{e}\]

Orthogonal Correction Approach

\[\mathbf{E}_f=S_f\mathbf{e}\]

Over-Relaxed Approach

\[\mathbf{E}_f=\frac{S_f}{\cos\theta}\mathbf{e}\]

Treatment of the Cross-Diffusion Term

The corss diffusion is treated as a source term in a deferred correction manner by computing its value using the current gradient field.

Gradient Computation

Average gradient over element C:

\[\nabla\phi_C=\frac{1}{V_C}\sum_{f}\phi_f\mathbf{S}_f\]

Gradient at the face f:

\[\nabla\phi_f=g_C\nabla\phi_C+g_F\nabla\phi_F\]

Algebraic Equation for Non-orthogonal Meshes

\[\begin{aligned} \sum_{f}(\mathbf{J}_f^{\phi,D})\cdot\mathbf{S}_f & = \sum_f\Gamma_f^\phi gDiff_f(\phi_C-\phi_F)+ \sum_{f}(-(\Gamma^\phi\nabla\phi)_f\cdot\mathbf{T}_f) \\ & = (\sum_fFluxC_f)\phi_C+\sum_fFluxF_f\phi_F+\sum_f(FluxV_f) \end{aligned}\]

where $gDiff_f=\frac{E_f}{d_{CF}}$ is the gemetric diffusion coefficient.

The final form of discretized diffusion equation is:

\[a_C\phi_C+\sum_F a_F\phi_F=b_C\]

where

\[a_F=FluxF_f=-\Gamma_f^\phi gDiff_f \\ a_c=\sum_f FluxC_f=-\sum_f FluxF_f=\sum_f \Gamma_f^\phi gDiff_f \\ b_C=Q_C^\phi V_C-\sum_f FluxV_f=Q_C^\phi+\sum_f(-(\Gamma^\phi\nabla\phi)_f\cdot\mathbf{T}_f)\]

Boundary conditions

Dirichlet boundary condition

Neumann boundary condition

Mixed bonundary condition

Skewness

skewness correction:

\[\phi_f=\phi_{f'}+(\nabla\phi)_{f'}\cdot\mathbf d_{f'f}\]

Anisotropic diffusion

Under-relaxation of the iterative solution process

\[\phi_C=\phi_C^*+\lambda^\phi\left(\frac{-\sum_F a_F\phi_F+b_C}{a_C}-\phi_C^*\right)\]

thus,

\[\frac{a_C}{\lambda^\phi}\phi_C+\sum_F a_F\phi_F=b_C+\frac{(1-\lambda^\phi)a_C}{\lambda^\phi}\phi_C^*\]