[Stefanie.S.W.]

Hello,

i am running a case with multiphaseEulerFoam (V6), where oil and water mixed with a surfactant collapse together in a cylinder at time 0. (based on the dam break tutorial)
It runs super stable, but when oil and the water solution are seperated again they start oscillating in a wave motion. As if a force from outside is applied. This is not the case, nothing is mooving, the only driving force is gravity.

I got the tip that the therms for pressure and velocity might not converge, so I tried to improve fvScheme and fvSolution. I got help from this great presentation: http://www.wolfdynamics.com/wiki/fvm_crash_intro.pdf
But now unfortunately my simulation crashes within 0.15 sek.

(all simulations run decomposed, 6 cores, simple decomposition)


Do you know which settings in fvSchemes and fvSolutions are wrong?

Can you recommend 2nd Order schemes to me?
Do you have additional tips for me?

(dicts are following)

Thank you a lot for your help!! Unfortunately I can not share the whole geometry (46926 cells, no errors). I hope that it is enough information! Please let me know if I forgot anything.

fvScheme: default settings are in the comments

FoamFile
{
    version     2.0;
    format      ascii;
    class       dictionary;
    location    "system";
    object      fvSchemes;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

ddtSchemes
{
    default   Euler;
}

gradSchemes
{
   default     cellLimited Gauss linear 1;
   grad(U)     cellLimited Gauss linear 1;//default         Gauss linear;
}

divSchemes
{
    "div\(phi,alpha.*\)"    Gauss linear;//vanLeer;
    "div\(phir,alpha.*,alpha.*\)"   Gauss linear;//vanLeer;

    "div\(alphaPhi.*,U.*\)" Gauss limitedLinearV 1;
    div(Rc)                 Gauss linear;
    "div\(phi.*,U.*\)"      Gauss limitedLinearV 1;
    //div(alphaPhi.water,U.water)
}

laplacianSchemes
{
    default         Gauss linear limited 1;//corrected;
}

interpolationSchemes
{
    default         linear;
}

snGradSchemes
{
    default         limited 1;//corrected;
}

//fluxRequired (ALWAYS DEACTIVATED)
//{
//    default         no;
//    p;
//    pcorr;
//}
// ************************************************************************* //

fvSolutions

FoamFile
{
    version     2.0;
    format      ascii;
    class       dictionary;
    location    "system";
    object      fvSolution;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

solvers
{
    "alpha.*"
    {
        nAlphaSubCycles 3; 
    }

    p_rgh
    {
        solver          GAMG;
        tolerance       1e-7;
        relTol          0.01; // 0.05;
        smoother       GaussSeidel;
    }

    p_rghFinal
    {
        solver          PCG;
        preconditioner
        {
            preconditioner  GAMG;
            tolerance       1e-7;
            relTol          0;
            nVcycles        2;
            smoother        GaussSeidel;
        }
        tolerance       1e-7;
        relTol          0;
        maxIter         20;
    }

    "pcorr.*"
    {
        $p_rghFinal;
        tolerance       1e-5;
        relTol          0;
    }

    U
    {
        solver          smoothSolver;
        smoother        GaussSeidel;
        tolerance       1e-8;
        relTol          0;//0.1;
        nSweeps         1;
    }
    
    "U.*" //22.09.2019
    {
        solver          smoothSolver;
        smoother        GaussSeidel;
        tolerance       1e-8;
        relTol          0;//0.1;
        nSweeps         1;
    }

    UFinal
    {
        $U;
        tolerance       1e-7;
        relTol          0;
    }
}

PIMPLE
{
    nCorrectors              3;
    nNonOrthogonalCorrectors 0;
}

relaxationFactors
{
    "U.*"           1;
}

 // ************************************************************************* //

Transport Properties:

FoamFile
{
    version     2.0;
    format      ascii;
    class       dictionary;
    location    "constant";
    object      transportProperties;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

phases
(
     water
     {
         nu         1e-06;
         kappa      1e-06;
         Cp         4195;
         rho        1000;

         diameterModel constant;
         constantCoeffs
         {
             d               5e-4; //13.08.2019 1e-3
         }
     }

     oil 
     {
         nu         9.6e-06;
         kappa      1e-06;
         Cp         4195;
         rho        805;

         diameterModel constant;
         constantCoeffs
         {
             d             1e-4;
         }
     }

     air
     {
         nu         1.48e-05;
         kappa      2.63e-2;
         Cp         1007;
         rho        1;

         diameterModel constant;
         constantCoeffs
         {
             d              2.5e-3;
         }
     }
     lm //Soap or Surfactant
     {
         nu         1.5e-05;
         kappa      1e-06;
         Cp         4195;
         rho        1060;

         diameterModel constant;
         constantCoeffs
         {
             d               4e-5; //2nd try: 5e-5
         }
     }
);

sigmas
(
    (air water)     0.07
    (air oil)       0.03
    (air lm)        0.017

);

interfaceCompression
(
    (air water)     1
    (air oil)       1 
    (air lm)        1
    (water oil)     0
    (water lm)      0
    (lm oil)        0
);

virtualMass
(
    (water oil)     0.2 //Change 19.08.2019
    (oil water)     0.2 //Change 19.08.2019
    (air oil)       0.4 //Change 23.08.2019
    (air water)     0.3 //Change 23.08.2019
    (lm water)      0.5
    (lm oil)        0.5
);

drag
(
    (air water)
    {
        type blended;

        air
        {
            type SchillerNaumann;
            residualPhaseFraction 0;
            residualSlip 0;
        }

        water
        {
            type SchillerNaumann;
            residualPhaseFraction 0;
            residualSlip 0;
        }

        residualPhaseFraction 1e-3;
        residualSlip 1e-3;
    }

    (air oil)
    {
        type blended;

        air
        {
            type SchillerNaumann;
            residualPhaseFraction 0;
            residualSlip 0;
        }

        oil
        {
            type SchillerNaumann;
            residualPhaseFraction 0;
            residualSlip 0;
        }

        residualPhaseFraction 1e-3;
        residualSlip 1e-3;
    }

    (water oil)
    {
        type blended;

        water
        {
            type SchillerNaumann;
            residualPhaseFraction 0;
            residualSlip 0;
        }

        oil
        {
            type SchillerNaumann;
            residualPhaseFraction 0;
            residualSlip 0;
        }

        residualPhaseFraction 1e-3;
        residualSlip 1e-3;
    }
    (air lm)
    {
        type blended;

        air
        {
            type SchillerNaumann;
            residualPhaseFraction 0;
            residualSlip 0;
        }

        lm
        {
            type SchillerNaumann;
            residualPhaseFraction 0;
            residualSlip 0;
        }

        residualPhaseFraction 1e-3;
        residualSlip 1e-3;
    }

    (lm oil)
    {
        type blended;

        lm
        {
            type SchillerNaumann;
            residualPhaseFraction 0;
            residualSlip 0;
        }

        oil
        {
            type SchillerNaumann;
            residualPhaseFraction 0;
            residualSlip 0;
        }

        residualPhaseFraction 1e-3;
        residualSlip 1e-3;
    }
    (lm water)
    {
        type blended;

        lm
        {
            type SchillerNaumann;
            residualPhaseFraction 0;
            residualSlip 0;
        }

        water
        {
            type SchillerNaumann;
            residualPhaseFraction 0;
            residualSlip 0;
        }

        residualPhaseFraction 1e-3;
        residualSlip 1e-3;
    }

 
);

// This is a dummy to support the Smagorinsky model
transportModel  Newtonian;
nu              nu [ 0 2 -1 0 0 0 0 ] 0;
// ************************************************************************* //

Turbulence Properties: LES

\*---------------------------------------------------------------------------*/
FoamFile
{
    version     2.0;
    format      ascii;
    class       dictionary;
    location    "constant";
    object      turbulenceProperties;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

simulationType  LES;

LES  
{
LESModel    Smagorinsky;
turbulence  on;
delta      maxDeltaxyz;
maxDeltaxyzCoeffs
{
    deltaCoeff 2;
}
}

// ************************************************************************* //

Exampe for Velocity:

FoamFile
{
    version     2.0;
    format      ascii;
    class       volVectorField;
    location    "0";
    object      U.water;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

dimensions      [0 1 -1 0 0 0 0];

internalField   uniform (0 0 0);

boundaryField
{
   
    wall
    {
        type               fixedValue;
    value           uniform (0 0 0);
    }
    atmosphere
    {
        type               pressureInletOutletVelocity;
        value              $internalField;
        phi                phi.water;
    }
}


// ************************************************************************* //

p_rgh

FoamFile
{
    version     2.0;
    format      ascii;
    class       volScalarField;
    object      p_rgh;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

dimensions      [1 -1 -2 0 0 0 0];

internalField   uniform 0;

boundaryField
{
    wall
    {
        type            fixedFluxPressure;
        value           uniform 0;
    }

    atmosphere
    {
        type            totalPressure;
        p0              uniform 0;
        U               U.air;
        phi             phi.air;
    }

    defaultFaces
    {
        type            empty;
    }
}

// ************************************************************************* //

nut

FoamFile
{
    version     2.0;
    format      ascii;
    class       volScalarField;
    location    "0";
    object      nut;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

dimensions      [0 2 -1 0 0 0 0];

internalField   uniform 0;

boundaryField
{
    wall
    {
        type            nutkWallFunction;
        value           uniform 0;
    }
    atmosphere
    {
        type            calculated;
        value           uniform 0;
    }
   
    defaultFaces
    {
        type            empty;
    }
}


// ************************************************************************* //

[Stefanie.S.W.]

With the default settings the residuals looked like this:

Courant Number mean: 0.0221092 max: 0.998021
deltaT = 0.000116101
Time = 0.234055

MULES: Solving for alpha.water
water volume fraction, min, max = 0.262989 -5.51866e-06 0.976194
MULES: Solving for alpha.oil
oil volume fraction, min, max = 0.066948 -8.3735e-08 1
MULES: Solving for alpha.air
air volume fraction, min, max = 0.661929 -5.72097e-05 1
MULES: Solving for alpha.lm
lm volume fraction, min, max = 0.00813367 -1.28815e-07 0.0348452
Phase-sum volume fraction, min, max = 1 0.999995 1.00001
MULES: Solving for alpha.water
water volume fraction, min, max = 0.262989 -3.14167e-06 0.976192
MULES: Solving for alpha.oil
oil volume fraction, min, max = 0.066948 -2.08892e-08 1
MULES: Solving for alpha.air
air volume fraction, min, max = 0.661929 -3.78788e-05 1
MULES: Solving for alpha.lm
lm volume fraction, min, max = 0.00813367 -1.2732e-07 0.0348463
Phase-sum volume fraction, min, max = 1 0.999995 1.00001
MULES: Solving for alpha.water
water volume fraction, min, max = 0.262989 -2.65234e-06 0.976189
MULES: Solving for alpha.oil
oil volume fraction, min, max = 0.066948 -2.08932e-08 1
MULES: Solving for alpha.air
air volume fraction, min, max = 0.661929 -1.53512e-05 1
MULES: Solving for alpha.lm
lm volume fraction, min, max = 0.00813367 -1.25778e-07 0.0348474
Phase-sum volume fraction, min, max = 1 0.999995 1.00001
GAMG:  Solving for p_rgh, Initial residual = 0.00172755, Final residual = 1.08339e-05, No Iterations 3
time step continuity errors : sum local = 9.22853e-07, global = 6.84195e-09, cumulative = 7.93007e-05
GAMG:  Solving for p_rgh, Initial residual = 0.00052862, Final residual = 2.28243e-06, No Iterations 4
time step continuity errors : sum local = 1.93883e-07, global = 2.0161e-08, cumulative = 7.93209e-05
GAMGPCG:  Solving for p_rgh, Initial residual = 0.000145525, Final residual = 2.23163e-07, No Iterations 20
time step continuity errors : sum local = 1.89153e-08, global = -8.59243e-10, cumulative = 7.932e-05
ExecutionTime = 1009.12 s

with the changed settings in fvScheme and fvSolution I got the following error:

[4] #0  Foam::error::printStack(Foam::Ostream&) at ??:?
[4] #1  Foam::sigFpe::sigHandler(int) at ??:?
[4] #2  ? in "/lib/x86_64-linux-gnu/libc.so.6"
[4] #3  Foam::GAMGSolver::scale(Foam::Field<double>&, Foam::Field<double>&, Foam::lduMatrix const&, Foam::FieldField<Foam::Field, double> const&, Foam::UPtrList<Foam::lduInterfaceField const> const&, Foam::Field<double> const&, unsigned char) const at ??:?
[4] #4  Foam::GAMGSolver::Vcycle(Foam::PtrList<Foam::lduMatrix::smoother> const&, Foam::Field<double>&, Foam::Field<double> const&, Foam::Field<double>&, Foam::Field<double>&, Foam::Field<double>&, Foam::Field<double>&, Foam::Field<double>&, Foam::PtrList<Foam::Field<double> >&, Foam::PtrList<Foam::Field<double> >&, unsigned char) const at ??:?
[4] #5  Foam::GAMGSolver::solve(Foam::Field<double>&, Foam::Field<double> const&, unsigned char) const at ??:?
[4] #6  Foam::fvMatrix<double>::solveSegregated(Foam::dictionary const&) at ??:?
[4] #7  Foam::fvMatrix<double>::solve(Foam::dictionary const&) in "/opt/openfoam6/platforms/linux64GccDPInt32Opt/bin/multiphaseEulerFoam"
[4] #8  ? in "/opt/openfoam6/platforms/linux64GccDPInt32Opt/bin/multiphaseEulerFoam"
[4] #9  __libc_start_main in "/lib/x86_64-linux-gnu/libc.so.6"
[4] #10  ? in "/opt/openfoam6/platforms/linux64GccDPInt32Opt/bin/multiphaseEulerFoam"

[Stefanie.S.W.]

Hello again,

attached is a picture of the oscillating behavior:

The movement starts "out of nowhere".

Here is my latest try to correct fvSchemes (simulation runs very fast and stable, less than 10 iterations needed per timestep) :

ddtSchemes
 {
     default    CrankNicolson 0;
 }

gradSchemes
 {
     default   cellLimited Gauss linear 1;
     grad(U)   cellLimited Gauss linear 1;
 }

divSchemes
 {
    default                          none;
    "div\(phi.*,U.*\)"                 Gauss linearUpwindV grad(U);
    "div\(phi,alpha.*\)"               Gauss linearUpwind default;
    "div\(phir,alpha.*,alpha.*\)"      Gauss linearUpwind default;
     "div\(alphaPhi.*,U.*\)"           Gauss linearUpwind default;
    div((nuEff*dev(T(grad(U)))))      Gauss linear;
    div(Rc)                 Gauss linear;
 }


laplacianSchemes
 {
   default   Gauss linear limited 1;
  }

interpolationSchemes
 {
  default   linear;
 }

snGradSchemes
 {
    default  limited 1;
 }

Do you have an idea what could cause the oscillating behavior?

[BlnPhoenix]

Could the movement at the interface be due to momentum conservation of the riseing phases?

[BlnPhoenix]

Also for a LES type of simulation your grid looks underresolved, this can be a potential problem.
Try running it in laminar first.

[arjun]

Originally Posted by Stefanie.S.W.

Hello again,

attached is a picture of the oscillating behavior:

The movement starts "out of nowhere".


Do you have an idea what could cause the oscillating behavior?

I have worked with this type of issue in my Software called Wildkatze and one of the reasons for this behaviour is not proper balancing of gravitational force.



It requires very careful balancing of this force. It took me considerable efforts to get this right in Wildkatze (and at some expense of stability too).



A simple way to simulate this behaviour is to have two fluids separated and where heavy fluid is in bottom. (just like your case) and it shall produce 0 movement.

[Stefanie.S.W.]

Thank you a lot for the tips, BlnPhoenix and Arjun!


I am running a laminar test at the moment. Also I am working on refining my mesh.

Maybe I also need different divSchemes for different phases? And "stronger" solver?