From 35565437f5ffb61f90cf166e9707c19988a9bab5 Mon Sep 17 00:00:00 2001 From: Henry Weller Date: Wed, 5 Feb 2020 10:56:12 +0000 Subject: [PATCH] reactingMultiphaseEulerFoam::multiphaseSystem: Updated implicitPhasePressure handling to improve stability and allow larger time-steps --- .../multiphaseSystem/multiphaseSystem.C | 697 +++++++++--------- .../multiphaseSystem/multiphaseSystem.H | 8 +- .../twoPhaseSystem/twoPhaseSystem.C | 8 +- .../laminar/mixerVessel2D/system/fvSolution | 1 + 4 files changed, 353 insertions(+), 361 deletions(-) diff --git a/applications/solvers/multiphase/reactingEulerFoam/reactingMultiphaseEulerFoam/multiphaseSystem/multiphaseSystem.C b/applications/solvers/multiphase/reactingEulerFoam/reactingMultiphaseEulerFoam/multiphaseSystem/multiphaseSystem.C index 66754f19b4..662aef5d8b 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 | Website: https://openfoam.org - \\ / A nd | Copyright (C) 2013-2019 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2013-2020 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -68,314 +68,6 @@ void Foam::multiphaseSystem::calcAlphas() } -void Foam::multiphaseSystem::solveAlphas -( - const PtrList& rAUs, - const PtrList& rAUfs -) -{ - forAll(phases(), phasei) - { - phases()[phasei].correctBoundaryConditions(); - } - - // Calculate the void fraction - volScalarField alphaVoid - ( - IOobject - ( - "alphaVoid", - mesh_.time().timeName(), - mesh_ - ), - mesh_, - dimensionedScalar(dimless, 1) - ); - forAll(stationaryPhases(), stationaryPhasei) - { - alphaVoid -= stationaryPhases()[stationaryPhasei]; - } - - // Generate face-alphas - PtrList alphafs(phases().size()); - forAll(phases(), phasei) - { - phaseModel& phase = phases()[phasei]; - alphafs.set - ( - phasei, - new surfaceScalarField - ( - IOobject::groupName("alphaf", phase.name()), - upwind(mesh_, phi_).interpolate(phase) - ) - ); - } - - // Create correction fluxes - PtrList alphaPhiCorrs(phases().size()); - forAll(stationaryPhases(), stationaryPhasei) - { - phaseModel& phase = stationaryPhases()[stationaryPhasei]; - - alphaPhiCorrs.set - ( - phase.index(), - new surfaceScalarField - ( - IOobject::groupName("alphaPhiCorr", phase.name()), - - upwind(mesh_, phi_).flux(phase) - ) - ); - } - - PtrList DbyAs; - if (implicitPhasePressure() && (rAUs.size() || rAUfs.size())) - { - DbyAs = this->DByAfs(rAUs, rAUfs); - } - - PtrList alphaDbyAs(phases().size()); - - forAll(movingPhases(), movingPhasei) - { - phaseModel& phase = movingPhases()[movingPhasei]; - volScalarField& alpha = phase; - - alphaPhiCorrs.set - ( - phase.index(), - new surfaceScalarField - ( - IOobject::groupName("alphaPhiCorr", phase.name()), - fvc::flux(phi_, alpha, "div(phi," + alpha.name() + ')') - ) - ); - - surfaceScalarField& alphaPhiCorr = alphaPhiCorrs[phase.index()]; - - forAll(phases(), phasei) - { - phaseModel& phase2 = phases()[phasei]; - volScalarField& alpha2 = phase2; - - if (&phase2 == &phase) continue; - - surfaceScalarField phir(phase.phi() - phase2.phi()); - - cAlphaTable::const_iterator cAlpha - ( - cAlphas_.find(phasePairKey(phase.name(), phase2.name())) - ); - - if (cAlpha != cAlphas_.end()) - { - surfaceScalarField phic - ( - (mag(phi_) + mag(phir))/mesh_.magSf() - ); - - phir += min(cAlpha()*phic, max(phic))*nHatf(alpha, alpha2); - } - - word phirScheme - ( - "div(phir," + alpha2.name() + ',' + alpha.name() + ')' - ); - - alphaPhiCorr += fvc::flux - ( - -fvc::flux(-phir, alpha2, phirScheme), - alpha, - phirScheme - ); - } - - if (implicitPhasePressure() && (rAUs.size() || rAUfs.size())) - { - alphaDbyAs.set - ( - phase.index(), - fvc::interpolate(max(alpha, scalar(0))) - *fvc::interpolate(max(1 - alpha, scalar(0))) - *DbyAs[phase.index()] - ); - - alphaPhiCorr += - alphaDbyAs[phase.index()] - *fvc::snGrad(alpha, "bounded")*mesh_.magSf(); - } - - phase.correctInflowOutflow(alphaPhiCorr); - - MULES::limit - ( - geometricOneField(), - alpha, - phi_, - alphaPhiCorr, - zeroField(), - zeroField(), - min(alphaVoid.primitiveField(), phase.alphaMax())(), - zeroField(), - true - ); - } - - // Limit the flux sums, fixing those of the stationary phases - labelHashSet fixedAlphaPhiCorrs; - forAll(stationaryPhases(), stationaryPhasei) - { - fixedAlphaPhiCorrs.insert(stationaryPhases()[stationaryPhasei].index()); - } - MULES::limitSum(alphafs, alphaPhiCorrs, fixedAlphaPhiCorrs); - - // Solve for the moving phase alphas - forAll(movingPhases(), movingPhasei) - { - phaseModel& phase = movingPhases()[movingPhasei]; - volScalarField& alpha = phase; - - surfaceScalarField& alphaPhi = alphaPhiCorrs[phase.index()]; - alphaPhi += upwind(mesh_, phi_).flux(phase); - phase.correctInflowOutflow(alphaPhi); - - volScalarField::Internal Sp - ( - IOobject - ( - "Sp", - mesh_.time().timeName(), - mesh_ - ), - mesh_, - dimensionedScalar(dimless/dimTime, 0) - ); - - volScalarField::Internal Su - ( - "Su", - min(alpha, scalar(1))*fvc::div(fvc::absolute(phi_, phase.U())) - ); - - if (phase.divU().valid()) - { - const scalarField& dgdt = phase.divU()(); - - forAll(dgdt, celli) - { - if (dgdt[celli] > 0.0) - { - Sp[celli] -= dgdt[celli]; - Su[celli] += dgdt[celli]; - } - else if (dgdt[celli] < 0.0) - { - Sp[celli] += - dgdt[celli] - *(1 - alpha[celli])/max(alpha[celli], 1e-4); - } - } - } - - forAll(phases(), phasej) - { - const phaseModel& phase2 = phases()[phasej]; - const volScalarField& alpha2 = phase2; - - if (&phase2 == &phase) continue; - - if (phase2.divU().valid()) - { - const scalarField& dgdt2 = phase2.divU()(); - - forAll(dgdt2, celli) - { - if (dgdt2[celli] < 0.0) - { - Sp[celli] += - dgdt2[celli] - *(1 - alpha2[celli])/max(alpha2[celli], 1e-4); - - Su[celli] -= - dgdt2[celli] - *alpha[celli]/max(alpha2[celli], 1e-4); - } - else if (dgdt2[celli] > 0.0) - { - Sp[celli] -= dgdt2[celli]; - } - } - } - } - - MULES::explicitSolve - ( - geometricOneField(), - alpha, - alphaPhi, - Sp, - Su - ); - - phase.alphaPhiRef() = alphaPhi; - - if (alphaDbyAs.set(phase.index())) - { - fvScalarMatrix alphaEqn - ( - fvm::ddt(alpha) - fvc::ddt(alpha) - - fvm::laplacian(alphaDbyAs[phase.index()], alpha, "bounded") - ); - - alphaEqn.solve(); - - phase.alphaPhiRef() += alphaEqn.flux(); - } - } - - // Report the phase fractions and the phase fraction sum - forAll(phases(), phasei) - { - phaseModel& phase = phases()[phasei]; - - Info<< phase.name() << " fraction, min, max = " - << phase.weightedAverage(mesh_.V()).value() - << ' ' << min(phase).value() - << ' ' << max(phase).value() - << endl; - } - - volScalarField sumAlphaMoving - ( - IOobject - ( - "sumAlphaMoving", - mesh_.time().timeName(), - mesh_ - ), - mesh_, - dimensionedScalar(dimless, 0) - ); - forAll(movingPhases(), movingPhasei) - { - sumAlphaMoving += movingPhases()[movingPhasei]; - } - - Info<< "Phase-sum volume fraction, min, max = " - << (sumAlphaMoving + 1 - alphaVoid)().weightedAverage(mesh_.V()).value() - << ' ' << min(sumAlphaMoving + 1 - alphaVoid).value() - << ' ' << max(sumAlphaMoving + 1 - alphaVoid).value() - << endl; - - // Correct the sum of the phase fractions to avoid drift - forAll(movingPhases(), movingPhasei) - { - movingPhases()[movingPhasei] *= alphaVoid/sumAlphaMoving; - } -} - - Foam::tmp Foam::multiphaseSystem::nHatfv ( const volScalarField& alpha1, @@ -660,49 +352,53 @@ void Foam::multiphaseSystem::solve const PtrList& rAUfs ) { - const Time& runTime = mesh_.time(); - const dictionary& alphaControls = mesh_.solverDict("alpha"); + label nAlphaSubCycles(alphaControls.lookup