[GV871]

I know it's a long shot, but I have similar problem in my simulations.
Can you share how have you solved your problem ?


Thanks

Originally Posted by Robin.Kamenicky

Dear Foamers,

I have a problem with computational stability of boiling case with reactingTwoPhaseEulerFoam. It diverges and floating point error occurs.

Description: I would like to simulate subcooled boiling of water on the attached geometry. Preferably to simulate all heat transfer stages (vapour blanket, nucleate boiling, natural convection).

Current case: Heat transfer is from the patch called solid_wall. The temperature of the wall is elevated highly beyond the saturation temperature of water at given pressure (Currently constant temperature, no all heat transfer stages needs to be currently observed). Fluid domain is a cylinder tank full of water with outlet on the top (patch called fluid_top). The fluid is not agitated.
The geometry is modelled as wedge (axisymmetric). At the outlet, reverse flow of the liquid phase is allowed.
Due to nature of the problem, one can expect abrupt changes and so stability issues.

The case setup is based on the tutorial $FOAM_TUTORIAL/multiphase/reactingTwoPhaseFlow/RAS/wallBoiling. The main differences are that no inlet is present in my case and very high temperature of the heated wall is used.

Computational Attempts:

• run it with PISO algorithm. Computation is more stable but crashes anyway.
• run it with PIMPLE algorithm. The crash speeds up. I am not able to get initial residuals of p_rgh to 1E-4 within 250 nOuterCorrectors (maximum I’ve run), 2 nCorrectors and relaxation factors mentioned in fvSolution. The p_rgh residuals decrease very slowly.
• Initialised k.*, epsilon.*, nut.* with previously run computation with exactly same setting but lower temperature of the heated wall (375K). Problem persists.
• Finer mesh. The simulation crashed even faster.
• K-epsilon turbulence model for liquid and laminar for gas. Problem persists.
• Smaller temperature difference between the heated wall and the fluid. The computation seems to be more stable.
• Many other attempts with changes such as increased minIter for p and h numerical solvers, various relaxation factors, various values of k, epsilon, nut, small initial velocity in the internalField (U.*). Various CFL number, deltaT and maxDeltaT (This helps)

controDict

application     reactingTwoPhaseEulerFoam;

startFrom       latestTime;

startTime       0;

stopAt          endTime;

endTime         40;

deltaT          1e-7;

writeControl    adjustableRunTime;

writeInterval   0.002;

purgeWrite      0;

writeFormat     ascii;

writePrecision  9;

writeCompression uncompressed;

timeFormat      general;

timePrecision   6;

runTimeModifiable yes;

adjustTimeStep  yes;

maxCo           0.01;

maxDeltaT       0.0001;

fvSchemes

ddtSchemes
{
    default         Euler;
}

gradSchemes
{
    default         Gauss linear;
}

divSchemes
{
    default                         none;

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

    "div\(alphaRhoPhi.*,U.*\)"      Gauss limitedLinearV 1;
    "div\(phi.*,U.*\)"              Gauss limitedLinearV 1;

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

    "div\(alphaRhoPhi.*,(k|epsilon).*\)"  Gauss upwind;
    "div\(phim,(k|epsilon)m\)"      Gauss upwind;

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

laplacianSchemes
{
    default         Gauss linear uncorrected;
}

interpolationSchemes
{
    default         linear;
}

snGradSchemes
{
    default         uncorrected;
}

fluxRequired
{
    default         no;
}

wallDist
{
    method          meshWave;
    nRequired       yes;
}

fvSolution

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

    p_rgh
    {
        solver          GAMG;
        smoother        DIC;
        tolerance       1e-8;
        relTol          0.01;
        maxIter         100;
        minIter         2;
    }

    p_rghFinal
    {
        $p_rgh;
        relTol          0;
    }

    "U.*"
    {
        solver          smoothSolver;
        smoother        symGaussSeidel;
        tolerance       1e-6;
        relTol          0;
        minIter         1;
    }

    "(e|h).*"
    {
        solver          smoothSolver;
        smoother        symGaussSeidel;
        tolerance       1e-12;
        relTol          0.001;
        minIter         1;
        maxIter         20;
    }

    "(k|epsilon|Theta).*"
    {
        solver          smoothSolver;
        smoother        symGaussSeidel;
        tolerance       1e-8;
        relTol          0;
        minIter         1;
    }

    Yi
    {
        solver          smoothSolver;
        smoother        symGaussSeidel;
        tolerance       1e-6;
        relTol          0;
        minIter         1;
        residualAlpha   1e-8;
    }
}

PIMPLE
{
    momnetumPredictor   yes;
    nOuterCorrectors    250;
    nCorrectors         2;
    nNonOrthogonalCorrectors 0;
    nEnergyCorrectors   2;
    faceMomentum        yes;
    
    residualControl
        {   
            p_rgh
            {
                relTol 0;
                tolerance 1e-4;
            }
        } 
}

relaxationFactors
{
    fields
    {
        iDmdt           0.1;
        p_rgh           0.1;
        p_rghFinal      0.6;
    }

    equations
    {
        ".*"            1;
        "h.*"           0.1;
        "U.*|k.*|epsilon.*|alpha.*"               0.1;
    }
}

phaseProperties

type    thermalPhaseChangeTwoPhaseSystem;

phases (gas liquid);

volatile    "water";

massTransfer on;

gas
{
    type            multiComponentPhaseModel;
    diameterModel   isothermal;
    constantCoeffs
    {
        d               0.00045;
    }

    isothermalCoeffs
    {
        d0               0.00045;
        p0              1e5;
    }
    Sc              0.7;

    residualAlpha   1e-4;
}

liquid
{
    type            multiComponentPhaseModel;
    diameterModel   constant;
    constantCoeffs
    {
        d               0.00045;
    }
    Sc              0.7;

    residualAlpha   1e-4;
}

blending
{
    default
    {
        type            linear;
        continuousPhase liquid;
        minFullyContinuousAlpha.gas 0.7;
        minPartlyContinuousAlpha.gas 0.5;
        minFullyContinuousAlpha.liquid 0.7;
        minPartlyContinuousAlpha.liquid 0.5;
    }

    heatTransfer
    {
        type            linear;
        continuousPhase liquid;
        minFullyContinuousAlpha.gas 1;
        minPartlyContinuousAlpha.gas 0;
        minFullyContinuousAlpha.liquid 1;
        minPartlyContinuousAlpha.liquid 0;
    }

    massTransfer
    {
        type            linear;
        continuousPhase liquid;
        minFullyContinuousAlpha.gas 1;
        minPartlyContinuousAlpha.gas 0;
        minFullyContinuousAlpha.liquid 1;
        minPartlyContinuousAlpha.liquid 0;
    }
}

surfaceTension
(
    (gas and liquid)
    {
        type            constant;
        sigma           0.07;
    }
);

saturationModel
{
    type polynomial;
    C<8>
    (
        308.0422
        0.0015096
        -1.61589e-8
        1.114106e-13
        -4.52216e-19
        1.05192e-24
        -1.2953e-30
        6.5365e-37
    );
};

aspectRatio
(
    (gas in liquid)
    {
        type            constant;
        E0              1.0;
    }

    (liquid in gas)
    {
        type            constant;
        E0              1.0;
    }
);

drag
(
    (gas in liquid)
    {
        type            SchillerNaumann;
        residualRe      1e-3;
        swarmCorrection
        {
            type        none;
        }
    }

    (liquid in gas)
    {
        type            SchillerNaumann;
        residualRe      1e-3;
        swarmCorrection
        {
            type        none;
        }
    }
);

virtualMass
(
    (gas in liquid)
    {
        type            constantCoefficient;
        Cvm             0.5;
    }

    (liquid in gas)
    {
        type            constantCoefficient;
        Cvm             0.5;
    }
);

interfaceComposition
();

heatTransfer.gas
(
    (gas in liquid)
    {
        type spherical;
        residualAlpha 1e-3;
    }

    (liquid in gas)
    {
        type RanzMarshall;
        residualAlpha 1e-3;
    }
);

heatTransfer.liquid
(
    (gas in liquid)
    {
        type RanzMarshall;
        residualAlpha 1e-3;
    }

    (liquid in gas)
    {
        type spherical;
        residualAlpha 1e-3;
    }
);

massTransfer.gas
();

massTransfer.liquid
();

lift
();

wallLubrication
(
    (gas in liquid)
    {
        type            Antal;
        Cw1             -0.01;
        Cw2             0.05;
        Cwc             10.0;
        Cwd             6.8;
        p               1.7;
    }
);

turbulentDispersion
(
    (gas in liquid)
    {
        type                Burns;
        sigma               0.7;
        Ctd                 1.0;
        residualAlpha       1e-3;
    }
);

// Minimum allowable pressure
pMin            10000;

thermophysicalProperties.gas

thermoType
{
    type            heRhoThermo;
    mixture         multiComponentMixture;
    transport       const;
    thermo          hRefConst;
    equationOfState perfectGas;
    specie          specie;
    energy          sensibleEnthalpy;
}

dpdt yes;

species
(
    water
);

inertSpecie water;

/* chemistryReader foamChemistryReader; */

/* foamChemistryFile "$FOAM_CASE/constant/reactions.gas"; */

water
{
    specie
    {
        molWeight       18.0153;
    }
    equationOfState
    {
        rho        1;
    }

    thermodynamics
    {
        Hf          0;
        Cp          12078.4;
        Tref        373.55;
        Href        2675500;
    }
    transport
    {
        mu          1.2256e-5;
        Pr          2.289;
    }
}

thermophysicalProperties.liquid

thermoType
{
    type            heRhoThermo;
    mixture         multiComponentMixture;
    transport       const;
    thermo          hRefConst;
    equationOfState perfectFluid;
    specie          specie;
    energy          sensibleEnthalpy;
}

dpdt yes;

species
(
    water
);

inertSpecie water;

"(mixture|H2O|water)"
{
    specie
    {
        molWeight       18.0153;
    }
    equationOfState
    {
        R           3000;
        rho0        959;
        rho        959;
    }
    thermodynamics
    {
        Hf          0;
        Cp          4195;
        Tref        373.55;
        Href        417500;
    }
    transport
    {
        mu          2.8291e-4;
        Pr          2.289;
    }
}

p_rgh

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

internalField   uniform 100000;

boundaryField
{
    fluid_bottom
    {
        type            fixedFluxPressure;
    }
    fluid_top
    {
        type            prghTotalPressure;
        U               U.liquid;
        p0               uniform 100000;
        value           uniform 100000;
    }
    tank_wall
    {
        type            fixedFluxPressure;
    }
    solid_wall
    {
        type            fixedFluxPressure;
    }
    right
    {
        type            wedge;
    }
    left
    {
        type            wedge;
    }
}

Part of log

gas min/max T 307.702909 - 2016.44048
liquid min/max T 359.969529 - 2389.91521
Tf.gasAndLiquid: min = 359.983643, mean = 367.376591, max = 403.907099
iDmdt.gasAndLiquid: min = -2137.69475, mean = -92.2414632, max = -2.62873604e-23, integral = -0.000462413491
wDmdt.gasAndLiquid: min = 0, mean = 1.94094026, max = 22.9273172, integral = 2.7677863e-07
dmdt.gasAndLiquid: min = -2114.76743, mean = -90.300523, max = -2.62873604e-23, integral = -0.000462136713
smoothSolver:  Solving for h.gas, Initial residual = 0.000705340909, Final residual = 1.52402318e-13, No Iterations 1
smoothSolver:  Solving for h.liquid, Initial residual = 0.0354795666, Final residual = 9.77502319e-18, No Iterations 1
gas min/max T 307.726027 - 1868.00206
liquid min/max T 359.969529 - 2538.90673
Tf.gasAndLiquid: min = 359.983643, mean = 367.376596, max = 403.907099
iDmdt.gasAndLiquid: min = -2167.11375, mean = -92.7952508, max = -2.63187519e-23, integral = -0.000411679688
wDmdt.gasAndLiquid: min = 0, mean = 1.94114191, max = 22.9373728, integral = 2.76781142e-07
dmdt.gasAndLiquid: min = -2144.17637, mean = -90.8541089, max = -2.63187519e-23, integral = -0.000411402907
GAMG:  Solving for p_rgh, Initial residual = 0.232598235, Final residual = 7.74806493e-17, No Iterations 2
GAMG:  Solving for p_rgh, Initial residual = 0.246949424, Final residual = 1.06193819e-16, No Iterations 2
PIMPLE: iteration 2
MULES: Solving for alpha.gas
MULES: Solving for alpha.gas
MULES: Solving for alpha.gas
alpha.gas volume fraction = 0.000147084939  Min(alpha1) = 5.96787709e-101  Max(alpha1) = 0.247168508
Constructing face momentum equations
smoothSolver:  Solving for h.gas, Initial residual = 0.668877348, Final residual = 7.52677685e-06, No Iterations 2
smoothSolver:  Solving for h.liquid, Initial residual = 0.0363311226, Final residual = 1.12267408e-17, No Iterations 1
gas min/max T -10269.4458 - 2197.07427
liquid min/max T 359.969529 - 2702.79741
[6] #0  Foam::error::printStack(Foam::Ostream&)Tf.gasAndLiquid: min = 359.983642, mean = 367.375326, max = 403.881676
 at ??:?
[6] #1  Foam::sigFpe::sigHandler(int) at ??:?
[6] #2  ? in "/lib/x86_64-linux-gnu/libc.so.6"
[6] #3  Foam::sqrt(Foam::Field<double>&, Foam::UList<double> const&) at ??:?
[6] #4  Foam::tmp<Foam::GeometricField<double, Foam::fvPatchField, Foam::volMesh> > Foam::sqrt<Foam::fvPatchField, Foam::volMesh>(Foam::tmp<Foam::GeometricField<double, Foam::fvPatchField, Foam::volMesh> > const&) at ??:?
[6] #5  Foam::heatTransferModels::RanzMarshall::K(double) const at ??:?
[6] #6  Foam::BlendedInterfacialModel<Foam::heatTransferModel>::K(double) const at ??:?
[6] #7  Foam::ThermalPhaseChangePhaseSystem<Foam::MomentumTransferPhaseSystem<Foam::twoPhaseSystem> >::correctThermo() at ??:?
[6] #8  ? at ??:?
[6] #9  __libc_start_main in "/lib/x86_64-linux-gnu/libc.so.6"
[6] #10  ? at ??:?
[MAE-ROW-RK:25098] *** Process received signal ***
[MAE-ROW-RK:25098] Signal: Floating point exception (8)
[MAE-ROW-RK:25098] Signal code:  (-6)
[MAE-ROW-RK:25098] Failing at address: 0x3e80000620a
[MAE-ROW-RK:25098] [ 0] /lib/x86_64-linux-gnu/libc.so.6(+0x354b0)[0x7f3cc8e714b0]
[MAE-ROW-RK:25098] [ 1] /lib/x86_64-linux-gnu/libc.so.6(gsignal+0x38)[0x7f3cc8e71428]
[MAE-ROW-RK:25098] [ 2] /lib/x86_64-linux-gnu/libc.so.6(+0x354b0)[0x7f3cc8e714b0]
[MAE-ROW-RK:25098] [ 3] /home/robin/OpenFOAM/OpenFOAM-5.0/platforms/linux64GccDPInt32Opt/lib/libOpenFOAM.so(_ZN4Foam4sqrtERNS_5FieldIdEERKNS_5UListIdEE+0x28)[0x7f3cca329f58]
[MAE-ROW-RK:25098] [ 4] /home/robin/OpenFOAM/OpenFOAM-5.0/platforms/linux64GccDPInt32Opt/lib/libreactingTwoPhaseSystem.so(_ZN4Foam4sqrtINS_12fvPatchFieldENS_7volMeshEEENS_3tmpINS_14GeometricFieldIdT_T0_EEEERKS8_+0x10a)[0x7f3cce32d20a]
[MAE-ROW-RK:25098] [ 5] /home/robin/OpenFOAM/OpenFOAM-5.0/platforms/linux64GccDPInt32Opt/lib/libreactingEulerianInterfacialModels.so(_ZNK4Foam18heatTransferModels12RanzMarshall1KEd+0x7d)[0x7f3ccdf79a2d]
[MAE-ROW-RK:25098] [ 6] /home/robin/OpenFOAM/OpenFOAM-5.0/platforms/linux64GccDPInt32Opt/lib/libreactingTwoPhaseSystem.so(_ZNK4Foam23BlendedInterfacialModelINS_17heatTransferModelEE1KEd+0x410)[0x7f3cce393fa0]
[MAE-ROW-RK:25098] [ 7] /home/robin/OpenFOAM/OpenFOAM-5.0/platforms/linux64GccDPInt32Opt/lib/libreactingTwoPhaseSystem.so(_ZN4Foam29ThermalPhaseChangePhaseSystemINS_27MomentumTransferPhaseSystemINS_14twoPhaseSystemEEEE13correctThermoEv+0xa0e)[0x7f3cce394ebe]
[MAE-ROW-RK:25098] [ 8] reactingTwoPhaseEulerFoam[0x436585]
[MAE-ROW-RK:25098] [ 9] /lib/x86_64-linux-gnu/libc.so.6(__libc_start_main+0xf0)[0x7f3cc8e5c830]
[MAE-ROW-RK:25098] [10] reactingTwoPhaseEulerFoam[0x44bb69]
[MAE-ROW-RK:25098] *** End of error message ***
--------------------------------------------------------------------------
mpirun noticed that process rank 6 with PID 25098 on node MAE-ROW-RK exited on signal 8 (Floating point exception).
--------------------------------------------------------------------------

Very high epsilon values are usually observed in the vicinity of the heated wall shortly before the crashes.

I enclose part of the case and a few pictures. Any recommendation/suggestion is welcomed.

Thanks,
Robin