Files
openfoam/applications/solvers/multiphase/interPhaseChangeFoam/alphaEqn.H
Sergio Ferraris 8170f2ad92 INT: Org integration of VOF, Euler phase solvers and models.
Integration of VOF MULES new interfaces. Update of VOF solvers and all instances
of MULES in the code.
Integration of reactingTwoPhaseEuler and reactingMultiphaseEuler solvers and sub-models
Updating reactingEuler tutorials accordingly (most of them tested)

New eRefConst thermo used in tutorials. Some modifications at thermo specie level
affecting mostly eThermo. hThermo mostly unaffected

New chtMultiRegionTwoPhaseEulerFoam solver for quenching and tutorial.

Phases sub-models for reactingTwoPhaseEuler and reactingMultiphaseEuler were moved
to src/phaseSystemModels/reactingEulerFoam in order to be used by BC for
chtMultiRegionTwoPhaseEulerFoam.

Update of interCondensatingEvaporatingFoam solver.
2019-06-07 09:38:35 +01:00

124 lines
3.2 KiB
C

{
word alphaScheme("div(phi,alpha)");
word alpharScheme("div(phirb,alpha)");
surfaceScalarField phir("phir", phic*interface.nHatf());
Pair<tmp<volScalarField>> vDotAlphal =
mixture->vDotAlphal();
const volScalarField& vDotcAlphal = vDotAlphal[0]();
const volScalarField& vDotvAlphal = vDotAlphal[1]();
const volScalarField vDotvmcAlphal(vDotvAlphal - vDotcAlphal);
tmp<surfaceScalarField> talphaPhi;
if (MULESCorr)
{
fvScalarMatrix alpha1Eqn
(
fv::EulerDdtScheme<scalar>(mesh).fvmDdt(alpha1)
+ fv::gaussConvectionScheme<scalar>
(
mesh,
phi,
upwind<scalar>(mesh, phi)
).fvmDiv(phi, alpha1)
- fvm::Sp(divU, alpha1)
==
fvm::Sp(vDotvmcAlphal, alpha1)
+ vDotcAlphal
);
alpha1Eqn.solve();
Info<< "Phase-1 volume fraction = "
<< alpha1.weightedAverage(mesh.Vsc()).value()
<< " Min(" << alpha1.name() << ") = " << min(alpha1).value()
<< " Max(" << alpha1.name() << ") = " << max(alpha1).value()
<< endl;
talphaPhi = alpha1Eqn.flux();
}
volScalarField alpha10("alpha10", alpha1);
for (int aCorr=0; aCorr<nAlphaCorr; aCorr++)
{
tmp<surfaceScalarField> talphaPhiCorr
(
fvc::flux
(
phi,
alpha1,
alphaScheme
)
+ fvc::flux
(
-fvc::flux(-phir, alpha2, alpharScheme),
alpha1,
alpharScheme
)
);
if (MULESCorr)
{
talphaPhiCorr.ref() -= talphaPhi();
volScalarField alpha100("alpha100", alpha10);
alpha10 = alpha1;
MULES::correct
(
geometricOneField(),
alpha1,
talphaPhi(),
talphaPhiCorr.ref(),
vDotvmcAlphal,
(
divU*(alpha10 - alpha100)
- vDotvmcAlphal*alpha10
)(),
oneField(),
zeroField()
);
// Under-relax the correction for all but the 1st corrector
if (aCorr == 0)
{
talphaPhi.ref() += talphaPhiCorr();
}
else
{
alpha1 = 0.5*alpha1 + 0.5*alpha10;
talphaPhi.ref() += 0.5*talphaPhiCorr();
}
}
else
{
MULES::explicitSolve
(
geometricOneField(),
alpha1,
phi,
talphaPhiCorr.ref(),
vDotvmcAlphal,
(divU*alpha1 + vDotcAlphal)(),
oneField(),
zeroField()
);
talphaPhi = talphaPhiCorr;
}
alpha2 = 1.0 - alpha1;
}
rhoPhi = talphaPhi()*(rho1 - rho2) + phi*rho2;
Info<< "Liquid phase volume fraction = "
<< alpha1.weightedAverage(mesh.V()).value()
<< " Min(" << alpha1.name() << ") = " << min(alpha1).value()
<< " Max(" << alpha1.name() << ") = " << max(alpha1).value()
<< endl;
}