From 17f2d5f189b2562ec2acbcb90c91d14b296a3251 Mon Sep 17 00:00:00 2001
From: Henry <Henry>
Date: Thu, 12 Sep 2013 16:58:39 +0100
Subject: [PATCH] Completed update of ddtCorr to use Uf for DyM solvers

---
 .../combustion/coldEngineFoam/coldEngineFoam.C  |  1 +
 .../solvers/combustion/engineFoam/engineFoam.C  |  1 +
 .../solvers/combustion/engineFoam/pEqn.H        | 17 +++++++++++++----
 .../rhoPimpleFoam/rhoPimpleDyMFoam/pEqn.H       | 11 +++++++++++
 .../rhoPimpleDyMFoam/rhoPimpleDyMFoam.C         |  3 +++
 .../compressible/sonicFoam/sonicDyMFoam/pEqn.H  | 14 ++++++++++++--
 .../sonicFoam/sonicDyMFoam/sonicDyMFoam.C       |  6 ------
 .../lagrangian/sprayFoam/sprayEngineFoam/pEqn.H | 17 +++++++++++++----
 .../sprayFoam/sprayEngineFoam/sprayEngineFoam.C |  1 +
 9 files changed, 55 insertions(+), 16 deletions(-)

diff --git a/applications/solvers/combustion/coldEngineFoam/coldEngineFoam.C b/applications/solvers/combustion/coldEngineFoam/coldEngineFoam.C
index bd5e2c3b99d..43eb616e173 100644
--- a/applications/solvers/combustion/coldEngineFoam/coldEngineFoam.C
+++ b/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"
diff --git a/applications/solvers/combustion/engineFoam/engineFoam.C b/applications/solvers/combustion/engineFoam/engineFoam.C
index 592abbd3e67..a35e22dc58d 100644
--- a/applications/solvers/combustion/engineFoam/engineFoam.C
+++ b/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"
diff --git a/applications/solvers/combustion/engineFoam/pEqn.H b/applications/solvers/combustion/engineFoam/pEqn.H
index 384ecbee598..2f6333cc976 100644
--- a/applications/solvers/combustion/engineFoam/pEqn.H
+++ b/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);
diff --git a/applications/solvers/compressible/rhoPimpleFoam/rhoPimpleDyMFoam/pEqn.H b/applications/solvers/compressible/rhoPimpleFoam/rhoPimpleDyMFoam/pEqn.H
index 40cf1094579..14ba97d1c2d 100644
--- a/applications/solvers/compressible/rhoPimpleFoam/rhoPimpleDyMFoam/pEqn.H
+++ b/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);
diff --git a/applications/solvers/compressible/rhoPimpleFoam/rhoPimpleDyMFoam/rhoPimpleDyMFoam.C b/applications/solvers/compressible/rhoPimpleFoam/rhoPimpleDyMFoam/rhoPimpleDyMFoam.C
index aae32851863..e633ebee602 100644
--- a/applications/solvers/compressible/rhoPimpleFoam/rhoPimpleDyMFoam/rhoPimpleDyMFoam.C
+++ b/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)));
 
diff --git a/applications/solvers/compressible/sonicFoam/sonicDyMFoam/pEqn.H b/applications/solvers/compressible/sonicFoam/sonicDyMFoam/pEqn.H
index f8f99dafa54..65afd1aeb88 100644
--- a/applications/solvers/compressible/sonicFoam/sonicDyMFoam/pEqn.H
+++ b/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());
+}
diff --git a/applications/solvers/compressible/sonicFoam/sonicDyMFoam/sonicDyMFoam.C b/applications/solvers/compressible/sonicFoam/sonicDyMFoam/sonicDyMFoam.C
index 11ec1de915e..62d4716fb5a 100644
--- a/applications/solvers/compressible/sonicFoam/sonicDyMFoam/sonicDyMFoam.C
+++ b/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
diff --git a/applications/solvers/lagrangian/sprayFoam/sprayEngineFoam/pEqn.H b/applications/solvers/lagrangian/sprayFoam/sprayEngineFoam/pEqn.H
index 62c3e0376e3..b10c24aa0b7 100644
--- a/applications/solvers/lagrangian/sprayFoam/sprayEngineFoam/pEqn.H
+++ b/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);
diff --git a/applications/solvers/lagrangian/sprayFoam/sprayEngineFoam/sprayEngineFoam.C b/applications/solvers/lagrangian/sprayFoam/sprayEngineFoam/sprayEngineFoam.C
index baafc8d5e3f..b91a3ad79de 100644
--- a/applications/solvers/lagrangian/sprayFoam/sprayEngineFoam/sprayEngineFoam.C
+++ b/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"
-- 
GitLab