From e3903fdb358817a970d974ff1df43e6785dca7b7 Mon Sep 17 00:00:00 2001 From: Will Bainbridge Date: Thu, 27 Feb 2020 10:24:32 +0000 Subject: [PATCH] reacting*EulerFoam: Separate continuity error from mass transfer continuityError is now just the transport inconsistency. Mass sources, whether as a result of fvOptions or phase-transfer/change processes, are not included. --- .../MomentumTransferPhaseSystem.C | 24 +++-------- .../OneResistanceHeatTransferPhaseSystem.C | 26 ++++-------- .../PhaseTransferPhaseSystem.C | 8 +--- .../ThermalPhaseChangePhaseSystem.C | 3 -- .../TwoResistanceHeatTransferPhaseSystem.C | 42 ++++++++----------- .../sizeGroup/shapeModels/fractal/fractal.C | 17 -------- .../MovingPhaseModel/MovingPhaseModel.C | 17 +++----- .../MovingPhaseModel/MovingPhaseModel.H | 9 ++-- .../phaseModel/phaseModel/phaseModel.C | 6 ++- .../phaseModel/phaseModel/phaseModel.H | 3 ++ .../phaseSystems/phaseSystem/phaseSystem.C | 35 ++++++++++++++++ .../phaseSystems/phaseSystem/phaseSystem.H | 8 ++-- .../populationBalanceModel.C | 16 +------ .../reactingMultiphaseEulerFoam/EEqns.H | 1 + .../reactingMultiphaseEulerFoam.C | 1 + .../reactingTwoPhaseEulerFoam/EEqns.H | 1 + .../reactingTwoPhaseEulerFoam.C | 1 + .../twoPhaseSystem/diameterModels/IATE/IATE.C | 12 ++++-- 18 files changed, 101 insertions(+), 129 deletions(-) diff --git a/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/PhaseSystems/MomentumTransferPhaseSystem/MomentumTransferPhaseSystem.C b/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/PhaseSystems/MomentumTransferPhaseSystem/MomentumTransferPhaseSystem.C index 9df985ec07..3d1a9a40ff 100644 --- a/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/PhaseSystems/MomentumTransferPhaseSystem/MomentumTransferPhaseSystem.C +++ b/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/PhaseSystems/MomentumTransferPhaseSystem/MomentumTransferPhaseSystem.C @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | Website: https://openfoam.org - \\ / A nd | Copyright (C) 2015-2019 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2015-2020 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -135,33 +135,19 @@ void Foam::MomentumTransferPhaseSystem::addDmdtUfs phaseModel& phase1 = this->phases()[pair.phase1().name()]; phaseModel& phase2 = this->phases()[pair.phase2().name()]; - // Note that the phase UEqn contains a continuity error term, which - // implicitly adds a mass transfer term of fvm::Sp(dmdt, U). These - // additions remove the part of this term corresponding to the supplied - // mass transfer rate; I.e.: - // - // DDt(alpha1*rho1*U1) + ... = dmdt21*U2 + dmdt12*U1 - // DDt(alpha1*rho1*U1) - contErr + ... = dmdt21*U2 + dmdt12*U1 - contErr - // DDt(alpha1*rho1*U1) - contErr + ... = dmdt21*U2 + dmdt12*U1 - dmdt*U1 - // DDt(alpha1*rho1*U1) - contErr + ... = dmdt21*U2 - dmdt21*U1 - // - // Where contErr is the continuity error associated with this mass - // transfer rate, dmdt. + const volScalarField dmdtf21(posPart(dmdtf)); + const volScalarField dmdtf12(negPart(dmdtf)); if (!phase1.stationary()) { - const volScalarField dmdtf21(posPart(dmdtf)); - *eqns[phase1.name()] += - dmdtf21*phase2.U() - fvm::Sp(dmdtf21, phase1.URef()); + dmdtf21*phase2.U() + fvm::Sp(dmdtf12, phase1.URef()); } if (!phase2.stationary()) { - const volScalarField dmdtf12(negPart(dmdtf)); - *eqns[phase2.name()] -= - dmdtf12*phase1.U() - fvm::Sp(dmdtf12, phase2.URef()); + dmdtf12*phase1.U() + fvm::Sp(dmdtf21, phase2.URef()); } } } diff --git a/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/PhaseSystems/OneResistanceHeatTransferPhaseSystem/OneResistanceHeatTransferPhaseSystem.C b/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/PhaseSystems/OneResistanceHeatTransferPhaseSystem/OneResistanceHeatTransferPhaseSystem.C index 7828a3b963..d80a0c9077 100644 --- a/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/PhaseSystems/OneResistanceHeatTransferPhaseSystem/OneResistanceHeatTransferPhaseSystem.C +++ b/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/PhaseSystems/OneResistanceHeatTransferPhaseSystem/OneResistanceHeatTransferPhaseSystem.C @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | Website: https://openfoam.org - \\ / A nd | Copyright (C) 2015-2019 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2015-2020 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -55,17 +55,13 @@ void Foam::OneResistanceHeatTransferPhaseSystem::addDmdtHefs const volScalarField K1(phase1.K()); const volScalarField K2(phase2.K()); - // Note that the phase EEqn contains a continuity error term. See - // MomentumTransferPhaseSystem::addDmdtUfs for an explanation of the - // fvm::Sp terms below. - // Transfer of energy from bulk to bulk - *eqns[phase1.name()] += dmdtf21*he2 - fvm::Sp(dmdtf21, he1); - *eqns[phase2.name()] -= dmdtf12*he1 - fvm::Sp(dmdtf12, he2); + *eqns[phase1.name()] += dmdtf21*he2 + fvm::Sp(dmdtf12, he1); + *eqns[phase2.name()] -= dmdtf12*he1 + fvm::Sp(dmdtf21, he2); // Transfer of kinetic energy - *eqns[phase1.name()] += dmdtf21*(K2 - K1); - *eqns[phase2.name()] -= dmdtf12*(K1 - K2); + *eqns[phase1.name()] += dmdtf21*K2 + dmdtf12*K1; + *eqns[phase2.name()] -= dmdtf12*K1 + dmdtf21*K2; } } @@ -91,10 +87,6 @@ void Foam::OneResistanceHeatTransferPhaseSystem::addDmidtHef const volScalarField K1(phase1.K()); const volScalarField K2(phase2.K()); - // Note that the phase EEqn contains a continuity error term. See - // MomentumTransferPhaseSystem::addDmdtUfs for an explanation of the - // fvm::Sp terms below. - forAllConstIter(HashPtrTable, *dmidtfIter(), dmidtfJter) { const word& member = dmidtfJter.key(); @@ -135,12 +127,12 @@ void Foam::OneResistanceHeatTransferPhaseSystem::addDmidtHef } // Transfer of energy from bulk to bulk - *eqns[phase1.name()] += dmidtf21*hei2 - fvm::Sp(dmidtf21, he1); - *eqns[phase2.name()] -= dmidtf12*hei1 - fvm::Sp(dmidtf12, he2); + *eqns[phase1.name()] += dmidtf21*hei2 + fvm::Sp(dmidtf12, he1); + *eqns[phase2.name()] -= dmidtf12*hei1 + fvm::Sp(dmidtf21, he2); // Transfer of kinetic energy - *eqns[phase1.name()] += dmidtf21*(K2 - K1); - *eqns[phase2.name()] -= dmidtf12*(K1 - K2); + *eqns[phase1.name()] += dmidtf21*K2 + dmidtf12*K1; + *eqns[phase2.name()] -= dmidtf12*K1 + dmidtf21*K2; } } } diff --git a/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/PhaseSystems/PhaseTransferPhaseSystem/PhaseTransferPhaseSystem.C b/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/PhaseSystems/PhaseTransferPhaseSystem/PhaseTransferPhaseSystem.C index 922f278a83..a63138d249 100644 --- a/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/PhaseSystems/PhaseTransferPhaseSystem/PhaseTransferPhaseSystem.C +++ b/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/PhaseSystems/PhaseTransferPhaseSystem/PhaseTransferPhaseSystem.C @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | Website: https://openfoam.org - \\ / A nd | Copyright (C) 2018-2019 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2018-2020 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -92,9 +92,6 @@ void Foam::PhaseTransferPhaseSystem::addDmdtYfs const phaseModel& phase1 = pair.phase1(); const phaseModel& phase2 = pair.phase2(); - // Note that the phase YiEqn does not contain a continuity error term, - // so the transfers below are complete. - forAll(phase1.Y(), Yi1) { const volScalarField& Y1 = phase1.Y()[Yi1]; @@ -122,9 +119,6 @@ void Foam::PhaseTransferPhaseSystem::addDmidtYf const phaseModel& phase1 = pair.phase1(); const phaseModel& phase2 = pair.phase2(); - // Note that the phase YiEqn does not contain a continuity error term, - // so the transfers below are complete. - forAllConstIter(HashPtrTable, *dmidtfIter(), dmidtfJter) { const word& member = dmidtfJter.key(); diff --git a/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/PhaseSystems/ThermalPhaseChangePhaseSystem/ThermalPhaseChangePhaseSystem.C b/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/PhaseSystems/ThermalPhaseChangePhaseSystem/ThermalPhaseChangePhaseSystem.C index 75b8d93375..21e52385e9 100644 --- a/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/PhaseSystems/ThermalPhaseChangePhaseSystem/ThermalPhaseChangePhaseSystem.C +++ b/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/PhaseSystems/ThermalPhaseChangePhaseSystem/ThermalPhaseChangePhaseSystem.C @@ -447,9 +447,6 @@ Foam::ThermalPhaseChangePhaseSystem::specieTransfer() const const volScalarField& Y1 = phase1.Y(volatile_); const volScalarField& Y2 = phase2.Y(volatile_); - // Note that the phase YiEqn does not contain a continuity error - // term, so these additions represent the entire mass transfer - const volScalarField dmdtf(this->totalDmdtf(pair)); *eqns[Y1.name()] += dmdtf; diff --git a/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/PhaseSystems/TwoResistanceHeatTransferPhaseSystem/TwoResistanceHeatTransferPhaseSystem.C b/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/PhaseSystems/TwoResistanceHeatTransferPhaseSystem/TwoResistanceHeatTransferPhaseSystem.C index 62351536bb..8095f689de 100644 --- a/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/PhaseSystems/TwoResistanceHeatTransferPhaseSystem/TwoResistanceHeatTransferPhaseSystem.C +++ b/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/PhaseSystems/TwoResistanceHeatTransferPhaseSystem/TwoResistanceHeatTransferPhaseSystem.C @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | Website: https://openfoam.org - \\ / A nd | Copyright (C) 2015-2019 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2015-2020 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -80,10 +80,6 @@ void Foam::TwoResistanceHeatTransferPhaseSystem::addDmdtHefs const volScalarField K1(phase1.K()); const volScalarField K2(phase2.K()); - // Note that the phase EEqn contains a continuity error term. See - // MomentumTransferPhaseSystem::addDmdtU for an explanation of the - // fvm::Sp terms below. - if (heatTransferModels_.found(key)) { // Assume a thermally-coupled mass transfer process. Calculate @@ -94,8 +90,8 @@ void Foam::TwoResistanceHeatTransferPhaseSystem::addDmdtHefs const volScalarField& Tf(*Tf_[key]); const volScalarField hef1(thermo1.he(thermo1.p(), Tf)); const volScalarField hef2(thermo2.he(thermo2.p(), Tf)); - *eqns[phase1.name()] += dmdtf*hef1 - fvm::Sp(dmdtf, he1); - *eqns[phase2.name()] -= dmdtf*hef2 - fvm::Sp(dmdtf, he2); + *eqns[phase1.name()] += dmdtf*hef1; + *eqns[phase2.name()] -= dmdtf*hef2; // Latent heat contribution const volScalarField L(hef2 + thermo2.hc() - hef1 - thermo1.hc()); @@ -106,8 +102,8 @@ void Foam::TwoResistanceHeatTransferPhaseSystem::addDmdtHefs *eqns[phase2.name()] += (1 - H1Fac)*dmdtf*L; // Transfer of kinetic energy - *eqns[phase1.name()] += dmdtf21*(K2 - K1); - *eqns[phase2.name()] -= dmdtf12*(K1 - K2); + *eqns[phase1.name()] += dmdtf21*K2 + dmdtf12*K1; + *eqns[phase2.name()] -= dmdtf12*K1 + dmdtf21*K2; } else { @@ -116,12 +112,12 @@ void Foam::TwoResistanceHeatTransferPhaseSystem::addDmdtHefs // change and therefore no interface state or latent heat... // Transfer of energy from bulk to bulk - *eqns[phase1.name()] += dmdtf21*he2 - fvm::Sp(dmdtf21, he1); - *eqns[phase2.name()] -= dmdtf12*he1 - fvm::Sp(dmdtf12, he2); + *eqns[phase1.name()] += dmdtf21*he2 + fvm::Sp(dmdtf12, he1); + *eqns[phase2.name()] -= dmdtf12*he1 + fvm::Sp(dmdtf21, he2); // Transfer of kinetic energy - *eqns[phase1.name()] += dmdtf21*(K2 - K1); - *eqns[phase2.name()] -= dmdtf12*(K1 - K2); + *eqns[phase1.name()] += dmdtf21*K2 + dmdtf12*K1; + *eqns[phase2.name()] -= dmdtf12*K1 + dmdtf21*K2; } } } @@ -150,10 +146,6 @@ void Foam::TwoResistanceHeatTransferPhaseSystem::addDmidtHef const volScalarField K1(phase1.K()); const volScalarField K2(phase2.K()); - // Note that the phase EEqn contains a continuity error term. See - // MomentumTransferPhaseSystem::addDmdtU for an explanation of the - // fvm::Sp terms below. - if (heatTransferModels_.found(key)) { // Assume a thermally-coupled mass transfer process. Calculate @@ -233,16 +225,16 @@ void Foam::TwoResistanceHeatTransferPhaseSystem::addDmidtHef const volScalarField Li(hefi2 + hci2 - hefi1 - hci1); // Transfer of energy from the interface into the bulk - *eqns[phase1.name()] += dmidtf*hefi1 - fvm::Sp(dmidtf, he1); - *eqns[phase2.name()] -= dmidtf*hefi2 - fvm::Sp(dmidtf, he2); + *eqns[phase1.name()] += dmidtf*hefi1; + *eqns[phase2.name()] -= dmidtf*hefi2; // Latent heat contribution *eqns[phase1.name()] += H1Fac*dmidtf*Li; *eqns[phase2.name()] += (1 - H1Fac)*dmidtf*Li; // Transfer of kinetic energy - *eqns[phase1.name()] += dmidtf21*(K2 - K1); - *eqns[phase2.name()] -= dmidtf12*(K1 - K2); + *eqns[phase1.name()] += dmidtf21*K2 + dmidtf12*K1; + *eqns[phase2.name()] -= dmidtf12*K1 + dmidtf21*K2; } } else @@ -297,12 +289,12 @@ void Foam::TwoResistanceHeatTransferPhaseSystem::addDmidtHef } // Transfer of energy from bulk to bulk - *eqns[phase1.name()] += dmidtf21*hei2 - fvm::Sp(dmidtf21, he1); - *eqns[phase2.name()] -= dmidtf12*hei1 - fvm::Sp(dmidtf12, he2); + *eqns[phase1.name()] += dmidtf21*hei2 + fvm::Sp(dmidtf12, he1); + *eqns[phase2.name()] -= dmidtf12*hei1 + fvm::Sp(dmidtf21, he2); // Transfer of kinetic energy - *eqns[phase1.name()] += dmidtf21*(K2 - K1); - *eqns[phase2.name()] -= dmidtf12*(K1 - K2); + *eqns[phase1.name()] += dmidtf21*K2 + dmidtf12*K1; + *eqns[phase2.name()] -= dmidtf12*K1 + dmidtf21*K2; } } } diff --git a/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/diameterModels/velocityGroup/sizeGroup/shapeModels/fractal/fractal.C b/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/diameterModels/velocityGroup/sizeGroup/shapeModels/fractal/fractal.C index 3b971c0e06..eaa5ef2301 100644 --- a/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/diameterModels/velocityGroup/sizeGroup/shapeModels/fractal/fractal.C +++ b/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/diameterModels/velocityGroup/sizeGroup/shapeModels/fractal/fractal.C @@ -215,23 +215,6 @@ void Foam::diameterModels::shapeModels::fractal::correct() + alpha*rho*fi*fvm::ddt(kappa_) + fvm::div(fAlphaRhoPhi, kappa_) + fvm::SuSp(-phase.continuityError()*fi, kappa_) - + fvm::SuSp - ( - ( - fi.VelocityGroup().dmdt() - - // Temporary term pending update of continuityError - - ( - phase.fluid().fvOptions() - ( - alpha, - const_cast(rho) - ) & rho - ) - ) - *fi, - kappa_ - ) == - sinteringModel_->R() + fvc::Su(Su_*rho, kappa_) diff --git a/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/phaseModel/MovingPhaseModel/MovingPhaseModel.C b/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/phaseModel/MovingPhaseModel/MovingPhaseModel.C index edc13703b7..7bb6c95fa5 100644 --- a/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/phaseModel/MovingPhaseModel/MovingPhaseModel.C +++ b/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/phaseModel/MovingPhaseModel/MovingPhaseModel.C @@ -203,13 +203,14 @@ Foam::MovingPhaseModel::~MovingPhaseModel() // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // template -void Foam::MovingPhaseModel::correctContinuityError() +void Foam::MovingPhaseModel::correctContinuityError +( + const volScalarField& source +) { volScalarField& rho = this->thermoRef().rho(); - continuityError_ = - fvc::ddt(*this, rho) + fvc::div(alphaRhoPhi_) - - (this->fluid().fvOptions()(*this, rho)&rho); + continuityError_ = fvc::ddt(*this, rho) + fvc::div(alphaRhoPhi_) - source; } @@ -218,7 +219,6 @@ void Foam::MovingPhaseModel::correct() { BasePhaseModel::correct(); this->fluid().MRF().correctBoundaryVelocity(U_); - correctContinuityError(); } @@ -247,13 +247,6 @@ void Foam::MovingPhaseModel::correctKinematics() } -template -void Foam::MovingPhaseModel::correctThermo() -{ - correctContinuityError(); -} - - template void Foam::MovingPhaseModel::correctTurbulence() { diff --git a/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/phaseModel/MovingPhaseModel/MovingPhaseModel.H b/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/phaseModel/MovingPhaseModel/MovingPhaseModel.H index d74a72ecb6..b5ec8f9a6b 100644 --- a/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/phaseModel/MovingPhaseModel/MovingPhaseModel.H +++ b/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/phaseModel/MovingPhaseModel/MovingPhaseModel.H @@ -103,9 +103,6 @@ private: //- Calculate and return the flux field tmp phi(const volVectorField& U) const; - //- Correct the continuity error - void correctContinuityError(); - public: @@ -128,12 +125,12 @@ public: //- Correct the phase properties other than the thermo and turbulence virtual void correct(); + //- Correct the continuity error + virtual void correctContinuityError(const volScalarField& source); + //- Correct the kinematics virtual void correctKinematics(); - //- Correct the thermodynamics - virtual void correctThermo(); - //- Correct the turbulence virtual void correctTurbulence(); diff --git a/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/phaseModel/phaseModel/phaseModel.C b/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/phaseModel/phaseModel/phaseModel.C index 5ab92a8986..ddac7038c4 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 | Website: https://openfoam.org - \\ / A nd | Copyright (C) 2015-2019 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2015-2020 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -143,6 +143,10 @@ void Foam::phaseModel::correct() } +void Foam::phaseModel::correctContinuityError(const volScalarField& source) +{} + + void Foam::phaseModel::correctKinematics() {} diff --git a/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/phaseModel/phaseModel/phaseModel.H b/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/phaseModel/phaseModel/phaseModel.H index e676fabc8b..2fd4611e39 100644 --- a/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/phaseModel/phaseModel/phaseModel.H +++ b/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/phaseModel/phaseModel/phaseModel.H @@ -186,6 +186,9 @@ public: //- Correct the phase properties virtual void correct(); + //- Correct the continuity error + virtual void correctContinuityError(const volScalarField& source); + //- Correct the kinematics virtual void correctKinematics(); diff --git a/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/phaseSystem/phaseSystem.C b/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/phaseSystem/phaseSystem.C index 7df7108328..30d8642f9e 100644 --- a/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/phaseSystem/phaseSystem.C +++ b/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/phaseSystem/phaseSystem.C @@ -450,6 +450,41 @@ void Foam::phaseSystem::correct() } +void Foam::phaseSystem::correctContinuityError() +{ + const PtrList dmdts = this->dmdts(); + + forAll(movingPhaseModels_, movingPhasei) + { + phaseModel& phase = movingPhaseModels_[movingPhasei]; + const volScalarField& alpha = phase; + volScalarField& rho = phase.thermoRef().rho(); + + volScalarField source + ( + volScalarField::New + ( + IOobject::groupName("source", phase.name()), + mesh_, + dimensionedScalar(dimDensity/dimTime, 0) + ) + ); + + if (fvOptions().appliesToField(rho.name())) + { + source += fvOptions()(alpha, rho)ρ + } + + if (dmdts.set(phase.index())) + { + source += dmdts[phase.index()]; + } + + phase.correctContinuityError(source); + } +} + + void Foam::phaseSystem::correctKinematics() { bool updateDpdt = false; diff --git a/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/phaseSystem/phaseSystem.H b/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/phaseSystem/phaseSystem.H index a90e067ffe..13a6ae96ba 100644 --- a/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/phaseSystem/phaseSystem.H +++ b/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/phaseSystem/phaseSystem.H @@ -187,10 +187,7 @@ protected: ) const; //- Generate pairs - void generatePairs - ( - const dictTable& modelDicts - ); + void generatePairs(const dictTable& modelDicts); //- Generate pairs and sub-model tables template @@ -523,6 +520,9 @@ public: //- Correct the fluid properties other than those listed below virtual void correct(); + //- Correct the continuity errors + virtual void correctContinuityError(); + //- Correct the kinematics virtual void correctKinematics(); diff --git a/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/populationBalanceModel/populationBalanceModel/populationBalanceModel.C b/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/populationBalanceModel/populationBalanceModel/populationBalanceModel.C index 0c58d3c13d..dd176dafaa 100644 --- a/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/populationBalanceModel/populationBalanceModel/populationBalanceModel.C +++ b/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/populationBalanceModel/populationBalanceModel/populationBalanceModel.C @@ -1304,21 +1304,7 @@ void Foam::diameterModels::populationBalanceModel::solve() phase.alphaRhoPhi(), fi ) - + fvm::SuSp - ( - fi.VelocityGroup().dmdt() - - phase.continuityError() - - // Temporary term pending update of continuityError - - ( - phase.fluid().fvOptions() - ( - alpha, - const_cast(rho) - ) & rho - ), - fi - ) + + fvm::SuSp(-phase.continuityError(), fi) == fvc::Su(Su_[i]*rho, fi) - fvm::SuSp(SuSp_[i]*rho, fi) diff --git a/applications/solvers/multiphase/reactingEulerFoam/reactingMultiphaseEulerFoam/EEqns.H b/applications/solvers/multiphase/reactingEulerFoam/reactingMultiphaseEulerFoam/EEqns.H index 6e8d5de9a0..dbaf73a7b4 100644 --- a/applications/solvers/multiphase/reactingEulerFoam/reactingMultiphaseEulerFoam/EEqns.H +++ b/applications/solvers/multiphase/reactingEulerFoam/reactingMultiphaseEulerFoam/EEqns.H @@ -33,6 +33,7 @@ for (int Ecorr=0; Ecorr