Pressure Correction Equation

The ability (of one) to solve eqs. (4.292) – (4.294) to obtain the velocity field relies on the accuracy of the pressure field, which is unknown a priori. If we assume that the pressure field is p* – which can be obtained from the previous time step, iteration, or guessing – the corresponding velocity field can be obtained by solving eqs. (4.292) – (4.294).

 $a_{e}u_{e}^{*}=\sum{a_{nb}u_{nb}^{*}+b+A_{e}(p_{P}^{*}-p_{E}^{*})}$ (1)
 $a_{n}v_{n}^{*}=\sum{a_{nb}v_{nb}^{*}+b+A_{n}(p_{P}^{*}-p_{N}^{*})}$ (2)
 $a_{t}w_{t}^{*}=\sum{a_{nb}w_{nb}^{*}+b+A_{t}(p_{P}^{*}-p_{T}^{*})}$ (3)

Since the velocity field obtained from eqs. (4.295) – (4.297) is based on a guessed pressure field, it may not satisfy the continuity equation. An algorithm based on pressure correction (Patankar, 1980) can be used to obtain the correct pressure field and consequently, the velocity field. If the correct pressure can be obtained by

 $\begin{matrix}{}\\\end{matrix}p=p^{*}+{p}'$ (4)

where p' is pressure correction, the new velocity field based on the corrected pressure field will be

 $\begin{matrix}{}\\\end{matrix}u=u^{*}+{u}',\text{ }v=v^{*}+{v}',\text{ }w=w^{*}+{w}'$ (5)

where u', v', and w' are velocity corrections. The corrected velocity field satisfies eqs(4.292) – (4.294), i.e.,

$a_{e}(u_{e}^{*}+{u}'_{e})=\sum{a_{nb}(u_{nb}^{*}+{u}'_{nb})+b+A_{e}[(p_{P}^{*}+{p}'_{P})-(p_{E}^{*}+{p}'_{E})}]$

$a_{n}(v_{n}^{*}+{v}'_{n})=\sum{a_{nb}(v_{nb}^{*}+{v}'_{nb})+b+A_{n}[(p_{P}^{*}+{p}'_{P})-(p_{N}^{*}+{p}'_{N})}]$

$a_{t}(w_{t}^{*}+{w}'_{t})=\sum{a_{nb}(w_{nb}^{*}+{w}'_{nb})+b+A_{t}[(p_{P}^{*}+{p}'_{P})-(p_{T}^{*}+{p}'_{T})}]$

Subtracting eqs. (4.295) – (4.297) from the above three equations yields the following equations for the velocity corrections:

 $a_{e}{u}'_{e}=\sum{a_{nb}{u}'_{nb}+A_{e}({p}'_{P}-{p}'_{E})}$ (6)
 $a_{n}{v}'_{n}=\sum{a_{nb}{v}'_{nb}+A_{n}({p}'_{P}-{p}'_{N})}$ (7)
 $a_{t}{w}'_{t}=\sum{a_{nb}{w}'_{nb}+A_{t}({p}'_{P}-{p}'_{T})}$ (8)

It can be seen from eqs. (4.300) – (4.302) that the correction of the velocity can be from (1) the difference on pressure correction at two neighbor points in the main grid (the second terms), and (2) velocity correction from all neighbor points (the first terms). The solution of the velocity correction from eqs. (4.300) – (4.302) will be very complicated because the velocity correction at neighbor points will be related to pressure corrections at more grid points, and eventually velocity correction at any point will be related to the pressure corrections in the entire computational domain. Assuming the contribution of the latter is negligible, which is equivalent to assuming anb = 0 in eqs. (4.300) – (4.302), the velocity correction becomes

 $\begin{matrix}{}\\\end{matrix}{u}'_{e}=d_{e}({p}'_{P}-{p}'_{E}),\text{ }{v}'_{n}=d_{n}({p}'_{P}-{p}'_{N}),\text{ }{w}'_{t}=d_{t}({p}'_{P}-{p}'_{T})$ (9)

where

 $\begin{matrix}{}\\\end{matrix}d_{e}=A_{e}/a_{e},\text{ }d_{n}=A_{n}/a_{n},\text{ }d_{t}=A_{t}/a_{t}$ (10)

Substituting eq. (4.303) into eq. (4.299), we have

 $u_{e}=u_{e}^{*}+d_{e}({p}'_{P}-{p}'_{E})$ (11)
 $v_{n}=v_{n}^{*}+d_{n}({p}'_{P}-{p}'_{N})$ (12)
 $w_{t}=w_{t}^{*}+d_{t}({p}'_{P}-{p}'_{T})$ (13)

which can be used to obtain the corrected velocity field from the pressure correction. The assumption of neglecting the first terms in eqs. (4.300) – (4.302) significantly simplifies the equations in order to obtain the corrected velocity field. If, however, this simplification results in unacceptable error in the final velocity field, it may not be used. As the iteration progresses, v * , gradually approaches the true velocity, the velocity correction at each point approaches zero. After the converged velocity is obtained, the velocity correction from all the neighbor points will be zero. Therefore, neglecting the first terms in eqs. (4.300) – (4.302) will not affect the ultimate velocity field. It will, however, affect how quickly a converged solution can be obtained. We now turn our attention to the algebraic equation for the pressure correction. Since the starred velocity field (u*, v*, and w*) does not satisfy the continuity equation, we hope that the corrected velocity will satisfy the continuity equation. For a three-dimensional flow, the continuity equation is

 $\frac{\partial \rho }{\partial t}+\frac{\partial (\rho u)}{\partial x}+\frac{\partial (\rho v)}{\partial y}+\frac{\partial (\rho w)}{\partial z}=0$ (14)

Integrating eq. (4.308) over the main control volume [see Fig. 4.23(a) for the two-dimensional projection of the three-dimensional control volume] and using a fully-implicit scheme, we obtain

\begin{align} & \frac{\rho _{P}-\rho _{P}^{0}}{\Delta t}\Delta x\Delta y\Delta z+[(\rho u)_{e}-(\rho u)_{w}]\Delta y\Delta z \\ & +[(\rho v)_{n}-(\rho v)_{s}]\Delta x\Delta z+[(\rho w)_{t}-(\rho w)_{b}]\Delta x\Delta y=0 \\ \end{align}

Substituting eqs. (4.305) – (4.307) into the above equation, the following algebraic equation for pressure correction is obtained:

 $\begin{matrix}{}\\\end{matrix}a_{P}{p}'_{P}=a_{E}{p}'_{E}+a_{W}{p}'_{W}+a_{N}{p}'_{N}+a_{S}{p}'_{S}+a_{T}{p}'_{T}+a_{B}{p}'_{B}+b$ (15)

where

 $\begin{matrix}{}\\\end{matrix}a_{E}=\rho _{e}d_{e}\Delta y\Delta z,\text{ }a_{W}=\rho _{w}d_{w}\Delta y\Delta z$ (16)
 $\begin{matrix}{}\\\end{matrix}a_{N}=\rho _{n}d_{n}\Delta z\Delta x,\text{ }a_{S}=\rho _{s}d_{s}\Delta z\Delta x$ (17)
 $\begin{matrix}{}\\\end{matrix}a_{T}=\rho _{t}d_{t}\Delta x\Delta y,\text{ }a_{B}=\rho _{b}d_{b}\Delta x\Delta y$ (18)
 $\begin{matrix}{}\\\end{matrix}a_{P}=a_{E}+a_{W}+a_{N}+a_{S}+a_{T}+a_{B}$ (19)
 \begin{align} & b=\frac{\rho _{P}^{0}-\rho _{P}}{\Delta t}\Delta x\Delta y\Delta z+[(\rho u^{*})_{w}-(\rho u^{*})_{e}]\Delta y\Delta z \\ & +[(\rho v^{*})_{s}-(\rho v^{*})_{n}]\Delta x\Delta z+[(\rho w^{*})_{b}-(\rho w^{*})_{t}]\Delta x\Delta y \\ \end{align} (20)
Figure 1: Control volume for continuity equation at boundary.

If the starred velocity terms (u*, v*, and w*) already satisfy the continuity equation, we have b = 0 from eq. (4.314). Therefore, b defined in eq. (4.314) can be viewed as the negative value of the “residual” mass for grid point P. If the continuity equation is satisfied, such residual mass should be zero. Practically, the maximum absolute value of b in the entire computational domain as well as the absolute value of the summation of b in the entire domain are used as criteria for the convergence. In order to apply eq(4.309) to obtain the pressure correction, the boundary conditions for pressure correction have to be specified. If the pressure at the boundary is given, no correction for the pressure at the boundary is needed so p' = 0 at the boundary. If the boundary velocity is given (see Fig. 1), the velocity correction for the boundary velocity (ue in Fig. 1) will be zero. It follows from eq. (4.303) that de for the control volume shown in Fig. 1 has to be zero. Therefore, one can set aE = 0 [see eq. (4.310)] if ue in Fig. 1 is given.