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 Community New Posts Updated Threads Search

Like Tree3Likes

Reply
 
LinkBack Thread Tools Search this Thread Display Modes
Old   June 27, 2018, 09:21
Default reactingTwoPhaseEulerFoam - stability of subcooled boiling simulation
  #1
Member
 
Robin Kamenicky
Join Date: Mar 2016
Posts: 74
Rep Power: 11
Robin.Kamenicky is on a distinguished road
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
Attached Images
File Type: png geometry.png (16.7 KB, 103 views)
Kummi likes this.
Robin.Kamenicky is offline   Reply With Quote

Old   June 28, 2018, 20:27
Default
  #2
New Member
 
Richard
Join Date: Apr 2018
Posts: 24
Rep Power: 8
Qinh is on a distinguished road
Hi Robin,

I am not very familiar with this solver, but I have two questions:

(1) Can the reactingTwoPhaseEulerFoam simulate the subcooled boiling flow? Because I have seen many guys developed their own solvers based on twoPhaseEulerFoam to simulate this kind of problems.

(2) Have you analysize the source code of the reactingTwoPhaseEulerFoam to figure out the models implemented in the solver? Do the models suitable for your simulation?

Best regards,
Qinhao
Qinh is offline   Reply With Quote

Old   June 29, 2018, 11:43
Default
  #3
Member
 
Robin Kamenicky
Join Date: Mar 2016
Posts: 74
Rep Power: 11
Robin.Kamenicky is on a distinguished road
Hi Qinh,

Firstly, thank you for the reply.

Despite it should be done earlier, I have been going through the source code just now. I will understand soon its capabilities much better.

Certainly, I have been reading phd of Henrik Rusche and paper from Chalmers university called "Description of reactingTwoPhaseEulerFoam solver with a focus on mass transfer modelling terms".
I am going to have a look further for boiling with twoPhaseEulerFoam.

Cheers,
Robin
Robin.Kamenicky is offline   Reply With Quote

Old   July 6, 2018, 08:46
Default
  #4
New Member
 
Richard
Join Date: Apr 2018
Posts: 24
Rep Power: 8
Qinh is on a distinguished road
Hello Robin,

Have you simulated subcooled boiling successfully using reactingTwoPhaseEulerFoam or other Foam? I hope you can give some hints on solving the same problems.

Best regards,
Qinhao
Qinh is offline   Reply With Quote

Old   July 8, 2018, 16:43
Default
  #5
Member
 
Robin Kamenicky
Join Date: Mar 2016
Posts: 74
Rep Power: 11
Robin.Kamenicky is on a distinguished road
Hi Qinhao,

well, I have been able to run subcooled boiling with reactingTwoPhaseEulerFoam. However, I have not done any qualitative assessment. Moreover, I still make some changes to adequatly simulate exact case.

Do you want some general hints or discuss your case?

Cheers,
Robin
Robin.Kamenicky is offline   Reply With Quote

Old   July 8, 2018, 21:25
Default
  #6
New Member
 
Richard
Join Date: Apr 2018
Posts: 24
Rep Power: 8
Qinh is on a distinguished road
Hello Robin,

I have tried to use reactingTwoPhaseEulerFoam to simulate subcooled boiling flow under 4.5MPa, just like the experiment conducted by Bartolomej for subcooled boiling heat transfer to water (CFD modelling of subcooled boiling—Concept, validation and application to fuel assembly design). However, the temperature and the pressure increases above the prescribed saturation values and no vapour is formed. I think it may due to the thermal physical properties. For example, in the tutorial wallBoiling/constant/thermophysicalProperties.liquid, the Pr of water is 2.289, but it should be 1.75. I don't know how to set them. Do you have some hints? Thanks a lot!

Best regards,
Qinhao
Qinh is offline   Reply With Quote

Old   July 9, 2018, 05:03
Default
  #7
Member
 
Robin Kamenicky
Join Date: Mar 2016
Posts: 74
Rep Power: 11
Robin.Kamenicky is on a distinguished road
Hi Qinhao,

just a quick question. Have you changed the saturationModel in constant/phaseProperties? The polynom in wallBoiling tutorial does not seem to be valid approximately above 0.5 MPa.

My case is at atmospheric pressure. I will have a further look later.

Best regards,
Robin
Robin.Kamenicky is offline   Reply With Quote

Old   July 9, 2018, 06:57
Default
  #8
New Member
 
Richard
Join Date: Apr 2018
Posts: 24
Rep Power: 8
Qinh is on a distinguished road
Hi Robin,

Yes, the unphysical simulation results may due to the improper saturation model settings. But how to modify it to fit different pressure? I have not find relevent statement about the model.

Beat regards,
Qinhao
Qinh is offline   Reply With Quote

Old   July 9, 2018, 09:35
Default
  #9
Member
 
Robin Kamenicky
Join Date: Mar 2016
Posts: 74
Rep Power: 11
Robin.Kamenicky is on a distinguished road
Hi Qinhao,

if you use as saturationModel the polynomial one then it is defined as

Code:
Polynomial equation for the saturation vapour temperature in terms of the vapour pressure (in Pa).

T_sat = sum_i C_i*p^i

where p is the pressure in Pa and C are the coefficients.
(viz the OF code itself).

I have just calculated a few values and found out that the polynom given in wallBoiling tutorial gives nonsense for the pressure above 0.5 MPa.

You can either make other polynom in the pressure range you need or have a look on other models provided in OF.

Let me know how it goes.

Best regards,
Robin
Robin.Kamenicky is offline   Reply With Quote

Old   July 10, 2018, 02:40
Default
  #10
New Member
 
Richard
Join Date: Apr 2018
Posts: 24
Rep Power: 8
Qinh is on a distinguished road
Hi Robin,

You are right, I have found that the default polynomial is suitable for pressure between 5000 Pa to 0.14 MPa. I fit a new fomula for pressure between 2.3MPa to 7.4MPa, and the simulation result for a 4.5MPa seems resonable. Additionally, I still puzzled about the therma dynamics in the constant/thermalphysicalProperties, what do the Href mean? The Tref may mean the satruation temperature.

The figure shows the fitting formula: https://www.dropbox.com/s/ijvyr1qlbzeuznm/fit.tif?dl=0 and https://www.dropbox.com/s/5m5j5pla6b...omial.tif?dl=0

Best regards,
Qinhao

Last edited by Qinh; July 10, 2018 at 05:35.
Qinh is offline   Reply With Quote

Old   July 10, 2018, 12:42
Default
  #11
Member
 
Robin Kamenicky
Join Date: Mar 2016
Posts: 74
Rep Power: 11
Robin.Kamenicky is on a distinguished road
Hi Qinhao,

for having a look at Hf, Cp, Tref, Href it is necessary to have a look at definition of thermoType dictionary an for options of thermo and equationOfState.

In the case of wallBoiling/constant/thermophysicalProperties.gas
  • thermo is hRefConst
  • equationOfState is perfectGas

In hRefConstThermoI.H can be see:
Hf is equal to Hc. Hc is chemical enthalpy.
Href and Tref are values used to calculate Ha and Hs. Ha is Absolute enthalpy and Hs is Sensible enthalpy.
Ha:
Code:
return Cp_*(T-Tref_) + Href_ + Hf_ + EquationOfState::H(p, T);
Hs:
Code:
return Cp_*(T-Tref_) + Href_ + EquationOfState::H(p, T);
Cp is the one defined in constant/thermophysicalProperties.gas

H(p. T) is zero because we use EquationOfState perfectGas viz. perfectGasI.H

The code mean that Ha=Hs for our case. Also, if current temperature T is equal to Tref, Hs = Hf = Href.

Overall, the values Cp, Tref, Href must correspond same state of specie.

According to ThermalPhaseChangePhaseSystem.C the saturation temperature is taken from saturationModel (not the value Tref. in thermophysicalProperties) as could be expected. Hence the values of Tref and Href in thermophysicalProperties do not have to be defined at the saturation point.

Cheers,
Robin
SHUBHAM9595 and Qinh like this.
Robin.Kamenicky is offline   Reply With Quote

Old   July 11, 2018, 22:02
Default
  #12
New Member
 
Richard
Join Date: Apr 2018
Posts: 24
Rep Power: 8
Qinh is on a distinguished road
Hello Robin,

I have a silly question.

Quote:
Overall, the values Cp, Tref, Href must correspond same state of specie.
I can find the specific heat capacity Cp easily in the widely used look-up table, but how can I define the Tref and Href ?

Best regards,
Qinhao
Qinh is offline   Reply With Quote

Old   July 12, 2018, 06:34
Default
  #13
New Member
 
Richard
Join Date: Apr 2018
Posts: 24
Rep Power: 8
Qinh is on a distinguished road
Hello Robin,

I learned that Tref is the saturation temperature at the system pressure. Href is the enthalpy at Tref and system pressure. They all can be referred in the look-up table.

Best regards,
Qinhao
Qinh is offline   Reply With Quote

Old   July 12, 2018, 09:19
Default
  #14
Member
 
Robin Kamenicky
Join Date: Mar 2016
Posts: 74
Rep Power: 11
Robin.Kamenicky is on a distinguished road
Hi Qinhao,

yes, you are right. You can get those values from any thermodynamic property table.

Only to clarify, as I have already mentioned, the data do not have to be necessarily for saturation point. However, it is sensible to define state point at which you compute (the setup mentioned previously means that cp is constant), hence the saturation point for you.

Cheers,
Robin
Robin.Kamenicky is offline   Reply With Quote

Old   July 13, 2018, 06:50
Default
  #15
New Member
 
Richard
Join Date: Apr 2018
Posts: 24
Rep Power: 8
Qinh is on a distinguished road
Hello Robin,

I have a another question, and I hope you can give me some hints.

The code will easily crash when I simulate a 3D case, so I have to simulate a simplfied 2D case -- transfer the pipe to a wedge type geometry. I tied to create a wedge using ICEM. However, the angle between the faces is suggested to be no more than 5 degrees, which will lead to a mesh of low quality. How to solve the problem? Or how did you create your geometry at the very beginning?

Thank you in advance!

Best regards,
Qinhao
Qinh is offline   Reply With Quote

Old   July 14, 2018, 13:01
Default
  #16
Member
 
Robin Kamenicky
Join Date: Mar 2016
Posts: 74
Rep Power: 11
Robin.Kamenicky is on a distinguished road
Hi Qinhao,

I have also started my computation with wedge geometry. However, I switched to simplified 3D geometry. I have not got a problem with the mesh quality.

In the case of wedge BC, some level of cell skewness is inevitable. Surely, you should keep it as low as possible. Keep the mesh quality high, as usually you would do, especially in the plane your mesh is straddling. To be honest, I do not know how the wedge BC treat the third direction. I have read somewhere that the skewness is not a big issue in this case (for the direction perpendicular to the coordinate plane your mesh is straddling).

Overall, if you are not able to make the mesh better in the third direction, carry on in the computation and just be aware of it. If you find some new information in this topic let me know or somebody else might add some note.

Anyway, have you run checkMesh utility and it has given you some warning?

Cheers,
Robin
Robin.Kamenicky is offline   Reply With Quote

Old   July 14, 2018, 21:15
Default
  #17
New Member
 
Richard
Join Date: Apr 2018
Posts: 24
Rep Power: 8
Qinh is on a distinguished road
Hello Robin,

I create the mesh of the wedge by using ICEM, and the quality is as shown in the attached figure. https://www.dropbox.com/s/tv9g2mnnpv...ality.JPG?dl=0 I think the quality is low. Fortunately, there is no warning when I run checkMesh. So it seems that the mesh can be adopted in the simulation.

Best regards,
Qinhao
Qinh is offline   Reply With Quote

Old   July 15, 2018, 15:16
Default
  #18
Member
 
Robin Kamenicky
Join Date: Mar 2016
Posts: 74
Rep Power: 11
Robin.Kamenicky is on a distinguished road
Hi Qinhao,

I was not able to load the dropBox page.

Anyway, if the checkMesh shows you Mesh OK, than go ahead. Surelly, the checkMesh ok does not automatically mean that the mesh is perfect. According to your non-ortagonality value you might want to setup nNonOrthogonalCorrectors.

Cheers,
Robin
Robin.Kamenicky is offline   Reply With Quote

Old   July 15, 2018, 20:32
Default
  #19
New Member
 
Richard
Join Date: Apr 2018
Posts: 24
Rep Power: 8
Qinh is on a distinguished road
Hello Robin,

I'm sorry about the link, and the website is attached again. https://www.dropbox.com/s/tv9g2mnnpv...ality.JPG?dl=0
Just as you said, I set the nNonOrthogonalCorrectors as 2, and the calculation seems reasonable.

Best regards,
Qinhao
Qinh is offline   Reply With Quote

Old   July 18, 2018, 14:45
Default
  #20
Member
 
Robin Kamenicky
Join Date: Mar 2016
Posts: 74
Rep Power: 11
Robin.Kamenicky is on a distinguished road
Hi Qinhao,

from the little I know about your case, it is difficult to judge the mesh but at the first glance it looks fine.

Cheers,
Robin
Robin.Kamenicky is offline   Reply With Quote

Reply

Tags
boiling, multiphase flow, reactingtwophaseeulerfoam, stability


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 20:54
Two-phase Eulerian model for nucleate subcooled boiling Edy OpenFOAM 5 February 11, 2017 16:21
Wall boiling - subcooled boiling kskong CFX 2 March 9, 2011 20:02
nucleate boiling simulation in CFX Anil CFX 3 August 25, 2010 14:18
Boiling simulation using Fluent Jake Lee FLUENT 0 February 10, 2005 02:09


All times are GMT -4. The time now is 01:08.