Skip to content
GitLab
Commit 58312bf8 authored by sergio's avatar sergio Committed by Andrew Heather
Browse files

INT: Back-tracking multiphaseSystem

parent bb9225f3
Hide whitespace changes
Inline Side-by-side
src/phaseSystemModels/reactingEulerFoam/reactingMultiphaseEulerFoam/multiphaseSystem/multiphaseSystem.C
View file @ 58312bf8
......@@ -62,9 +62,9 @@ void Foam::multiphaseSystem::calcAlphas()
scalar level = 0.0;
alphas_ == 0.0;
for (const phaseModel& phase : phases())
forAll(phases(), i)
{
alphas_ += level * phase;
alphas_ += level*phases()[i];
level += 1.0;
}
}
......@@ -72,9 +72,9 @@ void Foam::multiphaseSystem::calcAlphas()
void Foam::multiphaseSystem::solveAlphas()
{
for (phaseModel& phase : phases())
forAll(phases(), phasei)
{
phase.correctBoundaryConditions();
phases()[phasei].correctBoundaryConditions();
}
// Calculate the void fraction
......@@ -89,9 +89,9 @@ void Foam::multiphaseSystem::solveAlphas()
mesh_,
dimensionedScalar("one", dimless, 1)
);
for (const phaseModel& phase : stationaryPhases())
forAll(stationaryPhases(), stationaryPhasei)
{
alphaVoid -= phase;
alphaVoid -= stationaryPhases()[stationaryPhasei];
}
// Generate face-alphas
......@@ -112,8 +112,10 @@ void Foam::multiphaseSystem::solveAlphas()
// Create correction fluxes
PtrList<surfaceScalarField> alphaPhiCorrs(phases().size());
for (const phaseModel& phase : stationaryPhases())
forAll(stationaryPhases(), stationaryPhasei)
{
phaseModel& phase = stationaryPhases()[stationaryPhasei];
alphaPhiCorrs.set
(
phase.index(),
......@@ -124,9 +126,10 @@ void Foam::multiphaseSystem::solveAlphas()
)
);
}
for (const phaseModel& phase : movingPhases())
forAll(movingPhases(), movingPhasei)
{
const volScalarField& alpha = phase;
phaseModel& phase = movingPhases()[movingPhasei];
volScalarField& alpha = phase;
alphaPhiCorrs.set
(
......@@ -140,18 +143,21 @@ void Foam::multiphaseSystem::solveAlphas()
surfaceScalarField& alphaPhiCorr = alphaPhiCorrs[phase.index()];
for (const phaseModel& phase2 : phases())
forAll(phases(), phasei)
{
const volScalarField& alpha2 = phase2;
phaseModel& phase2 = phases()[phasei];
volScalarField& alpha2 = phase2;
if (&phase2 == &phase) continue;
surfaceScalarField phir(phase.phi() - phase2.phi());
const auto cAlpha =
cAlphas_.cfind(phasePairKey(phase.name(), phase2.name()));
cAlphaTable::const_iterator cAlpha
(
cAlphas_.find(phasePairKey(phase.name(), phase2.name()))
);
if (cAlpha.found())
if (cAlpha != cAlphas_.end())
{
surfaceScalarField phic
(
......@@ -161,7 +167,7 @@ void Foam::multiphaseSystem::solveAlphas()
phir += min(cAlpha()*phic, max(phic))*nHatf(phase, phase2);
}
const word phirScheme
word phirScheme
(
"div(phir," + alpha2.name() + ',' + alpha.name() + ')'
);
......@@ -192,15 +198,16 @@ void Foam::multiphaseSystem::solveAlphas()
// Limit the flux sums, fixing those of the stationary phases
labelHashSet fixedAlphaPhiCorrs;
for (const phaseModel& phase : stationaryPhases())
forAll(stationaryPhases(), stationaryPhasei)
{
fixedAlphaPhiCorrs.insert(phase.index());
fixedAlphaPhiCorrs.insert(stationaryPhases()[stationaryPhasei].index());
}
MULES::limitSum(alphafs, alphaPhiCorrs, fixedAlphaPhiCorrs);
// Solve for the moving phase alphas
for (phaseModel& phase : movingPhases())
forAll(movingPhases(), movingPhasei)
{
phaseModel& phase = movingPhases()[movingPhasei];
volScalarField& alpha = phase;
surfaceScalarField& alphaPhi = alphaPhiCorrs[phase.index()];
......@@ -216,7 +223,7 @@ void Foam::multiphaseSystem::solveAlphas()
mesh_
),
mesh_,
dimensionedScalar(dimless/dimTime)
dimensionedScalar("zero", dimless/dimTime, 0)
);
volScalarField::Internal Su
......@@ -245,8 +252,9 @@ void Foam::multiphaseSystem::solveAlphas()
}
}
for (const phaseModel& phase2 : phases())
forAll(phases(), phasej)
{
const phaseModel& phase2 = phases()[phasej];
const volScalarField& alpha2 = phase2;
if (&phase2 == &phase) continue;
......@@ -288,8 +296,10 @@ void Foam::multiphaseSystem::solveAlphas()
}
// Report the phase fractions and the phase fraction sum
for (phaseModel& phase : phases())
forAll(phases(), phasei)
{
phaseModel& phase = phases()[phasei];
Info<< phase.name() << " fraction, min, max = "
<< phase.weightedAverage(mesh_.V()).value()
<< ' ' << min(phase).value()
......@@ -306,11 +316,11 @@ void Foam::multiphaseSystem::solveAlphas()
mesh_
),
mesh_,
dimensionedScalar(dimless)
dimensionedScalar("zero", dimless, 0)
);
for (const phaseModel& phase : movingPhases())
forAll(movingPhases(), movingPhasei)
{
sumAlphaMoving += phase;
sumAlphaMoving += movingPhases()[movingPhasei];
}
Info<< "Phase-sum volume fraction, min, max = "
......@@ -320,9 +330,9 @@ void Foam::multiphaseSystem::solveAlphas()
<< endl;
// Correct the sum of the phase fractions to avoid drift
for (phaseModel& phase : movingPhases())
forAll(movingPhases(), movingPhasei)
{
phase *= alphaVoid/sumAlphaMoving;
movingPhases()[movingPhasei] *= alphaVoid/sumAlphaMoving;
}
}
......@@ -397,11 +407,12 @@ void Foam::multiphaseSystem::correctContactAngle
/mesh_.magSf().boundaryField()[patchi]
);
const auto tp =
acap.thetaProps()
.cfind(phasePairKey(phase1.name(), phase2.name()));
alphaContactAngleFvPatchScalarField::thetaPropsTable::
const_iterator tp =
acap.thetaProps()
.find(phasePairKey(phase1.name(), phase2.name()));
if (!tp.found())
if (tp == acap.thetaProps().end())
{
FatalErrorInFunction
<< "Cannot find interface "
......@@ -411,7 +422,7 @@ void Foam::multiphaseSystem::correctContactAngle
<< exit(FatalError);
}
const bool matched = (tp.key().first() == phase1.name());
bool matched = (tp.key().first() == phase1.name());
scalar theta0 = convertToRad*tp().theta0(matched);
scalarField theta(boundary[patchi].size(), theta0);
......@@ -510,7 +521,7 @@ Foam::multiphaseSystem::multiphaseSystem
IOobject::AUTO_WRITE
),
mesh,
dimensionedScalar(dimless)
dimensionedScalar("zero", dimless, 0)
),
cAlphas_(lookup("interfaceCompression")),
......@@ -518,17 +529,23 @@ Foam::multiphaseSystem::multiphaseSystem
deltaN_
(
"deltaN",
1e-8/cbrt(average(mesh_.V()))
1e-8/pow(average(mesh_.V()), 1.0/3.0)
)
{
for (const phaseModel& phase : phases())
forAll(phases(), phasei)
{
const volScalarField& alpha = phase;
mesh_.setFluxRequired(alpha.name());
volScalarField& alphai = phases()[phasei];
mesh_.setFluxRequired(alphai.name());
}
}
// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
Foam::multiphaseSystem::~multiphaseSystem()
{}
// * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * * //
Foam::tmp<Foam::surfaceScalarField> Foam::multiphaseSystem::surfaceTension
......@@ -542,21 +559,23 @@ Foam::tmp<Foam::surfaceScalarField> Foam::multiphaseSystem::surfaceTension
(
"surfaceTension",
mesh_,
dimensionedScalar(dimensionSet(1, -2, -2, 0, 0))
dimensionedScalar("zero", dimensionSet(1, -2, -2, 0, 0), 0)
)
);
tSurfaceTension.ref().setOriented();
for (const phaseModel& phase2 : phases())
forAll(phases(), phasej)
{
const phaseModel& phase2 = phases()[phasej];
if (&phase2 != &phase1)
{
phasePairKey key12(phase1.name(), phase2.name());
const auto cAlpha = cAlphas_.cfind(key12);
cAlphaTable::const_iterator cAlpha(cAlphas_.find(key12));
if (cAlpha.found())
if (cAlpha != cAlphas_.end())
{
tSurfaceTension.ref() +=
fvc::interpolate(sigma(key12)*K(phase1, phase2))
......@@ -567,6 +586,7 @@ Foam::tmp<Foam::surfaceScalarField> Foam::multiphaseSystem::surfaceTension
}
}
}
tSurfaceTension->setOriented();
return tSurfaceTension;
......@@ -582,16 +602,16 @@ Foam::multiphaseSystem::nearInterface() const
(
"nearInterface",
mesh_,
dimensionedScalar(dimless)
dimensionedScalar("zero", dimless, 0)
)
);
for (const phaseModel& phase : phases())
forAll(phases(), phasei)
{
tnearInt.ref() = max
(
tnearInt(),
pos0(phase - 0.01)*pos0(0.99 - phase)
pos0(phases()[phasei] - 0.01)*pos0(0.99 - phases()[phasei])
);
}
......@@ -621,9 +641,9 @@ void Foam::multiphaseSystem::solve()
List<volScalarField*> alphaPtrs(phases().size());
PtrList<surfaceScalarField> alphaPhiSums(phases().size());
label phasei = 0;
for (phaseModel& phase : phases())
forAll(phases(), phasei)
{
phaseModel& phase = phases()[phasei];
volScalarField& alpha = phase;
alphaPtrs[phasei] = &alpha;
......@@ -640,11 +660,9 @@ void Foam::multiphaseSystem::solve()
mesh_
),
mesh_,
dimensionedScalar(dimensionSet(0, 3, -1, 0, 0))
dimensionedScalar("zero", dimensionSet(0, 3, -1, 0, 0), 0)
)
);
++phasei;
}
for
......@@ -659,18 +677,15 @@ void Foam::multiphaseSystem::solve()
{
solveAlphas();
phasei = 0;
for (const phaseModel& phase : phases())
forAll(phases(), phasei)
{
alphaPhiSums[phasei] += phase.alphaPhi();
++phasei;
alphaPhiSums[phasei] += phases()[phasei].alphaPhi();
}
}
phasei = 0;
for (phaseModel& phase : phases())
forAll(phases(), phasei)
{
phaseModel& phase = phases()[phasei];
if (phase.stationary()) continue;
phase.alphaPhiRef() = alphaPhiSums[phasei]/nAlphaSubCycles;
......@@ -681,14 +696,14 @@ void Foam::multiphaseSystem::solve()
solveAlphas();
}
for (phaseModel& phase : phases())
forAll(phases(), phasei)
{
phaseModel& phase = phases()[phasei];
if (phase.stationary()) continue;
phase.alphaRhoPhiRef() =
fvc::interpolate(phase.rho())*phase.alphaPhi();
// Ensure the phase-fractions are bounded
phase.clip(0, 1);
}
......
src/phaseSystemModels/reactingEulerFoam/reactingMultiphaseEulerFoam/multiphaseSystem/multiphaseSystem.H
View file @ 58312bf8
......@@ -153,7 +153,7 @@ public:
//- Destructor
virtual ~multiphaseSystem() = default;
virtual ~multiphaseSystem();
// Selectors
......
Supports Markdown
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