diff --git a/applications/solvers/DNS/dnsFoam/dnsFoam.C b/applications/solvers/DNS/dnsFoam/dnsFoam.C
index 08073398a7..c29d451091 100644
--- a/applications/solvers/DNS/dnsFoam/dnsFoam.C
+++ b/applications/solvers/DNS/dnsFoam/dnsFoam.C
@@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
- \\ / A nd | Copyright (C) 2011-2012 OpenFOAM Foundation
+ \\ / A nd | Copyright (C) 2011-2013 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@@ -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..11ec1de915 100644
--- a/applications/solvers/compressible/sonicFoam/sonicDyMFoam/sonicDyMFoam.C
+++ b/applications/solvers/compressible/sonicFoam/sonicDyMFoam/sonicDyMFoam.C
@@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
- \\ / A nd | Copyright (C) 2011-2012 OpenFOAM Foundation
+ \\ / A nd | Copyright (C) 2011-2013 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@@ -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..2a0b08074a 100644
--- a/applications/solvers/compressible/sonicFoam/sonicLiquidFoam/sonicLiquidFoam.C
+++ b/applications/solvers/compressible/sonicFoam/sonicLiquidFoam/sonicLiquidFoam.C
@@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
- \\ / A nd | Copyright (C) 2011-2012 OpenFOAM Foundation
+ \\ / A nd | Copyright (C) 2011-2013 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@@ -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..945aab6ba3 100644
--- a/applications/solvers/electromagnetics/mhdFoam/mhdFoam.C
+++ b/applications/solvers/electromagnetics/mhdFoam/mhdFoam.C
@@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
- \\ / A nd | Copyright (C) 2011 OpenFOAM Foundation
+ \\ / A nd | Copyright (C) 2011-2013 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@@ -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..fc5b24ec6a 100644
--- a/applications/solvers/incompressible/pisoFoam/pisoFoam.C
+++ b/applications/solvers/incompressible/pisoFoam/pisoFoam.C
@@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
- \\ / A nd | Copyright (C) 2011-2012 OpenFOAM Foundation
+ \\ / A nd | Copyright (C) 2011-2013 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@@ -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..536d090d04 100644
--- a/applications/solvers/incompressible/shallowWaterFoam/shallowWaterFoam.C
+++ b/applications/solvers/incompressible/shallowWaterFoam/shallowWaterFoam.C
@@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
- \\ / A nd | Copyright (C) 2011-2012 OpenFOAM Foundation
+ \\ / A nd | Copyright (C) 2011-2013 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@@ -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..39b8a069bb 100644
--- a/applications/test/LduMatrix/LduMatrixTest3.C
+++ b/applications/test/LduMatrix/LduMatrixTest3.C
@@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
- \\ / A nd | Copyright (C) 2011-2012 OpenFOAM Foundation
+ \\ / A nd | Copyright (C) 2011-2013 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@@ -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..d817318a8a 100644
--- a/applications/test/PisoFoam/PisoFoam.C
+++ b/applications/test/PisoFoam/PisoFoam.C
@@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
- \\ / A nd | Copyright (C) 2011-2012 OpenFOAM Foundation
+ \\ / A nd | Copyright (C) 2011-2013 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@@ -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/OpenFOAM/db/Time/Time.C b/src/OpenFOAM/db/Time/Time.C
index cb5ca01bff..94c8b9cac4 100644
--- a/src/OpenFOAM/db/Time/Time.C
+++ b/src/OpenFOAM/db/Time/Time.C
@@ -129,6 +129,8 @@ void Foam::Time::adjustDeltaT()
}
}
}
+
+ functionObjects_.adjustTimeStep();
}
@@ -957,17 +959,25 @@ void Foam::Time::setEndTime(const scalar endTime)
}
-void Foam::Time::setDeltaT(const dimensionedScalar& deltaT)
+void Foam::Time::setDeltaT
+(
+ const dimensionedScalar& deltaT,
+ const bool bAdjustDeltaT
+)
{
- setDeltaT(deltaT.value());
+ setDeltaT(deltaT.value(), bAdjustDeltaT);
}
-void Foam::Time::setDeltaT(const scalar deltaT)
+void Foam::Time::setDeltaT(const scalar deltaT, const bool bAdjustDeltaT)
{
deltaT_ = deltaT;
deltaTchanged_ = true;
- adjustDeltaT();
+
+ if (bAdjustDeltaT)
+ {
+ adjustDeltaT();
+ }
}
diff --git a/src/OpenFOAM/db/Time/Time.H b/src/OpenFOAM/db/Time/Time.H
index 695f666c2e..0d7c83818f 100644
--- a/src/OpenFOAM/db/Time/Time.H
+++ b/src/OpenFOAM/db/Time/Time.H
@@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
- \\ / A nd | Copyright (C) 2011-2012 OpenFOAM Foundation
+ \\ / A nd | Copyright (C) 2011-2013 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@@ -527,10 +527,18 @@ public:
virtual void setEndTime(const scalar);
//- Reset time step
- virtual void setDeltaT(const dimensionedScalar&);
+ virtual void setDeltaT
+ (
+ const dimensionedScalar&,
+ const bool adjustDeltaT = true
+ );
//- Reset time step
- virtual void setDeltaT(const scalar);
+ virtual void setDeltaT
+ (
+ const scalar,
+ const bool adjustDeltaT = true
+ );
//- Set time to sub-cycle for the given number of steps
virtual TimeState subCycle(const label nSubCycles);
diff --git a/src/OpenFOAM/db/functionObjects/OutputFilterFunctionObject/OutputFilterFunctionObject.C b/src/OpenFOAM/db/functionObjects/OutputFilterFunctionObject/OutputFilterFunctionObject.C
index 7c6b95ec98..019a4f54aa 100644
--- a/src/OpenFOAM/db/functionObjects/OutputFilterFunctionObject/OutputFilterFunctionObject.C
+++ b/src/OpenFOAM/db/functionObjects/OutputFilterFunctionObject/OutputFilterFunctionObject.C
@@ -40,6 +40,7 @@ void Foam::OutputFilterFunctionObject::readDict()
dict_.readIfPresent("storeFilter", storeFilter_);
dict_.readIfPresent("timeStart", timeStart_);
dict_.readIfPresent("timeEnd", timeEnd_);
+ dict_.readIfPresent("nStepsToStartTimeChange", nStepsToStartTimeChange_);
}
@@ -214,6 +215,52 @@ bool Foam::OutputFilterFunctionObject::timeSet()
}
+template
+bool Foam::OutputFilterFunctionObject::adjustTimeStep()
+{
+ if
+ (
+ active()
+ && outputControl_.outputControl()
+ == outputFilterOutputControl::ocAdjustableTime
+ )
+ {
+ const label outputTimeIndex = outputControl_.outputTimeLastDump();
+ const scalar writeInterval = outputControl_.writeInterval();
+
+ scalar timeToNextWrite = max
+ (
+ 0.0,
+ (outputTimeIndex + 1)*writeInterval
+ - (time_.value() - time_.startTime().value())
+ );
+
+ scalar deltaT = time_.deltaTValue();
+
+ scalar nSteps = timeToNextWrite/deltaT - SMALL;
+
+ // function objects modify deltaT inside nStepsToStartTimeChange range
+ // NOTE: Potential problem if two function objects dump inside the same
+ //interval
+ if (nSteps < nStepsToStartTimeChange_)
+ {
+ label nStepsToNextWrite = label(nSteps) + 1;
+
+ scalar newDeltaT = timeToNextWrite/nStepsToNextWrite;
+
+ //Adjust time step
+ if (newDeltaT < deltaT)
+ {
+ deltaT = max(newDeltaT, 0.2*deltaT);
+ const_cast