diff --git a/applications/solvers/multiphase/reactingEulerFoam/Allwclean b/applications/solvers/multiphase/reactingEulerFoam/Allwclean
index 2fe4eb69c0e1a0f2a1959bc866580af9971686aa..4ac3e56e70a6c5b6ce27eaafdf2595162746148f 100755
--- a/applications/solvers/multiphase/reactingEulerFoam/Allwclean
+++ b/applications/solvers/multiphase/reactingEulerFoam/Allwclean
@@ -3,6 +3,5 @@ cd "${0%/*}" || exit    # Run from this directory
 
 wclean reactingTwoPhaseEulerFoam
 wclean reactingMultiphaseEulerFoam
-wclean libso functionObjects
 
 #------------------------------------------------------------------------------
diff --git a/applications/solvers/multiphase/reactingEulerFoam/Allwmake b/applications/solvers/multiphase/reactingEulerFoam/Allwmake
index 22b8bfef5ce1dc9b924b11849fe4c2c8b331a18f..afb6570f3afccad85b3fbdd5080f7e0343c72ed1 100755
--- a/applications/solvers/multiphase/reactingEulerFoam/Allwmake
+++ b/applications/solvers/multiphase/reactingEulerFoam/Allwmake
@@ -16,6 +16,5 @@ esac
 
 reactingTwoPhaseEulerFoam/Allwmake $targetType $*
 reactingMultiphaseEulerFoam/Allwmake $targetType $*
-wmake $targetType functionObjects
 
 #------------------------------------------------------------------------------
diff --git a/src/Allwmake b/src/Allwmake
index 812b206ecbfb6a699863635518ba16684e76b6c3..c2baebc1cae14eefe1630fb340a44e8948f29bd1 100755
--- a/src/Allwmake
+++ b/src/Allwmake
@@ -100,8 +100,8 @@ wmake $targetType engine
 
 conversion/Allwmake $targetType $*
 
-phaseSystemModels/Allwmake $targetType $*
 functionObjects/Allwmake $targetType $*
+phaseSystemModels/Allwmake $targetType $*
 
 wmake $targetType lumpedPointMotion
 wmake $targetType sixDoFRigidBodyMotion
diff --git a/src/functionObjects/field/Make/options b/src/functionObjects/field/Make/options
index 03f219f2604201beea9b231148617c8ef13b4342..c24acfb63201d0ff95ac4e860474a2a9bb544fe1 100644
--- a/src/functionObjects/field/Make/options
+++ b/src/functionObjects/field/Make/options
@@ -38,6 +38,4 @@ LIB_LIBS = \
     -lcompressibleTurbulenceModels \
     -lchemistryModel \
     -lreactionThermophysicalModels \
-    -lpairPatchAgglomeration \
-    -ltwoPhaseReactingTurbulenceModels \
-    -lreactingPhaseSystem
+    -lpairPatchAgglomeration
diff --git a/src/functionObjects/field/heatTransferCoeff/heatTransferCoeff.C b/src/functionObjects/field/heatTransferCoeff/heatTransferCoeff.C
index 9f66f0a86c6eb17eeaab813c873afae19871f99c..d9b48b62d17dfce81e9ec7d3bab0c1373319ce7f 100644
--- a/src/functionObjects/field/heatTransferCoeff/heatTransferCoeff.C
+++ b/src/functionObjects/field/heatTransferCoeff/heatTransferCoeff.C
@@ -5,7 +5,7 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2017 OpenCFD Ltd.
+    Copyright (C) 2017-2020 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -47,7 +47,7 @@ bool Foam::functionObjects::heatTransferCoeff::calc()
 {
     volScalarField& htc = mesh_.lookupObjectRef<volScalarField>(resultName_);
 
-    htcModelPtr_->calc(htc);
+    htcModelPtr_->calc(htc, htcModelPtr_->q());
 
     return true;
 }
@@ -63,14 +63,13 @@ Foam::functionObjects::heatTransferCoeff::heatTransferCoeff
 )
 :
     fieldExpression(name, runTime, dict),
-    htcModelPtr_()
+    htcModelPtr_(nullptr)
 {
     read(dict);
 
     setResultName(typeName, name + ":htc:" + htcModelPtr_->type());
 
-    volScalarField* heatTransferCoeffPtr
-    (
+    volScalarField* heatTransferCoeffPtr =
         new volScalarField
         (
             IOobject
@@ -83,8 +82,7 @@ Foam::functionObjects::heatTransferCoeff::heatTransferCoeff
             ),
             mesh_,
             dimensionedScalar(dimPower/dimArea/dimTemperature, Zero)
-        )
-    );
+        );
 
     mesh_.objectRegistry::store(heatTransferCoeffPtr);
 }
diff --git a/src/functionObjects/field/heatTransferCoeff/heatTransferCoeff.H b/src/functionObjects/field/heatTransferCoeff/heatTransferCoeff.H
index 01c42773305bdee832a34ab6fce756fdadd31b1f..08cf4258fa73724c9387057f59e14c06e645ebb8 100644
--- a/src/functionObjects/field/heatTransferCoeff/heatTransferCoeff.H
+++ b/src/functionObjects/field/heatTransferCoeff/heatTransferCoeff.H
@@ -122,6 +122,7 @@ SourceFiles
 namespace Foam
 {
 
+// Forward Declarations
 class heatTransferCoeffModel;
 
 namespace functionObjects
@@ -135,16 +136,19 @@ class heatTransferCoeff
 :
     public fieldExpression
 {
-
-private:
-
-    // Private data
+    // Private Data
 
         //- Heat transfer coefficient model
         autoPtr<heatTransferCoeffModel> htcModelPtr_;
 
 
-    // Private Member Functions
+protected:
+
+    // Protected Member Functions
+
+        //- Calculate the heat transfer coefficient field
+        //  \return true on success
+        virtual bool calc();
 
         //- No copy construct
         heatTransferCoeff(const heatTransferCoeff&) = delete;
@@ -153,13 +157,6 @@ private:
         void operator=(const heatTransferCoeff&) = delete;
 
 
-protected:
-
-        //- Calculate the heat transfer coefficient field and return true
-        // if successful
-        virtual bool calc();
-
-
 public:
 
     //- Runtime type information
@@ -185,7 +182,7 @@ public:
     // Member Functions
 
         //- Read the heatTransferCoeff data
-        virtual bool read(const dictionary&);
+        virtual bool read(const dictionary& dict);
 };
 
 
diff --git a/src/functionObjects/field/heatTransferCoeff/heatTransferCoeffModels/ReynoldsAnalogy/ReynoldsAnalogy.C b/src/functionObjects/field/heatTransferCoeff/heatTransferCoeffModels/ReynoldsAnalogy/ReynoldsAnalogy.C
index 3cd1f32642828098d68c9133e8bd6c4fab312a78..a32e632c384d234bcfebcd8e5d0f571ac47ee90d 100644
--- a/src/functionObjects/field/heatTransferCoeff/heatTransferCoeffModels/ReynoldsAnalogy/ReynoldsAnalogy.C
+++ b/src/functionObjects/field/heatTransferCoeff/heatTransferCoeffModels/ReynoldsAnalogy/ReynoldsAnalogy.C
@@ -243,7 +243,11 @@ bool Foam::heatTransferCoeffModels::ReynoldsAnalogy::read
 }
 
 
-void Foam::heatTransferCoeffModels::ReynoldsAnalogy::htc(volScalarField& htc)
+void Foam::heatTransferCoeffModels::ReynoldsAnalogy::htc
+(
+    volScalarField& htc,
+    const FieldField<Field, scalar>& q
+)
 {
     const FieldField<Field, scalar> CfBf(Cf());
     const scalar magU = mag(URef_);
diff --git a/src/functionObjects/field/heatTransferCoeff/heatTransferCoeffModels/ReynoldsAnalogy/ReynoldsAnalogy.H b/src/functionObjects/field/heatTransferCoeff/heatTransferCoeffModels/ReynoldsAnalogy/ReynoldsAnalogy.H
index 3e9b463622dae95eb6e91bc148c4fca462677346..1f227dd77696f57972266452667ed80b2623666e 100644
--- a/src/functionObjects/field/heatTransferCoeff/heatTransferCoeffModels/ReynoldsAnalogy/ReynoldsAnalogy.H
+++ b/src/functionObjects/field/heatTransferCoeff/heatTransferCoeffModels/ReynoldsAnalogy/ReynoldsAnalogy.H
@@ -131,7 +131,11 @@ protected:
         tmp<FieldField<Field, scalar>> Cf() const;
 
         //- Set the heat transfer coefficient
-        virtual void htc(volScalarField& htc);
+        virtual void htc
+        (
+            volScalarField& htc,
+            const FieldField<Field, scalar>& q
+        );
 
 
         //- No copy construct
diff --git a/src/functionObjects/field/heatTransferCoeff/heatTransferCoeffModels/fixedReferenceTemperature/fixedReferenceTemperature.C b/src/functionObjects/field/heatTransferCoeff/heatTransferCoeffModels/fixedReferenceTemperature/fixedReferenceTemperature.C
index d2bd89459aa596930477cb296e4e6b87733ba922..e7f8b7a72ea54539502b0aadd6b4e360bb23a40d 100644
--- a/src/functionObjects/field/heatTransferCoeff/heatTransferCoeffModels/fixedReferenceTemperature/fixedReferenceTemperature.C
+++ b/src/functionObjects/field/heatTransferCoeff/heatTransferCoeffModels/fixedReferenceTemperature/fixedReferenceTemperature.C
@@ -81,10 +81,11 @@ bool Foam::heatTransferCoeffModels::fixedReferenceTemperature::read
 
 void Foam::heatTransferCoeffModels::fixedReferenceTemperature::htc
 (
-    volScalarField& htc
+    volScalarField& htc,
+    const FieldField<Field, scalar>& q
 )
 {
-    const FieldField<Field, scalar> qBf(q());
+    //const FieldField<Field, scalar> qBf(q());
     const volScalarField& T = mesh_.lookupObject<volScalarField>(TName_);
     const volScalarField::Boundary& Tbf = T.boundaryField();
     const scalar eps = ROOTVSMALL;
@@ -92,7 +93,7 @@ void Foam::heatTransferCoeffModels::fixedReferenceTemperature::htc
     volScalarField::Boundary& htcBf = htc.boundaryFieldRef();
     for (const label patchi : patchSet_)
     {
-        htcBf[patchi] = qBf[patchi]/(TRef_ - Tbf[patchi] + eps);
+        htcBf[patchi] = q[patchi]/(TRef_ - Tbf[patchi] + eps);
     }
 }
 
diff --git a/src/functionObjects/field/heatTransferCoeff/heatTransferCoeffModels/fixedReferenceTemperature/fixedReferenceTemperature.H b/src/functionObjects/field/heatTransferCoeff/heatTransferCoeffModels/fixedReferenceTemperature/fixedReferenceTemperature.H
index e5f17ea8da0547b0404c10733f7bf936ee8003e0..589c48b7a3f080af7065eba619abe56800ebea4f 100644
--- a/src/functionObjects/field/heatTransferCoeff/heatTransferCoeffModels/fixedReferenceTemperature/fixedReferenceTemperature.H
+++ b/src/functionObjects/field/heatTransferCoeff/heatTransferCoeffModels/fixedReferenceTemperature/fixedReferenceTemperature.H
@@ -94,7 +94,11 @@ protected:
     // Protected Member Functions
 
         //- Set the heat transfer coefficient
-        virtual void htc(volScalarField& htc);
+        virtual void htc
+        (
+            volScalarField& htc,
+            const FieldField<Field, scalar>& q
+        );
 
         //- No copy construct
         fixedReferenceTemperature(const fixedReferenceTemperature&) = delete;
diff --git a/src/functionObjects/field/heatTransferCoeff/heatTransferCoeffModels/heatTransferCoeffModel/heatTransferCoeffModel.C b/src/functionObjects/field/heatTransferCoeff/heatTransferCoeffModels/heatTransferCoeffModel/heatTransferCoeffModel.C
index 0e078a8d784f11bb210b4bfb963b72041a564978..4fefa38596b064163cfc53e1b64f138e17c57301 100644
--- a/src/functionObjects/field/heatTransferCoeff/heatTransferCoeffModels/heatTransferCoeffModel/heatTransferCoeffModel.C
+++ b/src/functionObjects/field/heatTransferCoeff/heatTransferCoeffModels/heatTransferCoeffModel/heatTransferCoeffModel.C
@@ -30,7 +30,6 @@ License
 #include "fluidThermo.H"
 #include "turbulentTransportModel.H"
 #include "turbulentFluidThermoModel.H"
-#include "phaseSystem.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
@@ -91,32 +90,10 @@ Foam::heatTransferCoeffModel::q() const
             q[patchi] = alphabf[patchi]*hebf[patchi].snGrad();
         }
     }
-    else if (mesh_.foundObject<phaseSystem>("phaseProperties"))
-    {
-        const phaseSystem& fluid =
-        (
-            mesh_.lookupObject<phaseSystem>("phaseProperties")
-        );
-
-        for (label patchi : patchSet_)
-        {
-            forAll(fluid.phases(), phasei)
-            {
-                const phaseModel& phase = fluid.phases()[phasei];
-                const fvPatchScalarField& alpha =
-                    phase.boundaryField()[patchi];
-                const volScalarField& he = phase.thermo().he();
-                const volScalarField::Boundary& hebf = he.boundaryField();
-
-                q[patchi] +=
-                    alpha*phase.alphaEff(patchi)()*hebf[patchi].snGrad();
-            }
-        }
-    }
     else
     {
         FatalErrorInFunction
-            << "Unable to find a valid thermo model to evaluate q"
+            << "Unable to find a valid thermo model to evaluate q" << nl
             << exit(FatalError);
     }
 
@@ -158,11 +135,7 @@ Foam::heatTransferCoeffModel::heatTransferCoeffModel
 
 bool Foam::heatTransferCoeffModel::read(const dictionary& dict)
 {
-    patchSet_ =
-        mesh_.boundaryMesh().patchSet
-        (
-            dict.get<wordRes>("patches")
-        );
+    patchSet_ = mesh_.boundaryMesh().patchSet(dict.get<wordRes>("patches"));
 
     dict.readIfPresent("qr", qrName_);
 
@@ -170,9 +143,13 @@ bool Foam::heatTransferCoeffModel::read(const dictionary& dict)
 }
 
 
-bool Foam::heatTransferCoeffModel::calc(volScalarField& result)
+bool Foam::heatTransferCoeffModel::calc
+(
+    volScalarField& result,
+    const FieldField<Field, scalar>& q
+)
 {
-    htc(result);
+    htc(result, q);
 
     return true;
 }
diff --git a/src/functionObjects/field/heatTransferCoeff/heatTransferCoeffModels/heatTransferCoeffModel/heatTransferCoeffModel.H b/src/functionObjects/field/heatTransferCoeff/heatTransferCoeffModels/heatTransferCoeffModel/heatTransferCoeffModel.H
index d9e9ec13b6f8cf5db686ddd4a78c6bc0c851624b..ebeaf24da60c6101791e0e7c93c722fd75e83e9b 100644
--- a/src/functionObjects/field/heatTransferCoeff/heatTransferCoeffModels/heatTransferCoeffModel/heatTransferCoeffModel.H
+++ b/src/functionObjects/field/heatTransferCoeff/heatTransferCoeffModels/heatTransferCoeffModel/heatTransferCoeffModel.H
@@ -47,7 +47,6 @@ SourceFiles
 #include "dictionary.H"
 #include "HashSet.H"
 #include "volFields.H"
-
 #include "runTimeSelectionTables.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@@ -66,7 +65,7 @@ class heatTransferCoeffModel
 {
 protected:
 
-    // Protected Data
+    // Public Data
 
         //- Mesh reference
         const fvMesh& mesh_;
@@ -77,16 +76,20 @@ protected:
         //- Temperature name
         const word TName_;
 
-        //- Name of radiative heat flux name, default = qr
+        //- Name of radiative heat flux (default = qr)
         word qrName_;
 
 
-    // Protected Member Functions
+protected:
 
-        tmp<FieldField<Field, scalar>> q() const;
+    // Protected Member Functions
 
         //- Set the heat transfer coefficient
-        virtual void htc(volScalarField& htc) = 0;
+        virtual void htc
+        (
+            volScalarField& htc,
+            const FieldField<Field, scalar>& q
+        ) = 0;
 
         //- No copy construct
         heatTransferCoeffModel(const heatTransferCoeffModel&) = delete;
@@ -145,10 +148,42 @@ public:
 
     // Member Functions
 
+        //- The mesh reference
+        const fvMesh& mesh() const
+        {
+            return mesh_;
+        }
+
+        //- Wall patches to process
+        const labelHashSet& patchSet() const
+        {
+            return patchSet_;
+        }
+
+        //- Temperature name
+        const word& TName() const
+        {
+            return TName_;
+        }
+
+        //- Name of radiative heat flux
+        const word& qrName() const
+        {
+            return qrName_;
+        }
+
+
         //- Read from dictionary
         virtual bool read(const dictionary& dict);
 
-        virtual bool calc(volScalarField& result);
+        virtual bool calc
+        (
+            volScalarField& result,
+            const FieldField<Field, scalar>& q
+        );
+
+        //- Return q boundary fields
+        tmp<FieldField<Field, scalar>> q() const;
 };
 
 
diff --git a/src/functionObjects/field/heatTransferCoeff/heatTransferCoeffModels/localReferenceTemperature/localReferenceTemperature.C b/src/functionObjects/field/heatTransferCoeff/heatTransferCoeffModels/localReferenceTemperature/localReferenceTemperature.C
index b26f01502eb7c5e96287984b2ede80bbe227ef13..1492d33c4a63b1918941f28ec146ce0d16c2fb0c 100644
--- a/src/functionObjects/field/heatTransferCoeff/heatTransferCoeffModels/localReferenceTemperature/localReferenceTemperature.C
+++ b/src/functionObjects/field/heatTransferCoeff/heatTransferCoeffModels/localReferenceTemperature/localReferenceTemperature.C
@@ -74,10 +74,11 @@ bool Foam::heatTransferCoeffModels::localReferenceTemperature::read
 
 void Foam::heatTransferCoeffModels::localReferenceTemperature::htc
 (
-    volScalarField& htc
+    volScalarField& htc,
+    const FieldField<Field, scalar>& q
 )
 {
-    const FieldField<Field, scalar> qBf(q());
+    ///const FieldField<Field, scalar> qBf(q());
     const volScalarField& T = mesh_.lookupObject<volScalarField>(TName_);
     const volScalarField::Boundary& Tbf = T.boundaryField();
     const scalar eps = ROOTVSMALL;
@@ -87,7 +88,7 @@ void Foam::heatTransferCoeffModels::localReferenceTemperature::htc
     for (const label patchi : patchSet_)
     {
         const scalarField Tc(Tbf[patchi].patchInternalField());
-        htcBf[patchi] = qBf[patchi]/(Tc - Tbf[patchi] + eps);
+        htcBf[patchi] = q[patchi]/(Tc - Tbf[patchi] + eps);
     }
 }
 
diff --git a/src/functionObjects/field/heatTransferCoeff/heatTransferCoeffModels/localReferenceTemperature/localReferenceTemperature.H b/src/functionObjects/field/heatTransferCoeff/heatTransferCoeffModels/localReferenceTemperature/localReferenceTemperature.H
index f7578eec4d98b09280ace3eedce0fb4147b1ad46..7175b2091664630c1566615e11f4a3f208ed3379 100644
--- a/src/functionObjects/field/heatTransferCoeff/heatTransferCoeffModels/localReferenceTemperature/localReferenceTemperature.H
+++ b/src/functionObjects/field/heatTransferCoeff/heatTransferCoeffModels/localReferenceTemperature/localReferenceTemperature.H
@@ -81,7 +81,16 @@ class localReferenceTemperature
 :
     public heatTransferCoeffModel
 {
-    // Private Member Functions
+protected:
+
+    // Protected Member Functions
+
+        //- Set the heat transfer coefficient
+        virtual void htc
+        (
+            volScalarField& htc,
+            const FieldField<Field, scalar>& q
+        );
 
         //- No copy construct
         localReferenceTemperature(const localReferenceTemperature&) = delete;
@@ -90,14 +99,6 @@ class localReferenceTemperature
         void operator=(const localReferenceTemperature&) = delete;
 
 
-protected:
-
-    // Protected Member Functions
-
-        //- Set the heat transfer coefficient
-        virtual void htc(volScalarField& htc);
-
-
 public:
 
     //- Runtime type information
diff --git a/src/phaseSystemModels/Allwclean b/src/phaseSystemModels/Allwclean
index ab4d65ca0355012639caaaac67fb7badde8bea01..cc8f941470c33b473e1a9edcfaa31965b31e6fe8 100755
--- a/src/phaseSystemModels/Allwclean
+++ b/src/phaseSystemModels/Allwclean
@@ -11,5 +11,6 @@ wclean libso reactingEulerFoam/reactingMultiphaseEulerFoam/multiphaseCompressibl
 
 wclean libso reactingEulerFoam/reactingTwoPhaseEulerFoam/twoPhaseSystem
 wclean libso reactingEulerFoam/reactingTwoPhaseEulerFoam/twoPhaseCompressibleTurbulenceModels
+wclean libso  reactingEulerFoam/functionObjects
 
 #------------------------------------------------------------------------------
diff --git a/src/phaseSystemModels/Allwmake b/src/phaseSystemModels/Allwmake
index 686604a701fcdf2f2aea66116ec7b4f736e64608..65eb00e8ee7bdfeebf1c81f950a84905a7309a8c 100755
--- a/src/phaseSystemModels/Allwmake
+++ b/src/phaseSystemModels/Allwmake
@@ -27,5 +27,6 @@ wmake $targetType  reactingEulerFoam/reactingTwoPhaseEulerFoam/twoPhaseSystem
 wmake $targetType  reactingEulerFoam/reactingTwoPhaseEulerFoam/twoPhaseCompressibleTurbulenceModels
 
 wmake $targetType  reactingEulerFoam/phaseSystems
+wmake $targetType  reactingEulerFoam/functionObjects
 
 #------------------------------------------------------------------------------
diff --git a/applications/solvers/multiphase/reactingEulerFoam/functionObjects/Make/files b/src/phaseSystemModels/reactingEulerFoam/functionObjects/Make/files
similarity index 72%
rename from applications/solvers/multiphase/reactingEulerFoam/functionObjects/Make/files
rename to src/phaseSystemModels/reactingEulerFoam/functionObjects/Make/files
index ca7d02162dbd5f5feb601c39c9fae101549aa696..567c9619ea0ba90bf3f7aef88cb3150f079ee723 100644
--- a/applications/solvers/multiphase/reactingEulerFoam/functionObjects/Make/files
+++ b/src/phaseSystemModels/reactingEulerFoam/functionObjects/Make/files
@@ -1,4 +1,5 @@
 sizeDistribution/sizeDistribution.C
 phaseForces/phaseForces.C
+reactingEulerHtcModel/reactingEulerHtcModel.C
 
 LIB = $(FOAM_LIBBIN)/libreactingEulerFoamFunctionObjects
diff --git a/applications/solvers/multiphase/reactingEulerFoam/functionObjects/Make/options b/src/phaseSystemModels/reactingEulerFoam/functionObjects/Make/options
similarity index 100%
rename from applications/solvers/multiphase/reactingEulerFoam/functionObjects/Make/options
rename to src/phaseSystemModels/reactingEulerFoam/functionObjects/Make/options
diff --git a/applications/solvers/multiphase/reactingEulerFoam/functionObjects/phaseForces/phaseForces.C b/src/phaseSystemModels/reactingEulerFoam/functionObjects/phaseForces/phaseForces.C
similarity index 100%
rename from applications/solvers/multiphase/reactingEulerFoam/functionObjects/phaseForces/phaseForces.C
rename to src/phaseSystemModels/reactingEulerFoam/functionObjects/phaseForces/phaseForces.C
diff --git a/applications/solvers/multiphase/reactingEulerFoam/functionObjects/phaseForces/phaseForces.H b/src/phaseSystemModels/reactingEulerFoam/functionObjects/phaseForces/phaseForces.H
similarity index 100%
rename from applications/solvers/multiphase/reactingEulerFoam/functionObjects/phaseForces/phaseForces.H
rename to src/phaseSystemModels/reactingEulerFoam/functionObjects/phaseForces/phaseForces.H
diff --git a/src/phaseSystemModels/reactingEulerFoam/functionObjects/reactingEulerHtcModel/reactingEulerHtcModel.C b/src/phaseSystemModels/reactingEulerFoam/functionObjects/reactingEulerHtcModel/reactingEulerHtcModel.C
new file mode 100644
index 0000000000000000000000000000000000000000..85eb2500ba5410ab683d9121457318eaf3c9aafd
--- /dev/null
+++ b/src/phaseSystemModels/reactingEulerFoam/functionObjects/reactingEulerHtcModel/reactingEulerHtcModel.C
@@ -0,0 +1,178 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2020 OpenCFD Ltd.
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "reactingEulerHtcModel.H"
+#include "phaseSystem.H"
+#include "addToRunTimeSelectionTable.H"
+#include "dictionary.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace functionObjects
+{
+    defineTypeNameAndDebug(reactingEulerHtcModel, 0);
+    addToRunTimeSelectionTable
+    (
+        functionObject,
+        reactingEulerHtcModel,
+        dictionary
+    );
+}
+}
+
+
+// * * * * * * * * * * * * Protected Member Functions  * * * * * * * * * * * //
+
+Foam::tmp<Foam::FieldField<Foam::Field, Foam::scalar>>
+Foam::functionObjects::reactingEulerHtcModel::q() const
+{
+    const fvMesh& mesh = htcModelPtr_->mesh();
+
+    const volScalarField& T =
+        mesh.lookupObject<volScalarField>(htcModelPtr_->TName());
+
+    const volScalarField::Boundary& Tbf = T.boundaryField();
+
+    auto tq = tmp<FieldField<Field, scalar>>::New(Tbf.size());
+    auto& q = tq.ref();
+
+    forAll(q, patchi)
+    {
+        q.set(patchi, new Field<scalar>(Tbf[patchi].size(), Zero));
+    }
+
+    const auto* fluidPtr =
+        mesh.cfindObject<phaseSystem>("phaseProperties");
+
+    if (!fluidPtr)
+    {
+        FatalErrorInFunction
+            << "Unable to find a valid phaseSystem to evaluate q" << nl
+            << exit(FatalError);
+    }
+
+    const phaseSystem& fluid = *fluidPtr;
+
+    for (const label patchi : htcModelPtr_->patchSet())
+    {
+        for (const phaseModel& phase : fluid.phases())
+        {
+            const fvPatchScalarField& alpha = phase.boundaryField()[patchi];
+            const volScalarField& he = phase.thermo().he();
+            const volScalarField::Boundary& hebf = he.boundaryField();
+
+            q[patchi] +=
+                alpha*phase.alphaEff(patchi)()*hebf[patchi].snGrad();
+        }
+    }
+
+    // Add radiative heat flux contribution if present
+
+    const volScalarField* qrPtr =
+        mesh.cfindObject<volScalarField>(htcModelPtr_->qrName());
+
+    if (qrPtr)
+    {
+        const volScalarField::Boundary& qrbf = qrPtr->boundaryField();
+
+        for (const label patchi : htcModelPtr_->patchSet())
+        {
+            q[patchi] += qrbf[patchi];
+        }
+    }
+
+    return tq;
+}
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+bool Foam::functionObjects::reactingEulerHtcModel::calc()
+{
+    auto& htc =
+        htcModelPtr_->mesh().lookupObjectRef<volScalarField>(resultName_);
+
+    htcModelPtr_->calc(htc, q());
+
+    return true;
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::functionObjects::reactingEulerHtcModel::reactingEulerHtcModel
+(
+    const word& name,
+    const Time& runTime,
+    const dictionary& dict
+)
+:
+    fieldExpression(name, runTime, dict),
+    htcModelPtr_(nullptr)
+{
+    read(dict);
+
+    setResultName(typeName, name + ":htc:" + htcModelPtr_->type());
+
+    volScalarField* htcPtr =
+        new volScalarField
+        (
+            IOobject
+            (
+                resultName_,
+                mesh_.time().timeName(),
+                mesh_,
+                IOobject::NO_READ,
+                IOobject::NO_WRITE
+            ),
+            mesh_,
+            dimensionedScalar(dimPower/dimArea/dimTemperature, Zero)
+        );
+
+    mesh_.objectRegistry::store(htcPtr);
+}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+bool Foam::functionObjects::reactingEulerHtcModel::read(const dictionary& dict)
+{
+    if (fieldExpression::read(dict))
+    {
+        htcModelPtr_ = heatTransferCoeffModel::New(dict, mesh_, fieldName_);
+
+        htcModelPtr_->read(dict);
+
+        return true;
+    }
+
+    return false;
+}
+
+
+// ************************************************************************* //
diff --git a/src/phaseSystemModels/reactingEulerFoam/functionObjects/reactingEulerHtcModel/reactingEulerHtcModel.H b/src/phaseSystemModels/reactingEulerFoam/functionObjects/reactingEulerHtcModel/reactingEulerHtcModel.H
new file mode 100644
index 0000000000000000000000000000000000000000..7c1536bf7d8097f9301ed25228043170ee991cf1
--- /dev/null
+++ b/src/phaseSystemModels/reactingEulerFoam/functionObjects/reactingEulerHtcModel/reactingEulerHtcModel.H
@@ -0,0 +1,122 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2020 OpenCFD Ltd.
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::reactingEulerHtcModel
+
+Description
+    A heat transfer coefficient for reactingEuler solvers
+
+SourceFiles
+    reactingEulerHtcModel.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef reactingEulerHtcModel_H
+#define reactingEulerHtcModel_H
+
+#include "HashSet.H"
+#include "volFields.H"
+#include "fieldExpression.H"
+#include "runTimeSelectionTables.H"
+#include "heatTransferCoeffModel.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace functionObjects
+{
+
+/*---------------------------------------------------------------------------*\
+                   Class reactingEulerHtcModel Declaration
+\*---------------------------------------------------------------------------*/
+
+class reactingEulerHtcModel
+:
+    public fieldExpression
+{
+    // Private Data
+
+        //- Heat transfer coefficient model
+        autoPtr<heatTransferCoeffModel> htcModelPtr_;
+
+
+protected:
+
+    // Protected Member Functions
+
+        //- Calculate the heat transfer coefficient field
+        //  \return true on success
+        virtual bool calc();
+
+        //- Calculate heat flux
+        tmp<FieldField<Field, scalar>> q() const;
+
+        //- No copy construct
+        reactingEulerHtcModel(const reactingEulerHtcModel&) = delete;
+
+        //- No copy assignment
+        void operator=(const reactingEulerHtcModel&) = delete;
+
+
+public:
+
+    //- Runtime type information
+    TypeName("reactingEulerHtcModel");
+
+
+    // Constructors
+
+        //- Construct from components
+        reactingEulerHtcModel
+        (
+            const word& name,
+            const Time& runTime,
+            const dictionary& dict
+        );
+
+
+    //- Destructor
+    virtual ~reactingEulerHtcModel() = default;
+
+
+    // Member Functions
+
+        //- Read the heatTransferCoeff data
+        virtual bool read(const dictionary& dict);
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace functionObjects
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/applications/solvers/multiphase/reactingEulerFoam/functionObjects/sizeDistribution/sizeDistribution.C b/src/phaseSystemModels/reactingEulerFoam/functionObjects/sizeDistribution/sizeDistribution.C
similarity index 100%
rename from applications/solvers/multiphase/reactingEulerFoam/functionObjects/sizeDistribution/sizeDistribution.C
rename to src/phaseSystemModels/reactingEulerFoam/functionObjects/sizeDistribution/sizeDistribution.C
diff --git a/applications/solvers/multiphase/reactingEulerFoam/functionObjects/sizeDistribution/sizeDistribution.H b/src/phaseSystemModels/reactingEulerFoam/functionObjects/sizeDistribution/sizeDistribution.H
similarity index 100%
rename from applications/solvers/multiphase/reactingEulerFoam/functionObjects/sizeDistribution/sizeDistribution.H
rename to src/phaseSystemModels/reactingEulerFoam/functionObjects/sizeDistribution/sizeDistribution.H
diff --git a/tutorials/heatTransfer/buoyantPimpleFoam/hotRoom/system/controlDict b/tutorials/heatTransfer/buoyantPimpleFoam/hotRoom/system/controlDict
index eb8dd8a8e6405de9a2fda71346753983372f1b68..13683da17deedf67a6991cf80f5b3135af6f74bf 100644
--- a/tutorials/heatTransfer/buoyantPimpleFoam/hotRoom/system/controlDict
+++ b/tutorials/heatTransfer/buoyantPimpleFoam/hotRoom/system/controlDict
@@ -49,4 +49,19 @@ adjustTimeStep  no;
 
 maxCo           0.5;
 
+functions
+{
+    htc
+    {
+        type            heatTransferCoeff;
+        libs            (fieldFunctionObjects);
+        field           T;
+        writeControl    outputTime;
+        writeInterval   1;
+        htcModel        fixedReferenceTemperature;
+        patches         (ceiling);
+        TRef            373;
+    }
+}
+
 // ************************************************************************* //
diff --git a/tutorials/heatTransfer/chtMultiRegionTwoPhaseEulerFoam/solidQuenching2D/system/controlDict b/tutorials/heatTransfer/chtMultiRegionTwoPhaseEulerFoam/solidQuenching2D/system/controlDict
index 9434d196ca6415a0bdf64953cc9176f042e0ccf6..62866a741d662b5d6e34539b7d1a7ede0cbddc9f 100644
--- a/tutorials/heatTransfer/chtMultiRegionTwoPhaseEulerFoam/solidQuenching2D/system/controlDict
+++ b/tutorials/heatTransfer/chtMultiRegionTwoPhaseEulerFoam/solidQuenching2D/system/controlDict
@@ -63,8 +63,8 @@ functions
 {
     htc
     {
-        type            heatTransferCoeff;
-        libs            (fieldFunctionObjects);
+        type            reactingEulerHtcModel;
+        libs            (reactingEulerFoamFunctionObjects);
         region          water;
         field           T.liquid;
         writeControl    outputTime;
@@ -94,7 +94,7 @@ functions
 
         operation       areaAverage;
 
-        fields          ( heatTransferCoeff(T.liquid) );
+        fields          ( reactingEulerHtcModel(T.liquid) );
     }
 }