Merge pull request #125 from ParticulateFlow/release

Release 21.11

applications/solvers/cfdemSolverMultiphase/multiphaseMixture/multiphaseMixture.C

@@ -55,14 +55,21 @@ Foam::multiphaseMixture::calcNu() const
 {
     PtrDictionary<phase>::const_iterator iter = phases_.begin();
 
+    // 1/nu
+    tmp<volScalarField> tnuInv = iter()/iter().nu();
+    volScalarField& nuInv = tnuInv.ref();
+
+    // nu
     tmp<volScalarField> tnu = iter()*iter().nu();
     volScalarField& nu = tnu.ref();
 
     for (++iter; iter != phases_.end(); ++iter)
     {
-        nu += iter()*iter().nu();
+        nuInv += iter()/iter().nu();
     }
 
+    nu = 1/nuInv;
+    
     return tnu;
 }
 

applications/solvers/cfdemSolverMultiphaseScalar/Allwclean

@@ -0,0 +1,8 @@
+#!/bin/sh
+cd ${0%/*} || exit 1    # Run from this directory
+set -x
+
+wclean libso multiphaseMixture
+wclean
+
+#------------------------------------------------------------------------------

applications/solvers/cfdemSolverMultiphaseScalar/Allwmake

@@ -0,0 +1,12 @@
+#!/bin/sh
+cd ${0%/*} || exit 1    # Run from this directory
+
+# Parse arguments for library compilation
+targetType=libso
+. $WM_PROJECT_DIR/wmake/scripts/AllwmakeParseArguments
+set -x
+
+wmake $targetType multiphaseMixture
+wmake
+
+#------------------------------------------------------------------------------

applications/solvers/cfdemSolverMultiphaseScalar/CEqn.H

@@ -0,0 +1,22 @@
+        // get mixture properties
+        Cs = mixture.Cs();
+        diffusionCorrection = mixture.diffusionCorrection();
+        Deff = particleCloud.diffCoeffM().diffCoeff();
+
+        // get scalar source from DEM
+        particleCloud.massContributions(Sm);
+        particleCloud.massCoefficients(Smi);
+
+        fvScalarMatrix CEqn
+        (
+            fvm::ddt(voidfraction,C)
+          + fvm::div(phi,C)
+          - fvm::laplacian(Deff*voidfraction,C)
+          + fvm::div(fvc::interpolate(Deff*voidfraction)*diffusionCorrection*mesh.magSf(), C)
+          ==
+            Sm + fvm::Sp(Smi,C)
+        );
+
+        CEqn.relax();
+        fvOptions.constrain(CEqn);
+        CEqn.solve();

applications/solvers/cfdemSolverMultiphaseScalar/EEqn.H

@@ -0,0 +1,22 @@
+        // get mixture properties
+        Cp = mixture.Cp();
+        kf = mixture.kf();
+
+        // get scalar source from DEM
+        particleCloud.energyContributions(Qsource);
+        particleCloud.energyCoefficients(QCoeff);
+
+        fvScalarMatrix EEqn
+        (
+            rho*Cp*(fvm::ddt(voidfraction,T)
+          + fvm::div(phi,T))
+          - fvm::laplacian(thCond*voidfraction,T)
+          ==
+            Qsource + fvm::Sp(QCoeff,T)
+        );
+
+
+        EEqn.relax();
+        fvOptions.constrain(EEqn);
+        EEqn.solve();
+

applications/solvers/cfdemSolverMultiphaseScalar/Make/files

@@ -0,0 +1,3 @@
+cfdemSolverMultiphaseScalar.C
+
+EXE = $(CFDEM_APP_DIR)/cfdemSolverMultiphaseScalar

applications/solvers/cfdemSolverMultiphaseScalar/Make/options

@@ -0,0 +1,35 @@
+FOAM_VERSION_MAJOR := $(word 1,$(subst ., ,$(WM_PROJECT_VERSION)))
+PFLAGS+= -DOPENFOAM_VERSION_MAJOR=$(FOAM_VERSION_MAJOR)
+
+include $(CFDEM_ADD_LIBS_DIR)/additionalLibs
+
+EXE_INC = \
+     $(PFLAGS) \
+    -I$(CFDEM_OFVERSION_DIR) \
+    -ImultiphaseMixture/lnInclude \
+    -I$(LIB_SRC)/transportModels \
+    -I$(LIB_SRC)/transportModels/incompressible/lnInclude \
+    -I$(LIB_SRC)/transportModels/interfaceProperties/lnInclude \
+    -I$(LIB_SRC)/TurbulenceModels/turbulenceModels/lnInclude \
+    -I$(LIB_SRC)/TurbulenceModels/incompressible/lnInclude \
+    -I$(LIB_SRC)/finiteVolume/lnInclude \
+    -I$(LIB_SRC)/meshTools/lnInclude \
+    -I$(LIB_SRC)/sampling/lnInclude \
+    -I$(CFDEM_SRC_DIR)/lagrangian/cfdemParticle/lnInclude \
+    -I$(CFDEM_SRC_DIR)/lagrangian/cfdemParticle/cfdTools \
+    -Wno-deprecated-copy
+
+EXE_LIBS = \
+    -L$(CFDEM_LIB_DIR)\
+    -lcfdemMultiphaseInterFoamScalar \
+    -linterfaceProperties \
+    -lincompressibleTransportModels \
+    -lturbulenceModels \
+    -lincompressibleTurbulenceModels \
+    -lfiniteVolume \
+    -lfvOptions \
+    -lmeshTools \
+    -lsampling \
+    -l$(CFDEM_LIB_NAME) \
+     $(CFDEM_ADD_LIB_PATHS) \
+     $(CFDEM_ADD_LIBS)

applications/solvers/cfdemSolverMultiphaseScalar/UEqn.H

@@ -0,0 +1,61 @@
+const surfaceScalarField& rhoPhi(mixture.rhoPhi());
+
+volScalarField muEff = rho*(turbulence->nu() + turbulence->nut());
+
+if (modelType == "A")
+    muEff *= voidfraction;
+
+fvVectorMatrix UEqn
+(
+    fvm::ddt(rhoEps, U) - fvm::Sp(fvc::ddt(rhoEps),U)
+  + fvm::div(rhoPhi, U) - fvm::Sp(fvc::div(rhoPhi),U)
+  //+ particleCloud.divVoidfractionTau(U, voidfraction)
+  - fvm::laplacian(muEff, U) - fvc::div(muEff*dev2(fvc::grad(U)().T()))
+  ==
+    fvOptions(rho, U)
+  - fvm::Sp(Ksl,U)
+);
+
+UEqn.relax();
+
+fvOptions.constrain(UEqn);
+
+if (pimple.momentumPredictor() && (modelType=="B" || modelType=="Bfull"))
+{
+    solve
+    (
+        UEqn
+     ==
+        fvc::reconstruct
+        (
+            (- ghf*fvc::snGrad(rho) - fvc::snGrad(p_rgh)) * mesh.magSf()
+        )
+     +
+        fvc::reconstruct
+        (
+            mixture.surfaceTensionForce() * mesh.magSf()
+        ) * voidfraction
+     + Ksl*Us
+    );
+
+    fvOptions.correct(U);
+}
+else if (pimple.momentumPredictor())
+{
+    solve
+    (
+        UEqn
+     ==
+        fvc::reconstruct
+        (
+            (
+                mixture.surfaceTensionForce()
+              - ghf*fvc::snGrad(rho)
+              - fvc::snGrad(p_rgh)
+            ) * mesh.magSf()
+        ) * voidfraction
+     + Ksl*Us
+    );
+
+    fvOptions.correct(U);
+}

applications/solvers/cfdemSolverMultiphaseScalar/additionalChecks.H

@@ -0,0 +1,17 @@
+// Additional solver-specific checks
+
+// Useful if one wants to e.g. initialize floating particles using the Archimedes model
+if (particleCloud.couplingProperties().found("unrestrictedForceModelSelection"))
+{
+    Warning << "Using unrestrictedForceModelSelection, results may be incorrect!" << endl;
+} else
+{
+    #include "checkModelType.H"
+}
+
+word modelType = particleCloud.modelType();
+
+if(!particleCloud.couplingProperties().found("useDDTvoidfraction"))
+{
+    Warning << "Suppressing ddt(voidfraction) is not recommended with this solver as it may generate incorrect results!" << endl;
+}

applications/solvers/cfdemSolverMultiphaseScalar/alphaCourantNo.H

@@ -0,0 +1,21 @@
+scalar alphaCoNum = 0.0;
+scalar meanAlphaCoNum = 0.0;
+
+if (mesh.nInternalFaces())
+{
+    scalarField sumPhi
+    (
+        mixture.nearInterface()().primitiveField()
+       *fvc::surfaceSum(mag(phi))().primitiveField()
+    );
+
+    alphaCoNum = 0.5*gMax(sumPhi/mesh.V().field())*runTime.deltaTValue();
+
+    meanAlphaCoNum =
+        0.5*(gSum(sumPhi)/gSum(mesh.V().field()))*runTime.deltaTValue();
+}
+
+Info<< "Interface Courant Number mean: " << meanAlphaCoNum
+    << " max: " << alphaCoNum << endl;
+
+// ************************************************************************* //

applications/solvers/cfdemSolverMultiphaseScalar/cfdemSolverMultiphaseScalar.C

@@ -0,0 +1,164 @@
+/*---------------------------------------------------------------------------*\
+License
+
+    This 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.
+
+    This code 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 this code.  If not, see <http://www.gnu.org/licenses/>.
+
+    Copyright (C) 2018- Mathias Vångö, JKU Linz, Austria
+
+Application
+    cfdemSolverMultiphaseScalar
+
+Description
+    CFD-DEM solver for n incompressible fluids which captures the interfaces and
+    includes surface-tension and contact-angle effects for each phase. It is based
+    on the OpenFOAM(R)-4.x solver multiphaseInterFoam but extended to incorporate
+    DEM functionalities from the open-source DEM code LIGGGHTS.
+
+    Turbulence modelling is generic, i.e. laminar, RAS or LES may be selected.
+
+\*---------------------------------------------------------------------------*/
+
+#include "fvCFD.H"
+#include "multiphaseMixture.H"
+#include "turbulentTransportModel.H"
+#include "pimpleControl.H"
+#include "fvOptions.H"
+#include "CorrectPhi.H"
+
+#include "cfdemCloudEnergy.H"
+#include "implicitCouple.H"
+#include "clockModel.H"
+#include "smoothingModel.H"
+#include "forceModel.H"
+#include "thermCondModel.H"
+#include "diffCoeffModel.H"
+#include "energyModel.H"
+#include "massTransferModel.H"
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+int main(int argc, char *argv[])
+{
+#if OPENFOAM_VERSION_MAJOR >= 6
+    FatalError << "cfdemSolverMultiphase requires OpenFOAM 4.x or 5.x to work properly" << exit(FatalError);
+#endif
+
+    #include "postProcess.H"
+    #include "setRootCase.H"
+    #include "createTime.H"
+    #include "createMesh.H"
+    #include "createControl.H"
+    #include "initContinuityErrs.H"
+    #include "createFields.H"
+    #include "createFvOptions.H"
+    #include "correctPhi.H"
+    #include "CourantNo.H"
+
+    turbulence->validate();
+
+    // create cfdemCloud
+    cfdemCloudEnergy particleCloud(mesh);
+
+    #include "additionalChecks.H"
+
+    // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+    Info<< "\nStarting time loop\n" << endl;
+
+    while (runTime.loop())
+    {
+        #include "CourantNo.H"
+        #include "alphaCourantNo.H"
+
+        particleCloud.clockM().start(1,"Global");
+
+        Info<< "Time = " << runTime.timeName() << nl << endl;
+
+        particleCloud.clockM().start(2,"Coupling");
+        bool hasEvolved = particleCloud.evolve(voidfraction,Us,U);
+
+        if(hasEvolved)
+        {
+            particleCloud.smoothingM().smoothen(particleCloud.forceM(0).impParticleForces());
+        }
+
+        Info << "update Ksl.internalField()" << endl;
+        Ksl = particleCloud.momCoupleM(0).impMomSource();
+        Ksl.correctBoundaryConditions();
+
+       //Force Checks
+       vector fTotal(0,0,0);
+       vector fImpTotal = sum(mesh.V()*Ksl.internalField()*(Us.internalField()-U.internalField())).value();
+       reduce(fImpTotal, sumOp<vector>());
+       Info << "TotalForceExp: " << fTotal << endl;
+       Info << "TotalForceImp: " << fImpTotal << endl;
+
+        #include "solverDebugInfo.H"
+        particleCloud.clockM().stop("Coupling");
+
+        particleCloud.clockM().start(26,"Flow");
+
+        if(particleCloud.solveFlow())
+        {
+            mixture.solve();
+            rho = mixture.rho();
+            rhoEps = rho * voidfraction;
+
+            #include "EEqn.H"
+            #include "CEqn.H"
+
+             // --- Pressure-velocity PIMPLE corrector loop
+            while (pimple.loop())
+            {
+                #include "UEqn.H"
+
+                // --- Pressure corrector loop
+                while (pimple.correct())
+                {
+                    #include "pEqn.H"
+                }
+
+                if (pimple.turbCorr())
+                {
+                    turbulence->correct();
+                }
+            }
+        }
+        else
+        {
+            Info << "skipping flow solution." << endl;
+        }
+
+        particleCloud.clockM().start(31,"postFlow");
+        particleCloud.postFlow();
+        particleCloud.clockM().stop("postFlow");
+
+        runTime.write();
+
+        Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
+            << "  ClockTime = " << runTime.elapsedClockTime() << " s"
+            << nl << endl;
+
+        particleCloud.clockM().stop("Flow");
+        particleCloud.clockM().stop("Global");
+    }
+
+    Info<< "End\n" << endl;
+
+    return 0;
+}
+
+
+// ************************************************************************* //

applications/solvers/cfdemSolverMultiphaseScalar/correctPhi.H

@@ -0,0 +1,11 @@
+CorrectPhi
+(
+    U,
+    phi,
+    p_rgh,
+    dimensionedScalar("rAUf", dimTime/rho.dimensions(), 1),
+    geometricZeroField(),
+    pimple
+);
+
+#include "continuityErrs.H"

applications/solvers/cfdemSolverMultiphaseScalar/createFields.H

@@ -0,0 +1,342 @@
+//===============================
+// particle interaction modelling
+//===============================
+
+Info<< "\nReading momentum exchange field Ksl\n" << endl;
+volScalarField Ksl
+(
+    IOobject
+    (
+        "Ksl",
+        runTime.timeName(),
+        mesh,
+        IOobject::MUST_READ,
+        IOobject::AUTO_WRITE
+        ),
+    mesh
+    //dimensionedScalar("0", dimensionSet(1, -3, -1, 0, 0), 1.0)
+);
+
+Info<< "\nReading voidfraction field voidfraction = (Vgas/Vparticle)\n" << endl;
+volScalarField voidfraction
+(
+    IOobject
+    (
+        "voidfraction",
+        runTime.timeName(),
+        mesh,
+        IOobject::MUST_READ,
+        IOobject::AUTO_WRITE
+        ),
+    mesh
+);
+voidfraction.oldTime();
+
+Info<< "Reading particle velocity field Us\n" << endl;
+volVectorField Us
+(
+    IOobject
+    (
+        "Us",
+        runTime.timeName(),
+        mesh,
+        IOobject::MUST_READ,
+        IOobject::AUTO_WRITE
+        ),
+    mesh
+);
+
+Info<< "Reading field p_rgh\n" << endl;
+volScalarField p_rgh
+(
+    IOobject
+    (
+        "p_rgh",
+        runTime.timeName(),
+        mesh,
+        IOobject::MUST_READ,
+        IOobject::AUTO_WRITE
+    ),
+    mesh
+);
+
+Info<< "Reading field U\n" << endl;
+volVectorField U
+(
+    IOobject
+    (
+        "U",
+        runTime.timeName(),
+        mesh,
+        IOobject::MUST_READ,
+        IOobject::AUTO_WRITE
+    ),
+    mesh
+);
+
+Info<< "Reading/calculating face flux field phi\n" << endl;
+surfaceScalarField phi
+ (
+     IOobject
+     (
+        "phi",
+        runTime.timeName(),
+         mesh,
+        IOobject::READ_IF_PRESENT,
+        IOobject::AUTO_WRITE
+     ),
+     linearInterpolate(U*voidfraction) & mesh.Sf()
+ );
+
+multiphaseMixture mixture(U, phi, voidfraction);
+
+// Need to store rho for ddt(rho, U)
+volScalarField rho
+(
+    IOobject
+    (
+        "rho",
+        runTime.timeName(),
+        mesh,
+        IOobject::READ_IF_PRESENT,
+        IOobject::AUTO_WRITE
+    ),
+    mixture.rho()
+);
+rho.oldTime();
+
+//========================
+// scalar field modelling
+//========================
+Info<< "Reading/creating thermal fields\n" << endl;
+volScalarField T
+(
+    IOobject
+    (
+        "T",
+        runTime.timeName(),
+        mesh,
+        IOobject::MUST_READ,
+        IOobject::AUTO_WRITE
+    ),
+    mesh
+);
+
+volScalarField Qsource
+(
+    IOobject
+    (
+        "Qsource",
+        runTime.timeName(),
+        mesh,
+        IOobject::NO_READ,
+        IOobject::AUTO_WRITE
+    ),
+    mesh,
+    dimensionedScalar("zero", dimensionSet(1,-1,-3,0,0,0,0), 0.0)
+);
+
+volScalarField QCoeff
+(
+    IOobject
+    (
+        "Qsource",
+        runTime.timeName(),
+        mesh,
+        IOobject::NO_READ,
+        IOobject::AUTO_WRITE
+    ),
+    mesh,
+    dimensionedScalar("zero", dimensionSet(1,-1,-3,-1,0,0,0), 0.0)
+);
+
+volScalarField Cp
+(
+    IOobject
+    (
+        "Cp",
+        runTime.timeName(),
+        mesh,
+        IOobject::READ_IF_PRESENT,
+        IOobject::AUTO_WRITE
+    ),
+    mixture.Cp()
+);
+
+volScalarField kf
+(
+    IOobject
+    (
+        "kf",
+        runTime.timeName(),
+        mesh,
+        IOobject::READ_IF_PRESENT,
+        IOobject::AUTO_WRITE
+    ),
+    mixture.kf()
+);
+
+volScalarField thCond
+(
+    IOobject
+    (
+        "thCond",
+        runTime.timeName(),
+        mesh,
+        IOobject::READ_IF_PRESENT,
+        IOobject::AUTO_WRITE
+    ),
+    mesh,
+    dimensionedScalar("zero", dimensionSet(1,1,-3,-1,0,0,0), 0.0),
+    "zeroGradient"
+);
+
+Info<< "Reading/creating concentration fields\n" << endl;
+volScalarField C
+(
+    IOobject
+    (
+        "C",
+        runTime.timeName(),
+        mesh,
+        IOobject::MUST_READ,
+        IOobject::AUTO_WRITE
+    ),
+    mesh
+);
+
+volScalarField Sm
+(
+    IOobject
+    (
+        "Sm",
+        runTime.timeName(),
+        mesh,
+        IOobject::NO_READ,
+        IOobject::AUTO_WRITE
+    ),
+    mesh,
+    dimensionedScalar("zero", dimensionSet(1,-3,-1,0,0,0,0), 0.0)
+);
+
+volScalarField Smi
+(
+    IOobject
+    (
+        "Smi",
+        runTime.timeName(),
+        mesh,
+        IOobject::NO_READ,
+        IOobject::AUTO_WRITE
+    ),
+    mesh,
+    dimensionedScalar("zero", dimensionSet(0,0,-1,0,0,0,0), 0.0)
+);
+
+
+volScalarField D
+(
+    IOobject
+    (
+        "D",
+        runTime.timeName(),
+        mesh,
+        IOobject::READ_IF_PRESENT,
+        IOobject::AUTO_WRITE
+    ),
+    mixture.D()
+);
+
+volScalarField Deff
+(
+    IOobject
+    (
+        "Deff",
+        runTime.timeName(),
+        mesh,
+        IOobject::READ_IF_PRESENT,
+        IOobject::AUTO_WRITE
+    ),
+    mesh,
+    dimensionedScalar("zero", dimensionSet(0,2,-1,0,0,0,0), 0.0)
+);
+
+volScalarField Cs
+(
+    IOobject
+    (
+        "Cs",
+        runTime.timeName(),
+        mesh,
+        IOobject::READ_IF_PRESENT,
+        IOobject::AUTO_WRITE
+    ),
+    mixture.Cs()
+);
+
+surfaceScalarField diffusionCorrection
+(
+    IOobject
+    (
+        "diffusionCorrection",
+        runTime.timeName(),
+        mesh,
+        IOobject::NO_READ,
+        IOobject::AUTO_WRITE
+    ),
+    mixture.diffusionCorrection()
+);
+
+//========================
+
+volScalarField rhoEps ("rhoEps", rho * voidfraction);
+
+// Construct incompressible turbulence model
+autoPtr<incompressible::turbulenceModel> turbulence
+(
+    incompressible::turbulenceModel::New(U, phi, mixture)
+);
+
+
+#include "readGravitationalAcceleration.H"
+#include "readhRef.H"
+#include "gh.H"
+
+
+volScalarField p
+(
+    IOobject
+    (
+        "p",
+        runTime.timeName(),
+        mesh,
+        IOobject::NO_READ,
+        IOobject::AUTO_WRITE
+    ),
+    p_rgh + rho*gh
+);
+
+label pRefCell = 0;
+scalar pRefValue = 0.0;
+setRefCell
+(
+    p,
+    p_rgh,
+    pimple.dict(),
+    pRefCell,
+    pRefValue
+);
+
+if (p_rgh.needReference())
+{
+    p += dimensionedScalar
+    (
+        "p",
+        p.dimensions(),
+        pRefValue - getRefCellValue(p, pRefCell)
+    );
+}
+
+mesh.setFluxRequired(p_rgh.name());
+
+

applications/solvers/cfdemSolverMultiphaseScalar/multiphaseMixture/Make/files

@@ -0,0 +1,5 @@
+phase/phase.C
+alphaContactAngle/alphaContactAngleFvPatchScalarField.C
+multiphaseMixture.C
+
+LIB = $(CFDEM_LIB_DIR)/libcfdemMultiphaseInterFoamScalar

applications/solvers/cfdemSolverMultiphaseScalar/multiphaseMixture/Make/options

@@ -0,0 +1,18 @@
+FOAM_VERSION_MAJOR := $(word 1,$(subst ., ,$(WM_PROJECT_VERSION)))
+PFLAGS+= -DOPENFOAM_VERSION_MAJOR=$(FOAM_VERSION_MAJOR)
+
+EXE_INC = \
+     $(PFLAGS) \
+    -IalphaContactAngle \
+    -I$(LIB_SRC)/transportModels \
+    -I$(LIB_SRC)/transportModels/incompressible/lnInclude \
+    -I$(LIB_SRC)/transportModels/interfaceProperties/lnInclude \
+    -I$(LIB_SRC)/finiteVolume/lnInclude \
+    -I$(LIB_SRC)/meshTools/lnInclude \
+    -Wno-deprecated-copy
+
+LIB_LIBS = \
+    -linterfaceProperties \
+    -lincompressibleTransportModels \
+    -lfiniteVolume \
+    -lmeshTools

applications/solvers/cfdemSolverMultiphaseScalar/multiphaseMixture/alphaContactAngle/alphaContactAngleFvPatchScalarField.C

@@ -0,0 +1,146 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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 "alphaContactAngleFvPatchScalarField.H"
+#include "addToRunTimeSelectionTable.H"
+#include "fvPatchFieldMapper.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+alphaContactAngleFvPatchScalarField::interfaceThetaProps::interfaceThetaProps
+(
+    Istream& is
+)
+:
+    theta0_(readScalar(is)),
+    uTheta_(readScalar(is)),
+    thetaA_(readScalar(is)),
+    thetaR_(readScalar(is))
+{}
+
+
+Istream& operator>>
+(
+    Istream& is,
+    alphaContactAngleFvPatchScalarField::interfaceThetaProps& tp
+)
+{
+    is >> tp.theta0_ >> tp.uTheta_ >> tp.thetaA_ >> tp.thetaR_;
+    return is;
+}
+
+
+Ostream& operator<<
+(
+    Ostream& os,
+    const alphaContactAngleFvPatchScalarField::interfaceThetaProps& tp
+)
+{
+    os  << tp.theta0_ << token::SPACE
+        << tp.uTheta_ << token::SPACE
+        << tp.thetaA_ << token::SPACE
+        << tp.thetaR_;
+
+    return os;
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+alphaContactAngleFvPatchScalarField::alphaContactAngleFvPatchScalarField
+(
+    const fvPatch& p,
+    const DimensionedField<scalar, volMesh>& iF
+)
+:
+    zeroGradientFvPatchScalarField(p, iF)
+{}
+
+
+alphaContactAngleFvPatchScalarField::alphaContactAngleFvPatchScalarField
+(
+    const alphaContactAngleFvPatchScalarField& gcpsf,
+    const fvPatch& p,
+    const DimensionedField<scalar, volMesh>& iF,
+    const fvPatchFieldMapper& mapper
+)
+:
+    zeroGradientFvPatchScalarField(gcpsf, p, iF, mapper),
+    thetaProps_(gcpsf.thetaProps_)
+{}
+
+
+alphaContactAngleFvPatchScalarField::alphaContactAngleFvPatchScalarField
+(
+    const fvPatch& p,
+    const DimensionedField<scalar, volMesh>& iF,
+    const dictionary& dict
+)
+:
+    zeroGradientFvPatchScalarField(p, iF),
+    thetaProps_(dict.lookup("thetaProperties"))
+{
+    evaluate();
+}
+
+
+alphaContactAngleFvPatchScalarField::alphaContactAngleFvPatchScalarField
+(
+    const alphaContactAngleFvPatchScalarField& gcpsf,
+    const DimensionedField<scalar, volMesh>& iF
+)
+:
+    zeroGradientFvPatchScalarField(gcpsf, iF),
+    thetaProps_(gcpsf.thetaProps_)
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+void alphaContactAngleFvPatchScalarField::write(Ostream& os) const
+{
+    fvPatchScalarField::write(os);
+    os.writeKeyword("thetaProperties")
+        << thetaProps_ << token::END_STATEMENT << nl;
+    writeEntry("value", os);
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+makePatchTypeField
+(
+    fvPatchScalarField,
+    alphaContactAngleFvPatchScalarField
+);
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// ************************************************************************* //

applications/solvers/cfdemSolverMultiphaseScalar/multiphaseMixture/alphaContactAngle/alphaContactAngleFvPatchScalarField.H

@@ -0,0 +1,215 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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::alphaContactAngleFvPatchScalarField
+
+Description
+    Contact-angle boundary condition for multi-phase interface-capturing
+    simulations.  Used in conjuction with multiphaseMixture.
+
+SourceFiles
+    alphaContactAngleFvPatchScalarField.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef alphaContactAngleFvPatchScalarField_H
+#define alphaContactAngleFvPatchScalarField_H
+
+#include "zeroGradientFvPatchFields.H"
+#include "multiphaseMixture.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                           Class alphaContactAngleFvPatch Declaration
+\*---------------------------------------------------------------------------*/
+
+class alphaContactAngleFvPatchScalarField
+:
+    public zeroGradientFvPatchScalarField
+{
+public:
+
+    class interfaceThetaProps
+    {
+        //- Equilibrium contact angle
+        scalar theta0_;
+
+        //- Dynamic contact angle velocity scale
+        scalar uTheta_;
+
+        //- Limiting advancing contact angle
+        scalar thetaA_;
+
+        //- Limiting receeding contact angle
+        scalar thetaR_;
+
+
+    public:
+
+        // Constructors
+            interfaceThetaProps()
+            {}
+
+            interfaceThetaProps(Istream&);
+
+
+        // Member functions
+
+            //- Return the equilibrium contact angle theta0
+            scalar theta0(bool matched=true) const
+            {
+                if (matched) return theta0_;
+                else return 180.0 - theta0_;
+            }
+
+            //- Return the dynamic contact angle velocity scale
+            scalar uTheta() const
+            {
+                return uTheta_;
+            }
+
+            //- Return the limiting advancing contact angle
+            scalar thetaA(bool matched=true) const
+            {
+                if (matched) return thetaA_;
+                else return 180.0 - thetaA_;
+            }
+
+            //- Return the limiting receeding contact angle
+            scalar thetaR(bool matched=true) const
+            {
+                if (matched) return thetaR_;
+                else return 180.0 - thetaR_;
+            }
+
+
+        // IO functions
+
+            friend Istream& operator>>(Istream&, interfaceThetaProps&);
+            friend Ostream& operator<<(Ostream&, const interfaceThetaProps&);
+    };
+
+    typedef HashTable
+    <
+        interfaceThetaProps,
+        multiphaseMixture::interfacePair,
+        multiphaseMixture::interfacePair::hash
+    > thetaPropsTable;
+
+
+private:
+
+    // Private data
+
+        thetaPropsTable thetaProps_;
+
+
+public:
+
+    //- Runtime type information
+    TypeName("alphaContactAngle");
+
+
+    // Constructors
+
+        //- Construct from patch and internal field
+        alphaContactAngleFvPatchScalarField
+        (
+            const fvPatch&,
+            const DimensionedField<scalar, volMesh>&
+        );
+
+        //- Construct from patch, internal field and dictionary
+        alphaContactAngleFvPatchScalarField
+        (
+            const fvPatch&,
+            const DimensionedField<scalar, volMesh>&,
+            const dictionary&
+        );
+
+        //- Construct by mapping given alphaContactAngleFvPatchScalarField
+        //  onto a new patch
+        alphaContactAngleFvPatchScalarField
+        (
+            const alphaContactAngleFvPatchScalarField&,
+            const fvPatch&,
+            const DimensionedField<scalar, volMesh>&,
+            const fvPatchFieldMapper&
+        );
+
+        //- Construct and return a clone
+        virtual tmp<fvPatchScalarField> clone() const
+        {
+            return tmp<fvPatchScalarField>
+            (
+                new alphaContactAngleFvPatchScalarField(*this)
+            );
+        }
+
+        //- Construct as copy setting internal field reference
+        alphaContactAngleFvPatchScalarField
+        (
+            const alphaContactAngleFvPatchScalarField&,
+            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 alphaContactAngleFvPatchScalarField(*this, iF)
+            );
+        }
+
+
+    // Member functions
+
+        //- Return the contact angle properties
+        const thetaPropsTable& thetaProps() const
+        {
+            return thetaProps_;
+        }
+
+        //- Write
+        virtual void write(Ostream&) const;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

applications/solvers/cfdemSolverMultiphaseScalar/multiphaseMixture/multiphaseMixture.C

@@ -0,0 +1,929 @@
+/*---------------------------------------------------------------------------*\
+License
+
+    This 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.
+
+    This code 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 this code.  If not, see <http://www.gnu.org/licenses/>.
+
+    Copyright (C) 2018- Mathias Vångö, JKU Linz, Austria
+
+\*---------------------------------------------------------------------------*/
+
+#include "multiphaseMixture.H"
+#include "alphaContactAngleFvPatchScalarField.H"
+#include "Time.H"
+#include "subCycle.H"
+#include "MULES.H"
+#include "surfaceInterpolate.H"
+#include "fvcGrad.H"
+#include "fvcSnGrad.H"
+#include "fvcDiv.H"
+#include "fvcFlux.H"
+
+// * * * * * * * * * * * * * * * Static Member Data  * * * * * * * * * * * * //
+
+const Foam::scalar Foam::multiphaseMixture::convertToRad =
+    Foam::constant::mathematical::pi/180.0;
+
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+void Foam::multiphaseMixture::calcAlphas()
+{
+    scalar level = 0.0;
+    alphas_ == 0.0;
+
+    forAllIter(PtrDictionary<phase>, phases_, iter)
+    {
+        alphas_ += level*iter();
+        level += 1.0;
+    }
+}
+
+
+Foam::tmp<Foam::volScalarField>
+Foam::multiphaseMixture::calcNu() const
+{
+    PtrDictionary<phase>::const_iterator iter = phases_.begin();
+
+    // 1/nu
+    tmp<volScalarField> tnuInv = iter()/iter().nu();
+    volScalarField& nuInv = tnuInv.ref();
+
+    // nu
+    tmp<volScalarField> tnu = iter()*iter().nu();
+    volScalarField& nu = tnu.ref();
+
+    for (++iter; iter != phases_.end(); ++iter)
+    {
+        nuInv += iter()/iter().nu();
+    }
+
+    nu = 1/nuInv;
+
+    return tnu;
+}
+
+Foam::tmp<Foam::surfaceScalarField>
+Foam::multiphaseMixture::calcStf() const
+{
+    tmp<surfaceScalarField> tstf
+    (
+        new surfaceScalarField
+        (
+            IOobject
+            (
+                "stf",
+                mesh_.time().timeName(),
+                mesh_
+            ),
+            mesh_,
+            dimensionedScalar
+            (
+                "stf",
+                dimensionSet(1, -2, -2, 0, 0),
+                0.0
+            )
+        )
+    );
+
+    surfaceScalarField& stf = tstf.ref();
+
+    forAllConstIter(PtrDictionary<phase>, phases_, iter1)
+    {
+        const phase& alpha1 = iter1();
+
+        PtrDictionary<phase>::const_iterator iter2 = iter1;
+        ++iter2;
+
+        for (; iter2 != phases_.end(); ++iter2)
+        {
+            const phase& alpha2 = iter2();
+
+            sigmaTable::const_iterator sigma =
+                sigmas_.find(interfacePair(alpha1, alpha2));
+
+            if (sigma == sigmas_.end())
+            {
+                FatalErrorInFunction
+                    << "Cannot find interface " << interfacePair(alpha1, alpha2)
+                    << " in list of sigma values"
+                    << exit(FatalError);
+            }
+
+            stf += dimensionedScalar("sigma", dimSigma_, sigma())
+               *fvc::interpolate(K(alpha1, alpha2))*
+                (
+                    fvc::interpolate(alpha2)*fvc::snGrad(alpha1)
+                  - fvc::interpolate(alpha1)*fvc::snGrad(alpha2)
+                );
+        }
+    }
+
+    return tstf;
+}
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::multiphaseMixture::multiphaseMixture
+(
+    const volVectorField& U,
+    const surfaceScalarField& phi,
+    const volScalarField& voidfraction
+)
+:
+    IOdictionary
+    (
+        IOobject
+        (
+            "transportProperties",
+            U.time().constant(),
+            U.db(),
+            IOobject::MUST_READ_IF_MODIFIED,
+            IOobject::NO_WRITE
+        )
+    ),
+
+    phases_(lookup("phases"), phase::iNew(U, phi)),
+
+    mesh_(U.mesh()),
+    U_(U),
+    phi_(phi),
+    voidfraction_(voidfraction),
+    rhoPhi_
+    (
+        IOobject
+        (
+            "rhoPhi",
+            mesh_.time().timeName(),
+            mesh_,
+            IOobject::NO_READ,
+            IOobject::NO_WRITE
+        ),
+        mesh_,
+        dimensionedScalar("rhoPhi", dimMass/dimTime, 0.0)
+    ),
+    surfaceTensionForce_
+    (
+        IOobject
+        (
+            "surfaceTensionForce",
+            mesh_.time().timeName(),
+            mesh_,
+            IOobject::NO_READ,
+            IOobject::AUTO_WRITE
+        ),
+        mesh_,
+        dimensionedScalar("surfaceTensionForce", dimensionSet(1, -2, -2, 0, 0), 0.0)
+    ),
+    alphas_
+    (
+        IOobject
+        (
+            "alphas",
+            mesh_.time().timeName(),
+            mesh_,
+            IOobject::NO_READ,
+            IOobject::AUTO_WRITE
+        ),
+        mesh_,
+        dimensionedScalar("alphas", dimless, 0.0)
+    ),
+
+    nu_
+    (
+        IOobject
+        (
+            "nu",
+            mesh_.time().timeName(),
+            mesh_,
+            IOobject::NO_READ,
+            IOobject::AUTO_WRITE
+        ),
+        calcNu()
+    ),
+
+    sigmas_(lookup("sigmas")),
+    dimSigma_(1, 0, -2, 0, 0),
+    deltaN_
+    (
+        "deltaN",
+        1e-8/pow(average(mesh_.V()), 1.0/3.0)
+    )
+{
+    calcAlphas();
+    alphas_.write();
+    surfaceTensionForce_ = calcStf();
+
+}
+
+
+// * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
+
+Foam::tmp<Foam::volScalarField>
+Foam::multiphaseMixture::rho() const
+{
+    PtrDictionary<phase>::const_iterator iter = phases_.begin();
+
+    tmp<volScalarField> trho = iter()*iter().rho();
+    volScalarField& rho = trho.ref();
+
+    for (++iter; iter != phases_.end(); ++iter)
+    {
+        rho += iter()*iter().rho();
+    }
+
+    return trho;
+}
+
+
+Foam::tmp<Foam::scalarField>
+Foam::multiphaseMixture::rho(const label patchi) const
+{
+    PtrDictionary<phase>::const_iterator iter = phases_.begin();
+
+    tmp<scalarField> trho = iter().boundaryField()[patchi]*iter().rho().value();
+    scalarField& rho = trho.ref();
+
+    for (++iter; iter != phases_.end(); ++iter)
+    {
+        rho += iter().boundaryField()[patchi]*iter().rho().value();
+    }
+
+    return trho;
+}
+
+
+Foam::tmp<Foam::volScalarField>
+Foam::multiphaseMixture::mu() const
+{
+    Info << "In multiphasemixture mu()" << endl;
+    return rho()*nu();
+//    PtrDictionary<phase>::const_iterator iter = phases_.begin();
+
+//    tmp<volScalarField> tmu = iter()*iter().rho()*iter().nu();
+//    volScalarField& mu = tmu.ref();
+
+//    for (++iter; iter != phases_.end(); ++iter)
+//    {
+//        mu += iter()*iter().rho()*iter().nu();
+//    }
+
+//    return tmu;
+}
+
+
+Foam::tmp<Foam::scalarField>
+Foam::multiphaseMixture::mu(const label patchi) const
+{
+    PtrDictionary<phase>::const_iterator iter = phases_.begin();
+
+    tmp<scalarField> tmu =
+        iter().boundaryField()[patchi]
+       *iter().rho().value()
+       *iter().nu(patchi);
+    scalarField& mu = tmu.ref();
+
+    for (++iter; iter != phases_.end(); ++iter)
+    {
+        mu +=
+            iter().boundaryField()[patchi]
+           *iter().rho().value()
+           *iter().nu(patchi);
+    }
+
+    return tmu;
+}
+
+
+Foam::tmp<Foam::surfaceScalarField>
+Foam::multiphaseMixture::muf() const
+{
+
+    return nuf()*fvc::interpolate(rho());
+//    PtrDictionary<phase>::const_iterator iter = phases_.begin();
+
+//    tmp<surfaceScalarField> tmuf =
+//            fvc::interpolate(iter())*iter().rho()*fvc::interpolate(iter().nu());
+//    surfaceScalarField& muf = tmuf.ref();
+
+//    for (++iter; iter != phases_.end(); ++iter)
+//    {
+//        muf +=
+//                fvc::interpolate(iter())*iter().rho()*fvc::interpolate(iter().nu());
+//    }
+
+//    return tmuf;
+}
+
+
+Foam::tmp<Foam::volScalarField>
+Foam::multiphaseMixture::nu() const
+{
+    return nu_;
+}
+
+Foam::tmp<Foam::scalarField>
+Foam::multiphaseMixture::nu(const label patchi) const
+{
+    //return nu_.boundaryField()[patchi];
+    PtrDictionary<phase>::const_iterator iter = phases_.begin();
+
+    tmp<scalarField> tnu =
+        iter().boundaryField()[patchi]
+       *iter().nu(patchi);
+    scalarField& nu = tnu.ref();
+
+    for (++iter; iter != phases_.end(); ++iter)
+    {
+        nu +=
+            iter().boundaryField()[patchi]
+           *iter().nu(patchi);
+    }
+
+    return tnu;
+}
+
+
+Foam::tmp<Foam::surfaceScalarField>
+Foam::multiphaseMixture::nuf() const
+{
+    //return muf()/fvc::interpolate(rho());
+    PtrDictionary<phase>::const_iterator iter = phases_.begin();
+
+    tmp<surfaceScalarField> tnuf =
+        fvc::interpolate(iter())*fvc::interpolate(iter().nu());
+    surfaceScalarField& nuf = tnuf.ref();
+
+    for (++iter; iter != phases_.end(); ++iter)
+    {
+        nuf +=
+            fvc::interpolate(iter())*fvc::interpolate(iter().nu());
+    }
+
+    return tnuf;
+}
+
+Foam::tmp<Foam::volScalarField>
+Foam::multiphaseMixture::Cp() const
+{
+    PtrDictionary<phase>::const_iterator iter = phases_.begin();
+
+    // rho*Cp
+    tmp<volScalarField> trhoCp = iter()*iter().Cp()*iter().rho();
+    volScalarField& rhoCp = trhoCp.ref();
+
+    // Cp
+    tmp<volScalarField> tCp = iter()*iter().Cp();
+    volScalarField& Cp = tCp.ref();
+
+    for (++iter; iter != phases_.end(); ++iter)
+    {
+        rhoCp += iter()*iter().Cp()*iter().rho();
+    }
+    Cp = rhoCp/rho();
+    return tCp;
+}
+
+Foam::tmp<Foam::volScalarField>
+Foam::multiphaseMixture::kf() const
+{
+    PtrDictionary<phase>::const_iterator iter = phases_.begin();
+
+    // rho*Cp/kf
+    tmp<volScalarField> trhoCpkf = iter()*iter().rho()*iter().Cp()/iter().kf();
+    volScalarField& rhoCpkf = trhoCpkf.ref();
+
+    // kf
+    tmp<volScalarField> tkf = iter()*iter().kf();
+    volScalarField& kf = tkf.ref();
+
+    for (++iter; iter != phases_.end(); ++iter)
+    {
+      rhoCpkf += iter()*iter().rho()*iter().Cp()/iter().kf();
+    }
+
+    kf = rho()*Cp()/rhoCpkf;
+    return tkf;
+}
+
+Foam::tmp<Foam::volScalarField>
+Foam::multiphaseMixture::D() const
+{
+    PtrDictionary<phase>::const_iterator iter = phases_.begin();
+
+    // 1/D
+    tmp<volScalarField> tDInv = iter()/iter().D();
+    volScalarField& DInv = tDInv.ref();
+
+    // D
+    tmp<volScalarField> tD = iter()*iter().D();
+    volScalarField& D = tD.ref();
+
+    for (++iter; iter != phases_.end(); ++iter)
+    {
+      DInv += iter()/iter().D();
+    }
+
+    D = 1/DInv;
+    return tD;
+}
+
+Foam::tmp<Foam::volScalarField>
+Foam::multiphaseMixture::Cs() const
+{
+    PtrDictionary<phase>::const_iterator iter = phases_.begin();
+
+    // Cs
+    tmp<volScalarField> tCs = iter()*iter().Cs();
+    volScalarField& Cs = tCs.ref();
+
+    for (++iter; iter != phases_.end(); ++iter)
+    {
+        Cs += iter()*iter().Cs();
+    }
+
+    return tCs;
+}
+
+Foam::tmp<Foam::surfaceScalarField>
+Foam::multiphaseMixture::diffusionCorrection() const
+{
+
+    surfaceScalarField numerator
+    (
+        IOobject
+        (
+            "numerator",
+            mesh_.time().timeName(),
+            mesh_,
+            IOobject::NO_READ,
+            IOobject::NO_WRITE
+        ),
+        mesh_,
+        dimensionedScalar("zero", dimless/dimLength, 0.0)
+    );
+
+    surfaceScalarField denominator
+    (
+        IOobject
+        (
+            "denominator",
+            mesh_.time().timeName(),
+            mesh_,
+            IOobject::NO_READ,
+            IOobject::NO_WRITE
+        ),
+        mesh_,
+        dimensionedScalar("zero", dimless, 0.0)
+    );
+
+    PtrDictionary<phase>::const_iterator iter = phases_.begin();
+    const phase& alpha1 = iter();
+
+    for (++iter; iter != phases_.end(); ++iter)
+    {
+        const phase& alpha2 = iter();
+        scalar He = alpha1.Cs().value() / (alpha2.Cs().value() + SMALL);
+
+        numerator   += (1/He - 1) * fvc::snGrad(alpha2);
+        denominator += fvc::interpolate(alpha2) * (1/He - 1);
+    }
+
+    tmp<surfaceScalarField> correction = numerator / (denominator + 1 + SMALL);
+
+    /*
+    PtrDictionary<phase>::const_iterator iter = phases_.begin();
+
+    const phase& alphaL = iter();
+    ++iter;
+    const phase& alphaG = iter();
+    scalar He = alphaG.Cs().value() / (alphaL.Cs().value() + SMALL);
+
+    surfaceScalarField gradAlphaL = fvc::snGrad(alphaL);
+    surfaceScalarField surfAlphaL = fvc::interpolate(alphaL);
+
+    tmp<surfaceScalarField> correction = (1-He)/(surfAlphaL + He*(1-surfAlphaL) + 10*SMALL) * gradAlphaL;
+    */
+    return correction;
+}
+
+void Foam::multiphaseMixture::solve()
+{
+    correct();
+
+    const Time& runTime = mesh_.time();
+
+    volScalarField& alpha = phases_.first();
+
+    const dictionary& alphaControls = mesh_.solverDict("alpha");
+    label nAlphaSubCycles(readLabel(alphaControls.lookup("nAlphaSubCycles")));
+    scalar cAlpha(readScalar(alphaControls.lookup("cAlpha")));
+
+    if (nAlphaSubCycles > 1)
+    {
+        surfaceScalarField rhoPhiSum
+        (
+            IOobject
+            (
+                "rhoPhiSum",
+                runTime.timeName(),
+                mesh_
+            ),
+            mesh_,
+            dimensionedScalar("0", rhoPhi_.dimensions(), 0)
+        );
+
+        dimensionedScalar totalDeltaT = runTime.deltaT();
+
+        for
+        (
+            subCycle<volScalarField> alphaSubCycle(alpha, nAlphaSubCycles);
+            !(++alphaSubCycle).end();
+        )
+        {
+            FatalError << "Sub-cycling of the alpha equation not yet implemented!!" << abort(FatalError);
+            solveAlphas(cAlpha);
+            rhoPhiSum += (runTime.deltaT()/totalDeltaT)*rhoPhi_;
+        }
+
+        rhoPhi_ = rhoPhiSum;
+    }
+    else
+    {
+        solveAlphas(cAlpha);
+    }
+
+    // Update the mixture kinematic viscosity
+    nu_ = calcNu();
+
+    surfaceTensionForce_ = calcStf();
+}
+
+
+void Foam::multiphaseMixture::correct()
+{
+    forAllIter(PtrDictionary<phase>, phases_, iter)
+    {
+        iter().correct();
+    }
+}
+
+
+Foam::tmp<Foam::surfaceVectorField> Foam::multiphaseMixture::nHatfv
+(
+    const volScalarField& alpha1,
+    const volScalarField& alpha2
+) const
+{
+    /*
+    // Cell gradient of alpha
+    volVectorField gradAlpha =
+        alpha2*fvc::grad(alpha1) - alpha1*fvc::grad(alpha2);
+
+    // Interpolated face-gradient of alpha
+    surfaceVectorField gradAlphaf = fvc::interpolate(gradAlpha);
+    */
+
+    surfaceVectorField gradAlphaf
+    (
+        fvc::interpolate(alpha2)*fvc::interpolate(fvc::grad(alpha1))
+      - fvc::interpolate(alpha1)*fvc::interpolate(fvc::grad(alpha2))
+    );
+
+    // Face unit interface normal
+    return gradAlphaf/(mag(gradAlphaf) + deltaN_);
+}
+
+
+Foam::tmp<Foam::surfaceScalarField> Foam::multiphaseMixture::nHatf
+(
+    const volScalarField& alpha1,
+    const volScalarField& alpha2
+) const
+{
+    // Face unit interface normal flux
+    return nHatfv(alpha1, alpha2) & mesh_.Sf();
+}
+
+
+// Correction for the boundary condition on the unit normal nHat on
+// walls to produce the correct contact angle.
+
+// The dynamic contact angle is calculated from the component of the
+// velocity on the direction of the interface, parallel to the wall.
+
+void Foam::multiphaseMixture::correctContactAngle
+(
+    const phase& alpha1,
+    const phase& alpha2,
+    surfaceVectorField::Boundary& nHatb
+) const
+{
+    const volScalarField::Boundary& gbf
+        = alpha1.boundaryField();
+
+    const fvBoundaryMesh& boundary = mesh_.boundary();
+
+    forAll(boundary, patchi)
+    {
+        if (isA<alphaContactAngleFvPatchScalarField>(gbf[patchi]))
+        {
+            const alphaContactAngleFvPatchScalarField& acap =
+                refCast<const alphaContactAngleFvPatchScalarField>(gbf[patchi]);
+
+            vectorField& nHatPatch = nHatb[patchi];
+
+            vectorField AfHatPatch
+            (
+                mesh_.Sf().boundaryField()[patchi]
+               /mesh_.magSf().boundaryField()[patchi]
+            );
+
+            alphaContactAngleFvPatchScalarField::thetaPropsTable::
+                const_iterator tp =
+                acap.thetaProps().find(interfacePair(alpha1, alpha2));
+
+            if (tp == acap.thetaProps().end())
+            {
+                FatalErrorInFunction
+                    << "Cannot find interface " << interfacePair(alpha1, alpha2)
+                    << "\n    in table of theta properties for patch "
+                    << acap.patch().name()
+                    << exit(FatalError);
+            }
+
+            bool matched = (tp.key().first() == alpha1.name());
+
+            scalar theta0 = convertToRad*tp().theta0(matched);
+            scalarField theta(boundary[patchi].size(), theta0);
+
+            scalar uTheta = tp().uTheta();
+
+            // Calculate the dynamic contact angle if required
+            if (uTheta > SMALL)
+            {
+                scalar thetaA = convertToRad*tp().thetaA(matched);
+                scalar thetaR = convertToRad*tp().thetaR(matched);
+
+                // Calculated the component of the velocity parallel to the wall
+                vectorField Uwall
+                (
+                    U_.boundaryField()[patchi].patchInternalField()
+                  - U_.boundaryField()[patchi]
+                );
+                Uwall -= (AfHatPatch & Uwall)*AfHatPatch;
+
+                // Find the direction of the interface parallel to the wall
+                vectorField nWall
+                (
+                    nHatPatch - (AfHatPatch & nHatPatch)*AfHatPatch
+                );
+
+                // Normalise nWall
+                nWall /= (mag(nWall) + SMALL);
+
+                // Calculate Uwall resolved normal to the interface parallel to
+                // the interface
+                scalarField uwall(nWall & Uwall);
+
+                theta += (thetaA - thetaR)*tanh(uwall/uTheta);
+            }
+
+
+            // Reset nHatPatch to correspond to the contact angle
+
+            scalarField a12(nHatPatch & AfHatPatch);
+
+            scalarField b1(cos(theta));
+
+            scalarField b2(nHatPatch.size());
+
+            forAll(b2, facei)
+            {
+                b2[facei] = cos(acos(a12[facei]) - theta[facei]);
+            }
+
+            scalarField det(1.0 - a12*a12);
+
+            scalarField a((b1 - a12*b2)/det);
+            scalarField b((b2 - a12*b1)/det);
+
+            nHatPatch = a*AfHatPatch + b*nHatPatch;
+
+            nHatPatch /= (mag(nHatPatch) + deltaN_.value());
+        }
+    }
+}
+
+
+Foam::tmp<Foam::volScalarField> Foam::multiphaseMixture::K
+(
+    const phase& alpha1,
+    const phase& alpha2
+) const
+{
+    tmp<surfaceVectorField> tnHatfv = nHatfv(alpha1, alpha2);
+
+    correctContactAngle(alpha1, alpha2, tnHatfv.ref().boundaryFieldRef());
+
+    // Simple expression for curvature
+    return -fvc::div(tnHatfv & mesh_.Sf());
+}
+
+
+Foam::tmp<Foam::volScalarField>
+Foam::multiphaseMixture::nearInterface() const
+{
+    tmp<volScalarField> tnearInt
+    (
+        new volScalarField
+        (
+            IOobject
+            (
+                "nearInterface",
+                mesh_.time().timeName(),
+                mesh_
+            ),
+            mesh_,
+            dimensionedScalar("nearInterface", dimless, 0.0)
+        )
+    );
+
+    forAllConstIter(PtrDictionary<phase>, phases_, iter)
+    {
+        tnearInt.ref() = max(tnearInt(), pos(iter() - 0.01)*pos(0.99 - iter()));
+    }
+
+    return tnearInt;
+}
+
+
+void Foam::multiphaseMixture::solveAlphas
+(
+    const scalar cAlpha
+)
+{
+    static label nSolves=-1;
+    nSolves++;
+
+    word alphaScheme("div(phi,alpha)");
+    word alpharScheme("div(phirb,alpha)");
+
+    surfaceScalarField phic(mag(phi_/mesh_.magSf()));
+    phic = min(cAlpha*phic, max(phic));
+
+    PtrList<surfaceScalarField> alphaPhiCorrs(phases_.size());
+    int phasei = 0;
+
+    forAllIter(PtrDictionary<phase>, phases_, iter)
+    {
+        phase& alpha = iter();
+
+        alphaPhiCorrs.set
+        (
+            phasei,
+            new surfaceScalarField
+            (
+                "phi" + alpha.name() + "Corr",
+                fvc::flux
+                (
+                    phi_,
+                    alpha,
+                    alphaScheme
+                )
+            )
+        );
+
+        surfaceScalarField& alphaPhiCorr = alphaPhiCorrs[phasei];
+
+        forAllIter(PtrDictionary<phase>, phases_, iter2)
+        {
+            phase& alpha2 = iter2();
+
+            if (&alpha2 == &alpha) continue;
+
+            surfaceScalarField phir(phic*nHatf(alpha, alpha2));
+
+            alphaPhiCorr += fvc::flux
+            (
+                -fvc::flux(-phir, alpha2, alpharScheme),
+                alpha,
+                alpharScheme
+            );
+        }
+
+        MULES::limit
+        (
+            1.0/mesh_.time().deltaT().value(),
+            voidfraction_,
+            alpha,
+            phi_,
+            alphaPhiCorr,
+            zeroField(),
+            zeroField(),
+#if OPENFOAM_VERSION_MAJOR < 6
+            1,
+            0,
+#else
+            oneField(),
+            zeroField(),
+#endif
+            true
+        );
+
+        phasei++;
+    }
+
+    MULES::limitSum(alphaPhiCorrs);
+
+    rhoPhi_ = dimensionedScalar("0", dimensionSet(1, 0, -1, 0, 0), 0);
+
+    volScalarField sumAlpha
+    (
+        IOobject
+        (
+            "sumAlpha",
+            mesh_.time().timeName(),
+            mesh_
+        ),
+        mesh_,
+        dimensionedScalar("sumAlpha", dimless, 0)
+    );
+
+    phasei = 0;
+
+    forAllIter(PtrDictionary<phase>, phases_, iter)
+    {
+        phase& alpha = iter();
+
+        surfaceScalarField& alphaPhi = alphaPhiCorrs[phasei];
+        alphaPhi += upwind<scalar>(mesh_, phi_).flux(alpha);
+
+        MULES::explicitSolve
+        (
+            voidfraction_,
+            alpha,
+            alphaPhi,
+            zeroField(),
+            zeroField()
+        );
+
+        rhoPhi_ += alphaPhi*alpha.rho();
+
+        Info<< alpha.name() << " volume fraction, min, max = "
+            << alpha.weightedAverage(mesh_.V()).value()
+            << ' ' << min(alpha).value()
+            << ' ' << max(alpha).value()
+            << endl;
+
+        sumAlpha += alpha;
+
+        phasei++;
+    }
+
+    Info<< "Phase-sum volume fraction, min, max = "
+        << sumAlpha.weightedAverage(mesh_.V()).value()
+        << ' ' << min(sumAlpha).value()
+        << ' ' << max(sumAlpha).value()
+        << endl;
+
+    calcAlphas();
+}
+
+
+bool Foam::multiphaseMixture::read()
+{
+    if (transportModel::read())
+    {
+        bool readOK = true;
+
+        PtrList<entry> phaseData(lookup("phases"));
+        label phasei = 0;
+
+        forAllIter(PtrDictionary<phase>, phases_, iter)
+        {
+            readOK &= iter().read(phaseData[phasei++].dict());
+        }
+
+        lookup("sigmas") >> sigmas_;
+
+        return readOK;
+    }
+    else
+    {
+        return false;
+    }
+}
+
+
+// ************************************************************************* //

applications/solvers/cfdemSolverMultiphaseScalar/multiphaseMixture/multiphaseMixture.H

@@ -0,0 +1,299 @@
+/*---------------------------------------------------------------------------*\
+License
+
+    This 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.
+
+    This code 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 this code.  If not, see <http://www.gnu.org/licenses/>.
+
+    Copyright (C) 2018- Mathias Vångö, JKU Linz, Austria
+
+Class
+    multiphaseMixture
+
+Description
+    This class is based on the OpenFOAM(R) Foam::multiphaseMixture class,
+    which is an incompressible multi-phase mixture with built in solution
+    for the phase fractions with interface compression for interface-capturing.
+    It has been extended to include the void fraction in the volume fraction
+    transport equations.
+
+    Derived from transportModel so that it can be unsed in conjunction with
+    the incompressible turbulence models.
+
+    Surface tension and contact-angle is handled for the interface
+    between each phase-pair.
+
+SourceFiles
+    multiphaseMixture.C
+\*---------------------------------------------------------------------------*/
+
+#ifndef multiphaseMixture_H
+#define multiphaseMixture_H
+
+#include "incompressible/transportModel/transportModel.H"
+#include "IOdictionary.H"
+#include "phase.H"
+#include "PtrDictionary.H"
+#include "volFields.H"
+#include "surfaceFields.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                      Class multiphaseMixture Declaration
+\*---------------------------------------------------------------------------*/
+
+class multiphaseMixture
+:
+    public IOdictionary,
+    public transportModel
+{
+public:
+
+    class interfacePair
+    :
+        public Pair<word>
+    {
+    public:
+
+        class hash
+        :
+            public Hash<interfacePair>
+        {
+        public:
+
+            hash()
+            {}
+
+            label operator()(const interfacePair& key) const
+            {
+                return word::hash()(key.first()) + word::hash()(key.second());
+            }
+        };
+
+
+        // Constructors
+
+            interfacePair()
+            {}
+
+            interfacePair(const word& alpha1Name, const word& alpha2Name)
+            :
+                Pair<word>(alpha1Name, alpha2Name)
+            {}
+
+            interfacePair(const phase& alpha1, const phase& alpha2)
+            :
+                Pair<word>(alpha1.name(), alpha2.name())
+            {}
+
+
+        // Friend Operators
+
+            friend bool operator==
+            (
+                const interfacePair& a,
+                const interfacePair& b
+            )
+            {
+                return
+                (
+                    ((a.first() == b.first()) && (a.second() == b.second()))
+                 || ((a.first() == b.second()) && (a.second() == b.first()))
+                );
+            }
+
+            friend bool operator!=
+            (
+                const interfacePair& a,
+                const interfacePair& b
+            )
+            {
+                return (!(a == b));
+            }
+    };
+
+
+private:
+
+    // Private data
+
+        //- Dictionary of phases
+        PtrDictionary<phase> phases_;
+
+        const fvMesh& mesh_;
+        const volVectorField& U_;
+        const surfaceScalarField& phi_;
+        const volScalarField& voidfraction_;
+        surfaceScalarField rhoPhi_;
+        surfaceScalarField surfaceTensionForce_;
+        volScalarField alphas_;
+
+        volScalarField nu_;
+
+        typedef HashTable<scalar, interfacePair, interfacePair::hash>
+            sigmaTable;
+
+        sigmaTable sigmas_;
+        dimensionSet dimSigma_;
+
+        //- Stabilisation for normalisation of the interface normal
+        const dimensionedScalar deltaN_;
+
+        //- Conversion factor for degrees into radians
+        static const scalar convertToRad;
+
+
+    // Private member functions
+
+        void calcAlphas();
+
+        tmp<volScalarField> calcNu() const;
+
+        void solveAlphas(const scalar cAlpha);
+
+        tmp<surfaceVectorField> nHatfv
+        (
+            const volScalarField& alpha1,
+            const volScalarField& alpha2
+        ) const;
+
+        tmp<surfaceScalarField> nHatf
+        (
+            const volScalarField& alpha1,
+            const volScalarField& alpha2
+        ) const;
+
+        void correctContactAngle
+        (
+            const phase& alpha1,
+            const phase& alpha2,
+            surfaceVectorField::Boundary& nHatb
+        ) const;
+
+        tmp<volScalarField> K(const phase& alpha1, const phase& alpha2) const;
+        tmp<surfaceScalarField> calcStf() const;
+
+public:
+
+    // Constructors
+
+        //- Construct from components
+        multiphaseMixture
+        (
+            const volVectorField& U,
+            const surfaceScalarField& phi,
+            const volScalarField& voidfraction
+        );
+
+
+    //- Destructor
+    virtual ~multiphaseMixture()
+    {}
+
+
+    // Member Functions
+
+        //- Return the phases
+        const PtrDictionary<phase>& phases() const
+        {
+            return phases_;
+        }
+
+        //- Return the velocity
+        const volVectorField& U() const
+        {
+            return U_;
+        }
+
+        //- Return the volumetric flux
+        const surfaceScalarField& phi() const
+        {
+            return phi_;
+        }
+
+        const surfaceScalarField& rhoPhi() const
+        {
+            return rhoPhi_;
+        }
+
+        //- Return the mixture density
+        tmp<volScalarField> rho() const;
+
+        //- Return the mixture density for patch
+        tmp<scalarField> rho(const label patchi) const;
+
+        //- Return the dynamic laminar viscosity
+        tmp<volScalarField> mu() const;
+
+        //- Return the dynamic laminar viscosity for patch
+        tmp<scalarField> mu(const label patchi) const;
+
+        //- Return the face-interpolated dynamic laminar viscosity
+        tmp<surfaceScalarField> muf() const;
+
+        //- Return the kinematic laminar viscosity
+        tmp<volScalarField> nu() const;
+
+        //- Return the laminar viscosity for patch
+        tmp<scalarField> nu(const label patchi) const;
+
+        //- Return the face-interpolated dynamic laminar viscosity
+        tmp<surfaceScalarField> nuf() const;
+
+        //- Return the heat capacity
+        tmp<volScalarField> Cp() const;
+
+        //- Return the thermal conductivity
+        tmp<volScalarField> kf() const;
+
+        //- Return the diffusion coefficient
+        tmp<volScalarField> D() const;
+
+        //- Return the solubility
+        tmp<volScalarField> Cs() const;
+
+        //- Return the diffusion correction term
+        tmp<surfaceScalarField> diffusionCorrection() const;
+
+        tmp<surfaceScalarField> surfaceTensionForce() const
+        {
+            return surfaceTensionForce_;
+        }
+
+        //- Indicator of the proximity of the interface
+        //  Field values are 1 near and 0 away for the interface.
+        tmp<volScalarField> nearInterface() const;
+
+        //- Solve for the mixture phase-fractions
+        void solve();
+
+        //- Correct the mixture properties
+        void correct();
+
+        //- Read base transportProperties dictionary
+        bool read();
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

applications/solvers/cfdemSolverMultiphaseScalar/multiphaseMixture/phase/phase.C

@@ -0,0 +1,107 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2011-2015 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 "phase.H"
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::phase::phase
+(
+    const word& phaseName,
+    const dictionary& phaseDict,
+    const volVectorField& U,
+    const surfaceScalarField& phi
+)
+:
+    volScalarField
+    (
+        IOobject
+        (
+            IOobject::groupName("alpha", phaseName),
+            U.mesh().time().timeName(),
+            U.mesh(),
+            IOobject::MUST_READ,
+            IOobject::AUTO_WRITE
+        ),
+        U.mesh()
+    ),
+    name_(phaseName),
+    phaseDict_(phaseDict),
+    nuModel_
+    (
+        viscosityModel::New
+        (
+            IOobject::groupName("nu", phaseName),
+            phaseDict_,
+            U,
+            phi
+        )
+    ),
+    rho_("rho", dimDensity, phaseDict_),
+    Cp_("Cp", (dimSpecificHeatCapacity), phaseDict_),
+    kf_("kf", (dimPower/dimLength/dimTemperature), phaseDict_),
+    D_("D", dimViscosity, phaseDict_),
+    Cs_("Cs", dimDensity, phaseDict_)
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+Foam::autoPtr<Foam::phase> Foam::phase::clone() const
+{
+    NotImplemented;
+    return autoPtr<phase>(NULL);
+}
+
+
+void Foam::phase::correct()
+{
+    nuModel_->correct();
+}
+
+
+bool Foam::phase::read(const dictionary& phaseDict)
+{
+    phaseDict_ = phaseDict;
+
+    phaseDict_.lookup("Cp") >> Cp_;
+    phaseDict_.lookup("kf") >> kf_;
+    phaseDict_.lookup("D") >> D_;
+    phaseDict_.lookup("Cs") >> Cs_;
+
+    if (nuModel_->read(phaseDict_))
+    {
+        phaseDict_.lookup("rho") >> rho_;
+
+        return true;
+    }
+    else
+    {
+        return false;
+    }
+}
+
+
+// ************************************************************************* //

applications/solvers/cfdemSolverMultiphaseScalar/multiphaseMixture/phase/phase.H

@@ -0,0 +1,190 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2011-2015 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::phase
+
+Description
+    Single incompressible phase derived from the phase-fraction.
+    Used as part of the multiPhaseMixture for interface-capturing multi-phase
+    simulations.
+
+SourceFiles
+    phase.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef phase_H
+#define phase_H
+
+#include "volFields.H"
+#include "dictionaryEntry.H"
+#include "incompressible/viscosityModels/viscosityModel/viscosityModel.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                           Class phase Declaration
+\*---------------------------------------------------------------------------*/
+
+class phase
+:
+    public volScalarField
+{
+    // Private data
+
+        word name_;
+        dictionary phaseDict_;
+        autoPtr<viscosityModel> nuModel_;
+        dimensionedScalar rho_;
+        dimensionedScalar Cp_;
+        dimensionedScalar kf_;
+        dimensionedScalar D_;
+        dimensionedScalar Cs_;
+
+public:
+
+    // Constructors
+
+        //- Construct from components
+        phase
+        (
+            const word& name,
+            const dictionary& phaseDict,
+            const volVectorField& U,
+            const surfaceScalarField& phi
+        );
+
+        //- Return clone
+        autoPtr<phase> clone() const;
+
+        //- Return a pointer to a new phase created on freestore
+        //  from Istream
+        class iNew
+        {
+            const volVectorField& U_;
+            const surfaceScalarField& phi_;
+
+        public:
+
+            iNew
+            (
+                const volVectorField& U,
+                const surfaceScalarField& phi
+            )
+            :
+                U_(U),
+                phi_(phi)
+            {}
+
+            autoPtr<phase> operator()(Istream& is) const
+            {
+                dictionaryEntry ent(dictionary::null, is);
+                return autoPtr<phase>(new phase(ent.keyword(), ent, U_, phi_));
+            }
+        };
+
+
+    // Member Functions
+
+        const word& name() const
+        {
+            return name_;
+        }
+
+        const word& keyword() const
+        {
+            return name();
+        }
+
+        //- Return const-access to phase1 viscosityModel
+        const viscosityModel& nuModel() const
+        {
+            return nuModel_();
+        }
+
+        //- Return the kinematic laminar viscosity
+        tmp<volScalarField> nu() const
+        {
+            return nuModel_->nu();
+        }
+
+        //- Return the laminar viscosity for patch
+        tmp<scalarField> nu(const label patchi) const
+        {
+            return nuModel_->nu(patchi);
+        }
+
+        //- Return const-access to phase1 density
+        const dimensionedScalar& rho() const
+        {
+            return rho_;
+        }
+
+        //- Return const-access to phase1 heat capacity
+        const dimensionedScalar& Cp() const
+        {
+            return Cp_;
+        }
+
+         //- Return const-access to phase1 thermal conductivity
+        const dimensionedScalar& kf() const
+        {
+            return kf_;
+        }
+
+         //- Return const-access to phase1 diffusion coefficient
+        const dimensionedScalar& D() const
+        {
+            return D_;
+        }
+
+        //- Return const-access to phase1 solubility
+        const dimensionedScalar& Cs() const
+        {
+            return Cs_;
+        }
+
+        //- Correct the phase properties
+        void correct();
+
+        //-Inherit read from volScalarField
+        using volScalarField::read;
+
+        //- Read base transportProperties dictionary
+        bool read(const dictionary& phaseDict);
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

applications/solvers/cfdemSolverMultiphaseScalar/pEqn.H

@@ -0,0 +1,73 @@
+{
+    volScalarField rAU("rAU", 1.0/UEqn.A());
+    surfaceScalarField rAUepsf("rAUepsf", fvc::interpolate(rAU*voidfraction));
+    surfaceScalarField rAUepsSqf("rAUepsSqf", fvc::interpolate(rAU*voidfraction*voidfraction));
+    volVectorField Ueps("Ueps", U * voidfraction);
+
+    volVectorField HbyA(constrainHbyA(rAU*UEqn.H(), U, p_rgh));
+
+    surfaceScalarField phiHbyA
+    (
+        "phiHbyA",
+        fvc::flux(HbyA*voidfraction)
+      + fvc::interpolate(voidfraction*rho*rAU)*fvc::ddtCorr(U, phi)
+    );
+
+    adjustPhi(phiHbyA, U, p_rgh);
+
+    if (modelType == "A")
+        rAUepsf = rAUepsSqf;
+
+    surfaceScalarField phig (-ghf*fvc::snGrad(rho)*rAUepsf*mesh.magSf());
+
+    surfaceScalarField phiSt (mixture.surfaceTensionForce()*rAUepsSqf*mesh.magSf());
+
+    surfaceScalarField phiS (fvc::flux(voidfraction*Us*Ksl*rAU));
+
+    phiHbyA += phig + phiSt + phiS;
+
+    // Update the pressure BCs to ensure flux consistency
+    constrainPressure(p_rgh, Ueps, phiHbyA, rAUepsf);
+
+    while (pimple.correctNonOrthogonal())
+    {
+        fvScalarMatrix p_rghEqn
+        (
+            fvm::laplacian(rAUepsf, p_rgh) == particleCloud.ddtVoidfraction() + fvc::div(phiHbyA)
+        );
+
+        p_rghEqn.setReference(pRefCell, getRefCellValue(p_rgh, pRefCell));
+
+        p_rghEqn.solve(mesh.solver(p_rgh.select(pimple.finalInnerIter())));
+
+        if (pimple.finalNonOrthogonalIter())
+        {
+            phi = phiHbyA - p_rghEqn.flux();
+
+            p_rgh.relax();
+
+            if (modelType == "A")
+                U = HbyA + voidfraction*rAU*fvc::reconstruct((phig-p_rghEqn.flux()+phiSt)/rAUepsf) + rAU*Us*Ksl;
+            else
+                U = HbyA + rAU*fvc::reconstruct((phig-p_rghEqn.flux()+phiSt)/rAUepsf) + rAU*Us*Ksl;
+
+            U.correctBoundaryConditions();
+            fvOptions.correct(U);
+        }
+    }
+
+    #include "continuityErrs.H"
+
+    p == p_rgh + rho*gh;
+
+    if (p_rgh.needReference())
+    {
+        p += dimensionedScalar
+        (
+            "p",
+            p.dimensions(),
+            pRefValue - getRefCellValue(p, pRefCell)
+        );
+        p_rgh = p - rho*gh;
+    }
+}

applications/solvers/cfdemSolverPiso/UEqn.H

@@ -1,8 +1,11 @@
+particleCloud.otherForces(fOther);
+
 fvVectorMatrix UEqn
 (
     fvm::ddt(voidfraction,U) - fvm::Sp(fvc::ddt(voidfraction),U)
   + fvm::div(phi,U) - fvm::Sp(fvc::div(phi),U)
   + particleCloud.divVoidfractionTau(U, voidfraction)
+  - fOther/rho
   ==
     fvOptions(U)
   - fvm::Sp(Ksl/rho,U)

applications/solvers/cfdemSolverPiso/createFields.H

@@ -46,6 +46,21 @@
         //dimensionedScalar("0", dimensionSet(1, -3, -1, 0, 0), 1.0)
     );
 
+    Info<< "\nCreating body force field\n" << endl;
+    volVectorField fOther
+    (
+        IOobject
+        (
+            "fOther",
+            runTime.timeName(),
+            mesh,
+            IOobject::NO_READ,
+            IOobject::NO_WRITE
+        ),
+        mesh,
+        dimensionedVector("zero", dimensionSet(1,-2,-2,0,0,0,0), vector::zero)
+    );
+
     Info<< "\nReading voidfraction field voidfraction = (Vgas/Vparticle)\n" << endl;
     volScalarField voidfraction
     (

applications/solvers/cfdemSolverPisoFreeStreaming/Make/files

@@ -0,0 +1,3 @@
+cfdemSolverPisoFreeStreaming.C
+
+EXE=$(CFDEM_APP_DIR)/cfdemSolverPisoFreeStreaming

applications/solvers/cfdemSolverPisoFreeStreaming/Make/options

@@ -0,0 +1,28 @@
+include $(CFDEM_ADD_LIBS_DIR)/additionalLibs
+
+EXE_INC = \
+    -I$(CFDEM_OFVERSION_DIR) \
+    -I$(LIB_SRC)/finiteVolume/lnInclude \
+    -I$(LIB_SRC)/meshTools/lnInclude \
+    -I$(LIB_SRC)/TurbulenceModels/turbulenceModels/lnInclude \
+    -I$(LIB_SRC)/TurbulenceModels/incompressible/lnInclude \
+    -I$(LIB_SRC)/transportModels \
+    -I$(LIB_SRC)/transportModels/incompressible/singlePhaseTransportModel \
+    -I$(FOAM_SOLVERS)/incompressible/pisoFoam \
+    -I$(CFDEM_SRC_DIR)/lagrangian/cfdemParticle/lnInclude \
+    -I$(CFDEM_SRC_DIR)/lagrangian/cfdemParticle/cfdTools \
+    -I$(CFDEM_SRC_DIR)/lagrangian/cfdemParticle/derived/cfdemCloudRec \
+    -I$(LIB_SRC)/sampling/lnInclude
+
+EXE_LIBS = \
+    -L$(CFDEM_LIB_DIR)\
+    -lturbulenceModels \
+    -lincompressibleTurbulenceModels \
+    -lincompressibleTransportModels \
+    -lfiniteVolume \
+    -lmeshTools \
+    -lfvOptions \
+    -lsampling \
+    -l$(CFDEM_LIB_NAME) \
+     $(CFDEM_ADD_LIB_PATHS) \
+     $(CFDEM_ADD_LIBS)

applications/solvers/cfdemSolverPisoFreeStreaming/cfdemSolverPisoFreeStreaming.C

@@ -0,0 +1,126 @@
+/*---------------------------------------------------------------------------*\
+    CFDEMcoupling - Open Source CFD-DEM coupling
+
+    CFDEMcoupling is part of the CFDEMproject
+    www.cfdem.com
+                                Christoph Goniva, christoph.goniva@cfdem.com
+                                Copyright (C) 1991-2009 OpenCFD Ltd.
+                                Copyright (C) 2009-2012 JKU, Linz
+                                Copyright (C) 2012-     DCS Computing GmbH,Linz
+-------------------------------------------------------------------------------
+License
+    This file is part of CFDEMcoupling.
+
+    CFDEMcoupling 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.
+
+    CFDEMcoupling 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 CFDEMcoupling.  If not, see <http://www.gnu.org/licenses/>.
+
+Application
+    cfdemSolverPisoFreeStreaming
+
+Description
+    Transient solver for incompressible flow.
+    Turbulence modelling is generic, i.e. laminar, RAS or LES may be selected.
+    The code is an evolution of the solver pisoFoam in OpenFOAM(R) 1.6,
+    where additional functionality for CFD-DEM coupling is added.
+    the particles follow the fluid velocity
+\*---------------------------------------------------------------------------*/
+
+#include "fvCFD.H"
+#include "singlePhaseTransportModel.H"
+#include "turbulentTransportModel.H"
+#include "pisoControl.H"
+#include "fvOptions.H"
+
+#include "cfdemCloudRec.H"
+
+#include "cfdemCloud.H"
+#include "implicitCouple.H"
+#include "clockModel.H"
+#include "smoothingModel.H"
+#include "forceModel.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+int main(int argc, char *argv[])
+{
+    #include "postProcess.H"
+    #include "setRootCase.H"
+    #include "createTime.H"
+    #include "createMesh.H"
+    #include "createControl.H"
+    #include "createFields.H"
+    #include "createFvOptions.H"
+    #include "initContinuityErrs.H"
+
+    cfdemCloudRec<cfdemCloud> particleCloud(mesh);
+    // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+    Info<< "\nStarting time loop\n" << endl;
+    while (runTime.loop())
+    {
+
+        particleCloud.clockM().start(1,"Global");
+
+        Info<< "Time = " << runTime.timeName() << nl << endl;
+
+        #include "CourantNo.H"
+
+        // do particle stuff
+        particleCloud.clockM().start(2,"Coupling");
+
+        particleCloud.evolve(voidfraction,Us,U);
+
+
+        particleCloud.clockM().stop("Coupling");
+
+        particleCloud.clockM().start(26,"Flow");
+
+        if(particleCloud.solveFlow())
+        {
+            // Pressure-velocity PISO corrector
+            {
+                // Momentum predictor
+                 #include "UEqn.H"
+
+                // --- PISO loop
+
+                while (piso.correct())
+                {
+                    #include "pEqn.H"
+                }
+            }
+
+            laminarTransport.correct();
+            turbulence->correct();
+        }
+        else
+        {
+            Info << "skipping flow solution." << endl;
+        }
+
+        runTime.write();
+
+        Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
+            << "  ClockTime = " << runTime.elapsedClockTime() << " s"
+            << nl << endl;
+
+        particleCloud.clockM().stop("Flow");
+        particleCloud.clockM().stop("Global");
+    }
+
+    Info<< "End\n" << endl;
+
+    return 0;
+}
+
+
+// ************************************************************************* //

applications/solvers/cfdemSolverPisoFreeStreaming/createFields.H

@@ -0,0 +1,110 @@
+    Info<< "Reading field p\n" << endl;
+    volScalarField p
+    (
+        IOobject
+        (
+            "p",
+            runTime.timeName(),
+            mesh,
+            IOobject::MUST_READ,
+            IOobject::AUTO_WRITE
+        ),
+        mesh
+    );
+
+    Info<< "Reading physical velocity field U" << endl;
+    Info<< "Note: only if voidfraction at boundary is 1, U is superficial velocity!!!\n" << endl;
+    volVectorField U
+    (
+        IOobject
+        (
+            "U",
+            runTime.timeName(),
+            mesh,
+            IOobject::MUST_READ,
+            IOobject::AUTO_WRITE
+        ),
+        mesh
+    );
+
+//===============================
+// particle interaction modelling
+//===============================
+
+    Info<< "\nReading voidfraction field voidfraction = (Vgas/Vparticle)\n" << endl;
+    volScalarField voidfraction
+    (
+        IOobject
+        (
+            "voidfraction",
+            runTime.timeName(),
+            mesh,
+            IOobject::MUST_READ,
+            IOobject::AUTO_WRITE
+        ),
+        mesh
+    );
+
+    Info<< "\nCreating density field rho\n" << endl;
+    volScalarField rho
+    (
+        IOobject
+        (
+            "rho",
+            runTime.timeName(),
+            mesh,
+            IOobject::MUST_READ,
+            IOobject::AUTO_WRITE
+        ),
+        mesh
+    );
+
+    Info<< "Reading particle velocity field Us\n" << endl;
+    volVectorField Us
+    (
+        IOobject
+        (
+            "Us",
+            runTime.timeName(),
+            mesh,
+            IOobject::MUST_READ,
+            IOobject::AUTO_WRITE
+        ),
+        mesh
+    );
+
+//===============================
+
+//#   include "createPhi.H"
+#ifndef createPhi_H
+#define createPhi_H
+Info<< "Reading/calculating face flux field phi\n" << endl;
+surfaceScalarField phi
+(
+    IOobject
+    (
+        "phi",
+        runTime.timeName(),
+         mesh,
+        IOobject::READ_IF_PRESENT,
+        IOobject::AUTO_WRITE
+    ),
+    linearInterpolate(U) & mesh.Sf()
+);
+#endif
+
+
+
+    label pRefCell = 0;
+    scalar pRefValue = 0.0;
+    setRefCell(p, mesh.solutionDict().subDict("PISO"), pRefCell, pRefValue);
+
+
+    singlePhaseTransportModel laminarTransport(U, phi);
+
+    autoPtr<incompressible::turbulenceModel> turbulence
+    (
+        incompressible::turbulenceModel::New(U, phi, laminarTransport)
+    );
+
+#include "createMRF.H"

applications/solvers/cfdemSolverPisoScalar/createFields.H

@@ -46,6 +46,21 @@
         //dimensionedScalar("0", dimensionSet(0, 0, -1, 0, 0), 1.0)
     );
 
+    Info<< "\nCreating body force field\n" << endl;
+    volVectorField fOther
+    (
+        IOobject
+        (
+            "fOther",
+            runTime.timeName(),
+            mesh,
+            IOobject::NO_READ,
+            IOobject::NO_WRITE
+        ),
+        mesh,
+        dimensionedVector("zero", dimensionSet(1,-2,-2,0,0,0,0), vector::zero)
+    );
+
     Info<< "\nReading voidfraction field voidfraction = (Vgas/Vparticle)\n" << endl;
     volScalarField voidfraction
     (

applications/solvers/cfdemSolverRhoPimple/EEqn.H

@@ -6,7 +6,18 @@
     particleCloud.energyContributions(Qsource);
     particleCloud.energyCoefficients(QCoeff);
 
-    addSource = fvc::ddt(rhoeps, K) + fvc::div(phi, K)
+    Cpv = he.name() == "e" ? thermo.Cv() : thermo.Cp();
+
+// For implict T terms in the energy/enthalpy transport equation, use
+// (he_n+1 - he_n) / (T_n+1 - T_n) = Cpv to eliminate T_n+1 with he_n+1.
+// This formula is valid for ideal gases with e=e(T) and h=h(T). For
+// incompressible fluids, e=e(T) holds, too, but enthalpy would need correction
+// terms accounting for pressure variations.
+
+    fvScalarMatrix EEqn
+    (
+        fvm::ddt(rhoeps, he) + fvm::div(phi, he)
+      + fvc::ddt(rhoeps, K) + fvc::div(phi, K)
       + (
             he.name() == "e"
           ? fvc::div
@@ -16,25 +27,14 @@
                 "div(phiv,p)"
             )
           : -dpdt
-        );
-
-    Cpv = he.name() == "e" ? thermo.Cv() : thermo.Cp();
-
-    // correct source for the thermodynamic reference temperature
-    dimensionedScalar Tref("Tref", dimTemperature, T[0]-he[0]/(Cpv[0]+SMALL));
-    Qsource += QCoeff*Tref;
-
-    fvScalarMatrix EEqn
-    (
-        fvm::ddt(rhoeps, he) + fvm::div(phi, he)
-      + addSource
-     // net heat transfer from particles to fluid
+        )
       - Qsource
+      - QCoeff*T
       - fvm::Sp(QCoeff/Cpv, he)
-     // thermal conduction of the fluid with effective conductivity
+      + QCoeff/Cpv*he
+      - fvc::laplacian(voidfraction*thCond,T)
       - fvm::laplacian(voidfraction*thCond/Cpv,he)
-      // + particle-fluid energy transfer due to work
-      // + fluid energy dissipation due to shearing
+      + fvc::laplacian(voidfraction*thCond/Cpv,he)
      ==
         fvOptions(rho, he)
     );
@@ -50,7 +50,7 @@
 
     thermo.correct();
 
-    Info<< "T max/min : " << max(T).value() << " " << min(T).value() << endl;
+    Info<< "T max/min/ave : " << max(T).value() << " " << min(T).value() << " " << average(T).value() << endl;
 
     particleCloud.clockM().start(31,"energySolve");
     particleCloud.solve();

applications/solvers/cfdemSolverRhoPimple/cfdemSolverRhoPimple.C

@@ -32,6 +32,9 @@ Description
 #include "turbulentFluidThermoModel.H"
 #include "bound.H"
 #include "pimpleControl.H"
+#if OPENFOAM_VERSION_MAJOR >= 5
+#include "pressureControl.H"
+#endif
 #include "fvOptions.H"
 #include "localEulerDdtScheme.H"
 #include "fvcSmooth.H"
@@ -69,16 +72,19 @@ int main(int argc, char *argv[])
     #include "checkModelType.H"
 
     turbulence->validate();
-    //#include "compressibleCourantNo.H"
-    //#include "setInitialDeltaT.H"
+
+    #include "compressibleCourantNo.H"
+    #include "setInitialDeltaT.H"
 
     // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
     Info<< "\nStarting time loop\n" << endl;
+    bool firstStep = true;
 
     while (runTime.run())
     {
         #include "readTimeControls.H"
+
         #include "compressibleCourantNo.H"
         #include "setDeltaT.H"
 
@@ -90,6 +96,7 @@ int main(int argc, char *argv[])
 
         // do particle stuff
         particleCloud.clockM().start(2,"Coupling");
+
         bool hasEvolved = particleCloud.evolve(voidfraction,Us,U);
 
         if(hasEvolved && smoothenForces)
@@ -109,23 +116,35 @@ int main(int argc, char *argv[])
         Info << "TotalForceImp: " << fImpTotal << endl;
 
         #include "solverDebugInfo.H"
+
         particleCloud.clockM().stop("Coupling");
 
         particleCloud.clockM().start(26,"Flow");
 
 #if OPENFOAM_VERSION_MAJOR < 6
         if (pimple.nCorrPIMPLE() <= 1)
-#else
-        if (pimple.nCorrPimple() <= 1)
-#endif
         {
             #include "rhoEqn.H"
         }
+        rhoeps = rho*voidfraction;
+#endif
 
-        volScalarField rhoeps("rhoeps",rho*voidfraction);
         // --- Pressure-velocity PIMPLE corrector loop
         while (pimple.loop())
         {
+#if OPENFOAM_VERSION_MAJOR >= 6
+            if (pimple.firstIter())
+            {
+                #include "rhoEqn.H"
+                if (firstStep)
+                {
+                    rhoeps.oldTime() = rho.oldTime()*voidfraction.oldTime();
+                    firstStep = false;
+                }
+                rhoeps = rho*voidfraction;
+            }
+#endif
+
             #include "UEqn.H"
             #include "EEqn.H"
 
@@ -134,7 +153,6 @@ int main(int argc, char *argv[])
             {
                 // besides this pEqn, OF offers a "pimple consistent"-option
                 #include "pEqn.H"
-                rhoeps=rho*voidfraction;
             }
 
             if (pimple.turbCorr())
@@ -149,7 +167,6 @@ int main(int argc, char *argv[])
 
         runTime.write();
 
-
         Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
             << "  ClockTime = " << runTime.elapsedClockTime() << " s"
             << nl << endl;

applications/solvers/cfdemSolverRhoPimple/createFields.H

@@ -51,20 +51,49 @@ Info<< "Reading thermophysical properties\n" << endl;
         mesh
     );
 
-    volScalarField addSource
+    volScalarField rhoeps("rhoeps", rho*voidfraction);
+    rhoeps.oldTime(); // switch on saving old time
+
+    Info<< "Reading/calculating face flux field phi\n" << endl;
+    surfaceScalarField phi
     (
         IOobject
         (
-           "addSource",
-           runTime.timeName(),
-           mesh,
-           IOobject::READ_IF_PRESENT,
-           IOobject::AUTO_WRITE
-        ),
-        mesh,
-        dimensionedScalar("zero", dimensionSet(1,-1,-3,0,0,0,0), 0.0)
+            "phi",
+            runTime.timeName(),
+            mesh,
+            IOobject::READ_IF_PRESENT,
+            IOobject::AUTO_WRITE
+         ),
+         linearInterpolate(rho*U*voidfraction) & mesh.Sf()
     );
 
+#if OPENFOAM_VERSION_MAJOR < 5
+    dimensionedScalar rhoMax
+    (
+        dimensionedScalar::lookupOrDefault
+        (
+            "rhoMax",
+            pimple.dict(),
+            dimDensity,
+            GREAT
+        )
+    );
+
+    dimensionedScalar rhoMin
+    (
+        dimensionedScalar::lookupOrDefault
+        (
+            "rhoMin",
+            pimple.dict(),
+            dimDensity,
+            0
+        )
+    );
+#else
+    pressureControl pressureControl(p, rho, pimple.dict(), false);
+#endif
+
     Info<< "\nCreating fluid-particle heat flux field\n" << endl;
     volScalarField Qsource
     (
@@ -141,42 +170,6 @@ Info<< "Reading thermophysical properties\n" << endl;
         dimensionedVector("zero", dimensionSet(1,-2,-2,0,0,0,0), vector::zero)
     );
 
-    Info<< "Reading/calculating face flux field phi\n" << endl;
-    surfaceScalarField phi
-    (
-        IOobject
-        (
-            "phi",
-            runTime.timeName(),
-            mesh,
-            IOobject::READ_IF_PRESENT,
-            IOobject::AUTO_WRITE
-         ),
-         linearInterpolate(rho*U*voidfraction) & mesh.Sf()
-    );
-
-    dimensionedScalar rhoMax
-    (
-        dimensionedScalar::lookupOrDefault
-        (
-            "rhoMax",
-            pimple.dict(),
-            dimDensity,
-            GREAT
-        )
-    );
-
-    dimensionedScalar rhoMin
-    (
-        dimensionedScalar::lookupOrDefault
-        (
-            "rhoMin",
-            pimple.dict(),
-            dimDensity,
-            0
-        )
-    );
-
     bool smoothenForces
     (
         pimple.dict().lookupOrDefault<bool>

applications/solvers/cfdemSolverRhoPimple/pEqn.H

@@ -1,14 +1,19 @@
 rho = thermo.rho();
+#if OPENFOAM_VERSION_MAJOR < 5
 rho = max(rho, rhoMin);
 rho = min(rho, rhoMax);
 rho.relax();
+rhoeps = rho*voidfraction;
+#else
+rhoeps = rho*voidfraction;
+
+// Thermodynamic density needs to be updated by psi*d(p) after the
+// pressure solution
+const volScalarField psip0(psi*p);
+#endif
 
 volScalarField rAU(1.0/UEqn.A());
-surfaceScalarField rhorAUf("rhorAUf", fvc::interpolate(rhoeps*rAU));
-if (modelType=="A")
-{
-    rhorAUf *= fvc::interpolate(voidfraction);
-}
+surfaceScalarField rhorAUf("rhorAUf", (modelType=="A")?fvc::interpolate(voidfraction*rhoeps*rAU):fvc::interpolate(rhoeps*rAU));
 volVectorField HbyA(constrainHbyA(rAU*UEqn.H(), U, p));
 
 surfaceScalarField phiUs("phiUs", fvc::interpolate(rhoeps*rAU*Ksl*Us)& mesh.Sf());
@@ -18,30 +23,40 @@ if (pimple.nCorrPISO() <= 1)
     tUEqn.clear();
 }
 
+surfaceScalarField phiHbyA
+(
+    "phiHbyA",
+    fvc::interpolate(rhoeps)*fvc::flux(HbyA)
+  + rhorAUf*fvc::ddtCorr(rhoeps, U, phi)
+);
+
 if (pimple.transonic())
 {
 //  transonic version not implemented yet
 }
 else
 {
-    surfaceScalarField phiHbyA
-    (
-        "phiHbyA",
-        (
-            fvc::flux(rhoeps*HbyA)
-    //      + rhorAUf*fvc::ddtCorr(rho, U, phi)
-        )
-    );
-
     // flux without pressure gradient contribution
     phi = phiHbyA + phiUs;
 
     // Update the pressure BCs to ensure flux consistency
     constrainPressure(p, rhoeps, U, phi, rhorAUf);
 
+#if OPENFOAM_VERSION_MAJOR >= 5
+    fvScalarMatrix pDDtEqn
+    (
+        fvc::ddt(rhoeps)
+      + psi*voidfraction*correction(fvm::ddt(p))
+      + fvc::div(phi)
+     ==
+        fvOptions(psi, p, rho.name())
+    );
+#endif
+
     while (pimple.correctNonOrthogonal())
     {
         // Pressure corrector
+#if OPENFOAM_VERSION_MAJOR < 5
         fvScalarMatrix pEqn
         (
             fvm::ddt(psi*voidfraction, p)
@@ -50,6 +65,9 @@ else
           ==
             fvOptions(psi, p, rho.name())
         );
+#else
+        fvScalarMatrix pEqn(pDDtEqn - fvm::laplacian(rhorAUf, p));
+#endif
 
         pEqn.solve(mesh.solver(p.select(pimple.finalInnerIter())));
 
@@ -60,19 +78,18 @@ else
     }
 }
 
+// Thermodynamic density update
+#if OPENFOAM_VERSION_MAJOR >= 5
+thermo.correctRho(psi*p - psip0);
+#endif
+
 #include "rhoEqn.H"
 #include "compressibleContinuityErrsPU.H"
 
 // Explicitly relax pressure for momentum corrector
 p.relax();
-
-// Recalculate density from the relaxed pressure
-rho = thermo.rho();
-rho = max(rho, rhoMin);
-rho = min(rho, rhoMax);
-rho.relax();
-Info<< "rho max/min : " << max(rho).value()
-    << " " << min(rho).value() << endl;
+Info<< "p max/min/ave : " << max(p).value()
+    << " " << min(p).value() << " " << average(p).value() << endl;
 
 if (modelType=="A")
 {
@@ -86,6 +103,24 @@ U.correctBoundaryConditions();
 fvOptions.correct(U);
 K = 0.5*magSqr(U);
 
+// Recalculate density from the relaxed pressure
+#if OPENFOAM_VERSION_MAJOR >= 5
+if (pressureControl.limit(p))
+{
+    p.correctBoundaryConditions();
+}
+rho = thermo.rho();
+#else
+rho = thermo.rho();
+rho = max(rho, rhoMin);
+rho = min(rho, rhoMax);
+rho.relax();
+#endif
+rhoeps = rho*voidfraction;
+
+Info<< "rho max/min/ave : " << max(rho).value()
+    << " " << min(rho).value() << " " << average(rho).value() << endl;
+
 if (thermo.dpdt())
 {
     dpdt = fvc::ddt(voidfraction,p);

applications/solvers/cfdemSolverRhoPimpleChem/EEqn.H

@@ -8,9 +8,11 @@ particleCloud.energyCoefficients(QCoeff);
 
 Cpv = he.name() == "e" ? thermo.Cv() : thermo.Cp();
 
-// correct source for the thermodynamic reference temperature
-// dimensionedScalar Tref("Tref", dimTemperature, T[0]-he[0]/(Cpv[0]+SMALL));
-// Qsource += QCoeff*Tref;
+// For implict T terms in the energy/enthalpy transport equation, use
+// (he_n+1 - he_n) / (T_n+1 - T_n) = Cpv to eliminate T_n+1 with he_n+1.
+// This formula is valid for ideal gases with e=e(T) and h=h(T). For
+// incompressible fluids, e=e(T) holds, too, but enthalpy would need correction
+// terms accounting for pressure variations.
 
 fvScalarMatrix EEqn
 (
@@ -26,13 +28,13 @@ fvScalarMatrix EEqn
            )
           : -dpdt
        )
-    // net heat transfer from particles to fluid
       - Qsource
+      - QCoeff*T
       - fvm::Sp(QCoeff/Cpv, he)
-    // thermal conduction of the fluid with effective conductivity
+      + QCoeff/Cpv*he
+      - fvc::laplacian(voidfraction*thCond,T)
       - fvm::laplacian(voidfraction*thCond/Cpv,he)
-    // + particle-fluid energy transfer due to work
-    // + fluid energy dissipation due to shearing
+      + fvc::laplacian(voidfraction*thCond/Cpv,he)
      ==
 //      + combustion->Sh()
        fvOptions(rho, he)

applications/solvers/cfdemSolverRhoPimpleChem/Make/options

@@ -6,7 +6,6 @@ PFLAGS+= -Dcompre
 
 EXE_INC = \
      $(PFLAGS) \
-    -I../. \
     -I$(CFDEM_OFVERSION_DIR) \
     -I$(LIB_SRC)/finiteVolume/cfdTools \
     -I$(LIB_SRC)/finiteVolume/lnInclude \
@@ -27,7 +26,6 @@ EXE_INC = \
     -I$(LIB_SRC)/regionModels/surfaceFilmModels/lnInclude \
     -I$(LIB_SRC)/ODE/lnInclude \
     -I$(LIB_SRC)/combustionModels/lnInclude \
-    -I$(FOAM_SOLVERS)/combustion/reactingFoam \
     -Wno-deprecated-copy
 
 

applications/solvers/cfdemSolverRhoPimpleChem/cfdemSolverRhoPimpleChem.C

@@ -57,6 +57,8 @@ Description
 
 int main(int argc, char *argv[])
 {
+    #include "postProcess.H"
+
     #include "setRootCase.H"
     #include "createTime.H"
     #include "createMesh.H"
@@ -64,9 +66,10 @@ int main(int argc, char *argv[])
     #include "createTimeControls.H"
     #include "createRDeltaT.H"
 
-    #include "createFields.H"
-    #include "createFvOptions.H"
     #include "initContinuityErrs.H"
+    #include "createFields.H"
+    #include "createFieldRefs.H"
+    #include "createFvOptions.H"
 
     // create cfdemCloud
     #include "readGravitationalAcceleration.H"

applications/solvers/cfdemSolverRhoPimpleChem/createFields.H

@@ -24,7 +24,7 @@
     volScalarField W(thermo.W());
 #endif
 
-    bool propagateInertSpecie = true;
+    Switch propagateInertSpecie(true);
 
     const word inertSpecie(thermo.lookup("inertSpecie"));
 
@@ -40,9 +40,9 @@
             << exit(FatalError);
     }
 
+    Info<< "inert will be bounded in [" << inertLowerBound << "," << inertUpperBound << "]" << endl;
+
     volScalarField& p = thermo.p();
-    const volScalarField& T = thermo.T();
-    const volScalarField& psi = thermo.psi();
 
     multivariateSurfaceInterpolationScheme<scalar>::fieldTable fields;
 

applications/solvers/rcfdemSolverBase/createFields.H

@@ -37,12 +37,20 @@
             "URec",
             runTime.timeName(),
             mesh,
-            IOobject::MUST_READ,
-            IOobject::AUTO_WRITE
+            IOobject::READ_IF_PRESENT,
+            IOobject::NO_WRITE
         ),
-        mesh
+        mesh,
+        dimensionedVector("URec", dimensionSet(0, 1, -1, 0, 0), vector::zero)
     );
 
+    Switch updateURec(false);
+    if (URec.headerOk())
+    {
+        updateURec = true;
+        URec.writeOpt() = IOobject::AUTO_WRITE;
+    }
+
     volScalarField voidfractionRec
     (
         IOobject
@@ -50,12 +58,20 @@
             "voidfractionRec",
             runTime.timeName(),
             mesh,
-            IOobject::MUST_READ,
-            IOobject::AUTO_WRITE
+            IOobject::READ_IF_PRESENT,
+            IOobject::NO_WRITE
         ),
-        mesh
+        mesh,
+        dimensionedScalar("voidfractionRec", dimensionSet(0, 0, 0, 0, 0), 1.0)
     );
 
+    Switch updateVoidfractionRec(false);
+    if (voidfractionRec.headerOk())
+    {
+        updateVoidfractionRec = true;
+        voidfractionRec.writeOpt() = IOobject::AUTO_WRITE;
+    }
+
     volVectorField UsRec
     (
         IOobject
@@ -63,12 +79,20 @@
             "UsRec",
             runTime.timeName(),
             mesh,
-            IOobject::MUST_READ,
+            IOobject::READ_IF_PRESENT,
             IOobject::AUTO_WRITE
         ),
-        mesh
+        mesh,
+        dimensionedVector("URec", dimensionSet(0, 1, -1, 0, 0), vector::zero)
     );
 
+    Switch updateUsRec(false);
+    if (UsRec.headerOk())
+    {
+        updateUsRec = true;
+        UsRec.writeOpt() = IOobject::AUTO_WRITE;
+    }
+
     // calculated fields
     Info << "\nCreating fields subject to calculation\n" << endl;
     volScalarField voidfraction
@@ -78,7 +102,7 @@
             "voidfraction",
             runTime.timeName(),
             mesh,
-            IOobject::NO_READ,
+            IOobject::READ_IF_PRESENT,
             IOobject::AUTO_WRITE
         ),
         voidfractionRec
@@ -91,7 +115,7 @@
             "Us",
             runTime.timeName(),
             mesh,
-            IOobject::NO_READ,
+            IOobject::READ_IF_PRESENT,
             IOobject::AUTO_WRITE
         ),
         UsRec
@@ -111,11 +135,18 @@
             runTime.timeName(),
             mesh,
             IOobject::READ_IF_PRESENT,
-            IOobject::AUTO_WRITE
+            IOobject::NO_WRITE
         ),
         linearInterpolate(URec*voidfractionRec) & mesh.Sf()
     );
-    phiRec.write();
+
+    Switch updatePhiRec(false);
+    if (phiRec.headerOk())
+    {
+        updatePhiRec = true;
+        phiRec.writeOpt() = IOobject::AUTO_WRITE;
+        phiRec.write();
+    }
 
     singlePhaseTransportModel laminarTransport(URec, phiRec);
 
@@ -123,3 +154,40 @@
     (
         incompressible::turbulenceModel::New(URec, phiRec, laminarTransport)
     );
+
+    IOdictionary recDict
+    (
+        IOobject
+        (
+            "recProperties",
+            runTime.constant(),
+            mesh,
+            IOobject::MUST_READ,
+            IOobject::NO_WRITE
+        )
+    );
+
+    word voidfractionFieldName(recDict.lookupOrDefault<word>("voidfractionFieldName","voidfraction"));
+    word UFieldName(recDict.lookupOrDefault<word>("UFieldName","U"));
+    word UsFieldName(recDict.lookupOrDefault<word>("UsFieldName","Us"));
+    word fluxFieldName(recDict.lookupOrDefault<word>("fluxFieldName","phi"));
+
+
+   // place to put weight functions
+    IOdictionary weightDict
+    (
+        IOobject
+        (
+            "weightDict",
+            runTime.constant(),
+            mesh,
+            IOobject::READ_IF_PRESENT,
+            IOobject::NO_WRITE
+        )
+    );
+    if (!weightDict.headerOk())
+    {
+        weightDict.add("weights",scalarList(1,1.0));
+    }
+    scalarList weights(weightDict.lookup("weights"));
+    Info << "database initial weights: " << weights << endl;

applications/solvers/rcfdemSolverBase/rcfdemSolverBase.C

@@ -43,6 +43,7 @@ Rules
 #include "cfdemCloudRec.H"
 #include "recBase.H"
 #include "recModel.H"
+#include "recPath.H"
 
 #include "cfdemCloud.H"
 #include "clockModel.H"
@@ -57,7 +58,8 @@ int main(int argc, char *argv[])
     #include "createMesh.H"
     #include "createControl.H"
     #include "createFields.H"
-    #include "createFvOptions.H"
+
+    #include "readGravitationalAcceleration.H"
 
     cfdemCloudRec<cfdemCloud> particleCloud(mesh);
     recBase recurrenceBase(mesh);
@@ -67,8 +69,8 @@ int main(int argc, char *argv[])
     Info << "\nCalculating particle trajectories based on recurrence statistics\n" << endl;
 
     label recTimeIndex = 0;
-    scalar recTimeStep = recurrenceBase.recM().recTimeStep();
-    scalar startTime = runTime.startTime().value();
+    label stepCounter = 0;
+    label recTimeStep2CFDTimeStep = recurrenceBase.recM().recTimeStep2CFDTimeStep();
 
     while (runTime.run())
     {
@@ -85,12 +87,16 @@ int main(int argc, char *argv[])
 
         particleCloud.clockM().stop("Coupling");
 
+        stepCounter++;
 
-        if ( runTime.timeOutputValue() - startTime - (recTimeIndex+1)*recTimeStep + 1.0e-5 > 0.0 )
+        if (stepCounter == recTimeStep2CFDTimeStep)
         {
+            Info << "updating recurrence fields at time " << runTime.timeName() << "with recTimeIndex = " << recTimeIndex << nl << endl;
             recurrenceBase.updateRecFields();
-            #include "readFields.H"
+            #include "updateFields.H"
             recTimeIndex++;
+            stepCounter = 0;
+            recTimeStep2CFDTimeStep = recurrenceBase.recM().recTimeStep2CFDTimeStep();
         }
 
         particleCloud.clockM().start(27,"Output");
@@ -102,7 +108,6 @@ int main(int argc, char *argv[])
         Info << "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
              << "  ClockTime = " << runTime.elapsedClockTime() << " s"
              << nl << endl;
-
     }
 
     Info << "End\n" << endl;

applications/solvers/rcfdemSolverBase/updateFields.H

@@ -0,0 +1,29 @@
+scalarList wList(weightDict.lookupOrDefault("weights",scalarList(1,1.0)));
+
+recurrenceBase.recP().updateIntervalWeights(wList);
+
+if(recurrenceBase.recM().endOfPath())
+{
+    recurrenceBase.extendPath();
+}
+
+// update fields where necessary
+if (updateVoidfractionRec)
+{
+    recurrenceBase.recM().exportVolScalarField(voidfractionFieldName,voidfractionRec);
+}
+
+if (updateURec)
+{
+    recurrenceBase.recM().exportVolVectorField(UFieldName,URec);
+}
+
+if (updateUsRec)
+{
+    recurrenceBase.recM().exportVolVectorField(UsFieldName,UsRec);
+}
+
+if (updatePhiRec)
+{
+    recurrenceBase.recM().exportSurfaceScalarField(fluxFieldName,phiRec);
+}

applications/solvers/rcfdemSolverCoupledHeattransfer/TEqImp.H

@@ -10,6 +10,7 @@
     // main contribution due to gas expansion, not due to transport of kinetic energy
     // fvc::ddt(rhoeps, K) + fvc::div(phiRec, K)
 
+    // assuming constant Cv such that e = Cv * T
     fvScalarMatrix TEqn =
     (
         fvm::ddt(rhoeps, T)

applications/solvers/rcfdemSolverCoupledHeattransfer/rcfdemSolverCoupledHeattransfer.C

@@ -67,8 +67,8 @@ int main(int argc, char *argv[])
     Info << "\nCalculating particle trajectories based on recurrence statistics\n" << endl;
 
     label recTimeIndex = 0;
-    scalar recTimeStep = recurrenceBase.recM().recTimeStep();
-    scalar startTime = runTime.startTime().value();
+    label stepCounter = 0;
+    label recTimeStep2CFDTimeStep = recurrenceBase.recM().recTimeStep2CFDTimeStep();
 
     // control coupling behavior in case of substepping
     // assumes constant timestep size
@@ -101,12 +101,16 @@ int main(int argc, char *argv[])
         #include "TEqImp.H"
         particleCloud.clockM().stop("Flow");
 
+        stepCounter++;
+
         particleCloud.clockM().start(32,"ReadFields");
-        if ( runTime.timeOutputValue() - startTime - (recTimeIndex+1)*recTimeStep + 1.0e-5 > 0.0 )
+        if (stepCounter == recTimeStep2CFDTimeStep)
         {
             recurrenceBase.updateRecFields();
             #include "updateFields.H"
             recTimeIndex++;
+            stepCounter = 0;
+            recTimeStep2CFDTimeStep = recurrenceBase.recM().recTimeStep2CFDTimeStep();
         }
         particleCloud.clockM().stop("ReadFields");
 

applications/solvers/rcfdemSolverForcedTracers/Make/files

@@ -0,0 +1,3 @@
+rcfdemSolverForcedTracers.C
+
+EXE=$(CFDEM_APP_DIR)/rcfdemSolverForcedTracers

applications/solvers/rcfdemSolverForcedTracers/Make/options

@@ -0,0 +1,27 @@
+include $(CFDEM_ADD_LIBS_DIR)/additionalLibs
+
+EXE_INC = \
+    -I$(CFDEM_OFVERSION_DIR) \
+    -I$(LIB_SRC)/finiteVolume/lnInclude \
+    -I$(LIB_SRC)/meshTools/lnInclude \
+    -I$(LIB_SRC)/TurbulenceModels/turbulenceModels/lnInclude \
+    -I$(LIB_SRC)/TurbulenceModels/incompressible/lnInclude \
+    -I$(LIB_SRC)/transportModels \
+    -I$(LIB_SRC)/transportModels/incompressible/singlePhaseTransportModel \
+    -I$(CFDEM_SRC_DIR)/lagrangian/cfdemParticle/lnInclude \
+    -I$(CFDEM_SRC_DIR)/lagrangian/cfdemParticle/cfdTools \
+    -I$(CFDEM_SRC_DIR)/recurrence/lnInclude \
+    -I$(CFDEM_SRC_DIR)/lagrangian/cfdemParticle/derived/cfdemCloudRec \
+
+EXE_LIBS = \
+    -L$(CFDEM_LIB_DIR)\
+    -lrecurrence \
+    -lturbulenceModels \
+    -lincompressibleTurbulenceModels \
+    -lincompressibleTransportModels \
+    -lfiniteVolume \
+    -lmeshTools \
+    -lfvOptions \
+    -l$(CFDEM_LIB_NAME) \
+     $(CFDEM_ADD_LIB_PATHS) \
+$(CFDEM_ADD_LIBS)

applications/solvers/rcfdemSolverForcedTracers/createFields.H

@@ -0,0 +1,113 @@
+    // dummy fields
+    Info << "\nCreating dummy density field\n" << endl;
+
+    volScalarField rho
+    (
+        IOobject
+        (
+            "rho",
+            runTime.timeName(),
+            mesh,
+            IOobject::NO_READ,
+            IOobject::NO_WRITE
+        ),
+        mesh,
+        dimensionedScalar("rho", dimensionSet(1, -3, 0, 0, 0), 1.0)
+    );
+
+    // particle fields
+    Info << "\nCreating voidfraction and particle velocity fields\n" << endl;
+
+    volScalarField voidfraction
+    (
+        IOobject
+        (
+            "voidfraction",
+            runTime.timeName(),
+            mesh,
+            IOobject::MUST_READ,
+            IOobject::AUTO_WRITE
+        ),
+        mesh
+    );
+
+    volVectorField Us
+    (
+        IOobject
+        (
+            "Us",
+            runTime.timeName(),
+            mesh,
+            IOobject::MUST_READ,
+            IOobject::AUTO_WRITE
+        ),
+        mesh
+    );
+
+    // recurrence fields
+    Info << "\nCreating recurrence fields.\n" << endl;
+
+    volScalarField pRec
+    (
+        IOobject
+        (
+            "pRec",
+            runTime.timeName(),
+            mesh,
+            IOobject::READ_IF_PRESENT,
+            IOobject::NO_WRITE
+        ),
+        mesh,
+        dimensionedScalar("p", dimensionSet(1, 2, -2, 0, 0), 1.0)
+    );
+
+    volScalarField kRec
+    (
+        IOobject
+        (
+            "kRec",
+            runTime.timeName(),
+            mesh,
+            IOobject::READ_IF_PRESENT,
+            IOobject::NO_WRITE
+        ),
+        mesh,
+        dimensionedScalar("k", dimensionSet(0, 2, -2, 0, 0), 0.0)
+    );
+
+    volVectorField URec
+    (
+        IOobject
+        (
+            "URec",
+            runTime.timeName(),
+            mesh,
+            IOobject::MUST_READ,
+            IOobject::AUTO_WRITE
+        ),
+        mesh
+    );
+
+//===============================
+
+    Info << "Calculating face flux field phi\n" << endl;
+    surfaceScalarField phiRec
+    (
+        IOobject
+        (
+            "phiRec",
+            runTime.timeName(),
+            mesh,
+            IOobject::READ_IF_PRESENT,
+            IOobject::AUTO_WRITE
+        ),
+        linearInterpolate(URec*voidfraction) & mesh.Sf()
+    );
+    phiRec.write();
+
+    singlePhaseTransportModel laminarTransport(URec, phiRec);
+
+    autoPtr<incompressible::turbulenceModel> turbulence
+    (
+        incompressible::turbulenceModel::New(URec, phiRec, laminarTransport)
+    );

applications/solvers/rcfdemSolverForcedTracers/rcfdemSolverForcedTracers.C

@@ -0,0 +1,114 @@
+/*---------------------------------------------------------------------------*\
+    CFDEMcoupling academic - Open Source CFD-DEM coupling
+
+    Contributing authors:
+    Thomas Lichtenegger
+    Copyright (C) 2015- Johannes Kepler University, Linz
+-------------------------------------------------------------------------------
+License
+    This file is part of CFDEMcoupling academic.
+
+    CFDEMcoupling academic 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.
+
+    CFDEMcoupling academic 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 CFDEMcoupling academic.  If not, see <http://www.gnu.org/licenses/>.
+
+Application
+    rcfdemSolverForcedTracers
+
+Description
+    Moves tracers according to the activated force models on pressure and velocity
+    fields provided by a recurrence process
+
+\*---------------------------------------------------------------------------*/
+
+#include "fvCFD.H"
+#include "singlePhaseTransportModel.H"
+#include "turbulentTransportModel.H"
+#include "fvOptions.H"
+
+#include "recBase.H"
+#include "recModel.H"
+
+#include "cfdemCloud.H"
+#include "clockModel.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+int main(int argc, char *argv[])
+{
+    #include "postProcess.H"
+    #include "setRootCase.H"
+    #include "createTime.H"
+    #include "createMesh.H"
+    #include "createControl.H"
+    #include "createFields.H"
+
+    cfdemCloud particleCloud(mesh);
+    recBase recurrenceBase(mesh);
+
+    const IOdictionary& recProps = mesh.lookupObject<IOdictionary>("recProperties");
+    bool useRecP(recProps.lookupOrDefault<bool>("useRecP",false));
+    bool useRecK(recProps.lookupOrDefault<bool>("useRecK",false));
+
+    // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+    Info << "\nCalculating particle trajectories based on recurrence statistics\n" << endl;
+
+    label recTimeIndex = 0;
+    label stepCounter = 0;
+    label recTimeStep2CFDTimeStep = recurrenceBase.recM().recTimeStep2CFDTimeStep();
+
+    while (runTime.run())
+    {
+        runTime++;
+
+        // do stuff (every lagrangian time step)
+        particleCloud.clockM().start(1,"Global");
+
+        Info << "Time = " << runTime.timeName() << nl << endl;
+
+        particleCloud.clockM().start(2,"Coupling");
+
+        particleCloud.evolve(voidfraction,Us,URec);
+
+        particleCloud.clockM().stop("Coupling");
+
+        stepCounter++;
+
+        if (stepCounter == recTimeStep2CFDTimeStep)
+        {
+            recurrenceBase.updateRecFields();
+            #include "updateFields.H"
+            recTimeIndex++;
+            stepCounter = 0;
+            recTimeStep2CFDTimeStep = recurrenceBase.recM().recTimeStep2CFDTimeStep();
+        }
+
+        particleCloud.clockM().start(27,"Output");
+        runTime.write();
+        particleCloud.clockM().stop("Output");
+
+        particleCloud.clockM().stop("Global");
+
+        Info << "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
+             << "  ClockTime = " << runTime.elapsedClockTime() << " s"
+             << nl << endl;
+
+    }
+
+    Info << "End\n" << endl;
+
+    return 0;
+}
+
+
+// ************************************************************************* //

applications/solvers/rcfdemSolverForcedTracers/updateFields.H

@@ -0,0 +1,13 @@
+recurrenceBase.recM().exportVolVectorField("U",URec);
+
+if (useRecP)
+{
+    recurrenceBase.recM().exportVolScalarField("p",pRec);
+}
+
+if (useRecK)
+{
+    recurrenceBase.recM().exportVolScalarField("k",kRec);
+// in case database contains the velocity variance instead of k, do
+// kRec *= 0.5;
+}

applications/solvers/rcfdemSolverRhoSteadyPimple/EEqn.H

@@ -22,19 +22,33 @@
 
     Cpv = he.name() == "e" ? thermo.Cv() : thermo.Cp();
 
+// For implict T terms in the energy/enthalpy transport equation, use
+// (he_n+1 - he_n) / (T_n+1 - T_n) = Cpv to eliminate T_n+1 with he_n+1.
+// This formula is valid for ideal gases with e=e(T) and h=h(T). For
+// incompressible fluids, e=e(T) holds, too, but enthalpy would need correction
+// terms accounting for pressure variations.
 
     fvScalarMatrix EEqn
     (
         fvm::div(phi, he)
       + addSource
       - Qsource
+      - QCoeff*T
       - fvm::Sp(QCoeff/Cpv, he)
- //     - fvm::laplacian(voidfractionRec*kf/Cpv,he)
+      + QCoeff/Cpv*he
+      - fvc::laplacian(voidfractionRec*thCond,T)
       - fvm::laplacian(voidfractionRec*thCond/Cpv,he)
+      + fvc::laplacian(voidfractionRec*thCond/Cpv,he)
      ==
         fvOptions(rho, he)
     );
 
+    if (transientEEqn)
+    {
+        EEqn += fvm::ddt(rho,voidfractionRec,he);
+    }
+
+
     EEqn.relax();
 
     fvOptions.constrain(EEqn);

applications/solvers/rcfdemSolverRhoSteadyPimple/createFields.H

@@ -168,6 +168,8 @@ Info<< "Reading thermophysical properties\n" << endl;
          linearInterpolate(rho*U*voidfraction) & mesh.Sf()
     );
 
+    Switch transientEEqn(pimple.dict().lookupOrDefault<bool>("transientEEqn",false));
+
     dimensionedScalar rhoMax
     (
         dimensionedScalar::lookupOrDefault

applications/solvers/rcfdemSolverRhoSteadyPimple/rcfdemSolverRhoSteadyPimple.C

@@ -63,7 +63,6 @@ int main(int argc, char *argv[])
     #include "createControl.H"
     #include "createTimeControls.H"
     #include "createRDeltaT.H"
-    #include "initContinuityErrs.H"
     #include "createFields.H"
     #include "createFieldRefs.H"
     #include "createFvOptions.H"
@@ -82,13 +81,13 @@ int main(int argc, char *argv[])
     // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
     label recTimeIndex = 0;
-    scalar recTimeStep = recurrenceBase.recM().recTimeStep();
-    scalar startTime = runTime.startTime().value();
+    label stepCounter = 0;
+    label recTimeStep2CFDTimeStep = recurrenceBase.recM().recTimeStep2CFDTimeStep();
 
     const IOdictionary& couplingProps = particleCloud.couplingProperties();
     label nEveryFlow(couplingProps.lookupOrDefault<label>("nEveryFlow",1));
     Info << "Solving flow equations every " << nEveryFlow << " steps.\n" << endl;
-    label stepcounter = 0;
+    label totalStepCounter = 0;
 
     Info<< "\nStarting time loop\n" << endl;
 
@@ -132,7 +131,7 @@ int main(int argc, char *argv[])
 
         particleCloud.clockM().start(26,"Flow");
         volScalarField rhoeps("rhoeps",rho*voidfractionRec);
-        if (stepcounter%nEveryFlow==0)
+        if (totalStepCounter%nEveryFlow==0)
         {
             while (pimple.loop())
             {
@@ -165,7 +164,7 @@ int main(int argc, char *argv[])
                 }
             }
         }
-        stepcounter++;
+        totalStepCounter++;
         particleCloud.clockM().stop("Flow");
 
         particleCloud.clockM().start(31,"postFlow");
@@ -173,11 +172,16 @@ int main(int argc, char *argv[])
         particleCloud.clockM().stop("postFlow");
 
         particleCloud.clockM().start(32,"ReadFields");
-        if ( runTime.timeOutputValue() - startTime - (recTimeIndex+1)*recTimeStep + 1.0e-5 > 0.0 )
+
+        stepCounter++;
+
+        if (stepCounter == recTimeStep2CFDTimeStep)
         {
             recurrenceBase.updateRecFields();
             #include "updateFields.H"
             recTimeIndex++;
+            stepCounter = 0;
+            recTimeStep2CFDTimeStep = recurrenceBase.recM().recTimeStep2CFDTimeStep();
         }
         particleCloud.clockM().stop("ReadFields");
 

applications/solvers/rcfdemSolverRhoSteadyPimpleChem/EEqn.H

@@ -0,0 +1,67 @@
+// contributions to internal energy equation can be found in
+// Crowe et al.: "Multiphase flows with droplets and particles", CRC Press 1998
+{
+    // dim he = J / kg
+    volScalarField& he = thermo.he();
+    particleCloud.energyContributions(Qsource);
+    particleCloud.energyCoefficients(QCoeff);
+
+    addSource =
+        (
+            he.name() == "e"
+          ?
+            fvc::div(phi, K) +
+            fvc::div
+            (
+                fvc::absolute(phi/fvc::interpolate(rho), voidfractionRec*U),
+                p,
+                "div(phiv,p)"
+            )
+          : fvc::div(phi, K)
+        );
+
+    Cpv = he.name() == "e" ? thermo.Cv() : thermo.Cp();
+
+// For implict T terms in the energy/enthalpy transport equation, use
+// (he_n+1 - he_n) / (T_n+1 - T_n) = Cpv to eliminate T_n+1 with he_n+1.
+// This formula is valid for ideal gases with e=e(T) and h=h(T). For
+// incompressible fluids, e=e(T) holds, too, but enthalpy would need correction
+// terms accounting for pressure variations.
+
+    fvScalarMatrix EEqn
+    (
+      fvm::div(phi, he)
+      + addSource
+      - Qsource
+      - QCoeff*T
+      - fvm::Sp(QCoeff/Cpv, he)
+      + QCoeff/Cpv*he
+      - fvc::laplacian(voidfractionRec*thCond,T)
+      - fvm::laplacian(voidfractionRec*thCond/Cpv,he)
+      + fvc::laplacian(voidfractionRec*thCond/Cpv,he)
+     ==
+        fvOptions(rho, he)
+    );
+
+    if (transientEEqn)
+    {
+        EEqn += fvm::ddt(rho,voidfractionRec,he);
+    }
+
+    EEqn.relax();
+
+    fvOptions.constrain(EEqn);
+
+    EEqn.solve();
+
+    fvOptions.correct(he);
+
+    thermo.correct();
+
+    Info<< "T max/min : " << max(T).value() << " " << min(T).value() << endl;
+
+
+    particleCloud.clockM().start(31,"energySolve");
+    particleCloud.solve();
+    particleCloud.clockM().stop("energySolve");
+}

applications/solvers/rcfdemSolverRhoSteadyPimpleChem/Make/files

@@ -0,0 +1,3 @@
+rcfdemSolverRhoSteadyPimpleChem.C
+
+EXE=$(CFDEM_APP_DIR)/rcfdemSolverRhoSteadyPimpleChem

applications/solvers/rcfdemSolverRhoSteadyPimpleChem/Make/options

@@ -0,0 +1,52 @@
+include $(CFDEM_ADD_LIBS_DIR)/additionalLibs
+
+FOAM_VERSION_MAJOR := $(word 1,$(subst ., ,$(WM_PROJECT_VERSION)))
+PFLAGS+= -DOPENFOAM_VERSION_MAJOR=$(FOAM_VERSION_MAJOR)
+PFLAGS+= -Dcompre
+
+EXE_INC = \
+     $(PFLAGS) \
+    -I$(CFDEM_OFVERSION_DIR) \
+    -I$(LIB_SRC)/transportModels/compressible/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
+    -I$(LIB_SRC)/TurbulenceModels/turbulenceModels/lnInclude \
+    -I$(LIB_SRC)/TurbulenceModels/compressible/lnInclude \
+    -I$(LIB_SRC)/finiteVolume/cfdTools \
+    -I$(LIB_SRC)/finiteVolume/lnInclude \
+    -I$(LIB_SRC)/meshTools/lnInclude \
+    -I$(LIB_SRC)/sampling/lnInclude \
+    -I$(LIB_SRC)/fvOptions/lnInclude \
+    -I$(CFDEM_SRC_DIR)/lagrangian/cfdemParticle/lnInclude \
+    -I$(CFDEM_SRC_DIR)/lagrangian/cfdemParticle/cfdTools \
+    -I$(LIB_SRC)/thermophysicalModels/specie/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/reactionThermo/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/chemistryModel/lnInclude \
+    -I$(LIB_SRC)/regionModels/regionModel/lnInclude \
+    -I$(LIB_SRC)/regionModels/surfaceFilmModels/lnInclude \
+    -I$(LIB_SRC)/ODE/lnInclude \
+    -I$(LIB_SRC)/combustionModels/lnInclude \
+    -I$(CFDEM_SRC_DIR)/recurrence/lnInclude \
+
+EXE_LIBS = \
+    -L$(CFDEM_LIB_DIR)\
+    -lrecurrence \
+    -lcompressibleTransportModels \
+    -lfluidThermophysicalModels \
+    -lspecie \
+    -lturbulenceModels \
+    -lcompressibleTurbulenceModels \
+    -lfiniteVolume \
+    -lmeshTools \
+    -lsampling \
+    -lfvOptions \
+    -l$(CFDEM_LIB_COMP_NAME) \
+     $(CFDEM_ADD_LIB_PATHS) \
+     $(CFDEM_ADD_LIBS) \
+    -lreactionThermophysicalModels \
+    -lchemistryModel \
+    -lradiationModels \
+    -lregionModels \
+    -lsurfaceFilmModels \
+    -lODE \
+    -lcombustionModels

applications/solvers/rcfdemSolverRhoSteadyPimpleChem/UEqn.H

@@ -0,0 +1,35 @@
+// Solve the Momentum equation
+particleCloud.otherForces(fOther);
+
+fvVectorMatrix UEqn
+(
+    fvm::div(phi, U)
+  + particleCloud.divVoidfractionTau(U, voidfractionRec)
+  + fvm::Sp(Ksl,U)
+  - fOther
+ ==
+    fvOptions(rho, U)
+);
+
+if (totalStepCounter%nEveryFlow==0)
+{
+    UEqn.relax();
+
+    fvOptions.constrain(UEqn);
+
+    if (modelType=="B" || modelType=="Bfull")
+    {
+        solve(UEqn == -fvc::grad(p)+ Ksl*UsRec);
+    }
+    else
+    {
+        solve(UEqn == -voidfractionRec*fvc::grad(p)+ Ksl*UsRec);
+    }
+
+
+    #include "limitU.H"
+
+    fvOptions.correct(U);
+
+    K = 0.5*magSqr(U);
+}

applications/solvers/rcfdemSolverRhoSteadyPimpleChem/YEqn.H

@@ -0,0 +1,81 @@
+particleCloud.clockM().start(29,"Y");
+
+tmp<fv::convectionScheme<scalar> > mvConvection
+(
+    fv::convectionScheme<scalar>::New
+    (
+        mesh,
+        fields,
+        phi,
+        mesh.divScheme("div(phi,Yi_h)")
+    )
+);
+
+{
+    combustion->correct();
+#if OPENFOAM_VERSION_MAJOR < 5
+    dQ = combustion->dQ();
+#else
+    Qdot = combustion->Qdot();
+#endif
+    label inertIndex = -1;
+    volScalarField Yt(0.0*Y[0]);
+
+    forAll(Y, i)
+    {
+        if (Y[i].name() == inertSpecie) inertIndex = i;
+        if (Y[i].name() != inertSpecie || propagateInertSpecie)
+        {
+            volScalarField& Yi = Y[i];
+
+            fvScalarMatrix YiEqn
+            (
+                mvConvection->fvmDiv(phi, Yi)
+              - fvm::laplacian(voidfractionRec*turbulence->muEff(), Yi)
+              ==
+                combustion->R(Yi)
+              + particleCloud.chemistryM(0).Smi(i)*p/p.prevIter()
+              + fvOptions(rho, Yi)
+            );
+
+            YiEqn.relax();
+
+            fvOptions.constrain(YiEqn);
+
+            YiEqn.solve(mesh.solver("Yi"));
+
+            Yi.relax();
+
+            fvOptions.correct(Yi);
+
+            Yi.max(0.0);
+            if (Y[i].name() != inertSpecie) Yt += Yi;
+        }
+    }
+
+    if (inertIndex!=-1)
+    {
+        Y[inertIndex].max(inertLowerBound);
+        Y[inertIndex].min(inertUpperBound);
+    }
+
+    if (propagateInertSpecie)
+    {
+        if (inertIndex!=-1) Yt /= (1-Y[inertIndex] + ROOTVSMALL);
+        forAll(Y,i)
+        {
+            if (i!=inertIndex)
+            {
+                volScalarField& Yi = Y[i];
+                Yi = Yi/(Yt+ROOTVSMALL);
+            }
+        }
+    }
+    else
+    {
+        Y[inertIndex] = scalar(1) - Yt;
+        Y[inertIndex].max(0.0);
+    }
+}
+
+particleCloud.clockM().stop("Y");

applications/solvers/rcfdemSolverRhoSteadyPimpleChem/createFieldRefs.H

@@ -0,0 +1,2 @@
+const volScalarField& T = thermo.T();
+const volScalarField& psi = thermo.psi();

applications/solvers/rcfdemSolverRhoSteadyPimpleChem/createFields.H

@@ -0,0 +1,420 @@
+Info<< "Reading thermophysical properties\n" << endl;
+
+#if OPENFOAM_VERSION_MAJOR < 6
+    Info<< "Creating combustion model\n" << endl;
+    autoPtr<combustionModels::rhoCombustionModel> combustion
+    (
+        combustionModels::rhoCombustionModel::New(mesh)
+    );
+    rhoReactionThermo& thermo = combustion->thermo();
+#else
+    Info<< "Reading thermophysical properties\n" << endl;
+    autoPtr<rhoReactionThermo> pThermo(rhoReactionThermo::New(mesh));
+    rhoReactionThermo& thermo = pThermo();
+#endif
+    thermo.validate(args.executable(), "h", "e");
+
+    basicSpecieMixture& composition = thermo.composition();
+    PtrList<volScalarField>& Y = composition.Y();
+
+    // read molecular weight
+#if OPENFOAM_VERSION_MAJOR < 6
+    volScalarField W(composition.W());
+#else
+    volScalarField W(thermo.W());
+#endif
+
+    Switch propagateInertSpecie(thermo.lookupOrDefault<bool>("propagateInertSpecie",true));
+
+    const word inertSpecie(thermo.lookupOrDefault<word>("inertSpecie","none"));
+
+    const scalar inertLowerBound(thermo.lookupOrDefault<scalar>("inertLowerBound",0.0));
+
+    const scalar inertUpperBound(thermo.lookupOrDefault<scalar>("inertUpperBound",1.0));
+
+    if (!composition.contains(inertSpecie) && inertSpecie != "none")
+    {
+        FatalErrorIn(args.executable())
+            << "Specified inert specie '" << inertSpecie << "' not found in "
+            << "species list. Available species:" << composition.species()
+            << exit(FatalError);
+    }
+
+    Info<< "inert will be bounded in [" << inertLowerBound << "," << inertUpperBound << "]" << endl;
+
+    volScalarField& p = thermo.p();
+
+    multivariateSurfaceInterpolationScheme<scalar>::fieldTable fields;
+
+    forAll(Y, i)
+    {
+        fields.add(Y[i]);
+    }
+    fields.add(thermo.he());
+
+    Info<< "Reading field rho\n" << endl;
+    volScalarField rho
+    (
+        IOobject
+        (
+            "rho",
+            runTime.timeName(),
+            mesh,
+            IOobject::READ_IF_PRESENT,
+            IOobject::AUTO_WRITE
+        ),
+        thermo.rho()
+    );
+
+
+    Info<< "Reading field U\n" << endl;
+    volVectorField U
+    (
+        IOobject
+        (
+            "U",
+            runTime.timeName(),
+            mesh,
+            IOobject::MUST_READ,
+            IOobject::AUTO_WRITE
+        ),
+        mesh
+    );
+
+    Info<< "\nReading voidfraction field voidfraction = (Vgas/Vparticle)\n" << endl;
+    volScalarField voidfraction
+    (
+        IOobject
+        (
+            "voidfraction",
+            runTime.timeName(),
+            mesh,
+            IOobject::MUST_READ,
+            IOobject::AUTO_WRITE
+        ),
+        mesh
+    );
+
+    volScalarField voidfractionRec
+    (
+        IOobject
+        (
+            "voidfractionRec",
+            runTime.timeName(),
+            mesh,
+            IOobject::NO_READ,
+            IOobject::AUTO_WRITE
+        ),
+        voidfraction
+    );
+
+    volScalarField addSource
+    (
+        IOobject
+        (
+           "addSource",
+           runTime.timeName(),
+           mesh,
+           IOobject::READ_IF_PRESENT,
+           IOobject::AUTO_WRITE
+        ),
+        mesh,
+        dimensionedScalar("zero", dimensionSet(1,-1,-3,0,0,0,0), 0.0)
+    );
+
+    Info<< "\nCreating fluid-particle heat flux field\n" << endl;
+    volScalarField Qsource
+    (
+        IOobject
+        (
+            "Qsource",
+            runTime.timeName(),
+            mesh,
+            IOobject::NO_READ,
+            IOobject::AUTO_WRITE
+        ),
+        mesh,
+        dimensionedScalar("zero", dimensionSet(1,-1,-3,0,0,0,0), 0.0)
+    );
+
+    Info<< "\nCreating fluid-particle heat flux coefficient field\n" << endl;
+    volScalarField QCoeff
+    (
+        IOobject
+        (
+            "QCoeff",
+            runTime.timeName(),
+            mesh,
+            IOobject::NO_READ,
+            IOobject::AUTO_WRITE
+        ),
+        mesh,
+        dimensionedScalar("zero", dimensionSet(1,-1,-3,-1,0,0,0), 0.0)
+    );
+
+    Info<< "\nCreating fluid thermal conduction field\n" << endl;
+    volScalarField QFluidCond
+    (
+        IOobject
+        (
+            "QFluidCond",
+            runTime.timeName(),
+            mesh,
+            IOobject::NO_READ,
+            IOobject::NO_WRITE
+        ),
+        mesh,
+        dimensionedScalar("zero", dimensionSet(1,-1,-3,0,0,0,0), 0.0)
+    );
+
+
+    Info<< "\nCreating thermal conductivity field\n" << endl;
+    volScalarField thCond
+    (
+        IOobject
+        (
+            "thCond",
+            runTime.timeName(),
+            mesh,
+            IOobject::READ_IF_PRESENT,
+            IOobject::AUTO_WRITE
+        ),
+        mesh,
+        dimensionedScalar("zero", dimensionSet(1,1,-3,-1,0,0,0), 0.0),
+        "zeroGradient"
+    );
+
+    Info<< "\nCreating heat capacity field\n" << endl;
+    volScalarField Cpv
+    (
+        IOobject
+        (
+            "Cpv",
+            runTime.timeName(),
+            mesh,
+            IOobject::NO_READ,
+            IOobject::AUTO_WRITE
+        ),
+        mesh,
+        dimensionedScalar("zero", dimensionSet(0,2,-2,-1,0,0,0), 0.0)
+    );
+
+    Info<< "\nCreating body force field\n" << endl;
+    volVectorField fOther
+    (
+        IOobject
+        (
+            "fOther",
+            runTime.timeName(),
+            mesh,
+            IOobject::NO_READ,
+            IOobject::NO_WRITE
+        ),
+        mesh,
+        dimensionedVector("zero", dimensionSet(1,-2,-2,0,0,0,0), vector::zero)
+    );
+
+    Info<< "Reading/calculating face flux field phi\n" << endl;
+    surfaceScalarField phi
+    (
+        IOobject
+        (
+            "phi",
+            runTime.timeName(),
+            mesh,
+            IOobject::READ_IF_PRESENT,
+            IOobject::AUTO_WRITE
+         ),
+         linearInterpolate(rho*U*voidfraction) & mesh.Sf()
+    );
+
+    Switch transientEEqn(pimple.dict().lookupOrDefault<bool>("transientEEqn",false));
+
+    dimensionedScalar rhoMax
+    (
+        dimensionedScalar::lookupOrDefault
+        (
+            "rhoMax",
+            pimple.dict(),
+            dimDensity,
+            GREAT
+        )
+    );
+
+    dimensionedScalar rhoMin
+    (
+        dimensionedScalar::lookupOrDefault
+        (
+            "rhoMin",
+            pimple.dict(),
+            dimDensity,
+            0
+        )
+    );
+
+    dimensionedScalar pMax
+    (
+        dimensionedScalar::lookupOrDefault
+        (
+            "pMax",
+            pimple.dict(),
+            dimPressure,
+            GREAT
+        )
+    );
+
+    dimensionedScalar pMin
+    (
+        dimensionedScalar::lookupOrDefault
+        (
+            "pMin",
+            pimple.dict(),
+            dimPressure,
+            -GREAT
+        )
+    );
+
+    dimensionedScalar UMax
+    (
+        dimensionedScalar::lookupOrDefault
+        (
+            "UMax",
+            pimple.dict(),
+            dimVelocity,
+            -1.0
+        )
+    );
+
+    Info<< "Creating turbulence model\n" << endl;
+    autoPtr<compressible::turbulenceModel> turbulence
+    (
+        compressible::turbulenceModel::New
+        (
+            rho,
+            U,
+            phi,
+            thermo
+        )
+    );
+
+#if OPENFOAM_VERSION_MAJOR >= 6
+    Info<< "Creating combustion model\n" << endl;
+    autoPtr<CombustionModel<rhoReactionThermo>> combustion
+    (
+        CombustionModel<rhoReactionThermo>::New(thermo, turbulence())
+    );
+#endif
+
+    label pRefCell = 0;
+    scalar pRefValue = 0.0;
+    setRefCell(p, pimple.dict(), pRefCell, pRefValue);
+
+    mesh.setFluxRequired(p.name());
+
+    Info<< "Creating field dpdt\n" << endl;
+    volScalarField dpdt
+    (
+        IOobject
+        (
+            "dpdt",
+            runTime.timeName(),
+            mesh
+        ),
+        mesh,
+        dimensionedScalar("dpdt", p.dimensions()/dimTime, 0)
+    );
+
+    Info<< "Creating field kinetic energy K\n" << endl;
+    volScalarField K("K", 0.5*magSqr(U));
+
+#if OPENFOAM_VERSION_MAJOR < 5
+    volScalarField dQ
+    (
+        IOobject
+        (
+            "dQ",
+            runTime.timeName(),
+            mesh,
+            IOobject::NO_READ,
+            IOobject::AUTO_WRITE
+        ),
+        mesh,
+        dimensionedScalar("dQ", dimEnergy/dimTime, 0.0)
+    );
+#else
+    volScalarField Qdot
+    (
+        IOobject
+        (
+            "Qdot",
+            runTime.timeName(),
+            mesh,
+            IOobject::NO_READ,
+            IOobject::AUTO_WRITE
+        ),
+        mesh,
+        dimensionedScalar("Qdot", dimEnergy/dimVolume/dimTime, 0.0)
+    );
+#endif
+
+    Info<< "\nReading momentum exchange field Ksl\n" << endl;
+    volScalarField Ksl
+    (
+        IOobject
+        (
+            "Ksl",
+            runTime.timeName(),
+            mesh,
+            IOobject::READ_IF_PRESENT,
+            IOobject::AUTO_WRITE
+        ),
+        mesh,
+        dimensionedScalar("0", dimensionSet(1, -3, -1, 0, 0), 0.0)
+    );
+
+
+    Info<< "Reading particle velocity field Us\n" << endl;
+    volVectorField Us
+    (
+        IOobject
+        (
+            "Us",
+            runTime.timeName(),
+            mesh,
+            IOobject::MUST_READ,
+            IOobject::AUTO_WRITE
+        ),
+        mesh
+    );
+
+    volScalarField molarConc
+    (
+        IOobject
+        (
+            "molarConc",
+            runTime.timeName(),
+            mesh,
+            IOobject::NO_READ,
+            IOobject::AUTO_WRITE
+        ),
+        mesh,
+        dimensionedScalar("zero",dimensionSet(0, -3, 0, 0, 1),0)
+    );
+
+    volVectorField UsRec
+    (
+        IOobject
+        (
+            "UsRec",
+            runTime.timeName(),
+            mesh,
+            IOobject::NO_READ,
+            IOobject::AUTO_WRITE
+        ),
+        Us
+    );
+
+
+    dimensionedScalar kf("0", dimensionSet(1, 1, -3, -1, 0, 0, 0), 0.026);
+
+//===============================

applications/solvers/rcfdemSolverRhoSteadyPimpleChem/limitP.H

@@ -0,0 +1,2 @@
+p = max(p, pMin);
+p = min(p, pMax);

applications/solvers/rcfdemSolverRhoSteadyPimpleChem/limitU.H

@@ -0,0 +1,11 @@
+if (UMax.value() > 0)
+{
+    forAll(U,cellI)
+    {
+        scalar mU(mag(U[cellI]));
+        if (mU > UMax.value())
+        {
+            U[cellI] *= UMax.value() / mU;
+        }
+    }
+}

applications/solvers/rcfdemSolverRhoSteadyPimpleChem/molConc.H

@@ -0,0 +1,12 @@
+{
+    molarConc = 0.0 * molarConc;
+    forAll(Y, i)
+    {
+        volScalarField& Yi = Y[i];
+        dimensionedScalar mi("mi",dimensionSet(1, 0, 0, 0, -1),composition.W(i));
+        mi /= 1000.0; // g to kg
+        molarConc += rho * Yi / mi;
+    }
+}
+
+// ************************************************************************* //

applications/solvers/rcfdemSolverRhoSteadyPimpleChem/monitorMass.H

@@ -0,0 +1,9 @@
+{
+    m=gSum(rhoeps*1.0*rhoeps.mesh().V());
+    if(counter==0) m0=m;
+    counter++;
+    Info << "\ncurrent gas mass = " << m << "\n" << endl;
+    Info << "\ncurrent added gas mass = " << m-m0 << "\n" << endl;
+
+    QFluidCond = fvc::laplacian(voidfractionRec*thCond,T);
+}

applications/solvers/rcfdemSolverRhoSteadyPimpleChem/pEqn.H

@@ -0,0 +1,96 @@
+rho = thermo.rho();
+rho = max(rho, rhoMin);
+rho = min(rho, rhoMax);
+rho.relax();
+
+if (totalStepCounter%nEveryFlow==0)
+{
+
+volScalarField rAU(1.0/UEqn.A());
+surfaceScalarField rhorAUf("rhorAUf", fvc::interpolate(rhoeps*rAU));
+if (modelType=="A")
+{
+    rhorAUf *= fvc::interpolate(voidfractionRec);
+}
+volVectorField HbyA(constrainHbyA(rAU*UEqn.H(), U, p));
+
+surfaceScalarField phiUs("phiUs", fvc::interpolate(rhoeps*rAU*Ksl*UsRec)& mesh.Sf());
+
+
+if (pimple.transonic())
+{
+//  transonic version not implemented yet
+}
+else
+{
+    surfaceScalarField phiHbyA
+    (
+        "phiHbyA",
+        (
+            fvc::flux(rhoeps*HbyA)
+        )
+    );
+
+    // flux without pressure gradient contribution
+    phi = phiHbyA + phiUs;
+
+    // Update the pressure BCs to ensure flux consistency
+    constrainPressure(p, rhoeps, U, phi, rhorAUf);
+
+    volScalarField SmbyP(particleCloud.chemistryM(0).Sm() / p);
+
+    while (pimple.correctNonOrthogonal())
+    {
+        // Pressure corrector
+        fvScalarMatrix pEqn
+        (
+            fvc::div(phi)
+          - fvm::laplacian(rhorAUf, p)
+          ==
+            fvm::Sp(SmbyP, p)
+          + fvOptions(psi, p, rho.name())
+        );
+
+        pEqn.setReference(pRefCell, pRefValue);
+
+        pEqn.solve();
+
+        if (pimple.finalNonOrthogonalIter())
+        {
+            phi += pEqn.flux();
+        }
+    }
+}
+
+#include "rhoEqn.H"
+#include "compressibleContinuityErrsPU.H"
+
+// Explicitly relax pressure for momentum corrector
+p.relax();
+
+#include "limitP.H"
+
+// Recalculate density from the relaxed pressure
+rho = thermo.rho();
+rho = max(rho, rhoMin);
+rho = min(rho, rhoMax);
+rho.relax();
+Info<< "rho max/min : " << max(rho).value()
+    << " " << min(rho).value() << endl;
+
+if (modelType=="A")
+{
+    U = HbyA - rAU*(voidfractionRec*fvc::grad(p)-Ksl*UsRec);
+}
+else
+{
+    U = HbyA - rAU*(fvc::grad(p)-Ksl*UsRec);
+}
+
+#include "limitU.H"
+
+U.correctBoundaryConditions();
+fvOptions.correct(U);
+K = 0.5*magSqr(U);
+
+}

applications/solvers/rcfdemSolverRhoSteadyPimpleChem/rcfdemSolverRhoSteadyPimpleChem.C

@@ -0,0 +1,222 @@
+/*---------------------------------------------------------------------------*\
+License
+
+    This 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.
+
+    This code 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 this code.  If not, see <http://www.gnu.org/licenses/>.
+
+    Copyright (C) 2015- Thomas Lichtenegger, JKU Linz, Austria
+
+Application
+    rcfdemSolverRhoSteadyPimpleChem
+
+Description
+    Transient (DEM) + steady-state (CFD) solver for compressible flow using the
+    flexible PIMPLE (PISO-SIMPLE) algorithm. Particle-motion is obtained from
+    a recurrence process.
+    Turbulence modelling is generic, i.e. laminar, RAS or LES may be selected.
+    The code is an evolution of the solver rhoPimpleFoam in OpenFOAM(R) 4.x,
+    where additional functionality for CFD-DEM coupling is added.
+\*---------------------------------------------------------------------------*/
+
+#include "fvCFD.H"
+//#include "psiThermo.H"
+#include "turbulentFluidThermoModel.H"
+#if OPENFOAM_VERSION_MAJOR < 6
+#include "rhoCombustionModel.H"
+#else
+#include "rhoReactionThermo.H"
+#include "CombustionModel.H"
+#endif
+#include "bound.H"
+#include "pimpleControl.H"
+#include "fvOptions.H"
+#include "localEulerDdtScheme.H"
+#include "fvcSmooth.H"
+
+#include "cfdemCloudRec.H"
+#include "recBase.H"
+#include "recModel.H"
+#include "recPath.H"
+
+#include "cfdemCloudEnergy.H"
+#include "implicitCouple.H"
+#include "clockModel.H"
+#include "smoothingModel.H"
+#include "forceModel.H"
+#include "thermCondModel.H"
+#include "energyModel.H"
+#include "chemistryModel.H"
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+int main(int argc, char *argv[])
+{
+    #include "postProcess.H"
+
+    #include "setRootCase.H"
+    #include "createTime.H"
+    #include "createMesh.H"
+    #include "createControl.H"
+    #include "createTimeControls.H"
+    #include "createRDeltaT.H"
+
+    #include "initContinuityErrs.H"
+    #include "createFields.H"
+    #include "createFieldRefs.H"
+    #include "createFvOptions.H"
+
+    // create cfdemCloud
+    //#include "readGravitationalAcceleration.H"
+    cfdemCloudRec<cfdemCloudEnergy> particleCloud(mesh);
+    #include "checkModelType.H"
+    recBase recurrenceBase(mesh);
+    #include "updateFields.H"
+
+    turbulence->validate();
+    //#include "compressibleCourantNo.H"
+    //#include "setInitialDeltaT.H"
+
+    // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+    label recTimeIndex = 0;
+    label stepCounter = 0;
+    label recTimeStep2CFDTimeStep = recurrenceBase.recM().recTimeStep2CFDTimeStep();
+
+    const IOdictionary& couplingProps = particleCloud.couplingProperties();
+    label nEveryFlow(couplingProps.lookupOrDefault<label>("nEveryFlow",1));
+    Info << "Solving flow equations for U and p every " << nEveryFlow << " steps.\n" << endl;
+    label totalStepCounter = 0;
+
+    Info<< "\nStarting time loop\n" << endl;
+
+    scalar m(0.0);
+    scalar m0(0.0);
+    label counter(0);
+    p.storePrevIter();
+
+    while (runTime.run())
+    {
+        #include "readTimeControls.H"
+        #include "compressibleCourantNo.H"
+        #include "setDeltaT.H"
+
+        runTime++;
+
+        particleCloud.clockM().start(1,"Global");
+
+        Info<< "Time = " << runTime.timeName() << nl << endl;
+
+        // do particle stuff
+        particleCloud.clockM().start(2,"Coupling");
+        bool hasEvolved = particleCloud.evolve(voidfraction,Us,U);
+
+        //voidfraction = voidfractionRec;
+        //Us = UsRec;
+
+        if(hasEvolved)
+        {
+            particleCloud.smoothingM().smoothen(particleCloud.forceM(0).impParticleForces());
+        }
+
+        Info << "update Ksl.internalField()" << endl;
+        Ksl = particleCloud.momCoupleM(0).impMomSource();
+        Ksl.correctBoundaryConditions();
+
+        //Force Checks
+        vector fTotal(0,0,0);
+        vector fImpTotal = sum(mesh.V()*Ksl.primitiveFieldRef()*(Us.primitiveFieldRef()-U.primitiveFieldRef()));
+        reduce(fImpTotal, sumOp<vector>());
+        Info << "TotalForceExp: " << fTotal << endl;
+        Info << "TotalForceImp: " << fImpTotal << endl;
+
+        #include "solverDebugInfo.H"
+        particleCloud.clockM().stop("Coupling");
+
+        particleCloud.clockM().start(26,"Flow");
+        volScalarField rhoeps("rhoeps",rho*voidfractionRec);
+
+        while (pimple.loop())
+        {
+            // if needed, perform drag update here
+#if OPENFOAM_VERSION_MAJOR < 6
+            if (pimple.nCorrPIMPLE() <= 1)
+#else
+            if (pimple.nCorrPimple() <= 1)
+#endif
+            {
+                #include "rhoEqn.H"
+            }
+
+            // --- Pressure-velocity PIMPLE corrector loop
+
+
+            #include "UEqn.H"
+            #include "EEqn.H"
+
+            // --- Pressure corrector loop
+            while (pimple.correct())
+            {
+                // besides this pEqn, OF offers a "pimple consistent"-option
+                #include "molConc.H"
+                #include "pEqn.H"
+                rhoeps=rho*voidfractionRec;
+            }
+            #include "YEqn.H"
+
+            if (pimple.turbCorr())
+            {
+                turbulence->correct();
+            }
+        }
+
+        #include "monitorMass.H"
+
+        totalStepCounter++;
+        particleCloud.clockM().stop("Flow");
+
+        particleCloud.clockM().start(31,"postFlow");
+        particleCloud.postFlow();
+        particleCloud.clockM().stop("postFlow");
+
+        particleCloud.clockM().start(32,"ReadFields");
+        stepCounter++;
+
+        if (stepCounter == recTimeStep2CFDTimeStep)
+        {
+            recurrenceBase.updateRecFields();
+            #include "updateFields.H"
+            recTimeIndex++;
+            stepCounter = 0;
+            recTimeStep2CFDTimeStep = recurrenceBase.recM().recTimeStep2CFDTimeStep();
+        }
+        particleCloud.clockM().stop("ReadFields");
+
+        runTime.write();
+
+
+        Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
+            << "  ClockTime = " << runTime.elapsedClockTime() << " s"
+            << nl << endl;
+
+
+        particleCloud.clockM().stop("Global");
+    }
+
+    Info<< "End\n" << endl;
+
+    return 0;
+}
+
+
+// ************************************************************************* //

applications/solvers/rcfdemSolverRhoSteadyPimpleChem/rhoEqn.H

@@ -0,0 +1,21 @@
+{
+/*
+    fvScalarMatrix rhoEqn
+    (
+        //fvm::ddt(voidfraction,rho)
+      //+ 
+        fvc::div(phi)
+      ==
+        particleCloud.chemistryM(0).Sm()
+      + fvOptions(rho)
+    );
+
+    fvOptions.constrain(rhoEqn);
+
+    rhoEqn.solve();
+
+    fvOptions.correct(rho);
+*/
+}
+
+// ************************************************************************* //

applications/solvers/rcfdemSolverRhoSteadyPimpleChem/updateFields.H

@@ -1,4 +1,8 @@
+// is it neccessary to extend recurrence path?
+if(recurrenceBase.recM().endOfPath())
+{
+    recurrenceBase.extendPath();
+}
+
 recurrenceBase.recM().exportVolScalarField("voidfraction",voidfractionRec);
-recurrenceBase.recM().exportVolVectorField("U",URec);
 recurrenceBase.recM().exportVolVectorField("Us",UsRec);
-recurrenceBase.recM().exportSurfaceScalarField("phi",phiRec);

applications/solvers/rctfSpeciesTransport/CEq.H

@@ -0,0 +1,19 @@
+
+volScalarField alphaEff("alphaEff", turbulence->nu()/Sc + alphat);
+
+CEqn =
+(
+    fvm::ddt(C)
+    + fvm::div(phiRec, C)
+    - fvm::laplacian(alphaEff, C)
+    ==
+    fvOptions(C)
+);
+
+CEqn.relax(relaxCoeff);
+
+fvOptions.constrain(CEqn);
+
+CEqn.solve();
+
+fvOptions.correct(C);

applications/solvers/rctfSpeciesTransport/Make/files

@@ -0,0 +1,3 @@
+rctfSpeciesTransport.C
+
+EXE=$(CFDEM_APP_DIR)/rctfSpeciesTransport

applications/solvers/rctfSpeciesTransport/Make/options

@@ -0,0 +1,27 @@
+include $(CFDEM_ADD_LIBS_DIR)/additionalLibs
+
+EXE_INC = \
+    -I$(CFDEM_OFVERSION_DIR) \
+    -I$(LIB_SRC)/finiteVolume/lnInclude \
+    -I$(LIB_SRC)/meshTools/lnInclude \
+    -I$(LIB_SRC)/TurbulenceModels/turbulenceModels/lnInclude \
+    -I$(LIB_SRC)/TurbulenceModels/incompressible/lnInclude \
+    -I$(LIB_SRC)/transportModels \
+    -I$(LIB_SRC)/transportModels/incompressible/singlePhaseTransportModel \
+    -I$(CFDEM_SRC_DIR)/lagrangian/cfdemParticle/lnInclude \
+    -I$(CFDEM_SRC_DIR)/lagrangian/cfdemParticle/cfdTools \
+    -I$(CFDEM_SRC_DIR)/recurrence/lnInclude \
+    -I$(CFDEM_SRC_DIR)/lagrangian/cfdemParticle/derived/cfdemCloudRec \
+
+EXE_LIBS = \
+    -L$(CFDEM_LIB_DIR)\
+    -lrecurrence \
+    -lturbulenceModels \
+    -lincompressibleTurbulenceModels \
+    -lincompressibleTransportModels \
+    -lfiniteVolume \
+    -lmeshTools \
+    -lfvOptions \
+    -l$(CFDEM_LIB_NAME) \
+     $(CFDEM_ADD_LIB_PATHS) \
+     $(CFDEM_ADD_LIBS)

applications/solvers/rctfSpeciesTransport/continuityErrCalc.H

@@ -0,0 +1,17 @@
+// calculate the continuity error according to phiRec
+
+{
+    volScalarField contErr(fvc::div(phiRec));
+
+    scalar sumLocalContErr = runTime.deltaTValue()*
+        mag(contErr)().weightedAverage(mesh.V()).value();
+
+    scalar globalContErr = runTime.deltaTValue()*
+        contErr.weightedAverage(mesh.V()).value();
+    cumulativeContErr += globalContErr;
+
+    Info<< "time step continuity errors : sum local = " << sumLocalContErr
+        << ", global = " << globalContErr
+        << ", cumulative = " << cumulativeContErr
+        << endl;
+}

applications/solvers/rctfSpeciesTransport/createFields.H

@@ -0,0 +1,280 @@
+//creating the fields according to the recurrence dictionary
+
+IOdictionary recProperties_
+(
+    IOobject
+    (
+        "recProperties0",
+        runTime.constant(),
+        mesh,
+        IOobject::MUST_READ,
+        IOobject::NO_WRITE
+    )
+);
+
+List<wordList> fieldsDict_(recProperties_.lookup("fieldsPairs"));
+
+wordList fieldNames(fieldsDict_.size());
+
+for(int i = 0; i < fieldsDict_.size(); i++)
+{
+
+    fieldNames[i]= fieldsDict_[i][0];
+}
+
+Info<< "\n list of the fields: \n" << fieldNames << endl;
+
+//reading coherent velocity field name
+label k = findIndex(fieldNames,"coh_velocity");
+
+if (k < 0)
+{
+    FatalError <<"\n No field is defiened for the coherent velocity\n" << abort(FatalError);
+}
+const word Ucoh_pair = fieldsDict_[k][1];
+
+    volVectorField UcohRec
+    (
+        IOobject
+        (
+            Ucoh_pair,
+            runTime.timeName(),
+            mesh,
+            IOobject::READ_IF_PRESENT,
+            IOobject::AUTO_WRITE
+        ),
+        mesh,
+        dimensionedVector("zero",dimensionSet(0, 1, -1, 0, 0),vector::zero)
+    );
+
+//reading incoherent velocity field name
+k = findIndex(fieldNames,"inc_velocity");
+
+if (k < 0)
+{
+    FatalError <<"\n No field is defiened for the incoherent velocity\n" << abort(FatalError);
+}
+const word Uinc_pair = fieldsDict_[k][1];
+
+    volVectorField UincRec
+    (
+        IOobject
+        (
+            Uinc_pair,
+            runTime.timeName(),
+            mesh,
+            IOobject::READ_IF_PRESENT,
+            IOobject::AUTO_WRITE
+        ),
+        mesh,
+        dimensionedVector("zero",dimensionSet(0, 1, -1, 0, 0),vector::zero)
+    );
+
+//reading coherent turb kinetic energy field name
+k = findIndex(fieldNames,"kSGS_coh");
+if (k < 0)
+{
+    FatalError <<"\n No field is defiened for the coherent subgrid-scale turbulent kinetic energy\n" << abort(FatalError);
+}
+const word kSGScoh_pair = fieldsDict_[k][1];
+
+    volScalarField kcohRec
+    (
+        IOobject
+        (
+            kSGScoh_pair,
+            runTime.timeName(),
+            mesh,
+            IOobject::READ_IF_PRESENT,
+            IOobject::AUTO_WRITE
+        ),
+        mesh,
+        dimensionedScalar("zero",dimensionSet(0, 2, -2, 0, 0),0.0)
+     );
+
+//reading incoherent turb kinetic energy field name
+k = findIndex(fieldNames,"kSGS_inc");
+if (k < 0)
+{
+    FatalError <<"\n No field is defiened for the coherent subgrid-scale turbulent kinetic energy\n" << abort(FatalError);
+}
+const word kSGSinc_pair = fieldsDict_[k][1];
+
+    volScalarField kincRec
+    (
+        IOobject
+        (
+            kSGSinc_pair,
+            runTime.timeName(),
+            mesh,
+            IOobject::READ_IF_PRESENT,
+            IOobject::AUTO_WRITE
+        ),
+        mesh,
+        dimensionedScalar("zero",dimensionSet(0, 2, -2, 0, 0),0.0)
+     );
+
+// calculated fields
+Info<< "\nCreating cell volume field\n" << endl;
+
+volScalarField delta
+(
+    IOobject
+    (
+        "delta",
+        runTime.timeName(),
+        mesh,
+        IOobject::NO_READ,
+        IOobject::NO_WRITE
+    ),
+   mesh,
+   dimensionedScalar("delta", dimLength, 0.0)
+);
+
+delta.primitiveFieldRef()=pow(mesh.V(),1.0/3.0);
+delta.write();
+
+
+volVectorField URec
+(
+    IOobject
+    (
+        "URec",
+        runTime.timeName(),
+        mesh,
+        IOobject::MUST_READ,
+        IOobject::AUTO_WRITE
+    ),
+    mesh
+);
+Info<< "\nCreating turb kinetic energy field\n" << endl;
+
+volScalarField kRec
+(
+    IOobject
+    (
+        "kRec",
+        runTime.timeName(),
+        mesh,
+        IOobject::NO_READ,
+        IOobject::AUTO_WRITE
+    ),
+    kcohRec+kincRec
+);
+
+// check if there is any negative values
+forAll(kRec, cellI)
+{
+    if (kRec[cellI] < SMALL)
+    {
+        kRec[cellI] = 0.0;
+    }
+
+}
+
+const fvPatchList& patches = mesh.boundary();
+forAll(patches, patchI)
+{
+        kRec.boundaryFieldRef()[patchI] = 0.0;
+}
+
+
+kRec.write();
+
+Info<< "\nCreating turb viscosity field\n" << endl;
+volScalarField nutRec
+(
+    IOobject
+    (
+        "nutRec",
+        runTime.timeName(),
+        mesh,
+        IOobject::NO_READ,
+        IOobject::AUTO_WRITE
+    ),
+    sqrt(kRec)*delta*0.094
+);
+
+nutRec.write();
+
+
+Info<< "Calculating face flux field phiRec\n" << endl;
+surfaceScalarField phiRec
+ (
+     IOobject
+     (
+        "phiRec",
+        runTime.timeName(),
+        mesh,
+        IOobject::NO_READ,
+        IOobject::AUTO_WRITE
+     ),
+     linearInterpolate(URec) & mesh.Sf()
+ );
+
+    phiRec.write();
+
+    singlePhaseTransportModel laminarTransport(URec, phiRec);
+
+    autoPtr<incompressible::turbulenceModel> turbulence
+    (
+        incompressible::turbulenceModel::New(URec, phiRec, laminarTransport)
+    );
+  
+    dimensionedScalar Sc("Sc", dimless, laminarTransport);
+    dimensionedScalar Sct("Sct", dimless, laminarTransport);
+  
+volScalarField alphat
+(
+    IOobject
+    (
+        "alphat",
+        runTime.timeName(),
+        mesh,
+        IOobject::NO_READ,
+        IOobject::AUTO_WRITE
+    ),
+    nutRec/Sct
+);
+
+ // create the scalar field
+ Info<< "Creating scalar transport field\n" << endl;
+
+    volScalarField C
+    (
+        IOobject
+        (
+            "C",
+            runTime.timeName(),
+            mesh,
+            IOobject::MUST_READ,
+            IOobject::AUTO_WRITE
+        ),
+        mesh
+    );
+    
+    fvScalarMatrix CEqn(C, dimless*dimVolume/(dimTime));
+
+
+ Info<< "reading clockProperties\n" << endl;
+
+    IOdictionary clockProperties
+    (
+            IOobject
+            (
+                "clockProperties",
+                mesh.time().constant(),
+                mesh,
+                IOobject::MUST_READ,
+                IOobject::NO_WRITE
+            )
+    );
+
+    autoPtr<clockModel> myClock
+    (
+            clockModel::New
+            (
+                clockProperties,
+                mesh.time()
+            )
+    );

applications/solvers/rctfSpeciesTransport/createRecBase.H

@@ -0,0 +1,66 @@
+// check which recProperties dicts are present, read them in and construct a PtrList of recBases
+// names for dicts can be "recProperties" or "recPropertiesN" where N in {0, 1, ...} 
+
+#include "error.H"
+
+word dictName = "recProperties";
+wordList recPropertiesList(0);
+PtrList <recBase> recBases(0);
+label maxDictNumber = 100;
+
+{
+    IOdictionary recPropDict
+    (
+        IOobject
+        (
+            dictName,
+            mesh.time().constant(),
+            mesh,
+            IOobject::READ_IF_PRESENT,
+            IOobject::NO_WRITE
+        )
+    );
+
+    if (recPropDict.headerOk())
+    {
+        recPropertiesList.append(dictName);
+    }
+}
+
+
+for (label counter = 0; counter < maxDictNumber; counter++)
+{
+    word dictNameIter = dictName + Foam::name(counter);
+    IOdictionary recPropDict
+    (
+        IOobject
+        (
+            dictNameIter,
+            mesh.time().constant(),
+            mesh,
+            IOobject::READ_IF_PRESENT,
+            IOobject::NO_WRITE
+        )
+    );
+
+    if (recPropDict.headerOk())
+    {
+        recPropertiesList.append(dictNameIter);
+    }
+}
+
+if (recPropertiesList.size() == 0)
+{
+    FatalError << "no recProperties dicts found" << endl;
+}
+
+else
+{
+    Info << "found " << recPropertiesList.size() << " dicts with names " << recPropertiesList << endl;
+}
+
+
+for (label counter = 0; counter < recPropertiesList.size(); counter++)
+{
+    recBases.append( new recBase(mesh, recPropertiesList[counter]));
+}

applications/solvers/rctfSpeciesTransport/rctfSpeciesTransport.C

@@ -0,0 +1,132 @@
+/*---------------------------------------------------------------------------*\
+    CFDEMcoupling academic - Open Source CFD-DEM coupling
+
+    Contributing authors:
+    Thomas Lichtenegger, Gerhard Holzinger, Sanaz Abbasi
+    Copyright (C) 2015- Johannes Kepler University, Linz
+-------------------------------------------------------------------------------
+License
+    This file is part of CFDEMcoupling academic.
+
+    CFDEMcoupling academic 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.
+
+    CFDEMcoupling academic 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 CFDEMcoupling academic.  If not, see <http://www.gnu.org/licenses/>.
+
+Application
+    Turbulent Transport Recurrence Solver for modal decomposition
+
+Description
+    Solves a transport equation for a passive scalar on a single-phase solution
+    for a solver based on recurrence statistics
+
+Rules
+    Solution data to compute the recurrence statistics from, needs to
+    reside in $CASE_ROOT/dataBase(0...N)
+    Time step data in the first dataBase needs to be evenly spaced in time
+    A list of indices for the corresponding incoherent fields to coherent ones
+    should be provided.
+
+\*---------------------------------------------------------------------------*/
+
+#include "fvCFD.H"
+#include "singlePhaseTransportModel.H"
+#include "turbulentTransportModel.H"
+#include "fvOptions.H"
+
+#include "recBase.H"
+#include "recModel.H"
+
+#include "clockModel.H"
+
+#include "objectRegistry.H"
+#include "VectorSpace.H"
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+int main(int argc, char *argv[])
+{
+    #include "postProcess.H"
+    #include "setRootCase.H"
+    #include "createTime.H"
+    #include "createMesh.H"
+    #include "createControl.H"
+    #include "initContinuityErrs.H"
+    #include "createFields.H"
+    #include "createFvOptions.H"
+    scalar relaxCoeff(0.0);
+
+    //create recBases according to a list of recProperties
+    #include "createRecBase.H"
+    // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+
+    label recTimeIndex(0);
+    label currTimeIndex(0);
+
+    scalar recTimeStep_=recBases[0].recM().recTimeStep();
+    labelPairList incPairTimeIndex_(0);
+
+    IFstream pairFile("incIndexPairList");
+    pairFile >> incPairTimeIndex_;
+
+    while (runTime.run())
+    {
+
+        myClock().start(1,"Global");
+        runTime++;
+
+        myClock().start(11,"Total");
+
+        Info<< "Time = " << runTime.timeName() << nl << endl;
+
+        myClock().start(2,"fieldUpdate");
+
+        if ( runTime.timeOutputValue() - (recTimeIndex+1)*recTimeStep_ + 1.0e-5 > 0.0 )
+        {
+            Info<< "Updating fields at run time " << runTime.timeOutputValue()
+                << " corresponding to recurrence time " << (recTimeIndex+1)*recTimeStep_ << ".\n" << endl;
+            recBases[0].updateRecFields();
+            #include "readFields.H"
+
+            recTimeIndex++;
+        }
+
+        myClock().stop("fieldUpdate");
+
+        #include "continuityErrCalc.H"
+
+        myClock().start(3,"speciesEqn");
+        #include "CEq.H"
+        myClock().stop("speciesEqn");
+
+        myClock().stop("Total");
+
+        runTime.write();
+
+        Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
+            << "  ClockTime = " << runTime.elapsedClockTime() << " s"
+            << nl << endl;
+
+        myClock().stop("Global");
+
+    }
+
+
+    myClock().evalPar();
+    myClock().normHist();
+
+    Info<< "End\n" << endl;
+
+    return 0;
+}
+
+
+// ************************************************************************* //

applications/solvers/rctfSpeciesTransport/readFields.H

@@ -0,0 +1,40 @@
+
+currTimeIndex = recBases[0].recM().currentTimeIndex();
+
+Info << "current Time Index = " << currTimeIndex << endl;
+
+recBases[0].recM().exportVolVectorField(Ucoh_pair,UcohRec);
+recBases[0].recM().exportVolScalarField(kSGScoh_pair,kcohRec);
+
+
+label incTimeIndex = incPairTimeIndex_[currTimeIndex][1];
+
+Info << "  incoherent pair Time Index = " << incTimeIndex  << endl;
+
+UincRec = recBases[1].recM().exportVolVectorField(Uinc_pair,incTimeIndex);
+kincRec = recBases[1].recM().exportVolScalarField(kSGSinc_pair,incTimeIndex);
+
+kRec = kcohRec+kincRec;
+
+forAll(kRec, cellI)
+{
+    if (kRec[cellI] < SMALL)
+    {
+        kRec[cellI] = 0.0;
+    }
+}
+
+const fvPatchList& patches = mesh.boundary();
+forAll(patches, patchI)
+{
+        kRec.boundaryFieldRef()[patchI] = 0.0;
+}
+
+URec = UcohRec + UincRec;
+phiRec = linearInterpolate(URec) & mesh.Sf();
+
+nutRec = sqrt(kRec)*delta*0.094;
+
+alphat = nutRec/Sct;
+alphat.correctBoundaryConditions();
+

applications/solvers/recSolverTurbTransport/createFields.H

@@ -147,7 +147,6 @@ surfaceScalarField phiRec
     
     fvScalarMatrix TEqn(T, dimless*dimVolume/(dimTime));
 
-    scalar relaxCoeff(0.0);
 
  Info<< "reading clockProperties\n" << endl;
 

applications/solvers/recSolverTurbTransport/recSolverTurbTransport.C

@@ -1,6 +1,6 @@
 /*---------------------------------------------------------------------------*\
     CFDEMcoupling academic - Open Source CFD-DEM coupling
-    
+
     Contributing authors:
     Thomas Lichtenegger, Gerhard Holzinger
     Copyright (C) 2015- Johannes Kepler University, Linz
@@ -29,9 +29,9 @@ Description
     for a solver based on recurrence statistics
 
 Rules
-	Solution data to compute the recurrence statistics from, needs to 
-		reside in $CASE_ROOT/dataBase
-	Time step data in dataBase needs to be evenly spaced in time
+    Solution data to compute the recurrence statistics from, needs to
+    reside in $CASE_ROOT/dataBase
+    Time step data in dataBase needs to be evenly spaced in time
 
 \*---------------------------------------------------------------------------*/
 
@@ -55,34 +55,40 @@ int main(int argc, char *argv[])
     #include "createControl.H"
     #include "createFields.H"
     #include "createFvOptions.H"
-  
+    scalar relaxCoeff(0.0);
+
     recBase recurrenceBase(mesh);
 
     // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
     Info<< "\nCalculating particle trajectories based on recurrence statistics\n" << endl;
-    
+
     label recTimeIndex(0);
-    scalar recTimeStep_=recurrenceBase.recM().recTimeStep();
-    
+    label stepCounter = 0;
+    label recTimeStep2CFDTimeStep = recurrenceBase.recM().recTimeStep2CFDTimeStep();
+
     while (runTime.run())
     {
 
         myClock().start(1,"Global");
 
         runTime++;
-        
+
         Info<< "Time = " << runTime.timeName() << nl << endl;
 
         myClock().start(2,"fieldUpdate");
-        
-        if ( runTime.timeOutputValue() - (recTimeIndex+1)*recTimeStep_ + 1.0e-5 > 0.0 )
+
+        stepCounter++;
+
+        if (stepCounter == recTimeStep2CFDTimeStep)
         {
-	    Info << "Updating fields at run time " << runTime.timeOutputValue()
-	        << " corresponding to recurrence time " << (recTimeIndex+1)*recTimeStep_ << ".\n" << endl;
+            Info<< "Updating fields at run time " << runTime.timeOutputValue()
+                << " with recTimeIndex " << recTimeIndex << ".\n" << endl;
             recurrenceBase.updateRecFields();
             #include "readFields.H"
             recTimeIndex++;
+            stepCounter = 0;
+            recTimeStep2CFDTimeStep = recurrenceBase.recM().recTimeStep2CFDTimeStep();
         }
 
         myClock().stop("fieldUpdate");
@@ -92,15 +98,15 @@ int main(int argc, char *argv[])
         myClock().stop("speciesEqn");
 
         runTime.write();
-        
+
         Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
-        << "  ClockTime = " << runTime.elapsedClockTime() << " s"
-        << nl << endl;
-        
+            << "  ClockTime = " << runTime.elapsedClockTime() << " s"
+            << nl << endl;
+
         myClock().stop("Global");
 
     }
-    
+
     myClock().evalPar();
     myClock().normHist();
 

applications/solvers/rtfmSolverSpecies/createFields.H

@@ -132,4 +132,3 @@
 
     T.write();
 
-    scalar relaxCoeff(0.0);

applications/solvers/rtfmSolverSpecies/rtfmSolverSpecies.C

@@ -58,16 +58,18 @@ int main(int argc, char *argv[])
     #include "createControl.H"
     #include "createFields.H"
     #include "createFvOptions.H"
+    scalar relaxCoeff(0.0);
 
     cfdemCloudRec<cfdemCloud> particleCloud(mesh);
     recBase recurrenceBase(mesh);
 
     // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-    Info << "\nCalculating particle trajectories based on recurrence statistics\n" << endl;
+    Info<< "\nCalculating particle trajectories based on recurrence statistics\n" << endl;
 
     label recTimeIndex(0);
-    scalar recTimeStep_=recurrenceBase.recM().recTimeStep();
+    label stepCounter = 0;
+    label recTimeStep2CFDTimeStep = recurrenceBase.recM().recTimeStep2CFDTimeStep();
 
     while (runTime.run())
     {
@@ -76,20 +78,23 @@ int main(int argc, char *argv[])
         // do stuff (every lagrangian time step)
         particleCloud.clockM().start(1,"Global");
 
-        Info << "Time = " << runTime.timeName() << nl << endl;
+        Info<< "Time = " << runTime.timeName() << nl << endl;
 
         particleCloud.clockM().start(2,"Flow");
         #include "TEq.H"
         particleCloud.clockM().stop("Flow");
 
+        stepCounter++;
 
-        if ( runTime.timeOutputValue() - (recTimeIndex+1)*recTimeStep_ + 1.0e-5 > 0.0 )
+        if (stepCounter == recTimeStep2CFDTimeStep)
         {
-            Info << "Updating fields at run time " << runTime.timeOutputValue()
-                 << " corresponding to recurrence time " << (recTimeIndex+1)*recTimeStep_ << ".\n" << endl;
+            Info<< "Updating fields at run time " << runTime.timeOutputValue()
+                << " corresponding to recTimeIndex " << recTimeIndex << ".\n" << endl;
             recurrenceBase.updateRecFields();
             #include "readFields.H"
             recTimeIndex++;
+            stepCounter = 0;
+            recTimeStep2CFDTimeStep = recurrenceBase.recM().recTimeStep2CFDTimeStep();
         }
 
         particleCloud.clockM().start(27,"Output");
@@ -98,12 +103,12 @@ int main(int argc, char *argv[])
 
         particleCloud.clockM().stop("Global");
 
-        Info << "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
-             << "  ClockTime = " << runTime.elapsedClockTime() << " s"
-             << nl << endl;
+        Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
+            << "  ClockTime = " << runTime.elapsedClockTime() << " s"
+            << nl << endl;
     }
 
-    Info << "End\n" << endl;
+    Info<< "End\n" << endl;
 
     return 0;
 }

applications/utilities/displacementField/Make/files

@@ -0,0 +1,3 @@
+displacementField.C
+
+EXE = $(CFDEM_APP_DIR)/displacementField

applications/utilities/displacementField/Make/options

@@ -0,0 +1,7 @@
+EXE_INC = \
+    -I$(LIB_SRC)/finiteVolume/lnInclude \
+    -I$(LIB_SRC)/meshTools/lnInclude
+
+EXE_LIBS = \
+    -lfiniteVolume \
+    -lmeshTools

applications/utilities/displacementField/displacementField.C

@@ -0,0 +1,495 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     | Website:  https://openfoam.org
+    \\  /    A nd           | Copyright (C) 2011-2018 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/>.
+
+Application
+    displacementField
+
+\*---------------------------------------------------------------------------*/
+
+#include "fvCFD.H"
+#include "vectorList.H"
+#include <sstream>
+#include <string>
+#include <sys/stat.h>
+#include <unistd.h>
+#include <set>
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+void findPairs(labelList &, labelList &, labelPairList &);
+void findPairsUnordered(labelList &, labelList &, labelPairList &);
+void fillEmptyCells(fvMesh &, label, label, labelList &, volVectorField &, volVectorField &, scalarList &, volVectorField &, volVectorField &, bool, scalar);
+void nearestNeighborCells(fvMesh &, label, label, label, labelList &, labelList &);
+void normalizeFields(labelList &, volVectorField &, volVectorField &);
+void readDump(std::string, labelList &, vectorList &);
+scalar weightFun(scalar);
+
+int main(int argc, char *argv[])
+{
+    argList::addOption
+    (
+        "totalProcs",
+        "label",
+        "total number of parallel processes, defaults to 1"
+    );
+    argList::addOption
+    (
+        "thisProc",
+        "label",
+        "number of current process, defaults to 0"
+    );
+
+
+    #include "setRootCase.H"
+    #include "createTime.H"
+    #include "createMesh.H"
+
+    const label thisProc = args.optionLookupOrDefault("thisProc", 0);
+    const label totalProcs = args.optionLookupOrDefault("totalProcs", 1);
+
+    Info << "This is number " << thisProc << " of " << totalProcs << " processes." << endl;
+
+
+    // user-defined input for each case
+    IOdictionary displacementProperties
+    (
+        IOobject
+        (
+            "displacementProperties",
+            mesh.time().constant(),
+            mesh,
+            IOobject::MUST_READ,
+            IOobject::NO_WRITE
+        )
+    );
+
+    label dumpIndexStart(readLabel(displacementProperties.lookup("dumpIndexStart")));
+    label dumpIndexEnd(readLabel(displacementProperties.lookup("dumpIndexEnd")));
+    label dumpIndexInputIncrement(readLabel(displacementProperties.lookup("dumpIndexInputIncrement")));
+    label dumpIndexDisplacementIncrement(readLabel(displacementProperties.lookup("dumpIndexDisplacementIncrement")));
+    label nNeighMin(readLabel(displacementProperties.lookup("nNeighMin")));
+    label maxSearchLayers(displacementProperties.lookupOrDefault<label>("maxSearchLayers",0));
+    scalar timePerInputStep(readScalar(displacementProperties.lookup("timePerInputStep")));
+    scalar timePerDisplacementStep(readScalar(displacementProperties.lookup("timePerDisplacementStep")));
+    scalar startTime(readScalar(displacementProperties.lookup("startTime")));
+    std::string filepath=string(displacementProperties.lookup("filepath"));
+    std::string fileext=string(displacementProperties.lookupOrDefault<string>("fileextension",""));
+    bool interpolate=bool(displacementProperties.lookupOrDefault<bool>("fillEmptyCells",true));
+    bool averageMode=bool(displacementProperties.lookupOrDefault<bool>("averageMode",false));
+    
+    volVectorField defaultUs
+    (
+        IOobject
+        (
+            "defaultUDisp",
+            runTime.timeName(),
+            mesh,
+            IOobject::READ_IF_PRESENT,
+            IOobject::NO_WRITE
+        ),
+        mesh,
+        dimensionedVector("zero", dimensionSet(0,1,-1,0,0), vector::zero)
+    );
+
+    volVectorField defaultUsDirectedVariance
+    (
+        IOobject
+        (
+            "defaultUDispDirectedVariance",
+            runTime.timeName(),
+            mesh,
+            IOobject::READ_IF_PRESENT,
+            IOobject::NO_WRITE
+        ),
+        mesh,
+        dimensionedVector("zero", dimensionSet(0,1,-1,0,0), vector::zero)
+    );
+
+    scalar xmin=scalar(displacementProperties.lookupOrDefault<scalar>("xmin",-1e10));
+    scalar xmax=scalar(displacementProperties.lookupOrDefault<scalar>("xmax",1e10));
+    scalar ymin=scalar(displacementProperties.lookupOrDefault<scalar>("ymin",-1e10));
+    scalar ymax=scalar(displacementProperties.lookupOrDefault<scalar>("ymax",1e10));
+    scalar zmin=scalar(displacementProperties.lookupOrDefault<scalar>("zmin",-1e10));
+    scalar zmax=scalar(displacementProperties.lookupOrDefault<scalar>("zmax",1e10));
+    scalarList boundaries(6);
+    boundaries[0]=xmin;
+    boundaries[1]=xmax;
+    boundaries[2]=ymin;
+    boundaries[3]=ymax;
+    boundaries[4]=zmin;
+    boundaries[5]=zmax;
+
+    vectorList probePoints=vectorList(displacementProperties.lookupOrDefault<vectorList>("probePoints",vectorList(0)));
+    bool monitorProbes = false;
+    if (probePoints.size()>0) monitorProbes = true;
+    #include "OFstream.H"
+    OFstream monitoringDataFile("monitoringData.txt");
+    if (monitorProbes)
+    {
+        monitoringDataFile << "# monitoring data file" << endl;
+        monitoringDataFile << "# format: time nPerCell[p1] UDisp[p1] UDispDirectedVariance[p1] nPerCell[p2] ... " << endl;
+        for(label p=0;p<probePoints.size();p++)
+        {
+            vector pos = probePoints[p];
+            monitoringDataFile << "# point[" << p << "] = " << pos << endl;
+        }
+    }
+
+
+    label dumpIndex1 = dumpIndexStart + thisProc * dumpIndexInputIncrement;
+    label dumpIndex2 = dumpIndex1 + dumpIndexDisplacementIncrement;
+
+    volVectorField Us
+    (
+        IOobject
+        (
+            "UDisp",
+            runTime.timeName(),
+            mesh,
+            IOobject::NO_READ,
+            IOobject::NO_WRITE
+        ),
+        mesh,
+        dimensionedVector("zero", dimensionSet(0,1,-1,0,0), vector::zero)
+    );
+
+    volVectorField UsDirectedVariance
+    (
+        IOobject
+        (
+            "UDispDirectedVariance",
+            runTime.timeName(),
+            mesh,
+            IOobject::NO_READ,
+            IOobject::NO_WRITE
+        ),
+        mesh,
+        dimensionedVector("zero", dimensionSet(0,1,-1,0,0), vector::zero)
+    );
+
+    labelList particlesInCell(mesh.nCells(), 0);
+
+    scalar currTime=startTime + thisProc * timePerInputStep;
+    label timeIndex=thisProc;
+
+    while(true)
+    {
+        runTime.setTime(currTime,timeIndex);
+        // read dump files and check which particle indices are present in both
+        labelList indices1, indices2;
+        vectorList positions1, positions2;
+
+        std::stringstream ss;
+        ss << filepath << dumpIndex1 << fileext;
+        std::string filename1 = ss.str();
+        ss.str("");
+        ss << filepath << dumpIndex2 << fileext;
+        std::string filename2 = ss.str();
+
+        if (access( filename1.c_str(), F_OK ) == -1 || access( filename2.c_str(), F_OK ) == -1 || dumpIndex2 > dumpIndexEnd)
+        {
+            if (averageMode)
+            {
+                normalizeFields(particlesInCell, Us, UsDirectedVariance);
+                fillEmptyCells(mesh,nNeighMin,maxSearchLayers,particlesInCell,Us,UsDirectedVariance,boundaries,defaultUs,defaultUsDirectedVariance,interpolate,timePerDisplacementStep);
+
+                Us /= timePerDisplacementStep;
+                UsDirectedVariance /= timePerDisplacementStep;
+                Us.write();
+                UsDirectedVariance.write();
+            }
+            break;
+        }
+
+        Info << "\nReading" << endl;
+        Info << "\t" << filename1 << endl;
+        Info << "\t" << filename2 << endl;
+        Info << "corresponding to time = " << currTime << "." << endl;
+
+        readDump(filename1, indices1, positions1);
+        readDump(filename2, indices2, positions2);
+
+        labelPairList pairs;
+        findPairs(indices1,indices2,pairs);
+
+        // average particle displacements and their variance
+        Info << "Binning particle displacements on mesh." << endl;
+        vector position, displacement;
+        label line1, line2;
+        label cellI;
+
+        if (!averageMode)
+        {
+            Us *= 0.0;
+            UsDirectedVariance *= 0.0;
+            particlesInCell.clear();
+            particlesInCell.setSize(mesh.nCells(), 0);
+        }
+
+        for (label partI = 0; partI < pairs.size(); partI++)
+        {
+            line1 = pairs[partI].first();
+            line2 = pairs[partI].second();
+            position = positions1[line1];
+            displacement = positions2[line2] - positions1[line1];
+            cellI = mesh.findCell(position);
+            if (cellI < 0)  continue;
+            particlesInCell[cellI] += 1;
+            Us[cellI] += displacement;
+
+            for (label comp=0;comp<3;comp++)
+            {
+                UsDirectedVariance[cellI].component(comp) += displacement.component(comp)*displacement.component(comp);
+            }
+        }
+
+        if (!averageMode)
+        {
+            normalizeFields(particlesInCell, Us, UsDirectedVariance);
+            fillEmptyCells(mesh,nNeighMin,maxSearchLayers,particlesInCell,Us,UsDirectedVariance,boundaries,defaultUs,defaultUsDirectedVariance,interpolate,timePerDisplacementStep);
+
+            Us /= timePerDisplacementStep;
+            UsDirectedVariance /= timePerDisplacementStep;
+            Us.write();
+            UsDirectedVariance.write();
+        }
+
+        if (averageMode && monitorProbes)
+        {
+            monitoringDataFile << currTime << " ";
+            for(label p=0;p<probePoints.size();p++)
+            {
+                vector pos = probePoints[p];
+                label cellP = mesh.findCell(pos);
+                monitoringDataFile << " " << particlesInCell[cellP] << " " << Us[cellP]/timePerDisplacementStep << " " << UsDirectedVariance[cellP]/(timePerDisplacementStep*timePerDisplacementStep);
+            }
+            monitoringDataFile << endl;
+        }
+
+        dumpIndex1 += dumpIndexInputIncrement*totalProcs;
+        dumpIndex2 += dumpIndexInputIncrement*totalProcs;
+        currTime += timePerInputStep*totalProcs;
+        timeIndex += totalProcs;
+    }
+    return 0;
+}
+
+void readDump(std::string filename, labelList &indices, vectorList &positions)
+{
+    #include <fstream>
+
+    const label leadingLines = 9;
+    label lineCounter = 0;
+    label partIndex;
+    scalar x, y, z;
+
+    indices.clear();
+    positions.clear();
+
+    std::ifstream file(filename);
+    std::string str;
+    while (std::getline(file, str))
+    {
+        if (lineCounter >= leadingLines)
+        {
+            sscanf(str.c_str(), "%d %lf %lf %lf", &partIndex, &x, &y, &z);
+            indices.append(partIndex);
+            positions.append(vector(x,y,z));
+        }
+        lineCounter++;
+    }
+}
+
+void findPairs(labelList &indices1, labelList &indices2, labelPairList &pairs)
+{
+    // remove all entries from first list if they are not present in second list
+    // this assumes ordered entries
+
+    if (indices2.size() == 0) return;
+
+    for (label i=0;i<indices1.size();i++)
+    {
+        label j1 = -1;
+        label j2 = indices2.size();
+        label jmid = 0;
+        label index1 = indices1[i];
+        while(true)
+        {
+            jmid = (j1+j2)/2;
+            if (indices2[jmid] > index1) j2 = jmid;
+            else if (indices2[jmid] < index1) j1 = jmid;
+            else
+            {
+                pairs.append(labelPair(i,jmid));
+                break;
+            }
+            if (j2-j1 == 1) break;
+        }
+    }
+    Info << "findPairs: " << pairs.size() << " pairs found." << endl;
+}
+
+void findPairsUnordered(labelList &indices1, labelList &indices2, labelPairList &pairs)
+{
+    // remove all entries from first list if they are not present in second list
+
+    for (label i=0;i<indices1.size();i++)
+    {
+        for (label j=0;j<indices2.size();j++)
+        {
+            if (indices1[i] == indices2[j])
+            {
+                pairs.append(labelPair(i,j));
+                break;
+            }
+        }
+    }
+    Info << "findPairs: " << pairs.size() << " pairs found." << endl;
+}
+
+void fillEmptyCells(fvMesh &mesh, label nNeighMin, label maxSearchLayers, labelList &particlesInCell, volVectorField &Us, volVectorField& UsDirectedVariance,scalarList& boundaries, volVectorField &defaultUs, volVectorField &defaultUsDirectedVariance, bool interpolate, scalar dt)
+{
+    labelList neighborsWithValues;
+    scalar neighborSqrDistance;
+    scalar weight;
+    scalar weightSum;
+    scalarList weights;
+
+    Info << "Filling empty cells." << endl;
+    forAll(mesh.C(), cellI)
+    {
+        if (particlesInCell[cellI] > 0) continue;
+
+        vector position = mesh.C()[cellI];
+        label outsideBox = 0;
+        if (position.x() < boundaries[0] || position.x() > boundaries[1]) outsideBox++;
+        if (position.y() < boundaries[2] || position.y() > boundaries[3]) outsideBox++;
+        if (position.z() < boundaries[4] || position.z() > boundaries[5]) outsideBox++;
+
+        if (outsideBox > 0 || !interpolate)
+        {
+            Us[cellI] = defaultUs[cellI]*dt;
+            UsDirectedVariance[cellI] = defaultUsDirectedVariance[cellI]*dt;
+            continue;
+        }
+
+        nearestNeighborCells(mesh, cellI, nNeighMin, maxSearchLayers, particlesInCell, neighborsWithValues);
+        weightSum = 0.0;
+        weights.clear();
+        for (label neighI=0; neighI<neighborsWithValues.size(); neighI++)
+        {
+            neighborSqrDistance = magSqr(mesh.C()[cellI] - mesh.C()[neighborsWithValues[neighI]]);
+            weight = weightFun(neighborSqrDistance);
+            weights.append(weight);
+            weightSum += weight;
+        }
+        for (label neighI=0; neighI<neighborsWithValues.size(); neighI++)
+        {
+            weight = weights[neighI]/weightSum;
+            Us[cellI] += weight*Us[neighborsWithValues[neighI]];
+            UsDirectedVariance[cellI] += weight*UsDirectedVariance[neighborsWithValues[neighI]];
+        }
+        
+        if (neighborsWithValues.size() == 0)
+        {
+            Us[cellI] = defaultUs[cellI]*dt;
+            UsDirectedVariance[cellI] = defaultUsDirectedVariance[cellI]*dt;
+        }
+
+        // make sure no particles are placed outside of domain
+        // TODO: correct following implementation (meshSearch) and test it
+/*
+        vector shiftedPosition = position + dt * Us[cellI];
+        label cellJ = mesh.findCell(shiftedPosition);
+        if (cellJ < 0)
+        {
+            label cellK = mesh.findNearestCellWalk(shiftedPosition,cellI);
+            Us[cellI] = (mesh.C()[cellI] - mesh.C()[cellK]) / dt;
+        }
+*/
+    }
+}
+
+void nearestNeighborCells(fvMesh &mesh, label refCell, label nNeighMin, label maxSearchLayers, labelList &particlesInCell, labelList &neighborsWithValues)
+{
+    label numSearchLayers = 0;
+    std::set<label> neighbors;
+    std::set<label> newNeighbors;
+    std::set<label> recentNeighbors;
+
+    neighbors.insert(refCell);
+    recentNeighbors.insert(refCell);
+
+    neighborsWithValues.clear();
+
+    while(neighborsWithValues.size() < nNeighMin)
+    {
+        for (std::set<label>::iterator it=recentNeighbors.begin(); it!=recentNeighbors.end(); ++it)
+        {
+            labelList adjacent = mesh.cellCells()[*it];
+            label adj;
+            for (label j=0; j<adjacent.size(); j++)
+            {
+                adj = adjacent[j];
+                std::set<label>::iterator it2 = neighbors.find(adj);
+                if (it2 == neighbors.end())
+                {
+                    newNeighbors.insert(adj);
+                    neighbors.insert(adj);
+                    if (particlesInCell[adj] > 0) neighborsWithValues.append(adj);
+                }
+            }
+        }
+        
+        numSearchLayers++;
+        if (numSearchLayers > maxSearchLayers && maxSearchLayers > 0) return;
+
+        if (newNeighbors.size() == 0) return;
+        recentNeighbors.clear();
+        recentNeighbors = newNeighbors;
+        newNeighbors.clear();
+    }
+}
+
+void normalizeFields(labelList& particlesInCell, volVectorField& Us, volVectorField & UsDirectedVariance)
+{
+    for (label cellJ = 0; cellJ<particlesInCell.size(); cellJ++)
+    {
+        if (particlesInCell[cellJ] > 0)
+        {
+            Us[cellJ] /= particlesInCell[cellJ];
+            UsDirectedVariance[cellJ] /= particlesInCell[cellJ];
+            for (label comp=0;comp<3;comp++)
+            {
+                UsDirectedVariance[cellJ].component(comp) -= Us[cellJ].component(comp)*Us[cellJ].component(comp);
+                if (UsDirectedVariance[cellJ].component(comp) > 0) UsDirectedVariance[cellJ].component(comp) = Foam::sqrt(UsDirectedVariance[cellJ].component(comp));
+            }
+        }
+    }
+}
+
+scalar weightFun(scalar distSqr)
+{
+    // inverse distance weighting, order 2
+    return 1.0/distSqr;
+}
+// ************************************************************************* //

applications/utilities/rBaseMirror/rBaseMirror/Make/files

@@ -0,0 +1,3 @@
+rBaseMirror.C
+
+EXE=$(CFDEM_APP_DIR)/rBaseMirror

applications/utilities/rBaseMirror/rBaseMirror/createFields.H

@@ -13,5 +13,7 @@
     vector refPoint(mirrorProperties.lookup("refPoint"));
     vector refDirection(mirrorProperties.lookup("refDirection"));
 
-    word fieldName(mirrorProperties.lookup("fieldName"));
     word dataBaseName(mirrorProperties.lookup("dataBaseName"));
+
+    const wordList volScalarFieldNames(mirrorProperties.lookup("volScalarFields"));
+    const wordList volVectorFieldNames(mirrorProperties.lookup("volVectorFields"));

applications/utilities/rBaseMirror/rBaseMirror/rBaseMirror.C

@@ -66,16 +66,8 @@ int main(int argc, char *argv[])
     instantList timeDirs(recTime.times());
     recTime.setTime(timeDirs[0],0);
 
-    #include "readFields.H"
-
-    Info << fieldName << endl;
-
-    volScalarField transformedField = origField;
-
     scalar t;
 
-    label shiftedTimeI = 0;
-
     // check number of time directories
     label shift = 0;
     forAll(timeDirs, timeI)
@@ -95,37 +87,94 @@ int main(int argc, char *argv[])
     label cellI_transformed = -1;
     forAll(timeDirs, timeI)
     {
-
        recTime.setTime(timeDirs[timeI], timeI);
        t = recTime.value();
        if(t < startTime) continue;
        if(t > endTime) continue;
        Info << "time = " << t << ", time index = " << timeI << endl;
 
-       #include "readFields.H"
-
-       forAll(transformedField, cellI)
+       // volScalarFields
+       for (int sf = 0; sf < volScalarFieldNames.size(); sf++)
        {
-           vector position = mesh.C()[cellI];
-           vector transformedPosition = 2 * ((refPoint - position) & refDirection) * refDirection / (refDirection & refDirection) + position;
-           cellI_transformed = mesh.findCell(transformedPosition);
-           if(cellI_transformed < 0)
+           word fieldName = volScalarFieldNames[sf];
+           volScalarField origField
+           (
+               IOobject
+               (
+                    fieldName,
+                    recTime.timePath(),
+                    mesh,
+                    IOobject::MUST_READ,
+                    IOobject::NO_WRITE
+               ),
+                mesh
+           );
+
+           volScalarField transformedField = origField;
+
+           forAll(transformedField, cellI)
            {
-               Info << "Couldn't find transformed cell. Stopping." << endl;
-               return 0;
+               vector position = mesh.C()[cellI];
+               vector transformedPosition = 2 * ((refPoint - position) & refDirection) * refDirection / (refDirection & refDirection) + position;
+               cellI_transformed = mesh.findCell(transformedPosition);
+               if(cellI_transformed < 0)
+               {
+                   Info << "Couldn't find transformed cell. Stopping." << endl;
+                   return 0;
+               }
+
+               scalar value = origField[cellI_transformed];
+               scalar transformedValue = value;
+
+               transformedField[cellI] = transformedValue;
            }
 
-           scalar value = origField[cellI_transformed];
-           scalar transformedValue = value;
-
-           transformedField[cellI] = transformedValue;
+           runTime.setTime(recTime.value() + origTimeRange + dt, timeI + shift);
+           Info << "creating transformed field " << fieldName << " for time = " << recTime.value() + origTimeRange + dt << endl;
+           transformedField.write();
        }
 
-       shiftedTimeI = timeI + shift;
-       t = recTime.value() + origTimeRange + dt;
-       runTime.setTime(t, shiftedTimeI);
-       Info << "creating transformed fields for time = " << t << ", time index = " << shiftedTimeI << endl;
-       transformedField.write();
+
+       // volVectorFields
+       for (int vf = 0; vf < volVectorFieldNames.size(); vf++)
+       {
+           word fieldName = volVectorFieldNames[vf];
+           volVectorField origField
+           (
+               IOobject
+               (
+                    fieldName,
+                    recTime.timePath(),
+                    mesh,
+                    IOobject::MUST_READ,
+                    IOobject::NO_WRITE
+               ),
+                mesh
+           );
+
+           volVectorField transformedField = origField;
+
+           forAll(transformedField, cellI)
+           {
+               vector position = mesh.C()[cellI];
+               vector transformedPosition = 2 * ((refPoint - position) & refDirection) * refDirection / (refDirection & refDirection) + position;
+               cellI_transformed = mesh.findCell(transformedPosition);
+               if(cellI_transformed < 0)
+               {
+                   Info << "Couldn't find transformed cell. Stopping." << endl;
+                   return 0;
+               }
+
+               vector value = origField[cellI_transformed];
+               vector transformedValue = -2 * (value & refDirection) * refDirection / (refDirection & refDirection) + value;
+
+               transformedField[cellI] = transformedValue;
+           }
+
+           runTime.setTime(recTime.value() + origTimeRange + dt, timeI + shift);
+           Info << "creating transformed field " << fieldName << " for time = " << recTime.value() + origTimeRange + dt << endl;
+           transformedField.write();
+        }
     }
 
     Info << "\nEnd" << endl;

applications/utilities/rBaseMirror/rBaseMirrorScalar/Make/files

@@ -1,3 +0,0 @@
-rBaseMirrorScalar.C
-
-EXE=$(CFDEM_APP_DIR)/rBaseMirrorScalar

applications/utilities/rBaseMirror/rBaseMirrorScalar/createFields.H

@@ -1,17 +0,0 @@
-    IOdictionary mirrorProperties
-    (
-        IOobject
-        (
-            "mirrorProperties",
-            mesh.time().constant(),
-            mesh,
-            IOobject::MUST_READ,
-            IOobject::NO_WRITE
-        )
-    );
-
-    vector refPoint(mirrorProperties.lookup("refPoint"));
-    vector refDirection(mirrorProperties.lookup("refDirection"));
-
-    word fieldName(mirrorProperties.lookup("fieldName"));
-    word dataBaseName(mirrorProperties.lookup("dataBaseName"));

applications/utilities/rBaseMirror/rBaseMirrorScalar/readFields.H

@@ -1,13 +0,0 @@
-
-    volScalarField origField
-    (
-        IOobject
-        (
-            fieldName,
-            recTime.timePath(),
-            mesh,
-            IOobject::MUST_READ,
-            IOobject::NO_WRITE
-        ),
-        mesh
-    );

applications/utilities/rBaseMirror/rBaseMirrorVec/Make/files

@@ -1,3 +0,0 @@
-rBaseMirrorVec.C
-
-EXE=$(CFDEM_APP_DIR)/rBaseMirrorVec

applications/utilities/rBaseMirror/rBaseMirrorVec/Make/options

@@ -1,17 +0,0 @@
-include $(CFDEM_ADD_LIBS_DIR)/additionalLibs
-
-
-EXE_INC = \
-     $(PFLAGS) \
-    -I$(LIB_SRC)/finiteVolume/cfdTools \
-    -I$(LIB_SRC)/finiteVolume/lnInclude \
-    -I$(LIB_SRC)/meshTools/lnInclude \
-    -I$(LIB_SRC)/sampling/lnInclude \
-    -I$(LIB_SRC)/fvOptions/lnInclude \
-    -Wno-deprecated-copy
-
-EXE_LIBS = \
-    -lfiniteVolume \
-    -lmeshTools \
-    -lsampling \
-    -lfvOptions

applications/utilities/rBaseMirror/rBaseMirrorVec/rBaseMirrorVec.C

@@ -1,134 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 1991-2009 OpenCFD Ltd.
-     \\/     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 2 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, write to the Free Software Foundation,
-    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
-
-Application
-    rBaseMirror
-
-Description
-    Read time series and extend it by mirrored fields if geometry possesses
-    the same symmetry
-
-\*---------------------------------------------------------------------------*/
-
-#include "fvCFD.H"
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-//  Main program:
-
-int main(int argc, char *argv[])
-{
-    argList::noParallel();
-    timeSelector::addOptions();
-
-    #include "setRootCase.H"
-    #include "createTime.H"
-
-   // read in start and end time from controlDict
-
-    scalar startTime=runTime.startTime().value();
-    scalar endTime=runTime.endTime().value();
-    scalar origTimeRange = endTime - startTime;
-
-    Info << "start time = " << runTime.startTime() << endl;
-    Info << "end time = " << runTime.endTime() << endl;
-
-    // check which time directories are present
-    //instantList timeDirs = timeSelector::select0(runTime, args);
-    //runTime.setTime(timeDirs[0], 0);
-
-    #include "createMesh.H"
-
-    #include "createFields.H"
-
-    Foam::Time recTime(fileName(dataBaseName), "", "../system", "../constant", false);
-    instantList timeDirs(recTime.times());
-    recTime.setTime(timeDirs[0],0);
-
-    #include "readFields.H"
-
-    Info << fieldName << endl;
-
-    volVectorField transformedField = origField;
-
-    scalar t;
-
-    label shiftedTimeI = 0;
-
-    // check number of time directories
-    label shift = 0;
-    forAll(timeDirs, timeI)
-    {
-
-       if (recTime.timeName() == "constant") continue;
-       recTime.setTime(timeDirs[timeI], timeI);
-       t = recTime.value();
-       if(t < startTime) continue;
-       if(t > endTime) continue;
-       shift++;
-    }
-
-    scalar dt = origTimeRange / (shift - 1.0);
-    recTime.setEndTime(startTime + 2 * origTimeRange + dt);
-
-    label cellI_transformed = -1;
-    forAll(timeDirs, timeI)
-    {
-       recTime.setTime(timeDirs[timeI], timeI);
-       t = recTime.value();
-       if(t < startTime) continue;
-       if(t > endTime) continue;
-       Info << "time = " << t << ", time index = " << timeI << endl;
-
-       #include "readFields.H"
-
-       forAll(transformedField, cellI)
-       {
-           vector position = mesh.C()[cellI];
-           vector transformedPosition = 2 * ((refPoint - position) & refDirection) * refDirection / (refDirection & refDirection) + position;
-           cellI_transformed = mesh.findCell(transformedPosition);
-           if(cellI_transformed < 0)
-           {
-               Info << "Couldn't find transformed cell. Stopping." << endl;
-               return 0;
-           }
-
-           vector value = origField[cellI_transformed];
-           vector transformedValue = -2 * (value & refDirection) * refDirection / (refDirection & refDirection) + value;
-
-           transformedField[cellI] = transformedValue;
-       }
-
-       shiftedTimeI = timeI + shift;
-       t = recTime.value() + origTimeRange + dt;
-       runTime.setTime(t, shiftedTimeI);
-       Info << "creating transformed fields for time = " << t << ", time index = " << shiftedTimeI << endl;
-       transformedField.write();
-    }
-
-    Info << "\nEnd" << endl;
-
-    return 0;
-}
-
-// ************************************************************************* //

applications/utilities/rBaseMirror/rBaseMirrorVec/readFields.H

@@ -1,14 +0,0 @@
-
-volVectorField origField
-    (
-        IOobject
-        (
-            fieldName,
-            recTime.timePath(),
-            mesh,
-            IOobject::MUST_READ,
-            IOobject::NO_WRITE
-        ),
-        mesh
-    );
-

doc/CFDEMcoupling_models.txt

@@ -37,6 +37,7 @@ The "chemistryModels"_chemistryModel.html keyword is used to specify a list of
 models used for chemical reaction calculations.
 
 "diffusionCoefficients"_chemistryModel_diffusionCoefficients.html,
+initMultiLayers,
 "massTransferCoeff"_chemistryModel_massTransferCoeff.html,
 "off"_chemistryModel_noChemistry.html,
 "reactantPerParticle"_chemistryModel_reactantPerParticle.html,
@@ -65,18 +66,28 @@ that performs the data exchange between the DEM code and the CFD code.
 "twoWayOne2One"_dataExchangeModel_twoWayOne2One.html :tb(c=2,ea=c)
 
 
-6.6 Energy models :h4
+6.6 Diffusion coefficient models :h4
+
+The {diffCoeffModel} keyword entry specifies the model for the diffusion
+coefficient of dissolved spieces in the fluid phase in the presence of particles.
+
+SyamlalDiffCoeff,
+off :tb(c=2,ea=c)
+
+
+6.7 Energy models :h4
 
 The {energyModels} keyword specifies a list of energy models used for e.g.
 compressible, reacting flows.
 
-heatTransferGranConduction,
 heatTransferGunn,
+heatTransferInterGrain,
 heatTransferRanzMarshall,
-reactionHeat :tb(c=2,ea=c)
+reactionHeat,
+wallHeatTransferYagi :tb(c=2,ea=c)
 
 
-6.7 Force models :h4
+6.8 Force models :h4
 
 The "forceModels"_forceModel.html keyword specifies a list of models that exert
 a force on each DEM particle.
@@ -86,15 +97,15 @@ a force on each DEM particle.
 "BeetstraDrag"_forceModel_BeetstraDrag.html,
 BeetstraDragPoly,
 "DiFeliceDrag"_forceModel_DiFeliceDrag.html,
-"dSauter"_forceModel_dSauter.html,
 Fines,
 "GidaspowDrag"_forceModel_GidaspowDrag.html,
 "KochHillDrag"_forceModel_KochHillDrag.html,
 "LaEuScalarTemp"_forceModel_LaEuScalarTemp.html,
 "MeiLift"_forceModel_MeiLift.html,
+"ParmarBassetForce"_forceModel_ParmarBassetForce.html,
 "SchillerNaumannDrag"_forceModel_SchillerNaumannDrag.html,
 "ShirgaonkarIB"_forceModel_ShirgaonkarIB.html,
-dSauter,
+"dSauter"_forceModel_dSauter.html,
 deactivateForce,
 directedDiffusiveRelaxation,
 evaluateFluctuations,
@@ -112,12 +123,15 @@ particleDeformation,
 "pdCorrelation"_forceModel_pdCorrelation.html,
 potentialRelaxation,
 "surfaceTensionForce"_forceModel_surfaceTensionForce.html,
+terminalVelocity,
+turbulentDispersion,
+turbulentVelocityFluctuations,
 "virtualMassForce"_forceModel_virtualMassForce.html,
 "viscForce"_forceModel_viscForce.html,
 "volWeightedAverage"_forceModel_volWeightedAverage.html :tb(c=2,ea=c)
 
 
-6.7.1 Force sub-models :h5
+6.8.1 Force sub-models :h5
 
 The "forceSubModels"_forceSubModel.html keyword specifies a list
 of models that hold settings for a force model.
@@ -127,7 +141,7 @@ ScaleForce,
 scaleForceBoundary :tb(c=2,ea=c)
 
 
-6.8 LIGGGHTS command models :h4,link(lcm)
+6.9 LIGGGHTS command models :h4,link(lcm)
 
 The "liggghtsCommandModels"_liggghtsCommandModel.html keyword specifies a list
 of models that execute LIGGGHTS commands within a CFD run.
@@ -138,7 +152,7 @@ of models that execute LIGGGHTS commands within a CFD run.
 "writeLiggghts"_liggghtsCommandModel_writeLiggghts.html :tb(c=2,ea=c)
 
 
-6.9 Locate models :h4
+6.10 Locate models :h4
 
 The "locateModel"_locateModel.html keyword entry specifies the model used to
 search the CFD mesh for the CFD cell corresponding to a given position.
@@ -150,7 +164,15 @@ search the CFD mesh for the CFD cell corresponding to a given position.
 "turboEngine"_locateModel_turboEngineSearch.html :tb(c=2,ea=c)
 
 
-6.10 Mesh motion models :h4
+6.11 Mass transfer models :h4
+
+The {massTransferModels} keyword specifies a list of mass transfer models used
+for evaluating species transfer between particles and fluids.
+
+massTransferGunn :tb(c=2,ea=c)
+
+
+6.12 Mesh motion models :h4
 
 The "meshMotionModel"_meshMotionModel.html keyword entry specifies the model
 used to manipulate the CFD mesh according to the DEM mesh motion.
@@ -158,7 +180,7 @@ used to manipulate the CFD mesh according to the DEM mesh motion.
 "noMeshMotion"_meshMotionModel_noMeshMotion.html :tb(c=2,ea=c)
 
 
-6.11 Momentum coupling models :h4
+6.13 Momentum coupling models :h4
 
 The "momCoupleModels"_momCoupleModel.html keyword specifies a list of models
 used for momentum exchange between DEM and CFD simulation
@@ -169,17 +191,18 @@ deactivateCouple,
 "off"_momCoupleModel_noCouple.html :tb(c=2,ea=c)
 
 
-6.12 Other force models :h4
+6.14 Other force models :h4
 
 The {otherForceModels} keyword specifies a list of models that exert a force on
 each DEM particle.
 
 expParticleForces,
 gravity,
+secondaryPhaseInducedBuoyancy,
 weightSecondaryPhase :tb(c=2,ea=c)
 
 
-6.13 Probe models :h4
+6.15 Probe models :h4
 
 The "probeModel"_probeModel.html keyword entry specifies the probing features in
 CFDEMcoupling simulations.
@@ -188,7 +211,7 @@ CFDEMcoupling simulations.
 "particleProbe"_probeModel_particleProbe.html :tb(c=2,ea=c)
 
 
-6.14 Region models :h4
+6.16 Region models :h4
 
 The "regionModel"_regionModel.html keyword entry specifies the model used to
 select a certain region for coupled simulations.
@@ -196,17 +219,18 @@ select a certain region for coupled simulations.
 "allRegion"_regionModel_allRegion.html :tb(c=2,ea=c)
 
 
-6.15 Smoothing models :h4
+6.17 Smoothing models :h4
 
 The "smoothingModel"_smoothingModel.html keyword entry specifies the model for
 smoothing the exchange fields.
 
+"constDiffAndTemporalSmoothing"_smoothingModel_constDiffAndTemporalSmoothing.html,
 "constDiffSmoothing"_smoothingModel_constDiffSmoothing.html,
 "off"_smoothingModel_noSmoothing.html,
 "temporalSmoothing"_smoothingModel_temporalSmoothing.html :tb(c=2,ea=c)
 
 
-6.16 Thermal conductivity models :h4
+6.18 Thermal conductivity models :h4
 
 The {thermCondModel} keyword entry specifies the model for the thermal
 conductivity of the fluid phase in the presence of particles.
@@ -216,7 +240,7 @@ ZehnerSchluenderThermCond,
 off :tb(c=2,ea=c)
 
 
-6.17 Void fraction models :h4
+6.19 Void fraction models :h4
 
 The "voidFractionModel"_voidFractionModel.html keyword entry specifies the model
 accounting for the volume of the particles in the CFD domain.
@@ -229,3 +253,4 @@ accounting for the volume of the particles in the CFD domain.
 off,
 trilinear :tb(c=2,ea=c)
 
+

doc/CFDEMcoupling_solvers.txt

@@ -11,6 +11,7 @@ This section lists all CFDEMcoupling solvers alphabetically.
 
 "cfdemSolverIB"_cfdemSolverIB.html,
 "cfdemSolverMultiphase"_cfdemSolverMultiphase.html,
+"cfdemSolverMultiphaseScalar"_cfdemSolverMultiphaseScalar.html,
 "cfdemSolverPiso"_cfdemSolverPiso.html,
 "cfdemSolverPisoScalar"_cfdemSolverPisoScalar.html,
 "cfdemSolverRhoPimple"_cfdemSolverRhoPimple.html,

doc/cfdemSolverMultiphaseScalar.txt

@@ -0,0 +1,68 @@
+<!-- HTML_ONLY -->
+<HEAD>
+  <META CHARSET="utf-8">
+</HEAD>
+<!-- END_HTML_ONLY -->
+"CFDEMproject Website"_lws - "Main Page"_main :c
+
+:link(lws,http://www.cfdem.com)
+:link(main,CFDEMcoupling_Manual.html)
+
+:line
+
+cfdemSolverMultiphaseScalar command :h3
+
+[Description:]
+
+<!-- HTML_ONLY -->
+"cfdemSolverMultiphaseScalar" is a coupled CFD-DEM solver using the CFDEMcoupling framework. Based on the OpenFOAM solver multiphaseInterFoam&reg;(*) it has functionality to simulate several fluids using the Volume of Fluid approach, coupled with the DEM code LIGGGHTS for solid particles. Additionally, it enable evaluation of temperature and dissolved species concentration fields.
+
+<!-- END_HTML_ONLY -->
+
+<!-- RST
+
+"cfdemSolverMultiphaseScalar" is a coupled CFD-DEM solver using the CFDEMcoupling framework. Based on the OpenFOAM solver multiphaseInterFoam\ |reg|\ (*) it has functionality to simulate several fluids using the Volume of Fluid approach, coupled with the DEM code LIGGGHTS for solid particles. Additionally, it enable evaluation of temperature and dissolved species concentration fields.
+
+
+.. |reg|    unicode:: U+000AE .. REGISTERED SIGN
+
+END_RST -->
+
+For more details, see "Vångö et al. (2018)"_#Vångö2018 and "Nijssen et al. (2021)"_#Nijssen2021.
+
+:line
+
+:link(Vångö2018)
+[(Vångö2018)] M. Vångö, S. Pirker, T. Lichtenegger. (2018):
+"Unresolved CFD-DEM modeling of multiphase flow in densely packed particle beds",
+Applied Mathematical Modelling
+
+:link(Nijssen2021)
+[(Nijssen2021)] T.M.J. Nijssen, J.A.M. Kuipers, J. van der Stel, A.T. Adema, K.A. Buist. (2021):
+"article title here",
+journal name here
+
+:line
+
+<!-- HTML_ONLY -->
+NOTE:
+(*) This offering is not approved or endorsed by OpenCFD Limited, producer and
+distributor of the OpenFOAM software via www.openfoam.com, and owner of the
+OPENFOAM&reg; and OpenCFD&reg; trade marks.
+OPENFOAM&reg; is a registered trade mark of OpenCFD Limited, producer and
+distributor of the OpenFOAM software via www.openfoam.com.
+<!-- END_HTML_ONLY -->
+
+<!-- RST
+
+.. note::
+
+   (*) This offering is not approved or endorsed by OpenCFD Limited, producer
+   and distributor of the OpenFOAM software via www.openfoam.com, and owner of
+   the OPENFOAM\ |reg| and OpenCFD\ |reg| trade marks.
+   OPENFOAM\ |reg| is a registered trade mark of OpenCFD Limited, producer and
+   distributor of the OpenFOAM software via www.openfoam.com.
+
+.. |reg|    unicode:: U+000AE .. REGISTERED SIGN
+
+END_RST -->

doc/energyModel_wallHeatTransferYagi.txt

@@ -0,0 +1,83 @@
+"CFDEMproject Website"_lws - "Main Page"_main :c
+
+:link(lws,http://www.cfdem.com)
+:link(main,CFDEMcoupling_Manual.html)
+
+:line
+
+energyModel wallHeatTransferYagi command :h3
+
+[Syntax:]
+
+Defined in "couplingProperties"_CFDEMcoupling_dicts.html#couplingProperties
+dictionary.
+
+energyModels
+(
+    wallHeatTransferYagi
+);
+wallHeatTransferYagiProps
+\{
+    QWallFluidName         "QWallFluid";
+    QWallFluidCoeffNAme    "QWallFluidCoeff";
+    wallTempName           "wallTemp";
+    tempFieldName          "T";
+    voidfractionFieldName  "voidfraction";
+    voidfractionMax        scalar1;
+    maxSource              scalar2;
+    velFieldName           "U";
+    densityFieldName       "rho";
+    implicit               switch1;
+    verbose                switch2;
+\} :pre
+
+{QWallFluid} = name of the wall-fluid heat transfer rate field :ulb,l
+{QwallFluidCoeff} = name of the wall-fluid heat transfer coefficient field :l
+{wallTemp} = name of the wall temperature field :l
+{T} = name of the fluid temperature field :l
+{voidfraction} = name of the finite volume void fraction field :l 
+{scalar1} = maximum void fraction for a cell to be considered packed :l 
+{scalar2} = (optional, default 1e30) maximum allowed source term :l
+{U} = name of the finite volume fluid velocity field :l
+{rho} = name of the fluid density field :l
+{switch1} = (optional, default true) activate to treat the fluid temperature implicitly :l
+{switch2} = (optional, default false) activate to write additional fields, and write cell-based data to the terminal :l
+:ule
+
+[Examples:]
+
+energyModels
+(
+    wallHeatTransferYagi
+);
+wallHeatTransferYagiProps
+\{
+    QWallFluidName         "QWallFluid";
+    QWallFluidCoeffNAme    "QWallFluidCoeff";
+    wallTempName           "wallTemp";
+    tempFieldName          "T";
+    voidfractionFieldName  "voidfraction";
+    voidfractionMax        0.5;
+    velFieldName           "U";
+    densityFieldName       "rho";
+    implicit               true;
+\} :pre
+
+[Description:]
+
+The energy model performs calculation of wall-to-bed heat transfer for a
+packed bed, based on the correlation of "Yagi and Wakao (1959)"_#Yagi1959.
+
+:link(Yagi1959)
+[(Yagi1959)] S. Yagi, N. Wakao. (1959):
+"Heat and Mass Transfer from Wall to Fluid in Packed Beds",
+AIChE Journal
+
+[Restrictions:]
+
+IMPORTANT NOTE: Model not validated!
+
+[Related commands:]
+
+energyModel
+

doc/forceModel_MeiLift.txt

@@ -19,19 +19,23 @@ forceModels
 MeiLiftProps
 \{
     velFieldName        "U";
-    useSecondOrderTerms;
-    treatForceExplicit  switch1;
-    verbose             switch2;
-    interpolation       switch3;
-    scalarViscosity     switch4;
+    useShearInduced     switch1;
+    useSpinInduced      switch2;
+    combineShearSpin    switch3;
+    treatForceExplicit  switch4;
+    verbose             switch5;
+    interpolation       switch6;
+    scalarViscosity     switch7;
 \} :pre
 
 {U} = name of the finite volume fluid velocity field :ulb,l
-{useSecondOrderTerms} = (optional, default false) switch to activate second order terms in the lift force model :l
-{switch1} = (optional, default false) sub model switch, see "forceSubModel"_forceSubModel.html for details :l
-{switch2} = (optional, default false) switch to activate the report of per-particle quantities to the screen :l
-{switch3} = (optional, default false) switch to activate tri-linear interpolation of the flow quantities at the particle position :l
+{useShearInduced} = (optional, default true) switch to activate shear-induced (Saffman) lift :l
+{useSpinInduced} = (optional, default false) switch to activate spin-induced (Magnus) lift :l
+{combineShearSpin} = (optional, default false) switch to use the correlation for combined shear- and spin-induced lift by "Loth and Dorgan (2009)"_#LothDorgan2009 :l
 {switch4} = (optional, default false) sub model switch, see "forceSubModel"_forceSubModel.html for details :l
+{switch5} = (optional, default false) switch to activate the report of per-particle quantities to the screen :l
+{switch6} = (optional, default false) switch to activate tri-linear interpolation of the flow quantities at the particle position :l
+{switch7} = (optional, default false) sub model switch, see "forceSubModel"_forceSubModel.html for details :l
 :ule
 
 [Examples:]
@@ -43,7 +47,6 @@ forceModels
 MeiLiftProps
 \{
     velFieldName        "U";
-    useSecondOrderTerms;
     interpolation       true;
     verbose             true;
 \} :pre
@@ -52,9 +55,12 @@ MeiLiftProps
 
 The force model performs the calculation of forces (e.g. fluid-particle
 interaction forces) acting on each DEM particle. The {MeiLift} model calculates
-the lift force for each particle based on Loth and Dorgan (2009). In case the
-keyword "useSecondOrderTerms" is not specified, this lift force model uses the
-expression of McLaughlin (1991, J. Fluid Mech. 224:261-274).
+the lift force for each particle based on "Loth and Dorgan (2009)"_#LothDorgan2009.
+
+:link(LothDorgan2009)
+[(LothDorgan2009)] E. Loth, A.J. Dorgan. (2009):
+"An equation of motion for particles of finite Reynolds number and size",
+Environmental Fluid Mechanics
 
 [Restrictions:]
 

doc/forceModel_ParmarBassetForce.txt

@@ -0,0 +1,88 @@
+"CFDEMproject Website"_lws - "Main Page"_main :c
+
+:link(lws,http://www.cfdem.com)
+:link(main,CFDEMcoupling_Manual.html)
+
+:line
+
+forceModel ParmarBassetForce command :h3
+
+[Syntax:]
+
+Defined in "couplingProperties"_CFDEMcoupling_dicts.html#couplingProperties
+dictionary.
+
+forceModels
+(
+    ParmarBassetForce
+);
+ParmarBassetForceProps
+\{
+    velFieldName        "U";
+    UsFieldName         "Us";
+    nIntegral           scalar1;
+    discretisationOrder scalar2;
+    treatForceExplicit  switch1;
+    interpolation       switch2;
+    smoothingModel      "smoothingModel";
+\} :pre
+
+{U} = name of the finite volume fluid velocity field :ulb,l
+{Us} = name of the finite volume cell averaged particle velocity field :l
+{scalar1} = number of timesteps considered in the near history :l
+{scalar2} = (1 or 2) discretisation order of the far history differential equations :l
+{switch1} = (optional, default true) sub model switch, see "forceSubModel"_forceSubModel.html for details :l
+{switch2} = (optional, default false) sub model switch, see "forceSubModel"_forceSubModel.html for details :l
+{smoothingModel} = name of smoothing model for the dU/dt field :l
+:ule
+
+[Examples:]
+
+forceModels
+(
+    ParmarBassetForce
+);
+ParmarBassetForceProps
+\{
+    velFieldName        "U";
+    USFieldName         "Us"
+    nIntegral           10;
+    discretisationOrder 1;
+    smoothingModel      constDiffAndTemporalSmoothing;
+    constDiffAndTemporalSmoothingProps
+    \{
+        lowerLimit          1e-8;
+        upperLimit          1e8;
+        smoothingLength     0.1;
+        smoothingStrength   0.001;
+        correctBoundary     true;
+    \}
+\} :pre
+
+[Description:]
+
+The force model performs the calculation of forces (e.g. fluid-particle
+interaction forces) acting on each DEM particle. The {ParmarBassetForce} model
+calculates the Basset history force for each particle, based on the method by
+"Parmar et al. (2018)"_#Parmar2018.
+
+For more detail, see "Nijssen et al. (2020)"_#Nijssen2020.
+
+:link(Parmar2018)
+[(Parmar2018)] M. Parmar, S. Annamalai, S. Balachandar, A. Prosperetti. (2018):
+"Differential Formulation of the Viscous History Force on a Particle for Efficient and Accurate Computation",
+Journal of Fluid Mechanics
+
+:link(Nijssen2020)
+[(Nijssen2020)] T.M.J. Nijssen, J.A.M. Kuipers, J. van der Stel, A.T. Adema, K.A. Buist. (2020):
+"Complete liquid-solid momentum coupling for unresolved CFD-DEM simulations",
+International Journal of Multiphase Flow
+
+[Restrictions:]
+
+IMPORTANT NOTE: Model not validated!
+
+[Related commands:]
+
+"forceModel"_forceModel.html
+

doc/forceModel_virtualMassForce.txt

@@ -19,19 +19,29 @@ forceModels
 virtualMassForceProps
 \{
     velFieldName "U";
+    voidfractionFieldName "voidfraction";
+    UsFieldName "Us";
     phiFieldName "phi";
     splitUrelCalculation switch1;
+    useUs switch2;
+    useFelderhof switch3;
     Cadd scalar1;
-    treatForceExplicit switch2;
-    interpolation switch3;
+    treatForceExplicit switch4;
+    interpolation switch5;
+    smoothingModel "smoothingModel";
 \} :pre
 
 {U} = name of the finite volume fluid velocity field :ulb,l
+{voidfraction} = name of the finite volume void fraction field :l
+{Us} = name of the finite volume cell averaged particle velocity field :l
 {phi} = name of the finite volume flux field :l
 {switch1} = (optional, default false) indicator to split calculation of Urel between CFDEM and LIGGGHTS :l
+{switch2} = (optional, default false) indicator to use the cell-averaged particle velocity Us for calculation of Urel :l
+{switch3} = (optional, default false) indicator to use the correlation by Felderhof for the virtual mass coefficient :l
 {scalar1} = (optional, default 0.5) virtual mass coefficient :l
-{switch2} = (optional, default true) sub model switch, see "forceSubModel"_forceSubModel.html for details :l
-{switch3} = (optional, default false) sub model switch, see "forceSubModel"_forceSubModel.html for details :l
+{switch4} = (optional, default true) sub model switch, see "forceSubModel"_forceSubModel.html for details :l
+{switch5} = (optional, default false) sub model switch, see "forceSubModel"_forceSubModel.html for details :l
+{smoothingModel} = name of smoothing model for the DU/Dt field :l
 :ule
 
 [Examples:]
@@ -43,7 +53,10 @@ forceModels
 virtualMassForceProps
 \{
     velFieldName "U";
+    voidfractionFieldName "voidfraction"
+    USFieldName "Us"
     phiFieldName "phi";
+    smoothingModel off;
 \} :pre
 
 [Description:]
@@ -52,6 +65,15 @@ The force model performs the calculation of forces (e.g. fluid-particle
 interaction forces) acting on each DEM particle. The {virtualMassForce} model
 calculates the virtual mass force for each particle.
 
+For more detail, see "Nijssen et al. (2020)"_#Nijssen2020.
+
+:link(Nijssen2020)
+[(Nijssen2020)] T.M.J. Nijssen, J.A.M. Kuipers, J. van der Stel, A.T. Adema, K.A. Buist. (2020):
+"Complete liquid-solid momentum coupling for unresolved CFD-DEM simulations",
+International Journal of Multiphase Flow
+
+[Restrictions:]
+
 IMPORTANT NOTE: Model not validated!
 
 [Related commands:]

doc/smoothingModel_constDiffAndTemporalSmoothing.txt

@@ -0,0 +1,77 @@
+"CFDEMproject Website"_lws - "Main Page"_main :c
+
+:link(lws,http://www.cfdem.com)
+:link(main,CFDEMcoupling_Manual.html)
+
+:line
+
+smoothingModel constDiffAndTemporalSmoothing command :h3
+
+[Syntax:]
+
+Defined in dictionary depending on the application.
+
+smoothingModel constDiffAndTemporalSmoothing;
+constDiffAndTemporalSmoothingProps
+\{
+    lowerLimit                    number1;
+    upperLimit                    number2;
+    smoothingLength               lengthScale;
+    smoothingLengthReferenceField lengthScaleRefField;
+    smoothingStrength             smoothingStrength;
+    correctBoundary               switch1;
+    verbose;
+\} :pre
+
+{number1} = scalar fields will be bound to this lower value :ulb,l
+{number2} = scalar fields will be bound to this upper value :l
+{lengthScale} = length scale over which the exchange fields will be smoothed out :l
+{lengthScaleRefField} = length scale over which reference fields (e.g., the average particle velocity) will be smoothed out. Should be always larger than lengthScale. If not specified, will be equal to lengthScale. :l
+{smoothingStrength} = control parameter gamma for the smoothing, lower value yields stronger smoothing (gamma = 1 results in an equal contribution from the un-smoothed and smoothed fields)  :l
+{correctBoundary} = (optional, default false) activate to use purely temporal smoothing on the boundary field, avoids interpolation errors near the domain boundary :l
+{verbose} = (optional, default false) flag for debugging output :l
+:ule
+
+[Examples:]
+
+constDiffAndTemporalSmoothingProps
+\{
+    lowerLimit       0.1;
+    upperLimit       1e10;
+    smoothingLength  1500e-6;
+    smoothingLengthReferenceField 9000e-6;
+    referenceField   "p";
+    gamma            1.0;
+\} :pre
+
+[Description:]
+
+The {constDiffAndTemporalSmoothing} model is a smoothing model that combines the
+spacial smoothing of "constDiffSmoothing"_smoothingModel_constDiffSmoothing.html and
+the temporal smoothing of "temporalSmoothing"_smoothingModel_temporalSmoothing.html for
+the relaxation of a desired quantity. This model can be used to filter out
+high-frequency fluctuations (e.g. numerical noise) controlled via the smoothing length
+and the temporal smoothing strength parameter gamma.
+
+For more details, see "Vångö et al. (2018)"_#Vångö2018 and "Nijssen et al. (2020)"_#Nijssen2020.
+
+:link(Vångö2018)
+[(Vångö2018)] M. Vångö, S. Pirker, T. Lichtenegger. (2018):
+"Unresolved CFD-DEM modeling of multiphase flow in densely packed particle beds",
+Applied Mathematical Modelling
+
+:link(Nijssen2020)
+[(Nijssen2020)] T.M.J. Nijssen, J.A.M. Kuipers, J. van der Stel, A.T. Adema, K.A. Buist. (2020):
+"Complete liquid-solid momentum coupling for unresolved CFD-DEM simulations",
+International Journal of Multiphase Flow
+
+[Restrictions:]
+
+This model is tested in a limited number of flow situations.
+
+[Related commands:]
+
+"smoothingModel"_smoothingModel.html
+"constDiffSmoothing"_smoothingModel_constDiffSmoothing.html
+"temporalSmoothing"_smoothingModel_temporalSmoothing.html
+

etc/bashrc

@@ -17,7 +17,7 @@
 #------------------------------------------------------------------------------
 
 export CFDEM_PROJECT=CFDEM
-export CFDEM_VERSION=21.03
+export CFDEM_VERSION=21.11
 
 ################################################################################
 # USER EDITABLE PART: Changes made here may be lost with the next upgrade
@@ -221,7 +221,7 @@ alias cfdemCompCFDEMuti='bash $CFDEM_PROJECT_DIR/etc/compileCFDEMcoupling_uti.sh
 alias cfdemTestTUT='bash $CFDEM_PROJECT_DIR/etc/testTutorials.sh'
 
 #- shortcut to visualize the clock model data
-alias vizClock='python $CFDEM_UT_DIR/vizClock/matPlot.py'
+alias vizClock='python2 $CFDEM_UT_DIR/vizClock/matPlot.py'
 
 #- recursive touch of current directory
 alias touchRec='find ./* -exec touch {} \;'