multiphaseEulerFoam: Support for frozen flow switch
multiphaseEulerFoam now supports the "frozenFlow" setting in system/fvSolution/PIMPLE. With this switch on, the pressure-velocity system is not solved. Energy and species transport are simulated using fixed velocity and pressure fields. Non-transport effects on the phase fraction, such as mass transfer or fvOption sources, are also included. Note that this setting does not enforce conservation. In general, any processes that alter the phase fraction have to be set up carefully in order to avoid generating an inconsistent set of phase fractions. This switch has been enabled in two zero-dimensional test cases which are designed to operate at a constant pressure. The "frozenFlow" switch now enforces this constant pressure directly. Previously the pressure equation was being solved, but the system was not well posed, and the cases were failing to run as a result.
This commit is contained in:
Will Bainbridge
@ -101,51 +101,69 @@ int main(int argc, char *argv[])
|
||||
runTime++;
|
||||
Info<< "Time = " << runTime.timeName() << nl << endl;
|
||||
|
||||
// --- Pressure-velocity PIMPLE corrector loop
|
||||
while (pimple.loop())
|
||||
if (pimple.frozenFlow())
|
||||
{
|
||||
if (pimple.firstPimpleIter() || moveMeshOuterCorrectors)
|
||||
fluid.solve(rAUs, rAUfs);
|
||||
fluid.correct();
|
||||
fluid.correctContinuityError();
|
||||
|
||||
#include "YEqns.H"
|
||||
#include "EEqns.H"
|
||||
#include "pEqnComps.H"
|
||||
|
||||
forAll(phases, phasei)
|
||||
{
|
||||
phases[phasei].divU(-pEqnComps[phasei] & p_rgh);
|
||||
}
|
||||
}
|
||||
else
|
||||
{
|
||||
// --- Pressure-velocity PIMPLE corrector loop
|
||||
while (pimple.loop())
|
||||
{
|
||||
mesh.update();
|
||||
|
||||
if (mesh.changing())
|
||||
if (pimple.firstPimpleIter() || moveMeshOuterCorrectors)
|
||||
{
|
||||
gh = (g & mesh.C()) - ghRef;
|
||||
ghf = (g & mesh.Cf()) - ghRef;
|
||||
mesh.update();
|
||||
|
||||
fluid.meshUpdate();
|
||||
|
||||
if (checkMeshCourantNo)
|
||||
if (mesh.changing())
|
||||
{
|
||||
#include "meshCourantNo.H"
|
||||
gh = (g & mesh.C()) - ghRef;
|
||||
ghf = (g & mesh.Cf()) - ghRef;
|
||||
|
||||
fluid.meshUpdate();
|
||||
|
||||
if (checkMeshCourantNo)
|
||||
{
|
||||
#include "meshCourantNo.H"
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
fluid.solve(rAUs, rAUfs);
|
||||
fluid.correct();
|
||||
fluid.correctContinuityError();
|
||||
fluid.solve(rAUs, rAUfs);
|
||||
fluid.correct();
|
||||
fluid.correctContinuityError();
|
||||
|
||||
#include "YEqns.H"
|
||||
#include "YEqns.H"
|
||||
|
||||
if (faceMomentum)
|
||||
{
|
||||
#include "pUf/UEqns.H"
|
||||
#include "EEqns.H"
|
||||
#include "pUf/pEqn.H"
|
||||
}
|
||||
else
|
||||
{
|
||||
#include "pU/UEqns.H"
|
||||
#include "EEqns.H"
|
||||
#include "pU/pEqn.H"
|
||||
}
|
||||
if (faceMomentum)
|
||||
{
|
||||
#include "pUf/UEqns.H"
|
||||
#include "EEqns.H"
|
||||
#include "pUf/pEqn.H"
|
||||
}
|
||||
else
|
||||
{
|
||||
#include "pU/UEqns.H"
|
||||
#include "EEqns.H"
|
||||
#include "pU/pEqn.H"
|
||||
}
|
||||
|
||||
fluid.correctKinematics();
|
||||
fluid.correctKinematics();
|
||||
|
||||
if (pimple.turbCorr())
|
||||
{
|
||||
fluid.correctTurbulence();
|
||||
if (pimple.turbCorr())
|
||||
{
|
||||
fluid.correctTurbulence();
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
@ -0,0 +1,75 @@
|
||||
PtrList<fvScalarMatrix> pEqnComps(phases.size());
|
||||
|
||||
{
|
||||
PtrList<volScalarField> dmdts(fluid.dmdts());
|
||||
|
||||
forAll(phases, phasei)
|
||||
{
|
||||
phaseModel& phase = phases[phasei];
|
||||
const volScalarField& alpha = phase;
|
||||
volScalarField& rho = phase.thermoRef().rho();
|
||||
|
||||
pEqnComps.set(phasei, new fvScalarMatrix(p_rgh, dimVolume/dimTime));
|
||||
fvScalarMatrix& pEqnComp = pEqnComps[phasei];
|
||||
|
||||
// Density variation
|
||||
if (!phase.isochoric() || !phase.pure())
|
||||
{
|
||||
pEqnComp +=
|
||||
(
|
||||
fvc::ddt(alpha, rho) + fvc::div(phase.alphaRhoPhi())
|
||||
- fvc::Sp(fvc::ddt(alpha) + fvc::div(phase.alphaPhi()), rho)
|
||||
)/rho;
|
||||
}
|
||||
|
||||
// Mesh dilatation correction
|
||||
if (mesh.moving())
|
||||
{
|
||||
pEqnComp += fvc::div(mesh.phi())*alpha;
|
||||
}
|
||||
|
||||
// Compressibility
|
||||
if (!phase.incompressible())
|
||||
{
|
||||
if (pimple.transonic())
|
||||
{
|
||||
const surfaceScalarField phid
|
||||
(
|
||||
IOobject::groupName("phid", phase.name()),
|
||||
fvc::interpolate(phase.thermo().psi())*phase.phi()
|
||||
);
|
||||
|
||||
pEqnComp +=
|
||||
correction
|
||||
(
|
||||
(alpha/rho)*
|
||||
(
|
||||
phase.thermo().psi()*fvm::ddt(p_rgh)
|
||||
+ fvm::div(phid, p_rgh)
|
||||
- fvm::Sp(fvc::div(phid), p_rgh)
|
||||
)
|
||||
);
|
||||
|
||||
pEqnComps[phasei].relax();
|
||||
}
|
||||
else
|
||||
{
|
||||
pEqnComp +=
|
||||
(alpha*phase.thermo().psi()/rho)
|
||||
*correction(fvm::ddt(p_rgh));
|
||||
}
|
||||
}
|
||||
|
||||
// Option sources
|
||||
if (fvOptions.appliesToField(rho.name()))
|
||||
{
|
||||
pEqnComp -= (fvOptions(alpha, rho) & rho)/rho;
|
||||
}
|
||||
|
||||
// Mass transfer
|
||||
if (dmdts.set(phasei))
|
||||
{
|
||||
pEqnComp -= dmdts[phasei]/rho;
|
||||
}
|
||||
}
|
||||
}
|
||||
@ -231,77 +231,7 @@ while (pimple.correct())
|
||||
}
|
||||
|
||||
// Compressible pressure equations
|
||||
PtrList<fvScalarMatrix> pEqnComps(phases.size());
|
||||
PtrList<volScalarField> dmdts(fluid.dmdts());
|
||||
forAll(phases, phasei)
|
||||
{
|
||||
phaseModel& phase = phases[phasei];
|
||||
const volScalarField& alpha = phase;
|
||||
volScalarField& rho = phase.thermoRef().rho();
|
||||
|
||||
pEqnComps.set(phasei, new fvScalarMatrix(p_rgh, dimVolume/dimTime));
|
||||
fvScalarMatrix& pEqnComp = pEqnComps[phasei];
|
||||
|
||||
// Density variation
|
||||
if (!phase.isochoric() || !phase.pure())
|
||||
{
|
||||
pEqnComp +=
|
||||
(
|
||||
fvc::ddt(alpha, rho) + fvc::div(phase.alphaRhoPhi())
|
||||
- fvc::Sp(fvc::ddt(alpha) + fvc::div(phase.alphaPhi()), rho)
|
||||
)/rho;
|
||||
}
|
||||
|
||||
// Mesh dilatation correction
|
||||
if (mesh.moving())
|
||||
{
|
||||
pEqnComp += fvc::div(mesh.phi())*alpha;
|
||||
}
|
||||
|
||||
// Compressibility
|
||||
if (!phase.incompressible())
|
||||
{
|
||||
if (pimple.transonic())
|
||||
{
|
||||
const surfaceScalarField phid
|
||||
(
|
||||
IOobject::groupName("phid", phase.name()),
|
||||
fvc::interpolate(phase.thermo().psi())*phase.phi()
|
||||
);
|
||||
|
||||
pEqnComp +=
|
||||
correction
|
||||
(
|
||||
(alpha/rho)*
|
||||
(
|
||||
phase.thermo().psi()*fvm::ddt(p_rgh)
|
||||
+ fvm::div(phid, p_rgh)
|
||||
- fvm::Sp(fvc::div(phid), p_rgh)
|
||||
)
|
||||
);
|
||||
|
||||
pEqnComps[phasei].relax();
|
||||
}
|
||||
else
|
||||
{
|
||||
pEqnComp +=
|
||||
(alpha*phase.thermo().psi()/rho)
|
||||
*correction(fvm::ddt(p_rgh));
|
||||
}
|
||||
}
|
||||
|
||||
// Option sources
|
||||
if (fvOptions.appliesToField(rho.name()))
|
||||
{
|
||||
pEqnComp -= (fvOptions(alpha, rho) & rho)/rho;
|
||||
}
|
||||
|
||||
// Mass transfer
|
||||
if (dmdts.set(phasei))
|
||||
{
|
||||
pEqnComp -= dmdts[phasei]/rho;
|
||||
}
|
||||
}
|
||||
#include "pEqnComps.H"
|
||||
|
||||
// Cache p prior to solve for density update
|
||||
volScalarField p_rgh_0(p_rgh);
|
||||
|
||||
@ -222,82 +222,7 @@ while (pimple.correct())
|
||||
}
|
||||
|
||||
// Compressible pressure equations
|
||||
PtrList<fvScalarMatrix> pEqnComps(phases.size());
|
||||
PtrList<volScalarField> dmdts(fluid.dmdts());
|
||||
forAll(phases, phasei)
|
||||
{
|
||||
phaseModel& phase = phases[phasei];
|
||||
const volScalarField& alpha = phase;
|
||||
volScalarField& rho = phase.thermoRef().rho();
|
||||
|
||||
pEqnComps.set(phasei, new fvScalarMatrix(p_rgh, dimVolume/dimTime));
|
||||
fvScalarMatrix& pEqnComp = pEqnComps[phasei];
|
||||
|
||||
// Density variation
|
||||
if (!phase.isochoric() || !phase.pure())
|
||||
{
|
||||
pEqnComp +=
|
||||
(
|
||||
fvc::ddt(alpha, rho) + fvc::div(phase.alphaRhoPhi())
|
||||
- fvc::Sp(fvc::ddt(alpha) + fvc::div(phase.alphaPhi()), rho)
|
||||
)/rho;
|
||||
}
|
||||
|
||||
// Mesh dilatation correction
|
||||
if (mesh.moving())
|
||||
{
|
||||
pEqnComp += fvc::div(mesh.phi())*alpha;
|
||||
}
|
||||
|
||||
// Compressibility
|
||||
if (!phase.incompressible())
|
||||
{
|
||||
if (pimple.transonic())
|
||||
{
|
||||
const surfaceScalarField phid
|
||||
(
|
||||
IOobject::groupName("phid", phase.name()),
|
||||
fvc::interpolate(phase.thermo().psi())*phase.phi()
|
||||
);
|
||||
|
||||
pEqnComp +=
|
||||
correction
|
||||
(
|
||||
(alpha/rho)*
|
||||
(
|
||||
phase.thermo().psi()*fvm::ddt(p_rgh)
|
||||
+ fvm::div(phid, p_rgh)
|
||||
- fvm::Sp(fvc::div(phid), p_rgh)
|
||||
)
|
||||
);
|
||||
|
||||
deleteDemandDrivenData
|
||||
(
|
||||
pEqnComps[phasei].faceFluxCorrectionPtr()
|
||||
);
|
||||
|
||||
pEqnComps[phasei].relax();
|
||||
}
|
||||
else
|
||||
{
|
||||
pEqnComp +=
|
||||
(alpha*phase.thermo().psi()/rho)
|
||||
*correction(fvm::ddt(p_rgh));
|
||||
}
|
||||
}
|
||||
|
||||
// Option sources
|
||||
if (fvOptions.appliesToField(rho.name()))
|
||||
{
|
||||
pEqnComp -= (fvOptions(alpha, rho) & rho)/rho;
|
||||
}
|
||||
|
||||
// Mass transfer
|
||||
if (dmdts.set(phasei))
|
||||
{
|
||||
pEqnComp -= dmdts[phasei]/rho;
|
||||
}
|
||||
}
|
||||
#include "pEqnComps.H"
|
||||
|
||||
// Cache p prior to solve for density update
|
||||
volScalarField p_rgh_0(p_rgh);
|
||||
|
||||
@ -60,6 +60,8 @@ PIMPLE
|
||||
|
||||
pRefCell 0;
|
||||
pRefValue 1e5;
|
||||
|
||||
frozenFlow yes;
|
||||
}
|
||||
|
||||
relaxationFactors
|
||||
|
||||
@ -60,6 +60,8 @@ PIMPLE
|
||||
|
||||
pRefCell 0;
|
||||
pRefValue 1e5;
|
||||
|
||||
frozenFlow yes;
|
||||
}
|
||||
|
||||
relaxationFactors
|
||||
|
||||
Reference in New Issue
Block a user