PengRobinsonGas: Updated mixing-rule for critical pressure

Resolves bug-report http://www.openfoam.org/mantisbt/view.php?id=1329

src/thermophysicalModels/specie/equationOfState/PengRobinsonGas/PengRobinsonGas.C

@@ -33,6 +33,8 @@ Foam::PengRobinsonGas<Specie>::PengRobinsonGas(Istream& is)
 :
     Specie(is),
     Tc_(readScalar(is)),
+    Vc_(readScalar(is)),
+    Zc_(readScalar(is)),
     Pc_(readScalar(is)),
     omega_(readScalar(is))
 {
@@ -48,9 +50,13 @@ Foam::PengRobinsonGas<Specie>::PengRobinsonGas
 :
     Specie(dict),
     Tc_(readScalar(dict.subDict("equationOfState").lookup("Tc"))),
+    Vc_(readScalar(dict.subDict("equationOfState").lookup("Vc"))),
+    Zc_(1.0),
     Pc_(readScalar(dict.subDict("equationOfState").lookup("Pc"))),
     omega_(readScalar(dict.subDict("equationOfState").lookup("omega")))
-{}
+{
+    Zc_ = Pc_*Vc_/(specie::RR*Tc_);
+}
 
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
@@ -74,6 +80,8 @@ Foam::Ostream& Foam::operator<<
 {
     os  << static_cast<const Specie&>(pg)
         << token::SPACE << pg.Tc_
+        << token::SPACE << pg.Vc_
+        << token::SPACE << pg.Zc_
         << token::SPACE << pg.Pc_
         << token::SPACE << pg.omega_;
 

src/thermophysicalModels/specie/equationOfState/PengRobinsonGas/PengRobinsonGas.H

@@ -98,10 +98,16 @@ class PengRobinsonGas
         //- Critical Temperature [K]
         scalar Tc_;
 
+        //- Critical volume [m^3/kmol]
+        scalar Vc_;
+
+        //- Critical compression factor [-]
+        scalar Zc_;
+
         //- Critical Pressure [Pa]
         scalar Pc_;
 
-        //- Accentric factor [-]
+        //- Acentric factor [-]
         scalar omega_;
 
 
@@ -114,6 +120,8 @@ public:
         (
             const Specie& sp,
             const scalar& Tc,
+            const scalar& Vc,
+            const scalar& Zc,
             const scalar& Pc,
             const scalar& omega
         );
@@ -163,7 +171,7 @@ public:
             //- Return compressibility rho/p [s^2/m^2]
             inline scalar psi(scalar p, scalar T) const;
 
-            //- Return compression factor []
+            //- Return compression factor [-]
             inline scalar Z(scalar p, scalar T) const;
 
             //- Return (cp - cv) [J/(kmol K]

src/thermophysicalModels/specie/equationOfState/PengRobinsonGas/PengRobinsonGasI.H

@@ -21,7 +21,6 @@ License
     You should have received a copy of the GNU General Public License
     along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
 
-
 \*---------------------------------------------------------------------------*/
 
 #include "PengRobinsonGas.H"
@@ -34,12 +33,16 @@ inline Foam::PengRobinsonGas<Specie>::PengRobinsonGas
 (
     const Specie& sp,
     const scalar& Tc,
+    const scalar& Vc,
+    const scalar& Zc,
     const scalar& Pc,
     const scalar& omega
 )
 :
     Specie(sp),
     Tc_(Tc),
+    Vc_(Vc),
+    Zc_(Zc),
     Pc_(Pc),
     omega_(omega)
 {}
@@ -55,9 +58,11 @@ inline Foam::PengRobinsonGas<Specie>::PengRobinsonGas
 )
 :
     Specie(name, pg),
-    Tc_(pg.Tc),
-    Pc_(pg.Pc),
-    omega_(pg.omega)
+    Tc_(pg.Tc_),
+    Pc_(pg.Pc_),
+    Vc_(pg.Vc_),
+    Zc_(pg.Zc_),
+    omega_(pg.omega_)
 {}
 
 
@@ -96,6 +101,7 @@ Foam::PengRobinsonGas<Specie>::New
     );
 }
 
+
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
 template<class Specie>
@@ -214,7 +220,9 @@ inline void Foam::PengRobinsonGas<Specie>::operator+=
     scalar molr2 = pg.nMoles()/this->nMoles();
 
     Tc_ = molr1*Tc_ + molr2*pg.Tc_;
-    Pc_ = molr1*Pc_ + molr2*pg.Pc_;
+    Vc_ = molr1*Vc_ + molr2*pg.Vc_;
+    Zc_ = molr1*Zc_ + molr2*pg.Zc_;
+    Pc_ = specie::RR*Zc_*Tc_/Vc_;
     omega_ = molr1*omega_ + molr2*pg.omega_;
 }
 
@@ -233,7 +241,9 @@ inline void Foam::PengRobinsonGas<Specie>::operator-=
     scalar molr2 = pg.nMoles()/this->nMoles();
 
     Tc_ = molr1*Tc_ - molr2*pg.Tc_;
-    Pc_ = molr1*Pc_ - molr2*pg.Pc_;
+    Vc_ = molr1*Vc_ - molr2*pg.Vc_;
+    Zc_ = molr1*Zc_ - molr2*pg.Zc_;
+    Pc_ = specie::RR*Zc_*Tc_/Vc_;
     omega_ = molr1*omega_ - molr2*pg.omega_;
 }
 
@@ -259,12 +269,18 @@ Foam::PengRobinsonGas<Specie> Foam::operator+
     scalar molr1 = pg1.nMoles()/nMoles;
     scalar molr2 = pg2.nMoles()/nMoles;
 
+    const scalar Tc = molr1*pg1.Tc_ + molr2*pg2.Tc_;
+    const scalar Vc = molr1*pg1.Vc_ + molr2*pg2.Vc_;
+    const scalar Zc = molr1*pg1.Zc_ + molr2*pg2.Zc_;
+
     return PengRobinsonGas<Specie>
     (
         static_cast<const Specie&>(pg1)
       + static_cast<const Specie&>(pg2),
-        molr1*pg1.Tc_ + molr2*pg2.Tc_,
-        molr1*pg1.Pc_ + molr2*pg2.Pc_,
+        Tc,
+        Vc,
+        Zc,
+        specie::RR*Zc*Tc/Vc,
         molr1*pg1.omega_ + molr2*pg2.omega_
     );
 }
@@ -281,12 +297,18 @@ Foam::PengRobinsonGas<Specie> Foam::operator-
     scalar molr1 = pg1.nMoles()/nMoles;
     scalar molr2 = pg2.nMoles()/nMoles;
 
+    const scalar Tc = molr1*pg1.Tc_ + molr2*pg2.Tc_;
+    const scalar Vc = molr1*pg1.Vc_ + molr2*pg2.Vc_;
+    const scalar Zc = molr1*pg1.Zc_ + molr2*pg2.Zc_;
+
     return PengRobinsonGas<Specie>
     (
         static_cast<const Specie&>(pg1)
       - static_cast<const Specie&>(pg2),
-        molr1*pg1.Tc_ - molr2*pg2.Tc_,
-        molr1*pg1.Pc_ - molr2*pg2.Pc_,
+        Tc,
+        Vc,
+        Zc,
+        specie::RR*Zc*Tc/Vc,
         molr1*pg1.omega_ - molr2*pg2.omega_
     );
 }
@@ -303,6 +325,8 @@ Foam::PengRobinsonGas<Specie> Foam::operator*
     (
         s*static_cast<const Specie&>(pg),
         pg.Tc_,
+        pg.Vc_,
+        pg.Zc_,
         pg.Pc_,
         pg.omega_
     );