[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