From 9f54506fbf7432dd7fe58824dbee8f9a88299417 Mon Sep 17 00:00:00 2001 From: Henry Weller Date: Thu, 4 Jan 2018 17:23:00 +0000 Subject: [PATCH] reactingEulerFoam: improvements to population balance modeling Removed possibility for the user to specify a driftRate in the constantDrift model which is independent of a fvOptions mass source. The driftRate must be calculated from/be consistent with the mass source in order to yield a particle number conserving result. Made calculation of the over-all Sauter mean diameter of an entire population balance conditional on more than one velocityGroup being present. This diameter field is for post-processing purposes only and would be redundant in case of one velocityGroup being used. Solution control is extended to allow for solution of the population balance equation at the last PIMPLE loop only, using an optional switch. This can be beneficial in terms of simulation time as well as coupling between the population balance based diameter calculation and the rest of the equation system. Patch contributed by Institute of Fluid Dynamics, Helmholtz-Zentrum Dresden - Rossendorf (HZDR) and VTT Technical Research Centre of Finland Ltd. --- .../driftModels/constantDrift/constantDrift.C | 21 +- .../driftModels/constantDrift/constantDrift.H | 16 +- .../populationBalanceModel.C | 208 ++++++++++-------- .../populationBalanceModel.H | 6 +- .../constant/phaseProperties | 3 - .../constant/phaseProperties | 3 - .../system/fvSolution | 7 +- 7 files changed, 134 insertions(+), 130 deletions(-) diff --git a/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/populationBalanceModel/driftModels/constantDrift/constantDrift.C b/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/populationBalanceModel/driftModels/constantDrift/constantDrift.C index 86d5a28b4..43c461648 100644 --- a/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/populationBalanceModel/driftModels/constantDrift/constantDrift.C +++ b/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/populationBalanceModel/driftModels/constantDrift/constantDrift.C @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | Copyright (C) 2017 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2017-2018 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -51,8 +51,6 @@ Foam::diameterModels::driftModels::constantDrift::constantDrift ) : driftModel(popBal, dict), - calculateFromFvOptions_(dict.lookup("calculateFromFvOptions")), - rate_("rate", dimVolume/dimTime, dict.lookupOrDefault("rate", 1.0)), N_ ( IOobject @@ -64,12 +62,7 @@ Foam::diameterModels::driftModels::constantDrift::constantDrift popBal.mesh(), dimensionedScalar("Sui", inv(dimVolume), Zero) ) -{ - if (calculateFromFvOptions_ == false) - { - rate_.read(dict); - } -} +{} // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // @@ -97,15 +90,7 @@ void Foam::diameterModels::driftModels::constantDrift::driftRate phaseModel& phase = const_cast(fi.phase()); volScalarField& rho = phase.thermo().rho(); - if (calculateFromFvOptions_ == true) - { - driftRate +=(popBal_.fluid().fvOptions()(phase, rho)&rho) - /(N_*rho); - } - else - { - driftRate += rate_; - } + driftRate += (popBal_.fluid().fvOptions()(phase, rho)&rho)/(N_*rho); } diff --git a/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/populationBalanceModel/driftModels/constantDrift/constantDrift.H b/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/populationBalanceModel/driftModels/constantDrift/constantDrift.H index b9d2cead0..1c7f8badd 100644 --- a/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/populationBalanceModel/driftModels/constantDrift/constantDrift.H +++ b/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/populationBalanceModel/driftModels/constantDrift/constantDrift.H @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | Copyright (C) 2017 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2017-2018 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -25,10 +25,10 @@ Class Foam::diameterModels::driftModels::constant Description - Constant growth rate. Used for verification and validation of the drift - formulation implemented in the populationBalanceModel class. Rate is either - given in the corresponding dictionary or calculated from fvOptions mass - source, depending on Switch calculateFromFvOptions. + Constant drift rate within all classes. Used for verification and + validation of the drift formulation implemented in the + populationBalanceModel class. Rate is calculated from fvOptions mass source. + SourceFiles constant.C @@ -59,12 +59,6 @@ class constantDrift { // Private data - //- If true, calculate drift rate from fvOptions mass source - Switch calculateFromFvOptions_; - - //- Optional driftRate - dimensionedScalar rate_; - //- Total number concentration volScalarField N_; diff --git a/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/populationBalanceModel/populationBalanceModel/populationBalanceModel.C b/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/populationBalanceModel/populationBalanceModel/populationBalanceModel.C index 581beb975..2559e7cb2 100644 --- a/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/populationBalanceModel/populationBalanceModel/populationBalanceModel.C +++ b/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/populationBalanceModel/populationBalanceModel/populationBalanceModel.C @@ -969,6 +969,7 @@ Foam::diameterModels::populationBalanceModel::populationBalanceModel ( fluid.subDict("populationBalanceCoeffs").subDict(name_) ), + pimple_(mesh_.lookupObject("solutionControl")), continuousPhase_ ( mesh_.lookupObject @@ -1044,24 +1045,7 @@ Foam::diameterModels::populationBalanceModel::populationBalanceModel Zero ) ), - d_ - ( - IOobject - ( - IOobject::groupName("d", this->name()), - fluid.time().timeName(), - fluid.mesh(), - IOobject::NO_READ, - IOobject::AUTO_WRITE - ), - fluid.mesh(), - dimensionedScalar - ( - IOobject::groupName("d", this->name()), - dimLength, - Zero - ) - ) + d_() { this->registerVelocityGroups(); @@ -1083,6 +1067,30 @@ Foam::diameterModels::populationBalanceModel::populationBalanceModel dimensionedScalar("Sui", inv(dimTime), Zero) ) ); + + d_.reset + ( + new volScalarField + ( + IOobject + ( + IOobject::groupName("d", this->name()), + fluid.time().timeName(), + fluid.mesh(), + IOobject::NO_READ, + IOobject::AUTO_WRITE + ), + fluid.mesh(), + dimensionedScalar + ( + IOobject::groupName("d", this->name()), + dimLength, + Zero + ) + ) + ); + + } if (coalescence_.size() != 0) @@ -1201,92 +1209,110 @@ bool Foam::diameterModels::populationBalanceModel::writeData(Ostream& os) const void Foam::diameterModels::populationBalanceModel::solve() { - const dictionary& populationBalanceControls = mesh_.solverDict(name_); - label nCorr(readLabel(populationBalanceControls.lookup("nCorr"))); - scalar tolerance(readScalar(populationBalanceControls.lookup("tolerance"))); + const dictionary& solutionControls = mesh_.solverDict(name_); + bool solveOnFinalIterOnly + ( + solutionControls.lookupOrDefault + ( + "solveOnFinalIterOnly", + false + ) + ); - calcAlphas(); - calcVelocity(); - - if (nCorr > 0) + if (!solveOnFinalIterOnly || pimple_.finalIter()) { - preSolve(); - } + calcAlphas(); + calcVelocity(); - int iCorr = 0; - scalar initialResidual = 0; - scalar maxInitialResidual = 1; + label nCorr(readLabel(solutionControls.lookup("nCorr"))); + scalar tolerance + ( + readScalar(solutionControls.lookup("tolerance")) + ); - while - ( - maxInitialResidual > tolerance - && - ++iCorr <= nCorr - ) - { - Info<< "populationBalance " - << this->name() - << ": Iteration " - << iCorr - << endl; - - sources(); - - dmdt(); - - forAll(sizeGroups_, i) + if (nCorr > 0) { - sizeGroup& fi = *sizeGroups_[i]; - const phaseModel& phase = fi.phase(); - const volScalarField& alpha = phase; - const dimensionedScalar& residualAlpha = phase.residualAlpha(); - const volScalarField& rho = phase.thermo().rho(); + preSolve(); + } - fvScalarMatrix sizeGroupEqn - ( - fvm::ddt(alpha, rho, fi) - + fi.VelocityGroup().mvConvection()->fvmDiv + int iCorr = 0; + scalar initialResidual = 0; + scalar maxInitialResidual = 1; + + while + ( + maxInitialResidual > tolerance + && + ++iCorr <= nCorr + ) + { + Info<< "populationBalance " + << this->name() + << ": Iteration " + << iCorr + << endl; + + sources(); + + dmdt(); + + forAll(sizeGroups_, i) + { + sizeGroup& fi = *sizeGroups_[i]; + const phaseModel& phase = fi.phase(); + const volScalarField& alpha = phase; + const dimensionedScalar& residualAlpha = phase.residualAlpha(); + const volScalarField& rho = phase.thermo().rho(); + + fvScalarMatrix sizeGroupEqn ( - phase.alphaRhoPhi(), - fi - ) - - fvm::Sp + fvm::ddt(alpha, rho, fi) + + fi.VelocityGroup().mvConvection()->fvmDiv + ( + phase.alphaRhoPhi(), + fi + ) + - fvm::Sp + ( + fvc::ddt(alpha, rho) + fvc::div(phase.alphaRhoPhi()) + - fi.VelocityGroup().dmdt(), + fi + ) + == + Su_[i]*rho + - fvm::SuSp(SuSp_[i]*rho, fi) + + fvc::ddt(residualAlpha*rho, fi) + - fvm::ddt(residualAlpha*rho, fi) + ); + + sizeGroupEqn.relax ( - fvc::ddt(alpha, rho) + fvc::div(phase.alphaRhoPhi()) - - fi.VelocityGroup().dmdt(), - fi - ) - == - Su_[i]*rho - - fvm::SuSp(SuSp_[i]*rho, fi) - + fvc::ddt(residualAlpha*rho, fi) - - fvm::ddt(residualAlpha*rho, fi) - ); + fi.mesh().equationRelaxationFactor("f") + ); - sizeGroupEqn.relax - ( - fi.mesh().equationRelaxationFactor("f") - ); + initialResidual = sizeGroupEqn.solve().initialResidual(); - initialResidual = sizeGroupEqn.solve().initialResidual(); + maxInitialResidual = max + ( + initialResidual, + maxInitialResidual + ); + } + } - maxInitialResidual = max - ( - initialResidual, - maxInitialResidual - ); + if (nCorr > 0) + { + forAll(velocityGroups_, i) + { + velocityGroups_[i]->postSolve(); + } + } + + if (velocityGroups_.size() > 1) + { + d_() = dsm(); } } - - if (nCorr > 0) - { - forAll(velocityGroups_, i) - { - velocityGroups_[i]->postSolve(); - } - } - - d_ = dsm(); } // ************************************************************************* // diff --git a/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/populationBalanceModel/populationBalanceModel/populationBalanceModel.H b/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/populationBalanceModel/populationBalanceModel/populationBalanceModel.H index dfab349e6..47fc910b3 100644 --- a/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/populationBalanceModel/populationBalanceModel/populationBalanceModel.H +++ b/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/populationBalanceModel/populationBalanceModel/populationBalanceModel.H @@ -136,6 +136,7 @@ SourceFiles #include "sizeGroup.H" #include "phasePair.H" +#include "pimpleControl.H" // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // @@ -184,6 +185,9 @@ class populationBalanceModel //- Dictionary dictionary dict_; + //- Reference to pimpleControl + const pimpleControl& pimple_; + //- Continuous phase const phaseModel& continuousPhase_; @@ -245,7 +249,7 @@ class populationBalanceModel volVectorField U_; //- Mean Sauter diameter - volScalarField d_; + autoPtr d_; // Private member functions diff --git a/tutorials/multiphase/reactingMultiphaseEulerFoam/laminar/bubbleColumnPolydisperse/constant/phaseProperties b/tutorials/multiphase/reactingMultiphaseEulerFoam/laminar/bubbleColumnPolydisperse/constant/phaseProperties index a45fc0ee8..3386194a3 100644 --- a/tutorials/multiphase/reactingMultiphaseEulerFoam/laminar/bubbleColumnPolydisperse/constant/phaseProperties +++ b/tutorials/multiphase/reactingMultiphaseEulerFoam/laminar/bubbleColumnPolydisperse/constant/phaseProperties @@ -122,9 +122,6 @@ populationBalanceCoeffs ( densityChange{} ); - - zeroethOrderModels - (); } } diff --git a/tutorials/multiphase/reactingTwoPhaseEulerFoam/RAS/bubbleColumnPolydisperse/constant/phaseProperties b/tutorials/multiphase/reactingTwoPhaseEulerFoam/RAS/bubbleColumnPolydisperse/constant/phaseProperties index 7628b2951..30a8e888e 100644 --- a/tutorials/multiphase/reactingTwoPhaseEulerFoam/RAS/bubbleColumnPolydisperse/constant/phaseProperties +++ b/tutorials/multiphase/reactingTwoPhaseEulerFoam/RAS/bubbleColumnPolydisperse/constant/phaseProperties @@ -104,9 +104,6 @@ populationBalanceCoeffs ( densityChange{} ); - - zeroethOrderModels - (); } } diff --git a/tutorials/multiphase/reactingTwoPhaseEulerFoam/RAS/bubbleColumnPolydisperse/system/fvSolution b/tutorials/multiphase/reactingTwoPhaseEulerFoam/RAS/bubbleColumnPolydisperse/system/fvSolution index 05221ed7f..17023600f 100644 --- a/tutorials/multiphase/reactingTwoPhaseEulerFoam/RAS/bubbleColumnPolydisperse/system/fvSolution +++ b/tutorials/multiphase/reactingTwoPhaseEulerFoam/RAS/bubbleColumnPolydisperse/system/fvSolution @@ -25,9 +25,10 @@ solvers bubbles { - nCorr 1; - tolerance 1e-4; - renormalize false; + nCorr 1; + tolerance 1e-4; + renormalize false; + solveOnFinalIterOnly true; } p_rgh