[mercator]

Hi community,

as the title says, I am trying to understand Rhie-Chow interpolation, and I am unfortunately quite confused and in need of some advice.

(I have read http://www.cfd-online.com/Forums/mai...rpolation.html, but I am still not sure how it is actually implemented.)

I realize these are very basic questions, sorry for this...

1. If using SIMPLE, which step is actually modified by RC? Normally I do:
a) Solve convection-diffusion to obtain u^*.
b) Calculate divergence of u^*, and solve pressure Poisson eqn.
c) Calculate u^n+1 by adding pressure gradients.
My understanding is that RC *only* affects (b). Is that correct or not?

2. To derive the actual matrix used in (b), I use the same approach as for
SIMPLE (assume rho = delta_t = 1 here):
- start with u^n+1 = u^* - nabla P
- use continuity eqn. on both sides: div u^n+1 = div u^* - div nabla P
- require div u^n+1 = 0 => div nabla P = div u^*.
Now, in the book of Versteeg, formula (11.83) details the face velocity
for RC interpolation as (e.g.)
u^e = (u_P + u_E) / 2 + 1/2 (1/a_P + 1/a_E) (p_P - p_E)
- 1/4 1/a_P (p_W - p_E) - 1/4 1/a_E (p_P - p_EE)
So my understanding, again, is that we actually use this to discretize
the div u^* term on the rhs of the pressure equation. Is that so?

3. When doing so, the "rhs" suddenly contains many new pressure coefficients
in the linear system. I have now two choices: Use the old pressure values
(p^* in the SIMPLE lingo), which is essentially deferred correction, or
use the new values (p' = p as above). When using deferred correction,
I have the problem that the resulting velocity field after the first iteration
won't actually be divergence-free. When coupling p into the system, I
have a wide stencil which I don't really like when using multigrid methods.
So how is this done in practice?

4. And finally, to complete my confusion, I have tried to understand how
OpenFOAM does it from
http://www.tfd.chalmers.se/~hani/kur...7/rhiechow.pdf
They claim that OpenFOAM does not "actually" do RC (which I get since
no wide p stencil is used anywhere). But the paper is quite technical,
and the author never actually points out *where* OpenFOAM differs
from RC. It appears to me that they temporarily use face-based velocities
and then interpolate back, which would make it actually a MAC approach.
But is that true?

Again sorry for the lengthy post. Any insight from the experts would be greatly appreciated!

Cheers

[sbaffini]

This is how the momentum interpolation is done in OpenFOAM:

http://dx.doi.org/10.1016/j.compfluid.2011.11.003

[michujo]

Hi, I'm facing a similar problem when trying to implement Rhie-Chow in my little code here. I also don't like at all the wide stencil, especially when you approach the walls and you have to code variations of that expression depending on the cell type.

As I understand it, the interpolated face velocity will be used to calculate the convective fluxes across the east, west, north and south cell faces (Fe,Fw,Fn and Fs following notation from Versteeg&Malalasekera).

The convective fluxes will be used to compute the discretization coefficients aP, aE, ... based on your differentiation scheme (upwind, QUICK, ...), both for solving the convection-diffusion equation for u* and v*, and the corresponding aP_u and aP_v. I do not think you have to include the whole expression for the calculation of the pressure correction, you just use the old values of pressure in the calculation of ue, uw, un, us. Anyway you will repeat the calculation in the SIMPLE loop so you will eventually converge to the right values. What do you think?

Cheers,
Michujo.

[mercator]

Originally Posted by sbaffini

This is how the momentum interpolation is done in OpenFOAM:

http://dx.doi.org/10.1016/j.compfluid.2011.11.003

Thanks for the link. Unfortunately our uni does not have access to this journal, but I will try to use the library service. Amyway, this could take a few weeks.

Any insights on the other questions?

[mercator]

Originally Posted by michujo

Hi, I'm facing a similar problem when trying to implement Rhie-Chow in my little code here. I also don't like at all the wide stencil, especially when you approach the walls and you have to code variations of that expression depending on the cell type.

As I understand it, the interpolated face velocity will be used to calculate the convective fluxes across the east, west, north and south cell faces (Fe,Fw,Fn and Fs following notation from Versteeg&Malalasekera).

The convective fluxes will be used to compute the discretization coefficients aP, aE, ... based on your differentiation scheme (upwind, QUICK, ...), both for solving the convection-diffusion equation for u* and v*, and the corresponding aP_u and aP_v. I do not think you have to include the whole expression for the calculation of the pressure correction, you just use the old values of pressure in the calculation of ue, uw, un, us. Anyway you will repeat the calculation in the SIMPLE loop so you will eventually converge to the right values. What do you think?

Cheers,
Michujo.

Hi Michujo,

you're right, for ultimately convergence it shouldn't matter whether we use old or new values. Maybe it affects the speed of convergence though..

You could also be right that convective fluxes could be affected, but it seems to me that this is not necessary (?). After all, RC is needed to avoid checkerboard problems when projecting, so I guess one could get away with only doing it during projection. I am not sure though.

[michujo]

Hi. Yes I agree.

However I am not sure how to calculate the convective fluxes for the calculation of the discretization coefficients in the solution of u* and v*. If I do not need to use Rhie-Chow interpolation here, how shall I calculate the face velocities then? Do you think a linear interpolation between adjacent cells will suffice?

So far I have been trying to apply Rhie-Chow in the calculation of the convective fluxes and subsequent solution of u* and v*, but I could be completely wrong. I'll keep reading about it.

Cheers,
Michujo.

[mercator]

Originally Posted by michujo

Hi. Yes I agree.

However I am not sure how to calculate the convective fluxes for the calculation of the discretization coefficients in the solution of u* and v*. If I do not need to use Rhie-Chow interpolation here, how shall I calculate the face velocities then? Do you think a linear interpolation between adjacent cells will suffice?

That's what I do right now. Apart from the checkerboard problem (which I think has to do with the projection) that should work fine.

[RodriguezFatz]

Originally Posted by sbaffini

This is how the momentum interpolation is done in OpenFOAM:

http://dx.doi.org/10.1016/j.compfluid.2011.11.003

Thank you for the link, I just read it. I wish all papers were as good as this one. If you know the authors tell them they are great!

[michujo]

Hi, I finally got my little code to work.

1) Lid-driven cavity flow, 2D, incompressible, Reynolds=20 based on the top wall velocity.

2) Collocated grid, 15x15, uniform.

3) Upwind discretization scheme.

I tried two different methods to interpolate the cell face velocity: linear interpolation between adjacent cells and Rhie&Chow.
I attach two pictures showing the pressure fields. It can be seen that the pressure field obtained using Rhie&Chow is much smoother than the one obtained using linear interpolation.

Coming back to the original question:

1) I applied Rhie&Chow to calculate the cell face velocities and therefore the convective fluxes, which I used to solve the discretized u* and v* equations.

2) I did not explicitly introduce the expression for the cell face velocities in the pressure correction equation, which would have led to a wider stencil. The additional pressure terms appear implicitly in the mass imbalance term through the interpolated cell face velocities. I guess that a explicit substitution of the cell face velocities expressions in the pressure correction equation would probably improve convergence though, but so far this worked fine.

I also noticed that Rhie&Chow requires less than 1/3 of the SIMPLE iterations needed using linear interpolation for convergence. Anyone noticed something similar?

Cheers,
Michujo.

[RodriguezFatz]

What criterion did take for "convergence"? I am asking because the absolute value of residuals can change when you change the equations. Thus, the residual can appear lower without having a better solution.

[mercator]

Originally Posted by michujo

Hi, I finally got my little code to work.

/snip

Cheers,
Michujo.

Very interesting! I think I will try something similar then. Maybe I was confused because the textbook claimed that the continuity eqn was modified, which should not affect convection (pure momentum eqn). Otoh maybe you could do both as you indicated, at the expense of having wider stencils for the pressure eqn.

Just to be sure, so I understand: in the incompressible case the convective flux at the east face is simply F_e = u_e. And that is the velocity you replace by the RC term, using the guess pressure p^* as a corrector. Right?

Cheers

[michujo]

Hi.

1) For convergence I just track the residuals calculated as the difference between the old and new values of u, v and p within the SIMPLE loop. I take the maximum of those three and accept convergence when this is lower than my tolerance (I took 1e-5). Do you think I could implement a better procedure?

2) Hi, yes I agree with you. However I tried to avoid the wider stencil in the pressure correction equation. I'm assuming that the influence of the pressure field is implicitly accounted for in the mass imbalance term calculated with Rhie&Chow.

Yes, I first calculate the cell face velocities with Rhie&Chow (u_east, u_west, v_north and v_south) and with them the convective fluxes as you wrote. Yes, I use p* as the corrector term in the R&C expression.

Maybe it is better to use R&C expression explicitly in the pressure correction equation as well, but so far it seems to be working...

Thanks for your comments.

Cheers,
Michujo.

[RodriguezFatz]

Dear all,
I also busied myself with these implementation details. Unfortunately all explanations

I found (here and in books) weren't satisfying and I hope someone here can shed light on this:

1) In "Ferziger / Peric" Computational Methods for Fluid Dynamics (3rd Edition), the authors explain how modification of the Poisson-equation for the pressure prevents pressure oszillations in colocated arrangements. (page 196, et sqq.). As I understand it, the "outer" derivative (which stems from continity eq.) of the Poisson equation is evaluated using face values of pressure derivatives of the current cell. Thus, only "dx" and "dy" derivatives occour, instead of "2dx", "2dy". Now they argue, that they "introduced an inconsistency in the treatment of the pressure gradient in the momentum and pressure equations" (page 200). I don't get it, can someone tell me, what the inconsinstency is about?

Philipp.

[michujo]

Hi, I think the inconsistency he is talking about is the fourth order pressure term that is added to the momentum equations when interpolating values at the cell faces (where variables are not stored in a collocated arrangement), which is essentially the Rhie&Chow interpolation. Being this extra term fourth order, its effect on the accuracy is negligible for second or third order discretization schemes, yet it damps out pressure oscillations with a 2*deltax pattern (checkerboarding).

Anyone else has a better explanation? Thanks.

Cheers,
Michujo.

[RodriguezFatz]

Thank you for your answer, michujo! I don't see how this can be called "inconsistent". For me, "inconsistent" in this context would mean to treat the same term differently in the two equations. However, I think this isn't done here, right?

[michujo]

Hi, I just read Ferziger again, more carefully this time. As I understood it, the inconsistency consists in using a central difference approximation using points 2*deltax apart in Eq. 7.120 and, on the other hand, using a central difference approximation using points separated deltax in Eq.7.126 (and then interpolating the values of the fluxes at the cell faces). The difference between both approaches is the fourth order term.

I personally find it hard understanding Ferziger's book compared to other textbooks, it's so dense...

Cheers,
Michujo.

[RodriguezFatz]

Yes, that this is what he states. But these are just different approximations of the same equation, of course there is a difference. I am just wondering why he calls that "inconsistent".

Additionally a second thing I am thinking about:
After deriving the "dx"-based pressure equation, he states that
"Use of this eliminates the oscillation in the pressure field".
If this is already sufficient, why all the "Rhie Chow" stuff for the face velocities? Is this (Rhie Chow) made for a second kind of oscillations?