diff --git a/applications/solvers/cfdemSolverPisoScalar/cfdemSolverPisoScalar.C b/applications/solvers/cfdemSolverPisoScalar/cfdemSolverPisoScalar.C index 788c8f02..9fc59b59 100644 --- a/applications/solvers/cfdemSolverPisoScalar/cfdemSolverPisoScalar.C +++ b/applications/solvers/cfdemSolverPisoScalar/cfdemSolverPisoScalar.C @@ -25,7 +25,7 @@ License along with CFDEMcoupling. If not, see . Application - cfdemSolverPiso + cfdemSolverPisoScalar Description Transient solver for incompressible flow. diff --git a/src/lagrangian/cfdemParticle/Make/files b/src/lagrangian/cfdemParticle/Make/files index 6ddd9cc6..7cc1c491 100644 --- a/src/lagrangian/cfdemParticle/Make/files +++ b/src/lagrangian/cfdemParticle/Make/files @@ -111,6 +111,7 @@ $(voidFractionModels)/dividedVoidFractionMS/dividedVoidFractionMS.C $(voidFractionModels)/bigParticleVoidFraction/bigParticleVoidFraction.C $(voidFractionModels)/GaussVoidFraction/GaussVoidFraction.C $(voidFractionModels)/IBVoidFraction/IBVoidFraction.C +$(voidFractionModels)/trilinearVoidFraction/trilinearVoidFraction.C $(locateModels)/locateModel/locateModel.C $(locateModels)/locateModel/newLocateModel.C diff --git a/src/lagrangian/cfdemParticle/cfdTools/setupProbeModel.H b/src/lagrangian/cfdemParticle/cfdTools/setupProbeModel.H index d94133eb..cf15a61d 100644 --- a/src/lagrangian/cfdemParticle/cfdTools/setupProbeModel.H +++ b/src/lagrangian/cfdemParticle/cfdTools/setupProbeModel.H @@ -1,3 +1,2 @@ - //set probeModel parameters for this force model - particleCloud_.probeM().setOutputFile(); - particleCloud_.probeM().setCounter(); +//set probeModel parameters for this force model +if (probeIt_) { particleCloud_.probeM().setOutputFile(typeName+".logDat"); } diff --git a/src/lagrangian/cfdemParticle/cfdemCloud/cfdemCloud.C b/src/lagrangian/cfdemParticle/cfdemCloud/cfdemCloud.C index 359b51ac..79f2a03f 100644 --- a/src/lagrangian/cfdemParticle/cfdemCloud/cfdemCloud.C +++ b/src/lagrangian/cfdemParticle/cfdemCloud/cfdemCloud.C @@ -126,6 +126,7 @@ cfdemCloud::cfdemCloud mesh, dimensionedScalar("zero", dimensionSet(0,0,-1,0,0), 0) // 1/s ), + checkPeriodicCells_(false), turbulence_ ( mesh.lookupObject @@ -196,7 +197,7 @@ cfdemCloud::cfdemCloud clockModel::New ( couplingProperties_, - *this + mesh.time() ) ), smoothingModel_ @@ -225,6 +226,14 @@ cfdemCloud::cfdemCloud solveFlow_=Switch(couplingProperties_.lookup("solveFlow")); if (couplingProperties_.found("imExSplitFactor")) imExSplitFactor_ = readScalar(couplingProperties_.lookup("imExSplitFactor")); + + if(imExSplitFactor_ > 1.0) + FatalError << "You have set imExSplitFactor > 1 in your couplingProperties. Must be <= 1." + << abort(FatalError); + if(imExSplitFactor_ < 0.0) + FatalError << "You have set imExSplitFactor < 0 in your couplingProperties. Must be >= 0." + << abort(FatalError); + if (couplingProperties_.found("treatVoidCellsAsExplicitForce")) treatVoidCellsAsExplicitForce_ = readBool(couplingProperties_.lookup("treatVoidCellsAsExplicitForce")); if (couplingProperties_.found("verbose")) verbose_=true; @@ -289,7 +298,37 @@ cfdemCloud::cfdemCloud } dataExchangeM().setCG(); - if (!cgOK_ && cg_ > 1) FatalError<< "at least one of your models is not fit for cg !!!"<< abort(FatalError); + Switch cgWarnOnly_(couplingProperties_.lookupOrDefault("cgWarnOnly", true)); + if (!cgOK_ && cg_ > 1) + { + if (cgWarnOnly_) + Warning << "at least one of your models is not fit for cg !!!" << endl; + else + FatalError << "at least one of your models is not fit for cg !!!" << abort(FatalError); + } + + // check if simulation is a fully periodic box + const polyBoundaryMesh& patches = mesh_.boundaryMesh(); + int nPatchesCyclic(0); + int nPatchesNonCyclic(0); + forAll(patches, patchI) + { + const polyPatch& pp = patches[patchI]; + if (isA(pp) || isA(pp)) + ++nPatchesCyclic; + else if (!isA(pp)) + ++nPatchesNonCyclic; + } + + if (nPatchesNonCyclic == 0) + { + checkPeriodicCells_ = true; + } + else if (nPatchesCyclic > 0 && nPatchesNonCyclic > 0) + { + if (verbose_) Info << "nPatchesNonCyclic=" << nPatchesNonCyclic << ", nPatchesCyclic=" << nPatchesCyclic << endl; + Warning << "Periodic handing is disabled because the domain is not fully periodic!\n" << endl; + } } // * * * * * * * * * * * * * * * * Destructors * * * * * * * * * * * * * * // @@ -412,12 +451,15 @@ void cfdemCloud::setVectorAverages() ); if(verbose_) Info << "setVectorAverage done." << endl; } + // * * * * * * * * * * * * * * * public Member Functions * * * * * * * * * * * * * // + void cfdemCloud::checkCG(bool ok) { if(!cgOK_) return; if(!ok) cgOK_ = ok; } + void cfdemCloud::setPos(double**& pos) { for(int index = 0;index < numberOfParticles(); ++index) @@ -427,40 +469,32 @@ void cfdemCloud::setPos(double**& pos) } } } + // * * * * * * * * * * * * * * * ACCESS * * * * * * * * * * * * * // -label cfdemCloud::particleCell(int index) +label cfdemCloud::particleCell(int index) const { - label cellI = cellIDs()[index][0]; - return cellI; + return cellIDs()[index][0]; } -vector cfdemCloud::position(int index) +vector cfdemCloud::position(int index) const { - vector pos; - for(int i=0;i<3;i++) pos[i] = positions()[index][i]; - return pos; + return vector(positions()[index][0],positions()[index][1],positions()[index][2]); } -vector cfdemCloud::velocity(int index) +vector cfdemCloud::velocity(int index) const { - vector vel; - for(int i=0;i<3;i++) vel[i] = velocities()[index][i]; - return vel; + return vector(velocities()[index][0],velocities()[index][1],velocities()[index][2]); } -vector cfdemCloud::expForce(int index) +vector cfdemCloud::expForce(int index) const { - vector force; - for(int i=0;i<3;i++) force[i] = DEMForces()[index][i]; - return force; + return vector(DEMForces()[index][0],DEMForces()[index][1],DEMForces()[index][2]); } -vector cfdemCloud::fluidVel(int index) +vector cfdemCloud::fluidVel(int index) const { - vector vel; - for(int i=0;i<3;i++) vel[i] = fluidVels()[index][i]; - return vel; + return vector(fluidVels()[index][0],fluidVels()[index][1],fluidVels()[index][2]); } const forceModel& cfdemCloud::forceM(int i) @@ -468,43 +502,31 @@ const forceModel& cfdemCloud::forceM(int i) return forceModel_[i]; } -int cfdemCloud::nrForceModels() +label cfdemCloud::nrForceModels() const { return forceModels_.size(); } -int cfdemCloud::nrMomCoupleModels() +label cfdemCloud::nrMomCoupleModels() const { return momCoupleModels_.size(); } -scalar cfdemCloud::voidfraction(int index) +scalar cfdemCloud::voidfraction(int index) const { return voidfractions()[index][0]; } -label cfdemCloud::liggghtsCommandModelIndex(word name) +label cfdemCloud::liggghtsCommandModelIndex(word name) const { - int index=-1; forAll(liggghtsCommandModelList_,i) { if(liggghtsCommand()[i]().name() == name) { - index = i; - break; + return i; } } - return index; -} - -std::vector< std::vector >* cfdemCloud::getVprobe() -{ - return probeModel_->getVprobe(); -} - -std::vector< std::vector >* cfdemCloud::getSprobe() -{ - return probeModel_->getSprobe(); + return -1; } // * * * * * * * * * * * * * * * WRITE * * * * * * * * * * * * * // @@ -675,29 +697,6 @@ bool cfdemCloud::reAllocArrays() return false; } -bool cfdemCloud::reAllocArrays(int nP, bool forceRealloc) -{ - if( (numberOfParticlesChanged_ && !arraysReallocated_) || forceRealloc) - { - // get arrays of new length - dataExchangeM().allocateArray(positions_,0.,3,nP); - dataExchangeM().allocateArray(velocities_,0.,3,nP); - dataExchangeM().allocateArray(fluidVel_,0.,3,nP); - dataExchangeM().allocateArray(impForces_,0.,3,nP); - dataExchangeM().allocateArray(expForces_,0.,3,nP); - dataExchangeM().allocateArray(DEMForces_,0.,3,nP); - dataExchangeM().allocateArray(Cds_,0.,1,nP); - dataExchangeM().allocateArray(radii_,0.,1,nP); - dataExchangeM().allocateArray(voidfractions_,1.,voidFractionM().maxCellsPerParticle(),nP); - dataExchangeM().allocateArray(cellIDs_,-1,voidFractionM().maxCellsPerParticle(),nP); - dataExchangeM().allocateArray(particleWeights_,0.,voidFractionM().maxCellsPerParticle(),nP); - dataExchangeM().allocateArray(particleVolumes_,0.,voidFractionM().maxCellsPerParticle(),nP); - arraysReallocated_ = true; - return true; - } - return false; -} - tmp cfdemCloud::divVoidfractionTau(volVectorField& U,volScalarField& voidfraction) const { return diff --git a/src/lagrangian/cfdemParticle/cfdemCloud/cfdemCloud.H b/src/lagrangian/cfdemParticle/cfdemCloud/cfdemCloud.H index c199cd2b..bf014ddf 100644 --- a/src/lagrangian/cfdemParticle/cfdemCloud/cfdemCloud.H +++ b/src/lagrangian/cfdemParticle/cfdemCloud/cfdemCloud.H @@ -49,7 +49,7 @@ SourceFiles #include "fvCFD.H" #include "IFstream.H" -#include +#include "turbulenceModel.H" // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // @@ -160,6 +160,8 @@ protected: mutable volScalarField ddtVoidfraction_; + mutable Switch checkPeriodicCells_; + const turbulenceModel& turbulence_; autoPtr* forceModel_; @@ -228,41 +230,41 @@ public: void setPos(double **&); - word modelType(){ return modelType_; } + const word& modelType() const { return modelType_; } - label particleCell(int); + label particleCell(int) const; - vector position(int); + vector position(int) const; - vector velocity(int); + vector velocity(int) const; - vector expForce(int); + vector expForce(int) const; - vector fluidVel(int); + vector fluidVel(int) const; virtual const forceModel& forceM(int); - virtual int nrForceModels(); + virtual label nrForceModels() const; - virtual int nrMomCoupleModels(); + virtual label nrMomCoupleModels() const; - scalar voidfraction(int); + scalar voidfraction(int) const; - label liggghtsCommandModelIndex(word); + label liggghtsCommandModelIndex(word) const; - inline void setCG(double) const; + inline void setCG(double); - inline const scalar& cg() const; + inline scalar cg() const; - inline const bool& impDEMdrag() const; + inline bool impDEMdrag() const; - inline const bool& impDEMdragAcc() const; + inline bool impDEMdragAcc() const; - inline const scalar& imExSplitFactor() const; + inline scalar imExSplitFactor() const; - inline const bool& treatVoidCellsAsExplicitForce() const; + inline bool treatVoidCellsAsExplicitForce() const; - inline const bool& ignore() const; + inline bool ignore() const; inline const fvMesh& mesh() const; @@ -300,13 +302,13 @@ public: inline double ** particleWeights() const; - virtual inline label body(int); + virtual inline label body(int) const; - virtual inline double particleVolume(int); + virtual inline double particleVolume(int) const; - inline scalar radius(int); + inline scalar radius(int) const; - virtual inline double d(int); + virtual inline double d(int) const; inline scalar d32(bool recalc=true); virtual inline double dMin() {return -1;} @@ -322,11 +324,11 @@ public: //access to the particle's rotation and torque data virtual inline double ** DEMTorques() const {return NULL;} virtual inline double ** omegaArray() const {return NULL;} - virtual vector omega(int) const {return Foam::vector(0,0,0);} + virtual vector omega(int) const {return vector(0,0,0);} //access to the particles' orientation information virtual inline double ** exArray() const {return NULL;} - virtual vector ex(int) const {return Foam::vector(0,0,0);} + virtual vector ex(int) const {return vector(0,0,0);} //Detector if SRF module is enable or not virtual inline bool SRFOn(){return false;} @@ -339,7 +341,7 @@ public: inline bool arraysReallocated() const; - inline const wordList& forceModels(); + inline const wordList& forceModels() const; inline const voidFractionModel& voidFractionM() const; @@ -347,11 +349,11 @@ public: inline const momCoupleModel& momCoupleM(int) const; - inline const dataExchangeModel& dataExchangeM() const; + inline dataExchangeModel& dataExchangeM(); inline const IOModel& IOM() const; - inline const probeModel& probeM() const; + inline probeModel& probeM(); inline const averagingModel& averagingM() const; @@ -376,13 +378,10 @@ public: virtual bool reAllocArrays(); - virtual bool reAllocArrays(int nP, bool forceRealloc); //force number of particles during reallocation - - // IO - void writeScalarFieldToTerminal(double**&); + void writeScalarFieldToTerminal(double**&) const; - void writeVectorFieldToTerminal(double**&); + void writeVectorFieldToTerminal(double**&) const; // functions tmp divVoidfractionTau(volVectorField& ,volScalarField&) const; @@ -397,11 +396,9 @@ public: void resetArray(double**&,int,int,double resetVal=0.); - std::vector< std::vector >* getVprobe(); - - std::vector< std::vector >* getSprobe(); - void otherForces(volVectorField&); + + bool checkPeriodicCells() { return checkPeriodicCells_; } }; diff --git a/src/lagrangian/cfdemParticle/cfdemCloud/cfdemCloudI.H b/src/lagrangian/cfdemParticle/cfdemCloud/cfdemCloudI.H index 84cd29c5..2fbe1744 100644 --- a/src/lagrangian/cfdemParticle/cfdemCloud/cfdemCloudI.H +++ b/src/lagrangian/cfdemParticle/cfdemCloud/cfdemCloudI.H @@ -44,38 +44,38 @@ namespace Foam { // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // -inline void cfdemCloud::setCG(double cg) const +inline void cfdemCloud::setCG(double cg) { cg_ = cg; Info << "cg is set to: " << cg_ << endl; } -inline const bool& cfdemCloud::impDEMdrag() const +inline bool cfdemCloud::impDEMdrag() const { return impDEMdrag_; } -inline const bool& cfdemCloud::impDEMdragAcc() const +inline bool cfdemCloud::impDEMdragAcc() const { return impDEMdragAcc_; } -inline const scalar& cfdemCloud::imExSplitFactor() const +inline scalar cfdemCloud::imExSplitFactor() const { return imExSplitFactor_; } -inline const bool& cfdemCloud::treatVoidCellsAsExplicitForce() const +inline bool cfdemCloud::treatVoidCellsAsExplicitForce() const { return treatVoidCellsAsExplicitForce_; } -inline const scalar& cfdemCloud::cg() const +inline scalar cfdemCloud::cg() const { return cg_; } -inline const bool& cfdemCloud::ignore() const +inline bool cfdemCloud::ignore() const { return ignore_; } @@ -177,24 +177,24 @@ inline double ** cfdemCloud::particleWeights() const return particleWeights_; } -inline label cfdemCloud::body(int index) +inline label cfdemCloud::body(int index) const { return index; } -inline double cfdemCloud::particleVolume(int index) +inline double cfdemCloud::particleVolume(int index) const { return particleV_[index][0]; } -inline scalar cfdemCloud::radius(int index) +inline scalar cfdemCloud::radius(int index) const { return radii_[index][0]; } -inline double cfdemCloud::d(int index) +inline double cfdemCloud::d(int index) const { - return 2*radii_[index][0]; + return 2.*radii_[index][0]; } inline double cfdemCloud::d32(bool recalc) @@ -237,7 +237,7 @@ inline bool cfdemCloud::arraysReallocated() const return arraysReallocated_; } -inline const wordList& cfdemCloud::forceModels() +inline const wordList& cfdemCloud::forceModels() const { return forceModels_; } @@ -252,9 +252,9 @@ inline const momCoupleModel& cfdemCloud::momCoupleM(int i) const return momCoupleModel_[i]; } -inline const dataExchangeModel& cfdemCloud::dataExchangeM() const +inline dataExchangeModel& cfdemCloud::dataExchangeM() { - return dataExchangeModel_; + return dataExchangeModel_(); } inline const IOModel& cfdemCloud::IOM() const @@ -262,9 +262,9 @@ inline const IOModel& cfdemCloud::IOM() const return IOModel_; } -inline const probeModel& cfdemCloud::probeM() const +inline probeModel& cfdemCloud::probeM() { - return probeModel_; + return probeModel_(); } inline const voidFractionModel& cfdemCloud::voidFractionM() const diff --git a/src/lagrangian/cfdemParticle/cfdemCloud/cfdemCloudIO.C b/src/lagrangian/cfdemParticle/cfdemCloud/cfdemCloudIO.C index 51744a79..716c0eab 100644 --- a/src/lagrangian/cfdemParticle/cfdemCloud/cfdemCloudIO.C +++ b/src/lagrangian/cfdemParticle/cfdemCloud/cfdemCloudIO.C @@ -48,7 +48,7 @@ namespace Foam // * * * * * * * * * * * * * * * IO * * * * * * * * * * * * * // -void cfdemCloud::writeScalarFieldToTerminal(double**& array) +void cfdemCloud::writeScalarFieldToTerminal(double**& array) const { // init double array for (int i=0; i* energyModel_; - + autoPtr thermCondModel_; - + autoPtr chemistryModel_; - + void calcEnergyContributions(); void speciesExecute(); - + public: friend class energyModel; @@ -89,23 +89,23 @@ public: const energyModel& energyM(int); - + const thermCondModel& thermCondM(); - + const chemistryModel& chemistryM(); - int nrEnergyModels(); + label nrEnergyModels() const; + + inline const wordList& energyModels() const; + + bool implicitEnergyModel() const; - inline const wordList& energyModels(); - - bool& implicitEnergyModel(); - void energyContributions(volScalarField&); - + void energyCoefficients(volScalarField&); - + bool evolve(volScalarField&,volVectorField&,volVectorField&); - + void postFlow(); }; diff --git a/src/lagrangian/cfdemParticle/derived/cfdemCloudEnergy/cfdemCloudEnergyI.H b/src/lagrangian/cfdemParticle/derived/cfdemCloudEnergy/cfdemCloudEnergyI.H index 14231a4b..12983d18 100644 --- a/src/lagrangian/cfdemParticle/derived/cfdemCloudEnergy/cfdemCloudEnergyI.H +++ b/src/lagrangian/cfdemParticle/derived/cfdemCloudEnergy/cfdemCloudEnergyI.H @@ -25,7 +25,7 @@ namespace Foam { // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // -inline const wordList& cfdemCloudEnergy::energyModels() +inline const wordList& cfdemCloudEnergy::energyModels() const { return energyModels_; } diff --git a/src/lagrangian/cfdemParticle/derived/cfdemCloudIB/cfdemCloudIB.C b/src/lagrangian/cfdemParticle/derived/cfdemCloudIB/cfdemCloudIB.C index 3f578ef5..8fec93d2 100644 --- a/src/lagrangian/cfdemParticle/derived/cfdemCloudIB/cfdemCloudIB.C +++ b/src/lagrangian/cfdemParticle/derived/cfdemCloudIB/cfdemCloudIB.C @@ -91,7 +91,7 @@ bool cfdemCloudIB::reAllocArrays() if(cfdemCloud::reAllocArrays()) { // get arrays of new length - dataExchangeM().allocateArray(angularVelocities_,0,3); + dataExchangeM().allocateArray(angularVelocities_,0.,3); return true; } return false; @@ -224,11 +224,9 @@ void cfdemCloudIB::calcVelocityCorrection } -vector cfdemCloudIB::angularVelocity(int index) +vector cfdemCloudIB::angularVelocity(int index) const { - vector vel; - for(int i=0;i<3;i++) vel[i] = angularVelocities_[index][i]; - return vel; + return vector(angularVelocities_[index][0],angularVelocities_[index][1],angularVelocities_[index][2]); } // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // diff --git a/src/lagrangian/cfdemParticle/derived/cfdemCloudIB/cfdemCloudIB.H b/src/lagrangian/cfdemParticle/derived/cfdemCloudIB/cfdemCloudIB.H index 1aaffa5a..d1db244d 100644 --- a/src/lagrangian/cfdemParticle/derived/cfdemCloudIB/cfdemCloudIB.H +++ b/src/lagrangian/cfdemParticle/derived/cfdemCloudIB/cfdemCloudIB.H @@ -90,7 +90,7 @@ public: void calcVelocityCorrection(volScalarField&,volVectorField&,volScalarField&,volScalarField&); // this could be moved to an IB mom couple model // Access - vector angularVelocity(int); + vector angularVelocity(int) const; inline double ** angularVelocities() const { diff --git a/src/lagrangian/cfdemParticle/derived/cfdemCloudMS/cfdemCloudMS.C b/src/lagrangian/cfdemParticle/derived/cfdemCloudMS/cfdemCloudMS.C index 1b34cad0..40e5ddde 100644 --- a/src/lagrangian/cfdemParticle/derived/cfdemCloudMS/cfdemCloudMS.C +++ b/src/lagrangian/cfdemParticle/derived/cfdemCloudMS/cfdemCloudMS.C @@ -128,14 +128,14 @@ void cfdemCloudMS::getDEMdata() //- save clump volume and mass double **typeDH(NULL); - dataExchangeM().allocateArray(typeDH,-1,1,nClumpTypes()+1); + dataExchangeM().allocateArray(typeDH,-1.,1,nClumpTypes()+1); if(manDHdev_) // use manually defined dH { for(int k = 1;k <= nClumpTypes(); k++) typeDH[k][0]=dHbyV_[k-1]*typeVol_[k]; } else // calc dH from volAeqivalent shpere - { + { for(int k = 1;k <= nClumpTypes(); k++) typeDH[k][0]=pow(typeVol_[k]*1.9099,1./3.); // 6/pi=1.9099 // calc a hydraulic diameter as d of vol equal sphere } @@ -144,7 +144,7 @@ void cfdemCloudMS::getDEMdata() for(int ind = 0;ind < numberOfClumps(); ind++) { ct=clumpType()[0][ind]; - clumpVol_[ind][0] = typeVol_[ct]; + clumpVol_[ind][0] = typeVol_[ct]; clumpDH_[ind][0]=typeDH[ct][0]; //Info << "ct=" << ct << endl; //Info << "clumpVol()[ind][0]=" << clumpVol()[ind][0] << endl; @@ -192,18 +192,18 @@ bool cfdemCloudMS::reAllocArrays() if(cfdemCloud::reAllocArrays()) { // get arrays of new length - dataExchangeM().allocateArray(positionsCM_,0,3,"nbodies"); - dataExchangeM().allocateArray(velocitiesCM_,0,3,"nbodies"); + dataExchangeM().allocateArray(positionsCM_,0.,3,"nbodies"); + dataExchangeM().allocateArray(velocitiesCM_,0.,3,"nbodies"); dataExchangeM().allocateArray(cellIDsCM_,-1,1,"nbodies"); dataExchangeM().allocateArray(bodies_,0,1); dataExchangeM().allocateArray(nrigids_,0,1,"nbodies"); dataExchangeM().allocateArray(clumpType_,0,1,"nbodies"); - dataExchangeM().allocateArray(clumpVol_,0,1,"nbodies"); + dataExchangeM().allocateArray(clumpVol_,0.,1,"nbodies"); dataExchangeM().allocateArray(clumpDH_,1.,1,"nbodies"); - dataExchangeM().allocateArray(clumpWeights_,1,1,"nbodies"); - dataExchangeM().allocateArray(impForcesCM_,0,3,"nbodies"); - dataExchangeM().allocateArray(expForcesCM_,0,3,"nbodies"); - dataExchangeM().allocateArray(DEMForcesCM_,0,3,"nbodies"); + dataExchangeM().allocateArray(clumpWeights_,1.,1,"nbodies"); + dataExchangeM().allocateArray(impForcesCM_,0.,3,"nbodies"); + dataExchangeM().allocateArray(expForcesCM_,0.,3,"nbodies"); + dataExchangeM().allocateArray(DEMForcesCM_,0.,3,"nbodies"); return true; } return false; diff --git a/src/lagrangian/cfdemParticle/derived/cfdemCloudMS/cfdemCloudMS.H b/src/lagrangian/cfdemParticle/derived/cfdemCloudMS/cfdemCloudMS.H index 9888149b..df3f8352 100644 --- a/src/lagrangian/cfdemParticle/derived/cfdemCloudMS/cfdemCloudMS.H +++ b/src/lagrangian/cfdemParticle/derived/cfdemCloudMS/cfdemCloudMS.H @@ -117,9 +117,9 @@ public: // Member Functions // Access - inline label body(int); + inline label body(int) const; - inline double particleVolume(int); + inline double particleVolume(int) const; inline vector positionCM(int); @@ -127,9 +127,9 @@ public: inline label cellIDCM(int); - inline label nrigid(int); + inline label nrigid(int) const; - inline int nrForceModels(); + inline label nrForceModels() const; inline double **& positionsCM() const; diff --git a/src/lagrangian/cfdemParticle/derived/cfdemCloudMS/cfdemCloudMSI.H b/src/lagrangian/cfdemParticle/derived/cfdemCloudMS/cfdemCloudMSI.H index ecb76590..10abefb8 100644 --- a/src/lagrangian/cfdemParticle/derived/cfdemCloudMS/cfdemCloudMSI.H +++ b/src/lagrangian/cfdemParticle/derived/cfdemCloudMS/cfdemCloudMSI.H @@ -36,12 +36,12 @@ namespace Foam // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // -inline label cfdemCloudMS::body(int index) +inline label cfdemCloudMS::body(int index) const { return bodies_[0][index]-1; } -inline double cfdemCloudMS::particleVolume(int index) +inline double cfdemCloudMS::particleVolume(int index) const { int ind = body(index); // particle to clump ID @@ -50,7 +50,7 @@ inline double cfdemCloudMS::particleVolume(int index) Vp=clumpVol_[ind][0]; else Warning << "ind=" << ind << endl; - + int nR(nrigid(ind)); if(nR>0) Vp/=nR; return Vp; @@ -75,12 +75,12 @@ inline label cfdemCloudMS::cellIDCM(int index) return cellIDsCM_[index][0]; } -inline label cfdemCloudMS::nrigid(int index) +inline label cfdemCloudMS::nrigid(int index) const { return nrigids_[0][index]; } -inline int cfdemCloudMS::nrForceModels() +inline label cfdemCloudMS::nrForceModels() const { return forceModels_.size(); } diff --git a/src/lagrangian/cfdemParticle/subModels/averagingModel/averagingModel/averagingModel.C b/src/lagrangian/cfdemParticle/subModels/averagingModel/averagingModel/averagingModel.C index b0d8b4bc..1f159643 100644 --- a/src/lagrangian/cfdemParticle/subModels/averagingModel/averagingModel/averagingModel.C +++ b/src/lagrangian/cfdemParticle/subModels/averagingModel/averagingModel/averagingModel.C @@ -262,7 +262,7 @@ void averagingModel::setDSauter for(int index=0; index< particleCloud_.numberOfParticles(); index++) { if(myParticleType!=0) //in case a particle type is specified, only consider particles of the right type - if(myParticleType != particleCloud_.particleType(index)) continue; + if(myParticleType != particleCloud_.particleType(index)) continue; radius = particleCloud_.radii()[index][0] / scale_; //the primary particle diameter radiusPow2 = radius*radius; @@ -313,7 +313,7 @@ void averagingModel::resetWeightFields() const } -void Foam::averagingModel::undoWeightFields(double**const& mask) const +void averagingModel::undoWeightFields(double**const& mask) const { for(int index=0; index< particleCloud_.numberOfParticles(); index++) { @@ -326,41 +326,22 @@ void Foam::averagingModel::undoWeightFields(double**const& mask) const } } -tmp Foam::averagingModel::UsInterp() const +tmp averagingModel::UsInterp() const { - tmp tsource - ( - new volVectorField - ( - IOobject - ( - "Us_averagingModel", - particleCloud_.mesh().time().timeName(), - particleCloud_.mesh(), - IOobject::NO_READ, - IOobject::NO_WRITE - ), - particleCloud_.mesh(), - dimensionedVector - ( - "zero", - dimensionSet(0, 1, -1, 0, 0), - vector::zero - ) - ) - ); - if (particleCloud_.dataExchangeM().couplingStep() > 1) { - tsource.ref() = (1 - particleCloud_.dataExchangeM().timeStepFraction()) * UsPrev_ - + particleCloud_.dataExchangeM().timeStepFraction() * UsNext_; + return tmp + ( + new volVectorField("Us_averagingModel", (1. - particleCloud_.dataExchangeM().timeStepFraction()) * UsPrev_ + particleCloud_.dataExchangeM().timeStepFraction() * UsNext_) + ); } else { - tsource.ref() = UsNext_; + return tmp + ( + new volVectorField("Us_averagingModel", UsNext_) + ); } - - return tsource; } // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // diff --git a/src/lagrangian/cfdemParticle/subModels/chemistryModel/species/species.C b/src/lagrangian/cfdemParticle/subModels/chemistryModel/species/species.C index fc27c319..e0aaa6f5 100644 --- a/src/lagrangian/cfdemParticle/subModels/chemistryModel/species/species.C +++ b/src/lagrangian/cfdemParticle/subModels/chemistryModel/species/species.C @@ -147,11 +147,11 @@ species::species species::~species() { - delete partTemp_; - delete partRho_; + particleCloud_.dataExchangeM().destroy(partTemp_,1); + particleCloud_.dataExchangeM().destroy(partRho_,1); - for (int i=0; i=0) + if (cellI >= 0) { if(interpolation_) { @@ -213,12 +213,12 @@ void species::execute() partRho_[index][0]=rhofluid; for (int i=0; i=0 && index < 2) + if(particleCloud_.verbose() && index >= 0 && index < 2) { /*for(int i =0; i #include "clockModel.H" #include - +#include // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // @@ -48,19 +48,18 @@ defineRunTimeSelectionTable(clockModel, dictionary); // * * * * * * * * * * * * * public Member Functions * * * * * * * * * * * * // -void Foam::clockModel::start(int pos) const +void clockModel::start(int pos) const { start(pos,""); - return; } -void Foam::clockModel::start(int pos, const std::string& ident) const +void clockModel::start(int pos, const std::string& ident) const { - if (particleCloud_.mesh().time().value() > startTime_) + if (time_.value() > startTime_) { if (pos >= n_) // alternatively one fixed size? { - n_ = 2*n_; + n_ += 32; deltaT_.resize(n_,0); identifier_.resize(n_,""); nOfRuns_.resize(n_,0); @@ -73,14 +72,13 @@ void Foam::clockModel::start(int pos, const std::string& ident) const parent_[pos]=curParent_; curParent_ = pos; nOfRuns_[pos] += 1; - deltaT_[pos]-=std::clock(); + deltaT_[pos] -= std::clock(); } - return; } -void Foam::clockModel::stop() const +void clockModel::stop() const { - if (particleCloud_.mesh().time().value() > startTime_) + if (time_.value() > startTime_) { deltaT_[curParent_] += std::clock(); curLev_ -= 1; @@ -93,17 +91,16 @@ void Foam::clockModel::stop() const curParent_ = -1; } } - return; } -void Foam::clockModel::stop(const std::string& ident) const +void clockModel::stop(const std::string& ident) const { - if (particleCloud_.mesh().time().value() > startTime_) + if (time_.value() > startTime_) { deltaT_[curParent_] += std::clock(); if (curParent_ > 0 && identifier_[curParent_].compare(ident) != 0) { - Pout<<"Warning: stop identifier did not equal start identifier! "<= 0) @@ -115,10 +112,9 @@ void Foam::clockModel::stop(const std::string& ident) const curParent_ = -1; } } - return; } -std::string Foam::clockModel::eval() const +std::string clockModel::eval() const { std::ostringstream strs("Measurements in CPU-seconds:\n"); strs << "Name\tdeltaT\tnOfRuns\tlevel\tparentNr\tparentName\n"; @@ -149,7 +145,7 @@ std::string Foam::clockModel::eval() const return strs.str(); } -void Foam::clockModel::evalFile() const +void clockModel::evalFile() const { std::ofstream outFile; std::string fileName(path_/"timeEval.txt"); @@ -159,10 +155,10 @@ void Foam::clockModel::evalFile() const outFile.close(); } -void Foam::clockModel::evalPar() const +void clockModel::evalPar() const { int myrank, numprocs; - MPI_Comm_rank(MPI_COMM_WORLD,&myrank); + MPI_Comm_rank(MPI_COMM_WORLD, &myrank); MPI_Comm_size(MPI_COMM_WORLD, &numprocs); std::ofstream outFile; @@ -227,29 +223,27 @@ void Foam::clockModel::evalPar() const outFile << strs.str(); outFile.close(); } - - return; } -void Foam::clockModel::initElems() +void clockModel::initElems() { //init elems - for (int i = 0;i < n_; i++) + for (int i = 0; i < n_; ++i) { deltaT_[i] = 0; - identifier_[i] = ""; + identifier_[i].clear(); nOfRuns_[i] = 0; level_[i] = -1; parent_[i] = -2; } } -std::vector Foam::clockModel::calcShift() const +std::vector clockModel::calcShift() const { std::vector shifts(n_, 0); - for (int i=1; i Foam::clockModel::calcShift() const return shifts; } -void Foam::clockModel::normHist() const +void clockModel::normHist() const { int myrank, numprocs; MPI_Comm_rank(MPI_COMM_WORLD,&myrank); @@ -300,28 +294,28 @@ void Foam::clockModel::normHist() const Info << "===========================" << endl; getRAMUsage(); - return; } -void Foam::clockModel::plotHist(double buffIn,const std::string& identifier,int numprocs,int myrank) const +void clockModel::plotHist(double buffIn,const std::string& identifier,int numprocs,int myrank) const { double* globalTime_all = NULL; + if (myrank == 0) globalTime_all = new double[numprocs]; MPI_Gather(&buffIn, 1, MPI_DOUBLE, globalTime_all, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD); if (myrank == 0) for (int j=0; j(SwapMem)/1024.0; //kB -> MB - double RssMB = static_cast(RssMem)/1024.0; + double SwapMB = SwapMem/1024.0; //kB -> MB + double RssMB = RssMem/1024.0; inFile.close(); // set up communication between Procs and plot Stuff Info << " RAM USAGE HISTOGRAM in MB" << endl; plotHist(RssMB,"RSS memory used",numprocs,myrank); - if (SwapMem > 0) + + if (SwapMem > 0) { plotHist(SwapMB,"WARNING: Swap",numprocs,myrank); } @@ -398,26 +392,23 @@ void Foam::clockModel::getRAMUsage() const //Pout << "SWAP Memory used: " << SwapMem <<"MB\n"; //Pout << "Rss Memory used: " << RssMem <<"MB\n"; - - return; } // * * * * * * * * * * * * * private Member Functions * * * * * * * * * * * * // // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // // Construct from components -Foam::clockModel::clockModel +clockModel::clockModel ( const dictionary& dict, - cfdemCloud& sm + const Time& time ) : dict_(dict), - particleCloud_(sm), + time_(time), path_("clockData"), - startTime_(sm.mesh().time().startTime().value()+sm.mesh().time().deltaT().value()+SMALL), // delay start of measurement by deltaT - //startTime_(0), //no delay - n_(30), + startTime_(time.startTime().value()+time.deltaT().value()+SMALL), // delay start of measurement by deltaT + n_(32), deltaT_(n_), identifier_(n_), nOfRuns_(n_), @@ -432,7 +423,7 @@ Foam::clockModel::clockModel // * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * // -Foam::clockModel::~clockModel() +clockModel::~clockModel() {} diff --git a/src/lagrangian/cfdemParticle/subModels/clockModel/clockModel/clockModel.H b/src/lagrangian/cfdemParticle/subModels/clockModel/clockModel/clockModel.H index 942eae1b..ac9888c4 100644 --- a/src/lagrangian/cfdemParticle/subModels/clockModel/clockModel/clockModel.H +++ b/src/lagrangian/cfdemParticle/subModels/clockModel/clockModel/clockModel.H @@ -42,10 +42,8 @@ SourceFiles #define START(x) start(__COUNTER__,x) #include "fvCFD.H" -#include "cfdemCloud.H" -#include "dataExchangeModel.H" - #include + // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // namespace Foam @@ -61,7 +59,7 @@ protected: // Protected data const dictionary& dict_; - cfdemCloud& particleCloud_; + const Time& time_; fileName path_; scalar startTime_; @@ -88,9 +86,9 @@ public: dictionary, ( const dictionary& dict, - cfdemCloud& sm + const Time& time ), - (dict,sm) + (dict,time) ); @@ -100,7 +98,7 @@ public: clockModel ( const dictionary& dict, - cfdemCloud& sm + const Time& time ); @@ -114,7 +112,7 @@ public: static autoPtr New ( const dictionary& dict, - cfdemCloud& sm + const Time& time ); @@ -129,9 +127,9 @@ public: virtual void evalPar() const; void initElems(); std::vector calcShift() const; //detects empty indices in vector, when times are evaluated - void Hist() const; //calc Histogram - virtual void normHist() const; //calc normalized Histogram - void plotHist(double,const std::string&,int,int) const; //plot histogramm to terminal + void Hist() const; //calc Histogram + virtual void normHist() const; //calc normalized Histogram + void plotHist(double,const std::string&,int,int) const; //plot histogramm to terminal void getRAMUsage() const; }; diff --git a/src/lagrangian/cfdemParticle/subModels/clockModel/clockModel/newClockModel.C b/src/lagrangian/cfdemParticle/subModels/clockModel/clockModel/newClockModel.C index f789525a..c3816d23 100644 --- a/src/lagrangian/cfdemParticle/subModels/clockModel/clockModel/newClockModel.C +++ b/src/lagrangian/cfdemParticle/subModels/clockModel/clockModel/newClockModel.C @@ -44,7 +44,7 @@ namespace Foam autoPtr clockModel::New ( const dictionary& dict, - cfdemCloud& sm + const Time& time ) { word clockModelType @@ -52,7 +52,7 @@ autoPtr clockModel::New dict.lookup("clockModel") ); - Info<< "Selecting clockModel " + Info << "Selecting clockModel " << clockModelType << endl; @@ -73,7 +73,7 @@ autoPtr clockModel::New << abort(FatalError); } - return autoPtr(cstrIter()(dict,sm)); + return autoPtr(cstrIter()(dict,time)); } diff --git a/src/lagrangian/cfdemParticle/subModels/clockModel/noClock/noClock.C b/src/lagrangian/cfdemParticle/subModels/clockModel/noClock/noClock.C index 550724fb..b224b02c 100644 --- a/src/lagrangian/cfdemParticle/subModels/clockModel/noClock/noClock.C +++ b/src/lagrangian/cfdemParticle/subModels/clockModel/noClock/noClock.C @@ -56,10 +56,10 @@ addToRunTimeSelectionTable noClock::noClock ( const dictionary& dict, - cfdemCloud& sm + const Time& time ) : - clockModel(dict,sm) + clockModel(dict,time) { initElems(); } diff --git a/src/lagrangian/cfdemParticle/subModels/clockModel/noClock/noClock.H b/src/lagrangian/cfdemParticle/subModels/clockModel/noClock/noClock.H index 7af54865..11797501 100644 --- a/src/lagrangian/cfdemParticle/subModels/clockModel/noClock/noClock.H +++ b/src/lagrangian/cfdemParticle/subModels/clockModel/noClock/noClock.H @@ -47,7 +47,7 @@ namespace Foam { /*---------------------------------------------------------------------------*\ - Class noDrag Declaration + Class noClock Declaration \*---------------------------------------------------------------------------*/ class noClock @@ -67,7 +67,7 @@ public: noClock ( const dictionary& dict, - cfdemCloud& sm + const Time& time ); // Destructor diff --git a/src/lagrangian/cfdemParticle/subModels/clockModel/standardClock/standardClock.C b/src/lagrangian/cfdemParticle/subModels/clockModel/standardClock/standardClock.C index bda6e6ce..25c11b14 100644 --- a/src/lagrangian/cfdemParticle/subModels/clockModel/standardClock/standardClock.C +++ b/src/lagrangian/cfdemParticle/subModels/clockModel/standardClock/standardClock.C @@ -30,7 +30,6 @@ Description \*---------------------------------------------------------------------------*/ #include "error.H" -#include "IOModel.H" #include "standardClock.H" #include "addToRunTimeSelectionTable.H" @@ -57,12 +56,14 @@ addToRunTimeSelectionTable standardClock::standardClock ( const dictionary& dict, - cfdemCloud& sm + const Time& time ) : - clockModel(dict,sm) + clockModel(dict,time) { - path_=particleCloud_.IOM().createTimeDir(path_); + path_ = path_/time_.timeName(); + mkDir(path_,0777); + initElems(); } diff --git a/src/lagrangian/cfdemParticle/subModels/clockModel/standardClock/standardClock.H b/src/lagrangian/cfdemParticle/subModels/clockModel/standardClock/standardClock.H index 3b62490c..b4bf2a35 100644 --- a/src/lagrangian/cfdemParticle/subModels/clockModel/standardClock/standardClock.H +++ b/src/lagrangian/cfdemParticle/subModels/clockModel/standardClock/standardClock.H @@ -47,7 +47,7 @@ namespace Foam { /*---------------------------------------------------------------------------*\ - Class noDrag Declaration + Class standardClock Declaration \*---------------------------------------------------------------------------*/ class standardClock @@ -67,7 +67,7 @@ public: standardClock ( const dictionary& dict, - cfdemCloud& sm + const Time& time ); // Destructor diff --git a/src/lagrangian/cfdemParticle/subModels/dataExchangeModel/dataExchangeModel/dataExchangeModel.C b/src/lagrangian/cfdemParticle/subModels/dataExchangeModel/dataExchangeModel/dataExchangeModel.C index dfe718c1..df0ed4c4 100755 --- a/src/lagrangian/cfdemParticle/subModels/dataExchangeModel/dataExchangeModel/dataExchangeModel.C +++ b/src/lagrangian/cfdemParticle/subModels/dataExchangeModel/dataExchangeModel/dataExchangeModel.C @@ -45,7 +45,7 @@ defineRunTimeSelectionTable(dataExchangeModel, dictionary); // * * * * * * * * * * * * * * protected Member Functions * * * * * * * * * * * * * // -void Foam::dataExchangeModel::setNumberOfParticles(int numberOfParticles) const +void dataExchangeModel::setNumberOfParticles(int numberOfParticles) const { particleCloud_.setNumberOfParticles(numberOfParticles); } @@ -55,7 +55,7 @@ void Foam::dataExchangeModel::setNumberOfParticles(int numberOfParticles) const //==== // double ** -void Foam::dataExchangeModel::allocateArray +void dataExchangeModel::allocateArray ( double**& array, double initVal, @@ -77,7 +77,7 @@ void Foam::dataExchangeModel::allocateArray } } -void Foam::dataExchangeModel::allocateArray +void dataExchangeModel::allocateArray ( double**& array, double initVal, @@ -92,7 +92,7 @@ void Foam::dataExchangeModel::allocateArray allocateArray(array,initVal,width,len); } -void Foam::dataExchangeModel::destroy(double** array,int /*len*/) const +void dataExchangeModel::destroy(double** array,int /*len*/) const { if (array == NULL) return; @@ -103,7 +103,7 @@ void Foam::dataExchangeModel::destroy(double** array,int /*len*/) const //==== // int ** -void Foam::dataExchangeModel::allocateArray +void dataExchangeModel::allocateArray ( int**& array, int initVal, @@ -125,7 +125,7 @@ void Foam::dataExchangeModel::allocateArray } } -void Foam::dataExchangeModel::allocateArray +void dataExchangeModel::allocateArray ( int**& array, int initVal, @@ -140,7 +140,7 @@ void Foam::dataExchangeModel::allocateArray allocateArray(array,initVal,width,len); } -void Foam::dataExchangeModel::destroy(int** array,int /*len*/) const +void dataExchangeModel::destroy(int** array,int /*len*/) const { if (array == NULL) return; @@ -152,7 +152,7 @@ void Foam::dataExchangeModel::destroy(int** array,int /*len*/) const //==== // int * -void Foam::dataExchangeModel::allocateArray +void dataExchangeModel::allocateArray ( int*& array, int initVal, @@ -166,16 +166,15 @@ void Foam::dataExchangeModel::allocateArray array[i] = initVal; } -void Foam::dataExchangeModel::destroy(int* array) const +void dataExchangeModel::destroy(int* array) const { - if (array == NULL) return; delete [] array; } //==== //==== // double * -void Foam::dataExchangeModel::allocateArray +void dataExchangeModel::allocateArray ( double*& array, double initVal, @@ -189,16 +188,14 @@ void Foam::dataExchangeModel::allocateArray array[i] = initVal; } -void Foam::dataExchangeModel::destroy(double* array) const +void dataExchangeModel::destroy(double* array) const { - if (array == NULL) return; - delete [] array; } //==== -bool Foam::dataExchangeModel::couple(int i) const +bool dataExchangeModel::couple(int i) const { bool coupleNow = false; if (doCoupleNow()) @@ -209,7 +206,7 @@ bool Foam::dataExchangeModel::couple(int i) const return coupleNow; } -scalar Foam::dataExchangeModel::timeStepFraction() const +scalar dataExchangeModel::timeStepFraction() const { //return fraction between previous coupling TS and actual TS //scalar DEMtime = DEMts_ * couplingInterval_; @@ -220,25 +217,25 @@ scalar Foam::dataExchangeModel::timeStepFraction() const return frac; } -int Foam::dataExchangeModel::getNumberOfParticles() const +int dataExchangeModel::getNumberOfParticles() const { Warning << "ask for nr of particles - which is not supported for this dataExchange model" << endl; return -1; } -int Foam::dataExchangeModel::getNumberOfClumps() const +int dataExchangeModel::getNumberOfClumps() const { Warning << "ask for nr of clumps - which is not supported for this dataExchange model" << endl; return -1; } -int Foam::dataExchangeModel::getNumberOfTypes() const +int dataExchangeModel::getNumberOfTypes() const { Warning << "ask for nr of types - which is not supported for this dataExchange model" << endl; return -1; } -double* Foam::dataExchangeModel::getTypeVol() const +double* dataExchangeModel::getTypeVol() const { Warning << "ask for type volume - which is not supported for this dataExchange model" << endl; return NULL; diff --git a/src/lagrangian/cfdemParticle/subModels/dataExchangeModel/dataExchangeModel/dataExchangeModel.H b/src/lagrangian/cfdemParticle/subModels/dataExchangeModel/dataExchangeModel/dataExchangeModel.H index a1c7b78c..e6098fb7 100755 --- a/src/lagrangian/cfdemParticle/subModels/dataExchangeModel/dataExchangeModel/dataExchangeModel.H +++ b/src/lagrangian/cfdemParticle/subModels/dataExchangeModel/dataExchangeModel/dataExchangeModel.H @@ -245,16 +245,16 @@ public: for (int i=0;iforce->cg()); } + void setCG() { particleCloud_.setCG(lmp->force->cg()); } }; diff --git a/src/lagrangian/cfdemParticle/subModels/dataExchangeModel/twoWayMany2Many/twoWayMany2Many.C b/src/lagrangian/cfdemParticle/subModels/dataExchangeModel/twoWayMany2Many/twoWayMany2Many.C index 5dcb7f3c..74d50569 100644 --- a/src/lagrangian/cfdemParticle/subModels/dataExchangeModel/twoWayMany2Many/twoWayMany2Many.C +++ b/src/lagrangian/cfdemParticle/subModels/dataExchangeModel/twoWayMany2Many/twoWayMany2Many.C @@ -62,7 +62,7 @@ twoWayMany2Many::twoWayMany2Many : dataExchangeModel(dict,sm), propsDict_(dict.subDict(typeName + "Props")), - pbm_(sm.mesh().boundaryMesh()), + pbm_(sm.mesh().boundaryMesh()), pData_(sm.mesh().globalData()), procPatches_(pData_.processorPatches()), procPatchIndices_(pData_.processorPatchIndices()), @@ -280,7 +280,7 @@ void twoWayMany2Many::giveData //============ // double ** -void Foam::twoWayMany2Many::allocateArray +void twoWayMany2Many::allocateArray ( double**& array, double initVal, @@ -295,7 +295,7 @@ void Foam::twoWayMany2Many::allocateArray array[i][j] = initVal; } -void Foam::twoWayMany2Many::allocateArray +void twoWayMany2Many::allocateArray ( double**& array, double initVal, @@ -310,14 +310,14 @@ void Foam::twoWayMany2Many::allocateArray array[i][j] = initVal; } -void inline Foam::twoWayMany2Many::destroy(double** array,int len) const +void inline twoWayMany2Many::destroy(double** array,int len) const { lmp->memory->destroy(array); } //============ // int ** -void Foam::twoWayMany2Many::allocateArray +void twoWayMany2Many::allocateArray ( int**& array, int initVal, @@ -332,7 +332,7 @@ void Foam::twoWayMany2Many::allocateArray array[i][j] = initVal; } -void Foam::twoWayMany2Many::allocateArray +void twoWayMany2Many::allocateArray ( int**& array, int initVal, @@ -347,14 +347,14 @@ void Foam::twoWayMany2Many::allocateArray array[i][j] = initVal; } -void inline Foam::twoWayMany2Many::destroy(int** array,int len) const +void inline twoWayMany2Many::destroy(int** array,int len) const { lmp->memory->destroy(array); } //============ // double * -void Foam::twoWayMany2Many::allocateArray(double*& array, double initVal, int length) const +void twoWayMany2Many::allocateArray(double*& array, double initVal, int length) const { int len = max(length,1); lmp->memory->grow(array, len, "m2m:dbl*"); @@ -362,14 +362,14 @@ void Foam::twoWayMany2Many::allocateArray(double*& array, double initVal, int le array[i] = initVal; } -void inline Foam::twoWayMany2Many::destroy(double* array) const +void inline twoWayMany2Many::destroy(double* array) const { lmp->memory->destroy(array); } //============== // int * -void Foam::twoWayMany2Many::allocateArray(int*& array, int initVal, int length) const +void twoWayMany2Many::allocateArray(int*& array, int initVal, int length) const { int len = max(length,1); lmp->memory->grow(array, len, "m2m:int*"); @@ -377,14 +377,14 @@ void Foam::twoWayMany2Many::allocateArray(int*& array, int initVal, int length) array[i] = initVal; } -void inline Foam::twoWayMany2Many::destroy(int* array) const +void inline twoWayMany2Many::destroy(int* array) const { lmp->memory->destroy(array); } //============== -bool Foam::twoWayMany2Many::couple(int i) const +bool twoWayMany2Many::couple(int i) const { bool coupleNow = false; if (i==0) @@ -443,19 +443,19 @@ bool Foam::twoWayMany2Many::couple(int i) const return coupleNow; } -int Foam::twoWayMany2Many::getNumberOfParticles() const +int twoWayMany2Many::getNumberOfParticles() const { return liggghts_get_maxtag(lmp); } -int Foam::twoWayMany2Many::getNumberOfClumps() const +int twoWayMany2Many::getNumberOfClumps() const { Warning << "Foam::twoWayMany2Many::getNumberOfClumps() - changes necessary here" << endl; //return liggghts_get_maxtag_ms(lmp); return 1; } -void Foam::twoWayMany2Many::syncIDs() const +void twoWayMany2Many::syncIDs() const { particleCloud_.clockM().start(5,"recv_DEM_ids"); @@ -575,7 +575,7 @@ void Foam::twoWayMany2Many::syncIDs() const particleCloud_.clockM().stop("setup_Comm"); } -void Foam::twoWayMany2Many::locateParticle(int* id_lammpsSync, bool id_lammps_alloc_flag) const +void twoWayMany2Many::locateParticle(int* id_lammpsSync, bool id_lammps_alloc_flag) const { #if defined(version21) @@ -669,7 +669,7 @@ void Foam::twoWayMany2Many::locateParticle(int* id_lammpsSync, bool id_lammps_al } particleCloud_.clockM().stop("locate_Stage1"); particleCloud_.clockM().start(8,"locate_Stage2"); - + PstreamBuffers pBufs(Pstream::nonBlocking); forAll(particleTransferID, i) @@ -802,7 +802,7 @@ void Foam::twoWayMany2Many::locateParticle(int* id_lammpsSync, bool id_lammps_al if (firstRun_) { int* id_foam_nowhere_all; - Foam::dataExchangeModel::allocateArray(id_foam_nowhere_all,1,nlocal_foam_lostAll); + dataExchangeModel::allocateArray(id_foam_nowhere_all,1,nlocal_foam_lostAll); MPI_Allreduce(id_foamLostAll, id_foam_nowhere_all, nlocal_foam_lostAll, MPI_INT, MPI_MIN, MPI_COMM_WORLD); int i = 0; @@ -820,7 +820,7 @@ void Foam::twoWayMany2Many::locateParticle(int* id_lammpsSync, bool id_lammps_al for (int k=0;k<3;k++) id_lammpsVec_[j*3+k] = id_lammpsVec_[(nlocal_lammps_-1)*3+k]; - + nlocal_lammps_ -= 1; break; } @@ -828,7 +828,7 @@ void Foam::twoWayMany2Many::locateParticle(int* id_lammpsSync, bool id_lammps_al } i++; } - Foam::dataExchangeModel::destroy(id_foam_nowhere_all); + dataExchangeModel::destroy(id_foam_nowhere_all); id_foam_nowhere_all = NULL; if (id_lammps_alloc_flag) destroy(id_lammps_); id_lammps_ = NULL; diff --git a/src/lagrangian/cfdemParticle/subModels/dataExchangeModel/twoWayMany2Many/twoWayMany2Many.H b/src/lagrangian/cfdemParticle/subModels/dataExchangeModel/twoWayMany2Many/twoWayMany2Many.H index 76d5de73..96d8fc25 100644 --- a/src/lagrangian/cfdemParticle/subModels/dataExchangeModel/twoWayMany2Many/twoWayMany2Many.H +++ b/src/lagrangian/cfdemParticle/subModels/dataExchangeModel/twoWayMany2Many/twoWayMany2Many.H @@ -114,7 +114,7 @@ private: mutable int *cellID_foam_; mutable double *pos_foam_; - const polyBoundaryMesh& pbm_; + const polyBoundaryMesh& pbm_; const globalMeshData& pData_; const labelList& procPatches_; const labelList& procPatchIndices_; @@ -209,7 +209,7 @@ public: void syncIDs() const; void locateParticle(int*, bool) const; word myType() const { return typeName; } - void setCG() const { particleCloud_.setCG(lmp->force->cg()); } + void setCG() { particleCloud_.setCG(lmp->force->cg()); } }; diff --git a/src/lagrangian/cfdemParticle/subModels/energyModel/heatTransferGunn/heatTransferGunn.C b/src/lagrangian/cfdemParticle/subModels/energyModel/heatTransferGunn/heatTransferGunn.C index 4959ca76..69175e97 100644 --- a/src/lagrangian/cfdemParticle/subModels/energyModel/heatTransferGunn/heatTransferGunn.C +++ b/src/lagrangian/cfdemParticle/subModels/energyModel/heatTransferGunn/heatTransferGunn.C @@ -134,19 +134,19 @@ heatTransferGunn::heatTransferGunn if (calcPartTempField_) { 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; + 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 (verbose_) { ReField_.writeOpt() = IOobject::AUTO_WRITE; - NuField_.writeOpt() = IOobject::AUTO_WRITE; - ReField_.write(); - NuField_.write(); + NuField_.writeOpt() = IOobject::AUTO_WRITE; + ReField_.write(); + NuField_.write(); } } @@ -155,10 +155,10 @@ heatTransferGunn::heatTransferGunn heatTransferGunn::~heatTransferGunn() { - delete partTemp_; - delete partHeatFlux_; - delete partRe_; - delete partNu_; + particleCloud_.dataExchangeM().destroy(partTemp_,1); + particleCloud_.dataExchangeM().destroy(partHeatFlux_,1); + particleCloud_.dataExchangeM().destroy(partRe_,1); + particleCloud_.dataExchangeM().destroy(partNu_,1); } // * * * * * * * * * * * * * * * private Member Functions * * * * * * * * * * * * * // @@ -168,7 +168,7 @@ 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(verbose_) { particleCloud_.dataExchangeM().allocateArray(partRe_,initVal,1); @@ -188,10 +188,10 @@ void heatTransferGunn::calcEnergyContribution() // get DEM data particleCloud_.dataExchangeM().getData(partTempName_,"scalar-atom",partTemp_); - + if(calcPartTempField_) - { - partTempField_.primitiveFieldRef() = 0.0; + { + partTempField_.primitiveFieldRef() = 0.0; particleCloud_.averagingM().resetWeightFields(); particleCloud_.averagingM().setScalarAverage ( @@ -202,10 +202,10 @@ void heatTransferGunn::calcEnergyContribution() NULL ); - volScalarField sumTp (particleCloud_.averagingM().UsWeightField() * partTempField_); - dimensionedScalar aveTemp("aveTemp",dimensionSet(0,0,0,1,0,0,0), gSum(sumTp) / particleCloud_.numberOfParticles()); - partRelTempField_ = (partTempField_ - aveTemp) / (aveTemp - partRefTemp_); - Info << "heatTransferGunn: average part. temp = " << aveTemp.value() << endl; + volScalarField sumTp (particleCloud_.averagingM().UsWeightField() * partTempField_); + dimensionedScalar aveTemp("aveTemp",dimensionSet(0,0,0,1,0,0,0), gSum(sumTp) / particleCloud_.numberOfParticles()); + partRelTempField_ = (partTempField_ - aveTemp) / (aveTemp - partRefTemp_); + Info << "heatTransferGunn: average part. temp = " << aveTemp.value() << endl; } #ifdef compre @@ -213,7 +213,7 @@ void heatTransferGunn::calcEnergyContribution() #else const volScalarField mufField = particleCloud_.turbulence().nu()*rho_; #endif - + // calc La based heat flux scalar voidfraction(1); @@ -251,9 +251,9 @@ void heatTransferGunn::calcEnergyContribution() Ufluid = U_[cellI]; Tfluid = tempField_[cellI]; } - + if (voidfraction < 0.01) - voidfraction = 0.01; + voidfraction = 0.01; // calc relative velocity Us = particleCloud_.velocity(index); @@ -264,7 +264,7 @@ void heatTransferGunn::calcEnergyContribution() Pr = max(SMALL, Cp_ * muf / kf0_); Nup = Nusselt(voidfraction, Rep, Pr); - + scalar h = kf0_ * Nup / ds; scalar As = ds * ds * M_PI; // surface area of sphere @@ -272,14 +272,14 @@ void heatTransferGunn::calcEnergyContribution() // calc convective heat flux [W] heatFlux(index, h, As, Tfluid); heatFluxCoeff(index, h, As); - - if(verbose_) - { - partRe_[index][0] = Rep; - partNu_[index][0] = Nup; - } - - if(particleCloud_.verbose() && index >=0 && index <2) + + if(verbose_) + { + partRe_[index][0] = Rep; + partNu_[index][0] = Nup; + } + + if(particleCloud_.verbose() && index >=0 && index <2) { Info << "partHeatFlux = " << partHeatFlux_[index][0] << endl; Info << "magUr = " << magUr << endl; @@ -308,27 +308,27 @@ void heatTransferGunn::calcEnergyContribution() if(verbose_) { ReField_.primitiveFieldRef() = 0.0; - NuField_.primitiveFieldRef() = 0.0; - particleCloud_.averagingM().resetWeightFields(); + NuField_.primitiveFieldRef() = 0.0; + particleCloud_.averagingM().resetWeightFields(); particleCloud_.averagingM().setScalarAverage ( ReField_, partRe_, particleCloud_.particleWeights(), - particleCloud_.averagingM().UsWeightField(), + particleCloud_.averagingM().UsWeightField(), NULL ); - particleCloud_.averagingM().resetWeightFields(); + particleCloud_.averagingM().resetWeightFields(); particleCloud_.averagingM().setScalarAverage ( NuField_, partNu_, particleCloud_.particleWeights(), - particleCloud_.averagingM().UsWeightField(), + particleCloud_.averagingM().UsWeightField(), NULL ); } - + // limit source term forAll(QPartFluid_,cellI) { @@ -336,13 +336,13 @@ void heatTransferGunn::calcEnergyContribution() if(mag(EuFieldInCell) > maxSource_ ) { - Pout << "limiting source term\n" << endl ; + Pout << "limiting source term\n" << endl ; QPartFluid_[cellI] = sign(EuFieldInCell) * maxSource_; } } - + QPartFluid_.correctBoundaryConditions(); - + giveData(0); } @@ -354,13 +354,10 @@ void heatTransferGunn::addEnergyContribution(volScalarField& Qsource) const scalar heatTransferGunn::Nusselt(scalar voidfraction, scalar Rep, scalar Pr) const { - scalar Nup(0.0); - Nup = (7 - 10 * voidfraction + 5 * voidfraction * voidfraction) * + return (7 - 10 * voidfraction + 5 * voidfraction * voidfraction) * (1 + 0.7 * Foam::pow(Rep,0.2) * Foam::pow(Pr,0.33)) + - (1.33 - 2.4 * voidfraction + 1.2 * voidfraction * voidfraction) * - Foam::pow(Rep,0.7) * Foam::pow(Pr,0.33); - - return Nup; + (1.33 - 2.4 * voidfraction + 1.2 * voidfraction * voidfraction) * + Foam::pow(Rep,0.7) * Foam::pow(Pr,0.33); } void heatTransferGunn::heatFlux(label index, scalar h, scalar As, scalar Tfluid) @@ -370,7 +367,7 @@ void heatTransferGunn::heatFlux(label index, scalar h, scalar As, scalar Tfluid) void heatTransferGunn::heatFluxCoeff(label index, scalar h, scalar As) { - //no heat transfer coefficient in explicit model + //no heat transfer coefficient in explicit model } void heatTransferGunn::giveData(int call) @@ -379,7 +376,7 @@ void heatTransferGunn::giveData(int call) { Info << "total convective particle-fluid heat flux [W] (Eulerian) = " << gSum(QPartFluid_*1.0*QPartFluid_.mesh().V()) << endl; - particleCloud_.dataExchangeM().giveData(partHeatFluxName_,"scalar-atom", partHeatFlux_); + particleCloud_.dataExchangeM().giveData(partHeatFluxName_,"scalar-atom", partHeatFlux_); } } // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // diff --git a/src/lagrangian/cfdemParticle/subModels/energyModel/heatTransferGunnImplicit/heatTransferGunnImplicit.C b/src/lagrangian/cfdemParticle/subModels/energyModel/heatTransferGunnImplicit/heatTransferGunnImplicit.C index a204f6ca..79941b51 100644 --- a/src/lagrangian/cfdemParticle/subModels/energyModel/heatTransferGunnImplicit/heatTransferGunnImplicit.C +++ b/src/lagrangian/cfdemParticle/subModels/energyModel/heatTransferGunnImplicit/heatTransferGunnImplicit.C @@ -58,7 +58,7 @@ heatTransferGunnImplicit::heatTransferGunnImplicit partHeatFluxCoeff_(NULL) { allocateMyArrays(); - + // no limiting necessary for implicit heat transfer maxSource_ = 1e30; } @@ -68,7 +68,7 @@ heatTransferGunnImplicit::heatTransferGunnImplicit heatTransferGunnImplicit::~heatTransferGunnImplicit() { - delete partHeatFluxCoeff_; + particleCloud_.dataExchangeM().destroy(partHeatFluxCoeff_,1); } // * * * * * * * * * * * * * * * private Member Functions * * * * * * * * * * * * * // @@ -87,7 +87,7 @@ void heatTransferGunnImplicit::calcEnergyContribution() heatTransferGunn::calcEnergyContribution(); QPartFluidCoeff_.primitiveFieldRef() = 0.0; - + particleCloud_.averagingM().setScalarSum ( QPartFluidCoeff_, @@ -97,7 +97,7 @@ void heatTransferGunnImplicit::calcEnergyContribution() ); QPartFluidCoeff_.primitiveFieldRef() /= -QPartFluidCoeff_.mesh().V(); - + // QPartFluidCoeff_.correctBoundaryConditions(); } @@ -123,7 +123,7 @@ void heatTransferGunnImplicit::giveData(int call) { //Info << "total convective particle-fluid heat flux [W] (Eulerian) = " << gSum(QPartFluid_*1.0*QPartFluid_.mesh().V()) << endl; - particleCloud_.dataExchangeM().giveData(partHeatFluxName_,"scalar-atom", partHeatFlux_); + particleCloud_.dataExchangeM().giveData(partHeatFluxName_,"scalar-atom", partHeatFlux_); } } @@ -133,7 +133,7 @@ void heatTransferGunnImplicit::postFlow() scalar Tfluid(0.0); scalar Tpart(0.0); interpolationCellPoint TInterpolator_(tempField_); - + for(int index = 0;index < particleCloud_.numberOfParticles(); ++index) { cellI = particleCloud_.cellIDs()[index][0]; @@ -146,12 +146,12 @@ void heatTransferGunnImplicit::postFlow() } else Tfluid = tempField_[cellI]; - + Tpart = partTemp_[index][0]; partHeatFlux_[index][0] = (Tfluid - Tpart) * partHeatFluxCoeff_[index][0]; } } - + giveData(1); } // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // diff --git a/src/lagrangian/cfdemParticle/subModels/energyModel/reactionHeat/reactionHeat.C b/src/lagrangian/cfdemParticle/subModels/energyModel/reactionHeat/reactionHeat.C index a3cf1635..5ca1c6c6 100644 --- a/src/lagrangian/cfdemParticle/subModels/energyModel/reactionHeat/reactionHeat.C +++ b/src/lagrangian/cfdemParticle/subModels/energyModel/reactionHeat/reactionHeat.C @@ -77,7 +77,7 @@ reactionHeat::reactionHeat reactionHeat::~reactionHeat() { - delete reactionHeat_; + particleCloud_.dataExchangeM().destroy(reactionHeat_,1); } // * * * * * * * * * * * * * * * private Member Functions * * * * * * * * * * * * * // diff --git a/src/lagrangian/cfdemParticle/subModels/forceModel/Archimedes/Archimedes.C b/src/lagrangian/cfdemParticle/subModels/forceModel/Archimedes/Archimedes.C index 50f5a7e4..c218844d 100644 --- a/src/lagrangian/cfdemParticle/subModels/forceModel/Archimedes/Archimedes.C +++ b/src/lagrangian/cfdemParticle/subModels/forceModel/Archimedes/Archimedes.C @@ -68,10 +68,10 @@ Archimedes::Archimedes { //Append the field names to be probed - particleCloud_.probeM().initialize(typeName, "archimedesF.logDat"); + particleCloud_.probeM().initialize(typeName, typeName+".logDat"); particleCloud_.probeM().vectorFields_.append("archimedesForce"); //first entry must the be the force particleCloud_.probeM().scalarFields_.append("Vp"); - particleCloud_.probeM().writeHeader(); + particleCloud_.probeM().writeHeader(); if (propsDict_.found("twoDimensional")) @@ -84,18 +84,20 @@ Archimedes::Archimedes setForceSubModels(propsDict_); // define switches which can be read from dict - forceSubM(0).setSwitchesList(0,true); // activate treatExplicit switch - forceSubM(0).setSwitchesList(1,true); // activate treatForceDEM switch + forceSubM(0).setSwitchesList(SW_TREAT_FORCE_EXPLICIT,true); // activate treatExplicit switch + forceSubM(0).setSwitchesList(SW_TREAT_FORCE_DEM,true); // activate treatForceDEM switch forceSubM(0).readSwitches(); - if (modelType_=="A" || modelType_=="Bfull"){ - if(!forceSubM(0).switches()[1]) // treatDEM != true + if (modelType_=="A" || modelType_=="Bfull") + { + if(!forceSubM(0).switches()[SW_TREAT_FORCE_DEM]) // treatDEM != true { Warning << "Usually model type A and Bfull need Archimedes only on DEM side only (treatForceDEM=true)! are you sure about your settings?" << endl; } } - if (modelType_=="B"){ - if(forceSubM(0).switches()[1]) // treatDEM = true + else if (modelType_=="B") + { + if(forceSubM(0).switches()[SW_TREAT_FORCE_DEM]) // treatDEM = true { Warning << "Usually model type B needs Archimedes on CFD and DEM side (treatForceDEM=false)! are you sure about your settings?" << endl; } @@ -119,12 +121,12 @@ void Archimedes::setForce() const #include "setupProbeModel.H" - for(int index = 0;index < particleCloud_.numberOfParticles(); ++index) + for(int index = 0; index < particleCloud_.numberOfParticles(); ++index) { //if(mask[index][0]) //{ label cellI = particleCloud_.cellIDs()[index][0]; - force=vector::zero; + force = vector::zero; if (cellI > -1) // particle Found { @@ -133,7 +135,9 @@ void Archimedes::setForce() const scalar r = particleCloud_.radius(index); force = -g_.value()*forceSubM(0).rhoField()[cellI]*r*r*M_PI; // circle area Warning << "Archimedes::setForce() : this functionality is not tested!" << endl; - }else{ + } + else + { force = -g_.value()*forceSubM(0).rhoField()[cellI]*particleCloud_.particleVolume(index); } diff --git a/src/lagrangian/cfdemParticle/subModels/forceModel/ArchimedesIB/ArchimedesIB.C b/src/lagrangian/cfdemParticle/subModels/forceModel/ArchimedesIB/ArchimedesIB.C index 46c016b6..4bafc636 100644 --- a/src/lagrangian/cfdemParticle/subModels/forceModel/ArchimedesIB/ArchimedesIB.C +++ b/src/lagrangian/cfdemParticle/subModels/forceModel/ArchimedesIB/ArchimedesIB.C @@ -69,9 +69,9 @@ ArchimedesIB::ArchimedesIB g_(sm.mesh().lookupObject (gravityFieldName_)) { //Append the field names to be probed - particleCloud_.probeM().initialize(typeName, "archimedesIBF.logDat"); + particleCloud_.probeM().initialize(typeName, typeName+".logDat"); particleCloud_.probeM().vectorFields_.append("archimedesIBForce"); //first entry must the be the force - particleCloud_.probeM().writeHeader(); + particleCloud_.probeM().writeHeader(); if (propsDict_.found("twoDimensional")) { @@ -83,12 +83,12 @@ ArchimedesIB::ArchimedesIB setForceSubModels(propsDict_); // define switches which can be read from dict - forceSubM(0).setSwitchesList(0,true); // activate treatExplicit switch + forceSubM(0).setSwitchesList(SW_TREAT_FORCE_EXPLICIT,true); // activate treatExplicit switch // read those switches defined above, if provided in dict forceSubM(0).readSwitches(); - forceSubM(0).setSwitches(1,true); // treatDEM = true + forceSubM(0).setSwitches(SW_TREAT_FORCE_DEM,true); // treatDEM = true Info << "accounting for Archimedes only on DEM side!" << endl; particleCloud_.checkCG(true); @@ -109,7 +109,7 @@ void ArchimedesIB::setForce() const #include "setupProbeModel.H" - for(int index = 0;index < particleCloud_.numberOfParticles(); ++index) + for(int index = 0; index < particleCloud_.numberOfParticles(); ++index) { //if(mask[index][0]) //{ @@ -120,8 +120,8 @@ void ArchimedesIB::setForce() const if (cellI > -1) // particle Found { //force += -g_.value()*forceSubM(0).rhoField()[cellI]*forceSubM(0).rhoField().mesh().V()[cellI]*(1-particleCloud_.voidfractions()[index][subCell]);//mod by alice - force += -g_.value()*forceSubM(0).rhoField()[cellI]*particleCloud_.mesh().V()[cellI]*(1-voidfractions_[cellI]);//mod by alice - } + force += -g_.value()*forceSubM(0).rhoField()[cellI]*particleCloud_.mesh().V()[cellI]*(1-voidfractions_[cellI]);//mod by alice + } } //Set value fields and write the probe diff --git a/src/lagrangian/cfdemParticle/subModels/forceModel/BeetstraDrag/BeetstraDrag.C b/src/lagrangian/cfdemParticle/subModels/forceModel/BeetstraDrag/BeetstraDrag.C index 779a7a42..58677ad3 100644 --- a/src/lagrangian/cfdemParticle/subModels/forceModel/BeetstraDrag/BeetstraDrag.C +++ b/src/lagrangian/cfdemParticle/subModels/forceModel/BeetstraDrag/BeetstraDrag.C @@ -59,7 +59,7 @@ BeetstraDrag::BeetstraDrag scaleDrag_(1.) { //Append the field names to be probed - particleCloud_.probeM().initialize(typeName, "BeetstraDrag.logDat"); + particleCloud_.probeM().initialize(typeName, typeName+".logDat"); particleCloud_.probeM().vectorFields_.append("dragForce"); //first entry must be the force particleCloud_.probeM().vectorFields_.append("Urel"); particleCloud_.probeM().scalarFields_.append("Rep"); @@ -70,11 +70,11 @@ BeetstraDrag::BeetstraDrag // init force sub model setForceSubModels(propsDict_); // define switches which can be read from dict - forceSubM(0).setSwitchesList(0,true); // activate treatExplicit switch - forceSubM(0).setSwitchesList(2,true); // activate implDEM switch - forceSubM(0).setSwitchesList(3,true); // activate search for verbose switch - forceSubM(0).setSwitchesList(4,true); // activate search for interpolate switch - forceSubM(0).setSwitchesList(8,true); // activate scalarViscosity switch + forceSubM(0).setSwitchesList(SW_TREAT_FORCE_EXPLICIT,true); // activate treatExplicit switch + forceSubM(0).setSwitchesList(SW_IMPL_FORCE_DEM,true); // activate implDEM switch + forceSubM(0).setSwitchesList(SW_VERBOSE,true); // activate search for verbose switch + forceSubM(0).setSwitchesList(SW_INTERPOLATION,true); // activate search for interpolate switch + forceSubM(0).setSwitchesList(SW_SCALAR_VISCOSITY,true); // activate scalarViscosity switch forceSubM(0).readSwitches(); particleCloud_.checkCG(true); @@ -97,8 +97,11 @@ BeetstraDrag::~BeetstraDrag() void BeetstraDrag::setForce() const { if (scaleDia_ > 1) + { Info << "Beetstra using scale = " << scaleDia_ << endl; - else if (particleCloud_.cg() > 1){ + } + else if (particleCloud_.cg() > 1) + { scaleDia_=particleCloud_.cg(); Info << "Beetstra using scale from liggghts cg = " << scaleDia_ << endl; } @@ -125,13 +128,13 @@ void BeetstraDrag::setForce() const vector dragExplicit(0,0,0); scalar dragCoefficient(0); - + interpolationCellPoint voidfractionInterpolator_(voidfraction_); interpolationCellPoint UInterpolator_(U_); #include "setupProbeModel.H" - for(int index = 0;index < particleCloud_.numberOfParticles(); ++index) + for(int index = 0; index < particleCloud_.numberOfParticles(); ++index) { cellI = particleCloud_.cellIDs()[index][0]; drag = vector(0,0,0); @@ -163,7 +166,7 @@ void BeetstraDrag::setForce() const Ur = Ufluid-Us; magUr = mag(Ur); ds = 2*particleCloud_.radius(index); - ds_scaled = ds/scaleDia_; + ds_scaled = ds/scaleDia_; rho = rhoField[cellI]; nuf = nufField[cellI]; @@ -171,7 +174,7 @@ void BeetstraDrag::setForce() const 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)) + diff --git a/src/lagrangian/cfdemParticle/subModels/forceModel/BeetstraDrag/BeetstraDrag.H b/src/lagrangian/cfdemParticle/subModels/forceModel/BeetstraDrag/BeetstraDrag.H index a36ab894..f6e079ab 100644 --- a/src/lagrangian/cfdemParticle/subModels/forceModel/BeetstraDrag/BeetstraDrag.H +++ b/src/lagrangian/cfdemParticle/subModels/forceModel/BeetstraDrag/BeetstraDrag.H @@ -51,7 +51,7 @@ private: const volScalarField& voidfraction_; - word UsFieldName_; + word UsFieldName_; const volVectorField& UsField_; diff --git a/src/lagrangian/cfdemParticle/subModels/forceModel/DiFeliceDrag/DiFeliceDrag.C b/src/lagrangian/cfdemParticle/subModels/forceModel/DiFeliceDrag/DiFeliceDrag.C index 29bc0ef7..c15da1e1 100644 --- a/src/lagrangian/cfdemParticle/subModels/forceModel/DiFeliceDrag/DiFeliceDrag.C +++ b/src/lagrangian/cfdemParticle/subModels/forceModel/DiFeliceDrag/DiFeliceDrag.C @@ -34,8 +34,6 @@ Description #include "DiFeliceDrag.H" #include "addToRunTimeSelectionTable.H" -//#include - // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // namespace Foam @@ -74,12 +72,12 @@ DiFeliceDrag::DiFeliceDrag scaleDrag_(1.) { //Append the field names to be probed - particleCloud_.probeM().initialize(typeName, "diFeliceDrag.logDat"); - particleCloud_.probeM().vectorFields_.append("dragForce"); //first entry must the be the force - particleCloud_.probeM().vectorFields_.append("Urel"); //other are debug - particleCloud_.probeM().scalarFields_.append("Rep"); //other are debug - particleCloud_.probeM().scalarFields_.append("Cd"); //other are debug - particleCloud_.probeM().scalarFields_.append("voidfraction"); //other are debug + particleCloud_.probeM().initialize(typeName, typeName+".logDat"); + particleCloud_.probeM().vectorFields_.append("dragForce"); // first entry must the be the force + particleCloud_.probeM().vectorFields_.append("Urel"); // other are debug + particleCloud_.probeM().scalarFields_.append("Rep"); // other are debug + particleCloud_.probeM().scalarFields_.append("Cd"); // other are debug + particleCloud_.probeM().scalarFields_.append("voidfraction"); // other are debug particleCloud_.probeM().writeHeader(); particleCloud_.checkCG(true); @@ -92,11 +90,11 @@ DiFeliceDrag::DiFeliceDrag setForceSubModels(propsDict_); // define switches which can be read from dict - forceSubM(0).setSwitchesList(0,true); // activate treatExplicit switch - forceSubM(0).setSwitchesList(2,true); // activate implDEM switch - forceSubM(0).setSwitchesList(3,true); // activate search for verbose switch - forceSubM(0).setSwitchesList(4,true); // activate search for interpolate switch - forceSubM(0).setSwitchesList(8,true); // activate scalarViscosity switch + forceSubM(0).setSwitchesList(SW_TREAT_FORCE_EXPLICIT,true); // activate treatExplicit switch + forceSubM(0).setSwitchesList(SW_IMPL_FORCE_DEM,true); // activate implDEM switch + forceSubM(0).setSwitchesList(SW_VERBOSE,true); // activate search for verbose switch + forceSubM(0).setSwitchesList(SW_INTERPOLATION,true); // activate search for interpolate switch + forceSubM(0).setSwitchesList(SW_SCALAR_VISCOSITY,true); // activate scalarViscosity switch // read those switches defined above, if provided in dict forceSubM(0).readSwitches(); @@ -114,15 +112,18 @@ DiFeliceDrag::~DiFeliceDrag() void DiFeliceDrag::setForce() const { if (scaleDia_ > 1) + { Info << "DiFeliceDrag using scale = " << scaleDia_ << endl; - else if (particleCloud_.cg() > 1){ + } + else if (particleCloud_.cg() > 1) + { scaleDia_=particleCloud_.cg(); Info << "DiFeliceDrag using scale from liggghts cg = " << scaleDia_ << endl; } const volScalarField& nufField = forceSubM(0).nuField(); const volScalarField& rhoField = forceSubM(0).rhoField(); - + vector position(0,0,0); scalar voidfraction(1); vector Ufluid(0,0,0); @@ -144,7 +145,7 @@ void DiFeliceDrag::setForce() const #include "setupProbeModel.H" - for(int index = 0;index < particleCloud_.numberOfParticles(); index++) + for(int index = 0; index < particleCloud_.numberOfParticles(); ++index) { cellI = particleCloud_.cellIDs()[index][0]; drag = vector(0,0,0); @@ -159,7 +160,8 @@ void DiFeliceDrag::setForce() const position = particleCloud_.position(index); voidfraction = voidfractionInterpolator_.interpolate(position,cellI); Ufluid = UInterpolator_.interpolate(position,cellI); - }else + } + else { voidfraction = voidfraction_[cellI]; Ufluid = U_[cellI]; @@ -189,8 +191,8 @@ void DiFeliceDrag::setForce() const // calc particle's drag dragCoefficient = 0.125*Cd*rho *M_PI - *ds*ds - *scaleDia_ + *ds*ds + *scaleDia_ *pow(voidfraction,(2-Xi))*magUr *scaleDrag_; if (modelType_=="B") diff --git a/src/lagrangian/cfdemParticle/subModels/forceModel/DiFeliceDrag/DiFeliceDrag.H b/src/lagrangian/cfdemParticle/subModels/forceModel/DiFeliceDrag/DiFeliceDrag.H index a043b28e..b92e8c0f 100644 --- a/src/lagrangian/cfdemParticle/subModels/forceModel/DiFeliceDrag/DiFeliceDrag.H +++ b/src/lagrangian/cfdemParticle/subModels/forceModel/DiFeliceDrag/DiFeliceDrag.H @@ -68,7 +68,7 @@ private: const volScalarField& voidfraction_; - word UsFieldName_; + word UsFieldName_; const volVectorField& UsField_; // the average particle velocity field (for implicit/expliti force split) diff --git a/src/lagrangian/cfdemParticle/subModels/forceModel/Fines/ErgunStatFines.C b/src/lagrangian/cfdemParticle/subModels/forceModel/Fines/ErgunStatFines.C index 17fffd4e..e229c15e 100644 --- a/src/lagrangian/cfdemParticle/subModels/forceModel/Fines/ErgunStatFines.C +++ b/src/lagrangian/cfdemParticle/subModels/forceModel/Fines/ErgunStatFines.C @@ -66,7 +66,7 @@ ErgunStatFines::ErgunStatFines switchingVoidfraction_(0.8) { //Append the field names to be probed - particleCloud_.probeM().initialize(typeName, "ErgunStatFines.logDat"); + particleCloud_.probeM().initialize(typeName, typeName+".logDat"); particleCloud_.probeM().vectorFields_.append("dragForce"); //first entry must be the force particleCloud_.probeM().vectorFields_.append("Urel"); particleCloud_.probeM().scalarFields_.append("Rep"); @@ -77,26 +77,26 @@ ErgunStatFines::ErgunStatFines // init force sub model setForceSubModels(propsDict_); // define switches which can be read from dict - forceSubM(0).setSwitchesList(0,true); // activate treatExplicit switch - forceSubM(0).setSwitchesList(2,true); // activate implDEM switch - forceSubM(0).setSwitchesList(3,true); // activate search for verbose switch - forceSubM(0).setSwitchesList(4,true); // activate search for interpolate switch - forceSubM(0).setSwitchesList(8,true); // activate scalarViscosity switch + forceSubM(0).setSwitchesList(SW_TREAT_FORCE_EXPLICIT,true); // activate treatExplicit switch + forceSubM(0).setSwitchesList(SW_IMPL_FORCE_DEM,true); // activate implDEM switch + forceSubM(0).setSwitchesList(SW_VERBOSE,true); // activate search for verbose switch + forceSubM(0).setSwitchesList(SW_INTERPOLATION,true); // activate search for interpolate switch + forceSubM(0).setSwitchesList(SW_SCALAR_VISCOSITY,true); // activate scalarViscosity switch forceSubM(0).readSwitches(); particleCloud_.checkCG(true); if (propsDict_.found("scale")) - scaleDia_=scalar(readScalar(propsDict_.lookup("scale"))); + scaleDia_ = scalar(readScalar(propsDict_.lookup("scale"))); if (propsDict_.found("scaleDrag")) - scaleDrag_=scalar(readScalar(propsDict_.lookup("scaleDrag"))); + scaleDrag_ = scalar(readScalar(propsDict_.lookup("scaleDrag"))); if (propsDict_.found("switchingVoidfraction")) switchingVoidfraction_ = readScalar(propsDict_.lookup("switchingVoidfraction")); dictionary SauterDict(dict.subDict("dSauterProps")); if (SauterDict.found("scaleDist")) - scaleDist_=scalar(readScalar(SauterDict.lookup("scaleDist"))); - + scaleDist_ = scalar(readScalar(SauterDict.lookup("scaleDist"))); + } @@ -117,7 +117,9 @@ scalar ErgunStatFines::dSauter(label cellI) const void ErgunStatFines::setForce() const { if (scaleDia_ > 1) + { Info << "ErgunStatFines using scale = " << scaleDia_ << endl; + } else if (particleCloud_.cg() > 1) { scaleDia_=particleCloud_.cg(); @@ -131,7 +133,7 @@ void ErgunStatFines::setForce() const scalar voidfraction(1); vector Ufluid(0,0,0); vector drag(0,0,0); - label cellI=0; + label cellI = 0; vector Us(0,0,0); vector Ur(0,0,0); @@ -144,13 +146,13 @@ void ErgunStatFines::setForce() const scalar alphaPartEff(0); scalar CdMagUrLag(0); //Cd of the very particle - scalar betaP(0); //momentum exchange of the very particle + scalar betaP(0); //momentum exchange of the very particle vector dragExplicit(0,0,0); scalar dragCoefficient(0); scalar scaleDia3 = scaleDia_*scaleDia_*scaleDia_; - + interpolationCellPoint voidfractionInterpolator_(voidfraction_); interpolationCellPoint UInterpolator_(U_); @@ -158,15 +160,15 @@ void ErgunStatFines::setForce() const if(forceSubM(0).verbose()) Info << "Entering force loop of ErgunStatFines.\n" << endl; - - for(int index = 0;index < particleCloud_.numberOfParticles(); ++index) - { + + for(int index = 0; index < particleCloud_.numberOfParticles(); ++index) + { cellI = particleCloud_.cellIDs()[index][0]; drag = vector(0,0,0); dragExplicit = vector(0,0,0); betaP = 0; - Ufluid =vector(0,0,0); - voidfraction=0; + Ufluid = vector(0,0,0); + voidfraction = 0; dragCoefficient = 0; if (cellI > -1) // particle found @@ -174,34 +176,34 @@ void ErgunStatFines::setForce() const if( forceSubM(0).interpolation() ) { - position = particleCloud_.position(index); + position = particleCloud_.position(index); voidfraction = voidfractionInterpolator_.interpolate(position,cellI); - Ufluid = UInterpolator_.interpolate(position,cellI); + Ufluid = UInterpolator_.interpolate(position,cellI); } else { - voidfraction = voidfraction_[cellI]; + voidfraction = voidfraction_[cellI]; Ufluid = U_[cellI]; } // ensure voidfraction to be meaningful // problems could arise from interpolation or empty cells - if(voidfraction>0.999) + if(voidfraction > 0.999) voidfraction = 0.999; - else if(voidfraction<0.05) + else if(voidfraction < 0.05) voidfraction = 0.05; Us = particleCloud_.velocity(index); Ur = Ufluid-Us; magUr = mag(Ur); - dSauterMix = dSauterMix_[cellI]; - ds = 2*particleCloud_.radius(index); + dSauterMix = dSauterMix_[cellI]; + ds = 2*particleCloud_.radius(index); rho = rhoField[cellI]; nuf = nufField[cellI]; - Rep=0.0; - alphaPartEff = 1.0 - voidfraction + alphaSt_[cellI] + SMALL; + Rep = 0.0; + alphaPartEff = 1.0 - voidfraction + alphaSt_[cellI] + SMALL; // calc particle's drag coefficient (i.e., Force per unit slip velocity and per m³ PARTICLE) if(voidfraction > switchingVoidfraction_) //dilute, no static hold-up present @@ -210,7 +212,7 @@ void ErgunStatFines::setForce() const CdMagUrLag = (24.0*nuf/(dSauterMix*voidfraction)) //1/magUr missing here, but compensated in expression for betaP! *(scalar(1.0)+0.15*Foam::pow(Rep, 0.687)); - betaP = 0.75* alphaPartEff * ( + betaP = 0.75* alphaPartEff * ( rho*voidfraction*CdMagUrLag / (dSauterMix*Foam::pow(voidfraction,2.65)) @@ -224,12 +226,12 @@ void ErgunStatFines::setForce() const (1.75 * magUr * rho * alphaPartEff) /((dSauterMix*phi_)); } - + // calc particle's drag betaP /= (1-alphaPartEff); dragCoefficient = M_PI/6 * ds/scaleDia_ * ds/scaleDia_ * dSauter(cellI) * voidfraction / (1 - voidfraction) * betaP * scaleDrag_; dragCoefficient *= scaleDia3; - if (modelType_=="B") + if (modelType_ == "B") dragCoefficient /= voidfraction; drag = dragCoefficient * Ur; @@ -237,7 +239,7 @@ void ErgunStatFines::setForce() const // explicitCorr forceSubM(0).explicitCorr(drag,dragExplicit,dragCoefficient,Ufluid,U_[cellI],Us,UsField_[cellI],forceSubM(0).verbose()); - if(forceSubM(0).verbose() && index >=0 && index <2) + if(forceSubM(0).verbose() && index >= 0 && index < 2) { Pout << "cellI = " << cellI << endl; Pout << "index = " << index << endl; @@ -271,7 +273,7 @@ void ErgunStatFines::setForce() const forceSubM(0).partToArray(index,drag,dragExplicit,Ufluid,dragCoefficient); }// end loop particles - + if(forceSubM(0).verbose()) Pout << "Leaving force loop of ErgunStatFines.\n" << endl; } diff --git a/src/lagrangian/cfdemParticle/subModels/forceModel/Fines/ErgunStatFines.H b/src/lagrangian/cfdemParticle/subModels/forceModel/Fines/ErgunStatFines.H index 1b33e30e..00bb186a 100644 --- a/src/lagrangian/cfdemParticle/subModels/forceModel/Fines/ErgunStatFines.H +++ b/src/lagrangian/cfdemParticle/subModels/forceModel/Fines/ErgunStatFines.H @@ -50,18 +50,18 @@ private: word voidfractionFieldName_; const volScalarField& voidfraction_; - + const volScalarField& dSauter_; - + const volScalarField& dSauterMix_; - + const volScalarField& alphaP_; - + const volScalarField& alphaSt_; const scalar phi_; - word UsFieldName_; + word UsFieldName_; const volVectorField& UsField_; // the average particle velocity field diff --git a/src/lagrangian/cfdemParticle/subModels/forceModel/Fines/FanningDynFines.C b/src/lagrangian/cfdemParticle/subModels/forceModel/Fines/FanningDynFines.C index 997cfdf4..29cbdd36 100644 --- a/src/lagrangian/cfdemParticle/subModels/forceModel/Fines/FanningDynFines.C +++ b/src/lagrangian/cfdemParticle/subModels/forceModel/Fines/FanningDynFines.C @@ -65,16 +65,16 @@ FanningDynFines::FanningDynFines // init force sub model setForceSubModels(propsDict_); // define switches which can be read from dict - forceSubM(0).setSwitchesList(0,true); // activate treatExplicit switch - forceSubM(0).setSwitchesList(3,true); // activate search for verbose switch + forceSubM(0).setSwitchesList(SW_TREAT_FORCE_EXPLICIT,true); // activate treatExplicit switch + forceSubM(0).setSwitchesList(SW_VERBOSE,true); // activate search for verbose switch forceSubM(0).readSwitches(); - forceSubM(0).setSwitches(0,true); + forceSubM(0).setSwitches(SW_TREAT_FORCE_EXPLICIT,true); particleCloud_.checkCG(true); if (propsDict_.found("scale")) - scaleDia_=scalar(readScalar(propsDict_.lookup("scale"))); + scaleDia_ = scalar(readScalar(propsDict_.lookup("scale"))); if (propsDict_.found("scaleDrag")) - scaleDrag_=scalar(readScalar(propsDict_.lookup("scaleDrag"))); + scaleDrag_ = scalar(readScalar(propsDict_.lookup("scaleDrag"))); } @@ -91,18 +91,20 @@ void FanningDynFines::setForce() const { if(forceSubM(0).verbose()) Info << "Entering force loop of FanningDynFines.\n" << endl; - + if (scaleDia_ > 1) + { Info << "FanningDynFines using scale = " << scaleDia_ << endl; + } else if (particleCloud_.cg() > 1) { scaleDia_=particleCloud_.cg(); Info << "FanningDynFines using scale from liggghts cg = " << scaleDia_ << endl; } - + vector UDyn(0,0,0); vector drag(0,0,0); - label cellI=0; + label cellI = 0; vector Us(0,0,0); vector Ur(0,0,0); @@ -111,11 +113,11 @@ void FanningDynFines::setForce() const scalar scaleDia3 = scaleDia_*scaleDia_*scaleDia_; scalar dragCoefficient(0); - + #include "setupProbeModel.H" - for(int index = 0;index < particleCloud_.numberOfParticles(); ++index) + for(int index = 0; index < particleCloud_.numberOfParticles(); ++index) { cellI = particleCloud_.cellIDs()[index][0]; drag = vector(0,0,0); @@ -127,14 +129,14 @@ void FanningDynFines::setForce() const UDyn = UDyn_[cellI]; Us = UsField_[cellI]; Ur = UDyn-Us; - ds = 2*particleCloud_.radius(index); - ds_scaled = ds/scaleDia_; - + ds = 2 * particleCloud_.radius(index); + ds_scaled = ds/scaleDia_; + dragCoefficient = FanningCoeff_[cellI]; // calc particle's drag dragCoefficient *= M_PI/6 * ds_scaled * ds_scaled / alphaP_[cellI] * dSauter_[cellI] * scaleDia3 * scaleDrag_; - if (modelType_=="B") + if (modelType_ == "B") dragCoefficient /= voidfraction_[cellI]; drag = dragCoefficient * Ur; @@ -143,8 +145,8 @@ void FanningDynFines::setForce() const // write particle based data to global array forceSubM(0).partToArray(index,drag,vector::zero); } - - if(forceSubM(0).verbose()) + + if (forceSubM(0).verbose()) Info << "Leaving force loop of FanningDynFines.\n" << endl; } // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // diff --git a/src/lagrangian/cfdemParticle/subModels/forceModel/Fines/FanningDynFines.H b/src/lagrangian/cfdemParticle/subModels/forceModel/Fines/FanningDynFines.H index 195de6c6..090ece60 100644 --- a/src/lagrangian/cfdemParticle/subModels/forceModel/Fines/FanningDynFines.H +++ b/src/lagrangian/cfdemParticle/subModels/forceModel/Fines/FanningDynFines.H @@ -41,7 +41,7 @@ class FanningDynFines public forceModel { private: - + dictionary propsDict_; word velFieldName_; @@ -52,22 +52,22 @@ private: const volScalarField& voidfraction_; - word UsFieldName_; + word UsFieldName_; const volVectorField& UsField_; - + const volVectorField& UDyn_; - + const volScalarField& FanningCoeff_; - + const volScalarField& alphaP_; - + const volScalarField& dSauter_; mutable scalar scaleDia_; mutable scalar scaleDrag_; - + public: diff --git a/src/lagrangian/cfdemParticle/subModels/forceModel/Fines/Fines.H b/src/lagrangian/cfdemParticle/subModels/forceModel/Fines/Fines.H index a6ff8d71..06153a63 100644 --- a/src/lagrangian/cfdemParticle/subModels/forceModel/Fines/Fines.H +++ b/src/lagrangian/cfdemParticle/subModels/forceModel/Fines/Fines.H @@ -41,7 +41,7 @@ class Fines public forceModel { private: - mutable FinesFields finesFields_; + mutable FinesFields finesFields_; public: diff --git a/src/lagrangian/cfdemParticle/subModels/forceModel/Fines/FinesFields.C b/src/lagrangian/cfdemParticle/subModels/forceModel/Fines/FinesFields.C index 8f6b3b30..1147bbe6 100644 --- a/src/lagrangian/cfdemParticle/subModels/forceModel/Fines/FinesFields.C +++ b/src/lagrangian/cfdemParticle/subModels/forceModel/Fines/FinesFields.C @@ -56,7 +56,7 @@ FinesFields::FinesFields p_(sm.mesh().lookupObject (pFieldName_)), rhoGFieldName_(propsDict_.lookupOrDefault("rhoGFieldName","rho")), rhoG_(sm.mesh().lookupObject (rhoGFieldName_)), - dSauter_(sm.mesh().lookupObject ("dSauter")), + dSauter_(sm.mesh().lookupObject ("dSauter")), alphaG_ ( IOobject ( @@ -126,7 +126,7 @@ FinesFields::FinesFields ), sm.mesh(), dimensionedScalar("zero", dimensionSet(0,1,0,0,0), 0), - "zeroGradient" + "zeroGradient" ), DragCoeff_ ( IOobject @@ -151,7 +151,7 @@ FinesFields::FinesFields ), sm.mesh(), dimensionedScalar("zero", dimensionSet(0,1,0,0,0), 0), - "zeroGradient" + "zeroGradient" ), FanningCoeff_ ( IOobject @@ -199,8 +199,8 @@ FinesFields::FinesFields IOobject::NO_WRITE ), sm.mesh(), - dimensionedScalar("zero", dimensionSet(1,0,-1,0,0), 0) - //dimensionedVector("zero", dimensionSet(1,-2,-1,0,0), vector::zero) + dimensionedScalar("zero", dimensionSet(1,0,-1,0,0), 0) + //dimensionedVector("zero", dimensionSet(1,-2,-1,0,0), vector::zero) ), uDyn_ ( IOobject @@ -243,12 +243,12 @@ FinesFields::FinesFields if (propsDict_.found("rhoFine")) rhoFine_.value()=readScalar(propsDict_.lookup ("rhoFine")); else - FatalError <<"Please specify rhoFine.\n" << abort(FatalError); + FatalError <<"Please specify rhoFine.\n" << abort(FatalError); if (propsDict_.found("nuAve")) nuAve_.value()=readScalar(propsDict_.lookup ("nuAve")); if (propsDict_.found("alphaDynMax")) alphaDynMax_=readScalar(propsDict_.lookup ("alphaDynMax")); - + if(verbose_) { alphaG_.writeOpt() = IOobject::AUTO_WRITE; @@ -305,143 +305,144 @@ void FinesFields::update() } -void FinesFields::calcSource() +void FinesFields::calcSource() { - Sds_.primitiveFieldRef()=0; + Sds_.primitiveFieldRef() = 0; deltaAlpha_.primitiveFieldRef() = 0.0; scalar f(0.0); scalar critpore(0.0); scalar dmean(0.0); scalar d1(0.0); scalar d2(0.0); + forAll(Sds_,cellI) { // calculate everything in units auf dSauter critpore = nCrit_*dFine_.value()/dSauter_[cellI]; - // pore size from hydraulic radius - dmean = 2 * (1 - alphaP_[cellI]) / ( (1 + poresizeWidth_*poresizeWidth_/3) * 3 * alphaP_[cellI] ); - // Sweeney and Martin, Acta Materialia 51 (2003): ratio of hydraulic to pore throat radius - dmean /= ratioHydraulicPore_; - d1 = dmean * (1 - poresizeWidth_); - d2 = dmean * (1 + poresizeWidth_); - + // pore size from hydraulic radius + dmean = 2 * (1 - alphaP_[cellI]) / ( (1 + poresizeWidth_*poresizeWidth_/3) * 3 * alphaP_[cellI] ); + // Sweeney and Martin, Acta Materialia 51 (2003): ratio of hydraulic to pore throat radius + dmean /= ratioHydraulicPore_; + d1 = dmean * (1 - poresizeWidth_); + d2 = dmean * (1 + poresizeWidth_); + f = (critpore*critpore*critpore - d1 * d1 * d1) / (d2 * d2 * d2 - d1 * d1 * d1); - if (f<0) - { - f=0.0; - } - else if (f>1.0) - { - f=1.0; + if (f < 0) + { + f = 0.0; } - - // at this point, voidfraction is still calculated from the true particle sizes - deltaAlpha_[cellI] = f * (alphaMax_ - alphaP_[cellI]) - alphaSt_[cellI]; - // too much volume occupied: release it (50% per time step) - if (deltaAlpha_[cellI] < 0.0) - { - Sds_[cellI] = 0.5*deltaAlpha_[cellI]; - } - // volume too occupy available: deposit at most 80% of dyn hold up - else if (depRate_ * deltaAlpha_[cellI] > 0.8 * alphaDyn_[cellI]) - { - Sds_[cellI] = 0.8 * alphaDyn_[cellI]; - } - else - { - Sds_[cellI] = depRate_ * deltaAlpha_[cellI]; - } - } + else if (f > 1.0) + { + f = 1.0; + } + + // at this point, voidfraction is still calculated from the true particle sizes + deltaAlpha_[cellI] = f * (alphaMax_ - alphaP_[cellI]) - alphaSt_[cellI]; + // too much volume occupied: release it (50% per time step) + if (deltaAlpha_[cellI] < 0.0) + { + Sds_[cellI] = 0.5*deltaAlpha_[cellI]; + } + // volume too occupy available: deposit at most 80% of dyn hold up + else if (depRate_ * deltaAlpha_[cellI] > 0.8 * alphaDyn_[cellI]) + { + Sds_[cellI] = 0.8 * alphaDyn_[cellI]; + } + else + { + Sds_[cellI] = depRate_ * deltaAlpha_[cellI]; + } + } } -void FinesFields::integrateFields() +void FinesFields::integrateFields() { - + surfaceScalarField phiSt(linearInterpolate(UsField_) & particleCloud_.mesh().Sf()); surfaceScalarField phiDyn(linearInterpolate(uDyn_) & particleCloud_.mesh().Sf()); - + fvScalarMatrix alphaStEqn ( fvm::ddt(alphaSt_) - + fvm::div(phiSt,alphaSt_) - == - Sds_ + + fvm::div(phiSt,alphaSt_) + == + Sds_ ); fvScalarMatrix alphaDynEqn ( fvm::ddt(alphaDyn_) - + fvm::div(phiDyn,alphaDyn_) - - fvm::laplacian(diffCoeff_,alphaDyn_) - == - -Sds_ + + fvm::div(phiDyn,alphaDyn_) + - fvm::laplacian(diffCoeff_,alphaDyn_) + == + -Sds_ ); alphaStEqn.solve(); alphaDynEqn.solve(); - + if(smoothing_) particleCloud_.smoothingM().smoothen(alphaDyn_); - + // limit hold-ups, should be done more elegantly - + scalar alphaStErr(0.0); scalar alphaDynErr1(0.0); scalar alphaDynErr2(0.0); forAll(alphaSt_, cellI) { if (alphaSt_[cellI] < 0.0) - { - alphaStErr += alphaSt_[cellI] * particleCloud_.mesh().V()[cellI]; - alphaSt_[cellI] = 0.0; - } - - if (alphaDyn_[cellI] < 0.0) - { - alphaDynErr1 += alphaDyn_[cellI] * particleCloud_.mesh().V()[cellI]; - alphaDyn_[cellI] = 0.0; - } - else if (alphaDyn_[cellI] > alphaDynMax_) - { - alphaDynErr2 += (alphaDyn_[cellI] - alphaDynMax_) * particleCloud_.mesh().V()[cellI]; - alphaDyn_[cellI] = alphaDynMax_; - } + { + alphaStErr += alphaSt_[cellI] * particleCloud_.mesh().V()[cellI]; + alphaSt_[cellI] = 0.0; + } + + if (alphaDyn_[cellI] < 0.0) + { + alphaDynErr1 += alphaDyn_[cellI] * particleCloud_.mesh().V()[cellI]; + alphaDyn_[cellI] = 0.0; + } + else if (alphaDyn_[cellI] > alphaDynMax_) + { + alphaDynErr2 += (alphaDyn_[cellI] - alphaDynMax_) * particleCloud_.mesh().V()[cellI]; + alphaDyn_[cellI] = alphaDynMax_; + } } - + if (verbose_) { Sout << "[" << Pstream::myProcNo() << "] " << "amount of alphaSt added because of positivity requirement: " << -alphaStErr << endl; Sout << "[" << Pstream::myProcNo() << "] " << "amount of alphaDyn added because of positivity requirement: " << -alphaDynErr1 << endl; Sout << "[" << Pstream::myProcNo() << "] " << "amount of alphaDyn removed because of max. value: " << -alphaDynErr2 << endl; } - + alphaSt_.correctBoundaryConditions(); alphaDyn_.correctBoundaryConditions(); - + massFluxDyn_ = rhoFine_ * fvc::interpolate(alphaDyn_) * phiDyn; } -void FinesFields::updateAlphaG() +void FinesFields::updateAlphaG() { - alphaG_ = max(voidfraction_ - alphaSt_ - alphaDyn_, critVoidfraction_); + alphaG_ = max(voidfraction_ - alphaSt_ - alphaDyn_, critVoidfraction_); } -void FinesFields::updateAlphaP() +void FinesFields::updateAlphaP() { alphaP_ = 1.0 - voidfraction_ + SMALL; } -void FinesFields::updateDHydMix() +void FinesFields::updateDHydMix() { forAll(dHydMix_,cellI) { scalar aPSt = alphaP_[cellI] + alphaSt_[cellI]; - if(aPSt < SMALL || aPSt > 1 - SMALL) - dHydMix_[cellI] = SMALL; - else - dHydMix_[cellI] = 2*(1 - aPSt) / (3*aPSt ) * dSauterMix_[cellI]; + if(aPSt < SMALL || aPSt > 1 - SMALL) + dHydMix_[cellI] = SMALL; + else + dHydMix_[cellI] = 2*(1 - aPSt) / (3*aPSt ) * dSauterMix_[cellI]; } dHydMix_.correctBoundaryConditions(); } @@ -453,54 +454,54 @@ void FinesFields::updateDragCoeff() volScalarField Ref = dFine_ * alphaG_ / nuAve_ * mag(U_ - uDyn_); scalar Cd(0.0); scalar Ref1(0.0); - + // calculate drag coefficient for cells forAll(DragCoeff_,cellI) { Ref1 = Ref[cellI]; if(Ref1 <= SMALL) - Cd = 24.0 / SMALL; + Cd = 24.0 / SMALL; else if(Ref1 <= 1.0) - Cd = 24.0 / Ref1; - else if(Ref1 <= 1000) - Cd = 24 * (1.0 + 0.15 * Foam::pow(Ref1,0.687) ) / Ref1; - else - Cd = 0.44; - DragCoeff_[cellI] = Cd * beta[cellI]; + Cd = 24.0 / Ref1; + else if(Ref1 <= 1000) + Cd = 24 * (1.0 + 0.15 * Foam::pow(Ref1,0.687) ) / Ref1; + else + Cd = 0.44; + DragCoeff_[cellI] = Cd * beta[cellI]; } - + // calculate drag coefficient for faces forAll(DragCoeff_.boundaryField(), patchI) forAll(DragCoeff_.boundaryField()[patchI], faceI) { Ref1 = Ref.boundaryField()[patchI][faceI]; if(Ref1 <= SMALL) - Cd = 24.0 / SMALL; + Cd = 24.0 / SMALL; else if(Ref1 <= 1.0) - Cd = 24.0 / Ref1; - else if(Ref1 <= 1000) - Cd = 24 * (1.0 + 0.15 * Foam::pow(Ref1,0.687) ) / Ref1; - else - Cd = 0.44; - DragCoeff_.boundaryFieldRef()[patchI][faceI] = Cd * beta.boundaryFieldRef()[patchI][faceI]; + Cd = 24.0 / Ref1; + else if(Ref1 <= 1000) + Cd = 24 * (1.0 + 0.15 * Foam::pow(Ref1,0.687) ) / Ref1; + else + Cd = 0.44; + DragCoeff_.boundaryFieldRef()[patchI][faceI] = Cd * beta.boundaryFieldRef()[patchI][faceI]; } - + DragCoeff_ = max( DragCoeff_, dimensionedScalar("SMALL", dimensionSet(1,-3,-1,0,0), SMALL) ); } -void FinesFields::updateDSauter() +void FinesFields::updateDSauter() { forAll(dSauterMix_,cellI) { scalar aP = alphaP_[cellI]; - scalar aSt = alphaSt_[cellI]; - if(aSt < SMALL) - dSauterMix_[cellI] = dSauter_[cellI]; - else if(aP < SMALL) - dSauterMix_[cellI] = dFine_.value(); - else - dSauterMix_[cellI] = (aP + aSt) / (aP / dSauter_[cellI] + aSt / dFine_.value() ); + scalar aSt = alphaSt_[cellI]; + if(aSt < SMALL) + dSauterMix_[cellI] = dSauter_[cellI]; + else if(aP < SMALL) + dSauterMix_[cellI] = dFine_.value(); + else + dSauterMix_[cellI] = (aP + aSt) / (aP / dSauter_[cellI] + aSt / dFine_.value() ); } dSauterMix_.correctBoundaryConditions(); } @@ -514,7 +515,7 @@ void FinesFields::updateFanningCoeff() } -void FinesFields::updateFroude() +void FinesFields::updateFroude() { // seems like different authors use different conventions for the Froude number // Chen et al. (1994) define it in terms of a superficial velocity, @@ -534,27 +535,27 @@ void FinesFields::updateUDyn() volScalarField denom = FanningCoeff_ + DragCoeff_; uDyn_ = num / denom; - + // limit uDyn for stability reasons forAll(uDyn_,cellI) { scalar mU(mag(U_[cellI])); - scalar muDyn(mag(uDyn_[cellI])); + scalar muDyn(mag(uDyn_[cellI])); if(muDyn > mU && muDyn > SMALL) - { - uDyn_[cellI] *= mU / muDyn; - } + { + uDyn_[cellI] *= mU / muDyn; + } } - + forAll(uDyn_.boundaryField(), patchI) forAll(uDyn_.boundaryField()[patchI], faceI) { - scalar mU(mag(U_.boundaryField()[patchI][faceI])); - scalar muDyn(mag(uDyn_.boundaryField()[patchI][faceI])); + scalar mU(mag(U_.boundaryField()[patchI][faceI])); + scalar muDyn(mag(uDyn_.boundaryField()[patchI][faceI])); if(muDyn > mU && muDyn > SMALL) - { - uDyn_.boundaryFieldRef()[patchI][faceI] *= mU / muDyn; - } + { + uDyn_.boundaryFieldRef()[patchI][faceI] *= mU / muDyn; + } } } diff --git a/src/lagrangian/cfdemParticle/subModels/forceModel/Fines/FinesFields.H b/src/lagrangian/cfdemParticle/subModels/forceModel/Fines/FinesFields.H index 75ae75b9..3098ff2f 100644 --- a/src/lagrangian/cfdemParticle/subModels/forceModel/Fines/FinesFields.H +++ b/src/lagrangian/cfdemParticle/subModels/forceModel/Fines/FinesFields.H @@ -39,11 +39,11 @@ class FinesFields { private: cfdemCloud& particleCloud_; - + dictionary propsDict_; - + bool smoothing_; - + bool verbose_; word velFieldName_; @@ -54,95 +54,95 @@ private: volScalarField& voidfraction_; - word UsFieldName_; + word UsFieldName_; const volVectorField& UsField_; - + word pFieldName_; const volScalarField& p_; - + word rhoGFieldName_; const volScalarField& rhoG_; - + const volScalarField& dSauter_; - + volScalarField alphaG_; - + volScalarField alphaDyn_; - + volScalarField alphaP_; volScalarField alphaSt_; volScalarField deltaAlpha_; - + volScalarField dHydMix_; - + volScalarField DragCoeff_; - + volScalarField dSauterMix_; - + volScalarField FanningCoeff_; - + volScalarField Froude_; - + volScalarField Sds_; - + //volVectorField massFluxDyn_; surfaceScalarField massFluxDyn_; - + volVectorField uDyn_; dimensionedScalar dFine_; - + dimensionedScalar diffCoeff_; - + dimensionedScalar nuAve_; - + dimensionedScalar rhoFine_; - + const dimensionedVector g_; - + scalar alphaDynMax_; - + scalar alphaMax_; - + scalar critVoidfraction_; - + scalar depRate_; - + scalar exponent_; - + scalar nCrit_; - + scalar poresizeWidth_; - + scalar prefactor_; - + scalar ratioHydraulicPore_; - + void calcSource(); - + void integrateFields(); - + void updateAlphaG(); - + void updateAlphaP(); - + void updateDHydMix(); - + void updateDragCoeff(); - + void updateDSauter(); - + void updateFanningCoeff(); - + void updateFroude(); - + void updateUDyn(); - + public: //- Runtime type information @@ -165,7 +165,7 @@ public: // Member Functions - + void update(); }; diff --git a/src/lagrangian/cfdemParticle/subModels/forceModel/GidaspowDrag/GidaspowDrag.C b/src/lagrangian/cfdemParticle/subModels/forceModel/GidaspowDrag/GidaspowDrag.C index eaaeb94a..57bb5d3c 100644 --- a/src/lagrangian/cfdemParticle/subModels/forceModel/GidaspowDrag/GidaspowDrag.C +++ b/src/lagrangian/cfdemParticle/subModels/forceModel/GidaspowDrag/GidaspowDrag.C @@ -75,8 +75,8 @@ GidaspowDrag::GidaspowDrag switchingVoidfraction_(0.8) { //Append the field names to be probed - particleCloud_.probeM().initialize(typeName, "gidaspowDrag.logDat"); - particleCloud_.probeM().vectorFields_.append("dragForce"); //first entry must be the force + particleCloud_.probeM().initialize(typeName, typeName+".logDat"); + particleCloud_.probeM().vectorFields_.append("dragForce"); // first entry must be the force particleCloud_.probeM().vectorFields_.append("Urel"); particleCloud_.probeM().scalarFields_.append("Rep"); particleCloud_.probeM().scalarFields_.append("betaP"); @@ -86,18 +86,18 @@ GidaspowDrag::GidaspowDrag // init force sub model setForceSubModels(propsDict_); // define switches which can be read from dict - forceSubM(0).setSwitchesList(0,true); // activate treatExplicit switch - forceSubM(0).setSwitchesList(2,true); // activate implDEM switch - forceSubM(0).setSwitchesList(3,true); // activate search for verbose switch - forceSubM(0).setSwitchesList(4,true); // activate search for interpolate switch - forceSubM(0).setSwitchesList(8,true); // activate scalarViscosity switch + forceSubM(0).setSwitchesList(SW_TREAT_FORCE_EXPLICIT,true); // activate treatExplicit switch + forceSubM(0).setSwitchesList(SW_IMPL_FORCE_DEM,true); // activate implDEM switch + forceSubM(0).setSwitchesList(SW_VERBOSE,true); // activate search for verbose switch + forceSubM(0).setSwitchesList(SW_INTERPOLATION,true); // activate search for interpolate switch + forceSubM(0).setSwitchesList(SW_SCALAR_VISCOSITY,true); // activate scalarViscosity switch forceSubM(0).readSwitches(); particleCloud_.checkCG(true); if (propsDict_.found("scale")) - scaleDia_=scalar(readScalar(propsDict_.lookup("scale"))); + scaleDia_ = scalar(readScalar(propsDict_.lookup("scale"))); if (propsDict_.found("scaleDrag")) - scaleDrag_=scalar(readScalar(propsDict_.lookup("scaleDrag"))); + scaleDrag_ = scalar(readScalar(propsDict_.lookup("scaleDrag"))); if (propsDict_.found("switchingVoidfraction")) switchingVoidfraction_ = readScalar(propsDict_.lookup("switchingVoidfraction")); @@ -114,10 +114,13 @@ GidaspowDrag::~GidaspowDrag() void GidaspowDrag::setForce() const { - if (scaleDia_ > 1) + if (scaleDia_ > 1.) + { Info << "Gidaspow using scale = " << scaleDia_ << endl; - else if (particleCloud_.cg() > 1){ - scaleDia_=particleCloud_.cg(); + } + else if (particleCloud_.cg() > 1.) + { + scaleDia_ = particleCloud_.cg(); Info << "Gidaspow using scale from liggghts cg = " << scaleDia_ << endl; } @@ -128,7 +131,7 @@ void GidaspowDrag::setForce() const scalar voidfraction(1); vector Ufluid(0,0,0); vector drag(0,0,0); - label cellI=0; + label cellI = 0; vector Us(0,0,0); vector Ur(0,0,0); @@ -145,13 +148,13 @@ void GidaspowDrag::setForce() const vector dragExplicit(0,0,0); scalar dragCoefficient(0); - + interpolationCellPoint voidfractionInterpolator_(voidfraction_); interpolationCellPoint UInterpolator_(U_); #include "setupProbeModel.H" - for(int index = 0;index < particleCloud_.numberOfParticles(); ++index) + for(int index = 0; index < particleCloud_.numberOfParticles(); ++index) { //if(mask[index][0]) //{ @@ -160,26 +163,26 @@ void GidaspowDrag::setForce() const dragExplicit = vector(0,0,0); betaP = 0; Vs = 0; - Ufluid =vector(0,0,0); - voidfraction=0; + Ufluid = vector(0,0,0); + voidfraction = 0; dragCoefficient = 0; if (cellI > -1) // particle Found { - if( forceSubM(0).interpolation() ) + if ( forceSubM(0).interpolation() ) { - position = particleCloud_.position(index); + 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.0) voidfraction = 1.0; + else if (voidfraction < 0.1) voidfraction = 0.10; } else { - voidfraction = voidfraction_[cellI]; + voidfraction = voidfraction_[cellI]; Ufluid = U_[cellI]; } @@ -190,14 +193,14 @@ void GidaspowDrag::setForce() const rho = rhoField[cellI]; nuf = nufField[cellI]; - Rep=0.0; + Rep = 0.0; localPhiP = 1.0f-voidfraction+SMALL; Vs = ds*ds*ds*M_PI/6; // calc particle's drag coefficient (i.e., Force per unit slip velocity and per m³ PARTICLE) - if(voidfraction > switchingVoidfraction_) //dilute + if (voidfraction > switchingVoidfraction_) //dilute { - Rep=ds/scaleDia_*voidfraction*magUr/nuf; + Rep = ds/scaleDia_*voidfraction*magUr/nuf; CdMagUrLag = (24.0*nuf/(ds/scaleDia_*voidfraction)) //1/magUr missing here, but compensated in expression for betaP! *(scalar(1.0)+0.15*Foam::pow(Rep, 0.687)); @@ -218,7 +221,7 @@ void GidaspowDrag::setForce() const // calc particle's drag dragCoefficient = Vs*betaP*scaleDrag_; - if (modelType_=="B") + if (modelType_ == "B") dragCoefficient /= voidfraction; drag = dragCoefficient * Ur; @@ -226,7 +229,7 @@ void GidaspowDrag::setForce() const // explicitCorr forceSubM(0).explicitCorr(drag,dragExplicit,dragCoefficient,Ufluid,U_[cellI],Us,UsField_[cellI],forceSubM(0).verbose()); - if(forceSubM(0).verbose() && index >=0 && index <2) + if (forceSubM(0).verbose() && index >= 0 && index < 2) { Pout << "cellI = " << cellI << endl; Pout << "index = " << index << endl; @@ -244,7 +247,7 @@ void GidaspowDrag::setForce() const } //Set value fields and write the probe - if(probeIt_) + if (probeIt_) { #include "setupProbeModelfields.H" vValues.append(drag); //first entry must the be the force diff --git a/src/lagrangian/cfdemParticle/subModels/forceModel/GidaspowDrag/GidaspowDrag.H b/src/lagrangian/cfdemParticle/subModels/forceModel/GidaspowDrag/GidaspowDrag.H index 61c7d191..76102bcb 100644 --- a/src/lagrangian/cfdemParticle/subModels/forceModel/GidaspowDrag/GidaspowDrag.H +++ b/src/lagrangian/cfdemParticle/subModels/forceModel/GidaspowDrag/GidaspowDrag.H @@ -74,7 +74,7 @@ private: const scalar phi_; - word UsFieldName_; + word UsFieldName_; const volVectorField& UsField_; // the average particle velocity field diff --git a/src/lagrangian/cfdemParticle/subModels/forceModel/KochHillDrag/KochHillDrag.C b/src/lagrangian/cfdemParticle/subModels/forceModel/KochHillDrag/KochHillDrag.C index 8a7ec04f..7368d942 100644 --- a/src/lagrangian/cfdemParticle/subModels/forceModel/KochHillDrag/KochHillDrag.C +++ b/src/lagrangian/cfdemParticle/subModels/forceModel/KochHillDrag/KochHillDrag.C @@ -74,7 +74,7 @@ KochHillDrag::KochHillDrag scaleDrag_(1.) { //Append the field names to be probed - particleCloud_.probeM().initialize(typeName, "kochHillDrag.logDat"); + particleCloud_.probeM().initialize(typeName, typeName+".logDat"); particleCloud_.probeM().vectorFields_.append("dragForce"); //first entry must the be the force particleCloud_.probeM().vectorFields_.append("Urel"); //other are debug particleCloud_.probeM().scalarFields_.append("Rep"); //other are debug @@ -86,12 +86,12 @@ KochHillDrag::KochHillDrag setForceSubModels(propsDict_); // define switches which can be read from dict - forceSubM(0).setSwitchesList(0,true); // activate treatExplicit switch - forceSubM(0).setSwitchesList(2,true); // activate implDEM switch - forceSubM(0).setSwitchesList(3,true); // activate search for verbose switch - forceSubM(0).setSwitchesList(4,true); // activate search for interpolate switch - forceSubM(0).setSwitchesList(7,true); // activate implForceDEMacc switch - forceSubM(0).setSwitchesList(8,true); // activate scalarViscosity switch + forceSubM(0).setSwitchesList(SW_TREAT_FORCE_EXPLICIT,true); // activate treatExplicit switch + forceSubM(0).setSwitchesList(SW_IMPL_FORCE_DEM,true); // activate implDEM switch + forceSubM(0).setSwitchesList(SW_VERBOSE,true); // activate search for verbose switch + forceSubM(0).setSwitchesList(SW_INTERPOLATION,true); // activate search for interpolate switch + forceSubM(0).setSwitchesList(SW_IMPL_FORCE_DEM_ACCUMULATED,true); // activate implForceDEMacc switch + forceSubM(0).setSwitchesList(SW_SCALAR_VISCOSITY,true); // activate scalarViscosity switch // read those switches defined above, if provided in dict forceSubM(0).readSwitches(); @@ -99,9 +99,9 @@ KochHillDrag::KochHillDrag particleCloud_.checkCG(true); if (propsDict_.found("scale")) - scaleDia_=scalar(readScalar(propsDict_.lookup("scale"))); + scaleDia_ = scalar(readScalar(propsDict_.lookup("scale"))); if (propsDict_.found("scaleDrag")) - scaleDrag_=scalar(readScalar(propsDict_.lookup("scaleDrag"))); + scaleDrag_ = scalar(readScalar(propsDict_.lookup("scaleDrag"))); } @@ -134,7 +134,7 @@ void KochHillDrag::setForce() const vector drag(0,0,0); vector dragExplicit(0,0,0); scalar dragCoefficient(0); - label cellI=0; + label cellI = 0; vector Us(0,0,0); vector Ur(0,0,0); @@ -154,16 +154,16 @@ void KochHillDrag::setForce() const #include "setupProbeModel.H" - for (int index=0; index -1) // particle Found { @@ -174,8 +174,8 @@ void KochHillDrag::setForce() const Ufluid = UInterpolator_.interpolate(position,cellI); //Ensure interpolated void fraction to be meaningful // Info << " --> voidfraction: " << voidfraction << endl; - if (voidfraction > 1.00) voidfraction = 1.00; - if (voidfraction < 0.40) voidfraction = 0.40; + if (voidfraction > 1.0) voidfraction = 1.0; + else if (voidfraction < 0.4) voidfraction = 0.4; } else { @@ -227,10 +227,10 @@ void KochHillDrag::setForce() const // calc particle's drag dragCoefficient = Vs*betaP*scaleDrag_; - if (modelType_=="B") + if (modelType_ == "B") dragCoefficient /= voidfraction; - if (forceSubM(0).switches()[7]) // implForceDEMaccumulated=true + if (forceSubM(0).switches()[SW_IMPL_FORCE_DEM_ACCUMULATED]) // implForceDEMaccumulated=true { //get drag from the particle itself for (int j=0; j<3; j++) drag[j] = particleCloud_.fAccs()[index][j]/couplingInterval; @@ -270,7 +270,7 @@ void KochHillDrag::setForce() const sValues.append(betaP); sValues.append(voidfraction); particleCloud_.probeM().writeProbe(index, sValues, vValues); - } + } } // write particle based data to global array diff --git a/src/lagrangian/cfdemParticle/subModels/forceModel/KochHillDrag/KochHillDrag.H b/src/lagrangian/cfdemParticle/subModels/forceModel/KochHillDrag/KochHillDrag.H index 95d6712e..42c921f3 100644 --- a/src/lagrangian/cfdemParticle/subModels/forceModel/KochHillDrag/KochHillDrag.H +++ b/src/lagrangian/cfdemParticle/subModels/forceModel/KochHillDrag/KochHillDrag.H @@ -29,7 +29,7 @@ Description and OpenFOAM(R). Note: this code is not part of OpenFOAM(R) (see DISCLAIMER). Koch, Hill drag law - based on Koch Hill 2001,"Inertial effects in suspensions and porous-media + based on Koch Hill 2001,"Inertial effects in suspensions and porous-media flows", Annual Review of fluid mechanics. including interpolation of the velocity to the exact position including drag coefficient for implicit drag for DEM @@ -72,9 +72,9 @@ private: const volScalarField& voidfraction_; - word UsFieldName_; + word UsFieldName_; - const volVectorField& UsField_; + const volVectorField& UsField_; mutable scalar scaleDia_; diff --git a/src/lagrangian/cfdemParticle/subModels/forceModel/KochHillRWDrag/KochHillRWDrag.C b/src/lagrangian/cfdemParticle/subModels/forceModel/KochHillRWDrag/KochHillRWDrag.C index 083ae54d..616069c8 100755 --- a/src/lagrangian/cfdemParticle/subModels/forceModel/KochHillRWDrag/KochHillRWDrag.C +++ b/src/lagrangian/cfdemParticle/subModels/forceModel/KochHillRWDrag/KochHillRWDrag.C @@ -79,16 +79,16 @@ KochHillRWDrag::KochHillRWDrag RanGen_(label(0)) { - if (propsDict_.found("verbose")) verbose_=true; - if (propsDict_.found("interpolation")) interpolation_=true; - if (propsDict_.found("randomTauE")) randomTauE_=true; + if (propsDict_.found("verbose")) verbose_ = true; + if (propsDict_.found("interpolation")) interpolation_ = true; + if (propsDict_.found("randomTauE")) randomTauE_ = true; // init force sub model setForceSubModels(propsDict_); // define switches which can be read from dict - forceSubM(0).setSwitchesList(0,true); // activate treatExplicit switch - forceSubM(0).setSwitchesList(2,true); // activate implDEM switch + forceSubM(0).setSwitchesList(SW_TREAT_FORCE_EXPLICIT,true); // activate treatExplicit switch + forceSubM(0).setSwitchesList(SW_IMPL_FORCE_DEM,true); // activate implDEM switch // read those switches defined above, if provided in dict forceSubM(0).readSwitches(); @@ -120,8 +120,8 @@ KochHillRWDrag::KochHillRWDrag KochHillRWDrag::~KochHillRWDrag() { - delete partTime_; - delete partUfluct_; + particleCloud_.dataExchangeM().destroy(partTime_, 1); + particleCloud_.dataExchangeM().destroy(partUfluct_, 3); } // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // @@ -152,7 +152,7 @@ void KochHillRWDrag::setForce() const vector drag(0,0,0); vector dragExplicit(0,0,0); scalar dragCoefficient(0); - label cellI=0; + label cellI = 0; vector Us(0,0,0); vector Ur(0,0,0); @@ -195,17 +195,17 @@ void KochHillRWDrag::setForce() const //Info << "RW-TEST: We are in setForce() at t = " << t << endl; // TEST-Output - for (int index=0; index voidfraction: " << voidfraction << endl; - if (voidfraction > 1.00) voidfraction = 1.00; - if (voidfraction < 0.40) voidfraction = 0.40; + if (voidfraction > 1.0) voidfraction = 1.0; + else if (voidfraction < 0.4) voidfraction = 0.4; } else { diff --git a/src/lagrangian/cfdemParticle/subModels/forceModel/KochHillRWDrag/KochHillRWDrag.H b/src/lagrangian/cfdemParticle/subModels/forceModel/KochHillRWDrag/KochHillRWDrag.H index ba80cfb1..16d3f463 100755 --- a/src/lagrangian/cfdemParticle/subModels/forceModel/KochHillRWDrag/KochHillRWDrag.H +++ b/src/lagrangian/cfdemParticle/subModels/forceModel/KochHillRWDrag/KochHillRWDrag.H @@ -29,7 +29,7 @@ Description and OpenFOAM(R). Note: this code is not part of OpenFOAM(R) (see DISCLAIMER). Koch, Hill drag law - based on Koch Hill 2001,"Inertial effects in suspensions and porous-media + based on Koch Hill 2001,"Inertial effects in suspensions and porous-media flows", Annual Review of fluid mechanics. including interpolation of the velocity to the exact position including drag coefficient for implicit drag for DEM @@ -75,7 +75,7 @@ private: const volScalarField& voidfraction_; - word UsFieldName_; + word UsFieldName_; const volVectorField& UsField_; // the average particle velocity field diff --git a/src/lagrangian/cfdemParticle/subModels/forceModel/LaEuScalarTemp/LaEuScalarTemp.C b/src/lagrangian/cfdemParticle/subModels/forceModel/LaEuScalarTemp/LaEuScalarTemp.C index 09cb8c46..bb8813fc 100644 --- a/src/lagrangian/cfdemParticle/subModels/forceModel/LaEuScalarTemp/LaEuScalarTemp.C +++ b/src/lagrangian/cfdemParticle/subModels/forceModel/LaEuScalarTemp/LaEuScalarTemp.C @@ -89,9 +89,9 @@ LaEuScalarTemp::LaEuScalarTemp setForceSubModels(propsDict_); // define switches which can be read from dict - forceSubM(0).setSwitchesList(3,true); // activate search for verbose switch - forceSubM(0).setSwitchesList(4,true); // activate search for interpolate switch - forceSubM(0).setSwitchesList(8,true); // activate scalarViscosity switch + forceSubM(0).setSwitchesList(SW_VERBOSE,true); // activate search for verbose switch + forceSubM(0).setSwitchesList(SW_INTERPOLATION,true); // activate search for interpolate switch + forceSubM(0).setSwitchesList(SW_SCALAR_VISCOSITY,true); // activate scalarViscosity switch // read those switches defined above, if provided in dict forceSubM(0).readSwitches(); @@ -105,15 +105,15 @@ LaEuScalarTemp::LaEuScalarTemp LaEuScalarTemp::~LaEuScalarTemp() { - delete partTemp_; - delete partHeatFlux_; + particleCloud_.dataExchangeM().destroy(partTemp_,1); + particleCloud_.dataExchangeM().destroy(partHeatFlux_,1); } // * * * * * * * * * * * * * * * private Member Functions * * * * * * * * * * * * * // void LaEuScalarTemp::allocateMyArrays() const { // get memory for 2d arrays - double initVal=0.0; + double initVal = 0.0; particleCloud_.dataExchangeM().allocateArray(partTemp_,initVal,1); // field/initVal/with/lenghtFromLigghts particleCloud_.dataExchangeM().allocateArray(partHeatFlux_,initVal,1); } @@ -156,14 +156,14 @@ void LaEuScalarTemp::manipulateScalarField(volScalarField& EuField) const interpolationCellPoint UInterpolator_(U_); interpolationCellPoint TInterpolator_(tempField_); - for(int index = 0;index < particleCloud_.numberOfParticles(); ++index) + for(int index = 0; index < particleCloud_.numberOfParticles(); ++index) { //if(particleCloud_.regionM().inRegion()[index][0]) //{ cellI = particleCloud_.cellIDs()[index][0]; - if(cellI >= 0) + if (cellI >= 0) { - if(forceSubM(0).interpolation()) + if (forceSubM(0).interpolation()) { vector position = particleCloud_.position(index); voidfraction = voidfractionInterpolator_.interpolate(position,cellI); @@ -206,7 +206,7 @@ void LaEuScalarTemp::manipulateScalarField(volScalarField& EuField) const partHeatFlux_[index][0] = partHeatFlux; - if(forceSubM(0).verbose() && index >=0 && index <2) + if(forceSubM(0).verbose() && index >= 0 && index < 2) { Info << "partHeatFlux = " << partHeatFlux << endl; Info << "magUr = " << magUr << endl; diff --git a/src/lagrangian/cfdemParticle/subModels/forceModel/LaEuScalarTemp/LaEuScalarTemp.H b/src/lagrangian/cfdemParticle/subModels/forceModel/LaEuScalarTemp/LaEuScalarTemp.H index 01e6354f..c9de04f9 100644 --- a/src/lagrangian/cfdemParticle/subModels/forceModel/LaEuScalarTemp/LaEuScalarTemp.H +++ b/src/lagrangian/cfdemParticle/subModels/forceModel/LaEuScalarTemp/LaEuScalarTemp.H @@ -71,7 +71,7 @@ private: word voidfractionFieldName_; - const volScalarField& voidfraction_; // ref to voidfraction field + const volScalarField& voidfraction_; // ref to voidfraction field scalar maxSource_; // max (limited) value of src field @@ -81,7 +81,7 @@ private: word partTempName_; - mutable double **partTemp_; // Lagrangian array + mutable double **partTemp_; // Lagrangian array word partHeatFluxName_; diff --git a/src/lagrangian/cfdemParticle/subModels/forceModel/MeiLift/MeiLift.C b/src/lagrangian/cfdemParticle/subModels/forceModel/MeiLift/MeiLift.C index d8f688d6..0dca7e3b 100644 --- a/src/lagrangian/cfdemParticle/subModels/forceModel/MeiLift/MeiLift.C +++ b/src/lagrangian/cfdemParticle/subModels/forceModel/MeiLift/MeiLift.C @@ -73,16 +73,16 @@ MeiLift::MeiLift // init force sub model setForceSubModels(propsDict_); // define switches which can be read from dict - forceSubM(0).setSwitchesList(0,true); // activate treatExplicit switch - forceSubM(0).setSwitchesList(3,true); // activate search for verbose switch - forceSubM(0).setSwitchesList(4,true); // activate search for interpolate switch - forceSubM(0).setSwitchesList(8,true); // activate scalarViscosity switch + forceSubM(0).setSwitchesList(SW_TREAT_FORCE_EXPLICIT,true); // activate treatExplicit switch + forceSubM(0).setSwitchesList(SW_VERBOSE,true); // activate search for verbose switch + forceSubM(0).setSwitchesList(SW_INTERPOLATION,true); // activate search for interpolate switch + forceSubM(0).setSwitchesList(SW_SCALAR_VISCOSITY,true); // activate scalarViscosity switch forceSubM(0).readSwitches(); particleCloud_.checkCG(false); //Append the field names to be probed - particleCloud_.probeM().initialize(typeName, "meiLift.logDat"); + particleCloud_.probeM().initialize(typeName, typeName+".logDat"); particleCloud_.probeM().vectorFields_.append("liftForce"); //first entry must the be the force particleCloud_.probeM().vectorFields_.append("Urel"); //other are debug particleCloud_.probeM().vectorFields_.append("vorticity"); //other are debug @@ -134,11 +134,11 @@ void MeiLift::setForce() const #include "setupProbeModel.H" - for(int index = 0;index < particleCloud_.numberOfParticles(); index++) + for(int index = 0; index < particleCloud_.numberOfParticles(); ++index) { //if(mask[index][0]) //{ - lift = vector::zero; + lift = vector::zero; label cellI = particleCloud_.cellIDs()[index][0]; if (cellI > -1) // particle Found @@ -147,10 +147,10 @@ void MeiLift::setForce() const if( forceSubM(0).interpolation() ) { - position = particleCloud_.position(index); - Ur = UInterpolator_.interpolate(position,cellI) - - Us; - vorticity = VorticityInterpolator_.interpolate(position,cellI); + position = particleCloud_.position(index); + Ur = UInterpolator_.interpolate(position,cellI) + - Us; + vorticity = VorticityInterpolator_.interpolate(position,cellI); } else { @@ -159,7 +159,7 @@ void MeiLift::setForce() const vorticity=vorticityField[cellI]; } - magUr = mag(Ur); + magUr = mag(Ur); magVorticity = mag(vorticity); if (magUr > 0 && magVorticity > 0) @@ -170,26 +170,26 @@ void MeiLift::setForce() const // calc dimensionless properties Rep = ds*magUr/nuf; - Rew = magVorticity*ds*ds/nuf; + Rew = magVorticity*ds*ds/nuf; - alphaStar = magVorticity*ds/magUr/2.0; - epsilon = sqrt(2.0*alphaStar /Rep ); - omega_star=2.0*alphaStar; + alphaStar = magVorticity*ds/magUr/2.0; + epsilon = sqrt(2.0*alphaStar /Rep ); + omega_star = 2.0*alphaStar; //Basic model for the correction to the Saffman lift //Based on McLaughlin (1991) if(epsilon < 0.1) { - J_star = -140 *epsilon*epsilon*epsilon*epsilon*epsilon + J_star = -140 *epsilon*epsilon*epsilon*epsilon*epsilon *log( 1./(epsilon*epsilon+SMALL) ); } else if(epsilon > 20) { - J_star = 1.0-0.287/(epsilon*epsilon+SMALL); + J_star = 1.0-0.287/(epsilon*epsilon+SMALL); } else { - J_star = 0.3 + J_star = 0.3 *( 1.0 +tanh( 2.5 * log10(epsilon+0.191) ) ) @@ -197,11 +197,11 @@ void MeiLift::setForce() const +tanh( 6.0 * (epsilon-0.32) ) ); } - Cl=J_star*4.11*epsilon; //multiply McLaughlin's correction to the basic Saffman model + Cl = J_star * 4.11 * epsilon; //multiply McLaughlin's correction to the basic Saffman model - //Second order terms given by Loth and Dorgan 2009 + //Second order terms given by Loth and Dorgan 2009 if(useSecondOrderTerms_) - { + { Omega_eq = omega_star/2.0*(1.0-0.0075*Rew)*(1.0-0.062*sqrt(Rep)-0.001*Rep); Cl_star=1.0-(0.675+0.15*(1.0+tanh(0.28*(omega_star/2.0-2.0))))*tanh(0.18*sqrt(Rep)); Cl += Omega_eq*Cl_star; @@ -209,21 +209,21 @@ void MeiLift::setForce() const lift = 0.125*M_PI *rho - *Cl + *Cl *magUr*Ur^vorticity/magVorticity *ds*ds; - if (modelType_=="B") + if (modelType_ == "B") { voidfraction = particleCloud_.voidfraction(index); lift /= voidfraction; } } - //********************************** + //********************************** //SAMPLING AND VERBOSE OUTOUT - if( forceSubM(0).verbose() ) - { + if ( forceSubM(0).verbose() ) + { Pout << "index = " << index << endl; Pout << "Us = " << Us << endl; Pout << "Ur = " << Ur << endl; @@ -240,7 +240,7 @@ void MeiLift::setForce() const } //Set value fields and write the probe - if(probeIt_) + if (probeIt_) { #include "setupProbeModelfields.H" vValues.append(lift); //first entry must the be the force @@ -252,7 +252,7 @@ void MeiLift::setForce() const particleCloud_.probeM().writeProbe(index, sValues, vValues); } // END OF SAMPLING AND VERBOSE OUTOUT - //********************************** + //********************************** } // write particle based data to global array diff --git a/src/lagrangian/cfdemParticle/subModels/forceModel/SchillerNaumannDrag/SchillerNaumannDrag.C b/src/lagrangian/cfdemParticle/subModels/forceModel/SchillerNaumannDrag/SchillerNaumannDrag.C index 3fc632cd..0c89a4aa 100644 --- a/src/lagrangian/cfdemParticle/subModels/forceModel/SchillerNaumannDrag/SchillerNaumannDrag.C +++ b/src/lagrangian/cfdemParticle/subModels/forceModel/SchillerNaumannDrag/SchillerNaumannDrag.C @@ -69,7 +69,7 @@ SchillerNaumannDrag::SchillerNaumannDrag U_(sm.mesh().lookupObject (velFieldName_)) { //Append the field names to be probed - particleCloud_.probeM().initialize(typeName, "schillerNaumannDrag.logDat"); + particleCloud_.probeM().initialize(typeName, typeName+".logDat"); particleCloud_.probeM().vectorFields_.append("dragForce"); //first entry must the be the force particleCloud_.probeM().vectorFields_.append("Urel"); //other are debug particleCloud_.probeM().scalarFields_.append("Rep"); //other are debug @@ -82,7 +82,7 @@ SchillerNaumannDrag::SchillerNaumannDrag setForceSubModels(propsDict_); // define switches which can be read from dict - forceSubM(0).setSwitchesList(0,true); // activate treatExplicit switch + forceSubM(0).setSwitchesList(SW_TREAT_FORCE_EXPLICIT,true); // activate treatExplicit switch // read those switches defined above, if provided in dict forceSubM(0).readSwitches(); diff --git a/src/lagrangian/cfdemParticle/subModels/forceModel/ShirgaonkarIB/ShirgaonkarIB.C b/src/lagrangian/cfdemParticle/subModels/forceModel/ShirgaonkarIB/ShirgaonkarIB.C index 59569a45..d0002a81 100644 --- a/src/lagrangian/cfdemParticle/subModels/forceModel/ShirgaonkarIB/ShirgaonkarIB.C +++ b/src/lagrangian/cfdemParticle/subModels/forceModel/ShirgaonkarIB/ShirgaonkarIB.C @@ -72,7 +72,7 @@ ShirgaonkarIB::ShirgaonkarIB p_(sm.mesh().lookupObject (pressureFieldName_)) { //Append the field names to be probed - particleCloud_.probeM().initialize(typeName, "shirgaonkarIB.logDat"); + particleCloud_.probeM().initialize(typeName, typeName+".logDat"); particleCloud_.probeM().vectorFields_.append("dragForce"); //first entry must the be the force particleCloud_.probeM().writeHeader(); @@ -89,7 +89,7 @@ ShirgaonkarIB::ShirgaonkarIB setForceSubModels(propsDict_); // define switches which can be read from dict - forceSubM(0).setSwitchesList(0,true); // activate treatExplicit switch + forceSubM(0).setSwitchesList(SW_TREAT_FORCE_EXPLICIT,true); // activate treatExplicit switch // read those switches defined above, if provided in dict forceSubM(0).readSwitches(); diff --git a/src/lagrangian/cfdemParticle/subModels/forceModel/dSauter/dSauter.C b/src/lagrangian/cfdemParticle/subModels/forceModel/dSauter/dSauter.C index 0ccdfc56..081742c3 100644 --- a/src/lagrangian/cfdemParticle/subModels/forceModel/dSauter/dSauter.C +++ b/src/lagrangian/cfdemParticle/subModels/forceModel/dSauter/dSauter.C @@ -92,20 +92,20 @@ dSauter::dSauter ), sm.mesh(), dimensionedScalar("zero", dimensionSet(0,1,0,0,0), 0), - "zeroGradient" + "zeroGradient" ) { allocateMyArrays(); dSauter_.write(); - + // init force sub model setForceSubModels(propsDict_); if (propsDict_.found("scaleCG")) - scaleDia_=scalar(readScalar(propsDict_.lookup("scaleCG"))); + scaleDia_ = scalar(readScalar(propsDict_.lookup("scaleCG"))); if (propsDict_.found("scaleDist")) - scaleDiaDist_=scalar(readScalar(propsDict_.lookup("scaleDist"))); + scaleDiaDist_ = scalar(readScalar(propsDict_.lookup("scaleDist"))); } @@ -113,27 +113,31 @@ dSauter::dSauter dSauter::~dSauter() { - delete d2_; - delete d3_; + particleCloud_.dataExchangeM().destroy(d2_,1); + particleCloud_.dataExchangeM().destroy(d3_,1); } // * * * * * * * * * * * * * * * private Member Functions * * * * * * * * * * * * * // + void dSauter::allocateMyArrays() const { // get memory for 2d arrays - double initVal=0.0; + double initVal = 0.0; particleCloud_.dataExchangeM().allocateArray(d2_,initVal,1); // field/initVal/with/lenghtFromLigghts particleCloud_.dataExchangeM().allocateArray(d3_,initVal,1); } + // * * * * * * * * * * * * * * * public Member Functions * * * * * * * * * * * * * // void dSauter::setForce() const { if (scaleDia_ > 1) + { Info << "dSauter using scaleCG = " << scaleDia_ << endl; + } else if (particleCloud_.cg() > 1) { - scaleDia_=particleCloud_.cg(); + scaleDia_ = particleCloud_.cg(); Info << "dSauter using scaleCG from liggghts cg = " << scaleDia_ << endl; } @@ -142,21 +146,21 @@ void dSauter::setForce() const label cellI=0; scalar ds(0); scalar scale = scaleDiaDist_/scaleDia_; - - for(int index = 0;index < particleCloud_.numberOfParticles(); ++index) + + for(int index = 0; index < particleCloud_.numberOfParticles(); ++index) { cellI = particleCloud_.cellIDs()[index][0]; - if(cellI >= 0) + if (cellI >= 0) { - ds = particleCloud_.d(index); - d2_[index][0] = ds*ds; - d3_[index][0] = ds*ds*ds; - } + ds = particleCloud_.d(index); + d2_[index][0] = ds*ds; + d3_[index][0] = ds*ds*ds; + } } - + d2Field_.primitiveFieldRef() = 0.0; d3Field_.primitiveFieldRef() = 0.0; - + particleCloud_.averagingM().setScalarSum ( d2Field_, @@ -164,7 +168,7 @@ void dSauter::setForce() const particleCloud_.particleWeights(), NULL ); - + particleCloud_.averagingM().setScalarSum ( d3Field_, @@ -172,19 +176,19 @@ void dSauter::setForce() const particleCloud_.particleWeights(), NULL ); - + forAll(dSauter_,cellI) { - if(d2Field_[cellI] > ROOTVSMALL) + if (d2Field_[cellI] > ROOTVSMALL) { - dSauter_[cellI] = d3Field_[cellI] / d2Field_[cellI] * scale; + dSauter_[cellI] = d3Field_[cellI] / d2Field_[cellI] * scale; } else - { - dSauter_[cellI] = SMALL; - } + { + dSauter_[cellI] = SMALL; + } } - + dSauter_.correctBoundaryConditions(); } diff --git a/src/lagrangian/cfdemParticle/subModels/forceModel/forceModel/forceModel.C b/src/lagrangian/cfdemParticle/subModels/forceModel/forceModel/forceModel.C index 3619330c..0b593924 100644 --- a/src/lagrangian/cfdemParticle/subModels/forceModel/forceModel/forceModel.C +++ b/src/lagrangian/cfdemParticle/subModels/forceModel/forceModel/forceModel.C @@ -88,14 +88,16 @@ forceModel::forceModel probeIt_(sm.probeM().active()), requiresEx_(false), forceSubModels_(0), - forceSubModel_(new autoPtr[nrForceSubModels()]) + forceSubModel_(NULL) {} // * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * // forceModel::~forceModel() -{} +{ + delete [] forceSubModel_; +} // * * * * * * * * * * * * * * * * Member Fct * * * * * * * * * * * * * * * // /*tmp forceModel::provideScalarField() diff --git a/src/lagrangian/cfdemParticle/subModels/forceModel/forceModel/forceModel.H b/src/lagrangian/cfdemParticle/subModels/forceModel/forceModel/forceModel.H index 99b7733c..204652eb 100644 --- a/src/lagrangian/cfdemParticle/subModels/forceModel/forceModel/forceModel.H +++ b/src/lagrangian/cfdemParticle/subModels/forceModel/forceModel/forceModel.H @@ -131,7 +131,7 @@ public: virtual void manipulateScalarField(volScalarField&) const; // Access - word modelType() { return modelType_; } + const word& modelType() const { return modelType_; } inline volVectorField& impParticleForces() const { return impParticleForces_; } @@ -147,7 +147,7 @@ public: inline double ** fluidVel() const { return particleCloud_.fluidVel_; } - inline const bool& coupleForce() const { return coupleForce_; } + inline bool coupleForce() const { return coupleForce_; } virtual inline bool& requiresEx() { return requiresEx_; } @@ -155,11 +155,11 @@ public: void treatVoidCells() const; - inline const wordList& forceSubModels(){ return forceSubModels_; } + inline const wordList& forceSubModels() { return forceSubModels_; } - inline const forceSubModel& forceSubM(int i) const { return forceSubModel_[i]; } + inline forceSubModel& forceSubM(int i) const { return forceSubModel_[i](); } - inline int nrForceSubModels(){ return forceSubModels_.size(); } + inline int nrForceSubModels() const { return forceSubModels_.size(); } void setForceSubModels(dictionary& dict); }; diff --git a/src/lagrangian/cfdemParticle/subModels/forceModel/forceSubModels/ScaleForce/ScaleForce.C b/src/lagrangian/cfdemParticle/subModels/forceModel/forceSubModels/ScaleForce/ScaleForce.C index 942d3a59..89527e65 100644 --- a/src/lagrangian/cfdemParticle/subModels/forceModel/forceSubModels/ScaleForce/ScaleForce.C +++ b/src/lagrangian/cfdemParticle/subModels/forceModel/forceSubModels/ScaleForce/ScaleForce.C @@ -62,7 +62,7 @@ ScaleForce::~ScaleForce() // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // void ScaleForce::partToArray ( - label& index, + label index, vector& dragTot, const vector& dragEx, const vector& Ufluid, diff --git a/src/lagrangian/cfdemParticle/subModels/forceModel/forceSubModels/ScaleForce/ScaleForce.H b/src/lagrangian/cfdemParticle/subModels/forceModel/forceSubModels/ScaleForce/ScaleForce.H index cd44fee5..47715b85 100644 --- a/src/lagrangian/cfdemParticle/subModels/forceModel/forceSubModels/ScaleForce/ScaleForce.H +++ b/src/lagrangian/cfdemParticle/subModels/forceModel/forceSubModels/ScaleForce/ScaleForce.H @@ -41,9 +41,9 @@ class ScaleForce public forceSubModel { private: - + dictionary propsDict_; - + word scaleFieldName_; const volScalarField& scaleField_; @@ -70,8 +70,8 @@ public: // Member Functions - - void partToArray(label&, vector&, const vector&, const vector& Ufluid=vector::zero, scalar Cd=scalar(0)) const; + + void partToArray(label, vector&, const vector&, const vector& Ufluid=vector::zero, scalar Cd=scalar(0)) const; word myType() const{return typeName; }; }; diff --git a/src/lagrangian/cfdemParticle/subModels/forceModel/forceSubModels/forceSubModel/forceSubModel.C b/src/lagrangian/cfdemParticle/subModels/forceModel/forceSubModels/forceSubModel/forceSubModel.C index ba235213..fe3e98df 100644 --- a/src/lagrangian/cfdemParticle/subModels/forceModel/forceSubModels/forceSubModel/forceSubModel.C +++ b/src/lagrangian/cfdemParticle/subModels/forceModel/forceSubModels/forceSubModel/forceSubModel.C @@ -59,7 +59,7 @@ forceSubModel::forceSubModel dict_(dict), particleCloud_(sm), forceModel_(fm), - nrDefaultSwitches_(9), // !!! + nrDefaultSwitches_(SW_MAX), switchesNameList_(nrDefaultSwitches_), switchesList_(nrDefaultSwitches_), switches_(nrDefaultSwitches_), @@ -106,16 +106,15 @@ forceSubModel::forceSubModel rho_(sm.mesh().lookupObject (densityFieldName_)) { // init standard switch list - int iCounter(0); - switchesNameList_[iCounter]="treatForceExplicit"; iCounter++; //0 - switchesNameList_[iCounter]="treatForceDEM";iCounter++; //1 - switchesNameList_[iCounter]="implForceDEM";iCounter++; //2 - switchesNameList_[iCounter]="verbose";iCounter++; //3 - switchesNameList_[iCounter]="interpolation";iCounter++; //4 - switchesNameList_[iCounter]="useFilteredDragModel";iCounter++; //5 - switchesNameList_[iCounter]="useParcelSizeDependentFilteredDrag";iCounter++; //6 - switchesNameList_[iCounter]="implForceDEMaccumulated";iCounter++; //7 - switchesNameList_[iCounter]="scalarViscosity";iCounter++; //8 + switchesNameList_[SW_TREAT_FORCE_EXPLICIT] = "treatForceExplicit"; + switchesNameList_[SW_TREAT_FORCE_DEM] = "treatForceDEM"; + switchesNameList_[SW_IMPL_FORCE_DEM] = "implForceDEM"; + switchesNameList_[SW_VERBOSE] = "verbose"; + switchesNameList_[SW_INTERPOLATION] = "interpolation"; + switchesNameList_[SW_FILTERED_DRAG_MODEL] ="useFilteredDragModel"; + switchesNameList_[SW_PARCEL_SIZE_DEPENDENT_FILTERED_DRAG] = "useParcelSizeDependentFilteredDrag"; + switchesNameList_[SW_IMPL_FORCE_DEM_ACCUMULATED] = "implForceDEMaccumulated"; + switchesNameList_[SW_SCALAR_VISCOSITY] = "scalarViscosity"; // sanity check of what is defined above if(switchesNameList_.size() != nrDefaultSwitches_) @@ -131,7 +130,7 @@ forceSubModel::~forceSubModel() // * * * * * * * * * * * * * * * * Member Fct * * * * * * * * * * * * * * * // void forceSubModel::partToArray ( - label& index, + label index, vector& dragTot, const vector& dragEx, const vector& Ufluid, @@ -139,17 +138,17 @@ void forceSubModel::partToArray ) const { // forces for CFD - if(!switches_[1])// !treatDEM + if(!switches_[SW_TREAT_FORCE_DEM])// !treatDEM { - if(switches_[0]) // treatExplicit + if(switches_[SW_TREAT_FORCE_EXPLICIT]) // treatExplicit { for(int j=0;j<3;j++) myForceM().expForces()[index][j] += dragTot[j]; - } + } else //implicit treatment, taking explicit force contribution into account { - for(int j=0;j<3;j++) - { + for(int j=0;j<3;j++) + { myForceM().impForces()[index][j] += dragTot[j] - dragEx[j]; //only consider implicit part! myForceM().expForces()[index][j] += dragEx[j]; } @@ -157,16 +156,16 @@ void forceSubModel::partToArray } // forces for DEM - if(switches_[2]) // implForceDEM + if(switches_[SW_IMPL_FORCE_DEM]) // implForceDEM { for(int j=0;j<3;j++) myForceM().fluidVel()[index][j]=Ufluid[j]; - myForceM().Cds()[index][0]=Cd; + myForceM().Cds()[index][0] = Cd; } else { - for(int j=0;j<3;j++) + for(int j=0;j<3;j++) myForceM().DEMForces()[index][j] += dragTot[j]; } } @@ -183,7 +182,7 @@ void forceSubModel::explicitCorr vector& Us, const vector& UsCell, bool verbose, - label index + label index ) const { dragExplicit=vector::zero; @@ -191,7 +190,7 @@ void forceSubModel::explicitCorr // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // -void forceSubModel::readSwitches() const +void forceSubModel::readSwitches() { Info << "\nreading switches for forceSubModel:" << myType() << endl; forAll(switchesNameList_,i) @@ -200,23 +199,23 @@ void forceSubModel::readSwitches() const { Info << " looking for " << switchesNameList_[i] << " ..." << endl; if (dict_.found(switchesNameList_[i])) - switches_[i]=Switch(dict_.lookup(switchesNameList_[i])); - + switches_[i] = Switch(dict_.lookup(switchesNameList_[i])); + Info << "\t" << switchesNameList_[i] << " = " << switches_[i] << endl; - } + } } Info << endl; - if(switches_[2]) // implForceDEM=true + if(switches_[SW_IMPL_FORCE_DEM]) // implForceDEM=true { // communicate implForceDEM to particleCloud - particleCloud_.impDEMdrag_=true; + particleCloud_.impDEMdrag_ = true; // do sanity check // This can work if the accumulator is used, but is explicitely applied on the CFD side // Sanity check is therefore not necessary here /* - if(switches_[0]) // treatExplicit=true + if(switches_[SW_TREAT_FORCE_EXPLICIT]) // treatExplicit=true { FatalError << "Please check your settings, treatExplicit together with implForceDEM does not work!." << abort(FatalError); @@ -224,24 +223,24 @@ void forceSubModel::readSwitches() const */ } - if(switches_[7]) // implForceDEMaccumulated=true + if(switches_[SW_IMPL_FORCE_DEM_ACCUMULATED]) // implForceDEMaccumulated=true { // sanity check for implForceDEMaccumulated - if(!switches_[2]) //implForceDEM=false + if(!switches_[SW_IMPL_FORCE_DEM]) //implForceDEM=false { Warning<< "please check your settings, implForceDEMaccumulated without implForceDEM does not work! (using implForceDEMaccumulated=false)" << endl; - switches_[3]=false; + switches_[SW_VERBOSE] = false; }else { - particleCloud_.impDEMdragAcc_=true; + particleCloud_.impDEMdragAcc_ = true; } } - if(switches_[8]) // scalarViscosity=true + if(switches_[SW_SCALAR_VISCOSITY]) // scalarViscosity=true { Info << "Using a constant viscosity for this force model." << endl; dimensionedScalar nu0_("nu", dimensionSet(0, 2, -1, 0, 0), dict_.lookup("nu")); - nu_=volScalarField + nu_ = volScalarField ( IOobject ( @@ -259,7 +258,7 @@ void forceSubModel::readSwitches() const // look for old nomenclature if (dict_.found("treatExplicit") || dict_.found("treatDEM") || dict_.found("implDEM")) FatalError<< "You are using an old nomenclature for force model settings, please have a look at the forceSubModel doc." << abort(FatalError); - + // look for old nomenclature if (dict_.found("verbose")) Warning<< "Please make sure you use the new nomenclature for verbose force model settings, please have a look at the forceSubModel doc." << endl; @@ -274,7 +273,7 @@ const volScalarField& forceSubModel::nuField() const nu_=particleCloud_.turbulence().mu() / rho_; return nu_; #else - if(switches_[8]) // scalarViscosity=true + if(switches_[SW_SCALAR_VISCOSITY]) // scalarViscosity=true return nu_; else return particleCloud_.turbulence().nu(); @@ -287,9 +286,9 @@ const volScalarField& forceSubModel::muField() const return particleCloud_.turbulence().mu(); #else // passing the ref to nu*rho will not work->generate a mu_ field like nu_ - FatalError<< "implementation not complete!" << abort(FatalError); + FatalError << "implementation not complete!" << abort(FatalError); - if(switches_[8]) // scalarViscosity=true + if(switches_[SW_SCALAR_VISCOSITY]) // scalarViscosity=true return nu_*rho_; else return particleCloud_.turbulence().nu()*rho_; diff --git a/src/lagrangian/cfdemParticle/subModels/forceModel/forceSubModels/forceSubModel/forceSubModel.H b/src/lagrangian/cfdemParticle/subModels/forceModel/forceSubModels/forceSubModel/forceSubModel.H index f64173bf..591c525a 100644 --- a/src/lagrangian/cfdemParticle/subModels/forceModel/forceSubModels/forceSubModel/forceSubModel.H +++ b/src/lagrangian/cfdemParticle/subModels/forceModel/forceSubModels/forceSubModel/forceSubModel.H @@ -47,6 +47,19 @@ SourceFiles namespace Foam { +enum { + SW_TREAT_FORCE_EXPLICIT = 0, + SW_TREAT_FORCE_DEM, + SW_IMPL_FORCE_DEM, + SW_VERBOSE, + SW_INTERPOLATION, + SW_FILTERED_DRAG_MODEL, + SW_PARCEL_SIZE_DEPENDENT_FILTERED_DRAG, + SW_IMPL_FORCE_DEM_ACCUMULATED, + SW_SCALAR_VISCOSITY, + SW_MAX +}; + /*---------------------------------------------------------------------------*\ Class forceSubModel Declaration \*---------------------------------------------------------------------------*/ @@ -67,9 +80,9 @@ protected: wordList switchesNameList_; // names of switches available - mutable List switchesList_; // switches which are requested in dict + List switchesList_; // switches which are requested from dict - mutable List switches_; + List switches_; // switch status mutable volScalarField nu_; @@ -130,35 +143,33 @@ public: // Member Functions - virtual void partToArray(label&, vector&, const vector&, const vector& Ufluid=vector::zero, scalar Cd=scalar(0)) const; + virtual void partToArray(label, vector&, const vector&, const vector& Ufluid=vector::zero, scalar Cd=scalar(0)) const; virtual void explicitCorr(vector&, vector&, scalar&, vector&, const vector&, vector&, const vector&, bool,label index=100) const; // Access - inline bool verbose() const { return switches_[3]; } + inline bool verbose() const { return switches_[SW_VERBOSE]; } - inline bool interpolation() const { return switches_[4]; } + inline bool interpolation() const { return switches_[SW_INTERPOLATION]; } - inline bool useFilteredDragModel() const { return switches_[5]; } + inline bool useFilteredDragModel() const { return switches_[SW_FILTERED_DRAG_MODEL]; } - inline bool useParcelSizeDependentFilteredDrag() const { return switches_[6]; } + inline bool useParcelSizeDependentFilteredDrag() const { return switches_[SW_PARCEL_SIZE_DEPENDENT_FILTERED_DRAG]; } - virtual word myType() const=0; + virtual word myType() const = 0; - inline forceModel& myForceM() const { return forceModel_;} + inline forceModel& myForceM() const { return forceModel_; } - inline const List& switches() const { return switches_;} + inline const List& switches() const { return switches_; } - inline const wordList& switchesNameList() const { return switchesNameList_;} + inline const wordList& switchesNameList() const { return switchesNameList_; } - void setSwitchesList(label i, bool v) const { switchesList_[i] = v;} + void setSwitchesList(label i, Switch v) { switchesList_[i] = v; } - void setSwitches(label i, Switch v) const { switches_[i] = v;} + void setSwitches(label i, Switch v) { switches_[i] = v; } - virtual void readSwitches() const; - - const label& nrDefaultSwitches() const {return nrDefaultSwitches_;} + virtual void readSwitches(); const volScalarField& nuField() const; diff --git a/src/lagrangian/cfdemParticle/subModels/forceModel/forceSubModels/scaleForceBoundary/scaleForceBoundary.C b/src/lagrangian/cfdemParticle/subModels/forceModel/forceSubModels/scaleForceBoundary/scaleForceBoundary.C index 0270881b..aa2b0639 100644 --- a/src/lagrangian/cfdemParticle/subModels/forceModel/forceSubModels/scaleForceBoundary/scaleForceBoundary.C +++ b/src/lagrangian/cfdemParticle/subModels/forceModel/forceSubModels/scaleForceBoundary/scaleForceBoundary.C @@ -56,35 +56,35 @@ scaleForceBoundary::scaleForceBoundary if (propsDict_.found("x1") && propsDict_.found("x2")) { coordinateInner_=readScalar(propsDict_.lookup ("x1")); - coordinateOuter_=readScalar(propsDict_.lookup ("x2")); - dim_ = 0; + coordinateOuter_=readScalar(propsDict_.lookup ("x2")); + dim_ = 0; Info << "scaleForceBoundary: Limiting force in x direction." << endl; } else if (propsDict_.found("y1") && propsDict_.found("y2")) { coordinateInner_=readScalar(propsDict_.lookup ("y1")); - coordinateOuter_=readScalar(propsDict_.lookup ("y2")); - dim_ = 1; + coordinateOuter_=readScalar(propsDict_.lookup ("y2")); + dim_ = 1; Info << "scaleForceBoundary: Limiting force in y direction." << endl; } else if (propsDict_.found("z1") && propsDict_.found("z2")) { coordinateInner_=readScalar(propsDict_.lookup ("z1")); - coordinateOuter_=readScalar(propsDict_.lookup ("z2")); - dim_ = 2; + coordinateOuter_=readScalar(propsDict_.lookup ("z2")); + dim_ = 2; Info << "scaleForceBoundary: Limiting force in z direction." << endl; } - + if (propsDict_.found("outerValue")) outerVal_=readScalar(propsDict_.lookup ("outerValue")); - + if(coordinateOuter_ > coordinateInner_) orientation_ = 1; else orientation_ = -1; - + dist_ = fabs(coordinateOuter_ - coordinateInner_); - + } @@ -97,7 +97,7 @@ scaleForceBoundary::~scaleForceBoundary() // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // void scaleForceBoundary::partToArray ( - label& index, + label index, vector& dragTot, const vector& dragEx, const vector& Ufluid, diff --git a/src/lagrangian/cfdemParticle/subModels/forceModel/forceSubModels/scaleForceBoundary/scaleForceBoundary.H b/src/lagrangian/cfdemParticle/subModels/forceModel/forceSubModels/scaleForceBoundary/scaleForceBoundary.H index 807c860f..8a4f65cd 100644 --- a/src/lagrangian/cfdemParticle/subModels/forceModel/forceSubModels/scaleForceBoundary/scaleForceBoundary.H +++ b/src/lagrangian/cfdemParticle/subModels/forceModel/forceSubModels/scaleForceBoundary/scaleForceBoundary.H @@ -45,19 +45,19 @@ class scaleForceBoundary public forceSubModel { private: - + dictionary propsDict_; - + scalar coordinateInner_; - + scalar coordinateOuter_; - + scalar dist_; - + scalar outerVal_; - + int dim_; - + int orientation_; public: @@ -82,10 +82,10 @@ public: // Member Functions - - void partToArray(label&, vector&, const vector&, const vector& Ufluid=vector::zero, scalar Cd=scalar(0)) const; - word myType() const{return typeName; }; + void partToArray(label, vector&, const vector&, const vector& Ufluid=vector::zero, scalar Cd=scalar(0)) const; + + word myType() const {return typeName; } }; diff --git a/src/lagrangian/cfdemParticle/subModels/forceModel/gradPForce/gradPForce.C b/src/lagrangian/cfdemParticle/subModels/forceModel/gradPForce/gradPForce.C index 1885294c..586579b7 100644 --- a/src/lagrangian/cfdemParticle/subModels/forceModel/gradPForce/gradPForce.C +++ b/src/lagrangian/cfdemParticle/subModels/forceModel/gradPForce/gradPForce.C @@ -74,9 +74,9 @@ gradPForce::gradPForce setForceSubModels(propsDict_); // define switches which can be read from dict - forceSubM(0).setSwitchesList(0,true); // activate treatExplicit switch - forceSubM(0).setSwitchesList(1,true); // activate treatForceDEM switch - forceSubM(0).setSwitchesList(4,true); // activate search for interpolate switch + forceSubM(0).setSwitchesList(SW_TREAT_FORCE_EXPLICIT,true); // activate treatExplicit switch + forceSubM(0).setSwitchesList(SW_TREAT_FORCE_DEM,true); // activate treatForceDEM switch + forceSubM(0).setSwitchesList(SW_INTERPOLATION,true); // activate search for interpolate switch // read those switches defined above, if provided in dict forceSubM(0).readSwitches(); @@ -84,26 +84,30 @@ gradPForce::gradPForce if (modelType_ == "B") { FatalError <<"using model gradPForce with model type B is not valid\n" << abort(FatalError); - }else if (modelType_ == "Bfull") + } + else if (modelType_ == "Bfull") { - if(forceSubM(0).switches()[1]) + if(forceSubM(0).switches()[SW_TREAT_FORCE_DEM]) { Info << "Using treatForceDEM false!" << endl; - forceSubM(0).setSwitches(1,false); // treatForceDEM = false + forceSubM(0).setSwitches(SW_TREAT_FORCE_DEM,false); // treatForceDEM = false } - }else // modelType_=="A" + } + else // modelType_=="A" { - if(!forceSubM(0).switches()[1]) + if(!forceSubM(0).switches()[SW_TREAT_FORCE_DEM]) { Info << "Using treatForceDEM true!" << endl; - forceSubM(0).setSwitches(1,true); // treatForceDEM = true + forceSubM(0).setSwitches(SW_TREAT_FORCE_DEM,true); // treatForceDEM = true } } - if (propsDict_.found("useU")) useU_=true; - if (propsDict_.found("useAddedMass")) + if (propsDict_.found("useU")) + useU_ = true; + + if (propsDict_.found("useAddedMass")) { - addedMassCoeff_ = readScalar(propsDict_.lookup("useAddedMass")); + addedMassCoeff_ = readScalar(propsDict_.lookup("useAddedMass")); Info << "gradP will also include added mass with coefficient: " << addedMassCoeff_ << endl; Info << "WARNING: use fix nve/sphere/addedMass in LIGGGHTS input script to correctly account for added mass effects!" << endl; } @@ -113,7 +117,7 @@ gradPForce::gradPForce particleCloud_.checkCG(true); - particleCloud_.probeM().initialize(typeName, "gradP.logDat"); + particleCloud_.probeM().initialize(typeName, typeName+".logDat"); particleCloud_.probeM().vectorFields_.append("gradPForce"); //first entry must the be the force particleCloud_.probeM().scalarFields_.append("Vs"); particleCloud_.probeM().scalarFields_.append("rho"); @@ -152,21 +156,22 @@ void gradPForce::setForce() const #include "setupProbeModel.H" - for(int index = 0;index < particleCloud_.numberOfParticles(); index++) + for(int index = 0; index < particleCloud_.numberOfParticles(); ++index) { //if(mask[index][0]) //{ - force=vector(0,0,0); + force = vector(0,0,0); cellI = particleCloud_.cellIDs()[index][0]; if (cellI > -1) // particle Found { position = particleCloud_.position(index); - if(forceSubM(0).interpolation()) // use intepolated values for alpha (normally off!!!) + if (forceSubM(0).interpolation()) // use intepolated values for alpha (normally off!!!) { gradP = gradPInterpolator_.interpolate(position,cellI); - }else + } + else { gradP = gradPField[cellI]; } @@ -180,7 +185,7 @@ void gradPForce::setForce() const else force = -Vs*gradP*(1.0+addedMassCoeff_); - if(forceSubM(0).verbose() && index >=0 && index <2) + if (forceSubM(0).verbose() && index >= 0 && index < 2) { Info << "index = " << index << endl; Info << "gradP = " << gradP << endl; @@ -188,7 +193,7 @@ void gradPForce::setForce() const } //Set value fields and write the probe - if(probeIt_) + if (probeIt_) { #include "setupProbeModelfields.H" vValues.append(force); //first entry must the be the force diff --git a/src/lagrangian/cfdemParticle/subModels/forceModel/granKineticEnergy/granKineticEnergy.C b/src/lagrangian/cfdemParticle/subModels/forceModel/granKineticEnergy/granKineticEnergy.C index 350b20cb..072ea2de 100644 --- a/src/lagrangian/cfdemParticle/subModels/forceModel/granKineticEnergy/granKineticEnergy.C +++ b/src/lagrangian/cfdemParticle/subModels/forceModel/granKineticEnergy/granKineticEnergy.C @@ -67,12 +67,12 @@ granKineticEnergy::granKineticEnergy ), sm.mesh(), dimensionedScalar("zero", dimensionSet(0,2,-2,0,0), 0), - "zeroGradient" + "zeroGradient" ) { allocateMyArrays(); granKineticEnergy_.write(); - + // init force sub model setForceSubModels(propsDict_); @@ -82,14 +82,14 @@ granKineticEnergy::granKineticEnergy granKineticEnergy::~granKineticEnergy() { - delete vfluc_; + particleCloud_.dataExchangeM().destroy(vfluc_,1); } // * * * * * * * * * * * * * * * private Member Functions * * * * * * * * * * * * * // void granKineticEnergy::allocateMyArrays() const { // get memory for 2d arrays - double initVal=0.0; + double initVal = 0.0; particleCloud_.dataExchangeM().allocateArray(vfluc_,initVal,1); } // * * * * * * * * * * * * * * * public Member Functions * * * * * * * * * * * * * // @@ -98,20 +98,20 @@ void granKineticEnergy::setForce() const { allocateMyArrays(); - label cellI=0; + label cellI = 0; vector velfluc(0,0,0); - - for(int index = 0;index < particleCloud_.numberOfParticles(); ++index) + + for(int index = 0; index < particleCloud_.numberOfParticles(); ++index) { cellI = particleCloud_.cellIDs()[index][0]; if(cellI >= 0) { - velfluc = particleCloud_.velocity(index) - UsField_[cellI]; - vfluc_[index][0] = magSqr(velfluc); - } + velfluc = particleCloud_.velocity(index) - UsField_[cellI]; + vfluc_[index][0] = magSqr(velfluc); + } } - + granKineticEnergy_.primitiveFieldRef() = 0.0; particleCloud_.averagingM().resetWeightFields(); @@ -120,10 +120,10 @@ void granKineticEnergy::setForce() const granKineticEnergy_, vfluc_, particleCloud_.particleWeights(), - particleCloud_.averagingM().UsWeightField(), + particleCloud_.averagingM().UsWeightField(), NULL ); - + granKineticEnergy_ *= 0.5; } diff --git a/src/lagrangian/cfdemParticle/subModels/forceModel/interface/interface.C b/src/lagrangian/cfdemParticle/subModels/forceModel/interface/interface.C index 3507f08d..ad77200e 100644 --- a/src/lagrangian/cfdemParticle/subModels/forceModel/interface/interface.C +++ b/src/lagrangian/cfdemParticle/subModels/forceModel/interface/interface.C @@ -83,7 +83,7 @@ interface::interface setForceSubModels(propsDict_); // define switches which can be read from dict - forceSubM(0).setSwitchesList(0,true); // activate treatExplicit switch + forceSubM(0).setSwitchesList(SW_TREAT_FORCE_EXPLICIT,true); // activate treatExplicit switch // read those switches defined above, if provided in dict forceSubM(0).readSwitches(); diff --git a/src/lagrangian/cfdemParticle/subModels/forceModel/noDrag/noDrag.C b/src/lagrangian/cfdemParticle/subModels/forceModel/noDrag/noDrag.C index 0a403a02..f431f890 100644 --- a/src/lagrangian/cfdemParticle/subModels/forceModel/noDrag/noDrag.C +++ b/src/lagrangian/cfdemParticle/subModels/forceModel/noDrag/noDrag.C @@ -72,7 +72,7 @@ noDrag::noDrag setForceSubModels(propsDict_); // define switches which can be read from dict - forceSubM(0).setSwitchesList(0,true); // activate treatExplicit + forceSubM(0).setSwitchesList(SW_TREAT_FORCE_EXPLICIT,true); // activate treatExplicit // read those switches defined above, if provided in dict forceSubM(0).readSwitches(); @@ -98,7 +98,7 @@ void noDrag::setForce() const // Do nothing Info << "noDrag::setForce" << endl; label cellI=0; - bool treatExplicit=forceSubM(0).switches()[0]; + bool treatExplicit = forceSubM(0).switches()[SW_TREAT_FORCE_EXPLICIT]; for(int index = 0;index < particleCloud_.numberOfParticles(); ++index) { cellI = particleCloud_.cellIDs()[index][0]; @@ -125,7 +125,7 @@ void noDrag::setForce() const } } //========================== - } + } } } diff --git a/src/lagrangian/cfdemParticle/subModels/forceModel/virtualMassForce/virtualMassForce.C b/src/lagrangian/cfdemParticle/subModels/forceModel/virtualMassForce/virtualMassForce.C index 5443d948..2d33a4fb 100644 --- a/src/lagrangian/cfdemParticle/subModels/forceModel/virtualMassForce/virtualMassForce.C +++ b/src/lagrangian/cfdemParticle/subModels/forceModel/virtualMassForce/virtualMassForce.C @@ -83,28 +83,30 @@ virtualMassForce::virtualMassForce // init force sub model setForceSubModels(propsDict_); // define switches which can be read from dict - forceSubM(0).setSwitchesList(0,true); // activate treatExplicit switch - forceSubM(0).setSwitchesList(4,true); // activate search for interpolate switch + forceSubM(0).setSwitchesList(SW_TREAT_FORCE_EXPLICIT,true); // activate treatExplicit switch + forceSubM(0).setSwitchesList(SW_INTERPOLATION,true); // activate search for interpolate switch forceSubM(0).readSwitches(); //Extra switches/settings if(propsDict_.found("splitUrelCalculation")) { - splitUrelCalculation_ = readBool(propsDict_.lookup("splitUrelCalculation")); + splitUrelCalculation_ = readBool(propsDict_.lookup("splitUrelCalculation")); if(splitUrelCalculation_) + { Info << "Virtual mass model: will split the Urel calculation\n"; Info << "WARNING: be sure that LIGGGHTS integration takes ddtv_p implicitly into account! \n"; + } } if(propsDict_.found("Cadd")) { - Cadd_ = readScalar(propsDict_.lookup("Cadd")); + Cadd_ = readScalar(propsDict_.lookup("Cadd")); Info << "Virtual mass model: using non-standard Cadd = " << Cadd_ << endl; } particleCloud_.checkCG(true); //Append the field names to be probed - particleCloud_.probeM().initialize(typeName, "virtualMass.logDat"); + particleCloud_.probeM().initialize(typeName, typeName+".logDat"); particleCloud_.probeM().vectorFields_.append("virtualMassForce"); //first entry must the be the force particleCloud_.probeM().vectorFields_.append("Urel"); particleCloud_.probeM().vectorFields_.append("UrelOld"); @@ -119,7 +121,7 @@ virtualMassForce::virtualMassForce virtualMassForce::~virtualMassForce() { - delete UrelOld_; + particleCloud_.dataExchangeM().destroy(UrelOld_,3); } @@ -146,7 +148,7 @@ void virtualMassForce::setForce() const #include "setupProbeModel.H" - bool haveUrelOld_(false); + bool haveUrelOld_(false); for(int index = 0;index < particleCloud_.numberOfParticles(); index++) { @@ -156,19 +158,19 @@ void virtualMassForce::setForce() const if (cellI > -1) // particle Found { - if(forceSubM(0).interpolation()) + if(forceSubM(0).interpolation()) { - position = particleCloud_.position(index); - Ufluid = UInterpolator_.interpolate(position,cellI); + position = particleCloud_.position(index); + Ufluid = UInterpolator_.interpolate(position,cellI); } else { Ufluid = U_[cellI]; } - + if(splitUrelCalculation_) //if split, just use total derivative of fluid velocity - if(forceSubM(0).interpolation()) + if(forceSubM(0).interpolation()) { DDtU = DDtUInterpolator_.interpolate(position,cellI); } @@ -182,7 +184,7 @@ void virtualMassForce::setForce() const Ur = Ufluid - Us; } - + //Check of particle was on this CPU the last step if(UrelOld_[index][0]==NOTONCPU) //use 1. element to indicate that particle was on this CPU the last time step haveUrelOld_ = false; @@ -229,7 +231,7 @@ void virtualMassForce::setForce() const } else //particle not on this CPU UrelOld_[index][0]=NOTONCPU; - + // write particle based data to global array forceSubM(0).partToArray(index,virtualMassForce,vector::zero); } diff --git a/src/lagrangian/cfdemParticle/subModels/forceModel/viscForce/viscForce.C b/src/lagrangian/cfdemParticle/subModels/forceModel/viscForce/viscForce.C index 571915ab..30762ce8 100644 --- a/src/lagrangian/cfdemParticle/subModels/forceModel/viscForce/viscForce.C +++ b/src/lagrangian/cfdemParticle/subModels/forceModel/viscForce/viscForce.C @@ -71,10 +71,10 @@ viscForce::viscForce setForceSubModels(propsDict_); // define switches which can be read from dict - forceSubM(0).setSwitchesList(0,true); // activate treatExplicit switch - forceSubM(0).setSwitchesList(1,true); // activate treatForceDEM switch - forceSubM(0).setSwitchesList(4,true); // activate search for interpolate switch - forceSubM(0).setSwitchesList(8,true); // activate scalarViscosity switch + forceSubM(0).setSwitchesList(SW_TREAT_FORCE_EXPLICIT,true); // activate treatExplicit switch + forceSubM(0).setSwitchesList(SW_TREAT_FORCE_DEM,true); // activate treatForceDEM switch + forceSubM(0).setSwitchesList(SW_INTERPOLATION,true); // activate search for interpolate switch + forceSubM(0).setSwitchesList(SW_SCALAR_VISCOSITY,true); // activate scalarViscosity switch // read those switches defined above, if provided in dict forceSubM(0).readSwitches(); @@ -84,22 +84,22 @@ viscForce::viscForce FatalError <<"using model viscForce with model type B is not valid\n" << abort(FatalError); }else if (modelType_ == "Bfull") { - if(forceSubM(0).switches()[1]) + if(forceSubM(0).switches()[SW_TREAT_FORCE_DEM]) { Info << "Using treatForceDEM false!" << endl; - forceSubM(0).setSwitches(1,false); // treatForceDEM = false + forceSubM(0).setSwitches(SW_TREAT_FORCE_DEM,false); // treatForceDEM = false } }else // modelType_=="A" { - if(!forceSubM(0).switches()[1]) + if(!forceSubM(0).switches()[SW_TREAT_FORCE_DEM]) { Info << "Using treatForceDEM true!" << endl; - forceSubM(0).setSwitches(1,true); // treatForceDEM = true + forceSubM(0).setSwitches(SW_TREAT_FORCE_DEM,true); // treatForceDEM = true } } - if (propsDict_.found("useAddedMass")) + if (propsDict_.found("useAddedMass")) { addedMassCoeff_ = readScalar(propsDict_.lookup("useAddedMass")); Info << "viscForce will also include added mass with coefficient: " << addedMassCoeff_ << endl; @@ -109,7 +109,7 @@ viscForce::viscForce particleCloud_.checkCG(true); //Append the field names to be probed - particleCloud_.probeM().initialize(typeName, "visc.logDat"); + particleCloud_.probeM().initialize(typeName, typeName+".logDat"); particleCloud_.probeM().vectorFields_.append("viscForce"); //first entry must the be the force particleCloud_.probeM().scalarFields_.append("Vs"); particleCloud_.probeM().writeHeader(); @@ -160,7 +160,7 @@ void viscForce::setForce() const Vs = particleCloud_.particleVolume(index); - // calc the contribution of the deviatoric stress + // calc the contribution of the deviatoric stress // to the generalized buoyancy force force = -Vs*divTau*(1.0+addedMassCoeff_); diff --git a/src/lagrangian/cfdemParticle/subModels/forceModelMS/DiFeliceDragMS/DiFeliceDragMS.C b/src/lagrangian/cfdemParticle/subModels/forceModelMS/DiFeliceDragMS/DiFeliceDragMS.C index 27b4f278..7cb38d76 100644 --- a/src/lagrangian/cfdemParticle/subModels/forceModelMS/DiFeliceDragMS/DiFeliceDragMS.C +++ b/src/lagrangian/cfdemParticle/subModels/forceModelMS/DiFeliceDragMS/DiFeliceDragMS.C @@ -75,21 +75,21 @@ DiFeliceDragMS::DiFeliceDragMS //dH_(readScalar(propsDict_.lookup("hydraulicDiameter"))) { //Append the field names to be probed - particleCloud_.probeM().initialize(typeName, "diFeliceDrag.logDat"); - particleCloud_.probeM().vectorFields_.append("dragForce"); //first entry must the be the force - particleCloud_.probeM().vectorFields_.append("Urel"); //other are debug - particleCloud_.probeM().scalarFields_.append("Rep"); //other are debug - particleCloud_.probeM().scalarFields_.append("Cd"); //other are debug - particleCloud_.probeM().scalarFields_.append("voidfraction"); //other are debug + particleCloud_.probeM().initialize(typeName, typeName+".logDat"); + particleCloud_.probeM().vectorFields_.append("dragForce"); // first entry must the be the force + particleCloud_.probeM().vectorFields_.append("Urel"); // other are debug + particleCloud_.probeM().scalarFields_.append("Rep"); // other are debug + particleCloud_.probeM().scalarFields_.append("Cd"); // other are debug + particleCloud_.probeM().scalarFields_.append("voidfraction"); // other are debug particleCloud_.probeM().writeHeader(); // init force sub model setForceSubModels(propsDict_); // define switches which can be read from dict - forceSubM(0).setSwitchesList(0,true); // activate treatExplicit switch - forceSubM(0).setSwitchesList(3,true); // activate search for verbose switch - forceSubM(0).setSwitchesList(4,true); // activate search for interpolate switch + forceSubM(0).setSwitchesList(SW_TREAT_FORCE_EXPLICIT,true); // activate treatExplicit switch + forceSubM(0).setSwitchesList(SW_VERBOSE,true); // activate search for verbose switch + forceSubM(0).setSwitchesList(SW_INTERPOLATION,true); // activate search for interpolate switch // read those switches defined above, if provided in dict forceSubM(0).readSwitches(); @@ -98,11 +98,12 @@ DiFeliceDragMS::DiFeliceDragMS { Warning << " interpolation is commented for this force model - it seems to be unstable with AMI!" << endl; } + if (propsDict_.found("splitImplicitExplicit")) { Info << "will split implicit / explicit force contributions." << endl; splitImplicitExplicit_ = true; - if(!forceSubM(0).interpolation()) + if(!forceSubM(0).interpolation()) Info << "WARNING: will only consider fluctuating particle velocity in implicit / explicit force split!" << endl; } particleCloud_.checkCG(false); @@ -124,11 +125,11 @@ void DiFeliceDragMS::setForce() const const volScalarField& nufField = forceSubM(0).nuField(); const volScalarField& rhoField = forceSubM(0).rhoField(); - vector position(0,0,0); + //vector position(0,0,0); scalar voidfraction(1); vector Ufluid(0,0,0); vector drag(0,0,0); - label cellI=0; + label cellI = 0; vector Us(0,0,0); vector Ur(0,0,0); scalar ds(0); @@ -138,17 +139,17 @@ void DiFeliceDragMS::setForce() const scalar Rep(0); scalar Cd(0); - vector UfluidFluct(0,0,0); + vector UfluidFluct(0,0,0); vector UsFluct(0,0,0); vector dragExplicit(0,0,0); - scalar dragCoefficient(0); + scalar dragCoefficient(0); //interpolationCellPoint voidfractionInterpolator_(voidfraction_); //interpolationCellPoint UInterpolator_(U_); #include "setupProbeModel.H" - for(int index = 0;index < cloudRefMS().numberOfClumps(); index++) + for(int index = 0; index < cloudRefMS().numberOfClumps(); ++index) { //if(mask[index][0]) // would have to be transformed from body ID to particle ID @@ -214,7 +215,7 @@ void DiFeliceDragMS::setForce() const } } - if(forceSubM(0).verbose() && index >=0 && index <10) + if(forceSubM(0).verbose() && index >= 0 && index < 10) { Pout << "index = " << index << endl; Pout << "Us = " << Us << endl; @@ -246,8 +247,13 @@ void DiFeliceDragMS::setForce() const particleCloud_.probeM().writeProbe(index, sValues, vValues); } } + // set force on bodies - if(forceSubM(0).switches()[0]) for(int j=0;j<3;j++) cloudRefMS().expForcesCM()[index][j] += drag[j]; + if (forceSubM(0).switches()[SW_TREAT_FORCE_EXPLICIT]) + { + for(int j=0;j<3;j++) + cloudRefMS().expForcesCM()[index][j] += drag[j]; + } else //implicit treatment, taking explicit force contribution into account { for(int j=0;j<3;j++) @@ -256,7 +262,9 @@ void DiFeliceDragMS::setForce() const cloudRefMS().expForcesCM()[index][j] += dragExplicit[j]; } } - for(int j=0;j<3;j++) cloudRefMS().DEMForcesCM()[index][j] += drag[j]; + + for(int j=0;j<3;j++) + cloudRefMS().DEMForcesCM()[index][j] += drag[j]; //} } diff --git a/src/lagrangian/cfdemParticle/subModels/forceModelMS/DiFeliceDragMS/DiFeliceDragMS.H b/src/lagrangian/cfdemParticle/subModels/forceModelMS/DiFeliceDragMS/DiFeliceDragMS.H index 27bc6275..4a0efaf6 100644 --- a/src/lagrangian/cfdemParticle/subModels/forceModelMS/DiFeliceDragMS/DiFeliceDragMS.H +++ b/src/lagrangian/cfdemParticle/subModels/forceModelMS/DiFeliceDragMS/DiFeliceDragMS.H @@ -71,7 +71,7 @@ private: bool splitImplicitExplicit_; // use splitting of implicit and explict force contribution - word UsFieldName_; + word UsFieldName_; const volVectorField& UsField_; // the average particle velocity field (for implicit/expliti force split) diff --git a/src/lagrangian/cfdemParticle/subModels/momCoupleModel/explicitCouple/explicitCouple.C b/src/lagrangian/cfdemParticle/subModels/momCoupleModel/explicitCouple/explicitCouple.C index a5512441..190b6c85 100644 --- a/src/lagrangian/cfdemParticle/subModels/momCoupleModel/explicitCouple/explicitCouple.C +++ b/src/lagrangian/cfdemParticle/subModels/momCoupleModel/explicitCouple/explicitCouple.C @@ -108,38 +108,15 @@ explicitCouple::~explicitCouple() // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // tmp explicitCouple::expMomSource() const { - tmp tsource - ( - new volVectorField - ( - IOobject - ( - "f_explicitCouple", - particleCloud_.mesh().time().timeName(), - particleCloud_.mesh(), - IOobject::NO_READ, - IOobject::NO_WRITE - ), - particleCloud_.mesh(), - dimensionedVector - ( - "zero", - dimensionSet(1, -2, -2, 0, 0), // N/m3 - vector::zero - ), - "zeroGradient" - ) - ); - scalar tsf = particleCloud_.dataExchangeM().timeStepFraction(); - if(1-tsf < 1e-4) //tsf==1 + if (1. - tsf < 1e-4) //tsf==1 { // calc fNext forAll(fNext_,cellI) { fNext_[cellI] = arrayToField(cellI); - + // limiter for (int i=0;i<3;i++) { @@ -147,21 +124,27 @@ tmp explicitCouple::expMomSource() const if (magF > fLimit_[i]) fNext_[cellI][i] *= fLimit_[i]/magF; } } - tsource.ref() = fPrev_; - }else - { - tsource.ref() = (1 - tsf) * fPrev_ + tsf * fNext_; + return tmp + ( + new volVectorField("f_explicitCouple", fPrev_) + ); + } + else + { + return tmp + ( + new volVectorField("f_explicitCouple", (1. - tsf) * fPrev_ + tsf * fNext_) + ); } - return tsource; } -void Foam::explicitCouple::resetMomSourceField() const +void explicitCouple::resetMomSourceField() const { fPrev_.ref() = fNext_.ref(); fNext_.primitiveFieldRef() = vector::zero; } -inline vector Foam::explicitCouple::arrayToField(label cellI) const +inline vector explicitCouple::arrayToField(label cellI) const { return particleCloud_.forceM(0).expParticleForces()[cellI] / particleCloud_.mesh().V()[cellI]; } diff --git a/src/lagrangian/cfdemParticle/subModels/momCoupleModel/implicitCouple/implicitCouple.C b/src/lagrangian/cfdemParticle/subModels/momCoupleModel/implicitCouple/implicitCouple.C index 377fc137..124b7e9b 100644 --- a/src/lagrangian/cfdemParticle/subModels/momCoupleModel/implicitCouple/implicitCouple.C +++ b/src/lagrangian/cfdemParticle/subModels/momCoupleModel/implicitCouple/implicitCouple.C @@ -122,40 +122,19 @@ implicitCouple::~implicitCouple() tmp implicitCouple::impMomSource() const { - tmp tsource - ( - new volScalarField - ( - IOobject - ( - "Ksl_implicitCouple", - particleCloud_.mesh().time().timeName(), - particleCloud_.mesh(), - IOobject::NO_READ, - IOobject::NO_WRITE - ), - particleCloud_.mesh(), - dimensionedScalar - ( - "zero", - dimensionSet(1, -3, -1, 0, 0), // N/m3 / m/s - 0 - ) - ) - ); - scalar tsf = particleCloud_.dataExchangeM().timeStepFraction(); // calc Ksl - scalar Ur; - if(1-tsf < 1e-4) //tsf==1 + if (1. - tsf < 1e-4) //tsf==1 { + scalar Ur; + forAll(KslNext_,cellI) { Ur = mag(U_[cellI] - Us_[cellI]); - if(Ur > SMALL && alpha_[cellI] < maxAlpha_) //momentum exchange switched off if alpha too big + if (Ur > SMALL && alpha_[cellI] < maxAlpha_) //momentum exchange switched off if alpha too big { KslNext_[cellI] = mag(particleCloud_.forceM(0).impParticleForces()[cellI]) / Ur @@ -166,16 +145,21 @@ tmp implicitCouple::impMomSource() const // limiter if (KslNext_[cellI] > KslLimit_) KslNext_[cellI] = KslLimit_; } - tsource.ref() = KslPrev_; - }else - { - tsource.ref() = (1 - tsf) * KslPrev_ + tsf * KslNext_; + return tmp + ( + new volScalarField("Ksl_implicitCouple", KslPrev_) + ); + } + else + { + return tmp + ( + new volScalarField("Ksl_implicitCouple", (1. - tsf) * KslPrev_ + tsf * KslNext_) + ); } - - return tsource; } -void Foam::implicitCouple::resetMomSourceField() const +void implicitCouple::resetMomSourceField() const { KslPrev_.ref() = KslNext_.ref(); KslNext_.primitiveFieldRef() = 0; diff --git a/src/lagrangian/cfdemParticle/subModels/momCoupleModel/noCouple/noCouple.C b/src/lagrangian/cfdemParticle/subModels/momCoupleModel/noCouple/noCouple.C index 789907df..ffd629ec 100644 --- a/src/lagrangian/cfdemParticle/subModels/momCoupleModel/noCouple/noCouple.C +++ b/src/lagrangian/cfdemParticle/subModels/momCoupleModel/noCouple/noCouple.C @@ -72,7 +72,7 @@ noCouple::~noCouple() // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // -void Foam::noCouple::resetMomSourceField() const +void noCouple::resetMomSourceField() const { FatalError<<"the solver calls for resetMomSourceField() although you use IB method where this is not needed!\n" <<", check your solver! - PANIC -\n"; diff --git a/src/lagrangian/cfdemParticle/subModels/probeModel/particleProbe/particleProbe.C b/src/lagrangian/cfdemParticle/subModels/probeModel/particleProbe/particleProbe.C index c9bc6826..ca5bab6e 100644 --- a/src/lagrangian/cfdemParticle/subModels/probeModel/particleProbe/particleProbe.C +++ b/src/lagrangian/cfdemParticle/subModels/probeModel/particleProbe/particleProbe.C @@ -87,12 +87,12 @@ particleProbe::particleProbe if (propsDict_.found("verbose")) verbose_=true; if (propsDict_.found("verboseToFile")) verboseToFile_=true; - if (propsDict_.found("printEvery")) printEvery_= readScalar(propsDict_.lookup("printEvery")); - if (propsDict_.found("sampleAll")) sampleAll_=true; - if (propsDict_.found("probeDebug")) probeDebug_=true; - if (propsDict_.found("includePosition")) includePosition_=true; + if (propsDict_.found("printEvery")) printEvery_ = readScalar(propsDict_.lookup("printEvery")); + if (propsDict_.found("sampleAll")) sampleAll_ = true; + if (propsDict_.found("probeDebug")) probeDebug_ = true; + if (propsDict_.found("includePosition")) includePosition_ = true; - if (propsDict_.found("writePrecision")) writePrecision_= readScalar(propsDict_.lookup("writePrecision")); + if (propsDict_.found("writePrecision")) writePrecision_ = readScalar(propsDict_.lookup("writePrecision")); } @@ -100,25 +100,36 @@ particleProbe::particleProbe particleProbe::~particleProbe() { - clearProbes(); + forAll(sPtrList_, i) + delete sPtrList_[i]; } // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // -void particleProbe::setOutputFile() const +void particleProbe::setOutputFile(const word& logFileName) { - //set the current item ID - if(currItemId_== itemCounter_) - currItemId_=1; - else - currItemId_+=1; - sPtr = sPtrList_[currItemId_-1]; //set the pointer to the output file from list - probeIndex_=currItemId_-1; + if (itemCounter_ > 0 && verboseToFile_) + { + bool foundFile = false; + forAll(itemsToSample_, i) + { + if (itemsToSample_[i] == logFileName) + { + probeIndex_ = i; + foundFile = true; + } + } + + if(!foundFile) + FatalError << "particleProbe::setOutputFile for logFileName " << logFileName << " : " << "File not found" << abort(FatalError); + currItemId_ = probeIndex_ + 1; + setCounter(); + } } -void particleProbe::initialize(word typeName, word logFileName) const +void particleProbe::initialize(const word& modelName, const word& logFileName) { //update the list of items to be sampled ++itemCounter_; @@ -126,7 +137,7 @@ void particleProbe::initialize(word typeName, word logFileName) const // init environment //propsDict_ = particleCloud_.couplingProperties().subDict(typeName + "Props"); - name_ = typeName; + name_ = modelName; if (verboseToFile_) { @@ -136,8 +147,8 @@ void particleProbe::initialize(word typeName, word logFileName) const MPI_Comm_rank(MPI_COMM_WORLD, &rank_); //open a separate file for each processor - char* filecurrent_ = new char[strlen(logFileName.c_str()) + 4]; //reserve 4 chars for processor name - sprintf(filecurrent_,"%s%s%d", logFileName.c_str(), ".", rank_); + char* filecurrent_ = new char[logFileName.length() + 1 + 4 + 1]; //reserve 4 chars for processor name + sprintf(filecurrent_,"%s.%d", logFileName.c_str(), rank_); Info << "particleProbe for model " << name_ << " will write to file " << filecurrent_ << endl; @@ -171,120 +182,44 @@ void particleProbe::initialize(word typeName, word logFileName) const scalarFields_.clear(); vectorFields_.clear(); } - - return; } + void particleProbe::writeHeader() const { + if (verboseToFile_) + { + *sPtr << "#processor: " << rank_ << endl; + *sPtr << "#index time " << " "; + *sPtr << "|| vectorData: " << " "; - if(verboseToFile_ ) - { - *sPtr<<"#processor: " << rank_ << endl; - *sPtr<<"#index time " << " "; + forAll(vectorFields_, iter) + { + if (!probeDebug_ && iter > 0) break; + *sPtr << vectorFields_(iter) << " "; + } - - *sPtr<<"|| vectorData: " << " "; - - forAll(vectorFields_, iter) - { - if(!probeDebug_ && iter>0) break; - *sPtr << vectorFields_(iter) << " "; - } - - if(probeDebug_) - { - *sPtr<<"|| scalarData: " << " "; + if (probeDebug_) + { + *sPtr << "|| scalarData: " << " "; forAll(scalarFields_, iter) { - *sPtr << scalarFields_(iter) << " "; + *sPtr << scalarFields_(iter) << " "; } - } - - if(includePosition_) *sPtr<<" || position" << endl; - else *sPtr << endl; - } - -} - -void particleProbe::clearProbes() const -{ - for (unsigned int i=0; i sValues, Field vValues) const -{ - int vSize_=vProbes_.size(); - int sSize_=sProbes_.size(); - - //check if the particle already has an allocated vector. If not, create it. It should be only called at the beginning. - while(index >= vSize_) - { - std::vector particleVector_; - vProbes_.push_back(particleVector_); - vSize_=vProbes_.size(); - } - - while(index >= sSize_) - { - std::vector particleScalar_; - sProbes_.push_back(particleScalar_); - sSize_=sProbes_.size(); - } - - //register vector probes on the corresponding vector - forAll(vValues, iter) - { - int ProbeSize_=vProbes_[index].size(); - - if(probeIndex_ sValues, Field vValues) const + +void particleProbe::writeProbe(int index, Field sValues, Field vValues) { - updateProbes(index,sValues,vValues); //update probe vectors - - - if(printNow_ && checkIDForPrint(index) && verboseToFile_) + if (printNow_ && verboseToFile_ && checkIDForPrint(index)) { + sPtr = sPtrList_[probeIndex_]; //set the pointer to the output file from list + //index and time *sPtr << setprecision(IOstream::defaultPrecision()+7); *sPtr << index << tab << particleCloud_.mesh().time().value() << " "; @@ -298,7 +233,6 @@ void particleProbe::writeProbe(int index, Field sValues, Field v *sPtr << vValues[iter][0] << " "; *sPtr << vValues[iter][1] << " "; *sPtr << vValues[iter][2] << " "; - } //scalarFields @@ -324,42 +258,40 @@ void particleProbe::writeProbe(int index, Field sValues, Field v *sPtr << endl; } } - - return; } + bool particleProbe::checkIDForPrint(int index) const { - bool sampleThisId_ = false; - if(sampleAll_) sampleThisId_ = true; - else - { - forAll(particleIDsToSample_, iSample) - { - if(index==particleIDsToSample_[iSample]) sampleThisId_ = true; - } - } - return sampleThisId_; + if(sampleAll_) + { + return true; + } + else + { + forAll(particleIDsToSample_, iSample) + { + if (index == particleIDsToSample_[iSample]) return true; + } + } + return false; } -void particleProbe::setCounter() const +void particleProbe::setCounter() { - - //reset or increment counter for printing to file + //reset or increment counter for printing to file //Do only if called by first item in the list of items! - if(currItemId_==1) + if (currItemId_ == 1) { - printCounter_++; - if( printCounter_ >= printEvery_ ) + ++printCounter_; + + if (printCounter_ >= printEvery_) { - printCounter_=0; + printCounter_ = 0; printNow_ = true; } else printNow_ = false; - } - return; - } // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // diff --git a/src/lagrangian/cfdemParticle/subModels/probeModel/particleProbe/particleProbe.H b/src/lagrangian/cfdemParticle/subModels/probeModel/particleProbe/particleProbe.H index 3155a16e..e840f790 100644 --- a/src/lagrangian/cfdemParticle/subModels/probeModel/particleProbe/particleProbe.H +++ b/src/lagrangian/cfdemParticle/subModels/probeModel/particleProbe/particleProbe.H @@ -63,7 +63,7 @@ private: dictionary propsDict_; - mutable word name_; + word name_; cfdemCloud& particleCloud_; @@ -75,7 +75,7 @@ private: word dirName_; - mutable int rank_; + int rank_; mutable OFstream* sPtr; @@ -89,25 +89,23 @@ private: const labelList particleIDsToSample_; - mutable wordList itemsToSample_; + wordList itemsToSample_; - mutable List sPtrList_; + List sPtrList_; - mutable int itemCounter_; + int itemCounter_; - mutable int currItemId_; + int currItemId_; - mutable int printCounter_; + int printCounter_; - mutable bool printNow_; + bool printNow_; - mutable std::vector< std::vector > sProbes_; + std::vector probeName_; - mutable std::vector< std::vector > vProbes_; + int probeIndex_; - mutable std::vector probeName_; - - mutable int probeIndex_; + void setCounter(); public: @@ -131,17 +129,11 @@ public: ~particleProbe(); // Member Functions - void updateProbes(int index, Field sValues, Field vValues) const; - void initialize(word typeName, word logFileName) const; - void setOutputFile() const; + void initialize(const word& modelName, const word& logFileName); + void setOutputFile(const word& logFileName); void writeHeader() const; - void writeProbe(int index, Field sValues, Field vValues) const; + void writeProbe(int index, Field sValues, Field vValues); bool checkIDForPrint(int) const; - void setCounter() const; - void clearProbes() const; - std::vector< std::vector >* getVprobe() { return &vProbes_; } - std::vector< std::vector >* getSprobe() { return &sProbes_; } - }; diff --git a/src/lagrangian/cfdemParticle/subModels/probeModel/probeModel/probeModel.H b/src/lagrangian/cfdemParticle/subModels/probeModel/probeModel/probeModel.H index d30bd808..a7c811f0 100644 --- a/src/lagrangian/cfdemParticle/subModels/probeModel/probeModel/probeModel.H +++ b/src/lagrangian/cfdemParticle/subModels/probeModel/probeModel/probeModel.H @@ -133,16 +133,12 @@ public: // Member Functions - virtual void initialize(word typeName, word logFileName) const {} - virtual void setOutputFile() const {} + virtual void initialize(const word& modelName, const word& logFileName) {} + virtual void setOutputFile(const word& logFileName) {} virtual void writeHeader() const {} - virtual void writeProbe(int index, Field sValues, Field vValues) const {} + virtual void writeProbe(int index, Field sValues, Field vValues) {} virtual bool checkIDForPrint(int) const { return false; } - virtual void setCounter() const {} virtual bool active() const { return true; } - virtual std::vector< std::vector >* getVprobe() { return NULL; } - virtual std::vector< std::vector >* getSprobe() { return NULL; } - virtual void clearProbes() const {} // Access diff --git a/src/lagrangian/cfdemParticle/subModels/regionModel/regionModel/regionModel.C b/src/lagrangian/cfdemParticle/subModels/regionModel/regionModel/regionModel.C index d9a6ed6d..b1770d07 100644 --- a/src/lagrangian/cfdemParticle/subModels/regionModel/regionModel/regionModel.C +++ b/src/lagrangian/cfdemParticle/subModels/regionModel/regionModel/regionModel.C @@ -46,13 +46,13 @@ defineRunTimeSelectionTable(regionModel, dictionary); // * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // -void Foam::regionModel::reAllocArrays() const +void regionModel::reAllocArrays() const { if(particleCloud_.numberOfParticlesChanged()) { // get arrays of new length - particleCloud_.dataExchangeM().allocateArray(inRegion_,1,1); - particleCloud_.dataExchangeM().allocateArray(outRegion_,1,1); + particleCloud_.dataExchangeM().allocateArray(inRegion_,1.,1); + particleCloud_.dataExchangeM().allocateArray(outRegion_,1.,1); } } @@ -83,8 +83,8 @@ regionModel::regionModel regionModel::~regionModel() { - free(inRegion_); - free(outRegion_); + particleCloud_.dataExchangeM().destroy(inRegion_,1); + particleCloud_.dataExchangeM().destroy(outRegion_,1); } diff --git a/src/lagrangian/cfdemParticle/subModels/smoothingModel/constDiffSmoothing/constDiffSmoothing.C b/src/lagrangian/cfdemParticle/subModels/smoothingModel/constDiffSmoothing/constDiffSmoothing.C index 1af57abe..dc739cc2 100644 --- a/src/lagrangian/cfdemParticle/subModels/smoothingModel/constDiffSmoothing/constDiffSmoothing.C +++ b/src/lagrangian/cfdemParticle/subModels/smoothingModel/constDiffSmoothing/constDiffSmoothing.C @@ -6,7 +6,7 @@ Christoph Goniva, christoph.goniva@cfdem.com Copyright 2009-2012 JKU Linz Copyright 2012- DCS Computing GmbH, Linz - Copyright (C) 2013- Graz University of + Copyright (C) 2013- Graz University of Technology, IPPT ------------------------------------------------------------------------------- License @@ -71,10 +71,10 @@ constDiffSmoothing::constDiffSmoothing verbose_(false) { - if(propsDict_.found("verbose")) + if(propsDict_.found("verbose")) verbose_ = true; - if(propsDict_.found("smoothingLengthReferenceField")) + if(propsDict_.found("smoothingLengthReferenceField")) smoothingLengthReferenceField_.value() = double(readScalar(propsDict_.lookup("smoothingLengthReferenceField"))); checkFields(sSmoothField_); @@ -94,7 +94,7 @@ bool constDiffSmoothing::doSmoothing() const } -void Foam::constDiffSmoothing::smoothen(volScalarField& fieldSrc) const +void constDiffSmoothing::smoothen(volScalarField& fieldSrc) const { // Create scalar smooth field from virgin scalar smooth field template volScalarField sSmoothField = sSmoothField_; @@ -120,11 +120,11 @@ void Foam::constDiffSmoothing::smoothen(volScalarField& fieldSrc) const forAll(sSmoothField,cellI) { sSmoothField[cellI]=max(lowerLimit_,min(upperLimit_,sSmoothField[cellI])); - } + } // get data from working sSmoothField - will copy only values at new time fieldSrc=sSmoothField; - fieldSrc.correctBoundaryConditions(); + fieldSrc.correctBoundaryConditions(); if(verbose_) { @@ -135,7 +135,7 @@ void Foam::constDiffSmoothing::smoothen(volScalarField& fieldSrc) const } // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // -void Foam::constDiffSmoothing::smoothen(volVectorField& fieldSrc) const +void constDiffSmoothing::smoothen(volVectorField& fieldSrc) const { // Create scalar smooth field from virgin scalar smooth field template volVectorField vSmoothField = vSmoothField_; @@ -159,7 +159,7 @@ void Foam::constDiffSmoothing::smoothen(volVectorField& fieldSrc) const // get data from working vSmoothField fieldSrc=vSmoothField; - fieldSrc.correctBoundaryConditions(); + fieldSrc.correctBoundaryConditions(); if(verbose_) { @@ -170,7 +170,7 @@ void Foam::constDiffSmoothing::smoothen(volVectorField& fieldSrc) const } // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // -void Foam::constDiffSmoothing::smoothenReferenceField(volVectorField& fieldSrc) const +void constDiffSmoothing::smoothenReferenceField(volVectorField& fieldSrc) const { // Create scalar smooth field from virgin scalar smooth field template volVectorField vSmoothField = vSmoothField_; @@ -185,7 +185,7 @@ void Foam::constDiffSmoothing::smoothenReferenceField(volVectorField& fieldSrc) double sourceStrength = 1e5; //large number to keep reference values constant dimensionedScalar deltaT = vSmoothField.mesh().time().deltaT(); - DT_.value() = smoothingLengthReferenceField_.value() + DT_.value() = smoothingLengthReferenceField_.value() * smoothingLengthReferenceField_.value() / deltaT.value(); tmp NLarge @@ -194,7 +194,7 @@ void Foam::constDiffSmoothing::smoothenReferenceField(volVectorField& fieldSrc) ( IOobject ( - "xxx", + "NLarge", particleCloud_.mesh().time().timeName(), particleCloud_.mesh(), IOobject::NO_READ, @@ -218,14 +218,14 @@ void Foam::constDiffSmoothing::smoothenReferenceField(volVectorField& fieldSrc) ( fvm::ddt(vSmoothField) -fvm::laplacian( DT_, vSmoothField) - == + == NLarge() / deltaT * vSmoothField.oldTime() //add source to keep cell values constant -fvm::Sp( NLarge() / deltaT, vSmoothField) //add sink to keep cell values constant ); // get data from working vSmoothField fieldSrc=vSmoothField; - fieldSrc.correctBoundaryConditions(); + fieldSrc.correctBoundaryConditions(); if(verbose_) { diff --git a/src/lagrangian/cfdemParticle/subModels/smoothingModel/noSmoothing/noSmoothing.C b/src/lagrangian/cfdemParticle/subModels/smoothingModel/noSmoothing/noSmoothing.C index 01a73ab1..61f5bcfd 100644 --- a/src/lagrangian/cfdemParticle/subModels/smoothingModel/noSmoothing/noSmoothing.C +++ b/src/lagrangian/cfdemParticle/subModels/smoothingModel/noSmoothing/noSmoothing.C @@ -70,13 +70,13 @@ noSmoothing::~noSmoothing() // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // -void Foam::noSmoothing::smoothen(volScalarField& field) const +void noSmoothing::smoothen(volScalarField& field) const {} -void Foam::noSmoothing::smoothen(volVectorField& field) const +void noSmoothing::smoothen(volVectorField& field) const {} -void Foam::noSmoothing::smoothenReferenceField(volVectorField& field) const +void noSmoothing::smoothenReferenceField(volVectorField& field) const {} // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // diff --git a/src/lagrangian/cfdemParticle/subModels/smoothingModel/smoothingModel/smoothingModel.C b/src/lagrangian/cfdemParticle/subModels/smoothingModel/smoothingModel/smoothingModel.C index b906a160..a98d4bde 100644 --- a/src/lagrangian/cfdemParticle/subModels/smoothingModel/smoothingModel/smoothingModel.C +++ b/src/lagrangian/cfdemParticle/subModels/smoothingModel/smoothingModel/smoothingModel.C @@ -50,7 +50,7 @@ defineRunTimeSelectionTable(smoothingModel, dictionary); // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // // Construct from components -Foam::smoothingModel::smoothingModel +smoothingModel::smoothingModel ( const dictionary& dict, cfdemCloud& sm @@ -59,7 +59,7 @@ Foam::smoothingModel::smoothingModel dict_(dict), particleCloud_(sm), vSmoothField_ - ( + ( IOobject ( "vSmoothField", @@ -72,7 +72,7 @@ Foam::smoothingModel::smoothingModel dimensionedVector("zero", dimensionSet(0,0,0,0,0), vector::zero) ), sSmoothField_ - ( + ( IOobject ( "sSmoothField", @@ -89,7 +89,7 @@ Foam::smoothingModel::smoothingModel // * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * // -Foam::smoothingModel::~smoothingModel() +smoothingModel::~smoothingModel() {} @@ -110,7 +110,7 @@ void smoothingModel::checkFields(volVectorField& vSmoothField_) const { // currently it is detected if field was auto generated or defined // improvement would be changing the type here automatically - forAll(vSmoothField_.boundaryField(),patchI) + forAll(vSmoothField_.boundaryField(),patchI) if(vSmoothField_.boundaryField()[patchI].type()=="calculated") FatalError <<"Vector field:"<< vSmoothField_.name() << " must be defined.\n" << abort(FatalError); diff --git a/src/lagrangian/cfdemParticle/subModels/voidFractionModel/GaussVoidFraction/GaussVoidFraction.C b/src/lagrangian/cfdemParticle/subModels/voidFractionModel/GaussVoidFraction/GaussVoidFraction.C index 9a929d67..231e7e97 100644 --- a/src/lagrangian/cfdemParticle/subModels/voidFractionModel/GaussVoidFraction/GaussVoidFraction.C +++ b/src/lagrangian/cfdemParticle/subModels/voidFractionModel/GaussVoidFraction/GaussVoidFraction.C @@ -137,7 +137,7 @@ void GaussVoidFraction::setvoidFraction(double** const& mask,double**& voidfract if (hashSetLength > maxCellsPerParticle_) { FatalError<< "big particle algo found more cells ("<< hashSetLength - <<") than storage is prepered ("< 0) { diff --git a/src/lagrangian/cfdemParticle/subModels/voidFractionModel/IBVoidFraction/IBVoidFraction.C b/src/lagrangian/cfdemParticle/subModels/voidFractionModel/IBVoidFraction/IBVoidFraction.C index 2b60cfac..5e016a19 100755 --- a/src/lagrangian/cfdemParticle/subModels/voidFractionModel/IBVoidFraction/IBVoidFraction.C +++ b/src/lagrangian/cfdemParticle/subModels/voidFractionModel/IBVoidFraction/IBVoidFraction.C @@ -36,8 +36,6 @@ Description #include "locateModel.H" #include "dataExchangeModel.H" -#include - // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // namespace Foam @@ -68,17 +66,15 @@ IBVoidFraction::IBVoidFraction propsDict_(dict.subDict(typeName + "Props")), alphaMin_(readScalar(propsDict_.lookup("alphaMin"))), alphaLimited_(0), - scaleUpVol_(readScalar(propsDict_.lookup("scaleUpVol"))), - checkPeriodicCells_(false) + scaleUpVol_(readScalar(propsDict_.lookup("scaleUpVol"))) { - Info << "\n\n W A R N I N G - do not use in combination with differentialRegion model! \n\n" << endl; - //Info << "\n\n W A R N I N G - this model does not yet work properly! \n\n" << endl; - maxCellsPerParticle_=readLabel(propsDict_.lookup("maxCellsPerParticle")); + Info << "\n\n W A R N I N G - do not use in combination with differentialRegion model!\n\n" << endl; + maxCellsPerParticle_ = readLabel(propsDict_.lookup("maxCellsPerParticle")); - if(scaleUpVol_ < 1){ FatalError<< "scaleUpVol shloud be > 1."<< abort(FatalError); } - if(alphaMin_ > 1 || alphaMin_ < 0.01){ FatalError<< "alphaMin shloud be > 1 and < 0.01." << abort(FatalError); } - - if(propsDict_.found("checkPeriodicCells")) checkPeriodicCells_=true; + if (scaleUpVol_ < 1.0) + FatalError << "scaleUpVol shloud be > 1." << abort(FatalError); + if (alphaMin_ > 1.0 || alphaMin_ < 0.01) + FatalError << "alphaMin shloud be > 1 and < 0.01." << abort(FatalError); } @@ -96,25 +92,26 @@ void IBVoidFraction::setvoidFraction(double** const& mask,double**& voidfraction reAllocArrays(); - voidfractionNext_.ref()=1; + voidfractionNext_.ref() = 1.0; - for(int index=0; index< particleCloud_.numberOfParticles(); index++) + for (int index=0; index < particleCloud_.numberOfParticles(); index++) { //if(mask[index][0]) //{ //reset - for(int subcell=0;subcell= 0) { @@ -122,92 +119,66 @@ void IBVoidFraction::setvoidFraction(double** const& mask,double**& voidfraction //compute the voidfraction for the cell "particleCentreCellID vector cellCentrePosition = particleCloud_.mesh().C()[particleCenterCellID]; - scalar centreDist=mag(cellCentrePosition-positionCenter); + scalar fc = pointInParticle(index, positionCenter, cellCentrePosition); + vector minPeriodicParticlePos = positionCenter; - vector minPeriodicParticlePos=positionCenter; - if(checkPeriodicCells_) //consider minimal distance to all periodic images of this particle + if (particleCloud_.checkPeriodicCells()) //consider minimal distance to all periodic images of this particle { - centreDist = minPeriodicDistance(cellCentrePosition, positionCenter, globalBb, - minPeriodicParticlePos); + fc = minPeriodicDistance(index,cellCentrePosition, positionCenter, globalBb, minPeriodicParticlePos); } - if(centreDist + 0.5*sqrt(3.0)*pow(particleCloud_.mesh().V()[particleCenterCellID],0.33333) < radius) + scalar centreDist = mag(cellCentrePosition - minPeriodicParticlePos); + scalar corona = 0.5 * sqrt(3.0) * pow(particleCloud_.mesh().V()[particleCenterCellID], 1./3.); + vector coronaPoint = cellCentrePosition + (cellCentrePosition - minPeriodicParticlePos) * (corona / centreDist); + + if (pointInParticle(index, minPeriodicParticlePos, coronaPoint) < 0.0) { - voidfractionNext_[particleCenterCellID] = 0; + voidfractionNext_[particleCenterCellID] = 0.0; } else { const labelList& vertices = particleCloud_.mesh().cellPoints()[particleCenterCellID]; + const double ratio = 0.125; + forAll(vertices, i) { vector vertexPosition = particleCloud_.mesh().points()[vertices[i]]; - scalar centreVertexDist = mag(vertexPosition-positionCenter); - if(checkPeriodicCells_) //consider minimal distance to all periodic images of this particle + scalar fv = pointInParticle(index, positionCenter, vertexPosition); + + if (particleCloud_.checkPeriodicCells()) //consider minimal distance to all periodic images of this particle { - centreVertexDist = minPeriodicDistance(vertexPosition, positionCenter, globalBb, - minPeriodicParticlePos); + fv = minPeriodicDistance(index, vertexPosition, positionCenter, globalBb, minPeriodicParticlePos); } - if(centreDistradius) + else if (fc < 0.0 && fv > 0.0) { //compute lambda - scalar a = (vertexPosition - cellCentrePosition) - & (vertexPosition - cellCentrePosition); - scalar b = 2. * (vertexPosition - cellCentrePosition) - & (cellCentrePosition-minPeriodicParticlePos); - scalar c = ( (cellCentrePosition-minPeriodicParticlePos) - &(cellCentrePosition-minPeriodicParticlePos) - ) - - radius*radius; - - scalar lambda = 0.; - - if (b*b-4.*a*c>=0.) lambda = (-b+sqrt(b*b-4.*a*c))/(2.*a); - if (lambda > 0 && lambda <=1) voidfractionNext_[particleCenterCellID]-=lambda*0.125; - else - { - lambda = (-b-sqrt(b*b-4.*a*c))/(2.*a); - if (lambda > 0. && lambda <=1.) voidfractionNext_[particleCenterCellID]-=lambda*0.125; - } + scalar lambda = segmentParticleIntersection(index, minPeriodicParticlePos, cellCentrePosition, vertexPosition); + voidfractionNext_[particleCenterCellID] -= ratio * lambda; } - else if(centreDist>radius && centreVertexDist 0.0 && fv < 0.0) { - //compute another lambda too - scalar a = (vertexPosition - cellCentrePosition) - & (vertexPosition - cellCentrePosition); - scalar b = 2.* (vertexPosition - cellCentrePosition) - & (cellCentrePosition-minPeriodicParticlePos); - scalar c = ( (cellCentrePosition-minPeriodicParticlePos) - &(cellCentrePosition-minPeriodicParticlePos) - ) - - radius*radius; - scalar lambda = 0.; - - if(b*b-4.*a*c>=0.) lambda = (-b+sqrt(b*b-4.*a*c))/(2.*a); - if(lambda > 0. && lambda <=1.) voidfractionNext_[particleCenterCellID]-=(1.-lambda)*0.125; - else - { - lambda = (-b-sqrt(b*b-4.*a*c))/(2.*a); - if (lambda > 0. && lambda <=1.) voidfractionNext_[particleCenterCellID]-=(1.-lambda)*0.125; - } + //compute lambda + scalar lambda = segmentParticleIntersection(index, minPeriodicParticlePos, vertexPosition, cellCentrePosition); + voidfractionNext_[particleCenterCellID] -= ratio * lambda; } } } //end particle partially overlapping with cell //generating list with cell and subcells - buildLabelHashSet(radius, minPeriodicParticlePos, particleCenterCellID, hashSett, true); + buildLabelHashSet(index, minPeriodicParticlePos, particleCenterCellID, hashSett, true); //Add cells of periodic particle images on same processor - if(checkPeriodicCells_) + if (particleCloud_.checkPeriodicCells()) { int doPeriodicImage[3]; - for (int iDir=0; iDir<3; iDir++) + for (int iDir=0; iDir < 3; iDir++) { - doPeriodicImage[iDir]= 0; + doPeriodicImage[iDir] = 0; if ((minPeriodicParticlePos[iDir]+radius) > globalBb.max()[iDir]) { doPeriodicImage[iDir] = -1; @@ -238,10 +209,10 @@ void IBVoidFraction::setvoidFraction(double** const& mask,double**& voidfraction copyCounter++; } //y-direction - int currCopyCounter=copyCounter; + int currCopyCounter = copyCounter; if (doPeriodicImage[1] != 0) { - for(int yDirCop=0; yDirCop<=currCopyCounter; yDirCop++) + for (int yDirCop=0; yDirCop <= currCopyCounter; yDirCop++) { particlePosList.append( particlePosList[yDirCop] + vector( @@ -254,10 +225,10 @@ void IBVoidFraction::setvoidFraction(double** const& mask,double**& voidfraction } } //z-direction - currCopyCounter=copyCounter; + currCopyCounter = copyCounter; if (doPeriodicImage[2] != 0) { - for(int zDirCop=0; zDirCop<=currCopyCounter; zDirCop++) + for (int zDirCop=0; zDirCop <= currCopyCounter; zDirCop++) { particlePosList.append( particlePosList[zDirCop] + vector( @@ -273,7 +244,7 @@ void IBVoidFraction::setvoidFraction(double** const& mask,double**& voidfraction //add the nearest cell labels particleLabelList.append(particleCenterCellID); - for(int iPeriodicImage=1; iPeriodicImage<=copyCounter; iPeriodicImage++) + for (int iPeriodicImage=1; iPeriodicImage <= copyCounter; iPeriodicImage++) { label partCellId = particleCloud_.mesh().findNearestCell(particlePosList[iPeriodicImage]); particleLabelList.append(partCellId); @@ -287,17 +258,17 @@ void IBVoidFraction::setvoidFraction(double** const& mask,double**& voidfraction label hashSetLength = hashSett.size(); if (hashSetLength > maxCellsPerParticle_) { - FatalError<< "big particle algo found more cells ("<< hashSetLength - <<") than storage is prepared ("< 0) { cellsPerParticle_[index][0]=hashSetLength; hashSett.erase(particleCenterCellID); - for(label i=0;i= 0) { + if (voidfractionNext_[cellID] < 0.0) + voidfractionNext_[cellID] = 0.0; voidfractions[index][subcell] = voidfractionNext_[cellID]; } else @@ -325,14 +298,15 @@ void IBVoidFraction::setvoidFraction(double** const& mask,double**& voidfraction void IBVoidFraction::buildLabelHashSet ( - const scalar radius, + int index, const vector position, const label cellID, labelHashSet& hashSett, bool initialInsert //initial insertion of own cell )const { - if(initialInsert) hashSett.insert(cellID); + if(initialInsert) + hashSett.insert(cellID); const labelList& nc = particleCloud_.mesh().cellCells()[cellID]; forAll(nc,i) @@ -341,115 +315,99 @@ void IBVoidFraction::buildLabelHashSet vector cellCentrePosition = particleCloud_.mesh().C()[neighbor]; scalar centreDist = mag(cellCentrePosition-position); - if(!hashSett.found(neighbor) && centreDist + 0.5*sqrt(3.0)*pow(particleCloud_.mesh().V()[neighbor],0.33333) < radius) + scalar fc = pointInParticle(index, position, cellCentrePosition); + scalar corona = 0.5 * sqrt(3.0) * pow(particleCloud_.mesh().V()[neighbor], 1./3.); + vector coronaPoint = cellCentrePosition + (cellCentrePosition - position) * (corona / centreDist); + + if (!hashSett.found(neighbor) && pointInParticle(index, position, coronaPoint) < 0.0) { voidfractionNext_[neighbor] = 0.; - buildLabelHashSet(radius,position,neighbor,hashSett,true); + buildLabelHashSet(index,position,neighbor,hashSett,true); } - else if(!hashSett.found(neighbor) && centreDist < radius + sqrt(3.0)*pow(particleCloud_.mesh().V()[neighbor],0.33333)) + else if(!hashSett.found(neighbor)) { scalar scale = 1.; const labelList& vertexPoints = particleCloud_.mesh().cellPoints()[neighbor]; + const scalar ratio = 0.125; forAll(vertexPoints, j) { vector vertexPosition = particleCloud_.mesh().points()[vertexPoints[j]]; - scalar vertexDist = mag(vertexPosition - position); + scalar fv = pointInParticle(index, position, vertexPosition); - if (centreDist < radius) + if (fc < 0.0) { - if (vertexDist < radius) + if (fv < 0.0) { - scale -= 0.125; + scale -= ratio; } else { - scalar a = (vertexPosition - cellCentrePosition)&(vertexPosition - cellCentrePosition); - scalar b = 2.* (vertexPosition - cellCentrePosition)&(cellCentrePosition-position); - scalar c = ((cellCentrePosition-position)&(cellCentrePosition-position))-radius*radius; - scalar lambda = 0.; - - if(b*b-4.*a*c >= 0.) lambda = (-b+sqrt(b*b-4.*a*c))/(2.*a); - if (lambda > 0. && lambda <=1.) - { - scale -=lambda * 0.125; - } - else - { - lambda = (-b-sqrt(b*b-4.*a*c))/(2.*a); - if (lambda > 0. && lambda <= 1.) scale -=lambda * 0.125; - } + scalar lambda = segmentParticleIntersection(index, position, cellCentrePosition, vertexPosition); + scale -= lambda * ratio; } } - else if (vertexDist < radius) + else if (fv < 0.0) { - scalar a = (vertexPosition - cellCentrePosition)&(vertexPosition - cellCentrePosition); - scalar b = 2.* (vertexPosition - cellCentrePosition)&(cellCentrePosition-position); - scalar c = ((cellCentrePosition-position)&(cellCentrePosition-position))-radius*radius; - scalar lambda = 0.; - - if (b*b-4.*a*c >= 0.) lambda = (-b+sqrt(b*b-4.*a*c))/(2.*a); - if (lambda > 0. && lambda <= 1.) - { - scale -=(1.-lambda) * 0.125; - } - else - { - lambda = (-b-sqrt(b*b-4.*a*c))/(2.*a); - if (lambda > 0. && lambda <= 1.) scale -= (1.-lambda) * 0.125; - } + scalar lambda = segmentParticleIntersection(index, position, vertexPosition, cellCentrePosition); + scale -= lambda * ratio; } } - if(voidfractionNext_[neighbor] == 1.0) + if (scale < 0.0) + scale = 0.0; + + if (voidfractionNext_[neighbor] == 1.0) { voidfractionNext_[neighbor] = scale; } else { voidfractionNext_[neighbor] -= (1.0-scale); - if(voidfractionNext_[neighbor] < 0.) voidfractionNext_[neighbor] = 0.0; + if (voidfractionNext_[neighbor] < 0.) + voidfractionNext_[neighbor] = 0.0; } - if(!(scale == 1.0)) buildLabelHashSet(radius,position,neighbor,hashSett, true); + + if (!(scale == 1.0)) + buildLabelHashSet(index,position,neighbor,hashSett, true); } } } -//Function to determine minimal distance of point -//to one of the periodic images of a particle -double IBVoidFraction::minPeriodicDistance(vector cellCentrePosition, - vector positionCenter, - boundBox globalBb, - vector& minPeriodicPos)const +double IBVoidFraction::segmentParticleIntersection(int index, vector positionCenter, vector pointInside, vector pointOutside) const { - double centreDist = 999e32; - vector positionCenterPeriodic; + const scalar radius = particleCloud_.radius(index); + const scalar a = (pointOutside - pointInside) & (pointOutside - pointInside); + const scalar b = 2.*(pointOutside - pointInside) & (pointInside - positionCenter); + const scalar c = ((pointInside - positionCenter) & (pointInside - positionCenter)) - radius*radius; + const scalar D = b*b - 4.0*a*c; + const scalar eps = 1e-12; + scalar lambda_ = 0.0; + scalar lambda = 0.0; - for (int xDir=-1; xDir<=1; xDir++) + if (D >= 0.0) { - positionCenterPeriodic[0] = positionCenter[0] - + static_cast(xDir) - * (globalBb.max()[0] - globalBb.min()[0]); - for (int yDir=-1; yDir<=1; yDir++) + scalar sqrtD = sqrt(D); + lambda_ = (-b + sqrtD)/(2.0*a); + + if (lambda_ >= -eps && lambda_ <= 1.0+eps) { - positionCenterPeriodic[1] = positionCenter[1] - + static_cast(yDir) - * (globalBb.max()[1] - globalBb.min()[1]); - for (int zDir=-1; zDir<=1; zDir++) - { - positionCenterPeriodic[2] = positionCenter[2] - + static_cast(zDir) - * (globalBb.max()[2] - globalBb.min()[2]); - if (mag(cellCentrePosition-positionCenterPeriodic) < centreDist) - { - centreDist = mag(cellCentrePosition-positionCenterPeriodic); - minPeriodicPos = positionCenterPeriodic; - } - } + lambda = lambda_; + } + else + { + lambda_ = (-b - sqrtD)/(2.0*a); + if (lambda_ >= -eps && lambda_ <= 1.0+eps) + lambda = lambda_; } } - return centreDist; + if (lambda < 0.0) + return 0.0; + else if (lambda > 1.0) + return 1.0; + + return lambda; } // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // diff --git a/src/lagrangian/cfdemParticle/subModels/voidFractionModel/IBVoidFraction/IBVoidFraction.H b/src/lagrangian/cfdemParticle/subModels/voidFractionModel/IBVoidFraction/IBVoidFraction.H index 78620719..e2072025 100755 --- a/src/lagrangian/cfdemParticle/subModels/voidFractionModel/IBVoidFraction/IBVoidFraction.H +++ b/src/lagrangian/cfdemParticle/subModels/voidFractionModel/IBVoidFraction/IBVoidFraction.H @@ -28,7 +28,7 @@ Description This code is designed to realize coupled CFD-DEM simulations using LIGGGHTS and OpenFOAM(R). Note: this code is not part of OpenFOAM(R) (see DISCLAIMER). - void fraction model for the smooth representation of spheres with + void fraction model for the smooth representation of spheres with radius > cell length. contribution from Alice Hager @@ -64,11 +64,9 @@ private: const scalar alphaMin_; //NP min value of voidFraction - mutable bool alphaLimited_; + mutable bool alphaLimited_; const scalar scaleUpVol_; //NP scaling radius, keeping volume of particle - - mutable bool checkPeriodicCells_; public: @@ -95,22 +93,14 @@ public: void buildLabelHashSet ( - const scalar radius, + int index, const vector position, const label cellID, labelHashSet& hashSett, bool initialInsert ) const; - - double minPeriodicDistance - ( - vector cellCentrePosition, - vector positionCenter, - boundBox globalBb, - vector& minPeriodicPos - ) const; - - + + virtual double segmentParticleIntersection(int index, vector positionCenter, vector pointInside, vector pointOutside) const; }; diff --git a/src/lagrangian/cfdemParticle/subModels/voidFractionModel/bigParticleVoidFraction/bigParticleVoidFraction.C b/src/lagrangian/cfdemParticle/subModels/voidFractionModel/bigParticleVoidFraction/bigParticleVoidFraction.C index 179a3b6a..63405dca 100644 --- a/src/lagrangian/cfdemParticle/subModels/voidFractionModel/bigParticleVoidFraction/bigParticleVoidFraction.C +++ b/src/lagrangian/cfdemParticle/subModels/voidFractionModel/bigParticleVoidFraction/bigParticleVoidFraction.C @@ -132,7 +132,7 @@ void bigParticleVoidFraction::setvoidFraction(double** const& mask,double**& voi if (hashSetLength > maxCellsPerParticle_) { FatalError<< "big particle algo found more cells ("<< hashSetLength - <<") than storage is prepered ("< 0) { diff --git a/src/lagrangian/cfdemParticle/subModels/voidFractionModel/centreVoidFraction/centreVoidFraction.C b/src/lagrangian/cfdemParticle/subModels/voidFractionModel/centreVoidFraction/centreVoidFraction.C index f04cfc07..a9237811 100644 --- a/src/lagrangian/cfdemParticle/subModels/voidFractionModel/centreVoidFraction/centreVoidFraction.C +++ b/src/lagrangian/cfdemParticle/subModels/voidFractionModel/centreVoidFraction/centreVoidFraction.C @@ -140,6 +140,7 @@ void centreVoidFraction::setvoidFraction(double** const& mask,double**& voidfrac } //} } + voidfractionNext_.correctBoundaryConditions(); // bring voidfraction from Eulerian Field to particle array for(int index=0; index< particleCloud_.numberOfParticles(); index++) diff --git a/src/lagrangian/cfdemParticle/subModels/voidFractionModel/dividedVoidFraction/dividedVoidFraction.C b/src/lagrangian/cfdemParticle/subModels/voidFractionModel/dividedVoidFraction/dividedVoidFraction.C index d892d465..3066be22 100644 --- a/src/lagrangian/cfdemParticle/subModels/voidFractionModel/dividedVoidFraction/dividedVoidFraction.C +++ b/src/lagrangian/cfdemParticle/subModels/voidFractionModel/dividedVoidFraction/dividedVoidFraction.C @@ -35,9 +35,7 @@ Description #include "addToRunTimeSelectionTable.H" #include "locateModel.H" #include "dataExchangeModel.H" -#include -//#include // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // namespace Foam @@ -67,29 +65,95 @@ dividedVoidFraction::dividedVoidFraction voidFractionModel(dict,sm), propsDict_(dict.subDict(typeName + "Props")), verbose_(false), + procBoundaryCorrection_(propsDict_.lookupOrDefault("procBoundaryCorrection", false)), alphaMin_(readScalar(propsDict_.lookup("alphaMin"))), alphaLimited_(0), tooMuch_(0.0), interpolation_(false), cfdemUseOnly_(false) { - maxCellsPerParticle_ = 29; + maxCellsPerParticle_ = numberOfMarkerPoints; - if(alphaMin_ > 1 || alphaMin_ < 0.01){ Warning << "alphaMin should be < 1 and > 0.01 !!!" << endl; } - if (propsDict_.found("interpolation")){ - interpolation_=true; + if (alphaMin_ > 1.0 || alphaMin_ < 0.01) + Warning << "alphaMin should be < 1 and > 0.01 !!!" << endl; + + if (propsDict_.found("interpolation")) + { + interpolation_ = true; Warning << "interpolation for dividedVoidFraction does not yet work correctly!" << endl; - Info << "Using interpolated voidfraction field - do not use this in combination with interpolation in drag model!"<< endl; + Info << "Using interpolated voidfraction field - do not use this in combination with interpolation in drag model!" << endl; } checkWeightNporosity(propsDict_); - if (propsDict_.found("verbose")) verbose_=true; + if (propsDict_.found("verbose")) + verbose_ = true; if (propsDict_.found("cfdemUseOnly")) { cfdemUseOnly_ = readBool(propsDict_.lookup("cfdemUseOnly")); } + + if (procBoundaryCorrection_) + { + if (!(particleCloud_.locateM().type() == "engineIB")) + { + FatalError << typeName << ": You are requesting procBoundaryCorrection, this requires the use of engineIB!\n" + << abort(FatalError); + } + } + else + { + if (particleCloud_.locateM().type() == "engineIB") + { + FatalError << typeName << ": You are using engineIB, this requires using procBoundaryCorrection=true!\n" + << abort(FatalError); + } + } + + // generate marker points on unit sphere + label m = 0; + offsets[m][0] = offsets[m][1] = offsets[m][2] = 0.0; + ++m; + + // for 2 different radii + scalar r1 = cbrt(1.0/numberOfMarkerPoints); + scalar r2 = cbrt(15.0/numberOfMarkerPoints); + scalar r[] = { 0.75 * (r2*r2*r2*r2 - r1*r1*r1*r1)/(r2*r2*r2 - r1*r1*r1), + 0.75 * (1.0 - r2*r2*r2*r2)/(1.0 - r2*r2*r2) }; + + for (label ir = 0; ir < 2; ++ir) + { + // try 8 subpoints derived from spherical coordinates + for (scalar zeta = M_PI_4; zeta < constant::mathematical::twoPi; zeta += constant::mathematical::piByTwo) + { + for (scalar theta = M_PI_4; theta < constant::mathematical::pi; theta += constant::mathematical::piByTwo) + { + offsets[m][0] = r[ir] * Foam::sin(theta) * Foam::cos(zeta); + offsets[m][1] = r[ir] * Foam::sin(theta) * Foam::sin(zeta); + offsets[m][2] = r[ir] * Foam::cos(theta); + ++m; + } + } + // try 2 more subpoints for each coordinate direction (6 total) + for (label j = -1; j <= 1; j += 2) + { + offsets[m][0] = r[ir] * j; + offsets[m][1] = 0.; + offsets[m][2] = 0.; + ++m; + + offsets[m][0] = 0.; + offsets[m][1] = r[ir] * j; + offsets[m][2] = 0.; + ++m; + + offsets[m][0] = 0.; + offsets[m][1] = 0.; + offsets[m][2] = r[ir] * j; + ++m; + } + } } @@ -103,7 +167,7 @@ dividedVoidFraction::~dividedVoidFraction() void dividedVoidFraction::setvoidFraction(double** const& mask,double**& voidfractions,double**& particleWeights,double**& particleVolumes, double**& particleV) const { - if(cfdemUseOnly_) + if (cfdemUseOnly_) reAllocArrays(particleCloud_.numberOfParticles()); else reAllocArrays(); @@ -115,19 +179,20 @@ void dividedVoidFraction::setvoidFraction(double** const& mask,double**& voidfra scalar cellVol(0.); scalar scaleVol = weight(); scalar scaleRadius = pow(porosity(),1./3.); + const boundBox& globalBb = particleCloud_.mesh().bounds(); - for(int index=0; index< particleCloud_.numberOfParticles(); index++) + for (int index=0; index < particleCloud_.numberOfParticles(); index++) { - if(!checkParticleType(index)) continue; //skip this particle if not correct type + if (!checkParticleType(index)) continue; //skip this particle if not correct type //if(mask[index][0]) //{ // reset - for(int subcell=0;subcell= 0) // particel centre is in domain { cellVol = particleCloud_.mesh().V()[cellID]; - // for 2 different radii - const scalar delta_r = 0.293976*radius; - for (scalar r = 0.623926*radius; r < radius; r+=delta_r) + if (procBoundaryCorrection_) { - // try 8 subpoint derived from spherical coordinates - for (scalar zeta=M_PI_4; zeta<2.*M_PI; zeta+=M_PI_2) - { - for (scalar theta=M_PI_4; theta(j); - offset[1] = 0.; - offset[2] = 0.; - #include "setWeightedSource.H" //NP set source terms at position+offset + offset = radius * offsets[0]; + #include "setWeightedSource.H" // set source terms at position+offset + } - offset[0] = 0.; - offset[1] = r * static_cast(j); - offset[2] = 0.; - #include "setWeightedSource.H" //NP set source terms at position+offset + for (label i = 1; i < numberOfMarkerPoints; ++i) + { + offset = radius * offsets[i]; + #include "setWeightedSource.H" // set source terms at position+offset + } - offset[0] = 0.; - offset[1] = 0.; - offset[2] = r * static_cast(j); - #include "setWeightedSource.H" //NP set source terms at position+offset - } - }// end loop radiivoidfractions - - if(cellsSet>29 || cellsSet<0) + if (cellsSet > maxCellsPerParticle_ || cellsSet < 0) { Info << "ERROR cellsSet =" << cellsSet << endl; } - // set source for particle center; source 1/nPts+weight of all subpoints that have not been found - scalar centreWeight = 1./nPoints*(nPoints-cellsSet); - - // update voidfraction for each particle read - scalar newAlpha = voidfractionNext_[cellID]- volume*centreWeight/cellVol; - if(newAlpha > alphaMin_) voidfractionNext_[cellID] = newAlpha; - else + if (!procBoundaryCorrection_) { - voidfractionNext_[cellID] = alphaMin_; - tooMuch_ += (alphaMin_-newAlpha) * cellVol; - } + // set source for particle center; source 1/nPts+weight of all subpoints that have not been found + scalar centreWeight = 1./nPoints*(nPoints-cellsSet); - // store cellweight for each particle --- this should be done for subpoints as well!! - particleWeights[index][0] += centreWeight; - - // store particleVolume for each particle - particleVolumes[index][0] += volume*centreWeight; - particleV[index][0] += volume*centreWeight; - - /*//OUTPUT - if (index==0 && verbose_) - { - Info << "centre cellID = " << cellID << endl; - Info << "cellsPerParticle_=" << cellsPerParticle_[index][0] << endl; - - for(int i=0;i alphaMin_) { - if(i==0)Info << "cellids, voidfractions, particleWeights, : \n"; - Info << particleCloud_.cellIDs()[index][i] << " ," << endl; - Info << voidfractions[index][i] << " ," << endl; - Info << particleWeights[index][i] << " ," << endl; - } - }*/ + voidfractionNext_[cellID] = newAlpha; + } + else + { + voidfractionNext_[cellID] = alphaMin_; + tooMuch_ += (alphaMin_-newAlpha) * cellVol; + } + // store cellweight for each particle --- this should be done for subpoints as well!! + particleWeights[index][0] += centreWeight; + + // store particleVolume for each particle + particleVolumes[index][0] += volume*centreWeight; + particleV[index][0] += volume*centreWeight; + } }// end if in cell //}// end if in mask }// end loop all particles + voidfractionNext_.correctBoundaryConditions(); // reset counter of lost volume if (verbose_) Pout << "Total particle volume neglected: " << tooMuch_ << endl; @@ -234,46 +275,17 @@ void dividedVoidFraction::setvoidFraction(double** const& mask,double**& voidfra // bring voidfraction from Eulerian Field to particle array //interpolationCellPoint voidfractionInterpolator_(voidfractionNext_); //scalar voidfractionAtPos(0); - for(int index=0; index< particleCloud_.numberOfParticles(); index++) + for(int index=0; index < particleCloud_.numberOfParticles(); index++) { - /*if(interpolation_) { - label cellI = particleCloud_.cellIDs()[index][0]; - if(cellI >= 0) - { - position = particleCloud_.position(index); - voidfractionAtPos=voidfractionInterpolator_.interpolate(position,cellI); - }else{ - voidfractionAtPos=-1; - } - - for(int subcell=0;subcell= 0) - { - if(voidfractionAtPos > 0) - voidfractions[index][subcell] = voidfractionAtPos; - else - voidfractions[index][subcell] = voidfractionNext_[cellID]; - } - else - { - voidfractions[index][subcell] = -1.; - } - } - } - else*/ - { - for(int subcell=0;subcell= 0) + if (cellID >= 0) { voidfractions[index][subcell] = voidfractionNext_[cellID]; - } + } else { voidfractions[index][subcell] = -1.; diff --git a/src/lagrangian/cfdemParticle/subModels/voidFractionModel/dividedVoidFraction/dividedVoidFraction.H b/src/lagrangian/cfdemParticle/subModels/voidFractionModel/dividedVoidFraction/dividedVoidFraction.H index 1511ac8e..9f4d91a4 100644 --- a/src/lagrangian/cfdemParticle/subModels/voidFractionModel/dividedVoidFraction/dividedVoidFraction.H +++ b/src/lagrangian/cfdemParticle/subModels/voidFractionModel/dividedVoidFraction/dividedVoidFraction.H @@ -57,13 +57,17 @@ class dividedVoidFraction { private: + static const int numberOfMarkerPoints = 29; + dictionary propsDict_; bool verbose_; + Switch procBoundaryCorrection_; + const scalar alphaMin_; // min value of voidFraction - mutable bool alphaLimited_; + mutable bool alphaLimited_; mutable scalar tooMuch_; // particle volume which is lost due to voidFraction limitation @@ -71,6 +75,8 @@ private: bool cfdemUseOnly_; + vector offsets[numberOfMarkerPoints]; + virtual inline scalar Vp(int index, scalar radius, scalar scaleVol) const { return 4.188790205*radius*radius*radius*scaleVol; //4/3*pi=4.188790205 diff --git a/src/lagrangian/cfdemParticle/subModels/voidFractionModel/dividedVoidFraction/setWeightedSource.H b/src/lagrangian/cfdemParticle/subModels/voidFractionModel/dividedVoidFraction/setWeightedSource.H index 36d472c9..8a2ebaa5 100755 --- a/src/lagrangian/cfdemParticle/subModels/voidFractionModel/dividedVoidFraction/setWeightedSource.H +++ b/src/lagrangian/cfdemParticle/subModels/voidFractionModel/dividedVoidFraction/setWeightedSource.H @@ -3,6 +3,18 @@ { // locate subPoint vector subPosition = position + offset; + + if (particleCloud_.checkPeriodicCells()) + { + for (int iDir=0; iDir<3; iDir++) + { + if (subPosition[iDir] > globalBb.max()[iDir]) + subPosition[iDir] -= globalBb.max()[iDir]-globalBb.min()[iDir]; + else if (subPosition[iDir] < globalBb.min()[iDir]) + subPosition[iDir] += globalBb.max()[iDir]-globalBb.min()[iDir]; + } + } + label partCellId = particleCloud_.locateM().findSingleCell(subPosition,cellID); //NP fprintf(lmp->screen,"cellID=%d, partCellId=%d\n",static_cast(cellID),static_cast(partCellId)); diff --git a/src/lagrangian/cfdemParticle/subModels/voidFractionModel/trilinearVoidFraction/trilinearVoidFraction.C b/src/lagrangian/cfdemParticle/subModels/voidFractionModel/trilinearVoidFraction/trilinearVoidFraction.C new file mode 100644 index 00000000..b58e6ad7 --- /dev/null +++ b/src/lagrangian/cfdemParticle/subModels/voidFractionModel/trilinearVoidFraction/trilinearVoidFraction.C @@ -0,0 +1,392 @@ +/*---------------------------------------------------------------------------*\ + CFDEMcoupling - Open Source CFD-DEM coupling + + CFDEMcoupling is part of the CFDEMproject + www.cfdem.com + Christoph Goniva, christoph.goniva@cfdem.com + Copyright 2012- DCS Computing GmbH, Linz +------------------------------------------------------------------------------- +License + This file is part of CFDEMcoupling. + + CFDEMcoupling 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. + + CFDEMcoupling 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 CFDEMcoupling; if not, write to the Free Software Foundation, + Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA + +Description + This code is designed to realize coupled CFD-DEM simulations using LIGGGHTS + and OpenFOAM(R). Note: this code is not part of OpenFOAM(R) (see DISCLAIMER). +\*---------------------------------------------------------------------------*/ + +#include "error.H" + +#include "trilinearVoidFraction.H" +#include "addToRunTimeSelectionTable.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ + +// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * // + +defineTypeNameAndDebug(trilinearVoidFraction, 0); + +addToRunTimeSelectionTable +( + voidFractionModel, + trilinearVoidFraction, + dictionary +); + + +// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // + +// Construct from components +trilinearVoidFraction::trilinearVoidFraction +( + const dictionary& dict, + cfdemCloud& sm +) +: + voidFractionModel(dict,sm), + propsDict_(dict.subDict(typeName + "Props")), + alphaMin_(readScalar(propsDict_.lookup("alphaMin"))), + bb_(particleCloud_.mesh().points(),false), + cellVol_(particleCloud_.mesh().V()[0]), + cellLength_(pow(cellVol_,1./3.)), + nCellXYZ_ + ( + round((bb_.max()[0] - bb_.min()[0]) / cellLength_), + round((bb_.max()[1] - bb_.min()[1]) / cellLength_), + round((bb_.max()[2] - bb_.min()[2]) / cellLength_) + ) +{ + maxCellsPerParticle_ = 8; + checkWeightNporosity(propsDict_); + if (porosity() != 1.) FatalError << "porosity not used in trilinearVoidFraction" << abort(FatalError); + + Warning << "trilinearVoidFraction model is not yet complete and does not work near boundaries" << endl; +} + + +// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * // + +trilinearVoidFraction::~trilinearVoidFraction() +{} + + +// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // + +void trilinearVoidFraction::setvoidFraction(double** const& mask,double**& voidfractions,double**& particleWeights,double**& particleVolumes,double**& particleV) const +{ + reAllocArrays(); + + scalar radius(-1.); + scalar volume(0.); + const scalar scaleVol = weight(); + const scalar fourPiByThree = 4.0*constant::mathematical::pi/3.0; + + vector partPos(0.,0.,0.); + vector pt(0.,0.,0.); + vector posShift(0.,0.,0.); + vector offsetCell(0.,0.,0.); + vector offsetOrigin(0.,0.,0.); + + label i000(0); + label i100(0); + label i110(0); + label i101(0); + label i111(0); + label i010(0); + label i011(0); + label i001(0); + + scalar C000(0.); + scalar C100(0.); + scalar C110(0.); + scalar C101(0.); + scalar C111(0.); + scalar C010(0.); + scalar C011(0.); + scalar C001(0.); + + scalar x(0.); + scalar y(0.); + scalar z(0.); + + scalar a(0.); + scalar b(0.); + scalar c(0.); + + for(int index = 0; index < particleCloud_.numberOfParticles(); ++index) + { + // reset + cellsPerParticle_[index][0] = 8; + //TODO do we need to set particleVolumes, particleV? + // === + + label cellI = particleCloud_.cellIDs()[index][0]; + + if (cellI >= 0) // particel centre is in domain + { + radius = particleCloud_.radius(index); + volume = fourPiByThree * radius * radius * radius * scaleVol; + + // store volume for each particle + particleVolumes[index][0] = volume; + particleV[index][0] = volume; + + // find a,b,c + partPos = particleCloud_.position(index); + offsetCell = partPos - particleCloud_.mesh().C()[cellI]; + a = offsetCell[0]; + b = offsetCell[1]; + c = offsetCell[2]; + + // find "origin" index for mapping + if ( a > 0.) + { + if (b > 0.) + { + if (c > 0.) //FNE + i000 = cellI; + else //BNE + i000 = cellI - nCellXYZ_[0] * nCellXYZ_[1]; + } + else + { + if (c > 0.) //FSE + i000 = cellI - nCellXYZ_[0]; + else //BSE + i000 = cellI - nCellXYZ_[0] - nCellXYZ_[0] * nCellXYZ_[1]; + } + } + else + { + if (b > 0.) + { + if (c > 0.) //FNW + i000 = cellI - 1; + else //BNW + i000 = cellI - 1 - nCellXYZ_[0] * nCellXYZ_[1]; + } + else + { + if (c > 0.) //FSW + i000 = cellI - 1 - nCellXYZ_[0]; + else //BSW + i000 = cellI - 1 - nCellXYZ_[0] - nCellXYZ_[0] * nCellXYZ_[1]; + } + } + + // check boundaries + // TODO different handling for periodic and processor boundaries + pt = particleCloud_.mesh().C()[cellI]; + posShift = vector::zero; + if (a > 0.) + { + pt += vector(cellLength_,0.,0.); + if (pt[0] > bb_.max()[0]) + { + --i000; + posShift[0] = -a; + } + } + else + { + pt -= vector(cellLength_,0.,0.); + if (pt[0] < bb_.min()[0]) + { + ++i000; + posShift[0] = -a; + } + } + if (b > 0.) + { + pt += vector(0.,cellLength_,0.); + if (pt[1] > bb_.max()[1]) + { + i000 -= nCellXYZ_[0]; + posShift[1] = -b; + } + } + else + { + pt -= vector(0.,cellLength_,0.); + if (pt[1] < bb_.min()[1]) + { + i000 += nCellXYZ_[0]; + posShift[1] = -b; + } + } + if (c > 0.) + { + pt += vector(0.,0.,cellLength_); + if (pt[2] > bb_.max()[2]) + { + i000 -= nCellXYZ_[0] * nCellXYZ_[1]; + posShift[2] = -c; + } + } + else + { + pt -= vector(0.,0.,cellLength_); + if (pt[2] < bb_.min()[2]) + { + i000 += nCellXYZ_[0] * nCellXYZ_[1]; + posShift[2] = -c; + } + } + + // define other 7 indices + i100 = i000 + 1; + i110 = i100 + nCellXYZ_[0]; + i010 = i000 + nCellXYZ_[0]; + i001 = i000 + nCellXYZ_[0] * nCellXYZ_[1]; + i101 = i100 + nCellXYZ_[0] * nCellXYZ_[1]; + i111 = i110 + nCellXYZ_[0] * nCellXYZ_[1]; + i011 = i010 + nCellXYZ_[0] * nCellXYZ_[1]; + + // find x,y,z + // TODO changes needed here when generalized for quader cells + offsetOrigin = particleCloud_.mesh().C()[i000] - (partPos + posShift); + x = mag(offsetOrigin[0]) / cellLength_; + y = mag(offsetOrigin[1]) / cellLength_; + z = mag(offsetOrigin[2]) / cellLength_; + + // calculate the mapping coeffs + C000 = (1 - x) * (1 - y) * (1 - z); + C100 = x * (1 - y) * (1 - z); + C110 = x * y * (1 - z); + C010 = (1 - x) * y * (1 - z); + C001 = (1 - x) * (1 - y) * z; + C101 = x * (1 - y) * z; + C111 = x * y * z; + C011 = (1 - x) * y * z; + + // set weights + particleWeights[index][0] = C000; + particleWeights[index][1] = C100; + particleWeights[index][2] = C110; + particleWeights[index][3] = C010; + particleWeights[index][4] = C001; + particleWeights[index][5] = C101; + particleWeights[index][6] = C111; + particleWeights[index][7] = C011; + + // set cellIDs + particleCloud_.cellIDs()[index][0] = i000; + particleCloud_.cellIDs()[index][1] = i100; + particleCloud_.cellIDs()[index][2] = i110; + particleCloud_.cellIDs()[index][3] = i010; + particleCloud_.cellIDs()[index][4] = i001; + particleCloud_.cellIDs()[index][5] = i101; + particleCloud_.cellIDs()[index][6] = i111; + particleCloud_.cellIDs()[index][7] = i011; + + //distribute volume + // TODO use different cell volume when generalized for quader cells + voidfractionNext_[i000] -= volume*C000 / cellVol_; + voidfractionNext_[i100] -= volume*C100 / cellVol_; + voidfractionNext_[i010] -= volume*C010 / cellVol_; + voidfractionNext_[i001] -= volume*C001 / cellVol_; + voidfractionNext_[i101] -= volume*C101 / cellVol_; + voidfractionNext_[i011] -= volume*C011 / cellVol_; + voidfractionNext_[i110] -= volume*C110 / cellVol_; + voidfractionNext_[i111] -= volume*C111 / cellVol_; + + // debugging + /*Pout << "cellI=" << cellI << endl; + Pout << "a=" << a << endl; + Pout << "b=" << b << endl; + Pout << "c=" << c << endl; + Pout << "x=" << x << endl; + Pout << "y=" << y << endl; + Pout << "z=" << z << endl; + Pout << "i000=" << i000 << endl; + Pout << "i100=" << i100 << endl; + Pout << "i010=" << i010 << endl; + Pout << "i001=" << i001 << endl; + Pout << "i101=" << i101 << endl; + Pout << "i011=" << i011 << endl; + Pout << "i110=" << i110 << endl; + Pout << "i111=" << i111 << endl; + + Pout << "C000=" << C000 << endl; + Pout << "C100=" << C100 << endl; + Pout << "C010=" << C010 << endl; + Pout << "C001=" << C001 << endl; + Pout << "C101=" << C101 << endl; + Pout << "C011=" << C011 << endl; + Pout << "C110=" << C110 << endl; + Pout << "C111=" << C111 << endl; + Pout << "sum(Cijk)=" << C000+C100+C010+C001+C101+C011+C110+C111 << endl;*/ + + /*voidfractionNext_[i000]=0.999; + voidfractionNext_[i100]=0.999; + voidfractionNext_[i010]=0.999; + voidfractionNext_[i001]=0.999; + voidfractionNext_[i101]=0.999; + voidfractionNext_[i011]=0.999; + voidfractionNext_[i110]=0.999; + voidfractionNext_[i111]=0.999;*/ + + // limit volumefraction + // TODO implement limiter for all 8 indices + /*if(voidfractionNext_[cellI] < alphaMin_ ) + { + voidfractionNext_[cellI] = alphaMin_; + alphaLimited = true; + } + if(index==0 && alphaLimited) Info<<"alpha limited to" <= 0) + { + voidfractions[index][0] = voidfractionNext_[i000]; + voidfractions[index][1] = voidfractionNext_[i100]; + voidfractions[index][2] = voidfractionNext_[i110]; + voidfractions[index][3] = voidfractionNext_[i010]; + voidfractions[index][4] = voidfractionNext_[i001]; + voidfractions[index][5] = voidfractionNext_[i101]; + voidfractions[index][6] = voidfractionNext_[i111]; + voidfractions[index][7] = voidfractionNext_[i011]; + } + else + { + for (int i = 0; i < 8; ++i) + voidfractions[index][i] = -1.; + } + } +} + + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +} // End namespace Foam + +// ************************************************************************* // diff --git a/src/lagrangian/cfdemParticle/subModels/voidFractionModel/trilinearVoidFraction/trilinearVoidFraction.H b/src/lagrangian/cfdemParticle/subModels/voidFractionModel/trilinearVoidFraction/trilinearVoidFraction.H new file mode 100644 index 00000000..08a3bdd7 --- /dev/null +++ b/src/lagrangian/cfdemParticle/subModels/voidFractionModel/trilinearVoidFraction/trilinearVoidFraction.H @@ -0,0 +1,103 @@ +/*---------------------------------------------------------------------------*\ + CFDEMcoupling - Open Source CFD-DEM coupling + + CFDEMcoupling is part of the CFDEMproject + www.cfdem.com + Christoph Goniva, christoph.goniva@cfdem.com + Copyright 2012- DCS Computing GmbH, Linz +------------------------------------------------------------------------------- +License + This file is part of CFDEMcoupling. + + CFDEMcoupling 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. + + CFDEMcoupling 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 CFDEMcoupling; if not, write to the Free Software Foundation, + Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA + +Description + This code is designed to realize coupled CFD-DEM simulations using LIGGGHTS + and OpenFOAM(R). Note: this code is not part of OpenFOAM(R) (see DISCLAIMER). + +Class + trilinearVoidFraction + +SourceFiles + trilinearVoidFraction.C + +\*---------------------------------------------------------------------------*/ + +#ifndef trilinearVoidFraction_H +#define trilinearVoidFraction_H + +#include "voidFractionModel.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ + +/*---------------------------------------------------------------------------*\ + Class trilinearVoidFraction Declaration +\*---------------------------------------------------------------------------*/ + +class trilinearVoidFraction +: + public voidFractionModel +{ + +private: + dictionary propsDict_; + + const scalar alphaMin_; + + boundBox bb_; + + scalar cellVol_; + + scalar cellLength_; + + vector nCellXYZ_; + +public: + + //- Runtime type information + TypeName("trilinear"); + + + // Constructors + + //- Construct from components + trilinearVoidFraction + ( + const dictionary& dict, + cfdemCloud& sm + ); + + // Destructor + + ~trilinearVoidFraction(); + + + // Member Functions + void setvoidFraction(double** const& ,double**&, double**&, double**&, double**&) const; +}; + + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +} // End namespace Foam + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#endif + +// ************************************************************************* // diff --git a/src/lagrangian/cfdemParticle/subModels/voidFractionModel/voidFractionModel/newVoidFractionModel.C b/src/lagrangian/cfdemParticle/subModels/voidFractionModel/voidFractionModel/newVoidFractionModel.C index 1e4125a2..cf97f414 100644 --- a/src/lagrangian/cfdemParticle/subModels/voidFractionModel/voidFractionModel/newVoidFractionModel.C +++ b/src/lagrangian/cfdemParticle/subModels/voidFractionModel/voidFractionModel/newVoidFractionModel.C @@ -32,7 +32,6 @@ Description #include "error.H" #include "voidFractionModel.H" -#include "centreVoidFraction.H" // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // diff --git a/src/lagrangian/cfdemParticle/subModels/voidFractionModel/voidFractionModel/voidFractionModel.C b/src/lagrangian/cfdemParticle/subModels/voidFractionModel/voidFractionModel/voidFractionModel.C index b5ddbddf..38cb8cdd 100644 --- a/src/lagrangian/cfdemParticle/subModels/voidFractionModel/voidFractionModel/voidFractionModel.C +++ b/src/lagrangian/cfdemParticle/subModels/voidFractionModel/voidFractionModel/voidFractionModel.C @@ -100,70 +100,38 @@ voidFractionModel::~voidFractionModel() } // * * * * * * * * * * * * * * public Member Functions * * * * * * * * * * * * * // -tmp Foam::voidFractionModel::voidFractionInterp() const +tmp voidFractionModel::voidFractionInterp() const { - tmp tsource - ( - new volScalarField - ( - IOobject - ( - "alpha_voidFractionModel", - particleCloud_.mesh().time().timeName(), - particleCloud_.mesh(), - IOobject::NO_READ, - IOobject::NO_WRITE - ), - particleCloud_.mesh(), - dimensionedScalar - ( - "zero", - dimensionSet(0, 0, 0, 0, 0), - 0 - ) - ) - ); - scalar tsf = particleCloud_.dataExchangeM().timeStepFraction(); - if(1-tsf < 1e-4 && particleCloud_.dataExchangeM().couplingStep() > 1) //tsf==1 + + if (1. - tsf < 1e-4 && particleCloud_.dataExchangeM().couplingStep() > 1) //tsf==1 { - tsource.ref() = voidfractionPrev_; + return tmp + ( + new volScalarField("alpha_voidFractionModel", voidfractionPrev_) + ); } else { - tsource.ref() = (1 - tsf) * voidfractionPrev_ + tsf * voidfractionNext_; + return tmp + ( + new volScalarField("alpha_voidFractionModel", (1. - tsf) * voidfractionPrev_ + tsf * voidfractionNext_) + ); } - return tsource; } -void Foam::voidFractionModel::resetVoidFractions() const +void voidFractionModel::resetVoidFractions() const { voidfractionPrev_.ref() = voidfractionNext_.ref(); - voidfractionNext_.ref() = 1; + voidfractionNext_.ref() = 1.; } -/*void Foam::voidFractionModel::undoVoidFractions(double**const& mask) const -{ - voidfractionPrev_.internalField() = voidfractionNext_.internalField(); - - for(int index=0; index< particleCloud_.numberOfParticles(); index++) - { - if(mask[index][0]) - { - // undo voidfraction cause by particle - label cellI = particleCloud_.cellIDs()[index][0]; - scalar cellVolume=voidfractionNext_.mesh().V()[cellI]; - voidfractionNext_[cellI] += particleCloud_.particleVolumes()[index][0]/cellVolume; - } - } -}*/ - int** const& voidFractionModel::cellsPerParticle() const { return cellsPerParticle_; } -int Foam::voidFractionModel::maxCellsPerParticle() const +int voidFractionModel::maxCellsPerParticle() const { return maxCellsPerParticle_; } @@ -186,6 +154,60 @@ void voidFractionModel::reAllocArrays(int nP) const } } +scalar voidFractionModel::pointInParticle(int index, const vector& positionCenter, const vector& point, double scale) const +{ + const scalar radius = particleCloud_.radius(index); + + if(radius > SMALL) + { + scalar pointDistSq = magSqr(point - positionCenter); + return pointDistSq / (scale*scale * radius*radius) - 1.0; + } + else + { + return 0.; + } +} + +//Function to determine minimal distance of point +//to one of the periodic images of a particle +scalar voidFractionModel::minPeriodicDistance(int index, + const vector& cellCentrePosition, + const vector& positionCenter, + const boundBox& globalBb, + vector& minPeriodicPos) const +{ + scalar f = VGREAT; + vector positionCenterPeriodic; + + for(label xDir=-1; xDir<=1; ++xDir) + { + positionCenterPeriodic[0] = positionCenter[0] + + static_cast(xDir) + * (globalBb.max()[0]-globalBb.min()[0]); + for(label yDir=-1; yDir<=1; ++yDir) + { + positionCenterPeriodic[1] = positionCenter[1] + + static_cast(yDir) + * (globalBb.max()[1]-globalBb.min()[1]); + for(label zDir=-1; zDir<=1; ++zDir) + { + positionCenterPeriodic[2] = positionCenter[2] + + static_cast(zDir) + * (globalBb.max()[2]-globalBb.min()[2]); + + if(pointInParticle(index, positionCenterPeriodic, cellCentrePosition) < f) + { + f = pointInParticle(index, positionCenterPeriodic, cellCentrePosition); + minPeriodicPos = positionCenterPeriodic; + } + } + } + } + + return f; +} + // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // } // End namespace Foam diff --git a/src/lagrangian/cfdemParticle/subModels/voidFractionModel/voidFractionModel/voidFractionModel.H b/src/lagrangian/cfdemParticle/subModels/voidFractionModel/voidFractionModel/voidFractionModel.H index 84d88667..877b49bf 100644 --- a/src/lagrangian/cfdemParticle/subModels/voidFractionModel/voidFractionModel/voidFractionModel.H +++ b/src/lagrangian/cfdemParticle/subModels/voidFractionModel/voidFractionModel/voidFractionModel.H @@ -155,6 +155,17 @@ public: virtual void setParticleType(label type) const {} virtual bool checkParticleType(label) const { return true; } //consider all particles by default + + virtual scalar minPeriodicDistance + ( + int index, + const vector& cellCentrePosition, + const vector& positionCenter, + const boundBox& globalBb, + vector& minPeriodicPos + ) const; + + virtual scalar pointInParticle(int index, const vector& positionCenter, const vector& point, double scale=1.0) const; };