Finite Volume Method for the Compressible Viscous Gas Flows. SIMPLE-TS Derivation

This information is useful for people interested in derivation and understanding of finite volume method for calculation of unsteady, viscous, compressible, and heat-conductive flows at all speeds. Files with the symbolical derivation of numerical equations, examples, and C++ parallel code are also available. All considerations are straightforwardly applicable to unsteady/steady, viscous, incompressible, and isothermal flows.

Please cite my papers if you find the information useful:

K. Shterev and S. Stefanov, Pressure based finite volume method for calculation of compressible viscous gas flows, Journal of Computational Physics 229 (2010) pp. 461-480, doi:10.1016/j.jcp.2009.09.042 – the accepted manuscript can be downloaded from here, the paper in it`s final mode is available here. IF 3.023

K. S. Shterev and S. K. Stefanov, A Parallelization of Finite Volume Method for Calculation of Gas Microflows by Domain Decomposition Methods, 7th International Conference-Large-ScaleScientific Computations, Sozopol, Bulgaria, June 04-08,2009. Lecture Notes in Computer Science Volume 5910, 2010, DOI:10.1007/978 -3-642-12535 -5, SJR 0.295.

Kiril S. Shterev, GPU implementation of algorithm SIMPLE-TS for calculation of unsteady, viscous, compressible and heat-conductive gas flows, URL: https://arxiv.org/abs/1802.04243, 2018.

Here is placed derivation of the numerical algorithm SIMPLE-TS (Semi Implicit Method for Pressure Linked Equations – Time Step). SIMPLE-TS is a finite volume method for the calculation of unsteady, viscous, compressible, and heat-conductive flows. It can be reduced straightforwardly to calculate unsteady, viscous, incompressible, and isothermal flows. SIMPLE-TS is part of SIMPLE-like algorithms. A basic idea in the algorithm SIMPLE-TS is to keep the equation of state maximally satisfied within the entire iterative process. The main contribution is the substitution of density in the unsteady term in the pressure equation with pressure and temperature. This way we turned the numerical pressure equation into a stable that does not need a relaxation coefficient to ensure convergence. Furthermore, SIMPLE-TS is five times faster than SIMPLE and slightly faster than PISO.

Resources (available downloads) on this page:
Derivation and rearrangement of terms of equations are complicated and could be a source of errors. To reduce mistakes, and make derivation faster, and easier we advise you to use appropriate software to do symbolic calculations. You can download the derivation of the numerical equations of the algorithm SIMPLE-TS using a first-order upwind scheme for the approximation of convective terms with Mathematica from here:
SIMPLE_TS_2D_staggered_mesh_First_order_upwind_scheme_approximation_of_convective_terms.nb
SIMPLE_TS_2D_staggered_mesh_First_order_upwind_scheme_approximation_of_convective_terms.pdf

A second order total variation diminishing (TVD) scheme can be used to approximate convective terms. A good explanation of implementation of TVD schemes in finite volume methods is textbook of Versteeg and Malalasekera, Chapter 5 The finite volume method for convection-diffusion problems, . You can download the derivation of the numerical equations of the algorithm SIMPLE-TS using second order TVD scheme for the approximation of convective terms with Mathematica from here:
SIMPLE_TS_2D_staggered_mesh_Second_order_TVD_schemes_approximation_of_convective_terms.nb
SIMPLE_TS_2D_staggered_mesh_Second_order_TVD_schemes_approximation_of_convective_terms.pdf
You can download the program that calculates and plots results of the first time step of one-dimensional unsteady, isothermal pressure-driven flow in a duct from here: One_dimensional_flow_in_a_duct.m

Resources (available downloads) on the example page are:

parallel C++ source code with a first-order approximation of convective terms;
unsteady supersonic, compressible, viscous, heat-conductive fluid flow past a confined square in a micro-channel – Mach number 2.43 and Knudsen number 0.00283 (Reynolds number 1415);
unsteady subsonic, compressible, viscous, heat-conductive fluid flow past a confined square in a micro-channel – Mach number 0.1 and Knudsen number 0.00194 (Reynolds number 85);
increasing velocity at the channel inflow from Mach number 2.43 to Mach number 4.86, for Knudsen number 0.05;
Rayleigh-Bènard flow of a rarefied gas;

Derivation and Considerations

SIMPLE-TS calculates numerically Navier-Stokes equation. SIMPLE-TS can be implemented to calculate unsteady/steady, viscous, compressible, and head-conductive fluid and unsteady/steady, viscous, incompressible fluid. Seeking more generality we consider a unified two-dimensional system of equations describing the time evolution of a unsteady, viscous, compressible, heat-conducting gas flow.

(1) $\begin{equation*}\frac{\partial \rho }{\partial t} + \frac{\partial (\rho u)}{\partial x} + \frac{\partial (\rho v)}{\partial y} = 0\end{equation*}$

(2) $\begin{equation*}\begin{split}\frac{\partial(\rho u)}{\partial t} + \frac{\partial(\rho u u)}{\partial x} + \frac{\partial(\rho v u)}{\partial y} =& \rho g_x - A \frac{\partial p}{\partial x} + B\left[\frac{\partial}{\partial x}\left(\Gamma\frac{\partial u}{\partial x}\right) + \frac{\partial}{\partial y}\left(\Gamma \frac{\partial u}{\partial y}\right)\right] \\& + B \left\{\frac{\partial}{\partial x}\left(\Gamma \frac{\partial u}{\partial x}\right) + \frac{\partial}{\partial y}\left(\Gamma \frac{\partial v}{\partial x}\right) - \frac{2}{3} \frac{\partial}{\partial x}\left[\Gamma\left(\frac{\partial u}{\partial x} + \frac{\partial v}{\partial y}\right)\right]\right\}\end{split}\end{equation*}$

(3) $\begin{equation*}\begin{split}\frac{\partial(\rho v)}{\partial t} + \frac{\partial(\rho u v)}{\partial x} + \frac{\partial(\rho v v)}{\partial y} =& \rho g_y - A \frac{\partial p}{\partial y} + B\left[\frac{\partial}{\partial x}\left(\Gamma \frac{\partial v}{\partial x}\right) + \frac{\partial}{\partial y}\left(\Gamma \frac{\partial v}{\partial y}\right)\right] \\& + B \left\{\frac{\partial}{\partial y}\left(\Gamma \frac{\partial v}{\partial y}\right) + \frac{\partial}{\partial x}\left(\Gamma \frac{\partial u}{\partial y}\right) - \frac{2}{3} \frac{\partial}{\partial y}\left[\Gamma\left(\frac{\partial u}{\partial x} + \frac{\partial v}{\partial y}\right)\right]\right\}\end{split}\end{equation*}$

(4) $\begin{equation*}\frac{\partial(\rho T)}{\partial t} + \frac{\partial(\rho u T)}{\partial x} + \frac{\partial(\rho v T)}{\partial y} = C^{T1}\left[\frac{\partial}{\partial x}\left(\Gamma^{\lambda}\frac{\partial T}{\partial x}\right) + \frac{\partial}{\partial y}\left(\Gamma^{\lambda}\frac{\partial T}{\partial y}\right)\right] + C^{T2} \Gamma \Phi + C^{T3}\frac{Dp}{Dt}\end{equation*}$

(5) $\begin{equation*}p=\rho T\end{equation*}$

where:

(6) $\begin{equation*}\Phi=2\left[\left(\frac{\partial u}{\partial x}\right)^2 + \left(\frac{\partial v}{\partial y}\right)^2\right] + \left(\frac{\partial v}{\partial x} + \frac{\partial u}{\partial y}\right)^2 - \frac{2}{3}\left(\frac{\partial u}{\partial x} + \frac{\partial v}{\partial y}\right)^2\end{equation*}$

$u$ is the horizontal component of velocity, $v$ is the vertical component of velocity, $p$ is pressure, $T$ is temperature, $\rho$ is density, $t$ is time, $x$ and $y$ are coordinates of a Cartesian coordinate system. The reference quantities used to scale the system of equations (1) – (5) will be defined later, separately for each considered problem. Parameters $A$ , $B$ , $g_x$ , $g_y$ , $C^{T1}$ , $C^{T2}$ , $C^{T3}$ and diffusion coefficients $\Gamma$ and $\Gamma^\lambda$ , given in Eqs. (1) – (5), depend on the gas model and the equation non-dimensional form.

The discretization of the equation system is accomplished by using a backward staggered velocity grid, in which all dependent field variables (pressure, temperature, and density) are calculated at a cell center, and all flow variables (velocity components) are calculated at the surfaces of a cell (Fig. 1). The variable’s value at a defined point is considered an average value of the variable in the control volume. Using this basic assumption, we interpreted and converted some of the terms in the equations to a more suitable form to apply the finite volume method. As a result, it reduces the number of middle points and their interpolations, increasing the numerical scheme’s stability and accuracy. The important test for such stability and accuracy improvement is mesh convergence. When the scheme contains not sufficiently well-implemented interpolations in middle points, you will not be able to use small spacial steps. At the same time, the coarser mesh will obtain correct results that correspond to its accuracy. Note that “not sufficiently well” is hard to determine precisely, and such mesh convergence tests are time-consuming.
In general, when you use some approach, as in our case finite volume method, you have to follow the basic assumptions of the approach. Messing it partly with other approaches such as finite difference or finite element methods could lead to erroneous results. For example, the finite volume method integrates over a control volume. As a result numerical scheme is multiplied by a volume. On the other hand, finite difference uses Tailor series, and spacial steps appear in a denominator. So, if you use a mesh where both approaches are easily applicable as uniform Cartesian mesh, numerical expressions obtained for each scheme will be different. Most likely, in this case, you have to multiply finite difference numerical expression by control volume to obtain the same expression as the finite volume method. Here we followed the basic assumption of finite volume method, as they are presented in Patankar’s textbook .

For interpolation between two neighbour points a piecewise-linear profile is used. In that way we approximate the derivations of second order . For the interpolation of convective terms a first order upwind scheme is used.In Fig. 1 the following grid variables are denoted:

$x^{f}_{i}$ – the left frontier x-coordinate of the control volume of node $i$
$y^{f}_{j}$ – bottom frontier y-coordinate of the control volume of node $j$
$\Delta x_i = x^{f}_{i+1} - x^{f}_{i}$ – step on OX
$\Delta y_j = y^{f}_{j+1} - y^{f}_{j}$ – step on OY
$x^{v}_{i}$ – x-coordinate of the centre of the control volume of node $i$ , $x^{v}_{i} = x^{f}_{i} + 0.5\Delta x_i$
$y^{v}_{j}$ – y-coordinate of the centre of the control volume of node $j$ , $y^{v}_{j} = y^{f}_{j} + 0.5\Delta y_j$
$\phi_{i,j}$ – field variables defined at point $(x^{v}_{i}, y^{v}_{j})$
$u_{i,j}$ – the horizontal component of velocity $(x^{f}_{i}, y^{v}_{j})$
$v_{i,j}$ – the vertical component of velocity defined at point $(x^{v}_{i}, y^{f}_{j})$
$F^{x}_{i,j}$ – the convective mass flux through the surface between control volumes $(i-1,j)$ and $(i,j)$ (in horizontal direction)

(7) $\begin{equation*} F^{x}_{i,j} = \rho^{u}_{i,j} u_{i,j} \Delta y_j \end{equation*}$

$F^{y}_{i,j}$ – the convective mass flux through the surface between control volumes $(i,j-1)$ and $(i,j)$ (in vertical direction)

(8) $\begin{equation*} F^{y}_{i,j} = \rho^{v}_{i,j} v_{i,j} \Delta x_i, \end{equation*}$

where:

(9) $\begin{equation*} \rho^{u}_{i,j} = \begin{cases} 0.5(\rho_{i-1,j} + \rho_{i,j})& \text{if $u_{i,j} = 0$},\\ \rho_{i-1,j}& \text{if $u_{i,j} > 0$},\\ \rho_{i,j}& \text{if $u_{i,j} < 0$}. \end{cases} \end{equation*}$

(10) $\begin{equation*} \rho^{v}_{i,j} = \begin{cases} 0.5(\rho_{i,j-1} + \rho_{i,j})& \text{if $v_{i,j} = 0$},\\ \rho_{i,j-1}& \text{if $v_{i,j} > 0$},\\ \rho_{i,j}& \text{if $v_{i,j} < 0$}. \end{cases} \end{equation*}$

Density $(\rho)$ is computed at the middle points, i.e. $\rho^{u}_{i,j}$ at point $(x^{f}_{i}, y^{v}_{j})$ and $\rho^{v}_{i,j}$ at point $(x^{v}_{i}, y^{f}_{j})$ , by using upwind scheme and . The upwind calculation of density is of first order accuracy, but it brings more stability of the calculations, when a supersonic fluid flow is considered.

Finite volume discretization of the momentum equations

For brevity, here a finite volume discretization of the momentum equations is given for the y-momentum (3) only. The equation is written in the following convenient form:

(11) $\begin{equation*} \begin{split} \underbrace{\frac{\partial (\rho v)}{\partial t}}_{Unsteady\ term} + \underbrace{\frac{\partial (\rho u v)}{\partial x} + \frac{\partial (\rho v v)}{\partial y}}_{Convective\ terms} =& - \underbrace{A \frac{\partial p}{\partial y}}_{Pressure\ term} + \underbrace{B\left[\frac{\partial}{\partial x}\left(\Gamma \frac{\partial v}{\partial x}\right) + \frac{\partial}{\partial y}\left(\Gamma \frac{\partial v}{\partial y}\right)\right]}_{Diffusion\ terms}\\ & + \underbrace{\rho g_y + B \left\{\frac{\partial}{\partial y}\left(\Gamma \frac{\partial v}{\partial y}\right) + \frac{\partial}{\partial x}\left(\Gamma \frac{\partial u}{\partial y}\right) - \frac{2}{3} \frac{\partial}{\partial y}\left[\Gamma\left(\frac{\partial u}{\partial x} + \frac{\partial v}{\partial y}\right)\right]\right\}}_{Source\ term\ (S^v)} \end{split}\end{equation*}$

Integration of convective terms

The integration of (11) is accomplished for $v$ over a control volume (Fig. 1), as given in : integrating from $x^f_i$ to $x^f_{i+1}$ for OX, from $y^v_{j-1}$ to $y^f_j$ for OY and integrating over the volume from $x^f_i$ to $x^f_{i+1}$ for OX and from $y^f_j$ to $y^v_j$ for OY, see Fig. 2.1. Two volumes are integrated and the obtained expressions are summed. In the middle nodes, where $v$ is not defined, we have used a upwind interpolation scheme, instead the approach proposed in .

**Figure 2.1.** Vertical velocity control volume

We used upwind scheme when integrated convective terms over the volumes:

(12) $\begin{equation*} \begin{split} &\int_{x^f_i}^{x^f_{i+1}}\!\!\!\int_{y^f_j}^{y^v_{j}} \frac{\partial (\rho v v)}{\partial y} dy dx = \int_{x^f_i}^{x^f_{i+1}} \left[ (\rho v v)_{y^v_j} - (\rho v v)_{y^f_j} \right] dx \\ &= \left[ (\rho v v)_{x^f_i \text{ to } x^f_{i+1}, y^v_j} - (\rho v v)_{x^f_i \text{ to } x^f_{i+1}, y^f_j} \right] \Delta x_i \\ &\approx \left[ max\left(0, \overline{F}^y_{i,j} \right) v_{i,j} - max\left(0, -\overline{F}^y_{i,j} \right) v_{i,j+1} \right] - F^y_{i,j} v_{i,j} , \end{split} \end{equation*}$

where $\overline{F}^y_{i,j}$ is a vertical flux defined at a point $(x^v_i, y^v_j)$ . The standard definition is: $\overline{F}^y_{i,j} = 0.5 (F^y_{i,j} + F^y_{i,j+1})$ . In according to Finite Volume Mothod $\overline{F}^y_{i,j}$ can be defined as: $\overline{F}^y_{i,j} = 0.5 \rho_{i,j} (v_{i,j} + v_{i,j+1})$ , because $\rho$ is defined at point $(x^v_i, y^v_j)$ and doesn’t need interpolation that will impove stability and accuracy of the numerical scheme.

(13) $\begin{equation*} \begin{split} &\int_{x^f_i}^{x^f_{i+1}}\!\!\!\int_{y^v_{j-1}}^{y^f_j} \frac{\partial (\rho v v)}{\partial y} dy dx = \int_{x^f_i}^{x^f_{i+1}} \left[ (\rho v v)_{y^f_j} - (\rho v v)_{y^v_{j-1}} \right] dx \\ &\approx F^y_{i,j} v_{i,j} - \left[ max\left(0, \overline{F}^y_{i,{j-1}} \right) v_{i,{j-1}} - max\left(0, -\overline{F}^y_{i,{j-1}} \right) v_{i,j} \right]. \end{split} \end{equation*}$

Integration of convective term $\dfrac{\partial (\rho u v)}{\partial x}$ :

(14) $\begin{equation*} \begin{split} &\int_{y^f_j}^{y^v_j}\!\!\!\int_{x^f_i}^{x^f_{i+1}} \frac{\partial (\rho u v)}{\partial x} dx dy = \int_{y^f_j}^{y^v_j} \left[ (\rho u v)_{x^f_{i+1}} - (\rho u v)_{x^f_i} \right] dy \\ &= \left[ (\rho u v)_{x^f_{i+1}, y^v_j \text{ to } y^f_j} - (\rho u v)_{x^f_i, y^v_j \text{ to } y^f_j} \right] \frac{\Delta y_j}{2} \\ &\approx \frac{1}{2} \left[ max\left(0, F^x_{{i+1},j} \right) v_{i,j} - max\left(0, -F^x_{{i+1},j} \right) v_{{i+1},j} \right] \\ &\ \ \ \ - \frac{1}{2} \left[ max\left(0, F^x_{i,j} \right) v_{{i-1},j} - max\left(0, -F^x_{i,j} \right) v_{i,j} \right] \end{split} \end{equation*}$

(15) $\begin{equation*} \begin{split} &\int_{y^v_{j-1}}^{y^f_j}\!\!\!\int_{x^f_i}^{x^f_{i+1}} \frac{\partial (\rho u v)}{\partial x} dx dy = \int_{y^v_{j-1}}^{y^f_j} \left[ (\rho u v)_{x^f_{i+1}} - (\rho u v)_{x^f_i} \right] dy \\ &\approx \frac{1}{2} \left[ max\left(0, F^x_{{i+1},{j-1}} \right) v_{i,j} - max\left(0, -F^x_{{i+1},{j-1}} \right) v_{{i+1},j} \right] \\ &\ \ \ \ - \frac{1}{2} \left[ max\left(0, F^x_{i,{j-1}} \right) v_{{i-1},j} - max\left(0, -F^x_{i,{j-1}} \right) v_{i,j} \right] \end{split} \end{equation*}$

Integration of convective terms over the control volume for $v$ Fig. 2.1 is:

(16) $\begin{equation*} \begin{split} &\int_{x^f_i}^{x^f_{i+1}}\!\!\!\int_{y^f_{j-1}}^{y^f_j} \left[ \frac{\partial (\rho u v)}{\partial x} - \frac{\partial (\rho v v)}{\partial y} \right] dy dx \\ &= \int_{x^f_i}^{x^f_{i+1}}\!\!\!\int_{y^f_j}^{y^v_{j}} \frac{\partial (\rho v v)}{\partial y} dy dx + \int_{x^f_i}^{x^f_{i+1}}\!\!\!\int_{y^v_{j-1}}^{y^f_j} \frac{\partial (\rho v v)}{\partial y} dy dx \\ &\ \ \ + \int_{y^f_j}^{y^v_j}\!\!\!\int_{x^f_i}^{x^f_{i+1}} \frac{\partial (\rho u v)}{\partial x} dy dx + \int_{y^v_{j-1}}^{y^f_j}\!\!\!\int_{x^f_i}^{x^f_{i+1}} \frac{\partial (\rho u v)}{\partial x} dy dx \\ &\approx \left[ \begin{array}{l} max\left(0, \overline{F}^y_{i,j} \right) + max\left(0, -\overline{F}^y_{i,{j-1}} \right) \\ + \frac{1}{2} max\left(0, F^x_{{i+1},j} \right) + \frac{1}{2} max\left(0, -F^x_{i,j} \right) \\ + \frac{1}{2} max\left(0, F^x_{{i+1},{j-1}} \right) + \frac{1}{2} max\left(0, -F^x_{i,{j-1}} \right) \end{array} \right] v_{i,j} \\ &\ \ \ - \left[ \begin{array}{l} max\left(0, -\overline{F}^y_{i,j} \right) v_{i,j+1} + max\left(0, \overline{F}^y_{i,{j-1}} \right) v_{i,{j-1}} \\ + \frac{1}{2} \left[ max\left(0, -F^x_{{i+1},j} \right) + max\left(0, -F^x_{{i+1},{j-1}} \right) \right] v_{{i+1},j} \\ + \frac{1}{2} \left[ max\left(0, F^x_{i,{j-1}} \right) + max\left(0, F^x_{i,{j}} \right) \right] v_{{i-1},j} \end{array} \right] \\ &= a^{vc}_0 v_{i,j} - (a^{vc}_1 v_{i-1,j} + a^{vc}_2 v_{i+1,j} + a^{vc}_3 v_{i,j-1} + a^{vc}_4 v_{i,j+1}) \end{split} \end{equation*}$

where the superscripts $vc$ mean vertical velocity ( $v$ ) and convection ( $c$ ), subscripts 0, 1, 2, 3, and 4 denote coefficients of velocites $v_{i,j}$ , $v_{i-1,j}$ , $v_{i+1,j}$ , $v_{i,j-1}$ , and $v_{i,j+1}$ , respectively. Derived coefficients are eqaul to:

(17) $\begin{equation*} \begin{split} a^{vc}_0 =& max\left(0, \overline{F}^y_{i,j} \right) + max\left(0, -\overline{F}^y_{i,{j-1}} \right) \\ & + \frac{1}{2} max\left(0, F^x_{{i+1},j} \right) + \frac{1}{2} max\left(0, -F^x_{i,j} \right) \\ & + \frac{1}{2} max\left(0, F^x_{{i+1},{j-1}} \right) + \frac{1}{2} max\left(0, -F^x_{i,{j-1}} \right)\\ =& a^{vc}_1 + a^{vc}_2 + a^{vc}_3 + a^{vc}_4 + \frac{1}{2} (F^x_{{i+1},{j}} - F^x_{{i},{j}} + F^x_{{i+1},{j-1}} - F^x_{{i},{j-1}}) + \overline{F}^y_{{i},{j}} - \overline{F}^y_{{i},{j-1}} \\ a^{vc}_1 =& \frac{1}{2} \left[ max\left(0, F^x_{i,j} \right) + max\left(0, F^x_{i,{j-1}} \right) \right], \ \ \ \ \ \ \ \ \ \ a^{vc}_3 = max\left(0, \overline{F}^y_{i,{j-1}} \right) \\ a^{vc}_2 =& \frac{1}{2} \left[ max\left(0, -F^x_{{i+1},j} \right) + max\left(0, -F^x_{{i+1},{j-1}} \right) \right], a^{vc}_4 = max\left(0, -\overline{F}^y_{i,j} \right) \end{split} \end{equation*}$

Integration of diffusion terms

We presented the integration of the diffusion terms of the equation (11) over $v$ control volume, Fig. 2.1, term by term. We used linear approximation:

(18) $\begin{equation*} \begin{split} &\int_{y^v_{j-1}}^{y^v_j}\!\!\int_{x^f_i}^{x^f_{i+1}} \frac{\partial}{\partial x} \left(\Gamma \frac{\partial v}{\partial x} \right) dx dy = \int_{y^v_{j-1}}^{y^v_j}\left[ \left(\Gamma \frac{\partial v}{\partial x} \right)_{x^f_{i+1}} - \left(\Gamma \frac{\partial v}{\partial x} \right)_{x^f_{i}} \right] dy \\ &= \left[\Gamma\rvert_{x^f_{i+1}} \left(\frac{\partial v}{\partial x} \right)_{x^f_{i+1}} - \Gamma\rvert_{x^f_{i}} \left(\frac{\partial v}{\partial x} \right)_{x^f_{i}} \right] \frac{\Delta y_j + \Delta y_{j-1}}{2} \\ &\approx \left[\Gamma\rvert_{x^f_{i+1}} \left(\frac{v_{i+1,j} - v_{i,j}}{\frac{1}{2} (\Delta x_{i+1} + \Delta x_{i})} \right) - \Gamma\rvert_{x^f_{i}} \left(\frac{v_{i,j} - v_{i-1,j}}{\frac{1}{2} (\Delta x_{i} + \Delta x_{i-1})} \right) \right] \frac{\Delta y_j + \Delta y_{j-1}}{2} \\ &= \left[\Gamma\rvert_{x^f_{i+1}} \left(\frac{v_{i+1,j} - v_{i,j}}{\Delta x_{i+1} + \Delta x_{i}}\right) - \Gamma\rvert_{x^f_{i}} \left(\frac{v_{i,j} - v_{i-1,j}}{\Delta x_{i} + \Delta x_{i-1}}\right) \right] (\Delta y_j + \Delta y_{j-1}) \end{split} \end{equation*}$

where $\Gamma\rvert_{x^f_{i}}$ and $\Gamma\rvert_{x^f_{i+1}}$ are diffusion coefficients on the surfaces of the control volume $x^f_{i}$ and $x^f_{i+1}$ , respectively. The diffusion $\Gamma$ is defined at nodes, where are defined scalar variables, ( $x^v_i,y^v_j$ ). We assume that $\Gamma$ varies continuously between neighbour control volumes and we use billinear interpolation of the diffusion coefficients on the control volume surfaces to calculate $\Gamma$ values between nodes (in middle nodes) with coordinates $(x^f_i, y^f_j)$ for $\Gamma\rvert_{x^f_{i}}$ . The bilinear interpolation if between for neighbour nodes $(x^v_{i-1},y^v_{j-1})$ , $(x^v_{i},y^v_{j-1})$ , $(x^v_{i},y^v_{j-1})$ , and $(x^v_{i},y^v_{j})$ .

(19) $\begin{equation*} \begin{split} &\int_{x^f_i}^{x^f_{i+1}}\!\!\int_{y^v_{j-1}}^{y^v_j} \frac{\partial}{\partial y} \left(\Gamma \frac{\partial v}{\partial y} \right) dy dx = \int_{x^f_i}^{x^f_{i+1}} \left[\left(\Gamma \frac{\partial v}{\partial y} \right)_{y^v_j} - \left(\Gamma \frac{\partial v}{\partial y} \right)_{y^v_{j-1}}\right] dx \\ &= \left[\Gamma\rvert_{y^v_j} \left(\frac{\partial v}{\partial y} \right)_{y^v_j} - \Gamma\rvert_{y^v_{j-1}} \left(\frac{\partial v}{\partial y} \right)_{y^v_{j-1}}\right] \Delta x_i \\ &\approx \left[\Gamma_{i,j} \left(\frac{v_{i,j+1} - v_{i,j}}{\Delta y_j} \right) - \Gamma_{i,j-1} \left(\frac{v_{i,j} - v_{i,j-1}}{\Delta y_{j-1}}\right)\right] \Delta x_i \end{split} \end{equation*}$

Put all together of integrated diffusion terms over $v$ control volume, Fig. 2.1, using derived expressions:

(20) $\begin{equation*} \begin{split} &\int_{x^f_i}^{x^f_{i+1}}\!\!\int_{y^v_{j-1}}^{y^v_j} B \left[ \dfrac{\partial}{\partial x} \left(\Gamma \dfrac{\partial v}{\partial x} \right) + \dfrac{\partial}{\partial y} \left(\Gamma \dfrac{\partial v}{\partial y} \right) \right] dy dx \\ &\approx B \left[\Gamma\rvert_{x^f_{i+1}} \left(\frac{v_{i+1,j} - v_{i,j}}{\Delta x_{i+1} + \Delta x_{i}}\right) - \Gamma\rvert_{x^f_{i}} \left(\frac{v_{i,j} - v_{i-1,j}}{\Delta x_{i} + \Delta x_{i-1}}\right) \right] (\Delta y_j + \Delta y_{j-1}) \\ &\ \ \ \ + B \left[\Gamma_{i,j} \left(\dfrac{v_{i,j+1} - v_{i,j}}{\Delta y_j} \right) - \Gamma_{i,j-1} \left(\dfrac{v_{i,j} - v_{i,j-1}}{\Delta y_{j-1}}\right)\right] \Delta x_i \\ &= -v_{i,j} B \left(\begin{array}{l} \Gamma\rvert_{x^f_{i+1}} \dfrac{\Delta y_j + \Delta y_{j-1}}{\Delta x_{i+1} + \Delta x_{i}} - \Gamma\rvert_{x^f_{i}} \dfrac{\Delta y_j + \Delta y_{j-1}}{\Delta x_{i} + \Delta x_{i-1}} \\ + \Gamma_{i,j} \frac{\Delta x_i}{\Delta y_j} + \Gamma_{i,j-1} \frac{\Delta x_i}{\Delta y_{j-1}} \end{array} \right) \\ &\ \ \ \ + v_{i-1,j} B \Gamma\rvert_{x^f_{i}} \dfrac{\Delta y_j + \Delta y_{j-1}}{\Delta x_{i} + \Delta x_{i-1}} + v_{i+1,j} B \Gamma\rvert_{x^f_{i+1}} \dfrac{\Delta y_j + \Delta y_{j-1}}{\Delta x_{i+1} + \Delta x_{i}} \\ &\ \ \ \ + v_{i,j-1} B \Gamma_{i,j-1} \dfrac{\Delta x_i}{\Delta y_{j-1}} + v_{i,j+1} B \Gamma_{i,j} \dfrac{\Delta x_i}{\Delta y_j} \\ &= -(D^{vx}_{i,j} + D^{vx}_{i+1,j} + D^{vy}_{i,j} + D^{vy}_{i,j+1}) v_{i,j} \\ &\ \ \ \ + D^{vx}_{i,j} v_{i-1,j} + D^{vx}_{i+1,j} v_{i+1,j} + D^{vy}_{i,j} v_{i,j-1} + D^{vy}_{i,j+1} v_{i,j+1} \end{split} \end{equation*}$

where:

(21) $\begin{equation*} D^{vx}_{i,j} = B \Gamma\rvert_{x^f_{i}} \frac{\Delta y_j + \Delta y_{j-1}}{\Delta x_{i} + \Delta x_{i-1}},\ \ D^{vy}_{i,j} = B \Gamma_{i,j-1} \frac{\Delta x_i}{\Delta y_{j-1}} \end{equation*}$

With $D$ we denote the diffusion conductance at cell face. To calculate $D^{vx}_{i,j}$ we caclute $\Gamma\rvert_{x^f_{i}}$ using bilinear interpolation as we already descibed. To calculate $D^{vy}_{i,j}$ no interpolation is needed as $\Gamma$ is defined at nodes $(x^v_i,y^v_j)$ .

Integration of unsteady term

(22) $\begin{equation*} \begin{split} &\int_{x^f_i}^{x^f_{i+1}}\!\!\!\int_{y^f_{j}}^{y^v_j} \frac{\partial (\rho v)}{\partial t} dy dx \approx \frac{\partial (\rho v)}{\partial t} \frac{\Delta y_j}{2} \Delta x_i = \frac{\rho_{i,j} v_{i,j} - \rho^{n-1}_{i,j} v^{n-1}_{i,j}}{\Delta t} \frac{\Delta y_j}{2} \Delta x_i \end{split} \end{equation*}$

(23) $\begin{equation*} \begin{split} \int_{x^f_i}^{x^f_{i+1}}\!\!\!\int_{y^v_{j-1}}^{y^f_j} \frac{\partial (\rho v)}{\partial t} dy dx \approx \frac{\partial (\rho v)}{\partial t} \frac{\Delta y_{j-1}}{2} \Delta x_i = \frac{\rho_{i,j-1} v_{i,j} - \rho^{n-1}_{i,j-1} v^{n-1}_{i,j}}{\Delta t} \frac{\Delta y_{j-1}}{2} \Delta x_i \end{split} \end{equation*}$

The superscript $n-1$ denotes variables from previous time step, $\Delta t$ is a time step.
The integral overall $v$ control volume is:

(24) $\begin{equation*} \begin{split} &\int_{x^f_i}^{x^f_{i+1}}\!\!\!\int_{y^v_{j-1}}^{y^v_j} \frac{\partial (\rho v)}{\partial t} dy dx \\ &\approx \frac{\rho_{i,j} v_{i,j} - \rho^{n-1}_{i,j} v^{n-1}_{i,j}}{\Delta t} \frac{\Delta y_j}{2} \Delta x_i + \frac{\rho_{i,j-1} v_{i,j} - \rho^{n-1}_{i,j-1} v^{n-1}_{i,j}}{\Delta t} \frac{\Delta y_{j-1}}{2} \Delta x_i \\ &= (\rho_{i,j} \Delta y_{j} + \rho_{i,j-1} \Delta y_{j-1}) \frac{\Delta x_i}{2 \Delta t} v_{i,j} - (\rho^{n-1}_{i,j} \Delta y_{j} + \rho^{n-1}_{i,j-1} \Delta y_{j-1}) \frac{\Delta x_i}{2 \Delta t} v^{n-1}_{i,j} \end{split} \end{equation*}$

Integration of source term

The finite-difference representation of the source term $S^v$ (11) is expressed as:

(25) $\begin{equation*} S^v = S^v_c + S^v_p v_{i,j}, \end{equation*}$

where the coefficient $S^v_p$ is always less than or equal to zero; otherwise we could obtain an instability and divergence of the computational process, or physically unrealistic solutions .

Integration of the terms of the source term:

(26) $\begin{equation*} \begin{split} &\int_{y^v_{j-1}}^{y^v_j}\!\!\int_{x^f_i}^{x^f_{i+1}} \frac{\partial}{\partial x} \left(\Gamma \frac{\partial u}{\partial y} \right) dx dy = \int_{y^v_{j-1}}^{y^v_j} \left[ \left(\Gamma \frac{\partial u}{\partial y} \right)_{x^f_{i+1}} - \left(\Gamma \frac{\partial u}{\partial y} \right)_{x^f_{i}} \right] dy \\ &\ \ \approx \Gamma\rvert_{x^f_{i+1}} (u_{i+1,j} - u_{i+1,j-1}) - \Gamma\rvert_{x^f_{i}} (u_{i,j} - u_{i,j-1}) \\ &\int_{x^f_i}^{x^f_{i+1}}\!\!\int_{y^v_{j-1}}^{y^v_j} \frac{\partial}{\partial y} \left(\Gamma \frac{\partial u}{\partial x} \right) dy dx = \int_{x^f_i}^{x^f_{i+1}} \left[ \left(\Gamma \frac{\partial u}{\partial x} \right)_{y^v_{j}} - \left(\Gamma \frac{\partial u}{\partial x} \right)_{y^v_{j-1}} \right] dx \\ &\ \ \approx \Gamma_{i,j} (u_{i+1,j} - u_{i,j}) - \Gamma_{i,j-1} (u_{i+1,j-1} - u_{i,j-1}) \\ &\int_{y^v_{j-1}}^{y^v_j}\!\!\int_{x^f_i}^{x^f_{i+1}} (\rho g_y) dx dy = g_y \int_{y^v_{j-1}}^{y^v_j}\!\!\int_{x^f_i}^{x^f_{i+1}} \rho \ dx dy \\ &\ \ \approx g_y \left( \rho_{i,j} \frac{\Delta y_j}{2} - \rho_{i,j-1} \frac{\Delta y_{j-1}}{2} \right) \Delta x_i \end{split} \end{equation*}$

The integration of term $\frac{\partial}{\partial y}\left(\Gamma \frac{\partial v}{\partial y}\right)$ is presented in (19).
After some additional rearrangements the source term takes the form:

(27) $\begin{equation*} \begin{split} S^v_p = - \frac{B}{3} \left(\Gamma_{i,j-1} \frac{\Delta x_i}{\Delta y_{j-1}} + \Gamma_{i,j} \frac{\Delta x_i}{\Delta y_j}\right) = -\frac{1}{3}\left(D^{vy}_{i,j} + D^{vy}_{i,j+1}\right) \end{split} \end{equation*}$

(28) $\begin{equation*} \begin{split} S^v_c =& B \left[\frac{1}{3} \Gamma_{i,j} \frac{\Delta x_i}{\Delta y_j} v_{i,j+1} - \Gamma\!\!\mid_{x^f_i}\left(u_{i,j} - u_{i,j-1}\right) + \frac{1}{3} \Gamma_{i,j-1} \frac{\Delta x_i}{\Delta y_{j-1}} v_{i,j-1} \right. \\ & \ \ \ + \frac{2}{3} \Gamma_{i,j-1}\left(u_{i+1,j-1} - u_{i,j-1}\right) + \Gamma\!\!\mid_{x^f_{i+1}}\left(u_{i+1,j} - u_{i+1,j-1}\right) \\ & \left. \ \ \ - \frac{2}{3} \Gamma_{i,j}\left(u_{i+1,j} - u_{i,j}\right)\right]\\ &+ g_y\left(\rho_{i,j}\frac{\Delta y_j}{2} + \rho_{i,j-1}\frac{\Delta y_{j-1}}{2}\right)\Delta x_i \end{split} \end{equation*}$

You can see that $S^v_p$ is always less or equal to zero ( $S^v_p = 0,\ \text{if}\ B = 0$ ). It is worth noting that when considering an extended hydrodynamic system of equations high order terms such as thermal stresses etc. must enter in the source therm. In this case, it is important to transform the source term in a such way so that the condition (40) is fulfilled.

The numerical equation for vertical velocity

After rearrangement the numerical equation for $v$ reads:

(29) $\begin{equation*} v_{i,j} = \frac{a^v_1 v_{i-1,j} + a^v_2 v_{i+1,j} + a^v_3 v_{i,j-1} + a^v_4 v_{i,j+1} + b^v}{a^v_0} - d^v_{i,j} (p_{i,j} - p_{i,j-1}), \end{equation*}$

where the coefficients are:

(30) $\begin{equation*}\begin{split} a^v_0 =& a^v_1 + a^v_2 + a^v_3 + a^v_4 + 0.5 \left(F^x_{i+1,j} - F^x_{i,j} + F^x_{i+1,j-1} - F^x_{i,j-1}\right) \\ & + \overline{F}^y_{i,j} - \overline{F}^y_{i,j-1} - S^v_p + \left(\rho_{i,j} \Delta y_j + \rho_{i,j-1} \Delta y_{j-1} \right) \frac{\Delta x_i}{2 \Delta t} \\ a^v_1 =& a^{vc}_1 + D^{vx}_{i,j} \\ a^v_2 =& a^{vc}_2 + D^{vx}_{i+1,j} \\ a^v_3 =& a^{vc}_3 + D^{vy}_{i,j} \\ a^v_4 =& a^{vc}_4 + D^{vy}_{i,j+1} \\ b^v =& \left(\rho^{n-1}_{i,j} \Delta y_j + \rho^{n-1}_{i,j-1} \Delta y_{j-1}\right)\frac{\Delta x_i}{2 \Delta t}v^{n-1}_{i,j} + S^v_c,\ d^v_{i,j} = \frac{A}{a^v_0}\Delta x_i, \end{split}\end{equation*}$

(31) $\begin{equation*}\begin{split} a^{vc}_0 =& a^{vc}_1 + a^{vc}_2 + a^{vc}_3 + a^{vc}_4 \\ & + 0.5(F^x_{i+1,j} - F^x_{i,j} + F^x_{i+1,j-1} - F^x_{i,j-1}) + \overline{F}^y_{i,j} - \overline{F}^y_{i,j-1} \\ a^{vc}_1 =& 0.5\left[max\left(0, F^x_{i,j}\right)+max\left(0,F^x_{i,j-1}\right)\right] \\ a^{vc}_2 =& 0.5\left[max\left(0, -F^x_{i+1,j}\right) + max\left(0, -F^x_{i+1,j-1}\right)\right] \\ a^{vc}_3 =& max\left(0, \overline{F}^y_{i,j-1}\right)\\ a^{vc}_4 =& max\left(0, -\overline{F}^y_{i,j}\right). \end{split}\end{equation*}$

Derivation and rearrangement of terms of equations are complicated and could be a source of errors. To reduce mistakes, and make derivation faster, and easier we advise you to use appropriate software to do symbolic calculations. You can download the derivation of the numerical equations of the algorithm SIMPLE-TS using a first-order upwind scheme for the approximation of convective terms with Mathematica from here:
SIMPLE_TS_2D_staggered_mesh_First_order_upwind_scheme_approximation_of_convective_terms.nb
SIMPLE_TS_2D_staggered_mesh_First_order_upwind_scheme_approximation_of_convective_terms.pdf

A second order total variation diminishing (TVD) scheme can be used to approximate convective terms. A good explanation of implementation of TVD schemes in finite volume methods is textook of Versteeg and Malalasekera, Chapter 5 The finite volume method for convection-diffusion problems, . You can download derivation of the numerical equations of the algorithm SIMPLE-TS using second order TVD scheme for approximation of convective terms with Mathematica from here:
SIMPLE_TS_2D_staggered_mesh_Second_order_TVD_schemes_approximation_of_convective_terms.nb
SIMPLE_TS_2D_staggered_mesh_Second_order_TVD_schemes_approximation_of_convective_terms.pdf

Finite volume discretization of the pressure equation

After integration over the control volume (from $x^f_i$ to $x^f_{i+1}$ along x-axis and from $y^f_j$ to $y^f_{j+1}$ along y-axis (Fig. 1)), the continuity equation (1) takes the form:

(32) $\begin{equation*} \frac{\partial \rho}{\partial t}\Delta x_i \Delta y_j + \left(\rho^u_{i+1,j}u_{i+1,j} - \rho^u_{i,j}u_{i,j}\right)\Delta y_j + \left(\rho^v_{i,j+1}v_{i,j+1} - \rho^v_{i,j}v_{i,j}\right)\Delta x_i = 0 \end{equation*}$

Here pseudo velocities are used in the same way to SIMPLER. Therefore, the numerical equations for $u$ and $v$ (29) can be written in form:

(33) $\begin{equation*} u_{i,j} = \hat{u}_{i,j} - d^u_{i,j}\left(p_{i,j} - p_{i-1,j}\right) \end{equation*}$

(34) $\begin{equation*} v_{i,j} = \hat{v}_{i,j} - d^v_{i,j}\left(p_{i,j} - p_{i,j-1}\right) \end{equation*}$

where $\hat{u}_{i,j}$ and $\hat{v}_{i,j}$ are pseudo velocities:

(35) $\begin{equation*} \hat{u}_{i,j} = \frac{a^u_1 u_{i-1,j} + a^u_2 u_{i+1,j} + a^u_3 u_{i,j-1} + a^u_4 u_{i,j+1} + b^u}{a^u_0} \end{equation*}$

(36) $\begin{equation*} \hat{v}_{i,j} = \frac{a^v_1 v_{i-1,j} + a^v_2 v_{i+1,j} + a^v_3 v_{i,j-1} + a^v_4 v_{i,j+1} + b^v}{a^v_0} \end{equation*}$

An important part of our scheme is obtained by the substitution of density with pressure in the unsteady term in equation (32) (a similar idea was used in PISO ):

(37) $\begin{equation*} \frac{\partial \rho}{\partial t}\Delta x_i \Delta y_j = \frac{\partial \left(\frac{p}{T} \right)}{\partial t}\Delta x_i \Delta y_j = \left(\frac{p_{i,j}}{T_{i,j}} - \frac{p^{n-1}_{i,j}}{T^{n-1}_{i,j}}\right)\frac{\Delta x_i \Delta y_j}{\Delta t} \end{equation*}$

We multiply by $\Delta t$ to became numerical equation of pressure more stable for small sime steps and after some rearrangements we obtained the numerical equation for pressure as follows:

(38) $\begin{equation*} a^p_0 p_{i,j} = \left(a^{px}_{i,j}p_{i-1,j} + a^{px}_{i+1,j}p_{i+1,j} + a^{py}_{i,j}p_{i,j-1} + a^{py}_{i,j+1}p_{i,j+1}\right)\Delta t + b^p, \end{equation*}$

where:

(39) $\begin{equation*}\begin{split} &a^p_0 = \frac{1}{T_{i,j}}\Delta x_i \Delta y_j + \left(a^{px}_{i,j} + a^{px}_{i+1,j} + a^{py}_{i,j} + a^{py}_{i,j+1}\right)\Delta t \\ &b^p = \frac{p^{n-1}_{i,j}}{T^{n-1}_{i,j}}\Delta x_i \Delta y_j - \left(b^{px}_{i+1,j} - b^{px}_{i,j} + b^{py}_{i,j+1} - b^{py}_{i,j}\right)\Delta t \\ &a^{px}_{i,j} = \rho^u_{i,j}d^u_{i,j}\Delta y_j,\ a^{py}_{i,j} = \rho^v_{i,j}d^v_{i,j}\Delta x_i \\ &b^{px}_{i,j} = \rho^u_{i,j}\hat{u}_{i,j}\Delta y_j,\ b^{py}_{i,j} = \rho^v_{i,j}\hat{v}_{i,j}\Delta x_i \end{split}\end{equation*}$

Scarborough has shown that a sufficient condition for the convergence of an iterative method can be expressed in terms of the values of the coefficients of the numerical equations as follows:

(40) $\begin{equation*}\begin{split} \frac{\left|a^{px}_{i,j} \Delta t\right| + \left|a^{px}_{i+1,j} \Delta t \right| + \left|a^{py}_{i,j} \Delta t \right| + \left|a^{py}_{i,j+1} \Delta t \right|}{\left|a^p_0\right|} \begin{cases} \leq 1& \text{at all nodes},\\ < 1& \text{at one node at least} \end{cases} \end{split}\end{equation*}$

For SIMPLE-TS the validity of condition (40) can be proved easily. The coefficients (39) are positive, therefore, after the substitution of (39) in (40) the condition (40) takes the form:

(41) $\begin{equation*} \frac{\left(a^{px}_{i,j} + a^{px}_{i+1,j} + a^{py}_{i,j} + a^{py}_{i,j+1}\right)\Delta t}{\frac{1}{T_{i,j}}\Delta x_i \Delta y_j + \left(a^{px}_{i,j} + a^{px}_{i+1,j} + a^{py}_{i,j} + a^{py}_{i,j+1}\right)\Delta t} < 1. \end{equation*}$

A similar substitution done for methods containing term $\rho - \rho^{n-1}$ (like SIMPLE) in (40) gives the expression:

(42) $\begin{equation*} \frac{\left(a^{px}_{i,j} + a^{px}_{i+1,j} + a^{py}_{i,j} + a^{py}_{i,j+1}\right)\Delta t}{\left(a^{px}_{i,j} + a^{px}_{i+1,j} + a^{py}_{i,j} + a^{py}_{i,j+1}\right)\Delta t} = 1 \end{equation*}$

and to make the iterative process converging one has to use under-relaxation coefficients appropriately included in the iteration procedure.

The importance and difference of numerical pressure equation used in SIMPLE-TS and other SIMPLE-like methods are discussed in detail using a simplified test case of One-dimensional unsteady, isothermal pressure-driven flow in a duct which is presented in one of the following sections.

Finite volume discretization of the energy equation

In the SIMPLE-TS algorithm the temperature is calculated from the energy equation (4). The term $Dp/Dt$ in the right-hand side in the equation is inconvenient to apply finite volume method as it is because in contains time derivative of pressure and we have to do interpolations in middle point of pressure and velocities where they are not defined. This can be avoided by transforming the term $Dp/Dt$ using continuity equation (1) and equation of state (5) as follows:

(43) $\begin{equation*} \begin{split} &\rho \frac{DT}{Dt} - C^{T3} \frac{Dp}{Dt} = \rho \frac{DT}{Dt} - C^{T3} \left(\frac{\partial (\rho T)}{\partial t} + u \frac{\partial (\rho T)}{\partial x} + v \frac{\partial (\rho T)}{\partial y} \right) \\ &= \rho \frac{DT}{Dt} - \rho C^{T3} \frac{DT}{Dt} - T C^{T3} \left(\frac{\partial \rho}{\partial t} + u \frac{\partial \rho}{\partial x} + v \frac{\partial \rho}{\partial y} \right) \\ &= \rho (1 - C^{T3}) \frac{DT}{Dt} - T C^{T3} \left[\begin{array}{l} \frac{\partial \rho}{\partial t} + u \frac{\partial \rho}{\partial x} + v \frac{\partial \rho}{\partial y} + \left(\frac{\partial (\rho u)}{\partial x} + \frac{\partial (\rho v)}{\partial y} \right) \\ - \left(\frac{\partial (\rho u)}{\partial x} + \frac{\partial (\rho v)}{\partial y} \right) \end{array} \right] \\ &= \rho (1 - C^{T3}) \frac{DT}{Dt} - T C^{T3} \left[\begin{array}{l} \underbrace{\left(\frac{\partial \rho}{\partial t} + \frac{\partial (\rho u)}{\partial x} + \frac{\partial (\rho v)}{\partial y} \right)}_{\text{=0, since it is continuity equatio}\text{n.}} + u \frac{\partial \rho}{\partial x} + v \frac{\partial \rho}{\partial y} \\ - \left(u \frac{\partial \rho}{\partial x} + v \frac{\partial \rho}{\partial y} \right) - \rho \left( \frac{\partial u}{\partial x} + \frac{\partial v}{\partial y} \right) \end{array} \right] \\ &= \rho (1 - C^{T3}) \frac{DT}{Dt} + \rho T C^{T3} \left( \frac{\partial u}{\partial x} + \frac{\partial v}{\partial y} \right) \\ &= \rho (1 - C^{T3}) \frac{DT}{Dt} + p C^{T3} \left( \frac{\partial u}{\partial x} + \frac{\partial v}{\partial y} \right) \end{split} \end{equation*}$

The energy equation is transformed into convenient to apply finite volume method:

(44) $\begin{equation*}\begin{split} &\underbrace{\frac{\partial (\rho T)}{\partial t}}_{Unsteady\ term} + \underbrace{\frac{\partial (\rho u T)}{\partial x} + \frac{\partial (\rho v T)}{\partial y}}_{Convective\ terms} \\ &= \underbrace{C^{T1}_f\left[\frac{\partial}{\partial x}\left(\Gamma^\lambda \frac{\partial T}{\partial x} \right) + \frac{\partial}{\partial y}\left(\Gamma^\lambda \frac{\partial T}{\partial y} \right) \right]}_{Diffusion\ terms} + \underbrace{C^{T2}_f \Gamma \Phi + C^{T3}_f p \left(\frac{\partial u}{\partial x} + \frac{\partial v}{\partial y}\right)}_{Source\ term\left(S^T\right)}, \end{split}\end{equation*}$

where:

(45) $\begin{equation*} C^{T1}_f = \frac{C^{T1}}{1 - C^{T3}},\ C^{T2}_f = \frac{C^{T2}}{1 - C^{T3}},\ C^{T3}_f = \frac{C^{T3}}{C^{T3} - 1}. \end{equation*}$

The integration of energy equation is over a control volume for field variables: from $x^f_i$ to $x^f_{i+1}$ along x-axis and from $y^f_j$ to $y^f_{j+1}$ along y-axis.

Integration of convective terms

We used first-order upwind scheme for the integration of convective terms (44). It worth nothing to implement second order total variation diminishing (TVD) schemes.

(46) $\begin{equation*} \begin{split} &\int_{y^f_{j}}^{y^f_{j+1}}\!\!\int_{x^f_i}^{x^f_{i+1}} \frac{\partial (\rho u T)}{\partial x} dx dy = \int_{y^f_{j}}^{y^f_{j+1}} \left[ (\rho u T )_{x^f_{i+1}} - (\rho u T )_{x^f_{i}} \right] dy \\ &\approx \left[ max\left(0,F^x_{i+1,j} \right) T_{i,j} - max\left(0,-F^x_{i+1,j} \right) T_{i+1,j}\right] \\ &\ \ \ - \left[ max\left(0,F^x_{i,j} \right) T_{i-1,j} - max\left(0,-F^x_{i,j} \right) T_{i,j}\right] \end{split} \end{equation*}$

(47) $\begin{equation*} \begin{split} &\int_{x^f_i}^{x^f_{i+1}}\!\!\int_{y^f_{j}}^{y^f_{j+1}} \frac{\partial (\rho v T)}{\partial y} dy dx = \int_{x^f_i}^{x^f_{i+1}} \left[ (\rho v T )_{y^f_{j+1}} - (\rho v T )_{y^f_{j}} \right] dx \\ &\approx \left[ max\left(0,F^y_{i,j+1} \right) T_{i,j} - max\left(0,-F^y_{i,j+1} \right) T_{i,j+1}\right] \\ &\ \ \ - \left[ max\left(0,F^y_{i,j} \right) T_{i,j-1} - max\left(0,-F^y_{i,j} \right) T_{i,j}\right] \end{split} \end{equation*}$

From (46) and (47) follows:

(48) $\begin{equation*} \begin{split} &\int_{y^f_{j}}^{y^f_{j+1}}\!\!\int_{x^f_i}^{x^f_{i+1}} \left[ \frac{\partial (\rho u T)}{\partial x} + \frac{\partial (\rho v T)}{\partial y} \right] dx dy \\ &\approx \left[ max\left(0,F^x_{i+1,j} \right) T_{i,j} - max\left(0,-F^x_{i+1,j} \right) T_{i+1,j}\right] \\ &\ \ \ - \left[ max\left(0,F^x_{i,j} \right) T_{i-1,j} - max\left(0,-F^x_{i,j} \right) T_{i,j}\right] \\ &\ \ \ \left[ max\left(0,F^y_{i,j+1} \right) T_{i,j} - max\left(0,-F^y_{i,j+1} \right) T_{i,j+1}\right] \\ &\ \ \ - \left[ max\left(0,F^y_{i,j} \right) T_{i,j-1} - max\left(0,-F^y_{i,j} \right) T_{i,j}\right] \\ &= \left[ max\left(0,-F^x_{i,j} \right) + max\left(0,F^x_{i+1,j} \right) + max\left(0,-F^y_{i,j} \right) + max\left(0,F^y_{i,j+1} \right) \right] T_{i,j} \\ &\ \ \ - \left[ \begin{array}{l} max\left(0,F^x_{i,j} \right) T_{i-1,j} + max\left(0,-F^x_{i+1,j} \right) T_{i+1,j} \\ + max\left(0,F^y_{i,j} \right) T_{i,j-1} + max\left(0,-F^y_{i,j+1} \right) T_{i,j+1} \end{array} \right] \\ &= a^{Tc}_0 T_{i,j} - (a^{Tc}_1 T_{i-1,j} + a^{Tc}_2 T_{i+1,j} + a^{Tc}_3 T_{i,j-1} + a^{Tc}_4 T_{i,j+1}) \end{split} \end{equation*}$

where the superscipt $Tc$ denotes temperature ( $T$ ) and convection ( $c$ ), subscripts 0, 1, 2, 3, and 4 denote coefficients for $T_{i,j}$ , $T_{i-1,j}$ , $T_{i+1,j}$ , $T_{i,j-1}$ , and $T_{i,j+1}$ , respectively. Derived coefficients are:

(49) $\begin{equation*} \begin{split} a^{Tc}_0 &= max\left(0,-F^x_{i,j} \right) + max\left(0,F^x_{i+1,j} \right) + max\left(0,-F^y_{i,j} \right) + max\left(0,F^y_{i,j+1} \right) \\ &= a^{Tc}_1 + a^{Tc}_2 + a^{Tc}_3 + a^{Tc}_4 + F^x_{i+1,j} - F^x_{i,j} + F^y_{i,j+1} - F^y_{i,j} \\ a^{Tc}_1 &= max\left(0,F^x_{i,j} \right),\ \ a^{Tc}_2 = max\left(0,-F^x_{i+1,j} \right) \\ a^{Tc}_3 &= max\left(0,F^y_{i,j} \right),\ a^{Tc}_4 = max\left(0,-F^y_{i,j+1} \right) \end{split} \end{equation*}$

Integration of diffusion terms

Integration of diffusion terms in equation (44) over temperature control volume:

(50) $\begin{equation*} \begin{split} &\int_{y^f_{j}}^{y^f_{j+1}}\!\!\int_{x^f_i}^{x^f_{i+1}} C^{T1}_f \frac{\partial}{\partial x} \left( \Gamma^\lambda \frac{\partial T}{\partial x} \right) dx dy \\ &= C^{T1}_f \int_{y^f_{j}}^{y^f_{j+1}} \left[ \left( \Gamma^\lambda \frac{\partial T}{\partial x} \right)_{x^f_{i+1}} - \left( \Gamma^\lambda \frac{\partial T}{\partial x} \right)_{x^f_{i}} \right] dy \\ &\approx C^{T1}_f \left( \Gamma^\lambda\rvert_{x^f_{i+1}} \frac{T_{i+1,j} - T_{i,j}}{\frac{1}{2} (\Delta x_{i+1} + \Delta x_{i})} - \Gamma^\lambda\rvert_{x^f_{i}} \frac{T_{i,j} - T_{i-1,j}}{\frac{1}{2} (\Delta x_{i} + \Delta x_{i-1})} \right) \Delta y_j \end{split} \end{equation*}$

(51) $\begin{equation*} \begin{split} &\int_{x^f_i}^{x^f_{i+1}}\!\!\int_{y^f_{j}}^{y^f_{j+1}} C^{T1}_f \frac{\partial}{\partial y} \left( \Gamma^\lambda \frac{\partial T}{\partial y} \right) dy dx \\ &= C^{T1}_f \int_{x^f_i}^{x^f_{i+1}} \left[ \left( \Gamma^\lambda \frac{\partial T}{\partial y} \right)_{y^f_{j+1}} - \left( \Gamma^\lambda \frac{\partial T}{\partial y} \right)_{y^f_{j}} \right] dx \\ &\approx C^{T1}_f \left( \Gamma^\lambda\rvert_{y^f_{j+1}} \frac{T_{i,j+1} - T_{i,j}}{\frac{1}{2} (\Delta y_{j+1} + \Delta y_{j})} - \Gamma^\lambda\rvert_{y^f_{j}} \frac{T_{i,j} - T_{i,j-1}}{\frac{1}{2} (\Delta y_{j} + \Delta y_{j-1})} \right) \Delta x_i \end{split} \end{equation*}$

where $\Gamma^\lambda$ is defined at nodes of scalar variable. A harmonic average between two neighboring nodes is used to calculate $\Gamma^\lambda\rvert_{x^f_{i}}$ and $\Gamma^\lambda\rvert_{y^f_{j}}$ as it is used by Patabkar :

(52) $\begin{equation*} \Gamma ^\lambda \!\!\mid_{x^f_i} = \frac{(\Delta x_{i-1} + \Delta x_i)\Gamma^\lambda_{i-1,j} \Gamma^\lambda_{i,j}}{\Delta x_{i-1} \Gamma^\lambda_{i,j} + \Delta x_i \Gamma^\lambda_{i-1,j}},\ \Gamma ^\lambda \!\!\mid_{y^f_j} = \frac{(\Delta y_{j-1} + \Delta y_j)\Gamma^\lambda_{i,j-1} \Gamma^\lambda_{i,j}}{\Delta y_{j-1} \Gamma^\lambda_{i,j} + \Delta y_j \Gamma^\lambda_{i,j-1}} \end{equation*}$

The integral of diffusion terms over-temperature control volume is:

(53) $\begin{equation*} \begin{split} &\int_{y^f_{j}}^{y^f_{j+1}}\!\!\int_{x^f_i}^{x^f_{i+1}} C^{T1}_f \left[ \frac{\partial}{\partial x} \left( \Gamma^\lambda \frac{\partial T}{\partial x} \right) + \frac{\partial}{\partial y} \left( \Gamma^\lambda \frac{\partial T}{\partial y} \right) \right] dx dy \\ &\approx C^{T1}_f \left( \Gamma^\lambda\rvert_{x^f_{i+1}} \frac{T_{i+1,j} - T_{i,j}}{\frac{1}{2} (\Delta x_{i+1} + \Delta x_{i})} - \Gamma^\lambda\rvert_{x^f_{i}} \frac{T_{i,j} - T_{i-1,j}}{\frac{1}{2} (\Delta x_{i} + \Delta x_{i-1})} \right) \Delta y_j \\ &\ \ \ + C^{T1}_f \left( \Gamma^\lambda\rvert_{y^f_{j+1}} \frac{T_{i,j+1} - T_{i,j}}{\frac{1}{2} (\Delta y_{j+1} + \Delta y_{j})} - \Gamma^\lambda\rvert_{y^f_{j}} \frac{T_{i,j} - T_{i,j-1}}{\frac{1}{2} (\Delta y_{j} + \Delta y_{j-1})} \right) \Delta x_i \\ &= - (D^{Tx}_{i,j} + D^{Tx}_{i+1,j} + D^{Ty}_{i,j} + D^{Ty}_{i,j+1}) T_{i,j} \\ &\ \ \ + D^{Tx}_{i,j} T_{i-1,j} + D^{Tx}_{i+1,j} T_{i+1,j} + D^{Ty}_{i,j} T_{i,j-1} + D^{Ty}_{i,j+1} T_{i,j+1} \end{split} \end{equation*}$

where the diffusion coefficients are:

(54) $\begin{equation*} D^{Tx}_{i,j} = C^{T1}_f \Gamma ^\lambda \!\!\mid_{x^f_i} \frac{\Delta y_j}{0.5 (\Delta x_i + \Delta x_{i-1})},\ D^{Ty}_{i,j} = C^{T1}_f \Gamma ^\lambda \!\!\mid_{y^f_j} \frac{\Delta x_i}{0.5 (\Delta y_j + \Delta y_{j-1})} \end{equation*}$

Integration of unsteady term

(55) $\begin{equation*} \begin{split} \int_{y^f_{j}}^{y^f_{j+1}}\!\!\int_{x^f_i}^{x^f_{i+1}} \frac{\partial (\rho T)}{\partial t} dxdy = \frac{\partial (\rho T)}{\partial t} \Delta x_i \Delta y_j = \frac{\rho_{i,j} T_{i,j} - \rho^{n-1}_{i,j} T^{n-1}_{i,j}}{\Delta t} \Delta x_i \Delta y_j \end{split} \end{equation*}$

Integration of source term

The finite-difference representation of the source term $S^T$ is expressed in same way as $S^v$ as follows:

(56) $\begin{equation*} S^T = S^T_c + S^T_p T_{i,j}, \end{equation*}$

where:

(57) $\begin{equation*} \begin{split} S^T_p =& 0 \\ S^T_{c_{i,j}} =& C^{T2}_f \Gamma _{i,j} \left\{ 2 \left[\left(\frac{u_{i+1,j} - u_{i,j}}{\Delta x_i}\right)^2 + \left(\frac{v_{i,j+1} - v_{i,j}}{\Delta y_j}\right)^2 \right] \right. \\ & \ \ \ \ \ \ \ \ \ \ \ \ \ + \left(\frac{v(x^f_{i+1},y^v_j) - v(x^f_i,y^v_j)}{\Delta x_i} + \frac{u(x^v_i,y^f_{j+1}) - u(x^v_i,y^f_j)}{\Delta y_j}\right)^2 \\ & \left. \ \ \ \ \ \ \ \ \ \ \ \ \ -\frac{2}{3}\left(\frac{u_{i+1,j} - u_{i,j}}{\Delta x_i} + \frac{v_{i,j+1} - v_{i,j}}{\Delta y_j}\right)^2 \right\} \Delta x_i \Delta y_j \\ &+ C^{T3}_f p_{i,j} \left(\frac{u_{i+1,j} - u_{i,j}}{\Delta x_i} + \frac{v_{i,j+1} - v_{i,j}}{\Delta y_j}\right) \Delta x_i \Delta y_j \end{split} \end{equation*}$

To interpolate velocities $u(x^v_i,y^f_{j+1}), u(x^v_i,y^f_j), v(x^f_{i+1},y^v_j)$ and $v(x^f_i,y^v_j)$ a bi-linear interpolation between four neighboring nodes is used for each one.

The numerical equation for temperature

The discrete equation for temperature is:

(58) $\begin{equation*}\begin{split} a^T_0 T_{i,j} =& \Delta t \left(a^T_1 T_{i-1,j} + a^T_2 T_{i+1,j} + a^T_3 T_{i,j-1} + a^T_4 T_{i,j+1} + S^T_{c_{i,j}}\right) \\ & + \rho^{n-1}_{i,j} T^{n-1}_{i,j} \Delta x_i \Delta y_j, \end{split}\end{equation*}$

where:

(59) $\begin{equation*}\begin{split} a^T_0 =& \Delta t \left(a^T_1 + a^T_2 + a^T_3 + a^T_4 + F^x_{i+1,j} - F^x_{i,j} + F^y_{i,j+1} - F^y_{i,j}\right) \\ & + \rho_{i,j} \Delta x_i \Delta y_j \\ a^T_1 =& max\left(0, F^x_{i,j}\right) + D^{Tx}_{i,j},\ \ \ \ \ \ \ a^T_3 = max\left(0, F^y_{i,j}\right) + D^{Ty}_{i,j} \\ a^T_2 =& max\left(0, -F^x_{i+1,j}\right) + D^{Tx}_{i+1,j}, a^T_4 = max\left(0, -F^y_{i,j+1}\right) + D^{Ty}_{i,j+1}. \\\end{split}\end{equation*}$

Here we multiply by time step ( $\Delta t$ ) to increase numerical stability when use small time step in calculations, as we did in numerical equation for $v$ .

Algorithm SIMPLE-TS

Using the results of the previous section we are ready to construct the general algorithm of the new method. The SIMPLE-TS algorithm for solving unsteady, compressible, viscous and heat-conducting gas flows consists of the following steps:
The SIMPLE-TS algorithm
Step 1. Initialize $u$ , $v$ , $p$ , $T$ and compute $\rho$ using equation (5). Set time step $\Delta t$ .
Start loop 1:
    $\bullet$ Step 2. Set the initial condition for the calculated time step: $t := t + \Delta t$ , $u^{n-1} = u$ , $v^{n-1} = v$ , $p^{n-1} = p$ , $T^{n-1} = T$ , $\rho^{n-1} = \rho$ .
    $\bullet$ Start loop 2 (calculating a state for a new time step):
       $\bullet\bullet$ Step 3. Calculate convective and diffusion fluxes: $F^x$ (7), $F^y$ (8), $D^{ux}$ , $D^{uy}$ , $D^{vx}$ (21), $D^{vy}$ (21), $D^{Tx}$ (54), $D^{Ty}$ (54).
       $\bullet\bullet$ Step 4. Calculate pseudo velocities $\hat u$ (35), $\hat v$ (36), coefficients $d^u$ and $d^v$ (30), and the coefficients for pressure equation $a^{px}$ , $a^{py}$ , $b^{px}$ , $b^{py}$ (39).
       $\bullet\bullet$ Start loop 3:
          $\bullet\bullet\bullet$ Step 5. Solve the coupled equations for energy (58) and pressure (38).
       $\bullet\bullet$ Stop loop 3. In most cases two iterations are sufficient. In this is discussed in details.
       $\bullet\bullet$ Step 6. Compute $\rho$ (5), using $p$ and $T$ calculated within step 5, as well as velocities: $u$ (33) and $v$ (34).
       $\bullet\bullet$ Convergence of loop 2: Check for convergence of the iteration process for the current time step by using iteration criteria (60). If the iteration criteria are not satisfied go to step 3, otherwise continue.
    $\bullet$ Convergence of loop 1: If the final time is not reached go to step 2, otherwise stop the calculation.
The following iteration criteria must be satisfied:

(60) $\begin{equation*}\begin{split} &\left|p - p^{old}\right|_{max} < \epsilon^p,\ \left|T - T^{old}\right|_{max} < \epsilon^T \\ &\left|u - u^{old}\right|_{max} < \epsilon^u,\ \left|v - v^{old}\right|_{max} < \epsilon^v\end{split}\end{equation*}$

where $p$ , $T$ , $u$ and $v$ are values obtained from the current iteration, $p^{old}$ , $T^{old}$ , $u^{old}$ and $v^{old}$ are values obtained from the previous iteration, $\epsilon^p$ , $\epsilon^T$ , $\epsilon^u$ and $\epsilon^v$ are iteration criteria for $p$ , $T$ , $u$ and $v$ , respectively. When under-relaxation coefficients are included in the iteration scheme the convergence criteria take the following form:

(61) $\begin{equation*}\begin{split} &\frac{\left|p-p^{old}\right|_{max}}{\alpha_p}<\epsilon^p, \frac{\left|T-T^{old}\right|_{max}}{\alpha_T}<\epsilon^T,\\ &\frac{\left|u-u^{old}\right|_{max}}{\alpha_u}<\epsilon^u, \frac{\left|v-v^{old}\right|_{max}}{\alpha_v}<\epsilon^v\end{split}\end{equation*}$

where $\alpha_p$ , $\alpha_T$ , $\alpha_u$ and $\alpha_v$ are under-relaxation coefficients for $p$ , $T$ , $u$ and $v$ , respectively.
A basic idea in the algorithm SIMPLE-TS is to keep the equation of state maximally satisfied within the entire iterative process. This idea is accomplished by the calculation of pressure and temperature in loop 3 and the substitution of density with pressure into the unsteady term using equation of state. The loop 3 reduces the computational time for reaching a given accuracy (see Table. ??). A second benefit from loop 3 is that the mass in the entire computational volume is conserved within a given limits before completing the loop 2, since the loop 3 is run until satisfying criteria $\epsilon^p$ and $\epsilon^T$ .
The velocities, pressure, and temperature fields are calculated at the end of loop 2 using the values of the pseudo velocities obtained in step 3, and coefficients (step 4), which depend directly on the initial fields, are set at the beginning of loop 2.
The interpolation, applied to the middle points as explained in the previous section, is important for a successful implementation of the algorithm SIMPLE-TS.

Convert algorithm SIMPLE to SIMPLE-TS

Converting SIMPLE-like algorithm to SIMPLE-TS could be simple task. The most important is to is to keep the equation of state maximally satisfied within the entire iterative process that it is the main idea implemented in SIMPLE-TS. Here is given an example of converting algorithm SIMPLE to SIMPLE-TS. The following changes have to be accomplished:
1. The most important step is the substitution of density into the unsteady term of the pressure equation using the equation of state. In this way, the coefficients for the pressure equations, used in SIMPLE-TS, satisfy a sufficient criterion for the convergence of the iterative process (41) while the SIMPLE coefficients do not and SIMPLE needs under-relaxation coefficients to ensure convergence.
2. The calculation of energy and pressure equations in one loop (loop 3 of SIMPLE-TS), after calculation of pseudo velocities and before calculation of velocities.
3. The calculation of pseudo velocities.
4. The calculation of the absolute pressure instead of the pressure correction.
Of course the interpolation, applied to the middle points as explained in the previous section, is important for a successful implementation of the new algorithm and should be followed.

One-dimensional unsteady, isothermal pressure driven flow in a duct

To stress on the importance of the substitution of density into the unsteady term of the pressure equation using the equation of state let us consider a simplified calculation of one-dimensional unsteady isothermal pressure driven flow in a duct. The example is simple but useful to highlight the meaning of the replace of density with pressure using the equation of state (37). The diffusion coefficient $\Gamma$ depends on temperature and for this considerations is considered to a constant. The system of equations describing the one-dimensional, isothermal flow without external field is:

(62) $\begin{equation*}\frac{\partial \rho }{\partial t} + \frac{\partial (\rho u)}{\partial x} = 0\end{equation*}$

(63) $\begin{equation*}\frac{\partial(\rho u)}{\partial t} + \frac{\partial(\rho u u)}{\partial x} = -A \frac{\partial p}{\partial x} + \frac{4}{3}B\Gamma\frac{\partial^2 u}{\partial x^2}\end{equation*}$

The problem is calculated by using numerical pressure equation (32) with parameters $A = 1$ and $D^{ux} = 1$ . The inlet and outlet pressures are fixed at ( $p_{in} = 3$ ) and ( $p_{out} = 1$ ), at both inlet and outlet the gradient of the velocity is zero ( $\partial u / \partial x = 0$ ). Initially the fluid is at rest and the pressure is equal to the outlet pressure. The step on OX is $\Delta x = 0.1$ , the time step is $\Delta t = 0.1$ and the length of the duct is 10. At zero time the inlet is opened suddenly and the fluid at pressure $p_{in}$ enters the duct volume. Let us analyse what happens within the first iteration of the both SIMPLE and SIMPLE-TS algorithms.

**Figure 2.** One-dimensional staggered grid

Fig. 1 shows one-dimensional staggered mesh of discussed problem.

The numerical system of equations is:

(64) $\begin{equation*} u_{i} = \hat{u}_{i} - d^u_{i}\left(p_{i} - p_{i-1}\right) \end{equation*}$

(65) $\begin{equation*} a^p_{0_i} p_{i} = (a^{px}_{i} p_{i-1} + a^{px}_{i+1} p_{i+1}) \Delta t + b^p_i, \end{equation*}$

where the pseudo velocity ( $\hat{u}_{i}$ ) and the corresponding coefficients are:

(66) $\begin{equation*} \hat{u}_i = \frac{a^u_1 u_{i-1} + a^u_2 u_{i+1} + b^u}{a^u_{0_i}}, \end{equation*}$

(67) $\begin{equation*} \begin{split} a^u_1 &= a^{uc}_1 + D^{ux}_i, \ \ \ \ \ \ \ \ \ a^u_2 = a^{uc}_2 + D^{ux}_{i+1},\\ a^{uc}_1 &= max\left(0, \overline{F}^x_{i-1} \right), \ \ a^{uc}_2 = max\left(0, -\overline{F}^x_{i} \right),\\ b^u &= \frac{(\rho^{n-1}_{i-1} + \rho^{n-1}_{i}) u^{n-1}_{i} \Delta x 0.5}{\Delta t} + S^u_c,\\ a^u_{0_i} &= a^u_1 + a^u_2 + \overline{F}^x_{i} - \overline{F}^x_{i-1} - S^u_p + \frac{(\rho_{i-1} + \rho_{i}) \Delta x 0.5}{\Delta t},\\ S^u_c &= \frac{B}{3} \left(\frac{u_{i+1} \Gamma_i + u_{i-1} \Gamma_{i-1}}{\Delta x} \right),\ \ \ \ S^u_p = -\frac{D^{ux}_i + D^{ux}_{i+1}}{3},\\ \overline{F}^x_{i} &= 0.5 (F^x_{i} + F^x_{i+1}),\ \ F^x_i = u_i \rho^u_i,\ \ D^{ux}_i = \dfrac{B \Gamma_{i-1}}{\Delta x},\ \ d^u_i = \dfrac{A}{a^u_{0_i}} \end{split} \end{equation*}$

The coefficients with overline ( $\overline{a}^p_{0_i}$ and $\overline{b}^p_i$ ) are used in SIMPLE-like methods SIMPLE, SIMPLER, SIMPLEST-ANL and SIMPLEC. Note that the density is not substituted into the unsteady term of the pressure equation using the equation of state:

(68) $\begin{equation*} \begin{split} \overline{a}^p_{0_i} &= \left(a^{px}_{i} + a^{px}_{i+1} \right)\Delta t \\ \overline{b}^p_i &= -(\rho_{i} - \rho^{n-1}_{i}) \Delta x - \left(b^{px}_{i+1} - b^{px}_{i} \right)\Delta t \\ \end{split} \end{equation*}$

The coefficients used in the pressure equation that corresponds to SIMPLE-TS are:

(69) $\begin{equation*} \begin{split} a^p_{0_i} &= \frac{1}{T_{i}}\Delta x + \left(a^{px}_{i} + a^{px}_{i+1} \right)\Delta t,\ \ b^p_i = \frac{p^{n-1}_{i}}{T^{n-1}_{i}}\Delta x - \left(b^{px}_{i+1} - b^{px}_{i} \right)\Delta t \\ a^{px}_{i} &= \rho^u_{i}d^u_{i},\ \ b^{px}_{i} = \rho^u_{i}\hat{u}_{i} \end{split} \end{equation*}$

Note that the density in unsteady term in the pressure equation is substituted with pressure using equation of state. This is an important new step in SIMPLE-TS.

Let us do the first iteration manually to show where the idea of substitution of density using the equation of state came from. All computations are performed following the steps of the algorithm SIMPLE-TS. As we consider isothermal flow the temperature is denoted as $T$ .

Consider in details the first iteration of the SIMPLE-like methods with pressure coefficients (68) (without substitution of density):

Step 1. Initialize the arrays according to the initial and the boundary conditions of one-dimensional, pressure-driven, isothermal flow.
Initially, the velocity field is zero in all nodes:

(70) $\begin{equation*}u^{n-1}_i = 0,\ \text{for}\ i = 1, 2, ... , N_x + 1\end{equation*}$

The pressure field is:

(71) $\begin{equation*}\begin{split}p^{n-1}_1 =& p_{in}\\p^{n-1}_i =& p_{out},\ \text{for}\ i = 2, 3, ... , N_x\end{split}\end{equation*}$

The flow is isothermal, therefore the density is equal to the pressure field (note that here it is considered non-dimensional case):

(72) $\begin{equation*}\begin{split}\rho^{n-1} = p^{n-1}\end{split}\end{equation*}$

Start loop 1
Step 2. The initial guess used of the calculated time step is equal to the previous time step: $u = u^{n-1}$ , $p = p^{n-1}$ , and $\rho = \rho^{n-1}$ .
Start loop 2
Step 3. Calculate convective and diffusion fluxes.
The velocity filed in the computational domain is zero (70), therefore $F^x = 0$ .
As the problem is isothermal then $B = const$ and $\Gamma=const$ , therefore the diffusion fluxes are equal to constant: $D^{ux}_i = D^{ux} = const$ .
The pseudo velocities are equal to zero $\hat{u} = 0$ and coefficient $a^u_{0_i} = \dfrac{8}{3} D^{ux} + (\rho_{i-1} + \rho_{i}) \dfrac{\Delta x}{2\Delta t}$ .
The internal cells (Fig. 2) are calculated for $i = 2, 3, ... , N_x$ using the coefficients:

(73) $\begin{equation*}\begin{split}a^u_{0_2} &= \dfrac{8}{3} D^{ux} + (\rho_{in} + \rho_{out}) \dfrac{\Delta x}{2\Delta t}= \dfrac{8}{3} D^{ux} + (p_{in} + p_{out}) \dfrac{\Delta x}{2\Delta t} \\a^u_{0_3} &= a^u_{0_4} = ... = a^u_{0_{N_x}} = \dfrac{8}{3} D^{ux} + (\rho_{out} + \rho_{out}) \dfrac{\Delta x}{2\Delta t} = \dfrac{8}{3} D^{ux} + p_{out} \dfrac{\Delta x}{\Delta t} \\d^u_2 &= \dfrac{A}{a^u_{0_2}} \\d^u_3 &= d^u_4 = ... = d^u_{N_x} = \dfrac{A}{a^u_{0_3}}\\a^{px}_2 &= \rho^u_2 d^u_2 = \dfrac{1}{2} (\rho_1 + \rho_2) \dfrac{A}{a^u_{0_2}} = \dfrac{1}{2} (p_{in} + p_{out}) \dfrac{A}{a^u_{0_2}} \\a^{px}_3 &= a^{px}_4 = ... = a^{px}_{N_x} = p_{out}\dfrac{A}{a^u_{0_3}} \\b^{px} &= 0 \\\end{split}\end{equation*}$

The pressure equation coefficients (68) are:
$\overline{a}^p_{0_2} = (a^{px}_2 + a^{px}_3) \Delta t = \left(\dfrac{1}{2} (p_{in} + p{out}) \dfrac{A}{a^u_{0_2}} + p_{out}\dfrac{A}{a^u_{0_3}} \right) \Delta t$
$\overline{a}^p_{0_3} = \overline{a}^p_{0_4} = ... = \overline{a}^p_{0_{N_x-1}} = (a^{px}_3 + a^{px}_4) \Delta t = 2 p_{out}\dfrac{A}{a^u_{0_3}} \Delta t$
$\overline{b}^p = 0$ .
As $\rho - \rho^{n-1} = 0$ , the numerical pressure equation is:

(74) $\begin{equation*}\begin{split}p_i = \dfrac{(a^{px}_{i} p_{i-1} + a^{px}_{i+1} p_{i+1})\Delta t}{\overline{a}^p_{0_i}} = \dfrac{a^{px}_{i} p_{i-1} + a^{px}_{i+1} p_{i+1}}{a^{px}_{i} + a^{px}_{i+1}}\end{split}\end{equation*}$

After the substitution of the coefficients (73) in (74) we obtained the numerical equations for the pressure in cells:

(75) $\begin{equation*}\begin{split}p_2 = \dfrac{a^{px}_{2} p_{in} + a^{px}_{3} p_3}{a^{px}_{2} + a^{px}_{3}}, \end{split}\end{equation*}$

where $a^{px}_{2} \neq a^{px}_{3}$

(76) $\begin{equation*}\begin{split}p_i = \dfrac{a^{px}_{i} p_i + a^{px}_{i+1} p_{i+1}}{a^{px}_{i} + a^{px}_{i}} = \dfrac{p_i + p_{i+1}}{2}, \ \ \ \text{for}\ i = 3, 4, ... , N_x - 1\end{split}\end{equation*}$

Equation (75) shows that the pressure in cell 2 depends on the pressure in the neighbouring cells calculated within the considered time step and does not depend on time step, pressure, density, or velocity evaluated in the previous time step. The obtained pressure along the control volumes $i = 3, 4, ... , N_x - 1$ , (76) after the first iteration varies linearly from $p_2$ to $p_{out}$ . The other important test that conforms incorrect behaviour is when time step approaches zero in equation (76). The obtained solution of the pressure in the channel is linear interpolation of inlet and outlet pressure and does not depends on the initial pressure in the channel. As a result when the pressure increases at the duct inlet the disturbance propagates immediately to the end of the duct (Fig. 3). This behavior of fluid corresponds to incompressible fluid while here the model is compressible fluid. In compressible fluid disturbances propagate with a speed of sound that is a finite value and has to be implemented in the algorithm. The important question that arise to the experienced one is “How did SIMPLE-like methods obtain correct solutions of compressible fluid flows?”. This question will be clarified at the end of considerations about first iteration of SIMPLE-TS (next section).

Consider in details the first iteration of the SIMPLE-TS method with pressure coefficients (69) (with substitution of density with pressure):

Step 1. Initialize the arrays according to the initial and the boundary conditions of one-dimensional, pressure-driven, isothermal flow.
Initially, the velocity field is zero in all nodes:

(77) $\begin{equation*}u^{n-1}_i = 0,\ \text{for}\ i = 1, 2, ... , N_x + 1\end{equation*}$

The pressure field is:

(78) $\begin{equation*}\begin{split}p^{n-1}_1 =& p_{in}\\p^{n-1}_i =& p_{out},\ \text{for}\ i = 2, 3, ... , N_x\end{split}\end{equation*}$

The flow is isothermal, therefore the density is equal to the pressure field (note that here it is considered non-dimensional case):

(79) $\begin{equation*}\begin{split}\rho^{n-1} = p^{n-1}\end{split}\end{equation*}$

(80) $\begin{equation*}\begin{split}a^u_{0_2} &= \dfrac{8}{3} D^{ux} + (\rho_{in} + \rho_{out}) \dfrac{\Delta x}{2\Delta t}= \dfrac{8}{3} D^{ux} + (p_{in} + p_{out}) \dfrac{\Delta x}{2\Delta t} \\a^u_{0_3} &= a^u_{0_4} = ... = a^u_{0_{N_x}} = \dfrac{8}{3} D^{ux} + (\rho_{out} + \rho_{out}) \dfrac{\Delta x}{2\Delta t} = \dfrac{8}{3} D^{ux} + p_{out} \dfrac{\Delta x}{\Delta t} \\d^u_2 &= \dfrac{A}{a^u_{0_2}} \\d^u_3 &= d^u_4 = ... = d^u_{N_x} = \dfrac{A}{a^u_{0_3}}\\a^{px}_2 &= \rho^u_2 d^u_2 = \dfrac{1}{2} (\rho_1 + \rho_2) \dfrac{A}{a^u_{0_2}} = \dfrac{1}{2} (p_{in} + p_{out}) \dfrac{A}{a^u_{0_2}} \\a^{px}_3 &= a^{px}_4 = ... = a^{px}_{N_x} = p_{out}\dfrac{A}{a^u_{0_3}} \\b^{px} &= 0\end{split}\end{equation*}$

Up to here, the procedure is the same as SIMPLE-like methods.
The pressure equation coefficients with sunstitution of density with pressure (69) are:
$a^p_{0_2} = \frac{1}{T}\Delta x + (a^{px}_2 + a^{px}_3) \Delta t$
$a^p_{0_3} = {a}^p_{0_4} = ... = {a}^p_{0_{N_x-1}} = \frac{1}{T}\Delta x + 2 a^{px}_3 \Delta t$
$b^p_i = \frac{p^{n-1}_i}{T}\Delta x$ for $i=1,...,Nx$ .
The numerical pressure equation is:

(81) $\begin{equation*}\begin{split}p_i = \dfrac{(a^{px}_{i} p_{i-1} + a^{px}_{i+1} p_{i+1})\Delta t + b^p_i}{a^p_{0_i}} = \dfrac{(a^{px}_{i} p_{i-1} + a^{px}_{i+1} p_{i+1})\Delta t + \frac{p^{n-1}_i}{T}\Delta x}{\frac{1}{T}\Delta x + (a^{px}_{i} + a^{px}_{i+1})\Delta t}\end{split}\end{equation*}$

After the substitution of the coefficients (80) in (81) we obtained the numerical equations for the pressure in cells and taking into account that accordind to initial boundary conditions $p^{n-1}_1=p_{in}$ and $p^{n-1}_i=p_{out}$ for $i=2,3,...,N_x$ :

(82) $\begin{equation*}\begin{split}p_2 = \dfrac{(a^{px}_{2} p_{in} + a^{px}_{3} p_3) \Delta t + \frac{p_{out}}{T}\Delta x}{\frac{1}{T}\Delta x + (a^{px}_{2} + a^{px}_{3}) \Delta t} \end{split}\end{equation*}$

(83) $\begin{equation*}\begin{split}p_i = \dfrac{a^{px}_{3} ( p_{3} + p_4) \Delta t + \frac{p_{out}}{T}\Delta x}{\frac{1}{T}\Delta x + 2 a^{px}_{3} \Delta t}, \ \ \ \text{for}\ i = 3, 4, ... , N_x - 1\end{split}\end{equation*}$

The solution of the pressure in the channel depends on pressure on previous time step and does not propagates immediately to the end of the channel, see Fig. 3. The behavior of the numerical solution corresponds to the compressible fluid model that it should be. The other important test is when time step approaches zero in equation (83). The obtained solutiuon of the pressure in the channel is equal to the inituial pressure in the channel, i.e. $p_i=p_{out}$ for $i=2,3,…,N_x$ . The behavior of the pressure equation used in SIMPLE-TS method corresponds to the behavior of a compressible fluid.

Here we will answer the question “How did SIMPLE-like methods obtain correct solutions of compressible fluid flows?” from the previous section when we discussed the first iteration of the SIMPLE-like methods. Well having in mind the presented considerations we can set up a SIMPLE-like method to behave as close as possible to the SIMPLE-TS method. SIMPLE-like methods have to use under-relaxation coefficients to reduce the spread of the corrections of the pressure equation (disturbances in fluid) and have to use one iteration for calculation of pressure equation. In this way, the SIMPLE-like method could “update” changes in density that satisfy the equation of state indirectly in the iteration process, that it is actually the main idea implemented in SIMPLE-TS. So, in this way, SIMPLE-like methods can be set up to calculate compressible fluid models. The price/disadvantage of under-relaxation coefficients is that they slow down SIMPLE-like methods approximately from 5 to 10 times compared to SIMPLE-TS.

Considerations and numerical results of a first time step of a modelling of one-dimensional unsteady, isothermal pressure driven flow in a duct

Here are placed obtained pressure profiles at first iteration of unsteady, isothermal pressure-driven compressible flow in a duct. It was shown that in the case without substitution of density (in the pressure equation exist term $\rho - \rho^{n-1}$ ) the exact pressure profile is linear demonstrating a typical incompressible behavior of the fluid while in the case of substitution (37) the solution takes into account the compressibility of the fluid. The pressure profiles shown on Fig. 3 are obtained after the first iteration of the loop 2. In the case without substitution of density with pressure the equation of state is not satisfied within the iteration and a propagation of a disturbance from inlet to outlet is generated in each iteration making difficult the convergence of the iterative process. To ensure convergence of the algorithm in the case without substitution of density one must use under-relaxation coefficients and one iteration of the loop 3. The pressure profile of the exact solution of the pressure equation, when density is substituted, has a form of a compression wave (Fig. 3). Already after the first iteration the pressure profile is close to the final solution which is sought for within the given time step. Thus, the algorithm SIMPLE-TS maximally satisfies the equation of state during the iterative process and this leads to a better approximation of the solution after each iteration. In this way the number of iterations, and correspondingly the computational time, is reduced.
The one-dimensional pressure-driven case is important to show the importance of the idea of substitution of density with pressure with the equation of state in the pressure equation.

**Figure 3.** Pressure profiles of the exact solution of the pressure equation calculated in loop 3, obtained with and without substitution of the density into the unsteady term of the pressure equation.

You can download the program that calculates and plots results of the first time step of one-dimensional unsteady, isothermal pressure-driven flow in a duct from here: One_dimensional_flow_in_a_duct.m

References