Minor changes in rcfdemSolverRhoSteadyPimpleChem wrt to implicit source treatment and in energyModel/reactionHeat to allow smoothing.

applications/solvers/rcfdemSolverRhoSteadyPimpleChem/EEqn.H

@@ -59,4 +59,6 @@
     thermo.correct();
 
     Info<< "T max/min : " << max(T).value() << " " << min(T).value() << endl;
+
+    QFluidCond = fvc::laplacian(voidfractionRec*thCond,T);
 }

applications/solvers/rcfdemSolverRhoSteadyPimpleChem/YEqn.H

@@ -20,6 +20,7 @@ tmp<fv::convectionScheme<scalar> > mvConvection
 #endif
     label inertIndex = -1;
     volScalarField Yt(0.0*Y[0]);
+    Sm *= 0.0;
 
     forAll(Y, i)
     {
@@ -28,13 +29,17 @@ tmp<fv::convectionScheme<scalar> > mvConvection
         {
             volScalarField& Yi = Y[i];
 
+            volScalarField sourceField(particleCloud.chemistryM(0).Smi(i));
+            volScalarField Smi0(neg(sourceField)*sourceField/(Yi + Yismall));
+            volScalarField Smi1(pos0(sourceField)*sourceField);
             fvScalarMatrix YiEqn
             (
                 mvConvection->fvmDiv(phi, Yi)
               - fvm::laplacian(voidfractionRec*turbulence->muEff(), Yi)
               ==
                 combustion->R(Yi)
-              + particleCloud.chemistryM(0).Smi(i)*p/p.prevIter()
+                + fvm::Sp(Smi0,Yi)
+                + Smi1
               + fvOptions(rho, Yi)
             );
 
@@ -48,8 +53,11 @@ tmp<fv::convectionScheme<scalar> > mvConvection
 
             fvOptions.correct(Yi);
 
+            #include "monitorMassSinks.H"
             Yi.max(0.0);
             if (Y[i].name() != inertSpecie) Yt += Yi;
+            #include "monitorMassSources.H"
+            Sm += Smi0*Yi+Smi1;
         }
     }
 

applications/solvers/rcfdemSolverRhoSteadyPimpleChem/createFields.H

@@ -133,6 +133,20 @@ Info<< "Reading thermophysical properties\n" << endl;
         dimensionedScalar("zero", dimensionSet(1,-1,-3,0,0,0,0), 0.0)
     );
 
+    volScalarField Sm
+    (
+        IOobject
+        (
+            "Sm",
+            runTime.timeName(),
+            mesh,
+            IOobject::NO_READ,
+            IOobject::NO_WRITE
+         ),
+        mesh,
+        dimensionedScalar("zero",dimMass/(dimVol*dimTime),0.0)
+    );
+
     Info<< "\nCreating fluid-particle heat flux field\n" << endl;
     volScalarField Qsource
     (

applications/solvers/rcfdemSolverRhoSteadyPimpleChem/monitorMass.H

@@ -1,9 +0,0 @@
-{
-    m=gSum(rhoeps*1.0*rhoeps.mesh().V());
-    if(counter==0) m0=m;
-    counter++;
-    Info << "\ncurrent gas mass = " << m << "\n" << endl;
-    Info << "\ncurrent added gas mass = " << m-m0 << "\n" << endl;
-
-    QFluidCond = fvc::laplacian(voidfractionRec*thCond,T);
-}

applications/solvers/rcfdemSolverRhoSteadyPimpleChem/monitorMassSinks.H

@@ -0,0 +1,12 @@
+{
+    scalar massSink = 0.0;
+    forAll(Yi,cellI)
+    {
+        if (Yi[cellI] <= 0.0)
+        {
+            massSink += rhoeps[cellI]*Yi[cellI]*Yi.mesh().V()[cellI];
+        }
+    }
+    reduce(massSink, sumOp<scalar>());
+    Info << Y[i].name() << ": mass sink = " << massSink << endl;
+}

applications/solvers/rcfdemSolverRhoSteadyPimpleChem/monitorMassSource.H

@@ -0,0 +1,9 @@
+{
+    scalar sourceStrength = 0.0;
+    forAll(p,cellI)
+    {
+        sourceStrength += (Sm0[cellI]*p[cellI]+Sm1[cellI])*p.mesh().V()[cellI];
+    }
+    reduce(sourceStrength, sumOp<scalar>());
+    Info << "total mass source strength = " << sourceStrength << endl;
+}

applications/solvers/rcfdemSolverRhoSteadyPimpleChem/monitorMassSources.H

@@ -0,0 +1,9 @@
+{
+    scalar sourceStrength = 0.0;
+    forAll(Yi,cellI)
+    {
+        sourceStrength += (Smi0[cellI]*Yi[cellI]+Smi1[cellI])*Yi.mesh().V()[cellI];
+    }
+    reduce(sourceStrength, sumOp<scalar>());
+    Info << Y[i].name() << ": source strength = " << sourceStrength << endl;
+}

applications/solvers/rcfdemSolverRhoSteadyPimpleChem/pEqn.H

@@ -37,7 +37,8 @@ else
     // Update the pressure BCs to ensure flux consistency
     constrainPressure(p, rhoeps, U, phi, rhorAUf);
 
-    volScalarField SmbyP(particleCloud.chemistryM(0).Sm() / p);
+    volScalarField Sm0(neg(Sm)*Sm/(p + psmall));
+    volScalarField Sm1(pos0(Sm)*Sm);
 
     while (pimple.correctNonOrthogonal())
     {
@@ -47,7 +48,8 @@ else
             fvc::div(phi)
           - fvm::laplacian(rhorAUf, p)
           ==
-            fvm::Sp(SmbyP, p)
+            fvm::Sp(Sm0,p)
+          + Sm1
           + fvOptions(psi, p, rho.name())
         );
 
@@ -60,6 +62,7 @@ else
             phi += pEqn.flux();
         }
     }
+    #include "monitorMassSource.H"
 }
 
 #include "rhoEqn.H"

applications/solvers/rcfdemSolverRhoSteadyPimpleChem/rcfdemSolverRhoSteadyPimpleChem.C

@@ -76,16 +76,12 @@ int main(int argc, char *argv[])
     #include "createFieldRefs.H"
     #include "createFvOptions.H"
 
-    // create cfdemCloud
-    //#include "readGravitationalAcceleration.H"
     cfdemCloudRec<cfdemCloudEnergy> particleCloud(mesh);
     #include "checkModelType.H"
     recBase recurrenceBase(mesh);
     #include "updateFields.H"
 
     turbulence->validate();
-    //#include "compressibleCourantNo.H"
-    //#include "setInitialDeltaT.H"
 
     // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -104,6 +100,8 @@ int main(int argc, char *argv[])
     scalar m0(0.0);
     label counter(0);
     p.storePrevIter();
+    const dimensionedScalar psmall("psmall", dimPressure, small);
+    const dimensionedScalar Yismall("Yismall", dimless, small);
 
     while (runTime.run())
     {
@@ -121,9 +119,6 @@ int main(int argc, char *argv[])
         particleCloud.clockM().start(2,"Coupling");
         bool hasEvolved = particleCloud.evolve(voidfraction,Us,U);
 
-        //voidfraction = voidfractionRec;
-        //Us = UsRec;
-
         if(hasEvolved)
         {
             particleCloud.smoothingM().smoothen(particleCloud.forceM(0).impParticleForces());
@@ -160,8 +155,8 @@ int main(int argc, char *argv[])
 
             // --- Pressure-velocity PIMPLE corrector loop
 
-
             #include "UEqn.H"
+            #include "YEqn.H"
             #include "EEqn.H"
 
             // --- Pressure corrector loop
@@ -172,7 +167,6 @@ int main(int argc, char *argv[])
                 #include "pEqn.H"
                 rhoeps=rho*voidfractionRec;
             }
-            #include "YEqn.H"
 
             if (pimple.turbCorr())
             {
@@ -180,8 +174,6 @@ int main(int argc, char *argv[])
             }
         }
 
-        #include "monitorMass.H"
-
         totalStepCounter++;
         particleCloud.clockM().stop("Flow");
 

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

@@ -22,6 +22,7 @@ License
 #include "reactionHeat.H"
 #include "addToRunTimeSelectionTable.H"
 #include "dataExchangeModel.H"
+#include "smoothingModel.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -48,6 +49,7 @@ reactionHeat::reactionHeat
     interpolation_(propsDict_.lookupOrDefault<bool>("interpolation",false)),
     verbose_(propsDict_.lookupOrDefault<bool>("verbose",false)),
     execution_(true),
+    smoothen_(propsDict_.lookupOrDefault<bool>("smoothen",false)),
     mesh_(sm.mesh()),
     maxSource_(1e30),
     reactionHeatName_(propsDict_.lookupOrDefault<word>("reactionHeatName","reactionHeat")),
@@ -123,6 +125,10 @@ void reactionHeat::calcEnergyContribution()
     );
 
     reactionHeatField_.primitiveFieldRef() /= (reactionHeatField_.mesh().V() * Nevery_ * couplingTimestep_);
+    if (smoothen_)
+    {
+        particleCloud_.smoothingM().smoothen(reactionHeatField_);
+    }
 
     forAll(reactionHeatField_,cellI)
     {
@@ -134,11 +140,8 @@ void reactionHeat::calcEnergyContribution()
         }
     }
 
-    if (verbose_)
+    Info << "reaction heat per unit time = "
-    {
+                 << gSum(reactionHeatField_*1.0*reactionHeatField_.mesh().V()) << endl;
-        Info << "reaction heat per unit time = "
-                     << gSum(reactionHeatField_*1.0*reactionHeatField_.mesh().V()) << endl;
-    }
 
     reactionHeatField_.correctBoundaryConditions();
 }

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

@@ -52,6 +52,8 @@ protected:
 
     bool execution_;
 
+    bool smoothen_;
+
     const fvMesh& mesh_;
 
     scalar maxSource_;

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

@@ -69,7 +69,8 @@ smoothingModel::smoothingModel
             IOobject::NO_WRITE
         ),
         particleCloud_.mesh(),
-        dimensionedVector("zero", dimensionSet(0,0,0,0,0), vector::zero)
+        dimensionedVector("zero", dimensionSet(0,0,0,0,0), vector::zero),
+        "zeroGradient"
     ),
     sSmoothField_
     (
@@ -82,7 +83,8 @@ smoothingModel::smoothingModel
             IOobject::NO_WRITE
         ),
         particleCloud_.mesh(),
-        dimensionedScalar("zero", dimensionSet(0,0,0,0,0), 0)
+        dimensionedScalar("zero", dimensionSet(0,0,0,0,0), 0),
+        "zeroGradient"
     )
 {}