diff --git a/applications/solvers/compressible/rhoCentralFoam/createFields.H b/applications/solvers/compressible/rhoCentralFoam/createFields.H
index 1622dd3501d6230fa52be76444c6026848c06f15..c1618d6346a2b8c02dda50ec480b560ab4f043f5 100644
--- a/applications/solvers/compressible/rhoCentralFoam/createFields.H
+++ b/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();
diff --git a/applications/solvers/compressible/rhoCentralFoam/directionInterpolate.H b/applications/solvers/compressible/rhoCentralFoam/directionInterpolate.H
new file mode 100644
index 0000000000000000000000000000000000000000..443bab5f1f0456208692ae94e11b964da0baf09f
--- /dev/null
+++ b/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;
+}
+
+}
diff --git a/applications/solvers/compressible/rhoCentralFoam/rhoCentralDyMFoam/rhoCentralDyMFoam.C b/applications/solvers/compressible/rhoCentralFoam/rhoCentralDyMFoam/rhoCentralDyMFoam.C
index 4f8f95d11956f275f703a50a3c72e70dabd57a6e..4416a7cafafa1a651dec9fe2b5c8e963caf6d2e5 100644
--- a/applications/solvers/compressible/rhoCentralFoam/rhoCentralDyMFoam/rhoCentralDyMFoam.C
+++ b/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
- (
- "rho_pos",
- fvc::interpolate(rho, pos, "reconstruct(rho)")
- );
- surfaceScalarField rho_neg
- (
- "rho_neg",
- fvc::interpolate(rho, neg, "reconstruct(rho)")
- );
+ surfaceScalarField rho_pos(interpolate(rho, pos));
+ surfaceScalarField rho_neg(interpolate(rho, neg));
- surfaceVectorField rhoU_pos
- (
- "rhoU_pos",
- fvc::interpolate(rhoU, pos, "reconstruct(U)")
- );
- surfaceVectorField rhoU_neg
- (
- "rhoU_neg",
- fvc::interpolate(rhoU, neg, "reconstruct(U)")
- );
+ surfaceVectorField rhoU_pos(interpolate(rhoU, pos, U.name()));
+ surfaceVectorField rhoU_neg(interpolate(rhoU, neg, U.name()));
- volScalarField rPsi(1.0/psi);
- surfaceScalarField rPsi_pos
- (
- "rPsi_pos",
- fvc::interpolate(rPsi, pos, "reconstruct(T)")
- );
- surfaceScalarField rPsi_neg
- (
- "rPsi_neg",
- fvc::interpolate(rPsi, neg, "reconstruct(T)")
- );
+ volScalarField rPsi("rPsi", 1.0/psi);
+ surfaceScalarField rPsi_pos(interpolate(rPsi, pos, T.name()));
+ surfaceScalarField rPsi_neg(interpolate(rPsi, neg, T.name()));
- surfaceScalarField e_pos
- (
- "e_pos",
- fvc::interpolate(e, pos, "reconstruct(T)")
- );
- surfaceScalarField e_neg
- (
- "e_neg",
- fvc::interpolate(e, neg, "reconstruct(T)")
- );
+ surfaceScalarField e_pos(interpolate(e, pos, T.name()));
+ surfaceScalarField e_neg(interpolate(e, neg, T.name()));
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
diff --git a/applications/solvers/compressible/rhoCentralFoam/rhoCentralFoam.C b/applications/solvers/compressible/rhoCentralFoam/rhoCentralFoam.C
index 2ec5f476171c2c62ca08c44d661ba8b2bd572cb1..93a88c783b2297534f791582c4d0681b4fd13c9f 100644
--- a/applications/solvers/compressible/rhoCentralFoam/rhoCentralFoam.C
+++ b/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
- (
- "rho_pos",
- fvc::interpolate(rho, pos, "reconstruct(rho)")
- );
- surfaceScalarField rho_neg
- (
- "rho_neg",
- fvc::interpolate(rho, neg, "reconstruct(rho)")
- );
+ surfaceScalarField rho_pos(interpolate(rho, pos));
+ surfaceScalarField rho_neg(interpolate(rho, neg));
- surfaceVectorField rhoU_pos
- (
- "rhoU_pos",
- fvc::interpolate(rhoU, pos, "reconstruct(U)")
- );
- surfaceVectorField rhoU_neg
- (
- "rhoU_neg",
- fvc::interpolate(rhoU, neg, "reconstruct(U)")
- );
+ surfaceVectorField rhoU_pos(interpolate(rhoU, pos, U.name()));
+ surfaceVectorField rhoU_neg(interpolate(rhoU, neg, U.name()));
- volScalarField rPsi(1.0/psi);
- surfaceScalarField rPsi_pos
- (
- "rPsi_pos",
- fvc::interpolate(rPsi, pos, "reconstruct(T)")
- );
- surfaceScalarField rPsi_neg
- (
- "rPsi_neg",
- fvc::interpolate(rPsi, neg, "reconstruct(T)")
- );
+ volScalarField rPsi("rPsi", 1.0/psi);
+ surfaceScalarField rPsi_pos(interpolate(rPsi, pos, T.name()));
+ surfaceScalarField rPsi_neg(interpolate(rPsi, neg, T.name()));
- surfaceScalarField e_pos
- (
- "e_pos",
- fvc::interpolate(e, pos, "reconstruct(T)")
- );
- surfaceScalarField e_neg
- (
- "e_neg",
- fvc::interpolate(e, neg, "reconstruct(T)")
- );
+ surfaceScalarField e_pos(interpolate(e, pos, T.name()));
+ surfaceScalarField e_neg(interpolate(e, neg, T.name()));
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