Conservation of mass (compressible flow)

Let us denote the velocity calculated from the momentum equation by \boldsymbol{v}^*. Thus:

\boldsymbol{v}^*_P + \boldsymbol{H}_P (\boldsymbol{v}^*) = \boldsymbol{B}_P^{(m-1)} - D_P \nabla p_P^{(m-1)}

(649)

and

\displaystyle \boldsymbol{v}^*_f	\displaystyle = \overline{\boldsymbol{v} }^*_f - \overline{D}_f \left [ \frac{p... ... \overline{\nabla p}_f^{(m-1)} \cdot \boldsymbol{j}_F \right ] \boldsymbol{j}_f	(650)
	\displaystyle = \overline{\boldsymbol{v} }^*_f + \overline{D}_f (\overline{\nabla p}_f^{(m-1)} - \nabla p_f^{(m-1)})	(651)

The resulting mass flow \dot{m}^*_f through face f corresponds to:

\dot{m}^*_f = \overrightarrow{\rho }_f^{(m-1)} A_f \boldsymbol{v}^*_f \cdot \boldsymbol{n}_f,

(652)

where \overrightarrow{\rho }_f^{(m-1)} is obtained from

\overrightarrow{\rho } ^{(m-1)} = \frac{\overrightarrow{p} ^{(m-1)} }{r \overline{T} ^{(m-1)} }

(653)

while at the element centers we have (e.g. at P):

\rho_P^{(m-1)}=\frac{p_P ^{(m-1)}}{r T_P ^{(m-1)} }.

(654)

Now, mass conservation requires:

\frac{\partial \rho }{\partial t} + (\rho v_j),_j =0,

(655)

which after solving the conservation of momentum equation in iteration m of increment n reads (integrated form):

V_P \frac{\rho_P ^{(m-1)} - \rho _P ^{n-1}}{\Delta t} + \sum_{f}^{} \dot{m}^*_f=0,

(656)

which is probably not satisfied. Therefore, a correction of the density, pressure, velocity and mass flow is proposed in the form:

\displaystyle \rho^{(m)}	\displaystyle = \rho ^{(m-1)} + \rho '	(657)
\displaystyle p^{(m)}	\displaystyle = p ^{(m-1)} + p'	(658)
\displaystyle \boldsymbol{v}^{(m)}	\displaystyle = \boldsymbol{v}^* + \boldsymbol{v} '	(659)
\displaystyle \dot{m}_f^{(m)}	\displaystyle = \dot{m}_f^* + \dot{m}_f ' .	(660)

Assume that \boldsymbol{v},p also satisfy the conservation of momentum equation just solved, i.e.

\boldsymbol{v}_P ^{(m)} + \boldsymbol{H}_P (\boldsymbol{v} ^{(m)} ) = B_P ^{(m-1)} - D_P \nabla p ^{(m)} _P

(661)

and

\boldsymbol{v}_f ^{(m)} = \overline{\boldsymbol{v} }_f ^{(m)} + \overline{D}_f \left [ \overline{\nabla p}_f ^{(m)} - \nabla p_f ^{(m)} \right ],

(662)

leading to (subtract Equation () from Equation ())

\boldsymbol{v}_P ' + \boldsymbol{H}_P (\boldsymbol{v} ') = - D_P \nabla p'_P

(663)

and (subtract Equation () from Equation ())

\boldsymbol{v}_f' = \overline{\boldsymbol{v} }_f' + \overline{D}_f (\overline{\nabla p}_f' - \nabla p_f' ).

(664)

With the corrected values the conservation of mass reads:

V_P \frac{\rho_P ^{(m-1)} + \rho_P' - \rho_P^{n-1}}{\Delta t} + \sum_{f}^{}\dot{m}_f ^{(m)} = 0.

(665)

Now, \dot{m}_f ^{(m)} satisfies:

\displaystyle \dot{m}_f^{(m)}	\displaystyle = (\overrightarrow{\rho }_f ^{(m-1)} + \overrightarrow{\rho} _f') A_f (\boldsymbol{v}_f^* + \boldsymbol{v}_f ') \cdot \boldsymbol{n}_f	(666)
	\displaystyle \approx \dot{m}_f^* + \overrightarrow{\rho }_f ^{(m-1)} A_f \bold... ...n}_f + \overrightarrow{\rho }_f' A_f \boldsymbol{v}_f^* \cdot \boldsymbol{n}_f,	(667)

neglecting second order terms. Since

\overrightarrow{p}_f ^{(m-1)} = \overrightarrow{\rho }_f ^{(m-1)} r \overline{T} _f ^{(m-1)}

(668)

we also have (the temperature T is kept constant)

\overrightarrow{p}_f' = \overrightarrow{\rho }_f' r \overline{T} _f ^{(m-1)}

(669)

and the last term can be replaced by

\overrightarrow{\rho }_f' A_f \boldsymbol{v}_f^* \cdot \boldsymbol{n}_f = \frac{\overrightarrow{p}_f' }{r \overline{T} _f ^{(m-1)} } \frac{\dot{m}_f^*}{\overrightarrow{\rho } _f ^{(m-1)}}.

(670)

The second term in Equation () amounts to:

\displaystyle \overrightarrow{\rho }_f ^{(m-1)} A_f \boldsymbol{v}_f' \cdot \boldsymbol{n}_f	\displaystyle = \overrightarrow{\rho }_f ^{(m-1)} A_f \overline{\boldsymbol{v} }_f' \cdot \boldsymbol{n}_f
	\displaystyle + \overrightarrow{\rho }_f ^{(m-1)} A_f \overline{D}_f ( \overline{\nabla p}_f' - \nabla p_f') \cdot \boldsymbol{n}_f.	(671)

Since

\boldsymbol{v}_P' + \boldsymbol{H}_P (\boldsymbol{v}') = -D_P \nabla p'_P

(672)

we have

\boldsymbol{v}_F' + \boldsymbol{H}_F (\boldsymbol{v}') = -D_F \nabla p'_F

(673)

and (taking the mean)

\overline{\boldsymbol{v} }_f' + \overline{\boldsymbol{H} }_f (\boldsymbol{v}') = - \overline{D}_f \overline{\nabla p}_f',

(674)

and the second term now reads:

\displaystyle \overrightarrow{\rho }_f ^{(m-1)} A_f \boldsymbol{v}_f' \cdot \boldsymbol{n}_f	\displaystyle = - \overrightarrow{\rho }_f ^{(m-1)} A_f \overline{D}_f \nabla p _f' \cdot \boldsymbol{n}_f
	\displaystyle - \overrightarrow{\rho }_f ^{(m-1)} A_f \overline{\boldsymbol{H} }_f (\boldsymbol{v}') \cdot \boldsymbol{n}_f	(675)

leading to the folowing form for the conservation of mass:

	\displaystyle V_P \frac{\rho_P ^{(m-1)} + \rho_P' - \rho_P^{n-1}}{\Delta t} + \sum_{f}^{} \left \{ \dot{m}_f^* \right .
	\displaystyle - \overrightarrow{\rho }_f ^{(m-1)} A_f \overline{D}_f \nabla p_f... ...-1)} A_f \overline{\boldsymbol{H} }_f ( \boldsymbol{v}') \cdot \boldsymbol{n}_f
	\displaystyle + \left . \frac{\overrightarrow{p}_f' }{r \overline{T} _f ^{(m-1)} } \frac{\dot{m}_f^*}{\overrightarrow{\rho }_f ^{(m-1)} } \right \} = 0	(676)

	\displaystyle \frac{V_P p_P'}{r T_P ^{(m-1)} \Delta t} + \sum_{f}^{} \overright... ...rac{\dot{m}_f^*}{r \overline{T} _f ^{(m-1)} \overrightarrow{\rho }_f ^{(m-1)} }
	\displaystyle - \sum_{f}^{} \overrightarrow{\rho }_f ^{(m-1)} A_f \overline{D}_... ...m-1)} A_f \overline{\boldsymbol{H} }_f (\boldsymbol{v}') \cdot \boldsymbol{n}_f
	\displaystyle = V_P \frac{\rho_P ^{n-1} - \rho_P ^{(m-1)} }{\Delta t} - \sum_{f}^{} \dot{m}_f^*	(677)

The first term is transient-like, the second convection-like and the third diffusion-like. The fourth term containing the \boldsymbol{H}-terms is usually neglected (this corresponds to the so-called SIMPLE algorithm - Semi-Implicit Method for Pressure-Linked Equations [66]).

For the convective interpolation of the pressure correction the upwind difference scheme is always used.

The normal gradient in the diffusion term can be written as:

\nabla p_f' \cdot \boldsymbol{n}_f \approx \frac{p_F'-p_P'}{l_{PF}} + \nabla p_f' \cdot (\boldsymbol{n}_f - \boldsymbol{j}_f).

(678)

Therefore, Equation () has to be solved in an iterative way: in the first subiteration only the first term in Equation () is taken into account, in the following subiterations the first term in Equation () is based on the actual subiteration, the second term is calculated based on the p' solution from the previous subiteration. This can be expressed as follows:

(\nabla p_f' \cdot \boldsymbol{n}_f)^s \approx \frac{p_F'^s-p_P'^s}{l_{PF}} + \nabla p_f'^{s-1} \cdot (\boldsymbol{n}_f - \boldsymbol{j}_f),

(679)

where s is the subiteration number. Usually, only 2 subiterations are needed.

The solution of the conservation of mass yields p' at the element centers. Also here, underrelaxation is applied by taking only about 20 % of p' into account. This is typical for the SIMPLE algorithm and is due to the neglection of the \boldsymbol{H}-terms hereby increasing the risk of divergence. After determining p' the other corrections satisfy:

\boldsymbol{v}_P' = - D_P \nabla p_P' - \boldsymbol{H}_P(\boldsymbol{v}'),

(680)

\dot{m}_f' = \frac{\overrightarrow{p}_f' \dot{m}_f^* }{\overrightarrow{p}_f ^{(m-1)} } - \overrightarrow{\rho }_f ^{(m-1)} A_f \overline{D}_f \nabla p_f' \cdot \boldsymbol{n}_f - \overrightarrow{\rho }_f ^{(m-1)} A_f \overline{\boldsymbol{H} }_f (\boldsymbol{v}') \cdot \boldsymbol{n}_f,

(681)

(the \overline{\boldsymbol{H}}_f-terms are usually neglected)

\boldsymbol{v}_P ^{(m)} = \boldsymbol{v}_P^* + \boldsymbol{v}_P',

(682)

\dot{m}_f ^{(m)} = \dot{m}_f^* + \dot{m}_f',

(683)

and

p_P ^{(m)} = p_P ^{(m-1)} + 0.2 p_P'

(684)

\overline{p}_f ^{(m)} and \overline{\boldsymbol{v} }_f ^{(m)} are obtained from p_P ^{(m)} and \boldsymbol{v} _P ^{(m)} with the usual procedures.

Also here the boundary conditions are of utmost importance. The convection boundary conditions amount to:

Inlet: if no pressure is given (subsonic flow) the face correction is set equal to the element center correction: \overrightarrow{p}_f' \approx p_P'. If the pressure is given (supersonic flow) we have \overrightarrow{p}_f'=0.
Outlet: the same rules apply as for the inlet. If the pressure is given (subsonic flow) the correction is zero, if it is not given the element center value is taken.
wall or sliding conditions: since \dot{m}_f^*=0 the convection correction term is zero.

The diffusion boundary conditions satisfy:

Inlet and pressure given: since at the inlet the temperature is usually given, so is the density (gas equation). At an inlet also the velocity is specified, so the mass flow is fixed. Therefore, the total correction to the mass flow is zero. Since the pressure correction is zero, so is the convective mass flow correction. Therefore, the diffusive mass flow correction must also be zero.
Inlet and no pressure given: the term for the inlet face is neglected. This boils down to the assumption that the pressure correction does not change in normal inlet direction.
Outlet and pressure given: the gradient of the pressure correction is calculated between the outlet face and the element center. Since the correction at the face is zero Equation () is replaced by:

\nabla p_f' \cdot \boldsymbol{n}_f \approx \frac{-p_P'}{l_{Pf}} + \nabla p_f' \cdot (\boldsymbol{n}_f - \boldsymbol{j}_f). (685)
Outlet and no pressure given: same procedure as for the inlet without given pressure.
wall or sliding conditions: the mass flow is zero, so also the total correction is zero. Since the convective correction is zero (see above), the diffusive correction must also be zero.

In the SIMPLE algorithm the \boldsymbol{H}-terms were neglected. This leads to slower convergence and the need to introduce unterrelaxation in order to avoid divergence. Since the publication of the SIMPLE scheme by Patankar and Spalding [66] several other schemes have been proposed to improve the convergence. One of them is the SIMPLEC scheme, which stands for SIMPLE-Consistent [87]. Instead of neglecting the \boldsymbol{H}-terms it is assumed that the velocity correction in P is the weighted mean of the correction in its neighbors. The weighting coefficients are taken from the momentum equation. Mathematically this can be written as:

\boldsymbol{v}_P' \approx \frac{\sum_{F}^{} a_F \boldsymbol{v}_F' }{\sum_{F}^{} a_F }.

(686)

Now, the left hand side of Equation () can be written as:

\displaystyle \boldsymbol{v}_P' + \boldsymbol{H}_P (\boldsymbol{v}')	\displaystyle = \boldsymbol{v}_P' + \sum_{F}^{} \frac{a_F \boldsymbol{v}_F' }{a_P}
	\displaystyle = \boldsymbol{v}_P' + \boldsymbol{v}_P' \frac{\sum_{F}^{} a_F }{a_P}
	\displaystyle = \boldsymbol{v}_P' + \boldsymbol{v}_P' H_P(1)
	\displaystyle = \boldsymbol{v}_P' (1 + H_P(1)).	(687)

Now, Equation () can be replaced by:

\boldsymbol{v}_P' = -D_P^c \nabla p'_P,

(688)

where

D_P^c = \frac{D_P}{1+H_P(1)}.

(689)

In a similar way, Equation () amounts to:

\overline{\boldsymbol{v} }_f' = - \overline{D}_f^c \overline{\nabla p}_f',

(690)

and Equation () is replaced by:

\overrightarrow{\rho }_f ^{(m-1)} A_f \boldsymbol{v}_f' \cdot \boldsymbol{n}_f = - \overrightarrow{\rho }_f ^{(m-1)} A_f \overline{D}_f^c \nabla p _f' \cdot \boldsymbol{n}_f .

(691)

Therefore, in the SIMPLEC scheme, the fourth term in () is taken into account by replacing \overline{D}_f in the third term by \overline{D}_f^c. The SIMPLEC scheme usually converges faster, especially since no underrelaxation is needed any more.