From 1c2093c8b3db054c7add334fc7f2098829c46890 Mon Sep 17 00:00:00 2001
From: Henry Weller <http://cfd.direct>
Date: Tue, 17 Jan 2017 22:43:47 +0000
Subject: [PATCH] Multi-phase solvers: Improved handling of inflow/outflow BCs
in MULES
Avoids slight phase-fraction unboundedness at entertainment BCs and improved
robustness.
Additionally the phase-fractions in the multi-phase (rather than two-phase)
solvers are adjusted to avoid the slow growth of inconsistency ("drift") caused
by solving for all of the phase-fractions rather than deriving one from the
others.
---
.../multiphaseSystem/multiphaseSystem.C | 74 +++---
.../multiphaseSystem/phaseModel/phaseModel.C | 20 +-
.../multiphaseSystem/phaseModel/phaseModel.H | 5 +-
.../multiphaseMixture/multiphaseMixture.C | 10 +-
.../phaseModel/phaseModel/phaseModel.C | 20 +-
.../phaseModel/phaseModel/phaseModel.H | 5 +-
.../multiphaseSystem/multiphaseSystem.C | 50 ++--
.../twoPhaseSystem/twoPhaseSystem.C | 61 ++---
.../twoPhaseSystem/phaseModel/phaseModel.C | 18 +-
.../twoPhaseSystem/phaseModel/phaseModel.H | 6 +-
.../twoPhaseSystem/twoPhaseSystem.C | 13 +-
src/finiteVolume/Make/files | 1 -
.../solvers/MULES/CMULESTemplates.C | 50 ++--
.../fvMatrices/solvers/MULES/IMULES.C | 51 ----
.../fvMatrices/solvers/MULES/IMULES.H | 94 -------
.../solvers/MULES/IMULESTemplates.C | 236 ------------------
.../fvMatrices/solvers/MULES/MULESTemplates.C | 92 +++----
.../laminar/damBreak4phase/system/fvSolution | 2 +-
.../RAS/testTubeMixer/system/fvSolution | 2 +-
.../LES/nozzleFlow2D/system/fvSolution | 2 +-
.../damBreak/damBreak/system/fvSolution | 2 +-
.../laminar/damBreak4phase/system/fvSolution | 2 +-
.../damBreak4phaseFine/system/fvSolution | 2 +-
.../laminar/mixerVessel2D/system/fvSolution | 2 +-
24 files changed, 214 insertions(+), 606 deletions(-)
delete mode 100644 src/finiteVolume/fvMatrices/solvers/MULES/IMULES.C
delete mode 100644 src/finiteVolume/fvMatrices/solvers/MULES/IMULES.H
delete mode 100644 src/finiteVolume/fvMatrices/solvers/MULES/IMULESTemplates.C
diff --git a/applications/solvers/multiphase/multiphaseEulerFoam/multiphaseSystem/multiphaseSystem.C b/applications/solvers/multiphase/multiphaseEulerFoam/multiphaseSystem/multiphaseSystem.C
index be4def509ec..798e9c2eb6c 100644
--- a/applications/solvers/multiphase/multiphaseEulerFoam/multiphaseSystem/multiphaseSystem.C
+++ b/applications/solvers/multiphase/multiphaseEulerFoam/multiphaseSystem/multiphaseSystem.C
@@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
- \\ / A nd | Copyright (C) 2011-2016 OpenFOAM Foundation
+ \\ / A nd | Copyright (C) 2011-2017 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@@ -64,11 +64,8 @@ void Foam::multiphaseSystem::solveAlphas()
forAllIter(PtrDictionary<phaseModel>, phases_, iter)
{
- phaseModel& phase1 = iter();
- volScalarField& alpha1 = phase1;
-
- phase1.alphaPhi() =
- dimensionedScalar("0", dimensionSet(0, 3, -1, 0, 0), 0);
+ phaseModel& phase = iter();
+ volScalarField& alpha1 = phase;
alphaPhiCorrs.set
(
@@ -79,7 +76,7 @@ void Foam::multiphaseSystem::solveAlphas()
fvc::flux
(
phi_,
- phase1,
+ phase,
"div(phi," + alpha1.name() + ')'
)
)
@@ -92,13 +89,13 @@ void Foam::multiphaseSystem::solveAlphas()
phaseModel& phase2 = iter2();
volScalarField& alpha2 = phase2;
- if (&phase2 == &phase1) continue;
+ if (&phase2 == &phase) continue;
- surfaceScalarField phir(phase1.phi() - phase2.phi());
+ surfaceScalarField phir(phase.phi() - phase2.phi());
scalarCoeffSymmTable::const_iterator cAlpha
(
- cAlphas_.find(interfacePair(phase1, phase2))
+ cAlphas_.find(interfacePair(phase, phase2))
);
if (cAlpha != cAlphas_.end())
@@ -108,7 +105,7 @@ void Foam::multiphaseSystem::solveAlphas()
(mag(phi_) + mag(phir))/mesh_.magSf()
);
- phir += min(cAlpha()*phic, max(phic))*nHatf(phase1, phase2);
+ phir += min(cAlpha()*phic, max(phic))*nHatf(phase, phase2);
}
word phirScheme
@@ -119,39 +116,18 @@ void Foam::multiphaseSystem::solveAlphas()
alphaPhiCorr += fvc::flux
(
-fvc::flux(-phir, phase2, phirScheme),
- phase1,
+ phase,
phirScheme
);
}
- surfaceScalarField::Boundary& alphaPhiCorrBf =
- alphaPhiCorr.boundaryFieldRef();
-
- // Ensure that the flux at inflow BCs is preserved
- forAll(alphaPhiCorrBf, patchi)
- {
- fvsPatchScalarField& alphaPhiCorrp = alphaPhiCorrBf[patchi];
-
- if (!alphaPhiCorrp.coupled())
- {
- const scalarField& phi1p = phase1.phi().boundaryField()[patchi];
- const scalarField& alpha1p = alpha1.boundaryField()[patchi];
-
- forAll(alphaPhiCorrp, facei)
- {
- if (phi1p[facei] < 0)
- {
- alphaPhiCorrp[facei] = alpha1p[facei]*phi1p[facei];
- }
- }
- }
- }
+ phase.correctInflowOutflow(alphaPhiCorr);
MULES::limit
(
1.0/mesh_.time().deltaT().value(),
geometricOneField(),
- phase1,
+ phase,
phi_,
alphaPhiCorr,
zeroField(),
@@ -182,29 +158,30 @@ void Foam::multiphaseSystem::solveAlphas()
forAllIter(PtrDictionary<phaseModel>, phases_, iter)
{
- phaseModel& phase1 = iter();
+ phaseModel& phase = iter();
surfaceScalarField& alphaPhi = alphaPhiCorrs[phasei];
- alphaPhi += upwind<scalar>(mesh_, phi_).flux(phase1);
+ alphaPhi += upwind<scalar>(mesh_, phi_).flux(phase);
+ phase.correctInflowOutflow(alphaPhi);
MULES::explicitSolve
(
geometricOneField(),
- phase1,
+ phase,
alphaPhi,
zeroField(),
zeroField()
);
- phase1.alphaPhi() += alphaPhi;
+ phase.alphaPhi() = alphaPhi;
- Info<< phase1.name() << " volume fraction, min, max = "
- << phase1.weightedAverage(mesh_.V()).value()
- << ' ' << min(phase1).value()
- << ' ' << max(phase1).value()
+ Info<< phase.name() << " volume fraction, min, max = "
+ << phase.weightedAverage(mesh_.V()).value()
+ << ' ' << min(phase).value()
+ << ' ' << max(phase).value()
<< endl;
- sumAlpha += phase1;
+ sumAlpha += phase;
phasei++;
}
@@ -215,6 +192,15 @@ void Foam::multiphaseSystem::solveAlphas()
<< ' ' << max(sumAlpha).value()
<< endl;
+ // Correct the sum of the phase-fractions to avoid 'drift'
+ volScalarField sumCorr(1.0 - sumAlpha);
+ forAllIter(PtrDictionary<phaseModel>, phases_, iter)
+ {
+ phaseModel& phase = iter();
+ volScalarField& alpha = phase;
+ alpha += alpha*sumCorr;
+ }
+
calcAlphas();
}
diff --git a/applications/solvers/multiphase/multiphaseEulerFoam/multiphaseSystem/phaseModel/phaseModel.C b/applications/solvers/multiphase/multiphaseEulerFoam/multiphaseSystem/phaseModel/phaseModel.C
index 4392e62ac9b..1751e404a47 100644
--- a/applications/solvers/multiphase/multiphaseEulerFoam/multiphaseSystem/phaseModel/phaseModel.C
+++ b/applications/solvers/multiphase/multiphaseEulerFoam/multiphaseSystem/phaseModel/phaseModel.C
@@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
- \\ / A nd | Copyright (C) 2011-2016 OpenFOAM Foundation
+ \\ / A nd | Copyright (C) 2011-2017 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@@ -236,6 +236,24 @@ bool Foam::phaseModel::read(const dictionary& phaseDict)
}
+void Foam::phaseModel::correctInflowOutflow(surfaceScalarField& alphaPhi) const
+{
+ surfaceScalarField::Boundary& alphaPhiBf = alphaPhi.boundaryFieldRef();
+ const volScalarField::Boundary& alphaBf = boundaryField();
+ const surfaceScalarField::Boundary& phiBf = phi().boundaryField();
+
+ forAll(alphaPhiBf, patchi)
+ {
+ fvsPatchScalarField& alphaPhip = alphaPhiBf[patchi];
+
+ if (!alphaPhip.coupled())
+ {
+ alphaPhip = phiBf[patchi]*alphaBf[patchi];
+ }
+ }
+}
+
+
Foam::tmp<Foam::volScalarField> Foam::phaseModel::d() const
{
return dPtr_().d();
diff --git a/applications/solvers/multiphase/multiphaseEulerFoam/multiphaseSystem/phaseModel/phaseModel.H b/applications/solvers/multiphase/multiphaseEulerFoam/multiphaseSystem/phaseModel/phaseModel.H
index d1f1157b9a1..90e48183cae 100644
--- a/applications/solvers/multiphase/multiphaseEulerFoam/multiphaseSystem/phaseModel/phaseModel.H
+++ b/applications/solvers/multiphase/multiphaseEulerFoam/multiphaseSystem/phaseModel/phaseModel.H
@@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
- \\ / A nd | Copyright (C) 2011-2015 OpenFOAM Foundation
+ \\ / A nd | Copyright (C) 2011-2017 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@@ -208,6 +208,9 @@ public:
return alphaPhi_;
}
+ //- Ensure that the flux at inflow/outflow BCs is preserved
+ void correctInflowOutflow(surfaceScalarField& alphaPhi) const;
+
//- Correct the phase properties
void correct();
diff --git a/applications/solvers/multiphase/multiphaseInterFoam/multiphaseMixture/multiphaseMixture.C b/applications/solvers/multiphase/multiphaseInterFoam/multiphaseMixture/multiphaseMixture.C
index 5dd59dead15..21b5edefe45 100644
--- a/applications/solvers/multiphase/multiphaseInterFoam/multiphaseMixture/multiphaseMixture.C
+++ b/applications/solvers/multiphase/multiphaseInterFoam/multiphaseMixture/multiphaseMixture.C
@@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
- \\ / A nd | Copyright (C) 2011-2016 OpenFOAM Foundation
+ \\ / A nd | Copyright (C) 2011-2017 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@@ -681,6 +681,14 @@ void Foam::multiphaseMixture::solveAlphas
<< ' ' << max(sumAlpha).value()
<< endl;
+ // Correct the sum of the phase-fractions to avoid 'drift'
+ volScalarField sumCorr(1.0 - sumAlpha);
+ forAllIter(PtrDictionary<phase>, phases_, iter)
+ {
+ phase& alpha = iter();
+ alpha += alpha*sumCorr;
+ }
+
calcAlphas();
}
diff --git a/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/phaseModel/phaseModel/phaseModel.C b/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/phaseModel/phaseModel/phaseModel.C
index 9fa7c0f5898..da10824230f 100644
--- a/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/phaseModel/phaseModel/phaseModel.C
+++ b/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/phaseModel/phaseModel/phaseModel.C
@@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
- \\ / A nd | Copyright (C) 2015-2016 OpenFOAM Foundation
+ \\ / A nd | Copyright (C) 2015-2017 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@@ -165,6 +165,24 @@ bool Foam::phaseModel::compressible() const
}
+void Foam::phaseModel::correctInflowOutflow(surfaceScalarField& alphaPhi) const
+{
+ surfaceScalarField::Boundary& alphaPhiBf = alphaPhi.boundaryFieldRef();
+ const volScalarField::Boundary& alphaBf = boundaryField();
+ const surfaceScalarField::Boundary& phiBf = phi()().boundaryField();
+
+ forAll(alphaPhiBf, patchi)
+ {
+ fvsPatchScalarField& alphaPhip = alphaPhiBf[patchi];
+
+ if (!alphaPhip.coupled())
+ {
+ alphaPhip = phiBf[patchi]*alphaBf[patchi];
+ }
+ }
+}
+
+
const Foam::tmp<Foam::volScalarField>& Foam::phaseModel::divU() const
{
NotImplemented;
diff --git a/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/phaseModel/phaseModel/phaseModel.H b/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/phaseModel/phaseModel/phaseModel.H
index 719ef80a8c9..41fb67e715d 100644
--- a/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/phaseModel/phaseModel/phaseModel.H
+++ b/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/phaseModel/phaseModel/phaseModel.H
@@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
- \\ / A nd | Copyright (C) 2015 OpenFOAM Foundation
+ \\ / A nd | Copyright (C) 2015-2017 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@@ -283,6 +283,9 @@ public:
//- Access the mass flux of the phase
virtual surfaceScalarField& alphaRhoPhi() = 0;
+ //- Ensure that the flux at inflow/outflow BCs is preserved
+ void correctInflowOutflow(surfaceScalarField& alphaPhi) const;
+
// Transport
diff --git a/applications/solvers/multiphase/reactingEulerFoam/reactingMultiphaseEulerFoam/multiphaseSystem/multiphaseSystem.C b/applications/solvers/multiphase/reactingEulerFoam/reactingMultiphaseEulerFoam/multiphaseSystem/multiphaseSystem.C
index ca37c2245b4..805d99844e3 100644
--- a/applications/solvers/multiphase/reactingEulerFoam/reactingMultiphaseEulerFoam/multiphaseSystem/multiphaseSystem.C
+++ b/applications/solvers/multiphase/reactingEulerFoam/reactingMultiphaseEulerFoam/multiphaseSystem/multiphaseSystem.C
@@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
- \\ / A nd | Copyright (C) 2013-2016 OpenFOAM Foundation
+ \\ / A nd | Copyright (C) 2013-2017 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@@ -71,15 +71,17 @@ void Foam::multiphaseSystem::solveAlphas()
{
bool LTS = fv::localEulerDdt::enabled(mesh_);
+ forAll(phases(), phasei)
+ {
+ phases()[phasei].correctBoundaryConditions();
+ }
+
PtrList<surfaceScalarField> alphaPhiCorrs(phases().size());
forAll(phases(), phasei)
{
phaseModel& phase = phases()[phasei];
volScalarField& alpha1 = phase;
- phase.alphaPhi() =
- dimensionedScalar("0", dimensionSet(0, 3, -1, 0, 0), 0);
-
alphaPhiCorrs.set
(
phasei,
@@ -134,28 +136,7 @@ void Foam::multiphaseSystem::solveAlphas()
);
}
- surfaceScalarField::Boundary& alphaPhiCorrBf =
- alphaPhiCorr.boundaryFieldRef();
-
- // Ensure that the flux at inflow BCs is preserved
- forAll(alphaPhiCorr.boundaryField(), patchi)
- {
- fvsPatchScalarField& alphaPhiCorrp = alphaPhiCorrBf[patchi];
-
- if (!alphaPhiCorrp.coupled())
- {
- const scalarField& phi1p = phase.phi().boundaryField()[patchi];
- const scalarField& alpha1p = alpha1.boundaryField()[patchi];
-
- forAll(alphaPhiCorrp, facei)
- {
- if (phi1p[facei] < 0)
- {
- alphaPhiCorrp[facei] = alpha1p[facei]*phi1p[facei];
- }
- }
- }
- }
+ phase.correctInflowOutflow(alphaPhiCorr);
if (LTS)
{
@@ -215,8 +196,9 @@ void Foam::multiphaseSystem::solveAlphas()
phaseModel& phase = phases()[phasei];
volScalarField& alpha = phase;
- surfaceScalarField& alphaPhic = alphaPhiCorrs[phasei];
- alphaPhic += upwind<scalar>(mesh_, phi_).flux(phase);
+ surfaceScalarField& alphaPhi = alphaPhiCorrs[phasei];
+ alphaPhi += upwind<scalar>(mesh_, phi_).flux(phase);
+ phase.correctInflowOutflow(alphaPhi);
volScalarField::Internal Sp
(
@@ -298,12 +280,12 @@ void Foam::multiphaseSystem::solveAlphas()
(
geometricOneField(),
alpha,
- alphaPhic,
+ alphaPhi,
Sp,
Su
);
- phase.alphaPhi() += alphaPhic;
+ phase.alphaPhi() = alphaPhi;
Info<< phase.name() << " volume fraction, min, max = "
<< phase.weightedAverage(mesh_.V()).value()
@@ -319,6 +301,14 @@ void Foam::multiphaseSystem::solveAlphas()
<< ' ' << min(sumAlpha).value()
<< ' ' << max(sumAlpha).value()
<< endl;
+
+ // Correct the sum of the phase-fractions to avoid 'drift'
+ volScalarField sumCorr(1.0 - sumAlpha);
+ forAll(phases(), phasei)
+ {
+ volScalarField& alpha = phases()[phasei];
+ alpha += alpha*sumCorr;
+ }
}
diff --git a/applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/twoPhaseSystem/twoPhaseSystem.C b/applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/twoPhaseSystem/twoPhaseSystem.C
index 8115f05cd0b..62bbe6af40d 100644
--- a/applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/twoPhaseSystem/twoPhaseSystem.C
+++ b/applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/twoPhaseSystem/twoPhaseSystem.C
@@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
- \\ / A nd | Copyright (C) 2013-2016 OpenFOAM Foundation
+ \\ / A nd | Copyright (C) 2013-2017 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@@ -201,7 +201,6 @@ void Foam::twoPhaseSystem::solve()
word alphaScheme("div(phi," + alpha1.name() + ')');
word alpharScheme("div(phir," + alpha1.name() + ')');
- const surfaceScalarField& phi = this->phi();
const surfaceScalarField& phi1 = phase1_.phi();
const surfaceScalarField& phi2 = phase2_.phi();
@@ -236,9 +235,7 @@ void Foam::twoPhaseSystem::solve()
}
alpha1.correctBoundaryConditions();
- surfaceScalarField alpha1f(fvc::interpolate(max(alpha1, scalar(0))));
- surfaceScalarField phic("phic", phi);
surfaceScalarField phir("phir", phi1 - phi2);
tmp<surfaceScalarField> alphaDbyA;
@@ -279,7 +276,7 @@ void Foam::twoPhaseSystem::solve()
),
// Divergence term is handled explicitly to be
// consistent with the explicit transport solution
- fvc::div(phi)*min(alpha1, scalar(1))
+ fvc::div(phi_)*min(alpha1, scalar(1))
);
if (tdgdt.valid())
@@ -300,11 +297,11 @@ void Foam::twoPhaseSystem::solve()
}
}
- surfaceScalarField alphaPhic1
+ surfaceScalarField alphaPhi1
(
fvc::flux
(
- phic,
+ phi_,
alpha1,
alphaScheme
)
@@ -316,28 +313,7 @@ void Foam::twoPhaseSystem::solve()
)
);
- surfaceScalarField::Boundary& alphaPhic1Bf =
- alphaPhic1.boundaryFieldRef();
-
- // Ensure that the flux at inflow BCs is preserved
- forAll(alphaPhic1Bf, patchi)
- {
- fvsPatchScalarField& alphaPhic1p = alphaPhic1Bf[patchi];
-
- if (!alphaPhic1p.coupled())
- {
- const scalarField& phi1p = phi1.boundaryField()[patchi];
- const scalarField& alpha1p = alpha1.boundaryField()[patchi];
-
- forAll(alphaPhic1p, facei)
- {
- if (phi1p[facei] < 0)
- {
- alphaPhic1p[facei] = alpha1p[facei]*phi1p[facei];
- }
- }
- }
- }
+ phase1_.correctInflowOutflow(alphaPhi1);
if (nAlphaSubCycles > 1)
{
@@ -355,14 +331,14 @@ void Foam::twoPhaseSystem::solve()
!(++alphaSubCycle).end();
)
{
- surfaceScalarField alphaPhic10(alphaPhic1);
+ surfaceScalarField alphaPhi10(alphaPhi1);
MULES::explicitSolve
(
geometricOneField(),
alpha1,
- phi,
- alphaPhic10,
+ phi_,
+ alphaPhi10,
(alphaSubCycle.index()*Sp)(),
(Su - (alphaSubCycle.index() - 1)*Sp*alpha1)(),
phase1_.alphaMax(),
@@ -371,11 +347,11 @@ void Foam::twoPhaseSystem::solve()
if (alphaSubCycle.index() == 1)
{
- phase1_.alphaPhi() = alphaPhic10;
+ phase1_.alphaPhi() = alphaPhi10;
}
else
{
- phase1_.alphaPhi() += alphaPhic10;
+ phase1_.alphaPhi() += alphaPhi10;
}
}
@@ -387,15 +363,15 @@ void Foam::twoPhaseSystem::solve()
(
geometricOneField(),
alpha1,
- phi,
- alphaPhic1,
+ phi_,
+ alphaPhi1,
Sp,
Su,
phase1_.alphaMax(),
0
);
- phase1_.alphaPhi() = alphaPhic1;
+ phase1_.alphaPhi() = alphaPhi1;
}
if (alphaDbyA.valid())
@@ -415,8 +391,8 @@ void Foam::twoPhaseSystem::solve()
phase1_.alphaRhoPhi() =
fvc::interpolate(phase1_.rho())*phase1_.alphaPhi();
- phase2_.alphaPhi() = phi - phase1_.alphaPhi();
- alpha2 = scalar(1) - alpha1;
+ phase2_.alphaPhi() = phi_ - phase1_.alphaPhi();
+ phase2_.correctInflowOutflow(phase2_.alphaPhi());
phase2_.alphaRhoPhi() =
fvc::interpolate(phase2_.rho())*phase2_.alphaPhi();
@@ -425,6 +401,13 @@ void Foam::twoPhaseSystem::solve()
<< " Min(alpha1) = " << min(alpha1).value()
<< " Max(alpha1) = " << max(alpha1).value()
<< endl;
+
+ // Ensure the phase-fractions are bounded
+ alpha1.max(0);
+ alpha1.min(1);
+
+ // Update the phase-fraction of the other phase
+ alpha2 = scalar(1) - alpha1;
}
}
diff --git a/applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseSystem/phaseModel/phaseModel.C b/applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseSystem/phaseModel/phaseModel.C
index cd60d2cbb00..dab711fb2f0 100644
--- a/applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseSystem/phaseModel/phaseModel.C
+++ b/applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseSystem/phaseModel/phaseModel.C
@@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
- \\ / A nd | Copyright (C) 2011-2016 OpenFOAM Foundation
+ \\ / A nd | Copyright (C) 2011-2017 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@@ -248,27 +248,19 @@ bool Foam::phaseModel::read(const dictionary& phaseProperties)
}
-void Foam::phaseModel::correctInflowFlux(surfaceScalarField& alphaPhi) const
+void Foam::phaseModel::correctInflowOutflow(surfaceScalarField& alphaPhi) const
{
surfaceScalarField::Boundary& alphaPhiBf = alphaPhi.boundaryFieldRef();
+ const volScalarField::Boundary& alphaBf = boundaryField();
+ const surfaceScalarField::Boundary& phiBf = phi().boundaryField();
- // Ensure that the flux at inflow BCs is preserved
forAll(alphaPhiBf, patchi)
{
fvsPatchScalarField& alphaPhip = alphaPhiBf[patchi];
if (!alphaPhip.coupled())
{
- const scalarField& phip = phi().boundaryField()[patchi];
- const scalarField& alphap = boundaryField()[patchi];
-
- forAll(alphaPhip, facei)
- {
- if (phip[facei] < SMALL)
- {
- alphaPhip[facei] = alphap[facei]*phip[facei];
- }
- }
+ alphaPhip = phiBf[patchi]*alphaBf[patchi];
}
}
}
diff --git a/applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseSystem/phaseModel/phaseModel.H b/applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseSystem/phaseModel/phaseModel.H
index a5ae107c909..894907fe7cc 100644
--- a/applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseSystem/phaseModel/phaseModel.H
+++ b/applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseSystem/phaseModel/phaseModel.H
@@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
- \\ / A nd | Copyright (C) 2011-2016 OpenFOAM Foundation
+ \\ / A nd | Copyright (C) 2011-2017 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@@ -319,8 +319,8 @@ public:
return alphaRhoPhi_;
}
- //- Ensure that the flux at inflow BCs is preserved
- void correctInflowFlux(surfaceScalarField& alphaPhi) const;
+ //- Ensure that the flux at inflow/outflow BCs is preserved
+ void correctInflowOutflow(surfaceScalarField& alphaPhi) const;
//- Correct the phase properties
// other than the thermodynamics and turbulence
diff --git a/applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseSystem/twoPhaseSystem.C b/applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseSystem/twoPhaseSystem.C
index 808fe255f9e..c6861dcf2c0 100644
--- a/applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseSystem/twoPhaseSystem.C
+++ b/applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseSystem/twoPhaseSystem.C
@@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
- \\ / A nd | Copyright (C) 2013-2016 OpenFOAM Foundation
+ \\ / A nd | Copyright (C) 2013-2017 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@@ -444,7 +444,7 @@ void Foam::twoPhaseSystem::solve()
)
);
- phase1_.correctInflowFlux(alphaPhic1);
+ phase1_.correctInflowOutflow(alphaPhic1);
if (nAlphaSubCycles > 1)
{
@@ -515,8 +515,7 @@ void Foam::twoPhaseSystem::solve()
fvc::interpolate(phase1_.rho())*phase1_.alphaPhi();
phase2_.alphaPhi() = phi_ - phase1_.alphaPhi();
- alpha2 = scalar(1) - alpha1;
- phase2_.correctInflowFlux(phase2_.alphaPhi());
+ phase2_.correctInflowOutflow(phase2_.alphaPhi());
phase2_.alphaRhoPhi() =
fvc::interpolate(phase2_.rho())*phase2_.alphaPhi();
@@ -525,6 +524,12 @@ void Foam::twoPhaseSystem::solve()
<< " Min(" << alpha1.name() << ") = " << min(alpha1).value()
<< " Max(" << alpha1.name() << ") = " << max(alpha1).value()
<< endl;
+
+ // Ensure the phase-fractions are bounded
+ alpha1.max(0);
+ alpha1.min(1);
+
+ alpha2 = scalar(1) - alpha1;
}
}
diff --git a/src/finiteVolume/Make/files b/src/finiteVolume/Make/files
index 20ef3e4f134..d79bbff4e82 100644
--- a/src/finiteVolume/Make/files
+++ b/src/finiteVolume/Make/files
@@ -235,7 +235,6 @@ fvMatrices/fvMatrices.C
fvMatrices/fvScalarMatrix/fvScalarMatrix.C
fvMatrices/solvers/MULES/MULES.C
fvMatrices/solvers/MULES/CMULES.C
-fvMatrices/solvers/MULES/IMULES.C
fvMatrices/solvers/GAMGSymSolver/GAMGAgglomerations/faceAreaPairGAMGAgglomeration/faceAreaPairGAMGAgglomeration.C
interpolation = interpolation/interpolation
diff --git a/src/finiteVolume/fvMatrices/solvers/MULES/CMULESTemplates.C b/src/finiteVolume/fvMatrices/solvers/MULES/CMULESTemplates.C
index 41fc3b1194d..34ee43bede3 100644
--- a/src/finiteVolume/fvMatrices/solvers/MULES/CMULESTemplates.C
+++ b/src/finiteVolume/fvMatrices/solvers/MULES/CMULESTemplates.C
@@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
- \\ / A nd | Copyright (C) 2013-2016 OpenFOAM Foundation
+ \\ / A nd | Copyright (C) 2013-2017 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@@ -151,17 +151,17 @@ void Foam::MULES::limiterCorr
const dictionary& MULEScontrols = mesh.solverDict(psi.name());
- label nLimiterIter
+ const label nLimiterIter
(
readLabel(MULEScontrols.lookup("nLimiterIter"))
);
- scalar smoothLimiter
+ const scalar smoothLimiter
(
MULEScontrols.lookupOrDefault<scalar>("smoothLimiter", 0)
);
- scalar extremaCoeff
+ const scalar extremaCoeff
(
MULEScontrols.lookupOrDefault<scalar>("extremaCoeff", 0)
);
@@ -207,8 +207,8 @@ void Foam::MULES::limiterCorr
forAll(phiCorrIf, facei)
{
- label own = owner[facei];
- label nei = neighb[facei];
+ const label own = owner[facei];
+ const label nei = neighb[facei];
psiMaxn[own] = max(psiMaxn[own], psiIf[nei]);
psiMinn[own] = min(psiMinn[own], psiIf[nei]);
@@ -216,9 +216,9 @@ void Foam::MULES::limiterCorr
psiMaxn[nei] = max(psiMaxn[nei], psiIf[own]);
psiMinn[nei] = min(psiMinn[nei], psiIf[own]);
- scalar phiCorrf = phiCorrIf[facei];
+ const scalar phiCorrf = phiCorrIf[facei];
- if (phiCorrf > 0.0)
+ if (phiCorrf > 0)
{
sumPhip[own] += phiCorrf;
mSumPhim[nei] += phiCorrf;
@@ -253,7 +253,7 @@ void Foam::MULES::limiterCorr
{
forAll(phiCorrPf, pFacei)
{
- label pfCelli = pFaceCells[pFacei];
+ const label pfCelli = pFaceCells[pFacei];
psiMaxn[pfCelli] = max(psiMaxn[pfCelli], psiPf[pFacei]);
psiMinn[pfCelli] = min(psiMinn[pfCelli], psiPf[pFacei]);
@@ -262,11 +262,11 @@ void Foam::MULES::limiterCorr
forAll(phiCorrPf, pFacei)
{
- label pfCelli = pFaceCells[pFacei];
+ const label pfCelli = pFaceCells[pFacei];
- scalar phiCorrf = phiCorrPf[pFacei];
+ const scalar phiCorrf = phiCorrPf[pFacei];
- if (phiCorrf > 0.0)
+ if (phiCorrf > 0)
{
sumPhip[pfCelli] += phiCorrf;
}
@@ -309,17 +309,17 @@ void Foam::MULES::limiterCorr
for (int j=0; j<nLimiterIter; j++)
{
- sumlPhip = 0.0;
- mSumlPhim = 0.0;
+ sumlPhip = 0;
+ mSumlPhim = 0;
forAll(lambdaIf, facei)
{
- label own = owner[facei];
- label nei = neighb[facei];
+ const label own = owner[facei];
+ const label nei = neighb[facei];
- scalar lambdaPhiCorrf = lambdaIf[facei]*phiCorrIf[facei];
+ const scalar lambdaPhiCorrf = lambdaIf[facei]*phiCorrIf[facei];
- if (lambdaPhiCorrf > 0.0)
+ if (lambdaPhiCorrf > 0)
{
sumlPhip[own] += lambdaPhiCorrf;
mSumlPhim[nei] += lambdaPhiCorrf;
@@ -344,7 +344,7 @@ void Foam::MULES::limiterCorr
scalar lambdaPhiCorrf = lambdaPf[pFacei]*phiCorrfPf[pFacei];
- if (lambdaPhiCorrf > 0.0)
+ if (lambdaPhiCorrf > 0)
{
sumlPhip[pfCelli] += lambdaPhiCorrf;
}
@@ -379,7 +379,7 @@ void Foam::MULES::limiterCorr
forAll(lambdaIf, facei)
{
- if (phiCorrIf[facei] > 0.0)
+ if (phiCorrIf[facei] > 0)
{
lambdaIf[facei] = min
(
@@ -415,9 +415,9 @@ void Foam::MULES::limiterCorr
forAll(lambdaPf, pFacei)
{
- label pfCelli = pFaceCells[pFacei];
+ const label pfCelli = pFaceCells[pFacei];
- if (phiCorrfPf[pFacei] > 0.0)
+ if (phiCorrfPf[pFacei] > 0)
{
lambdaPf[pFacei] =
min(lambdaPf[pFacei], lambdap[pfCelli]);
@@ -438,11 +438,11 @@ void Foam::MULES::limiterCorr
forAll(lambdaPf, pFacei)
{
// Limit outlet faces only
- if (phiPf[pFacei] > SMALL*SMALL)
+ if ((phiPf[pFacei] + phiCorrfPf[pFacei]) > SMALL*SMALL)
{
- label pfCelli = pFaceCells[pFacei];
+ const label pfCelli = pFaceCells[pFacei];
- if (phiCorrfPf[pFacei] > 0.0)
+ if (phiCorrfPf[pFacei] > 0)
{
lambdaPf[pFacei] =
min(lambdaPf[pFacei], lambdap[pfCelli]);
diff --git a/src/finiteVolume/fvMatrices/solvers/MULES/IMULES.C b/src/finiteVolume/fvMatrices/solvers/MULES/IMULES.C
deleted file mode 100644
index 2bc88decb55..00000000000
--- a/src/finiteVolume/fvMatrices/solvers/MULES/IMULES.C
+++ /dev/null
@@ -1,51 +0,0 @@
-/*---------------------------------------------------------------------------*\
- ========= |
- \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
- \\ / O peration |
- \\ / A nd | Copyright (C) 2013 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 "IMULES.H"
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-void Foam::MULES::implicitSolve
-(
- volScalarField& psi,
- const surfaceScalarField& phi,
- surfaceScalarField& phiPsi,
- const scalar psiMax,
- const scalar psiMin
-)
-{
- implicitSolve
- (
- geometricOneField(),
- psi,
- phi,
- phiPsi,
- zeroField(), zeroField(),
- psiMax, psiMin
- );
-}
-
-
-// ************************************************************************* //
diff --git a/src/finiteVolume/fvMatrices/solvers/MULES/IMULES.H b/src/finiteVolume/fvMatrices/solvers/MULES/IMULES.H
deleted file mode 100644
index fdef6b283b9..00000000000
--- a/src/finiteVolume/fvMatrices/solvers/MULES/IMULES.H
+++ /dev/null
@@ -1,94 +0,0 @@
-/*---------------------------------------------------------------------------*\
- ========= |
- \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
- \\ / O peration |
- \\ / A nd | Copyright (C) 2013-2016 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/>.
-
-Global
- IMULES
-
-Description
- IMULES: Multidimensional universal limiter for implicit solution.
-
- Solve a convective-only transport equation using an explicit universal
- multi-dimensional limiter applied to an implicit formulation requiring
- iteration to guarantee boundedness. The number of iterations required
- to obtain boundedness increases with the Courant number of the simulation.
-
- It may be more efficient to use CMULES.
-
-SourceFiles
- IMULES.C
- IMULESTemplates.C
-
-\*---------------------------------------------------------------------------*/
-
-#ifndef IMULES_H
-#define IMULES_H
-
-#include "MULES.H"
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-namespace Foam
-{
-namespace MULES
-{
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-template<class RhoType, class SpType, class SuType>
-void implicitSolve
-(
- const RhoType& rho,
- volScalarField& gamma,
- const surfaceScalarField& phi,
- surfaceScalarField& phiCorr,
- const SpType& Sp,
- const SuType& Su,
- const scalar psiMax,
- const scalar psiMin
-);
-
-void implicitSolve
-(
- volScalarField& gamma,
- const surfaceScalarField& phi,
- surfaceScalarField& phiCorr,
- const scalar psiMax,
- const scalar psiMin
-);
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-} // End namespace MULES
-} // End namespace Foam
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-#ifdef NoRepository
- #include "IMULESTemplates.C"
-#endif
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-#endif
-
-// ************************************************************************* //
diff --git a/src/finiteVolume/fvMatrices/solvers/MULES/IMULESTemplates.C b/src/finiteVolume/fvMatrices/solvers/MULES/IMULESTemplates.C
deleted file mode 100644
index c1ecbb927fa..00000000000
--- a/src/finiteVolume/fvMatrices/solvers/MULES/IMULESTemplates.C
+++ /dev/null
@@ -1,236 +0,0 @@
-/*---------------------------------------------------------------------------*\
- ========= |
- \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
- \\ / O peration |
- \\ / A nd | Copyright (C) 2011-2016 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 "IMULES.H"
-#include "gaussConvectionScheme.H"
-#include "surfaceInterpolate.H"
-#include "fvmDdt.H"
-#include "fvmSup.H"
-#include "fvcDiv.H"
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-namespace Foam
-{
-namespace MULES
-{
- template<class RhoType>
- inline tmp<surfaceScalarField> interpolate(const RhoType& rho)
- {
- NotImplemented;
- return tmp<surfaceScalarField>(nullptr);
- }
-
- template<>
- inline tmp<surfaceScalarField> interpolate(const volScalarField& rho)
- {
- return fvc::interpolate(rho);
- }
-}
-}
-
-
-template<class RhoType, class SpType, class SuType>
-void Foam::MULES::implicitSolve
-(
- const RhoType& rho,
- volScalarField& psi,
- const surfaceScalarField& phi,
- surfaceScalarField& phiPsi,
- 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 maxIter
- (
- readLabel(MULEScontrols.lookup("maxIter"))
- );
-
- scalar maxUnboundedness
- (
- readScalar(MULEScontrols.lookup("maxUnboundedness"))
- );
-
- scalar CoCoeff
- (
- readScalar(MULEScontrols.lookup("CoCoeff"))
- );
-
- scalarField allCoLambda(mesh.nFaces());
-
- {
- slicedSurfaceScalarField CoLambda
- (
- IOobject
- (
- "CoLambda",
- mesh.time().timeName(),
- mesh,
- IOobject::NO_READ,
- IOobject::NO_WRITE,
- false
- ),
- mesh,
- dimless,
- allCoLambda,
- false // Use slices for the couples
- );
-
- if (phi.dimensions() == dimDensity*dimVelocity*dimArea)
- {
- tmp<surfaceScalarField> Cof =
- mesh.time().deltaT()*mesh.surfaceInterpolation::deltaCoeffs()
- *mag(phi/interpolate(rho))/mesh.magSf();
-
- CoLambda == 1.0/max(CoCoeff*Cof, scalar(1));
- }
- else
- {
- tmp<surfaceScalarField> Cof =
- mesh.time().deltaT()*mesh.surfaceInterpolation::deltaCoeffs()
- *mag(phi)/mesh.magSf();
-
- CoLambda == 1.0/max(CoCoeff*Cof, scalar(1));
- }
- }
-
- scalarField allLambda(allCoLambda);
- //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
- );
-
- linear<scalar> CDs(mesh);
- upwind<scalar> UDs(mesh, phi);
- //fv::uncorrectedSnGrad<scalar> snGrads(mesh);
-
- fvScalarMatrix psiConvectionDiffusion
- (
- fvm::ddt(rho, psi)
- + fv::gaussConvectionScheme<scalar>(mesh, phi, UDs).fvmDiv(phi, psi)
- //- fv::gaussLaplacianScheme<scalar, scalar>(mesh, CDs, snGrads)
- //.fvmLaplacian(Dpsif, psi)
- - fvm::Sp(Sp, psi)
- - Su
- );
-
- surfaceScalarField phiBD(psiConvectionDiffusion.flux());
-
- surfaceScalarField& phiCorr = phiPsi;
- phiCorr -= phiBD;
-
- for (label i=0; i<maxIter; i++)
- {
- if (i != 0 && i < 4)
- {
- allLambda = allCoLambda;
- }
-
- limiter
- (
- allLambda,
- 1.0/mesh.time().deltaTValue(),
- rho,
- psi,
- phiBD,
- phiCorr,
- Sp,
- Su,
- psiMax,
- psiMin
- );
-
- solve
- (
- psiConvectionDiffusion + fvc::div(lambda*phiCorr),
- MULEScontrols
- );
-
- scalar maxPsiM1 = gMax(psi.primitiveField()) - 1.0;
- scalar minPsi = gMin(psi.primitiveField());
-
- scalar unboundedness = max(max(maxPsiM1, 0.0), -min(minPsi, 0.0));
-
- if (unboundedness < maxUnboundedness)
- {
- break;
- }
- else
- {
- Info<< "MULES: max(" << psi.name() << " - 1) = " << maxPsiM1
- << " min(" << psi.name() << ") = " << minPsi << endl;
-
- phiBD = psiConvectionDiffusion.flux();
-
- /*
- word gammaScheme("div(phi,gamma)");
- word gammarScheme("div(phirb,gamma)");
-
- const surfaceScalarField& phir =
- mesh.lookupObject<surfaceScalarField>("phir");
-
- phiCorr =
- fvc::flux
- (
- phi,
- psi,
- gammaScheme
- )
- + fvc::flux
- (
- -fvc::flux(-phir, scalar(1) - psi, gammarScheme),
- psi,
- gammarScheme
- )
- - phiBD;
- */
- }
- }
-
- phiPsi = psiConvectionDiffusion.flux() + lambda*phiCorr;
-}
-
-
-// ************************************************************************* //
diff --git a/src/finiteVolume/fvMatrices/solvers/MULES/MULESTemplates.C b/src/finiteVolume/fvMatrices/solvers/MULES/MULESTemplates.C
index 3272d81d9f9..743cb8eae5c 100644
--- a/src/finiteVolume/fvMatrices/solvers/MULES/MULESTemplates.C
+++ b/src/finiteVolume/fvMatrices/solvers/MULES/MULESTemplates.C
@@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
- \\ / A nd | Copyright (C) 2011-2016 OpenFOAM Foundation
+ \\ / A nd | Copyright (C) 2011-2017 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@@ -185,12 +185,12 @@ void Foam::MULES::limiter
const dictionary& MULEScontrols = mesh.solverDict(psi.name());
- label nLimiterIter
+ const label nLimiterIter
(
MULEScontrols.lookupOrDefault<label>("nLimiterIter", 3)
);
- scalar smoothLimiter
+ const scalar smoothLimiter
(
MULEScontrols.lookupOrDefault<scalar>("smoothLimiter", 0)
);
@@ -228,8 +228,7 @@ void Foam::MULES::limiter
);
scalarField& lambdaIf = lambda;
- surfaceScalarField::Boundary& lambdaBf =
- lambda.boundaryFieldRef();
+ surfaceScalarField::Boundary& lambdaBf = lambda.boundaryFieldRef();
scalarField psiMaxn(psiIf.size(), psiMin);
scalarField psiMinn(psiIf.size(), psiMax);
@@ -241,8 +240,8 @@ void Foam::MULES::limiter
forAll(phiCorrIf, facei)
{
- label own = owner[facei];
- label nei = neighb[facei];
+ const label own = owner[facei];
+ const label nei = neighb[facei];
psiMaxn[own] = max(psiMaxn[own], psiIf[nei]);
psiMinn[own] = min(psiMinn[own], psiIf[nei]);
@@ -253,9 +252,9 @@ void Foam::MULES::limiter
sumPhiBD[own] += phiBDIf[facei];
sumPhiBD[nei] -= phiBDIf[facei];
- scalar phiCorrf = phiCorrIf[facei];
+ const scalar phiCorrf = phiCorrIf[facei];
- if (phiCorrf > 0.0)
+ if (phiCorrf > 0)
{
sumPhip[own] += phiCorrf;
mSumPhim[nei] += phiCorrf;
@@ -281,7 +280,7 @@ void Foam::MULES::limiter
forAll(phiCorrPf, pFacei)
{
- label pfCelli = pFaceCells[pFacei];
+ const label pfCelli = pFaceCells[pFacei];
psiMaxn[pfCelli] = max(psiMaxn[pfCelli], psiPNf[pFacei]);
psiMinn[pfCelli] = min(psiMinn[pfCelli], psiPNf[pFacei]);
@@ -291,7 +290,7 @@ void Foam::MULES::limiter
{
forAll(phiCorrPf, pFacei)
{
- label pfCelli = pFaceCells[pFacei];
+ const label pfCelli = pFaceCells[pFacei];
psiMaxn[pfCelli] = max(psiMaxn[pfCelli], psiPf[pFacei]);
psiMinn[pfCelli] = min(psiMinn[pfCelli], psiPf[pFacei]);
@@ -300,13 +299,13 @@ void Foam::MULES::limiter
forAll(phiCorrPf, pFacei)
{
- label pfCelli = pFaceCells[pFacei];
+ const label pfCelli = pFaceCells[pFacei];
sumPhiBD[pfCelli] += phiBDPf[pFacei];
- scalar phiCorrf = phiCorrPf[pFacei];
+ const scalar phiCorrf = phiCorrPf[pFacei];
- if (phiCorrf > 0.0)
+ if (phiCorrf > 0)
{
sumPhip[pfCelli] += phiCorrf;
}
@@ -376,17 +375,17 @@ void Foam::MULES::limiter
for (int j=0; j<nLimiterIter; j++)
{
- sumlPhip = 0.0;
- mSumlPhim = 0.0;
+ sumlPhip = 0;
+ mSumlPhim = 0;
forAll(lambdaIf, facei)
{
- label own = owner[facei];
- label nei = neighb[facei];
+ const label own = owner[facei];
+ const label nei = neighb[facei];
scalar lambdaPhiCorrf = lambdaIf[facei]*phiCorrIf[facei];
- if (lambdaPhiCorrf > 0.0)
+ if (lambdaPhiCorrf > 0)
{
sumlPhip[own] += lambdaPhiCorrf;
mSumlPhim[nei] += lambdaPhiCorrf;
@@ -407,11 +406,11 @@ void Foam::MULES::limiter
forAll(lambdaPf, pFacei)
{
- label pfCelli = pFaceCells[pFacei];
+ const label pfCelli = pFaceCells[pFacei];
+ const scalar lambdaPhiCorrf =
+ lambdaPf[pFacei]*phiCorrfPf[pFacei];
- scalar lambdaPhiCorrf = lambdaPf[pFacei]*phiCorrfPf[pFacei];
-
- if (lambdaPhiCorrf > 0.0)
+ if (lambdaPhiCorrf > 0)
{
sumlPhip[pfCelli] += lambdaPhiCorrf;
}
@@ -446,7 +445,7 @@ void Foam::MULES::limiter
forAll(lambdaIf, facei)
{
- if (phiCorrIf[facei] > 0.0)
+ if (phiCorrIf[facei] > 0)
{
lambdaIf[facei] = min
(
@@ -464,7 +463,6 @@ void Foam::MULES::limiter
}
}
-
forAll(lambdaBf, patchi)
{
fvsPatchScalarField& lambdaPf = lambdaBf[patchi];
@@ -482,9 +480,9 @@ void Foam::MULES::limiter
forAll(lambdaPf, pFacei)
{
- label pfCelli = pFaceCells[pFacei];
+ const label pfCelli = pFaceCells[pFacei];
- if (phiCorrfPf[pFacei] > 0.0)
+ if (phiCorrfPf[pFacei] > 0)
{
lambdaPf[pFacei] =
min(lambdaPf[pFacei], lambdap[pfCelli]);
@@ -496,33 +494,6 @@ void Foam::MULES::limiter
}
}
}
- else
- {
- const labelList& pFaceCells =
- mesh.boundary()[patchi].faceCells();
- const scalarField& phiBDPf = phiBDBf[patchi];
- const scalarField& phiCorrPf = phiCorrBf[patchi];
-
- forAll(lambdaPf, pFacei)
- {
- // Limit outlet faces only
- if ((phiBDPf[pFacei] + phiCorrPf[pFacei]) > SMALL*SMALL)
- {
- 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>());
@@ -549,6 +520,19 @@ void Foam::MULES::limit
surfaceScalarField phiBD(upwind<scalar>(psi.mesh(), phi).flux(psi));
+ surfaceScalarField::Boundary& phiBDBf = phiBD.boundaryFieldRef();
+ const surfaceScalarField::Boundary& phiPsiBf = phiPsi.boundaryField();
+
+ forAll(phiBDBf, patchi)
+ {
+ fvsPatchScalarField& phiBDPf = phiBDBf[patchi];
+
+ if (!phiBDPf.coupled())
+ {
+ phiBDPf = phiPsiBf[patchi];
+ }
+ }
+
surfaceScalarField& phiCorr = phiPsi;
phiCorr -= phiBD;
diff --git a/tutorials/multiphase/compressibleMultiphaseInterFoam/laminar/damBreak4phase/system/fvSolution b/tutorials/multiphase/compressibleMultiphaseInterFoam/laminar/damBreak4phase/system/fvSolution
index ab14e5ae1c2..7cf09c590a5 100644
--- a/tutorials/multiphase/compressibleMultiphaseInterFoam/laminar/damBreak4phase/system/fvSolution
+++ b/tutorials/multiphase/compressibleMultiphaseInterFoam/laminar/damBreak4phase/system/fvSolution
@@ -20,7 +20,7 @@ solvers
"alpha.*"
{
nAlphaSubCycles 4;
- cAlpha 2;
+ cAlpha 1;
}
pcorr
diff --git a/tutorials/multiphase/interDyMFoam/RAS/testTubeMixer/system/fvSolution b/tutorials/multiphase/interDyMFoam/RAS/testTubeMixer/system/fvSolution
index 7ebbdec3bd8..de2c8b8553c 100644
--- a/tutorials/multiphase/interDyMFoam/RAS/testTubeMixer/system/fvSolution
+++ b/tutorials/multiphase/interDyMFoam/RAS/testTubeMixer/system/fvSolution
@@ -21,7 +21,7 @@ solvers
{
nAlphaCorr 1;
nAlphaSubCycles 3;
- cAlpha 1.5;
+ cAlpha 1;
}
pcorr
diff --git a/tutorials/multiphase/interFoam/LES/nozzleFlow2D/system/fvSolution b/tutorials/multiphase/interFoam/LES/nozzleFlow2D/system/fvSolution
index ca68f313d5d..c095e44786c 100644
--- a/tutorials/multiphase/interFoam/LES/nozzleFlow2D/system/fvSolution
+++ b/tutorials/multiphase/interFoam/LES/nozzleFlow2D/system/fvSolution
@@ -21,7 +21,7 @@ solvers
{
nAlphaCorr 1;
nAlphaSubCycles 4;
- cAlpha 2;
+ cAlpha 1;
}
pcorr
diff --git a/tutorials/multiphase/interFoam/laminar/damBreak/damBreak/system/fvSolution b/tutorials/multiphase/interFoam/laminar/damBreak/damBreak/system/fvSolution
index fa206f8f198..81406469091 100644
--- a/tutorials/multiphase/interFoam/laminar/damBreak/damBreak/system/fvSolution
+++ b/tutorials/multiphase/interFoam/laminar/damBreak/damBreak/system/fvSolution
@@ -24,7 +24,7 @@ solvers
cAlpha 1;
MULESCorr yes;
- nLimiterIter 3;
+ nLimiterIter 5;
solver smoothSolver;
smoother symGaussSeidel;
diff --git a/tutorials/multiphase/multiphaseInterFoam/laminar/damBreak4phase/system/fvSolution b/tutorials/multiphase/multiphaseInterFoam/laminar/damBreak4phase/system/fvSolution
index fc08f6f97d3..b962745e2aa 100644
--- a/tutorials/multiphase/multiphaseInterFoam/laminar/damBreak4phase/system/fvSolution
+++ b/tutorials/multiphase/multiphaseInterFoam/laminar/damBreak4phase/system/fvSolution
@@ -20,7 +20,7 @@ solvers
"alpha.*"
{
nAlphaSubCycles 4;
- cAlpha 2;
+ cAlpha 1;
}
pcorr
diff --git a/tutorials/multiphase/multiphaseInterFoam/laminar/damBreak4phaseFine/system/fvSolution b/tutorials/multiphase/multiphaseInterFoam/laminar/damBreak4phaseFine/system/fvSolution
index a9e02cb6005..a271af11498 100644
--- a/tutorials/multiphase/multiphaseInterFoam/laminar/damBreak4phaseFine/system/fvSolution
+++ b/tutorials/multiphase/multiphaseInterFoam/laminar/damBreak4phaseFine/system/fvSolution
@@ -20,7 +20,7 @@ solvers
"alpha.*"
{
nAlphaSubCycles 4;
- cAlpha 2;
+ cAlpha 1;
}
pcorr
diff --git a/tutorials/multiphase/multiphaseInterFoam/laminar/mixerVessel2D/system/fvSolution b/tutorials/multiphase/multiphaseInterFoam/laminar/mixerVessel2D/system/fvSolution
index 26a94ec389a..64bf5c641c3 100644
--- a/tutorials/multiphase/multiphaseInterFoam/laminar/mixerVessel2D/system/fvSolution
+++ b/tutorials/multiphase/multiphaseInterFoam/laminar/mixerVessel2D/system/fvSolution
@@ -21,7 +21,7 @@ solvers
{
nAlphaCorr 4;
nAlphaSubCycles 4;
- cAlpha 2;
+ cAlpha 1;
}
pcorr
--
GitLab