Minor modifications for rcfdemSolverCoupledHeattransfer to improve stability.

applications/solvers/rcfdemSolverCoupledHeattransfer/TEqn.H

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

applications/solvers/rcfdemSolverCoupledHeattransfer/createFields.H

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

applications/solvers/rcfdemSolverCoupledHeattransfer/rcfdemSolverCoupledHeattransfer.C

@@ -98,7 +98,7 @@ int main(int argc, char *argv[])
 
         particleCloud.clockM().start(26,"Flow");
         #include "updateRho.H"
-        #include "TEqImp.H"
+        #include "TEqn.H"
         particleCloud.clockM().stop("Flow");
 
         stepCounter++;

applications/solvers/rcfdemSolverCoupledHeattransfer/updateRho.H

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