Rhie-Chow interpolation

Before continuing with the conservation of mass equation a new type of interpolation has to be introduced: the Rhie-Chow interpolation. So far we encountered for the primary variables such as velocity, static pressure or static temperature at internal faces the mean interpolation (e.g. \overline{\boldsymbol{v}_f } ) and the convective interpolation (e.g. \overrightarrow{\boldsymbol{v}_f }) . For external faces the extrapolation is dictated by the kind of boundary condition (e.g. outflow: constant extrapolation or wall/sliding condition: linear extrapolation) and these facial values are denoted without any extra sign.

The need for the Rhie-Chow interpolation comes from the fact that for a long time CFD calculations on collocated grids (i.e. grids for which the velocity, pressure, temperature etc. are unknown at the same positions) were not successfull due to the occurrence of checkerboard-type solutions. Therefore, people resorted to staggered grids, in whch the pressure was calculated at positions halfway the velocity positions [61]. Especially in three dimensions this leads to hopelessly complicated data structures. The use of the Rhie-Chow interpolation, however, solved this problem.

The Rhie-Chow interpolation (only defined on internal faces) starts from a velocity/pressure field satisfying the conservation of momentum equation. From the last section this equation can be written as follows at element center P:

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

The terms correspond to:

  1. A velocity term for element P at the current iteration. Notice that the equation has been divided by the diagonal coefficient a_P of the left hand side of the equation system
  2. \boldsymbol{H}_P(\boldsymbol{v}^{(m)}): velocity contributions from the element neighbors at the current iteration. Here \boldsymbol{H} is to be understood as a linear function.
  3. \boldsymbol{B}_P^{(m-1)}: the right hand side terms calculated in the previous iterations
  4. D_P \nabla p^{(m-1)}_P: the pressure gradient; D_P corresponds to V_P/a_P where V_P is the volume of element P. Notice that Gauss' theorem was not applied to this term.

More generically, this equation can be written for a velocity field \boldsymbol{v} and pressure field p satisfying the momentum equations as follows:

\boldsymbol{v}_P + \boldsymbol{H}_P(\boldsymbol{v}) = \boldsymbol{B}_P - D_P \nabla p_P (636)

and for element F:

\boldsymbol{v}_F + \boldsymbol{H}_F(\boldsymbol{v}) = \boldsymbol{B}_F - D_F \nabla p_F. (637)

Now we are looking for the velocity v_f at the face f satisfying the momentum equation as well, i.e.

\boldsymbol{v}_f + \boldsymbol{H}_f(\boldsymbol{v}) = \boldsymbol{B}_f - D_f \nabla p_f, (638)

where

\displaystyle \boldsymbol{H}_f \displaystyle = \overline{\boldsymbol{H} }_f = (\boldsymbol{H}_P + \boldsymbol{H}_F)/2 (639)
\displaystyle \boldsymbol{B}_f \displaystyle = \overline{\boldsymbol{B} }_f = (\boldsymbol{B}_P + \boldsymbol{B}_F)/2 (640)
\displaystyle D_f \displaystyle = \overline{D }_f = (D_P + D_F)/2. (641)

Hence:

\boldsymbol{v}_f + \overline{\boldsymbol{H} }_f (\boldsymbol{v}) = \overline{\boldsymbol{B} }_f - \overline{D}_f \nabla p _f. (642)

On the other hand. one can also just take the mean of equation ([*]) and ([*]) resulting in:

\overline{\boldsymbol{v}}_f + \overline{\boldsymbol{H} }_f (\boldsymbol{v}) = \overline{\boldsymbol{B} }_f - \frac{1}{2}(D_P \nabla p_P + D_F \nabla p_F), (643)

which is within second order accurary equivalent to [61]:

\overline{\boldsymbol{v}}_f + \overline{\boldsymbol{H} }_f (\boldsymbol{v}) = \overline{\boldsymbol{B} }_f - \overline{D}_f \overline{\nabla p}_f. (644)

Subtracting Equation ([*]) from Equation ([*]) yields:

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

Recall that (Equation ([*])) \nabla p_f is the corrected form of \overline{\nabla p}_f obtained by enforcing the short range gradient based on elements P and F, i.e.

\nabla p_f \cdot \boldsymbol{e}_{PF} = (p_F-p_P)/l_{PF}. (646)

This is done by defining

\nabla p_f = \overline{\nabla p}_f + \left [ \frac{p_F-p_P}{l_{PF}} - \overline{\nabla p}_f \cdot \boldsymbol{j}_f \right ] \boldsymbol{j}_f. (647)

Substituting Equation ([*]) into Equation ([*]) yields

\boldsymbol{v}_f = \overline{\boldsymbol{v} }_f - \overline{D}_f \left [ \frac{p_F-p_P}{l_{PF}} - \overline{\nabla p}_f \cdot \boldsymbol{j}_f \right ] \boldsymbol{j}_f. (648)

Whereas \overline{\boldsymbol{v} }_f is obtained through mean interpolation \boldsymbol{v}_f is obtained by Rhie-Chow interpolation. It is a kind of improved value for the velocity by enforcing the short-range correctness of the pressure gradient.