diff --git a/applications/solvers/rcfdemSolverRhoSteadyPimpleChem/EEqn.H b/applications/solvers/rcfdemSolverRhoSteadyPimpleChem/EEqn.H index c524365d..27fe0cfd 100644 --- a/applications/solvers/rcfdemSolverRhoSteadyPimpleChem/EEqn.H +++ b/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); } diff --git a/applications/solvers/rcfdemSolverRhoSteadyPimpleChem/YEqn.H b/applications/solvers/rcfdemSolverRhoSteadyPimpleChem/YEqn.H index c5cbaca0..74b72b18 100644 --- a/applications/solvers/rcfdemSolverRhoSteadyPimpleChem/YEqn.H +++ b/applications/solvers/rcfdemSolverRhoSteadyPimpleChem/YEqn.H @@ -20,6 +20,7 @@ tmp > mvConvection #endif label inertIndex = -1; volScalarField Yt(0.0*Y[0]); + Sm *= 0.0; forAll(Y, i) { @@ -28,13 +29,17 @@ tmp > 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 > 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; } } diff --git a/applications/solvers/rcfdemSolverRhoSteadyPimpleChem/createFields.H b/applications/solvers/rcfdemSolverRhoSteadyPimpleChem/createFields.H index 2a96af04..af6c1ba1 100644 --- a/applications/solvers/rcfdemSolverRhoSteadyPimpleChem/createFields.H +++ b/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 ( diff --git a/applications/solvers/rcfdemSolverRhoSteadyPimpleChem/monitorMass.H b/applications/solvers/rcfdemSolverRhoSteadyPimpleChem/monitorMass.H deleted file mode 100644 index 40b6d36f..00000000 --- a/applications/solvers/rcfdemSolverRhoSteadyPimpleChem/monitorMass.H +++ /dev/null @@ -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); -} diff --git a/applications/solvers/rcfdemSolverRhoSteadyPimpleChem/monitorMassSinks.H b/applications/solvers/rcfdemSolverRhoSteadyPimpleChem/monitorMassSinks.H new file mode 100644 index 00000000..3ef2ac9a --- /dev/null +++ b/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()); + Info << Y[i].name() << ": mass sink = " << massSink << endl; +} diff --git a/applications/solvers/rcfdemSolverRhoSteadyPimpleChem/monitorMassSource.H b/applications/solvers/rcfdemSolverRhoSteadyPimpleChem/monitorMassSource.H new file mode 100644 index 00000000..3f907e79 --- /dev/null +++ b/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()); + Info << "total mass source strength = " << sourceStrength << endl; +} diff --git a/applications/solvers/rcfdemSolverRhoSteadyPimpleChem/monitorMassSources.H b/applications/solvers/rcfdemSolverRhoSteadyPimpleChem/monitorMassSources.H new file mode 100644 index 00000000..a2350462 --- /dev/null +++ b/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()); + Info << Y[i].name() << ": source strength = " << sourceStrength << endl; +} diff --git a/applications/solvers/rcfdemSolverRhoSteadyPimpleChem/pEqn.H b/applications/solvers/rcfdemSolverRhoSteadyPimpleChem/pEqn.H index b4439732..936a96c0 100644 --- a/applications/solvers/rcfdemSolverRhoSteadyPimpleChem/pEqn.H +++ b/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" diff --git a/applications/solvers/rcfdemSolverRhoSteadyPimpleChem/rcfdemSolverRhoSteadyPimpleChem.C b/applications/solvers/rcfdemSolverRhoSteadyPimpleChem/rcfdemSolverRhoSteadyPimpleChem.C index ef6420eb..9e5c0e02 100644 --- a/applications/solvers/rcfdemSolverRhoSteadyPimpleChem/rcfdemSolverRhoSteadyPimpleChem.C +++ b/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 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"); diff --git a/src/lagrangian/cfdemParticle/subModels/energyModel/reactionHeat/reactionHeat.C b/src/lagrangian/cfdemParticle/subModels/energyModel/reactionHeat/reactionHeat.C index 14f5c876..3ac05a94 100644 --- a/src/lagrangian/cfdemParticle/subModels/energyModel/reactionHeat/reactionHeat.C +++ b/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("interpolation",false)), verbose_(propsDict_.lookupOrDefault("verbose",false)), execution_(true), + smoothen_(propsDict_.lookupOrDefault("smoothen",false)), mesh_(sm.mesh()), maxSource_(1e30), reactionHeatName_(propsDict_.lookupOrDefault("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(); } diff --git a/src/lagrangian/cfdemParticle/subModels/energyModel/reactionHeat/reactionHeat.H b/src/lagrangian/cfdemParticle/subModels/energyModel/reactionHeat/reactionHeat.H index 73a4e2fa..db2f9ae9 100644 --- a/src/lagrangian/cfdemParticle/subModels/energyModel/reactionHeat/reactionHeat.H +++ b/src/lagrangian/cfdemParticle/subModels/energyModel/reactionHeat/reactionHeat.H @@ -52,6 +52,8 @@ protected: bool execution_; + bool smoothen_; + const fvMesh& mesh_; scalar maxSource_; diff --git a/src/lagrangian/cfdemParticle/subModels/smoothingModel/smoothingModel/smoothingModel.C b/src/lagrangian/cfdemParticle/subModels/smoothingModel/smoothingModel/smoothingModel.C index a98d4bde..61422602 100644 --- a/src/lagrangian/cfdemParticle/subModels/smoothingModel/smoothingModel/smoothingModel.C +++ b/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" ) {}