rhoCentralFoam: Corrected the BCs for positive (outgoing) fluxes

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

applications/solvers/compressible/rhoCentralFoam/createFields.H

@@ -8,6 +8,7 @@ psiThermo& thermo = pThermo();
 
 volScalarField& p = thermo.p();
 volScalarField& e = thermo.he();
+const volScalarField& T = thermo.T();
 const volScalarField& psi = thermo.psi();
 const volScalarField& mu = thermo.mu();
 

applications/solvers/compressible/rhoCentralFoam/directionInterpolate.H

@@ -0,0 +1,45 @@
+namespace Foam
+{
+
+//- Interpolate field vf according to direction dir
+template<class Type>
+tmp<GeometricField<Type, fvsPatchField, surfaceMesh> > interpolate
+(
+    const GeometricField<Type, fvPatchField, volMesh>& vf,
+    const surfaceScalarField& dir,
+    const word& reconFieldName = word::null
+)
+{
+    tmp<GeometricField<Type, fvsPatchField, surfaceMesh> > tsf
+    (
+        fvc::interpolate
+        (
+            vf,
+            dir,
+            "reconstruct("
+          + (reconFieldName != word::null ? reconFieldName : vf.name())
+          + ')'
+        )
+    );
+
+    GeometricField<Type, fvsPatchField, surfaceMesh>& sf = tsf();
+
+    sf.rename(vf.name() + '_' + dir.name());
+
+    // Correct BCs of the positive (outgoing) fluxes
+    if (dir[0] > 0)
+    {
+        forAll(sf.boundaryField(), patchi)
+        {
+            if (!sf.boundaryField()[patchi].coupled())
+            {
+                sf.boundaryField()[patchi] =
+                    vf.boundaryField()[patchi].patchInternalField();
+            }
+        }
+    }
+
+    return tsf;
+}
+
+}

applications/solvers/compressible/rhoCentralFoam/rhoCentralDyMFoam/rhoCentralDyMFoam.C

@@ -26,7 +26,7 @@ Application
 
 Description
     Density-based compressible flow solver based on central-upwind schemes of
-    Kurganov and Tadmor
+    Kurganov and Tadmor with support for mesh-motion and topology changes
 
 \*---------------------------------------------------------------------------*/
 
@@ -36,6 +36,7 @@ Description
 #include "turbulentFluidThermoModel.H"
 #include "zeroGradientFvPatchFields.H"
 #include "fixedRhoFvPatchScalarField.H"
+#include "directionInterpolate.H"
 #include "motionSolver.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@@ -72,52 +73,20 @@ int main(int argc, char *argv[])
         // Do any mesh changes
         mesh.update();
 
-        // --- upwind interpolation of primitive fields on faces
+        // --- Directed interpolation of primitive fields onto faces
 
-        surfaceScalarField rho_pos
+        surfaceScalarField rho_pos(interpolate(rho, pos));
-        (
+        surfaceScalarField rho_neg(interpolate(rho, neg));
-            "rho_pos",
-            fvc::interpolate(rho, pos, "reconstruct(rho)")
-        );
-        surfaceScalarField rho_neg
-        (
-            "rho_neg",
-            fvc::interpolate(rho, neg, "reconstruct(rho)")
-        );
 
-        surfaceVectorField rhoU_pos
+        surfaceVectorField rhoU_pos(interpolate(rhoU, pos, U.name()));
-        (
+        surfaceVectorField rhoU_neg(interpolate(rhoU, neg, U.name()));
-            "rhoU_pos",
-            fvc::interpolate(rhoU, pos, "reconstruct(U)")
-        );
-        surfaceVectorField rhoU_neg
-        (
-            "rhoU_neg",
-            fvc::interpolate(rhoU, neg, "reconstruct(U)")
-        );
 
-        volScalarField rPsi(1.0/psi);
+        volScalarField rPsi("rPsi", 1.0/psi);
-        surfaceScalarField rPsi_pos
+        surfaceScalarField rPsi_pos(interpolate(rPsi, pos, T.name()));
-        (
+        surfaceScalarField rPsi_neg(interpolate(rPsi, neg, T.name()));
-            "rPsi_pos",
-            fvc::interpolate(rPsi, pos, "reconstruct(T)")
-        );
-        surfaceScalarField rPsi_neg
-        (
-            "rPsi_neg",
-            fvc::interpolate(rPsi, neg, "reconstruct(T)")
-        );
 
-        surfaceScalarField e_pos
+        surfaceScalarField e_pos(interpolate(e, pos, T.name()));
-        (
+        surfaceScalarField e_neg(interpolate(e, neg, T.name()));
-            "e_pos",
-            fvc::interpolate(e, pos, "reconstruct(T)")
-        );
-        surfaceScalarField e_neg
-        (
-            "e_neg",
-            fvc::interpolate(e, neg, "reconstruct(T)")
-        );
 
         surfaceVectorField U_pos("U_pos", rhoU_pos/rho_pos);
         surfaceVectorField U_neg("U_neg", rhoU_neg/rho_neg);
@@ -135,16 +104,16 @@ int main(int argc, char *argv[])
             phiv_neg -= mesh.phi();
         }
 
-        volScalarField c(sqrt(thermo.Cp()/thermo.Cv()*rPsi));
+        volScalarField c("c", sqrt(thermo.Cp()/thermo.Cv()*rPsi));
         surfaceScalarField cSf_pos
         (
             "cSf_pos",
-            fvc::interpolate(c, pos, "reconstruct(T)")*mesh.magSf()
+            interpolate(c, pos, T.name())*mesh.magSf()
         );
         surfaceScalarField cSf_neg
         (
             "cSf_neg",
-            fvc::interpolate(c, neg, "reconstruct(T)")*mesh.magSf()
+            interpolate(c, neg, T.name())*mesh.magSf()
         );
 
         surfaceScalarField ap
@@ -206,7 +175,7 @@ int main(int argc, char *argv[])
             phiEp += mesh.phi()*(a_pos*p_pos + a_neg*p_neg);
         }
 
-        volScalarField muEff(turbulence->muEff());
+        volScalarField muEff("muEff", turbulence->muEff());
         volTensorField tauMC("tauMC", muEff*dev2(Foam::T(fvc::grad(U))));
 
         // --- Solve density

applications/solvers/compressible/rhoCentralFoam/rhoCentralFoam.C

@@ -35,6 +35,7 @@ Description
 #include "turbulentFluidThermoModel.H"
 #include "zeroGradientFvPatchFields.H"
 #include "fixedRhoFvPatchScalarField.H"
+#include "directionInterpolate.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -61,52 +62,20 @@ int main(int argc, char *argv[])
 
     while (runTime.run())
     {
-        // --- upwind interpolation of primitive fields on faces
+        // --- Directed interpolation of primitive fields onto faces
 
-        surfaceScalarField rho_pos
+        surfaceScalarField rho_pos(interpolate(rho, pos));
-        (
+        surfaceScalarField rho_neg(interpolate(rho, neg));
-            "rho_pos",
-            fvc::interpolate(rho, pos, "reconstruct(rho)")
-        );
-        surfaceScalarField rho_neg
-        (
-            "rho_neg",
-            fvc::interpolate(rho, neg, "reconstruct(rho)")
-        );
 
-        surfaceVectorField rhoU_pos
+        surfaceVectorField rhoU_pos(interpolate(rhoU, pos, U.name()));
-        (
+        surfaceVectorField rhoU_neg(interpolate(rhoU, neg, U.name()));
-            "rhoU_pos",
-            fvc::interpolate(rhoU, pos, "reconstruct(U)")
-        );
-        surfaceVectorField rhoU_neg
-        (
-            "rhoU_neg",
-            fvc::interpolate(rhoU, neg, "reconstruct(U)")
-        );
 
-        volScalarField rPsi(1.0/psi);
+        volScalarField rPsi("rPsi", 1.0/psi);
-        surfaceScalarField rPsi_pos
+        surfaceScalarField rPsi_pos(interpolate(rPsi, pos, T.name()));
-        (
+        surfaceScalarField rPsi_neg(interpolate(rPsi, neg, T.name()));
-            "rPsi_pos",
-            fvc::interpolate(rPsi, pos, "reconstruct(T)")
-        );
-        surfaceScalarField rPsi_neg
-        (
-            "rPsi_neg",
-            fvc::interpolate(rPsi, neg, "reconstruct(T)")
-        );
 
-        surfaceScalarField e_pos
+        surfaceScalarField e_pos(interpolate(e, pos, T.name()));
-        (
+        surfaceScalarField e_neg(interpolate(e, neg, T.name()));
-            "e_pos",
-            fvc::interpolate(e, pos, "reconstruct(T)")
-        );
-        surfaceScalarField e_neg
-        (
-            "e_neg",
-            fvc::interpolate(e, neg, "reconstruct(T)")
-        );
 
         surfaceVectorField U_pos("U_pos", rhoU_pos/rho_pos);
         surfaceVectorField U_neg("U_neg", rhoU_neg/rho_neg);
@@ -117,16 +86,16 @@ int main(int argc, char *argv[])
         surfaceScalarField phiv_pos("phiv_pos", U_pos & mesh.Sf());
         surfaceScalarField phiv_neg("phiv_neg", U_neg & mesh.Sf());
 
-        volScalarField c(sqrt(thermo.Cp()/thermo.Cv()*rPsi));
+        volScalarField c("c", sqrt(thermo.Cp()/thermo.Cv()*rPsi));
         surfaceScalarField cSf_pos
         (
             "cSf_pos",
-            fvc::interpolate(c, pos, "reconstruct(T)")*mesh.magSf()
+            interpolate(c, pos, T.name())*mesh.magSf()
         );
         surfaceScalarField cSf_neg
         (
             "cSf_neg",
-            fvc::interpolate(c, neg, "reconstruct(T)")*mesh.magSf()
+            interpolate(c, neg, T.name())*mesh.magSf()
         );
 
         surfaceScalarField ap
@@ -188,7 +157,7 @@ int main(int argc, char *argv[])
           + aSf*p_pos - aSf*p_neg
         );
 
-        volScalarField muEff(turbulence->muEff());
+        volScalarField muEff("muEff", turbulence->muEff());
         volTensorField tauMC("tauMC", muEff*dev2(Foam::T(fvc::grad(U))));
 
         // --- Solve density