Updates and corrections for the polynomial thermo

Added missing operators, and converted polys to be molar-intensive
to enable mixing using existing molar-weighted strategy

src/thermophysicalModels/specie/equationOfState/icoPolynomial/icoPolynomial.C

@@ -39,7 +39,9 @@ icoPolynomial<PolySize>::icoPolynomial(Istream& is)
 :
     specie(is),
     rhoPolynomial_("rhoPolynomial", is)
-{}
+{
+    rhoPolynomial_ *= this->W();
+}
 
 
 // * * * * * * * * * * * * * * * Ostream Operator  * * * * * * * * * * * * * //
@@ -47,7 +49,8 @@ icoPolynomial<PolySize>::icoPolynomial(Istream& is)
 template<int PolySize>
 Ostream& operator<<(Ostream& os, const icoPolynomial<PolySize>& ip)
 {
-    os  << static_cast<const specie&>(ip);
+    os  << static_cast<const specie&>(ip) << tab
+        << "rhoPolynomial" << tab << ip.rhoPolynomial_/ip.W();
 
     os.check
     (

src/thermophysicalModels/specie/equationOfState/icoPolynomial/icoPolynomial.H

@@ -99,7 +99,8 @@ class icoPolynomial
 {
     // Private data
 
-        //- Density [kg/m^3]
+        //- Density
+        //  Note: input in [kg/m3], but internally uses [kg/m3/kmol]
         Polynomial<PolySize> rhoPolynomial_;
 
 
@@ -117,6 +118,9 @@ public:
         //- Construct from Istream
         icoPolynomial(Istream&);
 
+        //- Construct as copy
+        inline icoPolynomial(const icoPolynomial&);
+
         //- Construct as named copy
         inline icoPolynomial(const word& name, const icoPolynomial&);
 
@@ -141,6 +145,7 @@ public:
 
     // Member operators
 
+        inline icoPolynomial& operator=(const icoPolynomial&);
         inline void operator+=(const icoPolynomial&);
         inline void operator-=(const icoPolynomial&);
 

src/thermophysicalModels/specie/equationOfState/icoPolynomial/icoPolynomialI.H

@@ -42,6 +42,17 @@ inline Foam::icoPolynomial<PolySize>::icoPolynomial
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
+template<int PolySize>
+inline Foam::icoPolynomial<PolySize>::icoPolynomial
+(
+    const icoPolynomial<PolySize>& ip
+)
+:
+    specie(ip),
+    rhoPolynomial_(ip.rhoPolynomial_)
+{}
+
+
 template<int PolySize>
 inline Foam::icoPolynomial<PolySize>::icoPolynomial
 (
@@ -78,7 +89,7 @@ Foam::icoPolynomial<PolySize>::New(Istream& is)
 template<int PolySize>
 inline Foam::scalar Foam::icoPolynomial<PolySize>::rho(scalar, scalar T) const
 {
-    return rhoPolynomial_.evaluate(T);
+    return rhoPolynomial_.evaluate(T)/this->W();
 }
 
 
@@ -98,6 +109,20 @@ inline Foam::scalar Foam::icoPolynomial<PolySize>::Z(scalar, scalar) const
 
 // * * * * * * * * * * * * * * * Member Operators  * * * * * * * * * * * * * //
 
+template<int PolySize>
+inline Foam::icoPolynomial<PolySize>& Foam::icoPolynomial<PolySize>::operator=
+(
+    const icoPolynomial<PolySize>& ip
+)
+{
+    specie::operator=(ip);
+
+    rhoPolynomial_ = ip.rhoPolynomial_;
+
+    return *this;
+}
+
+
 template<int PolySize>
 inline void Foam::icoPolynomial<PolySize>::operator+=
 (
@@ -122,6 +147,7 @@ inline void Foam::icoPolynomial<PolySize>::operator-=
 )
 {
     scalar molr1 = this->nMoles();
+
     specie::operator-=(ip);
 
     molr1 /= this->nMoles();
@@ -147,15 +173,15 @@ Foam::icoPolynomial<PolySize> Foam::operator+
     const icoPolynomial<PolySize>& ip2
 )
 {
-    scalar mol1 = ip1.nMoles();
+    scalar nMoles = ip1.nMoles() + ip2.nMoles();
-    scalar mol2 = ip2.nMoles();
+    scalar molr1 = ip1.nMoles()/nMoles;
-    scalar nMoles = mol1 + mol2;
+    scalar molr2 = ip2.nMoles()/nMoles;
 
     return icoPolynomial<PolySize>
     (
         static_cast<const specie&>(ip1)
       + static_cast<const specie&>(ip2),
-        (mol1/nMoles)*ip1.rhoPolynomial_ + (mol2/nMoles)*ip2.rhoPolynomial_
+        molr1*ip1.rhoPolynomial_ + molr2*ip2.rhoPolynomial_
     );
 }
 
@@ -167,15 +193,15 @@ Foam::icoPolynomial<PolySize> Foam::operator-
     const icoPolynomial<PolySize>& ip2
 )
 {
-    scalar mol1 = ip1.nMoles();
+    scalar nMoles = ip1.nMoles() + ip2.nMoles();
-    scalar mol2 = ip2.nMoles();
+    scalar molr1 = ip1.nMoles()/nMoles;
-    scalar nMoles = mol1 + mol2;
+    scalar molr2 = ip2.nMoles()/nMoles;
 
     return icoPolynomial<PolySize>
     (
         static_cast<const specie&>(ip1)
       - static_cast<const specie&>(ip2),
-        (mol1/nMoles)*ip1.rhoPolynomial_ - (mol2/nMoles)*ip2.rhoPolynomial_
+        molr1*ip1.rhoPolynomial_ - molr2*ip2.rhoPolynomial_
     );
 }
 

src/thermophysicalModels/specie/thermo/hPolynomial/hPolynomialThermo.C

@@ -39,11 +39,21 @@ Foam::hPolynomialThermo<EquationOfState, PolySize>::hPolynomialThermo
     Hf_(readScalar(is)),
     Sf_(readScalar(is)),
     cpPolynomial_("cpPolynomial", is),
-    dhPolynomial_(cpPolynomial_.integrate()),
+    hPolynomial_(),
-    dsPolynomial_(cpPolynomial_.integrateMinus1())
+    sPolynomial_()
 {
-    // Offset dh poly so that it is relative to the enthalpy at Tstd
+    Hf_ *= this->W();
-    dhPolynomial_[0] -= dhPolynomial_.evaluate(specie::Tstd);
+    Sf_ *= this->W();
+    cpPolynomial_ *= this->W();
+
+    hPolynomial_ = cpPolynomial_.integrate();
+    sPolynomial_ = cpPolynomial_.integrateMinus1();
+
+    // Offset h poly so that it is relative to the enthalpy at Tstd
+    hPolynomial_[0] += Hf_ - hPolynomial_.evaluate(specie::Tstd);
+
+    // Offset s poly so that it is relative to the entropy at Tstd
+    sPolynomial_[0] += Sf_ - sPolynomial_.evaluate(specie::Tstd);
 }
 
 
@@ -57,15 +67,17 @@ Foam::Ostream& Foam::operator<<
 )
 {
     os  << static_cast<const EquationOfState&>(pt) << tab
-        << pt.Hf_ << tab
+        << pt.Hf_/pt.W() << tab
         << pt.Sf_ << tab
-        << pt.cpPolynomial_ << tab
+        << "cpPolynomial" << tab << pt.cpPolynomial_/pt.W();
-        << pt.dhPolynomial_ << tab
-        << pt.dsPolynomial;
 
     os.check
     (
-        "operator<<(Ostream& os, const hPolynomialThermo<EquationOfState>& pt)"
+        "operator<<"
+        "("
+            "Ostream&, "
+            "const hPolynomialThermo<EquationOfState, PolySize>&"
+        ")"
     );
 
     return os;

src/thermophysicalModels/specie/thermo/hPolynomial/hPolynomialThermo.H

@@ -100,20 +100,22 @@ class hPolynomialThermo
 {
     // Private data
 
-        //- Heat of formation [J/kg]
+        //- Heat of formation
+        //  Note: input in [J/kg], but internally uses [J/kmol]
         scalar Hf_;
 
-        //- Standard entropy [J/(kg.K)]
+        //- Standard entropy
+        //  Note: input in [J/kg/K], but internally uses [J/kmol/K]
         scalar Sf_;
 
         //- Specific heat at constant pressure [J/(kg.K)]
         Polynomial<PolySize> cpPolynomial_;
 
         //- Enthalpy - derived from cp [J/kg] - relative to Tstd
-        typename Polynomial<PolySize>::intPolyType dhPolynomial_;
+        typename Polynomial<PolySize>::intPolyType hPolynomial_;
 
-        //- Entropy - derived from cp [J/(kg.K)]
+        //- Entropy - derived from cp [J/(kg.K)] - relative to Tstd
-        Polynomial<PolySize> dsPolynomial_;
+        Polynomial<PolySize> sPolynomial_;
 
 
     // Private member functions
@@ -125,8 +127,8 @@ class hPolynomialThermo
             const scalar Hf,
             const scalar Sf,
             const Polynomial<PolySize>& cpPoly,
-            const typename Polynomial<PolySize>::intPolyType& dhPoly,
+            const typename Polynomial<PolySize>::intPolyType& hPoly,
-            const Polynomial<PolySize>& dsPoly
+            const Polynomial<PolySize>& sPoly
         );
 
 
@@ -137,6 +139,9 @@ public:
         //- Construct from dictionary
         hPolynomialThermo(Istream& is);
 
+        //- Construct as copy
+        inline hPolynomialThermo(const hPolynomialThermo&);
+
         //- Construct as a named copy
         inline hPolynomialThermo(const word&, const hPolynomialThermo&);
 
@@ -161,8 +166,10 @@ public:
 
     // Member operators
 
+        inline hPolynomialThermo& operator=(const hPolynomialThermo&);
         inline void operator+=(const hPolynomialThermo&);
         inline void operator-=(const hPolynomialThermo&);
+        inline void operator*=(const scalar);
 
 
     // Friend operators

src/thermophysicalModels/specie/thermo/hPolynomial/hPolynomialThermoI.H

@@ -35,21 +35,36 @@ inline Foam::hPolynomialThermo<EquationOfState, PolySize>::hPolynomialThermo
     const scalar Hf,
     const scalar Sf,
     const Polynomial<PolySize>& cpPoly,
-    const typename Polynomial<PolySize>::intPolyType& dhPoly,
+    const typename Polynomial<PolySize>::intPolyType& hPoly,
-    const Polynomial<PolySize>& dsPoly
+    const Polynomial<PolySize>& sPoly
 )
 :
     EquationOfState(pt),
     Hf_(Hf),
     Sf_(Sf),
     cpPolynomial_(cpPoly),
-    dhPolynomial_(dhPoly),
+    hPolynomial_(hPoly),
-    dsPolynomial_(dsPoly)
+    sPolynomial_(sPoly)
 {}
 
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
+template<class EquationOfState, int PolySize>
+inline Foam::hPolynomialThermo<EquationOfState, PolySize>::hPolynomialThermo
+(
+    const hPolynomialThermo& pt
+)
+:
+    EquationOfState(pt),
+    Hf_(pt.Hf_),
+    Sf_(pt.Sf_),
+    cpPolynomial_(pt.cpPolynomial_),
+    hPolynomial_(pt.hPolynomial_),
+    sPolynomial_(pt.sPolynomial_)
+{}
+
+
 template<class EquationOfState, int PolySize>
 inline Foam::hPolynomialThermo<EquationOfState, PolySize>::hPolynomialThermo
 (
@@ -61,8 +76,8 @@ inline Foam::hPolynomialThermo<EquationOfState, PolySize>::hPolynomialThermo
     Hf_(pt.Hf_),
     Sf_(pt.Sf_),
     cpPolynomial_(pt.cpPolynomial_),
-    dhPolynomial_(pt.dhPolynomial_),
+    hPolynomial_(pt.hPolynomial_),
-    dsPolynomial_(pt.dsPolynomial_)
+    sPolynomial_(pt.sPolynomial_)
 {}
 
 
@@ -74,7 +89,7 @@ inline Foam::scalar Foam::hPolynomialThermo<EquationOfState, PolySize>::cp
     const scalar T
 ) const
 {
-    return cpPolynomial_.evaluate(T)*this->W();
+    return cpPolynomial_.evaluate(T);
 }
 
 
@@ -84,7 +99,7 @@ inline Foam::scalar Foam::hPolynomialThermo<EquationOfState, PolySize>::h
     const scalar T
 ) const
 {
-    return (dhPolynomial_.evaluate(T) + Hf_)*this->W();
+    return hPolynomial_.evaluate(T);
 }
 
 
@@ -94,7 +109,7 @@ inline Foam::scalar Foam::hPolynomialThermo<EquationOfState, PolySize>::hs
     const scalar T
 ) const
 {
-    return dhPolynomial_.evaluate(T)*this->W();
+    return h(T) - hc();
 }
 
 
@@ -102,7 +117,7 @@ template<class EquationOfState, int PolySize>
 inline Foam::scalar Foam::hPolynomialThermo<EquationOfState, PolySize>::hc()
 const
 {
-    return Hf_*this->W();
+    return Hf_;
 }
 
 
@@ -112,12 +127,31 @@ inline Foam::scalar Foam::hPolynomialThermo<EquationOfState, PolySize>::s
     const scalar T
 ) const
 {
-    return (dsPolynomial_.evaluate(T) + Sf_)*this->W();
+    return sPolynomial_.evaluate(T);
 }
 
 
 // * * * * * * * * * * * * * * * Member Operators  * * * * * * * * * * * * * //
 
+template<class EquationOfState, int PolySize>
+inline Foam::hPolynomialThermo<EquationOfState, PolySize>&
+Foam::hPolynomialThermo<EquationOfState, PolySize>::operator=
+(
+    const hPolynomialThermo<EquationOfState, PolySize>& pt
+)
+{
+    EquationOfState::operator=(pt);
+
+    Hf_ = pt.Hf_;
+    Sf_ = pt.Sf_;
+    cpPolynomial_ = pt.cpPolynomial_;
+    hPolynomial_ = pt.hPolynomial_;
+    sPolynomial_ = pt.sPolynomial_;
+
+    return *this;
+}
+
+
 template<class EquationOfState, int PolySize>
 inline void Foam::hPolynomialThermo<EquationOfState, PolySize>::operator+=
 (
@@ -134,8 +168,8 @@ inline void Foam::hPolynomialThermo<EquationOfState, PolySize>::operator+=
     Hf_ = molr1*Hf_ + molr2*pt.Hf_;
     Sf_ = molr1*Sf_ + molr2*pt.Sf_;
     cpPolynomial_ = molr1*cpPolynomial_ + molr2*pt.cpPolynomial_;
-    dhPolynomial_ = molr1*dhPolynomial_ + molr2*pt.dhPolynomial_;
+    hPolynomial_ = molr1*hPolynomial_ + molr2*pt.hPolynomial_;
-    dsPolynomial_ = molr1*dsPolynomial_ + molr2*pt.dsPolynomial_;
+    sPolynomial_ = molr1*sPolynomial_ + molr2*pt.sPolynomial_;
 }
 
 
@@ -155,8 +189,18 @@ inline void Foam::hPolynomialThermo<EquationOfState, PolySize>::operator-=
     Hf_ = molr1*Hf_ - molr2*pt.Hf_;
     Sf_ = molr1*Sf_ - molr2*pt.Sf_;
     cpPolynomial_ = molr1*cpPolynomial_ - molr2*pt.cpPolynomial_;
-    dhPolynomial_ = molr1*dhPolynomial_ - molr2*pt.dhPolynomial_;
+    hPolynomial_ = molr1*hPolynomial_ - molr2*pt.hPolynomial_;
-    dsPolynomial_ = molr1*dsPolynomial_ - molr2*pt.dsPolynomial_;
+    sPolynomial_ = molr1*sPolynomial_ - molr2*pt.sPolynomial_;
+}
+
+
+template<class EquationOfState, int PolySize>
+inline void Foam::hPolynomialThermo<EquationOfState, PolySize>::operator*=
+(
+    const scalar s
+)
+{
+    EquationOfState::operator*=(s);
 }
 
 
@@ -169,22 +213,20 @@ inline Foam::hPolynomialThermo<EquationOfState, PolySize> Foam::operator+
     const hPolynomialThermo<EquationOfState, PolySize>& pt2
 )
 {
-    EquationOfState eofs
+    EquationOfState eofs = pt1;
-    (
+    eofs += pt2;
-        static_cast<const EquationOfState&>(pt1)
-      + static_cast<const EquationOfState&>(pt2)
-    );
 
     scalar molr1 = pt1.nMoles()/eofs.nMoles();
     scalar molr2 = pt2.nMoles()/eofs.nMoles();
+
     return hPolynomialThermo<EquationOfState, PolySize>
     (
         eofs,
         molr1*pt1.Hf_ + molr2*pt2.Hf_,
         molr1*pt1.Sf_ + molr2*pt2.Sf_,
         molr1*pt1.cpPolynomial_ + molr2*pt2.cpPolynomial_,
-        molr1*pt1.dhPolynomial_ + molr2*pt2.dhPolynomial_,
+        molr1*pt1.hPolynomial_ + molr2*pt2.hPolynomial_,
-        molr1*pt1.dsPolynomial_ + molr2*pt2.dsPolynomial_
+        molr1*pt1.sPolynomial_ + molr2*pt2.sPolynomial_
     );
 }
 
@@ -196,22 +238,20 @@ inline Foam::hPolynomialThermo<EquationOfState, PolySize> Foam::operator-
     const hPolynomialThermo<EquationOfState, PolySize>& pt2
 )
 {
-    EquationOfState eofs
+    EquationOfState eofs = pt1;
-    (
+    eofs -= pt2;
-        static_cast<const EquationOfState&>(pt1)
-      - static_cast<const EquationOfState&>(pt2)
-    );
 
     scalar molr1 = pt1.nMoles()/eofs.nMoles();
     scalar molr2 = pt2.nMoles()/eofs.nMoles();
+
     return hPolynomialThermo<EquationOfState, PolySize>
     (
         eofs,
         molr1*pt1.Hf_ - molr2*pt2.Hf_,
         molr1*pt1.Sf_ - molr2*pt2.Sf_,
         molr1*pt1.cpPolynomial_ - molr2*pt2.cpPolynomial_,
-        molr1*pt1.dhPolynomial_ - molr2*pt2.dhPolynomial_,
+        molr1*pt1.hPolynomial_ - molr2*pt2.hPolynomial_,
-        molr1*pt1.dsPolynomial_ - molr2*pt2.dsPolynomial_
+        molr1*pt1.sPolynomial_ - molr2*pt2.sPolynomial_
     );
 }
 
@@ -229,8 +269,8 @@ inline Foam::hPolynomialThermo<EquationOfState, PolySize> Foam::operator*
         pt.Hf_,
         pt.Sf_,
         pt.cpPolynomial_,
-        pt.dhPolynomial_,
+        pt.hPolynomial_,
-        pt.dsPolynomial_
+        pt.sPolynomial_
     );
 }
 

src/thermophysicalModels/specie/transport/polynomial/polynomialTransport.C

@@ -35,7 +35,10 @@ Foam::polynomialTransport<Thermo, PolySize>::polynomialTransport(Istream& is)
     Thermo(is),
     muPolynomial_("muPolynomial", is),
     kappaPolynomial_("kappaPolynomial", is)
-{}
+{
+    muPolynomial_ *= this->W();
+    kappaPolynomial_ *= this->W();
+}
 
 
 // * * * * * * * * * * * * * * * IOstream Operators  * * * * * * * * * * * * //
@@ -47,7 +50,9 @@ Foam::Ostream& Foam::operator<<
     const polynomialTransport<Thermo, PolySize>& pt
 )
 {
-    os << static_cast<const Thermo&>(pt);
+    os  << static_cast<const Thermo&>(pt) << tab
+        << "muPolynomial" << tab << pt.muPolynomial_/pt.W() << tab
+        << "kappaPolynomial" << tab << pt.kappaPolynomial_/pt.W();
 
     os.check
     (

src/thermophysicalModels/specie/transport/polynomial/polynomialTransport.H

@@ -96,9 +96,11 @@ class polynomialTransport
     // Private data
 
         //- Dynamic viscosity
+        //  Note: input in [Pa.s], but internally uses [Pa.s/kmol]
         Polynomial<PolySize> muPolynomial_;
 
         //- Thermal conductivity
+        //  Note: input in [W/m/K], but internally uses [W/m/K/kmol]
         Polynomial<PolySize> kappaPolynomial_;
 
 
@@ -117,6 +119,9 @@ public:
 
     // Constructors
 
+        //- Construct copy
+        inline polynomialTransport(const polynomialTransport&);
+
         //- Construct as named copy
         inline polynomialTransport(const word&, const polynomialTransport&);
 
@@ -148,6 +153,9 @@ public:
     // Member operators
 
         inline polynomialTransport& operator=(const polynomialTransport&);
+        inline void operator+=(const polynomialTransport&);
+        inline void operator-=(const polynomialTransport&);
+        inline void operator*=(const scalar);
 
 
     // Friend operators

src/thermophysicalModels/specie/transport/polynomial/polynomialTransportI.H

@@ -28,6 +28,18 @@ License
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
+template<class Thermo, int PolySize>
+inline Foam::polynomialTransport<Thermo, PolySize>::polynomialTransport
+(
+    const polynomialTransport& pt
+)
+:
+    Thermo(pt),
+    muPolynomial_(pt.muPolynomial_),
+    kappaPolynomial_(pt.kappaPolynomial_)
+{}
+
+
 template<class Thermo, int PolySize>
 inline Foam::polynomialTransport<Thermo, PolySize>::polynomialTransport
 (
@@ -85,7 +97,7 @@ inline Foam::scalar Foam::polynomialTransport<Thermo, PolySize>::mu
     const scalar T
 ) const
 {
-    return muPolynomial_.evaluate(T);
+    return muPolynomial_.evaluate(T)/this->W();
 }
 
 
@@ -95,7 +107,7 @@ inline Foam::scalar Foam::polynomialTransport<Thermo, PolySize>::kappa
     const scalar T
 ) const
 {
-    return kappaPolynomial_.evaluate(T);
+    return kappaPolynomial_.evaluate(T)/this->W();
 }
 
 
@@ -132,6 +144,52 @@ Foam::polynomialTransport<Thermo, PolySize>::operator=
 }
 
 
+template<class Thermo, int PolySize>
+inline void Foam::polynomialTransport<Thermo, PolySize>::operator+=
+(
+    const polynomialTransport<Thermo, PolySize>& pt
+)
+{
+    scalar molr1 = this->nMoles();
+
+    Thermo::operator+=(pt);
+
+    molr1 /= this->nMoles();
+    scalar molr2 = pt.nMoles()/this->nMoles();
+
+    muPolynomial_ = molr1*muPolynomial_ + molr2*pt.muPolynomial_;
+    kappaPolynomial_ = molr1*kappaPolynomial_ + molr2*pt.kappaPolynomial_;
+}
+
+
+template<class Thermo, int PolySize>
+inline void Foam::polynomialTransport<Thermo, PolySize>::operator-=
+(
+    const polynomialTransport<Thermo, PolySize>& pt
+)
+{
+    scalar molr1 = this->nMoles();
+
+    Thermo::operator-=(pt);
+
+    molr1 /= this->nMoles();
+    scalar molr2 = pt.nMoles()/this->nMoles();
+
+    muPolynomial_ = molr1*muPolynomial_ - molr2*pt.muPolynomial_;
+    kappaPolynomial_ = molr1*kappaPolynomial_ - molr2*pt.kappaPolynomial_;
+}
+
+
+template<class Thermo, int PolySize>
+inline void Foam::polynomialTransport<Thermo, PolySize>::operator*=
+(
+    const scalar s
+)
+{
+    Thermo::operator*=(s);
+}
+
+
 // * * * * * * * * * * * * * * * Friend Operators  * * * * * * * * * * * * * //
 
 template<class Thermo, int PolySize>
@@ -146,7 +204,6 @@ inline Foam::polynomialTransport<Thermo, PolySize> Foam::operator+
         static_cast<const Thermo&>(pt1) + static_cast<const Thermo&>(pt2)
     );
 
-
     scalar molr1 = pt1.nMoles()/t.nMoles();
     scalar molr2 = pt2.nMoles()/t.nMoles();