From d13c9810e86f277923bbd4af02d0ea4c1880e032 Mon Sep 17 00:00:00 2001
From: Henry <Henry>
Date: Fri, 3 May 2013 15:46:29 +0100
Subject: [PATCH] VoF Solvers: Relocate the correction, sub-cycling and
compressions controls from PIMPLE to the alpha1 sub-dict MULES: Add support
for explicit correction interPhaseChangeFoam: Add option for explicit MULES
or as correction to an upwind solution Deprecate implicit form of MULES
---
.../compressibleInterFoam/alphaEqnsSubCycle.H | 4 +-
.../compressibleInterDyMFoam/readControls.H | 12 -
.../compressibleInterFoam.C | 4 +-
.../multiphase/compressibleInterFoam/pEqn.H | 6 +-
.../compressibleInterFoam/readControls.H | 13 -
.../readTwoPhaseEulerFoamControls.H | 4 +-
.../interFoam/LTSInterFoam/alphaEqnSubCycle.H | 3 +-
.../interFoam/LTSInterFoam/setrDeltaT.H | 5 +-
.../multiphase/interFoam/alphaEqnSubCycle.H | 3 +-
.../interMixingFoam/alphaEqnsSubCycle.H | 4 +-
.../interPhaseChangeFoam/alphaEqn.H | 124 ++---
.../interPhaseChangeFoam/alphaEqnSubCycle.H | 14 +-
.../interPhaseChangeFoam.C | 2 +-
.../multiphaseSystem/multiphaseSystem.C | 5 +-
.../multiphaseMixture/multiphaseMixture.C | 11 +-
.../twoLiquidMixingFoam/alphaEqnSubCycle.H | 3 +-
.../readTwoPhaseEulerFoamControls.H | 3 +-
.../fvMatrices/solvers/MULES/MULES.C | 19 +
.../fvMatrices/solvers/MULES/MULES.H | 61 +++
.../fvMatrices/solvers/MULES/MULESTemplates.C | 426 ++++++++++++++++++
.../interfaceProperties/interfaceProperties.C | 2 +-
.../cavitatingBullet/system/controlDict | 32 +-
.../cavitatingBullet/system/fvSchemes | 6 +-
.../cavitatingBullet/system/fvSolution | 26 +-
24 files changed, 635 insertions(+), 157 deletions(-)
delete mode 100644 applications/solvers/multiphase/compressibleInterFoam/readControls.H
diff --git a/applications/solvers/multiphase/compressibleInterFoam/alphaEqnsSubCycle.H b/applications/solvers/multiphase/compressibleInterFoam/alphaEqnsSubCycle.H
index be34ab52ed0..accc3197ad9 100644
--- a/applications/solvers/multiphase/compressibleInterFoam/alphaEqnsSubCycle.H
+++ b/applications/solvers/multiphase/compressibleInterFoam/alphaEqnsSubCycle.H
@@ -1,7 +1,5 @@
{
- label nAlphaCorr(readLabel(pimple.dict().lookup("nAlphaCorr")));
-
- label nAlphaSubCycles(readLabel(pimple.dict().lookup("nAlphaSubCycles")));
+ #include "alphaControls.H"
surfaceScalarField phic(mag(phi/mesh.magSf()));
phic = min(interface.cAlpha()*phic, max(phic));
diff --git a/applications/solvers/multiphase/compressibleInterFoam/compressibleInterDyMFoam/readControls.H b/applications/solvers/multiphase/compressibleInterFoam/compressibleInterDyMFoam/readControls.H
index 1a718434a81..cf509804604 100644
--- a/applications/solvers/multiphase/compressibleInterFoam/compressibleInterDyMFoam/readControls.H
+++ b/applications/solvers/multiphase/compressibleInterFoam/compressibleInterDyMFoam/readControls.H
@@ -1,17 +1,5 @@
#include "readTimeControls.H"
- label nAlphaCorr(readLabel(pimple.dict().lookup("nAlphaCorr")));
-
- label nAlphaSubCycles(readLabel(pimple.dict().lookup("nAlphaSubCycles")));
-
- if (nAlphaSubCycles > 1 && pimple.nCorrPIMPLE() != 1)
- {
- FatalErrorIn(args.executable())
- << "Sub-cycling alpha is only allowed for PISO, "
- "i.e. when the number of outer-correctors = 1"
- << exit(FatalError);
- }
-
bool correctPhi =
pimple.dict().lookupOrDefault<Switch>("correctPhi", true);
diff --git a/applications/solvers/multiphase/compressibleInterFoam/compressibleInterFoam.C b/applications/solvers/multiphase/compressibleInterFoam/compressibleInterFoam.C
index a42d9119804..ed5faa7c86a 100644
--- a/applications/solvers/multiphase/compressibleInterFoam/compressibleInterFoam.C
+++ b/applications/solvers/multiphase/compressibleInterFoam/compressibleInterFoam.C
@@ -56,7 +56,7 @@ int main(int argc, char *argv[])
pimpleControl pimple(mesh);
- #include "readControls.H"
+ #include "readTimeControls.H"
#include "initContinuityErrs.H"
#include "createFields.H"
#include "CourantNo.H"
@@ -68,7 +68,7 @@ int main(int argc, char *argv[])
while (runTime.run())
{
- #include "readControls.H"
+ #include "readTimeControls.H"
#include "CourantNo.H"
#include "setDeltaT.H"
diff --git a/applications/solvers/multiphase/compressibleInterFoam/pEqn.H b/applications/solvers/multiphase/compressibleInterFoam/pEqn.H
index 73babb08f04..13193788219 100644
--- a/applications/solvers/multiphase/compressibleInterFoam/pEqn.H
+++ b/applications/solvers/multiphase/compressibleInterFoam/pEqn.H
@@ -85,8 +85,8 @@
if (pimple.finalNonOrthogonalIter())
{
- //p = max(p_rgh + (alpha1*rho1 + alpha2*rho2)*gh, pMin);
- //p_rgh = p - (alpha1*rho1 + alpha2*rho2)*gh;
+ p = max(p_rgh + (alpha1*rho1 + alpha2*rho2)*gh, pMin);
+ p_rgh = p - (alpha1*rho1 + alpha2*rho2)*gh;
dgdt =
(
@@ -102,7 +102,7 @@
}
}
- p = max(p_rgh + (alpha1*rho1 + alpha2*rho2)*gh, pMin);
+ // p = max(p_rgh + (alpha1*rho1 + alpha2*rho2)*gh, pMin);
// Update densities from change in p_rgh
rho1 += psi1*(p_rgh - p_rgh_0);
diff --git a/applications/solvers/multiphase/compressibleInterFoam/readControls.H b/applications/solvers/multiphase/compressibleInterFoam/readControls.H
deleted file mode 100644
index 87a055d641f..00000000000
--- a/applications/solvers/multiphase/compressibleInterFoam/readControls.H
+++ /dev/null
@@ -1,13 +0,0 @@
- #include "readTimeControls.H"
-
- label nAlphaCorr(readLabel(pimple.dict().lookup("nAlphaCorr")));
-
- label nAlphaSubCycles(readLabel(pimple.dict().lookup("nAlphaSubCycles")));
-
- if (nAlphaSubCycles > 1 && pimple.nCorrPIMPLE() != 1)
- {
- FatalErrorIn(args.executable())
- << "Sub-cycling alpha is only allowed for PISO operation, "
- "i.e. when the number of outer-correctors = 1"
- << exit(FatalError);
- }
diff --git a/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/readTwoPhaseEulerFoamControls.H b/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/readTwoPhaseEulerFoamControls.H
index cb8b4efc3bb..9913595cf29 100644
--- a/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/readTwoPhaseEulerFoamControls.H
+++ b/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/readTwoPhaseEulerFoamControls.H
@@ -1,4 +1,2 @@
#include "readTimeControls.H"
-
- int nAlphaCorr(readInt(pimple.dict().lookup("nAlphaCorr")));
- int nAlphaSubCycles(readInt(pimple.dict().lookup("nAlphaSubCycles")));
+ #include "alphaControls.H"
diff --git a/applications/solvers/multiphase/interFoam/LTSInterFoam/alphaEqnSubCycle.H b/applications/solvers/multiphase/interFoam/LTSInterFoam/alphaEqnSubCycle.H
index 57c78027a47..428876fd229 100644
--- a/applications/solvers/multiphase/interFoam/LTSInterFoam/alphaEqnSubCycle.H
+++ b/applications/solvers/multiphase/interFoam/LTSInterFoam/alphaEqnSubCycle.H
@@ -1,5 +1,4 @@
-label nAlphaCorr(readLabel(pimple.dict().lookup("nAlphaCorr")));
-label nAlphaSubCycles(readLabel(pimple.dict().lookup("nAlphaSubCycles")));
+#include "alphaControls.H"
if (nAlphaSubCycles > 1)
{
diff --git a/applications/solvers/multiphase/interFoam/LTSInterFoam/setrDeltaT.H b/applications/solvers/multiphase/interFoam/LTSInterFoam/setrDeltaT.H
index 8525a8bf673..f419153506c 100644
--- a/applications/solvers/multiphase/interFoam/LTSInterFoam/setrDeltaT.H
+++ b/applications/solvers/multiphase/interFoam/LTSInterFoam/setrDeltaT.H
@@ -130,10 +130,7 @@
<< ", " << gMax(1/rDeltaT.internalField()) << endl;
}
- label nAlphaSubCycles
- (
- readLabel(pimpleDict.lookup("nAlphaSubCycles"))
- );
+ #include "alphaControls.H"
rSubDeltaT = rDeltaT*nAlphaSubCycles;
}
diff --git a/applications/solvers/multiphase/interFoam/alphaEqnSubCycle.H b/applications/solvers/multiphase/interFoam/alphaEqnSubCycle.H
index 57c78027a47..428876fd229 100644
--- a/applications/solvers/multiphase/interFoam/alphaEqnSubCycle.H
+++ b/applications/solvers/multiphase/interFoam/alphaEqnSubCycle.H
@@ -1,5 +1,4 @@
-label nAlphaCorr(readLabel(pimple.dict().lookup("nAlphaCorr")));
-label nAlphaSubCycles(readLabel(pimple.dict().lookup("nAlphaSubCycles")));
+#include "alphaControls.H"
if (nAlphaSubCycles > 1)
{
diff --git a/applications/solvers/multiphase/interFoam/interMixingFoam/alphaEqnsSubCycle.H b/applications/solvers/multiphase/interFoam/interMixingFoam/alphaEqnsSubCycle.H
index 97a09ce017d..9e3f4f6bd47 100644
--- a/applications/solvers/multiphase/interFoam/interMixingFoam/alphaEqnsSubCycle.H
+++ b/applications/solvers/multiphase/interFoam/interMixingFoam/alphaEqnsSubCycle.H
@@ -1,6 +1,4 @@
-const dictionary& pimpleDict = pimple.dict();
-label nAlphaCorr(readLabel(pimpleDict.lookup("nAlphaCorr")));
-label nAlphaSubCycles(readLabel(pimpleDict.lookup("nAlphaSubCycles")));
+#include "alphaControls.H"
if (nAlphaSubCycles > 1)
{
diff --git a/applications/solvers/multiphase/interPhaseChangeFoam/alphaEqn.H b/applications/solvers/multiphase/interPhaseChangeFoam/alphaEqn.H
index 75015542b55..9e5461614be 100644
--- a/applications/solvers/multiphase/interPhaseChangeFoam/alphaEqn.H
+++ b/applications/solvers/multiphase/interPhaseChangeFoam/alphaEqn.H
@@ -4,9 +4,41 @@
surfaceScalarField phir("phir", phic*interface.nHatf());
- for (int gCorr=0; gCorr<nAlphaCorr; gCorr++)
+ Pair<tmp<volScalarField> > vDotAlphal =
+ twoPhaseProperties->vDotAlphal();
+ const volScalarField& vDotcAlphal = vDotAlphal[0]();
+ const volScalarField& vDotvAlphal = vDotAlphal[1]();
+ const volScalarField vDotvmcAlphal(vDotvAlphal - vDotcAlphal);
+
+ tmp<surfaceScalarField> tphiAlpha;
+
+ if (MULESCorr)
{
- surfaceScalarField phiAlpha
+ fvScalarMatrix alpha1Eqn
+ (
+ fvm::ddt(alpha1)
+ + fvm::div(phi, alpha1, "UD") - fvm::Sp(divU, alpha1)
+ ==
+ fvm::Sp(vDotvmcAlphal, alpha1)
+ + vDotcAlphal
+ );
+
+ alpha1Eqn.solve();
+
+ Info<< "Phase-1 volume fraction = "
+ << alpha1.weightedAverage(mesh.Vsc()).value()
+ << " Min(alpha1) = " << min(alpha1).value()
+ << " Max(alpha1) = " << max(alpha1).value()
+ << endl;
+
+ tphiAlpha = alpha1Eqn.flux();
+ }
+
+ volScalarField alpha10("alpha10", alpha1);
+
+ for (int aCorr=0; aCorr<nAlphaCorr; aCorr++)
+ {
+ tmp<surfaceScalarField> tphiAlphaCorr
(
fvc::flux
(
@@ -22,66 +54,52 @@
)
);
- Pair<tmp<volScalarField> > vDotAlphal =
- twoPhaseProperties->vDotAlphal();
- const volScalarField& vDotcAlphal = vDotAlphal[0]();
- const volScalarField& vDotvAlphal = vDotAlphal[1]();
+ if (MULESCorr)
+ {
+ tphiAlphaCorr() -= tphiAlpha();
- volScalarField Sp
- (
- IOobject
+
+ volScalarField alpha100("alpha100", alpha10);
+ alpha10 = alpha1;
+
+ MULES::correct
(
- "Sp",
- runTime.timeName(),
- mesh
- ),
- vDotvAlphal - vDotcAlphal
- );
+ geometricOneField(),
+ alpha1,
+ tphiAlphaCorr(),
+ vDotvmcAlphal,
+ (
+ divU*(alpha10 - alpha100)
+ - vDotvmcAlphal*alpha10
+ )(),
+ 1,
+ 0
+ );
- volScalarField Su
- (
- IOobject
+ tphiAlpha() += tphiAlphaCorr();
+ }
+ else
+ {
+ MULES::explicitSolve
(
- "Su",
- runTime.timeName(),
- mesh
- ),
- // Divergence term is handled explicitly to be
- // consistent with the explicit transport solution
- divU*alpha1
- + vDotcAlphal
- );
+ geometricOneField(),
+ alpha1,
+ phi,
+ tphiAlphaCorr(),
+ vDotvmcAlphal,
+ (divU*alpha1 + vDotcAlphal)(),
+ 1,
+ 0
+ );
- //MULES::explicitSolve
- //(
- // geometricOneField(),
- // alpha1,
- // phi,
- // phiAlpha,
- // Sp,
- // Su,
- // 1,
- // 0
- //);
-
- MULES::implicitSolve
- (
- geometricOneField(),
- alpha1,
- phi,
- phiAlpha,
- Sp,
- Su,
- 1,
- 0
- );
+ tphiAlpha = tphiAlphaCorr;
+ }
alpha2 = 1.0 - alpha1;
- rhoPhi +=
- (runTime.deltaT()/totalDeltaT)
- *(phiAlpha*(rho1 - rho2) + phi*rho2);
}
+ rhoPhi = tphiAlpha()*(rho1 - rho2) + phi*rho2;
+
Info<< "Liquid phase volume fraction = "
<< alpha1.weightedAverage(mesh.V()).value()
<< " Min(alpha1) = " << min(alpha1).value()
diff --git a/applications/solvers/multiphase/interPhaseChangeFoam/alphaEqnSubCycle.H b/applications/solvers/multiphase/interPhaseChangeFoam/alphaEqnSubCycle.H
index 5e7fd27f752..a4338b907f8 100644
--- a/applications/solvers/multiphase/interPhaseChangeFoam/alphaEqnSubCycle.H
+++ b/applications/solvers/multiphase/interPhaseChangeFoam/alphaEqnSubCycle.H
@@ -11,21 +11,18 @@ surfaceScalarField rhoPhi
);
{
- const dictionary& pimpleDict = pimple.dict();
-
- label nAlphaCorr(readLabel(pimpleDict.lookup("nAlphaCorr")));
-
- label nAlphaSubCycles(readLabel(pimpleDict.lookup("nAlphaSubCycles")));
+ #include "alphaControls.H"
surfaceScalarField phic(mag(phi/mesh.magSf()));
phic = min(interface.cAlpha()*phic, max(phic));
volScalarField divU(fvc::div(phi));
- dimensionedScalar totalDeltaT = runTime.deltaT();
-
if (nAlphaSubCycles > 1)
{
+ dimensionedScalar totalDeltaT = runTime.deltaT();
+ surfaceScalarField rhoPhiSum("rhoPhiSum", rhoPhi);
+
for
(
subCycle<volScalarField> alphaSubCycle(alpha1, nAlphaSubCycles);
@@ -33,7 +30,10 @@ surfaceScalarField rhoPhi
)
{
#include "alphaEqn.H"
+ rhoPhiSum += (runTime.deltaT()/totalDeltaT)*rhoPhi;
}
+
+ rhoPhi = rhoPhiSum;
}
else
{
diff --git a/applications/solvers/multiphase/interPhaseChangeFoam/interPhaseChangeFoam.C b/applications/solvers/multiphase/interPhaseChangeFoam/interPhaseChangeFoam.C
index edfc34dfb83..bc6af055f0c 100644
--- a/applications/solvers/multiphase/interPhaseChangeFoam/interPhaseChangeFoam.C
+++ b/applications/solvers/multiphase/interPhaseChangeFoam/interPhaseChangeFoam.C
@@ -41,7 +41,7 @@ Description
\*---------------------------------------------------------------------------*/
#include "fvCFD.H"
-#include "IMULES.H"
+#include "MULES.H"
#include "subCycle.H"
#include "interfaceProperties.H"
#include "phaseChangeTwoPhaseMixture.H"
diff --git a/applications/solvers/multiphase/multiphaseEulerFoam/multiphaseSystem/multiphaseSystem.C b/applications/solvers/multiphase/multiphaseEulerFoam/multiphaseSystem/multiphaseSystem.C
index 0cde87042d9..a19b6a3d9b5 100644
--- a/applications/solvers/multiphase/multiphaseEulerFoam/multiphaseSystem/multiphaseSystem.C
+++ b/applications/solvers/multiphase/multiphaseEulerFoam/multiphaseSystem/multiphaseSystem.C
@@ -812,9 +812,8 @@ void Foam::multiphaseSystem::solve()
const Time& runTime = mesh_.time();
- const dictionary& pimpleDict = mesh_.solutionDict().subDict("PIMPLE");
-
- label nAlphaSubCycles(readLabel(pimpleDict.lookup("nAlphaSubCycles")));
+ const dictionary& alphaControls = mesh_.solverDict(phases_.first().name());
+ label nAlphaSubCycles(readLabel(alphaControls.lookup("nAlphaSubCycles")));
if (nAlphaSubCycles > 1)
{
diff --git a/applications/solvers/multiphase/multiphaseInterFoam/multiphaseMixture/multiphaseMixture.C b/applications/solvers/multiphase/multiphaseInterFoam/multiphaseMixture/multiphaseMixture.C
index b3c4ecb1115..76b4dab8dd8 100644
--- a/applications/solvers/multiphase/multiphaseInterFoam/multiphaseMixture/multiphaseMixture.C
+++ b/applications/solvers/multiphase/multiphaseInterFoam/multiphaseMixture/multiphaseMixture.C
@@ -249,15 +249,12 @@ void Foam::multiphaseMixture::solve()
const Time& runTime = mesh_.time();
- const dictionary& pimpleDict = mesh_.solutionDict().subDict("PIMPLE");
-
- label nAlphaSubCycles(readLabel(pimpleDict.lookup("nAlphaSubCycles")));
-
- scalar cAlpha(readScalar(pimpleDict.lookup("cAlpha")));
-
-
volScalarField& alpha = phases_.first();
+ const dictionary& alphaControls = mesh_.solverDict(alpha.name());
+ label nAlphaSubCycles(readLabel(alphaControls.lookup("nAlphaSubCycles")));
+ scalar cAlpha(readScalar(alphaControls.lookup("cAlpha")));
+
if (nAlphaSubCycles > 1)
{
surfaceScalarField rhoPhiSum(0.0*rhoPhi_);
diff --git a/applications/solvers/multiphase/twoLiquidMixingFoam/alphaEqnSubCycle.H b/applications/solvers/multiphase/twoLiquidMixingFoam/alphaEqnSubCycle.H
index 57c78027a47..428876fd229 100644
--- a/applications/solvers/multiphase/twoLiquidMixingFoam/alphaEqnSubCycle.H
+++ b/applications/solvers/multiphase/twoLiquidMixingFoam/alphaEqnSubCycle.H
@@ -1,5 +1,4 @@
-label nAlphaCorr(readLabel(pimple.dict().lookup("nAlphaCorr")));
-label nAlphaSubCycles(readLabel(pimple.dict().lookup("nAlphaSubCycles")));
+#include "alphaControls.H"
if (nAlphaSubCycles > 1)
{
diff --git a/applications/solvers/multiphase/twoPhaseEulerFoam/readTwoPhaseEulerFoamControls.H b/applications/solvers/multiphase/twoPhaseEulerFoam/readTwoPhaseEulerFoamControls.H
index 1abe3ef10d4..813fdd8f103 100644
--- a/applications/solvers/multiphase/twoPhaseEulerFoam/readTwoPhaseEulerFoamControls.H
+++ b/applications/solvers/multiphase/twoPhaseEulerFoam/readTwoPhaseEulerFoamControls.H
@@ -1,5 +1,4 @@
#include "readTimeControls.H"
+ #include "alphaControls.H"
- int nAlphaCorr(readInt(pimple.dict().lookup("nAlphaCorr")));
- int nAlphaSubCycles(readInt(pimple.dict().lookup("nAlphaSubCycles")));
Switch correctAlpha(pimple.dict().lookup("correctAlpha"));
diff --git a/src/finiteVolume/fvMatrices/solvers/MULES/MULES.C b/src/finiteVolume/fvMatrices/solvers/MULES/MULES.C
index 62ca855e785..ff5994e670c 100644
--- a/src/finiteVolume/fvMatrices/solvers/MULES/MULES.C
+++ b/src/finiteVolume/fvMatrices/solvers/MULES/MULES.C
@@ -48,6 +48,25 @@ void Foam::MULES::explicitSolve
}
+void Foam::MULES::correct
+(
+ volScalarField& psi,
+ surfaceScalarField& phiPsiCorr,
+ const scalar psiMax,
+ const scalar psiMin
+)
+{
+ correct
+ (
+ geometricOneField(),
+ psi,
+ phiPsiCorr,
+ zeroField(), zeroField(),
+ psiMax, psiMin
+ );
+}
+
+
void Foam::MULES::limitSum(UPtrList<scalarField>& phiPsiCorrs)
{
forAll(phiPsiCorrs[0], facei)
diff --git a/src/finiteVolume/fvMatrices/solvers/MULES/MULES.H b/src/finiteVolume/fvMatrices/solvers/MULES/MULES.H
index 2307c289e40..c1d824b9503 100644
--- a/src/finiteVolume/fvMatrices/solvers/MULES/MULES.H
+++ b/src/finiteVolume/fvMatrices/solvers/MULES/MULES.H
@@ -92,6 +92,38 @@ void explicitSolve
const scalar psiMin
);
+
+template<class RhoType, class SpType, class SuType>
+void correct
+(
+ const RhoType& rho,
+ volScalarField& psi,
+ const surfaceScalarField& phiCorr,
+ const SpType& Sp,
+ const SuType& Su
+);
+
+template<class RhoType, class SpType, class SuType>
+void correct
+(
+ const RhoType& rho,
+ volScalarField& psi,
+ surfaceScalarField& phiCorr,
+ const SpType& Sp,
+ const SuType& Su,
+ const scalar psiMax,
+ const scalar psiMin
+);
+
+void correct
+(
+ volScalarField& psi,
+ surfaceScalarField& phiCorr,
+ const scalar psiMax,
+ const scalar psiMin
+);
+
+
template<class RhoType, class SpType, class SuType>
void limiter
(
@@ -122,6 +154,35 @@ void limit
const bool returnCorr
);
+
+template<class RhoType, class SpType, class SuType>
+void limiterCorr
+(
+ scalarField& allLambda,
+ const RhoType& rho,
+ const volScalarField& psi,
+ const surfaceScalarField& phiCorr,
+ const SpType& Sp,
+ const SuType& Su,
+ const scalar psiMax,
+ const scalar psiMin,
+ const label nLimiterIter
+);
+
+template<class RhoType, class SpType, class SuType>
+void limitCorr
+(
+ const RhoType& rho,
+ const volScalarField& psi,
+ surfaceScalarField& phiCorr,
+ const SpType& Sp,
+ const SuType& Su,
+ const scalar psiMax,
+ const scalar psiMin,
+ const label nLimiterIter
+);
+
+
void limitSum(UPtrList<scalarField>& phiPsiCorrs);
template<class SurfaceScalarFieldList>
diff --git a/src/finiteVolume/fvMatrices/solvers/MULES/MULESTemplates.C b/src/finiteVolume/fvMatrices/solvers/MULES/MULESTemplates.C
index a3cb9f17e9d..e0718f923a9 100644
--- a/src/finiteVolume/fvMatrices/solvers/MULES/MULESTemplates.C
+++ b/src/finiteVolume/fvMatrices/solvers/MULES/MULESTemplates.C
@@ -96,6 +96,74 @@ void Foam::MULES::explicitSolve
}
+template<class RhoType, class SpType, class SuType>
+void Foam::MULES::correct
+(
+ const RhoType& rho,
+ volScalarField& psi,
+ const surfaceScalarField& phiCorr,
+ const SpType& Sp,
+ const SuType& Su
+)
+{
+ Info<< "MULES: Correcting " << psi.name() << endl;
+
+ const fvMesh& mesh = psi.mesh();
+
+ const scalar deltaT = mesh.time().deltaTValue();
+
+ scalarField psiIf(psi.size(), 0);
+ fvc::surfaceIntegrate(psiIf, phiCorr);
+
+ if (mesh.moving())
+ {
+ psi.internalField() =
+ (
+ rho.field()*psi.internalField()/deltaT
+ + Su.field()
+ - psiIf
+ )/(rho.field()/deltaT - Sp.field());
+ }
+ else
+ {
+ psi.internalField() =
+ (
+ rho.field()*psi.internalField()/deltaT
+ + Su.field()
+ - psiIf
+ )/(rho.field()/deltaT - Sp.field());
+ }
+
+ psi.correctBoundaryConditions();
+}
+
+
+template<class RhoType, class SpType, class SuType>
+void Foam::MULES::correct
+(
+ const RhoType& rho,
+ volScalarField& psi,
+ surfaceScalarField& phiCorr,
+ const SpType& Sp,
+ const SuType& Su,
+ const scalar psiMax,
+ const scalar psiMin
+)
+{
+ const fvMesh& mesh = psi.mesh();
+
+ const dictionary& MULEScontrols = mesh.solverDict(psi.name());
+
+ label nLimiterIter
+ (
+ readLabel(MULEScontrols.lookup("nLimiterIter"))
+ );
+
+ limitCorr(rho, psi, phiCorr, Sp, Su, psiMax, psiMin, nLimiterIter);
+ correct(rho, psi, phiCorr, Sp, Su);
+}
+
+
template<class RhoType, class SpType, class SuType>
void Foam::MULES::limiter
(
@@ -511,6 +579,364 @@ void Foam::MULES::limit
}
+template<class RhoType, class SpType, class SuType>
+void Foam::MULES::limiterCorr
+(
+ scalarField& allLambda,
+ const RhoType& rho,
+ const volScalarField& psi,
+ const surfaceScalarField& phiCorr,
+ const SpType& Sp,
+ const SuType& Su,
+ const scalar psiMax,
+ const scalar psiMin,
+ const label nLimiterIter
+)
+{
+ const scalarField& psiIf = psi;
+ const volScalarField::GeometricBoundaryField& psiBf = psi.boundaryField();
+
+ const fvMesh& mesh = psi.mesh();
+
+ const labelUList& owner = mesh.owner();
+ const labelUList& neighb = mesh.neighbour();
+ tmp<volScalarField::DimensionedInternalField> tVsc = mesh.Vsc();
+ const scalarField& V = tVsc();
+ const scalar deltaT = mesh.time().deltaTValue();
+
+ const scalarField& phiCorrIf = phiCorr;
+ const surfaceScalarField::GeometricBoundaryField& phiCorrBf =
+ phiCorr.boundaryField();
+
+ slicedSurfaceScalarField lambda
+ (
+ IOobject
+ (
+ "lambda",
+ mesh.time().timeName(),
+ mesh,
+ IOobject::NO_READ,
+ IOobject::NO_WRITE,
+ false
+ ),
+ mesh,
+ dimless,
+ allLambda,
+ false // Use slices for the couples
+ );
+
+ scalarField& lambdaIf = lambda;
+ surfaceScalarField::GeometricBoundaryField& lambdaBf =
+ lambda.boundaryField();
+
+ scalarField psiMaxn(psiIf.size(), psiMin);
+ scalarField psiMinn(psiIf.size(), psiMax);
+
+ scalarField sumPhip(psiIf.size(), VSMALL);
+ scalarField mSumPhim(psiIf.size(), VSMALL);
+
+ forAll(phiCorrIf, facei)
+ {
+ label own = owner[facei];
+ label nei = neighb[facei];
+
+ psiMaxn[own] = max(psiMaxn[own], psiIf[nei]);
+ psiMinn[own] = min(psiMinn[own], psiIf[nei]);
+
+ psiMaxn[nei] = max(psiMaxn[nei], psiIf[own]);
+ psiMinn[nei] = min(psiMinn[nei], psiIf[own]);
+
+ scalar phiCorrf = phiCorrIf[facei];
+
+ if (phiCorrf > 0.0)
+ {
+ sumPhip[own] += phiCorrf;
+ mSumPhim[nei] += phiCorrf;
+ }
+ else
+ {
+ mSumPhim[own] -= phiCorrf;
+ sumPhip[nei] -= phiCorrf;
+ }
+ }
+
+ forAll(phiCorrBf, patchi)
+ {
+ const fvPatchScalarField& psiPf = psiBf[patchi];
+ const scalarField& phiCorrPf = phiCorrBf[patchi];
+
+ const labelList& pFaceCells = mesh.boundary()[patchi].faceCells();
+
+ if (psiPf.coupled())
+ {
+ const scalarField psiPNf(psiPf.patchNeighbourField());
+
+ forAll(phiCorrPf, pFacei)
+ {
+ label pfCelli = pFaceCells[pFacei];
+
+ psiMaxn[pfCelli] = max(psiMaxn[pfCelli], psiPNf[pFacei]);
+ psiMinn[pfCelli] = min(psiMinn[pfCelli], psiPNf[pFacei]);
+ }
+ }
+ else
+ {
+ forAll(phiCorrPf, pFacei)
+ {
+ label pfCelli = pFaceCells[pFacei];
+
+ psiMaxn[pfCelli] = max(psiMaxn[pfCelli], psiPf[pFacei]);
+ psiMinn[pfCelli] = min(psiMinn[pfCelli], psiPf[pFacei]);
+ }
+ }
+
+ forAll(phiCorrPf, pFacei)
+ {
+ label pfCelli = pFaceCells[pFacei];
+
+ scalar phiCorrf = phiCorrPf[pFacei];
+
+ if (phiCorrf > 0.0)
+ {
+ sumPhip[pfCelli] += phiCorrf;
+ }
+ else
+ {
+ mSumPhim[pfCelli] -= phiCorrf;
+ }
+ }
+ }
+
+ psiMaxn = min(psiMaxn, psiMax);
+ psiMinn = max(psiMinn, psiMin);
+
+ // scalar smooth = 0.5;
+ // psiMaxn = min((1.0 - smooth)*psiIf + smooth*psiMaxn, psiMax);
+ // psiMinn = max((1.0 - smooth)*psiIf + smooth*psiMinn, psiMin);
+
+ psiMaxn =
+ V
+ *(
+ (rho.field()/deltaT - Sp.field())*psiMaxn
+ - Su.field()
+ - rho.field()*psi.internalField()/deltaT
+ );
+
+ psiMinn =
+ V
+ *(
+ Su.field()
+ - (rho.field()/deltaT - Sp.field())*psiMinn
+ + rho.field()*psi.internalField()/deltaT
+ );
+
+ scalarField sumlPhip(psiIf.size());
+ scalarField mSumlPhim(psiIf.size());
+
+ for (int j=0; j<nLimiterIter; j++)
+ {
+ sumlPhip = 0.0;
+ mSumlPhim = 0.0;
+
+ forAll(lambdaIf, facei)
+ {
+ label own = owner[facei];
+ label nei = neighb[facei];
+
+ scalar lambdaPhiCorrf = lambdaIf[facei]*phiCorrIf[facei];
+
+ if (lambdaPhiCorrf > 0.0)
+ {
+ sumlPhip[own] += lambdaPhiCorrf;
+ mSumlPhim[nei] += lambdaPhiCorrf;
+ }
+ else
+ {
+ mSumlPhim[own] -= lambdaPhiCorrf;
+ sumlPhip[nei] -= lambdaPhiCorrf;
+ }
+ }
+
+ forAll(lambdaBf, patchi)
+ {
+ scalarField& lambdaPf = lambdaBf[patchi];
+ const scalarField& phiCorrfPf = phiCorrBf[patchi];
+
+ const labelList& pFaceCells = mesh.boundary()[patchi].faceCells();
+
+ forAll(lambdaPf, pFacei)
+ {
+ label pfCelli = pFaceCells[pFacei];
+
+ scalar lambdaPhiCorrf = lambdaPf[pFacei]*phiCorrfPf[pFacei];
+
+ if (lambdaPhiCorrf > 0.0)
+ {
+ sumlPhip[pfCelli] += lambdaPhiCorrf;
+ }
+ else
+ {
+ mSumlPhim[pfCelli] -= lambdaPhiCorrf;
+ }
+ }
+ }
+
+ forAll(sumlPhip, celli)
+ {
+ sumlPhip[celli] =
+ max(min
+ (
+ (sumlPhip[celli] + psiMaxn[celli])/mSumPhim[celli],
+ 1.0), 0.0
+ );
+
+ mSumlPhim[celli] =
+ max(min
+ (
+ (mSumlPhim[celli] + psiMinn[celli])/sumPhip[celli],
+ 1.0), 0.0
+ );
+ }
+
+ const scalarField& lambdam = sumlPhip;
+ const scalarField& lambdap = mSumlPhim;
+
+ forAll(lambdaIf, facei)
+ {
+ if (phiCorrIf[facei] > 0.0)
+ {
+ lambdaIf[facei] = min
+ (
+ lambdaIf[facei],
+ min(lambdap[owner[facei]], lambdam[neighb[facei]])
+ );
+ }
+ else
+ {
+ lambdaIf[facei] = min
+ (
+ lambdaIf[facei],
+ min(lambdam[owner[facei]], lambdap[neighb[facei]])
+ );
+ }
+ }
+
+
+ forAll(lambdaBf, patchi)
+ {
+ fvsPatchScalarField& lambdaPf = lambdaBf[patchi];
+ const scalarField& phiCorrfPf = phiCorrBf[patchi];
+ const fvPatchScalarField& psiPf = psiBf[patchi];
+
+ if (isA<wedgeFvPatch>(mesh.boundary()[patchi]))
+ {
+ lambdaPf = 0;
+ }
+ else if (psiPf.coupled())
+ {
+ const labelList& pFaceCells =
+ mesh.boundary()[patchi].faceCells();
+
+ forAll(lambdaPf, pFacei)
+ {
+ label pfCelli = pFaceCells[pFacei];
+
+ if (phiCorrfPf[pFacei] > 0.0)
+ {
+ lambdaPf[pFacei] =
+ min(lambdaPf[pFacei], lambdap[pfCelli]);
+ }
+ else
+ {
+ lambdaPf[pFacei] =
+ min(lambdaPf[pFacei], lambdam[pfCelli]);
+ }
+ }
+ }
+ else
+ {
+ const labelList& pFaceCells =
+ mesh.boundary()[patchi].faceCells();
+ const scalarField& phiCorrPf = phiCorrBf[patchi];
+
+ forAll(lambdaPf, pFacei)
+ {
+ // Limit outlet faces only
+ if (phiCorrPf[pFacei] > 0)
+ {
+ label pfCelli = pFaceCells[pFacei];
+
+ if (phiCorrfPf[pFacei] > 0.0)
+ {
+ lambdaPf[pFacei] =
+ min(lambdaPf[pFacei], lambdap[pfCelli]);
+ }
+ else
+ {
+ lambdaPf[pFacei] =
+ min(lambdaPf[pFacei], lambdam[pfCelli]);
+ }
+ }
+ }
+ }
+ }
+
+ syncTools::syncFaceList(mesh, allLambda, minEqOp<scalar>());
+ }
+}
+
+
+template<class RhoType, class SpType, class SuType>
+void Foam::MULES::limitCorr
+(
+ const RhoType& rho,
+ const volScalarField& psi,
+ surfaceScalarField& phiCorr,
+ const SpType& Sp,
+ const SuType& Su,
+ const scalar psiMax,
+ const scalar psiMin,
+ const label nLimiterIter
+)
+{
+ const fvMesh& mesh = psi.mesh();
+
+ scalarField allLambda(mesh.nFaces(), 1.0);
+
+ slicedSurfaceScalarField lambda
+ (
+ IOobject
+ (
+ "lambda",
+ mesh.time().timeName(),
+ mesh,
+ IOobject::NO_READ,
+ IOobject::NO_WRITE,
+ false
+ ),
+ mesh,
+ dimless,
+ allLambda,
+ false // Use slices for the couples
+ );
+
+ limiterCorr
+ (
+ allLambda,
+ rho,
+ psi,
+ phiCorr,
+ Sp,
+ Su,
+ psiMax,
+ psiMin,
+ nLimiterIter
+ );
+
+ phiCorr *= lambda;
+}
+
+
template<class SurfaceScalarFieldList>
void Foam::MULES::limitSum(SurfaceScalarFieldList& phiPsiCorrs)
{
diff --git a/src/transportModels/interfaceProperties/interfaceProperties.C b/src/transportModels/interfaceProperties/interfaceProperties.C
index 63a1340ee92..ec9ace6ef94 100644
--- a/src/transportModels/interfaceProperties/interfaceProperties.C
+++ b/src/transportModels/interfaceProperties/interfaceProperties.C
@@ -159,7 +159,7 @@ Foam::interfaceProperties::interfaceProperties
(
readScalar
(
- alpha1.mesh().solutionDict().subDict("PIMPLE").lookup("cAlpha")
+ alpha1.mesh().solverDict(alpha1.name()).lookup("cAlpha")
)
),
sigma_(dict.lookup("sigma")),
diff --git a/tutorials/multiphase/interPhaseChangeFoam/cavitatingBullet/system/controlDict b/tutorials/multiphase/interPhaseChangeFoam/cavitatingBullet/system/controlDict
index f190eaf2ff8..d1a84670894 100644
--- a/tutorials/multiphase/interPhaseChangeFoam/cavitatingBullet/system/controlDict
+++ b/tutorials/multiphase/interPhaseChangeFoam/cavitatingBullet/system/controlDict
@@ -14,37 +14,37 @@ FoamFile
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-application interPhaseChangeFoam;
+application interPhaseChangeFoam;
-startFrom latestTime;
+startFrom latestTime;
-startTime 0;
+startTime 0;
-stopAt endTime;
+stopAt endTime;
-endTime 0.05;
+endTime 0.05;
-deltaT 1e-8;
+deltaT 1e-8;
-writeControl adjustableRunTime;
+writeControl adjustableRunTime;
-writeInterval 0.001;
+writeInterval 0.001;
-purgeWrite 0;
+purgeWrite 0;
-writeFormat ascii;
+writeFormat ascii;
-writePrecision 6;
+writePrecision 6;
-writeCompression uncompressed;
+writeCompression uncompressed;
-timeFormat general;
+timeFormat general;
-runTimeModifiable yes;
+runTimeModifiable yes;
-adjustTimeStep on;
+adjustTimeStep on;
-maxCo 1.0;
+maxCo 5;
// ************************************************************************* //
diff --git a/tutorials/multiphase/interPhaseChangeFoam/cavitatingBullet/system/fvSchemes b/tutorials/multiphase/interPhaseChangeFoam/cavitatingBullet/system/fvSchemes
index 9143f603f4b..8ddd1fe903d 100644
--- a/tutorials/multiphase/interPhaseChangeFoam/cavitatingBullet/system/fvSchemes
+++ b/tutorials/multiphase/interPhaseChangeFoam/cavitatingBullet/system/fvSchemes
@@ -32,6 +32,7 @@ divSchemes
div(phi,k) Gauss linearUpwind grad(k);
div(phi,alpha) Gauss vanLeer;
div(phirb,alpha) Gauss interfaceCompression;
+ UD Gauss upwind;
div((muEff*dev(T(grad(U))))) Gauss linear;
}
@@ -47,10 +48,7 @@ laplacianSchemes
snGradSchemes
{
- default none;
- snGrad(pd) limited corrected 0.5;
- snGrad(rho) limited corrected 0.5;
- snGrad(alpha1) limited corrected 0.5;
+ default limited corrected 0.5;
}
fluxRequired
diff --git a/tutorials/multiphase/interPhaseChangeFoam/cavitatingBullet/system/fvSolution b/tutorials/multiphase/interPhaseChangeFoam/cavitatingBullet/system/fvSolution
index bfbd5df2715..a1f8f8ffd5e 100644
--- a/tutorials/multiphase/interPhaseChangeFoam/cavitatingBullet/system/fvSolution
+++ b/tutorials/multiphase/interPhaseChangeFoam/cavitatingBullet/system/fvSolution
@@ -18,22 +18,24 @@ solvers
{
alpha1
{
- maxUnboundedness 1e-5;
- CoCoeff 2;
- maxIter 5;
- nLimiterIter 2;
+ cAlpha 0;
+ nAlphaCorr 2;
+ nAlphaSubCycles 1;
- solver PBiCG;
- preconditioner DILU;
- tolerance 1e-12;
- relTol 0.1;
+ MULESCorr yes;
+ nLimiterIter 5;
+
+ solver PBiCG;
+ preconditioner DILU;
+ tolerance 1e-8;
+ relTol 0;
};
- U
+ "U.*"
{
solver PBiCG;
preconditioner DILU;
- tolerance 1e-06;
+ tolerance 1e-6;
relTol 0;
};
@@ -96,10 +98,6 @@ PIMPLE
nOuterCorrectors 1;
nCorrectors 3;
nNonOrthogonalCorrectors 0;
-
- cAlpha 0;
- nAlphaCorr 1;
- nAlphaSubCycles 1;
}
--
GitLab