[Rachit]

Hi everyone,

I have been trying to write a MATLAB code for 2D flow of the SIMPLE algorithm for incompressible flows on a collocated grid. I do understand that the Rhie Chow interpolation is used against the normal interpolation to calculate the velocity on the faces, in order to avoid checkerboarding effect. However, from what I gathered from the literature, there is no reference to where is the Rhie Chow interpolation used. In the calculation face velocities in the discretization of the continuity equation? Or is it used to calculate the face velocities in the calculation of convective flux in the momentum equation as well? If that is so then what is the point of using schemes like Upwind, Central Difference, QUICK, etc?

I would be extremely grateful, if someone could help me out on this.

Thanks in advance,
Rachit

[dsanthosh2]

you have to use Rhie and Chow interpolation in the continuity equation. Not in the momemtume equation.

[RodriguezFatz]

Originally Posted by dsanthosh2

you have to use Rhie and Chow interpolation in the continuity equation. Not in the momemtume equation.

But how do you get the flux at the cell faces for the convective term in the momentum euqation?

Originally Posted by Rachit

Or is it used to calculate the face velocities in the calculation of convective flux in the momentum equation as well? If that is so then what is the point of using schemes like Upwind, Central Difference, QUICK, etc?

I think it is used in both cases, continuity and momentum. Remember, that you don't really directly solve the continuity equation in schemes like SIMPLE...
The momentum term is not linear. Thus, you have to interpolate the flux on the cell face (as the "transporting" medium) and you have to interpolate the velocity on the cell faces (as the "independent variable"). If you have "rho*u*u", then "rho*u" is the flux through the face, and one single "u" is the independent variable. That's the one, you create your linear system of equations for. And this is also the one, you can decide the way of interpolation for (as you described).

If you have "rho*v*u" and you are creating the momentum equation for "u", then "rho*v" is the flux through the face, and "u" is the one you solve with your matrix...

Did that become clearer?

[andy_]

Originally Posted by Rachit

Hi everyone,

I have been trying to write a MATLAB code for 2D flow of the SIMPLE algorithm for incompressible flows on a collocated grid. I do understand that the Rhie Chow interpolation is used against the normal interpolation to calculate the velocity on the faces, in order to avoid checkerboarding effect. However, from what I gathered from the literature, there is no reference to where is the Rhie Chow interpolation used. In the calculation face velocities in the discretization of the continuity equation? Or is it used to calculate the face velocities in the calculation of convective flux in the momentum equation as well? If that is so then what is the point of using schemes like Upwind, Central Difference, QUICK, etc?

I would be extremely grateful, if someone could help me out on this.

Thanks in advance,
Rachit

Calling the addition of pressure gradient terms to the mass flux "interpolation" is rather misleading because there are no pressure terms in the mass flux. Nonetheless adding pressure smoothing terms to the mass flux is a requirement for many schemes.

The pressure smoothing terms are usually needed in the mass fluxes in the continuity equation to oppose pressure/velocity decoupling. For the mass fluxes in the other transport equations the pressure smoothing terms are optional but are usually included so that what is a mass flux is consistent and various numerical fixes that rely on the mass fluxes in a cell summing to zero hold.

The flux of momentum involves squared velocities. One velocity can be involved in determining the mass flux and the other velocity that is being transported by the mass flux can be evaluated using upwind, central difference, QUICK or whatever.

[RodriguezFatz]

Originally Posted by andy_

One velocity can be involved in determining the mass flux and the other velocity that is being transported by the mass flux can be evaluated using upwind, central difference, QUICK or whatever.

That's what I was trying to explain. Better read andy's post - it's easier to understand.

[andy_]

Originally Posted by RodriguezFatz

That's what I was trying to explain. Better read andy's post - it's easier to understand.

Had I seen that you had posted I would not have replied. My technique of putting posts to reply to in tabs and then getting on with other things while replying perhaps needs modifying.

[Rachit]

Dear Santo, Rodriguez and Andy,

Thank you for your replies. I shall work it out this way in my code. Thanks once again.

Regards,
Rachit

[Rachit]

Hi,

Looks like my confusion pertaining the Rhie Chow interpolation hasn't ended yet. From what's there in Versteeg and Malalasakera, to calculate the velocity at face e, I will need to the pressure at EE. What if my EE is outside the domain? Is there any way find the value of PEE without using the concept of ghost cells?

Also, to calculate the pressure smoothing term I need to aP which is the coefficient of variable u at point P. However, the coefficient is itself made up of convective flux terms (something for which I am using the Rhie Chow interpolation to calculate). What do I do then? Do I take the flux values from the previous step or one has to do something else?

Thanks in advance,
Rachit

[andy_]

> What if my EE is outside the domain?

What to do about the pressure smoothing on the boundary is up to you. For example, on a solid boundary the mass flow is zero but the pressure gradients will generally be non zero. So do you use a zero mass flux or include the pressure smoothing term? If the former, global mass conservation will be physically reasonable but there will be a jump caused by smoothing on one face but not the other and in the presence of strong pressure gradients from body forces like swirl it can be large and unreasonable. If the latter, you will have unreasonable mass fluxes on the solution boundary and issues such as coming up with ways to evaluate pressure gradients at the boundary.

> Do I take the flux values from the previous step or one has to do something else?

Storing the mass fluxes along with the solution variables is a common thing to do for schemes using Rhie and Chow type pressure smoothing.

[FMDenaro]

Justo to give you my opinion ...

First, start considering everywhere in the domain the Hodge decomposition

Vn+1 + Grad phi = V*

1) in the interior, after computing V* somehow, you will search for a "pressure" gradient such that it ensures Div Vn+1 = 0 -> Div Grad phi = Div V*
2) on the boundaries, you only have to prescribe something about Vn+1 and use the Hodge decomposition to ensures the correct mass flux by solving the modified pressure equation.

On co-located grids, you can use the approximate projection method (APM) to avoid checkerboard modes in the pressure, the RC interpolation can be applied on the flux defined by the Hodge decomposition.
However, the RC interpolation has some degrading in the accuracy...

If you want, some more deatils about possible remedy for checkerboard modes are addressed at:
http://onlinelibrary.wiley.com/doi/1....1368/abstract
http://onlinelibrary.wiley.com/doi/1....1368/abstract

[malaboss]

Hi,
I am actually also trying to understand how is done the Rhie Chow interpolation in OpenFOAM and there are several things that still remain unclear. I understand that we interpolate pressure on the faces of the cell and then use this face pressure to construct the central difference scheme for pressure, which couples the adjacent pressures and avoid the checkerboard effect. This is refered as the pressure smoothing by Philipp and Andy.

1) I don't get were it is formulated in OpenFOAM. I see some descriptions of the RC interpolation with the usage of the H and A Matrix, but I don't really understand the link between theses matrices and the actual RC interpolation.

2) In the literature, some improvements of the RC interpolation are mentioned. Then, I am not sure which version has been implemented in OpenFOAM. Do you have some info about it ?

3) In order to better understand the effect of this anti-chekboarding interpolation, do you know if it is possible to disable the Rhie Chow interpolation in the OpenFOAM and how to do it ?

Thanks in advance for your time and for your help.

[sbaffini]

In my opinion, a very detailed and clear set of notes on unstructured FV methods for CFD is the one by Prof. Murthy, former Fluent developer:

https://engineering.purdue.edu/ME608...ass-notes.html

where you can find most of the things you need about a CFD solver.

For what concerns OpenFOAM specifically, consider the following article:

Tukovic, Jasak: A moving mesh finite volume interface tracking method for surface tension dominated interfacial fluid flow. Comp. Fluids 55 (2012), pp. 70-84

Here, the OpenFOAM implementation and modifications are discussed at length.

However, consider that most of RC has nothing to do with the pressure equation, where it simply reduces to use an Approximate Projection Method (APM), as discussed by F.M. Denaro above. More roughly speaking this just means using a compact, straight, Laplacian operator in the pressure equation. Moreover, this is actually your only option for unstructured FV codes.

Where RC really comes into play is in the determination of the mass flux, which in order to be consistent with the above pressure equation has to include a pressure part. Still, this part is only a perturbation vanishing with the grid spacing going to 0.

[malaboss]

Hi and thanks for the answer.
I found in Tukovic and Jasak's paper what I was looking for.
I will sum it up to make sure that I really get it and to save time to other people interested in it :

OpenFOAM uses collocated grid in order to handle complex geometries. One known problem about the collocated grids is the checkboarder effect, that is to say, the decoupling between adjacent cells for pressure and velocity which will not be able to damp possible oscillations of the fields, leading to an unphysical result.
So why do we have a decoupling between adjacent cells ?

Let's consider the 3 same cubic cells (W P and E separated by faces which locations are denoted by w and e) and a 1D case, and a constant density case.
We do not want to calculate grad(P) with a 1st order scheme but a central difference scheme.

Momentum equation applied to the 3 cells

u_P=(H/A)_P - grad(p)_P/A_P
u_W=(H/A)_W - grad(p)_W/A_W
u_E=(H/A)_E - grad(p)_E/A_E

I do not discretize grad(p) for the moment because this is where lies the RC interpolation

Continuity equation with FVM applied to the cell P.

rho ue Se - rho uw Sw=0
We need to know the velocity at the faces e and w, but we have it at the center.
So we approximate ue by (uE+uP)/2 and uw by (uW+uP)/2
As a result
ue=1/2 * [(H/A)_E + (H/A)_P - grad(p)_P/A_P - grad(p)_E/A_E ]
uw=1/2 * [(H/A)_W + (H/A)_P - grad(p)_P/A_P - grad(p)_W/A_W ]
And the continuity equation reduces to
1/2 * [(H/A)_E + (H/A)_P - grad(p)_P/A_P - grad(p)_E/A_E ] - 1/2 * [(H/A)_W + (H/A)_P - grad(p)_P/A_P - grad(p)_W/A_W ] = 0

This is an approximation for

1/2 * [(H/A)_e - grad(p)_e/A_e ] - 1/2 * [(H/A)_w - grad(p)_w/A_w ] = 0

The goal of the RC interpolation is to involve Pressure at adjacent cells in the continuity equation.

What is not RC interpolation
Calculate grap(P) at cell centers or neighboring cells first.
Then the continuity equations without the H/A terms reduces to
1/2 * [- grad(p)_P - grad(p)_E ] - 1/2 * [-grad(p)_P - grad(p)_W ]
= 1/(4*deltaX) * [- (P_E-P_W)- (P_P-P_EE) +(P_E-P_W) +(P_P-P_WW) ]
= 1/(4*deltaX) * [P_EE -P_WW ]
Here only P_EE and P_WW are coupled leading to unphysical oscillations

What is RC interpolation
Calculate grap(P) at cell faces directly
Then the continuity equations without the H/A terms reduces to
1/2 * [- grad(p)_e ] - 1/2 * [-grad(p)_w ]
= 1/(4*deltaX) * [- (P_P-P_E) - (P_P-P_W)]
= 1/(4*deltaX) * [P_E - 2P_P +P_W ]
Here, adjacent cells are coupled !

Remark
In the continuity equation, the term (H/A)_P will also disappear. I don't know which treatment is done by RC interpolation for this.
1/2 * [(H/A)_E + (H/A)_P - grad(p)_P/A_P - grad(p)_E/A_E ] - 1/2 * [(H/A)_W + (H/A)_P - grad(p)_P/A_P - grad(p)_W/A_W ] = 0

Answer to question 1 : How is formulated the RC interpolation in OpenFOAM ?
The particular gradient calculation at the faces for p is done in the formulation of the laplacian of p which calculates a flux of gradient p, in the pressure correction equation. This laplacian of p is equivalent to div(snGrad(p)). snGrad(P) is the Rhie Chow interpolation

Answer to question 2 : Has the last improvement of RC interpolation been implemented in OpenFOAM ?
The basic Rhie chow interpolation seems implemented this way. If the implementation of rhie chow interpolation is really what I described then there is room for improvement.

Answer to question 3 : Can we remove the Rhie Chow interpolation from OpenFOAM ?
Yes, we can remove RC interpolation from OpenFOAM by just replacing laplacian(p) by div(grad(p))

Hope that helps someone !

[sbaffini]

I am not an expert in OpenFOAM, but i would consider the following two points:

1) Improvements in the RC procedure, as i understand them, have been around for a while and i would be impressed if the OpenFOAM formultation would still be based on the original one. I'm referring here to the time-step/under-relaxation dependence. I don't know of other relevant improvements. These are also clearly mentioned in the Tukovic paper.

2) OpenFOAM works with classes for implicit (fvm) and explicit (fvc) operators which are clearly distinguished. With the former you can build up a discretization matrix, with the latter you can only work on the source term side. This is fundamental to understand and also clarifies what approach you can actually use for the Poisson equation.

In OpenFOAM, according to my knowledge, you can only use the laplacian as an impicit discretization and not the extended one. There are different reasons for this, like the inherent difficulty to track the extended stencil for each cell, the possible bad properties of such a wide band matrix, the nearly impossible task to obtain the discretization coefficients for such a discrete operator. I am not aware of any unstructured code doing differently.

When looked in this perspective, RC is just a justification for such a compact Laplacian and the resulting derivation of the correct mass fluxes to use in the remaining equations.