Merge branch 'feature-sub-cooling-ext' into 'develop'

Additional sub-cooling heat transfer correlations for liquid H2

See merge request Development/openfoam!509

applications/solvers/multiphase/compressibleInterFoam/twoPhaseMixtureThermo/twoPhaseMixtureThermo.H

@@ -6,6 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2013-2017 OpenFOAM Foundation
+    Copyright (C) 2021 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -203,6 +204,18 @@ public:
                 const label patchi
             ) const;
 
+            //- Heat capacity using pressure and temperature
+            virtual tmp<scalarField> Cp
+            (
+                const scalarField& p,
+                const scalarField& T,
+                const labelList& cells
+            ) const
+            {
+                NotImplemented;
+                return tmp<scalarField>::New(p);
+            }
+
             //- Heat capacity at constant volume [J/kg/K]
             virtual tmp<volScalarField> Cv() const;
 
@@ -214,6 +227,18 @@ public:
                 const label patchi
             ) const;
 
+            //- Density from pressure and temperature
+            virtual tmp<scalarField> rhoEoS
+            (
+                const scalarField& p,
+                const scalarField& T,
+                const labelList& cells
+            ) const
+            {
+                NotImplemented;
+                return tmp<scalarField>::New(p);
+            }
+
             //- Gamma = Cp/Cv []
             virtual tmp<volScalarField> gamma() const;
 

applications/solvers/multiphase/compressibleMultiphaseInterFoam/multiphaseMixtureThermo/multiphaseMixtureThermo.H

@@ -318,6 +318,18 @@ public:
                 const label patchi
             ) const;
 
+            //- Heat capacity using pressure and temperature
+            virtual tmp<scalarField> Cp
+            (
+                const scalarField& p,
+                const scalarField& T,
+                const labelList& cells
+            ) const
+            {
+                NotImplemented;
+                return tmp<scalarField>::New(p);
+            }
+
             //- Heat capacity at constant volume [J/kg/K]
             virtual tmp<volScalarField> Cv() const;
 
@@ -329,6 +341,18 @@ public:
                 const label patchi
             ) const;
 
+            //- Density from pressure and temperature
+            virtual tmp<scalarField> rhoEoS
+            (
+                const scalarField& p,
+                const scalarField& T,
+                const labelList& cells
+            ) const
+            {
+                NotImplemented;
+                return tmp<scalarField>::New(p);
+            }
+
             //- Gamma = Cp/Cv []
             virtual tmp<volScalarField> gamma() const;
 

applications/solvers/multiphase/interCondensatingEvaporatingFoam/temperaturePhaseChangeTwoPhaseMixtures/twoPhaseMixtureEThermo/twoPhaseMixtureEThermo.H

@@ -5,7 +5,7 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2016-2020 OpenCFD Ltd.
+    Copyright (C) 2016-2021 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -173,6 +173,18 @@ public:
             const label patchi
         ) const;
 
+        //- Heat capacity using pressure and temperature
+        virtual tmp<scalarField> Cp
+        (
+            const scalarField& p,
+            const scalarField& T,
+            const labelList& cells
+        ) const
+        {
+            NotImplemented;
+            return tmp<scalarField>::New(p);
+        }
+
         //- Return Cv of the mixture
         virtual tmp<volScalarField> Cv() const;
 
@@ -184,6 +196,18 @@ public:
             const label patchI
         ) const;
 
+        //- Density from pressure and temperature
+        virtual tmp<scalarField> rhoEoS
+        (
+            const scalarField& p,
+            const scalarField& T,
+            const labelList& cells
+        ) const
+        {
+            NotImplemented;
+            return tmp<scalarField>::New(p);
+        }
+
         //- Gamma = Cp/Cv []
         virtual tmp<volScalarField> gamma() const;
 

src/phaseSystemModels/multiphaseInter/phasesSystem/phaseSystem/phaseSystem.C

@@ -496,6 +496,18 @@ Foam::tmp<Foam::scalarField> Foam::phaseSystem::Cv
 }
 
 
+Foam::tmp<Foam::scalarField> Foam::phaseSystem::rhoEoS
+(
+    const scalarField& p,
+    const scalarField& T,
+    const labelList& cells
+) const
+{
+    NotImplemented;
+    return nullptr;
+}
+
+
 Foam::tmp<Foam::volScalarField> Foam::phaseSystem::gamma() const
 {
     auto iter = phaseModels_.cbegin();

src/phaseSystemModels/multiphaseInter/phasesSystem/phaseSystem/phaseSystem.H

@@ -351,6 +351,18 @@ public:
             const label patchi
         ) const;
 
+        //- Heat capacity using pressure and temperature
+        virtual tmp<scalarField> Cp
+        (
+            const scalarField& p,
+            const scalarField& T,
+            const labelList& cells
+        ) const
+        {
+            NotImplemented;
+            return tmp<scalarField>::New(p);
+        }
+
         //- Return Cv of the mixture
         virtual tmp<volScalarField> Cv() const;
 
@@ -362,6 +374,14 @@ public:
             const label patchI
         ) const;
 
+        //- Density from pressure and temperature
+        virtual tmp<scalarField> rhoEoS
+        (
+            const scalarField& p,
+            const scalarField& T,
+            const labelList& cells
+        ) const;
+
         //- Gamma = Cp/Cv []
         virtual tmp<volScalarField> gamma() const;
 

src/phaseSystemModels/reactingEuler/multiphaseSystem/Make/files

@@ -210,10 +210,12 @@ $(CHFModels)/Zuber/Zuber.C
 CHFSubCoolModels = $(wallBoilingSubModels)/CHFSubCoolModels
 $(CHFSubCoolModels)/CHFSubCoolModel/CHFSubCoolModel.C
 $(CHFSubCoolModels)/HuaXu/HuaXu.C
+$(CHFSubCoolModels)/Tatsumoto/Tatsumoto.C
 
 filmBoiling = $(wallBoilingSubModels)/filmBoilingModels
 $(filmBoiling)/filmBoilingModel/filmBoilingModel.C
 $(filmBoiling)/Bromley/Bromley.C
+$(filmBoiling)/BreenWestwater/BreenWestwater.C
 
 Leidenfrost = $(wallBoilingSubModels)/LeidenfrostModels
 $(Leidenfrost)/LeidenfrostModel/LeidenfrostModel.C
@@ -226,6 +228,12 @@ $(MHFModels)/Jeschar/Jeschar.C
 TDNBModels = $(wallBoilingSubModels)/TDNBModels
 $(TDNBModels)/TDNBModel/TDNBModel.C
 $(TDNBModels)/Schroeder/Schroeder.C
+$(TDNBModels)/Shirai/Shirai.C
+
+nucleateFluxModels = $(wallBoilingSubModels)/nucleateFluxModels
+$(nucleateFluxModels)/nucleateFluxModel/nucleateFluxModel.C
+$(nucleateFluxModels)/Kutadeladze/Kutadeladze.C
+$(nucleateFluxModels)/exponential/exponential.C
 
 /* Turbulence */
 turbulence/multiphaseCompressibleTurbulenceModels.C

src/phaseSystemModels/reactingEuler/multiphaseSystem/derivedFvPatchFields/alphaContactAngle/alphaContactAngleFvPatchScalarField.H

@@ -6,6 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2018 OpenFOAM Foundation
+    Copyright (C) 2021 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -28,7 +29,40 @@ Class
 
 Description
     Contact-angle boundary condition for multi-phase interface-capturing
-    simulations.  Used in conjunction with multiphaseSystem.
+    simulations.  Used in conjunction with \c multiphaseSystem.
+
+Usage
+    Example of the boundary condition specification:
+    \verbatim
+    <patch>
+    {
+        // Mandatory entries
+        type                      alphaContactAngle;
+        thetaProperties
+        (
+            (<phase1> <phase2>) <scalar1> <scalar2> <scalar3> <scalar4>
+            (<phase3> <phase2>) <scalar1> <scalar2> <scalar3> <scalar4>
+            ...
+        );
+
+        // Inherited entries
+        ...
+    }
+    \endverbatim
+
+    where the entries mean:
+    \table
+      Property  | Description                             | Type | Reqd | Deflt
+      type      | Type name: alphaContactAngle            | word | yes  | -
+      thetaProperties | Contact-angle properties          | dict | yes  | -
+      \<scalar1\>  | Equilibrium contact angle            | scalar | yes  |-
+      \<scalar2\>  | Dynamic contact angle velocity scale | scalar | yes  |-
+      \<scalar3\>  | Limiting advancing contact angle     | scalar | yes  |-
+      \<scalar4\>  | Limiting receding contact angle      | scalar | yes  |-
+    \endtable
+
+    The inherited entries are elaborated in:
+     - \link zeroGradientFvPatchFields.H \endlink
 
 SourceFiles
     alphaContactAngleFvPatchScalarField.C
@@ -47,7 +81,7 @@ namespace Foam
 {
 
 /*---------------------------------------------------------------------------*\
-                           Class alphaContactAngleFvPatch Declaration
+                    Class alphaContactAngleFvPatch Declaration
 \*---------------------------------------------------------------------------*/
 
 class alphaContactAngleFvPatchScalarField
@@ -74,13 +108,16 @@ public:
     public:
 
         // Constructors
+
+            //- Default construct
             interfaceThetaProps()
             {}
 
+            //- Construct from Istream
             interfaceThetaProps(Istream&);
 
 
-        // Member functions
+        // Member Functions
 
             //- Return the equilibrium contact angle theta0
             scalar theta0(bool matched=true) const
@@ -90,7 +127,7 @@ public:
             }
 
             //- Return the dynamic contact angle velocity scale
-            scalar uTheta() const
+            scalar uTheta() const noexcept
             {
                 return uTheta_;
             }
@@ -110,7 +147,7 @@ public:
             }
 
 
-        // IO functions
+        // IOstream operators
 
             friend Istream& operator>>(Istream&, interfaceThetaProps&);
             friend Ostream& operator<<(Ostream&, const interfaceThetaProps&);
@@ -126,8 +163,9 @@ public:
 
 private:
 
-    // Private data
+    // Private Data
 
+        //- Interface properties
         thetaPropsTable thetaProps_;
 
 
@@ -155,7 +193,7 @@ public:
         );
 
         //- Construct by mapping given alphaContactAngleFvPatchScalarField
-        //  onto a new patch
+        //- onto a new patch
         alphaContactAngleFvPatchScalarField
         (
             const alphaContactAngleFvPatchScalarField&,
@@ -193,10 +231,10 @@ public:
         }
 
 
-    // Member functions
+    // Member Functions
 
         //- Return the contact angle properties
-        const thetaPropsTable& thetaProps() const
+        const thetaPropsTable& thetaProps() const noexcept
         {
             return thetaProps_;
         }

src/phaseSystemModels/reactingEuler/multiphaseSystem/derivedFvPatchFields/alphatFixedDmdtWallBoilingWallFunction/alphatFixedDmdtWallBoilingWallFunctionFvPatchScalarField.C

@@ -65,7 +65,7 @@ alphatFixedDmdtWallBoilingWallFunctionFvPatchScalarField
 )
 :
     alphatPhaseChangeJayatillekeWallFunctionFvPatchScalarField(p, iF, dict),
-    vaporPhaseName_(dict.lookup("vaporPhase")),
+    vaporPhaseName_(dict.get<word>("vaporPhase")),
     relax_(dict.getOrDefault<scalar>("relax", 1)),
     fixedDmdt_(dict.getOrDefault<scalar>("fixedDmdt", 0)),
     L_(dict.getOrDefault<scalar>("L", 0))
@@ -88,6 +88,8 @@ alphatFixedDmdtWallBoilingWallFunctionFvPatchScalarField
         iF,
         mapper
     ),
+    vaporPhaseName_(psf.vaporPhaseName_),
+    relax_(psf.relax_),
     fixedDmdt_(psf.fixedDmdt_),
     L_(psf.L_)
 {}
@@ -100,6 +102,7 @@ alphatFixedDmdtWallBoilingWallFunctionFvPatchScalarField
 )
 :
     alphatPhaseChangeJayatillekeWallFunctionFvPatchScalarField(psf),
+    vaporPhaseName_(psf.vaporPhaseName_),
     relax_(psf.relax_),
     fixedDmdt_(psf.fixedDmdt_),
     L_(psf.L_)
@@ -114,6 +117,7 @@ alphatFixedDmdtWallBoilingWallFunctionFvPatchScalarField
 )
 :
     alphatPhaseChangeJayatillekeWallFunctionFvPatchScalarField(psf, iF),
+    vaporPhaseName_(psf.vaporPhaseName_),
     relax_(psf.relax_),
     fixedDmdt_(psf.fixedDmdt_),
     L_(psf.L_)
@@ -129,10 +133,8 @@ activePhasePair(const phasePairKey& phasePair) const
     {
         return true;
     }
-    else
-    {
-        return false;
-    }
+
+    return false;
 }
 
 
@@ -143,14 +145,12 @@ dmdt(const phasePairKey& phasePair) const
     {
         return dmdt_;
     }
-    else
-    {
-        FatalErrorInFunction
-            << " dmdt requested for invalid phasePair!"
-            << abort(FatalError);
 
-        return mDotL_;
-    }
+    FatalErrorInFunction
+        << " dmdt requested for invalid phasePair!"
+        << abort(FatalError);
+
+    return mDotL_;
 }
 
 
@@ -161,14 +161,12 @@ mDotL(const phasePairKey& phasePair) const
     {
         return mDotL_;
     }
-    else
-    {
-        FatalErrorInFunction
-            << " mDotL requested for invalid phasePair!"
-            << abort(FatalError);
 
-        return mDotL_;
-    }
+    FatalErrorInFunction
+        << " mDotL requested for invalid phasePair!"
+        << abort(FatalError);
+
+    return mDotL_;
 }
 
 
@@ -179,7 +177,9 @@ void alphatFixedDmdtWallBoilingWallFunctionFvPatchScalarField::updateCoeffs()
         return;
     }
 
-    dmdt_ = (1 - relax_)*dmdt_ + relax_*fixedDmdt_;
+    dmdt_ *= (scalar(1) - relax_);
+    dmdt_ += relax_*fixedDmdt_;
+
     mDotL_ = dmdt_*L_;
 
     operator==(calcAlphat(*this));

src/phaseSystemModels/reactingEuler/multiphaseSystem/derivedFvPatchFields/alphatFixedDmdtWallBoilingWallFunction/alphatFixedDmdtWallBoilingWallFunctionFvPatchScalarField.H

@@ -6,6 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2015-2018 OpenFOAM Foundation
+    Copyright (C) 2021 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -28,9 +29,44 @@ Class
         alphatFixedDmdtWallBoilingWallFunctionFvPatchScalarField
 
 Description
-    A simple alphatPhaseChangeJayatillekeWallFunctionFvPatchScalarField with
+    A simple \c alphatPhaseChangeJayatillekeWallFunctionFvPatchScalarField with
     a fixed volumetric phase-change mass flux.
 
+Usage
+    Example of the boundary condition specification:
+    \verbatim
+    <patch>
+    {
+        // Mandatory entries
+        type         compressible::alphatFixedDmdtWallBoilingWallFunction;
+        vaporPhase   <word>;
+
+        // Optional entries
+        relax        <scalar>;
+        fixedDmdt    <scalar>;
+        L            <scalar>;
+
+        // Inherited entries
+        ...
+    }
+    \endverbatim
+
+    where the entries mean:
+    \table
+      Property  | Description                             | Type | Reqd | Deflt
+      type      | Type name:                             <!--
+      -->  compressible::alphatFixedDmdtWallBoilingWallFunction <!--
+      -->                                                 | word | yes  | -
+      vaporPhase | Name of the vapor phase                | word | yes  | -
+      relax      | Relaxation factor for dmdt           | scalar | no   | 1.0
+      fixedDmdt  | Volumetric phase-change mass flux in near wall cells <!--
+      -->                                               | scalar | no   | 0.0
+      L          | Latent heat                          | scalar | no   | 0.0
+    \endtable
+
+    The inherited entries are elaborated in:
+     -\link alphatPhaseChangeJayatillekeWallFunctionFvPatchScalarField.H\endlink
+
 See also
     Foam::compressible::
         alphatPhaseChangeJayatillekeWallFunctionFvPatchScalarField
@@ -60,12 +96,12 @@ class alphatFixedDmdtWallBoilingWallFunctionFvPatchScalarField
 :
     public alphatPhaseChangeJayatillekeWallFunctionFvPatchScalarField
 {
-    // Private data
+    // Private Data
 
-        //- name on the phase
+        //- Name of the vapor phase
         word vaporPhaseName_;
 
-        //- dmdt relaxationFactor
+        //- Relaxation factor for dmdt
         scalar relax_;
 
         //- Volumetric phase-change mass flux in near wall cells
@@ -99,8 +135,8 @@ public:
         );
 
         //- Construct by mapping given
-        //  alphatFixedDmdtWallBoilingWallFunctionFvPatchScalarField
-        //  onto a new patch
+        //- alphatFixedDmdtWallBoilingWallFunctionFvPatchScalarField
+        //- onto a new patch
         alphatFixedDmdtWallBoilingWallFunctionFvPatchScalarField
         (
             const alphatFixedDmdtWallBoilingWallFunctionFvPatchScalarField&,
@@ -151,7 +187,7 @@ public:
         }
 
 
-    // Member functions
+    // Member Functions
 
         //- Is there phase change mass transfer for this phasePair
         virtual bool activePhasePair(const phasePairKey&) const;
@@ -162,7 +198,8 @@ public:
         //- Return the rate of phase-change for specific phase pair
         virtual const scalarField& mDotL(const phasePairKey&) const;
 
-        // Evaluation functions
+
+        // Evaluation
 
             //- Update the coefficients associated with the patch field
             virtual void updateCoeffs();

src/phaseSystemModels/reactingEuler/multiphaseSystem/derivedFvPatchFields/alphatPhaseChangeJayatillekeWallFunction/alphatPhaseChangeJayatillekeWallFunctionFvPatchScalarField.C

@@ -45,14 +45,12 @@ namespace compressible
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
-scalar alphatPhaseChangeJayatillekeWallFunctionFvPatchScalarField::maxExp_
-    = 50.0;
 scalar alphatPhaseChangeJayatillekeWallFunctionFvPatchScalarField::tolerance_
     = 0.01;
 label alphatPhaseChangeJayatillekeWallFunctionFvPatchScalarField::maxIters_
     = 10;
 
-// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+// * * * * * * * * * * * * * Protected Member Functions  * * * * * * * * * * * //
 
 void alphatPhaseChangeJayatillekeWallFunctionFvPatchScalarField::checkType()
 {
@@ -83,18 +81,18 @@ alphatPhaseChangeJayatillekeWallFunctionFvPatchScalarField::yPlusTherm
     const scalarField& Prat
 ) const
 {
-    tmp<scalarField> typsf(new scalarField(this->size()));
-    scalarField& ypsf = typsf.ref();
+    auto typsf = tmp<scalarField>::New(this->size());
+    auto& ypsf = typsf.ref();
 
     forAll(ypsf, facei)
     {
         scalar ypt = 11.0;
 
-        for (int i=0; i<maxIters_; i++)
+        for (int i = 0; i < maxIters_; ++i)
         {
-            scalar f = ypt - (log(E_*ypt)/kappa_ + P[facei])/Prat[facei];
-            scalar df = 1 - 1.0/(ypt*kappa_*Prat[facei]);
-            scalar yptNew = ypt - f/df;
+            const scalar f = ypt - (log(E_*ypt)/kappa_ + P[facei])/Prat[facei];
+            const scalar df = 1.0 - 1.0/(ypt*kappa_*Prat[facei]);
+            const scalar yptNew = ypt - f/df;
 
             if (yptNew < VSMALL)
             {
@@ -116,26 +114,23 @@ alphatPhaseChangeJayatillekeWallFunctionFvPatchScalarField::yPlusTherm
     return typsf;
 }
 
+
 tmp<scalarField>
 alphatPhaseChangeJayatillekeWallFunctionFvPatchScalarField::calcAlphat
 (
     const scalarField& prevAlphat
 ) const
 {
-
     // Lookup the fluid model
     const phaseSystem& fluid =
         db().lookupObject<phaseSystem>("phaseProperties");
 
-    const phaseModel& phase
-    (
-        fluid.phases()[internalField().group()]
-    );
+    const phaseModel& phase = fluid.phases()[internalField().group()];
 
     const label patchi = patch().index();
 
     // Retrieve turbulence properties from model
-    const phaseCompressibleTurbulenceModel& turbModel =
+    const auto& turbModel =
         db().lookupObject<phaseCompressibleTurbulenceModel>
         (
             IOobject::groupName(turbulenceModel::propertiesName, phase.name())
@@ -163,11 +158,6 @@ alphatPhaseChangeJayatillekeWallFunctionFvPatchScalarField::calcAlphat
     const fvPatchScalarField& hew =
         phase.thermo().he().boundaryField()[patchi];
 
-    const fvPatchScalarField& Tw =
-        phase.thermo().T().boundaryField()[patchi];
-
-    scalarField Tp(Tw.patchInternalField());
-
     // Heat flux [W/m2] - lagging alphatw
     const scalarField qDot
     (
@@ -178,18 +168,19 @@ alphatPhaseChangeJayatillekeWallFunctionFvPatchScalarField::calcAlphat
 
     scalarField yPlus(uTau*y/(muw/rhow));
 
-    scalarField Pr(muw/alphaw);
+    const scalarField Pr(muw/alphaw);
 
     // Molecular-to-turbulent Prandtl number ratio
-    scalarField Prat(Pr/Prt_);
+    const scalarField Prat(Pr/Prt_);
 
     // Thermal sublayer thickness
-    scalarField P(this->Psmooth(Prat));
+    const scalarField P(this->Psmooth(Prat));
 
-    scalarField yPlusTherm(this->yPlusTherm(P, Prat));
+    tmp<scalarField> tyPlusTherm = this->yPlusTherm(P, Prat);
+    const scalarField& yPlusTherm = tyPlusTherm.cref();
 
-    tmp<scalarField> talphatConv(new scalarField(this->size()));
-    scalarField& alphatConv = talphatConv.ref();
+    auto talphatConv = tmp<scalarField>::New(this->size());
+    auto& alphatConv = talphatConv.ref();
 
     // Populate boundary values
     forAll(alphatConv, facei)
@@ -198,19 +189,20 @@ alphatPhaseChangeJayatillekeWallFunctionFvPatchScalarField::calcAlphat
         scalar alphaEff = 0.0;
         if (yPlus[facei] < yPlusTherm[facei])
         {
-            scalar A = qDot[facei]*rhow[facei]*uTau[facei]*y[facei];
-            scalar B = qDot[facei]*Pr[facei]*yPlus[facei];
-            scalar C = Pr[facei]*0.5*rhow[facei]*uTau[facei]*sqr(magUp[facei]);
+            const scalar A = qDot[facei]*rhow[facei]*uTau[facei]*y[facei];
+            const scalar B = qDot[facei]*Pr[facei]*yPlus[facei];
+            const scalar C =
+                Pr[facei]*0.5*rhow[facei]*uTau[facei]*sqr(magUp[facei]);
             alphaEff = A/(B + C + VSMALL);
         }
         else
         {
-            scalar A = qDot[facei]*rhow[facei]*uTau[facei]*y[facei];
-            scalar B =
+            const scalar A = qDot[facei]*rhow[facei]*uTau[facei]*y[facei];
+            const scalar B =
                 qDot[facei]*Prt_*(1.0/kappa_*log(E_*yPlus[facei]) + P[facei]);
-            scalar magUc =
+            const scalar magUc =
                 uTau[facei]/kappa_*log(E_*yPlusTherm[facei]) - mag(Uw[facei]);
-            scalar C =
+            const scalar C =
                 0.5*rhow[facei]*uTau[facei]
                *(Prt_*sqr(magUp[facei]) + (Pr[facei] - Prt_)*sqr(magUc));
             alphaEff = A/(B + C + VSMALL);

src/phaseSystemModels/reactingEuler/multiphaseSystem/derivedFvPatchFields/alphatPhaseChangeJayatillekeWallFunction/alphatPhaseChangeJayatillekeWallFunctionFvPatchScalarField.H

@@ -6,6 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2015-2018 OpenFOAM Foundation
+    Copyright (C) 2021 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -33,26 +34,39 @@ Description
     the Eulerian multiphase solvers.
 
 Usage
-    \table
-        Property     | Description             | Required    | Default value
-        Prt          | Turbulent Prandtl number | no         | 0.85
-        Cmu          | Model coefficient       | no          | 0.09
-        kappa        | von Karman constant     | no          | 0.41
-        E            | Model coefficient       | no          | 9.8
-    \endtable
-
     Example of the boundary condition specification:
     \verbatim
     <patchName>
     {
+        // Mandatory entries
         type            alphatPhaseChangeJayatillekeWallFunction;
-        Prt             0.85;
-        kappa           0.41;
-        E               9.8;
-        value           uniform 0; // optional value entry
+
+        // Optional entries
+        Prt             <scalar>;
+        Cmu             <scalar>;
+        kappa           <scalar>;
+        E               <scalar>;
+
+        // Inherited entries
+        ...
     }
     \endverbatim
 
+    where the entries mean:
+    \table
+      Property  | Description                             | Type | Reqd | Deflt
+      type      | Type name:                             <!--
+      -->  compressible::alphatPhaseChangeJayatillekeWallFunction <!--
+      -->                                                 | word | yes  | -
+      Prt       | Turbulent Prandtl number              | scalar | no   | 0.85
+      Cmu       | Empirical model coefficient           | scalar | no   | 0.09
+      kappa     | Von Karman constant                   | scalar | no   | 0.41
+      E         | Wall roughness parameter              | scalar | no   | 9.8
+    \endtable
+
+    The inherited entries are elaborated in:
+      - \link alphatPhaseChangeWallFunctionFvPatchScalarField.H \endlink
+
 See also
     Foam::compressible::alphatPhaseChangeWallFunctionFvPatchScalarField
 
@@ -81,27 +95,28 @@ class alphatPhaseChangeJayatillekeWallFunctionFvPatchScalarField
 :
     public alphatPhaseChangeWallFunctionFvPatchScalarField
 {
-
 protected:
 
-    // Protected data
+    // Protected Data
 
         //- Turbulent Prandtl number
         scalar Prt_;
 
-        //- Cmu coefficient
+        //- Empirical model coefficient
         scalar Cmu_;
 
         //- Von Karman constant
         scalar kappa_;
 
-        //- E coefficient
+        //- Wall roughness parameter
         scalar E_;
 
         // Solution parameters
 
-            static scalar maxExp_;
+            //- Absolute tolerance
             static scalar tolerance_;
+
+            //- Maximum number of iterations
             static label maxIters_;
 
 
@@ -151,8 +166,8 @@ public:
         );
 
         //- Construct by mapping given
-        //  alphatPhaseChangeJayatillekeWallFunctionFvPatchScalarField
-        //  onto a new patch
+        //- alphatPhaseChangeJayatillekeWallFunctionFvPatchScalarField
+        //- onto a new patch
         alphatPhaseChangeJayatillekeWallFunctionFvPatchScalarField
         (
             const alphatPhaseChangeJayatillekeWallFunctionFvPatchScalarField&,
@@ -203,9 +218,9 @@ public:
         }
 
 
-    // Member functions
+    // Member Functions
 
-        // Evaluation functions
+        // Evaluation
 
             //- Update the coefficients associated with the patch field
             virtual void updateCoeffs();

src/phaseSystemModels/reactingEuler/multiphaseSystem/derivedFvPatchFields/alphatPhaseChangeWallFunction/alphatPhaseChangeWallFunctionFvPatchScalarField.H

@@ -6,6 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2015-2018 OpenFOAM Foundation
+    Copyright (C) 2021 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -29,9 +30,33 @@ Class
 Description
     Abstract base-class for all alphatWallFunctions supporting phase-change.
 
+Usage
+    Example of the boundary condition specification:
+    \verbatim
+    <patchName>
+    {
+        // Optional entries
+        dmdt        <scalarField>;
+        mDotL       <scalarField>;
+
+        // Inherited entries
+        ...
+    }
+    \endverbatim
+
+    where the entries mean:
+    \table
+      Property  | Description                            | Type | Reqd | Deflt
+      dmdt      | Rate of phase-change            | scalarField | no   | 0.0
+      mDotL     | Latent heat of the phase-change | scalarField | no   | 0.0
+    \endtable
+
+    The inherited entries are elaborated in:
+      - \link fixedValueFvPatchFields.H \endlink
+
 See also
-    Foam::fixedValueFvPatchScalarField
-    Foam::alphatWallFunctionFvPatchScalarField
+  - Foam::fixedValueFvPatchScalarField
+  - Foam::alphatWallFunctionFvPatchScalarField
 
 SourceFiles
     alphatPhaseChangeWallFunctionFvPatchScalarField.C
@@ -94,8 +119,8 @@ public:
         );
 
         //- Construct by mapping given
-        //  alphatPhaseChangeWallFunctionFvPatchScalarField
-        //  onto a new patch
+        //- alphatPhaseChangeWallFunctionFvPatchScalarField
+        //- onto a new patch
         alphatPhaseChangeWallFunctionFvPatchScalarField
         (
             const alphatPhaseChangeWallFunctionFvPatchScalarField&,
@@ -165,9 +190,9 @@ public:
         }
 
 
-        // Evaluation Functions
+        // Evaluation
 
-        //- Update the coefficients associated with the patch field
+            //- Update the coefficients associated with the patch field
             virtual void updateCoeffs() = 0;
 
 

src/phaseSystemModels/reactingEuler/multiphaseSystem/derivedFvPatchFields/alphatWallBoilingWallFunction/alphatWallBoilingWallFunctionFvPatchScalarField.C

@@ -85,13 +85,16 @@ alphatWallBoilingWallFunctionFvPatchScalarField
     nucleationSiteModel_(nullptr),
     departureDiamModel_(nullptr),
     departureFreqModel_(nullptr),
+    nucleatingModel_(nullptr),
     filmBoilingModel_(nullptr),
     LeidenfrostModel_(nullptr),
     CHFModel_(nullptr),
     CHFSoobModel_(nullptr),
     MHFModel_(nullptr),
     TDNBModel_(nullptr),
-    wp_(1)
+    wp_(1),
+    liquidTatYplus_(false),
+    regimeTypes_(p.size(), -1)
 {
     AbyV_ = this->patch().magSf();
     forAll(AbyV_, facei)
@@ -113,7 +116,7 @@ alphatWallBoilingWallFunctionFvPatchScalarField
     alphatPhaseChangeJayatillekeWallFunctionFvPatchScalarField(p, iF, dict),
     otherPhaseName_(dict.get<word>("otherPhase")),
     phaseType_(phaseTypeNames_.get("phaseType", dict)),
-    relax_(Function1<scalar>::New("relax", dict, &db())),
+    relax_(Function1<scalar>::New("relax", dict)),
     AbyV_(p.size(), 0),
     alphatConv_(p.size(), 0),
     dDep_(p.size(), 1e-5),
@@ -123,13 +126,16 @@ alphatWallBoilingWallFunctionFvPatchScalarField
     nucleationSiteModel_(nullptr),
     departureDiamModel_(nullptr),
     departureFreqModel_(nullptr),
+    nucleatingModel_(nullptr),
     filmBoilingModel_(nullptr),
     LeidenfrostModel_(nullptr),
     CHFModel_(nullptr),
     CHFSoobModel_(nullptr),
     MHFModel_(nullptr),
     TDNBModel_(nullptr),
-    wp_(1)
+    wp_(1),
+    liquidTatYplus_(dict.getOrDefault<bool>("liquidTatYplus", false)),
+    regimeTypes_(p.size(), -1)
 {
     // Check that otherPhaseName != this phase
     if (internalField().group() == otherPhaseName_)
@@ -158,23 +164,43 @@ alphatWallBoilingWallFunctionFvPatchScalarField
         }
         case liquidPhase:
         {
-            nucleationSiteModel_ =
-                wallBoilingModels::nucleationSiteModel::New
+            partitioningModel_ =
+                wallBoilingModels::partitioningModel::New
                 (
-                    dict.subDict("nucleationSiteModel")
+                    dict.subDict("partitioningModel")
                 );
 
-            departureDiamModel_ =
-                wallBoilingModels::departureDiameterModel::New
-                (
-                    dict.subDict("departureDiamModel")
-                );
+            // If nucleating model is specifed use it. Otherwise use
+            // RPI wall boiling model
+            const dictionary* nucleatingDict
+                = dict.findDict("nucleateFluxModel");
+
+            if (nucleatingDict)
+            {
+                nucleatingModel_ =
+                wallBoilingModels::nucleateFluxModel::New(*nucleatingDict);
+            }
+            else
+            {
+                nucleationSiteModel_ =
+                    wallBoilingModels::nucleationSiteModel::New
+                    (
+                        dict.subDict("nucleationSiteModel")
+                    );
+
+                departureDiamModel_ =
+                    wallBoilingModels::departureDiameterModel::New
+                    (
+                        dict.subDict("departureDiamModel")
+                    );
+
+                departureFreqModel_ =
+                    wallBoilingModels::departureFrequencyModel::New
+                    (
+                        dict.subDict("departureFreqModel")
+                    );
+            }
 
-            departureFreqModel_ =
-                wallBoilingModels::departureFrequencyModel::New
-                (
-                    dict.subDict("departureFreqModel")
-                );
 
             const dictionary* LeidenfrostDict =
                 dict.findDict("LeidenfrostModel");
@@ -284,13 +310,16 @@ alphatWallBoilingWallFunctionFvPatchScalarField
     partitioningModel_(psf.partitioningModel_),
     nucleationSiteModel_(psf.nucleationSiteModel_),
     departureDiamModel_(psf.departureDiamModel_),
+    nucleatingModel_(psf.nucleatingModel_),
     filmBoilingModel_(psf.filmBoilingModel_),
     LeidenfrostModel_(psf.LeidenfrostModel_),
     CHFModel_(psf.CHFModel_),
     CHFSoobModel_(psf.CHFSoobModel_),
     MHFModel_(psf.MHFModel_),
     TDNBModel_(psf.TDNBModel_),
-    wp_(psf.wp_)
+    wp_(psf.wp_),
+    liquidTatYplus_(psf.liquidTatYplus_),
+    regimeTypes_(psf.regimeTypes_)
 {}
 
 
@@ -312,13 +341,16 @@ alphatWallBoilingWallFunctionFvPatchScalarField
     partitioningModel_(psf.partitioningModel_),
     nucleationSiteModel_(psf.nucleationSiteModel_),
     departureDiamModel_(psf.departureDiamModel_),
+    nucleatingModel_(psf.nucleatingModel_),
     filmBoilingModel_(psf.filmBoilingModel_),
     LeidenfrostModel_(psf.LeidenfrostModel_),
     CHFModel_(psf.CHFModel_),
     CHFSoobModel_(psf.CHFSoobModel_),
     MHFModel_(psf.MHFModel_),
     TDNBModel_(psf.TDNBModel_),
-    wp_(psf.wp_)
+    wp_(psf.wp_),
+    liquidTatYplus_(psf.liquidTatYplus_),
+    regimeTypes_(psf.regimeTypes_)
 {}
 
 
@@ -341,13 +373,16 @@ alphatWallBoilingWallFunctionFvPatchScalarField
     partitioningModel_(psf.partitioningModel_),
     nucleationSiteModel_(psf.nucleationSiteModel_),
     departureDiamModel_(psf.departureDiamModel_),
+    nucleatingModel_(psf.nucleatingModel_),
     filmBoilingModel_(psf.filmBoilingModel_),
     LeidenfrostModel_(psf.LeidenfrostModel_),
     CHFModel_(psf.CHFModel_),
     CHFSoobModel_(psf.CHFSoobModel_),
     MHFModel_(psf.MHFModel_),
     TDNBModel_(psf.TDNBModel_),
-    wp_(psf.wp_)
+    wp_(psf.wp_),
+    liquidTatYplus_(psf.liquidTatYplus_),
+    regimeTypes_(psf.regimeTypes_)
 {}
 
 
@@ -360,10 +395,8 @@ activePhasePair(const phasePairKey& phasePair) const
     {
         return true;
     }
-    else
-    {
-        return false;
-    }
+
+    return false;
 }
 
 const scalarField& alphatWallBoilingWallFunctionFvPatchScalarField::
@@ -373,14 +406,12 @@ dmdt(const phasePairKey& phasePair) const
     {
         return dmdt_;
     }
-    else
-    {
-        FatalErrorInFunction
-            << " dmdt requested for invalid phasePair!"
-            << abort(FatalError);
 
-        return dmdt_;
-    }
+    FatalErrorInFunction
+        << " dmdt requested for invalid phasePair!"
+        << abort(FatalError);
+
+    return dmdt_;
 }
 
 const scalarField& alphatWallBoilingWallFunctionFvPatchScalarField::
@@ -390,19 +421,17 @@ mDotL(const phasePairKey& phasePair) const
     {
         return mDotL_;
     }
-    else
-    {
-        FatalErrorInFunction
-            << " mDotL requested for invalid phasePair!"
-            << abort(FatalError);
 
-        return mDotL_;
-    }
+    FatalErrorInFunction
+        << " mDotL requested for invalid phasePair!"
+        << abort(FatalError);
+
+    return mDotL_;
 }
 
+
 void alphatWallBoilingWallFunctionFvPatchScalarField::updateCoeffs()
 {
-
     if (updated())
     {
         return;
@@ -423,13 +452,13 @@ void alphatWallBoilingWallFunctionFvPatchScalarField::updateCoeffs()
             db().lookupObject<phaseSystem>("phaseProperties")
         );
 
-    const saturationModel& satModel =
+    const auto& satModel =
         db().lookupObject<saturationModel>("saturationModel");
 
     const label patchi = patch().index();
 
     const scalar t = this->db().time().timeOutputValue();
-    scalar relax = relax_->value(t);
+    const scalar relax = relax_->value(t);
 
     switch (phaseType_)
     {
@@ -440,11 +469,18 @@ void alphatWallBoilingWallFunctionFvPatchScalarField::updateCoeffs()
                 fluid.phases()[internalField().group()]
             );
 
+            const tmp<scalarField> talphaw = vapor.thermo().alpha(patchi);
+            const scalarField& alphaw = talphaw();
+
             const fvPatchScalarField& hewv =
                 vapor.thermo().he().boundaryField()[patchi];
 
             // Vapor Liquid phase fraction at the wall
-            const scalarField vaporw(vapor.boundaryField()[patchi]);
+            const scalarField vaporw
+            (
+                max(vapor.boundaryField()[patchi], scalar(1e-16))
+            );
+            const scalarField liquidw(1.0 - vaporw);
 
             // NOTE! Assumes 1-thisPhase for liquid fraction in
             // multiphase simulations
@@ -453,11 +489,6 @@ void alphatWallBoilingWallFunctionFvPatchScalarField::updateCoeffs()
                 partitioningModel_->fLiquid(1-vaporw)
             );
 
-            const tmp<scalarField> talphaw = vapor.thermo().alpha(patchi);
-            const scalarField& alphaw = talphaw();
-
-            const scalarField heSnGrad(max(hewv.snGrad(), scalar(1e-16)));
-
             // Convective thermal diffusivity for single phase
             const scalarField alphatv(calcAlphat(*this));
 
@@ -465,7 +496,7 @@ void alphatWallBoilingWallFunctionFvPatchScalarField::updateCoeffs()
             {
                 this->operator[](i) =
                 (
-                    (1 - fLiquid[i])*(alphatv[i])
+                    (1 - fLiquid[i])*(alphatv[i] + alphaw[i])
                    /max(vaporw[i], scalar(1e-8))
                 );
             }
@@ -479,6 +510,9 @@ void alphatWallBoilingWallFunctionFvPatchScalarField::updateCoeffs()
 
                 const scalarField qEff(vaporw*(*this + alphaw)*hewv.snGrad());
 
+                 Info<< "  qEffVap: " << gMin(qEff) << " - "
+                     << gMax(qEff) << endl;
+
                 scalar Qeff = gSum(qEff*patch().magSf());
                 Info<< " Effective heat transfer rate to vapor:" << Qeff
                     << nl << endl;
@@ -488,27 +522,30 @@ void alphatWallBoilingWallFunctionFvPatchScalarField::updateCoeffs()
         case liquidPhase:
         {
             // Check that nucleationSiteModel has been constructed
-            if (!nucleationSiteModel_)
+            if (!nucleatingModel_)
             {
-                FatalErrorInFunction
-                    << "nucleationSiteModel has not been constructed!"
-                    << abort(FatalError);
-            }
+                if (!nucleationSiteModel_)
+                {
+                    FatalErrorInFunction
+                        << "nucleationSiteModel has not been constructed!"
+                        << abort(FatalError);
+                }
 
-            // Check that departureDiameterModel has been constructed
-            if (!departureDiamModel_)
-            {
-                FatalErrorInFunction
-                    << "departureDiameterModel has not been constructed!"
-                    << abort(FatalError);
-            }
+                // Check that departureDiameterModel has been constructed
+                if (!departureDiamModel_)
+                {
+                    FatalErrorInFunction
+                        << "departureDiameterModel has not been constructed!"
+                        << abort(FatalError);
+                }
 
-            // Check that nucleationSiteModel has been constructed
-            if (!departureFreqModel_)
-            {
-                FatalErrorInFunction
-                    << "departureFrequencyModel has not been constructed!"
-                    << abort(FatalError);
+                // Check that nucleationSiteModel has been constructed
+                if (!departureFreqModel_)
+                {
+                    FatalErrorInFunction
+                        << "departureFrequencyModel has not been constructed!"
+                        << abort(FatalError);
+                }
             }
 
             const phaseModel& liquid
@@ -519,7 +556,7 @@ void alphatWallBoilingWallFunctionFvPatchScalarField::updateCoeffs()
             const phaseModel& vapor(fluid.phases()[otherPhaseName_]);
 
             // Retrieve turbulence properties from models
-            const phaseCompressibleTurbulenceModel& turbModel =
+            const auto& turbModel =
                 db().lookupObject<phaseCompressibleTurbulenceModel>
                 (
                     IOobject::groupName
@@ -528,7 +565,7 @@ void alphatWallBoilingWallFunctionFvPatchScalarField::updateCoeffs()
                         liquid.name()
                     )
                 );
-            const phaseCompressibleTurbulenceModel& vaporTurbModel =
+            const auto& vaporTurbModel =
                 db().lookupObject<phaseCompressibleTurbulenceModel>
                 (
                     IOobject::groupName
@@ -540,7 +577,7 @@ void alphatWallBoilingWallFunctionFvPatchScalarField::updateCoeffs()
 
             const tmp<scalarField> tnutw = turbModel.nut(patchi);
 
-            const scalar Cmu25(pow025(Cmu_));
+            const scalar Cmu25 = pow025(Cmu_);
 
             const scalarField& y = turbModel.y()[patchi];
 
@@ -593,7 +630,7 @@ void alphatWallBoilingWallFunctionFvPatchScalarField::updateCoeffs()
                 satModel.Tsat(liquid.thermo().p());
 
             const volScalarField& Tsat = tTsat();
-            const fvPatchScalarField& Tsatw(Tsat.boundaryField()[patchi]);
+            const fvPatchScalarField& Tsatw = Tsat.boundaryField()[patchi];
             const scalarField Tsatc(Tsatw.patchInternalField());
 
             const fvPatchScalarField& pw =
@@ -602,34 +639,40 @@ void alphatWallBoilingWallFunctionFvPatchScalarField::updateCoeffs()
             const fvPatchScalarField& hew =
                 liquid.thermo().he().boundaryField()[patchi];
 
-            const scalarField hw
+            const scalarField hwLiqSat
             (
                 liquid.thermo().he().member() == "e"
-              ? hew.patchInternalField() + pw/rhow.patchInternalField()
-              : hew.patchInternalField()
+              ? liquid.thermo().he(pw, Tsatc, patchi)
+                    + pw/rhow.patchInternalField()
+              : liquid.thermo().he(pw, Tsatc, patchi)
             );
 
             const scalarField L
             (
                 vapor.thermo().he().member() == "e"
-              ? vapor.thermo().he(pw, Tsatc, patchi) + pw/rhoVaporw - hw
-              : vapor.thermo().he(pw, Tsatc, patchi) - hw
+              ? vapor.thermo().he(pw, Tsatc, patchi) + pw/rhoVaporw - hwLiqSat
+              : vapor.thermo().he(pw, Tsatc, patchi) - hwLiqSat
             );
 
             // Liquid phase fraction at the wall
             const scalarField liquidw(liquid.boundaryField()[patchi]);
 
+            // Partition between phases
             const scalarField fLiquid(partitioningModel_->fLiquid(liquidw));
 
+            scalarField Tl(Tc);
             // Liquid temperature at y+=250 is estimated from logarithmic
             // thermal wall function (Koncar, Krepper & Egorov, 2005)
-            const scalarField Tplus_y250
-            (
-                Prt_*(log(E_*250)/kappa_ + P)
-            );
-            const scalarField Tplus(Prt_*(log(E_*yPlus)/kappa_ + P));
-            scalarField Tl(Tw - (Tplus_y250/Tplus)*(Tw - Tc));
-            Tl = max(Tc - 40, Tl);
+            if (liquidTatYplus_)
+            {
+                const scalarField Tplus_y250
+                (
+                    Prt_*(log(E_*250)/kappa_ + P)
+                );
+                const scalarField Tplus(Prt_*(log(E_*yPlus)/kappa_ + P));
+                scalarField Tl(Tw - (Tplus_y250/Tplus)*(Tw - Tc));
+                Tl = max(Tc - 40, Tl);
+            }
 
             // Film, transient boiling regimes
             scalarField Qtb(this->size(), 0);
@@ -661,6 +704,11 @@ void alphatWallBoilingWallFunctionFvPatchScalarField::updateCoeffs()
                     )
                 );
 
+                if (debug == 2)
+                {
+                    Info << "CHF : " << CHF << endl;
+                }
+
                 // Effect of sub-cooling to the CHF in saturated conditions
                 const scalarField CHFSubCool
                 (
@@ -675,6 +723,11 @@ void alphatWallBoilingWallFunctionFvPatchScalarField::updateCoeffs()
                     )
                 );
 
+                if (debug == 2)
+                {
+                    Info << "CHF Sub Cool factor : " << CHFSubCool << endl;
+                }
+
                 const scalarField CHFtotal(CHF*CHFSubCool);
 
                 tDNB =
@@ -688,6 +741,12 @@ void alphatWallBoilingWallFunctionFvPatchScalarField::updateCoeffs()
                         L
                     );
 
+                if (debug == 2)
+                {
+                    Info<< "Temperature departure from biling : "
+                        << tDNB << endl;
+                }
+
                 const scalarField MHF
                 (
                     MHFModel_->MHF
@@ -701,6 +760,11 @@ void alphatWallBoilingWallFunctionFvPatchScalarField::updateCoeffs()
                     )
                 );
 
+                if (debug == 2)
+                {
+                    Info<< "MHF : " << MHF << endl;
+                }
+
                 TLeiden =
                     LeidenfrostModel_->TLeid
                     (
@@ -712,6 +776,11 @@ void alphatWallBoilingWallFunctionFvPatchScalarField::updateCoeffs()
                         L
                     );
 
+                if (debug == 2)
+                {
+                    Info<< "Leidenfrost Temp : " << TLeiden << endl;
+                }
+
                 // htc for film boiling
                 htcFilmBoiling =
                     filmBoilingModel_->htcFilmBoil
@@ -724,8 +793,12 @@ void alphatWallBoilingWallFunctionFvPatchScalarField::updateCoeffs()
                         L
                     );
 
-                // htc for film transition boiling
+                if (debug == 2)
+                {
+                    Info<< "Htc film boiling : " << htcFilmBoiling << endl;
+                }
 
+                // htc for film transition boiling
                 // Indicator between CHF (phi = 0) and MHF (phi = 1)
                 const scalarField phi
                 (
@@ -744,11 +817,45 @@ void alphatWallBoilingWallFunctionFvPatchScalarField::updateCoeffs()
 
             }
 
+            // Convective heat transfer area for Sub-cool boiling
+            scalarField A1(this->size(), 0);
+            qq_ = Zero;
+            scalarField dmdtSubCooling(this->size(), 0);
 
-            // Sub-cool boiling Nucleation
-            const scalarField N
-            (
-                nucleationSiteModel_->N
+            if (nucleatingModel_)
+            {
+                dmdtSubCooling =
+                    nucleatingModel_->qNucleate
+                    (
+                        liquid,
+                        vapor,
+                        patchi,
+                        Tl,
+                        Tsatw,
+                        L
+                    )*AbyV_/L;
+
+
+                dmdtSubCooling *= fLiquid;
+            }
+            else
+            {
+                // Sub-cool boiling Nucleation
+                const scalarField N
+                (
+                    nucleationSiteModel_->N
+                    (
+                        liquid,
+                        vapor,
+                        patchi,
+                        Tl,
+                        Tsatw,
+                        L
+                    )
+                );
+
+                // Bubble departure diameter:
+                dDep_ = departureDiamModel_->dDeparture
                 (
                     liquid,
                     vapor,
@@ -756,44 +863,75 @@ void alphatWallBoilingWallFunctionFvPatchScalarField::updateCoeffs()
                     Tl,
                     Tsatw,
                     L
-                )
-            );
+                );
 
-            // Bubble departure diameter:
-            dDep_ = departureDiamModel_->dDeparture
-            (
-                liquid,
-                vapor,
-                patchi,
-                Tl,
-                Tsatw,
-                L
-            );
-
-            // Bubble departure frequency:
-            const scalarField fDep
-            (
-                departureFreqModel_->fDeparture
+                // Bubble departure frequency:
+                const scalarField fDep
                 (
-                    liquid,
-                    vapor,
-                    patchi,
-                    dDep_
-                )
-            );
+                    departureFreqModel_->fDeparture
+                    (
+                        liquid,
+                        vapor,
+                        patchi,
+                        dDep_
+                    )
+                );
+
+                scalarField Ja
+                (
+                    rhow*Cpw*(Tsatw - Tl)/(rhoVaporw*L)
+                );
+
+                scalarField Al
+                (
+                    fLiquid*4.8*exp(min(-Ja/80, log(VGREAT)))
+                );
+
+                scalarField A2
+                (
+                    min(pi*sqr(dDep_)*N*Al/4, scalar(1))
+                );
+
+                A1 = max(1 - A2, scalar(1e-4));
+
+                // Following Bowring(1962)
+                scalarField A2E
+                (
+                    min(pi*sqr(dDep_)*N*Al/4, scalar(5))
+                );
+
+                dmdtSubCooling =
+                (
+                    (1.0/6.0)*A2E*dDep_*rhoVaporw*fDep*AbyV_
+                );
+
+                scalarField hQ
+                (
+                    2*(alphaw*Cpw)*fDep
+                    *sqrt
+                    (
+                        (0.8/max(fDep, SMALL))/(pi*alphaw/rhow)
+                    )
+                );
+
+                // Quenching heat flux in Sub-cool boiling
+                qq_ =
+                    (
+                        (1 - relax)*qq_
+                      + relax*A2*hQ*max(Tw - Tl, scalar(0))
+                    );
+            }
 
             // Convective thermal diffusivity for single phase
             alphatConv_ = calcAlphat(alphatConv_);
 
-            // Convective heat transfer area for Sub-cool boiling
-            scalarField A1(this->size(), 0);
-
             const scalarField hewSn(hew.snGrad());
 
+            // AlphaEff for film regime
             scalarField alphaFilm(this->size(), 0);
 
             // Use to identify regimes per face
-            labelField regimeTypes(A1.size(), -1);
+            regimeTypes_ = -1;
 
             forAll(*this, i)
             {
@@ -803,63 +941,17 @@ void alphatWallBoilingWallFunctionFvPatchScalarField::updateCoeffs()
                     if (Tw[i] < tDNB[i])
                     {
                         // Sub-cool boiling
-                        regimeTypes[i] = regimeType::subcool;
+                        regimeTypes_[i] = regimeType::subcool;
 
-                        // Del Valle & Kenning (1985)
-                        const scalar Ja =
-                            rhow[i]*Cpw[i]*(Tsatw[i] - Tl[i])
-                            /(rhoVaporw[i]*L[i]);
-
-                        const scalar Al =
-                            fLiquid[i]*4.8*exp(min(-Ja/80, log(VGREAT)));
-
-                        const scalar A2
-                        (
-                            min(pi*sqr(dDep_[i])*N[i]*Al/4, scalar(1))
-                        );
-
-                        A1[i] = max(1 - A2, scalar(1e-4));
-
-                        // Following Bowring(1962)
-                        const scalar A2E
-                        (
-                            min
-                            (
-                                pi*sqr(dDep_[i])*N[i]*Al/4,
-                                scalar(5)
-                            )
-                        );
-
-                        // Volumetric mass source in the near wall cell due
-                        // to the wall boiling
                         dmdt_[i] =
-                            (
-                                (1 - relax)*dmdt_[i]
-                                + relax*(1.0/6.0)*A2E*dDep_[i]*rhoVaporw[i]
-                                * fDep[i]*AbyV_[i]
-                            );
+                        (
+                            (1 - relax)*dmdt_[i] + relax*dmdtSubCooling[i]
+                        );
 
                         // Volumetric source in the near wall cell due to
                         // the wall boiling
                         mDotL_[i] = dmdt_[i]*L[i];
 
-                        // Quenching heat transfer coefficient
-                        const scalar hQ
-                        (
-                            2*(alphaw[i]*Cpw[i])*fDep[i]
-                            *sqrt
-                            (
-                                (0.8/max(fDep[i], SMALL))/(pi*alphaw[i]/rhow[i])
-                            )
-                        );
-
-                        // Quenching heat flux in Sub-cool boiling
-                        qq_[i] =
-                            (
-                                (1 - relax)*qq_[i]
-                                + relax*A2*hQ*max(Tw[i] - Tl[i], scalar(0))
-                            );
-
                         this->operator[](i) =
                         (
                             max
@@ -867,17 +959,27 @@ void alphatWallBoilingWallFunctionFvPatchScalarField::updateCoeffs()
                                 A1[i]*alphatConv_[i]
                                 + (
                                     (qq_[i] + mDotL_[i]/AbyV_[i])
-                                    / max(hewSn[i], scalar(1e-16))
+                                  / max(hewSn[i], scalar(1e-16))
                                 )
                                 /max(liquidw[i], scalar(1e-8)),
-                                1e-8
+                                scalar(1e-8)
                             )
                         );
+
+                        if (debug == 2)
+                        {
+                            Info<< "Sub-cool boiling:  " <<  nl
+                                << "  fraction Liq: " << fLiquid[i] << nl
+                                << "  Heat flux: "
+                                << (qq_[i] + mDotL_[i]/AbyV_[i]) << nl
+                                << "  delta Tsub: " << (Tw[i] - Tsatw[i])
+                                << endl;
+                        }
                     }
                     else if (Tw[i] > tDNB[i] && Tw[i] < TLeiden[i])
                     {
                         // transient boiling
-                        regimeTypes[i] = regimeType::transient;
+                        regimeTypes_[i] = regimeType::transient;
 
                         // No convective heat transfer
                         alphatConv_[i] = 0.0;
@@ -885,44 +987,53 @@ void alphatWallBoilingWallFunctionFvPatchScalarField::updateCoeffs()
                         // transient boiling
                         dmdt_[i] =
                                 fLiquid[i]
-                                *(
+                               *(
                                     relax*Qtb[i]*AbyV_[i]/L[i]
-                                    + (1 - relax)*dmdt_[i]
+                                 + (1 - relax)*dmdt_[i]
                                 );
 
-                        mDotL_[i] = dmdt_[i]*L[i];
 
+                        mDotL_[i] = dmdt_[i]*L[i];
 
                         // No quenching flux
                         qq_[i] = 0.0;
 
                         this->operator[](i) =
-                        max
-                        (
+                            max
                             (
-                                mDotL_[i]/AbyV_[i]
-                                /max(hewSn[i], scalar(1e-16))
-                            )/max(liquidw[i], scalar(1e-8)),
-                            1e-8
-                        );
+                                (
+                                    mDotL_[i]/AbyV_[i]
+                                   /max(hewSn[i], scalar(1e-16))
+                                )/max(liquidw[i], scalar(1e-8)),
+                                scalar(1e-8)
+                            );
+
+                        if (debug == 2)
+                        {
+                            Info<< "Transient boiling:  " <<  nl
+                                << "  fraction Liq: " << fLiquid[i] << nl
+                                << "  Heat flux: " << Qtb[i] << nl
+                                << "  delta Tsub: " << (Tw[i] - Tsatw[i])
+                                << endl;
+                        }
 
                     }
                     else if (Tw[i] > TLeiden[i])
                     {
-                        regimeTypes[i] = regimeType::film; // film boiling
+                        regimeTypes_[i] = regimeType::film; // film boiling
 
                         // No convective heat transfer
                         alphatConv_[i] = 0.0;
 
                         // Film boiling
                         dmdt_[i] =
-                            fLiquid[i]
-                            *(
-                                relax*htcFilmBoiling[i]
-                                *max(Tw[i] - Tsatw[i], 0)
-                                *AbyV_[i]/L[i]
-                                + (1 - relax)*dmdt_[i]
-                            );
+                                fLiquid[i]
+                               *(
+                                    relax*htcFilmBoiling[i]
+                                    *max(Tw[i] - Tsatw[i], scalar(0))
+                                    *AbyV_[i]/L[i]
+                                  + (1 - relax)*dmdt_[i]
+                                );
 
 
                         mDotL_[i] = dmdt_[i]*L[i];
@@ -939,17 +1050,27 @@ void alphatWallBoilingWallFunctionFvPatchScalarField::updateCoeffs()
                         // alphaFilm in the coupled BC. We subtract
                         // alpha and divide by phase to get a net alphaFilm
                         this->operator[](i) =
-                            (
-                                alphaFilm[i]/max(liquidw[i], scalar(1e-8))
-                              - alphaw[i]
-                            );
+                        (
+                            alphaFilm[i]/max(liquidw[i], scalar(1e-8))
+                          - alphaw[i]
+                        );
+
+                        if (debug == 2)
+                        {
+                            Info<< "Film boiling:  " <<  nl
+                                << "  fraction Liq: " << fLiquid[i] << nl
+                                << "  Heat flux: "
+                                << htcFilmBoiling[i]*(Tw[i] - Tsatw[i]) << nl
+                                << "  delta Tsub: " << (Tw[i] - Tsatw[i])
+                                << endl;
+                        }
                     }
                 }
                 else
                 {
                     // Tw below Tsat. No boiling single phase convection
                     // Single phase
-                    regimeTypes[i] = regimeType::nonBoiling;
+                    regimeTypes_[i] = regimeType::nonBoiling;
 
                     qq_[i] = 0.0;
                     mDotL_[i] = 0.0;
@@ -962,7 +1083,7 @@ void alphatWallBoilingWallFunctionFvPatchScalarField::updateCoeffs()
                         (
                             fLiquid[i]*(alphatConv_[i])
                             /max(liquidw[i], scalar(1e-8)),
-                            1e-8
+                            scalar(1e-8)
                         )
                     );
                 }
@@ -972,11 +1093,15 @@ void alphatWallBoilingWallFunctionFvPatchScalarField::updateCoeffs()
             {
                 const scalarField qEff
                 (
-                    liquidw*(*this + alphaw)*hew.snGrad()
+                    fLiquid*liquidw*(*this + alphaw)*hew.snGrad()
                 );
 
                 Info<< "alphat for liquid:  " <<  nl << endl;
 
+                Info<< "  qEffLiq: " << gMin(qEff) << " - "
+                    << gMax(qEff) << endl;
+
+
                 Info<< "  alphatl: " << gMin((*this)) << " - "
                     << gMax((*this)) << endl;
 
@@ -990,7 +1115,7 @@ void alphatWallBoilingWallFunctionFvPatchScalarField::updateCoeffs()
                 Info<< " Effective heat transfer rate to liquid: " << Qeff
                     << endl << nl;
 
-                if (debug & 2)
+                if (debug == 2)
                 {
                     scalarField nSubCools(this->size(), 0);
                     scalarField nTransients(this->size(), 0);
@@ -1000,7 +1125,7 @@ void alphatWallBoilingWallFunctionFvPatchScalarField::updateCoeffs()
                     forAll(*this, i)
                     {
                         //faceRegimes[i] = regimeTypes[i];
-                        switch (regimeTypes[i])
+                        switch (regimeTypes_[i])
                         {
                             case regimeType::subcool:
                                 nSubCools[i] = 1;
@@ -1037,7 +1162,7 @@ void alphatWallBoilingWallFunctionFvPatchScalarField::updateCoeffs()
 
                     const scalarField qc
                     (
-                        nNonBoilings*fLiquid*A1*(alphatConv_ + alphaw)
+                        nNonBoilings*fLiquid*(alphatConv_ + alphaw)
                         *hew.snGrad()
                     );
 
@@ -1070,26 +1195,15 @@ void alphatWallBoilingWallFunctionFvPatchScalarField::updateCoeffs()
                     (
                         fLiquid*nSubCools*
                         (
-                            A1*alphatConv_*hew.snGrad()
-                            + qe() + qq()
+                            A1*alphatConv_*hew.snGrad() + qe() + qq()
                         )
                     );
 
-
                     scalar QsubCool = gSum(qSubCool*patch().magSf());
 
                     Info<< " Sub Cool boiling heat transfer: " << QsubCool
                         << endl;
 
-                    Info<< "  N: " << gMin(nSubCools*N) << " - "
-                        << gMax(nSubCools*N) << endl;
-                    Info<< "  dDep: " << gMin(nSubCools*dDep_) << " - "
-                        << gMax(nSubCools*dDep_) << endl;
-                    Info<< "  fDep: " << gMin(nSubCools*fDep) << " - "
-                        << gMax(nSubCools*fDep) << endl;
-                    Info<< "  A1: " << gMin(nSubCools*A1) << " - "
-                        << gMax(nSubCools*A1) << endl;
-
                     Info<< nl;
                 }
             }
@@ -1125,17 +1239,33 @@ void alphatWallBoilingWallFunctionFvPatchScalarField::write(Ostream& os) const
             break;
         case liquidPhase:
         {
-            os.beginBlock("nucleationSiteModel");
-            nucleationSiteModel_->write(os);
-            os.endBlock();
+            if (nucleationSiteModel_)
+            {
+                os.beginBlock("nucleationSiteModel");
+                nucleationSiteModel_->write(os);
+                os.endBlock();
+            }
 
-            os.beginBlock("departureDiamModel");
-            departureDiamModel_->write(os);
-            os.endBlock();
+            if (departureDiamModel_)
+            {
+                os.beginBlock("departureDiamModel");
+                departureDiamModel_->write(os);
+                os.endBlock();
+            }
 
-            os.beginBlock("departureFreqModel");
-            departureFreqModel_->write(os);
-            os.endBlock();
+            if (departureFreqModel_)
+            {
+                os.beginBlock("departureFreqModel");
+                departureFreqModel_->write(os);
+                os.endBlock();
+            }
+
+            if (nucleatingModel_)
+            {
+                os.beginBlock("nucleateFluxModel");
+                nucleatingModel_->write(os);
+                os.endBlock();
+            }
 
             if (filmBoilingModel_)
             {
@@ -1181,6 +1311,7 @@ void alphatWallBoilingWallFunctionFvPatchScalarField::write(Ostream& os) const
 
             os.writeEntry("K", K_);
             os.writeEntry("wp", wp_);
+            os.writeEntry("liquidTatYplus", liquidTatYplus_);
             break;
         }
     }
@@ -1191,6 +1322,7 @@ void alphatWallBoilingWallFunctionFvPatchScalarField::write(Ostream& os) const
     dDep_.writeEntry("dDep", os);
     qq_.writeEntry("qQuenching", os);
     alphatConv_.writeEntry("alphatConv", os);
+
     writeEntry("value", os);
 }
 

src/phaseSystemModels/reactingEuler/multiphaseSystem/derivedFvPatchFields/alphatWallBoilingWallFunction/alphatWallBoilingWallFunctionFvPatchScalarField.H

@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2015-2018 OpenFOAM Foundation
-    Copyright (C) 2018 OpenCFD Ltd
+    Copyright (C) 2018-2021 OpenCFD Ltd
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -31,27 +31,30 @@ Description
     A thermal wall function for simulation of boiling wall.
 
     This alpha wall function can handle the following regimes:
-        single phase
-        subcooled nucleate wall boiling
-        transitional boiling
-        film boiling.
+      - single phase
+      - subcooled nucleate wall boiling
+      - transitional boiling
+      - film boiling
 
     The wall function uses a partition method to transfer heat either
     to the liquid or vapor phase. At the moment, this function works
-    in a wall temperature fixed mode. i.e, there is no consideration
+    in a wall temperature fixed mode, i.e. there is no consideration
     for the sudden change of heat transfer coefficient (htc) after
     reaching TDBN (deviation from nucleate boiling temperature).
 
-     References:
+    References:
     \verbatim
-        Numerical simulation of immersion quenching process of an engine cylinder head
-        Vedanth Srinivasan, Kil-Min Moon, David Greif, De Ming Wang, Myung-hwan Kim
-        Applied Mathematical Modelling 34 (2010) 2111-2128
+        Srinivasan, V., Moon, K. M., Greif, D.,
+        Wang, D. M., & Kim, M. H. (2010).
+        Numerical simulation of immersion quenching
+        process of an engine cylinder head.
+        Applied Mathematical Modelling, 34(8), 2111-2128.
+        DOI:10.1016/j.apm.2009.10.023
     \endverbatim
 
 
     For the single phase non-boiling regime the standard
-    JayatillekeWallFunction is used.
+    \c JayatillekeWallFunction is used.
 
     For the sub-cool nucleate boiling regime the following runtime
     selectable submodels are used:
@@ -66,21 +69,23 @@ Description
 
     References:
     \verbatim
-        "On the modeling of multidimensional effects in boiling channels"
-        Kurul, N., Podowski, M.Z.,
-        ANS Proceedings, National Heat Transfer Conference,
-        Minneapolis, Minnesota, USA, July 28-31, 1991,
+        Kurul, N., & Podowski, M. Z. (1991).
+        On the modeling of multidimensional effects in boiling channels.
+        Proceedings of the 27th National Heat Transfer Conference.
+        Minneapolis, Minn, USA, July 28-31, 1991.
         ISBN: 0-89448-162-1, pp. 30-40
+
+        Peltola, J., & Pättikangas, T. (2012).
+        Development and validation of a boiling model
+        for OpenFOAM multiphase solver.
+        Proceedings of the CFD4NRS-4. p. 59.
+        Daejeon, Democratic People's Republic of Korea, September 10-12, 2012.
     \endverbatim
 
-    \verbatim
-        "Development and validation of a boiling model for OpenFOAM
-        multiphase solver"
-        Peltola, J., Pättikangas, T.J.H.,
-        CFD4NRS-4 Conference Proceedings, paper 59,
-        Daejeon, Korea, September 10-12 2012
-    \endverbatim
 
+    Alternatively a correlation can be used instead of the RPI wall boiling model.
+    If the keyword nucleatingModel a model is provided the BC uses it
+    instead of the RPI model.
 
     The transition boiling regime flux (TBF) is modelled following
     a temperature based linear interpolation between the critical heat flux
@@ -89,131 +94,180 @@ Description
     (TLeiden) a linear interpolation is used between CHF and MHF.
 
     Thus, the following models are required:
-        LeidenfrostModel
-        CHFModel
-        CHFSubCoolModel
-        MHFModel
-        TDNBModel
-        filmBoilingModel
+      - LeidenfrostModel
+      - CHFModel
+      - CHFSubCoolModel
+      - MHFModel
+      - TDNBModel
+      - filmBoilingModel
 
     The linear interpolation is as follows:
 
-        TBF = CHF*phi + (1 - phi)*MHF
+    \f[
+        TBF = CHF*\phi + (1 - \phi)*MHF
+    \f]
 
-        where phi:
+    with
+    \f[
+        \phi = w_p*(T_w - T_{DNB})/(T_{Leiden} - T_{DNB})
+    \f]
 
-            phi = wp*(Tw - TDNB)/(TLeiden - TDNB),
-
-        where:
-            wp model constant
-            Tw wall temperature
+    where:
+    \vartable
+        w_p | Model constant
+        T_w | Wall temperature [K]
+    \endvartable
 
 
-    The film boiling regime is applied when Tw is larger than TLeiden. In
-    this regime the corrlation from the filmBoilingModel is used for
-    calculating the cht from the wall.
+    The film boiling regime is applied when \f$T_w\f$ is larger than
+    \f$T_{Leiden}\f$. In this regime the correlation from the
+    \c filmBoilingModel is used for calculating the cht from the wall.
 
-    The filmBoilingModel is needed in the vapor field in order to calculate
+    The \c filmBoilingModel is needed in the vapor field in order to calculate
     the heat transfer to the vapor phase in film boiling regime.
 
 
 Usage
-    \table
-        Property     | Description             | Required    | Default value
-        phaseType    | 'vapor' or 'liquid'     | yes         |
-        relax        |wall boiling model relaxation| yes     |
-        Prt          | inherited from alphatPhaseChangeJayatillekeWallFunction
-        Cmu          | inherited from alphatPhaseChangeJayatillekeWallFunction
-        kappa        | inherited from alphatPhaseChangeJayatillekeWallFunction
-        E            | inherited from alphatPhaseChangeJayatillekeWallFunction
-        dmdt         | phase change mass flux  | no         |
-        value        | initial alphat value    | yes         |
-
-        if phaseType 'vapor':
-
-        partitioningModel|                     | yes        |
-        filmBoilingModel |                     | yes        |
-        LeidenfrostModel |                     | yes        |
-
-        if phaseType 'liquid':
-
-        partitioningModel|                     | yes        |
-        nucleationSiteModel|                     | yes      |
-        departureDiamModel|                     | yes       |
-        departureFreqModel|                     | yes       |
-        K                | bubbles area constant| no       | 4
-
-        LeidenfrostModel |                      | no       |
-        CHFModel         |                      | no       |
-        CHFSubCoolModel  |                      | no       |
-        MHFModel         |                      | no       |
-        TDNBModel        |                      | no       |
-        filmBoilingModel |                      | no       |
-        wp               |                      | no       | 1
-    \endtable
-
-    NOTE: Runtime selectabale submodels may require model specific entries
-
-    Example usage:
+    Example of the boundary condition specification:
     \verbatim
-    hotWall
+    <patchName>
     {
+        // Mandatory entries
         type            compressible::alphatWallBoilingWallFunction;
-        phaseType       liquid;
-        Prt             0.85;
-        Cmu             0.09;
-        kappa           0.41;
-        E               9.8;
-        relax           0.1;
-        dmdt            uniform 0;
+        phaseType       <word>;
+        otherPhase      <word>;
+        relax           <Function1<scalar>>;
+
         partitioningModel
         {
             type        Lavieville;
             alphaCrit   0.2;
         }
-        nucleationSiteModel
-        {
-            type        LemmertChawla;
-        }
-        departureDiamModel
-        {
-            type        TolubinskiKostanchuk;
-        }
-        departureFreqModel
-        {
-            type        Cole;
-        }
 
-        LeidenfrostModel
-        {
-            type        Spiegler;
-            Tcrit       647;
-        }
-        CHFModel
-        {
-            type        Zuber;
-        }
-        CHFSubCoolModel
-        {
-            type        HuaXu;
-            Kburn       0.5;
-        }
-        MHFModel
-        {
-            type        Jeschar;
-            Kmhf        1;
-        }
-        TDNBModel
-        {
-            type        Schroeder;
-        }
-        filmBoilingModel
-        {
-            type        Bromley;
-        }
-        value           uniform 0.01;
+        // Conditional entries
+
+            // Option-1: phaseType=vapor
+
+                // Optional entries
+                LeidenfrostModel
+                {
+                    type        Spiegler;
+                    Tcrit       647;
+                }
+
+                filmBoilingModel
+                {
+                    type        Bromley;
+                }
+
+
+            // Option-2: phaseType=liquid
+            nucleationSiteModel
+            {
+                type        LemmertChawla;
+            }
+
+            departureDiamModel
+            {
+                type        TolubinskiKostanchuk;
+            }
+
+            departureFreqModel
+            {
+                type        Cole;
+            }
+
+                // Optional entries
+                LeidenfrostModel
+                {
+                    type        Spiegler;
+                    Tcrit       647;
+                }
+
+                CHFModel
+                {
+                    type        Zuber;
+                }
+
+                CHFSubCoolModel
+                {
+                    type        HuaXu;
+                    Kburn       0.5;
+                }
+
+                MHFModel
+                {
+                    type        Jeschar;
+                    Kmhf        1;
+                }
+
+                TDNBModel
+                {
+                    type        Schroeder;
+                }
+
+                filmBoilingModel
+                {
+                    type        Bromley;
+                }
+
+                dDep    <scalarField>;
+                K       <scalar>;
+                wp      <scalar>;
+                qQuenching <scalarField>;
+
+
+        // Optional entries
+        alphatConv      <scalarField>;
+
+        //Inherited entries
+        ...
     \endverbatim
 
+    where the entries mean:
+    \table
+      Property  | Description                            | Type | Reqd | Deflt
+      type      | compressible::alphatWallBoilingWallFunction | word | yes | -
+      phaseType | Name of phase type                     | word | yes  | -
+      otherPhase | Name of other phase                   | word | yes  | -
+      relax     | Relaxation factor for dmdt | Function1\<scalar\> <!--
+                -->                                              | yes | -
+      alphatConv | Convective turbulent thermal diffusivity    <!--
+                -->                                 | scalarField | no | 0
+      partitioningModel | Run-time selected heat flux partitioning model <!--
+                -->                                      | dict | yes  | -
+    \endtable
+
+    Options for the \c phaseType and \c otherPhase entries:
+    \verbatim
+      vapor       | Vapor phase
+      liquid      | Liquid phase
+    \endverbatim
+
+    when \c phaseType=liquid:
+    \table
+      Property  | Description                             | Type | Reqd | Deflt
+      nucleationSiteModel | Nucleation site density model | dict | yes  | -
+      departureDiamModel  | Bubble departure diameter model     <!--
+                        -->                               | dict | yes  | -
+      departureFreqModel | Bubble departure frequency model | dict | yes | -
+      LeidenfrostModel   | Leidenfrost temperature model  | dict | no   | -
+      CHFModel | Critical heat flux model                 | dict | no   | -
+      CHFSubCoolModel | CHF sub-cool model                | dict | no   | -
+      MHFModel        | Minium heat flux model            | dict | no   | -
+      TDNBModel       | Departure from nulceate boiling model | dict | no | -
+      filmBoilingModel | Film boiling model               | dict | no   | -
+      K | Model constant for area of bubbles            | scalar | no   | 4.0
+      wp | Wetting parameter for transient boiling      | scalar | no   | 1.0
+    \endtable
+
+    The inherited entries are elaborated in:
+     -\link alphatPhaseChangeJayatillekeWallFunctionFvPatchScalarField.H\endlink
+
+Notes
+  - Runtime selectabale submodels may require model specific entries
+  - \c phaseType and \c otherPhase entries should be the opposite of each other.
+
 See also
     Foam::alphatPhaseChangeJayatillekeWallFunctionFvPatchField
 
@@ -232,6 +286,7 @@ SourceFiles
 #include "nucleationSiteModel.H"
 #include "departureDiameterModel.H"
 #include "departureFrequencyModel.H"
+#include "nucleateFluxModel.H"
 
 #include "LeidenfrostModel.H"
 #include "filmBoilingModel.H"
@@ -248,7 +303,7 @@ namespace compressible
 {
 
 /*---------------------------------------------------------------------------*\
-            Class alphatWallBoilingWallFunctionFvPatchScalarField Declaration
+        Class alphatWallBoilingWallFunctionFvPatchScalarField Declaration
 \*---------------------------------------------------------------------------*/
 
 class alphatWallBoilingWallFunctionFvPatchScalarField
@@ -257,7 +312,7 @@ class alphatWallBoilingWallFunctionFvPatchScalarField
 {
 public:
 
-    // Data types
+    // Public Enumerations
 
         //- Enumeration listing the possible operational modes
         enum phaseType
@@ -269,7 +324,7 @@ public:
 
 private:
 
-    // Private data
+    // Private Data
 
         //- Enumeration of regimes per face
         enum regimeType
@@ -280,16 +335,16 @@ private:
             nonBoiling
         };
 
-        //- name of the other phase (vapor/liquid phase)
+        //- Name of the other phase (vapor/liquid phase)
         word otherPhaseName_;
 
-        //- Heat source type names
+        //- Names of heat source types
         static const Enum<phaseType> phaseTypeNames_;
 
         //- Heat source type
         phaseType phaseType_;
 
-        //- dmdt relaxationFactor
+        //- Relaxation factor for dmdt
         autoPtr<Function1<scalar>> relax_;
 
         //- Patch face area by cell volume
@@ -325,6 +380,10 @@ private:
             autoPtr<wallBoilingModels::departureFrequencyModel>
                 departureFreqModel_;
 
+            //- Run-time sub-cooling heat flux correlatiom
+            autoPtr<wallBoilingModels::nucleateFluxModel>
+                nucleatingModel_;
+
 
         // Film boiling model
 
@@ -332,6 +391,7 @@ private:
             autoPtr<wallBoilingModels::filmBoilingModel>
                 filmBoilingModel_;
 
+
         // Transition boiling model
 
             //- Run-time selected for Leidenfrost temperature
@@ -353,6 +413,12 @@ private:
             //- Wetting parameter for transient boiling
             scalar wp_;
 
+            //- Use Liquid temperature at y+=250
+            bool liquidTatYplus_;
+
+            //- Face regime
+            labelField regimeTypes_;
+
 
 public:
 
@@ -378,8 +444,8 @@ public:
         );
 
         //- Construct by mapping given
-        //  alphatWallBoilingWallFunctionFvPatchScalarField
-        //  onto a new patch
+        //- alphatWallBoilingWallFunctionFvPatchScalarField
+        //- onto a new patch
         alphatWallBoilingWallFunctionFvPatchScalarField
         (
             const alphatWallBoilingWallFunctionFvPatchScalarField&,
@@ -423,7 +489,7 @@ public:
         }
 
 
-    // Member functions
+    // Member Functions
 
         using alphatPhaseChangeWallFunctionFvPatchScalarField::dmdt;
 
@@ -454,7 +520,13 @@ public:
             return mDotL_/AbyV_;
         }
 
-        // Evaluation functions
+        //- Return const reference to the face regime
+        const labelField& regimeTypes() const noexcept
+        {
+            return regimeTypes_;
+        }
+
+        // Evaluation
 
             //- Update the coefficients associated with the patch field
             virtual void updateCoeffs();

src/phaseSystemModels/reactingEuler/multiphaseSystem/derivedFvPatchFields/copiedFixedValue/copiedFixedValueFvPatchScalarField.H

@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2015-2018 OpenFOAM Foundation
-    Copyright (C) 2020 OpenCFD Ltd.
+    Copyright (C) 2020-2021 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -30,6 +30,30 @@ Class
 Description
     Copies the boundary values from a user specified field.
 
+Usage
+    Example of the boundary condition specification:
+    \verbatim
+    <patchName>
+    {
+        // Mandatory entries
+        type                copiedFixedValue;
+        sourceFieldName     <word>;
+
+        // Inherited entries
+        ...
+    }
+    \endverbatim
+
+    where the entries mean:
+    \table
+      Property  | Description                       | Type | Reqd | Deflt
+      type      | Type name: copiedFixedValue       | word | yes  | -
+      sourceFieldName | Name of the source field    | word | yes  | -
+    \endtable
+
+    The inherited entries are elaborated in:
+      - \link fixedValueFvPatchFields.H \endlink
+
 See also
     Foam::fixedValueFvPatchField
 
@@ -60,6 +84,7 @@ protected:
 
     // Protected Data
 
+        //- Name of the source field
         word sourceFieldName_;
 
 public:
@@ -86,8 +111,8 @@ public:
         );
 
         //- Construct by mapping given
-        //  copiedFixedValueFvPatchScalarField
-        //  onto a new patch
+        //- copiedFixedValueFvPatchScalarField
+        //- onto a new patch
         copiedFixedValueFvPatchScalarField
         (
             const copiedFixedValueFvPatchScalarField&,
@@ -133,7 +158,7 @@ public:
 
     // Member Functions
 
-        // Evaluation Functions
+        // Evaluation
 
             //- Update the coefficients associated with the patch field
             virtual void updateCoeffs();

src/phaseSystemModels/reactingEuler/multiphaseSystem/derivedFvPatchFields/fixedMultiPhaseHeatFlux/fixedMultiPhaseHeatFluxFvPatchScalarField.C

@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2015-2020 OpenFOAM Foundation
-    Copyright (C) 2020 OpenCFD Ltd.
+    Copyright (C) 2020-2021 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -144,7 +144,7 @@ void Foam::fixedMultiPhaseHeatFluxFvPatchScalarField::updateCoeffs()
 
         if (debug)
         {
-            scalarField q0(T.snGrad()*alpha*kappaEff);
+            const scalarField q0(T.snGrad()*alpha*kappaEff);
             Q += q0;
 
             Info<< patch().name() << " " << phase.name()
@@ -162,7 +162,7 @@ void Foam::fixedMultiPhaseHeatFluxFvPatchScalarField::updateCoeffs()
             << gSum(patch().magSf()*Q) << " W" << endl;
     }
 
-    operator==((1 - relax_)*Tp + relax_*max(Tmin_,(q_ + A)/(B)));
+    operator==((scalar(1) - relax_)*Tp + relax_*max(Tmin_,(q_ + A)/(B)));
 
     fixedValueFvPatchScalarField::updateCoeffs();
 }

src/phaseSystemModels/reactingEuler/multiphaseSystem/derivedFvPatchFields/fixedMultiPhaseHeatFlux/fixedMultiPhaseHeatFluxFvPatchScalarField.H

@@ -6,6 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2015-2020 OpenFOAM Foundation
+    Copyright (C) 2021 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -27,13 +28,44 @@ Class
     Foam::fixedMultiPhaseHeatFluxFvPatchScalarField
 
 Description
-    Calculates a wall temperature that produces the specified overall wall heat
-    flux across all the phases in an Eulerian multi-phase simulation.
+    Calculates a wall temperature that produces
+    the specified overall wall heat flux across
+    all the phases in an Eulerian multi-phase simulation.
 
-    Intended to be used with copiedFixedValue to ensure that phase wall
-    temperature are consistent:
-        - Set 'fixedMultiPhaseHeatFlux' boundary for one of the phases
-        - Use 'copiedFixedValue' for all the other phases.
+    Intended to be used with \c copiedFixedValue
+    to ensure that phase wall temperature are consistent:
+        - Set \c fixedMultiPhaseHeatFlux boundary for one of the phases
+        - Use \c copiedFixedValue for all the other phases.
+
+Usage
+    Example of the boundary condition specification:
+    \verbatim
+    <patchName>
+    {
+        // Mandatory entries
+        type            fixedMultiPhaseHeatFlux;
+        q               <scalarField>;
+
+        // Optional entries
+        relax           <scalar>;
+        Tmin            <scalar>;
+
+        // Inherited entries
+        ...
+    }
+    \endverbatim
+
+    where the entries mean:
+    \table
+      Property  | Description                         | Type | Reqd | Deflt
+      type      | Type name: fixedMultiPhaseHeatFlux  | word | yes  | -
+      q         | Heat power [W] or flux [W/m2]       | scalarField | yes | -
+      relax     | Relaxation factor                   | scalar | no | 1.0
+      Tmin      | Minimum temperature limit [K]       | scalar | no | 273.0
+    \endtable
+
+    The inherited entries are elaborated in:
+      - \link fixedValueFvPatchFields.H \endlink
 
 See also
     Foam::fixedValueFvPatchField
@@ -97,8 +129,8 @@ public:
         );
 
         //- Construct by mapping given
-        //  fixedMultiPhaseHeatFluxFvPatchScalarField
-        //  onto a new patch
+        //- fixedMultiPhaseHeatFluxFvPatchScalarField
+        //- onto a new patch
         fixedMultiPhaseHeatFluxFvPatchScalarField
         (
             const fixedMultiPhaseHeatFluxFvPatchScalarField&,
@@ -144,7 +176,7 @@ public:
 
     // Member Functions
 
-        // Mapping functions
+        // Mapping
 
             //- Map (and resize as needed) from self given a mapping object
             //  Used to update fields following mesh topology change
@@ -155,7 +187,7 @@ public:
             virtual void rmap(const fvPatchScalarField&, const labelList&);
 
 
-        // Evaluation Functions
+        // Evaluation
 
             //- Update the coefficients associated with the patch field
             virtual void updateCoeffs();

src/phaseSystemModels/reactingEuler/multiphaseSystem/derivedFvPatchFields/wallBoilingSubModels/CHFModels/CHFModel/CHFModel.H

@@ -5,7 +5,7 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2018-2020 OpenCFD Ltd
+    Copyright (C) 2018-2021 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -27,7 +27,8 @@ Class
     Foam::wallBoilingModels::CHFModel
 
 Description
-    Base class for nucleation site density models
+    Base class for critical heat flux (CHF)
+    correlation models for boiling flows.
 
 SourceFiles
     CHFModel.C
@@ -103,7 +104,10 @@ public:
             const scalarField& L
         ) const = 0;
 
-        virtual void write(Ostream& os) const;
+        // I-O
+
+            //- Write
+            virtual void write(Ostream& os) const;
 };
 
 

src/phaseSystemModels/reactingEuler/multiphaseSystem/derivedFvPatchFields/wallBoilingSubModels/CHFModels/Zuber/Zuber.C

@@ -5,7 +5,7 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2018-2020 OpenCFD Ltd
+    Copyright (C) 2018-2021 OpenCFD Ltd
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -63,12 +63,6 @@ Foam::wallBoilingModels::CHFModels::Zuber::Zuber
 {}
 
 
-// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
-
-Foam::wallBoilingModels::CHFModels::Zuber::~Zuber()
-{}
-
-
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
 Foam::tmp<Foam::scalarField>
@@ -82,12 +76,18 @@ Foam::wallBoilingModels::CHFModels::Zuber::CHF
     const scalarField& L
 ) const
 {
-    const uniformDimensionedVectorField& g =
+    const auto& g =
         liquid.mesh().time().lookupObject<uniformDimensionedVectorField>("g");
 
-    const scalarField rhoVapor(vapor.thermo().rho(patchi));
-    const scalarField rhoLiq(liquid.thermo().rho(patchi));
+    const labelUList& cells = liquid.mesh().boundary()[patchi].faceCells();
 
+    const scalarField& pw = liquid.thermo().p().boundaryField()[patchi];
+
+    tmp<scalarField> trhoVapor = vapor.thermo().rhoEoS(pw ,Tsatw, cells);
+    const scalarField& rhoVapor = trhoVapor.ref();
+
+    tmp<scalarField> trhoLiq = liquid.thermo().rhoEoS(pw, Tsatw, cells);
+    const scalarField& rhoLiq = trhoLiq.ref();
     const phasePairKey pair(liquid.name(), vapor.name());
     const scalarField sigma
     (

src/phaseSystemModels/reactingEuler/multiphaseSystem/derivedFvPatchFields/wallBoilingSubModels/CHFModels/Zuber/Zuber.H

@@ -5,7 +5,7 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2018 OpenCFD Ltd
+    Copyright (C) 2018-2021 OpenCFD Ltd
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -24,17 +24,40 @@ License
     along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
 
 Class
-    Foam::wallBoilingModels:CHFModels:::Zuber
+    Foam::wallBoilingModels::CHFModels::Zuber
 
 Description
-    Critical heat flux (CHF) correlation
+    A critical heat flux (CHF) correlation model
+    based on Zuber (1958) for boiling flows.
 
-    References:
+    Reference:
     \verbatim
-        N. Zuber, On the stability of boiling heat transfer,
-        Trans. ASME 80 (1958) 711
+        Zuber, N. (1958).
+        On the stability of boiling heat transfer.
+        Trans. Am. Soc. Mech. Engrs., 80.
+        URL:https://www.osti.gov/biblio/4326542
     \endverbatim
 
+Usage
+    Example of the model specification:
+    \verbatim
+    CHFModel
+    {
+        // Mandatory entries
+        type            Zuber;
+
+        // Optional entries
+        Cn              <scalar>;
+    }
+    \endverbatim
+
+    where the entries mean:
+    \table
+      Property  | Description                         | Type | Reqd | Deflt
+      type      | Type name: Zuber                    | word | yes  | -
+      Cn        | Model coefficient                 | scalar | no   | 0.131
+    \endtable
+
 SourceFiles
     Zuber.C
 
@@ -62,17 +85,27 @@ class Zuber
 :
     public CHFModel
 {
+    // Private Data
 
-    // Private data:
-
-        //- Coefficient constant
+        //- Model coefficient
         scalar Cn_;
 
+
+    // Private Member Functions
+
+        //- No copy construct
+        Zuber(const Zuber&) = delete;
+
+        //- No copy assignment
+        void operator=(const Zuber&) = delete;
+
+
 public:
 
     //- Runtime type information
     TypeName("Zuber");
 
+
     // Constructors
 
         //- Construct from a dictionary
@@ -80,7 +113,7 @@ public:
 
 
     //- Destructor
-    virtual ~Zuber();
+    virtual ~Zuber() = default;
 
 
     // Member Functions
@@ -96,7 +129,10 @@ public:
             const scalarField& L
         ) const;
 
-        virtual void write(Ostream& os) const;
+        // I-O
+
+            // Write
+            virtual void write(Ostream& os) const;
 };
 
 

src/phaseSystemModels/reactingEuler/multiphaseSystem/derivedFvPatchFields/wallBoilingSubModels/CHFSubCoolModels/CHFSubCoolModel/CHFSubCoolModel.H

@@ -5,7 +5,7 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2018-2020 OpenCFD Ltd
+    Copyright (C) 2018-2021 OpenCFD Ltd
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -27,7 +27,8 @@ Class
     Foam::CHFModels::CHFSubCoolModel
 
 Description
-    Base class for nucleation site density models
+    Base class for critical heat flux (CHF)
+    sub-cooling correlation models for boiling flows.
 
 SourceFiles
     CHFSubCoolModel.C
@@ -103,7 +104,10 @@ public:
             const scalarField& L
         ) const = 0;
 
-        virtual void write(Ostream& os) const;
+        // I-O
+
+            //- Write
+            virtual void write(Ostream& os) const;
 };
 
 

src/phaseSystemModels/reactingEuler/multiphaseSystem/derivedFvPatchFields/wallBoilingSubModels/CHFSubCoolModels/HuaXu/HuaXu.C

@@ -5,7 +5,7 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2018-2020 OpenCFD Ltd
+    Copyright (C) 2018-2021 OpenCFD Ltd
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -63,12 +63,6 @@ Foam::wallBoilingModels::CHFModels::HuaXu::HuaXu
 {}
 
 
-// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
-
-Foam::wallBoilingModels::CHFModels::HuaXu::~HuaXu()
-{}
-
-
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
 Foam::tmp<Foam::scalarField>
@@ -82,13 +76,22 @@ Foam::wallBoilingModels::CHFModels::HuaXu::CHFSubCool
     const scalarField& L
 ) const
 {
-    const uniformDimensionedVectorField& g =
+    const auto& g =
         liquid.mesh().time().lookupObject<uniformDimensionedVectorField>("g");
 
     const scalarField alphaLiq(liquid.alpha(patchi));
 
-    const scalarField rhoVapor(vapor.thermo().rho(patchi));
-    const scalarField rhoLiq(liquid.thermo().rho(patchi));
+    const labelUList& cells = liquid.mesh().boundary()[patchi].faceCells();
+
+    const scalarField& pw = liquid.thermo().p().boundaryField()[patchi];
+
+    tmp<scalarField> trhoVapor = vapor.thermo().rhoEoS(pw, Tsatw, cells);
+    const scalarField& rhoVapor = trhoVapor.ref();
+
+    tmp<scalarField> trhoLiq = liquid.thermo().rhoEoS(pw, Tsatw, cells);
+    const scalarField& rhoLiq = trhoLiq.ref();
+
+
     tmp<volScalarField> tCp = liquid.thermo().Cp();
     const volScalarField& Cp = tCp();
     const fvPatchScalarField& Cpw = Cp.boundaryField()[patchi];
@@ -105,7 +108,7 @@ Foam::wallBoilingModels::CHFModels::HuaXu::CHFSubCool
        /
         (
             alphaLiq
-          * pow(mag(g.value())*(rhoLiq-rhoVapor), 0.25)
+          * pow025(mag(g.value())*(rhoLiq-rhoVapor))
           * sqrt(rhoVapor)
         )
     );
@@ -115,8 +118,7 @@ Foam::wallBoilingModels::CHFModels::HuaXu::CHFSubCool
         rhoLiq*Cpw*max(Tsatw - Tl, scalar(0))/(rhoVapor*L)
     );
 
-    return
-        Kburn_*(1 + 0.345*Ja/pow(Pe, 0.25));
+    return Kburn_*(scalar(1) + 0.345*Ja/pow025(Pe));
 }
 
 

src/phaseSystemModels/reactingEuler/multiphaseSystem/derivedFvPatchFields/wallBoilingSubModels/CHFSubCoolModels/HuaXu/HuaXu.H

@@ -5,7 +5,7 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2018 OpenCFD Ltd
+    Copyright (C) 2018-2021 OpenCFD Ltd
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -24,19 +24,41 @@ License
     along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
 
 Class
-    Foam::wallBoilingModels:CHFModels:::HuaXu
+    Foam::wallBoilingModels::CHFModels::HuaXu
 
 Description
+    A critical heat flux (CHF) sub-cooling correlation model
+    based on Hua-Xu (2000) for boiling flows.
 
-    Critical heat flux for soob cool boiling flows.
-
-    References:
+    Reference:
     \verbatim
-        T.C. Hua, J.J. Xu, Quenching boiling in subcooled liquid nitrogen
-        for solidification of aqueous materials, Mater.
-        Sci. Eng. A 292 (2000) 169–172.
+        Hua, T. C., & Xu, J. J. (2000).
+        Quenching boiling in subcooled liquid nitrogen
+        for solidification of aqueous materials.
+        Materials Science and Engineering: A, 292(2), 169-172.
+        DOI:10.1016/S0921-5093(00)01004-2
     \endverbatim
 
+Usage
+    Example of the model specification:
+    \verbatim
+    CHFSubCoolModel
+    {
+        // Mandatory entries
+        type            HuaXu;
+
+        // Optional entries
+        Kburn           <scalar>;
+    }
+    \endverbatim
+
+    where the entries mean:
+    \table
+      Property  | Description                         | Type | Reqd | Deflt
+      type      | Type name: HuaXu                    | word | yes  | -
+      Kburn     | Burn out factor                   | scalar | no   | 1.5
+    \endtable
+
 SourceFiles
     HuaXu.C
 
@@ -64,17 +86,27 @@ class HuaXu
 :
     public CHFSubCoolModel
 {
-
-    // Private data:
+    // Private Data
 
         //- Burn out factor
         scalar Kburn_;
 
+
+    // Private Member Functions
+
+        //- No copy construct
+        HuaXu(const HuaXu&) = delete;
+
+        //- No copy assignment
+        void operator=(const HuaXu&) = delete;
+
+
 public:
 
     //- Runtime type information
     TypeName("HuaXu");
 
+
     // Constructors
 
         //- Construct from a dictionary
@@ -82,7 +114,7 @@ public:
 
 
     //- Destructor
-    virtual ~HuaXu();
+    virtual ~HuaXu() = default;
 
 
     // Member Functions
@@ -98,8 +130,10 @@ public:
             const scalarField& L
         ) const;
 
+        // I-O
 
-        virtual void write(Ostream& os) const;
+            //- Write
+            virtual void write(Ostream& os) const;
 };
 
 

src/phaseSystemModels/reactingEuler/multiphaseSystem/derivedFvPatchFields/wallBoilingSubModels/CHFSubCoolModels/Tatsumoto/Tatsumoto.C

@@ -0,0 +1,105 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2021 OpenCFD Ltd
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "Tatsumoto.H"
+#include "addToRunTimeSelectionTable.H"
+#include "uniformDimensionedFields.H"
+#include "phasePairKey.H"
+#include "phaseSystem.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace wallBoilingModels
+{
+namespace CHFModels
+{
+    defineTypeNameAndDebug(Tatsumoto, 0);
+    addToRunTimeSelectionTable
+    (
+        CHFSubCoolModel,
+        Tatsumoto,
+        dictionary
+    );
+}
+}
+}
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::wallBoilingModels::CHFModels::Tatsumoto::Tatsumoto
+(
+    const dictionary& dict
+)
+:
+    CHFSubCoolModel(),
+    K_(dict.getOrDefault<scalar>("K", 0.23))
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+Foam::tmp<Foam::scalarField>
+Foam::wallBoilingModels::CHFModels::Tatsumoto::CHFSubCool
+(
+    const phaseModel& liquid,
+    const phaseModel& vapor,
+    const label patchi,
+    const scalarField& Tl,
+    const scalarField& Tsatw,
+    const scalarField& L
+) const
+{
+    const labelUList& cells = liquid.mesh().boundary()[patchi].faceCells();
+    const scalarField& pw = liquid.thermo().p().boundaryField()[patchi];
+
+    tmp<scalarField> trhoVapor = vapor.thermo().rhoEoS(pw, Tsatw, cells);
+    const scalarField& rhoVapor = trhoVapor.ref();
+
+    tmp<scalarField> trhoLiq = liquid.thermo().rhoEoS(pw, Tsatw, cells);
+    const scalarField& rhoLiq = trhoLiq.ref();
+
+    tmp<scalarField> tCp = liquid.thermo().Cp(pw, Tsatw, cells);
+    const scalarField& Cp = tCp();
+
+    return
+        1 + K_*pow(rhoVapor/rhoLiq, 0.8)*Cp*max(Tsatw - Tl, scalar(0))/L;
+}
+
+
+void Foam::wallBoilingModels::CHFModels::Tatsumoto::write
+(
+    Ostream& os
+) const
+{
+    CHFSubCoolModel::write(os);
+    os.writeEntry("K", K_);
+}
+
+
+// ************************************************************************* //

src/phaseSystemModels/reactingEuler/multiphaseSystem/derivedFvPatchFields/wallBoilingSubModels/CHFSubCoolModels/Tatsumoto/Tatsumoto.H

@@ -0,0 +1,138 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2021 OpenCFD Ltd
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::wallBoilingModels::CHFModels::Tatsumoto
+
+Description
+    A critical heat flux (CHF) sub-cooling correlation model.
+
+Usage
+    Example of the model specification:
+    \verbatim
+    CHFSubCoolModel
+    {
+        // Mandatory entries
+        type            Tatsumoto;
+
+        // Optional entries
+        K               <scalar>;
+    }
+    \endverbatim
+
+    where the entries mean:
+    \table
+      Property  | Description                         | Type | Reqd | Deflt
+      type      | Type name: Tatsumoto                | word | yes  | -
+      K         | Model coefficient                   | scalar | no | 0.23
+    \endtable
+
+SourceFiles
+    Tatsumoto.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef Tatsumoto_H
+#define Tatsumoto_H
+
+#include "CHFSubCoolModel.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace wallBoilingModels
+{
+namespace CHFModels
+{
+
+/*---------------------------------------------------------------------------*\
+                         Class Tatsumoto Declaration
+\*---------------------------------------------------------------------------*/
+
+class Tatsumoto
+:
+    public CHFSubCoolModel
+{
+    // Private Data
+
+        //- Model coefficient
+        scalar K_;
+
+
+    // Private Member Functions
+
+        //- No copy construct
+        Tatsumoto(const Tatsumoto&) = delete;
+
+        //- No copy assignment
+        void operator=(const Tatsumoto&) = delete;
+
+
+public:
+
+    //- Runtime type information
+    TypeName("Tatsumoto");
+
+
+    // Constructors
+
+        //- Construct from a dictionary
+        Tatsumoto(const dictionary& dict);
+
+
+    //- Destructor
+    virtual ~Tatsumoto() = default;
+
+
+    // Member Functions
+
+        //- Calculate and return the nucleation-site density
+        virtual tmp<scalarField> CHFSubCool
+        (
+            const phaseModel& liquid,
+            const phaseModel& vapor,
+            const label patchi,
+            const scalarField& Tl,
+            const scalarField& Tsatw,
+            const scalarField& L
+        ) const;
+
+        //- Write
+        virtual void write(Ostream& os) const;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace CHFModels
+} // End namespace wallBoilingModels
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/phaseSystemModels/reactingEuler/multiphaseSystem/derivedFvPatchFields/wallBoilingSubModels/LeidenfrostModels/LeidenfrostModel/LeidenfrostModel.H

@@ -5,7 +5,7 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2018-2020 OpenCFD Ltd
+    Copyright (C) 2018-2021 OpenCFD Ltd
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -27,7 +27,7 @@ Class
     Foam::wallBoilingModels::LeidenfrostModel
 
 Description
-    Base class for nucleation site density models
+    Base class for Leidenfrost-effect models.
 
 SourceFiles
     LeidenfrostModel.C
@@ -92,7 +92,7 @@ public:
 
     // Member Functions
 
-        //- Calculate temperature
+        //- Calculate and return the Leidenfrost temperature
         virtual tmp<scalarField> TLeid
         (
             const phaseModel& liquid,
@@ -103,7 +103,10 @@ public:
             const scalarField& L
         ) const = 0;
 
-        virtual void write(Ostream& os) const;
+        // I-O
+
+            //- Write
+            virtual void write(Ostream& os) const;
 };
 
 

src/phaseSystemModels/reactingEuler/multiphaseSystem/derivedFvPatchFields/wallBoilingSubModels/LeidenfrostModels/Spiegler/Spiegler.C

@@ -5,7 +5,7 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2018-2020 OpenCFD Ltd
+    Copyright (C) 2018-2021 OpenCFD Ltd
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -61,12 +61,6 @@ Foam::wallBoilingModels::LeidenfrostModels::Spiegler::Spiegler
 {}
 
 
-// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
-
-Foam::wallBoilingModels::LeidenfrostModels::Spiegler::~Spiegler()
-{}
-
-
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
 Foam::tmp<Foam::scalarField>
@@ -80,13 +74,10 @@ Foam::wallBoilingModels::LeidenfrostModels::Spiegler::TLeid
     const scalarField& L
 ) const
 {
-    return tmp<scalarField>
+    return tmp<scalarField>::New
     (
-        new scalarField
-        (
-            liquid.thermo().p().boundaryField()[patchi].size(),
-            27*Tcrit_/32
-        )
+        liquid.thermo().p().boundaryField()[patchi].size(),
+        scalar(27)*Tcrit_/scalar(32)
     );
 }
 
@@ -100,4 +91,5 @@ void Foam::wallBoilingModels::LeidenfrostModels::Spiegler::write
     os.writeEntry("Tcrit", Tcrit_);
 }
 
+
 // ************************************************************************* //

src/phaseSystemModels/reactingEuler/multiphaseSystem/derivedFvPatchFields/wallBoilingSubModels/LeidenfrostModels/Spiegler/Spiegler.H

@@ -5,7 +5,7 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2018 OpenCFD Ltd
+    Copyright (C) 2018-2021 OpenCFD Ltd
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -24,18 +24,41 @@ License
     along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
 
 Class
-    Foam::wallBoilingModels:LeidenfrostModels:::Spiegler
+    Foam::wallBoilingModels::LeidenfrostModels::Spiegler
 
 Description
-    Leidenfrost temperature model.
+    A model for Leidenfrost effects based on
+    Spiegler et al. (1963) for boiling flows.
 
-    References:
+    Reference:
     \verbatim
-        SPIEGLER P., HOPENFELD J., SILBERBERG M., BUMPUS J. and NORMAN A.,
-        Onset of stable film boiling and the foam limit, International
-        Journal of Heat and Mass Transfer, 6,11, pp.987-989, 1963
+        Spiegler, P., Hopenfeld, J., Silberberg, M.,
+        Bumpus Jr, C. F., & Norman, A. (1963).
+        Onset of stable film boiling and the foam limit.
+        International Journal of Heat and Mass Transfer, 6(11), 987-989.
+        DOI:10.1016/0017-9310(63)90053-X
     \endverbatim
 
+Usage
+    Example of the model specification:
+    \verbatim
+    LeidenfrostModel
+    {
+        // Mandatory entries
+        type            Spiegler;
+
+        // Optional entries
+        Tcrit            <scalar>;
+    }
+    \endverbatim
+
+    where the entries mean:
+    \table
+      Property  | Description                         | Type | Reqd | Deflt
+      type      | Type name: Spiegler                 | word | yes  | -
+      Tcrit     | Critical temperature [K]          | scalar | no   | 374.0
+    \endtable
+
 SourceFiles
     Spiegler.C
 
@@ -63,11 +86,19 @@ class Spiegler
 :
     public LeidenfrostModel
 {
+    // Private Data
 
-private:
+        //- Critical temperature
+        scalar Tcrit_;
 
-    //- Critical temperature
-    scalar Tcrit_;
+
+    // Private Member Functions
+
+        //- No copy construct
+        Spiegler(const Spiegler&) = delete;
+
+        //- No copy assignment
+        void operator=(const Spiegler&) = delete;
 
 
 public:
@@ -75,6 +106,7 @@ public:
     //- Runtime type information
     TypeName("Spiegler");
 
+
     // Constructors
 
         //- Construct from a dictionary
@@ -82,12 +114,12 @@ public:
 
 
     //- Destructor
-    virtual ~Spiegler();
+    virtual ~Spiegler() = default;
 
 
     // Member Functions
 
-        //- Calculate and return the nucleation-site density
+        //- Calculate and return the Leidenfrost temperature
         virtual tmp<scalarField> TLeid
         (
             const phaseModel& liquid,
@@ -98,8 +130,10 @@ public:
             const scalarField& L
         ) const;
 
+        // I-O
 
-        virtual void write(Ostream& os) const;
+            //- Write
+            virtual void write(Ostream& os) const;
 };
 
 

src/phaseSystemModels/reactingEuler/multiphaseSystem/derivedFvPatchFields/wallBoilingSubModels/MHFModels/Jeschar/Jeschar.C

@@ -5,7 +5,7 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2018-2020 OpenCFD Ltd
+    Copyright (C) 2018-2021 OpenCFD Ltd
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -59,13 +59,7 @@ Foam::wallBoilingModels::CHFModels::Jeschar::Jeschar
 )
 :
     MHFModel(),
-    Kmhf_(dict.getOrDefault<scalar>("Kmhf", 1))
-{}
-
-
-// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
-
-Foam::wallBoilingModels::CHFModels::Jeschar::~Jeschar()
+    Kmhf_(dict.getOrDefault<scalar>("Kmhf", 0.16))
 {}
 
 
@@ -85,8 +79,15 @@ Foam::wallBoilingModels::CHFModels::Jeschar::MHF
     const uniformDimensionedVectorField& g =
         liquid.mesh().time().lookupObject<uniformDimensionedVectorField>("g");
 
-    const scalarField rhoVapor(vapor.thermo().rho(patchi));
-    const scalarField rhoLiq(liquid.thermo().rho(patchi));
+    const labelUList& cells = liquid.mesh().boundary()[patchi].faceCells();
+
+    const scalarField& pw = liquid.thermo().p().boundaryField()[patchi];
+
+    tmp<scalarField> trhoVapor = vapor.thermo().rhoEoS(pw, Tsatw, cells);
+    const scalarField& rhoVapor = trhoVapor.ref();
+
+    tmp<scalarField> trhoLiq = liquid.thermo().rhoEoS(pw, Tsatw, cells);
+    const scalarField& rhoLiq = trhoLiq.ref();
 
     const phasePairKey pair(liquid.name(), vapor.name());
     const scalarField sigma
@@ -95,10 +96,13 @@ Foam::wallBoilingModels::CHFModels::Jeschar::MHF
     );
 
     return
-        Kmhf_*0.09*rhoVapor*L
+        Kmhf_*rhoVapor*L
        *(
             pow(sigma/(mag(g.value())*(rhoLiq - rhoVapor)), 0.25)
-          * sqrt(mag(g.value())*(rhoLiq - rhoVapor)/(rhoLiq + rhoVapor))
+          * sqrt
+            (
+                mag(g.value())*(rhoLiq - rhoVapor)/(rhoLiq + rhoVapor + VSMALL)
+            )
         );
 }
 

src/phaseSystemModels/reactingEuler/multiphaseSystem/derivedFvPatchFields/wallBoilingSubModels/MHFModels/Jeschar/Jeschar.H

@@ -5,7 +5,7 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2018 OpenCFD Ltd
+    Copyright (C) 2018-2021 OpenCFD Ltd
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -24,18 +24,42 @@ License
     along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
 
 Class
-    Foam::wallBoilingModels:MHFModels:::Jeschar
+    Foam::wallBoilingModels::MHFModels::Jeschar
 
 Description
-    Minimum heat flux (MHF) model.
+    A model for minimum heat flux based on
+    Jeschar et al. (1992) for boiling flows.
 
-    References:
+    Reference:
     \verbatim
-        Jeschar, E. Specht, C. Kohler, Heat Transfer during Cooling of
-        Heated Metallic Objects with Evaporating Liquids,
-        Theory and Technology in Quenching, Springer, 1992. Chapter 4.
+        Jeschar, R., Specht, E., & Köhler, C. (1992).
+        Heat transfer during cooling of heated
+        metallic objects with evaporating liquids.
+        In Theory and Technology of Quenching (pp. 73-92).
+        Springer, Berlin, Heidelberg.
+        DOI:10.1007/978-3-662-01596-4_4
     \endverbatim
 
+Usage
+    Example of the model specification:
+    \verbatim
+    MHFModel
+    {
+        // Mandatory entries
+        type            Jeschar;
+
+        // Optional entries
+        Kmhf            <scalar>;
+    }
+    \endverbatim
+
+    where the entries mean:
+    \table
+      Property  | Description                         | Type | Reqd | Deflt
+      type      | Type name: Jeschar                  | word | yes  | -
+      Kmhf      | Burn out factor                   | scalar | no   | 1.0
+    \endtable
+
 SourceFiles
     Jeschar.C
 
@@ -63,17 +87,27 @@ class Jeschar
 :
     public MHFModel
 {
-
-    // Private data:
+    // Private Data
 
         //- Burn out factor
         scalar Kmhf_;
 
+
+    // Private Member Functions
+
+        //- No copy construct
+        Jeschar(const Jeschar&) = delete;
+
+        //- No copy assignment
+        void operator=(const Jeschar&) = delete;
+
+
 public:
 
     //- Runtime type information
     TypeName("Jeschar");
 
+
     // Constructors
 
         //- Construct from a dictionary
@@ -81,12 +115,12 @@ public:
 
 
     //- Destructor
-    virtual ~Jeschar();
+    virtual ~Jeschar() = default;
 
 
     // Member Functions
 
-        //- Calculate and return the nucleation-site density
+        //- Calculate and return the minimum heat flux
         virtual tmp<scalarField> MHF
         (
             const phaseModel& liquid,
@@ -97,8 +131,10 @@ public:
             const scalarField& L
         ) const;
 
+        // I-O
 
-        virtual void write(Ostream& os) const;
+            //- Write
+            virtual void write(Ostream& os) const;
 };
 
 

src/phaseSystemModels/reactingEuler/multiphaseSystem/derivedFvPatchFields/wallBoilingSubModels/MHFModels/MHFModel/MHFModel.H

@@ -5,7 +5,7 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2018-2020 OpenCFD Ltd
+    Copyright (C) 2018-2021 OpenCFD Ltd
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -27,7 +27,7 @@ Class
     Foam::MHFModels::MHFModel
 
 Description
-    Base class for nucleation site density models
+    Base class for minimum heat flux (MHF) models.
 
 SourceFiles
     MHFModel.C
@@ -92,7 +92,7 @@ public:
 
     // Member Functions
 
-        //- Calculate temperature
+        //- Calculate and return the minimum heat flux
         virtual tmp<scalarField> MHF
         (
             const phaseModel& liquid,
@@ -103,7 +103,10 @@ public:
             const scalarField& L
         ) const = 0;
 
-        virtual void write(Ostream& os) const;
+        // I-O
+
+            //- Write
+            virtual void write(Ostream& os) const;
 };
 
 

src/phaseSystemModels/reactingEuler/multiphaseSystem/derivedFvPatchFields/wallBoilingSubModels/TDNBModels/Schroeder/Schroeder.C

@@ -5,7 +5,7 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2018-2020 OpenCFD Ltd
+    Copyright (C) 2018-2021 OpenCFD Ltd
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -64,12 +64,6 @@ Foam::wallBoilingModels::TDNBModels::Schroeder::Schroeder
 {}
 
 
-// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
-
-Foam::wallBoilingModels::TDNBModels::Schroeder::~Schroeder()
-{}
-
-
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
 Foam::tmp<Foam::scalarField>
@@ -92,7 +86,7 @@ Foam::wallBoilingModels::TDNBModels::Schroeder::TDNB
         Tsatw
         /
         (
-            1 - log(2*kg_ + 1)*(R.value()*Tsatw)/(W*L)
+            scalar(1) - log(scalar(2)*kg_ + scalar(1))*(R.value()*Tsatw)/(W*L)
         );
 }
 
@@ -106,4 +100,5 @@ void Foam::wallBoilingModels::TDNBModels::Schroeder::write
     os.writeEntry("kg", kg_);
 }
 
+
 // ************************************************************************* //

src/phaseSystemModels/reactingEuler/multiphaseSystem/derivedFvPatchFields/wallBoilingSubModels/TDNBModels/Schroeder/Schroeder.H

@@ -5,7 +5,7 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2018 OpenCFD Ltd
+    Copyright (C) 2018-2021 OpenCFD Ltd
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -24,33 +24,45 @@ License
     along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
 
 Class
-    Foam::wallBoilingModels:TDNBModels:::Schroeder
+    Foam::wallBoilingModels::TDNBModels::Schroeder
 
 Description
-    Departure from nulceate boiling correlation.
+    A model for departure from nucleation boiling based on
+    Schroeder-Richter and Bartsch (1994) for boiling flows.
 
-    References:
+    Reference:
     \verbatim
-        Schroeder-Richter D. and Bartsch G. Analytical calculation of
-        DNB-superheating by a pos-tulated thermo-mechanical effect of
-        nucleate boiling.
-        International Journal of Multiphase Flow, 20(6):1143–1167, 1994.
+        Schroeder-Richter, D., & Bartsch, G. (1994).
+        Analytical calculation of DNB-superheating by a postulated
+        thermo-mechanical effect of nucleate boiling.
+        International journal of multiphase flow, 20(6), 1143-1167.
+        DOI:10.1016/0301-9322(94)90060-4
 
+        Theler, G., & Freis, D. (2011).
+        Theoretical critical heat flux prediction based on non-equilibrium
+        thermodynamics considerations of the subcooled boiling phenomenon.
+        Mecánica Computacional, 30(19), 1713-1732.
     \endverbatim
 
-
+Usage
+    Example of the model specification:
     \verbatim
-        THEORETICAL CRITICAL HEAT FLUX PREDICTION BASEDNON-EQUILIBRIUM
-        THERMODYNAMICS CONSIDERATIONSTHE SUBCOOLED BOILING PHENOMENON
+    TDNBModel
+    {
+        // Mandatory entries
+        type            Schroeder;
 
-        Germán Thelera and Daniel Freisba TECNA
-        Estudios y Proyectos de Ingenierı́a S.A.
-        Encarnación Ezcurra 365, C1107CLA Buenos Aires, Argentina
-        Westinghouse
-        Electric Germany GmbH
-        Dudenstraße 44, 68167 Mannheim, Germany
+        // Optional entries
+        kg              <scalar>;
+    }
     \endverbatim
 
+    where the entries mean:
+    \table
+      Property   | Description                       | Type | Reqd | Deflt
+      type       | Type name: Schroeder              | word | yes  | -
+      kg  | Isoentropic expansion factor for ideal gases | scalar | no  | 1.666
+    \endtable
 
 SourceFiles
     Schroeder.C
@@ -79,19 +91,29 @@ class Schroeder
 :
     public TDNBModel
 {
-
-    // Private data:
+    // Private Data
 
         //- Isoentropic expansion factor for ideal gases
         // 5/3 monoatomic
         // 7/5 diatomic
         scalar kg_;
 
+
+    // Private Member Functions
+
+        //- No copy construct
+        Schroeder(const Schroeder&) = delete;
+
+        //- No copy assignment
+        void operator=(const Schroeder&) = delete;
+
+
 public:
 
     //- Runtime type information
     TypeName("Schroeder");
 
+
     // Constructors
 
         //- Construct from a dictionary
@@ -99,12 +121,12 @@ public:
 
 
     //- Destructor
-    virtual ~Schroeder();
+    virtual ~Schroeder() = default;
 
 
     // Member Functions
 
-        //- Calculate and return the nucleation-site density
+        //- Calculate and return the departure from nulceate boiling correlation
         virtual tmp<scalarField> TDNB
         (
             const phaseModel& liquid,
@@ -115,8 +137,10 @@ public:
             const scalarField& L
         ) const;
 
+        // I-O
 
-        virtual void write(Ostream& os) const;
+            //- Write
+            virtual void write(Ostream& os) const;
 };
 
 

src/phaseSystemModels/reactingEuler/multiphaseSystem/derivedFvPatchFields/wallBoilingSubModels/TDNBModels/Shirai/Shirai.C

@@ -0,0 +1,98 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2021 OpenCFD Ltd
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "Shirai.H"
+#include "addToRunTimeSelectionTable.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace wallBoilingModels
+{
+namespace TDNBModels
+{
+    defineTypeNameAndDebug(Shirai, 0);
+    addToRunTimeSelectionTable
+    (
+        TDNBModel,
+        Shirai,
+        dictionary
+    );
+}
+}
+}
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::wallBoilingModels::TDNBModels::Shirai::Shirai
+(
+    const dictionary& dict
+)
+:
+    TDNBModel(),
+    Tc_(dict.get<scalar>("Tc")),
+    Pc_(dict.get<scalar>("Pc"))
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+Foam::tmp<Foam::scalarField>
+Foam::wallBoilingModels::TDNBModels::Shirai::TDNB
+(
+    const phaseModel& liquid,
+    const phaseModel& vapor,
+    const label patchi,
+    const scalarField& Tl,
+    const scalarField& Tsatw,
+    const scalarField& L
+) const
+{
+    tmp<scalarField> tp = liquid.thermo().p().boundaryField()[patchi];
+
+    const scalarField pRatio(max(min(tp/Pc_, scalar(1)), scalar(0)));
+
+    return
+    (
+        (0.8823*pow3(pRatio) - 1.8938*sqr(pRatio) + 1.4322*pRatio + 0.6289)*Tc_
+    );
+}
+
+
+void Foam::wallBoilingModels::TDNBModels::Shirai::write
+(
+    Ostream& os
+) const
+{
+    TDNBModel::write(os);
+    os.writeEntry("Tc", Tc_);
+    os.writeEntry("Pc", Pc_);
+}
+
+
+// ************************************************************************* //

src/phaseSystemModels/reactingEuler/multiphaseSystem/derivedFvPatchFields/wallBoilingSubModels/TDNBModels/Shirai/Shirai.H

@@ -0,0 +1,155 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2021 OpenCFD Ltd
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::wallBoilingModels::TDNBModels::Shirai
+
+Description
+    Temperature of departure from nulceate boiling correlation.
+
+    References:
+    \verbatim
+        Shirai, Y., Tatsumoto, H., Shiotsu, M., Hata, K.,
+        Kobayashi, H., Naruo, Y., & Inatani, Y. (2010).
+        Boiling heat transfer from a horizontal flat
+        plate in a pool of liquid hydrogen.
+        Cryogenics, 50(6-7), 410-416.
+        DOI:10.1016/j.cryogenics.2010.04.001
+    \endverbatim
+
+Usage
+    Example of the model specification:
+    \verbatim
+    TDNBModel
+    {
+        // Mandatory entries
+        type            Shirai;
+        Tc              <scalar>;
+        Pc              <scalar>;
+    }
+    \endverbatim
+
+    where the entries mean:
+    \table
+      Property   | Description                       | Type | Reqd | Deflt
+      type       | Type name: Shirai                 | word | yes  | -
+      Tc         | Critical temperature            | scalar | yes  | -
+      Pc         | Critical pressure               | scalar | yes  | -
+    \endtable
+
+Note
+  - Correlation based on fiting data from Fig 11 from above references.
+  - Suitable for liquid Helium, Nitrogen, Oxygen and Hydrogen.
+
+SourceFiles
+    Shirai.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef Shirai_H
+#define Shirai_H
+
+#include "TDNBModel.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace wallBoilingModels
+{
+namespace TDNBModels
+{
+
+/*---------------------------------------------------------------------------*\
+                         Class Shirai Declaration
+\*---------------------------------------------------------------------------*/
+
+class Shirai
+:
+    public TDNBModel
+{
+    // Private Data
+
+        //- Critical temperature
+        scalar Tc_;
+
+        //- Critical pressure
+        scalar Pc_;
+
+
+    // Private Member Functions
+
+        //- No copy construct
+        Shirai(const Shirai&) = delete;
+
+        //- No copy assignment
+        void operator=(const Shirai&) = delete;
+
+
+public:
+
+    //- Runtime type information
+    TypeName("Shirai");
+
+
+    // Constructors
+
+        //- Construct from a dictionary
+        Shirai(const dictionary& dict);
+
+
+    //- Destructor
+    virtual ~Shirai() = default;
+
+
+    // Member Functions
+
+        //- Calculate and return the nucleation-site density
+        virtual tmp<scalarField> TDNB
+        (
+            const phaseModel& liquid,
+            const phaseModel& vapor,
+            const label patchi,
+            const scalarField& Tl,
+            const scalarField& Tsatw,
+            const scalarField& L
+        ) const;
+
+        //- Write
+        virtual void write(Ostream& os) const;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace TDNBModels
+} // End namespace wallBoilingModels
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/phaseSystemModels/reactingEuler/multiphaseSystem/derivedFvPatchFields/wallBoilingSubModels/TDNBModels/TDNBModel/TDNBModel.H

@@ -5,7 +5,7 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2018-2020 OpenCFD Ltd
+    Copyright (C) 2018-2021 OpenCFD Ltd
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -27,7 +27,7 @@ Class
     Foam::wallBoilingModels::TDNBModel
 
 Description
-    Base class for nucleation site density models
+    Base class for departure from nucleation boiling models.
 
 SourceFiles
     TDNBModel.C
@@ -92,7 +92,7 @@ public:
 
     // Member Functions
 
-        //- Calculate temperature
+        //- Calculate and return the departure from nulceate boiling correlation
         virtual tmp<scalarField> TDNB
         (
             const phaseModel& liquid,
@@ -103,7 +103,10 @@ public:
             const scalarField& L
         ) const = 0;
 
-        virtual void write(Ostream& os) const;
+        // I-O
+
+            //- Write
+            virtual void write(Ostream& os) const;
 };
 
 

src/phaseSystemModels/reactingEuler/multiphaseSystem/derivedFvPatchFields/wallBoilingSubModels/departureDiameterModels/KocamustafaogullariIshii/KocamustafaogullariIshii.C

@@ -67,13 +67,6 @@ KocamustafaogullariIshii::KocamustafaogullariIshii
 {}
 
 
-// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
-
-Foam::wallBoilingModels::departureDiameterModels::
-KocamustafaogullariIshii::~KocamustafaogullariIshii()
-{}
-
-
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
 Foam::tmp<Foam::scalarField>
@@ -89,7 +82,7 @@ KocamustafaogullariIshii::dDeparture
 ) const
 {
     // Gravitational acceleration
-    const uniformDimensionedVectorField& g =
+    const auto& g =
         liquid.mesh().time().lookupObject<uniformDimensionedVectorField>("g");
 
     const scalarField rhoLiquid(liquid.thermo().rho(patchi));
@@ -97,10 +90,9 @@ KocamustafaogullariIshii::dDeparture
 
     const scalarField rhoM((rhoLiquid - rhoVapor)/rhoVapor);
 
-    const tmp<volScalarField>& tsigma
-    (
-        liquid.fluid().sigma(phasePairKey(liquid.name(), vapor.name()))
-    );
+    const tmp<volScalarField>& tsigma =
+        liquid.fluid().sigma(phasePairKey(liquid.name(), vapor.name()));
+
     const volScalarField& sigma = tsigma();
     const fvPatchScalarField& sigmaw = sigma.boundaryField()[patchi];
 

src/phaseSystemModels/reactingEuler/multiphaseSystem/derivedFvPatchFields/wallBoilingSubModels/departureDiameterModels/KocamustafaogullariIshii/KocamustafaogullariIshii.H

@@ -6,6 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2016-2018 OpenFOAM Foundation
+    Copyright (C) 2021 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -27,17 +28,35 @@ Class
     Foam::wallBoilingModels::departureDiameterModels::KocamustafaogullariIshii
 
 Description
-    A correlation for bubble departure diameter.
-
-    Requires model parameter 'phi': contact angle in degrees.
+    A correlation for bubble departure diameter modelling
+    based on Kocamustafaogullari-Ishii (1983) for boiling flows.
 
     Reference:
     \verbatim
         Kocamustafaogullari, G., & Ishii, M. (1983).
         Interfacial area and nucleation site density in boiling systems.
         International Journal of Heat and Mass Transfer, 26(9), 1377-1387.
+        DOI:10.1016/S0017-9310(83)80069-6
     \endverbatim
 
+Usage
+    Example of the model specification:
+    \verbatim
+    departureDiameterModel
+    {
+        // Mandatory entries
+        type            KocamustafaogullariIshii;
+        phi             <scalar>;
+    }
+    \endverbatim
+
+    where the entries mean:
+    \table
+      Property  | Description                         | Type | Reqd | Deflt
+      type      | Type name: KocamustafaogullariIshii | word | yes  | -
+      phi       | Contact angle [deg]                 | scalar | yes  | -
+    \endtable
+
 SourceFiles
     KocamustafaogullariIshii.C
 
@@ -58,19 +77,28 @@ namespace departureDiameterModels
 {
 
 /*---------------------------------------------------------------------------*\
-                         Class KocamustafaogullariIshii Declaration
+                  Class KocamustafaogullariIshii Declaration
 \*---------------------------------------------------------------------------*/
 
 class KocamustafaogullariIshii
 :
     public departureDiameterModel
 {
-    // Private data
+    // Private Data
 
         //- Contact angle
         scalar phi_;
 
 
+    // Private Member Functions
+
+        //- No copy construct
+        KocamustafaogullariIshii(const KocamustafaogullariIshii&) = delete;
+
+        //- No copy assignment
+        void operator=(const KocamustafaogullariIshii&) = delete;
+
+
 public:
 
     //- Runtime type information
@@ -83,7 +111,7 @@ public:
 
 
     //- Destructor
-    virtual ~KocamustafaogullariIshii();
+    virtual ~KocamustafaogullariIshii() = default;
 
 
     // Member Functions
@@ -99,7 +127,10 @@ public:
             const scalarField& L
         ) const;
 
-        virtual void write(Ostream& os) const;
+        // I-O
+
+            //- Write
+            virtual void write(Ostream& os) const;
 };
 
 

src/phaseSystemModels/reactingEuler/multiphaseSystem/derivedFvPatchFields/wallBoilingSubModels/departureDiameterModels/TolubinskiKostanchuk/TolubinskiKostanchuk.C

@@ -64,13 +64,6 @@ TolubinskiKostanchuk::TolubinskiKostanchuk
 {}
 
 
-// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
-
-Foam::wallBoilingModels::departureDiameterModels::
-TolubinskiKostanchuk::~TolubinskiKostanchuk()
-{}
-
-
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
 Foam::tmp<Foam::scalarField>
@@ -85,7 +78,7 @@ TolubinskiKostanchuk::dDeparture
     const scalarField& L
 ) const
 {
-    return max(min(dRef_*exp(-(Tsatw - Tl)/45), dMax_), dMin_);
+    return max(min(dRef_*exp(-(Tsatw - Tl)/scalar(45)), dMax_), dMin_);
 }
 
 

src/phaseSystemModels/reactingEuler/multiphaseSystem/derivedFvPatchFields/wallBoilingSubModels/departureDiameterModels/TolubinskiKostanchuk/TolubinskiKostanchuk.H

@@ -6,6 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2016-2018 OpenFOAM Foundation
+    Copyright (C) 2021 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -27,16 +28,41 @@ Class
     Foam::wallBoilingModels::departureDiameterModels::TolubinskiKostanchuk
 
 Description
-    Tolubinski-Kostanchuk correlation for bubble departure diameter.
+    A correlation for bubble departure diameter modelling
+    based on Tolubinski-Kostanchuk (1970) for boiling flows.
 
     Reference:
     \verbatim
         Tolubinsky, V. I., & Kostanchuk, D. M. (1970).
-        Vapour bubbles growth rate and heat transfer intensity at subcooled
-        water boiling.
-        In International Heat Transfer Conference 4 (Vol. 23). Begel House Inc.
+        Vapour bubbles growth rate and heat transfer
+        intensity at subcooled water boiling.
+        In International Heat Transfer Conference 4 (Vol. 23). Begel House Inc..
     \endverbatim
 
+Usage
+    Example of the model specification:
+    \verbatim
+    departureDiameterModel
+    {
+        // Mandatory entries
+        type            TolubinskiKostanchuk;
+
+        // Optional entries
+        dRef            <scalar>;
+        dMax            <scalar>;
+        dMin            <scalar>;
+    }
+    \endverbatim
+
+    where the entries mean:
+    \table
+      Property  | Description                         | Type | Reqd | Deflt
+      type      | Type name: TolubinskiKostanchuk     | word | yes  | -
+      dRef      | Coefficient of the temperature term | scalar | no | 6e-4
+      dMax      | Maximum diameter                    | scalar | no | 0.0014
+      dMin      | Minimum diameter                    | scalar | no | 1e-6
+    \endtable
+
 SourceFiles
     TolubinskiKostanchuk.C
 
@@ -64,8 +90,7 @@ class TolubinskiKostanchuk
 :
     public departureDiameterModel
 {
-
-    // Private data:
+    // Private Data
 
         //- Coefficient of the temperature term
         scalar dRef_;
@@ -76,6 +101,16 @@ class TolubinskiKostanchuk
         //- Minimum diameter
         scalar dMin_;
 
+
+    // Private Member Functions
+
+        //- No copy construct
+        TolubinskiKostanchuk(const TolubinskiKostanchuk&) = delete;
+
+        //- No copy assignment
+        void operator=(const TolubinskiKostanchuk&) = delete;
+
+
 public:
 
     //- Runtime type information
@@ -89,7 +124,7 @@ public:
 
 
     //- Destructor
-    virtual ~TolubinskiKostanchuk();
+    virtual ~TolubinskiKostanchuk() = default;
 
 
     // Member Functions
@@ -105,7 +140,10 @@ public:
             const scalarField& L
         ) const;
 
-        virtual void write(Ostream& os) const;
+        // I-O
+
+            // Write
+            virtual void write(Ostream& os) const;
 };
 
 

src/phaseSystemModels/reactingEuler/multiphaseSystem/derivedFvPatchFields/wallBoilingSubModels/departureDiameterModels/departureDiameterModel/departureDiameterModel.C

@@ -76,4 +76,5 @@ void Foam::wallBoilingModels::departureDiameterModel::write(Ostream& os) const
     os.writeEntry("type", this->type());
 }
 
+
 // ************************************************************************* //

src/phaseSystemModels/reactingEuler/multiphaseSystem/derivedFvPatchFields/wallBoilingSubModels/departureDiameterModels/departureDiameterModel/departureDiameterModel.H

@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2016-2018 OpenFOAM Foundation
-    Copyright (C) 2020 OpenCFD Ltd.
+    Copyright (C) 2020-2021 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -28,7 +28,7 @@ Class
     Foam::wallBoilingModels::departureDiameterModel
 
 Description
-    Base class for bubble departure diameter models
+    Base class for bubble departure diameter models for boiling flows.
 
 SourceFiles
     departureDiameterModel.C
@@ -104,7 +104,10 @@ public:
             const scalarField& L
         ) const = 0;
 
-        virtual void write(Ostream& os) const;
+        // I-O
+
+            //- Write
+            virtual void write(Ostream& os) const;
 };
 
 

src/phaseSystemModels/reactingEuler/multiphaseSystem/derivedFvPatchFields/wallBoilingSubModels/departureFrequencyModels/Cole/Cole.C

@@ -62,13 +62,6 @@ Cole::Cole(const dictionary& dict)
 {}
 
 
-// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
-
-Foam::wallBoilingModels::departureFrequencyModels::
-Cole::~Cole()
-{}
-
-
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
 Foam::tmp<Foam::scalarField>
@@ -82,7 +75,7 @@ Cole::fDeparture
 ) const
 {
     // Gravitational acceleration
-    const uniformDimensionedVectorField& g =
+    const auto& g =
         liquid.mesh().time().lookupObject<uniformDimensionedVectorField>("g");
 
     const scalarField rhoLiquid(liquid.thermo().rho(patchi));
@@ -90,9 +83,9 @@ Cole::fDeparture
 
     return sqrt
     (
-        4*mag(g).value()
+        scalar(4)*mag(g).value()
        *max(rhoLiquid - rhoVapor, scalar(0.1))
-       /(3*dDep*rhoLiquid)
+       /(scalar(3)*dDep*rhoLiquid)
     );
 }
 

src/phaseSystemModels/reactingEuler/multiphaseSystem/derivedFvPatchFields/wallBoilingSubModels/departureFrequencyModels/Cole/Cole.H

@@ -6,6 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2016-2018 OpenFOAM Foundation
+    Copyright (C) 2021 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -27,16 +28,34 @@ Class
     Foam::wallBoilingModels::departureFrequencyModels::Cole
 
 Description
-    Cole correlation for bubble departure frequency.
+    A correlation for bubble departure frequency modelling
+    based on Cole (1960) for boiling flows.
 
     Reference:
     \verbatim
         Cole, R. (1960).
-        A photographic study of pool boiling in the region of the critical heat
-        flux.
+        A photographic study of pool boiling
+        in the region of the critical heat flux.
         AIChE Journal, 6(4), 533-538.
+        DOI:10.1002/aic.690060405
     \endverbatim
 
+Usage
+    Example of the model specification:
+    \verbatim
+    departureFrequencyModel
+    {
+        // Mandatory entries
+        type            Cole;
+    }
+    \endverbatim
+
+    where the entries mean:
+    \table
+      Property  | Description                         | Type | Reqd | Deflt
+      type      | Type name: Cole                     | word | yes  | -
+    \endtable
+
 SourceFiles
     Cole.C
 
@@ -64,6 +83,14 @@ class Cole
 :
     public departureFrequencyModel
 {
+    // Private Member Functions
+
+        //- No copy construct
+        Cole(const Cole&) = delete;
+
+        //- No copy assignment
+        void operator=(const Cole&) = delete;
+
 
 public:
 
@@ -78,7 +105,7 @@ public:
 
 
     //- Destructor
-    virtual ~Cole();
+    virtual ~Cole() = default;
 
 
     // Member Functions

src/phaseSystemModels/reactingEuler/multiphaseSystem/derivedFvPatchFields/wallBoilingSubModels/departureFrequencyModels/departureFrequencyModel/departureFrequencyModel.H

@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2016-2018 OpenFOAM Foundation
-    Copyright (C) 2020 OpenCFD Ltd.
+    Copyright (C) 2020-2021 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -28,7 +28,7 @@ Class
     Foam::wallBoilingModels::departureFrequencyModel
 
 Description
-    Base class for bubble departure frequency models
+    Base class for bubble departure frequency models for boiling flows.
 
 SourceFiles
     departureFrequencyModel.C
@@ -102,7 +102,10 @@ public:
             const scalarField& dDep
         ) const = 0;
 
-        virtual void write(Ostream& os) const;
+        // I-O
+
+            //- Write
+            virtual void write(Ostream& os) const;
 };
 
 

src/phaseSystemModels/reactingEuler/multiphaseSystem/derivedFvPatchFields/wallBoilingSubModels/filmBoilingModels/BreenWestwater/BreenWestwater.C

@@ -0,0 +1,136 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2021 OpenCFD Ltd
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "BreenWestwater.H"
+#include "addToRunTimeSelectionTable.H"
+#include "uniformDimensionedFields.H"
+#include "phasePairKey.H"
+#include "phaseSystem.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace wallBoilingModels
+{
+namespace filmBoilingModels
+{
+    defineTypeNameAndDebug(BreenWestwater, 0);
+    addToRunTimeSelectionTable
+    (
+        filmBoilingModel,
+        BreenWestwater,
+        dictionary
+    );
+}
+}
+}
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::wallBoilingModels::filmBoilingModels::BreenWestwater::BreenWestwater
+(
+    const dictionary& dict
+)
+:
+    filmBoilingModel(),
+    Cn_(dict.getOrDefault<scalar>("Cn", 0.37))
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+Foam::tmp<Foam::scalarField>
+Foam::wallBoilingModels::filmBoilingModels::BreenWestwater::htcFilmBoil
+(
+    const phaseModel& liquid,
+    const phaseModel& vapor,
+    const label patchi,
+    const scalarField& Tl,
+    const scalarField& Tsatw,
+    const scalarField& L
+) const
+{
+    const fvPatchScalarField& Tw =
+        liquid.thermo().T().boundaryField()[patchi];
+
+    const auto& g =
+        liquid.mesh().time().lookupObject<uniformDimensionedVectorField>("g");
+
+    const labelUList& cells = liquid.mesh().boundary()[patchi].faceCells();
+    const scalarField& pw = liquid.thermo().p().boundaryField()[patchi];
+
+    tmp<scalarField> trhoVapor = vapor.thermo().rhoEoS(pw, Tsatw, cells);
+    const scalarField& rhoVapor = trhoVapor.ref();
+
+    tmp<scalarField> trhoLiq = liquid.thermo().rhoEoS(pw, Tsatw, cells);
+    const scalarField& rhoLiq = trhoLiq.ref();
+
+
+    const scalarField kappaLiquid(liquid.kappa(patchi));
+
+    tmp<scalarField> tCp = vapor.thermo().Cp(pw, Tsatw, cells);
+    const scalarField& CpVapor = tCp();
+
+    const scalarField nuLiquid(liquid.nu(patchi));
+
+    const scalarField Leff
+    (
+        L*sqr(1 + 0.34*CpVapor*max((Tw-Tsatw), scalar(0))/L)
+    );
+
+    const scalarField rhoDiff(rhoLiq - rhoVapor);
+
+    const phasePairKey pair(liquid.name(), vapor.name());
+    const scalarField sigma
+    (
+        liquid.fluid().sigma(pair)().boundaryField()[patchi]
+    );
+
+    return
+         Cn_
+        /pow(sigma/mag(g.value())/rhoDiff, 1/8)
+        /pow
+         (
+             nuLiquid*max((Tw-Tsatw), scalar(1e-3))
+            /(pow3(kappaLiquid)*rhoVapor*Leff*mag(g.value())*rhoDiff)
+            , 0.25
+         );
+}
+
+
+void Foam::wallBoilingModels::filmBoilingModels::BreenWestwater::write
+(
+    Ostream& os
+) const
+{
+    filmBoilingModel::write(os);
+    os.writeEntry("Cn", Cn_);
+}
+
+
+// ************************************************************************* //

src/phaseSystemModels/reactingEuler/multiphaseSystem/derivedFvPatchFields/wallBoilingSubModels/filmBoilingModels/BreenWestwater/BreenWestwater.H

@@ -0,0 +1,151 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2021 OpenCFD Ltd
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::wallBoilingModels::filmBoilingModels::BreenWestwater
+
+Description
+    Boiling film correlation.
+    A correlation for boiling film modelling
+    based on Breen & Westwater (1965) for boiling flows.
+
+    References:
+    \verbatim
+        Breen B.P. & Westwater J. W. (1965)
+        Effect of diameter of horizontal
+        tubes on film boiling heat tranfer.
+        Chem. Eng. Progr. 58. No 7.
+    \endverbatim
+
+Usage
+    Example of the model specification:
+    \verbatim
+    filmBoilingModel
+    {
+        // Mandatory entries
+        type            BreenWestwater;
+
+        // Optional entries
+        Cn              <scalar>;
+        an              <scalar>;
+        bn              <scalar>;
+        n               <scalar>;
+    }
+    \endverbatim
+
+    where the entries mean:
+    \table
+      Property   | Description                       | Type | Reqd | Deflt
+      type       | Type name: BreenWestwater         | word | yes  | -
+      Cn         | Model coefficient               | scalar | no   | 0.37
+    \endtable
+
+SourceFiles
+    BreenWestwater.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef BreenWestwater_H
+#define BreenWestwater_H
+
+#include "filmBoilingModel.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace wallBoilingModels
+{
+namespace filmBoilingModels
+{
+
+/*---------------------------------------------------------------------------*\
+                         Class BreenWestwater Declaration
+\*---------------------------------------------------------------------------*/
+
+class BreenWestwater
+:
+    public filmBoilingModel
+{
+    // Private Data
+
+        //- Model coefficient
+        scalar Cn_;
+
+
+    // Private Member Functions
+
+        //- No copy construct
+        BreenWestwater(const BreenWestwater&) = delete;
+
+        //- No copy assignment
+        void operator=(const BreenWestwater&) = delete;
+
+
+public:
+
+    //- Runtime type information
+    TypeName("BreenWestwater");
+
+
+    // Constructors
+
+        //- Construct from a dictionary
+        BreenWestwater(const dictionary& dict);
+
+
+    //- Destructor
+    virtual ~BreenWestwater() = default;
+
+
+    // Member Functions
+
+        //- Calculate and return the nucleation-site density
+        virtual tmp<scalarField> htcFilmBoil
+        (
+            const phaseModel& liquid,
+            const phaseModel& vapor,
+            const label patchi,
+            const scalarField& Tl,
+            const scalarField& Tsatw,
+            const scalarField& L
+        ) const;
+
+        //- Write
+        virtual void write(Ostream& os) const;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace filmBoilingModels
+} // End namespace wallBoilingModels
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/phaseSystemModels/reactingEuler/multiphaseSystem/derivedFvPatchFields/wallBoilingSubModels/filmBoilingModels/Bromley/Bromley.C

@@ -5,7 +5,7 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2018-2020 OpenCFD Ltd
+    Copyright (C) 2018-2021 OpenCFD Ltd
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -66,12 +66,6 @@ Foam::wallBoilingModels::filmBoilingModels::Bromley::Bromley
 {}
 
 
-// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
-
-Foam::wallBoilingModels::filmBoilingModels::Bromley::~Bromley()
-{}
-
-
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
 Foam::tmp<Foam::scalarField>
@@ -85,26 +79,28 @@ Foam::wallBoilingModels::filmBoilingModels::Bromley::htcFilmBoil
     const scalarField& L
 ) const
 {
-
     const fvPatchScalarField& Tw =
         liquid.thermo().T().boundaryField()[patchi];
-    const uniformDimensionedVectorField& g =
+    const auto& g =
         liquid.mesh().time().lookupObject<uniformDimensionedVectorField>("g");
 
-    const fvPatchScalarField& rhoVaporw
-    (
-        vapor.thermo().rho()().boundaryField()[patchi]
-    );
+    const labelUList& cells = liquid.mesh().boundary()[patchi].faceCells();
+
+    const scalarField& pw = liquid.thermo().p().boundaryField()[patchi];
+
+    tmp<scalarField> trhoVapor = vapor.thermo().rhoEoS(pw, Tsatw, cells);
+    const scalarField& rhoVapor = trhoVapor.ref();
+
+    tmp<scalarField> trhoLiq = liquid.thermo().rhoEoS(pw, Tsatw, cells);
+    const scalarField& rhoLiq = trhoLiq.ref();
+
 
-    const scalarField rhoLiq(liquid.thermo().rho(patchi));
     const scalarField kappaVapor(vapor.kappa(patchi));
 
-    tmp<volScalarField> tCp = vapor.thermo().Cp();
-    const volScalarField& Cp = tCp();
-    const scalarField& CpVapor = Cp.boundaryField()[patchi];
+    tmp<scalarField> tCp = vapor.thermo().Cp(pw, Tsatw, cells);
+    const scalarField& CpVapor = tCp();
 
     const scalarField muVapor(vapor.mu(patchi));
-    //const scalarField dbVapor(vapor.d()().boundaryField()[patchi]);
 
     const scalarField htcRad
     (
@@ -113,14 +109,13 @@ Foam::wallBoilingModels::filmBoilingModels::Bromley::htcFilmBoil
     );
 
     return
-        Cn_*pow
+        Cn_*pow025
         (
             pow3(kappaVapor)
-           *rhoVaporw*(rhoLiq - rhoVaporw)*mag(g.value())
-           *(L + 0.4*CpVapor*max((Tw-Tsatw), scalar(0)))
-           /(L_*muVapor*max((Tw-Tsatw), scalar(1e-4))),
-            0.25
-        ) + 0.75*htcRad;
+           *rhoVapor*(rhoLiq - rhoVapor)*mag(g.value())
+           *(L + scalar(0.4)*CpVapor*max((Tw-Tsatw), scalar(0)))
+           /(L_*muVapor*max((Tw-Tsatw), scalar(1e-4)))
+        ) + scalar(0.75)*htcRad;
 }
 
 

src/phaseSystemModels/reactingEuler/multiphaseSystem/derivedFvPatchFields/wallBoilingSubModels/filmBoilingModels/Bromley/Bromley.H

@@ -5,7 +5,7 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2018 OpenCFD Ltd
+    Copyright (C) 2018-2021 OpenCFD Ltd
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -24,17 +24,43 @@ License
     along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
 
 Class
-    Foam::wallBoilingModels:filmBoilingModels:::Bromley
+    Foam::wallBoilingModels::filmBoilingModels::Bromley
 
 Description
-    Boiling film correlation.
+    A correlation for boiling film modelling
+    based on Bromley (1950) for boiling flows.
 
-    References:
+    Reference:
     \verbatim
-        A. Bromley, Heat transfer in stable film boiling,
-        Chem. Eng. Prog. 58 (1950) 67–72.
+        Bromley, L. A. (1950).
+        Heat transfer in stable film boiling.
+        Chemical Engineering Progress, 46, 221-227.
     \endverbatim
 
+Usage
+    Example of the model specification:
+    \verbatim
+    filmBoilingModel
+    {
+        // Mandatory entries
+        type            Bromley;
+        L               <scalar>;
+
+        // Optional entries
+        Cn              <scalar>;
+        emissivity      <scalar>;
+    }
+    \endverbatim
+
+    where the entries mean:
+    \table
+      Property  | Description                         | Type | Reqd | Deflt
+      type      | Type name: Bromley                  | word | yes  | -
+      L         | Characteristic length scale       | scalar | yes  | -
+      Cn  | Coefficient for nucleation site density | scalar | no   | 0.62
+      emissivity | Wall emissivity                  | scalar | no   | 1.0
+    \endtable
+
 SourceFiles
     Bromley.C
 
@@ -62,8 +88,7 @@ class Bromley
 :
     public filmBoilingModel
 {
-
-    // Private data:
+    // Private Data
 
         //- Coefficient for nucleation site density
         scalar Cn_;
@@ -71,9 +96,19 @@ class Bromley
         //- Wall emissivity
         scalar emissivity_;
 
-        //- Characteristic lenght-scale
+        //- Characteristic length scale
         scalar L_;
 
+
+    // Private Member Functions
+
+        //- No copy construct
+        Bromley(const Bromley&) = delete;
+
+        //- No copy assignment
+        void operator=(const Bromley&) = delete;
+
+
 public:
 
     //- Runtime type information
@@ -86,12 +121,12 @@ public:
 
 
     //- Destructor
-    virtual ~Bromley();
+    virtual ~Bromley() = default;
 
 
     // Member Functions
 
-        //- Calculate and return the nucleation-site density
+        //- Calculate and return the film boiling correlation
         virtual tmp<scalarField> htcFilmBoil
         (
             const phaseModel& liquid,
@@ -102,8 +137,10 @@ public:
             const scalarField& L
         ) const;
 
+        // I-O
 
-        virtual void write(Ostream& os) const;
+            //- Write
+            virtual void write(Ostream& os) const;
 };
 
 

src/phaseSystemModels/reactingEuler/multiphaseSystem/derivedFvPatchFields/wallBoilingSubModels/filmBoilingModels/filmBoilingModel/filmBoilingModel.H

@@ -5,7 +5,7 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2018-2020 OpenCFD Ltd
+    Copyright (C) 2018-2021 OpenCFD Ltd
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -27,7 +27,7 @@ Class
     Foam::filmBoilingModels::filmBoilingModel
 
 Description
-    Base class for nucleation site density models
+    Base class for film boiling models.
 
 SourceFiles
     filmBoilingModel.C
@@ -92,7 +92,7 @@ public:
 
     // Member Functions
 
-        //- Calculate temperature
+        //- Calculate and return the film boiling correlation
         virtual tmp<scalarField> htcFilmBoil
         (
             const phaseModel& liquid,
@@ -103,7 +103,10 @@ public:
             const scalarField& L
         ) const = 0;
 
-        virtual void write(Ostream& os) const;
+        // I-O
+
+            //- Write
+            virtual void write(Ostream& os) const;
 };
 
 

src/phaseSystemModels/reactingEuler/multiphaseSystem/derivedFvPatchFields/wallBoilingSubModels/nucleateFluxModels/Kutadeladze/Kutadeladze.C

@@ -0,0 +1,134 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2021 OpenCFD Ltd
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "Kutadeladze.H"
+#include "addToRunTimeSelectionTable.H"
+#include "uniformDimensionedFields.H"
+#include "phasePairKey.H"
+#include "phaseSystem.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace wallBoilingModels
+{
+namespace nucleateFluxModels
+{
+    defineTypeNameAndDebug(Kutadeladze, 0);
+    addToRunTimeSelectionTable
+    (
+        nucleateFluxModel,
+        Kutadeladze,
+        dictionary
+    );
+}
+}
+}
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::wallBoilingModels::nucleateFluxModels::Kutadeladze::Kutadeladze
+(
+    const dictionary& dict
+)
+:
+    nucleateFluxModel(),
+    Cn_(dict.getOrDefault<scalar>("Cn", 5.66e-10)),
+    an_(dict.getOrDefault<scalar>("an", 2.5)),
+    bn_(dict.getOrDefault<scalar>("bn", 1)),
+    n_(dict.getOrDefault<scalar>("n", 1))
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+Foam::tmp<Foam::scalarField>
+Foam::wallBoilingModels::nucleateFluxModels::Kutadeladze::qNucleate
+(
+    const phaseModel& liquid,
+    const phaseModel& vapor,
+    const label patchi,
+    const scalarField& Tl,
+    const scalarField& Tsatw,
+    const scalarField& L
+) const
+{
+    const auto& p = liquid.mesh().lookupObject<volScalarField>("p");
+    const scalarField& pb = p.boundaryField()[patchi];
+
+    const labelUList& cells = liquid.mesh().boundary()[patchi].faceCells();
+
+    tmp<scalarField> trhoVapor = vapor.thermo().rhoEoS(pb, Tsatw, cells);
+    const scalarField& rhoVapor = trhoVapor.ref();
+
+    tmp<scalarField> trhoLiq = liquid.thermo().rhoEoS(pb, Tsatw, cells);
+    const scalarField& rhoLiq = trhoLiq.ref();
+
+    const phasePairKey pair(liquid.name(), vapor.name());
+    const scalarField sigma
+    (
+        liquid.fluid().sigma(pair)().boundaryField()[patchi]
+    );
+
+    const fvPatchScalarField& Tw =
+        liquid.thermo().T().boundaryField()[patchi];
+
+    const scalarField kappaLiquid(liquid.kappa(patchi));
+
+    tmp<scalarField> tCpliq = liquid.thermo().Cp(pb, Tsatw, cells);
+    const scalarField& Cpliquid = tCpliq();
+
+    const scalarField muLiquid(liquid.mu(patchi));
+
+    const scalarField deltaTsub
+    (
+        pow(max((Tw-Tsatw), scalar(0)), an_)
+    );
+
+    return
+        Cn_*kappaLiquid*pow(Cpliquid, 1.5)*pow(rhoLiq,1.28)*pow(pb,1.75)
+       *deltaTsub
+       /
+       (pow(muLiquid, 0.625)*pow(sigma,0.9)*pow(L,1.5)*pow(rhoVapor,1.5));
+}
+
+
+void Foam::wallBoilingModels::nucleateFluxModels::Kutadeladze::write
+(
+    Ostream& os
+) const
+{
+    nucleateFluxModel::write(os);
+    os.writeEntry("Cn", Cn_);
+    os.writeEntry("an", an_);
+    os.writeEntry("bn", bn_);
+    os.writeEntry("n", n_);
+}
+
+
+// ************************************************************************* //

src/phaseSystemModels/reactingEuler/multiphaseSystem/derivedFvPatchFields/wallBoilingSubModels/nucleateFluxModels/Kutadeladze/Kutadeladze.H

@@ -0,0 +1,160 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2021 OpenCFD Ltd
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::wallBoilingModels::nucleateFluxModels::Kutadeladze
+
+Description
+    Nucleate flux sub-cooling correlation
+
+    References:
+    \verbatim
+        Wang, L., Li, Y., Zhang, F., Xie, F., & Ma, Y. (2016).
+        Correlations for calculating heat transfer of hydrogen pool boiling.
+        International Journal of Hydrogen Energy, 41(38), 17118-17131.
+    \endverbatim
+
+Usage
+    Example of the model specification:
+    \verbatim
+    nucleateFluxModel
+    {
+        // Mandatory entries
+        type            Kutadeladze;
+
+        // Optional entries
+        Cn              <scalar>;
+        an              <scalar>;
+        bn              <scalar>;
+        n               <scalar>;
+    }
+    \endverbatim
+
+    where the entries mean:
+    \table
+      Property   | Description                       | Type | Reqd | Deflt
+      type       | Type name: Kutadeladze            | word | yes  | -
+      Cn         | Model coefficient               | scalar | no   | 5.66e-10
+      an      | Coefficient for deltaT sub-cooling | scalar | no   | 2.5
+      bn      | Exponent for n                     | scalar | no   | 1
+      n       | Nucleating boiling paramemeter     | scalar | no   | 1
+    \endtable
+
+SourceFiles
+    Kutadeladze.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef Kutadeladze_H
+#define Kutadeladze_H
+
+#include "nucleateFluxModel.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace wallBoilingModels
+{
+namespace nucleateFluxModels
+{
+
+/*---------------------------------------------------------------------------*\
+                         Class Kutadeladze Declaration
+\*---------------------------------------------------------------------------*/
+
+class Kutadeladze
+:
+    public nucleateFluxModel
+{
+    // Private Data
+
+        //- Model coefficient
+        scalar Cn_;
+
+        //- Coefficient for deltaT sub-cooling
+        scalar an_;
+
+        //- Exponent for n
+        scalar bn_;
+
+        //- Nucleating boiling paramemeter
+        scalar n_;
+
+
+    // Private Member Functions
+
+        //- No copy construct
+        Kutadeladze(const Kutadeladze&) = delete;
+
+        //- No copy assignment
+        void operator=(const Kutadeladze&) = delete;
+
+
+public:
+
+    //- Runtime type information
+    TypeName("Kutadeladze");
+
+
+    // Constructors
+
+        //- Construct from a dictionary
+        Kutadeladze(const dictionary& dict);
+
+
+    //- Destructor
+    virtual ~Kutadeladze() = default;
+
+
+    // Member Functions
+
+        //- Calculate and return the nucleation-site density
+        virtual tmp<scalarField> qNucleate
+        (
+            const phaseModel& liquid,
+            const phaseModel& vapor,
+            const label patchi,
+            const scalarField& Tl,
+            const scalarField& Tsatw,
+            const scalarField& L
+        ) const;
+
+        //- Write
+        virtual void write(Ostream& os) const;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace nucleateFluxModels
+} // End namespace wallBoilingModels
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/phaseSystemModels/reactingEuler/multiphaseSystem/derivedFvPatchFields/wallBoilingSubModels/nucleateFluxModels/exponential/exponential.C

@@ -0,0 +1,97 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2021 OpenCFD Ltd
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "exponential.H"
+#include "addToRunTimeSelectionTable.H"
+#include "uniformDimensionedFields.H"
+#include "phasePairKey.H"
+#include "phaseSystem.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace wallBoilingModels
+{
+namespace nucleateFluxModels
+{
+    defineTypeNameAndDebug(exponential, 0);
+    addToRunTimeSelectionTable
+    (
+        nucleateFluxModel,
+        exponential,
+        dictionary
+    );
+}
+}
+}
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::wallBoilingModels::nucleateFluxModels::exponential::exponential
+(
+    const dictionary& dict
+)
+:
+    nucleateFluxModel(),
+    a_(dict.getOrDefault<scalar>("a", 6309)),
+    b_(dict.getOrDefault<scalar>("b", 2.52))
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+Foam::tmp<Foam::scalarField>
+Foam::wallBoilingModels::nucleateFluxModels::exponential::qNucleate
+(
+    const phaseModel& liquid,
+    const phaseModel& vapor,
+    const label patchi,
+    const scalarField& Tl,
+    const scalarField& Tsatw,
+    const scalarField& L
+) const
+{
+    const fvPatchScalarField& Tw =
+        liquid.thermo().T().boundaryField()[patchi];
+
+    return a_*pow(max((Tw-Tsatw), scalar(0)), b_);
+}
+
+
+void Foam::wallBoilingModels::nucleateFluxModels::exponential::write
+(
+    Ostream& os
+) const
+{
+    nucleateFluxModel::write(os);
+    os.writeEntry("a", a_);
+    os.writeEntry("b", b_);
+}
+
+
+// ************************************************************************* //

src/phaseSystemModels/reactingEuler/multiphaseSystem/derivedFvPatchFields/wallBoilingSubModels/nucleateFluxModels/exponential/exponential.H

@@ -0,0 +1,158 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2021 OpenCFD Ltd
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::wallBoilingModels::nucleateFluxModels::exponential
+
+Description
+    Nucleate flux sub-cooling correlation
+
+    References:
+    \verbatim
+        Frost W. & Dzakowic G. S. (1967)
+        An extension of the methods of predicting incipient
+        boiling on commercially finished surfaces.
+        ASME AIChE Heat Transfer Conf 67-HT-61, Seattle.
+
+        Shirai, Y., Tatsumoto, H., Shiotsu, M., Hata, K.,
+        Kobayashi, H., Naruo, Y., & Inatani, Y. (2010).
+        Boiling heat transfer from a horizontal flat
+        plate in a pool of liquid hydrogen.
+        Cryogenics, 50(6-7), 410-416.
+        DOI:10.1016/j.cryogenics.2010.04.001
+    \endverbatim
+
+Usage
+    Example of the model specification:
+    \verbatim
+    nucleateFluxModel
+    {
+        // Mandatory entries
+        type            exponential;
+
+        // Optional entries
+        a               <scalar>;
+        b               <scalar>;
+    }
+    \endverbatim
+
+    where the entries mean:
+    \table
+      Property   | Description                       | Type | Reqd | Deflt
+      type       | Type name: exponential            | word | yes  | -
+      a          | Pre-factor coefficient          | scalar | no   | 6309
+      b          | Exponent coefficient            | scalar | no   | 2.52
+    \endtable
+
+SourceFiles
+    exponential.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef exponential_H
+#define exponential_H
+
+#include "nucleateFluxModel.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace wallBoilingModels
+{
+namespace nucleateFluxModels
+{
+
+/*---------------------------------------------------------------------------*\
+                         Class exponential Declaration
+\*---------------------------------------------------------------------------*/
+
+class exponential
+:
+    public nucleateFluxModel
+{
+    // Private Data
+
+        //- Pre-factor coefficient
+        scalar a_;
+
+        //- Exponent coefficient
+        scalar b_;
+
+
+    // Private Member Functions
+
+        //- No copy construct
+        exponential(const exponential&) = delete;
+
+        //- No copy assignment
+        void operator=(const exponential&) = delete;
+
+
+public:
+
+    //- Runtime type information
+    TypeName("exponential");
+
+
+    // Constructors
+
+        //- Construct from a dictionary
+        exponential(const dictionary& dict);
+
+
+    //- Destructor
+    virtual ~exponential() = default;
+
+
+    // Member Functions
+
+        //- Calculate and return the nucleation-site density
+        virtual tmp<scalarField> qNucleate
+        (
+            const phaseModel& liquid,
+            const phaseModel& vapor,
+            const label patchi,
+            const scalarField& Tl,
+            const scalarField& Tsatw,
+            const scalarField& L
+        ) const;
+
+        //- Write
+        virtual void write(Ostream& os) const;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace nucleateFluxModels
+} // End namespace wallBoilingModels
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/phaseSystemModels/reactingEuler/multiphaseSystem/derivedFvPatchFields/wallBoilingSubModels/nucleateFluxModels/nucleateFluxModel/nucleateFluxModel.C

@@ -0,0 +1,78 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2021 OpenCFD Ltd
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "nucleateFluxModel.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+    namespace wallBoilingModels
+    {
+        defineTypeNameAndDebug(nucleateFluxModel, 0);
+        defineRunTimeSelectionTable(nucleateFluxModel, dictionary);
+    }
+}
+
+// * * * * * * * * * * * * * * * * Selectors * * * * * * * * * * * * * * * * //
+
+Foam::autoPtr<Foam::wallBoilingModels::nucleateFluxModel>
+Foam::wallBoilingModels::nucleateFluxModel::New
+(
+    const dictionary& dict
+)
+{
+    const word modelType(dict.get<word>("type"));
+
+    Info<< "Selecting nucleateFluxModel: " << modelType << endl;
+
+    auto cstrIter = dictionaryConstructorTablePtr_->cfind(modelType);
+
+    if (!cstrIter.found())
+    {
+        FatalIOErrorInLookup
+        (
+            dict,
+            "nucleateFluxModel",
+            modelType,
+            *dictionaryConstructorTablePtr_
+        ) << abort(FatalIOError);
+    }
+
+    return cstrIter()(dict);
+}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+void Foam::wallBoilingModels::nucleateFluxModel::write(Ostream& os) const
+{
+    os.writeEntry("type", this->type());
+}
+
+
+// ************************************************************************* //

src/phaseSystemModels/reactingEuler/multiphaseSystem/derivedFvPatchFields/wallBoilingSubModels/nucleateFluxModels/nucleateFluxModel/nucleateFluxModel.H

@@ -0,0 +1,118 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2021 OpenCFD Ltd
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::wallBoilingModels::nucleateFluxModel
+
+Description
+    Base class for nucleation flux models
+
+SourceFiles
+    nucleateFluxModel.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef nucleateFluxModel_H
+#define nucleateFluxModel_H
+
+#include "volFields.H"
+#include "dictionary.H"
+#include "runTimeSelectionTables.H"
+#include "phaseModel.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace wallBoilingModels
+{
+/*---------------------------------------------------------------------------*\
+                         Class nucleateFluxModel Declaration
+\*---------------------------------------------------------------------------*/
+
+class nucleateFluxModel
+{
+public:
+
+    //- Runtime type information
+    TypeName("nucleateFluxModel");
+
+
+    //- Declare runtime construction
+    declareRunTimeSelectionTable
+    (
+        autoPtr,
+        nucleateFluxModel,
+        dictionary,
+        (
+            const dictionary& dict
+        ),
+        (dict)
+    );
+
+
+    // Generated Methods
+
+        //- Default construct
+        nucleateFluxModel() = default;
+
+        //- Destructor
+        virtual ~nucleateFluxModel() = default;
+
+
+    // Selectors
+
+        //- Select default constructed
+        static autoPtr<nucleateFluxModel> New(const dictionary& dict);
+
+
+    // Member Functions
+
+        //- Calculate nucleate heat flux
+        virtual tmp<scalarField> qNucleate
+        (
+            const phaseModel& liquid,
+            const phaseModel& vapor,
+            const label patchi,
+            const scalarField& Tl,
+            const scalarField& Tsatw,
+            const scalarField& L
+        ) const = 0;
+
+        //- Write
+        virtual void write(Ostream& os) const;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace wallBoilingModels
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/phaseSystemModels/reactingEuler/multiphaseSystem/derivedFvPatchFields/wallBoilingSubModels/nucleationSiteModels/LemmertChawla/LemmertChawla.C

@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2016-2018 OpenFOAM Foundation
-    Copyright (C) 2020 OpenCFD Ltd.
+    Copyright (C) 2020-2021 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -61,12 +61,6 @@ Foam::wallBoilingModels::nucleationSiteModels::LemmertChawla::LemmertChawla
 {}
 
 
-// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
-
-Foam::wallBoilingModels::nucleationSiteModels::LemmertChawla::~LemmertChawla()
-{}
-
-
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
 Foam::tmp<Foam::scalarField>
@@ -83,7 +77,17 @@ Foam::wallBoilingModels::nucleationSiteModels::LemmertChawla::N
     const fvPatchScalarField& Tw =
         liquid.thermo().T().boundaryField()[patchi];
 
-    return Cn_*9.922e5*pow(max((Tw - Tsatw)/10, scalar(0)), 1.805);
+    return Cn_*9.922e5*pow(max((Tw - Tsatw)/scalar(10), scalar(0)), 1.805);
+}
+
+
+void Foam::wallBoilingModels::nucleationSiteModels::LemmertChawla::write
+(
+    Ostream& os
+) const
+{
+    nucleationSiteModel::write(os);
+    os.writeEntry("Cn", Cn_);
 }
 
 

src/phaseSystemModels/reactingEuler/multiphaseSystem/derivedFvPatchFields/wallBoilingSubModels/nucleationSiteModels/LemmertChawla/LemmertChawla.H

@@ -6,6 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2016-2018 OpenFOAM Foundation
+    Copyright (C) 2021 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -27,10 +28,11 @@ Class
     Foam::wallBoilingModels::nucleationSiteModels::LemmertChawla
 
 Description
-    Lemmert & Chawla function for nucleation site density,
-    correlation by Egorov & Menter.
+    A model for nucleation site density based on
+    Lemmert-Chawla (1977) function and
+    Egorov-Menter (2004) correlation for boiling flows.
 
-    References:
+    Reference:
     \verbatim
         Lemmert, M., & Chawla, J. M. (1977).
         Influence of flow velocity on surface boiling heat transfer coefficient.
@@ -42,6 +44,26 @@ Description
         Technical Report ANSYS/TR-04-10, ANSYS Gmbh.
     \endverbatim
 
+Usage
+    Example of the model specification:
+    \verbatim
+    nucleationSiteModel
+    {
+        // Mandatory entries
+        type            LemmertChawla;
+
+        // Optional entries
+        Cn              <scalar>;
+    }
+    \endverbatim
+
+    where the entries mean:
+    \table
+      Property  | Description                         | Type | Reqd | Deflt
+      type      | Type name: LemmertChawla            | word | yes  | -
+      Cn  | Coefficient for nucleation site density | scalar | no   | 1.0
+    \endtable
+
 SourceFiles
     LemmertChawla.C
 
@@ -69,17 +91,27 @@ class LemmertChawla
 :
     public nucleationSiteModel
 {
-
-    // Private data:
+    // Private Data
 
         //- Coefficient for nucleation site density
         scalar Cn_;
 
+
+    // Private Member Functions
+
+        //- No copy construct
+        LemmertChawla(const LemmertChawla&) = delete;
+
+        //- No copy assignment
+        void operator=(const LemmertChawla&) = delete;
+
+
 public:
 
     //- Runtime type information
     TypeName("LemmertChawla");
 
+
     // Constructors
 
         //- Construct from a dictionary
@@ -87,7 +119,7 @@ public:
 
 
     //- Destructor
-    virtual ~LemmertChawla();
+    virtual ~LemmertChawla() = default;
 
 
     // Member Functions
@@ -102,6 +134,11 @@ public:
             const scalarField& Tsatw,
             const scalarField& L
         ) const;
+
+        // I-O
+
+            //- Write
+            virtual void write(Ostream& os) const;
 };
 
 

src/phaseSystemModels/reactingEuler/multiphaseSystem/derivedFvPatchFields/wallBoilingSubModels/nucleationSiteModels/nucleationSiteModel/nucleationSiteModel.H

@@ -6,6 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2016-2018 OpenFOAM Foundation
+    Copyright (C) 2021 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -27,7 +28,7 @@ Class
     Foam::wallBoilingModels::nucleationSiteModel
 
 Description
-    Base class for nucleation site density models
+    Base class for nucleation site density models.
 
 SourceFiles
     nucleationSiteModel.C
@@ -105,7 +106,10 @@ public:
             const scalarField& L
         ) const = 0;
 
-        virtual void write(Ostream& os) const;
+        // I-O
+
+            //- Write
+            virtual void write(Ostream& os) const;
 };
 
 

src/phaseSystemModels/reactingEuler/multiphaseSystem/derivedFvPatchFields/wallBoilingSubModels/partitioningModels/Lavieville/Lavieville.H

@@ -6,6 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2016-2018 OpenFOAM Foundation
+    Copyright (C) 2021 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -27,19 +28,35 @@ Class
     Foam::wallBoilingModels::partitioningModels::Lavieville
 
 Description
-    Lavieville wall heat flux partitioning model.
-
-    Model parameters:
-        alphaCrit: critical liquid fraction
+    A model for wall heat flux partitioning based on
+    Lavieville et al. (2006) for boiling flows.
 
     Reference:
     \verbatim
-        Lavieville, J., Quemerais, E., Mimouni, S., Boucker, M., &
-        Mechitoua, N. (2006).
-        NEPTUNE CFD V1. 0 theory manual.
+        Lavieville, J., Quemerais, E., Mimouni, S.,
+        Boucker, M., & Mechitoua, N. (2006).
+        NEPTUNE CFD V1.0 theory manual.
         NEPTUNE report Nept_2004_L1, 2(3).
     \endverbatim
 
+Usage
+    Example of the model specification:
+    \verbatim
+    partitioningModel
+    {
+        // Mandatory entries
+        type            Lavieville;
+        alphaCrit       <scalar>;
+    }
+    \endverbatim
+
+    where the entries mean:
+    \table
+      Property      | Description                       | Type | Reqd | Deflt
+      type          | Type name: Lavieville             | word | yes  | -
+      alphaCrit     | Critical liquid fraction        | scalar | yes  | -
+    \endtable
+
 SourceFiles
     Lavieville.C
 

src/phaseSystemModels/reactingEuler/multiphaseSystem/derivedFvPatchFields/wallBoilingSubModels/partitioningModels/cosine/cosine.C

@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2016-2019 OpenFOAM Foundation
-    Copyright (C) 2020 OpenCFD Ltd.
+    Copyright (C) 2020-2021 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -60,13 +60,6 @@ cosine::cosine(const dictionary& dict)
 {}
 
 
-// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
-
-Foam::wallBoilingModels::partitioningModels::
-cosine::~cosine()
-{}
-
-
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
 Foam::tmp<Foam::scalarField>
@@ -81,9 +74,9 @@ cosine::fLiquid
        *(
             neg(alphaLiquid0_ - alphaLiquid)
            *(
-                0.5
+                scalar(0.5)
                *(
-                    1 - cos
+                    scalar(1) - cos
                     (
                         constant::mathematical::pi
                        *(alphaLiquid - alphaLiquid0_)

src/phaseSystemModels/reactingEuler/multiphaseSystem/derivedFvPatchFields/wallBoilingSubModels/partitioningModels/cosine/cosine.H

@@ -6,6 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2016-2018 OpenFOAM Foundation
+    Copyright (C) 2021 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -27,20 +28,46 @@ Class
     Foam::wallBoilingModels::partitioningModels::cosine
 
 Description
-    Cosine wall heat flux partitioning model.
+    A cosine model for wall heat flux partitioning based on
+    Tentner et al. (2006) for boiling flows.
 
-    Proposed threshold liquid fractions:
+    Reference:
+    \verbatim
+        Tentner, A., Lo, S., Ioilev, A., Samigulin, M.,
+        Ustinenko, V., Melnikov, V., Kozlov, V. (2006).
+        Advances in computational fluid dynamics modeling
+        of two phase flow in a boiling water reactor fuel assembly.
+        In: Proc. Int. Conf. Nuclear Engineering ICONE-14.
+        Miami, Florida, USA, July 17–20.
+    \endverbatim
+
+Usage
+    Example of the model specification:
+    \verbatim
+    partitioningModel
+    {
+        // Mandatory entries
+        type            cosine;
+        alphaLiquid1    <scalar>;
+        alphaLiquid0    <scalar>;
+    }
+    \endverbatim
+
+    where the entries mean:
+    \table
+      Property      | Description                       | Type | Reqd | Deflt
+      type          | Type name: cosine                 | word | yes  | -
+      alphaLiquid1  | Model parameters for threshold liquid phase fraction <!--
+                    -->                               | scalar | yes  | -
+      alphaLiquid0  | Model parameters for threshold liquid phase fraction <!--
+                    -->                               | scalar | yes  | -
+    \endtable
+
+Note
+  - Proposed threshold liquid fractions:
       - alphaLiquid1 0.1
       - alphaLiquid0 0.05
 
-    \verbatim
-        Tentner, A., Lo, S., & Kozlov, V. (2006).
-        Advances in computational fluid dynamics modeling
-        of two-phase flow in boiling water reactor fuel assemblies.
-        In International Conference of Nuclear Engineering,
-        Miami, Florida, USA.
-    \endverbatim
-
 SourceFiles
     cosine.C
 
@@ -68,13 +95,24 @@ class cosine
 :
     public partitioningModel
 {
-    // Private data
+    // Private Data
 
-        // Model parameters, threshold liquid phase fractions
+        //- Model parameters for threshold liquid phase fraction
         scalar alphaLiquid1_;
+
+        //- Model parameters for threshold liquid phase fraction
         scalar alphaLiquid0_;
 
 
+    // Private Member Functions
+
+        //- No copy construct
+        cosine(const cosine&) = delete;
+
+        //- No copy assignment
+        void operator=(const cosine&) = delete;
+
+
 public:
 
     //- Runtime type information
@@ -88,7 +126,7 @@ public:
 
 
     //- Destructor
-    virtual ~cosine();
+    virtual ~cosine() = default;
 
 
     // Member Functions
@@ -96,7 +134,10 @@ public:
         //- Calculate and return the wall heat-flux partitioning
         virtual tmp<scalarField> fLiquid(const scalarField& alphaLiquid) const;
 
-        virtual void write(Ostream& os) const;
+        // I-O
+
+            //- Write
+            virtual void write(Ostream& os) const;
 };
 
 

src/phaseSystemModels/reactingEuler/multiphaseSystem/derivedFvPatchFields/wallBoilingSubModels/partitioningModels/linear/linear.C

@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2016-2019 OpenFOAM Foundation
-    Copyright (C) 2020 OpenCFD Ltd.
+    Copyright (C) 2020-2021 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -60,13 +60,6 @@ linear::linear(const dictionary& dict)
 {}
 
 
-// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
-
-Foam::wallBoilingModels::partitioningModels::
-linear::~linear()
-{}
-
-
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
 Foam::tmp<Foam::scalarField>

src/phaseSystemModels/reactingEuler/multiphaseSystem/derivedFvPatchFields/wallBoilingSubModels/partitioningModels/linear/linear.H

@@ -6,6 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2016-2018 OpenFOAM Foundation
+    Copyright (C) 2021 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -27,21 +28,47 @@ Class
     Foam::wallBoilingModels::partitioningModels::linear
 
 Description
-    Linear wall heat flux partitioning model.
-
-    Proposed threshold liquid fractions:
-      - alphaLiquid1 0.1
-      - alphaLiquid0 0.05
+    A linear model for wall heat flux partitioning based on
+    Ioilev et al. (2007) for boiling flows.
 
+    Reference:
     \verbatim
-        Ioilev, A., Samigulin, M., Ustinenko (2007).
+        Ioilev, A., Samigulin, M., Ustinenko, V.,
+        Kucherova, P., Tentner, A., Lo, S., & Splawski, A. (2007).
         Advances in the modeling of cladding heat transfer
         and critical heat flux in boiling water reactor fuel assemblies.
-        In Proc. 12th International Topical Meeting on
-        Nuclear Reactor Thermal Hydraulics (NURETH-12),
+        In Proc. 12th International Topical Meeting on Nuclear Reactor
+        Thermal Hydraulics (NURETH-12).
         Pittsburgh, Pennsylvania, USA.
     \endverbatim
 
+Usage
+    Example of the model specification:
+    \verbatim
+    partitioningModel
+    {
+        // Mandatory entries
+        type            linear;
+        alphaLiquid1    <scalar>;
+        alphaLiquid0    <scalar>;
+    }
+    \endverbatim
+
+    where the entries mean:
+    \table
+      Property      | Description                       | Type | Reqd | Deflt
+      type          | Type name: linear                 | word | yes  | -
+      alphaLiquid1  | Model parameters for threshold liquid phase fraction <!--
+                    -->                               | scalar | yes  | -
+      alphaLiquid0  | Model parameters for threshold liquid phase fraction <!--
+                    -->                               | scalar | yes  | -
+    \endtable
+
+Note
+  - Proposed threshold liquid fractions:
+      - alphaLiquid1 0.1
+      - alphaLiquid0 0.05
+
 SourceFiles
     linear.C
 
@@ -69,13 +96,24 @@ class linear
 :
     public partitioningModel
 {
-    // Private data
+    // Private Data
 
-        //- Model parameters, threshold liquid phase fractions
+        //- Model parameters for threshold liquid phase fraction
         scalar alphaLiquid1_;
+
+        //- Model parameters for threshold liquid phase fraction
         scalar alphaLiquid0_;
 
 
+    // Private Member Functions
+
+        //- No copy construct
+        linear(const linear&) = delete;
+
+        //- No copy assignment
+        void operator=(const linear&) = delete;
+
+
 public:
 
     //- Runtime type information
@@ -89,7 +127,7 @@ public:
 
 
     //- Destructor
-    virtual ~linear();
+    virtual ~linear() = default;
 
 
     // Member Functions
@@ -97,7 +135,10 @@ public:
         //- Calculate and return the wall heat-flux partitioning
         virtual tmp<scalarField> fLiquid(const scalarField& alphaLiquid) const;
 
-        virtual void write(Ostream& os) const;
+        // I-O
+
+            //- Write
+            virtual void write(Ostream& os) const;
 };
 
 

src/phaseSystemModels/reactingEuler/multiphaseSystem/derivedFvPatchFields/wallBoilingSubModels/partitioningModels/partitioningModel/partitioningModel.H

@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2016-2018 OpenFOAM Foundation
-    Copyright (C) 2020 OpenCFD Ltd.
+    Copyright (C) 2020-2021 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -28,7 +28,7 @@ Class
     Foam::wallBoilingModels::partitioningModel
 
 Description
-    Base class for wall heat flux partitioning models
+    Base class for wall heat flux partitioning models.
 
 SourceFiles
     partitioningModel.C
@@ -98,7 +98,10 @@ public:
             const scalarField& alphaLiquid
         ) const = 0;
 
-        virtual void write(Ostream& os) const;
+        // I-O
+
+            //- Write
+            virtual void write(Ostream& os) const;
 };
 
 

src/phaseSystemModels/reactingEuler/multiphaseSystem/derivedFvPatchFields/wallBoilingSubModels/partitioningModels/phaseFraction/phaseFraction.C

@@ -6,6 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2016-2018 OpenFOAM Foundation
+    Copyright (C) 2021 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -57,13 +58,6 @@ phaseFraction::phaseFraction(const dictionary& dict)
 {}
 
 
-// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
-
-Foam::wallBoilingModels::partitioningModels::
-phaseFraction::~phaseFraction()
-{}
-
-
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
 Foam::tmp<Foam::scalarField>

src/phaseSystemModels/reactingEuler/multiphaseSystem/derivedFvPatchFields/wallBoilingSubModels/partitioningModels/phaseFraction/phaseFraction.H

@@ -6,6 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2016-2018 OpenFOAM Foundation
+    Copyright (C) 2021 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -56,6 +57,14 @@ class phaseFraction
 :
     public partitioningModel
 {
+    // Private Member Functions
+
+        //- No copy construct
+        phaseFraction(const phaseFraction&) = delete;
+
+        //- No copy assignment
+        void operator=(const phaseFraction&) = delete;
+
 
 public:
 
@@ -70,7 +79,7 @@ public:
 
 
     //- Destructor
-    virtual ~phaseFraction();
+    virtual ~phaseFraction() = default;
 
 
     // Member Functions

src/thermophysicalModels/basic/basicThermo/basicThermo.H

@@ -443,6 +443,14 @@ public:
                 const label patchi
             ) const = 0;
 
+            //- Heat capacity using pressure and temperature [J/kg/K]
+            virtual tmp<scalarField> Cp
+            (
+                const scalarField& p,
+                const scalarField& T,
+                const labelList& cells
+            ) const = 0;
+
             //- Heat capacity at constant volume [J/kg/K]
             virtual tmp<volScalarField> Cv() const = 0;
 
@@ -454,6 +462,14 @@ public:
                 const label patchi
             ) const = 0;
 
+            //- Density from pressure and temperature from EoS
+            virtual tmp<scalarField> rhoEoS
+            (
+                const scalarField& p,
+                const scalarField& T,
+                const labelList& cells
+            ) const = 0;
+
             //- Gamma = Cp/Cv []
             virtual tmp<volScalarField> gamma() const = 0;
 

src/thermophysicalModels/basic/heThermo/heThermo.C

@@ -376,6 +376,28 @@ Foam::tmp<Foam::scalarField> Foam::heThermo<BasicThermo, MixtureType>::Cp
 }
 
 
+template<class BasicThermo, class MixtureType>
+Foam::tmp<Foam::scalarField>
+Foam::heThermo<BasicThermo, MixtureType>::Cp
+(
+    const scalarField& p,
+    const scalarField& T,
+    const labelList& cells
+) const
+{
+    auto tCp = tmp<scalarField>::New(T.size());
+    auto& Cp = tCp.ref();
+
+    forAll(cells, i)
+    {
+        const label celli = cells[i];
+        Cp[i] = this->cellMixture(celli).Cp(p[i], T[i]);
+    }
+
+    return tCp;
+}
+
+
 template<class BasicThermo, class MixtureType>
 Foam::tmp<Foam::volScalarField>
 Foam::heThermo<BasicThermo, MixtureType>::Cp() const
@@ -449,6 +471,28 @@ Foam::heThermo<BasicThermo, MixtureType>::Cv
 }
 
 
+template<class BasicThermo, class MixtureType>
+Foam::tmp<Foam::scalarField>
+Foam::heThermo<BasicThermo, MixtureType>::rhoEoS
+(
+    const scalarField& p,
+    const scalarField& T,
+    const labelList& cells
+) const
+{
+    auto tRho = tmp<scalarField>::New(T.size());
+    auto& rho = tRho.ref();
+
+    forAll(cells, i)
+    {
+        const label celli = cells[i];
+        rho[i] = this->cellMixture(celli).rho(p[i], T[i]);
+    }
+
+    return tRho;
+}
+
+
 template<class BasicThermo, class MixtureType>
 Foam::tmp<Foam::volScalarField>
 Foam::heThermo<BasicThermo, MixtureType>::Cv() const

src/thermophysicalModels/basic/heThermo/heThermo.H

@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2017 OpenFOAM Foundation
-    Copyright (C) 2015-2017 OpenCFD Ltd.
+    Copyright (C) 2015-2021 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -228,6 +228,14 @@ public:
                 const label patchi
             ) const;
 
+            //- Heat capacity using pressure and temperature
+            virtual tmp<scalarField> Cp
+            (
+                const scalarField& p,
+                const scalarField& T,
+                const labelList& cells
+            ) const;
+
             //- Heat capacity at constant pressure [J/kg/K]
             virtual tmp<volScalarField> Cp() const;
 
@@ -239,6 +247,14 @@ public:
                 const label patchi
             ) const;
 
+            //- Density from pressure and temperature
+            virtual tmp<scalarField> rhoEoS
+            (
+                const scalarField& p,
+                const scalarField& T,
+                const labelList& cells
+            ) const;
+
             //- Heat capacity at constant volume [J/kg/K]
             virtual tmp<volScalarField> Cv() const;
 

src/thermophysicalModels/specie/thermo/hPolynomial/hPolynomialThermo.C

@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2017 OpenFOAM Foundation
-    Copyright (C) 2020 OpenCFD Ltd.
+    Copyright (C) 2020-2021 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -41,17 +41,22 @@ Foam::hPolynomialThermo<EquationOfState, PolySize>::hPolynomialThermo
     Hf_(dict.subDict("thermodynamics").get<scalar>("Hf")),
     Sf_(dict.subDict("thermodynamics").get<scalar>("Sf")),
     CpCoeffs_(dict.subDict("thermodynamics").lookup(coeffsName("Cp"))),
+    Tref_(dict.subDict("thermodynamics").getOrDefault<scalar>("Tref", Tstd)),
+    Href_(dict.subDict("thermodynamics").getOrDefault<scalar>("Href", 0)),
+    Sref_(dict.subDict("thermodynamics").getOrDefault<scalar>("Sref", 0)),
+    Pref_(dict.subDict("thermodynamics").getOrDefault<scalar>("Pref", Pstd)),
     hCoeffs_(),
     sCoeffs_()
 {
     hCoeffs_ = CpCoeffs_.integral();
     sCoeffs_ = CpCoeffs_.integralMinus1();
+    // Offset h poly so that it is relative to the enthalpy at Tref
+    hCoeffs_[0] +=
+        Href_ - hCoeffs_.value(Tref_) - EquationOfState::H(Pstd, Tref_);
 
-    // Offset h poly so that it is relative to the enthalpy at Tstd
-    hCoeffs_[0] += Hf_ - hCoeffs_.value(Tstd);
-
-    // Offset s poly so that it is relative to the entropy at Tstd
-    sCoeffs_[0] += Sf_ - sCoeffs_.value(Tstd);
+    // Offset s poly so that it is relative to the entropy at Tref
+    sCoeffs_[0] +=
+        Sref_ - sCoeffs_.value(Tref_) - EquationOfState::S(Pstd, Tref_);
 }
 
 
@@ -70,6 +75,10 @@ void Foam::hPolynomialThermo<EquationOfState, PolySize>::write
         os.beginBlock("thermodynamics");
         os.writeEntry("Hf", Hf_);
         os.writeEntry("Sf", Sf_);
+        os.writeEntry("Tref", Tref_);
+        os.writeEntry("Hsref", Href_);
+        os.writeEntry("Sref", Sref_);
+        os.writeEntry("Pref", Pref_);
         os.writeEntry(coeffsName("Cp"), CpCoeffs_);
         os.endBlock();
     }

src/thermophysicalModels/specie/thermo/hPolynomial/hPolynomialThermo.H

@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2017 OpenFOAM Foundation
-    Copyright (C) 2020 OpenCFD Ltd.
+    Copyright (C) 2020-2021 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -141,6 +141,18 @@ class hPolynomialThermo
         //- Specific heat at constant pressure polynomial coeffs
         Polynomial<PolySize> CpCoeffs_;
 
+        //- Reference temperature
+        scalar Tref_;
+
+        //- Reference enthalphy
+        scalar Href_;
+
+        //- Reference entropy
+        scalar Sref_;
+
+        //- Reference pressure
+        scalar Pref_;
+
         //- Enthalpy polynomial coeffs - derived from cp [J/kg]
         //  NOTE: relative to Tstd
         typename Polynomial<PolySize>::intPolyType hCoeffs_;