zhyang-dev/OpenFOAM-12
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

compressibleVoF,multiphaseEuler: Renamed compressibility dilatation dgdt to vDot

Browse Source 
Currently in compressibleVoF vDot contains only the compressibility dilatation
effect whereas in multiphaseEuler the effect of sources are also included but
this will be refactored shortly so that the handling of mass sources and
compressibility is consistent between VoF and Euler-Euler solvers.

The previously hard-coded 1e-4 division stabilisation used when linearising vDot
for bounded semi-implicit solution of the phase-fractions is now an optional
user-input with keyword vDotResidualAlpha, e.g. in multiphaseEuler:

solvers
{
    "alpha.*"
    {
        nAlphaCorr          1;
        nAlphaSubCycles     2;
        vDotResidualAlpha   1e-6;
    }
    .
    .
    .
This commit is contained in:
Henry Weller
2023-11-03 13:19:52 +00:00
parent 7ec1e2f1a5
commit 0ed84ff137
15 changed files with 77 additions and 52 deletions
Download Patch File Download Diff File Expand all files Collapse all files

30
applications/modules/compressibleMultiphaseVoF/alphaPredictor.C
View File

@ -135,7 +135,7 @@ void Foam::solvers::compressibleMultiphaseVoF::alphaSolve
mesh
),
mesh,
dimensionedScalar(alpha.dgdt().dimensions(), 0)
dimensionedScalar(alpha.vDot().dimensions(), 0)
);
volScalarField::Internal Su
@ -152,18 +152,18 @@ void Foam::solvers::compressibleMultiphaseVoF::alphaSolve
);
{
const scalarField& dgdt = alpha.dgdt();
const scalarField& vDot = alpha.vDot();
forAll(dgdt, celli)
forAll(vDot, celli)
{
if (dgdt[celli] < 0.0 && alpha[celli] > 0.0)
if (vDot[celli] < 0.0 && alpha[celli] > 0.0)
{
Sp[celli] += dgdt[celli]*alpha[celli];
Su[celli] -= dgdt[celli]*alpha[celli];
Sp[celli] += vDot[celli]*alpha[celli];
Su[celli] -= vDot[celli]*alpha[celli];
}
else if (dgdt[celli] > 0.0 && alpha[celli] < 1.0)
else if (vDot[celli] > 0.0 && alpha[celli] < 1.0)
{
Sp[celli] -= dgdt[celli]*(1.0 - alpha[celli]);
Sp[celli] -= vDot[celli]*(1.0 - alpha[celli]);
}
}
}
@ -175,18 +175,18 @@ void Foam::solvers::compressibleMultiphaseVoF::alphaSolve
if (&alpha2 == &alpha) continue;
const scalarField& dgdt2 = alpha2.dgdt();
const scalarField& vDot2 = alpha2.vDot();
forAll(dgdt2, celli)
forAll(vDot2, celli)
{
if (dgdt2[celli] > 0.0 && alpha2[celli] < 1.0)
if (vDot2[celli] > 0.0 && alpha2[celli] < 1.0)
{
Sp[celli] -= dgdt2[celli]*(1.0 - alpha2[celli]);
Su[celli] += dgdt2[celli]*alpha[celli];
Sp[celli] -= vDot2[celli]*(1.0 - alpha2[celli]);
Su[celli] += vDot2[celli]*alpha[celli];
}
else if (dgdt2[celli] < 0.0 && alpha2[celli] > 0.0)
else if (vDot2[celli] < 0.0 && alpha2[celli] > 0.0)
{
Sp[celli] += dgdt2[celli]*alpha2[celli];
Sp[celli] += vDot2[celli]*alpha2[celli];
}
}
}

4
applications/modules/compressibleMultiphaseVoF/compressibleMultiphaseVoFMixture/compressibleVoFphase/compressibleVoFphase.C
View File

@ -47,11 +47,11 @@ Foam::compressibleVoFphase::compressibleVoFphase
mesh,
dimensionedScalar(dimless, 0)
),
dgdt_
vDot_
(
IOobject
(
IOobject::groupName("dgdt", name),
IOobject::groupName("vDot", name),
mesh.time().name(),
mesh,
IOobject::READ_IF_PRESENT,

10
applications/modules/compressibleMultiphaseVoF/compressibleMultiphaseVoFMixture/compressibleVoFphase/compressibleVoFphase.H
View File

@ -67,7 +67,7 @@ class compressibleVoFphase
volScalarField Alpha_;
//- Phase compressibility contribution
volScalarField::Internal dgdt_;
volScalarField::Internal vDot_;
public:
@ -140,15 +140,15 @@ public:
}
//- Return const-access to phase divergence
const volScalarField::Internal& dgdt() const
const volScalarField::Internal& vDot() const
{
return dgdt_;
return vDot_;
}
//- Return access to phase divergence
volScalarField::Internal& dgdt()
volScalarField::Internal& vDot()
{
return dgdt_;
return vDot_;
}
void correct

2
applications/modules/compressibleMultiphaseVoF/pressureCorrector.C
View File

@ -140,7 +140,7 @@ void Foam::solvers::compressibleMultiphaseVoF::pressureCorrector()
{
compressibleVoFphase& phase = phases[phasei];
phase.dgdt() =
phase.vDot() =
pos0(phase)
*(p_rghEqnComps[phasei] & p_rgh)/phase.thermo().rho();
}

24
applications/modules/compressibleVoF/alphaSuSp.C
View File

@ -31,10 +31,16 @@ License
void Foam::solvers::compressibleVoF::alphaSuSp
(
tmp<volScalarField::Internal>& tSu,
tmp<volScalarField::Internal>& tSp
tmp<volScalarField::Internal>& tSp,
const dictionary& alphaControls
)
{
const dimensionedScalar Szero(dgdt.dimensions(), 0);
const scalar vDotResidualAlpha
(
alphaControls.lookupOrDefault("vDotResidualAlpha", 1e-4)
);
const dimensionedScalar Szero(vDot.dimensions(), 0);
tSp = volScalarField::Internal::New("Sp", mesh, Szero);
tSu = volScalarField::Internal::New("Su", mesh, Szero);
@ -66,16 +72,18 @@ void Foam::solvers::compressibleVoF::alphaSuSp
Sp += alpha1ByRho2*alphaRho2Sup.Sp();
}
forAll(dgdt, celli)
forAll(vDot, celli)
{
if (dgdt[celli] > 0.0)
if (vDot[celli] > 0.0)
{
Sp[celli] -= dgdt[celli]/max(1.0 - alpha1[celli], 1e-4);
Su[celli] += dgdt[celli]/max(1.0 - alpha1[celli], 1e-4);
Sp[celli] -=
vDot[celli]/max(1.0 - alpha1[celli], vDotResidualAlpha);
Su[celli] +=
vDot[celli]/max(1.0 - alpha1[celli], vDotResidualAlpha);
}
else if (dgdt[celli] < 0.0)
else if (vDot[celli] < 0.0)
{
Sp[celli] += dgdt[celli]/max(alpha1[celli], 1e-4);
Sp[celli] += vDot[celli]/max(alpha1[celli], vDotResidualAlpha);
}
}
}

4
applications/modules/compressibleVoF/compressibleVoF.C
View File

@ -58,11 +58,11 @@ Foam::solvers::compressibleVoF::compressibleVoF(fvMesh& mesh)
p(mixture_.p()),
dgdt
vDot
(
IOobject
(
"dgdt",
"vDot",
runTime.name(),
mesh,
IOobject::READ_IF_PRESENT,

5
applications/modules/compressibleVoF/compressibleVoF.H
View File

@ -93,7 +93,7 @@ protected:
volScalarField& p;
//- Compressibility source
volScalarField::Internal dgdt;
volScalarField::Internal vDot;
// Pressure reference
@ -171,7 +171,8 @@ protected:
virtual void alphaSuSp
(
tmp<volScalarField::Internal>& Su,
tmp<volScalarField::Internal>& Sp
tmp<volScalarField::Internal>& Sp,
const dictionary& alphaControls
);
//- Return the momentum equation stress term

2
applications/modules/compressibleVoF/pressureCorrector.C
View File

@ -185,7 +185,7 @@ void Foam::solvers::compressibleVoF::pressureCorrector()
if (pimple.finalNonOrthogonalIter())
{
dgdt =
vDot =
(
alpha1*(p_rghEqnComp2 & p_rgh)
- alpha2*(p_rghEqnComp1 & p_rgh)

3
applications/modules/incompressibleDriftFlux/alphaSuSp.C
View File

@ -49,7 +49,8 @@ Foam::solvers::incompressibleDriftFlux::alphaPhi
void Foam::solvers::incompressibleDriftFlux::alphaSuSp
(
tmp<volScalarField::Internal>& tSu,
tmp<volScalarField::Internal>& tSp
tmp<volScalarField::Internal>& tSp,
const dictionary& alphaControls
)
{
if (!divergent()) return;

3
applications/modules/incompressibleDriftFlux/incompressibleDriftFlux.H
View File

@ -159,7 +159,8 @@ protected:
virtual void alphaSuSp
(
tmp<volScalarField::Internal>& Su,
tmp<volScalarField::Internal>& Sp
tmp<volScalarField::Internal>& Sp,
const dictionary& alphaControls
);
//- Correct the interface properties following mesh-change

3
applications/modules/incompressibleVoF/alphaSuSp.C
View File

@ -31,7 +31,8 @@ License
void Foam::solvers::incompressibleVoF::alphaSuSp
(
tmp<volScalarField::Internal>& tSu,
tmp<volScalarField::Internal>& tSp
tmp<volScalarField::Internal>& tSp,
const dictionary& alphaControls
)
{
if (!divergent()) return;

3
applications/modules/incompressibleVoF/incompressibleVoF.H
View File

@ -138,7 +138,8 @@ protected:
virtual void alphaSuSp
(
tmp<volScalarField::Internal>& Su,
tmp<volScalarField::Internal>& Sp
tmp<volScalarField::Internal>& Sp,
const dictionary& alphaControls
);
//- Return the momentum equation stress term

31
applications/modules/multiphaseEuler/phaseSystem/phaseSystem/phaseSystemSolve.C
View File

@ -64,6 +64,11 @@ void Foam::phaseSystem::solve(const PtrList<volScalarField>& rAs)
alphaControls.lookupOrDefault<Switch>("meanFluxReference", false)
);
const scalar vDotResidualAlpha
(
alphaControls.lookupOrDefault("vDotResidualAlpha", 1e-4)
);
// Optional reference phase which is not solved for
// but obtained from the sum of the other phases
phaseModel* referencePhasePtr = nullptr;
@ -172,11 +177,11 @@ void Foam::phaseSystem::solve(const PtrList<volScalarField>& rAs)
if (dilatation)
{
// Construct the dilatation rate source term
volScalarField::Internal dgdt
volScalarField::Internal vDot
(
volScalarField::Internal::New
(
"dgdt",
"vDot",
mesh_,
dimensionedScalar(dimless/dimTime, 0)
)
@ -191,12 +196,12 @@ void Foam::phaseSystem::solve(const PtrList<volScalarField>& rAs)
{
if (!phase.stationary() && phase.divU().valid())
{
dgdt += alpha2()*phase.divU()()();
vDot += alpha2()*phase.divU()()();
}
if (!phase2.stationary() && phase2.divU().valid())
{
dgdt -= alpha()*phase2.divU()()();
vDot -= alpha()*phase2.divU()()();
}
}
}
@ -204,16 +209,22 @@ void Foam::phaseSystem::solve(const PtrList<volScalarField>& rAs)
volScalarField::Internal& Sp = Sps[phasei];
volScalarField::Internal& Su = Sus[phasei];
forAll(dgdt, celli)
forAll(vDot, celli)
{
if (dgdt[celli] > 0)
if (vDot[celli] > 0)
{
Sp[celli] -= dgdt[celli]/max(1 - alpha[celli], 1e-4);
Su[celli] += dgdt[celli]/max(1 - alpha[celli], 1e-4);
Sp[celli] -=
vDot[celli]
/max(1 - alpha[celli], vDotResidualAlpha);
Su[celli] +=
vDot[celli]
/max(1 - alpha[celli], vDotResidualAlpha);
}
else if (dgdt[celli] < 0)
else if (vDot[celli] < 0)
{
Sp[celli] += dgdt[celli]/max(alpha[celli], 1e-4);
Sp[celli] +=
vDot[celli]
/max(alpha[celli], vDotResidualAlpha);
}
}
}

2
applications/modules/twoPhaseSolver/alphaPredictor.C
View File

@ -171,7 +171,7 @@ void Foam::solvers::twoPhaseSolver::alphaSolve
tmp<volScalarField::Internal> Su;
tmp<volScalarField::Internal> Sp;
alphaSuSp(Su, Sp);
alphaSuSp(Su, Sp, alphaControls);
if (MULESCorr)
{

3
applications/modules/twoPhaseSolver/twoPhaseSolver.H
View File

@ -128,7 +128,8 @@ protected:
virtual void alphaSuSp
(
tmp<volScalarField::Internal>& Su,
tmp<volScalarField::Internal>& Sp
tmp<volScalarField::Internal>& Sp,
const dictionary& alphaControls
) = 0;
//- Correct the interface properties following mesh-change