Completed update of ddtCorr to use Uf for DyM solvers

applications/solvers/combustion/coldEngineFoam/coldEngineFoam.C

@@ -48,6 +48,7 @@ int main(int argc, char *argv[])
     #include "createEngineMesh.H"
     #include "createFields.H"
     #include "createFvOptions.H"
+    #include "createRhoUf.H"
     #include "initContinuityErrs.H"
     #include "readEngineTimeControls.H"
     #include "compressibleCourantNo.H"

applications/solvers/combustion/engineFoam/engineFoam.C

@@ -72,6 +72,7 @@ int main(int argc, char *argv[])
     #include "readCombustionProperties.H"
     #include "createFields.H"
     #include "createFvOptions.H"
+    #include "createRhoUf.H"
     #include "initContinuityErrs.H"
     #include "readEngineTimeControls.H"
     #include "compressibleCourantNo.H"

applications/solvers/combustion/engineFoam/pEqn.H

@@ -15,12 +15,12 @@ if (pimple.transonic())
        *(
             (
                 (fvc::interpolate(rho*HbyA) & mesh.Sf())
-              //***HGW + rhorAUf*fvc::ddtCorr(rho, U, phi)
+              + rhorAUf*fvc::ddtCorr(rho, U, rhoUf)
             )/fvc::interpolate(rho)
-          - fvc::meshPhi(rho, U)
         )
     );
 
+    fvc::makeRelative(phid, psi, U);
     fvOptions.makeRelative(fvc::interpolate(psi), phid);
 
     while (pimple.correctNonOrthogonal())
@@ -51,11 +51,11 @@ else
         "phiHbyA",
         (
             (fvc::interpolate(rho*HbyA) & mesh.Sf())
-          //***HGW + rhorAUf*fvc::ddtCorr(rho, U, phi)
+          + rhorAUf*fvc::ddtCorr(rho, U, rhoUf)
         )
-      - fvc::interpolate(rho)*fvc::meshPhi(rho, U)
     );
 
+    fvc::makeRelative(phiHbyA, rho, U);
     fvOptions.makeRelative(fvc::interpolate(rho), phiHbyA);
 
     while (pimple.correctNonOrthogonal())
@@ -88,6 +88,15 @@ U.correctBoundaryConditions();
 fvOptions.correct(U);
 K = 0.5*magSqr(U);
 
+{
+    rhoUf = fvc::interpolate(rho*U);
+    surfaceVectorField n(mesh.Sf()/mesh.magSf());
+    rhoUf +=
+        mesh.Sf()
+       *(fvc::absolute(phi, rho, U) - (mesh.Sf() & rhoUf))
+       /sqr(mesh.magSf());
+}
+
 if (thermo.dpdt())
 {
     dpdt = fvc::ddt(p) - fvc::div(fvc::meshPhi(rho, U), p);

applications/solvers/compressible/rhoPimpleFoam/rhoPimpleDyMFoam/pEqn.H

@@ -26,6 +26,7 @@ if (pimple.transonic())
         )/fvc::interpolate(rho)
     );
 
+    fvc::makeRelative(phid, psi, U);
     fvOptions.makeRelative(fvc::interpolate(psi), phid);
 
     while (pimple.correctNonOrthogonal())
@@ -58,6 +59,7 @@ else
       + rhorAUf*fvc::ddtCorr(rho, U, rhoUf)
     );
 
+    fvc::makeRelative(phiHbyA, rho, U);
     fvOptions.makeRelative(fvc::interpolate(rho), phiHbyA);
 
     while (pimple.correctNonOrthogonal())
@@ -102,6 +104,15 @@ U.correctBoundaryConditions();
 fvOptions.correct(U);
 K = 0.5*magSqr(U);
 
+{
+    rhoUf = fvc::interpolate(rho*U);
+    surfaceVectorField n(mesh.Sf()/mesh.magSf());
+    rhoUf +=
+        mesh.Sf()
+       *(fvc::absolute(phi, rho, U) - (mesh.Sf() & rhoUf))
+       /sqr(mesh.magSf());
+}
+
 if (thermo.dpdt())
 {
     dpdt = fvc::ddt(p);

applications/solvers/compressible/rhoPimpleFoam/rhoPimpleDyMFoam/rhoPimpleDyMFoam.C

@@ -76,6 +76,9 @@ int main(int argc, char *argv[])
         Info<< "Time = " << runTime.timeName() << nl << endl;
 
         {
+            // Store momentum to set rhoUf for introduced faces.
+            volVectorField rhoU("rhoU", rho*U);
+
             // Store divrhoU from the previous time-step/mesh for the correctPhi
             volScalarField divrhoU(fvc::div(fvc::absolute(phi, rho, U)));
 

applications/solvers/compressible/sonicFoam/sonicDyMFoam/pEqn.H

@@ -13,12 +13,13 @@ surfaceScalarField phid
    *(
         (
             (fvc::interpolate(rho*HbyA) & mesh.Sf())
-          //***HGW + rhorAUf*fvc::ddtCorr(rho, U, phi)
+          + rhorAUf*fvc::ddtCorr(rho, U, rhoUf)
         )/fvc::interpolate(rho)
-      - fvc::meshPhi(rho, U)
     )
 );
 
+fvc::makeRelative(phid, psi, U);
+
 for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
 {
     fvScalarMatrix pEqn
@@ -37,3 +38,12 @@ for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
 
 U = HbyA - rAU*fvc::grad(p);
 U.correctBoundaryConditions();
+
+{
+    rhoUf = fvc::interpolate(rho*U);
+    surfaceVectorField n(mesh.Sf()/mesh.magSf());
+    rhoUf +=
+        mesh.Sf()
+       *(fvc::absolute(phi, rho, U) - (mesh.Sf() & rhoUf))
+       /sqr(mesh.magSf());
+}

applications/solvers/compressible/sonicFoam/sonicDyMFoam/sonicDyMFoam.C

@@ -64,12 +64,6 @@ int main(int argc, char *argv[])
 
         mesh.movePoints(motionPtr->newPoints());
 
-        // Calculate absolute flux from the mapped surface velocity
-        phi = mesh.Sf() & rhoUf;
-
-        // Make the flux relative to the mesh motion
-        fvc::makeRelative(phi, rho, U);
-
         #include "rhoEqn.H"
 
         // --- Pressure-velocity PIMPLE corrector loop

applications/solvers/lagrangian/sprayFoam/sprayEngineFoam/pEqn.H

@@ -15,12 +15,12 @@ if (pimple.transonic())
        *(
             (
                 (fvc::interpolate(rho*HbyA) & mesh.Sf())
-              //***HGW + rhorAUf*fvc::ddtCorr(rho, U, phi)
+              + rhorAUf*fvc::ddtCorr(rho, U, rhoUf)
             )/fvc::interpolate(rho)
-          - fvc::meshPhi(rho, U)
         )
     );
 
+    fvc::makeRelative(phid, psi, U);
     fvOptions.makeRelative(fvc::interpolate(psi), phid);
 
     while (pimple.correctNonOrthogonal())
@@ -52,11 +52,11 @@ else
         "phiHbyA",
         (
             (fvc::interpolate(HbyA) & mesh.Sf())
-          //***HGW + rhorAUf*fvc::ddtCorr(rho, U, phi)
+          + rhorAUf*fvc::ddtCorr(rho, U, rhoUf)
         )
-      - fvc::interpolate(rho)*fvc::meshPhi(rho, U)
     );
 
+    fvc::makeRelative(phiHbyA, rho, U);
     fvOptions.makeRelative(fvc::interpolate(rho), phiHbyA);
 
     while (pimple.correctNonOrthogonal())
@@ -90,6 +90,15 @@ U.correctBoundaryConditions();
 fvOptions.correct(U);
 K = 0.5*magSqr(U);
 
+{
+    rhoUf = fvc::interpolate(rho*U);
+    surfaceVectorField n(mesh.Sf()/mesh.magSf());
+    rhoUf +=
+        mesh.Sf()
+       *(fvc::absolute(phi, rho, U) - (mesh.Sf() & rhoUf))
+       /sqr(mesh.magSf());
+}
+
 if (thermo.dpdt())
 {
     dpdt = fvc::ddt(p) - fvc::div(fvc::meshPhi(rho, U), p);

applications/solvers/lagrangian/sprayFoam/sprayEngineFoam/sprayEngineFoam.C

@@ -52,6 +52,7 @@ int main(int argc, char *argv[])
     #include "readGravitationalAcceleration.H"
     #include "createFields.H"
     #include "createFvOptions.H"
+    #include "createRhoUf.H"
     #include "createClouds.H"
     #include "createRadiationModel.H"
     #include "initContinuityErrs.H"