multiphaseEulerFoam: Work in progress

applications/solvers/multiphase/multiphaseEulerFoam/CourantNo.H

@@ -39,6 +39,15 @@ if (mesh.nInternalFaces())
         fvc::surfaceSum(mag(phi))().internalField()
     );
 
+    forAllIter(PtrDictionary<phaseModel>, fluid.phases(), iter)
+    {
+        sumPhi = max
+        (
+            sumPhi,
+            fvc::surfaceSum(mag(iter().phi()))().internalField()
+        );
+    }
+
     CoNum = 0.5*gMax(sumPhi/mesh.V().field())*runTime.deltaTValue();
 
     meanCoNum =

applications/solvers/multiphase/multiphaseEulerFoam/CourantNos.H

@@ -1,12 +0,0 @@
-#   include "CourantNo.H"
-
-// {
-//     scalar UrCoNum = 0.5*gMax
-//     (
-//         fvc::surfaceSum(mag(phi1 - phi2))().internalField()/mesh.V().field()
-//     )*runTime.deltaTValue();
-
-//     Info<< "Max Ur Courant Number = " << UrCoNum << endl;
-
-//     CoNum = max(CoNum, UrCoNum);
-// }

applications/solvers/multiphase/multiphaseEulerFoam/createFields.H

@@ -42,7 +42,7 @@
         dimensionedScalar("phi", dimArea*dimVelocity, 0)
     );
 
-    multiphaseSystem fluid(mesh, phi);
+    multiphaseSystem fluid(U, phi);
 
     forAllIter(PtrDictionary<phaseModel>, fluid.phases(), iter)
     {
@@ -77,9 +77,8 @@
     scalar pRefValue = 0.0;
     setRefCell(p, mesh.solutionDict().subDict("PIMPLE"), pRefCell, pRefValue);
 
-    singlePhaseTransportModel laminarTransport(U, phi);
 
     autoPtr<incompressible::LESModel> sgsModel
     (
-        incompressible::LESModel::New(U, phi, laminarTransport)
+        incompressible::LESModel::New(U, phi, fluid)
     );

applications/solvers/multiphase/multiphaseEulerFoam/interfacialModels/Make/files

@@ -8,6 +8,7 @@ dragModels/Gibilaro/Gibilaro.C
 dragModels/WenYu/WenYu.C
 dragModels/SyamlalOBrien/SyamlalOBrien.C
 dragModels/blended/blended.C
+dragModels/interface/interface.C
 
 heatTransferModels/heatTransferModel/heatTransferModel.C
 heatTransferModels/heatTransferModel/newHeatTransferModel.C

applications/solvers/multiphase/multiphaseEulerFoam/interfacialModels/dragModels/interface/interface.C

@@ -0,0 +1,93 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2011 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 "interface.H"
+#include "addToRunTimeSelectionTable.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace dragModels
+{
+    defineTypeNameAndDebug(interface, 0);
+
+    addToRunTimeSelectionTable
+    (
+        dragModel,
+        interface,
+        dictionary
+    );
+}
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::dragModels::interface::interface
+(
+    const dictionary& interfaceDict,
+    const phaseModel& phase1,
+    const phaseModel& phase2
+)
+:
+    dragModel(interfaceDict, phase1, phase2)
+{}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+Foam::dragModels::interface::~interface()
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+Foam::tmp<Foam::volScalarField> Foam::dragModels::interface::K
+(
+    const volScalarField& Ur
+) const
+{
+    return tmp<volScalarField>
+    (
+        new volScalarField
+        (
+            IOobject
+            (
+                "K",
+                Ur.mesh().time().timeName(),
+                Ur.mesh(),
+                IOobject::NO_READ,
+                IOobject::NO_WRITE,
+                false
+            ),
+            Ur.mesh(),
+            dimensionedScalar("K", dimDensity/dimTime, 0)
+        )
+    );
+}
+
+
+// ************************************************************************* //

applications/solvers/multiphase/multiphaseEulerFoam/interfacialModels/dragModels/interface/interface.H

@@ -0,0 +1,92 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2011 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::dragModels::interface
+
+Description
+    Drag between phase separated by a VoF resolved interface.
+
+SourceFiles
+    interface.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef interface_H
+#define interface_H
+
+#include "dragModel.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace dragModels
+{
+
+/*---------------------------------------------------------------------------*\
+                           Class interface Declaration
+\*---------------------------------------------------------------------------*/
+
+class interface
+:
+    public dragModel
+{
+
+public:
+
+    //- Runtime type information
+    TypeName("interface");
+
+
+    // Constructors
+
+        //- Construct from components
+        interface
+        (
+            const dictionary& interfaceDict,
+            const phaseModel& phase1,
+            const phaseModel& phase2
+        );
+
+
+    //- Destructor
+    virtual ~interface();
+
+
+    // Member Functions
+
+        tmp<volScalarField> K(const volScalarField& Ur) const;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace dragModels
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

applications/solvers/multiphase/multiphaseEulerFoam/multiphaseEulerFoam.C

@@ -58,7 +58,7 @@ int main(int argc, char *argv[])
     #include "initContinuityErrs.H"
     #include "readTimeControls.H"
     #include "correctPhi.H"
-    #include "CourantNos.H"
+    #include "CourantNo.H"
     #include "setInitialDeltaT.H"
 
     // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@@ -68,7 +68,7 @@ int main(int argc, char *argv[])
     while (runTime.run())
     {
         #include "readTimeControls.H"
-        #include "CourantNos.H"
+        #include "CourantNo.H"
         #include "setDeltaT.H"
 
         runTime++;
@@ -80,6 +80,7 @@ int main(int argc, char *argv[])
             sgsModel->correct();
             fluid.solve();
             rho = fluid.rho();
+            #include "zonePhaseVolumes.H"
 
             //#include "interfacialCoeffs.H"
             //#include "TEqns.H"

applications/solvers/multiphase/multiphaseEulerFoam/multiphaseFixedFluxPressure/multiphaseFixedFluxPressureFvPatchScalarField.C

@@ -0,0 +1,181 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2011 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 "multiphaseFixedFluxPressureFvPatchScalarField.H"
+#include "fvPatchFieldMapper.H"
+#include "volFields.H"
+#include "surfaceFields.H"
+#include "addToRunTimeSelectionTable.H"
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::multiphaseFixedFluxPressureFvPatchScalarField::
+multiphaseFixedFluxPressureFvPatchScalarField
+(
+    const fvPatch& p,
+    const DimensionedField<scalar, volMesh>& iF
+)
+:
+    fixedGradientFvPatchScalarField(p, iF),
+    phi0Name_("phi0"),
+    phiName_("phi"),
+    rhoName_("rho")
+{}
+
+
+Foam::multiphaseFixedFluxPressureFvPatchScalarField::
+multiphaseFixedFluxPressureFvPatchScalarField
+(
+    const multiphaseFixedFluxPressureFvPatchScalarField& ptf,
+    const fvPatch& p,
+    const DimensionedField<scalar, volMesh>& iF,
+    const fvPatchFieldMapper& mapper
+)
+:
+    fixedGradientFvPatchScalarField(ptf, p, iF, mapper),
+    phi0Name_(ptf.phi0Name_),
+    phiName_(ptf.phiName_),
+    rhoName_(ptf.rhoName_)
+{}
+
+
+Foam::multiphaseFixedFluxPressureFvPatchScalarField::
+multiphaseFixedFluxPressureFvPatchScalarField
+(
+    const fvPatch& p,
+    const DimensionedField<scalar, volMesh>& iF,
+    const dictionary& dict
+)
+:
+    fixedGradientFvPatchScalarField(p, iF),
+    phi0Name_(dict.lookupOrDefault<word>("phi0", "phi0")),
+    phiName_(dict.lookupOrDefault<word>("phi", "phi")),
+    rhoName_(dict.lookupOrDefault<word>("rho", "rho"))
+{
+    if (dict.found("gradient"))
+    {
+        gradient() = scalarField("gradient", dict, p.size());
+        fixedGradientFvPatchScalarField::updateCoeffs();
+        fixedGradientFvPatchScalarField::evaluate();
+    }
+    else
+    {
+        fvPatchField<scalar>::operator=(patchInternalField());
+        gradient() = 0.0;
+    }
+}
+
+
+Foam::multiphaseFixedFluxPressureFvPatchScalarField::
+multiphaseFixedFluxPressureFvPatchScalarField
+(
+    const multiphaseFixedFluxPressureFvPatchScalarField& wbppsf
+)
+:
+    fixedGradientFvPatchScalarField(wbppsf),
+    phi0Name_(wbppsf.phi0Name_),
+    phiName_(wbppsf.phiName_),
+    rhoName_(wbppsf.rhoName_)
+{}
+
+
+Foam::multiphaseFixedFluxPressureFvPatchScalarField::
+multiphaseFixedFluxPressureFvPatchScalarField
+(
+    const multiphaseFixedFluxPressureFvPatchScalarField& wbppsf,
+    const DimensionedField<scalar, volMesh>& iF
+)
+:
+    fixedGradientFvPatchScalarField(wbppsf, iF),
+    phi0Name_(wbppsf.phi0Name_),
+    phiName_(wbppsf.phiName_),
+    rhoName_(wbppsf.rhoName_)
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+void Foam::multiphaseFixedFluxPressureFvPatchScalarField::updateCoeffs()
+{
+    if (updated())
+    {
+        return;
+    }
+
+    const surfaceScalarField& phi0 =
+        db().lookupObject<surfaceScalarField>(phi0Name_);
+
+    const surfaceScalarField& phi =
+        db().lookupObject<surfaceScalarField>(phiName_);
+
+    fvsPatchField<scalar> phi0p =
+        patch().patchField<surfaceScalarField, scalar>(phi0);
+
+    fvsPatchField<scalar> phip =
+        patch().patchField<surfaceScalarField, scalar>(phi);
+
+    if (phi.dimensions() == dimDensity*dimVelocity*dimArea)
+    {
+        const fvPatchField<scalar>& rhop =
+            patch().lookupPatchField<volScalarField, scalar>(rhoName_);
+
+        phip /= rhop;
+    }
+
+    const fvsPatchField<scalar>& Dpp =
+        patch().lookupPatchField<surfaceScalarField, scalar>("Dp");
+
+    gradient() = (phi0p - phip)/patch().magSf()/Dpp;
+
+    fixedGradientFvPatchScalarField::updateCoeffs();
+}
+
+
+void Foam::multiphaseFixedFluxPressureFvPatchScalarField::write
+(
+    Ostream& os
+) const
+{
+    fvPatchScalarField::write(os);
+    writeEntryIfDifferent<word>(os, "phi0", "phi0", phi0Name_);
+    writeEntryIfDifferent<word>(os, "phi", "phi", phiName_);
+    writeEntryIfDifferent<word>(os, "rho", "rho", rhoName_);
+    gradient().writeEntry("gradient", os);
+    writeEntry("value", os);
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+    makePatchTypeField
+    (
+        fvPatchScalarField,
+        multiphaseFixedFluxPressureFvPatchScalarField
+    );
+}
+
+// ************************************************************************* //

applications/solvers/multiphase/multiphaseEulerFoam/multiphaseFixedFluxPressure/multiphaseFixedFluxPressureFvPatchScalarField.H

@@ -0,0 +1,154 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2011 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::multiphaseFixedFluxPressureFvPatchScalarField
+
+Description
+    Foam::multiphaseFixedFluxPressureFvPatchScalarField
+
+SourceFiles
+    multiphaseFixedFluxPressureFvPatchScalarField.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef multiphaseFixedFluxPressureFvPatchScalarFields_H
+#define multiphaseFixedFluxPressureFvPatchScalarFields_H
+
+#include "fvPatchFields.H"
+#include "fixedGradientFvPatchFields.H"
+#include "Switch.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+              Class multiphaseFixedFluxPressureFvPatch Declaration
+\*---------------------------------------------------------------------------*/
+
+class multiphaseFixedFluxPressureFvPatchScalarField
+:
+    public fixedGradientFvPatchScalarField
+{
+    // Private data
+
+        //- Name of the predicted flux transporting the field
+        word phi0Name_;
+
+        //- Name of the flux transporting the field
+        word phiName_;
+
+        //- Name of the density field used to normalise the mass flux
+        //  if neccessary
+        word rhoName_;
+
+
+public:
+
+    //- Runtime type information
+    TypeName("multiphaseFixedFluxPressure");
+
+
+    // Constructors
+
+        //- Construct from patch and internal field
+        multiphaseFixedFluxPressureFvPatchScalarField
+        (
+            const fvPatch&,
+            const DimensionedField<scalar, volMesh>&
+        );
+
+        //- Construct from patch, internal field and dictionary
+        multiphaseFixedFluxPressureFvPatchScalarField
+        (
+            const fvPatch&,
+            const DimensionedField<scalar, volMesh>&,
+            const dictionary&
+        );
+
+        //- Construct by mapping given
+        //  multiphaseFixedFluxPressureFvPatchScalarField onto a new patch
+        multiphaseFixedFluxPressureFvPatchScalarField
+        (
+            const multiphaseFixedFluxPressureFvPatchScalarField&,
+            const fvPatch&,
+            const DimensionedField<scalar, volMesh>&,
+            const fvPatchFieldMapper&
+        );
+
+        //- Construct as copy
+        multiphaseFixedFluxPressureFvPatchScalarField
+        (
+            const multiphaseFixedFluxPressureFvPatchScalarField&
+        );
+
+        //- Construct and return a clone
+        virtual tmp<fvPatchScalarField> clone() const
+        {
+            return tmp<fvPatchScalarField>
+            (
+                new multiphaseFixedFluxPressureFvPatchScalarField(*this)
+            );
+        }
+
+        //- Construct as copy setting internal field reference
+        multiphaseFixedFluxPressureFvPatchScalarField
+        (
+            const multiphaseFixedFluxPressureFvPatchScalarField&,
+            const DimensionedField<scalar, volMesh>&
+        );
+
+        //- Construct and return a clone setting internal field reference
+        virtual tmp<fvPatchScalarField> clone
+        (
+            const DimensionedField<scalar, volMesh>& iF
+        ) const
+        {
+            return tmp<fvPatchScalarField>
+            (
+                new multiphaseFixedFluxPressureFvPatchScalarField(*this, iF)
+            );
+        }
+
+
+    // Member functions
+
+        //- Update the coefficients associated with the patch field
+        virtual void updateCoeffs();
+
+        //- Write
+        virtual void write(Ostream&) const;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

applications/solvers/multiphase/multiphaseEulerFoam/multiphaseSystem/multiphaseSystem.C

@@ -441,24 +441,15 @@ Foam::tmp<Foam::volScalarField> Foam::multiphaseSystem::K
 
 Foam::multiphaseSystem::multiphaseSystem
 (
-    const fvMesh& mesh,
+    const volVectorField& U,
     const surfaceScalarField& phi
 )
 :
-    IOdictionary
+    transportModel(U, phi),
-    (
-        IOobject
-        (
-            "transportProperties",
-            mesh.time().constant(),
-            mesh,
-            IOobject::MUST_READ_IF_MODIFIED,
-            IOobject::NO_WRITE
-        )
-    ),
-    phases_(lookup("phases"), phaseModel::iNew(mesh)),
 
-    mesh_(mesh),
+    phases_(lookup("phases"), phaseModel::iNew(U.mesh())),
+
+    mesh_(U.mesh()),
     phi_(phi),
 
     alphas_
@@ -514,7 +505,6 @@ Foam::multiphaseSystem::multiphaseSystem
 
 // * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
 
-
 Foam::tmp<Foam::volScalarField> Foam::multiphaseSystem::rho() const
 {
     PtrDictionary<phaseModel>::const_iterator iter = phases_.begin();
@@ -530,6 +520,21 @@ Foam::tmp<Foam::volScalarField> Foam::multiphaseSystem::rho() const
 }
 
 
+Foam::tmp<Foam::volScalarField> Foam::multiphaseSystem::nu() const
+{
+    PtrDictionary<phaseModel>::const_iterator iter = phases_.begin();
+
+    tmp<volScalarField> tnu = iter()*iter().nu();
+
+    for (++iter; iter != phases_.end(); ++iter)
+    {
+        tnu() += iter()*iter().nu();
+    }
+
+    return tnu;
+}
+
+
 Foam::tmp<Foam::volScalarField> Foam::multiphaseSystem::Cvm
 (
     const phaseModel& phase
@@ -770,19 +775,80 @@ void Foam::multiphaseSystem::solve()
 
     label nAlphaSubCycles(readLabel(pimpleDict.lookup("nAlphaSubCycles")));
 
-    volScalarField& alpha = phases_.first();
-
     if (nAlphaSubCycles > 1)
     {
         dimensionedScalar totalDeltaT = runTime.deltaT();
 
+        PtrList<volScalarField> alpha0s(phases_.size());
+        PtrList<surfaceScalarField> phiSums(phases_.size());
+
+        int phasei = 0;
+        forAllIter(PtrDictionary<phaseModel>, phases_, iter)
+        {
+            phaseModel& phase = iter();
+            volScalarField& alpha = phase;
+
+            alpha0s.set
+            (
+                phasei,
+                new volScalarField(alpha.oldTime())
+            );
+
+            phiSums.set
+            (
+                phasei,
+                new surfaceScalarField
+                (
+                    IOobject
+                    (
+                        "phiSum" + alpha.name(),
+                        runTime.timeName(),
+                        mesh_
+                    ),
+                    mesh_,
+                    dimensionedScalar("0", dimensionSet(0, 3, -1, 0, 0), 0)
+                )
+            );
+
+            phasei++;
+        }
+
         for
         (
-            subCycle<volScalarField> alphaSubCycle(alpha, nAlphaSubCycles);
+            subCycleTime alphaSubCycle
+            (
+                const_cast<Time&>(runTime),
+                nAlphaSubCycles
+            );
             !(++alphaSubCycle).end();
         )
         {
             solveAlphas();
+
+            int phasei = 0;
+            forAllIter(PtrDictionary<phaseModel>, phases_, iter)
+            {
+                phiSums[phasei] += (runTime.deltaT()/totalDeltaT)*iter().phi();
+                phasei++;
+            }
+        }
+
+        phasei = 0;
+        forAllIter(PtrDictionary<phaseModel>, phases_, iter)
+        {
+            phaseModel& phase = iter();
+            volScalarField& alpha = phase;
+
+            phase.phi() = phiSums[phasei];
+
+            // Correct the time index of the field to correspond to the global time
+            alpha.timeIndex() = runTime.timeIndex();
+
+            // Reset the old-time field value
+            alpha.oldTime() = alpha0s[phasei];
+            alpha.oldTime().timeIndex() = runTime.timeIndex();
+
+            phasei++;
         }
     }
     else

applications/solvers/multiphase/multiphaseEulerFoam/multiphaseSystem/multiphaseSystem.H

@@ -61,7 +61,7 @@ namespace Foam
 
 class multiphaseSystem
 :
-    public IOdictionary
+    public transportModel
 {
 
 public:
@@ -247,7 +247,7 @@ public:
         //- Construct from components
         multiphaseSystem
         (
-            const fvMesh& mesh,
+            const volVectorField& U,
             const surfaceScalarField& phi
         );
 
@@ -274,6 +274,9 @@ public:
         //- Return the mixture density
         tmp<volScalarField> rho() const;
 
+        //- Return the mixture laminar viscosity
+        tmp<volScalarField> nu() const;
+
         //- Return the virtual-mass coefficient for the given phase
         tmp<volScalarField> Cvm(const phaseModel& phase) const;
 
@@ -305,6 +308,10 @@ public:
         //- Solve for the mixture phase-fractions
         void solve();
 
+        //- Dummy correct
+        void correct()
+        {}
+
         //- Read base transportProperties dictionary
         bool read();
 };

applications/solvers/multiphase/multiphaseEulerFoam/pEqn.H

@@ -41,6 +41,20 @@
         phasei++;
     }
 
+    surfaceScalarField phi0
+    (
+        IOobject
+        (
+            "phi0",
+            runTime.timeName(),
+            mesh,
+            IOobject::NO_READ,
+            IOobject::AUTO_WRITE
+        ),
+        mesh,
+        dimensionedScalar("phi0", dimArea*dimVelocity, 0)
+    );
+
     phasei = 0;
     forAllIter(PtrDictionary<phaseModel>, fluid.phases(), iter)
     {
@@ -60,7 +74,8 @@
           + fvc::ddtPhiCorr(rAUs[phasei], alpha, phase.U(), phase.phi())
         );
         mrfZones.relativeFlux(phase.phi());
-        phase.phi() += rAlphaAUfs[phasei]*(g & mesh.Sf());
+        surfaceScalarField pphi0("pphi0", phase.phi());
+        pphi0 += rAlphaAUfs[phasei]*(g & mesh.Sf());
 
         multiphaseSystem::dragModelTable::const_iterator dmIter =
             fluid.dragModels().begin();
@@ -99,13 +114,16 @@
                     phasej++;
                 }
 
-                phase.phi() +=
+                pphi0 +=
                     fvc::interpolate
                     ((1.0/phase.rho())*rAUs[phasei]*(*dcIter()))
                    *phi0s[phasej];
             }
         }
 
+        phi0 += alphafs[phasei]*pphi0;
+        phase.phi() = pphi0;
+
         phasei++;
     }
 
@@ -132,13 +150,13 @@
         phasei++;
     }
     Dp = mag(Dp);
-    adjustPhi(phi, U, p);
+    adjustPhi(phi0, U, p);
 
     for(int nonOrth=0; nonOrth<=pimple.nNonOrthCorr(); nonOrth++)
     {
         fvScalarMatrix pEqnIncomp
         (
-            fvc::div(phi)
+            fvc::div(phi0)
           - fvm::laplacian(Dp, p)
         );
 
@@ -169,6 +187,7 @@
             }
 
             // dgdt =
+
             // (
             //     pos(alpha2)*(pEqnComp2 & p)/rho2
             //   - pos(alpha1)*(pEqnComp1 & p)/rho1

applications/solvers/multiphase/multiphaseEulerFoam/zonePhaseVolumes.H

@@ -0,0 +1,26 @@
+{
+    const scalarField& V = mesh.V();
+
+    forAll(mesh.cellZones(), czi)
+    {
+        const labelList& cellLabels = mesh.cellZones()[czi];
+
+        forAllConstIter(PtrDictionary<phaseModel>, fluid.phases(), iter)
+        {
+            const volScalarField& alpha = iter();
+            scalar phaseVolume = 0;
+
+            forAll(cellLabels, cli)
+            {
+                label celli = cellLabels[cli];
+                phaseVolume += alpha[celli]*V[celli];
+            }
+
+            reduce(phaseVolume, sumOp<scalar>());
+
+            Info<< alpha.name()
+                << " phase volume in zone " << mesh.cellZones()[czi].name()
+                << " = " << phaseVolume*1e6 << " ml " << endl;
+        }
+    }
+}