Specific heat and thermal conductivities for various species.

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

@@ -76,14 +76,15 @@ heatTransferGunn::heatTransferGunn
     partTempField_
     (   IOobject
         (
-            "particleTemp",
+            "partTemp",
             sm.mesh().time().timeName(),
             sm.mesh(),
             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

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

@@ -59,6 +59,7 @@ heatTransferGunnPartField::heatTransferGunnPartField
 	"zeroGradient"
     ),
     partRhoField_(sm.mesh().lookupObject<volScalarField>("partRho")),
+    typeCp_(propsDict_.lookupOrDefault<scalarList>("specificHeatCapacities",scalarList(1,-1.0))),
     partCp_(NULL),
     thermCondModel_
     (
@@ -70,11 +71,16 @@ heatTransferGunnPartField::heatTransferGunnPartField
     ),
     fvOptions(fv::options::New(sm.mesh()))
 {
-    if(!implicit_)
+    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);
+    }
+
     allocateMyArrays();
 }
 
@@ -102,6 +108,19 @@ void heatTransferGunnPartField::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
@@ -155,19 +174,20 @@ void heatTransferGunnPartField::postFlow()
 void heatTransferGunnPartField::solve()
 {
     Info << "patch types of partTemp boundary: " << partTempField_.boundaryField().types() << endl;
-    volScalarField Qsource = QPartFluidCoeff_*tempField_/partCpField_;
+    volScalarField Qsource = QPartFluidCoeff_*tempField_;
     volScalarField alphaP = 1.0 - voidfraction_;
-    volScalarField partRhoEff = alphaP*partRhoField_;
+    volScalarField partCpEff = alphaP*partRhoField_*partCpField_;
-//    volScalarField thCond = thermCondModel_().thermCond();
+    volScalarField thCondEff = alphaP*thermCondModel_().thermCond();
+//    thCondEff.correctBoundaryConditions();
 
 
     fvScalarMatrix partTEqn
     (
         Qsource
-      - fvm::Sp(QPartFluidCoeff_/partCpField_, partTempField_)
+      - fvm::Sp(QPartFluidCoeff_, partTempField_)
- //     - fvm::laplacian(alphaP*thCond/partCpField_,partTempField_)
+      - fvm::laplacian(thCondEff,partTempField_)
      ==
-        fvOptions(partRhoEff, partTempField_)
+        fvOptions(partCpEff, partTempField_)
     );
   // if transient add time derivative - need particle density and specific heat fields
   // if sources activated add sources

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

@@ -29,6 +29,7 @@ License
 #include "cfdemCloudEnergy.H"
 #include "heatTransferGunn.H"
 #include "fvOptions.H"
+#include "scalarList.H"
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 namespace Foam
@@ -48,6 +49,8 @@ private:
 
     const volScalarField& partRhoField_;
 
+    scalarList typeCp_;
+
     mutable double **partCp_;
 
     autoPtr<thermCondModel> thermCondModel_;

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

@@ -50,9 +50,24 @@ ZehnerSchluenderThermCond::ZehnerSchluenderThermCond
 :
     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_)),
-    ks0_(transportProperties_.lookup("ks")),
+    typeKs_(propsDict_.lookupOrDefault<scalarList>("thermalConductivities",scalarList(1,-1.0))),
+    partKs_(NULL),
     wallQFactorName_(propsDict_.lookupOrDefault<word>("wallQFactorName","wallQFactor")),
     wallQFactor_
     (   IOobject
@@ -68,6 +83,17 @@ ZehnerSchluenderThermCond::ZehnerSchluenderThermCond
     ),
     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;
@@ -79,9 +105,17 @@ ZehnerSchluenderThermCond::ZehnerSchluenderThermCond
 // * * * * * * * * * * * * * * * * 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
@@ -106,7 +140,8 @@ tmp<volScalarField> ZehnerSchluenderThermCond::thermCond() const
 
     volScalarField& svf = tvf.ref();
 
-    scalar A = ks0_.value()/kf0_.value();
+    calcPartKsField();
+    scalar A = 0.0;
     scalar B = 0.0;
     scalar C = 0.0;
     scalar k = 0.0;
@@ -116,12 +151,11 @@ tmp<volScalarField> ZehnerSchluenderThermCond::thermCond() const
 
     forAll(svf, cellI)
     {
-        // debugging
-        Pout << "calculating field in cell " << cellI << endl;
         voidfraction = voidfraction_[cellI];
-        if(voidfraction > 1.0 - SMALL) svf[cellI] = 0.0;
+        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);
             OnemBoA = 1.0 - B/A;
             C = (A - 1) / (OnemBoA * OnemBoA) * B/A * log(A/B) - (B - 1)/OnemBoA - 0.5 * (B + 1);
@@ -131,8 +165,6 @@ tmp<volScalarField> ZehnerSchluenderThermCond::thermCond() const
         }
     }
 
-    // debugging
-    Pout << "patch types of svf boundary: " << svf.boundaryField().types() << endl;
     svf.correctBoundaryConditions();
 
     // if a wallQFactor field is present, use it to scale heat transport through a patch
@@ -151,8 +183,35 @@ tmp<volScalarField> ZehnerSchluenderThermCond::thermDiff() const
     FatalError << "ZehnerSchluenderThermCond does not provide thermal diffusivity." << abort(FatalError);
 }
 
+void ZehnerSchluenderThermCond::calcPartKsField() const
+{
+    if (typeKs_.size() <= 1) return;
 
+    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

@@ -34,6 +34,7 @@ SourceFiles
 #define ZehnerSchluenderThermCond_H
 
 #include "thermCondModel.H"
+#include "scalarList.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -53,11 +54,15 @@ private:
   
     dictionary propsDict_;
 
+    mutable volScalarField partKsField_;
+  
     word voidfractionFieldName_;
   
     const volScalarField& voidfraction_;
     
-    dimensionedScalar ks0_;
+    scalarList typeKs_;
+
+    mutable double **partKs_;
     
     word wallQFactorName_;
     
@@ -66,6 +71,10 @@ private:
     
     bool hasWallQFactor_;
 
+    void allocateMyArrays() const;
+
+    void calcPartKsField() const;
+
 public:
 
     //- Runtime type information