Merge pull request #146 from ParticulateFlow/release

Release 23.02

.github/pull_request_template.md

@@ -0,0 +1,24 @@
+<!-- Please provide a general summary of your changes in the title above. -->
+
+## Description of proposed changes
+<!-- Describe your changes in detail. -->
+
+## Types of changes
+<!-- What types of changes does your code introduce? Put an `x` in all the boxes that apply. -->
+<!-- Please try to limit your pull request to one type, submit multiple pull requests if needed. -->
+- [ ] Bugfix
+- [ ] Feature
+- [ ] Refactoring (no functional changes, no api changes)
+- [ ] Build related changes
+- [ ] Documentation updates
+- [ ] Other (please describe): 
+
+## Checklist
+<!-- Go over all the following points, and put an `x` in all the boxes that apply. -->
+- [ ] Code compiles correctly (mandatory for bugfixes / features / refactoring / build process)
+- [ ] Tests for the changes have been added / updated (mandatory for bugfixes / features)
+- [ ] Documentation has been added / updated (mandatory for bugfixes / features)
+
+## Further comments
+<!-- If this is a relatively large or complex change, kick off the discussion by explaining
+why you chose the solution you did and what alternatives you considered, etc... -->

.gitignore

@@ -8,7 +8,7 @@ log.*
 *.swp
 *.swo
 
-**/linux*Gcc*/
+**/linux*cc*/
 **/.vscode
 
 lnInclude

applications/solvers/cfdemSolverIB/Make/options

@@ -9,12 +9,12 @@ EXE_INC = \
     -I$(LIB_SRC)/transportModels \
     -I$(LIB_SRC)/transportModels/incompressible/singlePhaseTransportModel \
     -I$(CFDEM_SRC_DIR)/lagrangian/cfdemParticle/lnInclude \
-    -I$(LIB_SRC)/sampling/lnInclude \
     -I$(LIB_SRC)/dynamicFvMesh/lnInclude \
     -I$(LIB_SRC)/dynamicMesh/lnInclude \
     -I$(LIB_SRC)/dynamicMesh/dynamicFvMesh/lnInclude \
     -I$(LIB_SRC)/dynamicMesh/dynamicMesh/lnInclude \
     -I$(LIB_SRC)/fvOptions/lnInclude \
+    -I$(LIB_SRC)/sampling/lnInclude \
     -Wno-deprecated-copy
 
 EXE_LIBS = \
@@ -27,6 +27,7 @@ EXE_LIBS = \
     -ldynamicFvMesh \
     -ldynamicMesh \
     -lfvOptions \
+    -lsampling \
     -l$(CFDEM_LIB_NAME) \
      $(CFDEM_ADD_LIB_PATHS) \
      $(CFDEM_ADD_LIBS)

applications/solvers/cfdemSolverIB/cfdemSolverIB.C

@@ -6,7 +6,8 @@
                                 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
+                                Copyright (C) 2012-2015 DCS Computing GmbH,Linz
+                                Copyright (C) 2015-     JKU, Linz
 -------------------------------------------------------------------------------
 License
     This file is part of CFDEMcoupling.
@@ -29,11 +30,14 @@ Application
 
 Description
     Transient solver for incompressible flow.
-    The code is an evolution of the solver pisoFoam in OpenFOAM(R) 1.6, 
+    The code is an evolution of the solver pisoFoam in OpenFOAM(R) 1.6,
     where additional functionality for CFD-DEM coupling using immersed body
     (fictitious domain) method is added.
 Contributions
     Alice Hager
+    Daniel Queteschiner
+    Thomas Lichtenegger
+    Achuth N. Balachandran Nair
 \*---------------------------------------------------------------------------*/
 
 
@@ -53,23 +57,21 @@ Contributions
 
 #include "cellSet.H"
 
+#include "fvOptions.H"  // added the fvOptions library
+
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 int main(int argc, char *argv[])
 {
+
     #include "setRootCase.H"
-
     #include "createTime.H"
-
     #include "createDynamicFvMesh.H"
-
     #include "createControl.H"
-
     #include "createTimeControls.H"
-
     #include "createFields.H"
-
     #include "initContinuityErrs.H"
+    #include "createFvOptions.H"
 
     // create cfdemCloud
     #include "readGravitationalAcceleration.H"
@@ -93,24 +95,31 @@ int main(int argc, char *argv[])
 
         // do particle stuff
         Info << "- evolve()" << endl;
-        particleCloud.evolve();
+        particleCloud.evolve(Us);
 
         // Pressure-velocity PISO corrector
         {
+            MRF.correctBoundaryVelocity(U);
+
             // Momentum predictor
 
             fvVectorMatrix UEqn
             (
-                fvm::ddt(voidfraction,U)
+                fvm::ddt(voidfraction,U) + MRF.DDt(U)
               + fvm::div(phi, U)
               + turbulence->divDevReff(U)
+             ==
+                fvOptions(U)
             );
 
             UEqn.relax();
 
+            fvOptions.constrain(UEqn);
+
             if (piso.momentumPredictor())
             {
                 solve(UEqn == -fvc::grad(p));
+                fvOptions.correct(U);
             }
 
             // --- PISO loop
@@ -126,6 +135,7 @@ int main(int argc, char *argv[])
 
                 adjustPhi(phi, U, p);
 
+
                 while (piso.correctNonOrthogonal())
                 {
                     // Pressure corrector
@@ -152,12 +162,15 @@ int main(int argc, char *argv[])
             }
         }
 
+        laminarTransport.correct();
         turbulence->correct();
 
         Info << "particleCloud.calcVelocityCorrection() " << endl;
         volScalarField voidfractionNext=mesh.lookupObject<volScalarField>("voidfractionNext");
         particleCloud.calcVelocityCorrection(p,U,phiIB,voidfractionNext);
 
+        fvOptions.correct(U);
+
         runTime.write();
 
         Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"

applications/solvers/cfdemSolverIB/createFields.H

@@ -26,21 +26,6 @@
         ),
         mesh
     );
-    //mod by alice
-    Info<< "Reading physical velocity field U" << endl;
-    Info<< "Note: only if voidfraction at boundary is 1, U is superficial velocity!!!\n" << endl;
-    volVectorField Us
-    (
-        IOobject
-        (
-            "Us",
-            runTime.timeName(),
-            mesh,
-            IOobject::MUST_READ,
-            IOobject::AUTO_WRITE
-        ),
-        mesh
-    );
 
 //========================
 // drag law modelling
@@ -76,9 +61,8 @@
         mesh
     );
 
-
     //mod by alice
-    Info<< "Reading field phiIB\n" << endl;
+    Info<< "Reading field voidfraction\n" << endl;
     volScalarField voidfraction
     (
         IOobject
@@ -91,6 +75,21 @@
         ),
         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"
@@ -126,3 +125,5 @@
     );
 
 //===========================
+
+#include "createMRF.H"

applications/solvers/cfdemSolverIBContinuousForcing/Make/files

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

applications/solvers/cfdemSolverIBContinuousForcing/Make/options

@@ -0,0 +1,34 @@
+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$(LIB_SRC)/dynamicFvMesh/lnInclude \
+    -I$(LIB_SRC)/dynamicMesh/lnInclude \
+    -I$(LIB_SRC)/dynamicMesh/dynamicFvMesh/lnInclude \
+    -I$(LIB_SRC)/dynamicMesh/dynamicMesh/lnInclude \
+    -I$(LIB_SRC)/fvOptions/lnInclude \
+    -I$(LIB_SRC)/sampling/lnInclude \
+    -Wno-deprecated-copy
+
+EXE_LIBS = \
+    -L$(CFDEM_LIB_DIR)\
+    -lturbulenceModels \
+    -lincompressibleTurbulenceModels \
+    -lincompressibleTransportModels \
+    -lfiniteVolume \
+    -lmeshTools \
+    -ldynamicFvMesh \
+    -ldynamicMesh \
+    -lfvOptions \
+    -lsampling \
+    -l$(CFDEM_LIB_NAME) \
+     $(CFDEM_ADD_LIB_PATHS) \
+     $(CFDEM_ADD_LIBS)
+

applications/solvers/cfdemSolverIBContinuousForcing/UEqn.H

@@ -0,0 +1,19 @@
+fvVectorMatrix UEqn
+(
+    fvm::ddt(voidfractionNext,U) + MRF.DDt(U)
+  + fvm::div(phi, U)
+  + turbulence->divDevReff(U)
+ ==
+    fvOptions(U)
+ +  (lambda*(1-voidfractionNext)/U.mesh().time().deltaT())*(fvc::Sp(1,Us)-fvm::Sp(1,U))
+);
+
+UEqn.relax();
+
+fvOptions.constrain(UEqn);
+
+if (piso.momentumPredictor())
+{
+    solve(UEqn == -fvc::grad(p));
+    fvOptions.correct(U);
+}

applications/solvers/cfdemSolverIBContinuousForcing/cfdemSolverIBContinuousForcing.C

@@ -0,0 +1,129 @@
+/*---------------------------------------------------------------------------*\
+    CFDEMcoupling - Open Source CFD-DEM coupling
+
+    CFDEMcoupling is part of the CFDEMproject
+    www.cfdem.com
+                                Copyright (C) 1991-2009 OpenCFD Ltd.
+                                Copyright (C) 2009-2012 JKU, Linz
+                                Copyright (C) 2012-2015 DCS Computing GmbH,Linz
+                                Copyright (C) 2015-     JKU, 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
+    cfdemSolverIBContinuousForcing
+
+Description
+    Transient solver for incompressible flow.
+    The code is an evolution of the solver pisoFoam in OpenFOAM(R) 1.6,
+    where additional functionality for CFD-DEM coupling using immersed body
+    (fictitious domain) method and a continuous forcing approach is added.
+Contributions
+    Alice Hager
+    Achuth N. Balachandran Nair
+\*---------------------------------------------------------------------------*/
+
+
+#include "fvCFD.H"
+#include "singlePhaseTransportModel.H"
+#include "turbulentTransportModel.H"
+#include "pisoControl.H"
+
+#include "cfdemCloudIB.H"
+#include "implicitCouple.H"
+
+#include "averagingModel.H"
+#include "regionModel.H"
+#include "voidFractionModel.H"
+
+#include "dynamicFvMesh.H"
+
+#include "cellSet.H"
+
+#include "fvOptions.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+int main(int argc, char *argv[])
+{
+
+    #include "setRootCase.H"
+    #include "createTime.H"
+    #include "createDynamicFvMesh.H"
+    #include "createControl.H"
+    #include "createTimeControls.H"
+    #include "createFields.H"
+    #include "initContinuityErrs.H"
+    #include "createFvOptions.H"
+
+    // create cfdemCloud
+    #include "readGravitationalAcceleration.H"
+    cfdemCloudIB particleCloud(mesh);
+
+    // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+    Info<< "\nStarting time loop\n" << endl;
+
+    while (runTime.loop())
+    {
+        Info<< "Time = " << runTime.timeName() << nl << endl;
+
+        //=== dyM ===================
+        interFace = mag(mesh.lookupObject<volScalarField>("voidfractionNext"));
+        mesh.update(); //dyM
+
+        #include "readTimeControls.H"
+        #include "CourantNo.H"
+        #include "setDeltaT.H"
+
+        // do particle stuff
+        Info << "- evolve()" << endl;
+        particleCloud.evolve(Us);
+
+        volScalarField voidfractionNext=mesh.lookupObject<volScalarField>("voidfractionNext");
+
+        // Pressure-velocity PISO corrector
+        {
+            MRF.correctBoundaryVelocity(U);
+
+            // Momentum predictor
+            #include "UEqn.H"
+
+            // --- PISO loop
+            while (piso.correct())
+            {
+                #include "pEqn.H"
+            }
+        }
+
+        laminarTransport.correct();
+        turbulence->correct();
+
+        runTime.write();
+
+        Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
+            << "  ClockTime = " << runTime.elapsedClockTime() << " s"
+            << nl << endl;
+    }
+
+    Info<< "End\n" << endl;
+
+    return 0;
+}
+
+
+// ************************************************************************* //

applications/solvers/cfdemSolverIBContinuousForcing/createFields.H

@@ -0,0 +1,143 @@
+    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
+    );
+
+    Info<< "Reading particle velocity field Us\n" << endl;
+    volVectorField Us
+    (
+        IOobject
+        (
+            "Us",
+            runTime.timeName(),
+            mesh,
+            IOobject::MUST_READ,
+            IOobject::AUTO_WRITE
+        ),
+        mesh
+    );
+
+    Info<< "Reading the penalization factor field lambda\n" << endl;
+    volScalarField lambda
+    (
+        IOobject
+        (
+            "lambda",
+            runTime.timeName(),
+            mesh,
+            IOobject::MUST_READ,
+            IOobject::AUTO_WRITE
+        ),
+        mesh
+    );
+
+//========================
+// drag law modelling
+//========================
+
+    Info<< "\nCreating dummy density field rho = 1\n" << endl;
+    volScalarField rho
+    (
+        IOobject
+        (
+            "rho",
+            runTime.timeName(),
+            mesh,
+            IOobject::READ_IF_PRESENT,
+            IOobject::AUTO_WRITE
+        ),
+        mesh,
+        dimensionedScalar("0", dimensionSet(1, -3, 0, 0, 0), 1.0)
+    );
+
+
+    Info<< "Reading field phiIB\n" << endl;
+    volScalarField phiIB
+    (
+        IOobject
+        (
+            "phiIB",
+            runTime.timeName(),
+            mesh,
+            IOobject::MUST_READ,
+            IOobject::AUTO_WRITE
+        ),
+        mesh
+    );
+
+    //mod by alice
+    Info<< "Reading field voidfraction\n" << endl;
+    volScalarField voidfraction
+    (
+        IOobject
+        (
+            "voidfraction",
+            runTime.timeName(),
+            mesh,
+            IOobject::MUST_READ,
+            IOobject::AUTO_WRITE
+        ),
+        mesh
+    );
+
+//========================
+
+#   include "createPhi.H"
+
+    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)
+    );
+
+//=== dyM ===================
+
+    Info<< "Reading field interFace\n" << endl;
+    volScalarField interFace
+    (
+        IOobject
+        (
+            "interFace",
+            runTime.timeName(),
+            mesh,
+            IOobject::READ_IF_PRESENT,
+            IOobject::AUTO_WRITE
+        ),
+        mesh,
+        //dimensionedScalar("0", dimensionSet(0, -1, 0, 0, 0), 0.0)
+        dimensionedScalar("0", dimensionSet(0, 0, 0, 0, 0), 0.0)
+    );
+
+//===========================
+
+#include "createMRF.H"

applications/solvers/cfdemSolverIBContinuousForcing/pEqn.H

@@ -0,0 +1,35 @@
+volScalarField rUA = 1.0/UEqn.A();
+surfaceScalarField rUAf(fvc::interpolate(rUA));
+
+U = rUA*UEqn.H();
+
+phi = (fvc::interpolate(U) & mesh.Sf())
+    + rUAf*fvc::ddtCorr(U, phi);            // Is there additional flux term due to the particle presence?
+
+MRF.makeRelative(phi);
+
+adjustPhi(phi, U, p);
+
+while (piso.correctNonOrthogonal())
+{
+    // Pressure corrector
+
+    fvScalarMatrix pEqn
+    (
+        fvm::laplacian(rUA, p) == fvc::div(phi) + particleCloud.ddtVoidfraction()
+    );
+
+    pEqn.setReference(pRefCell, pRefValue);
+
+    pEqn.solve(mesh.solver(p.select(piso.finalInnerIter())));
+
+    if (piso.finalNonOrthogonalIter())
+    {
+        phi -= pEqn.flux();
+    }
+}
+
+#include "continuityErrs.H"
+U -= rUA*fvc::grad(p); // should we add a pressure correction?
+U.correctBoundaryConditions();
+fvOptions.correct(U);

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 multiphaseMixtureScalar
+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 multiphaseMixtureScalar
+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) \
+    -ImultiphaseMixtureScalar/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 "multiphaseMixtureScalar.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()
+ );
+
+multiphaseMixtureScalar 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/multiphaseMixtureScalar/Make/files

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

applications/solvers/cfdemSolverMultiphaseScalar/multiphaseMixtureScalar/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/multiphaseMixtureScalar/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/multiphaseMixtureScalar/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 multiphaseMixtureScalar.
+
+SourceFiles
+    alphaContactAngleFvPatchScalarField.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef alphaContactAngleFvPatchScalarField_H
+#define alphaContactAngleFvPatchScalarField_H
+
+#include "zeroGradientFvPatchFields.H"
+#include "multiphaseMixtureScalar.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,
+        multiphaseMixtureScalar::interfacePair,
+        multiphaseMixtureScalar::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/multiphaseMixtureScalar/multiphaseMixtureScalar.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 "multiphaseMixtureScalar.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::multiphaseMixtureScalar::convertToRad =
+    Foam::constant::mathematical::pi/180.0;
+
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+void Foam::multiphaseMixtureScalar::calcAlphas()
+{
+    scalar level = 0.0;
+    alphas_ == 0.0;
+
+    forAllIter(PtrDictionary<phase>, phases_, iter)
+    {
+        alphas_ += level*iter();
+        level += 1.0;
+    }
+}
+
+
+Foam::tmp<Foam::volScalarField>
+Foam::multiphaseMixtureScalar::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::multiphaseMixtureScalar::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::multiphaseMixtureScalar::multiphaseMixtureScalar
+(
+    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::multiphaseMixtureScalar::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::multiphaseMixtureScalar::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::multiphaseMixtureScalar::mu() const
+{
+    Info << "In multiphaseMixtureScalar 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::multiphaseMixtureScalar::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::multiphaseMixtureScalar::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::multiphaseMixtureScalar::nu() const
+{
+    return nu_;
+}
+
+Foam::tmp<Foam::scalarField>
+Foam::multiphaseMixtureScalar::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::multiphaseMixtureScalar::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::multiphaseMixtureScalar::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::multiphaseMixtureScalar::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::multiphaseMixtureScalar::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::multiphaseMixtureScalar::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::multiphaseMixtureScalar::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::multiphaseMixtureScalar::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::multiphaseMixtureScalar::correct()
+{
+    forAllIter(PtrDictionary<phase>, phases_, iter)
+    {
+        iter().correct();
+    }
+}
+
+
+Foam::tmp<Foam::surfaceVectorField> Foam::multiphaseMixtureScalar::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::multiphaseMixtureScalar::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::multiphaseMixtureScalar::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::multiphaseMixtureScalar::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::multiphaseMixtureScalar::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::multiphaseMixtureScalar::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::multiphaseMixtureScalar::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/multiphaseMixtureScalar/multiphaseMixtureScalar.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
+    multiphaseMixtureScalar
+
+Description
+    This class is based on the OpenFOAM(R) Foam::multiphaseMixtureScalar 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
+    multiphaseMixtureScalar.C
+\*---------------------------------------------------------------------------*/
+
+#ifndef multiphaseMixtureScalar_H
+#define multiphaseMixtureScalar_H
+
+#include "incompressible/transportModel/transportModel.H"
+#include "IOdictionary.H"
+#include "phase.H"
+#include "PtrDictionary.H"
+#include "volFields.H"
+#include "surfaceFields.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                      Class multiphaseMixtureScalar Declaration
+\*---------------------------------------------------------------------------*/
+
+class multiphaseMixtureScalar
+:
+    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
+        multiphaseMixtureScalar
+        (
+            const volVectorField& U,
+            const surfaceScalarField& phi,
+            const volScalarField& voidfraction
+        );
+
+
+    //- Destructor
+    virtual ~multiphaseMixtureScalar()
+    {}
+
+
+    // 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/multiphaseMixtureScalar/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/multiphaseMixtureScalar/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 multiphaseMixtureScalar 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,29 @@
+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 \
+    -Wno-deprecated-copy
+
+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,9 +50,5 @@
 
     thermo.correct();
 
-    Info<< "T max/min : " << max(T).value() << " " << min(T).value() << endl;
-
-    particleCloud.clockM().start(31,"energySolve");
-    particleCloud.solve();
-    particleCloud.clockM().stop("energySolve");
+    Info<< "T max/min/ave : " << max(T).value() << " " << min(T).value() << " " << average(T).value() << endl;
 }

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)
@@ -53,8 +55,4 @@ fvScalarMatrix EEqn
     Info << "Cpv :" << max(Cpv).value() << " " << min(Cpv).value() << endl;
     Info<<  "T max/min : " << max(T).value() << " " << min(T).value() << endl;
     Info << "he max/min : " << max(he).value() << "  " << min(he).value() << endl;
-
-    particleCloud.clockM().start(31,"energySolve");
-    particleCloud.solve();
-    particleCloud.clockM().stop("energySolve");
 }

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/YEqn.H

@@ -60,13 +60,13 @@ tmp<fv::convectionScheme<scalar> > mvConvection
 
     if (propagateInertSpecie)
     {
-        if (inertIndex!=-1) Yt /= (1-Y[inertIndex] + VSMALL);
+        if (inertIndex!=-1) Yt /= (1-Y[inertIndex] + ROOTVSMALL);
         forAll(Y,i)
         {
             if (i!=inertIndex)
             {
                 volScalarField& Yi = Y[i];
-                Yi = Yi/(Yt+VSMALL);
+                Yi = Yi/(Yt+ROOTVSMALL);
             }
         }
     }

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;
 
@@ -303,18 +303,4 @@
         mesh,
         dimensionedScalar("zero",dimensionSet(0, -3, 0, 0, 1),0)
     );
-
-    volScalarField dSauter
-    (
-        IOobject
-        (
-            "dSauter",
-            runTime.timeName(),
-            mesh,
-            IOobject::READ_IF_PRESENT,
-            IOobject::AUTO_WRITE
-        ),
-        mesh,
-        dimensionedScalar("zero",dimensionSet(0, 1, 0, 0, 0,0,0),0)
-    );
 //===============================

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/TEqn.H

@@ -1,4 +1,4 @@
-    volScalarField rhoeps = rhoRec*voidfractionRec;
+    rhoeps = rhoRec*voidfractionRec;
 
     particleCloud.energyContributions(Qsource);
 
@@ -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)
@@ -22,10 +23,17 @@
         fvOptions(rhoeps, T)    // no fvOptions support yet
     );
 
-    fvOptions.constrain(TEqn); // no fvOptions support yet
+ //   fvOptions.constrain(TEqn); // no fvOptions support yet
+
+    TEqn.relax();
 
     TEqn.solve();
 
+    T = max(T, TMin);
+    T = min(T, TMax);
+
+    Info<< "T max/min/ave : " << max(T).value() << " " << min(T).value() << " " << average(T).value() << endl;
+
     particleCloud.clockM().start(31,"postFlow");
     counter++;
 

applications/solvers/rcfdemSolverCoupledHeattransfer/createFields.H

@@ -73,6 +73,9 @@
         dimensionedVector("zero", dimensionSet(0, 1, -1, 0, 0), vector::zero)
     );
 
+    volScalarField rhoeps("rhoeps", rhoRec*voidfractionRec);
+    rhoeps.oldTime(); // switch on saving old time
+
     // heat transfer fields
     Info << "\nCreating heat transfer fields.\n" << endl;
 
@@ -228,3 +231,25 @@
         )
     );
     weightDict.add("weights",scalarList(1,1.0));
+
+    dimensionedScalar TMax
+    (
+        dimensionedScalar::lookupOrDefault
+        (
+            "TMax",
+            transportProps,
+            dimTemperature,
+            GREAT
+        )
+    );
+
+    dimensionedScalar TMin
+    (
+        dimensionedScalar::lookupOrDefault
+        (
+            "TMin",
+            transportProps,
+            dimTemperature,
+            0.0
+        )
+    );

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
@@ -98,15 +98,19 @@ int main(int argc, char *argv[])
 
         particleCloud.clockM().start(26,"Flow");
         #include "updateRho.H"
-        #include "TEqImp.H"
+        #include "TEqn.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/rcfdemSolverCoupledHeattransfer/updateRho.H

@@ -1 +1,2 @@
-rhoRec = pRec / (T * R);
+dimensionedScalar Tave = T.weightedAverage(voidfractionRec);
+rhoRec = pRec / (Tave * R);

applications/solvers/rcfdemSolverForcedTracers/Make/files

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

applications/solvers/rcfdemSolverForcedTracers/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$(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 \
+    -Wno-deprecated-copy
+
+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);
@@ -46,9 +60,4 @@
     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/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,62 @@
+// 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;
+}

applications/solvers/rcfdemSolverRhoSteadyPimpleChem/Make/files

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

applications/solvers/rcfdemSolverRhoSteadyPimpleChem/Make/options

@@ -0,0 +1,53 @@
+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 \
+    -Wno-deprecated-copy
+
+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,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$(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 \
+    -Wno-deprecated-copy
+
+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,562 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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, scalarList &, volVectorField &, volVectorField &, scalarList &, volVectorField &, volVectorField &, bool, scalar);
+void nearestNeighborCells(fvMesh &, label, label, label, scalarList &, labelList &);
+void normalizeFields(scalarList &, volVectorField &, volVectorField &);
+void readDump(std::string, labelList &, scalarList &, vectorList &);
+scalar weightFun(scalar);
+label maxNumParticles = 1000000;
+scalar minVol = 1e-12;
+scalar Pi43 = 4.1888;
+label posIndex = -1;
+label posRad = -1;
+label posX = -1;
+label posY = -1;
+label posZ = -1;
+
+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));
+    posIndex = readLabel(displacementProperties.lookup("posIndex"));
+    posRad = readLabel(displacementProperties.lookup("posRad"));
+    posX = readLabel(displacementProperties.lookup("posX"));
+    posY = readLabel(displacementProperties.lookup("posY"));
+    posZ = readLabel(displacementProperties.lookup("posZ"));
+    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 defaultUsDirectedStdDev
+    (
+        IOobject
+        (
+            "defaultUDispDirectedStdDev",
+            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 UsDirectedStdDev
+    (
+        IOobject
+        (
+            "UDispDirectedStdDev",
+            runTime.timeName(),
+            mesh,
+            IOobject::NO_READ,
+            IOobject::NO_WRITE
+        ),
+        mesh,
+        dimensionedVector("zero", dimensionSet(0,1,-1,0,0), vector::zero)
+    );
+
+    scalarList particleVolInCell(mesh.nCells(), 0.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;
+        scalarList radii1, radii2;
+        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(particleVolInCell, Us, UsDirectedStdDev);
+                fillEmptyCells(mesh,nNeighMin,maxSearchLayers,particleVolInCell,Us,UsDirectedStdDev,boundaries,defaultUs,defaultUsDirectedStdDev,interpolate,timePerDisplacementStep);
+
+                Us /= timePerDisplacementStep;
+                UsDirectedStdDev /= timePerDisplacementStep;
+                Us.write();
+                UsDirectedStdDev.write();
+            }
+            break;
+        }
+
+        Info << "\nReading" << endl;
+        Info << "\t" << filename1 << endl;
+        Info << "\t" << filename2 << endl;
+        Info << "corresponding to time = " << currTime << "." << endl;
+
+        readDump(filename1, indices1, radii1, positions1);
+        readDump(filename2, indices2, radii2, positions2);
+
+        labelPairList pairs;
+        findPairs(indices1,indices2,pairs);
+
+        // average particle displacements and their variance
+        Info << "Binning particle displacements on mesh." << endl;
+        vector position, displacement;
+        scalar radius, volume;
+        label line1, line2;
+        label cellI;
+
+        if (!averageMode)
+        {
+            Us *= 0.0;
+            UsDirectedStdDev *= 0.0;
+            particleVolInCell.clear();
+            particleVolInCell.setSize(mesh.nCells(), 0);
+        }
+
+        for (label partI = 0; partI < pairs.size(); partI++)
+        {
+            line1 = pairs[partI].first();
+            line2 = pairs[partI].second();
+            position = positions1[line1];
+            cellI = mesh.findCell(position);
+            if (cellI < 0)  continue;
+            displacement = positions2[line2] - positions1[line1];
+            radius = radii1[line1];
+            volume = Pi43 * radius * radius * radius;
+            particleVolInCell[cellI] += volume;
+            Us[cellI] += displacement*volume;
+
+            for (label comp=0;comp<3;comp++)
+            {
+                UsDirectedStdDev[cellI].component(comp) += displacement.component(comp)*displacement.component(comp)*volume;
+            }
+        }
+
+        if (!averageMode)
+        {
+            normalizeFields(particleVolInCell, Us, UsDirectedStdDev);
+            fillEmptyCells(mesh,nNeighMin,maxSearchLayers,particleVolInCell,Us,UsDirectedStdDev,boundaries,defaultUs,defaultUsDirectedStdDev,interpolate,timePerDisplacementStep);
+
+            Us /= timePerDisplacementStep;
+            UsDirectedStdDev /= timePerDisplacementStep;
+            Us.write();
+            UsDirectedStdDev.write();
+        }
+
+        if (averageMode && monitorProbes)
+        {
+            monitoringDataFile << currTime << " ";
+            for(label p=0;p<probePoints.size();p++)
+            {
+                vector pos = probePoints[p];
+                label cellP = mesh.findCell(pos);
+                monitoringDataFile << " " << particleVolInCell[cellP] << " " << Us[cellP]/timePerDisplacementStep << " " << UsDirectedStdDev[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, scalarList &radii, vectorList &positions)
+{
+    #include <fstream>
+
+    const label leadingLines = 9;
+    label lineCounter = 0;
+    label partIndex;
+    scalar r = 1.0, x = 0.0, y = 0.0, z = 0.0;
+
+    indices.clear();
+    radii.clear();
+    positions.clear();
+
+    indices.setSize(maxNumParticles);
+    radii.setSize(maxNumParticles);
+    positions.setSize(maxNumParticles);
+
+    std::ifstream file(filename);
+    std::string str;
+    std::string word;
+    label wordcounter;
+    while (std::getline(file, str))
+    {
+        if (lineCounter >= leadingLines)
+        {
+            std::istringstream ss(str);
+            wordcounter = 0;
+            while (ss >> word)
+            {
+                if (wordcounter == posIndex)
+                {
+                    partIndex = stoi(word);
+                }
+                else if (wordcounter == posRad)
+                {
+                    r = stod(word);
+                }
+                else if (wordcounter == posX)
+                {
+                    x = stod(word);
+                }
+                else if (wordcounter == posY)
+                {
+                    y = stod(word);
+                }
+                else if (wordcounter == posZ)
+                {
+                    z = stod(word);
+                }
+                wordcounter++;
+            }
+//            sscanf(str.c_str(), "%d %lf %lf %lf", &partIndex, &x, &y, &z);
+            indices[lineCounter-leadingLines] = partIndex;
+            radii[lineCounter-leadingLines] = r;
+            positions[lineCounter-leadingLines] = vector(x,y,z);
+        }
+        lineCounter++;
+    }
+
+    label readLines = lineCounter - leadingLines;
+    indices.resize(readLines);
+    radii.resize(readLines);
+    positions.resize(readLines);
+}
+
+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
+
+    pairs.clear();
+    pairs.setSize(maxNumParticles);
+    label pairCounter = 0;
+
+    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[pairCounter]=labelPair(i,jmid);
+                pairCounter++;
+                break;
+            }
+            if (j2-j1 == 1) break;
+        }
+    }
+    pairs.resize(pairCounter);
+    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
+    pairs.clear();
+    pairs.setSize(maxNumParticles);
+    label pairCounter = 0;
+
+    for (label i=0;i<indices1.size();i++)
+    {
+        for (label j=0;j<indices2.size();j++)
+        {
+            if (indices1[i] == indices2[j])
+            {
+                pairs[pairCounter]=labelPair(i,j);
+                pairCounter++;
+                break;
+            }
+        }
+    }
+    pairs.resize(pairCounter);
+    Info << "findPairs: " << pairs.size() << " pairs found." << endl;
+}
+
+void fillEmptyCells(fvMesh &mesh, label nNeighMin, label maxSearchLayers, scalarList &particleVolInCell, volVectorField &Us, volVectorField& UsDirectedStdDev,scalarList& boundaries, volVectorField &defaultUs, volVectorField &defaultUsDirectedStdDev, bool interpolate, scalar dt)
+{
+    labelList neighborsWithValues;
+    scalar neighborSqrDistance;
+    scalar weight;
+    scalar weightSum;
+    scalarList weights;
+
+    Info << "Filling empty cells." << endl;
+    forAll(mesh.C(), cellI)
+    {
+        if (particleVolInCell[cellI] > minVol) 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;
+            UsDirectedStdDev[cellI] = defaultUsDirectedStdDev[cellI]*dt;
+            continue;
+        }
+
+        nearestNeighborCells(mesh, cellI, nNeighMin, maxSearchLayers, particleVolInCell, 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]];
+            UsDirectedStdDev[cellI] += weight*UsDirectedStdDev[neighborsWithValues[neighI]];
+        }
+
+        if (neighborsWithValues.size() == 0)
+        {
+            Us[cellI] = defaultUs[cellI]*dt;
+            UsDirectedStdDev[cellI] = defaultUsDirectedStdDev[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, scalarList &particleVolInCell, 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 (particleVolInCell[adj] > minVol) neighborsWithValues.append(adj);
+                }
+            }
+        }
+
+        numSearchLayers++;
+        if (numSearchLayers > maxSearchLayers && maxSearchLayers > 0) return;
+
+        if (newNeighbors.size() == 0) return;
+        recentNeighbors.clear();
+        recentNeighbors = newNeighbors;
+        newNeighbors.clear();
+    }
+}
+
+void normalizeFields(scalarList& particleVolInCell, volVectorField& Us, volVectorField & UsDirectedStdDev)
+{
+    for (label cellJ = 0; cellJ<particleVolInCell.size(); cellJ++)
+    {
+        if (particleVolInCell[cellJ] > minVol)
+        {
+            Us[cellJ] /= particleVolInCell[cellJ];
+            UsDirectedStdDev[cellJ] /= particleVolInCell[cellJ];
+            for (label comp=0;comp<3;comp++)
+            {
+                UsDirectedStdDev[cellJ].component(comp) -= Us[cellJ].component(comp)*Us[cellJ].component(comp);
+                if (UsDirectedStdDev[cellJ].component(comp) > 0) UsDirectedStdDev[cellJ].component(comp) = Foam::sqrt(UsDirectedStdDev[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"));