as I learned from the awesome book^{1}

If properly handled, source term improves solution stability. Consider a discretised conservation equation for a general variable \phi with source in element C written below:

\begin{equation} a_{C}\bold \phi_{C} + \sum_{F} a_{F}\phi_{F} = Q_{C}^{\phi}V_{C} \end{equation}

If the RHS is large compared to the rest of the terms, the rate of convergence will be affected. In such cases, the rate of convergence can be improved by source term linearisation.

**Source term linearisation**

Here, Q_{C} is expanded using Taylor series up to first order.

\begin{equation} Q_{C}^{\phi} = Q_{C}^{\phi^{*}}+\Big(\frac{\partial Q^{\phi} }{\partial \phi}\Big)_{C}^{*} (\phi_{C}-\phi_{C}^{*}) \end{equation}

Hence the RHS of (1) can be written as

\begin{split} Q_{C}^{\phi}V_{C} &= Q_{C}^{\phi^{*}}V_{C}+\Big(\frac{\partial Q^{\phi} }{\partial \phi}\Big)_{C}^{*} (\phi_{C}-\phi_{C}^{*})V_{C}\\ &= (\frac{\partial Q^{\phi} }{\partial \phi}\Big)_{C}^{*} V_{C}\phi_{C} +\Big[ Q_{C}^{\phi^{*}} - (\frac{\partial Q^{\phi} }{\partial \phi}\Big)_{C}^{*} \phi_{C}^{*} \Big]V_{C}\\ &= FluxC_{C} \phi_{C} + FluxV_{C} \end{split}

First term in the RHS is treated as implicit and second as explicit. Substituting for RHS in (1)

\begin{equation} [a_{C}-FluxC_{C}]\bold \phi_{C} + \sum_{F} a_{F}\phi_{F} = FluxV_{C} \end{equation}

FluxC_{C}should be -ve, else the system will loose diagonal dominance and Scarborough criterion shall not be fulfilled. This shall result in divergence. Note that, still a_{C} = -(a_{E}+a_{W}+a_{S}+a_{N}) pg # 215 , Moukalled et. al

Consider the case of radiation heat transfer. Q_{C} = A (T_{\inf}^{4} - T_{C}^{4}). If we take FluxC_{C}= 0 and FluxV_{C} = A (T_{\inf}^{4} - T_{C}^{4}) , for T_{C} > T_{\inf}, FluxV_{C} will be negative which may result in negative T_{C} value during iteration. If we employ source term linearisation, then FluxC_{C} = -4AT_{C}^{3} V_{C} and FluxV_{C} = A (T_{\inf}^{4} + 3T_{C}^{4}) V_{C} which is perfect because of negative FluxC_{C} and positive FluxV_{C}.

**Patankar’s implicit relaxation**

Consider a conservation equation of the following form

\begin{equation} a_{C}\bold \phi_{C} + \sum_{F} a_{F}\phi_{F} = b_{C} \end{equation}

Adding an subtracting previous value of variable,

\begin{split} \phi_{C} &= \frac{b_{C} - \sum a_{F}\phi_{F}}{a_{C}}\\ &=\phi_{C}^{*} + \frac{b_{C} - \sum a_{F}\phi_{F}}{a_{C}} - \phi_{C}^{*} \\ \end{split}

In order to relax, we apply \lambda to the second term in RHS , (which is the difference with previous value)

\phi_{C} =\phi_{C}^{*} + \lambda \Big(\frac{b_{C} - \sum a_{F}\phi_{F}}{a_{C}} - \phi_{C}^{*} \Big)\\

rearranging,

\frac{a_{C}}{\lambda}\phi_{C}+ \sum a_{F}\phi_{F} = b_{C}+ a_{C}( \frac{1}{\lambda}-1)\phi_{C}^{*}

which can be written in the original form as (4) by reconsidering,

\begin{split} \frac{a_{C}}{\lambda} &\rightarrow a_{C}\\ b_{C}+ a_{C}( \frac{1}{\lambda}-1)\phi_{C}^{*} &\rightarrow b_{C} \end{split}

It is clear that this method improves the diagonal dominance if \lambda <1

**E-Factor method**

In this method (4) is written in the following way:

\begin{split} a_{C}\phi_{C} &= b_{C} - \sum a_{F}\phi_{F}\\ &=\lambda( b_{C} - \sum a_{F}\phi_{F})+ (1-\lambda)a_{C}\phi_{C}^{*}\\ &=\frac{E}{1+E}( b_{C} - \sum a_{F}\phi_{F})+ \frac{1}{1+E}a_{C}\phi_{C}^{*} \text{, replacing }\lambda \text{ by} \frac{E}{1+E} \\ \end{split}

a_{C}\phi_{C}(1+\frac{1}{E}) +\sum a_{F}\phi_{F} =b_{C} +\frac{1}{E} a_{C}\phi_{C}^{*}

This under-relaxation effect can be interpreted in terms of some artificial transient time scale that advances \phi_{C} at each iteration. The time step \Delta t is proportional to \Delta t^{*} which is the characteristic time scale defined as

\Delta t^{*} = \frac{\rho_{C}V_{C}}{a_{C}}

Characteristic time scale is related to the time needed to diffuse and convect a change of \phi_{C} across the element. Hence E is like element CFL number. The following points rae to be noted:

- The \Delta t^{*} is proportional to element volume.
- Solution in smaller element advance slower compared to a bigger element.
- This is also a characteristic of Patankar’s implicit relaxation.
- E and \lambda are related by E = \frac{1}{1-\lambda}
- General values E: 4 – 10, \lambda : 0.75 – 0.9.

**False transient relaxation**

This can be viewed as a modification of Euler’s first order implicit transient method. Here, diagonal dominance id enhanced by adding pseudo transient term to diagonal coefficient and a pseudo previous time term to RHS for balance.Equation (4) is return as following:

\begin{equation} (a_{C}+a_{C_{0}})\bold \phi_{C} + \sum_{F} a_{F}\phi_{F} = b_{C}+a_{C_{0}}\phi_{C} ^{*}\end{equation}

where,

a_{C_{0}} = \frac{\rho_{C}V_{C}}{\Delta t}

where, \Delta t is a user supplied time step.

The following points shall be noted for under relaxation in general

- Different equations have Different URF
- Different parts of the domain may employ Different URF
- Different iteration may have Different URF value.