solver, tutorial +mixClasses changed

applications/solvers/compressible/realFluidPisoFoam/Make/files

@@ -0,0 +1,3 @@
+realFluidPisoFoam.C
+
+EXE = $(FOAM_APPBIN)/realFluidPisoFoam

applications/solvers/compressible/realFluidPisoFoam/Make/options

@@ -0,0 +1,14 @@
+EXE_INC = \
+    -I$(LIB_SRC)/thermophysicalModels/specie/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
+    -I$(LIB_SRC)/turbulenceModels/compressible/turbulenceModel \
+    -I$(LIB_SRC)/finiteVolume/lnInclude
+
+EXE_LIBS = \
+    -lfiniteVolume \
+    -lmeshTools \
+    -lbasicThermophysicalModels \
+    -lspecie \
+    -lcompressibleTurbulenceModel \
+    -lcompressibleRASModels \
+    -lcompressibleLESModels 

applications/solvers/compressible/realFluidPisoFoam/UEqn.H

@@ -0,0 +1,11 @@
+    fvVectorMatrix UEqn
+    (
+        fvm::ddt(rho, U)
+      + fvm::div(phi, U)
+      + turbulence->divDevRhoReff(U)
+    );
+
+    if (momentumPredictor)
+    {
+        solve(UEqn == -fvc::grad(p));
+    }

applications/solvers/compressible/realFluidPisoFoam/createFields.H

@@ -0,0 +1,61 @@
+    Info<< "Reading thermophysical properties\n" << endl;
+
+    autoPtr<basicPsiThermo> pThermo
+    (
+        basicPsiThermo::New(mesh)
+    );
+    basicPsiThermo& thermo = pThermo();
+
+    volScalarField& p = thermo.p();
+    volScalarField& h = thermo.h();
+    const volScalarField& psi = thermo.psi();
+    const volScalarField& drhodh = thermo.drhodh();
+
+    bool realFluid=mesh.solutionDict().subDict("PISO").lookupOrDefault<bool>("realFluid",false);
+
+    volScalarField rho
+    (
+        IOobject
+        (
+            "rho",
+            runTime.timeName(),
+            mesh,
+            IOobject::NO_READ,
+            IOobject::AUTO_WRITE
+        ),
+        thermo.rho()
+    );
+
+    Info<< "\nReading field U\n" << endl;
+    volVectorField U
+    (
+        IOobject
+        (
+            "U",
+            runTime.timeName(),
+            mesh,
+            IOobject::MUST_READ,
+            IOobject::AUTO_WRITE
+        ),
+        mesh
+    );
+
+#   include "compressibleCreatePhi.H"
+
+
+    Info<< "Creating turbulence model\n" << endl;
+    autoPtr<compressible::turbulenceModel> turbulence
+    (
+        compressible::turbulenceModel::New
+        (
+            rho,
+            U,
+            phi,
+            thermo
+        )
+    );
+
+
+    Info<< "Creating field DpDt\n" << endl;
+    volScalarField DpDt =
+        fvc::DDt(surfaceScalarField("phiU", phi/fvc::interpolate(rho)), p);

applications/solvers/compressible/realFluidPisoFoam/hEqn.H

@@ -0,0 +1,12 @@
+{
+    solve
+    (
+        fvm::ddt(rho, h)
+      + fvm::div(phi, h)
+      - fvm::laplacian(turbulence->alphaEff(), h)
+     ==
+        DpDt
+    );
+
+    thermo.correct();
+}

applications/solvers/compressible/realFluidPisoFoam/pEqn.H

@@ -0,0 +1,98 @@
+rho = thermo.rho();
+
+volScalarField rUA = 1.0/UEqn.A();
+U = rUA*UEqn.H();
+
+if (realFluid)
+{
+    phi =
+        fvc::interpolate(rho)*
+        (
+            (fvc::interpolate(U) & mesh.Sf())
+          + fvc::ddtPhiCorr(rUA, rho, U, phi)
+        );
+
+    for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
+    {
+        fvScalarMatrix pEqn
+        (
+            psi*fvm::ddt(p)
+          + drhodh*fvc::ddt(h)
+          + fvc::div(phi)
+          - fvm::laplacian(rho*rUA, p)
+        );
+
+        pEqn.solve();
+
+        if (nonOrth == nNonOrthCorr)
+        {
+            phi += pEqn.flux();
+        }
+    }
+}
+else
+{
+    if (transonic)
+    {
+        surfaceScalarField phid
+        (
+            "phid",
+            fvc::interpolate(psi)
+           *(
+                (fvc::interpolate(U) & mesh.Sf())
+               + fvc::ddtPhiCorr(rUA, rho, U, phi)
+            )
+        );
+
+        for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
+        {
+            fvScalarMatrix pEqn
+            (
+                fvm::ddt(psi, p)
+              + fvm::div(phid, p)
+              - fvm::laplacian(rho*rUA, p)
+            );
+
+            pEqn.solve();
+
+            if (nonOrth == nNonOrthCorr)
+            {
+                phi == pEqn.flux();
+            }
+        }
+    }
+    else
+    {
+        phi =
+        fvc::interpolate(rho)*
+        (
+            (fvc::interpolate(U) & mesh.Sf())
+          + fvc::ddtPhiCorr(rUA, rho, U, phi)
+        );
+
+        for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
+        {
+            fvScalarMatrix pEqn
+            (
+                fvm::ddt(psi, p)
+              + fvc::div(phi)
+              - fvm::laplacian(rho*rUA, p)
+            );
+
+            pEqn.solve();
+
+            if (nonOrth == nNonOrthCorr)
+            {
+                phi += pEqn.flux();
+            }
+        }
+    }
+}
+
+#include "rhoEqn.H"
+#include "compressibleContinuityErrs.H"
+
+U -= rUA*fvc::grad(p);
+U.correctBoundaryConditions();
+
+DpDt = fvc::DDt(surfaceScalarField("phiU", phi/fvc::interpolate(rho)), p);

applications/solvers/compressible/realFluidPisoFoam/realFluidPisoFoam.C

@@ -0,0 +1,97 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright held by original author
+     \\/     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 2 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, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+Application
+    rhoPisoFoam
+
+Description
+    Transient PISO solver for compressible, laminar or turbulent flow.
+   
+    CL: rhoPisoFoam with a changed pressure equation for non-perfect gas fluids
+    CL: see realFluid flag
+
+\*---------------------------------------------------------------------------*/
+
+#include "fvCFD.H"
+#include "basicPsiThermo.H"
+#include "turbulenceModel.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+int main(int argc, char *argv[])
+{
+    #include "setRootCase.H"
+
+    #include "createTime.H"
+    #include "createMesh.H"
+    #include "createFields.H"
+    #include "initContinuityErrs.H"
+    #include "readTimeControls.H"
+    #include "compressibleCourantNo.H"
+    #include "setInitialDeltaT.H"
+
+
+    // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+    Info<< "\nStarting time loop\n" << endl;
+
+    while (runTime.run())
+    {
+        #include "readTimeControls.H"
+        #include "readPISOControls.H"
+        #include "compressibleCourantNo.H"
+        #include "setDeltaT.H"
+
+        runTime++;
+
+        Info<< "Time = " << runTime.timeName() << nl << endl;
+
+        #include "rhoEqn.H"
+        #include "UEqn.H"
+
+        // --- PISO loop
+        for (int corr=1; corr<=nCorr; corr++)
+        {
+            #include "hEqn.H"
+            #include "pEqn.H"
+        }
+
+        turbulence->correct();
+
+        rho = thermo.rho();
+
+        runTime.write();
+
+        Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
+            << "  ClockTime = " << runTime.elapsedClockTime() << " s"
+            << nl << endl;
+    }
+
+    Info<< "End\n" << endl;
+
+    return 0;
+}
+
+
+// ************************************************************************* //

applications/solvers/compressible/rhoPisoFoam/createFields.H

@@ -9,9 +9,6 @@
     volScalarField& p = thermo.p();
     volScalarField& h = thermo.h();
     const volScalarField& psi = thermo.psi();
-    const volScalarField& drhodh = thermo.drhodh();
-
-    bool realFluid=mesh.solutionDict().subDict("PISO").lookupOrDefault<bool>("realFluid",false);
 
     volScalarField rho
     (

applications/solvers/compressible/rhoPisoFoam/pEqn.H

@@ -3,7 +3,36 @@ rho = thermo.rho();
 volScalarField rUA = 1.0/UEqn.A();
 U = rUA*UEqn.H();
 
-if (realFluid)
+if (transonic)
+{
+    surfaceScalarField phid
+    (
+        "phid",
+        fvc::interpolate(psi)
+       *(
+            (fvc::interpolate(U) & mesh.Sf())
+          + fvc::ddtPhiCorr(rUA, rho, U, phi)
+        )
+    );
+
+    for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
+    {
+        fvScalarMatrix pEqn
+        (
+            fvm::ddt(psi, p)
+          + fvm::div(phid, p)
+          - fvm::laplacian(rho*rUA, p)
+        );
+
+        pEqn.solve();
+
+        if (nonOrth == nNonOrthCorr)
+        {
+            phi == pEqn.flux();
+        }
+    }
+}
+else
 {
     phi =
         fvc::interpolate(rho)*
@@ -16,8 +45,7 @@ if (realFluid)
     {
         fvScalarMatrix pEqn
         (
-            psi*fvm::ddt(p)
+            fvm::ddt(psi, p)
-          + drhodh*fvc::ddt(h)
           + fvc::div(phi)
           - fvm::laplacian(rho*rUA, p)
         );
@@ -30,64 +58,6 @@ if (realFluid)
         }
     }
 }
-else
-{
-    if (transonic)
-    {
-        surfaceScalarField phid
-        (
-            "phid",
-            fvc::interpolate(psi)
-           *(
-                (fvc::interpolate(U) & mesh.Sf())
-               + fvc::ddtPhiCorr(rUA, rho, U, phi)
-            )
-        );
-
-        for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
-        {
-            fvScalarMatrix pEqn
-            (
-                fvm::ddt(psi, p)
-              + fvm::div(phid, p)
-              - fvm::laplacian(rho*rUA, p)
-            );
-
-            pEqn.solve();
-
-            if (nonOrth == nNonOrthCorr)
-            {
-                phi == pEqn.flux();
-            }
-        }
-    }
-    else
-    {
-        phi =
-        fvc::interpolate(rho)*
-        (
-            (fvc::interpolate(U) & mesh.Sf())
-          + fvc::ddtPhiCorr(rUA, rho, U, phi)
-        );
-
-        for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
-        {
-            fvScalarMatrix pEqn
-            (
-                fvm::ddt(psi, p)
-              + fvc::div(phi)
-              - fvm::laplacian(rho*rUA, p)
-            );
-
-            pEqn.solve();
-
-            if (nonOrth == nNonOrthCorr)
-            {
-                phi += pEqn.flux();
-            }
-        }
-    }
-}
 
 #include "rhoEqn.H"
 #include "compressibleContinuityErrs.H"

src/thermophysicalModels/specie/Make/files

@@ -12,12 +12,13 @@ $(reactions)/makeChemkinReactions.C
 $(reactions)/makeReactionThermoReactions.C
 $(reactions)/makeLangmuirHinshelwoodReactions.C
 
-$(equationOfState)/aungierRedlichKwong/aungierRedlichKwong.C
+$(equationOfState)/cubicEquationOfState/aungierRedlichKwong/aungierRedlichKwong.C
-$(equationOfState)/pengRobinson/pengRobinson.C
+$(equationOfState)/cubicEquationOfState/pengRobinson/pengRobinson.C
-$(equationOfState)/redlichKwong/redlichKwong.C
+$(equationOfState)/cubicEquationOfState/redlichKwong/redlichKwong.C
-$(equationOfState)/soaveRedlichKwong/soaveRedlichKwong.C
+$(equationOfState)/cubicEquationOfState/soaveRedlichKwong/soaveRedlichKwong.C
-$(equationOfState)/mixtureRedlichKwong/mixtureRedlichKwong.C
+$(equationOfState)/cubicEquationOfState/mixtures/mixtureRedlichKwong/mixtureRedlichKwong.C
-$(equationOfState)/mixturePengRobinson/mixturePengRobinson.C
+$(equationOfState)/cubicEquationOfState/mixtures/mixturePengRobinson/mixturePengRobinson.C
-$(equationOfState)/mixtureSoaveRedlichKwong/mixtureSoaveRedlichKwong.C
+$(equationOfState)/cubicEquationOfState/mixtures/mixtureSoaveRedlichKwong/mixtureSoaveRedlichKwong.C
+
 
 LIB = $(FOAM_LIBBIN)/libspecie

src/thermophysicalModels/specie/equationOfState/cubicEquationOfState/aungierRedlichKwong/aungierRedlichKwong.C

@@ -41,18 +41,6 @@ Germany
 namespace Foam
 {
 
-/* * * * * * * * * * * * * * * Private static data * * * * * * * * * * * * * */
-
-const scalar aungierRedlichKwong::rhoMin_
-(
-    debug::tolerances("aungierRedlichKwongRhoMin", 1e-3)
-);
-
-const scalar aungierRedlichKwong::rhoMax_
-(
-    debug::tolerances("aungierRedlichKwongRhoMax", 1500)
-);
-
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
 aungierRedlichKwong::aungierRedlichKwong(Istream& is)
@@ -62,22 +50,76 @@ aungierRedlichKwong::aungierRedlichKwong(Istream& is)
     Tcrit_(readScalar(is)),
     azentricFactor_(readScalar(is)),
     rhocrit_(readScalar(is)),    
-    a_(0.42747*pow(this->RR,2)*pow(Tcrit_,2)/pcrit_),
+    a0_(0.42747*pow(this->RR,2)*pow(Tcrit_,2)/pcrit_),
     b_(0.08664*this->RR*Tcrit_/pcrit_),
-    c_(this->RR*Tcrit_/(pcrit_+(a_/(this->W()/rhocrit_*(this->W()/rhocrit_+b_))))+b_-this->W()/rhocrit_),
+    c_(this->RR*Tcrit_/(pcrit_+(a0_/(this->W()/rhocrit_*(this->W()/rhocrit_+b_))))+b_-this->W()/rhocrit_),
     n_(0.4986+1.2735*azentricFactor_+0.4754*pow(azentricFactor_,2)),
-        // Starting GUESS for the density by ideal gas law
+    TSave(0.0), 
-        rhostd_(this->rho(Pstd,Tstd,Pstd/(Tstd*this->R())))
+    //CL: Only uses the default values
+    rhoMin_(1e-3),
+    rhoMax_(1500),
+    b2_(pow(b_,2)),
+    b3_(pow(b_,3)),
+    b4_(pow(b_,4)),
+    b5_(pow(b_,5)),
+    c2_(pow(c_,2)),
+    // Starting GUESS for the density by ideal gas law
+    rhostd_(this->rho(Pstd,Tstd,Pstd/(Tstd*this->R())))
 {
     is.check("aungierRedlichKwong::aungierRedlichKwong(Istream& is)");
 }
+//CL: Constructed needed in OpenFOAM 2.x.x
+//CL: Code works fine, but compiling problem in OpenFOAM 1.6.ext
+//CL:  because specie has no constructor using dict
+/*
+aungierRedlichKwong::aungierRedlichKwong(const dictionary& dict)
+:
+    specie(dict),
+    pcrit_(readScalar(dict.subDict("equationOfState").lookup("pCritical"))),
+    Tcrit_(readScalar(dict.subDict("equationOfState").lookup("TCritical"))),
+    azentricFactor_(readScalar(dict.subDict("equationOfState").lookup("azentricFactor"))),
+    rhocrit_(readScalar(dict.subDict("equationOfState").lookup("rhoCritical"))),
+    //CL: rhoMin and rhoMax are only used as boundaries for the bisection method (see rho function)
+    //CL: important: rhoMin and rhoMax are not used as boundary for the newton solver
+    //CL: therefore, rho can be larger than rhoMax and smaller than rhoMin
+    rhoMin_(dict.subDict("equationOfState").lookupOrDefault("rhoMin",1e-3)),
+    rhoMax_(dict.subDict("equationOfState").lookupOrDefault("rhoMax",1500)),
+    a0_(0.42747*pow(this->RR,2)*pow(Tcrit_,2)/pcrit_),
+    b_(0.08664*this->RR*Tcrit_/pcrit_),
+    c_(this->RR*Tcrit_/(pcrit_+(a0_/(this->W()/rhocrit_*(this->W()/rhocrit_+b_))))+b_-this->W()/rhocrit_),
+    n_(0.4986+1.2735*azentricFactor_+0.4754*pow(azentricFactor_,2)),
+    TSave(0.0), 
+    b2_(pow(b_,2)),
+    b3_(pow(b_,3)),
+    b4_(pow(b_,4)),
+    b5_(pow(b_,5)),
+    c2_(pow(c_,2)),
+    // Starting GUESS for the density by ideal gas law
+    rhostd_(this->rho(Pstd,Tstd,Pstd/(Tstd*this->R())))
+{}
 
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
+void Foam::aungierRedlichKwong::write(Ostream& os) const
+{
+    specie::write(os);
+
+    dictionary dict("equationOfState");
+    dict.add("pCritical", pcrit_);
+    dict.add("TCritical", Tcrit_);
+    dict.add("azentricFactor", azentricFactor_);
+    dict.add("rhoCritical", rhocrit_);
+    dict.add("rhoMin", rhoMin_);
+    dict.add("rhoMax", rhoMax_);
+
+    os  << indent << dict.dictName() << dict;
+}
+*/
 // * * * * * * * * * * * * * * * Ostream Operator  * * * * * * * * * * * * * //
 
 Ostream& operator<<(Ostream& os, const aungierRedlichKwong& ark)
 {
-    os  << static_cast<const specie&>(ark)<< tab
+    os  << static_cast<const specie&>(ark)<< token::SPACE 
         << ark.pcrit_ << tab<< ark.Tcrit_<< tab<<ark.azentricFactor_<< tab<<ark.rhocrit_;
 
     os.check("Ostream& operator<<(Ostream& os, const aungierRedlichKwong& st)");

src/thermophysicalModels/specie/equationOfState/cubicEquationOfState/aungierRedlichKwong/aungierRedlichKwong.H

@@ -65,6 +65,9 @@ class aungierRedlichKwong
     public specie
 
 {
+
+protected:
+
     // Private data
         scalar pcrit_;
         scalar Tcrit_; 
@@ -72,21 +75,37 @@ class aungierRedlichKwong
         scalar rhocrit_;
 
     //Aungier Redlich Kwong factors
-        mutable scalar a_;
+        mutable scalar a0_;
         mutable scalar b_;
         mutable scalar c_;
         mutable scalar n_;
 
+    //CL: pow of constants (b_, c_) used in the code e.g. b2_=b*b;
+        mutable scalar b2_;
+        mutable scalar b3_;
+        mutable scalar b4_;
+        mutable scalar b5_;
+        mutable scalar c2_;
+
     //Density @STD, initialise after a, b!
         scalar rhostd_;
 
-    //Static Private data
+    //CL: rhoMin and rhoMax are only used as boundaries for the bisection methode (see rho function)
+    scalar rhoMax_;
+    scalar rhoMin_;
 
-    //should be read from the fvSolution file where rhoMax and rhoMin values must be define ( for rhoSimpleFoam)
+    //CL: Variables to save the values of a, dadT and d2adT2 of the mixture
-    //HR: Don't know, yet. Let's make these static for starters
+    //CL: Variables must corrected for changing temperatures
-        static const scalar rhoMax_;
+    mutable scalar aSave;
-        static const scalar rhoMin_; 
+    mutable scalar daSave;
+    mutable scalar d2aSave;
 
+    //CL: save the temperature for which the save coefficients (amix,dadTmix,d2adT2mix) are correct
+    mutable scalar TSave;
+
+    //Protected functions
+    //CL: function updates the coefficients (aSave, daSave, d2aSave)
+    inline void updateModelCoefficients(const  scalar T) const; 
 
 public:
   
@@ -96,16 +115,15 @@ public:
         //- Construct from components
         inline aungierRedlichKwong
         (
-            const specie& sp, 
+            const specie& sp
-            scalar pcrit, 
-            scalar Tcrit,
-            scalar azentricFactor,
-            scalar rhocrit
         );
 
         //- Construct from Istream
         aungierRedlichKwong(Istream&);
 
+        //- Construct from dictionary
+        //aungierRedlichKwong(const dictionary& dict);
+
         //- Construct as named copy
         inline aungierRedlichKwong(const word& name, const aungierRedlichKwong&);
 
@@ -119,6 +137,36 @@ public:
 
 	inline scalar rhostd() const;
 
+        //CL: Model coefficient a(T)
+	inline scalar a(const  scalar T)const;
+
+        //CL: temperature deriviative of model coefficient a(T)
+	inline scalar dadT(const  scalar T)const;
+
+        //CL: second order temperature deriviative of model coefficient a(T)
+	inline scalar d2adT2(const  scalar T)const;
+
+        //Return Aungier Redlich Kwong factors
+	inline scalar a0()const;
+
+	inline scalar b()const;
+
+	inline scalar c()const;
+
+	inline scalar n()const;
+
+        //CL: return power of constants (b_, c_)
+	inline scalar b2()const;
+
+	inline scalar b3()const;
+
+	inline scalar b4()const;
+
+	inline scalar b5()const;
+
+	inline scalar c2()const;
+
+
         //CL: Equation of state
         inline scalar p(const scalar rho, const scalar T) const;
 
@@ -184,10 +232,31 @@ public:
         ) const;
 
 
+        // I-O
+
+            //- Write to Ostream
+            //void write(Ostream& os) const;
+
+
     // Member operators
 
+        inline void operator+=(const aungierRedlichKwong&);
+
+
     // Friend operators
 
+        inline friend aungierRedlichKwong operator+
+        (
+            const aungierRedlichKwong&,
+            const aungierRedlichKwong&
+        );
+
+        inline friend aungierRedlichKwong operator*
+        (
+            const scalar s,
+            const aungierRedlichKwong&
+        );
+
     // Ostream Operator
 
         friend Ostream& operator<<(Ostream&, const aungierRedlichKwong&);

src/thermophysicalModels/specie/equationOfState/cubicEquationOfState/aungierRedlichKwong/aungierRedlichKwongI.H

@@ -42,24 +42,10 @@ namespace Foam
 // Construct from components
 inline aungierRedlichKwong::aungierRedlichKwong
 (
-    const specie& sp,
+    const specie& sp
-    scalar pcrit, 
-    scalar Tcrit,
-    scalar azentricFactor,
-    scalar rhocrit
 )
 :
-    specie(sp),    
+    specie(sp)
-    pcrit_(pcrit),
-    Tcrit_(Tcrit),
-    azentricFactor_(azentricFactor),
-    rhocrit_(rhocrit),
-    a_(0.42747*pow(this->RR,2)*pow(Tcrit_,2)/pcrit_),
-    b_(0.08664*this->RR*Tcrit_/pcrit_),
-    c_(this->RR*Tcrit_/(pcrit_+(a_/(this->W()/rhocrit_*(this->W()/rhocrit_+b_))))+b_-this->W()/rhocrit_),
-    n_(0.4986+1.2735*azentricFactor_+0.4754*pow(azentricFactor_,2)),
-    // Starting GUESS for the density by ideal gas law
-        rhostd_(this->rho(Pstd,Tstd,Pstd/(Tstd*this->R())))
 {}
 
 
@@ -71,11 +57,16 @@ inline aungierRedlichKwong::aungierRedlichKwong(const word& name, const aungierR
     Tcrit_(pg.Tcrit_),
     azentricFactor_(pg.azentricFactor_),
     rhocrit_(pg.rhocrit_),    
-    a_(pg.a_),
+    a0_(pg.a0_),
     b_(pg.b_),
     c_(pg.c_),
     n_(pg.n_),
-    rhostd_(pg.rhostd_)
+    rhostd_(pg.rhostd_),
+    b2_(pg.b2_),
+    b3_(pg.b3_),
+    b4_(pg.b4_),
+    b5_(pg.b5_),
+    c2_(pg.c2_)
 {}
 
 
@@ -95,18 +86,128 @@ inline autoPtr<aungierRedlichKwong> aungierRedlichKwong::New(Istream& is)
 
 
 // * * * * * * * * * * * * *  Member Functions  * * * * * * * * * * * //
+
 inline scalar aungierRedlichKwong::rhostd()const
 {
     return rhostd_;
 }
 
+inline void aungierRedlichKwong::updateModelCoefficients(const scalar T)const
+{
+    aSave=a0_*pow(T/Tcrit_,-n());
+    daSave=-a0_*n()*pow(T/Tcrit_,-n())/T;
+    d2aSave=a0_*(n() *n()+n())/(T*T)*pow(T/Tcrit_,-n());
+
+    //CL: saving the temperature at which the coefficients are valid
+    TSave=T;
+}
+
+
+//CL: Model coefficient a(T)
+inline scalar aungierRedlichKwong::a(const scalar T)const
+{
+    //CL: check if a has already been calculated for this temperature
+    if(TSave==T)
+    { 
+        return aSave;  
+    }
+    //CL: If not, recalculate a(T), dadT(T) and d2adT2(T)
+    else
+    {
+        updateModelCoefficients(T);
+        return aSave;
+    }     
+}
+
+
+//CL: temperature deriviative of model coefficient a(T)
+inline scalar aungierRedlichKwong::dadT(const scalar T)const
+{
+    // check if a has already been calculated for this temperature
+    if(TSave==T)
+    {
+        return daSave;  
+    }
+    //CL: If not, recalculate a(T), dadT(T) and d2adT2(T)
+    else
+    {
+        updateModelCoefficients(T);
+        return daSave;
+    }    
+}
+
+
+//CL: second order temperature deriviative of model coefficient a(T)
+inline scalar aungierRedlichKwong::d2adT2(const scalar T)const
+{
+    // check if a has already been calculated for this temperature
+    if(TSave==T)
+    {
+        return d2aSave;  
+    }
+    //CL: If not, recalculate a(T), dadT(T) and d2adT2(T)
+    else
+    {
+        updateModelCoefficients(T);
+        return d2aSave;
+    }    
+}
+
+//Aungier Redlich Kwong factors
+inline scalar aungierRedlichKwong::a0()const
+{
+    return a0_;
+}
+
+
+inline scalar aungierRedlichKwong::b()const
+{
+    return b_;
+}
+
+inline scalar aungierRedlichKwong::c()const
+{
+    return c_;
+}
+
+inline scalar aungierRedlichKwong::n()const
+{
+    return n_;
+}
+
+//CL: pow of constants (b(), c()) used in the code e.g. b2_=b*b;
+inline scalar aungierRedlichKwong::b2()const
+{
+    return b2_;
+}
+
+inline scalar aungierRedlichKwong::b3()const
+{
+    return b3_;
+}
+
+inline scalar aungierRedlichKwong::b4()const
+{
+    return b4_;
+}
+
+inline scalar aungierRedlichKwong::b5()const
+{
+    return b5_;
+}
+
+inline scalar aungierRedlichKwong::c2()const
+{
+    return c2_;
+}
+
 
 //returns the pressure for a given density and temperature    
 inline scalar aungierRedlichKwong::p(const scalar rho,const scalar T) const
 {
     scalar Vm = this->W()/rho;
-    return this->RR*T/(Vm-b_+c_)
+    return this->RR*T/(Vm-b()+c())
-        -a_*pow(T,-n_)/(pow(Tcrit_,-n_)*Vm*(Vm+b_));
+               -a(T)/(Vm*(Vm+b()));
 }
 
 
@@ -114,9 +215,12 @@ inline scalar aungierRedlichKwong::p(const scalar rho,const scalar T) const
 //(molar values)
 inline scalar aungierRedlichKwong::dpdv(const scalar rho, const scalar T) const
 {
-    scalar Vm = this->W()/rho;
+   scalar Vm = this->W()/rho;
-    return  a_*pow(T,-n_)*pow(Tcrit_,n_)*(b_+2*Vm)/(pow(Vm,2)*pow((b_+Vm),2))
+    scalar Vm2 = Vm*Vm;
-        -this->RR*T/(pow((b_-Vm-c_),2));
+
+    return  (a(T)*(b3()-2*b2()*c()+b()*(c()+Vm)*(c()-3*Vm)+2*Vm*pow(c()+Vm,2))
+                -this->RR*T*Vm2*(b2()+2*b()*Vm+Vm2))
+                    /(Vm2*pow(b()-c()-Vm,2)*pow(b()+Vm,2));
 }
 
 
@@ -125,8 +229,7 @@ inline scalar aungierRedlichKwong::dpdv(const scalar rho, const scalar T) const
 inline scalar aungierRedlichKwong::dpdT(const scalar rho, const scalar T) const
 {
     scalar Vm = this->W()/rho;
-    return a_*n_*pow(T,-n_-1)*pow(Tcrit_,n_)/(Vm*(Vm+b_))
+    return this->RR/(Vm-b()+c())-dadT(T)/(Vm*(Vm+b()));
-        -this->RR/(b_-Vm-c_);
 }
 
 
@@ -156,8 +259,7 @@ inline scalar aungierRedlichKwong::integral_p_dv
 ) const
 {
     scalar Vm = this->W()/rho;
-    return -pow((T/Tcrit_),-n_)*(a_*log(Vm)/(b_)
+    return this->RR*T*log(-b()+c()+Vm)+a(T)*log(b()+Vm)/b()-a(T)*log(Vm)/b();
-        -a_*log(b_+Vm)/b_)+this->RR*T*log(Vm-b_+c_);
 }
 
 
@@ -170,8 +272,7 @@ inline scalar aungierRedlichKwong::integral_dpdT_dv
 ) const
 {
     scalar Vm = this->W()/rho;
-    return pow(T,-n_-1)*pow(Tcrit_,n_)*(a_*n_*log(Vm)/(b_) 
+    return this->RR*log(-b()+c()+Vm)+dadT(T)*log(b()+Vm)/b()-dadT(T)*log(Vm)/b();
-        -a_*n_*log(b_+Vm)/(b_))+ this->RR*log(-b_+c_+Vm);
 }
 
 
@@ -179,7 +280,7 @@ inline scalar aungierRedlichKwong::integral_dpdT_dv
 inline scalar aungierRedlichKwong::d2pdT2(const scalar rho,const scalar T) const
 {
     scalar Vm = this->W()/rho;
-    return -a_*n_*pow(T,-n_-2)*pow(Tcrit_,n_)*(n_+1)/(Vm*(b_+Vm));
+    return -d2adT2(T)/(Vm*(Vm+b()));
 }
 
 
@@ -187,9 +288,14 @@ inline scalar aungierRedlichKwong::d2pdT2(const scalar rho,const scalar T) const
 inline scalar aungierRedlichKwong::d2pdv2(const scalar rho,const scalar T) const
 {
     scalar Vm = this->W()/rho;
-    return -2*a_*pow(T,-n_)*pow(Tcrit_,n_)*(pow(b_,2)
+    scalar Vm2 = Vm*Vm;
-        +3*b_*Vm+3*pow(Vm,2))/(pow(Vm,3)*pow((b_+Vm),3))
+    scalar Vm3 = Vm*Vm*Vm;
-        -2*this->RR*T/(pow((b_-Vm-c_),3));
+
+    return -2*(a(T)*(b5()-3*b4()*c()+3*b3()*(c2()-c()*Vm-Vm2)
+               -b2()*(c()+Vm)*(c2()-7*c()*Vm+Vm2)+3*b()*Vm*pow(c()+Vm,2)*(2*Vm-c())
+                   -3*Vm2*pow(c()+Vm,3))+this->RR*T*Vm3*(b3()+3*b2()*Vm
+                        +3*b()*Vm2+Vm3))
+                            /(Vm3*pow(b()-c()-Vm,3)*pow(b()+Vm,3));
 }
 
 
@@ -211,8 +317,11 @@ inline scalar aungierRedlichKwong::d2vdT2(const scalar rho, const scalar T) cons
 inline scalar aungierRedlichKwong::d2pdvdT(const scalar rho, const scalar T) const
 {
     scalar Vm = this->W()/rho;
-    return -a_*n_*pow(T,-n_-1)*pow(Tcrit_,n_)*(b_+2*Vm)/(pow(Vm,2)*pow((b_+Vm),2))
+    scalar Vm2 = Vm*Vm;
-        -this->RR/(pow((b_-c_-Vm),2));
+
+    return (dadT(T)*(b3()-2*b2()*c()+b()*(c()+Vm)*(c()-3*Vm)+2*Vm*pow(c()+Vm,2))
+               -this->RR*Vm2*(b2()+2*b()*Vm+Vm2))
+                   /(Vm2*pow(b()-c()-Vm,2)*pow(b()+Vm,2));
 }
 
 
@@ -225,8 +334,7 @@ inline scalar aungierRedlichKwong::integral_d2pdT2_dv
 ) const 
 {
     scalar Vm = this->W()/rho;
-    return pow(T,-n_-2)*pow(Tcrit_,n_)*(a_*n_*(n_+1)*log(b_+Vm)
+    return d2adT2(T)*log(b()+Vm)/b()-d2adT2(T)*log(Vm)/b();
-        /b_-a_*n_*(1+n_)*log(Vm)/b_);
 }
 
 
@@ -327,7 +435,7 @@ inline scalar aungierRedlichKwong::rho(
        {
            FatalErrorIn
            (
-               "inline scalar redlichKwong::rho(const scalar p, const scalar T, const scalar rho0) const "
+               "inline scalar aungierRedlichKwong::rho(const scalar p, const scalar T, const scalar rho0) const "
            )   << "Maximum number of iterations exceeded"
                << abort(FatalError);
        }
@@ -362,8 +470,37 @@ inline scalar aungierRedlichKwong::Z( const scalar p, const scalar T,const scala
 
 // * * * * * * * * * * * * * * * Member Operators  * * * * * * * * * * * * * //
 
+
+inline void aungierRedlichKwong::operator+=(const aungierRedlichKwong& pr)
+{
+    specie::operator+=(pr);
+}
+
+
 // * * * * * * * * * * * * * * * Friend Operators  * * * * * * * * * * * * * //
 
+inline aungierRedlichKwong operator+
+(
+    const aungierRedlichKwong& pr1,
+    const aungierRedlichKwong& pr2
+)
+{
+    return aungierRedlichKwong
+    (
+        static_cast<const specie&>(pr1)
+      + static_cast<const specie&>(pr2)
+    );
+}
+
+inline aungierRedlichKwong operator*
+(
+    const scalar s,
+    const aungierRedlichKwong& pr
+)
+{
+    return aungierRedlichKwong(s*static_cast<const specie&>(pr));
+}
+
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 } // End namespace Foam

src/thermophysicalModels/specie/equationOfState/cubicEquationOfState/mixtures/README

@@ -0,0 +1,11 @@
+CL: Real gas mixture classes using mixture models.
+CL: other models are possible 
+
+
+CL: Mixtures based on a pseudo critical point approach
+--> mixtureAungierRedlichKwong
+
+CL: Mixtures using the van der waals mixing rule
+--> mixtureSoaveRedlichKwong
+--> mixturePengRobinson
+--> mixtureRedlichKwong

src/thermophysicalModels/specie/equationOfState/cubicEquationOfState/mixtures/mixtureAungierRedlichKwong/mixtureAungierRedlichKwong.C

@@ -23,7 +23,7 @@ License
     Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
 
 Description
-    Peng Robinson equation of state.
+    Mixture Aungier Redlich Kwong equation of state.
 
 Author
 Christian Lucas
@@ -33,7 +33,7 @@ Germany
 
 \*---------------------------------------------------------------------------*/
 
-#include "pengRobinson.H"
+#include "mixtureAungierRedlichKwong.H"
 #include "IOstreams.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@@ -41,44 +41,41 @@ Germany
 namespace Foam
 {
 
-/* * * * * * * * * * * * * * * Private static data * * * * * * * * * * * * * */
-
-const scalar pengRobinson::rhoMin_
-(
-    debug::tolerances("pengRobinsonRhoMin", 1e-3)
-);
-
-const scalar pengRobinson::rhoMax_
-(
-    debug::tolerances("pengRobinsonRhoMax", 1500)
-);
-
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
-pengRobinson::pengRobinson(Istream& is)
+mixtureAungierRedlichKwong::mixtureAungierRedlichKwong(Istream& is)
 :
-    specie(is),
+    aungierRedlichKwong(is),
-    pcrit_(readScalar(is)),
+    Vcrit_(this->W()/rhocrit_),
-    Tcrit_(readScalar(is)),
+    Zcrit_(pcrit_/((this->R()*Tcrit_*rhocrit_)))
-    azentricFactor_(readScalar(is)),
-    a0_(0.457235*pow(this->RR,2)*pow(Tcrit_,2)/pcrit_),
-    b_(0.077796*this->RR*Tcrit_/pcrit_), 
-    n_(0.37464+1.54226*azentricFactor_-0.26992*pow(azentricFactor_,2)),
-    TSave(0.0), 
-    rhostd_(this->rho(Pstd,Tstd,Pstd/(Tstd*this->R())))  	
 {
-    is.check("pengRobinson::pengRobinson(Istream& is)");
+    is.check("mixtureAungierRedlichKwong::mixtureAungierRedlichKwong(Istream& is)");
 }
+//CL: Constructed needed in OpenFOAM 2.x.x
+//CL: Code works fine, but compiling problem in OpenFOAM 1.6.ext
+//CL:  because specie has no constructor using dict
+/*
+mixtureAungierRedlichKwong::mixtureAungierRedlichKwong(const dictionary& dict)
+:
+    aungierRedlichKwong(dict),
+    Vcrit_(this->W()/rhocrit_),
+    Zcrit_(pcrit_/((this->R()*Tcrit_*rhocrit_)))
+{}
 
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
+void Foam::mixtureAungierRedlichKwong::write(Ostream& os) const
+{
+    aungierRedlichKwong::write(os);
+}
+*/
 // * * * * * * * * * * * * * * * Ostream Operator  * * * * * * * * * * * * * //
 
-Ostream& operator<<(Ostream& os, const pengRobinson& pr)
+Ostream& operator<<(Ostream& os, const mixtureAungierRedlichKwong& ark)
 {
-    os  << static_cast<const specie&>(pr)<< tab
+    os  << static_cast<const specie&>(ark)<< token::SPACE;
-        << pr.pcrit_ << tab<< pr.Tcrit_<< tab << pr.azentricFactor_;
 
-    os.check("Ostream& operator<<(Ostream& os, const pengRobinson& st)");
+    os.check("Ostream& operator<<(Ostream& os, const mixtureAungierRedlichKwong& st)");
     return os;
 }
 

src/thermophysicalModels/specie/equationOfState/cubicEquationOfState/mixtures/mixtureAungierRedlichKwong/mixtureAungierRedlichKwong.H

@@ -0,0 +1,164 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright held by original author
+     \\/     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 2 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, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+Class
+    Foam::mixtureAungierRedlichKwong
+
+Description
+    Mixture Aungier Redlich Kwong equation of state.
+
+    Mixture based on a pseudo critical point approach. 
+    
+    For further information, see:
+    BOOK Title: The Properties of Gases And Liquids, Fifth Edition, McGraw-Hill, Chapter 5
+    Authors: B. Poling; J. Prausnitz; J. O'Connell
+
+
+SourceFiles
+    mixtureAungierRedlichKwongI.H
+    mixtureAungierRedlichKwong.C
+
+Author
+Christian Lucas
+Institut für Thermodynamik
+Technische Universität Braunschweig 
+Germany
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef mixtureAungierRedlichKwong_H
+#define mixtureAungierRedlichKwong_H
+
+#include "specie.H"
+#include "aungierRedlichKwong.H"
+#include "autoPtr.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                           Class mixtureAungierRedlichKwong Declaration
+\*---------------------------------------------------------------------------*/
+
+class mixtureAungierRedlichKwong
+:
+    public aungierRedlichKwong
+
+{
+
+protected:
+
+    // data at the critical point
+    scalar Zcrit_;
+        
+    //CL. molar volume at the critical point
+    scalar Vcrit_;
+
+    //CL: save the concentrations of each component of the mixture
+    //CL: needs to be multiplied by this->W() to get the molar fractions
+    //mutable DynamicList<scalar> weigths;
+    
+    //CL: saves a pointer to the pure component classes of the mixture
+   // mutable DynamicList<mixtureAungierRedlichKwong*> mixtureComponents;
+
+
+public:
+  
+
+    // Constructors
+
+        //- Construct from components
+        inline mixtureAungierRedlichKwong
+        (
+            const aungierRedlichKwong& sp, 
+            scalar Tcrit,
+            scalar azentricFactor,
+            scalar Vcrit,
+            scalar Zcrit,
+            scalar rhoMin,
+            scalar rhoMax
+        );
+
+        //- Construct from Istream
+        mixtureAungierRedlichKwong(Istream&);
+
+        //- Construct from dictionary
+        //mixtureAungierRedlichKwong(const dictionary& dict);
+
+        //- Construct as named copy
+        inline mixtureAungierRedlichKwong(const word& name, const mixtureAungierRedlichKwong&);
+
+        //- Construct and return a clone
+        inline autoPtr<mixtureAungierRedlichKwong> clone() const;
+
+        // Selector from Istream
+        inline static autoPtr<mixtureAungierRedlichKwong> New(Istream& is);
+
+   // Member functions
+
+        // I-O
+
+            //- Write to Ostream
+            //void write(Ostream& os) const;
+
+
+  // Member operators
+
+        inline void operator+=(const mixtureAungierRedlichKwong&);
+
+  // Friend operators
+
+        inline friend mixtureAungierRedlichKwong operator+
+        (
+            const mixtureAungierRedlichKwong&,
+            const mixtureAungierRedlichKwong&
+        );
+
+        inline friend mixtureAungierRedlichKwong operator*
+        (
+            const scalar s,
+            const mixtureAungierRedlichKwong&
+        );
+
+    // Ostream Operator
+
+        friend Ostream& operator<<(Ostream&, const mixtureAungierRedlichKwong&);
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#include "mixtureAungierRedlichKwongI.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/thermophysicalModels/specie/equationOfState/cubicEquationOfState/mixtures/mixtureAungierRedlichKwong/mixtureAungierRedlichKwongI.H

@@ -0,0 +1,212 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright held by original author
+     \\/     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 2 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, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+Author
+Christian Lucas
+Institut für Thermodynamik
+Technische Universität Braunschweig 
+Germany
+
+\*---------------------------------------------------------------------------*/
+
+#include "mixtureAungierRedlichKwong.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+// Construct from components
+inline mixtureAungierRedlichKwong::mixtureAungierRedlichKwong
+(
+    const aungierRedlichKwong& sp,
+    scalar Tcrit,
+    scalar azentricFactor,
+    scalar Vcrit,
+    scalar Zcrit,
+    scalar rhoMin,
+    scalar rhoMax
+)
+:
+    aungierRedlichKwong(sp),    
+    Vcrit_(Vcrit),
+    Zcrit_(Zcrit)
+{
+    //CL: Saving critical data
+    Tcrit_=Tcrit;
+    pcrit_=Zcrit*this->RR*Tcrit/Vcrit;  
+    rhocrit_=this->W()/Vcrit_;    
+    azentricFactor_=azentricFactor;
+
+    //CL: calculating the aungier redlich kwong coefficience
+    a0_=0.42747*pow(this->RR,2)*pow(Tcrit_,2)/pcrit_;
+    b_=0.08664*this->RR*Tcrit_/pcrit_;
+    c_=this->RR*Tcrit_/(pcrit_+(a0_/(this->W()/rhocrit_*(this->W()/rhocrit_+b_))))+b_-this->W()/rhocrit_;
+    n_=0.4986+1.2735*azentricFactor_+0.4754*pow(azentricFactor_,2);
+
+    //CL: set rhoMax and rhoMin
+    rhoMin_=rhoMin;
+    rhoMax_=rhoMax;
+
+    //CL: 
+    b2_=pow(b_,2);
+    b3_=pow(b_,3);
+    b4_=pow(b_,4);
+    b5_=pow(b_,5);
+    c2_=pow(c_,2);
+
+    //CL: Starting GUESS for the density by ideal gas law
+    rhostd_=this->rho(Pstd,Tstd,Pstd/(Tstd*this->R()));
+}
+
+
+
+
+// Construct as named copy
+inline mixtureAungierRedlichKwong::mixtureAungierRedlichKwong(const word& name, const mixtureAungierRedlichKwong& pg)
+:
+    aungierRedlichKwong(name, pg),
+    Vcrit_(pg.Vcrit_),
+    Zcrit_(pg.Zcrit_)
+{}
+
+
+
+// Construct and return a clone
+inline autoPtr<mixtureAungierRedlichKwong> mixtureAungierRedlichKwong::clone() const
+{
+    return autoPtr<mixtureAungierRedlichKwong>(new mixtureAungierRedlichKwong(*this));
+}
+
+
+// Selector from Istream
+inline autoPtr<mixtureAungierRedlichKwong> mixtureAungierRedlichKwong::New(Istream& is)
+{
+    return autoPtr<mixtureAungierRedlichKwong>(new mixtureAungierRedlichKwong(is));
+}
+
+
+
+// * * * * * * * * * * * * * * * Member Operators  * * * * * * * * * * * * * //
+
+
+inline void mixtureAungierRedlichKwong::operator+=(const mixtureAungierRedlichKwong& pr)
+{
+    scalar molr1 = this->nMoles();
+
+    aungierRedlichKwong::operator+=(pr);
+
+    molr1 /= this->nMoles();
+    scalar molr2 = pr.nMoles()/this->nMoles();
+
+    //CL: calculating new critical point data
+    Tcrit_=molr1*Tcrit_+molr2*pr.Tcrit_;
+    Zcrit_=molr1*Zcrit_+molr2*pr.Zcrit_;
+    Vcrit_=molr1*Vcrit_+molr2*pr.Vcrit_;
+    pcrit_=Zcrit_*this->RR*Tcrit_/Vcrit_;
+    rhocrit_=this->W()/Vcrit_;
+    
+    //CL: calculating new azentric factor
+    azentricFactor_=molr1*azentricFactor_+molr2*pr.azentricFactor_;
+
+    //CL: set rhoMax and rhoMin
+    rhoMin_=min(rhoMin_, pr.rhoMin_);
+    rhoMax_=max(rhoMax_, pr.rhoMax_);
+
+    //CL: calculating the aungier redlich kwong coefficience
+    a0_=0.42747*pow(this->RR,2)*pow(Tcrit_,2)/pcrit_;
+    b_=0.08664*this->RR*Tcrit_/pcrit_;
+    c_=this->RR*Tcrit_/(pcrit_+(a0_/(this->W()/rhocrit_*(this->W()/rhocrit_+b_))))+b_-this->W()/rhocrit_;
+    n_=0.4986+1.2735*azentricFactor_+0.4754*pow(azentricFactor_,2);
+
+    //CL: 
+    b2_=pow(b_,2);
+    b3_=pow(b_,3);
+    b4_=pow(b_,4);
+    b5_=pow(b_,5);
+    c2_=pow(c_,2);
+
+    //CL: Starting GUESS for the density by ideal gas law
+    rhostd_=this->rho(Pstd,Tstd,Pstd/(Tstd*this->R()));
+}
+
+
+// * * * * * * * * * * * * * * * Friend Operators  * * * * * * * * * * * * * //
+
+inline mixtureAungierRedlichKwong operator+
+(
+    const mixtureAungierRedlichKwong& pr1,
+    const mixtureAungierRedlichKwong& pr2
+)
+{
+
+    aungierRedlichKwong ark
+    (
+        static_cast<const aungierRedlichKwong&>(pr1)
+       +static_cast<const aungierRedlichKwong&>(pr2)
+    );
+
+    return mixtureAungierRedlichKwong
+           (
+               ark,
+               pr1.nMoles()/ark.nMoles()*pr1.Tcrit_
+              +pr2.nMoles()/ark.nMoles()*pr2.Tcrit_,
+               pr1.nMoles()/ark.nMoles()*pr1.azentricFactor_
+              +pr2.nMoles()/ark.nMoles()*pr2.azentricFactor_,
+               pr1.nMoles()/ark.nMoles()*pr1.Vcrit_
+              +pr2.nMoles()/ark.nMoles()*pr2.Vcrit_,
+               pr1.nMoles()/ark.nMoles()*pr1.Zcrit_
+              +pr2.nMoles()/ark.nMoles()*pr2.Zcrit_,
+               min(pr1.rhoMin_,pr2.rhoMin_),
+               max(pr1.rhoMax_,pr2.rhoMax_)
+           );
+}
+
+
+inline mixtureAungierRedlichKwong operator*
+(
+    const scalar s,
+    const mixtureAungierRedlichKwong& pr
+)
+{
+    return mixtureAungierRedlichKwong
+           (
+               s*static_cast<const aungierRedlichKwong&>(pr),
+               pr.Tcrit_, 
+               pr.azentricFactor_, 
+               pr.Vcrit_, 
+               pr.Zcrit_, 
+               pr.rhoMin_, 
+               pr.rhoMax_
+           );
+}
+
+
+
+
+} // End namespace Foam
+
+// ************************************************************************* //

src/thermophysicalModels/specie/equationOfState/cubicEquationOfState/mixtures/mixturePengRobinson/mixturePengRobinson.C

@@ -0,0 +1,157 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright held by original author
+     \\/     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 2 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, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+Description
+    Peng Robionson equation of state for mixtures.
+
+Author
+Christian Lucas
+Institut für Thermodynamik
+Technische Universität Braunschweig 
+Germany
+
+\*---------------------------------------------------------------------------*/
+
+#include "mixturePengRobinson.H"
+#include "IOstreams.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+mixturePengRobinson::mixturePengRobinson(Istream& is)
+:
+    pengRobinson(is),
+    numOfComp(1),
+    singleComponent(1),
+    //CL: no real gas mixture correction when stream constructor is used
+    realMixtureCorr_(false)
+{ 
+    //CL: set size of weigths, mixtureComponents ... to 10,
+    //CL: when more mixture components are used
+    //CL: size of the DynamicLis increases automatically
+    weigths.setSize(10);
+    mixtureComponents.setSize(10);
+
+    //Save a pointer of this object in the mixtureComponents array
+    mixtureComponents[0]=this;
+
+    is.check("mixturePengRobinson::mixturePengRobinson(Istream& is)");
+}
+//CL: Constructed needed in OpenFOAM 2.x.x
+//CL: Code works fine, but compiling problem in OpenFOAM 1.6.ext
+//CL:  because specie has no constructor using dict
+/*
+mixturePengRobinson::mixturePengRobinson(const dictionary& dict)
+:
+    pengRobinson(dict),
+    numOfComp(1),
+    singleComponent(1),
+    realMixtureCorr_(dict.subDict("equationOfState").lookupOrDefault("realMixtureCorrection",false))
+{
+    //CL: set size of weigths, mixtureComponents ... to 10,
+    //CL: when more mixture components are used
+    //CL: size of the DynamicLis increases automatically
+    weigths.setSize(10);
+    mixtureComponents.setSize(10);
+
+    //Save a pointer of this object in the mixtureComponents array
+    mixtureComponents[0]=this;
+
+    if(realMixtureCorr_==true)
+    {
+        //CL: reads number of mixture components
+        //CL: this number is needed to calculate the number of correction coefficients
+        nCom_=dict.subDict("equationOfState").lookupOrDefault("numberOfComponents",1);
+
+            if (nCom_<1||nCom_==1)
+            {
+                FatalErrorIn
+                (
+                    "mixturePengRobinson::mixturePengRobinson(const dictionary& dict)"
+                )   << "mixture has only one component or number is not defined correctly, "
+                    << "recheck numberOfComponents in the  thermophysicalProperties dict of your case"
+                    << abort(FatalError);
+            }
+
+        label i;
+        label j;
+        label k=0;
+ 
+        // CL: saves the number of mixtureCorrectionCoefficients needed for this mixture
+        // CL: need to be set to 1
+        label sizeOfVector_=1;
+
+        //CL: size of the vector depends on the number of components
+        //CL: determine the size of the vector    
+        for (i=3; i<=nCom_;i++)
+        {   
+            sizeOfVector_=sizeOfVector_+(i-1);
+        }
+ 
+        //CL: setting the size 
+        realMixtureCorrCoef_.setSize(sizeOfVector_);
+       
+        //CL: Reading the real mixture correction coefficients
+        //CL: Naming convention e.g for 3 component mixture:
+        //CL: mixtureCorrectionCoefficient_12, mixtureCorrectionCoefficient_13, mixtureCorrectionCoefficient_23
+        for(i=1;i<=nCom_-1;i++)
+        {
+            for(j=i+1;j<=nCom_;j++)
+            {
+                realMixtureCorrCoef_[k]=dict.subDict("equationOfState")
+                    .lookupOrDefault("mixtureCorrectionCoefficient_" + Foam::name(i) + Foam::name(j),0.0);
+                k++;
+            }
+        }
+    }
+}
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+void Foam::mixturePengRobinson::write(Ostream& os) const
+{
+    pengRobinson::write(os);
+
+}
+*/
+// * * * * * * * * * * * * * * * Ostream Operator  * * * * * * * * * * * * * //
+
+Ostream& operator<<(Ostream& os, const mixturePengRobinson& pr)
+{
+    os  << static_cast<const pengRobinson&>(pr)<< token::SPACE ;
+
+    os.check("Ostream& operator<<(Ostream& os, const mixturePengRobinson& st)");
+    return os;
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// ************************************************************************* //

src/thermophysicalModels/specie/equationOfState/cubicEquationOfState/mixtures/mixturePengRobinson/mixturePengRobinson.H

@@ -83,15 +83,6 @@ protected:
     //CL: needs to be multiplied by this->W() to get the molar fractions
     mutable DynamicList<scalar> weigths;
     
-    //CL: stores the values of a(T) of all componentes
-    mutable DynamicList<scalar> aComponents;
-
-    //CL: stores the values of the derivatives da/dT of all componentes
-    mutable DynamicList<scalar> daComponents;
-
-    //CL: stores the values of the second order derivatives d2a/dT2 of all componentes
-    mutable DynamicList<scalar> d2aComponents;
-    
     //CL: saves a pointer to the pure component classes of the mixture
     mutable DynamicList<mixturePengRobinson*> mixtureComponents;
    
@@ -101,8 +92,23 @@ protected:
     //Protected functions
 
     //CL: function updates the coefficients (aSave, daSave, d2aSave, b) of the mixture 
+    //CL: this is the function with the mixing rule
     inline void updateModelCoefficients(const  scalar T) const; 
 
+    //CL: Variables used in real gas mixture correction
+
+        //CL: If true, the real gas mixture correction is used
+        mutable bool  realMixtureCorr_;
+ 
+        //CL: number of mixture components, needed to calculate the mixture correction factors needed
+        //CL: do not mistake this variable with numOfComp,
+        //CL: numOfComp is a counter to counts the number of components while the mixture is constructed
+        mutable label nCom_;
+ 
+        //CL: stores real mixture correction coefficients
+        mutable DynamicList<scalar> realMixtureCorrCoef_;
+
+
 public:
 
     // Constructors
@@ -119,7 +125,9 @@ public:
             scalar b,
             scalar Tcrit,
             scalar n,
-            scalar rhostd
+            scalar rhostd,
+            scalar rhoMin,
+            scalar rhoMax
         );
 
         //- Construct from components
@@ -129,12 +137,17 @@ public:
             const pengRobinson& pr,
             label numOfComp,
             DynamicList<scalar> weigths,
-            DynamicList<mixturePengRobinson*> mixtureComponents
+            DynamicList<mixturePengRobinson*> mixtureComponents,
+            scalar rhoMin,
+            scalar rhoMax
         );
 
         //- Construct from Istream
         mixturePengRobinson(Istream&);
 
+        //- Construct from dictionary
+        //mixturePengRobinson(const dictionary& dict);
+
         //- Construct as named copy
         inline mixturePengRobinson(const word& name, const mixturePengRobinson&);
 
@@ -146,6 +159,13 @@ public:
 
    // Member functions
 
+        //CL: mixture coefficience
+        inline scalar amix(const scalar T, const label i, const  label j) const;
+                 
+        inline scalar dadTmix(const scalar T, const label i, const  label j) const;
+
+        inline scalar d2adT2mix(const scalar T, const label i, const  label j) const;
+
         inline scalar p(const scalar rho, const scalar T) const;
 
         //first order derivatives 
@@ -209,6 +229,11 @@ public:
             const scalar rho0
         ) const;
 
+        // I-O
+
+            //- Write to Ostream
+            //void write(Ostream& os) const;
+
   // Member operators
 
         inline void operator+=(const mixturePengRobinson&);

src/thermophysicalModels/specie/equationOfState/cubicEquationOfState/mixtures/mixturePengRobinson/mixturePengRobinsonI.H

@@ -51,7 +51,9 @@ inline mixturePengRobinson::mixturePengRobinson
     scalar b,
     scalar Tcrit,
     scalar n,
-    scalar rhostd
+    scalar rhostd,
+    scalar rhoMin,
+    scalar rhoMax
 )
 : 
     pengRobinson(pr),
@@ -65,6 +67,8 @@ inline mixturePengRobinson::mixturePengRobinson
     Tcrit_=Tcrit; 
     n_=n;
     rhostd_=rhostd;
+    rhoMin_=rhoMin;
+    rhoMax_=rhoMax;
 }
 
 // Construct from components
@@ -74,7 +78,9 @@ inline mixturePengRobinson::mixturePengRobinson
     const pengRobinson& pr,
     label numOfComp,
     DynamicList<scalar> weigths,
-    DynamicList<mixturePengRobinson*> mixtureComponents
+    DynamicList<mixturePengRobinson*> mixtureComponents,
+    scalar rhoMin,
+    scalar rhoMax
 )
 : 
     pengRobinson(pr),
@@ -85,11 +91,8 @@ inline mixturePengRobinson::mixturePengRobinson
     
 {
     TSave=0.0;
-    //CL: following three DynamicList have no size in the new object
+    rhoMin_=rhoMin;
-    //CL: resize DynamicList
+    rhoMax_=rhoMax;
-    aComponents.resize(numOfComp);
-    daComponents.resize(numOfComp);
-    d2aComponents.resize(numOfComp);
     rhostd_=this->rho(Pstd, Tstd, Pstd*this->W()/(Tstd*this->R()));
 }
 
@@ -125,49 +128,47 @@ inline void mixturePengRobinson::updateModelCoefficients(const  scalar T) const
         // Checking if the mixture coefficient were already calculated for this temperature
         if(TSave!=T)
         {
-            label i,j;
+            label i,j,k;
  
             aSave=0;
             daSave=0;
             d2aSave=0;
             b_=0;
 
-            //CL: filling vector a, dadT and d2adT2
-            for (i=0;i<numOfComp;i++)
-            {
-                aComponents[i]=mixtureComponents[i]->a(T);
-                daComponents[i]=mixtureComponents[i]->dadT(T);
-                d2aComponents[i]=mixtureComponents[i]->d2adT2(T);
-            }
-
             for (i=0;i<numOfComp;i++)
             { 
                 //CL: getting a(T), dadT(T) and d2adT2(T) for the mixture 
                 //CL: using van der waals mixing rule
-                //CL: TODO: Include possibility to use coupling cofficients k_ij (see paper reference)
                 for (j=0;j<numOfComp;j++)
                 {
-                    aSave=aSave+weigths[i]*weigths[j]*pow(aComponents[i]*aComponents[j],0.5)*this->W()*this->W();
+
-   
+                    //CL: use mixture correction cofficients k_ij (see paper reference)
-                    // first and second order temperature derivative of the van der waals mixing rule for a(T) 
+                    if(mixtureComponents[0]->realMixtureCorr_==true)
-                    if(i==j)
                     {
-                        daSave=daSave+weigths[i]*weigths[i]*daComponents[i]*this->W()*this->W();
+                        // first and second order temperature derivative of the van der waals mixing rule for a(T) 
-                        d2aSave=d2aSave+weigths[i]*weigths[i]*d2aComponents[i]*this->W()*this->W();
+                        if(i==j)
+                        {
+                            aSave=aSave+amix(T,i,j);
+                            daSave=daSave+dadTmix(T,i,j);
+                            d2aSave=d2aSave+d2adT2mix(T,i,j);
+                        }
+                        else 
+                        {
+                            //CL: gives the position of the correction factor in the vector realMixtureCorrCoef_
+                            k=i+j-1;
+
+                            aSave=aSave+amix(T,i,j)*(1-mixtureComponents[0]->realMixtureCorrCoef_[k]);
+                            daSave=daSave+dadTmix(T,i,j)*(1-mixtureComponents[0]->realMixtureCorrCoef_[k]);
+                            d2aSave=d2aSave+d2adT2mix(T,i,j)*(1-mixtureComponents[0]->realMixtureCorrCoef_[k]);  
+                        }    
+
                     }
-                    else 
+                    else
                     {
-                        daSave=daSave+0.5*weigths[i]*weigths[j]*this->W()*this->W()
+                        aSave=aSave+amix(T,i,j);
-                            *(pow(aComponents[i]/aComponents[j],0.5)*daComponents[j]+pow(aComponents[j]/aComponents[i],0.5)*daComponents[i]);
+                        daSave=daSave+dadTmix(T,i,j);
-                   
+                        d2aSave=d2aSave+d2adT2mix(T,i,j);
-                        d2aSave=d2aSave+0.5*weigths[i]*weigths[j]*this->W()*this->W()*
+                    }           
-                             (
-                                 pow(aComponents[i]/aComponents[j],0.5)*d2aComponents[i]+pow(aComponents[j]/aComponents[i],0.5)*d2aComponents[j]
-                                     +0.5*pow(aComponents[i]*aComponents[j],-0.5)*(pow(daComponents[i],2)+pow(daComponents[j],2))
-                                         -0.5*daComponents[i]*daComponents[j]*(pow(aComponents[i]/aComponents[j],0.5)/aComponents[j]
-                                             +pow(aComponents[j]/aComponents[i],0.5)/aComponents[i])
-                             );   
-                    }              
                 }
 
                 //CL: getting b for the mixture 
@@ -177,11 +178,57 @@ inline void mixturePengRobinson::updateModelCoefficients(const  scalar T) const
                 //CL: saving the temperature at which the mixture coefficients are valid
                 TSave=T;
             }
+ 
+            b2_=b_*b_;
+            b3_=pow(b_,3);
+            b4_=pow(b_,4);
+            b5_=pow(b_,5);
         }
     }
-
 }
 
+//CL: returns part of  "a" for the mixture
+//CL: Important: a is a function of T
+//CL: TODO: RECHECK
+inline scalar mixturePengRobinson::amix(const scalar T, const label i, const label j) const
+{
+
+    return weigths[i]*weigths[j]*this->W()*this->W()*
+               pow(mixtureComponents[i]->a(T)*mixtureComponents[j]->a(T),0.5);
+}
+
+//CL: returns part of the first derivative of "a" for the mixture
+//CL: Important: a is a function of T
+//CL: TODO: RECHECK
+inline scalar mixturePengRobinson::dadTmix(const scalar T, const label i, const label j) const
+{
+
+    return 0.5*weigths[i]*weigths[j]*this->W()*this->W()
+               *(pow(mixtureComponents[i]->a(T)/mixtureComponents[j]->a(T),0.5)*mixtureComponents[j]->dadT(T)
+                   +pow(mixtureComponents[j]->a(T)/mixtureComponents[i]->a(T),0.5)*mixtureComponents[i]->dadT(T));
+}
+
+//CL: returns part of the second derivative of "a" for the mixture
+//CL: Important: a is a function of T
+//CL: TODO: RECHECK
+inline scalar mixturePengRobinson::d2adT2mix(const scalar T, const label i, const label j) const
+{
+
+    return  0.5*weigths[i]*weigths[j]*this->W()*this->W()*
+            (
+                pow(mixtureComponents[i]->a(T)/mixtureComponents[j]->a(T),0.5)*mixtureComponents[j]->d2adT2(T)
+                    +pow(mixtureComponents[j]->a(T)/mixtureComponents[i]->a(T),0.5)*mixtureComponents[i]->d2adT2(T)
+                        +pow(mixtureComponents[i]->a(T)*mixtureComponents[j]->a(T),-0.5)
+                            *mixtureComponents[i]->dadT(T)*mixtureComponents[j]->dadT(T)
+                                -0.5*pow(mixtureComponents[i]->a(T)/mixtureComponents[j]->a(T),0.5)*1/mixtureComponents[j]->a(T)
+                                    *pow(mixtureComponents[j]->dadT(T),2)
+                                        -0.5*pow(mixtureComponents[j]->a(T)/mixtureComponents[i]->a(T),0.5)
+                                            *1/mixtureComponents[i]->a(T)*pow(mixtureComponents[i]->dadT(T),2)
+            );
+}
+
+
+
 
 //returns the pressure for a given density and temperature    
 inline scalar mixturePengRobinson::p(const scalar rho, const scalar T) const
@@ -387,14 +434,15 @@ inline void mixturePengRobinson::operator+=(const mixturePengRobinson& pr)
     //CL: increase number of Components by 1
     numOfComp=numOfComp+1;
 
-    //CL: resize DynamicList so that they have a defined size in the mixture object
-    aComponents.resize(numOfComp);
-    daComponents.resize(numOfComp);
-    d2aComponents.resize(numOfComp);
-
     //CL: setting TSave=0--> results that all coefficients (a, da, d2a) are calculated in updateModelCoefficients
     TSave=0.0;
     singleComponent=0;
+
+    //CL:setting rho boundaries
+    rhoMin_=min(rhoMin_,pr.rhoMin_);
+    rhoMax_=max(rhoMax_,pr.rhoMax_);
+
+    //CL: calculating rho @ std
     rhostd_=this->rho(Pstd, Tstd, Pstd*this->W()/(Tstd*this->R()));
 }
 
@@ -424,7 +472,7 @@ inline mixturePengRobinson operator+
     mixtureComponents[pr1.numOfComp]=pr2.mixtureComponents[0];
 
     return mixturePengRobinson(static_cast<const pengRobinson&>(pr1)+static_cast<const pengRobinson&>(pr2),
-        pr1.numOfComp+1, weigths, mixtureComponents);
+        pr1.numOfComp+1, weigths, mixtureComponents, min(pr1.rhoMin_,pr2.rhoMin_), max(pr1.rhoMax_,pr2.rhoMax_));
 }
 
 
@@ -440,7 +488,7 @@ inline mixturePengRobinson operator*
     weigths[pr.numOfComp-1]=s*pr.nMoles(); 
 
    return mixturePengRobinson(s*static_cast<const pengRobinson&>(pr), pr.numOfComp, weigths,
-       pr.mixtureComponents,pr.a0_,pr.b_,pr.Tcrit_,pr.n_,pr.rhostd_);
+       pr.mixtureComponents,pr.a0_,pr.b_,pr.Tcrit_,pr.n_,pr.rhostd_, pr.rhoMin_, pr.rhoMax_);
 }
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

src/thermophysicalModels/specie/equationOfState/cubicEquationOfState/mixtures/mixtureRedlichKwong/mixtureRedlichKwong.C

@@ -0,0 +1,154 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright held by original author
+     \\/     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 2 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, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+Description
+    Redlich Kwong equation of state for mixtures.
+
+Author
+Christian Lucas
+Institut für Thermodynamik
+Technische Universität Braunschweig 
+Germany
+
+\*---------------------------------------------------------------------------*/
+
+#include "mixtureRedlichKwong.H"
+#include "IOstreams.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+mixtureRedlichKwong::mixtureRedlichKwong(Istream& is)
+:
+    redlichKwong(is),
+    numOfComp(1),
+    //CL: no real gas mixture correction when stream constructor is used
+    realMixtureCorr_(false)
+{ 
+    //CL: set size of weigths and mixtureComponents to 10
+    //CL: when more mixture components are used
+    //CL: size of the DynamicLis increases automatically
+    weigths.setSize(10);
+    mixtureComponents.setSize(10);
+
+    //Save a pointer of this object in the mixtureComponents array
+    mixtureComponents[0]=this;
+    is.check("mixtureRedlichKwong::mixtureRedlichKwong(Istream& is)");
+}
+//CL: Constructed needed in OpenFOAM 2.x.x
+//CL: Code works fine, but compiling problem in OpenFOAM 1.6.ext
+//CL:  because specie has no constructor using dict
+/*
+mixtureRedlichKwong::mixtureRedlichKwong(const dictionary& dict)
+:
+    redlichKwong(dict),
+    numOfComp(1),
+    realMixtureCorr_(dict.subDict("equationOfState").lookupOrDefault("realMixtureCorrection",false))
+{
+    //CL: set size of weigths and mixtureComponents to 10
+    //CL: when more mixture components are used
+    //CL: size of the DynamicLis increases automatically
+    weigths.setSize(10);
+    mixtureComponents.setSize(10);
+
+    //Save a pointer of this object in the mixtureComponents array
+    mixtureComponents[0]=this;
+
+    if(realMixtureCorr_==true)
+    {
+        //CL: reads number of mixture components
+        //CL: this number is needed to calculate the number of correction coefficients
+        nCom_=dict.subDict("equationOfState").lookupOrDefault("numberOfComponents",1);
+
+            if (nCom_<1||nCom_==1)
+            {
+                FatalErrorIn
+                (
+                    "mixtureRedlichKwong::mixtureRedlichKwong(const dictionary& dict)"
+                )   << "mixture has only one component or number is not defined correctly, "
+                    << "recheck numberOfComponents in the  thermophysicalProperties dict of your case"
+                    << abort(FatalError);
+            }
+
+        label i;
+        label j;
+        label k=0;
+      
+ 
+        // CL: saves the number of mixtureCorrectionCoefficients needed for this mixture
+        // CL: need to be set to 1
+        label sizeOfVector_=1;
+
+        //CL: size of the vector depends on the number of components
+        //CL: determine the size of the vector    
+        for (i=3; i<=nCom_;i++)
+        {   
+            sizeOfVector_=sizeOfVector_+(i-1);
+        }
+ 
+        //CL: setting the size 
+        realMixtureCorrCoef_.setSize(sizeOfVector_);
+       
+        //CL: Reading the real mixture correction coefficients
+        //CL: Naming convention e.g for 3 component mixture:
+        //CL: mixtureCorrectionCoefficient_12, mixtureCorrectionCoefficient_13, mixtureCorrectionCoefficient_23
+        for(i=1;i<=nCom_-1;i++)
+        {
+            for(j=i+1;j<=nCom_;j++)
+            {
+                realMixtureCorrCoef_[k]=dict.subDict("equationOfState")
+                    .lookupOrDefault("mixtureCorrectionCoefficient_" + Foam::name(i) + Foam::name(j),0.0);
+                k++;
+            }
+        }
+    }
+}
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+void Foam::mixtureRedlichKwong::write(Ostream& os) const
+{
+    redlichKwong::write(os);
+
+}
+*/
+// * * * * * * * * * * * * * * * Ostream Operator  * * * * * * * * * * * * * //
+
+Ostream& operator<<(Ostream& os, const mixtureRedlichKwong& rk)
+{
+    os  << static_cast<const redlichKwong&>(rk)<< token::SPACE ;
+
+    os.check("Ostream& operator<<(Ostream& os, const mixtureRedlichKwong& st)");
+    return os;
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// ************************************************************************* //

src/thermophysicalModels/specie/equationOfState/cubicEquationOfState/mixtures/mixtureRedlichKwong/mixtureRedlichKwong.H

@@ -76,7 +76,8 @@ class mixtureRedlichKwong
     
     // Private functions
 
-    //CL: function updates the model coefficients (a,b) of the mixture 
+    //CL: function updates the model coefficients (a,b) of the mixture
+    //CL: this is the function with the mixing rule 
     inline void updateModelCoefficients() const; 
 
 protected:
@@ -91,6 +92,20 @@ protected:
     //CL: counts the number of components
     mutable label numOfComp;
 
+
+    //CL: Variables used in real gas mixture correction
+
+        //CL: If true, the real gas mixture correction is used
+        mutable bool  realMixtureCorr_;
+
+        //CL: number of mixture components, needed to calculate the mixture correction factors needed
+        //CL: do not mistake this variable with numOfComp,
+        //CL: numOfComp is a counter to counts the number of components while the mixture is constructed
+        mutable label nCom_;
+
+        //CL: stores real mixture correction coefficients
+        mutable DynamicList<scalar> realMixtureCorrCoef_;
+
 public:
 
     // Constructors
@@ -105,7 +120,9 @@ public:
             DynamicList<mixtureRedlichKwong*> mixtureComponents,
             scalar a, 
             scalar b,
-            scalar rhostd
+            scalar rhostd,
+            scalar rhoMin,
+            scalar rhoMax
         );
 
         //- Construct from components
@@ -115,12 +132,17 @@ public:
             const redlichKwong& rK,
             label numOfComp,
             DynamicList<scalar> weigths,
-            DynamicList<mixtureRedlichKwong*> mixtureComponents
+            DynamicList<mixtureRedlichKwong*> mixtureComponents,
+            scalar rhoMin,
+            scalar rhoMax
         );
 
         //- Construct from Istream
         mixtureRedlichKwong(Istream&);
 
+        //- Construct from dictionary
+        //mixtureRedlichKwong(const dictionary& dict);
+
         //- Construct as named copy
         inline mixtureRedlichKwong(const word& name, const mixtureRedlichKwong&);
 
@@ -130,6 +152,12 @@ public:
         // Selector from Istream
         inline static autoPtr<mixtureRedlichKwong> New(Istream& is);
 
+        // I-O
+
+            //- Write to Ostream
+            //void write(Ostream& os) const;
+
+
   // Member operators
 
         inline void operator+=(const mixtureRedlichKwong&);

src/thermophysicalModels/specie/equationOfState/cubicEquationOfState/mixtures/mixtureRedlichKwong/mixtureRedlichKwongI.H

@@ -49,7 +49,9 @@ inline mixtureRedlichKwong::mixtureRedlichKwong
     DynamicList<mixtureRedlichKwong*> mixtureComponents,
     scalar a,
     scalar b,
-    scalar rhostd
+    scalar rhostd,
+    scalar rhoMin,
+    scalar rhoMax
 )
 : 
     redlichKwong(rK),
@@ -59,7 +61,14 @@ inline mixtureRedlichKwong::mixtureRedlichKwong
 {
     a_=a;
     b_=b;
+    b2_=b*b;
+    b3_=pow(b,3);
+    b5_=pow(b,4);
     rhostd_=rhostd;
+
+    //CL:setting rho boundaries
+    rhoMin_=rhoMin;
+    rhoMax_=rhoMax;
 }
 
 // Construct from components
@@ -69,7 +78,9 @@ inline mixtureRedlichKwong::mixtureRedlichKwong
     const redlichKwong& rK,
     label numOfComp,
     DynamicList<scalar> weigths,
-    DynamicList<mixtureRedlichKwong*> mixtureComponents
+    DynamicList<mixtureRedlichKwong*> mixtureComponents,
+    scalar rhoMin,
+    scalar rhoMax
 )
 : 
     redlichKwong(rK),
@@ -79,6 +90,12 @@ inline mixtureRedlichKwong::mixtureRedlichKwong
 {
     //CL: update model coefficients
     updateModelCoefficients();
+
+    //CL:setting rho boundaries
+    rhoMin_=rhoMin;
+    rhoMax_=rhoMax;
+
+    //CL: calculating rho @ std
     rhostd_=this->rho(Pstd, Tstd, Pstd/(Tstd*this->R()));
 }
 
@@ -110,7 +127,7 @@ inline autoPtr<mixtureRedlichKwong> mixtureRedlichKwong::New(Istream& is)
 //CL: uses the van der waals mixing rule
 inline void mixtureRedlichKwong::updateModelCoefficients() const
 {  
-    label i,j;
+    label i,j,k;
       
     a_=0;
     b_=0;
@@ -119,16 +136,38 @@ inline void mixtureRedlichKwong::updateModelCoefficients() const
     {
         //CL: getting a for the mixture 
         //CL: using van der waals mixing rule
-        //CL: TODO: Include possibility to use coupling cofficients k_ij (see paper reference)
         for (j=0;j<numOfComp;j++)
         {
-            a_=a_+weigths[i]*weigths[j]*pow(mixtureComponents[i]->a()*mixtureComponents[j]->a(),0.5)*this->W()*this->W();
+            //CL: Important: a is not a function of T
+            if(mixtureComponents[0]->realMixtureCorr_==true)
+            {
+                if(j!=i)
+                {
+                    //CL: gives the position of the correction factor in the vector realMixtureCorrCoef_
+                    k=i+j-1;
+
+                    a_=a_+weigths[i]*weigths[j]*pow(mixtureComponents[i]->a()*mixtureComponents[j]->a(),0.5)
+                        *this->W()*this->W()*(1-mixtureComponents[0]->realMixtureCorrCoef_[k]);
+                }
+                else
+                {
+                    a_=a_+weigths[i]*weigths[j]*pow(mixtureComponents[i]->a()*mixtureComponents[j]->a(),0.5)*this->W()*this->W();
+                }
+            }
+            //CL: don't use mixture correction cofficients k_ij (see paper reference)
+            else
+            {
+                a_=a_+weigths[i]*weigths[j]*pow(mixtureComponents[i]->a()*mixtureComponents[j]->a(),0.5)*this->W()*this->W();    
+            } 
+            //CL: getting b for the mixture 
+            //CL: using van der waals mixing rule
+            b_=b_+weigths[i]*mixtureComponents[i]->b()*this->W();
         }
-  
-        //CL: getting b for the mixture 
-        //CL: using van der waals mixing rule
-        b_=b_+weigths[i]*mixtureComponents[i]->b()*this->W();
     }
+    
+    b2_=b_*b_;
+    b3_=pow(b_,3);
+    b5_=pow(b_,4);
 }
 
 // * * * * * * * * * * * * * * * Member Operators  * * * * * * * * * * * * * //
@@ -155,7 +194,13 @@ inline void mixtureRedlichKwong::operator+=(const mixtureRedlichKwong& rK)
 
     //CL: update model coefficients
     updateModelCoefficients();
-    rhostd_=this->rho(Pstd, Tstd, Pstd/(Tstd*this->R()));
+
+    //CL:setting rho boundaries
+    rhoMin_=min(rhoMin_,rK.rhoMin_);
+    rhoMax_=max(rhoMax_,rK.rhoMax_);
+
+    //CL: calculating rho @ std
+    rhostd_=this->rho(Pstd, Tstd, Pstd*this->W()/(Tstd*this->R()));
 }
 
 // * * * * * * * * * * * * * * * Friend Operators  * * * * * * * * * * * * * //
@@ -184,7 +229,7 @@ inline mixtureRedlichKwong operator+
     mixtureComponents[rK1.numOfComp]=rK2.mixtureComponents[0];
 
     return mixtureRedlichKwong(static_cast<const redlichKwong&>(rK1)+static_cast<const redlichKwong&>(rK2),
-        rK1.numOfComp+1, weigths, mixtureComponents);
+        rK1.numOfComp+1, weigths, mixtureComponents, min(rK1.rhoMin_,rK2.rhoMin_), max(rK1.rhoMax_,rK2.rhoMax_));
 }
 
 
@@ -198,7 +243,8 @@ inline mixtureRedlichKwong operator*
     DynamicList<scalar> weigths=rK.weigths;
     weigths[rK.numOfComp-1]=s*rK.nMoles(); 
 
-   return mixtureRedlichKwong(s*static_cast<const redlichKwong&>(rK), rK.numOfComp, weigths, rK.mixtureComponents,rK.a_,rK.b_,rK.rhostd_);
+   return mixtureRedlichKwong(s*static_cast<const redlichKwong&>(rK), rK.numOfComp, weigths, 
+              rK.mixtureComponents,rK.a_,rK.b_,rK.rhostd_, rK.rhoMin_, rK.rhoMax_);
 }
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

src/thermophysicalModels/specie/equationOfState/cubicEquationOfState/mixtures/mixtureSoaveRedlichKwong/mixtureSoaveRedlichKwong.C

@@ -0,0 +1,157 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright held by original author
+     \\/     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 2 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, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+Description
+    Soave Redlich Kwong equation of state for mixtures.
+
+Author
+Christian Lucas
+Institut für Thermodynamik
+Technische Universität Braunschweig 
+Germany
+
+\*---------------------------------------------------------------------------*/
+
+#include "mixtureSoaveRedlichKwong.H"
+#include "IOstreams.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+mixtureSoaveRedlichKwong::mixtureSoaveRedlichKwong(Istream& is)
+:
+    soaveRedlichKwong(is),
+    numOfComp(1),
+    singleComponent(1),
+    //CL: no real gas mixture correction when stream constructor is used
+    realMixtureCorr_(false)
+{ 
+    //CL: set size of weigths, mixtureComponents ... to 10,
+    //CL: when more mixture componentents are used
+    //CL: size of the DynamicLis increases automatically
+    weigths.setSize(10);
+    mixtureComponents.setSize(10);
+
+    //Save a pointer of this object in the mixtureComponents array
+    mixtureComponents[0]=this;
+    is.check("mixtureSoaveRedlichKwong::mixtureSoaveRedlichKwong(Istream& is)");
+}
+
+//CL: Constructed needed in OpenFOAM 2.x.x
+//CL: Code works fine, but compiling problem in OpenFOAM 1.6.ext
+//CL:  because specie has no constructor using dict
+/*
+mixtureSoaveRedlichKwong::mixtureSoaveRedlichKwong(const dictionary& dict)
+:
+    soaveRedlichKwong(dict),
+    numOfComp(1),
+    singleComponent(1),
+    realMixtureCorr_(dict.subDict("equationOfState").lookupOrDefault("realMixtureCorrection",false))
+{ 
+    //CL: set size of weigths, mixtureComponents ... to 10,
+    //CL: when more mixture componentents are used
+    //CL: size of the DynamicLis increases automatically
+    weigths.setSize(10);
+    mixtureComponents.setSize(10);
+
+    //Save a pointer of this object in the mixtureComponents array
+    mixtureComponents[0]=this;
+
+    if(realMixtureCorr_==true)
+    {
+        //CL: reads number of mixture components
+        //CL: this number is needed to calculate the number of correction coefficients
+        nCom_=dict.subDict("equationOfState").lookupOrDefault("numberOfComponents",1);
+
+            if (nCom_<1||nCom_==1)
+            {
+                FatalErrorIn
+                (
+                    "mixturePengRobinson::mixturePengRobinson(const dictionary& dict)"
+                )   << "mixture has only one component or number is not defined correctly, "
+                    << "recheck numberOfComponents in the  thermophysicalProperties dict of your case"
+                    << abort(FatalError);
+            }
+
+        label i;
+        label j;
+        label k=0;
+ 
+        // CL: saves the number of mixtureCorrectionCoefficients needed for this mixture
+        // CL: need to be set to 1
+        label sizeOfVector_=1;
+
+        //CL: size of the vector depends on the number of components
+        //CL: determine the size of the vector    
+        for (i=3; i<=nCom_;i++)
+        {   
+            sizeOfVector_=sizeOfVector_+(i-1);
+        }
+ 
+        //CL: setting the size 
+        realMixtureCorrCoef_.setSize(sizeOfVector_);
+       
+        //CL: Reading the real mixture correction coefficients
+        //CL: Naming convention e.g for 3 component mixture:
+        //CL: mixtureCorrectionCoefficient_12, mixtureCorrectionCoefficient_13, mixtureCorrectionCoefficient_23
+        for(i=1;i<=nCom_-1;i++)
+        {
+            for(j=i+1;j<=nCom_;j++)
+            {
+                realMixtureCorrCoef_[k]=dict.subDict("equationOfState")
+                    .lookupOrDefault("mixtureCorrectionCoefficient_" + Foam::name(i) + Foam::name(j),0.0);
+                k++;
+            }
+        }
+    }
+}
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+void Foam::mixtureSoaveRedlichKwong::write(Ostream& os) const
+{
+    soaveRedlichKwong::write(os);
+
+}
+*/
+// * * * * * * * * * * * * * * * Ostream Operator  * * * * * * * * * * * * * //
+
+Ostream& operator<<(Ostream& os, const mixtureSoaveRedlichKwong& srk)
+{
+    os  << static_cast<const soaveRedlichKwong&>(srk)<< token::SPACE ;
+
+    os.check("Ostream& operator<<(Ostream& os, const mixtureSoaveRedlichKwong& st)");
+    return os;
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// ************************************************************************* //

src/thermophysicalModels/specie/equationOfState/cubicEquationOfState/mixtures/mixtureSoaveRedlichKwong/mixtureSoaveRedlichKwong.H

@@ -83,15 +83,6 @@ protected:
     //CL: needs to be multiplied by this->W() to get the molar fractions
     mutable DynamicList<scalar> weigths;
     
-    //CL: stores the values of a(T) of all componentes
-    mutable DynamicList<scalar> aComponents;
-
-    //CL: stores the values of the derivatives of da/dT of all componentes
-    mutable DynamicList<scalar> daComponents;
-
-    //CL: stores the values of the second order derivatives of d2a/dT2 of all componentes
-    mutable DynamicList<scalar> d2aComponents;
-    
     //CL: saves a pointer to the pure component classes of the mixture
     mutable DynamicList<mixtureSoaveRedlichKwong*> mixtureComponents;
    
@@ -103,6 +94,20 @@ protected:
     //CL: function updates the coefficients (aSave, daSave, d2aSave, b) of the mixture 
     inline void updateModelCoefficients(const  scalar T) const; 
 
+    //CL: Variables used in real gas mixture correction
+
+        //CL: If true, the real gas mixture correction is used
+        mutable bool  realMixtureCorr_;
+ 
+        //CL: number of mixture components, needed to calculate the mixture correction factors needed
+        //CL: do not mistake this variable with numOfComp,
+        //CL: numOfComp is a counter to counts the number of components while the mixture is constructed
+        mutable label nCom_;
+ 
+        //CL: stores real mixture correction coefficients
+        mutable DynamicList<scalar> realMixtureCorrCoef_;
+
+
 public:
 
     // Constructors
@@ -119,7 +124,9 @@ public:
             scalar b,
             scalar Tcrit,
             scalar n,
-            scalar rhostd
+            scalar rhostd,
+            scalar rhoMin,
+            scalar rhoMax
         );
 
         //- Construct from components
@@ -129,12 +136,17 @@ public:
             const soaveRedlichKwong& srk,
             label numOfComp,
             DynamicList<scalar> weigths,
-            DynamicList<mixtureSoaveRedlichKwong*> mixtureComponents
+            DynamicList<mixtureSoaveRedlichKwong*> mixtureComponents,
+            scalar rhoMin,
+            scalar rhoMax
         );
 
         //- Construct from Istream
         mixtureSoaveRedlichKwong(Istream&);
 
+        //- Construct from dictionary
+        //mixtureSoaveRedlichKwong(const dictionary& dict);
+
         //- Construct as named copy
         inline mixtureSoaveRedlichKwong(const word& name, const mixtureSoaveRedlichKwong&);
 
@@ -146,6 +158,13 @@ public:
 
    // Member functions
 
+        //CL: mixture coefficience
+        inline scalar amix(const scalar T, const label i, const  label j) const;
+                 
+        inline scalar dadTmix(const scalar T, const label i, const  label j) const;
+
+        inline scalar d2adT2mix(const scalar T, const label i, const  label j) const;
+
         inline scalar p(const scalar rho, const scalar T) const;
 
         //first order derivatives 
@@ -209,6 +228,11 @@ public:
             const scalar rho0
         ) const;
 
+        // I-O
+
+            //- Write to Ostream
+            //void write(Ostream& os) const;
+
   // Member operators
 
         inline void operator+=(const mixtureSoaveRedlichKwong&);

src/thermophysicalModels/specie/equationOfState/cubicEquationOfState/mixtures/mixtureSoaveRedlichKwong/mixtureSoaveRedlichKwongI.H

@@ -51,7 +51,9 @@ inline mixtureSoaveRedlichKwong::mixtureSoaveRedlichKwong
     scalar b,
     scalar Tcrit,
     scalar n,
-    scalar rhostd
+    scalar rhostd,
+    scalar rhoMin,
+    scalar rhoMax
 )
 : 
     soaveRedlichKwong(srk),
@@ -65,6 +67,8 @@ inline mixtureSoaveRedlichKwong::mixtureSoaveRedlichKwong
     Tcrit_=Tcrit; 
     n_=n;
     rhostd_=rhostd;
+    rhoMin_=rhoMin;
+    rhoMax_=rhoMax;
 }
 
 
@@ -75,7 +79,9 @@ inline mixtureSoaveRedlichKwong::mixtureSoaveRedlichKwong
     const soaveRedlichKwong& srk,
     label numOfComp,
     DynamicList<scalar> weigths,
-    DynamicList<mixtureSoaveRedlichKwong*> mixtureComponents
+    DynamicList<mixtureSoaveRedlichKwong*> mixtureComponents,
+    scalar rhoMin,
+    scalar rhoMax
 )
 : 
     soaveRedlichKwong(srk),
@@ -86,11 +92,8 @@ inline mixtureSoaveRedlichKwong::mixtureSoaveRedlichKwong
     
 {
     TSave=0.0;
-    //CL: following three DynamicList have no size in the new object
+    rhoMin_=rhoMin;
-    //CL: resize DynamicList
+    rhoMax_=rhoMax;
-    aComponents.resize(numOfComp);
-    daComponents.resize(numOfComp);
-    d2aComponents.resize(numOfComp);
     rhostd_=this->rho(Pstd, Tstd, Pstd*this->W()/(Tstd*this->R()));
 }
 
@@ -128,49 +131,47 @@ inline void mixtureSoaveRedlichKwong::updateModelCoefficients(const  scalar T) c
         // Checking if the mixture coefficient were already calculated for this temperature
         if(TSave!=T)
         {
-            label i,j;
+            label i,j,k;
  
             aSave=0;
             daSave=0;
             d2aSave=0;
             b_=0;
 
-            //CL: filling vector a, dadT and d2adT2
-            for (i=0;i<numOfComp;i++)
-            {
-                aComponents[i]=mixtureComponents[i]->a(T);
-                daComponents[i]=mixtureComponents[i]->dadT(T);
-                d2aComponents[i]=mixtureComponents[i]->d2adT2(T);
-            }
-
             for (i=0;i<numOfComp;i++)
             { 
                 //CL: getting a(T), dadT(T) and d2adT2(T) for the mixture 
                 //CL: using van der waals mixing rule
-                //CL: TODO: Include possibility to use coupling cofficients k_ij (see paper reference)
                 for (j=0;j<numOfComp;j++)
                 {
-                    aSave=aSave+weigths[i]*weigths[j]*pow(aComponents[i]*aComponents[j],0.5)*this->W()*this->W();
+
-   
+                    //CL: use mixture correction cofficients k_ij (see paper reference)
-                    // first and second order temperature derivative of the van der waals mixing rule for a(T) 
+                    if(mixtureComponents[0]->realMixtureCorr_==true)
-                    if(i==j)
                     {
-                        daSave=daSave+weigths[i]*weigths[i]*daComponents[i]*this->W()*this->W();
+                        // first and second order temperature derivative of the van der waals mixing rule for a(T) 
-                        d2aSave=d2aSave+weigths[i]*weigths[i]*d2aComponents[i]*this->W()*this->W();
+                        if(i==j)
+                        {
+                            aSave=aSave+amix(T,i,j);
+                            daSave=daSave+dadTmix(T,i,j);
+                            d2aSave=d2aSave+d2adT2mix(T,i,j);
+                        }
+                        else 
+                        {
+                            //CL: gives the position of the correction factor in the vector realMixtureCorrCoef_
+                            k=i+j-1;
+
+                            aSave=aSave+amix(T,i,j)*(1-mixtureComponents[0]->realMixtureCorrCoef_[k]);
+                            daSave=daSave+dadTmix(T,i,j)*(1-mixtureComponents[0]->realMixtureCorrCoef_[k]);
+                            d2aSave=d2aSave+d2adT2mix(T,i,j)*(1-mixtureComponents[0]->realMixtureCorrCoef_[k]);  
+                        }    
+
                     }
-                    else 
+                    else
                     {
-                        daSave=daSave+0.5*weigths[i]*weigths[j]*this->W()*this->W()
+                        aSave=aSave+amix(T,i,j);
-                            *(pow(aComponents[i]/aComponents[j],0.5)*daComponents[j]+pow(aComponents[j]/aComponents[i],0.5)*daComponents[i]);
+                        daSave=daSave+dadTmix(T,i,j);
-                   
+                        d2aSave=d2aSave+d2adT2mix(T,i,j);
-                        d2aSave=d2aSave+0.5*weigths[i]*weigths[j]*this->W()*this->W()*
+                    }           
-                             (
-                                 pow(aComponents[i]/aComponents[j],0.5)*d2aComponents[i]+pow(aComponents[j]/aComponents[i],0.5)*d2aComponents[j]
-                                     +0.5*pow(aComponents[i]*aComponents[j],-0.5)*(pow(daComponents[i],2)+pow(daComponents[j],2))
-                                         -0.5*daComponents[i]*daComponents[j]*(pow(aComponents[i]/aComponents[j],0.5)/aComponents[j]
-                                             +pow(aComponents[j]/aComponents[i],0.5)/aComponents[i])
-                             );   
-                    }              
                 }
 
                 //CL: getting b for the mixture 
@@ -180,10 +181,53 @@ inline void mixtureSoaveRedlichKwong::updateModelCoefficients(const  scalar T) c
                 //CL: saving the temperature at which the mixture coefficients are valid
                 TSave=T;
             }
+ 
+            b2_=b_*b_;
+            b3_=pow(b_,3);
+            b5_=pow(b_,5);
         }
     }
 }
 
+//CL: returns part of  "a" for the mixture
+//CL: Important: a is a function of T
+//CL: TODO: RECHECK
+inline scalar mixtureSoaveRedlichKwong::amix(const scalar T, const label i, const label j) const
+{
+
+    return weigths[i]*weigths[j]*this->W()*this->W()*
+               pow(mixtureComponents[i]->a(T)*mixtureComponents[j]->a(T),0.5);
+}
+
+//CL: returns part of the first derivative of "a" for the mixture
+//CL: Important: a is a function of T
+//CL: TODO: RECHECK
+inline scalar mixtureSoaveRedlichKwong::dadTmix(const scalar T, const label i, const label j) const
+{
+
+    return 0.5*weigths[i]*weigths[j]*this->W()*this->W()
+               *(pow(mixtureComponents[i]->a(T)/mixtureComponents[j]->a(T),0.5)*mixtureComponents[j]->dadT(T)
+                   +pow(mixtureComponents[j]->a(T)/mixtureComponents[i]->a(T),0.5)*mixtureComponents[i]->dadT(T));
+}
+
+//CL: returns part of the second derivative of "a" for the mixture
+//CL: Important: a is a function of T
+//CL: TODO: RECHECK
+inline scalar mixtureSoaveRedlichKwong::d2adT2mix(const scalar T, const label i, const label j) const
+{
+
+    return  0.5*weigths[i]*weigths[j]*this->W()*this->W()*
+            (
+                pow(mixtureComponents[i]->a(T)/mixtureComponents[j]->a(T),0.5)*mixtureComponents[j]->d2adT2(T)
+                    +pow(mixtureComponents[j]->a(T)/mixtureComponents[i]->a(T),0.5)*mixtureComponents[i]->d2adT2(T)
+                        +pow(mixtureComponents[i]->a(T)*mixtureComponents[j]->a(T),-0.5)
+                            *mixtureComponents[i]->dadT(T)*mixtureComponents[j]->dadT(T)
+                                -0.5*pow(mixtureComponents[i]->a(T)/mixtureComponents[j]->a(T),0.5)*1/mixtureComponents[j]->a(T)
+                                    *pow(mixtureComponents[j]->dadT(T),2)
+                                        -0.5*pow(mixtureComponents[j]->a(T)/mixtureComponents[i]->a(T),0.5)
+                                            *1/mixtureComponents[i]->a(T)*pow(mixtureComponents[i]->dadT(T),2)
+            );
+}
 
 //returns the pressure for a given density and temperature    
 inline scalar mixtureSoaveRedlichKwong::p(const scalar rho, const scalar T) const
@@ -390,14 +434,15 @@ inline void mixtureSoaveRedlichKwong::operator+=(const mixtureSoaveRedlichKwong&
     //CL: increase number of Components by 1
     numOfComp=numOfComp+1;
 
-    //CL: resize DynamicList so that they have a defined size in the mixture object
-    aComponents.resize(numOfComp);
-    daComponents.resize(numOfComp);
-    d2aComponents.resize(numOfComp);
-
     //CL: setting TSave=0--> results that all coefficients (a, da, d2a) are calculated in updateModelCoefficients
     TSave=0.0;
     singleComponent=0;
+
+    //CL:setting rho boundaries
+    rhoMin_=min(rhoMin_,srk.rhoMin_);
+    rhoMax_=max(rhoMax_,srk.rhoMax_);
+
+    //CL: calculating rho @ std
     rhostd_=this->rho(Pstd, Tstd, Pstd*this->W()/(Tstd*this->R()));
 }
 
@@ -428,7 +473,7 @@ inline mixtureSoaveRedlichKwong operator+
     mixtureComponents[srk1.numOfComp]=srk2.mixtureComponents[0];
 
     return mixtureSoaveRedlichKwong(static_cast<const soaveRedlichKwong&>(srk1)+static_cast<const soaveRedlichKwong&>(srk2),
-        srk1.numOfComp+1, weigths, mixtureComponents);
+        srk1.numOfComp+1, weigths, mixtureComponents, min(srk1.rhoMin_,srk2.rhoMin_), max(srk1.rhoMax_,srk2.rhoMax_));
 }
 
 
@@ -444,7 +489,7 @@ inline mixtureSoaveRedlichKwong operator*
     weigths[srk.numOfComp-1]=s*srk.nMoles(); 
 
    return mixtureSoaveRedlichKwong(s*static_cast<const soaveRedlichKwong&>(srk), srk.numOfComp, weigths,
-       srk.mixtureComponents,srk.a0_,srk.b_,srk.Tcrit_,srk.n_,srk.rhostd_);
+       srk.mixtureComponents,srk.a0_,srk.b_,srk.Tcrit_,srk.n_,srk.rhostd_, srk.rhoMin_, srk.rhoMax_);
 }
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

src/thermophysicalModels/specie/equationOfState/cubicEquationOfState/pengRobinson/pengRobinson.C

@@ -0,0 +1,125 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright held by original author
+     \\/     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 2 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, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+Description
+    Peng Robinson equation of state.
+
+Author
+Christian Lucas
+Institut für Thermodynamik
+Technische Universität Braunschweig 
+Germany
+
+\*---------------------------------------------------------------------------*/
+
+#include "pengRobinson.H"
+#include "IOstreams.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+pengRobinson::pengRobinson(Istream& is)
+:
+    specie(is),
+    pcrit_(readScalar(is)),
+    Tcrit_(readScalar(is)),
+    azentricFactor_(readScalar(is)),
+    a0_(0.457235*pow(this->RR,2)*pow(Tcrit_,2)/pcrit_),
+    b_(0.077796*this->RR*Tcrit_/pcrit_), 
+    n_(0.37464+1.54226*azentricFactor_-0.26992*pow(azentricFactor_,2)),
+    TSave(0.0), 
+    rhostd_(this->rho(Pstd,Tstd,Pstd/(Tstd*this->R()))),
+    //CL: Only uses the default values
+    rhoMin_(1e-3),
+    rhoMax_(1500),
+    b2_(pow(b_,2)),
+    b3_(pow(b_,3)),
+    b4_(pow(b_,4)),
+    b5_(pow(b_,5))  	
+{
+    is.check("pengRobinson::pengRobinson(Istream& is)");
+}
+//CL: Constructed needed in OpenFOAM 2.x.x
+//CL: Code works fine, but compiling problem in OpenFOAM 1.6.ext
+//CL:  because specie has no constructor using dict
+/*
+pengRobinson::pengRobinson(const dictionary& dict)
+:
+    specie(dict),
+    pcrit_(readScalar(dict.subDict("equationOfState").lookup("pCritical"))),
+    Tcrit_(readScalar(dict.subDict("equationOfState").lookup("TCritical"))),
+    azentricFactor_(readScalar(dict.subDict("equationOfState").lookup("azentricFactor"))),
+    //CL: rhoMin and rhoMax are only used as boundaries for the bisection methode (see rho function)
+    //CL: important: rhoMin and rhoMax are not used as boundary for the newton solver
+    //CL: therefore, rho can be larger than rhoMax and smaller than rhoMin
+    rhoMin_(dict.subDict("equationOfState").lookupOrDefault("rhoMin",1e-3)),
+    rhoMax_(dict.subDict("equationOfState").lookupOrDefault("rhoMax",1500)),
+    a0_(0.457235*pow(this->RR,2)*pow(Tcrit_,2)/pcrit_),
+    b_(0.077796*this->RR*Tcrit_/pcrit_), 
+    n_(0.37464+1.54226*azentricFactor_-0.26992*pow(azentricFactor_,2)),
+    TSave(0.0), 
+    b2_(pow(b_,2)),
+    b3_(pow(b_,3)),
+    b4_(pow(b_,4)),
+    b5_(pow(b_,5)),
+    rhostd_(this->rho(Pstd,Tstd,Pstd/(Tstd*this->R()))) 
+{}
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+void Foam::pengRobinson::write(Ostream& os) const
+{
+    specie::write(os);
+
+    dictionary dict("equationOfState");
+    dict.add("pCritical", pcrit_);
+    dict.add("TCritical", Tcrit_);
+    dict.add("azentricFactor", azentricFactor_);
+    dict.add("rhoMin", rhoMin_);
+    dict.add("rhoMax", rhoMax_);
+
+    os  << indent << dict.dictName() << dict;
+}
+*/
+// * * * * * * * * * * * * * * * Ostream Operator  * * * * * * * * * * * * * //
+
+Ostream& operator<<(Ostream& os, const pengRobinson& pr)
+{
+    os  << static_cast<const specie&>(pr)<< token::SPACE 
+        << pr.pcrit_ << tab<< pr.Tcrit_<< tab << pr.azentricFactor_;
+
+    os.check("Ostream& operator<<(Ostream& os, const pengRobinson& st)");
+    return os;
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// ************************************************************************* //

src/thermophysicalModels/specie/equationOfState/cubicEquationOfState/pengRobinson/pengRobinson.H

@@ -68,12 +68,6 @@ class pengRobinson
     public specie
 
 {
-    // Static Private data
-
-    //should be read from the fvSolution file where rhoMax and rhoMin values must be defined (for rhoSimpleFoam)
-    //HR: Don't know, yet. Let's make these static for starters
-    static const scalar rhoMax_;
-    static const scalar rhoMin_;
 
 protected:
     // Protected data
@@ -86,6 +80,15 @@ protected:
     scalar a0_;
     mutable scalar b_;
 
+    //CL: pow of constants b_ used in the code e.g. b2_=b*b;
+    mutable scalar b2_;
+    mutable scalar b3_;
+    mutable scalar b4_;
+    mutable scalar b5_;
+
+    //CL: rhoMin and rhoMax are only used as boundaries for the bisection methode (see rho function)
+    scalar rhoMax_;
+    scalar rhoMin_;
 
     //CL: Variables to save the values of a, dadT and d2adT2 of the mixture
     //CL: Variables must corrected for changing temperatures
@@ -99,7 +102,6 @@ protected:
     //- Density @STD, initialise after a0, b!
     mutable scalar rhostd_;
 
-
     //Protected functions
     //CL: function updates the coefficients (aSave, daSave, d2aSave)
     inline void updateModelCoefficients(const  scalar T) const; 
@@ -118,6 +120,9 @@ public:
         //- Construct from Istream
         pengRobinson(Istream&);
 
+        //- Construct from dictionary
+        //pengRobinson(const dictionary& dict);
+
         //- Construct as named copy
         inline pengRobinson(const word& name,const pengRobinson&);
 
@@ -139,12 +144,22 @@ public:
         //CL: second order temperature deriviative of model coefficient a(T)
 	inline scalar d2adT2(const  scalar T)const;
 
+        //Return Peng Robinson factors
 	inline scalar a0()const;
 
 	inline scalar b()const;
 
 	inline scalar n()const;
 
+        //CL: return power of constants b_
+	inline scalar b2()const;
+
+	inline scalar b3()const;
+
+	inline scalar b4()const;
+
+	inline scalar b5()const;
+
         //CL: Equation of state
         inline scalar p(const scalar rho, const scalar T) const;
 
@@ -214,6 +229,12 @@ public:
         ) const;
 
 
+        // I-O
+
+            //- Write to Ostream
+            //void write(Ostream& os) const;
+
+
     // Member operators
 
         inline void operator+=(const pengRobinson&);

src/thermophysicalModels/specie/equationOfState/cubicEquationOfState/pengRobinson/pengRobinsonI.H

@@ -61,7 +61,11 @@ inline pengRobinson::pengRobinson(const word& name, const pengRobinson& pr)
     b_(pr.b_),
     n_(pr.n_),
     rhostd_(pr.rhostd_),
-    TSave(0) 
+    TSave(0),
+    b2_(pr.b2_),
+    b3_(pr.b3_),
+    b4_(pr.b4_),
+    b5_(pr.b5_)
 {}
 
 
@@ -166,13 +170,33 @@ inline scalar pengRobinson::n()const
     return n_;
 }
 
+//CL: pow of constant b() used in the code e.g. b2_=b*b;
+inline scalar pengRobinson::b2()const
+{
+    return b2_;
+}
+
+inline scalar pengRobinson::b3()const
+{
+    return b3_;
+}
+
+inline scalar pengRobinson::b4()const
+{
+    return b4_;
+}
+
+inline scalar pengRobinson::b5()const
+{
+    return b5_;
+}
 
 //returns the pressure for a given density and temperature    
 inline scalar pengRobinson::p(const scalar rho,const scalar T) const
 {
     scalar Vm = this->W()/rho;
      
-    return this->RR*T/(Vm-b_)-a(T)/(pow(Vm,2)+2*b_*Vm-pow(b_,2));
+    return this->RR*T/(Vm-b())-a(T)/(pow(Vm,2)+2*b()*Vm-b2());
 }
 
 
@@ -181,19 +205,21 @@ inline scalar pengRobinson::p(const scalar rho,const scalar T) const
 inline scalar pengRobinson::dpdv(const scalar rho, const scalar T) const
 {
     scalar Vm = this->W()/rho;
+    scalar Vm2 = Vm*Vm;
+    scalar Vm3 = Vm*Vm*Vm;
 
     return(
               2*a(T)*
               (
-                  pow(b_,3)-pow(b_,2)*Vm-b_*pow(Vm,2)+pow(Vm,3)
+                  b3()-b2()*Vm-b()*Vm2+Vm3
               )
               -this->RR*T*
               (
-                  pow(b_,4)-4*pow(b_,3)*Vm+2*pow(b_,2)*pow(Vm,2)
+                  b4()-4*b3()*Vm+2*b2()*Vm2
-                      +4*b_*pow(Vm,3)+pow(Vm,4)
+                      +4*b()*Vm3+pow(Vm,4)
               )
           )
-          /(pow(b_-Vm,2)*pow(pow(b_,2)-2*b_*Vm-pow(Vm,2),2));
+          /(pow(b()-Vm,2)*pow(b2()-2*b()*Vm-Vm2,2));
 }
 
 
@@ -203,7 +229,7 @@ inline scalar pengRobinson::dpdT(const scalar rho, const scalar T) const
 {
     scalar Vm = this->W()/rho;
 
-    return this->RR/(Vm-b_)-dadT(T)/(pow(Vm,2)+2*b_*Vm-pow(b_,2));
+    return this->RR/(Vm-b())-dadT(T)/(pow(Vm,2)+2*b()*Vm-b2());
 }
 
 
@@ -229,7 +255,8 @@ inline scalar pengRobinson::integral_p_dv(const scalar rho,const scalar T) const
 {
     scalar Vm = this->W()/rho;
 
-    return -pow(2,0.5)*a(T)*log(b_*(1-pow(2,0.5))+Vm)/(4*b_)+this->RR*T*log(Vm-b_)+pow(2,0.5)*a(T)*log(b_*(pow(2,0.5)+1)+Vm)/(4*b_);
+    return -pow(2,0.5)*a(T)*log(b()*(1-pow(2,0.5))+Vm)/(4*b())+this->RR*T*log(Vm-b())
+               +pow(2,0.5)*a(T)*log(b()*(pow(2,0.5)+1)+Vm)/(4*b());
 }
 
 
@@ -239,7 +266,8 @@ inline scalar pengRobinson::integral_dpdT_dv(const scalar rho,const scalar T) co
 {
     scalar Vm = this->W()/rho;
 
-    return -pow(2,0.5)*dadT(T)*log(b_*(1-pow(2,0.5))+Vm)/(4*b_)+this->RR*log(Vm-b_)+pow(2,0.5)*dadT(T)*log(b_*(pow(2,0.5)+1)+Vm)/(4*b_);
+    return -pow(2,0.5)*dadT(T)*log(b()*(1-pow(2,0.5))+Vm)/(4*b())
+               +this->RR*log(Vm-b())+pow(2,0.5)*dadT(T)*log(b()*(pow(2,0.5)+1)+Vm)/(4*b());
 }
 
 
@@ -248,7 +276,7 @@ inline scalar pengRobinson::d2pdT2(const scalar rho,const scalar T) const
 {
     scalar Vm = this->W()/rho;
 
-    return -d2adT2(T)/(pow(Vm,2)+2*b_*Vm-pow(b_,2));
+    return -d2adT2(T)/(pow(Vm,2)+2*b()*Vm-b2());
 }
 
 
@@ -256,20 +284,24 @@ inline scalar pengRobinson::d2pdT2(const scalar rho,const scalar T) const
 inline scalar pengRobinson::d2pdv2(const scalar rho,const scalar T) const
 {
     scalar Vm = this->W()/rho;
+    scalar Vm2 = Vm*Vm;
+    scalar Vm3 = Vm*Vm*Vm;
+    scalar Vm4 = Vm*Vm*Vm*Vm;
+    scalar Vm5 = Vm*Vm*Vm*Vm*Vm;
 
     return 2*
            (
                a(T)*
                (
-                   5*pow(b_,5)-9*pow(b_,4)*Vm+4*pow(b_,2)*pow(Vm,3)+3*b_*pow(Vm,4)-3*pow(Vm,5)
+                   5*b5()-9*b4()*Vm+4*b2()*Vm3+3*b()*Vm4-3*Vm5
                )
                -this->RR*T*
                (
-                   pow(b_,6)-6*pow(b_,5)*Vm+9*pow(b_,4)*pow(Vm,2)+4*pow(b_,3)*pow(Vm,3)
+                   pow(b(),6)-6*b5()*Vm+9*b4()*Vm2+4*b3()*Vm3
-                       -9*pow(b_,2)*pow(Vm,4)-6*b_*pow(Vm,5)-pow(Vm,6)
+                       -9*b2()*Vm4-6*b()*Vm5-pow(Vm,6)
                )
            )
-           /(pow(b_-Vm,3)*pow(pow(b_,2)-2*b_*Vm-pow(Vm,2),3));
+           /(pow(b()-Vm,3)*pow(b2()-2*b()*Vm-Vm2,3));
 }
 
 
@@ -291,19 +323,21 @@ inline scalar pengRobinson::d2vdT2(const scalar rho, const scalar T) const
 inline scalar pengRobinson::d2pdvdT(const scalar rho, const scalar T) const
 {
     scalar Vm = this->W()/rho;
+    scalar Vm2 = Vm*Vm;
+    scalar Vm3 = Vm*Vm*Vm;
   
     return(
               2*dadT(T)*
               (
-                  pow(b_,3)-pow(b_,2)*Vm-b_*pow(Vm,2)+pow(Vm,3)
+                  b3()-b2()*Vm-b()*Vm2+Vm3
               )
               -this->RR*
               (
-                  pow(b_,4)-4*pow(b_,3)*Vm+2*pow(b_,2)*pow(Vm,2)
+                  b4()-4*b3()*Vm+2*b2()*Vm2
-                      +4*b_*pow(Vm,3)+pow(Vm,4)
+                      +4*b()*Vm3+pow(Vm,4)
               )
           )
-          /(pow(b_-Vm,2)*pow(pow(b_,2)-2*b_*Vm-pow(Vm,2),2));
+          /(pow(b()-Vm,2)*pow(b2()-2*b()*Vm-Vm2,2));
 }
 
 
@@ -313,7 +347,7 @@ inline scalar pengRobinson::integral_d2pdT2_dv(const scalar rho,const scalar T)
 {
     scalar Vm = this->W()/rho;
 
-    return pow(2,0.5)*d2adT2(T)*log(b_*(pow(2,0.5)+1)+Vm)/(4*b_)-pow(2,0.5)*d2adT2(T)*log(b_*(1-pow(2,0.5))+Vm)/(4*b_);
+    return pow(2,0.5)*d2adT2(T)*log(b()*(pow(2,0.5)+1)+Vm)/(4*b())-pow(2,0.5)*d2adT2(T)*log(b()*(1-pow(2,0.5))+Vm)/(4*b());
 }
 
 
@@ -342,7 +376,7 @@ inline scalar pengRobinson::rho(
     const scalar rho0
 ) const
 {
- 
+
    scalar molarVolumePrevIteration;
    scalar molarVolume;
    label iter=0;
@@ -415,7 +449,7 @@ inline scalar pengRobinson::rho(
        {
            FatalErrorIn
            (
-               "inline scalar redlichKwong::rho(const scalar p, const scalar T, const scalar rho0) const "
+               "inline scalar pengRobinson::rho(const scalar p, const scalar T, const scalar rho0) const "
            )   << "Maximum number of iterations exceeded"
                << abort(FatalError);
        }

src/thermophysicalModels/specie/equationOfState/cubicEquationOfState/redlichKwong/redlichKwong.C

@@ -41,18 +41,6 @@ Germany
 namespace Foam
 {
 
-/* * * * * * * * * * * * * * * Private static data * * * * * * * * * * * * * */
-
-const scalar redlichKwong::rhoMin_
-(
-    debug::tolerances("redlichKwongRhoMin", 1e-3)
-);
-
-const scalar redlichKwong::rhoMax_
-(
-    debug::tolerances("redlichKwongRhoMax", 1500)
-);
-
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
 redlichKwong::redlichKwong(Istream& is)
@@ -62,25 +50,67 @@ redlichKwong::redlichKwong(Istream& is)
     Tcrit_(readScalar(is)),
     a_(0.42748*pow(this->RR,2)*pow(Tcrit_,2.5)/pcrit_),
     b_(0.08664*this->RR*Tcrit_/pcrit_),
-        // Starting GUESS for the density by ideal gas law
+    //CL: Only uses the default values
-        rhostd_(this->rho(Pstd, Tstd, Pstd/(Tstd*this->R())))
+    rhoMin_(1e-3),
+    rhoMax_(1500),
+    b2_(pow(b_,2)),
+    b3_(pow(b_,3)),
+    b5_(pow(b_,5)),
+    // Starting GUESS for the density by ideal gas law
+    rhostd_(this->rho(Pstd, Tstd, Pstd/(Tstd*this->R())))
 { 
     is.check("redlichKwong::redlichKwong(Istream& is)");
 }
+//CL: Constructed needed in OpenFOAM 2.x.x
+//CL: Code works fine, but compiling problem in OpenFOAM 1.6.ext
+//CL:  because specie has no constructor using dict
+/*
+redlichKwong::redlichKwong(const dictionary& dict)
+:
+    specie(dict),
+    pcrit_(readScalar(dict.subDict("equationOfState").lookup("pCritical"))),
+    Tcrit_(readScalar(dict.subDict("equationOfState").lookup("TCritical"))),
+    //CL: rhoMin and rhoMax are only used as boundaries for the bisection methode (see rho function)
+    //CL: important: rhoMin and rhoMax are not used as boundary for the newton solver
+    //CL: therefore, rho can be larger than rhoMax and smaller than rhoMin
+    rhoMin_(dict.subDict("equationOfState").lookupOrDefault("rhoMin",1e-3)),
+    rhoMax_(dict.subDict("equationOfState").lookupOrDefault("rhoMax",1500)),
+    a_(0.42748*pow(this->RR,2)*pow(Tcrit_,2.5)/pcrit_),
+    b_(0.08664*this->RR*Tcrit_/pcrit_),
+    b2_(pow(b_,2)),
+    b3_(pow(b_,3)),
+    b5_(pow(b_,5)),
+    // Starting GUESS for the density by ideal gas law
+    rhostd_(this->rho(Pstd, Tstd, Pstd/(Tstd*this->R()))) 
+{}
 
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+void Foam::redlichKwong::write(Ostream& os) const
+{
+    specie::write(os);
+
+    dictionary dict("equationOfState");
+    dict.add("pCritical", pcrit_);
+    dict.add("TCritical", Tcrit_);
+    dict.add("rhoMin", rhoMin_);
+    dict.add("rhoMax", rhoMax_);
+
+    os  << indent << dict.dictName() << dict;
+}
+*/
 
 // * * * * * * * * * * * * * * * Ostream Operator  * * * * * * * * * * * * * //
 
 Ostream& operator<<(Ostream& os, const redlichKwong& rk)
 {
-    os  << static_cast<const specie&>(rk)<< tab
+    os  << static_cast<const specie&>(rk)<< token::SPACE 
         << rk.pcrit_ << tab<< rk.Tcrit_;
 
     os.check("Ostream& operator<<(Ostream& os, const redlichKwong& st)");
     return os;
 }
 
-
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 } // End namespace Foam

src/thermophysicalModels/specie/equationOfState/cubicEquationOfState/redlichKwong/redlichKwong.H

@@ -59,22 +59,27 @@ class redlichKwong
 :
     public specie
 {
-    // Static Private data
 
-    //should be read from the fvSolution file where rhoMax and rhoMin values must be defined (for rhoSimpleFoam)
-    //HR: Don't know, yet. Let's make these static for starters
-    static const scalar rhoMax_;
-    static const scalar rhoMin_;
 
 protected:
+
     // Protected data
     scalar pcrit_;
     scalar Tcrit_; 
 
+    //CL: rhoMin and rhoMax are only used as boundaries for the bisection methode (see rho function)
+    scalar rhoMax_;
+    scalar rhoMin_;
+
     //-Redlich Kwong factors
     mutable scalar a_;
     mutable scalar b_;
 
+    //CL: pow of constants b_ used in the code e.g. b2_=b*b;
+    mutable scalar b2_;
+    mutable scalar b3_;
+    mutable scalar b5_;
+
     //- Density @STD, initialise after a, b!
     mutable scalar rhostd_;
 
@@ -91,6 +96,9 @@ public:
         //- Construct from Istream
         redlichKwong(Istream&);
 
+        //- Construct from dictionary
+        //redlichKwong(const dictionary& dict);
+
         //- Construct as named copy
         inline redlichKwong(const word& name, const redlichKwong&);
 
@@ -104,12 +112,20 @@ public:
 
         inline scalar rhostd() const;
 
+        //Return Redlich Kwong factors
         inline scalar a() const;
 
         inline scalar b() const;
 
         inline scalar p(const scalar rho, const scalar T) const;
 
+        //CL: return power of constants b_
+	inline scalar b2()const;
+
+	inline scalar b3()const;
+
+	inline scalar b5()const;
+
 
         //first order derivatives 
         inline scalar dpdv(const scalar rho, const  scalar T) const;
@@ -173,6 +189,12 @@ public:
         ) const;
 
 
+        // I-O
+
+            //- Write to Ostream
+            //void write(Ostream& os) const;
+
+
     // Member operators
 
         inline void operator+=(const redlichKwong&);

src/thermophysicalModels/specie/equationOfState/cubicEquationOfState/redlichKwong/redlichKwongI.H

@@ -57,7 +57,10 @@ inline redlichKwong::redlichKwong(const word& name, const redlichKwong& rk)
     Tcrit_(rk.Tcrit_),
     a_(rk.a_),
     b_(rk.b_),
-    rhostd_(rk.rhostd_)
+    rhostd_(rk.rhostd_),
+    b2_(rk.b2_),
+    b3_(rk.b3_),
+    b5_(rk.b5_)
 {}
 
 
@@ -94,6 +97,22 @@ inline scalar redlichKwong::b()const
     return b_;
 }
 
+inline scalar redlichKwong::b2()const
+{
+    return b2_;
+}
+
+inline scalar redlichKwong::b3()const
+{
+    return b3_;
+}
+
+inline scalar redlichKwong::b5()const
+{
+    return b5_;
+}
+
+
 
 //returns the pressure for a given density and temperature    
 inline scalar redlichKwong::p(const scalar rho, const scalar T) const
@@ -108,9 +127,11 @@ inline scalar redlichKwong::p(const scalar rho, const scalar T) const
 inline scalar redlichKwong::dpdv(const scalar rho, const scalar T) const
 {
     scalar Vm = this->W()/rho;
-    return (a_*(pow(b_,3) - 3*b_*pow(Vm,2) + 2*pow(Vm,3))
+    scalar Vm2 = Vm*Vm;
-    - this->RR*pow(T,1.5)*pow(Vm,2)*(pow(b_,2) + 2*b_*Vm + pow(Vm,2)))
+
-        /(sqrt(T)*pow(Vm,2)*pow((b_ + Vm),2)*pow( (b_ - Vm),2));
+    return (a_*(b3() - 3*b_*Vm2 + 2*pow(Vm,3))
+    - this->RR*pow(T,1.5)*Vm2*(b2() + 2*b_*Vm + Vm2))
+        /(sqrt(T)*Vm2*pow((b_ + Vm),2)*pow( (b_ - Vm),2));
 }
 
 
@@ -182,19 +203,22 @@ inline scalar redlichKwong::d2pdT2(const scalar rho, const scalar T) const
 inline scalar redlichKwong::d2pdv2(const scalar rho, const scalar T) const
 {
     scalar Vm = this->W()/rho;
+    scalar Vm2 = Vm*Vm;
+    scalar Vm3 = Vm*Vm*Vm;
+
     return
     (
         2*(
             a_*(
-                pow(b_,5)-3*pow(b_,3)*pow(Vm,2)
+                b5()-3*b3()*Vm2
-                - pow(b_,2)*pow(Vm,3)
+                - b2()*Vm3
                 + 6*b_*pow(Vm,4)-3*pow(Vm,5)
             )
-            + this->RR*pow(T,1.5)*pow(Vm,3)*(pow(b_,3)
+            + this->RR*pow(T,1.5)*Vm3*(b3()
-            + 3*pow(b_,2)*Vm
+            + 3*b2()*Vm
-            + 3*b_*pow(Vm,2)+pow(Vm,3))
+            + 3*b_*Vm2+Vm3)
         )
-        /(sqrt(T)*pow(Vm,3)*pow((b_ + Vm),3)*pow(Vm-b_,3))
+        /(sqrt(T)*Vm3*pow((b_ + Vm),3)*pow(Vm-b_,3))
     );
 }
 
@@ -217,13 +241,14 @@ inline scalar redlichKwong::d2vdT2(const scalar rho, const scalar T) const
 inline scalar redlichKwong::d2pdvdT(const scalar rho, const scalar T) const
 {
     scalar Vm = this->W()/rho;
+    scalar Vm2 = Vm*Vm;
     
     return
     -(0.5*(
-        a_*(pow(b_,3) - 3*b_*pow(Vm,2) + 2*pow(Vm,3))
+        a_*(b3() - 3*b_*Vm2 + 2*pow(Vm,3))
-      + 2*this->RR*pow(T,1.5)*pow(Vm,2)*(pow(b_,2) + 2*b_*Vm + pow(Vm,2))
+      + 2*this->RR*pow(T,1.5)*Vm2*(b2() + 2*b_*Vm + Vm2)
     ))
-    /(pow(T,1.5)*pow(Vm,2)*pow(b_ + Vm,2)*pow(b_ - Vm,2));
+    /(pow(T,1.5)*Vm2*pow(b_ + Vm,2)*pow(b_ - Vm,2));
 }
 
 
@@ -266,7 +291,7 @@ inline scalar redlichKwong::isothermalCompressiblity
 }
 
 
-//- Return density [kg/m^3]on 
+//- Return density [kg/m^3]
 inline scalar redlichKwong::rho
 (
     const scalar p,
@@ -292,6 +317,7 @@ inline scalar redlichKwong::rho
        label i=0;
        do
        {
+           //CL: modified Newton solver
            molarVolume=molarVolumePrevIteration
                -(
                    (this->p((this->W()/molarVolumePrevIteration),T) - p)
@@ -302,8 +328,7 @@ inline scalar redlichKwong::rho
            if (i>8) 
            {
                //CL: using bisection methode as backup,
-               //CL: solution must be between rho=0.001 to rho=1500;
+               //CL: solution must be between rhoMin_ to rhoMax
-               //CL: if not, change rhoMax_ and rhoMin_
                for(i=0; i<200; i++)
                {
                    scalar f1 = this->p(rho1,T) - p;

src/thermophysicalModels/specie/equationOfState/cubicEquationOfState/soaveRedlichKwong/soaveRedlichKwong.C

@@ -41,17 +41,6 @@ Germany
 namespace Foam
 {
 
-/* * * * * * * * * * * * * * * Private static data * * * * * * * * * * * * * */
-const scalar soaveRedlichKwong::rhoMin_
-(
-    debug::tolerances("soaveRedlichKwongRhoMin", 1e-3)
-);
-
-const scalar soaveRedlichKwong::rhoMax_
-(
-    debug::tolerances("soaveRedlichKwongRhoMax", 1500)
-);
-
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
 soaveRedlichKwong::soaveRedlichKwong(Istream& is)
@@ -63,18 +52,65 @@ soaveRedlichKwong::soaveRedlichKwong(Istream& is)
     a0_(0.42747*pow(this->RR,2)*pow(Tcrit_,2)/(pcrit_)),
     b_(0.08664*this->RR*Tcrit_/pcrit_),
     n_(0.48508+1.55171*azentricFactor_-0.15613*pow(azentricFactor_,2)),
-        // Starting GUESS for the density by ideal gas law
+    TSave(0.0), 
-        rhostd_(this->rho(Pstd,Tstd,Pstd/(Tstd*this->R())))
+    //CL: Only uses the default values
+    rhoMin_(1e-3),
+    rhoMax_(1500),
+    b2_(pow(b_,2)),
+    b3_(pow(b_,3)),
+    b5_(pow(b_,5)),
+    // Starting GUESS for the density by ideal gas law
+    rhostd_(this->rho(Pstd,Tstd,Pstd/(Tstd*this->R())))
 {
     is.check("soaveRedlichKwong::soaveRedlichKwong(Istream& is)"); 	
 }
+//CL: Constructed needed in OpenFOAM 2.x.x
+//CL: Code works fine, but compiling problem in OpenFOAM 1.6.ext
+//CL:  because specie has no constructor using dict
+/*
+soaveRedlichKwong::soaveRedlichKwong(const dictionary& dict)
+:
+    specie(dict),
+    pcrit_(readScalar(dict.subDict("equationOfState").lookup("pCritical"))),
+    Tcrit_(readScalar(dict.subDict("equationOfState").lookup("TCritical"))),
+    azentricFactor_(readScalar(dict.subDict("equationOfState").lookup("azentricFactor"))),
+    //CL: rhoMin and rhoMax are only used as boundaries for the bisection methode (see rho function)
+    //CL: important: rhoMin and rhoMax are not used as boundary for the newton solver
+    //CL: therefore, rho can be larger than rhoMax and smaller than rhoMin
+    rhoMin_(dict.subDict("equationOfState").lookupOrDefault("rhoMin",1e-3)),
+    rhoMax_(dict.subDict("equationOfState").lookupOrDefault("rhoMax",1500)),
+    a0_(0.42747*pow(this->RR,2)*pow(Tcrit_,2)/(pcrit_)),
+    b_(0.08664*this->RR*Tcrit_/pcrit_),
+    n_(0.48508+1.55171*azentricFactor_-0.15613*pow(azentricFactor_,2)),
+    TSave(0.0), 
+    b2_(pow(b_,2)),
+    b3_(pow(b_,3)),
+    b5_(pow(b_,5)),
+    // Starting GUESS for the density by ideal gas law
+    rhostd_(this->rho(Pstd,Tstd,Pstd/(Tstd*this->R())))
+{}
 
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
+void Foam::soaveRedlichKwong::write(Ostream& os) const
+{
+    specie::write(os);
+
+    dictionary dict("equationOfState");
+    dict.add("pCritical", pcrit_);
+    dict.add("TCritical", Tcrit_);
+    dict.add("azentricFactor", azentricFactor_);
+    dict.add("rhoMin", rhoMin_);
+    dict.add("rhoMax", rhoMax_);
+
+    os  << indent << dict.dictName() << dict;
+}
+*/
 // * * * * * * * * * * * * * * * Ostream Operator  * * * * * * * * * * * * * //
 
 Ostream& operator<<(Ostream& os, const soaveRedlichKwong& srk)
 {
-    os  << static_cast<const specie&>(srk)<< tab
+    os  << static_cast<const specie&>(srk)<< token::SPACE 
         << srk.pcrit_ << tab<< srk.Tcrit_<<tab<<srk.azentricFactor_;
 
     os.check("Ostream& operator<<(Ostream& os, const soaveRedlichKwong& st)");

src/thermophysicalModels/specie/equationOfState/cubicEquationOfState/soaveRedlichKwong/soaveRedlichKwong.H

@@ -66,13 +66,6 @@ class soaveRedlichKwong
     public specie
 
 {
-    // Static Private data
-
-    //should be read from the fvSolution file where rhoMax and rhoMin values must be defined (for rhoSimpleFoam)
-    //HR: Don't know, yet. Let's make these static for starters
-    static const scalar rhoMax_;
-    static const scalar rhoMin_;
-
 
 protected:
     // Protected data
@@ -85,6 +78,14 @@ protected:
     scalar a0_;
     mutable scalar b_;
 
+    //CL: pow of constants b_ used in the code e.g. b2_=b*b;
+    mutable scalar b2_;
+    mutable scalar b3_;
+    mutable scalar b5_;
+
+    //CL: rhoMin and rhoMax are only used as boundaries for the bisection methode (see rho function)
+    scalar rhoMax_;
+    scalar rhoMin_;
 
     //CL: Variables to save the values of a, dadT and d2adT2 of the mixture
     //CL: Variables must corrected for changing temperatures
@@ -117,6 +118,9 @@ public:
         //- Construct from Istream
         soaveRedlichKwong(Istream&);
 
+        //- Construct from dictionary
+        //soaveRedlichKwong(const dictionary& dict);
+
         //- Construct as named copy
         inline soaveRedlichKwong(const word& name,const soaveRedlichKwong&);
 
@@ -138,12 +142,20 @@ public:
         //CL: second order temperature deriviative of model coefficient a(T)
 	inline scalar d2adT2(const  scalar T)const;
 
+        //Return Soave Redlich Kwong factors
 	inline scalar a0()const;
 
 	inline scalar b()const;
 
 	inline scalar n()const;
 
+        //CL: return power of constants b_
+	inline scalar b2()const;
+
+	inline scalar b3()const;
+
+	inline scalar b5()const;
+
         //CL: Equation of state
         inline scalar p(const scalar rho, const scalar T) const;
 
@@ -209,6 +221,12 @@ public:
         ) const;
 
 
+        // I-O
+
+            //- Write to Ostream
+            //void write(Ostream& os) const;
+
+
     // Member operators
 
 

src/thermophysicalModels/specie/equationOfState/cubicEquationOfState/soaveRedlichKwong/soaveRedlichKwongI.H

@@ -62,7 +62,10 @@ inline soaveRedlichKwong::soaveRedlichKwong(const word& name,const  soaveRedlich
     b_(srk.b_),
     n_(srk.n_),
     rhostd_(srk.rhostd_),
-    TSave(0) 
+    TSave(0),
+    b2_(srk.b2_),
+    b3_(srk.b3_),
+    b5_(srk.b5_) 
 {}
 
 
@@ -167,6 +170,22 @@ inline scalar soaveRedlichKwong::n()const
     return n_;
 }
 
+//CL: pow of constant b() used in the code e.g. b2_=b*b;
+inline scalar soaveRedlichKwong::b2()const
+{
+    return b2_;
+}
+
+inline scalar soaveRedlichKwong::b3()const
+{
+    return b3_;
+}
+
+inline scalar soaveRedlichKwong::b5()const
+{
+    return b5_;
+}
+
 
 //returns the pressure for a given density and temperature    
 inline scalar soaveRedlichKwong::p(const scalar rho,const scalar T) const
@@ -186,13 +205,14 @@ inline scalar soaveRedlichKwong::p(const scalar rho,const scalar T) const
 inline scalar soaveRedlichKwong::dpdv(const scalar rho, const scalar T) const
 {
     scalar Vm = this->W()/rho;
+    scalar Vm2 = Vm*Vm;
 
     return 
     (
-        a(T)*(pow(b_,3)-3*b_*pow(Vm,2)+2*pow(Vm,3))
+        a(T)*(b3()-3*b_*Vm2+2*pow(Vm,3))
-        -this->RR*T*pow(Vm,2)*(pow(b_,2)+2*b_*Vm+pow(Vm,2))
+        -this->RR*T*Vm2*(b2()+2*b_*Vm+Vm2)
     )
-    /(pow(Vm,2)*pow(b_+Vm,2)*pow(b_-Vm,2));   
+    /(Vm2*pow(b_+Vm,2)*pow(b_-Vm,2));   
 }
 
 
@@ -260,20 +280,22 @@ inline scalar soaveRedlichKwong::d2pdT2(const scalar rho,const scalar T) const
 inline scalar soaveRedlichKwong::d2pdv2(const scalar rho,const scalar T) const
 {
     scalar Vm = this->W()/rho;
+    scalar Vm2 = Vm*Vm;
+    scalar Vm3 = Vm*Vm*Vm;
 
     return 
     2*
     (
         a(T)*
         (
-            pow(b_,5)-3*pow(b_,3)*pow(Vm,2)-pow(b_,2)*pow(Vm,3)+6*b_*pow(Vm,4)-3*pow(Vm,5)
+            b5()-3*b3()*Vm2-b2()*Vm3+6*b_*pow(Vm,4)-3*pow(Vm,5)
         )
-        +this->RR*T*pow(Vm,3)*
+        +this->RR*T*Vm3*
         (
-            pow(b_,3)+3*pow(b_,2)*Vm+3*b_*pow(Vm,2)+pow(Vm,3)
+            b3()+3*b2()*Vm+3*b_*Vm2+Vm3
         )
     )
-    /(pow(Vm,3)*pow(b_+Vm,3)*pow(Vm-b_,3));   
+    /(Vm3*pow(b_+Vm,3)*pow(Vm-b_,3));   
 }
 
 
@@ -295,13 +317,14 @@ inline scalar soaveRedlichKwong::d2vdT2(const scalar rho, const scalar T) const
 inline scalar soaveRedlichKwong::d2pdvdT(const scalar rho, const scalar T) const
 {
     scalar Vm = this->W()/rho;
+    scalar Vm2 = Vm*Vm;
 
     return 
     (
-        dadT(T)*(pow(b_,3)-3*b_*pow(Vm,2)+2*pow(Vm,3))
+        dadT(T)*(b3()-3*b_*Vm2+2*pow(Vm,3))
-        -this->RR*pow(Vm,2)*(pow(b_,2)+2*b_*Vm+pow(Vm,2))
+        -this->RR*Vm2*(b2()+2*b_*Vm+Vm2)
     )
-    /(pow(Vm,2)*pow(b_+Vm,2)*pow(b_-Vm,2)); 
+    /(Vm2*pow(b_+Vm,2)*pow(b_-Vm,2)); 
 }
 
 
@@ -411,7 +434,7 @@ inline scalar soaveRedlichKwong::rho(
        {
            FatalErrorIn
            (
-               "inline scalar redlichKwong::rho(const scalar p, const scalar T, const scalar rho0) const "
+               "inline scalar soaveRedlichKwong::rho(const scalar p, const scalar T, const scalar rho0) const "
            )   << "Maximum number of iterations exceeded"
                << abort(FatalError);
        }

src/thermophysicalModels/specie/equationOfState/mixturePengRobinson/mixturePengRobinson.C

@@ -1,83 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     |
-    \\  /    A nd           | Copyright held by original author
-     \\/     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 2 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, write to the Free Software Foundation,
-    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
-
-Description
-    Peng Robionson equation of state for mixtures.
-
-Author
-Christian Lucas
-Institut für Thermodynamik
-Technische Universität Braunschweig 
-Germany
-
-\*---------------------------------------------------------------------------*/
-
-#include "mixturePengRobinson.H"
-#include "IOstreams.H"
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-namespace Foam
-{
-
-// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
-
-mixturePengRobinson::mixturePengRobinson(Istream& is)
-:
-    pengRobinson(is),
-    numOfComp(1),
-    singleComponent(1)
-{ 
-    //CL: set size of weigths, mixtureComponents ... to 10,
-    //CL: when more mixture componentents are used
-    //CL: size of the DynamicLis increases automatically
-    weigths.setSize(10);
-    mixtureComponents.setSize(10);
-    aComponents.setSize(10);
-    daComponents.setSize(10);
-    d2aComponents.setSize(10);
-
-    //Save a pointer of this object in the mixtureComponents array
-    mixtureComponents[0]=this;
-
-    is.check("mixturePengRobinson::mixturePengRobinson(Istream& is)");
-}
-
-
-// * * * * * * * * * * * * * * * Ostream Operator  * * * * * * * * * * * * * //
-
-Ostream& operator<<(Ostream& os, const mixturePengRobinson& pr)
-{
-    os  << static_cast<const pengRobinson&>(pr)<< tab;
-
-    os.check("Ostream& operator<<(Ostream& os, const mixturePengRobinson& st)");
-    return os;
-}
-
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-} // End namespace Foam
-
-// ************************************************************************* //

src/thermophysicalModels/specie/equationOfState/mixtureRedlichKwong/mixtureRedlichKwong.C

@@ -1,78 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     |
-    \\  /    A nd           | Copyright held by original author
-     \\/     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 2 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, write to the Free Software Foundation,
-    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
-
-Description
-    Redlich Kwong equation of state for mixtures.
-
-Author
-Christian Lucas
-Institut für Thermodynamik
-Technische Universität Braunschweig 
-Germany
-
-\*---------------------------------------------------------------------------*/
-
-#include "mixtureRedlichKwong.H"
-#include "IOstreams.H"
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-namespace Foam
-{
-
-// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
-
-mixtureRedlichKwong::mixtureRedlichKwong(Istream& is)
-:
-    redlichKwong(is),
-    numOfComp(1)
-{ 
-    //CL: set size of weigths and mixtureComponents to 10
-    //CL: when more mixture componentents are used
-    //CL: size of the DynamicLis increases automatically
-    weigths.setSize(10);
-    mixtureComponents.setSize(10);
-
-    //Save a pointer of this object in the mixtureComponents array
-    mixtureComponents[0]=this;
-    is.check("mixtureRedlichKwong::mixtureRedlichKwong(Istream& is)");
-}
-
-
-// * * * * * * * * * * * * * * * Ostream Operator  * * * * * * * * * * * * * //
-
-Ostream& operator<<(Ostream& os, const mixtureRedlichKwong& rk)
-{
-    os  << static_cast<const redlichKwong&>(rk)<< tab;
-
-    os.check("Ostream& operator<<(Ostream& os, const mixtureRedlichKwong& st)");
-    return os;
-}
-
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-} // End namespace Foam
-
-// ************************************************************************* //

src/thermophysicalModels/specie/equationOfState/mixtureSoaveRedlichKwong/mixtureSoaveRedlichKwong.C

@@ -1,82 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     |
-    \\  /    A nd           | Copyright held by original author
-     \\/     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 2 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, write to the Free Software Foundation,
-    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
-
-Description
-    Soave Redlich Kwong equation of state for mixtures.
-
-Author
-Christian Lucas
-Institut für Thermodynamik
-Technische Universität Braunschweig 
-Germany
-
-\*---------------------------------------------------------------------------*/
-
-#include "mixtureSoaveRedlichKwong.H"
-#include "IOstreams.H"
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-namespace Foam
-{
-
-// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
-
-mixtureSoaveRedlichKwong::mixtureSoaveRedlichKwong(Istream& is)
-:
-    soaveRedlichKwong(is),
-    numOfComp(1),
-    singleComponent(1)
-{ 
-    //CL: set size of weigths, mixtureComponents ... to 10,
-    //CL: when more mixture componentents are used
-    //CL: size of the DynamicLis increases automatically
-    weigths.setSize(10);
-    mixtureComponents.setSize(10);
-    aComponents.setSize(10);
-    daComponents.setSize(10);
-    d2aComponents.setSize(10);
-
-    //Save a pointer of this object in the mixtureComponents array
-    mixtureComponents[0]=this;
-    is.check("mixtureSoaveRedlichKwong::mixtureSoaveRedlichKwong(Istream& is)");
-}
-
-
-// * * * * * * * * * * * * * * * Ostream Operator  * * * * * * * * * * * * * //
-
-Ostream& operator<<(Ostream& os, const mixtureSoaveRedlichKwong& srk)
-{
-    os  << static_cast<const soaveRedlichKwong&>(srk)<< tab;
-
-    os.check("Ostream& operator<<(Ostream& os, const mixtureSoaveRedlichKwong& st)");
-    return os;
-}
-
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-} // End namespace Foam
-
-// ************************************************************************* //

tutorials/compressible/realFluidPisoFoam/ras/backStep/0.org/T

@@ -21,25 +21,25 @@ internalField uniform 333.15;
 
 boundaryField
 {
-     INLET
+    inlet
     {
         type            fixedValue;
         value           uniform 333.15;
     }
-    OUTLET
+    outlet
     {
         type            zeroGradient;
     }
 
-    WALL_1
+    upperWall
     {
         type            zeroGradient;
     }
-    WALL_2
+    lowerWall
     {
         type            zeroGradient;
     }
-    frontAndBackPlanes
+    frontAndBack
     {
         type            empty;
     }

tutorials/compressible/realFluidPisoFoam/ras/backStep/0.org/U

@@ -17,32 +17,31 @@ FoamFile
 
 dimensions      [0 1 -1 0 0 0 0];
 
-internalField   uniform (5 0 0);
+internalField   uniform (0 0 0);
 
 boundaryField
 {
-     INLET
+    inlet
     {
-        type            inletOutlet;
+        type            pressureInletVelocity;
-        inletValue      uniform (5 0 0);
+        value           uniform (3.5 0 0);
-        value           uniform (5 0 0);
     }
-    OUTLET
+    outlet
     {
         type            zeroGradient;
     }
 
-    WALL_1
+    upperWall
     {
         type            fixedValue;
         value           uniform (0 0 0);
     }
-    WALL_2
+    lowerWall
     {
         type            fixedValue;
         value           uniform (0 0 0);
     }
-    frontAndBackPlanes
+    frontAndBack
     {
         type            empty;
     }

tutorials/compressible/realFluidPisoFoam/ras/backStep/0.org/epsilon

@@ -1,8 +1,8 @@
 /*--------------------------------*- C++ -*----------------------------------*\
 | =========                 |                                                 |
-| \\      /  F ield         | OpenFOAM Extend Project: Open source CFD        |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
-|  \\    /   O peration     | Version:  1.6-ext                               |
+|  \\    /   O peration     | Version:  2.1.1                                 |
-|   \\  /    A nd           | Web:      www.extend-project.de                 |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
 |    \\/     M anipulation  |                                                 |
 \*---------------------------------------------------------------------------*/
 FoamFile
@@ -21,7 +21,7 @@ internalField   uniform 1120;
 
 boundaryField
 {
-    WALL_1
+    upperWall
     {
         type            compressible::epsilonWallFunction;
         Cmu             0.09;
@@ -29,7 +29,7 @@ boundaryField
         E               9.8;
         value           uniform 1120;
     }
-    WALL_2
+    lowerWall
     {
         type            compressible::epsilonWallFunction;
         Cmu             0.09;
@@ -37,18 +37,18 @@ boundaryField
         E               9.8;
         value           uniform 1120;
     }
-    INLET
+    inlet
     {
         type            fixedValue;
         value           uniform 1120;
     }
-    OUTLET
+    outlet
     {
         type            inletOutlet;
         inletValue      uniform 1120;
         value           uniform 1120;
     }
-    frontAndBackPlanes
+    frontAndBack
     {
         type            empty;
     }

tutorials/compressible/realFluidPisoFoam/ras/backStep/0.org/k

@@ -1,8 +1,8 @@
 /*--------------------------------*- C++ -*----------------------------------*\
 | =========                 |                                                 |
-| \\      /  F ield         | OpenFOAM Extend Project: Open source CFD        |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
-|  \\    /   O peration     | Version:  1.6-ext                               |
+|  \\    /   O peration     | Version:  2.1.1                                 |
-|   \\  /    A nd           | Web:      www.extend-project.de                 |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
 |    \\/     M anipulation  |                                                 |
 \*---------------------------------------------------------------------------*/
 FoamFile
@@ -21,27 +21,29 @@ internalField   uniform 5;
 
 boundaryField
 {
-    WALL_1
+    upperWall
     {
         type            compressible::kqRWallFunction;
         value           uniform 5;
     }
-    WALL_2
+    lowerWall
     {
         type            compressible::kqRWallFunction;
         value           uniform 5;
     }
-    INLET
+    inlet
     {
         type            turbulentIntensityKineticEnergyInlet;
         intensity       0.01;
+        U               U;
+        phi             phi;
         value           uniform 5;
     }
-    OUTLET
+    outlet
     {
         type            zeroGradient;
     }
-    frontAndBackPlanes
+    frontAndBack
     {
         type            empty;
     }

tutorials/compressible/realFluidPisoFoam/ras/backStep/0.org/p

@@ -21,30 +21,30 @@ internalField   uniform 80e5;
 
 boundaryField
 {
-     INLET
+    inlet
     {
         type            totalPressure;
         rho             rho;
         psi             none;
         gamma           1.4;
-        p0              uniform 8.01e+06;
+        p0              uniform 80.01e5;
-        value           uniform 8.01e+06;
+        value           uniform 80.01e5;
     }
-    OUTLET
+    outlet
     {
         type            fixedValue;
         value           uniform 8e+06;
     }
 
-    WALL_1
+    upperWall
     {
         type            zeroGradient;
     }
-    WALL_2
+    lowerWall
     {
         type            zeroGradient;
     }
-    frontAndBackPlanes
+    frontAndBack
     {
         type            empty;
     }

tutorials/compressible/realFluidPisoFoam/ras/backStep/0/T

@@ -11,35 +11,35 @@ FoamFile
     format      ascii;
     class       volScalarField;
     location    "0";
-    object      p;
+    object      T;
 }
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-dimensions      [1 -1 -2 0 0 0 0];
+dimensions      [0 0 0 1 0 0 0];
 
-internalField   uniform 80e5;
+internalField uniform 333.15;
 
 boundaryField
 {
-     INLET
+    inlet
-    {
-   type            zeroGradient;
-    }
-    OUTLET
     {
         type            fixedValue;
-        value           uniform 8e+06;
+        value           uniform 333.15;
+    }
+    outlet
+    {
+        type            zeroGradient;
     }
 
-    WALL_1
+    upperWall
     {
         type            zeroGradient;
     }
-    WALL_2
+    lowerWall
     {
         type            zeroGradient;
     }
-    frontAndBackPlanes
+    frontAndBack
     {
         type            empty;
     }

tutorials/compressible/realFluidPisoFoam/ras/backStep/0/U

@@ -17,33 +17,31 @@ FoamFile
 
 dimensions      [0 1 -1 0 0 0 0];
 
-internalField   uniform (4.8 0 0);
+internalField   uniform (0 0 0);
 
 boundaryField
 {
-     INLET
+    inlet
     {
         type            pressureInletVelocity;
-        value           uniform (4.8 0 0);
+        value           uniform (3.5 0 0);
     }
-    OUTLET
+    outlet
     {
-        type            inletOutlet;
+        type            zeroGradient;
-        inletValue      uniform (0 0 0);
-        value           uniform (4.8 0 0);
     }
 
-    WALL_1
+    upperWall
     {
         type            fixedValue;
         value           uniform (0 0 0);
     }
-    WALL_2
+    lowerWall
     {
         type            fixedValue;
         value           uniform (0 0 0);
     }
-    frontAndBackPlanes
+    frontAndBack
     {
         type            empty;
     }

tutorials/compressible/realFluidPisoFoam/ras/backStep/0/epsilon

@@ -1,8 +1,8 @@
 /*--------------------------------*- C++ -*----------------------------------*\
 | =========                 |                                                 |
-| \\      /  F ield         | OpenFOAM Extend Project: Open source CFD        |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
-|  \\    /   O peration     | Version:  1.6-ext                               |
+|  \\    /   O peration     | Version:  2.1.1                                 |
-|   \\  /    A nd           | Web:      www.extend-project.de                 |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
 |    \\/     M anipulation  |                                                 |
 \*---------------------------------------------------------------------------*/
 FoamFile
@@ -21,7 +21,7 @@ internalField   uniform 1120;
 
 boundaryField
 {
-    WALL_1
+    upperWall
     {
         type            compressible::epsilonWallFunction;
         Cmu             0.09;
@@ -29,7 +29,7 @@ boundaryField
         E               9.8;
         value           uniform 1120;
     }
-    WALL_2
+    lowerWall
     {
         type            compressible::epsilonWallFunction;
         Cmu             0.09;
@@ -37,18 +37,18 @@ boundaryField
         E               9.8;
         value           uniform 1120;
     }
-    INLET
+    inlet
     {
         type            fixedValue;
         value           uniform 1120;
     }
-    OUTLET
+    outlet
     {
         type            inletOutlet;
         inletValue      uniform 1120;
         value           uniform 1120;
     }
-    frontAndBackPlanes
+    frontAndBack
     {
         type            empty;
     }

tutorials/compressible/realFluidPisoFoam/ras/backStep/0/k

@@ -1,8 +1,8 @@
 /*--------------------------------*- C++ -*----------------------------------*\
 | =========                 |                                                 |
-| \\      /  F ield         | OpenFOAM Extend Project: Open source CFD        |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
-|  \\    /   O peration     | Version:  1.6-ext                               |
+|  \\    /   O peration     | Version:  2.1.1                                 |
-|   \\  /    A nd           | Web:      www.extend-project.de                 |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
 |    \\/     M anipulation  |                                                 |
 \*---------------------------------------------------------------------------*/
 FoamFile
@@ -21,27 +21,29 @@ internalField   uniform 5;
 
 boundaryField
 {
-    WALL_1
+    upperWall
     {
         type            compressible::kqRWallFunction;
         value           uniform 5;
     }
-    WALL_2
+    lowerWall
     {
         type            compressible::kqRWallFunction;
         value           uniform 5;
     }
-    INLET
+    inlet
     {
         type            turbulentIntensityKineticEnergyInlet;
         intensity       0.01;
+        U               U;
+        phi             phi;
         value           uniform 5;
     }
-    OUTLET
+    outlet
     {
         type            zeroGradient;
     }
-    frontAndBackPlanes
+    frontAndBack
     {
         type            empty;
     }

tutorials/compressible/realFluidPisoFoam/ras/backStep/0/p

@@ -11,37 +11,40 @@ FoamFile
     format      ascii;
     class       volScalarField;
     location    "0";
-    object      T;
+    object      p;
 }
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-dimensions      [0 0 0 1 0 0 0];
+dimensions      [1 -1 -2 0 0 0 0];
 
-internalField uniform 533.15;
+internalField   uniform 80e5;
 
 boundaryField
 {
-     INLET
+    inlet
+    {
+        type            totalPressure;
+        rho             rho;
+        psi             none;
+        gamma           1.4;
+        p0              uniform 80.01e5;
+        value           uniform 80.01e5;
+    }
+    outlet
     {
         type            fixedValue;
-        value           uniform 533.15;
+        value           uniform 8e+06;
-    }
-    OUTLET
-    {
-        type            inletOutlet;
-        inletValue      uniform 533.15;
-        value           uniform 533.15;
     }
 
-    WALL_1
+    upperWall
     {
         type            zeroGradient;
     }
-    WALL_2
+    lowerWall
     {
         type            zeroGradient;
     }
-    frontAndBackPlanes
+    frontAndBack
     {
         type            empty;
     }

tutorials/compressible/realFluidPisoFoam/ras/backStep/Allrun

@@ -0,0 +1,9 @@
+#!/bin/sh
+cd ${0%/*} || exit 1    # run from this directory
+
+
+# Source tutorial run functions
+. $WM_PROJECT_DIR/bin/tools/RunFunctions
+
+runApplication blockMesh
+runApplication `getApplication`

tutorials/compressible/realFluidPisoFoam/ras/backStep/constant/polyMesh/blockMeshDict

@@ -0,0 +1,173 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  2.1.1                                 |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    object      blockMeshDict;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+convertToMeters 0.001;
+
+vertices
+(
+    (-20.6 0 -0.5)
+    (-20.6 3 -0.5)
+    (-20.6 12.7 -0.5)
+    (-20.6 25.4 -0.5)
+    (0 -25.4 -0.5)
+    (0 -5 -0.5)
+    (0 0 -0.5)
+    (0 3 -0.5)
+    (0 12.7 -0.5)
+    (0 25.4 -0.5)
+    (206 -25.4 -0.5)
+    (206 -8.5 -0.5)
+    (206 0 -0.5)
+    (206 6.5 -0.5)
+    (206 17 -0.5)
+    (206 25.4 -0.5)
+    (290 -16.6 -0.5)
+    (290 -6.3 -0.5)
+    (290 0 -0.5)
+    (290 4.5 -0.5)
+    (290 11 -0.5)
+    (290 16.6 -0.5)
+    (-20.6 0 0.5)
+    (-20.6 3 0.5)
+    (-20.6 12.7 0.5)
+    (-20.6 25.4 0.5)
+    (0 -25.4 0.5)
+    (0 -5 0.5)
+    (0 0 0.5)
+    (0 3 0.5)
+    (0 12.7 0.5)
+    (0 25.4 0.5)
+    (206 -25.4 0.5)
+    (206 -8.5 0.5)
+    (206 0 0.5)
+    (206 6.5 0.5)
+    (206 17 0.5)
+    (206 25.4 0.5)
+    (290 -16.6 0.5)
+    (290 -6.3 0.5)
+    (290 0 0.5)
+    (290 4.5 0.5)
+    (290 11 0.5)
+    (290 16.6 0.5)
+);
+
+blocks
+(
+    hex (0 6 7 1 22 28 29 23) (18 7 1) simpleGrading (0.5 1.8 1)
+    hex (1 7 8 2 23 29 30 24) (18 10 1) simpleGrading (0.5 4 1)
+    hex (2 8 9 3 24 30 31 25) (18 13 1) simpleGrading (0.5 0.25 1)
+    hex (4 10 11 5 26 32 33 27) (180 18 1) simpleGrading (4 1 1)
+    hex (5 11 12 6 27 33 34 28) (180 9 1) edgeGrading (4 4 4 4 0.5 1 1 0.5 1 1 1 1)
+    hex (6 12 13 7 28 34 35 29) (180 7 1) edgeGrading (4 4 4 4 1.8 1 1 1.8 1 1 1 1)
+    hex (7 13 14 8 29 35 36 30) (180 10 1) edgeGrading (4 4 4 4 4 1 1 4 1 1 1 1)
+    hex (8 14 15 9 30 36 37 31) (180 13 1) simpleGrading (4 0.25 1)
+    hex (10 16 17 11 32 38 39 33) (25 18 1) simpleGrading (2.5 1 1)
+    hex (11 17 18 12 33 39 40 34) (25 9 1) simpleGrading (2.5 1 1)
+    hex (12 18 19 13 34 40 41 35) (25 7 1) simpleGrading (2.5 1 1)
+    hex (13 19 20 14 35 41 42 36) (25 10 1) simpleGrading (2.5 1 1)
+    hex (14 20 21 15 36 42 43 37) (25 13 1) simpleGrading (2.5 0.25 1)
+);
+
+edges
+(
+);
+
+boundary
+(
+    inlet
+    {
+        type patch;
+        faces
+        (
+            (0 22 23 1)
+            (1 23 24 2)
+            (2 24 25 3)
+        );
+    }
+    outlet
+    {
+        type patch;
+        faces
+        (
+            (16 17 39 38)
+            (17 18 40 39)
+            (18 19 41 40)
+            (19 20 42 41)
+            (20 21 43 42)
+        );
+    }
+    upperWall
+    {
+        type wall;
+        faces
+        (
+            (3 25 31 9)
+            (9 31 37 15)
+            (15 37 43 21)
+        );
+    }
+    lowerWall
+    {
+        type wall;
+        faces
+        (
+            (0 6 28 22)
+            (6 5 27 28)
+            (5 4 26 27)
+            (4 10 32 26)
+            (10 16 38 32)
+        );
+    }
+    frontAndBack
+    {
+        type empty;
+        faces
+        (
+            (22 28 29 23)
+            (23 29 30 24)
+            (24 30 31 25)
+            (26 32 33 27)
+            (27 33 34 28)
+            (28 34 35 29)
+            (29 35 36 30)
+            (30 36 37 31)
+            (32 38 39 33)
+            (33 39 40 34)
+            (34 40 41 35)
+            (35 41 42 36)
+            (36 42 43 37)
+            (0 1 7 6)
+            (1 2 8 7)
+            (2 3 9 8)
+            (4 5 11 10)
+            (5 6 12 11)
+            (6 7 13 12)
+            (7 8 14 13)
+            (8 9 15 14)
+            (10 11 17 16)
+            (11 12 18 17)
+            (12 13 19 18)
+            (13 14 20 19)
+            (14 15 21 20)
+        );
+    }
+);
+
+mergePatchPairs
+(
+);
+
+// ************************************************************************* //

tutorials/compressible/realFluidPisoFoam/ras/backStep/constant/polyMesh/boundary

@@ -1,8 +1,8 @@
 /*--------------------------------*- C++ -*----------------------------------*\
 | =========                 |                                                 |
-| \\      /  F ield         | OpenFOAM Extend Project: Open source CFD        |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
-|  \\    /   O peration     | Version:  1.6-ext                               |
+|  \\    /   O peration     | Version:  2.1.1                                 |
-|   \\  /    A nd           | Web:      www.extend-project.de                 |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
 |    \\/     M anipulation  |                                                 |
 \*---------------------------------------------------------------------------*/
 FoamFile
@@ -17,35 +17,35 @@ FoamFile
 
 5
 (
-    WALL_1
+    inlet
-    {
-        type            wall;
-        nFaces          99;
-        startFace       3644;
-    }
-    WALL_2
-    {
-        type            wall;
-        nFaces          99;
-        startFace       3743;
-    }
-    INLET
     {
         type            patch;
-        nFaces          19;
+        nFaces          30;
-        startFace       3842;
+        startFace       24170;
     }
-    OUTLET
+    outlet
     {
         type            patch;
-        nFaces          19;
+        nFaces          57;
-        startFace       3861;
+        startFace       24200;
     }
-    frontAndBackPlanes
+    upperWall
+    {
+        type            wall;
+        nFaces          223;
+        startFace       24257;
+    }
+    lowerWall
+    {
+        type            wall;
+        nFaces          250;
+        startFace       24480;
+    }
+    frontAndBack
     {
         type            empty;
-        nFaces          3762;
+        nFaces          24450;
-        startFace       3880;
+        startFace       24730;
     }
 )
 

tutorials/compressible/realFluidPisoFoam/ras/backStep/system/controlDict

@@ -15,12 +15,10 @@ FoamFile
 }
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-application     rhoPisoFoam;
+application     realFluidPisoFoam;
 
 startFrom       latestTime;
 
-//startTime       0;
-
 stopAt          endTime;
 
 endTime         0.5;
@@ -29,7 +27,7 @@ deltaT          1e-5;
 
 writeControl    runTime;
 
-writeInterval  1e-2;
+writeInterval   1e-2;
 
 purgeWrite      0;
 
@@ -45,7 +43,7 @@ timePrecision   10;
 
 adjustTimeStep  yes;
 
-maxCo          0.5;
+maxCo           0.5;
 
 maxDeltaT       1e-2;
 

tutorials/compressible/realFluidPisoFoam/ras/backStep/system/fvSchemes

@@ -28,7 +28,7 @@ gradSchemes
 
 divSchemes
 {
-    default         none;
+    default         Gauss linear;
     div(phi,U)      Gauss limitedLinearV 1;
     div(phid,p)     Gauss limitedLinear 1;
     div(phiU,p)     Gauss linear;

tutorials/compressible/realFluidPisoFoam/ras/cavity_IAPWS97/0.org/T

@@ -0,0 +1,39 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  2.1.x                                 |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       volScalarField;
+    object      T;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [0 0 0 1 0 0 0];
+
+internalField   uniform 300;
+
+boundaryField
+{
+    movingWall
+    {
+        type            zeroGradient;
+    }
+
+    fixedWalls
+    {
+        type            zeroGradient;
+    }
+
+    frontAndBack
+    {
+        type            empty;
+    }
+}
+
+// ************************************************************************* //

tutorials/compressible/realFluidPisoFoam/ras/cavity_IAPWS97/0.org/U

@@ -0,0 +1,43 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  2.1.x                                 |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       volVectorField;
+    object      U;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [0 1 -1 0 0 0 0];
+
+internalField   uniform (0 0 0);
+
+boundaryField
+{
+    movingWall
+    {
+        type            fixedValue;
+        // Field Value
+        value           uniform (1 0 0);
+    }
+
+    fixedWalls
+    {
+        type            fixedValue;
+        // Field Value
+        value           uniform (0 0 0);
+    }
+
+    frontAndBack
+    {
+        type            empty;
+    }
+}
+
+// ************************************************************************* //

tutorials/compressible/realFluidPisoFoam/ras/cavity_IAPWS97/0.org/epsilon

@@ -0,0 +1,43 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  2.1.x                                 |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       volScalarField;
+    location    "0";
+    object      epsilon;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [ 0 2 -3 0 0 0 0 ];
+
+internalField   uniform 0.000765;
+
+boundaryField
+{
+    movingWall
+    {
+        type            compressible::epsilonWallFunction;
+        value           uniform 0;
+    }
+
+    fixedWalls
+    {
+        type            compressible::epsilonWallFunction;
+        value           uniform 0;
+    }
+
+    frontAndBack
+    {
+        type            empty;
+    }
+}
+
+
+// ************************************************************************* //

tutorials/compressible/realFluidPisoFoam/ras/cavity_IAPWS97/0.org/k

@@ -0,0 +1,41 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  2.1.x                                 |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       volScalarField;
+    location    "0";
+    object      k;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [0 2 -2 0 0 0 0];
+
+internalField   uniform 0.00325;
+
+boundaryField
+{
+    movingWall
+    {
+        type            compressible::kqRWallFunction;
+        value           uniform 0.00325;
+    }
+    fixedWalls
+    {
+        type            compressible::kqRWallFunction;
+        value           uniform 0.00325;
+    }
+    frontAndBack
+    {
+        type            empty;
+    }
+}
+
+
+// ************************************************************************* //

tutorials/compressible/realFluidPisoFoam/ras/cavity_IAPWS97/0.org/omega

@@ -0,0 +1,47 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  2.1.x                                 |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       volScalarField;
+    location    "0";
+    object      omega;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [0 0 -1 0 0 0 0];
+
+internalField   uniform 2.6;
+
+boundaryField
+{
+    movingWall
+    {
+        type            compressible::omegaWallFunction;
+        Cmu             0.09;
+        kappa           0.41;
+        E               9.8;
+        value           uniform 2.6;
+    }
+    fixedWalls
+    {
+        type            compressible::omegaWallFunction;
+        Cmu             0.09;
+        kappa           0.41;
+        E               9.8;
+        value           uniform 2.6;
+    }
+    frontAndBack
+    {
+        type            empty;
+    }
+}
+
+
+// ************************************************************************* //

tutorials/compressible/realFluidPisoFoam/ras/cavity_IAPWS97/0.org/p

@@ -0,0 +1,39 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  2.1.x                                 |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       volScalarField;
+    object      p;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [1 -1 -2 0 0 0 0];
+
+internalField   uniform 100000;
+
+boundaryField
+{
+    movingWall
+    {
+        type            zeroGradient;
+    }
+
+    fixedWalls
+    {
+        type            zeroGradient;
+    }
+
+    frontAndBack
+    {
+        type            empty;
+    }
+}
+
+// ************************************************************************* //

tutorials/compressible/realFluidPisoFoam/ras/cavity_IAPWS97/0/T

@@ -0,0 +1,39 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  2.1.x                                 |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       volScalarField;
+    object      T;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [0 0 0 1 0 0 0];
+
+internalField   uniform 300;
+
+boundaryField
+{
+    movingWall
+    {
+        type            zeroGradient;
+    }
+
+    fixedWalls
+    {
+        type            zeroGradient;
+    }
+
+    frontAndBack
+    {
+        type            empty;
+    }
+}
+
+// ************************************************************************* //

tutorials/compressible/realFluidPisoFoam/ras/cavity_IAPWS97/0/U

@@ -0,0 +1,43 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  2.1.x                                 |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       volVectorField;
+    object      U;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [0 1 -1 0 0 0 0];
+
+internalField   uniform (0 0 0);
+
+boundaryField
+{
+    movingWall
+    {
+        type            fixedValue;
+        // Field Value
+        value           uniform (1 0 0);
+    }
+
+    fixedWalls
+    {
+        type            fixedValue;
+        // Field Value
+        value           uniform (0 0 0);
+    }
+
+    frontAndBack
+    {
+        type            empty;
+    }
+}
+
+// ************************************************************************* //

tutorials/compressible/realFluidPisoFoam/ras/cavity_IAPWS97/0/epsilon

@@ -0,0 +1,43 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  2.1.x                                 |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       volScalarField;
+    location    "0";
+    object      epsilon;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [ 0 2 -3 0 0 0 0 ];
+
+internalField   uniform 0.000765;
+
+boundaryField
+{
+    movingWall
+    {
+        type            compressible::epsilonWallFunction;
+        value           uniform 0;
+    }
+
+    fixedWalls
+    {
+        type            compressible::epsilonWallFunction;
+        value           uniform 0;
+    }
+
+    frontAndBack
+    {
+        type            empty;
+    }
+}
+
+
+// ************************************************************************* //

tutorials/compressible/realFluidPisoFoam/ras/cavity_IAPWS97/0/k

@@ -0,0 +1,41 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  2.1.x                                 |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       volScalarField;
+    location    "0";
+    object      k;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [0 2 -2 0 0 0 0];
+
+internalField   uniform 0.00325;
+
+boundaryField
+{
+    movingWall
+    {
+        type            compressible::kqRWallFunction;
+        value           uniform 0.00325;
+    }
+    fixedWalls
+    {
+        type            compressible::kqRWallFunction;
+        value           uniform 0.00325;
+    }
+    frontAndBack
+    {
+        type            empty;
+    }
+}
+
+
+// ************************************************************************* //

tutorials/compressible/realFluidPisoFoam/ras/cavity_IAPWS97/0/omega

@@ -0,0 +1,47 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  2.1.x                                 |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       volScalarField;
+    location    "0";
+    object      omega;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [0 0 -1 0 0 0 0];
+
+internalField   uniform 2.6;
+
+boundaryField
+{
+    movingWall
+    {
+        type            compressible::omegaWallFunction;
+        Cmu             0.09;
+        kappa           0.41;
+        E               9.8;
+        value           uniform 2.6;
+    }
+    fixedWalls
+    {
+        type            compressible::omegaWallFunction;
+        Cmu             0.09;
+        kappa           0.41;
+        E               9.8;
+        value           uniform 2.6;
+    }
+    frontAndBack
+    {
+        type            empty;
+    }
+}
+
+
+// ************************************************************************* //

tutorials/compressible/realFluidPisoFoam/ras/cavity_IAPWS97/0/p

@@ -0,0 +1,39 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  2.1.x                                 |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       volScalarField;
+    object      p;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [1 -1 -2 0 0 0 0];
+
+internalField   uniform 100000;
+
+boundaryField
+{
+    movingWall
+    {
+        type            zeroGradient;
+    }
+
+    fixedWalls
+    {
+        type            zeroGradient;
+    }
+
+    frontAndBack
+    {
+        type            empty;
+    }
+}
+
+// ************************************************************************* //

tutorials/compressible/realFluidPisoFoam/ras/cavity_IAPWS97/Allrun

@@ -0,0 +1,9 @@
+#!/bin/sh
+cd ${0%/*} || exit 1    # run from this directory
+
+
+# Source tutorial run functions
+. $WM_PROJECT_DIR/bin/tools/RunFunctions
+
+runApplication blockMesh
+runApplication `getApplication`

tutorials/compressible/realFluidPisoFoam/ras/cavity_IAPWS97/README

@@ -14,7 +14,7 @@ Important files in this case:
 
 1. fvSolution 
 
-PISO subdict --> see the realFluid --> this must be true (see modification in the pEqn of rhoPisoFoam)
+PISO subdict --> see the realFluid --> this must be true (see modification in the pEqn of realFluidPisoFoam)
 
 2. controlDict: 
 

tutorials/compressible/realFluidPisoFoam/ras/cavity_IAPWS97/constant/polyMesh/blockMeshDict

@@ -0,0 +1,75 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  2.1.x                                 |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    object      blockMeshDict;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+convertToMeters 0.1;
+
+vertices
+(
+    (0 0 0)
+    (1 0 0)
+    (1 1 0)
+    (0 1 0)
+    (0 0 0.1)
+    (1 0 0.1)
+    (1 1 0.1)
+    (0 1 0.1)
+);
+
+blocks
+(
+    hex (0 1 2 3 4 5 6 7) (20 20 1) simpleGrading (1 1 1)
+);
+
+edges
+(
+);
+
+boundary
+(
+    movingWall
+    {
+        type wall;
+        faces
+        (
+            (3 7 6 2)
+        );
+    }
+    fixedWalls
+    {
+        type wall;
+        faces
+        (
+            (0 4 7 3)
+            (2 6 5 1)
+            (1 5 4 0)
+        );
+    }
+    frontAndBack
+    {
+        type empty;
+        faces
+        (
+            (0 3 2 1)
+            (4 5 6 7)
+        );
+    }
+);
+
+mergePatchPairs
+(
+);
+
+// ************************************************************************* //

tutorials/compressible/realFluidPisoFoam/ras/cavity_IAPWS97/constant/polyMesh/boundary

@@ -0,0 +1,40 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  2.1.1                                 |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       polyBoundaryMesh;
+    location    "constant/polyMesh";
+    object      boundary;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+3
+(
+    movingWall
+    {
+        type            wall;
+        nFaces          20;
+        startFace       760;
+    }
+    fixedWalls
+    {
+        type            wall;
+        nFaces          60;
+        startFace       780;
+    }
+    frontAndBack
+    {
+        type            empty;
+        nFaces          800;
+        startFace       840;
+    }
+)
+
+// ************************************************************************* //

tutorials/compressible/realFluidPisoFoam/ras/cavity_IAPWS97/system/controlDict

@@ -15,21 +15,20 @@ FoamFile
 }
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-application     rhoPisoFoam;
+application     realFluidPisoFoam;
 
 startFrom       latestTime;
 
-//startTime       0;
-
 stopAt          endTime;
 
-endTime         2.5;
+endTime         0.5;
 
 deltaT          1e-5;
 
 writeControl    runTime;
 
-writeInterval  1e-2;
+writeInterval   0.1ls
+;
 
 purgeWrite      0;
 
@@ -45,12 +44,13 @@ timePrecision   10;
 
 adjustTimeStep  yes;
 
-maxCo          0.5;
+maxCo           0.5;
 
 maxDeltaT       1e-3;
 
 runTimeModifiable yes;
 
+//CL: load libraries need for the water properties
 libs
 (
     "libIAPWSThermo.so"

tutorials/compressible/realFluidPisoFoam/ras/cavity_IAPWS97/system/fvSchemes

@@ -28,7 +28,7 @@ gradSchemes
 
 divSchemes
 {
-    default         none;
+    default         Gauss linear;
     div(phi,U)      Gauss limitedLinearV 0.5;
     div(phid,p)     Gauss limitedLinear  0.5;
     div(phiU,p)     Gauss limitedLinear  0.5;

tutorials/compressible/rhoPisoFoam/ras/IAPWS97_Pipe/constant/polyMesh/cellZones

@@ -1,21 +0,0 @@
-/*--------------------------------*- C++ -*----------------------------------*\
-| =========                 |                                                 |
-| \\      /  F ield         | OpenFOAM Extend Project: Open source CFD        |
-|  \\    /   O peration     | Version:  1.6-ext                               |
-|   \\  /    A nd           | Web:      www.extend-project.de                 |
-|    \\/     M anipulation  |                                                 |
-\*---------------------------------------------------------------------------*/
-FoamFile
-{
-    version     2.0;
-    format      ascii;
-    class       regIOobject;
-    location    "constant/polyMesh";
-    object      cellZones;
-}
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-0
-()
-
-// ************************************************************************* //

tutorials/compressible/rhoPisoFoam/ras/IAPWS97_Pipe/constant/polyMesh/faceZones

@@ -1,21 +0,0 @@
-/*--------------------------------*- C++ -*----------------------------------*\
-| =========                 |                                                 |
-| \\      /  F ield         | OpenFOAM Extend Project: Open source CFD        |
-|  \\    /   O peration     | Version:  1.6-ext                               |
-|   \\  /    A nd           | Web:      www.extend-project.de                 |
-|    \\/     M anipulation  |                                                 |
-\*---------------------------------------------------------------------------*/
-FoamFile
-{
-    version     2.0;
-    format      ascii;
-    class       regIOobject;
-    location    "constant/polyMesh";
-    object      faceZones;
-}
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-0
-()
-
-// ************************************************************************* //

tutorials/compressible/rhoPisoFoam/ras/IAPWS97_Pipe/constant/polyMesh/pointZones

@@ -1,21 +0,0 @@
-/*--------------------------------*- C++ -*----------------------------------*\
-| =========                 |                                                 |
-| \\      /  F ield         | OpenFOAM Extend Project: Open source CFD        |
-|  \\    /   O peration     | Version:  1.6-ext                               |
-|   \\  /    A nd           | Web:      www.extend-project.de                 |
-|    \\/     M anipulation  |                                                 |
-\*---------------------------------------------------------------------------*/
-FoamFile
-{
-    version     2.0;
-    format      ascii;
-    class       regIOobject;
-    location    "constant/polyMesh";
-    object      pointZones;
-}
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-0
-()
-
-// ************************************************************************* //

tutorials/compressible/rhoPisoFoam/ras/IAPWS97_Pipe/constant/polyMesh/zoneToPatchName

@@ -1,40 +0,0 @@
-/*--------------------------------*- C++ -*----------------------------------*\
-| =========                 |                                                 |
-| \\      /  F ield         | OpenFOAM Extend Project: Open source CFD        |
-|  \\    /   O peration     | Version:  1.6-ext                               |
-|   \\  /    A nd           | Web:      www.extend-project.de                 |
-|    \\/     M anipulation  |                                                 |
-\*---------------------------------------------------------------------------*/
-FoamFile
-{
-    version     2.0;
-    format      ascii;
-    class       wordList;
-    location    "constant/polyMesh";
-    object      zoneToPatchName;
-}
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-
-16
-(
-unknown
-unknown
-unknown
-unknown
-unknown
-unknown
-unknown
-unknown
-unknown
-unknown
-unknown
-int_SOLID
-WALL_1
-WALL_2
-INLET
-OUTLET
-)
-
-
-// ************************************************************************* //

tutorials/compressible/rhoPisoFoam/ras/pipe/constant/polyMesh/boundary

@@ -1,52 +0,0 @@
-/*--------------------------------*- C++ -*----------------------------------*\
-| =========                 |                                                 |
-| \\      /  F ield         | OpenFOAM Extend Project: Open source CFD        |
-|  \\    /   O peration     | Version:  1.6-ext                               |
-|   \\  /    A nd           | Web:      www.extend-project.de                 |
-|    \\/     M anipulation  |                                                 |
-\*---------------------------------------------------------------------------*/
-FoamFile
-{
-    version     2.0;
-    format      ascii;
-    class       polyBoundaryMesh;
-    location    "constant/polyMesh";
-    object      boundary;
-}
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-5
-(
-    WALL_1
-    {
-        type            wall;
-        nFaces          99;
-        startFace       3644;
-    }
-    WALL_2
-    {
-        type            wall;
-        nFaces          99;
-        startFace       3743;
-    }
-    INLET
-    {
-        type            patch;
-        nFaces          19;
-        startFace       3842;
-    }
-    OUTLET
-    {
-        type            patch;
-        nFaces          19;
-        startFace       3861;
-    }
-    frontAndBackPlanes
-    {
-        type            empty;
-        nFaces          3762;
-        startFace       3880;
-    }
-)
-
-// ************************************************************************* //

tutorials/compressible/rhoPisoFoam/ras/pipe/constant/polyMesh/cellZones

@@ -1,21 +0,0 @@
-/*--------------------------------*- C++ -*----------------------------------*\
-| =========                 |                                                 |
-| \\      /  F ield         | OpenFOAM Extend Project: Open source CFD        |
-|  \\    /   O peration     | Version:  1.6-ext                               |
-|   \\  /    A nd           | Web:      www.extend-project.de                 |
-|    \\/     M anipulation  |                                                 |
-\*---------------------------------------------------------------------------*/
-FoamFile
-{
-    version     2.0;
-    format      ascii;
-    class       regIOobject;
-    location    "constant/polyMesh";
-    object      cellZones;
-}
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-0
-()
-
-// ************************************************************************* //

tutorials/compressible/rhoPisoFoam/ras/pipe/constant/polyMesh/faceZones

@@ -1,21 +0,0 @@
-/*--------------------------------*- C++ -*----------------------------------*\
-| =========                 |                                                 |
-| \\      /  F ield         | OpenFOAM Extend Project: Open source CFD        |
-|  \\    /   O peration     | Version:  1.6-ext                               |
-|   \\  /    A nd           | Web:      www.extend-project.de                 |
-|    \\/     M anipulation  |                                                 |
-\*---------------------------------------------------------------------------*/
-FoamFile
-{
-    version     2.0;
-    format      ascii;
-    class       regIOobject;
-    location    "constant/polyMesh";
-    object      faceZones;
-}
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-0
-()
-
-// ************************************************************************* //