zhyang-dev/OpenFOAM-6
1
0
Fork 0
mirror of https://github.com/OpenFOAM/OpenFOAM-6.git synced 2025-12-08 06:57:46 +00:00
Code Issues Actions Packages Projects Releases Wiki Activity

interFoam family of solvers: Improved Crank-Nicolson implementation

Browse Source 
Fewer limiter iterations are now required to obtain sufficient boundedness and
restart is more consistent.
This commit is contained in:
Henry Weller
2017-03-06 17:34:10 +00:00
parent 5546fe444d
commit 13aacde1c7
9 changed files with 70 additions and 76 deletions
Download Patch File Download Diff File Expand all files Collapse all files

84
applications/solvers/multiphase/VoF/alphaEqn.H
View File

@ -2,46 +2,57 @@
word alphaScheme("div(phi,alpha)");
word alpharScheme("div(phirb,alpha)");
tmp<fv::ddtScheme<scalar>> ddtAlpha
(
fv::ddtScheme<scalar>::New
(
mesh,
mesh.ddtScheme("ddt(alpha)")
)
);
// Set the off-centering coefficient according to ddt scheme
scalar ocCoeff = 0;
if
(
isType<fv::EulerDdtScheme<scalar>>(ddtAlpha())
|| isType<fv::localEulerDdtScheme<scalar>>(ddtAlpha())
)
{
ocCoeff = 0;
}
else if (isType<fv::CrankNicolsonDdtScheme<scalar>>(ddtAlpha()))
{
if (nAlphaSubCycles > 1)
tmp<fv::ddtScheme<scalar>> tddtAlpha
(
fv::ddtScheme<scalar>::New
(
mesh,
mesh.ddtScheme("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
<< "Sub-cycling is not supported "
"with the CrankNicolson ddt scheme"
<< "Only Euler and CrankNicolson ddt schemes are supported"
<< exit(FatalError);
}
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);
// Standard face-flux compression coefficient
@ -136,8 +147,8 @@
(
fvc::flux
(
phi,
alpha1,
phiCN(),
cnCoeff*alpha1 + (1.0 - cnCoeff)*alpha1.oldTime(),
alphaScheme
)
+ fvc::flux
@ -148,13 +159,6 @@
)
);
// Calculate the Crank-Nicolson off-centred alpha flux
if (ocCoeff > 0)
{
talphaPhiUn =
cnCoeff*talphaPhiUn + (1.0 - cnCoeff)*alphaPhi.oldTime();
}
if (MULESCorr)
{
tmp<surfaceScalarField> talphaPhiCorr(talphaPhiUn() - alphaPhi);

20
applications/solvers/multiphase/VoF/createAlphaFluxes.H Normal file
View File

@ -0,0 +1,20 @@
IOobject alphaPhiHeader
(
"alphaPhi",
runTime.timeName(),
mesh,
IOobject::READ_IF_PRESENT,
IOobject::AUTO_WRITE
);
const bool alphaRestart = alphaPhiHeader.headerOk();
// MULES flux from previous time-step
surfaceScalarField alphaPhi
(
alphaPhiHeader,
phi*fvc::interpolate(alpha1)
);
// MULES Correction
tmp<surfaceScalarField> talphaPhiCorr0;

1
applications/solvers/multiphase/compressibleInterFoam/compressibleInterDyMFoam/compressibleInterDyMFoam.C
View File

@ -64,6 +64,7 @@ int main(int argc, char *argv[])
#include "initContinuityErrs.H"
#include "createControl.H"
#include "createFields.H"
#include "createAlphaFluxes.H"
#include "createFvOptions.H"
#include "createUf.H"
#include "createControls.H"

1
applications/solvers/multiphase/compressibleInterFoam/compressibleInterFoam.C
View File

@ -61,6 +61,7 @@ int main(int argc, char *argv[])
#include "createControl.H"
#include "createTimeControls.H"
#include "createFields.H"
#include "createAlphaFluxes.H"
#include "createFvOptions.H"
volScalarField& p = mixture.p();

17
applications/solvers/multiphase/compressibleInterFoam/createFields.H
View File

@ -101,21 +101,4 @@ autoPtr<compressible::turbulenceModel> turbulence
Info<< "Creating field kinetic energy K\n" << endl;
volScalarField K("K", 0.5*magSqr(U));
// MULES flux from previous time-step
surfaceScalarField alphaPhi
(
IOobject
(
"alphaPhi",
runTime.timeName(),
mesh,
IOobject::READ_IF_PRESENT,
IOobject::AUTO_WRITE
),
phi*fvc::interpolate(alpha1)
);
// MULES Correction
tmp<surfaceScalarField> talphaPhiCorr0;
#include "createMRF.H"

17
applications/solvers/multiphase/interFoam/createFields.H
View File

@ -121,21 +121,4 @@ if (p_rgh.needReference())
mesh.setFluxRequired(p_rgh.name());
mesh.setFluxRequired(alpha1.name());
// MULES flux from previous time-step
surfaceScalarField alphaPhi
(
IOobject
(
"alphaPhi",
runTime.timeName(),
mesh,
IOobject::READ_IF_PRESENT,
IOobject::AUTO_WRITE
),
phi*fvc::interpolate(alpha1)
);
// MULES Correction
tmp<surfaceScalarField> talphaPhiCorr0;
#include "createMRF.H"

1
applications/solvers/multiphase/interFoam/interDyMFoam/interDyMFoam.C
View File

@ -60,6 +60,7 @@ int main(int argc, char *argv[])
#include "createTimeControls.H"
#include "createDyMControls.H"
#include "createFields.H"
#include "createAlphaFluxes.H"
#include "createFvOptions.H"
volScalarField rAU

1
applications/solvers/multiphase/interFoam/interFoam.C
View File

@ -63,6 +63,7 @@ int main(int argc, char *argv[])
#include "createTimeControls.H"
#include "initContinuityErrs.H"
#include "createFields.H"
#include "createAlphaFluxes.H"
#include "createFvOptions.H"
#include "correctPhi.H"

4
src/finiteVolume/finiteVolume/ddtSchemes/CrankNicolsonDdtScheme/CrankNicolsonDdtScheme.C
View File

@ -198,7 +198,7 @@ scalar CrankNicolsonDdtScheme<Type>::coef_
const DDt0Field<GeoField>& ddt0
) const
{
if (mesh().time().timeIndex() - ddt0.startTimeIndex() > 0)
if (mesh().time().timeIndex() > ddt0.startTimeIndex())
{
return 1 + ocCoeff_;
}
@ -216,7 +216,7 @@ scalar CrankNicolsonDdtScheme<Type>::coef0_
const DDt0Field<GeoField>& ddt0
) const
{
if (mesh().time().timeIndex() - ddt0.startTimeIndex() > 1)
if (mesh().time().timeIndex() > ddt0.startTimeIndex() + 1)
{
return 1 + ocCoeff_;
}