diff --git a/applications/solvers/modules/multiphaseEuler/functionObjects/wallBoilingProperties/wallBoilingProperties.C b/applications/solvers/modules/multiphaseEuler/functionObjects/wallBoilingProperties/wallBoilingProperties.C
index 635b459766..152c5e541a 100644
--- a/applications/solvers/modules/multiphaseEuler/functionObjects/wallBoilingProperties/wallBoilingProperties.C
+++ b/applications/solvers/modules/multiphaseEuler/functionObjects/wallBoilingProperties/wallBoilingProperties.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2022 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2022-2023 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -165,10 +165,8 @@ bool Foam::functionObjects::wallBoilingProperties::write()
                 if (isA<alphatWallBoilingWallFunction>(alphatBf[patchi]))
                 {
                     const alphatWallBoilingWallFunction& alphatw =
-                        refCast
-                        <
-                            const alphatWallBoilingWallFunction
-                        >(alphatBf[patchi]);
+                        refCast<const alphatWallBoilingWallFunction>
+                        (alphatBf[patchi]);
 
                     dDeparture.boundaryFieldRef()[patchi] =
                         alphatw.dDeparture();
diff --git a/applications/solvers/modules/multiphaseEuler/multiphaseThermophysicalTransportModels/Make/files b/applications/solvers/modules/multiphaseEuler/multiphaseThermophysicalTransportModels/Make/files
index 01bf80bb84..301e729592 100644
--- a/applications/solvers/modules/multiphaseEuler/multiphaseThermophysicalTransportModels/Make/files
+++ b/applications/solvers/modules/multiphaseEuler/multiphaseThermophysicalTransportModels/Make/files
@@ -22,7 +22,7 @@ derivedFvPatchFields/wallBoilingSubModels/departureFrequencyModels/Kocamustafaog
 
 derivedFvPatchFields/alphatPhaseChangeWallFunctionBase/alphatPhaseChangeWallFunctionBase.C
 derivedFvPatchFields/alphatWallBoilingWallFunction/alphatWallBoilingWallFunctionFvPatchScalarField.C
-derivedFvPatchFields/fixedMultiphaseHeatFlux/fixedMultiphaseHeatFluxFvPatchScalarField.C
+derivedFvPatchFields/uniformFixedMultiphaseHeatFlux/uniformFixedMultiphaseHeatFluxFvPatchScalarField.C
 derivedFvPatchFields/coupledMultiphaseTemperature/coupledMultiphaseTemperatureFvPatchScalarField.C
 
 LIB = $(FOAM_LIBBIN)/libmultiphaseThermophysicalTransportModels
diff --git a/applications/solvers/modules/multiphaseEuler/multiphaseThermophysicalTransportModels/derivedFvPatchFields/alphatPhaseChangeWallFunctionBase/alphatPhaseChangeWallFunctionBase.C b/applications/solvers/modules/multiphaseEuler/multiphaseThermophysicalTransportModels/derivedFvPatchFields/alphatPhaseChangeWallFunctionBase/alphatPhaseChangeWallFunctionBase.C
index fb0ce53a43..b1ca3104ca 100644
--- a/applications/solvers/modules/multiphaseEuler/multiphaseThermophysicalTransportModels/derivedFvPatchFields/alphatPhaseChangeWallFunctionBase/alphatPhaseChangeWallFunctionBase.C
+++ b/applications/solvers/modules/multiphaseEuler/multiphaseThermophysicalTransportModels/derivedFvPatchFields/alphatPhaseChangeWallFunctionBase/alphatPhaseChangeWallFunctionBase.C
@@ -94,23 +94,6 @@ bool alphatPhaseChangeWallFunctionBase::activeInterface
 }
 
 
-const scalarField& alphatPhaseChangeWallFunctionBase::dmdtf
-(
-    const phaseInterface& interface
-) const
-{
-    if (!activeInterface(interface))
-    {
-        FatalErrorInFunction
-            << "Phase change mass transfer rate requested for interface on "
-            << "which there is no phase change "
-            << abort(FatalError);
-    }
-
-    return dmdtf();
-}
-
-
 void alphatPhaseChangeWallFunctionBase::write(Ostream& os) const
 {
     writeEntry(os, "otherPhase", otherPhaseName_);
diff --git a/applications/solvers/modules/multiphaseEuler/multiphaseThermophysicalTransportModels/derivedFvPatchFields/alphatPhaseChangeWallFunctionBase/alphatPhaseChangeWallFunctionBase.H b/applications/solvers/modules/multiphaseEuler/multiphaseThermophysicalTransportModels/derivedFvPatchFields/alphatPhaseChangeWallFunctionBase/alphatPhaseChangeWallFunctionBase.H
index c0beeaf51f..3fcf98a142 100644
--- a/applications/solvers/modules/multiphaseEuler/multiphaseThermophysicalTransportModels/derivedFvPatchFields/alphatPhaseChangeWallFunctionBase/alphatPhaseChangeWallFunctionBase.H
+++ b/applications/solvers/modules/multiphaseEuler/multiphaseThermophysicalTransportModels/derivedFvPatchFields/alphatPhaseChangeWallFunctionBase/alphatPhaseChangeWallFunctionBase.H
@@ -96,9 +96,6 @@ public:
         //- Return the rate of phase-change
         virtual const scalarField& dmdtf() const = 0;
 
-        //- Return the rate of phase-change for the given interface
-        const scalarField& dmdtf(const phaseInterface&) const;
-
         //- Write
         virtual void write(Ostream&) const;
 };
diff --git a/applications/solvers/modules/multiphaseEuler/multiphaseThermophysicalTransportModels/derivedFvPatchFields/alphatWallBoilingWallFunction/alphatWallBoilingWallFunctionFvPatchScalarField.C b/applications/solvers/modules/multiphaseEuler/multiphaseThermophysicalTransportModels/derivedFvPatchFields/alphatWallBoilingWallFunction/alphatWallBoilingWallFunctionFvPatchScalarField.C
index c5b982092e..c2f34f7e54 100644
--- a/applications/solvers/modules/multiphaseEuler/multiphaseThermophysicalTransportModels/derivedFvPatchFields/alphatWallBoilingWallFunction/alphatWallBoilingWallFunctionFvPatchScalarField.C
+++ b/applications/solvers/modules/multiphaseEuler/multiphaseThermophysicalTransportModels/derivedFvPatchFields/alphatWallBoilingWallFunction/alphatWallBoilingWallFunctionFvPatchScalarField.C
@@ -32,6 +32,12 @@ License
 #include "phaseCompressibleMomentumTransportModel.H"
 #include "interfaceSaturationTemperatureModel.H"
 #include "rhoMulticomponentThermo.H"
+
+#include "fixedValueFvPatchFields.H"
+#include "zeroGradientFvPatchFields.H"
+#include "fixedGradientFvPatchFields.H"
+#include "mixedFvPatchFields.H"
+
 #include "addToRunTimeSelectionTable.H"
 
 using namespace Foam::constant::mathematical;
@@ -63,6 +69,464 @@ namespace Foam
 namespace compressible
 {
 
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+struct alphatWallBoilingWallFunctionFvPatchScalarField::properties
+{
+    // Data
+
+        //- Wall function field
+        const alphatWallBoilingWallFunctionFvPatchScalarField& field;
+
+        //- Phase
+        const phaseModel& phase;
+
+        //- Other phase
+        const phaseModel& otherPhase;
+
+        //- Volume fraction
+        const scalarField& alphaw;
+
+        //- Other volume fraction
+        const scalarField& otherAlphaw;
+
+        //- Interface
+        const phaseInterface interface;
+
+        //- Phase thermophysical transport model
+        const fluidThermophysicalTransportModel& ttm;
+
+        //- Phase convective turbulent thermal diffusivity
+        const scalarField alphatConv;
+
+
+    //- Constructor
+    properties
+    (
+        const alphatWallBoilingWallFunctionFvPatchScalarField& field,
+        const phaseModel& phase,
+        const phaseModel& otherPhase
+    )
+    :
+        field(field),
+        phase(phase),
+        otherPhase(otherPhase),
+        alphaw(phase.boundaryField()[patchi()]),
+        otherAlphaw(otherPhase.boundaryField()[patchi()]),
+        interface(phase, otherPhase),
+        ttm
+        (
+            field.db().lookupType<fluidThermophysicalTransportModel>
+            (
+                phase.name()
+            )
+        ),
+        alphatConv
+        (
+            alphatJayatillekeWallFunctionFvPatchScalarField::alphat
+            (
+                ttm,
+                field.Prt_,
+                patchi()
+            )
+        )
+    {}
+
+
+    // Member Fuctions
+
+        //- Patch
+        inline const fvPatch& patch() const
+        {
+            return field.patch();
+        }
+
+        //- Patch index
+        inline label patchi() const
+        {
+            return patch().index();
+        }
+};
+
+
+struct alphatWallBoilingWallFunctionFvPatchScalarField::boilingLiquidProperties
+:
+    public alphatWallBoilingWallFunctionFvPatchScalarField::properties
+{
+    // Data
+
+        //- Name of the volatile specie
+        const word volatileSpecie;
+
+        //- Patch area by neighbouring cell volume ratio
+        const scalarField AbyV;
+
+        //- Liquid density
+        const tmp<scalarField> trhoLiquidw;
+        const scalarField& rhoLiquidw;
+
+        //- Vapour density
+        const tmp<scalarField> trhoVapourw;
+        const scalarField& rhoVapourw;
+
+        //- Liquid heat capacity
+        const scalarField& Cpw;
+
+        //- Liquid laminar kinematic viscosity
+        const tmp<scalarField> tnuw;
+        const scalarField& nuw;
+
+        //- Liquid laminar thermal diffusivity
+        const scalarField kappaByCp;
+
+        //- Liquid viscosity wall function
+        const nutWallFunctionFvPatchScalarField& nutw;
+
+        //- Dimensionless wall distance
+        const scalarField yPlus;
+
+        //- Smoothing function
+        const scalarField P;
+
+        //- Cell temperature
+        const scalarField Tc;
+
+        //- Saturation temperature
+        const scalarField Tsat;
+
+        //- Latent heat
+        const scalarField L;
+
+
+    //- Constructor
+    boilingLiquidProperties
+    (
+        const alphatWallBoilingWallFunctionFvPatchScalarField& field,
+        const phaseModel& liquid,
+        const phaseModel& vapour
+    )
+    :
+        properties(field, liquid, vapour),
+        volatileSpecie
+        (
+            liquid.fluid().lookupOrDefault<word>("volatile", "none")
+        ),
+        AbyV
+        (
+            patch().magSf()
+           /scalarField
+            (
+                patch().boundaryMesh().mesh().V(),
+                patch().faceCells()
+            )
+        ),
+        trhoLiquidw(liquid.thermo().rho(patchi())),
+        rhoLiquidw(trhoLiquidw()),
+        trhoVapourw(vapour.thermo().rho(patchi())),
+        rhoVapourw(trhoVapourw()),
+        Cpw(liquid.thermo().Cp().boundaryField()[patchi()]),
+        tnuw(liquid.thermo().nu(patchi())),
+        nuw(tnuw()),
+        kappaByCp
+        (
+            liquid.thermo().kappa().boundaryField()[patchi()]
+           /liquid.thermo().Cp().boundaryField()[patchi()]
+        ),
+        nutw
+        (
+            nutWallFunctionFvPatchScalarField::nutw
+            (
+                ttm.momentumTransport(),
+                patchi()
+            )
+        ),
+        yPlus
+        (
+            pow025(nutw.Cmu())
+           *sqrt(ttm.momentumTransport().k()().boundaryField()[patchi()])
+           *ttm.momentumTransport().y()[patchi()]
+           /nuw
+        ),
+        P
+        (
+            alphatJayatillekeWallFunctionFvPatchScalarField::P
+            (
+                rhoLiquidw*nuw/kappaByCp/field.Prt_
+            )
+        ),
+        Tc
+        (
+            liquid.thermo().T().boundaryField()[patchi()].patchInternalField()
+        ),
+        Tsat
+        (
+            liquid.fluid().lookupInterfacialModel
+            <
+                interfaceSaturationTemperatureModel
+            >
+            (interface)
+           .Tsat(liquid.thermo().p())()
+           .boundaryField()[patchi()]
+        ),
+        L
+        (
+            volatileSpecie != "none"
+         ? -refCast<const heatTransferPhaseSystem>(liquid.fluid())
+           .Li
+            (
+                interface,
+                volatileSpecie,
+                scalarField(patch().size(), +1),
+                Tsat,
+                patch().faceCells(),
+                heatTransferPhaseSystem::latentHeatScheme::upwind
+            )
+         : -refCast<const heatTransferPhaseSystem>(liquid.fluid())
+           .L
+            (
+                interface,
+                scalarField(patch().size(), +1),
+                Tsat,
+                patch().faceCells(),
+                heatTransferPhaseSystem::latentHeatScheme::upwind
+            )
+        )
+    {}
+};
+
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+tmp<scalarField>
+alphatWallBoilingWallFunctionFvPatchScalarField::calcBoiling
+(
+    const boilingLiquidProperties& props,
+    const scalarField& Tw,
+    scalarField& dDep,
+    scalarField& fDep,
+    scalarField& N,
+    scalarField& qq,
+    scalarField& qe,
+    scalarField& dmdtf
+) const
+{
+    scalarField Tl;
+    if (!useLiquidTemperatureWallFunction_)
+    {
+        Tl = props.Tc;
+    }
+    else
+    {
+        // Liquid temperature at y+=250 is estimated from the
+        // logarithmic thermal wall function of Koncar, Krepper
+        // & Egorov (2005)
+        const scalarField TyPlus250
+        (
+            Prt_
+           *(
+                log(props.nutw.E()*250)
+               /props.nutw.kappa()
+              + props.P
+           )
+        );
+
+        const scalarField TyPlus
+        (
+            Prt_
+           *(
+                log(props.nutw.E()*max(props.yPlus, scalar(11)))
+               /props.nutw.kappa()
+              + props.P
+            )
+        );
+
+        Tl = Tw - (TyPlus250/TyPlus)*(Tw - props.Tc);
+    }
+
+    // Bubble departure diameter
+    dDep =
+        departureDiameterModel_->dDeparture
+        (
+            props.phase,
+            props.otherPhase,
+            patch().index(),
+            Tl,
+            props.Tsat,
+            props.L
+        );
+
+    // Bubble departure frequency
+    fDep =
+        departureFrequencyModel_->fDeparture
+        (
+            props.phase,
+            props.otherPhase,
+            patch().index(),
+            Tl,
+            props.Tsat,
+            props.L,
+            dDep
+        );
+
+    // Nucleation site density
+    N =
+        nucleationSiteModel_->N
+        (
+            props.phase,
+            props.otherPhase,
+            patch().index(),
+            Tl,
+            props.Tsat,
+            props.L,
+            dDep,
+            fDep
+        );
+
+    // Del Valle & Kenning (1985)
+    const scalarField Ja
+    (
+        props.rhoLiquidw
+       *props.Cpw
+       *(props.Tsat - Tl)
+       /(props.rhoVapourw*props.L)
+    );
+
+    const scalarField Al
+    (
+        fLiquid_*4.8*exp(min(-Ja/80, log(vGreat)))
+    );
+
+    scalarField A2(min(pi*sqr(dDep)*N*Al/4, scalar(1)));
+    const scalarField A1(max(1 - A2, scalar(1e-4)));
+    scalarField A2E(min(pi*sqr(dDep)*N*Al/4, scalar(5)));
+
+    if (props.volatileSpecie != "none" && !props.phase.pure())
+    {
+        const scalarField& Yvolatile =
+            props.phase
+           .Y(props.volatileSpecie)
+           .boundaryField()[patch().index()];
+        A2E *= Yvolatile;
+        A2 *= Yvolatile;
+    }
+
+    // Volumetric mass source in the near wall cell due to the
+    // wall boiling
+    dmdtf = (1.0/6.0)*A2E*dDep*props.rhoVapourw*fDep*props.AbyV;
+
+    // Quenching heat transfer coefficient
+    const scalarField hQ
+    (
+        2*props.kappaByCp*props.Cpw*fDep
+       *sqrt((tau_/max(fDep, small))/(pi*props.kappaByCp/props.rhoLiquidw))
+    );
+
+    // Quenching heat flux
+    qq = A2*hQ*max(Tw - Tl, scalar(0));
+
+    // Evaporation heat flux
+    qe = dmdtf*props.L/props.AbyV;
+
+    // Return total sum of convective, quenching and evaporative heat fluxes
+    const scalarField gradTw
+    (
+        patch().deltaCoeffs()*max(Tw - props.Tc, small*props.Tc)
+    );
+    return A1*props.alphatConv*props.Cpw*gradTw + qq_ + qe_;
+}
+
+
+tmp<scalarField>
+alphatWallBoilingWallFunctionFvPatchScalarField::calcBoiling
+(
+    const boilingLiquidProperties& props,
+    const scalarField& Tw
+) const
+{
+    scalarField dDep(dDep_);
+    scalarField fDep(fDep_);
+    scalarField N(N_);
+    scalarField qq(qq_);
+    scalarField qe(qe_);
+    scalarField dmdtf(dmdtf_);
+
+    return calcBoiling(props, Tw, dDep, fDep, N, qq, qe, dmdtf);
+}
+
+
+tmp<scalarField>
+alphatWallBoilingWallFunctionFvPatchScalarField::evaluateBoiling
+(
+    const boilingLiquidProperties& props,
+    const scalarField& Tw
+)
+{
+    return calcBoiling(props, Tw, dDep_, fDep_, N_, qq_, qe_, dmdtf_);
+}
+
+
+const fvPatchScalarField&
+alphatWallBoilingWallFunctionFvPatchScalarField::getTemperaturePatchField
+(
+    const boilingLiquidProperties& props,
+    scalarField& isFixed,
+    scalarField& h,
+    scalarField& hTaPlusQa
+) const
+{
+    isFixed.setSize(patch().size());
+    h.setSize(patch().size());
+    hTaPlusQa.setSize(patch().size());
+
+    const fvPatchScalarField& Tw =
+        props.phase.thermo().T().boundaryField()[patch().index()];
+
+    if (isA<fixedValueFvPatchScalarField>(Tw))
+    {
+        isFixed = 1;
+        h = rootVGreat;
+        hTaPlusQa = rootVGreat*Tw;
+    }
+    else if (isA<zeroGradientFvPatchScalarField>(Tw))
+    {
+        isFixed = 0;
+        h = 0;
+        hTaPlusQa = 0;
+    }
+    else if (isA<fixedGradientFvPatchScalarField>(Tw))
+    {
+        const fixedGradientFvPatchScalarField& Twm =
+            refCast<const fixedGradientFvPatchScalarField>(Tw);
+
+        isFixed = 0;
+        h = 0;
+        hTaPlusQa = (*this)*props.Cpw*Twm.gradient();
+    }
+    else if (isA<mixedFvPatchScalarField>(Tw))
+    {
+        const mixedFvPatchScalarField& Twm =
+            refCast<const mixedFvPatchScalarField>(Tw);
+
+        isFixed = pos(Twm.valueFraction() - 1 + rootSmall);
+        h =
+            Twm.valueFraction()
+           /max(1 - Twm.valueFraction(), rootVSmall)
+           *(*this)*props.Cpw*patch().deltaCoeffs();
+        hTaPlusQa =
+            h*Twm.refValue()
+          + (*this)*props.Cpw*Twm.refGrad();
+    }
+    else
+    {
+        FatalErrorInFunction
+            << "Temperature boundary condition type not recognised"
+            << exit(FatalError);
+    }
+
+    return Tw;
+}
+
+
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
 alphatWallBoilingWallFunctionFvPatchScalarField::
@@ -77,7 +541,8 @@ alphatWallBoilingWallFunctionFvPatchScalarField
 
     phaseType_(liquidPhase),
     useLiquidTemperatureWallFunction_(true),
-    relax_(1),
+    tolerance_(rootSmall),
+
     Prt_(0.85),
     tau_(0.8),
 
@@ -112,7 +577,8 @@ alphatWallBoilingWallFunctionFvPatchScalarField
     (
         dict.lookupOrDefault<Switch>("useLiquidTemperatureWallFunction", true)
     ),
-    relax_(dict.lookupOrDefault<scalar>("relax", 1)),
+    tolerance_(dict.lookupOrDefault<scalar>("tolerance", rootSmall)),
+
     Prt_(dict.lookupOrDefault<scalar>("Prt", 0.85)),
     tau_(dict.lookupOrDefault<scalar>("bubbleWaitingTimeRatio", 0.8)),
 
@@ -179,9 +645,9 @@ alphatWallBoilingWallFunctionFvPatchScalarField
         {
             qq_ = scalarField("qQuenching", dict, p.size());
         }
-        if (dict.found("qEvaporation"))
+        if (dict.found("qEvaporative"))
         {
-            qe_ = scalarField("qEvaporation", dict, p.size());
+            qe_ = scalarField("qEvaporative", dict, p.size());
         }
         if (dict.found("dmdtf"))
         {
@@ -205,7 +671,8 @@ alphatWallBoilingWallFunctionFvPatchScalarField
 
     phaseType_(psf.phaseType_),
     useLiquidTemperatureWallFunction_(psf.useLiquidTemperatureWallFunction_),
-    relax_(psf.relax_),
+    tolerance_(psf.tolerance_),
+
     Prt_(psf.Prt_),
     tau_(psf.tau_),
 
@@ -236,7 +703,8 @@ alphatWallBoilingWallFunctionFvPatchScalarField
 
     phaseType_(psf.phaseType_),
     useLiquidTemperatureWallFunction_(psf.useLiquidTemperatureWallFunction_),
-    relax_(psf.relax_),
+    tolerance_(psf.tolerance_),
+
     Prt_(psf.Prt_),
     tau_(psf.tau_),
 
@@ -322,51 +790,29 @@ void alphatWallBoilingWallFunctionFvPatchScalarField::updateCoeffs()
         return;
     }
 
-    // Lookup the fluid model
+    // Lookup the fluid model and the phases
     const phaseSystem& fluid =
         db().lookupObject<phaseSystem>(phaseSystem::propertiesName);
 
-    const word volatileSpecie(fluid.lookupOrDefault<word>("volatile", "none"));
-
-    const label patchi = patch().index();
-
     switch (phaseType_)
     {
         case vaporPhase:
         {
             const phaseModel& vapor = fluid.phases()[internalField().group()];
+            const phaseModel& liquid = fluid.phases()[otherPhaseName_];
 
-            // Vapor thermophysical transport model
-            const fluidThermophysicalTransportModel& vaporTtm =
-                db().lookupType<fluidThermophysicalTransportModel>
-                (
-                    vapor.name()
-                );
-
-            // Vapor phase fraction at the wall
-            const scalarField& vaporw = vapor.boundaryField()[patchi];
+            // Construct boiling properties
+            const properties props(*this, vapor, liquid);
 
             // Partitioning. Note: Assumes that there is only only one liquid
             // phase and all 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()
-                )
-            );
+            fLiquid_ = partitioningModel_->fLiquid(props.otherAlphaw);
 
             operator==
             (
-                alphatConv*(vaporw/(1 - liquidw + small) )
-               *(1 - fLiquid_)/max(vaporw, scalar(1e-8))
+                (1 - fLiquid_)
+               /max(1 - props.otherAlphaw, rootSmall)
+               *props.alphatConv
             );
 
             break;
@@ -377,316 +823,100 @@ void alphatWallBoilingWallFunctionFvPatchScalarField::updateCoeffs()
             const phaseModel& liquid = fluid.phases()[internalField().group()];
             const phaseModel& vapor = fluid.phases()[otherPhaseName_];
 
-            const phaseInterface interface(vapor, liquid);
-
-            // Liquid thermophysical and momentum transport models
-            const fluidThermophysicalTransportModel& liquidTtm =
-                db().lookupType<fluidThermophysicalTransportModel>
-                (
-                    liquid.name()
-                );
-            const compressibleMomentumTransportModel& liquidMtm =
-                liquidTtm.momentumTransport();
-
-            // Convective thermal diffusivity
-            const scalarField alphatConv
-            (
-                alphatJayatillekeWallFunctionFvPatchScalarField::alphat
-                (
-                    liquidTtm,
-                    Prt_,
-                    patch().index()
-                )
-            );
-
-            // Quit if no saturation temperature model exists for this
-            // interface. Saturation modelling is considered to be the
-            // fundamental enabling model for thermal mass transfers.
+            // Boiling is enabled by the presence of saturation temperature
+            // modelling. This is consistent with interfacial thermal phase
+            // changes.
             if
             (
                !fluid.foundInterfacialModel
                 <
                     interfaceSaturationTemperatureModel
-                >
-                (interface)
+                >(phaseInterface(liquid, vapor))
             )
             {
-                Info<< "Saturation model for interface " << interface.name()
+                // Construct non-boiling properties
+                const properties props(*this, liquid, vapor);
+
+                Info<< "Saturation model for interface "
+                    << props.interface.name()
                     << " not found. Wall boiling disabled." << endl;
 
-                operator==(alphatConv);
-
-                break;
+                operator==(props.alphatConv);
             }
-
-            // Lookup and calculate turbulent and thermal properties
-            const nutWallFunctionFvPatchScalarField& nutw =
-                nutWallFunctionFvPatchScalarField::nutw(liquidMtm, patchi);
-
-            const scalar Cmu25(pow025(nutw.Cmu()));
-
-            const scalarField& y = liquidMtm.y()[patchi];
-
-            const tmp<scalarField> tnuw = liquidMtm.nu(patchi);
-            const scalarField& nuw = tnuw();
-
-            const tmp<scalarField> talphaw
-            (
-                liquid.thermo().kappa().boundaryField()[patchi]
-               /liquid.thermo().Cp().boundaryField()[patchi]
-            );
-            const scalarField& alphaw = talphaw();
-
-            const tmp<volScalarField> tk = liquidMtm.k();
-            const volScalarField& k = tk();
-            const fvPatchScalarField& kw = k.boundaryField()[patchi];
-
-            const fvPatchVectorField& Uw =
-                liquidMtm.U().boundaryField()[patchi];
-            const scalarField magUp(mag(Uw.patchInternalField() - Uw));
-            const scalarField magGradUw(mag(Uw.snGrad()));
-
-            const tmp<scalarField> trhoLiquidw = liquid.thermo().rho(patchi);
-            const scalarField rhoLiquidw = trhoLiquidw();
-
-            const tmp<scalarField> trhoVaporw = vapor.thermo().rho(patchi);
-            const scalarField rhoVaporw = trhoVaporw();
-
-            const fvPatchScalarField& hew =
-                liquid.thermo().he().boundaryField()[patchi];
-
-            const fvPatchScalarField& Tw =
-                liquid.thermo().T().boundaryField()[patchi];
-
-            const scalarField Tc(Tw.patchInternalField());
-
-            const scalarField uTau(Cmu25*sqrt(kw));
-
-            const scalarField yPlus(uTau*y/nuw);
-
-            const scalarField Pr(rhoLiquidw*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
-                )
-            );
-
-            const scalarField Cpw(liquid.thermo().Cp(Tw, patchi));
-
-            // Saturation temperature
-            const interfaceSaturationTemperatureModel& satModel =
-                fluid.lookupInterfacialModel
-                <
-                    interfaceSaturationTemperatureModel
-                >(interface);
-            const tmp<volScalarField> tTsat =
-                satModel.Tsat(liquid.thermo().p());
-            const volScalarField& Tsat = tTsat();
-            const fvPatchScalarField& Tsatw(Tsat.boundaryField()[patchi]);
-
-            // Latent heat
-            const scalarField L
-            (
-                volatileSpecie != "none"
-              ? -refCast<const heatTransferPhaseSystem>(fluid)
-                .Li
-                 (
-                     interface,
-                     volatileSpecie,
-                     dmdtf_,
-                     Tsat,
-                     patch().faceCells(),
-                     heatTransferPhaseSystem::latentHeatScheme::upwind
-                 )
-              : -refCast<const heatTransferPhaseSystem>(fluid)
-                .L
-                (
-                     interface,
-                     dmdtf_,
-                     Tsat,
-                     patch().faceCells(),
-                     heatTransferPhaseSystem::latentHeatScheme::upwind
-                 )
-            );
-
-            // Liquid phase fraction at the wall
-            const scalarField liquidw(liquid.boundaryField()[patchi]);
-
-            // Partitioning
-            fLiquid_ = partitioningModel_->fLiquid(liquidw);
-
-            // Iterative solution for the wall temperature
-            label maxIter(10);
-            for (label i=0; i<maxIter; i++)
+            else
             {
-                scalarField Tl(Tc);
+                // Construct boiling properties
+                const boilingLiquidProperties props(*this, liquid, vapor);
 
-                if (useLiquidTemperatureWallFunction_)
+                // Partitioning. Note: Assumes that there is only only one
+                // liquid phase and all other phases are vapor.
+                fLiquid_ = partitioningModel_->fLiquid(props.alphaw);
+
+                // Get the temperature boundary condition and extract its
+                // physical parameters
+                scalarField TwpfIsFixed, TwpfH, TwpfHTaPlusQa;
+                const fvPatchScalarField& Twpf =
+                    getTemperaturePatchField
+                    (
+                        props,
+                        TwpfIsFixed,
+                        TwpfH,
+                        TwpfHTaPlusQa
+                    );
+
+                // Define the residual. This should be monotonic in Tw.
+                auto R = [&](const scalarField& Tw)
                 {
-                    // Liquid temperature at y+=250 is estimated from the
-                    // logarithmic thermal wall function of Koncar, Krepper
-                    // & Egorov (2005)
-                    const scalarField TyPlus250
-                    (
-                        Prt_*(log(nutw.E()*250)/nutw.kappa() + P)
-                    );
+                    return calcBoiling(props, Tw) - TwpfHTaPlusQa + TwpfH*Tw;
+                };
 
-                    const scalarField TyPlus
-                    (
-                        Prt_
-                       *(
-                            log(nutw.E()*max(yPlus, scalar(11)))
-                           /nutw.kappa()
-                          + P
-                        )
-                    );
-
-                    Tl = Tw - (TyPlus250/TyPlus)*(Tw - Tc);
-                }
-
-                // Bubble departure diameter
-                dDep_ =
-                    departureDiameterModel_->dDeparture
-                    (
-                        liquid,
-                        vapor,
-                        patchi,
-                        Tl,
-                        Tsatw,
-                        L
-                    );
-
-                // Bubble departure frequency
-                fDep_ =
-                    departureFrequencyModel_->fDeparture
-                    (
-                        liquid,
-                        vapor,
-                        patchi,
-                        Tl,
-                        Tsatw,
-                        L,
-                        dDep_
-                    );
-
-                // Nucleation site density
-                N_ =
-                    nucleationSiteModel_->N
-                    (
-                        liquid,
-                        vapor,
-                        patchi,
-                        Tl,
-                        Tsatw,
-                        L,
-                        dDep_,
-                        fDep_
-                    );
-
-                // Del Valle & Kenning (1985)
-                const scalarField Ja
+                // Solve using interval bisection. Boiling cannot occur below
+                // the saturation temperature, so that is taken to be the lower
+                // bound. The upper bound is harder to define. For now, take
+                // twice the current superheat as the maximum. The solution is
+                // likely to be below this value. If it is not, then the
+                // iteration will converge to this upper limit, creating a new
+                // superheat of twice the current value. Subsequent time-steps
+                // will then double this superheat again and again until it
+                // does eventually bound the solution.
+                const scalarField isBoiling(neg(R(props.Tsat)));
+                scalarField Tw0(props.Tsat);
+                scalarField Tw1
                 (
-                    rhoLiquidw*Cpw*(Tsatw - Tl)/(rhoVaporw*L)
-                );
-
-                const scalarField Al
-                (
-                    fLiquid_*4.8*exp(min(-Ja/80, log(vGreat)))
-                );
-
-                scalarField A2(min(pi*sqr(dDep_)*N_*Al/4, scalar(1)));
-                const scalarField A1(max(1 - A2, scalar(1e-4)));
-                scalarField A2E(min(pi*sqr(dDep_)*N_*Al/4, scalar(5)));
-
-                if (volatileSpecie != "none" && !liquid.pure())
-                {
-                    const volScalarField& Yvolatile =
-                        liquid.Y(volatileSpecie);
-                    A2E *= Yvolatile.boundaryField()[patchi];
-                    A2 *= Yvolatile.boundaryField()[patchi];
-                }
-
-                // Patch area by neighbouring cell volume ratio
-                const scalarField AbyV
-                (
-                    patch().magSf()
-                   /scalarField
+                    max
                     (
-                        patch().boundaryMesh().mesh().V(),
-                        patch().faceCells()
+                        Twpf + (Twpf - props.Tsat),
+                        props.Tsat*(1 + sqrt(tolerance_))
                     )
                 );
-
-                // Volumetric mass source in the near wall cell due to the
-                // wall boiling
-                dmdtf_ =
-                    (1 - relax_)*dmdtf_
-                  + relax_*(1.0/6.0)*A2E*dDep_*rhoVaporw*fDep_*AbyV;
-
-                // Quenching heat transfer coefficient
-                const scalarField hQ
+                scalar e =
+                    gMax((1 - TwpfIsFixed)*isBoiling*(Tw1 - Tw0)/(Tw0 + Tw1));
+                for (; e > tolerance_; e /= 2)
+                {
+                    const scalarField TwM((Tw0 + Tw1)/2);
+                    const scalarField rM(R(TwM));
+                    Tw0 = pos(rM)*Tw0 + neg0(rM)*TwM;
+                    Tw1 = pos(rM)*TwM + neg0(rM)*Tw1;
+                }
+                const scalarField Tw
                 (
-                    2*alphaw*Cpw*fDep_
-                   *sqrt((tau_/max(fDep_, small))/(pi*alphaw/rhoLiquidw))
+                    TwpfIsFixed*Twpf + (1 - TwpfIsFixed)*(Tw0 + Tw1)/2
                 );
 
-                // Quenching heat flux
-                qq_ =
-                    (1 - relax_)*qq_
-                  + relax_*(A2*hQ*max(Tw - Tl, scalar(0)));
+                // Use solution to re-evaluate the boiling and set the thermal
+                // diffusivity to recover the calculated heat flux
+                const scalarField gradTw
+                (
+                    patch().deltaCoeffs()*max(Tw - props.Tc, small*props.Tc)
+                );
 
-                // Evaporation heat flux
-                qe_ = dmdtf_*L/AbyV;
+                const scalarField q(evaluateBoiling(props, Tw));
 
-                // Set an effective thermal diffusivity that corresponds to the
-                // calculated convective, quenching and evaporative heat fluxes
                 operator==
                 (
-                    (
-                        A1*alphatConv
-                      + (qq_ + qe_)/max(hew.snGrad(), scalar(1e-16))
-                    )
-                   /max(liquidw, scalar(1e-8))
+                    isBoiling*q/props.Cpw/gradTw/max(props.alphaw, rootSmall)
+                  + (1 - isBoiling)*props.alphatConv
                 );
-
-                // Evaluate the temperature condition and estimate the
-                // remaining residual error
-                const scalarField TsupPrev(max((Tw - Tsatw), scalar(0)));
-                const_cast<fvPatchScalarField&>(Tw).evaluate();
-                const scalarField TsupNew(max((Tw - Tsatw), scalar(0)));
-                const scalar maxErr(gMax(mag(TsupPrev - TsupNew)));
-
-                if (maxErr < 1e-1)
-                {
-                    if (i > 0)
-                    {
-                        Info<< "Wall boiling wall function iterations: "
-                            << i + 1 << endl;
-                    }
-
-                    break;
-                }
-
-                if (i == maxIter - 1)
-                {
-                    Info<< "Maximum number of wall boiling wall function "
-                        << "iterations (" << maxIter << ") reached." << endl
-                        << "Maximum change in wall temperature on last "
-                        << "iteration: " << maxErr << endl;
-                }
             }
 
             break;
@@ -710,7 +940,9 @@ void alphatWallBoilingWallFunctionFvPatchScalarField::write(Ostream& os) const
         "useLiquidTemperatureWallFunction",
         useLiquidTemperatureWallFunction_
     );
-    writeEntry(os, "relax", relax_);
+    writeEntry(os, "tolerance", tolerance_);
+
+    // Parameters
     writeEntry(os, "Prt", Prt_);
     writeEntry(os, "bubbleWaitingTimeRatio", tau_);
 
diff --git a/applications/solvers/modules/multiphaseEuler/multiphaseThermophysicalTransportModels/derivedFvPatchFields/alphatWallBoilingWallFunction/alphatWallBoilingWallFunctionFvPatchScalarField.H b/applications/solvers/modules/multiphaseEuler/multiphaseThermophysicalTransportModels/derivedFvPatchFields/alphatWallBoilingWallFunction/alphatWallBoilingWallFunctionFvPatchScalarField.H
index 73d625a62e..ba31b4a841 100644
--- a/applications/solvers/modules/multiphaseEuler/multiphaseThermophysicalTransportModels/derivedFvPatchFields/alphatWallBoilingWallFunction/alphatWallBoilingWallFunctionFvPatchScalarField.H
+++ b/applications/solvers/modules/multiphaseEuler/multiphaseThermophysicalTransportModels/derivedFvPatchFields/alphatWallBoilingWallFunction/alphatWallBoilingWallFunctionFvPatchScalarField.H
@@ -70,7 +70,7 @@ Usage
         phaseType      | 'vapor' or 'liquid'      | yes         |
         useLiquidTemperatureWallFunction | \\
         Use wall function to calculate liquid temperature? | no | yes
-        relax          | relaxation factor        | no          | 1
+        tolerance      | solution tolerance       | no          | rootSmall
         dmdt           | phase change mass flux   | no          | uniform 0
         Prt            | turbulent Prandtl number | no          | 0.85
     \endtable
@@ -95,9 +95,8 @@ Usage
     \verbatim
     hotWall
     {
-        type            compressible::alphatWallBoiling2WallFunction;
+        type            compressible::alphatWallBoilingWallFunction;
         phaseType       liquid;
-        relax           0.1;
         dmdt            uniform 0;
         Prt             0.85;
         partitioningModel
@@ -168,6 +167,15 @@ public:
 
 private:
 
+    // Private classes
+
+        //- Struct to hold cached properties for one of the phases
+        struct properties;
+
+        //- Struct to hold cached properties for the liquid when it is boiling
+        struct boilingLiquidProperties;
+
+
     // Private Data
 
         // Controls
@@ -178,8 +186,11 @@ private:
             //- Estimate liquid temperature using logarithmic wall function?
             Switch useLiquidTemperatureWallFunction_;
 
-            //- Relaxation factor
-            scalar relax_;
+            //- Solution tolerance
+            scalar tolerance_;
+
+
+        // Parameters
 
             //- Turbulent Prandtl number
             scalar Prt_;
@@ -231,6 +242,57 @@ private:
             scalarField dmdtf_;
 
 
+        // Private Member Functions
+
+            //- Calculate the boiling for the given wall temperature. Return
+            //  the total sum of all heat fluxes. Set the properties passed by
+            //  non-const reference. Used by the functions below.
+            tmp<scalarField> calcBoiling
+            (
+                const boilingLiquidProperties& props,
+                const scalarField& Tw,
+                scalarField& dDep,
+                scalarField& fDep,
+                scalarField& N,
+                scalarField& qq,
+                scalarField& qe,
+                scalarField& dmdtf
+            ) const;
+
+            //- Calculate the boiling for the given wall temperature. Return
+            //  the total sum of all heat fluxes. Use this to solve the balance
+            //  between the heat fluxes specified by the boiling models and the
+            //  temperature boundary condition without changing the stored
+            //  boiling state.
+            tmp<scalarField> calcBoiling
+            (
+                const boilingLiquidProperties& props,
+                const scalarField& Tw
+            ) const;
+
+            //- Calculate the boiling for the given wall temperature. Return
+            //  the total sum of all heat fluxes. Also set the stored boiling
+            //  state. Use this after solving with the final wall temperature
+            //  to set the boiling state.
+            tmp<scalarField> evaluateBoiling
+            (
+                const boilingLiquidProperties& props,
+                const scalarField& Tw
+            );
+
+            //- Get the temperature patch field and the parameters associated
+            //  with its boundary condition that are necessary for
+            //  approximately evaluating the boundary condition's heat flux at
+            //  a given wall temperature.
+            const fvPatchScalarField& getTemperaturePatchField
+            (
+                const boilingLiquidProperties& props,
+                scalarField& isFixed,
+                scalarField& h,
+                scalarField& hTaPlusQa
+            ) const;
+
+
 public:
 
     //- Runtime type information
diff --git a/applications/solvers/modules/multiphaseEuler/multiphaseThermophysicalTransportModels/derivedFvPatchFields/fixedMultiphaseHeatFlux/fixedMultiphaseHeatFluxFvPatchScalarField.C b/applications/solvers/modules/multiphaseEuler/multiphaseThermophysicalTransportModels/derivedFvPatchFields/fixedMultiphaseHeatFlux/fixedMultiphaseHeatFluxFvPatchScalarField.C
deleted file mode 100644
index 0734ddd859..0000000000
--- a/applications/solvers/modules/multiphaseEuler/multiphaseThermophysicalTransportModels/derivedFvPatchFields/fixedMultiphaseHeatFlux/fixedMultiphaseHeatFluxFvPatchScalarField.C
+++ /dev/null
@@ -1,214 +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 "fixedMultiphaseHeatFluxFvPatchScalarField.H"
-#include "fvPatchFieldMapper.H"
-#include "addToRunTimeSelectionTable.H"
-
-#include "phaseSystem.H"
-#include "compressibleMomentumTransportModels.H"
-#include "phaseCompressibleMomentumTransportModel.H"
-
-// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
-
-Foam::fixedMultiphaseHeatFluxFvPatchScalarField::
-fixedMultiphaseHeatFluxFvPatchScalarField
-(
-    const fvPatch& p,
-    const DimensionedField<scalar, volMesh>& iF
-)
-:
-    fixedValueFvPatchScalarField(p, iF),
-    q_(p.size(), 0.0),
-    relax_(1.0),
-    Tmin_(0.0)
-{}
-
-
-Foam::fixedMultiphaseHeatFluxFvPatchScalarField::
-fixedMultiphaseHeatFluxFvPatchScalarField
-(
-    const fvPatch& p,
-    const DimensionedField<scalar, volMesh>& iF,
-    const dictionary& dict
-)
-:
-    fixedValueFvPatchScalarField(p, iF, dict),
-    q_("q", dict, p.size()),
-    relax_(dict.lookupOrDefault<scalar>("relax", 1.0)),
-    Tmin_(dict.lookupOrDefault<scalar>("Tmin", 273))
-{}
-
-
-Foam::fixedMultiphaseHeatFluxFvPatchScalarField::
-fixedMultiphaseHeatFluxFvPatchScalarField
-(
-    const fixedMultiphaseHeatFluxFvPatchScalarField& psf,
-    const fvPatch& p,
-    const DimensionedField<scalar, volMesh>& iF,
-    const fvPatchFieldMapper& mapper
-)
-:
-    fixedValueFvPatchScalarField(psf, p, iF, mapper),
-    q_(mapper(psf.q_)),
-    relax_(psf.relax_),
-    Tmin_(psf.Tmin_)
-{}
-
-
-Foam::fixedMultiphaseHeatFluxFvPatchScalarField::
-fixedMultiphaseHeatFluxFvPatchScalarField
-(
-    const fixedMultiphaseHeatFluxFvPatchScalarField& psf,
-    const DimensionedField<scalar, volMesh>& iF
-)
-:
-    fixedValueFvPatchScalarField(psf, iF),
-    q_(psf.q_),
-    relax_(psf.relax_),
-    Tmin_(psf.Tmin_)
-{}
-
-
-
-// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
-
-void Foam::fixedMultiphaseHeatFluxFvPatchScalarField::updateCoeffs()
-{
-    if (updated())
-    {
-        return;
-    }
-
-    // Lookup the fluid model
-    const phaseSystem& fluid =
-            db().lookupObject<phaseSystem>(phaseSystem::propertiesName);
-
-    const scalarField& Tp = *this;
-
-    scalarField A(Tp.size(), scalar(0));
-    scalarField B(Tp.size(), scalar(0));
-    scalarField Q(Tp.size(), scalar(0));
-
-    forAll(fluid.phases(), phasei)
-    {
-        const phaseModel& phase = fluid.phases()[phasei];
-        const fluidThermo& thermo = phase.thermo();
-
-        const fvPatchScalarField& alpha =
-            phase.boundaryField()[patch().index()];
-
-        const fvPatchScalarField& T =
-            thermo.T().boundaryField()[patch().index()];
-
-        const scalarField kappaEff(phase.kappaEff(patch().index()));
-
-        if (debug)
-        {
-            scalarField q0(T.snGrad()*alpha*kappaEff);
-            Q += q0;
-
-            Info<< patch().name() << " " << phase.name()
-                << ": Heat flux " << gMin(q0) << " - " << gMax(q0) << endl;
-        }
-
-        A += T.patchInternalField()*alpha*kappaEff*patch().deltaCoeffs();
-        B += alpha*kappaEff*patch().deltaCoeffs();
-    }
-
-    if (debug)
-    {
-        Info<< patch().name() << " " << ": overall heat flux "
-            << gMin(Q) << " - " << gMax(Q) << " W/m2, power: "
-            << gSum(patch().magSf()*Q) << " W" << endl;
-    }
-
-    operator==((1 - relax_)*Tp + relax_*max(Tmin_,(q_ + A)/(B)));
-
-    fixedValueFvPatchScalarField::updateCoeffs();
-}
-
-
-void Foam::fixedMultiphaseHeatFluxFvPatchScalarField::autoMap
-(
-    const fvPatchFieldMapper& m
-)
-{
-    fixedValueFvPatchScalarField::autoMap(m);
-    m(q_, q_);
-}
-
-
-void Foam::fixedMultiphaseHeatFluxFvPatchScalarField::rmap
-(
-    const fvPatchScalarField& ptf,
-    const labelList& addr
-)
-{
-    fixedValueFvPatchScalarField::rmap(ptf, addr);
-
-    const fixedMultiphaseHeatFluxFvPatchScalarField& mptf =
-        refCast<const fixedMultiphaseHeatFluxFvPatchScalarField>(ptf);
-
-    q_.rmap(mptf.q_, addr);
-}
-
-
-void Foam::fixedMultiphaseHeatFluxFvPatchScalarField::reset
-(
-    const fvPatchScalarField& ptf
-)
-{
-    fixedValueFvPatchScalarField::reset(ptf);
-
-    const fixedMultiphaseHeatFluxFvPatchScalarField& mptf =
-        refCast<const fixedMultiphaseHeatFluxFvPatchScalarField>(ptf);
-
-    q_.reset(mptf.q_);
-}
-
-
-void Foam::fixedMultiphaseHeatFluxFvPatchScalarField::write(Ostream& os) const
-{
-    fvPatchField<scalar>::write(os);
-    writeEntry(os, "relax", relax_);
-    writeEntry(os, "q", q_);
-    writeEntry(os, "value", *this);
-}
-
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-namespace Foam
-{
-    makePatchTypeField
-    (
-        fvPatchScalarField,
-        fixedMultiphaseHeatFluxFvPatchScalarField
-    );
-}
-
-
-// ************************************************************************* //
diff --git a/applications/solvers/modules/multiphaseEuler/multiphaseThermophysicalTransportModels/derivedFvPatchFields/uniformFixedMultiphaseHeatFlux/uniformFixedMultiphaseHeatFluxFvPatchScalarField.C b/applications/solvers/modules/multiphaseEuler/multiphaseThermophysicalTransportModels/derivedFvPatchFields/uniformFixedMultiphaseHeatFlux/uniformFixedMultiphaseHeatFluxFvPatchScalarField.C
new file mode 100644
index 0000000000..e132410e01
--- /dev/null
+++ b/applications/solvers/modules/multiphaseEuler/multiphaseThermophysicalTransportModels/derivedFvPatchFields/uniformFixedMultiphaseHeatFlux/uniformFixedMultiphaseHeatFluxFvPatchScalarField.C
@@ -0,0 +1,198 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     | Website:  https://openfoam.org
+    \\  /    A nd           | Copyright (C) 2015-2023 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 "uniformFixedMultiphaseHeatFluxFvPatchScalarField.H"
+#include "fvPatchFieldMapper.H"
+#include "phaseSystem.H"
+#include "addToRunTimeSelectionTable.H"
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::uniformFixedMultiphaseHeatFluxFvPatchScalarField::
+uniformFixedMultiphaseHeatFluxFvPatchScalarField
+(
+    const fvPatch& p,
+    const DimensionedField<scalar, volMesh>& iF
+)
+:
+    mixedFvPatchScalarField(p, iF),
+    q_(nullptr),
+    relax_(1)
+{}
+
+
+Foam::uniformFixedMultiphaseHeatFluxFvPatchScalarField::
+uniformFixedMultiphaseHeatFluxFvPatchScalarField
+(
+    const fvPatch& p,
+    const DimensionedField<scalar, volMesh>& iF,
+    const dictionary& dict
+)
+:
+    mixedFvPatchScalarField(p, iF, dict, false),
+    q_(Function1<scalar>::New("q", dict)),
+    relax_(dict.lookupOrDefault<scalar>("relax", 1))
+{
+    valueFraction() = 1;
+    refValue() = patchInternalField();
+    refGrad() = Zero;
+
+    operator==(patchInternalField());
+}
+
+
+Foam::uniformFixedMultiphaseHeatFluxFvPatchScalarField::
+uniformFixedMultiphaseHeatFluxFvPatchScalarField
+(
+    const uniformFixedMultiphaseHeatFluxFvPatchScalarField& psf,
+    const fvPatch& p,
+    const DimensionedField<scalar, volMesh>& iF,
+    const fvPatchFieldMapper& mapper
+)
+:
+    mixedFvPatchScalarField(psf, p, iF, mapper),
+    q_(psf.q_, false),
+    relax_(psf.relax_)
+{}
+
+
+Foam::uniformFixedMultiphaseHeatFluxFvPatchScalarField::
+uniformFixedMultiphaseHeatFluxFvPatchScalarField
+(
+    const uniformFixedMultiphaseHeatFluxFvPatchScalarField& psf,
+    const DimensionedField<scalar, volMesh>& iF
+)
+:
+    mixedFvPatchScalarField(psf, iF),
+    q_(psf.q_, false),
+    relax_(psf.relax_)
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+void Foam::uniformFixedMultiphaseHeatFluxFvPatchScalarField::updateCoeffs()
+{
+    if (updated())
+    {
+        return;
+    }
+
+    const label patchi = patch().index();
+
+    const scalar q = q_->value(db().time().userTimeValue());
+
+    const phaseSystem& fluid =
+        db().lookupObject<phaseSystem>(phaseSystem::propertiesName);
+
+    const phaseModel& thisPhase = fluid.phases()[internalField().group()];
+
+    // Sums of alpha*kappaEff
+    scalarField sumAlphaKappaEff(patch().size(), rootVSmall);
+    scalarField sumNotThisAlphaKappaEff(patch().size(), rootVSmall);
+    scalarField sumNotThisAlphaKappaEffT(patch().size(), rootVSmall);
+
+    // Contributions from phases other than this one
+    forAll(fluid.phases(), phasei)
+    {
+        const phaseModel& phase = fluid.phases()[phasei];
+
+        if (&phase != &thisPhase)
+        {
+            const scalarField& alpha = phase.boundaryField()[patchi];
+            const scalarField kappaEff(phase.kappaEff(patchi));
+            const scalarField alphaKappaEff(alpha*kappaEff);
+            const fvPatchScalarField& T =
+                phase.thermo().T().boundaryField()[patchi];
+
+            sumAlphaKappaEff += alphaKappaEff;
+            sumNotThisAlphaKappaEff += alphaKappaEff;
+            sumNotThisAlphaKappaEffT += alphaKappaEff*T.patchInternalField();
+        }
+    }
+
+    // Contribution from this phase, and stabilisation
+    const scalarField& alpha = thisPhase.boundaryField()[patchi];
+    const scalarField kappaEff(thisPhase.kappaEff(patchi));
+    const scalarField alphaKappaEff(alpha*kappaEff);
+    const fvPatchScalarField& T =
+        thisPhase.thermo().T().boundaryField()[patchi];
+
+    sumAlphaKappaEff += alphaKappaEff;
+    sumNotThisAlphaKappaEff =
+        max
+        (
+            sumNotThisAlphaKappaEff,
+            rootSmall*kappaEff
+        );
+    sumNotThisAlphaKappaEffT =
+        max
+        (
+            sumNotThisAlphaKappaEffT,
+            rootSmall*kappaEff*T.patchInternalField()
+        );
+
+    // Mixed parameters
+    valueFraction() = sumNotThisAlphaKappaEff/sumAlphaKappaEff;
+    refValue() = sumNotThisAlphaKappaEffT/sumNotThisAlphaKappaEff;
+    refGrad() = q/max(alpha, rootSmall)/kappaEff;
+
+    // Modify mixed parameters for under-relaxation
+    if (relax_ != 1)
+    {
+        const scalarField f(valueFraction());
+        valueFraction() = 1 - relax_*(1 - f);
+        refValue() = (f*relax_*refValue() + (1 - relax_)*T)/valueFraction();
+        //refGrad() = refGrad(); // No change
+    }
+
+    mixedFvPatchScalarField::updateCoeffs();
+}
+
+
+void Foam::uniformFixedMultiphaseHeatFluxFvPatchScalarField::write
+(
+    Ostream& os
+) const
+{
+    mixedFvPatchScalarField::write(os);
+    writeEntry(os, q_());
+    writeEntry(os, "relax", relax_);
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+    makePatchTypeField
+    (
+        fvPatchScalarField,
+        uniformFixedMultiphaseHeatFluxFvPatchScalarField
+    );
+}
+
+
+// ************************************************************************* //
diff --git a/applications/solvers/modules/multiphaseEuler/multiphaseThermophysicalTransportModels/derivedFvPatchFields/fixedMultiphaseHeatFlux/fixedMultiphaseHeatFluxFvPatchScalarField.H b/applications/solvers/modules/multiphaseEuler/multiphaseThermophysicalTransportModels/derivedFvPatchFields/uniformFixedMultiphaseHeatFlux/uniformFixedMultiphaseHeatFluxFvPatchScalarField.H
similarity index 58%
rename from applications/solvers/modules/multiphaseEuler/multiphaseThermophysicalTransportModels/derivedFvPatchFields/fixedMultiphaseHeatFlux/fixedMultiphaseHeatFluxFvPatchScalarField.H
rename to applications/solvers/modules/multiphaseEuler/multiphaseThermophysicalTransportModels/derivedFvPatchFields/uniformFixedMultiphaseHeatFlux/uniformFixedMultiphaseHeatFluxFvPatchScalarField.H
index f2b4aeb48a..5a93875cd9 100644
--- a/applications/solvers/modules/multiphaseEuler/multiphaseThermophysicalTransportModels/derivedFvPatchFields/fixedMultiphaseHeatFlux/fixedMultiphaseHeatFluxFvPatchScalarField.H
+++ b/applications/solvers/modules/multiphaseEuler/multiphaseThermophysicalTransportModels/derivedFvPatchFields/uniformFixedMultiphaseHeatFlux/uniformFixedMultiphaseHeatFluxFvPatchScalarField.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
@@ -22,29 +22,44 @@ License
     along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
 
 Class
-    Foam::fixedMultiphaseHeatFluxFvPatchScalarField
+    Foam::uniformFixedMultiphaseHeatFluxFvPatchScalarField
 
 Description
-    Calculates a wall temperature that produces the specified overall wall heat
-    flux across all the phases in an Eulerian multi-phase simulation.
+    Uniform fixed heat flux boundary condition for Eulerian multi-phase cases.
+    Constructs a mixed constraint which portions the heat flux between the
+    phases in such a way as to keep the boundary temperature uniform across all
+    phases. The heat flux can be specified as a time-varying function, and an
+    under-relaxation factor can be supplied if this is necessary to maintain
+    stability.
 
-    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.
+Usage
+    \table
+        Property     | Description                | Required | Default value
+        q            | Heat flux [w/m^2]          | yes      |
+        relax        | Relaxation factor          | no       | 1
+    \endtable
 
-See also
-    Foam::fixedValueFvPatchField
+    Example of the boundary condition specification:
+    \verbatim
+    <patchName>
+    {
+        type            uniformFixedMultiphaseHeatFlux;
+        q               1000;
+        relax           0.3;
+        value           $internalField;
+    }
+    \endverbatim
 
 SourceFiles
-    fixedMultiphaseHeatFluxFvPatchScalarField.C
+    uniformFixedMultiphaseHeatFluxFvPatchScalarField.C
 
 \*---------------------------------------------------------------------------*/
 
-#ifndef fixedMultiphaseHeatFluxFvPatchScalarField_H
-#define fixedMultiphaseHeatFluxFvPatchScalarField_H
+#ifndef uniformFixedMultiphaseHeatFluxFvPatchScalarField_H
+#define uniformFixedMultiphaseHeatFluxFvPatchScalarField_H
 
-#include "fixedValueFvPatchFields.H"
+#include "mixedFvPatchFields.H"
+#include "Function1.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -52,42 +67,39 @@ namespace Foam
 {
 
 /*---------------------------------------------------------------------------*\
-          Class fixedMultiphaseHeatFluxFvPatchScalarField Declaration
+          Class uniformFixedMultiphaseHeatFluxFvPatchScalarField Declaration
 \*---------------------------------------------------------------------------*/
 
-class fixedMultiphaseHeatFluxFvPatchScalarField
+class uniformFixedMultiphaseHeatFluxFvPatchScalarField
 :
-    public fixedValueFvPatchScalarField
+    public mixedFvPatchScalarField
 {
     // Private Data
 
-        //- Heat power [W] or flux [W/m^2]
-        scalarField q_;
+        //- Heat flux [W/m^2]
+        autoPtr<Function1<scalar>> q_;
 
         //- Relaxation factor
         scalar relax_;
 
-        //- Minimum temperature limit [K]
-        scalar Tmin_;
-
 
 public:
 
     //- Runtime type information
-    TypeName("fixedMultiphaseHeatFlux");
+    TypeName("uniformFixedMultiphaseHeatFlux");
 
 
     // Constructors
 
         //- Construct from patch and internal field
-        fixedMultiphaseHeatFluxFvPatchScalarField
+        uniformFixedMultiphaseHeatFluxFvPatchScalarField
         (
             const fvPatch&,
             const DimensionedField<scalar, volMesh>&
         );
 
         //- Construct from patch, internal field and dictionary
-        fixedMultiphaseHeatFluxFvPatchScalarField
+        uniformFixedMultiphaseHeatFluxFvPatchScalarField
         (
             const fvPatch&,
             const DimensionedField<scalar, volMesh>&,
@@ -95,26 +107,26 @@ public:
         );
 
         //- Construct by mapping given
-        //  fixedMultiphaseHeatFluxFvPatchScalarField
+        //  uniformFixedMultiphaseHeatFluxFvPatchScalarField
         //  onto a new patch
-        fixedMultiphaseHeatFluxFvPatchScalarField
+        uniformFixedMultiphaseHeatFluxFvPatchScalarField
         (
-            const fixedMultiphaseHeatFluxFvPatchScalarField&,
+            const uniformFixedMultiphaseHeatFluxFvPatchScalarField&,
             const fvPatch&,
             const DimensionedField<scalar, volMesh>&,
             const fvPatchFieldMapper&
         );
 
         //- Disallow copy without setting internal field reference
-        fixedMultiphaseHeatFluxFvPatchScalarField
+        uniformFixedMultiphaseHeatFluxFvPatchScalarField
         (
-            const fixedMultiphaseHeatFluxFvPatchScalarField&
+            const uniformFixedMultiphaseHeatFluxFvPatchScalarField&
         ) = delete;
 
         //- Copy constructor setting internal field reference
-        fixedMultiphaseHeatFluxFvPatchScalarField
+        uniformFixedMultiphaseHeatFluxFvPatchScalarField
         (
-            const fixedMultiphaseHeatFluxFvPatchScalarField&,
+            const uniformFixedMultiphaseHeatFluxFvPatchScalarField&,
             const DimensionedField<scalar, volMesh>&
         );
 
@@ -126,28 +138,13 @@ public:
         {
             return tmp<fvPatchScalarField>
             (
-                new fixedMultiphaseHeatFluxFvPatchScalarField(*this, iF)
+                new uniformFixedMultiphaseHeatFluxFvPatchScalarField(*this, iF)
             );
         }
 
 
     // Member Functions
 
-        // Mapping functions
-
-            //- Map (and resize as needed) from self given a mapping object
-            //  Used to update fields following mesh topology change
-            virtual void autoMap(const fvPatchFieldMapper&);
-
-            //- Reverse map the given fvPatchField onto this fvPatchField
-            //  Used to reconstruct fields
-            virtual void rmap(const fvPatchScalarField&, const labelList&);
-
-            //- Reset the fvPatchField to the given fvPatchField
-            //  Used for mesh to mesh mapping
-            virtual void reset(const fvPatchScalarField&);
-
-
         // Evaluation functions
 
             //- Update the coefficients associated with the patch field
diff --git a/applications/solvers/modules/multiphaseEuler/phaseSystems/PhaseSystems/ThermalPhaseChangePhaseSystem/ThermalPhaseChangePhaseSystem.C b/applications/solvers/modules/multiphaseEuler/phaseSystems/PhaseSystems/ThermalPhaseChangePhaseSystem/ThermalPhaseChangePhaseSystem.C
index 484a5b272e..81b22e2b97 100644
--- a/applications/solvers/modules/multiphaseEuler/phaseSystems/PhaseSystems/ThermalPhaseChangePhaseSystem/ThermalPhaseChangePhaseSystem.C
+++ b/applications/solvers/modules/multiphaseEuler/phaseSystems/PhaseSystems/ThermalPhaseChangePhaseSystem/ThermalPhaseChangePhaseSystem.C
@@ -624,7 +624,9 @@ Foam::ThermalPhaseChangePhaseSystem<BasePhaseSystem>::correctInterfaceThermo()
 
         // Nucleation mass transfer update
         {
-            typedef compressible::alphatPhaseChangeWallFunctionBase alphatwType;
+            typedef
+                compressible::alphatPhaseChangeWallFunctionBase
+                alphatWallFunction;
 
             volScalarField& nDmdtf(*this->nDmdtfs_[interface]);
             nDmdtf = Zero;
@@ -649,10 +651,10 @@ Foam::ThermalPhaseChangePhaseSystem<BasePhaseSystem>::correctInterfaceThermo()
                     const fvPatchScalarField& alphatp =
                         alphat.boundaryField()[patchi];
 
-                    if (!isA<alphatwType>(alphatp)) continue;
+                    if (!isA<alphatWallFunction>(alphatp)) continue;
 
-                    const alphatwType& alphatw =
-                        refCast<const alphatwType>(alphatp);
+                    const alphatWallFunction& alphatw =
+                        refCast<const alphatWallFunction>(alphatp);
 
                     if (!alphatw.activeInterface(interface)) continue;
 
@@ -667,7 +669,7 @@ Foam::ThermalPhaseChangePhaseSystem<BasePhaseSystem>::correctInterfaceThermo()
                     nDmdtfp =
                         scalarField(nDmdtfp)
                       - (interfaceIter.index() == 0 ? +1 : -1)
-                       *alphatw.dmdtf(interface);
+                       *alphatw.dmdtf();
                 }
             }
 
diff --git a/applications/solvers/modules/multiphaseEuler/phaseSystems/populationBalanceModel/nucleationModels/wallBoiling/wallBoiling.C b/applications/solvers/modules/multiphaseEuler/phaseSystems/populationBalanceModel/nucleationModels/wallBoiling/wallBoiling.C
index 49a5c27d5d..4da9569618 100644
--- a/applications/solvers/modules/multiphaseEuler/phaseSystems/populationBalanceModel/nucleationModels/wallBoiling/wallBoiling.C
+++ b/applications/solvers/modules/multiphaseEuler/phaseSystems/populationBalanceModel/nucleationModels/wallBoiling/wallBoiling.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2018-2022 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2018-2023 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -78,10 +78,7 @@ void Foam::diameterModels::nucleationModels::wallBoiling::precompute()
 
     forAll(alphatBf, patchi)
     {
-        if
-        (
-            isA<alphatWallBoilingWallFunction>(alphatBf[patchi])
-        )
+        if (isA<alphatWallBoilingWallFunction>(alphatBf[patchi]))
         {
             const alphatWallBoilingWallFunction& alphatw =
                 refCast<const alphatWallBoilingWallFunction>(alphatBf[patchi]);
diff --git a/tutorials/modules/CHT/wallBoiling/0/fluid/alphat.gas b/tutorials/modules/CHT/wallBoiling/0/fluid/alphat.gas
index 248d4e028c..0bc8ec3dad 100644
--- a/tutorials/modules/CHT/wallBoiling/0/fluid/alphat.gas
+++ b/tutorials/modules/CHT/wallBoiling/0/fluid/alphat.gas
@@ -35,13 +35,10 @@ boundaryField
         otherPhase      liquid;
         phaseType       vapor;
         Prt             0.85;
-        Cmu             0.09;
-        kappa           0.41;
-        E               9.8;
         partitioningModel
         {
-            type        Lavieville;
-            alphaCrit   0.2;
+            type            Lavieville;
+            alphaCrit       0.2;
         }
         value           uniform 0;
     }
diff --git a/tutorials/modules/CHT/wallBoiling/0/fluid/alphat.liquid b/tutorials/modules/CHT/wallBoiling/0/fluid/alphat.liquid
index ed9dd4682f..ac08737bb0 100644
--- a/tutorials/modules/CHT/wallBoiling/0/fluid/alphat.liquid
+++ b/tutorials/modules/CHT/wallBoiling/0/fluid/alphat.liquid
@@ -32,31 +32,27 @@ boundaryField
     wall
     {
         type            compressible::alphatWallBoilingWallFunction;
-        Prt             0.85;
-        Cmu             0.09;
-        kappa           0.41;
-        E               9.8;
-        relax           1.0;
         otherPhase      gas;
         phaseType       liquid;
+        Prt             0.85;
         partitioningModel
         {
-            type        Lavieville;
-            alphaCrit   0.2;
+            type            Lavieville;
+            alphaCrit       0.2;
         }
         nucleationSiteModel
         {
-            type        LemmertChawla;
-            Cn          1;
-            NRef        30000000;
-            deltaTRef   10;
+            type            LemmertChawla;
+            Cn              1;
+            NRef            30000000;
+            deltaTRef       10;
         }
         departureDiamModel
         {
-            type        TolubinskiKostanchuk;
-            dRef        0.00024;
-            dMax        0.0014;
-            dMin        1e-06;
+            type            TolubinskiKostanchuk;
+            dRef            0.00024;
+            dMax            0.0014;
+            dMin            1e-06;
         }
         departureFreqModel
         {
diff --git a/tutorials/modules/CHT/wallBoiling/Allrun b/tutorials/modules/CHT/wallBoiling/Allrun
index 4ff27e983a..71a8c62fae 100755
--- a/tutorials/modules/CHT/wallBoiling/Allrun
+++ b/tutorials/modules/CHT/wallBoiling/Allrun
@@ -7,7 +7,7 @@ cd ${0%/*} || exit 1    # run from this directory
 runApplication blockMesh
 runApplication extrudeMesh
 runApplication splitMeshRegions -cellZones -overwrite
-runApplication foamDictionary constant/fluid/polyMesh/boundary -entry entry0/inlet/neighbourRegion -set "fluid"
+
 paraFoam -region fluid -touch
 paraFoam -region solid -touch
 
@@ -16,6 +16,7 @@ runApplication decomposePar -allRegions
 runParallel $(getApplication)
 
 runApplication reconstructPar -latestTime -allRegions
+
 runApplication foamPostProcess -latestTime -region fluid -func "
     graphCell
     (
@@ -24,9 +25,6 @@ runApplication foamPostProcess -latestTime -region fluid -func "
         end=(3.4901 0.0096 0),
         fields=(alpha.gas T.liquid T.gas)
     )"
-
-./validation/createGraphs
-
 runApplication -append foamPostProcess  -region fluid -latestTime -func "
     patchSurface
     (
@@ -37,6 +35,10 @@ runApplication -append foamPostProcess  -region fluid -latestTime -func "
         fields=(dDeparture.liquid fDeparture.liquid nucleationSiteDensity.liquid wetFraction.liquid qQuenching.liquid qEvaporative.liquid)
     )"
 
-./validation/createWallBoilingPropertiesGraphs
+if ! isTest "$@"
+then
+    ./validation/createGraphs
+    ./validation/createWallBoilingPropertiesGraphs
+fi
 
 #------------------------------------------------------------------------------
diff --git a/tutorials/modules/multiphaseEuler/wallBoilingIATE/0/T.gas b/tutorials/modules/multiphaseEuler/wallBoilingIATE/0/T.gas
index d148bc370c..ce14b165f4 100644
--- a/tutorials/modules/multiphaseEuler/wallBoilingIATE/0/T.gas
+++ b/tutorials/modules/multiphaseEuler/wallBoilingIATE/0/T.gas
@@ -34,9 +34,19 @@ boundaryField
     }
     wall
     {
-        type            mappedValue;
-        neighbourPatch  wall;
-        field           T.liquid;
+        type            uniformFixedMultiphaseHeatFlux;
+        q
+        {
+            type            scale;
+            value           73890;
+            scale
+            {
+                type            linearRamp;
+                start           0.5;
+                duration        0.01;
+            }
+        }
+        relax           0.3;
         value           $internalField;
     }
     front
diff --git a/tutorials/modules/multiphaseEuler/wallBoilingIATE/0/T.liquid b/tutorials/modules/multiphaseEuler/wallBoilingIATE/0/T.liquid
index 05ce63a7b9..8fd3dfaae5 100644
--- a/tutorials/modules/multiphaseEuler/wallBoilingIATE/0/T.liquid
+++ b/tutorials/modules/multiphaseEuler/wallBoilingIATE/0/T.liquid
@@ -34,10 +34,19 @@ boundaryField
     }
     wall
     {
-        type            fixedMultiphaseHeatFlux;
-        relax           0.6;
-        q               uniform 0;
-        phase           "liquid";
+        type            uniformFixedMultiphaseHeatFlux;
+        q
+        {
+            type            scale;
+            value           73890;
+            scale
+            {
+                type            linearRamp;
+                start           0.5;
+                duration        0.01;
+            }
+        }
+        relax           0.3;
         value           $internalField;
     }
     front
diff --git a/tutorials/modules/multiphaseEuler/wallBoilingIATE/0/alphat.gas b/tutorials/modules/multiphaseEuler/wallBoilingIATE/0/alphat.gas
index fd4ce93354..3fbb4b5bab 100644
--- a/tutorials/modules/multiphaseEuler/wallBoilingIATE/0/alphat.gas
+++ b/tutorials/modules/multiphaseEuler/wallBoilingIATE/0/alphat.gas
@@ -36,13 +36,10 @@ boundaryField
         otherPhase      liquid;
         phaseType       vapor;
         Prt             0.85;
-        Cmu             0.09;
-        kappa           0.41;
-        E               9.8;
         partitioningModel
         {
-            type        Lavieville;
-            alphaCrit   0.2;
+            type            Lavieville;
+            alphaCrit       0.2;
         }
         value           uniform 0;
     }
diff --git a/tutorials/modules/multiphaseEuler/wallBoilingIATE/0/alphat.liquid b/tutorials/modules/multiphaseEuler/wallBoilingIATE/0/alphat.liquid
index 3e00ea7cd1..ddbb10e0ee 100644
--- a/tutorials/modules/multiphaseEuler/wallBoilingIATE/0/alphat.liquid
+++ b/tutorials/modules/multiphaseEuler/wallBoilingIATE/0/alphat.liquid
@@ -33,31 +33,27 @@ boundaryField
     wall
     {
         type            compressible::alphatWallBoilingWallFunction;
-        Prt             0.85;
-        Cmu             0.09;
-        kappa           0.41;
-        E               9.8;
-        relax           0.6;
         otherPhase      gas;
         phaseType       liquid;
+        Prt             0.85;
         partitioningModel
         {
-            type        Lavieville;
-            alphaCrit   0.2;
+            type            Lavieville;
+            alphaCrit       0.2;
         }
         nucleationSiteModel
         {
-            type        LemmertChawla;
-            Cn          1;
-            NRef        30000000;
-            deltaTRef   10;
+            type            LemmertChawla;
+            Cn              1;
+            NRef            30000000;
+            deltaTRef       10;
         }
         departureDiamModel
         {
-            type        TolubinskiKostanchuk;
-            dRef        0.00024;
-            dMax        0.0014;
-            dMin        1e-06;
+            type            TolubinskiKostanchuk;
+            dRef            0.00024;
+            dMax            0.0014;
+            dMin            1e-06;
         }
         departureFreqModel
         {
diff --git a/tutorials/modules/multiphaseEuler/wallBoilingIATE/Allrun b/tutorials/modules/multiphaseEuler/wallBoilingIATE/Allrun
index bd96ed4575..02f2a6c4c2 100755
--- a/tutorials/modules/multiphaseEuler/wallBoilingIATE/Allrun
+++ b/tutorials/modules/multiphaseEuler/wallBoilingIATE/Allrun
@@ -10,16 +10,8 @@ runApplication decomposePar
 
 runParallel $(getApplication)
 
-if ! isTest "$@"
-then
-    runApplication -a foamDictionary system/controlDict -entry endTime -set 4
-    runApplication -a foamDictionary system/controlDict -entry startTime -set 0.5
-    runParallel -a foamDictionary 0.5/T.liquid -entry boundaryField/wall/q -set "uniform 73890"
-    runParallel -a foamDictionary 0.5/U.liquid -entry boundaryField/inlet/type -set "fixedValue"
-    runParallel -a $(getApplication)
-fi
-
 runApplication reconstructPar -latestTime
+
 runApplication foamPostProcess -latestTime -func "
     graphCell
     (
@@ -28,9 +20,6 @@ runApplication foamPostProcess -latestTime -func "
         end=(3.4901 0.0096 0),
         fields=(alpha.gas T.liquid T.gas d.gas)
     )"
-
-./validation/createGraphs
-
 runApplication -append foamPostProcess -latestTime -func "
     patchSurface
     (
@@ -41,6 +30,10 @@ runApplication -append foamPostProcess -latestTime -func "
         fields=(dDeparture.liquid fDeparture.liquid nucleationSiteDensity.liquid wetFraction.liquid qQuenching.liquid qEvaporative.liquid)
     )"
 
-./validation/createWallBoilingPropertiesGraphs
+if ! isTest "$@"
+then
+    ./validation/createGraphs
+    ./validation/createWallBoilingPropertiesGraphs
+fi
 
 #------------------------------------------------------------------------------
diff --git a/tutorials/modules/multiphaseEuler/wallBoilingIATE/system/controlDict.orig b/tutorials/modules/multiphaseEuler/wallBoilingIATE/system/controlDict
similarity index 99%
rename from tutorials/modules/multiphaseEuler/wallBoilingIATE/system/controlDict.orig
rename to tutorials/modules/multiphaseEuler/wallBoilingIATE/system/controlDict
index e3da7f558f..3f5180bdfd 100644
--- a/tutorials/modules/multiphaseEuler/wallBoilingIATE/system/controlDict.orig
+++ b/tutorials/modules/multiphaseEuler/wallBoilingIATE/system/controlDict
@@ -24,7 +24,7 @@ startTime       0;
 
 stopAt          endTime;
 
-endTime         0.5;
+endTime         4;
 
 deltaT          0.0001;
 
diff --git a/tutorials/modules/multiphaseEuler/wallBoilingPolydisperse/0/T.gas b/tutorials/modules/multiphaseEuler/wallBoilingPolydisperse/0/T.gas
index d148bc370c..ce14b165f4 100644
--- a/tutorials/modules/multiphaseEuler/wallBoilingPolydisperse/0/T.gas
+++ b/tutorials/modules/multiphaseEuler/wallBoilingPolydisperse/0/T.gas
@@ -34,9 +34,19 @@ boundaryField
     }
     wall
     {
-        type            mappedValue;
-        neighbourPatch  wall;
-        field           T.liquid;
+        type            uniformFixedMultiphaseHeatFlux;
+        q
+        {
+            type            scale;
+            value           73890;
+            scale
+            {
+                type            linearRamp;
+                start           0.5;
+                duration        0.01;
+            }
+        }
+        relax           0.3;
         value           $internalField;
     }
     front
diff --git a/tutorials/modules/multiphaseEuler/wallBoilingPolydisperse/0/T.liquid b/tutorials/modules/multiphaseEuler/wallBoilingPolydisperse/0/T.liquid
index 05ce63a7b9..8fd3dfaae5 100644
--- a/tutorials/modules/multiphaseEuler/wallBoilingPolydisperse/0/T.liquid
+++ b/tutorials/modules/multiphaseEuler/wallBoilingPolydisperse/0/T.liquid
@@ -34,10 +34,19 @@ boundaryField
     }
     wall
     {
-        type            fixedMultiphaseHeatFlux;
-        relax           0.6;
-        q               uniform 0;
-        phase           "liquid";
+        type            uniformFixedMultiphaseHeatFlux;
+        q
+        {
+            type            scale;
+            value           73890;
+            scale
+            {
+                type            linearRamp;
+                start           0.5;
+                duration        0.01;
+            }
+        }
+        relax           0.3;
         value           $internalField;
     }
     front
diff --git a/tutorials/modules/multiphaseEuler/wallBoilingPolydisperse/0/alphat.gas b/tutorials/modules/multiphaseEuler/wallBoilingPolydisperse/0/alphat.gas
index fd4ce93354..3fbb4b5bab 100644
--- a/tutorials/modules/multiphaseEuler/wallBoilingPolydisperse/0/alphat.gas
+++ b/tutorials/modules/multiphaseEuler/wallBoilingPolydisperse/0/alphat.gas
@@ -36,13 +36,10 @@ boundaryField
         otherPhase      liquid;
         phaseType       vapor;
         Prt             0.85;
-        Cmu             0.09;
-        kappa           0.41;
-        E               9.8;
         partitioningModel
         {
-            type        Lavieville;
-            alphaCrit   0.2;
+            type            Lavieville;
+            alphaCrit       0.2;
         }
         value           uniform 0;
     }
diff --git a/tutorials/modules/multiphaseEuler/wallBoilingPolydisperse/0/alphat.liquid b/tutorials/modules/multiphaseEuler/wallBoilingPolydisperse/0/alphat.liquid
index 3e00ea7cd1..ddbb10e0ee 100644
--- a/tutorials/modules/multiphaseEuler/wallBoilingPolydisperse/0/alphat.liquid
+++ b/tutorials/modules/multiphaseEuler/wallBoilingPolydisperse/0/alphat.liquid
@@ -33,31 +33,27 @@ boundaryField
     wall
     {
         type            compressible::alphatWallBoilingWallFunction;
-        Prt             0.85;
-        Cmu             0.09;
-        kappa           0.41;
-        E               9.8;
-        relax           0.6;
         otherPhase      gas;
         phaseType       liquid;
+        Prt             0.85;
         partitioningModel
         {
-            type        Lavieville;
-            alphaCrit   0.2;
+            type            Lavieville;
+            alphaCrit       0.2;
         }
         nucleationSiteModel
         {
-            type        LemmertChawla;
-            Cn          1;
-            NRef        30000000;
-            deltaTRef   10;
+            type            LemmertChawla;
+            Cn              1;
+            NRef            30000000;
+            deltaTRef       10;
         }
         departureDiamModel
         {
-            type        TolubinskiKostanchuk;
-            dRef        0.00024;
-            dMax        0.0014;
-            dMin        1e-06;
+            type            TolubinskiKostanchuk;
+            dRef            0.00024;
+            dMax            0.0014;
+            dMin            1e-06;
         }
         departureFreqModel
         {
diff --git a/tutorials/modules/multiphaseEuler/wallBoilingPolydisperse/Allrun b/tutorials/modules/multiphaseEuler/wallBoilingPolydisperse/Allrun
index 59ef3762d6..26ba676605 100755
--- a/tutorials/modules/multiphaseEuler/wallBoilingPolydisperse/Allrun
+++ b/tutorials/modules/multiphaseEuler/wallBoilingPolydisperse/Allrun
@@ -11,16 +11,8 @@ runApplication decomposePar
 
 runParallel $(getApplication)
 
-if ! isTest "$@"
-then
-    runApplication -a foamDictionary system/controlDict -entry endTime -set 4
-    runApplication -a foamDictionary system/controlDict -entry startTime -set 0.5
-    runParallel -a foamDictionary 0.5/T.liquid -entry boundaryField/wall/q -set "uniform 73890"
-    runParallel -a foamDictionary 0.5/U.liquid -entry boundaryField/inlet/type -set "fixedValue"
-    runParallel -a $(getApplication)
-fi
-
 runApplication reconstructPar -latestTime
+
 runApplication foamPostProcess -latestTime -func "
     graphCell
     (
@@ -29,9 +21,6 @@ runApplication foamPostProcess -latestTime -func "
         end=(3.4901 0.0096 0),
         fields=(alpha.gas T.liquid T.gas d.gas)
     )"
-
-./validation/createGraphs
-
 runApplication -append foamPostProcess -latestTime -func "
     patchSurface
     (
@@ -42,6 +31,10 @@ runApplication -append foamPostProcess -latestTime -func "
         fields=(dDeparture.liquid fDeparture.liquid nucleationSiteDensity.liquid wetFraction.liquid qQuenching.liquid qEvaporative.liquid)
     )"
 
-./validation/createWallBoilingPropertiesGraphs
+if ! isTest "$@"
+then
+    ./validation/createGraphs
+    ./validation/createWallBoilingPropertiesGraphs
+fi
 
 #------------------------------------------------------------------------------
diff --git a/tutorials/modules/multiphaseEuler/wallBoilingPolydisperse/system/controlDict.orig b/tutorials/modules/multiphaseEuler/wallBoilingPolydisperse/system/controlDict
similarity index 99%
rename from tutorials/modules/multiphaseEuler/wallBoilingPolydisperse/system/controlDict.orig
rename to tutorials/modules/multiphaseEuler/wallBoilingPolydisperse/system/controlDict
index aa211b5811..af87f5da21 100644
--- a/tutorials/modules/multiphaseEuler/wallBoilingPolydisperse/system/controlDict.orig
+++ b/tutorials/modules/multiphaseEuler/wallBoilingPolydisperse/system/controlDict
@@ -24,7 +24,7 @@ startTime       0;
 
 stopAt          endTime;
 
-endTime         0.5;
+endTime         4;
 
 deltaT          0.0001;
 
diff --git a/tutorials/modules/multiphaseEuler/wallBoilingPolydisperseTwoGroups/0/T.gas b/tutorials/modules/multiphaseEuler/wallBoilingPolydisperseTwoGroups/0/T.gas
index d148bc370c..ce14b165f4 100644
--- a/tutorials/modules/multiphaseEuler/wallBoilingPolydisperseTwoGroups/0/T.gas
+++ b/tutorials/modules/multiphaseEuler/wallBoilingPolydisperseTwoGroups/0/T.gas
@@ -34,9 +34,19 @@ boundaryField
     }
     wall
     {
-        type            mappedValue;
-        neighbourPatch  wall;
-        field           T.liquid;
+        type            uniformFixedMultiphaseHeatFlux;
+        q
+        {
+            type            scale;
+            value           73890;
+            scale
+            {
+                type            linearRamp;
+                start           0.5;
+                duration        0.01;
+            }
+        }
+        relax           0.3;
         value           $internalField;
     }
     front
diff --git a/tutorials/modules/multiphaseEuler/wallBoilingPolydisperseTwoGroups/0/T.gas2 b/tutorials/modules/multiphaseEuler/wallBoilingPolydisperseTwoGroups/0/T.gas2
index ff9a4b7859..a3e4b002bb 100644
--- a/tutorials/modules/multiphaseEuler/wallBoilingPolydisperseTwoGroups/0/T.gas2
+++ b/tutorials/modules/multiphaseEuler/wallBoilingPolydisperseTwoGroups/0/T.gas2
@@ -10,7 +10,7 @@ FoamFile
     format      ascii;
     class       volScalarField;
     location    "5";
-    object      T.gas2;
+    object      T.gas;
 }
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -34,9 +34,19 @@ boundaryField
     }
     wall
     {
-        type            mappedValue;
-        neighbourPatch  wall;
-        field           T.liquid;
+        type            uniformFixedMultiphaseHeatFlux;
+        q
+        {
+            type            scale;
+            value           73890;
+            scale
+            {
+                type            linearRamp;
+                start           0.5;
+                duration        0.01;
+            }
+        }
+        relax           0.3;
         value           $internalField;
     }
     front
diff --git a/tutorials/modules/multiphaseEuler/wallBoilingPolydisperseTwoGroups/0/T.liquid b/tutorials/modules/multiphaseEuler/wallBoilingPolydisperseTwoGroups/0/T.liquid
index 05ce63a7b9..8fd3dfaae5 100644
--- a/tutorials/modules/multiphaseEuler/wallBoilingPolydisperseTwoGroups/0/T.liquid
+++ b/tutorials/modules/multiphaseEuler/wallBoilingPolydisperseTwoGroups/0/T.liquid
@@ -34,10 +34,19 @@ boundaryField
     }
     wall
     {
-        type            fixedMultiphaseHeatFlux;
-        relax           0.6;
-        q               uniform 0;
-        phase           "liquid";
+        type            uniformFixedMultiphaseHeatFlux;
+        q
+        {
+            type            scale;
+            value           73890;
+            scale
+            {
+                type            linearRamp;
+                start           0.5;
+                duration        0.01;
+            }
+        }
+        relax           0.3;
         value           $internalField;
     }
     front
diff --git a/tutorials/modules/multiphaseEuler/wallBoilingPolydisperseTwoGroups/0/alphat.gas b/tutorials/modules/multiphaseEuler/wallBoilingPolydisperseTwoGroups/0/alphat.gas
index fd4ce93354..3fbb4b5bab 100644
--- a/tutorials/modules/multiphaseEuler/wallBoilingPolydisperseTwoGroups/0/alphat.gas
+++ b/tutorials/modules/multiphaseEuler/wallBoilingPolydisperseTwoGroups/0/alphat.gas
@@ -36,13 +36,10 @@ boundaryField
         otherPhase      liquid;
         phaseType       vapor;
         Prt             0.85;
-        Cmu             0.09;
-        kappa           0.41;
-        E               9.8;
         partitioningModel
         {
-            type        Lavieville;
-            alphaCrit   0.2;
+            type            Lavieville;
+            alphaCrit       0.2;
         }
         value           uniform 0;
     }
diff --git a/tutorials/modules/multiphaseEuler/wallBoilingPolydisperseTwoGroups/0/alphat.gas2 b/tutorials/modules/multiphaseEuler/wallBoilingPolydisperseTwoGroups/0/alphat.gas2
index fd4ce93354..3fbb4b5bab 100644
--- a/tutorials/modules/multiphaseEuler/wallBoilingPolydisperseTwoGroups/0/alphat.gas2
+++ b/tutorials/modules/multiphaseEuler/wallBoilingPolydisperseTwoGroups/0/alphat.gas2
@@ -36,13 +36,10 @@ boundaryField
         otherPhase      liquid;
         phaseType       vapor;
         Prt             0.85;
-        Cmu             0.09;
-        kappa           0.41;
-        E               9.8;
         partitioningModel
         {
-            type        Lavieville;
-            alphaCrit   0.2;
+            type            Lavieville;
+            alphaCrit       0.2;
         }
         value           uniform 0;
     }
diff --git a/tutorials/modules/multiphaseEuler/wallBoilingPolydisperseTwoGroups/0/alphat.liquid b/tutorials/modules/multiphaseEuler/wallBoilingPolydisperseTwoGroups/0/alphat.liquid
index 3e00ea7cd1..ddbb10e0ee 100644
--- a/tutorials/modules/multiphaseEuler/wallBoilingPolydisperseTwoGroups/0/alphat.liquid
+++ b/tutorials/modules/multiphaseEuler/wallBoilingPolydisperseTwoGroups/0/alphat.liquid
@@ -33,31 +33,27 @@ boundaryField
     wall
     {
         type            compressible::alphatWallBoilingWallFunction;
-        Prt             0.85;
-        Cmu             0.09;
-        kappa           0.41;
-        E               9.8;
-        relax           0.6;
         otherPhase      gas;
         phaseType       liquid;
+        Prt             0.85;
         partitioningModel
         {
-            type        Lavieville;
-            alphaCrit   0.2;
+            type            Lavieville;
+            alphaCrit       0.2;
         }
         nucleationSiteModel
         {
-            type        LemmertChawla;
-            Cn          1;
-            NRef        30000000;
-            deltaTRef   10;
+            type            LemmertChawla;
+            Cn              1;
+            NRef            30000000;
+            deltaTRef       10;
         }
         departureDiamModel
         {
-            type        TolubinskiKostanchuk;
-            dRef        0.00024;
-            dMax        0.0014;
-            dMin        1e-06;
+            type            TolubinskiKostanchuk;
+            dRef            0.00024;
+            dMax            0.0014;
+            dMin            1e-06;
         }
         departureFreqModel
         {
diff --git a/tutorials/modules/multiphaseEuler/wallBoilingPolydisperseTwoGroups/Allrun b/tutorials/modules/multiphaseEuler/wallBoilingPolydisperseTwoGroups/Allrun
index 9935a96442..b268255cb8 100755
--- a/tutorials/modules/multiphaseEuler/wallBoilingPolydisperseTwoGroups/Allrun
+++ b/tutorials/modules/multiphaseEuler/wallBoilingPolydisperseTwoGroups/Allrun
@@ -11,16 +11,8 @@ runApplication decomposePar
 
 runParallel $(getApplication)
 
-if ! isTest "$@"
-then
-    runApplication -a foamDictionary system/controlDict -entry endTime -set 4
-    runApplication -a foamDictionary system/controlDict -entry startTime -set 0.5
-    runParallel -a foamDictionary 0.5/T.liquid -entry boundaryField/wall/q -set "uniform 73890"
-    runParallel -a foamDictionary 0.5/U.liquid -entry boundaryField/inlet/type -set "fixedValue"
-    runParallel -a $(getApplication)
-fi
-
 runApplication reconstructPar -latestTime
+
 runApplication foamPostProcess -latestTime -func "
     graphCell
     (
@@ -29,9 +21,6 @@ runApplication foamPostProcess -latestTime -func "
         end=(3.4901 0.0096 0),
         fields=(alpha.gas alpha.gas2 alpha.liquid T.liquid T.gas d.bubbles)
     )"
-
-./validation/createGraphs
-
 runApplication -append foamPostProcess -latestTime -func "
     patchSurface
     (
@@ -42,6 +31,10 @@ runApplication -append foamPostProcess -latestTime -func "
         fields=(dDeparture.liquid fDeparture.liquid nucleationSiteDensity.liquid wetFraction.liquid qQuenching.liquid qEvaporative.liquid)
     )"
 
-./validation/createWallBoilingPropertiesGraphs
+if ! isTest "$@"
+then
+    ./validation/createGraphs
+    ./validation/createWallBoilingPropertiesGraphs
+fi
 
 #------------------------------------------------------------------------------
diff --git a/tutorials/modules/multiphaseEuler/wallBoilingPolydisperseTwoGroups/system/controlDict.orig b/tutorials/modules/multiphaseEuler/wallBoilingPolydisperseTwoGroups/system/controlDict
similarity index 99%
rename from tutorials/modules/multiphaseEuler/wallBoilingPolydisperseTwoGroups/system/controlDict.orig
rename to tutorials/modules/multiphaseEuler/wallBoilingPolydisperseTwoGroups/system/controlDict
index d593c66f8d..562fdcf3f6 100644
--- a/tutorials/modules/multiphaseEuler/wallBoilingPolydisperseTwoGroups/system/controlDict.orig
+++ b/tutorials/modules/multiphaseEuler/wallBoilingPolydisperseTwoGroups/system/controlDict
@@ -24,7 +24,7 @@ startTime       0;
 
 stopAt          endTime;
 
-endTime         0.5;
+endTime         4;
 
 deltaT          0.0001;