Merge pull request #56 from ParticulateFlow/feature/cfdemSolverRhoSimple

Feature/cfdem solver rho simple

applications/solvers/cfdemSolverRhoPimple/EEqn.H

@@ -39,7 +39,7 @@
      ==
         fvOptions(rho, he)
     );
-    
+
 
     EEqn.relax();
 
@@ -53,7 +53,8 @@
 
     Info<< "T max/min : " << max(T).value() << " " << min(T).value() << endl;
 
-    particleCloud.clockM().start(31,"postFlow");
+    particleCloud.clockM().start(31,"energySolve");
-    particleCloud.postFlow();
+    particleCloud.solve();
-    particleCloud.clockM().stop("postFlow");
+    particleCloud.clockM().stop("energySolve");
+
 }

applications/solvers/cfdemSolverRhoPimple/cfdemSolverRhoPimple.C

@@ -139,6 +139,10 @@ int main(int argc, char *argv[])
             }
         }
 
+        particleCloud.clockM().start(31,"postFlow");
+        particleCloud.postFlow();
+        particleCloud.clockM().stop("postFlow");
+
         runTime.write();
 
 

applications/solvers/cfdemSolverRhoSimple/EEqn.H

@@ -0,0 +1,57 @@
+// contributions to internal energy equation can be found in
+// Crowe et al.: "Multiphase flows with droplets and particles", CRC Press 1998
+{
+    // dim he = J / kg
+    volScalarField& he = thermo.he();
+    particleCloud.energyContributions(Qsource);
+    particleCloud.energyCoefficients(QCoeff);
+
+    //thDiff=particleCloud.thermCondM().thermDiff();   
+    thCond=particleCloud.thermCondM().thermCond(); 
+
+    addSource = 
+        (
+            he.name() == "e"
+          ? 
+            fvc::div(phi, K) +
+            fvc::div
+            (
+                fvc::absolute(phi/fvc::interpolate(rho), voidfraction*U),
+                p,
+                "div(phiv,p)"
+            )
+          : fvc::div(phi, K)
+        );
+
+    Cpv = he.name() == "e" ? thermo.Cv() : thermo.Cp();
+
+
+    fvScalarMatrix EEqn
+    (
+        fvm::div(phi, he)
+      + addSource
+      - Qsource
+      - fvm::Sp(QCoeff/Cpv, he)
+      - fvm::laplacian(voidfraction*thCond/Cpv,he)
+     ==
+        fvOptions(rho, he)
+    );
+
+
+    EEqn.relax();
+
+    fvOptions.constrain(EEqn);
+
+    EEqn.solve();
+
+    fvOptions.correct(he);
+
+    thermo.correct();
+
+    Info<< "T max/min : " << max(T).value() << " " << min(T).value() << endl;
+
+
+    particleCloud.clockM().start(31,"energySolve");
+    particleCloud.solve();
+    particleCloud.clockM().stop("energySolve");
+}

applications/solvers/cfdemSolverRhoSimple/Make/files

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

applications/solvers/cfdemSolverRhoSimple/Make/options

@@ -0,0 +1,32 @@
+include $(CFDEM_ADD_LIBS_DIR)/additionalLibs
+
+PFLAGS+= -Dcompre
+
+EXE_INC = \
+     $(PFLAGS) \
+    -I$(CFDEM_OFVERSION_DIR) \
+    -I$(LIB_SRC)/transportModels/compressible/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
+    -I$(LIB_SRC)/TurbulenceModels/turbulenceModels/lnInclude \
+    -I$(LIB_SRC)/TurbulenceModels/compressible/lnInclude \
+    -I$(LIB_SRC)/finiteVolume/cfdTools \
+    -I$(LIB_SRC)/finiteVolume/lnInclude \
+    -I$(LIB_SRC)/meshTools/lnInclude \
+    -I$(LIB_SRC)/sampling/lnInclude \
+    -I$(CFDEM_SRC_DIR)/lagrangian/cfdemParticle/lnInclude \
+    -I$(CFDEM_SRC_DIR)/lagrangian/cfdemParticle/cfdTools \
+
+EXE_LIBS = \
+    -L$(CFDEM_LIB_DIR)\
+    -lcompressibleTransportModels \
+    -lfluidThermophysicalModels \
+    -lspecie \
+    -lturbulenceModels \
+    -lcompressibleTurbulenceModels \
+    -lfiniteVolume \
+    -lmeshTools \
+    -lsampling \
+    -lfvOptions \
+    -l$(CFDEM_LIB_COMP_NAME) \
+     $(CFDEM_ADD_LIB_PATHS) \
+     $(CFDEM_ADD_LIBS)

applications/solvers/cfdemSolverRhoSimple/UEqn.H

@@ -0,0 +1,30 @@
+// Solve the Momentum equation
+particleCloud.otherForces(fOther);
+
+tmp<fvVectorMatrix> tUEqn
+(
+    fvm::div(phi, U)
+  + particleCloud.divVoidfractionTau(U, voidfraction)
+  + fvm::Sp(Ksl,U)
+  - fOther
+ ==
+    fvOptions(rho, U)
+);
+fvVectorMatrix& UEqn = tUEqn.ref();
+
+UEqn.relax();
+
+fvOptions.constrain(UEqn);
+
+if (modelType=="B" || modelType=="Bfull")
+{
+    solve(UEqn == -fvc::grad(p)+ Ksl*Us);
+}
+else
+{
+    solve(UEqn == -voidfraction*fvc::grad(p)+ Ksl*Us);
+}
+
+fvOptions.correct(U);
+
+K = 0.5*magSqr(U);

applications/solvers/cfdemSolverRhoSimple/cfdemSolverRhoSimple.C

@@ -0,0 +1,142 @@
+/*---------------------------------------------------------------------------*\
+License
+
+    This is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    This code is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with this code.  If not, see <http://www.gnu.org/licenses/>.
+
+    Copyright (C) 2015- Thomas Lichtenegger, JKU Linz, Austria
+
+Application
+    cfdemSolverRhoPimple
+
+Description
+    Transient solver for compressible flow using the flexible PIMPLE (PISO-SIMPLE)
+    algorithm.
+    Turbulence modelling is generic, i.e. laminar, RAS or LES may be selected.
+    The code is an evolution of the solver rhoPimpleFoam in OpenFOAM(R) 4.x,
+    where additional functionality for CFD-DEM coupling is added.
+\*---------------------------------------------------------------------------*/
+
+#include "fvCFD.H"
+#include "psiThermo.H"
+#include "turbulentFluidThermoModel.H"
+#include "bound.H"
+#include "simpleControl.H"
+#include "fvOptions.H"
+#include "localEulerDdtScheme.H"
+#include "fvcSmooth.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 "createTimeControls.H"
+    #include "createRDeltaT.H"
+    #include "initContinuityErrs.H"
+    #include "createFields.H"
+    #include "createFieldRefs.H"
+    #include "createFvOptions.H"
+
+    // create cfdemCloud
+    #include "readGravitationalAcceleration.H"
+    cfdemCloudEnergy particleCloud(mesh);
+    #include "checkModelType.H"
+
+    turbulence->validate();
+
+    // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+    Info<< "\nStarting time loop\n" << endl;
+
+    while (simple.loop())
+    {
+        particleCloud.clockM().start(1,"Global");
+
+        Info<< "Time = " << runTime.timeName() << nl << endl;
+
+        // do particle stuff
+        particleCloud.clockM().start(2,"Coupling");
+        bool hasEvolved = particleCloud.evolve(voidfraction,Us,U);
+
+        if(hasEvolved)
+        {
+            particleCloud.smoothingM().smoothen(particleCloud.forceM(0).impParticleForces());
+        }
+
+        Info << "update Ksl.internalField()" << endl;
+        Ksl = particleCloud.momCoupleM(0).impMomSource();
+        Ksl.correctBoundaryConditions();
+
+        //Force Checks
+        vector fTotal(0,0,0);
+        vector fImpTotal = sum(mesh.V()*Ksl.primitiveFieldRef()*(Us.primitiveFieldRef()-U.primitiveFieldRef()));
+        reduce(fImpTotal, sumOp<vector>());
+        Info << "TotalForceExp: " << fTotal << endl;
+        Info << "TotalForceImp: " << fImpTotal << endl;
+
+        #include "solverDebugInfo.H"
+        particleCloud.clockM().stop("Coupling");
+
+        particleCloud.clockM().start(26,"Flow");
+
+        volScalarField rhoeps("rhoeps",rho*voidfraction);
+        // Pressure-velocity SIMPLE corrector
+
+        #include "UEqn.H"
+
+
+        // besides this pEqn, OF offers a "simple consistent"-option
+        #include "pEqn.H"
+        rhoeps=rho*voidfraction;
+
+        #include "EEqn.H"
+
+        turbulence->correct();
+
+        particleCloud.clockM().start(32,"postFlow");
+        if(hasEvolved) particleCloud.postFlow();
+        particleCloud.clockM().stop("postFlow");
+
+        runTime.write();
+
+
+        Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
+            << "  ClockTime = " << runTime.elapsedClockTime() << " s"
+            << nl << endl;
+
+        particleCloud.clockM().stop("Flow");
+        particleCloud.clockM().stop("Global");
+    }
+
+    Info<< "End\n" << endl;
+
+    return 0;
+}
+
+
+// ************************************************************************* //

applications/solvers/cfdemSolverRhoSimple/createFieldRefs.H

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

applications/solvers/cfdemSolverRhoSimple/createFields.H

@@ -0,0 +1,242 @@
+Info<< "Reading thermophysical properties\n" << endl;
+
+    autoPtr<psiThermo> pThermo
+    (
+        psiThermo::New(mesh)
+    );
+    psiThermo& thermo = pThermo();
+    thermo.validate(args.executable(), "h", "e");
+    volScalarField& p = thermo.p();
+
+    Info<< "Reading field rho\n" << endl;
+    volScalarField rho
+    (
+        IOobject
+        (
+            "rho",
+            runTime.timeName(),
+            mesh,
+            IOobject::READ_IF_PRESENT,
+            IOobject::AUTO_WRITE
+        ),
+        thermo.rho()
+    );
+
+
+    Info<< "Reading field U\n" << endl;
+    volVectorField U
+    (
+        IOobject
+        (
+            "U",
+            runTime.timeName(),
+            mesh,
+            IOobject::MUST_READ,
+            IOobject::AUTO_WRITE
+        ),
+        mesh
+    );
+
+    Info<< "\nReading voidfraction field voidfraction = (Vgas/Vparticle)\n" << endl;
+    volScalarField voidfraction
+    (
+        IOobject
+        (
+            "voidfraction",
+            runTime.timeName(),
+            mesh,
+            IOobject::MUST_READ,
+            IOobject::AUTO_WRITE
+        ),
+        mesh
+    );
+
+    volScalarField addSource
+    (
+        IOobject
+        (
+           "addSource",
+           runTime.timeName(),
+           mesh,
+           IOobject::MUST_READ,
+           IOobject::AUTO_WRITE
+        ),
+        mesh
+    );
+
+    Info<< "\nCreating fluid-particle heat flux field\n" << endl;
+    volScalarField Qsource
+    (
+        IOobject
+        (
+            "Qsource",
+            runTime.timeName(),
+            mesh,
+            IOobject::NO_READ,
+            IOobject::AUTO_WRITE
+        ),
+        mesh,
+        dimensionedScalar("zero", dimensionSet(1,-1,-3,0,0,0,0), 0.0)
+    );
+
+    Info<< "\nCreating fluid-particle heat flux coefficient field\n" << endl;
+    volScalarField QCoeff
+    (
+        IOobject
+        (
+            "QCoeff",
+            runTime.timeName(),
+            mesh,
+            IOobject::NO_READ,
+            IOobject::AUTO_WRITE
+        ),
+        mesh,
+        dimensionedScalar("zero", dimensionSet(1,-1,-3,-1,0,0,0), 0.0)
+    );
+
+    Info<< "\nCreating thermal conductivity field\n" << endl;
+    volScalarField thCond
+    (
+        IOobject
+        (
+            "thCond",
+            runTime.timeName(),
+            mesh,
+            IOobject::NO_READ,
+            IOobject::AUTO_WRITE
+        ),
+        mesh,
+        dimensionedScalar("zero", dimensionSet(1,1,-3,-1,0,0,0), 0.0)
+    );
+
+    Info<< "\nCreating heat capacity field\n" << endl;
+    volScalarField Cpv
+    (
+        IOobject
+        (
+            "Cpv",
+            runTime.timeName(),
+            mesh,
+            IOobject::NO_READ,
+            IOobject::AUTO_WRITE
+        ),
+        mesh,
+        dimensionedScalar("zero", dimensionSet(0,2,-2,-1,0,0,0), 0.0)
+    );
+
+    Info<< "\nCreating body force field\n" << endl;
+    volVectorField fOther
+    (
+        IOobject
+        (
+            "fOther",
+            runTime.timeName(),
+            mesh,
+            IOobject::NO_READ,
+            IOobject::NO_WRITE
+        ),
+        mesh,
+        dimensionedVector("zero", dimensionSet(1,-2,-2,0,0,0,0), vector::zero)
+    );
+
+    Info<< "Reading/calculating face flux field phi\n" << endl;
+    surfaceScalarField phi
+    (
+        IOobject
+        (
+            "phi",
+            runTime.timeName(),
+            mesh,
+            IOobject::READ_IF_PRESENT,
+            IOobject::AUTO_WRITE
+         ),
+         linearInterpolate(rho*U*voidfraction) & mesh.Sf()
+    );
+
+    dimensionedScalar rhoMax
+    (
+        dimensionedScalar::lookupOrDefault
+        (
+            "rhoMax",
+            simple.dict(),
+            dimDensity,
+            GREAT
+        )
+    );
+
+    dimensionedScalar rhoMin
+    (
+        dimensionedScalar::lookupOrDefault
+        (
+            "rhoMin",
+            simple.dict(),
+            dimDensity,
+            0
+        )
+    );
+
+    Info<< "Creating turbulence model\n" << endl;
+    autoPtr<compressible::turbulenceModel> turbulence
+    (
+        compressible::turbulenceModel::New
+        (
+            rho,
+            U,
+            phi,
+            thermo
+        )
+    );
+
+    label pRefCell = 0;
+    scalar pRefValue = 0.0;
+    setRefCell(p, simple.dict(), pRefCell, pRefValue);
+
+    mesh.setFluxRequired(p.name());
+
+    Info<< "Creating field dpdt\n" << endl;
+    volScalarField dpdt
+    (
+        IOobject
+        (
+            "dpdt",
+            runTime.timeName(),
+            mesh
+        ),
+        mesh,
+        dimensionedScalar("dpdt", p.dimensions()/dimTime, 0)
+    );
+
+    Info<< "Creating field kinetic energy K\n" << endl;
+    volScalarField K("K", 0.5*magSqr(U));
+
+    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<< "Reading particle velocity field Us\n" << endl;
+    volVectorField Us
+    (
+        IOobject
+        (
+            "Us",
+            runTime.timeName(),
+            mesh,
+            IOobject::MUST_READ,
+            IOobject::AUTO_WRITE
+        ),
+        mesh
+    );
+
+//===============================

applications/solvers/cfdemSolverRhoSimple/pEqn.H

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

etc/solver-list.txt

@@ -1,5 +1,6 @@
 cfdemSolverPisoMS/dir
 cfdemSolverPiso/dir
 cfdemSolverRhoPimple/dir
+cfdemSolverRhoSimple/dir
 cfdemSolverIB/dir
 cfdemSolverPisoScalar/dir

src/lagrangian/cfdemParticle/Make/files

@@ -32,12 +32,13 @@ $(cfdTools)/newGlobal.C
 $(energyModels)/energyModel/energyModel.C
 $(energyModels)/energyModel/newEnergyModel.C
 $(energyModels)/heatTransferGunn/heatTransferGunn.C
-$(energyModels)/heatTransferGunnImplicit/heatTransferGunnImplicit.C
+$(energyModels)/heatTransferGunnPartField/heatTransferGunnPartField.C
 $(energyModels)/reactionHeat/reactionHeat.C
 
 $(thermCondModels)/thermCondModel/thermCondModel.C
 $(thermCondModels)/thermCondModel/newThermCondModel.C
 $(thermCondModels)/SyamlalThermCond/SyamlalThermCond.C
+$(thermCondModels)/ZehnerSchluenderThermCond/ZehnerSchluenderThermCond.C
 $(thermCondModels)/noTherm/noThermCond.C
 
 $(forceModels)/forceModel/forceModel.C
@@ -63,6 +64,7 @@ $(forceModels)/particleCellVolume/particleCellVolume.C
 $(forceModels)/fieldTimeAverage/fieldTimeAverage.C
 $(forceModels)/volWeightedAverage/volWeightedAverage.C
 $(forceModels)/BeetstraDrag/BeetstraDrag.C
+$(forceModels)/BeetstraDragPoly/BeetstraDragPoly.C
 $(forceModels)/dSauter/dSauter.C
 $(forceModels)/Fines/Fines.C
 $(forceModels)/Fines/FinesFields.C

src/lagrangian/cfdemParticle/cfdemCloud/cfdemCloud.C

@@ -83,6 +83,9 @@ cfdemCloud::cfdemCloud
     ignore_(false),
     allowCFDsubTimestep_(true),
     limitDEMForces_(false),
+    getParticleDensities_(couplingProperties_.lookupOrDefault<bool>("getParticleDensities",false)),
+    getParticleEffVolFactors_(couplingProperties_.lookupOrDefault<bool>("getParticleEffVolFactors",false)),
+    getParticleTypes_(couplingProperties_.lookupOrDefault<bool>("getParticleTypes",false)),
     modelType_(couplingProperties_.lookup("modelType")),
     positions_(NULL),
     velocities_(NULL),
@@ -95,6 +98,9 @@ cfdemCloud::cfdemCloud
     radii_(NULL),
     voidfractions_(NULL),
     cellIDs_(NULL),
+    particleDensities_(NULL),
+    particleEffVolFactors_(NULL),
+    particleTypes_(NULL),
     particleWeights_(NULL),
     particleVolumes_(NULL),
     particleV_(NULL),
@@ -127,6 +133,19 @@ cfdemCloud::cfdemCloud
         mesh,
         dimensionedScalar("zero", dimensionSet(0,0,-1,0,0), 0)  // 1/s
     ),
+    particleDensityField_
+    (
+        IOobject
+        (
+            "partRho",
+            mesh.time().timeName(),
+            mesh,
+            IOobject::READ_IF_PRESENT,
+            IOobject::NO_WRITE
+        ),
+        mesh,
+        dimensionedScalar("zero", dimensionSet(1,-3,0,0,0), 0.0)
+    ),
     checkPeriodicCells_(false),
     turbulence_
     (
@@ -358,6 +377,9 @@ cfdemCloud::~cfdemCloud()
     dataExchangeM().destroy(particleWeights_,1);
     dataExchangeM().destroy(particleVolumes_,1);
     dataExchangeM().destroy(particleV_,1);
+    if(getParticleDensities_) dataExchangeM().destroy(particleDensities_,1);
+    if(getParticleEffVolFactors_) dataExchangeM().destroy(particleEffVolFactors_,1);
+    if(getParticleTypes_) dataExchangeM().destroy(particleTypes_,1);
 }
 
 // * * * * * * * * * * * * * * * private Member Functions  * * * * * * * * * * * * * //
@@ -369,6 +391,10 @@ void cfdemCloud::getDEMdata()
 
     if(impDEMdragAcc_)
         dataExchangeM().getData("dragAcc","vector-atom",fAcc_); // array is used twice - might be necessary to clean it first
+
+    if(getParticleDensities_) dataExchangeM().getData("density","scalar-atom",particleDensities_);
+    if(getParticleEffVolFactors_) dataExchangeM().getData("effvolfactor","scalar-atom",particleEffVolFactors_);
+    if(getParticleTypes_) dataExchangeM().getData("type","scalar-atom",particleTypes_);
 }
 
 void cfdemCloud::giveDEMdata()
@@ -446,6 +472,25 @@ void cfdemCloud::setParticleForceField()
     );
 }
 
+void cfdemCloud::setScalarAverages()
+{
+    if(!getParticleDensities_) return;
+    if(verbose_) Info << "- setScalarAverage" << endl;
+
+    particleDensityField_.primitiveFieldRef() = 0.0;
+    averagingM().resetWeightFields();
+    averagingM().setScalarAverage
+    (
+        particleDensityField_,
+        particleDensities_,
+        particleWeights_,
+        averagingM().UsWeightField(),
+        NULL
+    );
+
+    if(verbose_) Info << "setScalarAverage done." << endl;
+}
+
 void cfdemCloud::setVectorAverages()
 {
     if(verbose_) Info << "- setVectorAverage(Us,velocities_,weights_)" << endl;
@@ -598,7 +643,8 @@ bool cfdemCloud::evolve
             clockM().stop("setvoidFraction");
 
             // set average particles velocity field
-            clockM().start(20,"setVectorAverage");
+            clockM().start(20,"setAverages");
+            setScalarAverages();
             setVectorAverages();
 
 
@@ -612,7 +658,7 @@ bool cfdemCloud::evolve
             if(!treatVoidCellsAsExplicitForce())
                 smoothingM().smoothenReferenceField(averagingM().UsNext());
 
-            clockM().stop("setVectorAverage");
+            clockM().stop("setAverages");
         }
 
         //============================================
@@ -703,6 +749,9 @@ bool cfdemCloud::reAllocArrays()
         dataExchangeM().allocateArray(particleWeights_,0.,voidFractionM().maxCellsPerParticle());
         dataExchangeM().allocateArray(particleVolumes_,0.,voidFractionM().maxCellsPerParticle());
         dataExchangeM().allocateArray(particleV_,0.,1);
+        if(getParticleDensities_) dataExchangeM().allocateArray(particleDensities_,0.,1);
+        if(getParticleEffVolFactors_) dataExchangeM().allocateArray(particleEffVolFactors_,0.,1);
+        if(getParticleTypes_) dataExchangeM().allocateArray(particleTypes_,0,1);
         arraysReallocated_ = true;
         return true;
     }

src/lagrangian/cfdemParticle/cfdemCloud/cfdemCloud.H

@@ -96,6 +96,12 @@ protected:
 
     bool limitDEMForces_;
 
+    bool getParticleDensities_;
+
+    bool getParticleEffVolFactors_;
+
+    bool getParticleTypes_;
+
     scalar maxDEMForce_;
 
     const word modelType_;
@@ -122,6 +128,12 @@ protected:
 
     mutable int **cellIDs_;
 
+    mutable double **particleDensities_;
+
+    mutable double **particleEffVolFactors_;
+
+    mutable int **particleTypes_;
+
     mutable double **particleWeights_;
 
     mutable double **particleVolumes_;
@@ -162,6 +174,8 @@ protected:
 
     mutable volScalarField ddtVoidfraction_;
 
+    mutable volScalarField particleDensityField_;
+
     mutable Switch checkPeriodicCells_;
 
     const turbulenceModel& turbulence_;
@@ -205,6 +219,8 @@ protected:
 
     virtual void setParticleForceField();
 
+    virtual void setScalarAverages();
+
     virtual void setVectorAverages();
 
 public:
@@ -323,9 +339,17 @@ public:
        virtual inline int maxType() {return -1;}
        virtual inline bool multipleTypesDMax() {return false;}
        virtual inline bool multipleTypesDMin() {return false;}
-       virtual inline double ** particleDensity() const {return NULL;}
+
-       virtual inline int ** particleTypes() const {return NULL;}
+       inline bool getParticleDensities() const;
-       virtual label particleType(label index) const {return -1;}
+       virtual inline double ** particleDensities() const;
+       virtual inline double particleDensity(label index) const;
+
+       inline bool getParticleEffVolFactors() const;
+       virtual inline double particleEffVolFactor(label index) const;
+
+       inline bool getParticleTypes() const;
+       virtual inline int ** particleTypes() const;
+       virtual inline label particleType(label index) const;
 
        //access to the particle's rotation and torque data
        virtual inline double ** DEMTorques() const {return NULL;}

src/lagrangian/cfdemParticle/cfdemCloud/cfdemCloudI.H

@@ -216,6 +216,58 @@ inline double cfdemCloud::d32(bool recalc)
     return d32_;
 }
 
+inline bool cfdemCloud::getParticleDensities() const
+{
+    return getParticleDensities_;
+}
+
+inline double ** cfdemCloud::particleDensities() const
+{
+    return particleDensities_;
+}
+
+inline double cfdemCloud::particleDensity(label index) const
+{
+    if(!getParticleDensities_) return -1.0;
+    else
+    {
+        return particleDensities_[index][0];
+    }
+}
+
+inline bool cfdemCloud::getParticleEffVolFactors() const
+{
+    return getParticleEffVolFactors_;
+}
+
+inline double cfdemCloud::particleEffVolFactor(label index) const
+{
+    if(!getParticleEffVolFactors_) return -1.0;
+    else
+    {
+        return particleEffVolFactors_[index][0];
+    }
+}
+
+inline bool cfdemCloud::getParticleTypes() const
+{
+    return getParticleTypes_;
+}
+
+inline int ** cfdemCloud::particleTypes() const
+{
+    return particleTypes_;
+}
+
+inline label cfdemCloud::particleType(label index) const
+{
+    if(!getParticleDensities_) return -1;
+    else
+    {
+        return particleTypes_[index][0];
+    }
+}
+
 inline int cfdemCloud::numberOfParticles() const
 {
     return numberOfParticles_;

src/lagrangian/cfdemParticle/derived/cfdemCloudEnergy/cfdemCloudEnergy.C

@@ -169,6 +169,12 @@ void cfdemCloudEnergy::postFlow()
         energyModel_[i]().postFlow();
 }
 
+void cfdemCloudEnergy::solve()
+{
+    for (int i=0;i<nrEnergyModels();i++)
+        energyModel_[i]().solve();
+}
+
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 } // End namespace Foam

src/lagrangian/cfdemParticle/derived/cfdemCloudEnergy/cfdemCloudEnergy.H

@@ -108,6 +108,8 @@ public:
 
        void postFlow();
 
+       void solve();
+
 };
 
 

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

@@ -104,6 +104,8 @@ public:
 	virtual void calcEnergyContribution() = 0;
 	
 	virtual void postFlow() {}
+
+	virtual void solve() {}
 	
 	scalar Cp() const;
 

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

@@ -43,8 +43,10 @@ heatTransferGunn::heatTransferGunn
 :
     energyModel(dict,sm),
     propsDict_(dict.subDict(typeName + "Props")),
+    expNusselt_(propsDict_.lookupOrDefault<bool>("expNusselt",false)),
     interpolation_(propsDict_.lookupOrDefault<bool>("interpolation",false)),
     verbose_(propsDict_.lookupOrDefault<bool>("verbose",false)),
+    implicit_(propsDict_.lookupOrDefault<bool>("implicit",true)),
     QPartFluidName_(propsDict_.lookupOrDefault<word>("QPartFluidName","QPartFluid")),
     QPartFluid_
     (   IOobject
@@ -58,17 +60,31 @@ heatTransferGunn::heatTransferGunn
         sm.mesh(),
         dimensionedScalar("zero", dimensionSet(1,-1,-3,0,0,0,0), 0.0)
     ),
+    QPartFluidCoeffName_(propsDict_.lookupOrDefault<word>("QPartFluidCoeffName","QPartFluidCoeff")),
+    QPartFluidCoeff_
+    (   IOobject
+        (
+            QPartFluidCoeffName_,
+            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)
+    ),
     partTempField_
     (   IOobject
         (
-            "particleTemp",
+            "partTemp",
             sm.mesh().time().timeName(),
             sm.mesh(),
-            IOobject::NO_READ,
+            IOobject::READ_IF_PRESENT,
             IOobject::NO_WRITE
         ),
         sm.mesh(),
-        dimensionedScalar("zero", dimensionSet(0,0,0,1,0,0,0), 0.0)
+        dimensionedScalar("zero", dimensionSet(0,0,0,1,0,0,0), 0.0),
+        "zeroGradient"
     ),
     partRelTempField_
     (   IOobject
@@ -108,6 +124,8 @@ heatTransferGunn::heatTransferGunn
     ),
     partRefTemp_("partRefTemp", dimensionSet(0,0,0,1,0,0,0), 0.0),
     calcPartTempField_(propsDict_.lookupOrDefault<bool>("calcPartTempField",false)),
+    calcPartTempAve_(propsDict_.lookupOrDefault<bool>("calcPartTempAve",false)),
+    partTempAve_(0.0),
     tempFieldName_(propsDict_.lookupOrDefault<word>("tempFieldName","T")),
     tempField_(sm.mesh().lookupObject<volScalarField> (tempFieldName_)),
     voidfractionFieldName_(propsDict_.lookupOrDefault<word>("voidfractionFieldName","voidfraction")),
@@ -121,32 +139,47 @@ heatTransferGunn::heatTransferGunn
     partTemp_(NULL),
     partHeatFluxName_(propsDict_.lookup("partHeatFluxName")),
     partHeatFlux_(NULL),
+    partHeatFluxCoeff_(NULL),
     partRe_(NULL),
     partNu_(NULL)
 {
-     allocateMyArrays();
+    allocateMyArrays();
 
     if (propsDict_.found("maxSource"))
     {
         maxSource_=readScalar(propsDict_.lookup ("maxSource"));
         Info << "limiting eulerian source field to: " << maxSource_ << endl;
     }
+
     if (calcPartTempField_)
     {
+        calcPartTempAve_ = true;
         if (propsDict_.found("partRefTemp"))
+        {
             partRefTemp_.value()=readScalar(propsDict_.lookup ("partRefTemp"));
+        }
         partTempField_.writeOpt() = IOobject::AUTO_WRITE;
         partRelTempField_.writeOpt() = IOobject::AUTO_WRITE;
         partTempField_.write();
         partRelTempField_.write();
         Info <<  "Particle temperature field activated." << endl;
     }
+
+    if (!implicit_)
+    {
+        QPartFluidCoeff_.writeOpt() = IOobject::NO_WRITE;
+    }
+
     if (verbose_)
     {
         ReField_.writeOpt() = IOobject::AUTO_WRITE;
         NuField_.writeOpt() = IOobject::AUTO_WRITE;
         ReField_.write();
         NuField_.write();
+        if (expNusselt_)
+        {
+          FatalError <<"Cannot read and create NuField at the same time!\n" << abort(FatalError);
+        }
     }
 }
 
@@ -159,6 +192,10 @@ heatTransferGunn::~heatTransferGunn()
     particleCloud_.dataExchangeM().destroy(partHeatFlux_,1);
     particleCloud_.dataExchangeM().destroy(partRe_,1);
     particleCloud_.dataExchangeM().destroy(partNu_,1);
+    if (implicit_)
+    {
+        particleCloud_.dataExchangeM().destroy(partHeatFluxCoeff_,1);
+    }
 }
 
 // * * * * * * * * * * * * * * * private Member Functions  * * * * * * * * * * * * * //
@@ -168,6 +205,10 @@ void heatTransferGunn::allocateMyArrays() const
     double initVal=0.0;
     particleCloud_.dataExchangeM().allocateArray(partTemp_,initVal,1);  // field/initVal/with/lenghtFromLigghts
     particleCloud_.dataExchangeM().allocateArray(partHeatFlux_,initVal,1);
+    if(implicit_)
+    {
+        particleCloud_.dataExchangeM().allocateArray(partHeatFluxCoeff_,initVal,1);
+    }
 
     if(verbose_)
     {
@@ -227,6 +268,7 @@ void heatTransferGunn::calcEnergyContribution()
     scalar Rep(0);
     scalar Pr(0);
     scalar Nup(0);
+    scalar Tsum(0.0);
 
 
     interpolationCellPoint<scalar> voidfractionInterpolator_(voidfraction_);
@@ -265,13 +307,12 @@ void heatTransferGunn::calcEnergyContribution()
 
                 Nup = Nusselt(voidfraction, Rep, Pr);
 
-
+                Tsum += partTemp_[index][0];
                 scalar h = kf0_ * Nup / ds;
                 scalar As = ds * ds * M_PI; // surface area of sphere
 
                 // calc convective heat flux [W]
                 heatFlux(index, h, As, Tfluid);
-                heatFluxCoeff(index, h, As);
 
                 if(verbose_)
                 {
@@ -295,6 +336,16 @@ void heatTransferGunn::calcEnergyContribution()
             }
     }
 
+    // gather particle temperature sums and obtain average
+    if(calcPartTempAve_)
+    {
+        reduce(Tsum, sumOp<scalar>());
+        partTempAve_ = Tsum / particleCloud_.numberOfParticles();
+        Info << "mean particle temperature = " << partTempAve_ << endl;
+    }
+
+    if(calcPartTempField_) partTempField();
+
     particleCloud_.averagingM().setScalarSum
     (
         QPartFluid_,
@@ -305,6 +356,21 @@ void heatTransferGunn::calcEnergyContribution()
 
     QPartFluid_.primitiveFieldRef() /= -QPartFluid_.mesh().V();
 
+    if(implicit_)
+    {
+        QPartFluidCoeff_.primitiveFieldRef() = 0.0;
+
+        particleCloud_.averagingM().setScalarSum
+        (
+            QPartFluidCoeff_,
+            partHeatFluxCoeff_,
+            particleCloud_.particleWeights(),
+            NULL
+        );
+
+        QPartFluidCoeff_.primitiveFieldRef() /= -QPartFluidCoeff_.mesh().V();
+    }
+
     if(verbose_)
     {
         ReField_.primitiveFieldRef() = 0.0;
@@ -329,22 +395,22 @@ void heatTransferGunn::calcEnergyContribution()
         );
     }
 
-    // limit source term
+    // limit source term in explicit treatment
-    forAll(QPartFluid_,cellI)
+    if(!implicit_)
     {
-        scalar EuFieldInCell = QPartFluid_[cellI];
+        forAll(QPartFluid_,cellI)
-
-        if(mag(EuFieldInCell) > maxSource_ )
         {
-             Pout << "limiting source term\n"  << endl  ;
+            scalar EuFieldInCell = QPartFluid_[cellI];
-             QPartFluid_[cellI] = sign(EuFieldInCell) * maxSource_;
+
+            if(mag(EuFieldInCell) > maxSource_ )
+            {
+                 Pout << "limiting source term\n"  << endl  ;
+                 QPartFluid_[cellI] = sign(EuFieldInCell) * maxSource_;
+            }
         }
     }
 
     QPartFluid_.correctBoundaryConditions();
-
-    giveData(0);
-
 }
 
 void heatTransferGunn::addEnergyContribution(volScalarField& Qsource) const
@@ -352,6 +418,14 @@ void heatTransferGunn::addEnergyContribution(volScalarField& Qsource) const
     Qsource += QPartFluid_;
 }
 
+void heatTransferGunn::addEnergyCoefficient(volScalarField& Qsource) const
+{
+    if(implicit_)
+    {
+        Qsource += QPartFluidCoeff_;
+    }
+}
+
 scalar heatTransferGunn::Nusselt(scalar voidfraction, scalar Rep, scalar Pr) const
 {
     return (7 - 10 * voidfraction + 5 * voidfraction * voidfraction) *
@@ -362,22 +436,77 @@ scalar heatTransferGunn::Nusselt(scalar voidfraction, scalar Rep, scalar Pr) con
 
 void heatTransferGunn::heatFlux(label index, scalar h, scalar As, scalar Tfluid)
 {
-    partHeatFlux_[index][0] = h * As * (Tfluid - partTemp_[index][0]);
+    scalar hAs = h * As;
-}
+    partHeatFlux_[index][0] = - hAs * partTemp_[index][0];
-
+    if(!implicit_)
-void heatTransferGunn::heatFluxCoeff(label index, scalar h, scalar As)
-{
-    //no heat transfer coefficient in explicit model
-}
-
-void heatTransferGunn::giveData(int call)
-{
-    if(call == 0)
     {
-        Info << "total convective particle-fluid heat flux [W] (Eulerian) = " << gSum(QPartFluid_*1.0*QPartFluid_.mesh().V()) << endl;
+        partHeatFlux_[index][0] += hAs * Tfluid;
-
-        particleCloud_.dataExchangeM().giveData(partHeatFluxName_,"scalar-atom", partHeatFlux_);
     }
+    else
+    {
+        partHeatFluxCoeff_[index][0] = hAs;
+    }
+}
+
+void heatTransferGunn::giveData()
+{
+    Info << "total convective particle-fluid heat flux [W] (Eulerian) = " << gSum(QPartFluid_*1.0*QPartFluid_.mesh().V()) << endl;
+
+    particleCloud_.dataExchangeM().giveData(partHeatFluxName_,"scalar-atom", partHeatFlux_);
+}
+
+void heatTransferGunn::postFlow()
+{
+    if(implicit_)
+    {
+        label cellI;
+        scalar Tfluid(0.0);
+        scalar Tpart(0.0);
+        interpolationCellPoint<scalar> TInterpolator_(tempField_);
+
+        for(int index = 0;index < particleCloud_.numberOfParticles(); ++index)
+        {
+                cellI = particleCloud_.cellIDs()[index][0];
+                if(cellI >= 0)
+                {
+                    if(interpolation_)
+                    {
+                        vector position = particleCloud_.position(index);
+                        Tfluid = TInterpolator_.interpolate(position,cellI);
+                    }
+                    else
+                    {
+                        Tfluid = tempField_[cellI];
+                    }
+
+                    Tpart = partTemp_[index][0];
+                    partHeatFlux_[index][0] = (Tfluid - Tpart) * partHeatFluxCoeff_[index][0];
+                }
+        }
+    }
+    giveData();
+}
+
+scalar heatTransferGunn::aveTpart() const
+{
+    return partTempAve_;
+}
+
+void heatTransferGunn::partTempField()
+{
+    partTempField_.primitiveFieldRef() = 0.0;
+    particleCloud_.averagingM().resetWeightFields();
+    particleCloud_.averagingM().setScalarAverage
+    (
+        partTempField_,
+        partTemp_,
+        particleCloud_.particleWeights(),
+        particleCloud_.averagingM().UsWeightField(),
+        NULL
+    );
+
+    dimensionedScalar aveTemp("aveTemp",dimensionSet(0,0,0,1,0,0,0), partTempAve_);
+    partRelTempField_ = (partTempField_ - aveTemp) / (aveTemp - partRefTemp_);
 }
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 

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

@@ -15,7 +15,7 @@ License
     along with this code.  If not, see <http://www.gnu.org/licenses/>.
 
     Copyright (C) 2015- Thomas Lichtenegger, JKU Linz, Austria
-    
+
     Description
     Correlation for Nusselt number according to
     Gunn, D. J. International Journal of Heat and Mass Transfer 21.4 (1978)
@@ -44,27 +44,39 @@ class heatTransferGunn
 protected:
 
     dictionary propsDict_;
-    
+
+    bool expNusselt_;
+
     bool interpolation_;
-    
+
     bool verbose_;
-    
+
+    bool implicit_;
+
     word QPartFluidName_;
-        
+
     volScalarField QPartFluid_;
-    
+
-    volScalarField partTempField_;
+    word QPartFluidCoeffName_;
-    
+
+    volScalarField QPartFluidCoeff_;
+
+    mutable volScalarField partTempField_;
+
     volScalarField partRelTempField_;
-    
+
     volScalarField ReField_;
-    
+
     volScalarField NuField_;
-    
+
     dimensionedScalar partRefTemp_;
-    
+
     bool calcPartTempField_;
 
+    bool calcPartTempAve_;
+
+    scalar partTempAve_;
+
     word tempFieldName_;
 
     const volScalarField& tempField_;            // ref to temperature field
@@ -78,33 +90,35 @@ protected:
     word velFieldName_;
 
     const volVectorField& U_;
-    
+
     word densityFieldName_;
-    
+
     const volScalarField& rho_;
 
     word partTempName_;
 
-    mutable double **partTemp_;                          // Lagrangian array
+    mutable double **partTemp_;
 
     word partHeatFluxName_;
 
-    mutable double **partHeatFlux_;              // Lagrangian array
+    mutable double **partHeatFlux_;
-    
+
+    mutable double **partHeatFluxCoeff_;
+
     mutable double **partRe_;
-    
+
     mutable double **partNu_;
 
     void allocateMyArrays() const;
-    
+
+    void partTempField();
+
     scalar Nusselt(scalar, scalar, scalar) const;
-    
+
-    virtual void giveData(int);
+    virtual void giveData();
-    
+
     virtual void heatFlux(label, scalar, scalar, scalar);
-    
+
-    virtual void heatFluxCoeff(label, scalar, scalar);
-       
 public:
 
     //- Runtime type information
@@ -128,10 +142,14 @@ public:
     // Member Functions
 
         void addEnergyContribution(volScalarField&) const;
-	
+
-	void addEnergyCoefficient(volScalarField&) const {}
+	void addEnergyCoefficient(volScalarField&) const;
-	
+
 	void calcEnergyContribution();
+
+	void postFlow();
+
+        scalar aveTpart() const;
 };
 
 

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

@@ -0,0 +1,244 @@
+/*---------------------------------------------------------------------------*\
+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 "heatTransferGunnPartField.H"
+#include "addToRunTimeSelectionTable.H"
+#include "thermCondModel.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+defineTypeNameAndDebug(heatTransferGunnPartField, 0);
+
+addToRunTimeSelectionTable(energyModel, heatTransferGunnPartField, dictionary);
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+// Construct from components
+heatTransferGunnPartField::heatTransferGunnPartField
+(
+    const dictionary& dict,
+    cfdemCloudEnergy& sm
+)
+:
+    heatTransferGunn(dict,sm),
+    partCpField_
+    (   
+        IOobject
+        (
+            "partCp",
+            sm.mesh().time().timeName(),
+            sm.mesh(),
+            IOobject::READ_IF_PRESENT,
+            IOobject::NO_WRITE
+        ),
+        sm.mesh(),
+        dimensionedScalar("zero", dimensionSet(0,2,-2,-1,0,0,0), 0.0),
+        "zeroGradient"
+    ),
+    partRhoField_(sm.mesh().lookupObject<volScalarField>("partRho")),
+    typeCp_(propsDict_.lookupOrDefault<scalarList>("specificHeatCapacities",scalarList(1,-1.0))),
+    partCp_(NULL),
+    pTMax_(dimensionedScalar("pTMax",dimensionSet(0,0,0,1,0,0,0), -1.0)),
+    pTMin_(dimensionedScalar("pTMin",dimensionSet(0,0,0,1,0,0,0), -1.0)),
+    thermCondModel_
+    (
+        thermCondModel::New
+        (
+            propsDict_,
+            sm
+        )
+    ),
+    fvOptions(fv::options::New(sm.mesh()))
+{
+    if (!implicit_)
+    {
+        FatalError << "heatTransferGunnPartField requires implicit heat transfer treatment." << abort(FatalError);
+    }
+
+    if (typeCp_[0] < 0.0)
+    {
+        FatalError << "heatTransferGunnPartField: provide list of specific heat capacities." << abort(FatalError);
+    }
+    
+    if (propsDict_.found("pTMax"))
+    {
+        pTMax_.value()=scalar(readScalar(propsDict_.lookup("pTMax"))); 
+    }
+    
+    if (propsDict_.found("pTMin"))
+    {
+        pTMin_.value()=scalar(readScalar(propsDict_.lookup("pTMin"))); 
+    }
+
+    partTempField_.writeOpt() = IOobject::AUTO_WRITE;
+
+    allocateMyArrays();
+}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+heatTransferGunnPartField::~heatTransferGunnPartField()
+{
+    particleCloud_.dataExchangeM().destroy(partCp_,1);
+}
+
+// * * * * * * * * * * * * * * * private Member Functions  * * * * * * * * * * * * * //
+
+void heatTransferGunnPartField::allocateMyArrays() const
+{
+    double initVal=0.0;
+    particleCloud_.dataExchangeM().allocateArray(partCp_,initVal,1);
+}
+// * * * * * * * * * * * * * * * * Member Fct  * * * * * * * * * * * * * * * //
+void heatTransferGunnPartField::calcEnergyContribution()
+{
+    allocateMyArrays();
+    heatTransferGunn::calcEnergyContribution();
+
+    // if heat sources in particles present, pull them here
+
+    // loop over all particles to fill partCp_ based on type
+    label cellI=0;
+    label partType = 0;
+    for(int index = 0;index < particleCloud_.numberOfParticles(); ++index)
+    {
+            cellI = particleCloud_.cellIDs()[index][0];
+            if(cellI >= 0)
+            {
+                partType = particleCloud_.particleType(index);
+                // LIGGGGHTS counts types 1, 2, ..., C++ array starts at 0
+                partCp_[index][0] = typeCp_[partType - 1];
+            }
+    }
+
+    partCpField_.primitiveFieldRef() = 0.0;
+    particleCloud_.averagingM().resetWeightFields();
+    particleCloud_.averagingM().setScalarAverage
+    (
+        partCpField_,
+        partCp_,
+        particleCloud_.particleWeights(),
+        particleCloud_.averagingM().UsWeightField(),
+        NULL
+    );
+
+}
+
+void heatTransferGunnPartField::addEnergyContribution(volScalarField& Qsource) const
+{
+    Qsource -= QPartFluidCoeff_*partTempField_;
+}
+
+void heatTransferGunnPartField::giveData()
+{
+    particleCloud_.dataExchangeM().giveData(partTempName_,"scalar-atom", partTemp_);
+}
+
+void heatTransferGunnPartField::postFlow()
+{
+    label cellI;
+    scalar Tpart(0.0);
+    interpolationCellPoint<scalar> partTInterpolator_(partTempField_);
+
+    particleCloud_.dataExchangeM().allocateArray(partTemp_,0.0,1);
+
+    for(int index = 0;index < particleCloud_.numberOfParticles(); ++index)
+    {
+            cellI = particleCloud_.cellIDs()[index][0];
+            if(cellI >= 0)
+            {
+                if(interpolation_)
+                {
+                    vector position = particleCloud_.position(index);
+                    Tpart = partTInterpolator_.interpolate(position,cellI);
+                }
+                else
+                {
+                    Tpart = partTempField_[cellI];
+                }
+                partTemp_[index][0] = Tpart;
+            }
+    }
+
+    giveData();
+}
+
+void heatTransferGunnPartField::solve()
+{
+    // for some weird reason, the particle-fluid heat transfer fields were defined with a negative sign
+
+    volScalarField alphaP = 1.0 - voidfraction_;
+    volScalarField correctedQPartFluidCoeff(QPartFluidCoeff_);
+    // no heattransfer in empty cells -- for numerical stability, add small amount
+    forAll(correctedQPartFluidCoeff,cellI)
+    {
+        if (-correctedQPartFluidCoeff[cellI] < SMALL) correctedQPartFluidCoeff[cellI] = -SMALL;
+    }
+
+    volScalarField Qsource = correctedQPartFluidCoeff*tempField_;
+    volScalarField partCpEff = alphaP*partRhoField_*partCpField_;
+    volScalarField thCondEff = alphaP*thermCondModel_().thermCond();
+//    thCondEff.correctBoundaryConditions();
+
+
+
+    fvScalarMatrix partTEqn
+    (
+     //   fvm::ddt(partCpEff, partTempField_)
+     // + Qsource
+        Qsource
+      - fvm::Sp(correctedQPartFluidCoeff, partTempField_)
+      - fvm::laplacian(thCondEff,partTempField_)
+     ==
+        fvOptions(partCpEff, partTempField_)
+    );
+  // if transient add time derivative - need particle density and specific heat fields
+  // if sources activated add sources
+  // if convection activated add convection
+
+    partTEqn.relax();
+
+    fvOptions.constrain(partTEqn);
+
+    partTEqn.solve();
+
+    partTempField_.relax();
+
+    fvOptions.correct(partTempField_);
+
+    if (pTMax_.value() > 0.0) partTempField_ = min(partTempField_, pTMax_);
+    if (pTMin_.value() > 0.0) partTempField_ = max(partTempField_, pTMin_);
+
+
+    Info<< "partT max/min : " << max(partTempField_).value() << " " << min(partTempField_).value() << endl;
+
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// ************************************************************************* //

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

@@ -0,0 +1,108 @@
+/*---------------------------------------------------------------------------*\
+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 Nusselt number according to
+    Gunn, D. J. International Journal of Heat and Mass Transfer 21.4 (1978)
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef heatTransferGunnPartField_H
+#define heatTransferGunnPartField_H
+
+#include "fvCFD.H"
+#include "cfdemCloudEnergy.H"
+#include "heatTransferGunn.H"
+#include "fvOptions.H"
+#include "scalarList.H"
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+class thermCondModel;
+/*---------------------------------------------------------------------------*\
+                           Class heatTransferGunnPartField Declaration
+\*---------------------------------------------------------------------------*/
+
+class heatTransferGunnPartField
+:
+    public heatTransferGunn
+{
+private:
+
+    volScalarField partCpField_;
+
+    const volScalarField& partRhoField_;
+
+    scalarList typeCp_;
+
+    mutable double **partCp_;
+
+    dimensionedScalar pTMax_;
+
+    dimensionedScalar pTMin_;
+
+    autoPtr<thermCondModel> thermCondModel_;
+
+    fv::options& fvOptions;
+
+    void allocateMyArrays() const;
+
+    void giveData();
+
+public:
+
+    //- Runtime type information
+    TypeName("heatTransferGunnPartField");
+
+    // Constructors
+
+        //- Construct from components
+        heatTransferGunnPartField
+        (
+            const dictionary& dict,
+            cfdemCloudEnergy& sm
+        );
+
+
+    // Destructor
+
+        virtual ~heatTransferGunnPartField();
+
+
+    // Member Functions
+        void addEnergyContribution(volScalarField&) const;
+
+	void calcEnergyContribution();
+
+	void postFlow();
+
+        void solve();
+
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

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

@@ -12,6 +12,7 @@ License
     along with this code.  If not, see <http://www.gnu.org/licenses/>.
 
     Copyright (C) 2015- Thomas Lichtenegger, JKU Linz, Austria
+    Contributing author: Paul Kieckhefen, TU Hamburg, Germany
 
 \*---------------------------------------------------------------------------*/
 
@@ -49,14 +50,26 @@ BeetstraDrag::BeetstraDrag
 :
     forceModel(dict,sm),
     propsDict_(dict.subDict(typeName + "Props")),
+    multiTypes_(false),
     velFieldName_(propsDict_.lookup("velFieldName")),
     U_(sm.mesh().lookupObject<volVectorField> (velFieldName_)),
     voidfractionFieldName_(propsDict_.lookup("voidfractionFieldName")),
     voidfraction_(sm.mesh().lookupObject<volScalarField> (voidfractionFieldName_)),
+    minVoidfraction_(propsDict_.lookupOrDefault<scalar>("minVoidfraction",0.1)),
     UsFieldName_(propsDict_.lookup("granVelFieldName")),
     UsField_(sm.mesh().lookupObject<volVectorField> (UsFieldName_)),
     scaleDia_(1.),
-    scaleDrag_(1.)
+    typeCG_(propsDict_.lookupOrDefault<scalarList>("coarseGrainingFactors",scalarList(1,1.0))),
+    scaleDrag_(1.),
+    rhoP_(0.),
+    rho_(0.),
+    Lc2_(0.),
+    dPrim_(0.),
+    nuf_(0.),
+    g_(9.81),
+    k_(0.05),
+    useGC_(false),
+    usePC_(false)
 {
     //Append the field names to be probed
     particleCloud_.probeM().initialize(typeName, typeName+".logDat");
@@ -78,10 +91,42 @@ BeetstraDrag::BeetstraDrag
     forceSubM(0).readSwitches();
 
     particleCloud_.checkCG(true);
-    if (propsDict_.found("scale"))
+    if (propsDict_.found("scale") && typeCG_.size()==1)
+    {
+        // if "scale" is specified and there's only one single type, use "scale"
         scaleDia_=scalar(readScalar(propsDict_.lookup("scale")));
+        typeCG_[0] = scaleDia_;
+    }
     if (propsDict_.found("scaleDrag"))
+    {
         scaleDrag_=scalar(readScalar(propsDict_.lookup("scaleDrag")));
+    }
+
+    if (typeCG_.size()>1) multiTypes_ = true;
+
+    if (propsDict_.found("useFilteredDragModel"))
+    {
+        useGC_ = true;
+        g_=propsDict_.lookupOrDefault<scalar>("g", 9.81);
+        dPrim_=scalar(readScalar(propsDict_.lookup("dPrim")));
+        rhoP_=scalar(readScalar(propsDict_.lookup("rhoP")));
+        rho_=scalar(readScalar(propsDict_.lookup("rho")));
+        nuf_=scalar(readScalar(propsDict_.lookup("nuf")));
+	    scalar ut = terminalVelocity(1., dPrim_, nuf_, rho_, rhoP_, g_);
+	    scalar Frp = ut*ut/g_/dPrim_;
+        Lc2_ = ut*ut/g_*pow(Frp, -.6666667); // n is hardcoded as -2/3
+        Info << "using grid coarsening correction with Lc2 = " << Lc2_ << " and ut = " << ut << " and Frp = " << Frp<< endl;
+
+        if (propsDict_.found("useParcelSizeDependentFilteredDrag"))
+        {
+            usePC_ = true;
+            if (propsDict_.found("k"))
+                k_=scalar(readScalar(propsDict_.lookup("k")));
+            Info << "using particle coarsening correction with k = " << k_ << endl;
+        }
+    }
+
+
 
 }
 
@@ -96,9 +141,9 @@ BeetstraDrag::~BeetstraDrag()
 
 void BeetstraDrag::setForce() const
 {
-    if (scaleDia_ > 1)
+    if (typeCG_.size()>1 || typeCG_[0] > 1)
     {
-        Info << "Beetstra using scale = " << scaleDia_ << endl;
+        Info << "Beetstra using scale = " << typeCG_ << endl;
     }
     else if (particleCloud_.cg() > 1)
     {
@@ -119,12 +164,18 @@ void BeetstraDrag::setForce() const
     vector Ur(0,0,0);
     scalar ds(0);
     scalar ds_scaled(0);
-    scalar scaleDia3 = scaleDia_*scaleDia_*scaleDia_;
+    scalar dSauter(0);
+    scalar scaleDia3 = typeCG_[0]*typeCG_[0]*typeCG_[0];
     scalar nuf(0);
     scalar rho(0);
     scalar magUr(0);
     scalar Rep(0);
     scalar localPhiP(0);
+    scalar GCcorr(1.);
+    scalar PCcorr(1.);
+
+    scalar cg = typeCG_[0];
+    label partType = 1;
 
     vector dragExplicit(0,0,0);
     scalar dragCoefficient(0);
@@ -146,43 +197,70 @@ void BeetstraDrag::setForce() const
             if (cellI > -1) // particle found
             {
 
-                if( forceSubM(0).interpolation() )
+                if ( forceSubM(0).interpolation() )
                 {
                     position     = particleCloud_.position(index);
                     voidfraction = voidfractionInterpolator_.interpolate(position,cellI);
                     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 > 1.00) voidfraction = 1.0;
-                    if(voidfraction<0.10) voidfraction = 0.10;
+                    if (voidfraction < minVoidfraction_) voidfraction = minVoidfraction_;
                 }
                 else
                 {
                     voidfraction = voidfraction_[cellI];
                     Ufluid = U_[cellI];
                 }
+                // in case a fines phase is present, void fraction needs to be adapted
+                adaptVoidfraction(voidfraction, cellI);
+
+                if (multiTypes_)
+                {
+                    partType = particleCloud_.particleType(index);
+                    cg = typeCG_[partType - 1];
+                    scaleDia3 = cg*cg*cg;
+                }
 
                 Us = particleCloud_.velocity(index);
                 Ur = Ufluid-Us;
                 magUr = mag(Ur);
                 ds = 2*particleCloud_.radius(index);
-                ds_scaled = ds/scaleDia_;
+                ds_scaled = ds/cg;
+                dSauter = meanSauterDiameter(ds_scaled, cellI);
+
                 rho = rhoField[cellI];
                 nuf = nufField[cellI];
 
-                Rep=0.0;
                 localPhiP = 1.0f-voidfraction+SMALL;
 
                 // calc particle's drag coefficient (i.e., Force per unit slip velocity and Stokes drag)
 
-                Rep=ds_scaled*voidfraction*magUr/nuf+SMALL;
+                Rep=dSauter*voidfraction*magUr/nuf + SMALL;
-                dragCoefficient = 10.0*localPhiP/(voidfraction*voidfraction) +
+
-                                  voidfraction*voidfraction*(1.0+1.5*Foam::sqrt(localPhiP)) +
+                dragCoefficient = F(voidfraction, Rep)
-                                  0.413*Rep/(24*voidfraction*voidfraction)*(1.0/voidfraction+3*voidfraction*localPhiP+8.4*Foam::pow(Rep,-0.343))/
+                                   *3*M_PI*nuf*rho*voidfraction
-                                  (1+Foam::pow(10,3*localPhiP)*Foam::pow(Rep,-0.5*(1+4*localPhiP)));
+                                   *effDiameter(ds_scaled, cellI, index)
+                                   *scaleDia3*scaleDrag_;
+
+                // calculate filtering corrections
+                if (useGC_)
+                {
+                    GCcorr = 1.-h(localPhiP)
+                                /(  a(localPhiP)
+                                    *Lc2_
+                                    /std::cbrt(U_.mesh().V()[cellI])
+                                  + 1.
+                                 );
+                    if (usePC_)
+                    {
+                        PCcorr = Foam::exp(k_*(1.-cg));
+                    }
+                }
+
+                // apply filtering corrections
+                dragCoefficient *= GCcorr*PCcorr;
 
-                // calc particle's drag
-                dragCoefficient *= 3*M_PI*ds_scaled*nuf*rho*voidfraction*scaleDia3*scaleDrag_;
                 if (modelType_=="B")
                     dragCoefficient /= voidfraction;
 
@@ -198,11 +276,13 @@ void BeetstraDrag::setForce() const
                     Pout << "Us = " << Us << endl;
                     Pout << "Ur = " << Ur << endl;
                     Pout << "ds = " << ds << endl;
-                    Pout << "ds/scale = " << ds/scaleDia_ << endl;
+                    Pout << "ds/scale = " << ds/cg << endl;
                     Pout << "rho = " << rho << endl;
                     Pout << "nuf = " << nuf << endl;
                     Pout << "voidfraction = " << voidfraction << endl;
                     Pout << "Rep = " << Rep << endl;
+                    Pout << "GCcorr = " << GCcorr << endl;
+		            Pout << "PCcorr = " << PCcorr << endl;
                     Pout << "drag = " << drag << endl;
                 }
 
@@ -223,7 +303,149 @@ void BeetstraDrag::setForce() const
     }
 }
 
+/*********************************************************
+ * "Drag Force of Intermediate Reynolds Number Flow Past *
+ *  Mono- and Bidisperse Arrays of Spheres", eq. 16      *
+ *  R Beetstra, M. A. van der Hoef, JAM Kuipers          *
+ *  AIChE Journal 53(2) (2007)                           *
+ *********************************************************/
+double BeetstraDrag::F(double voidfraction, double Rep) const
+{
+    double localPhiP = max(SMALL,min(1.-SMALL,1.-voidfraction));
+    return   10.0*localPhiP/(voidfraction*voidfraction)
+           + voidfraction*voidfraction*(1.0+1.5*Foam::sqrt(localPhiP))
+           + 0.413*Rep/(24*voidfraction*voidfraction)
+                       *(1.0/voidfraction
+                          +3*voidfraction*localPhiP
+                          +8.4*Foam::pow(Rep,-0.343)
+                        )
+                        /(1+Foam::pow(10,3*localPhiP)
+                            *Foam::pow(Rep,-0.5*(1+4*localPhiP))
+                         );
 
+}
+ 
+
+/*********************************************************
+ * "A drag model for filtered Euler-Lagange simulations  *
+ *  of clustered gas-particle suspension",               *
+ *  S. Radl, S. Sundaresan,                              *
+ *  Chemical Engineering Science 117 (2014).             *
+ *********************************************************/
+double BeetstraDrag::a(double phiP) const
+{
+    double a0m = 0.;
+    double a1m = 0.;
+    double a2m = 0.;
+    double a3m = 0.;
+    double phipam = 0.;
+    if      (phiP < 0.016)
+    {
+      a0m = 21.51;
+    }
+    else if (phiP < 0.100)
+    {
+      a0m =  1.96; a1m =  29.40; a2m =  164.91; a3m = -1923.;
+    }
+    else if (phiP < 0.180)
+    {
+      a0m =  4.63; a1m =   4.68; a2m = -412.04; a3m =  2254.; phipam = 0.10;
+    }
+    else if (phiP < 0.250)
+    {
+      a0m =  3.52; a1m = -17.99; a2m =  128.80; a3m =  -603.; phipam = 0.18;
+    }
+    else if (phiP < 0.400)
+    {
+      a0m =  2.68; a1m =  -8.20; a2m =    2.18; a3m = 112.33; phipam = 0.25;
+    }
+    else
+    {
+      a0m =  1.79;
+    }
+    const double phiP_minus_phipam = phiP-phipam;
+    return a0m + a1m*phiP_minus_phipam + a2m*phiP_minus_phipam*phiP_minus_phipam
+           + a3m*phiP_minus_phipam*phiP_minus_phipam*phiP_minus_phipam;
+}
+
+double BeetstraDrag::h(double phiP) const
+{
+  double h0m = 0.;
+  double h1m = 0.;
+  double h2m = 0.;
+  double h3m = 0.;
+  double phiphm = 0.;
+
+  if      (phiP < 0.03)
+  {
+                 h1m = 7.97;
+  }
+  else if (phiP < 0.08)
+  {
+    h0m = 0.239; h1m =  4.640; h2m =   -4.410; h3m =  253.630; phiphm = 0.03;
+  }
+  else if (phiP < 0.12)
+  {
+    h0m = 0.492; h1m =  6.100; h2m =   33.630; h3m = -789.600; phiphm = 0.08;
+  }
+  else if (phiP < 0.18)
+  {
+    h0m = 0.739; h1m =  5.010; h2m =  -61.100; h3m =  310.800; phiphm = 0.12;
+  }
+  else if (phiP < 0.34)
+  {
+    h0m = 0.887; h1m =  1.030; h2m =   -5.170; h3m =    5.990; phiphm = 0.18;
+  }
+  else if (phiP < 0.48)
+  {
+    h0m = 0.943; h1m = -0.170; h2m =   -2.290; h3m =   -9.120; phiphm = 0.34;
+  }
+  else if (phiP < 0.55)
+  {
+    h0m = 0.850; h1m = -1.350; h2m =   -6.130; h3m = -132.600; phiphm = 0.48;
+  }
+  else
+  {
+    h0m = 0.680; h1m = -2.340; h2m = -225.200;                 phiphm = 0.55;
+  }
+  const double phiP_minus_phiphm = phiP-phiphm;
+  return h0m + h1m*phiP_minus_phiphm + h2m*phiP_minus_phiphm*phiP_minus_phiphm
+         + h3m*phiP_minus_phiphm*phiP_minus_phiphm*phiP_minus_phiphm;
+}
+
+double BeetstraDrag::terminalVelocity(double voidfraction, double dp, double nuf, double rhof, double rhop, double g) const
+{
+    scalar u0(dp*dp*fabs(rhof-rhop)*g/18./rhof/nuf);
+    scalar Re(u0*dp/nuf);
+    scalar res(1.);
+    scalar u(u0);
+    scalar Fi(0);
+    scalar CdSt(0);
+    Info << "uo: " << u0<<endl;
+    int i = 0;
+    while ((res > 1.e-6) && (i<100))
+    {
+        Info << "Iteration " << i;
+        u0 = u;
+        Info << ", u0 = " << u0;
+        CdSt = 24/Re;
+        Info << ", CdSt = " << CdSt;
+        Fi = F(voidfraction, Re);
+        Info << ", F = ";
+        u = sqrt(1.333333333*fabs(rhof-rhop)*g*dp
+                  /(CdSt*voidfraction*Fi*rhof)
+            );
+        Info << ", u = " << u;
+        Re = fabs(u)*dp/nuf*voidfraction;
+        res = fabs((u-u0)/u);
+        Info << "Res: " << res << endl;
+        i++;
+    }
+    if (res >1.e-6)
+        FatalError << "Terminal velocity calculation diverged!" << endl;
+
+    return u;
+}
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 } // End namespace Foam

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

@@ -40,9 +40,11 @@ class BeetstraDrag
 :
     public forceModel
 {
-private:
+protected:
     dictionary propsDict_;
 
+    bool multiTypes_;
+
     word velFieldName_;
 
     const volVectorField& U_;
@@ -51,14 +53,52 @@ private:
 
     const volScalarField& voidfraction_;
 
+    const scalar minVoidfraction_;
+
     word UsFieldName_;
 
     const volVectorField& UsField_;
 
     mutable scalar scaleDia_;
 
+    scalarList typeCG_;
+
     mutable scalar scaleDrag_;
 
+    mutable scalar rhoP_;
+
+    mutable scalar rho_;
+
+    mutable scalar Lc2_;
+
+    mutable scalar dPrim_;
+
+    mutable scalar nuf_;
+
+    mutable scalar g_;
+
+    mutable scalar k_;
+
+    bool useGC_;
+
+    bool usePC_;
+
+    virtual void adaptVoidfraction(double&, label) const {}
+
+    virtual scalar effDiameter(double d, label cellI, label index) const {return d;}
+    
+    virtual scalar meanSauterDiameter(double d, label cellI) const {return d;}
+
+    double F(double, double) const;
+
+    double terminalVelocity(double, double, double, double, double, double) const; 
+
+    double a(double) const;
+
+    double h(double) const;
+
+
+
 public:
 
     //- Runtime type information

src/lagrangian/cfdemParticle/subModels/forceModel/BeetstraDragPoly/BeetstraDragPoly.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 "BeetstraDragPoly.H"
+#include "addToRunTimeSelectionTable.H"
+#include "averagingModel.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+defineTypeNameAndDebug(BeetstraDragPoly, 0);
+
+addToRunTimeSelectionTable
+(
+    forceModel,
+    BeetstraDragPoly,
+    dictionary
+);
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+// Construct from components
+BeetstraDragPoly::BeetstraDragPoly
+(
+    const dictionary& dict,
+    cfdemCloud& sm
+)
+:
+    BeetstraDrag(dict,sm),
+    fines_(propsDict_.lookupOrDefault<bool>("fines",false)),
+    dFine_(1.0)
+{
+    // if fines are present, take mixture dSauter, otherwise normal dSauter
+    if (fines_)
+    {
+        dFine_ = readScalar(propsDict_.lookup("dFine"));
+        volScalarField& alphaP(const_cast<volScalarField&>(sm.mesh().lookupObject<volScalarField> ("alphaP")));
+        alphaP_.set(&alphaP);
+        volScalarField& alphaSt(const_cast<volScalarField&>(sm.mesh().lookupObject<volScalarField> ("alphaSt")));
+        alphaSt_.set(&alphaSt);
+        volScalarField& dSauter(const_cast<volScalarField&>(sm.mesh().lookupObject<volScalarField> ("dSauterMix")));
+        dSauter_.set(&dSauter);
+    }
+    else
+    {
+        volScalarField& dSauter(const_cast<volScalarField&>(sm.mesh().lookupObject<volScalarField> ("dSauter")));
+        dSauter_.set(&dSauter);
+    }
+}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+BeetstraDragPoly::~BeetstraDragPoly()
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+void BeetstraDragPoly::adaptVoidfraction(double& voidfraction, label cellI) const
+{
+    voidfraction -= alphaSt_()[cellI];
+    if (voidfraction < minVoidfraction_) voidfraction = minVoidfraction_;
+}
+
+scalar BeetstraDragPoly::effDiameter(double d, label cellI, label index) const
+{
+    scalar dS = dSauter_()[cellI];
+    scalar effD = d*d / dS + 0.064*d*d*d*d / (dS*dS*dS);
+    
+    if (fines_)
+    {
+        scalar fineCorr = dFine_*dFine_ / dS + 0.064*dFine_*dFine_*dFine_*dFine_ / (dS*dS*dS);
+        fineCorr *= d*d*d / (dFine_*dFine_*dFine_) * alphaSt_()[cellI] / alphaP_()[cellI];
+        effD += fineCorr;
+    }
+    
+    if (particleCloud_.getParticleEffVolFactors()) 
+    {
+        scalar effVolFac = particleCloud_.particleEffVolFactor(index);
+        effD *= effVolFac;
+    }
+
+    return effD;
+}
+
+scalar BeetstraDragPoly::meanSauterDiameter(double d, label cellI) const
+{
+    return dSauter_()[cellI];
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// ************************************************************************* //

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

@@ -0,0 +1,93 @@
+/*---------------------------------------------------------------------------*\
+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
+    drag law for monodisperse systems according to
+    Beetstra et al. AIChE J 53.2 (2007)
+
+SourceFiles
+    BeetstraDragPoly.C
+\*---------------------------------------------------------------------------*/
+
+#ifndef BeetstraDragPoly_H
+#define BeetstraDragPoly_H
+
+#include "BeetstraDrag.H"
+#include "interpolationCellPoint.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                           Class BeetstraDragPoly Declaration
+\*---------------------------------------------------------------------------*/
+
+class BeetstraDragPoly
+:
+    public BeetstraDrag
+{
+protected:
+
+    const bool fines_;
+    
+    scalar dFine_;
+
+    autoPtr<volScalarField> alphaP_;
+    
+    autoPtr<volScalarField> alphaSt_;
+
+    autoPtr<volScalarField> dSauter_;
+
+    void adaptVoidfraction(double&, label) const;
+
+    scalar effDiameter(double, label, label) const;
+    
+    scalar meanSauterDiameter(double, label) const;
+
+public:
+
+    //- Runtime type information
+    TypeName("BeetstraDragPoly");
+
+
+    // Constructors
+
+        //- Construct from components
+        BeetstraDragPoly
+        (
+            const dictionary& dict,
+            cfdemCloud& sm
+        );
+
+    // Destructor
+
+        ~BeetstraDragPoly();
+
+
+    // Member Functions
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

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

@@ -53,10 +53,10 @@ dSauter::dSauter
 :
     forceModel(dict,sm),
     propsDict_(dict.subDict(typeName + "Props")),
+    multiTypes_(false),
     d2_(NULL),
     d3_(NULL),
-    scaleDia_(1.0),
+    typeCG_(propsDict_.lookupOrDefault<scalarList>("coarseGrainingFactors",scalarList(1,1.0))),
-    scaleDiaDist_(1.0),
     d2Field_
     (   IOobject
         (
@@ -95,17 +95,13 @@ dSauter::dSauter
         "zeroGradient"
     )
 {
+    if (typeCG_.size()>1) multiTypes_ = true;
     allocateMyArrays();
     dSauter_.write();
 
 
     // init force sub model
     setForceSubModels(propsDict_);
-
-    if (propsDict_.found("scaleCG"))
-        scaleDia_ = scalar(readScalar(propsDict_.lookup("scaleCG")));
-    if (propsDict_.found("scaleDist"))
-        scaleDiaDist_ = scalar(readScalar(propsDict_.lookup("scaleDist")));
 }
 
 
@@ -131,30 +127,37 @@ void dSauter::allocateMyArrays() const
 
 void dSauter::setForce() const
 {
-    if (scaleDia_ > 1)
+    if (typeCG_.size()>1 || typeCG_[0] > 1.0)
     {
-        Info << "dSauter using scaleCG = " << scaleDia_ << endl;
+        Info << "dSauter using CG factor(s) = " << typeCG_ << endl;
-    }
-    else if (particleCloud_.cg() > 1)
-    {
-        scaleDia_ = particleCloud_.cg();
-        Info << "dSauter using scaleCG from liggghts cg = " << scaleDia_ << endl;
     }
 
     allocateMyArrays();
 
-    label cellI=0;
+    label cellI = 0;
-    scalar ds(0);
+    label partType = 1;
-    scalar scale = scaleDiaDist_/scaleDia_;
+    scalar cg = typeCG_[0];
+    scalar ds = 0.0;
+    scalar effVolFac = 1.0;
+
 
     for(int index = 0; index < particleCloud_.numberOfParticles(); ++index)
     {
         cellI = particleCloud_.cellIDs()[index][0];
         if (cellI >= 0)
         {
+            if (particleCloud_.getParticleEffVolFactors()) 
+            {
+                effVolFac = particleCloud_.particleEffVolFactor(index);
+            }
+            if (multiTypes_) 
+            {
+                partType = particleCloud_.particleType(index);
+                cg = typeCG_[partType - 1];
+            }
             ds = particleCloud_.d(index);
-            d2_[index][0] = ds*ds;
+            d2_[index][0] = ds*ds*effVolFac*cg;
-            d3_[index][0] = ds*ds*ds;
+            d3_[index][0] = ds*ds*ds*effVolFac;
         }
     }
 
@@ -181,7 +184,7 @@ void dSauter::setForce() const
     {
         if (d2Field_[cellI] > ROOTVSMALL)
         {
-            dSauter_[cellI] = d3Field_[cellI] / d2Field_[cellI] * scale;
+            dSauter_[cellI] = d3Field_[cellI] / d2Field_[cellI];
         }
         else
         {

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

@@ -43,13 +43,13 @@ private:
 
     dictionary propsDict_;
 
+    bool multiTypes_;
+
     mutable double **d2_;
     
     mutable double **d3_;
 
-    mutable scalar scaleDia_;      // scaling due to coarse graining
+    scalarList typeCG_;
-
-    mutable scalar scaleDiaDist_;  // scaling due to modified particle size distribution
     
     mutable volScalarField d2Field_;
     

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

@@ -0,0 +1,225 @@
+/*---------------------------------------------------------------------------*\
+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 "ZehnerSchluenderThermCond.H"
+#include "addToRunTimeSelectionTable.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+defineTypeNameAndDebug(ZehnerSchluenderThermCond, 0);
+
+addToRunTimeSelectionTable
+(
+    thermCondModel,
+    ZehnerSchluenderThermCond,
+    dictionary
+);
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+// Construct from components
+ZehnerSchluenderThermCond::ZehnerSchluenderThermCond
+(
+    const dictionary& dict,
+    cfdemCloud& sm
+)
+:
+    thermCondModel(dict,sm),
+    propsDict_(dict.subDict(typeName + "Props")),
+    partKsField_
+    (
+        IOobject
+        (
+            "partKs",
+            sm.mesh().time().timeName(),
+            sm.mesh(),
+            IOobject::READ_IF_PRESENT,
+            IOobject::NO_WRITE
+        ),
+        sm.mesh(),
+        dimensionedScalar("one", dimensionSet(1, 1, -3, -1,0,0,0), 1.0),
+	"zeroGradient"
+    ),
+    voidfractionFieldName_(propsDict_.lookupOrDefault<word>("voidfractionFieldName","voidfraction")),
+    voidfraction_(sm.mesh().lookupObject<volScalarField> (voidfractionFieldName_)),
+    typeKs_(propsDict_.lookupOrDefault<scalarList>("thermalConductivities",scalarList(1,-1.0))),
+    partKs_(NULL),
+    wallQFactorName_(propsDict_.lookupOrDefault<word>("wallQFactorName","wallQFactor")),
+    wallQFactor_
+    (   IOobject
+        (
+            wallQFactorName_,
+            sm.mesh().time().timeName(),
+            sm.mesh(),
+            IOobject::READ_IF_PRESENT,
+            IOobject::AUTO_WRITE
+        ),
+        sm.mesh(),
+        dimensionedScalar("zero", dimensionSet(0,0,0,0,0,0,0), 1.0)
+    ),
+    hasWallQFactor_(false)
+{
+    if (typeKs_[0] < 0.0)
+    {
+        FatalError << "ZehnerSchluenderThermCond: provide list of thermal conductivities." << abort(FatalError);
+    }
+
+    if (typeKs_.size() > 1) allocateMyArrays();
+    else
+    {
+        partKsField_ *= typeKs_[0];
+    }
+
+    if (wallQFactor_.headerOk())
+    {
+        Info << "Found field for scaling wall heat flux.\n" << endl;
+        hasWallQFactor_ = true;
+    }
+}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+ZehnerSchluenderThermCond::~ZehnerSchluenderThermCond()
+{
+    if (typeKs_.size() > 1) particleCloud_.dataExchangeM().destroy(partKs_,1);
+}
+
+
+void ZehnerSchluenderThermCond::allocateMyArrays() const
+{
+    double initVal=0.0;
+    particleCloud_.dataExchangeM().allocateArray(partKs_,initVal,1);
+}
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+tmp<volScalarField> ZehnerSchluenderThermCond::thermCond() const
+{
+    tmp<volScalarField> tvf
+    (
+        new volScalarField
+        (
+            IOobject
+            (
+                "tmpThCond",
+                voidfraction_.instance(),
+                voidfraction_.mesh(),
+                IOobject::NO_READ,
+                IOobject::NO_WRITE
+            ),
+            voidfraction_.mesh(),
+            dimensionedScalar("zero", dimensionSet(1,1,-3,-1,0,0,0), 0.0),
+            "zeroGradient"
+        )
+    );
+
+    volScalarField& svf = tvf.ref();
+
+    calcPartKsField();
+    scalar A = 0.0;
+    scalar B = 0.0;
+    scalar C = 0.0;
+    scalar k = 0.0;
+    scalar InvOnemBoA = 0.0;
+    scalar voidfraction = 0.0;
+    scalar w = 7.26e-3;
+
+    forAll(svf, cellI)
+    {
+        voidfraction = voidfraction_[cellI];
+        if(voidfraction > 1.0 - SMALL || partKsField_[cellI] < SMALL) svf[cellI] = 0.0;
+        else
+        {
+            A = partKsField_[cellI]/kf0_.value();
+            B = 1.25 * Foam::pow((1 - voidfraction) / voidfraction, 1.11);
+            InvOnemBoA = 1.0/(1.0 - B/A);
+            C = (A - 1) * InvOnemBoA * InvOnemBoA * B/A * log(A/B) - (B - 1) * InvOnemBoA - 0.5 * (B + 1);
+            C *= 2.0 * InvOnemBoA;
+            k = Foam::sqrt(1 - voidfraction) * (w * A + (1 - w) * C) * kf0_.value();
+            svf[cellI] = k / (1 - voidfraction);
+        }
+    }
+
+    svf.correctBoundaryConditions();
+
+    // if a wallQFactor field is present, use it to scale heat transport through a patch
+    if (hasWallQFactor_)
+    {
+        wallQFactor_.correctBoundaryConditions();
+        forAll(wallQFactor_.boundaryField(), patchi)
+            svf.boundaryFieldRef()[patchi] *= wallQFactor_.boundaryField()[patchi];
+    }
+
+    return tvf;
+}
+
+tmp<volScalarField> ZehnerSchluenderThermCond::thermDiff() const
+{
+    FatalError << "ZehnerSchluenderThermCond does not provide thermal diffusivity." << abort(FatalError);
+    return tmp<volScalarField>(NULL);
+}
+
+void ZehnerSchluenderThermCond::calcPartKsField() const
+{
+    if (typeKs_.size() <= 1) return;
+
+    if (!particleCloud_.getParticleTypes())
+    {
+        FatalError << "ZehnerSchluenderThermCond needs data for more than one type, but types are not communicated." << abort(FatalError);
+    }
+
+    allocateMyArrays();
+    label cellI=0;
+    label partType = 0;
+    for(int index = 0;index < particleCloud_.numberOfParticles(); ++index)
+    {
+            cellI = particleCloud_.cellIDs()[index][0];
+            if(cellI >= 0)
+            {
+                partType = particleCloud_.particleType(index);
+                // LIGGGGHTS counts types 1, 2, ..., C++ array starts at 0
+                partKs_[index][0] = typeKs_[partType - 1];
+            }
+    }
+
+    partKsField_.primitiveFieldRef() = 0.0;
+    particleCloud_.averagingM().resetWeightFields();
+    particleCloud_.averagingM().setScalarAverage
+    (
+        partKsField_,
+        partKs_,
+        particleCloud_.particleWeights(),
+        particleCloud_.averagingM().UsWeightField(),
+        NULL
+    );
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// ************************************************************************* //

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

@@ -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
+
+
+Description
+    thermal conductivity of PARTICLE phase in presence of a fluid according to
+    Zehner and Schluender (1970)
+
+Class
+    ZehnerSchluenderThermCond
+
+SourceFiles
+    ZehnerSchluenderThermCond.C
+
+\*---------------------------------------------------------------------------*/
+
+
+#ifndef ZehnerSchluenderThermCond_H
+#define ZehnerSchluenderThermCond_H
+
+#include "thermCondModel.H"
+#include "scalarList.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                           Class ZehnerSchluenderThermCond Declaration
+\*---------------------------------------------------------------------------*/
+
+class ZehnerSchluenderThermCond
+:
+    public thermCondModel
+{
+  
+private:
+  
+    dictionary propsDict_;
+
+    mutable volScalarField partKsField_;
+  
+    word voidfractionFieldName_;
+  
+    const volScalarField& voidfraction_;
+    
+    scalarList typeKs_;
+
+    mutable double **partKs_;
+    
+    word wallQFactorName_;
+    
+    // ratio of half-cell-size and near-wall film
+    mutable volScalarField wallQFactor_;
+    
+    bool hasWallQFactor_;
+
+    void allocateMyArrays() const;
+
+    void calcPartKsField() const;
+
+public:
+
+    //- Runtime type information
+    TypeName("ZehnerSchluenderThermCond");
+
+
+    // Constructors
+
+        //- Construct from components
+        ZehnerSchluenderThermCond
+        (
+            const dictionary& dict,
+            cfdemCloud& sm
+        );
+
+    // Destructor
+
+        ~ZehnerSchluenderThermCond();
+
+
+    // Member Functions
+	
+	tmp<volScalarField> thermCond() const;
+	
+	tmp<volScalarField> thermDiff() const;
+       
+       
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/lagrangian/cfdemParticle/subModels/voidFractionModel/GaussVoidFraction/GaussVoidFraction.C

@@ -95,7 +95,7 @@ void GaussVoidFraction::setvoidFraction(double** const& mask,double**& voidfract
 
     scalar radius(-1);
     scalar volume(0);
-    scalar scaleVol= weight();
+    scalar scaleVol = weight();
     scalar scaleRadius = pow(porosity(),1./3.);
 
     for(int index=0; index< particleCloud_.numberOfParticles(); index++)
@@ -115,6 +115,7 @@ void GaussVoidFraction::setvoidFraction(double** const& mask,double**& voidfract
             label particleCenterCellID=particleCloud_.cellIDs()[index][0];
 
             radius = particleCloud_.radius(index);
+            if (multiWeights_) scaleVol = weight(index);
             volume = constant::mathematical::fourPiByThree*radius*radius*radius*scaleVol;
             radius *= scaleRadius;
 

src/lagrangian/cfdemParticle/subModels/voidFractionModel/bigParticleVoidFraction/bigParticleVoidFraction.C

@@ -94,7 +94,7 @@ void bigParticleVoidFraction::setvoidFraction(double** const& mask,double**& voi
 
     scalar radius(-1);
     scalar volume(0);
-    scalar scaleVol= weight();
+    scalar scaleVol = weight();
     scalar scaleRadius = pow(porosity(),1./3.);
 
     for(int index=0; index< particleCloud_.numberOfParticles(); index++)
@@ -113,6 +113,7 @@ void bigParticleVoidFraction::setvoidFraction(double** const& mask,double**& voi
             //collecting data
             label particleCenterCellID=particleCloud_.cellIDs()[index][0];
             radius =  particleCloud_.radius(index);
+            if (multiWeights_) scaleVol = weight(index);
             volume = constant::mathematical::fourPiByThree*radius*radius*radius*scaleVol;
             radius *= scaleRadius;
             vector positionCenter=particleCloud_.position(index);

src/lagrangian/cfdemParticle/subModels/voidFractionModel/centreVoidFraction/centreVoidFraction.C

@@ -85,7 +85,7 @@ void centreVoidFraction::setvoidFraction(double** const& mask,double**& voidfrac
     scalar radius(-1);
     scalar volume(0);
     scalar cellVol(0);
-    scalar scaleVol= weight();
+    scalar scaleVol = weight();
 
     for(int index=0; index< particleCloud_.numberOfParticles(); index++)
     {
@@ -99,6 +99,7 @@ void centreVoidFraction::setvoidFraction(double** const& mask,double**& voidfrac
 
             if (cellI >= 0)  // particel centre is in domain
             {
+                if (multiWeights_) scaleVol = weight(index);
                 cellVol = voidfractionNext_.mesh().V()[cellI];
                 radius = particleCloud_.radius(index);
                 volume = constant::mathematical::fourPiByThree*radius*radius*radius*scaleVol;

src/lagrangian/cfdemParticle/subModels/voidFractionModel/dividedVoidFraction/dividedVoidFraction.C

@@ -200,6 +200,7 @@ void dividedVoidFraction::setvoidFraction(double** const& mask,double**& voidfra
             position = particleCloud_.position(index);
             cellID = particleCloud_.cellIDs()[index][0];
             radius = particleCloud_.radius(index);
+            if (multiWeights_) scaleVol = weight(index);
             volume = Vp(index,radius,scaleVol);
             radius *= scaleRadius;
             cellVol = 0;

src/lagrangian/cfdemParticle/subModels/voidFractionModel/trilinearVoidFraction/trilinearVoidFraction.C

@@ -94,7 +94,7 @@ void trilinearVoidFraction::setvoidFraction(double** const& mask,double**& voidf
 
     scalar radius(-1.);
     scalar volume(0.);
-    const scalar scaleVol = weight();
+    scalar scaleVol = weight();
 
     vector partPos(0.,0.,0.);
     vector pt(0.,0.,0.);
@@ -140,6 +140,7 @@ void trilinearVoidFraction::setvoidFraction(double** const& mask,double**& voidf
         if (cellI >= 0)  // particel centre is in domain
         {
             radius = particleCloud_.radius(index);
+            if (multiWeights_) scaleVol = weight(index);
             volume = constant::mathematical::fourPiByThree * radius * radius * radius * scaleVol;
 
             // store volume for each particle

src/lagrangian/cfdemParticle/subModels/voidFractionModel/voidFractionModel/voidFractionModel.C

@@ -57,6 +57,7 @@ voidFractionModel::voidFractionModel
 :
     dict_(dict),
     particleCloud_(sm),
+    multiWeights_(false),
     voidfractionPrev_
     (   IOobject
         (
@@ -89,6 +90,7 @@ voidFractionModel::voidFractionModel
     porosity_(1.)
 {
     particleCloud_.dataExchangeM().allocateArray(cellsPerParticle_,1,1);
+    if (particleCloud_.getParticleEffVolFactors()) multiWeights_ = true;
 }
 
 

src/lagrangian/cfdemParticle/subModels/voidFractionModel/voidFractionModel/voidFractionModel.H

@@ -63,6 +63,8 @@ protected:
 
         cfdemCloud& particleCloud_;
 
+        bool multiWeights_;
+
         mutable volScalarField voidfractionPrev_;
 
         mutable volScalarField voidfractionNext_;
@@ -132,6 +134,11 @@ public:
 
         inline scalar weight()const { return weight_; }
 
+        inline scalar weight(label index) const
+        {
+            return particleCloud_.particleEffVolFactor(index);
+        }
+
         inline scalar porosity()const { return porosity_; }
 
         inline void checkWeightNporosity(dictionary& propsDict) const

src/lagrangian/cfdemParticleComp/Make/files

@@ -31,12 +31,13 @@ $(cfdTools)/newGlobal.C
 $(energyModels)/energyModel/energyModel.C
 $(energyModels)/energyModel/newEnergyModel.C
 $(energyModels)/heatTransferGunn/heatTransferGunn.C
-$(energyModels)/heatTransferGunnImplicit/heatTransferGunnImplicit.C
+$(energyModels)/heatTransferGunnPartField/heatTransferGunnPartField.C
 $(energyModels)/reactionHeat/reactionHeat.C
 
 $(thermCondModels)/thermCondModel/thermCondModel.C
 $(thermCondModels)/thermCondModel/newThermCondModel.C
 $(thermCondModels)/SyamlalThermCond/SyamlalThermCond.C
+$(thermCondModels)/ZehnerSchluenderThermCond/ZehnerSchluenderThermCond.C
 $(thermCondModels)/noTherm/noThermCond.C
 
 $(forceModels)/forceModel/forceModel.C
@@ -59,6 +60,7 @@ $(forceModels)/viscForce/viscForce.C
 $(forceModels)/MeiLift/MeiLift.C
 $(forceModels)/particleCellVolume/particleCellVolume.C
 $(forceModels)/BeetstraDrag/BeetstraDrag.C
+$(forceModels)/BeetstraDragPoly/BeetstraDragPoly.C
 $(forceModels)/dSauter/dSauter.C
 $(forceModels)/Fines/Fines.C
 $(forceModels)/Fines/FinesFields.C

tutorials/cfdemSolverPiso/ErgunTestCG/CFD/constant/couplingProperties

@@ -33,7 +33,7 @@ couplingInterval 100;
 
 voidFractionModel divided;//centre;//
 
-locateModel engine;//turboEngine;//
+locateModel engine;//turboEngineM2M;//
 
 meshMotionModel noMeshMotion;
 
@@ -49,7 +49,7 @@ averagingModel dense;//dilute;//
 
 clockModel off;//standardClock;//
 
-smoothingModel off;// constDiffSmoothing; //
+smoothingModel off;// localPSizeDiffSmoothing;// constDiffSmoothing; //
 
 forceModels
 (
@@ -61,6 +61,7 @@ forceModels
     GidaspowDrag
     //Archimedes
     //volWeightedAverage
+    //totalMomentumExchange
     particleCellVolume
 );
 
@@ -74,6 +75,14 @@ turbulenceModelType "turbulenceProperties";
 //===========================================================================//
 // sub-model properties
 
+localPSizeDiffSmoothingProps
+{
+    lowerLimit 0.1;
+    upperLimit 1e10;
+    dSmoothingLength 1.5e-3;
+    Csmoothing    1.0; 
+}
+
 constDiffSmoothingProps
 {
     lowerLimit 0.1;
@@ -123,7 +132,13 @@ volWeightedAverageProps
     lowerThreshold 0;
     verbose true;
 }
-
+totalMomentumExchangeProps
+{
+    implicitMomExFieldName "Ksl";
+    explicitMomExFieldName "none";
+    fluidVelFieldName "U";
+    granVelFieldName "Us";
+}
 GidaspowDragProps
 {
     verbose true;
@@ -147,12 +162,15 @@ KochHillDragProps
 BeetstraDragProps
 {
     velFieldName "U";
-    gravityFieldName "g";
-    rhoParticle     2000.;
     voidfractionFieldName "voidfraction";
+    granVelFieldName "Us";
     interpolation true;
-    useFilteredDragModel ;
+//    useFilteredDragModel;
-    useParcelSizeDependentFilteredDrag ;
+//    useParcelSizeDependentFilteredDrag;
+    g       9.81;
+    rhoP    7000.;
+    rho     10.;
+    nuf     1.5e-4;
     k       0.05;
     aLimit  0.0;
 //    verbose true;