Finite volume method for unsteady flow

The conservation equation for the transport of a scalar in unsteady flow has the general form as ^[2]

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

${\displaystyle \rho }$ is density and ${\displaystyle \phi }$ is conservative form of all fluid flow,
${\displaystyle \Gamma }$ is the Diffusion coefficient and ${\displaystyle S}$ is the Source term. ${\displaystyle \operatorname {div} \left(\rho \phi \upsilon \right)}$ is Net rate of flow of ${\displaystyle \phi }$ out of fluid element(convection),
${\displaystyle \operatorname {div} \left(\Gamma \operatorname {grad} \phi \right)}$ is Rate of increase of ${\displaystyle \phi }$ due to diffusion,
${\displaystyle S_{\phi }}$ is Rate of increase of ${\displaystyle \phi }$ due to sources.

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

The first term of the equation reflects the unsteadiness of the flow and is absent in case of steady flows. The finite volume integration of the governing equation is carried out over a control volume and also over a finite time step ∆t.

${\displaystyle \int \limits _{cv}\!\!\!\int _{t}^{t+\Delta t}\left({\frac {\partial \rho \phi }{\partial t}}\,\mathrm {d} t\right)\,\mathrm {d} V+\int _{t}^{t+\Delta t}\!\!\!\int \limits _{A}\left(n.{\rho \phi u}\,\mathrm {d} A\right)\,\mathrm {d} t=\int _{t}^{t+\Delta t}\!\!\!\int \limits _{A}\left(n\cdot \left(\Gamma \operatorname {grad} \phi \right)\,\mathrm {d} A\right)\,\mathrm {d} t+\int _{t}^{t+\Delta t}\!\!\!\int \limits _{cv}S_{\phi }\,\mathrm {d} V\,\mathrm {d} t}$

The control volume integration of the steady part of the equation is similar to the steady state governing equation's integration. We need to focus on the integration of the unsteady component of the equation. To get a feel of the integration technique, we refer to the one-dimensional unsteady heat conduction equation.^[3]

${\displaystyle \rho c{\frac {\partial T}{\partial t}}={\frac {\partial }{\partial x}}{\frac {k\partial T}{\partial x}}+S}$

${\displaystyle \int _{t}^{t+\Delta t}\!\!\!\int \limits _{cv}\rho c{\frac {\partial T}{\partial t}}\,\mathrm {d} V\,\mathrm {d} t=\int _{t}^{t+\Delta t}\!\!\!\int \limits _{cv}{\frac {\partial }{\partial x}}{\frac {k\partial T}{\partial x}}\,\mathrm {d} V\,\mathrm {d} t+\int _{t}^{t+\Delta t}\!\!\!\int \limits _{cv}S\,\mathrm {d} V\,\mathrm {d} t}$

${\displaystyle \int _{e}^{w}\!\!\!\int _{t}^{t+\Delta t}\left(\rho c{\frac {\partial T}{\partial t}}\,\mathrm {d} t\right)\,\mathrm {d} V=\int _{t}^{t+\Delta t}\left[\left(kA{\frac {\partial T}{\partial x}}\right)_{e}-\left(kA{\frac {\partial T}{\partial x}}\right)_{w}\right]\,\mathrm {d} t+\int _{t}^{t+\Delta t}{\bar {S}}\Delta V\,\mathrm {d} t}$

Now, holding the assumption of the temperature at the node being prevalent in the entire control volume, the left side of the equation can be written as ^[4]

${\displaystyle \int \limits _{cv}\!\!\!\int _{t}^{t+\Delta t}\left(\rho c{\frac {\partial T}{\partial t}}\,\mathrm {d} t\right)\,\mathrm {d} V=\rho c\left(T_{P}-{T_{P}}^{O}\right)\Delta V}$

By using a first order backward differencing scheme, we can write the right hand side of the equation as

${\displaystyle \rho c\left(T_{P}-{T_{P}}^{0}\right)\Delta V=\int _{t}^{t+\Delta t}\left[\left(K_{e}A{\frac {T_{E}-T_{P}}{\delta x_{PE}}}\right)-\left(K_{w}A{\frac {T_{P}-T_{W}}{\delta x_{WP}}}\right)\right]\,\mathrm {d} t+\int _{t}^{t+\Delta t}{\bar {S}}\Delta V\,\mathrm {d} t}$

Now to evaluate the right hand side of the equation we use a weighting parameter ${\displaystyle \theta }$ between 0 and 1, and we write the integration of ${\displaystyle T_{P}}$

${\displaystyle I_{T}=\int _{t}^{t+\Delta t}T_{P}\,\mathrm {d} t=\left[\theta T_{P}+\left(1-\theta \right){T_{P}}^{0}\right]\Delta t}$

Now, the exact form of the final discretised equation depends on the value of ${\displaystyle \Theta }$. As the variance of ${\displaystyle \Theta }$ is 0< ${\displaystyle \Theta }$ <1, the scheme to be used to calculate ${\displaystyle T_{P}}$ depends on the value of the ${\displaystyle \Theta }$ Thus\\ ${\displaystyle \rho c{\frac {\left(T_{P}-{T_{P}}^{0}\right)}{\Delta t}}\Delta x=\theta [\left(K_{e}{\frac {T_{E}-T_{P}}{\delta x_{PE}}}\right)-\left(K_{w}{\frac {T_{P}-T_{W}}{\delta x_{WP}}}\right)]+(1-\theta )[\left(K_{e}{\frac {T_{E}-T_{P}}{\delta x_{PE}}}\right)-\left(K_{w}{\frac {T_{P}-T_{W}}{\delta x_{WP}}}\right)]+{\bar {S}}\Delta x}$

1. Explicit Scheme in the explicit scheme the source term is linearised as ${\displaystyle b=S_{u}+{S_{P}}{T_{P}}^{0}}$. We substitute ${\displaystyle \theta =0}$ to get the explicit discretisation i.e.:^[5]

${\displaystyle a_{P}T_{P}=a_{w}{T_{w}}^{0}+a_{e}{T_{e}}^{0}+\left[{a_{P}}^{0}-\left(a_{w}+a_{e}-S_{P}\right)\right]{T_{P}}^{0}+S_{u}}$

where ${\displaystyle a_{P}={a_{P}}^{0}}$. One thing worth noting is that the right side contains values at the old time step and hence the left side can be calculated by forward matching in time. The scheme is based on backward differencing and its Taylor series truncation error is first order with respect to time. All coefficients need to be positive. For constant k and uniform grid spacing, ${\displaystyle \delta x_{PE}=\delta x_{WP}=\Delta x}$ this condition may be written as

${\displaystyle \rho c{\frac {\Delta x}{\Delta t}}>{\frac {2K}{\Delta x}}}$

This inequality sets a stringent condition on the maximum time step that can be used and represents a serious limitation on the scheme. It becomes very expensive to improve the spatial accuracy because the maximum possible time step needs to be reduced as the square of ${\displaystyle \Delta x}$ ^[6]

2. Crank-Nicolson scheme : the Crank-Nicolson method results from setting ${\displaystyle \theta ={\frac {1}{2}}}$. The discretised unsteady heat conduction equation becomes

${\displaystyle a_{P}T_{P}=a_{E}\left[{\frac {T_{E}+{T_{E}}^{0}}{2}}\right]+a_{W}\left[{\frac {T_{W}+{T_{W}}^{0}}{2}}\right]+\left[{a_{P}}^{0}-{\frac {a_{E}}{2}}-{\frac {a_{W}}{2}}\right]{T_{P}}^{0}+b}$

Where ${\displaystyle a_{P}={\frac {a_{W}+a_{E}}{2}}+{a_{P}}^{0}-{\frac {S_{P}}{2}}}$

Since more than one unknown value of T at the new time level is present in equation the method is implicit and simultaneous equations for all node points need to be solved at each time step. Although schemes with ${\displaystyle {\frac {1}{2}}<\theta <1}$ including the Crank-Nicolson scheme, are unconditionally stable for all values of the time step it is more important to ensure that all coefficients are positive for physically realistic and bounded results. This is the case if the coefficient of ${\displaystyle {T_{P}}^{0}}$ satisfies the following condition

${\displaystyle {a_{P}}^{0}=\left[{\frac {a_{E}+a_{W}}{2}}\right]}$

which leads to

${\displaystyle \Delta t<\rho c{\frac {\Delta x^{2}}{K}}}$

the Crank-Nicolson is based on central differencing and hence is second order accurate in time. The overall accuracy of a computation depends also on the spatial differencing practice, so the Crank-Nicolson scheme is normally used in conjunction with spatial central differencing

3. Fully implicit Scheme when the value of Ѳ is set to 1 we get the fully implicit scheme. The discretised equation is: ^[7]

${\displaystyle a_{P}T_{P}=a_{W}T_{W}+a_{E}T_{E}+{a_{P}}^{0}{T_{P}}^{0}+S_{u}}$

${\displaystyle a_{P}={a_{P}}^{0}+a_{W}+a_{E}-S_{P}}$

Both sides of the equation contain temperatures at the new time step, and a system of algebraic equations must be solved at each time level. The time marching procedure starts with a given initial field of temperatures ${\displaystyle T^{0}}$. The system of equations is solved after selecting time step ${\displaystyle \Delta t}$. Next the solution ${\displaystyle T}$ is assigned to ${\displaystyle T^{0}}$ and the procedure is repeated to progress the solution by a further time step. It can be seen that all coefficients are positive, which makes the implicit scheme unconditionally stable for any size of time step. Since the accuracy of the scheme is only first-order in time, small time steps are needed to ensure the accuracy of results. The implicit method is recommended for general purpose transient calculations because of its robustness and unconditional stability