06893a0bc6
Much of the VoF functionality, particularly relating to momentum solution, is independent of the number of phases and it is useful to hold this generic VoF data and functionality in an abstract base-class and derive twoPhaseVoFSolver and multiphaseVoFSolver from it, adding two-phase and multiphase functionality respectively.
402 lines
11 KiB
C++
402 lines
11 KiB
C++
/*---------------------------------------------------------------------------*\
|
|
========= |
|
|
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
|
|
\\ / O peration | Website: https://openfoam.org
|
|
\\ / A nd | Copyright (C) 2023 OpenFOAM Foundation
|
|
\\/ M anipulation |
|
|
-------------------------------------------------------------------------------
|
|
License
|
|
This file is part of OpenFOAM.
|
|
|
|
OpenFOAM is free software: you can redistribute it and/or modify it
|
|
under the terms of the GNU General Public License as published by
|
|
the Free Software Foundation, either version 3 of the License, or
|
|
(at your option) any later version.
|
|
|
|
OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
|
|
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
|
|
FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
|
|
for more details.
|
|
|
|
You should have received a copy of the GNU General Public License
|
|
along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
|
|
|
|
\*---------------------------------------------------------------------------*/
|
|
|
|
#include "twoPhaseVoFSolver.H"
|
|
#include "subCycle.H"
|
|
#include "interfaceCompression.H"
|
|
#include "CMULES.H"
|
|
#include "CrankNicolsonDdtScheme.H"
|
|
#include "fvcFlux.H"
|
|
#include "fvmSup.H"
|
|
|
|
// * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * * //
|
|
|
|
void Foam::solvers::twoPhaseVoFSolver::alphaSolve
|
|
(
|
|
const dictionary& alphaControls
|
|
)
|
|
{
|
|
const label nAlphaSubCycles(alphaControls.lookup<label>("nAlphaSubCycles"));
|
|
|
|
const label nAlphaCorr(alphaControls.lookup<label>("nAlphaCorr"));
|
|
|
|
const bool MULESCorr
|
|
(
|
|
alphaControls.lookupOrDefault<Switch>("MULESCorr", false)
|
|
);
|
|
|
|
// Apply the compression correction from the previous iteration
|
|
// Improves efficiency for steady-simulations but can only be applied
|
|
// once the alpha field is reasonably steady, i.e. fully developed
|
|
const bool alphaApplyPrevCorr
|
|
(
|
|
alphaControls.lookupOrDefault<Switch>("alphaApplyPrevCorr", false)
|
|
);
|
|
|
|
const word alphaScheme(mesh.schemes().div(divAlphaName)[1].wordToken());
|
|
|
|
ITstream compressionScheme
|
|
(
|
|
compressionSchemes.found(alphaScheme)
|
|
? mesh.schemes().div(divAlphaName)
|
|
: ITstream
|
|
(
|
|
divAlphaName,
|
|
tokenList
|
|
{
|
|
word("Gauss"),
|
|
word("interfaceCompression"),
|
|
alphaScheme,
|
|
alphaControls.lookup<scalar>("cAlpha")
|
|
}
|
|
)
|
|
);
|
|
|
|
|
|
// Set the off-centering coefficient according to ddt scheme
|
|
scalar ocCoeff = 0;
|
|
{
|
|
tmp<fv::ddtScheme<scalar>> tddtAlpha
|
|
(
|
|
fv::ddtScheme<scalar>::New
|
|
(
|
|
mesh,
|
|
mesh.schemes().ddt("ddt(alpha)")
|
|
)
|
|
);
|
|
const fv::ddtScheme<scalar>& ddtAlpha = tddtAlpha();
|
|
|
|
if
|
|
(
|
|
isType<fv::EulerDdtScheme<scalar>>(ddtAlpha)
|
|
|| isType<fv::localEulerDdtScheme<scalar>>(ddtAlpha)
|
|
)
|
|
{
|
|
ocCoeff = 0;
|
|
}
|
|
else if (isType<fv::CrankNicolsonDdtScheme<scalar>>(ddtAlpha))
|
|
{
|
|
if (nAlphaSubCycles > 1)
|
|
{
|
|
FatalErrorInFunction
|
|
<< "Sub-cycling is not supported "
|
|
"with the CrankNicolson ddt scheme"
|
|
<< exit(FatalError);
|
|
}
|
|
|
|
if
|
|
(
|
|
alphaRestart
|
|
|| mesh.time().timeIndex() > mesh.time().startTimeIndex() + 1
|
|
)
|
|
{
|
|
ocCoeff =
|
|
refCast<const fv::CrankNicolsonDdtScheme<scalar>>(ddtAlpha)
|
|
.ocCoeff();
|
|
}
|
|
}
|
|
else
|
|
{
|
|
FatalErrorInFunction
|
|
<< "Only Euler and CrankNicolson ddt schemes are supported"
|
|
<< exit(FatalError);
|
|
}
|
|
}
|
|
|
|
// Set the time blending factor, 1 for Euler
|
|
scalar cnCoeff = 1.0/(1.0 + ocCoeff);
|
|
|
|
tmp<surfaceScalarField> phiCN(phi);
|
|
|
|
// Calculate the Crank-Nicolson off-centred volumetric flux
|
|
if (ocCoeff > 0)
|
|
{
|
|
phiCN = surfaceScalarField::New
|
|
(
|
|
"phiCN",
|
|
cnCoeff*phi + (1.0 - cnCoeff)*phi.oldTime()
|
|
);
|
|
}
|
|
|
|
tmp<volScalarField> divU;
|
|
|
|
if (divergent())
|
|
{
|
|
divU =
|
|
(
|
|
mesh.moving()
|
|
? fvc::div(phiCN + mesh.phi())
|
|
: fvc::div(phiCN)
|
|
);
|
|
}
|
|
|
|
tmp<volScalarField::Internal> Su;
|
|
tmp<volScalarField::Internal> Sp;
|
|
|
|
alphaSuSp(Su, Sp);
|
|
|
|
if (MULESCorr)
|
|
{
|
|
fvScalarMatrix alpha1Eqn
|
|
(
|
|
(
|
|
LTS
|
|
? fv::localEulerDdtScheme<scalar>(mesh).fvmDdt(alpha1)
|
|
: fv::EulerDdtScheme<scalar>(mesh).fvmDdt(alpha1)
|
|
)
|
|
+ fv::gaussConvectionScheme<scalar>
|
|
(
|
|
mesh,
|
|
phiCN,
|
|
upwind<scalar>(mesh, phiCN)
|
|
).fvmDiv(phiCN, alpha1)
|
|
);
|
|
|
|
if (divU.valid())
|
|
{
|
|
alpha1Eqn -= Su() + fvm::Sp(Sp() + divU(), alpha1);
|
|
}
|
|
|
|
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;
|
|
|
|
tmp<surfaceScalarField> talphaPhi1UD(alpha1Eqn.flux());
|
|
alphaPhi1 = talphaPhi1UD();
|
|
|
|
if (alphaApplyPrevCorr && talphaPhi1Corr0.valid())
|
|
{
|
|
Info<< "Applying the previous iteration compression flux" << endl;
|
|
MULES::correct
|
|
(
|
|
geometricOneField(),
|
|
alpha1,
|
|
alphaPhi1,
|
|
talphaPhi1Corr0.ref(),
|
|
oneField(),
|
|
zeroField()
|
|
);
|
|
|
|
alphaPhi1 += talphaPhi1Corr0();
|
|
}
|
|
|
|
// Cache the upwind-flux
|
|
talphaPhi1Corr0 = talphaPhi1UD;
|
|
|
|
alpha2 = 1.0 - alpha1;
|
|
|
|
interface.correct();
|
|
}
|
|
|
|
for (int aCorr=0; aCorr<nAlphaCorr; aCorr++)
|
|
{
|
|
// Split operator
|
|
tmp<surfaceScalarField> talphaPhi1Un
|
|
(
|
|
fvc::flux
|
|
(
|
|
phiCN(),
|
|
(cnCoeff*alpha1 + (1.0 - cnCoeff)*alpha1.oldTime())(),
|
|
compressionScheme.rewind()
|
|
)
|
|
);
|
|
|
|
if (MULESCorr)
|
|
{
|
|
tmp<surfaceScalarField> talphaPhi1Corr(talphaPhi1Un() - alphaPhi1);
|
|
volScalarField alpha10("alpha10", alpha1);
|
|
|
|
if (divU.valid())
|
|
{
|
|
MULES::correct
|
|
(
|
|
geometricOneField(),
|
|
alpha1,
|
|
talphaPhi1Un(),
|
|
talphaPhi1Corr.ref(),
|
|
Sp(),
|
|
(-Sp()*alpha1)(),
|
|
oneField(),
|
|
zeroField()
|
|
);
|
|
}
|
|
else
|
|
{
|
|
MULES::correct
|
|
(
|
|
geometricOneField(),
|
|
alpha1,
|
|
talphaPhi1Un(),
|
|
talphaPhi1Corr.ref(),
|
|
oneField(),
|
|
zeroField()
|
|
);
|
|
}
|
|
|
|
// Under-relax the correction for all but the 1st corrector
|
|
if (aCorr == 0)
|
|
{
|
|
alphaPhi1 += talphaPhi1Corr();
|
|
}
|
|
else
|
|
{
|
|
alpha1 = 0.5*alpha1 + 0.5*alpha10;
|
|
alphaPhi1 += 0.5*talphaPhi1Corr();
|
|
}
|
|
}
|
|
else
|
|
{
|
|
alphaPhi1 = talphaPhi1Un;
|
|
|
|
if (divU.valid())
|
|
{
|
|
MULES::explicitSolve
|
|
(
|
|
geometricOneField(),
|
|
alpha1,
|
|
phiCN,
|
|
alphaPhi1,
|
|
Sp(),
|
|
(Su() + divU()*min(alpha1(), scalar(1)))(),
|
|
oneField(),
|
|
zeroField()
|
|
);
|
|
}
|
|
else
|
|
{
|
|
MULES::explicitSolve
|
|
(
|
|
geometricOneField(),
|
|
alpha1,
|
|
phiCN,
|
|
alphaPhi1,
|
|
oneField(),
|
|
zeroField()
|
|
);
|
|
}
|
|
}
|
|
|
|
alpha2 = 1.0 - alpha1;
|
|
|
|
// Correct only the mixture interface for the interface compression flux
|
|
interface.correct();
|
|
}
|
|
|
|
if (alphaApplyPrevCorr && MULESCorr)
|
|
{
|
|
talphaPhi1Corr0 = alphaPhi1 - talphaPhi1Corr0;
|
|
talphaPhi1Corr0.ref().rename("alphaPhi1Corr0");
|
|
}
|
|
else
|
|
{
|
|
talphaPhi1Corr0.clear();
|
|
}
|
|
|
|
if
|
|
(
|
|
word(mesh.schemes().ddt("ddt(rho,U)"))
|
|
!= fv::EulerDdtScheme<vector>::typeName
|
|
&& word(mesh.schemes().ddt("ddt(rho,U)"))
|
|
!= fv::localEulerDdtScheme<vector>::typeName
|
|
)
|
|
{
|
|
if (ocCoeff > 0)
|
|
{
|
|
// Calculate the end-of-time-step alpha flux
|
|
alphaPhi1 =
|
|
(alphaPhi1 - (1.0 - cnCoeff)*alphaPhi1.oldTime())/cnCoeff;
|
|
}
|
|
}
|
|
|
|
Info<< "Phase-1 volume fraction = "
|
|
<< alpha1.weightedAverage(mesh.Vsc()).value()
|
|
<< " Min(" << alpha1.name() << ") = " << min(alpha1).value()
|
|
<< " Max(" << alpha1.name() << ") = " << max(alpha1).value()
|
|
<< endl;
|
|
}
|
|
|
|
|
|
void Foam::solvers::twoPhaseVoFSolver::alphaPredictor()
|
|
{
|
|
const dictionary& alphaControls = mesh.solution().solverDict(alpha1.name());
|
|
|
|
const label nAlphaSubCycles(alphaControls.lookup<label>("nAlphaSubCycles"));
|
|
|
|
|
|
if (nAlphaSubCycles > 1)
|
|
{
|
|
dimensionedScalar totalDeltaT = runTime.deltaT();
|
|
tmp<volScalarField> trSubDeltaT;
|
|
|
|
if (LTS)
|
|
{
|
|
trSubDeltaT =
|
|
fv::localEulerDdt::localRSubDeltaT(mesh, nAlphaSubCycles);
|
|
}
|
|
|
|
// Create a temporary alphaPhi1 to accumulate the sub-cycled alphaPhi1
|
|
tmp<surfaceScalarField> talphaPhi1
|
|
(
|
|
surfaceScalarField::New
|
|
(
|
|
"alphaPhi1",
|
|
mesh,
|
|
dimensionedScalar(alphaPhi1.dimensions(), 0)
|
|
)
|
|
);
|
|
|
|
List<volScalarField*> alphaPtrs({&alpha1, &alpha2});
|
|
|
|
for
|
|
(
|
|
subCycle<volScalarField, subCycleFields> alphaSubCycle
|
|
(
|
|
alphaPtrs,
|
|
nAlphaSubCycles
|
|
);
|
|
!(++alphaSubCycle).end();
|
|
)
|
|
{
|
|
alphaSolve(alphaControls);
|
|
talphaPhi1.ref() += (runTime.deltaT()/totalDeltaT)*alphaPhi1;
|
|
}
|
|
|
|
alphaPhi1 = talphaPhi1();
|
|
}
|
|
else
|
|
{
|
|
alphaSolve(alphaControls);
|
|
}
|
|
|
|
mixture.correct();
|
|
}
|
|
|
|
|
|
// ************************************************************************* //
|