reacting*EulerFoam: Various consistency improvements

Mass transfer rates now have a more comprehensive naming convention.
"dmdt" means a bulk/mixture transfer, whilst "dmidt" is for a
specie-specific transfer. "dmdt" implies a transfer into a phase, whilst
"dmdtf" means a transfer across an interface. Tables or lists of
transfers are denoted by pluralising the name with the suffix "s"; e.g.,
"dmdtfs". All registered mass transfer rate fields have names which
include the name of the sub-model or phase system which generated them.

The phaseTransfer models have been changed so that the mixture and the
specie-specific mass transfers are independent. This simplifies the
naming convention required for registering the resulting mass transfers
and reduces the amount of logic necessary in the phase system.

The inheritance pattern of the alphat wall functions has been altered so
that the code and parameters relating to phase change are reused, and so
that the base (the Jayatilleke wall function) more closely resembles the
library implementation. This should make it easier to remove it when the
library function is generalised enough to use it directly.

The phaseSystem::zero*Field construction functions have been removed as
their behaviour regarding registration was not clear, and in most
instances of their usage the GeometriField<...>::New methods are
similarly convenient.

applications/solvers/multiphase/reactingEulerFoam/derivedFvPatchFields/Make/files

@@ -20,9 +20,9 @@ wallBoilingSubModels/departureFrequencyModels/departureFrequencyModel/newDepartu
 wallBoilingSubModels/departureFrequencyModels/Cole/Cole.C
 wallBoilingSubModels/departureFrequencyModels/KocamustafaogullariIshiiDepartureFrequency/KocamustafaogullariIshiiDepartureFrequency.C
 
+alphatPhaseJayatillekeWallFunction/alphatPhaseJayatillekeWallFunctionFvPatchScalarField.C
 alphatPhaseChangeWallFunction/alphatPhaseChangeWallFunctionFvPatchScalarField.C
-alphatPhaseChangeJayatillekeWallFunction/alphatPhaseChangeJayatillekeWallFunctionFvPatchScalarField.C
+alphatFixedDmdtfWallBoilingWallFunction/alphatFixedDmdtfWallBoilingWallFunctionFvPatchScalarField.C
-alphatFixedDmdtWallBoilingWallFunction/alphatFixedDmdtWallBoilingWallFunctionFvPatchScalarField.C
 alphatWallBoilingWallFunction/alphatWallBoilingWallFunctionFvPatchScalarField.C
 copiedFixedValue/copiedFixedValueFvPatchScalarField.C
 fixedMultiPhaseHeatFlux/fixedMultiPhaseHeatFluxFvPatchScalarField.C

applications/solvers/multiphase/reactingEulerFoam/derivedFvPatchFields/alphatFixedDmdtWallBoilingWallFunction/alphatFixedDmdtWallBoilingWallFunctionFvPatchScalarField.C

@@ -1,215 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2015-2019 OpenFOAM Foundation
-     \\/     M anipulation  |
--------------------------------------------------------------------------------
-License
-    This file is part of OpenFOAM.
-
-    OpenFOAM is free software: you can redistribute it and/or modify it
-    under the terms of the GNU General Public License as published by
-    the Free Software Foundation, either version 3 of the License, or
-    (at your option) any later version.
-
-    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
-    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
-    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
-    for more details.
-
-    You should have received a copy of the GNU General Public License
-    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
-
-\*---------------------------------------------------------------------------*/
-
-#include "alphatFixedDmdtWallBoilingWallFunctionFvPatchScalarField.H"
-#include "fvPatchFieldMapper.H"
-#include "addToRunTimeSelectionTable.H"
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-namespace Foam
-{
-namespace compressible
-{
-
-// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
-
-alphatFixedDmdtWallBoilingWallFunctionFvPatchScalarField::
-alphatFixedDmdtWallBoilingWallFunctionFvPatchScalarField
-(
-    const fvPatch& p,
-    const DimensionedField<scalar, volMesh>& iF
-)
-:
-    alphatPhaseChangeJayatillekeWallFunctionFvPatchScalarField(p, iF),
-    vaporPhaseName_("vapor"),
-    relax_(1.0),
-    fixedDmdt_(0.0),
-    L_(0.0)
-{}
-
-
-alphatFixedDmdtWallBoilingWallFunctionFvPatchScalarField::
-alphatFixedDmdtWallBoilingWallFunctionFvPatchScalarField
-(
-    const fvPatch& p,
-    const DimensionedField<scalar, volMesh>& iF,
-    const dictionary& dict
-)
-:
-    alphatPhaseChangeJayatillekeWallFunctionFvPatchScalarField(p, iF, dict),
-    vaporPhaseName_(dict.lookup("vaporPhase")),
-    relax_(dict.lookupOrDefault<scalar>("relax", 1.0)),
-    fixedDmdt_(dict.lookupOrDefault<scalar>("fixedDmdt", 0.0)),
-    L_(dict.lookupOrDefault<scalar>("L", 0.0))
-{}
-
-
-alphatFixedDmdtWallBoilingWallFunctionFvPatchScalarField::
-alphatFixedDmdtWallBoilingWallFunctionFvPatchScalarField
-(
-    const alphatFixedDmdtWallBoilingWallFunctionFvPatchScalarField& psf,
-    const fvPatch& p,
-    const DimensionedField<scalar, volMesh>& iF,
-    const fvPatchFieldMapper& mapper
-)
-:
-    alphatPhaseChangeJayatillekeWallFunctionFvPatchScalarField
-    (
-        psf,
-        p,
-        iF,
-        mapper
-    ),
-    fixedDmdt_(psf.fixedDmdt_),
-    L_(psf.L_)
-{}
-
-
-alphatFixedDmdtWallBoilingWallFunctionFvPatchScalarField::
-alphatFixedDmdtWallBoilingWallFunctionFvPatchScalarField
-(
-    const alphatFixedDmdtWallBoilingWallFunctionFvPatchScalarField& psf
-)
-:
-    alphatPhaseChangeJayatillekeWallFunctionFvPatchScalarField(psf),
-    relax_(psf.relax_),
-    fixedDmdt_(psf.fixedDmdt_),
-    L_(psf.L_)
-{}
-
-
-alphatFixedDmdtWallBoilingWallFunctionFvPatchScalarField::
-alphatFixedDmdtWallBoilingWallFunctionFvPatchScalarField
-(
-    const alphatFixedDmdtWallBoilingWallFunctionFvPatchScalarField& psf,
-    const DimensionedField<scalar, volMesh>& iF
-)
-:
-    alphatPhaseChangeJayatillekeWallFunctionFvPatchScalarField(psf, iF),
-    relax_(psf.relax_),
-    fixedDmdt_(psf.fixedDmdt_),
-    L_(psf.L_)
-{}
-
-
-// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
-
-bool alphatFixedDmdtWallBoilingWallFunctionFvPatchScalarField::
-activePhasePair(const phasePairKey& phasePair) const
-{
-    if (phasePair == phasePairKey(vaporPhaseName_, internalField().group()))
-    {
-        return true;
-    }
-    else
-    {
-        return false;
-    }
-}
-
-
-const scalarField& alphatFixedDmdtWallBoilingWallFunctionFvPatchScalarField::
-dmdt(const phasePairKey& phasePair) const
-{
-    if (activePhasePair(phasePair))
-    {
-        return dmdt_;
-    }
-    else
-    {
-        FatalErrorInFunction
-            << " dmdt requested for invalid phasePair!"
-            << abort(FatalError);
-
-        return mDotL_;
-    }
-}
-
-
-const scalarField& alphatFixedDmdtWallBoilingWallFunctionFvPatchScalarField::
-mDotL(const phasePairKey& phasePair) const
-{
-    if (activePhasePair(phasePair))
-    {
-        return mDotL_;
-    }
-    else
-    {
-        FatalErrorInFunction
-            << " mDotL requested for invalid phasePair!"
-            << abort(FatalError);
-
-        return mDotL_;
-    }
-}
-
-
-void alphatFixedDmdtWallBoilingWallFunctionFvPatchScalarField::updateCoeffs()
-{
-    if (updated())
-    {
-        return;
-    }
-
-    dmdt_ = (1 - relax_)*dmdt_ + relax_*fixedDmdt_;
-    mDotL_ = dmdt_*L_;
-
-    operator==(calcAlphat(*this));
-
-    fixedValueFvPatchScalarField::updateCoeffs();
-}
-
-
-void alphatFixedDmdtWallBoilingWallFunctionFvPatchScalarField::write
-(
-    Ostream& os
-) const
-{
-    fvPatchField<scalar>::write(os);
-    writeEntry(os, "vaporPhase", vaporPhaseName_);
-    writeEntry(os, "relax", relax_);
-    writeEntry(os, "fixedDmdt", fixedDmdt_);
-    writeEntry(os, "L", L_);
-    writeEntry(os, "dmdt", dmdt_);
-    writeEntry(os, "value", *this);
-}
-
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-makePatchTypeField
-(
-    fvPatchScalarField,
-    alphatFixedDmdtWallBoilingWallFunctionFvPatchScalarField
-);
-
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-} // End namespace compressible
-} // End namespace Foam
-
-// ************************************************************************* //

applications/solvers/multiphase/reactingEulerFoam/derivedFvPatchFields/alphatFixedDmdtfWallBoilingWallFunction/alphatFixedDmdtfWallBoilingWallFunctionFvPatchScalarField.C

@@ -0,0 +1,156 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     | Website:  https://openfoam.org
+    \\  /    A nd           | Copyright (C) 2015-2019 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "alphatFixedDmdtfWallBoilingWallFunctionFvPatchScalarField.H"
+#include "fvPatchFieldMapper.H"
+#include "addToRunTimeSelectionTable.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace compressible
+{
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+alphatFixedDmdtfWallBoilingWallFunctionFvPatchScalarField::
+alphatFixedDmdtfWallBoilingWallFunctionFvPatchScalarField
+(
+    const fvPatch& p,
+    const DimensionedField<scalar, volMesh>& iF
+)
+:
+    alphatPhaseChangeWallFunctionFvPatchScalarField(p, iF),
+    fixedDmdtf_(0),
+    L_(0)
+{}
+
+
+alphatFixedDmdtfWallBoilingWallFunctionFvPatchScalarField::
+alphatFixedDmdtfWallBoilingWallFunctionFvPatchScalarField
+(
+    const fvPatch& p,
+    const DimensionedField<scalar, volMesh>& iF,
+    const dictionary& dict
+)
+:
+    alphatPhaseChangeWallFunctionFvPatchScalarField(p, iF, dict),
+    fixedDmdtf_(dict.lookupOrDefault<scalar>("fixedDmdtf", 0)),
+    L_(dict.lookupOrDefault<scalar>("L", 0))
+{}
+
+
+alphatFixedDmdtfWallBoilingWallFunctionFvPatchScalarField::
+alphatFixedDmdtfWallBoilingWallFunctionFvPatchScalarField
+(
+    const alphatFixedDmdtfWallBoilingWallFunctionFvPatchScalarField& psf,
+    const fvPatch& p,
+    const DimensionedField<scalar, volMesh>& iF,
+    const fvPatchFieldMapper& mapper
+)
+:
+    alphatPhaseChangeWallFunctionFvPatchScalarField
+    (
+        psf,
+        p,
+        iF,
+        mapper
+    ),
+    fixedDmdtf_(psf.fixedDmdtf_),
+    L_(psf.L_)
+{}
+
+
+alphatFixedDmdtfWallBoilingWallFunctionFvPatchScalarField::
+alphatFixedDmdtfWallBoilingWallFunctionFvPatchScalarField
+(
+    const alphatFixedDmdtfWallBoilingWallFunctionFvPatchScalarField& psf
+)
+:
+    alphatPhaseChangeWallFunctionFvPatchScalarField(psf),
+    fixedDmdtf_(psf.fixedDmdtf_),
+    L_(psf.L_)
+{}
+
+
+alphatFixedDmdtfWallBoilingWallFunctionFvPatchScalarField::
+alphatFixedDmdtfWallBoilingWallFunctionFvPatchScalarField
+(
+    const alphatFixedDmdtfWallBoilingWallFunctionFvPatchScalarField& psf,
+    const DimensionedField<scalar, volMesh>& iF
+)
+:
+    alphatPhaseChangeWallFunctionFvPatchScalarField(psf, iF),
+    fixedDmdtf_(psf.fixedDmdtf_),
+    L_(psf.L_)
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+void alphatFixedDmdtfWallBoilingWallFunctionFvPatchScalarField::updateCoeffs()
+{
+    if (updated())
+    {
+        return;
+    }
+
+    dmdtf_ = (1 - relax_)*dmdtf_ + relax_*fixedDmdtf_;
+    dmdtLf_ = dmdtf_*L_;
+
+    operator==(calcAlphat(*this));
+
+    fixedValueFvPatchScalarField::updateCoeffs();
+}
+
+
+void alphatFixedDmdtfWallBoilingWallFunctionFvPatchScalarField::write
+(
+    Ostream& os
+) const
+{
+    alphatPhaseChangeWallFunctionFvPatchScalarField::write(os);
+
+    writeEntry(os, "fixedDmdtf", fixedDmdtf_);
+    writeEntry(os, "L", L_);
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+makePatchTypeField
+(
+    fvPatchScalarField,
+    alphatFixedDmdtfWallBoilingWallFunctionFvPatchScalarField
+);
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace compressible
+} // End namespace Foam
+
+// ************************************************************************* //

applications/solvers/multiphase/reactingEulerFoam/derivedFvPatchFields/alphatFixedDmdtfWallBoilingWallFunction/alphatFixedDmdtfWallBoilingWallFunctionFvPatchScalarField.H

@@ -23,25 +23,24 @@ License
 
 Class
     Foam::compressible::
-        alphatFixedDmdtWallBoilingWallFunctionFvPatchScalarField
+        alphatFixedDmdtfWallBoilingWallFunctionFvPatchScalarField
 
 Description
-    A simple alphatPhaseChangeJayatillekeWallFunctionFvPatchScalarField with
+    A simple alphatPhaseChangeWallFunctionFvPatchScalarField with
     a fixed volumetric phase-change mass flux.
 
 See also
-    Foam::compressible::
+    Foam::compressible::alphatPhaseChangeWallFunctionFvPatchScalarField
-        alphatPhaseChangeJayatillekeWallFunctionFvPatchScalarField
 
 SourceFiles
-    alphatFixedDmdtWallBoilingWallFunctionFvPatchScalarField.C
+    alphatFixedDmdtfWallBoilingWallFunctionFvPatchScalarField.C
 
 \*---------------------------------------------------------------------------*/
 
-#ifndef alphatFixedDmdtWallBoilingWallFunctionFvPatchScalarField_H
+#ifndef alphatFixedDmdtfWallBoilingWallFunctionFvPatchScalarField_H
-#define alphatFixedDmdtWallBoilingWallFunctionFvPatchScalarField_H
+#define alphatFixedDmdtfWallBoilingWallFunctionFvPatchScalarField_H
 
-#include "alphatPhaseChangeJayatillekeWallFunctionFvPatchScalarField.H"
+#include "alphatPhaseChangeWallFunctionFvPatchScalarField.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -51,23 +50,17 @@ namespace compressible
 {
 
 /*---------------------------------------------------------------------------*\
-   Class alphatFixedDmdtWallBoilingWallFunctionFvPatchScalarField Declaration
+   Class alphatFixedDmdtfWallBoilingWallFunctionFvPatchScalarField Declaration
 \*---------------------------------------------------------------------------*/
 
-class alphatFixedDmdtWallBoilingWallFunctionFvPatchScalarField
+class alphatFixedDmdtfWallBoilingWallFunctionFvPatchScalarField
 :
-    public alphatPhaseChangeJayatillekeWallFunctionFvPatchScalarField
+    public alphatPhaseChangeWallFunctionFvPatchScalarField
 {
     // Private Data
 
-        //- name on the phase
-        word vaporPhaseName_;
-
-        //- dmdt relaxationFactor
-        scalar relax_;
-
         //- Volumetric phase-change mass flux in near wall cells
-        scalar fixedDmdt_;
+        scalar fixedDmdtf_;
 
         //- Latent heat
         scalar L_;
@@ -76,20 +69,20 @@ class alphatFixedDmdtWallBoilingWallFunctionFvPatchScalarField
 public:
 
     //- Runtime type information
-    TypeName("compressible::alphatFixedDmdtWallBoilingWallFunction");
+    TypeName("compressible::alphatFixedDmdtfWallBoilingWallFunction");
 
 
     // Constructors
 
         //- Construct from patch and internal field
-        alphatFixedDmdtWallBoilingWallFunctionFvPatchScalarField
+        alphatFixedDmdtfWallBoilingWallFunctionFvPatchScalarField
         (
             const fvPatch&,
             const DimensionedField<scalar, volMesh>&
         );
 
         //- Construct from patch, internal field and dictionary
-        alphatFixedDmdtWallBoilingWallFunctionFvPatchScalarField
+        alphatFixedDmdtfWallBoilingWallFunctionFvPatchScalarField
         (
             const fvPatch&,
             const DimensionedField<scalar, volMesh>&,
@@ -97,20 +90,20 @@ public:
         );
 
         //- Construct by mapping given
-        //  alphatFixedDmdtWallBoilingWallFunctionFvPatchScalarField
+        //  alphatFixedDmdtfWallBoilingWallFunctionFvPatchScalarField
         //  onto a new patch
-        alphatFixedDmdtWallBoilingWallFunctionFvPatchScalarField
+        alphatFixedDmdtfWallBoilingWallFunctionFvPatchScalarField
         (
-            const alphatFixedDmdtWallBoilingWallFunctionFvPatchScalarField&,
+            const alphatFixedDmdtfWallBoilingWallFunctionFvPatchScalarField&,
             const fvPatch&,
             const DimensionedField<scalar, volMesh>&,
             const fvPatchFieldMapper&
         );
 
         //- Copy constructor
-        alphatFixedDmdtWallBoilingWallFunctionFvPatchScalarField
+        alphatFixedDmdtfWallBoilingWallFunctionFvPatchScalarField
         (
-            const alphatFixedDmdtWallBoilingWallFunctionFvPatchScalarField&
+            const alphatFixedDmdtfWallBoilingWallFunctionFvPatchScalarField&
         );
 
         //- Construct and return a clone
@@ -118,7 +111,7 @@ public:
         {
             return tmp<fvPatchScalarField>
             (
-                new alphatFixedDmdtWallBoilingWallFunctionFvPatchScalarField
+                new alphatFixedDmdtfWallBoilingWallFunctionFvPatchScalarField
                 (
                     *this
                 )
@@ -126,9 +119,9 @@ public:
         }
 
         //- Copy constructor setting internal field reference
-        alphatFixedDmdtWallBoilingWallFunctionFvPatchScalarField
+        alphatFixedDmdtfWallBoilingWallFunctionFvPatchScalarField
         (
-            const alphatFixedDmdtWallBoilingWallFunctionFvPatchScalarField&,
+            const alphatFixedDmdtfWallBoilingWallFunctionFvPatchScalarField&,
             const DimensionedField<scalar, volMesh>&
         );
 
@@ -140,7 +133,7 @@ public:
         {
             return tmp<fvPatchScalarField>
             (
-                new alphatFixedDmdtWallBoilingWallFunctionFvPatchScalarField
+                new alphatFixedDmdtfWallBoilingWallFunctionFvPatchScalarField
                 (
                     *this,
                     iF
@@ -151,15 +144,6 @@ public:
 
     // Member Functions
 
-        //- Is there phase change mass transfer for this phasePair
-        virtual bool activePhasePair(const phasePairKey&) const;
-
-        //- Return the rate of phase-change for specific phase pair
-        virtual const scalarField& dmdt(const phasePairKey&) const;
-
-        //- Return the rate of phase-change for specific phase pair
-        virtual const scalarField& mDotL(const phasePairKey&) const;
-
         // Evaluation functions
 
             //- Update the coefficients associated with the patch field

applications/solvers/multiphase/reactingEulerFoam/derivedFvPatchFields/alphatPhaseChangeWallFunction/alphatPhaseChangeWallFunctionFvPatchScalarField.C

@@ -48,9 +48,11 @@ alphatPhaseChangeWallFunctionFvPatchScalarField
     const DimensionedField<scalar, volMesh>& iF
 )
 :
-    fixedValueFvPatchScalarField(p, iF),
+    alphatPhaseJayatillekeWallFunctionFvPatchScalarField(p, iF),
-    dmdt_(p.size(), 0),
+    otherPhaseName_(word::null),
-    mDotL_(p.size(), 0)
+    relax_(1),
+    dmdtf_(p.size(), 0),
+    dmdtLf_(p.size(), 0)
 {}
 
 
@@ -62,18 +64,31 @@ alphatPhaseChangeWallFunctionFvPatchScalarField
     const dictionary& dict
 )
 :
-    fixedValueFvPatchScalarField(p, iF, dict),
+    alphatPhaseJayatillekeWallFunctionFvPatchScalarField(p, iF, dict),
-    dmdt_(p.size(), 0),
+    otherPhaseName_(dict.lookup("otherPhase")),
-    mDotL_(p.size(), 0)
+    relax_(dict.lookupOrDefault<scalar>("relax", 1)),
+    dmdtf_(p.size(), 0),
+    dmdtLf_(p.size(), 0)
 {
-    if (dict.found("dmdt"))
+    // Check that otherPhaseName != this phase
+    if (internalField().group() == otherPhaseName_)
     {
-        dmdt_ = scalarField("dmdt", dict, p.size());
+        FatalErrorInFunction
+            << "otherPhase should be the name of the vapor phase that "
+            << "corresponds to the liquid base or vice versa" << nl
+            << "This phase: " << internalField().group() << nl
+            << "otherPhase: " << otherPhaseName_
+            << abort(FatalError);
     }
 
-    if (dict.found("mDotL"))
+    if (dict.found("dmdtf"))
     {
-        dmdt_ = scalarField("mDotL", dict, p.size());
+        dmdtf_ = scalarField("dmdtf", dict, p.size());
+    }
+
+    if (dict.found("dmdtLf"))
+    {
+        dmdtf_ = scalarField("dmdtLf", dict, p.size());
     }
 }
 
@@ -87,9 +102,11 @@ alphatPhaseChangeWallFunctionFvPatchScalarField
     const fvPatchFieldMapper& mapper
 )
 :
-    fixedValueFvPatchScalarField(ptf, p, iF, mapper),
+    alphatPhaseJayatillekeWallFunctionFvPatchScalarField(ptf, p, iF, mapper),
-    dmdt_(mapper(ptf.dmdt_)),
+    otherPhaseName_(ptf.otherPhaseName_),
-    mDotL_(mapper(ptf.mDotL_))
+    relax_(ptf.relax_),
+    dmdtf_(mapper(ptf.dmdtf_)),
+    dmdtLf_(mapper(ptf.dmdtLf_))
 {}
 
 
@@ -99,9 +116,11 @@ alphatPhaseChangeWallFunctionFvPatchScalarField
     const alphatPhaseChangeWallFunctionFvPatchScalarField& awfpsf
 )
 :
-    fixedValueFvPatchScalarField(awfpsf),
+    alphatPhaseJayatillekeWallFunctionFvPatchScalarField(awfpsf),
-    dmdt_(awfpsf.dmdt_),
+    otherPhaseName_(awfpsf.otherPhaseName_),
-    mDotL_(awfpsf.mDotL_)
+    relax_(awfpsf.relax_),
+    dmdtf_(awfpsf.dmdtf_),
+    dmdtLf_(awfpsf.dmdtLf_)
 {}
 
 
@@ -112,21 +131,88 @@ alphatPhaseChangeWallFunctionFvPatchScalarField
     const DimensionedField<scalar, volMesh>& iF
 )
 :
-    fixedValueFvPatchScalarField(awfpsf, iF),
+    alphatPhaseJayatillekeWallFunctionFvPatchScalarField(awfpsf, iF),
-    dmdt_(awfpsf.dmdt_),
+    otherPhaseName_(awfpsf.otherPhaseName_),
-    mDotL_(awfpsf.mDotL_)
+    relax_(awfpsf.relax_),
+    dmdtf_(awfpsf.dmdtf_),
+    dmdtLf_(awfpsf.dmdtLf_)
 {}
 
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
-void alphatPhaseChangeWallFunctionFvPatchScalarField::
+bool alphatPhaseChangeWallFunctionFvPatchScalarField::
-write(Ostream& os) const
+activePhasePair(const phasePairKey& phasePair) const
 {
-    fvPatchField<scalar>::write(os);
+    if (phasePair == phasePairKey(otherPhaseName_, internalField().group()))
-    writeEntry(os, "dmdt", dmdt_);
+    {
-    writeEntry(os, "mDotL", mDotL_);
+        return true;
-    writeEntry(os, "value", *this);
+    }
+    else
+    {
+        return false;
+    }
+}
+
+
+const scalarField&
+alphatPhaseChangeWallFunctionFvPatchScalarField::dmdtf() const
+{
+    return dmdtf_;
+}
+
+
+const scalarField& alphatPhaseChangeWallFunctionFvPatchScalarField::
+dmdtf(const phasePairKey& phasePair) const
+{
+    if (activePhasePair(phasePair))
+    {
+        return dmdtf_;
+    }
+    else
+    {
+        FatalErrorInFunction
+            << " dmdtf requested for invalid phasePair!"
+            << abort(FatalError);
+
+        return dmdtLf_;
+    }
+}
+
+
+const scalarField&
+alphatPhaseChangeWallFunctionFvPatchScalarField::dmdtLf() const
+{
+    return dmdtLf_;
+}
+
+
+const scalarField& alphatPhaseChangeWallFunctionFvPatchScalarField::
+dmdtLf(const phasePairKey& phasePair) const
+{
+    if (activePhasePair(phasePair))
+    {
+        return dmdtLf_;
+    }
+    else
+    {
+        FatalErrorInFunction
+            << " dmdtLf requested for invalid phasePair!"
+            << abort(FatalError);
+
+        return dmdtLf_;
+    }
+}
+
+
+void alphatPhaseChangeWallFunctionFvPatchScalarField::write(Ostream& os) const
+{
+    alphatPhaseJayatillekeWallFunctionFvPatchScalarField::write(os);
+
+    writeEntry(os, "otherPhase", otherPhaseName_);
+    writeEntry(os, "relax", relax_);
+    writeEntry(os, "dmdtf", dmdtf_);
+    writeEntry(os, "dmdtLf", dmdtLf_);
 }
 
 

applications/solvers/multiphase/reactingEulerFoam/derivedFvPatchFields/alphatPhaseChangeWallFunction/alphatPhaseChangeWallFunctionFvPatchScalarField.H

@@ -28,8 +28,7 @@ Description
     Abstract base-class for all alphatWallFunctions supporting phase-change.
 
 See also
-    Foam::fixedValueFvPatchScalarField
+    Foam::alphatPhaseJayatillekeWallFunctionFvPatchScalarField
-    Foam::alphatWallFunctionFvPatchScalarField
 
 SourceFiles
     alphatPhaseChangeWallFunctionFvPatchScalarField.C
@@ -39,7 +38,7 @@ SourceFiles
 #ifndef alphatPhaseChangeWallFunctionFvPatchScalarField_H
 #define alphatPhaseChangeWallFunctionFvPatchScalarField_H
 
-#include "fixedValueFvPatchFields.H"
+#include "alphatPhaseJayatillekeWallFunctionFvPatchScalarField.H"
 #include "phasePairKey.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@@ -55,17 +54,23 @@ namespace compressible
 
 class alphatPhaseChangeWallFunctionFvPatchScalarField
 :
-    public fixedValueFvPatchScalarField
+    public alphatPhaseJayatillekeWallFunctionFvPatchScalarField
 {
 protected:
 
     // Protected data
 
+        //- Name of the other phase
+        const word otherPhaseName_;
+
+        //- Relaxation factor
+        const scalar relax_;
+
         //- Rate of phase-change
-        scalarField dmdt_;
+        scalarField dmdtf_;
 
         //- Latent heat of the phase-change
-        scalarField mDotL_;
+        scalarField dmdtLf_;
 
 
 public:
@@ -118,47 +123,22 @@ public:
 
     // Member Functions
 
-        //- Return the rate of phase-change
-        virtual const scalarField& dmdt() const
-        {
-            return dmdt_;
-        }
-
-        //- Return the enthalpy source due to phase-change
-        virtual const scalarField& mDotL() const
-        {
-            return mDotL_;
-        }
 
         //- Is there phase change mass transfer for this phasePair
-        virtual bool activePhasePair(const phasePairKey&) const
+        bool activePhasePair(const phasePairKey&) const;
-        {
+
-            return false;
+        //- Return the rate of phase-change
-        }
+        const scalarField& dmdtf() const;
 
         //- Return the rate of phase-change for specific phase pair
-        virtual const scalarField& dmdt(const phasePairKey&) const
+        const scalarField& dmdtf(const phasePairKey&) const;
-        {
+
-            return dmdt_;
+        //- Return the enthalpy source due to phase-change
-        }
+        const scalarField& dmdtLf() const;
 
         //- Return the rate of phase-change for specific phase pair
-        virtual const scalarField& mDotL(const phasePairKey&) const
+        const scalarField& dmdtLf(const phasePairKey&) const;
-        {
-            return mDotL_;
-        }
 
-        //- Return the rate of phase-change for specific phase
-        virtual scalarField dmdt(const word&) const
-        {
-            return dmdt_;
-        }
-
-        //- Return the enthalpy source due to phase-change for specific phase
-        virtual scalarField mDotL(const word&) const
-        {
-            return mDotL_;
-        }
 
         // Evaluation functions
 

applications/solvers/multiphase/reactingEulerFoam/derivedFvPatchFields/alphatPhaseJayatillekeWallFunction/alphatPhaseJayatillekeWallFunctionFvPatchScalarField.C

@@ -23,11 +23,9 @@ License
 
 \*---------------------------------------------------------------------------*/
 
-#include "alphatPhaseChangeJayatillekeWallFunctionFvPatchScalarField.H"
+#include "alphatPhaseJayatillekeWallFunctionFvPatchScalarField.H"
 #include "phaseSystem.H"
-#include "compressibleTurbulenceModel.H"
+#include "phaseCompressibleTurbulenceModel.H"
-#include "ThermalDiffusivity.H"
-#include "PhaseCompressibleTurbulenceModel.H"
 #include "addToRunTimeSelectionTable.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@@ -39,27 +37,23 @@ namespace compressible
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
-scalar alphatPhaseChangeJayatillekeWallFunctionFvPatchScalarField::maxExp_
+scalar alphatPhaseJayatillekeWallFunctionFvPatchScalarField::maxExp_  = 50.0;
-    = 50.0;
+scalar alphatPhaseJayatillekeWallFunctionFvPatchScalarField::tolerance_ = 0.01;
-scalar alphatPhaseChangeJayatillekeWallFunctionFvPatchScalarField::tolerance_
+label alphatPhaseJayatillekeWallFunctionFvPatchScalarField::maxIters_ = 10;
-    = 0.01;
-label alphatPhaseChangeJayatillekeWallFunctionFvPatchScalarField::maxIters_
-    = 10;
 
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
-tmp<scalarField>
+tmp<scalarField> alphatPhaseJayatillekeWallFunctionFvPatchScalarField::Psmooth
-alphatPhaseChangeJayatillekeWallFunctionFvPatchScalarField::Psmooth
 (
     const scalarField& Prat
 ) const
 {
-    return 9.24*(pow(Prat, 0.75) - 1)*(1 + 0.28*exp(-0.007*Prat));
+    return 9.24*(pow(Prat, 0.75) - 1.0)*(1.0 + 0.28*exp(-0.007*Prat));
 }
 
 
 tmp<scalarField>
-alphatPhaseChangeJayatillekeWallFunctionFvPatchScalarField::yPlusTherm
+alphatPhaseJayatillekeWallFunctionFvPatchScalarField::yPlusTherm
 (
     const nutWallFunctionFvPatchScalarField& nutw,
     const scalarField& P,
@@ -83,6 +77,7 @@ alphatPhaseChangeJayatillekeWallFunctionFvPatchScalarField::yPlusTherm
             if (yptNew < vSmall)
             {
                 ypsf[facei] = 0;
+                break;
             }
             else if (mag(yptNew - ypt) < tolerance_)
             {
@@ -100,8 +95,76 @@ alphatPhaseChangeJayatillekeWallFunctionFvPatchScalarField::yPlusTherm
     return typsf;
 }
 
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+alphatPhaseJayatillekeWallFunctionFvPatchScalarField::
+alphatPhaseJayatillekeWallFunctionFvPatchScalarField
+(
+    const fvPatch& p,
+    const DimensionedField<scalar, volMesh>& iF
+)
+:
+    fixedValueFvPatchScalarField(p, iF),
+    Prt_(0.85)
+{}
+
+
+alphatPhaseJayatillekeWallFunctionFvPatchScalarField::
+alphatPhaseJayatillekeWallFunctionFvPatchScalarField
+(
+    const alphatPhaseJayatillekeWallFunctionFvPatchScalarField& ptf,
+    const fvPatch& p,
+    const DimensionedField<scalar, volMesh>& iF,
+    const fvPatchFieldMapper& mapper
+)
+:
+    fixedValueFvPatchScalarField(ptf, p, iF, mapper),
+    Prt_(ptf.Prt_)
+{}
+
+
+alphatPhaseJayatillekeWallFunctionFvPatchScalarField::
+alphatPhaseJayatillekeWallFunctionFvPatchScalarField
+(
+    const fvPatch& p,
+    const DimensionedField<scalar, volMesh>& iF,
+    const dictionary& dict
+)
+:
+    fixedValueFvPatchScalarField(p, iF, dict),
+    Prt_(dict.lookupOrDefault<scalar>("Prt", 0.85))
+{}
+
+
+
+alphatPhaseJayatillekeWallFunctionFvPatchScalarField::
+alphatPhaseJayatillekeWallFunctionFvPatchScalarField
+(
+    const alphatPhaseJayatillekeWallFunctionFvPatchScalarField& awfpsf
+)
+:
+    fixedValueFvPatchScalarField(awfpsf),
+    Prt_(awfpsf.Prt_)
+{}
+
+
+alphatPhaseJayatillekeWallFunctionFvPatchScalarField::
+alphatPhaseJayatillekeWallFunctionFvPatchScalarField
+(
+    const alphatPhaseJayatillekeWallFunctionFvPatchScalarField& awfpsf,
+    const DimensionedField<scalar, volMesh>& iF
+)
+:
+    fixedValueFvPatchScalarField(awfpsf, iF),
+    Prt_(awfpsf.Prt_)
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
 tmp<scalarField>
-alphatPhaseChangeJayatillekeWallFunctionFvPatchScalarField::calcAlphat
+alphatPhaseJayatillekeWallFunctionFvPatchScalarField::calcAlphat
 (
     const scalarField& prevAlphat
 ) const
@@ -221,73 +284,7 @@ alphatPhaseChangeJayatillekeWallFunctionFvPatchScalarField::calcAlphat
 }
 
 
-// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+void alphatPhaseJayatillekeWallFunctionFvPatchScalarField::updateCoeffs()
-
-alphatPhaseChangeJayatillekeWallFunctionFvPatchScalarField::
-alphatPhaseChangeJayatillekeWallFunctionFvPatchScalarField
-(
-    const fvPatch& p,
-    const DimensionedField<scalar, volMesh>& iF
-)
-:
-    alphatPhaseChangeWallFunctionFvPatchScalarField(p, iF),
-    Prt_(0.85)
-{}
-
-
-alphatPhaseChangeJayatillekeWallFunctionFvPatchScalarField::
-alphatPhaseChangeJayatillekeWallFunctionFvPatchScalarField
-(
-    const fvPatch& p,
-    const DimensionedField<scalar, volMesh>& iF,
-    const dictionary& dict
-)
-:
-    alphatPhaseChangeWallFunctionFvPatchScalarField(p, iF, dict),
-    Prt_(dict.lookupOrDefault<scalar>("Prt", 0.85))
-{}
-
-
-alphatPhaseChangeJayatillekeWallFunctionFvPatchScalarField::
-alphatPhaseChangeJayatillekeWallFunctionFvPatchScalarField
-(
-    const alphatPhaseChangeJayatillekeWallFunctionFvPatchScalarField& ptf,
-    const fvPatch& p,
-    const DimensionedField<scalar, volMesh>& iF,
-    const fvPatchFieldMapper& mapper
-)
-:
-    alphatPhaseChangeWallFunctionFvPatchScalarField(ptf, p, iF, mapper),
-    Prt_(ptf.Prt_)
-{}
-
-
-alphatPhaseChangeJayatillekeWallFunctionFvPatchScalarField::
-alphatPhaseChangeJayatillekeWallFunctionFvPatchScalarField
-(
-    const alphatPhaseChangeJayatillekeWallFunctionFvPatchScalarField& awfpsf
-)
-:
-    alphatPhaseChangeWallFunctionFvPatchScalarField(awfpsf),
-    Prt_(awfpsf.Prt_)
-{}
-
-
-alphatPhaseChangeJayatillekeWallFunctionFvPatchScalarField::
-alphatPhaseChangeJayatillekeWallFunctionFvPatchScalarField
-(
-    const alphatPhaseChangeJayatillekeWallFunctionFvPatchScalarField& awfpsf,
-    const DimensionedField<scalar, volMesh>& iF
-)
-:
-    alphatPhaseChangeWallFunctionFvPatchScalarField(awfpsf, iF),
-    Prt_(awfpsf.Prt_)
-{}
-
-
-// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
-
-void alphatPhaseChangeJayatillekeWallFunctionFvPatchScalarField::updateCoeffs()
 {
     if (updated())
     {
@@ -300,14 +297,13 @@ void alphatPhaseChangeJayatillekeWallFunctionFvPatchScalarField::updateCoeffs()
 }
 
 
-void alphatPhaseChangeJayatillekeWallFunctionFvPatchScalarField::write
+void alphatPhaseJayatillekeWallFunctionFvPatchScalarField::write
 (
     Ostream& os
 ) const
 {
     fvPatchField<scalar>::write(os);
     writeEntry(os, "Prt", Prt_);
-    writeEntry(os, "dmdt", dmdt_);
     writeEntry(os, "value", *this);
 }
 
@@ -317,10 +313,9 @@ void alphatPhaseChangeJayatillekeWallFunctionFvPatchScalarField::write
 makePatchTypeField
 (
     fvPatchScalarField,
-    alphatPhaseChangeJayatillekeWallFunctionFvPatchScalarField
+    alphatPhaseJayatillekeWallFunctionFvPatchScalarField
 );
 
-
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 } // End namespace compressible

applications/solvers/multiphase/reactingEulerFoam/derivedFvPatchFields/alphatPhaseJayatillekeWallFunction/alphatPhaseJayatillekeWallFunctionFvPatchScalarField.H

@@ -23,7 +23,7 @@ License
 
 Class
     Foam::compressible::
-        alphatPhaseChangeJayatillekeWallFunctionFvPatchScalarField
+        alphatPhaseJayatillekeWallFunctionFvPatchScalarField
 
 Description
     This boundary condition provides a thermal wall function for turbulent
@@ -43,7 +43,7 @@ Usage
     \verbatim
     <patchName>
     {
-        type            alphatPhaseChangeJayatillekeWallFunction;
+        type            alphatPhaseJayatillekeWallFunction;
         Prt             0.85;
         kappa           0.41;
         E               9.8;
@@ -55,14 +55,14 @@ See also
     Foam::compressible::alphatPhaseChangeWallFunctionFvPatchScalarField
 
 SourceFiles
-    alphatPhaseChangeJayatillekeWallFunctionFvPatchScalarField.C
+    alphatPhaseJayatillekeWallFunctionFvPatchScalarField.C
 
 \*---------------------------------------------------------------------------*/
 
-#ifndef compressible_alphatPhaseChangeJayatillekeWallFunctionFvPatchScalarField_H
+#ifndef compressible_alphatPhaseJayatillekeWallFunctionFvPatchScalarField_H
-#define compressible_alphatPhaseChangeJayatillekeWallFunctionFvPatchScalarField_H
+#define compressible_alphatPhaseJayatillekeWallFunctionFvPatchScalarField_H
 
-#include "alphatPhaseChangeWallFunctionFvPatchScalarField.H"
+#include "fixedValueFvPatchFields.H"
 #include "nutWallFunctionFvPatchScalarField.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@@ -73,12 +73,12 @@ namespace compressible
 {
 
 /*---------------------------------------------------------------------------*\
- Class alphatPhaseChangeJayatillekeWallFunctionFvPatchScalarField Declaration
+ Class alphatPhaseJayatillekeWallFunctionFvPatchScalarField Declaration
 \*---------------------------------------------------------------------------*/
 
-class alphatPhaseChangeJayatillekeWallFunctionFvPatchScalarField
+class alphatPhaseJayatillekeWallFunctionFvPatchScalarField
 :
-    public alphatPhaseChangeWallFunctionFvPatchScalarField
+    public fixedValueFvPatchScalarField
 {
 
 protected:
@@ -108,30 +108,24 @@ protected:
             const scalarField& Prat
         ) const;
 
-        //- Update turbulent thermal diffusivity
-        tmp<scalarField> calcAlphat
-        (
-            const scalarField& prevAlphat
-        ) const;
-
 
 public:
 
     //- Runtime type information
-    TypeName("compressible::alphatPhaseChangeJayatillekeWallFunction");
+    TypeName("compressible::alphatPhaseJayatillekeWallFunction");
 
 
     // Constructors
 
         //- Construct from patch and internal field
-        alphatPhaseChangeJayatillekeWallFunctionFvPatchScalarField
+        alphatPhaseJayatillekeWallFunctionFvPatchScalarField
         (
             const fvPatch&,
             const DimensionedField<scalar, volMesh>&
         );
 
         //- Construct from patch, internal field and dictionary
-        alphatPhaseChangeJayatillekeWallFunctionFvPatchScalarField
+        alphatPhaseJayatillekeWallFunctionFvPatchScalarField
         (
             const fvPatch&,
             const DimensionedField<scalar, volMesh>&,
@@ -139,20 +133,20 @@ public:
         );
 
         //- Construct by mapping given
-        //  alphatPhaseChangeJayatillekeWallFunctionFvPatchScalarField
+        //  alphatPhaseJayatillekeWallFunctionFvPatchScalarField
         //  onto a new patch
-        alphatPhaseChangeJayatillekeWallFunctionFvPatchScalarField
+        alphatPhaseJayatillekeWallFunctionFvPatchScalarField
         (
-            const alphatPhaseChangeJayatillekeWallFunctionFvPatchScalarField&,
+            const alphatPhaseJayatillekeWallFunctionFvPatchScalarField&,
             const fvPatch&,
             const DimensionedField<scalar, volMesh>&,
             const fvPatchFieldMapper&
         );
 
         //- Copy constructor
-        alphatPhaseChangeJayatillekeWallFunctionFvPatchScalarField
+        alphatPhaseJayatillekeWallFunctionFvPatchScalarField
         (
-            const alphatPhaseChangeJayatillekeWallFunctionFvPatchScalarField&
+            const alphatPhaseJayatillekeWallFunctionFvPatchScalarField&
         );
 
         //- Construct and return a clone
@@ -160,7 +154,7 @@ public:
         {
             return tmp<fvPatchScalarField>
             (
-                new alphatPhaseChangeJayatillekeWallFunctionFvPatchScalarField
+                new alphatPhaseJayatillekeWallFunctionFvPatchScalarField
                 (
                     *this
                 )
@@ -168,9 +162,9 @@ public:
         }
 
         //- Copy constructor setting internal field reference
-        alphatPhaseChangeJayatillekeWallFunctionFvPatchScalarField
+        alphatPhaseJayatillekeWallFunctionFvPatchScalarField
         (
-            const alphatPhaseChangeJayatillekeWallFunctionFvPatchScalarField&,
+            const alphatPhaseJayatillekeWallFunctionFvPatchScalarField&,
             const DimensionedField<scalar, volMesh>&
         );
 
@@ -182,7 +176,7 @@ public:
         {
             return tmp<fvPatchScalarField>
             (
-                new alphatPhaseChangeJayatillekeWallFunctionFvPatchScalarField
+                new alphatPhaseJayatillekeWallFunctionFvPatchScalarField
                 (
                     *this,
                     iF
@@ -195,6 +189,9 @@ public:
 
         // Evaluation functions
 
+            //- Evaluate the turbulent thermal diffusivity
+            tmp<scalarField> calcAlphat(const scalarField& prevAlphat) const;
+
             //- Update the coefficients associated with the patch field
             virtual void updateCoeffs();
 

applications/solvers/multiphase/reactingEulerFoam/derivedFvPatchFields/alphatWallBoilingWallFunction/alphatWallBoilingWallFunctionFvPatchScalarField.C

@@ -73,10 +73,8 @@ alphatWallBoilingWallFunctionFvPatchScalarField
     const DimensionedField<scalar, volMesh>& iF
 )
 :
-    alphatPhaseChangeJayatillekeWallFunctionFvPatchScalarField(p, iF),
+    alphatPhaseChangeWallFunctionFvPatchScalarField(p, iF),
-    otherPhaseName_("vapor"),
     phaseType_(liquidPhase),
-    relax_(0.1),
     AbyV_(p.size(), 0),
     alphatConv_(p.size(), 0),
     dDep_(p.size(), 1e-5),
@@ -103,10 +101,8 @@ alphatWallBoilingWallFunctionFvPatchScalarField
     const dictionary& dict
 )
 :
-    alphatPhaseChangeJayatillekeWallFunctionFvPatchScalarField(p, iF, dict),
+    alphatPhaseChangeWallFunctionFvPatchScalarField(p, iF, dict),
-    otherPhaseName_(dict.lookup("otherPhase")),
     phaseType_(phaseTypeNames_.read(dict.lookup("phaseType"))),
-    relax_(dict.lookupOrDefault<scalar>("relax", 0.1)),
     AbyV_(p.size(), 0),
     alphatConv_(p.size(), 0),
     dDep_(p.size(), 1e-5),
@@ -116,18 +112,6 @@ alphatWallBoilingWallFunctionFvPatchScalarField
     departureDiamModel_(nullptr),
     departureFreqModel_(nullptr)
 {
-
-    // Check that otherPhaseName != this phase
-    if (internalField().group() == otherPhaseName_)
-    {
-        FatalErrorInFunction
-            << "otherPhase should be the name of the vapor phase that "
-            << "corresponds to the liquid base of vice versa" << nl
-            << "This phase: " << internalField().group() << nl
-            << "otherPhase: " << otherPhaseName_
-            << abort(FatalError);
-    }
-
     switch (phaseType_)
     {
         case vaporPhase:
@@ -138,7 +122,7 @@ alphatWallBoilingWallFunctionFvPatchScalarField
                     dict.subDict("partitioningModel")
                 );
 
-            dmdt_ = 0;
+            dmdtf_ = 0;
 
             break;
         }
@@ -205,16 +189,14 @@ alphatWallBoilingWallFunctionFvPatchScalarField
     const fvPatchFieldMapper& mapper
 )
 :
-    alphatPhaseChangeJayatillekeWallFunctionFvPatchScalarField
+    alphatPhaseChangeWallFunctionFvPatchScalarField
     (
         psf,
         p,
         iF,
         mapper
     ),
-    otherPhaseName_(psf.otherPhaseName_),
     phaseType_(psf.phaseType_),
-    relax_(psf.relax_),
     AbyV_(psf.AbyV_),
     alphatConv_(mapper(psf.alphatConv_)),
     dDep_(mapper(psf.dDep_)),
@@ -232,10 +214,8 @@ alphatWallBoilingWallFunctionFvPatchScalarField
     const alphatWallBoilingWallFunctionFvPatchScalarField& psf
 )
 :
-    alphatPhaseChangeJayatillekeWallFunctionFvPatchScalarField(psf),
+    alphatPhaseChangeWallFunctionFvPatchScalarField(psf),
-    otherPhaseName_(psf.otherPhaseName_),
     phaseType_(psf.phaseType_),
-    relax_(psf.relax_),
     AbyV_(psf.AbyV_),
     alphatConv_(psf.alphatConv_),
     dDep_(psf.dDep_),
@@ -254,10 +234,8 @@ alphatWallBoilingWallFunctionFvPatchScalarField
     const DimensionedField<scalar, volMesh>& iF
 )
 :
-    alphatPhaseChangeJayatillekeWallFunctionFvPatchScalarField(psf, iF),
+    alphatPhaseChangeWallFunctionFvPatchScalarField(psf, iF),
-    otherPhaseName_(psf.otherPhaseName_),
     phaseType_(psf.phaseType_),
-    relax_(psf.relax_),
     AbyV_(psf.AbyV_),
     alphatConv_(psf.alphatConv_),
     dDep_(psf.dDep_),
@@ -271,53 +249,6 @@ alphatWallBoilingWallFunctionFvPatchScalarField
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
-bool alphatWallBoilingWallFunctionFvPatchScalarField::
-activePhasePair(const phasePairKey& phasePair) const
-{
-    if (phasePair == phasePairKey(otherPhaseName_, internalField().group()))
-    {
-        return true;
-    }
-    else
-    {
-        return false;
-    }
-}
-
-const scalarField& alphatWallBoilingWallFunctionFvPatchScalarField::
-dmdt(const phasePairKey& phasePair) const
-{
-    if (activePhasePair(phasePair))
-    {
-        return dmdt_;
-    }
-    else
-    {
-        FatalErrorInFunction
-            << " dmdt requested for invalid phasePair!"
-            << abort(FatalError);
-
-        return dmdt_;
-    }
-}
-
-const scalarField& alphatWallBoilingWallFunctionFvPatchScalarField::
-mDotL(const phasePairKey& phasePair) const
-{
-    if (activePhasePair(phasePair))
-    {
-        return mDotL_;
-    }
-    else
-    {
-        FatalErrorInFunction
-            << " mDotL requested for invalid phasePair!"
-            << abort(FatalError);
-
-        return mDotL_;
-    }
-}
-
 void alphatWallBoilingWallFunctionFvPatchScalarField::updateCoeffs()
 {
     if (updated())
@@ -594,13 +525,13 @@ void alphatWallBoilingWallFunctionFvPatchScalarField::updateCoeffs()
 
                 // Volumetric mass source in the near wall cell due to the
                 // wall boiling
-                dmdt_ =
+                dmdtf_ =
-                    (1 - relax_)*dmdt_
+                    (1 - relax_)*dmdtf_
                   + relax_*(1.0/6.0)*A2E*dDep_*rhoVaporw*fDep*AbyV_;
 
                 // Volumetric source in the near wall cell due to the wall
                 // boiling
-                mDotL_ = dmdt_*L;
+                dmdtLf_ = dmdtf_*L;
 
                 // Quenching heat transfer coefficient
                 const scalarField hQ
@@ -657,8 +588,8 @@ void alphatWallBoilingWallFunctionFvPatchScalarField::updateCoeffs()
                     Info<< "  A2: " << gMin(A2) << " - " << gMax(A2) << endl;
                     Info<< "  A2E: " << gMin(A2E) << " - "
                         << gMax(A2E) << endl;
-                    Info<< "  dmdtW: " << gMin(dmdt_) << " - "
+                    Info<< "  dmdtW: " << gMin(dmdtf_) << " - "
-                        << gMax(dmdt_) << endl;
+                        << gMax(dmdtf_) << endl;
                     Info<< "  qc: " << gMin(qc) << " - " << gMax(qc) << endl;
                     Info<< "  qq: " << gMin(fLiquid*qq()) << " - "
                         << gMax(fLiquid*qq()) << endl;
@@ -706,11 +637,12 @@ void alphatWallBoilingWallFunctionFvPatchScalarField::updateCoeffs()
 
 void alphatWallBoilingWallFunctionFvPatchScalarField::write(Ostream& os) const
 {
-    fvPatchField<scalar>::write(os);
+    alphatPhaseChangeWallFunctionFvPatchScalarField::write(os);
 
     writeEntry(os, "phaseType", phaseTypeNames_[phaseType_]);
-
+    writeEntry(os, "alphatConv", alphatConv_);
-    writeEntry(os, "relax", relax_);
+    writeEntry(os, "dDep", dDep_);
+    writeEntry(os, "qQuenching", qq_);
 
     switch (phaseType_)
     {
@@ -747,13 +679,6 @@ void alphatWallBoilingWallFunctionFvPatchScalarField::write(Ostream& os) const
             break;
         }
     }
-
-    writeEntry(os, "otherPhase", otherPhaseName_);
-    writeEntry(os, "dmdt", dmdt_);
-    writeEntry(os, "dDep", dDep_);
-    writeEntry(os, "qQuenching", qq_);
-    writeEntry(os, "alphatConv", alphatConv_);
-    writeEntry(os, "value", *this);
 }
 
 

applications/solvers/multiphase/reactingEulerFoam/derivedFvPatchFields/alphatWallBoilingWallFunction/alphatWallBoilingWallFunctionFvPatchScalarField.H

@@ -131,7 +131,7 @@ SourceFiles
 #ifndef compressible_alphatWallBoilingWallFunctionFvPatchScalarField_H
 #define compressible_alphatWallBoilingWallFunctionFvPatchScalarField_H
 
-#include "alphatPhaseChangeJayatillekeWallFunctionFvPatchScalarField.H"
+#include "alphatPhaseChangeWallFunctionFvPatchScalarField.H"
 #include "partitioningModel.H"
 #include "nucleationSiteModel.H"
 #include "departureDiameterModel.H"
@@ -150,7 +150,7 @@ namespace compressible
 
 class alphatWallBoilingWallFunctionFvPatchScalarField
 :
-    public alphatPhaseChangeJayatillekeWallFunctionFvPatchScalarField
+    public alphatPhaseChangeWallFunctionFvPatchScalarField
 {
 public:
 
@@ -163,23 +163,17 @@ public:
             liquidPhase
         };
 
+        //- Heat source type names
+        static const NamedEnum<phaseType, 2> phaseTypeNames_;
+
 
 private:
 
     // Private Data
 
-        //- name of the other phase (vapor/liquid phase)
-        word otherPhaseName_;
-
-        //- Heat source type names
-        static const NamedEnum<phaseType, 2> phaseTypeNames_;
-
         //- Heat source type
         phaseType phaseType_;
 
-        //- dmdt relaxationFactor
-        scalar relax_;
-
         //- Patch face area by cell volume
         scalarField AbyV_;
 
@@ -280,17 +274,6 @@ public:
 
     // Member Functions
 
-        using alphatPhaseChangeWallFunctionFvPatchScalarField::dmdt;
-
-        //- Is there phase change mass transfer for this phasePair
-        virtual bool activePhasePair(const phasePairKey&) const;
-
-        //- Return the rate of phase-change for specific phase pair
-        virtual const scalarField& dmdt(const phasePairKey&) const;
-
-        //- Return the rate of phase-change for specific phase pair
-        virtual const scalarField& mDotL(const phasePairKey&) const;
-
         //- Return the departure diameter field
         const scalarField& dDeparture() const
         {
@@ -306,9 +289,10 @@ public:
         //- Return the evaporation surface heat flux [W/m^2]
         tmp<scalarField> qe() const
         {
-            return mDotL_/AbyV_;
+            return dmdtLf_/AbyV_;
         }
 
+
         // Evaluation functions
 
             //- Update the coefficients associated with the patch field

applications/solvers/multiphase/reactingEulerFoam/interfacialCompositionModels/interfaceCompositionModels/interfaceCompositionModel/interfaceCompositionModel.C

@@ -135,20 +135,20 @@ Foam::tmp<Foam::volScalarField> Foam::interfaceCompositionModel::L
 }
 
 
-void Foam::interfaceCompositionModel::addMDotL
+void Foam::interfaceCompositionModel::addDmdtL
 (
     const volScalarField& K,
     const volScalarField& Tf,
-    volScalarField& mDotL,
+    volScalarField& dmdtL,
-    volScalarField& mDotLPrime
+    volScalarField& dmdtLPrime
 ) const
 {
     forAllConstIter(hashedWordList, species_, iter)
     {
         const volScalarField rhoKDL(thermo_.rho()*K*D(*iter)*L(*iter, Tf));
 
-        mDotL += rhoKDL*dY(*iter, Tf);
+        dmdtL += rhoKDL*dY(*iter, Tf);
-        mDotLPrime += rhoKDL*YfPrime(*iter, Tf);
+        dmdtLPrime += rhoKDL*YfPrime(*iter, Tf);
     }
 }
 

applications/solvers/multiphase/reactingEulerFoam/interfacialCompositionModels/interfaceCompositionModels/interfaceCompositionModel/interfaceCompositionModel.H

@@ -182,12 +182,12 @@ public:
             ) const;
 
             //- Add latent heat flow rate to total
-            virtual void addMDotL
+            virtual void addDmdtL
             (
                 const volScalarField& K,
                 const volScalarField& Tf,
-                volScalarField& mDotL,
+                volScalarField& dmdtL,
-                volScalarField& mDotLPrime
+                volScalarField& dmdtLPrime
             ) const;
 
 

applications/solvers/multiphase/reactingEulerFoam/interfacialModels/phaseTransferModels/deposition/deposition.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
@@ -63,8 +63,14 @@ Foam::phaseTransferModels::deposition::~deposition()
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
+bool Foam::phaseTransferModels::deposition::mixture() const
+{
+    return true;
+}
+
+
 Foam::tmp<Foam::volScalarField>
-Foam::phaseTransferModels::deposition::dmdt() const
+Foam::phaseTransferModels::deposition::dmdtf() const
 {
     const phaseModel* dropletPtr = nullptr;
     scalar sign = 1;

applications/solvers/multiphase/reactingEulerFoam/interfacialModels/phaseTransferModels/deposition/deposition.H

@@ -92,8 +92,11 @@ public:
 
     // Member Functions
 
-        //- The mass transfer rate
+        //- Whether or not there is a mass transfer rate for the mixture
-        virtual tmp<volScalarField> dmdt() const;
+        virtual bool mixture() const;
+
+        //- The mass transfer rate for the mixture
+        virtual tmp<volScalarField> dmdtf() const;
 };
 
 

applications/solvers/multiphase/reactingEulerFoam/interfacialModels/phaseTransferModels/phaseTransferModel/phaseTransferModel.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
@@ -39,6 +39,8 @@ namespace Foam
 const Foam::dimensionSet Foam::phaseTransferModel::dimDmdt =
     Foam::dimDensity/Foam::dimTime;
 
+const Foam::hashedWordList Foam::phaseTransferModel::noSpecies_ =
+    Foam::hashedWordList();
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
@@ -48,6 +50,15 @@ Foam::phaseTransferModel::phaseTransferModel
     const phasePair& pair
 )
 :
+    regIOobject
+    (
+        IOobject
+        (
+            IOobject::groupName(typeName, pair.name()),
+            pair.phase1().mesh().time().timeName(),
+            pair.phase1().mesh()
+        )
+    ),
     pair_(pair)
 {}
 
@@ -58,4 +69,37 @@ Foam::phaseTransferModel::~phaseTransferModel()
 {}
 
 
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+bool Foam::phaseTransferModel::mixture() const
+{
+    return false;
+}
+
+
+Foam::tmp<Foam::volScalarField> Foam::phaseTransferModel::dmdtf() const
+{
+    return tmp<volScalarField>(nullptr);
+}
+
+
+const Foam::hashedWordList& Foam::phaseTransferModel::species() const
+{
+    return noSpecies_;
+}
+
+
+Foam::HashPtrTable<Foam::volScalarField>
+Foam::phaseTransferModel::dmidtf() const
+{
+    return HashPtrTable<volScalarField>();
+}
+
+
+bool Foam::phaseTransferModel::writeData(Ostream& os) const
+{
+    return os.good();
+}
+
+
 // ************************************************************************* //

applications/solvers/multiphase/reactingEulerFoam/interfacialModels/phaseTransferModels/phaseTransferModel/phaseTransferModel.H

@@ -38,6 +38,8 @@ SourceFiles
 
 #include "volFields.H"
 #include "dictionary.H"
+#include "HashPtrTable.H"
+#include "hashedWordList.H"
 #include "runTimeSelectionTables.H"
 
 namespace Foam
@@ -50,6 +52,8 @@ class phasePair;
 \*---------------------------------------------------------------------------*/
 
 class phaseTransferModel
+:
+    public regIOobject
 {
 protected:
 
@@ -85,6 +89,9 @@ public:
         //- Mass transfer rate dimensions
         static const dimensionSet dimDmdt;
 
+        //- Empty species list
+        static const hashedWordList noSpecies_;
+
 
     // Constructors
 
@@ -111,29 +118,20 @@ public:
 
     // Member Functions
 
-        //- Return the active species for the current pair,
+        //- Whether or not there is a mass transfer rate for the mixture
-        //  by default returns empty list
+        virtual bool mixture() const;
-        virtual const wordList activeSpecies() const
-        {
-            return wordList();
-        };
 
-        //- Do possible corrections before dmdt calculations
+        //- The mass transfer rate for the mixture
-        virtual void correct() {};
+        virtual tmp<volScalarField> dmdtf() const;
 
-        //- The mass transfer rate
+        //- The list of individual species that are transferred
-        virtual tmp<volScalarField> dmdt() const = 0;
+        virtual const hashedWordList& species() const;
 
-        //- The mass transfer rate for the given species, this needs to be
+        //- The mass transfer rate for individual species
-        //  implemented only if the model uses non uniform species transfer
+        virtual HashPtrTable<volScalarField> dmidtf() const;
-        virtual tmp<volScalarField> speciesDmdt
+
-        (
+        //- Dummy write for regIOobject
-            const word speciesName
+        bool writeData(Ostream& os) const;
-        ) const
-        {
-            NotImplemented;
-            return dmdt();
-        };
 };
 
 

applications/solvers/multiphase/reactingEulerFoam/interfacialModels/phaseTransferModels/reactionDriven/reactionDriven.C

@@ -74,44 +74,33 @@ Foam::phaseTransferModels::reactionDriven::~reactionDriven()
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
-Foam::tmp<Foam::volScalarField>
+const Foam::hashedWordList&
-Foam::phaseTransferModels::reactionDriven::dmdt() const
+Foam::phaseTransferModels::reactionDriven::species() const
 {
-    volScalarField dmdt
+    return species_;
-    (
-        volScalarField::New
-        (
-            "reactionDrivendmdt",
-            pair_.phase1().mesh(),
-            dimensionedScalar(Foam::phaseTransferModel::dimDmdt, 0)
-        )
-    );
-
-    forAllConstIter
-    (
-        wordList,
-        species_,
-        sIter
-    )
-    {
-        dmdt += speciesDmdt(*sIter);
-    }
-
-    return dmdt;
 }
 
 
-Foam::tmp<Foam::volScalarField>
+Foam::HashPtrTable<Foam::volScalarField>
-Foam::phaseTransferModels::reactionDriven::speciesDmdt
+Foam::phaseTransferModels::reactionDriven::dmidtf() const
-(
-    const word speciesName
-) const
 {
-    // Reaction rate query requires non const Y
+    HashPtrTable<volScalarField> result;
-    volScalarField& Y =
-        const_cast<volScalarField&>(reactingPhase_.Y(speciesName));
 
-    return sign_*reactingPhase_*reactingPhase_.R(Y) & Y;
+    forAll(species_, i)
+    {
+        const word name = species_[i];
+
+        volScalarField& Y =
+            const_cast<volScalarField&>(reactingPhase_.Y(name));
+
+        result.set
+        (
+            species_[i],
+            (sign_*reactingPhase_*reactingPhase_.R(Y) & Y).ptr()
+        );
+    }
+
+    return result;
 };
 
 

applications/solvers/multiphase/reactingEulerFoam/interfacialModels/phaseTransferModels/reactionDriven/reactionDriven.H

@@ -74,7 +74,7 @@ private:
         const scalar sign_;
 
         //- List of species changing phase
-        wordList species_;
+        const hashedWordList species_;
 
 
 public:
@@ -99,20 +99,11 @@ public:
 
     // Member Functions
 
-        //- Return the active species for the current pair
+        //- The list of individual species that are transferred
-        virtual const wordList activeSpecies() const
+        virtual const hashedWordList& species() const;
-        {
-            return species_;
-        };
 
-        //- The mass transfer rate
+        //- The mass transfer rate for individual species
-        virtual tmp<volScalarField> dmdt() const;
+        virtual HashPtrTable<volScalarField> dmidtf() const;
-
-        //- The mass transfer rate for the given species
-        virtual tmp<volScalarField> speciesDmdt
-        (
-            const word speciesName
-        ) const;
 };
 
 

applications/solvers/multiphase/reactingEulerFoam/phaseSystems/BlendedInterfacialModel/BlendedInterfacialModel.C

@@ -31,46 +31,97 @@ License
 
 namespace Foam
 {
+namespace blendedInterfacialModel
+{
 
-// * * * * * * * * * * * * * Static Member Functions * * * * * * * * * * * * //
+template<class GeoField>
+inline tmp<GeoField> interpolate(tmp<volScalarField> f);
 
 template<>
-inline tmp<Foam::volScalarField>
+inline tmp<Foam::volScalarField> interpolate(tmp<volScalarField> f)
-blendedInterfacialModel::interpolate(tmp<volScalarField> f)
 {
     return f;
 }
 
-
 template<>
-inline tmp<Foam::surfaceScalarField>
+inline tmp<Foam::surfaceScalarField> interpolate(tmp<volScalarField> f)
-blendedInterfacialModel::interpolate(tmp<volScalarField> f)
 {
     return fvc::interpolate(f);
 }
 
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+} // End namespace blendedInterfacialModel
-
 } // End namespace Foam
 
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
 template<class ModelType>
-template<class GeoField>
+template<template<class> class PatchField, class GeoMesh>
-void Foam::BlendedInterfacialModel<ModelType>::correctFixedFluxBCs
+void Foam::BlendedInterfacialModel<ModelType>::calculateBlendingCoeffs
 (
-    GeoField& field
+    tmp<GeometricField<scalar, PatchField, GeoMesh>>& f1,
+    tmp<GeometricField<scalar, PatchField, GeoMesh>>& f2,
+    const bool subtract
 ) const
 {
-    typename GeoField::Boundary& fieldBf = field.boundaryFieldRef();
+    typedef GeometricField<scalar, PatchField, GeoMesh> scalarGeoField;
 
-    forAll(phase1_.phi()().boundaryField(), patchi)
+    if (model_.valid() && subtract)
+    {
+        FatalErrorInFunction
+            << "Cannot treat an interfacial model with no distinction between "
+            << "continuous and dispersed phases as signed"
+            << exit(FatalError);
+    }
+
+    if (model_.valid() || model1In2_.valid())
+    {
+        f1 =
+            blendedInterfacialModel::interpolate<scalarGeoField>
+            (
+                blending_.f1(phase1_, phase2_)
+            );
+    }
+
+    if (model_.valid() || model2In1_.valid())
+    {
+        f2 =
+            (subtract ? -1 : +1)
+           *blendedInterfacialModel::interpolate<scalarGeoField>
+            (
+                blending_.f2(phase1_, phase2_)
+            );
+    }
+}
+
+
+template<class ModelType>
+template<class Type, template<class> class PatchField, class GeoMesh>
+void Foam::BlendedInterfacialModel<ModelType>::correctFixedFluxBCs
+(
+    GeometricField<Type, PatchField, GeoMesh>& field
+) const
+{
+    typedef GeometricField<Type, PatchField, GeoMesh> typeGeoField;
+
+    typename typeGeoField::Boundary& fieldBf = field.boundaryFieldRef();
+
+    forAll(fieldBf, patchi)
     {
         if
         (
-            isA<fixedValueFvsPatchScalarField>
             (
-                phase1_.phi()().boundaryField()[patchi]
+                !phase1_.stationary()
+             && isA<fixedValueFvsPatchScalarField>
+                (
+                    phase1_.phi()().boundaryField()[patchi]
+                )
+            )
+         || (
+                !phase2_.stationary()
+             && isA<fixedValueFvsPatchScalarField>
+                (
+                    phase2_.phi()().boundaryField()[patchi]
+                )
             )
         )
         {
@@ -103,73 +154,30 @@ Foam::BlendedInterfacialModel<ModelType>::evaluate
     typedef GeometricField<Type, PatchField, GeoMesh> typeGeoField;
 
     tmp<scalarGeoField> f1, f2;
+    calculateBlendingCoeffs(f1, f2, subtract);
 
-    if (model_.valid() || model1In2_.valid())
+    tmp<typeGeoField> x =
-    {
+        typeGeoField::New
-        f1 =
-            blendedInterfacialModel::interpolate<scalarGeoField>
-            (
-                blending_.f1(phase1_, phase2_)
-            );
-    }
-
-    if (model_.valid() || model2In1_.valid())
-    {
-        f2 =
-            blendedInterfacialModel::interpolate<scalarGeoField>
-            (
-                blending_.f2(phase1_, phase2_)
-            );
-    }
-
-    tmp<typeGeoField> x
-    (
-        new typeGeoField
         (
-            IOobject
+            ModelType::typeName + ":"
-            (
+          + IOobject::groupName(name, phasePair(phase1_, phase2_).name()),
-                ModelType::typeName + ":" + name,
-                phase1_.mesh().time().timeName(),
-                phase1_.mesh(),
-                IOobject::NO_READ,
-                IOobject::NO_WRITE,
-                false
-            ),
             phase1_.mesh(),
-            dimensioned<Type>("zero", dims, Zero)
+            dimensioned<Type>(dims, Zero)
-        )
+        );
-    );
 
     if (model_.valid())
     {
-        if (subtract)
+        x.ref() += (scalar(1) - f1() - f2())*(model_().*method)(args ...);
-        {
-            FatalErrorInFunction
-                << "Cannot treat an interfacial model with no distinction "
-                << "between continuous and dispersed phases as signed"
-                << exit(FatalError);
-        }
-
-        x.ref() += (model_().*method)(args ...)*(scalar(1) - f1() - f2());
     }
 
     if (model1In2_.valid())
     {
-        x.ref() += (model1In2_().*method)(args ...)*f1;
+        x.ref() += f1*(model1In2_().*method)(args ...);
     }
 
     if (model2In1_.valid())
     {
-        tmp<typeGeoField> dx = (model2In1_().*method)(args ...)*f2;
+        x.ref() += f2*(model2In1_().*method)(args ...);
-
-        if (subtract)
-        {
-            x.ref() -= dx;
-        }
-        else
-        {
-            x.ref() += dx;
-        }
     }
 
     if
@@ -185,6 +193,96 @@ Foam::BlendedInterfacialModel<ModelType>::evaluate
 }
 
 
+template<class ModelType>
+template
+<
+    class Type,
+    template<class> class PatchField,
+    class GeoMesh,
+    class ... Args
+>
+Foam::HashPtrTable<Foam::GeometricField<Type, PatchField, GeoMesh>>
+Foam::BlendedInterfacialModel<ModelType>::evaluate
+(
+    HashPtrTable<GeometricField<Type, PatchField, GeoMesh>>
+    (ModelType::*method)(Args ...) const,
+    const word& name,
+    const dimensionSet& dims,
+    const bool subtract,
+    Args ... args
+) const
+{
+    typedef GeometricField<scalar, PatchField, GeoMesh> scalarGeoField;
+    typedef GeometricField<Type, PatchField, GeoMesh> typeGeoField;
+
+    tmp<scalarGeoField> f1, f2;
+    calculateBlendingCoeffs(f1, f2, subtract);
+
+    HashPtrTable<typeGeoField> xs;
+
+    auto addToXs = [&]
+    (
+        const scalarGeoField& f,
+        const HashPtrTable<typeGeoField>& dxs
+    )
+    {
+        forAllConstIter(typename HashPtrTable<typeGeoField>, dxs, dxIter)
+        {
+            if (xs.found(dxIter.key()))
+            {
+                *xs[dxIter.key()] += f**dxIter();
+            }
+            else
+            {
+                xs.insert
+                (
+                    dxIter.key(),
+                    typeGeoField::New
+                    (
+                        ModelType::typeName + ':'
+                      + IOobject::groupName
+                        (
+                            IOobject::groupName(name, dxIter.key()),
+                            phasePair(phase1_, phase2_).name()
+                        ),
+                        f**dxIter()
+                    ).ptr()
+                );
+            }
+        }
+    };
+
+    if (model_.valid())
+    {
+        addToXs(scalar(1) - f1() - f2(), (model_().*method)(args ...));
+    }
+
+    if (model1In2_.valid())
+    {
+        addToXs(f1, (model1In2_().*method)(args ...));
+    }
+
+    if (model2In1_.valid())
+    {
+        addToXs(f2, (model1In2_().*method)(args ...));
+    }
+
+    if
+    (
+        correctFixedFluxBCs_
+     && (model_.valid() || model1In2_.valid() || model2In1_.valid())
+    )
+    {
+        forAllIter(typename HashPtrTable<typeGeoField>, xs, xIter)
+        {
+            correctFixedFluxBCs(*xIter());
+        }
+    }
+
+    return xs;
+}
+
+
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
 template<class ModelType>
@@ -290,46 +388,6 @@ Foam::BlendedInterfacialModel<ModelType>::~BlendedInterfacialModel()
 
 // * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
 
-template<class ModelType>
-void Foam::BlendedInterfacialModel<ModelType>::correct()
-{
-    if (model1In2_.valid())
-    {
-        model1In2_->correct();
-    }
-    if (model2In1_.valid())
-    {
-        model2In1_->correct();
-    }
-    if (model_.valid())
-    {
-        model_->correct();
-    }
-}
-
-
-template<class ModelType>
-const Foam::wordList Foam::BlendedInterfacialModel<ModelType>::
-activeSpecies() const
-{
-    wordList activeSpeciesList;
-    if (model1In2_.valid())
-    {
-        activeSpeciesList.append(model1In2_->activeSpecies());
-    }
-    if (model2In1_.valid())
-    {
-        activeSpeciesList.append(model2In1_->activeSpecies());
-    }
-    if (model_.valid())
-    {
-        activeSpeciesList.append(model_->activeSpecies());
-    }
-
-    return activeSpeciesList;
-}
-
-
 template<class ModelType>
 Foam::tmp<Foam::volScalarField>
 Foam::BlendedInterfacialModel<ModelType>::K() const
@@ -384,29 +442,50 @@ Foam::BlendedInterfacialModel<ModelType>::D() const
 
 
 template<class ModelType>
-Foam::tmp<Foam::volScalarField>
+bool Foam::BlendedInterfacialModel<ModelType>::mixture() const
-Foam::BlendedInterfacialModel<ModelType>::dmdt() const
 {
-    return evaluate(&ModelType::dmdt, "dmdt", ModelType::dimDmdt, false);
+    return
+        (model1In2_.valid() && model1In2_->mixture())
+     || (model2In1_.valid() && model2In1_->mixture())
+     || (model_.valid() && model_->mixture());
 }
 
 
 template<class ModelType>
 Foam::tmp<Foam::volScalarField>
-Foam::BlendedInterfacialModel<ModelType>::speciesDmdt
+Foam::BlendedInterfacialModel<ModelType>::dmdtf() const
-(
-    const word speciesName
-) const
 {
-    return
+    return evaluate(&ModelType::dmdtf, "dmdtf", ModelType::dimDmdt, true);
-        evaluate
+}
-        (
+
-            &ModelType::speciesDmdt,
+
-            "speciesDmdt",
+template<class ModelType>
-            ModelType::dimDmdt,
+Foam::hashedWordList Foam::BlendedInterfacialModel<ModelType>::species() const
-            false,
+{
-            speciesName
+    wordList species;
-        );
+
+    if (model1In2_.valid())
+    {
+        species.append(model1In2_->species());
+    }
+    if (model2In1_.valid())
+    {
+        species.append(model2In1_->species());
+    }
+    if (model_.valid())
+    {
+        species.append(model_->species());
+    }
+
+    return hashedWordList(move(species));
+}
+
+
+template<class ModelType>
+Foam::HashPtrTable<Foam::volScalarField>
+Foam::BlendedInterfacialModel<ModelType>::dmidtf() const
+{
+    return evaluate(&ModelType::dmidtf, "dmidtf", ModelType::dimDmdt, true);
 }
 
 

applications/solvers/multiphase/reactingEulerFoam/phaseSystems/BlendedInterfacialModel/BlendedInterfacialModel.H

@@ -37,6 +37,8 @@ SourceFiles
 #include "blendingMethod.H"
 #include "phasePair.H"
 #include "orderedPhasePair.H"
+#include "HashPtrTable.H"
+#include "hashedWordList.H"
 #include "geometricZeroField.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@@ -44,21 +46,6 @@ SourceFiles
 namespace Foam
 {
 
-/*---------------------------------------------------------------------------*\
-                    Class blendedInterfacialModel Declaration
-\*---------------------------------------------------------------------------*/
-
-class blendedInterfacialModel
-{
-    public:
-
-        //- Convenience function to interpolate blending values. Needs to be
-        //  specialised, so can't sit in the templated class.
-        template<class GeoField>
-        static inline tmp<GeoField> interpolate(tmp<volScalarField> f);
-};
-
-
 /*---------------------------------------------------------------------------*\
                     Class BlendedInterfacialModel Declaration
 \*---------------------------------------------------------------------------*/
@@ -94,9 +81,21 @@ class BlendedInterfacialModel
 
     // Private Member Functions
 
+        //- Calculate the blending coefficients
+        template<template<class> class PatchField, class GeoMesh>
+        void calculateBlendingCoeffs
+        (
+            tmp<GeometricField<scalar, PatchField, GeoMesh>>& f1,
+            tmp<GeometricField<scalar, PatchField, GeoMesh>>& f2,
+            const bool subtract
+        ) const;
+
         //- Correct coeff/value on fixed flux boundary conditions
-        template<class GeoField>
+        template<class Type, template<class> class PatchField, class GeoMesh>
-        void correctFixedFluxBCs(GeoField& field) const;
+        void correctFixedFluxBCs
+        (
+            GeometricField<Type, PatchField, GeoMesh>& field
+        ) const;
 
         //- Return the blended coeff/value
         template
@@ -116,6 +115,24 @@ class BlendedInterfacialModel
             Args ... args
         ) const;
 
+        //- Return the blended coeff/value
+        template
+        <
+            class Type,
+            template<class> class PatchField,
+            class GeoMesh,
+            class ... Args
+        >
+        HashPtrTable<GeometricField<Type, PatchField, GeoMesh>> evaluate
+        (
+            HashPtrTable<GeometricField<Type, PatchField, GeoMesh>>
+            (ModelType::*method)(Args ...) const,
+            const word& name,
+            const dimensionSet& dims,
+            const bool subtract,
+            Args ... args
+        ) const;
+
 
 public:
 
@@ -162,12 +179,6 @@ public:
 
     // Member Functions
 
-        //- Correct the underlying models
-        void correct();
-
-        //- Return the list of active species
-        const wordList activeSpecies() const;
-
         //- Return the blended force coefficient
         tmp<volScalarField> K() const;
 
@@ -187,11 +198,17 @@ public:
         //- Return the blended diffusivity
         tmp<volScalarField> D() const;
 
-        //- Return the blended mass transfer rate
+        //- Return the list of individual species that are transferred
-        tmp<volScalarField> dmdt() const;
+        bool mixture() const;
 
-        //- Return the blended species transfer rate
+        //- Return the blended mass transfer rate
-        tmp<volScalarField> speciesDmdt(const word speciesName) const;
+        tmp<volScalarField> dmdtf() const;
+
+        //- Return the list of individual species that are transferred
+        hashedWordList species() const;
+
+        //- Return the blended mass transfer rates for individual species
+        HashPtrTable<volScalarField> dmidtf() const;
 
         //- Dummy write for regIOobject
         bool writeData(Ostream& os) const;

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

@@ -32,13 +32,12 @@ License
 // * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * * //
 
 template<class BasePhaseSystem>
-Foam::autoPtr<Foam::phaseSystem::dmdtTable>
+Foam::autoPtr<Foam::phaseSystem::dmdtfTable>
 Foam::InterfaceCompositionPhaseChangePhaseSystem<BasePhaseSystem>::
-iDmdts() const
+dmdtfs() const
 {
-    autoPtr<phaseSystem::dmdtTable> iDmdtsPtr(new phaseSystem::dmdtTable);
+    autoPtr<phaseSystem::dmdtfTable> dmdtPtr(new phaseSystem::dmdtfTable);
-
+    phaseSystem::dmdtfTable& dmdtfs = dmdtPtr();
-    phaseSystem::dmdtTable& iDmdts = iDmdtsPtr();
 
     forAllConstIter
     (
@@ -66,43 +65,43 @@ iDmdts() const
             (
                 hashedWordList,
                 compositionModel.species(),
-                memberIter
+                specieIter
             )
             {
-                const word& member = *memberIter;
+                const word& specie = *specieIter;
 
-                const label dmidtSign = pairIter.index() == 0 ? +1 : -1;
+                tmp<volScalarField> dmidtf
-
-                tmp<volScalarField> dmidt
                 (
-                    *(*iDmdtSu_[pair])[member]
+                    (pairIter.index() == 0 ? +1 : -1)
-                  + *(*iDmdtSp_[pair])[member]*phase.Y(member)
+                   *(
+                        *(*dmidtfSus_[pair])[specie]
+                      + *(*dmidtfSps_[pair])[specie]*phase.Y(specie)
+                    )
                 );
 
-                if (iDmdts.found(pair))
+                if (dmdtfs.found(pair))
                 {
-                    *iDmdts[pair] += dmidtSign*dmidt;
+                    *dmdtfs[pair] += dmidtf;
                 }
                 else
                 {
-                    iDmdts.insert(pair, (dmidtSign*dmidt).ptr());
+                    dmdtfs.insert(pair, dmidtf.ptr());
                 }
             }
         }
     }
 
-    return iDmdtsPtr;
+    return dmdtPtr;
 }
 
 
 template<class BasePhaseSystem>
-Foam::autoPtr<Foam::phaseSystem::dmidtTable>
+Foam::autoPtr<Foam::phaseSystem::dmidtfTable>
 Foam::InterfaceCompositionPhaseChangePhaseSystem<BasePhaseSystem>::
-iDmidts() const
+dmidtfs() const
 {
-    autoPtr<phaseSystem::dmidtTable> iDmidtsPtr(new phaseSystem::dmidtTable);
+    autoPtr<phaseSystem::dmidtfTable> dmidtfsPtr(new phaseSystem::dmidtfTable);
-
+    phaseSystem::dmidtfTable& dmidtfs = dmidtfsPtr();
-    phaseSystem::dmidtTable& iDmidts = iDmidtsPtr();
 
     forAllConstIter
     (
@@ -114,9 +113,9 @@ iDmidts() const
         const phasePair& pair =
             this->phasePairs_[interfaceCompositionModelIter.key()];
 
-        if (!iDmidts.found(pair))
+        if (!dmidtfs.found(pair))
         {
-            iDmidts.insert(pair, new HashPtrTable<volScalarField>());
+            dmidtfs.insert(pair, new HashPtrTable<volScalarField>());
         }
 
         forAllConstIter(phasePair, pair, pairIter)
@@ -135,32 +134,34 @@ iDmidts() const
             (
                 hashedWordList,
                 compositionModel.species(),
-                memberIter
+                specieIter
             )
             {
-                const word& member = *memberIter;
+                const word& specie = *specieIter;
 
-                const label dmidtSign = pairIter.index() == 0 ? +1 : -1;
 
-                tmp<volScalarField> dmidt
+                tmp<volScalarField> dmidtf
                 (
-                    *(*iDmdtSu_[pair])[member]
+                    (pairIter.index() == 0 ? +1 : -1)
-                  + *(*iDmdtSp_[pair])[member]*phase.Y(member)
+                   *(
+                        *(*dmidtfSus_[pair])[specie]
+                      + *(*dmidtfSps_[pair])[specie]*phase.Y(specie)
+                    )
                 );
 
-                if (iDmidts[pair]->found(member))
+                if (dmidtfs[pair]->found(specie))
                 {
-                    *(*iDmidts[pair])[member] += dmidtSign*dmidt;
+                    *(*dmidtfs[pair])[specie] += dmidtf;
                 }
                 else
                 {
-                    iDmidts[pair]->insert(member, (dmidtSign*dmidt).ptr());
+                    dmidtfs[pair]->insert(specie, dmidtf.ptr());
                 }
             }
         }
     }
 
-    return iDmidtsPtr;
+    return dmidtfsPtr;
 }
 
 
@@ -255,8 +256,8 @@ InterfaceCompositionPhaseChangePhaseSystem
         const phasePair& pair =
             this->phasePairs_[interfaceCompositionModelIter.key()];
 
-        iDmdtSu_.insert(pair, new HashPtrTable<volScalarField>());
+        dmidtfSus_.insert(pair, new HashPtrTable<volScalarField>());
-        iDmdtSp_.insert(pair, new HashPtrTable<volScalarField>());
+        dmidtfSps_.insert(pair, new HashPtrTable<volScalarField>());
 
         forAllConstIter(phasePair, pair, pairIter)
         {
@@ -272,23 +273,57 @@ InterfaceCompositionPhaseChangePhaseSystem
             (
                 hashedWordList,
                 compositionModel.species(),
-                memberIter
+                specieIter
             )
             {
-                const word& member = *memberIter;
+                const word& specie = *specieIter;
 
-                iDmdtSu_[pair]->insert
+                dmidtfSus_[pair]->insert
                 (
-                    member,
+                    specie,
-                    zeroVolField<scalar>(pair, "iDmdtSu", dimDensity/dimTime)
+                    new volScalarField
-                   .ptr()
+                    (
+                        IOobject
+                        (
+                            IOobject::groupName
+                            (
+                                IOobject::groupName
+                                (
+                                    "interfaceCompositionPhaseChange:dmidtfSu",
+                                    specie
+                                ),
+                                pair.name()
+                            ),
+                            this->mesh().time().timeName(),
+                            this->mesh()
+                        ),
+                        this->mesh(),
+                        dimensionedScalar(dimDensity/dimTime, 0)
+                    )
                 );
 
-                iDmdtSp_[pair]->insert
+                dmidtfSps_[pair]->insert
                 (
-                    member,
+                    specie,
-                    zeroVolField<scalar>(pair, "iDmdtSp", dimDensity/dimTime)
+                    new volScalarField
-                   .ptr()
+                    (
+                        IOobject
+                        (
+                            IOobject::groupName
+                            (
+                                IOobject::groupName
+                                (
+                                    "interfaceCompositionPhaseChange:dmidtfSp",
+                                    specie
+                                ),
+                                pair.name()
+                            ),
+                            this->mesh().time().timeName(),
+                            this->mesh()
+                        ),
+                        this->mesh(),
+                        dimensionedScalar(dimDensity/dimTime, 0)
+                    )
                 );
             }
         }
@@ -306,24 +341,71 @@ Foam::InterfaceCompositionPhaseChangePhaseSystem<BasePhaseSystem>::
 
 // * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
 
+template<class BasePhaseSystem>
+Foam::tmp<Foam::volScalarField>
+Foam::InterfaceCompositionPhaseChangePhaseSystem<BasePhaseSystem>::dmdtf
+(
+    const phasePairKey& key
+) const
+{
+    tmp<volScalarField> tDmdtf = BasePhaseSystem::dmdtf(key);
+
+    if (interfaceCompositionModels_.found(key))
+    {
+        const phasePair& pair = this->phasePairs_[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)
+                    );
+            }
+        }
+    }
+
+    return tDmdtf;
+}
+
+
 template<class BasePhaseSystem>
 Foam::PtrList<Foam::volScalarField>
 Foam::InterfaceCompositionPhaseChangePhaseSystem<BasePhaseSystem>::dmdts() const
 {
     PtrList<volScalarField> dmdts(BasePhaseSystem::dmdts());
 
-    autoPtr<phaseSystem::dmdtTable> iDmdtsPtr = this->iDmdts();
+    autoPtr<phaseSystem::dmdtfTable> dmdtfsPtr = this->dmdtfs();
+    const phaseSystem::dmdtfTable& dmdtfs = dmdtfsPtr();
 
-    const phaseSystem::dmdtTable& iDmdts = iDmdtsPtr();
+    forAllConstIter(phaseSystem::dmdtfTable, dmdtfs, dmdtfsIter)
-
-    forAllConstIter(phaseSystem::dmdtTable, iDmdts, iDmdtIter)
     {
-        const phasePair& pair = this->phasePairs_[iDmdtIter.key()];
+        const phasePair& pair = this->phasePairs_[dmdtfsIter.key()];
         const phaseModel& phase = pair.phase1();
         const phaseModel& otherPhase = pair.phase2();
 
-        addField(phase, "dmdt", *iDmdtIter(), dmdts);
+        addField(phase, "dmdt", *dmdtfsIter(), dmdts);
-        addField(otherPhase, "dmdt", - *iDmdtIter(), dmdts);
+        addField(otherPhase, "dmdt", - *dmdtfsIter(), dmdts);
     }
 
     return dmdts;
@@ -340,7 +422,7 @@ momentumTransfer()
 
     phaseSystem::momentumTransferTable& eqns = eqnsPtr();
 
-    this->addDmdtU(iDmdts(), eqns);
+    this->addDmdtUfs(dmdtfs(), eqns);
 
     return eqnsPtr;
 }
@@ -356,7 +438,7 @@ momentumTransferf()
 
     phaseSystem::momentumTransferTable& eqns = eqnsPtr();
 
-    this->addDmdtU(iDmdts(), eqns);
+    this->addDmdtUfs(dmdtfs(), eqns);
 
     return eqnsPtr;
 }
@@ -372,7 +454,7 @@ heatTransfer() const
 
     phaseSystem::heatTransferTable& eqns = eqnsPtr();
 
-    this->addDmidtHe(iDmidts(), eqns);
+    this->addDmidtHef(dmidtfs(), eqns);
 
     return eqnsPtr;
 }
@@ -390,7 +472,7 @@ specieTransfer() const
 
     // Explicit
     /*
-    this->addDmidtY(iDmidts(), eqns);
+    this->addDmidtYf(dmidtfs(), eqns);
     */
 
     // Semi-implicit
@@ -421,22 +503,22 @@ specieTransfer() const
             (
                 hashedWordList,
                 compositionModel.species(),
-                memberIter
+                specieIter
             )
             {
-                const word& member = *memberIter;
+                const word& specie = *specieIter;
 
                 // Implicit transport through this phase
-                *eqns[phase.Y(member).name()] +=
+                *eqns[phase.Y(specie).name()] +=
-                    *(*iDmdtSu_[pair])[member]
+                    *(*dmidtfSus_[pair])[specie]
-                  + fvm::Sp(*(*iDmdtSp_[pair])[member], phase.Y(member));
+                  + fvm::Sp(*(*dmidtfSps_[pair])[specie], phase.Y(specie));
 
                 // Explicit transport out of the other phase
-                if (eqns.found(IOobject::groupName(member, otherPhase.name())))
+                if (eqns.found(IOobject::groupName(specie, otherPhase.name())))
                 {
-                    *eqns[otherPhase.Y(member).name()] -=
+                    *eqns[otherPhase.Y(specie).name()] -=
-                        *(*iDmdtSu_[pair])[member]
+                        *(*dmidtfSus_[pair])[specie]
-                      + *(*iDmdtSp_[pair])[member]*phase.Y(member);
+                      + *(*dmidtfSps_[pair])[specie]*phase.Y(specie);
                 }
             }
         }
@@ -486,16 +568,16 @@ correct()
             (
                 hashedWordList,
                 compositionModel.species(),
-                memberIter
+                specieIter
             )
             {
-                const word& member = *memberIter;
+                const word& specie = *specieIter;
 
-                const volScalarField KD(K*compositionModel.D(member));
+                const volScalarField KD(K*compositionModel.D(specie));
-                const volScalarField Yf(compositionModel.Yf(member, Tf));
+                const volScalarField Yf(compositionModel.Yf(specie, Tf));
 
-                *(*iDmdtSu_[pair])[member] = phase.rho()*KD*Yf;
+                *(*dmidtfSus_[pair])[specie] = phase.rho()*KD*Yf;
-                *(*iDmdtSp_[pair])[member] = - phase.rho()*KD;
+                *(*dmidtfSps_[pair])[specie] = - phase.rho()*KD;
             }
         }
     }
@@ -546,40 +628,40 @@ correctInterfaceThermo()
 
         for (label i = 0; i < nInterfaceCorrectors_; ++ i)
         {
-            tmp<volScalarField> mDotL =
+            tmp<volScalarField> dmdtLf =
-                zeroVolField<scalar>
+                volScalarField::New
                 (
-                    pair,
+                    IOobject::groupName("dmdtLf", pair.name()),
-                    "mDotL",
+                    this->mesh(),
-                    dimEnergy/dimVolume/dimTime
+                    dimensionedScalar(dimEnergy/dimVolume/dimTime, 0)
                 );
-            tmp<volScalarField> mDotLPrime =
+            tmp<volScalarField> dmdtLfPrime =
-                zeroVolField<scalar>
+                volScalarField::New
                 (
-                    pair,
+                    IOobject::groupName("dmdtLfPrime", pair.name()),
-                    "mDotLPrime",
+                    this->mesh(),
-                    mDotL().dimensions()/dimTemperature
+                    dimensionedScalar(dmdtLf().dimensions()/dimTemperature, 0)
                 );
 
             // Add latent heats from forward and backward models
             if (this->interfaceCompositionModels_[pair].first().valid())
             {
-                this->interfaceCompositionModels_[pair].first()->addMDotL
+                this->interfaceCompositionModels_[pair].first()->addDmdtL
                 (
                     diffusiveMassTransferModels_[pair].first()->K(),
                     Tf,
-                    mDotL.ref(),
+                    dmdtLf.ref(),
-                    mDotLPrime.ref()
+                    dmdtLfPrime.ref()
                 );
             }
             if (this->interfaceCompositionModels_[pair].second().valid())
             {
-                this->interfaceCompositionModels_[pair].second()->addMDotL
+                this->interfaceCompositionModels_[pair].second()->addDmdtL
                 (
                   - diffusiveMassTransferModels_[pair].second()->K(),
                     Tf,
-                    mDotL.ref(),
+                    dmdtLf.ref(),
-                    mDotLPrime.ref()
+                    dmdtLfPrime.ref()
                 );
             }
 
@@ -589,10 +671,10 @@ correctInterfaceThermo()
                 (
                     H1*(Tf - pair.phase1().thermo().T())
                   + H2*(Tf - pair.phase2().thermo().T())
-                  - mDotL
+                  - dmdtLf
                 )
                /(
-                    max(H1 + H2 - mDotLPrime, HSmall)
+                    max(H1 + H2 - dmdtLfPrime, HSmall)
                 );
 
             Tf.correctBoundaryConditions();

applications/solvers/multiphase/reactingEulerFoam/phaseSystems/PhaseSystems/InterfaceCompositionPhaseChangePhaseSystem/InterfaceCompositionPhaseChangePhaseSystem.H

@@ -59,8 +59,6 @@ class InterfaceCompositionPhaseChangePhaseSystem
 :
     public BasePhaseSystem
 {
-private:
-
     // Private typedefs
 
         typedef HashTable
@@ -92,19 +90,19 @@ private:
             interfaceCompositionModelTable interfaceCompositionModels_;
 
             //- The explicit part of the interfacial mass transfer rates
-            phaseSystem::dmidtTable iDmdtSu_;
+            phaseSystem::dmidtfTable dmidtfSus_;
 
             //- The implicit part of the interfacial mass transfer rates
-            phaseSystem::dmidtTable iDmdtSp_;
+            phaseSystem::dmidtfTable dmidtfSps_;
 
 
     // Private member functions
 
-        //- Return the bulk mass transfers across each interface
+        //- Return mass transfers across each interface
-        autoPtr<phaseSystem::dmdtTable> iDmdts() const;
+        autoPtr<phaseSystem::dmdtfTable> dmdtfs() const;
 
-        //- Return the species mass transfers across each interface
+        //- Return species mass transfers across each interface
-        autoPtr<phaseSystem::dmidtTable> iDmidts() const;
+        autoPtr<phaseSystem::dmidtfTable> dmidtfs() const;
 
 
 public:
@@ -121,6 +119,9 @@ public:
 
     // Member Functions
 
+        //- Return the mass transfer rate for an interface
+        virtual tmp<volScalarField> dmdtf(const phasePairKey& key) const;
+
         //- Return the mass transfer rates for each phase
         virtual PtrList<volScalarField> dmdts() const;
 

applications/solvers/multiphase/reactingEulerFoam/phaseSystems/PhaseSystems/MomentumTransferPhaseSystem/MomentumTransferPhaseSystem.C

@@ -119,18 +119,18 @@ Foam::MomentumTransferPhaseSystem<BasePhaseSystem>::Vm
 // * * * * * * * * * * * * Protected Member Functions * * * * * * * * * * * //
 
 template<class BasePhaseSystem>
-void Foam::MomentumTransferPhaseSystem<BasePhaseSystem>::addDmdtU
+void Foam::MomentumTransferPhaseSystem<BasePhaseSystem>::addDmdtUfs
 (
-    const phaseSystem::dmdtTable& dmdts,
+    const phaseSystem::dmdtfTable& dmdtfs,
     phaseSystem::momentumTransferTable& eqns
 )
 {
-    forAllConstIter(phaseSystem::dmdtTable, dmdts, dmdtIter)
+    forAllConstIter(phaseSystem::dmdtfTable, dmdtfs, dmdtfIter)
     {
-        const phasePairKey& key = dmdtIter.key();
+        const phasePairKey& key = dmdtfIter.key();
         const phasePair& pair(this->phasePairs_[key]);
 
-        const volScalarField dmdt(Pair<word>::compare(pair, key)**dmdtIter());
+        const volScalarField dmdtf(Pair<word>::compare(pair, key)**dmdtfIter());
 
         phaseModel& phase1 = this->phases()[pair.phase1().name()];
         phaseModel& phase2 = this->phases()[pair.phase2().name()];
@@ -150,18 +150,18 @@ void Foam::MomentumTransferPhaseSystem<BasePhaseSystem>::addDmdtU
 
         if (!phase1.stationary())
         {
-            const volScalarField dmdt21(posPart(dmdt));
+            const volScalarField dmdtf21(posPart(dmdtf));
 
             *eqns[phase1.name()] +=
-                dmdt21*phase2.U() - fvm::Sp(dmdt21, phase1.URef());
+                dmdtf21*phase2.U() - fvm::Sp(dmdtf21, phase1.URef());
         }
 
         if (!phase2.stationary())
         {
-            const volScalarField dmdt12(negPart(dmdt));
+            const volScalarField dmdtf12(negPart(dmdtf));
 
             *eqns[phase2.name()] -=
-                dmdt12*phase1.U() - fvm::Sp(dmdt12, phase2.URef());
+                dmdtf12*phase1.U() - fvm::Sp(dmdtf12, phase2.URef());
         }
     }
 }
@@ -220,13 +220,33 @@ MomentumTransferPhaseSystem
         Kds_.insert
         (
             pair,
-            zeroVolField<scalar>(pair, "Kd", dragModel::dimK).ptr()
+            new volScalarField
+            (
+                IOobject
+                (
+                    IOobject::groupName("Kd", pair.name()),
+                    this->mesh().time().timeName(),
+                    this->mesh()
+                ),
+                this->mesh(),
+                dimensionedScalar(dragModel::dimK, 0)
+            )
         );
 
         Kdfs_.insert
         (
             pair,
-            zeroSurfaceField<scalar>(pair, "Kdf", dragModel::dimK).ptr()
+            new surfaceScalarField
+            (
+                IOobject
+                (
+                    IOobject::groupName("Kdf", pair.name()),
+                    this->mesh().time().timeName(),
+                    this->mesh()
+                ),
+                this->mesh(),
+                dimensionedScalar(dragModel::dimK, 0)
+            )
         );
     }
 
@@ -242,13 +262,33 @@ MomentumTransferPhaseSystem
         Vms_.insert
         (
             pair,
-            zeroVolField<scalar>(pair, "Vm", virtualMassModel::dimK).ptr()
+            new volScalarField
+            (
+                IOobject
+                (
+                    IOobject::groupName("Vm", pair.name()),
+                    this->mesh().time().timeName(),
+                    this->mesh()
+                ),
+                this->mesh(),
+                dimensionedScalar(virtualMassModel::dimK, 0)
+            )
         );
 
         Vmfs_.insert
         (
             pair,
-            zeroSurfaceField<scalar>(pair, "Vmf", virtualMassModel::dimK).ptr()
+            new surfaceScalarField
+            (
+                IOobject
+                (
+                    IOobject::groupName("Vmf", pair.name()),
+                    this->mesh().time().timeName(),
+                    this->mesh()
+                ),
+                this->mesh(),
+                dimensionedScalar(virtualMassModel::dimK, 0)
+            )
         );
     }
 }

applications/solvers/multiphase/reactingEulerFoam/phaseSystems/PhaseSystems/MomentumTransferPhaseSystem/MomentumTransferPhaseSystem.H

@@ -192,9 +192,9 @@ protected:
     // Protected Member Functions
 
         //- Add momentum transfer terms which result from bulk mass transfers
-        void addDmdtU
+        void addDmdtUfs
         (
-            const phaseSystem::dmdtTable& dmdts,
+            const phaseSystem::dmdtfTable& dmdtfs,
             phaseSystem::momentumTransferTable& eqns
         );
 

applications/solvers/multiphase/reactingEulerFoam/phaseSystems/PhaseSystems/OneResistanceHeatTransferPhaseSystem/OneResistanceHeatTransferPhaseSystem.C

@@ -32,20 +32,19 @@ License
 // * * * * * * * * * * * * Protected Member Functions * * * * * * * * * * * //
 
 template<class BasePhaseSystem>
-void Foam::OneResistanceHeatTransferPhaseSystem<BasePhaseSystem>::addDmdtHe
+void Foam::OneResistanceHeatTransferPhaseSystem<BasePhaseSystem>::addDmdtHefs
 (
-    const phaseSystem::dmdtTable& dmdts,
+    const phaseSystem::dmdtfTable& dmdtfs,
     phaseSystem::heatTransferTable& eqns
 ) const
 {
-    forAllConstIter(phaseSystem::dmdtTable, dmdts, dmdtIter)
+    forAllConstIter(phaseSystem::dmdtfTable, dmdtfs, dmdtfIter)
     {
-        const phasePairKey& key = dmdtIter.key();
+        const phasePairKey& key = dmdtfIter.key();
         const phasePair& pair(this->phasePairs_[key]);
-
+        const volScalarField dmdtf(Pair<word>::compare(pair, key)**dmdtfIter());
-        const volScalarField dmdt(Pair<word>::compare(pair, key)**dmdtIter());
+        const volScalarField dmdtf21(posPart(dmdtf));
-        const volScalarField dmdt21(posPart(dmdt));
+        const volScalarField dmdtf12(negPart(dmdtf));
-        const volScalarField dmdt12(negPart(dmdt));
 
         const phaseModel& phase1 = pair.phase1();
         const phaseModel& phase2 = pair.phase2();
@@ -57,30 +56,30 @@ void Foam::OneResistanceHeatTransferPhaseSystem<BasePhaseSystem>::addDmdtHe
         const volScalarField K2(phase2.K());
 
         // Note that the phase EEqn contains a continuity error term. See
-        // MomentumTransferPhaseSystem::addDmdtU for an explanation of the
+        // MomentumTransferPhaseSystem::addDmdtUfs for an explanation of the
         // fvm::Sp terms below.
 
         // Transfer of energy from bulk to bulk
-        *eqns[phase1.name()] += dmdt21*he2 - fvm::Sp(dmdt21, he1);
+        *eqns[phase1.name()] += dmdtf21*he2 - fvm::Sp(dmdtf21, he1);
-        *eqns[phase2.name()] -= dmdt12*he1 - fvm::Sp(dmdt12, he2);
+        *eqns[phase2.name()] -= dmdtf12*he1 - fvm::Sp(dmdtf12, he2);
 
         // Transfer of kinetic energy
-        *eqns[phase1.name()] += dmdt21*(K2 - K1);
+        *eqns[phase1.name()] += dmdtf21*(K2 - K1);
-        *eqns[phase2.name()] -= dmdt12*(K1 - K2);
+        *eqns[phase2.name()] -= dmdtf12*(K1 - K2);
     }
 }
 
 
 template<class BasePhaseSystem>
-void Foam::OneResistanceHeatTransferPhaseSystem<BasePhaseSystem>::addDmidtHe
+void Foam::OneResistanceHeatTransferPhaseSystem<BasePhaseSystem>::addDmidtHef
 (
-    const phaseSystem::dmidtTable& dmidts,
+    const phaseSystem::dmidtfTable& dmidtfs,
     phaseSystem::heatTransferTable& eqns
 ) const
 {
-    forAllConstIter(phaseSystem::dmidtTable, dmidts, dmidtIter)
+    forAllConstIter(phaseSystem::dmidtfTable, dmidtfs, dmidtfIter)
     {
-        const phasePairKey& key = dmidtIter.key();
+        const phasePairKey& key = dmidtfIter.key();
         const phasePair& pair(this->phasePairs_[key]);
 
         const phaseModel& phase1 = pair.phase1();
@@ -93,19 +92,19 @@ void Foam::OneResistanceHeatTransferPhaseSystem<BasePhaseSystem>::addDmidtHe
         const volScalarField K2(phase2.K());
 
         // Note that the phase EEqn contains a continuity error term. See
-        // MomentumTransferPhaseSystem::addDmdtU for an explanation of the
+        // MomentumTransferPhaseSystem::addDmdtUfs for an explanation of the
         // fvm::Sp terms below.
 
-        forAllConstIter(HashPtrTable<volScalarField>, *dmidtIter(), dmidtJter)
+        forAllConstIter(HashPtrTable<volScalarField>, *dmidtfIter(), dmidtfJter)
         {
-            const word& member = dmidtJter.key();
+            const word& member = dmidtfJter.key();
 
-            const volScalarField dmidt
+            const volScalarField dmidtf
             (
-                Pair<word>::compare(pair, key)**dmidtJter()
+                Pair<word>::compare(pair, key)**dmidtfJter()
             );
-            const volScalarField dmidt21(posPart(dmidt));
+            const volScalarField dmidtf21(posPart(dmidtf));
-            const volScalarField dmidt12(negPart(dmidt));
+            const volScalarField dmidtf12(negPart(dmidtf));
 
             // Create the energies for the transferring specie
             volScalarField hei1(he1);
@@ -136,12 +135,12 @@ void Foam::OneResistanceHeatTransferPhaseSystem<BasePhaseSystem>::addDmidtHe
             }
 
             // Transfer of energy from bulk to bulk
-            *eqns[phase1.name()] += dmidt21*hei2 - fvm::Sp(dmidt21, he1);
+            *eqns[phase1.name()] += dmidtf21*hei2 - fvm::Sp(dmidtf21, he1);
-            *eqns[phase2.name()] -= dmidt12*hei1 - fvm::Sp(dmidt12, he2);
+            *eqns[phase2.name()] -= dmidtf12*hei1 - fvm::Sp(dmidtf12, he2);
 
             // Transfer of kinetic energy
-            *eqns[phase1.name()] += dmidt21*(K2 - K1);
+            *eqns[phase1.name()] += dmidtf21*(K2 - K1);
-            *eqns[phase2.name()] -= dmidt12*(K1 - K2);
+            *eqns[phase2.name()] -= dmidtf12*(K1 - K2);
         }
     }
 }

applications/solvers/multiphase/reactingEulerFoam/phaseSystems/PhaseSystems/OneResistanceHeatTransferPhaseSystem/OneResistanceHeatTransferPhaseSystem.H

@@ -82,16 +82,16 @@ protected:
     // Protected Member Functions
 
         //- Add energy transfer terms which result from bulk mass transfers
-        void addDmdtHe
+        void addDmdtHefs
         (
-            const phaseSystem::dmdtTable& dmdts,
+            const phaseSystem::dmdtfTable& dmdtfs,
             phaseSystem::heatTransferTable& eqns
         ) const;
 
         //- Add energy transfer terms which result from specie mass transfers
-        void addDmidtHe
+        void addDmidtHef
         (
-            const phaseSystem::dmidtTable& dmidts,
+            const phaseSystem::dmidtfTable& dmidtfs,
             phaseSystem::heatTransferTable& eqns
         ) const;
 

applications/solvers/multiphase/reactingEulerFoam/phaseSystems/PhaseSystems/PhaseTransferPhaseSystem/PhaseTransferPhaseSystem.C

@@ -30,66 +30,14 @@ License
 // * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * * //
 
 template<class BasePhaseSystem>
-void Foam::PhaseTransferPhaseSystem<BasePhaseSystem>::addToSpeciesDmdt
+Foam::autoPtr<Foam::phaseSystem::dmdtfTable>
-(
+Foam::PhaseTransferPhaseSystem<BasePhaseSystem>::totalDmdtfs() const
-    const phasePair& pair,
-    const wordList species
-)
 {
-    if (species.empty())
+    autoPtr<phaseSystem::dmdtfTable> totalDmdtfsPtr
-    {
-        return;
-    }
-
-    if (!specieDmdt_.found(pair))
-    {
-        specieDmdt_.insert(pair, new HashPtrTable<volScalarField>());
-    }
-
-    forAllConstIter
     (
-        wordList,
+        new phaseSystem::dmdtfTable
-        species,
+    );
-        sIter
+    phaseSystem::dmdtfTable& totalDmdtfs = totalDmdtfsPtr();
-    )
-    {
-        const word specieName(*sIter);
-        if (!specieDmdt_[pair]->found(specieName))
-        {
-            specieDmdt_[pair]->insert
-            (
-                specieName,
-                new volScalarField
-                (
-                    IOobject
-                    (
-                        IOobject::groupName
-                        (
-                            specieName+":dmdt",
-                            pair.name()
-                        ),
-                        this->mesh().time().timeName(),
-                        this->mesh(),
-                        IOobject::NO_READ,
-                        IOobject::AUTO_WRITE
-                    ),
-                    this->mesh(),
-                    dimensionedScalar(Foam::phaseTransferModel::dimDmdt, 0)
-                )
-            );
-        }
-    }
-}
-
-
-template<class BasePhaseSystem>
-Foam::autoPtr<Foam::phaseSystem::dmdtTable>
-Foam::PhaseTransferPhaseSystem<BasePhaseSystem>::
-sumDmdts() const
-{
-    autoPtr<phaseSystem::dmdtTable> sumDmdtsPtr(new phaseSystem::dmdtTable);
-
-    phaseSystem::dmdtTable& sumDmdts = sumDmdtsPtr();
 
     forAllConstIter
     (
@@ -101,52 +49,45 @@ sumDmdts() const
         const phasePair& pair =
             this->phasePairs_[phaseTransferModelIter.key()];
 
-        sumDmdts.insert
+        totalDmdtfs.insert(pair, phaseSystem::dmdtf(pair).ptr());
-        (
-            pair,
-            zeroVolField<scalar>(pair, "dmdt", dimDensity/dimTime).ptr()
-        );
 
-        if (specieDmdt_.found(pair))
+        if (phaseTransferModelIter()->mixture())
         {
-            forAllConstIter
+            *totalDmdtfs[pair] += *dmdtfs_[pair];
-            (
-                HashPtrTable<volScalarField>,
-                *specieDmdt_[pair],
-                sDmdtIter
-            )
-            {
-                *sumDmdts[pair] +=
-                    phaseTransferModelIter()->speciesDmdt(sDmdtIter.key());
-            }
         }
-        else
+
+        forAllConstIter
+        (
+            HashPtrTable<volScalarField>,
+            *dmidtfs_[pair],
+            dmidtfIter
+        )
         {
-            *sumDmdts[pair] += *rDmdt_[pair];
+            *totalDmdtfs[pair] += *dmidtfIter();
         }
     }
 
-    return sumDmdtsPtr;
+    return totalDmdtfsPtr;
 }
 
 
 // * * * * * * * * * * * * Protected Member Functions * * * * * * * * * * * //
 
 template<class BasePhaseSystem>
-void Foam::PhaseTransferPhaseSystem<BasePhaseSystem>::addDmdtY
+void Foam::PhaseTransferPhaseSystem<BasePhaseSystem>::addDmdtYfs
 (
-    const phaseSystem::dmdtTable& dmdts,
+    const phaseSystem::dmdtfTable& dmdtfs,
     phaseSystem::specieTransferTable& eqns
 ) const
 {
-    forAllConstIter(phaseSystem::dmdtTable, dmdts, dmdtIter)
+    forAllConstIter(phaseSystem::dmdtfTable, dmdtfs, dmdtfIter)
     {
-        const phasePairKey& key = dmdtIter.key();
+        const phasePairKey& key = dmdtfIter.key();
         const phasePair& pair(this->phasePairs_[key]);
 
-        const volScalarField dmdt(Pair<word>::compare(pair, key)**dmdtIter());
+        const volScalarField dmdtf(Pair<word>::compare(pair, key)**dmdtfIter());
-        const volScalarField dmdt12(negPart(dmdt));
+        const volScalarField dmdtf12(negPart(dmdtf));
-        const volScalarField dmdt21(posPart(dmdt));
+        const volScalarField dmdtf21(posPart(dmdtf));
 
         const phaseModel& phase1 = pair.phase1();
         const phaseModel& phase2 = pair.phase2();
@@ -159,23 +100,23 @@ void Foam::PhaseTransferPhaseSystem<BasePhaseSystem>::addDmdtY
             const volScalarField& Y1 = phase1.Y()[Yi1];
             const volScalarField& Y2 = phase2.Y(Y1.member());
 
-            *eqns[Y1.name()] += dmdt21*Y2 + fvm::Sp(dmdt12, Y1);
+            *eqns[Y1.name()] += dmdtf21*Y2 + fvm::Sp(dmdtf12, Y1);
-            *eqns[Y2.name()] -= dmdt12*Y1 + fvm::Sp(dmdt21, Y2);
+            *eqns[Y2.name()] -= dmdtf12*Y1 + fvm::Sp(dmdtf21, Y2);
         }
     }
 }
 
 
 template<class BasePhaseSystem>
-void Foam::PhaseTransferPhaseSystem<BasePhaseSystem>::addDmidtY
+void Foam::PhaseTransferPhaseSystem<BasePhaseSystem>::addDmidtYf
 (
-    const phaseSystem::dmidtTable& dmidts,
+    const phaseSystem::dmidtfTable& dmidtfs,
     phaseSystem::specieTransferTable& eqns
 ) const
 {
-    forAllConstIter(phaseSystem::dmidtTable, dmidts, dmidtIter)
+    forAllConstIter(phaseSystem::dmidtfTable, dmidtfs, dmidtfIter)
     {
-        const phasePairKey& key = dmidtIter.key();
+        const phasePairKey& key = dmidtfIter.key();
         const phasePair& pair(this->phasePairs_[key]);
 
         const phaseModel& phase1 = pair.phase1();
@@ -184,25 +125,25 @@ void Foam::PhaseTransferPhaseSystem<BasePhaseSystem>::addDmidtY
         // Note that the phase YiEqn does not contain a continuity error term,
         // so the transfers below are complete.
 
-        forAllConstIter(HashPtrTable<volScalarField>, *dmidtIter(), dmidtJter)
+        forAllConstIter(HashPtrTable<volScalarField>, *dmidtfIter(), dmidtfJter)
         {
-            const word& member = dmidtJter.key();
+            const word& member = dmidtfJter.key();
 
-            const volScalarField dmidt
+            const volScalarField dmidtf
             (
-                Pair<word>::compare(pair, key)**dmidtJter()
+                Pair<word>::compare(pair, key)**dmidtfJter()
             );
 
             if (!phase1.pure())
             {
                 const volScalarField& Y1 = phase1.Y(member);
-                *eqns[Y1.name()] += dmidt;
+                *eqns[Y1.name()] += dmidtf;
             }
 
             if (!phase2.pure())
             {
                 const volScalarField& Y2 = phase2.Y(member);
-                *eqns[Y2.name()] -= dmidt;
+                *eqns[Y2.name()] -= dmidtf;
             }
         }
     }
@@ -235,14 +176,59 @@ Foam::PhaseTransferPhaseSystem<BasePhaseSystem>::PhaseTransferPhaseSystem
     {
         const phasePair& pair = this->phasePairs_[phaseTransferModelIter.key()];
 
-        addToSpeciesDmdt(pair, phaseTransferModelIter()->activeSpecies());
+        if (phaseTransferModelIter()->mixture())
-
-        if (!specieDmdt_.found(pair))
         {
-            this->rDmdt_.insert
+            dmdtfs_.insert
             (
-                phaseTransferModelIter.key(),
+                pair,
-                zeroVolField<scalar>(pair, "rDmdt", dimDensity/dimTime).ptr()
+                new volScalarField
+                (
+                    IOobject
+                    (
+                        IOobject::groupName
+                        (
+                            "phaseTransfer:dmdtf",
+                            pair.name()
+                        ),
+                        this->mesh().time().timeName(),
+                        this->mesh()
+                    ),
+                    this->mesh(),
+                    dimensionedScalar(dimDensity/dimTime, 0)
+                )
+            );
+        }
+
+        dmidtfs_.insert(pair, new HashPtrTable<volScalarField>());
+
+        const hashedWordList species(phaseTransferModelIter()->species());
+
+        forAllConstIter(hashedWordList, species, specieIter)
+        {
+            const word& specie = *specieIter;
+
+            dmidtfs_[pair]->insert
+            (
+                specie,
+                new volScalarField
+                (
+                    IOobject
+                    (
+                        IOobject::groupName
+                        (
+                            IOobject::groupName
+                            (
+                                "phaseTransfer:dmidtf",
+                                specie
+                            ),
+                            pair.name()
+                        ),
+                        this->mesh().time().timeName(),
+                        this->mesh()
+                    ),
+                    this->mesh(),
+                    dimensionedScalar(dimDensity/dimTime, 0)
+                )
             );
         }
     }
@@ -252,30 +238,60 @@ Foam::PhaseTransferPhaseSystem<BasePhaseSystem>::PhaseTransferPhaseSystem
 // * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
 
 template<class BasePhaseSystem>
-Foam::PhaseTransferPhaseSystem<BasePhaseSystem>::
+Foam::PhaseTransferPhaseSystem<BasePhaseSystem>::~PhaseTransferPhaseSystem()
-~PhaseTransferPhaseSystem()
 {}
 
 
 // * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
 
+template<class BasePhaseSystem>
+Foam::tmp<Foam::volScalarField>
+Foam::PhaseTransferPhaseSystem<BasePhaseSystem>::dmdtf
+(
+    const phasePairKey& key
+) const
+{
+    tmp<volScalarField> tDmdtf = BasePhaseSystem::dmdtf(key);
+
+    if (phaseTransferModels_.found(key))
+    {
+        const label dmdtfSign(Pair<word>::compare(this->phasePairs_[key], key));
+
+        if (phaseTransferModels_[key]->mixture())
+        {
+            tDmdtf.ref() += dmdtfSign**dmdtfs_[key];
+        }
+
+        forAllConstIter
+        (
+            HashPtrTable<volScalarField>,
+            *dmidtfs_[key],
+            dmidtfIter
+        )
+        {
+            tDmdtf.ref() += dmdtfSign**dmidtfIter();
+        }
+    }
+
+    return tDmdtf;
+}
+
+
 template<class BasePhaseSystem>
 Foam::PtrList<Foam::volScalarField>
 Foam::PhaseTransferPhaseSystem<BasePhaseSystem>::dmdts() const
 {
     PtrList<volScalarField> dmdts(BasePhaseSystem::dmdts());
 
-    autoPtr<phaseSystem::dmdtTable> sumDmdtsPtr = this->sumDmdts();
+    autoPtr<phaseSystem::dmdtfTable> totalDmdtfsPtr = this->totalDmdtfs();
+    const phaseSystem::dmdtfTable& totalDmdtfs = totalDmdtfsPtr();
 
-    const phaseSystem::dmdtTable& sumDmdts = sumDmdtsPtr();
+    forAllConstIter(phaseSystem::dmdtfTable, totalDmdtfs, totalDmdtfIter)
-
-    forAllConstIter(phaseSystem::dmdtTable, sumDmdts, sumDmdtIter)
     {
-        const phasePair& pair = this->phasePairs_[sumDmdtIter.key()];
+        const phasePair& pair = this->phasePairs_[totalDmdtfIter.key()];
-
-        addField(pair.phase1(), "dmdt", *sumDmdtIter(), dmdts);
-        addField(pair.phase2(), "dmdt", - *sumDmdtIter(), dmdts);
 
+        addField(pair.phase1(), "dmdt", *totalDmdtfIter(), dmdts);
+        addField(pair.phase2(), "dmdt", - *totalDmdtfIter(), dmdts);
     }
 
     return dmdts;
@@ -291,7 +307,7 @@ Foam::PhaseTransferPhaseSystem<BasePhaseSystem>::momentumTransfer()
 
     phaseSystem::momentumTransferTable& eqns = eqnsPtr();
 
-    this->addDmdtU(sumDmdts(), eqns);
+    this->addDmdtUfs(totalDmdtfs(), eqns);
 
     return eqnsPtr;
 }
@@ -306,7 +322,7 @@ Foam::PhaseTransferPhaseSystem<BasePhaseSystem>::momentumTransferf()
 
     phaseSystem::momentumTransferTable& eqns = eqnsPtr();
 
-    this->addDmdtU(sumDmdts(), eqns);
+    this->addDmdtUfs(totalDmdtfs(), eqns);
 
     return eqnsPtr;
 }
@@ -321,9 +337,8 @@ Foam::PhaseTransferPhaseSystem<BasePhaseSystem>::heatTransfer() const
 
     phaseSystem::heatTransferTable& eqns = eqnsPtr();
 
-    this->addDmdtHe(rDmdt_, eqns);
+    this->addDmdtHefs(dmdtfs_, eqns);
-    this->addDmidtHe(specieDmdt_, eqns);
+    this->addDmidtHef(dmidtfs_, eqns);
-
 
     return eqnsPtr;
 }
@@ -357,8 +372,8 @@ Foam::PhaseTransferPhaseSystem<BasePhaseSystem>::specieTransfer() const
         }
     }
 
-    this->addDmdtY(rDmdt_, eqns);
+    this->addDmdtYfs(dmdtfs_, eqns);
-    this->addDmidtY(specieDmdt_, eqns);
+    this->addDmidtYf(dmidtfs_, eqns);
 
     return eqnsPtr;
 }
@@ -369,6 +384,7 @@ void Foam::PhaseTransferPhaseSystem<BasePhaseSystem>::correct()
 {
     BasePhaseSystem::correct();
 
+    // Reset all the mass transfer rates to zero
     forAllConstIter
     (
         phaseTransferModelTable,
@@ -376,25 +392,24 @@ void Foam::PhaseTransferPhaseSystem<BasePhaseSystem>::correct()
         phaseTransferModelIter
     )
     {
-        if (specieDmdt_.found(phaseTransferModelIter.key()))
+        const phasePair& pair = this->phasePairs_[phaseTransferModelIter.key()];
+
+        if (phaseTransferModelIter()->mixture())
         {
-            forAllIter
+            *dmdtfs_[pair] = Zero;
-            (
-                HashPtrTable<volScalarField>,
-                *specieDmdt_(phaseTransferModelIter.key()),
-                sDmdtIter
-            )
-            {
-                (*sDmdtIter)->ref() = Zero;
-            }
         }
-        else
+
+        const hashedWordList species(phaseTransferModelIter()->species());
+
+        forAllConstIter(hashedWordList, species, specieIter)
         {
-            *rDmdt_[phaseTransferModelIter.key()] =
+            const word& specie = *specieIter;
-                dimensionedScalar(dimDensity/dimTime, 0);
+
+            *(*dmidtfs_[pair])[specie] = Zero;
         }
     }
 
+    // Evaluate the models and sum the results into the mass transfer tables
     forAllIter
     (
         phaseTransferModelTable,
@@ -402,25 +417,21 @@ void Foam::PhaseTransferPhaseSystem<BasePhaseSystem>::correct()
         phaseTransferModelIter
     )
     {
-        phaseTransferModelIter()->correct();
+        const phasePair& pair = this->phasePairs_[phaseTransferModelIter.key()];
 
-        if (specieDmdt_.found(phaseTransferModelIter.key()))
+        if (phaseTransferModelIter()->mixture())
         {
-            forAllIter
+            *dmdtfs_[pair] += phaseTransferModelIter()->dmdtf();
-            (
-                HashPtrTable<volScalarField>,
-                *specieDmdt_(phaseTransferModelIter.key()),
-                sDmdtIter
-            )
-            {
-                (*sDmdtIter)->ref() +=
-                    phaseTransferModelIter()->speciesDmdt(sDmdtIter.key());
-            }
         }
-        else
+
+        const HashPtrTable<volScalarField> dmidtf
+        (
+            phaseTransferModelIter()->dmidtf()
+        );
+
+        forAllConstIter(HashPtrTable<volScalarField>, dmidtf, dmidtfIter)
         {
-            *rDmdt_[phaseTransferModelIter.key()] +=
+            *(*dmidtfs_[pair])[dmidtfIter.key()] += *dmidtfIter();
-                phaseTransferModelIter()->dmdt();
         }
     }
 }

applications/solvers/multiphase/reactingEulerFoam/phaseSystems/PhaseSystems/PhaseTransferPhaseSystem/PhaseTransferPhaseSystem.H

@@ -76,22 +76,16 @@ private:
             phaseTransferModelTable phaseTransferModels_;
 
             //- Bulk mass transfer rates
-            phaseSystem::dmdtTable rDmdt_;
+            phaseSystem::dmdtfTable dmdtfs_;
 
             //- Specie mass transfer rates
-            phaseSystem::dmidtTable specieDmdt_;
+            phaseSystem::dmidtfTable dmidtfs_;
+
 
     // Private member functions
 
-        //- Add given species to the speciesDmdt for the given pair
+        //- Return the summed mass transfers across each interface
-        virtual void addToSpeciesDmdt
+        autoPtr<phaseSystem::dmdtfTable> totalDmdtfs() const;
-        (
-            const phasePair& key,
-            const wordList species
-        );
-
-        //- Return the sum of specie transfers across each interface
-        autoPtr<phaseSystem::dmdtTable> sumDmdts() const;
 
 
 protected:
@@ -99,16 +93,16 @@ protected:
     // Protected member functions
 
         //- Add specie transfer terms which result from bulk mass transfers
-        void addDmdtY
+        void addDmdtYfs
         (
-            const phaseSystem::dmdtTable& dmdts,
+            const phaseSystem::dmdtfTable& dmdtfs,
             phaseSystem::specieTransferTable& eqns
         ) const;
 
         //- Add specie transfer terms which result from specie mass transfers
-        void addDmidtY
+        void addDmidtYf
         (
-            const phaseSystem::dmidtTable& dmidts,
+            const phaseSystem::dmidtfTable& dmidtfs,
             phaseSystem::specieTransferTable& eqns
         ) const;
 
@@ -127,6 +121,9 @@ public:
 
     // Member Functions
 
+        //- Return the mass transfer rate for an interface
+        virtual tmp<volScalarField> dmdtf(const phasePairKey& key) const;
+
         //- Return the mass transfer rates for each phase
         virtual PtrList<volScalarField> dmdts() const;
 

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

@@ -25,27 +25,6 @@ License
 
 #include "PopulationBalancePhaseSystem.H"
 
-// * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * * //
-
-template<class BasePhaseSystem>
-Foam::tmp<Foam::volScalarField>
-Foam::PopulationBalancePhaseSystem<BasePhaseSystem>::pDmdt
-(
-    const phasePairKey& key
-) const
-{
-    if (!pDmdt_.found(key))
-    {
-        const phasePair& pair = this->phasePairs_[key];
-        return zeroVolField<scalar>(pair, "pDmdt", dimDensity/dimTime);
-    }
-
-    const scalar pDmdtSign(Pair<word>::compare(pDmdt_.find(key).key(), key));
-
-    return pDmdtSign**pDmdt_[key];
-}
-
-
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
 template<class BasePhaseSystem>
@@ -60,7 +39,7 @@ PopulationBalancePhaseSystem
     populationBalances_
     (
         this->lookup("populationBalances"),
-        diameterModels::populationBalanceModel::iNew(*this, pDmdt_)
+        diameterModels::populationBalanceModel::iNew(*this, dmdtfs_)
     )
 {
     forAll(populationBalances_, i)
@@ -88,7 +67,7 @@ PopulationBalancePhaseSystem
                 );
             }
 
-            pDmdt_.insert
+            dmdtfs_.insert
             (
                 key,
                 new volScalarField
@@ -97,7 +76,7 @@ PopulationBalancePhaseSystem
                     (
                         IOobject::groupName
                         (
-                            "pDmdt",
+                            "populationBalance:dmdtf",
                             this->phasePairs_[key]->name()
                         ),
                         this->mesh().time().timeName(),
@@ -124,16 +103,36 @@ Foam::PopulationBalancePhaseSystem<BasePhaseSystem>::
 
 // * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
 
+template<class BasePhaseSystem>
+Foam::tmp<Foam::volScalarField>
+Foam::PopulationBalancePhaseSystem<BasePhaseSystem>::dmdtf
+(
+    const phasePairKey& key
+) const
+{
+    tmp<volScalarField> tDmdt = BasePhaseSystem::dmdtf(key);
+
+    if (!dmdtfs_.found(key))
+    {
+        const label dmdtSign(Pair<word>::compare(this->phasePairs_[key], key));
+
+        tDmdt.ref() += dmdtSign**dmdtfs_[key];
+    }
+
+    return tDmdt;
+}
+
+
 template<class BasePhaseSystem>
 Foam::PtrList<Foam::volScalarField>
 Foam::PopulationBalancePhaseSystem<BasePhaseSystem>::dmdts() const
 {
     PtrList<volScalarField> dmdts(BasePhaseSystem::dmdts());
 
-    forAllConstIter(phaseSystem::dmdtTable, pDmdt_, pDmdtIter)
+    forAllConstIter(phaseSystem::dmdtfTable, dmdtfs_, dmdtfIter)
     {
-        const phasePair& pair = this->phasePairs_[pDmdtIter.key()];
+        const phasePair& pair = this->phasePairs_[dmdtfIter.key()];
-        const volScalarField& pDmdt = *pDmdtIter();
+        const volScalarField& pDmdt = *dmdtfIter();
 
         addField(pair.phase1(), "dmdt", pDmdt, dmdts);
         addField(pair.phase2(), "dmdt", - pDmdt, dmdts);
@@ -152,7 +151,7 @@ Foam::PopulationBalancePhaseSystem<BasePhaseSystem>::momentumTransfer()
 
     phaseSystem::momentumTransferTable& eqns = eqnsPtr();
 
-    this->addDmdtU(pDmdt_, eqns);
+    this->addDmdtUfs(dmdtfs_, eqns);
 
     return eqnsPtr;
 }
@@ -167,7 +166,7 @@ Foam::PopulationBalancePhaseSystem<BasePhaseSystem>::momentumTransferf()
 
     phaseSystem::momentumTransferTable& eqns = eqnsPtr();
 
-    this->addDmdtU(pDmdt_, eqns);
+    this->addDmdtUfs(dmdtfs_, eqns);
 
     return eqnsPtr;
 }
@@ -182,7 +181,7 @@ Foam::PopulationBalancePhaseSystem<BasePhaseSystem>::heatTransfer() const
 
     phaseSystem::heatTransferTable& eqns = eqnsPtr();
 
-    this->addDmdtHe(pDmdt_, eqns);
+    this->addDmdtHefs(dmdtfs_, eqns);
 
     return eqnsPtr;
 }
@@ -197,7 +196,7 @@ Foam::PopulationBalancePhaseSystem<BasePhaseSystem>::specieTransfer() const
 
     phaseSystem::specieTransferTable& eqns = eqnsPtr();
 
-    this->addDmdtY(pDmdt_, eqns);
+    this->addDmdtYfs(dmdtfs_, eqns);
 
     return eqnsPtr;
 }

applications/solvers/multiphase/reactingEulerFoam/phaseSystems/PhaseSystems/PopulationBalancePhaseSystem/PopulationBalancePhaseSystem.H

@@ -55,21 +55,13 @@ class PopulationBalancePhaseSystem
 :
     public BasePhaseSystem
 {
-protected:
+    // Private data
-
-    // Protected data
 
         //- Population balances
         PtrList<diameterModels::populationBalanceModel> populationBalances_;
 
         //- Mass transfer rates
-        phaseSystem::dmdtTable pDmdt_;
+        phaseSystem::dmdtfTable dmdtfs_;
-
-
-    // Protected member functions
-
-        //- Return the population balance induced mass transfer rate
-        virtual tmp<volScalarField> pDmdt(const phasePairKey& key) const;
 
 
 public:
@@ -86,6 +78,9 @@ public:
 
     // Member Functions
 
+        //- Return the mass transfer rate for an interface
+        virtual tmp<volScalarField> dmdtf(const phasePairKey& key) const;
+
         //- Return the mass transfer rates for each phase
         virtual PtrList<volScalarField> dmdts() const;
 

applications/solvers/multiphase/reactingEulerFoam/phaseSystems/PhaseSystems/ThermalPhaseChangePhaseSystem/ThermalPhaseChangePhaseSystem.C

@@ -32,39 +32,28 @@ License
 
 template<class BasePhaseSystem>
 Foam::tmp<Foam::volScalarField>
-Foam::ThermalPhaseChangePhaseSystem<BasePhaseSystem>::iDmdt
+Foam::ThermalPhaseChangePhaseSystem<BasePhaseSystem>::totalDmdtf
 (
     const phasePairKey& key
 ) const
 {
-    if (!iDmdt_.found(key))
+    tmp<volScalarField> tTotalDmdtf = phaseSystem::dmdtf(key);
+    volScalarField& totalDmdtf = tTotalDmdtf.ref();
+
+    const scalar tTotalDmdtfSign =
+        Pair<word>::compare(this->phasePairs_[key], key);
+
+    if (dmdtfs_.found(key))
     {
-        const phasePair& pair = this->phasePairs_[key];
+        totalDmdtf += *dmdtfs_[key];
-        return zeroVolField<scalar>(pair, "iDmdt", dimDensity/dimTime);
     }
 
-    const scalar dmdtSign(Pair<word>::compare(iDmdt_.find(key).key(), key));
+    if (nDmdtfs_.found(key))
-
-    return dmdtSign**iDmdt_[key];
-}
-
-
-template<class BasePhaseSystem>
-Foam::tmp<Foam::volScalarField>
-Foam::ThermalPhaseChangePhaseSystem<BasePhaseSystem>::wDmdt
-(
-    const phasePairKey& key
-) const
-{
-    if (!wDmdt_.found(key))
     {
-        const phasePair& pair = this->phasePairs_[key];
+        totalDmdtf += *nDmdtfs_[key];
-        return zeroVolField<scalar>(pair, "wDmdt", dimDensity/dimTime);
     }
 
-    const scalar dmdtSign(Pair<word>::compare(wDmdt_.find(key).key(), key));
+    return tTotalDmdtfSign*tTotalDmdtf;
-
-    return dmdtSign**wDmdt_[key];
 }
 
 
@@ -85,7 +74,6 @@ ThermalPhaseChangePhaseSystem
     ),
     phaseChange_(this->lookup("phaseChange"))
 {
-
     forAllConstIter
     (
         phaseSystem::phasePairTable,
@@ -101,14 +89,18 @@ ThermalPhaseChangePhaseSystem
         }
 
         // Initially assume no mass transfer
-        iDmdt_.insert
+        dmdtfs_.insert
         (
             pair,
             new volScalarField
             (
                 IOobject
                 (
-                    IOobject::groupName("iDmdt", pair.name()),
+                    IOobject::groupName
+                    (
+                        "thermalPhaseChange:dmdtf",
+                        pair.name()
+                    ),
                     this->mesh().time().timeName(),
                     this->mesh(),
                     IOobject::READ_IF_PRESENT,
@@ -120,14 +112,18 @@ ThermalPhaseChangePhaseSystem
         );
 
         // Initially assume no mass transfer
-        wDmdt_.insert
+        nDmdtfs_.insert
         (
             pair,
             new volScalarField
             (
                 IOobject
                 (
-                    IOobject::groupName("wDmdt", pair.name()),
+                    IOobject::groupName
+                    (
+                        "thermalPhaseChange:nucleation:dmdtf",
+                        pair.name()
+                    ),
                     this->mesh().time().timeName(),
                     this->mesh(),
                     IOobject::READ_IF_PRESENT,
@@ -139,14 +135,18 @@ ThermalPhaseChangePhaseSystem
         );
 
         // Initially assume no mass transfer
-        wMDotL_.insert
+        nDmdtLfs_.insert
         (
             pair,
             new volScalarField
             (
                 IOobject
                 (
-                    IOobject::groupName("wMDotL", pair.name()),
+                    IOobject::groupName
+                    (
+                        "thermalPhaseChange:nucleation:dmdtLf",
+                        pair.name()
+                    ),
                     this->mesh().time().timeName(),
                     this->mesh(),
                     IOobject::READ_IF_PRESENT,
@@ -178,28 +178,39 @@ Foam::ThermalPhaseChangePhaseSystem<BasePhaseSystem>::saturation() const
 }
 
 
+template<class BasePhaseSystem>
+Foam::tmp<Foam::volScalarField>
+Foam::ThermalPhaseChangePhaseSystem<BasePhaseSystem>::dmdtf
+(
+    const phasePairKey& key
+) const
+{
+    return BasePhaseSystem::dmdtf(key) + this->totalDmdtf(key);
+}
+
+
 template<class BasePhaseSystem>
 Foam::PtrList<Foam::volScalarField>
 Foam::ThermalPhaseChangePhaseSystem<BasePhaseSystem>::dmdts() const
 {
     PtrList<volScalarField> dmdts(BasePhaseSystem::dmdts());
 
-    forAllConstIter(iDmdtTable, iDmdt_, iDmdtIter)
+    forAllConstIter(phaseSystem::dmdtfTable, dmdtfs_, dmdtfIter)
     {
-        const phasePair& pair = this->phasePairs_[iDmdtIter.key()];
+        const phasePair& pair = this->phasePairs_[dmdtfIter.key()];
-        const volScalarField& iDmdt = *iDmdtIter();
+        const volScalarField& dmdtf = *dmdtfIter();
 
-        addField(pair.phase1(), "dmdt", iDmdt, dmdts);
+        addField(pair.phase1(), "dmdt", dmdtf, dmdts);
-        addField(pair.phase2(), "dmdt", - iDmdt, dmdts);
+        addField(pair.phase2(), "dmdt", - dmdtf, dmdts);
     }
 
-    forAllConstIter(wDmdtTable, wDmdt_, wDmdtIter)
+    forAllConstIter(phaseSystem::dmdtfTable, nDmdtfs_, nDmdtfIter)
     {
-        const phasePair& pair = this->phasePairs_[wDmdtIter.key()];
+        const phasePair& pair = this->phasePairs_[nDmdtfIter.key()];
-        const volScalarField& wDmdt = *wDmdtIter();
+        const volScalarField& nDmdtf = *nDmdtfIter();
 
-        addField(pair.phase1(), "dmdt", wDmdt, dmdts);
+        addField(pair.phase1(), "dmdt", nDmdtf, dmdts);
-        addField(pair.phase2(), "dmdt", - wDmdt, dmdts);
+        addField(pair.phase2(), "dmdt", - nDmdtf, dmdts);
     }
 
     return dmdts;
@@ -216,8 +227,8 @@ momentumTransfer()
 
     phaseSystem::momentumTransferTable& eqns = eqnsPtr();
 
-    this->addDmdtU(iDmdt_, eqns);
+    this->addDmdtUfs(dmdtfs_, eqns);
-    this->addDmdtU(wDmdt_, eqns);
+    this->addDmdtUfs(nDmdtfs_, eqns);
 
     return eqnsPtr;
 }
@@ -233,8 +244,8 @@ momentumTransferf()
 
     phaseSystem::momentumTransferTable& eqns = eqnsPtr();
 
-    this->addDmdtU(iDmdt_, eqns);
+    this->addDmdtUfs(dmdtfs_, eqns);
-    this->addDmdtU(wDmdt_, eqns);
+    this->addDmdtUfs(nDmdtfs_, eqns);
 
     return eqnsPtr;
 }
@@ -302,35 +313,35 @@ Foam::ThermalPhaseChangePhaseSystem<BasePhaseSystem>::heatTransfer() const
         const volScalarField K1(phase1.K());
         const volScalarField K2(phase2.K());
 
-        const volScalarField dmdt(this->iDmdt(pair) + this->wDmdt(pair));
+        const volScalarField dmdtf(this->totalDmdtf(pair));
-        const volScalarField dmdt21(posPart(dmdt));
+        const volScalarField dmdtf21(posPart(dmdtf));
-        const volScalarField dmdt12(negPart(dmdt));
+        const volScalarField dmdtf12(negPart(dmdtf));
 
-        *eqns[phase1.name()] += - fvm::Sp(dmdt21, he1) + dmdt21*(K2 - K1);
+        *eqns[phase1.name()] += - fvm::Sp(dmdtf21, he1) + dmdtf21*(K2 - K1);
 
-        *eqns[phase2.name()] -= - fvm::Sp(dmdt12, he2) + dmdt12*(K1 - K2);
+        *eqns[phase2.name()] -= - fvm::Sp(dmdtf12, he2) + dmdtf12*(K1 - K2);
 
         if (this->heatTransferModels_.found(phasePairIter.key()))
         {
             const volScalarField& Tf(*this->Tf_[pair]);
 
             *eqns[phase1.name()] +=
-                dmdt21*phase1.thermo().he(phase1.thermo().p(), Tf);
+                dmdtf21*phase1.thermo().he(phase1.thermo().p(), Tf);
 
             *eqns[phase2.name()] -=
-                dmdt12*phase2.thermo().he(phase2.thermo().p(), Tf);
+                dmdtf12*phase2.thermo().he(phase2.thermo().p(), Tf);
         }
         else
         {
-            *eqns[phase1.name()] += dmdt21*he2;
+            *eqns[phase1.name()] += dmdtf21*he2;
 
-            *eqns[phase2.name()] -= dmdt12*he1;
+            *eqns[phase2.name()] -= dmdtf12*he1;
         }
 
-        if (this->wMDotL_.found(phasePairIter.key()))
+        if (this->nDmdtLfs_.found(phasePairIter.key()))
         {
-            *eqns[phase1.name()] += negPart(*this->wMDotL_[pair]);
+            *eqns[phase1.name()] += negPart(*this->nDmdtLfs_[pair]);
-            *eqns[phase2.name()] -= posPart(*this->wMDotL_[pair]);
+            *eqns[phase2.name()] -= posPart(*this->nDmdtLfs_[pair]);
 
             if
             (
@@ -341,13 +352,13 @@ Foam::ThermalPhaseChangePhaseSystem<BasePhaseSystem>::heatTransfer() const
                 if (phase1.thermo().he().member() == "e")
                 {
                     *eqns[phase1.name()] +=
-                        phase1.thermo().p()*dmdt/phase1.thermo().rho();
+                        phase1.thermo().p()*dmdtf/phase1.thermo().rho();
                 }
 
                 if (phase2.thermo().he().member() == "e")
                 {
                     *eqns[phase2.name()] -=
-                        phase2.thermo().p()*dmdt/phase2.thermo().rho();
+                        phase2.thermo().p()*dmdtf/phase2.thermo().rho();
                 }
             }
         }
@@ -413,10 +424,10 @@ Foam::ThermalPhaseChangePhaseSystem<BasePhaseSystem>::specieTransfer() const
             // Note that the phase YiEqn does not contain a continuity error
             // term, so these additions represent the entire mass transfer
 
-            const volScalarField dmdt(this->iDmdt(pair) + this->wDmdt(pair));
+            const volScalarField dmdtf(this->totalDmdtf(pair));
 
-            *eqns[name] += dmdt;
+            *eqns[name] += dmdtf;
-            *eqns[otherName] -= dmdt;
+            *eqns[otherName] -= dmdtf;
         }
     }
 
@@ -454,7 +465,7 @@ Foam::ThermalPhaseChangePhaseSystem<BasePhaseSystem>::correctInterfaceThermo()
 
         const volScalarField& p(phase1.thermo().p());
 
-        volScalarField& iDmdt(*this->iDmdt_[pair]);
+        volScalarField& dmdtf(*this->dmdtfs_[pair]);
         volScalarField& Tf(*this->Tf_[pair]);
 
         const volScalarField Tsat(saturationModel_->Tsat(phase1.thermo().p()));
@@ -488,22 +499,22 @@ Foam::ThermalPhaseChangePhaseSystem<BasePhaseSystem>::correctInterfaceThermo()
 
         volScalarField L
         (
-            (neg0(iDmdt)*hf2 + pos(iDmdt)*h2)
+            (neg0(dmdtf)*hf2 + pos(dmdtf)*h2)
-          - (pos0(iDmdt)*hf1 + neg(iDmdt)*h1)
+          - (pos0(dmdtf)*hf1 + neg(dmdtf)*h1)
         );
 
-        volScalarField iDmdtNew(iDmdt);
+        volScalarField dmdtfNew(dmdtf);
 
         if (phaseChange_)
         {
             volScalarField H1(heatTransferModelIter().first()->K(0));
             volScalarField H2(heatTransferModelIter().second()->K(0));
 
-            iDmdtNew = (H1*(Tsat - T1) + H2*(Tsat - T2))/L;
+            dmdtfNew = (H1*(Tsat - T1) + H2*(Tsat - T2))/L;
         }
         else
         {
-            iDmdtNew == dimensionedScalar(iDmdt.dimensions(), 0);
+            dmdtfNew == dimensionedScalar(dmdtf.dimensions(), 0);
         }
 
         volScalarField H1(heatTransferModelIter().first()->K());
@@ -513,7 +524,7 @@ Foam::ThermalPhaseChangePhaseSystem<BasePhaseSystem>::correctInterfaceThermo()
         H1.max(small);
         H2.max(small);
 
-        Tf = (H1*T1 + H2*T2 + iDmdtNew*L)/(H1 + H2);
+        Tf = (H1*T1 + H2*T2 + dmdtfNew*L)/(H1 + H2);
 
         Info<< "Tf." << pair.name()
             << ": min = " << min(Tf.primitiveField())
@@ -521,23 +532,24 @@ Foam::ThermalPhaseChangePhaseSystem<BasePhaseSystem>::correctInterfaceThermo()
             << ", max = " << max(Tf.primitiveField())
             << endl;
 
-        scalar iDmdtRelax(this->mesh().fieldRelaxationFactor("iDmdt"));
+        scalar dmdtfRelax =
-        iDmdt = (1 - iDmdtRelax)*iDmdt + iDmdtRelax*iDmdtNew;
+            this->mesh().fieldRelaxationFactor(dmdtf.member());
+        dmdtf = (1 - dmdtfRelax)*dmdtf + dmdtfRelax*dmdtfNew;
 
         if (phaseChange_)
         {
-            Info<< "iDmdt." << pair.name()
+            Info<< dmdtf.name()
-                << ": min = " << min(iDmdt.primitiveField())
+                << ": min = " << min(dmdtf.primitiveField())
-                << ", mean = " << average(iDmdt.primitiveField())
+                << ", mean = " << average(dmdtf.primitiveField())
-                << ", max = " << max(iDmdt.primitiveField())
+                << ", max = " << max(dmdtf.primitiveField())
-                << ", integral = " << fvc::domainIntegrate(iDmdt).value()
+                << ", integral = " << fvc::domainIntegrate(dmdtf).value()
                 << endl;
         }
 
-        volScalarField& wDmdt(*this->wDmdt_[pair]);
+        volScalarField& nDmdtf(*this->nDmdtfs_[pair]);
-        volScalarField& wMDotL(*this->wMDotL_[pair]);
+        volScalarField& nDmdtLf(*this->nDmdtLfs_[pair]);
-        wDmdt = Zero;
+        nDmdtf = Zero;
-        wMDotL = Zero;
+        nDmdtLf = Zero;
 
         bool wallBoilingActive = false;
 
@@ -550,14 +562,14 @@ Foam::ThermalPhaseChangePhaseSystem<BasePhaseSystem>::correctInterfaceThermo()
             (
                 phase.mesh().foundObject<volScalarField>
                 (
-                    "alphat." +  phase.name()
+                    IOobject::groupName("alphat", phase.name())
                 )
             )
             {
                 const volScalarField& alphat =
                     phase.mesh().lookupObject<volScalarField>
                     (
-                        "alphat." +  phase.name()
+                        IOobject::groupName("alphat", phase.name())
                     );
 
                 const fvPatchList& patches = this->mesh().boundary();
@@ -585,22 +597,21 @@ Foam::ThermalPhaseChangePhaseSystem<BasePhaseSystem>::correctInterfaceThermo()
                         {
                             wallBoilingActive = true;
 
-                            const scalarField& patchDmdt =
+                            const scalarField& patchDmdtf = PCpatch.dmdtf(key);
-                                PCpatch.dmdt(key);
+                            const scalarField& patchDmdtLf =
-                            const scalarField& patchMDotL =
+                                PCpatch.dmdtLf(key);
-                                PCpatch.mDotL(key);
 
                             const scalar sign
                             (
                                 Pair<word>::compare(pair, key)
                             );
 
-                            forAll(patchDmdt, facei)
+                            forAll(patchDmdtf, facei)
                             {
                                 const label faceCelli =
                                     currPatch.faceCells()[facei];
-                                wDmdt[faceCelli] -= sign*patchDmdt[facei];
+                                nDmdtf[faceCelli] -= sign*patchDmdtf[facei];
-                                wMDotL[faceCelli] -= sign*patchMDotL[facei];
+                                nDmdtLf[faceCelli] -= sign*patchDmdtLf[facei];
                             }
                         }
                     }
@@ -610,11 +621,11 @@ Foam::ThermalPhaseChangePhaseSystem<BasePhaseSystem>::correctInterfaceThermo()
 
         if (wallBoilingActive)
         {
-            Info<< "wDmdt." << pair.name()
+            Info<< nDmdtf.name()
-                << ": min = " << min(wDmdt.primitiveField())
+                << ": min = " << min(nDmdtf.primitiveField())
-                << ", mean = " << average(wDmdt.primitiveField())
+                << ", mean = " << average(nDmdtf.primitiveField())
-                << ", max = " << max(wDmdt.primitiveField())
+                << ", max = " << max(nDmdtf.primitiveField())
-                << ", integral = " << fvc::domainIntegrate(wDmdt).value()
+                << ", integral = " << fvc::domainIntegrate(nDmdtf).value()
                 << endl;
         }
     }

applications/solvers/multiphase/reactingEulerFoam/phaseSystems/PhaseSystems/ThermalPhaseChangePhaseSystem/ThermalPhaseChangePhaseSystem.H

@@ -55,21 +55,7 @@ class ThermalPhaseChangePhaseSystem
 :
     public BasePhaseSystem
 {
-protected:
+    // Private data
-
-    // Protected typedefs
-
-        typedef HashPtrTable<volScalarField, phasePairKey, phasePairKey::hash>
-            iDmdtTable;
-
-        typedef HashPtrTable<volScalarField, phasePairKey, phasePairKey::hash>
-            wDmdtTable;
-
-        typedef HashPtrTable<volScalarField, phasePairKey, phasePairKey::hash>
-            wMDotLTable;
-
-
-    // Protected data
 
         //- Name of the volatile specie
         word volatile_;
@@ -80,23 +66,20 @@ protected:
         // Phase change enabled
         Switch phaseChange_;
 
-        //- Interfacial Mass transfer rate
+        //- Mass transfer rates
-        iDmdtTable iDmdt_;
+        phaseSystem::dmdtfTable dmdtfs_;
 
-        //- Boundary Mass transfer rate
+        //- Nucleate Mass transfer rates
-        wDmdtTable wDmdt_;
+        phaseSystem::dmdtfTable nDmdtfs_;
 
-        //- Boundary thermal energy transfer rate
+        //- Nucleate thermal energy transfer rates
-        wMDotLTable wMDotL_;
+        phaseSystem::dmdtfTable nDmdtLfs_;
 
 
-    // Protected member functions
+    // Private member functions
 
-        //- Return the interfacial mass transfer rate for a pair
+        //- Return the total mass transfer rate for a pair
-        tmp<volScalarField> iDmdt(const phasePairKey& key) const;
+        tmp<volScalarField> totalDmdtf(const phasePairKey& key) const;
-
-        //- Return the boundary mass transfer rate for a pair
-        tmp<volScalarField> wDmdt(const phasePairKey& key) const;
 
 
 public:
@@ -116,6 +99,9 @@ public:
         //- Return the saturationModel
         const saturationModel& saturation() const;
 
+        //- Return the mass transfer rate for an interface
+        virtual tmp<volScalarField> dmdtf(const phasePairKey& key) const;
+
         //- Return the mass transfer rates for each phase
         virtual PtrList<volScalarField> dmdts() const;
 

applications/solvers/multiphase/reactingEulerFoam/phaseSystems/PhaseSystems/TwoResistanceHeatTransferPhaseSystem/TwoResistanceHeatTransferPhaseSystem.C

@@ -32,9 +32,9 @@ License
 // * * * * * * * * * * * * Protected Member Functions * * * * * * * * * * * //
 
 template<class BasePhaseSystem>
-void Foam::TwoResistanceHeatTransferPhaseSystem<BasePhaseSystem>::addDmdtHe
+void Foam::TwoResistanceHeatTransferPhaseSystem<BasePhaseSystem>::addDmdtHefs
 (
-    const phaseSystem::dmdtTable& dmdts,
+    const phaseSystem::dmdtfTable& dmdtfs,
     phaseSystem::heatTransferTable& eqns
 ) const
 {
@@ -62,14 +62,14 @@ void Foam::TwoResistanceHeatTransferPhaseSystem<BasePhaseSystem>::addDmdtHe
     // latent heat mDot*L terms are added below (as well as other standard
     // terms representing sensible energy and kinetic energy differences).
 
-    forAllConstIter(phaseSystem::dmdtTable, dmdts, dmdtIter)
+    forAllConstIter(phaseSystem::dmdtfTable, dmdtfs, dmdtfIter)
     {
-        const phasePairKey& key = dmdtIter.key();
+        const phasePairKey& key = dmdtfIter.key();
         const phasePair& pair(this->phasePairs_[key]);
 
-        const volScalarField dmdt(Pair<word>::compare(pair, key)**dmdtIter());
+        const volScalarField dmdtf(Pair<word>::compare(pair, key)**dmdtfIter());
-        const volScalarField dmdt21(posPart(dmdt));
+        const volScalarField dmdtf21(posPart(dmdtf));
-        const volScalarField dmdt12(negPart(dmdt));
+        const volScalarField dmdtf12(negPart(dmdtf));
 
         const phaseModel& phase1 = pair.phase1();
         const phaseModel& phase2 = pair.phase2();
@@ -94,20 +94,20 @@ void Foam::TwoResistanceHeatTransferPhaseSystem<BasePhaseSystem>::addDmdtHe
             const volScalarField& Tf(*Tf_[key]);
             const volScalarField hef1(thermo1.he(thermo1.p(), Tf));
             const volScalarField hef2(thermo2.he(thermo2.p(), Tf));
-            *eqns[phase1.name()] += dmdt*hef1 - fvm::Sp(dmdt, he1);
+            *eqns[phase1.name()] += dmdtf*hef1 - fvm::Sp(dmdtf, he1);
-            *eqns[phase2.name()] -= dmdt*hef2 - fvm::Sp(dmdt, he2);
+            *eqns[phase2.name()] -= dmdtf*hef2 - fvm::Sp(dmdtf, he2);
 
             // Latent heat contribution
             const volScalarField L(hef2 + thermo2.hc() - hef1 - thermo1.hc());
             const volScalarField H1(heatTransferModels_[key].first()->K());
             const volScalarField H2(heatTransferModels_[key].second()->K());
             const volScalarField H1Fac(H1/(H1 + H2));
-            *eqns[phase1.name()] += H1Fac*dmdt*L;
+            *eqns[phase1.name()] += H1Fac*dmdtf*L;
-            *eqns[phase2.name()] += (1 - H1Fac)*dmdt*L;
+            *eqns[phase2.name()] += (1 - H1Fac)*dmdtf*L;
 
             // Transfer of kinetic energy
-            *eqns[phase1.name()] += dmdt21*(K2 - K1);
+            *eqns[phase1.name()] += dmdtf21*(K2 - K1);
-            *eqns[phase2.name()] -= dmdt12*(K1 - K2);
+            *eqns[phase2.name()] -= dmdtf12*(K1 - K2);
         }
         else
         {
@@ -116,29 +116,29 @@ void Foam::TwoResistanceHeatTransferPhaseSystem<BasePhaseSystem>::addDmdtHe
             // change and therefore no interface state or latent heat...
 
             // Transfer of energy from bulk to bulk
-            *eqns[phase1.name()] += dmdt21*he2 - fvm::Sp(dmdt21, he1);
+            *eqns[phase1.name()] += dmdtf21*he2 - fvm::Sp(dmdtf21, he1);
-            *eqns[phase2.name()] -= dmdt12*he1 - fvm::Sp(dmdt12, he2);
+            *eqns[phase2.name()] -= dmdtf12*he1 - fvm::Sp(dmdtf12, he2);
 
             // Transfer of kinetic energy
-            *eqns[phase1.name()] += dmdt21*(K2 - K1);
+            *eqns[phase1.name()] += dmdtf21*(K2 - K1);
-            *eqns[phase2.name()] -= dmdt12*(K1 - K2);
+            *eqns[phase2.name()] -= dmdtf12*(K1 - K2);
         }
     }
 }
 
 
 template<class BasePhaseSystem>
-void Foam::TwoResistanceHeatTransferPhaseSystem<BasePhaseSystem>::addDmidtHe
+void Foam::TwoResistanceHeatTransferPhaseSystem<BasePhaseSystem>::addDmidtHef
 (
-    const phaseSystem::dmidtTable& dmidts,
+    const phaseSystem::dmidtfTable& dmidtfs,
     phaseSystem::heatTransferTable& eqns
 ) const
 {
     // See the addDmdtHe method for a description of the latent heat terms
 
-    forAllConstIter(phaseSystem::dmidtTable, dmidts, dmidtIter)
+    forAllConstIter(phaseSystem::dmidtfTable, dmidtfs, dmidtfIter)
     {
-        const phasePairKey& key = dmidtIter.key();
+        const phasePairKey& key = dmidtfIter.key();
         const phasePair& pair(this->phasePairs_[key]);
 
         const phaseModel& phase1 = pair.phase1();
@@ -176,18 +176,18 @@ void Foam::TwoResistanceHeatTransferPhaseSystem<BasePhaseSystem>::addDmidtHe
             forAllConstIter
             (
                 HashPtrTable<volScalarField>,
-                *dmidtIter(),
+                *dmidtfIter(),
-                dmidtJter
+                dmidtfJter
             )
             {
-                const word& member = dmidtJter.key();
+                const word& member = dmidtfJter.key();
 
-                const volScalarField dmidt
+                const volScalarField dmidtf
                 (
-                    Pair<word>::compare(pair, key)**dmidtJter()
+                    Pair<word>::compare(pair, key)**dmidtfJter()
                 );
-                const volScalarField dmidt21(posPart(dmidt));
+                const volScalarField dmidtf21(posPart(dmidtf));
-                const volScalarField dmidt12(negPart(dmidt));
+                const volScalarField dmidtf12(negPart(dmidtf));
 
                 // Create the interface energies for the transferring specie
                 volScalarField hefi1(hef1), hci1(hc1);
@@ -233,16 +233,16 @@ void Foam::TwoResistanceHeatTransferPhaseSystem<BasePhaseSystem>::addDmidtHe
                 const volScalarField Li(hefi2 + hci2 - hefi1 - hci1);
 
                 // Transfer of energy from the interface into the bulk
-                *eqns[phase1.name()] += dmidt*hefi1 - fvm::Sp(dmidt, he1);
+                *eqns[phase1.name()] += dmidtf*hefi1 - fvm::Sp(dmidtf, he1);
-                *eqns[phase2.name()] -= dmidt*hefi2 - fvm::Sp(dmidt, he2);
+                *eqns[phase2.name()] -= dmidtf*hefi2 - fvm::Sp(dmidtf, he2);
 
                 // Latent heat contribution
-                *eqns[phase1.name()] += H1Fac*dmidt*Li;
+                *eqns[phase1.name()] += H1Fac*dmidtf*Li;
-                *eqns[phase2.name()] += (1 - H1Fac)*dmidt*Li;
+                *eqns[phase2.name()] += (1 - H1Fac)*dmidtf*Li;
 
                 // Transfer of kinetic energy
-                *eqns[phase1.name()] += dmidt21*(K2 - K1);
+                *eqns[phase1.name()] += dmidtf21*(K2 - K1);
-                *eqns[phase2.name()] -= dmidt12*(K1 - K2);
+                *eqns[phase2.name()] -= dmidtf12*(K1 - K2);
             }
         }
         else
@@ -255,18 +255,18 @@ void Foam::TwoResistanceHeatTransferPhaseSystem<BasePhaseSystem>::addDmidtHe
             forAllConstIter
             (
                 HashPtrTable<volScalarField>,
-                *dmidtIter(),
+                *dmidtfIter(),
-                dmidtJter
+                dmidtfJter
             )
             {
-                const word& member = dmidtJter.key();
+                const word& member = dmidtfJter.key();
 
-                const volScalarField dmidt
+                const volScalarField dmidtf
                 (
-                    Pair<word>::compare(pair, key)**dmidtJter()
+                    Pair<word>::compare(pair, key)**dmidtfJter()
                 );
-                const volScalarField dmidt21(posPart(dmidt));
+                const volScalarField dmidtf21(posPart(dmidtf));
-                const volScalarField dmidt12(negPart(dmidt));
+                const volScalarField dmidtf12(negPart(dmidtf));
 
                 // Create the energies for the transferring specie
                 volScalarField hei1(he1);
@@ -297,12 +297,12 @@ void Foam::TwoResistanceHeatTransferPhaseSystem<BasePhaseSystem>::addDmidtHe
                 }
 
                 // Transfer of energy from bulk to bulk
-                *eqns[phase1.name()] += dmidt21*hei2 - fvm::Sp(dmidt21, he1);
+                *eqns[phase1.name()] += dmidtf21*hei2 - fvm::Sp(dmidtf21, he1);
-                *eqns[phase2.name()] -= dmidt12*hei1 - fvm::Sp(dmidt12, he2);
+                *eqns[phase2.name()] -= dmidtf12*hei1 - fvm::Sp(dmidtf12, he2);
 
                 // Transfer of kinetic energy
-                *eqns[phase1.name()] += dmidt21*(K2 - K1);
+                *eqns[phase1.name()] += dmidtf21*(K2 - K1);
-                *eqns[phase2.name()] -= dmidt12*(K1 - K2);
+                *eqns[phase2.name()] -= dmidtf12*(K1 - K2);
             }
         }
     }

applications/solvers/multiphase/reactingEulerFoam/phaseSystems/PhaseSystems/TwoResistanceHeatTransferPhaseSystem/TwoResistanceHeatTransferPhaseSystem.H

@@ -90,16 +90,16 @@ protected:
     // Protected Member Functions
 
         //- Add energy transfer terms which result from bulk mass transfers
-        void addDmdtHe
+        void addDmdtHefs
         (
-            const phaseSystem::dmdtTable& dmdts,
+            const phaseSystem::dmdtfTable& dmdtfs,
             phaseSystem::heatTransferTable& eqns
         ) const;
 
         //- Add energy transfer terms which result from specie mass transfers
-        void addDmidtHe
+        void addDmidtHef
         (
-            const phaseSystem::dmidtTable& dmidts,
+            const phaseSystem::dmidtfTable& dmidtfs,
             phaseSystem::heatTransferTable& eqns
         ) const;
 

applications/solvers/multiphase/reactingEulerFoam/phaseSystems/phaseModel/StationaryPhaseModel/StationaryPhaseModel.C

@@ -83,7 +83,12 @@ template<class BasePhaseModel>
 Foam::tmp<Foam::volVectorField>
 Foam::StationaryPhaseModel<BasePhaseModel>::U() const
 {
-    return zeroVolField<vector>(*this, "U", dimVelocity);
+    return volVectorField::New
+    (
+        IOobject::groupName("U", this->name()),
+        this->mesh(),
+        dimensionedVector(dimVelocity, Zero)
+    );
 }
 
 
@@ -103,7 +108,12 @@ template<class BasePhaseModel>
 Foam::tmp<Foam::surfaceScalarField>
 Foam::StationaryPhaseModel<BasePhaseModel>::phi() const
 {
-    return zeroSurfaceField<scalar>(*this, "phi", dimVolume/dimTime);
+    return surfaceScalarField::New
+    (
+        IOobject::groupName("phi", this->name()),
+        this->mesh(),
+        dimensionedScalar(dimVolume/dimTime, 0)
+    );
 }
 
 
@@ -123,7 +133,12 @@ template<class BasePhaseModel>
 Foam::tmp<Foam::surfaceScalarField>
 Foam::StationaryPhaseModel<BasePhaseModel>::alphaPhi() const
 {
-    return zeroSurfaceField<scalar>(*this, "alphaPhi", dimVolume/dimTime);
+    return surfaceScalarField::New
+    (
+        IOobject::groupName("alphaPhi", this->name()),
+        this->mesh(),
+        dimensionedScalar(dimVolume/dimTime, 0)
+    );
 }
 
 
@@ -143,7 +158,12 @@ template<class BasePhaseModel>
 Foam::tmp<Foam::surfaceScalarField>
 Foam::StationaryPhaseModel<BasePhaseModel>::alphaRhoPhi() const
 {
-    return zeroSurfaceField<scalar>(*this, "alphaRhoPhi", dimMass/dimTime);
+    return surfaceScalarField::New
+    (
+        IOobject::groupName("alphaRhoPhi", this->name()),
+        this->mesh(),
+        dimensionedScalar(dimMass/dimTime, 0)
+    );
 }
 
 
@@ -163,7 +183,12 @@ template<class BasePhaseModel>
 Foam::tmp<Foam::volVectorField>
 Foam::StationaryPhaseModel<BasePhaseModel>::DUDt() const
 {
-    return zeroVolField<vector>(*this, "DUDt", dimVelocity/dimTime);
+    return volVectorField::New
+    (
+        IOobject::groupName("DUDt", this->name()),
+        this->mesh(),
+        dimensionedVector(dimVelocity/dimTime, Zero)
+    );
 }
 
 
@@ -171,7 +196,12 @@ template<class BasePhaseModel>
 Foam::tmp<Foam::surfaceScalarField>
 Foam::StationaryPhaseModel<BasePhaseModel>::DUDtf() const
 {
-    return zeroSurfaceField<scalar>(*this, "DUDtf", dimVolume/sqr(dimTime));
+    return surfaceScalarField::New
+    (
+        IOobject::groupName("DUDtf", this->name()),
+        this->mesh(),
+        dimensionedScalar(dimVolume/sqr(dimTime), 0)
+    );
 }
 
 
@@ -179,7 +209,12 @@ template<class BasePhaseModel>
 Foam::tmp<Foam::volScalarField>
 Foam::StationaryPhaseModel<BasePhaseModel>::continuityError() const
 {
-    return zeroVolField<scalar>(*this, "contErr", dimDensity/dimTime);
+    return volScalarField::New
+    (
+        IOobject::groupName("continuityError", this->name()),
+        this->mesh(),
+        dimensionedScalar(dimDensity/dimTime, 0)
+    );
 }
 
 
@@ -187,7 +222,12 @@ template<class BasePhaseModel>
 Foam::tmp<Foam::volScalarField>
 Foam::StationaryPhaseModel<BasePhaseModel>::continuityErrorFlow() const
 {
-    return zeroVolField<scalar>(*this, "contErrFlow", dimDensity/dimTime);
+    return volScalarField::New
+    (
+        IOobject::groupName("continuityErrorFlow", this->name()),
+        this->mesh(),
+        dimensionedScalar(dimDensity/dimTime, 0)
+    );
 }
 
 
@@ -195,7 +235,12 @@ template<class BasePhaseModel>
 Foam::tmp<Foam::volScalarField>
 Foam::StationaryPhaseModel<BasePhaseModel>::continuityErrorSources() const
 {
-    return zeroVolField<scalar>(*this, "contErrSources", dimDensity/dimTime);
+    return volScalarField::New
+    (
+        IOobject::groupName("continuityErrorSources", this->name()),
+        this->mesh(),
+        dimensionedScalar(dimDensity/dimTime, 0)
+    );
 }
 
 
@@ -203,7 +248,12 @@ template<class BasePhaseModel>
 Foam::tmp<Foam::volScalarField>
 Foam::StationaryPhaseModel<BasePhaseModel>::K() const
 {
-    return zeroVolField<scalar>(*this, "K", sqr(dimVelocity));
+    return volScalarField::New
+    (
+        IOobject::groupName("K", this->name()),
+        this->mesh(),
+        dimensionedScalar(sqr(dimVelocity), 0)
+    );
 }
 
 
@@ -231,7 +281,12 @@ template<class BasePhaseModel>
 Foam::tmp<Foam::volScalarField>
 Foam::StationaryPhaseModel<BasePhaseModel>::mut() const
 {
-    return zeroVolField<scalar>(*this, "mut", dimDynamicViscosity);
+    return volScalarField::New
+    (
+        IOobject::groupName("mut", this->name()),
+        this->mesh(),
+        dimensionedScalar(dimDynamicViscosity, 0)
+    );
 }
 
 
@@ -247,7 +302,12 @@ template<class BasePhaseModel>
 Foam::tmp<Foam::volScalarField>
 Foam::StationaryPhaseModel<BasePhaseModel>::nut() const
 {
-    return zeroVolField<scalar>(*this, "nut", dimViscosity);
+    return volScalarField::New
+    (
+        IOobject::groupName("nut", this->name()),
+        this->mesh(),
+        dimensionedScalar(dimViscosity, 0)
+    );
 }
 
 
@@ -295,7 +355,12 @@ template<class BasePhaseModel>
 Foam::tmp<Foam::volScalarField>
 Foam::StationaryPhaseModel<BasePhaseModel>::k() const
 {
-    return zeroVolField<scalar>(*this, "k", sqr(dimVelocity));
+    return volScalarField::New
+    (
+        IOobject::groupName("k", this->name()),
+        this->mesh(),
+        dimensionedScalar(sqr(dimVelocity), 0)
+    );
 }
 
 
@@ -303,7 +368,12 @@ template<class BasePhaseModel>
 Foam::tmp<Foam::volScalarField>
 Foam::StationaryPhaseModel<BasePhaseModel>::pPrime() const
 {
-    return zeroVolField<scalar>(*this, "pPrime", dimPressure);
+    return volScalarField::New
+    (
+        IOobject::groupName("pPrime", this->name()),
+        this->mesh(),
+        dimensionedScalar(dimPressure, 0)
+    );
 }
 
 

applications/solvers/multiphase/reactingEulerFoam/phaseSystems/phaseModel/phaseModel/phaseModel.H

@@ -222,12 +222,6 @@ public:
 
         // Thermo
 
-            //- Return whether the phase is isothermal
-            virtual bool isothermal() const = 0;
-
-            //- Return the enthalpy equation
-            virtual tmp<fvScalarMatrix> heEqn() = 0;
-
             //- Return the thermophysical model
             virtual const rhoThermo& thermo() const = 0;
 
@@ -237,6 +231,12 @@ public:
             //- Return the density field
             virtual tmp<volScalarField> rho() const = 0;
 
+            //- Return whether the phase is isothermal
+            virtual bool isothermal() const = 0;
+
+            //- Return the enthalpy equation
+            virtual tmp<fvScalarMatrix> heEqn() = 0;
+
 
         // Species
 

applications/solvers/multiphase/reactingEulerFoam/phaseSystems/phaseSystem/phaseSystem.C

@@ -302,7 +302,7 @@ Foam::phaseSystem::E(const phasePairKey& key) const
         return volScalarField::New
         (
             aspectRatioModel::typeName + ":E",
-            this->mesh_,
+            mesh_,
             dimensionedScalar(dimless, 1)
         );
     }
@@ -321,7 +321,7 @@ Foam::phaseSystem::sigma(const phasePairKey& key) const
         return volScalarField::New
         (
             surfaceTensionModel::typeName + ":sigma",
-            this->mesh_,
+            mesh_,
             dimensionedScalar(surfaceTensionModel::dimSigma, 0)
         );
     }
@@ -339,17 +339,36 @@ Foam::phaseSystem::sigma(const phasePairKey& key, label patchi) const
     {
         return tmp<scalarField>
         (
-            new scalarField(this->mesh_.boundary()[patchi].size(), 0)
+            new scalarField(mesh_.boundary()[patchi].size(), 0)
         );
     }
 }
 
 
+Foam::tmp<Foam::volScalarField> Foam::phaseSystem::dmdtf
+(
+    const phasePairKey& key
+) const
+{
+    const phasePair pair
+    (
+        phaseModels_[key.first()],
+        phaseModels_[key.second()]
+    );
+
+    return volScalarField::New
+    (
+        IOobject::groupName("dmdtf", pair.name()),
+        mesh(),
+        dimensionedScalar(dimDensity/dimTime, 0)
+    );
+}
+
+
+
 Foam::PtrList<Foam::volScalarField> Foam::phaseSystem::dmdts() const
 {
-    PtrList<volScalarField> dmdts(this->phaseModels_.size());
+    return PtrList<volScalarField>(phaseModels_.size());
-
-    return dmdts;
 }
 
 

applications/solvers/multiphase/reactingEulerFoam/phaseSystems/phaseSystem/phaseSystem.H

@@ -89,7 +89,7 @@ public:
 
         typedef
             HashPtrTable<volScalarField, phasePairKey, phasePairKey::hash>
-            dmdtTable;
+            dmdtfTable;
 
         typedef
             HashPtrTable
@@ -98,7 +98,7 @@ public:
                 phasePairKey,
                 phasePairKey::hash
             >
-            dmidtTable;
+            dmidtfTable;
 
 
 
@@ -415,6 +415,9 @@ public:
                 label patchi
             ) const;
 
+            //- Return the mass transfer rate for an interface
+            virtual tmp<volScalarField> dmdtf(const phasePairKey& key) const;
+
             //- Return the mass transfer rates for each phase
             virtual PtrList<volScalarField> dmdts() const;
 

applications/solvers/multiphase/reactingEulerFoam/phaseSystems/phaseSystem/phaseSystemTemplates.C

@@ -444,40 +444,6 @@ namespace Foam
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-template<class Type, class Group>
-inline tmp<GeometricField<Type, fvPatchField, volMesh>> zeroVolField
-(
-    const Group& group,
-    const word& name,
-    const dimensionSet& dims
-)
-{
-    return GeometricField<Type, fvPatchField, volMesh>::New
-    (
-        IOobject::groupName(name, group.name()),
-        group.mesh(),
-        dimensioned<Type>(dims, pTraits<Type>::zero)
-    );
-}
-
-
-template<class Type, class Group>
-inline tmp<GeometricField<Type, fvsPatchField, surfaceMesh>> zeroSurfaceField
-(
-    const Group& group,
-    const word& name,
-    const dimensionSet& dims
-)
-{
-    return GeometricField<Type, fvsPatchField, surfaceMesh>::New
-    (
-        IOobject::groupName(name, group.name()),
-        group.mesh(),
-        dimensioned<Type>(dims, pTraits<Type>::zero)
-    );
-}
-
-
 template<class GeoField, class Group>
 inline void addField
 (

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

@@ -55,8 +55,8 @@ Foam::diameterModels::driftModels::phaseChange::phaseChange
     pairKeys_(dict.lookup("pairs")),
     numberWeighted_(dict.lookupOrDefault<Switch>("numberWeighted", false)),
     W_(pairKeys_.size()),
-    specieName_(dict.lookupOrDefault("specie", word())),
+    dmdtfName_(dict.lookup("dmdtf")),
-    dmdtName_(specieName_ == word::null ? "iDmdt" : specieName_+":dmdt")
+    specieName_(dict.lookupOrDefault("specie", word()))
 {
     const phaseSystem& fluid = popBal_.fluid();
 
@@ -144,20 +144,28 @@ void Foam::diameterModels::driftModels::phaseChange::addToDriftRate
 
         if (pair.contains(vg.phase()))
         {
-            const volScalarField& iDmdt =
+            const volScalarField& dmidtf =
                 popBal_.mesh().lookupObject<volScalarField>
                 (
-                    IOobject::groupName(dmdtName_, pair.name())
+                    IOobject::groupName
+                    (
+                        IOobject::groupName
+                        (
+                            dmdtfName_,
+                            specieName_
+                        ),
+                        pair.name()
+                    )
                 );
 
-            const scalar dmdtSign =
+            const scalar dmidtfSign =
                 vg.phase().name() == pair.first() ? +1 : -1;
 
             const sizeGroup& fi = popBal_.sizeGroups()[i];
 
             tmp<volScalarField> dDriftRate
             (
-                dmdtSign*iDmdt/(fi.phase().rho()*W_[k])
+                dmidtfSign*dmidtf/(fi.phase().rho()*W_[k])
             );
 
             if (!numberWeighted_)

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

@@ -64,17 +64,17 @@ class phaseChange
 
         //- Distribute phase change mass flux between sizeGroups based on the
         //  number concentration, rather than the interfacial area
-        Switch numberWeighted_;
+        const Switch numberWeighted_;
 
         //- Weighting with which the phase change mass flux is distributed
         PtrList<volScalarField> W_;
 
+        //- Name of the mass transfer rate
+        const word dmdtfName_;
+
         //- Name of the specie which changes phase, empty by default
         const word specieName_;
 
-        //- The name of the relevant dmdt term for lookup
-        const word dmdtName_;
-
 
 public:
 

applications/solvers/multiphase/reactingEulerFoam/phaseSystems/populationBalanceModel/nucleationModels/reactionDriven/reactionDriven.C

@@ -91,18 +91,8 @@ reactionDriven
             phasePair(velGroup_.phase(), reactingPhase_)
         ]
     ),
-    specie_(dict.lookup("specie")),
+    dmdtfName_(dict.lookup("dmdtf")),
-    nDmdt_
+    specieName_(dict.lookup("specie"))
-    (
-        IOobject
-        (
-            IOobject::groupName(type() + ":N", pair_.name()),
-            popBal.time().timeName(),
-            popBal.mesh()
-        ),
-        popBal.mesh(),
-        dimensionedScalar(dimDensity/dimTime, Zero)
-    )
 {
     if
     (
@@ -121,20 +111,6 @@ reactionDriven
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
-void Foam::diameterModels::nucleationModels::reactionDriven::correct()
-{
-    const volScalarField& speciesDmdt =
-        popBal_.mesh().lookupObject<volScalarField>
-        (
-            IOobject::groupName(specie_+":dmdt", pair_.name())
-        );
-
-    const scalar sign = velGroup_.phase().name() == pair_.first() ? +1 : -1;
-
-    nDmdt_ = sign*speciesDmdt;
-}
-
-
 void
 Foam::diameterModels::nucleationModels::reactionDriven::addToNucleationRate
 (
@@ -145,8 +121,25 @@ Foam::diameterModels::nucleationModels::reactionDriven::addToNucleationRate
     const sizeGroup& fi = popBal_.sizeGroups()[i];
     const volScalarField& rho = fi.phase().rho();
 
+    const volScalarField& dmidtf =
+        popBal_.mesh().lookupObject<volScalarField>
+        (
+            IOobject::groupName
+            (
+                IOobject::groupName
+                (
+                    dmdtfName_,
+                    specieName_
+                ),
+                pair_.name()
+            )
+        );
+
+    const scalar dmidtfSign =
+        velGroup_.phase().name() == pair_.first() ? +1 : -1;
+
     nucleationRate +=
-        popBal_.eta(i, pi/6.0*pow3(dNuc_))*nDmdt_/rho/fi.x();
+        popBal_.eta(i, pi/6.0*pow3(dNuc_))*dmidtfSign*dmidtf/rho/fi.x();
 }
 
 

applications/solvers/multiphase/reactingEulerFoam/phaseSystems/populationBalanceModel/nucleationModels/reactionDriven/reactionDriven.H

@@ -58,7 +58,7 @@ class reactionDriven
     // Private Data
 
         //- Sphere-equivalent nucleation diameter
-        dimensionedScalar dNuc_;
+        const dimensionedScalar dNuc_;
 
         //- Velocity group in which the nucleation occurs
         const velocityGroup& velGroup_;
@@ -67,13 +67,13 @@ class reactionDriven
         const phaseModel& reactingPhase_;
 
         //- PhasePair between which the nucleation occurs
-        phasePair pair_;
+        const phasePair pair_;
 
-        //- Names of the specie which nucleates
+        //- Name of the mass transfer rate
-        word specie_;
+        const word dmdtfName_;
 
-        //- Nucleation mass transfer rate
+        //- Name of the specie which nucleates
-        volScalarField nDmdt_;
+        const word specieName_;
 
 
 public:
@@ -97,9 +97,6 @@ public:
 
     // Member Functions
 
-        //- Correct diameter independent expressions
-        virtual void correct();
-
         //- Add to nucleationRate
         virtual void addToNucleationRate
         (

applications/solvers/multiphase/reactingEulerFoam/phaseSystems/populationBalanceModel/nucleationModels/wallBoiling/wallBoiling.C

@@ -169,7 +169,7 @@ Foam::diameterModels::nucleationModels::wallBoiling::addToNucleationRate
             const alphatWallBoilingWallFunction& alphatw =
                 refCast<const alphatWallBoilingWallFunction>(alphatBf[patchi]);
 
-            const scalarField& dmdt = alphatw.dmdt();
+            const scalarField& dmdt = alphatw.dmdtf();
             const scalarField& dDep = alphatw.dDeparture();
 
             const labelList& faceCells = alphatw.patch().faceCells();

applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/twoPhaseSystem/diameterModels/IATE/IATE.C

@@ -106,6 +106,7 @@ Foam::tmp<Foam::volScalarField> Foam::diameterModels::IATE::dsm() const
     return max(6/max(kappai_, 6/dMax_), dMin_);
 }
 
+
 void Foam::diameterModels::IATE::correct()
 {
     volScalarField alphaAv

applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/twoPhaseSystem/diameterModels/IATE/IATEsources/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) 2016-2018 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2016-2019 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -55,8 +55,9 @@ Foam::diameterModels::IATEsources::phaseChange::phaseChange
 )
 :
     IATEsource(iate),
-    pairName_(dict.lookup("pairName")),
+    otherPhaseName_(dict.lookup("otherPhase")),
-    iDmdtPtr_(nullptr)
+    dmdtfName_(dict.lookup("dmdtf")),
+    specieName_(dict.lookupOrDefault<word>("specie", word::null))
 {}
 
 
@@ -69,23 +70,26 @@ Foam::diameterModels::IATEsources::phaseChange::R
     volScalarField& kappai
 ) const
 {
-    if (!iDmdtPtr_)
+    const phasePair& pair =
-    {
+        phase().fluid().phasePairs()
-        iDmdtPtr_ = &alphai.mesh().lookupObject<volScalarField>
+        [
+            phasePairKey(phase().name(), otherPhaseName_)
+        ];
+
+    const volScalarField& dmdtf =
+        alphai.mesh().lookupObject<volScalarField>
         (
-            IOobject::groupName("iDmdt", pairName_)
+            IOobject::groupName
+            (
+                IOobject::groupName(dmdtfName_, specieName_),
+                pair.name()
+            )
         );
-    }
 
-    const volScalarField& iDmdt = *iDmdtPtr_;
+    const label dmdtfSign =
+        phase().name() == pair.first() ? +1 : -1;
 
-    return -fvm::SuSp
+    return -fvm::SuSp(dmdtfSign*dmdtf/(3*alphai*phase().rho()), kappai);
-    (
-        (1.0/3.0)
-       *iDmdt()
-       /(alphai()*phase().rho()()),
-        kappai
-    );
 }
 
 

applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/twoPhaseSystem/diameterModels/IATE/IATEsources/phaseChange/phaseChange.H

@@ -56,11 +56,15 @@ class phaseChange
 {
     // Private Data
 
-        //- Phase pair name
+        //- Name of the other phase
-        word pairName_;
+        const word otherPhaseName_;
+
+        //- Name of the interfacial mass transfer rate
+        const word dmdtfName_;
+
+        //- Name of the specie which is transferring (optional)
+        const word specieName_;
 
-        //- Pointer to the corresponding iDmdt field
-        mutable const volScalarField* iDmdtPtr_;
 
 public:
 

applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/twoPhaseSystem/diameterModels/IATE/IATEsources/wallBoiling/wallBoiling.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2016-2018 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2016-2019 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -117,7 +117,7 @@ Foam::diameterModels::IATEsources::wallBoiling::R
             const alphatWallBoilingWallFunction& alphatw =
                 refCast<const alphatWallBoilingWallFunction>(alphatBf[patchi]);
 
-            const scalarField& dmdt = alphatw.dmdt();
+            const scalarField& dmdt = alphatw.dmdtf();
             const scalarField& dDep = alphatw.dDeparture();
 
             const labelList& faceCells = alphatw.patch().faceCells();

tutorials/multiphase/reactingMultiphaseEulerFoam/RAS/wallBoilingPolydisperseTwoGroups/constant/phaseProperties

@@ -135,10 +135,13 @@ populationBalanceCoeffs
         (
             phaseChange
             {
-                pairs ( (gas and liquid) (gas2 and liquid) );
+                pairs ((gas and liquid) (gas2 and liquid));
+                dmdtf thermalPhaseChange:dmdtf;
             }
 
-            densityChange{}
+            densityChange
+            {
+            }
         );
 
         nucleationModels

tutorials/multiphase/reactingMultiphaseEulerFoam/RAS/wallBoilingPolydisperseTwoGroups/system/fvSolution

@@ -108,7 +108,7 @@ relaxationFactors
 {
     fields
     {
-        iDmdt           0.8;
+        thermalPhaseChange:dmdtf 0.8;
     }
 
     equations

tutorials/multiphase/reactingTwoPhaseEulerFoam/RAS/bubblePipe/0/alphat.gas

@@ -33,7 +33,7 @@ boundaryField
     }
     wall
     {
-        type            compressible::alphatPhaseChangeJayatillekeWallFunction;
+        type            compressible::alphatPhaseJayatillekeWallFunction;
         Prt             0.85;
         Cmu             0.09;
         kappa           0.41;

tutorials/multiphase/reactingTwoPhaseEulerFoam/RAS/bubblePipe/0/alphat.liquid

@@ -33,7 +33,7 @@ boundaryField
     }
     wall
     {
-        type            compressible::alphatPhaseChangeJayatillekeWallFunction;
+        type            compressible::alphatPhaseJayatillekeWallFunction;
         Prt             0.85;
         Cmu             0.09;
         kappa           0.41;

tutorials/multiphase/reactingTwoPhaseEulerFoam/RAS/bubblePipe/system/fvSolution

@@ -81,7 +81,7 @@ relaxationFactors
 {
     fields
     {
-        iDmdt           1.0;
+        thermalPhaseChange:dmdtf 1.0;
     }
 
     equations

tutorials/multiphase/reactingTwoPhaseEulerFoam/RAS/wallBoiling/0/alphat.liquid

@@ -39,7 +39,6 @@ boundaryField
         kappa           0.41;
         E               9.8;
         relax           0.2;
-        dmdt            uniform 0;
         otherPhase      gas;
         phaseType       liquid;
         partitioningModel

tutorials/multiphase/reactingTwoPhaseEulerFoam/RAS/wallBoiling/system/fvSolution

@@ -108,7 +108,7 @@ relaxationFactors
 {
     fields
     {
-        iDmdt           1.0;
+        thermalPhaseChange:dmdtf 1.0;
     }
 
     equations

tutorials/multiphase/reactingTwoPhaseEulerFoam/RAS/wallBoilingIATE/0/alphat.liquid

@@ -39,7 +39,6 @@ boundaryField
         kappa           0.41;
         E               9.8;
         relax           0.2;
-        dmdt            uniform 0;
         otherPhase      gas;
         phaseType       liquid;
         partitioningModel

tutorials/multiphase/reactingTwoPhaseEulerFoam/RAS/wallBoilingIATE/constant/phaseProperties

@@ -54,7 +54,8 @@ gas
 
             phaseChange
             {
-                pairName gasAndLiquid;
+                otherPhase  liquid;
+                dmdtf       thermalPhaseChange:dmdtf;
             }
 
             wallBoiling

tutorials/multiphase/reactingTwoPhaseEulerFoam/RAS/wallBoilingIATE/system/fvSolution

@@ -118,7 +118,7 @@ relaxationFactors
 {
     fields
     {
-        iDmdt           0.8;
+        thermalPhaseChange:dmdtf 1.0;
     }
 
     equations

tutorials/multiphase/reactingTwoPhaseEulerFoam/RAS/wallBoilingPolydisperse/0/alphat.liquid

@@ -39,7 +39,6 @@ boundaryField
         kappa           0.41;
         E               9.8;
         relax           0.2;
-        dmdt            uniform 0;
         otherPhase      gas;
         phaseType       liquid;
         partitioningModel

tutorials/multiphase/reactingTwoPhaseEulerFoam/RAS/wallBoilingPolydisperse/constant/phaseProperties

@@ -111,10 +111,13 @@ populationBalanceCoeffs
         (
             phaseChange
             {
-                pairs ( (gas and liquid) );
+                pairs   ((gas and liquid));
+                dmdtf   thermalPhaseChange:dmdtf;
             }
 
-            densityChange{}
+            densityChange
+            {
+            }
         );
 
         nucleationModels

tutorials/multiphase/reactingTwoPhaseEulerFoam/RAS/wallBoilingPolydisperse/system/fvSolution

@@ -108,7 +108,7 @@ relaxationFactors
 {
     fields
     {
-        iDmdt           0.8;
+        thermalPhaseChange:dmdtf 1.0;
     }
 
     equations

tutorials/multiphase/reactingTwoPhaseEulerFoam/laminar/steamInjection/system/fvSolution

@@ -89,7 +89,7 @@ relaxationFactors
 {
     fields
     {
-        iDmdt           1;
+        thermalPhaseChange:dmdtf 1.0;
     }
 
     equations

tutorials/multiphase/reactingTwoPhaseEulerFoam/laminar/titaniaSynthesis/constant/phaseProperties

@@ -123,6 +123,7 @@ populationBalanceCoeffs
                 nucleationDiameter  4E-10;
                 velocityGroup       particles;
                 reactingPhase       vapor;
+                dmdtf               phaseTransfer:dmidtf;
                 specie              TiO2;
             }
         );

tutorials/multiphase/reactingTwoPhaseEulerFoam/laminar/titaniaSynthesisSurface/constant/phaseProperties

@@ -123,9 +123,10 @@ populationBalanceCoeffs
         (
             phaseChange
             {
-                pairs ( (particles and vapor) );
+                pairs               ((particles and vapor));
-                specie TiO2_s;
+                dmdtf               phaseTransfer:dmidtf;
-                surfaceGrowthType ParkRogak;
+                specie              TiO2_s;
+                surfaceGrowthType   ParkRogak;
             }
         );
 
@@ -136,6 +137,7 @@ populationBalanceCoeffs
                 nucleationDiameter  4E-10;
                 velocityGroup       particles;
                 reactingPhase       vapor;
+                dmdtf               phaseTransfer:dmidtf;
                 specie              TiO2;
             }
         );