diff --git a/applications/solvers/DNS/dnsFoam/dnsFoam.C b/applications/solvers/DNS/dnsFoam/dnsFoam.C
index 08073398a7..d5a804841f 100644
--- a/applications/solvers/DNS/dnsFoam/dnsFoam.C
+++ b/applications/solvers/DNS/dnsFoam/dnsFoam.C
@@ -85,7 +85,8 @@ int main(int argc, char *argv[])
for (int corr=1; corr<=1; corr++)
{
- volScalarField rAU("Dp", 1.0/UEqn.A());
+ volScalarField rAU(1.0/UEqn.A());
+ surfaceScalarField Dp("Dp", fvc::interpolate(rAU));
volVectorField HbyA("HbyA", U);
HbyA = rAU*UEqn.H();
@@ -93,12 +94,12 @@ int main(int argc, char *argv[])
(
"phiHbyA",
(fvc::interpolate(HbyA) & mesh.Sf())
- + fvc::ddtPhiCorr(rAU, U, phi)
+ + Dp*fvc::ddtCorr(U, phi)
);
fvScalarMatrix pEqn
(
- fvm::laplacian(rAU, p) == fvc::div(phiHbyA)
+ fvm::laplacian(Dp, p) == fvc::div(phiHbyA)
);
pEqn.solve();
diff --git a/applications/solvers/combustion/PDRFoam/pEqn.H b/applications/solvers/combustion/PDRFoam/pEqn.H
index e1c4274f33..3053d4c68c 100644
--- a/applications/solvers/combustion/PDRFoam/pEqn.H
+++ b/applications/solvers/combustion/PDRFoam/pEqn.H
@@ -1,6 +1,7 @@
rho = thermo.rho();
volScalarField rAU(1.0/UEqn.A());
+
volVectorField HbyA("HbyA", U);
HbyA = invA & UEqn.H();
@@ -11,9 +12,9 @@ if (pimple.transonic())
"phid",
fvc::interpolate(psi)
*(
- (fvc::interpolate(HbyA) & mesh.Sf())
- + fvc::ddtPhiCorr(rAU, rho, U, phi)
- )
+ (fvc::interpolate(rho*HbyA) & mesh.Sf())
+ + fvc::interpolate(rho*rAU)*fvc::ddtCorr(rho, U, phi)
+ )/fvc::interpolate(rho)
);
while (pimple.correctNonOrthogonal())
@@ -38,10 +39,9 @@ else
surfaceScalarField phiHbyA
(
"phiHbyA",
- fvc::interpolate(rho)*
(
- (fvc::interpolate(HbyA) & mesh.Sf())
- + fvc::ddtPhiCorr(rAU, rho, U, phi)
+ (fvc::interpolate(rho*HbyA) & mesh.Sf())
+ + fvc::interpolate(rho*rAU)*fvc::ddtCorr(rho, U, phi)
)
);
diff --git a/applications/solvers/combustion/XiFoam/pEqn.H b/applications/solvers/combustion/XiFoam/pEqn.H
index c1effcd60d..83a44ab7c2 100644
--- a/applications/solvers/combustion/XiFoam/pEqn.H
+++ b/applications/solvers/combustion/XiFoam/pEqn.H
@@ -1,6 +1,8 @@
rho = thermo.rho();
volScalarField rAU(1.0/UEqn.A());
+surfaceScalarField Dp("Dp", fvc::interpolate(rho*rAU));
+
volVectorField HbyA("HbyA", U);
HbyA = rAU*UEqn.H();
@@ -11,9 +13,9 @@ if (pimple.transonic())
"phid",
fvc::interpolate(psi)
*(
- (fvc::interpolate(HbyA) & mesh.Sf())
- + fvc::ddtPhiCorr(rAU, rho, U, phi)
- )
+ (fvc::interpolate(rho*HbyA) & mesh.Sf())
+ + Dp*fvc::ddtCorr(rho, U, phi)
+ )/fvc::interpolate(rho)
);
fvOptions.makeRelative(fvc::interpolate(psi), phid);
@@ -24,7 +26,7 @@ if (pimple.transonic())
(
fvm::ddt(psi, p)
+ fvm::div(phid, p)
- - fvm::laplacian(rho*rAU, p)
+ - fvm::laplacian(Dp, p)
==
fvOptions(psi, p, rho.name())
);
@@ -44,10 +46,9 @@ else
surfaceScalarField phiHbyA
(
"phiHbyA",
- fvc::interpolate(rho)
- *(
- (fvc::interpolate(HbyA) & mesh.Sf())
- + fvc::ddtPhiCorr(rAU, rho, U, phi)
+ (
+ (fvc::interpolate(rho*HbyA) & mesh.Sf())
+ + Dp*fvc::ddtCorr(rho, U, phi)
)
);
@@ -59,7 +60,7 @@ else
(
fvm::ddt(psi, p)
+ fvc::div(phiHbyA)
- - fvm::laplacian(rho*rAU, p)
+ - fvm::laplacian(Dp, p)
==
fvOptions(psi, p, rho.name())
);
diff --git a/applications/solvers/combustion/engineFoam/pEqn.H b/applications/solvers/combustion/engineFoam/pEqn.H
index 39939333c2..a92e7285ab 100644
--- a/applications/solvers/combustion/engineFoam/pEqn.H
+++ b/applications/solvers/combustion/engineFoam/pEqn.H
@@ -1,6 +1,8 @@
rho = thermo.rho();
volScalarField rAU(1.0/UEqn.A());
+surfaceScalarField Dp("Dp", fvc::interpolate(rho*rAU));
+
volVectorField HbyA("HbyA", U);
HbyA = rAU*UEqn.H();
@@ -10,7 +12,13 @@ if (pimple.transonic())
(
"phid",
fvc::interpolate(psi)
- *((fvc::interpolate(HbyA) & mesh.Sf()) - fvc::meshPhi(rho, U))
+ *(
+ (
+ (fvc::interpolate(rho*HbyA) & mesh.Sf())
+ //***HGW + Dp*fvc::ddtCorr(rho, U, phi)
+ )/fvc::interpolate(rho)
+ - fvc::meshPhi(rho, U)
+ )
);
fvOptions.makeRelative(fvc::interpolate(psi), phid);
@@ -41,8 +49,11 @@ else
surfaceScalarField phiHbyA
(
"phiHbyA",
- fvc::interpolate(rho)
- *((fvc::interpolate(HbyA) & mesh.Sf()) - fvc::meshPhi(rho, U))
+ (
+ (fvc::interpolate(rho*HbyA) & mesh.Sf())
+ //***HGW + Dp*fvc::ddtCorr(rho, U, phi)
+ )
+ - fvc::interpolate(rho)*fvc::meshPhi(rho, U)
);
fvOptions.makeRelative(fvc::interpolate(rho), phiHbyA);
diff --git a/applications/solvers/combustion/fireFoam/pEqn.H b/applications/solvers/combustion/fireFoam/pEqn.H
index 9953e83007..c98ac61f3b 100644
--- a/applications/solvers/combustion/fireFoam/pEqn.H
+++ b/applications/solvers/combustion/fireFoam/pEqn.H
@@ -1,22 +1,22 @@
rho = thermo.rho();
volScalarField rAU(1.0/UEqn.A());
-surfaceScalarField rhorAUf("Dp", fvc::interpolate(rho*rAU));
+surfaceScalarField Dp("Dp", fvc::interpolate(rho*rAU));
+
volVectorField HbyA("HbyA", U);
HbyA = rAU*UEqn.H();
phi.boundaryField() =
fvc::interpolate(rho.boundaryField()*U.boundaryField())
& mesh.Sf().boundaryField();
-surfaceScalarField phig(-rhorAUf*ghf*fvc::snGrad(rho)*mesh.magSf());
+surfaceScalarField phig(-Dp*ghf*fvc::snGrad(rho)*mesh.magSf());
surfaceScalarField phiHbyA
(
"phiHbyA",
- fvc::interpolate(rho)
- *(
- (fvc::interpolate(HbyA) & mesh.Sf())
- + fvc::ddtPhiCorr(rAU, rho, U, phi)
+ (
+ (fvc::interpolate(rho*HbyA) & mesh.Sf())
+ + Dp*fvc::ddtCorr(rho, U, phi)
)
+ phig
);
@@ -30,7 +30,7 @@ while (pimple.correctNonOrthogonal())
fvc::ddt(psi, rho)*gh
+ fvc::div(phiHbyA)
+ fvm::ddt(psi, p_rgh)
- - fvm::laplacian(rhorAUf, p_rgh)
+ - fvm::laplacian(Dp, p_rgh)
==
parcels.Srho()
+ surfaceFilm.Srho()
@@ -44,7 +44,7 @@ while (pimple.correctNonOrthogonal())
if (pimple.finalNonOrthogonalIter())
{
phi = phiHbyA + p_rghEqn.flux();
- U = HbyA + rAU*fvc::reconstruct((p_rghEqn.flux() + phig)/rhorAUf);
+ U = HbyA + rAU*fvc::reconstruct((p_rghEqn.flux() + phig)/Dp);
U.correctBoundaryConditions();
fvOptions.correct(U);
}
diff --git a/applications/solvers/combustion/reactingFoam/pEqn.H b/applications/solvers/combustion/reactingFoam/pEqn.H
index 0a8a6ce9d0..0dc5d422db 100644
--- a/applications/solvers/combustion/reactingFoam/pEqn.H
+++ b/applications/solvers/combustion/reactingFoam/pEqn.H
@@ -1,6 +1,8 @@
rho = thermo.rho();
volScalarField rAU(1.0/UEqn.A());
+surfaceScalarField Dp("Dp", fvc::interpolate(rho*rAU));
+
volVectorField HbyA("HbyA", U);
HbyA = rAU*UEqn.H();
@@ -11,9 +13,9 @@ if (pimple.transonic())
"phid",
fvc::interpolate(psi)
*(
- (fvc::interpolate(HbyA) & mesh.Sf())
- + fvc::ddtPhiCorr(rAU, rho, U, phi)
- )
+ (fvc::interpolate(rho*HbyA) & mesh.Sf())
+ + Dp*fvc::ddtCorr(rho, U, phi)
+ )/fvc::interpolate(rho)
);
fvOptions.makeRelative(fvc::interpolate(psi), phid);
@@ -44,10 +46,9 @@ else
surfaceScalarField phiHbyA
(
"phiHbyA",
- fvc::interpolate(rho)
- *(
- (fvc::interpolate(HbyA) & mesh.Sf())
- + fvc::ddtPhiCorr(rAU, rho, U, phi)
+ (
+ (fvc::interpolate(rho*HbyA) & mesh.Sf())
+ + Dp*fvc::ddtCorr(rho, U, phi)
)
);
diff --git a/applications/solvers/combustion/reactingFoam/rhoReactingBuoyantFoam/pEqn.H b/applications/solvers/combustion/reactingFoam/rhoReactingBuoyantFoam/pEqn.H
index 7f312e5030..8602725bdb 100644
--- a/applications/solvers/combustion/reactingFoam/rhoReactingBuoyantFoam/pEqn.H
+++ b/applications/solvers/combustion/reactingFoam/rhoReactingBuoyantFoam/pEqn.H
@@ -6,20 +6,19 @@
thermo.rho() -= psi*p;
volScalarField rAU(1.0/UEqn.A());
- surfaceScalarField rhorAUf("Dp", fvc::interpolate(rho*rAU));
+ surfaceScalarField Dp("Dp", fvc::interpolate(rho*rAU));
volVectorField HbyA("HbyA", U);
HbyA = rAU*UEqn.H();
- surfaceScalarField phig(-rhorAUf*ghf*fvc::snGrad(rho)*mesh.magSf());
+ surfaceScalarField phig(-Dp*ghf*fvc::snGrad(rho)*mesh.magSf());
surfaceScalarField phiHbyA
(
"phiHbyA",
- fvc::interpolate(rho)
- *(
- (fvc::interpolate(U) & mesh.Sf())
- + fvc::ddtPhiCorr(rAU, rho, U, phi)
+ (
+ (fvc::interpolate(rho*HbyA) & mesh.Sf())
+ + Dp*fvc::ddtCorr(rho, U, phi)
)
+ phig
);
@@ -39,7 +38,7 @@
fvScalarMatrix p_rghEqn
(
p_rghDDtEqn
- - fvm::laplacian(rhorAUf, p_rgh)
+ - fvm::laplacian(Dp, p_rgh)
);
fvOptions.constrain(p_rghEqn);
@@ -56,7 +55,7 @@
// Correct the momentum source with the pressure gradient flux
// calculated from the relaxed pressure
- U = HbyA + rAU*fvc::reconstruct((phig + p_rghEqn.flux())/rhorAUf);
+ U = HbyA + rAU*fvc::reconstruct((phig + p_rghEqn.flux())/Dp);
U.correctBoundaryConditions();
fvOptions.correct(U);
K = 0.5*magSqr(U);
diff --git a/applications/solvers/combustion/reactingFoam/rhoReactingFoam/pEqn.H b/applications/solvers/combustion/reactingFoam/rhoReactingFoam/pEqn.H
index 461e1ca8fb..3c92ce8119 100644
--- a/applications/solvers/combustion/reactingFoam/rhoReactingFoam/pEqn.H
+++ b/applications/solvers/combustion/reactingFoam/rhoReactingFoam/pEqn.H
@@ -6,6 +6,8 @@
thermo.rho() -= psi*p;
volScalarField rAU(1.0/UEqn.A());
+ surfaceScalarField Dp("Dp", fvc::interpolate(rho*rAU));
+
volVectorField HbyA("HbyA", U);
HbyA = rAU*UEqn.H();
@@ -14,8 +16,10 @@
surfaceScalarField phiHbyA
(
"phiHbyA",
- (fvc::interpolate(HbyA) & mesh.Sf())
- + fvc::ddtPhiCorr(rAU, rho, U, phi)
+ (
+ (fvc::interpolate(rho*HbyA) & mesh.Sf())
+ + Dp*fvc::ddtCorr(rho, U, phi)
+ )/fvc::interpolate(rho)
);
fvOptions.makeRelative(phiHbyA);
@@ -35,7 +39,7 @@
fvScalarMatrix pEqn
(
pDDtEqn
- - fvm::laplacian(rho*rAU, p)
+ - fvm::laplacian(Dp, p)
==
fvOptions(psi, p, rho.name())
);
@@ -55,10 +59,9 @@
surfaceScalarField phiHbyA
(
"phiHbyA",
- fvc::interpolate(rho)
- *(
- (fvc::interpolate(HbyA) & mesh.Sf())
- + fvc::ddtPhiCorr(rAU, rho, U, phi)
+ (
+ (fvc::interpolate(rho*HbyA) & mesh.Sf())
+ + Dp*fvc::ddtCorr(rho, U, phi)
)
);
@@ -77,7 +80,7 @@
fvScalarMatrix pEqn
(
pDDtEqn
- - fvm::laplacian(rho*rAU, p)
+ - fvm::laplacian(Dp, p)
);
fvOptions.constrain(pEqn);
diff --git a/applications/solvers/compressible/rhoPimpleFoam/pEqn.H b/applications/solvers/compressible/rhoPimpleFoam/pEqn.H
index fcc40881a1..373bafcd3b 100644
--- a/applications/solvers/compressible/rhoPimpleFoam/pEqn.H
+++ b/applications/solvers/compressible/rhoPimpleFoam/pEqn.H
@@ -4,6 +4,8 @@ rho = min(rho, rhoMax);
rho.relax();
volScalarField rAU(1.0/UEqn().A());
+surfaceScalarField Dp("Dp", fvc::interpolate(rho*rAU));
+
volVectorField HbyA("HbyA", U);
HbyA = rAU*UEqn().H();
@@ -19,15 +21,13 @@ if (pimple.transonic())
"phid",
fvc::interpolate(psi)
*(
- (fvc::interpolate(HbyA) & mesh.Sf())
- + fvc::ddtPhiCorr(rAU, rho, U, phi)
- )
+ (fvc::interpolate(rho*HbyA) & mesh.Sf())
+ + Dp*fvc::ddtCorr(rho, U, phi)
+ )/fvc::interpolate(rho)
);
fvOptions.makeRelative(fvc::interpolate(psi), phid);
- volScalarField Dp("Dp", rho*rAU);
-
while (pimple.correctNonOrthogonal())
{
fvScalarMatrix pEqn
@@ -54,10 +54,9 @@ else
surfaceScalarField phiHbyA
(
"phiHbyA",
- fvc::interpolate(rho)
- *(
- (fvc::interpolate(HbyA) & mesh.Sf())
- + fvc::ddtPhiCorr(rAU, rho, U, phi)
+ (
+ (fvc::interpolate(rho*HbyA) & mesh.Sf())
+ + Dp*fvc::ddtCorr(rho, U, phi)
)
);
diff --git a/applications/solvers/compressible/rhoPimpleFoam/rhoPimpleDyMFoam/pEqn.H b/applications/solvers/compressible/rhoPimpleFoam/rhoPimpleDyMFoam/pEqn.H
index c1ceb88c17..d2681663e6 100644
--- a/applications/solvers/compressible/rhoPimpleFoam/rhoPimpleDyMFoam/pEqn.H
+++ b/applications/solvers/compressible/rhoPimpleFoam/rhoPimpleDyMFoam/pEqn.H
@@ -4,6 +4,8 @@ rho = min(rho, rhoMax);
rho.relax();
volScalarField rAU(1.0/UEqn().A());
+surfaceScalarField Dp("Dp", fvc::interpolate(rho*rAU));
+
volVectorField HbyA("HbyA", U);
HbyA = rAU*UEqn().H();
@@ -19,15 +21,13 @@ if (pimple.transonic())
"phid",
fvc::interpolate(psi)
*(
- (fvc::interpolate(HbyA) & mesh.Sf())
- + fvc::ddtPhiCorr(rAU, rho, U, phiAbs)
- )
+ (fvc::interpolate(rho*HbyA) & mesh.Sf())
+ + Dp*fvc::ddtCorr(rho, U, phiAbs)
+ )/fvc::interpolate(rho)
);
fvOptions.makeRelative(fvc::interpolate(psi), phid);
- volScalarField Dp("Dp", rho*rAU);
-
while (pimple.correctNonOrthogonal())
{
fvScalarMatrix pEqn
@@ -54,18 +54,12 @@ else
surfaceScalarField phiHbyA
(
"phiHbyA",
- fvc::interpolate(rho)
- *(
- (fvc::interpolate(HbyA) & mesh.Sf())
- - fvc::meshPhi(rho, U)
- + fvc::ddtPhiCorr(rAU, rho, U, phiAbs)
- )
+ (fvc::interpolate(rho*HbyA) & mesh.Sf())
+ + Dp*fvc::ddtCorr(rho, U, phiAbs)
);
fvOptions.makeRelative(fvc::interpolate(rho), phiHbyA);
- volScalarField Dp("Dp", rho*rAU);
-
while (pimple.correctNonOrthogonal())
{
// Pressure corrector
diff --git a/applications/solvers/compressible/rhoPimpleFoam/rhoPimpleDyMFoam/rhoPimpleDyMFoam.C b/applications/solvers/compressible/rhoPimpleFoam/rhoPimpleDyMFoam/rhoPimpleDyMFoam.C
index fb74f2b39a..dac6980917 100644
--- a/applications/solvers/compressible/rhoPimpleFoam/rhoPimpleDyMFoam/rhoPimpleDyMFoam.C
+++ b/applications/solvers/compressible/rhoPimpleFoam/rhoPimpleDyMFoam/rhoPimpleDyMFoam.C
@@ -59,7 +59,7 @@ int main(int argc, char *argv[])
#include "CourantNo.H"
#include "setInitialDeltaT.H"
- // Create old-time absolute flux for ddtPhiCorr
+ // Create old-time absolute flux for ddtCorr
surfaceScalarField phiAbs("phiAbs", phi);
@@ -75,7 +75,7 @@ int main(int argc, char *argv[])
// Make the fluxes absolute before mesh-motion
fvc::makeAbsolute(phi, rho, U);
- // Update absolute flux for ddtPhiCorr
+ // Update absolute flux for ddtCorr
phiAbs = phi;
#include "setDeltaT.H"
diff --git a/applications/solvers/compressible/rhoPimpleFoam/rhoPimplecFoam/pEqn.H b/applications/solvers/compressible/rhoPimpleFoam/rhoPimplecFoam/pEqn.H
index de09cb53e2..f8837f0aa4 100644
--- a/applications/solvers/compressible/rhoPimpleFoam/rhoPimplecFoam/pEqn.H
+++ b/applications/solvers/compressible/rhoPimpleFoam/rhoPimplecFoam/pEqn.H
@@ -20,9 +20,9 @@ if (pimple.transonic())
"phid",
fvc::interpolate(psi)
*(
- (fvc::interpolate(HbyA) & mesh.Sf())
- + fvc::ddtPhiCorr(rAU, rho, U, phi)
- )
+ (fvc::interpolate(rho*HbyA) & mesh.Sf())
+ + fvc::interpolate(rho*rAU)*fvc::ddtCorr(rho, U, phi)
+ )/fvc::interpolate(rho)
);
fvOptions.makeRelative(fvc::interpolate(psi), phid);
@@ -64,10 +64,9 @@ else
surfaceScalarField phiHbyA
(
"phiHbyA",
- fvc::interpolate(rho)
- *(
- (fvc::interpolate(HbyA) & mesh.Sf())
- + fvc::ddtPhiCorr(rAU, rho, U, phi)
+ (
+ (fvc::interpolate(rho*HbyA) & mesh.Sf())
+ + fvc::interpolate(rho*rAU)*fvc::ddtCorr(rho, U, phi)
)
);
diff --git a/applications/solvers/compressible/rhoSimpleFoam/pEqn.H b/applications/solvers/compressible/rhoSimpleFoam/pEqn.H
index 364d26a50b..444f5e996a 100644
--- a/applications/solvers/compressible/rhoSimpleFoam/pEqn.H
+++ b/applications/solvers/compressible/rhoSimpleFoam/pEqn.H
@@ -12,7 +12,9 @@
surfaceScalarField phid
(
"phid",
- fvc::interpolate(psi)*(fvc::interpolate(HbyA) & mesh.Sf())
+ fvc::interpolate(psi)
+ *(fvc::interpolate(rho*HbyA) & mesh.Sf())
+ /fvc::interpolate(rho)
);
fvOptions.makeRelative(fvc::interpolate(psi), phid);
@@ -47,7 +49,7 @@
surfaceScalarField phiHbyA
(
"phiHbyA",
- fvc::interpolate(rho)*(fvc::interpolate(HbyA) & mesh.Sf())
+ fvc::interpolate(rho*HbyA) & mesh.Sf()
);
fvOptions.makeRelative(fvc::interpolate(rho), phiHbyA);
diff --git a/applications/solvers/compressible/rhoSimpleFoam/rhoSimplecFoam/pEqn.H b/applications/solvers/compressible/rhoSimpleFoam/rhoSimplecFoam/pEqn.H
index 353593b7f4..d61a80ffa8 100644
--- a/applications/solvers/compressible/rhoSimpleFoam/rhoSimplecFoam/pEqn.H
+++ b/applications/solvers/compressible/rhoSimpleFoam/rhoSimplecFoam/pEqn.H
@@ -13,7 +13,9 @@ if (simple.transonic())
surfaceScalarField phid
(
"phid",
- fvc::interpolate(psi)*(fvc::interpolate(HbyA) & mesh.Sf())
+ fvc::interpolate(psi)
+ *(fvc::interpolate(rho*HbyA) & mesh.Sf())
+ /fvc::interpolate(rho)
);
surfaceScalarField phic
@@ -57,7 +59,7 @@ else
surfaceScalarField phiHbyA
(
"phiHbyA",
- fvc::interpolate(rho)*(fvc::interpolate(HbyA) & mesh.Sf())
+ fvc::interpolate(rho*HbyA) & mesh.Sf()
);
closedVolume = adjustPhi(phiHbyA, U, p);
diff --git a/applications/solvers/compressible/sonicFoam/pEqn.H b/applications/solvers/compressible/sonicFoam/pEqn.H
index 4d700b3a41..2832bf421a 100644
--- a/applications/solvers/compressible/sonicFoam/pEqn.H
+++ b/applications/solvers/compressible/sonicFoam/pEqn.H
@@ -1,22 +1,21 @@
rho = thermo.rho();
volScalarField rAU(1.0/UEqn.A());
+surfaceScalarField Dp("Dp", fvc::interpolate(rho*rAU));
+
volVectorField HbyA("HbyA", U);
HbyA = rAU*UEqn.H();
surfaceScalarField phid
(
"phid",
- fvc::interpolate(psi)
- *(
- (fvc::interpolate(HbyA) & mesh.Sf())
- + fvc::ddtPhiCorr(rAU, rho, U, phi)
- )
+ fvc::interpolate(psi)*
+ (
+ (mesh.Sf() & fvc::interpolate(rho*HbyA))
+ + Dp*fvc::ddtCorr(rho, U, phi)
+ )/fvc::interpolate(rho)
);
-
-volScalarField Dp("Dp", rho*rAU);
-
// Non-orthogonal pressure corrector loop
while (pimple.correctNonOrthogonal())
{
diff --git a/applications/solvers/compressible/sonicFoam/sonicDyMFoam/createRhoUf.H b/applications/solvers/compressible/sonicFoam/sonicDyMFoam/createRhoUf.H
new file mode 100644
index 0000000000..bb05d10f66
--- /dev/null
+++ b/applications/solvers/compressible/sonicFoam/sonicDyMFoam/createRhoUf.H
@@ -0,0 +1,56 @@
+/*---------------------------------------------------------------------------*\
+ ========= |
+ \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
+ \\ / O peration |
+ \\ / A nd | Copyright (C) 2013 OpenFOAM Foundation
+ \\/ M anipulation |
+-------------------------------------------------------------------------------
+License
+ This file is part of OpenFOAM.
+
+ OpenFOAM is free software: you can redistribute it and/or modify it
+ under the terms of the GNU General Public License as published by
+ the Free Software Foundation, either version 3 of the License, or
+ (at your option) any later version.
+
+ OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+ ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+ FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
+ for more details.
+
+ You should have received a copy of the GNU General Public License
+ along with OpenFOAM. If not, see .
+
+Global
+ createUf
+
+Description
+ Creates and initialises the velocity velocity field Uf.
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef createUf_H
+#define createUf_H
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+Info<< "Reading/calculating face velocity rhoUf\n" << endl;
+
+surfaceVectorField rhoUf
+(
+ IOobject
+ (
+ "rhoUf",
+ runTime.timeName(),
+ mesh,
+ IOobject::READ_IF_PRESENT,
+ IOobject::AUTO_WRITE
+ ),
+ linearInterpolate(rho*U)
+);
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/applications/solvers/compressible/sonicFoam/sonicDyMFoam/pEqn.H b/applications/solvers/compressible/sonicFoam/sonicDyMFoam/pEqn.H
index 21dc48614e..91becd8aa1 100644
--- a/applications/solvers/compressible/sonicFoam/sonicDyMFoam/pEqn.H
+++ b/applications/solvers/compressible/sonicFoam/sonicDyMFoam/pEqn.H
@@ -1,6 +1,8 @@
rho = thermo.rho();
volScalarField rAU(1.0/UEqn.A());
+surfaceScalarField Dp("Dp", fvc::interpolate(rho*rAU));
+
volVectorField HbyA("HbyA", U);
HbyA = rAU*UEqn.H();
@@ -9,13 +11,14 @@ surfaceScalarField phid
"phid",
fvc::interpolate(psi)
*(
- (fvc::interpolate(HbyA) & mesh.Sf())
+ (
+ (fvc::interpolate(rho*HbyA) & mesh.Sf())
+ //***HGW + Dp*fvc::ddtCorr(rho, U, phi)
+ )/fvc::interpolate(rho)
- fvc::meshPhi(rho, U)
)
);
-volScalarField Dp("Dp", rho*rAU);
-
for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
{
fvScalarMatrix pEqn
diff --git a/applications/solvers/compressible/sonicFoam/sonicDyMFoam/sonicDyMFoam.C b/applications/solvers/compressible/sonicFoam/sonicDyMFoam/sonicDyMFoam.C
index 2f7ac9962a..2504f19766 100644
--- a/applications/solvers/compressible/sonicFoam/sonicDyMFoam/sonicDyMFoam.C
+++ b/applications/solvers/compressible/sonicFoam/sonicDyMFoam/sonicDyMFoam.C
@@ -44,6 +44,7 @@ int main(int argc, char *argv[])
#include "createTime.H"
#include "createMesh.H"
#include "createFields.H"
+ #include "createRhoUf.H"
#include "initContinuityErrs.H"
pimpleControl pimple(mesh);
@@ -63,6 +64,12 @@ 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/compressible/sonicFoam/sonicLiquidFoam/sonicLiquidFoam.C b/applications/solvers/compressible/sonicFoam/sonicLiquidFoam/sonicLiquidFoam.C
index 1ac1ae4ae3..cc7f734afc 100644
--- a/applications/solvers/compressible/sonicFoam/sonicLiquidFoam/sonicLiquidFoam.C
+++ b/applications/solvers/compressible/sonicFoam/sonicLiquidFoam/sonicLiquidFoam.C
@@ -76,6 +76,8 @@ int main(int argc, char *argv[])
while (pimple.correct())
{
volScalarField rAU(1.0/UEqn.A());
+ surfaceScalarField Dp("Dp", fvc::interpolate(rho*rAU));
+
U = rAU*UEqn.H();
surfaceScalarField phid
@@ -83,13 +85,12 @@ int main(int argc, char *argv[])
"phid",
psi
*(
- (fvc::interpolate(U) & mesh.Sf())
- + fvc::ddtPhiCorr(rAU, rho, U, phi)
- )
+ (fvc::interpolate(rho*U) & mesh.Sf())
+ + Dp*fvc::ddtCorr(rho, U, phi)
+ )/fvc::interpolate(rho)
);
phi = (rhoO/psi)*phid;
- volScalarField Dp("Dp", rho*rAU);
fvScalarMatrix pEqn
(
diff --git a/applications/solvers/electromagnetics/mhdFoam/mhdFoam.C b/applications/solvers/electromagnetics/mhdFoam/mhdFoam.C
index 51f04820eb..d920e41000 100644
--- a/applications/solvers/electromagnetics/mhdFoam/mhdFoam.C
+++ b/applications/solvers/electromagnetics/mhdFoam/mhdFoam.C
@@ -93,6 +93,8 @@ int main(int argc, char *argv[])
for (int corr=0; corr turbulence
- (
- incompressible::turbulenceModel::New(U, phi, laminarTransport)
- );
diff --git a/applications/solvers/incompressible/pimpleFoam/pimpleDyMFoam/createUf.H b/applications/solvers/incompressible/pimpleFoam/pimpleDyMFoam/createUf.H
new file mode 100644
index 0000000000..aab92ce0e2
--- /dev/null
+++ b/applications/solvers/incompressible/pimpleFoam/pimpleDyMFoam/createUf.H
@@ -0,0 +1,56 @@
+/*---------------------------------------------------------------------------*\
+ ========= |
+ \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
+ \\ / O peration |
+ \\ / A nd | Copyright (C) 2013 OpenFOAM Foundation
+ \\/ M anipulation |
+-------------------------------------------------------------------------------
+License
+ This file is part of OpenFOAM.
+
+ OpenFOAM is free software: you can redistribute it and/or modify it
+ under the terms of the GNU General Public License as published by
+ the Free Software Foundation, either version 3 of the License, or
+ (at your option) any later version.
+
+ OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+ ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+ FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
+ for more details.
+
+ You should have received a copy of the GNU General Public License
+ along with OpenFOAM. If not, see .
+
+Global
+ createUf
+
+Description
+ Creates and initialises the velocity velocity field Uf.
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef createUf_H
+#define createUf_H
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+Info<< "Reading/calculating face velocity Uf\n" << endl;
+
+surfaceVectorField Uf
+(
+ IOobject
+ (
+ "Uf",
+ runTime.timeName(),
+ mesh,
+ IOobject::READ_IF_PRESENT,
+ IOobject::AUTO_WRITE
+ ),
+ linearInterpolate(U)
+);
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/applications/solvers/incompressible/pimpleFoam/pimpleDyMFoam/pEqn.H b/applications/solvers/incompressible/pimpleFoam/pimpleDyMFoam/pEqn.H
index c9c0ce3c58..827daa4ce8 100644
--- a/applications/solvers/incompressible/pimpleFoam/pimpleDyMFoam/pEqn.H
+++ b/applications/solvers/incompressible/pimpleFoam/pimpleDyMFoam/pEqn.H
@@ -1,3 +1,5 @@
+surfaceScalarField Dp("Dp", fvc::interpolate(rAU));
+
volVectorField HbyA("HbyA", U);
HbyA = rAU*UEqn().H();
@@ -10,7 +12,7 @@ surfaceScalarField phiHbyA
(
"phiHbyA",
(fvc::interpolate(HbyA) & mesh.Sf())
- + fvc::ddtPhiCorr(rAU, U, phiAbs)
+ + Dp*fvc::ddtCorr(U, Uf)
);
if (p.needReference())
@@ -24,7 +26,7 @@ while (pimple.correctNonOrthogonal())
{
fvScalarMatrix pEqn
(
- fvm::laplacian(rAU, p) == fvc::div(phiHbyA)
+ fvm::laplacian(Dp, p) == fvc::div(phiHbyA)
);
pEqn.setReference(pRefCell, pRefValue);
@@ -42,9 +44,15 @@ while (pimple.correctNonOrthogonal())
// Explicitly relax pressure for momentum corrector
p.relax();
-// Make the fluxes relative to the mesh motion
-fvc::makeRelative(phi, U);
-
U = HbyA - rAU*fvc::grad(p);
U.correctBoundaryConditions();
fvOptions.correct(U);
+
+{
+ Uf = fvc::interpolate(U);
+ surfaceVectorField n(mesh.Sf()/mesh.magSf());
+ Uf += mesh.Sf()*(phi - (mesh.Sf() & Uf))/sqr(mesh.magSf());
+}
+
+// Make the fluxes relative to the mesh motion
+fvc::makeRelative(phi, U);
diff --git a/applications/solvers/incompressible/pimpleFoam/pimpleDyMFoam/pimpleDyMFoam.C b/applications/solvers/incompressible/pimpleFoam/pimpleDyMFoam/pimpleDyMFoam.C
index 883a20f6cf..12b0b438a3 100644
--- a/applications/solvers/incompressible/pimpleFoam/pimpleDyMFoam/pimpleDyMFoam.C
+++ b/applications/solvers/incompressible/pimpleFoam/pimpleDyMFoam/pimpleDyMFoam.C
@@ -51,15 +51,13 @@ int main(int argc, char *argv[])
pimpleControl pimple(mesh);
#include "createFields.H"
+ #include "createUf.H"
#include "createFvOptions.H"
#include "readTimeControls.H"
#include "createPcorrTypes.H"
#include "CourantNo.H"
#include "setInitialDeltaT.H"
- // Create old-time absolute flux for ddtPhiCorr
- surfaceScalarField phiAbs("phiAbs", phi);
-
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
Info<< "\nStarting time loop\n" << endl;
@@ -69,12 +67,6 @@ int main(int argc, char *argv[])
#include "readControls.H"
#include "CourantNo.H"
- // Make the fluxes absolute
- fvc::makeAbsolute(phi, U);
-
- // Update absolute flux for ddtPhiCorr
- phiAbs = phi;
-
#include "setDeltaT.H"
runTime++;
@@ -83,12 +75,15 @@ int main(int argc, char *argv[])
mesh.update();
+ // Calculate absolute flux from the mapped surface velocity
+ phi = mesh.Sf() & Uf;
+
if (mesh.changing() && correctPhi)
{
#include "correctPhi.H"
}
- // Make the fluxes relative to the mesh motion
+ // Make the flux relative to the mesh motion
fvc::makeRelative(phi, U);
if (mesh.changing() && checkMeshCourantNo)
diff --git a/applications/solvers/incompressible/pisoFoam/pisoFoam.C b/applications/solvers/incompressible/pisoFoam/pisoFoam.C
index 1e45897e75..ecbb02ff20 100644
--- a/applications/solvers/incompressible/pisoFoam/pisoFoam.C
+++ b/applications/solvers/incompressible/pisoFoam/pisoFoam.C
@@ -87,7 +87,7 @@ int main(int argc, char *argv[])
(
"phiHbyA",
(fvc::interpolate(HbyA) & mesh.Sf())
- + fvc::ddtPhiCorr(rAU, U, phi)
+ + fvc::interpolate(rAU)*fvc::ddtCorr(U, phi)
);
adjustPhi(phiHbyA, U, p);
diff --git a/applications/solvers/incompressible/shallowWaterFoam/shallowWaterFoam.C b/applications/solvers/incompressible/shallowWaterFoam/shallowWaterFoam.C
index f279f91eec..2107e789af 100644
--- a/applications/solvers/incompressible/shallowWaterFoam/shallowWaterFoam.C
+++ b/applications/solvers/incompressible/shallowWaterFoam/shallowWaterFoam.C
@@ -111,7 +111,7 @@ int main(int argc, char *argv[])
(
"phiHbyA",
(fvc::interpolate(HbyA) & mesh.Sf())
- + fvc::ddtPhiCorr(rAU, h, hU, phi)
+ + fvc::interpolate(rAU)*fvc::ddtCorr(h, hU, phi)
- phih0
);
diff --git a/applications/solvers/lagrangian/coalChemistryFoam/pEqn.H b/applications/solvers/lagrangian/coalChemistryFoam/pEqn.H
index 82f1e9ff63..200ea289d8 100644
--- a/applications/solvers/lagrangian/coalChemistryFoam/pEqn.H
+++ b/applications/solvers/lagrangian/coalChemistryFoam/pEqn.H
@@ -1,6 +1,8 @@
rho = thermo.rho();
volScalarField rAU(1.0/UEqn.A());
+surfaceScalarField Dp("Dp", fvc::interpolate(rho*rAU));
+
volVectorField HbyA("HbyA", U);
HbyA = rAU*UEqn.H();
@@ -11,9 +13,9 @@ if (pimple.transonic())
"phid",
fvc::interpolate(psi)
*(
- (fvc::interpolate(HbyA) & mesh.Sf())
- + fvc::ddtPhiCorr(rAU, rho, U, phi)
- )
+ (fvc::interpolate(rho*HbyA) & mesh.Sf())
+ + Dp*fvc::ddtCorr(rho, U, phi)
+ )/fvc::interpolate(rho)
);
fvOptions.makeRelative(fvc::interpolate(psi), phid);
@@ -24,7 +26,7 @@ if (pimple.transonic())
(
fvm::ddt(psi, p)
+ fvm::div(phid, p)
- - fvm::laplacian(rho*rAU, p)
+ - fvm::laplacian(Dp, p)
==
coalParcels.Srho()
+ fvOptions(psi, p, rho.name())
@@ -45,10 +47,9 @@ else
surfaceScalarField phiHbyA
(
"phiHbyA",
- fvc::interpolate(rho)
- *(
- (fvc::interpolate(HbyA) & mesh.Sf())
- + fvc::ddtPhiCorr(rAU, rho, U, phi)
+ (
+ (fvc::interpolate(rho*HbyA) & mesh.Sf())
+ + Dp*fvc::ddtCorr(rho, U, phi)
)
);
@@ -60,7 +61,7 @@ else
(
fvm::ddt(psi, p)
+ fvc::div(phiHbyA)
- - fvm::laplacian(rho*rAU, p)
+ - fvm::laplacian(Dp, p)
==
coalParcels.Srho()
+ fvOptions(psi, p, rho.name())
diff --git a/applications/solvers/lagrangian/reactingParcelFilmFoam/pEqn.H b/applications/solvers/lagrangian/reactingParcelFilmFoam/pEqn.H
index 7df99b9b94..1e504462b6 100644
--- a/applications/solvers/lagrangian/reactingParcelFilmFoam/pEqn.H
+++ b/applications/solvers/lagrangian/reactingParcelFilmFoam/pEqn.H
@@ -1,19 +1,19 @@
rho = thermo.rho();
volScalarField rAU(1.0/UEqn.A());
-surfaceScalarField rhorAUf("Dp", fvc::interpolate(rho*rAU));
+surfaceScalarField Dp("Dp", fvc::interpolate(rho*rAU));
+
volVectorField HbyA("HbyA", U);
HbyA = rAU*UEqn.H();
-surfaceScalarField phig(-rhorAUf*ghf*fvc::snGrad(rho)*mesh.magSf());
+surfaceScalarField phig(-Dp*ghf*fvc::snGrad(rho)*mesh.magSf());
surfaceScalarField phiHbyA
(
"phiHbyA",
- fvc::interpolate(rho)
- *(
- (fvc::interpolate(HbyA) & mesh.Sf())
- + fvc::ddtPhiCorr(rAU, rho, U, phi)
+ (
+ (fvc::interpolate(rho*HbyA) & mesh.Sf())
+ + Dp*fvc::ddtCorr(rho, U, phi)
)
+ phig
);
@@ -27,7 +27,7 @@ while (pimple.correctNonOrthogonal())
fvc::ddt(psi, rho)*gh
+ fvc::div(phiHbyA)
+ fvm::ddt(psi, p_rgh)
- - fvm::laplacian(rhorAUf, p_rgh)
+ - fvm::laplacian(Dp, p_rgh)
==
parcels.Srho()
+ surfaceFilm.Srho()
@@ -41,7 +41,7 @@ while (pimple.correctNonOrthogonal())
if (pimple.finalNonOrthogonalIter())
{
phi = phiHbyA + p_rghEqn.flux();
- U = HbyA + rAU*fvc::reconstruct((p_rghEqn.flux() + phig)/rhorAUf);
+ U = HbyA + rAU*fvc::reconstruct((p_rghEqn.flux() + phig)/Dp);
U.correctBoundaryConditions();
fvOptions.correct(U);
}
diff --git a/applications/solvers/lagrangian/reactingParcelFoam/LTSReactingParcelFoam/LTSReactingParcelFoam.C b/applications/solvers/lagrangian/reactingParcelFoam/LTSReactingParcelFoam/LTSReactingParcelFoam.C
index 34aecb65b4..7481ae79cd 100644
--- a/applications/solvers/lagrangian/reactingParcelFoam/LTSReactingParcelFoam/LTSReactingParcelFoam.C
+++ b/applications/solvers/lagrangian/reactingParcelFoam/LTSReactingParcelFoam/LTSReactingParcelFoam.C
@@ -30,9 +30,6 @@ Description
parcels and porous media, including run-time selectable finitite volume
options, e.g. sources, constraints
- Note: ddtPhiCorr not used here when porous zones are active
- - not well defined for porous calculations
-
\*---------------------------------------------------------------------------*/
#include "fvCFD.H"
diff --git a/applications/solvers/lagrangian/reactingParcelFoam/pEqn.H b/applications/solvers/lagrangian/reactingParcelFoam/pEqn.H
index 6dd16f0bdb..420bc1ef0c 100644
--- a/applications/solvers/lagrangian/reactingParcelFoam/pEqn.H
+++ b/applications/solvers/lagrangian/reactingParcelFoam/pEqn.H
@@ -6,16 +6,17 @@
thermo.rho() -= psi*p;
volScalarField rAU(1.0/UEqn.A());
+ surfaceScalarField Dp("Dp", fvc::interpolate(rho*rAU));
+
volVectorField HbyA("HbyA", U);
HbyA = rAU*UEqn.H();
surfaceScalarField phiHbyA
(
"phiHbyA",
- fvc::interpolate(rho)
- *(
- (fvc::interpolate(HbyA) & mesh.Sf())
- + fvc::ddtPhiCorr(rAU, rho, U, phi)
+ (
+ (fvc::interpolate(rho*HbyA) & mesh.Sf())
+ + Dp*fvc::ddtCorr(rho, U, phi)
)
);
@@ -35,7 +36,7 @@
fvScalarMatrix pEqn
(
pDDtEqn
- - fvm::laplacian(rho*rAU, p)
+ - fvm::laplacian(Dp, p)
);
fvOptions.constrain(pEqn);
diff --git a/applications/solvers/lagrangian/sprayFoam/pEqn.H b/applications/solvers/lagrangian/sprayFoam/pEqn.H
index af65d3016b..f62a17286c 100644
--- a/applications/solvers/lagrangian/sprayFoam/pEqn.H
+++ b/applications/solvers/lagrangian/sprayFoam/pEqn.H
@@ -1,6 +1,8 @@
rho = thermo.rho();
volScalarField rAU(1.0/UEqn.A());
+surfaceScalarField Dp("Dp", fvc::interpolate(rho*rAU));
+
volVectorField HbyA("HbyA", U);
HbyA = rAU*UEqn.H();
@@ -11,9 +13,9 @@ if (pimple.transonic())
"phid",
fvc::interpolate(psi)
*(
- (fvc::interpolate(HbyA) & mesh.Sf())
- + fvc::ddtPhiCorr(rAU, rho, U, phi)
- )
+ (fvc::interpolate(rho*HbyA) & mesh.Sf())
+ + Dp*fvc::ddtCorr(rho, U, phi)
+ )/fvc::interpolate(rho)
);
fvOptions.makeRelative(fvc::interpolate(psi), phid);
@@ -24,7 +26,7 @@ if (pimple.transonic())
(
fvm::ddt(psi, p)
+ fvm::div(phid, p)
- - fvm::laplacian(rho*rAU, p)
+ - fvm::laplacian(Dp, p)
==
parcels.Srho()
+ fvOptions(psi, p, rho.name())
@@ -45,10 +47,9 @@ else
surfaceScalarField phiHbyA
(
"phiHbyA",
- fvc::interpolate(rho)
- *(
- (fvc::interpolate(HbyA) & mesh.Sf())
- + fvc::ddtPhiCorr(rAU, rho, U, phi)
+ (
+ (fvc::interpolate(rho*HbyA) & mesh.Sf())
+ + Dp*fvc::ddtCorr(rho, U, phi)
)
);
@@ -60,7 +61,7 @@ else
(
fvm::ddt(psi, p)
+ fvc::div(phiHbyA)
- - fvm::laplacian(rho*rAU, p)
+ - fvm::laplacian(Dp, p)
==
parcels.Srho()
+ fvOptions(psi, p, rho.name())
diff --git a/applications/solvers/lagrangian/sprayFoam/sprayEngineFoam/pEqn.H b/applications/solvers/lagrangian/sprayFoam/sprayEngineFoam/pEqn.H
index cf7acbd602..b7de7e2f79 100644
--- a/applications/solvers/lagrangian/sprayFoam/sprayEngineFoam/pEqn.H
+++ b/applications/solvers/lagrangian/sprayFoam/sprayEngineFoam/pEqn.H
@@ -1,6 +1,8 @@
rho = thermo.rho();
volScalarField rAU(1.0/UEqn.A());
+surfaceScalarField Dp("Dp", fvc::interpolate(rho*rAU));
+
volVectorField HbyA("HbyA", U);
HbyA = rAU*UEqn.H();
@@ -11,8 +13,11 @@ if (pimple.transonic())
"phid",
fvc::interpolate(psi)
*(
- ((fvc::interpolate(HbyA) & mesh.Sf()) - fvc::meshPhi(rho, U))
- + fvc::ddtPhiCorr(rAU, rho, U, phi)
+ (
+ (fvc::interpolate(rho*HbyA) & mesh.Sf())
+ //***HGW + Dp*fvc::ddtCorr(rho, U, phi)
+ )/fvc::interpolate(rho)
+ - fvc::meshPhi(rho, U)
)
);
@@ -24,7 +29,7 @@ if (pimple.transonic())
(
fvm::ddt(psi, p)
+ fvm::div(phid, p)
- - fvm::laplacian(rho*rAU, p)
+ - fvm::laplacian(Dp, p)
==
parcels.Srho()
+ fvOptions(psi, p, rho.name())
@@ -45,11 +50,11 @@ else
surfaceScalarField phiHbyA
(
"phiHbyA",
- fvc::interpolate(rho)
- *(
- ((fvc::interpolate(HbyA) & mesh.Sf()) - fvc::meshPhi(rho, U))
- + fvc::ddtPhiCorr(rAU, rho, U, phi)
+ (
+ (fvc::interpolate(HbyA) & mesh.Sf())
+ //***HGW + Dp*fvc::ddtCorr(rho, U, phi)
)
+ - fvc::interpolate(rho)*fvc::meshPhi(rho, U)
);
fvOptions.makeRelative(fvc::interpolate(rho), phiHbyA);
@@ -60,7 +65,7 @@ else
(
fvm::ddt(psi, p)
+ fvc::div(phiHbyA)
- - fvm::laplacian(rho*rAU, p)
+ - fvm::laplacian(Dp, p)
==
parcels.Srho()
+ fvOptions(psi, p, rho.name())
diff --git a/applications/solvers/multiphase/cavitatingFoam/cavitatingDyMFoam/pEqn.H b/applications/solvers/multiphase/cavitatingFoam/cavitatingDyMFoam/pEqn.H
index 5ad7fdeb66..f869961a3c 100644
--- a/applications/solvers/multiphase/cavitatingFoam/cavitatingDyMFoam/pEqn.H
+++ b/applications/solvers/multiphase/cavitatingFoam/cavitatingDyMFoam/pEqn.H
@@ -12,16 +12,16 @@
surfaceScalarField rhof("rhof", fvc::interpolate(rho));
volScalarField rAU(1.0/UEqn.A());
- surfaceScalarField rAUf("Dp", rhof*fvc::interpolate(rAU));
+ surfaceScalarField Dp("Dp", fvc::interpolate(rho*rAU));
volVectorField HbyA("HbyA", U);
HbyA = rAU*UEqn.H();
phiv = (fvc::interpolate(HbyA) & mesh.Sf())
- + fvc::ddtPhiCorr(rAU, rho, U, phivAbs);
+ + Dp*fvc::ddtCorr(U, phivAbs);
fvc::makeRelative(phiv, U);
- surfaceScalarField phiGradp(rAUf*mesh.magSf()*fvc::snGrad(p));
+ surfaceScalarField phiGradp(Dp*mesh.magSf()*fvc::snGrad(p));
phiv -= phiGradp/rhof;
@@ -35,7 +35,7 @@
+ psi*correction(fvm::ddt(p))
+ fvc::div(phiv, rho)
+ fvc::div(phiGradp)
- - fvm::laplacian(rAUf, p)
+ - fvm::laplacian(Dp, p)
);
pEqn.solve(mesh.solver(p.select(pimple.finalInnerIter())));
diff --git a/applications/solvers/multiphase/cavitatingFoam/pEqn.H b/applications/solvers/multiphase/cavitatingFoam/pEqn.H
index da43e67ce7..e9d059aadb 100644
--- a/applications/solvers/multiphase/cavitatingFoam/pEqn.H
+++ b/applications/solvers/multiphase/cavitatingFoam/pEqn.H
@@ -12,15 +12,15 @@
surfaceScalarField rhof("rhof", fvc::interpolate(rho));
volScalarField rAU(1.0/UEqn.A());
- surfaceScalarField rAUf("Dp", rhof*fvc::interpolate(rAU));
+ surfaceScalarField Dp("Dp", fvc::interpolate(rho*rAU));
volVectorField HbyA("HbyA", U);
HbyA = rAU*UEqn.H();
phiv = (fvc::interpolate(HbyA) & mesh.Sf())
- + fvc::ddtPhiCorr(rAU, rho, U, phiv);
+ + Dp*fvc::ddtCorr(U, phiv);
- surfaceScalarField phiGradp(rAUf*mesh.magSf()*fvc::snGrad(p));
+ surfaceScalarField phiGradp(Dp*mesh.magSf()*fvc::snGrad(p));
phiv -= phiGradp/rhof;
@@ -32,7 +32,7 @@
- (rhol0 + (psil - psiv)*pSat)*fvc::ddt(alphav) - pSat*fvc::ddt(psi)
+ fvc::div(phiv, rho)
+ fvc::div(phiGradp)
- - fvm::laplacian(rAUf, p)
+ - fvm::laplacian(Dp, p)
);
pEqn.solve(mesh.solver(p.select(pimple.finalInnerIter())));
diff --git a/applications/solvers/multiphase/compressibleInterFoam/pEqn.H b/applications/solvers/multiphase/compressibleInterFoam/pEqn.H
index 1319378821..adc380978c 100644
--- a/applications/solvers/multiphase/compressibleInterFoam/pEqn.H
+++ b/applications/solvers/multiphase/compressibleInterFoam/pEqn.H
@@ -1,6 +1,6 @@
{
volScalarField rAU("rAU", 1.0/UEqn.A());
- surfaceScalarField rAUf("Dp", fvc::interpolate(rAU));
+ surfaceScalarField Dp("Dp", fvc::interpolate(rAU));
volVectorField HbyA("HbyA", U);
HbyA = rAU*UEqn.H();
@@ -9,7 +9,7 @@
(
"phiHbyA",
(fvc::interpolate(HbyA) & mesh.Sf())
- + fvc::ddtPhiCorr(rAU, rho, U, phi)
+ + fvc::interpolate(rho*rAU)*fvc::ddtCorr(U, phi)
);
phi = phiHbyA;
@@ -18,7 +18,7 @@
(
fvc::interpolate(interface.sigmaK())*fvc::snGrad(alpha1)
- ghf*fvc::snGrad(rho)
- )*rAUf*mesh.magSf()
+ )*Dp*mesh.magSf()
);
phiHbyA += phig;
@@ -70,7 +70,7 @@
fvScalarMatrix p_rghEqnIncomp
(
fvc::div(phiHbyA)
- - fvm::laplacian(rAUf, p_rgh)
+ - fvm::laplacian(Dp, p_rgh)
);
solve
@@ -97,7 +97,7 @@
phi = phiHbyA + p_rghEqnIncomp.flux();
U = HbyA
- + rAU*fvc::reconstruct((phig + p_rghEqnIncomp.flux())/rAUf);
+ + rAU*fvc::reconstruct((phig + p_rghEqnIncomp.flux())/Dp);
U.correctBoundaryConditions();
}
}
diff --git a/applications/solvers/multiphase/interFoam/MRFInterFoam/pEqn.H b/applications/solvers/multiphase/interFoam/MRFInterFoam/pEqn.H
index 5d2acc8b68..be0e723493 100644
--- a/applications/solvers/multiphase/interFoam/MRFInterFoam/pEqn.H
+++ b/applications/solvers/multiphase/interFoam/MRFInterFoam/pEqn.H
@@ -1,6 +1,6 @@
{
volScalarField rAU("rAU", 1.0/UEqn.A());
- surfaceScalarField rAUf("Dp", fvc::interpolate(rAU));
+ surfaceScalarField Dp("Dp", fvc::interpolate(rAU));
volVectorField HbyA("HbyA", U);
HbyA = rAU*UEqn.H();
@@ -9,7 +9,7 @@
(
"phiHbyA",
(fvc::interpolate(HbyA) & mesh.Sf())
- + fvc::ddtPhiCorr(rAU, rho, U, phi)
+ + fvc::interpolate(rho*rAU)*fvc::ddtCorr(U, phi)
);
adjustPhi(phiHbyA, U, p_rgh);
mrfZones.makeRelative(phiHbyA);
@@ -20,7 +20,7 @@
(
fvc::interpolate(interface.sigmaK())*fvc::snGrad(alpha1)
- ghf*fvc::snGrad(rho)
- )*rAUf*mesh.magSf()
+ )*Dp*mesh.magSf()
);
phiHbyA += phig;
@@ -29,7 +29,7 @@
{
fvScalarMatrix p_rghEqn
(
- fvm::laplacian(rAUf, p_rgh) == fvc::div(phiHbyA)
+ fvm::laplacian(Dp, p_rgh) == fvc::div(phiHbyA)
);
p_rghEqn.setReference(pRefCell, getRefCellValue(p_rgh, pRefCell));
@@ -40,7 +40,7 @@
{
phi = phiHbyA - p_rghEqn.flux();
- U = HbyA + rAU*fvc::reconstruct((phig - p_rghEqn.flux())/rAUf);
+ U = HbyA + rAU*fvc::reconstruct((phig - p_rghEqn.flux())/Dp);
U.correctBoundaryConditions();
fvOptions.correct(U);
}
diff --git a/applications/solvers/multiphase/interFoam/interDyMFoam/correctPhi.H b/applications/solvers/multiphase/interFoam/interDyMFoam/correctPhi.H
index 75560b6337..deb4edc388 100644
--- a/applications/solvers/multiphase/interFoam/interDyMFoam/correctPhi.H
+++ b/applications/solvers/multiphase/interFoam/interDyMFoam/correctPhi.H
@@ -1,20 +1,27 @@
+if (mesh.changing())
{
- #include "continuityErrs.H"
-
- wordList pcorrTypes
- (
- p_rgh.boundaryField().size(),
- zeroGradientFvPatchScalarField::typeName
- );
-
- forAll (p_rgh.boundaryField(), i)
+ forAll(U.boundaryField(), patchI)
{
- if (p_rgh.boundaryField()[i].fixesValue())
+ if (U.boundaryField()[patchI].fixesValue())
{
- pcorrTypes[i] = fixedValueFvPatchScalarField::typeName;
+ U.boundaryField()[patchI].initEvaluate();
}
}
+ forAll(U.boundaryField(), patchI)
+ {
+ if (U.boundaryField()[patchI].fixesValue())
+ {
+ U.boundaryField()[patchI].evaluate();
+
+ phi.boundaryField()[patchI] =
+ U.boundaryField()[patchI]
+ & mesh.Sf().boundaryField()[patchI];
+ }
+ }
+}
+
+{
volScalarField pcorr
(
IOobject
@@ -30,17 +37,13 @@
pcorrTypes
);
- dimensionedScalar rAUf("(1|A(U))", dimTime/rho.dimensions(), 1.0);
-
- adjustPhi(phi, U, pcorr);
-
- fvc::makeAbsolute(phi, U);
+ dimensionedScalar Dp("Dp", dimTime/rho.dimensions(), 1.0);
while (pimple.correctNonOrthogonal())
{
fvScalarMatrix pcorrEqn
(
- fvm::laplacian(rAUf, pcorr) == fvc::div(phi)
+ fvm::laplacian(Dp, pcorr) == fvc::div(phi)
);
pcorrEqn.setReference(pRefCell, pRefValue);
@@ -49,9 +52,6 @@
if (pimple.finalNonOrthogonalIter())
{
phi -= pcorrEqn.flux();
- phiAbs = phi;
- phiAbs.oldTime() = phi;
- fvc::makeRelative(phi, U);
}
}
}
diff --git a/applications/solvers/multiphase/interFoam/interDyMFoam/createPcorrTypes.H b/applications/solvers/multiphase/interFoam/interDyMFoam/createPcorrTypes.H
new file mode 100644
index 0000000000..dfd4afb49b
--- /dev/null
+++ b/applications/solvers/multiphase/interFoam/interDyMFoam/createPcorrTypes.H
@@ -0,0 +1,13 @@
+ wordList pcorrTypes
+ (
+ p_rgh.boundaryField().size(),
+ zeroGradientFvPatchScalarField::typeName
+ );
+
+ for (label i=0; i.
+
+Global
+ createUf
+
+Description
+ Creates and initialises the velocity velocity field Uf.
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef createUf_H
+#define createUf_H
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+Info<< "Reading/calculating face velocity Uf\n" << endl;
+
+surfaceVectorField Uf
+(
+ IOobject
+ (
+ "Uf",
+ runTime.timeName(),
+ mesh,
+ IOobject::READ_IF_PRESENT,
+ IOobject::AUTO_WRITE
+ ),
+ linearInterpolate(U)
+);
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/applications/solvers/multiphase/interFoam/interDyMFoam/interDyMFoam.C b/applications/solvers/multiphase/interFoam/interDyMFoam/interDyMFoam.C
index 0c22dd01b7..75e4c2f8d4 100644
--- a/applications/solvers/multiphase/interFoam/interDyMFoam/interDyMFoam.C
+++ b/applications/solvers/multiphase/interFoam/interDyMFoam/interDyMFoam.C
@@ -50,15 +50,13 @@ int main(int argc, char *argv[])
#include "createTime.H"
#include "createDynamicFvMesh.H"
#include "initContinuityErrs.H"
- #include "createFields.H"
pimpleControl pimple(mesh);
+ #include "createFields.H"
+ #include "createUf.H"
#include "readTimeControls.H"
-
- surfaceScalarField phiAbs("phiAbs", phi);
- fvc::makeAbsolute(phiAbs, U);
-
+ #include "createPcorrTypes.H"
#include "correctPhi.H"
#include "CourantNo.H"
#include "setInitialDeltaT.H"
@@ -80,21 +78,7 @@ int main(int argc, char *argv[])
scalar timeBeforeMeshUpdate = runTime.elapsedCpuTime();
- {
- // Ensure old-time U exists for mapping
- U.oldTime();
-
- // Calculate the relative velocity used to map the relative flux phi
- volVectorField Urel("Urel", U);
-
- if (mesh.moving())
- {
- Urel -= fvc::reconstruct(fvc::meshPhi(U));
- }
-
- // Do any mesh changes
- mesh.update();
- }
+ mesh.update();
if (mesh.changing())
{
@@ -108,7 +92,13 @@ int main(int argc, char *argv[])
if (mesh.changing() && correctPhi)
{
+ // Calculate absolute flux from the mapped surface velocity
+ phi = mesh.Sf() & Uf;
+
#include "correctPhi.H"
+
+ // Make the flux relative to the mesh motion
+ fvc::makeRelative(phi, U);
}
if (mesh.changing() && checkMeshCourantNo)
diff --git a/applications/solvers/multiphase/interFoam/interDyMFoam/pEqn.H b/applications/solvers/multiphase/interFoam/interDyMFoam/pEqn.H
index 6fa4330997..51190711c6 100644
--- a/applications/solvers/multiphase/interFoam/interDyMFoam/pEqn.H
+++ b/applications/solvers/multiphase/interFoam/interDyMFoam/pEqn.H
@@ -1,6 +1,6 @@
{
volScalarField rAU("rAU", 1.0/UEqn.A());
- surfaceScalarField rAUf("Dp", fvc::interpolate(rAU));
+ surfaceScalarField Dp("Dp", fvc::interpolate(rAU));
volVectorField HbyA("HbyA", U);
HbyA = rAU*UEqn.H();
@@ -9,7 +9,7 @@
(
"phiHbyA",
(fvc::interpolate(HbyA) & mesh.Sf())
- + fvc::ddtPhiCorr(rAU, rho, U, phiAbs)
+ + fvc::interpolate(rho*rAU)*fvc::ddtCorr(U, Uf)
);
if (p_rgh.needReference())
@@ -19,14 +19,14 @@
fvc::makeAbsolute(phiHbyA, U);
}
- phiAbs = phiHbyA;
+ surfaceScalarField phiAbs("phiAbs", phiHbyA);
surfaceScalarField phig
(
(
fvc::interpolate(interface.sigmaK())*fvc::snGrad(alpha1)
- ghf*fvc::snGrad(rho)
- )*rAUf*mesh.magSf()
+ )*Dp*mesh.magSf()
);
phiHbyA += phig;
@@ -35,7 +35,7 @@
{
fvScalarMatrix p_rghEqn
(
- fvm::laplacian(rAUf, p_rgh) == fvc::div(phiHbyA)
+ fvm::laplacian(Dp, p_rgh) == fvc::div(phiHbyA)
);
p_rghEqn.setReference(pRefCell, getRefCellValue(p_rgh, pRefCell));
@@ -46,7 +46,7 @@
{
phi = phiHbyA - p_rghEqn.flux();
- U = HbyA + rAU*fvc::reconstruct((phig - p_rghEqn.flux())/rAUf);
+ U = HbyA + rAU*fvc::reconstruct((phig - p_rghEqn.flux())/Dp);
U.correctBoundaryConditions();
fvOptions.correct(U);
}
@@ -54,7 +54,11 @@
#include "continuityErrs.H"
- phiAbs = phi;
+ {
+ Uf = fvc::interpolate(U);
+ surfaceVectorField n(mesh.Sf()/mesh.magSf());
+ Uf += mesh.Sf()*(phi - (mesh.Sf() & Uf))/sqr(mesh.magSf());
+ }
// Make the fluxes relative to the mesh motion
fvc::makeRelative(phi, U);
diff --git a/applications/solvers/multiphase/interFoam/pEqn.H b/applications/solvers/multiphase/interFoam/pEqn.H
index 0ec531f0dd..d8651bfa81 100644
--- a/applications/solvers/multiphase/interFoam/pEqn.H
+++ b/applications/solvers/multiphase/interFoam/pEqn.H
@@ -1,6 +1,6 @@
{
volScalarField rAU("rAU", 1.0/UEqn.A());
- surfaceScalarField rAUf("Dp", fvc::interpolate(rAU));
+ surfaceScalarField Dp("Dp", fvc::interpolate(rAU));
volVectorField HbyA("HbyA", U);
HbyA = rAU*UEqn.H();
@@ -9,7 +9,7 @@
(
"phiHbyA",
(fvc::interpolate(HbyA) & mesh.Sf())
- + fvc::ddtPhiCorr(rAU, rho, U, phi)
+ + fvc::interpolate(rho*rAU)*fvc::ddtCorr(U, phi)
);
adjustPhi(phiHbyA, U, p_rgh);
@@ -20,7 +20,7 @@
(
fvc::interpolate(interface.sigmaK())*fvc::snGrad(alpha1)
- ghf*fvc::snGrad(rho)
- )*rAUf*mesh.magSf()
+ )*Dp*mesh.magSf()
);
phiHbyA += phig;
@@ -29,7 +29,7 @@
{
fvScalarMatrix p_rghEqn
(
- fvm::laplacian(rAUf, p_rgh) == fvc::div(phiHbyA)
+ fvm::laplacian(Dp, p_rgh) == fvc::div(phiHbyA)
);
p_rghEqn.setReference(pRefCell, getRefCellValue(p_rgh, pRefCell));
@@ -40,7 +40,7 @@
{
phi = phiHbyA - p_rghEqn.flux();
- U = HbyA + rAU*fvc::reconstruct((phig - p_rghEqn.flux())/rAUf);
+ U = HbyA + rAU*fvc::reconstruct((phig - p_rghEqn.flux())/Dp);
U.correctBoundaryConditions();
fvOptions.correct(U);
}
diff --git a/applications/solvers/multiphase/interPhaseChangeFoam/interPhaseChangeDyMFoam/interPhaseChangeDyMFoam.C b/applications/solvers/multiphase/interPhaseChangeFoam/interPhaseChangeDyMFoam/interPhaseChangeDyMFoam.C
index da11697426..1246ee203d 100644
--- a/applications/solvers/multiphase/interPhaseChangeFoam/interPhaseChangeDyMFoam/interPhaseChangeDyMFoam.C
+++ b/applications/solvers/multiphase/interPhaseChangeFoam/interPhaseChangeDyMFoam/interPhaseChangeDyMFoam.C
@@ -60,14 +60,13 @@ int main(int argc, char *argv[])
#include "createDynamicFvMesh.H"
#include "readGravitationalAcceleration.H"
#include "initContinuityErrs.H"
- #include "createFields.H"
- #include "readTimeControls.H"
pimpleControl pimple(mesh);
- surfaceScalarField phiAbs("phiAbs", phi);
- fvc::makeAbsolute(phiAbs, U);
-
+ #include "createFields.H"
+ #include "../interFoam/interDyMFoam/createUf.H"
+ #include "readTimeControls.H"
+ #include "../interFoam/interDyMFoam/createPcorrTypes.H"
#include "../interFoam/interDyMFoam/correctPhi.H"
#include "CourantNo.H"
#include "setInitialDeltaT.H"
@@ -88,21 +87,7 @@ int main(int argc, char *argv[])
scalar timeBeforeMeshUpdate = runTime.elapsedCpuTime();
- {
- // Ensure old-time U exists for mapping
- U.oldTime();
-
- // Calculate the relative velocity used to map the relative flux phi
- volVectorField Urel("Urel", U);
-
- if (mesh.moving())
- {
- Urel -= fvc::reconstruct(fvc::meshPhi(U));
- }
-
- // Do any mesh changes
- mesh.update();
- }
+ mesh.update();
if (mesh.changing())
{
@@ -116,7 +101,13 @@ int main(int argc, char *argv[])
if (mesh.changing() && correctPhi)
{
+ // Calculate absolute flux from the mapped surface velocity
+ phi = mesh.Sf() & Uf;
+
#include "../interFoam/interDyMFoam/correctPhi.H"
+
+ // Make the flux relative to the mesh motion
+ fvc::makeRelative(phi, U);
}
if (mesh.changing() && checkMeshCourantNo)
diff --git a/applications/solvers/multiphase/interPhaseChangeFoam/interPhaseChangeDyMFoam/pEqn.H b/applications/solvers/multiphase/interPhaseChangeFoam/interPhaseChangeDyMFoam/pEqn.H
index afb6478b9d..1169de5c54 100644
--- a/applications/solvers/multiphase/interPhaseChangeFoam/interPhaseChangeDyMFoam/pEqn.H
+++ b/applications/solvers/multiphase/interPhaseChangeFoam/interPhaseChangeDyMFoam/pEqn.H
@@ -1,6 +1,6 @@
{
volScalarField rAU("rAU", 1.0/UEqn.A());
- surfaceScalarField rAUf("Dp", fvc::interpolate(rAU));
+ surfaceScalarField Dp("Dp", fvc::interpolate(rAU));
volVectorField HbyA("HbyA", U);
HbyA = rAU*UEqn.H();
@@ -9,7 +9,7 @@
(
"phiHbyA",
(fvc::interpolate(HbyA) & mesh.Sf())
- + fvc::ddtPhiCorr(rAU, rho, U, phiAbs)
+ + fvc::interpolate(rho*rAU)*fvc::ddtCorr(U, Uf)
);
if (p_rgh.needReference())
@@ -19,14 +19,14 @@
fvc::makeAbsolute(phiHbyA, U);
}
- phiAbs = phiHbyA;
+ surfaceScalarField phiAbs("phiAbs", phiHbyA);
surfaceScalarField phig
(
(
fvc::interpolate(interface.sigmaK())*fvc::snGrad(alpha1)
- ghf*fvc::snGrad(rho)
- )*rAUf*mesh.magSf()
+ )*Dp*mesh.magSf()
);
phiHbyA += phig;
@@ -39,7 +39,7 @@
{
fvScalarMatrix p_rghEqn
(
- fvc::div(phiHbyA) - fvm::laplacian(rAUf, p_rgh)
+ fvc::div(phiHbyA) - fvm::laplacian(Dp, p_rgh)
- (vDotvP - vDotcP)*(pSat - rho*gh) + fvm::Sp(vDotvP - vDotcP, p_rgh)
);
@@ -51,13 +51,17 @@
{
phi = phiHbyA + p_rghEqn.flux();
- U = HbyA + rAU*fvc::reconstruct((phig + p_rghEqn.flux())/rAUf);
+ U = HbyA + rAU*fvc::reconstruct((phig + p_rghEqn.flux())/Dp);
U.correctBoundaryConditions();
fvOptions.correct(U);
}
}
- phiAbs = phi;
+ {
+ Uf = fvc::interpolate(U);
+ surfaceVectorField n(mesh.Sf()/mesh.magSf());
+ Uf += mesh.Sf()*(phi - (mesh.Sf() & Uf))/sqr(mesh.magSf());
+ }
// Make the fluxes relative to the mesh motion
fvc::makeRelative(phi, U);
diff --git a/applications/solvers/multiphase/interPhaseChangeFoam/pEqn.H b/applications/solvers/multiphase/interPhaseChangeFoam/pEqn.H
index 7d9b71669a..75f06420ee 100644
--- a/applications/solvers/multiphase/interPhaseChangeFoam/pEqn.H
+++ b/applications/solvers/multiphase/interPhaseChangeFoam/pEqn.H
@@ -1,6 +1,6 @@
{
volScalarField rAU("rAU", 1.0/UEqn.A());
- surfaceScalarField rAUf("Dp", fvc::interpolate(rAU));
+ surfaceScalarField Dp("Dp", fvc::interpolate(rAU));
volVectorField HbyA("HbyA", U);
HbyA = rAU*UEqn.H();
@@ -9,7 +9,7 @@
(
"phiHbyA",
(fvc::interpolate(HbyA) & mesh.Sf())
- + fvc::ddtPhiCorr(rAU, rho, U, phi)
+ + fvc::interpolate(rho*rAU)*fvc::ddtCorr(U, phi)
);
adjustPhi(phiHbyA, U, p_rgh);
phi = phiHbyA;
@@ -19,7 +19,7 @@
(
fvc::interpolate(interface.sigmaK())*fvc::snGrad(alpha1)
- ghf*fvc::snGrad(rho)
- )*rAUf*mesh.magSf()
+ )*Dp*mesh.magSf()
);
phiHbyA += phig;
@@ -32,7 +32,7 @@
{
fvScalarMatrix p_rghEqn
(
- fvc::div(phiHbyA) - fvm::laplacian(rAUf, p_rgh)
+ fvc::div(phiHbyA) - fvm::laplacian(Dp, p_rgh)
- (vDotvP - vDotcP)*(pSat - rho*gh) + fvm::Sp(vDotvP - vDotcP, p_rgh)
);
@@ -44,7 +44,7 @@
{
phi = phiHbyA + p_rghEqn.flux();
- U = HbyA + rAU*fvc::reconstruct((phig + p_rghEqn.flux())/rAUf);
+ U = HbyA + rAU*fvc::reconstruct((phig + p_rghEqn.flux())/Dp);
U.correctBoundaryConditions();
fvOptions.correct(U);
}
diff --git a/applications/solvers/multiphase/multiphaseEulerFoam/pEqn.H b/applications/solvers/multiphase/multiphaseEulerFoam/pEqn.H
index 8677157de4..00d41d3bb4 100644
--- a/applications/solvers/multiphase/multiphaseEulerFoam/pEqn.H
+++ b/applications/solvers/multiphase/multiphaseEulerFoam/pEqn.H
@@ -94,7 +94,7 @@
phiHbyAs[phasei] =
(
(fvc::interpolate(HbyAs[phasei]) & mesh.Sf())
- + fvc::ddtPhiCorr(rAUs[phasei], alpha, phase.U(), phase.phi())
+ + rAlphaAUfs[phasei]*fvc::ddtCorr(phase.U(), phase.phi())
);
mrfZones.makeRelative(phiHbyAs[phasei]);
diff --git a/applications/solvers/multiphase/multiphaseInterFoam/MRFMultiphaseInterFoam/pEqn.H b/applications/solvers/multiphase/multiphaseInterFoam/MRFMultiphaseInterFoam/pEqn.H
index e7f4c3ac8e..afa9d6058e 100644
--- a/applications/solvers/multiphase/multiphaseInterFoam/MRFMultiphaseInterFoam/pEqn.H
+++ b/applications/solvers/multiphase/multiphaseInterFoam/MRFMultiphaseInterFoam/pEqn.H
@@ -1,6 +1,6 @@
{
volScalarField rAU("rAU", 1.0/UEqn.A());
- surfaceScalarField rAUf("Dp", fvc::interpolate(rAU));
+ surfaceScalarField Dp("Dp", fvc::interpolate(rAU));
volVectorField HbyA("HbyA", U);
HbyA = rAU*UEqn.H();
@@ -9,7 +9,7 @@
(
"phiHbyA",
(fvc::interpolate(HbyA) & mesh.Sf())
- + fvc::ddtPhiCorr(rAU, rho, U, phi)
+ + fvc::interpolate(rho*rAU)*fvc::ddtCorr(U, phi)
);
adjustPhi(phiHbyA, U, p_rgh);
mrfZones.makeRelative(phiHbyA);
@@ -17,7 +17,7 @@
surfaceScalarField phig
(
- - ghf*fvc::snGrad(rho)*rAUf*mesh.magSf()
+ - ghf*fvc::snGrad(rho)*Dp*mesh.magSf()
);
phiHbyA += phig;
@@ -26,7 +26,7 @@
{
fvScalarMatrix p_rghEqn
(
- fvm::laplacian(rAUf, p_rgh) == fvc::div(phiHbyA)
+ fvm::laplacian(Dp, p_rgh) == fvc::div(phiHbyA)
);
p_rghEqn.setReference(pRefCell, getRefCellValue(p_rgh, pRefCell));
@@ -37,7 +37,7 @@
{
phi = phiHbyA - p_rghEqn.flux();
- U = HbyA + rAU*fvc::reconstruct((phig - p_rghEqn.flux())/rAUf);
+ U = HbyA + rAU*fvc::reconstruct((phig - p_rghEqn.flux())/Dp);
U.correctBoundaryConditions();
}
}
diff --git a/applications/solvers/multiphase/multiphaseInterFoam/pEqn.H b/applications/solvers/multiphase/multiphaseInterFoam/pEqn.H
index 82bf3dbd59..d092aab3e2 100644
--- a/applications/solvers/multiphase/multiphaseInterFoam/pEqn.H
+++ b/applications/solvers/multiphase/multiphaseInterFoam/pEqn.H
@@ -1,6 +1,6 @@
{
volScalarField rAU("rAU", 1.0/UEqn.A());
- surfaceScalarField rAUf("Dp", fvc::interpolate(rAU));
+ surfaceScalarField Dp("Dp", fvc::interpolate(rAU));
volVectorField HbyA("HbyA", U);
HbyA = rAU*UEqn.H();
@@ -9,7 +9,7 @@
(
"phiHbyA",
(fvc::interpolate(HbyA) & mesh.Sf())
- + fvc::ddtPhiCorr(rAU, rho, U, phi)
+ + fvc::interpolate(rho*rAU)*fvc::ddtCorr(U, phi)
);
adjustPhi(phiHbyA, U, p_rgh);
phi = phiHbyA;
@@ -19,7 +19,7 @@
(
mixture.surfaceTensionForce()
- ghf*fvc::snGrad(rho)
- )*rAUf*mesh.magSf()
+ )*Dp*mesh.magSf()
);
phiHbyA += phig;
@@ -28,7 +28,7 @@
{
fvScalarMatrix p_rghEqn
(
- fvm::laplacian(rAUf, p_rgh) == fvc::div(phiHbyA)
+ fvm::laplacian(Dp, p_rgh) == fvc::div(phiHbyA)
);
p_rghEqn.setReference(pRefCell, getRefCellValue(p_rgh, pRefCell));
@@ -39,7 +39,7 @@
{
phi = phiHbyA - p_rghEqn.flux();
- U = HbyA + rAU*fvc::reconstruct((phig - p_rghEqn.flux())/rAUf);
+ U = HbyA + rAU*fvc::reconstruct((phig - p_rghEqn.flux())/Dp);
U.correctBoundaryConditions();
}
}
diff --git a/applications/solvers/multiphase/potentialFreeSurfaceFoam/pEqn.H b/applications/solvers/multiphase/potentialFreeSurfaceFoam/pEqn.H
index 024598673f..9d5f21a94d 100644
--- a/applications/solvers/multiphase/potentialFreeSurfaceFoam/pEqn.H
+++ b/applications/solvers/multiphase/potentialFreeSurfaceFoam/pEqn.H
@@ -13,7 +13,7 @@ surfaceScalarField phiHbyA
(
"phiHbyA",
(fvc::interpolate(HbyA) & mesh.Sf())
- + fvc::ddtPhiCorr(rAU, U, phi)
+ + rAUf*fvc::ddtCorr(U, phi)
);
adjustPhi(phiHbyA, U, p_gh);
diff --git a/applications/solvers/multiphase/settlingFoam/pEqn.H b/applications/solvers/multiphase/settlingFoam/pEqn.H
index 6d61ec8730..6bc6e9a6c1 100644
--- a/applications/solvers/multiphase/settlingFoam/pEqn.H
+++ b/applications/solvers/multiphase/settlingFoam/pEqn.H
@@ -1,7 +1,6 @@
{
volScalarField rAU("rAU", 1.0/UEqn.A());
-
- surfaceScalarField rAUf("Dp", fvc::interpolate(rho*rAU));
+ surfaceScalarField Dp("Dp", fvc::interpolate(rho*rAU));
volVectorField HbyA("HbyA", U);
HbyA = rAU*UEqn.H();
@@ -9,17 +8,16 @@
surfaceScalarField phiHbyA
(
"phiHbyA",
- fvc::interpolate(rho)
- *(
- (fvc::interpolate(HbyA) & mesh.Sf())
- + fvc::ddtPhiCorr(rAU, rho, U, phi)
+ (
+ (fvc::interpolate(rho*HbyA) & mesh.Sf())
+ + Dp*fvc::ddtCorr(rho, U, phi)
)
);
phi = phiHbyA;
surfaceScalarField phig
(
- - ghf*fvc::snGrad(rho)*rAUf*mesh.magSf()
+ - ghf*fvc::snGrad(rho)*Dp*mesh.magSf()
);
phiHbyA += phig;
@@ -28,7 +26,7 @@
{
fvScalarMatrix p_rghEqn
(
- fvm::laplacian(rAUf, p_rgh) == fvc::ddt(rho) + fvc::div(phiHbyA)
+ fvm::laplacian(Dp, p_rgh) == fvc::ddt(rho) + fvc::div(phiHbyA)
);
p_rghEqn.setReference(pRefCell, getRefCellValue(p_rgh, pRefCell));
@@ -39,7 +37,7 @@
{
phi = phiHbyA - p_rghEqn.flux();
- U = HbyA + rAU*fvc::reconstruct((phig - p_rghEqn.flux())/rAUf);
+ U = HbyA + rAU*fvc::reconstruct((phig - p_rghEqn.flux())/Dp);
U.correctBoundaryConditions();
}
}
diff --git a/applications/solvers/multiphase/twoLiquidMixingFoam/pEqn.H b/applications/solvers/multiphase/twoLiquidMixingFoam/pEqn.H
index e090522de6..94495f4f30 100644
--- a/applications/solvers/multiphase/twoLiquidMixingFoam/pEqn.H
+++ b/applications/solvers/multiphase/twoLiquidMixingFoam/pEqn.H
@@ -1,6 +1,6 @@
{
volScalarField rAU("rAU", 1.0/UEqn.A());
- surfaceScalarField rAUf("Dp", fvc::interpolate(rAU));
+ surfaceScalarField Dp("Dp", fvc::interpolate(rAU));
volVectorField HbyA("HbyA", U);
HbyA = rAU*UEqn.H();
@@ -9,14 +9,14 @@
(
"phiHbyA",
(fvc::interpolate(HbyA) & mesh.Sf())
- + fvc::ddtPhiCorr(rAU, rho, U, phi)
+ + fvc::interpolate(rho*rAU)*fvc::ddtCorr(U, phi)
);
adjustPhi(phiHbyA, U, p_rgh);
phi = phiHbyA;
surfaceScalarField phig
(
- - ghf*fvc::snGrad(rho)*rAUf*mesh.magSf()
+ - ghf*fvc::snGrad(rho)*Dp*mesh.magSf()
);
phiHbyA += phig;
@@ -25,7 +25,7 @@
{
fvScalarMatrix p_rghEqn
(
- fvm::laplacian(rAUf, p_rgh) == fvc::div(phiHbyA)
+ fvm::laplacian(Dp, p_rgh) == fvc::div(phiHbyA)
);
p_rghEqn.setReference(pRefCell, getRefCellValue(p_rgh, pRefCell));
@@ -36,7 +36,7 @@
{
phi = phiHbyA - p_rghEqn.flux();
- U = HbyA + rAU*fvc::reconstruct((phig - p_rghEqn.flux())/rAUf);
+ U = HbyA + rAU*fvc::reconstruct((phig - p_rghEqn.flux())/Dp);
U.correctBoundaryConditions();
}
}
diff --git a/applications/solvers/multiphase/twoPhaseEulerFoam/pEqn.H b/applications/solvers/multiphase/twoPhaseEulerFoam/pEqn.H
index 21a38ea5a6..6485ee0323 100644
--- a/applications/solvers/multiphase/twoPhaseEulerFoam/pEqn.H
+++ b/applications/solvers/multiphase/twoPhaseEulerFoam/pEqn.H
@@ -47,14 +47,14 @@
(
IOobject::groupName("phiHbyA", phase1.name()),
(fvc::interpolate(HbyA1) & mesh.Sf())
- + fvc::ddtPhiCorr(rAU1, alpha1, U1, phi1)
+ + rAlphaAU1f*fvc::ddtCorr(U1, phi1)
);
surfaceScalarField phiHbyA2
(
IOobject::groupName("phiHbyA", phase2.name()),
(fvc::interpolate(HbyA2) & mesh.Sf())
- + fvc::ddtPhiCorr(rAU2, alpha2, U2, phi2)
+ + rAlphaAU2f*fvc::ddtCorr(U2, phi2)
);
phiHbyA1 +=
diff --git a/applications/test/LduMatrix/LduMatrixTest3.C b/applications/test/LduMatrix/LduMatrixTest3.C
index 596559e6a7..e9cce12310 100644
--- a/applications/test/LduMatrix/LduMatrixTest3.C
+++ b/applications/test/LduMatrix/LduMatrixTest3.C
@@ -105,7 +105,7 @@ int main(int argc, char *argv[])
U = rAU*UEqn.H();
phi = (fvc::interpolate(U) & mesh.Sf())
- + fvc::ddtPhiCorr(rAU, U, phi);
+ + fvc::ddtCorr(rAU, U, phi);
adjustPhi(phi, U, p);
diff --git a/applications/test/PisoFoam/PisoFoam.C b/applications/test/PisoFoam/PisoFoam.C
index bbd28e6f88..ecf0df2f50 100644
--- a/applications/test/PisoFoam/PisoFoam.C
+++ b/applications/test/PisoFoam/PisoFoam.C
@@ -87,7 +87,7 @@ int main(int argc, char *argv[])
(
"phiHbyA",
(fvc::interpolate(HbyA) & mesh.Sf())
- + fvc::ddtPhiCorr(rAU, U, phi)
+ + fvc::ddtCorr(rAU, U, phi)
);
adjustPhi(phiHbyA, U, p);
diff --git a/applications/test/RhoPimpleFoam/pEqn.H b/applications/test/RhoPimpleFoam/pEqn.H
index fcc40881a1..465e81d0e9 100644
--- a/applications/test/RhoPimpleFoam/pEqn.H
+++ b/applications/test/RhoPimpleFoam/pEqn.H
@@ -20,7 +20,7 @@ if (pimple.transonic())
fvc::interpolate(psi)
*(
(fvc::interpolate(HbyA) & mesh.Sf())
- + fvc::ddtPhiCorr(rAU, rho, U, phi)
+ + fvc::ddtCorr(rAU, rho, U, phi)
)
);
@@ -57,7 +57,7 @@ else
fvc::interpolate(rho)
*(
(fvc::interpolate(HbyA) & mesh.Sf())
- + fvc::ddtPhiCorr(rAU, rho, U, phi)
+ + fvc::ddtCorr(rAU, rho, U, phi)
)
);
diff --git a/src/finiteVolume/finiteVolume/ddtSchemes/CoEulerDdtScheme/CoEulerDdtScheme.C b/src/finiteVolume/finiteVolume/ddtSchemes/CoEulerDdtScheme/CoEulerDdtScheme.C
index d6ef3989cd..fab1e295d3 100644
--- a/src/finiteVolume/finiteVolume/ddtSchemes/CoEulerDdtScheme/CoEulerDdtScheme.C
+++ b/src/finiteVolume/finiteVolume/ddtSchemes/CoEulerDdtScheme/CoEulerDdtScheme.C
@@ -478,55 +478,163 @@ CoEulerDdtScheme::fvmDdt
}
+template
+tmp::fluxFieldType>
+CoEulerDdtScheme::fvcDdtUfCorr
+(
+ const GeometricField& U,
+ const GeometricField& Uf
+)
+{
+ IOobject ddtIOobject
+ (
+ "ddtCorr(" + U.name() + ',' + Uf.name() + ')',
+ mesh().time().timeName(),
+ mesh()
+ );
+
+ const surfaceScalarField rDeltaT(fvc::interpolate(CorDeltaT()));
+
+ fluxFieldType phiCorr
+ (
+ mesh().Sf() & (Uf.oldTime() - fvc::interpolate(U.oldTime()))
+ );
+
+ return tmp
+ (
+ new fluxFieldType
+ (
+ ddtIOobject,
+ this->fvcDdtPhiCoeff
+ (
+ U.oldTime(),
+ (mesh().Sf() & Uf.oldTime()),
+ phiCorr
+ )
+ *rDeltaT*phiCorr
+ )
+ );
+}
+
+
template
tmp::fluxFieldType>
CoEulerDdtScheme::fvcDdtPhiCorr
(
- const volScalarField& rA,
const GeometricField& U,
const fluxFieldType& phi
)
{
IOobject ddtIOobject
(
- "ddtPhiCorr(" + rA.name() + ',' + U.name() + ',' + phi.name() + ')',
+ "ddtCorr(" + U.name() + ',' + phi.name() + ')',
mesh().time().timeName(),
mesh()
);
- if (mesh().moving())
+ const surfaceScalarField rDeltaT(fvc::interpolate(CorDeltaT()));
+
+ fluxFieldType phiCorr
+ (
+ phi.oldTime() - (mesh().Sf() & fvc::interpolate(U.oldTime()))
+ );
+
+ return tmp
+ (
+ new fluxFieldType
+ (
+ ddtIOobject,
+ this->fvcDdtPhiCoeff(U.oldTime(), phi.oldTime(), phiCorr)
+ *rDeltaT*phiCorr
+ )
+ );
+}
+
+
+template
+tmp::fluxFieldType>
+CoEulerDdtScheme::fvcDdtUfCorr
+(
+ const volScalarField& rho,
+ const GeometricField& U,
+ const GeometricField& Uf
+)
+{
+ IOobject ddtIOobject
+ (
+ "ddtCorr(" + rho.name() + ',' + U.name() + ',' + Uf.name() + ')',
+ mesh().time().timeName(),
+ mesh()
+ );
+
+ const surfaceScalarField rDeltaT(fvc::interpolate(CorDeltaT()));
+
+ if
+ (
+ U.dimensions() == dimVelocity
+ && Uf.dimensions() == dimDensity*dimVelocity
+ )
{
+ GeometricField rhoU0
+ (
+ rho.oldTime()*U.oldTime()
+ );
+
+ fluxFieldType phiCorr
+ (
+ mesh().Sf() & (Uf.oldTime() - fvc::interpolate(rhoU0))
+ );
+
return tmp
(
new fluxFieldType
(
ddtIOobject,
- mesh(),
- dimensioned::type>
+ this->fvcDdtPhiCoeff
(
- "0",
- rA.dimensions()*phi.dimensions()/dimTime,
- pTraits::type>::zero
+ rhoU0,
+ mesh().Sf() & Uf.oldTime(),
+ phiCorr
)
+ *rDeltaT*phiCorr
+ )
+ );
+ }
+ else if
+ (
+ U.dimensions() == dimDensity*dimVelocity
+ && Uf.dimensions() == dimDensity*dimVelocity
+ )
+ {
+ fluxFieldType phiCorr
+ (
+ mesh().Sf() & (Uf.oldTime() - fvc::interpolate(U.oldTime()))
+ );
+
+ return tmp
+ (
+ new fluxFieldType
+ (
+ ddtIOobject,
+ this->fvcDdtPhiCoeff
+ (
+ U.oldTime(),
+ mesh().Sf() & Uf.oldTime(),
+ phiCorr
+ )
+ *rDeltaT*phiCorr
)
);
}
else
{
- const volScalarField rDeltaT(CorDeltaT());
-
- return tmp
+ FatalErrorIn
(
- new fluxFieldType
- (
- ddtIOobject,
- this->fvcDdtPhiCoeff(U.oldTime(), phi.oldTime())*
- (
- fvc::interpolate(rDeltaT*rA)*phi.oldTime()
- - (fvc::interpolate(rDeltaT*rA*U.oldTime()) & mesh().Sf())
- )
- )
- );
+ "CoEulerDdtScheme::fvcDdtPhiCorr"
+ ) << "dimensions of Uf are not correct"
+ << abort(FatalError);
+
+ return fluxFieldType::null();
}
}
@@ -535,122 +643,76 @@ template
tmp::fluxFieldType>
CoEulerDdtScheme::fvcDdtPhiCorr
(
- const volScalarField& rA,
const volScalarField& rho,
const GeometricField& U,
const fluxFieldType& phi
)
{
+ dimensionedScalar rDeltaT = 1.0/mesh().time().deltaT();
+
IOobject ddtIOobject
(
- "ddtPhiCorr("
- + rA.name() + ',' + rho.name() + ',' + U.name() + ',' + phi.name() + ')',
+ "ddtCorr(" + rho.name() + ',' + U.name() + ',' + phi.name() + ')',
mesh().time().timeName(),
mesh()
);
- if (mesh().moving())
+ if
+ (
+ U.dimensions() == dimVelocity
+ && phi.dimensions() == rho.dimensions()*dimVelocity*dimArea
+ )
{
+ GeometricField rhoU0
+ (
+ rho.oldTime()*U.oldTime()
+ );
+
+ fluxFieldType phiCorr
+ (
+ phi.oldTime() - (mesh().Sf() & fvc::interpolate(rhoU0))
+ );
+
return tmp
(
new fluxFieldType
(
ddtIOobject,
- mesh(),
- dimensioned::type>
- (
- "0",
- rA.dimensions()*rho.dimensions()*phi.dimensions()/dimTime,
- pTraits::type>::zero
- )
+ this->fvcDdtPhiCoeff(rhoU0, phi.oldTime(), phiCorr)
+ *rDeltaT*phiCorr
+ )
+ );
+ }
+ else if
+ (
+ U.dimensions() == rho.dimensions()*dimVelocity
+ && phi.dimensions() == rho.dimensions()*dimVelocity*dimArea
+ )
+ {
+ fluxFieldType phiCorr
+ (
+ phi.oldTime() - (mesh().Sf() & fvc::interpolate(U.oldTime()))
+ );
+
+ return tmp
+ (
+ new fluxFieldType
+ (
+ ddtIOobject,
+ this->fvcDdtPhiCoeff(U.oldTime(), phi.oldTime(), phiCorr)
+ *rDeltaT*phiCorr
)
);
}
else
{
- const volScalarField rDeltaT(CorDeltaT());
+ FatalErrorIn
+ (
+ "CoEulerDdtScheme::fvcDdtPhiCorr"
+ ) << "dimensions of phi are not correct"
+ << abort(FatalError);
- if
- (
- U.dimensions() == dimVelocity
- && phi.dimensions() == dimVelocity*dimArea
- )
- {
- return tmp
- (
- new fluxFieldType
- (
- ddtIOobject,
- this->fvcDdtPhiCoeff(U.oldTime(), phi.oldTime())
- *(
- fvc::interpolate(rDeltaT*rA*rho.oldTime())*phi.oldTime()
- - (fvc::interpolate(rDeltaT*rA*rho.oldTime()*U.oldTime())
- & mesh().Sf())
- )
- )
- );
- }
- else if
- (
- U.dimensions() == dimVelocity
- && phi.dimensions() == dimDensity*dimVelocity*dimArea
- )
- {
- return tmp
- (
- new fluxFieldType
- (
- ddtIOobject,
- this->fvcDdtPhiCoeff
- (
- U.oldTime(),
- phi.oldTime()/fvc::interpolate(rho.oldTime())
- )
- *(
- fvc::interpolate(rDeltaT*rA*rho.oldTime())
- *phi.oldTime()/fvc::interpolate(rho.oldTime())
- - (
- fvc::interpolate
- (
- rDeltaT*rA*rho.oldTime()*U.oldTime()
- ) & mesh().Sf()
- )
- )
- )
- );
- }
- else if
- (
- U.dimensions() == dimDensity*dimVelocity
- && phi.dimensions() == dimDensity*dimVelocity*dimArea
- )
- {
- return tmp
- (
- new fluxFieldType
- (
- ddtIOobject,
- this->fvcDdtPhiCoeff
- (rho.oldTime(), U.oldTime(), phi.oldTime())
- * (
- fvc::interpolate(rDeltaT*rA)*phi.oldTime()
- - (
- fvc::interpolate(rDeltaT*rA*U.oldTime())&mesh().Sf()
- )
- )
- )
- );
- }
- else
- {
- FatalErrorIn
- (
- "CoEulerDdtScheme::fvcDdtPhiCorr"
- ) << "dimensions of phi are not correct"
- << abort(FatalError);
-
- return fluxFieldType::null();
- }
+ return fluxFieldType::null();
}
}
diff --git a/src/finiteVolume/finiteVolume/ddtSchemes/CoEulerDdtScheme/CoEulerDdtScheme.H b/src/finiteVolume/finiteVolume/ddtSchemes/CoEulerDdtScheme/CoEulerDdtScheme.H
index 52401609fa..a59a780967 100644
--- a/src/finiteVolume/finiteVolume/ddtSchemes/CoEulerDdtScheme/CoEulerDdtScheme.H
+++ b/src/finiteVolume/finiteVolume/ddtSchemes/CoEulerDdtScheme/CoEulerDdtScheme.H
@@ -158,16 +158,27 @@ public:
typedef typename ddtScheme::fluxFieldType fluxFieldType;
- tmp fvcDdtPhiCorr
+ tmp fvcDdtUfCorr
(
- const volScalarField& rA,
const GeometricField& U,
- const fluxFieldType& phi
+ const GeometricField& Uf
+ );
+
+ tmp fvcDdtPhiCorr
+ (
+ const GeometricField& U,
+ const fluxFieldType& phi
+ );
+
+ tmp fvcDdtUfCorr
+ (
+ const volScalarField& rho,
+ const GeometricField& U,
+ const GeometricField& Uf
);
tmp fvcDdtPhiCorr
(
- const volScalarField& rA,
const volScalarField& rho,
const GeometricField& U,
const fluxFieldType& phi
@@ -180,10 +191,17 @@ public:
};
+template<>
+tmp CoEulerDdtScheme::fvcDdtUfCorr
+(
+ const GeometricField& U,
+ const GeometricField& Uf
+);
+
+
template<>
tmp CoEulerDdtScheme::fvcDdtPhiCorr
(
- const volScalarField& rA,
const volScalarField& U,
const surfaceScalarField& phi
);
@@ -192,7 +210,6 @@ tmp CoEulerDdtScheme::fvcDdtPhiCorr
template<>
tmp CoEulerDdtScheme::fvcDdtPhiCorr
(
- const volScalarField& rA,
const volScalarField& rho,
const volScalarField& U,
const surfaceScalarField& phi
diff --git a/src/finiteVolume/finiteVolume/ddtSchemes/CrankNicolsonDdtScheme/CrankNicolsonDdtScheme.C b/src/finiteVolume/finiteVolume/ddtSchemes/CrankNicolsonDdtScheme/CrankNicolsonDdtScheme.C
index 53d106a75d..c0622c4635 100644
--- a/src/finiteVolume/finiteVolume/ddtSchemes/CrankNicolsonDdtScheme/CrankNicolsonDdtScheme.C
+++ b/src/finiteVolume/finiteVolume/ddtSchemes/CrankNicolsonDdtScheme/CrankNicolsonDdtScheme.C
@@ -646,7 +646,6 @@ CrankNicolsonDdtScheme::fvmDdt
fvMatrix& fvm = tfvm();
scalar rDtCoef = rDtCoef_(ddt0).value();
-
fvm.diag() = rDtCoef*mesh().V();
vf.oldTime().oldTime();
@@ -875,90 +874,197 @@ CrankNicolsonDdtScheme::fvmDdt
}
-
template
tmp::fluxFieldType>
-CrankNicolsonDdtScheme::fvcDdtPhiCorr
+CrankNicolsonDdtScheme::fvcDdtUfCorr
(
- const volScalarField& rA,
const GeometricField& U,
- const fluxFieldType& phi
+ const GeometricField& Uf
)
{
- DDt0Field >& dUdt0 =
+ DDt0Field >& ddt0 =
ddt0_ >
(
"ddt0(" + U.name() + ')',
U.dimensions()
);
- DDt0Field& dphidt0 =
- ddt0_
- (
- "ddt0(" + phi.name() + ')',
- phi.dimensions()
- );
+ dimensionedScalar rDeltaT = rDtCoef_(ddt0);
IOobject ddtIOobject
(
- "ddtPhiCorr(" + rA.name() + ',' + U.name() + ',' + phi.name() + ')',
+ "ddtCorr(" + U.name() + ',' + Uf.name() + ')',
mesh().time().timeName(),
mesh()
);
- dimensionedScalar rDtCoef = rDtCoef_(dUdt0);
+ fluxFieldType phiCorr
+ (
+ mesh().Sf() & (Uf.oldTime() - fvc::interpolate(U.oldTime()))
+ );
- if (mesh().moving())
+ return tmp
+ (
+ new fluxFieldType
+ (
+ ddtIOobject,
+ this->fvcDdtPhiCoeff
+ (
+ U.oldTime(),
+ mesh().Sf() & Uf.oldTime(),
+ phiCorr
+ )
+ *rDeltaT*phiCorr
+ )
+ );
+}
+
+
+template
+tmp::fluxFieldType>
+CrankNicolsonDdtScheme::fvcDdtPhiCorr
+(
+ const GeometricField& U,
+ const fluxFieldType& phi
+)
+{
+ DDt0Field >& ddt0 =
+ ddt0_ >
+ (
+ "ddt0(" + U.name() + ')',
+ U.dimensions()
+ );
+
+ dimensionedScalar rDeltaT = rDtCoef_(ddt0);
+
+ IOobject ddtIOobject
+ (
+ "ddtCorr(" + U.name() + ',' + phi.name() + ')',
+ mesh().time().timeName(),
+ mesh()
+ );
+
+ fluxFieldType phiCorr
+ (
+ phi.oldTime() - (mesh().Sf() & fvc::interpolate(U.oldTime()))
+ );
+
+ return tmp
+ (
+ new fluxFieldType
+ (
+ ddtIOobject,
+ this->fvcDdtPhiCoeff(U.oldTime(), phi.oldTime(), phiCorr)
+ *rDeltaT*phiCorr
+ )
+ );
+}
+
+
+template
+tmp::fluxFieldType>
+CrankNicolsonDdtScheme::fvcDdtUfCorr
+(
+ const volScalarField& rho,
+ const GeometricField& U,
+ const GeometricField& Uf
+)
+{
+ IOobject ddtIOobject
+ (
+ "ddtCorr(" + rho.name() + ',' + U.name() + ',' + Uf.name() + ')',
+ mesh().time().timeName(),
+ mesh()
+ );
+
+ if
+ (
+ U.dimensions() == dimVelocity
+ && Uf.dimensions() == rho.dimensions()*dimVelocity
+ )
{
- return tmp
+ DDt0Field >& ddt0 =
+ ddt0_ >
+ (
+ "ddt0(" + rho.name() + ',' + U.name() + ')',
+ U.dimensions()
+ );
+
+ dimensionedScalar rDeltaT = rDtCoef_(ddt0);
+
+ GeometricField rhoU0
+ (
+ rho.oldTime()*U.oldTime()
+ );
+
+ fluxFieldType phiCorr
+ (
+ mesh().Sf() & (Uf.oldTime() - fvc::interpolate(rhoU0))
+ );
+
+ tmp ddtCorr
(
new fluxFieldType
(
ddtIOobject,
- mesh(),
- dimensioned::type>
+ this->fvcDdtPhiCoeff
(
- "0",
- rDtCoef.dimensions()*rA.dimensions()*phi.dimensions(),
- pTraits::type>::zero
+ rhoU0,
+ mesh().Sf() & Uf.oldTime(),
+ phiCorr
)
+ *rDeltaT*phiCorr
)
);
+
+ return ddtCorr;
+ }
+ else if
+ (
+ U.dimensions() == rho.dimensions()*dimVelocity
+ && Uf.dimensions() == rho.dimensions()*dimVelocity
+ )
+ {
+ DDt0Field >& ddt0 =
+ ddt0_ >
+ (
+ "ddt0(" + U.name() + ')',
+ U.dimensions()
+ );
+
+ dimensionedScalar rDeltaT = rDtCoef_(ddt0);
+
+ fluxFieldType phiCorr
+ (
+ mesh().Sf() & (Uf.oldTime() - fvc::interpolate(U.oldTime()))
+ );
+
+ tmp ddtCorr
+ (
+ new fluxFieldType
+ (
+ ddtIOobject,
+ this->fvcDdtPhiCoeff
+ (
+ U.oldTime(),
+ mesh().Sf() & Uf.oldTime(),
+ phiCorr
+ )
+ *rDeltaT*phiCorr
+ )
+ );
+
+ return ddtCorr;
}
else
{
- if (evaluate(dUdt0))
- {
- dUdt0 =
- rDtCoef0_(dUdt0)*(U.oldTime() - U.oldTime().oldTime())
- - offCentre_(dUdt0());
- }
-
- if (evaluate(dphidt0))
- {
- dphidt0 =
- rDtCoef0_(dphidt0)*(phi.oldTime() - phi.oldTime().oldTime())
- - offCentre_(dphidt0());
- }
-
- return tmp
+ FatalErrorIn
(
- new fluxFieldType
- (
- ddtIOobject,
- this->fvcDdtPhiCoeff(U.oldTime(), phi.oldTime())
- *fvc::interpolate(rA)
- *(
- (rDtCoef*phi.oldTime() + offCentre_(dphidt0()))
- - (
- fvc::interpolate
- (
- (rDtCoef*U.oldTime() + offCentre_(dUdt0()))
- ) & mesh().Sf()
- )
- )
- )
- );
+ "CrankNicolsonDdtScheme::fvcDdtPhiCorr"
+ ) << "dimensions of Uf are not correct"
+ << abort(FatalError);
+
+ return fluxFieldType::null();
}
}
@@ -967,203 +1073,96 @@ template
tmp::fluxFieldType>
CrankNicolsonDdtScheme::fvcDdtPhiCorr
(
- const volScalarField& rA,
const volScalarField& rho,
const GeometricField& U,
const fluxFieldType& phi
)
{
- DDt0Field >& dUdt0 =
- ddt0_ >
- (
- "ddt0(" + U.name() + ')',
- U.dimensions()
- );
-
- DDt0Field& dphidt0 =
- ddt0_
- (
- "ddt0(" + phi.name() + ')',
- U.dimensions()*dimArea
- );
-
IOobject ddtIOobject
(
- "ddtPhiCorr("
- + rA.name() + ',' + rho.name() + ',' + U.name() + ',' + phi.name() + ')',
+ "ddtCorr(" + rho.name() + ',' + U.name() + ',' + phi.name() + ')',
mesh().time().timeName(),
mesh()
);
- dimensionedScalar rDtCoef = rDtCoef_(dUdt0);
-
- if (mesh().moving())
+ if
+ (
+ U.dimensions() == dimVelocity
+ && phi.dimensions() == rho.dimensions()*dimVelocity*dimArea
+ )
{
- return tmp
+ DDt0Field >& ddt0 =
+ ddt0_ >
+ (
+ "ddt0(" + rho.name() + ',' + U.name() + ')',
+ U.dimensions()
+ );
+
+ dimensionedScalar rDeltaT = rDtCoef_(ddt0);
+
+ GeometricField rhoU0
+ (
+ rho.oldTime()*U.oldTime()
+ );
+
+ fluxFieldType phiCorr
+ (
+ phi.oldTime() - (mesh().Sf() & fvc::interpolate(rhoU0))
+ );
+
+ tmp ddtCorr
(
new fluxFieldType
(
ddtIOobject,
- mesh(),
- dimensioned::type>
- (
- "0",
- rA.dimensions()*rho.dimensions()*phi.dimensions()/dimTime,
- pTraits::type>::zero
- )
+ this->fvcDdtPhiCoeff(rhoU0, phi.oldTime(), phiCorr)
+ *rDeltaT*phiCorr
)
);
+
+ return ddtCorr;
+ }
+ else if
+ (
+ U.dimensions() == rho.dimensions()*dimVelocity
+ && phi.dimensions() == rho.dimensions()*dimVelocity*dimArea
+ )
+ {
+ DDt0Field >& ddt0 =
+ ddt0_ >
+ (
+ "ddt0(" + U.name() + ')',
+ U.dimensions()
+ );
+
+ dimensionedScalar rDeltaT = rDtCoef_(ddt0);
+
+ fluxFieldType phiCorr
+ (
+ phi.oldTime() - (mesh().Sf() & fvc::interpolate(U.oldTime()))
+ );
+
+ tmp ddtCorr
+ (
+ new fluxFieldType
+ (
+ ddtIOobject,
+ this->fvcDdtPhiCoeff(U.oldTime(), phi.oldTime(), phiCorr)
+ *rDeltaT*phiCorr
+ )
+ );
+
+ return ddtCorr;
}
else
{
- if
+ FatalErrorIn
(
- U.dimensions() == dimVelocity
- && phi.dimensions() == dimVelocity*dimArea
- )
- {
- if (evaluate(dUdt0))
- {
- dUdt0 = rDtCoef0_(dUdt0)*
- (
- U.oldTime() - U.oldTime().oldTime()
- ) - offCentre_(dUdt0());
- }
+ "CrankNicolsonDdtScheme::fvcDdtPhiCorr"
+ ) << "dimensions of phi are not correct"
+ << abort(FatalError);
- if (evaluate(dphidt0))
- {
- dphidt0 = rDtCoef0_(dphidt0)*
- (
- phi.oldTime()
- - fvc::interpolate(rho.oldTime().oldTime()/rho.oldTime())
- *phi.oldTime().oldTime()
- ) - offCentre_(dphidt0());
- }
-
- return tmp
- (
- new fluxFieldType
- (
- ddtIOobject,
- this->fvcDdtPhiCoeff(U.oldTime(), phi.oldTime())*
- (
- fvc::interpolate(rA*rho.oldTime())
- *(rDtCoef*phi.oldTime() + offCentre_(dphidt0()))
- - (
- fvc::interpolate
- (
- rA*rho.oldTime()
- *(rDtCoef*U.oldTime() + offCentre_(dUdt0()))
- ) & mesh().Sf()
- )
- )
- )
- );
- }
- else if
- (
- U.dimensions() == dimVelocity
- && phi.dimensions() == rho.dimensions()*dimVelocity*dimArea
- )
- {
- if (evaluate(dUdt0))
- {
- dUdt0 = rDtCoef0_(dUdt0)*
- (
- U.oldTime() - U.oldTime().oldTime()
- ) - offCentre_(dUdt0());
- }
-
- if (evaluate(dphidt0))
- {
- dphidt0 = rDtCoef0_(dphidt0)*
- (
- phi.oldTime()
- /fvc::interpolate(rho.oldTime())
- - phi.oldTime().oldTime()
- /fvc::interpolate(rho.oldTime().oldTime())
- ) - offCentre_(dphidt0());
- }
-
- return tmp
- (
- new fluxFieldType
- (
- ddtIOobject,
- this->fvcDdtPhiCoeff
- (
- U.oldTime(),
- phi.oldTime()/fvc::interpolate(rho.oldTime())
- )
- *(
- fvc::interpolate(rA*rho.oldTime())
- *(
- rDtCoef*phi.oldTime()/fvc::interpolate(rho.oldTime())
- + offCentre_(dphidt0())
- )
- - (
- fvc::interpolate
- (
- rA*rho.oldTime()
- *(rDtCoef*U.oldTime() + offCentre_(dUdt0()))
- ) & mesh().Sf()
- )
- )
- )
- );
- }
- else if
- (
- U.dimensions() == rho.dimensions()*dimVelocity
- && phi.dimensions() == rho.dimensions()*dimVelocity*dimArea
- )
- {
- if (evaluate(dUdt0))
- {
- dUdt0 = rDtCoef0_(dUdt0)*
- (
- U.oldTime() - U.oldTime().oldTime()
- ) - offCentre_(dUdt0());
- }
-
- if (evaluate(dphidt0))
- {
- dphidt0 = rDtCoef0_(dphidt0)*
- (
- phi.oldTime() - phi.oldTime().oldTime()
- ) - offCentre_(dphidt0());
- }
-
- return tmp
- (
- new fluxFieldType
- (
- ddtIOobject,
- this->fvcDdtPhiCoeff
- (rho.oldTime(), U.oldTime(), phi.oldTime())
- * (
- fvc::interpolate(rA)
- *(rDtCoef*phi.oldTime() + offCentre_(dphidt0()))
- - (
- fvc::interpolate
- (
- rA*(rDtCoef*U.oldTime() + offCentre_(dUdt0()))
- ) & mesh().Sf()
- )
- )
- )
- );
- }
- else
- {
- FatalErrorIn
- (
- "CrankNicolsonDdtScheme::fvcDdtPhiCorr"
- ) << "dimensions of phi are not correct"
- << abort(FatalError);
-
- return fluxFieldType::null();
- }
+ return fluxFieldType::null();
}
}
diff --git a/src/finiteVolume/finiteVolume/ddtSchemes/CrankNicolsonDdtScheme/CrankNicolsonDdtScheme.H b/src/finiteVolume/finiteVolume/ddtSchemes/CrankNicolsonDdtScheme/CrankNicolsonDdtScheme.H
index 315072306c..cee3d96fbd 100644
--- a/src/finiteVolume/finiteVolume/ddtSchemes/CrankNicolsonDdtScheme/CrankNicolsonDdtScheme.H
+++ b/src/finiteVolume/finiteVolume/ddtSchemes/CrankNicolsonDdtScheme/CrankNicolsonDdtScheme.H
@@ -229,16 +229,27 @@ public:
typedef typename ddtScheme::fluxFieldType fluxFieldType;
- tmp fvcDdtPhiCorr
+ tmp fvcDdtUfCorr
(
- const volScalarField& rA,
const GeometricField& U,
- const fluxFieldType& phi
+ const GeometricField& Uf
+ );
+
+ tmp fvcDdtPhiCorr
+ (
+ const GeometricField& U,
+ const fluxFieldType& phi
+ );
+
+ tmp fvcDdtUfCorr
+ (
+ const volScalarField& rho,
+ const GeometricField& U,
+ const GeometricField& Uf
);
tmp fvcDdtPhiCorr
(
- const volScalarField& rA,
const volScalarField& rho,
const GeometricField& U,
const fluxFieldType& phi
@@ -252,10 +263,17 @@ public:
};
+template<>
+tmp CrankNicolsonDdtScheme::fvcDdtUfCorr
+(
+ const GeometricField& U,
+ const GeometricField& Uf
+);
+
+
template<>
tmp CrankNicolsonDdtScheme::fvcDdtPhiCorr
(
- const volScalarField& rA,
const volScalarField& U,
const surfaceScalarField& phi
);
@@ -264,7 +282,6 @@ tmp CrankNicolsonDdtScheme::fvcDdtPhiCorr
template<>
tmp CrankNicolsonDdtScheme::fvcDdtPhiCorr
(
- const volScalarField& rA,
const volScalarField& rho,
const volScalarField& U,
const surfaceScalarField& phi
diff --git a/src/finiteVolume/finiteVolume/ddtSchemes/EulerDdtScheme/EulerDdtScheme.C b/src/finiteVolume/finiteVolume/ddtSchemes/EulerDdtScheme/EulerDdtScheme.C
index d48a634491..ba137df9e6 100644
--- a/src/finiteVolume/finiteVolume/ddtSchemes/EulerDdtScheme/EulerDdtScheme.C
+++ b/src/finiteVolume/finiteVolume/ddtSchemes/EulerDdtScheme/EulerDdtScheme.C
@@ -370,34 +370,38 @@ EulerDdtScheme::fvmDdt
template
tmp::fluxFieldType>
-EulerDdtScheme::fvcDdtPhiCorr
+EulerDdtScheme::fvcDdtUfCorr
(
- const volScalarField& rA,
const GeometricField& U,
- const fluxFieldType& phiAbs
+ const GeometricField& Uf
)
{
dimensionedScalar rDeltaT = 1.0/mesh().time().deltaT();
IOobject ddtIOobject
(
- "ddtPhiCorr(" + rA.name() + ',' + U.name() + ',' + phiAbs.name() + ')',
+ "ddtCorr(" + U.name() + ',' + Uf.name() + ')',
mesh().time().timeName(),
mesh()
);
- tmp phiCorr =
- phiAbs.oldTime() - (fvc::interpolate(U.oldTime()) & mesh().Sf());
-
- phiCorr().boundaryField() = pTraits::type>::zero;
+ fluxFieldType phiCorr
+ (
+ mesh().Sf() & (Uf.oldTime() - fvc::interpolate(U.oldTime()))
+ );
return tmp
(
new fluxFieldType
(
ddtIOobject,
- this->fvcDdtPhiCoeff(U.oldTime(), phiAbs.oldTime(), phiCorr())
- * fvc::interpolate(rDeltaT*rA)*phiCorr
+ this->fvcDdtPhiCoeff
+ (
+ U.oldTime(),
+ mesh().Sf() & Uf.oldTime(),
+ phiCorr
+ )
+ *rDeltaT*phiCorr
)
);
}
@@ -407,21 +411,50 @@ template
tmp::fluxFieldType>
EulerDdtScheme::fvcDdtPhiCorr
(
- const volScalarField& rA,
- const volScalarField& rho,
const GeometricField& U,
- const fluxFieldType& phiAbs
+ const fluxFieldType& phi
)
{
dimensionedScalar rDeltaT = 1.0/mesh().time().deltaT();
IOobject ddtIOobject
(
- "ddtPhiCorr("
- + rA.name() + ','
- + rho.name() + ','
- + U.name() + ','
- + phiAbs.name() + ')',
+ "ddtCorr(" + U.name() + ',' + phi.name() + ')',
+ mesh().time().timeName(),
+ mesh()
+ );
+
+ fluxFieldType phiCorr
+ (
+ phi.oldTime() - (mesh().Sf() & fvc::interpolate(U.oldTime()))
+ );
+
+ return tmp
+ (
+ new fluxFieldType
+ (
+ ddtIOobject,
+ this->fvcDdtPhiCoeff(U.oldTime(), phi.oldTime(), phiCorr)
+ *rDeltaT*phiCorr
+ )
+ );
+}
+
+
+template
+tmp::fluxFieldType>
+EulerDdtScheme::fvcDdtUfCorr
+(
+ const volScalarField& rho,
+ const GeometricField& U,
+ const GeometricField& Uf
+)
+{
+ dimensionedScalar rDeltaT = 1.0/mesh().time().deltaT();
+
+ IOobject ddtIOobject
+ (
+ "ddtCorr(" + rho.name() + ',' + U.name() + ',' + Uf.name() + ')',
mesh().time().timeName(),
mesh()
);
@@ -429,95 +462,144 @@ EulerDdtScheme::fvcDdtPhiCorr
if
(
U.dimensions() == dimVelocity
- && phiAbs.dimensions() == dimVelocity*dimArea
+ && Uf.dimensions() == rho.dimensions()*dimVelocity
)
{
- tmp ddtPhiCorr
+ GeometricField rhoU0
(
- new fluxFieldType
- (
- ddtIOobject,
- rDeltaT
- * this->fvcDdtPhiCoeff(U.oldTime(), phiAbs.oldTime())
- * (
- fvc::interpolate(rA*rho.oldTime())*phiAbs.oldTime()
- - (
- fvc::interpolate(rA*rho.oldTime()*U.oldTime())
- & mesh().Sf()
- )
- )
- )
+ rho.oldTime()*U.oldTime()
);
- ddtPhiCorr().boundaryField() = pTraits::type>::zero;
+ fluxFieldType phiCorr
+ (
+ mesh().Sf() & (Uf.oldTime() - fvc::interpolate(rhoU0))
+ );
- return ddtPhiCorr;
- }
- else if
- (
- U.dimensions() == dimVelocity
- && phiAbs.dimensions() == rho.dimensions()*dimVelocity*dimArea
- )
- {
- tmp ddtPhiCorr
+ return tmp
(
new fluxFieldType
(
ddtIOobject,
- rDeltaT
- * this->fvcDdtPhiCoeff
+ this->fvcDdtPhiCoeff
(
- U.oldTime(),
- phiAbs.oldTime()/fvc::interpolate(rho.oldTime())
- )
- * (
- fvc::interpolate(rA*rho.oldTime())
- * phiAbs.oldTime()/fvc::interpolate(rho.oldTime())
- - (
- fvc::interpolate
- (
- rA*rho.oldTime()*U.oldTime()
- ) & mesh().Sf()
- )
+ rhoU0,
+ mesh().Sf() & Uf.oldTime(),
+ phiCorr
)
+ *rDeltaT*phiCorr
)
);
-
- ddtPhiCorr().boundaryField() = pTraits::type>::zero;
-
- return ddtPhiCorr;
}
else if
(
U.dimensions() == rho.dimensions()*dimVelocity
- && phiAbs.dimensions() == rho.dimensions()*dimVelocity*dimArea
+ && Uf.dimensions() == rho.dimensions()*dimVelocity
)
{
- tmp ddtPhiCorr
+ fluxFieldType phiCorr
+ (
+ mesh().Sf() & (Uf.oldTime() - fvc::interpolate(U.oldTime()))
+ );
+
+ return tmp
(
new fluxFieldType
(
ddtIOobject,
- rDeltaT
- * this->fvcDdtPhiCoeff
- (rho.oldTime(), U.oldTime(), phiAbs.oldTime())
- * (
- fvc::interpolate(rA)*phiAbs.oldTime()
- - (fvc::interpolate(rA*U.oldTime()) & mesh().Sf())
+ this->fvcDdtPhiCoeff
+ (
+ U.oldTime(),
+ mesh().Sf() & Uf.oldTime(),
+ phiCorr
)
+ *rDeltaT*phiCorr
)
);
-
- ddtPhiCorr().boundaryField() = pTraits::type>::zero;
-
- return ddtPhiCorr;
}
else
{
FatalErrorIn
(
"EulerDdtScheme::fvcDdtPhiCorr"
- ) << "dimensions of phiAbs are not correct"
+ ) << "dimensions of Uf are not correct"
+ << abort(FatalError);
+
+ return fluxFieldType::null();
+ }
+}
+
+
+template
+tmp::fluxFieldType>
+EulerDdtScheme::fvcDdtPhiCorr
+(
+ const volScalarField& rho,
+ const GeometricField& U,
+ const fluxFieldType& phi
+)
+{
+ dimensionedScalar rDeltaT = 1.0/mesh().time().deltaT();
+
+ IOobject ddtIOobject
+ (
+ "ddtCorr(" + rho.name() + ',' + U.name() + ',' + phi.name() + ')',
+ mesh().time().timeName(),
+ mesh()
+ );
+
+ if
+ (
+ U.dimensions() == dimVelocity
+ && phi.dimensions() == rho.dimensions()*dimVelocity*dimArea
+ )
+ {
+ GeometricField rhoU0
+ (
+ rho.oldTime()*U.oldTime()
+ );
+
+ fluxFieldType phiCorr
+ (
+ phi.oldTime() - (mesh().Sf() & fvc::interpolate(rhoU0))
+ );
+
+ return tmp
+ (
+ new fluxFieldType
+ (
+ ddtIOobject,
+ this->fvcDdtPhiCoeff(rhoU0, phi.oldTime(), phiCorr)
+ *rDeltaT*phiCorr
+ )
+ );
+ }
+ else if
+ (
+ U.dimensions() == rho.dimensions()*dimVelocity
+ && phi.dimensions() == rho.dimensions()*dimVelocity*dimArea
+ )
+ {
+ fluxFieldType phiCorr
+ (
+ phi.oldTime() - (mesh().Sf() & fvc::interpolate(U.oldTime()))
+ );
+
+ return tmp
+ (
+ new fluxFieldType
+ (
+ ddtIOobject,
+ this->fvcDdtPhiCoeff(U.oldTime(), phi.oldTime(), phiCorr)
+ *rDeltaT*phiCorr
+ )
+ );
+ }
+ else
+ {
+ FatalErrorIn
+ (
+ "EulerDdtScheme::fvcDdtPhiCorr"
+ ) << "dimensions of phi are not correct"
<< abort(FatalError);
return fluxFieldType::null();
diff --git a/src/finiteVolume/finiteVolume/ddtSchemes/EulerDdtScheme/EulerDdtScheme.H b/src/finiteVolume/finiteVolume/ddtSchemes/EulerDdtScheme/EulerDdtScheme.H
index 27085e7dbf..1e6f0fec13 100644
--- a/src/finiteVolume/finiteVolume/ddtSchemes/EulerDdtScheme/EulerDdtScheme.H
+++ b/src/finiteVolume/finiteVolume/ddtSchemes/EulerDdtScheme/EulerDdtScheme.H
@@ -136,16 +136,27 @@ public:
typedef typename ddtScheme::fluxFieldType fluxFieldType;
- tmp fvcDdtPhiCorr
+ tmp fvcDdtUfCorr
(
- const volScalarField& rA,
const GeometricField& U,
- const fluxFieldType& phi
+ const GeometricField& Uf
+ );
+
+ tmp fvcDdtPhiCorr
+ (
+ const GeometricField& U,
+ const fluxFieldType& phi
+ );
+
+ tmp fvcDdtUfCorr
+ (
+ const volScalarField& rho,
+ const GeometricField& U,
+ const GeometricField& Uf
);
tmp fvcDdtPhiCorr
(
- const volScalarField& rA,
const volScalarField& rho,
const GeometricField& U,
const fluxFieldType& phi
@@ -158,10 +169,17 @@ public:
};
+template<>
+tmp EulerDdtScheme::fvcDdtUfCorr
+(
+ const GeometricField& U,
+ const GeometricField& Uf
+);
+
+
template<>
tmp EulerDdtScheme::fvcDdtPhiCorr
(
- const volScalarField& rA,
const volScalarField& U,
const surfaceScalarField& phi
);
@@ -170,7 +188,6 @@ tmp EulerDdtScheme::fvcDdtPhiCorr
template<>
tmp EulerDdtScheme::fvcDdtPhiCorr
(
- const volScalarField& rA,
const volScalarField& rho,
const volScalarField& U,
const surfaceScalarField& phi
diff --git a/src/finiteVolume/finiteVolume/ddtSchemes/SLTSDdtScheme/SLTSDdtScheme.C b/src/finiteVolume/finiteVolume/ddtSchemes/SLTSDdtScheme/SLTSDdtScheme.C
index dd531848e6..95a786c5ac 100644
--- a/src/finiteVolume/finiteVolume/ddtSchemes/SLTSDdtScheme/SLTSDdtScheme.C
+++ b/src/finiteVolume/finiteVolume/ddtSchemes/SLTSDdtScheme/SLTSDdtScheme.C
@@ -481,55 +481,163 @@ SLTSDdtScheme::fvmDdt
}
+template
+tmp::fluxFieldType>
+SLTSDdtScheme::fvcDdtUfCorr
+(
+ const GeometricField& U,
+ const GeometricField& Uf
+)
+{
+ IOobject ddtIOobject
+ (
+ "ddtCorr(" + U.name() + ',' + Uf.name() + ')',
+ mesh().time().timeName(),
+ mesh()
+ );
+
+ const surfaceScalarField rDeltaT(fvc::interpolate(SLrDeltaT()));
+
+ fluxFieldType phiCorr
+ (
+ mesh().Sf() & (Uf.oldTime() - fvc::interpolate(U.oldTime()))
+ );
+
+ return tmp
+ (
+ new fluxFieldType
+ (
+ ddtIOobject,
+ this->fvcDdtPhiCoeff
+ (
+ U.oldTime(),
+ (mesh().Sf() & Uf.oldTime()),
+ phiCorr
+ )
+ *rDeltaT*phiCorr
+ )
+ );
+}
+
+
template
tmp::fluxFieldType>
SLTSDdtScheme::fvcDdtPhiCorr
(
- const volScalarField& rA,
const GeometricField& U,
const fluxFieldType& phi
)
{
IOobject ddtIOobject
(
- "ddtPhiCorr(" + rA.name() + ',' + U.name() + ',' + phi.name() + ')',
+ "ddtCorr(" + U.name() + ',' + phi.name() + ')',
mesh().time().timeName(),
mesh()
);
- if (mesh().moving())
+ const surfaceScalarField rDeltaT(fvc::interpolate(SLrDeltaT()));
+
+ fluxFieldType phiCorr
+ (
+ phi.oldTime() - (mesh().Sf() & fvc::interpolate(U.oldTime()))
+ );
+
+ return tmp
+ (
+ new fluxFieldType
+ (
+ ddtIOobject,
+ this->fvcDdtPhiCoeff(U.oldTime(), phi.oldTime(), phiCorr)
+ *rDeltaT*phiCorr
+ )
+ );
+}
+
+
+template
+tmp::fluxFieldType>
+SLTSDdtScheme::fvcDdtUfCorr
+(
+ const volScalarField& rho,
+ const GeometricField& U,
+ const GeometricField& Uf
+)
+{
+ IOobject ddtIOobject
+ (
+ "ddtCorr(" + rho.name() + ',' + U.name() + ',' + Uf.name() + ')',
+ mesh().time().timeName(),
+ mesh()
+ );
+
+ const surfaceScalarField rDeltaT(fvc::interpolate(SLrDeltaT()));
+
+ if
+ (
+ U.dimensions() == dimVelocity
+ && Uf.dimensions() == dimDensity*dimVelocity
+ )
{
+ GeometricField rhoU0
+ (
+ rho.oldTime()*U.oldTime()
+ );
+
+ fluxFieldType phiCorr
+ (
+ mesh().Sf() & (Uf.oldTime() - fvc::interpolate(rhoU0))
+ );
+
return tmp
(
new fluxFieldType
(
ddtIOobject,
- mesh(),
- dimensioned::type>
+ this->fvcDdtPhiCoeff
(
- "0",
- rA.dimensions()*phi.dimensions()/dimTime,
- pTraits::type>::zero
+ rhoU0,
+ mesh().Sf() & Uf.oldTime(),
+ phiCorr
)
+ *rDeltaT*phiCorr
+ )
+ );
+ }
+ else if
+ (
+ U.dimensions() == dimDensity*dimVelocity
+ && Uf.dimensions() == dimDensity*dimVelocity
+ )
+ {
+ fluxFieldType phiCorr
+ (
+ mesh().Sf() & (Uf.oldTime() - fvc::interpolate(U.oldTime()))
+ );
+
+ return tmp
+ (
+ new fluxFieldType
+ (
+ ddtIOobject,
+ this->fvcDdtPhiCoeff
+ (
+ U.oldTime(),
+ mesh().Sf() & Uf.oldTime(),
+ phiCorr
+ )
+ *rDeltaT*phiCorr
)
);
}
else
{
- const volScalarField rDeltaT(SLrDeltaT());
-
- return tmp
+ FatalErrorIn
(
- new fluxFieldType
- (
- ddtIOobject,
- this->fvcDdtPhiCoeff(U.oldTime(), phi.oldTime())*
- (
- fvc::interpolate(rDeltaT*rA)*phi.oldTime()
- - (fvc::interpolate(rDeltaT*rA*U.oldTime()) & mesh().Sf())
- )
- )
- );
+ "SLTSDdtScheme::fvcDdtPhiCorr"
+ ) << "dimensions of Uf are not correct"
+ << abort(FatalError);
+
+ return fluxFieldType::null();
}
}
@@ -538,122 +646,76 @@ template
tmp::fluxFieldType>
SLTSDdtScheme::fvcDdtPhiCorr
(
- const volScalarField& rA,
const volScalarField& rho,
const GeometricField& U,
const fluxFieldType& phi
)
{
+ dimensionedScalar rDeltaT = 1.0/mesh().time().deltaT();
+
IOobject ddtIOobject
(
- "ddtPhiCorr("
- + rA.name() + ',' + rho.name() + ',' + U.name() + ',' + phi.name() + ')',
+ "ddtCorr(" + rho.name() + ',' + U.name() + ',' + phi.name() + ')',
mesh().time().timeName(),
mesh()
);
- if (mesh().moving())
+ if
+ (
+ U.dimensions() == dimVelocity
+ && phi.dimensions() == rho.dimensions()*dimVelocity*dimArea
+ )
{
+ GeometricField rhoU0
+ (
+ rho.oldTime()*U.oldTime()
+ );
+
+ fluxFieldType phiCorr
+ (
+ phi.oldTime() - (mesh().Sf() & fvc::interpolate(rhoU0))
+ );
+
return tmp
(
new fluxFieldType
(
ddtIOobject,
- mesh(),
- dimensioned::type>
- (
- "0",
- rA.dimensions()*rho.dimensions()*phi.dimensions()/dimTime,
- pTraits::type>::zero
- )
+ this->fvcDdtPhiCoeff(rhoU0, phi.oldTime(), phiCorr)
+ *rDeltaT*phiCorr
+ )
+ );
+ }
+ else if
+ (
+ U.dimensions() == rho.dimensions()*dimVelocity
+ && phi.dimensions() == rho.dimensions()*dimVelocity*dimArea
+ )
+ {
+ fluxFieldType phiCorr
+ (
+ phi.oldTime() - (mesh().Sf() & fvc::interpolate(U.oldTime()))
+ );
+
+ return tmp
+ (
+ new fluxFieldType
+ (
+ ddtIOobject,
+ this->fvcDdtPhiCoeff(U.oldTime(), phi.oldTime(), phiCorr)
+ *rDeltaT*phiCorr
)
);
}
else
{
- const volScalarField rDeltaT(SLrDeltaT());
+ FatalErrorIn
+ (
+ "SLTSDdtScheme::fvcDdtPhiCorr"
+ ) << "dimensions of phi are not correct"
+ << abort(FatalError);
- if
- (
- U.dimensions() == dimVelocity
- && phi.dimensions() == dimVelocity*dimArea
- )
- {
- return tmp
- (
- new fluxFieldType
- (
- ddtIOobject,
- this->fvcDdtPhiCoeff(U.oldTime(), phi.oldTime())
- *(
- fvc::interpolate(rDeltaT*rA*rho.oldTime())*phi.oldTime()
- - (fvc::interpolate(rDeltaT*rA*rho.oldTime()*U.oldTime())
- & mesh().Sf())
- )
- )
- );
- }
- else if
- (
- U.dimensions() == dimVelocity
- && phi.dimensions() == dimDensity*dimVelocity*dimArea
- )
- {
- return tmp
- (
- new fluxFieldType
- (
- ddtIOobject,
- this->fvcDdtPhiCoeff
- (
- U.oldTime(),
- phi.oldTime()/fvc::interpolate(rho.oldTime())
- )
- *(
- fvc::interpolate(rDeltaT*rA*rho.oldTime())
- *phi.oldTime()/fvc::interpolate(rho.oldTime())
- - (
- fvc::interpolate
- (
- rDeltaT*rA*rho.oldTime()*U.oldTime()
- ) & mesh().Sf()
- )
- )
- )
- );
- }
- else if
- (
- U.dimensions() == dimDensity*dimVelocity
- && phi.dimensions() == dimDensity*dimVelocity*dimArea
- )
- {
- return tmp
- (
- new fluxFieldType
- (
- ddtIOobject,
- this->fvcDdtPhiCoeff
- (rho.oldTime(), U.oldTime(), phi.oldTime())
- * (
- fvc::interpolate(rDeltaT*rA)*phi.oldTime()
- - (
- fvc::interpolate(rDeltaT*rA*U.oldTime())&mesh().Sf()
- )
- )
- )
- );
- }
- else
- {
- FatalErrorIn
- (
- "SLTSDdtScheme::fvcDdtPhiCorr"
- ) << "dimensions of phi are not correct"
- << abort(FatalError);
-
- return fluxFieldType::null();
- }
+ return fluxFieldType::null();
}
}
diff --git a/src/finiteVolume/finiteVolume/ddtSchemes/SLTSDdtScheme/SLTSDdtScheme.H b/src/finiteVolume/finiteVolume/ddtSchemes/SLTSDdtScheme/SLTSDdtScheme.H
index 7e1793a75e..dad42a1553 100644
--- a/src/finiteVolume/finiteVolume/ddtSchemes/SLTSDdtScheme/SLTSDdtScheme.H
+++ b/src/finiteVolume/finiteVolume/ddtSchemes/SLTSDdtScheme/SLTSDdtScheme.H
@@ -159,16 +159,27 @@ public:
typedef typename ddtScheme::fluxFieldType fluxFieldType;
- tmp fvcDdtPhiCorr
+ tmp fvcDdtUfCorr
(
- const volScalarField& rA,
const GeometricField& U,
- const fluxFieldType& phi
+ const GeometricField& Uf
+ );
+
+ tmp fvcDdtPhiCorr
+ (
+ const GeometricField& U,
+ const fluxFieldType& phi
+ );
+
+ tmp fvcDdtUfCorr
+ (
+ const volScalarField& rho,
+ const GeometricField& U,
+ const GeometricField& Uf
);
tmp