diff --git a/applications/solvers/cfdemSolverRhoPimple/EEqn.H b/applications/solvers/cfdemSolverRhoPimple/EEqn.H index 606d9fe2..0f1e6aff 100644 --- a/applications/solvers/cfdemSolverRhoPimple/EEqn.H +++ b/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.postFlow(); - particleCloud.clockM().stop("postFlow"); + particleCloud.clockM().start(31,"energySolve"); + particleCloud.solve(); + particleCloud.clockM().stop("energySolve"); + } diff --git a/applications/solvers/cfdemSolverRhoPimple/cfdemSolverRhoPimple.C b/applications/solvers/cfdemSolverRhoPimple/cfdemSolverRhoPimple.C index b9cf2fcb..39fa85cc 100644 --- a/applications/solvers/cfdemSolverRhoPimple/cfdemSolverRhoPimple.C +++ b/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(); diff --git a/applications/solvers/cfdemSolverRhoSimple/EEqn.H b/applications/solvers/cfdemSolverRhoSimple/EEqn.H new file mode 100644 index 00000000..3dbb5c3f --- /dev/null +++ b/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"); +} diff --git a/applications/solvers/cfdemSolverRhoSimple/Make/files b/applications/solvers/cfdemSolverRhoSimple/Make/files new file mode 100644 index 00000000..d9db0744 --- /dev/null +++ b/applications/solvers/cfdemSolverRhoSimple/Make/files @@ -0,0 +1,3 @@ +cfdemSolverRhoSimple.C + +EXE=$(CFDEM_APP_DIR)/cfdemSolverRhoSimple diff --git a/applications/solvers/cfdemSolverRhoSimple/Make/options b/applications/solvers/cfdemSolverRhoSimple/Make/options new file mode 100644 index 00000000..0377ece5 --- /dev/null +++ b/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) diff --git a/applications/solvers/cfdemSolverRhoSimple/UEqn.H b/applications/solvers/cfdemSolverRhoSimple/UEqn.H new file mode 100644 index 00000000..8fbd30f2 --- /dev/null +++ b/applications/solvers/cfdemSolverRhoSimple/UEqn.H @@ -0,0 +1,30 @@ +// Solve the Momentum equation +particleCloud.otherForces(fOther); + +tmp 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); diff --git a/applications/solvers/cfdemSolverRhoSimple/cfdemSolverRhoSimple.C b/applications/solvers/cfdemSolverRhoSimple/cfdemSolverRhoSimple.C new file mode 100644 index 00000000..4e1821fc --- /dev/null +++ b/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 . + + 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()); + 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; +} + + +// ************************************************************************* // diff --git a/applications/solvers/cfdemSolverRhoSimple/createFieldRefs.H b/applications/solvers/cfdemSolverRhoSimple/createFieldRefs.H new file mode 100644 index 00000000..5842906a --- /dev/null +++ b/applications/solvers/cfdemSolverRhoSimple/createFieldRefs.H @@ -0,0 +1,2 @@ +const volScalarField& T = thermo.T(); +const volScalarField& psi = thermo.psi(); diff --git a/applications/solvers/cfdemSolverRhoSimple/createFields.H b/applications/solvers/cfdemSolverRhoSimple/createFields.H new file mode 100644 index 00000000..a9225229 --- /dev/null +++ b/applications/solvers/cfdemSolverRhoSimple/createFields.H @@ -0,0 +1,242 @@ +Info<< "Reading thermophysical properties\n" << endl; + + autoPtr 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 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 + ); + +//=============================== diff --git a/applications/solvers/cfdemSolverRhoSimple/pEqn.H b/applications/solvers/cfdemSolverRhoSimple/pEqn.H new file mode 100644 index 00000000..5b2dbcbd --- /dev/null +++ b/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); \ No newline at end of file diff --git a/etc/solver-list.txt b/etc/solver-list.txt index 9d9911fa..60c0eb39 100644 --- a/etc/solver-list.txt +++ b/etc/solver-list.txt @@ -1,5 +1,6 @@ cfdemSolverPisoMS/dir cfdemSolverPiso/dir cfdemSolverRhoPimple/dir +cfdemSolverRhoSimple/dir cfdemSolverIB/dir cfdemSolverPisoScalar/dir diff --git a/src/lagrangian/cfdemParticle/Make/files b/src/lagrangian/cfdemParticle/Make/files index 7cc1c491..93e0999c 100644 --- a/src/lagrangian/cfdemParticle/Make/files +++ b/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 diff --git a/src/lagrangian/cfdemParticle/cfdemCloud/cfdemCloud.C b/src/lagrangian/cfdemParticle/cfdemCloud/cfdemCloud.C index 89e2cead..9f5d8955 100644 --- a/src/lagrangian/cfdemParticle/cfdemCloud/cfdemCloud.C +++ b/src/lagrangian/cfdemParticle/cfdemCloud/cfdemCloud.C @@ -83,6 +83,9 @@ cfdemCloud::cfdemCloud ignore_(false), allowCFDsubTimestep_(true), limitDEMForces_(false), + getParticleDensities_(couplingProperties_.lookupOrDefault("getParticleDensities",false)), + getParticleEffVolFactors_(couplingProperties_.lookupOrDefault("getParticleEffVolFactors",false)), + getParticleTypes_(couplingProperties_.lookupOrDefault("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; } diff --git a/src/lagrangian/cfdemParticle/cfdemCloud/cfdemCloud.H b/src/lagrangian/cfdemParticle/cfdemCloud/cfdemCloud.H index 0f74cc1f..6f46acd1 100644 --- a/src/lagrangian/cfdemParticle/cfdemCloud/cfdemCloud.H +++ b/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;} - virtual label particleType(label index) const {return -1;} + + inline bool getParticleDensities() const; + 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;} diff --git a/src/lagrangian/cfdemParticle/cfdemCloud/cfdemCloudI.H b/src/lagrangian/cfdemParticle/cfdemCloud/cfdemCloudI.H index 2fbe1744..02b3a7dc 100644 --- a/src/lagrangian/cfdemParticle/cfdemCloud/cfdemCloudI.H +++ b/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_; diff --git a/src/lagrangian/cfdemParticle/derived/cfdemCloudEnergy/cfdemCloudEnergy.C b/src/lagrangian/cfdemParticle/derived/cfdemCloudEnergy/cfdemCloudEnergy.C index 1ed2727b..4a2f44d3 100644 --- a/src/lagrangian/cfdemParticle/derived/cfdemCloudEnergy/cfdemCloudEnergy.C +++ b/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("expNusselt",false)), interpolation_(propsDict_.lookupOrDefault("interpolation",false)), verbose_(propsDict_.lookupOrDefault("verbose",false)), + implicit_(propsDict_.lookupOrDefault("implicit",true)), QPartFluidName_(propsDict_.lookupOrDefault("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("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("calcPartTempField",false)), + calcPartTempAve_(propsDict_.lookupOrDefault("calcPartTempAve",false)), + partTempAve_(0.0), tempFieldName_(propsDict_.lookupOrDefault("tempFieldName","T")), tempField_(sm.mesh().lookupObject (tempFieldName_)), voidfractionFieldName_(propsDict_.lookupOrDefault("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 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()); + 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 - forAll(QPartFluid_,cellI) + // limit source term in explicit treatment + if(!implicit_) { - scalar EuFieldInCell = QPartFluid_[cellI]; - - if(mag(EuFieldInCell) > maxSource_ ) + forAll(QPartFluid_,cellI) { - Pout << "limiting source term\n" << endl ; - QPartFluid_[cellI] = sign(EuFieldInCell) * maxSource_; + scalar EuFieldInCell = QPartFluid_[cellI]; + + 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]); -} - -void heatTransferGunn::heatFluxCoeff(label index, scalar h, scalar As) -{ - //no heat transfer coefficient in explicit model -} - -void heatTransferGunn::giveData(int call) -{ - if(call == 0) + scalar hAs = h * As; + partHeatFlux_[index][0] = - hAs * partTemp_[index][0]; + if(!implicit_) { - Info << "total convective particle-fluid heat flux [W] (Eulerian) = " << gSum(QPartFluid_*1.0*QPartFluid_.mesh().V()) << endl; - - particleCloud_.dataExchangeM().giveData(partHeatFluxName_,"scalar-atom", partHeatFlux_); + partHeatFlux_[index][0] += hAs * Tfluid; } + 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 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_); } // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // diff --git a/src/lagrangian/cfdemParticle/subModels/energyModel/heatTransferGunn/heatTransferGunn.H b/src/lagrangian/cfdemParticle/subModels/energyModel/heatTransferGunn/heatTransferGunn.H index 4862e5ed..aa6facf7 100644 --- a/src/lagrangian/cfdemParticle/subModels/energyModel/heatTransferGunn/heatTransferGunn.H +++ b/src/lagrangian/cfdemParticle/subModels/energyModel/heatTransferGunn/heatTransferGunn.H @@ -15,7 +15,7 @@ License along with this code. If not, see . 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; }; diff --git a/src/lagrangian/cfdemParticle/subModels/energyModel/heatTransferGunnPartField/heatTransferGunnPartField.C b/src/lagrangian/cfdemParticle/subModels/energyModel/heatTransferGunnPartField/heatTransferGunnPartField.C new file mode 100644 index 00000000..2d8af0e6 --- /dev/null +++ b/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 . + + 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("partRho")), + typeCp_(propsDict_.lookupOrDefault("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 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 + +// ************************************************************************* // diff --git a/src/lagrangian/cfdemParticle/subModels/energyModel/heatTransferGunnPartField/heatTransferGunnPartField.H b/src/lagrangian/cfdemParticle/subModels/energyModel/heatTransferGunnPartField/heatTransferGunnPartField.H new file mode 100644 index 00000000..296c9a49 --- /dev/null +++ b/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 . + + 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_; + + 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 + +// ************************************************************************* // diff --git a/src/lagrangian/cfdemParticle/subModels/forceModel/BeetstraDrag/BeetstraDrag.C b/src/lagrangian/cfdemParticle/subModels/forceModel/BeetstraDrag/BeetstraDrag.C index b37eaba7..59c39254 100644 --- a/src/lagrangian/cfdemParticle/subModels/forceModel/BeetstraDrag/BeetstraDrag.C +++ b/src/lagrangian/cfdemParticle/subModels/forceModel/BeetstraDrag/BeetstraDrag.C @@ -12,6 +12,7 @@ License along with this code. If not, see . 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 (velFieldName_)), voidfractionFieldName_(propsDict_.lookup("voidfractionFieldName")), voidfraction_(sm.mesh().lookupObject (voidfractionFieldName_)), + minVoidfraction_(propsDict_.lookupOrDefault("minVoidfraction",0.1)), UsFieldName_(propsDict_.lookup("granVelFieldName")), UsField_(sm.mesh().lookupObject (UsFieldName_)), scaleDia_(1.), - scaleDrag_(1.) + typeCG_(propsDict_.lookupOrDefault("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("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<0.10) voidfraction = 0.10; + if (voidfraction > 1.00) voidfraction = 1.0; + 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; - dragCoefficient = 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))); + Rep=dSauter*voidfraction*magUr/nuf + SMALL; + + dragCoefficient = F(voidfraction, Rep) + *3*M_PI*nuf*rho*voidfraction + *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< 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 diff --git a/src/lagrangian/cfdemParticle/subModels/forceModel/BeetstraDrag/BeetstraDrag.H b/src/lagrangian/cfdemParticle/subModels/forceModel/BeetstraDrag/BeetstraDrag.H index f6e079ab..c2128855 100644 --- a/src/lagrangian/cfdemParticle/subModels/forceModel/BeetstraDrag/BeetstraDrag.H +++ b/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 diff --git a/src/lagrangian/cfdemParticle/subModels/forceModel/BeetstraDragPoly/BeetstraDragPoly.C b/src/lagrangian/cfdemParticle/subModels/forceModel/BeetstraDragPoly/BeetstraDragPoly.C new file mode 100644 index 00000000..0d886657 --- /dev/null +++ b/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 . + + 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("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(sm.mesh().lookupObject ("alphaP"))); + alphaP_.set(&alphaP); + volScalarField& alphaSt(const_cast(sm.mesh().lookupObject ("alphaSt"))); + alphaSt_.set(&alphaSt); + volScalarField& dSauter(const_cast(sm.mesh().lookupObject ("dSauterMix"))); + dSauter_.set(&dSauter); + } + else + { + volScalarField& dSauter(const_cast(sm.mesh().lookupObject ("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 + +// ************************************************************************* // diff --git a/src/lagrangian/cfdemParticle/subModels/forceModel/BeetstraDragPoly/BeetstraDragPoly.H b/src/lagrangian/cfdemParticle/subModels/forceModel/BeetstraDragPoly/BeetstraDragPoly.H new file mode 100644 index 00000000..88965010 --- /dev/null +++ b/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 . + + 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 alphaP_; + + autoPtr alphaSt_; + + autoPtr 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 + +// ************************************************************************* // diff --git a/src/lagrangian/cfdemParticle/subModels/forceModel/dSauter/dSauter.C b/src/lagrangian/cfdemParticle/subModels/forceModel/dSauter/dSauter.C index 081742c3..9f1a555b 100644 --- a/src/lagrangian/cfdemParticle/subModels/forceModel/dSauter/dSauter.C +++ b/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), - scaleDiaDist_(1.0), + typeCG_(propsDict_.lookupOrDefault("coarseGrainingFactors",scalarList(1,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; - } - else if (particleCloud_.cg() > 1) - { - scaleDia_ = particleCloud_.cg(); - Info << "dSauter using scaleCG from liggghts cg = " << scaleDia_ << endl; + Info << "dSauter using CG factor(s) = " << typeCG_ << endl; } allocateMyArrays(); - label cellI=0; - scalar ds(0); - scalar scale = scaleDiaDist_/scaleDia_; + label cellI = 0; + label partType = 1; + 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; - d3_[index][0] = ds*ds*ds; + d2_[index][0] = ds*ds*effVolFac*cg; + 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 { diff --git a/src/lagrangian/cfdemParticle/subModels/forceModel/dSauter/dSauter.H b/src/lagrangian/cfdemParticle/subModels/forceModel/dSauter/dSauter.H index 7a81fab5..f1e1ca6e 100644 --- a/src/lagrangian/cfdemParticle/subModels/forceModel/dSauter/dSauter.H +++ b/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 - - mutable scalar scaleDiaDist_; // scaling due to modified particle size distribution + scalarList typeCG_; mutable volScalarField d2Field_; diff --git a/src/lagrangian/cfdemParticle/subModels/thermCondModel/ZehnerSchluenderThermCond/ZehnerSchluenderThermCond.C b/src/lagrangian/cfdemParticle/subModels/thermCondModel/ZehnerSchluenderThermCond/ZehnerSchluenderThermCond.C new file mode 100644 index 00000000..7e4ca613 --- /dev/null +++ b/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 . + + 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("voidfractionFieldName","voidfraction")), + voidfraction_(sm.mesh().lookupObject (voidfractionFieldName_)), + typeKs_(propsDict_.lookupOrDefault("thermalConductivities",scalarList(1,-1.0))), + partKs_(NULL), + wallQFactorName_(propsDict_.lookupOrDefault("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 ZehnerSchluenderThermCond::thermCond() const +{ + tmp 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 ZehnerSchluenderThermCond::thermDiff() const +{ + FatalError << "ZehnerSchluenderThermCond does not provide thermal diffusivity." << abort(FatalError); + return tmp(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 + +// ************************************************************************* // diff --git a/src/lagrangian/cfdemParticle/subModels/thermCondModel/ZehnerSchluenderThermCond/ZehnerSchluenderThermCond.H b/src/lagrangian/cfdemParticle/subModels/thermCondModel/ZehnerSchluenderThermCond/ZehnerSchluenderThermCond.H new file mode 100644 index 00000000..b3cafd7c --- /dev/null +++ b/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 . + + 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 thermCond() const; + + tmp thermDiff() const; + + +}; + + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +} // End namespace Foam + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#endif + +// ************************************************************************* // diff --git a/src/lagrangian/cfdemParticle/subModels/voidFractionModel/GaussVoidFraction/GaussVoidFraction.C b/src/lagrangian/cfdemParticle/subModels/voidFractionModel/GaussVoidFraction/GaussVoidFraction.C index 3964b074..3db3e053 100644 --- a/src/lagrangian/cfdemParticle/subModels/voidFractionModel/GaussVoidFraction/GaussVoidFraction.C +++ b/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; diff --git a/src/lagrangian/cfdemParticle/subModels/voidFractionModel/bigParticleVoidFraction/bigParticleVoidFraction.C b/src/lagrangian/cfdemParticle/subModels/voidFractionModel/bigParticleVoidFraction/bigParticleVoidFraction.C index d176668f..e7dc28bc 100644 --- a/src/lagrangian/cfdemParticle/subModels/voidFractionModel/bigParticleVoidFraction/bigParticleVoidFraction.C +++ b/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); diff --git a/src/lagrangian/cfdemParticle/subModels/voidFractionModel/centreVoidFraction/centreVoidFraction.C b/src/lagrangian/cfdemParticle/subModels/voidFractionModel/centreVoidFraction/centreVoidFraction.C index 0f1fb2bf..6aa96bce 100644 --- a/src/lagrangian/cfdemParticle/subModels/voidFractionModel/centreVoidFraction/centreVoidFraction.C +++ b/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; diff --git a/src/lagrangian/cfdemParticle/subModels/voidFractionModel/dividedVoidFraction/dividedVoidFraction.C b/src/lagrangian/cfdemParticle/subModels/voidFractionModel/dividedVoidFraction/dividedVoidFraction.C index 3066be22..6648d3c7 100644 --- a/src/lagrangian/cfdemParticle/subModels/voidFractionModel/dividedVoidFraction/dividedVoidFraction.C +++ b/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; diff --git a/src/lagrangian/cfdemParticle/subModels/voidFractionModel/trilinearVoidFraction/trilinearVoidFraction.C b/src/lagrangian/cfdemParticle/subModels/voidFractionModel/trilinearVoidFraction/trilinearVoidFraction.C index 5e395e31..8413532a 100644 --- a/src/lagrangian/cfdemParticle/subModels/voidFractionModel/trilinearVoidFraction/trilinearVoidFraction.C +++ b/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 diff --git a/src/lagrangian/cfdemParticle/subModels/voidFractionModel/voidFractionModel/voidFractionModel.C b/src/lagrangian/cfdemParticle/subModels/voidFractionModel/voidFractionModel/voidFractionModel.C index 3368fc7a..424b8eef 100644 --- a/src/lagrangian/cfdemParticle/subModels/voidFractionModel/voidFractionModel/voidFractionModel.C +++ b/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; } diff --git a/src/lagrangian/cfdemParticle/subModels/voidFractionModel/voidFractionModel/voidFractionModel.H b/src/lagrangian/cfdemParticle/subModels/voidFractionModel/voidFractionModel/voidFractionModel.H index 877b49bf..9a760b31 100644 --- a/src/lagrangian/cfdemParticle/subModels/voidFractionModel/voidFractionModel/voidFractionModel.H +++ b/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 diff --git a/src/lagrangian/cfdemParticleComp/Make/files b/src/lagrangian/cfdemParticleComp/Make/files index b119a3b3..c586a4fb 100644 --- a/src/lagrangian/cfdemParticleComp/Make/files +++ b/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 diff --git a/tutorials/cfdemSolverPiso/ErgunTestCG/CFD/constant/couplingProperties b/tutorials/cfdemSolverPiso/ErgunTestCG/CFD/constant/couplingProperties index fd901353..1216be50 100644 --- a/tutorials/cfdemSolverPiso/ErgunTestCG/CFD/constant/couplingProperties +++ b/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 ; - useParcelSizeDependentFilteredDrag ; +// useFilteredDragModel; +// useParcelSizeDependentFilteredDrag; + g 9.81; + rhoP 7000.; + rho 10.; + nuf 1.5e-4; k 0.05; aLimit 0.0; // verbose true;