reactingtwoPhaseEulerFoam: Wall boiling model refinements

Patch contributed by Juho Peltola, VTT.

Resolves patch request https://bugs.openfoam.org/view.php?id=2521

applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/twoPhaseCompressibleTurbulenceModels/derivedFvPatchFields/alphatWallBoilingWallFunction/alphatWallBoilingWallFunctionFvPatchScalarField.C

@@ -86,7 +86,8 @@ alphatWallBoilingWallFunctionFvPatchScalarField
     relax_(0.5),
     AbyV_(p.size(), 0),
     alphatConv_(p.size(), 0),
-    dDep_(p.size(), 0),
+    dDep_(p.size(), 1e-5),
+    qq_(p.size(), 0),
     partitioningModel_(nullptr),
     nucleationSiteModel_(nullptr),
     departureDiamModel_(nullptr),
@@ -114,7 +115,8 @@ alphatWallBoilingWallFunctionFvPatchScalarField
     relax_(dict.lookupOrDefault<scalar>("relax", 0.5)),
     AbyV_(p.size(), 0),
     alphatConv_(p.size(), 0),
-    dDep_(p.size(), 0),
+    dDep_(p.size(), 1e-5),
+    qq_(p.size(), 0),
     partitioningModel_(nullptr),
     nucleationSiteModel_(nullptr),
     departureDiamModel_(nullptr),
@@ -165,6 +167,11 @@ alphatWallBoilingWallFunctionFvPatchScalarField
                 dDep_ = scalarField("dDep", dict, p.size());
             }
 
+            if (dict.found("qQuenching"))
+            {
+                qq_ = scalarField("qQuenching", dict, p.size());
+            }
+
             break;
         }
     }
@@ -203,6 +210,7 @@ alphatWallBoilingWallFunctionFvPatchScalarField
     AbyV_(psf.AbyV_),
     alphatConv_(psf.alphatConv_, mapper),
     dDep_(psf.dDep_, mapper),
+    qq_(psf.qq_, mapper),
     partitioningModel_(psf.partitioningModel_),
     nucleationSiteModel_(psf.nucleationSiteModel_),
     departureDiamModel_(psf.departureDiamModel_),
@@ -221,6 +229,7 @@ alphatWallBoilingWallFunctionFvPatchScalarField
     AbyV_(psf.AbyV_),
     alphatConv_(psf.alphatConv_),
     dDep_(psf.dDep_),
+    qq_(psf.qq_),
     partitioningModel_(psf.partitioningModel_),
     nucleationSiteModel_(psf.nucleationSiteModel_),
     departureDiamModel_(psf.departureDiamModel_),
@@ -240,6 +249,7 @@ alphatWallBoilingWallFunctionFvPatchScalarField
     AbyV_(psf.AbyV_),
     alphatConv_(psf.alphatConv_),
     dDep_(psf.dDep_),
+    qq_(psf.qq_),
     partitioningModel_(psf.partitioningModel_),
     nucleationSiteModel_(psf.nucleationSiteModel_),
     departureDiamModel_(psf.departureDiamModel_),
@@ -436,14 +446,7 @@ void alphatWallBoilingWallFunctionFvPatchScalarField::updateCoeffs()
 
                 const scalarField Tplus(Prt_*(log(E_*yPlus)/kappa_ + P));
                 scalarField Tl(Tw - (Tplus_y250/Tplus)*(Tw - Tc));
-                Tl = max(Tc - 40, min(Tc, Tl));
+                Tl = max(Tc - 40, Tl);
-                const scalarField Tsub(max(Tsatw - Tl, scalar(0)));
-
-                // Wall heat flux partitioning
-                const scalarField fLiquid
-                (
-                    partitioningModel_->fLiquid(liquidw)
-                );
 
                 // Nucleation site density:
                 const scalarField N
@@ -453,7 +456,9 @@ void alphatWallBoilingWallFunctionFvPatchScalarField::updateCoeffs()
                         liquid,
                         vapor,
                         patchi,
-                        Tsatw
+                        Tl,
+                        Tsatw,
+                        L
                     )
                 );
 
@@ -463,7 +468,9 @@ void alphatWallBoilingWallFunctionFvPatchScalarField::updateCoeffs()
                     liquid,
                     vapor,
                     patchi,
-                    Tsub
+                    Tl,
+                    Tsatw,
+                    L
                 );
 
                 // Bubble departure frequency:
@@ -481,19 +488,24 @@ void alphatWallBoilingWallFunctionFvPatchScalarField::updateCoeffs()
                 // Area fractions:
 
                 // Del Valle & Kenning (1985)
-                const scalarField Ja(rhoLiquidw*Cpw*Tsub/(rhoVaporw*L));
+                const scalarField Ja
-                const scalarField Al(fLiquid*4.8*exp(-Ja/80));
+                (
+                    rhoLiquidw*Cpw*(Tsatw - Tl)/(rhoVaporw*L)
+                );
+                const scalarField Al
+                (
+                    fLiquid*4.8*exp( min(-Ja/80,log(VGREAT)))
+                );
 
                 const scalarField A2(min(pi*sqr(dDep_)*N*Al/4, scalar(1)));
                 const scalarField A1(max(1 - A2, scalar(1e-4)));
                 const scalarField A2E(min(pi*sqr(dDep_)*N*Al/4, scalar(5)));
 
-                // Wall evaporation heat flux [kg/s3 = J/m2s]
-                const scalarField Qe((1.0/6.0)*A2E*dDep_*rhoVaporw*fDep*L);
-
                 // Volumetric mass source in the near wall cell due to the
                 // wall boiling
-                dmdt_ = (1 - relax_)*dmdt_ + relax_*Qe*AbyV_/L;
+                dmdt_ =
+                    (1 - relax_)*dmdt_
+                  + relax_*(1.0/6.0)*A2E*dDep_*rhoVaporw*fDep*AbyV_;
 
                 // Volumetric source in the near wall cell due to the wall
                 // boiling
@@ -506,7 +518,7 @@ void alphatWallBoilingWallFunctionFvPatchScalarField::updateCoeffs()
                 );
 
                 // Quenching heat flux
-                const scalarField Qq(A2*hQ*max(Tw - Tl, scalar(0)));
+                qq_ = (A2*hQ*max(Tw - Tl, scalar(0)));
 
                 // Effective thermal diffusivity that corresponds to the
                 // calculated convective, quenching and evaporative heat fluxes
@@ -515,7 +527,7 @@ void alphatWallBoilingWallFunctionFvPatchScalarField::updateCoeffs()
                 (
                     (
                         A1*alphatConv_
-                      + (Qq + Qe)/max(hew.snGrad(), scalar(1e-16))
+                      + (qq_ + qe())/max(hew.snGrad(), scalar(1e-16))
                     )
                    /max(liquidw, scalar(1e-8))
                 );
@@ -538,12 +550,12 @@ void alphatWallBoilingWallFunctionFvPatchScalarField::updateCoeffs()
 
                 if (debug)
                 {
-                    const scalarField Qc
+                    const scalarField qc
                     (
                         fLiquid*A1*(alphatConv_ + alphaw)*hew.snGrad()
                     );
 
-                    const scalarField QEff
+                    const scalarField qEff
                     (
                         liquidw*(*this + alphaw)*hew.snGrad()
                     );
@@ -562,13 +574,13 @@ void alphatWallBoilingWallFunctionFvPatchScalarField::updateCoeffs()
                         << gMax(A2E) << endl;
                     Info<< "  dmdtW: " << gMin(dmdt_) << " - "
                         << gMax(dmdt_) << endl;
-                    Info<< "  Qc: " << gMin(Qc) << " - " << gMax(Qc) << endl;
+                    Info<< "  qc: " << gMin(qc) << " - " << gMax(qc) << endl;
-                    Info<< "  Qq: " << gMin(fLiquid*Qq) << " - "
+                    Info<< "  qq: " << gMin(fLiquid*qq()) << " - "
-                        << gMax(fLiquid*Qq) << endl;
+                        << gMax(fLiquid*qq()) << endl;
-                    Info<< "  Qe: " << gMin(fLiquid*Qe) << " - "
+                    Info<< "  qe: " << gMin(fLiquid*qe()) << " - "
-                        << gMax(fLiquid*Qe) << endl;
+                        << gMax(fLiquid*qe()) << endl;
-                    Info<< "  QEff: " << gMin(QEff) << " - "
+                    Info<< "  qEff: " << gMin(qEff) << " - "
-                        << gMax(QEff) << endl;
+                        << gMax(qEff) << endl;
                     Info<< "  alphat: " << gMin(*this) << " - "
                         << gMax(*this) << endl;
                     Info<< "  alphatConv: " << gMin(alphatConv_)
@@ -636,6 +648,7 @@ void alphatWallBoilingWallFunctionFvPatchScalarField::write(Ostream& os) const
 
     dmdt_.writeEntry("dmdt", os);
     dDep_.writeEntry("dDep", os);
+    qq_.writeEntry("qQuenching", os);
     alphatConv_.writeEntry("alphatConv", os);
     writeEntry("value", os);
 }

applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/twoPhaseCompressibleTurbulenceModels/derivedFvPatchFields/alphatWallBoilingWallFunction/alphatWallBoilingWallFunctionFvPatchScalarField.H

@@ -181,6 +181,9 @@ private:
         //- Departure diameter field
         scalarField dDep_;
 
+        //- Quenching surface heat flux
+        scalarField qq_;
+
         //- Run-time selected heat flux partitioning model
         autoPtr<wallBoilingModels::partitioningModel>
             partitioningModel_;
@@ -275,6 +278,17 @@ public:
             return dDep_;
         }
 
+        //- Return the quenching surface heat flux [W/m2]
+        const scalarField& qq() const
+        {
+            return qq_;
+        }
+
+        //- Return the evaporation surface heat flux [W/m2]
+        tmp<scalarField> qe() const
+        {
+            return mDotL_/AbyV_;
+        }
 
         // Evaluation functions
 

applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/twoPhaseCompressibleTurbulenceModels/derivedFvPatchFields/wallBoilingSubModels/departureDiameterModels/KocamustafaogullariIshii/KocamustafaogullariIshii.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2016 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2016-2017 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -80,7 +80,9 @@ KocamustafaogullariIshii::dDeparture
     const phaseModel& liquid,
     const phaseModel& vapor,
     const label patchi,
-    const scalarField& Tsub
+    const scalarField& Tl,
+    const scalarField& Tsatw,
+    const scalarField& L
 ) const
 {
     // Gravitational acceleration

applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/twoPhaseCompressibleTurbulenceModels/derivedFvPatchFields/wallBoilingSubModels/departureDiameterModels/KocamustafaogullariIshii/KocamustafaogullariIshii.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2016 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2016-2017 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -92,7 +92,9 @@ public:
             const phaseModel& liquid,
             const phaseModel& vapor,
             const label patchi,
-            const scalarField& Tsub
+            const scalarField& Tl,
+            const scalarField& Tsatw,
+            const scalarField& L
         ) const;
 
         virtual void write(Ostream& os) const;

applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/twoPhaseCompressibleTurbulenceModels/derivedFvPatchFields/wallBoilingSubModels/departureDiameterModels/TolubinskiKostanchuk/TolubinskiKostanchuk.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2016 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2016-2017 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -54,7 +54,10 @@ TolubinskiKostanchuk::TolubinskiKostanchuk
     const dictionary& dict
 )
 :
-    departureDiameterModel()
+    departureDiameterModel(),
+    dRef_(dict.lookupOrDefault<scalar>("dRef", 6e-4)),
+    dMax_(dict.lookupOrDefault<scalar>("dMax", 0.0014)),
+    dMin_(dict.lookupOrDefault<scalar>("dMin", 1e-6))
 {}
 
 
@@ -74,10 +77,22 @@ TolubinskiKostanchuk::dDeparture
     const phaseModel& liquid,
     const phaseModel& vapor,
     const label patchi,
-    const scalarField& Tsub
+    const scalarField& Tl,
+    const scalarField& Tsatw,
+    const scalarField& L
 ) const
 {
-    return max(min(0.0006*exp(-Tsub/45), scalar(0.0014)), scalar(1e-6));
+    return max(min(dRef_*exp(-(Tsatw-Tl)/45), dMax_), dMin_);
+}
+
+
+void Foam::wallBoilingModels::departureDiameterModels::
+TolubinskiKostanchuk::write(Ostream& os) const
+{
+    departureDiameterModel::write(os);
+    os.writeKeyword("dRef") << dRef_ << token::END_STATEMENT << nl;
+    os.writeKeyword("dMax") << dMax_ << token::END_STATEMENT << nl;
+    os.writeKeyword("dMin") << dMin_ << token::END_STATEMENT << nl;
 }
 
 

applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/twoPhaseCompressibleTurbulenceModels/derivedFvPatchFields/wallBoilingSubModels/departureDiameterModels/TolubinskiKostanchuk/TolubinskiKostanchuk.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2016 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2016-2017 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -63,6 +63,17 @@ class TolubinskiKostanchuk
     public departureDiameterModel
 {
 
+    // Private data:
+
+        //- Coefficient of the temperature term
+        scalar dRef_;
+
+        //- Maximum diameter
+        scalar dMax_;
+
+        //- Minimum diameter
+        scalar dMin_;
+
 public:
 
     //- Runtime type information
@@ -87,8 +98,12 @@ public:
             const phaseModel& liquid,
             const phaseModel& vapor,
             const label patchi,
-            const scalarField& Tsub
+            const scalarField& Tl,
+            const scalarField& Tsatw,
+            const scalarField& L
         ) const;
+
+        virtual void write(Ostream& os) const;
 };
 
 

applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/twoPhaseCompressibleTurbulenceModels/derivedFvPatchFields/wallBoilingSubModels/departureDiameterModels/departureDiameterModel/departureDiameterModel.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2016 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2016-2017 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -107,7 +107,9 @@ public:
             const phaseModel& liquid,
             const phaseModel& vapor,
             const label patchi,
-            const scalarField& Tsub
+            const scalarField& Tl,
+            const scalarField& Tsatw,
+            const scalarField& L
         ) const = 0;
 
         virtual void write(Ostream& os) const;

applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/twoPhaseCompressibleTurbulenceModels/derivedFvPatchFields/wallBoilingSubModels/nucleationSiteModels/LemmertChawla/LemmertChawla.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2015-2016 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2016-2017 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -53,7 +53,8 @@ Foam::wallBoilingModels::nucleationSiteModels::LemmertChawla::LemmertChawla
     const dictionary& dict
 )
 :
-    nucleationSiteModel()
+    nucleationSiteModel(),
+    Cn_(dict.lookupOrDefault<scalar>("Cn", 1))
 {}
 
 
@@ -71,13 +72,15 @@ Foam::wallBoilingModels::nucleationSiteModels::LemmertChawla::N
     const phaseModel& liquid,
     const phaseModel& vapor,
     const label patchi,
-    const fvPatchScalarField& Tsatw
+    const scalarField& Tl,
+    const scalarField& Tsatw,
+    const scalarField& L
 ) const
 {
     const fvPatchScalarField& Tw =
         liquid.thermo().T().boundaryField()[patchi];
 
-    return 0.8*9.922e5*pow(max((Tw - Tsatw)/10, scalar(0)), 1.805);
+    return Cn_*9.922e5*pow(max((Tw - Tsatw)/10, scalar(0)), 1.805);
 }
 
 

applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/twoPhaseCompressibleTurbulenceModels/derivedFvPatchFields/wallBoilingSubModels/nucleationSiteModels/LemmertChawla/LemmertChawla.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2016 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2016-2017 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -68,6 +68,11 @@ class LemmertChawla
     public nucleationSiteModel
 {
 
+    // Private data:
+
+        //- Coefficient for nucleation site density
+        scalar Cn_;
+
 public:
 
     //- Runtime type information
@@ -91,7 +96,9 @@ public:
             const phaseModel& liquid,
             const phaseModel& vapor,
             const label patchi,
-            const fvPatchScalarField& Tsatw
+            const scalarField& Tl,
+            const scalarField& Tsatw,
+            const scalarField& L
         ) const;
 };
 

applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/twoPhaseCompressibleTurbulenceModels/derivedFvPatchFields/wallBoilingSubModels/nucleationSiteModels/nucleationSiteModel/nucleationSiteModel.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2016 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2016-2017 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -107,7 +107,9 @@ public:
             const phaseModel& liquid,
             const phaseModel& vapor,
             const label patchi,
-            const fvPatchScalarField& Tsatw
+            const scalarField& Tl,
+            const scalarField& Tsatw,
+            const scalarField& L
         ) const = 0;
 
         virtual void write(Ostream& os) const;