From 01b4bf55efe52826d483aedd405188687ffaa763 Mon Sep 17 00:00:00 2001 From: Will Bainbridge Date: Thu, 10 Nov 2022 13:49:13 +0000 Subject: [PATCH] multiphaseEuler: Improved mass transfer linearisation w.r.t pressure Mass transfer rates are now updated following a change in the pressure if the mass transfer modelling provides a pressure coefficient. This means that pimple correctors can be used to improve the behaviour of mass transfer processes that coupled closely to the pressure field. --- .../multiphaseEuler/cellPressureCorrector.C | 15 ++++++++++++++- .../multiphaseEuler/compressibilityEqns.C | 9 +++++---- .../multiphaseEuler/facePressureCorrector.C | 15 ++++++++++++++- .../multiphaseEuler/multiphaseEuler.H | 6 +++++- .../multiphaseEuler/pressureCorrector.C | 5 ++++- 5 files changed, 42 insertions(+), 8 deletions(-) diff --git a/applications/solvers/modules/fluid/multiphaseEuler/multiphaseEuler/cellPressureCorrector.C b/applications/solvers/modules/fluid/multiphaseEuler/multiphaseEuler/cellPressureCorrector.C index 13be34880a..e7a5866d80 100644 --- a/applications/solvers/modules/fluid/multiphaseEuler/multiphaseEuler/cellPressureCorrector.C +++ b/applications/solvers/modules/fluid/multiphaseEuler/multiphaseEuler/cellPressureCorrector.C @@ -102,6 +102,10 @@ void Foam::solvers::multiphaseEuler::cellPressureCorrector() // Explicit force fluxes PtrList phiFs(fluid.phiFs(rAUs)); + // Mass transfer rates + PtrList dmdts(fluid.dmdts()); + PtrList d2mdtdps(fluid.d2mdtdps()); + // --- Pressure corrector loop while (pimple.correct()) { @@ -282,7 +286,7 @@ void Foam::solvers::multiphaseEuler::cellPressureCorrector() } // Compressible pressure equations - PtrList pEqnComps(compressibilityEqns()); + PtrList pEqnComps(compressibilityEqns(dmdts, d2mdtdps)); // Cache p prior to solve for density update volScalarField p_rgh_0(p_rgh); @@ -407,6 +411,15 @@ void Foam::solvers::multiphaseEuler::cellPressureCorrector() phase.thermoRef().rho() += phase.thermo().psi()*(p_rgh - p_rgh_0); } + // Update mass transfer rates for change in p_rgh + forAll(phases, phasei) + { + if (dmdts.set(phasei) && d2mdtdps.set(phasei)) + { + dmdts[phasei] += d2mdtdps[phasei]*(p_rgh - p_rgh_0); + } + } + // Correct p_rgh for consistency with p and the updated densities rho = fluid.rho(); p_rgh = p - rho*buoyancy.gh; diff --git a/applications/solvers/modules/fluid/multiphaseEuler/multiphaseEuler/compressibilityEqns.C b/applications/solvers/modules/fluid/multiphaseEuler/multiphaseEuler/compressibilityEqns.C index 5ef3e1ad0f..68e730ed63 100644 --- a/applications/solvers/modules/fluid/multiphaseEuler/multiphaseEuler/compressibilityEqns.C +++ b/applications/solvers/modules/fluid/multiphaseEuler/multiphaseEuler/compressibilityEqns.C @@ -34,13 +34,14 @@ License // * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * * // Foam::PtrList -Foam::solvers::multiphaseEuler::compressibilityEqns() const +Foam::solvers::multiphaseEuler::compressibilityEqns +( + const PtrList& dmdts, + const PtrList& d2mdtdps +) const { PtrList pEqnComps(phases.size()); - PtrList dmdts(fluid.dmdts()); - PtrList d2mdtdps(fluid.d2mdtdps()); - forAll(phases, phasei) { phaseModel& phase = phases[phasei]; diff --git a/applications/solvers/modules/fluid/multiphaseEuler/multiphaseEuler/facePressureCorrector.C b/applications/solvers/modules/fluid/multiphaseEuler/multiphaseEuler/facePressureCorrector.C index 61c0a9acfe..3671416d51 100644 --- a/applications/solvers/modules/fluid/multiphaseEuler/multiphaseEuler/facePressureCorrector.C +++ b/applications/solvers/modules/fluid/multiphaseEuler/multiphaseEuler/facePressureCorrector.C @@ -113,6 +113,10 @@ void Foam::solvers::multiphaseEuler::facePressureCorrector() // Explicit force fluxes PtrList phiFfs(fluid.phiFfs(rAUfs)); + // Mass transfer rates + PtrList dmdts(fluid.dmdts()); + PtrList d2mdtdps(fluid.d2mdtdps()); + // --- Pressure corrector loop while (pimple.correct()) { @@ -269,7 +273,7 @@ void Foam::solvers::multiphaseEuler::facePressureCorrector() } // Compressible pressure equations - PtrList pEqnComps(compressibilityEqns()); + PtrList pEqnComps(compressibilityEqns(dmdts, d2mdtdps)); // Cache p prior to solve for density update volScalarField p_rgh_0(p_rgh); @@ -381,6 +385,15 @@ void Foam::solvers::multiphaseEuler::facePressureCorrector() phase.thermoRef().rho() += phase.thermo().psi()*(p_rgh - p_rgh_0); } + // Update mass transfer rates for change in p_rgh + forAll(phases, phasei) + { + if (dmdts.set(phasei) && d2mdtdps.set(phasei)) + { + dmdts[phasei] += d2mdtdps[phasei]*(p_rgh - p_rgh_0); + } + } + // Correct p_rgh for consistency with p and the updated densities rho = fluid.rho(); p_rgh = p - rho*buoyancy.gh; diff --git a/applications/solvers/modules/fluid/multiphaseEuler/multiphaseEuler/multiphaseEuler.H b/applications/solvers/modules/fluid/multiphaseEuler/multiphaseEuler/multiphaseEuler.H index 117be21c6a..c4504ba887 100644 --- a/applications/solvers/modules/fluid/multiphaseEuler/multiphaseEuler/multiphaseEuler.H +++ b/applications/solvers/modules/fluid/multiphaseEuler/multiphaseEuler/multiphaseEuler.H @@ -185,7 +185,11 @@ private: void facePressureCorrector(); //- Return the list of pressure equation compressibility contributions - PtrList compressibilityEqns() const; + PtrList compressibilityEqns + ( + const PtrList& dmdts, + const PtrList& d2mdtdps + ) const; public: diff --git a/applications/solvers/modules/fluid/multiphaseEuler/multiphaseEuler/pressureCorrector.C b/applications/solvers/modules/fluid/multiphaseEuler/multiphaseEuler/pressureCorrector.C index 7079c46131..71ac2fe2f1 100644 --- a/applications/solvers/modules/fluid/multiphaseEuler/multiphaseEuler/pressureCorrector.C +++ b/applications/solvers/modules/fluid/multiphaseEuler/multiphaseEuler/pressureCorrector.C @@ -44,7 +44,10 @@ void Foam::solvers::multiphaseEuler::pressureCorrector() } else { - PtrList pEqnComps(compressibilityEqns()); + PtrList pEqnComps + ( + compressibilityEqns(fluid.dmdts(), fluid.d2mdtdps()) + ); forAll(phases, phasei) {