[dduque]

Hi,

I had this idea to model a Rayleigh-Taylor instability for the simplest possible case. I needed a funny gravity field that would pull up on one of the phases, and down on the other. I hacked the interFoam solver minimally, and I think I have got correct results. My first trial was actually succesful. UEqn.H was changed to:

    MRF.correctBoundaryVelocity(U);

    fvVectorMatrix UEqn
    (
        fvm::ddt(rho, U) + fvm::div(rhoPhi, U)
      + MRF.DDt(rho, U)
      + turbulence->divDevRhoReff(rho, U)
     ==
        fvOptions(rho, U)
            // fake gravity field
            + (2*alpha1-1) * rho * g
    );

    UEqn.relax();


    fvOptions.constrain(UEqn);

    if (pimple.momentumPredictor())
    {
        solve
        (
            UEqn
         ==
            fvc::reconstruct
            (
             (
                    mixture.surfaceTensionForce()
                  - ghf*fvc::snGrad(rho)
//                  - fvc::snGrad(p_rgh)
                  - fvc::snGrad(p)
              ) * mesh.magSf()
            )
         );

        fvOptions.correct(U);
    }

Not much of a change, I introduce the funny gravity and make sure it's p, not p_rgh, the one whose gradient is evaluated.

But then I thought adding the gravity on the right-hand side of the equation would be more correct:

    MRF.correctBoundaryVelocity(U);

    fvVectorMatrix UEqn
    (
        fvm::ddt(rho, U) + fvm::div(rhoPhi, U)
      + MRF.DDt(rho, U)
      + turbulence->divDevRhoReff(rho, U)
     ==
        fvOptions(rho, U)
    );

    UEqn.relax();

    fvOptions.constrain(UEqn);

    if (pimple.momentumPredictor())
    {
        solve
        (
            UEqn
         ==
            fvc::reconstruct
            (
             (
                    mixture.surfaceTensionForce()
                  - ghf*fvc::snGrad(rho)
//                  - fvc::snGrad(p_rgh)
                  - fvc::snGrad(p)
              ) * mesh.magSf()
            )
            // fake gravity field
            + (2*alpha1-1) * rho * g
         );

        fvOptions.correct(U);
    }

Well, no, with the latter version the system does not evolve at all. I wonder why. More importantly, I wonder if I am doing the right thing. Results are ok, but they look somewhat squearish when compared with a Lagrangian (SPH-like) simulation
comparison_step_60.jpg

Thanks in advance for any help!

[Phicau]

Hi Daniel,

probably you have not switched on the momentum predictor in fvSolution, it is usually switched off by default.

Best,

Pablo

[dduque]

Originally Posted by Phicau

Hi Daniel,

probably you have not switched on the momentum predictor in fvSolution, it is usually switched off by default.

Best,

Pablo

Thanks, Pablo.

You are right, it would seem that " if (pimple.momentumPredictor()) " would require the momentum predictor step. However, I have tried enabling it and now the velocity seems to evolve ... but not the pressure, which does not change, or the color field.

Aside from " momentumPredictor yes;" at the PIMPLE section, I had to add the following to fvSolutions.

  UFinal
    {
        $U;
        solver          smoothSolver;
        smoother        symGaussSeidel;
        tolerance       1e-06;
        relTol          0;
    }

What I don't get is, that equation contains the standard rho ( grad p + g ) term, so it should be solved no matter what. Right ??

[It's actually grad( p - rho g z ) + rho z grad(rho), but it's the same ]

[Phicau]

Hi again,

not really, you should always include it in the first one (to assemble the system of equations properly), and if you want a momentum predictor, also in the second equation.

Momentum predictor just solves an intermediate step and updates U. Then it applies that value to obtain the independent terms (explicit) and solve the PPE. Otherwise, the U used is the one from the beginning of the time step.

Best,

Pablo

[dduque]

Originally Posted by Phicau

you should always include it in the first one (to assemble the system of equations properly), and if you want a momentum predictor, also in the second equation.

Thanks, that seems to work. I find little difference with the previous results, but that may be ok. To check incompressibility, I am running foamCalc div phi, and the resulting field is between about -3e-4 and 3e-4. I guess that's small. Perhaps all this could go in another thread.