diff --git a/applications/solvers/multiphase/VoF/alphaEqn.H b/applications/solvers/multiphase/VoF/alphaEqn.H
index 570442a304914253f103a22e9740298073298879..abae00fb38541d4830e4a768e4deda155b240f90 100644
--- a/applications/solvers/multiphase/VoF/alphaEqn.H
+++ b/applications/solvers/multiphase/VoF/alphaEqn.H
@@ -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);
diff --git a/applications/solvers/multiphase/VoF/createAlphaFluxes.H b/applications/solvers/multiphase/VoF/createAlphaFluxes.H
new file mode 100644
index 0000000000000000000000000000000000000000..e75c2fe0b0d3fdef29e851caf9c8f96371ec26e0
--- /dev/null
+++ b/applications/solvers/multiphase/VoF/createAlphaFluxes.H
@@ -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;
diff --git a/applications/solvers/multiphase/compressibleInterFoam/compressibleInterDyMFoam/compressibleInterDyMFoam.C b/applications/solvers/multiphase/compressibleInterFoam/compressibleInterDyMFoam/compressibleInterDyMFoam.C
index f82665a49e44292e193317e6d7d4804ceafc6b8d..896af439af9fd80ea6e7789f5b4d0bef5df09742 100644
--- a/applications/solvers/multiphase/compressibleInterFoam/compressibleInterDyMFoam/compressibleInterDyMFoam.C
+++ b/applications/solvers/multiphase/compressibleInterFoam/compressibleInterDyMFoam/compressibleInterDyMFoam.C
@@ -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"
diff --git a/applications/solvers/multiphase/compressibleInterFoam/compressibleInterFoam.C b/applications/solvers/multiphase/compressibleInterFoam/compressibleInterFoam.C
index f6ec5d20cd2bca42b0bc60b6a1447dbab94fa08c..ce6446d0e1da6ff4e8d0ec4f5fe99c0b8653029e 100644
--- a/applications/solvers/multiphase/compressibleInterFoam/compressibleInterFoam.C
+++ b/applications/solvers/multiphase/compressibleInterFoam/compressibleInterFoam.C
@@ -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();
diff --git a/applications/solvers/multiphase/compressibleInterFoam/createFields.H b/applications/solvers/multiphase/compressibleInterFoam/createFields.H
index 9bac4d825922be2cae8a78f1f656001b49136c99..46168e6b87f7bfb44bba6c1313f37bf2e72b1787 100644
--- a/applications/solvers/multiphase/compressibleInterFoam/createFields.H
+++ b/applications/solvers/multiphase/compressibleInterFoam/createFields.H
@@ -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"
diff --git a/applications/solvers/multiphase/interFoam/createFields.H b/applications/solvers/multiphase/interFoam/createFields.H
index 4a82afbd294aa3a314ab5708e837f182f86c929d..dca977548e9809dee679b16a0b225b3a1b44382f 100644
--- a/applications/solvers/multiphase/interFoam/createFields.H
+++ b/applications/solvers/multiphase/interFoam/createFields.H
@@ -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"
diff --git a/applications/solvers/multiphase/interFoam/interDyMFoam/interDyMFoam.C b/applications/solvers/multiphase/interFoam/interDyMFoam/interDyMFoam.C
index df826479f0b10cbdec3b04362ba9bfc4ae6177cb..6a6da9d2d2c4746a43d92a47ff6cc26dce47a0c7 100644
--- a/applications/solvers/multiphase/interFoam/interDyMFoam/interDyMFoam.C
+++ b/applications/solvers/multiphase/interFoam/interDyMFoam/interDyMFoam.C
@@ -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
diff --git a/applications/solvers/multiphase/interFoam/interFoam.C b/applications/solvers/multiphase/interFoam/interFoam.C
index 93d70dea5e7520d17fe66ae1b5294a7b3c5dac31..7b755b6d192a71b5a30222b87fe76298bf5708f7 100644
--- a/applications/solvers/multiphase/interFoam/interFoam.C
+++ b/applications/solvers/multiphase/interFoam/interFoam.C
@@ -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"
diff --git a/src/finiteVolume/finiteVolume/ddtSchemes/CrankNicolsonDdtScheme/CrankNicolsonDdtScheme.C b/src/finiteVolume/finiteVolume/ddtSchemes/CrankNicolsonDdtScheme/CrankNicolsonDdtScheme.C
index 956898e4bd8c63bd3c9d57259569fa0814863ee8..0c7fe053b3bc6a0fbb5dd2e542ce3ced53f0ee53 100644
--- a/src/finiteVolume/finiteVolume/ddtSchemes/CrankNicolsonDdtScheme/CrankNicolsonDdtScheme.C
+++ b/src/finiteVolume/finiteVolume/ddtSchemes/CrankNicolsonDdtScheme/CrankNicolsonDdtScheme.C
@@ -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_;
}