Is the Courant number actually 10000, or is this the value you read by OpenFOAM? If so, are you using an unstructured grid? Does the mesh pass checkMesh (skewness?)?

If you use Euler or CN schemes, you're already using an implicit scheme.

Best,

Dear all,

I am very interested in this topic. Actually, I also met the same problem in my unsteady computation. My unsteady computation is unstable if I don't use a relaxation factor.

To my point of view, it is not correct to use relaxation factor in an unsteady simulation. In OpenFoam, the momentum equations are probably (I didn't read the code very deeply) solved by

A . U = H - grad (p) + rhok . g -------------- (Predictor)

I printed the A values at different relaxation factor and they are different. That means the equation you solve becomes

A' . U = H - grad (p) + rhok . g + R(f) -------------- (Predictor)

where R(f) is related to relaxation factor. These A values are also used in pressure correction step. So the solutions will be different no matter how many sub-iterations you use in the pressure correction step.

The equation becomes more explicit with increasing relaxation factor. I think the solution will be the same only when the time step is small enough.

Best regards,
Y. J.

... use coefficients without relaxation!
 

Hi,

thanks for you post and you are right. That is the reason why I suggested to use coefficients without relaxation for the pressure equation in post #26 in this thread you are reading here: http://www.cfd-online.com/Forums/ope...tml#post247123. You have to do that in order to achieve an unsteady, implicit solver stable with arbitrarily large time steps as my thread title implies. However, on reaching convergence, this does not play a role, at least if you do not omit some parts of the coefficients.

So, your "problem" is already solved.

Best regards, Ulf.

Hello,

elegant solutions to this problem have been published (literature search helps ;)), which do not rely on modifications of the iterative procedure (they are not related to it, so why should you change that?), but simply fix the problem at its root, the Rhie-Chow interpolation formula. :D

For example, following Zhang and Zhao (Proceedings of HT-FED04, 2004 ASME Heat Transfer/Fluid Engineering Summer Conference, July 11-15, 2004, Charlotte, North Carolina, USA), you can write (for what possible on a forum ;)) the face velocity as:

u_(i+1/2) = H + B_(i+1/2) DeltaP

where

H = f*H_(i+1) + (1-f)*H_i

and f the interpolation factor. At this point, you can think to apply relaxation as follows:

u_(i+1/2) = URF*(H + B_(i+1/2) * DeltaP) + (1-URF)*U0_(i+1/2)

with U0_(i+1/2) value of U_(i+1/2) at the previous iteration. Now, with some algebra, substituting the definition of H, you find

u_(i+1/2) = u#_(i+1/2) + URF*DeltaX*B_(i+1/2)*DeltaP + (1-URF)*(u0_(i+1/2) - u0#_(i+1/2))

with

u#_(i+1/2) = f*u_(i+1) + (1-f)*u_i

u0#_(i+1/2) = f*u0_(i+1) + (1-f)*u0_i

which ensures the face velocity converges to the correct value, independently from the applied relaxation.

P.S. The dependency on URF is not related to the time step. It is always there, if you use the standard Rhie-Chow formula. Such a formula introduces also a dependency on time step, and the same paper cited above shows how to deal with it (and some treatment should be implemented in ddtPhiCorr, if my understanding is correct).

Best,

Hi,
I enjoyed reading this discussion and I'd like to ask if there is any final conclusion
on the implementation?

I am also interested in this topic.
How to deal with extremely high CFL using pimpleDyMFoam?