populationBalanceModel: phaseChange: Improvements

A number of improvements have been made to the population balance phase
change drift model.

- The model now checks the ordering of the phase pairs and changes the
sign of the drift rate accordingly.

- The phase change mass flux and weights are calculated for each
velocity group, so the drift rate and phase change mass flux should be
consistent for each velocity group.

- By default the phase change mass flux is distributed between the size
groups based on the interfacial area of each group. For backward
compatibility number weighting can be enabled with a new
"numberWeighted" option.

The model now requires the user to provide a list of phase pairs in the
usual parenthesised form, rather than using the name. For example:

    phaseChange
    {
        pairs ((gas and liquid));
    }

Patch contributed by Juho Peltola, VTT.

applications/solvers/multiphase/reactingEulerFoam/phaseSystems/populationBalanceModel/driftModels/phaseChange/phaseChange.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2018 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2018-2019 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -52,57 +52,71 @@ Foam::diameterModels::driftModels::phaseChange::phaseChange
 )
 :
     driftModel(popBal, dict),
-    pairNames_(dict.lookup("pairNames")),
-    iDmdt_
-    (
-        IOobject
+    pairKeys_(dict.lookup("pairs")),
+    numberWeighted_(dict.lookupOrDefault<Switch>("numberWeighted", false)),
+    W_(pairKeys_.size())
+{
+    const phaseSystem& fluid = popBal_.fluid();
+
+    forAll(pairKeys_, i)
+    {
+        const phasePair& pair = fluid.phasePairs()[pairKeys_[i]];
+
+        W_.set
         (
-            "iDmdt",
-            popBal.time().timeName(),
-            popBal.mesh()
-        ),
-        popBal.mesh(),
-        dimensionedScalar(dimDensity/dimTime, Zero)
-    ),
-    N_
-    (
-        IOobject
-        (
-            "N",
-            popBal.mesh().time().timeName(),
-            popBal.mesh()
-        ),
-        popBal.mesh(),
-        dimensionedScalar(inv(dimVolume), Zero)
-    )
-{}
+            i,
+            new volScalarField
+            (
+                IOobject
+                (
+                    IOobject::groupName(type() + ":W", pair.name()),
+                    popBal_.mesh().time().timeName(),
+                    popBal_.mesh()
+                ),
+                popBal_.mesh(),
+                dimensionedScalar
+                (
+                    inv(numberWeighted_ ? dimVolume : dimLength),
+                    Zero
+                )
+            )
+        );
+    }
+}
 
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
-
 void Foam::diameterModels::driftModels::phaseChange::correct()
 {
-    iDmdt_ = Zero;
+    const phaseSystem& fluid = popBal_.fluid();
 
-    forAll(pairNames_, i)
+    forAll(pairKeys_, i)
     {
-        const word& pairName = pairNames_[i];
-
-        iDmdt_ +=
-            popBal_.mesh().lookupObject<volScalarField>
-            (
-                IOobject::groupName("iDmdt", pairName)
-            );
+        W_[i] = Zero;
     }
 
-    N_ = Zero;
-
-    forAll(popBal_.sizeGroups(), i)
+    forAll(pairKeys_, k)
     {
-        const sizeGroup& fi = popBal_.sizeGroups()[i];
+        if (fluid.phasePairs().found(pairKeys_[k]))
+        {
+            const phasePair& pair = fluid.phasePairs()[pairKeys_[k]];
 
-        N_ += fi*max(fi.phase(), small)/fi.x();
+            forAll(popBal_.velocityGroups(), j)
+            {
+                const velocityGroup& vgj = popBal_.velocityGroups()[j];
+                if (pair.contains(vgj.phase()))
+                {
+                    forAll(vgj.sizeGroups(), i)
+                    {
+                        const sizeGroup& fi = vgj.sizeGroups()[i];
+                        W_[k] +=
+                            fi*max(fi.phase(), small)
+                           /(numberWeighted_ ? fi.x() : fi.d());
+                    }
+                }
+            }
+        }
     }
 }
 
@@ -113,9 +127,39 @@ void Foam::diameterModels::driftModels::phaseChange::addToDriftRate
     const label i
 )
 {
-    const sizeGroup& fi = popBal_.sizeGroups()[i];
+    const velocityGroup& vg = popBal_.sizeGroups()[i].VelocityGroup();
 
-    driftRate += iDmdt_/(N_*fi.phase().rho());
+    forAll(pairKeys_, k)
+    {
+        const phasePair& pair =
+                popBal_.fluid().phasePairs()[pairKeys_[k]];
+
+        if (pair.contains(vg.phase()))
+        {
+            const volScalarField& iDmdt =
+                popBal_.mesh().lookupObject<volScalarField>
+                (
+                    IOobject::groupName("iDmdt", pair.name())
+                );
+
+            const scalar iDmdtSign =
+                vg.phase().name() == pair.first() ? +1 : -1;
+
+            const sizeGroup& fi = popBal_.sizeGroups()[i];
+
+            tmp<volScalarField> dDriftRate
+            (
+                iDmdtSign*iDmdt/(fi.phase().rho()*W_[k])
+            );
+
+            if (!numberWeighted_)
+            {
+                dDriftRate.ref() *= fi.x()/fi.d();
+            }
+
+            driftRate += dDriftRate;
+        }
+    }
 }
 
 

applications/solvers/multiphase/reactingEulerFoam/phaseSystems/populationBalanceModel/driftModels/phaseChange/phaseChange.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2018 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2018-2019 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -25,10 +25,12 @@ Class
     Foam::diameterModels::driftModels::phaseChange
 
 Description
-    Drift induced by interfacial phaseChange.
+    Drift induced by interfacial phase change. By default phase change mass
+    flux is distributed between sizeGroups of each velocityGroup with phase
+    change based on interfacial area of each size group.
 
 SourceFiles
-    isothermal.C
+    phaseChange.C
 
 \*---------------------------------------------------------------------------*/
 
@@ -56,15 +58,17 @@ class phaseChange
 {
     // Private data
 
-        //- Names of unordered phasePairs between which phaseChange occurs, i.e.
-        //  "(gasIAndLiquid gasIIAndLiquid)"
-        List<word> pairNames_;
+        //- PhasePairs between which phaseChange occurs, e.g.,
+        //  "((gasI and liquid) (gasII and liquid))"
+        List<phasePairKey> pairKeys_;
 
-        //- Total mass transfer rate due to phaseChange
-        volScalarField iDmdt_;
+        //- Distribute phase change mass flux between sizeGroups based on the
+        //  number concentration, rather than the interfacial area
+        Switch numberWeighted_;
+
+        //- Weighting with which the phase change mass flux is distributed
+        PtrList<volScalarField> W_;
 
-        //- Total number concentration
-        volScalarField N_;
 
 public:
 
@@ -73,6 +77,7 @@ public:
 
     // Constructor
 
+        //- Construct from a population balance model and a dictionary
         phaseChange
         (
             const populationBalanceModel& popBal,

tutorials/multiphase/reactingTwoPhaseEulerFoam/RAS/wallBoilingPolyDisperse/constant/phaseProperties

@@ -99,7 +99,7 @@ populationBalanceCoeffs
         (
             phaseChange
             {
-                pairNames (gasAndLiquid);
+                pairs    ((gas and liquid));
             }
 
             densityChange{}

tutorials/multiphase/reactingTwoPhaseEulerFoam/RAS/wallBoilingPolyDisperse/system/fvSolution

@@ -28,7 +28,7 @@ solvers
         nCorr                   1;
         tolerance               1e-4;
         renormalizeOnRestart    true;
-        renormalize             false;
+        renormalize             true;
         solveOnFinalIterOnly    true;
     }
 
@@ -108,12 +108,13 @@ relaxationFactors
 {
     fields
     {
-        iDmdt           0.1;
+        iDmdt           0.2;
     }
 
     equations
     {
         ".*"            1;
+        "e\..*"         0.2;
     }
 }