From 19f27d279d0b6305ae5e055a25121ff43cad4d74 Mon Sep 17 00:00:00 2001 From: Paul Kieckhefen Date: Wed, 30 May 2018 17:41:48 +0200 Subject: [PATCH] add postprocessing model for momentum-diameter correlation. --- .../forceModel/pdCorrelation/pdCorrelation.C | 334 ++++++++++++++++++ .../forceModel/pdCorrelation/pdCorrelation.H | 99 ++++++ 2 files changed, 433 insertions(+) create mode 100644 src/lagrangian/cfdemParticle/subModels/forceModel/pdCorrelation/pdCorrelation.C create mode 100644 src/lagrangian/cfdemParticle/subModels/forceModel/pdCorrelation/pdCorrelation.H diff --git a/src/lagrangian/cfdemParticle/subModels/forceModel/pdCorrelation/pdCorrelation.C b/src/lagrangian/cfdemParticle/subModels/forceModel/pdCorrelation/pdCorrelation.C new file mode 100644 index 00000000..44c428a3 --- /dev/null +++ b/src/lagrangian/cfdemParticle/subModels/forceModel/pdCorrelation/pdCorrelation.C @@ -0,0 +1,334 @@ +/*---------------------------------------------------------------------------*\ +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 . + +Description + calculates the correlation between particle momentum and diameter + +SourceFiles + pdCorrelation.C + +Contributing Author + 2018 Paul Kieckhefen, TUHH +\*---------------------------------------------------------------------------*/ + +#include "error.H" + +#include "pdCorrelation.H" +#include "addToRunTimeSelectionTable.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ + +// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * // + +defineTypeNameAndDebug(pdCorrelation, 0); + +addToRunTimeSelectionTable +( + forceModel, + pdCorrelation, + dictionary +); + +const dimensionSet dimMomentum(dimForce*dimTime); + +// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // + +// Construct from components +pdCorrelation::pdCorrelation +( + const dictionary& dict, + cfdemCloud& sm +) +: + forceModel(dict,sm), + propsDict_(dict.subDict(typeName + "Props")), + d_(nullptr), + p_(nullptr), + d2_(nullptr), + pd_(nullptr), + cg3_(nullptr), + dField_ + ( IOobject + ( + "dField", + sm.mesh().time().timeName(), + sm.mesh(), + IOobject::NO_READ, + IOobject::AUTO_WRITE + ), + sm.mesh(), + dimensionedScalar("zero", dimLength, 0), + "zeroGradient" + ), + pField_ + ( IOobject + ( + "pField", + sm.mesh().time().timeName(), + sm.mesh(), + IOobject::NO_READ, + IOobject::NO_WRITE + ), + sm.mesh(), + dimensionedVector("zero", dimMomentum, vector::zero) + ), + d2Field_ + ( IOobject + ( + "d2Field", + sm.mesh().time().timeName(), + sm.mesh(), + IOobject::NO_READ, + IOobject::NO_WRITE + ), + sm.mesh(), + dimensionedScalar("zero", dimLength*dimLength, 0) + ), + pdField_ + ( IOobject + ( + "pdCorrelation", + sm.mesh().time().timeName(), + sm.mesh(), + IOobject::NO_READ, + IOobject::AUTO_WRITE + ), + sm.mesh(), + dimensionedVector("zero", dimMomentum*dimLength, vector::zero), + "zeroGradient" + ), + cg3FieldPtr_(), + typeCG_ + ( + propsDict_.lookupOrDefault + ( + "coarseGrainingFactors", + scalarList(1, sm.cg()) + ) + ), + particleDensities_ + ( + propsDict_.lookupOrDefault + ( + "particleDensities", + scalarList(1, -1.) + ) + ), + constantCG_(typeCG_.size() > 0), + CG_(!(constantCG_ && abs(typeCG_[0] - 1.) < SMALL)) +{ + if ((particleDensities_[0] < 0) && !particleCloud_.getParticleDensities()) + { + FatalError<< "Please set the particle density either in LIGGGHTS" + << "or the pdCorrelationPropsDict." + << "In the first case, set getParticleDensities in the couplingProperties." + << abort(FatalError); + } + + allocateMyArrays(); + if (CG_) + { + cg3FieldPtr_.set + ( new volScalarField + ( IOobject + ( + "CG3Field", + sm.mesh().time().timeName(), + sm.mesh(), + IOobject::NO_READ, + IOobject::NO_WRITE + ), + sm.mesh(), + dimensionedScalar("zero", dimless, 0) + ) + ); + } + + dField_.write(); + pdField_.write(); + + + // init force sub model + setForceSubModels(propsDict_); +} + +// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * // + +pdCorrelation::~pdCorrelation() +{ + if (!constantCG_) + { + particleCloud_.dataExchangeM().destroy(cg3_, 1); + } + particleCloud_.dataExchangeM().destroy(d_, 1); + particleCloud_.dataExchangeM().destroy(p_, 3); + particleCloud_.dataExchangeM().destroy(d2_, 1); + particleCloud_.dataExchangeM().destroy(pd_, 3); +} + +// * * * * * * * * * * * * * * * private Member Functions * * * * * * * * * * * * * // +void pdCorrelation::allocateMyArrays() const +{ + // get memory for 2d arrays + double initVal = 0.0; + particleCloud_.dataExchangeM().allocateArray(d_, initVal, 1); + particleCloud_.dataExchangeM().allocateArray(p_, initVal, 3); + particleCloud_.dataExchangeM().allocateArray(d2_, initVal, 1); + particleCloud_.dataExchangeM().allocateArray(pd_, initVal, 3); + if (CG_) + { + particleCloud_.dataExchangeM().allocateArray(cg3_, initVal, 1); + } +} +// * * * * * * * * * * * * * * * public Member Functions * * * * * * * * * * * * * // + +void pdCorrelation::setForce() const +{ + const fvMesh& mesh = particleCloud_.mesh(); + + if (!mesh.write()) return; // skip if it's not write time + + allocateMyArrays(); + + const Switch densityFromList + ( + particleCloud_.getParticleTypes() + && particleCloud_.getParticleDensities() + ); + + const scalar piOverSix = constant::mathematical::pi / 6.; + + label celli(0); + scalar dp(0.); + scalar cg(typeCG_[0]); + scalar rhop(particleDensities_[0]); + label typei(0); + scalar particleV(0.); + for(int pi = 0; pi < particleCloud_.numberOfParticles(); ++pi) + { + celli = particleCloud_.cellIDs()[pi][0]; + if (celli < 0) continue; + + dp = particleCloud_.d(pi); + if (particleCloud_.getParticleTypes()) + { + typei = particleCloud_.particleTypes()[pi][0] - 1; + rhop = densityFromList ? particleDensities_[typei] + : particleCloud_.particleDensity(pi); + cg = constantCG_ ? typeCG_[typei] : typeCG_[0]; + } + + particleV = piOverSix * rhop * dp * dp * dp; + d_ [pi][0] = cg * cg * dp; + d2_[pi][0] = cg * dp * dp; + p_ [pi][0] = particleV * particleCloud_.velocities()[pi][0]; + p_ [pi][1] = particleV * particleCloud_.velocities()[pi][1]; + p_ [pi][2] = particleV * particleCloud_.velocities()[pi][2]; + pd_[pi][0] = p_[pi][0] * dp / cg; + pd_[pi][1] = p_[pi][1] * dp / cg; + pd_[pi][2] = p_[pi][2] * dp / cg; + + if (CG_) + { + cg3_[pi][0] = cg * cg * cg; + } + } + + dField_.primitiveFieldRef() = 0.0; + pField_.primitiveFieldRef() = vector::zero; + d2Field_.primitiveFieldRef() = 0.0; + pdField_.primitiveFieldRef() = vector::zero; + + particleCloud_.averagingM().setScalarSum + ( + dField_, + d_, + particleCloud_.particleWeights(), + nullptr + ); + particleCloud_.averagingM().setVectorSum + ( + pField_, + p_, + particleCloud_.particleWeights(), + nullptr + ); + particleCloud_.averagingM().setScalarSum + ( + d2Field_, + d2_, + particleCloud_.particleWeights(), + nullptr + ); + particleCloud_.averagingM().setVectorSum + ( + pdField_, + pd_, + particleCloud_.particleWeights(), + nullptr + ); + + if (CG_) + { + volScalarField& cg3Field_ = cg3FieldPtr_(); + particleCloud_.averagingM().setScalarSum + ( + cg3Field_, + cg3_, + particleCloud_.particleWeights(), + nullptr + ); + + scalar oneOverCG3(1.); + forAll(cg3Field_, celli) + { + if (cg3Field_[celli] < 1.) continue; + oneOverCG3 = 1./cg3Field_[celli]; + oneOverCG3 = oneOverCG3 * oneOverCG3 * oneOverCG3; + + dField_ [celli] *= oneOverCG3; + pField_ [celli] *= oneOverCG3; + d2Field_[celli] *= oneOverCG3; + pdField_[celli] *= oneOverCG3; + } + } + + scalar denom(0.); + forAll(dField_, celli) + { + if (dField_[celli] < SMALL) + { + pdField_[celli] = vector::zero; + continue; + } + denom = d2Field_[celli] - dField_[celli] * dField_[celli]; + if (denom < SMALL) + { + pdField_[celli] = vector::zero; + continue; + } + + pdField_[celli] -= dField_[celli] * pField_[celli]; + pdField_[celli] /= denom; + } +} + + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +} // End namespace Foam + +// ************************************************************************* // diff --git a/src/lagrangian/cfdemParticle/subModels/forceModel/pdCorrelation/pdCorrelation.H b/src/lagrangian/cfdemParticle/subModels/forceModel/pdCorrelation/pdCorrelation.H new file mode 100644 index 00000000..817e10af --- /dev/null +++ b/src/lagrangian/cfdemParticle/subModels/forceModel/pdCorrelation/pdCorrelation.H @@ -0,0 +1,99 @@ +/*---------------------------------------------------------------------------*\ +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 + calculates the granular kinetic energy + +SourceFiles + pdCorrelation.C +\*---------------------------------------------------------------------------*/ + +#ifndef pdCorrelation_H +#define pdCorrelation_H + +#include "forceModel.H" +#include "averagingModel.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ + +/*---------------------------------------------------------------------------*\ + Class pdCorrelation Declaration +\*---------------------------------------------------------------------------*/ + +class pdCorrelation +: + public forceModel +{ +private: + + dictionary propsDict_; + + mutable double **d_; + mutable double **p_; + mutable double **d2_; + mutable double **pd_; + mutable double **cg3_; + + mutable volScalarField dField_; + mutable volVectorField pField_; + mutable volScalarField d2Field_; + mutable volVectorField pdField_; + + mutable autoPtr cg3FieldPtr_; + + const scalarList typeCG_; + const scalarList particleDensities_; + const Switch constantCG_; + const Switch CG_; + + void allocateMyArrays() const; + +public: + + //- Runtime type information + TypeName("pdCorrelation"); + + // Constructors + + //- Construct from components + pdCorrelation + ( + const dictionary& dict, + cfdemCloud& sm + ); + + // Destructor + + ~pdCorrelation(); + + + // Member Functions + void setForce() const; + +}; + + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +} // End namespace Foam + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#endif + +// ************************************************************************* //