CFD Online Logo CFD Online URL
www.cfd-online.com
[Sponsors]
Home > Forums > Software User Forums > OpenFOAM > OpenFOAM Running, Solving & CFD

reactingTwoPhaseEulerFoam - stability of subcooled boiling simulation

Register Blogs Members List Search Today's Posts Mark Forums Read

Like Tree3Likes

Reply
 
LinkBack Thread Tools Search this Thread Display Modes
Old   September 20, 2021, 11:15
Default
  #21
New Member
 
Join Date: Mar 2020
Posts: 10
Rep Power: 4
GV871 is on a distinguished road
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





Quote:
Originally Posted by Robin.Kamenicky View Post
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
Code:
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
Code:
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
Code:
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
Code:
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
Code:

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
Code:
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
Code:
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
Code:
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
GV871 is offline   Reply With Quote

Reply

Tags
boiling, multiphase flow, reactingtwophaseeulerfoam, stability

Thread Tools Search this Thread
Search this Thread:

Advanced Search
Display Modes

Posting Rules
You may not post new threads
You may not post replies
You may not post attachments
You may not edit your posts

BB code is On
Smilies are On
[IMG] code is On
HTML code is Off
Trackbacks are Off
Pingbacks are On
Refbacks are On


Similar Threads
Thread Thread Starter Forum Replies Last Post
Fluent simulation subcooled boiling at low pressure cscfx Fluent Multiphase 1 August 9, 2017 21:54
Two-phase Eulerian model for nucleate subcooled boiling Edy OpenFOAM 5 February 11, 2017 17:21
Wall boiling - subcooled boiling kskong CFX 2 March 9, 2011 21:02
nucleate boiling simulation in CFX Anil CFX 3 August 25, 2010 15:18
Boiling simulation using Fluent Jake Lee FLUENT 0 February 10, 2005 03:09


All times are GMT -4. The time now is 17:58.