Merge branch 'integration-foundation' into 'develop'

INT: openfoam.org code integration

See merge request Development/openfoam!312

applications/utilities/postProcessing/lagrangian/particleTracks/particleTracks.C

@@ -69,7 +69,7 @@ int main(int argc, char *argv[])
 
     // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-    fileName vtkPath(runTime.path()/"VTK");
+    const fileName vtkPath(runTime.rootPath()/runTime.globalCaseName()/"VTK");
     mkDir(vtkPath);
 
     Info<< "Scanning times to determine track data for cloud " << cloudName

applications/utilities/postProcessing/lagrangian/steadyParticleTracks/steadyParticleTracks.C

@@ -137,7 +137,7 @@ int main(int argc, char *argv[])
 
     // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-    fileName vtkPath(runTime.path()/"VTK");
+    const fileName vtkPath(runTime.rootPath()/runTime.globalCaseName()/"VTK");
     mkDir(vtkPath);
 
     forAll(timeDirs, timeI)
@@ -145,7 +145,7 @@ int main(int argc, char *argv[])
         runTime.setTime(timeDirs[timeI], timeI);
         Info<< "Time = " << runTime.timeName() << endl;
 
-        fileName vtkTimePath(runTime.path()/"VTK"/runTime.timeName());
+        const fileName vtkTimePath(vtkPath/runTime.timeName());
         mkDir(vtkTimePath);
 
         Info<< "    Reading particle positions" << endl;

src/finiteVolume/cfdTools/general/porosityModel/porosityModel/porosityModel.C

@@ -42,12 +42,13 @@ namespace Foam
 
 void Foam::porosityModel::adjustNegativeResistance(dimensionedVector& resist)
 {
-    scalar maxCmpt = max(0, cmptMax(resist.value()));
+    scalar maxCmpt = cmptMax(resist.value());
 
     if (maxCmpt < 0)
     {
         FatalErrorInFunction
-            << "Negative resistances are invalid, resistance = " << resist
+            << "Cannot have all resistances set to negative, resistance = "
+            << resist
             << exit(FatalError);
     }
     else

src/lagrangian/intermediate/clouds/Templates/MPPICCloud/MPPICCloud.C

@@ -202,6 +202,13 @@ void Foam::MPPICCloud<CloudType>::motion
 
     if (dampingModel_->active())
     {
+        if (this->mesh().moving())
+        {
+            FatalErrorInFunction
+                << "MPPIC damping modelling does not support moving meshes."
+                << exit(FatalError);
+        }
+
         // update averages
         td.updateAverages(cloud);
 
@@ -226,6 +233,13 @@ void Foam::MPPICCloud<CloudType>::motion
 
     if (packingModel_->active())
     {
+        if (this->mesh().moving())
+        {
+            FatalErrorInFunction
+                << "MPPIC packing modelling does not support moving meshes."
+                << exit(FatalError);
+        }
+
         // same procedure as for damping
         td.updateAverages(cloud);
         packingModel_->cacheFields(true);

src/phaseSystemModels/reactingEulerFoam/interfacialModels/dragModels/segregated/segregated.C

@@ -109,7 +109,7 @@ Foam::tmp<Foam::volScalarField> Foam::dragModels::segregated::K() const
     L.primitiveFieldRef() = cbrt(mesh.V());
     L.correctBoundaryConditions();
 
-    volScalarField I
+    const volScalarField I
     (
         alpha1
        /max
@@ -118,7 +118,7 @@ Foam::tmp<Foam::volScalarField> Foam::dragModels::segregated::K() const
             pair_.phase1().residualAlpha() + pair_.phase2().residualAlpha()
         )
     );
-    volScalarField magGradI
+    const volScalarField magGradI
     (
         max
         (
@@ -127,28 +127,36 @@ Foam::tmp<Foam::volScalarField> Foam::dragModels::segregated::K() const
         )
     );
 
-    volScalarField muI
+    const volScalarField muI
     (
         rho1*nu1*rho2*nu2
        /(rho1*nu1 + rho2*nu2)
     );
-    volScalarField muAlphaI
+
+    const volScalarField limitedAlpha1
     (
-        alpha1*rho1*nu1*alpha2*rho2*nu2
-       /(
-           max(alpha1, pair_.phase1().residualAlpha())*rho1*nu1
-         + max(alpha2, pair_.phase2().residualAlpha())*rho2*nu2
-        )
+        max(alpha1, pair_.phase1().residualAlpha())
     );
 
-    volScalarField ReI
+    const volScalarField limitedAlpha2
+    (
+        max(alpha2, pair_.phase2().residualAlpha())
+    );
+
+    const volScalarField muAlphaI
+    (
+        alpha1*rho1*nu1*alpha2*rho2*nu2
+       /(limitedAlpha1*rho1*nu1 + limitedAlpha2*rho2*nu2)
+    );
+
+    const volScalarField ReI
     (
         pair_.rho()
        *pair_.magUr()
-       /(magGradI*muI)
+       /(magGradI*limitedAlpha1*limitedAlpha2*muI)
     );
 
-    volScalarField lambda(m_*ReI + n_*muAlphaI/muI);
+    const volScalarField lambda(m_*ReI + n_*muAlphaI/muI);
 
     return lambda*sqr(magGradI)*muI;
 }

src/phaseSystemModels/reactingEulerFoam/phaseSystems/reactionThermo/hRefConstThermos.C

@@ -33,6 +33,7 @@ License
 
 #include "specie.H"
 #include "perfectGas.H"
+#include "rPolynomial.H"
 #include "perfectFluid.H"
 #include "rhoConst.H"
 
@@ -84,6 +85,19 @@ constTransport
     >
 > constRefFluidHThermoPhysics;
 
+typedef
+constTransport
+<
+    species::thermo
+    <
+        hRefConstThermo
+        <
+            rPolynomial<specie>
+        >,
+        sensibleEnthalpy
+    >
+> constRefrPolFluidHThermoPhysics;
+
 typedef
 constTransport
 <
@@ -110,6 +124,19 @@ constTransport
     >
 > constRefFluidEThermoPhysics;
 
+typedef
+constTransport
+<
+    species::thermo
+    <
+        eRefConstThermo
+        <
+            rPolynomial<specie>
+        >,
+        sensibleInternalEnergy
+    >
+> constRefrPolFluidEThermoPhysics;
+
 typedef
 constTransport
 <
@@ -151,6 +178,18 @@ makeThermos
     specie
 );
 
+makeThermos
+(
+    rhoThermo,
+    heRhoThermo,
+    pureMixture,
+    constTransport,
+    sensibleEnthalpy,
+    hRefConstThermo,
+    rPolynomial,
+    specie
+);
+
 makeThermos
 (
     rhoThermo,
@@ -190,6 +229,18 @@ makeThermos
     specie
 );
 
+makeThermos
+(
+    rhoThermo,
+    heRhoThermo,
+    pureMixture,
+    constTransport,
+    sensibleInternalEnergy,
+    eRefConstThermo,
+    rPolynomial,
+    specie
+);
+
 makeThermos
 (
     rhoThermo,
@@ -235,6 +286,15 @@ makeThermoPhysicsReactionThermos
     constRefFluidEThermoPhysics
 );
 
+makeThermoPhysicsReactionThermos
+(
+    rhoThermo,
+    rhoReactionThermo,
+    heRhoThermo,
+    multiComponentMixture,
+    constRefrPolFluidEThermoPhysics
+);
+
 makeThermoPhysicsReactionThermos
 (
     rhoThermo,
@@ -265,6 +325,15 @@ makeThermoPhysicsReactionThermos
     constRefFluidHThermoPhysics
 );
 
+makeThermoPhysicsReactionThermos
+(
+    rhoThermo,
+    rhoReactionThermo,
+    heRhoThermo,
+    multiComponentMixture,
+    constRefrPolFluidHThermoPhysics
+);
+
 makeThermoPhysicsReactionThermos
 (
     rhoThermo,

src/rigidBodyDynamics/Make/files

@@ -33,6 +33,7 @@ restraints/linearDamper/linearDamper.C
 restraints/linearAxialAngularSpring/linearAxialAngularSpring.C
 restraints/sphericalAngularDamper/sphericalAngularDamper.C
 restraints/prescribedRotation/prescribedRotation.C
+restraints/externalForce/externalForce.C
 restraints/softWall/softWall.C
 
 rigidBodyModel/rigidBodyModel.C

src/rigidBodyDynamics/restraints/externalForce/externalForce.C

@@ -0,0 +1,130 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2019 OpenFOAM Foundation
+-------------------------------------------------------------------------------
+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 "externalForce.H"
+#include "rigidBodyModel.H"
+#include "rigidBodyModelState.H"
+#include "OneConstant.H"
+#include "addToRunTimeSelectionTable.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace RBD
+{
+namespace restraints
+{
+    defineTypeNameAndDebug(externalForce, 0);
+
+    addToRunTimeSelectionTable
+    (
+        restraint,
+        externalForce,
+        dictionary
+    );
+}
+}
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::RBD::restraints::externalForce::externalForce
+(
+    const word& name,
+    const dictionary& dict,
+    const rigidBodyModel& model
+)
+:
+    restraint(name, dict, model),
+    externalForce_(nullptr)
+{
+    read(dict);
+}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+Foam::RBD::restraints::externalForce::~externalForce()
+{}
+
+
+// * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
+
+void Foam::RBD::restraints::externalForce::restrain
+(
+    scalarField& tau,
+    Field<spatialVector>& fx,
+    const rigidBodyModelState& state
+) const
+{
+    const vector force = externalForce_().value(state.t());
+    const vector moment(location_ ^ force);
+
+    if (model_.debug)
+    {
+        Info<< " location " << location_
+            << " force " << force
+            << " moment " << moment
+            << endl;
+    }
+
+    // Accumulate the force for the restrained body
+    fx[bodyIndex_] += spatialVector(moment, force);
+}
+
+
+bool Foam::RBD::restraints::externalForce::read
+(
+    const dictionary& dict
+)
+{
+    restraint::read(dict);
+
+    coeffs_.lookup("location") >> location_;
+
+    externalForce_ = Function1<vector>::New("force", coeffs_);
+
+    return true;
+}
+
+
+void Foam::RBD::restraints::externalForce::write
+(
+    Ostream& os
+) const
+{
+    restraint::write(os);
+
+    os.writeEntry("location", location_);
+
+    externalForce_().writeData(os);
+}
+
+
+// ************************************************************************* //

src/rigidBodyDynamics/restraints/externalForce/externalForce.H

@@ -0,0 +1,146 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2019 OpenFOAM Foundation
+-------------------------------------------------------------------------------
+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::RBD::restraints::externalForce
+
+Group
+    grpRigidBodyDynamicsRestraints
+
+Description
+    Time-dependent external force restraint using Function1.
+
+Usage
+    Example applying a constant force to the floatingObject:
+    \verbatim
+    restraints
+    {
+        force
+        {
+            type        externalForce;
+            body        floatingObject;
+            location    (0 0 0);
+            force       (100 0 0);
+        }
+    }
+    \endverbatim
+
+SourceFiles
+    externalForce.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef RBD_restraints_externalForce_H
+#define RBD_restraints_externalForce_H
+
+#include "rigidBodyRestraint.H"
+#include "Function1.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace RBD
+{
+namespace restraints
+{
+
+/*---------------------------------------------------------------------------*\
+                          Class externalForce Declaration
+\*---------------------------------------------------------------------------*/
+
+class externalForce
+:
+    public restraint
+{
+    // Private data
+
+        //- External force (N)
+        autoPtr<Function1<vector>> externalForce_;
+
+        //- Point of application of the force
+        vector location_;
+
+
+public:
+
+    //- Runtime type information
+    TypeName("externalForce");
+
+
+    // Constructors
+
+        //- Construct from components
+        externalForce
+        (
+            const word& name,
+            const dictionary& dict,
+            const rigidBodyModel& model
+        );
+
+        //- Construct and return a clone
+        virtual autoPtr<restraint> clone() const
+        {
+            return autoPtr<restraint>
+            (
+                new externalForce(*this)
+            );
+        }
+
+
+    //- Destructor
+    virtual ~externalForce();
+
+
+    // Member Functions
+
+        //- Accumulate the retraint internal joint forces into the tau field and
+        //  external forces into the fx field
+        virtual void restrain
+        (
+            scalarField& tau,
+            Field<spatialVector>& fx,
+            const rigidBodyModelState& state
+        ) const;
+
+        //- Update properties from given dictionary
+        virtual bool read(const dictionary& dict);
+
+        //- Write
+        virtual void write(Ostream&) const;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace restraints
+} // End namespace RBD
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/rigidBodyDynamics/restraints/linearAxialAngularSpring/linearAxialAngularSpring.C

@@ -77,7 +77,8 @@ Foam::RBD::restraints::linearAxialAngularSpring::~linearAxialAngularSpring()
 void Foam::RBD::restraints::linearAxialAngularSpring::restrain
 (
     scalarField& tau,
-    Field<spatialVector>& fx
+    Field<spatialVector>& fx,
+    const rigidBodyModelState& state
 ) const
 {
     vector refDir = rotationTensor(vector(1, 0, 0), axis_) & vector(0, 1, 0);

src/rigidBodyDynamics/restraints/linearAxialAngularSpring/linearAxialAngularSpring.H

@@ -111,7 +111,8 @@ public:
         virtual void restrain
         (
             scalarField& tau,
-            Field<spatialVector>& fx
+            Field<spatialVector>& fx,
+            const rigidBodyModelState& state
         ) const;
 
         //- Update properties from given dictionary

src/rigidBodyDynamics/restraints/linearDamper/linearDamper.C

@@ -76,7 +76,8 @@ Foam::RBD::restraints::linearDamper::~linearDamper()
 void Foam::RBD::restraints::linearDamper::restrain
 (
     scalarField& tau,
-    Field<spatialVector>& fx
+    Field<spatialVector>& fx,
+    const rigidBodyModelState& state
 ) const
 {
     vector force = -coeff_*model_.v(model_.master(bodyID_)).l();

src/rigidBodyDynamics/restraints/linearDamper/linearDamper.H

@@ -102,7 +102,8 @@ public:
         virtual void restrain
         (
             scalarField& tau,
-            Field<spatialVector>& fx
+            Field<spatialVector>& fx,
+            const rigidBodyModelState& state
         ) const;
 
         //- Update properties from given dictionary

src/rigidBodyDynamics/restraints/linearSpring/linearSpring.C

@@ -76,7 +76,8 @@ Foam::RBD::restraints::linearSpring::~linearSpring()
 void Foam::RBD::restraints::linearSpring::restrain
 (
     scalarField& tau,
-    Field<spatialVector>& fx
+    Field<spatialVector>& fx,
+    const rigidBodyModelState& state
 ) const
 {
     point attachmentPt = bodyPoint(refAttachmentPt_);

src/rigidBodyDynamics/restraints/linearSpring/linearSpring.H

@@ -115,7 +115,8 @@ public:
         virtual void restrain
         (
             scalarField& tau,
-            Field<spatialVector>& fx
+            Field<spatialVector>& fx,
+            const rigidBodyModelState& state
         ) const;
 
         //- Update properties from given dictionary

src/rigidBodyDynamics/restraints/prescribedRotation/prescribedRotation.C

@@ -78,7 +78,8 @@ Foam::RBD::restraints::prescribedRotation::~prescribedRotation()
 void Foam::RBD::restraints::prescribedRotation::restrain
 (
     scalarField& tau,
-    Field<spatialVector>& fx
+    Field<spatialVector>& fx,
+    const rigidBodyModelState& state
 ) const
 {
     vector refDir = rotationTensor(vector(1, 0, 0), axis_) & vector(0, 1, 0);

src/rigidBodyDynamics/restraints/prescribedRotation/prescribedRotation.H

@@ -123,7 +123,8 @@ public:
         virtual void restrain
         (
             scalarField& tau,
-            Field<spatialVector>& fx
+            Field<spatialVector>& fx,
+            const rigidBodyModelState& state
         ) const;
 
         //- Update properties from given dictionary

src/rigidBodyDynamics/restraints/restraint/rigidBodyRestraint.H

@@ -50,6 +50,7 @@ SourceFiles
 #include "point.H"
 #include "scalarField.H"
 #include "runTimeSelectionTables.H"
+#include "rigidBodyModelState.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -164,7 +165,8 @@ public:
         virtual void restrain
         (
             scalarField& tau,
-            Field<spatialVector>& fx
+            Field<spatialVector>& fx,
+            const rigidBodyModelState& state
         ) const = 0;
 
         //- Update properties from given dictionary

src/rigidBodyDynamics/restraints/softWall/softWall.C

@@ -77,7 +77,8 @@ Foam::RBD::restraints::softWall::~softWall()
 void Foam::RBD::restraints::softWall::restrain
 (
     scalarField& tau,
-    Field<spatialVector>& fx
+    Field<spatialVector>& fx,
+    const rigidBodyModelState& state
 ) const
 {
 

src/rigidBodyDynamics/restraints/softWall/softWall.H

@@ -121,7 +121,8 @@ public:
         virtual void restrain
         (
             scalarField& tau,
-            Field<spatialVector>& fx
+            Field<spatialVector>& fx,
+            const rigidBodyModelState& state
         ) const;
 
         //- Update properties from given dictionary

src/rigidBodyDynamics/restraints/sphericalAngularDamper/sphericalAngularDamper.C

@@ -76,7 +76,8 @@ Foam::RBD::restraints::sphericalAngularDamper::~sphericalAngularDamper()
 void Foam::RBD::restraints::sphericalAngularDamper::restrain
 (
     scalarField& tau,
-    Field<spatialVector>& fx
+    Field<spatialVector>& fx,
+    const rigidBodyModelState& state
 ) const
 {
     vector moment = -coeff_*model_.v(model_.master(bodyID_)).w();

src/rigidBodyDynamics/restraints/sphericalAngularDamper/sphericalAngularDamper.H

@@ -103,7 +103,8 @@ public:
         virtual void restrain
         (
             scalarField& tau,
-            Field<spatialVector>& fx
+            Field<spatialVector>& fx,
+            const rigidBodyModelState& state
         ) const;
 
         //- Update properties from given dictionary

src/rigidBodyDynamics/rigidBodyModel/forwardDynamics.C

@@ -34,7 +34,8 @@ License
 void Foam::RBD::rigidBodyModel::applyRestraints
 (
     scalarField& tau,
-    Field<spatialVector>& fx
+    Field<spatialVector>& fx,
+    const rigidBodyModelState& state
 ) const
 {
     if (restraints_.empty())
@@ -47,7 +48,7 @@ void Foam::RBD::rigidBodyModel::applyRestraints
         DebugInfo << "Restraint " << restraints_[ri].name();
 
         // Accumulate the restraint forces
-        restraints_[ri].restrain(tau, fx);
+        restraints_[ri].restrain(tau, fx, state);
     }
 }
 

src/rigidBodyDynamics/rigidBodyModel/rigidBodyModel.H

@@ -340,7 +340,12 @@ public:
 
         //- Apply the restraints and accumulate the internal joint forces
         //  into the tau field and external forces into the fx field
-        void applyRestraints(scalarField& tau, Field<spatialVector>& fx) const;
+        void applyRestraints
+        (
+            scalarField& tau,
+            Field<spatialVector>& fx,
+            const rigidBodyModelState& state
+        ) const;
 
         //- Calculate the joint acceleration qDdot from the joint state q,
         //  velocity qDot, internal force tau (in the joint frame) and

src/rigidBodyDynamics/rigidBodySolvers/CrankNicolson/CrankNicolson.C

@@ -74,7 +74,7 @@ void Foam::RBD::rigidBodySolvers::CrankNicolson::solve
     // Accumulate the restraint forces
     scalarField rtau(tau);
     Field<spatialVector> rfx(fx);
-    model_.applyRestraints(rtau, rfx);
+    model_.applyRestraints(rtau, rfx, state());
 
     // Calculate the accelerations for the given state and forces
     model_.forwardDynamics(state(), rtau, rfx);

src/rigidBodyDynamics/rigidBodySolvers/Newmark/Newmark.C

@@ -81,7 +81,7 @@ void Foam::RBD::rigidBodySolvers::Newmark::solve
     // Accumulate the restraint forces
     scalarField rtau(tau);
     Field<spatialVector> rfx(fx);
-    model_.applyRestraints(rtau, rfx);
+    model_.applyRestraints(rtau, rfx, state());
 
     // Calculate the accelerations for the given state and forces
     model_.forwardDynamics(state(), rtau, rfx);

src/rigidBodyDynamics/rigidBodySolvers/symplectic/symplectic.C

@@ -83,7 +83,7 @@ void Foam::RBD::rigidBodySolvers::symplectic::solve
     // Accumulate the restraint forces
     scalarField rtau(tau);
     Field<spatialVector> rfx(fx);
-    model_.applyRestraints(rtau, rfx);
+    model_.applyRestraints(rtau, rfx, state());
 
     // Calculate the body acceleration for the given state
     // and restraint forces

src/thermophysicalModels/barotropicCompressibilityModel/Chung/Chung.C

@@ -54,6 +54,12 @@ Foam::compressibilityModels::Chung::Chung
 )
 :
     barotropicCompressibilityModel(compressibilityProperties, gamma, psiName),
+    pSat_
+    (
+        "pSat",
+        dimPressure,
+        compressibilityProperties_
+    ),
     psiv_
     (
         "psiv",
@@ -66,12 +72,6 @@ Foam::compressibilityModels::Chung::Chung
         dimCompressibility,
         compressibilityProperties_
     ),
-    rhovSat_
-    (
-        "rhovSat",
-        dimDensity,
-        compressibilityProperties_
-    ),
     rholSat_
     (
         "rholSat",
@@ -91,8 +91,8 @@ void Foam::compressibilityModels::Chung::correct()
     (
         sqrt
         (
-            (rhovSat_/psiv_)
-           /((scalar(1) - gamma_)*rhovSat_/psiv_ + gamma_*rholSat_/psil_)
+            pSat_
+           /((scalar(1) - gamma_)*pSat_ + gamma_*rholSat_/psil_)
         )
     );
 
@@ -111,9 +111,9 @@ bool Foam::compressibilityModels::Chung::read
 {
     barotropicCompressibilityModel::read(compressibilityProperties);
 
+    compressibilityProperties_.readEntry("pSat", pSat_);
     compressibilityProperties_.readEntry("psiv", psiv_);
     compressibilityProperties_.readEntry("psil", psil_);
-    compressibilityProperties_.readEntry("rhovSat", rhovSat_);
     compressibilityProperties_.readEntry("rholSat", rholSat_);
 
     return true;

src/thermophysicalModels/barotropicCompressibilityModel/Chung/Chung.H

@@ -57,10 +57,10 @@ class Chung
 {
     // Private data
 
+        dimensionedScalar pSat_;
         dimensionedScalar psiv_;
         dimensionedScalar psil_;
 
-        dimensionedScalar rhovSat_;
         dimensionedScalar rholSat_;
 
 

src/thermophysicalModels/barotropicCompressibilityModel/Wallis/Wallis.C

@@ -54,6 +54,12 @@ Foam::compressibilityModels::Wallis::Wallis
 )
 :
     barotropicCompressibilityModel(compressibilityProperties, gamma, psiName),
+    pSat_
+    (
+       "pSat",
+       dimPressure,
+       compressibilityProperties_
+    ),
     psiv_
     (
         "psiv",
@@ -66,12 +72,7 @@ Foam::compressibilityModels::Wallis::Wallis
         dimCompressibility,
         compressibilityProperties_
     ),
-    rhovSat_
-    (
-        "rhovSat",
-        dimDensity,
-        compressibilityProperties_
-    ),
+    rhovSat_(psiv_*pSat_),
     rholSat_
     (
         "rholSat",
@@ -89,7 +90,7 @@ void Foam::compressibilityModels::Wallis::correct()
 {
     psi_ =
         (gamma_*rhovSat_ + (scalar(1) - gamma_)*rholSat_)
-       *(gamma_*psiv_/rhovSat_ + (scalar(1) - gamma_)*psil_/rholSat_);
+       *(gamma_/pSat_ + (scalar(1) - gamma_)*psil_/rholSat_);
 }
 
 
@@ -100,9 +101,10 @@ bool Foam::compressibilityModels::Wallis::read
 {
     barotropicCompressibilityModel::read(compressibilityProperties);
 
+    compressibilityProperties_.readEntry("pSat", pSat_);
     compressibilityProperties_.readEntry("psiv", psiv_);
     compressibilityProperties_.readEntry("psil", psil_);
-    compressibilityProperties_.readEntry("rhovSat", rhovSat_);
+    rhovSat_ = psiv_*pSat_;
     compressibilityProperties_.readEntry("rholSat", rholSat_);
 
     return true;

src/thermophysicalModels/barotropicCompressibilityModel/Wallis/Wallis.H

@@ -57,6 +57,7 @@ class Wallis
 {
     // Private data
 
+        dimensionedScalar pSat_;
         dimensionedScalar psiv_;
         dimensionedScalar psil_;
 

src/thermophysicalModels/basic/derivedFvPatchFields/gradientEnergy/gradientEnergyFvPatchScalarField.C

@@ -119,6 +119,13 @@ void Foam::gradientEnergyFvPatchScalarField::updateCoeffs()
 }
 
 
+void Foam::gradientEnergyFvPatchScalarField::write(Ostream& os) const
+{
+    fixedGradientFvPatchScalarField::write(os);
+    os.writeEntry("value", *this);
+}
+
+
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 namespace Foam

src/thermophysicalModels/basic/derivedFvPatchFields/gradientEnergy/gradientEnergyFvPatchScalarField.H

@@ -158,6 +158,9 @@ public:
 
             //- Update the coefficients associated with the patch field
             virtual void updateCoeffs();
+
+        //- Write
+        virtual void write(Ostream&) const;
 };
 
 

src/thermophysicalModels/basic/rhoThermo/rhoThermos.C

@@ -33,6 +33,7 @@ License
 #include "incompressiblePerfectGas.H"
 #include "Boussinesq.H"
 #include "rhoConst.H"
+#include "rPolynomial.H"
 #include "perfectFluid.H"
 #include "PengRobinsonGas.H"
 #include "adiabaticPerfectFluid.H"
@@ -122,6 +123,18 @@ makeThermos
     specie
 );
 
+makeThermos
+(
+    rhoThermo,
+    heRhoThermo,
+    pureMixture,
+    constTransport,
+    sensibleEnthalpy,
+    hConstThermo,
+    rPolynomial,
+    specie
+);
+
 makeThermos
 (
     rhoThermo,
@@ -146,6 +159,18 @@ makeThermos
     specie
 );
 
+makeThermos
+(
+    rhoThermo,
+    heRhoThermo,
+    pureMixture,
+    polynomialTransport,
+    sensibleEnthalpy,
+    hPolynomialThermo,
+    rPolynomial,
+    specie
+);
+
 makeThermos
 (
     rhoThermo,
@@ -328,6 +353,18 @@ makeThermos
     specie
 );
 
+makeThermos
+(
+    rhoThermo,
+    heRhoThermo,
+    pureMixture,
+    constTransport,
+    sensibleInternalEnergy,
+    hConstThermo,
+    rPolynomial,
+    specie
+);
+
 makeThermos
 (
     rhoThermo,
@@ -340,6 +377,18 @@ makeThermos
     specie
 );
 
+makeThermos
+(
+    rhoThermo,
+    heRhoThermo,
+    pureMixture,
+    constTransport,
+    sensibleInternalEnergy,
+    eConstThermo,
+    rPolynomial,
+    specie
+);
+
 makeThermos
 (
     rhoThermo,

src/thermophysicalModels/reactionThermo/mixtures/SpecieMixture/SpecieMixture.C

@@ -60,6 +60,13 @@ Foam::scalar Foam::SpecieMixture<MixtureType>::W
 }
 
 
+template<class MixtureType>
+Foam::scalar Foam::SpecieMixture<MixtureType>::Hc(const label speciei) const
+{
+    return this->getLocalThermo(speciei).Hc();
+}
+
+
 template<class MixtureType>
 Foam::scalar Foam::SpecieMixture<MixtureType>::Cp
 (
@@ -84,6 +91,18 @@ Foam::scalar Foam::SpecieMixture<MixtureType>::Cv
 }
 
 
+template<class MixtureType>
+Foam::scalar Foam::SpecieMixture<MixtureType>::HE
+(
+    const label speciei,
+    const scalar p,
+    const scalar T
+) const
+{
+    return this->getLocalThermo(speciei).HE(p, T);
+}
+
+
 template<class MixtureType>
 Foam::scalar Foam::SpecieMixture<MixtureType>::Ha
 (
@@ -108,16 +127,6 @@ Foam::scalar Foam::SpecieMixture<MixtureType>::Hs
 }
 
 
-template<class MixtureType>
-Foam::scalar Foam::SpecieMixture<MixtureType>::Hc
-(
-    const label speciei
-) const
-{
-    return this->getLocalThermo(speciei).Hc();
-}
-
-
 template<class MixtureType>
 Foam::scalar Foam::SpecieMixture<MixtureType>::S
 (

src/thermophysicalModels/reactionThermo/mixtures/SpecieMixture/SpecieMixture.H

@@ -87,6 +87,9 @@ public:
             //- Molecular weight of the given specie [kg/kmol]
             virtual scalar W(const label speciei) const;
 
+            //- Chemical enthalpy [J/kg]
+            virtual scalar Hc(const label speciei) const;
+
 
         // Per specie thermo properties
 
@@ -106,6 +109,14 @@ public:
                 const scalar T
             ) const;
 
+            //- Enthalpy/Internal energy [J/kg]
+            virtual scalar HE
+            (
+                const label speciei,
+                const scalar p,
+                const scalar T
+            ) const;
+
             //- Absolute enthalpy [J/kg]
             virtual scalar Ha
             (
@@ -122,8 +133,6 @@ public:
                 const scalar T
             ) const;
 
-            //- Chemical enthalpy [J/kg]
-            virtual scalar Hc(const label speciei) const;
 
             //- Entropy [J/(kg K)]
             virtual scalar S

src/thermophysicalModels/reactionThermo/mixtures/basicSpecieMixture/basicSpecieMixture.H

@@ -89,6 +89,9 @@ public:
             //- Molecular weight of the given specie [kg/kmol]
             virtual scalar W(const label speciei) const = 0;
 
+            //- Chemical enthalpy [J/kg]
+            virtual scalar Hc(const label speciei) const = 0;
+
 
         // Per specie thermo properties
 
@@ -108,6 +111,14 @@ public:
                 const scalar T
             ) const = 0;
 
+            //- Enthalpy/Internal energy [J/kg]
+            virtual scalar HE
+            (
+                const label speciei,
+                const scalar p,
+                const scalar T
+            ) const = 0;
+
             //- Absolute enthalpy [J/kg]
             virtual scalar Ha
             (
@@ -124,8 +135,6 @@ public:
                 const scalar T
             ) const = 0;
 
-            //- Chemical enthalpy [J/kg]
-            virtual scalar Hc(const label speciei) const = 0;
 
             //- Entropy [J/(kg K)]
             virtual scalar S

src/thermophysicalModels/reactionThermo/rhoReactionThermo/rhoReactionThermos.C

@@ -39,6 +39,7 @@ License
 #include "sensibleEnthalpy.H"
 #include "thermo.H"
 #include "rhoConst.H"
+#include "rPolynomial.H"
 #include "perfectFluid.H"
 #include "adiabaticPerfectFluid.H"
 #include "Boussinesq.H"

src/thermophysicalModels/specie/equationOfState/rPolynomial/rPolynomial.C

@@ -0,0 +1,65 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2019 OpenFOAM Foundation
+-------------------------------------------------------------------------------
+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 "rPolynomial.H"
+#include "IOstreams.H"
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+template<class Specie>
+Foam::rPolynomial<Specie>::rPolynomial(const dictionary& dict)
+:
+    Specie(dict),
+    C_(dict.subDict("equationOfState").lookup("C"))
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+template<class Specie>
+void Foam::rPolynomial<Specie>::write(Ostream& os) const
+{
+    Specie::write(os);
+
+    dictionary dict("equationOfState");
+    dict.add("C", C_);
+
+    os  << indent << dict.dictName() << dict;
+}
+
+
+// * * * * * * * * * * * * * * * Ostream Operator  * * * * * * * * * * * * * //
+
+template<class Specie>
+Foam::Ostream& Foam::operator<<(Ostream& os, const rPolynomial<Specie>& pf)
+{
+    pf.write(os);
+    return os;
+}
+
+
+// ************************************************************************* //

src/thermophysicalModels/specie/equationOfState/rPolynomial/rPolynomial.H

@@ -0,0 +1,270 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2019 OpenFOAM Foundation
+-------------------------------------------------------------------------------
+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::rPolynomial
+
+Description
+    Reciprocal polynomial equation of state for liquids and solids
+
+    \f[
+        1/\rho = C_0 + C_1 T + C_2 T^2 - C_3 p - C_4 p T
+    \f]
+
+    This polynomial for the reciprocal of the density provides a much better fit
+    than the equivalent polynomial for the density and has the advantage that it
+    support coefficient mixing to support liquid and solid mixtures in an
+    efficient manner.
+
+Usage
+    \table
+        Property     | Description
+        C            | Density polynomial coefficients
+    \endtable
+
+    Example of the specification of the equation of state for pure water:
+    \verbatim
+    equationOfState
+    {
+        C (0.001278 -2.1055e-06 3.9689e-09 4.3772e-13 -2.0225e-16);
+    }
+    \endverbatim
+    Note: This fit is based on the small amount of data which is freely
+    available for the range 20-65degC and 1-100bar.
+
+SourceFiles
+    rPolynomialI.H
+    rPolynomial.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef rPolynomial_H
+#define rPolynomial_H
+
+#include "autoPtr.H"
+#include "VectorSpace.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// Forward declaration of friend functions and operators
+
+template<class Specie> class rPolynomial;
+
+template<class Specie>
+inline rPolynomial<Specie> operator+
+(
+    const rPolynomial<Specie>&,
+    const rPolynomial<Specie>&
+);
+
+template<class Specie>
+inline rPolynomial<Specie> operator*
+(
+    const scalar,
+    const rPolynomial<Specie>&
+);
+
+template<class Specie>
+inline rPolynomial<Specie> operator==
+(
+    const rPolynomial<Specie>&,
+    const rPolynomial<Specie>&
+);
+
+template<class Specie>
+Ostream& operator<<
+(
+    Ostream&,
+    const rPolynomial<Specie>&
+);
+
+
+/*---------------------------------------------------------------------------*\
+                           Class rPolynomial Declaration
+\*---------------------------------------------------------------------------*/
+
+template<class Specie>
+class rPolynomial
+:
+    public Specie
+{
+    // Private Data
+
+        class coeffList
+        :
+            public VectorSpace<coeffList, scalar, 5>
+        {
+            public:
+
+            // Constructors
+
+                //- Construct null
+                inline coeffList()
+                {}
+
+                //- Construct from Istream
+                inline coeffList(Istream& is)
+                :
+                    VectorSpace<coeffList, scalar, 5>(is)
+                {}
+        };
+
+
+        //- Density coefficients
+        coeffList C_;
+
+
+public:
+
+    // Constructors
+
+        //- Construct from components
+        inline rPolynomial
+        (
+            const Specie& sp,
+            const coeffList& coeffs
+        );
+
+        //- Construct from dictionary
+        rPolynomial(const dictionary& dict);
+
+        //- Construct as named copy
+        inline rPolynomial(const word& name, const rPolynomial&);
+
+        //- Construct and return a clone
+        inline autoPtr<rPolynomial> clone() const;
+
+        // Selector from dictionary
+        inline static autoPtr<rPolynomial> New(const dictionary& dict);
+
+
+    // Member Functions
+
+        //- Return the instantiated type name
+        static word typeName()
+        {
+            return "rPolynomial<" + word(Specie::typeName_()) + '>';
+        }
+
+
+        // Fundamental properties
+
+            //- Is the equation of state is incompressible i.e. rho != f(p)
+            static const bool incompressible = false;
+
+            //- Is the equation of state is isochoric i.e. rho = const
+            static const bool isochoric = false;
+
+            //- Return density [kg/m^3]
+            inline scalar rho(scalar p, scalar T) const;
+
+            //- Return enthalpy departure [J/kg]
+            inline scalar H(const scalar p, const scalar T) const;
+
+            //- Return Cp departure [J/(kg K]
+            inline scalar Cp(scalar p, scalar T) const;
+
+            //- Return internal energy departure [J/kg]
+            inline scalar E(const scalar p, const scalar T) const;
+
+            //- Return Cv departure [J/(kg K]
+            inline scalar Cv(scalar p, scalar T) const;
+
+            //- Return entropy [J/kg/K]
+            inline scalar S(const scalar p, const scalar T) const;
+
+            //- Return compressibility [s^2/m^2]
+            inline scalar psi(scalar p, scalar T) const;
+
+            //- Return compression factor []
+            inline scalar Z(scalar p, scalar T) const;
+
+            //- Return (Cp - Cv) [J/(kg K]
+            inline scalar CpMCv(scalar p, scalar T) const;
+
+
+        // IO
+
+            //- Write to Ostream
+            void write(Ostream& os) const;
+
+
+    // Member Operators
+
+        inline void operator+=(const rPolynomial&);
+        inline void operator*=(const scalar);
+
+
+    // Friend operators
+
+        friend rPolynomial operator+ <Specie>
+        (
+            const rPolynomial&,
+            const rPolynomial&
+        );
+
+        friend rPolynomial operator* <Specie>
+        (
+            const scalar s,
+            const rPolynomial&
+        );
+
+        friend rPolynomial operator== <Specie>
+        (
+            const rPolynomial&,
+            const rPolynomial&
+        );
+
+
+    // Ostream Operator
+
+        friend Ostream& operator<< <Specie>
+        (
+            Ostream&,
+            const rPolynomial&
+        );
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#include "rPolynomialI.H"
+
+#ifdef NoRepository
+    #include "rPolynomial.C"
+#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/thermophysicalModels/specie/equationOfState/rPolynomial/rPolynomialI.H

@@ -0,0 +1,235 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2019 OpenFOAM Foundation
+-------------------------------------------------------------------------------
+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 "rPolynomial.H"
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+template<class Specie>
+inline Foam::rPolynomial<Specie>::rPolynomial
+(
+    const Specie& sp,
+    const coeffList& coeffs
+)
+:
+    Specie(sp),
+    C_(coeffs)
+{}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+template<class Specie>
+inline Foam::rPolynomial<Specie>::rPolynomial
+(
+    const word& name,
+    const rPolynomial<Specie>& rp
+)
+:
+    Specie(name, rp),
+    C_(rp.C_)
+{}
+
+
+template<class Specie>
+inline Foam::autoPtr<Foam::rPolynomial<Specie>>
+Foam::rPolynomial<Specie>::clone() const
+{
+    return autoPtr<rPolynomial<Specie>>::New(*this);
+}
+
+
+template<class Specie>
+inline Foam::autoPtr<Foam::rPolynomial<Specie>>
+Foam::rPolynomial<Specie>::New
+(
+    const dictionary& dict
+)
+{
+    return autoPtr<rPolynomial<Specie>>::New(dict);
+}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+template<class Specie>
+inline Foam::scalar Foam::rPolynomial<Specie>::rho(scalar p, scalar T) const
+{
+    return 1/(C_[0] + (C_[1] + C_[2]*T - C_[4]*p)*T - C_[3]*p);
+}
+
+
+template<class Specie>
+inline Foam::scalar Foam::rPolynomial<Specie>::H(scalar p, scalar T) const
+{
+    return 0;
+}
+
+
+template<class Specie>
+inline Foam::scalar Foam::rPolynomial<Specie>::Cp(scalar p, scalar T) const
+{
+    return 0;
+}
+
+
+template<class Specie>
+inline Foam::scalar Foam::rPolynomial<Specie>::E(scalar p, scalar T) const
+{
+    return 0;
+}
+
+
+template<class Specie>
+inline Foam::scalar Foam::rPolynomial<Specie>::Cv(scalar p, scalar T) const
+{
+    return 0;
+}
+
+
+template<class Specie>
+inline Foam::scalar Foam::rPolynomial<Specie>::S(scalar p, scalar T) const
+{
+    return 0;
+}
+
+
+template<class Specie>
+inline Foam::scalar Foam::rPolynomial<Specie>::psi(scalar p, scalar T) const
+{
+    return sqr(rho(p, T))*(C_[3] + C_[4]*T);
+}
+
+
+template<class Specie>
+inline Foam::scalar Foam::rPolynomial<Specie>::Z(scalar p, scalar T) const
+{
+    return 1;
+}
+
+
+template<class Specie>
+inline Foam::scalar Foam::rPolynomial<Specie>::CpMCv(scalar p, scalar T) const
+{
+    return 0;
+}
+
+
+// * * * * * * * * * * * * * * * Member Operators  * * * * * * * * * * * * * //
+
+template<class Specie>
+inline void Foam::rPolynomial<Specie>::operator+=
+(
+    const rPolynomial<Specie>& rp
+)
+{
+    const scalar Y1 = this->Y();
+    Specie::operator+=(rp);
+
+    if (mag(this->Y()) > SMALL)
+    {
+        C_ = (Y1*C_ + rp.Y()*rp.C_)/this->Y();
+    }
+}
+
+
+template<class Specie>
+inline void Foam::rPolynomial<Specie>::operator*=(const scalar s)
+{
+    Specie::operator*=(s);
+}
+
+
+// * * * * * * * * * * * * * * * Friend Operators  * * * * * * * * * * * * * //
+
+template<class Specie>
+inline Foam::rPolynomial<Specie> Foam::operator+
+(
+    const rPolynomial<Specie>& rp1,
+    const rPolynomial<Specie>& rp2
+)
+{
+    Specie sp
+    (
+        static_cast<const Specie&>(rp1)
+      + static_cast<const Specie&>(rp2)
+    );
+
+    if (mag(sp.Y()) < SMALL)
+    {
+        return rPolynomial<Specie>
+        (
+            sp,
+            rp1.C_
+        );
+    }
+
+    return rPolynomial<Specie>
+    (
+        sp,
+        (rp1.Y()*rp1.C_ + rp2.Y()*rp2.C_)/sp.Y()
+    );
+}
+
+
+template<class Specie>
+inline Foam::rPolynomial<Specie> Foam::operator*
+(
+    const scalar s,
+    const rPolynomial<Specie>& rp
+)
+{
+    return rPolynomial<Specie>
+    (
+        s*static_cast<const Specie&>(rp),
+        rp.C_
+    );
+}
+
+
+template<class Specie>
+inline Foam::rPolynomial<Specie> Foam::operator==
+(
+    const rPolynomial<Specie>& rp1,
+    const rPolynomial<Specie>& rp2
+)
+{
+    Specie sp
+    (
+        static_cast<const Specie&>(rp1)
+     == static_cast<const Specie&>(rp2)
+    );
+
+    return rPolynomial<Specie>
+    (
+        sp,
+        (rp1.Y()*rp1.C_ - rp2.Y()*rp2.C_)/sp.Y()
+    );
+}
+
+
+// ************************************************************************* //

src/thermophysicalModels/specie/include/thermoPhysicsTypes.H

@@ -38,6 +38,7 @@ Description
 #include "specie.H"
 #include "perfectGas.H"
 #include "incompressiblePerfectGas.H"
+#include "rPolynomial.H"
 #include "perfectFluid.H"
 #include "adiabaticPerfectFluid.H"
 #include "rhoConst.H"
@@ -155,6 +156,20 @@ namespace Foam
     >
     constFluidHThermoPhysics;
 
+    typedef
+    constTransport
+    <
+        species::thermo
+        <
+            hConstThermo
+            <
+                rPolynomial<specie>
+            >,
+            sensibleEnthalpy
+        >
+    >
+    constrPolFluidHThermoPhysics;
+
     typedef
     constTransport
     <
@@ -280,6 +295,20 @@ namespace Foam
     >
     constFluidEThermoPhysics;
 
+    typedef
+    constTransport
+    <
+        species::thermo
+        <
+            eConstThermo
+            <
+                perfectFluid<specie>
+            >,
+            sensibleInternalEnergy
+        >
+    >
+    constrPolFluidEThermoPhysics;
+
     typedef
     constTransport
     <