changed realGas EqnOfState + nasaPoly + realGasSpecieThermo

src/thermophysicalModels/basic/Make/files

@@ -8,6 +8,7 @@ psiThermo/basicPsiThermo/newBasicPsiThermo.C
 psiThermo/hPsiThermo/hPsiThermos.C
 psiThermo/hsPsiThermo/hsPsiThermos.C
 psiThermo/ePsiThermo/ePsiThermos.C
+psiThermo/realGasHThermo/realGasHThermos.C
 
 rhoThermo/basicRhoThermo/basicRhoThermo.C
 rhoThermo/basicRhoThermo/newBasicRhoThermo.C
@@ -28,15 +29,4 @@ derivedFvPatchFields/temperatureDirectedInletOutletVelocity/temperatureDirectedI
 derivedFvPatchFields/isentropicTotalPressure/isentropicTotalPressureFvPatchScalarField.C
 derivedFvPatchFields/isentropicTotalTemperature/isentropicTotalTemperatureFvPatchScalarField.C
 
-
-psiThermo/realGasHThermo/realGasHThermos.C
-
-derivedFvPatchFieldsRealGas/fixedEnthalpyRealFluids/fixedEnthalpyRealFluids.C
-derivedFvPatchFieldsRealGas/fixedInternalEnergyRealFluids/fixedInternalEnergyRealFluids.C
-derivedFvPatchFieldsRealGas/gradientEnthalpyRealFluids/gradientEnthalpyRealFluids.C
-derivedFvPatchFieldsRealGas/gradientInternalEnergyRealFluids/gradientInternalEnergyRealFluids.C
-derivedFvPatchFieldsRealGas/mixedEnthalpyRealFluids/mixedEnthalpyRealFluids.C
-derivedFvPatchFieldsRealGas/mixedInternalEnergyRealFluids/mixedInternalEnergyRealFluids.C
-derivedFvPatchFieldsRealGas/realFluidWallHeatTransfer/realFluidWallHeatTransferFvPatchScalarField.C
-
 LIB = $(FOAM_LIBBIN)/libbasicThermophysicalModels

src/thermophysicalModels/basic/basicThermo/basicThermo.C

@@ -33,13 +33,6 @@ License
 #include "fixedInternalEnergyFvPatchScalarField.H"
 #include "gradientInternalEnergyFvPatchScalarField.H"
 #include "mixedInternalEnergyFvPatchScalarField.H"
-#include "zeroGradientFvPatchFields.H"
-#include "fixedEnthalpyRealFluids.H"
-#include "gradientEnthalpyRealFluids.H"
-#include "mixedEnthalpyRealFluids.H"
-#include "fixedInternalEnergyRealFluids.H"
-#include "gradientInternalEnergyRealFluids.H"
-#include "mixedInternalEnergyRealFluids.H"
 
 /* * * * * * * * * * * * * * Private Static Data * * * * * * * * * * * * * * */
 
@@ -149,106 +142,6 @@ void Foam::basicThermo::eBoundaryCorrection(volScalarField& e)
     }
 }
 
-//real Gas Boundary Conditions
-
-Foam::wordList Foam::basicThermo::hRealBoundaryTypes()
-{
-    const volScalarField::GeometricBoundaryField& tbf = T_.boundaryField();
-
-    wordList hbt = tbf.types();
-
-    forAll(tbf, patchi)
-    {
-        if (isA<fixedValueFvPatchScalarField>(tbf[patchi]))
-        {
-            hbt[patchi] = fixedEnthalpyRealFluids::typeName;
-        }
-        else if
-        (
-            isA<zeroGradientFvPatchScalarField>(tbf[patchi])
-         || isA<fixedGradientFvPatchScalarField>(tbf[patchi])
-        )
-        {
-            hbt[patchi] = gradientEnthalpyRealFluids::typeName;
-        }
-        else if (isA<mixedFvPatchScalarField>(tbf[patchi]))
-        {
-            hbt[patchi] = mixedEnthalpyRealFluids::typeName;
-        }
-    }
-
-    return hbt;
-}
-
-
-void Foam::basicThermo::hRealBoundaryCorrection(volScalarField& h)
-{
-    volScalarField::GeometricBoundaryField& hbf = h.boundaryField();
-
-    forAll(hbf, patchi)
-    {
-        if (isA<gradientEnthalpyRealFluids>(hbf[patchi]))
-        {
-            refCast<gradientEnthalpyRealFluids>(hbf[patchi]).gradient()
-                = hbf[patchi].fvPatchField::snGrad();
-        }
-        else if (isA<mixedEnthalpyRealFluids>(hbf[patchi]))
-        {
-            refCast<mixedEnthalpyRealFluids>(hbf[patchi]).refGrad()
-                = hbf[patchi].fvPatchField::snGrad();
-        }
-    }
-}
-
-
-Foam::wordList Foam::basicThermo::eRealBoundaryTypes()
-{
-    const volScalarField::GeometricBoundaryField& tbf = T_.boundaryField();
-
-    wordList ebt = tbf.types();
-
-    forAll(tbf, patchi)
-    {
-        if (isA<fixedValueFvPatchScalarField>(tbf[patchi]))
-        {
-            ebt[patchi] = fixedInternalEnergyRealFluids::typeName;
-        }
-        else if
-        (
-            isA<zeroGradientFvPatchScalarField>(tbf[patchi])
-         || isA<fixedGradientFvPatchScalarField>(tbf[patchi])
-        )
-        {
-            ebt[patchi] = gradientInternalEnergyRealFluids::typeName;
-        }
-        else if (isA<mixedFvPatchScalarField>(tbf[patchi]))
-        {
-            ebt[patchi] = mixedInternalEnergyRealFluids::typeName;
-        }
-    }
-
-    return ebt;
-}
-
-
-void Foam::basicThermo::eRealBoundaryCorrection(volScalarField& e)
-{
-    volScalarField::GeometricBoundaryField& ebf = e.boundaryField();
-
-    forAll(ebf, patchi)
-    {
-        if (isA<gradientInternalEnergyRealFluids>(ebf[patchi]))
-        {
-            refCast<gradientInternalEnergyRealFluids>(ebf[patchi])
-                .gradient() = ebf[patchi].fvPatchField::snGrad();
-        }
-        else if (isA<mixedInternalEnergyRealFluids>(ebf[patchi]))
-        {
-            refCast<mixedInternalEnergyRealFluids>(ebf[patchi])
-                .refGrad() = ebf[patchi].fvPatchField::snGrad();
-        }
-    }
-}
 
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //

src/thermophysicalModels/basic/derivedFvPatchFields/fixedEnthalpy/fixedEnthalpyFvPatchScalarField.C

@@ -103,12 +103,16 @@ void Foam::fixedEnthalpyFvPatchScalarField::updateCoeffs()
         const_cast<fvPatchScalarField&>(thermo.T().boundaryField()[patchi]);
     Tw.evaluate();
 
+    fvPatchScalarField& pw =
+        const_cast<fvPatchScalarField&>(thermo.p().boundaryField()[patchi]);
+    pw.evaluate();
+
     if
     (
         dimensionedInternalField().name() == db().mangleFileName("h")
     )
     {
-        operator==(thermo.h(Tw, patchi));
+        operator==(pw,thermo.h(Tw, patchi));
     }
     else if
     (
@@ -118,7 +122,7 @@ void Foam::fixedEnthalpyFvPatchScalarField::updateCoeffs()
     }
     else
     {
-        operator==(thermo.hs(Tw, patchi));
+        operator==(pw,thermo.hs(Tw, patchi));
     }
 
     fixedValueFvPatchScalarField::updateCoeffs();

src/thermophysicalModels/basic/derivedFvPatchFields/fixedInternalEnergy/fixedInternalEnergyFvPatchScalarField.C

@@ -107,7 +107,12 @@ void fixedInternalEnergyFvPatchScalarField::updateCoeffs()
     fvPatchScalarField& Tw =
         const_cast<fvPatchScalarField&>(thermo.T().boundaryField()[patchi]);
     Tw.evaluate();
-    operator==(thermo.e(Tw, patchi));
+
+    fvPatchScalarField& pw = 
+        const_cast<fvPatchScalarField&>(thermo.p().boundaryField()[patchi]);
+    pw.evaluate();
+
+    operator==(thermo.e(pw,Tw, patchi));
 
     fixedValueFvPatchScalarField::updateCoeffs();
 }

src/thermophysicalModels/basic/derivedFvPatchFields/gradientEnthalpy/gradientEnthalpyFvPatchScalarField.C

@@ -104,6 +104,11 @@ void Foam::gradientEnthalpyFvPatchScalarField::updateCoeffs()
 
     Tw.evaluate();
 
+    fvPatchScalarField& pw =
+        const_cast<fvPatchScalarField&>(thermo.p().boundaryField()[patchi]);
+
+    pw.evaluate();
+
     if
     (
         dimensionedInternalField().name() == db().mangleFileName("h")
@@ -112,8 +117,8 @@ void Foam::gradientEnthalpyFvPatchScalarField::updateCoeffs()
         gradient() = thermo.Cp(Tw, patchi)*Tw.snGrad()
           + patch().deltaCoeffs()*
             (
-                thermo.h(Tw, patchi)
+                thermo.h(pw, Tw, patchi)
-              - thermo.h(Tw, patch().faceCells())
+              - thermo.h(pw, Tw, patch().faceCells())
             );
     }
     else if
@@ -127,8 +132,8 @@ void Foam::gradientEnthalpyFvPatchScalarField::updateCoeffs()
         gradient() = thermo.Cp(Tw, patchi)*Tw.snGrad()
         + patch().deltaCoeffs()*
         (
-            thermo.hs(Tw, patchi)
+            thermo.hs(pw, Tw, patchi)
-          - thermo.hs(Tw, patch().faceCells())
+          - thermo.hs(pw, Tw, patch().faceCells())
         );
     }
 

src/thermophysicalModels/basic/derivedFvPatchFields/gradientInternalEnergy/gradientInternalEnergyFvPatchScalarField.C

@@ -114,11 +114,15 @@ void gradientInternalEnergyFvPatchScalarField::updateCoeffs()
 
     Tw.evaluate();
 
-    gradient() = thermo.Cv(Tw, patchi)*Tw.snGrad()
+    fvPatchScalarField& pw =
+        const_cast<fvPatchScalarField&>(thermo.p().boundaryField()[patchi]);
+    pw.evaluate();
+
+    gradient() = thermo.Cv(pw,Tw, patchi)*Tw.snGrad()
       + patch().deltaCoeffs()*
         (
-            thermo.e(Tw, patchi)
+            thermo.e(pw,Tw, patchi)
-          - thermo.e(Tw, patch().faceCells())
+          - thermo.e(pw,Tw, patch().faceCells())
         );
 
     fixedGradientFvPatchScalarField::updateCoeffs();

src/thermophysicalModels/basic/derivedFvPatchFields/mixedEnthalpy/mixedEnthalpyFvPatchScalarField.C

@@ -1,25 +1,26 @@
 /*---------------------------------------------------------------------------*\
   =========                 |
-  \\      /  F ield         | foam-extend: Open Source CFD
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     | Version:     3.2
+   \\    /   O peration     |
-    \\  /    A nd           | Web:         http://www.foam-extend.org
+    \\  /    A nd           | Copyright held by original author
-     \\/     M anipulation  | For copyright notice see file Copyright
+     \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
-    This file is part of foam-extend.
+    This file is part of OpenFOAM.
 
-    foam-extend is free software: you can redistribute it and/or modify it
+    OpenFOAM is free software; you can redistribute it and/or modify it
     under the terms of the GNU General Public License as published by the
-    Free Software Foundation, either version 3 of the License, or (at your
+    Free Software Foundation; either version 2 of the License, or (at your
     option) any later version.
 
-    foam-extend is distributed in the hope that it will be useful, but
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
-    WITHOUT ANY WARRANTY; without even the implied warranty of
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
-    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
-    General Public License for more details.
+    for more details.
 
     You should have received a copy of the GNU General Public License
-    along with foam-extend.  If not, see <http://www.gnu.org/licenses/>.
+    along with OpenFOAM; if not, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
 
 \*---------------------------------------------------------------------------*/
 
@@ -110,6 +111,10 @@ void Foam::mixedEnthalpyFvPatchScalarField::updateCoeffs()
 
     Tw.evaluate();
 
+    fvPatchScalarField& pw =
+        const_cast<fvPatchScalarField&>(thermo.p().boundaryField()[patchi]);
+    pw.evaluate();
+
     valueFraction() = Tw.valueFraction();
 
     if
@@ -121,8 +126,8 @@ void Foam::mixedEnthalpyFvPatchScalarField::updateCoeffs()
         refGrad() = thermo.Cp(Tw, patchi)*Tw.refGrad()
           + patch().deltaCoeffs()*
             (
-                thermo.h(Tw, patchi)
+                thermo.h(pw, Tw, patchi)
-              - thermo.h(Tw, patch().faceCells())
+              - thermo.h(pw, Tw, patch().faceCells())
             );
     }
     else if
@@ -137,8 +142,8 @@ void Foam::mixedEnthalpyFvPatchScalarField::updateCoeffs()
         refGrad() = thermo.Cp(Tw, patchi)*Tw.refGrad()
           + patch().deltaCoeffs()*
             (
-                thermo.hs(Tw, patchi)
+                thermo.hs(pw, Tw, patchi)
-              - thermo.hs(Tw, patch().faceCells())
+              - thermo.hs(pw, Tw, patch().faceCells())
             );
     }
 

src/thermophysicalModels/basic/derivedFvPatchFields/mixedInternalEnergy/mixedInternalEnergyFvPatchScalarField.C

@@ -115,13 +115,17 @@ void mixedInternalEnergyFvPatchScalarField::updateCoeffs()
 
     Tw.evaluate();
 
+    fvPatchScalarField& pw =
+        const_cast<fvPatchScalarField&>(thermo.p().boundaryField()[patchi]);
+    pw.evaluate();
+
     valueFraction() = Tw.valueFraction();
-    refValue() = thermo.e(Tw.refValue(), patchi);
+    refValue() = thermo.e(pw,Tw.refValue(), patchi);
-    refGrad() = thermo.Cv(Tw, patchi)*Tw.refGrad()
+    refGrad() = thermo.Cv(pw,Tw, patchi)*Tw.refGrad()
       + patch().deltaCoeffs()*
         (
-            thermo.e(Tw, patchi)
+            thermo.e(pw,Tw, patchi)
-          - thermo.e(Tw, patch().faceCells())
+          - thermo.e(pw,Tw, patch().faceCells())
         );
 
     mixedFvPatchScalarField::updateCoeffs();

src/thermophysicalModels/basic/derivedFvPatchFields/wallHeatTransfer/wallHeatTransferFvPatchScalarField.C

@@ -157,7 +157,12 @@ void Foam::wallHeatTransferFvPatchScalarField::updateCoeffs()
     const label patchi = patch().index();
 
     const scalarField& Tw = thermo.T().boundaryField()[patchi];
-    scalarField Cpw = thermo.Cp(Tw, patchi);
+
+    fvPatchScalarField& pw =
+        const_cast<fvPatchScalarField&>(thermo.p().boundaryField()[patchi]);
+    pw.evaluate();
+
+    scalarField Cpw = thermo.Cp(pw,Tw, patchi);
 
     valueFraction() =
         1.0/

src/thermophysicalModels/basic/derivedFvPatchFieldsRealGas/fixedEnthalpyRealFluids/fixedEnthalpyRealFluids.C

@@ -1,147 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     |
-    \\  /    A nd           | 
-     \\/     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
-
-Modified by
-Christian Lucas
-Institut für Thermodynamik
-Technische Universität Braunschweig 
-Germany
-
-
-\*---------------------------------------------------------------------------*/
-
-#include "fixedEnthalpyRealFluids.H"
-#include "addToRunTimeSelectionTable.H"
-#include "fvPatchFieldMapper.H"
-#include "volFields.H"
-#include "basicPsiThermo.H"
-
-// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
-
-Foam::fixedEnthalpyRealFluids::fixedEnthalpyRealFluids
-(
-    const fvPatch& p,
-    const DimensionedField<scalar, volMesh>& iF
-)
-:
-    fixedValueFvPatchScalarField(p, iF)
-{}
-
-
-Foam::fixedEnthalpyRealFluids::fixedEnthalpyRealFluids
-(
-    const fixedEnthalpyRealFluids& ptf,
-    const fvPatch& p,
-    const DimensionedField<scalar, volMesh>& iF,
-    const fvPatchFieldMapper& mapper
-)
-:
-    fixedValueFvPatchScalarField(ptf, p, iF, mapper)
-{}
-
-
-Foam::fixedEnthalpyRealFluids::fixedEnthalpyRealFluids
-(
-    const fvPatch& p,
-    const DimensionedField<scalar, volMesh>& iF,
-    const dictionary& dict
-)
-:
-    fixedValueFvPatchScalarField(p, iF, dict)
-{}
-
-
-Foam::fixedEnthalpyRealFluids::fixedEnthalpyRealFluids
-(
-    const fixedEnthalpyRealFluids& tppsf
-)
-:
-    fixedValueFvPatchScalarField(tppsf)
-{}
-
-
-Foam::fixedEnthalpyRealFluids::fixedEnthalpyRealFluids
-(
-    const fixedEnthalpyRealFluids& tppsf,
-    const DimensionedField<scalar, volMesh>& iF
-)
-:
-    fixedValueFvPatchScalarField(tppsf, iF)
-{}
-
-
-// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
-
-void Foam::fixedEnthalpyRealFluids::updateCoeffs()
-{
-    if (updated())
-    {
-        return;
-    }
-
-    const basicPsiThermo& thermo = db().lookupObject<basicPsiThermo>
-    (
-        "thermophysicalProperties"
-    );
-
-    const label patchi = patch().index();
-
-    fvPatchScalarField& Tw =
-        const_cast<fvPatchScalarField&>(thermo.T().boundaryField()[patchi]);
-    Tw.evaluate();
-
-    fvPatchScalarField& pw =
-        const_cast<fvPatchScalarField&>(thermo.p().boundaryField()[patchi]);
-    pw.evaluate();
-
-
-    if (dimensionedInternalField().name() == "h")
-    {
-   
-        operator==(thermo.h(pw,Tw, patchi));
-    }
-    else
-    {
-        operator==(thermo.hs(Tw, patchi));
-    }
-
-    fixedValueFvPatchScalarField::updateCoeffs();
-
-    
-}
-
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-namespace Foam
-{
-    makePatchTypeField
-    (
-        fvPatchScalarField,
-        fixedEnthalpyRealFluids
-    );
-}
-
-
-// ************************************************************************* //

src/thermophysicalModels/basic/derivedFvPatchFieldsRealGas/fixedEnthalpyRealFluids/fixedEnthalpyRealFluids.H

@@ -1,146 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     |
-    \\  /    A nd           | 
-     \\/     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::fixedEnthalpyRealFluids
-
-Description
-    A fixed boundary condition for enthalpy
-
-SourceFiles
-    fixedEnthalpyRealFluids.C
-
-Modified by
-Christian Lucas
-Institut für Thermodynamik
-Technische Universität Braunschweig 
-Germany
-
-\*---------------------------------------------------------------------------*/
-
-#ifndef fixedEnthalpyRealFluids_H
-#define fixedEnthalpyRealFluids_H
-
-#include "fixedValueFvPatchFields.H"
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-namespace Foam
-{
-
-/*---------------------------------------------------------------------------*\
-                Class fixedEnthalpyRealFluids Declaration
-\*---------------------------------------------------------------------------*/
-
-class fixedEnthalpyRealFluids
-:
-    public fixedValueFvPatchScalarField
-{
-
-public:
-
-    //- Runtime type information
-    TypeName("fixedEnthalpy");
-
-
-    // Constructors
-
-        //- Construct from patch and internal field
-        fixedEnthalpyRealFluids
-        (
-            const fvPatch&,
-            const DimensionedField<scalar, volMesh>&
-        );
-
-        //- Construct from patch, internal field and dictionary
-        fixedEnthalpyRealFluids
-        (
-            const fvPatch&,
-            const DimensionedField<scalar, volMesh>&,
-            const dictionary&
-        );
-
-        //- Construct by mapping given fixedEnthalpyFvPatchScalarField
-        // onto a new patch
-        fixedEnthalpyRealFluids
-        (
-            const fixedEnthalpyRealFluids&,
-            const fvPatch&,
-            const DimensionedField<scalar, volMesh>&,
-            const fvPatchFieldMapper&
-        );
-
-        //- Construct as copy
-        fixedEnthalpyRealFluids
-        (
-            const fixedEnthalpyRealFluids&
-        );
-
-        //- Construct and return a clone
-        virtual tmp<fvPatchScalarField> clone() const
-        {
-            return tmp<fvPatchScalarField>
-            (
-                new fixedEnthalpyRealFluids(*this)
-            );
-        }
-
-        //- Construct as copy setting internal field reference
-        fixedEnthalpyRealFluids
-        (
-            const fixedEnthalpyRealFluids&,
-            const DimensionedField<scalar, volMesh>&
-        );
-
-        //- Construct and return a clone setting internal field reference
-        virtual tmp<fvPatchScalarField> clone
-        (
-            const DimensionedField<scalar, volMesh>& iF
-        ) const
-        {
-            return tmp<fvPatchScalarField>
-            (
-                new fixedEnthalpyRealFluids(*this, iF)
-            );
-        }
-
-
-    // Member functions
-
-        // Evaluation functions
-
-            //- Update the coefficients associated with the patch field
-            virtual void updateCoeffs();
-};
-
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-} // End namespace Foam
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-#endif
-
-// ************************************************************************* //

src/thermophysicalModels/basic/derivedFvPatchFieldsRealGas/fixedInternalEnergyRealFluids/fixedInternalEnergyRealFluids.C

@@ -1,137 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     |
-    \\  /    A nd           | 
-     \\/     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
-
-Modified by
-Christian Lucas
-Institut für Thermodynamik
-Technische Universität Braunschweig 
-Germany
-
-\*---------------------------------------------------------------------------*/
-
-#include "fixedInternalEnergyRealFluids.H"
-#include "addToRunTimeSelectionTable.H"
-#include "fvPatchFieldMapper.H"
-#include "volFields.H"
-#include "basicPsiThermo.H"
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-namespace Foam
-{
-
-// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
-
-fixedInternalEnergyRealFluids::fixedInternalEnergyRealFluids
-(
-    const fvPatch& p,
-    const DimensionedField<scalar, volMesh>& iF
-)
-:
-    fixedValueFvPatchScalarField(p, iF)
-{}
-
-
-fixedInternalEnergyRealFluids::fixedInternalEnergyRealFluids
-(
-    const fixedInternalEnergyRealFluids& ptf,
-    const fvPatch& p,
-    const DimensionedField<scalar, volMesh>& iF,
-    const fvPatchFieldMapper& mapper
-)
-:
-    fixedValueFvPatchScalarField(ptf, p, iF, mapper)
-{}
-
-
-fixedInternalEnergyRealFluids::fixedInternalEnergyRealFluids
-(
-    const fvPatch& p,
-    const DimensionedField<scalar, volMesh>& iF,
-    const dictionary& dict
-)
-:
-    fixedValueFvPatchScalarField(p, iF, dict)
-{}
-
-
-fixedInternalEnergyRealFluids::fixedInternalEnergyRealFluids
-(
-    const fixedInternalEnergyRealFluids& tppsf
-)
-:
-    fixedValueFvPatchScalarField(tppsf)
-{}
-
-
-fixedInternalEnergyRealFluids::fixedInternalEnergyRealFluids
-(
-    const fixedInternalEnergyRealFluids& tppsf,
-    const DimensionedField<scalar, volMesh>& iF
-)
-:
-    fixedValueFvPatchScalarField(tppsf, iF)
-{}
-
-
-// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
-
-void fixedInternalEnergyRealFluids::updateCoeffs()
-{
-    if (updated())
-    {
-        return;
-    }
-
-    const basicPsiThermo& thermo = db().lookupObject<basicPsiThermo>
-    (
-        "thermophysicalProperties"
-    );
-
-    
-    const label patchi = patch().index();
-
-    fvPatchScalarField& Tw = 
-        const_cast<fvPatchScalarField&>(thermo.T().boundaryField()[patchi]);
-    Tw.evaluate();
-
-    fvPatchScalarField& pw =
-        const_cast<fvPatchScalarField&>(thermo.p().boundaryField()[patchi]);
-    pw.evaluate();
-
-    operator==(thermo.e(pw,Tw, patchi));
-
-    fixedValueFvPatchScalarField::updateCoeffs();
-}
-
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-makePatchTypeField(fvPatchScalarField, fixedInternalEnergyRealFluids);
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-} // End namespace Foam
-
-// ************************************************************************* //

src/thermophysicalModels/basic/derivedFvPatchFieldsRealGas/fixedInternalEnergyRealFluids/fixedInternalEnergyRealFluids.H

@@ -1,146 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     |
-    \\  /    A nd           | 
-     \\/     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::fixedInternalEnergyRealFluids
-
-Description
-    A fixed boundary condition for internal energy
-
-SourceFiles
-    fixedInternalEnergyRealFluids.C
-
-Modified by
-Christian Lucas
-Institut für Thermodynamik
-Technische Universität Braunschweig 
-Germany
-
-\*---------------------------------------------------------------------------*/
-
-#ifndef fixedInternalEnergyRealFluids_H
-#define fixedInternalEnergyRealFluids_H
-
-#include "fixedValueFvPatchFields.H"
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-namespace Foam
-{
-
-/*---------------------------------------------------------------------------*\
-                Class fixedInternalEnergyRealFluids Declaration
-\*---------------------------------------------------------------------------*/
-
-class fixedInternalEnergyRealFluids
-:
-    public fixedValueFvPatchScalarField
-{
-
-public:
-
-    //- Runtime type information
-    TypeName("fixedInternalEnergy");
-
-
-    // Constructors
-
-        //- Construct from patch and internal field
-        fixedInternalEnergyRealFluids
-        (
-            const fvPatch&,
-            const DimensionedField<scalar, volMesh>&
-        );
-
-        //- Construct from patch, internal field and dictionary
-        fixedInternalEnergyRealFluids
-        (
-            const fvPatch&,
-            const DimensionedField<scalar, volMesh>&,
-            const dictionary&
-        );
-
-        //- Construct by mapping given fixedInternalEnergyRealFluids
-        // onto a new patch
-        fixedInternalEnergyRealFluids
-        (
-            const fixedInternalEnergyRealFluids&,
-            const fvPatch&,
-            const DimensionedField<scalar, volMesh>&,
-            const fvPatchFieldMapper&
-        );
-
-        //- Construct as copy
-        fixedInternalEnergyRealFluids
-        (
-            const fixedInternalEnergyRealFluids&
-        );
-
-        //- Construct and return a clone
-        virtual tmp<fvPatchScalarField> clone() const
-        {
-            return tmp<fvPatchScalarField>
-            (
-                new fixedInternalEnergyRealFluids(*this)
-            );
-        }
-
-        //- Construct as copy setting internal field reference
-        fixedInternalEnergyRealFluids
-        (
-            const fixedInternalEnergyRealFluids&,
-            const DimensionedField<scalar, volMesh>&
-        );
-
-        //- Construct and return a clone setting internal field reference
-        virtual tmp<fvPatchScalarField> clone
-        (
-            const DimensionedField<scalar, volMesh>& iF
-        ) const
-        {
-            return tmp<fvPatchScalarField>
-            (
-                new fixedInternalEnergyRealFluids(*this, iF)
-            );
-        }
-
-
-    // Member functions
-
-        // Evaluation functions
-
-            //- Update the coefficients associated with the patch field
-            virtual void updateCoeffs();
-};
-
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-} // End namespace Foam
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-#endif
-
-// ************************************************************************* //

src/thermophysicalModels/basic/derivedFvPatchFieldsRealGas/gradientEnthalpyRealFluids/gradientEnthalpyRealFluids.C

@@ -1,154 +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
-
-Modified by
-Christian Lucas
-Institut für Thermodynamik
-Technische Universität Braunschweig 
-Germany
-
-\*---------------------------------------------------------------------------*/
-
-#include "gradientEnthalpyRealFluids.H"
-#include "addToRunTimeSelectionTable.H"
-#include "fvPatchFieldMapper.H"
-#include "volFields.H"
-#include "basicPsiThermo.H"
-
-// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
-
-Foam::gradientEnthalpyRealFluids::gradientEnthalpyRealFluids
-(
-    const fvPatch& p,
-    const DimensionedField<scalar, volMesh>& iF
-)
-:
-    fixedGradientFvPatchScalarField(p, iF)
-{}
-
-
-Foam::gradientEnthalpyRealFluids::gradientEnthalpyRealFluids
-(
-    const gradientEnthalpyRealFluids& ptf,
-    const fvPatch& p,
-    const DimensionedField<scalar, volMesh>& iF,
-    const fvPatchFieldMapper& mapper
-)
-:
-    fixedGradientFvPatchScalarField(ptf, p, iF, mapper)
-{}
-
-
-Foam::gradientEnthalpyRealFluids::gradientEnthalpyRealFluids
-(
-    const fvPatch& p,
-    const DimensionedField<scalar, volMesh>& iF,
-    const dictionary& dict
-)
-:
-    fixedGradientFvPatchScalarField(p, iF, dict)
-{}
-
-
-Foam::gradientEnthalpyRealFluids::gradientEnthalpyRealFluids
-(
-    const gradientEnthalpyRealFluids& tppsf
-)
-:
-    fixedGradientFvPatchScalarField(tppsf)
-{}
-
-
-Foam::gradientEnthalpyRealFluids::gradientEnthalpyRealFluids
-(
-    const gradientEnthalpyRealFluids& tppsf,
-    const DimensionedField<scalar, volMesh>& iF
-)
-:
-    fixedGradientFvPatchScalarField(tppsf, iF)
-{}
-
-
-// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
-
-void Foam::gradientEnthalpyRealFluids::updateCoeffs()
-{
-    if (updated())
-    {
-        return;
-    }
-
-    const basicPsiThermo& thermo = db().lookupObject<basicPsiThermo>
-    (
-        "thermophysicalProperties"
-    );
-
-    const label patchi = patch().index();
-
-    fvPatchScalarField& pw =
-        const_cast<fvPatchScalarField&>(thermo.p().boundaryField()[patchi]);
-
-    pw.evaluate();
-
-    fvPatchScalarField& Tw =
-        const_cast<fvPatchScalarField&>(thermo.T().boundaryField()[patchi]);
-
-    Tw.evaluate();
-
-    if (dimensionedInternalField().name() == "h")
-    {
-        gradient() = thermo.Cp(pw,Tw, patchi)*Tw.snGrad()
-        + patch().deltaCoeffs()*
-        (
-            thermo.h(pw,Tw, patchi)
-          - thermo.h(pw,Tw, patch().faceCells())
-        );
-    }
- /*   else
-    {
-        gradient() = thermo.Cp(Tw, patchi)*Tw.snGrad()
-        + patch().deltaCoeffs()*
-        (
-            thermo.hs(Tw, patchi)
-          - thermo.hs(Tw, patch().faceCells())
-        );
-    }
-*/
-    fixedGradientFvPatchScalarField::updateCoeffs();
-}
-
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-namespace Foam
-{
-    makePatchTypeField
-    (
-        fvPatchScalarField,
-        gradientEnthalpyRealFluids
-    );
-}
-
-
-// ************************************************************************* //

src/thermophysicalModels/basic/derivedFvPatchFieldsRealGas/gradientEnthalpyRealFluids/gradientEnthalpyRealFluids.H

@@ -1,146 +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
-
-Class
-    Foam::gradientEnthalpyRealFluids
-
-Description
-    Gradient boundary condition for enthalpy
-
-SourceFiles
-    gradientEnthalpyRealFluids.C
-
-Modified by
-Christian Lucas
-Institut für Thermodynamik
-Technische Universität Braunschweig 
-Germany
-
-\*---------------------------------------------------------------------------*/
-
-#ifndef gradientEnthalpyRealFluids_H
-#define gradientEnthalpyRealFluids_H
-
-#include "fixedGradientFvPatchFields.H"
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-namespace Foam
-{
-
-/*---------------------------------------------------------------------------*\
-               Class gradientEnthalpyRealFluids Declaration
-\*---------------------------------------------------------------------------*/
-
-class gradientEnthalpyRealFluids
-:
-    public fixedGradientFvPatchScalarField
-{
-
-public:
-
-    //- Runtime type information
-    TypeName("gradientEnthalpy");
-
-
-    // Constructors
-
-        //- Construct from patch and internal field
-        gradientEnthalpyRealFluids
-        (
-            const fvPatch&,
-            const DimensionedField<scalar, volMesh>&
-        );
-
-        //- Construct from patch, internal field and dictionary
-        gradientEnthalpyRealFluids
-        (
-            const fvPatch&,
-            const DimensionedField<scalar, volMesh>&,
-            const dictionary&
-        );
-
-        //- Construct by mapping given gradientEnthalpyRealFluids
-        // onto a new patch
-        gradientEnthalpyRealFluids
-        (
-            const gradientEnthalpyRealFluids&,
-            const fvPatch&,
-            const DimensionedField<scalar, volMesh>&,
-            const fvPatchFieldMapper&
-        );
-
-        //- Construct as copy
-        gradientEnthalpyRealFluids
-        (
-            const gradientEnthalpyRealFluids&
-        );
-
-        //- Construct and return a clone
-        virtual tmp<fvPatchScalarField> clone() const
-        {
-            return tmp<fvPatchScalarField>
-            (
-                new gradientEnthalpyRealFluids(*this)
-            );
-        }
-
-        //- Construct as copy setting internal field reference
-        gradientEnthalpyRealFluids
-        (
-            const gradientEnthalpyRealFluids&,
-            const DimensionedField<scalar, volMesh>&
-        );
-
-        //- Construct and return a clone setting internal field reference
-        virtual tmp<fvPatchScalarField> clone
-        (
-            const DimensionedField<scalar, volMesh>& iF
-        ) const
-        {
-            return tmp<fvPatchScalarField>
-            (
-                new gradientEnthalpyRealFluids(*this, iF)
-            );
-        }
-
-
-    // Member functions
-
-        // Evaluation functions
-
-            //- Update the coefficients associated with the patch field
-            virtual void updateCoeffs();
-};
-
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-} // End namespace Foam
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-#endif
-
-// ************************************************************************* //

src/thermophysicalModels/basic/derivedFvPatchFieldsRealGas/gradientInternalEnergyRealFluids/gradientInternalEnergyRealFluids.C

@@ -1,146 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     |
-    \\  /    A nd           | 
-     \\/     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
-
-Modified by
-Christian Lucas
-Institut für Thermodynamik
-Technische Universität Braunschweig 
-Germany
-
-\*---------------------------------------------------------------------------*/
-
-#include "gradientInternalEnergyRealFluids.H"
-#include "addToRunTimeSelectionTable.H"
-#include "fvPatchFieldMapper.H"
-#include "volFields.H"
-#include "basicPsiThermo.H"
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-namespace Foam
-{
-
-// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
-
-gradientInternalEnergyRealFluids::gradientInternalEnergyRealFluids
-(
-    const fvPatch& p,
-    const DimensionedField<scalar, volMesh>& iF
-)
-:
-    fixedGradientFvPatchScalarField(p, iF)
-{}
-
-
-gradientInternalEnergyRealFluids::gradientInternalEnergyRealFluids
-(
-    const gradientInternalEnergyRealFluids& ptf,
-    const fvPatch& p,
-    const DimensionedField<scalar, volMesh>& iF,
-    const fvPatchFieldMapper& mapper
-)
-:
-    fixedGradientFvPatchScalarField(ptf, p, iF, mapper)
-{}
-
-
-gradientInternalEnergyRealFluids::gradientInternalEnergyRealFluids
-(
-    const fvPatch& p,
-    const DimensionedField<scalar, volMesh>& iF,
-    const dictionary& dict
-)
-:
-    fixedGradientFvPatchScalarField(p, iF, dict)
-{}
-
-
-gradientInternalEnergyRealFluids::gradientInternalEnergyRealFluids
-(
-    const gradientInternalEnergyRealFluids& tppsf
-)
-:
-    fixedGradientFvPatchScalarField(tppsf)
-{}
-
-
-gradientInternalEnergyRealFluids::gradientInternalEnergyRealFluids
-(
-    const gradientInternalEnergyRealFluids& tppsf,
-    const DimensionedField<scalar, volMesh>& iF
-)
-:
-    fixedGradientFvPatchScalarField(tppsf, iF)
-{}
-
-
-// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
-
-void gradientInternalEnergyRealFluids::updateCoeffs()
-{
-    if (updated())
-    {
-        return;
-    }
-
-    const basicPsiThermo& thermo = db().lookupObject<basicPsiThermo>
-    (
-        "thermophysicalProperties"
-    );
-
-    
-    const label patchi = patch().index();
-
-    fvPatchScalarField& Tw = 
-        const_cast<fvPatchScalarField&>(thermo.T().boundaryField()[patchi]);
-
-    Tw.evaluate();
-
-    fvPatchScalarField& pw =
-        const_cast<fvPatchScalarField&>(thermo.p().boundaryField()[patchi]);
-
-    pw.evaluate();
-
-
-
-    gradient() = thermo.Cv(pw,Tw, patchi)*Tw.snGrad()
-      + patch().deltaCoeffs()*
-        (
-            thermo.e(pw,Tw, patchi)
-          - thermo.e(pw,Tw, patch().faceCells())
-        );
-
-    fixedGradientFvPatchScalarField::updateCoeffs();
-}
-
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-makePatchTypeField(fvPatchScalarField, gradientInternalEnergyRealFluids);
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-} // End namespace Foam
-
-// ************************************************************************* //

src/thermophysicalModels/basic/derivedFvPatchFieldsRealGas/gradientInternalEnergyRealFluids/gradientInternalEnergyRealFluids.H

@@ -1,146 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     |
-    \\  /    A nd           | 
-     \\/     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::gradientInternalEnergyRealFluids
-
-Description
-    Gradient boundary condition for internal energy
-
-SourceFiles
-    gradientInternalEnergyRealFluids.C
-
-Modified by
-Christian Lucas
-Institut für Thermodynamik
-Technische Universität Braunschweig 
-Germany
-
-\*---------------------------------------------------------------------------*/
-
-#ifndef gradientInternalEnergyRealFluids_H
-#define gradientInternalEnergyRealFluids_H
-
-#include "fixedGradientFvPatchFields.H"
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-namespace Foam
-{
-
-/*---------------------------------------------------------------------------*\
-               Class gradientInternalEnergyRealFluids Declaration
-\*---------------------------------------------------------------------------*/
-
-class gradientInternalEnergyRealFluids
-:
-    public fixedGradientFvPatchScalarField
-{
-
-public:
-
-    //- Runtime type information
-    TypeName("gradientInternalEnergy");
-
-
-    // Constructors
-
-        //- Construct from patch and internal field
-        gradientInternalEnergyRealFluids
-        (
-            const fvPatch&,
-            const DimensionedField<scalar, volMesh>&
-        );
-
-        //- Construct from patch, internal field and dictionary
-        gradientInternalEnergyRealFluids
-        (
-            const fvPatch&,
-            const DimensionedField<scalar, volMesh>&,
-            const dictionary&
-        );
-
-        //- Construct by mapping given gradientInternalEnergyRealFluids
-        // onto a new patch
-        gradientInternalEnergyRealFluids
-        (
-            const gradientInternalEnergyRealFluids&,
-            const fvPatch&,
-            const DimensionedField<scalar, volMesh>&,
-            const fvPatchFieldMapper&
-        );
-
-        //- Construct as copy
-        gradientInternalEnergyRealFluids
-        (
-            const gradientInternalEnergyRealFluids&
-        );
-
-        //- Construct and return a clone
-        virtual tmp<fvPatchScalarField> clone() const
-        {
-            return tmp<fvPatchScalarField>
-            (
-                new gradientInternalEnergyRealFluids(*this)
-            );
-        }
-
-        //- Construct as copy setting internal field reference
-        gradientInternalEnergyRealFluids
-        (
-            const gradientInternalEnergyRealFluids&,
-            const DimensionedField<scalar, volMesh>&
-        );
-
-        //- Construct and return a clone setting internal field reference
-        virtual tmp<fvPatchScalarField> clone
-        (
-            const DimensionedField<scalar, volMesh>& iF
-        ) const
-        {
-            return tmp<fvPatchScalarField>
-            (
-                new gradientInternalEnergyRealFluids(*this, iF)
-            );
-        }
-
-
-    // Member functions
-
-        // Evaluation functions
-
-            //- Update the coefficients associated with the patch field
-            virtual void updateCoeffs();
-};
-
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-} // End namespace Foam
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-#endif
-
-// ************************************************************************* //

src/thermophysicalModels/basic/derivedFvPatchFieldsRealGas/mixedEnthalpyRealFluids/mixedEnthalpyRealFluids.C

@@ -1,164 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     |
-    \\  /    A nd           | 
-     \\/     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
-
-Modified by
-Christian Lucas
-Institut für Thermodynamik
-Technische Universität Braunschweig 
-Germany
-
-\*---------------------------------------------------------------------------*/
-
-#include "mixedEnthalpyRealFluids.H"
-#include "addToRunTimeSelectionTable.H"
-#include "fvPatchFieldMapper.H"
-#include "volFields.H"
-#include "basicPsiThermo.H"
-
-// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
-
-Foam::mixedEnthalpyRealFluids::mixedEnthalpyRealFluids
-(
-    const fvPatch& p,
-    const DimensionedField<scalar, volMesh>& iF
-)
-:
-    mixedFvPatchScalarField(p, iF)
-{
-    valueFraction() = 0.0;
-    refValue() = 0.0;
-    refGrad() = 0.0;
-}
-
-
-Foam::mixedEnthalpyRealFluids::mixedEnthalpyRealFluids
-(
-    const mixedEnthalpyRealFluids& ptf,
-    const fvPatch& p,
-    const DimensionedField<scalar, volMesh>& iF,
-    const fvPatchFieldMapper& mapper
-)
-:
-    mixedFvPatchScalarField(ptf, p, iF, mapper)
-{}
-
-
-Foam::mixedEnthalpyRealFluids::mixedEnthalpyRealFluids
-(
-    const fvPatch& p,
-    const DimensionedField<scalar, volMesh>& iF,
-    const dictionary& dict
-)
-:
-    mixedFvPatchScalarField(p, iF, dict)
-{}
-
-
-Foam::mixedEnthalpyRealFluids::mixedEnthalpyRealFluids
-(
-    const mixedEnthalpyRealFluids& tppsf
-)
-:
-    mixedFvPatchScalarField(tppsf)
-{}
-
-
-Foam::mixedEnthalpyRealFluids::mixedEnthalpyRealFluids
-(
-    const mixedEnthalpyRealFluids& tppsf,
-    const DimensionedField<scalar, volMesh>& iF
-)
-:
-    mixedFvPatchScalarField(tppsf, iF)
-{}
-
-
-// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
-
-void Foam::mixedEnthalpyRealFluids::updateCoeffs()
-{
-    if (updated())
-    {
-        return;
-    }
-
-    const basicPsiThermo& thermo = db().lookupObject<basicPsiThermo>
-    (
-        "thermophysicalProperties"
-    );
-
-
-    const label patchi = patch().index();
-
-    mixedFvPatchScalarField& Tw = refCast<mixedFvPatchScalarField>
-    (
-        const_cast<fvPatchScalarField&>(thermo.T().boundaryField()[patchi])
-    );
-
-    Tw.evaluate();
-
-    fvPatchScalarField& pw =
-        const_cast<fvPatchScalarField&>(thermo.p().boundaryField()[patchi]);
-    pw.evaluate();
-
-    valueFraction() = Tw.valueFraction();
-
-    if (dimensionedInternalField().name() == "h")
-    {
-        refValue() = thermo.h(pw,Tw.refValue(), patchi);
-        refGrad() = thermo.Cp(pw,Tw, patchi)*Tw.refGrad()
-        + patch().deltaCoeffs()*
-         (
-            thermo.h(pw,Tw, patchi)
-          - thermo.h(pw,Tw, patch().faceCells())
-         );
-    }
-  /*  else
-    {
-        refValue() = thermo.hs(Tw.refValue(), patchi);
-        refGrad() = thermo.CpBC(pw,Tw, patchi)*Tw.refGrad()
-        + patch().deltaCoeffs()*
-         (
-            thermo.hs(Tw, patchi)
-          - thermo.hs(Tw, patch().faceCells())
-         );
-    }*/
-
-    mixedFvPatchScalarField::updateCoeffs();
-}
-
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-namespace Foam
-{
-    makePatchTypeField
-    (
-        fvPatchScalarField,
-        mixedEnthalpyRealFluids
-    );
-}
-
-
-// ************************************************************************* //

src/thermophysicalModels/basic/derivedFvPatchFieldsRealGas/mixedEnthalpyRealFluids/mixedEnthalpyRealFluids.H

@@ -1,146 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     |
-    \\  /    A nd           | 
-     \\/     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::mixedEnthalpyRealFluids
-
-Description
-    Mixed boundary conditions for enthalpy
-
-SourceFiles
-    mixedEnthalpyRealFluids.C
-
-Modified by
-Christian Lucas
-Institut für Thermodynamik
-Technische Universität Braunschweig 
-Germany
-
-\*---------------------------------------------------------------------------*/
-
-#ifndef mixedEnthalpyRealFluids_H
-#define mixedEnthalpyRealFluids_H
-
-#include "mixedFvPatchFields.H"
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-namespace Foam
-{
-
-/*---------------------------------------------------------------------------*\
-               Class mixedEnthalpyRealFluids Declaration
-\*---------------------------------------------------------------------------*/
-
-class mixedEnthalpyRealFluids
-:
-    public mixedFvPatchScalarField
-{
-
-public:
-
-    //- Runtime type information
-    TypeName("mixedEnthalpy");
-
-
-    // Constructors
-
-        //- Construct from patch and internal field
-        mixedEnthalpyRealFluids
-        (
-            const fvPatch&,
-            const DimensionedField<scalar, volMesh>&
-        );
-
-        //- Construct from patch, internal field and dictionary
-        mixedEnthalpyRealFluids
-        (
-            const fvPatch&,
-            const DimensionedField<scalar, volMesh>&,
-            const dictionary&
-        );
-
-        //- Construct by mapping given mixedEnthalpyRealFluids
-        // onto a new patch
-        mixedEnthalpyRealFluids
-        (
-            const mixedEnthalpyRealFluids&,
-            const fvPatch&,
-            const DimensionedField<scalar, volMesh>&,
-            const fvPatchFieldMapper&
-        );
-
-        //- Construct as copy
-        mixedEnthalpyRealFluids
-        (
-            const mixedEnthalpyRealFluids&
-        );
-
-        //- Construct and return a clone
-        virtual tmp<fvPatchScalarField> clone() const
-        {
-            return tmp<fvPatchScalarField>
-            (
-                new mixedEnthalpyRealFluids(*this)
-            );
-        }
-
-        //- Construct as copy setting internal field reference
-        mixedEnthalpyRealFluids
-        (
-            const mixedEnthalpyRealFluids&,
-            const DimensionedField<scalar, volMesh>&
-        );
-
-        //- Construct and return a clone setting internal field reference
-        virtual tmp<fvPatchScalarField> clone
-        (
-            const DimensionedField<scalar, volMesh>& iF
-        ) const
-        {
-            return tmp<fvPatchScalarField>
-            (
-                new mixedEnthalpyRealFluids(*this, iF)
-            );
-        }
-
-
-    // Member functions
-
-        // Evaluation functions
-
-            //- Update the coefficients associated with the patch field
-            virtual void updateCoeffs();
-};
-
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-} // End namespace Foam
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-#endif
-
-// ************************************************************************* //

src/thermophysicalModels/basic/derivedFvPatchFieldsRealGas/mixedInternalEnergyRealFluids/mixedInternalEnergyRealFluids.C

@@ -1,151 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     |
-    \\  /    A nd           | 
-     \\/     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
-
-Modified by
-Christian Lucas
-Institut für Thermodynamik
-Technische Universität Braunschweig 
-Germany
-
-\*---------------------------------------------------------------------------*/
-
-#include "mixedInternalEnergyRealFluids.H"
-#include "addToRunTimeSelectionTable.H"
-#include "fvPatchFieldMapper.H"
-#include "volFields.H"
-#include "basicPsiThermo.H"
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-namespace Foam
-{
-
-// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
-
-mixedInternalEnergyRealFluids::mixedInternalEnergyRealFluids
-(
-    const fvPatch& p,
-    const DimensionedField<scalar, volMesh>& iF
-)
-:
-    mixedFvPatchScalarField(p, iF)
-{
-    valueFraction() = 0.0;
-    refValue() = 0.0;
-    refGrad() = 0.0;
-}
-
-
-mixedInternalEnergyRealFluids::mixedInternalEnergyRealFluids
-(
-    const mixedInternalEnergyRealFluids& ptf,
-    const fvPatch& p,
-    const DimensionedField<scalar, volMesh>& iF,
-    const fvPatchFieldMapper& mapper
-)
-:
-    mixedFvPatchScalarField(ptf, p, iF, mapper)
-{}
-
-
-mixedInternalEnergyRealFluids::mixedInternalEnergyRealFluids
-(
-    const fvPatch& p,
-    const DimensionedField<scalar, volMesh>& iF,
-    const dictionary& dict
-)
-:
-    mixedFvPatchScalarField(p, iF, dict)
-{}
-
-
-mixedInternalEnergyRealFluids::mixedInternalEnergyRealFluids
-(
-    const mixedInternalEnergyRealFluids& tppsf
-)
-:
-    mixedFvPatchScalarField(tppsf)
-{}
-
-
-mixedInternalEnergyRealFluids::mixedInternalEnergyRealFluids
-(
-    const mixedInternalEnergyRealFluids& tppsf,
-    const DimensionedField<scalar, volMesh>& iF
-)
-:
-    mixedFvPatchScalarField(tppsf, iF)
-{}
-
-
-// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
-
-void mixedInternalEnergyRealFluids::updateCoeffs()
-{
-    if (updated())
-    {
-        return;
-    }
-
-    const basicPsiThermo& thermo = db().lookupObject<basicPsiThermo>
-    (
-        "thermophysicalProperties"
-    );
-
-    
-    const label patchi = patch().index();
-
-    mixedFvPatchScalarField& Tw = refCast<mixedFvPatchScalarField>
-    (
-        const_cast<fvPatchScalarField&>(thermo.T().boundaryField()[patchi])
-    );
-
-    Tw.evaluate();
-
-    fvPatchScalarField& pw =
-        const_cast<fvPatchScalarField&>(thermo.p().boundaryField()[patchi]);
-    pw.evaluate();
-
-    valueFraction() = Tw.valueFraction();
-    refValue() = thermo.e(pw,Tw.refValue(), patchi);
-    refGrad() = thermo.Cv(pw,Tw, patchi)*Tw.refGrad()
-      + patch().deltaCoeffs()*
-        (
-            thermo.e(pw,Tw, patchi)
-          - thermo.e(pw,Tw, patch().faceCells())
-        );
-
-    mixedFvPatchScalarField::updateCoeffs();
-}
-
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-makePatchTypeField(fvPatchScalarField, mixedInternalEnergyRealFluids);
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-} // End namespace Foam
-
-// ************************************************************************* //

src/thermophysicalModels/basic/derivedFvPatchFieldsRealGas/mixedInternalEnergyRealFluids/mixedInternalEnergyRealFluids.H

@@ -1,146 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     |
-    \\  /    A nd           | 
-     \\/     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::mixedInternalEnergyRealFluids
-
-Description
-    Mixed boundary conditions for internal energy
-
-SourceFiles
-    mixedInternalEnergyRealFluids.C
-
-Modified by
-Christian Lucas
-Institut für Thermodynamik
-Technische Universität Braunschweig 
-Germany
-
-\*---------------------------------------------------------------------------*/
-
-#ifndef mixedInternalEnergyRealFluids_H
-#define mixedInternalEnergyRealFluids_H
-
-#include "mixedFvPatchFields.H"
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-namespace Foam
-{
-
-/*---------------------------------------------------------------------------*\
-               Class mixedInternalEnergyRealFluids Declaration
-\*---------------------------------------------------------------------------*/
-
-class mixedInternalEnergyRealFluids
-:
-    public mixedFvPatchScalarField
-{
-
-public:
-
-    //- Runtime type information
-    TypeName("mixedInternalEnergy");
-
-
-    // Constructors
-
-        //- Construct from patch and internal field
-        mixedInternalEnergyRealFluids
-        (
-            const fvPatch&,
-            const DimensionedField<scalar, volMesh>&
-        );
-
-        //- Construct from patch, internal field and dictionary
-        mixedInternalEnergyRealFluids
-        (
-            const fvPatch&,
-            const DimensionedField<scalar, volMesh>&,
-            const dictionary&
-        );
-
-        //- Construct by mapping given mixedInternalEnergyRealFluids
-        // onto a new patch
-        mixedInternalEnergyRealFluids
-        (
-            const mixedInternalEnergyRealFluids&,
-            const fvPatch&,
-            const DimensionedField<scalar, volMesh>&,
-            const fvPatchFieldMapper&
-        );
-
-        //- Construct as copy
-        mixedInternalEnergyRealFluids
-        (
-            const mixedInternalEnergyRealFluids&
-        );
-
-        //- Construct and return a clone
-        virtual tmp<fvPatchScalarField> clone() const
-        {
-            return tmp<fvPatchScalarField>
-            (
-                new mixedInternalEnergyRealFluids(*this)
-            );
-        }
-
-        //- Construct as copy setting internal field reference
-        mixedInternalEnergyRealFluids
-        (
-            const mixedInternalEnergyRealFluids&,
-            const DimensionedField<scalar, volMesh>&
-        );
-
-        //- Construct and return a clone setting internal field reference
-        virtual tmp<fvPatchScalarField> clone
-        (
-            const DimensionedField<scalar, volMesh>& iF
-        ) const
-        {
-            return tmp<fvPatchScalarField>
-            (
-                new mixedInternalEnergyRealFluids(*this, iF)
-            );
-        }
-
-
-    // Member functions
-
-        // Evaluation functions
-
-            //- Update the coefficients associated with the patch field
-            virtual void updateCoeffs();
-};
-
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-} // End namespace Foam
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-#endif
-
-// ************************************************************************* //

src/thermophysicalModels/basic/derivedFvPatchFieldsRealGas/realFluidWallHeatTransfer/realFluidWallHeatTransferFvPatchScalarField.C

@@ -1,198 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     |
-    \\  /    A nd           | 
-     \\/     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
-
-Modified by
-Christian Lucas
-Institut für Thermodynamik
-Technische Universität Braunschweig 
-Germany
-
-\*---------------------------------------------------------------------------*/
-
-#include "realFluidWallHeatTransferFvPatchScalarField.H"
-#include "addToRunTimeSelectionTable.H"
-#include "fvPatchFieldMapper.H"
-#include "volFields.H"
-#include "basicPsiThermo.H"
-
-// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
-
-Foam::realFluidWallHeatTransferFvPatchScalarField::realFluidWallHeatTransferFvPatchScalarField
-(
-    const fvPatch& p,
-    const DimensionedField<scalar, volMesh>& iF
-)
-:
-    mixedFvPatchScalarField(p, iF),
-    Tinf_(p.size(), 0.0),
-    alphaWall_(p.size(), 0.0)
-{
-    refValue() = 0.0;
-    refGrad() = 0.0;
-    valueFraction() = 0.0;
-}
-
-
-Foam::realFluidWallHeatTransferFvPatchScalarField::realFluidWallHeatTransferFvPatchScalarField
-(
-    const realFluidWallHeatTransferFvPatchScalarField& ptf,
-    const fvPatch& p,
-    const DimensionedField<scalar, volMesh>& iF,
-    const fvPatchFieldMapper& mapper
-)
-:
-    mixedFvPatchScalarField(ptf, p, iF, mapper),
-    Tinf_(ptf.Tinf_, mapper),
-    alphaWall_(ptf.alphaWall_, mapper)
-{}
-
-
-Foam::realFluidWallHeatTransferFvPatchScalarField::realFluidWallHeatTransferFvPatchScalarField
-(
-    const fvPatch& p,
-    const DimensionedField<scalar, volMesh>& iF,
-    const dictionary& dict
-)
-:
-    mixedFvPatchScalarField(p, iF),
-    Tinf_("Tinf", dict, p.size()),
-    alphaWall_("alphaWall", dict, p.size())
-{
-    refValue() = Tinf_;
-    refGrad() = 0.0;
-    valueFraction() = 0.0;
-
-    if (dict.found("value"))
-    {
-        fvPatchField<scalar>::operator=
-        (
-            scalarField("value", dict, p.size())
-        );
-    }
-    else
-    {
-        evaluate();
-    }
-}
-
-
-Foam::realFluidWallHeatTransferFvPatchScalarField::realFluidWallHeatTransferFvPatchScalarField
-(
-    const realFluidWallHeatTransferFvPatchScalarField& tppsf
-)
-:
-    mixedFvPatchScalarField(tppsf),
-    Tinf_(tppsf.Tinf_),
-    alphaWall_(tppsf.alphaWall_)
-{}
-
-
-Foam::realFluidWallHeatTransferFvPatchScalarField::realFluidWallHeatTransferFvPatchScalarField
-(
-    const realFluidWallHeatTransferFvPatchScalarField& tppsf,
-    const DimensionedField<scalar, volMesh>& iF
-)
-:
-    mixedFvPatchScalarField(tppsf, iF),
-    Tinf_(tppsf.Tinf_),
-    alphaWall_(tppsf.alphaWall_)
-{}
-
-
-// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
-
-void Foam::realFluidWallHeatTransferFvPatchScalarField::autoMap
-(
-    const fvPatchFieldMapper& m
-)
-{
-    scalarField::autoMap(m);
-    Tinf_.autoMap(m);
-    alphaWall_.autoMap(m);
-}
-
-
-void Foam::realFluidWallHeatTransferFvPatchScalarField::rmap
-(
-    const fvPatchScalarField& ptf,
-    const labelList& addr
-)
-{
-    mixedFvPatchScalarField::rmap(ptf, addr);
-
-    const realFluidWallHeatTransferFvPatchScalarField& tiptf =
-        refCast<const realFluidWallHeatTransferFvPatchScalarField>(ptf);
-
-    Tinf_.rmap(tiptf.Tinf_, addr);
-    alphaWall_.rmap(tiptf.alphaWall_, addr);
-}
-
-
-void Foam::realFluidWallHeatTransferFvPatchScalarField::updateCoeffs()
-{
-    if (updated())
-    {
-        return;
-    }
-
-    const basicPsiThermo& thermo = db().lookupObject<basicPsiThermo>
-    (
-        "thermophysicalProperties"
-    );
-
-    const label patchi = patch().index();
-
-    const scalarField& Tw = thermo.T().boundaryField()[patchi];
-    const scalarField& pw = thermo.p().boundaryField()[patchi];
-    scalarField Cpw = thermo.Cp(pw,Tw, patchi);
-
-    valueFraction() =
-        1.0/
-        (
-            1.0
-          + Cpw*thermo.alpha().boundaryField()[patchi]
-           *patch().deltaCoeffs()/alphaWall_
-        );
-
-    mixedFvPatchScalarField::updateCoeffs();
-}
-
-
-void Foam::realFluidWallHeatTransferFvPatchScalarField::write(Ostream& os) const
-{
-    fvPatchScalarField::write(os);
-    Tinf_.writeEntry("Tinf", os);
-    alphaWall_.writeEntry("alphaWall", os);
-    writeEntry("value", os);
-}
-
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-namespace Foam
-{
-    makePatchTypeField(fvPatchScalarField, realFluidWallHeatTransferFvPatchScalarField);
-}
-
-// ************************************************************************* //

src/thermophysicalModels/basic/derivedFvPatchFieldsRealGas/realFluidWallHeatTransfer/realFluidWallHeatTransferFvPatchScalarField.H

@@ -1,201 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     |
-    \\  /    A nd           | 
-     \\/     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::realFluidWallHeatTransferFvPatchScalarField
-
-Description
-    Enthalpy boundary conditions for wall heat transfer
-
-SourceFiles
-    realFluidWallHeatTransferFvPatchScalarField.C
-
-Modified by
-Christian Lucas
-Institut für Thermodynamik
-Technische Universität Braunschweig 
-Germany
-
-\*---------------------------------------------------------------------------*/
-
-#ifndef realFluidWallHeatTransferFvPatchScalarField_H
-#define realFluidWallHeatTransferFvPatchScalarField_H
-
-#include "mixedFvPatchFields.H"
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-namespace Foam
-{
-
-/*---------------------------------------------------------------------------*\
-                  Class wallHeatTransferFvPatch Declaration
-\*---------------------------------------------------------------------------*/
-
-class realFluidWallHeatTransferFvPatchScalarField
-:
-    public mixedFvPatchScalarField
-{
-    // Private data
-
-        //- Tinf
-        scalarField Tinf_;
-
-        //- alphaWall
-        scalarField alphaWall_;
-
-
-public:
-
-    //- Runtime type information
-    TypeName("realFluidWallHeatTransfer");
-
-
-    // Constructors
-
-        //- Construct from patch and internal field
-        realFluidWallHeatTransferFvPatchScalarField
-        (
-            const fvPatch&,
-            const DimensionedField<scalar, volMesh>&
-        );
-
-        //- Construct from patch, internal field and dictionary
-        realFluidWallHeatTransferFvPatchScalarField
-        (
-            const fvPatch&,
-            const DimensionedField<scalar, volMesh>&,
-            const dictionary&
-        );
-
-        //- Construct by mapping given realFluidWallHeatTransferFvPatchScalarField
-        // onto a new patch
-        realFluidWallHeatTransferFvPatchScalarField
-        (
-            const realFluidWallHeatTransferFvPatchScalarField&,
-            const fvPatch&,
-            const DimensionedField<scalar, volMesh>&,
-            const fvPatchFieldMapper&
-        );
-
-        //- Construct as copy
-        realFluidWallHeatTransferFvPatchScalarField
-        (
-            const realFluidWallHeatTransferFvPatchScalarField&
-        );
-
-        //- Construct and return a clone
-        virtual tmp<fvPatchScalarField> clone() const
-        {
-            return tmp<fvPatchScalarField>
-            (
-                new realFluidWallHeatTransferFvPatchScalarField(*this)
-            );
-        }
-
-        //- Construct as copy setting internal field reference
-        realFluidWallHeatTransferFvPatchScalarField
-        (
-            const realFluidWallHeatTransferFvPatchScalarField&,
-            const DimensionedField<scalar, volMesh>&
-        );
-
-        //- Construct and return a clone setting internal field reference
-        virtual tmp<fvPatchScalarField> clone
-        (
-            const DimensionedField<scalar, volMesh>& iF
-        ) const
-        {
-            return tmp<fvPatchScalarField>
-            (
-                new realFluidWallHeatTransferFvPatchScalarField(*this, iF)
-            );
-        }
-
-
-    // Member functions
-
-        // Access
-
-            //- Return Tinf
-            const scalarField& Tinf() const
-            {
-                return Tinf_;
-            }
-
-            //- Return reference to Tinf to allow adjustment
-            scalarField& Tinf()
-            {
-                return Tinf_;
-            }
-
-            //- Return alphaWall
-            const scalarField& alphaWall() const
-            {
-                return alphaWall_;
-            }
-
-            //- Return reference to alphaWall to allow adjustment
-            scalarField& alphaWall()
-            {
-                return alphaWall_;
-            }
-
-
-        // Mapping functions
-
-            //- Map (and resize as needed) from self given a mapping object
-            virtual void autoMap
-            (
-                const fvPatchFieldMapper&
-            );
-
-            //- Reverse map the given fvPatchField onto this fvPatchField
-            virtual void rmap
-            (
-                const fvPatchScalarField&,
-                const labelList&
-            );
-
-
-        // Evaluation functions
-
-            //- Update the coefficients associated with the patch field
-            virtual void updateCoeffs();
-
-
-        //- Write
-        virtual void write(Ostream&) const;
-};
-
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-} // End namespace Foam
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-#endif
-
-// ************************************************************************* //

src/thermophysicalModels/basic/mixtures/basicMixture/basicMixtures.C

@@ -47,6 +47,13 @@ Description
 
 #include "addToRunTimeSelectionTable.H"
 
+#include "redlichKwong.H"
+#include "pengRobinson.H"
+#include "aungierRedlichKwong.H"
+#include "soaveRedlichKwong.H"
+#include "realGasSpecieThermo.H"
+#include "nasaHeatCapacityPolynomial.H"
+
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 namespace Foam
@@ -103,7 +110,17 @@ makeBasicMixturePhys
 makeBasicRealFluidMixture
 (
     pureMixture,
-    realGasSutherlandTransport,
+    sutherlandTransport,
+    realGasSpecieThermo,
+    nasaHeatCapacityPolynomial,
+    redlichKwong
+);
+
+
+makeBasicRealFluidMixture
+(
+    pureMixture,
+    sutherlandTransport,
     realGasSpecieThermo,
     nasaHeatCapacityPolynomial,
     pengRobinson
@@ -112,7 +129,7 @@ makeBasicRealFluidMixture
 makeBasicRealFluidMixture
 (
     pureMixture,
-    realGasSutherlandTransport,
+    sutherlandTransport,
     realGasSpecieThermo,
     nasaHeatCapacityPolynomial,
     aungierRedlichKwong
@@ -122,7 +139,7 @@ makeBasicRealFluidMixture
 makeBasicRealFluidMixture
 (
     pureMixture,
-    realGasSutherlandTransport,
+    sutherlandTransport,
     realGasSpecieThermo,
     nasaHeatCapacityPolynomial,
     soaveRedlichKwong
@@ -133,7 +150,7 @@ makeBasicRealFluidMixture
 makeBasicRealFluidMixture
 (
     pureMixture,
-    realGasConstTransport,
+    constTransport,
     realGasSpecieThermo,
     nasaHeatCapacityPolynomial,
     redlichKwong
@@ -143,7 +160,7 @@ makeBasicRealFluidMixture
 makeBasicRealFluidMixture
 (
     pureMixture,
-    realGasConstTransport,
+    constTransport,
     realGasSpecieThermo,
     nasaHeatCapacityPolynomial,
     pengRobinson
@@ -153,7 +170,7 @@ makeBasicRealFluidMixture
 makeBasicRealFluidMixture
 (
     pureMixture,
-    realGasConstTransport,
+    constTransport,
     realGasSpecieThermo,
     nasaHeatCapacityPolynomial,
     aungierRedlichKwong
@@ -163,7 +180,7 @@ makeBasicRealFluidMixture
 makeBasicRealFluidMixture
 (
     pureMixture,
-    realGasConstTransport,
+    constTransport,
     realGasSpecieThermo,
     nasaHeatCapacityPolynomial,
     soaveRedlichKwong

src/thermophysicalModels/basic/psiThermo/realGasHThermo/realGasHThermo.C

@@ -129,7 +129,7 @@ Foam::realGasHThermo<MixtureType>::realGasHThermo(const fvMesh& mesh)
         ),
         mesh,
         dimensionSet(0, 2, -2, 0, 0),
-        this->hRealBoundaryTypes()
+        this->hBoundaryTypes()
     ),
     
     rho_
@@ -184,7 +184,7 @@ Foam::realGasHThermo<MixtureType>::realGasHThermo(const fvMesh& mesh)
             h((this->rho_.boundaryField()[patchi])  ,  this->T_.boundaryField()[patchi], patchi);                    
     }
 
-    hRealBoundaryCorrection(h_);
+    hBoundaryCorrection(h_);
     calculate();
 
     // Switch on saving old time

src/thermophysicalModels/basic/psiThermo/realGasHThermo/realGasHThermos.C

@@ -42,10 +42,10 @@ Germany
 #include "soaveRedlichKwong.H"
 #include "nasaHeatCapacityPolynomial.H"
 #include "realGasSpecieThermo.H"
+#include "constTransport.H"
+#include "sutherlandTransport.H"
 
 #include "pureMixture.H"
-#include "realGasSutherlandTransport.H"
-#include "realGasConstTransport.H"
 #include "realGasHThermo.H"
 
 
@@ -60,7 +60,47 @@ makeBasicRealGasThermo
 (
     realGasHThermo,
     pureMixture,
-    realGasSutherlandTransport,
+    sutherlandTransport,
+    realGasSpecieThermo,
+    nasaHeatCapacityPolynomial,
+    pengRobinson
+);
+
+makeBasicRealGasThermo
+(
+    realGasHThermo,
+    pureMixture,
+    sutherlandTransport,
+    realGasSpecieThermo,
+    nasaHeatCapacityPolynomial,
+    aungierRedlichKwong
+);
+
+makeBasicRealGasThermo
+(
+    realGasHThermo,
+    pureMixture,
+    sutherlandTransport,
+    realGasSpecieThermo,
+    nasaHeatCapacityPolynomial,
+    redlichKwong
+);
+
+makeBasicRealGasThermo
+(
+    realGasHThermo,
+    pureMixture,
+    sutherlandTransport,
+    realGasSpecieThermo,
+    nasaHeatCapacityPolynomial,
+    soaveRedlichKwong
+);
+
+makeBasicRealGasThermo
+(
+    realGasHThermo,
+    pureMixture,
+    constTransport,
     realGasSpecieThermo,
     nasaHeatCapacityPolynomial,
     pengRobinson
@@ -73,89 +113,32 @@ makeBasicRealGasThermo
 (
     realGasHThermo,
     pureMixture,
-    realGasSutherlandTransport,
+    constTransport,
     realGasSpecieThermo,
     nasaHeatCapacityPolynomial,
     aungierRedlichKwong
 );
 
-
-
-
 makeBasicRealGasThermo
 (
     realGasHThermo,
     pureMixture,
-    realGasSutherlandTransport,
+    constTransport,
     realGasSpecieThermo,
     nasaHeatCapacityPolynomial,
     redlichKwong
 );
 
-
-
-
 makeBasicRealGasThermo
 (
     realGasHThermo,
     pureMixture,
-    realGasSutherlandTransport,
+    constTransport,
     realGasSpecieThermo,
     nasaHeatCapacityPolynomial,
     soaveRedlichKwong
 );
 
-makeBasicRealGasThermo
-(
-    realGasHThermo,
-    pureMixture,
-    realGasConstTransport,
-    realGasSpecieThermo,
-    nasaHeatCapacityPolynomial,
-    pengRobinson
-);
-
-
-
-
-makeBasicRealGasThermo
-(
-    realGasHThermo,
-    pureMixture,
-    realGasConstTransport,
-    realGasSpecieThermo,
-    nasaHeatCapacityPolynomial,
-    aungierRedlichKwong
-);
-
-
-
-
-makeBasicRealGasThermo
-(
-    realGasHThermo,
-    pureMixture,
-    realGasConstTransport,
-    realGasSpecieThermo,
-    nasaHeatCapacityPolynomial,
-    redlichKwong
-);
-
-
-
-
-makeBasicRealGasThermo
-(
-    realGasHThermo,
-    pureMixture,
-    realGasConstTransport,
-    realGasSpecieThermo,
-    nasaHeatCapacityPolynomial,
-    soaveRedlichKwong
-);
-
-
-
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 } // End namespace Foam

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

@@ -43,6 +43,19 @@ 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)
@@ -50,13 +63,16 @@ aungierRedlichKwong::aungierRedlichKwong(Istream& is)
     specie(is),
     pcrit_(readScalar(is)),
     Tcrit_(readScalar(is)),
-	azentricFactor_(readScalar(is)),
+    azentricFactor_(readScalar(is)),
-	rhocrit_(readScalar(is))
+    rhocrit_(readScalar(is)),    
+    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*this->W()/(Tstd*this->R())))
 {
     is.check("aungierRedlichKwong::aungierRedlichKwong(Istream& is)");
-    rhostd_=this->rho(Pstd,Tstd,Pstd*this->W()/(Tstd*this->R()));
-    rhoMax_=1500;	
-    rhoMin_=0.001;
 }
 
 
@@ -65,7 +81,8 @@ aungierRedlichKwong::aungierRedlichKwong(Istream& is)
 Ostream& operator<<(Ostream& os, const aungierRedlichKwong& pg)
 {
     os  << static_cast<const specie&>(pg)<< tab
-        << pg.pcrit_ << tab<< pg.Tcrit_<<tab<<pg.azentricFactor_<<tab<<pg.rhocrit_;
+        << pg.pcrit_ << tab<< pg.Tcrit_<< tab<<pg.azentricFactor_<< tab<<pg.rhocrit_;
+
     os.check("Ostream& operator<<(Ostream& os, const aungierRedlichKwong& st)");
     return os;
 }

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

@@ -45,8 +45,6 @@ Germany
 
 #include "specie.H"
 #include "autoPtr.H"
-#include "word.H"
-#include "scalar.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -62,6 +60,26 @@ class aungierRedlichKwong
     public specie
 
 {
+    // Private data
+        scalar pcrit_;
+        scalar Tcrit_; 
+        scalar azentricFactor_;
+        scalar rhocrit_;
+
+    //Aungier Redlich Kwong factors
+        scalar a_;
+        scalar b_;
+        scalar c_;
+        scalar n_;
+
+    //Density @STD, initialise after a, b!
+        scalar rhostd_;
+
+    // Static Private data
+
+    //should be read from the fvSolution file where rhoMax and rhoMin values must be define ( for rhoSimpleFoam)
+        static const scalar rhoMax_;
+        static const scalar rhoMin_; 
 
 
 public:
@@ -70,19 +88,20 @@ public:
     // Constructors
 
         //- Construct from components
-        inline aungierRedlichKwong(
+        inline aungierRedlichKwong
-		const specie& sp, 
+        (
-		scalar pcrit, 
+            const specie& sp, 
-		scalar Tcrit,
+            scalar pcrit, 
-		scalar azentricFactor,
+            scalar Tcrit,
-	        scalar rhocrit
+            scalar azentricFactor,
-		);
+            scalar rhocrit
+        );
 
         //- Construct from Istream
         aungierRedlichKwong(Istream&);
 
         //- Construct as named copy
-        inline aungierRedlichKwong(const word& name, aungierRedlichKwong&);
+        inline aungierRedlichKwong(const word& name, const aungierRedlichKwong&);
 
         //- Construct and return a clone
         inline autoPtr<aungierRedlichKwong> clone() const;
@@ -90,52 +109,72 @@ public:
         // Selector from Istream
         inline static autoPtr<aungierRedlichKwong> New(Istream& is);
 
-        //Member Variabels
+   // Member functions
-	    scalar pcrit_;
+
- 	    scalar Tcrit_; 
+	inline scalar rhostd() const;
-            scalar rhostd_;
-	    scalar azentricFactor_;
- 	    scalar rhocrit_;
-	    scalar rhoMax_;	//should be read from the fvSolution file where rhoMax and rhoMin values must be define ( for rhoSimpleFoam)
-    	    scalar rhoMin_;	
 
-    // Member functions
-        inline scalar pcrit() const;
-        inline scalar Tcrit() const;
-	inline scalar azentricFactor() const;
-	inline scalar rhostd()const;
-        inline scalar rhocrit() const;
         inline scalar p(const scalar rho, const scalar T) const;
 
-	//-Redlich Kwong factors
+        //first order derivatives 
-     	inline scalar a() const;
+        inline scalar dpdv(const scalar rho, const  scalar T) const;
-	inline scalar b() const;
+
-	inline scalar c() const;
-    inline scalar n() const;
-        //derivatives 
-        inline scalar dpdv(const scalar rho,const  scalar T) const;
         inline scalar dpdT(const scalar rho, const scalar T) const;
-        inline scalar dvdT(const scalar rho,const  scalar T) const;
+
+        inline scalar dvdT(const scalar rho, const  scalar T) const;
+
         inline scalar dvdp(const scalar rho, const scalar T) const; 
-        inline scalar isobarExpCoef(const scalar rho,const  scalar T) const; 
-        inline scalar isothermalCompressiblity(const scalar rho,const  scalar T) const; 
-        inline scalar integral_d2pdT2_dv(const scalar rho,const  scalar T) const ; // Used for cv
-        inline scalar d2pdv2(const scalar rho,const  scalar T) const;    // not Used At The Moment
-	inline scalar d2pdT2(const scalar rho,const  scalar T) const;    // not Used At The Moment
-        inline scalar d2pdvdT(const scalar rho,const  scalar T) const;    // not Used At The Moment
-	inline scalar d2vdT2(const scalar rho,const  scalar T) const;    // not Used At The Moment
-        inline scalar integral_p_dv(const scalar rho,const  scalar T) const;     //Used for internal Energy
-        inline scalar integral_dpdT_dv(const scalar rho,const  scalar T) const;  //Used for Entropy
 
-        //- Return density [kg/m^3] // rho0 is the starting point of the newton solver used to calculate rho
+        inline scalar isobarExpCoef(const scalar rho, const  scalar T) const; 
-        inline scalar rho(const scalar p,const  scalar T,const scalar rho0) const;
-        inline scalar rho(const scalar p,const  scalar T) const;
 
-        //- Return compressibility drho/dp [s^2/m^2]   
+        inline scalar isothermalCompressiblity
+        (
+            const scalar rho,
+            const scalar T
+        ) const; 
+        
+        // Used for cv
+        inline scalar integral_d2pdT2_dv
+        (
+            const scalar rho,
+            const scalar T
+        ) const ; 
+
+        // second order derivatives, not Used At The Moment
+        inline scalar d2pdv2(const scalar rho, const  scalar T) const;    
+
+	inline scalar d2pdT2(const scalar rho,   scalar T) const;    
+
+        inline scalar d2pdvdT(const scalar rho, const  scalar T) const;    
+
+	inline scalar d2vdT2(const scalar rho, const  scalar T) const;    
+
+        //Used for internal Energy
+        inline scalar integral_p_dv(const scalar rho, const  scalar T) const;
+
+        //Used for Entropy
+        inline scalar integral_dpdT_dv(const scalar rho, const  scalar T) const;  
+
+        //- Return density [kg/m^3] 
+        // rho0 is the starting point of the newton solver used to calculate rho
+        inline scalar rho
+	(
+            const scalar p,
+            const scalar T,
+            const scalar rho0
+        ) const;
+ 
+        inline scalar rho(const scalar p, const  scalar T) const;
+
+        //- Return compressibility drho/dp at T=constant [s^2/m^2]
         inline scalar psi(const scalar rho, const scalar T) const;
 
-        //- Return compression factor []	// rho0 is the starting point of the newton solver used to calculate rho
+        //- Return compression factor  []
-        inline scalar Z(const scalar p,const  scalar T,const scalar rho0) const;
+        inline scalar Z
+        (
+            const scalar p,
+            const  scalar T, 
+            const scalar rho0
+        ) const;
 
 
     // Member operators

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

@@ -40,230 +40,9 @@ namespace Foam
 
 
 
-
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
-
-// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
-
-inline scalar aungierRedlichKwong::pcrit()const
-{
-return pcrit_;
-}
-
-
-
-inline scalar aungierRedlichKwong::Tcrit()const
-{
-return Tcrit_;
-}
-
-
-inline scalar aungierRedlichKwong::rhostd()const
-{
-return rhostd_;
-}
-
-inline scalar aungierRedlichKwong::rhocrit()const
-{
-return rhocrit_;
-}
-
-// Returns the Azentric Factor (Acentric Factor)
-inline scalar aungierRedlichKwong::azentricFactor() const
-{
-return azentricFactor_;
-}
-
-//returns the pressure for a given density and temperature    
-inline scalar aungierRedlichKwong::p(const scalar rho,const scalar T) const
-{
-
-return this->RR*T/((this->W()/rho)-this->b()+this->c())
--this->a()*pow(T,-this->n())/(pow(this->Tcrit(),-this->n())*(this->W()/rho)*((this->W()/rho)+this->b()));
-}
-
-
-
-// Factor a of the redlich Kwong equation of state  
-//(molar values)
-inline scalar aungierRedlichKwong::a() const
-{
-return 0.42747*pow(this->RR,2)*pow(this->Tcrit(),2)/this->pcrit();
-}
-
-
-
-// Faktor b of the redlich Kwong equation of state 
-//(molar values)
-inline scalar aungierRedlichKwong::b() const
-{
-return 0.08664*this->RR*this->Tcrit()/this->pcrit();
-}
-// Factor c of the redlich Kwong equation of state 
-//(molar values)
-inline scalar aungierRedlichKwong::c() const
-{
-return this->RR*this->Tcrit()/(this->pcrit()+(this->a()/(this->W()/this->rhocrit()*(this->W()/this->rhocrit()+this->b()))))+this->b()-this->W()/this->rhocrit();
-}
-// Factor n of the redlich Kwong equation of state 
-inline scalar aungierRedlichKwong::n() const
-{
-return 0.4986+1.2735*this->azentricFactor()+0.4754*pow(this->azentricFactor(),2);
-}
-
-
-//* * * * * * * * * * * * * Derivatives * * * * * * * * * * * //
-
-
-
-//Real deviative dp/dv at constant temperature 
-//(molar values)
-inline scalar aungierRedlichKwong::dpdv(const scalar rho, const scalar T) const
-{
-
-
- return  this->a()*pow(T,-this->n())*pow(this->Tcrit(),this->n())*(this->b()+2*(this->W()/rho))/(pow((this->W()/rho),2)*pow((this->b()+(this->W()/rho)),2))-this->RR*T/(pow((this->b()-(this->W()/rho)-this->c()),2));
-}
-
-
-
-
-//Real deviative dp/dT at constant molar volume 
-//(molar values)
-inline scalar aungierRedlichKwong::dpdT(const scalar rho, const scalar T) const
-{
-return this->a()*this->n()*pow(T,-this->n()-1)*pow(this->Tcrit(),this->n())/((this->W()/rho)*((this->W()/rho)+this->b()))-
-this->RR/(this->b()-(this->W()/rho)-this->c());
-}
-
-
-
-
-//Real  deviative dv/dT at constant pressure 
-//using implicit differentiation
-//(molar values)
-inline scalar aungierRedlichKwong::dvdT(const scalar rho,const scalar T) const
-{
-return (-1)*this->dpdT(rho,T)/this->dpdv(rho,T);
-}
-
-
-
-
-//Real  deviative dv/dp at constant temperature 
-//(molar values)
-inline scalar aungierRedlichKwong::dvdp(const scalar rho,const scalar T)  const 
-{
-return 1/this->dpdv(rho,T);
-}
-
-
-
-
-
-
-//needed to calculate the internal energy
-//(molar values)
-inline scalar aungierRedlichKwong::integral_p_dv(const scalar rho,const scalar T) const
-{
-return -pow((T/this->Tcrit()),-this->n())*(this->a()*log((this->W()/rho))/(this->b())-this->a()*log(this->b()+(this->W()/rho))/this->b())+this->RR*T*log((this->W()/rho)-this->b()+this->c());
-}
-
-
-
-//needed to calculate the entropy
-//(molar values)
-inline scalar aungierRedlichKwong::integral_dpdT_dv(const scalar rho,const scalar T) const
-{
-return pow(T,-this->n()-1)*pow(this->Tcrit(),this->n())*
-(this->a()*this->n()*log((this->W()/rho))/(this->b()) 
--this->a()*this->n()*log(this->b()+(this->W()/rho))/(this->b()))
-+ this->RR*log(-this->b()+this->c()+(this->W()/rho));
-}
-
-
-
-
-//* * * * * * * * * * * * * second order Derivative based functions * * * * * * * * * * * //
-
-
-
-//(molar values)
-inline scalar aungierRedlichKwong::d2pdT2(const scalar rho,const scalar T) const
-{
-return -this->a()*this->n()*pow(T,-this->n()-2)*pow(this->Tcrit(),this->n())*(this->n()+1)/((this->W()/rho)*(this->b()+(this->W()/rho)));
-}
-
-
-
-
-//(molar values)
-inline scalar aungierRedlichKwong::d2pdv2(const scalar rho,const scalar T) const
-{
- return -2*this->a()*pow(T,-this->n())*pow(this->Tcrit(),this->n())*(pow(this->b(),2)+3*this->b()*(this->W()/rho)+3*pow((this->W()/rho),2))/(pow((this->W()/rho),3)*pow((this->b()+(this->W()/rho)),3))-2*this->RR*T/(pow((this->b()-(this->W()/rho)-this->c()),3));
-}
-
-
-
-
-//(molar values)
-//using second order implicit differentiation 
-inline scalar aungierRedlichKwong::d2vdT2(const scalar rho, const scalar T) const
-{
-return -(pow(this->dpdT(rho,T),2)*this->d2pdv2(rho,T)+   pow(this->dpdv(rho,T),2)  *this->d2pdT2(rho,T)- 2*this->dpdv(rho,T)*this->dpdT(rho,T)*this->d2pdvdT(rho,T))/(   pow(this->dpdv(rho,T),3));
-}
-
-
-
-
-//(molar values)
-inline scalar aungierRedlichKwong::d2pdvdT(const scalar rho, const scalar T) const
-{
-return -this->a()*this->n()*pow(T,-this->n()-1)*pow(this->Tcrit(),this->n())*(this->b()+2*(this->W()/rho))/(pow((this->W()/rho),2)*pow((this->b()+(this->W()/rho)),2))-this->RR/(pow((this->b()-this->c()-(this->W()/rho)),2));
-}
-
-
-
-
-// the result of this intergal is needed for the nasa based cp polynomial 
-//(molar values)
-inline scalar aungierRedlichKwong::integral_d2pdT2_dv(const scalar rho,const scalar T) const 
-{
-return pow(T,-this->n()-2)*pow(this->Tcrit(),this->n())*(this->a()*this->n()*(this->n()+1)*log(this->b()+(this->W()/rho))/this->b()-this->a()*this->n()*(1+this->n())*log((this->W()/rho))/this->b());
-}
-
-
-
-
-
-//* * * * * * * * * * * * * thermodynamic properties * * * * * * * * * * * //
-
-
-//Isobar expansion Coefficent beta = 1/v (dv/dt) at constant p
-//(molar values)
-inline scalar aungierRedlichKwong::isobarExpCoef(const scalar rho,const scalar T)  const
-{
- return this->dvdT(rho, T)*rho/this->W();
-}
-
-//isothemal compressiblity kappa 
-//(molar values)
-inline scalar aungierRedlichKwong::isothermalCompressiblity(const scalar rho,const scalar T) const
-{
-return this->isobarExpCoef(rho, T)/this->dpdT(rho, T);
-//also possible : return -this->dvdp(rho,T)*rho/this->W();
-}
-
-
-
-
-
-
-
-
-
 
 // Construct from components
-// Starting GUESS for the density by ideal gas law
 inline aungierRedlichKwong::aungierRedlichKwong
 (
     const specie& sp,
@@ -276,26 +55,32 @@ inline aungierRedlichKwong::aungierRedlichKwong
     specie(sp),    
     pcrit_(pcrit),
     Tcrit_(Tcrit),
-	azentricFactor_(azentricFactor),
+    azentricFactor_(azentricFactor),
-	rhocrit_(rhocrit)
+    rhocrit_(rhocrit),
-{
+    a_(0.42747*pow(this->RR,2)*pow(Tcrit_,2)/pcrit_),
-	rhostd_=this->rho(Pstd,Tstd,Pstd*this->W()/(Tstd*this->R()));	// Starting GUESS for the density by ideal gas law
+    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*this->W()/(Tstd*this->R())))
+{}
 
 
-// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
-
 // Construct as named copy
-
+inline aungierRedlichKwong::aungierRedlichKwong(const word& name, const aungierRedlichKwong& pg)
-inline aungierRedlichKwong::aungierRedlichKwong(const word& name,  aungierRedlichKwong& pg)
 :
     specie(name, pg),
     pcrit_(pg.pcrit_),
     Tcrit_(pg.Tcrit_),
-	azentricFactor_(pg.azentricFactor_)
+    azentricFactor_(pg.azentricFactor_),
-{
+    rhocrit_(pg.rhocrit_),    
-  pg.rhostd_=this->rho(Pstd,Tstd,  (Pstd*this->W()/(Tstd*this->R())));	// Starting GUESS for the density by ideal gas law
+    a_(pg.a_),
-}
+    b_(pg.b_),
+    c_(pg.c_),
+    n_(pg.n_),
+    rhostd_(pg.rhostd_)
+{}
+
 
 
 // Construct and return a clone
@@ -312,105 +97,267 @@ inline autoPtr<aungierRedlichKwong> aungierRedlichKwong::New(Istream& is)
 }
 
 
-// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+// * * * * * * * * * * * * *  Member Functions  * * * * * * * * * * * //
+inline scalar aungierRedlichKwong::rhostd()const
+{
+    return rhostd_;
+}
+
+//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_)
+        -a_*pow(T,-n_)/(pow(Tcrit_,-n_)*Vm*(Vm+b_));
+}
+
+
+//Real deviative dp/dv at constant temperature 
+//(molar values)
+inline scalar aungierRedlichKwong::dpdv(const scalar rho, const scalar T) const
+{
+    scalar Vm = this->W()/rho;
+    return  a_*pow(T,-n_)*pow(Tcrit_,n_)*(b_+2*Vm)/(pow(Vm,2)*pow((b_+Vm),2))
+        -this->RR*T/(pow((b_-Vm-c_),2));
+}
+
+
+//Real deviative dp/dT at constant molar volume 
+//(molar values)
+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_))
+        -this->RR/(b_-Vm-c_);
+}
+
+
+//Real  deviative dv/dT at constant pressure 
+//using implicit differentiation
+//(molar values)
+inline scalar aungierRedlichKwong::dvdT(const scalar rho,const scalar T) const
+{
+    return (-1)*this->dpdT(rho,T)/this->dpdv(rho,T);
+}
+
+
+//Real  deviative dv/dp at constant temperature 
+//(molar values)
+inline scalar aungierRedlichKwong::dvdp(const scalar rho,const scalar T)  const 
+{
+    return 1/this->dpdv(rho,T);
+}
+
+
+//needed to calculate the internal energy
+//(molar values)
+inline scalar aungierRedlichKwong::integral_p_dv
+(
+    const scalar rho,
+    const scalar T
+) const
+{
+    scalar Vm = this->W()/rho;
+    return -pow((T/Tcrit_),-n_)*(a_*log(Vm)/(b_)
+        -a_*log(b_+Vm)/b_)+this->RR*T*log(Vm-b_+c_);
+}
+
+
+
+//needed to calculate the entropy
+//(molar values)
+inline scalar aungierRedlichKwong::integral_dpdT_dv
+( 
+    const scalar rho,
+    const scalar T
+) const
+{
+    scalar Vm = this->W()/rho;
+    return pow(T,-n_-1)*pow(Tcrit_,n_)*(a_*n_*log(Vm)/(b_) 
+        -a_*n_*log(b_+Vm)/(b_))+ this->RR*log(-b_+c_+Vm);
+}
+
+
+//(molar values)
+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));
+}
+
+
+
+
+//(molar values)
+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)
+        +3*b_*Vm+3*pow(Vm,2))/(pow(Vm,3)*pow((b_+Vm),3))
+        -2*this->RR*T/(pow((b_-Vm-c_),3));
+}
+
+
+
+
+//(molar values)
+//using second order implicit differentiation 
+inline scalar aungierRedlichKwong::d2vdT2(const scalar rho, const scalar T) const
+{
+    return 
+    -(
+          pow(this->dpdT(rho,T),2)*this->d2pdv2(rho,T)
+        + pow(this->dpdv(rho,T),2)*this->d2pdT2(rho,T)
+        - 2*this->dpdv(rho,T)*this->dpdT(rho,T)*this->d2pdvdT(rho,T)
+    )
+    /(pow(this->dpdv(rho,T),3));
+}
+
+
+
+
+//(molar values)
+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))
+        -this->RR/(pow((b_-c_-Vm),2));
+}
+
+
+
+
+// the result of this intergal is needed for the nasa based cp polynomial 
+//(molar values)
+inline scalar aungierRedlichKwong::integral_d2pdT2_dv
+(
+    const scalar rho,
+    const scalar T
+) const 
+{
+    scalar Vm = this->W()/rho;
+    return pow(T,-n_-2)*pow(Tcrit_,n_)*(a_*n_*(n_+1)*log(b_+Vm)
+        /b_-a_*n_*(1+n_)*log(Vm)/b_);
+}
+
+
+//Isobar expansion Coefficent beta = 1/v (dv/dt) at constant p
+//(molar values)
+inline scalar aungierRedlichKwong::isobarExpCoef(const scalar rho,const scalar T)  const
+{
+    return this->dvdT(rho, T)*rho/this->W();
+}
+
+
+//isothemal compressiblity kappa (not Thermal conductivity)
+//(molar values)
+inline scalar aungierRedlichKwong::isothermalCompressiblity(const scalar rho,const scalar T) const
+{
+    return this->isobarExpCoef(rho, T)/this->dpdT(rho, T);
+    //also possible : return -this->dvdp(rho,T)*rho/this->W();
+}
+
 
 //- Return density [kg/m^3]
-inline scalar aungierRedlichKwong::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
 {
  
-
    scalar molarVolumePrevIteration;
    scalar molarVolume;
    int iter=0;
    int maxIter_=400;
    scalar tol_=1e-8;
-   int i;
+   scalar rho1=rhoMax_;
-   scalar rho1=rhoMax_, rho2=rhoMin_,rho3, f1,f2,f3;
+   scalar rho2=rhoMin_;
-
-
-molarVolume=this->W()/rho0;	 
 
+   molarVolume=this->W()/rho0;	 
 
    do
    {
-   molarVolumePrevIteration= molarVolume;
+       molarVolumePrevIteration= molarVolume;
    
-   i=0;
+       label i=0;
-      	do
+       do
-   	{		
+       {
-  	molarVolume=molarVolumePrevIteration-((this->p((this->W()/molarVolumePrevIteration),T)-p)/(this->dpdv((this->W()/       	
+           molarVolume=molarVolumePrevIteration
-  	molarVolumePrevIteration),T)))/(pow(2,i)); 
+               -(
-	i++;
+                   (this->p((this->W()/molarVolumePrevIteration),T) - p)
-		if(i>8) 
+                   /(this->dpdv((this->W()/molarVolumePrevIteration),T))
-		{
+               )/pow(2,i);
 
-		//using bisection methode as backup, solution must be between rho=0.001 to rho=1500; 
+           i++;
-		  for(i=0;i<200;i++) 
+           if (i>8) 
-		  {
+           {
-                  	f1= (this->p(rho1,T)-p);
+               //using bisection methode as backup,
-                  	f2= (this->p(rho2,T)-p);
+               //solution must be between rho=0.001 to rho=1500;
-                  	rho3=(rho1+rho2)/2;
+               for(i=0; i<200; i++)
-		  	f3=(this->p(rho3,T)-p);
+               {
+                   scalar f1 = this->p(rho1,T) - p;
+                   scalar f2 = this->p(rho2,T) - p;
+                   scalar rho3 = (rho1 + rho2)/2;
+                   scalar f3 = this->p(rho3,T) - p;
 	
-			if ((f2<0 && f3>0)||(f2>0 &&f3<0))	
+                   if ((f2 < 0 && f3 > 0) || (f2 > 0 && f3 < 0))
-		  	{
+                   {
-		 		rho1=rho3;
+                       rho1=rho3;
-			}  
+                   }  
- 		  	else if ((f1<0 && f3>0)||(f1>0 &&f3<0))
+                   else if ((f1 < 0 && f3 > 0)||(f1 > 0 && f3 < 0))
-                  	{
+                   {
-		 		rho2=rho3;
+                       rho2=rho3;
-			}
+                   }		 
+                   else	
+                   {
+                       rho2=(rho2 + rho3)/2;
+                   }
 
-			else	
+                   if(mag(f3) < p*tol_)
-			{
+                   {	  	  
-                        	rho2=(rho2+rho3)/2;
+                       molarVolume=this->W()/rho3;
-			}
+                       molarVolumePrevIteration=this->W()/rho3;
-
+                       break; 
-			if(mag(f3)<p*tol_)
+                   }
-                  	{	  	  
+                   else
-				molarVolume=this->W()/rho3;
+                   {
-                 	 	molarVolumePrevIteration=this->W()/rho3;
+                       molarVolumePrevIteration=this->W()/rho3;
-				break; 
+                   }
-		  	}
+               }		  
-			else
+           }
-			{
+       }
-				molarVolumePrevIteration=this->W()/rho3;
+       while
-			}
+       (
-		    }		  
+           mag(this->p((this->W()/molarVolume),T) - p)
-  		}  
+         > mag(this->p((this->W()/molarVolumePrevIteration),T) - p)
-        }
+       );
-        while(mag(this->p((this->W()/molarVolume),T)-p)>mag(this->p((this->W()/molarVolumePrevIteration),T)-p));
-
-           if (iter++ > maxIter_)
-        {
-            FatalErrorIn
-            (
-                "inline scalar aungierRedlichKwong::rho(const scalar p,const scalar T,const scalar rho0) const "
-            )   << "Maximum number of iterations exceeded"
-                << abort(FatalError);
-        }
 
+       if (iter++ > maxIter_)
+       {
+           FatalErrorIn
+           (
+               "inline scalar redlichKwong::rho(const scalar p, const scalar T, const scalar rho0) const "
+           )   << "Maximum number of iterations exceeded"
+               << abort(FatalError);
+       }
    }
-    while(mag(molarVolumePrevIteration-molarVolume)>tol_*(this->W()/rho0));
+   while(mag(molarVolumePrevIteration-molarVolume) > tol_*(this->W()/rho0));
 
    return this->W()/molarVolume;
- 
 }
 
-//- Return density [kg/m^3]on 
+//- Return density [kg/m^3]
 inline scalar aungierRedlichKwong::rho(const scalar p,const scalar T) const
 {
-
+    //using perfect gas equation as starting point
-
+    return rho(p,T,p/(this->R()*T));
-scalar rho0=p/(this->R()*T); //using perfect gas equation as starting point
-return rho(p,T,rho0);
-
 }
 
-//- Return compressibility drho/dp [s^2/m^2]
+//- Return compressibility drho/dp at T=constant [s^2/m^2]
 inline scalar aungierRedlichKwong::psi(const scalar rho, const scalar T) const
 {
-   return -this->dvdp(rho,T)*pow(rho,2)/this->W();
+    return -this->dvdp(rho,T)*pow(rho,2)/this->W();
 }
 
 //- Return compression factor []

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

@@ -41,6 +41,18 @@ 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)
@@ -48,14 +60,13 @@ pengRobinson::pengRobinson(Istream& is)
     specie(is),
     pcrit_(readScalar(is)),
     Tcrit_(readScalar(is)),
-    azentricFactor_(readScalar(is))
+    azentricFactor_(readScalar(is)),
+    a_(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)), 
+    rhostd_(this->rho(Pstd,Tstd,Pstd*this->W()/(Tstd*this->R())))  	
 {
     is.check("pengRobinson::pengRobinson(Istream& is)");
-    rhostd_=this->rho(Pstd,Tstd,Pstd*this->W()/(Tstd*this->R()));
-    rhoMax_=1500;	
-    rhoMin_=0.001;
-
-
    
 }
 

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

@@ -47,8 +47,7 @@ Germany
 
 #include "specie.H"
 #include "autoPtr.H"
-#include "word.H"
+
-#include "scalar.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -64,6 +63,25 @@ class pengRobinson
     public specie
 
 {
+   //Member Variabels
+    scalar pcrit_;
+    scalar Tcrit_; 
+    scalar azentricFactor_;
+
+    //-Peng Robinson factors
+    scalar a_;
+    scalar b_;
+    scalar n_;
+
+    //- Density @STD, initialise after a, b!
+    scalar rhostd_;
+
+    // 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_;
 
 
 public:
@@ -83,7 +101,7 @@ public:
         pengRobinson(Istream&);
 
         //- Construct as named copy
-        inline pengRobinson(const word& name, pengRobinson&);
+        inline pengRobinson(const word& name,const pengRobinson&);
 
         //- Construct and return a clone
         inline autoPtr<pengRobinson> clone() const;
@@ -91,49 +109,78 @@ public:
         // Selector from Istream
         inline static autoPtr<pengRobinson> New(Istream& is);
 
-        //Member Variabels
+
-	    scalar pcrit_;
+	
- 	    scalar Tcrit_; 
-            scalar rhostd_;
-	    scalar azentricFactor_;
-	    scalar rhoMax_;	//should be read from the fvSolution file where rhoMax and rhoMin values must be define ( for rhoSimpleFoam)
-    	    scalar rhoMin_;	
 
     // Member functions
-        inline scalar pcrit() const;
-        inline scalar Tcrit() const;
-	inline scalar azentricFactor() const;
 	inline scalar rhostd()const;
+
         inline scalar p(const scalar rho, const scalar T) const;
 
-	//-Redlich Kwong factors
+        //first order derivatives 
-     	inline scalar a() const;
-	inline scalar b() const;
-	inline scalar n() const;
-
-        //derivatives 
         inline scalar dpdv(const scalar rho,const  scalar T) const;
-        inline scalar dpdT(const scalar rho, const scalar T) const;
-        inline scalar dvdT(const scalar rho,const  scalar T) const;
-        inline scalar dvdp(const scalar rho, const scalar T) const; 
-        inline scalar isobarExpCoef(const scalar rho,const  scalar T) const; 
-        inline scalar isothermalCompressiblity(const scalar rho,const  scalar T) const; 
-        inline scalar integral_d2pdT2_dv(const scalar rho,const  scalar T) const ; 	// Used for cv
-        inline scalar d2pdv2(const scalar rho,const  scalar T) const;    		// not Used At The Moment
-	inline scalar d2pdT2(const scalar rho,const  scalar T) const;    		// not Used At The Moment
-        inline scalar d2pdvdT(const scalar rho,const  scalar T) const;    		// not Used At The Moment
-	inline scalar d2vdT2(const scalar rho,const  scalar T) const;    		// not Used At The Moment
-        inline scalar integral_p_dv(const scalar rho,const  scalar T) const;     	//Used for internal Energy
-        inline scalar integral_dpdT_dv(const scalar rho,const  scalar T) const;  	//Used for Entropy
 
-        //- Return density [kg/m^3] // rho0 is the starting point of the newton solver used to calculate rho
+        inline scalar dpdT(const scalar rho, const scalar T) const;
-        inline scalar rho(const scalar p,const  scalar T,const scalar rho0) const;
+
+        inline scalar dvdT(const scalar rho,const  scalar T) const;
+
+        inline scalar dvdp(const scalar rho, const scalar T) const; 
+
+        inline scalar isobarExpCoef
+        (
+            const scalar rho,
+            const  scalar T
+        ) const; 
+
+        inline scalar isothermalCompressiblity
+        (
+            const scalar rho,
+            const  scalar T
+        ) const; 
+
+        // Used for cv
+        inline scalar integral_d2pdT2_dv
+        ( 
+            const scalar rho,
+            const  scalar T
+        ) const ; 
+	
+        // second order derivatives, not Used At The Moment
+        inline scalar d2pdv2(const scalar rho,const  scalar T) const;    		
+	
+        inline scalar d2pdT2(const scalar rho,const  scalar T) const;    		
+        
+        inline scalar d2pdvdT(const scalar rho,const  scalar T) const;    		
+	
+        inline scalar d2vdT2(const scalar rho,const  scalar T) const;    
+	
+        //Used for internal Energy
+        inline scalar integral_p_dv(const scalar rho,const  scalar T) const;
+        
+        //Used for Entropy
+        inline scalar integral_dpdT_dv(const scalar rho,const  scalar T) const;  	
+
+        //- Return density [kg/m^3] 
+        // rho0 is the starting point of the newton solver used to calculate rho
+        inline scalar rho
+        (
+            const scalar p,
+            const scalar T,
+            const scalar rho0
+        ) const;
+
         inline scalar rho(const scalar p,const  scalar T) const;
-        //- Return compressibility rho/p [s^2/m^2]   // rho0 is the starting point of the newton solver used to calculate rho
+
+        //- Return compressibility drho/dp at T=constant [s^2/m^2]
         inline scalar psi(const scalar rho, const scalar T) const;
 
         //- Return compression factor []
-        inline scalar Z(const scalar p,const  scalar T,const scalar rho0) const;
+        inline scalar Z
+        (
+            const scalar p,
+            const scalar T,
+            const scalar rho0
+        ) const;
 
 
     // Member operators

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

@@ -39,243 +39,7 @@ namespace Foam
 {
 
 
-
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
-
-
-// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
-
-inline scalar pengRobinson::pcrit()const
-{
-return pcrit_;
-}
-
-
-
-inline scalar pengRobinson::Tcrit()const
-{
-return Tcrit_;
-}
-
-
-inline scalar pengRobinson::rhostd()const
-{
-return rhostd_;
-}
-
-
-// Returns the Azentric Factor (Acentric Factor)
-inline scalar pengRobinson::azentricFactor() const
-{
-return azentricFactor_;
-}
-
-//returns the pressure for a given density and temperature    
-inline scalar pengRobinson::p(const scalar rho,const scalar T) const
-{
-return  this->RR*T/((this->W()/rho)-this->b())-(this->a()*pow((1+this->n()*(1-pow((T/this->Tcrit()),0.5))),2))/(pow((this->W()/rho),2)+2*this->b()*(this->W()/rho)-pow(this->b(),2));
-}
-
-
-
-// Factor a of the pengRobinson equation of state  
-//(molar values)
-inline scalar pengRobinson::a() const
-{
-return 0.457235*pow(this->RR,2)*pow(this->Tcrit(),2)/this->pcrit();
-}
-
-
-
-// Factor b of the pengRobinson equation of state 
-//(molar values)
-inline scalar pengRobinson::b() const
-{
-return 0.077796*this->RR*this->Tcrit()/this->pcrit(); 
-}
-
-// Factor nb of the pengRobinson equation of state 
-//(molar values)
-inline scalar pengRobinson::n() const
-{
-// this equation is only valid for asentricFactors <0.49. 
-return 0.37464+1.54226*this->azentricFactor()-0.26992*pow(this->azentricFactor(),2); 
-}
-
-
-//* * * * * * * * * * * * * Derivatives * * * * * * * * * * * //
-
-
-
-//Real deviative dp/dv at constant temperature 
-//(molar values)
-inline scalar pengRobinson::dpdv(const scalar rho, const scalar T) const
-{
-return -(4*this->a()*this->n()*this->Tcrit()*(this->b()-(this->W()/rho))*(pow(this->b(),2)-pow((this->W()/rho),2))*(this->n()+1)
-*pow((T/this->Tcrit()),0.5)
-+this->Tcrit()*(-2*this->a()*pow((this->n()+1),2)*(pow(this->b(),3)-pow(this->b(),2)*(this->W()/rho)
--this->b()*pow((this->W()/rho),2)+pow((this->W()/rho),3))
-+this->RR*T*(pow(this->b(),4)-4*pow(this->b(),3)*(this->W()/rho)+2*pow(this->b(),2)*pow((this->W()/rho),2)
-+4*this->b()*pow((this->W()/rho),3)+pow((this->W()/rho),4)))
--2*this->a()*pow(this->n(),2)*T*(pow(this->b(),3)-pow(this->b(),2)*(this->W()/rho)-this->b()*pow((this->W()/rho),2)+pow((this->W()/rho),3)))
-/(this->Tcrit()*pow((this->b()-(this->W()/rho)),2)*pow((pow(this->b(),2)-2*this->b()*(this->W()/rho)-pow((this->W()/rho),2)),2));
-}
-
-
-
-//Real deviative dp/dT at constant molar volume 
-//(molar values)
-inline scalar pengRobinson::dpdT(const scalar rho, const scalar T) const
-{
-return (-this->a()*this->n()*(this->n()+1)*pow((T/this->Tcrit()),0.5)/(T*(pow(this->b(),2)-2*this->b()*(this->W()/rho)-pow((this->W()/rho),2)))+this->a()*pow(this->n(),2)/(this->Tcrit()*(pow(this->b(),2)-2*this->b()*(this->W()/rho)-pow((this->W()/rho),2)))-this->RR/(this->b()-(this->W()/rho)));
-}
-
-
-
-
-//Real  deviative dv/dT at constant pressure 
-//by using implicit differentiation
-//(molar values)
-inline scalar pengRobinson::dvdT(const scalar rho,const scalar T) const
-{
-return (-1)*this->dpdT(rho,T)/this->dpdv(rho,T);
-}
-
-
-
-
-//(molar values)
-inline scalar pengRobinson::dvdp(const scalar rho,const scalar T)  const 
-{
-return 1/this->dpdv(rho,T);
-}
-
-
-
-
-
-//(molar values)
-//needed to calculate the internal energy
-inline scalar pengRobinson::integral_p_dv(const scalar rho,const scalar T) const
-{
-return  pow(2,0.5)*this->a()*(2*this->n()*this->Tcrit()*(this->n()+1)*pow(T/this->Tcrit(),0.5)
--this->Tcrit()*(pow(this->n(),2)+2*this->n()+1)-pow(this->n(),2)*T)*log(this->b()*(1-pow(2,0.5))+(this->W()/rho))/(4*this->b()*this->Tcrit())
-+this->RR*T*log((this->W()/rho)-this->b())
--pow(2,0.5)*this->a()*(2*this->n()*this->Tcrit()*(this->n()+1)*pow(T/this->Tcrit(),0.5)
--this->Tcrit()*(pow(this->n(),2)+2*this->n()+1)-pow(this->n(),2)*T)
-*log(this->b()*(pow(2,0.5)+1)+(this->W()/rho))/(4*this->b()*this->Tcrit());
-}
-
-
-//(molar values)
-//needed to calculate the entropy
-inline scalar pengRobinson::integral_dpdT_dv(const scalar rho,const scalar T) const
-{
-
-return (pow(2,0.5)*this->a()*this->n()*(this->n()+1)*pow(T/this->Tcrit(),0.5)/(4*this->b()*T)
--pow(2,0.5)*this->a()*pow(this->n(),2)/(4*this->b()*this->Tcrit()))
-*log(this->b()*(1-pow(2,0.5))+(this->W()/rho))
-+this->RR*log((this->W()/rho)-this->b())
-+(pow(2,0.5)*this->a()*pow(this->n(),2)/(4*this->b()*this->Tcrit())
--pow(2,0.5)*this->a()*this->n()*(this->n()+1)*pow(T/this->Tcrit(),0.5)/(4*this->b()*T))
-*log(this->b()*(pow(2,0.5)+1)+(this->W()/rho));
-
-}
-
-
-
-
-//* * * * * * * * * * * * * second order Derivative based functions * * * * * * * * * * * //
-
-
-
-//(molar values)
-inline scalar pengRobinson::d2pdT2(const scalar rho,const scalar T) const
-{
-return this->a()*this->n()*(this->n()+1)*pow((T/this->Tcrit()),0.5)/(2*pow(T,2)*(pow(this->b(),2)-2*this->b()*(this->W()/rho)-pow((this->W()/rho),2)));
-}
-
-
-
-
-//(molar values)
-inline scalar pengRobinson::d2pdv2(const scalar rho,const scalar T) const
-{
-return 2*(2*this->a()*this->n()*this->Tcrit()*(this->b()-(this->W()/rho))*(this->n()+1)*
-(5*pow(this->b(),4)-4*pow(this->b(),3)*(this->W()/rho)-4*pow(this->b(),2)*pow((this->W()/rho),2)
-+3*pow((this->W()/rho),4))*pow(T/this->Tcrit(),0.5)
-+this->Tcrit()*(this->RR*T*(pow(this->b(),6)-6*pow(this->b(),5)*(this->W()/rho)+9*pow(this->b(),4)*pow((this->W()/rho),2)
-+4*pow(this->b(),3)*pow((this->W()/rho),3)-9*pow(this->b(),2)*pow((this->W()/rho),4)-6*this->b()*pow((this->W()/rho),5)
--pow((this->W()/rho),6))
--this->a()*pow((this->n()+1),2)*(5*pow(this->b(),5)-9*pow(this->b(),4)*(this->W()/rho)
-+4*pow(this->b(),2)*pow((this->W()/rho),3)+3*this->b()*pow((this->W()/rho),4)-3*pow((this->W()/rho),5)))
--this->a()*pow(this->n(),2)*T*(5*pow(this->b(),5)-9*pow(this->b(),4)*(this->W()/rho)+4*pow(this->b(),2)*pow((this->W()/rho),3)
-+3*this->b()*pow((this->W()/rho),4)-3*pow((this->W()/rho),5)))
-/(this->Tcrit()*pow((pow(this->b(),2)-2*this->b()*(this->W()/rho)-pow((this->W()/rho),2)),3)*pow(((this->W()/rho)-this->b()),3));
-}
-
-
-
-
-//(molar values)
-//using second order implicit differentiation 
-inline scalar pengRobinson::d2vdT2(const scalar rho, const scalar T) const
-{
-return -(pow(this->dpdT(rho,T),2)*this->d2pdv2(rho,T)+   pow(this->dpdv(rho,T),2)  *this->d2pdT2(rho,T)- 2*this->dpdv(rho,T)*this->dpdT(rho,T)*this->d2pdvdT(rho,T))/(   pow(this->dpdv(rho,T),3));
-}
-
-
-
-
-//(molar values)
-inline scalar pengRobinson::d2pdvdT(const scalar rho, const scalar T) const
-{
-return -2*this->a()*this->n()*(this->b()+(this->W()/rho))*(this->n()+1)*pow(T/this->Tcrit(),0.5)
-/(T*pow((pow(this->b(),2)-2*this->b()*(this->W()/rho)-pow((this->W()/rho),2)),2))+(2*this->a()*pow(this->n(),2)*(this->b()+(this->W()/rho)))/(this->Tcrit()*pow(pow(this->b(),2)-2*this->b()*(this->W()/rho)-pow((this->W()/rho),2),2))-this->RR/(pow(this->b()-(this->W()/rho),2));
-}
-
-
-
-
-// the result of this intergal is needed for the nasa based cp polynomial 
-//(molar values)
-inline scalar pengRobinson::integral_d2pdT2_dv(const scalar rho,const scalar T) const 
-{
-return pow(2,0.5)*this->a()*this->n()*(this->n()+1)*pow(T/this->Tcrit(),0.5)
-*log(this->b()*(pow(2,0.5)+1)+(this->W()/rho))/(8*this->b()*pow(T,2))
--pow(2,0.5)*this->a()*this->n()*(this->n()+1)*pow(T/this->Tcrit(),0.5)
-*log(this->b()*(1-pow(2,0.5))+(this->W()/rho))/(8*this->b()*pow(T,2));
-}
-
-
-
-
-
-//* * * * * * * * * * * * * thermodynamic properties * * * * * * * * * * * //
-
-
-//Isobar expansion Coefficent beta = 1/v (dv/dt) at constant p
-//(molar values)
-inline scalar pengRobinson::isobarExpCoef(const scalar rho,const scalar T)  const
-{
- return this->dvdT(rho, T)*rho/this->W();
-}
-
-//isothemal compressiblity kappa 
-//(molar values)
-inline scalar pengRobinson::isothermalCompressiblity(const scalar rho,const scalar T) const
-{
-return this->isobarExpCoef(rho, T)/this->dpdT(rho, T);
-//also possible : return -this->dvdp(rho,T)*rho/this->W();
-}
-
-
-
-
-
-
-
-
-
 
 // Construct from components
 inline pengRobinson::pengRobinson
@@ -290,25 +54,25 @@ inline pengRobinson::pengRobinson
     specie(sp),    
     pcrit_(pcrit),
     Tcrit_(Tcrit),
-    azentricFactor_(azentricFactor)
+    azentricFactor_(azentricFactor),
-{
+    a_(0.457235*pow(this->RR,2)*pow(Tcrit_,2)/pcrit_),
-	rhostd_=this->rho(Pstd,Tstd,Pstd*this->W()/(Tstd*this->R()));	// Starting GUESS for the density by ideal gas law
+    b_(0.077796*this->RR*Tcrit_/pcrit_), 
-}
+    n_(0.37464+1.54226*azentricFactor_-0.26992*pow(azentricFactor_,2)), 
-
+    rhostd_(this->rho(Pstd,Tstd,Pstd*this->W()/(Tstd*this->R()))) 
-
+{}
-// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
 // Construct as named copy
-
+inline pengRobinson::pengRobinson(const word& name, const pengRobinson& pg)
-inline pengRobinson::pengRobinson(const word& name,  pengRobinson& pg)
 :
     specie(name, pg),
     pcrit_(pg.pcrit_),
     Tcrit_(pg.Tcrit_),
-    azentricFactor_(pg.azentricFactor_)
+    azentricFactor_(pg.azentricFactor_),
-{
+    a_(pg.a_),
-  pg.rhostd_=this->rho(Pstd,Tstd,  (Pstd*this->W()/(Tstd*this->R())));	// Starting GUESS for the density by ideal gas law
+    b_(pg.b_),
-}
+    n_(pg.n_),
+    rhostd_(pg.rhostd_)
+{}
 
 
 // Construct and return a clone
@@ -325,110 +89,333 @@ inline autoPtr<pengRobinson> pengRobinson::New(Istream& is)
 }
 
 
-// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+// * * * * * * * * * * * * *  Member Functions  * * * * * * * * * * * //
+
+inline scalar pengRobinson::rhostd()const
+{
+    return rhostd_;
+}
+
+
+//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_*pow((1+n_*(1-pow((T/Tcrit_),0.5))),2))
+        /(pow(Vm,2)+2*b_*Vm-pow(b_,2));
+}
+
+
+//Real deviative dp/dv at constant temperature 
+//(molar values)
+inline scalar pengRobinson::dpdv(const scalar rho, const scalar T) const
+{
+    scalar Vm = this->W()/rho;
+    return 
+    -(
+        4*a_*n_*Tcrit_*(b_-Vm)*(pow(b_,2)
+        -pow(Vm,2))*(n_+1)*pow((T/Tcrit_),0.5)
+        +Tcrit_*(-2*a_*pow((n_+1),2)*(pow(b_,3)-pow(b_,2)*Vm
+        -b_*pow(Vm,2)+pow(Vm,3))
+        +this->RR*T*(pow(b_,4)-4*pow(b_,3)*Vm
+        +2*pow(b_,2)*pow(Vm,2)
+        +4*b_*pow(Vm,3)+pow(Vm,4)))
+        -2*a_*pow(n_,2)*T*(pow(b_,3)-pow(b_,2)*Vm
+        -b_*pow(Vm,2)+pow(Vm,3))
+     )
+     /(Tcrit_*pow((b_-Vm),2)*pow((pow(b_,2)
+        -2*b_*Vm-pow(Vm,2)),2));
+}
+
+
+//Real deviative dp/dT at constant molar volume 
+//(molar values)
+inline scalar pengRobinson::dpdT(const scalar rho, const scalar T) const
+{
+    scalar Vm = this->W()/rho;
+    return 
+    (
+        -a_*n_*(n_+1)*pow((T/Tcrit_),0.5)
+            /(T*(pow(b_,2)
+        -2*b_*Vm-pow(Vm,2)))
+        +a_*pow(n_,2)/(Tcrit_*(pow(b_,2)
+        -2*b_*Vm-pow(Vm,2)))
+        -this->RR/(b_-Vm)
+    );
+}
+
+
+//Real  deviative dv/dT at constant pressure 
+//by using implicit differentiation
+//(molar values)
+inline scalar pengRobinson::dvdT(const scalar rho,const scalar T) const
+{
+    return (-1)*this->dpdT(rho,T)/this->dpdv(rho,T);
+}
+
+
+//(molar values)
+inline scalar pengRobinson::dvdp(const scalar rho,const scalar T)  const 
+{
+    return 1/this->dpdv(rho,T);
+}
+
+
+//(molar values)
+//needed to calculate the internal energy
+inline scalar pengRobinson::integral_p_dv(const scalar rho,const scalar T) const
+{
+    scalar Vm = this->W()/rho;
+    return  pow(2,0.5)*a_*(2*n_*Tcrit_*(n_+1)*pow(T/Tcrit_,0.5)
+        -Tcrit_*(pow(n_,2)+2*n_+1)-pow(n_,2)*T)
+            *log(b_*(1-pow(2,0.5))+Vm)/(4*b_*Tcrit_)
+        +this->RR*T*log(Vm-b_)
+        -pow(2,0.5)*a_*(2*n_*Tcrit_*(n_+1)*pow(T/Tcrit_,0.5)
+        -Tcrit_*(pow(n_,2)+2*n_+1)
+        -pow(n_,2)*T)*log(b_*(pow(2,0.5)+1)+Vm)
+            /(4*b_*Tcrit_);
+}
+
+
+//(molar values)
+//needed to calculate the entropy
+inline scalar pengRobinson::integral_dpdT_dv(const scalar rho,const scalar T) const
+{
+    scalar Vm = this->W()/rho;
+    return (pow(2,0.5)*a_*n_*(n_+1)*pow(T/Tcrit_,0.5)/(4*b_*T)
+        -pow(2,0.5)*a_*pow(n_,2)/(4*b_*Tcrit_))
+            *log(b_*(1-pow(2,0.5))+Vm)
+        +this->RR*log(Vm-b_)
+        +(pow(2,0.5)*a_*pow(n_,2)/(4*b_*Tcrit_)
+        -pow(2,0.5)*a_*n_*(n_+1)*pow(T/Tcrit_,0.5)/(4*b_*T))
+            *log(b_*(pow(2,0.5)+1)+Vm);
+}
+
+
+//(molar values)
+inline scalar pengRobinson::d2pdT2(const scalar rho,const scalar T) const
+{
+    scalar Vm = this->W()/rho;
+    return 
+    a_*n_*(n_+1)*pow((T/Tcrit_),0.5)
+        /(2*pow(T,2)*(pow(b_,2)
+    -2*b_*Vm-pow(Vm,2)));
+}
+
+
+//(molar values)
+inline scalar pengRobinson::d2pdv2(const scalar rho,const scalar T) const
+{
+    scalar Vm = this->W()/rho;
+    return 
+    2*(
+        2*a_*n_*Tcrit_*(b_-Vm)*(n_+1)*
+        (
+            5*pow(b_,4)-4*pow(b_,3)*Vm
+            -4*pow(b_,2)*pow(Vm,2)
+            +3*pow(Vm,4)
+        )
+        *pow(T/Tcrit_,0.5)
+        +Tcrit_*
+        (
+            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)
+                -9*pow(b_,2)*pow(Vm,4)
+                -6*b_*pow(Vm,5)
+                -pow(Vm,6)
+            )
+            -a_*pow((n_+1),2)*
+            (
+                 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)
+            )
+         )
+         -a_*pow(n_,2)*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)
+         )
+     )
+     /
+     (
+         Tcrit_*
+         pow((pow(b_,2)
+             -2*b_*Vm-pow(Vm,2)),3)
+         *pow((Vm-b_),3)
+     );
+}
+
+
+//(molar values)
+//using second order implicit differentiation 
+inline scalar pengRobinson::d2vdT2(const scalar rho, const scalar T) const
+{
+    return 
+    -(
+          pow(this->dpdT(rho,T),2)*this->d2pdv2(rho,T)
+        + pow(this->dpdv(rho,T),2)*this->d2pdT2(rho,T)
+        - 2*this->dpdv(rho,T)*this->dpdT(rho,T)*this->d2pdvdT(rho,T)
+    )
+    /(pow(this->dpdv(rho,T),3));
+}
+
+
+//(molar values)
+inline scalar pengRobinson::d2pdvdT(const scalar rho, const scalar T) const
+{
+    scalar Vm = this->W()/rho;
+    return -2*a_*n_*(b_+Vm)*(n_+1)*pow(T/Tcrit_,0.5)
+        /(T*pow((pow(b_,2)-2*b_*Vm-pow(Vm,2)),2))
+    +(2*a_*pow(n_,2)*(b_+Vm))
+        /(Tcrit_*pow(pow(b_,2)-2*b_*Vm-pow(Vm,2),2))
+    -this->RR/(pow(b_-Vm,2));
+}
+
+
+// the result of this intergal is needed for the nasa based cp polynomial 
+//(molar values)
+inline scalar pengRobinson::integral_d2pdT2_dv(const scalar rho,const scalar T) const 
+{
+    scalar Vm = this->W()/rho;
+    return pow(2,0.5)*a_*n_*(n_+1)*pow(T/Tcrit_,0.5)
+        *log(b_*(pow(2,0.5)+1)+Vm)/(8*b_*pow(T,2))
+    -pow(2,0.5)*a_*n_*(n_+1)*pow(T/Tcrit_,0.5)
+        *log(b_*(1-pow(2,0.5))+Vm)/(8*b_*pow(T,2));
+}
+
+
+//Isobar expansion Coefficent beta = 1/v (dv/dt) at constant p
+//(molar values)
+inline scalar pengRobinson::isobarExpCoef(const scalar rho,const scalar T)  const
+{
+    return this->dvdT(rho, T)*rho/this->W();
+}
+
+
+//isothemal compressiblity kappa (not Thermal conductivity)
+//(molar values)
+inline scalar pengRobinson::isothermalCompressiblity(const scalar rho,const scalar T) const
+{
+    return this->isobarExpCoef(rho, T)/this->dpdT(rho, T);
+    //also possible : return -this->dvdp(rho,T)*rho/this->W();
+}
+
 
 //- Return density [kg/m^3]
-// TO DO Include a max Iteration number loop and abort function 
+inline scalar pengRobinson::rho(
-inline scalar pengRobinson::rho(const scalar p,const scalar T,const scalar rho0) const
+    const scalar p,
+    const scalar T,
+    const scalar rho0
+) const
 {
  
    scalar molarVolumePrevIteration;
    scalar molarVolume;
    int iter=0;
-   int i;
    int maxIter_=400;
    scalar tol_=1e-8;
-   scalar rho1=rhoMax_, rho2=rhoMin_,rho3, f1,f2,f3;
+   scalar rho1=rhoMax_;
+   scalar rho2=rhoMin_;
 
-molarVolume=this->W()/rho0;	 
+   molarVolume=this->W()/rho0;	 
 
    do
    {
-   molarVolumePrevIteration= molarVolume;
+       molarVolumePrevIteration= molarVolume;
-   i=0;
-      	do
-   	{
    
-  	molarVolume=molarVolumePrevIteration-((this->p((this->W()/molarVolumePrevIteration),T)-p)/(this->dpdv((this->W()/       	
+       label i=0;
-  	molarVolumePrevIteration),T)))/(pow(2,i)); 
+       do
+       {
+           molarVolume=molarVolumePrevIteration
+               -(
+                   (this->p((this->W()/molarVolumePrevIteration),T) - p)
+                   /(this->dpdv((this->W()/molarVolumePrevIteration),T))
+               )/pow(2,i);
 
+           i++;
+           if (i>8) 
+           {
+               //using bisection methode as backup,
+               //solution must be between rho=0.001 to rho=1500;
+               for(i=0; i<200; i++)
+               {
+                   scalar f1 = this->p(rho1,T) - p;
+                   scalar f2 = this->p(rho2,T) - p;
+                   scalar rho3 = (rho1 + rho2)/2;
+                   scalar f3 = this->p(rho3,T) - p;
 	
-	i++;
+                   if ((f2 < 0 && f3 > 0) || (f2 > 0 && f3 < 0))
-		if(i>8) 
+                   {
-		{
+                       rho1=rho3;
-		  //using bisection methode as backup, solution must be between rho=0.001 to rho=1500; 
+                   }  
-		  for(i=0;i<200;i++)
+                   else if ((f1 < 0 && f3 > 0)||(f1 > 0 && f3 < 0))
-		  {
+                   {
+                       rho2=rho3;
+                   }		 
+                   else	
+                   {
+                       rho2=(rho2 + rho3)/2;
+                   }
 
-                  	f1= (this->p(rho1,T)-p);
+                   if(mag(f3) < p*tol_)
-                  	f2= (this->p(rho2,T)-p);
+                   {	  	  
-                  	rho3=(rho1+rho2)/2;
+                       molarVolume=this->W()/rho3;
-		  	f3=(this->p(rho3,T)-p);
+                       molarVolumePrevIteration=this->W()/rho3;
-	
+                       break; 
-			if ((f2<0 && f3>0)||(f2>0 &&f3<0))	
+                   }
-		  	{
+                   else
-		 		rho1=rho3;
+                   {
-			}  
+                       molarVolumePrevIteration=this->W()/rho3;
- 		  	else if ((f1<0 && f3>0)||(f1>0 &&f3<0))
+                   }
-                  	{
+               }		  
-		 		rho2=rho3;
+           }
-			}
+       }
-		 
+       while
-			else	
+       (
-			{
+           mag(this->p((this->W()/molarVolume),T) - p)
-                        	rho2=(rho2+rho3)/2;
+         > mag(this->p((this->W()/molarVolumePrevIteration),T) - p)
-			}
+       );
-
-			if(mag(f3)<p*tol_)
-                  	{	  	  
-				molarVolume=this->W()/rho3;
-                 	 	molarVolumePrevIteration=this->W()/rho3;
-				break; 
-		  	}
-			else
-			{
-				molarVolumePrevIteration=this->W()/rho3;
-			}
-		    }		  
-  		}
-
-
-        }
-        while(mag(this->p((this->W()/molarVolume),T)-p)>mag(this->p((this->W()/molarVolumePrevIteration),T)-p));
-
-           if (iter++ > maxIter_)
-        {
-            FatalErrorIn
-            (
-                "inline scalar pengRobinson::rho(const scalar p,const scalar T,const scalar rho0) const "
-            )   << "Maximum number of iterations exceeded"
-                << abort(FatalError);
-        }
 
+       if (iter++ > maxIter_)
+       {
+           FatalErrorIn
+           (
+               "inline scalar redlichKwong::rho(const scalar p, const scalar T, const scalar rho0) const "
+           )   << "Maximum number of iterations exceeded"
+               << abort(FatalError);
+       }
    }
-   while(mag(molarVolumePrevIteration-molarVolume)>tol_*(this->W()/rho0));
+   while(mag(molarVolumePrevIteration-molarVolume) > tol_*(this->W()/rho0));
-
-
-
-
 
    return this->W()/molarVolume;
- 
 }
 
+
 //- Return density [kg/m^3]on 
 inline scalar pengRobinson::rho(const scalar p,const scalar T) const
 {
-	scalar rho0=p/(this->R()*T); //using perfect gas equation as starting point
+    // using perfect gas equation as starting point
-	return rho(p,T,rho0); 
+    return rho(p,T,p/(this->R()*T)); 
 }
 
-//- Return compressibility drho/dp [s^2/m^2] 
+
+//- Return compressibility drho/dp at T=constant [s^2/m^2] 
 inline scalar pengRobinson::psi(const scalar rho, const scalar T) const
 {
     return -this->dvdp(rho,T)*pow(rho,2)/this->W();
 }
 
+
 //- Return compression factor []
 inline scalar pengRobinson::Z( const scalar p, const scalar T,const scalar rho0) const
 {

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

@@ -63,8 +63,8 @@ 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
+        // Starting GUESS for the density by ideal gas law
-    rhostd_(this->rho(Pstd, Tstd, Pstd*this->W()/(Tstd*this->R())))
+        rhostd_(this->rho(Pstd, Tstd, Pstd*this->W()/(Tstd*this->R())))
 { 
     is.check("redlichKwong::redlichKwong(Istream& is)");
 }

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

@@ -53,7 +53,7 @@ namespace Foam
 {
 
 /*---------------------------------------------------------------------------*\
-                           Class perfectGas Declaration
+                           Class redlichKwong Declaration
 \*---------------------------------------------------------------------------*/
 
 class redlichKwong
@@ -62,27 +62,22 @@ class redlichKwong
 {
     // Private data
 
-        scalar pcrit_;
+    scalar pcrit_;
+    scalar Tcrit_; 
 
- 	    scalar Tcrit_; 
+    //-Redlich Kwong factors
-
+    scalar a_;
-	//-Redlich Kwong factors
+    scalar b_;
-
-	    scalar a_;
-
-        scalar b_;
 
     //- Density @STD, initialise after a, b!
-
+    scalar rhostd_;
-        scalar rhostd_;
 
     // Static Private data
 
-	//should be read from the fvSolution file where rhoMax and rhoMin values must be defined (for rhoSimpleFoam)
+    //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 rhoMax_;
-
+    static const scalar rhoMin_;
-        static const scalar rhoMin_;
 
 
 public:
@@ -96,7 +91,7 @@ public:
             const specie& sp, 
             scalar pcrit, 
             scalar Tcrit
-		);
+        );
 
         //- Construct from Istream
         redlichKwong(Istream&);
@@ -117,7 +112,7 @@ public:
         inline scalar p(const scalar rho, const scalar T) const;
 
 
-        // Derivatives
+        //first order derivatives 
         inline scalar dpdv(const scalar rho, const  scalar T) const;
 
         inline scalar dpdT(const scalar rho, const scalar T) const;
@@ -147,7 +142,7 @@ public:
         // Used for Entropy
         inline scalar integral_dpdT_dv(const scalar rho, const scalar T) const;
 
-        // Not Used At The Moment
+        // second order derivatives, not Used At The Moment
         inline scalar d2pdv2(const scalar rho, const scalar T) const;
 
         inline scalar d2pdT2(const scalar rho, const scalar T) const;
@@ -161,13 +156,13 @@ public:
         inline scalar rho
         (
             const scalar p,
-            const  scalar T,
+            const scalar T,
             const scalar rho0
         ) const;
 
         inline scalar rho(const scalar p, const scalar T) const;
 
-        //- Return compressibility drho/dp [s^2/m^2]  
+        //- Return compressibility drho/dp at T=constant [s^2/m^2]  
         inline scalar psi(const scalar rho, const scalar T) const;
 
         //- Return compression factor []

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

@@ -252,7 +252,7 @@ inline scalar redlichKwong::isobarExpCoef
 }
 
 
-//isothemal compressiblity kappa 
+//isothemal compressiblity kappa (not Thermal conductivity)
 //(molar values)
 inline scalar redlichKwong::isothermalCompressiblity
 (
@@ -360,11 +360,11 @@ inline scalar redlichKwong::rho
 inline scalar redlichKwong::rho(const scalar p, const scalar T) const
 {
     // using perfect gas equation as starting point
-	return rho(p,T,p/(this->R()*T)); 
+    return rho(p,T,p/(this->R()*T));
 }
 
 
-//- Return compressibility drho/dp [s^2/m^2]
+//- Return compressibility drho/dp at T=constant [s^2/m^2]
 inline scalar redlichKwong::psi(const scalar rho, const scalar T) const
 {
     return -this->dvdp(rho,T)*pow(rho,2)/this->W();

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

@@ -41,6 +41,17 @@ 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)
@@ -48,12 +59,14 @@ soaveRedlichKwong::soaveRedlichKwong(Istream& is)
     specie(is),
     pcrit_(readScalar(is)),
     Tcrit_(readScalar(is)),
-	azentricFactor_(readScalar(is))
+    azentricFactor_(readScalar(is)),
+    a_(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
+        rhostd_(this->rho(Pstd,Tstd,Pstd*this->W()/(Tstd*this->R())))
 {
     is.check("soaveRedlichKwong::soaveRedlichKwong(Istream& is)"); 	
-    rhostd_=this->rho(Pstd,Tstd,Pstd*this->W()/(Tstd*this->R()));
-    rhoMax_=1500;	
-    rhoMin_=0.001;    	
 }
 
 

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

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 1991-2009 OpenCFD Ltd.
+    \\  /    A nd           | 
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -45,8 +45,7 @@ Germany
 
 #include "specie.H"
 #include "autoPtr.H"
-#include "word.H"
+
-#include "scalar.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -54,7 +53,7 @@ namespace Foam
 {
 
 /*---------------------------------------------------------------------------*\
-                           Class perfectGas Declaration
+                           Class Soave Readlich Kwong Declaration
 \*---------------------------------------------------------------------------*/
 
 class soaveRedlichKwong
@@ -62,7 +61,26 @@ class soaveRedlichKwong
     public specie
 
 {
+    // Private data
 
+    scalar pcrit_;
+    scalar Tcrit_; 
+    scalar azentricFactor_;
+
+    //-Soave Redlich Kwong factors
+    scalar a_;
+    scalar b_;
+    scalar n_;
+
+    //- Density @STD, initialise after a, b!
+    scalar rhostd_;
+
+    // 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_;
 
 public:
   
@@ -81,7 +99,7 @@ public:
         soaveRedlichKwong(Istream&);
 
         //- Construct as named copy
-        inline soaveRedlichKwong(const word& name, soaveRedlichKwong&);
+        inline soaveRedlichKwong(const word& name,const soaveRedlichKwong&);
 
         //- Construct and return a clone
         inline autoPtr<soaveRedlichKwong> clone() const;
@@ -89,50 +107,72 @@ public:
         // Selector from Istream
         inline static autoPtr<soaveRedlichKwong> New(Istream& is);
   
-        //Member Variabels
-	    scalar pcrit_;
- 	    scalar Tcrit_; 
-            scalar rhostd_;
-	    scalar azentricFactor_;
-	    scalar rhoMax_;	//should be read from the fvSolution file where rhoMax and rhoMin values must be define ( for rhoSimpleFoam)
-    	    scalar rhoMin_;
-
-
     // Member functions
-        inline scalar pcrit() const;
-        inline scalar Tcrit() const;
-	inline scalar azentricFactor() const;
 	inline scalar rhostd()const;
+
         inline scalar p(const scalar rho, const scalar T) const;
 
-	//-Redlich Kwong factors
-     	inline scalar a() const;
-	inline scalar b() const;
-	inline scalar n() const;
 
-        //derivatives 
+        //first order derivatives  
         inline scalar dpdv(const scalar rho,const  scalar T) const;
-        inline scalar dpdT(const scalar rho, const scalar T) const;
-        inline scalar dvdT(const scalar rho,const  scalar T) const;
-        inline scalar dvdp(const scalar rho, const scalar T) const; 
-        inline scalar isobarExpCoef(const scalar rho,const  scalar T) const; 
-        inline scalar isothermalCompressiblity(const scalar rho,const  scalar T) const; 
-        inline scalar integral_d2pdT2_dv(const scalar rho,const  scalar T) const ; // Used for cv
-        inline scalar d2pdv2(const scalar rho,const  scalar T) const;    // not Used At The Moment
-	inline scalar d2pdT2(const scalar rho,const  scalar T) const;    // not Used At The Moment
-        inline scalar d2pdvdT(const scalar rho,const  scalar T) const;    // not Used At The Moment
-	inline scalar d2vdT2(const scalar rho,const  scalar T) const;    // not Used At The Moment
-        inline scalar integral_p_dv(const scalar rho,const  scalar T) const;     //Used for internal Energy
-        inline scalar integral_dpdT_dv(const scalar rho,const  scalar T) const;  //Used for Entropy
 
-        //- Return density [kg/m^3] // rho0 is the starting point of the newton solver used to calculate rho
+        inline scalar dpdT(const scalar rho, const scalar T) const;
-        inline scalar rho(const scalar p,const  scalar T,const scalar rho0) const;
+
+        inline scalar dvdT(const scalar rho,const  scalar T) const;
+
+        inline scalar dvdp(const scalar rho, const scalar T) const; 
+
+        inline scalar isobarExpCoef
+        (
+            const scalar rho,
+            const scalar T
+        ) const; 
+
+        inline scalar isothermalCompressiblity
+        (
+            const scalar rho,
+            const scalar T
+        ) const; 
+
+        // Used for cv
+        inline scalar integral_d2pdT2_dv(const scalar rho,const  scalar T) const ; 
+
+        // second order derivatives, not Used At The Moment
+        inline scalar d2pdv2(const scalar rho,const  scalar T) const;    
+
+	inline scalar d2pdT2(const scalar rho,const  scalar T) const;    
+
+        inline scalar d2pdvdT(const scalar rho,const  scalar T) const;    
+
+	inline scalar d2vdT2(const scalar rho,const  scalar T) const;    
+
+        //Used for internal Energy
+        inline scalar integral_p_dv(const scalar rho,const  scalar T) const;    
+
+        //Used for Entropy
+        inline scalar integral_dpdT_dv(const scalar rho,const  scalar T) const; 
+
+        //- Return density [kg/m^3] 
+        // rho0 is the starting point of the newton solver used to calculate rho
+        inline scalar rho
+        (
+            const scalar p,
+            const scalar T,
+            const scalar rho0
+        ) const;
+
         inline scalar rho(const scalar p,const  scalar T) const;
-        //- Return compressibility rho/p [s^2/m^2]   
+
+        //- Return compressibility drho/dp at T=constant [s^2/m^2] 
         inline scalar psi(const scalar rho, const scalar T) const;
 
         //- Return compression factor []
-        inline scalar Z(const scalar p,const  scalar T,const scalar rho0) const;
+        inline scalar Z
+        (
+            const scalar p,
+            const scalar T,
+            const scalar rho0
+        ) const;
 
 
     // Member operators

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

@@ -38,225 +38,7 @@ Germany
 namespace Foam
 {
 
-
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
-
-
-
-// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
-
-inline scalar soaveRedlichKwong::pcrit()const
-{
-return pcrit_;
-}
-
-
-
-inline scalar soaveRedlichKwong::Tcrit()const
-{
-return Tcrit_;
-}
-
-
-inline scalar soaveRedlichKwong::rhostd()const
-{
-return rhostd_;
-}
-
-
-// Returns the Azentric Factor (Acentric Factor)
-inline scalar soaveRedlichKwong::azentricFactor() const
-{
-return azentricFactor_;
-}
-
-//returns the pressure for a given density and temperature    
-inline scalar soaveRedlichKwong::p(const scalar rho,const scalar T) const
-{
-return  (this->RR*T/((this->W()/rho)-this->b())-this->a()*pow((1+this->n()*(1-pow((T/this->Tcrit()),0.5))),2)/((this->W()/rho)*((this->W()/rho)+this->b())));
-}
-
-
-
-// Faktor a of the soave redlich Kwong equation of state  
-//(molar values)
-inline scalar soaveRedlichKwong::a() const
-{
-return 0.42747*pow(this->RR,2)*pow(this->Tcrit(),2)/(this->pcrit());
-}
-
-
-
-// Faktor b of the soave redlich Kwong equation of state 
-//(molar values)
-inline scalar soaveRedlichKwong::b() const
-{
-return 0.08664*this->RR*this->Tcrit()/this->pcrit(); 
-}
-
-// Faktor n of the soave redlich Kwong equation of state  
-//(molar values)
-inline scalar soaveRedlichKwong::n() const
-{
-return 0.48508+1.55171*this->azentricFactor()-0.15613*pow(this->azentricFactor(),2);
-}
-
-
-//* * * * * * * * * * * * * Derivatives * * * * * * * * * * * //
-
-
-
-//Real deviative dp/dv at constant temperature 
-//(molar values)
-inline scalar soaveRedlichKwong::dpdv(const scalar rho, const scalar T) const
-{
- return -(2*this->a()*this->n()*this->Tcrit()*(this->b()-(this->W()/rho))*(pow(this->b(),2)+this->b()*(this->W()/rho)-2*pow((this->W()/rho),2))*(this->n()+1)*pow((T/this->Tcrit()),0.5)+this->Tcrit()*(this->RR*T*pow((this->W()/rho),2)*(pow(this->b(),2)+2*this->b()*(this->W()/rho)+pow((this->W()/rho),2))-this->a()*(pow(this->b(),3)-3*this->b()*pow((this->W()/rho),2)+2*pow((this->W()/rho),3))*pow((this->n()+1),2))-this->a()*pow(this->n(),2)*T*(pow(this->b(),3)-3*this->b()*pow((this->W()/rho),2)+2*pow((this->W()/rho),3)))/(pow((this->W()/rho),2)*this->Tcrit()*pow((this->b()+(this->W()/rho)),2)*pow((this->b()-(this->W()/rho)),2));
-}
-
-
-
-
-//Real deviative dp/dT at constant molar volume 
-//(molar values)
-inline scalar soaveRedlichKwong::dpdT(const scalar rho, const scalar T) const
-{
-return this->a()*this->n()*(this->n()+1)*pow((T/this->Tcrit()),0.5)/(T*(this->W()/rho)*(this->b()+(this->W()/rho)))-this->a()*pow(this->n(),2)/((this->W()/rho)*this->Tcrit()*(this->b()+(this->W()/rho)))-this->RR/(this->b()-(this->W()/rho));
-}
-
-
-
-
-//Real  deviative dv/dT at constant pressure
-//using implicit differentiation
-// (molar values)
-inline scalar soaveRedlichKwong::dvdT(const scalar rho,const scalar T) const
-{
-return (-1)*this->dpdT(rho,T)/this->dpdv(rho,T);
-}
-
-
-
-
-//Real  deviative dv/dp at constant temperature 
-//(molar values)
-inline scalar soaveRedlichKwong::dvdp(const scalar rho,const scalar T)  const 
-{
-return 1/this->dpdv(rho,T);
-}
-
-
-
-//needed to calculate the internal energy 
-//(molar values)
-inline scalar soaveRedlichKwong::integral_p_dv(const scalar rho,const scalar T) const
-{
-return this->RR*T*log((this->W()/rho)-this->b())
--(this->a()*(2*this->n()*this->Tcrit()*(this->n()+1)*pow(T/this->Tcrit(),0.5)-this->Tcrit()*(pow(this->n(),2)+2*this->n()+1)-pow(this->n(),2)*T)*log(this->b()+(this->W()/rho)))/(this->b()*this->Tcrit())
-+this->a()*(2*this->n()*this->Tcrit()*(this->n()+1)*pow(T/this->Tcrit(),0.5)-this->Tcrit()*(pow(this->n(),2)+2*this->n()+1)-pow(this->n(),2)*T)*log((this->W()/rho))/(this->b()*this->Tcrit());
-}
-
-
-
-//needed to calculate the entropy
-//(molar values)
-inline scalar soaveRedlichKwong::integral_dpdT_dv(const scalar rho,const scalar T) const
-{
-return this->RR*log((this->W()/rho)-this->b())+(pow(this->n(),2)*this->a()/(this->b()*this->Tcrit())-this->a()*this->n()*(this->n()+1)*pow((T/this->Tcrit()),0.5)/(this->b()*T))*log(this->b()+(this->W()/rho))+(this->a()*this->n()*(this->n()+1)*pow((T/this->Tcrit()),0.5)/(this->b()*T)-this->a()*pow(this->n(),2)/(this->b()*this->Tcrit()))*log((this->W()/rho));
-}
-
-
-
-
-//* * * * * * * * * * * * * second order Derivative based functions * * * * * * * * * * * //
-
-
-
-//(molar values)
-inline scalar soaveRedlichKwong::d2pdT2(const scalar rho,const scalar T) const
-{
-return -this->a()*this->n()*(this->n()+1)*pow(T/this->Tcrit(),0.5)/(2*pow(T,2)*(this->W()/rho)*(this->b()+(this->W()/rho)));
-}
-
-
-
-
-//(molar values) 
-inline scalar soaveRedlichKwong::d2pdv2(const scalar rho,const scalar T) const
-{
- return 2*(2*this->a()*this->n()*this->Tcrit()*(this->b()-(this->W()/rho))
-    *(pow(this->b(),4)+pow(this->b(),3)*(this->W()/rho)-2*pow(this->b(),2)*pow((this->W()/rho),2)
-	-3*this->b()*pow((this->W()/rho),3)+3*pow((this->W()/rho),4))*(this->n()+1)*pow(T/this->Tcrit(),0.5)
-	-this->Tcrit()*(this->a()*(pow(this->b(),5)-3*pow(this->b(),3)*pow((this->W()/rho),2)-pow(this->b(),2)*pow((this->W()/rho),3)
-	+6*this->b()*pow((this->W()/rho),4)-3*pow((this->W()/rho),5))*pow((this->n()+1),2)
-	+this->RR*T*pow((this->W()/rho),3)*(pow(this->b(),3)+3*pow(this->b(),2)*(this->W()/rho)+3*this->b()*pow((this->W()/rho),2)+pow((this->W()/rho),3)))
-	-this->a()*pow(this->n(),2)*T*(pow(this->b(),5)-3*pow(this->b(),3)*pow((this->W()/rho),2)-pow(this->b(),2)*pow((this->W()/rho),3)
-	+6*this->b()*pow((this->W()/rho),4)-3*pow((this->W()/rho),5)))
-	/(pow((this->W()/rho),3)*this->Tcrit()*pow(((this->W()/rho)+this->b()),3)*pow((this->b()-(this->W()/rho)),3));
-}
-
-
-
-
-//(molar values)
-// using second Order implicit differentiation 
-inline scalar soaveRedlichKwong::d2vdT2(const scalar rho, const scalar T) const
-{
-return -(pow(this->dpdT(rho,T),2)*this->d2pdv2(rho,T)+   pow(this->dpdv(rho,T),2)  *this->d2pdT2(rho,T)- 2*this->dpdv(rho,T)*this->dpdT(rho,T)*this->d2pdvdT(rho,T))/(   pow(this->dpdv(rho,T),3));
-}
-
-
-
-
-//(molar values)
-inline scalar soaveRedlichKwong::d2pdvdT(const scalar rho, const scalar T) const
-{
-return -this->a()*this->n()*(this->b()+2*(this->W()/rho))*(this->n()+1)*pow(T/this->Tcrit(),0.5)
-/(T*pow((this->W()/rho),2)*pow((this->b()+(this->W()/rho)),2))
--(this->RR*pow((this->W()/rho),2)*this->Tcrit()*(pow(this->b(),2)+2*this->b()*(this->W()/rho)+pow((this->W()/rho),2))
-  -this->a()*pow(this->n(),2)*(pow(this->b(),3)-3*this->b()*pow((this->W()/rho),2)+2*pow((this->W()/rho),3)))
-  /(pow((this->W()/rho),2)*this->Tcrit()*pow((this->b()+(this->W()/rho)),2)*pow((this->b()-(this->W()/rho)),2));
-}
-
-
-
-
-// the result of this intergal is needed for the nasa based cp polynomial 
-//(molar values)
-inline scalar soaveRedlichKwong::integral_d2pdT2_dv(const scalar rho,const scalar T) const 
-{
-return this->a()*this->n()*(this->n()+1)*pow(T/this->Tcrit(),0.5)*log(this->b()+(this->W()/rho))
-/(2*this->b()*pow(T,2))
--this->a()*this->n()*(this->n()+1)*pow(T/this->Tcrit(),0.5)*log((this->W()/rho))/(2*this->b()*pow(T,2));
-}
-
-
-
-
-
-//* * * * * * * * * * * * * thermodynamic properties * * * * * * * * * * * //
-
-
-//Isobar expansion Coefficent beta = 1/v (dv/dt) at constant p
-//(molar values)
-inline scalar soaveRedlichKwong::isobarExpCoef(const scalar rho,const scalar T)  const
-{
-return this->dvdT(rho, T)*rho/this->W();
-}
-
-//isothemal compressiblity kappa 
-//(molar values)
-inline scalar soaveRedlichKwong::isothermalCompressiblity(const scalar rho,const scalar T) const
-{
-return this->isobarExpCoef(rho, T)/this->dpdT(rho, T);
-}
-
-
-
-
-
-
-
-
-
 
 // Construct from components
 inline soaveRedlichKwong::soaveRedlichKwong
@@ -264,31 +46,36 @@ inline soaveRedlichKwong::soaveRedlichKwong
     const specie& sp,
     scalar pcrit, 
     scalar Tcrit,
-	scalar azentricFactor
+    scalar azentricFactor
     
 )
 :
     specie(sp),    
     pcrit_(pcrit),
     Tcrit_(Tcrit),
-	azentricFactor_(azentricFactor)
+    azentricFactor_(azentricFactor),
-{
+    a_(0.42747*pow(this->RR,2)*pow(Tcrit_,2)/(pcrit_)),
-	rhostd_=this->rho(Pstd,Tstd,Pstd*this->W()/(Tstd*this->R()));
+    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
+        rhostd_(this->rho(Pstd,Tstd,Pstd*this->W()/(Tstd*this->R())))
+{}
+
 
 
-// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
 // Construct as named copy
-inline soaveRedlichKwong::soaveRedlichKwong(const word& name,  soaveRedlichKwong& pg)
+inline soaveRedlichKwong::soaveRedlichKwong(const word& name,const  soaveRedlichKwong& pg)
 :
     specie(name, pg),
     pcrit_(pg.pcrit_),
     Tcrit_(pg.Tcrit_),
-	azentricFactor_(pg.azentricFactor_)
+    azentricFactor_(pg.azentricFactor_),
-{
+    a_(pg.a_),
-  pg.rhostd_=this->rho(Pstd,Tstd,  (Pstd*this->W()/(Tstd*this->R())));	
+    b_(pg.b_),
-}
+    n_(pg.n_),
+    rhostd_(pg.rhostd_)
+{}
 
 
 // Construct and return a clone
@@ -305,10 +92,247 @@ inline autoPtr<soaveRedlichKwong> soaveRedlichKwong::New(Istream& is)
 }
 
 
-// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+// * * * * * * * * * * * * *  Member Functions  * * * * * * * * * * * //
 
-//- Return density [kg/m^3]on 
+inline scalar soaveRedlichKwong::rhostd()const
-inline scalar soaveRedlichKwong::rho(const scalar p,const scalar T,const scalar rho0) const
+{
+    return rhostd_;
+}
+
+
+//returns the pressure for a given density and temperature    
+inline scalar soaveRedlichKwong::p(const scalar rho,const scalar T) const
+{
+    scalar Vm = this->W()/rho;
+    return 
+    (
+        this->RR*T/(Vm-b_)
+        -a_*pow((1+n_*(1-pow((T/Tcrit_),0.5))),2)
+            /(Vm*(Vm+b_))
+    );
+}
+
+
+//Real deviative dp/dv at constant temperature 
+//(molar values)
+inline scalar soaveRedlichKwong::dpdv(const scalar rho, const scalar T) const
+{
+    scalar Vm = this->W()/rho;
+    return 
+    -(
+        2*a_*n_*Tcrit_*(b_-Vm)*(pow(b_,2)
+        +b_*Vm-2*pow(Vm,2))*(n_+1)*pow((T/Tcrit_),0.5)
+        +Tcrit_*
+        (
+            this->RR*T*pow(Vm,2)*(pow(b_,2)
+                +2*b_*Vm+pow(Vm,2))
+        -a_*
+            (pow(b_,3)-3*b_*pow(Vm,2)+
+            2*pow(Vm,3))
+        *pow((n_+1),2)
+        )
+        -a_*pow(n_,2)*T*(pow(b_,3)
+        -3*b_*pow(Vm,2)+2*pow(Vm,3))
+     )
+     /
+     (
+         pow(Vm,2)*Tcrit_*pow((b_+Vm),2)
+             *pow((b_-Vm),2)
+     );
+}
+
+
+//Real deviative dp/dT at constant molar volume 
+//(molar values)
+inline scalar soaveRedlichKwong::dpdT(const scalar rho, const scalar T) const
+{
+    scalar Vm = this->W()/rho;
+    return a_*n_*(n_+1)*pow((T/Tcrit_),0.5)
+        /(T*Vm*(b_+Vm))
+    -a_*pow(n_,2)
+        /(Vm*Tcrit_*(b_+Vm))
+    -this->RR/(b_-Vm);
+}
+
+
+//Real  deviative dv/dT at constant pressure
+//using implicit differentiation
+// (molar values)
+inline scalar soaveRedlichKwong::dvdT(const scalar rho,const scalar T) const
+{
+    return (-1)*this->dpdT(rho,T)/this->dpdv(rho,T);
+}
+
+
+//Real  deviative dv/dp at constant temperature 
+//(molar values)
+inline scalar soaveRedlichKwong::dvdp(const scalar rho,const scalar T)  const 
+{
+    return 1/this->dpdv(rho,T);
+}
+
+
+//needed to calculate the internal energy 
+//(molar values)
+inline scalar soaveRedlichKwong::integral_p_dv(const scalar rho,const scalar T) const
+{
+    scalar Vm = this->W()/rho;
+    return this->RR*T*log(Vm-b_)
+    -(
+        a_*(2*n_*Tcrit_*(n_+1)*pow(T/Tcrit_,0.5)
+        -Tcrit_*(pow(n_,2)+2*n_+1)
+        -pow(n_,2)*T)*log(b_+Vm)
+     )
+     /(b_*Tcrit_)
+     +a_*
+     (
+         2*n_*Tcrit_*(n_+1)*pow(T/Tcrit_,0.5)
+         -Tcrit_*(pow(n_,2)+2*n_+1)-pow(n_,2)*T
+     )
+     *log(Vm)/(b_*Tcrit_);
+}
+
+
+//needed to calculate the entropy
+//(molar values)
+inline scalar soaveRedlichKwong::integral_dpdT_dv(const scalar rho,const scalar T) const
+{
+    scalar Vm = this->W()/rho;
+    return this->RR*log(Vm-b_)
+    +(
+        pow(n_,2)*a_/(b_*Tcrit_)
+        -a_*n_*(n_+1)*pow((T/Tcrit_),0.5)
+            /(b_*T)
+    )*log(b_+Vm)
+    +(a_*n_*(n_+1)*pow((T/Tcrit_),0.5)
+        /(b_*T)
+    -a_*pow(n_,2)/(b_*Tcrit_))*log(Vm);
+}
+
+
+//(molar values)
+inline scalar soaveRedlichKwong::d2pdT2(const scalar rho,const scalar T) const
+{
+    scalar Vm = this->W()/rho;
+    return -a_*n_*(n_+1)*pow(T/Tcrit_,0.5)
+        /(2*pow(T,2)*Vm*(b_+Vm));
+}
+
+
+//(molar values) 
+inline scalar soaveRedlichKwong::d2pdv2(const scalar rho,const scalar T) const
+{
+    scalar Vm = this->W()/rho;
+    return 
+    2*(
+        2*a_*n_*Tcrit_*(b_-Vm)*
+        (
+            pow(b_,4)+pow(b_,3)*Vm
+            -2*pow(b_,2)*pow(Vm,2)
+	    -3*b_*pow(Vm,3)
+            +3*pow(Vm,4)
+        )
+        *(n_+1)*pow(T/Tcrit_,0.5)
+	-Tcrit_*
+        (
+            a_*
+            (
+                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)
+            )*pow((n_+1),2)
+	    +this->RR*T*pow(Vm,3)*
+            (
+                pow(b_,3)
+                +3*pow(b_,2)*Vm
+                +3*b_*pow(Vm,2)
+                +pow(Vm,3)
+            )
+        )
+	-a_*pow(n_,2)*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)
+        )
+    )
+    /(pow(Vm,3)*Tcrit_*pow((Vm+b_),3)
+        *pow((b_-Vm),3));
+}
+
+
+//(molar values)
+// using second Order implicit differentiation 
+inline scalar soaveRedlichKwong::d2vdT2(const scalar rho, const scalar T) const
+{
+    return 
+    -(
+          pow(this->dpdT(rho,T),2)*this->d2pdv2(rho,T)
+        + pow(this->dpdv(rho,T),2)*this->d2pdT2(rho,T)
+        - 2*this->dpdv(rho,T)*this->dpdT(rho,T)*this->d2pdvdT(rho,T)
+    )
+    /(pow(this->dpdv(rho,T),3));
+}
+
+
+//(molar values)
+inline scalar soaveRedlichKwong::d2pdvdT(const scalar rho, const scalar T) const
+{
+    scalar Vm = this->W()/rho;
+    return -a_*n_*(b_+2*Vm)*(n_+1)*pow(T/Tcrit_,0.5)
+        /(T*pow(Vm,2)*pow((b_+Vm),2))
+    -(
+        this->RR*pow(Vm,2)*Tcrit_*(pow(b_,2)
+        +2*b_*Vm+pow(Vm,2))
+        -a_*pow(n_,2)*
+        (
+            pow(b_,3)
+            -3*b_*pow(Vm,2)
+            +2*pow(Vm,3)
+        )
+    )
+  /(pow(Vm,2)*Tcrit_*pow((b_+Vm),2)*pow((b_-Vm),2));
+}
+
+
+// the result of this intergal is needed for the nasa based cp polynomial 
+//(molar values)
+inline scalar soaveRedlichKwong::integral_d2pdT2_dv(const scalar rho,const scalar T) const 
+{    
+    scalar Vm = this->W()/rho;
+    return a_*n_*(n_+1)*pow(T/Tcrit_,0.5)*log(b_+Vm)
+        /(2*b_*pow(T,2))
+    -a_*n_*(n_+1)*pow(T/Tcrit_,0.5)*log(Vm)
+        /(2*b_*pow(T,2));
+}
+
+
+//Isobar expansion Coefficent beta = 1/v (dv/dt) at constant p
+//(molar values)
+inline scalar soaveRedlichKwong::isobarExpCoef(const scalar rho,const scalar T)  const
+{
+    return this->dvdT(rho, T)*rho/this->W();
+}
+
+
+//isothemal compressiblity kappa (not Thermal conductivity)
+//(molar values)
+inline scalar soaveRedlichKwong::isothermalCompressiblity(const scalar rho,const scalar T) const
+{
+    return this->isobarExpCoef(rho, T)/this->dpdT(rho, T);
+}
+
+
+//- Return density [kg/m^3]
+inline scalar soaveRedlichKwong::rho(
+    const scalar p,
+    const scalar T,
+    const scalar rho0
+) const
 {
  
    scalar molarVolumePrevIteration;
@@ -316,101 +340,104 @@ inline scalar soaveRedlichKwong::rho(const scalar p,const scalar T,const scalar
    int iter=0;
    int maxIter_=400;
    scalar tol_=1e-8;
-   int i;
+   scalar rho1=rhoMax_;
-   scalar rho1=rhoMax_, rho2=rhoMin_,rho3, f1,f2,f3;
+   scalar rho2=rhoMin_;
 
-molarVolume=this->W()/rho0;	 
+   molarVolume=this->W()/rho0;	 
 
    do
    {
-   molarVolumePrevIteration= molarVolume;
+       molarVolumePrevIteration= molarVolume;
    
+       label i=0;
+       do
+       {
+           molarVolume=molarVolumePrevIteration
+               -(
+                   (this->p((this->W()/molarVolumePrevIteration),T) - p)
+                   /(this->dpdv((this->W()/molarVolumePrevIteration),T))
+               )/pow(2,i);
 
-   i=0;
+           i++;
-      	do
+           if (i>8) 
-   	{
+           {
-  	molarVolume=molarVolumePrevIteration-((this->p((this->W()/molarVolumePrevIteration),T)-p)/(this->dpdv((this->W()/       	
+               //using bisection methode as backup,
-  	molarVolumePrevIteration),T)))/(pow(2,i)); 
+               //solution must be between rho=0.001 to rho=1500;
+               for(i=0; i<200; i++)
+               {
+                   scalar f1 = this->p(rho1,T) - p;
+                   scalar f2 = this->p(rho2,T) - p;
+                   scalar rho3 = (rho1 + rho2)/2;
+                   scalar f3 = this->p(rho3,T) - p;
 	
-	i++;
+                   if ((f2 < 0 && f3 > 0) || (f2 > 0 && f3 < 0))
-		if(i>8) 
+                   {
-		{
+                       rho1=rho3;
-		  //using bisection methode as backup, solution must be between rho=0.001 to rho=1500; 
+                   }  
-		  for(i=0;i<200;i++)
+                   else if ((f1 < 0 && f3 > 0)||(f1 > 0 && f3 < 0))
-		  {
+                   {
-                  	f1= (this->p(rho1,T)-p);
+                       rho2=rho3;
-                  	f2= (this->p(rho2,T)-p);
+                   }		 
-                  	rho3=(rho1+rho2)/2;
+                   else	
-		  	f3=(this->p(rho3,T)-p);
+                   {
+                       rho2=(rho2 + rho3)/2;
+                   }
 
-			if ((f2<0 && f3>0)||(f2>0 &&f3<0))	
+                   if(mag(f3) < p*tol_)
-		  	{
+                   {	  	  
-		 		rho1=rho3;
+                       molarVolume=this->W()/rho3;
-			}  
+                       molarVolumePrevIteration=this->W()/rho3;
- 		  	else if ((f1<0 && f3>0)||(f1>0 &&f3<0))
+                       break; 
-                  	{
+                   }
-		 		rho2=rho3;
+                   else
-			}
+                   {
-		 
+                       molarVolumePrevIteration=this->W()/rho3;
-			else	
+                   }
-			{
+               }		  
-                        	rho2=(rho2+rho3)/2;
+           }
-			}
+       }
-
+       while
-			if(mag(f3)<p*tol_)
+       (
-                  	{	  	  
+           mag(this->p((this->W()/molarVolume),T) - p)
-				molarVolume=this->W()/rho3;
+         > mag(this->p((this->W()/molarVolumePrevIteration),T) - p)
-                 	 	molarVolumePrevIteration=this->W()/rho3;
+       );
-				break; 
-		  	}
-			
-			else
-			{
-				molarVolumePrevIteration=this->W()/rho3;
-			}
-		    }		  
-  		}
-        }
-        while(mag(this->p((this->W()/molarVolume),T)-p)>mag(this->p((this->W()/molarVolumePrevIteration),T)-p));
-
-
-           if (iter++ > maxIter_)
-        {
-            FatalErrorIn
-            (
-               "inline scalar soaveRedlichKwong::rho(const scalar p,const scalar T,const scalar rho0) const "
-            )   << "Maximum number of iterations exceeded"
-                << abort(FatalError);
-        }
 
+       if (iter++ > maxIter_)
+       {
+           FatalErrorIn
+           (
+               "inline scalar redlichKwong::rho(const scalar p, const scalar T, const scalar rho0) const "
+           )   << "Maximum number of iterations exceeded"
+               << abort(FatalError);
+       }
    }
-   while(mag(molarVolumePrevIteration-molarVolume)>tol_*(this->W()/rho0));
+   while(mag(molarVolumePrevIteration-molarVolume) > tol_*(this->W()/rho0));
 
    return this->W()/molarVolume;
-
 }
 
 //- Return density [kg/m^3]on 
 inline scalar soaveRedlichKwong::rho(const scalar p,const scalar T) const
 {
 
-	scalar rho0=p/(this->R()*T); //using perfect gas equation as starting point
+    //using perfect gas equation as starting point
-	return rho(p,T,rho0);
+    return rho(p,T,p/(this->R()*T));
 }
 
-//- Return compressibility drho/dp [s^2/m^2]
+
+//- Return compressibility drho/dp at T=constant [s^2/m^2]
 inline scalar soaveRedlichKwong::psi(const scalar rho, const scalar T) const
 {
-  return -this->dvdp(rho,T)*pow(rho,2)/this->W();
+    return -this->dvdp(rho,T)*pow(rho,2)/this->W();
 }
 
+
 //- Return compression factor []
 inline scalar soaveRedlichKwong::Z( const scalar p, const scalar T,const scalar rho0) const
 {
-        return (p*this->rho(p,T,rho0))/(this->R()*T);
+    return (p*this->rho(p,T,rho0))/(this->R()*T);
 }
 
-
 // * * * * * * * * * * * * * * * Member Operators  * * * * * * * * * * * * * //
 
 inline void soaveRedlichKwong::operator+=(const soaveRedlichKwong& pg)

src/thermophysicalModels/specie/thermo/nasaHeatCapacityPolynomial/nasaHeatCapacityPolynomial.C

@@ -46,10 +46,16 @@ Foam::nasaHeatCapacityPolynomial<equationOfState>::nasaHeatCapacityPolynomial(Is
     a4_(readScalar(is)),
     a5_(readScalar(is)),
     a6_(readScalar(is)),
-    a7_(readScalar(is))
+    a7_(readScalar(is)),
+    //values for some need terms at std
+        e0_std(e0(this->Tstd)),
+	s0_std(s0(this->Tstd)),
+        integral_p_dv_std(this->integral_p_dv(this->rhostd(),this->Tstd)),
+        integral_dpdT_dv_std(this->integral_dpdT_dv(this->rhostd(),this->Tstd)),
+    // cp @ STD (needed to limit cp for stability
+        cp_std(this->cp_nonLimited(this->rhostd(),this->Tstd)) 	
 {
     is.check("nasaHeatCapacityPolynomial::nasaHeatCapacityPolynomial(Istream& is)"); 	
-    cp_std=this->cp_nonLimited(this->rhostd(),this->Tstd); 	// cp @ STD (needed to limit cp for stability
 }
 
 
@@ -64,7 +70,6 @@ Foam::Ostream& Foam::operator<<
 {
     os  << static_cast<const equationOfState&>(ct) << tab
         << ct.a1_ << tab<< ct.a2_ << tab << ct.a3_ << tab << ct.a4_ << tab << ct.a5_ << tab << ct.a6_ << tab << ct.a7_ ;
-    cout<<"nasa polynomal start"<<nl<<endl;
     os.check("Ostream& operator<<(Ostream& os, const nasaHeatCapacityPolynomial& ct)");
     return os;
 }

src/thermophysicalModels/specie/thermo/nasaHeatCapacityPolynomial/nasaHeatCapacityPolynomial.H

@@ -115,6 +115,10 @@ class nasaHeatCapacityPolynomial
 	scalar a5_;
 	scalar a6_;
         scalar a7_;
+	scalar e0_std;
+	scalar s0_std;
+	scalar integral_p_dv_std;
+	scalar integral_dpdT_dv_std;
 	scalar cp_std;
 
     // Private member functions
@@ -124,22 +128,20 @@ class nasaHeatCapacityPolynomial
         (
             const equationOfState& st,
             const scalar a1,
-			const scalar a2,
+	    const scalar a2,
-			const scalar a3,
+            const scalar a3,
-			const scalar a4,
+	    const scalar a4,
-		    const scalar a5,
+	    const scalar a5,
-		    const scalar a6,
+	    const scalar a6,
             const scalar a7
-			);
+	);
 
 
 public:
 
 //Variable
 
-
     // Constructors
-
         //- Construct from Istream
         nasaHeatCapacityPolynomial(Istream&);
 
@@ -152,49 +154,41 @@ public:
         //- Selector from Istream
         inline static autoPtr<nasaHeatCapacityPolynomial> New(Istream& is);
 
-
     // Member Functions
         
- 	//useful functions
+	//- perfect Gas Enthalpy [J/kmol]
-           //used to calculate the internal energy
+        inline scalar h0(const scalar T) const;
-	   //ideal Gas Enthalpy
-            inline scalar h0(const scalar T) const;
 	    
+        //- perfect Gas Entropy [J/(kmol K)]
+        inline scalar s0(const scalar T) const;
 
-           // used to calculate the entropy
+        //- perfect Gas internal Energy  [J/kmol]
-           //ideal Gas Entropy
+        inline scalar e0(const scalar T) const;
-            inline scalar s0(const scalar T) const;
 
+        //- perfect gas Heat capacity at constant pressure [J/(kmol K)]
+        inline scalar cv0(const scalar T) const;
 
-	   //used to calculate the internal energy
+        //- perfect gas Heat capacity at constant pressure [J/(kmol K)]
-           //ideal Gas Enthalpy 
+        inline scalar cp0(const scalar T) const;
-           inline scalar e0(const scalar T) const;
 
-        // Fundamental properties
+        //- Limited Heat capacity at constant pressure [J/(kmol K)]
+        inline scalar cp(const scalar rho, const scalar T) const;
 
-             //- ideal gas Heat capacity at constant pressure [J/(kmol K)]
+        //- non Limited Heat capacity at constant pressure [J/(kmol K)]
-            inline scalar cv0(const scalar T) const;
+	inline scalar cp_nonLimited(const scalar rho, const scalar T) const;
 
-             //- ideal gas Heat capacity at constant pressure [J/(kmol K)]
+        //- Heat capacity at constant pressure [J/(kmol K)]
-            inline scalar cp0(const scalar T) const;
+        inline scalar cv(const scalar rho, const scalar T) const;
 
-            //- Heat capacity at constant pressure [J/(kmol K)]
+        //- Enthalpy [J/kmol]
-            inline scalar cp(const scalar rho, const scalar T) const;
+        inline scalar h(const scalar rho, const scalar T) const;
 
-            //- Heat capacity at constant pressure [J/(kmol K)]
+        //- Entropy [J/(kmol K)]
-	    inline scalar cp_nonLimited(const scalar rho, const scalar T) const;
+        inline scalar s(const scalar rho,const scalar T) const;
 	    
-            //- Heat capacity at constant pressure [J/(kmol K)]
+        //- Internal Energy [J/kmol]
-            inline scalar cv(const scalar rho, const scalar T) const;
+        inline scalar e(const scalar rho, const scalar T) const;
 
-            //- Enthalpy [J/kmol]
-            inline scalar h(const scalar rho, const scalar T) const;
-
-            //- Entropy [J/(kmol K)]
-            inline scalar s(const scalar rho,const scalar T) const;
-	    
-            //- Internal Energy [J/kmol]
-            inline scalar e(const scalar rho, const scalar T) const;
 
 
     // Member operators

src/thermophysicalModels/specie/thermo/nasaHeatCapacityPolynomial/nasaHeatCapacityPolynomialI.H

@@ -33,7 +33,7 @@ Germany
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
 
-
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
 template<class equationOfState>
 inline Foam::nasaHeatCapacityPolynomial<equationOfState>::nasaHeatCapacityPolynomial
@@ -55,12 +55,16 @@ inline Foam::nasaHeatCapacityPolynomial<equationOfState>::nasaHeatCapacityPolyno
 	    a4_(a4),
 	    a5_(a5),
 	    a6_(a6),
-            a7_(a7)
+            a7_(a7),
-{
+            e0_std(e0(this->Tstd)),
-}
+	    s0_std(s0(this->Tstd)),
+            integral_p_dv_std(this->integral_p_dv(this->rhostd(),this->Tstd)),
+            integral_dpdT_dv_std(this->integral_dpdT_dv(this->rhostd(),this->Tstd)),
+            cp_std(this->cp_nonLimited(this->rhostd(),this->Tstd)) 
+{}
+
 
 
-// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
 template<class equationOfState>
 inline Foam::nasaHeatCapacityPolynomial<equationOfState>::nasaHeatCapacityPolynomial
@@ -76,9 +80,13 @@ inline Foam::nasaHeatCapacityPolynomial<equationOfState>::nasaHeatCapacityPolyno
 	    a4_(ct.a4_),
 	    a5_(ct.a5_),
 	    a6_(ct.a6_),
-            a7_(ct.a7_)
+            a7_(ct.a7_),
-{
+            e0_std(ct.e0_std),
-}
+            s0_std(ct.s0_std),
+            integral_p_dv_std(ct.integral_p_dv_std),
+            integral_dpdT_dv_std(ct.integral_dpdT_dv_std),
+            cp_std(ct.cp_std)
+{}
 
 
 template<class equationOfState>
@@ -102,12 +110,7 @@ Foam::nasaHeatCapacityPolynomial<equationOfState>::New(Istream& is)
     );
 }
 
-
+ // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
-
-
-
-
- // * * * * * * * * * * * * * * *useful functions* * * * * * * * * * * * * //
 
 //used to calculate the internal energy
 //perfect gas enthalpy 
@@ -117,7 +120,17 @@ inline Foam::scalar Foam::nasaHeatCapacityPolynomial<equationOfState>::h0
     const scalar T
 ) const
 {
-return this->RR*T*(-this->a1_*pow(T,-2)+this->a2_*log(T)/T+this->a3_+0.5*this->a4_*T+(this->a5_*pow(T,2))/3+(this->a6_*pow(T,3))/4+(this->a7_*pow(T,4))/5);
+    return 
+    this->RR*T*
+    (
+        -this->a1_*pow(T,-2)
+        +this->a2_*log(T)/T
+        +this->a3_
+        +0.5*this->a4_*T
+        +(this->a5_*pow(T,2))/3
+        +(this->a6_*pow(T,3))/4
+        +(this->a7_*pow(T,4))/5
+    );
 }
 
 
@@ -129,7 +142,7 @@ inline Foam::scalar Foam::nasaHeatCapacityPolynomial<equationOfState>::e0
     const scalar T
 ) const
 {
-return this->h0(T) - this->RR*T;
+    return this->h0(T) - this->RR*T;
 }
 	    
 // used to calculate the entropy
@@ -140,11 +153,19 @@ inline Foam::scalar Foam::nasaHeatCapacityPolynomial<equationOfState>::s0
     const scalar T
 ) const
 {
-return this->RR*( this->a1_*(-1)/(2*pow(T,2))-this->a2_/T+this->a3_*log(T)+this->a4_*T+(this->a5_*pow(T,2))/2+(this->a6_*pow(T,3))/3+(this->a7_*pow(T,4))/4);
+     return this->RR*
+     ( 
+          this->a1_*(-1)/(2*pow(T,2))
+          -this->a2_/T+this->a3_*log(T)
+          +this->a4_*T
+          +(this->a5_*pow(T,2))/2
+          +(this->a6_*pow(T,3))/3
+          +(this->a7_*pow(T,4))/4
+     );
 }
 
 
-// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
 
 //perfect gas cp
 template<class equationOfState>
@@ -153,8 +174,16 @@ inline Foam::scalar Foam::nasaHeatCapacityPolynomial<equationOfState>::cp0
     const scalar T
 ) const
 {
-
+    return this->RR*
-   return this->RR*(this->a1_*1/pow(T,2)+this->a2_*1/T+this->a3_+this->a4_*T+this->a5_*pow(T,2)+this->a6_*pow(T,3)+this->a7_*pow(T,4));
+    (
+        this->a1_*1/pow(T,2)
+        +this->a2_*1/T
+        +this->a3_
+        +this->a4_*T
+        +this->a5_*pow(T,2)
+        +this->a6_*pow(T,3)
+        +this->a7_*pow(T,4)
+    );
 }
 
 
@@ -166,7 +195,7 @@ inline Foam::scalar Foam::nasaHeatCapacityPolynomial<equationOfState>::cv0
 ) const
 {
 
-   return this->cp0(T)-this->RR;
+    return this->cp0(T)-this->RR;
 }
 
 
@@ -180,14 +209,25 @@ inline Foam::scalar Foam::nasaHeatCapacityPolynomial<equationOfState>::cp
      const scalar T 
 ) const
 {
-// Problem --> dpdv(rho,T)  is =0 at some points within the vapour dome. To increase stability, (dp/dv) has to be limited 
+    // Problem --> dpdv(rho,T)  is =0 at some points within the vapour dome. To increase stability, (dp/dv) has to be limited 
-// cp can be negative within the vapor dome. To avoid this nonphysical result, the absolute value is used.
+    // cp can be negative within the vapor dome. To avoid this nonphysical result, the absolute value is used.
-// within the vapourdome and at the critical point, cp increases to very high values --> infinity, 
+    // within the vapourdome and at the critical point, cp increases to very high values --> infinity, 
-// this would decrease the stability, so cp will be limited to 20 time the cp @ STD    
+    // this would decrease the stability, so cp will be limited to 20 time the cp @ STD    
- return  min(cp_std*20,fabs(this->cv(rho,T)-T*pow((this->dpdT(rho, T)),2)/min(this->dpdv(rho, T),-1))); 
+    return  
+    min
+    (
+        cp_std*20,
+        fabs
+        (
+            this->cv(rho,T)
+            -T*pow((this->dpdT(rho, T)),2)
+                /min(this->dpdv(rho, T),-1)
+        )
+    ); 
 }
 
-// this function is needed to get cp @ STD (without the limit imposed in the function above), which in turn is needed to limit the cp in the function above
+// this function is needed to get cp @ STD (without the limit imposed in the function above),
+// which in turn is needed to limit the cp in the function above
 template<class equationOfState>
 inline Foam::scalar Foam::nasaHeatCapacityPolynomial<equationOfState>::cp_nonLimited
 (
@@ -195,7 +235,7 @@ inline Foam::scalar Foam::nasaHeatCapacityPolynomial<equationOfState>::cp_nonLim
      const scalar T 
 ) const
 {
- return  fabs(this->cv(rho,T)-T*pow((this->dpdT(rho, T)),2)/min(this->dpdv(rho, T),-1));    
+     return  fabs(this->cv(rho,T)-T*pow((this->dpdT(rho, T)),2)/min(this->dpdv(rho, T),-1));    
 }
 
 
@@ -209,8 +249,7 @@ inline Foam::scalar Foam::nasaHeatCapacityPolynomial<equationOfState>::cv
     
 ) const
 {
-
+    return this->cv0(T)+T*this->integral_d2pdT2_dv(rho, T);
-   return this->cv0(T)+T*this->integral_d2pdT2_dv(rho, T);
 }
 
 
@@ -224,12 +263,12 @@ inline Foam::scalar Foam::nasaHeatCapacityPolynomial<equationOfState>::h
     const scalar T   
 ) const
 {
-  return this->e(rho,T)+this->p(rho,T)/rho*this->W()-this->Pstd/this->rhostd()*this->W();
+    return this->e(rho,T)+this->p(rho,T)/rho*this->W()-this->Pstd/this->rhostd()*this->W();
 }
 
 
 // function to calculate real gas internal energy
-// important assumption used is that the internal Energie is 0 at STD conditions.
+// important assumption used: internal Energie is 0 at STD conditions.
 // equation: du= cv0 dT +[T*dp/dT -p]dv
 template<class equationOfState>
 inline Foam::scalar Foam::nasaHeatCapacityPolynomial<equationOfState>::e
@@ -238,19 +277,20 @@ inline Foam::scalar Foam::nasaHeatCapacityPolynomial<equationOfState>::e
     const scalar T    
 ) const
 {
-return (	-this->Tstd*this->integral_dpdT_dv(this->rhostd(),this->Tstd)
+    return 
-			+this->integral_p_dv(this->rhostd(),this->Tstd)
+    (
-				+this->e0(T)-this->e0(this->Tstd)
+        -this->Tstd*integral_dpdT_dv_std
-					+T*this->integral_dpdT_dv(rho,T)
+	+integral_p_dv_std
-						-this->integral_p_dv(rho,T)
+	+this->e0(T)-e0_std
-	);
+	+T*this->integral_dpdT_dv(rho,T)
- 
+	-this->integral_p_dv(rho,T)
+    );
 }
 
 
 
 //function to calculate real gas entropy
-// important assumption used is that the Entropy is 0 at STD conditions.
+// important assumption used: the Entropy is 0 at STD conditions.
 // equation: ds= cv0/T * dT +  dp/dT *dv
 // --> integral cv0/T dT = s0(T) -s0(Tstd) - R*ln(T/Tstd)  --> due to s0(T)-s0(Tstd)=integral cp0/T dT
 template<class equationOfState>
@@ -261,9 +301,15 @@ inline Foam::scalar Foam::nasaHeatCapacityPolynomial<equationOfState>::s
     
 ) const
 {
-  return -this->integral_dpdT_dv(this->rhostd(),this->Tstd)+(this->s0(T)-this->s0(this->Tstd))-this->RR*log(T/this->Tstd)+ this->integral_dpdT_dv(rho,T);
+    return -integral_dpdT_dv_std
+    +(this->s0(T)-s0_std)
+    -this->RR*log(T/this->Tstd)
+    + this->integral_dpdT_dv(rho,T);
 }
 
+
+
+
 // * * * * * * * * * * * * * * * Member Operators  * * * * * * * * * * * * * //
 
 
@@ -283,11 +329,11 @@ inline void Foam::nasaHeatCapacityPolynomial<equationOfState>::operator+=
 
     a1_ = molr1*a1_ + molr2*ct.a1_;
     a2_ = molr1*a2_ + molr2*ct.a2_;
-	a3_ = molr1*a3_ + molr2*ct.a3_;
+    a3_ = molr1*a3_ + molr2*ct.a3_;
-	a4_ = molr1*a4_ + molr2*ct.a4_;
+    a4_ = molr1*a4_ + molr2*ct.a4_;
-	a5_ = molr1*a5_ + molr2*ct.a5_;
+    a5_ = molr1*a5_ + molr2*ct.a5_;
-	a6_ = molr1*a6_ + molr2*ct.a6_;
+    a6_ = molr1*a6_ + molr2*ct.a6_;
-	a7_ = molr1*a7_ + molr2*ct.a7_;
+    a7_ = molr1*a7_ + molr2*ct.a7_;
 }
 
 
@@ -306,11 +352,11 @@ inline void Foam::nasaHeatCapacityPolynomial<equationOfState>::operator-=
 
     a1_ = molr1*a1_ - molr2*ct.a1_;
     a2_ = molr1*a2_ - molr2*ct.a2_;
-	a3_ = molr1*a3_ - molr2*ct.a3_;
+    a3_ = molr1*a3_ - molr2*ct.a3_;
-	a4_ = molr1*a4_ - molr2*ct.a4_;
+    a4_ = molr1*a4_ - molr2*ct.a4_;
-	a5_ = molr1*a5_ - molr2*ct.a5_;
+    a5_ = molr1*a5_ - molr2*ct.a5_;
-	a6_ = molr1*a6_ - molr2*ct.a6_;
+    a6_ = molr1*a6_ - molr2*ct.a6_;
-	a7_ = molr1*a7_ - molr2*ct.a7_;
+    a7_ = molr1*a7_ - molr2*ct.a7_;
 }
 
 
@@ -336,15 +382,15 @@ inline Foam::nasaHeatCapacityPolynomial<equationOfState> Foam::operator+
       + ct2.nMoles()/eofs.nMoles()*ct2.a1_,
         ct1.nMoles()/eofs.nMoles()*ct1.a2_
       + ct2.nMoles()/eofs.nMoles()*ct2.a2_,
-	    ct1.nMoles()/eofs.nMoles()*ct1.a3_
+	ct1.nMoles()/eofs.nMoles()*ct1.a3_
       + ct2.nMoles()/eofs.nMoles()*ct2.a3_,
         ct1.nMoles()/eofs.nMoles()*ct1.a4_
       + ct2.nMoles()/eofs.nMoles()*ct2.a4_,
-	    ct1.nMoles()/eofs.nMoles()*ct1.a5_
+	ct1.nMoles()/eofs.nMoles()*ct1.a5_
       + ct2.nMoles()/eofs.nMoles()*ct2.a5_,
         ct1.nMoles()/eofs.nMoles()*ct1.a6_
       + ct2.nMoles()/eofs.nMoles()*ct2.a6_,    
-	    ct1.nMoles()/eofs.nMoles()*ct1.a7_
+	ct1.nMoles()/eofs.nMoles()*ct1.a7_
       + ct2.nMoles()/eofs.nMoles()*ct2.a7_
     );
 }
@@ -370,15 +416,15 @@ inline Foam::nasaHeatCapacityPolynomial<equationOfState> Foam::operator-
       - ct2.nMoles()/eofs.nMoles()*ct2.a1_,
         ct1.nMoles()/eofs.nMoles()*ct1.a2_
       - ct2.nMoles()/eofs.nMoles()*ct2.a2_,
-	    ct1.nMoles()/eofs.nMoles()*ct1.a3_
+	ct1.nMoles()/eofs.nMoles()*ct1.a3_
       - ct2.nMoles()/eofs.nMoles()*ct2.a3_,
         ct1.nMoles()/eofs.nMoles()*ct1.a4_
       - ct2.nMoles()/eofs.nMoles()*ct2.a4_,
-	    ct1.nMoles()/eofs.nMoles()*ct1.a5_
+	ct1.nMoles()/eofs.nMoles()*ct1.a5_
       - ct2.nMoles()/eofs.nMoles()*ct2.a5_,
         ct1.nMoles()/eofs.nMoles()*ct1.a6_
       - ct2.nMoles()/eofs.nMoles()*ct2.a6_,    
-	    ct1.nMoles()/eofs.nMoles()*ct1.a7_
+	ct1.nMoles()/eofs.nMoles()*ct1.a7_
       - ct2.nMoles()/eofs.nMoles()*ct2.a7_
     );
 }
@@ -396,11 +442,11 @@ inline Foam::nasaHeatCapacityPolynomial<equationOfState> Foam::operator*
         s*static_cast<const equationOfState&>(ct),
         ct.a1_,
         ct.a2_,
-		ct.a3_,
+	ct.a3_,
         ct.a4_,
-		ct.a5_,
+	ct.a5_,
         ct.a6_,
-		ct.a7_
+	ct.a7_
     );
 }
 

src/thermophysicalModels/specie/thermo/realGasSpecieThermo/realGasSpecieThermo.H

@@ -205,8 +205,6 @@ public:
                 //- Enthalpy [J/kg]
                 inline scalar H(const scalar rho, const scalar T) const;
 
-
-
                 //- Sensible enthalpy [J/kg]
                // inline scalar Hs(const scalar T) const;
 
@@ -225,9 +223,13 @@ public:
                 //- Helmholtz free energy [J/kg]
                 inline scalar A(const scalar rho, const scalar T) const;
 
+                //Other variables
 
+	        //- Return compressibility drho/dp at h=constant [s^2/m^2] 
+                inline scalar psiH(const scalar rho, const scalar T) const;
 
-
+                //- Return compressibility drho/dp at e=constant [s^2/m^2]     
+                inline scalar psiE(const scalar rho, const scalar T) const;
 
 
         // Energy->temperature  inversion functions

src/thermophysicalModels/specie/thermo/realGasSpecieThermo/realGasSpecieThermoI.H

@@ -47,7 +47,7 @@ inline Foam::realGasSpecieThermo<thermo>::realGasSpecieThermo
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
 
-// using a stabilizing newton solver 
+// using two one dimensional newton solvers in a row
 template<class thermo>
 inline void Foam::realGasSpecieThermo<thermo>::T
 (
@@ -59,47 +59,59 @@ inline void Foam::realGasSpecieThermo<thermo>::T
     scalar (realGasSpecieThermo<thermo>::*dFdT)(const scalar,const scalar) const
 ) const
 {
-
+    scalar Test ;
-   
-    scalar Test = T0;
     scalar Tnew = T0;
+    scalar rhoOld;
     scalar rho=rho0;
     scalar Ttol = T0*tol_;
-    int    iter = 0;
+    scalar rhotol=rho0*tol_;
-int i;
+    int iter = 0;
+    int i;
 
     do
     {
         Test = Tnew;   
-        
+	rhoOld=rho;
-	rho=this->rho(p,Test,rho);
+	rho=this->rho(p,Test,rhoOld);
-
+        i=0;
-
-   i=0;
-;
       	do
    	{   
-        Tnew = Test - ((this->*F)(rho,Test) - f)/(this->*dFdT)(rho,Test)/(pow(2,i)); 
+            // using a stabilizing newton solver 
-	i++;				// if the solve is diverging, the time step is reduced until the solver converges	
+            // if the solve is diverging, the time step is reduced until the solver converges
-         }while((i<20)&&((mag((this->*F)(rho,Tnew) - f) > mag((this->*F)(rho,Test) - f))));	
+            Tnew = Test - ((this->*F)(rho,Test) - f)/(this->*dFdT)(rho,Test)/(pow(2,i)); 
+	    i++;					
+        }while
+        (
+            (i<20)
+            &&
+            ((
+                mag((this->*F)(rho,Tnew) - f) 
+                > 
+                mag((this->*F)(rho,Test) - f)
+            ))
+        );	
 
         if (iter++ > maxIter_)
         {
             FatalErrorIn
             (
-                "realGasSpecieThermo<thermo>::T(scalar f, scalar T0, "
+                "realGasSpecieThermo<thermo>::T(scalar f, scalar T0, scalar p, scalar rho0, "
                 "scalar (realGasSpecieThermo<thermo>::*F)(const scalar) const, "
                 "scalar (realGasSpecieThermo<thermo>::*dFdT)(const scalar) const"
                 ") const"
             )   << "Maximum number of iterations exceeded"
                 << abort(FatalError);
         }
+    } while
+    // both fields must converge
+    ( 
+        (mag(mag(Tnew) - mag(Test)) > Ttol)
+        ||
+        (mag(mag(rho) - mag(rhoOld)) > rhotol)
+    );		
 
-
+    rho0=rho;
-    } while (mag(mag(Tnew) - mag(Test)) > Ttol);		
+    T0=Tnew;
-
-  rho0=rho;
-  T0=Tnew;
 }
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
@@ -184,8 +196,52 @@ inline Foam::scalar Foam::realGasSpecieThermo<thermo>::A(const scalar rho,  cons
     return this->a(rho, T)/this->W();
 }
 
+//- Return compressibility drho/dp at h=constant [s^2/m^2]
+//- using Bridgeman's Table
+template<class thermo>
+inline Foam::scalar Foam::realGasSpecieThermo<thermo>::psiH
+(
+    const scalar rho,
+    const scalar T
     
+) const
+{
+    scalar cp=this->cp(rho,T);
+    scalar beta=this->isobarExpCoef(rho,T);
+    return  
+    -(
+        T*pow(beta,2)/cp
+        -beta/cp
+        -this->isothermalCompressiblity(rho,T)*rho/this->W()
+    )*this->W();
+}
 
+//- Return compressibility drho/dp at e=constant [s^2/m^2]
+//- using Bridgeman's Table
+template<class thermo>
+inline Foam::scalar Foam::realGasSpecieThermo<thermo>::psiE
+(
+    const scalar rho,
+    const scalar T
+    
+) const
+{
+    scalar Vm = this->W()/rho;
+    scalar cp=this->cp(rho,T);
+    scalar beta=this->isobarExpCoef(rho,T);
+    return  
+    -(
+        (
+	    T*pow(beta,2)*Vm
+	    -this->isothermalCompressiblity(rho,T)*cp
+	)
+	/
+	(
+            cp*Vm
+            -beta*this->p(rho,T)*pow(Vm,2)
+	)
+    )*this->W();
+}
 
 template<class thermo>
 inline void Foam::realGasSpecieThermo<thermo>::TH

src/thermophysicalModels/specie/transport/const/constTransport.H

@@ -137,6 +137,14 @@ public:
         //- Thermal diffusivity for enthalpy [kg/ms]
         inline scalar alpha(const scalar T) const;
 
+
+        //- Thermal conductivity [W/mK]
+        inline scalar kappa(const scalar rho,const scalar T) const;
+
+        //- Thermal diffusivity for enthalpy [kg/ms]
+        inline scalar alpha(const scalar rho,const scalar T) const;
+
+
         // Species diffusivity
         //inline scalar D(const scalar T) const;
 

src/thermophysicalModels/specie/transport/const/constTransportI.H

@@ -116,7 +116,26 @@ inline scalar constTransport<thermo>::alpha(const scalar T) const
     return Cp_*mu(T)*rPr/CpBar;
 }
 
+// Thermal conductivity [W/mK]
+template<class thermo>
+inline scalar constTransport<thermo>::kappa(const scalar rho,const scalar T) const
+{
+    return this->Cp(rho,T)*mu(T)*rPr;
+}
 
+
+// Thermal diffusivity for enthalpy [kg/ms]
+template<class thermo>
+inline scalar constTransport<thermo>::alpha(const scalar rho,const scalar T) const
+{
+    scalar Cp_ = this->Cp(rho,T);
+
+    scalar deltaT = T - specie::Tstd;
+    scalar CpBar =
+        (deltaT*(this->H(rho,T) - this->H(this->rhostd(),specie::Tstd)) + Cp_)/(sqr(deltaT) + 1);
+
+    return Cp_*mu(T)*rPr/CpBar;
+}
 // * * * * * * * * * * * * * * * Member Operators  * * * * * * * * * * * * * //
 
 template<class thermo>

src/thermophysicalModels/specie/transport/realGasConstTransport/realGasConstTransport.C

@@ -1,79 +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
-
-    Modified version of constTransport, modified to be used with realGases
-
-    Constant properties Transport package.  Templated ito a given
-    thermodynamics package (needed for thermal conductivity).
-
-Modified by
-Christian Lucas
-Institut für Thermodynamik
-Technische Universität Braunschweig 
-Germany
-
-\*---------------------------------------------------------------------------*/
-
-#include "realGasConstTransport.H"
-#include "IOstreams.H"
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-namespace Foam
-{
-
-// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
-
-template<class thermo>
-realGasConstTransport<thermo>::realGasConstTransport(Istream& is)
-:
-    thermo(is),
-    Mu(readScalar(is)),
-    rPr(1.0/readScalar(is))
-{
-    is.check("realGasConstTransport::realGasConstTransport(Istream& is)");
-}
-
-
-// * * * * * * * * * * * * * * * IOstream Operators  * * * * * * * * * * * * //
-
-template<class thermo>
-Ostream& operator<<(Ostream& os, const realGasConstTransport<thermo>& ct)
-{
-    operator<<(os, static_cast<const thermo&>(ct));
-    os << tab << ct.Mu << tab << 1.0/ct.rPr;
-
-    os.check("Ostream& operator<<(Ostream& os, const realGasConstTransport& ct)");
-
-    return os;
-}
-
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-} // End namespace Foam
-
-// ************************************************************************* //

src/thermophysicalModels/specie/transport/realGasConstTransport/realGasConstTransport.H

@@ -1,204 +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
-
-Class
-    Foam::realGasConstTransport
-
-Description
-    Constant properties Transport package.
-    Templated into a given thermodynamics package (needed for thermal
-    conductivity).
-
-SourceFiles
-    realGasConstTransportI.H
-    realGasConstTransport.C
-
-Modified by
-Christian Lucas
-Institut für Thermodynamik
-Technische Universität Braunschweig 
-Germany
-\*---------------------------------------------------------------------------*/
-
-#ifndef realGasConstTransport_H
-#define realGasConstTransport_H
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-namespace Foam
-{
-
-// Forward declaration of friend functions and operators
-
-template<class thermo> class realGasConstTransport;
-
-template<class thermo>
-inline realGasConstTransport<thermo> operator+
-(
-    const realGasConstTransport<thermo>&,
-    const realGasConstTransport<thermo>&
-);
-
-template<class thermo>
-inline realGasConstTransport<thermo> operator-
-(
-    const realGasConstTransport<thermo>&,
-    const realGasConstTransport<thermo>&
-);
-
-template<class thermo>
-inline realGasConstTransport<thermo> operator*
-(
-    const scalar,
-    const realGasConstTransport<thermo>&
-);
-
-template<class thermo>
-inline realGasConstTransport<thermo> operator==
-(
-    const realGasConstTransport<thermo>&,
-    const realGasConstTransport<thermo>&
-);
-
-template<class thermo>
-Ostream& operator<<
-(
-    Ostream&,
-    const realGasConstTransport<thermo>&
-);
-
-
-/*---------------------------------------------------------------------------*\
-                           Class realGasConstTransport Declaration
-\*---------------------------------------------------------------------------*/
-
-template<class thermo>
-class realGasConstTransport
-:
-    public thermo
-{
-    // Private data
-
-        //- Constant viscosity and reciprocal Prandtl Number.
-        scalar Mu, rPr;
-
-
-    // Private member functions
-
-        //- Construct from components
-        inline realGasConstTransport
-        (
-            const thermo& t,
-            const scalar nu,
-            const scalar Pr
-        );
-
-
-public:
-
-    // Constructors
-
-        //- Construct as named copy
-        inline realGasConstTransport(const word&, const realGasConstTransport&);
-
-        //- Construct from Istream
-        realGasConstTransport(Istream&);
-
-
-    // Member functions
-
-        //- Dynamic viscosity [kg/ms]
-        inline scalar mu(const scalar T) const;
-
-        //- Thermal conductivity [W/mK]
-        inline scalar kappa(const scalar rho,const scalar T) const;
-
-        //- Thermal diffusivity for enthalpy [kg/ms]
-        inline scalar alpha(const scalar rho,const scalar T) const;
-
-
-
-
-    // Member operators
-
-        inline realGasConstTransport& operator=
-        (
-            const realGasConstTransport&
-        );
-
-
-    // Friend operators
-
-        friend realGasConstTransport operator+ <thermo>
-        (
-            const realGasConstTransport&,
-            const realGasConstTransport&
-        );
-
-        friend realGasConstTransport operator- <thermo>
-        (
-            const realGasConstTransport&,
-            const realGasConstTransport&
-        );
-
-        friend realGasConstTransport operator* <thermo>
-        (
-            const scalar,
-            const realGasConstTransport&
-        );
-
-        friend realGasConstTransport operator== <thermo>
-        (
-            const realGasConstTransport&,
-            const realGasConstTransport&
-        );
-
-
-    // Ostream Operator
-
-        friend Ostream& operator<< <thermo>
-        (
-            Ostream&,
-            const realGasConstTransport&
-        );
-};
-
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-} // End namespace Foam
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-#include "realGasConstTransportI.H"
-
-#ifdef NoRepository
-#   include "realGasConstTransport.C"
-#endif
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-#endif
-
-// ************************************************************************* //

src/thermophysicalModels/specie/transport/realGasConstTransport/realGasConstTransportI.H

@@ -1,199 +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
-
-Modified by
-Christian Lucas
-Institut für Thermodynamik
-Technische Universität Braunschweig 
-Germany
-
-\*---------------------------------------------------------------------------*/
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-namespace Foam
-{
-
-// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
-
-// Construct from components
-template<class thermo>
-inline realGasConstTransport<thermo>::realGasConstTransport
-(
-    const thermo& t,
-    const scalar mu,
-    const scalar Pr
-)
-:
-    thermo(t),
-    Mu(mu),
-    rPr(1.0/Pr)
-{}
-
-
-// Construct as named copy
-template<class thermo>
-inline realGasConstTransport<thermo>::realGasConstTransport
-(
-    const word& name,
-    const realGasConstTransport& ct
-)
-:
-    thermo(name, ct),
-    Mu(ct.Mu),
-    rPr(ct.rPr)
-{}
-
-
-// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
-
-// Dynamic viscosity [kg/ms]
-template<class thermo>
-inline scalar realGasConstTransport<thermo>::mu(const scalar T) const
-{
-    return Mu;
-}
-
-
-// Thermal conductivity [W/mK]
-template<class thermo>
-inline scalar realGasConstTransport<thermo>::kappa(const scalar rho,const scalar T) const
-{
-    return this->Cp(rho,T)*mu(T)*rPr;
-}
-
-
-// Thermal diffusivity for enthalpy [kg/ms]
-template<class thermo>
-inline scalar realGasConstTransport<thermo>::alpha(const scalar rho,const scalar T) const
-{
-    scalar Cp_ = this->Cp(rho,T);
-
-    scalar deltaT = T - specie::Tstd;
-    scalar CpBar =
-        (deltaT*(this->H(rho,T) - this->H(this->rhostd(),specie::Tstd)) + Cp_)/(sqr(deltaT) + 1);
-
-    return Cp_*mu(T)*rPr/CpBar;
-}
-
-
-// * * * * * * * * * * * * * * * Member Operators  * * * * * * * * * * * * * //
-
-template<class thermo>
-inline realGasConstTransport<thermo>& realGasConstTransport<thermo>::operator=
-(
-    const realGasConstTransport<thermo>& ct
-)
-{
-    thermo::operator=(ct);
-
-    Mu = ct.Mu;
-    rPr = ct.rPr;
-
-    return *this;
-}
-
-
-// * * * * * * * * * * * * * * * Friend Operators  * * * * * * * * * * * * * //
-
-template<class thermo>
-inline realGasConstTransport<thermo> operator+
-(
-    const realGasConstTransport<thermo>& ct1,
-    const realGasConstTransport<thermo>& ct2
-)
-{
-    thermo t
-    (
-        static_cast<const thermo&>(ct1) + static_cast<const thermo&>(ct2)
-    );
-
-    scalar molr1 = ct1.nMoles()/t.nMoles();
-    scalar molr2 = ct2.nMoles()/t.nMoles();
-
-    return realGasConstTransport<thermo>
-    (
-        t,
-        molr1*ct1.Mu + molr2*ct2.Mu,
-        molr1*ct1.rPr + molr2*ct2.rPr
-    );
-}
-
-
-template<class thermo>
-inline realGasConstTransport<thermo> operator-
-(
-    const realGasConstTransport<thermo>& ct1,
-    const realGasConstTransport<thermo>& ct2
-)
-{
-    thermo t
-    (
-        static_cast<const thermo&>(ct1) - static_cast<const thermo&>(ct2)
-    );
-
-    scalar molr1 = ct1.nMoles()/t.nMoles();
-    scalar molr2 = ct2.nMoles()/t.nMoles();
-
-    return realGasConstTransport<thermo>
-    (
-        t,
-        molr1*ct1.Mu - molr2*ct2.Mu,
-        molr1*ct1.rPr - molr2*ct2.rPr
-    );
-}
-
-
-template<class thermo>
-inline realGasConstTransport<thermo> operator*
-(
-    const scalar s,
-    const realGasConstTransport<thermo>& ct
-)
-{
-    return realGasConstTransport<thermo>
-    (
-        s*static_cast<const thermo&>(ct),
-        ct.Mu,
-        ct.rPr
-    );
-}
-
-
-template<class thermo>
-inline realGasConstTransport<thermo> operator==
-(
-    const realGasConstTransport<thermo>& ct1,
-    const realGasConstTransport<thermo>& ct2
-)
-{
-    return ct2 - ct1;
-}
-
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-} // End namespace Foam
-
-// ************************************************************************* //

src/thermophysicalModels/specie/transport/realGasSutherlandTransport/realGasSutherlandTransport.C

@@ -1,74 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     |
-    \\  /    A nd           | 
-     \\/     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
-
-Modified by
-Christian Lucas
-Institut für Thermodynamik
-Technische Universität Braunschweig 
-Germany
-
-\*---------------------------------------------------------------------------*/
-
-#include "realGasSutherlandTransport.H"
-#include "IOstreams.H"
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-namespace Foam
-{
-
-// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
-
-template<class thermo>
-realGasSutherlandTransport<thermo>::realGasSutherlandTransport(Istream& is)
-:
-    thermo(is),
-    As(readScalar(is)),
-    Ts(readScalar(is))
-{
-    is.check("realGasSutherlandTransport<thermo>::realGasSutherlandTransport(Istream&)");
-}
-
-
-// * * * * * * * * * * * * * * * IOstream Operators  * * * * * * * * * * * * //
-
-template<class thermo>
-Ostream& operator<<(Ostream& os, const realGasSutherlandTransport<thermo>& st)
-{
-    os << static_cast<const thermo&>(st) << tab << st.As << tab << st.Ts;
-    cout<<"sutherland start"<<nl<<endl;
-    os.check
-    (
-        "Ostream& operator<<(Ostream&, const realGasSutherlandTransport<thermo>&)"
-    );
-
-    return os;
-}
-
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-} // End namespace Foam
-
-// ************************************************************************* //

src/thermophysicalModels/specie/transport/realGasSutherlandTransport/realGasSutherlandTransport.H

@@ -1,237 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     |
-    \\  /    A nd           | 
-     \\/     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::realGasSutherlandTransport
-
-Description
-
-    Copy of the sutherlandTransport class, modified to be used with real gas properties
-
-    Transport package using Sutherland's formula.
-
-    Templated into a given thermodynamics package (needed for thermal
-    conductivity).
-
-    Dynamic viscosity [kg/m.s]
-    @f[
-        \mu = A_s \frac{\sqrt{T}}{1 + T_s / T}
-    @f]
-
-SourceFiles
-    realGasSutherlandTransportI.H
-    realGasSutherlandTransport.C
-
-Modified by
-Christian Lucas
-Institut für Thermodynamik
-Technische Universität Braunschweig 
-Germany
-
-\*---------------------------------------------------------------------------*/
-
-#ifndef realGasSutherlandTransport_H
-#define realGasSutherlandTransport_H
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-namespace Foam
-{
-
-// Forward declaration of friend functions and operators
-
-template<class thermo> class realGasSutherlandTransport;
-
-template<class thermo>
-inline realGasSutherlandTransport<thermo> operator+
-(
-    const realGasSutherlandTransport<thermo>&,
-    const realGasSutherlandTransport<thermo>&
-);
-
-template<class thermo>
-inline realGasSutherlandTransport<thermo> operator-
-(
-    const realGasSutherlandTransport<thermo>&,
-    const realGasSutherlandTransport<thermo>&
-);
-
-template<class thermo>
-inline realGasSutherlandTransport<thermo> operator*
-(
-    const scalar,
-    const realGasSutherlandTransport<thermo>&
-);
-
-template<class thermo>
-inline realGasSutherlandTransport<thermo> operator==
-(
-    const realGasSutherlandTransport<thermo>&,
-    const realGasSutherlandTransport<thermo>&
-);
-
-template<class thermo>
-Ostream& operator<<
-(
-    Ostream&,
-    const realGasSutherlandTransport<thermo>&
-);
-
-
-/*---------------------------------------------------------------------------*\
-                           Class realGasSutherlandTransport Declaration
-\*---------------------------------------------------------------------------*/
-
-template<class thermo>
-class realGasSutherlandTransport
-:
-    public thermo
-{
-    // Private data
-
-        // Sutherland's coefficients
-        scalar As, Ts;
-
-
-    // Private member functions
-
-        //- Calculate the Sutherland coefficients
-        //  given two viscosities and temperatures
-        inline void calcCoeffs
-        (
-            const scalar mu1, const scalar T1,
-            const scalar mu2, const scalar T2
-        );
-
-
-public:
-
-    // Constructors
-
-        //- Construct from components
-        inline realGasSutherlandTransport
-        (
-            const thermo& t,
-            const scalar as,
-            const scalar ts
-        );
-
-        //- Construct from two viscosities
-        inline realGasSutherlandTransport
-        (
-            const thermo& t,
-            const scalar mu1, const scalar T1,
-            const scalar mu2, const scalar T2
-        );
-
-        //- Construct as named copy
-        inline realGasSutherlandTransport(const word&, const realGasSutherlandTransport&);
-
-        //- Construct from Istream
-        realGasSutherlandTransport(Istream&);
-
-        //- Construct and return a clone
-        inline autoPtr<realGasSutherlandTransport> clone() const;
-
-        // Selector from Istream
-        inline static autoPtr<realGasSutherlandTransport> New(Istream& is);
-
-
-    // Member functions 
-
-
-        //- Dynamic viscosity [kg/ms]
-        inline scalar mu(const scalar T) const;
-
-
-        //- Thermal conductivity [W/mK]
-        inline scalar kappa(const scalar rho,const scalar T) const;
-
-        //- Thermal diffusivity for enthalpy [kg/ms]
-        inline scalar alpha(const scalar rho,const scalar T) const;
-
-
-
-    // Member operators
-
-        inline realGasSutherlandTransport& operator=
-        (
-            const realGasSutherlandTransport&
-        );
-
-
-    // Friend operators
-
-        friend realGasSutherlandTransport operator+ <thermo>
-        (
-            const realGasSutherlandTransport&,
-            const realGasSutherlandTransport&
-        );
-
-        friend realGasSutherlandTransport operator- <thermo>
-        (
-            const realGasSutherlandTransport&,
-            const realGasSutherlandTransport&
-        );
-
-        friend realGasSutherlandTransport operator* <thermo>
-        (
-            const scalar,
-            const realGasSutherlandTransport&
-        );
-
-        friend realGasSutherlandTransport operator== <thermo>
-        (
-            const realGasSutherlandTransport&,
-            const realGasSutherlandTransport&
-        );
-
-
-    // Ostream Operator
-
-        friend Ostream& operator<< <thermo>
-        (
-            Ostream&,
-            const realGasSutherlandTransport&
-        );
-};
-
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-} // End namespace Foam
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-#include "realGasSutherlandTransportI.H"
-
-#ifdef NoRepository
-#   include "realGasSutherlandTransport.C"
-#endif
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-#endif
-
-// ************************************************************************* //

src/thermophysicalModels/specie/transport/realGasSutherlandTransport/realGasSutherlandTransportI.H

@@ -1,262 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     |
-    \\  /    A nd           | 
-     \\/     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
-
-Modified by
-Christian Lucas
-Institut für Thermodynamik
-Technische Universität Braunschweig 
-Germany
-
-
-\*---------------------------------------------------------------------------*/
-
-#include "specie.H"
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-namespace Foam
-{
-
-// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
-
-template<class thermo>
-inline void realGasSutherlandTransport<thermo>::calcCoeffs
-(
-    const scalar mu1, const scalar T1,
-    const scalar mu2, const scalar T2
-)
-{
-    scalar rootT1 = sqrt(T1);
-    scalar mu1rootT2 = mu1*sqrt(T2);
-    scalar mu2rootT1 = mu2*rootT1;
-
-    Ts = (mu2rootT1 - mu1rootT2)/(mu1rootT2/T1 - mu2rootT1/T2);
-
-    As = mu1*(1.0 + Ts/T1)/rootT1;
-}
-
-
-// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
-
-// Construct from components
-template<class thermo>
-inline realGasSutherlandTransport<thermo>::realGasSutherlandTransport
-(
-    const thermo& t,
-    const scalar as,
-    const scalar ts
-)
-:
-    thermo(t),
-    As(as),
-    Ts(ts)
-{}
-
-
-// Construct from components
-template<class thermo>
-inline realGasSutherlandTransport<thermo>::realGasSutherlandTransport
-(
-    const thermo& t,
-    const scalar mu1, const scalar T1,
-    const scalar mu2, const scalar T2
-)
-:
-    thermo(t)
-{
-    calcCoeffs(mu1, T1, mu2, T2);
-}
-
-
-//- Construct as named copy
-template<class thermo>
-inline realGasSutherlandTransport<thermo>::realGasSutherlandTransport
-(
-    const word& name,
-    const realGasSutherlandTransport& st
-)
-:
-    thermo(name, st),
-    As(st.As),
-    Ts(st.Ts)
-{}
-
-
-// Construct and return a clone
-template<class thermo>
-inline autoPtr<realGasSutherlandTransport<thermo> > realGasSutherlandTransport<thermo>::clone
-() const
-{
-    return autoPtr<realGasSutherlandTransport<thermo> >
-    (
-        new realGasSutherlandTransport<thermo>(*this)
-    );
-}
-
-
-// Selector from Istream
-template<class thermo>
-inline autoPtr<realGasSutherlandTransport<thermo> > realGasSutherlandTransport<thermo>::New
-(
-    Istream& is
-)
-{
-    return autoPtr<realGasSutherlandTransport<thermo> >
-    (
-        new realGasSutherlandTransport<thermo>(is)
-    );
-}
-
-
-// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
-
-// Dynamic viscosity [kg/ms]
-template<class thermo>
-inline scalar realGasSutherlandTransport<thermo>::mu(const scalar T) const
-{
-    return As*::sqrt(T)/(1.0 + Ts/T);
-}
-
-
-// Thermal conductivity [W/mK]
-template<class thermo>
-inline scalar realGasSutherlandTransport<thermo>::kappa(const scalar rho, const scalar T) const
-{
-    scalar Cv_ = this->Cv(rho,T);
-    return mu(T)*Cv_*(1.32 + 1.77*this->R()/Cv_);
-}
-
-
-// Thermal diffusivity for enthalpy [kg/ms]
-template<class thermo>
-inline scalar realGasSutherlandTransport<thermo>::alpha(const scalar rho, const scalar T) const
-{
-    scalar Cv_ = this->Cv(rho,T);
-    scalar Cp_ = this->Cp(rho,T);
-
-    scalar deltaT = T - specie::Tstd;
-    scalar CpBar =
-        (deltaT*(this->H(rho, T) - this->H(this->rhostd(),specie::Tstd)) + Cp_)/(sqr(deltaT) + 1);
-
-    return mu(T)*Cv_*(1.32 + 1.77*this->R()/Cv_)/CpBar;
-}
-// * * * * * * * * * * * * * * * Member Operators  * * * * * * * * * * * * * //
-
-template<class thermo>
-inline realGasSutherlandTransport<thermo>& realGasSutherlandTransport<thermo>::operator=
-(
-    const realGasSutherlandTransport<thermo>& st
-)
-{
-    thermo::operator=(st);
-
-    As = st.As;
-    Ts = st.Ts;
-
-    return *this;
-}
-
-
-// * * * * * * * * * * * * * * * Friend Operators  * * * * * * * * * * * * * //
-
-template<class thermo>
-inline realGasSutherlandTransport<thermo> operator+
-(
-    const realGasSutherlandTransport<thermo>& st1,
-    const realGasSutherlandTransport<thermo>& st2
-)
-{
-    thermo t
-    (
-        static_cast<const thermo&>(st1) + static_cast<const thermo&>(st2)
-    );
-
-    scalar molr1 = st1.nMoles()/t.nMoles();
-    scalar molr2 = st2.nMoles()/t.nMoles();
-
-    return realGasSutherlandTransport<thermo>
-    (
-        t,
-        molr1*st1.As + molr2*st2.As,
-        molr1*st1.Ts + molr2*st2.Ts
-    );
-}
-
-
-template<class thermo>
-inline realGasSutherlandTransport<thermo> operator-
-(
-    const realGasSutherlandTransport<thermo>& st1,
-    const realGasSutherlandTransport<thermo>& st2
-)
-{
-    thermo t
-    (
-        static_cast<const thermo&>(st1) - static_cast<const thermo&>(st2)
-    );
-
-    scalar molr1 = st1.nMoles()/t.nMoles();
-    scalar molr2 = st2.nMoles()/t.nMoles();
-
-    return realGasSutherlandTransport<thermo>
-    (
-        t,
-        molr1*st1.As - molr2*st2.As,
-        molr1*st1.Ts - molr2*st2.Ts
-    );
-}
-
-
-template<class thermo>
-inline realGasSutherlandTransport<thermo> operator*
-(
-    const scalar s,
-    const realGasSutherlandTransport<thermo>& st
-)
-{
-    return realGasSutherlandTransport<thermo>
-    (
-        s*static_cast<const thermo&>(st),
-        st.As,
-        st.Ts
-    );
-}
-
-
-template<class thermo>
-inline realGasSutherlandTransport<thermo> operator==
-(
-    const realGasSutherlandTransport<thermo>& st1,
-    const realGasSutherlandTransport<thermo>& st2
-)
-{
-    return st2 - st1;
-}
-
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-} // End namespace Foam
-
-// ************************************************************************* //

src/thermophysicalModels/specie/transport/sutherland/sutherlandTransport.H

@@ -159,6 +159,14 @@ public:
         //- Thermal diffusivity for enthalpy [kg/ms]
         inline scalar alpha(const scalar T) const;
 
+
+        //- Thermal conductivity [W/mK]
+        inline scalar kappa(const scalar rho,const scalar T) const;
+
+        //- Thermal diffusivity for enthalpy [kg/ms]
+        inline scalar alpha(const scalar rho,const scalar T) const;
+
+
         // Species diffusivity
         //inline scalar D(const scalar T) const;
 

src/thermophysicalModels/specie/transport/sutherland/sutherlandTransportI.H

@@ -155,7 +155,28 @@ inline scalar sutherlandTransport<thermo>::alpha(const scalar T) const
     return mu(T)*Cv_*(1.32 + 1.77*this->R()/Cv_)/CpBar;
 }
 
+// Thermal conductivity [W/mK]
+template<class thermo>
+inline scalar sutherlandTransport<thermo>::kappa(const scalar rho, const scalar T) const
+{
+    scalar Cv_ = this->Cv(rho,T);
+    return mu(T)*Cv_*(1.32 + 1.77*this->R()/Cv_);
+}
 
+
+// Thermal diffusivity for enthalpy [kg/ms]
+template<class thermo>
+inline scalar sutherlandTransport<thermo>::alpha(const scalar rho, const scalar T) const
+{
+    scalar Cv_ = this->Cv(rho,T);
+    scalar Cp_ = this->Cp(rho,T);
+
+    scalar deltaT = T - specie::Tstd;
+    scalar CpBar =
+        (deltaT*(this->H(rho, T) - this->H(this->rhostd(),specie::Tstd)) + Cp_)/(sqr(deltaT) + 1);
+
+    return mu(T)*Cv_*(1.32 + 1.77*this->R()/Cv_)/CpBar;
+}
 // * * * * * * * * * * * * * * * Member Operators  * * * * * * * * * * * * * //
 
 template<class thermo>

tutorials/compressible/rhoPisoFoam/ras/pipe/0/alphat

@@ -1,53 +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       volScalarField;
-    location    "0";
-    object      alphat;
-}
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-dimensions      [1 -1 -1 0 0 0 0];
-
-internalField   uniform 0;
-
-boundaryField
-{
-    WALL_1
-    {
-        type            alphatWallFunction;
-        Prt             0.85;
-        value           uniform 0;
-    }
-    WALL_2
-    {
-        type            alphatWallFunction;
-        Prt             0.85;
-        value           uniform 0;
-    }
-    INLET
-    {
-        type            calculated;
-        value           uniform 0;
-    }
-    OUTLET
-    {
-        type            calculated;
-        value           uniform 0;
-    }
-    frontAndBackPlanes
-    {
-        type            empty;
-    }
-}
-
-
-// ************************************************************************* //

tutorials/compressible/rhoPisoFoam/ras/pipe/0/mut

@@ -1,57 +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       volScalarField;
-    location    "0";
-    object      mut;
-}
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-dimensions      [1 -1 -1 0 0 0 0];
-
-internalField   uniform 0;
-
-boundaryField
-{
-    WALL_1
-    {
-        type            mutWallFunction;
-        Cmu             0.09;
-        kappa           0.41;
-        E               9.8;
-        value           uniform 0;
-    }
-    WALL_2
-    {
-        type            mutWallFunction;
-        Cmu             0.09;
-        kappa           0.41;
-        E               9.8;
-        value           uniform 0;
-    }
-    INLET
-    {
-        type            calculated;
-        value           uniform 0;
-    }
-    OUTLET
-    {
-        type            calculated;
-        value           uniform 0;
-    }
-    frontAndBackPlanes
-    {
-        type            empty;
-    }
-}
-
-
-// ************************************************************************* //

tutorials/compressible/rhoPisoFoam/ras/pipe/constant/thermophysicalProperties

@@ -22,7 +22,7 @@ thermoType      realGasHThermo<pureMixture<realGasSutherlandTransport<realGasSpe
 mixture         CO2 1 44.01 73.773e5 304.13 0.22394 467.6 49436.5054 -626.411601 5.30172524 0.002503813816 -0.0000002127308728 -0.000000000768998878 2.849677801E-13 1.4792e-06 116;
 
 //thermoType      realGasHThermo<pureMixture<realGasSutherlandTransport<realGasSpecieThermo<nasaHeatCapacityPolynomial<redlichKwong>>>>>;
-//mixture         CO2 1 44.01 73.773e5 304.13 0.22394 49436.5054 -626.411601 5.30172524 0.002503813816 -0.0000002127308728 -0.000000000768998878 2.849677801e-13 1.4792e-06 116;
+//mixture         CO2 1 44.01 73.773e5 304.13  49436.5054 -626.411601 5.30172524 0.002503813816 -0.0000002127308728 -0.000000000768998878 2.849677801e-13 1.4792e-06 116;
 
 //thermoType      realGasHThermo<pureMixture<realGasSutherlandTransport<realGasSpecieThermo<nasaHeatCapacityPolynomial<pengRobinson>>>>>;
 //mixture        CO2 1 44.01 73.773e5 304.13 0.22394 49436.5054 -626.411601 5.30172524 0.002503813816 -0.0000002127308728 -0.000000000768998878 2.849677801e-13 1.4792e-06 116;
@@ -34,7 +34,7 @@ mixture         CO2 1 44.01 73.773e5 304.13 0.22394 467.6 49436.5054 -626.411601
 //mixture         CO2 1 44.01 73.773e5 304.13 0.22394 467.6 49436.5054 -626.411601 5.30172524 0.002503813816 -0.0000002127308728 -0.000000000768998878 2.849677801E-13 1e-06 0.7;
 
 //thermoType      realGasHThermo<pureMixture<realGasConstTransport<realGasSpecieThermo<nasaHeatCapacityPolynomial<redlichKwong>>>>>;
-//mixture         CO2 1 44.01 73.773e5 304.13 0.22394 49436.5054 -626.411601 5.30172524 0.002503813816 -0.0000002127308728 -0.000000000768998878 2.849677801e-13 1e-06 0.7;
+//mixture         CO2 1 44.01 73.773e5 304.13  49436.5054 -626.411601 5.30172524 0.002503813816 -0.0000002127308728 -0.000000000768998878 2.849677801e-13 1e-06 0.7;
 
 //thermoType      realGasHThermo<pureMixture<realGasConstTransport<realGasSpecieThermo<nasaHeatCapacityPolynomial<pengRobinson>>>>>;
 //mixture        CO2 1 44.01 73.773e5 304.13 0.22394 49436.5054 -626.411601 5.30172524 0.002503813816 -0.0000002127308728 -0.000000000768998878 2.849677801e-13 1e-06 0.7;