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. For Studying Finite-volume method for unsteady flow there is some governing equations >

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}\left(\rho \phi \upsilon\right) = \operatorname{div}\left(\Gamma \operatorname{grad} \phi\right) + 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}\left(\rho \phi \upsilon\right)$$ is Net rate of flow of $$ \phi $$ out of fluid element(convection),

$$\operatorname{div}\left(\Gamma \operatorname{grad} \phi\right) $$ 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} \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.

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

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

$$\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(k A \frac{\partial T} {\partial x}\right)_e - \left(k A \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

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

$$ \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 $$ \theta $$ between 0 and 1, and we write the integration of $$ T_P $$

$$ 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 $$ \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 $$ Thus\\ $$ \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 $$

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 + \left[ {a_P}^0 - \left( a_w + a_e - S_P \right)\right] {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 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, $$ \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 $$

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

$$ 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 $$ 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 $$ \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 $$ {T_P}^0$$ satisfies the following condition

$$ {a_P}^0 = \left[ \frac {a_E + a_W} {2} \right]$$

which leads to

$$ \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:

$$ a_P T_P = a_W T_W + a_E T_E + {a_P}^0 {T_P}^0 + S_u $$

$$ 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  $$ T^0 $$. The system of equations is solved after selecting time step $$ \Delta t $$. Next the solution $$ T $$ is assigned to $$ 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