realGasEThermo + realGas_pEqn 4 rhoPisoFoam

applications/solvers/compressible/rhoPisoFoam/createFields.H

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

applications/solvers/compressible/rhoPisoFoam/pEqn.H

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

src/thermophysicalModels/basic/Make/files

@@ -9,6 +9,7 @@ psiThermo/hPsiThermo/hPsiThermos.C
 psiThermo/hsPsiThermo/hsPsiThermos.C
 psiThermo/ePsiThermo/ePsiThermos.C
 psiThermo/realGasHThermo/realGasHThermos.C
+psiThermo/realGasEThermo/realGasEThermos.C
 
 rhoThermo/basicRhoThermo/basicRhoThermo.C
 rhoThermo/basicRhoThermo/newBasicRhoThermo.C

src/thermophysicalModels/basic/psiThermo/basicPsiThermo/basicPsiThermo.C

@@ -48,4 +48,18 @@ Foam::basicPsiThermo::~basicPsiThermo()
 {}
 
 
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+const Foam::volScalarField& Foam::basicPsiThermo::drhodh() const
+{
+    notImplemented("basicPsiThermo::drhodh()");
+    return const_cast<volScalarField&>(volScalarField::null());
+}
+
+const Foam::volScalarField& Foam::basicPsiThermo::drhode() const
+{
+    notImplemented("basicPsiThermo::drhode()");
+    return const_cast<volScalarField&>(volScalarField::null());
+}
+
 // ************************************************************************* //

src/thermophysicalModels/basic/psiThermo/basicPsiThermo/basicPsiThermo.H

@@ -108,6 +108,12 @@ public:
             {
                 return p_*psi();
             }
+
+            //CL: drhodh needed for pressure equation of the real gas solver 
+            virtual const volScalarField& drhodh() const;
+
+            //CL: drhode needed for pressure equation of the real gas solver 
+            virtual const volScalarField& drhode() const;
 };
 
 

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

@@ -43,6 +43,7 @@ void Foam::realGasHThermo<MixtureType>::calculate()
     scalarField& TCells = this->T_.internalField();
     scalarField& rhoCells= this->rho_.internalField();
     scalarField& psiCells = this->psi_.internalField();
+    scalarField& drhodhCells = this->drhodh_.internalField();
     scalarField& muCells = this->mu_.internalField();
     scalarField& alphaCells = this->alpha_.internalField();
 
@@ -54,6 +55,7 @@ void Foam::realGasHThermo<MixtureType>::calculate()
 
         mixture_.TH(hCells[celli], TCells[celli], pCells[celli], rhoCells[celli]);
         psiCells[celli]=mixture_.psi(rhoCells[celli], TCells[celli]);
+        drhodhCells[celli]=mixture_.drhodH(rhoCells[celli], TCells[celli]);
         muCells[celli] = mixture_.mu(TCells[celli]);
         alphaCells[celli] = mixture_.alpha(rhoCells[celli], TCells[celli]);
     }
@@ -64,9 +66,9 @@ void Foam::realGasHThermo<MixtureType>::calculate()
         fvPatchScalarField& pp = this->p_.boundaryField()[patchi];
         fvPatchScalarField& pT = this->T_.boundaryField()[patchi];
         fvPatchScalarField& ppsi = this->psi_.boundaryField()[patchi];
+        fvPatchScalarField& pdrhodh = this->drhodh_.boundaryField()[patchi];
         fvPatchScalarField& prho = this->rho_.boundaryField()[patchi];
         fvPatchScalarField& ph = h_.boundaryField()[patchi];
-
         fvPatchScalarField& pmu = this->mu_.boundaryField()[patchi];
         fvPatchScalarField& palpha = this->alpha_.boundaryField()[patchi];
 
@@ -81,8 +83,8 @@ void Foam::realGasHThermo<MixtureType>::calculate()
 
                 prho[facei] = mixture_.rho(pp[facei], pT[facei],prho[facei]);
                 ppsi[facei]=mixture_.psi(prho[facei],pT[facei]);
+                pdrhodh[facei]=mixture_.drhodH(prho[facei],pT[facei]);
 		ph[facei] = mixture_.H(prho[facei], pT[facei]);
-
                 pmu[facei] = mixture_.mu(pT[facei]);
                 palpha[facei] = mixture_.alpha(prho[facei],pT[facei]);
             }
@@ -93,9 +95,11 @@ void Foam::realGasHThermo<MixtureType>::calculate()
             {
                 const typename MixtureType::thermoType& mixture_ =
                     this->patchFaceMixture(patchi, facei);
+
                 mixture_.TH(ph[facei], pT[facei],pp[facei],prho[facei]);
                 pmu[facei] = mixture_.mu(pT[facei]);
 		ppsi[facei]=mixture_.psi(prho[facei],pT[facei]);
+		pdrhodh[facei]=mixture_.drhodH(prho[facei],pT[facei]);
                 palpha[facei] = mixture_.alpha(prho[facei],pT[facei]);
             }
         }
@@ -121,11 +125,11 @@ Foam::realGasHThermo<MixtureType>::realGasHThermo(const fvMesh& mesh)
             IOobject::NO_READ,
             IOobject::NO_WRITE
         ),
-        mesh,
+        mesh, 
         dimensionSet(0, 2, -2, 0, 0),
         this->hBoundaryTypes()
-    ),
+    ),   
-    
+
     rho_
     (
         IOobject
@@ -138,7 +142,23 @@ Foam::realGasHThermo<MixtureType>::realGasHThermo(const fvMesh& mesh)
         ),
 	mesh,
 	dimDensity
+    ),
+
+    drhodh_
+    (
+        IOobject
+        (
+            "drhodh",
+            mesh.time().timeName(),
+            mesh,
+            IOobject::NO_READ,
+            IOobject::NO_WRITE
+        ),
+	mesh,
+	dimensionSet(1, -5, 2, 0, 0)
     )
+
+
 {
 
     scalarField& hCells = h_.internalField();
@@ -150,15 +170,13 @@ Foam::realGasHThermo<MixtureType>::realGasHThermo(const fvMesh& mesh)
     forAll(rhoCells, celli)
     {
      rhoCells[celli]=this->cellMixture(celli).rho(pCells[celli],TCells[celli]);
-     
     }
 
 
     forAll(rho_.boundaryField(), patchi)
     {
         rho_.boundaryField()[patchi] ==
-            rho(this->T_.boundaryField()[patchi], patchi);
+            rho(this->T_.boundaryField()[patchi], patchi); 
-   
     }
 
 
@@ -218,7 +236,7 @@ Foam::tmp<Foam::scalarField> Foam::realGasHThermo<MixtureType>::h
     const labelList& cells
 ) const
 {
-    //CL: needing the pressure of the internal field to calculate the realGas enthalpy
+    //CL: need the pressure of the internal field to calculate the realGas enthalpy
     //CL: this is done this way to assure compatibility to old OF Thermo-Versions 
     const scalarField& pCells = this->p_.internalField();
 
@@ -240,7 +258,7 @@ Foam::tmp<Foam::scalarField> Foam::realGasHThermo<MixtureType>::h
     const label patchi
 ) const
 {
-    //CL: needing the pressure at the patch to calculate the realGas enthalpy
+    //CL: need the pressure at the patch to calculate the realGas enthalpy
     //CL: this is done this way to assure compatibility to old OF Thermo-Versions 
     const fvPatchScalarField& pp = this->p_.boundaryField()[patchi];
 
@@ -262,7 +280,7 @@ Foam::tmp<Foam::scalarField> Foam::realGasHThermo<MixtureType>::rho
     const label patchi
 ) const
 {
-    //CL: needing the pressure at the patch to calculate the realGas enthalpy
+    //CL: need the pressure at the patch to calculate the realGas enthalpy
     //CL: this is done this way to assure compatibility to old OF Thermo-Versions 
     const fvPatchScalarField& pp = this->p_.boundaryField()[patchi];
 
@@ -284,7 +302,7 @@ Foam::tmp<Foam::scalarField> Foam::realGasHThermo<MixtureType>::Cp
     const label patchi
 ) const
 {
-    //CL: needing the pressure at the patch to calculate the realGas enthalpy
+    //CL: need the pressure at the patch to calculate the realGas enthalpy
     //CL: this is done this way to assure compatibility to old OF Thermo-Versions 
     const fvPatchScalarField& pp = this->p_.boundaryField()[patchi];
 
@@ -352,7 +370,7 @@ Foam::tmp<Foam::scalarField> Foam::realGasHThermo<MixtureType>::Cv
     const label patchi
 ) const
 {
-    //CL: needing the pressure at the patch to calculate the realGas enthalpy
+    //CL: need the pressure at the patch to calculate the realGas enthalpy
     //CL: this is done this way to assure compatibility to old OF Thermo-Versions 
     const fvPatchScalarField& pp = this->p_.boundaryField()[patchi];
 
@@ -374,7 +392,7 @@ Foam::tmp<Foam::volScalarField> Foam::realGasHThermo<MixtureType>::rho()  const
 
     const fvMesh& mesh = this->T_.mesh();
 
-    //CL: create an rho Field, which will be return
+    //CL: create a rho Field, which will be return
     //CL: the problem is that this function is "const", 
     //CL: so a new variabel is needed
       tmp<volScalarField> trho

src/thermophysicalModels/basic/psiThermo/realGasHThermo/realGasHThermo.H

@@ -53,7 +53,7 @@ namespace Foam
 {
 
 /*---------------------------------------------------------------------------*\
-                         Class hPsiThermo Declaration
+                         Class realGasHThermo Declaration
 \*---------------------------------------------------------------------------*/
 
 template<class MixtureType>
@@ -70,6 +70,8 @@ class realGasHThermo
         //- DensityField
         volScalarField rho_;
 
+        //- drhodh_Field
+        volScalarField drhodh_;
 
     // Private member functions
 
@@ -128,6 +130,12 @@ public:
                 return h_;
             }
 
+            //CL: drhodh needed for pressure equation of the real gas solver 
+            virtual const volScalarField& drhodh() const
+            {
+                return drhodh_;
+            }
+
         // Fields derived from thermodynamic state variables
 
             //- Enthalpy for cell-set [J/kg]

src/thermophysicalModels/specie/thermo/realGasSpecieThermo/realGasSpecieThermo.C

@@ -28,10 +28,8 @@ Institut für Thermodynamik
 Technische Universität Braunschweig 
 Germany
 
-
 \*---------------------------------------------------------------------------*/
 
-
 #include "realGasSpecieThermo.H"
 #include "IOstreams.H"
 

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

@@ -40,8 +40,6 @@ SourceFiles
     realGasSpecieThermoI.H
     realGasSpecieThermo.C
 
-
-
 Author
 Christian Lucas
 Institut für Thermodynamik
@@ -100,7 +98,7 @@ Ostream& operator<<
 
 
 /*---------------------------------------------------------------------------*\
-                           Class specieThermo Declaration
+                           Class realGasSpecieThermo Declaration
 \*---------------------------------------------------------------------------*/
 
 template<class thermo>
@@ -120,7 +118,7 @@ class realGasSpecieThermo
     // Private member functions
 
         // return the temperature corresponding to the value of the
-        //  thermodynamic property f, given the function f = F(T) and dF(T)/dT
+        //  thermodynamic property f, given the function f = F(rho,T) and dF(rho,T)/dT
         inline void T
         (
             scalar f,
@@ -131,12 +129,8 @@ class realGasSpecieThermo
             scalar (realGasSpecieThermo::*dFdT)(const scalar,const scalar) const
         ) const;
 
-  
-
 public:
 
-
-
         // Constructors
 
         //- construct from components
@@ -154,35 +148,17 @@ public:
         // Fundamaental properties
         // (These functions must be provided in derived types)
 
-            // Heat capacity at constant pressure [J/(kmol K)]
-            //scalar cp(const scalar) const;
-
-            // Enthalpy [J/kmol]
-            //scalar h(const scalar) const;
-
             // Sensible enthalpy [J/kmol]
             //virtual scalar hs(const scalar) const;
 
             // Chemical enthalpy [J/kmol]
             //virtual scalar hc(const scalar) const;
 
-            // Entropy [J/(kmol K)]
-            //virtual scalar s(const scalar) const;
-
-
         // Calculate and return derived properties
         // (These functions need not provided in derived types)
 
-            // Mole specific properties
+                //CL: isentropic expansion factor "gamma" (heat capacity ratio for perfect gas) 
-
+                inline scalar gamma(const scalar T, const scalar rho) const;
-                //- Heat capacity at constant volume [J/(kmol K)]
-                //inline scalar cv(const scalar T, const scalar rho) const;
-
-                //- gamma = cp/cv []
-                //inline scalar gamma(const scalar T, const scalar rho) const;
-
-                // Internal energy [J/kmol]
-                // inline scalar e(const scalar rho, const scalar T) const;
 
                 //- Sensible internal energy [J/kmol]
                 // inline scalar es(const scalar p, const scalar rho) const;
@@ -223,7 +199,7 @@ public:
                 //- Helmholtz free energy [J/kg]
                 inline scalar A(const scalar rho, const scalar T) const;
 
-                //Other variables
+                //CL: Other variables
 
 	        //- Return compressibility drho/dp at h=constant [s^2/m^2] 
                 inline scalar psiH(const scalar rho, const scalar T) const;
@@ -231,14 +207,17 @@ public:
                 //- Return compressibility drho/dp at e=constant [s^2/m^2]     
                 inline scalar psiE(const scalar rho, const scalar T) const;
 
+	        //- Return compressibility drho/dH at p=constant 
+                inline scalar drhodH(const scalar rho, const scalar T) const;
+
+	        //- Return compressibility drho/dE at p=constant 
+                inline scalar drhodE(const scalar rho, const scalar T) const;
 
         // Energy->temperature  inversion functions
 
             //- Temperature from Enthalpy given an initial temperature T0
             inline void TH(const scalar H,  scalar &T0,const scalar p, scalar &psi0) const;
 
-
-
             //- Temperature from internal energy given an initial temperature T0
             inline void TE(const scalar E, scalar &T0,const scalar p, scalar &psi0) const;
 

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

@@ -43,11 +43,10 @@ inline Foam::realGasSpecieThermo<thermo>::realGasSpecieThermo
 {}
 
 
-
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
 
-// using two one dimensional newton solvers in a row
+//CL: using two one dimensional newton solvers in a row
 template<class thermo>
 inline void Foam::realGasSpecieThermo<thermo>::T
 (
@@ -76,8 +75,8 @@ inline void Foam::realGasSpecieThermo<thermo>::T
         i=0;
       	do
    	{   
-            // using a stabilizing newton solver 
+            //CL: using a stabilizing newton solver 
-            // if the solve is diverging, the time step is reduced until the solver converges
+            //CL: if the solve is diverging, the time step is reduced until the solver converges
             Tnew = Test - ((this->*F)(rho,Test) - f)/(this->*dFdT)(rho,Test)/(pow(2,i)); 
 	    i++;					
         }while
@@ -103,7 +102,7 @@ inline void Foam::realGasSpecieThermo<thermo>::T
                 << abort(FatalError);
         }
     } while
-    // both fields must converge
+    //CL: both fields must converge
     ( 
         (mag(mag(Tnew) - mag(Test)) > Ttol)
         ||
@@ -130,18 +129,24 @@ inline Foam::realGasSpecieThermo<thermo>::realGasSpecieThermo
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
 
+template<class thermo>
+inline Foam::scalar Foam::realGasSpecieThermo<thermo>::gamma(const scalar rho,  const scalar T ) const
+{
+    return -1/(rho*this->p(rho,T))*this->cp(rho,T)/this->cv(rho,T)*this->dpdv(rho,T);
+}
+
 
 template<class thermo>
 inline Foam::scalar Foam::realGasSpecieThermo<thermo>::g(const scalar rho,  const scalar T ) const
 {
-    return this->h(rho, this->pReturn(rho,T)) - T*this->s(rho, this->pReturn(rho,T));
+    return this->h(rho, this->p(rho,T)) - T*this->s(rho, this->p(rho,T));
 }
 
 
 template<class thermo>
 inline Foam::scalar Foam::realGasSpecieThermo<thermo>::a(const scalar rho, const scalar T ) const
 {
-    return this->e(rho,this->pReturn(rho,T)) - T*this->s(rho,this->pReturn(rho,T));
+    return this->e(rho,this->p(rho,T)) - T*this->s(rho,this->p(rho,T));
 }
 
 
@@ -166,9 +171,6 @@ inline Foam::scalar Foam::realGasSpecieThermo<thermo>::H(const scalar rho, const
 }
 
 
-
-
-
 template<class thermo>
 inline Foam::scalar Foam::realGasSpecieThermo<thermo>::S(const scalar rho, const scalar T) const
 {
@@ -196,8 +198,8 @@ 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]
+//CL:- Return compressibility drho/dp at h=constant [s^2/m^2]
-//- using Bridgeman's Table
+//CL:- using Bridgeman's Table
 template<class thermo>
 inline Foam::scalar Foam::realGasSpecieThermo<thermo>::psiH
 (
@@ -206,17 +208,17 @@ inline Foam::scalar Foam::realGasSpecieThermo<thermo>::psiH
     
 ) const
 {
-
     scalar beta=this->isobarExpCoef(rho,T);
+
     return  
     -(
-        (T*pow(beta,2)-beta)/this->cp(rho,T)
+        (T*beta*beta-beta)/this->Cp(rho,T)
-        -this->isothermalCompressiblity(rho,T)*rho/this->W()
+        -this->isothermalCompressiblity(rho,T)*rho
-    )*this->W();
+    );
 }
 
-//- Return compressibility drho/dp at e=constant [s^2/m^2]
+//CL:- Return compressibility drho/dp at e=constant [s^2/m^2]
-//- using Bridgeman's Table
+//CL:- using Bridgeman's Table
 template<class thermo>
 inline Foam::scalar Foam::realGasSpecieThermo<thermo>::psiE
 (
@@ -225,21 +227,50 @@ inline Foam::scalar Foam::realGasSpecieThermo<thermo>::psiE
     
 ) const
 {
-    scalar Vm = this->W()/rho;
+    scalar V = 1/rho;
-    scalar cp=this->cp(rho,T);
+    scalar cp=this->Cp(rho,T);
     scalar beta=this->isobarExpCoef(rho,T);
+
     return  
     -(
         (
-	    T*pow(beta,2)*Vm
+	    T*pow(beta,2)*V
 	    -this->isothermalCompressiblity(rho,T)*cp
 	)
 	/
 	(
-            cp*Vm
+            cp*V
-            -beta*this->p(rho,T)*pow(Vm,2)
+            -beta*this->p(rho,T)*pow(V,2)
 	)
-    )*this->W();
+    );
+}
+
+//CL:- Returns drho/dH at p=constant 
+//CL:- using Bridgeman's Table
+template<class thermo>
+inline Foam::scalar Foam::realGasSpecieThermo<thermo>::drhodH
+(
+    const scalar rho,
+    const scalar T
+    
+) const
+{
+    return -(rho*this->isobarExpCoef(rho,T))/this->Cp(rho,T);
+}
+
+//CL:- Returns drho/dE at p=constant 
+//CL:- using Bridgeman's Table
+template<class thermo>
+inline Foam::scalar Foam::realGasSpecieThermo<thermo>::drhodE
+(
+    const scalar rho,
+    const scalar T
+    
+) const
+{
+    scalar beta=this->isobarExpCoef(rho,T);
+
+    return -(rho*beta)/(this->Cp(rho,T)-beta*this->p(rho,T)/rho);
 }
 
 template<class thermo>

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

@@ -116,6 +116,7 @@ inline scalar constTransport<thermo>::alpha(const scalar T) const
     return Cp_*mu(T)*rPr/CpBar;
 }
 
+// CL: for real gas thermo
 // Thermal conductivity [W/mK]
 template<class thermo>
 inline scalar constTransport<thermo>::kappa(const scalar rho,const scalar T) const
@@ -123,7 +124,7 @@ inline scalar constTransport<thermo>::kappa(const scalar rho,const scalar T) con
     return this->Cp(rho,T)*mu(T)*rPr;
 }
 
-
+// CL: for real gas thermo
 // Thermal diffusivity for enthalpy [kg/ms]
 template<class thermo>
 inline scalar constTransport<thermo>::alpha(const scalar rho,const scalar T) const

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

@@ -155,6 +155,7 @@ inline scalar sutherlandTransport<thermo>::alpha(const scalar T) const
     return mu(T)*Cv_*(1.32 + 1.77*this->R()/Cv_)/CpBar;
 }
 
+// CL: for real gas thermo
 // Thermal conductivity [W/mK]
 template<class thermo>
 inline scalar sutherlandTransport<thermo>::kappa(const scalar rho, const scalar T) const
@@ -163,7 +164,7 @@ inline scalar sutherlandTransport<thermo>::kappa(const scalar rho, const scalar
     return mu(T)*Cv_*(1.32 + 1.77*this->R()/Cv_);
 }
 
-
+// CL: for real gas thermo
 // Thermal diffusivity for enthalpy [kg/ms]
 template<class thermo>
 inline scalar sutherlandTransport<thermo>::alpha(const scalar rho, const scalar T) const

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

@@ -18,28 +18,28 @@ FoamFile
 
 
 
-thermoType      realGasHThermo<pureMixture<realGasSutherlandTransport<realGasSpecieThermo<nasaHeatCapacityPolynomial<aungierRedlichKwong>>>>>;
+thermoType      realGasHThermo<pureMixture<sutherlandTransport<realGasSpecieThermo<nasaHeatCapacityPolynomial<aungierRedlichKwong>>>>>;
 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>>>>>;
+//thermoType      realGasHThermo<pureMixture<sutherlandTransport<realGasSpecieThermo<nasaHeatCapacityPolynomial<redlichKwong>>>>>;
-//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;
-
-//thermoType      realGasHThermo<pureMixture<realGasSutherlandTransport<realGasSpecieThermo<nasaHeatCapacityPolynomial<soaveRedlichKwong>>>>>;
 //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;
 
-//thermoType      realGasHThermo<pureMixture<realGasConstTransport<realGasSpecieThermo<nasaHeatCapacityPolynomial<aungierRedlichKwong>>>>>;
+//thermoType      realGasHThermo<pureMixture<sutherlandTransport<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;
+
+//thermoType      realGasHThermo<pureMixture<sutherlandTransport<realGasSpecieThermo<nasaHeatCapacityPolynomial<soaveRedlichKwong>>>>>;
+//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;
+
+//thermoType      realGasHThermo<pureMixture<constTransport<realGasSpecieThermo<nasaHeatCapacityPolynomial<aungierRedlichKwong>>>>>;
 //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>>>>>;
+//thermoType      realGasHThermo<pureMixture<constTransport<realGasSpecieThermo<nasaHeatCapacityPolynomial<redlichKwong>>>>>;
-//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;
+//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;
 
-//thermoType      realGasHThermo<pureMixture<realGasConstTransport<realGasSpecieThermo<nasaHeatCapacityPolynomial<pengRobinson>>>>>;
+//thermoType      realGasHThermo<pureMixture<constTransport<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;
 
-//thermoType      realGasHThermo<pureMixture<realGasConstTransport<realGasSpecieThermo<nasaHeatCapacityPolynomial<soaveRedlichKwong>>>>>;
+//thermoType      realGasHThermo<pureMixture<constTransport<realGasSpecieThermo<nasaHeatCapacityPolynomial<soaveRedlichKwong>>>>>;
 //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;
 
 
@@ -52,7 +52,7 @@ mixture         CO2 1 44.01 73.773e5 304.13 0.22394 467.6 49436.5054 -626.411601
 //   77.773e5 				--> critical pressure
 //   304.13 				--> critical temperatur
 //   0.22394 				--> acentric factor
-//   467.6     				--> critical density (only for aungier redlich kwong
+//   467.6     				--> critical density (only for aungier redlich kwong)
 //   49436.5054 --> 2.849677801e-13 	--> 7 heat capacity polynomial coefficent's
 //   ....				--> two coefficent's for sutherland model or for the constTransport model
 // *********************************************************************************************************************** //

tutorials/compressible/rhoPisoFoam/ras/pipe/system/controlDict

@@ -41,13 +41,13 @@ writeCompression uncompressed;
 
 timeFormat      general;
 
-timePrecision   20;
+timePrecision   10;
 
 adjustTimeStep  yes;
 
-maxCo          0.05;
+maxCo          0.5;
 
-maxDeltaT       1e-4;
+maxDeltaT       1e-2;
 
 runTimeModifiable yes;
 

tutorials/compressible/rhoPisoFoam/ras/pipe/system/fvSchemes

@@ -1,8 +1,8 @@
 /*--------------------------------*- C++ -*----------------------------------*\
 | =========                 |                                                 |
-| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+| \\      /  F ield         | OpenFOAM Extend Project: Open Source CFD        |
-|  \\    /   O peration     | Version:  1.6                                   |
+|  \\    /   O peration     | Version:  1.6-ext                               |
-|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|   \\  /    A nd           | Web:      www.extend-project.de                 |
 |    \\/     M anipulation  |                                                 |
 \*---------------------------------------------------------------------------*/
 FoamFile
@@ -22,24 +22,25 @@ ddtSchemes
 
 gradSchemes
 {
- default         cellMDLimited Gauss linear 0.333;
+    default         Gauss linear;
- grad(p)         cellMDLimited Gauss linear 0.333;
+    grad(p)         Gauss linear;
 }
 
 divSchemes
 {
     default         none;
-    div(phi,U)      Gauss limitedLinearV 0.5;
+    div(phi,U)      Gauss limitedLinearV 1;
-    div(phid,p)     Gauss limitedLinear  0.5;
+    div(phid,p)     Gauss limitedLinear 1;
-    div(phiU,p)     Gauss limitedLinear  0.5;
+    div(phiU,p)     Gauss linear;
-    div(phi,h)      Gauss limitedLinear  0.5;
+    div(phi,h)      Gauss limitedLinear 1;
-    div(phi,k)      Gauss limitedLinear  0.5;
+    div(phi,k)      Gauss limitedLinear 1;
-    div(phi,epsilon) Gauss limitedLinear 0.5;
+    div(phi,epsilon) Gauss limitedLinear 1;
-    div(phi,R)      Gauss limitedLinear 0.5;
+    div(phi,R)      Gauss limitedLinear 1;
-    div(phi,omega)  Gauss limitedLinear  0.5;;
+    div(phi,omega)  Gauss limitedLinear 1;
-    div(U)          Gauss limitedLinear  0.5;
+    div((rho*R))    Gauss linear;
+    div(R)          Gauss linear;
+    div(U)          Gauss linear;
     div((muEff*dev2(grad(U).T()))) Gauss linear;
-
 }
 
 laplacianSchemes
@@ -53,7 +54,6 @@ laplacianSchemes
     laplacian(DomegaEff,omega) Gauss linear corrected;
     laplacian((rho*(1|A(U))),p) Gauss linear corrected;
     laplacian(alphaEff,h) Gauss linear corrected;
-
 }
 
 interpolationSchemes
@@ -72,5 +72,4 @@ fluxRequired
     p               ;
 }
 
-
 // ************************************************************************* //

tutorials/compressible/rhoPisoFoam/ras/pipe/system/fvSolution

@@ -35,8 +35,6 @@ p
         mergeLevels     1;
 }
 
-
-
     U
     {
         solver          smoothSolver;
@@ -90,10 +88,10 @@ p
 
 PISO
 {
+    realFluid true;
     nNonOrthogonalCorrectors 0;
-    nCorrectors     6;
+    nCorrectors     2;
     momentumPredictor yes;
-   // transonic true;
 }
 
 

tutorials/compressible/rhoPisoFoam/ras/pipe/system/sampleDict

@@ -1,44 +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       dictionary;
-    location    "system";
-    object      sampleDict;
-}
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-interpolationScheme cellPoint;
-
-setFormat     raw;
-
-sets
-(
-    leftPatch
-    {
-
-        type            uniform;
-        axis            y;
-        start           (0.49 0 0.0);
-        end             (0.49 0.01 0);
-        nPoints         100;
-    }
-);
-
-fields
-(
-    U p T
-);
-
-surfaces
-(
-);
-
-// ************************************************************************* //