merge TU/e changes into PFM master branch

applications/solvers/cfdemSolverMultiphase/multiphaseMixture/multiphaseMixture.C

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

applications/solvers/cfdemSolverMultiphaseScalar/Allwclean

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

applications/solvers/cfdemSolverMultiphaseScalar/Allwmake

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

applications/solvers/cfdemSolverMultiphaseScalar/EEqn.H

@@ -0,0 +1,26 @@
+        // get mixture properties
+        Cp  = mixture.Cp();
+	kf  = mixture.kf();
+        keff=particleCloud.thermCondM().thermCond(); 
+
+        // get scalar source from DEM        
+        particleCloud.energyContributions(Qsource);
+        particleCloud.energyCoefficients(QCoeff);
+
+        fvScalarMatrix EEqn
+        (
+	 rho*Cp*(fvm::ddt(voidfraction,T)
+		 + fvm::div(phi,T))
+	 - fvc::div(keff*voidfraction*fvc::grad(T))
+         ==
+	 Qsource + fvm::Sp(QCoeff,T)
+        );
+
+
+        EEqn.relax();
+        fvOptions.constrain(EEqn);
+        EEqn.solve();
+
+        particleCloud.clockM().start(31,"postFlow");
+        particleCloud.postFlow();
+        particleCloud.clockM().stop("postFlow");

applications/solvers/cfdemSolverMultiphaseScalar/Make/files

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

applications/solvers/cfdemSolverMultiphaseScalar/Make/options

@@ -0,0 +1,30 @@
+include $(CFDEM_ADD_LIBS_DIR)/additionalLibs
+
+EXE_INC = \
+	-I$(CFDEM_OFVERSION_DIR) \
+    -ImultiphaseMixture/lnInclude \
+    -I$(LIB_SRC)/transportModels \
+    -I$(LIB_SRC)/transportModels/incompressible/lnInclude \
+    -I$(LIB_SRC)/transportModels/interfaceProperties/lnInclude \
+    -I$(LIB_SRC)/TurbulenceModels/turbulenceModels/lnInclude \
+    -I$(LIB_SRC)/TurbulenceModels/incompressible/lnInclude \
+    -I$(LIB_SRC)/finiteVolume/lnInclude \
+    -I$(LIB_SRC)/meshTools/lnInclude \
+    -I$(LIB_SRC)/sampling/lnInclude \
+    -I$(CFDEM_SRC_DIR)/lagrangian/cfdemParticle/lnInclude \
+    -I$(CFDEM_SRC_DIR)/lagrangian/cfdemParticle/cfdTools \
+
+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,153 @@
+/*---------------------------------------------------------------------------*\
+License
+
+    This is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    This code is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with this code.  If not, see <http://www.gnu.org/licenses/>.
+
+    Copyright (C) 2018- Mathias Vångö, JKU Linz, Austria
+
+Application
+    cfdemSolverMultiphaseScalar
+
+Description
+    CFD-DEM solver for n incompressible fluids which captures the interfaces and
+    includes surface-tension and contact-angle effects for each phase. It is based
+    on the OpenFOAM(R)-4.x solver multiphaseInterFoam but extended to incorporate
+    DEM functionalities from the open-source DEM code LIGGGHTS.
+
+    Turbulence modelling is generic, i.e. laminar, RAS or LES may be selected.
+
+\*---------------------------------------------------------------------------*/
+
+#include "fvCFD.H"
+#include "multiphaseMixture.H"
+#include "turbulentTransportModel.H"
+#include "pimpleControl.H"
+#include "fvOptions.H"
+#include "CorrectPhi.H"
+
+#include "cfdemCloudEnergy.H"
+#include "implicitCouple.H"
+#include "clockModel.H"
+#include "smoothingModel.H"
+#include "forceModel.H"
+#include "thermCondModel.H"
+#include "energyModel.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"
+    #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"
+            
+             // --- 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;
+        }
+
+        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,246 @@
+//===============================
+// particle interaction modelling
+//===============================
+
+Info<< "\nReading momentum exchange field Ksl\n" << endl;
+volScalarField Ksl
+(
+    IOobject
+    (
+        "Ksl",
+        runTime.timeName(),
+        mesh,
+        IOobject::MUST_READ,
+        IOobject::AUTO_WRITE
+        ),
+    mesh
+    //dimensionedScalar("0", dimensionSet(1, -3, -1, 0, 0), 1.0)
+);
+
+Info<< "\nReading voidfraction field voidfraction = (Vgas/Vparticle)\n" << endl;
+volScalarField voidfraction
+(
+    IOobject
+    (
+        "voidfraction",
+        runTime.timeName(),
+        mesh,
+        IOobject::MUST_READ,
+        IOobject::AUTO_WRITE
+        ),
+    mesh
+);
+voidfraction.oldTime();
+
+Info<< "Reading particle velocity field Us\n" << endl;
+volVectorField Us
+(
+    IOobject
+    (
+        "Us",
+        runTime.timeName(),
+        mesh,
+        IOobject::MUST_READ,
+        IOobject::AUTO_WRITE
+        ),
+    mesh
+);
+
+Info<< "Reading field p_rgh\n" << endl;
+volScalarField p_rgh
+(
+    IOobject
+    (
+        "p_rgh",
+        runTime.timeName(),
+        mesh,
+        IOobject::MUST_READ,
+        IOobject::AUTO_WRITE
+    ),
+    mesh
+);
+
+Info<< "Reading field U\n" << endl;
+volVectorField U
+(
+    IOobject
+    (
+        "U",
+        runTime.timeName(),
+        mesh,
+        IOobject::MUST_READ,
+        IOobject::AUTO_WRITE
+    ),
+    mesh
+);
+
+Info<< "Reading/calculating face flux field phi\n" << endl;
+surfaceScalarField phi
+ (
+     IOobject
+     (
+        "phi",
+        runTime.timeName(),
+         mesh,
+        IOobject::READ_IF_PRESENT,
+        IOobject::AUTO_WRITE
+     ),
+     linearInterpolate(U*voidfraction) & mesh.Sf()
+ );
+
+multiphaseMixture mixture(U, phi, voidfraction);
+
+// Need to store rho for ddt(rho, U)
+volScalarField rho
+(
+    IOobject
+    (
+        "rho",
+        runTime.timeName(),
+        mesh,
+        IOobject::READ_IF_PRESENT,
+        IOobject::AUTO_WRITE
+    ),
+    mixture.rho()
+);
+rho.oldTime();
+
+//========================
+// scalar field modelling
+//========================
+Info<< "Reading/creating thermal fields\n" << endl;
+volScalarField T
+(
+    IOobject
+    (
+        "T",
+        runTime.timeName(),
+        mesh,
+        IOobject::MUST_READ,
+        IOobject::AUTO_WRITE
+    ),
+    mesh
+);
+
+volScalarField Qsource
+(
+    IOobject
+    (
+        "Qsource",
+        runTime.timeName(),
+        mesh,
+        IOobject::NO_READ,
+        IOobject::AUTO_WRITE
+    ),
+    mesh,
+    dimensionedScalar("zero", dimensionSet(1,-1,-3,0,0,0,0), 0.0)
+);
+
+volScalarField QCoeff
+(
+    IOobject
+    (
+        "Qsource",
+        runTime.timeName(),
+        mesh,
+        IOobject::NO_READ,
+        IOobject::AUTO_WRITE
+    ),
+    mesh,
+    dimensionedScalar("zero", dimensionSet(1,-1,-3,-1,0,0,0), 0.0)
+);
+
+volScalarField Cp
+(
+    IOobject
+    (
+        "Cp",
+        runTime.timeName(),
+        mesh,
+        IOobject::READ_IF_PRESENT,
+        IOobject::AUTO_WRITE
+    ),
+    mixture.Cp()
+);
+
+volScalarField kf
+(
+    IOobject
+    (
+        "kf",
+        runTime.timeName(),
+        mesh,
+        IOobject::READ_IF_PRESENT,
+        IOobject::AUTO_WRITE
+    ),
+    mixture.kf()
+);
+
+volScalarField keff
+(
+    IOobject
+    (
+        "keff",
+        runTime.timeName(),
+        mesh,
+        IOobject::NO_READ,
+        IOobject::AUTO_WRITE
+    ),
+    mesh,
+    dimensionedScalar("zero", dimensionSet(1,1,-3,-1,0,0,0), 0.0)
+);
+
+
+//========================
+
+volScalarField rhoEps ("rhoEps", rho * voidfraction);
+
+// Construct incompressible turbulence model
+autoPtr<incompressible::turbulenceModel> turbulence
+(
+    incompressible::turbulenceModel::New(U, phi, mixture)
+);
+
+
+#include "readGravitationalAcceleration.H"
+#include "readhRef.H"
+#include "gh.H"
+
+
+volScalarField p
+(
+    IOobject
+    (
+        "p",
+        runTime.timeName(),
+        mesh,
+        IOobject::NO_READ,
+        IOobject::AUTO_WRITE
+    ),
+    p_rgh + rho*gh
+);
+
+label pRefCell = 0;
+scalar pRefValue = 0.0;
+setRefCell
+(
+    p,
+    p_rgh,
+    pimple.dict(),
+    pRefCell,
+    pRefValue
+);
+
+if (p_rgh.needReference())
+{
+    p += dimensionedScalar
+    (
+        "p",
+        p.dimensions(),
+        pRefValue - getRefCellValue(p, pRefCell)
+    );
+}
+
+mesh.setFluxRequired(p_rgh.name());
+
+

applications/solvers/cfdemSolverMultiphaseScalar/multiphaseMixture/Make/files

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

applications/solvers/cfdemSolverMultiphaseScalar/multiphaseMixture/Make/options

@@ -0,0 +1,13 @@
+EXE_INC = \
+    -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
+
+LIB_LIBS = \
+    -linterfaceProperties \
+    -lincompressibleTransportModels \
+    -lfiniteVolume \
+    -lmeshTools

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

@@ -0,0 +1,146 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "alphaContactAngleFvPatchScalarField.H"
+#include "addToRunTimeSelectionTable.H"
+#include "fvPatchFieldMapper.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+alphaContactAngleFvPatchScalarField::interfaceThetaProps::interfaceThetaProps
+(
+    Istream& is
+)
+:
+    theta0_(readScalar(is)),
+    uTheta_(readScalar(is)),
+    thetaA_(readScalar(is)),
+    thetaR_(readScalar(is))
+{}
+
+
+Istream& operator>>
+(
+    Istream& is,
+    alphaContactAngleFvPatchScalarField::interfaceThetaProps& tp
+)
+{
+    is >> tp.theta0_ >> tp.uTheta_ >> tp.thetaA_ >> tp.thetaR_;
+    return is;
+}
+
+
+Ostream& operator<<
+(
+    Ostream& os,
+    const alphaContactAngleFvPatchScalarField::interfaceThetaProps& tp
+)
+{
+    os  << tp.theta0_ << token::SPACE
+        << tp.uTheta_ << token::SPACE
+        << tp.thetaA_ << token::SPACE
+        << tp.thetaR_;
+
+    return os;
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+alphaContactAngleFvPatchScalarField::alphaContactAngleFvPatchScalarField
+(
+    const fvPatch& p,
+    const DimensionedField<scalar, volMesh>& iF
+)
+:
+    zeroGradientFvPatchScalarField(p, iF)
+{}
+
+
+alphaContactAngleFvPatchScalarField::alphaContactAngleFvPatchScalarField
+(
+    const alphaContactAngleFvPatchScalarField& gcpsf,
+    const fvPatch& p,
+    const DimensionedField<scalar, volMesh>& iF,
+    const fvPatchFieldMapper& mapper
+)
+:
+    zeroGradientFvPatchScalarField(gcpsf, p, iF, mapper),
+    thetaProps_(gcpsf.thetaProps_)
+{}
+
+
+alphaContactAngleFvPatchScalarField::alphaContactAngleFvPatchScalarField
+(
+    const fvPatch& p,
+    const DimensionedField<scalar, volMesh>& iF,
+    const dictionary& dict
+)
+:
+    zeroGradientFvPatchScalarField(p, iF),
+    thetaProps_(dict.lookup("thetaProperties"))
+{
+    evaluate();
+}
+
+
+alphaContactAngleFvPatchScalarField::alphaContactAngleFvPatchScalarField
+(
+    const alphaContactAngleFvPatchScalarField& gcpsf,
+    const DimensionedField<scalar, volMesh>& iF
+)
+:
+    zeroGradientFvPatchScalarField(gcpsf, iF),
+    thetaProps_(gcpsf.thetaProps_)
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+void alphaContactAngleFvPatchScalarField::write(Ostream& os) const
+{
+    fvPatchScalarField::write(os);
+    os.writeKeyword("thetaProperties")
+        << thetaProps_ << token::END_STATEMENT << nl;
+    writeEntry("value", os);
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+makePatchTypeField
+(
+    fvPatchScalarField,
+    alphaContactAngleFvPatchScalarField
+);
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// ************************************************************************* //

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

@@ -0,0 +1,215 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::alphaContactAngleFvPatchScalarField
+
+Description
+    Contact-angle boundary condition for multi-phase interface-capturing
+    simulations.  Used in conjuction with multiphaseMixture.
+
+SourceFiles
+    alphaContactAngleFvPatchScalarField.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef alphaContactAngleFvPatchScalarField_H
+#define alphaContactAngleFvPatchScalarField_H
+
+#include "zeroGradientFvPatchFields.H"
+#include "multiphaseMixture.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                           Class alphaContactAngleFvPatch Declaration
+\*---------------------------------------------------------------------------*/
+
+class alphaContactAngleFvPatchScalarField
+:
+    public zeroGradientFvPatchScalarField
+{
+public:
+
+    class interfaceThetaProps
+    {
+        //- Equilibrium contact angle
+        scalar theta0_;
+
+        //- Dynamic contact angle velocity scale
+        scalar uTheta_;
+
+        //- Limiting advancing contact angle
+        scalar thetaA_;
+
+        //- Limiting receeding contact angle
+        scalar thetaR_;
+
+
+    public:
+
+        // Constructors
+            interfaceThetaProps()
+            {}
+
+            interfaceThetaProps(Istream&);
+
+
+        // Member functions
+
+            //- Return the equilibrium contact angle theta0
+            scalar theta0(bool matched=true) const
+            {
+                if (matched) return theta0_;
+                else return 180.0 - theta0_;
+            }
+
+            //- Return the dynamic contact angle velocity scale
+            scalar uTheta() const
+            {
+                return uTheta_;
+            }
+
+            //- Return the limiting advancing contact angle
+            scalar thetaA(bool matched=true) const
+            {
+                if (matched) return thetaA_;
+                else return 180.0 - thetaA_;
+            }
+
+            //- Return the limiting receeding contact angle
+            scalar thetaR(bool matched=true) const
+            {
+                if (matched) return thetaR_;
+                else return 180.0 - thetaR_;
+            }
+
+
+        // IO functions
+
+            friend Istream& operator>>(Istream&, interfaceThetaProps&);
+            friend Ostream& operator<<(Ostream&, const interfaceThetaProps&);
+    };
+
+    typedef HashTable
+    <
+        interfaceThetaProps,
+        multiphaseMixture::interfacePair,
+        multiphaseMixture::interfacePair::hash
+    > thetaPropsTable;
+
+
+private:
+
+    // Private data
+
+        thetaPropsTable thetaProps_;
+
+
+public:
+
+    //- Runtime type information
+    TypeName("alphaContactAngle");
+
+
+    // Constructors
+
+        //- Construct from patch and internal field
+        alphaContactAngleFvPatchScalarField
+        (
+            const fvPatch&,
+            const DimensionedField<scalar, volMesh>&
+        );
+
+        //- Construct from patch, internal field and dictionary
+        alphaContactAngleFvPatchScalarField
+        (
+            const fvPatch&,
+            const DimensionedField<scalar, volMesh>&,
+            const dictionary&
+        );
+
+        //- Construct by mapping given alphaContactAngleFvPatchScalarField
+        //  onto a new patch
+        alphaContactAngleFvPatchScalarField
+        (
+            const alphaContactAngleFvPatchScalarField&,
+            const fvPatch&,
+            const DimensionedField<scalar, volMesh>&,
+            const fvPatchFieldMapper&
+        );
+
+        //- Construct and return a clone
+        virtual tmp<fvPatchScalarField> clone() const
+        {
+            return tmp<fvPatchScalarField>
+            (
+                new alphaContactAngleFvPatchScalarField(*this)
+            );
+        }
+
+        //- Construct as copy setting internal field reference
+        alphaContactAngleFvPatchScalarField
+        (
+            const alphaContactAngleFvPatchScalarField&,
+            const DimensionedField<scalar, volMesh>&
+        );
+
+        //- Construct and return a clone setting internal field reference
+        virtual tmp<fvPatchScalarField> clone
+        (
+            const DimensionedField<scalar, volMesh>& iF
+        ) const
+        {
+            return tmp<fvPatchScalarField>
+            (
+                new alphaContactAngleFvPatchScalarField(*this, iF)
+            );
+        }
+
+
+    // Member functions
+
+        //- Return the contact angle properties
+        const thetaPropsTable& thetaProps() const
+        {
+            return thetaProps_;
+        }
+
+        //- Write
+        virtual void write(Ostream&) const;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

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

@@ -0,0 +1 @@
+../alphaContactAngle/alphaContactAngleFvPatchScalarField.C

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

@@ -0,0 +1 @@
+../alphaContactAngle/alphaContactAngleFvPatchScalarField.H

applications/solvers/cfdemSolverMultiphaseScalar/multiphaseMixture/lnInclude/multiphaseMixture.C

@@ -0,0 +1 @@
+../multiphaseMixture.C

applications/solvers/cfdemSolverMultiphaseScalar/multiphaseMixture/lnInclude/multiphaseMixture.H

@@ -0,0 +1 @@
+../multiphaseMixture.H

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

@@ -0,0 +1 @@
+../phase/phase.C

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

@@ -0,0 +1 @@
+../phase/phase.H

applications/solvers/cfdemSolverMultiphaseScalar/multiphaseMixture/multiphaseMixture.C

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

applications/solvers/cfdemSolverMultiphaseScalar/multiphaseMixture/multiphaseMixture.H

@@ -0,0 +1,290 @@
+/*---------------------------------------------------------------------------*\
+License
+
+    This is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    This code is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with this code.  If not, see <http://www.gnu.org/licenses/>.
+
+    Copyright (C) 2018- Mathias Vångö, JKU Linz, Austria
+
+Class
+    multiphaseMixture
+
+Description
+    This class is based on the OpenFOAM(R) Foam::multiphaseMixture class,
+    which is an incompressible multi-phase mixture with built in solution
+    for the phase fractions with interface compression for interface-capturing.
+    It has been extended to include the void fraction in the volume fraction
+    transport equations.
+
+    Derived from transportModel so that it can be unsed in conjunction with
+    the incompressible turbulence models.
+
+    Surface tension and contact-angle is handled for the interface
+    between each phase-pair.
+
+SourceFiles
+    multiphaseMixture.C
+\*---------------------------------------------------------------------------*/
+
+#ifndef multiphaseMixture_H
+#define multiphaseMixture_H
+
+#include "incompressible/transportModel/transportModel.H"
+#include "IOdictionary.H"
+#include "phase.H"
+#include "PtrDictionary.H"
+#include "volFields.H"
+#include "surfaceFields.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                      Class multiphaseMixture Declaration
+\*---------------------------------------------------------------------------*/
+
+class multiphaseMixture
+:
+    public IOdictionary,
+    public transportModel
+{
+public:
+
+    class interfacePair
+    :
+        public Pair<word>
+    {
+    public:
+
+        class hash
+        :
+            public Hash<interfacePair>
+        {
+        public:
+
+            hash()
+            {}
+
+            label operator()(const interfacePair& key) const
+            {
+                return word::hash()(key.first()) + word::hash()(key.second());
+            }
+        };
+
+
+        // Constructors
+
+            interfacePair()
+            {}
+
+            interfacePair(const word& alpha1Name, const word& alpha2Name)
+            :
+                Pair<word>(alpha1Name, alpha2Name)
+            {}
+
+            interfacePair(const phase& alpha1, const phase& alpha2)
+            :
+                Pair<word>(alpha1.name(), alpha2.name())
+            {}
+
+
+        // Friend Operators
+
+            friend bool operator==
+            (
+                const interfacePair& a,
+                const interfacePair& b
+            )
+            {
+                return
+                (
+                    ((a.first() == b.first()) && (a.second() == b.second()))
+                 || ((a.first() == b.second()) && (a.second() == b.first()))
+                );
+            }
+
+            friend bool operator!=
+            (
+                const interfacePair& a,
+                const interfacePair& b
+            )
+            {
+                return (!(a == b));
+            }
+    };
+
+
+private:
+
+    // Private data
+
+        //- Dictionary of phases
+        PtrDictionary<phase> phases_;
+
+        const fvMesh& mesh_;
+        const volVectorField& U_;
+        const surfaceScalarField& phi_;
+        const volScalarField& voidfraction_;
+        surfaceScalarField rhoPhi_;
+        surfaceScalarField surfaceTensionForce_;
+        volScalarField alphas_;
+
+        volScalarField nu_;
+
+        typedef HashTable<scalar, interfacePair, interfacePair::hash>
+            sigmaTable;
+
+        sigmaTable sigmas_;
+        dimensionSet dimSigma_;
+
+        //- Stabilisation for normalisation of the interface normal
+        const dimensionedScalar deltaN_;
+
+        //- Conversion factor for degrees into radians
+        static const scalar convertToRad;
+
+
+    // Private member functions
+
+        void calcAlphas();
+        
+        tmp<volScalarField> calcNu() const;
+
+        void solveAlphas(const scalar cAlpha);
+
+        tmp<surfaceVectorField> nHatfv
+        (
+            const volScalarField& alpha1,
+            const volScalarField& alpha2
+        ) const;
+
+        tmp<surfaceScalarField> nHatf
+        (
+            const volScalarField& alpha1,
+            const volScalarField& alpha2
+        ) const;
+
+        void correctContactAngle
+        (
+            const phase& alpha1,
+            const phase& alpha2,
+            surfaceVectorField::Boundary& nHatb
+        ) const;
+
+        tmp<volScalarField> K(const phase& alpha1, const phase& alpha2) const;
+        tmp<surfaceScalarField> calcStf() const;
+
+public:
+
+    // Constructors
+
+        //- Construct from components
+        multiphaseMixture
+        (
+            const volVectorField& U,
+            const surfaceScalarField& phi,
+            const volScalarField& voidfraction
+        );
+
+
+    //- Destructor
+    virtual ~multiphaseMixture()
+    {}
+
+
+    // Member Functions
+
+        //- Return the phases
+        const PtrDictionary<phase>& phases() const
+        {
+            return phases_;
+        }
+
+        //- Return the velocity
+        const volVectorField& U() const
+        {
+            return U_;
+        }
+
+        //- Return the volumetric flux
+        const surfaceScalarField& phi() const
+        {
+            return phi_;
+        }
+
+        const surfaceScalarField& rhoPhi() const
+        {
+            return rhoPhi_;
+        }
+
+        //- Return the mixture density
+        tmp<volScalarField> rho() const;
+
+        //- Return the mixture density for patch
+        tmp<scalarField> rho(const label patchi) const;
+
+        //- Return the dynamic laminar viscosity
+        tmp<volScalarField> mu() const;
+
+        //- Return the dynamic laminar viscosity for patch
+        tmp<scalarField> mu(const label patchi) const;
+
+        //- Return the face-interpolated dynamic laminar viscosity
+        tmp<surfaceScalarField> muf() const;
+
+        //- Return the kinematic laminar viscosity
+        tmp<volScalarField> nu() const;
+
+        //- Return the laminar viscosity for patch
+        tmp<scalarField> nu(const label patchi) const;
+
+        //- Return the face-interpolated dynamic laminar viscosity
+        tmp<surfaceScalarField> nuf() const;
+
+        //- Return the heat capacity
+        tmp<volScalarField> Cp() const;
+
+        //- Return the thermal conductivity
+        tmp<volScalarField> kf() const;
+
+        tmp<surfaceScalarField> surfaceTensionForce() const
+        {
+            return surfaceTensionForce_;
+        }
+        
+        //- Indicator of the proximity of the interface
+        //  Field values are 1 near and 0 away for the interface.
+        tmp<volScalarField> nearInterface() const;
+
+        //- Solve for the mixture phase-fractions
+        void solve();
+
+        //- Correct the mixture properties
+        void correct();
+
+        //- Read base transportProperties dictionary
+        bool read();
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

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

@@ -0,0 +1,103 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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_)
+{}
+
+
+// * * * * * * * * * * * * * * * 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_;
+
+    if (nuModel_->read(phaseDict_))
+    {
+        phaseDict_.lookup("rho") >> rho_;
+
+        return true;
+    }
+    else
+    {
+        return false;
+    }
+}
+
+
+// ************************************************************************* //

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

@@ -0,0 +1,176 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2011-2015 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::phase
+
+Description
+    Single incompressible phase derived from the phase-fraction.
+    Used as part of the multiPhaseMixture for interface-capturing multi-phase
+    simulations.
+
+SourceFiles
+    phase.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef phase_H
+#define phase_H
+
+#include "volFields.H"
+#include "dictionaryEntry.H"
+#include "incompressible/viscosityModels/viscosityModel/viscosityModel.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                           Class phase Declaration
+\*---------------------------------------------------------------------------*/
+
+class phase
+:
+    public volScalarField
+{
+    // Private data
+
+        word name_;
+        dictionary phaseDict_;
+        autoPtr<viscosityModel> nuModel_;
+        dimensionedScalar rho_;
+        dimensionedScalar Cp_;
+        dimensionedScalar kf_;
+
+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_;
+        }
+
+        //- 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/pEqn.H

@@ -21,7 +21,7 @@ phi = phiHbyA + rAUf*(fvc::interpolate(Ksl/rho) * phiS);
 MRF.makeRelative(phi);
 
 // adjustment of phi (only in cases w.o. p boundary conditions) not tested yet
-adjustPhi(phi, U, p);
+//adjustPhi(phi, U, p);
 
 // Update the pressure BCs to ensure flux consistency
 constrainPressure(p, Uvoidfraction, phiHbyA, rAUvoidfraction, MRF);

etc/bashrc

@@ -218,7 +218,7 @@ alias cfdemCompCFDEMuti='bash $CFDEM_PROJECT_DIR/etc/compileCFDEMcoupling_uti.sh
 alias cfdemTestTUT='bash $CFDEM_PROJECT_DIR/etc/testTutorials.sh'
 
 #- shortcut to visualize the clock model data
-alias vizClock='python $CFDEM_UT_DIR/vizClock/matPlot.py'
+alias vizClock='python2 $CFDEM_UT_DIR/vizClock/matPlot.py'
 
 #- recursive touch of current directory
 alias touchRec='find ./* -exec touch {} \;'

etc/solver-list.txt

@@ -6,3 +6,4 @@ cfdemSolverIB/dir
 cfdemSolverPisoScalar/dir
 cfdemSolverRhoPimpleChem/dir
 cfdemSolverMultiphase/dir
+cfdemSolverMultiphaseScalar/dir

src/finiteVolume/fields/fvPatchFields/derived/uniformFixedValueTube/uniformFixedValueTubeFvPatchField.C

@@ -48,7 +48,8 @@ uniformFixedValueTubeFvPatchField<Type>::uniformFixedValueTubeFvPatchField
     densityFieldName_("rho"),
     tubeLength_(-1.),
     tubeDiameter_(-1.),
-    zeta_(-1)
+    zeta_(-1),
+    p0_(-1)
 {
     Info << "error, wrong constructor1!" << endl;
 }
@@ -71,7 +72,8 @@ uniformFixedValueTubeFvPatchField<Type>::uniformFixedValueTubeFvPatchField
     densityFieldName_("rho"),
     tubeLength_(ptf.tubeLength_),
     tubeDiameter_(ptf.tubeDiameter_),
-    zeta_(ptf.zeta_)
+    zeta_(ptf.zeta_),
+    p0_(ptf.p0_)
 {
     const scalar t = this->db().time().timeOutputValue();
     fvPatchField<Type>::operator==(uniformValue_->value(t));
@@ -89,13 +91,14 @@ uniformFixedValueTubeFvPatchField<Type>::uniformFixedValueTubeFvPatchField
 :
     fixedValueFvPatchField<Type>(p, iF),
     uniformValue_(Function1<Type>::New("uniformValue", dict)),
-	pName_("p"), //JOKER
+	pName_("p_rgh"), //JOKER
 	phiName_("phi"), //JOKER
     velocityFieldName_("U"),
     densityFieldName_("rho"),
     tubeLength_(readScalar(dict.lookup("tubeLength"))),
     tubeDiameter_(readScalar(dict.lookup("tubeDiameter"))),
-    zeta_(readScalar(dict.lookup("zeta")))
+    zeta_(readScalar(dict.lookup("zeta"))),
+    p0_(readScalar(dict.lookup("p0")))
 {
     const scalar t = this->db().time().timeOutputValue();
     fvPatchField<Type>::operator==(uniformValue_->value(t));
@@ -110,13 +113,14 @@ uniformFixedValueTubeFvPatchField<Type>::uniformFixedValueTubeFvPatchField
 :
     fixedValueFvPatchField<Type>(ptf),
     uniformValue_(ptf.uniformValue_().clone().ptr()),
-	pName_("p"), //JOKER
+	pName_("p_rgh"), //JOKER
 	phiName_("phi"), //JOKER
     velocityFieldName_("U"),
     densityFieldName_("rho"),
     tubeLength_(ptf.tubeLength_),
     tubeDiameter_(ptf.tubeDiameter_),
-    zeta_(ptf.zeta_)
+    zeta_(ptf.zeta_),
+    p0_(ptf.p0_)
 {
     const scalar t = this->db().time().timeOutputValue();
     fvPatchField<Type>::operator==(uniformValue_->value(t));
@@ -132,13 +136,14 @@ uniformFixedValueTubeFvPatchField<Type>::uniformFixedValueTubeFvPatchField
 :
     fixedValueFvPatchField<Type>(ptf, iF),
     uniformValue_(ptf.uniformValue_().clone().ptr()),
-	pName_("p"), //JOKER
+	pName_("p_rgh"), //JOKER
 	phiName_("phi"), //JOKER
     velocityFieldName_("U"),
     densityFieldName_("rho"),
     tubeLength_(ptf.tubeLength_),
     tubeDiameter_(ptf.tubeDiameter_),
-    zeta_(ptf.zeta_)
+    zeta_(ptf.zeta_),
+    p0_(ptf.p0_)
 {
     const scalar t = this->db().time().timeOutputValue();
     fvPatchField<Type>::operator==(uniformValue_->value(t));
@@ -185,7 +190,7 @@ void uniformFixedValueTubeFvPatchField<Type>::updateCoeffs()
     // calc pressure drop
 //    fvPatchField<Type>::operator==(0.5*zeta*density*mag(velocity)*mag(velocity)*uniformValue_->value(t));
 	//dP = zeta * l/d * rho u^2 / 2
-	fvPatchField<Type>::operator==((pressure+(0.5*zeta_*tubeLength_/tubeDiameter_*density*phip/this->patch().magSf()*phip/this->patch().magSf()-pressure)*uRelFact)*uniformValue_->value(t));
+	fvPatchField<Type>::operator==((pressure+(0.5*zeta_*tubeLength_/tubeDiameter_*density*phip/this->patch().magSf()*phip/this->patch().magSf()+p0_-pressure)*uRelFact)*uniformValue_->value(t));
 	//fvPatchField<Type>::operator==((100.*uniformValue_->value(t)));
     fixedValueFvPatchField<Type>::updateCoeffs();
 }
@@ -200,6 +205,7 @@ void uniformFixedValueTubeFvPatchField<Type>::write(Ostream& os) const
     os.writeKeyword("tubeLength") << tubeLength_ << token::END_STATEMENT << nl;
     os.writeKeyword("tubeDiameter") << tubeDiameter_ << token::END_STATEMENT << nl;
     os.writeKeyword("zeta") << zeta_ << token::END_STATEMENT << nl;
+    os.writeKeyword("p0") << p0_ << token::END_STATEMENT << nl;
 }
 
 

src/finiteVolume/fields/fvPatchFields/derived/uniformFixedValueTube/uniformFixedValueTubeFvPatchField.H

@@ -74,6 +74,8 @@ class uniformFixedValueTubeFvPatchField
 
         scalar zeta_;
 
+        scalar p0_;
+
 public:
 
     //- Runtime type information

src/lagrangian/cfdemParticle/Make/files

@@ -41,6 +41,8 @@ $(energyModels)/energyModel/energyModel.C
 $(energyModels)/energyModel/newEnergyModel.C
 $(energyModels)/heatTransferGunn/heatTransferGunn.C
 $(energyModels)/heatTransferGunnPartField/heatTransferGunnPartField.C
+$(energyModels)/wallHeatTransferYagi/wallHeatTransferYagi.C
+$(energyModels)/wallHeatTransferYagiImplicit/wallHeatTransferYagiImplicit.C
 $(energyModels)/reactionHeat/reactionHeat.C
 
 $(thermCondModels)/thermCondModel/thermCondModel.C
@@ -65,6 +67,7 @@ $(forceModels)/KochHillDrag/KochHillDrag.C
 $(forceModels)/KochHillRWDrag/KochHillRWDrag.C
 $(forceModels)/LaEuScalarTemp/LaEuScalarTemp.C
 $(forceModels)/virtualMassForce/virtualMassForce.C
+$(forceModels)/ParmarBassetForce/ParmarBassetForce.C
 $(forceModels)/gradPForce/gradPForce.C
 $(forceModels)/viscForce/viscForce.C
 $(forceModels)/MeiLift/MeiLift.C
@@ -174,5 +177,6 @@ $(smoothingModels)/smoothingModel/newSmoothingModel.C
 $(smoothingModels)/noSmoothing/noSmoothing.C
 $(smoothingModels)/constDiffSmoothing/constDiffSmoothing.C
 $(smoothingModels)/temporalSmoothing/temporalSmoothing.C
+$(smoothingModels)/constDiffAndTemporalSmoothing/constDiffAndTemporalSmoothing.C
 
 LIB = $(CFDEM_LIB_DIR)/lib$(CFDEM_LIB_NAME)

src/lagrangian/cfdemParticle/cfdemCloud/cfdemCloud.C

@@ -87,6 +87,7 @@ cfdemCloud::cfdemCloud
     getParticleDensities_(couplingProperties_.lookupOrDefault<bool>("getParticleDensities",false)),
     getParticleEffVolFactors_(couplingProperties_.lookupOrDefault<bool>("getParticleEffVolFactors",false)),
     getParticleTypes_(couplingProperties_.lookupOrDefault<bool>("getParticleTypes",false)),
+    getParticleAngVels_(couplingProperties_.lookupOrDefault<bool>("getParticleAngVels",false)),
     maxDEMForce_(0.),
     modelType_(couplingProperties_.lookup("modelType")),
     positions_(NULL),
@@ -103,6 +104,7 @@ cfdemCloud::cfdemCloud
     particleDensities_(NULL),
     particleEffVolFactors_(NULL),
     particleTypes_(NULL),
+    particleAngVels_(NULL),
     particleWeights_(NULL),
     particleVolumes_(NULL),
     particleV_(NULL),
@@ -390,6 +392,7 @@ cfdemCloud::~cfdemCloud()
     if(getParticleDensities_) dataExchangeM().destroy(particleDensities_,1);
     if(getParticleEffVolFactors_) dataExchangeM().destroy(particleEffVolFactors_,1);
     if(getParticleTypes_) dataExchangeM().destroy(particleTypes_,1);
+    if(getParticleAngVels_) dataExchangeM().destroy(particleAngVels_,1);
 }
 
 // * * * * * * * * * * * * * * * private Member Functions  * * * * * * * * * * * * * //
@@ -405,6 +408,7 @@ void cfdemCloud::getDEMdata()
     if(getParticleDensities_) dataExchangeM().getData("density","scalar-atom",particleDensities_);
     if(getParticleEffVolFactors_) dataExchangeM().getData("effvolfactor","scalar-atom",particleEffVolFactors_);
     if(getParticleTypes_) dataExchangeM().getData("type","scalar-atom",particleTypes_);
+    if(getParticleAngVels_) dataExchangeM().getData("omega","vector-atom",particleAngVels_);
 }
 
 void cfdemCloud::giveDEMdata()
@@ -448,7 +452,12 @@ void cfdemCloud::setForces()
     resetArray(expForces_,numberOfParticles(),3);
     resetArray(DEMForces_,numberOfParticles(),3);
     resetArray(Cds_,numberOfParticles(),1);
-    for (int i=0;i<cfdemCloud::nrForceModels();i++) cfdemCloud::forceM(i).setForce();
+    for (int i=0;i<cfdemCloud::nrForceModels();i++)
+    {
+        clockM().start(30+i,forceModels_[i]);
+        cfdemCloud::forceM(i).setForce();
+        clockM().stop(forceModels_[i]);
+    }
 
     if (limitDEMForces_)
     {
@@ -669,6 +678,7 @@ bool cfdemCloud::evolve
         //CHECK JUST TIME-INTERPOLATE ALREADY SMOOTHENED VOIDFRACTIONNEXT AND UsNEXT FIELD
         //      IMPLICIT FORCE CONTRIBUTION AND SOLVER USE EXACTLY THE SAME AVERAGED
         //      QUANTITIES AT THE GRID!
+
         const scalar timeStepFrac = dataExchangeM().timeStepFraction();
         int old_precision = Info().precision(10);
         Info << "\n timeStepFraction() = " << timeStepFrac << endl;
@@ -756,6 +766,7 @@ bool cfdemCloud::reAllocArrays()
         if(getParticleDensities_) dataExchangeM().allocateArray(particleDensities_,0.,1);
         if(getParticleEffVolFactors_) dataExchangeM().allocateArray(particleEffVolFactors_,0.,1);
         if(getParticleTypes_) dataExchangeM().allocateArray(particleTypes_,0,1);
+        if(getParticleAngVels_) dataExchangeM().allocateArray(particleAngVels_,0.,3);
         arraysReallocated_ = true;
         return true;
     }

src/lagrangian/cfdemParticle/cfdemCloud/cfdemCloud.H

@@ -102,6 +102,10 @@ protected:
 
     const bool getParticleTypes_;
 
+    bool getParticleDensities_;
+
+    bool getParticleAngVels_;
+
     scalar maxDEMForce_;
 
     const word modelType_;
@@ -134,6 +138,8 @@ protected:
 
     int **particleTypes_;
 
+    double **particleAngVels_;
+
     double **particleWeights_;
 
     double **particleVolumes_;
@@ -351,7 +357,11 @@ public:
        virtual inline int ** particleTypes() const;
        virtual inline label particleType(label index) const;
 
-       //access to the particle's rotation and torque data
+       inline bool getParticleAngVels() const;
+       virtual inline double ** particleAngVels() const;
+       virtual vector particleAngVel(label index) const; 
+
+       //access to the particle's and torque data
        virtual inline double ** DEMTorques() const {return NULL;}
        virtual inline double ** omegaArray() const {return NULL;}
        virtual vector omega(int) const {return vector::zero;}

src/lagrangian/cfdemParticle/cfdemCloud/cfdemCloudI.H

@@ -268,6 +268,27 @@ inline label cfdemCloud::particleType(label index) const
     }
 }
 
+inline bool cfdemCloud::getParticleAngVels() const
+{
+    return getParticleAngVels_;
+}
+	  
+inline double ** cfdemCloud::particleAngVels() const
+{
+    return particleAngVels_;
+}
+
+inline vector cfdemCloud::particleAngVel(label index) const
+{
+    if(!getParticleAngVels_) return vector::zero;
+    else
+    {
+    	vector angVel;
+    	for(int i=0;i<3;i++) angVel[i] = particleAngVels_[index][i];
+        return angVel;
+    }
+}
+
 inline int cfdemCloud::numberOfParticles() const
 {
     return numberOfParticles_;

src/lagrangian/cfdemParticle/subModels/energyModel/heatTransferGunn/heatTransferGunn.C

@@ -146,6 +146,11 @@ heatTransferGunn::heatTransferGunn
     partNu_(NULL),
     scaleDia_(1.),
     typeCG_(propsDict_.lookupOrDefault<scalarList>("coarseGrainingFactors",scalarList(1,1.0)))
+    multiphase_(propsDict_.lookupOrDefault<bool>("multiphase",false)),
+    kfFieldName_(propsDict_.lookupOrDefault<word>("kfFieldName",voidfractionFieldName_)), // use voidfractionField as dummy to prevent lookup error when not using multiphase
+    kfField_(sm.mesh().lookupObject<volScalarField> (kfFieldName_)),
+    CpFieldName_(propsDict_.lookupOrDefault<word>("CpFieldName",voidfractionFieldName_)), // use voidfractionField as dummy to prevent lookup error when not using multiphase
+    CpField_(sm.mesh().lookupObject<volScalarField> (CpFieldName_))
 {
     allocateMyArrays();
 
@@ -341,7 +346,15 @@ void heatTransferGunn::calcEnergyContribution()
                 ds = 2.*particleCloud_.radius(index);
                 ds_scaled = ds/cg;
                 muf = mufField[cellI];
+
                 Rep = ds_scaled * magUr * voidfraction * rho_[cellI]/ muf;
+
+                if(multiphase_)
+                {
+                	kf0_ = kfField_[cellI];
+                	Cp_  = CpField_[cellI];
+                }
+
                 Pr = max(SMALL, Cp_ * muf / kf0_);
 
                 Nup = Nusselt(voidfraction, Rep, Pr);

src/lagrangian/cfdemParticle/subModels/energyModel/heatTransferGunn/heatTransferGunn.H

@@ -116,6 +116,16 @@ protected:
 
     scalarList typeCG_;
 
+    bool multiphase_;
+
+    word kfFieldName_;
+
+    const volScalarField& kfField_;
+
+    word CpFieldName_;
+
+    const volScalarField& CpField_;
+
     void allocateMyArrays() const;
 
     void partTempField();

src/lagrangian/cfdemParticle/subModels/energyModel/wallHeatTransferYagi/wallHeatTransferYagi.C

@@ -0,0 +1,329 @@
+/*---------------------------------------------------------------------------*\
+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
+
+\*---------------------------------------------------------------------------*/
+
+#include "error.H"
+#include "wallHeatTransferYagi.H"
+#include "addToRunTimeSelectionTable.H"
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+  // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+  defineTypeNameAndDebug(wallHeatTransferYagi, 0);
+
+  addToRunTimeSelectionTable(energyModel, wallHeatTransferYagi, dictionary);
+
+  // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+
+  wallHeatTransferYagi::wallHeatTransferYagi
+  (
+    const dictionary& dict,
+    cfdemCloudEnergy& sm
+  )
+  :
+  energyModel(dict,sm),
+  propsDict_(dict.subDict(typeName + "Props")),
+  interpolation_(propsDict_.lookupOrDefault<bool>("interpolation",false)),
+  verbose_(propsDict_.lookupOrDefault<bool>("verbose",false)),
+  QWallFluidName_(propsDict_.lookupOrDefault<word>("QWallFluidName","QWallFluid")),
+  QWallFluid_
+  (   IOobject
+      (
+          QWallFluidName_,
+          sm.mesh().time().timeName(),
+          sm.mesh(),
+          IOobject::READ_IF_PRESENT,
+          IOobject::AUTO_WRITE
+      ),
+      sm.mesh(),
+      dimensionedScalar("zero", dimensionSet(1,-1,-3,0,0,0,0), 0.0)
+  ),
+  wallTempName_(propsDict_.lookup("wallTempName")),
+  wallTemp_
+  (   IOobject
+    (
+      wallTempName_,
+      sm.mesh().time().timeName(),
+      sm.mesh(),
+      IOobject::READ_IF_PRESENT,
+      IOobject::AUTO_WRITE
+    ),
+    sm.mesh(),
+    dimensionedScalar("zero", dimensionSet(0,0,0,1,0,0,0), 0.0)
+  ),
+  ReField_
+  (   IOobject
+    (
+      "ReField",
+      sm.mesh().time().timeName(),
+      sm.mesh(),
+      IOobject::NO_READ,
+      IOobject::NO_WRITE
+    ),
+    sm.mesh(),
+    dimensionedScalar("zero", dimensionSet(0,0,0,0,0,0,0), 0.0)
+  ),
+  PrField_
+  (   IOobject
+    (
+      "PrField",
+      sm.mesh().time().timeName(),
+      sm.mesh(),
+      IOobject::NO_READ,
+      IOobject::NO_WRITE
+    ),
+    sm.mesh(),
+    dimensionedScalar("zero", dimensionSet(0,0,0,0,0,0,0), 0.0)
+  ),
+  tempFieldName_(propsDict_.lookupOrDefault<word>("tempFieldName","T")),
+  tempField_(sm.mesh().lookupObject<volScalarField> (tempFieldName_)),
+  voidfractionFieldName_(propsDict_.lookupOrDefault<word>("voidfractionFieldName","voidfraction")),
+  voidfraction_(sm.mesh().lookupObject<volScalarField> (voidfractionFieldName_)),
+  voidfractionMax_(readScalar(propsDict_.lookup("voidfractionMax"))),
+  maxSource_(1e30),
+  velFieldName_(propsDict_.lookupOrDefault<word>("velFieldName","U")),
+  U_(sm.mesh().lookupObject<volVectorField> (velFieldName_)),
+  UsFieldName_(propsDict_.lookup("granVelFieldName")),
+  Us_(sm.mesh().lookupObject<volVectorField> (UsFieldName_)),
+  densityFieldName_(propsDict_.lookupOrDefault<word>("densityFieldName","rho")),
+  rho_(sm.mesh().lookupObject<volScalarField> (densityFieldName_)),
+  partRe_(NULL),
+  multiphase_(propsDict_.lookupOrDefault<bool>("multiphase",false)),
+  kfFieldName_(propsDict_.lookupOrDefault<word>("kfFieldName",voidfractionFieldName_)), // use voidfractionField as dummy to prevent lookup error when not using multiphase
+  kfField_(sm.mesh().lookupObject<volScalarField> (kfFieldName_)),
+  CpFieldName_(propsDict_.lookupOrDefault<word>("CpFieldName",voidfractionFieldName_)), // use voidfractionField as dummy to prevent lookup error when not using multiphase
+  CpField_(sm.mesh().lookupObject<volScalarField> (CpFieldName_))
+  {
+    allocateMyArrays();
+
+    if (propsDict_.found("maxSource"))
+    {
+        maxSource_=readScalar(propsDict_.lookup ("maxSource"));
+        Info << "limiting wall source field to: " << maxSource_ << endl;
+    }
+
+    if (verbose_)
+    {
+      ReField_.writeOpt() = IOobject::AUTO_WRITE;
+      PrField_.writeOpt() = IOobject::AUTO_WRITE;
+      ReField_.write();
+      PrField_.write();
+    }
+    
+    // currently it is detected if field was auto generated or defined
+    // improvement would be changing the type here automatically
+    forAll(wallTemp_.boundaryField(),patchI)
+        if(wallTemp_.boundaryField()[patchI].type()=="calculated")
+            FatalError <<"Scalar field:"<< wallTemp_.name() << " must be defined.\n" << abort(FatalError);
+
+    wallTemp_.writeOpt() = IOobject::AUTO_WRITE;
+  }
+
+  // * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+  wallHeatTransferYagi::~wallHeatTransferYagi()
+  {
+    particleCloud_.dataExchangeM().destroy(partRe_,1);
+  }
+
+  // * * * * * * * * * * * * * * * private Member Functions  * * * * * * * * * * * * * //
+  void wallHeatTransferYagi::allocateMyArrays() const
+  {
+    // get memory for 2d arrays
+    double initVal=0.0;
+
+    if(verbose_)
+    {
+      particleCloud_.dataExchangeM().allocateArray(partRe_,initVal,1);
+    }
+  }
+
+  // * * * * * * * * * * * * * * * * Member Fct  * * * * * * * * * * * * * * * //
+
+  void wallHeatTransferYagi::calcEnergyContribution()
+  {
+    allocateMyArrays();
+
+    // reset Scalar field
+    QWallFluid_.primitiveFieldRef() = 0.0;
+
+    #ifdef compre
+    const volScalarField mufField = particleCloud_.turbulence().mu();
+    #else
+    const volScalarField mufField = particleCloud_.turbulence().nu()*rho_;
+    #endif
+
+    interpolationCellPoint<scalar> voidfractionInterpolator_(voidfraction_);
+    interpolationCellPoint<vector> UInterpolator_(U_);
+    interpolationCellPoint<scalar> TInterpolator_(tempField_);
+
+    // calculate Rep
+    for(int index = 0;index < particleCloud_.numberOfParticles(); ++index)
+    {
+      label cellI = particleCloud_.cellIDs()[index][0];
+      if(cellI >= 0)
+      {
+        scalar voidfraction;
+        vector Ufluid;
+
+        if(interpolation_)
+        {
+          vector position = particleCloud_.position(index);
+          voidfraction = voidfractionInterpolator_.interpolate(position,cellI);
+          Ufluid = UInterpolator_.interpolate(position,cellI);
+        }
+        else
+        {
+          voidfraction = voidfraction_[cellI];
+          Ufluid = U_[cellI];
+        }
+
+        if (voidfraction < 0.01)
+          voidfraction = 0.01;
+
+        vector Us = particleCloud_.velocity(index);
+        scalar magUr = mag(Ufluid - Us);
+        scalar ds = 2.*particleCloud_.radius(index);
+        scalar muf = mufField[cellI];
+
+        scalar Rep = ds * magUr * voidfraction * rho_[cellI]/ muf;
+        partRe_[index][0] = Rep;
+      }
+    }
+
+    // calculate Rep field
+    particleCloud_.averagingM().resetWeightFields();
+    particleCloud_.averagingM().setScalarAverage
+    (
+      ReField_,
+      partRe_,
+      particleCloud_.particleWeights(),
+      particleCloud_.averagingM().UsWeightField(),
+      NULL
+    );
+
+    // calculate Pr field
+    PrField_ = CpField_ * mufField / kfField_;
+
+    const fvPatchList& patches = U_.mesh().boundary();
+
+    // calculate flux
+    forAll(patches, patchi)
+    {
+      const fvPatch& curPatch = patches[patchi];
+
+      if (wallTemp_.boundaryField().types()[patchi] == "fixedValue")
+      {
+        if(tempField_.boundaryField().types()[patchi] == "zeroGradient")
+        {
+
+          forAll(curPatch, facei)
+          {							
+            label faceCelli = curPatch.faceCells()[facei];
+
+            // calculate Urel
+            scalar magG = mag(U_[faceCelli]-Us_[faceCelli])*voidfraction_[faceCelli]*rho_[faceCelli];
+
+            // calculate H
+            scalar H;
+            if (voidfraction_[faceCelli]<=voidfractionMax_)
+              H = 0.2087 * (pow(ReField_[faceCelli]+SMALL,-0.20)) * CpField_[faceCelli] * magG / (pow(PrField_[faceCelli],2/3) + SMALL);
+            else
+              H = 0;
+
+            // get delta T (wall-fluid)
+            scalar Twall  = wallTemp_.boundaryField()[patchi][facei];
+            scalar Tfluid = tempField_[faceCelli];
+            scalar deltaT = Twall - Tfluid;
+
+            // get area
+            scalar area = curPatch.magSf()[facei];
+
+            // calculate heat flux
+            heatFlux(faceCelli, H, area, Twall, Tfluid);
+            heatFluxCoeff(faceCelli, H, area);
+
+            if(verbose_ && facei >=0 && facei <2)
+            {
+              Info << "####################" << endl;
+              Info << "cellID: " << faceCelli << endl;
+              Info << "G : " << magG << endl;
+              Info << "Re: " << ReField_[faceCelli] << endl;
+              Info << "Pr: " << PrField_[faceCelli] << endl;
+              Info << "Cp: " << CpField_[faceCelli] << endl;
+              Info << "kf: " << kfField_[faceCelli] << endl;
+              Info << "H : " << H << endl;
+              Info << "Twall: " << Twall << endl;
+              Info << "Tfluid: " << Tfluid << endl;
+              Info << "dT: " << deltaT << endl;
+              Info << "q: " << H*deltaT << endl;
+              Info << "area: " << area << endl;
+              Info << "Q:" << H*deltaT*area << endl;
+            }
+          }		    
+        }
+        else
+        {
+            FatalError << "wallHeatTransferYagi requires zeroGradient BC for temperature field" << endl;
+        }		    
+      }	    
+    }
+
+    QWallFluid_.primitiveFieldRef() /= QWallFluid_.mesh().V();
+
+    // limit source term
+    forAll(QWallFluid_,cellI)
+    {
+        scalar EuFieldInCell = QWallFluid_[cellI];
+
+        if(mag(EuFieldInCell) > maxSource_ )
+        {
+             Pout << "limiting source term"  << endl  ;
+             QWallFluid_[cellI] = sign(EuFieldInCell) * maxSource_;
+        }
+    }
+
+    QWallFluid_.correctBoundaryConditions();
+
+  }
+
+  void wallHeatTransferYagi::addEnergyContribution(volScalarField& Qsource) const
+  {
+    Qsource += QWallFluid_;
+  }
+
+  void wallHeatTransferYagi::heatFlux(label faceCelli, scalar H, scalar area, scalar Twall, scalar Tfluid)
+  {
+    QWallFluid_[faceCelli] += H * area * (Twall - Tfluid);
+  }
+
+  void wallHeatTransferYagi::heatFluxCoeff(label faceCelli, scalar H, scalar area)
+  {
+    //no heat transfer coefficient in explicit model
+  }
+
+  // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// ************************************************************************* //

src/lagrangian/cfdemParticle/subModels/energyModel/wallHeatTransferYagi/wallHeatTransferYagi.H

@@ -0,0 +1,144 @@
+/*---------------------------------------------------------------------------*\
+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
+
+    Description
+    Correlation for wall-to-bed heat transfer according to
+    Yagi, S. A.I.Ch.E. Journal 5.1 (1959)
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef wallHeatTransferYagi_H
+#define wallHeatTransferYagi_H
+
+#include "fvCFD.H"
+#include "cfdemCloudEnergy.H"
+#include "energyModel.H"
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                           Class wallHeatTransferYagi Declaration
+\*---------------------------------------------------------------------------*/
+
+class wallHeatTransferYagi
+:
+    public energyModel
+{
+protected:
+
+    dictionary propsDict_;
+
+    bool interpolation_;
+
+    bool verbose_;
+
+    word QWallFluidName_;
+        
+    volScalarField QWallFluid_;
+
+    word wallTempName_;
+
+    volScalarField wallTemp_;
+
+    volScalarField ReField_;
+
+    volScalarField PrField_;
+
+    word tempFieldName_;
+
+    const volScalarField& tempField_;
+
+    word voidfractionFieldName_;
+
+    const volScalarField& voidfraction_;
+
+    scalar voidfractionMax_;
+
+    scalar maxSource_;                           // max (limited) value of src field
+
+    word velFieldName_;
+
+    const volVectorField& U_;
+
+    word UsFieldName_;
+
+    const volVectorField& Us_;
+
+    word densityFieldName_;
+
+    const volScalarField& rho_;
+
+    mutable double **partRe_;
+
+    bool multiphase_;
+
+    word kfFieldName_;
+
+    const volScalarField& kfField_;
+
+    word CpFieldName_;
+
+    const volScalarField& CpField_;
+
+    void allocateMyArrays() const;
+
+    virtual void heatFlux(label, scalar, scalar, scalar, scalar);
+
+    virtual void heatFluxCoeff(label, scalar, scalar);
+
+public:
+
+    //- Runtime type information
+    TypeName("wallHeatTransferYagi");
+
+    // Constructors
+
+        //- Construct from components
+        wallHeatTransferYagi
+        (
+            const dictionary& dict,
+            cfdemCloudEnergy& sm
+        );
+
+
+    // Destructor
+
+        virtual ~wallHeatTransferYagi();
+
+
+    // Member Functions
+
+    void addEnergyContribution(volScalarField&) const;
+
+	void addEnergyCoefficient(volScalarField&) const {}
+
+	void calcEnergyContribution();
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/lagrangian/cfdemParticle/subModels/energyModel/wallHeatTransferYagiImplicit/wallHeatTransferYagiImplicit.C

@@ -0,0 +1,116 @@
+/*---------------------------------------------------------------------------*\
+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
+
+\*---------------------------------------------------------------------------*/
+
+#include "error.H"
+#include "wallHeatTransferYagiImplicit.H"
+#include "addToRunTimeSelectionTable.H"
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+defineTypeNameAndDebug(wallHeatTransferYagiImplicit, 0);
+
+addToRunTimeSelectionTable(energyModel, wallHeatTransferYagiImplicit, dictionary);
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+// Construct from components
+wallHeatTransferYagiImplicit::wallHeatTransferYagiImplicit
+(
+    const dictionary& dict,
+    cfdemCloudEnergy& sm
+)
+:
+    wallHeatTransferYagi(dict,sm),
+    QWallFluidCoeffName_(propsDict_.lookupOrDefault<word>("QWallFluidCoeffName","QWallFluidCoeff")),
+    QWallFluidCoeff_
+    (   IOobject
+        (
+            QWallFluidCoeffName_,
+            sm.mesh().time().timeName(),
+            sm.mesh(),
+            IOobject::READ_IF_PRESENT,
+            IOobject::AUTO_WRITE
+        ),
+        sm.mesh(),
+        dimensionedScalar("zero", dimensionSet(1,-1,-3,-1,0,0,0), 0.0)
+    ),
+    partHeatFluxCoeff_(NULL)
+{
+    allocateMyArrays();
+
+    // no limiting necessary for implicit heat transfer
+    maxSource_ = 1e30;
+}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+wallHeatTransferYagiImplicit::~wallHeatTransferYagiImplicit()
+{
+    particleCloud_.dataExchangeM().destroy(partHeatFluxCoeff_,1);
+}
+
+// * * * * * * * * * * * * * * * private Member Functions  * * * * * * * * * * * * * //
+void wallHeatTransferYagiImplicit::allocateMyArrays() const
+{
+//    wallHeatTransferYagi::allocateMyArrays();
+    double initVal=0.0;
+    particleCloud_.dataExchangeM().allocateArray(partHeatFluxCoeff_,initVal,1);
+}
+// * * * * * * * * * * * * * * * * Member Fct  * * * * * * * * * * * * * * * //
+
+void wallHeatTransferYagiImplicit::calcEnergyContribution()
+{
+    allocateMyArrays();
+
+    QWallFluidCoeff_.primitiveFieldRef() = 0.0;
+
+    wallHeatTransferYagi::calcEnergyContribution();
+
+    QWallFluidCoeff_.primitiveFieldRef() /= QWallFluidCoeff_.mesh().V();
+
+//    QWallFluidCoeff_.correctBoundaryConditions();
+
+}
+
+void wallHeatTransferYagiImplicit::addEnergyCoefficient(volScalarField& Qsource) const
+{
+    Qsource += QWallFluidCoeff_;
+}
+
+void wallHeatTransferYagiImplicit::heatFlux(label faceCelli, scalar H, scalar area, scalar Twall, scalar Tfluid)
+{
+    QWallFluid_[faceCelli] += H * area * Twall;
+}
+
+void wallHeatTransferYagiImplicit::heatFluxCoeff(label faceCelli, scalar H, scalar area)
+{
+    QWallFluidCoeff_[faceCelli] -= H * area;
+}
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// ************************************************************************* //

src/lagrangian/cfdemParticle/subModels/energyModel/wallHeatTransferYagiImplicit/wallHeatTransferYagiImplicit.H

@@ -0,0 +1,95 @@
+/*---------------------------------------------------------------------------*\
+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
+    
+    Description
+    Correlation for wall-to-bed heat transfer according to
+    Yagi, S. A.I.Ch.E. Journal 5.1 (1959)
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef wallHeatTransferYagiImplicit_H
+#define wallHeatTransferYagiImplicit_H
+
+#include "fvCFD.H"
+#include "cfdemCloudEnergy.H"
+#include "wallHeatTransferYagi.H"
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                           Class wallHeatTransferYagiImplicit Declaration
+\*---------------------------------------------------------------------------*/
+
+class wallHeatTransferYagiImplicit
+:
+    public wallHeatTransferYagi
+{
+private:
+    
+    word QWallFluidCoeffName_;
+        
+    volScalarField QWallFluidCoeff_;
+
+    mutable double **partHeatFluxCoeff_;              // Lagrangian array
+
+    void allocateMyArrays() const;
+
+    virtual void heatFlux(label, scalar, scalar, scalar, scalar);
+
+    virtual void heatFluxCoeff(label, scalar, scalar);
+       
+public:
+
+    //- Runtime type information
+    TypeName("wallHeatTransferYagiImplicit");
+
+    // Constructors
+
+        //- Construct from components
+        wallHeatTransferYagiImplicit
+        (
+            const dictionary& dict,
+            cfdemCloudEnergy& sm
+        );
+
+
+    // Destructor
+
+        virtual ~wallHeatTransferYagiImplicit();
+
+
+    // Member Functions
+	
+	void addEnergyCoefficient(volScalarField&) const;
+	
+	void calcEnergyContribution();
+
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/lagrangian/cfdemParticle/subModels/forceModel/BeetstraDrag/BeetstraDrag.C

@@ -70,14 +70,21 @@ BeetstraDrag::BeetstraDrag
     k_(0.05),
     useGC_(false),
     usePC_(false)
+    polydisperse_(propsDict_.lookupOrDefault<bool>("polydisperse",false))
 {
+
+    if(polydisperse_)
+	Info << "Drag model: using polydisperse correction factor \n";
+
+
     //Append the field names to be probed
     particleCloud_.probeM().initialize(typeName, typeName+".logDat");
     particleCloud_.probeM().vectorFields_.append("dragForce"); //first entry must  be the force
     particleCloud_.probeM().vectorFields_.append("Urel");
     particleCloud_.probeM().scalarFields_.append("Rep");
-    particleCloud_.probeM().scalarFields_.append("betaP");
     particleCloud_.probeM().scalarFields_.append("voidfraction");
+    particleCloud_.probeM().scalarFields_.append("Fdrag");
+    particleCloud_.probeM().scalarFields_.append("y_polydisperse");
     particleCloud_.probeM().writeHeader();
 
     // init force sub model
@@ -159,7 +166,7 @@ void BeetstraDrag::setForce() const
     scalar voidfraction(1);
     vector Ufluid(0,0,0);
     vector drag(0,0,0);
-    label cellI=0;
+    label  cellI = 0;
 
     vector Us(0,0,0);
     vector Ur(0,0,0);
@@ -177,12 +184,14 @@ void BeetstraDrag::setForce() const
 
     scalar cg = typeCG_[0];
     label partType = 1;
+    scalar yi(0);
 
     vector dragExplicit(0,0,0);
     scalar dragCoefficient(0);
 
     interpolationCellPoint<scalar> voidfractionInterpolator_(voidfraction_);
     interpolationCellPoint<vector> UInterpolator_(U_);
+    interpolationCellPoint<scalar> dSauterInterpolator_(dSauterField_);
 
     #include "setupProbeModel.H"
 
@@ -205,6 +214,7 @@ void BeetstraDrag::setForce() const
                     Ufluid       = UInterpolator_.interpolate(position,cellI);
                     //Ensure interpolated void fraction to be meaningful
                     // Info << " --> voidfraction: " << voidfraction << endl;
+
                     if (voidfraction > 1.00) voidfraction = 1.0;
                     if (voidfraction < minVoidfraction_) voidfraction = minVoidfraction_;
                 }
@@ -223,6 +233,7 @@ void BeetstraDrag::setForce() const
                     scaleDia3 = cg*cg*cg;
                 }
 
+<<<<<<< HEAD
                 Us = particleCloud_.velocity(index);
                 Ur = Ufluid-Us;
                 magUr = mag(Ur);
@@ -243,6 +254,13 @@ void BeetstraDrag::setForce() const
                                    *3*M_PI*nuf*rho*voidfraction
                                    *effDiameter(ds_scaled, cellI, index)
                                    *scaleDia3*scaleDrag_;
+		if (polydisperse_)
+                {
+                    // Beetstra et. al (2007), Eq. (21)
+                    yi = ds_scaled/dSauter;
+                    dragCoefficient *= (voidfraction*yi + 0.064*voidfraction*yi*yi*yi);
+                }
+		
 
                 // calculate filtering corrections
                 if (useGC_)
@@ -274,6 +292,7 @@ void BeetstraDrag::setForce() const
                 {
                     Pout << "cellI = " << cellI << endl;
                     Pout << "index = " << index << endl;
+                    Pout << "Ufluid = " << Ufluid << endl;
                     Pout << "Us = " << Us << endl;
                     Pout << "Ur = " << Ur << endl;
                     Pout << "ds = " << ds << endl;
@@ -295,6 +314,8 @@ void BeetstraDrag::setForce() const
                     vValues.append(Ur);
                     sValues.append(Rep);
                     sValues.append(voidfraction);
+                    sValues.append(Fdrag);
+                    sValues.append(yi);
                     particleCloud_.probeM().writeProbe(index, sValues, vValues);
                 }
             }

src/lagrangian/cfdemParticle/subModels/forceModel/BeetstraDrag/BeetstraDrag.H

@@ -102,7 +102,7 @@ protected:
 
     double h(double) const;
 
-
+    mutable bool polydisperse_;
 
 public:
 

src/lagrangian/cfdemParticle/subModels/forceModel/MeiLift/MeiLift.C

@@ -81,9 +81,29 @@ MeiLift::MeiLift
     propsDict_(dict.subDict(typeName + "Props")),
     velFieldName_(propsDict_.lookup("velFieldName")),
     U_(sm.mesh().lookupObject<volVectorField> (velFieldName_)),
-    useSecondOrderTerms_(false)
+	useShearInduced_(propsDict_.lookupOrDefault<bool>("useShearInduced",true)),
+	useSpinInduced_(propsDict_.lookupOrDefault<bool>("useSpinInduced",false)),
+	combineShearSpin_(propsDict_.lookupOrDefault<bool>("combineShearSpin",false))
 {
-    if (propsDict_.found("useSecondOrderTerms")) useSecondOrderTerms_=true;
+	// read switches
+
+	if(useShearInduced_)
+	    Info << "Lift model: including shear-induced term.\n";
+
+	if(useSpinInduced_)
+    {
+	    Info << "Lift model: including spin-induced term.\n";
+	    Info << "Make sure to use a rolling friction model in LIGGGHTS!\n";
+	    if(!dict.lookupOrDefault<bool>("getParticleAngVels",false))
+	    	FatalError << "Lift model: useSpinInduced=true requires getParticleAngVels=true in couplingProperties" << abort(FatalError);
+	}
+
+	if(combineShearSpin_)
+	{
+	    Info << "Lift model: combining shear- and spin-terms by assuming equilibrium spin-rate.\n";
+	    if(!useShearInduced_ || !useSpinInduced_)
+	        FatalError << "Shear- and spin-induced lift must be activated in order to combine." << abort(FatalError);
+	}
 
     // init force sub model
     setForceSubModels(propsDict_);
@@ -100,10 +120,15 @@ MeiLift::MeiLift
     particleCloud_.probeM().initialize(typeName, typeName+".logDat");
     particleCloud_.probeM().vectorFields_.append("liftForce"); //first entry must the be the force
     particleCloud_.probeM().vectorFields_.append("Urel");      //other are debug
-    particleCloud_.probeM().vectorFields_.append("vorticity"); //other are debug
-    particleCloud_.probeM().scalarFields_.append("Rep");       //other are debug
-    particleCloud_.probeM().scalarFields_.append("Rew");       //other are debug
-    particleCloud_.probeM().scalarFields_.append("J_star");    //other are debug
+    particleCloud_.probeM().vectorFields_.append("vorticity");
+    particleCloud_.probeM().vectorFields_.append("Ang_velocity");
+    particleCloud_.probeM().scalarFields_.append("Rep");
+    particleCloud_.probeM().scalarFields_.append("Rew");
+    particleCloud_.probeM().scalarFields_.append("J*");
+    particleCloud_.probeM().scalarFields_.append("Cl(shear)");
+    particleCloud_.probeM().scalarFields_.append("Cl*(spin)");
+    particleCloud_.probeM().scalarFields_.append("Omega_eq");
+    particleCloud_.probeM().scalarFields_.append("Cl(combined)");
     particleCloud_.probeM().writeHeader();
 }
 
@@ -122,27 +147,44 @@ void MeiLift::setForce() const
     const volScalarField& nufField = forceSubM(0).nuField();
     const volScalarField& rhoField = forceSubM(0).rhoField();
 
+    // vectors
     vector position(0,0,0);
     vector lift(0,0,0);
     vector Us(0,0,0);
     vector Ur(0,0,0);
+    vector Omega(0,0,0);
+    vector vorticity(0,0,0);
+
+    // properties
     scalar magUr(0);
     scalar magVorticity(0);
+    scalar magOmega(0);
     scalar ds(0);
     scalar nuf(0);
     scalar rho(0);
     scalar voidfraction(1);
+
+    // dimensionless groups
     scalar Rep(0);
     scalar Rew(0);
-    scalar Cl(0);
-    scalar Cl_star(0);
-    scalar J_star(0);
-    scalar Omega_eq(0);
-    scalar alphaStar(0);
-    scalar epsilonSqr(0.0);
-    scalar epsilon(0);
+
+    // shear induced
     scalar omega_star(0);
-    vector vorticity(0,0,0);
+    scalar Clshear(0);
+    scalar J_star(0);
+    scalar alphaStar(0);
+    scalar epsilonSqr(0);
+    scalar epsilon(0);
+
+    // spin induced
+    scalar Omega_star(0);
+    scalar Clspin_star(0);
+
+    // shear-spin combination
+    scalar Omega_eq(0);
+    scalar Clcombined(0);
+
+
     volVectorField vorticityField = fvc::curl(U_);
 
     interpolationCellPoint<vector> UInterpolator_(U_);
@@ -152,14 +194,14 @@ void MeiLift::setForce() const
 
     for (int index = 0; index < particleCloud_.numberOfParticles(); ++index)
     {
-        //if(mask[index][0])
-        //{
             lift = vector::zero;
             label cellI = particleCloud_.cellIDs()[index][0];
 
             if (cellI > -1) // particle Found
             {
+                // properties
             	Us = particleCloud_.velocity(index);
+            	Omega = particleCloud_.particleAngVel(index);
 
                 if (forceSubM(0).interpolation())
                 {
@@ -173,24 +215,24 @@ void MeiLift::setForce() const
                     vorticity = vorticityField[cellI];
                 }
 
-                magUr        = mag(Ur);
-                magVorticity = mag(vorticity);
-
-                if (magUr > 0 && magVorticity > 0)
-                {
                 ds  	= 2. * particleCloud_.radius(index);
                 nuf 	= nufField[cellI];
                 rho 	= rhoField[cellI];
 
-                    // calc dimensionless properties
+                magUr 	= mag(Ur);
                 Rep 	= ds*magUr/nuf;
-                    Rew = magVorticity*ds*ds/nuf;
 
+                // shear-induced lift
+                if (useShearInduced_)
+                {
+                	magVorticity 	= mag(vorticity);
+                	Rew 			= magVorticity*ds*ds/nuf;
                 	omega_star 		= magVorticity * ds / magUr;
                     alphaStar 		= 0.5 * omega_star;
                     epsilonSqr 		= omega_star / Rep;
                     epsilon 		= sqrt(epsilonSqr);
 
+
                     //Basic model for the correction to the Saffman lift
                     //McLaughlin (1991), Mei (1992), Loth and Dorgan (2009)
                     //J_star = 0.443 * J
@@ -212,35 +254,76 @@ void MeiLift::setForce() const
                                  * (1.0   + tanh(2.5 * (log10(epsilon) + 0.191)))
                                  * (0.667 + tanh(6.0 * (      epsilon  - 0.32 )));
                     }
-                    //Loth and Dorgan (2009), Eq (31): Saffman lift-coefficient: ClSaff = 12.92 / pi * epsilon ~ 4.11 * epsilon
-                    //Loth and Dorgan (2009), Eq (32)
-                    Cl = J_star * 4.11 * epsilon; //multiply correction to the basic Saffman model
 
-                    //Second order terms given by Loth and Dorgan (2009)
-                    if (useSecondOrderTerms_)
-                    {
-                        scalar sqrtRep = sqrt(Rep);
-                        //Loth and Dorgan (2009), Eq (34)
-                        Cl_star = 1.0 - (0.675 + 0.15 * (1.0 + tanh(0.28 * (alphaStar - 2.0)))) * tanh(0.18 * sqrtRep);
-                        //Loth and Dorgan (2009), Eq (38)
-                        Omega_eq = alphaStar * (1.0 - 0.0075 * Rew) * (1.0 - 0.062 * sqrtRep - 0.001 * Rep);
-                        //Loth and Dorgan (2009), Eq (39)
-                        Cl += Omega_eq * Cl_star;
+                    //Loth and Dorgan (2009), Eq (31), Eq (32)
+                    Clshear = J_star * 4.11 * epsilon; //multiply correction to the basic Saffman model
                 }
 
+                if (useSpinInduced_)
+                {
+                	magOmega 		= mag(Omega);
+                	Omega_star 		= magOmega * ds / magUr;
+
+                	//Loth and Dorgan (2009), Eq (34)
+                	Clspin_star		= 1.0 - (0.675 + 0.15 * (1.0 + tanh(0.28 * (Omega_star - 2.0)))) * tanh(0.18 * sqrt(Rep));
+                }
+
+                if (combineShearSpin_)
+                {
+                    //Loth and Dorgan (2009), Eq (38)
+                    Omega_eq 		= alphaStar * (1.0 - 0.0075 * Rew) * (1.0 - 0.062 * sqrt(Rep) - 0.001 * Rep);
+                    //Loth and Dorgan (2009), Eq (39)
+                    Clcombined 		= Clshear + Clspin_star * Omega_eq;
+
+                    if (magUr>0.0 && magVorticity>0.0)
+                    {
                     	//Loth and Dorgan (2009), Eq (27)
+                         // please note: in Loth and Dorgan Vrel = Vp-Uf, here Urel = Uf-Vp. Hence the reversed cross products.
                     	lift = 0.125 * constant::mathematical::pi
                     			* rho
-                           * Cl
-                           * magUr * Ur ^ vorticity / magVorticity
+								* Clcombined
+								* magUr * magUr
+								* (Ur ^ vorticity) / mag(Ur ^ vorticity) // force direction
 								* ds * ds;
+                    }
+                }
+                else
+                {
+                	//Loth and Dorgan (2009), Eq (36)
+                    // please note: in Loth and Dorgan Vrel = Vp-Uf, here Urel = Uf-Vp. Hence the reversed cross products.
+                	if (useShearInduced_)
+                	{
+                        if (magUr>0.0 && magVorticity>0.0)
+                        {
+                        	//Loth and Dorgan (2009), Eq (27)
+                        	lift += 0.125 * constant::mathematical::pi
+                        			* rho
+									* Clshear
+									* magUr * magUr
+                                    * (Ur ^ vorticity) / mag(Ur ^ vorticity) // force direction
+									* ds * ds;
+                        }
+                	}
+                	if (useSpinInduced_)
+                	{
+                        if (magUr>0.0 && magOmega>0.0)
+                        {
+                        	//Loth and Dorgan (2009), Eq (33)
+                        	lift += 0.125 * constant::mathematical::pi
+                        			* rho
+									* Clspin_star
+									* (Ur ^ Omega)
+									* ds * ds * ds;
+                        }
+                	}
+                }
 
                 if (modelType_ == "B")
                 {
                     voidfraction = particleCloud_.voidfraction(index);
                     lift /= voidfraction;
                 }
-                }
+
 
                 //**********************************
                 //SAMPLING AND VERBOSE OUTOUT
@@ -250,6 +333,7 @@ void MeiLift::setForce() const
                     Pout << "Us = " << Us << endl;
                     Pout << "Ur = " << Ur << endl;
                     Pout << "vorticity = " << vorticity << endl;
+                    Pout << "Omega = " << Omega << endl;
                     Pout << "ds = " << ds << endl;
                     Pout << "rho = " << rho << endl;
                     Pout << "nuf = " << nuf << endl;
@@ -258,6 +342,10 @@ void MeiLift::setForce() const
                     Pout << "alphaStar = " << alphaStar << endl;
                     Pout << "epsilon = " << epsilon << endl;
                     Pout << "J_star = " << J_star << endl;
+                    Pout << "Omega_eq = " << Omega_eq << endl;
+                    Pout << "Clshear = " <<  Clshear<< endl;
+                    Pout << "Clspin_star = " << Clspin_star << endl;
+                    Pout << "Clcombined = " << Clcombined << endl;
                     Pout << "lift = " << lift << endl;
                 }
 
@@ -268,9 +356,14 @@ void MeiLift::setForce() const
                     vValues.append(lift);   //first entry must the be the force
                     vValues.append(Ur);
                     vValues.append(vorticity);
+                    vValues.append(Omega);
                     sValues.append(Rep);
                     sValues.append(Rew);
                     sValues.append(J_star);
+                    sValues.append(Clshear);
+                    sValues.append(Clspin_star);
+                    sValues.append(Omega_eq);
+                    sValues.append(Clcombined);
                     particleCloud_.probeM().writeProbe(index, sValues, vValues);
                 }
                 // END OF SAMPLING AND VERBOSE OUTOUT
@@ -279,7 +372,6 @@ void MeiLift::setForce() const
             }
             // write particle based data to global array
             forceSubM(0).partToArray(index,lift,vector::zero);
-        //}
     }
 
 }

src/lagrangian/cfdemParticle/subModels/forceModel/MeiLift/MeiLift.H

@@ -83,7 +83,11 @@ private:
 
     const volVectorField& U_;
 
-    bool useSecondOrderTerms_;
+    bool useShearInduced_;
+
+    bool useSpinInduced_;
+
+    bool combineShearSpin_;
 
 public:
 

src/lagrangian/cfdemParticle/subModels/forceModel/ParmarBassetForce/ParmarBassetForce.C

@@ -0,0 +1,626 @@
+/*---------------------------------------------------------------------------*\
+    CFDEMcoupling - Open Source CFD-DEM coupling
+
+    CFDEMcoupling is part of the CFDEMproject
+    www.cfdem.com
+                                Christoph Goniva, christoph.goniva@cfdem.com
+                                Copyright 2009-2012 JKU Linz
+                                Copyright 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, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+Description
+    This code is designed to realize coupled CFD-DEM simulations using LIGGGHTS
+    and OpenFOAM(R). Note: this code is not part of OpenFOAM(R) (see DISCLAIMER).
+\*---------------------------------------------------------------------------*/
+
+#include "error.H"
+
+#include "ParmarBassetForce.H"
+#include "addToRunTimeSelectionTable.H"
+#include "smoothingModel.H"
+#include "constDiffSmoothing.H"
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+#define NOTONCPU 9999
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+defineTypeNameAndDebug(ParmarBassetForce, 0);
+
+addToRunTimeSelectionTable
+(
+    forceModel,
+    ParmarBassetForce,
+    dictionary
+);
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+// Construct from components
+ParmarBassetForce::ParmarBassetForce
+(
+    const dictionary& dict,
+    cfdemCloud& sm
+)
+:
+    forceModel(dict,sm),
+    propsDict_(dict.subDict(typeName + "Props")),
+    velFieldName_(propsDict_.lookup("velFieldName")),
+    U_(sm.mesh().lookupObject<volVectorField> (velFieldName_)),
+	UsFieldName_(propsDict_.lookup("granVelFieldName")),
+	Us_(sm.mesh().lookupObject<volVectorField> (UsFieldName_)),
+	nInt_(readScalar(propsDict_.lookup("nIntegral"))),
+	discOrder_(readScalar(propsDict_.lookup("discretisationOrder"))),
+	nHist_(nInt_+discOrder_+1),
+	ddtUrelHist_(nHist_,NULL), 		               // UrelHist_[ndt in past][particle ID][dim]
+	rHist_(nHist_,NULL),                           // rHist_[ndt in past][particle ID][0]
+	FHist_(2,List<double**>(2*discOrder_+1,NULL)), // FHist_[k={1,2}-1][ndt in past][particle ID][dim]
+	gH0_(NULL),
+	tRef_(NULL),
+	mRef_(NULL),
+	lRef_(NULL),
+    Urel_
+    (   IOobject
+        (
+            "Urel",
+            sm.mesh().time().timeName(),
+            sm.mesh(),
+            IOobject::NO_READ,
+            IOobject::AUTO_WRITE
+        ),
+        sm.mesh(),
+        dimensionedVector("zero", dimensionSet(0,1,-1,0,0,0,0), vector(0,0,0))
+    ),
+    ddtUrel_
+    (   IOobject
+        (
+            "ddtUrel",
+            sm.mesh().time().timeName(),
+            sm.mesh(),
+            IOobject::NO_READ,
+            IOobject::AUTO_WRITE
+        ),
+        sm.mesh(),
+        dimensionedVector("zero", dimensionSet(0,1,-2,0,0,0,0), vector(0,0,0))
+    ),
+    smoothingModel_
+    (
+        smoothingModel::New
+        (
+            propsDict_,
+            sm
+        )
+    )
+{
+
+    // init force sub model
+    setForceSubModels(propsDict_);
+    // define switches which can be read from dict
+    forceSubM(0).setSwitchesList(SW_TREAT_FORCE_EXPLICIT,true); // activate treatExplicit switch
+    forceSubM(0).setSwitchesList(SW_INTERPOLATION,true); // activate search for interpolate switch
+    forceSubM(0).readSwitches();
+
+    //Extra switches/settings
+    particleCloud_.checkCG(true);
+
+    if (discOrder_!=1 || discOrder_!=2)
+    	FatalError << "Parmar Basset Force: Discretisation order > 2 not implemented!" << endl;
+
+    //Append the field names to be probed
+    particleCloud_.probeM().initialize(typeName, typeName+".logDat");
+    particleCloud_.probeM().vectorFields_.append("ParmarBassetForce"); //first entry must the be the force
+    particleCloud_.probeM().vectorFields_.append("Urel");
+    particleCloud_.probeM().vectorFields_.append("ddtUrel");
+    //
+    particleCloud_.probeM().vectorFields_.append("UrelNoSmooth");
+    particleCloud_.probeM().vectorFields_.append("ddtUrelNoSmooth");
+    //
+    particleCloud_.probeM().vectorFields_.append("Fshort");
+    particleCloud_.probeM().vectorFields_.append("Flong1");
+    particleCloud_.probeM().vectorFields_.append("Flong2");
+    particleCloud_.probeM().scalarFields_.append("ReRef");
+    particleCloud_.probeM().scalarFields_.append("tRef");
+    particleCloud_.probeM().scalarFields_.append("mRef");
+    particleCloud_.probeM().scalarFields_.append("lRef");
+    particleCloud_.probeM().scalarFields_.append("r");
+    particleCloud_.probeM().scalarFields_.append("t0");
+    particleCloud_.probeM().scalarFields_.append("nKnown");
+    particleCloud_.probeM().writeHeader();
+
+    Urel_ = Us_ - U_;
+}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+ParmarBassetForce::~ParmarBassetForce()
+{
+	particleCloud_.dataExchangeM().destroy(gH0_,    1);
+    particleCloud_.dataExchangeM().destroy(tRef_,   1);
+    particleCloud_.dataExchangeM().destroy(mRef_,   1);
+    particleCloud_.dataExchangeM().destroy(lRef_,   1);
+
+    for (int i=0; i<nHist_; i++)
+    {
+    	particleCloud_.dataExchangeM().destroy(ddtUrelHist_[i],3);
+    	particleCloud_.dataExchangeM().destroy(rHist_      [i],1);
+    }
+
+    for (int i=0; i<3; i++)
+    	for (int k=0; k<2; k++)
+    		particleCloud_.dataExchangeM().destroy(FHist_[k][i],3);
+
+}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+void ParmarBassetForce::setForce() const
+{
+	// allocate arrays
+	if(particleCloud_.numberOfParticlesChanged())
+    	reAllocArrays();
+
+    vector position(0,0,0);
+    vector Urel(0,0,0);
+    vector ddtUrel(0,0,0);
+
+    scalar t0min    = 0.00;
+
+    const volScalarField& nufField = forceSubM(0).nuField();
+    const volScalarField& rhoField = forceSubM(0).rhoField();
+
+    double c[2][4][3]   = {{{0.004045,0.011788,-0.851409},
+    		                {0.184350,0.009672,-0.425408},
+			                {0.000746,0.014390,-1.262288},
+			                {0.021708,0.016611,-0.867352}},
+
+    		               {{0.536783,0.005453,-1.430163},
+    				        {1.683212,0.003598,-0.773872},
+					        {0.316350,0.004367,-1.423649},
+					        {0.420423,0.000677,-0.670551}}
+    					   };
+
+    double chi[2][4][2] = {{{1.464787,  0.422501},
+    		                {0.785093, -0.078301},
+			                {0.585753, -0.118737},
+			                {-0.077407,-0.016665}},
+
+    		               {{1.001435, -0.003566},
+    				        {0.415345, -0.136697},
+					        {0.536055, -0.205041},
+					        {-0.032671,-0.021047}}
+                           };
+
+    #include "setupProbeModel.H"
+
+    Urel_ = Us_ - U_;
+    ddtUrel_ = fvc::ddt(Us_) - fvc::ddt(U_) - (Us_ & fvc::grad(U_));
+
+    //
+    volVectorField UrelNoSmooth_ = Urel_;
+    volVectorField ddtUrelNoSmooth_ = ddtUrel_;
+    //
+
+    smoothingM().smoothen(Urel_);
+    smoothingM().smoothen(ddtUrel_);
+
+    interpolationCellPoint<vector> UrelInterpolator_(Urel_);
+    interpolationCellPoint<vector> ddtUrelInterpolator_(ddtUrel_);
+
+    //
+    interpolationCellPoint<vector> UrelNoSmoothInterpolator_(UrelNoSmooth_);
+    interpolationCellPoint<vector> ddtUrelNoSmoothInterpolator_(ddtUrelNoSmooth_);
+    //
+
+    for(int index = 0;index <  particleCloud_.numberOfParticles(); index++)
+    {
+            vector ParmarBassetForce(0,0,0);
+            vector Fshort(0,0,0);
+            label  cellI = particleCloud_.cellIDs()[index][0];
+
+            if (cellI > -1) // particle Found
+            {
+            	// particle position (m)
+            	position = particleCloud_.position(index);
+
+            	// relative velocity (m/s)
+            	if(forceSubM(0).interpolation())
+            		Urel = UrelInterpolator_.interpolate(position,cellI);
+            	else
+            		Urel = Urel_[cellI];
+
+            	// acceleration (m/s2)
+            	if(forceSubM(0).interpolation())
+            		ddtUrel = ddtUrelInterpolator_.interpolate(position,cellI);
+            	else
+            		ddtUrel = ddtUrel_[cellI];
+
+            	//********* set/get reference point *********//
+
+                scalar rs = particleCloud_.radius(index);
+                scalar Re = 2.*rs*mag(Urel)/nufField[cellI];
+                scalar gH = (0.75 + 0.105*Re)/(Re+SMALL); // Eq. 3.3
+                scalar r;
+
+                if (gH0_[index][0]!=NOTONCPU)
+				{
+                	r = pow(gH0_[index][0]/gH,1.5); // Eq. 3.4
+
+					if (r<.25 || r>2.)
+					{
+						gH0_[index][0]            = NOTONCPU; //reset reference
+            			ddtUrelHist_[0][index][0] = NOTONCPU; //reset ddtU history (only component used for checking nKnown)
+					}
+
+				}
+
+                if (gH0_[index][0]==NOTONCPU)
+                {
+                	// calculate new reference values
+                	scalar tRef = (rs*rs) / nufField[cellI] * (gH*gH) * 4.3354; // (256/pi)^(1/3) = 4.3354 (Eq 3.1)
+                	scalar Vs   = rs*rs*rs*M_PI*4/3;
+                	scalar mRef = Vs*rhoField[cellI] * gH * 5.2863;             // 9/(2*sqrt(pi))*(256/pi)^(1/6) = 5.2863 (Eq. 3.2)
+
+                	 gH0_[index][0] = gH;
+                	tRef_[index][0] = tRef;
+                	mRef_[index][0] = mRef;
+                	lRef_[index][0] = rs;
+					r = 1;
+                }
+
+            	scalar mps    = lRef_[index][0]/tRef_[index][0]; // m/s
+            	scalar mpss   = mps/tRef_[index][0]; // m/s^2
+            	scalar newton = mpss*mRef_[index][0]; // kg*m/s^2
+            	scalar dt     = U_.mesh().time().deltaT().value() / tRef_[index][0];  // dim.less
+            	scalar t0     = nInt_*dt; // dim.less
+
+            	//********* update histories *********//
+
+            	// non-dimensionlise
+            	Urel    /= mps;
+            	ddtUrel /= mpss;
+
+            	// update ddtUrel  and r history
+            	update_ddtUrelHist(ddtUrel,index); // add current dim.less ddtUrel to history
+            	update_rHist(r,index); // add current r to history
+
+            	// warning and reset for too small t0
+            	if (t0<t0min)
+            	{
+            		Pout << "ParmarBassetForce WARNING: t0 = " << t0 << " at ID = " << index <<  endl;
+                	gH0_[index][0]            = NOTONCPU; //reset reference
+                	ddtUrelHist_[0][index][0] = NOTONCPU; //reset ddtU history (only component used for checking nKnown)
+            	}
+
+            	// check length of known history
+            	int nKnown = 0;
+            	for (int j=0; j<nHist_; j++) // loop over past times
+            	{
+            		if (ddtUrelHist_[j][index][0] == NOTONCPU)
+            			break;
+            		else
+            			nKnown++;
+            	}
+
+            	//********* short term force computing (full integral) *********//
+
+            	// number of values to be used for short term force computing
+            	int nShort = min(nKnown,nInt_+1);
+
+            	// int_0^dt K(r,xi) dxi * ddtU(t) dxi (singularity treated by assuming constant acceleration)
+            	if (nShort>0)
+            	{
+            		for (int i=0; i<3; i++) // loop over dimensions
+            			Fshort[i] = -calculateK0(r,dt) * ddtUrelHist_[0][index][i];
+            	}
+
+            	// int_dt^t0 K(r,xi) * ddtU(t-xi) dxi (trapezoid rule)
+            	if (nShort>2)
+            	{
+            		for (int j=1; j<nShort; j++)
+            		{
+            			scalar xi = j*dt;
+            			scalar K  = pow((pow(xi,.25) + rHist_[j][index][0]*xi),-2.); // Eq. 3.4
+
+            			for (int i=0; i<3; i++) // loop over dimensions
+            				Fshort[i] -= trapWeight(j,nShort) * K * ddtUrelHist_[j][index][i] * dt;
+            		}
+            	}
+
+            	//********* long term force computing (differential form) *********//
+
+            	// update F1, F2 history
+            	update_FHist(vector(0., 0., 0.),vector(0., 0., 0.),index);
+
+            	// initialise ddtUrel(t0) and Flong(:) as 0 and r(t0) as 1
+            	if (nKnown == nInt_)
+            	{
+            		for (int j=nInt_; j<nHist_; j++) // loop over past times
+            		{
+            			rHist_[j][index][0] = 1.;
+            			for (int i=0; i<3; i++) // loop over dimensions
+            				ddtUrelHist_[j][index][i] = 0.0;
+            		}
+
+					for (int j=0; j<2*discOrder_+1; j++) // loop over past times
+						for (int i=0; i<3; i++) // loop over dimensions
+							for (int k=0; k<2; k++) // loop over F1, F2				
+								FHist_[k][j][index][i] = 0.0;
+            		nKnown = nHist_;
+            	}
+
+            	// solve ODEs
+            	if (nKnown == nHist_)
+            	{
+            		for (int k=0; k<2; k++) // loop over F1, F2
+            		{
+            			//calculate coefficients
+            			double C[4];
+            			calculateCoeffs(k,t0,rHist_[nInt_][index][0],c,chi,C);
+
+            			// solve Eq. 3.20
+            			solveFlongODE(k,C,dt,index);
+            		}
+            	}
+
+            	//********* total force *********//
+
+            	// sum and convert to N
+            	for (int i=0; i<3; i++) // loop over dimensions
+            	{
+            		ParmarBassetForce[i] = Fshort[i];
+            		for (int k=0; k<2; k++) // loop over F1, F2
+            			ParmarBassetForce[i] += FHist_[k][0][index][i];
+            	}
+            	ParmarBassetForce *= newton;
+
+            	if( forceSubM(0).verbose() ) //&& index>100 && index < 105)
+            	{
+            	}
+
+            	// Set value fields and write the probe
+            	if(probeIt_)
+            	{
+            		scalar ReRef = 0.75/(gH0_[index][0]-0.105);
+
+            		vector Flong1(0,0,0);
+            		vector Flong2(0,0,0);
+
+            		for (int i=0; i<3; i++) // loop over dimensions
+            		{
+                    	Flong1[i] = FHist_[0][0][index][i];
+                    	Flong2[i] = FHist_[1][0][index][i];
+                    }
+
+            		//
+                	// relative velocity (m/s)
+                	vector UrelNoSmooth;
+                	vector ddtUrelNoSmooth;
+
+            		if(forceSubM(0).interpolation())
+                		UrelNoSmooth = UrelNoSmoothInterpolator_.interpolate(position,cellI);
+                	else
+                		UrelNoSmooth = UrelNoSmooth_[cellI];
+
+                	// acceleration (m/s2)
+                	if(forceSubM(0).interpolation())
+                		ddtUrelNoSmooth = ddtUrelNoSmoothInterpolator_.interpolate(position,cellI);
+                	else
+                		ddtUrelNoSmooth = ddtUrelNoSmooth_[cellI];
+
+                	UrelNoSmooth /= mps;
+                	ddtUrelNoSmooth /= mpss;
+            		//
+
+                    #include "setupProbeModelfields.H"
+                    vValues.append(ParmarBassetForce);           //first entry must the be the force
+                    vValues.append(Urel);
+                    vValues.append(ddtUrel);
+                    //
+                    vValues.append(UrelNoSmooth);
+                    vValues.append(ddtUrelNoSmooth);
+                    //
+                    vValues.append(Fshort);
+                    vValues.append(Flong1);
+                    vValues.append(Flong2);
+                    sValues.append(ReRef);
+                    sValues.append(tRef_[index][0]);
+                    sValues.append(mRef_[index][0]);
+                    sValues.append(lRef_[index][0]);
+                    sValues.append(r);
+                    sValues.append(t0);
+                    sValues.append(nKnown);
+                    particleCloud_.probeM().writeProbe(index, sValues, vValues);
+                }
+            }
+            else
+            {
+            	// not on CPU
+            	gH0_[index][0]            = NOTONCPU; //reset reference
+            	ddtUrelHist_[0][index][0] = NOTONCPU; //reset ddtU history (only component used for checking nKnown)
+            }
+
+            // write particle based data to global array
+            forceSubM(0).partToArray(index,ParmarBassetForce,vector::zero);
+    }
+}
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+void Foam::ParmarBassetForce::reAllocArrays() const
+{
+	Pout << "ParmarBassetForce::reAllocArrays..." << endl;
+
+	particleCloud_.dataExchangeM().allocateArray(gH0_,     NOTONCPU,1);
+    particleCloud_.dataExchangeM().allocateArray(tRef_,    NOTONCPU,1);
+    particleCloud_.dataExchangeM().allocateArray(mRef_,    NOTONCPU,1);
+    particleCloud_.dataExchangeM().allocateArray(lRef_,    NOTONCPU,1);
+
+    for (int i=0; i<nHist_; i++)
+    {
+    	particleCloud_.dataExchangeM().allocateArray(ddtUrelHist_[i],NOTONCPU,3);
+    	particleCloud_.dataExchangeM().allocateArray(rHist_      [i],NOTONCPU,1);
+    }
+
+    for (int i=0; i<2*discOrder_+1; i++)
+        for (int k=0; k<2; k++)
+        	particleCloud_.dataExchangeM().allocateArray(FHist_[k][i],0.0,3);
+}
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+scalar Foam::ParmarBassetForce::calculateK0(scalar r, scalar dt) const
+{
+	scalar gamma = pow(r,0.333)*pow(dt,0.25);
+
+	/*
+	scalar K0 = 2./(9.*pow(r,0.666)) *
+			(
+					6.*gamma*gamma/(gamma*gamma*gamma+1)
+					+ 2.*sqrt(3.)*atan(sqrt(3.)*gamma/(2.-gamma))
+					+ log((gamma*gamma-gamma+1.)/(gamma+1.)/(gamma+1.))
+			); // int_0^dt K(r,xi) dxi
+	*/
+	scalar K0 = 2./(9.*pow(r,0.666)) *
+			(
+			  -   .37224 * gamma
+			  + 12.1587  * gamma*gamma
+			  -  6.4884  * gamma*gamma*gamma
+			); // third order approximation to the eq. above
+
+	return K0;
+}
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+scalar Foam::ParmarBassetForce::trapWeight(int i, int n) const
+{
+	if ( (i==1) | (i==(n-1)) )
+		return 0.5;
+	else
+		return 1.0;
+}
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+void Foam::ParmarBassetForce::update_ddtUrelHist(vector ddtUrel, int index) const
+{
+	for (int i=0; i<3; i++) // loop over dimensions
+	{
+		for (int j=nHist_-1; j>0; j--) // loop over past times
+			ddtUrelHist_[j][index][i] = ddtUrelHist_[j-1][index][i];
+
+		ddtUrelHist_[0][index][i] = ddtUrel[i];
+	}
+}
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+void Foam::ParmarBassetForce::update_rHist(scalar r, int index) const
+{
+	for (int j=nHist_-1; j>0; j--) // loop over past times
+		rHist_[j][index][0] = rHist_[j-1][index][0];
+
+	rHist_[0][index][0] = r;
+}
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+void Foam::ParmarBassetForce::update_FHist(vector F1, vector F2, int index) const
+{
+	for (int i=0; i<3; i++) // loop over dimensions
+	{
+		for (int k=0; k<2; k++) // loop over F1, F2
+			for (int j=2*discOrder_; j>0; j--) // loop over past times
+				FHist_[k][j][index][i] = FHist_[k][j-1][index][i];
+
+		FHist_[0][0][index][i] = F1[i];
+		FHist_[1][0][index][i] = F2[i];
+	}
+}
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+void Foam::ParmarBassetForce::calculateCoeffs(int k, scalar t0, scalar r, double c[2][4][3], double chi[2][4][2], double C[4]) const
+{
+	for (int j=0; j<4; j++) // loop over terms
+	{
+		C[j]  = (c[k][j][0] * pow((t0 + c[k][j][1]),c[k][j][2]));    //t0 dependency
+		C[j] *=	(1 + chi[k][j][0]*(r-1) + chi[k][j][1]*(r-1)*(r-1)); // r dependency
+	}
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+void Foam::ParmarBassetForce::solveFlongODE(int k, double C[4], scalar dt, int index) const
+{
+	if (discOrder_==1)
+	{
+		for (int i=0; i<3; i++) // loop over dimensions
+		{
+			FHist_[k][0][index][i] =
+					(
+						(     C[1]/dt+2/(dt*dt)) * FHist_[k][1][index][i]
+					  - (             1/(dt*dt)) * FHist_[k][2][index][i]
+		              - (C[2]+C[3]/dt        ) * ddtUrelHist_[nInt_  ][index][i]
+					  + (     C[3]/dt        ) * ddtUrelHist_[nInt_+1][index][i]
+					) / (C[0]+C[1]/dt+1/(dt*dt)); // Eq. 3.20 using first order temporal discretisation
+		}
+	}
+	if (discOrder_==2)
+	{
+		for (int i=0; i<3; i++) // loop over dimensions
+		{
+			FHist_[k][0][index][i] =
+					(
+						(       4*C[1]/(2*dt) + 24/(4*dt*dt)) * FHist_[k][1][index][i]
+					  - (         C[1]/(2*dt) + 22/(4*dt*dt)) * FHist_[k][2][index][i]
+					  + (                        8/(4*dt*dt)) * FHist_[k][3][index][i]
+					  - (                        1/(4*dt*dt)) * FHist_[k][4][index][i]
+
+		              - (C[2] + 3*C[3]/(2*dt)               ) * ddtUrelHist_[nInt_  ][index][i]
+					  + (       4*C[3]/(2*dt)               ) * ddtUrelHist_[nInt_+1][index][i]
+					  - (         C[3]/(2*dt)               ) * ddtUrelHist_[nInt_+2][index][i]
+					) / (C[0] + 3*C[1]/(2*dt) +  9/(4*dt*dt)); // Eq. 3.20 using second order temporal discretisation
+		}
+	}
+}
+
+void Foam::ParmarBassetForce::rescaleHist(scalar tScale, scalar mScale, scalar lScale, scalar rScale, int index) const
+{
+	for (int i=0; i<3; i++) // loop over dimensions
+    {
+    	// rescale ddtU history
+    	for (int j=0; j<nHist_; j++) // loop over past times
+    		if (ddtUrelHist_[j][index][i] != NOTONCPU)
+    			ddtUrelHist_[j][index][i] *= lScale/(tScale*tScale);
+
+    	// rescale F1, F2 history
+    	for (int k=0; k<2; k++) // loop over F1, F2
+    		for (int j=0; j<(2*discOrder_+1); j++) // loop over past times
+    			FHist_[k][j][index][i] *= mScale*lScale/(tScale*tScale);
+    }
+    // rescale r history
+    for (int j=0; j<nHist_; j++) // loop over past times
+    	if (rHist_[j][index][0] != NOTONCPU)
+    		rHist_[j][index][0] /= pow(rScale,1.5);
+}
+
+} // End namespace Foam
+
+// ************************************************************************* //
+

src/lagrangian/cfdemParticle/subModels/forceModel/ParmarBassetForce/ParmarBassetForce.H

@@ -0,0 +1,155 @@
+/*---------------------------------------------------------------------------*\
+    CFDEMcoupling - Open Source CFD-DEM coupling
+
+    CFDEMcoupling is part of the CFDEMproject
+    www.cfdem.com
+                                Christoph Goniva, christoph.goniva@cfdem.com
+                                Copyright 2009-2012 JKU Linz
+                                Copyright 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, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+Description
+    This code is designed to realize coupled CFD-DEM simulations using LIGGGHTS
+    and OpenFOAM(R). Note: this code is not part of OpenFOAM(R) (see DISCLAIMER).
+
+Class
+    ParmarBassetForce
+
+SourceFiles
+    ParmarBassetForce.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef ParmarBassetForce_H
+#define ParmarBassetForce_H
+
+#include "forceModel.H"
+#include "dataExchangeModel.H"
+#include "interpolationCellPoint.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                           Class ParmarBassetForce Declaration
+\*---------------------------------------------------------------------------*/
+
+class ParmarBassetForce
+:
+    public forceModel
+{
+private:
+    dictionary propsDict_;
+
+    word velFieldName_;
+
+    const volVectorField& U_;
+
+    word UsFieldName_;
+
+    const volVectorField& Us_;
+
+    scalar nInt_;                          //no of timesteps to solve full integral
+
+    scalar discOrder_;                     //ODE discretisation order
+
+    scalar nHist_;                         //no of timesteps to save ddtUrel history for
+
+    mutable List<double**> ddtUrelHist_;
+
+    mutable List<double**> rHist_;
+
+    mutable List<List<double**>> FHist_;
+
+    mutable double** gH0_;
+
+    mutable double** tRef_;
+
+    mutable double** mRef_;
+
+    mutable double** lRef_;
+
+    mutable volVectorField Urel_;
+
+    mutable volVectorField ddtUrel_;
+
+    autoPtr<smoothingModel> smoothingModel_;
+
+
+public:
+
+    //- Runtime type information
+    TypeName("ParmarBassetForce");
+
+
+    // Constructors
+
+        //- Construct from components
+        ParmarBassetForce
+        (
+            const dictionary& dict,
+            cfdemCloud& sm
+        );
+
+    // Destructor
+
+        ~ParmarBassetForce();
+
+
+    // Member Functions
+        void setForce() const;
+
+
+        void reAllocArrays() const;
+
+        inline const smoothingModel& smoothingM() const
+        {
+            return smoothingModel_;
+        }
+
+        scalar calculateK0(scalar r, scalar dt) const;
+
+        scalar trapWeight(int i, int n) const;
+
+        void update_ddtUrelHist(vector ddtUrel, int index) const;
+
+        void update_rHist(scalar r, int index) const;
+
+        void update_FHist(vector F1, vector F2, int index) const;
+
+        void calculateCoeffs(int k, scalar t0, scalar r, double c[2][4][3], double chi[2][4][2], double C[4]) const;
+
+        void solveFlongODE(int k, double C[4], scalar dt, int index) const;
+
+        void rescaleHist(scalar tScale, scalar mScale, scalar lScale, scalar rScale, int index) const;
+
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/lagrangian/cfdemParticle/subModels/forceModel/virtualMassForce/virtualMassForce.C

@@ -67,11 +67,37 @@ virtualMassForce::virtualMassForce
     propsDict_(dict.subDict(typeName + "Props")),
     velFieldName_(propsDict_.lookup("velFieldName")),
     U_(sm.mesh().lookupObject<volVectorField> (velFieldName_)),
+    voidfractionFieldName_(propsDict_.lookup("voidfractionFieldName")),
+    voidfraction_(sm.mesh().lookupObject<volScalarField> (voidfractionFieldName_)),
+    UsFieldName_(propsDict_.lookup("granVelFieldName")),
+    Us_(sm.mesh().lookupObject<volVectorField> (UsFieldName_)),
     phiFieldName_(propsDict_.lookup("phiFieldName")),
     phi_(sm.mesh().lookupObject<surfaceScalarField> (phiFieldName_)),
     UrelOld_(NULL),
     splitUrelCalculation_(propsDict_.lookupOrDefault<bool>("splitUrelCalculation",false)),
-    Cadd_(0.5)
+    useUs_(propsDict_.lookupOrDefault<bool>("useUs",false)),
+    useFelderhof_(propsDict_.lookupOrDefault<bool>("useFelderhof",false)),
+    Cadd_(0.5),
+    DDtUrel_
+    (   IOobject
+        (
+            "DDtUrel",
+            sm.mesh().time().timeName(),
+            sm.mesh(),
+            IOobject::NO_READ,
+            IOobject::AUTO_WRITE
+        ),
+        sm.mesh(),
+        dimensionedVector("zero", dimensionSet(0,1,-2,0,0,0,0), vector(0,0,0))
+    ),
+    smoothingModel_
+    (
+        smoothingModel::New
+        (
+            propsDict_,
+            sm
+        )
+    )
 {
 
     if (particleCloud_.dataExchangeM().maxNumberOfParticles() > 0)
@@ -98,17 +124,39 @@ virtualMassForce::virtualMassForce
         Cadd_ = readScalar(propsDict_.lookup("Cadd"));
         Info << "Virtual mass model: using non-standard Cadd = " << Cadd_ << endl;
     }
+    if(useUs_)
+    {
+        Info << "Virtual mass model: using averaged Us \n";
+        Info << "WARNING: ignoring virtual mass of relative particle motion/collisions \n";
+        
+        if(splitUrelCalculation_)
+        {
+            FatalError << "Virtual mass model: useUs=true requires splitUrelCalculation_=false" << abort(FatalError);
+        }
+    }
+    if(useFelderhof_)
+    {
+        Info << "Virtual mass model: using Cadd correlation by Felderhof \n";
+        Info << "WARNING: ignoring user-set Cadd \n";
+
+	    if(!dict.lookupOrDefault<bool>("getParticleDensities",false))
+	    	FatalError << "Virtual mass model: useFelderhof=true requires getParticleDensities=true in couplingProperties" << abort(FatalError);
+    }
 
     particleCloud_.checkCG(true);
 
     //Append the field names to be probed
     particleCloud_.probeM().initialize(typeName, typeName+".logDat");
     particleCloud_.probeM().vectorFields_.append("virtualMassForce"); //first entry must the be the force
-    particleCloud_.probeM().vectorFields_.append("Urel");
-    particleCloud_.probeM().vectorFields_.append("UrelOld");
-    particleCloud_.probeM().vectorFields_.append("ddtUrel");
+    particleCloud_.probeM().vectorFields_.append("acceleration");
+    //
+    particleCloud_.probeM().vectorFields_.append("accelerationNoSmooth");
+    //
+    particleCloud_.probeM().scalarFields_.append("Cadd");
     particleCloud_.probeM().scalarFields_.append("Vs");
     particleCloud_.probeM().scalarFields_.append("rho");
+    particleCloud_.probeM().scalarFields_.append("voidfraction");
+    particleCloud_.probeM().scalarFields_.append("specificGravity");
     particleCloud_.probeM().writeHeader();
 }
 
@@ -125,22 +173,28 @@ virtualMassForce::~virtualMassForce()
 
 void virtualMassForce::setForce() const
 {
+	if (!useUs_ || !splitUrelCalculation_ )
 		reAllocArrays();
 
     scalar dt = U_.mesh().time().deltaT().value();
 
-    vector position(0,0,0);
-    vector Ufluid(0,0,0);
-    vector Ur(0,0,0);
-    vector DDtU(0,0,0);
-
-    //Compute extra vfields in case it is needed
-    volVectorField DDtU_(0.0*U_/U_.mesh().time().deltaT());
+    //acceleration field
     if(splitUrelCalculation_)
-        DDtU_ = fvc::ddt(U_) + fvc::div(phi_, U_); //Total Derivative of fluid velocity
+        DDtUrel_ = fvc::ddt(U_) + fvc::div(phi_, U_); //Total Derivative of fluid velocity
+    else if(useUs_)
+        DDtUrel_ = fvc::ddt(U_) + fvc::div(phi_, U_) - fvc::ddt(Us_); //Total Derivative of fluid velocity minus average particle velocity
 
-    interpolationCellPoint<vector> UInterpolator_(   U_);
-    interpolationCellPoint<vector> DDtUInterpolator_(DDtU_);
+    //
+        volVectorField DDtUrelNoSmooth_ = DDtUrel_;
+        interpolationCellPoint<vector> DDtUrelNoSmoothInterpolator_(DDtUrelNoSmooth_);
+    //
+    
+    // smoothen
+    smoothingM().smoothen(DDtUrel_);
+
+    interpolationCellPoint<vector> UInterpolator_(U_);
+    interpolationCellPoint<vector> DDtUrelInterpolator_(DDtUrel_);
+    interpolationCellPoint<scalar> voidfractionInterpolator_(voidfraction_);
 
     #include "setupProbeModel.H"
 
@@ -149,37 +203,42 @@ void virtualMassForce::setForce() const
     for(int index = 0;index <  particleCloud_.numberOfParticles(); index++)
     {
             vector virtualMassForce(0,0,0);
+            vector position(0,0,0);
+            vector DDtUrel(0,0,0);
+
+            scalar voidfraction(1);
+            scalar sg(1);
+
             label cellI = particleCloud_.cellIDs()[index][0];
 
             if (cellI > -1) // particle Found
             {
-
+                // particle position
                 if(forceSubM(0).interpolation())
-                {
                     position = particleCloud_.position(index);
-                    Ufluid      = UInterpolator_.interpolate(position,cellI);
-                }
-                else
+
+                //********* acceleration value *********//
+                if(splitUrelCalculation_ || useUs_ )
                 {
-                    Ufluid = U_[cellI];
-                }
-
-
-                if(splitUrelCalculation_)  //if split, just use total derivative of fluid velocity
+                    // DDtUrel from acceleration field 
                     if(forceSubM(0).interpolation())
-                {
-                    DDtU = DDtUInterpolator_.interpolate(position,cellI);
+                		DDtUrel = DDtUrelInterpolator_.interpolate(position,cellI);
+                	else
+                		DDtUrel = DDtUrel_[cellI];
                 }
                 else
                 {
-                    DDtU = DDtU_[cellI];
-                }
-                else
-                {
-                    vector Us = particleCloud_.velocity(index);
-                    Ur = Ufluid - Us;
-                }
+                    // DDtUrel from UrelOld
+                    vector Ufluid(0,0,0);
 
+                    // relative velocity
+                    if(forceSubM(0).interpolation())
+                        Ufluid      = UInterpolator_.interpolate(position,cellI);
+                    else
+                        Ufluid = U_[cellI];
+
+                    vector Us   = particleCloud_.velocity(index);
+                    vector Urel = Ufluid - Us;
 
                     //Check of particle was on this CPU the last step
                 	if(UrelOld_[index][0]==NOTONCPU) //use 1. element to indicate that particle was on this CPU the last time step
@@ -187,24 +246,51 @@ void virtualMassForce::setForce() const
                 	else
                 		haveUrelOld_ = true;
 
-                vector UrelOld(0.,0.,0.);
-                vector ddtUrel(0.,0.,0.);
+                	vector UrelOld(0,0,0);
+
                 	for(int j=0;j<3;j++)
                 	{
                 		UrelOld[j]         = UrelOld_[index][j];
-                    UrelOld_[index][j] = Ur[j];
+                		UrelOld_[index][j] = Urel[j];
                 	}
                 	if(haveUrelOld_ ) //only compute force if we have old (relative) velocity
-                    ddtUrel = (Ur-UrelOld)/dt;
-
-                if(splitUrelCalculation_) //we can always compute the total derivative in case we split
-                    ddtUrel = DDtU;
+                		DDtUrel = (Urel-UrelOld)/dt;
+                }
 
+                //********* Cadd value *********//
+                scalar rho  = forceSubM(0).rhoField()[cellI];
                 scalar ds = 2*particleCloud_.radius(index);
                 scalar Vs = ds*ds*ds*M_PI/6;
-                scalar rho = forceSubM(0).rhoField()[cellI];
                 
-                virtualMassForce = Cadd_ * rho * Vs * ddtUrel;
+                scalar Cadd;
+                
+                if (useFelderhof_)
+                {                    
+                    scalar epsilons(0);
+                    scalar logsg(0);
+
+                	if(forceSubM(0).interpolation())
+		                voidfraction = voidfractionInterpolator_.interpolate(position,cellI);
+		            else
+		                voidfraction = voidfraction_[cellI];
+
+                    sg       = particleCloud_.particleDensity(index) / rho;
+                    logsg    = log(sg);
+                    epsilons = 1-voidfraction;
+
+                    Cadd = 0.5
+                    		+ ( 0.047*logsg + 0.13)*epsilons
+							+ (-0.066*logsg - 0.58)*epsilons*epsilons
+							+ (               1.42)*epsilons*epsilons*epsilons;
+                }
+                else
+                {
+                    // use predefined value
+                    Cadd = Cadd_;
+                }
+
+                //********* calculate force *********//
+                virtualMassForce = Cadd_ * rho * Vs * DDtUrel;
 
                 if( forceSubM(0).verbose() ) //&& index>100 && index < 105)
                 {
@@ -212,21 +298,33 @@ void virtualMassForce::setForce() const
                     Pout << "position = " << particleCloud_.position(index) << endl;
                 }
 
-                //Set value fields and write the probe
+                // Set value fields and write the probe
                 if(probeIt_)
                 {
+                    //
+                    vector DDtUrelNoSmooth(0,0,0);
+                    if(forceSubM(0).interpolation())
+                		DDtUrelNoSmooth = DDtUrelInterpolator_.interpolate(position,cellI);
+                	else
+                		DDtUrelNoSmooth = DDtUrelNoSmooth_[cellI];
+                    //
                     #include "setupProbeModelfields.H"
                     vValues.append(virtualMassForce);           //first entry must the be the force
-                    vValues.append(Ur);
-                    vValues.append(UrelOld);
-                    vValues.append(ddtUrel);
+                    vValues.append(DDtUrel);
+                    //
+                    vValues.append(DDtUrelNoSmooth);
+                    //
+                    sValues.append(Cadd);
                     sValues.append(Vs);
                     sValues.append(rho);
+                    sValues.append(voidfraction);
+                    sValues.append(sg);
                     particleCloud_.probeM().writeProbe(index, sValues, vValues);
                 }
             }
             else    //particle not on this CPU
-                UrelOld_[index][0]=NOTONCPU;
+                if (!useUs_ || !splitUrelCalculation_ )
+                    UrelOld_[index][0] = NOTONCPU;
 
             // write particle based data to global array
             forceSubM(0).partToArray(index,virtualMassForce,vector::zero);

src/lagrangian/cfdemParticle/subModels/forceModel/virtualMassForce/virtualMassForce.H

@@ -67,6 +67,14 @@ private:
 
     const volVectorField& U_;
 
+    word voidfractionFieldName_;
+
+    const volScalarField& voidfraction_;
+
+    word UsFieldName_;
+
+    const volVectorField& Us_;
+
     word phiFieldName_;
 
     const surfaceScalarField& phi_;
@@ -76,7 +84,15 @@ private:
     const bool splitUrelCalculation_;      //indicator to split calculation of Urel between CFDEM and LIGGGHTS
                                            //requires the integration fix to take dv/dt into account!
 
-    mutable double Cadd_;                  //virtual mass constant
+    mutable bool useUs_;				   //indicator to use averaged Us instead of v_p
+
+    mutable bool useFelderhof_;			   //indicator to use Cadd values from Felderhof
+
+    mutable double Cadd_;                  //custom Cadd value
+
+    mutable volVectorField DDtUrel_;
+
+    autoPtr<smoothingModel> smoothingModel_;
 
 public:
 
@@ -102,6 +118,11 @@ public:
         void setForce() const;
 
         void reAllocArrays() const;
+
+        inline const smoothingModel& smoothingM() const
+        {
+            return smoothingModel_;
+        }
 };
 
 

src/lagrangian/cfdemParticle/subModels/liggghtsCommandModel/liggghtsCommandModel/liggghtsCommandModel.C

@@ -147,7 +147,7 @@ void liggghtsCommandModel::checkTimeSettings(const dictionary& propsDict)
                 // if this makes troubles try floor((startTime_+SMALL)/.. as above
                 firstCouplingStep_ = floor((startTime_+SMALL-simStartTime)/DEMts/couplingInterval);
                 lastCouplingStep_ = floor((endTime_+SMALL-simStartTime)/DEMts/couplingInterval);
-                couplingStepInterval_ = floor(timeInterval_+SMALL/DEMts/couplingInterval);
+                couplingStepInterval_ = floor((timeInterval_+SMALL)/DEMts/couplingInterval);
             }
             else      //runEveryCouplingStep  or writeStep
             {
@@ -279,6 +279,11 @@ void liggghtsCommandModel::parseCommandList(wordList& commandList,labelList& lab
         else if (add=="dot")    add = ".";
         else if (add=="dotdot") add = "..";
         else if (add=="slash")  add = "/";
+	else if (add=="dollar") add = "$";
+	else if (add=="curlyOpen") add = "{";
+	else if (add=="curlyClose") add = "}";
+	else if (add=="squareOpen") add = "[";
+	else if (add=="squareClose") add = "]";
         else if (add=="noBlanks")  // no blanks after the following words
         {
             add = "";

src/lagrangian/cfdemParticle/subModels/smoothingModel/constDiffAndTemporalSmoothing/constDiffAndTemporalSmoothing.C

@@ -0,0 +1,265 @@
+/*---------------------------------------------------------------------------*\
+    CFDEMcoupling - Open Source CFD-DEM coupling
+
+    CFDEMcoupling is part of the CFDEMproject
+    www.cfdem.com
+                                Christoph Goniva, christoph.goniva@cfdem.com
+                                Copyright 2009-2012 JKU Linz
+                                Copyright 2012-     DCS Computing GmbH, Linz
+                                Copyright (C) 2013-     Graz University of
+                                                        Technology, IPPT
+-------------------------------------------------------------------------------
+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, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+Description
+    This code is designed to realize coupled CFD-DEM simulations using LIGGGHTS
+    and OpenFOAM(R). Note: this code is not part of OpenFOAM(R) (see DISCLAIMER).
+\*---------------------------------------------------------------------------*/
+
+#include "error.H"
+
+#include "constDiffAndTemporalSmoothing.H"
+#include "addToRunTimeSelectionTable.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+defineTypeNameAndDebug(constDiffAndTemporalSmoothing, 0);
+
+addToRunTimeSelectionTable
+(
+    smoothingModel,
+    constDiffAndTemporalSmoothing,
+    dictionary
+);
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+// Construct from components
+constDiffAndTemporalSmoothing::constDiffAndTemporalSmoothing
+(
+    const dictionary& dict,
+    cfdemCloud& sm
+)
+:
+    smoothingModel(dict,sm),
+    propsDict_(dict.subDict(typeName + "Props")),
+    lowerLimit_(readScalar(propsDict_.lookup("lowerLimit"))),
+    upperLimit_(readScalar(propsDict_.lookup("upperLimit"))),
+    smoothingLength_(dimensionedScalar("smoothingLength",dimensionSet(0,1,0,0,0,0,0), readScalar(propsDict_.lookup("smoothingLength")))),
+    smoothingLengthReferenceField_(dimensionedScalar("smoothingLengthReferenceField",dimensionSet(0,1,0,0,0,0,0), readScalar(propsDict_.lookup("smoothingLength")))),
+    DT_("DT", dimensionSet(0,2,-1,0,0), 0.),
+	gamma_(readScalar(propsDict_.lookup("smoothingStrength"))),
+	correctBoundary_(propsDict_.lookupOrDefault<bool>("correctBoundary",false)),
+    verbose_(false)
+{
+
+    if(propsDict_.found("verbose"))
+        verbose_ = true;
+
+    if(propsDict_.found("smoothingLengthReferenceField"))
+       smoothingLengthReferenceField_.value() = double(readScalar(propsDict_.lookup("smoothingLengthReferenceField")));
+
+    checkFields(sSmoothField_);
+    checkFields(vSmoothField_);
+}
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+constDiffAndTemporalSmoothing::~constDiffAndTemporalSmoothing()
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+bool constDiffAndTemporalSmoothing::doSmoothing() const
+{
+    return true;
+}
+
+
+void constDiffAndTemporalSmoothing::smoothen(volScalarField& fieldSrc) const
+{
+    // Create scalar smooth field from virgin scalar smooth field template
+    volScalarField sSmoothField = sSmoothField_;
+
+    sSmoothField.dimensions().reset(fieldSrc.dimensions());
+    sSmoothField.ref()=fieldSrc.internalField();
+    sSmoothField.correctBoundaryConditions();
+    sSmoothField.oldTime().dimensions().reset(fieldSrc.dimensions());
+    sSmoothField.oldTime()=fieldSrc;
+    sSmoothField.oldTime().correctBoundaryConditions();
+
+    double deltaT = sSmoothField.mesh().time().deltaTValue();
+    DT_.value() = smoothingLength_.value() * smoothingLength_.value() / deltaT;
+
+    // do smoothing
+    solve
+    (
+        fvm::ddt(sSmoothField)
+       -fvm::laplacian(DT_, sSmoothField)
+    );
+
+    // bound sSmoothField_
+    forAll(sSmoothField,cellI)
+    {
+        sSmoothField[cellI]=max(lowerLimit_,min(upperLimit_,sSmoothField[cellI]));
+    }
+
+    // get data from working sSmoothField - will copy only values at new time
+    fieldSrc=sSmoothField;
+    fieldSrc.correctBoundaryConditions();
+
+    if(verbose_)
+    {
+        Info << "min/max(fieldoldTime) (unsmoothed): " << min(sSmoothField.oldTime()) << tab << max(sSmoothField.oldTime()) << endl;
+        Info << "min/max(fieldSrc): " << min(fieldSrc) << tab << max(fieldSrc) << endl;
+        Info << "min/max(fieldSrc.oldTime): " << min(fieldSrc.oldTime()) << tab << max(fieldSrc.oldTime()) << endl;
+    }
+
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+void constDiffAndTemporalSmoothing::smoothen(volVectorField& fieldSrc) const
+{
+    // Create scalar smooth field from virgin scalar smooth field template
+    volVectorField vSmoothField = vSmoothField_;
+
+    vSmoothField.dimensions().reset(fieldSrc.dimensions());
+    vSmoothField.ref()=fieldSrc.internalField();
+    vSmoothField.correctBoundaryConditions();
+    vSmoothField.oldTime().dimensions().reset(fieldSrc.dimensions());
+    vSmoothField.oldTime()=fieldSrc;
+    vSmoothField.oldTime().correctBoundaryConditions();
+
+    dimensionedScalar deltaT = vSmoothField.mesh().time().deltaT();
+    DT_.value() = smoothingLength_.value() * smoothingLength_.value() / deltaT.value();
+
+    // do spacial smoothing
+    solve
+    (
+        fvm::ddt(vSmoothField)
+       -fvm::laplacian(DT_, vSmoothField)
+    );
+
+    // create fields for temporal smoothing
+    volVectorField refField = vSmoothField;
+
+    vSmoothField.ref()=fieldSrc.oldTime().internalField();
+    vSmoothField.correctBoundaryConditions();
+    vSmoothField.oldTime()=fieldSrc.oldTime();
+    vSmoothField.oldTime().correctBoundaryConditions();
+
+    // do temporal smoothing
+    solve
+    (
+        fvm::ddt(vSmoothField)
+        -
+        gamma_/deltaT * (refField - fvm::Sp(1.0,vSmoothField))
+    );
+
+    // create temporally smoothened boundary field
+    if (correctBoundary_)
+    	vSmoothField.boundaryFieldRef() = 1/(1+gamma_)*fieldSrc.oldTime().boundaryField() + gamma_/(1+gamma_)*fieldSrc.boundaryField();
+
+    // get data from working vSmoothField
+    fieldSrc=vSmoothField;
+    fieldSrc.correctBoundaryConditions();
+
+    if(verbose_)
+    {
+        Info << "min/max(fieldoldTime) (unsmoothed): " << min(vSmoothField.oldTime()) << tab << max(vSmoothField.oldTime()) << endl;
+        Info << "min/max(fieldSrc): " << min(fieldSrc) << tab << max(fieldSrc) << endl;
+        Info << "min/max(fieldSrc.oldTime): " << min(fieldSrc.oldTime()) << tab << max(fieldSrc.oldTime()) << endl;
+    }
+}
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+void constDiffAndTemporalSmoothing::smoothenReferenceField(volVectorField& fieldSrc) const
+{
+    // Create scalar smooth field from virgin scalar smooth field template
+    volVectorField vSmoothField = vSmoothField_;
+
+    vSmoothField.dimensions().reset(fieldSrc.dimensions());
+    vSmoothField.ref()=fieldSrc.internalField();
+    vSmoothField.correctBoundaryConditions();
+    vSmoothField.oldTime().dimensions().reset(fieldSrc.dimensions());
+    vSmoothField.oldTime()=fieldSrc;
+    vSmoothField.oldTime().correctBoundaryConditions();
+
+    double sourceStrength = 1e5; //large number to keep reference values constant
+
+    dimensionedScalar deltaT = vSmoothField.mesh().time().deltaT();
+    DT_.value() = smoothingLengthReferenceField_.value()
+                         * smoothingLengthReferenceField_.value() / deltaT.value();
+
+    tmp<volScalarField> NLarge
+    (
+        new volScalarField
+        (
+            IOobject
+            (
+                "NLarge",
+                particleCloud_.mesh().time().timeName(),
+                particleCloud_.mesh(),
+                IOobject::NO_READ,
+                IOobject::NO_WRITE
+            ),
+            particleCloud_.mesh(),
+            0.0
+        )
+    );
+
+
+    //loop over particles and map max particle diameter to Euler Grid
+    forAll(vSmoothField,cellI)
+    {
+        if ( mag(vSmoothField.oldTime().internalField()[cellI]) > 0.0f)  // have a vector in the OLD vSmoothField, so keep it!
+            NLarge.ref()[cellI] = sourceStrength;
+    }
+
+    // do the smoothing
+    solve
+    (
+        fvm::ddt(vSmoothField)
+       -fvm::laplacian( DT_, vSmoothField)
+       ==
+        NLarge() / deltaT * vSmoothField.oldTime()  //add source to keep cell values constant
+       -fvm::Sp( NLarge() / deltaT, vSmoothField)   //add sink to keep cell values constant
+    );
+
+    // get data from working vSmoothField
+    fieldSrc=vSmoothField;
+    fieldSrc.correctBoundaryConditions();
+
+    if(verbose_)
+    {
+        Info << "min/max(fieldoldTime) (unsmoothed): " << min(vSmoothField.oldTime()) << tab << max(vSmoothField.oldTime()) << endl;
+        Info << "min/max(fieldSrc): " << min(fieldSrc) << tab << max(fieldSrc) << endl;
+        Info << "min/max(fieldSrc.oldTime): " << min(fieldSrc.oldTime()) << tab << max(fieldSrc.oldTime()) << endl;
+    }
+
+}
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// ************************************************************************* //

src/lagrangian/cfdemParticle/subModels/smoothingModel/constDiffAndTemporalSmoothing/constDiffAndTemporalSmoothing.H

@@ -0,0 +1,114 @@
+/*---------------------------------------------------------------------------*\
+    CFDEMcoupling - Open Source CFD-DEM coupling
+
+    CFDEMcoupling is part of the CFDEMproject
+    www.cfdem.com
+                                Christoph Goniva, christoph.goniva@cfdem.com
+                                Copyright 2009-2012 JKU Linz
+                                Copyright 2012-     DCS Computing GmbH, Linz
+                                Copyright (C) 2013-     Graz University of  
+                                                        Technology, IPPT
+-------------------------------------------------------------------------------
+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, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+Description
+    This code is designed to realize coupled CFD-DEM simulations using LIGGGHTS
+    and OpenFOAM(R). Note: this code is not part of OpenFOAM(R) (see DISCLAIMER).
+
+Class
+    constDiffAndTemporalSmoothing
+
+SourceFiles
+    constDiffAndTemporalSmoothing.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef constDiffAndTemporalSmoothing_H
+#define constDiffAndTemporalSmoothing_H
+
+#include "smoothingModel.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                           Class constDiffAndTemporalSmoothing Declaration
+\*---------------------------------------------------------------------------*/
+
+class constDiffAndTemporalSmoothing
+:
+    public smoothingModel
+{
+
+private:
+
+    dictionary propsDict_;
+    scalar lowerLimit_;
+    scalar upperLimit_;
+    dimensionedScalar smoothingLength_;
+    dimensionedScalar smoothingLengthReferenceField_;
+    mutable dimensionedScalar DT_;
+    scalar gamma_;
+    bool correctBoundary_;
+    bool verbose_;
+
+public:
+
+    //- Runtime type information
+    TypeName("constDiffAndTemporalSmoothing");
+
+
+    // Constructors
+
+        //- Construct from components
+        constDiffAndTemporalSmoothing
+        (
+            const dictionary& dict,
+            cfdemCloud& sm
+        );
+
+    // Destructor
+
+        ~constDiffAndTemporalSmoothing();
+
+
+    // Member Functions
+        bool doSmoothing() const;
+
+        //void dSmoothing(volScalarField&) const;
+        
+        void smoothen(volScalarField&) const;
+
+        void smoothen(volVectorField&) const;
+
+        void smoothenReferenceField(volVectorField&) const;
+
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/lagrangian/cfdemParticle/subModels/thermCondModel/SyamlalThermCond/SyamlalThermCond.C

@@ -67,7 +67,12 @@ SyamlalThermCond::SyamlalThermCond
         sm.mesh(),
         dimensionedScalar("zero", dimensionSet(0,0,0,0,0,0,0), 1.0)
     ),
-    hasWallQFactor_(false)
+    hasWallQFactor_(false),
+    multiphase_(propsDict_.lookupOrDefault<bool>("multiphase",false)),
+    kfFieldName_(propsDict_.lookupOrDefault<word>("kfFieldName",voidfractionFieldName_)), // use voidfractionField as dummy to prevent lookup error when not using multiphase
+    kfField_(sm.mesh().lookupObject<volScalarField> (kfFieldName_)),
+    CpFieldName_(propsDict_.lookupOrDefault<word>("CpFieldName",voidfractionFieldName_)), // use voidfractionField as dummy to prevent lookup error when not using multiphase
+    CpField_(sm.mesh().lookupObject<volScalarField> (CpFieldName_))
 {
     if (wallQFactor_.headerOk())
     {
@@ -106,12 +111,17 @@ tmp<volScalarField> SyamlalThermCond::thermCond() const
 
     volScalarField& svf = tvf.ref();
 
-
     forAll(svf,cellI)
     {
-        if (1-voidfraction_[cellI] < SMALL) svf[cellI] = kf0_.value();
+        scalar kf;
+        if (multiphase_)
+	    kf = kfField_[cellI];
+        else
+	    kf = kf0_.value();
+      
+        if (1-voidfraction_[cellI] < SMALL) svf[cellI] = kf;
         else if (voidfraction_[cellI] < SMALL) svf[cellI] = 0.0;
-        else svf[cellI] = (1-sqrt(1-voidfraction_[cellI]+SMALL)) / (voidfraction_[cellI]) * kf0_.value();
+        else svf[cellI] = (1-sqrt(1-voidfraction_[cellI]+SMALL)) / (voidfraction_[cellI]) * kf;
     }
 
     // if a wallQFactor field is present, use it to scale heat transport through a patch
@@ -127,6 +137,9 @@ tmp<volScalarField> SyamlalThermCond::thermCond() const
 
 tmp<volScalarField> SyamlalThermCond::thermDiff() const
 {
+	if (multiphase_)
+		return thermCond()/(rho_*CpField_);
+	else
 		return thermCond()/(rho_*Cp_);
 }
 

src/lagrangian/cfdemParticle/subModels/thermCondModel/SyamlalThermCond/SyamlalThermCond.H

@@ -68,6 +68,16 @@ private:
 
     bool hasWallQFactor_;
 
+    bool multiphase_;
+
+    word kfFieldName_;
+
+    const volScalarField& kfField_;
+
+    word CpFieldName_;
+
+    const volScalarField& CpField_;
+
 public:
 
     //- Runtime type information

src/lagrangian/cfdemParticleComp/Make/files

@@ -40,6 +40,8 @@ $(energyModels)/energyModel/energyModel.C
 $(energyModels)/energyModel/newEnergyModel.C
 $(energyModels)/heatTransferGunn/heatTransferGunn.C
 $(energyModels)/heatTransferGunnPartField/heatTransferGunnPartField.C
+$(energyModels)/wallHeatTransferYagi/wallHeatTransferYagi.C
+$(energyModels)/wallHeatTransferYagiImplicit/wallHeatTransferYagiImplicit.C
 $(energyModels)/reactionHeat/reactionHeat.C
 
 $(thermCondModels)/thermCondModel/thermCondModel.C
@@ -63,6 +65,7 @@ $(forceModels)/ShirgaonkarIB/ShirgaonkarIB.C
 $(forceModels)/KochHillDrag/KochHillDrag.C
 $(forceModels)/LaEuScalarTemp/LaEuScalarTemp.C
 $(forceModels)/virtualMassForce/virtualMassForce.C
+$(forceModels)/ParmarBassetForce/ParmarBassetForce.C
 $(forceModels)/gradPForce/gradPForce.C
 $(forceModels)/viscForce/viscForce.C
 $(forceModels)/MeiLift/MeiLift.C