[hester]

Hello everyone,

I'm investigating a water air bubble column in a cuboid after Deen.
When using more complex modells for the interface momentum transport coefficients I get problems at the outlet area (see pictures below). With SchillerNaumann for drag coefficient everything worked fine, while with IshiiZuber for drag coefficient the simulation crashes.

This is my case setup:

//========= 1. Mesh =============================//

The cuboid measurements are 0.15 m x 0.15 m x 0.45 m. The inlet at the bottom is modelled as rectangular. The water surface is not included.
I generated my mesh using blockMesh. It is equidistant. CheckMesh is okay.

//========= 2. Boundary conditions ==================//

I am using the following boundary condtitions:

alpha.air

internalField   uniform 0;
boundaryField
{
    inlet
    {
        type            fixedValue;
        value           uniform 1;
    }
    outlet
    {
        type            inletOutlet;
        phi             phi.air;
        inletValue      uniform 0;
        value           uniform 0;
    }
    walls
    {
        type            zeroGradient;
    }
}

U.air

internalField   uniform (0 0 0);
boundaryField
{
    inlet
    {
        type            fixedValue;
        value           uniform (0 0 0.0784);
    }
    outlet
    {
        type            zeroGradient;
    }
    walls
    {
        type            fixedValue;
        value           uniform (0 0 0);
    }
}

U.water

internalField   uniform (0 0 0);
boundaryField
{
    inlet
    {
        type            fixedValue;
        value           uniform (0 0 0);
    }
    outlet
    {
        type            slip;
    }
    walls
    {
        type            fixedValue;
        value           uniform (0 0 0);
    }
}

p

internalField   uniform 100000;
boundaryField
{
    inlet
    {
        type            zeroGradient;
    }
    outlet
    {
        type            fixedValue;
        value           $internalField;
    }
    walls
    {
        type            zeroGradient;
    }
}

k.*

internalField   uniform 1e-4;
boundaryField
{
    inlet
    {
        type             turbulentIntensityKineticEnergyInlet;
        intensity       0.05;
        value            uniform 0;
    }
    outlet
    {
        type            inletOutlet;
        phi             phi.air;
        inletValue      $internalField;
        value           uniform 0;
    }
    walls
    {
        type            zeroGradient;
    }
}

epsilon.air

internalField   uniform 1e-5;
boundaryField
{
    inlet
    {
        type            turbulentMixingLengthDissipationRateInlet;
        mixingLength    1e-3;
        k               k.air;
        value           uniform 0;
    }
    outlet
    {
        type            inletOutlet;
        phi             phi.air;
        inletValue      $internalField;
        value           uniform 0;
    }
    walls
    {
        type            zeroGradient;
    }
}

epsilon.water

internalField   uniform 1e-5;
boundaryField
{
    inlet
    {
        type            fixedValue; 
        value           $internalField;
    }
    outlet
    {
        type            inletOutlet;
        phi             phi.water;
        inletValue      $internalField;
        value           uniform 0;
    }
    walls
    {
        type            zeroGradient;
    }
}

//========= 3. Phase Properties =====================//

I'm assuming a constant air bubble diameter of 4 mm. I only consider drag, lift and virtual mass force. Drag is modelled after Ishii and Zuber as describes here:
http://cfd-online.com/Forums/openfoa...eulerfoam.html
Lift and virtual mass coefficents are both constant with a value of 0.5.

//========= 4. Discretization =======================//

Used methods in fvScheme

ddtSchemes
{
    default         Euler;
}

gradSchemes
{
    default         Gauss linear;
}

divSchemes
{
    
    // from bubble column tutorial
    default                     Gauss linear;

    "div\(phi,alpha.*\)"        Gauss vanLeer;
    "div\(phir,alpha.*\)"       Gauss vanLeer

    "div\(alphaPhi.*,U.*\)"     Gauss limitedLinearV 1;
    "div\(phi.*,U.*\)"          Gauss limitedLinearV 1;
    "div\(alphaPhi.*,.*rho.*\)" Gauss linear;

    "div\(alphaPhi.*,(h|e).*\)" Gauss limitedLinear 1;
    "div\(alphaPhi.*,(K.*|p)\)" Gauss limitedLinear 1;

    "div\(alphaPhi.*,(k|epsilon).*\)"  Gauss limitedLinear 1;
    "div\(phim,(k|epsilon)m\)"    Gauss limitedLinear 1;

    "div\(\(\(alpha.*nuEff.*\)*dev2\(T\(grad\(U.*\)\)\)\)\)" Gauss linear;
}

laplacianSchemes
{
    default         Gauss linear uncorrected;
}

interpolationSchemes
{
    default         linear;
}

snGradSchemes
{
    default         uncorrected;
}

fluxRequired
{
    default         no;
    p               ;
}

//========= 5. Solution ===========================//

I only use GAMG as suggested by my professor.
fvSolution:

solvers
{
    alpha.air
    {
        nAlphaCorr      1;
        nAlphaSubCycles 2;
    }

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

    pFinal
    {
        $p;
        relTol          0;
    }

    
    "e.*"
    {
        solver          GAMG;
        smoother        DILU;
        nPreSweeps      0;
        nPostSweeps     2;
        nFinestSweeps   2;
        cacheAgglomeration true;
        nCellsInCoarsestLevel 160;
        agglomerator    faceAreaPair;
        mergeLevels     2;
        tolerance       1e-6;
        relTol          0.0001;
    }
    
    "(k|epsilon|Theta).*"
    {
        solver          GAMG;
        smoother        DILU;
        nPreSweeps      0;
        nPostSweeps     2;
        nFinestSweeps   2;
        cacheAgglomeration true;
        nCellsInCoarsestLevel 160;
        agglomerator    faceAreaPair;
        mergeLevels     2;
        tolerance       1e-6;
        relTol          0.0001;
    }
}

PIMPLE
{
    nOuterCorrectors 1;
    nCorrectors      2;
    nNonOrthogonalCorrectors 0;
    
    residualControl
    {
        "(U|k|epsilon)"
        {
            relTol          0.0001;
            tolerance       0.0001;
        }
    }

    turbOnFinalIterOnly off;
}

relaxationFactors
{
}

The used Co number is 0.3. The mixtureKEpsilon modell is applied.

//========= 6. Problem ===========================//

After air reaches the outlet area it didn't seem to leave. alpha.air gets big (still <1) in the upper cells and the solution eventually crashes shortly after 7.6 s flow. The flow field for U.water shows there is water inflow. I first set up my case using SchillerNaumann for drag and everything worked fine. But using IshiiZuber instead breaks my simulation. I already tried different things that I found on the forum:

• setting turbulence setting to laminar: same problem, turbulence modell is not the problem
• Co = 0.1 instead of 0.3: same problem, is just takes longer to crash
• nCorrectors = 8 in PISO algorithm: same problem, is just takes much longer to crash
• including the water surface: high velocities for air in area above water surface, simulation crashes
• using 1st order discretization: same problem, athough simulation doesn't crash

I hope someone knows what might be the problem or even has encountered the problem already and knows how to deal with it. I've been struggling with this case for a while now. Any help is greatly appreciated.

Regards,
hester

[hester]

Here are the source files of the used Ishii Zuber drag model.

[demichie]

Hi Hester,
I am wondering why you are using a "slip" boundary condition at the outlet for the liquid velocity and not a "zeroGradient". Have you tested the zeroGradient?

Ciao
Mattia

[hester]

Hello Mattia,

thank you for your answer. I tried zeroGradient for the liquid face at the outlet. I get no simulation problems but the results were wrong. I didn't observe the recirculation pattern in the liquid face. Also was the velocity compared with experimental results too low. So zeroGradient is not what I want here. I use slip because the the outlet boundary is where the water surface would be if I included it in my simulation. I don't want any water to leave the domain.

Regards,
hester

[BlnPhoenix]

Originally Posted by hester

Hello Mattia,

thank you for your answer. I tried zeroGradient for the liquid face at the outlet. I get no simulation problems but the results were wrong. I didn't observe the recirculation pattern in the liquid face. Also was the velocity compared with experimental results too low. So zeroGradient is not what I want here. I use slip because the the outlet boundary is where the water surface would be if I included it in my simulation. I don't want any water to leave the domain.

Regards,
hester

Hello Hester,

were you able to resolve your problem at the outlet region? I think the slip BC is OK to model water surface but i recognized you did not post your alpha.water file. What did you put there for outlet?

For the other approach, when i try to model water and air as two different parts in my domain via setFields i experience that even with a fine mesh the water/air interface (surface) becomes very fuzzy over time which i can not explain or resolve so far..