Commit 13866c28 authored by Henry Weller's avatar Henry Weller Committed by Andrew Heather
Browse files

rhoPimpleFoam: Updated transonic option to be consistent with sonicFoam

Improves stability on start-up and allows running at slightly larger time-steps.
parent 3861d8da
rho = thermo.rho();
volScalarField rAU(1.0/UEqn.A());
surfaceScalarField rhorAUf("rhorAUf", fvc::interpolate(rho*rAU));
volVectorField HbyA(constrainHbyA(rAU*UEqn.H(), U, p));
......@@ -10,7 +12,7 @@ if (pimple.nCorrPISO() <= 1)
surfaceScalarField phiHbyA
(
"phiHbyA",
fvc::flux(rho*HbyA)
fvc::interpolate(rho)*fvc::flux(HbyA)
+ rhorAUf*fvc::ddtCorr(rho, U, phi)
);
......@@ -26,7 +28,8 @@ if (pimple.transonic())
"phid",
(fvc::interpolate(psi)/fvc::interpolate(rho))*phiHbyA
);
phiHbyA -= fvc::interpolate(p)*phid;
phiHbyA -= fvc::interpolate(psi*p)*phiHbyA/fvc::interpolate(rho);
while (pimple.correctNonOrthogonal())
{
......@@ -84,9 +87,11 @@ U.correctBoundaryConditions();
fvOptions.correct(U);
K = 0.5*magSqr(U);
pressureControl.limit(p);
p.correctBoundaryConditions();
rho = thermo.rho();
if (pressureControl.limit(p))
{
p.correctBoundaryConditions();
rho = thermo.rho();
}
if (!pimple.transonic())
{
......
rho = thermo.rho();
volScalarField rAU(1.0/UEqn.A());
volScalarField rAtU(1.0/(1.0/rAU - UEqn.H1()));
volVectorField HbyA(constrainHbyA(rAU*UEqn.H(), U, p));
......@@ -11,7 +13,7 @@ surfaceScalarField phiHbyA
(
"phiHbyA",
(
fvc::flux(rho*HbyA)
fvc::interpolate(rho)*fvc::flux(HbyA)
+ fvc::interpolate(rho*rAU)*fvc::ddtCorr(rho, U, phi)
)
);
......@@ -33,7 +35,7 @@ if (pimple.transonic())
phiHbyA +=
fvc::interpolate(rho*(rAtU - rAU))*fvc::snGrad(p)*mesh.magSf()
- fvc::interpolate(p)*phid;
- fvc::interpolate(psi*p)*phiHbyA/fvc::interpolate(rho);
HbyA -= (rAU - rAtU)*fvc::grad(p);
......@@ -96,9 +98,11 @@ U.correctBoundaryConditions();
fvOptions.correct(U);
K = 0.5*magSqr(U);
pressureControl.limit(p);
p.correctBoundaryConditions();
rho = thermo.rho();
if (pressureControl.limit(p))
{
p.correctBoundaryConditions();
rho = thermo.rho();
}
if (!pimple.transonic())
{
......
{
volScalarField rAU(1.0/UEqn.A());
surfaceScalarField rhorAUf("rhorAUf", fvc::interpolate(rho*rAU));
volVectorField HbyA(constrainHbyA(rAU*UEqn.H(), U, p));
tUEqn.clear();
volScalarField rAU(1.0/UEqn.A());
surfaceScalarField rhorAUf("rhorAUf", fvc::interpolate(rho*rAU));
volVectorField HbyA(constrainHbyA(rAU*UEqn.H(), U, p));
tUEqn.clear();
bool closedVolume = false;
bool closedVolume = false;
surfaceScalarField phiHbyA("phiHbyA", fvc::interpolate(rho)*fvc::flux(HbyA));
MRF.makeRelative(fvc::interpolate(rho), phiHbyA);
surfaceScalarField phiHbyA("phiHbyA", fvc::flux(rho*HbyA));
MRF.makeRelative(fvc::interpolate(rho), phiHbyA);
// Update the pressure BCs to ensure flux consistency
constrainPressure(p, rho, U, phiHbyA, rhorAUf, MRF);
if (simple.transonic())
{
surfaceScalarField phid
(
"phid",
(fvc::interpolate(psi)/fvc::interpolate(rho))*phiHbyA
);
// Update the pressure BCs to ensure flux consistency
constrainPressure(p, rho, U, phiHbyA, rhorAUf, MRF);
phiHbyA -= fvc::interpolate(psi*p)*phiHbyA/fvc::interpolate(rho);
if (simple.transonic())
while (simple.correctNonOrthogonal())
{
surfaceScalarField phid
fvScalarMatrix pEqn
(
fvc::div(phiHbyA)
+ fvm::div(phid, p)
- fvm::laplacian(rhorAUf, p)
==
fvOptions(psi, p, rho.name())
);
// Relax the pressure equation to ensure diagonal-dominance
pEqn.relax();
pEqn.setReference
(
"phid",
(fvc::interpolate(psi)/fvc::interpolate(rho))*phiHbyA
pressureControl.refCell(),
pressureControl.refValue()
);
phiHbyA -= fvc::interpolate(p)*phid;
while (simple.correctNonOrthogonal())
pEqn.solve();
if (simple.finalNonOrthogonalIter())
{
fvScalarMatrix pEqn
(
fvc::div(phiHbyA)
+ fvm::div(phid, p)
- fvm::laplacian(rhorAUf, p)
==
fvOptions(psi, p, rho.name())
);
// Relax the pressure equation to ensure diagonal-dominance
pEqn.relax();
pEqn.setReference
(
pressureControl.refCell(),
pressureControl.refValue()
);
pEqn.solve();
if (simple.finalNonOrthogonalIter())
{
phi = phiHbyA + pEqn.flux();
}
phi = phiHbyA + pEqn.flux();
}
}
else
}
else
{
closedVolume = adjustPhi(phiHbyA, U, p);
while (simple.correctNonOrthogonal())
{
closedVolume = adjustPhi(phiHbyA, U, p);
fvScalarMatrix pEqn
(
fvc::div(phiHbyA)
- fvm::laplacian(rhorAUf, p)
==
fvOptions(psi, p, rho.name())
);
while (simple.correctNonOrthogonal())
pEqn.setReference
(
pressureControl.refCell(),
pressureControl.refValue()
);
pEqn.solve();
if (simple.finalNonOrthogonalIter())
{
fvScalarMatrix pEqn
(
fvc::div(phiHbyA)
- fvm::laplacian(rhorAUf, p)
==
fvOptions(psi, p, rho.name())
);
pEqn.setReference
(
pressureControl.refCell(),
pressureControl.refValue()
);
pEqn.solve();
if (simple.finalNonOrthogonalIter())
{
phi = phiHbyA + pEqn.flux();
}
phi = phiHbyA + pEqn.flux();
}
}
}
#include "incompressible/continuityErrs.H"
#include "incompressible/continuityErrs.H"
// Explicitly relax pressure for momentum corrector
p.relax();
// Explicitly relax pressure for momentum corrector
p.relax();
U = HbyA - rAU*fvc::grad(p);
U.correctBoundaryConditions();
fvOptions.correct(U);
U = HbyA - rAU*fvc::grad(p);
U.correctBoundaryConditions();
fvOptions.correct(U);
pressureControl.limit(p);
bool pLimited = pressureControl.limit(p);
// For closed-volume cases adjust the pressure and density levels
// to obey overall mass continuity
if (closedVolume)
{
p += (initialMass - fvc::domainIntegrate(psi*p))
/fvc::domainIntegrate(psi);
}
// For closed-volume cases adjust the pressure and density levels
// to obey overall mass continuity
if (closedVolume)
{
p += (initialMass - fvc::domainIntegrate(psi*p))
/fvc::domainIntegrate(psi);
}
if (pLimited || closedVolume)
{
p.correctBoundaryConditions();
}
rho = thermo.rho();
rho = thermo.rho();
if (!simple.transonic())
{
rho.relax();
}
if (!simple.transonic())
{
rho.relax();
}
rho = thermo.rho();
volScalarField rAU(1.0/UEqn.A());
volScalarField rAtU(1.0/(1.0/rAU - UEqn.H1()));
volVectorField HbyA(constrainHbyA(rAU*UEqn.H(), U, p));
......@@ -5,7 +7,7 @@ tUEqn.clear();
bool closedVolume = false;
surfaceScalarField phiHbyA("phiHbyA", fvc::flux(rho*HbyA));
surfaceScalarField phiHbyA("phiHbyA", fvc::interpolate(rho)*fvc::flux(HbyA));
MRF.makeRelative(fvc::interpolate(rho), phiHbyA);
volScalarField rhorAtU("rhorAtU", rho*rAtU);
......@@ -23,7 +25,7 @@ if (simple.transonic())
phiHbyA +=
fvc::interpolate(rho*(rAtU - rAU))*fvc::snGrad(p)*mesh.magSf()
- fvc::interpolate(p)*phid;
- fvc::interpolate(psi*p)*phiHbyA/fvc::interpolate(rho);
HbyA -= (rAU - rAtU)*fvc::grad(p);
......@@ -98,7 +100,7 @@ U = HbyA - rAtU*fvc::grad(p);
U.correctBoundaryConditions();
fvOptions.correct(U);
pressureControl.limit(p);
bool pLimited = pressureControl.limit(p);
// For closed-volume cases adjust the pressure and density levels
// to obey overall mass continuity
......@@ -108,9 +110,11 @@ if (closedVolume)
/fvc::domainIntegrate(psi);
}
p.correctBoundaryConditions();
if (pLimited || closedVolume)
{
p.correctBoundaryConditions();
}
// Recalculate density from the relaxed pressure
rho = thermo.rho();
if (!simple.transonic())
......
......@@ -207,14 +207,24 @@ Foam::pressureControl::pressureControl
// * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * * //
void Foam::pressureControl::limit(volScalarField& p) const
bool Foam::pressureControl::limit(volScalarField& p) const
{
Info<< "pressureControl: p max/min "
<< max(p).value() << " "
<< min(p).value() << endl;
scalar pMax = max(p).value();
scalar pMin = min(p).value();
p = max(p, pMin_);
p = min(p, pMax_);
if (pMax > pMax_.value() || pMin < pMin_.value())
{
Info<< "pressureControl: p max/min " << pMax << " " << pMin << endl;
p = max(p, pMin_);
p = min(p, pMax_);
return true;
}
else
{
return false;
}
}
......
......@@ -89,8 +89,8 @@ public:
//- Return the pressure reference level
inline scalar refValue() const;
//- Limit the pressure
void limit(volScalarField& p) const;
//- Limit the pressure if necessary and return true if so
bool limit(volScalarField& p) const;
};
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment