multiphaseEuler: Remove duplication of thermal wall functions

The standard Jayatilleke thermal wall function now permits evaluation
via static functions. The boiling wall function now uses these
functions, thereby removing the phase-Jayatilleke base class and
associated duplication of the Jayatilleke model.

applications/solvers/modules/multiphaseEuler/multiphaseThermophysicalTransportModels/Make/files

@@ -20,7 +20,6 @@ derivedFvPatchFields/wallBoilingSubModels/departureFrequencyModels/departureFreq
 derivedFvPatchFields/wallBoilingSubModels/departureFrequencyModels/Cole/Cole.C
 derivedFvPatchFields/wallBoilingSubModels/departureFrequencyModels/KocamustafaogullariIshiiDepartureFrequency/KocamustafaogullariIshiiDepartureFrequency.C
 
-derivedFvPatchFields/alphatPhaseJayatillekeWallFunction/alphatPhaseJayatillekeWallFunctionFvPatchScalarField.C
 derivedFvPatchFields/alphatPhaseChangeWallFunction/alphatPhaseChangeWallFunctionFvPatchScalarField.C
 derivedFvPatchFields/alphatFixedDmdtfWallBoilingWallFunction/alphatFixedDmdtfWallBoilingWallFunctionFvPatchScalarField.C
 derivedFvPatchFields/alphatWallBoilingWallFunction/alphatWallBoilingWallFunctionFvPatchScalarField.C

applications/solvers/modules/multiphaseEuler/multiphaseThermophysicalTransportModels/Make/options

@@ -14,6 +14,7 @@ EXE_INC = \
     -I$(LIB_SRC)/MomentumTransportModels/phaseCompressible/lnInclude \
     -I$(LIB_SRC)/ThermophysicalTransportModels/thermophysicalTransportModel/lnInclude \
     -I$(LIB_SRC)/ThermophysicalTransportModels/coupledThermophysicalTransportModels/lnInclude \
+    -I$(LIB_SRC)/ThermophysicalTransportModels/fluid/lnInclude \
     -I$(LIB_SRC)/finiteVolume/lnInclude \
     -I$(LIB_SRC)/meshTools/lnInclude
 

applications/solvers/modules/multiphaseEuler/multiphaseThermophysicalTransportModels/derivedFvPatchFields/alphatFixedDmdtfWallBoilingWallFunction/alphatFixedDmdtfWallBoilingWallFunctionFvPatchScalarField.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2015-2022 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2015-2023 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -24,6 +24,8 @@ License
 \*---------------------------------------------------------------------------*/
 
 #include "alphatFixedDmdtfWallBoilingWallFunctionFvPatchScalarField.H"
+#include "alphatJayatillekeWallFunctionFvPatchScalarField.H"
+#include "fluidThermophysicalTransportModel.H"
 #include "fvPatchFieldMapper.H"
 #include "addToRunTimeSelectionTable.H"
 
@@ -70,13 +72,7 @@ alphatFixedDmdtfWallBoilingWallFunctionFvPatchScalarField
     const fvPatchFieldMapper& mapper
 )
 :
-    alphatPhaseChangeWallFunctionFvPatchScalarField
+    alphatPhaseChangeWallFunctionFvPatchScalarField(psf, p, iF, mapper),
-    (
-        psf,
-        p,
-        iF,
-        mapper
-    ),
     fixedDmdtf_(psf.fixedDmdtf_)
 {}
 
@@ -104,7 +100,18 @@ void alphatFixedDmdtfWallBoilingWallFunctionFvPatchScalarField::updateCoeffs()
 
     dmdtf_ = (1 - relax_)*dmdtf_ + relax_*fixedDmdtf_;
 
-    operator==(calcAlphat(*this));
+    operator==
+    (
+        alphatJayatillekeWallFunctionFvPatchScalarField::alphat
+        (
+            db().lookupType<fluidThermophysicalTransportModel>
+            (
+                internalField().group()
+            ),
+            0.85,
+            patch().index()
+        )
+    );
 
     fixedValueFvPatchScalarField::updateCoeffs();
 }

applications/solvers/modules/multiphaseEuler/multiphaseThermophysicalTransportModels/derivedFvPatchFields/alphatPhaseChangeWallFunction/alphatPhaseChangeWallFunctionFvPatchScalarField.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2015-2022 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2015-2023 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -41,6 +41,7 @@ namespace compressible
 
 defineTypeNameAndDebug(alphatPhaseChangeWallFunctionFvPatchScalarField, 0);
 
+
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
 alphatPhaseChangeWallFunctionFvPatchScalarField::
@@ -50,7 +51,7 @@ alphatPhaseChangeWallFunctionFvPatchScalarField
     const DimensionedField<scalar, volMesh>& iF
 )
 :
-    alphatPhaseJayatillekeWallFunctionFvPatchScalarField(p, iF),
+    fixedValueFvPatchScalarField(p, iF),
     otherPhaseName_(word::null),
     relax_(1),
     dmdtf_(p.size(), 0)
@@ -65,7 +66,7 @@ alphatPhaseChangeWallFunctionFvPatchScalarField
     const dictionary& dict
 )
 :
-    alphatPhaseJayatillekeWallFunctionFvPatchScalarField(p, iF, dict),
+    fixedValueFvPatchScalarField(p, iF, dict),
     otherPhaseName_(dict.lookup("otherPhase")),
     relax_(dict.lookupOrDefault<scalar>("relax", 1)),
     dmdtf_(p.size(), 0)
@@ -97,7 +98,7 @@ alphatPhaseChangeWallFunctionFvPatchScalarField
     const fvPatchFieldMapper& mapper
 )
 :
-    alphatPhaseJayatillekeWallFunctionFvPatchScalarField(ptf, p, iF, mapper),
+    fixedValueFvPatchScalarField(ptf, p, iF, mapper),
     otherPhaseName_(ptf.otherPhaseName_),
     relax_(ptf.relax_),
     dmdtf_(mapper(ptf.dmdtf_))
@@ -111,7 +112,7 @@ alphatPhaseChangeWallFunctionFvPatchScalarField
     const DimensionedField<scalar, volMesh>& iF
 )
 :
-    alphatPhaseJayatillekeWallFunctionFvPatchScalarField(awfpsf, iF),
+    fixedValueFvPatchScalarField(awfpsf, iF),
     otherPhaseName_(awfpsf.otherPhaseName_),
     relax_(awfpsf.relax_),
     dmdtf_(awfpsf.dmdtf_)
@@ -120,23 +121,16 @@ alphatPhaseChangeWallFunctionFvPatchScalarField
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
-bool alphatPhaseChangeWallFunctionFvPatchScalarField::
+bool alphatPhaseChangeWallFunctionFvPatchScalarField::activeInterface
-activeInterface(const phaseInterface& interface) const
+(
+    const phaseInterface& interface
+) const
 {
     const phaseSystem& fluid = interface.fluid();
 
-    if
+    return
-    (
         interface.contains(fluid.phases()[internalField().group()])
-     && interface.contains(fluid.phases()[otherPhaseName_])
+     && interface.contains(fluid.phases()[otherPhaseName_]);
-    )
-    {
-        return true;
-    }
-    else
-    {
-        return false;
-    }
 }
 
 
@@ -147,22 +141,20 @@ alphatPhaseChangeWallFunctionFvPatchScalarField::dmdtf() const
 }
 
 
-const scalarField& alphatPhaseChangeWallFunctionFvPatchScalarField::
+const scalarField& alphatPhaseChangeWallFunctionFvPatchScalarField::dmdtf
-dmdtf(const phaseInterface& interface) const
+(
+    const phaseInterface& interface
+) const
 {
-    if (activeInterface(interface))
+    if (!activeInterface(interface))
-    {
-        return dmdtf_;
-    }
-    else
     {
         FatalErrorInFunction
             << "Phase change mass transfer rate requested for interface on "
             << "which there is no phase change "
             << abort(FatalError);
-
-        return dmdtf_;
     }
+
+    return dmdtf_;
 }
 
 
@@ -171,7 +163,7 @@ void alphatPhaseChangeWallFunctionFvPatchScalarField::autoMap
     const fvPatchFieldMapper& m
 )
 {
-    alphatPhaseJayatillekeWallFunctionFvPatchScalarField::autoMap(m);
+    fixedValueFvPatchScalarField::autoMap(m);
 
     m(dmdtf_, dmdtf_);
 }
@@ -183,7 +175,7 @@ void alphatPhaseChangeWallFunctionFvPatchScalarField::rmap
     const labelList& addr
 )
 {
-    alphatPhaseJayatillekeWallFunctionFvPatchScalarField::rmap(ptf, addr);
+    fixedValueFvPatchScalarField::rmap(ptf, addr);
 
     const alphatPhaseChangeWallFunctionFvPatchScalarField& tiptf =
         refCast<const alphatPhaseChangeWallFunctionFvPatchScalarField>(ptf);
@@ -197,7 +189,7 @@ void alphatPhaseChangeWallFunctionFvPatchScalarField::reset
     const fvPatchScalarField& ptf
 )
 {
-    alphatPhaseJayatillekeWallFunctionFvPatchScalarField::reset(ptf);
+    fixedValueFvPatchScalarField::reset(ptf);
 
     const alphatPhaseChangeWallFunctionFvPatchScalarField& tiptf =
         refCast<const alphatPhaseChangeWallFunctionFvPatchScalarField>(ptf);
@@ -208,7 +200,7 @@ void alphatPhaseChangeWallFunctionFvPatchScalarField::reset
 
 void alphatPhaseChangeWallFunctionFvPatchScalarField::write(Ostream& os) const
 {
-    alphatPhaseJayatillekeWallFunctionFvPatchScalarField::write(os);
+    fixedValueFvPatchScalarField::write(os);
 
     writeEntry(os, "otherPhase", otherPhaseName_);
     writeEntry(os, "relax", relax_);

applications/solvers/modules/multiphaseEuler/multiphaseThermophysicalTransportModels/derivedFvPatchFields/alphatPhaseChangeWallFunction/alphatPhaseChangeWallFunctionFvPatchScalarField.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2015-2022 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2015-2023 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -27,9 +27,6 @@ Class
 Description
     Abstract base-class for all alphatWallFunctions supporting phase-change.
 
-See also
-    Foam::alphatPhaseJayatillekeWallFunctionFvPatchScalarField
-
 SourceFiles
     alphatPhaseChangeWallFunctionFvPatchScalarField.C
 
@@ -38,7 +35,7 @@ SourceFiles
 #ifndef alphatPhaseChangeWallFunctionFvPatchScalarField_H
 #define alphatPhaseChangeWallFunctionFvPatchScalarField_H
 
-#include "alphatPhaseJayatillekeWallFunctionFvPatchScalarField.H"
+#include "fixedValueFvPatchFields.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -56,7 +53,7 @@ namespace compressible
 
 class alphatPhaseChangeWallFunctionFvPatchScalarField
 :
-    public alphatPhaseJayatillekeWallFunctionFvPatchScalarField
+    public fixedValueFvPatchScalarField
 {
 protected:
 

applications/solvers/modules/multiphaseEuler/multiphaseThermophysicalTransportModels/derivedFvPatchFields/alphatPhaseJayatillekeWallFunction/alphatPhaseJayatillekeWallFunctionFvPatchScalarField.C

@@ -1,315 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2015-2022 OpenFOAM Foundation
-     \\/     M anipulation  |
--------------------------------------------------------------------------------
-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 "alphatPhaseJayatillekeWallFunctionFvPatchScalarField.H"
-#include "phaseSystem.H"
-#include "phaseCompressibleMomentumTransportModel.H"
-#include "addToRunTimeSelectionTable.H"
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-namespace Foam
-{
-namespace compressible
-{
-
-// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
-
-scalar alphatPhaseJayatillekeWallFunctionFvPatchScalarField::maxExp_  = 50.0;
-scalar alphatPhaseJayatillekeWallFunctionFvPatchScalarField::tolerance_ = 0.01;
-label alphatPhaseJayatillekeWallFunctionFvPatchScalarField::maxIters_ = 10;
-
-// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
-
-tmp<scalarField> alphatPhaseJayatillekeWallFunctionFvPatchScalarField::Psmooth
-(
-    const scalarField& Prat
-) const
-{
-    return 9.24*(pow(Prat, 0.75) - 1.0)*(1.0 + 0.28*exp(-0.007*Prat));
-}
-
-
-tmp<scalarField>
-alphatPhaseJayatillekeWallFunctionFvPatchScalarField::yPlusTherm
-(
-    const nutWallFunctionFvPatchScalarField& nutw,
-    const scalarField& P,
-    const scalarField& Prat
-) const
-{
-    tmp<scalarField> typsf(new scalarField(this->size()));
-    scalarField& ypsf = typsf.ref();
-
-    forAll(ypsf, facei)
-    {
-        scalar ypt = 11.0;
-
-        for (int i=0; i<maxIters_; i++)
-        {
-            const scalar f =
-                ypt - (log(nutw.E()*ypt)/nutw.kappa() + P[facei])/Prat[facei];
-            const scalar df = 1 - 1.0/(ypt*nutw.kappa()*Prat[facei]);
-            const scalar yptNew = ypt - f/df;
-
-            if (yptNew < vSmall)
-            {
-                ypsf[facei] = 0;
-                break;
-            }
-            else if (mag(yptNew - ypt) < tolerance_)
-            {
-                ypsf[facei] = yptNew;
-            }
-            else
-            {
-                ypt = yptNew;
-            }
-        }
-
-        ypsf[facei] = ypt;
-    }
-
-    return typsf;
-}
-
-
-// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
-
-alphatPhaseJayatillekeWallFunctionFvPatchScalarField::
-alphatPhaseJayatillekeWallFunctionFvPatchScalarField
-(
-    const fvPatch& p,
-    const DimensionedField<scalar, volMesh>& iF
-)
-:
-    fixedValueFvPatchScalarField(p, iF),
-    Prt_(0.85)
-{}
-
-
-alphatPhaseJayatillekeWallFunctionFvPatchScalarField::
-alphatPhaseJayatillekeWallFunctionFvPatchScalarField
-(
-    const fvPatch& p,
-    const DimensionedField<scalar, volMesh>& iF,
-    const dictionary& dict
-)
-:
-    fixedValueFvPatchScalarField(p, iF, dict),
-    Prt_(dict.lookupOrDefault<scalar>("Prt", 0.85))
-{}
-
-
-alphatPhaseJayatillekeWallFunctionFvPatchScalarField::
-alphatPhaseJayatillekeWallFunctionFvPatchScalarField
-(
-    const alphatPhaseJayatillekeWallFunctionFvPatchScalarField& ptf,
-    const fvPatch& p,
-    const DimensionedField<scalar, volMesh>& iF,
-    const fvPatchFieldMapper& mapper
-)
-:
-    fixedValueFvPatchScalarField(ptf, p, iF, mapper),
-    Prt_(ptf.Prt_)
-{}
-
-
-alphatPhaseJayatillekeWallFunctionFvPatchScalarField::
-alphatPhaseJayatillekeWallFunctionFvPatchScalarField
-(
-    const alphatPhaseJayatillekeWallFunctionFvPatchScalarField& awfpsf,
-    const DimensionedField<scalar, volMesh>& iF
-)
-:
-    fixedValueFvPatchScalarField(awfpsf, iF),
-    Prt_(awfpsf.Prt_)
-{}
-
-
-// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
-
-tmp<scalarField>
-alphatPhaseJayatillekeWallFunctionFvPatchScalarField::calcAlphat
-(
-    const scalarField& prevAlphat
-) const
-{
-    // Lookup the fluid model
-    const phaseSystem& fluid =
-        db().lookupObject<phaseSystem>(phaseSystem::propertiesName);
-
-    const phaseModel& phase
-    (
-        fluid.phases()[internalField().group()]
-    );
-
-    const label patchi = patch().index();
-
-    // Retrieve turbulence properties from model
-    const phaseCompressible::momentumTransportModel& turbModel =
-        db().lookupType<phaseCompressible::momentumTransportModel>
-        (
-            phase.name()
-        );
-
-    const nutWallFunctionFvPatchScalarField& nutw =
-        nutWallFunctionFvPatchScalarField::nutw(turbModel, patchi);
-
-    const scalar Cmu25 = pow025(nutw.Cmu());
-
-    const scalarField& y = turbModel.y()[patchi];
-
-    const tmp<scalarField> tnuw = turbModel.nu(patchi);
-    const scalarField& nuw = tnuw();
-
-    const tmp<scalarField> talphaw
-    (
-        phase.thermo().kappa().boundaryField()[patchi]
-       /phase.thermo().Cp().boundaryField()[patchi]
-    );
-    const scalarField& alphaw = talphaw();
-
-    const tmp<volScalarField> tk = turbModel.k();
-    const volScalarField& k = tk();
-    const fvPatchScalarField& kw = k.boundaryField()[patchi];
-
-    const fvPatchVectorField& Uw = turbModel.U().boundaryField()[patchi];
-    const scalarField magUp(mag(Uw.patchInternalField() - Uw));
-    const scalarField magGradUw(mag(Uw.snGrad()));
-
-    const fvPatchScalarField& rhow = turbModel.rho().boundaryField()[patchi];
-    const fvPatchScalarField& hew =
-        phase.thermo().he().boundaryField()[patchi];
-
-    const fvPatchScalarField& Tw =
-        phase.thermo().T().boundaryField()[patchi];
-
-    const scalarField Tp(Tw.patchInternalField());
-
-    // Heat flux [W/m^2] - lagging alphatw
-    const scalarField qDot
-    (
-        (prevAlphat + alphaw)*hew.snGrad()
-    );
-
-    const scalarField uTau(Cmu25*sqrt(kw));
-
-    const scalarField yPlus(uTau*y/nuw);
-
-    const scalarField Pr(rhow*nuw/alphaw);
-
-    // Molecular-to-turbulent Prandtl number ratio
-    const scalarField Prat(Pr/Prt_);
-
-    // Thermal sublayer thickness
-    const scalarField P(this->Psmooth(Prat));
-
-    const scalarField yPlusTherm(this->yPlusTherm(nutw, P, Prat));
-
-    tmp<scalarField> talphatConv(new scalarField(this->size()));
-    scalarField& alphatConv = talphatConv.ref();
-
-    // Populate boundary values
-    forAll(alphatConv, facei)
-    {
-        // Evaluate new effective thermal diffusivity
-        scalar alphaEff = 0.0;
-        if (yPlus[facei] < yPlusTherm[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
-        {
-            const scalar A = qDot[facei]*rhow[facei]*uTau[facei]*y[facei];
-
-            const scalar B =
-                qDot[facei]*Prt_
-               *(1.0/nutw.kappa()*log(nutw.E()*yPlus[facei]) + P[facei]);
-
-            const scalar magUc =
-                uTau[facei]/nutw.kappa()
-               *log(nutw.E()*yPlusTherm[facei]) - mag(Uw[facei]);
-
-            const scalar C =
-                0.5*rhow[facei]*uTau[facei]
-               *(Prt_*sqr(magUp[facei]) + (Pr[facei] - Prt_)*sqr(magUc));
-
-            alphaEff = A/(B + C + vSmall);
-        }
-
-        // Update convective heat transfer turbulent thermal diffusivity
-        alphatConv[facei] = max(0.0, alphaEff - alphaw[facei]);
-    }
-
-    return talphatConv;
-}
-
-
-void alphatPhaseJayatillekeWallFunctionFvPatchScalarField::updateCoeffs()
-{
-    if (updated())
-    {
-        return;
-    }
-
-    operator==(calcAlphat(*this));
-
-    fixedValueFvPatchScalarField::updateCoeffs();
-}
-
-
-void alphatPhaseJayatillekeWallFunctionFvPatchScalarField::write
-(
-    Ostream& os
-) const
-{
-    fvPatchField<scalar>::write(os);
-    writeEntry(os, "Prt", Prt_);
-    writeEntry(os, "value", *this);
-}
-
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-makePatchTypeField
-(
-    fvPatchScalarField,
-    alphatPhaseJayatillekeWallFunctionFvPatchScalarField
-);
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-} // End namespace compressible
-} // End namespace Foam
-
-// ************************************************************************* //

applications/solvers/modules/multiphaseEuler/multiphaseThermophysicalTransportModels/derivedFvPatchFields/alphatPhaseJayatillekeWallFunction/alphatPhaseJayatillekeWallFunctionFvPatchScalarField.H

@@ -1,203 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2015-2022 OpenFOAM Foundation
-     \\/     M anipulation  |
--------------------------------------------------------------------------------
-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::compressible::
-        alphatPhaseJayatillekeWallFunctionFvPatchScalarField
-
-Description
-    This boundary condition provides a thermal wall function for turbulent
-    thermal diffusivity (usually\c alphat) based on the Jayatilleke model for
-    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>
-    {
-        type            alphatPhaseJayatillekeWallFunction;
-        Prt             0.85;
-        kappa           0.41;
-        E               9.8;
-        value           uniform 0; // optional value entry
-    }
-    \endverbatim
-
-See also
-    Foam::compressible::alphatPhaseChangeWallFunctionFvPatchScalarField
-
-SourceFiles
-    alphatPhaseJayatillekeWallFunctionFvPatchScalarField.C
-
-\*---------------------------------------------------------------------------*/
-
-#ifndef compressible_alphatPhaseJayatillekeWallFunctionFvPatchScalarField_H
-#define compressible_alphatPhaseJayatillekeWallFunctionFvPatchScalarField_H
-
-#include "fixedValueFvPatchFields.H"
-#include "nutWallFunctionFvPatchScalarField.H"
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-namespace Foam
-{
-namespace compressible
-{
-
-/*---------------------------------------------------------------------------*\
-    Class alphatPhaseJayatillekeWallFunctionFvPatchScalarField Declaration
-\*---------------------------------------------------------------------------*/
-
-class alphatPhaseJayatillekeWallFunctionFvPatchScalarField
-:
-    public fixedValueFvPatchScalarField
-{
-
-protected:
-
-    // Protected data
-
-        //- Turbulent Prandtl number
-        scalar Prt_;
-
-        // Solution parameters
-
-            static scalar maxExp_;
-            static scalar tolerance_;
-            static label maxIters_;
-
-
-    // Protected Member Functions
-
-        //- 'P' function
-        tmp<scalarField> Psmooth(const scalarField& Prat) const;
-
-        //- Calculate y+ at the edge of the thermal laminar sublayer
-        tmp<scalarField> yPlusTherm
-        (
-            const nutWallFunctionFvPatchScalarField& nutw,
-            const scalarField& P,
-            const scalarField& Prat
-        ) const;
-
-
-public:
-
-    //- Runtime type information
-    TypeName("compressible::alphatPhaseJayatillekeWallFunction");
-
-
-    // Constructors
-
-        //- Construct from patch and internal field
-        alphatPhaseJayatillekeWallFunctionFvPatchScalarField
-        (
-            const fvPatch&,
-            const DimensionedField<scalar, volMesh>&
-        );
-
-        //- Construct from patch, internal field and dictionary
-        alphatPhaseJayatillekeWallFunctionFvPatchScalarField
-        (
-            const fvPatch&,
-            const DimensionedField<scalar, volMesh>&,
-            const dictionary&
-        );
-
-        //- Construct by mapping given
-        //  alphatPhaseJayatillekeWallFunctionFvPatchScalarField
-        //  onto a new patch
-        alphatPhaseJayatillekeWallFunctionFvPatchScalarField
-        (
-            const alphatPhaseJayatillekeWallFunctionFvPatchScalarField&,
-            const fvPatch&,
-            const DimensionedField<scalar, volMesh>&,
-            const fvPatchFieldMapper&
-        );
-
-        //- Disallow copy without setting internal field reference
-        alphatPhaseJayatillekeWallFunctionFvPatchScalarField
-        (
-            const alphatPhaseJayatillekeWallFunctionFvPatchScalarField&
-        ) = delete;
-
-        //- Copy constructor setting internal field reference
-        alphatPhaseJayatillekeWallFunctionFvPatchScalarField
-        (
-            const alphatPhaseJayatillekeWallFunctionFvPatchScalarField&,
-            const DimensionedField<scalar, volMesh>&
-        );
-
-        //- Construct and return a clone setting internal field reference
-        virtual tmp<fvPatchScalarField> clone
-        (
-            const DimensionedField<scalar, volMesh>& iF
-        ) const
-        {
-            return tmp<fvPatchScalarField>
-            (
-                new alphatPhaseJayatillekeWallFunctionFvPatchScalarField
-                (
-                    *this,
-                    iF
-                )
-            );
-        }
-
-
-    // Member Functions
-
-        // Evaluation functions
-
-            //- Evaluate the turbulent thermal diffusivity
-            tmp<scalarField> calcAlphat(const scalarField& prevAlphat) const;
-
-            //- Update the coefficients associated with the patch field
-            virtual void updateCoeffs();
-
-
-        // I-O
-
-            //- Write
-            virtual void write(Ostream&) const;
-};
-
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-} // End namespace compressible
-} // End namespace Foam
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-#endif
-
-// ************************************************************************* //

applications/solvers/modules/multiphaseEuler/multiphaseThermophysicalTransportModels/derivedFvPatchFields/alphatWallBoilingWallFunction/alphatWallBoilingWallFunctionFvPatchScalarField.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2015-2022 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2015-2023 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -24,6 +24,8 @@ License
 \*---------------------------------------------------------------------------*/
 
 #include "alphatWallBoilingWallFunctionFvPatchScalarField.H"
+#include "alphatJayatillekeWallFunctionFvPatchScalarField.H"
+#include "fluidThermophysicalTransportModel.H"
 #include "phaseSystem.H"
 #include "heatTransferPhaseSystem.H"
 #include "compressibleMomentumTransportModels.H"
@@ -39,23 +41,19 @@ using namespace Foam::constant::mathematical;
 template<>
 const char* Foam::NamedEnum
 <
-    Foam::compressible::
+    Foam::compressible::alphatWallBoilingWallFunctionFvPatchScalarField::
-    alphatWallBoilingWallFunctionFvPatchScalarField::phaseType,
+        phaseType,
     2
->::names[] =
+>::names[] = {"vapor", "liquid"};
-{
-    "vapor",
-    "liquid"
-};
 
 const Foam::NamedEnum
 <
-    Foam::compressible::
+    Foam::compressible::alphatWallBoilingWallFunctionFvPatchScalarField::
-    alphatWallBoilingWallFunctionFvPatchScalarField::phaseType,
+        phaseType,
     2
 >
-Foam::compressible::
+Foam::compressible::alphatWallBoilingWallFunctionFvPatchScalarField::
-alphatWallBoilingWallFunctionFvPatchScalarField::phaseTypeNames_;
+    phaseTypeNames_;
 
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@@ -76,6 +74,7 @@ alphatWallBoilingWallFunctionFvPatchScalarField
 :
     alphatPhaseChangeWallFunctionFvPatchScalarField(p, iF),
     phaseType_(liquidPhase),
+    Prt_(0.85),
     AbyV_(p.size(), 0),
     alphatConv_(p.size(), 0),
     dDep_(p.size(), 1e-5),
@@ -109,6 +108,7 @@ alphatWallBoilingWallFunctionFvPatchScalarField
 :
     alphatPhaseChangeWallFunctionFvPatchScalarField(p, iF, dict),
     phaseType_(phaseTypeNames_.read(dict.lookup("phaseType"))),
+    Prt_(dict.lookupOrDefault<scalar>("Prt", 0.85)),
     AbyV_(p.size(), 0),
     alphatConv_(p.size(), 0),
     dDep_(p.size(), 1e-5),
@@ -231,14 +231,9 @@ alphatWallBoilingWallFunctionFvPatchScalarField
     const fvPatchFieldMapper& mapper
 )
 :
-    alphatPhaseChangeWallFunctionFvPatchScalarField
+    alphatPhaseChangeWallFunctionFvPatchScalarField(psf, p, iF, mapper),
-    (
-        psf,
-        p,
-        iF,
-        mapper
-    ),
     phaseType_(psf.phaseType_),
+    Prt_(psf.Prt_),
     AbyV_(mapper(psf.AbyV_)),
     alphatConv_(mapper(psf.alphatConv_)),
     dDep_(mapper(psf.dDep_)),
@@ -264,6 +259,7 @@ alphatWallBoilingWallFunctionFvPatchScalarField
 :
     alphatPhaseChangeWallFunctionFvPatchScalarField(psf, iF),
     phaseType_(psf.phaseType_),
+    Prt_(psf.Prt_),
     AbyV_(psf.AbyV_),
     alphatConv_(psf.alphatConv_),
     dDep_(psf.dDep_),
@@ -364,20 +360,36 @@ void alphatWallBoilingWallFunctionFvPatchScalarField::updateCoeffs()
         {
             const phaseModel& vapor = fluid.phases()[internalField().group()];
 
+            // Vapor thermophysical transport model
+            const fluidThermophysicalTransportModel& vaporTtm =
+                db().lookupType<fluidThermophysicalTransportModel>
+                (
+                    vapor.name()
+                );
+
             // Vapor phase fraction at the wall
             const scalarField& vaporw = vapor.boundaryField()[patchi];
 
-            // Partitioning
+            // Partitioning. Note: Assumes that there is only only one liquid
-            // NOTE! Assumes that there is only only one liquid phase and all
+            // phase and all other phases are vapor.
-            // other phases are vapor
-
             const phaseModel& liquid = fluid.phases()[otherPhaseName_];
             const scalarField& liquidw = liquid.boundaryField()[patchi];
             fLiquid_ = partitioningModel_->fLiquid(liquidw);
 
+            // Vapour thermal diffusivity
+            const scalarField alphatConv
+            (
+                alphatJayatillekeWallFunctionFvPatchScalarField::alphat
+                (
+                    vaporTtm,
+                    Prt_,
+                    patch().index()
+                )
+            );
+
             operator==
             (
-                calcAlphat(*this)*(vaporw/(1 - liquidw + small) )
+                alphatConv*(vaporw/(1 - liquidw + small) )
                *(1 - fLiquid_)/max(vaporw, scalar(1e-8))
             );
             break;
@@ -394,60 +406,58 @@ void alphatWallBoilingWallFunctionFvPatchScalarField::updateCoeffs()
                 fluid.foundInterfacialModel
                 <
                     interfaceSaturationTemperatureModel
-                >(interface)
+                >
+                (interface)
             )
             {
-                // Retrieve turbulence properties from models
+                // Liquid thermophysical and momentum transport models
-                const phaseCompressible::momentumTransportModel& turbModel
+                const fluidThermophysicalTransportModel& liquidTtm =
-                  = db().lookupType<phaseCompressible::momentumTransportModel>
+                    db().lookupType<fluidThermophysicalTransportModel>
                     (
                         liquid.name()
                     );
-                const phaseCompressible::momentumTransportModel& vaporTurbModel
+                const compressibleMomentumTransportModel& liquidMtm =
-                  = db().lookupType<phaseCompressible::momentumTransportModel>
+                    liquidTtm.momentumTransport();
-                    (
-                        vapor.name()
-                    );
 
                 const nutWallFunctionFvPatchScalarField& nutw =
-                    nutWallFunctionFvPatchScalarField::nutw(turbModel, patchi);
+                    nutWallFunctionFvPatchScalarField::nutw(liquidMtm, patchi);
 
                 const scalar Cmu25(pow025(nutw.Cmu()));
 
-                const scalarField& y = turbModel.y()[patchi];
+                const scalarField& y = liquidMtm.y()[patchi];
 
-                const tmp<scalarField> tnuw = turbModel.nu(patchi);
+                const tmp<scalarField> tnuw = liquidMtm.nu(patchi);
                 const scalarField& nuw = tnuw();
 
-                const rhoThermo& lThermo = liquid.thermo();
-
                 const tmp<scalarField> talphaw
                 (
-                    lThermo.kappa().boundaryField()[patchi]
+                    liquid.thermo().kappa().boundaryField()[patchi]
-                   /lThermo.Cp().boundaryField()[patchi]
+                   /liquid.thermo().Cp().boundaryField()[patchi]
                 );
                 const scalarField& alphaw = talphaw();
 
-                const tmp<volScalarField> tk = turbModel.k();
+                const tmp<volScalarField> tk = liquidMtm.k();
                 const volScalarField& k = tk();
                 const fvPatchScalarField& kw = k.boundaryField()[patchi];
 
                 const fvPatchVectorField& Uw =
-                    turbModel.U().boundaryField()[patchi];
+                    liquidMtm.U().boundaryField()[patchi];
                 const scalarField magUp(mag(Uw.patchInternalField() - Uw));
                 const scalarField magGradUw(mag(Uw.snGrad()));
 
-                const fvPatchScalarField& rhoLiquidw =
+                const tmp<scalarField> trhoLiquidw =
-                    turbModel.rho().boundaryField()[patchi];
+                    liquid.thermo().rho(patchi);
+                const scalarField rhoLiquidw = trhoLiquidw();
 
-                const fvPatchScalarField& rhoVaporw =
+                const tmp<scalarField> trhoVaporw =
-                    vaporTurbModel.rho().boundaryField()[patchi];
+                    vapor.thermo().rho(patchi);
+                const scalarField rhoVaporw = trhoVaporw();
 
                 const fvPatchScalarField& hew =
-                    lThermo.he().boundaryField()[patchi];
+                    liquid.thermo().he().boundaryField()[patchi];
 
                 const fvPatchScalarField& Tw =
-                    lThermo.T().boundaryField()[patchi];
+                    liquid.thermo().T().boundaryField()[patchi];
 
                 const scalarField Tc(Tw.patchInternalField());
 
@@ -461,11 +471,21 @@ void alphatWallBoilingWallFunctionFvPatchScalarField::updateCoeffs()
                 const scalarField Prat(Pr/Prt_);
 
                 // Thermal sublayer thickness
-                const scalarField P(this->Psmooth(Prat));
+                const scalarField P
+                (
+                    alphatJayatillekeWallFunctionFvPatchScalarField::P(Prat)
+                );
+                const scalarField yPlusTherm
+                (
+                    alphatJayatillekeWallFunctionFvPatchScalarField::yPlusTherm
+                    (
+                        nutw,
+                        P,
+                        Prat
+                    )
+                );
 
-                const scalarField yPlusTherm(this->yPlusTherm(nutw, P, Prat));
+                const scalarField Cpw(liquid.thermo().Cp(Tw, patchi));
-
-                const scalarField Cpw(lThermo.Cp(Tw, patchi));
 
                 // Saturation temperature
                 const interfaceSaturationTemperatureModel& satModel =
@@ -473,7 +493,8 @@ void alphatWallBoilingWallFunctionFvPatchScalarField::updateCoeffs()
                     <
                         interfaceSaturationTemperatureModel
                     >(interface);
-                const tmp<volScalarField> tTsat = satModel.Tsat(lThermo.p());
+                const tmp<volScalarField> tTsat =
+                    satModel.Tsat(liquid.thermo().p());
                 const volScalarField& Tsat = tTsat();
                 const fvPatchScalarField& Tsatw(Tsat.boundaryField()[patchi]);
 
@@ -493,7 +514,7 @@ void alphatWallBoilingWallFunctionFvPatchScalarField::updateCoeffs()
                      )
                   : -refCast<const heatTransferPhaseSystem>(fluid)
                     .L
-                     (
+                    (
                          interface,
                          dmdtf_,
                          Tsat,
@@ -509,7 +530,13 @@ void alphatWallBoilingWallFunctionFvPatchScalarField::updateCoeffs()
                 fLiquid_ = partitioningModel_->fLiquid(liquidw);
 
                 // Convective thermal diffusivity
-                alphatConv_ = calcAlphat(alphatConv_);
+                alphatConv_ =
+                    alphatJayatillekeWallFunctionFvPatchScalarField::alphat
+                    (
+                        liquidTtm,
+                        Prt_,
+                        patch().index()
+                    );
 
                 label maxIter(10);
                 for (label i=0; i<maxIter; i++)
@@ -787,6 +814,7 @@ makePatchTypeField
     alphatWallBoilingWallFunctionFvPatchScalarField
 );
 
+
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 } // End namespace compressible

applications/solvers/modules/multiphaseEuler/multiphaseThermophysicalTransportModels/derivedFvPatchFields/alphatWallBoilingWallFunction/alphatWallBoilingWallFunctionFvPatchScalarField.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2015-2022 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2015-2023 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -66,29 +66,25 @@ Description
 
 Usage
     \table
-        Property     | Description             | Required    | Default value
+        Property       | Description              | Required    | Default value
-        phaseType    | 'vapor' or 'liquid'     | yes         |
+        phaseType      | 'vapor' or 'liquid'      | yes         |
-        relax        | wall boiling model relaxation| yes     |
+        relax          | relaxation factor        | no          | 1
-        Prt   | inherited from alphatPhaseChangeJayatillekeWallFunction | yes |
+        dmdt           | phase change mass flux   | no          | uniform 0
-        Cmu   | inherited from alphatPhaseChangeJayatillekeWallFunction | yes |
+        Prt            | turbulent Prandtl number | no          | 0.85
-        kappa | inherited from alphatPhaseChangeJayatillekeWallFunction | yes |
-        E     | inherited from alphatPhaseChangeJayatillekeWallFunction | yes |
-        dmdt         | phase change mass flux  | yes         |
-        value        | initial alphat value    | yes         |
     \endtable
 
         if phaseType 'vapor':
     \table
-        partitioningModel|                     | yes         |
+        partitioningModel |                       | yes         |
     \endtable
 
         if phaseType 'liquid':
     \table
-        partitioningModel|                     | yes         |
+        partitioningModel |                       | yes         |
-        nucleationSiteModel|                   | yes         |
+        nucleationSiteModel |                     | yes         |
-        departureDiamModel|                    | yes         |
+        departureDiamModel |                      | yes         |
-        departureFreqModel|                    | yes         |
+        departureFreqModel |                      | yes         |
-        bubbleWaitingTime|                     | no          | yes
+        bubbleWaitingTimeRatio |                  | no          | 0.8
     \endtable
 
     NOTE: Runtime selectable submodels may require model specific entries
@@ -99,12 +95,9 @@ Usage
     {
         type            compressible::alphatWallBoiling2WallFunction;
         phaseType       liquid;
-        Prt             0.85;
-        Cmu             0.09;
-        kappa           0.41;
-        E               9.8;
         relax           0.1;
         dmdt            uniform 0;
+        Prt             0.85;
         partitioningModel
         {
             type        Lavieville;
@@ -125,9 +118,6 @@ Usage
         value           uniform 0.01;
     \endverbatim
 
-See also
-    Foam::alphatPhaseChangeJayatillekeWallFunctionFvPatchField
-
 SourceFiles
     alphatWallBoilingWallFunctionFvPatchScalarField.C
 
@@ -179,6 +169,9 @@ private:
         //- Heat source type
         phaseType phaseType_;
 
+        //- Turbulent Prandtl number
+        scalar Prt_;
+
         //- Patch face area by cell volume
         scalarField AbyV_;
 

src/ThermophysicalTransportModels/fluid/derivedFvPatchFields/alphatWallFunctions/alphatJayatillekeWallFunction/alphatJayatillekeWallFunctionFvPatchScalarField.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2011-2022 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2023 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -40,49 +40,6 @@ scalar alphatJayatillekeWallFunctionFvPatchScalarField::maxExp_ = 50.0;
 scalar alphatJayatillekeWallFunctionFvPatchScalarField::tolerance_ = 0.01;
 label alphatJayatillekeWallFunctionFvPatchScalarField::maxIters_ = 10;
 
-// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
-
-scalar alphatJayatillekeWallFunctionFvPatchScalarField::Psmooth
-(
-    const scalar Prat
-) const
-{
-    return 9.24*(pow(Prat, 0.75) - 1.0)*(1.0 + 0.28*exp(-0.007*Prat));
-}
-
-
-scalar alphatJayatillekeWallFunctionFvPatchScalarField::yPlusTherm
-(
-    const nutWallFunctionFvPatchScalarField& nutw,
-    const scalar P,
-    const scalar Prat
-) const
-{
-    scalar ypt = 11.0;
-
-    for (int i=0; i<maxIters_; i++)
-    {
-        scalar f = ypt - (log(nutw.E()*ypt)/nutw.kappa() + P)/Prat;
-        scalar df = 1.0 - 1.0/(ypt*nutw.kappa()*Prat);
-        scalar yptNew = ypt - f/df;
-
-        if (yptNew < vSmall)
-        {
-            return 0;
-        }
-        else if (mag(yptNew - ypt) < tolerance_)
-        {
-            return yptNew;
-        }
-        else
-        {
-            ypt = yptNew;
-        }
-     }
-
-    return ypt;
-}
-
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
@@ -139,21 +96,67 @@ alphatJayatillekeWallFunctionFvPatchScalarField
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
-void alphatJayatillekeWallFunctionFvPatchScalarField::updateCoeffs()
+tmp<scalarField> alphatJayatillekeWallFunctionFvPatchScalarField::P
+(
+    const scalarField& Prat
+)
 {
-    if (updated())
+    return 9.24*(pow(Prat, 0.75) - 1.0)*(1.0 + 0.28*exp(-0.007*Prat));
+}
+
+
+tmp<scalarField> alphatJayatillekeWallFunctionFvPatchScalarField::yPlusTherm
+(
+    const nutWallFunctionFvPatchScalarField& nutw,
+    const scalarField& P,
+    const scalarField& Prat
+)
+{
+    tmp<scalarField> typt(new scalarField(nutw.size()));
+    scalarField& ypt = typt.ref();
+
+    const scalar E = nutw.E();
+    const scalar kappa = nutw.kappa();
+
+    forAll(ypt, facei)
     {
-        return;
+        ypt[facei] = 11.0;
+
+        for (int i=0; i<maxIters_; i++)
+        {
+            const scalar f =
+                ypt[facei] - (log(E*ypt[facei])/kappa + P[facei])/Prat[facei];
+
+            const scalar df = 1.0 - 1.0/(ypt[facei]*nutw.kappa()*Prat[facei]);
+
+            const scalar dypt = - f/df;
+
+            ypt[facei] += dypt;
+
+            if (ypt[facei] < vSmall)
+            {
+                ypt[facei] = 0;
+                break;
+            }
+
+            if (mag(dypt) < tolerance_)
+            {
+                break;
+            }
+        }
     }
 
-    const label patchi = patch().index();
+    return typt;
+}
 
-    const fluidThermophysicalTransportModel& ttm =
-        db().lookupType<fluidThermophysicalTransportModel>
-        (
-            internalField().group()
-        );
 
+tmp<scalarField> alphatJayatillekeWallFunctionFvPatchScalarField::alphat
+(
+    const fluidThermophysicalTransportModel& ttm,
+    const scalar Prt,
+    const label patchi
+)
+{
     const compressibleMomentumTransportModel& turbModel =
         ttm.momentumTransport();
 
@@ -174,8 +177,6 @@ void alphatJayatillekeWallFunctionFvPatchScalarField::updateCoeffs()
     );
     const scalarField& alphaw = talphaw();
 
-    scalarField& alphatw = *this;
-
     const tmp<volScalarField> tk = turbModel.k();
     const volScalarField& k = tk();
 
@@ -189,34 +190,47 @@ void alphatJayatillekeWallFunctionFvPatchScalarField::updateCoeffs()
     // Enthalpy gradient
     const scalarField gradHew(hew.snGrad());
 
+    // Molecular Prandtl number
+    const scalarField Pr(rhow*nuw/alphaw);
+
+    // Molecular-to-turbulent Prandtl number ratio
+    const scalarField Prat(Pr/Prt);
+
+    // Thermal sublayer thickness
+    const scalarField P
+    (
+        alphatJayatillekeWallFunctionFvPatchScalarField::P(Prat)
+    );
+    const scalarField yPlusTherm
+    (
+        alphatJayatillekeWallFunctionFvPatchScalarField::yPlusTherm
+        (
+            nutw,
+            P,
+            Prat
+        )
+    );
+
     // Populate boundary values
+    tmp<scalarField> talphatw(new scalarField(nutw.size()));
+    scalarField& alphatw = talphatw.ref();
     forAll(alphatw, facei)
     {
-        label celli = patch().faceCells()[facei];
+        const label celli = nutw.patch().faceCells()[facei];
 
-        scalar uTau = Cmu25*sqrt(k[celli]);
+        const scalar uTau = Cmu25*sqrt(k[celli]);
 
-        scalar yPlus = uTau*y[facei]/nuw[facei];
+        const scalar yPlus = uTau*y[facei]/nuw[facei];
-
-        // Molecular Prandtl number
-        scalar Pr = rhow[facei]*nuw[facei]/alphaw[facei];
-
-        // Molecular-to-turbulent Prandtl number ratio
-        scalar Prat = Pr/Prt_;
-
-        // Thermal sublayer thickness
-        scalar P = Psmooth(Prat);
-        scalar yPlusTherm = this->yPlusTherm(nutw, P, Prat);
 
         // Evaluate new effective thermal diffusivity
         scalar alphaEff = 0.0;
-        if (yPlus < yPlusTherm)
+        if (yPlus < yPlusTherm[facei])
         {
             const scalar A = gradHew[facei]*rhow[facei]*uTau*y[facei];
 
-            const scalar B = gradHew[facei]*Pr*yPlus;
+            const scalar B = gradHew[facei]*Pr[facei]*yPlus;
 
-            const scalar C = Pr*0.5*rhow[facei]*uTau*sqr(magUp[facei]);
+            const scalar C = Pr[facei]*0.5*rhow[facei]*uTau*sqr(magUp[facei]);
 
             alphaEff = (A - C)/(B + sign(B)*rootVSmall);
         }
@@ -225,14 +239,16 @@ void alphatJayatillekeWallFunctionFvPatchScalarField::updateCoeffs()
             const scalar A = gradHew[facei]*rhow[facei]*uTau*y[facei];
 
             const scalar B =
-                gradHew[facei]*Prt_*(1.0/nutw.kappa()*log(nutw.E()*yPlus) + P);
+                gradHew[facei]*Prt
+               *(1.0/nutw.kappa()*log(nutw.E()*yPlus) + P[facei]);
 
             const scalar magUc =
-                uTau/nutw.kappa()*log(nutw.E()*yPlusTherm) - mag(Uw[facei]);
+                uTau/nutw.kappa()
+               *log(nutw.E()*yPlusTherm[facei]) - mag(Uw[facei]);
 
             const scalar C =
                 0.5*rhow[facei]*uTau
-               *(Prt_*sqr(magUp[facei]) + (Pr - Prt_)*sqr(magUc));
+               *(Prt*sqr(magUp[facei]) + (Pr[facei] - Prt)*sqr(magUc));
 
             alphaEff = (A - C)/(B + sign(B)*rootVSmall);
         }
@@ -246,6 +262,25 @@ void alphatJayatillekeWallFunctionFvPatchScalarField::updateCoeffs()
             min(max(alphaEff - alphaw[facei], alphatwMin), alphatwMax);
     }
 
+    return talphatw;
+}
+
+
+void alphatJayatillekeWallFunctionFvPatchScalarField::updateCoeffs()
+{
+    if (updated())
+    {
+        return;
+    }
+
+    const fluidThermophysicalTransportModel& ttm =
+        db().lookupType<fluidThermophysicalTransportModel>
+        (
+            internalField().group()
+        );
+
+    this->operator==(alphat(ttm, Prt_, patch().index()));
+
     fixedValueFvPatchField<scalar>::updateCoeffs();
 }
 

src/ThermophysicalTransportModels/fluid/derivedFvPatchFields/alphatWallFunctions/alphatJayatillekeWallFunction/alphatJayatillekeWallFunctionFvPatchScalarField.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2011-2022 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2023 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -32,9 +32,6 @@ 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:
@@ -43,12 +40,13 @@ Usage
     {
         type            alphatJayatillekeWallFunction;
         Prt             0.85;
-        kappa           0.41;
+        value           uniform 0;
-        E               9.8;
-        value           uniform 0; // optional value entry
     }
     \endverbatim
 
+    Note that other model constants (i.e., Cmu, kappa and E) are obtained from
+    the corresponding turbulent viscosity boundary condition.
+
 See also
     Foam::fixedValueFvPatchField
 
@@ -67,6 +65,9 @@ SourceFiles
 
 namespace Foam
 {
+
+class fluidThermophysicalTransportModel;
+
 namespace compressible
 {
 
@@ -91,20 +92,6 @@ class alphatJayatillekeWallFunctionFvPatchScalarField
             static label maxIters_;
 
 
-    // Private Member Functions
-
-        //- `P' function
-        scalar Psmooth(const scalar Prat) const;
-
-        //- Calculate y+ at the edge of the thermal laminar sublayer
-        scalar yPlusTherm
-        (
-            const nutWallFunctionFvPatchScalarField& nutw,
-            const scalar P,
-            const scalar Prat
-        ) const;
-
-
 public:
 
     //- Runtime type information
@@ -173,6 +160,25 @@ public:
 
         // Evaluation functions
 
+            //- Calculate the smoothing function
+            static tmp<scalarField> P(const scalarField& Prat);
+
+            //- Calculate y+ at the edge of the thermal laminar sublayer
+            static tmp<scalarField> yPlusTherm
+            (
+                const nutWallFunctionFvPatchScalarField& nutw,
+                const scalarField& P,
+                const scalarField& Prat
+            );
+
+            //- Calculate the turbulent thermal diffusivity
+            static tmp<scalarField> alphat
+            (
+                const fluidThermophysicalTransportModel& ttm,
+                const scalar Prt,
+                const label patchi
+            );
+
             //- Update the coefficients associated with the patch field
             virtual void updateCoeffs();
 

tutorials/modules/multiphaseEuler/Grossetete/0/alphat.gas

@@ -32,11 +32,8 @@ boundaryField
     }
     walls
     {
-        type            compressible::alphatPhaseJayatillekeWallFunction;
+        type            compressible::alphatJayatillekeWallFunction;
         Prt             0.85;
-        Cmu             0.09;
-        kappa           0.41;
-        E               9.8;
         value           uniform 0;
     }
     front

tutorials/modules/multiphaseEuler/Grossetete/0/alphat.liquid

@@ -32,11 +32,8 @@ boundaryField
     }
     walls
     {
-        type            compressible::alphatPhaseJayatillekeWallFunction;
+        type            compressible::alphatJayatillekeWallFunction;
         Prt             0.85;
-        Cmu             0.09;
-        kappa           0.41;
-        E               9.8;
         value           uniform 0;
     }
     front

tutorials/modules/multiphaseEuler/hydrofoil/0/alphat.gas

@@ -35,20 +35,14 @@ boundaryField
     }
     hydrofoil
     {
-        type            compressible::alphatPhaseJayatillekeWallFunction;
+        type            compressible::alphatJayatillekeWallFunction;
         Prt             0.85;
-        Cmu             0.09;
-        kappa           0.41;
-        E               9.8;
         value           uniform 0;
     }
     walls
     {
-        type            compressible::alphatPhaseJayatillekeWallFunction;
+        type            compressible::alphatJayatillekeWallFunction;
         Prt             0.85;
-        Cmu             0.09;
-        kappa           0.41;
-        E               9.8;
         value           uniform 0;
     }
 }

tutorials/modules/multiphaseEuler/hydrofoil/0/alphat.liquid

@@ -35,20 +35,14 @@ boundaryField
     }
     hydrofoil
     {
-        type            compressible::alphatPhaseJayatillekeWallFunction;
+        type            compressible::alphatJayatillekeWallFunction;
         Prt             0.85;
-        Cmu             0.09;
-        kappa           0.41;
-        E               9.8;
         value           uniform 0;
     }
     walls
     {
-        type            compressible::alphatPhaseJayatillekeWallFunction;
+        type            compressible::alphatJayatillekeWallFunction;
         Prt             0.85;
-        Cmu             0.09;
-        kappa           0.41;
-        E               9.8;
         value           uniform 0;
     }
 }