From 000876bbc64c3cf347cd14e142cc40c727ac4491 Mon Sep 17 00:00:00 2001 From: Will Bainbridge Date: Tue, 13 Apr 2021 15:47:38 +0100 Subject: [PATCH] multiphaseEulerFoam: Pressure implicit thermal mass transfer The phase system now have the ability to specify the derivative of mass transfer rates w.r.t. pressure. This permits implicit handling of pressure-coupled mass transfer processes. This implicit handling has been applied to the mass transfers that are modelled by the thermal phase system. This should result in significant stability improvements. The implicit handling can be toggled on or off by means of a "pressureImplicit" switch in constant/phaseProperties. It is on by default. Patch contributed by Juho Peltola, VTT. --- .../multiphaseEulerFoam/pEqnComps.H | 5 ++ .../ThermalPhaseChangePhaseSystem.C | 71 ++++++++++++++++++- .../ThermalPhaseChangePhaseSystem.H | 11 ++- .../phaseSystems/phaseSystem/phaseSystem.C | 6 ++ .../phaseSystems/phaseSystem/phaseSystem.H | 4 ++ 5 files changed, 95 insertions(+), 2 deletions(-) diff --git a/applications/solvers/multiphase/multiphaseEulerFoam/multiphaseEulerFoam/pEqnComps.H b/applications/solvers/multiphase/multiphaseEulerFoam/multiphaseEulerFoam/pEqnComps.H index 84c283e7a7..e7601ee5fb 100644 --- a/applications/solvers/multiphase/multiphaseEulerFoam/multiphaseEulerFoam/pEqnComps.H +++ b/applications/solvers/multiphase/multiphaseEulerFoam/multiphaseEulerFoam/pEqnComps.H @@ -2,6 +2,7 @@ PtrList pEqnComps(phases.size()); { PtrList dmdts(fluid.dmdts()); + PtrList d2mdtdps(fluid.d2mdtdps()); forAll(phases, phasei) { @@ -71,5 +72,9 @@ PtrList pEqnComps(phases.size()); { pEqnComp -= dmdts[phasei]/rho; } + if (d2mdtdps.set(phasei)) + { + pEqnComp -= correction(fvm::Sp(d2mdtdps[phasei]/rho, p_rgh)); + } } } diff --git a/applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/PhaseSystems/ThermalPhaseChangePhaseSystem/ThermalPhaseChangePhaseSystem.C b/applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/PhaseSystems/ThermalPhaseChangePhaseSystem/ThermalPhaseChangePhaseSystem.C index 272386b579..23a851c79b 100644 --- a/applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/PhaseSystems/ThermalPhaseChangePhaseSystem/ThermalPhaseChangePhaseSystem.C +++ b/applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/PhaseSystems/ThermalPhaseChangePhaseSystem/ThermalPhaseChangePhaseSystem.C @@ -68,7 +68,11 @@ ThermalPhaseChangePhaseSystem : BasePhaseSystem(mesh), volatile_(this->template lookupOrDefault("volatile", "none")), - dmdt0s_(this->phases().size()) + dmdt0s_(this->phases().size()), + pressureImplicit_ + ( + this->template lookupOrDefault("pressureImplicit", true) + ) { this->generatePairsAndSubModels ( @@ -137,6 +141,28 @@ ThermalPhaseChangePhaseSystem ) ); + d2mdtdpfs_.insert + ( + pair, + new volScalarField + ( + IOobject + ( + IOobject::groupName + ( + "thermalPhaseChange:d2mdtdpf", + pair.name() + ), + this->mesh().time().timeName(), + this->mesh(), + IOobject::NO_READ, + IOobject::AUTO_WRITE + ), + this->mesh(), + dimensionedScalar((dimDensity/dimTime)/dimPressure, 0) + ) + ); + Tfs_.insert ( pair, @@ -242,6 +268,25 @@ Foam::ThermalPhaseChangePhaseSystem::dmdts() const } +template +Foam::PtrList +Foam::ThermalPhaseChangePhaseSystem::d2mdtdps() const +{ + PtrList d2mdtdps(BasePhaseSystem::d2mdtdps()); + + forAllConstIter(phaseSystem::dmdtfTable, d2mdtdpfs_, d2mdtdpfIter) + { + const phasePair& pair = this->phasePairs_[d2mdtdpfIter.key()]; + const volScalarField& d2mdtdpf = *d2mdtdpfIter(); + + addField(pair.phase1(), "d2mdtdp", d2mdtdpf, d2mdtdps); + addField(pair.phase2(), "d2mdtdp", - d2mdtdpf, d2mdtdps); + } + + return d2mdtdps; +} + + template Foam::autoPtr Foam::ThermalPhaseChangePhaseSystem:: @@ -498,6 +543,30 @@ Foam::ThermalPhaseChangePhaseSystem::correctInterfaceThermo() + pos(dmdtfNew)*phase2.Y(volatile_); } + if (pressureImplicit_) + { + volScalarField& d2mdtdpf(*this->d2mdtdpfs_[pair]); + + const dimensionedScalar dp(rootSmall*thermo1.p().average()); + + const volScalarField dTsatdp + ( + ( + saturationModelIter()->Tsat(thermo1.p() + dp/2) + - saturationModelIter()->Tsat(thermo1.p() - dp/2) + )/dp + ); + + d2mdtdpf = (H1 + H2)*dTsatdp/L; + + if (volatile_ != "none") + { + d2mdtdpf *= + neg0(dmdtfNew)*phase1.Y(volatile_) + + pos(dmdtfNew)*phase2.Y(volatile_); + } + } + H1 = this->heatTransferModels_[pair].first()->K(); H2 = this->heatTransferModels_[pair].second()->K(); diff --git a/applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/PhaseSystems/ThermalPhaseChangePhaseSystem/ThermalPhaseChangePhaseSystem.H b/applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/PhaseSystems/ThermalPhaseChangePhaseSystem/ThermalPhaseChangePhaseSystem.H index 94e805fd5b..a594b8c963 100644 --- a/applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/PhaseSystems/ThermalPhaseChangePhaseSystem/ThermalPhaseChangePhaseSystem.H +++ b/applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/PhaseSystems/ThermalPhaseChangePhaseSystem/ThermalPhaseChangePhaseSystem.H @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | Website: https://openfoam.org - \\ / A nd | Copyright (C) 2015-2020 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2015-2021 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -101,6 +101,9 @@ class ThermalPhaseChangePhaseSystem //- Mass transfer rates phaseSystem::dmdtfTable dmdtfs_; + //- Mass transfer linearisation coeffs + phaseSystem::dmdtfTable d2mdtdpfs_; + //- Interface temperatures phaseSystem::dmdtfTable Tfs_; @@ -114,6 +117,9 @@ class ThermalPhaseChangePhaseSystem // function PtrList dmdt0s_; + //- Switch to control whether or not mass transfer rates are linearised + // in the pressure equation + Switch pressureImplicit_; // Private Member Functions @@ -144,6 +150,9 @@ public: //- Return the mass transfer rates for each phase virtual PtrList dmdts() const; + //- Return the mass transfer linearisation coeffs for each phase + virtual PtrList d2mdtdps() const; + //- Return the momentum transfer matrices for the cell-based algorithm virtual autoPtr momentumTransfer(); diff --git a/applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/phaseSystem/phaseSystem.C b/applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/phaseSystem/phaseSystem.C index 7bc08e3f57..8d05a2f768 100644 --- a/applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/phaseSystem/phaseSystem.C +++ b/applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/phaseSystem/phaseSystem.C @@ -671,6 +671,12 @@ Foam::PtrList Foam::phaseSystem::dmdts() const } +Foam::PtrList Foam::phaseSystem::d2mdtdps() const +{ + return PtrList(phaseModels_.size()); +} + + bool Foam::phaseSystem::incompressible() const { forAll(phaseModels_, phasei) diff --git a/applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/phaseSystem/phaseSystem.H b/applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/phaseSystem/phaseSystem.H index 80e662083b..8264d34b8c 100644 --- a/applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/phaseSystem/phaseSystem.H +++ b/applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/phaseSystem/phaseSystem.H @@ -544,6 +544,10 @@ public: //- Return the mass transfer rates for each phase virtual PtrList dmdts() const; + //- Return the mass transfer pressure implicit coefficients + // for each phase + virtual PtrList d2mdtdps() const; + //- Return incompressibility bool incompressible() const;