Merge branch 'master' of github.com:OpenFOAM/OpenFOAM-2.3.x

applications/solvers/compressible/rhoSimpleFoam/rhoSimplecFoam/pEqn.H

@@ -17,6 +17,8 @@ if (simple.transonic())
        *(fvc::interpolate(HbyA) & mesh.Sf())
     );
 
+    fvOptions.makeRelative(fvc::interpolate(psi), phid);
+
     surfaceScalarField phic
     (
         "phic",
@@ -61,6 +63,8 @@ else
         fvc::interpolate(rho*HbyA) & mesh.Sf()
     );
 
+    fvOptions.makeRelative(fvc::interpolate(rho), phiHbyA);
+
     closedVolume = adjustPhi(phiHbyA, U, p);
 
     phiHbyA += fvc::interpolate(rho*(rAtU - rAU))*fvc::snGrad(p)*mesh.magSf();

applications/solvers/multiphase/driftFluxFoam/Allwclean

@@ -2,7 +2,7 @@
 cd ${0%/*} || exit 1    # run from this directory
 set -x
 
-wclean libso viscosityModels
+wclean libso mixtureViscosityModels
 wclean libso relativeVelocityModels
 wclean
 

applications/solvers/multiphase/driftFluxFoam/Allwmake

@@ -2,7 +2,7 @@
 cd ${0%/*} || exit 1    # run from this directory
 set -x
 
-wmake libso viscosityModels
+wmake libso mixtureViscosityModels
 wmake libso relativeVelocityModels
 wmake
 

applications/solvers/multiphase/driftFluxFoam/Make/files

@@ -1,3 +1,4 @@
+incompressibleTwoPhaseInteractingMixture/incompressibleTwoPhaseInteractingMixture.C
 driftFluxFoam.C
 
 EXE = $(FOAM_APPBIN)/driftFluxFoam

applications/solvers/multiphase/driftFluxFoam/Make/options

@@ -1,4 +1,6 @@
 EXE_INC = \
+    -IincompressibleTwoPhaseInteractingMixture \
+    -ImixtureViscosityModels/lnInclude \
     -I$(LIB_SRC)/finiteVolume/lnInclude \
     -I$(LIB_SRC)/meshTools/lnInclude \
     -I$(LIB_SRC)/sampling/lnInclude \

applications/solvers/multiphase/driftFluxFoam/UEqn.H

@@ -4,7 +4,7 @@
     (
         fvm::ddt(rho, U)
       + fvm::div(rhoPhi, U)
-      + fvc::div(uRelModel.tau(), "div(phiUkm,Ukm)")
+      + fvc::div(UdmModel.tauDm())
       - fvm::laplacian(muEff, U)
       - fvc::div(muEff*dev2(T(fvc::grad(U))))
      ==

applications/solvers/multiphase/driftFluxFoam/alphaEqnSubCycle.H

@@ -13,7 +13,7 @@
 
     surfaceScalarField phir
     (
-        mesh.Sf() & fvc::interpolate(uRelModel.Udm())
+        mesh.Sf() & fvc::interpolate(UdmModel.Udm())
     );
 
     if (nAlphaSubCycles > 1)

applications/solvers/multiphase/driftFluxFoam/createFields.H

@@ -33,13 +33,13 @@
     // ~~~~~~~~~
 
     Info<< "Reading transportProperties\n" << endl;
-    incompressibleTwoPhaseMixture twoPhaseProperties(U, phi);
+    incompressibleTwoPhaseInteractingMixture twoPhaseProperties(U, phi);
 
     volScalarField& alpha1(twoPhaseProperties.alpha1());
     volScalarField& alpha2(twoPhaseProperties.alpha2());
 
-    const dimensionedScalar& rho1 = twoPhaseProperties.rho1();
+    const dimensionedScalar& rho1 = twoPhaseProperties.rhod();
-    const dimensionedScalar& rho2 = twoPhaseProperties.rho2();
+    const dimensionedScalar& rho2 = twoPhaseProperties.rhoc();
 
     IOdictionary transportProperties
     (
@@ -86,7 +86,7 @@
     // Relative Velocity
     // ~~~~~~~~~~~~~~~~~
 
-    autoPtr<relativeVelocityModel> uRelModelPtr
+    autoPtr<relativeVelocityModel> UdmModelPtr
     (
         relativeVelocityModel::New
         (
@@ -95,7 +95,7 @@
         )
     );
 
-    relativeVelocityModel& uRelModel(uRelModelPtr());
+    relativeVelocityModel& UdmModel(UdmModelPtr());
 
 
     // Turbulence

applications/solvers/multiphase/driftFluxFoam/driftFluxFoam.C

@@ -36,7 +36,7 @@ Description
 #include "fvCFD.H"
 #include "CMULES.H"
 #include "subCycle.H"
-#include "incompressibleTwoPhaseMixture.H"
+#include "incompressibleTwoPhaseInteractingMixture.H"
 #include "relativeVelocityModel.H"
 #include "nearWallDist.H"
 #include "wallFvPatch.H"
@@ -82,11 +82,12 @@ int main(int argc, char *argv[])
         {
             #include "alphaControls.H"
 
-            uRelModel.update();
+            UdmModel.correct();
 
             #include "alphaEqnSubCycle.H"
 
             twoPhaseProperties.correct();
+            Info<< average(twoPhaseProperties.mu()) << endl;
 
             #include "UEqn.H"
 

applications/solvers/multiphase/driftFluxFoam/incompressibleTwoPhaseInteractingMixture/incompressibleTwoPhaseInteractingMixture.C

@@ -0,0 +1,135 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2014 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 "incompressibleTwoPhaseInteractingMixture.H"
+#include "addToRunTimeSelectionTable.H"
+#include "surfaceFields.H"
+#include "fvc.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+    defineTypeNameAndDebug(incompressibleTwoPhaseInteractingMixture, 0);
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::incompressibleTwoPhaseInteractingMixture::
+incompressibleTwoPhaseInteractingMixture
+(
+    const volVectorField& U,
+    const surfaceScalarField& phi
+)
+:
+    IOdictionary
+    (
+        IOobject
+        (
+            "transportProperties",
+            U.time().constant(),
+            U.db(),
+            IOobject::MUST_READ_IF_MODIFIED,
+            IOobject::NO_WRITE
+        )
+    ),
+    twoPhaseMixture(U.mesh(), *this),
+
+    muModel_
+    (
+        mixtureViscosityModel::New
+        (
+            "mu",
+            subDict(phase1Name_),
+            U,
+            phi
+        )
+    ),
+
+    nucModel_
+    (
+        viscosityModel::New
+        (
+            "nuc",
+            subDict(phase2Name_),
+            U,
+            phi
+        )
+    ),
+
+    rhod_("rho", dimDensity, muModel_->viscosityProperties().lookup("rho")),
+    rhoc_("rho", dimDensity, nucModel_->viscosityProperties().lookup("rho")),
+
+    U_(U),
+    phi_(phi),
+
+    mu_
+    (
+        IOobject
+        (
+            "mu",
+            U_.time().timeName(),
+            U_.db()
+        ),
+        U_.mesh(),
+        dimensionedScalar("mu", dimensionSet(1, -1, -1, 0, 0), 0),
+        calculatedFvPatchScalarField::typeName
+    )
+{
+    correct();
+}
+
+
+// * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
+
+bool Foam::incompressibleTwoPhaseInteractingMixture::read()
+{
+    if (regIOobject::read())
+    {
+        if
+        (
+            muModel_().read(subDict(phase1Name_))
+         && nucModel_().read(subDict(phase2Name_))
+        )
+        {
+            muModel_->viscosityProperties().lookup("rho") >> rhod_;
+            nucModel_->viscosityProperties().lookup("rho") >> rhoc_;
+
+            return true;
+        }
+        else
+        {
+            return false;
+        }
+    }
+    else
+    {
+        return false;
+    }
+}
+
+
+// ************************************************************************* //

applications/solvers/multiphase/driftFluxFoam/incompressibleTwoPhaseInteractingMixture/incompressibleTwoPhaseInteractingMixture.H

@@ -0,0 +1,172 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2014 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/>.
+
+Class
+    Foam::incompressibleTwoPhaseInteractingMixture
+
+Description
+    A two-phase incompressible transportModel for interacting phases
+    requiring the direct evaluation of the mixture viscosity,
+    e.g. activated sludge or slurry.
+
+SourceFiles
+    incompressibleTwoPhaseInteractingMixture.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef incompressibleTwoPhaseInteractingMixture_H
+#define incompressibleTwoPhaseInteractingMixture_H
+
+#include "incompressible/transportModel/transportModel.H"
+#include "incompressible/viscosityModels/viscosityModel/viscosityModel.H"
+#include "mixtureViscosityModel.H"
+#include "twoPhaseMixture.H"
+#include "IOdictionary.H"
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+         Class incompressibleTwoPhaseInteractingMixture Declaration
+\*---------------------------------------------------------------------------*/
+
+class incompressibleTwoPhaseInteractingMixture
+:
+    public IOdictionary,
+    public transportModel,
+    public twoPhaseMixture
+{
+protected:
+
+    // Protected data
+
+        autoPtr<mixtureViscosityModel> muModel_;
+        autoPtr<viscosityModel> nucModel_;
+
+        dimensionedScalar rhod_;
+        dimensionedScalar rhoc_;
+
+        const volVectorField& U_;
+        const surfaceScalarField& phi_;
+
+        volScalarField mu_;
+
+
+public:
+
+    TypeName("incompressibleTwoPhaseInteractingMixture");
+
+
+    // Constructors
+
+        //- Construct from components
+        incompressibleTwoPhaseInteractingMixture
+        (
+            const volVectorField& U,
+            const surfaceScalarField& phi
+        );
+
+
+    //- Destructor
+    virtual ~incompressibleTwoPhaseInteractingMixture()
+    {}
+
+
+    // Member Functions
+
+        //- Return const-access to the mixture viscosityModel
+        const mixtureViscosityModel& muModel() const
+        {
+            return muModel_();
+        }
+
+        //- Return const-access to the continuous-phase viscosityModel
+        const viscosityModel& nucModel() const
+        {
+            return nucModel_();
+        }
+
+        //- Return const-access to the dispersed-phase density
+        const dimensionedScalar& rhod() const
+        {
+            return rhod_;
+        }
+
+        //- Return const-access to continuous-phase density
+        const dimensionedScalar& rhoc() const
+        {
+            return rhoc_;
+        };
+
+        //- Return const-access to the mixture velocity
+        const volVectorField& U() const
+        {
+            return U_;
+        }
+
+        //- Return the dynamic mixture viscosity
+        tmp<volScalarField> mu() const
+        {
+            return mu_;
+        }
+
+        //- Return the mixture viscosity
+        virtual tmp<volScalarField> nu() const
+        {
+            notImplemented("incompressibleTwoPhaseInteractingMixture::nu()");
+            return volScalarField::null();
+        }
+
+        //- Return the mixture viscosity for patch
+        virtual tmp<scalarField> nu(const label patchi) const
+        {
+            notImplemented
+            (
+                "incompressibleTwoPhaseInteractingMixture::nu(const label)"
+            );
+            return scalarField::null();
+        }
+
+        //- Correct the laminar viscosity
+        virtual void correct()
+        {
+            mu_ = muModel_->mu(rhoc_*nucModel_->nu());
+        }
+
+        //- Read base transportProperties dictionary
+        virtual bool read();
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

applications/solvers/multiphase/driftFluxFoam/kEpsilon.H

@@ -21,7 +21,7 @@ if (turbulence)
         Cmu*k/sigmak*(g & fvc::grad(rho))/(epsilon + epsilonMin)
     );
 
-    volScalarField muc(twoPhaseProperties.nuModel2().nu()*rho2);
+    volScalarField muc(twoPhaseProperties.nucModel().nu()*rho2);
 
     #include "wallFunctions.H"
 

applications/solvers/multiphase/driftFluxFoam/mixtureViscosityModels/BinghamPlastic/BinghamPlastic.C

@@ -25,20 +25,19 @@ License
 
 #include "BinghamPlastic.H"
 #include "addToRunTimeSelectionTable.H"
-#include "surfaceFields.H"
+#include "fvcGrad.H"
-#include "fvc.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
 namespace Foam
 {
-namespace viscosityModels
+namespace mixtureViscosityModels
 {
     defineTypeNameAndDebug(BinghamPlastic, 0);
 
     addToRunTimeSelectionTable
     (
-        viscosityModel,
+        mixtureViscosityModel,
         BinghamPlastic,
         dictionary
     );
@@ -46,56 +45,9 @@ namespace viscosityModels
 }
 
 
-// * * * * * * * * * * * * Protected Member Functions  * * * * * * * * * * * //
-
-Foam::tmp<Foam::volScalarField>
-Foam::viscosityModels::BinghamPlastic::correctionNu
-(
-    const dimensionedScalar& rhoc,
-    const dimensionedScalar& rhop,
-    const volScalarField& nuc
-) const
-{
-    volScalarField
-        tauy
-        (
-            yieldStressCoeff_
-           *(
-                pow
-                (
-                    scalar(10),
-                    yieldStressExponent_
-                   *(max(alpha_, scalar(0)) + yieldStressOffset_)
-                )
-              - pow
-                (
-                    scalar(10),
-                    yieldStressExponent_*yieldStressOffset_
-                )
-            )
-        );
-
-    volScalarField
-        nup
-        (
-            plastic::correctionNu(rhoc, rhop, nuc)
-        );
-
-    dimensionedScalar tauySmall("tauySmall", tauy.dimensions(), SMALL);
-
-    return
-        tauy
-       /(
-            mag(fvc::grad(U_))
-          + 1.0e-4*(tauy + tauySmall)/(nup + (rhoc/rhop)*nuc)
-        )
-      + nup;
-}
-
-
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
-Foam::viscosityModels::BinghamPlastic::BinghamPlastic
+Foam::mixtureViscosityModels::BinghamPlastic::BinghamPlastic
 (
     const word& name,
     const dictionary& viscosityProperties,
@@ -104,16 +56,72 @@ Foam::viscosityModels::BinghamPlastic::BinghamPlastic
 )
 :
     plastic(name, viscosityProperties, U, phi, typeName),
-    yieldStressCoeff_(plasticCoeffs_.lookup("yieldStressCoeff")),
+    yieldStressCoeff_
-    yieldStressExponent_(plasticCoeffs_.lookup("yieldStressExponent")),
+    (
-    yieldStressOffset_(plasticCoeffs_.lookup("yieldStressOffset")),
+        "BinghamCoeff",
+        dimensionSet(1, -1, -2, 0, 0),
+        plasticCoeffs_.lookup("BinghamCoeff")
+    ),
+    yieldStressExponent_
+    (
+        "BinghamExponent",
+        dimless,
+        plasticCoeffs_.lookup("BinghamExponent")
+    ),
+    yieldStressOffset_
+    (
+        "BinghamOffset",
+        dimless,
+        plasticCoeffs_.lookup("BinghamOffset")
+    ),
     U_(U)
 {}
 
 
 // * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
 
-bool Foam::viscosityModels::BinghamPlastic::read
+Foam::tmp<Foam::volScalarField>
+Foam::mixtureViscosityModels::BinghamPlastic::mu
+(
+    const volScalarField& muc
+) const
+{
+    volScalarField tauy
+    (
+        yieldStressCoeff_
+       *(
+            pow
+            (
+                scalar(10),
+                yieldStressExponent_
+               *(max(alpha_, scalar(0)) + yieldStressOffset_)
+            )
+          - pow
+            (
+                scalar(10),
+                yieldStressExponent_*yieldStressOffset_
+            )
+        )
+    );
+
+    volScalarField mup(plastic::mu(muc));
+
+    dimensionedScalar tauySmall("tauySmall", tauy.dimensions(), SMALL);
+
+    return min
+    (
+        tauy
+       /(
+            mag(fvc::grad(U_))
+          + 1.0e-4*(tauy + tauySmall)/mup
+        )
+      + mup,
+        muMax_
+    );
+}
+
+
+bool Foam::mixtureViscosityModels::BinghamPlastic::read
 (
     const dictionary& viscosityProperties
 )

applications/solvers/multiphase/driftFluxFoam/mixtureViscosityModels/BinghamPlastic/BinghamPlastic.H

@@ -22,7 +22,7 @@ License
     along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
 
 Class
-    Foam::viscosityModels::BinghamPlastic
+    Foam::mixtureViscosityModels::BinghamPlastic
 
 Description
      Viscosity correction model for Bingham plastics.
@@ -41,7 +41,7 @@ SourceFiles
 
 namespace Foam
 {
-namespace viscosityModels
+namespace mixtureViscosityModels
 {
 
 /*---------------------------------------------------------------------------*\
@@ -69,17 +69,6 @@ protected:
         const volVectorField& U_;
 
 
-    // Protected Member Functions
-
-        //- Calculate and return the laminar viscosity correction
-        virtual tmp<volScalarField> correctionNu
-        (
-            const dimensionedScalar& rhoc,
-            const dimensionedScalar& rhop,
-            const volScalarField& nuc
-        ) const;
-
-
 public:
 
     //- Runtime type information
@@ -105,6 +94,10 @@ public:
 
     // Member Functions
 
+        //- Return the mixture viscosity
+        //  given the viscosity of the continuous phase
+        tmp<volScalarField> mu(const volScalarField& muc) const;
+
         //- Read transportProperties dictionary
         bool read(const dictionary& viscosityProperties);
 };
@@ -112,7 +105,7 @@ public:
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-} // End namespace viscosityModels
+} // End namespace mixtureViscosityModels
 } // End namespace Foam
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

applications/solvers/multiphase/driftFluxFoam/mixtureViscosityModels/Make/files

@@ -0,0 +1,7 @@
+mixtureViscosityModel/mixtureViscosityModel.C
+mixtureViscosityModel/mixtureViscosityModelNew.C
+plastic/plastic.C
+BinghamPlastic/BinghamPlastic.C
+slurry/slurry.C
+
+LIB = $(FOAM_LIBBIN)/libdriftFluxTransportModels

applications/solvers/multiphase/driftFluxFoam/mixtureViscosityModels/Make/options

@@ -1,4 +1,5 @@
 EXE_INC = \
+    -I../incompressibleTwoPhaseInteractingMixture \
     -I$(LIB_SRC)/finiteVolume/lnInclude \
     -I$(LIB_SRC)/transportModels \
     -I$(LIB_SRC)/transportModels/twoPhaseMixture/lnInclude \

applications/solvers/multiphase/driftFluxFoam/mixtureViscosityModels/mixtureViscosityModel/mixtureViscosityModel.C

@@ -0,0 +1,65 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2014 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 "mixtureViscosityModel.H"
+#include "volFields.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+    defineTypeNameAndDebug(mixtureViscosityModel, 0);
+    defineRunTimeSelectionTable(mixtureViscosityModel, dictionary);
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::mixtureViscosityModel::mixtureViscosityModel
+(
+    const word& name,
+    const dictionary& viscosityProperties,
+    const volVectorField& U,
+    const surfaceScalarField& phi
+)
+:
+    name_(name),
+    viscosityProperties_(viscosityProperties),
+    U_(U),
+    phi_(phi)
+{}
+
+
+// * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
+
+bool Foam::mixtureViscosityModel::read(const dictionary& viscosityProperties)
+{
+    viscosityProperties_ = viscosityProperties;
+
+    return true;
+}
+
+
+// ************************************************************************* //

applications/solvers/multiphase/driftFluxFoam/mixtureViscosityModels/mixtureViscosityModel/mixtureViscosityModel.H

@@ -0,0 +1,165 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2014 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/>.
+
+Namespace
+    Foam::mixtureViscosityModels
+
+Description
+    A namespace for incompressible mixtureViscosityModel implementations.
+
+Class
+    Foam::mixtureViscosityModel
+
+Description
+    An abstract base class for incompressible mixtureViscosityModels.
+
+    The strain rate is defined by:
+
+        mag(symm(grad(U)))
+
+
+SourceFiles
+    mixtureViscosityModel.C
+    mixtureViscosityModelNew.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef mixtureViscosityModel_H
+#define mixtureViscosityModel_H
+
+#include "dictionary.H"
+#include "volFieldsFwd.H"
+#include "surfaceFieldsFwd.H"
+#include "dimensionedScalar.H"
+#include "runTimeSelectionTables.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                Class mixtureViscosityModel Declaration
+\*---------------------------------------------------------------------------*/
+
+class mixtureViscosityModel
+{
+
+protected:
+
+    // Protected data
+
+        word name_;
+        dictionary viscosityProperties_;
+
+        const volVectorField& U_;
+        const surfaceScalarField& phi_;
+
+
+    // Private Member Functions
+
+        //- Disallow copy construct
+        mixtureViscosityModel(const mixtureViscosityModel&);
+
+        //- Disallow default bitwise assignment
+        void operator=(const mixtureViscosityModel&);
+
+
+public:
+
+    //- Runtime type information
+    TypeName("mixtureViscosityModel");
+
+
+    // Declare run-time constructor selection table
+
+        declareRunTimeSelectionTable
+        (
+            autoPtr,
+            mixtureViscosityModel,
+            dictionary,
+            (
+                const word& name,
+                const dictionary& viscosityProperties,
+                const volVectorField& U,
+                const surfaceScalarField& phi
+            ),
+            (name, viscosityProperties, U, phi)
+        );
+
+
+    // Selectors
+
+        //- Return a reference to the selected viscosity model
+        static autoPtr<mixtureViscosityModel> New
+        (
+            const word& name,
+            const dictionary& viscosityProperties,
+            const volVectorField& U,
+            const surfaceScalarField& phi
+        );
+
+
+    // Constructors
+
+        //- Construct from components
+        mixtureViscosityModel
+        (
+            const word& name,
+            const dictionary& viscosityProperties,
+            const volVectorField& U,
+            const surfaceScalarField& phi
+        );
+
+
+    //- Destructor
+    virtual ~mixtureViscosityModel()
+    {}
+
+
+    // Member Functions
+
+        //- Return the phase transport properties dictionary
+        const dictionary& viscosityProperties() const
+        {
+            return viscosityProperties_;
+        }
+
+        //- Return the mixture viscosity
+        //  given the viscosity of the continuous phase
+        virtual tmp<volScalarField> mu(const volScalarField& muc) const = 0;
+
+        //- Read transportProperties dictionary
+        virtual bool read(const dictionary& viscosityProperties) = 0;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

applications/solvers/multiphase/driftFluxFoam/mixtureViscosityModels/mixtureViscosityModel/mixtureViscosityModelNew.C

@@ -0,0 +1,65 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2014 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 "mixtureViscosityModel.H"
+#include "volFields.H"
+#include "surfaceFields.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+Foam::autoPtr<Foam::mixtureViscosityModel> Foam::mixtureViscosityModel::New
+(
+    const word& name,
+    const dictionary& viscosityProperties,
+    const volVectorField& U,
+    const surfaceScalarField& phi
+)
+{
+    const word modelType(viscosityProperties.lookup("transportModel"));
+
+    Info<< "Selecting incompressible transport model " << modelType << endl;
+
+    dictionaryConstructorTable::iterator cstrIter =
+        dictionaryConstructorTablePtr_->find(modelType);
+
+    if (cstrIter == dictionaryConstructorTablePtr_->end())
+    {
+        FatalErrorIn
+        (
+            "mixtureViscosityModel::New(const volVectorField&, "
+            "const surfaceScalarField&)"
+        )   << "Unknown mixtureViscosityModel type "
+            << modelType << nl << nl
+            << "Valid mixtureViscosityModels are : " << endl
+            << dictionaryConstructorTablePtr_->sortedToc()
+            << exit(FatalError);
+    }
+
+    return autoPtr<mixtureViscosityModel>
+        (cstrIter()(name, viscosityProperties, U, phi));
+}
+
+
+// ************************************************************************* //

applications/solvers/multiphase/driftFluxFoam/mixtureViscosityModels/plastic/plastic.C

@@ -25,20 +25,18 @@ License
 
 #include "plastic.H"
 #include "addToRunTimeSelectionTable.H"
-#include "surfaceFields.H"
-#include "incompressibleTwoPhaseMixture.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
 namespace Foam
 {
-namespace viscosityModels
+namespace mixtureViscosityModels
 {
     defineTypeNameAndDebug(plastic, 0);
 
     addToRunTimeSelectionTable
     (
-        viscosityModel,
+        mixtureViscosityModel,
         plastic,
         dictionary
     );
@@ -46,90 +44,9 @@ namespace viscosityModels
 }
 
 
-// * * * * * * * * * * * * Protected Member Functions  * * * * * * * * * * * //
-
-Foam::tmp<Foam::volScalarField>
-Foam::viscosityModels::plastic::calcNu() const
-{
-    const incompressibleTwoPhaseMixture& twoPhaseProperties =
-        alpha_.mesh().lookupObject<incompressibleTwoPhaseMixture>
-        (
-            "transportProperties"
-        );
-
-    bool isThisIsPhase1(&twoPhaseProperties.nuModel1() == this);
-
-    dimensionedScalar
-        rhoc
-        (
-            isThisIsPhase1
-          ? twoPhaseProperties.rho2()
-          : twoPhaseProperties.rho1()
-        );
-
-    dimensionedScalar
-        rhop
-        (
-            isThisIsPhase1
-          ? twoPhaseProperties.rho1()
-          : twoPhaseProperties.rho2()
-        );
-
-    volScalarField
-        nuc
-        (
-            (
-                isThisIsPhase1
-              ? twoPhaseProperties.nuModel2()
-              : twoPhaseProperties.nuModel1()
-            ).nu()
-        );
-
-    volScalarField
-        nup
-        (
-            correctionNu(rhoc, rhop, nuc)
-        );
-
-    return
-        max
-        (
-            nuMin_,
-            min
-            (
-                nuMax_,
-                (
-                    nup + (rhoc/rhop)*nuc*alpha_
-                )
-            )
-        )
-       /max(alpha_, SMALL);
-}
-
-
-Foam::tmp<Foam::volScalarField>
-Foam::viscosityModels::plastic::correctionNu
-(
-    const dimensionedScalar& rhoc,
-    const dimensionedScalar& rhop,
-    const volScalarField& nuc
-) const
-{
-    return
-        plasticViscosityCoeff_
-       *(
-            pow
-            (
-                scalar(10),
-                plasticViscosityExponent_*alpha_
-            ) - scalar(1)
-        );
-}
-
-
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
-Foam::viscosityModels::plastic::plastic
+Foam::mixtureViscosityModels::plastic::plastic
 (
     const word& name,
     const dictionary& viscosityProperties,
@@ -138,18 +55,26 @@ Foam::viscosityModels::plastic::plastic
     const word modelName
 )
 :
-    viscosityModel(name, viscosityProperties, U, phi),
+    mixtureViscosityModel(name, viscosityProperties, U, phi),
     plasticCoeffs_(viscosityProperties.subDict(modelName + "Coeffs")),
     plasticViscosityCoeff_
     (
-        plasticCoeffs_.lookup("plasticViscosityCoeff")
+        "coeff",
+        dimensionSet(1, -1, -1, 0, 0),
+        plasticCoeffs_.lookup("coeff")
     ),
     plasticViscosityExponent_
     (
-        plasticCoeffs_.lookup("plasticViscosityExponent")
+        "exponent",
+        dimless,
+        plasticCoeffs_.lookup("exponent")
+    ),
+    muMax_
+    (
+        "muMax",
+        dimensionSet(1, -1, -1, 0, 0),
+        plasticCoeffs_.lookup("muMax")
     ),
-    nuMin_(plasticCoeffs_.lookup("nuMin")),
-    nuMax_(plasticCoeffs_.lookup("nuMax")),
     alpha_
     (
         U.mesh().lookupObject<volScalarField>
@@ -160,38 +85,43 @@ Foam::viscosityModels::plastic::plastic
                 viscosityProperties.dictName()
             )
         )
-    ),
-    nu_
-    (
-        IOobject
-        (
-            name,
-            U_.time().timeName(),
-            U_.db(),
-            IOobject::NO_READ,
-            IOobject::AUTO_WRITE
-        ),
-        U_.mesh(),
-        dimensionedScalar("nu", dimViscosity, 0)
     )
 {}
 
 
 // * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
 
-bool Foam::viscosityModels::plastic::read
+Foam::tmp<Foam::volScalarField>
+Foam::mixtureViscosityModels::plastic::mu(const volScalarField& muc) const
+{
+    return min
+    (
+        muc
+      + plasticViscosityCoeff_
+       *(
+            pow
+            (
+                scalar(10),
+                plasticViscosityExponent_*alpha_
+            ) - scalar(1)
+        ),
+        muMax_
+    );
+}
+
+
+bool Foam::mixtureViscosityModels::plastic::read
 (
     const dictionary& viscosityProperties
 )
 {
-    viscosityModel::read(viscosityProperties);
+    mixtureViscosityModel::read(viscosityProperties);
 
     plasticCoeffs_ = viscosityProperties.subDict(typeName + "Coeffs");
 
     plasticCoeffs_.lookup("k") >> plasticViscosityCoeff_;
     plasticCoeffs_.lookup("n") >> plasticViscosityExponent_;
-    plasticCoeffs_.lookup("nuMin") >> nuMin_;
+    plasticCoeffs_.lookup("muMax") >> muMax_;
-    plasticCoeffs_.lookup("nuMax") >> nuMax_;
 
     return true;
 }

applications/solvers/multiphase/driftFluxFoam/mixtureViscosityModels/plastic/plastic.H

@@ -22,7 +22,7 @@ License
     along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
 
 Class
-    Foam::viscosityModels::plastic
+    Foam::mixtureViscosityModels::plastic
 
 Description
      Viscosity correction model for a generic power-law plastic.
@@ -35,7 +35,7 @@ SourceFiles
 #ifndef plastic_H
 #define plastic_H
 
-#include "viscosityModel.H"
+#include "mixtureViscosityModel.H"
 #include "dimensionedScalar.H"
 #include "volFields.H"
 
@@ -44,9 +44,9 @@ SourceFiles
 namespace Foam
 {
 
-class incompressibleTwoPhaseMixture;
+class incompressibleTwoPhaseInteractingMixture;
 
-namespace viscosityModels
+namespace mixtureViscosityModels
 {
 
 /*---------------------------------------------------------------------------*\
@@ -55,7 +55,7 @@ namespace viscosityModels
 
 class plastic
 :
-    public viscosityModel
+    public mixtureViscosityModel
 {
 protected:
 
@@ -70,32 +70,12 @@ protected:
         //- Plastic viscosity exponent
         dimensionedScalar plasticViscosityExponent_;
 
-        //- Minimum viscosity
-        dimensionedScalar nuMin_;
-
         //- Maximum viscosity
-        dimensionedScalar nuMax_;
+        dimensionedScalar muMax_;
 
         //- Plastic phase fraction
         const volScalarField& alpha_;
 
-        //- Viscosity
-        volScalarField nu_;
-
-
-    // Protected Member Functions
-
-        //- Calculate and return the laminar viscosity
-        virtual tmp<volScalarField> calcNu() const;
-
-        //- Calculate and return the laminar viscosity correction
-        virtual tmp<volScalarField> correctionNu
-        (
-            const dimensionedScalar& rhoc,
-            const dimensionedScalar& rhop,
-            const volScalarField& nuc
-        ) const;
-
 
 public:
 
@@ -123,23 +103,9 @@ public:
 
     // Member Functions
 
-        //- Return the laminar viscosity
+        //- Return the mixture viscosity
-        tmp<volScalarField> nu() const
+        //  given the viscosity of the continuous phase
-        {
+        tmp<volScalarField> mu(const volScalarField& muc) const;
-            return nu_;
-        }
-
-        //- Return the laminar viscosity for patch
-        tmp<scalarField> nu(const label patchi) const
-        {
-            return nu_.boundaryField()[patchi];
-        }
-
-        //- Correct the laminar viscosity
-        void correct()
-        {
-            nu_ = calcNu();
-        }
 
         //- Read transportProperties dictionary
         bool read(const dictionary& viscosityProperties);
@@ -148,7 +114,7 @@ public:
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-} // End namespace viscosityModels
+} // End namespace mixtureViscosityModels
 } // End namespace Foam
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

applications/solvers/multiphase/driftFluxFoam/mixtureViscosityModels/slurry/slurry.C

@@ -0,0 +1,94 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2014 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 "slurry.H"
+#include "addToRunTimeSelectionTable.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace mixtureViscosityModels
+{
+    defineTypeNameAndDebug(slurry, 0);
+
+    addToRunTimeSelectionTable
+    (
+        mixtureViscosityModel,
+        slurry,
+        dictionary
+    );
+}
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::mixtureViscosityModels::slurry::slurry
+(
+    const word& name,
+    const dictionary& viscosityProperties,
+    const volVectorField& U,
+    const surfaceScalarField& phi,
+    const word modelName
+)
+:
+    mixtureViscosityModel(name, viscosityProperties, U, phi),
+    alpha_
+    (
+        U.mesh().lookupObject<volScalarField>
+        (
+            IOobject::groupName
+            (
+                viscosityProperties.lookupOrDefault<word>("alpha", "alpha"),
+                viscosityProperties.dictName()
+            )
+        )
+    )
+{}
+
+
+// * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
+
+Foam::tmp<Foam::volScalarField>
+Foam::mixtureViscosityModels::slurry::mu(const volScalarField& muc) const
+{
+    return
+    (
+        muc*(1.0 + 2.5*alpha_ + 10.05*sqr(alpha_) + 0.00273*exp(16.6*alpha_))
+    );
+}
+
+
+bool Foam::mixtureViscosityModels::slurry::read
+(
+    const dictionary& viscosityProperties
+)
+{
+    return true;
+}
+
+
+// ************************************************************************* //

applications/solvers/multiphase/driftFluxFoam/mixtureViscosityModels/slurry/slurry.H

@@ -0,0 +1,121 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2014 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/>.
+
+Class
+    Foam::mixtureViscosityModels::slurry
+
+Description
+     Thomas' viscosity correction for slurry.
+
+     References:
+     \verbatim
+         "Transport characteristics of suspension:
+          VIII. A note on the viscosity of Newtonian suspensions
+          of uniform spherical particles".
+          D.G. Thomas,
+          J. Colloid Sci. 20 (3), 1965, p267.
+    \endverbatim
+
+SourceFiles
+    slurry.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef slurry_H
+#define slurry_H
+
+#include "mixtureViscosityModel.H"
+#include "dimensionedScalar.H"
+#include "volFields.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+class incompressibleTwoPhaseInteractingMixture;
+
+namespace mixtureViscosityModels
+{
+
+/*---------------------------------------------------------------------------*\
+                           Class slurry Declaration
+\*---------------------------------------------------------------------------*/
+
+class slurry
+:
+    public mixtureViscosityModel
+{
+protected:
+
+    // Protected data
+
+        //- Slurry phase fraction
+        const volScalarField& alpha_;
+
+
+public:
+
+    //- Runtime type information
+    TypeName("slurry");
+
+
+    // Constructors
+
+        //- Construct from components
+        slurry
+        (
+            const word& name,
+            const dictionary& viscosityProperties,
+            const volVectorField& U,
+            const surfaceScalarField& phi,
+            const word modelName=typeName
+        );
+
+
+    //- Destructor
+    ~slurry()
+    {}
+
+
+    // Member Functions
+
+        //- Return the mixture viscosity
+        //  given the viscosity of the continuous phase
+        tmp<volScalarField> mu(const volScalarField& muc) const;
+
+        //- Read transportProperties dictionary
+        bool read(const dictionary& viscosityProperties);
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace mixtureViscosityModels
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

applications/solvers/multiphase/driftFluxFoam/relativeVelocityModels/Make/options

@@ -1,4 +1,6 @@
 EXE_INC = \
+    -I../incompressibleTwoPhaseInteractingMixture \
+    -I../mixtureViscosityModels/lnInclude \
     -I$(LIB_SRC)/finiteVolume/lnInclude \
     -I$(LIB_SRC)/transportModels \
     -I$(LIB_SRC)/transportModels/twoPhaseMixture/lnInclude \

applications/solvers/multiphase/driftFluxFoam/relativeVelocityModels/general/general.C

@@ -43,13 +43,13 @@ namespace relativeVelocityModels
 Foam::relativeVelocityModels::general::general
 (
     const dictionary& dict,
-    const incompressibleTwoPhaseMixture& mixture
+    const incompressibleTwoPhaseInteractingMixture& mixture
 )
 :
     relativeVelocityModel(dict, mixture),
-    a_(dict.lookup("a")),
+    a_("a", dimless, dict.lookup("a")),
-    a1_(dict.lookup("a1")),
+    a1_("a1", dimless, dict.lookup("a1")),
-    V0_(dict.lookup("V0")),
+    V0_("V0", dimVelocity, dict.lookup("V0")),
     residualAlpha_(dict.lookup("residualAlpha"))
 {}
 
@@ -62,16 +62,16 @@ Foam::relativeVelocityModels::general::~general()
 
 // * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
 
-Foam::tmp<Foam::volVectorField>
+void Foam::relativeVelocityModels::general::correct()
-Foam::relativeVelocityModels::general::Ur() const
 {
-    return
+    Udm_ =
-        V0_
+        (rhoc_/rho())
+       *V0_
        *(
-            exp(-a_*max(alphaD_ - residualAlpha_, scalar(0)))
+            exp(-a_*max(alphad_ - residualAlpha_, scalar(0)))
-          - exp(-a1_*max(alphaD_ - residualAlpha_, scalar(0)))
+          - exp(-a1_*max(alphad_ - residualAlpha_, scalar(0)))
         )
-       /max(alphaC_, residualAlpha_);
+       /max(alphac_, residualAlpha_);
 }
 
 

applications/solvers/multiphase/driftFluxFoam/relativeVelocityModels/general/general.H

@@ -79,7 +79,7 @@ public:
         general
         (
             const dictionary& dict,
-            const incompressibleTwoPhaseMixture& mixture
+            const incompressibleTwoPhaseInteractingMixture& mixture
         );
 
 
@@ -89,8 +89,8 @@ public:
 
     // Member Functions
 
-        //- Relative velocity
+        //- Update the diffusion velocity
-        virtual tmp<volVectorField> Ur() const;
+        virtual void correct();
 };
 
 

applications/solvers/multiphase/driftFluxFoam/relativeVelocityModels/relativeVelocityModel/relativeVelocityModel.C

@@ -39,69 +39,26 @@ namespace Foam
 Foam::relativeVelocityModel::relativeVelocityModel
 (
     const dictionary& dict,
-    const incompressibleTwoPhaseMixture& mixture
+    const incompressibleTwoPhaseInteractingMixture& mixture
 )
 :
     mixture_(mixture),
-
+    alphac_(mixture.alpha2()),
-    continuousPhaseName_(dict.lookup("continuousPhase")),
+    alphad_(mixture.alpha1()),
-
+    rhoc_(mixture.rhoc()),
-    alphaC_
+    rhod_(mixture.rhod()),
-    (
-        mixture.phase1Name() == continuousPhaseName_
-      ? mixture.alpha1()
-      : mixture.alpha2()
-    ),
-
-    alphaD_
-    (
-        mixture.phase1Name() == continuousPhaseName_
-      ? mixture.alpha2()
-      : mixture.alpha1()
-    ),
-
-    rhoC_
-    (
-        mixture.phase1Name() == continuousPhaseName_
-      ? mixture.rho1()
-      : mixture.rho2()
-    ),
-
-    rhoD_
-    (
-        mixture.phase1Name() == continuousPhaseName_
-      ? mixture.rho2()
-      : mixture.rho1()
-    ),
 
     Udm_
     (
         IOobject
         (
             "Udm",
-            alphaC_.time().timeName(),
+            alphac_.time().timeName(),
-            alphaC_.mesh()
+            alphac_.mesh()
         ),
-        alphaC_.mesh(),
+        alphac_.mesh(),
         dimensionedVector("Udm", dimVelocity, vector::zero),
         mixture.U().boundaryField().types()
-    ),
-
-    tau_
-    (
-        IOobject
-        (
-            "Udm",
-            alphaC_.time().timeName(),
-            alphaC_.mesh()
-        ),
-        alphaC_.mesh(),
-        dimensionedSymmTensor
-        (
-            "Udm",
-            sqr(dimVelocity)*dimDensity,
-            symmTensor::zero
-        )
     )
 {}
 
@@ -111,7 +68,7 @@ Foam::relativeVelocityModel::relativeVelocityModel
 Foam::autoPtr<Foam::relativeVelocityModel> Foam::relativeVelocityModel::New
 (
     const dictionary& dict,
-    const incompressibleTwoPhaseMixture& mixture
+    const incompressibleTwoPhaseInteractingMixture& mixture
 )
 {
     word modelType(dict.lookup(typeName));
@@ -156,19 +113,28 @@ Foam::relativeVelocityModel::~relativeVelocityModel()
 
 // * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
 
-void Foam::relativeVelocityModel::update()
+tmp<volScalarField> Foam::relativeVelocityModel::rho() const
 {
-    tmp<volVectorField> URel(Ur());
+    return alphac_*rhoc_ + alphad_*rhod_;
+}
 
-    tmp<volScalarField> betaC(alphaC_*rhoC_);
-    tmp<volScalarField> betaD(alphaD_*rhoD_);
-    tmp<volScalarField> rhoM(betaC() + betaD());
 
-    tmp<volVectorField> Udm = URel()*betaC()/rhoM;
+tmp<volSymmTensorField> Foam::relativeVelocityModel::tauDm() const
-    tmp<volVectorField> Ucm = Udm() - URel;
+{
+    volScalarField betac(alphac_*rhoc_);
+    volScalarField betad(alphad_*rhod_);
 
-    Udm_ = Udm();
+    // Calculate the relative velocity of the continuous phase w.r.t the mean
-    tau_ = betaD*sqr(Udm) + betaC*sqr(Ucm);
+    volVectorField Ucm(betad*Udm_/betac);
+
+    return tmp<volSymmTensorField>
+    (
+        new volSymmTensorField
+        (
+            "tauDm",
+            betad*sqr(Udm_) + betac*sqr(Ucm)
+        )
+    );
 }
 
 

applications/solvers/multiphase/driftFluxFoam/relativeVelocityModels/relativeVelocityModel/relativeVelocityModel.H

@@ -36,7 +36,7 @@ SourceFiles
 
 #include "fvCFD.H"
 #include "dictionary.H"
-#include "incompressibleTwoPhaseMixture.H"
+#include "incompressibleTwoPhaseInteractingMixture.H"
 #include "runTimeSelectionTables.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@@ -63,28 +63,25 @@ protected:
     // Protected data
 
         //- Mixture properties
-        const incompressibleTwoPhaseMixture& mixture_;
+        const incompressibleTwoPhaseInteractingMixture& mixture_;
 
         //- Name of the continuous phase
         const word continuousPhaseName_;
 
         //- Continuous phase fraction
-        const volScalarField& alphaC_;
+        const volScalarField& alphac_;
 
         //- Dispersed phase fraction
-        const volScalarField& alphaD_;
+        const volScalarField& alphad_;
 
         //- Continuous density
-        const dimensionedScalar& rhoC_;
+        const dimensionedScalar& rhoc_;
 
         //- Dispersed density
-        const dimensionedScalar& rhoD_;
+        const dimensionedScalar& rhod_;
 
         //- Dispersed diffusion velocity
-        volVectorField Udm_;
+        mutable volVectorField Udm_;
-
-        //- Stress
-        volSymmTensorField tau_;
 
 
 public:
@@ -98,7 +95,8 @@ public:
         autoPtr,
         relativeVelocityModel,
         dictionary,
-        (const dictionary& dict, const incompressibleTwoPhaseMixture& mixture),
+        (const dictionary& dict,
+         const incompressibleTwoPhaseInteractingMixture& mixture),
         (dict, mixture)
     );
 
@@ -109,7 +107,7 @@ public:
         relativeVelocityModel
         (
             const dictionary& dict,
-            const incompressibleTwoPhaseMixture& mixture
+            const incompressibleTwoPhaseInteractingMixture& mixture
         );
 
 
@@ -117,7 +115,7 @@ public:
     static autoPtr<relativeVelocityModel> New
     (
         const dictionary& dict,
-        const incompressibleTwoPhaseMixture& mixture
+        const incompressibleTwoPhaseInteractingMixture& mixture
     );
 
 
@@ -127,8 +125,8 @@ public:
 
     // Member Functions
 
-        //- Calculate the relative velocity of the dispersed phase
+        //- Return the mixture mean density
-        virtual tmp<volVectorField> Ur() const = 0;
+        tmp<volScalarField> rho() const;
 
         //- Return the diffusion velocity of the dispersed phase
         const volVectorField& Udm() const
@@ -137,13 +135,10 @@ public:
         }
 
         //- Return the stress tensor due to the phase transport
-        const volSymmTensorField& tau() const
+        tmp<volSymmTensorField> tauDm() const;
-        {
-            return tau_;
-        }
 
-        //- Update the stored diffusion velocity and stress
+        //- Update the diffusion velocity
-        void update();
+        virtual void correct() = 0;
 };
 
 

applications/solvers/multiphase/driftFluxFoam/relativeVelocityModels/simple/simple.C

@@ -43,13 +43,13 @@ namespace relativeVelocityModels
 Foam::relativeVelocityModels::simple::simple
 (
     const dictionary& dict,
-    const incompressibleTwoPhaseMixture& mixture
+    const incompressibleTwoPhaseInteractingMixture& mixture
 )
 :
     relativeVelocityModel(dict, mixture),
-    a_(dict.lookup("a")),
+    a_("a", dimless, dict.lookup("a")),
-    V0_(dict.lookup("V0")),
+    V0_("V0", dimVelocity, dict.lookup("V0")),
-    residualAlpha_(dict.lookup("residualAlpha"))
+    residualAlpha_("residualAlpha", dimless, dict.lookup("residualAlpha"))
 {}
 
 
@@ -61,13 +61,9 @@ Foam::relativeVelocityModels::simple::~simple()
 
 // * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
 
-Foam::tmp<Foam::volVectorField>
+void Foam::relativeVelocityModels::simple::correct()
-Foam::relativeVelocityModels::simple::Ur() const
 {
-    return
+    Udm_ = (rhoc_/rho())*V0_*pow(scalar(10), -a_*max(alphad_, scalar(0)));
-        V0_
-       *pow(scalar(10), -a_*max(alphaD_, scalar(0)))
-       /max(alphaC_, residualAlpha_);
 }
 
 

applications/solvers/multiphase/driftFluxFoam/relativeVelocityModels/simple/simple.H

@@ -76,7 +76,7 @@ public:
         simple
         (
             const dictionary& dict,
-            const incompressibleTwoPhaseMixture& mixture
+            const incompressibleTwoPhaseInteractingMixture& mixture
         );
 
 
@@ -86,8 +86,8 @@ public:
 
     // Member Functions
 
-        //- Relative velocity
+        //- Update the diffusion velocity
-        virtual tmp<volVectorField> Ur() const;
+        virtual void correct();
 };
 
 

applications/solvers/multiphase/driftFluxFoam/viscosityModels/Make/files

@@ -1,4 +0,0 @@
-plastic/plastic.C
-BinghamPlastic/BinghamPlastic.C
-
-LIB = $(FOAM_LIBBIN)/libdriftFluxTransportModels

applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseSystem/phaseModel/phaseModel.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2014 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -63,6 +63,7 @@ Foam::phaseModel::phaseModel
     (
         phaseProperties.subDict(name_)
     ),
+    alphaMax_(phaseDict_.lookupOrDefault("alphaMax", 1.0)),
     thermo_(rhoThermo::New(fluid.mesh(), name_)),
     U_
     (

applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseSystem/phaseModel/phaseModel.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2014 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -71,6 +71,9 @@ class phaseModel
 
         dictionary phaseDict_;
 
+        //- Optional maximum phase-fraction (e.g. packing limit)
+        scalar alphaMax_;
+
         //- Thermophysical properties
         autoPtr<rhoThermo> thermo_;
 
@@ -123,6 +126,13 @@ public:
         //- Return the other phase in this two-phase system
         const phaseModel& otherPhase() const;
 
+        //- Optional maximum phase-fraction (e.g. packing limit)
+        //  Defaults to 1
+        scalar alphaMax() const
+        {
+            return alphaMax_;
+        }
+
         //- Return the Sauter-mean diameter
         tmp<volScalarField> d() const;
 

applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseSystem/twoPhaseSystem.C

@@ -248,7 +248,7 @@ Foam::twoPhaseSystem::twoPhaseSystem
             pair2In1_
         )
     );
-    
+
     wallLubrication_.set
     (
         new BlendedInterfacialModel<wallLubricationModel>
@@ -264,7 +264,7 @@ Foam::twoPhaseSystem::twoPhaseSystem
             pair2In1_
         )
     );
-    
+
     turbulentDispersion_.set
     (
         new BlendedInterfacialModel<turbulentDispersionModel>
@@ -508,7 +508,7 @@ void Foam::twoPhaseSystem::solve()
                 alphaPhic1,
                 Sp,
                 Su,
-                1,
+                phase1_.alphaMax(),
                 0
             );
 

src/OpenFOAM/db/dictionary/functionEntries/codeStream/codeStream.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2014 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -286,7 +286,10 @@ Foam::functionEntries::codeStream::getFunction
     }
 
     bool haveLib = lib;
-    reduce(haveLib, andOp<bool>());
+    if (!regIOobject::masterOnlyReading)
+    {
+        reduce(haveLib, andOp<bool>());
+    }
 
     if (!haveLib)
     {

src/finiteVolume/fields/fvPatchFields/derived/pressureInletOutletVelocity/pressureInletOutletVelocityFvPatchVectorField.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2013 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2014 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -25,10 +25,10 @@ License
 
 #include "pressureInletOutletVelocityFvPatchVectorField.H"
 #include "addToRunTimeSelectionTable.H"
-#include "fvPatchFieldMapper.H"
 #include "volFields.H"
 #include "surfaceFields.H"
 
+
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
 Foam::pressureInletOutletVelocityFvPatchVectorField::
@@ -38,14 +38,12 @@ pressureInletOutletVelocityFvPatchVectorField
     const DimensionedField<vector, volMesh>& iF
 )
 :
-    mixedFvPatchVectorField(p, iF),
+    directionMixedFvPatchVectorField(p, iF),
-    phiName_("phi"),
+    phiName_("phi")
-    rhoName_("rho"),
-    applyTangentialVelocity_(false)
 {
-    refValue() = *this;
+    refValue() = vector::zero;
     refGrad() = vector::zero;
-    valueFraction() = 0.0;
+    valueFraction() = symmTensor::zero;
 }
 
 
@@ -58,12 +56,10 @@ pressureInletOutletVelocityFvPatchVectorField
     const fvPatchFieldMapper& mapper
 )
 :
-    mixedFvPatchVectorField(ptf, p, iF, mapper),
+    directionMixedFvPatchVectorField(ptf, p, iF, mapper),
-    phiName_(ptf.phiName_),
+    phiName_(ptf.phiName_)
-    rhoName_(ptf.rhoName_),
-    applyTangentialVelocity_(ptf.applyTangentialVelocity_)
 {
-    if (applyTangentialVelocity_)
+    if (ptf.tangentialVelocity_.size())
     {
         tangentialVelocity_ = mapper(ptf.tangentialVelocity_);
     }
@@ -78,26 +74,25 @@ pressureInletOutletVelocityFvPatchVectorField
     const dictionary& dict
 )
 :
-    mixedFvPatchVectorField(p, iF),
+    directionMixedFvPatchVectorField(p, iF),
-    phiName_(dict.lookupOrDefault<word>("phi", "phi")),
+    phiName_(dict.lookupOrDefault<word>("phi", "phi"))
-    rhoName_(dict.lookupOrDefault<word>("rho", "rho")),
-    applyTangentialVelocity_(false)
 {
     fvPatchVectorField::operator=(vectorField("value", dict, p.size()));
 
     if (dict.found("tangentialVelocity"))
     {
-        applyTangentialVelocity_ = true;
-
         setTangentialVelocity
         (
             vectorField("tangentialVelocity", dict, p.size())
         );
     }
+    else
+    {
+        refValue() = vector::zero;
+    }
 
-    refValue() = *this;
     refGrad() = vector::zero;
-    valueFraction() = 0.0;
+    valueFraction() = symmTensor::zero;
 }
 
 
@@ -107,11 +102,9 @@ pressureInletOutletVelocityFvPatchVectorField
     const pressureInletOutletVelocityFvPatchVectorField& pivpvf
 )
 :
-    mixedFvPatchVectorField(pivpvf),
+    directionMixedFvPatchVectorField(pivpvf),
     phiName_(pivpvf.phiName_),
-    rhoName_(pivpvf.rhoName_),
+    tangentialVelocity_(pivpvf.tangentialVelocity_)
-    tangentialVelocity_(pivpvf.tangentialVelocity_),
-    applyTangentialVelocity_(pivpvf.applyTangentialVelocity_)
 {}
 
 
@@ -122,11 +115,9 @@ pressureInletOutletVelocityFvPatchVectorField
     const DimensionedField<vector, volMesh>& iF
 )
 :
-    mixedFvPatchVectorField(pivpvf, iF),
+    directionMixedFvPatchVectorField(pivpvf, iF),
     phiName_(pivpvf.phiName_),
-    rhoName_(pivpvf.rhoName_),
+    tangentialVelocity_(pivpvf.tangentialVelocity_)
-    tangentialVelocity_(pivpvf.tangentialVelocity_),
-    applyTangentialVelocity_(pivpvf.applyTangentialVelocity_)
 {}
 
 
@@ -135,10 +126,9 @@ pressureInletOutletVelocityFvPatchVectorField
 void Foam::pressureInletOutletVelocityFvPatchVectorField::
 setTangentialVelocity(const vectorField& tangentialVelocity)
 {
-    applyTangentialVelocity_ = true;
     tangentialVelocity_ = tangentialVelocity;
-    vectorField n(patch().nf());
+    const vectorField n(patch().nf());
-    tangentialVelocity_ -= n*(n & tangentialVelocity_);
+    refValue() = tangentialVelocity_ - n*(n & tangentialVelocity_);
 }
 
 
@@ -147,8 +137,8 @@ void Foam::pressureInletOutletVelocityFvPatchVectorField::autoMap
     const fvPatchFieldMapper& m
 )
 {
-    mixedFvPatchVectorField::autoMap(m);
+    directionMixedFvPatchVectorField::autoMap(m);
-    if (applyTangentialVelocity_)
+    if (tangentialVelocity_.size())
     {
         tangentialVelocity_.autoMap(m);
     }
@@ -161,9 +151,9 @@ void Foam::pressureInletOutletVelocityFvPatchVectorField::rmap
     const labelList& addr
 )
 {
-    mixedFvPatchVectorField::rmap(ptf, addr);
+    directionMixedFvPatchVectorField::rmap(ptf, addr);
 
-    if (applyTangentialVelocity_)
+    if (tangentialVelocity_.size())
     {
         const pressureInletOutletVelocityFvPatchVectorField& tiptf =
             refCast<const pressureInletOutletVelocityFvPatchVectorField>(ptf);
@@ -180,67 +170,25 @@ void Foam::pressureInletOutletVelocityFvPatchVectorField::updateCoeffs()
         return;
     }
 
-    const surfaceScalarField& phi =
-        db().lookupObject<surfaceScalarField>(phiName_);
-
     const fvsPatchField<scalar>& phip =
-        patch().patchField<surfaceScalarField, scalar>(phi);
+        patch().lookupPatchField<surfaceScalarField, scalar>(phiName_);
 
-    vectorField n(patch().nf());
+    valueFraction() = neg(phip)*(I - sqr(patch().nf()));
-    const Field<scalar>& magSf = patch().magSf();
 
-    if (phi.dimensions() == dimVelocity*dimArea)
+    directionMixedFvPatchVectorField::updateCoeffs();
-    {
+    directionMixedFvPatchVectorField::evaluate();
-        refValue() = n*phip/magSf;
-    }
-    else if (phi.dimensions() == dimDensity*dimVelocity*dimArea)
-    {
-        const fvPatchField<scalar>& rhop =
-            patch().lookupPatchField<volScalarField, scalar>(rhoName_);
-
-        refValue() = n*phip/(rhop*magSf);
-    }
-    else
-    {
-        FatalErrorIn
-        (
-            "pressureInletOutletVelocityFvPatchVectorField::"
-            "updateCoeffs()"
-        )   << "dimensions of phi are not correct" << nl
-            << "    on patch " << this->patch().name()
-            << " of field " << this->dimensionedInternalField().name()
-            << " in file " << this->dimensionedInternalField().objectPath()
-            << exit(FatalError);
-    }
-
-    if (applyTangentialVelocity_)
-    {
-        // Adjust the tangential velocity to conserve kinetic energy
-        // of the entrained fluid
-        // scalarField magSqrUt(magSqr(tangentialVelocity_));
-        // scalarField scale
-        // (
-        //     sqrt(max(magSqrUt - magSqr(refValue()), scalar(0))/magSqrUt)
-        // );
-        // refValue() += scale*tangentialVelocity_;
-        refValue() += tangentialVelocity_;
-    }
-
-    valueFraction() = 1.0 - pos(phip);
-
-    mixedFvPatchVectorField::updateCoeffs();
 }
 
 
 void Foam::pressureInletOutletVelocityFvPatchVectorField::write
 (
     Ostream& os
-) const
+)
+const
 {
     fvPatchVectorField::write(os);
     writeEntryIfDifferent<word>(os, "phi", "phi", phiName_);
-    writeEntryIfDifferent<word>(os, "rho", "rho", rhoName_);
+    if (tangentialVelocity_.size())
-    if (applyTangentialVelocity_)
     {
         tangentialVelocity_.writeEntry("tangentialVelocity", os);
     }
@@ -255,7 +203,9 @@ void Foam::pressureInletOutletVelocityFvPatchVectorField::operator=
     const fvPatchField<vector>& pvf
 )
 {
-    fvPatchField<vector>::operator=(pvf);
+    tmp<vectorField> normalValue = transform(valueFraction(), refValue());
+    tmp<vectorField> transformGradValue = transform(I - valueFraction(), pvf);
+    fvPatchField<vector>::operator=(normalValue + transformGradValue);
 }
 
 

src/finiteVolume/fields/fvPatchFields/derived/pressureInletOutletVelocity/pressureInletOutletVelocityFvPatchVectorField.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2013 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2014 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -25,21 +25,21 @@ Class
     Foam::pressureInletOutletVelocityFvPatchVectorField
 
 Group
-    grpInletBoundaryConditions grpOutletBoundaryConditions
+    grpInletletBoundaryConditions grpOutletBoundaryConditions
 
 Description
     This velocity inlet/outlet boundary condition is applied to pressure
     boundaries where the pressure is specified.  A zero-gradient condition is
     applied for outflow (as defined by the flux); for inflow, the velocity is
-    obtained from the patch-face normal component of the internal-cell value and
+    obtained from the patch-face normal component of the internal-cell value.
-    the tangential patch velocity can be optionally specified.
+
+    The tangential patch velocity can be optionally specified.
 
     \heading Patch usage
 
     \table
         Property     | Description             | Required    | Default value
         phi          | flux field name         | no          | phi
-        rho          | density field name      | no          | rho
         tangentialVelocity | tangential velocity field | no  |
     \endtable
 
@@ -49,7 +49,6 @@ Description
     {
         type            pressureInletOutletVelocity;
         phi             phi;
-        rho             rho;
         tangentialVelocity uniform (0 0 0);
         value           uniform 0;
     }
@@ -70,7 +69,7 @@ SourceFiles
 #define pressureInletOutletVelocityFvPatchVectorField_H
 
 #include "fvPatchFields.H"
-#include "mixedFvPatchFields.H"
+#include "directionMixedFvPatchFields.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -78,29 +77,21 @@ namespace Foam
 {
 
 /*---------------------------------------------------------------------------*\
-    Class pressureInletOutletVelocityFvPatchVectorField Declaration
+       Class pressureInletOutletVelocityFvPatchVectorField Declaration
 \*---------------------------------------------------------------------------*/
 
 class pressureInletOutletVelocityFvPatchVectorField
 :
-    public mixedFvPatchVectorField
+    public directionMixedFvPatchVectorField
 {
     // Private data
 
         //- Flux field name
         word phiName_;
 
-        //- Density field name
-        word rhoName_;
-
-
-protected:
-
         //- Optional tangential velocity component
         vectorField tangentialVelocity_;
 
-        bool applyTangentialVelocity_;
-
 
 public:
 
@@ -117,17 +108,6 @@ public:
             const DimensionedField<vector, volMesh>&
         );
 
-        //- Construct by mapping given
-        //  pressureInletOutletVelocityFvPatchVectorField
-        //  onto a new patch
-        pressureInletOutletVelocityFvPatchVectorField
-        (
-            const pressureInletOutletVelocityFvPatchVectorField&,
-            const fvPatch&,
-            const DimensionedField<vector, volMesh>&,
-            const fvPatchFieldMapper&
-        );
-
         //- Construct from patch, internal field and dictionary
         pressureInletOutletVelocityFvPatchVectorField
         (
@@ -136,6 +116,16 @@ public:
             const dictionary&
         );
 
+        //- Construct by mapping given
+        //  pressureInletOutletVelocityFvPatchVectorField onto a new patch
+        pressureInletOutletVelocityFvPatchVectorField
+        (
+            const pressureInletOutletVelocityFvPatchVectorField&,
+            const fvPatch&,
+            const DimensionedField<vector, volMesh>&,
+            const fvPatchFieldMapper&
+        );
+
         //- Construct as copy
         pressureInletOutletVelocityFvPatchVectorField
         (
@@ -147,10 +137,7 @@ public:
         {
             return tmp<fvPatchVectorField>
             (
-                new pressureInletOutletVelocityFvPatchVectorField
+                new pressureInletOutletVelocityFvPatchVectorField(*this)
-                (
-                    *this
-                )
             );
         }
 
@@ -169,11 +156,7 @@ public:
         {
             return tmp<fvPatchVectorField>
             (
-                new pressureInletOutletVelocityFvPatchVectorField
+                new pressureInletOutletVelocityFvPatchVectorField(*this, iF)
-                (
-                    *this,
-                    iF
-                )
             );
         }
 
@@ -194,18 +177,6 @@ public:
                 return phiName_;
             }
 
-            //- Return the name of rho
-            const word& rhoName() const
-            {
-                return rhoName_;
-            }
-
-            //- Return reference to the name of rho to allow adjustment
-            word& rhoName()
-            {
-                return rhoName_;
-            }
-
             //- Return the tangential velocity
             const vectorField& tangentialVelocity() const
             {

src/finiteVolume/fields/fvPatchFields/derived/rotatingPressureInletOutletVelocity/rotatingPressureInletOutletVelocityFvPatchVectorField.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2014 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -33,17 +33,17 @@ License
 void Foam::rotatingPressureInletOutletVelocityFvPatchVectorField::
 calcTangentialVelocity()
 {
-    applyTangentialVelocity_ = true;
-
     const scalar t = this->db().time().timeOutputValue();
     vector om = omega_->value(t);
 
     vector axisHat = om/mag(om);
-    tangentialVelocity_ =
+    const vectorField tangentialVelocity
-        (-om) ^ (patch().Cf() - axisHat*(axisHat & patch().Cf()));
+    (
+        (-om) ^ (patch().Cf() - axisHat*(axisHat & patch().Cf()))
+    );
 
     const vectorField n(patch().nf());
-    tangentialVelocity_ -= n*(n & tangentialVelocity_);
+    refValue() = tangentialVelocity - n*(n & tangentialVelocity);
 }
 
 
@@ -72,7 +72,9 @@ rotatingPressureInletOutletVelocityFvPatchVectorField
 :
     pressureInletOutletVelocityFvPatchVectorField(ptf, p, iF, mapper),
     omega_(ptf.omega_().clone().ptr())
-{}
+{
+    calcTangentialVelocity();
+}
 
 
 Foam::rotatingPressureInletOutletVelocityFvPatchVectorField::
@@ -85,7 +87,9 @@ rotatingPressureInletOutletVelocityFvPatchVectorField
 :
     pressureInletOutletVelocityFvPatchVectorField(p, iF, dict),
     omega_(DataEntry<vector>::New("omega", dict))
-{}
+{
+    calcTangentialVelocity();
+}
 
 
 Foam::rotatingPressureInletOutletVelocityFvPatchVectorField::
@@ -96,7 +100,9 @@ rotatingPressureInletOutletVelocityFvPatchVectorField
 :
     pressureInletOutletVelocityFvPatchVectorField(rppvf),
     omega_(rppvf.omega_().clone().ptr())
-{}
+{
+    calcTangentialVelocity();
+}
 
 
 Foam::rotatingPressureInletOutletVelocityFvPatchVectorField::
@@ -108,24 +114,13 @@ rotatingPressureInletOutletVelocityFvPatchVectorField
 :
     pressureInletOutletVelocityFvPatchVectorField(rppvf, iF),
     omega_(rppvf.omega_().clone().ptr())
-{}
+{
+    calcTangentialVelocity();
+}
 
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
-void Foam::rotatingPressureInletOutletVelocityFvPatchVectorField::updateCoeffs()
-{
-    if (updated())
-    {
-        return;
-    }
-
-    calcTangentialVelocity();
-
-    pressureInletOutletVelocityFvPatchVectorField::updateCoeffs();
-}
-
-
 void Foam::rotatingPressureInletOutletVelocityFvPatchVectorField::write
 (
     Ostream& os

src/finiteVolume/fields/fvPatchFields/derived/rotatingPressureInletOutletVelocity/rotatingPressureInletOutletVelocityFvPatchVectorField.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2014 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -177,9 +177,6 @@ public:
 
     // Member functions
 
-        //- Update the coefficients associated with the patch field
-        virtual void updateCoeffs();
-
         //- Write
         virtual void write(Ostream&) const;
 };

src/thermophysicalModels/chemistryModel/chemistryModel/basicChemistryModel/basicChemistryModel.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2014 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -146,6 +146,13 @@ public:
                     const label i
                 ) = 0;
 
+                //- Return reaction rate of the specieI in reactionI
+                virtual tmp<DimensionedField<scalar, volMesh> > calculateRR
+                (
+                    const label reactionI,
+                    const label specieI
+                ) const = 0;
+
 
             // Chemistry solution
 

src/thermophysicalModels/chemistryModel/chemistryModel/chemistryModel/chemistryModel.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2014 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -626,6 +626,88 @@ Foam::label Foam::chemistryModel<CompType, ThermoType>::nEqns() const
 }
 
 
+template<class CompType, class ThermoType>
+Foam::tmp<Foam::DimensionedField<Foam::scalar, Foam::volMesh> >
+Foam::chemistryModel<CompType, ThermoType>::calculateRR
+(
+    const label reactionI,
+    const label specieI
+) const
+{
+    scalar pf, cf, pr, cr;
+    label lRef, rRef;
+
+    const volScalarField rho
+    (
+        IOobject
+        (
+            "rho",
+            this->time().timeName(),
+            this->mesh(),
+            IOobject::NO_READ,
+            IOobject::NO_WRITE,
+            false
+        ),
+        this->thermo().rho()
+    );
+
+    tmp<DimensionedField<scalar, volMesh> > tRR
+    (
+        new DimensionedField<scalar, volMesh>
+        (
+            IOobject
+            (
+                "RR",
+                this->mesh().time().timeName(),
+                this->mesh(),
+                IOobject::NO_READ,
+                IOobject::NO_WRITE
+            ),
+            this->mesh(),
+            dimensionedScalar("zero", dimMass/dimVolume/dimTime, 0.0)
+        )
+    );
+
+    DimensionedField<scalar, volMesh>& RR = tRR();
+
+    const scalarField& T = this->thermo().T();
+    const scalarField& p = this->thermo().p();
+
+    forAll(rho, celli)
+    {
+        const scalar rhoi = rho[celli];
+        const scalar Ti = T[celli];
+        const scalar pi = p[celli];
+
+        scalarField c(nSpecie_, 0.0);
+        for (label i=0; i<nSpecie_; i++)
+        {
+            const scalar Yi = Y_[i][celli];
+            c[i] = rhoi*Yi/specieThermo_[i].W();
+        }
+
+        const scalar w = omegaI
+        (
+            reactionI,
+            c,
+            Ti,
+            pi,
+            pf,
+            cf,
+            lRef,
+            pr,
+            cr,
+            rRef
+        );
+
+        RR[celli] = w*specieThermo_[specieI].W();
+
+    }
+
+    return tRR;
+}
+
+
 template<class CompType, class ThermoType>
 void Foam::chemistryModel<CompType, ThermoType>::calculate()
 {

src/thermophysicalModels/chemistryModel/chemistryModel/chemistryModel/chemistryModel.H

@@ -198,6 +198,13 @@ public:
                 const label i
             );
 
+            //- Return reaction rate of the specieI in reactionI
+            virtual tmp<DimensionedField<scalar, volMesh> > calculateRR
+            (
+                const label reactionI,
+                const label specieI
+            ) const;
+
             //- Solve the reaction system for the given time step
             //  and return the characteristic time
             virtual scalar solve(const scalar deltaT);

src/thermophysicalModels/solidChemistryModel/basicSolidChemistryModel/basicSolidChemistryModel.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2013 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2013-2014 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -84,4 +84,27 @@ Foam::basicSolidChemistryModel::RR(const label i)
 }
 
 
+Foam::tmp<Foam::DimensionedField<Foam::scalar, Foam::volMesh> >
+Foam::basicSolidChemistryModel::calculateRR
+(
+    const label reactionI,
+    const label specieI
+) const
+{
+    notImplemented
+    (
+        "Foam::DimensionedField<Foam::scalar, Foam::volMesh>&"
+        "basicSolidChemistryModel::calculateRR(const label)"
+    );
+
+    return dynamic_cast<tmp<DimensionedField<scalar, volMesh> >&>
+    (
+        const_cast<DimensionedField<scalar, volMesh>& >
+        (
+            DimensionedField<scalar, volMesh>::null()
+        )
+    );
+}
+
+
 // ************************************************************************* //

src/thermophysicalModels/solidChemistryModel/basicSolidChemistryModel/basicSolidChemistryModel.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2013-2013 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2013-2014 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -133,6 +133,13 @@ public:
             const label i
         ) const = 0;
 
+        //- Returns the reaction rate of the specieI in reactionI
+        virtual tmp<DimensionedField<scalar, volMesh> > calculateRR
+        (
+            const label reactionI,
+            const label specieI
+        ) const;
+
         //- Return sensible enthalpy for gas i [J/Kg]
         virtual tmp<volScalarField> gasHs
         (

src/thermophysicalModels/specie/equationOfState/PengRobinsonGas/PengRobinsonGas.C

@@ -2,16 +2,16 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 1991-2007 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 2014 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
 
-    OpenFOAM is free software; you can redistribute it and/or modify it
+    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
+    under the terms of the GNU General Public License as published by
-    Free Software Foundation; either version 2 of the License, or (at your
+    the Free Software Foundation, either version 3 of the License, or
-    option) any later version.
+    (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
@@ -19,26 +19,17 @@ License
     for more details.
 
     You should have received a copy of the GNU General Public License
-    along with OpenFOAM; if not, write to the Free Software Foundation,
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
-    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
-
-Description
-    Perfect gas equation of state.
 
 \*---------------------------------------------------------------------------*/
 
 #include "PengRobinsonGas.H"
 #include "IOstreams.H"
 
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-namespace Foam
-{
-
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
 template<class Specie>
-PengRobinsonGas<Specie>::PengRobinsonGas(Istream& is)
+Foam::PengRobinsonGas<Specie>::PengRobinsonGas(Istream& is)
 :
     Specie(is),
     Tc_(readScalar(is)),
@@ -62,83 +53,20 @@ Foam::PengRobinsonGas<Specie>::PengRobinsonGas
 {}
 
 
-// * * * * * * * * * * * * * * * Friend Operators  * * * * * * * * * * * * * //
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
-
-template<class Specie>
-PengRobinsonGas<Specie> operator+
-(
-    const PengRobinsonGas<Specie>& pg1,
-    const PengRobinsonGas<Specie>& pg2
-)
-{
-    scalar nMoles = pg1.nMoles() + pg2.nMoles();
-    scalar molr1 = pg1.nMoles()/nMoles;
-    scalar molr2 = pg2.nMoles()/nMoles;
-
-    return PengRobinsonGas<Specie>
-    (
-        static_cast<const Specie&>(pg1)
-      + static_cast<const Specie&>(pg2),
-        molr1*pg1.Tc_ + molr2*pg2.Tc_,
-        molr1*pg1.Pc_ + molr2*pg2.Pc_,
-        molr1*pg1.omega_ + molr2*pg2.omega_
-    );
-}
 
 
 template<class Specie>
-PengRobinsonGas<Specie> operator-
+void Foam::PengRobinsonGas<Specie>::write(Ostream& os) const
-(
-    const PengRobinsonGas<Specie>& pg1,
-    const PengRobinsonGas<Specie>& pg2
-)
 {
-    scalar nMoles = pg1.nMoles() + pg2.nMoles();
+    Specie::write(os);
-    scalar molr1 = pg1.nMoles()/nMoles;
-    scalar molr2 = pg2.nMoles()/nMoles;
-
-    return PengRobinsonGas<Specie>
-    (
-        static_cast<const Specie&>(pg1)
-      - static_cast<const Specie&>(pg2),
-        molr1*pg1.Tc_ - molr2*pg2.Tc_,
-        molr1*pg1.Pc_ - molr2*pg2.Pc_,
-        molr1*pg1.omega_ - molr2*pg2.omega_
-    );
 }
 
 
-template<class Specie>
-PengRobinsonGas<Specie> operator*
-(
-    const scalar s,
-    const PengRobinsonGas<Specie>& pg
-)
-{
-    return PengRobinsonGas<Specie>
-    (
-        s*static_cast<const Specie&>(pg),
-        pg.Tc_,
-        pg.Pc_,
-        pg.omega_
-    );
-}
-
-
-template<class Specie>
-PengRobinsonGas<Specie> operator==
-(
-    const PengRobinsonGas<Specie>& pg1,
-    const PengRobinsonGas<Specie>& pg2
-)
-{
-    return pg2 - pg1;
-}
-
 // * * * * * * * * * * * * * * * Ostream Operator  * * * * * * * * * * * * * //
 
 template<class Specie>
-Ostream& operator<<
+Foam::Ostream& Foam::operator<<
 (
     Ostream& os,
     const PengRobinsonGas<Specie>& pg
@@ -157,8 +85,4 @@ Ostream& operator<<
 }
 
 
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-} // End namespace Foam
-
 // ************************************************************************* //

src/thermophysicalModels/specie/equationOfState/PengRobinsonGas/PengRobinsonGas.H

@@ -93,10 +93,9 @@ class PengRobinsonGas
 :
     public Specie
 {
-private:
     // Private data
 
-         //- Critical Temperature [K]
+        //- Critical Temperature [K]
         scalar Tc_;
 
         //- Critical Pressure [Pa]

src/thermophysicalModels/specie/equationOfState/PengRobinsonGas/PengRobinsonGasI.H

@@ -27,13 +27,10 @@ License
 #include "PengRobinsonGas.H"
 #include "mathematicalConstants.H"
 
-namespace Foam
-{
-
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
 template<class Specie>
-inline PengRobinsonGas<Specie>::PengRobinsonGas
+inline Foam::PengRobinsonGas<Specie>::PengRobinsonGas
 (
     const Specie& sp,
     const scalar& Tc,
@@ -51,7 +48,7 @@ inline PengRobinsonGas<Specie>::PengRobinsonGas
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
 template<class Specie>
-inline PengRobinsonGas<Specie>::PengRobinsonGas
+inline Foam::PengRobinsonGas<Specie>::PengRobinsonGas
 (
     const word& name,
     const PengRobinsonGas& pg
@@ -65,7 +62,8 @@ inline PengRobinsonGas<Specie>::PengRobinsonGas
 
 
 template<class Specie>
-inline autoPtr<PengRobinsonGas <Specie> > PengRobinsonGas<Specie>::clone() const
+inline Foam::autoPtr<Foam::PengRobinsonGas <Specie> >
+Foam::PengRobinsonGas<Specie>::clone() const
 {
     return autoPtr<PengRobinsonGas<Specie> >
     (
@@ -75,7 +73,8 @@ inline autoPtr<PengRobinsonGas <Specie> > PengRobinsonGas<Specie>::clone() const
 
 
 template<class Specie>
-inline autoPtr<PengRobinsonGas<Specie> > PengRobinsonGas<Specie>::New
+inline Foam::autoPtr<Foam::PengRobinsonGas<Specie> >
+Foam::PengRobinsonGas<Specie>::New
 (
     Istream& is
 )
@@ -100,7 +99,11 @@ Foam::PengRobinsonGas<Specie>::New
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
 template<class Specie>
-inline scalar PengRobinsonGas<Specie>::rho(scalar p, scalar T) const
+inline Foam::scalar Foam::PengRobinsonGas<Specie>::rho
+(
+    scalar p,
+    scalar T
+) const
 {
     scalar z = Z(p, T);
     return p/(z*this->R()*T);
@@ -108,7 +111,11 @@ inline scalar PengRobinsonGas<Specie>::rho(scalar p, scalar T) const
 
 
 template<class Specie>
-inline scalar PengRobinsonGas<Specie>::psi(scalar p, scalar T) const
+inline Foam::scalar Foam::PengRobinsonGas<Specie>::psi
+(
+    scalar p,
+    scalar T
+) const
 {
     scalar z = Z(p, T);
     return 1.0/(z*this->R()*T);
@@ -116,7 +123,11 @@ inline scalar PengRobinsonGas<Specie>::psi(scalar p, scalar T) const
 
 
 template<class Specie>
-inline scalar PengRobinsonGas<Specie>::Z(scalar p, scalar T) const
+inline Foam::scalar Foam::PengRobinsonGas<Specie>::Z
+(
+    scalar p,
+    scalar T
+) const
 {
     scalar a = 0.45724*sqr(this->R())*sqr(Tc_)/Pc_;
     scalar b = 0.07780*this->R()*Tc_/Pc_;
@@ -191,7 +202,7 @@ inline Foam::scalar Foam::PengRobinsonGas<Specie>::cpMcv
 // * * * * * * * * * * * * * * * Member Operators  * * * * * * * * * * * * * //
 
 template<class Specie>
-inline void PengRobinsonGas<Specie>::operator+=
+inline void Foam::PengRobinsonGas<Specie>::operator+=
 (
     const PengRobinsonGas<Specie>& pg
 )
@@ -207,8 +218,9 @@ inline void PengRobinsonGas<Specie>::operator+=
     omega_ = molr1*omega_ + molr2*pg.omega_;
 }
 
+
 template<class Specie>
-inline void PengRobinsonGas<Specie>::operator-=
+inline void Foam::PengRobinsonGas<Specie>::operator-=
 (
     const PengRobinsonGas<Specie>& pg
 )
@@ -225,15 +237,87 @@ inline void PengRobinsonGas<Specie>::operator-=
     omega_ = molr1*omega_ - molr2*pg.omega_;
 }
 
+
 template<class Specie>
-inline void PengRobinsonGas<Specie>::operator*=(const scalar s)
+inline void Foam::PengRobinsonGas<Specie>::operator*=(const scalar s)
 {
      Specie::operator*=(s);
 }
 
 
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+// * * * * * * * * * * * * * * * Friend Operators  * * * * * * * * * * * * * //
+
+
+template<class Specie>
+Foam::PengRobinsonGas<Specie> Foam::operator+
+(
+    const PengRobinsonGas<Specie>& pg1,
+    const PengRobinsonGas<Specie>& pg2
+)
+{
+    scalar nMoles = pg1.nMoles() + pg2.nMoles();
+    scalar molr1 = pg1.nMoles()/nMoles;
+    scalar molr2 = pg2.nMoles()/nMoles;
+
+    return PengRobinsonGas<Specie>
+    (
+        static_cast<const Specie&>(pg1)
+      + static_cast<const Specie&>(pg2),
+        molr1*pg1.Tc_ + molr2*pg2.Tc_,
+        molr1*pg1.Pc_ + molr2*pg2.Pc_,
+        molr1*pg1.omega_ + molr2*pg2.omega_
+    );
+}
+
+
+template<class Specie>
+Foam::PengRobinsonGas<Specie> Foam::operator-
+(
+    const PengRobinsonGas<Specie>& pg1,
+    const PengRobinsonGas<Specie>& pg2
+)
+{
+    scalar nMoles = pg1.nMoles() + pg2.nMoles();
+    scalar molr1 = pg1.nMoles()/nMoles;
+    scalar molr2 = pg2.nMoles()/nMoles;
+
+    return PengRobinsonGas<Specie>
+    (
+        static_cast<const Specie&>(pg1)
+      - static_cast<const Specie&>(pg2),
+        molr1*pg1.Tc_ - molr2*pg2.Tc_,
+        molr1*pg1.Pc_ - molr2*pg2.Pc_,
+        molr1*pg1.omega_ - molr2*pg2.omega_
+    );
+}
+
+
+template<class Specie>
+Foam::PengRobinsonGas<Specie> Foam::operator*
+(
+    const scalar s,
+    const PengRobinsonGas<Specie>& pg
+)
+{
+    return PengRobinsonGas<Specie>
+    (
+        s*static_cast<const Specie&>(pg),
+        pg.Tc_,
+        pg.Pc_,
+        pg.omega_
+    );
+}
+
+
+template<class Specie>
+Foam::PengRobinsonGas<Specie> Foam::operator==
+(
+    const PengRobinsonGas<Specie>& pg1,
+    const PengRobinsonGas<Specie>& pg2
+)
+{
+    return pg2 - pg1;
+}
 
-} // End namespace Foam
 
 // ************************************************************************* //

tutorials/multiphase/driftFluxFoam/ras/dahl/constant/transportProperties

@@ -15,7 +15,7 @@ FoamFile
 }
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-phases          (sludge water);
+phases (sludge water);
 
 sludge
 {
@@ -23,37 +23,35 @@ sludge
 
     "(plastic|BinghamPlastic)Coeffs"
     {
-        plasticViscosityCoeff plasticViscosityCoeff [ 0 2 -1 0 0 0 0 ] 1.1595e-07;
+        coeff       0.00023143;
-        plasticViscosityExponent plasticViscosityExponent [ 0 0 0 0 0 0 0 ] 179.26;
+        exponent    179.26;
 
-        yieldStressCoeff yieldStressCoeff [ 0 2 -2 0 0 0 0 ] 2.1137e-07;
+        BinghamCoeff    0.00042189;
-        yieldStressExponent yieldStressExponent [ 0 0 0 0 0 0 0 ] 1050.8;
+        BinghamExponent 1050.8;
-        yieldStressOffset yieldStressOffset [ 0 0 0 0 0 0 0 ] 0;
+        BinghamOffset   0;
 
-        nuMin       nuMin [ 0 2 -1 0 0 0 0 ] 1e-10;
+        muMax       10;
-        nuMax       nuMax [ 0 2 -1 0 0 0 0 ] 5e-3;
     }
 
-    rho         rho [ 1 -3 0 0 0 0 0 ] 1996;
+    rho         1996;
 }
 
 water
 {
     transportModel  Newtonian;
 
-    nu      nu [ 0 2 -1 0 0 0 0 ] 1.7871e-06;
+    nu              1.7871e-06;
-    rho     rho [ 1 -3 0 0 0 0 0 ] 996;
+    rho             996;
 }
 
 relativeVelocityModel simple;
 
 "(simple|general)Coeffs"
 {
-    continuousPhase water;
+    V0              (0 -0.002198 0);
-    V0              V0 [ 0 1 -1 0 0 0 0 ] ( 0 -0.002198 0 );
+    a               285.84;
-    a               a [ 0 0 0 0 0 0 0 ] 285.84;
+    a1              0.1;
-    a1              a1 [ 0 0 0 0 0 0 0 ] 0.1;
+    residualAlpha   0;
-    residualAlpha   residualAlpha [ 0 0 0 0 0 0 0 ] 0;
 }
 
 

tutorials/multiphase/driftFluxFoam/ras/dahl/system/fvSchemes

@@ -30,7 +30,7 @@ divSchemes
     default             none;
 
     div(rhoPhi,U)       Gauss linearUpwind grad(U);
-    div(phiUkm,Ukm)     Gauss linear;
+    div(tauDm)          Gauss linear;
     "div\(phi,alpha.*\)" Gauss vanLeer;
     "div\(phirb,alpha.*\)" Gauss linear;
     div(rhoPhi,k)       Gauss limitedLinear 1;

tutorials/multiphase/driftFluxFoam/ras/mixerVessel2D/constant/transportProperties

@@ -15,7 +15,7 @@ FoamFile
 }
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-phases          (sludge water);
+phases (sludge water);
 
 sludge
 {
@@ -23,37 +23,35 @@ sludge
 
     "(plastic|BinghamPlastic)Coeffs"
     {
-        plasticViscosityCoeff plasticViscosityCoeff [ 0 2 -1 0 0 0 0 ] 1.1595e-07;
+        coeff       0.00023143;
-        plasticViscosityExponent plasticViscosityExponent [ 0 0 0 0 0 0 0 ] 179.26;
+        exponent    179.26;
 
-        yieldStressCoeff yieldStressCoeff [ 0 2 -2 0 0 0 0 ] 2.1137e-07;
+        BinghamCoeff    0.00042189;
-        yieldStressExponent yieldStressExponent [ 0 0 0 0 0 0 0 ] 1050.8;
+        BinghamExponent 1050.8;
-        yieldStressOffset yieldStressOffset [ 0 0 0 0 0 0 0 ] 0;
+        BinghamOffset   0;
 
-        nuMin       nuMin [ 0 2 -1 0 0 0 0 ] 1e-10;
+        muMax       10;
-        nuMax       nuMax [ 0 2 -1 0 0 0 0 ] 5e-3;
     }
 
-    rho         rho [ 1 -3 0 0 0 0 0 ] 1996;
+    rho         1996;
 }
 
 water
 {
     transportModel  Newtonian;
 
-    nu      nu [ 0 2 -1 0 0 0 0 ] 1.7871e-06;
+    nu              1.7871e-06;
-    rho     rho [ 1 -3 0 0 0 0 0 ] 996;
+    rho             996;
 }
 
 relativeVelocityModel simple;
 
 "(simple|general)Coeffs"
 {
-    continuousPhase water;
+    V0              (0 -0.002198 0);
-    V0              V0 [ 0 1 -1 0 0 0 0 ] ( 0 -0.002198 0 );
+    a               285.84;
-    a               a [ 0 0 0 0 0 0 0 ] 285.84;
+    a1              0.1;
-    a1              a1 [ 0 0 0 0 0 0 0 ] 0.1;
+    residualAlpha   0;
-    residualAlpha   residualAlpha [ 0 0 0 0 0 0 0 ] 0;
 }
 
 

tutorials/multiphase/driftFluxFoam/ras/mixerVessel2D/system/fvSchemes

@@ -30,7 +30,7 @@ divSchemes
     default             none;
 
     div(rhoPhi,U)       Gauss linearUpwind grad(U);
-    div(phiUkm,Ukm)     Gauss linear;
+    div(tauDm)          Gauss linear;
     "div\(phi,alpha.*\)" Gauss vanLeer;
     "div\(phirb,alpha.*\)" Gauss linear;
     div(rhoPhi,k)       Gauss limitedLinear 1;

tutorials/multiphase/driftFluxFoam/ras/tank3D/constant/transportProperties

@@ -15,7 +15,7 @@ FoamFile
 }
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-phases          (sludge water);
+phases (sludge water);
 
 sludge
 {
@@ -23,38 +23,35 @@ sludge
 
     "(plastic|BinghamPlastic)Coeffs"
     {
-        plasticViscosityCoeff plasticViscosityCoeff [ 0 2 -1 0 0 0 0 ] 2.2210e-07;
+        coeff       0.00023143;
-        plasticViscosityExponent plasticViscosityExponent [ 0 0 0 0 0 0 0 ] 0.17926;
+        exponent    0.17926;
 
-        yieldStressCoeff yieldStressCoeff [ 0 2 -2 0 0 0 0 ] 5.3233e-10;
+        BinghamCoeff    5.5469e-07;
-        yieldStressExponent yieldStressExponent [ 0 0 0 0 0 0 0 ] 95.25;
+        BinghamExponent 95.25;
-        yieldStressOffset yieldStressOffset [ 0 0 0 0 0 0 0 ] 0;
+        BinghamOffset   0;
 
-        nuMin       nuMin [ 0 2 -1 0 0 0 0 ] 1e-10;
+        muMax       10;
-        nuMax       nuMax [ 0 2 -1 0 0 0 0 ] 5e-3;
     }
 
-    rho         rho [ 1 -3 0 0 0 0 0 ] 1042;
+    rho         1042;
 }
 
 water
 {
     transportModel  Newtonian;
 
-    nu      nu [ 0 2 -1 0 0 0 0 ] 1.78e-06;
+    nu              1.78e-06;
-    rho     rho [ 1 -3 0 0 0 0 0 ] 1000;
+    rho             1000;
 }
 
 relativeVelocityModel simple;
 
 "(simple|general)Coeffs"
 {
-    continuousPhase water;
+    V0              (0 -0.002198 0);
-    V0              V0 [ 0 1 -1 0 0 0 0 ] ( 0 -0.002198 0 );
+    a               8.84;
-    a               a [ 0 0 0 0 0 0 0 ] 8.84;
+    a1              0.1;
-    a1              a1 [ 0 0 0 0 0 0 0 ] 0.1;
+    residualAlpha   0;
-    residualAlpha   residualAlpha [ 0 0 0 0 0 0 0 ] 0;
 }
 
-
 // ************************************************************************* //

tutorials/multiphase/driftFluxFoam/ras/tank3D/system/fvSchemes

@@ -30,7 +30,7 @@ divSchemes
     default         none;
 
     div(rhoPhi,U)       Gauss linearUpwind grad(U);
-    div(phiUkm,Ukm)     Gauss linear;
+    div(tauDm)          Gauss linear;
     "div\(phi,alpha.*\)" Gauss vanLeer;
     "div\(phirb,alpha.*\)" Gauss linear;
     div(rhoPhi,k)       Gauss limitedLinear 1;

tutorials/multiphase/interDyMFoam/ras/DTCHull/0.org/U

@@ -44,6 +44,7 @@ boundaryField
     atmosphere
     {
         type            pressureInletOutletVelocity;
+        tangentialVelocity $internalField;
         value           uniform (0 0 0);
     }
 

tutorials/multiphase/twoPhaseEulerFoam/RAS/fluidisedBed/constant/phaseProperties

@@ -24,6 +24,8 @@ particles
     {
         d               3e-4;
     }
+
+    alphaMax 0.62;
 }
 
 air

tutorials/multiphase/twoPhaseEulerFoam/laminar/fluidisedBed/constant/phaseProperties

@@ -24,6 +24,8 @@ particles
     {
         d               3e-4;
     }
+
+    alphaMax 0.62;
 }
 
 air