Skip to content

Commit

Permalink
multiphaseEuler::MomentumTransferPhaseSystem:alphaDByAf: multiphase c…
Browse files Browse the repository at this point in the history 
…onsistent replacement for DByAfs

In order that the phase-fractions sum to 1 it is necessary that the same
diffusivity is used for ALL phases in the implicitPhasePressure option.  This is
guaranteed by the new alphaDByAf function which returns a single
surfaceScalarField diffusivity to be used when forming the Laplacian term in the
implicit phase-fraction diffusion correction equation in phaseSystemSolve.

The phase-pressure and turbulent dispersion interface terms are summed over all
phases and interfaces in alphaDByAf to form a single diffusivity.
  • Loading branch information
Henry Weller committed Dec 1, 2022
1 parent e77a170 commit e8078ca
Show file tree
Hide file tree
Showing 4 changed files with 108 additions and 179 deletions.
228 changes: 87 additions & 141 deletions ...Euler/phaseSystems/PhaseSystems/MomentumTransferPhaseSystem/MomentumTransferPhaseSystem.C
View file
Expand Up @@ -79,6 +79,31 @@ void Foam::MomentumTransferPhaseSystem<BasePhaseSystem>::addDmdtUfs
}


template<class BasePhaseSystem>
void Foam::MomentumTransferPhaseSystem<BasePhaseSystem>::addTmpField
(
tmp<surfaceScalarField>& result,
const tmp<surfaceScalarField>& field
) const
{
if (result.valid())
{
result.ref() += field;
}
else
{
if (field.isTmp())
{
result = field;
}
else
{
result = field().clone();
}
}
}


// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //

template<class BasePhaseSystem>
Expand Down Expand Up @@ -863,168 +888,89 @@ implicitPhasePressure() const


template<class BasePhaseSystem>
Foam::PtrList<Foam::surfaceScalarField>
Foam::MomentumTransferPhaseSystem<BasePhaseSystem>::DByAfs
Foam::tmp<Foam::surfaceScalarField>
Foam::MomentumTransferPhaseSystem<BasePhaseSystem>::alphaDByAf
(
const PtrList<volScalarField>& rAUs,
const PtrList<surfaceScalarField>& rAUfs
) const
{
PtrList<surfaceScalarField> DByAfs(this->phaseModels_.size());
tmp<surfaceScalarField> alphaDByAf;

if (rAUfs.size())
// Add the phase pressure
forAll(this->phaseModels_, phasei)
{
// Add the phase pressure
forAll(this->phaseModels_, phasei)
{
const phaseModel& phase = this->phaseModels_[phasei];

const surfaceScalarField pPrimeByAf
(
rAUfs[phasei]*fvc::interpolate(phase.pPrime())
);

addField(phase, "DByAf", pPrimeByAf, DByAfs);

forAll(this->phaseModels_, phasej)
{
if (phasej != phasei)
{
const phaseModel& phase2 = this->phaseModels_[phasej];

addField
(
phase2,
"DByAf",
fvc::interpolate
(
phase2
/max(1 - phase, phase2.residualAlpha())
)*pPrimeByAf,
DByAfs
);
}
}
}
const phaseModel& phase = this->phaseModels_[phasei];

// Add the turbulent dispersion
forAllConstIter
addTmpField
(
turbulentDispersionModelTable,
turbulentDispersionModels_,
turbulentDispersionModelIter
)
{
const phaseInterface& interface =
turbulentDispersionModelIter()->interface();
alphaDByAf,
fvc::interpolate(max(phase, scalar(0)))
*fvc::interpolate(max(1 - phase, scalar(0))) // Optional
*(
rAUfs.size()

const surfaceScalarField Df
(
fvc::interpolate(turbulentDispersionModelIter()->D())
);
// Face-momentum form
? rAUfs[phasei]*fvc::interpolate(phase.pPrime())

const surfaceScalarField alpha12f
(
fvc::interpolate(interface.phase1() + interface.phase2())
);

addField
(
interface.phase1(),
"DByAf",
rAUfs[interface.phase1().index()]*Df
/max(alpha12f, interface.phase1().residualAlpha()),
DByAfs
);
addField
(
interface.phase2(),
"DByAf",
rAUfs[interface.phase2().index()]*Df
/max(alpha12f, interface.phase2().residualAlpha()),
DByAfs
);
}
// Cell-momentum form
: fvc::interpolate(rAUs[phasei]*phase.pPrime())
)
);
}
else
{
// Add the phase pressure
forAll(this->phaseModels_, phasei)
{
const phaseModel& phase = this->phaseModels_[phasei];

const surfaceScalarField pPrimeByAf
(
fvc::interpolate(rAUs[phasei]*phase.pPrime())
);

addField(phase, "DByAf", pPrimeByAf, DByAfs);

forAll(this->phaseModels_, phasej)
{
if (phasej != phasei)
{
const phaseModel& phase2 = this->phaseModels_[phasej];

addField
(
phase2,
"DByAf",
fvc::interpolate
(
phase2
/max(1 - phase, phase2.residualAlpha())
)*pPrimeByAf,
DByAfs
);
}
}
}
// Add the turbulent dispersion
forAllConstIter
(
turbulentDispersionModelTable,
turbulentDispersionModels_,
turbulentDispersionModelIter
)
{
const phaseInterface& interface =
turbulentDispersionModelIter()->interface();

// Add the turbulent dispersion
forAllConstIter
const surfaceScalarField alpha1f
(
turbulentDispersionModelTable,
turbulentDispersionModels_,
turbulentDispersionModelIter
)
{
const phaseInterface& interface =
turbulentDispersionModelIter()->interface();

const volScalarField D(turbulentDispersionModelIter()->D());
fvc::interpolate(max(interface.phase1(), scalar(0)))
);

const volScalarField alpha12
(
interface.phase1() + interface.phase2()
);
const surfaceScalarField alpha2f
(
fvc::interpolate(max(interface.phase2(), scalar(0)))
);

addField
(
interface.phase1(),
"DByAf",
fvc::interpolate
addTmpField
(
alphaDByAf,
alpha1f*alpha2f
/max(alpha1f + alpha2f, interface.phase1().residualAlpha())
*(
rAUfs.size()

// Face-momentum form
? max
(
rAUs[interface.phase1().index()]*D
/max(alpha12, interface.phase1().residualAlpha())
),
DByAfs
);
addField
(
interface.phase2(),
"DByAf",
fvc::interpolate
rAUfs[interface.phase1().index()],
rAUfs[interface.phase2().index()]
)
*fvc::interpolate(turbulentDispersionModelIter()->D())

// Cell-momentum form
: fvc::interpolate
(
rAUs[interface.phase2().index()]*D
/max(alpha12, interface.phase2().residualAlpha())
),
DByAfs
);
}
max
(
rAUs[interface.phase1().index()],
rAUs[interface.phase2().index()]
)
*turbulentDispersionModelIter()->D()
)
)
);
}

return DByAfs;
return alphaDByAf;
}


Expand Down
8 changes: 7 additions & 1 deletion ...Euler/phaseSystems/PhaseSystems/MomentumTransferPhaseSystem/MomentumTransferPhaseSystem.H
View file
Expand Up @@ -175,6 +175,12 @@ protected:
phaseSystem::momentumTransferTable& eqns
);

void addTmpField
(
tmp<surfaceScalarField>& result,
const tmp<surfaceScalarField>& field
) const;


public:

Expand Down Expand Up @@ -277,7 +283,7 @@ public:

//- Return the phase diffusivity
// divided by the momentum central coefficient
virtual PtrList<surfaceScalarField> DByAfs
virtual tmp<surfaceScalarField> alphaDByAf
(
const PtrList<volScalarField>& rAUs,
const PtrList<surfaceScalarField>& rAUfs
Expand Down
2 changes: 1 addition & 1 deletion applications/solvers/modules/multiphaseEuler/phaseSystems/phaseSystem/phaseSystem.H
View file
Expand Up @@ -557,7 +557,7 @@ public:

//- Return the phase diffusivity
// divided by the momentum central coefficient
virtual PtrList<surfaceScalarField> DByAfs
virtual tmp<surfaceScalarField> alphaDByAf
(
const PtrList<volScalarField>& rAUs,
const PtrList<surfaceScalarField>& rAUfs
Expand Down
49 changes: 13 additions & 36 deletions applications/solvers/modules/multiphaseEuler/phaseSystems/phaseSystem/phaseSystemSolve.C
View file
Expand Up @@ -100,25 +100,6 @@ void Foam::phaseSystem::solve
phases()[phasei].correctBoundaryConditions();
}

PtrList<surfaceScalarField> alphaPhiDbyA0s(phases().size());
if (implicitPhasePressure() && (rAUs.size() || rAUfs.size()))
{
const PtrList<surfaceScalarField> DByAfs(this->DByAfs(rAUs, rAUfs));

forAll(solvePhases, solvePhasei)
{
const phaseModel& phase = solvePhases[solvePhasei];
const volScalarField& alpha = phase;

alphaPhiDbyA0s.set
(
phase.index(),
DByAfs[phase.index()]
*fvc::snGrad(alpha, "bounded")*mesh_.magSf()
);
}
}

// Calculate the void fraction
volScalarField alphaVoid
(
Expand Down Expand Up @@ -269,6 +250,12 @@ void Foam::phaseSystem::solve
// Create correction fluxes
PtrList<surfaceScalarField> alphaPhis(phases().size());

tmp<surfaceScalarField> alphaDByAf;
if (implicitPhasePressure() && (rAUs.size() || rAUfs.size()))
{
alphaDByAf = this->alphaDByAf(rAUs, rAUfs);
}

forAll(movingPhases(), movingPhasei)
{
const phaseModel& phase = movingPhases()[movingPhasei];
Expand Down Expand Up @@ -383,12 +370,11 @@ void Foam::phaseSystem::solve
}
}

if (alphaPhiDbyA0s.set(phase.index()))
if (alphaDByAf.valid())
{
alphaPhi +=
fvc::interpolate(max(alpha, scalar(0)))
*fvc::interpolate(max(1 - alpha, scalar(0)))
*alphaPhiDbyA0s[phase.index()];
alphaDByAf()
*fvc::snGrad(alpha, "bounded")*mesh_.magSf();
}

phase.correctInflowOutflow(alphaPhi);
Expand Down Expand Up @@ -465,29 +451,20 @@ void Foam::phaseSystem::solve
}
}

if (implicitPhasePressure() && (rAUs.size() || rAUfs.size()))
if (alphaDByAf.valid())
{
const PtrList<surfaceScalarField> DByAfs
(
this->DByAfs(rAUs, rAUfs)
);
// Update alphaDByAf due to changes in alpha
alphaDByAf = this->alphaDByAf(rAUs, rAUfs);

forAll(solvePhases, solvePhasei)
{
phaseModel& phase = solvePhases[solvePhasei];
volScalarField& alpha = phase;

const surfaceScalarField alphaDbyA
(
fvc::interpolate(max(alpha, scalar(0)))
*fvc::interpolate(max(1 - alpha, scalar(0)))
*DByAfs[phase.index()]
);

fvScalarMatrix alphaEqn
(
fvm::ddt(alpha) - fvc::ddt(alpha)
- fvm::laplacian(alphaDbyA, alpha, "bounded")
- fvm::laplacian(alphaDByAf(), alpha, "bounded")
);

alphaEqn.solve();
Expand Down

0 comments on commit e8078ca

Please sign in to comment.