Merge branch 'master' into dicts

This commit is contained in:
andy
2010-10-12 12:13:39 +01:00
45 changed files with 1551 additions and 1517 deletions

View File

@ -85,15 +85,15 @@ int main(int argc, char *argv[])
for (int corr=1; corr<=1; corr++) for (int corr=1; corr<=1; corr++)
{ {
volScalarField rUA = 1.0/UEqn.A(); volScalarField rAU = 1.0/UEqn.A();
U = rUA*UEqn.H(); U = rAU*UEqn.H();
phi = (fvc::interpolate(U) & mesh.Sf()) phi = (fvc::interpolate(U) & mesh.Sf())
+ fvc::ddtPhiCorr(rUA, U, phi); + fvc::ddtPhiCorr(rAU, U, phi);
fvScalarMatrix pEqn fvScalarMatrix pEqn
( (
fvm::laplacian(rUA, p) == fvc::div(phi) fvm::laplacian(rAU, p) == fvc::div(phi)
); );
pEqn.solve(); pEqn.solve();
@ -102,7 +102,7 @@ int main(int argc, char *argv[])
#include "continuityErrs.H" #include "continuityErrs.H"
U -= rUA*fvc::grad(p); U -= rAU*fvc::grad(p);
U.correctBoundaryConditions(); U.correctBoundaryConditions();
} }

View File

@ -1,6 +1,6 @@
rho = thermo.rho(); rho = thermo.rho();
volScalarField rUA = 1.0/UEqn.A(); volScalarField rAU = 1.0/UEqn.A();
U = invA & UEqn.H(); U = invA & UEqn.H();
if (transonic) if (transonic)
@ -11,7 +11,7 @@ if (transonic)
fvc::interpolate(psi) fvc::interpolate(psi)
*( *(
(fvc::interpolate(U) & mesh.Sf()) (fvc::interpolate(U) & mesh.Sf())
+ fvc::ddtPhiCorr(rUA, rho, U, phi) + fvc::ddtPhiCorr(rAU, rho, U, phi)
) )
); );
@ -38,7 +38,7 @@ else
fvc::interpolate(rho)* fvc::interpolate(rho)*
( (
(fvc::interpolate(U) & mesh.Sf()) (fvc::interpolate(U) & mesh.Sf())
+ fvc::ddtPhiCorr(rUA, rho, U, phi) + fvc::ddtPhiCorr(rAU, rho, U, phi)
); );
for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++) for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)

View File

@ -1,7 +1,7 @@
rho = thermo.rho(); rho = thermo.rho();
volScalarField rUA = 1.0/UEqn.A(); volScalarField rAU = 1.0/UEqn.A();
U = rUA*UEqn.H(); U = rAU*UEqn.H();
if (transonic) if (transonic)
{ {
@ -11,7 +11,7 @@ if (transonic)
fvc::interpolate(psi) fvc::interpolate(psi)
*( *(
(fvc::interpolate(U) & mesh.Sf()) (fvc::interpolate(U) & mesh.Sf())
+ fvc::ddtPhiCorr(rUA, rho, U, phi) + fvc::ddtPhiCorr(rAU, rho, U, phi)
) )
); );
@ -21,7 +21,7 @@ if (transonic)
( (
fvm::ddt(psi, p) fvm::ddt(psi, p)
+ fvm::div(phid, p) + fvm::div(phid, p)
- fvm::laplacian(rho*rUA, p) - fvm::laplacian(rho*rAU, p)
); );
pEqn.solve(); pEqn.solve();
@ -38,7 +38,7 @@ else
fvc::interpolate(rho) fvc::interpolate(rho)
*( *(
(fvc::interpolate(U) & mesh.Sf()) (fvc::interpolate(U) & mesh.Sf())
+ fvc::ddtPhiCorr(rUA, rho, U, phi) + fvc::ddtPhiCorr(rAU, rho, U, phi)
); );
for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++) for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
@ -47,7 +47,7 @@ else
( (
fvm::ddt(psi, p) fvm::ddt(psi, p)
+ fvc::div(phi) + fvc::div(phi)
- fvm::laplacian(rho*rUA, p) - fvm::laplacian(rho*rAU, p)
); );
pEqn.solve(); pEqn.solve();
@ -62,7 +62,7 @@ else
#include "rhoEqn.H" #include "rhoEqn.H"
#include "compressibleContinuityErrs.H" #include "compressibleContinuityErrs.H"
U -= rUA*fvc::grad(p); U -= rAU*fvc::grad(p);
U.correctBoundaryConditions(); U.correctBoundaryConditions();
DpDt = fvc::DDt(surfaceScalarField("phiU", phi/fvc::interpolate(rho)), p); DpDt = fvc::DDt(surfaceScalarField("phiU", phi/fvc::interpolate(rho)), p);

View File

@ -1,7 +1,7 @@
rho = thermo.rho(); rho = thermo.rho();
volScalarField rUA = 1.0/UEqn.A(); volScalarField rAU = 1.0/UEqn.A();
U = rUA*UEqn.H(); U = rAU*UEqn.H();
if (transonic) if (transonic)
{ {
@ -11,7 +11,7 @@ if (transonic)
fvc::interpolate(psi) fvc::interpolate(psi)
*( *(
(fvc::interpolate(U) & mesh.Sf()) (fvc::interpolate(U) & mesh.Sf())
+ fvc::ddtPhiCorr(rUA, rho, U, phi) + fvc::ddtPhiCorr(rAU, rho, U, phi)
) )
); );
@ -21,7 +21,7 @@ if (transonic)
( (
fvm::ddt(psi, p) fvm::ddt(psi, p)
+ fvm::div(phid, p) + fvm::div(phid, p)
- fvm::laplacian(rho*rUA, p) - fvm::laplacian(rho*rAU, p)
== ==
Sevap Sevap
); );
@ -40,7 +40,7 @@ else
fvc::interpolate(rho) fvc::interpolate(rho)
*( *(
(fvc::interpolate(U) & mesh.Sf()) (fvc::interpolate(U) & mesh.Sf())
+ fvc::ddtPhiCorr(rUA, rho, U, phi) + fvc::ddtPhiCorr(rAU, rho, U, phi)
); );
for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++) for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
@ -49,7 +49,7 @@ else
( (
fvm::ddt(psi, p) fvm::ddt(psi, p)
+ fvc::div(phi) + fvc::div(phi)
- fvm::laplacian(rho*rUA, p) - fvm::laplacian(rho*rAU, p)
== ==
Sevap Sevap
); );
@ -66,7 +66,7 @@ else
#include "rhoEqn.H" #include "rhoEqn.H"
#include "compressibleContinuityErrs.H" #include "compressibleContinuityErrs.H"
U -= rUA*fvc::grad(p); U -= rAU*fvc::grad(p);
U.correctBoundaryConditions(); U.correctBoundaryConditions();
DpDt = fvc::DDt(surfaceScalarField("phiU", phi/fvc::interpolate(rho)), p); DpDt = fvc::DDt(surfaceScalarField("phiU", phi/fvc::interpolate(rho)), p);

View File

@ -1,7 +1,7 @@
rho = thermo.rho(); rho = thermo.rho();
volScalarField rUA = 1.0/UEqn.A(); volScalarField rAU = 1.0/UEqn.A();
U = rUA*UEqn.H(); U = rAU*UEqn.H();
if (transonic) if (transonic)
{ {
@ -18,7 +18,7 @@ if (transonic)
( (
fvm::ddt(psi, p) fvm::ddt(psi, p)
+ fvm::div(phid, p, "div(phid,p)") + fvm::div(phid, p, "div(phid,p)")
- fvm::laplacian(rho*rUA, p) - fvm::laplacian(rho*rAU, p)
); );
pEqn.solve(); pEqn.solve();
@ -40,7 +40,7 @@ else
( (
fvm::ddt(psi, p) fvm::ddt(psi, p)
+ fvc::div(phi) + fvc::div(phi)
- fvm::laplacian(rho*rUA, p) - fvm::laplacian(rho*rAU, p)
); );
pEqn.solve(); pEqn.solve();
@ -55,7 +55,7 @@ else
#include "rhoEqn.H" #include "rhoEqn.H"
#include "compressibleContinuityErrs.H" #include "compressibleContinuityErrs.H"
U -= rUA*fvc::grad(p); U -= rAU*fvc::grad(p);
U.correctBoundaryConditions(); U.correctBoundaryConditions();
DpDt = fvc::DDt(surfaceScalarField("phiU", phi/fvc::interpolate(rho)), p); DpDt = fvc::DDt(surfaceScalarField("phiU", phi/fvc::interpolate(rho)), p);

View File

@ -1,29 +1,29 @@
rho = thermo.rho(); rho = thermo.rho();
volScalarField rUA = 1.0/UEqn.A(); volScalarField rAU = 1.0/UEqn.A();
surfaceScalarField rhorUAf("(rho*(1|A(U)))", fvc::interpolate(rho*rUA)); surfaceScalarField rhorAUf("(rho*(1|A(U)))", fvc::interpolate(rho*rAU));
U = rUA*UEqn.H(); U = rAU*UEqn.H();
surfaceScalarField phiU surfaceScalarField phiU
( (
fvc::interpolate(rho) fvc::interpolate(rho)
*( *(
(fvc::interpolate(U) & mesh.Sf()) (fvc::interpolate(U) & mesh.Sf())
+ fvc::ddtPhiCorr(rUA, rho, U, phi) + fvc::ddtPhiCorr(rAU, rho, U, phi)
) )
); );
phi = phiU - rhorUAf*ghf*fvc::snGrad(rho)*mesh.magSf(); phi = phiU - rhorAUf*ghf*fvc::snGrad(rho)*mesh.magSf();
for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++) for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
{ {
surfaceScalarField rhorUAf = fvc::interpolate(rho*rUA); surfaceScalarField rhorAUf = fvc::interpolate(rho*rAU);
fvScalarMatrix p_rghEqn fvScalarMatrix p_rghEqn
( (
fvm::ddt(psi, p_rgh) + fvc::ddt(psi, rho)*gh fvm::ddt(psi, p_rgh) + fvc::ddt(psi, rho)*gh
+ fvc::div(phi) + fvc::div(phi)
- fvm::laplacian(rhorUAf, p_rgh) - fvm::laplacian(rhorAUf, p_rgh)
); );
p_rghEqn.solve p_rghEqn.solve
@ -52,7 +52,7 @@ p = p_rgh + rho*gh;
#include "rhoEqn.H" #include "rhoEqn.H"
#include "compressibleContinuityErrs.H" #include "compressibleContinuityErrs.H"
U += rUA*fvc::reconstruct((phi - phiU)/rhorUAf); U += rAU*fvc::reconstruct((phi - phiU)/rhorAUf);
U.correctBoundaryConditions(); U.correctBoundaryConditions();
DpDt = fvc::DDt(surfaceScalarField("phiU", phi/fvc::interpolate(rho)), p); DpDt = fvc::DDt(surfaceScalarField("phiU", phi/fvc::interpolate(rho)), p);

View File

@ -1,7 +1,7 @@
rho = thermo.rho(); rho = thermo.rho();
volScalarField rUA = 1.0/UEqn.A(); volScalarField rAU = 1.0/UEqn.A();
U = rUA*UEqn.H(); U = rAU*UEqn.H();
if (transonic) if (transonic)
{ {
@ -11,7 +11,7 @@ if (transonic)
fvc::interpolate(psi) fvc::interpolate(psi)
*( *(
(fvc::interpolate(U) & mesh.Sf()) (fvc::interpolate(U) & mesh.Sf())
+ fvc::ddtPhiCorr(rUA, rho, U, phi) + fvc::ddtPhiCorr(rAU, rho, U, phi)
) )
); );
@ -21,7 +21,7 @@ if (transonic)
( (
fvm::ddt(psi, p) fvm::ddt(psi, p)
+ fvm::div(phid, p) + fvm::div(phid, p)
- fvm::laplacian(rho*rUA, p) - fvm::laplacian(rho*rAU, p)
); );
pEqn.solve(); pEqn.solve();
@ -38,7 +38,7 @@ else
fvc::interpolate(rho) fvc::interpolate(rho)
*( *(
(fvc::interpolate(U) & mesh.Sf()) (fvc::interpolate(U) & mesh.Sf())
+ fvc::ddtPhiCorr(rUA, rho, U, phi) + fvc::ddtPhiCorr(rAU, rho, U, phi)
); );
for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++) for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
@ -47,7 +47,7 @@ else
( (
fvm::ddt(psi, p) fvm::ddt(psi, p)
+ fvc::div(phi) + fvc::div(phi)
- fvm::laplacian(rho*rUA, p) - fvm::laplacian(rho*rAU, p)
); );
pEqn.solve(); pEqn.solve();
@ -62,7 +62,7 @@ else
#include "rhoEqn.H" #include "rhoEqn.H"
#include "compressibleContinuityErrs.H" #include "compressibleContinuityErrs.H"
U -= rUA*fvc::grad(p); U -= rAU*fvc::grad(p);
U.correctBoundaryConditions(); U.correctBoundaryConditions();
DpDt = fvc::DDt(surfaceScalarField("phiU", phi/fvc::interpolate(rho)), p); DpDt = fvc::DDt(surfaceScalarField("phiU", phi/fvc::interpolate(rho)), p);

View File

@ -5,14 +5,14 @@
// pressure solution - done in 2 parts. Part 1: // pressure solution - done in 2 parts. Part 1:
thermo.rho() -= psi*p; thermo.rho() -= psi*p;
volScalarField rUA = 1.0/UEqn.A(); volScalarField rAU = 1.0/UEqn.A();
U = rUA*UEqn.H(); U = rAU*UEqn.H();
if (transonic) if (transonic)
{ {
surfaceScalarField phiv = surfaceScalarField phiv =
(fvc::interpolate(U) & mesh.Sf()) (fvc::interpolate(U) & mesh.Sf())
+ fvc::ddtPhiCorr(rUA, rho, U, phi); + fvc::ddtPhiCorr(rAU, rho, U, phi);
phi = fvc::interpolate(rho)*phiv; phi = fvc::interpolate(rho)*phiv;
@ -28,7 +28,7 @@
( (
fvc::ddt(rho) + fvc::div(phi) fvc::ddt(rho) + fvc::div(phi)
+ correction(fvm::ddt(psi, p) + fvm::div(phid, p)) + correction(fvm::ddt(psi, p) + fvm::div(phid, p))
- fvm::laplacian(rho*rUA, p) - fvm::laplacian(rho*rAU, p)
); );
pEqn.solve pEqn.solve
@ -58,7 +58,7 @@
fvc::interpolate(rho) fvc::interpolate(rho)
*( *(
(fvc::interpolate(U) & mesh.Sf()) (fvc::interpolate(U) & mesh.Sf())
+ fvc::ddtPhiCorr(rUA, rho, U, phi) + fvc::ddtPhiCorr(rAU, rho, U, phi)
); );
for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++) for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
@ -67,7 +67,7 @@
( (
fvc::ddt(rho) + psi*correction(fvm::ddt(p)) fvc::ddt(rho) + psi*correction(fvm::ddt(p))
+ fvc::div(phi) + fvc::div(phi)
- fvm::laplacian(rho*rUA, p) - fvm::laplacian(rho*rAU, p)
); );
pEqn.solve pEqn.solve
@ -98,7 +98,7 @@
#include "rhoEqn.H" #include "rhoEqn.H"
#include "compressibleContinuityErrs.H" #include "compressibleContinuityErrs.H"
U -= rUA*fvc::grad(p); U -= rAU*fvc::grad(p);
U.correctBoundaryConditions(); U.correctBoundaryConditions();
DpDt = fvc::DDt(surfaceScalarField("phiU", phi/fvc::interpolate(rho)), p); DpDt = fvc::DDt(surfaceScalarField("phiU", phi/fvc::interpolate(rho)), p);

View File

@ -9,7 +9,7 @@ tmp<fvVectorMatrix> UEqn
UEqn().relax(); UEqn().relax();
volScalarField rUA = 1.0/UEqn().A(); volScalarField rAU = 1.0/UEqn().A();
if (momentumPredictor) if (momentumPredictor)
{ {
@ -17,6 +17,6 @@ if (momentumPredictor)
} }
else else
{ {
U = rUA*(UEqn().H() - fvc::grad(p)); U = rAU*(UEqn().H() - fvc::grad(p));
U.correctBoundaryConditions(); U.correctBoundaryConditions();
} }

View File

@ -1,6 +1,6 @@
rho = thermo.rho(); rho = thermo.rho();
U = rUA*UEqn().H(); U = rAU*UEqn().H();
if (nCorr <= 1) if (nCorr <= 1)
{ {
@ -15,7 +15,7 @@ if (transonic)
fvc::interpolate(psi) fvc::interpolate(psi)
*( *(
(fvc::interpolate(U) & mesh.Sf()) (fvc::interpolate(U) & mesh.Sf())
+ fvc::ddtPhiCorr(rUA, rho, U, phi) + fvc::ddtPhiCorr(rAU, rho, U, phi)
) )
); );
@ -25,7 +25,7 @@ if (transonic)
( (
fvm::ddt(psi, p) fvm::ddt(psi, p)
+ fvm::div(phid, p) + fvm::div(phid, p)
- fvm::laplacian(rho*rUA, p) - fvm::laplacian(rho*rAU, p)
); );
pEqn.solve pEqn.solve
@ -55,7 +55,7 @@ else
fvc::interpolate(rho)* fvc::interpolate(rho)*
( (
(fvc::interpolate(U) & mesh.Sf()) (fvc::interpolate(U) & mesh.Sf())
+ fvc::ddtPhiCorr(rUA, rho, U, phi) + fvc::ddtPhiCorr(rAU, rho, U, phi)
); );
for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++) for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
@ -65,7 +65,7 @@ else
( (
fvm::ddt(psi, p) fvm::ddt(psi, p)
+ fvc::div(phi) + fvc::div(phi)
- fvm::laplacian(rho*rUA, p) - fvm::laplacian(rho*rAU, p)
); );
pEqn.solve pEqn.solve
@ -101,7 +101,7 @@ rho = thermo.rho();
Info<< "rho max/min : " << max(rho).value() Info<< "rho max/min : " << max(rho).value()
<< " " << min(rho).value() << endl; << " " << min(rho).value() << endl;
U -= rUA*fvc::grad(p); U -= rAU*fvc::grad(p);
U.correctBoundaryConditions(); U.correctBoundaryConditions();
DpDt = fvc::DDt(surfaceScalarField("phiU", phi/fvc::interpolate(rho)), p); DpDt = fvc::DDt(surfaceScalarField("phiU", phi/fvc::interpolate(rho)), p);

View File

@ -12,7 +12,7 @@ UEqn().relax();
mrfZones.addCoriolis(rho, UEqn()); mrfZones.addCoriolis(rho, UEqn());
pZones.addResistance(UEqn()); pZones.addResistance(UEqn());
volScalarField rUA = 1.0/UEqn().A(); volScalarField rAU = 1.0/UEqn().A();
if (momentumPredictor) if (momentumPredictor)
{ {
@ -20,6 +20,6 @@ if (momentumPredictor)
} }
else else
{ {
U = rUA*(UEqn().H() - fvc::grad(p)); U = rAU*(UEqn().H() - fvc::grad(p));
U.correctBoundaryConditions(); U.correctBoundaryConditions();
} }

View File

@ -1,7 +1,7 @@
rho = thermo.rho(); rho = thermo.rho();
volScalarField rUA = 1.0/UEqn().A(); volScalarField rAU = 1.0/UEqn().A();
U = rUA*UEqn().H(); U = rAU*UEqn().H();
if (nCorr <= 1) if (nCorr <= 1)
{ {
@ -16,7 +16,7 @@ if (transonic)
fvc::interpolate(psi) fvc::interpolate(psi)
*( *(
(fvc::interpolate(U) & mesh.Sf()) (fvc::interpolate(U) & mesh.Sf())
+ fvc::ddtPhiCorr(rUA, rho, U, phi) + fvc::ddtPhiCorr(rAU, rho, U, phi)
) )
); );
mrfZones.relativeFlux(fvc::interpolate(psi), phid); mrfZones.relativeFlux(fvc::interpolate(psi), phid);
@ -27,7 +27,7 @@ if (transonic)
( (
fvm::ddt(psi, p) fvm::ddt(psi, p)
+ fvm::div(phid, p) + fvm::div(phid, p)
- fvm::laplacian(rho*rUA, p) - fvm::laplacian(rho*rAU, p)
); );
pEqn.solve pEqn.solve
@ -57,7 +57,7 @@ else
fvc::interpolate(rho)* fvc::interpolate(rho)*
( (
(fvc::interpolate(U) & mesh.Sf()) (fvc::interpolate(U) & mesh.Sf())
//+ fvc::ddtPhiCorr(rUA, rho, U, phi) //+ fvc::ddtPhiCorr(rAU, rho, U, phi)
); );
mrfZones.relativeFlux(fvc::interpolate(rho), phi); mrfZones.relativeFlux(fvc::interpolate(rho), phi);
@ -68,7 +68,7 @@ else
( (
fvm::ddt(psi, p) fvm::ddt(psi, p)
+ fvc::div(phi) + fvc::div(phi)
- fvm::laplacian(rho*rUA, p) - fvm::laplacian(rho*rAU, p)
); );
pEqn.solve pEqn.solve
@ -109,7 +109,7 @@ else
<< " " << min(rho).value() << endl; << " " << min(rho).value() << endl;
} }
U -= rUA*fvc::grad(p); U -= rAU*fvc::grad(p);
U.correctBoundaryConditions(); U.correctBoundaryConditions();
DpDt = fvc::DDt(surfaceScalarField("phiU", phi/fvc::interpolate(rho)), p); DpDt = fvc::DDt(surfaceScalarField("phiU", phi/fvc::interpolate(rho)), p);

View File

@ -3,8 +3,8 @@ rho = max(rho, rhoMin);
rho = min(rho, rhoMax); rho = min(rho, rhoMax);
rho.relax(); rho.relax();
volScalarField rUA = 1.0/UEqn().A(); volScalarField rAU = 1.0/UEqn().A();
U = rUA*UEqn().H(); U = rAU*UEqn().H();
UEqn.clear(); UEqn.clear();
bool closedVolume = false; bool closedVolume = false;
@ -22,7 +22,7 @@ if (transonic)
fvScalarMatrix pEqn fvScalarMatrix pEqn
( (
fvm::div(phid, p) fvm::div(phid, p)
- fvm::laplacian(rho*rUA, p) - fvm::laplacian(rho*rAU, p)
); );
// Relax the pressure equation to ensure diagonal-dominance // Relax the pressure equation to ensure diagonal-dominance
@ -47,7 +47,7 @@ else
{ {
fvScalarMatrix pEqn fvScalarMatrix pEqn
( (
fvm::laplacian(rho*rUA, p) == fvc::div(phi) fvm::laplacian(rho*rAU, p) == fvc::div(phi)
); );
pEqn.setReference(pRefCell, pRefValue); pEqn.setReference(pRefCell, pRefValue);
@ -67,7 +67,7 @@ else
// Explicitly relax pressure for momentum corrector // Explicitly relax pressure for momentum corrector
p.relax(); p.relax();
U -= rUA*fvc::grad(p); U -= rAU*fvc::grad(p);
U.correctBoundaryConditions(); U.correctBoundaryConditions();
// For closed-volume cases adjust the pressure and density levels // For closed-volume cases adjust the pressure and density levels

View File

@ -1,7 +1,7 @@
rho = thermo.rho(); rho = thermo.rho();
volScalarField rUA = 1.0/UEqn.A(); volScalarField rAU = 1.0/UEqn.A();
U = rUA*UEqn.H(); U = rAU*UEqn.H();
surfaceScalarField phid surfaceScalarField phid
( (
@ -9,7 +9,7 @@ surfaceScalarField phid
fvc::interpolate(psi) fvc::interpolate(psi)
*( *(
(fvc::interpolate(U) & mesh.Sf()) (fvc::interpolate(U) & mesh.Sf())
+ fvc::ddtPhiCorr(rUA, rho, U, phi) + fvc::ddtPhiCorr(rAU, rho, U, phi)
) )
); );
@ -19,7 +19,7 @@ for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
( (
fvm::ddt(psi, p) fvm::ddt(psi, p)
+ fvm::div(phid, p) + fvm::div(phid, p)
- fvm::laplacian(rho*rUA, p) - fvm::laplacian(rho*rAU, p)
); );
pEqn.solve(); pEqn.solve();
@ -33,5 +33,5 @@ for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
#include "rhoEqn.H" #include "rhoEqn.H"
#include "compressibleContinuityErrs.H" #include "compressibleContinuityErrs.H"
U -= rUA*fvc::grad(p); U -= rAU*fvc::grad(p);
U.correctBoundaryConditions(); U.correctBoundaryConditions();

View File

@ -71,8 +71,8 @@ int main(int argc, char *argv[])
for (int corr=0; corr<nCorr; corr++) for (int corr=0; corr<nCorr; corr++)
{ {
volScalarField rUA = 1.0/UEqn.A(); volScalarField rAU = 1.0/UEqn.A();
U = rUA*UEqn.H(); U = rAU*UEqn.H();
surfaceScalarField phid surfaceScalarField phid
( (
@ -80,7 +80,7 @@ int main(int argc, char *argv[])
psi psi
*( *(
(fvc::interpolate(U) & mesh.Sf()) (fvc::interpolate(U) & mesh.Sf())
+ fvc::ddtPhiCorr(rUA, rho, U, phi) + fvc::ddtPhiCorr(rAU, rho, U, phi)
) )
); );
@ -91,7 +91,7 @@ int main(int argc, char *argv[])
fvm::ddt(psi, p) fvm::ddt(psi, p)
+ fvc::div(phi) + fvc::div(phi)
+ fvm::div(phid, p) + fvm::div(phid, p)
- fvm::laplacian(rho*rUA, p) - fvm::laplacian(rho*rAU, p)
); );
pEqn.solve(); pEqn.solve();
@ -100,7 +100,7 @@ int main(int argc, char *argv[])
#include "compressibleContinuityErrs.H" #include "compressibleContinuityErrs.H"
U -= rUA*fvc::grad(p); U -= rAU*fvc::grad(p);
U.correctBoundaryConditions(); U.correctBoundaryConditions();
} }

View File

@ -92,18 +92,18 @@ int main(int argc, char *argv[])
for (int corr=0; corr<nCorr; corr++) for (int corr=0; corr<nCorr; corr++)
{ {
volScalarField rUA = 1.0/UEqn.A(); volScalarField rAU = 1.0/UEqn.A();
U = rUA*UEqn.H(); U = rAU*UEqn.H();
phi = (fvc::interpolate(U) & mesh.Sf()) phi = (fvc::interpolate(U) & mesh.Sf())
+ fvc::ddtPhiCorr(rUA, U, phi); + fvc::ddtPhiCorr(rAU, U, phi);
for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++) for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
{ {
fvScalarMatrix pEqn fvScalarMatrix pEqn
( (
fvm::laplacian(rUA, p) == fvc::div(phi) fvm::laplacian(rAU, p) == fvc::div(phi)
); );
pEqn.setReference(pRefCell, pRefValue); pEqn.setReference(pRefCell, pRefValue);
@ -117,7 +117,7 @@ int main(int argc, char *argv[])
#include "continuityErrs.H" #include "continuityErrs.H"
U -= rUA*fvc::grad(p); U -= rAU*fvc::grad(p);
U.correctBoundaryConditions(); U.correctBoundaryConditions();
} }
} }

View File

@ -1,20 +1,20 @@
{ {
volScalarField rUA("rUA", 1.0/UEqn.A()); volScalarField rAU("rAU", 1.0/UEqn.A());
surfaceScalarField rUAf("(1|A(U))", fvc::interpolate(rUA)); surfaceScalarField rAUf("(1|A(U))", fvc::interpolate(rAU));
U = rUA*UEqn.H(); U = rAU*UEqn.H();
phi = (fvc::interpolate(U) & mesh.Sf()) phi = (fvc::interpolate(U) & mesh.Sf())
+ fvc::ddtPhiCorr(rUA, U, phi); + fvc::ddtPhiCorr(rAU, U, phi);
surfaceScalarField buoyancyPhi = rUAf*ghf*fvc::snGrad(rhok)*mesh.magSf(); surfaceScalarField buoyancyPhi = rAUf*ghf*fvc::snGrad(rhok)*mesh.magSf();
phi -= buoyancyPhi; phi -= buoyancyPhi;
for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++) for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
{ {
fvScalarMatrix p_rghEqn fvScalarMatrix p_rghEqn
( (
fvm::laplacian(rUAf, p_rgh) == fvc::div(phi) fvm::laplacian(rAUf, p_rgh) == fvc::div(phi)
); );
p_rghEqn.setReference(pRefCell, getRefCellValue(p_rgh, pRefCell)); p_rghEqn.setReference(pRefCell, getRefCellValue(p_rgh, pRefCell));
@ -44,7 +44,7 @@
// Correct the momentum source with the pressure gradient flux // Correct the momentum source with the pressure gradient flux
// calculated from the relaxed pressure // calculated from the relaxed pressure
U -= rUA*fvc::reconstruct((buoyancyPhi + p_rghEqn.flux())/rUAf); U -= rAU*fvc::reconstruct((buoyancyPhi + p_rghEqn.flux())/rAUf);
U.correctBoundaryConditions(); U.correctBoundaryConditions();
} }
} }

View File

@ -1,21 +1,21 @@
{ {
volScalarField rUA("rUA", 1.0/UEqn().A()); volScalarField rAU("rAU", 1.0/UEqn().A());
surfaceScalarField rUAf("(1|A(U))", fvc::interpolate(rUA)); surfaceScalarField rAUf("(1|A(U))", fvc::interpolate(rAU));
U = rUA*UEqn().H(); U = rAU*UEqn().H();
UEqn.clear(); UEqn.clear();
phi = fvc::interpolate(U) & mesh.Sf(); phi = fvc::interpolate(U) & mesh.Sf();
adjustPhi(phi, U, p_rgh); adjustPhi(phi, U, p_rgh);
surfaceScalarField buoyancyPhi = rUAf*ghf*fvc::snGrad(rhok)*mesh.magSf(); surfaceScalarField buoyancyPhi = rAUf*ghf*fvc::snGrad(rhok)*mesh.magSf();
phi -= buoyancyPhi; phi -= buoyancyPhi;
for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++) for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
{ {
fvScalarMatrix p_rghEqn fvScalarMatrix p_rghEqn
( (
fvm::laplacian(rUAf, p_rgh) == fvc::div(phi) fvm::laplacian(rAUf, p_rgh) == fvc::div(phi)
); );
p_rghEqn.setReference(pRefCell, getRefCellValue(p_rgh, pRefCell)); p_rghEqn.setReference(pRefCell, getRefCellValue(p_rgh, pRefCell));
@ -32,7 +32,7 @@
// Correct the momentum source with the pressure gradient flux // Correct the momentum source with the pressure gradient flux
// calculated from the relaxed pressure // calculated from the relaxed pressure
U -= rUA*fvc::reconstruct((buoyancyPhi + p_rghEqn.flux())/rUAf); U -= rAU*fvc::reconstruct((buoyancyPhi + p_rghEqn.flux())/rAUf);
U.correctBoundaryConditions(); U.correctBoundaryConditions();
} }
} }

View File

@ -5,18 +5,18 @@
// pressure solution - done in 2 parts. Part 1: // pressure solution - done in 2 parts. Part 1:
thermo.rho() -= psi*p_rgh; thermo.rho() -= psi*p_rgh;
volScalarField rUA = 1.0/UEqn.A(); volScalarField rAU = 1.0/UEqn.A();
surfaceScalarField rhorUAf("(rho*(1|A(U)))", fvc::interpolate(rho*rUA)); surfaceScalarField rhorAUf("(rho*(1|A(U)))", fvc::interpolate(rho*rAU));
U = rUA*UEqn.H(); U = rAU*UEqn.H();
phi = fvc::interpolate(rho)* phi = fvc::interpolate(rho)*
( (
(fvc::interpolate(U) & mesh.Sf()) (fvc::interpolate(U) & mesh.Sf())
+ fvc::ddtPhiCorr(rUA, rho, U, phi) + fvc::ddtPhiCorr(rAU, rho, U, phi)
); );
surfaceScalarField buoyancyPhi = -rhorUAf*ghf*fvc::snGrad(rho)*mesh.magSf(); surfaceScalarField buoyancyPhi = -rhorAUf*ghf*fvc::snGrad(rho)*mesh.magSf();
phi += buoyancyPhi; phi += buoyancyPhi;
for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++) for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
@ -25,7 +25,7 @@
( (
fvc::ddt(rho) + psi*correction(fvm::ddt(p_rgh)) fvc::ddt(rho) + psi*correction(fvm::ddt(p_rgh))
+ fvc::div(phi) + fvc::div(phi)
- fvm::laplacian(rhorUAf, p_rgh) - fvm::laplacian(rhorAUf, p_rgh)
); );
p_rghEqn.solve p_rghEqn.solve
@ -53,7 +53,7 @@
// Correct the momentum source with the pressure gradient flux // Correct the momentum source with the pressure gradient flux
// calculated from the relaxed pressure // calculated from the relaxed pressure
U += rUA*fvc::reconstruct((buoyancyPhi + p_rghEqn.flux())/rhorUAf); U += rAU*fvc::reconstruct((buoyancyPhi + p_rghEqn.flux())/rhorAUf);
U.correctBoundaryConditions(); U.correctBoundaryConditions();
} }
} }

View File

@ -2,23 +2,23 @@
rho = thermo.rho(); rho = thermo.rho();
rho.relax(); rho.relax();
volScalarField rUA = 1.0/UEqn().A(); volScalarField rAU = 1.0/UEqn().A();
surfaceScalarField rhorUAf("(rho*(1|A(U)))", fvc::interpolate(rho*rUA)); surfaceScalarField rhorAUf("(rho*(1|A(U)))", fvc::interpolate(rho*rAU));
U = rUA*UEqn().H(); U = rAU*UEqn().H();
UEqn.clear(); UEqn.clear();
phi = fvc::interpolate(rho)*(fvc::interpolate(U) & mesh.Sf()); phi = fvc::interpolate(rho)*(fvc::interpolate(U) & mesh.Sf());
bool closedVolume = adjustPhi(phi, U, p_rgh); bool closedVolume = adjustPhi(phi, U, p_rgh);
surfaceScalarField buoyancyPhi = rhorUAf*ghf*fvc::snGrad(rho)*mesh.magSf(); surfaceScalarField buoyancyPhi = rhorAUf*ghf*fvc::snGrad(rho)*mesh.magSf();
phi -= buoyancyPhi; phi -= buoyancyPhi;
for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++) for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
{ {
fvScalarMatrix p_rghEqn fvScalarMatrix p_rghEqn
( (
fvm::laplacian(rhorUAf, p_rgh) == fvc::div(phi) fvm::laplacian(rhorAUf, p_rgh) == fvc::div(phi)
); );
p_rghEqn.setReference(pRefCell, getRefCellValue(p_rgh, pRefCell)); p_rghEqn.setReference(pRefCell, getRefCellValue(p_rgh, pRefCell));
@ -34,7 +34,7 @@
// Correct the momentum source with the pressure gradient flux // Correct the momentum source with the pressure gradient flux
// calculated from the relaxed pressure // calculated from the relaxed pressure
U -= rUA*fvc::reconstruct((buoyancyPhi + p_rghEqn.flux())/rhorUAf); U -= rAU*fvc::reconstruct((buoyancyPhi + p_rghEqn.flux())/rhorAUf);
U.correctBoundaryConditions(); U.correctBoundaryConditions();
} }
} }

View File

@ -1,24 +1,24 @@
{ {
rho = thermo.rho(); rho = thermo.rho();
volScalarField rUA = 1.0/UEqn().A(); volScalarField rAU = 1.0/UEqn().A();
surfaceScalarField rhorUAf("(rho*(1|A(U)))", fvc::interpolate(rho*rUA)); surfaceScalarField rhorAUf("(rho*(1|A(U)))", fvc::interpolate(rho*rAU));
U = rUA*UEqn().H(); U = rAU*UEqn().H();
UEqn.clear(); UEqn.clear();
phi = fvc::interpolate(rho)*(fvc::interpolate(U) & mesh.Sf()); phi = fvc::interpolate(rho)*(fvc::interpolate(U) & mesh.Sf());
bool closedVolume = adjustPhi(phi, U, p); bool closedVolume = adjustPhi(phi, U, p);
surfaceScalarField buoyancyPhi = surfaceScalarField buoyancyPhi =
rhorUAf*fvc::interpolate(rho)*(g & mesh.Sf()); rhorAUf*fvc::interpolate(rho)*(g & mesh.Sf());
phi += buoyancyPhi; phi += buoyancyPhi;
for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++) for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
{ {
fvScalarMatrix pEqn fvScalarMatrix pEqn
( (
fvm::laplacian(rhorUAf, p) == fvc::div(phi) fvm::laplacian(rhorAUf, p) == fvc::div(phi)
); );
pEqn.setReference(pRefCell, pRefValue); pEqn.setReference(pRefCell, pRefValue);
@ -42,8 +42,8 @@
// Correct the momentum source with the pressure gradient flux // Correct the momentum source with the pressure gradient flux
// calculated from the relaxed pressure // calculated from the relaxed pressure
U += rUA*(rho*g - fvc::grad(p)); U += rAU*(rho*g - fvc::grad(p));
//U += rUA*fvc::reconstruct((buoyancyPhi - pEqn.flux())/rhorUAf); //U += rAU*fvc::reconstruct((buoyancyPhi - pEqn.flux())/rhorAUf);
U.correctBoundaryConditions(); U.correctBoundaryConditions();
} }
} }

View File

@ -5,16 +5,16 @@
rho = min(rho, rhoMax[i]); rho = min(rho, rhoMax[i]);
rho.relax(); rho.relax();
volScalarField rUA = 1.0/UEqn().A(); volScalarField rAU = 1.0/UEqn().A();
surfaceScalarField rhorUAf("(rho*(1|A(U)))", fvc::interpolate(rho*rUA)); surfaceScalarField rhorAUf("(rho*(1|A(U)))", fvc::interpolate(rho*rAU));
U = rUA*UEqn().H(); U = rAU*UEqn().H();
UEqn.clear(); UEqn.clear();
phi = fvc::interpolate(rho)*(fvc::interpolate(U) & mesh.Sf()); phi = fvc::interpolate(rho)*(fvc::interpolate(U) & mesh.Sf());
bool closedVolume = adjustPhi(phi, U, p_rgh); bool closedVolume = adjustPhi(phi, U, p_rgh);
surfaceScalarField buoyancyPhi = rhorUAf*ghf*fvc::snGrad(rho)*mesh.magSf(); surfaceScalarField buoyancyPhi = rhorAUf*ghf*fvc::snGrad(rho)*mesh.magSf();
phi -= buoyancyPhi; phi -= buoyancyPhi;
// Solve pressure // Solve pressure
@ -22,7 +22,7 @@
{ {
fvScalarMatrix p_rghEqn fvScalarMatrix p_rghEqn
( (
fvm::laplacian(rhorUAf, p_rgh) == fvc::div(phi) fvm::laplacian(rhorAUf, p_rgh) == fvc::div(phi)
); );
p_rghEqn.setReference(pRefCell, getRefCellValue(p_rgh, pRefCell)); p_rghEqn.setReference(pRefCell, getRefCellValue(p_rgh, pRefCell));
@ -39,7 +39,7 @@
// Correct the momentum source with the pressure gradient flux // Correct the momentum source with the pressure gradient flux
// calculated from the relaxed pressure // calculated from the relaxed pressure
U -= rUA*fvc::reconstruct((buoyancyPhi + p_rghEqn.flux())/rhorUAf); U -= rAU*fvc::reconstruct((buoyancyPhi + p_rghEqn.flux())/rhorAUf);
U.correctBoundaryConditions(); U.correctBoundaryConditions();
} }
} }

View File

@ -3,21 +3,21 @@
rho = thermo.rho(); rho = thermo.rho();
volScalarField rUA = 1.0/UEqn().A(); volScalarField rAU = 1.0/UEqn().A();
surfaceScalarField rhorUAf("(rho*(1|A(U)))", fvc::interpolate(rho*rUA)); surfaceScalarField rhorAUf("(rho*(1|A(U)))", fvc::interpolate(rho*rAU));
U = rUA*UEqn().H(); U = rAU*UEqn().H();
surfaceScalarField phiU surfaceScalarField phiU
( (
fvc::interpolate(rho) fvc::interpolate(rho)
*( *(
(fvc::interpolate(U) & mesh.Sf()) (fvc::interpolate(U) & mesh.Sf())
+ fvc::ddtPhiCorr(rUA, rho, U, phi) + fvc::ddtPhiCorr(rAU, rho, U, phi)
) )
); );
phi = phiU - rhorUAf*ghf*fvc::snGrad(rho)*mesh.magSf(); phi = phiU - rhorAUf*ghf*fvc::snGrad(rho)*mesh.magSf();
for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++) for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
{ {
@ -25,7 +25,7 @@
( (
fvm::ddt(psi, p_rgh) + fvc::ddt(psi, rho)*gh fvm::ddt(psi, p_rgh) + fvc::ddt(psi, rho)*gh
+ fvc::div(phi) + fvc::div(phi)
- fvm::laplacian(rhorUAf, p_rgh) - fvm::laplacian(rhorAUf, p_rgh)
); );
p_rghEqn.solve p_rghEqn.solve
@ -50,7 +50,7 @@
} }
// Correct velocity field // Correct velocity field
U += rUA*fvc::reconstruct((phi - phiU)/rhorUAf); U += rAU*fvc::reconstruct((phi - phiU)/rhorAUf);
U.correctBoundaryConditions(); U.correctBoundaryConditions();
p = p_rgh + rho*gh; p = p_rgh + rho*gh;

View File

@ -77,13 +77,13 @@ int main(int argc, char *argv[])
// --- PISO loop // --- PISO loop
volScalarField rUA = 1.0/UEqn.A(); volScalarField rAU = 1.0/UEqn.A();
for (int corr=0; corr<nCorr; corr++) for (int corr=0; corr<nCorr; corr++)
{ {
U = rUA*UEqn.H(); U = rAU*UEqn.H();
phi = (fvc::interpolate(U) & mesh.Sf()) phi = (fvc::interpolate(U) & mesh.Sf())
+ fvc::ddtPhiCorr(rUA, U, phi); + fvc::ddtPhiCorr(rAU, U, phi);
adjustPhi(phi, U, p); adjustPhi(phi, U, p);
@ -91,7 +91,7 @@ int main(int argc, char *argv[])
{ {
fvScalarMatrix pEqn fvScalarMatrix pEqn
( (
fvm::laplacian(rUA, p) == fvc::div(phi) fvm::laplacian(rAU, p) == fvc::div(phi)
); );
pEqn.setReference(pRefCell, pRefValue); pEqn.setReference(pRefCell, pRefValue);
@ -113,7 +113,7 @@ int main(int argc, char *argv[])
#include "continuityErrs.H" #include "continuityErrs.H"
U -= rUA*fvc::grad(p); U -= rAU*fvc::grad(p);
U.correctBoundaryConditions(); U.correctBoundaryConditions();
} }
@ -127,9 +127,9 @@ int main(int argc, char *argv[])
// Calculate the pressure gradient increment needed to // Calculate the pressure gradient increment needed to
// adjust the average flow-rate to the correct value // adjust the average flow-rate to the correct value
dimensionedScalar gragPplus = dimensionedScalar gragPplus =
(magUbar - magUbarStar)/rUA.weightedAverage(mesh.V()); (magUbar - magUbarStar)/rAU.weightedAverage(mesh.V());
U += flowDirection*rUA*gragPplus; U += flowDirection*rAU*gragPplus;
gradP += gragPplus; gradP += gragPplus;

View File

@ -66,11 +66,11 @@ int main(int argc, char *argv[])
for (int corr=0; corr<nCorr; corr++) for (int corr=0; corr<nCorr; corr++)
{ {
volScalarField rUA = 1.0/UEqn.A(); volScalarField rAU = 1.0/UEqn.A();
U = rUA*UEqn.H(); U = rAU*UEqn.H();
phi = (fvc::interpolate(U) & mesh.Sf()) phi = (fvc::interpolate(U) & mesh.Sf())
+ fvc::ddtPhiCorr(rUA, U, phi); + fvc::ddtPhiCorr(rAU, U, phi);
adjustPhi(phi, U, p); adjustPhi(phi, U, p);
@ -78,7 +78,7 @@ int main(int argc, char *argv[])
{ {
fvScalarMatrix pEqn fvScalarMatrix pEqn
( (
fvm::laplacian(rUA, p) == fvc::div(phi) fvm::laplacian(rAU, p) == fvc::div(phi)
); );
pEqn.setReference(pRefCell, pRefValue); pEqn.setReference(pRefCell, pRefValue);
@ -92,7 +92,7 @@ int main(int argc, char *argv[])
#include "continuityErrs.H" #include "continuityErrs.H"
U -= rUA*fvc::grad(p); U -= rAU*fvc::grad(p);
U.correctBoundaryConditions(); U.correctBoundaryConditions();
} }

View File

@ -69,11 +69,11 @@ int main(int argc, char *argv[])
for (int corr=0; corr<nCorr; corr++) for (int corr=0; corr<nCorr; corr++)
{ {
volScalarField rUA = 1.0/UEqn.A(); volScalarField rAU = 1.0/UEqn.A();
U = rUA*UEqn.H(); U = rAU*UEqn.H();
phi = (fvc::interpolate(U) & mesh.Sf()) phi = (fvc::interpolate(U) & mesh.Sf())
+ fvc::ddtPhiCorr(rUA, U, phi); + fvc::ddtPhiCorr(rAU, U, phi);
adjustPhi(phi, U, p); adjustPhi(phi, U, p);
@ -81,7 +81,7 @@ int main(int argc, char *argv[])
{ {
fvScalarMatrix pEqn fvScalarMatrix pEqn
( (
fvm::laplacian(rUA, p) == fvc::div(phi) fvm::laplacian(rAU, p) == fvc::div(phi)
); );
pEqn.setReference(pRefCell, pRefValue); pEqn.setReference(pRefCell, pRefValue);
@ -95,7 +95,7 @@ int main(int argc, char *argv[])
#include "continuityErrs.H" #include "continuityErrs.H"
U -= rUA*fvc::grad(p); U -= rAU*fvc::grad(p);
U.correctBoundaryConditions(); U.correctBoundaryConditions();
} }

View File

@ -0,0 +1,22 @@
// Solve the Momentum equation
tmp<fvVectorMatrix> UEqn
(
fvm::ddt(U)
+ fvm::div(phi, U)
+ turbulence->divDevReff(U)
);
UEqn().relax();
rAU = 1.0/UEqn().A();
if (momentumPredictor)
{
solve(UEqn() == -fvc::grad(p));
}
else
{
U = rAU*(UEqn().H() - fvc::grad(p));
U.correctBoundaryConditions();
}

View File

@ -79,11 +79,11 @@ int main(int argc, char *argv[])
for (int corr=0; corr<nCorr; corr++) for (int corr=0; corr<nCorr; corr++)
{ {
volScalarField rUA = 1.0/UEqn.A(); volScalarField rAU = 1.0/UEqn.A();
U = rUA*UEqn.H(); U = rAU*UEqn.H();
phi = (fvc::interpolate(U) & mesh.Sf()) phi = (fvc::interpolate(U) & mesh.Sf())
+ fvc::ddtPhiCorr(rUA, U, phi); + fvc::ddtPhiCorr(rAU, U, phi);
adjustPhi(phi, U, p); adjustPhi(phi, U, p);
@ -94,7 +94,7 @@ int main(int argc, char *argv[])
fvScalarMatrix pEqn fvScalarMatrix pEqn
( (
fvm::laplacian(rUA, p) == fvc::div(phi) fvm::laplacian(rAU, p) == fvc::div(phi)
); );
pEqn.setReference(pRefCell, pRefValue); pEqn.setReference(pRefCell, pRefValue);
@ -120,7 +120,7 @@ int main(int argc, char *argv[])
#include "continuityErrs.H" #include "continuityErrs.H"
U -= rUA*fvc::grad(p); U -= rAU*fvc::grad(p);
U.correctBoundaryConditions(); U.correctBoundaryConditions();
} }
} }

View File

@ -89,22 +89,22 @@ int main(int argc, char *argv[])
// --- PISO loop // --- PISO loop
for (int corr=0; corr<nCorr; corr++) for (int corr=0; corr<nCorr; corr++)
{ {
volScalarField rUA = 1.0/hUEqn.A(); volScalarField rAU = 1.0/hUEqn.A();
surfaceScalarField ghrUAf = magg*fvc::interpolate(h*rUA); surfaceScalarField ghrAUf = magg*fvc::interpolate(h*rAU);
surfaceScalarField phih0 = ghrUAf*mesh.magSf()*fvc::snGrad(h0); surfaceScalarField phih0 = ghrAUf*mesh.magSf()*fvc::snGrad(h0);
if (rotating) if (rotating)
{ {
hU = rUA*(hUEqn.H() - (F ^ hU)); hU = rAU*(hUEqn.H() - (F ^ hU));
} }
else else
{ {
hU = rUA*hUEqn.H(); hU = rAU*hUEqn.H();
} }
phi = (fvc::interpolate(hU) & mesh.Sf()) phi = (fvc::interpolate(hU) & mesh.Sf())
+ fvc::ddtPhiCorr(rUA, h, hU, phi) + fvc::ddtPhiCorr(rAU, h, hU, phi)
- phih0; - phih0;
for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++) for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
@ -113,7 +113,7 @@ int main(int argc, char *argv[])
( (
fvm::ddt(h) fvm::ddt(h)
+ fvc::div(phi) + fvc::div(phi)
- fvm::laplacian(ghrUAf, h) - fvm::laplacian(ghrAUf, h)
); );
if (ucorr < nOuterCorr-1 || corr < nCorr-1) if (ucorr < nOuterCorr-1 || corr < nCorr-1)
@ -131,7 +131,7 @@ int main(int argc, char *argv[])
} }
} }
hU -= rUA*h*magg*fvc::grad(h + h0); hU -= rAU*h*magg*fvc::grad(h + h0);
// Constrain the momentum to be in the geometry if 3D geometry // Constrain the momentum to be in the geometry if 3D geometry
if (mesh.nGeometricD() == 3) if (mesh.nGeometricD() == 3)

View File

@ -1,7 +1,7 @@
rho = thermo.rho(); rho = thermo.rho();
volScalarField rUA = 1.0/UEqn.A(); volScalarField rAU = 1.0/UEqn.A();
U = rUA*UEqn.H(); U = rAU*UEqn.H();
if (transonic) if (transonic)
{ {
@ -11,7 +11,7 @@ if (transonic)
fvc::interpolate(psi) fvc::interpolate(psi)
*( *(
(fvc::interpolate(U) & mesh.Sf()) (fvc::interpolate(U) & mesh.Sf())
+ fvc::ddtPhiCorr(rUA, rho, U, phi) + fvc::ddtPhiCorr(rAU, rho, U, phi)
) )
); );
@ -21,7 +21,7 @@ if (transonic)
( (
fvm::ddt(psi, p) fvm::ddt(psi, p)
+ fvm::div(phid, p) + fvm::div(phid, p)
- fvm::laplacian(rho*rUA, p) - fvm::laplacian(rho*rAU, p)
== ==
coalParcels.Srho() coalParcels.Srho()
); );
@ -53,7 +53,7 @@ else
fvc::interpolate(rho)* fvc::interpolate(rho)*
( (
(fvc::interpolate(U) & mesh.Sf()) (fvc::interpolate(U) & mesh.Sf())
+ fvc::ddtPhiCorr(rUA, rho, U, phi) + fvc::ddtPhiCorr(rAU, rho, U, phi)
); );
for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++) for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
@ -62,7 +62,7 @@ else
( (
fvm::ddt(psi, p) fvm::ddt(psi, p)
+ fvc::div(phi) + fvc::div(phi)
- fvm::laplacian(rho*rUA, p) - fvm::laplacian(rho*rAU, p)
== ==
coalParcels.Srho() coalParcels.Srho()
); );
@ -92,7 +92,7 @@ else
#include "rhoEqn.H" #include "rhoEqn.H"
#include "compressibleContinuityErrs.H" #include "compressibleContinuityErrs.H"
U -= rUA*fvc::grad(p); U -= rAU*fvc::grad(p);
U.correctBoundaryConditions(); U.correctBoundaryConditions();
DpDt = fvc::DDt(surfaceScalarField("phiU", phi/fvc::interpolate(rho)), p); DpDt = fvc::DDt(surfaceScalarField("phiU", phi/fvc::interpolate(rho)), p);

View File

@ -1,7 +1,7 @@
rho = thermo.rho(); rho = thermo.rho();
volScalarField rUA = 1.0/UEqn.A(); volScalarField rAU = 1.0/UEqn.A();
U = rUA*UEqn.H(); U = rAU*UEqn.H();
if (transonic) if (transonic)
{ {
@ -11,7 +11,7 @@ if (transonic)
fvc::interpolate(psi) fvc::interpolate(psi)
*( *(
(fvc::interpolate(U) & mesh.Sf()) (fvc::interpolate(U) & mesh.Sf())
+ fvc::ddtPhiCorr(rUA, rho, U, phi) + fvc::ddtPhiCorr(rAU, rho, U, phi)
) )
); );
@ -21,7 +21,7 @@ if (transonic)
( (
fvm::ddt(psi, p) fvm::ddt(psi, p)
+ fvm::div(phid, p) + fvm::div(phid, p)
- fvm::laplacian(rho*rUA, p) - fvm::laplacian(rho*rAU, p)
== ==
parcels.Srho() parcels.Srho()
+ surfaceFilm.Srho() + surfaceFilm.Srho()
@ -41,7 +41,7 @@ else
fvc::interpolate(rho) fvc::interpolate(rho)
*( *(
(fvc::interpolate(U) & mesh.Sf()) (fvc::interpolate(U) & mesh.Sf())
+ fvc::ddtPhiCorr(rUA, rho, U, phi) + fvc::ddtPhiCorr(rAU, rho, U, phi)
); );
for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++) for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
@ -50,7 +50,7 @@ else
( (
fvm::ddt(psi, p) fvm::ddt(psi, p)
+ fvc::div(phi) + fvc::div(phi)
- fvm::laplacian(rho*rUA, p) - fvm::laplacian(rho*rAU, p)
== ==
parcels.Srho() parcels.Srho()
+ surfaceFilm.Srho() + surfaceFilm.Srho()
@ -68,7 +68,7 @@ else
#include "rhoEqn.H" #include "rhoEqn.H"
#include "compressibleContinuityErrs.H" #include "compressibleContinuityErrs.H"
U -= rUA*fvc::grad(p); U -= rAU*fvc::grad(p);
U.correctBoundaryConditions(); U.correctBoundaryConditions();
DpDt = fvc::DDt(surfaceScalarField("phiU", phi/fvc::interpolate(rho)), p); DpDt = fvc::DDt(surfaceScalarField("phiU", phi/fvc::interpolate(rho)), p);

View File

@ -1,7 +1,7 @@
rho = thermo.rho(); rho = thermo.rho();
volScalarField rUA = 1.0/UEqn.A(); volScalarField rAU = 1.0/UEqn.A();
U = rUA*UEqn.H(); U = rAU*UEqn.H();
if (transonic) if (transonic)
{ {
@ -11,7 +11,7 @@ if (transonic)
fvc::interpolate(psi) fvc::interpolate(psi)
*( *(
(fvc::interpolate(U) & mesh.Sf()) (fvc::interpolate(U) & mesh.Sf())
+ fvc::ddtPhiCorr(rUA, rho, U, phi) + fvc::ddtPhiCorr(rAU, rho, U, phi)
) )
); );
@ -21,7 +21,7 @@ if (transonic)
( (
fvm::ddt(psi, p) fvm::ddt(psi, p)
+ fvm::div(phid, p) + fvm::div(phid, p)
- fvm::laplacian(rho*rUA, p) - fvm::laplacian(rho*rAU, p)
== ==
parcels.Srho() parcels.Srho()
); );
@ -40,7 +40,7 @@ else
fvc::interpolate(rho) fvc::interpolate(rho)
*( *(
(fvc::interpolate(U) & mesh.Sf()) (fvc::interpolate(U) & mesh.Sf())
+ fvc::ddtPhiCorr(rUA, rho, U, phi) + fvc::ddtPhiCorr(rAU, rho, U, phi)
); );
for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++) for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
@ -49,7 +49,7 @@ else
( (
fvm::ddt(psi, p) fvm::ddt(psi, p)
+ fvc::div(phi) + fvc::div(phi)
- fvm::laplacian(rho*rUA, p) - fvm::laplacian(rho*rAU, p)
== ==
parcels.Srho() parcels.Srho()
); );
@ -66,7 +66,7 @@ else
#include "rhoEqn.H" #include "rhoEqn.H"
#include "compressibleContinuityErrs.H" #include "compressibleContinuityErrs.H"
U -= rUA*fvc::grad(p); U -= rAU*fvc::grad(p);
U.correctBoundaryConditions(); U.correctBoundaryConditions();
DpDt = fvc::DDt(surfaceScalarField("phiU", phi/fvc::interpolate(rho)), p); DpDt = fvc::DDt(surfaceScalarField("phiU", phi/fvc::interpolate(rho)), p);

View File

@ -11,16 +11,16 @@
surfaceScalarField rhof = fvc::interpolate(rho, "rhof"); surfaceScalarField rhof = fvc::interpolate(rho, "rhof");
volScalarField rUA = 1.0/UEqn.A(); volScalarField rAU = 1.0/UEqn.A();
surfaceScalarField rUAf("rUAf", rhof*fvc::interpolate(rUA)); surfaceScalarField rAUf("rAUf", rhof*fvc::interpolate(rAU));
volVectorField HbyA = rUA*UEqn.H(); volVectorField HbyA = rAU*UEqn.H();
phiv = (fvc::interpolate(HbyA) & mesh.Sf()) phiv = (fvc::interpolate(HbyA) & mesh.Sf())
+ fvc::ddtPhiCorr(rUA, rho, U, phiv); + fvc::ddtPhiCorr(rAU, rho, U, phiv);
p.boundaryField().updateCoeffs(); p.boundaryField().updateCoeffs();
surfaceScalarField phiGradp = rUAf*mesh.magSf()*fvc::snGrad(p); surfaceScalarField phiGradp = rAUf*mesh.magSf()*fvc::snGrad(p);
phiv -= phiGradp/rhof; phiv -= phiGradp/rhof;
@ -34,7 +34,7 @@
- (rhol0 + (psil - psiv)*pSat)*fvc::ddt(gamma) - pSat*fvc::ddt(psi) - (rhol0 + (psil - psiv)*pSat)*fvc::ddt(gamma) - pSat*fvc::ddt(psi)
+ fvc::div(phiv, rho) + fvc::div(phiv, rho)
+ fvc::div(phiGradp) + fvc::div(phiGradp)
- fvm::laplacian(rUAf, p) - fvm::laplacian(rAUf, p)
); );
if (corr == nCorr-1 && nonOrth == nNonOrthCorr) if (corr == nCorr-1 && nonOrth == nNonOrthCorr)
@ -79,7 +79,7 @@
// Correct velocity // Correct velocity
U = HbyA - rUA*fvc::grad(p); U = HbyA - rAU*fvc::grad(p);
// Remove the swirl component of velocity for "wedge" cases // Remove the swirl component of velocity for "wedge" cases
if (piso.found("removeSwirl")) if (piso.found("removeSwirl"))

View File

@ -1,6 +1,6 @@
{ {
volScalarField rUA = 1.0/UEqn.A(); volScalarField rAU = 1.0/UEqn.A();
surfaceScalarField rUAf = fvc::interpolate(rUA); surfaceScalarField rAUf = fvc::interpolate(rAU);
tmp<fvScalarMatrix> p_rghEqnComp; tmp<fvScalarMatrix> p_rghEqnComp;
@ -24,27 +24,27 @@
} }
U = rUA*UEqn.H(); U = rAU*UEqn.H();
surfaceScalarField phiU surfaceScalarField phiU
( (
"phiU", "phiU",
(fvc::interpolate(U) & mesh.Sf()) (fvc::interpolate(U) & mesh.Sf())
+ fvc::ddtPhiCorr(rUA, rho, U, phi) + fvc::ddtPhiCorr(rAU, rho, U, phi)
); );
phi = phiU + phi = phiU +
( (
fvc::interpolate(interface.sigmaK())*fvc::snGrad(alpha1) fvc::interpolate(interface.sigmaK())*fvc::snGrad(alpha1)
- ghf*fvc::snGrad(rho) - ghf*fvc::snGrad(rho)
)*rUAf*mesh.magSf(); )*rAUf*mesh.magSf();
for(int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++) for(int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
{ {
fvScalarMatrix p_rghEqnIncomp fvScalarMatrix p_rghEqnIncomp
( (
fvc::div(phi) fvc::div(phi)
- fvm::laplacian(rUAf, p_rgh) - fvm::laplacian(rAUf, p_rgh)
); );
solve solve
@ -75,7 +75,7 @@
} }
} }
U += rUA*fvc::reconstruct((phi - phiU)/rUAf); U += rAU*fvc::reconstruct((phi - phiU)/rAUf);
U.correctBoundaryConditions(); U.correctBoundaryConditions();
p = max p = max

View File

@ -1,6 +1,6 @@
{ {
volScalarField rUA = 1.0/UEqn.A(); volScalarField rAU = 1.0/UEqn.A();
surfaceScalarField rUAf = fvc::interpolate(rUA); surfaceScalarField rAUf = fvc::interpolate(rAU);
tmp<fvScalarMatrix> p_rghEqnComp; tmp<fvScalarMatrix> p_rghEqnComp;
@ -24,27 +24,27 @@
} }
U = rUA*UEqn.H(); U = rAU*UEqn.H();
surfaceScalarField phiU surfaceScalarField phiU
( (
"phiU", "phiU",
(fvc::interpolate(U) & mesh.Sf()) (fvc::interpolate(U) & mesh.Sf())
+ fvc::ddtPhiCorr(rUA, rho, U, phi) + fvc::ddtPhiCorr(rAU, rho, U, phi)
); );
phi = phiU + phi = phiU +
( (
fvc::interpolate(interface.sigmaK())*fvc::snGrad(alpha1) fvc::interpolate(interface.sigmaK())*fvc::snGrad(alpha1)
- ghf*fvc::snGrad(rho) - ghf*fvc::snGrad(rho)
)*rUAf*mesh.magSf(); )*rAUf*mesh.magSf();
for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++) for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
{ {
fvScalarMatrix p_rghEqnIncomp fvScalarMatrix p_rghEqnIncomp
( (
fvc::div(phi) fvc::div(phi)
- fvm::laplacian(rUAf, p_rgh) - fvm::laplacian(rAUf, p_rgh)
); );
solve solve
@ -75,7 +75,7 @@
} }
} }
U += rUA*fvc::reconstruct((phi - phiU)/rUAf); U += rAU*fvc::reconstruct((phi - phiU)/rAUf);
U.correctBoundaryConditions(); U.correctBoundaryConditions();
p = max p = max

View File

@ -1,14 +1,14 @@
{ {
volScalarField rUA = 1.0/UEqn.A(); volScalarField rAU = 1.0/UEqn.A();
surfaceScalarField rUAf = fvc::interpolate(rUA); surfaceScalarField rAUf = fvc::interpolate(rAU);
U = rUA*UEqn.H(); U = rAU*UEqn.H();
surfaceScalarField phiU surfaceScalarField phiU
( (
"phiU", "phiU",
(fvc::interpolate(U) & mesh.Sf()) (fvc::interpolate(U) & mesh.Sf())
+ fvc::ddtPhiCorr(rUA, rho, U, phi) + fvc::ddtPhiCorr(rAU, rho, U, phi)
); );
adjustPhi(phiU, U, p_rgh); adjustPhi(phiU, U, p_rgh);
@ -17,7 +17,7 @@
( (
fvc::interpolate(interface.sigmaK())*fvc::snGrad(alpha1) fvc::interpolate(interface.sigmaK())*fvc::snGrad(alpha1)
- ghf*fvc::snGrad(rho) - ghf*fvc::snGrad(rho)
)*rUAf*mesh.magSf(); )*rAUf*mesh.magSf();
Pair<tmp<volScalarField> > vDotP = twoPhaseProperties->vDotP(); Pair<tmp<volScalarField> > vDotP = twoPhaseProperties->vDotP();
const volScalarField& vDotcP = vDotP[0](); const volScalarField& vDotcP = vDotP[0]();
@ -27,7 +27,7 @@
{ {
fvScalarMatrix p_rghEqn fvScalarMatrix p_rghEqn
( (
fvc::div(phi) - fvm::laplacian(rUAf, p_rgh) fvc::div(phi) - fvm::laplacian(rAUf, p_rgh)
- (vDotvP - vDotcP)*(pSat - rho*gh) + fvm::Sp(vDotvP - vDotcP, p_rgh) - (vDotvP - vDotcP)*(pSat - rho*gh) + fvm::Sp(vDotvP - vDotcP, p_rgh)
); );
@ -52,7 +52,7 @@
} }
} }
U += rUA*fvc::reconstruct((phi - phiU)/rUAf); U += rAU*fvc::reconstruct((phi - phiU)/rAUf);
U.correctBoundaryConditions(); U.correctBoundaryConditions();
#include "continuityErrs.H" #include "continuityErrs.H"

View File

@ -1,27 +1,27 @@
volScalarField rUA = 1.0/UEqn.A(); volScalarField rAU = 1.0/UEqn.A();
surfaceScalarField rUAf surfaceScalarField rAUf
( (
"(rho*(1|A(U)))", "(rho*(1|A(U)))",
fvc::interpolate(rho)*fvc::interpolate(rUA) fvc::interpolate(rho)*fvc::interpolate(rAU)
); );
U = rUA*UEqn.H(); U = rAU*UEqn.H();
phi = phi =
fvc::interpolate(rho) fvc::interpolate(rho)
*( *(
(fvc::interpolate(U) & mesh.Sf()) (fvc::interpolate(U) & mesh.Sf())
+ fvc::ddtPhiCorr(rUA, rho, U, phi) + fvc::ddtPhiCorr(rAU, rho, U, phi)
); );
surfaceScalarField phiU("phiU", phi); surfaceScalarField phiU("phiU", phi);
phi -= ghf*fvc::snGrad(rho)*rUAf*mesh.magSf(); phi -= ghf*fvc::snGrad(rho)*rAUf*mesh.magSf();
for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++) for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
{ {
fvScalarMatrix p_rghEqn fvScalarMatrix p_rghEqn
( (
fvm::laplacian(rUAf, p_rgh) == fvc::ddt(rho) + fvc::div(phi) fvm::laplacian(rAUf, p_rgh) == fvc::ddt(rho) + fvc::div(phi)
); );
p_rghEqn.setReference(pRefCell, getRefCellValue(p_rgh, pRefCell)); p_rghEqn.setReference(pRefCell, getRefCellValue(p_rgh, pRefCell));
@ -49,5 +49,5 @@ if (p_rgh.needReference())
#include "rhoEqn.H" #include "rhoEqn.H"
#include "compressibleContinuityErrs.H" #include "compressibleContinuityErrs.H"
U += rUA*fvc::reconstruct((phi - phiU)/rUAf); U += rAU*fvc::reconstruct((phi - phiU)/rAUf);
U.correctBoundaryConditions(); U.correctBoundaryConditions();

View File

@ -102,11 +102,11 @@ int main(int argc, char *argv[])
for (int corr=0; corr<nCorr; corr++) for (int corr=0; corr<nCorr; corr++)
{ {
volScalarField rUA = 1.0/UEqn.A(); volScalarField rAU = 1.0/UEqn.A();
U = rUA*UEqn.H(); U = rAU*UEqn.H();
phi = (fvc::interpolate(U) & mesh.Sf()) phi = (fvc::interpolate(U) & mesh.Sf())
+ fvc::ddtPhiCorr(rUA, U, phi); + fvc::ddtPhiCorr(rAU, U, phi);
adjustPhi(phi, U, p); adjustPhi(phi, U, p);
@ -114,7 +114,7 @@ int main(int argc, char *argv[])
{ {
fvScalarMatrix pEqn fvScalarMatrix pEqn
( (
fvm::laplacian(rUA, p) == fvc::div(phi) fvm::laplacian(rAU, p) == fvc::div(phi)
); );
pEqn.setReference(pRefCell, pRefValue); pEqn.setReference(pRefCell, pRefValue);
@ -128,7 +128,7 @@ int main(int argc, char *argv[])
# include "continuityErrs.H" # include "continuityErrs.H"
U -= rUA*fvc::grad(p); U -= rAU*fvc::grad(p);
U.correctBoundaryConditions(); U.correctBoundaryConditions();
} }

View File

@ -153,6 +153,12 @@ int main(int argc, char *argv[])
"internalFacesOnly", "internalFacesOnly",
"do not convert boundary faces" "do not convert boundary faces"
); );
argList::addBoolOption
(
"updateFields",
"update fields to include new patches:"
" NOTE: updated field values may need to be edited"
);
#include "setRootCase.H" #include "setRootCase.H"
#include "createTime.H" #include "createTime.H"
@ -235,39 +241,45 @@ int main(int argc, char *argv[])
IOobjectList objects(mesh, runTime.timeName()); IOobjectList objects(mesh, runTime.timeName());
// Read vol fields. // Read vol fields.
if (args.optionFound("updateFields"))
{
Info<< "Reading geometric fields" << nl << endl;
PtrList<volScalarField> vsFlds;
ReadFields(mesh, objects, vsFlds);
PtrList<volScalarField> vsFlds; PtrList<volVectorField> vvFlds;
ReadFields(mesh, objects, vsFlds); ReadFields(mesh, objects, vvFlds);
PtrList<volVectorField> vvFlds; PtrList<volSphericalTensorField> vstFlds;
ReadFields(mesh, objects, vvFlds); ReadFields(mesh, objects, vstFlds);
PtrList<volSphericalTensorField> vstFlds; PtrList<volSymmTensorField> vsymtFlds;
ReadFields(mesh, objects, vstFlds); ReadFields(mesh, objects, vsymtFlds);
PtrList<volSymmTensorField> vsymtFlds; PtrList<volTensorField> vtFlds;
ReadFields(mesh, objects, vsymtFlds); ReadFields(mesh, objects, vtFlds);
PtrList<volTensorField> vtFlds; // Read surface fields.
ReadFields(mesh, objects, vtFlds);
// Read surface fields. PtrList<surfaceScalarField> ssFlds;
ReadFields(mesh, objects, ssFlds);
PtrList<surfaceScalarField> ssFlds; PtrList<surfaceVectorField> svFlds;
ReadFields(mesh, objects, ssFlds); ReadFields(mesh, objects, svFlds);
PtrList<surfaceVectorField> svFlds; PtrList<surfaceSphericalTensorField> sstFlds;
ReadFields(mesh, objects, svFlds); ReadFields(mesh, objects, sstFlds);
PtrList<surfaceSphericalTensorField> sstFlds; PtrList<surfaceSymmTensorField> ssymtFlds;
ReadFields(mesh, objects, sstFlds); ReadFields(mesh, objects, ssymtFlds);
PtrList<surfaceSymmTensorField> ssymtFlds;
ReadFields(mesh, objects, ssymtFlds);
PtrList<surfaceTensorField> stFlds;
ReadFields(mesh, objects, stFlds);
PtrList<surfaceTensorField> stFlds;
ReadFields(mesh, objects, stFlds);
}
else
{
Info<< "Not updating geometric fields" << nl << endl;
}
// Mesh change container // Mesh change container
polyTopoChange meshMod(mesh); polyTopoChange meshMod(mesh);

View File

@ -161,13 +161,13 @@ Foam::pressureGradientExplicitSource::Su() const
void Foam::pressureGradientExplicitSource::update() void Foam::pressureGradientExplicitSource::update()
{ {
const volScalarField& rUA = const volScalarField& rAU =
mesh_.lookupObject<volScalarField>("(1|A(" + U_.name() + "))"); mesh_.lookupObject<volScalarField>("(1|A(" + U_.name() + "))");
// Integrate flow variables over cell set // Integrate flow variables over cell set
scalar volTot = 0.0; scalar volTot = 0.0;
scalar magUbarAve = 0.0; scalar magUbarAve = 0.0;
scalar rUAave = 0.0; scalar rAUave = 0.0;
forAllConstIter(cellSet, selectedCellSet_, iter) forAllConstIter(cellSet, selectedCellSet_, iter)
{ {
label cellI = iter.key(); label cellI = iter.key();
@ -176,27 +176,27 @@ void Foam::pressureGradientExplicitSource::update()
volTot += volCell; volTot += volCell;
magUbarAve += (flowDir_ & U_[cellI])*volCell; magUbarAve += (flowDir_ & U_[cellI])*volCell;
rUAave += rUA[cellI]*volCell; rAUave += rAU[cellI]*volCell;
} }
// Collect across all processors // Collect across all processors
reduce(volTot, sumOp<scalar>()); reduce(volTot, sumOp<scalar>());
reduce(magUbarAve, sumOp<scalar>()); reduce(magUbarAve, sumOp<scalar>());
reduce(rUAave, sumOp<scalar>()); reduce(rAUave, sumOp<scalar>());
// Volume averages // Volume averages
magUbarAve /= volTot; magUbarAve /= volTot;
rUAave /= volTot; rAUave /= volTot;
// Calculate the pressure gradient increment needed to adjust the average // Calculate the pressure gradient increment needed to adjust the average
// flow-rate to the desired value // flow-rate to the desired value
scalar gradPplus = (mag(Ubar_) - magUbarAve)/rUAave; scalar gradPplus = (mag(Ubar_) - magUbarAve)/rAUave;
// Apply correction to velocity field // Apply correction to velocity field
forAllConstIter(cellSet, selectedCellSet_, iter) forAllConstIter(cellSet, selectedCellSet_, iter)
{ {
label cellI = iter.key(); label cellI = iter.key();
U_[cellI] += flowDir_*rUA[cellI]*gradPplus; U_[cellI] += flowDir_*rAU[cellI]*gradPplus;
} }
// Update pressure gradient // Update pressure gradient

View File

@ -478,10 +478,10 @@ void Foam::surfaceFilmModels::kinematicSingleLayer::solveThickness
Info<< "kinematicSingleLayer::solveThickness()" << endl; Info<< "kinematicSingleLayer::solveThickness()" << endl;
} }
volScalarField rUA = 1.0/UEqn.A(); volScalarField rAU = 1.0/UEqn.A();
U_ = rUA*UEqn.H(); U_ = rAU*UEqn.H();
surfaceScalarField deltarUAf = fvc::interpolate(delta_*rUA); surfaceScalarField deltarAUf = fvc::interpolate(delta_*rAU);
surfaceScalarField rhof = fvc::interpolate(rho_); surfaceScalarField rhof = fvc::interpolate(rho_);
surfaceScalarField phiAdd surfaceScalarField phiAdd
@ -500,13 +500,13 @@ void Foam::surfaceFilmModels::kinematicSingleLayer::solveThickness
( (
"phid", "phid",
(fvc::interpolate(U_*rho_) & filmRegion_.Sf()) (fvc::interpolate(U_*rho_) & filmRegion_.Sf())
- deltarUAf*phiAdd*rhof - deltarAUf*phiAdd*rhof
); );
constrainFilmField(phid, 0.0); constrainFilmField(phid, 0.0);
surfaceScalarField ddrhorUAppf = surfaceScalarField ddrhorAUppf =
fvc::interpolate(delta_)*deltarUAf*rhof*fvc::interpolate(pp); fvc::interpolate(delta_)*deltarAUf*rhof*fvc::interpolate(pp);
// constrainFilmField(ddrhorUAppf, 0.0); // constrainFilmField(ddrhorAUppf, 0.0);
for (int nonOrth=0; nonOrth<=nNonOrthCorr_; nonOrth++) for (int nonOrth=0; nonOrth<=nNonOrthCorr_; nonOrth++)
{ {
@ -515,7 +515,7 @@ void Foam::surfaceFilmModels::kinematicSingleLayer::solveThickness
( (
fvm::ddt(rho_, delta_) fvm::ddt(rho_, delta_)
+ fvm::div(phid, delta_) + fvm::div(phid, delta_)
- fvm::laplacian(ddrhorUAppf, delta_) - fvm::laplacian(ddrhorAUppf, delta_)
== ==
rhoSp_ rhoSp_
); );
@ -537,7 +537,7 @@ void Foam::surfaceFilmModels::kinematicSingleLayer::solveThickness
delta_.max(0.0); delta_.max(0.0);
// Update U field // Update U field
U_ -= fvc::reconstruct(deltarUAf*phiAdd); U_ -= fvc::reconstruct(deltarAUf*phiAdd);
// Remove any patch-normal components of velocity // Remove any patch-normal components of velocity
U_ -= nHat_*(nHat_ & U_); U_ -= nHat_*(nHat_ & U_);

View File

@ -1,7 +1,7 @@
rho = thermo.rho(); rho = thermo.rho();
volScalarField rUA = 1.0/UEqn.A(); volScalarField rAU = 1.0/UEqn.A();
U = rUA*UEqn.H(); U = rAU*UEqn.H();
if (transonic) if (transonic)
{ {
@ -11,7 +11,7 @@ if (transonic)
fvc::interpolate(psi) fvc::interpolate(psi)
*( *(
(fvc::interpolate(U) & mesh.Sf()) (fvc::interpolate(U) & mesh.Sf())
+ fvc::ddtPhiCorr(rUA, rho, U, phi) + fvc::ddtPhiCorr(rAU, rho, U, phi)
) )
); );
@ -21,7 +21,7 @@ if (transonic)
( (
fvm::ddt(psi, p) fvm::ddt(psi, p)
+ fvm::div(phid, p) + fvm::div(phid, p)
- fvm::laplacian(rho*rUA, p) - fvm::laplacian(rho*rAU, p)
); );
pEqn.solve(); pEqn.solve();
@ -38,7 +38,7 @@ else
fvc::interpolate(rho) fvc::interpolate(rho)
*( *(
(fvc::interpolate(U) & mesh.Sf()) (fvc::interpolate(U) & mesh.Sf())
+ fvc::ddtPhiCorr(rUA, rho, U, phi) + fvc::ddtPhiCorr(rAU, rho, U, phi)
); );
for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++) for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
@ -47,7 +47,7 @@ else
( (
fvm::ddt(psi, p) fvm::ddt(psi, p)
+ fvc::div(phi) + fvc::div(phi)
- fvm::laplacian(rho*rUA, p) - fvm::laplacian(rho*rAU, p)
); );
pEqn.solve(); pEqn.solve();
@ -62,7 +62,7 @@ else
#include "rhoEqn.H" #include "rhoEqn.H"
#include "compressibleContinuityErrs.H" #include "compressibleContinuityErrs.H"
U -= rUA*fvc::grad(p); U -= rAU*fvc::grad(p);
U.correctBoundaryConditions(); U.correctBoundaryConditions();
DpDt = fvc::DDt(surfaceScalarField("phiU", phi/fvc::interpolate(rho)), p); DpDt = fvc::DDt(surfaceScalarField("phiU", phi/fvc::interpolate(rho)), p);

View File

@ -44,7 +44,7 @@ laplacianSchemes
laplacian(nuf,rhoU) Gauss linear corrected; laplacian(nuf,rhoU) Gauss linear corrected;
laplacian(muEff,U) Gauss linear corrected; laplacian(muEff,U) Gauss linear corrected;
laplacian(rrhoUAf,p) Gauss linear corrected; laplacian(rrhoUAf,p) Gauss linear corrected;
laplacian(rUAf,p) Gauss linear corrected; laplacian(rAUf,p) Gauss linear corrected;
laplacian(DkEff,k) Gauss linear corrected; laplacian(DkEff,k) Gauss linear corrected;
laplacian(1,p) Gauss linear corrected; laplacian(1,p) Gauss linear corrected;
} }

View File

@ -44,7 +44,7 @@ laplacianSchemes
laplacian(nuf,rhoU) Gauss linear corrected; laplacian(nuf,rhoU) Gauss linear corrected;
laplacian(muEff,U) Gauss linear corrected; laplacian(muEff,U) Gauss linear corrected;
laplacian(rrhoUAf,p) Gauss linear corrected; laplacian(rrhoUAf,p) Gauss linear corrected;
laplacian(rUAf,p) Gauss linear corrected; laplacian(rAUf,p) Gauss linear corrected;
laplacian(DkEff,k) Gauss linear corrected; laplacian(DkEff,k) Gauss linear corrected;
laplacian(1,p) Gauss linear corrected; laplacian(1,p) Gauss linear corrected;
} }

View File

@ -45,7 +45,7 @@ laplacianSchemes
laplacian(nuf,rhoU) Gauss linear corrected; laplacian(nuf,rhoU) Gauss linear corrected;
laplacian(muEff,U) Gauss linear corrected; laplacian(muEff,U) Gauss linear corrected;
laplacian(rrhoUAf,p) Gauss linear corrected; laplacian(rrhoUAf,p) Gauss linear corrected;
laplacian(rUAf,p) Gauss linear corrected; laplacian(rAUf,p) Gauss linear corrected;
laplacian(DomegaEff,omega) Gauss linear corrected; laplacian(DomegaEff,omega) Gauss linear corrected;
laplacian(DkEff,k) Gauss linear corrected; laplacian(DkEff,k) Gauss linear corrected;
laplacian(1,p) Gauss linear corrected; laplacian(1,p) Gauss linear corrected;