[Pedro24]

Hi ,

For each time step, I found that the order of the computations is a little bit strange.

1) It computes the new surface force term using the alpha field at the time "n"

       twoPhaseProperties.correct();

2) Compute the alpha field "n+1" using the velocity field U at the time "n"

        #include "alphaEqnSubCycle.H"

3) Compute the pressure-velocity field at the time "n+1" using the surface term force at the time "n"

        // --- Pressure-velocity PIMPLE corrector loop
        while (pimple.loop())
        {
            #include "UEqn.H"

            // --- Pressure corrector loop
            while (pimple.correct())
            {
                #include "pEqn.H"
            }

            if (pimple.turbCorr())
            {
                turbulence->correct();
            }
        }

Why don't we inverse the step 1 and 2. It will allows to use the surface term force at the time "n+1" in the navier-stokes equations ?

It seems more logical to update the surface force term (step 1) just after computing the new alpha field. In my mind, the algorithm should be :

1) update properties time "n"
2) compute implicitly velocity-pressure field "n+1"
3) compute alpha field "n+1" using the velocity-pressure previously computed

Am I wrong ?

Thank You,

Pierre

[Pedro24]

I confirm what I said before, this is an error in the "interFoam" solver (and the derivated solvers).

In the H. Rusche thesis, page 162, we can see the following solution procedure :

steps 1-4 : refer to the moving frame

step 5 : transport "alpha1"

step 6 : update properties (smooth gamma => curvature computation => surface force term)

step 7 and 8 : PISO-loop to compute pressure-velocity fields.


Maybe this error should be reported to the OpenFOAM Foundation to correct it for the following versions.


Pierre

[akidess]

How are step 5-8 different than what is already implemented?

[Pedro24]

In the current version of "interFoam", the step 6 is computed before the step 5.

First it computes the surface force term :

twoPhaseProperties.correct();

and then it computes the new alpha1 field :

#include "alphaEqnSubCycle.H"

So, the surface force term used in the PISO loop is computed using the old-time alpha field, which is different from the usual solution procedure.

[akidess]

The interface curvature for the surface force term is calculated by calling interface.correct() after the alphaEqnSubCycle.

[Pedro24]

Ok, I had inverted the "correct" function of twoPhaseMixture for the viscosity with the "correct" function of interfaceProperties for the curvature.

Thank you for your reply