2.3 Rhie - Chow Interpolation Method

  

Figure 2.3.1: One Dimensional Control Volume

If the one dimensional x momentum equation is considered

$$\dxdt{(\rho\,u)}{t} + \dxdt{(\rho u u)}{x} \; = \; \dxdt{\;}{x} \left( \mu \dxdt{u}{x} \right) - \dxdt{p}{x}$$ (2.3.1)

where u is the x component of velocity, then the pressure derivative term must be discretised over the control volume. If the control volume considered contains the node P with neighbouring nodes W and E as in Figure 2.3.1 then the contribution to the discretised equation from this term is p_e - p_w. Assuming that the faces e and w are midway between the nodes either side of them, then using linear interpolation for the face values of pressure

$$p_w - p_e \; = \; \frac{p_W + p_P}{2} - \frac{p_P + p_E}{2} \; = \; \frac{p_W + p_E}{2} \nomsort{p}{$p$ }{Pressure}$$ (2.3.2)

This indicates that the pressure differential term in the momentum equation leads to a discretised form which is a relationship between alternate nodal pressure values rather than adjacent ones. As a consequence the pressure differential is calculated over a coarser mesh than the one used for all other quantities. It also means that an oscillatory pressure field, in which, for example, the pressure at consecutive grid points follows the sequence 1,100,1,100,1 etc, would be treated as a uniform pressure field as alternate grid points have the same value. In higher dimensions the problem leads to the checkerboard effect where in each of the dimensions alternate grid points have the same pressure value but not adjacent nodes. This means that in three dimensions the pressure field could take 8 different values but still the discretised pressure differential term would give a zero pressure differential. Given any pressure field this checkerboard field can be added to it without changing the pressure gradient values. Thus if p is the correct pressure field and u the correct velocity field which together satisfy the discretised equations, then so does $p + \lambda \, p_c$ and u where p_c is the checkerboard pressure field. Any process has only one correct pressure field at any time. This means that any pressure obtained as a result of solving the momentum and continuity equations may predict the correct velocity whilst having an error in the pressure field equal to $\lambda\,p_c$.

Consider now the steady one dimensional continuity equation for an incompressible flow, which is

$$\frac{du}{dx} \; = \; 0$$ (2.3.3)

Integrating this over the control volume of Figure 2.3.1, and assuming that the faces lie midway between their nodes, the following equation is obtained
 
The problems of an equation which tests alternate nodal values, rather than adjacent ones, is again obtained. Hence, similar to the pressure, a checkerboard velocity field would satisfy the continuity equation.

The standard mechanism for handling the above problem is the use of a staggered mesh. The staggered grid for the velocity components was first suggested by Harlow and Welch[55] and forms the basis of the SIMPLE procedure[14] and of the SIVA procedure[13]. In the staggered grid the u velocity is calculated at the faces of the control volume which have their normal in the x coordinate direction. Similarly the v and w velocities are calculated at the faces of the control volume with normals in the y and z directions respectively. The calculation of the velocities at the control volume faces means that in equation (2.3.6) the values of u_e and u_w are obtained directly from the solution of the u component of velocity. This means that the discretised continuity equation now consists of a relationship between adjacent u values and so the checkerboard velocity pattern is no longer possible. Similarly the staggered grid leads to the pressure differential between adjacent nodes becoming the driving force of the velocity component solved for between these nodes.

This staggering of the velocity components is of huge benefit on structured meshes but the method does not extend easily to unstructured meshes. In the approach presented here all quantities are solved and stored at the element centroid. The face values of the velocity components have to be calculated from these element based values. This leads to the need to employ an alternative interpolation method, to those described earlier in this section, which does not suffer from the checkerboard effect. The Rhie - Chow interpolation method[44] offers one approach which satisfies these requirements.

The u momentum equation

$$\nfour \nfour \nfour \nfour \frac{\displaystyle \partial (\... ... u ) = -\dxdt{p}{x} + \text{div}(\mu \, \text{grad}(u) ) + S$$ (2.3.4)

can be discretised, using the techniques described earlier in this chapter, over the control volume about a node P to produce an equation which can be written in the form

$$a_P u_P + (\nabla_x p)_P \; = \; ( \sum a_{nb} u_{nb} )_P + S_P$$ (2.3.5)

where $\nabla_x p$ is the discretised contribution from the pressure gradient term. Similarly for the adjacent node A

$$a_A u_A + (\nabla_x p)_A \; = \; ( \sum a_{nb} u_{nb} )_A + S_A$$ (2.3.6)

From the conservation principle of the control volume formulation the u velocity at a point on the face between the nodes must also have a discretised momentum equation of the form

$$a_f u_f + (\nabla_x p)_f \; = \; ( \sum a_{nb} u_{nb} )_f + S_f$$ (2.3.7)

The Rhie - Chow interpolation method uses the equations (2.3.8) and (2.3.9) to approximate a solution of (2.3.10). It is assumed that the right hand side of (2.3.10) may be approximated by using a weighted linear interpolation of the corresponding terms in (2.3.8) and (2.3.9). Thus

$$\nfour \nfour \nfour a_f u_f + (\nabla_x p)_f \; = \; \overl... ... + S_f} \; = \; \overline{a_f u_f} + \overline{(\nabla_x p)_f}$$ (2.3.8)

where the overline in the above equation indicates a weighted linear interpolation of the variable. Assuming that $a_f \approx \overline{a_f}$ then

$$u_f \; = \; \overline{u_f} + \overline{d_f} ( \overline{\nabla_x p_f} - \nabla_x p_f )$$ (2.3.9)

where if $\alpha$ is the weighting factor used

$$\begin{split} \overline{u_f} \; &= \; \alpha u_P + ( 1-\alph... ...\alpha ) a_A \\ \overline{d_f} \; &= \; a^{-1}_f \end{split}$$ (2.3.10)

All that is now required to complete the Rhie - Chow interpolation method is to select the weighting factor $\alpha$. The obvious choice is that $\alpha$ should be equal to the distance of the node in element A to the face divided by the distance from the same node to the node in element P. The problem with this choice was highlighted when solving two dimensional flow in a duct. A very fine mesh was used near the inflow boundary in an attempt to accurately model the formation of the flow profile. In the region beyond this fine mesh a much coarser mesh was used, see Figure 2.3.2.

  

Figure 2.3.2: Mesh Used for Rhie - Chow Weighting Test Case

In this region the velocity on the centre line of the duct does not change and the pressure decreases at a constant rate. Using weighted linear interpolation leads to the centre line values of the velocity component along the duct shown in Figure 2.3.3. As a comparison a plot from Flow3d[56], Version 2.3.2, is also shown on the graph. Flow3d, like CWNN , uses Rhie and Chow interpolation and collocated variables but, unlike CWNN , it uses weighted averaging in the calculation of the convection fluxes.

  

Figure 2.3.3: Rhie - Chow Weighting Test Case, Velocity Along Duct

It can be seen that immediately upon changing to the coarser mesh the value of the velocity exhibits a physically unrealistic decrease. Examining the terms of the Rhie - Chow interpolation method for this case, assuming the decrease takes place after the profile has been fully formed, should give u_A = u_P = u_f. This means that the correction terms should be equal to zero. As the pressure gradient is constant, C, at this point

$$\begin{split} \left( \, \nabla_x \, p\,\right)_P \; & = \; V... ...ft( \, \nabla_x \, p\,\right)_A \; & = \; V_A \, C \end{split}$$ (2.3.11)

where V_P is the volume of element P. Hence

$$\overline{\left( \, \nabla_x \, p\,\right)_f} \; = \; C\, \left\{ \alpha\, V_P + (1-\alpha)\,V_A \right\}$$ (2.3.12)

The second correction term gives

$$\begin{split} \left( \, \nabla_x \, p\,\right)_f \; &= \; (p... ...\, A_f \\ &= \; (p_P + d_{AP}\, C - p_P ) \, A_f \end{split}$$ (2.3.13)

where d_AP is the distance between the nodes in elements A and P. As these nodes lie at the centroids of the elements, for the mesh used in this simulation

$$d_{AP} \; = \; \frac{ 0.5\,V_P + 0.5\, V_A}{A_f}$$ (2.3.14)

and so

$$\left( \, \nabla_x \, p\,\right)_f \; = \; C\, \left( \frac{V_P}{2}+\frac{V_A}{2} \right)$$ (2.3.15)

Unless the value of $\alpha$ is equal to 0.5 there will be a non zero contribution from the correction term, which is equal to equation (2.3.15) minus equation (2.3.18). This shows that the weighting factor used in the pressure gradient correction term should be 0.5 rather than distance weighted interpolation. For consistency the value of the weighting factor is set to 0.5 in all terms. This choice needs to examined further to ascertain the validity of the use of the factor 0.5 in each term. A search of relevant literature found only one reference to the choice of weighting factors in Rhie - Chow interpolation. In the CFX4[57], Version 4.1, manual it is stated that the code uses weighted interpolation on the velocity and 0.5 factors in the pressure gradient terms.