multiphaseEulerFoam: Added interface composition population balance system

A phase system has been added that combines interface composition phase
change modelling and population balance modelling. It can be selected in
constant/phaseProperties using the following setting:

    type    interfaceCompositionPhaseChangePopulationBalanceMultiphaseSystem;

In order to facilitate the operation of population balance models with
mass transfers generated by the interface composition system, the full
interfacial mass transfer rate is now being stored, as well as the
linearised specie transfers. This means the mass transfer rates
associated with the interface composition mechanism can be looked up and
utilised by sub models.

Based on a patch contributed by Institute of Fluid Dynamics,
Helmholtz-Zentrum Dresden - Rossendorf (HZDR)

applications/solvers/multiphase/multiphaseEulerFoam/multiphaseEulerFoam/multiphaseSystems/multiphaseSystems.C

@@ -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
@@ -119,6 +119,30 @@ namespace Foam
         populationBalanceMultiphaseSystem
     );
 
+    typedef
+        InterfaceCompositionPhaseChangePhaseSystem
+        <
+            PopulationBalancePhaseSystem
+            <
+                PhaseTransferPhaseSystem
+                <
+                    TwoResistanceHeatTransferPhaseSystem
+                    <
+                        MomentumTransferPhaseSystem<phaseSystem>
+                    >
+                >
+            >
+        >
+        interfaceCompositionPhaseChangePopulationBalanceMultiphaseSystem;
+
+    addNamedToRunTimeSelectionTable
+    (
+        phaseSystem,
+        interfaceCompositionPhaseChangePopulationBalanceMultiphaseSystem,
+        dictionary,
+        interfaceCompositionPhaseChangePopulationBalanceMultiphaseSystem
+    );
+
     typedef
         ThermalPhaseChangePhaseSystem
         <

applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/PhaseSystems/InterfaceCompositionPhaseChangePhaseSystem/InterfaceCompositionPhaseChangePhaseSystem.C

@@ -34,10 +34,13 @@ License
 template<class BasePhaseSystem>
 Foam::autoPtr<Foam::phaseSystem::dmdtfTable>
 Foam::InterfaceCompositionPhaseChangePhaseSystem<BasePhaseSystem>::
-dmdtfs() const
+totalDmdtfs() const
 {
-    autoPtr<phaseSystem::dmdtfTable> dmdtPtr(new phaseSystem::dmdtfTable);
-    phaseSystem::dmdtfTable& dmdtfs = dmdtPtr();
+    autoPtr<phaseSystem::dmdtfTable> totalDmdtfsPtr
+    (
+        new phaseSystem::dmdtfTable
+    );
+    phaseSystem::dmdtfTable& totalDmdtfs = totalDmdtfsPtr();
 
     forAllConstIter
     (
@@ -79,29 +82,32 @@ dmdtfs() const
                     )
                 );
 
-                if (dmdtfs.found(pair))
+                if (totalDmdtfs.found(pair))
                 {
-                    *dmdtfs[pair] += dmidtf;
+                    *totalDmdtfs[pair] += dmidtf;
                 }
                 else
                 {
-                    dmdtfs.insert(pair, dmidtf.ptr());
+                    totalDmdtfs.insert(pair, dmidtf.ptr());
                 }
             }
         }
     }
 
-    return dmdtPtr;
+    return totalDmdtfsPtr;
 }
 
 
 template<class BasePhaseSystem>
 Foam::autoPtr<Foam::phaseSystem::dmidtfTable>
 Foam::InterfaceCompositionPhaseChangePhaseSystem<BasePhaseSystem>::
-dmidtfs() const
+totalDmidtfs() const
 {
-    autoPtr<phaseSystem::dmidtfTable> dmidtfsPtr(new phaseSystem::dmidtfTable);
-    phaseSystem::dmidtfTable& dmidtfs = dmidtfsPtr();
+    autoPtr<phaseSystem::dmidtfTable> totalDmidtfsPtr
+    (
+        new phaseSystem::dmidtfTable
+    );
+    phaseSystem::dmidtfTable& totalDmidtfs = totalDmidtfsPtr();
 
     forAllConstIter
     (
@@ -113,9 +119,9 @@ dmidtfs() const
         const phasePair& pair =
             this->phasePairs_[interfaceCompositionModelIter.key()];
 
-        if (!dmidtfs.found(pair))
+        if (!totalDmidtfs.found(pair))
         {
-            dmidtfs.insert(pair, new HashPtrTable<volScalarField>());
+            totalDmidtfs.insert(pair, new HashPtrTable<volScalarField>());
         }
 
         forAllConstIter(phasePair, pair, pairIter)
@@ -149,19 +155,19 @@ dmidtfs() const
                     )
                 );
 
-                if (dmidtfs[pair]->found(specie))
+                if (totalDmidtfs[pair]->found(specie))
                 {
-                    *(*dmidtfs[pair])[specie] += dmidtf;
+                    *(*totalDmidtfs[pair])[specie] += dmidtf;
                 }
                 else
                 {
-                    dmidtfs[pair]->insert(specie, dmidtf.ptr());
+                    totalDmidtfs[pair]->insert(specie, dmidtf.ptr());
                 }
             }
         }
     }
 
-    return dmidtfsPtr;
+    return totalDmidtfsPtr;
 }
 
 
@@ -258,6 +264,28 @@ InterfaceCompositionPhaseChangePhaseSystem
         const phasePair& pair =
             this->phasePairs_[interfaceCompositionModelIter.key()];
 
+        dmdtfs_.insert
+        (
+            pair,
+            new volScalarField
+            (
+                IOobject
+                (
+                    IOobject::groupName
+                    (
+                        "interfaceCompositionPhaseChange:dmdtf",
+                        pair.name()
+                    ),
+                    this->mesh().time().timeName(),
+                    this->mesh(),
+                    IOobject::READ_IF_PRESENT,
+                    IOobject::AUTO_WRITE
+                ),
+                this->mesh(),
+                dimensionedScalar(dimDensity/dimTime, Zero)
+            )
+        );
+
         dmidtfSus_.insert(pair, new HashPtrTable<volScalarField>());
 
         dmidtfSps_.insert(pair, new HashPtrTable<volScalarField>());
@@ -374,39 +402,11 @@ Foam::InterfaceCompositionPhaseChangePhaseSystem<BasePhaseSystem>::dmdtf
 {
     tmp<volScalarField> tDmdtf = BasePhaseSystem::dmdtf(key);
 
-    if (interfaceCompositionModels_.found(key))
+    if (dmdtfs_.found(key))
     {
-        const phasePair& pair = this->phasePairs_[key];
+        const label dmdtSign(Pair<word>::compare(this->phasePairs_[key], key));
 
-        forAllConstIter(phasePair, pair, pairIter)
-        {
-            const autoPtr<interfaceCompositionModel>& compositionModelPtr =
-                interfaceCompositionModels_[key][pairIter.index()];
-
-            if (!compositionModelPtr.valid()) continue;
-
-            const interfaceCompositionModel& compositionModel =
-                compositionModelPtr();
-
-            const phaseModel& phase = *pairIter;
-
-            forAllConstIter
-            (
-                hashedWordList,
-                compositionModel.species(),
-                specieIter
-            )
-            {
-                const word& specie = *specieIter;
-
-                tDmdtf.ref() +=
-                    Pair<word>::compare(pair, key)
-                   *(
-                        *(*dmidtfSus_[pair])[specie]
-                      + *(*dmidtfSps_[pair])[specie]*phase.Y(specie)
-                    );
-            }
-        }
+        tDmdtf.ref() += dmdtSign**dmdtfs_[key];
     }
 
     return tDmdtf;
@@ -419,10 +419,7 @@ Foam::InterfaceCompositionPhaseChangePhaseSystem<BasePhaseSystem>::dmdts() const
 {
     PtrList<volScalarField> dmdts(BasePhaseSystem::dmdts());
 
-    autoPtr<phaseSystem::dmdtfTable> dmdtfsPtr = this->dmdtfs();
-    const phaseSystem::dmdtfTable& dmdtfs = dmdtfsPtr();
-
-    forAllConstIter(phaseSystem::dmdtfTable, dmdtfs, dmdtfsIter)
+    forAllConstIter(phaseSystem::dmdtfTable, dmdtfs_, dmdtfsIter)
     {
         const phasePair& pair = this->phasePairs_[dmdtfsIter.key()];
         const phaseModel& phase = pair.phase1();
@@ -446,7 +443,7 @@ momentumTransfer()
 
     phaseSystem::momentumTransferTable& eqns = eqnsPtr();
 
-    this->addDmdtUfs(dmdtfs(), eqns);
+    this->addDmdtUfs(dmdtfs_, eqns);
 
     return eqnsPtr;
 }
@@ -462,7 +459,7 @@ momentumTransferf()
 
     phaseSystem::momentumTransferTable& eqns = eqnsPtr();
 
-    this->addDmdtUfs(dmdtfs(), eqns);
+    this->addDmdtUfs(dmdtfs_, eqns);
 
     return eqnsPtr;
 }
@@ -480,7 +477,7 @@ heatTransfer() const
 
     this->addDmidtHefs
     (
-        dmidtfs(),
+        totalDmidtfs(),
         Tfs_,
         latentHeatScheme::symmetric,
         latentHeatTransfer::mass,
@@ -503,7 +500,7 @@ specieTransfer() const
 
     // Explicit
     /*
-    this->addDmidtYf(dmidtfs(), eqns);
+    this->addDmidtYf(totalDmidtfs(), eqns);
     */
 
     // Semi-implicit
@@ -612,6 +609,18 @@ correct()
             }
         }
     }
+
+    dmdtfs_ = totalDmdtfs();
+}
+
+
+template<class BasePhaseSystem>
+void Foam::InterfaceCompositionPhaseChangePhaseSystem<BasePhaseSystem>::
+correctSpecies()
+{
+    BasePhaseSystem::correctSpecies();
+
+    dmdtfs_ = totalDmdtfs();
 }
 
 

applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/PhaseSystems/InterfaceCompositionPhaseChangePhaseSystem/InterfaceCompositionPhaseChangePhaseSystem.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
@@ -93,10 +93,13 @@ class InterfaceCompositionPhaseChangePhaseSystem
             //- Interface composition models
             interfaceCompositionModelTable interfaceCompositionModels_;
 
-            //- The explicit part of the interfacial mass transfer rates
+            //- Interface mass transfer rates
+            phaseSystem::dmdtfTable dmdtfs_;
+
+            //- The explicit part of the interface specie transfer rates
             phaseSystem::dmidtfTable dmidtfSus_;
 
-            //- The implicit part of the interfacial mass transfer rates
+            //- The implicit part of the interface specie transfer rates
             phaseSystem::dmidtfTable dmidtfSps_;
 
             //- Interface temperatures
@@ -106,10 +109,10 @@ class InterfaceCompositionPhaseChangePhaseSystem
     // Private Member Functions
 
         //- Return mass transfers across each interface
-        autoPtr<phaseSystem::dmdtfTable> dmdtfs() const;
+        autoPtr<phaseSystem::dmdtfTable> totalDmdtfs() const;
 
         //- Return species mass transfers across each interface
-        autoPtr<phaseSystem::dmidtfTable> dmidtfs() const;
+        autoPtr<phaseSystem::dmidtfTable> totalDmidtfs() const;
 
 
 public:
@@ -148,6 +151,9 @@ public:
         //- Correct the fluid properties other than those listed below
         virtual void correct();
 
+        //- Correct the mass transfer rates for the new species mass fractions
+        virtual void correctSpecies();
+
         //- Correct the interface temperatures
         virtual void correctInterfaceThermo();
 

applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/PhaseSystems/PopulationBalancePhaseSystem/PopulationBalancePhaseSystem.C

@@ -113,16 +113,16 @@ Foam::PopulationBalancePhaseSystem<BasePhaseSystem>::dmdtf
     const phasePairKey& key
 ) const
 {
-    tmp<volScalarField> tDmdt = BasePhaseSystem::dmdtf(key);
+    tmp<volScalarField> tDmdtf = BasePhaseSystem::dmdtf(key);
 
-    if (!dmdtfs_.found(key))
+    if (dmdtfs_.found(key))
     {
         const label dmdtSign(Pair<word>::compare(this->phasePairs_[key], key));
 
-        tDmdt.ref() += dmdtSign**dmdtfs_[key];
+        tDmdtf.ref() += dmdtSign**dmdtfs_[key];
     }
 
-    return tDmdt;
+    return tDmdtf;
 }