Finite volume method for unsteady flow

Unsteady flows are characterized as flows in which the properties of the fluid are time dependent. It gets reflected in the governing equations as the time derivative of the properties are absent.

Governing Equation

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

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

$\rho$ is density and $\phi$ is conservative form of all fluid flow,
$\Gamma$ is the Diffusion coefficient and $S$ is the Source term. $\operatorname {div} (\rho \phi \upsilon )$ is 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),

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.

$\int \limits _{cv}\!\!\!\int _{t}^{t+\Delta t}({\frac {\partial \rho \phi }{\partial t}}\,dt)\,dV+\int _{t}^{t+\Delta t}\!\!\!\int \limits _{A}(n.{\rho \phi u}\,dA)\,dt=\int _{t}^{t+\Delta t}\!\!\!\int \limits _{A}(n.(\Gamma \operatorname {grad} \phi )\,dA)\,dt+\int _{t}^{t+\Delta t}\!\!\!\int \limits _{cv}S_{\phi }\,dV\,dt$

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.

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

$\int _{t}^{t+\Delta t}\!\!\!\int \limits _{cv}\rho c{\frac {\partial T}{\partial t}}\,dV\,dt=\int _{t}^{t+\Delta t}\!\!\!\int \limits _{cv}{\frac {\partial {\frac {k\partial T}{\partial x}}}{\partial x}}\,dV\,dt+\int _{t}^{t+\Delta t}\!\!\!\int \limits _{cv}S\,dV\,dt$

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

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

$\int \limits _{cv}\!\!\!\int _{t}^{t+\Delta t}(\rho c{\frac {\partial T}{\partial t}}\,dt)\,dV=\rho c(T_{P}-{T_{P}}^{O})\Delta V$

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

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

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

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

Now, the exact form of the final discretised equation depends on the value of $\Theta$. As the variance of $\Theta$ is 0< $\Theta$ <1, the scheme to be used to calculate $T_{P}$ depends on the value of the $\Theta$

Different Schemes

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

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

where $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 it’s Taylor series truncation error is first order with respect to time. All coefficients need to be positive. For constant k and uniform grid spacing, $\delta x_{PE}=\delta x_{WP}=\Delta x$ this condition may be written as

$\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 $\Delta x$