ENH: Update of alphaBoiling BC, Bromley and tutorial

applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionTwoPhaseEulerFoam/derivedFvPatchFields/turbulentTemperatureTwoPhaseRadCoupledMixed/turbulentTemperatureTwoPhaseRadCoupledMixedFvPatchScalarField.C

@@ -378,8 +378,11 @@ updateCoeffs()
         // TvNbr: vapour Tp
         scalarField c(TcNbr*KDeltaLiqNbr + TvNbr*KDeltaVapNbr);
 
-        valueFraction() = KDeltaNbr/(KDeltaNbr + KDelta);
+        //valueFraction() = KDeltaNbr/(KDeltaNbr + KDelta);
-        refValue() = c/KDeltaNbr;
+        //refValue() = c/KDeltaNbr;
+        scalarField KDeltaLiqVapNbr(KDeltaLiqNbr + KDeltaVapNbr);
+        valueFraction() = KDeltaLiqVapNbr/(KDeltaLiqVapNbr + KDelta);
+        refValue() = c/KDeltaLiqVapNbr;
         refGrad() = (qr + qrNbr)/kappa(Tp);
 
         if (debug)
@@ -393,7 +396,7 @@ updateCoeffs()
                     << " walltemperature "
                     << " min:" << gMin(Tp)
                     << " max:" << gMax(Tp)
-                    << " avg:" << gAverage(Tp)
+                    << " avg:" << gAverage(Tp) << nl
                     << endl;
         }
     }

applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionTwoPhaseEulerFoam/fluid/createFluidFields.H

@@ -21,7 +21,178 @@ const uniformDimensionedVectorField& g = meshObjects::gravity::New(runTime);
 
 PtrList<pimpleControl> pimpleFluid(fluidRegions.size());
 
-// Populate fluid field pointer lists
+
+//Debug Fields
+/*
+PtrList<volScalarField> faceRegimesFluid(fluidRegions.size());
+PtrList<volScalarField> qcFluid(fluidRegions.size());
+PtrList<volScalarField> qFilmFluid(fluidRegions.size());
+PtrList<volScalarField> htcFilmBoilingFluid(fluidRegions.size());
+PtrList<volScalarField> qtbFluid(fluidRegions.size());
+PtrList<volScalarField> qSubCoolFluid(fluidRegions.size());
+PtrList<volScalarField> CHFtotalFluid(fluidRegions.size());
+PtrList<volScalarField> TdnbFluid(fluidRegions.size());
+PtrList<volScalarField> phiFluid(fluidRegions.size());
+
+forAll(fluidRegions, i)
+{
+    faceRegimesFluid.set
+    (
+        i,
+        new volScalarField
+        (
+            IOobject
+            (
+                "faceRegimes",
+                runTime.timeName(),
+                fluidRegions[i],
+                IOobject::NO_READ,
+                IOobject::AUTO_WRITE
+            ),
+            fluidRegions[i],
+            dimensionedScalar(dimless, Zero)
+
+        )
+    );
+    qcFluid.set
+    (
+        i,
+        new volScalarField
+        (
+            IOobject
+            (
+                "qc",
+                runTime.timeName(),
+                fluidRegions[i],
+                IOobject::NO_READ,
+                IOobject::AUTO_WRITE
+            ),
+            fluidRegions[i],
+            dimensionedScalar(dimless, Zero)
+        )
+    );
+    qFilmFluid.set
+    (
+        i,
+        new volScalarField
+        (
+            IOobject
+            (
+                "qFilm",
+                runTime.timeName(),
+                fluidRegions[i],
+                IOobject::NO_READ,
+                IOobject::AUTO_WRITE
+            ),
+            fluidRegions[i],
+            dimensionedScalar(dimless, Zero)
+        )
+    );
+    htcFilmBoilingFluid.set
+    (
+        i,
+        new volScalarField
+        (
+            IOobject
+            (
+                "htcFilmBoiling",
+                runTime.timeName(),
+                fluidRegions[i],
+                IOobject::NO_READ,
+                IOobject::AUTO_WRITE
+            ),
+            fluidRegions[i],
+            dimensionedScalar(dimless, Zero)
+        )
+    );
+    qtbFluid.set
+    (
+        i,
+        new volScalarField
+        (
+            IOobject
+            (
+                "qtb",
+                runTime.timeName(),
+                fluidRegions[i],
+                IOobject::NO_READ,
+                IOobject::AUTO_WRITE
+            ),
+            fluidRegions[i],
+            dimensionedScalar(dimless, Zero)
+        )
+    );
+    qSubCoolFluid.set
+    (
+        i,
+        new volScalarField
+        (
+            IOobject
+            (
+                "qSubCool",
+                runTime.timeName(),
+                fluidRegions[i],
+                IOobject::NO_READ,
+                IOobject::AUTO_WRITE
+            ),
+            fluidRegions[i],
+            dimensionedScalar(dimless, Zero)
+        )
+    );
+    CHFtotalFluid.set
+    (
+        i,
+        new volScalarField
+        (
+            IOobject
+            (
+                "CHFtotal",
+                runTime.timeName(),
+                fluidRegions[i],
+                IOobject::NO_READ,
+                IOobject::AUTO_WRITE
+            ),
+            fluidRegions[i],
+            dimensionedScalar(dimless, Zero)
+        )
+    );
+     TdnbFluid.set
+    (
+        i,
+        new volScalarField
+        (
+            IOobject
+            (
+                "Tdnb",
+                runTime.timeName(),
+                fluidRegions[i],
+                IOobject::NO_READ,
+                IOobject::AUTO_WRITE
+            ),
+            fluidRegions[i],
+            dimensionedScalar(dimless, Zero)
+        )
+    );
+     phiFluid.set
+    (
+        i,
+        new volScalarField
+        (
+            IOobject
+            (
+                "phiTb",
+                runTime.timeName(),
+                fluidRegions[i],
+                IOobject::NO_READ,
+                IOobject::AUTO_WRITE
+            ),
+            fluidRegions[i],
+            dimensionedScalar(dimless, Zero)
+        )
+    );
+}
+*/
+
 forAll(fluidRegions, i)
 {
     Info<< "*** Reading fluid mesh thermophysical properties for region "

applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionTwoPhaseEulerFoam/fluid/setRegionFluidFields.H

@@ -56,3 +56,16 @@
     pimpleControl& pimple = pimpleFluid[i];
 
     const dimensionedScalar& pMin = pMinFluid[i];
+
+    // Debug fields
+/*
+    volScalarField& faceRegimes = faceRegimesFluid[i];
+    volScalarField& qc = qcFluid[i];
+    volScalarField& qFilm = qFilmFluid[i];
+    volScalarField& htcFilmBoiling = htcFilmBoilingFluid[i];
+    volScalarField& qtb = qtbFluid[i];
+    volScalarField& qSubCool = qSubCoolFluid[i];
+    volScalarField& CHFtotal = CHFtotalFluid[i];
+    volScalarField& Tdnb = TdnbFluid[i];
+    volScalarField& phiTb = phiFluid[i];
+*/

src/phaseSystemModels/reactingEulerFoam/derivedFvPatchFields/alphatWallBoilingWallFunction/alphatWallBoilingWallFunctionFvPatchScalarField.C

@@ -74,7 +74,7 @@ alphatWallBoilingWallFunctionFvPatchScalarField
     alphatPhaseChangeJayatillekeWallFunctionFvPatchScalarField(p, iF),
     otherPhaseName_("vapor"),
     phaseType_(liquidPhase),
-    relax_(0.5),
+    relax_(),
     AbyV_(p.size(), 0),
     alphatConv_(p.size(), 0),
     dDep_(p.size(), 1e-5),
@@ -112,7 +112,7 @@ alphatWallBoilingWallFunctionFvPatchScalarField
     alphatPhaseChangeJayatillekeWallFunctionFvPatchScalarField(p, iF, dict),
     otherPhaseName_(dict.lookup("otherPhase")),
     phaseType_(phaseTypeNames_.read(dict.lookup("phaseType"))),
-    relax_(dict.lookupOrDefault<scalar>("relax", 0.5)),
+    relax_(Function1<scalar>::New("relax", dict)),
     AbyV_(p.size(), 0),
     alphatConv_(p.size(), 0),
     dDep_(p.size(), 1e-5),
@@ -155,20 +155,6 @@ alphatWallBoilingWallFunctionFvPatchScalarField
             const dictionary* LeidenfrostDict =
                 dict.findDict("LeidenfrostModel");
 
-            if (LeidenfrostDict)
-            {
-                LeidenfrostModel_ =
-                    wallBoilingModels::LeidenfrostModel::New(*LeidenfrostDict);
-            }
-
-            const dictionary* filmDict = dict.findDict("filmBoilingModel");
-
-            if (filmDict)
-            {
-                filmBoilingModel_ =
-                    wallBoilingModels::filmBoilingModel::New(*filmDict);
-            }
-
             dmdt_ = 0;
 
             break;
@@ -304,7 +290,7 @@ alphatWallBoilingWallFunctionFvPatchScalarField
     ),
     otherPhaseName_(psf.otherPhaseName_),
     phaseType_(psf.phaseType_),
-    relax_(psf.relax_),
+    relax_(psf.relax_.clone()),
     AbyV_(psf.AbyV_),
     alphatConv_(psf.alphatConv_, mapper),
     dDep_(psf.dDep_, mapper),
@@ -332,7 +318,7 @@ alphatWallBoilingWallFunctionFvPatchScalarField
     alphatPhaseChangeJayatillekeWallFunctionFvPatchScalarField(psf),
     otherPhaseName_(psf.otherPhaseName_),
     phaseType_(psf.phaseType_),
-    relax_(psf.relax_),
+    relax_(psf.relax_.clone()),
     AbyV_(psf.AbyV_),
     alphatConv_(psf.alphatConv_),
     dDep_(psf.dDep_),
@@ -361,7 +347,7 @@ alphatWallBoilingWallFunctionFvPatchScalarField
     alphatPhaseChangeJayatillekeWallFunctionFvPatchScalarField(psf, iF),
     otherPhaseName_(psf.otherPhaseName_),
     phaseType_(psf.phaseType_),
-    relax_(psf.relax_),
+    relax_(psf.relax_.clone()),
     AbyV_(psf.AbyV_),
     alphatConv_(psf.alphatConv_),
     dDep_(psf.dDep_),
@@ -457,6 +443,9 @@ void alphatWallBoilingWallFunctionFvPatchScalarField::updateCoeffs()
 
     const label patchi = patch().index();
 
+    const scalar t = this->db().time().timeOutputValue();
+    scalar relax = relax_->value(t);
+
     switch (phaseType_)
     {
         case vaporPhase:
@@ -469,102 +458,9 @@ void alphatWallBoilingWallFunctionFvPatchScalarField::updateCoeffs()
             const fvPatchScalarField& hewv =
                 vapor.thermo().he().boundaryField()[patchi];
 
-            const phaseModel& liquid(fluid.phases()[otherPhaseName_]);
-
-            const phaseCompressibleTurbulenceModel& turbModel =
-                db().lookupObject<phaseCompressibleTurbulenceModel>
-                (
-                    IOobject::groupName
-                    (
-                        turbulenceModel::propertiesName,
-                        liquid.name()
-                    )
-                );
-
-            const phaseCompressibleTurbulenceModel& vaporTurbModel =
-                db().lookupObject<phaseCompressibleTurbulenceModel>
-                (
-                    IOobject::groupName
-                    (
-                        turbulenceModel::propertiesName,
-                        vapor.name()
-                    )
-                );
-
-            const fvPatchScalarField& rhoVaporw =
-                vaporTurbModel.rho().boundaryField()[patchi];
-
             // Vapor Liquid phase fraction at the wall
             const scalarField vaporw(vapor.boundaryField()[patchi]);
 
-            const fvPatchScalarField& Tw =
-                liquid.thermo().T().boundaryField()[patchi];
-            const scalarField Tc(Tw.patchInternalField());
-
-             // Saturation temperature
-            const tmp<volScalarField> tTsat =
-                satModel.Tsat(liquid.thermo().p());
-            const volScalarField& Tsat = tTsat();
-            const fvPatchScalarField& Tsatw(Tsat.boundaryField()[patchi]);
-            const scalarField Tsatc(Tsatw.patchInternalField());
-
-            const fvPatchScalarField& hewl =
-                liquid.thermo().he().boundaryField()[patchi];
-
-            const fvPatchScalarField& pw =
-                liquid.thermo().p().boundaryField()[patchi];
-
-            const fvPatchScalarField& rhow =
-                turbModel.rho().boundaryField()[patchi];
-
-            const scalarField hw
-            (
-                liquid.thermo().he().member() == "e"
-              ? hewl.patchInternalField() + pw/rhow.patchInternalField()
-              : hewl.patchInternalField()
-            );
-
-            const scalarField L
-            (
-                vapor.thermo().he().member() == "e"
-              ? vapor.thermo().he(pw, Tsatc, patchi) + pw/rhoVaporw - hw
-              : vapor.thermo().he(pw, Tsatc, patchi) - hw
-            );
-
-            // Film boiling models
-
-            scalarField htcFilmBoiling(this->size(), 0);
-            scalarField TLeiden(this->size(), GREAT);
-
-            if (filmBoilingModel_.valid() && LeidenfrostModel_.valid())
-            {
-                // htc for film boiling
-                htcFilmBoiling =
-                    filmBoilingModel_->htcFilmBoil
-                    (
-                        liquid,
-                        vapor,
-                        patchi,
-                        Tc,
-                        Tsatw,
-                        L
-                    );
-
-                // Leidenfrost Temperature
-                TLeiden =
-                    LeidenfrostModel_->TLeid
-                    (
-                        liquid,
-                        vapor,
-                        patchi,
-                        Tc,
-                        Tsatw,
-                        L
-                    );
-            }
-
-            const scalarField qFilm(htcFilmBoiling*max(Tw - Tsatw, scalar(0)));
-
             // NOTE! Assumes 1-thisPhase for liquid fraction in
             // multiphase simulations
             const scalarField fLiquid
@@ -581,24 +477,6 @@ void alphatWallBoilingWallFunctionFvPatchScalarField::updateCoeffs()
             const scalarField alphatv(calcAlphat(*this));
 
             forAll (*this, i)
-            {
-                if (Tw[i] > TLeiden[i])
-                {
-                    this->operator[](i) =
-                    (
-                        max
-                        (
-                            (1 - fLiquid[i])
-                           *(
-                               (qFilm[i]/heSnGrad[i])
-                              /max(vaporw[i], scalar(1e-8))
-                             - alphaw[i]
-                            ),
-                            -alphaw[i]
-                        )
-                    );
-                }
-                else
             {
                 this->operator[](i) =
                 (
@@ -606,7 +484,6 @@ void alphatWallBoilingWallFunctionFvPatchScalarField::updateCoeffs()
                    /max(vaporw[i], scalar(1e-8))
                 );
             }
-            }
 
             if (debug)
             {
@@ -759,11 +636,12 @@ void alphatWallBoilingWallFunctionFvPatchScalarField::updateCoeffs()
 
             const scalarField fLiquid(partitioningModel_->fLiquid(liquidw));
 
-            for (label i=0; i<10; i++)
-            {
             // 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_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);
@@ -872,12 +750,42 @@ void alphatWallBoilingWallFunctionFvPatchScalarField::updateCoeffs()
                         (
                             wp_*(Tw - tDNB)/(TLeiden - tDNB),
                             scalar(0)
-                            )
+                        ),
-                            , scalar(1)
+                        scalar(1)
                     )
                 );
 
                 Qtb = CHFtotal*(1 - phi) + phi*MHF;
+
+
+                if (debug & 2)
+                {
+
+                    scalarField& phiFluid =
+                        const_cast<fvPatchScalarField&>
+                        (
+                            patch().lookupPatchField
+                            <
+                                volScalarField,
+                                scalar
+                            >("phiTb")
+                        );
+
+                    phiFluid = phi;
+
+                    scalarField& CHFFluid =
+                        const_cast<fvPatchScalarField&>
+                        (
+                            patch().lookupPatchField
+                            <
+                                volScalarField,
+                                scalar
+                            >("CHFtotal")
+                        );
+
+                    CHFFluid = CHFtotal;
+                }
+
             }
 
 
@@ -941,43 +849,27 @@ void alphatWallBoilingWallFunctionFvPatchScalarField::updateCoeffs()
                         // Sub-cool boiling
                         regimeTypes[i] = regimeType::subcool;
 
-                            Tl = (Tw - (Tplus_y250/Tplus)*(Tw - Tc));
-                            Tl = max(Tc - 40, Tl);
-
-                            // Area fractions:
-
-                            /*
                         // Del Valle & Kenning (1985)
-                            const scalarField Ja
+                        const scalar Ja =
-                            (
+                            rhow[i]*Cpw[i]*(Tsatw[i] - Tl[i])
-                                rhoLiquidw*Cpw*(Tsatw - Tl)/(rhoVaporw*L)
+                            /(rhoVaporw[i]*L[i]);
-                            );
 
-                            const scalarField Al
+                        const scalar Al =
-                            (
+                            fLiquid[i]*4.8*exp(min(-Ja/80, log(VGREAT)));
-                                fLiquid*4.8*exp(min(-Ja/80, log(VGREAT)))
-                            );
-                            */
 
-                            // More simple method to calculate area affected by
-                            // bubbles
                         const scalar A2
                         (
-                                min
+                            min(pi*sqr(dDep_[i])*N[i]*Al/4, scalar(1))
-                                (
-                                    fLiquid[i]*pi*sqr(dDep_[i])*N[i]*K_/4,
-                                    scalar(1)
-                                )
                         );
 
-                            A1[i] = max(1 - A2, 0.0);
+                        A1[i] = max(1 - A2, scalar(1e-4));
 
                         // Following Bowring(1962)
                         const scalar A2E
                         (
                             min
                             (
-                                    fLiquid[i]*pi*sqr(dDep_[i])*N[i],
+                                pi*sqr(dDep_[i])*N[i]*Al/4,
                                 scalar(5)
                             )
                         );
@@ -986,8 +878,8 @@ void alphatWallBoilingWallFunctionFvPatchScalarField::updateCoeffs()
                         // to the wall boiling
                         dmdt_[i] =
                             (
-                                    (1 - relax_)*dmdt_[i]
+                                (1 - relax)*dmdt_[i]
-                                  + relax_*(1.0/6.0)*A2E*dDep_[i]*rhoVaporw[i]
+                                + relax*(1.0/6.0)*A2E*dDep_[i]*rhoVaporw[i]
                                 * fDep[i]*AbyV_[i]
                             );
 
@@ -1001,15 +893,15 @@ void alphatWallBoilingWallFunctionFvPatchScalarField::updateCoeffs()
                             2*(alphaw[i]*Cpw[i])*fDep[i]
                             *sqrt
                             (
-                                    (0.8/fDep[i])/(pi*alphaw[i]/rhow[i])
+                                (0.8/max(fDep[i], SMALL))/(pi*alphaw[i]/rhow[i])
                             )
                         );
 
                         // Quenching heat flux in Sub-cool boiling
                         qq_[i] =
                             (
-                                    (1 - relax_)*qq_[i]
+                                (1 - relax)*qq_[i]
-                                  + relax_*A2*hQ*max(Tw[i] - Tl[i], scalar(0))
+                                + relax*A2*hQ*max(Tw[i] - Tl[i], scalar(0))
                             );
 
                         this->operator[](i) =
@@ -1038,17 +930,17 @@ void alphatWallBoilingWallFunctionFvPatchScalarField::updateCoeffs()
                         dmdt_[i] =
                                 fLiquid[i]
                                 *(
-                                    relax_*Qtb[i]*AbyV_[i]/L[i]
+                                    relax*Qtb[i]*AbyV_[i]/L[i]
-                                  + (1 - relax_)*dmdt_[i]
+                                    + (1 - relax)*dmdt_[i]
                                 );
 
                         mDotL_[i] = dmdt_[i]*L[i];
 
+
                         // No quenching flux
                         qq_[i] = 0.0;
 
                         this->operator[](i) =
-                            (
                         max
                         (
                             (
@@ -1056,8 +948,8 @@ void alphatWallBoilingWallFunctionFvPatchScalarField::updateCoeffs()
                                 /max(hewSn[i], scalar(1e-16))
                             )/max(liquidw[i], scalar(1e-8)),
                             1e-8
-                                )
                         );
+
                     }
                     else if (Tw[i] > TLeiden[i])
                     {
@@ -1070,9 +962,10 @@ void alphatWallBoilingWallFunctionFvPatchScalarField::updateCoeffs()
                         dmdt_[i] =
                             fLiquid[i]
                             *(
-                                    relax_*htcFilmBoiling[i]
+                                relax*htcFilmBoiling[i]
-                                   *max(Tw[i] - Tsatw[i], 0)*AbyV_[i]/L[i]
+                                *max(Tw[i] - Tsatw[i], 0)
-                                 + (1 - relax_)*dmdt_[i]
+                                *AbyV_[i]/L[i]
+                                + (1 - relax)*dmdt_[i]
                             );
 
 
@@ -1087,12 +980,12 @@ void alphatWallBoilingWallFunctionFvPatchScalarField::updateCoeffs()
                         );
 
                         // alphat is added alphal and multiplied by phase
-                            // alphaFilm is lower than alphal. Then we substract
+                        // alphaFilm in the coupled BC. We substract
-                            // alpha and divide by phase to get alphaFilm
+                        // alpha and divide by phase to get a net alphaFilm
                         this->operator[](i) =
                             (
-                                (alphaFilm[i]/max(liquidw[i], scalar(1e-8))
+                                alphaFilm[i]/max(liquidw[i], scalar(1e-8))
-                              - alphaw[i])
+                              - alphaw[i]
                             );
                     }
                 }
@@ -1101,7 +994,7 @@ void alphatWallBoilingWallFunctionFvPatchScalarField::updateCoeffs()
                     // Tw below Tsat. No boiling single phase convection
                     // Single phase
                     regimeTypes[i] = regimeType::nonBoiling;
-                        A1[i] = 1.0;
+
                     qq_[i] = 0.0;
                     mDotL_[i] = 0.0;
                     dmdt_[i] = 0.0;
@@ -1111,22 +1004,14 @@ void alphatWallBoilingWallFunctionFvPatchScalarField::updateCoeffs()
                     (
                         max
                         (
-                                fLiquid[i]*A1[i]*(alphatConv_[i])
+                            fLiquid[i]*(alphatConv_[i])
                             /max(liquidw[i], scalar(1e-8)),
-                                1e-16
+                            1e-8
                         )
                     );
                 }
             }
 
-                scalarField TsupPrev(max((Tw - Tsatw), scalar(0)));
-
-                // NOTE: lagging Tw update.
-                //const_cast<fvPatchScalarField&>(Tw).evaluate();
-                scalarField TsupNew(max((Tw - Tsatw), scalar(0)));
-
-                scalar maxErr(max(mag(TsupPrev - TsupNew)));
-
             if (debug)
             {
                 const scalarField qEff
@@ -1151,6 +1036,18 @@ void alphatWallBoilingWallFunctionFvPatchScalarField::updateCoeffs()
 
                 if (debug & 2)
                 {
+                    /*
+                    scalarField& faceRegimes =
+                        const_cast<fvPatchScalarField&>
+                        (
+                            patch().lookupPatchField
+                            <
+                                volScalarField,
+                                scalar
+                            >("faceRegimes")
+                        );
+                    */
+
                     scalar nSubCool(0);
                     scalar nTransient(0);
                     scalar nFilm(0);
@@ -1163,6 +1060,7 @@ void alphatWallBoilingWallFunctionFvPatchScalarField::updateCoeffs()
 
                     forAll (*this, i)
                     {
+                        //faceRegimes[i] = regimeTypes[i];
                         switch (regimeTypes[i])
                         {
                             case regimeType::subcool:
@@ -1200,15 +1098,54 @@ void alphatWallBoilingWallFunctionFvPatchScalarField::updateCoeffs()
                         nNonBoilings*fLiquid*A1*(alphatConv_ + alphaw)
                         *hew.snGrad()
                     );
+                    /*
+                    scalarField& qcFluid =
+                        const_cast<fvPatchScalarField&>
+                        (
+                            patch().lookupPatchField
+                            <
+                                volScalarField,
+                                scalar
+                            >("qc")
+                        );
+
+                    qcFluid = qc;
+                    */
 
                     scalar Qc = gSum(qc*patch().magSf());
                     Info<< " Convective heat transfer: " << Qc << endl;
 
                     const scalarField qFilm
                     (
-                            relax_*fLiquid*nFilms*htcFilmBoiling*(Tw - Tsatw)
+                        relax*fLiquid*nFilms*htcFilmBoiling*(Tw - Tsatw)
+                    );
+                    /*
+                    scalarField& qFilmFluid =
+                        const_cast<fvPatchScalarField&>
+                        (
+                            patch().lookupPatchField
+                            <
+                                volScalarField,
+                                scalar
+                            >("qFilm")
                         );
 
+                    qFilmFluid = qFilm;
+                    */
+
+                    /*
+                    scalarField& htcFilm =
+                        const_cast<fvPatchScalarField&>
+                    (
+                        patch().lookupPatchField
+                        <
+                            volScalarField,
+                            scalar
+                        >("htcFilmBoiling")
+                    );
+
+                    htcFilm = htcFilmBoiling;
+                    */
                     scalar QFilm = gSum(qFilm*patch().magSf());
                     Info<< " Film boiling heat transfer: " << QFilm << endl;
 
@@ -1216,12 +1153,37 @@ void alphatWallBoilingWallFunctionFvPatchScalarField::updateCoeffs()
                         << gMin(nFilms*htcFilmBoiling)
                         << " - "
                         << gMax(nFilms*htcFilmBoiling) << endl;
+                    /*
+                    scalarField& qtbFluid =
+                        const_cast<fvPatchScalarField&>
+                        (
+                            patch().lookupPatchField
+                            <
+                                volScalarField,
+                                scalar
+                            >("qtb")
+                        );
+                    qtbFluid = fLiquid*Qtb;
+                    */
                     
                     scalar Qtbtot =
                         gSum(fLiquid*nTransients*Qtb*patch().magSf());
                     Info<< " Transient boiling heat transfer:" << Qtbtot
                         << endl;
 
+                    /*
+                    scalarField& TdnbFluid =
+                        const_cast<fvPatchScalarField&>
+                        (
+                            patch().lookupPatchField
+                            <
+                                volScalarField,
+                                scalar
+                            >("Tdnb")
+                        );
+                    TdnbFluid = tDNB;
+                    */
+
                     Info<< " TDNB: " << gMin(tDNB) << " - " << gMax(tDNB)
                         << endl;
 
@@ -1234,6 +1196,20 @@ void alphatWallBoilingWallFunctionFvPatchScalarField::updateCoeffs()
                         )
                     );
                     
+                    /*
+                    scalarField& qSubCoolFluid =
+                        const_cast<fvPatchScalarField&>
+                        (
+                            patch().lookupPatchField
+                            <
+                                volScalarField,
+                                scalar
+                            >("qSubCool")
+                        );
+
+                    qSubCoolFluid = qSubCool;
+                    */
+                    
                     scalar QsubCool = gSum(qSubCool*patch().magSf());
 
                     Info<< " Sub Cool boiling heat transfer: " << QsubCool
@@ -1249,24 +1225,10 @@ void alphatWallBoilingWallFunctionFvPatchScalarField::updateCoeffs()
                         << gMax(nSubCools*A1) << endl;
 
                     Info<< nl;
-
                 }
-
             }
-
-                if (maxErr < 1e-1)
-                {
-                    if (i > 0)
-                    {
-                        Info<< "Wall boiling wall function iterations: "
-                            << i + 1 << endl;
         }
         break;
-                }
-
-            }
-            break;
-        }
         default:
         {
             FatalErrorInFunction
@@ -1286,7 +1248,7 @@ void alphatWallBoilingWallFunctionFvPatchScalarField::write(Ostream& os) const
     os.writeKeyword("phaseType") << phaseTypeNames_[phaseType_]
         << token::END_STATEMENT << nl;
 
-    os.writeKeyword("relax") << relax_ << token::END_STATEMENT << nl;
+    relax_->writeData(os);
 
     switch (phaseType_)
     {

src/phaseSystemModels/reactingEulerFoam/derivedFvPatchFields/alphatWallBoilingWallFunction/alphatWallBoilingWallFunctionFvPatchScalarField.H

@@ -224,6 +224,8 @@ SourceFiles
 #ifndef compressible_alphatWallBoilingWallFunctionFvPatchScalarField_H
 #define compressible_alphatWallBoilingWallFunctionFvPatchScalarField_H
 
+#include "Function1.H"
+
 #include "alphatPhaseChangeJayatillekeWallFunctionFvPatchScalarField.H"
 #include "partitioningModel.H"
 #include "nucleationSiteModel.H"
@@ -287,7 +289,7 @@ private:
         phaseType phaseType_;
 
         //- dmdt relaxationFactor
-        scalar relax_;
+         autoPtr<Function1<scalar>> relax_;
 
         //- Patch face area by cell volume
         scalarField AbyV_;

src/phaseSystemModels/reactingEulerFoam/derivedFvPatchFields/wallBoilingSubModels/TDNBModels/Schroeder/Schroeder.C

@@ -58,7 +58,7 @@ Foam::wallBoilingModels::TDNBModels::Schroeder::Schroeder
 )
 :
     TDNBModel(),
-    kg_(dict.lookupOrDefault<scalar>("kg", 5/3))
+    kg_(dict.lookupOrDefault<scalar>("kg", 1.666))
 {}
 
 

src/phaseSystemModels/reactingEulerFoam/derivedFvPatchFields/wallBoilingSubModels/filmBoilingModels/Bromley/Bromley.C

@@ -26,6 +26,9 @@ License
 #include "Bromley.H"
 #include "addToRunTimeSelectionTable.H"
 #include "uniformDimensionedFields.H"
+#include "constants.H"
+
+using namespace Foam::constant;
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
@@ -55,7 +58,9 @@ Foam::wallBoilingModels::filmBoilingModels::Bromley::Bromley
 )
 :
     filmBoilingModel(),
-    Cn_(dict.lookupOrDefault<scalar>("Cn", 0.62))
+    Cn_(dict.lookupOrDefault<scalar>("Cn", 0.62)),
+    emissivity_(dict.lookupOrDefault<scalar>("emissivity", 1)),
+    L_(dict.get<scalar>("L"))
 {}
 
 
@@ -84,7 +89,11 @@ Foam::wallBoilingModels::filmBoilingModels::Bromley::htcFilmBoil
     const uniformDimensionedVectorField& g =
         liquid.mesh().time().lookupObject<uniformDimensionedVectorField>("g");
 
-    const scalarField rhoVapor(vapor.thermo().rho(patchi));
+    const fvPatchScalarField& rhoVaporw
+    (
+        vapor.thermo().rho()().boundaryField()[patchi]
+    );
+
     const scalarField rhoLiq(liquid.thermo().rho(patchi));
     const scalarField kappaVapor(vapor.kappa(patchi));
 
@@ -93,17 +102,23 @@ Foam::wallBoilingModels::filmBoilingModels::Bromley::htcFilmBoil
     const scalarField& CpVapor = Cp.boundaryField()[patchi];
 
     const scalarField muVapor(vapor.mu(patchi));
-    const scalarField dbVapor(vapor.d()().boundaryField()[patchi]);
+    //const scalarField dbVapor(vapor.d()().boundaryField()[patchi]);
+
+    const scalarField htcRad
+    (
+        emissivity_*physicoChemical::sigma.value()*(pow4(Tw) - pow4(Tsatw))
+      / max((Tw - Tsatw), scalar(1e-4))
+    );
 
     return
         Cn_*pow
         (
             pow3(kappaVapor)
-           *rhoVapor*(rhoLiq - rhoVapor)*mag(g.value())
+           *rhoVaporw*(rhoLiq - rhoVaporw)*mag(g.value())
            *(L + 0.4*CpVapor*max((Tw-Tsatw), scalar(0)))
-           /(dbVapor*muVapor*max((Tw-Tsatw), scalar(1e-4))),
+           /(L_*muVapor*max((Tw-Tsatw), scalar(1e-4))),
             0.25
-        );
+        ) + 0.75*htcRad;
 }
 
 
@@ -114,6 +129,7 @@ void Foam::wallBoilingModels::filmBoilingModels::Bromley::write
 {
     filmBoilingModel::write(os);
     os.writeKeyword("Cn") << Cn_ << token::END_STATEMENT << nl;
+    os.writeKeyword("emissivity") << emissivity_ << token::END_STATEMENT << nl;
 }
 
 

src/phaseSystemModels/reactingEulerFoam/derivedFvPatchFields/wallBoilingSubModels/filmBoilingModels/Bromley/Bromley.H

@@ -66,6 +66,12 @@ class Bromley
         //- Coefficient for nucleation site density
         scalar Cn_;
 
+        //- Wall emissivity
+        scalar emissivity_;
+
+        //- Characteristic lenght-scale
+        scalar L_;
+
 public:
 
     //- Runtime type information

tutorials/heatTransfer/chtMultiRegionTwoPhaseEulerFoam/solidQuenching2D/0.orig/water/alphat.gas

@@ -43,23 +43,13 @@ boundaryField
         Cmu             0.09;
         kappa           0.41;
         E               9.8;
+        relax           constant 1;
         partitioningModel
         {
             type        Lavieville;
             alphaCrit   0.2;
         }
-
+        value           uniform 1e-8;
-        LeidenfrostModel
-        {
-            type        Spiegler;
-            Tcrit       647;
-        }
-
-        filmBoilingModel
-        {
-            type        Bromley;
-        }
-        value           uniform 0;
     }
     maxY
     {

tutorials/heatTransfer/chtMultiRegionTwoPhaseEulerFoam/solidQuenching2D/0.orig/water/alphat.liquid

@@ -43,6 +43,7 @@ boundaryField
         Cmu             0.09;
         kappa           0.41;
         E               9.8;
+        relax           constant 1;
         dmdt            uniform 0;
         partitioningModel
         {
@@ -88,6 +89,7 @@ boundaryField
         filmBoilingModel
         {
             type        Bromley;
+            L           1e-3;
         }
         value           uniform 1e-8;
     }