multiphaseEuler: Updated to us the new phaseSolidThermophysicalTransportModel class

for thermophysical transport within stationary solid phases.  This provides a
consistent interface to heat transport within solids for single and now
multiphase solvers so that for example the wallHeatFlux functionObject can now
be used with multiphaseEuler, see tutorials/multiphaseEuler/boilingBed.
Also this development supports anisotropic thermal conductivity within the
stationary solid regions which was not possible previously.

The tutorials/multiphaseEuler/bed and tutorials/multiphaseEuler/boilingBed
tutorial cases have been updated for phaseSolidThermophysicalTransportModel by
changing the thermo type in physicalProperties.solid to heSolidThermo.  This
change will need to be made to all multiphaseEuler cases involving stationary
phases.

applications/modules/multiphaseEuler/cellPressureCorrector.C

@@ -498,7 +498,10 @@ void Foam::solvers::multiphaseEuler::cellPressureCorrector()
         forAll(phases, phasei)
         {
             phaseModel& phase = phases_[phasei];
-            phase.rho() += phase.thermo().psi()*(p_rgh - p_rgh_0);
+            if (!phase.incompressible())
+            {
+                phase.rho() += phase.fluidThermo().psi()*(p_rgh - p_rgh_0);
+            }
         }
 
         // Update mass transfer rates for change in p_rgh

applications/modules/multiphaseEuler/compressibilityEqns.C

@@ -75,7 +75,7 @@ Foam::solvers::multiphaseEuler::compressibilityEqns
                 const surfaceScalarField phid
                 (
                     IOobject::groupName("phid", phase.name()),
-                    fvc::interpolate(phase.thermo().psi())*phase.phi()
+                    fvc::interpolate(phase.fluidThermo().psi())*phase.phi()
                 );
 
                 pEqnComp +=
@@ -83,7 +83,7 @@ Foam::solvers::multiphaseEuler::compressibilityEqns
                     (
                         (alpha/rho)*
                         (
-                            phase.thermo().psi()*fvm::ddt(p_rgh)
+                            phase.fluidThermo().psi()*fvm::ddt(p_rgh)
                           + fvm::div(phid, p_rgh)
                           - fvm::Sp(fvc::div(phid), p_rgh)
                         )
@@ -94,7 +94,7 @@ Foam::solvers::multiphaseEuler::compressibilityEqns
             else
             {
                 pEqnComp +=
-                    (alpha*phase.thermo().psi()/rho)
+                    (alpha*phase.fluidThermo().psi()/rho)
                    *correction(fvm::ddt(p_rgh));
             }
         }

applications/modules/multiphaseEuler/facePressureCorrector.C

@@ -387,7 +387,10 @@ void Foam::solvers::multiphaseEuler::facePressureCorrector()
         forAll(phases, phasei)
         {
             phaseModel& phase = phases_[phasei];
-            phase.rho() += phase.thermo().psi()*(p_rgh - p_rgh_0);
+            if (!phase.incompressible())
+            {
+                phase.rho() += phase.fluidThermo().psi()*(p_rgh - p_rgh_0);
+            }
         }
 
         // Update mass transfer rates for change in p_rgh

applications/modules/multiphaseEuler/fvModels/interfaceTurbulenceDamping/interfaceTurbulenceDamping.C

@@ -145,13 +145,14 @@ void Foam::fv::interfaceTurbulenceDamping::addRhoSup
 
     volScalarField::Internal aSqrnu
     (
-        movingPhases[0]*sqr(movingPhases[0].thermo().nu()()())
+        movingPhases[0]*sqr(movingPhases[0].fluidThermo().nu()()())
     );
 
     for (label phasei=1; phasei<movingPhases.size(); phasei++)
     {
         aSqrnu +=
-            movingPhases[phasei]*sqr(movingPhases[phasei].thermo().nu()()());
+            movingPhases[phasei]
+           *sqr(movingPhases[phasei].fluidThermo().nu()()());
     }
 
     if (field.name() == "epsilon")
@@ -275,7 +276,7 @@ void Foam::fv::interfaceTurbulenceDamping::addSup
 
     const volScalarField::Internal aSqrnu
     (
-        alpha*sqr(phase_.thermo().nu()()())
+        alpha*sqr(phase_.fluidThermo().nu()()())
     );
 
     if (field.name() == IOobject::groupName("epsilon", phaseName_))

applications/modules/multiphaseEuler/interfacialModels/dragModels/AttouFerschneider/AttouFerschneider.C

@@ -58,7 +58,7 @@ Foam::dragModels::AttouFerschneider::KGasLiquid
     const volScalarField magURel(mag(gas.U() - liquid.U()));
 
     return
-        E2_*gas.thermo().mu()*sqr(oneMinusGas/solid.d())*sqr(cbrtR)
+        E2_*gas.fluidThermo().mu()*sqr(oneMinusGas/solid.d())*sqr(cbrtR)
        /max(gas, gas.residualAlpha())
       + E2_*gas.rho()*magURel*(1 - gas)/solid.d()*cbrtR;
 }
@@ -78,7 +78,7 @@ Foam::dragModels::AttouFerschneider::KGasSolid
     );
 
     return
-        E1_*gas.thermo().mu()*sqr(oneMinusGas/solid.d())*sqr(cbrtR)
+        E1_*gas.fluidThermo().mu()*sqr(oneMinusGas/solid.d())*sqr(cbrtR)
        /max(gas, gas.residualAlpha())
       + E2_*gas.rho()*mag(gas.U())*(1 - gas)/solid.d()*cbrtR;
 }
@@ -94,7 +94,7 @@ Foam::dragModels::AttouFerschneider::KLiquidSolid
     const phaseModel& gas = liquid.fluid().phases()[gasName_];
 
     return
-        E1_*liquid.thermo().mu()
+        E1_*liquid.fluidThermo().mu()
        *sqr(max(solid, solid.residualAlpha())/solid.d())
        /max(liquid, liquid.residualAlpha())
       + E2_*liquid.rho()*mag(gas.U())*solid/solid.d();

applications/modules/multiphaseEuler/interfacialModels/dragModels/IshiiZuber/IshiiZuber.C

@@ -65,8 +65,8 @@ Foam::dragModels::IshiiZuber::CdRe() const
     const volScalarField Re(interface_.Re());
     const volScalarField Eo(interface_.Eo());
 
-    const volScalarField mud(interface_.dispersed().thermo().mu());
+    const volScalarField mud(interface_.dispersed().fluidThermo().mu());
-    const volScalarField muc(interface_.continuous().thermo().mu());
+    const volScalarField muc(interface_.continuous().fluidThermo().mu());
 
     const volScalarField muStar((mud + 0.4*muc)/(mud + muc));
 

applications/modules/multiphaseEuler/interfacialModels/dragModels/aerosolDrag/aerosolDrag.C

@@ -71,7 +71,7 @@ Foam::dragModels::aerosolDrag::~aerosolDrag()
 Foam::tmp<Foam::volScalarField> Foam::dragModels::aerosolDrag::CdRe() const
 {
     const volScalarField& T = interface_.continuous().thermo().T();
-    const volScalarField& p = interface_.continuous().thermo().p();
+    const volScalarField& p = interface_.continuous().fluidThermo().p();
     tmp<volScalarField> td(interface_.dispersed().d());
     const volScalarField& d = td();
 

applications/modules/multiphaseEuler/interfacialModels/dragModels/dispersedDragModel/dispersedDragModel.C

@@ -62,7 +62,7 @@ Foam::tmp<Foam::volScalarField> Foam::dragModels::dispersedDragModel::Ki() const
        *CdRe()
        *swarmCorrection_->Cs()
        *interface_.continuous().rho()
-       *interface_.continuous().thermo().nu()
+       *interface_.continuous().fluidThermo().nu()
        /sqr(interface_.dispersed().d());
 }
 

applications/modules/multiphaseEuler/interfacialModels/dragModels/segregated/segregated.C

@@ -75,8 +75,8 @@ Foam::tmp<Foam::volScalarField> Foam::dragModels::segregated::K() const
     const volScalarField::Internal& rho1(interface_.phase1().rho());
     const volScalarField::Internal& rho2(interface_.phase2().rho());
 
-    tmp<volScalarField> tnu1(interface_.phase1().thermo().nu());
+    tmp<volScalarField> tnu1(interface_.phase1().fluidThermo().nu());
-    tmp<volScalarField> tnu2(interface_.phase2().thermo().nu());
+    tmp<volScalarField> tnu2(interface_.phase2().fluidThermo().nu());
 
     const volScalarField::Internal& nu1(tnu1());
     const volScalarField::Internal& nu2(tnu2());

applications/modules/multiphaseEuler/interfacialModels/dragModels/timeScaleFilteredDrag/timeScaleFilteredDrag.C

@@ -83,7 +83,7 @@ Foam::dragModels::timeScaleFilteredDrag::CdRe() const
        /0.75
        /swarmCorrection_->Cs()
        /interface_.continuous().rho()
-       /interface_.continuous().thermo().nu()
+       /interface_.continuous().fluidThermo().nu()
        /minRelaxTime_
     );
 

applications/modules/multiphaseEuler/interfacialModels/interfaceCompositionModels/interfaceCompositionModel/interfaceCompositionModel.H

@@ -71,7 +71,7 @@ class interfaceCompositionModel
         const rhoFluidMulticomponentThermo& thermo_;
 
         //- General thermo model for the other side of the interface
-        const rhoFluidThermo& otherThermo_;
+        const rhoThermo& otherThermo_;
 
 
 public:
@@ -133,7 +133,7 @@ public:
             inline const rhoFluidMulticomponentThermo& thermo() const;
 
             //- Return the other thermo
-            inline const rhoFluidThermo& otherThermo() const;
+            inline const rhoThermo& otherThermo() const;
 
             //- Return the other multicomponent thermo
             inline const rhoFluidMulticomponentThermo&

applications/modules/multiphaseEuler/interfacialModels/interfaceCompositionModels/interfaceCompositionModel/interfaceCompositionModelI.H

@@ -48,7 +48,7 @@ Foam::interfaceCompositionModel::thermo() const
 }
 
 
-const Foam::rhoFluidThermo& Foam::interfaceCompositionModel::otherThermo() const
+const Foam::rhoThermo& Foam::interfaceCompositionModel::otherThermo() const
 {
     return otherThermo_;
 }

applications/modules/multiphaseEuler/interfacialModels/liftModels/LegendreMagnaudet/LegendreMagnaudet.C

@@ -67,7 +67,7 @@ Foam::tmp<Foam::volScalarField> Foam::liftModels::LegendreMagnaudet::Cl() const
     const volScalarField Sr
     (
         sqr(interface_.dispersed().d())
-       /(Re*interface_.continuous().thermo().nu())
+       /(Re*interface_.continuous().fluidThermo().nu())
        *mag(fvc::grad(interface_.continuous().U()))
     );
 

applications/modules/multiphaseEuler/interfacialModels/liftModels/Moraga/Moraga.C

@@ -66,7 +66,7 @@ Foam::tmp<Foam::volScalarField> Foam::liftModels::Moraga::Cl() const
     volScalarField sqrSr
     (
         sqr(interface_.dispersed().d())
-       /interface_.continuous().thermo().nu()
+       /interface_.continuous().fluidThermo().nu()
        *mag(fvc::grad(interface_.continuous().U()))
     );
 

applications/modules/multiphaseEuler/interfacialModels/liftModels/SaffmanMei/SaffmanMei.C

@@ -69,7 +69,7 @@ Foam::tmp<Foam::volScalarField> Foam::liftModels::SaffmanMei::Cl() const
         mag(fvc::curl(interface_.continuous().U()))
        *sqr(interface_.dispersed().d())
        /(
-            interface_.continuous().thermo().nu()
+            interface_.continuous().fluidThermo().nu()
           + dimensionedScalar(dimViscosity, small)
         )
     );

applications/modules/multiphaseEuler/interfacialModels/phaseTransferModels/cavitation/cavitation.C

@@ -87,7 +87,7 @@ Foam::phaseTransferModels::cavitation::dmdtf() const
 
     const Pair<tmp<volScalarField::Internal>> coeffs(cavitation_->mDot12P());
 
-    const volScalarField::Internal& p = interface_.phase1().thermo().p();
+    const volScalarField::Internal& p = interface_.phase1().fluidThermo().p();
     const volScalarField::Internal pSat1(cavitation_->pSat1());
     const volScalarField::Internal pSat2(cavitation_->pSat2());
 

applications/modules/multiphaseEuler/multiphaseEuler.C

@@ -169,7 +169,7 @@ Foam::solvers::multiphaseEuler::multiphaseEuler(fvMesh& mesh)
 
     phi_(fluid_.phi()),
 
-    p_(phases_[0].thermo().p()),
+    p_(phases_[0].fluidThermo().p()),
 
     p_rgh(buoyancy.p_rgh),
 

applications/modules/multiphaseEuler/phaseSystem/Make/options

@@ -7,6 +7,7 @@ EXE_INC = \
     -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
     -I$(LIB_SRC)/thermophysicalModels/specie/lnInclude \
     -I$(LIB_SRC)/thermophysicalModels/multicomponentThermo/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/solidThermo/lnInclude \
     -I$(LIB_SRC)/thermophysicalModels/saturationModels/lnInclude \
     -I$(LIB_SRC)/twoPhaseModels/twoPhaseMixture/lnInclude \
     -I$(LIB_SRC)/twoPhaseModels/compressibleTwoPhases/lnInclude \
@@ -16,6 +17,8 @@ EXE_INC = \
     -I$(LIB_SRC)/MomentumTransportModels/phaseCompressible/lnInclude \
     -I$(LIB_SRC)/ThermophysicalTransportModels/thermophysicalTransportModel/lnInclude \
     -I$(LIB_SRC)/ThermophysicalTransportModels/fluid/lnInclude \
+    -I$(LIB_SRC)/ThermophysicalTransportModels/solid/lnInclude \
+    -I$(LIB_SRC)/ThermophysicalTransportModels/phaseSolid/lnInclude \
     -I$(LIB_SRC)/combustionModels/lnInclude
 
 LIB_LIBS = \
@@ -26,4 +29,5 @@ LIB_LIBS = \
     -lmultiphaseProperties \
     -lphaseCompressibleMomentumTransportModels \
     -lphaseFluidMulticomponentThermophysicalTransportModels \
+    -lphaseSolidThermophysicalTransportModels \
     -lcombustionModels

applications/modules/multiphaseEuler/phaseSystem/diameterModels/IATE/IATEsources/IATEsource/IATEsource.C

@@ -92,7 +92,7 @@ Foam::tmp<Foam::volScalarField> Foam::diameterModels::IATEsource::Ut() const
 
 Foam::tmp<Foam::volScalarField> Foam::diameterModels::IATEsource::Re() const
 {
-    return max(Ur()*phase().d()/otherPhase().thermo().nu(), scalar(1e-3));
+    return max(Ur()*phase().d()/otherPhase().fluidThermo().nu(), scalar(1e-3));
 }
 
 Foam::tmp<Foam::volScalarField> Foam::diameterModels::IATEsource::CD() const
@@ -118,7 +118,7 @@ Foam::tmp<Foam::volScalarField> Foam::diameterModels::IATEsource::Mo() const
         phase().db().lookupObject<uniformDimensionedVectorField>("g");
 
     return
-        mag(g)*pow4(otherPhase().thermo().nu())*sqr(otherPhase().rho())
+        mag(g)*pow4(otherPhase().fluidThermo().nu())*sqr(otherPhase().rho())
        *(otherPhase().rho() - phase().rho())
        /pow3(sigma());
 }

applications/modules/multiphaseEuler/phaseSystem/diameterModels/linearTsubDiameter/linearTsubDiameter.C

@@ -119,7 +119,7 @@ void Foam::diameterModels::linearTsub::correct()
 
         const volScalarField Tsub
         (
-            satModel.Tsat(liquid.thermo().p()) - liquid.thermo().T()
+            satModel.Tsat(liquid.fluidThermo().p()) - liquid.thermo().T()
         );
 
         d_ = max

applications/modules/multiphaseEuler/phaseSystem/phaseInterface/dispersedPhaseInterface/dispersedPhaseInterface.C

@@ -111,14 +111,14 @@ Foam::tmp<Foam::volVectorField> Foam::dispersedPhaseInterface::Ur() const
 
 Foam::tmp<Foam::volScalarField> Foam::dispersedPhaseInterface::Re() const
 {
-    return magUr()*dispersed().d()/continuous().thermo().nu();
+    return magUr()*dispersed().d()/continuous().fluidThermo().nu();
 }
 
 
 Foam::tmp<Foam::volScalarField> Foam::dispersedPhaseInterface::Pr() const
 {
     return
-         continuous().thermo().nu()
+         continuous().fluidThermo().nu()
         *continuous().thermo().Cp()
         *continuous().rho()
         /continuous().thermo().kappa();
@@ -148,10 +148,10 @@ Foam::tmp<Foam::volScalarField> Foam::dispersedPhaseInterface::Mo() const
 {
     return
         mag(g())
-       *continuous().thermo().nu()
+       *continuous().fluidThermo().nu()
        *pow3
         (
-            continuous().thermo().nu()
+            continuous().fluidThermo().nu()
            *continuous().rho()
            /sigma()
         );

applications/modules/multiphaseEuler/phaseSystem/phaseInterface/phaseInterface/phaseInterfaceI.H

@@ -41,13 +41,13 @@ inline const Foam::volScalarField& Foam::phaseInterface::alpha2() const
 
 inline const Foam::rhoFluidThermo& Foam::phaseInterface::thermo1() const
 {
-    return phase1().thermo();
+    return phase1().fluidThermo();
 }
 
 
 inline const Foam::rhoFluidThermo& Foam::phaseInterface::thermo2() const
 {
-    return phase2().thermo();
+    return phase2().fluidThermo();
 }
 
 

applications/modules/multiphaseEuler/phaseSystem/phaseModel/phaseModel.H

@@ -255,10 +255,16 @@ public:
         // Thermo
 
             //- Return the thermophysical model
-            virtual const rhoFluidThermo& thermo() const = 0;
+            virtual const rhoThermo& thermo() const = 0;
 
             //- Access the thermophysical model
-            virtual rhoFluidThermo& thermo() = 0;
+            virtual rhoThermo& thermo() = 0;
+
+            //- Return the thermophysical model
+            virtual const rhoFluidThermo& fluidThermo() const = 0;
+
+            //- Access the thermophysical model
+            virtual rhoFluidThermo& fluidThermo() = 0;
 
             //- Return the density field
             virtual const volScalarField& rho() const = 0;

applications/modules/multiphaseEuler/phaseSystem/phaseModels/AnisothermalPhaseModel/AnisothermalPhaseModel.C

@@ -72,7 +72,15 @@ Foam::AnisothermalPhaseModel<BasePhaseModel>::AnisothermalPhaseModel
 )
 :
     BasePhaseModel(fluid, phaseName, referencePhase, index),
-    g_(fluid.mesh().lookupObject<uniformDimensionedVectorField>("g"))
+    g_(fluid.mesh().lookupObject<uniformDimensionedVectorField>("g")),
+    thermophysicalTransport_
+    (
+        PhaseThermophysicalTransportModel
+        <
+            phaseCompressible::momentumTransportModel,
+            transportThermoModel
+        >::New(this->momentumTransport_, this->thermo_)
+    )
 {}
 
 
@@ -101,6 +109,48 @@ bool Foam::AnisothermalPhaseModel<BasePhaseModel>::isothermal() const
 }
 
 
+template<class BasePhaseModel>
+void Foam::AnisothermalPhaseModel<BasePhaseModel>::
+predictThermophysicalTransport()
+{
+    BasePhaseModel::predictThermophysicalTransport();
+    thermophysicalTransport_->predict();
+}
+
+
+template<class BasePhaseModel>
+void Foam::AnisothermalPhaseModel<BasePhaseModel>::
+correctThermophysicalTransport()
+{
+    BasePhaseModel::correctThermophysicalTransport();
+    thermophysicalTransport_->correct();
+}
+
+
+template<class BasePhaseModel>
+Foam::tmp<Foam::scalarField>
+Foam::AnisothermalPhaseModel<BasePhaseModel>::kappaEff(const label patchi) const
+{
+    return thermophysicalTransport_->kappaEff(patchi);
+}
+
+
+template<class BasePhaseModel>
+Foam::tmp<Foam::fvScalarMatrix>
+Foam::AnisothermalPhaseModel<BasePhaseModel>::divq(volScalarField& he) const
+{
+    return thermophysicalTransport_->divq(he);
+}
+
+
+template<class BasePhaseModel>
+Foam::tmp<Foam::fvScalarMatrix>
+Foam::AnisothermalPhaseModel<BasePhaseModel>::divj(volScalarField& Yi) const
+{
+    return thermophysicalTransport_->divj(Yi);
+}
+
+
 template<class BasePhaseModel>
 Foam::tmp<Foam::fvScalarMatrix>
 Foam::AnisothermalPhaseModel<BasePhaseModel>::heEqn()
@@ -145,9 +195,9 @@ Foam::AnisothermalPhaseModel<BasePhaseModel>::heEqn()
             fvc::div
             (
                 fvc::absolute(alphaRhoPhi, alpha, rho, U),
-                this->thermo().p()/rho
+                this->fluidThermo().p()/rho
             )
-          + (fvc::ddt(alpha) - contErr/rho)*this->thermo().p()
+          + (fvc::ddt(alpha) - contErr/rho)*this->fluidThermo().p()
         );
     }
     else if (this->thermo_->dpdt())

applications/modules/multiphaseEuler/phaseSystem/phaseModels/AnisothermalPhaseModel/AnisothermalPhaseModel.H

@@ -38,12 +38,34 @@ SourceFiles
 
 #include "phaseModel.H"
 #include "uniformDimensionedFields.H"
+#include "phaseCompressibleMomentumTransportModel.H"
+#include "PhaseThermophysicalTransportModel.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 namespace Foam
 {
 
+// Trait for converting the ThermoModel's thermo type to the thermo type needed
+// for the thermophysical transport model type; i.e., from rho-type thermo to
+// fluid-type thermo.
+
+template<class ThermoModel>
+struct MovingPhaseModelTransportThermoModel;
+
+template<>
+struct MovingPhaseModelTransportThermoModel<rhoFluidThermo>
+{
+    typedef fluidThermo type;
+};
+
+template<>
+struct MovingPhaseModelTransportThermoModel<rhoFluidMulticomponentThermo>
+{
+    typedef fluidMulticomponentThermo type;
+};
+
+
 /*---------------------------------------------------------------------------*\
                    Class AnisothermalPhaseModel Declaration
 \*---------------------------------------------------------------------------*/
@@ -58,6 +80,24 @@ class AnisothermalPhaseModel
         //- Gravitational acceleration
         const uniformDimensionedVectorField& g_;
 
+        //- Thermo type for the thermophysical transport model
+        typedef
+            typename MovingPhaseModelTransportThermoModel
+            <
+                typename BasePhaseModel::thermoModel
+            >::type
+            transportThermoModel;
+
+        //- Thermophysical transport model
+        autoPtr
+        <
+            PhaseThermophysicalTransportModel
+            <
+                phaseCompressible::momentumTransportModel,
+                transportThermoModel
+            >
+        > thermophysicalTransport_;
+
 
     // Private Member Functions
 
@@ -87,14 +127,31 @@ public:
 
     // Member Functions
 
-        //- Correct the thermodynamics
-        virtual void correctThermo();
-
         //- Return whether the phase is isothermal
         virtual bool isothermal() const;
 
+        //- Correct the thermodynamics
+        virtual void correctThermo();
+
+        //- Predict the energy transport e.g. alphat
+        virtual void predictThermophysicalTransport();
+
+        //- Correct the energy transport e.g. alphat
+        virtual void correctThermophysicalTransport();
+
+        //- Return the effective thermal conductivity on a patch
+        virtual tmp<scalarField> kappaEff(const label patchi) const;
+
+        //- Return the source term for the energy equation
+        virtual tmp<fvScalarMatrix> divq(volScalarField& he) const;
+
+        //- Return the source term for the given specie mass-fraction
+        //  equation
+        virtual tmp<fvScalarMatrix> divj(volScalarField& Yi) const;
+
         //- Return the enthalpy equation
         virtual tmp<fvScalarMatrix> heEqn();
+
 };
 
 

applications/modules/multiphaseEuler/phaseSystem/phaseModels/IsothermalPhaseModel/IsothermalPhaseModel.C

@@ -56,17 +56,17 @@ void Foam::IsothermalPhaseModel<BasePhaseModel>::correctThermo()
     BasePhaseModel::correctThermo();
 
     // Correct the thermo, but make sure that the temperature remains the same
-    tmp<volScalarField> TCopy
+    // tmp<volScalarField> TCopy
-    (
+    // (
-        volScalarField::New
+    //     volScalarField::New
-        (
+    //     (
-            this->thermo().T().name() + ":Copy",
+    //         this->thermo().T().name() + ":Copy",
-            this->thermo().T()
+    //         this->thermo().T()
-        )
+    //     )
-    );
+    // );
-    this->thermo_->he() = this->thermo().he(this->thermo().p(), TCopy);
+    // this->thermo_->he() = this->thermo().he(this->fluidThermo().p(), TCopy);
-    this->thermo_->correct();
+    // this->thermo_->correct();
-    this->thermo_->T() = TCopy;
+    // this->thermo_->T() = TCopy;
 }
 
 
@@ -77,6 +77,15 @@ bool Foam::IsothermalPhaseModel<BasePhaseModel>::isothermal() const
 }
 
 
+template<class BasePhaseModel>
+Foam::tmp<Foam::scalarField>
+Foam::IsothermalPhaseModel<BasePhaseModel>::kappaEff(const label patchi) const
+{
+    NotImplemented;
+    return this->thermo().kappa().boundaryField()[patchi];
+}
+
+
 template<class BasePhaseModel>
 Foam::tmp<Foam::fvScalarMatrix>
 Foam::IsothermalPhaseModel<BasePhaseModel>::heEqn()

applications/modules/multiphaseEuler/phaseSystem/phaseModels/IsothermalPhaseModel/IsothermalPhaseModel.H

@@ -78,6 +78,9 @@ public:
         //- Return whether the phase is isothermal
         virtual bool isothermal() const;
 
+        //- Return the effective thermal conductivity on a patch
+        virtual tmp<scalarField> kappaEff(const label patchi) const;
+
         //- Return the enthalpy equation
         virtual tmp<fvScalarMatrix> heEqn();
 };

applications/modules/multiphaseEuler/phaseSystem/phaseModels/MovingPhaseModel/MovingPhaseModel.C

@@ -175,14 +175,6 @@ Foam::MovingPhaseModel<BasePhaseModel>::MovingPhaseModel
             *this
         )
     ),
-    thermophysicalTransport_
-    (
-        PhaseThermophysicalTransportModel
-        <
-            phaseCompressible::momentumTransportModel,
-            transportThermoModel
-        >::New(momentumTransport_, this->thermo_)
-    ),
     continuityError_
     (
         IOobject
@@ -265,14 +257,6 @@ void Foam::MovingPhaseModel<BasePhaseModel>::predictMomentumTransport()
 }
 
 
-template<class BasePhaseModel>
-void Foam::MovingPhaseModel<BasePhaseModel>::predictThermophysicalTransport()
-{
-    BasePhaseModel::predictThermophysicalTransport();
-    thermophysicalTransport_->predict();
-}
-
-
 template<class BasePhaseModel>
 void Foam::MovingPhaseModel<BasePhaseModel>::correctMomentumTransport()
 {
@@ -281,14 +265,6 @@ void Foam::MovingPhaseModel<BasePhaseModel>::correctMomentumTransport()
 }
 
 
-template<class BasePhaseModel>
-void Foam::MovingPhaseModel<BasePhaseModel>::correctThermophysicalTransport()
-{
-    BasePhaseModel::correctThermophysicalTransport();
-    thermophysicalTransport_->correct();
-}
-
-
 template<class BasePhaseModel>
 void Foam::MovingPhaseModel<BasePhaseModel>::correctUf()
 {
@@ -581,28 +557,4 @@ Foam::MovingPhaseModel<BasePhaseModel>::pPrimef() const
 }
 
 
-template<class BasePhaseModel>
-Foam::tmp<Foam::scalarField>
-Foam::MovingPhaseModel<BasePhaseModel>::kappaEff(const label patchi) const
-{
-    return thermophysicalTransport_->kappaEff(patchi);
-}
-
-
-template<class BasePhaseModel>
-Foam::tmp<Foam::fvScalarMatrix>
-Foam::MovingPhaseModel<BasePhaseModel>::divq(volScalarField& he) const
-{
-    return thermophysicalTransport_->divq(he);
-}
-
-
-template<class BasePhaseModel>
-Foam::tmp<Foam::fvScalarMatrix>
-Foam::MovingPhaseModel<BasePhaseModel>::divj(volScalarField& Yi) const
-{
-    return thermophysicalTransport_->divj(Yi);
-}
-
-
 // ************************************************************************* //

applications/modules/multiphaseEuler/phaseSystem/phaseModels/MovingPhaseModel/MovingPhaseModel.H

@@ -46,7 +46,6 @@ SourceFiles
 #define MovingPhaseModel_H
 
 #include "phaseModel.H"
-#include "PhaseThermophysicalTransportModel.H"
 #include "phaseCompressibleMomentumTransportModel.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@@ -54,25 +53,6 @@ SourceFiles
 namespace Foam
 {
 
-// Trait for converting the ThermoModel's thermo type to the thermo type needed
-// for the thermophysical transport model type; i.e., from rho-type thermo to
-// fluid-type thermo.
-
-template<class ThermoModel>
-struct MovingPhaseModelTransportThermoModel;
-
-template<>
-struct MovingPhaseModelTransportThermoModel<rhoFluidThermo>
-{
-    typedef fluidThermo type;
-};
-
-template<>
-struct MovingPhaseModelTransportThermoModel<rhoFluidMulticomponentThermo>
-{
-    typedef fluidMulticomponentThermo type;
-};
-
 /*---------------------------------------------------------------------------*\
                       Class MovingPhaseModel Declaration
 \*---------------------------------------------------------------------------*/
@@ -84,17 +64,6 @@ class MovingPhaseModel
 {
 protected:
 
-    // Protected typedefs
-
-        //- Thermo type for the thermophysical transport model
-        typedef
-            typename MovingPhaseModelTransportThermoModel
-            <
-                typename BasePhaseModel::thermoModel
-            >::type
-            transportThermoModel;
-
-
     // Protected data
 
         //- Velocity field
@@ -118,16 +87,6 @@ protected:
         //- Turbulence model
         autoPtr<phaseCompressible::momentumTransportModel> momentumTransport_;
 
-        //- Thermophysical transport model
-        autoPtr
-        <
-            PhaseThermophysicalTransportModel
-            <
-                phaseCompressible::momentumTransportModel,
-                transportThermoModel
-            >
-        > thermophysicalTransport_;
-
         //- Continuity error
         volScalarField continuityError_;
 
@@ -175,15 +134,9 @@ public:
         //- Predict the momentumTransport
         virtual void predictMomentumTransport();
 
-        //- Predict the energy transport e.g. alphat
-        virtual void predictThermophysicalTransport();
-
         //- Correct the momentumTransport
         virtual void correctMomentumTransport();
 
-        //- Correct the energy transport e.g. alphat
-        virtual void correctThermophysicalTransport();
-
         //- Correct the face velocity for moving meshes
         virtual void correctUf();
 
@@ -277,19 +230,6 @@ public:
             //- Return the face-phase-pressure'
             // (derivative of phase-pressure w.r.t. phase-fraction)
             virtual tmp<surfaceScalarField> pPrimef() const;
-
-
-        // Thermophysical transport
-
-            //- Return the effective thermal conductivity on a patch
-            virtual tmp<scalarField> kappaEff(const label patchi) const;
-
-            //- Return the source term for the energy equation
-            virtual tmp<fvScalarMatrix> divq(volScalarField& he) const;
-
-            //- Return the source term for the given specie mass-fraction
-            //  equation
-            virtual tmp<fvScalarMatrix> divj(volScalarField& Yi) const;
 };
 
 

applications/modules/multiphaseEuler/phaseSystem/phaseModels/SolidThermalPhaseModel/SolidThermalPhaseModel.C

@@ -26,6 +26,7 @@ License
 #include "SolidThermalPhaseModel.H"
 #include "fvmDdt.H"
 #include "fvmSup.H"
+#include "fvcLaplacian.H"
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
@@ -38,7 +39,11 @@ Foam::SolidThermalPhaseModel<BasePhaseModel>::SolidThermalPhaseModel
     const label index
 )
 :
-    BasePhaseModel(fluid, phaseName, referencePhase, index)
+    BasePhaseModel(fluid, phaseName, referencePhase, index),
+    thermophysicalTransport_
+    (
+        phaseSolidThermophysicalTransportModel::New(*this, this->thermo_)
+    )
 {}
 
 
@@ -51,6 +56,13 @@ Foam::SolidThermalPhaseModel<BasePhaseModel>::~SolidThermalPhaseModel()
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
+template<class BasePhaseModel>
+bool Foam::SolidThermalPhaseModel<BasePhaseModel>::isothermal() const
+{
+    return false;
+}
+
+
 template<class BasePhaseModel>
 void Foam::SolidThermalPhaseModel<BasePhaseModel>::correctThermo()
 {
@@ -61,9 +73,36 @@ void Foam::SolidThermalPhaseModel<BasePhaseModel>::correctThermo()
 
 
 template<class BasePhaseModel>
-bool Foam::SolidThermalPhaseModel<BasePhaseModel>::isothermal() const
+void Foam::SolidThermalPhaseModel<BasePhaseModel>::
+predictThermophysicalTransport()
 {
-    return false;
+    BasePhaseModel::predictThermophysicalTransport();
+    thermophysicalTransport_->predict();
+}
+
+
+template<class BasePhaseModel>
+void Foam::SolidThermalPhaseModel<BasePhaseModel>::
+correctThermophysicalTransport()
+{
+    BasePhaseModel::correctThermophysicalTransport();
+    thermophysicalTransport_->correct();
+}
+
+
+template<class BasePhaseModel>
+Foam::tmp<Foam::scalarField>
+Foam::SolidThermalPhaseModel<BasePhaseModel>::kappaEff(const label patchi) const
+{
+    return thermophysicalTransport_->kappaEff(patchi);
+}
+
+
+template<class BasePhaseModel>
+Foam::tmp<Foam::fvScalarMatrix>
+Foam::SolidThermalPhaseModel<BasePhaseModel>::divq(volScalarField& he) const
+{
+    return thermophysicalTransport_->divq(he);
 }
 
 

applications/modules/multiphaseEuler/phaseSystem/phaseModels/SolidThermalPhaseModel/SolidThermalPhaseModel.H

@@ -38,6 +38,7 @@ SourceFiles
 #define SolidThermalPhaseModel_H
 
 #include "phaseModel.H"
+#include "phaseSolidThermophysicalTransportModel.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -53,6 +54,12 @@ class SolidThermalPhaseModel
 :
     public BasePhaseModel
 {
+    // Thermophysical transport
+
+        //- Pointer to the solid thermophysical transport model
+        autoPtr<phaseSolidThermophysicalTransportModel>
+            thermophysicalTransport_;
+
 
 public:
 
@@ -73,11 +80,23 @@ public:
 
     // Member Functions
 
+        //- Return whether the phase is isothermal
+        virtual bool isothermal() const;
+
         //- Correct the thermodynamics
         virtual void correctThermo();
 
-        //- Return whether the phase is isothermal
+        //- Predict the energy transport e.g. alphat
-        virtual bool isothermal() const;
+        virtual void predictThermophysicalTransport();
+
+        //- Correct the energy transport e.g. alphat
+        virtual void correctThermophysicalTransport();
+
+        //- Return the effective thermal conductivity on a patch
+        virtual tmp<scalarField> kappaEff(const label patchi) const;
+
+        //- Return the source term for the energy equation
+        virtual tmp<fvScalarMatrix> divq(volScalarField& he) const;
 
         //- Return the enthalpy equation
         virtual tmp<fvScalarMatrix> heEqn();

applications/modules/multiphaseEuler/phaseSystem/phaseModels/SolidThermoPhaseModel/SolidThermoPhaseModel.C

@@ -0,0 +1,170 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     | Website:  https://openfoam.org
+    \\  /    A nd           | Copyright (C) 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 "SolidThermoPhaseModel.H"
+#include "phaseSystem.H"
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+template<class BasePhaseModel, class ThermoModel>
+Foam::SolidThermoPhaseModel<BasePhaseModel, ThermoModel>::SolidThermoPhaseModel
+(
+    const phaseSystem& fluid,
+    const word& phaseName,
+    const bool referencePhase,
+    const label index
+)
+:
+    BasePhaseModel(fluid, phaseName, referencePhase, index),
+    viscosity(),
+    thermo_(ThermoModel::New(fluid.mesh(), this->name()))
+{
+    thermo_->validate
+    (
+        IOobject::groupName(phaseModel::typeName, this->name()),
+        "h",
+        "e"
+    );
+}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+template<class BasePhaseModel, class ThermoModel>
+Foam::SolidThermoPhaseModel<BasePhaseModel, ThermoModel>::
+~SolidThermoPhaseModel()
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+template<class BasePhaseModel, class ThermoModel>
+bool
+Foam::SolidThermoPhaseModel<BasePhaseModel, ThermoModel>::incompressible() const
+{
+    return thermo_().incompressible();
+}
+
+
+template<class BasePhaseModel, class ThermoModel>
+bool Foam::SolidThermoPhaseModel<BasePhaseModel, ThermoModel>::isochoric() const
+{
+    return thermo_().isochoric();
+}
+
+
+template<class BasePhaseModel, class ThermoModel>
+const Foam::rhoThermo&
+Foam::SolidThermoPhaseModel<BasePhaseModel, ThermoModel>::thermo() const
+{
+    return thermo_();
+}
+
+
+template<class BasePhaseModel, class ThermoModel>
+Foam::rhoThermo&
+Foam::SolidThermoPhaseModel<BasePhaseModel, ThermoModel>::thermo()
+{
+    return thermo_();
+}
+
+
+template<class BasePhaseModel, class ThermoModel>
+const Foam::rhoFluidThermo&
+Foam::SolidThermoPhaseModel<BasePhaseModel, ThermoModel>::fluidThermo() const
+{
+    NotImplemented;
+    return refCast<const rhoFluidThermo>(thermo_());
+}
+
+
+template<class BasePhaseModel, class ThermoModel>
+Foam::rhoFluidThermo&
+Foam::SolidThermoPhaseModel<BasePhaseModel, ThermoModel>::fluidThermo()
+{
+    NotImplemented;
+    return refCast<rhoFluidThermo>(thermo_());
+}
+
+
+template<class BasePhaseModel, class ThermoModel>
+const Foam::volScalarField&
+Foam::SolidThermoPhaseModel<BasePhaseModel, ThermoModel>::rho() const
+{
+    return thermo_->rho();
+}
+
+
+template<class BasePhaseModel, class ThermoModel>
+Foam::volScalarField&
+Foam::SolidThermoPhaseModel<BasePhaseModel, ThermoModel>::rho()
+{
+    return thermo_->rho();
+}
+
+
+template<class BasePhaseModel, class ThermoModel>
+Foam::tmp<Foam::volScalarField>
+Foam::SolidThermoPhaseModel<BasePhaseModel, ThermoModel>::mu() const
+{
+    NotImplemented;
+    return tmp<volScalarField>(nullptr);
+}
+
+
+template<class BasePhaseModel, class ThermoModel>
+Foam::tmp<Foam::scalarField>
+Foam::SolidThermoPhaseModel<BasePhaseModel, ThermoModel>::mu
+(
+    const label patchi
+) const
+{
+    NotImplemented;
+    return tmp<scalarField>(nullptr);
+}
+
+
+template<class BasePhaseModel, class ThermoModel>
+Foam::tmp<Foam::volScalarField>
+Foam::SolidThermoPhaseModel<BasePhaseModel, ThermoModel>::nu() const
+{
+    NotImplemented;
+    return tmp<volScalarField>(nullptr);
+}
+
+
+template<class BasePhaseModel, class ThermoModel>
+Foam::tmp<Foam::scalarField>
+Foam::SolidThermoPhaseModel<BasePhaseModel, ThermoModel>::nu
+(
+    const label patchi
+) const
+{
+    NotImplemented;
+    return tmp<scalarField>(nullptr);
+}
+
+
+// ************************************************************************* //

applications/modules/multiphaseEuler/phaseSystem/phaseModels/SolidThermoPhaseModel/SolidThermoPhaseModel.H

@@ -0,0 +1,149 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     | Website:  https://openfoam.org
+    \\  /    A nd           | Copyright (C) 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/>.
+
+Class
+    Foam::SolidThermoPhaseModel
+
+Description
+    Class which represents a solid phase with a thermodynamic model. Provides
+    access to the thermodynamic variables. Note that the thermo model itself is
+    not returned as this class could be substituted in the hierarchy for one
+    which mirrors the functionality, but does not include a thermo model; an
+    incompressible phase model, for example.
+
+SourceFiles
+    SolidThermoPhaseModel.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef SolidThermoPhaseModel_H
+#define SolidThermoPhaseModel_H
+
+#include "phaseModel.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+class rhoThermo;
+class rhoFluidThermo;
+
+/*---------------------------------------------------------------------------*\
+                      Class SolidThermoPhaseModel Declaration
+\*---------------------------------------------------------------------------*/
+
+template<class BasePhaseModel, class ThermoModel>
+class SolidThermoPhaseModel
+:
+    public BasePhaseModel,
+    public viscosity
+{
+protected:
+
+    // Protected data
+
+        //- Thermophysical model
+        autoPtr<ThermoModel> thermo_;
+
+
+public:
+
+    typedef ThermoModel thermoModel;
+
+
+    // Constructors
+
+        SolidThermoPhaseModel
+        (
+            const phaseSystem& fluid,
+            const word& phaseName,
+            const bool referencePhase,
+            const label index
+        );
+
+
+    //- Destructor
+    virtual ~SolidThermoPhaseModel();
+
+
+    // Member Functions
+
+        // Thermo
+
+            //- Return whether the phase is incompressible
+            virtual bool incompressible() const;
+
+            //- Return whether the phase is constant density
+            virtual bool isochoric() const;
+
+            //- Return the thermophysical model
+            virtual const rhoThermo& thermo() const;
+
+            //- Access the thermophysical model
+            virtual rhoThermo& thermo();
+
+            //- Return the thermophysical model
+            virtual const rhoFluidThermo& fluidThermo() const;
+
+            //- Access the thermophysical model
+            virtual rhoFluidThermo& fluidThermo();
+
+            //- Return the density field
+            virtual const volScalarField& rho() const;
+
+            //- Access the the density field
+            virtual volScalarField& rho();
+
+
+        // Transport
+
+            //- Return the laminar dynamic viscosity
+            virtual tmp<volScalarField> mu() const;
+
+            //- Return the laminar dynamic viscosity on a patch
+            virtual tmp<scalarField> mu(const label patchi) const;
+
+            //- Return the laminar kinematic viscosity
+            virtual tmp<volScalarField> nu() const;
+
+            //- Return the laminar kinematic viscosity on a patch
+            virtual tmp<scalarField> nu(const label patchi) const;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#ifdef NoRepository
+    #include "SolidThermoPhaseModel.C"
+#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

applications/modules/multiphaseEuler/phaseSystem/phaseModels/StationaryPhaseModel/StationaryPhaseModel.C

@@ -24,7 +24,6 @@ License
 \*---------------------------------------------------------------------------*/
 
 #include "StationaryPhaseModel.H"
-#include "fvcLaplacian.H"
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
@@ -335,14 +334,6 @@ void Foam::StationaryPhaseModel<BasePhaseModel>::divU
 }
 
 
-template<class BasePhaseModel>
-Foam::tmp<Foam::scalarField>
-Foam::StationaryPhaseModel<BasePhaseModel>::kappaEff(const label patchi) const
-{
-    return this->thermo().kappa().boundaryField()[patchi];
-}
-
-
 template<class BasePhaseModel>
 Foam::tmp<Foam::volScalarField>
 Foam::StationaryPhaseModel<BasePhaseModel>::k() const
@@ -368,28 +359,4 @@ Foam::StationaryPhaseModel<BasePhaseModel>::pPrimef() const
 }
 
 
-template<class BasePhaseModel>
-Foam::tmp<Foam::fvScalarMatrix>
-Foam::StationaryPhaseModel<BasePhaseModel>::divq(volScalarField& he) const
-{
-    const volScalarField& alpha = *this;
-
-    const surfaceScalarField alphaKappa
-    (
-        alpha.name() + '*' + this->thermo().kappa().name(),
-        fvc::interpolate(alpha)*fvc::interpolate(this->thermo().kappa())
-    );
-
-    // Return heat flux source as an implicit energy correction
-    // to the temperature gradient flux
-    return
-       -fvc::laplacian(alphaKappa, this->thermo().T())
-       -fvm::laplacianCorrection
-        (
-            alphaKappa/fvc::interpolate(this->thermo().Cpv()),
-            he
-        );
-}
-
-
 // ************************************************************************* //

applications/modules/multiphaseEuler/phaseSystem/phaseModels/StationaryPhaseModel/StationaryPhaseModel.H

@@ -160,21 +160,12 @@ public:
 
         // Momentum transport
 
-            //- Return the effective thermal conductivity on a patch
-            virtual tmp<scalarField> kappaEff(const label patchi) const;
-
             //- Return the turbulent kinetic energy
             virtual tmp<volScalarField> k() const;
 
             //- Return the face-phase-pressure'
             // (derivative of phase-pressure w.r.t. phase-fraction)
             virtual tmp<surfaceScalarField> pPrimef() const;
-
-
-        // Thermophysical transport
-
-            //- Return the source term for the energy equation
-            virtual tmp<fvScalarMatrix> divq(volScalarField& he) const;
 };
 
 

applications/modules/multiphaseEuler/phaseSystem/phaseModels/ThermoPhaseModel/ThermoPhaseModel.C

@@ -74,7 +74,7 @@ bool Foam::ThermoPhaseModel<BasePhaseModel, ThermoModel>::isochoric() const
 
 
 template<class BasePhaseModel, class ThermoModel>
-const Foam::rhoFluidThermo&
+const Foam::rhoThermo&
 Foam::ThermoPhaseModel<BasePhaseModel, ThermoModel>::thermo() const
 {
     return thermo_();
@@ -82,13 +82,29 @@ Foam::ThermoPhaseModel<BasePhaseModel, ThermoModel>::thermo() const
 
 
 template<class BasePhaseModel, class ThermoModel>
-Foam::rhoFluidThermo&
+Foam::rhoThermo&
 Foam::ThermoPhaseModel<BasePhaseModel, ThermoModel>::thermo()
 {
     return thermo_();
 }
 
 
+template<class BasePhaseModel, class ThermoModel>
+const Foam::rhoFluidThermo&
+Foam::ThermoPhaseModel<BasePhaseModel, ThermoModel>::fluidThermo() const
+{
+    return thermo_();
+}
+
+
+template<class BasePhaseModel, class ThermoModel>
+Foam::rhoFluidThermo&
+Foam::ThermoPhaseModel<BasePhaseModel, ThermoModel>::fluidThermo()
+{
+    return thermo_();
+}
+
+
 template<class BasePhaseModel, class ThermoModel>
 const Foam::volScalarField&
 Foam::ThermoPhaseModel<BasePhaseModel, ThermoModel>::rho() const

applications/modules/multiphaseEuler/phaseSystem/phaseModels/ThermoPhaseModel/ThermoPhaseModel.H

@@ -46,6 +46,7 @@ SourceFiles
 namespace Foam
 {
 
+class rhoThermo;
 class rhoFluidThermo;
 
 /*---------------------------------------------------------------------------*\
@@ -97,10 +98,16 @@ public:
             virtual bool isochoric() const;
 
             //- Return the thermophysical model
-            virtual const rhoFluidThermo& thermo() const;
+            virtual const rhoThermo& thermo() const;
 
             //- Access the thermophysical model
-            virtual rhoFluidThermo& thermo();
+            virtual rhoThermo& thermo();
+
+            //- Return the thermophysical model
+            virtual const rhoFluidThermo& fluidThermo() const;
+
+            //- Access the thermophysical model
+            virtual rhoFluidThermo& fluidThermo();
 
             //- Return the density field
             virtual const volScalarField& rho() const;

applications/modules/multiphaseEuler/phaseSystem/phaseModels/phaseModels.C

@@ -28,10 +28,13 @@ License
 #include "rhoFluidThermo.H"
 #include "rhoFluidMulticomponentThermo.H"
 
+#include "solidThermo.H"
+
 #include "combustionModel.H"
 
 #include "phaseModel.H"
 #include "ThermoPhaseModel.H"
+#include "SolidThermoPhaseModel.H"
 #include "IsothermalPhaseModel.H"
 #include "AnisothermalPhaseModel.H"
 #include "SolidThermalPhaseModel.H"
@@ -83,10 +86,10 @@ namespace Foam
                 <
                     StationaryPhaseModel
                     <
-                        ThermoPhaseModel
+                        SolidThermoPhaseModel
                         <
                             phaseModel,
-                            rhoFluidThermo
+                            solidThermo
                         >
                     >
                 >
@@ -139,10 +142,10 @@ namespace Foam
                 <
                     StationaryPhaseModel
                     <
-                        ThermoPhaseModel
+                        SolidThermoPhaseModel
                         <
                             phaseModel,
-                            rhoFluidThermo
+                            solidThermo
                         >
                     >
                 >
@@ -159,9 +162,9 @@ namespace Foam
     );
 
     typedef
-        AnisothermalPhaseModel
+        MulticomponentPhaseModel
         <
-            MulticomponentPhaseModel
+            AnisothermalPhaseModel
             <
                 InertPhaseModel
                 <
@@ -197,37 +200,9 @@ namespace Foam
     );
 
     typedef
-        IsothermalPhaseModel
+        MulticomponentPhaseModel
         <
-            MulticomponentPhaseModel
+            AnisothermalPhaseModel
-            <
-                InertPhaseModel
-                <
-                    MovingPhaseModel
-                    <
-                        ThermoPhaseModel
-                        <
-                            phaseModel,
-                            rhoFluidMulticomponentThermo
-                        >
-                    >
-                >
-            >
-        >
-        multicomponentIsothermalPhaseModel;
-
-    addNamedToRunTimeSelectionTable
-    (
-        phaseModel,
-        multicomponentIsothermalPhaseModel,
-        phaseSystem,
-        multicomponentIsothermalPhaseModel
-    );
-
-    typedef
-        AnisothermalPhaseModel
-        <
-            MulticomponentPhaseModel
             <
                 ReactingPhaseModel
                 <

applications/modules/multiphaseEuler/phaseSystem/phaseSystem/phaseSystem.C

@@ -606,7 +606,7 @@ void Foam::phaseSystem::correctKinematics()
     // Update the pressure time-derivative if required
     if (updateDpdt)
     {
-        dpdt_ = fvc::ddt(phaseModels_.begin()().thermo().p());
+        dpdt_ = fvc::ddt(phaseModels_.begin()().fluidThermo().p());
     }
 }
 
@@ -800,7 +800,7 @@ void Foam::phaseSystem::correctPhi
             phaseModel& phase = phases()[phasei];
             const volScalarField& alpha = phase;
 
-            psi += alpha*phase.thermo().psi()/phase.rho();
+            psi += alpha*phase.fluidThermo().psi()/phase.rho();
         }
 
         fv::correctPhi

applications/modules/multiphaseEuler/phaseSystems/HeatTransferPhaseSystem/HeatTransferPhaseSystem.C

@@ -47,8 +47,8 @@ void Foam::HeatTransferPhaseSystem<BasePhaseSystem>::addDmdtHefs
 
         const phaseModel& phase1 = interface.phase1();
         const phaseModel& phase2 = interface.phase2();
-        const rhoFluidThermo& thermo1 = phase1.thermo();
+        const rhoFluidThermo& thermo1 = phase1.fluidThermo();
-        const rhoFluidThermo& thermo2 = phase2.thermo();
+        const rhoFluidThermo& thermo2 = phase2.fluidThermo();
         const volScalarField& he1 = thermo1.he();
         const volScalarField& he2 = thermo2.he();
         const volScalarField hs1(thermo1.hs());
@@ -89,8 +89,8 @@ void Foam::HeatTransferPhaseSystem<BasePhaseSystem>::addDmidtHefs
 
         const phaseModel& phase1 = interface.phase1();
         const phaseModel& phase2 = interface.phase2();
-        const rhoFluidThermo& thermo1 = phase1.thermo();
+        const rhoFluidThermo& thermo1 = phase1.fluidThermo();
-        const rhoFluidThermo& thermo2 = phase2.thermo();
+        const rhoFluidThermo& thermo2 = phase2.fluidThermo();
         const rhoFluidMulticomponentThermo* mcThermoPtr1 =
             isA<rhoFluidMulticomponentThermo>(thermo1)
           ? &refCast<const rhoFluidMulticomponentThermo>(thermo1)
@@ -195,8 +195,8 @@ void Foam::HeatTransferPhaseSystem<BasePhaseSystem>::addDmdtHefsWithoutL
 
         const phaseModel& phase1 = interface.phase1();
         const phaseModel& phase2 = interface.phase2();
-        const rhoFluidThermo& thermo1 = phase1.thermo();
+        const rhoFluidThermo& thermo1 = phase1.fluidThermo();
-        const rhoFluidThermo& thermo2 = phase2.thermo();
+        const rhoFluidThermo& thermo2 = phase2.fluidThermo();
         const volScalarField& he1 = thermo1.he();
         const volScalarField& he2 = thermo2.he();
         const volScalarField K1(phase1.K());
@@ -305,8 +305,8 @@ void Foam::HeatTransferPhaseSystem<BasePhaseSystem>::addDmidtHefsWithoutL
 
         const phaseModel& phase1 = interface.phase1();
         const phaseModel& phase2 = interface.phase2();
-        const rhoFluidThermo& thermo1 = phase1.thermo();
+        const rhoFluidThermo& thermo1 = phase1.fluidThermo();
-        const rhoFluidThermo& thermo2 = phase2.thermo();
+        const rhoFluidThermo& thermo2 = phase2.fluidThermo();
         const rhoFluidMulticomponentThermo* mcThermoPtr1 =
             isA<rhoFluidMulticomponentThermo>(thermo1)
           ? &refCast<const rhoFluidMulticomponentThermo>(thermo1)
@@ -512,8 +512,8 @@ Foam::HeatTransferPhaseSystem<BasePhaseSystem>::L
     const latentHeatScheme scheme
 ) const
 {
-    const rhoFluidThermo& thermo1 = interface.phase1().thermo();
+    const rhoFluidThermo& thermo1 = interface.phase1().fluidThermo();
-    const rhoFluidThermo& thermo2 = interface.phase2().thermo();
+    const rhoFluidThermo& thermo2 = interface.phase2().fluidThermo();
 
     // Interface enthalpies
     const volScalarField haf1(thermo1.ha(thermo1.p(), Tf));
@@ -552,8 +552,8 @@ Foam::HeatTransferPhaseSystem<BasePhaseSystem>::L
     const latentHeatScheme scheme
 ) const
 {
-    const rhoFluidThermo& thermo1 = interface.phase1().thermo();
+    const rhoFluidThermo& thermo1 = interface.phase1().fluidThermo();
-    const rhoFluidThermo& thermo2 = interface.phase2().thermo();
+    const rhoFluidThermo& thermo2 = interface.phase2().fluidThermo();
 
     // Interface enthalpies
     const scalarField haf1(thermo1.ha(Tf, cells));
@@ -595,8 +595,8 @@ Foam::HeatTransferPhaseSystem<BasePhaseSystem>::Li
     const latentHeatScheme scheme
 ) const
 {
-    const rhoFluidThermo& thermo1 = interface.phase1().thermo();
+    const rhoFluidThermo& thermo1 = interface.phase1().fluidThermo();
-    const rhoFluidThermo& thermo2 = interface.phase2().thermo();
+    const rhoFluidThermo& thermo2 = interface.phase2().fluidThermo();
     const rhoFluidMulticomponentThermo* mcThermoPtr1 =
         isA<rhoFluidMulticomponentThermo>(thermo1)
       ? &refCast<const rhoFluidMulticomponentThermo>(thermo1)
@@ -668,8 +668,8 @@ Foam::HeatTransferPhaseSystem<BasePhaseSystem>::Li
     const latentHeatScheme scheme
 ) const
 {
-    const rhoFluidThermo& thermo1 = interface.phase1().thermo();
+    const rhoFluidThermo& thermo1 = interface.phase1().fluidThermo();
-    const rhoFluidThermo& thermo2 = interface.phase2().thermo();
+    const rhoFluidThermo& thermo2 = interface.phase2().fluidThermo();
     const rhoFluidMulticomponentThermo* mcThermoPtr1 =
         isA<rhoFluidMulticomponentThermo>(thermo1)
       ? &refCast<const rhoFluidMulticomponentThermo>(thermo1)

applications/modules/multiphaseEuler/phaseSystems/ThermalPhaseChangePhaseSystem/ThermalPhaseChangePhaseSystem.C

@@ -180,8 +180,8 @@ ThermalPhaseChangePhaseSystem
                     IOobject::AUTO_WRITE
                 ),
                 (
-                    interface.phase1().thermo().T()
+                    interface.phase1().fluidThermo().T()
-                  + interface.phase2().thermo().T()
+                  + interface.phase2().fluidThermo().T()
                 )/2
             )
         );
@@ -205,7 +205,7 @@ ThermalPhaseChangePhaseSystem
                 ),
                 saturationModels_[interface]->Tsat
                 (
-                    interface.phase1().thermo().p()
+                    interface.phase1().fluidThermo().p()
                 )
             )
         );
@@ -363,7 +363,7 @@ Foam::ThermalPhaseChangePhaseSystem<BasePhaseSystem>::heatTransfer() const
         Tns.insert
         (
             interface,
-            satModel.Tsat(interface.phase1().thermo().p()).ptr()
+            satModel.Tsat(interface.phase1().fluidThermo().p()).ptr()
         );
     }
 
@@ -457,7 +457,7 @@ Foam::ThermalPhaseChangePhaseSystem<BasePhaseSystem>::heatTransfer() const
                     fvm::Sp
                     (
                         dmdt0s_[phase.index()] - dmdts[phase.index()],
-                        phase.thermo().he()
+                        phase.fluidThermo().he()
                     );
             }
         }
@@ -562,8 +562,8 @@ Foam::ThermalPhaseChangePhaseSystem<BasePhaseSystem>::correctInterfaceThermo()
         const phaseInterface& interface = saturationModelIter()->interface();
         const phaseModel& phase1 = interface.phase1();
         const phaseModel& phase2 = interface.phase2();
-        const rhoFluidThermo& thermo1 = phase1.thermo();
+        const rhoFluidThermo& thermo1 = phase1.fluidThermo();
-        const rhoFluidThermo& thermo2 = phase2.thermo();
+        const rhoFluidThermo& thermo2 = phase2.fluidThermo();
         const volScalarField& T1(thermo1.T());
         const volScalarField& T2(thermo2.T());
 

applications/modules/multiphaseEuler/phaseSystems/TwoResistanceHeatTransferPhaseSystem/TwoResistanceHeatTransferPhaseSystem.C

@@ -59,8 +59,8 @@ void Foam::TwoResistanceHeatTransferPhaseSystem<BasePhaseSystem>::addDmdtHefs
 
         const phaseModel& phase1 = interface.phase1();
         const phaseModel& phase2 = interface.phase2();
-        const rhoFluidThermo& thermo1 = phase1.thermo();
+        const rhoThermo& thermo1 = phase1.thermo();
-        const rhoFluidThermo& thermo2 = phase2.thermo();
+        const rhoThermo& thermo2 = phase2.thermo();
 
         // Transfer coefficients
         const sidedBlendedHeatTransferModel& heatTransferModel =
@@ -124,8 +124,8 @@ void Foam::TwoResistanceHeatTransferPhaseSystem<BasePhaseSystem>::addDmidtHefs
 
         const phaseModel& phase1 = interface.phase1();
         const phaseModel& phase2 = interface.phase2();
-        const rhoFluidThermo& thermo1 = phase1.thermo();
+        const rhoThermo& thermo1 = phase1.thermo();
-        const rhoFluidThermo& thermo2 = phase2.thermo();
+        const rhoThermo& thermo2 = phase2.thermo();
 
         // Transfer coefficients
         const sidedBlendedHeatTransferModel& heatTransferModel =

applications/modules/multiphaseEuler/populationBalance/binaryBreakupModels/Liao/Liao.C

@@ -122,7 +122,7 @@ addToBinaryBreakupRate
     {
         const volScalarField tauShear
         (
-            BShear_*continuousPhase.thermo().mu()*shearStrainRate_
+            BShear_*continuousPhase.fluidThermo().mu()*shearStrainRate_
         );
 
         binaryBreakupRate +=
@@ -135,7 +135,7 @@ addToBinaryBreakupRate
         const volScalarField tauEddy
         (
             pos0(kolmogorovLengthScale_ - fj.dSph())
-           *BEddy_*continuousPhase.thermo().mu()*eddyStrainRate_
+           *BEddy_*continuousPhase.fluidThermo().mu()*eddyStrainRate_
         );
 
         binaryBreakupRate +=

applications/modules/multiphaseEuler/populationBalance/binaryBreakupModels/Liao/LiaoBase.C

@@ -93,7 +93,7 @@ void Foam::diameterModels::LiaoBase::precompute()
     kolmogorovLengthScale_ =
         pow025
         (
-            pow3(populationBalance_.continuousPhase().thermo().nu())
+            pow3(populationBalance_.continuousPhase().fluidThermo().nu())
            /populationBalance_.continuousTurbulence().epsilon()
         );
 
@@ -106,7 +106,7 @@ void Foam::diameterModels::LiaoBase::precompute()
            (
                populationBalance_.continuousPhase().rho()
               *populationBalance_.continuousTurbulence().epsilon()
-              /populationBalance_.continuousPhase().thermo().mu()
+              /populationBalance_.continuousPhase().fluidThermo().mu()
            );
 
     if (uTerminal_.empty())
@@ -119,7 +119,7 @@ void Foam::diameterModels::LiaoBase::precompute()
         (
             "nuc",
             dimViscosity,
-            gAverage(populationBalance_.continuousPhase().thermo().nu()())
+            gAverage(populationBalance_.continuousPhase().fluidThermo().nu()())
         );
 
         const dimensionedScalar rhoc

applications/modules/multiphaseEuler/populationBalance/binaryBreakupModels/LuoSvendsen/LuoSvendsen.C

@@ -156,7 +156,7 @@ void Foam::diameterModels::binaryBreakupModels::LuoSvendsen::precompute()
         (
             pow3
             (
-                popBal_.continuousPhase().thermo().nu()
+                popBal_.continuousPhase().fluidThermo().nu()
             )
            /popBal_.continuousTurbulence().epsilon()
         );

applications/modules/multiphaseEuler/populationBalance/breakupModels/Kusters/Kusters.C

@@ -76,7 +76,7 @@ void Foam::diameterModels::breakupModels::Kusters::setBreakupRate
         sqrt
         (
             4*popBal_.continuousTurbulence().epsilon()
-           /(15*pi*popBal_.continuousPhase().thermo().nu())
+           /(15*pi*popBal_.continuousPhase().fluidThermo().nu())
         )
        *exp
         (

applications/modules/multiphaseEuler/populationBalance/breakupModels/Laakkonen/Laakkonen.C

@@ -94,7 +94,7 @@ Foam::diameterModels::breakupModels::Laakkonen::setBreakupRate
                     continuousPhase.rho()*pow(fi.dSph(), 5.0/3.0)
                    *pow(popBal_.continuousTurbulence().epsilon(), 2.0/3.0)
                 )
-              + C3_*continuousPhase.thermo().mu()
+              + C3_*continuousPhase.fluidThermo().mu()
                /(
                     sqrt(continuousPhase.rho()*fi.phase().rho())
                    *cbrt(popBal_.continuousTurbulence().epsilon())

applications/modules/multiphaseEuler/populationBalance/coalescenceModels/AdachiStuartFokkink/AdachiStuartFokkink.C

@@ -79,7 +79,7 @@ addToCoalescenceRate
        *sqrt
         (
             0.3*pi*popBal_.continuousTurbulence().epsilon()
-           /popBal_.continuousPhase().thermo().nu()
+           /popBal_.continuousPhase().fluidThermo().nu()
         )
        *pow3(fi.d() + fj.d());
 }

applications/modules/multiphaseEuler/populationBalance/coalescenceModels/BrownianCollisions/BrownianCollisions.C

@@ -88,7 +88,7 @@ BrownianCollisions
 void Foam::diameterModels::coalescenceModels::BrownianCollisions::precompute()
 {
     const volScalarField& T = popBal_.continuousPhase().thermo().T();
-    const volScalarField& p = popBal_.continuousPhase().thermo().p();
+    const volScalarField& p = popBal_.continuousPhase().fluidThermo().p();
 
     lambda_ = k*T/(sqrt(2.0)*pi*p*sqr(sigma_));
 }
@@ -108,7 +108,7 @@ addToCoalescenceRate
 
     const volScalarField& T = popBal_.continuousPhase().thermo().T();
 
-    tmp<volScalarField> tmu(popBal_.continuousPhase().thermo().mu());
+    tmp<volScalarField> tmu(popBal_.continuousPhase().fluidThermo().mu());
     const volScalarField& mu = tmu();
 
     const volScalarField Cci

applications/modules/multiphaseEuler/populationBalance/coalescenceModels/CoulaloglouTavlarides/CoulaloglouTavlarides.C

@@ -82,7 +82,7 @@ addToCoalescenceRate
        *cbrt(popBal_.continuousTurbulence().epsilon())/(1 + popBal_.alphas())
        *exp
         (
-          - C2_*continuousPhase.thermo().mu()*continuousPhase.rho()
+          - C2_*continuousPhase.fluidThermo().mu()*continuousPhase.rho()
            *popBal_.continuousTurbulence().epsilon()
            /sqr(popBal_.sigmaWithContinuousPhase(fi.phase()))
            /pow3(1 + popBal_.alphas())

applications/modules/multiphaseEuler/populationBalance/coalescenceModels/LiaoCoalescence/LiaoCoalescence.C

@@ -231,7 +231,7 @@ addToCoalescenceRate
       + pos0(kolmogorovLengthScale_ - (fi.dSph() + fj.dSph()))
        *exp
         (
-          - 3*continuousPhase.thermo().mu()*dEq*eddyStrainRate_
+          - 3*continuousPhase.fluidThermo().mu()*dEq*eddyStrainRate_
            /(4*popBal_.sigmaWithContinuousPhase(fi.phase()))
            *log
             (

applications/modules/multiphaseEuler/populationBalance/coalescenceModels/turbulentShear/turbulentShear.C

@@ -76,7 +76,7 @@ addToCoalescenceRate
 
     const volScalarField& rho = popBal_.continuousPhase().rho();
     tmp<volScalarField> epsilon(popBal_.continuousTurbulence().epsilon());
-    tmp<volScalarField> mu(popBal_.continuousPhase().thermo().mu());
+    tmp<volScalarField> mu(popBal_.continuousPhase().fluidThermo().mu());
 
     coalescenceRate += C_*sqrt(epsilon*rho/mu)*pow3(fi.d() + fj.d());
 }

applications/modules/multiphaseEuler/thermophysicalTransportModels/derivedFvPatchFields/alphatWallBoilingWallFunction/alphatWallBoilingWallFunctionFvPatchScalarField.C

@@ -225,7 +225,7 @@ struct alphatWallBoilingWallFunctionFvPatchScalarField::boilingLiquidProperties
         trhoVapourw(vapour.thermo().rho(patchi())),
         rhoVapourw(trhoVapourw()),
         Cpw(liquid.thermo().Cp().boundaryField()[patchi()]),
-        tnuw(liquid.thermo().nu(patchi())),
+        tnuw(liquid.fluidThermo().nu(patchi())),
         nuw(tnuw()),
         kappaByCp
         (
@@ -265,7 +265,7 @@ struct alphatWallBoilingWallFunctionFvPatchScalarField::boilingLiquidProperties
                 interfaceSaturationTemperatureModel
             >
             (interface)
-           .Tsat(liquid.thermo().p())()
+           .Tsat(liquid.fluidThermo().p())()
            .boundaryField()[patchi()]
         ),
         L

applications/modules/multiphaseEuler/thermophysicalTransportModels/derivedFvPatchFields/coupledMultiphaseTemperature/coupledMultiphaseTemperatureFvPatchScalarField.C

@@ -52,7 +52,7 @@ void Foam::coupledMultiphaseTemperatureFvPatchScalarField::getThis
     forAll(fluid.phases(), phasei)
     {
         const phaseModel& phase = fluid.phases()[phasei];
-        const fluidThermo& thermo = phase.thermo();
+        const rhoThermo& thermo = phase.thermo();
 
         const fvPatchScalarField& Tw =
             thermo.T().boundaryField()[patch().index()];
@@ -105,7 +105,7 @@ void Foam::coupledMultiphaseTemperatureFvPatchScalarField::getNbr
     forAll(fluid.phases(), phasei)
     {
         const phaseModel& phase = fluid.phases()[phasei];
-        const fluidThermo& thermo = phase.thermo();
+        const rhoThermo& thermo = phase.thermo();
 
         const fvPatchScalarField& alpha =
             phase.boundaryField()[patch().index()];
@@ -143,7 +143,7 @@ void Foam::coupledMultiphaseTemperatureFvPatchScalarField::getNbr
     forAll(fluid.phases(), phasei)
     {
         const phaseModel& phase = fluid.phases()[phasei];
-        const fluidThermo& thermo = phase.thermo();
+        const rhoThermo& thermo = phase.thermo();
 
         const fvPatchScalarField& alpha =
             phase.boundaryField()[patch().index()];

applications/modules/multiphaseEuler/thermophysicalTransportModels/heatTransferModels/wallBoiling/wallBoilingHeatTransfer.C

@@ -237,9 +237,9 @@ Foam::heatTransferModels::wallBoilingHeatTransfer::K
     const phaseModel& vapour = fluid.phases()[vapourPhaseName_];
     const phaseModel& solid = interface_.dispersed();
 
-    const rhoFluidThermo& lThermo = liquid.thermo();
+    const rhoThermo& lThermo = liquid.thermo();
-    const rhoFluidThermo& vThermo = vapour.thermo();
+    const rhoThermo& vThermo = vapour.thermo();
-    const rhoFluidThermo& sThermo = solid.thermo();
+    const rhoThermo& sThermo = solid.thermo();
 
     // Estimate the surface temperature from the surrounding temperature and
     // heat transfer coefficients. Note that a lagged value of K is used in
@@ -263,12 +263,12 @@ Foam::heatTransferModels::wallBoilingHeatTransfer::K
             interfaceSaturationTemperatureModel
         >
         (phaseInterface(liquid, vapour))
-       .Tsat(liquid.thermo().p())()
+       .Tsat(liquid.fluidThermo().p())()
     );
 
     const volScalarField L
     (
-        vThermo.ha(lThermo.p(), Tsat) - lThermo.ha()
+        vThermo.ha(liquid.fluidThermo().p(), Tsat) - lThermo.ha()
     );
 
     // Wetted fraction

src/MomentumTransportModels/phaseCompressible/RAS/LaheyKEpsilon/LaheyKEpsilon.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2013-2022 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2013-2023 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -203,7 +203,7 @@ tmp<volScalarField> LaheyKEpsilon<BasicMomentumTransportModel>::bubbleG() const
         Cp_
        *(
             pow3(magUr)
-          + pow(drag.CdRe()*liquid.thermo().nu()/gas.d(), 4.0/3.0)
+          + pow(drag.CdRe()*liquid.fluidThermo().nu()/gas.d(), 4.0/3.0)
            *pow(magUr, 5.0/3.0)
         )
        *gas

src/MomentumTransportModels/phaseCompressible/RAS/continuousGasKEpsilon/continuousGasKEpsilon.C

@@ -132,7 +132,7 @@ void continuousGasKEpsilon<BasicMomentumTransportModel>::correctNut()
     volScalarField rhodv(gas.rho() + virtualMass.Cvm()*liquid.rho());
     volScalarField thetag
     (
-        (rhodv/(18*liquid.rho()*liquid.thermo().nu()))*sqr(gas.d())
+        (rhodv/(18*liquid.rho()*liquid.fluidThermo().nu()))*sqr(gas.d())
     );
     volScalarField expThetar
     (

src/MomentumTransportModels/phaseCompressible/RAS/mixtureKEpsilon/mixtureKEpsilon.C

@@ -472,7 +472,7 @@ mixtureKEpsilon<BasicMomentumTransportModel>::bubbleG() const
        *pos(alphap_ - gas)*liquid*liquid.rho()
        *(
             pow3(magUr)
-          + pow(drag.CdRe()*liquid.thermo().nu()/gas.d(), 4.0/3.0)
+          + pow(drag.CdRe()*liquid.fluidThermo().nu()/gas.d(), 4.0/3.0)
            *pow(magUr, 5.0/3.0)
         )
        *gas

src/ThermophysicalTransportModels/Allwmake

@@ -10,8 +10,12 @@ wmake $targetType solid
 
 wmake $targetType fluidThermo
 wmake $targetType fluidMulticomponentThermo
+
 wmake $targetType phaseFluidThermo
 wmake $targetType phaseFluidMulticomponentThermo
+
+wmake $targetType phaseSolid
+
 wmake $targetType coupledThermophysicalTransportModels
 
 #------------------------------------------------------------------------------

src/ThermophysicalTransportModels/phaseSolid/Make/files

@@ -0,0 +1,4 @@
+phaseSolidThermophysicalTransportModel/phaseSolidThermophysicalTransportModel.C
+phaseSolidThermophysicalTransportModels.C
+
+LIB = $(FOAM_LIBBIN)/libphaseSolidThermophysicalTransportModels

src/ThermophysicalTransportModels/phaseSolid/Make/options

@@ -0,0 +1,16 @@
+EXE_INC = \
+    -I$(LIB_SRC)/ThermophysicalTransportModels/thermophysicalTransportModel/lnInclude \
+    -I$(LIB_SRC)/ThermophysicalTransportModels/solid/lnInclude \
+    -I$(LIB_SRC)/physicalProperties/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/specie/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/solidThermo/lnInclude \
+    -I$(LIB_SRC)/finiteVolume/lnInclude \
+    -I$(LIB_SRC)/meshTools/lnInclude \
+
+LIB_LIBS = \
+    -lthermophysicalTransportModel \
+    -lsolidThermo \
+    -lspecie \
+    -lfiniteVolume \
+    -lmeshTools

src/ThermophysicalTransportModels/phaseSolid/phaseSolidThermophysicalTransportModel/phaseSolidThermophysicalTransportModel.C

@@ -0,0 +1,176 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     | Website:  https://openfoam.org
+    \\  /    A nd           | Copyright (C) 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 "phaseSolidThermophysicalTransportModel.H"
+#include "isotropic.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+    defineRunTimeSelectionTable
+    (
+        phaseSolidThermophysicalTransportModel,
+        dictionary
+    );
+}
+
+
+// * * * * * * * * * * * * * Protected Member Functions  * * * * * * * * * * //
+
+void Foam::phaseSolidThermophysicalTransportModel::printCoeffs
+(const word& type)
+{
+    if (printCoeffs_)
+    {
+        Info<< coeffDict_.dictName() << coeffDict_ << endl;
+    }
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::phaseSolidThermophysicalTransportModel::
+phaseSolidThermophysicalTransportModel
+(
+    const word& type,
+    const alphaField& alpha,
+    const solidThermo& thermo
+)
+:
+    thermophysicalTransportModel(thermo.mesh(), alpha.group()),
+    alpha_(alpha),
+    thermo_(thermo),
+    printCoeffs_(lookupOrDefault<Switch>("printCoeffs", false)),
+    coeffDict_(optionalSubDict(type + "Coeffs"))
+{}
+
+
+// * * * * * * * * * * * * * * * * * Selectors * * * * * * * * * * * * * * * //
+
+Foam::autoPtr<Foam::phaseSolidThermophysicalTransportModel>
+Foam::phaseSolidThermophysicalTransportModel::New
+(
+    const alphaField& alpha,
+    const solidThermo& thermo
+)
+{
+    typeIOobject<IOdictionary> header
+    (
+        IOobject::groupName
+        (
+            phaseSolidThermophysicalTransportModel::typeName,
+            alpha.group()
+        ),
+        thermo.mesh().time().constant(),
+        thermo.mesh(),
+        IOobject::MUST_READ,
+        IOobject::NO_WRITE,
+        false
+    );
+
+    if (header.headerOk())
+    {
+        const word modelType(IOdictionary(header).lookup("model"));
+
+        Info<< "Selecting solid thermophysical transport model "
+            << modelType << endl;
+
+        typename dictionaryConstructorTable::iterator cstrIter =
+        dictionaryConstructorTablePtr_->find(modelType);
+
+        if (cstrIter == dictionaryConstructorTablePtr_->end())
+        {
+            FatalErrorInFunction
+                << "Unknown solid thermophysical transport model "
+                << modelType << nl << nl
+                << "Available models:" << endl
+                << dictionaryConstructorTablePtr_->sortedToc()
+                << exit(FatalError);
+        }
+
+        return autoPtr<phaseSolidThermophysicalTransportModel>
+        (
+            cstrIter()(alpha, thermo)
+        );
+    }
+    else
+    {
+        Info<< "Selecting default solid thermophysical transport model "
+            << solidThermophysicalTransportModels::
+               isotropic<phaseSolidThermophysicalTransportModel>::typeName
+            << endl;
+
+        return autoPtr<phaseSolidThermophysicalTransportModel>
+        (
+            new solidThermophysicalTransportModels::
+            isotropic<phaseSolidThermophysicalTransportModel>(alpha, thermo)
+        );
+    }
+}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+Foam::tmp<Foam::volScalarField>
+Foam::phaseSolidThermophysicalTransportModel::kappa() const
+{
+    return thermo().kappa();
+}
+
+
+Foam::tmp<Foam::scalarField>
+Foam::phaseSolidThermophysicalTransportModel::kappa
+(
+    const label patchi
+) const
+{
+    return thermo().kappa().boundaryField()[patchi];
+}
+
+
+bool Foam::phaseSolidThermophysicalTransportModel::read()
+{
+    if (regIOobject::read())
+    {
+        coeffDict_ <<= optionalSubDict(type() + "Coeffs");
+        return true;
+    }
+    else
+    {
+        return false;
+    }
+}
+
+
+void Foam::phaseSolidThermophysicalTransportModel::predict()
+{}
+
+
+void Foam::phaseSolidThermophysicalTransportModel::correct()
+{}
+
+
+// ************************************************************************* //

src/ThermophysicalTransportModels/phaseSolid/phaseSolidThermophysicalTransportModel/phaseSolidThermophysicalTransportModel.H

@@ -0,0 +1,195 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     | Website:  https://openfoam.org
+    \\  /    A nd           | Copyright (C) 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/>.
+
+Class
+    Foam::phaseSolidThermophysicalTransportModel
+
+Description
+    Abstract base class for solid thermophysical transport models
+
+SourceFiles
+    phaseSolidThermophysicalTransportModel.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef phaseSolidThermophysicalTransportModel_H
+#define phaseSolidThermophysicalTransportModel_H
+
+#include "thermophysicalTransportModel.H"
+#include "solidThermo.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+              Class phaseSolidThermophysicalTransportModel Declaration
+\*---------------------------------------------------------------------------*/
+
+class phaseSolidThermophysicalTransportModel
+:
+    public thermophysicalTransportModel
+{
+public:
+
+    typedef volScalarField alphaField;
+
+
+protected:
+
+    // Protected data
+
+        const alphaField& alpha_;
+
+        //- Reference to the solid thermophysical properties
+        const solidThermo& thermo_;
+
+        //- Flag to print the model coeffs at run-time
+        Switch printCoeffs_;
+
+        //- Model coefficients dictionary
+        dictionary coeffDict_;
+
+
+    // Protected Member Functions
+
+        //- Print model coefficients
+        virtual void printCoeffs(const word& type);
+
+
+public:
+
+    // Declare run-time constructor selection table
+
+        declareRunTimeSelectionTable
+        (
+            autoPtr,
+            phaseSolidThermophysicalTransportModel,
+            dictionary,
+            (
+                const alphaField& alpha,
+                const solidThermo& thermo
+            ),
+            (alpha, thermo)
+        );
+
+
+    // Constructors
+
+        //- Construct from solid thermophysical properties
+        phaseSolidThermophysicalTransportModel
+        (
+            const word& type,
+            const alphaField& alpha,
+            const solidThermo& thermo
+        );
+
+
+    // Selectors
+
+        //- Return a reference to the selected thermophysical transport model
+        static autoPtr<phaseSolidThermophysicalTransportModel> New
+        (
+            const alphaField& alpha,
+            const solidThermo& thermo
+        );
+
+
+    //- Destructor
+    virtual ~phaseSolidThermophysicalTransportModel()
+    {}
+
+
+    // Member Functions
+
+        //- Read model coefficients if they have changed
+        virtual bool read() = 0;
+
+        //- Return the phase fraction field
+        const alphaField& alpha() const
+        {
+            return alpha_;
+        }
+
+        //- Access function to solid thermophysical properties
+        virtual const solidThermo& thermo() const
+        {
+            return thermo_;
+        }
+
+        //- Const access to the coefficients dictionary
+        virtual const dictionary& coeffDict() const
+        {
+            return coeffDict_;
+        }
+
+        //- Thermal conductivity [W/m/K]
+        virtual tmp<volScalarField> kappa() const;
+
+        //- Thermal conductivity for patch [W/m/K]
+        virtual tmp<scalarField> kappa(const label patchi) const;
+
+        //- Effective thermal conductivity
+        //  of mixture [W/m/K]
+        virtual tmp<volScalarField> kappaEff() const
+        {
+            return kappa();
+        }
+
+        //- Effective thermal conductivity
+        //  of mixture for patch [W/m/K]
+        virtual tmp<scalarField> kappaEff(const label patchi) const
+        {
+            return kappa(patchi);
+        }
+
+        //- Return the heat flux [W/m^2]
+        virtual tmp<surfaceScalarField> q() const = 0;
+
+        //- Return the patch heat flux correction [W/m^2]
+        //  For isotropic or patch-aligned thermal conductivity qCorr is null
+        virtual tmp<scalarField> qCorr(const label patchi) const = 0;
+
+        //- Return the source term for the energy equation
+        virtual tmp<fvScalarMatrix> divq(volScalarField& he) const = 0;
+
+        //- Predict the thermophysical transport coefficients if possible
+        //  without solving thermophysical transport model equations
+        virtual void predict() = 0;
+
+        //- Solve the thermophysical transport model equations
+        //  and correct the thermophysical transport coefficients
+        virtual void correct();
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/ThermophysicalTransportModels/phaseSolid/phaseSolidThermophysicalTransportModels.C

@@ -0,0 +1,71 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     | Website:  https://openfoam.org
+    \\  /    A nd           | Copyright (C) 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 "phaseSolidThermophysicalTransportModel.H"
+#include "isotropic.H"
+#include "anisotropic.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace solidThermophysicalTransportModels
+{
+    typedef isotropic<phaseSolidThermophysicalTransportModel>
+        isotropicSolidThermophysicalTransportModel;
+
+    defineNamedTemplateTypeNameAndDebug
+    (
+        isotropicSolidThermophysicalTransportModel,
+        0
+    );
+
+    addToRunTimeSelectionTable
+    (
+        phaseSolidThermophysicalTransportModel,
+        isotropicSolidThermophysicalTransportModel,
+        dictionary
+    );
+
+    typedef anisotropic<phaseSolidThermophysicalTransportModel>
+        anisotropicSolidThermophysicalTransportModel;
+
+    defineNamedTemplateTypeNameAndDebug
+    (
+        anisotropicSolidThermophysicalTransportModel,
+        0
+    );
+
+    addToRunTimeSelectionTable
+    (
+        phaseSolidThermophysicalTransportModel,
+        anisotropicSolidThermophysicalTransportModel,
+        dictionary
+    );
+}
+}
+
+
+// ************************************************************************* //

src/ThermophysicalTransportModels/solid/Make/files

@@ -1,5 +1,4 @@
 solidThermophysicalTransportModel/solidThermophysicalTransportModel.C
-isotropic/isotropic.C
+solidThermophysicalTransportModels.C
-anisotropic/anisotropic.C
 
 LIB = $(FOAM_LIBBIN)/libsolidThermophysicalTransportModels

src/ThermophysicalTransportModels/solid/anisotropic/anisotropic.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
@@ -29,33 +29,18 @@ License
 #include "fvcSnGrad.H"
 #include "addToRunTimeSelectionTable.H"
 
-// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
-
-namespace Foam
-{
-namespace solidThermophysicalTransportModels
-{
-    defineTypeNameAndDebug(anisotropic, 0);
-    addToRunTimeSelectionTable
-    (
-        solidThermophysicalTransportModel,
-        anisotropic,
-        dictionary
-    );
-}
-}
-
-
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
-void Foam::solidThermophysicalTransportModels::anisotropic::
+template<class SolidThermophysicalTransportModel>
+void Foam::solidThermophysicalTransportModels::
+anisotropic<SolidThermophysicalTransportModel>::
 setZonesPatchFaces() const
 {
     if (!zoneCoordinateSystems_.size()) return;
 
     // Find all the patch faces adjacent to zones
 
-    const fvMesh& mesh = thermo().mesh();
+    const fvMesh& mesh = this->thermo().mesh();
     const fvBoundaryMesh& patches = mesh.boundary();
     const labelList& own = mesh.faceOwner();
 
@@ -110,23 +95,30 @@ setZonesPatchFaces() const
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
-Foam::solidThermophysicalTransportModels::anisotropic::anisotropic
+template<class SolidThermophysicalTransportModel>
+Foam::solidThermophysicalTransportModels::
+anisotropic<SolidThermophysicalTransportModel>::anisotropic
 (
+    const alphaField& alpha,
     const solidThermo& thermo
 )
 :
-    solidThermophysicalTransportModel(typeName, thermo),
+    SolidThermophysicalTransportModel(typeName, alpha, thermo),
     UpdateableMeshObject(*this, thermo.mesh()),
-    coordinateSystem_(coordinateSystem::New(thermo.mesh(), coeffDict())),
+    coordinateSystem_(coordinateSystem::New(thermo.mesh(), this->coeffDict())),
     boundaryAligned_
     (
-        coeffDict().lookupOrDefault<Switch>("boundaryAligned", false)
+        this->coeffDict().template lookupOrDefault<Switch>
+        (
+            "boundaryAligned",
+            false
+        )
     ),
     aligned_(thermo.mesh().boundary().size(), true)
 {
-    if (coeffDict().found("zones"))
+    if (this->coeffDict().found("zones"))
     {
-        const dictionary& zonesDict(coeffDict().subDict("zones"));
+        const dictionary& zonesDict(this->coeffDict().subDict("zones"));
 
         Info<< "    Reading coordinate system for zones:" << endl;
 
@@ -216,15 +208,19 @@ Foam::solidThermophysicalTransportModels::anisotropic::anisotropic
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
-bool Foam::solidThermophysicalTransportModels::anisotropic::read()
+template<class SolidThermophysicalTransportModel>
+bool Foam::solidThermophysicalTransportModels::
+anisotropic<SolidThermophysicalTransportModel>::read()
 {
 
     return true;
 }
 
 
+template<class SolidThermophysicalTransportModel>
 Foam::tmp<Foam::volSymmTensorField>
-Foam::solidThermophysicalTransportModels::anisotropic::Kappa() const
+Foam::solidThermophysicalTransportModels::
+anisotropic<SolidThermophysicalTransportModel>::Kappa() const
 {
     const solidThermo& thermo = this->thermo();
     const fvMesh& mesh = thermo.mesh();
@@ -308,8 +304,10 @@ Foam::solidThermophysicalTransportModels::anisotropic::Kappa() const
 }
 
 
+template<class SolidThermophysicalTransportModel>
 Foam::tmp<Foam::symmTensorField>
-Foam::solidThermophysicalTransportModels::anisotropic::Kappa
+Foam::solidThermophysicalTransportModels::
+anisotropic<SolidThermophysicalTransportModel>::Kappa
 (
     const label patchi
 ) const
@@ -356,27 +354,33 @@ Foam::solidThermophysicalTransportModels::anisotropic::Kappa
 }
 
 
+template<class SolidThermophysicalTransportModel>
 Foam::tmp<Foam::volScalarField>
-Foam::solidThermophysicalTransportModels::anisotropic::kappa() const
+Foam::solidThermophysicalTransportModels::
+anisotropic<SolidThermophysicalTransportModel>::kappa() const
 {
     NotImplemented;
     return tmp<volScalarField>(nullptr);
 }
 
 
+template<class SolidThermophysicalTransportModel>
 Foam::tmp<Foam::scalarField>
-Foam::solidThermophysicalTransportModels::anisotropic::kappa
+Foam::solidThermophysicalTransportModels::
+anisotropic<SolidThermophysicalTransportModel>::kappa
 (
     const label patchi
 ) const
 {
-    const vectorField n(thermo().mesh().boundary()[patchi].nf());
+    const vectorField n(this->thermo().mesh().boundary()[patchi].nf());
     return n & Kappa(patchi) & n;
 }
 
 
+template<class SolidThermophysicalTransportModel>
 Foam::tmp<Foam::surfaceScalarField>
-Foam::solidThermophysicalTransportModels::anisotropic::q() const
+Foam::solidThermophysicalTransportModels::
+anisotropic<SolidThermophysicalTransportModel>::q() const
 {
     const solidThermo& thermo = this->thermo();
     const fvMesh& mesh = thermo.mesh();
@@ -384,25 +388,31 @@ Foam::solidThermophysicalTransportModels::anisotropic::q() const
     return surfaceScalarField::New
     (
         "q",
-       -fvm::laplacian(Kappa(), thermo.T())().flux()/mesh.magSf()
+       -fvm::laplacian(this->alpha()*Kappa(), thermo.T())().flux()/mesh.magSf()
     );
 }
 
 
+template<class SolidThermophysicalTransportModel>
 Foam::tmp<Foam::scalarField>
-Foam::solidThermophysicalTransportModels::anisotropic::qCorr
+Foam::solidThermophysicalTransportModels::
+anisotropic<SolidThermophysicalTransportModel>::qCorr
 (
     const label patchi
 ) const
 {
     if (!aligned_[patchi])
     {
-        tmp<volVectorField> gradT(fvc::grad(thermo().T()));
+        tmp<volVectorField> gradT(fvc::grad(this->thermo().T()));
 
-        const vectorField n(thermo().mesh().boundary()[patchi].nf());
+        const vectorField n(this->thermo().mesh().boundary()[patchi].nf());
         const vectorField nKappa(n & Kappa(patchi));
 
-        return -(nKappa - n*(nKappa & n)) & gradT().boundaryField()[patchi];
+        return
+           -(
+                this->alpha().boundaryField()[patchi]
+               *((nKappa - n*(nKappa & n)) & gradT().boundaryField()[patchi])
+            );
     }
     else
     {
@@ -411,8 +421,10 @@ Foam::solidThermophysicalTransportModels::anisotropic::qCorr
 }
 
 
+template<class SolidThermophysicalTransportModel>
 Foam::tmp<Foam::fvScalarMatrix>
-Foam::solidThermophysicalTransportModels::anisotropic::divq
+Foam::solidThermophysicalTransportModels::
+anisotropic<SolidThermophysicalTransportModel>::divq
 (
     volScalarField& e
 ) const
@@ -425,18 +437,21 @@ Foam::solidThermophysicalTransportModels::anisotropic::divq
     // Return heat flux source as an implicit energy correction
     // to the temperature gradient flux
     return
-       -fvc::laplacian(Kappa, thermo.T())
+       -fvc::laplacian(this->alpha()*Kappa, thermo.T())
        -fvm::laplacianCorrection
         (
-            (Sf & fvc::interpolate(Kappa/thermo.Cv()) & Sf)/sqr(magSf),
+            (Sf & fvc::interpolate(this->alpha()*Kappa/thermo.Cv()) & Sf)
+           /sqr(magSf),
             e
         );
 }
 
 
-void Foam::solidThermophysicalTransportModels::anisotropic::predict()
+template<class SolidThermophysicalTransportModel>
+void Foam::solidThermophysicalTransportModels::
+anisotropic<SolidThermophysicalTransportModel>::predict()
 {
-    solidThermophysicalTransportModel::predict();
+    SolidThermophysicalTransportModel::predict();
 
     // Recalculate zonesPatchFaces if they have been deleted
     // following mesh changes
@@ -447,13 +462,17 @@ void Foam::solidThermophysicalTransportModels::anisotropic::predict()
 }
 
 
-bool Foam::solidThermophysicalTransportModels::anisotropic::movePoints()
+template<class SolidThermophysicalTransportModel>
+bool Foam::solidThermophysicalTransportModels::
+anisotropic<SolidThermophysicalTransportModel>::movePoints()
 {
     return true;
 }
 
 
-void Foam::solidThermophysicalTransportModels::anisotropic::topoChange
+template<class SolidThermophysicalTransportModel>
+void Foam::solidThermophysicalTransportModels::
+anisotropic<SolidThermophysicalTransportModel>::topoChange
 (
     const polyTopoChangeMap& map
 )
@@ -463,7 +482,9 @@ void Foam::solidThermophysicalTransportModels::anisotropic::topoChange
 }
 
 
-void Foam::solidThermophysicalTransportModels::anisotropic::mapMesh
+template<class SolidThermophysicalTransportModel>
+void Foam::solidThermophysicalTransportModels::
+anisotropic<SolidThermophysicalTransportModel>::mapMesh
 (
     const polyMeshMap& map
 )
@@ -473,7 +494,9 @@ void Foam::solidThermophysicalTransportModels::anisotropic::mapMesh
 }
 
 
-void Foam::solidThermophysicalTransportModels::anisotropic::distribute
+template<class SolidThermophysicalTransportModel>
+void Foam::solidThermophysicalTransportModels::
+anisotropic<SolidThermophysicalTransportModel>::distribute
 (
     const polyDistributionMap& map
 )

src/ThermophysicalTransportModels/solid/anisotropic/anisotropic.H

@@ -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
@@ -114,9 +114,10 @@ namespace solidThermophysicalTransportModels
                            Class anisotropic Declaration
 \*---------------------------------------------------------------------------*/
 
+template<class SolidThermophysicalTransportModel>
 class anisotropic
 :
-    public solidThermophysicalTransportModel,
+    public SolidThermophysicalTransportModel,
     public UpdateableMeshObject<fvMesh>
 {
     // Private member data
@@ -150,6 +151,9 @@ class anisotropic
 
 public:
 
+    typedef typename SolidThermophysicalTransportModel::alphaField
+        alphaField;
+
     //- Runtime type information
     TypeName("anisotropic");
 
@@ -157,7 +161,11 @@ public:
     // Constructors
 
         //- Construct from solid thermophysical properties
-        anisotropic(const solidThermo& thermo);
+        anisotropic
+        (
+            const alphaField& alpha,
+            const solidThermo& thermo
+        );
 
 
     //- Destructor
@@ -213,6 +221,12 @@ public:
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
+#ifdef NoRepository
+    #include "anisotropic.C"
+#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
 #endif
 
 // ************************************************************************* //

src/ThermophysicalTransportModels/solid/isotropic/isotropic.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
@@ -29,31 +29,17 @@ License
 #include "fvcSnGrad.H"
 #include "addToRunTimeSelectionTable.H"
 
-// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
-
-namespace Foam
-{
-namespace solidThermophysicalTransportModels
-{
-    defineTypeNameAndDebug(isotropic, 0);
-    addToRunTimeSelectionTable
-    (
-        solidThermophysicalTransportModel,
-        isotropic,
-        dictionary
-    );
-}
-}
-
-
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
-Foam::solidThermophysicalTransportModels::isotropic::isotropic
+template<class SolidThermophysicalTransportModel>
+Foam::solidThermophysicalTransportModels::
+isotropic<SolidThermophysicalTransportModel>::isotropic
 (
+    const alphaField& alpha,
     const solidThermo& thermo
 )
 :
-    solidThermophysicalTransportModel(typeName, thermo)
+    SolidThermophysicalTransportModel(typeName, alpha, thermo)
 {
     if (!thermo.isotropic())
     {
@@ -67,32 +53,41 @@ Foam::solidThermophysicalTransportModels::isotropic::isotropic
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
+template<class SolidThermophysicalTransportModel>
 const Foam::dictionary&
-Foam::solidThermophysicalTransportModels::isotropic::coeffDict() const
+Foam::solidThermophysicalTransportModels::
+isotropic<SolidThermophysicalTransportModel>::coeffDict() const
 {
     return dictionary::null;
 }
 
 
-bool Foam::solidThermophysicalTransportModels::isotropic::read()
+template<class SolidThermophysicalTransportModel>
+bool Foam::solidThermophysicalTransportModels::
+isotropic<SolidThermophysicalTransportModel>::read()
 {
     return true;
 }
 
 
+template<class SolidThermophysicalTransportModel>
 Foam::tmp<Foam::surfaceScalarField>
-Foam::solidThermophysicalTransportModels::isotropic::q() const
+Foam::solidThermophysicalTransportModels::
+isotropic<SolidThermophysicalTransportModel>::q() const
 {
     return surfaceScalarField::New
     (
         "q",
-       -fvc::interpolate(kappa())*fvc::snGrad(thermo().T())
+       -fvc::interpolate(this->alpha()*this->kappa())
+       *fvc::snGrad(this->thermo().T())
     );
 }
 
 
+template<class SolidThermophysicalTransportModel>
 Foam::tmp<Foam::scalarField>
-Foam::solidThermophysicalTransportModels::isotropic::qCorr
+Foam::solidThermophysicalTransportModels::
+isotropic<SolidThermophysicalTransportModel>::qCorr
 (
     const label patchi
 ) const
@@ -101,8 +96,10 @@ Foam::solidThermophysicalTransportModels::isotropic::qCorr
 }
 
 
+template<class SolidThermophysicalTransportModel>
 Foam::tmp<Foam::fvScalarMatrix>
-Foam::solidThermophysicalTransportModels::isotropic::divq
+Foam::solidThermophysicalTransportModels::
+isotropic<SolidThermophysicalTransportModel>::divq
 (
     volScalarField& e
 ) const
@@ -112,14 +109,16 @@ Foam::solidThermophysicalTransportModels::isotropic::divq
     // Return heat flux source as an implicit energy correction
     // to the temperature gradient flux
     return
-       -fvc::laplacian(kappa(), thermo.T())
+       -fvc::laplacian(this->alpha()*this->kappa(), thermo.T())
-       -fvm::laplacianCorrection(kappa()/thermo.Cv(), e);
+       -fvm::laplacianCorrection(this->alpha()*this->kappa()/thermo.Cv(), e);
 }
 
 
-void Foam::solidThermophysicalTransportModels::isotropic::predict()
+template<class SolidThermophysicalTransportModel>
+void Foam::solidThermophysicalTransportModels::
+isotropic<SolidThermophysicalTransportModel>::predict()
 {
-    solidThermophysicalTransportModel::predict();
+    SolidThermophysicalTransportModel::predict();
 }
 
 

src/ThermophysicalTransportModels/solid/isotropic/isotropic.H

@@ -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
@@ -52,13 +52,17 @@ namespace solidThermophysicalTransportModels
                            Class isotropic Declaration
 \*---------------------------------------------------------------------------*/
 
+template<class SolidThermophysicalTransportModel>
 class isotropic
 :
-    public solidThermophysicalTransportModel
+    public SolidThermophysicalTransportModel
 {
 
 public:
 
+    typedef typename SolidThermophysicalTransportModel::alphaField
+        alphaField;
+
     //- Runtime type information
     TypeName("isotropic");
 
@@ -66,7 +70,11 @@ public:
     // Constructors
 
         //- Construct from solid thermophysical properties
-        isotropic(const solidThermo& thermo);
+        isotropic
+        (
+            const alphaField& alpha,
+            const solidThermo& thermo
+        );
 
 
     //- Destructor
@@ -103,6 +111,12 @@ public:
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
+#ifdef NoRepository
+    #include "isotropic.C"
+#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
 #endif
 
 // ************************************************************************* //

src/ThermophysicalTransportModels/solid/solidThermophysicalTransportModel/solidThermophysicalTransportModel.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
@@ -51,10 +51,12 @@ void Foam::solidThermophysicalTransportModel::printCoeffs
 Foam::solidThermophysicalTransportModel::solidThermophysicalTransportModel
 (
     const word& type,
+    const alphaField& alpha,
     const solidThermo& thermo
 )
 :
     thermophysicalTransportModel(thermo.mesh(), word::null),
+    alpha_(alpha),
     thermo_(thermo),
     printCoeffs_(lookupOrDefault<Switch>("printCoeffs", false)),
     coeffDict_(optionalSubDict(type + "Coeffs"))
@@ -98,18 +100,24 @@ Foam::solidThermophysicalTransportModel::New(const solidThermo& thermo)
 
         return autoPtr<solidThermophysicalTransportModel>
         (
-            cstrIter()(thermo)
+            cstrIter()(geometricOneField(), thermo)
         );
     }
     else
     {
         Info<< "Selecting default solid thermophysical transport model "
-            << solidThermophysicalTransportModels::isotropic::typeName
+            << solidThermophysicalTransportModels::
+               isotropic<solidThermophysicalTransportModel>::typeName
             << endl;
 
         return autoPtr<solidThermophysicalTransportModel>
         (
-            new solidThermophysicalTransportModels::isotropic(thermo)
+            new solidThermophysicalTransportModels::
+            isotropic<solidThermophysicalTransportModel>
+            (
+                geometricOneField(),
+                thermo
+            )
         );
     }
 }

src/ThermophysicalTransportModels/solid/solidThermophysicalTransportModel/solidThermophysicalTransportModel.H

@@ -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
@@ -51,10 +51,17 @@ class solidThermophysicalTransportModel
 :
     public thermophysicalTransportModel
 {
+public:
+
+    typedef geometricOneField alphaField;
+
+
 protected:
 
     // Protected data
 
+        alphaField alpha_;
+
         //- Reference to the solid thermophysical properties
         const solidThermo& thermo_;
 
@@ -81,9 +88,10 @@ public:
             solidThermophysicalTransportModel,
             dictionary,
             (
+                const alphaField& alpha,
                 const solidThermo& thermo
             ),
-            (thermo)
+            (alpha, thermo)
         );
 
 
@@ -93,6 +101,7 @@ public:
         solidThermophysicalTransportModel
         (
             const word& type,
+            const alphaField& alpha,
             const solidThermo& thermo
         );
 
@@ -116,6 +125,12 @@ public:
         //- Read model coefficients if they have changed
         virtual bool read() = 0;
 
+        //- Return the phase fraction field
+        const alphaField& alpha() const
+        {
+            return alpha_;
+        }
+
         //- Access function to solid thermophysical properties
         virtual const solidThermo& thermo() const
         {

src/ThermophysicalTransportModels/solid/solidThermophysicalTransportModels.C

@@ -0,0 +1,71 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     | Website:  https://openfoam.org
+    \\  /    A nd           | Copyright (C) 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 "solidThermophysicalTransportModel.H"
+#include "isotropic.H"
+#include "anisotropic.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace solidThermophysicalTransportModels
+{
+    typedef isotropic<solidThermophysicalTransportModel>
+        isotropicSolidThermophysicalTransportModel;
+
+    defineNamedTemplateTypeNameAndDebug
+    (
+        isotropicSolidThermophysicalTransportModel,
+        0
+    );
+
+    addToRunTimeSelectionTable
+    (
+        solidThermophysicalTransportModel,
+        isotropicSolidThermophysicalTransportModel,
+        dictionary
+    );
+
+    typedef anisotropic<solidThermophysicalTransportModel>
+        anisotropicSolidThermophysicalTransportModel;
+
+    defineNamedTemplateTypeNameAndDebug
+    (
+        anisotropicSolidThermophysicalTransportModel,
+        0
+    );
+
+    addToRunTimeSelectionTable
+    (
+        solidThermophysicalTransportModel,
+        anisotropicSolidThermophysicalTransportModel,
+        dictionary
+    );
+}
+}
+
+
+// ************************************************************************* //

tutorials/multiphaseEuler/bed/constant/physicalProperties.solid

@@ -16,10 +16,10 @@ FoamFile
 
 thermoType
 {
-    type            heRhoThermo;
+    type            heSolidThermo;
     mixture         pureMixture;
-    transport       const;
+    transport       constIsoSolid;
-    thermo          hConst;
+    thermo          eConst;
     equationOfState rhoConst;
     specie          specie;
     energy          sensibleInternalEnergy;
@@ -37,13 +37,12 @@ mixture
     }
     thermodynamics
     {
-        Cp          6000;
+        Cv          6000;
         Hf          0;
     }
     transport
     {
-        mu          0;
+        kappa       0;
-        Pr          1;
     }
 }
 

tutorials/multiphaseEuler/boilingBed/0/T.solid

@@ -30,7 +30,8 @@ boundaryField
     }
     wall
     {
-        type            zeroGradient;
+        type            fixedValue;
+        value           uniform 350;
     }
     front
     {

tutorials/multiphaseEuler/boilingBed/0/wallHeatFlux.solid

@@ -0,0 +1,49 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     | Website:  https://openfoam.org
+    \\  /    A nd           | Version:  dev
+     \\/     M anipulation  |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    format      ascii;
+    class       volScalarField;
+    location    "0";
+    object      wallHeatFlux.solid;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [1 0 -3 0 0 0 0];
+
+internalField   uniform 0;
+
+boundaryField
+{
+    inlet
+    {
+        type            calculated;
+        value           uniform 0;
+    }
+    outlet
+    {
+        type            calculated;
+        value           uniform 0;
+    }
+    front
+    {
+        type            wedge;
+    }
+    back
+    {
+        type            wedge;
+    }
+    wall
+    {
+        type            calculated;
+        value           uniform 0;
+    }
+}
+
+
+// ************************************************************************* //

tutorials/multiphaseEuler/boilingBed/constant/physicalProperties.solid

@@ -10,15 +10,15 @@ FoamFile
     format      ascii;
     class       dictionary;
     location    "constant";
-    object      thermophysicalProperties.solid;
+    object      physicalProperties.solid;
 }
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 thermoType
 {
-    type            heRhoThermo;
+    type            heSolidThermo;
     mixture         pureMixture;
-    transport       const;
+    transport       constIsoSolid;
     thermo          eConst;
     equationOfState rhoConst;
     specie          specie;
@@ -37,12 +37,11 @@ mixture
     }
     thermodynamics
     {
-        Hf          0;
         Cv          520;
+        Hf          0;
     }
     transport
     {
-        mu          1;
         kappa       16.7;
     }
 }

tutorials/multiphaseEuler/boilingBed/system/controlDict

@@ -164,6 +164,14 @@ functions
         operation       volIntegrate;
         fields          (continuityError.gas continuityError.liquid);
     }
+
+    writeWallheatfluxSolid
+    {
+        type            wallHeatFlux;
+        libs            ( "libmultiphaseEulerFoamFunctionObjects.so" );
+        writeControl    writeTime;
+        phase           solid;
+    }
 }