From ad92287afca99ae73da43022fe8021026507fec6 Mon Sep 17 00:00:00 2001 From: Henry Weller 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 be4def509..798e9c2eb 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, 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, phases_, iter) { - phaseModel& phase1 = iter(); + phaseModel& phase = iter(); surfaceScalarField& alphaPhi = alphaPhiCorrs[phasei]; - alphaPhi += upwind(mesh_, phi_).flux(phase1); + alphaPhi += upwind(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, 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 4392e62ac..1751e404a 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::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 d1f1157b9..90e48183c 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 5dd59dead..21b5edefe 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, 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 9fa7c0f58..da1082423 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::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 719ef80a8..41fb67e71 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 ca37c2245..805d99844 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 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(mesh_, phi_).flux(phase); + surfaceScalarField& alphaPhi = alphaPhiCorrs[phasei]; + alphaPhi += upwind(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 8115f05cd..62bbe6af4 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 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 cd60d2cbb..dab711fb2 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 a5ae107c9..894907fe7 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 808fe255f..c6861dcf2 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 20ef3e4f1..d79bbff4e 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 41fc3b119..34ee43bed 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("smoothLimiter", 0) ); - scalar extremaCoeff + const scalar extremaCoeff ( MULEScontrols.lookupOrDefault("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 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 2bc88decb..000000000 --- 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 . - -\*---------------------------------------------------------------------------*/ - -#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 fdef6b283..000000000 --- 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 . - -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 -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 c1ecbb927..000000000 --- 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 . - -\*---------------------------------------------------------------------------*/ - -#include "IMULES.H" -#include "gaussConvectionScheme.H" -#include "surfaceInterpolate.H" -#include "fvmDdt.H" -#include "fvmSup.H" -#include "fvcDiv.H" - -// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // - -namespace Foam -{ -namespace MULES -{ - template - inline tmp interpolate(const RhoType& rho) - { - NotImplemented; - return tmp(nullptr); - } - - template<> - inline tmp interpolate(const volScalarField& rho) - { - return fvc::interpolate(rho); - } -} -} - - -template -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 Cof = - mesh.time().deltaT()*mesh.surfaceInterpolation::deltaCoeffs() - *mag(phi/interpolate(rho))/mesh.magSf(); - - CoLambda == 1.0/max(CoCoeff*Cof, scalar(1)); - } - else - { - tmp 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 CDs(mesh); - upwind UDs(mesh, phi); - //fv::uncorrectedSnGrad snGrads(mesh); - - fvScalarMatrix psiConvectionDiffusion - ( - fvm::ddt(rho, psi) - + fv::gaussConvectionScheme(mesh, phi, UDs).fvmDiv(phi, psi) - //- fv::gaussLaplacianScheme(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("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 3272d81d9..743cb8eae 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