Improved steadyCompressibleMRFFoam

applications/solvers/compressible/steadyCompressibleMRFFoam/createFields.H

@@ -9,7 +9,8 @@
     volScalarField& p = thermo.p();
     volScalarField& h = thermo.h();
     const volScalarField& T = thermo.T();
-    const volScalarField& psi = thermo.psi();
+    volScalarField psis("psi", thermo.psi()/thermo.Cp()*thermo.Cv());
+    psis.oldTime();
 
     volScalarField rho
     (

applications/solvers/compressible/steadyCompressibleMRFFoam/hEqn.H

@@ -25,4 +25,5 @@
 
     // Bounding of enthalpy taken out
     thermo.correct();
+    psis = thermo.psi()/thermo.Cp()*thermo.Cv();
 }

applications/solvers/compressible/steadyCompressibleMRFFoam/iEqn.H

@@ -53,4 +53,5 @@
     i = h - 0.5*magSqr(Urot);
 
     thermo.correct();
+    psis = thermo.psi()/thermo.Cp()*thermo.Cv();
 }

applications/solvers/compressible/steadyCompressibleMRFFoam/pEqn.H

@@ -1,29 +1,23 @@
 {
     volScalarField rUA = 1.0/UEqn.A();
 
+    surfaceScalarField psisf = fvc::interpolate(psis);
+    surfaceScalarField rhof = fvc::interpolate(rho);
+
+    // Needs to be outside of loop since p is changing, but psi and rho are not.
+    surfaceScalarField rhoReff = rhof - psisf*fvc::interpolate(p);
+
     for (int corr = 0; corr < nCorr; corr++)
     {
         U = rUA*UEqn.H();
 
-        surfaceScalarField psif = fvc::interpolate(psi);
-        surfaceScalarField rhof = fvc::interpolate(rho);
+        phi = rhoReff*(fvc::interpolate(U) & mesh.Sf());
 
-        // Execute ddtPhiCorr before recalculating flux
-        // HJ, 27/Apr/2010
-        surfaceScalarField phid
-        (
-            "phid",
-            psif*(fvc::interpolate(U) & mesh.Sf())
-        );
+        surfaceScalarField phid("phid", psisf/rhoReff*phi);
 
-        // Make flux relative within the MRF zone
-        mrfZones.relativeFlux(psif, phid);
-
-        // Calculate phi for boundary conditions
-        phi = fvc::interpolate(rho*U) & mesh.Sf();
-
-        // Make flux relative within the MRF zone
-        mrfZones.relativeFlux(rhof, phi);
+        // Make fluxes relative within the MRF zone
+        mrfZones.relativeFlux(psisf, phid);
+        mrfZones.relativeFlux(rhoReff, phi);
 
         p.storePrevIter();
 
@@ -31,8 +25,9 @@
         {
             fvScalarMatrix pEqn
             (
-                fvm::ddt(psi, p)
+                fvm::ddt(psis, p)
               + fvm::div(phid, p)
+              + fvc::div(phi)
               - fvm::laplacian(rho*rUA, p)
             );
 
@@ -47,7 +42,7 @@
             // Calculate the flux
             if (nonOrth == nNonOrthCorr)
             {
-                phi = pEqn.flux();
+                phi += pEqn.flux();
             }
         }
 

src/finiteVolume/fields/fvPatchFields/derived/totalPressure/totalPressureFvPatchScalarField.C

@@ -254,6 +254,59 @@ void Foam::totalPressureFvPatchScalarField::updateCoeffs()
 }
 
 
+Foam::tmp<Foam::scalarField>
+Foam::totalPressureFvPatchScalarField::snGrad() const
+{
+    return tmp<scalarField>
+    (
+        new scalarField(this->size(), 0.0)
+    );
+}
+
+
+Foam::tmp<Foam::scalarField>
+Foam::totalPressureFvPatchScalarField::valueInternalCoeffs
+(
+    const tmp<scalarField>&
+) const
+{
+    return tmp<scalarField>
+    (
+        new scalarField(this->size(), 0.0)
+    );
+}
+
+
+Foam::tmp<Foam::scalarField>
+Foam::totalPressureFvPatchScalarField::valueBoundaryCoeffs
+(
+    const tmp<scalarField>&
+) const
+{
+    return *this;
+}
+
+
+Foam::tmp<Foam::scalarField>
+Foam::totalPressureFvPatchScalarField::gradientInternalCoeffs() const
+{
+    return tmp<scalarField>
+    (
+        new scalarField(this->size(), 0.0)
+    );
+}
+
+
+Foam::tmp<Foam::scalarField>
+Foam::totalPressureFvPatchScalarField::gradientBoundaryCoeffs() const
+{
+    return tmp<scalarField>
+    (
+        new scalarField(this->size(), 0.0)
+    );
+}
+
+
 void Foam::totalPressureFvPatchScalarField::write(Ostream& os) const
 {
     fvPatchScalarField::write(os);

src/finiteVolume/fields/fvPatchFields/derived/totalPressure/totalPressureFvPatchScalarField.H

@@ -207,6 +207,30 @@ public:
             //- Update the coefficients associated with the patch field
             virtual void updateCoeffs();
 
+            //- Return patch-normal gradient
+            virtual tmp<scalarField> snGrad() const;
+
+            //- Return the matrix diagonal coefficients corresponding to the
+            //  evaluation of the value of this patchField with given weights
+            virtual tmp<scalarField> valueInternalCoeffs
+            (
+                const tmp<scalarField>&
+            ) const;
+
+            //- Return the matrix source coefficients corresponding to the
+            //  evaluation of the value of this patchField with given weights
+            virtual tmp<scalarField> valueBoundaryCoeffs
+            (
+                const tmp<scalarField>&
+            ) const;
+
+            //- Return the matrix diagonal coefficients corresponding to the
+            //  evaluation of the gradient of this patchField
+            virtual tmp<scalarField> gradientInternalCoeffs() const;
+
+            //- Return the matrix source coefficients corresponding to the
+            //  evaluation of the gradient of this patchField
+            virtual tmp<scalarField> gradientBoundaryCoeffs() const;
 
         //- Write
         virtual void write(Ostream&) const;

src/thermophysicalModels/basic/derivedFvPatchFields/isentropicTotalPressure/isentropicTotalPressureFvPatchScalarField.C

@@ -184,6 +184,59 @@ void Foam::isentropicTotalPressureFvPatchScalarField::updateCoeffs
 }
 
 
+Foam::tmp<Foam::scalarField>
+Foam::isentropicTotalPressureFvPatchScalarField::snGrad() const
+{
+    return tmp<scalarField>
+    (
+        new scalarField(this->size(), 0.0)
+    );
+}
+
+
+Foam::tmp<Foam::scalarField>
+Foam::isentropicTotalPressureFvPatchScalarField::valueInternalCoeffs
+(
+    const tmp<scalarField>&
+) const
+{
+    return tmp<scalarField>
+    (
+        new scalarField(this->size(), 0.0)
+    );
+}
+
+
+Foam::tmp<Foam::scalarField>
+Foam::isentropicTotalPressureFvPatchScalarField::valueBoundaryCoeffs
+(
+    const tmp<scalarField>&
+) const
+{
+    return *this;
+}
+
+
+Foam::tmp<Foam::scalarField>
+Foam::isentropicTotalPressureFvPatchScalarField::gradientInternalCoeffs() const
+{
+    return tmp<scalarField>
+    (
+        new scalarField(this->size(), 0.0)
+    );
+}
+
+
+Foam::tmp<Foam::scalarField>
+Foam::isentropicTotalPressureFvPatchScalarField::gradientBoundaryCoeffs() const
+{
+    return tmp<scalarField>
+    (
+        new scalarField(this->size(), 0.0)
+    );
+}
+
+
 void Foam::isentropicTotalPressureFvPatchScalarField::write
 (
     Ostream& os

src/thermophysicalModels/basic/derivedFvPatchFields/isentropicTotalPressure/isentropicTotalPressureFvPatchScalarField.H

@@ -134,6 +134,31 @@ public:
 
         // Access
 
+            //- Return the name of the velocity field
+            const word& UName() const
+            {
+                return UName_;
+            }
+
+            //- Return reference to the name of the velocity field
+            //  to allow adjustment
+            word& UName()
+            {
+                return UName_;
+            }
+
+            //- Return the heat capacity ratio
+            scalar gamma() const
+            {
+                return gamma_;
+            }
+
+            //- Return reference to the heat capacity ratio to allow adjustment
+            scalar& gamma()
+            {
+                return gamma_;
+            }
+
             //- Return the total pressure
             const scalarField& p0() const
             {