[mars91]

Dear Foamers

I'm having a problem in with multiphaseInterFoam (OpenFOAM 4x) that is driving me crazy. I'm trying to simulate the flow around a cylinder in shallow water (2D simulation). The multi-phase simulation includes 3 fluids which are area air, water and sludge. I'm assuming all 3 phases are Newtonian, in-miscible and in compressible. The simulation is in laminar regime (U = 0.1 m/s) the diameter of the cylinder is 20 cm (located at 0,0,0), the bottom of the domain is located at -0.56 m, cause I want to capture the behavior of shallow water as well. The main dimensions of the fluid domain are -3 m the inlet, 3 m outlet, bottom -0.56 m, top 1 m. The overall result looks like the following.

1.jpg

Observations
1: The hydraulic jump in the back of the cylinder shouldn't look like that cause the simulation is in laminar regime.

Screenshot from 2020-10-02 16-10-16.png

2: There is an hydraulic jump right after the inlet BC which is not realistic (see bellow)

Screenshot from 2020-10-02 16-10-42.png

3: The pressure distribution displays extremely high values for a simple laminar case (velocity 0.1 m/s)

4.png

Here I also attached the velocity field

p_rgh5.jpg

I'm not completely sure what I'm doing wrong. I've been trying many different BC and meshing refinements without any success.

If you have any recommendation or any idea of how can I solve this problem I'll be deeply grateful

Here I attached my case setup:

Alpha.mud

 dimensions      [0 0 0 0 0 0 0];

internalField   uniform 0.0;

boundaryField

{
    inletWater
   {
        type            fixedValue;
	value		uniform 0;
    }

    inletAir
    {
        type            fixedValue;
	value		uniform 0;
    }

    inletSand
     {
        type            fixedValue;
	value		uniform 1;
    }
    outletWater
     {
	type            zeroGradient;
    }
    outletAir
     {
	type            zeroGradient;
    }
    outletSand
     {
	type            zeroGradient;
    }
    bottom
    {
        type            zeroGradient;
    }
    atmosphere
   {
        type            zeroGradient;
    }
    cylinder
    {
        type            zeroGradient;
    }

    front
    {
        type            empty;
    }
    back
    {
        type            empty;
    }

}

U

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

internalField   uniform (0.1 0 0);

boundaryField
{
    inletWater
    {
        type            fixedValue;
	value           uniform (0.1 0 0);
    }

    inletAir
    {
        type            fixedValue;
	value           uniform (0.1 0 0);
    }

    inletSand
    {
        type            fixedValue;
	value           uniform (0.1 0 0);
	}
    outletWater
    {
        type            inletOutlet;
        inletValue      uniform (0 0 0);
        value           $internalField;
    }
    outletAir
    {
        type            inletOutlet;
        inletValue      uniform (0 0 0);
        value           $internalField;
    }

    outletSand
   {
        type            inletOutlet;
        inletValue      uniform (0 0 0);
        value           $internalField;
    }

    bottom
    {
        type            movingWallVelocity;
	value           uniform (0.1 0 0);
    }
    atmosphere
    {
        type            fixedValue;
	value           uniform (0.1 0 0);
    }

    cylinder
    {
        type            noSlip;
    }

    front
    {
        type            empty;
    }
    back
    {
        type            empty;
    }

}

p_rgh

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

internalField   uniform 0;

boundaryField

{
    inletWater
   {
        type            zeroGradient;
    }

    inletAir
    {
        type            zeroGradient;
    }

    inletSand
     {
        type            zeroGradient;
    }
    outletWater
     {
        type            zeroGradient;
    }
    outletAir
     {
        type            zeroGradient;
    }
    outletSand
     {
        type            zeroGradient;
    }
    bottom
    {
        type            zeroGradient;
    }
    atmosphere
    {
        type            zeroGradient;   
    }

    cylinder
    {
        type            zeroGradient;
    }

    front
    {
        type            empty;
    }
    back
    {
        type            empty;
    }

}

lambda

dimensions      [0 0 0 0 0 0 0];

internalField   uniform 1.0;

boundaryField
{
    inletWater
    {
        type           zeroGradient;	
    }

    inletAir
     {
        type           zeroGradient;	
    }

    inletSand
     {
        type           zeroGradient;	
    }

    outletWater
    {
        type           zeroGradient;
    }
    outletAir
    {
        type           zeroGradient;
    }

    outletSand
   {
        type           zeroGradient;
    }

    bottom
     {
        type           zeroGradient;
    }
    atmosphere
    {
        type           zeroGradient;
    }

    cylinder
    {
        type           zeroGradient;
    }

    front
    {
        type            empty;
    }
    back
    {
        type            empty;
    }

}

fvSchemes

ddtSchemes
{
    default         Euler;
}

gradSchemes
{
        default         Gauss linear;
	grad(U)		Gauss linear; //Added for the solver
	grad(gamma)	Gauss linear; //Added for the solver
}


divSchemes
{

    div(phi,alpha)	Gauss vanLeer; //Added for solver
    div(phirb,alpha)	Gauss linear; //Added for solver
    div(((rho*nuEff)*dev2(T(grad(U)))))	Gauss linear; //Added for solver
    div(rhoPhi,U)  			Gauss upwind; 

laplacianSchemes
{
    default         Gauss linear corrected;
}

interpolationSchemes
{
    default  linear;
}

snGradSchemes
{
    default         corrected;
}

fvSolution

solvers
{
    "(gamma|alpha).*"
    {
        nAlphacCorr      1;
		nAlphadCorr      1;
		nAlphaSubCycles  1;
        cAlpha           1;


        MULESCorr       no;
        nLimiterIter    3;

        solver          PBiCG;
        preconditioner  DILU;
        tolerance       1e-15;
        relTol          0;
    }

    pcorr
    {
        solver          PCG;
        preconditioner  DIC;
        tolerance       1e-5;
        relTol          0;
    }

	pcorrFinal
	{
	solver		PCG;
	preconditioner	DIC;
	tolerance	1e-5;
	relTol		0;
	}


    p_rgh
    {
        solver          GAMG;
        tolerance       1e-06;
        relTol          0.0001;
        smoother        DIC;
        nPreSweeps      0;
        nPostSweeps     2;
        nFinestSweeps   2;
        cacheAgglomeration true;
        nCellsInCoarsestLevel 10;
        agglomerator    faceAreaPair;
        mergeLevels     1;
    }

    p_rghFinal
    {
        $p_rgh;
        tolerance       1e-6;
        relTol          0;
    }

    "(k|epsilon|omega)"
    {
        solver          PBiCG;
        preconditioner  DILU;
        tolerance       1e-6;
        relTol          0.001;
    }

    "(k|epsilon|omega)Final"
    {
        solver          PBiCG;
        preconditioner  DILU;
        tolerance       1e-6;
        relTol          0;
    }

   lambda
    {
    solver PBiCG;
    preconditioner DILU;
    tolerance 1e-05;
    relTol 0.1;
    }

    lambdaFinal
    {
    $lambda;
    relTol 0;
    }

}


PIMPLE
{
    momentumPredictor no;
    nCorrectors     	1;
	nOuterCorrectors    1;
    nNonOrthogonalCorrectors 0;
    pRefCell        0;
    pRefValue       0;
}

transportProperties

phases 
(
mud

{
    transportModel  Newtonian;
    nu              [0 2 -1 0 0 0 0] 2.36542e-04; //Kinematic viscosity of mud Type E
    rho             [1 -3 0 0 0 0 0] 1226.0;      //Density of mud Type E
}

water
{
    transportModel  Newtonian;
    nu              1.09e-06;
    rho             998.8;
}

air
{
    transportModel  Newtonian;
    nu              [0 2 -1 0 0 0 0] 1.48e-05;
    rho             [1 -3 0 0 0 0 0] 1;
}
);

sigmas //          [1 0 -2 0 0 0 0] 0.07;
(

	(air water) 0.07
	(air mud) 0.07
	(mud water) 0.07
);

Thanks a lot in advance

Regards

Marco

[simrego]

Hi!


I think it won't solve your problem, but.
For the div(rhoPhi,u) don't use upwind. Try with a higher order scheme. linearUpwind or linear. Corrected scheme is enough for you? If your orthogonality is higher than 60-70, try to add a slight limiter to it.
And you are checking the p_rgh as the pressure field while it is the p - rho*g*h field. Check the p field instead!
And are you sure your sigmas are correct?? I have no clue what is the surface tension between the air-mud and mud-water, but It is strange that they are all the same (maybe they are).
And why do you have only zeroGradient p_rgh boundary conditions? At the outlet or atmosphere (if the top is the atmosphere) you should define the pressure not a zeroGradient BC.
Also at the atmosphere why do you have a movingWallVelocity? At the atmosphere you usually don't define the velocity just leave as a zeroGradient or inletOutlet or whatever.


Finally your timeSteps are converged?

[mars91]

Thanks Simrego for your interest in this post.

Yes honestly there is no a reason why I used upwind for the divergence terms, I just didn't change the initial schemes from the base case I used to build mine. I tried with a higher schemes but the result didn't change

I was focusing on the results of the dynamic pressure cause the values are quite extreme for this case I think. Here I attached the p field. The tendency looks ok however I'm wondering why I have such a negative pressure in the air phase?

p.jpg

The sigmas I think are not a big problem here, I run the simulation without surface tension and it didn't change anything.

I'm pretty sure the problem might be in the definition of the BC. With the setting I presented I was trying to model a sort of open channel simulation, but it clearly didn't work. I'll try your recommendations.

The time step I'm using is 1e-4 sec, which for this physics I think is quite small. Not sure that replies to your last question tho.

Looking forward to continue with this discussion

Regards

Marco

[simrego]

Hi!


You have a negative pressure because you have a closed domain based on your boundary conditions, and you define the reference pressure (in closed domains you have to define it) as 0 Pa in the 0th cell (see in your fvSolution/PIMPLE subdict, as pRefValue and pRefCell). And this way you prescribe the pressure as 0 Pa in that cell which is somewhere in the domain (you can find it if you want). But since it is an incompressible case, it does not matter. In incompressible solvers only the pressure gradient is important which will be the same, you can only "shift" the pressure field. As soon as you remove that pRefCell and pRefValue entry probably you won'T be able to run the solver because you have no pressure definition in the boundary conditions or anywhere. That's why you usually set the outlet pressure as 0 in most of the incompressible cases. In that way for example on the inlet you immediately get the pressure drop. But you can set anything, you just shift the pressure field up or down. (Again, this is valid only in incompressible cases!!!)


But the pRefValue and pRefCell stuff is not good in your case because you should not have a closed domain, and after correcting your boundary conditions this negative pressure problem should vanish.