diff --git a/applications/solvers/multiphase/icoReactingMultiphaseInterFoam/TEqn.H b/applications/solvers/multiphase/icoReactingMultiphaseInterFoam/TEqn.H index 68744f7043b31aab583ad837596d7b466864a12d..c0405456173739cc1bbfe9e3f80a46cc89983b9d 100644 --- a/applications/solvers/multiphase/icoReactingMultiphaseInterFoam/TEqn.H +++ b/applications/solvers/multiphase/icoReactingMultiphaseInterFoam/TEqn.H @@ -1,7 +1,7 @@ { radiation->correct(); rhoCp = rho*fluid.Cp(); - + const surfaceScalarField rhoCpPhi(fvc::interpolate(fluid.Cp())*rhoPhi); const volScalarField kappaEff @@ -31,7 +31,7 @@ fvOptions.correct(T); fluid.correct(); - + Info<< "min/max(T) = " << min(T).value() << ", " << max(T).value() << endl; } diff --git a/applications/solvers/multiphase/icoReactingMultiphaseInterFoam/UEqn.H b/applications/solvers/multiphase/icoReactingMultiphaseInterFoam/UEqn.H index 4422f0bd578fa506395f09f037b26e48c320dae4..9d0a424d265403ac1010b61e25e8d8470f94e2fd 100644 --- a/applications/solvers/multiphase/icoReactingMultiphaseInterFoam/UEqn.H +++ b/applications/solvers/multiphase/icoReactingMultiphaseInterFoam/UEqn.H @@ -2,7 +2,7 @@ fvVectorMatrix UEqn ( fvm::ddt(rho, U) + fvm::div(rhoPhi, U) - //- fvm::Sp(fvc::ddt(rho) + fvc::div(rhoPhi), U) +//- fvm::Sp(fvc::ddt(rho) + fvc::div(rhoPhi), U) + turbulence->divDevRhoReff(U) == fvOptions(rho, U) diff --git a/applications/solvers/multiphase/icoReactingMultiphaseInterFoam/icoReactingMultiphaseInterFoam.C b/applications/solvers/multiphase/icoReactingMultiphaseInterFoam/icoReactingMultiphaseInterFoam.C index 0c42797f4ca493d098a05aa4fc01233de095215a..9bbe933d8c0b7d8e830235913fa2e9989705a3ef 100644 --- a/applications/solvers/multiphase/icoReactingMultiphaseInterFoam/icoReactingMultiphaseInterFoam.C +++ b/applications/solvers/multiphase/icoReactingMultiphaseInterFoam/icoReactingMultiphaseInterFoam.C @@ -5,7 +5,7 @@ \\ / A nd | www.openfoam.com \\/ M anipulation | ------------------------------------------------------------------------------- - Copyright (C) 2017 OpenCFD Ltd. + Copyright (C) 2017-2020 OpenCFD Ltd. ------------------------------------------------------------------------------- License This file is part of OpenFOAM. @@ -98,13 +98,13 @@ int main(int argc, char *argv[]) Info<< "Time = " << runTime.timeName() << nl << endl; + fluid.correctMassSources(T); fluid.solve(); rho = fluid.rho(); // --- Pressure-velocity PIMPLE corrector loop while (pimple.loop()) { - solve(fvm::ddt(rho) + fvc::div(rhoPhi)); #include "UEqn.H" #include "YEqns.H" #include "TEqn.H" diff --git a/applications/solvers/multiphase/icoReactingMultiphaseInterFoam/massTransferModels/Lee/Lee.C b/applications/solvers/multiphase/icoReactingMultiphaseInterFoam/massTransferModels/Lee/Lee.C index 7588b83c9dccb1860a5822445e7a2242a2ec4cd6..190b7fbb1d04fb74bcec85cbb289d1e1b7f8dd3f 100644 --- a/applications/solvers/multiphase/icoReactingMultiphaseInterFoam/massTransferModels/Lee/Lee.C +++ b/applications/solvers/multiphase/icoReactingMultiphaseInterFoam/massTransferModels/Lee/Lee.C @@ -50,6 +50,44 @@ Foam::meltingEvaporationModels::Lee::Lee template Foam::tmp Foam::meltingEvaporationModels::Lee::Kexp +( + const volScalarField& refValue +) +{ + { + const volScalarField from + ( + min(max(this->pair().from(), scalar(0)), scalar(1)) + ); + + const volScalarField coeff + ( + C_*from*this->pair().from().rho()*pos(from - alphaMin_) + *(refValue - Tactivate_) + /Tactivate_ + ); + + if (sign(C_.value()) > 0) + { + return + ( + coeff*pos(refValue - Tactivate_) + ); + } + else + { + return + ( + coeff*pos(Tactivate_ - refValue) + ); + } + } +} + + +template +Foam::tmp +Foam::meltingEvaporationModels::Lee::KSp ( label variable, const volScalarField& refValue @@ -61,31 +99,68 @@ Foam::meltingEvaporationModels::Lee::Kexp ( min(max(this->pair().from(), scalar(0)), scalar(1)) ); + + const volScalarField coeff + ( + C_*from*this->pair().from().rho()*pos(from - alphaMin_) + /Tactivate_ + ); if (sign(C_.value()) > 0) { return ( - C_ - * from - * this->pair().from().rho() - * (refValue.oldTime() - Tactivate_) - * pos(from - alphaMin_) - * pos(refValue.oldTime() - Tactivate_)/Tactivate_ + coeff*pos(refValue - Tactivate_) ); } else { return ( - -C_ - * from - * this->pair().from().rho() - * pos(from - alphaMin_) - * (Tactivate_ - refValue.oldTime()) - * pos(Tactivate_ - refValue.oldTime())/Tactivate_ + coeff*pos(Tactivate_ - refValue) ); + } + } + else + { + return tmp (); + } +} + +template +Foam::tmp +Foam::meltingEvaporationModels::Lee::KSu +( + label variable, + const volScalarField& refValue +) +{ + if (this->modelVariable_ == variable) + { + volScalarField from + ( + min(max(this->pair().from(), scalar(0)), scalar(1)) + ); + + const volScalarField coeff + ( + C_*from*this->pair().from().rho()*pos(from - alphaMin_) + ); + + if (sign(C_.value()) > 0) + { + return + ( + -coeff*pos(refValue - Tactivate_) + ); + } + else + { + return + ( + coeff*pos(Tactivate_ - refValue) + ); } } else @@ -103,4 +178,12 @@ Foam::meltingEvaporationModels::Lee::Tactivate() const } +template +bool +Foam::meltingEvaporationModels::Lee::includeDivU() +{ + return true; +} + + // ************************************************************************* // diff --git a/applications/solvers/multiphase/icoReactingMultiphaseInterFoam/massTransferModels/Lee/Lee.H b/applications/solvers/multiphase/icoReactingMultiphaseInterFoam/massTransferModels/Lee/Lee.H index 9743ddcfb2089987d9d8cfe03d04b908a170e085..97059d4f29981bcbc8e35682561d402f12672134 100644 --- a/applications/solvers/multiphase/icoReactingMultiphaseInterFoam/massTransferModels/Lee/Lee.H +++ b/applications/solvers/multiphase/icoReactingMultiphaseInterFoam/massTransferModels/Lee/Lee.H @@ -146,12 +146,29 @@ public: //- Explicit mass transfer coefficient virtual tmp Kexp ( - label variable, + const volScalarField& field + ); + + //- Implicit mass transfer coefficient + virtual tmp KSp + ( + label modelVariable, + const volScalarField& field + ); + + //- Explicit mass transfer coefficient + virtual tmp KSu + ( + label modelVariable, const volScalarField& field ); //- Return T transition between phases virtual const dimensionedScalar& Tactivate() const; + + //- Adds and substract alpha*div(U) as a source term + // for alpha, substituting div(U) = mDot(1/rho1 - 1/rho2) + virtual bool includeDivU(); }; diff --git a/applications/solvers/multiphase/icoReactingMultiphaseInterFoam/massTransferModels/interfaceCompositionModel/interfaceCompositionModel.C b/applications/solvers/multiphase/icoReactingMultiphaseInterFoam/massTransferModels/interfaceCompositionModel/interfaceCompositionModel.C index 442b484bcd453706fb215c60b9b99329e02b121d..8d54a861079e5b839c62d3c0f1049f456414f831 100644 --- a/applications/solvers/multiphase/icoReactingMultiphaseInterFoam/massTransferModels/interfaceCompositionModel/interfaceCompositionModel.C +++ b/applications/solvers/multiphase/icoReactingMultiphaseInterFoam/massTransferModels/interfaceCompositionModel/interfaceCompositionModel.C @@ -90,4 +90,10 @@ const Foam::word Foam::interfaceCompositionModel::variable() const } +bool Foam::interfaceCompositionModel::includeDivU() +{ + return true; +} + + // ************************************************************************* // diff --git a/applications/solvers/multiphase/icoReactingMultiphaseInterFoam/massTransferModels/interfaceCompositionModel/interfaceCompositionModel.H b/applications/solvers/multiphase/icoReactingMultiphaseInterFoam/massTransferModels/interfaceCompositionModel/interfaceCompositionModel.H index a4ccf813a232d885c5c0f50def13d06e5671047e..498290199b6cec26173a2fc1f8c4a7bb0a803ac9 100644 --- a/applications/solvers/multiphase/icoReactingMultiphaseInterFoam/massTransferModels/interfaceCompositionModel/interfaceCompositionModel.H +++ b/applications/solvers/multiphase/icoReactingMultiphaseInterFoam/massTransferModels/interfaceCompositionModel/interfaceCompositionModel.H @@ -67,7 +67,8 @@ public: { T, /* temperature based */ P, /* pressure based */ - Y /* mass fraction based */ + Y, /* mass fraction based */ + alpha /* alpha source */ }; static const Enum modelVariableNames; @@ -170,15 +171,33 @@ public: const volScalarField& Tf ) const = 0; - //- Explicit mass transfer coefficient + //- Explicit full mass transfer virtual tmp Kexp + ( + const volScalarField& field + ) = 0; + + //- Implicit mass transfer + virtual tmp KSp + ( + label modelVariable, + const volScalarField& field + ) = 0; + + //- Explicit mass transfer + virtual tmp KSu ( label modelVariable, const volScalarField& field ) = 0; + //- Reference value virtual const dimensionedScalar& Tactivate() const = 0; + + //- Adds and substract alpha*div(U) as a source term + // for alpha, substituting div(U) = mDot(1/rho1 - 1/rho2) + virtual bool includeDivU(); //- Returns the variable on which the model is based const word variable() const; diff --git a/applications/solvers/multiphase/icoReactingMultiphaseInterFoam/massTransferModels/kineticGasEvaporation/kineticGasEvaporation.C b/applications/solvers/multiphase/icoReactingMultiphaseInterFoam/massTransferModels/kineticGasEvaporation/kineticGasEvaporation.C index 53a398f44ff194268774e5ae386d47d40852b13a..3c544fdaeaa6923e8f7f79ac80ccab298d7fc1f1 100644 --- a/applications/solvers/multiphase/icoReactingMultiphaseInterFoam/massTransferModels/kineticGasEvaporation/kineticGasEvaporation.C +++ b/applications/solvers/multiphase/icoReactingMultiphaseInterFoam/massTransferModels/kineticGasEvaporation/kineticGasEvaporation.C @@ -88,9 +88,8 @@ Foam::meltingEvaporationModels::kineticGasEvaporation template Foam::tmp Foam::meltingEvaporationModels::kineticGasEvaporation -::Kexp(label variable, const volScalarField& field) +::Kexp(const volScalarField& field) { - if (this->modelVariable_ == variable) { const volScalarField& to = this->pair().to(); @@ -246,10 +245,29 @@ Foam::meltingEvaporationModels::kineticGasEvaporation return massFluxEvap*areaDensity*Nl*from; } - else - { - return tmp (); - } +} + +template +Foam::tmp +Foam::meltingEvaporationModels::kineticGasEvaporation::KSu +( + label variable, + const volScalarField& refValue +) +{ + return tmp (); +} + + +template +Foam::tmp +Foam::meltingEvaporationModels::kineticGasEvaporation::KSp +( + label variable, + const volScalarField& refValue +) +{ + return tmp (); } diff --git a/applications/solvers/multiphase/icoReactingMultiphaseInterFoam/massTransferModels/kineticGasEvaporation/kineticGasEvaporation.H b/applications/solvers/multiphase/icoReactingMultiphaseInterFoam/massTransferModels/kineticGasEvaporation/kineticGasEvaporation.H index 8f6fbace00b5d506532edd6bbd3a5466af92e1bf..4c85bf2e03ba66381b475386aa6fc4ef52163ea5 100644 --- a/applications/solvers/multiphase/icoReactingMultiphaseInterFoam/massTransferModels/kineticGasEvaporation/kineticGasEvaporation.H +++ b/applications/solvers/multiphase/icoReactingMultiphaseInterFoam/massTransferModels/kineticGasEvaporation/kineticGasEvaporation.H @@ -185,7 +185,20 @@ public: //- Explicit mass transfer coefficient virtual tmp Kexp ( - label variable, + const volScalarField& field + ); + + //- Implicit mass transfer coefficient + virtual tmp KSp + ( + label modelVariable, + const volScalarField& field + ); + + //- Explicit mass transfer coefficient + virtual tmp KSu + ( + label modelVariable, const volScalarField& field ); diff --git a/applications/solvers/multiphase/icoReactingMultiphaseInterFoam/pEqn.H b/applications/solvers/multiphase/icoReactingMultiphaseInterFoam/pEqn.H index 5dd79906202ee39662b03480b807caab48d77f30..463f59da9e47bb3962ba7a46c2acb58de94bc8d8 100644 --- a/applications/solvers/multiphase/icoReactingMultiphaseInterFoam/pEqn.H +++ b/applications/solvers/multiphase/icoReactingMultiphaseInterFoam/pEqn.H @@ -33,47 +33,9 @@ ( fvc::div(phiHbyA) - fvm::laplacian(rAUf, p_rgh) + + fluid.volTransfer(p_rgh) ); - - forAllConstIters(fluid.totalPhasePairs(), iter) - { - const phasePair& pair = iter()(); - - const phaseModel& phase1 = pair.phase1(); - const phaseModel& phase2 = pair.phase2(); - - const phasePairKey key12 - ( - phase1.name(), - phase2.name(), - true - ); - - // Mass transfer from phase2 to phase1 - tmp tdmdt12(fluid.dmdt(key12)); - const volScalarField& dmdt12 = tdmdt12(); - - const phasePairKey key21 - ( - phase2.name(), - phase1.name(), - true - ); - - // Mass transfer from phase1 to phase2 - tmp tdmdt21(fluid.dmdt(key21)); - const volScalarField& dmdt21 = tdmdt21(); - - const volScalarField dmdtNet(dmdt21 - dmdt12); - - p_rghEqn += - dmdtNet* - ( - - fluid.coeffs(phase1.name()) - + fluid.coeffs(phase2.name()) - ); - } - + p_rghEqn.setReference(pRefCell, pRefValue); p_rghEqn.solve(mesh.solver(p_rgh.select(pimple.finalInnerIter()))); diff --git a/applications/solvers/multiphase/icoReactingMultiphaseInterFoam/phasesSystem/MassTransferPhaseSystem/MassTransferPhaseSystem.C b/applications/solvers/multiphase/icoReactingMultiphaseInterFoam/phasesSystem/MassTransferPhaseSystem/MassTransferPhaseSystem.C index 59b04864798300fddf2d18daffbbbe9a22cbdf7e..73a763bc3bccd6601cae69b8a535b48fd6dba5cb 100644 --- a/applications/solvers/multiphase/icoReactingMultiphaseInterFoam/phasesSystem/MassTransferPhaseSystem/MassTransferPhaseSystem.C +++ b/applications/solvers/multiphase/icoReactingMultiphaseInterFoam/phasesSystem/MassTransferPhaseSystem/MassTransferPhaseSystem.C @@ -5,7 +5,7 @@ \\ / A nd | www.openfoam.com \\/ M anipulation | ------------------------------------------------------------------------------- - Copyright (C) 2017 OpenCFD Ltd. + Copyright (C) 2017-202 OpenCFD Ltd. ------------------------------------------------------------------------------- License This file is part of OpenFOAM. @@ -32,6 +32,7 @@ License #include "fvcDiv.H" #include "fvmSup.H" #include "fvMatrix.H" +#include "volFields.H" #include "fundamentalConstants.H" // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // @@ -152,9 +153,7 @@ Foam::MassTransferPhaseSystem::dmdt ( "dmdt", this->mesh().time().timeName(), - this->mesh(), - IOobject::NO_READ, - IOobject::NO_WRITE + this->mesh() ), this->mesh(), dimensionedScalar(dimDensity/dimTime, Zero) @@ -213,15 +212,45 @@ Foam::MassTransferPhaseSystem::heatTransfer ( "tdmdtYki", this->mesh().time().timeName(), - this->mesh(), - IOobject::NO_READ, - IOobject::NO_WRITE + this->mesh() ), this->mesh(), dimensionedScalar(dimDensity/dimTime, Zero) ) ); volScalarField& dmdtNetki = tdmdtNetki.ref(); + + tmp tSp + ( + new volScalarField + ( + IOobject + ( + "Sp", + this->mesh().time().timeName(), + this->mesh() + ), + this->mesh(), + dimensionedScalar(dimDensity/dimTime/dimTemperature, Zero) + ) + ); + volScalarField& Sp = tSp.ref(); + + tmp tSu + ( + new volScalarField + ( + IOobject + ( + "Su", + this->mesh().time().timeName(), + this->mesh() + ), + this->mesh(), + dimensionedScalar(dimDensity/dimTime, Zero) + ) + ); + volScalarField& Su = tSu.ref(); if (massTransferModels_.found(keyik)) @@ -229,14 +258,28 @@ Foam::MassTransferPhaseSystem::heatTransfer autoPtr& interfacePtr = massTransferModels_[keyik]; - // Explicit temperature mass transfer rate - tmp Kexp = - interfacePtr->Kexp(interfaceCompositionModel::T, T); - - if (Kexp.valid()) + dmdtNetki -= *dmdt_[keyik]; + + tmp KSp = + interfacePtr->KSp(interfaceCompositionModel::T, T); + + if (KSp.valid()) { - dmdtNetki -= Kexp.ref(); - *dmdt_[keyik] = Kexp.ref(); + Sp -= KSp.ref(); + } + + tmp KSu = + interfacePtr->KSu(interfaceCompositionModel::T, T); + + if (KSu.valid()) + { + Su -= KSu.ref(); + } + + // If linearization is not provided used full explicit + if (!KSp.valid() && !KSu.valid()) + { + Su -= *dmdt_[keyik]; } } @@ -246,37 +289,438 @@ Foam::MassTransferPhaseSystem::heatTransfer autoPtr& interfacePtr = massTransferModels_[keyki]; - // Explicit temperature mass transfer rate - const tmp Kexp = - interfacePtr->Kexp(interfaceCompositionModel::T, T); + dmdtNetki += *dmdt_[keyki]; - if (Kexp.valid()) + + tmp KSp = + interfacePtr->KSp(interfaceCompositionModel::T, T); + + if (KSp.valid()) { - dmdtNetki += Kexp.ref(); - *dmdt_[keyki] = Kexp.ref(); + Sp += KSp.ref(); + } + + tmp KSu = + interfacePtr->KSu(interfaceCompositionModel::T, T); + + if (KSu.valid()) + { + Su += KSu.ref(); + } + + // If linearization is not provided used full explicit + if (!KSp.valid() && !KSu.valid()) + { + Su += *dmdt_[keyki]; } - } - word keyikName(phasei.name() + phasek.name()); - word keykiName(phasek.name() + phasei.name()); + tmp L = calculateL(dmdtNetki, keyik, keyki, T); - eqn -= - ( - dmdtNetki - *( - calculateL(dmdtNetki, keyik, keyki, T) - - (phasek.Cp() - phasei.Cp()) - * constant::standard::Tstd - ) + //eqn -= dmdtNetki*L; + eqn -= fvm::Sp(Sp*L.ref(), T) + Su*L.ref(); + } + } + } + return tEqnPtr; +} + + +template +Foam::tmp +Foam::MassTransferPhaseSystem::volTransfer +( + const volScalarField& p +) +{ + tmp tEqnPtr + ( + new fvScalarMatrix(p, dimVolume/dimTime) + ); + + fvScalarMatrix& eqn = tEqnPtr.ref(); + + tmp tSp + ( + new volScalarField + ( + IOobject + ( + "Sp", + this->mesh().time().timeName(), + this->mesh() + ), + this->mesh(), + dimensionedScalar(dimless/dimTime/dimPressure, Zero) + ) + ); + volScalarField& Sp = tSp.ref(); + + tmp tSu + ( + new volScalarField + ( + IOobject + ( + "Su", + this->mesh().time().timeName(), + this->mesh() + ), + this->mesh(), + dimensionedScalar(dimless/dimTime, Zero) + ) + ); + volScalarField& Su = tSu.ref(); + + forAllConstIters(this->totalPhasePairs(), iter) + { + const phasePair& pair = iter()(); + + const phaseModel& phase1 = pair.phase1(); + const phaseModel& phase2 = pair.phase2(); + + const phasePairKey key12 + ( + phase1.name(), + phase2.name(), + true + ); + + if (massTransferModels_.found(key12)) + { + autoPtr& interfacePtr = + massTransferModels_[key12]; + + tmp KSp = + interfacePtr->KSp(interfaceCompositionModel::P, p); + + if (KSp.valid()) + { + Sp += KSp.ref(); + } + + tmp KSu = + interfacePtr->KSu(interfaceCompositionModel::P, p); + + if (KSu.valid()) + { + Su += KSu.ref(); + } + + // If linearization is not provided used full explicit + if (!KSp.valid() && !KSu.valid()) + { + Su -= + *dmdt_[key12] + *( + - this->coeffs(phase1.name()) + + this->coeffs(phase2.name()) + ); + } + } + + const phasePairKey key21 + ( + phase2.name(), + phase1.name(), + true + ); + + if (massTransferModels_.found(key21)) + { + autoPtr& interfacePtr = + massTransferModels_[key21]; + + tmp KSp = + interfacePtr->KSp(interfaceCompositionModel::P, p); + + if (KSp.valid()) + { + Sp += KSp.ref(); + } + + tmp KSu = + interfacePtr->KSu(interfaceCompositionModel::P, p); + + if (KSu.valid()) + { + Su += KSu.ref(); + } + + // If linearization is not provided used full explicit + if (!KSp.valid() && !KSu.valid()) + { + Su += + *dmdt_[key21] + *( + - this->coeffs(phase1.name()) + + this->coeffs(phase2.name()) ); } } + } + + eqn += fvm::Sp(Sp, p) + Su; return tEqnPtr; } +template +void Foam::MassTransferPhaseSystem::correctMassSources +( + const volScalarField& T +) +{ + forAllConstIters(this->phaseModels_, iteri) + { + const phaseModel& phasei = iteri()(); + + auto iterk = iteri; + + for (++iterk; iterk != this->phaseModels_.end(); ++iterk) + { + if (iteri()().name() != iterk()().name()) + { + const phaseModel& phasek = iterk()(); + + // Phase i to phase k + const phasePairKey keyik(phasei.name(), phasek.name(), true); + + // Phase k to phase i + const phasePairKey keyki(phasek.name(), phasei.name(), true); + + if (massTransferModels_.found(keyik)) + { + autoPtr& interfacePtr = + massTransferModels_[keyik]; + + tmp Kexp = interfacePtr->Kexp(T); + + *dmdt_[keyik] = Kexp.ref(); + + } + + if (massTransferModels_.found(keyki)) + { + autoPtr& interfacePtr = + massTransferModels_[keyki]; + + // Explicit temperature mass transfer rate + const tmp Kexp = interfacePtr->Kexp(T); + + *dmdt_[keyki] = Kexp.ref(); + } + } + } + } +} + + +template +void Foam::MassTransferPhaseSystem::alphaTransfer +( + SuSpTable& Su, + SuSpTable& Sp +) +{ + // This term adds and substract alpha*div(U) as a source term + // for alpha, substituting div(U) = mDot(1/rho1 - 1/rho2) + bool includeDivU(true); + + forAllConstIters(this->totalPhasePairs(), iter) + { + const phasePair& pair = iter()(); + + const phaseModel& phase1 = pair.phase1(); + const phaseModel& phase2 = pair.phase2(); + + const volScalarField& alpha1 = pair.phase1(); + const volScalarField& alpha2 = pair.phase2(); + + tmp tCoeffs1 = this->coeffs(phase1.name()); + const volScalarField& coeffs1 = tCoeffs1(); + + tmp tCoeffs2 = this->coeffs(phase2.name()); + const volScalarField& coeffs2 = tCoeffs2(); + + // Phase 1 to phase 2 + const phasePairKey key12 + ( + phase1.name(), + phase2.name(), + true + ); + + tmp tdmdt12(this->dmdt(key12)); + volScalarField& dmdt12 = tdmdt12.ref(); + + if (massTransferModels_.found(key12)) + { + autoPtr& interfacePtr = + massTransferModels_[key12]; + + tmp KSu = + interfacePtr->KSu(interfaceCompositionModel::alpha, phase1); + + if (KSu.valid()) + { + dmdt12 = KSu.ref(); + } + + includeDivU = interfacePtr->includeDivU(); + } + + + // Phase 2 to phase 1 + const phasePairKey key21 + ( + phase2.name(), + phase1.name(), + true + ); + + tmp tdmdt21(this->dmdt(key21)); + volScalarField& dmdt21 = tdmdt21.ref(); + + if (massTransferModels_.found(key21)) + { + autoPtr& interfacePtr = + massTransferModels_[key21]; + + tmp KSu = + interfacePtr->KSu(interfaceCompositionModel::alpha, phase2); + + if (KSu.valid()) + { + dmdt21 = KSu.ref(); + } + + includeDivU = interfacePtr->includeDivU(); + } + + volScalarField::Internal& SpPhase1 = Sp[phase1.name()]; + + volScalarField::Internal& SuPhase1 = Su[phase1.name()]; + + volScalarField::Internal& SpPhase2 = Sp[phase2.name()]; + + volScalarField::Internal& SuPhase2 = Su[phase2.name()]; + + const volScalarField dmdtNet(dmdt21 - dmdt12); + + const volScalarField coeffs12(coeffs1 - coeffs2); + + const surfaceScalarField& phi = this->phi(); + + if (includeDivU) + { + SuPhase1 += + fvc::div(phi)*min(max(alpha1, scalar(0)), scalar(1)); + + SuPhase2 += + fvc::div(phi)*min(max(alpha2, scalar(0)), scalar(1)); + } + + // NOTE: dmdtNet is distributed in terms = + // Source for phase 1 = + // dmdtNet/rho1 + // - alpha1*dmdtNet(1/rho1 - 1/rho2) + + forAll(dmdtNet, celli) + { + scalar dmdt21 = dmdtNet[celli]; + scalar coeffs12Cell = coeffs12[celli]; + + scalar alpha1Limited = max(min(alpha1[celli], 1.0), 0.0); + + // exp. + SuPhase1[celli] += coeffs1[celli]*dmdt21; + + if (includeDivU) + { + if (dmdt21 > 0) + { + if (coeffs12Cell > 0) + { + // imp + SpPhase1[celli] -= dmdt21*coeffs12Cell; + } + else if (coeffs12Cell < 0) + { + // exp + SuPhase1[celli] -= + dmdt21*coeffs12Cell*alpha1Limited; + } + } + else if (dmdt21 < 0) + { + if (coeffs12Cell > 0) + { + // exp + SuPhase1[celli] -= + dmdt21*coeffs12Cell*alpha1Limited; + } + else if (coeffs12Cell < 0) + { + // imp + SpPhase1[celli] -= dmdt21*coeffs12Cell; + } + } + } + } + + forAll(dmdtNet, celli) + { + scalar dmdt12 = -dmdtNet[celli]; + scalar coeffs21Cell = -coeffs12[celli]; + + scalar alpha2Limited = max(min(alpha2[celli], 1.0), 0.0); + + // exp + SuPhase2[celli] += coeffs2[celli]*dmdt12; + + if (includeDivU) + { + if (dmdt12 > 0) + { + if (coeffs21Cell > 0) + { + // imp + SpPhase2[celli] -= dmdt12*coeffs21Cell; + } + else if (coeffs21Cell < 0) + { + // exp + SuPhase2[celli] -= + dmdt12*coeffs21Cell*alpha2Limited; + } + } + else if (dmdt12 < 0) + { + if (coeffs21Cell > 0) + { + // exp + SuPhase2[celli] -= + coeffs21Cell*dmdt12*alpha2Limited; + } + else if (coeffs21Cell < 0) + { + // imp + SpPhase2[celli] -= dmdt12*coeffs21Cell; + } + } + } + + } + + // Update ddtAlphaMax + this->ddtAlphaMax_ = + max(gMax((dmdt21*coeffs1)()), gMax((dmdt12*coeffs2)())); + } +} + + template void Foam::MassTransferPhaseSystem::massSpeciesTransfer ( diff --git a/applications/solvers/multiphase/icoReactingMultiphaseInterFoam/phasesSystem/MassTransferPhaseSystem/MassTransferPhaseSystem.H b/applications/solvers/multiphase/icoReactingMultiphaseInterFoam/phasesSystem/MassTransferPhaseSystem/MassTransferPhaseSystem.H index fadfe2344ac1bd290d08f1d4796afe48a67dd417..0a6e77ab10db009b616c911eb95600ef136a1be5 100644 --- a/applications/solvers/multiphase/icoReactingMultiphaseInterFoam/phasesSystem/MassTransferPhaseSystem/MassTransferPhaseSystem.H +++ b/applications/solvers/multiphase/icoReactingMultiphaseInterFoam/phasesSystem/MassTransferPhaseSystem/MassTransferPhaseSystem.H @@ -5,7 +5,7 @@ \\ / A nd | www.openfoam.com \\/ M anipulation | ------------------------------------------------------------------------------- - Copyright (C) 2017 OpenCFD Ltd. + Copyright (C) 2017-2020 OpenCFD Ltd. ------------------------------------------------------------------------------- License This file is part of OpenFOAM. @@ -59,7 +59,7 @@ public: // Public typedefs - typedef + typedef HashTable < autoPtr, @@ -67,6 +67,9 @@ public: phasePairKey::hash > massTransferModelTable; + + + typedef HashTable SuSpTable; protected: @@ -123,9 +126,17 @@ public: // Mass transfer functions - //- Return the heat transfer matrix and fill dmdt for phases + //- Return the heat transfer matrix virtual tmp heatTransfer(const volScalarField& T); - + + //- Return the volumetric rate transfer matrix + virtual tmp volTransfer(const volScalarField& p); + + //- Correct/calculates mass sources dmdt for phases + virtual void correctMassSources(const volScalarField& T); + + //- Calculate mass transfer for alpha's + virtual void alphaTransfer(SuSpTable& Su, SuSpTable& Sp); //- Calculate mass transfer for species virtual void massSpeciesTransfer diff --git a/applications/solvers/multiphase/icoReactingMultiphaseInterFoam/phasesSystem/phaseModel/MultiComponentPhaseModel/MultiComponentPhaseModel.C b/applications/solvers/multiphase/icoReactingMultiphaseInterFoam/phasesSystem/phaseModel/MultiComponentPhaseModel/MultiComponentPhaseModel.C index 7bf399d3e5fc577a43f678d94c00227000c58004..d06b50380c1cbb7f48d6858f2f597054e3eb1f4d 100644 --- a/applications/solvers/multiphase/icoReactingMultiphaseInterFoam/phasesSystem/phaseModel/MultiComponentPhaseModel/MultiComponentPhaseModel.C +++ b/applications/solvers/multiphase/icoReactingMultiphaseInterFoam/phasesSystem/phaseModel/MultiComponentPhaseModel/MultiComponentPhaseModel.C @@ -90,9 +90,7 @@ MultiComponentPhaseModel ( IOobject::groupName("X" + species_[i], phaseName), fluid.mesh().time().timeName(), - fluid.mesh(), - IOobject::NO_READ, - IOobject::NO_WRITE + fluid.mesh() ), Y()[i] ) diff --git a/applications/solvers/multiphase/icoReactingMultiphaseInterFoam/phasesSystem/phaseSystem/multiphaseSystem/multiphaseSystem.C b/applications/solvers/multiphase/icoReactingMultiphaseInterFoam/phasesSystem/phaseSystem/multiphaseSystem/multiphaseSystem.C index 5fbf895fdeecef7f18e254388489ed9fee0778c2..a6dcf73aa751a1027d43ad5edaae620c41fd8309 100644 --- a/applications/solvers/multiphase/icoReactingMultiphaseInterFoam/phasesSystem/phaseSystem/multiphaseSystem/multiphaseSystem.C +++ b/applications/solvers/multiphase/icoReactingMultiphaseInterFoam/phasesSystem/phaseSystem/multiphaseSystem/multiphaseSystem.C @@ -5,7 +5,7 @@ \\ / A nd | www.openfoam.com \\/ M anipulation | ------------------------------------------------------------------------------- - Copyright (C) 2017-2019 OpenCFD Ltd. + Copyright (C) 2017-2020 OpenCFD Ltd. ------------------------------------------------------------------------------- License This file is part of OpenFOAM. @@ -120,166 +120,60 @@ Foam::multiphaseSystem::multiphaseSystem // * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * * // void Foam::multiphaseSystem::calculateSuSp() -{ - forAllConstIters(totalPhasePairs_, iter) - { - const phasePair& pair = iter()(); - - const phaseModel& phase1 = pair.phase1(); - const phaseModel& phase2 = pair.phase2(); - - const volScalarField& alpha1 = pair.phase1(); - const volScalarField& alpha2 = pair.phase2(); - - tmp tCoeffs1 = this->coeffs(phase1.name()); - const volScalarField& coeffs1 = tCoeffs1(); +{ + this->alphaTransfer(Su_, Sp_); +} - tmp tCoeffs2 = this->coeffs(phase2.name()); - const volScalarField& coeffs2 = tCoeffs2(); - // Phase 1 to phase 2 - const phasePairKey key12 +void Foam::multiphaseSystem::solve() +{ + const dictionary& alphaControls = mesh_.solverDict("alpha"); + label nAlphaSubCycles(alphaControls.get