diff --git a/src/lagrangian/cfdemParticle/cfdTools/debugInfo.H b/src/lagrangian/cfdemParticle/cfdTools/debugInfo.H index e2071fbc..285e291f 100755 --- a/src/lagrangian/cfdemParticle/cfdTools/debugInfo.H +++ b/src/lagrangian/cfdemParticle/cfdTools/debugInfo.H @@ -1,7 +1,7 @@ - +{ //========================================================================// - scalar countCell=0; // number of cells touched by particles - int points=0; // number of particles and sub-points + label countCell=0; // number of cells touched by particles + label points=0; // number of particles and sub-points scalar totalParticleWeights=0; // total weight of all particles and sub-points vector totalForce_array(0,0,0); // total force on particles based on particle array vector totalForce_field(0,0,0); // forceField of forceM(), used to calc Ksl @@ -93,3 +93,4 @@ Info <<"meanR_array = "<< meanR_array << endl; Info <<"=============================================================================" << endl; Info << endl; +} diff --git a/src/lagrangian/cfdemParticle/cfdTools/mathExtra.H b/src/lagrangian/cfdemParticle/cfdTools/mathExtra.H index fb45ecfd..ee946de0 100644 --- a/src/lagrangian/cfdemParticle/cfdTools/mathExtra.H +++ b/src/lagrangian/cfdemParticle/cfdTools/mathExtra.H @@ -27,7 +27,7 @@ License 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). - + Copyright of this contribution: Copyright 2014- TU Graz, IPPT ------------------------------------------------------------------------- */ @@ -36,15 +36,37 @@ Description #define CFDEM_MATH_EXTRA_H #include -//#include "math.h" #include #include #include #include +#include "mathematicalConstants.H" #define TOLERANCE_ORTHO 1e-10 -namespace MathExtra +namespace Foam +{ +namespace constant +{ +namespace mathematical +{ + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + + const scalar piByFour(0.25*pi); + const scalar fourPiByThree(4.0*pi/3.0); + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +} // End namespace mathematical +} // End namespace constant + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +} // End namespace Foam + + +namespace MathExtra { // inline void outerProduct(double *vec1, double *vec2, double **m); @@ -81,13 +103,13 @@ inline bool spheroidGeometry(double radius, double aspectRatio, //inputs //INPUT // radius ...volume-equivalent radius of the spheroid // aspectRatio ...major/minor aspect ratio - + //OUTPUT - // ai ... - // bi ... - // ei ... - // Le ... - + // ai ... + // bi ... + // ei ... + // Le ... + if(radius<=0.0) //avoid troubles in case radius is 0 or negative return false; @@ -119,7 +141,7 @@ inline bool spheroidGeometry2(double radius, double aspectRatio, //inputs //INPUT // radius ...volume-equivalent radius of the spheroid // aspectRatio ...major/minor aspect ratio - + //OUTPUT // XAe ...Eccentricity dependet parameter // YAe ...Eccentricity dependet parameter @@ -212,7 +234,7 @@ inline bool permutationDotDyadic( //Generate permutation tensor double permutationT[3][3][3]; permutationTensor(permutationT); - + //Step 1: compute dot prodcut of permutation tensor double permutationDotProd[3][3]; zeroize33(permutationDotProd); @@ -227,7 +249,7 @@ inline bool permutationDotDyadic( for(int iY=0; iY<3; iY++) for(int iZ=0; iZ<3; iZ++) tensor[iX][iY][iZ] = permutationDotProd[iX][iY] - * vector[iZ]; + * vector[iZ]; return true; } @@ -240,7 +262,7 @@ inline bool doubleDotTensor333Tensor33(double tensor333[3][3][3], ) { result[0]=0.0;result[1]=0.0;result[2]=0.0; - + for(int iX=0; iX<3; iX++) for(int iY=0; iY<3; iY++) for(int iZ=0; iZ<3; iZ++) diff --git a/src/lagrangian/cfdemParticle/subModels/forceModel/MeiLift/MeiLift.C b/src/lagrangian/cfdemParticle/subModels/forceModel/MeiLift/MeiLift.C index 0dca7e3b..c753e8c1 100644 --- a/src/lagrangian/cfdemParticle/subModels/forceModel/MeiLift/MeiLift.C +++ b/src/lagrangian/cfdemParticle/subModels/forceModel/MeiLift/MeiLift.C @@ -5,7 +5,8 @@ www.cfdem.com Christoph Goniva, christoph.goniva@cfdem.com Copyright 2009-2012 JKU Linz - Copyright 2012- DCS Computing GmbH, Linz + Copyright 2012-2015 DCS Computing GmbH, Linz + Copyright 2015- JKU Linz ------------------------------------------------------------------------------- License This file is part of CFDEMcoupling. @@ -27,6 +28,20 @@ License 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). + + This function is based on the derivation in R. Mei, + An approximate expression for shear lift force on a spherical particle at a + finite Reynolds number, + Int. J. Multiph. Flow 18 (1992) 145–147 + + The data for this functions is based on J.B. Mclaughlin, + Inertial migration of a small sphere in linear shear flows, + Journal of Fluid Mechanics. 224 (1991) 261-274. + + The second order terms are based on E. Loth and A. J. Dorgan, + An equation of motion for particles of finite Reynolds number and size, + Environ. Fluid Mech. 9 (2009) 187–206 + and can be added to the lift coefficient if desired \*---------------------------------------------------------------------------*/ #include "error.H" @@ -84,11 +99,11 @@ MeiLift::MeiLift //Append the field names to be probed 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 - particleCloud_.probeM().scalarFields_.append("Rep"); //other are debug - particleCloud_.probeM().scalarFields_.append("Rew"); //other are debug - particleCloud_.probeM().scalarFields_.append("J_star"); //other are debug + particleCloud_.probeM().vectorFields_.append("Urel"); //other are debug + particleCloud_.probeM().vectorFields_.append("vorticity"); //other are debug + particleCloud_.probeM().scalarFields_.append("Rep"); //other are debug + particleCloud_.probeM().scalarFields_.append("Rew"); //other are debug + particleCloud_.probeM().scalarFields_.append("J_star"); //other are debug particleCloud_.probeM().writeHeader(); } @@ -124,6 +139,7 @@ void MeiLift::setForce() const scalar J_star(0); scalar Omega_eq(0); scalar alphaStar(0); + scalar epsilonSqr(0.0); scalar epsilon(0); scalar omega_star(0); vector vorticity(0,0,0); @@ -134,7 +150,7 @@ 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]) //{ @@ -145,18 +161,16 @@ void MeiLift::setForce() const { Us = particleCloud_.velocity(index); - if( forceSubM(0).interpolation() ) + if (forceSubM(0).interpolation()) { position = particleCloud_.position(index); - Ur = UInterpolator_.interpolate(position,cellI) - - Us; + Ur = UInterpolator_.interpolate(position,cellI) - Us; vorticity = VorticityInterpolator_.interpolate(position,cellI); } else { - Ur = U_[cellI] - - Us; - vorticity=vorticityField[cellI]; + Ur = U_[cellI] - Us; + vorticity = vorticityField[cellI]; } magUr = mag(Ur); @@ -164,7 +178,7 @@ void MeiLift::setForce() const if (magUr > 0 && magVorticity > 0) { - ds = 2*particleCloud_.radius(index); + ds = 2. * particleCloud_.radius(index); nuf = nufField[cellI]; rho = rhoField[cellI]; @@ -172,46 +186,54 @@ void MeiLift::setForce() const Rep = ds*magUr/nuf; Rew = magVorticity*ds*ds/nuf; - alphaStar = magVorticity*ds/magUr/2.0; - epsilon = sqrt(2.0*alphaStar /Rep ); - omega_star = 2.0*alphaStar; + omega_star = magVorticity * ds / magUr; + alphaStar = 0.5 * omega_star; + epsilonSqr = omega_star / Rep; + epsilon = sqrt(epsilonSqr); //Basic model for the correction to the Saffman lift - //Based on McLaughlin (1991) - if(epsilon < 0.1) + //McLaughlin (1991), Mei (1992), Loth and Dorgan (2009) + //J_star = 0.443 * J + if (epsilon < 0.1) //epsilon << 1 { - J_star = -140 *epsilon*epsilon*epsilon*epsilon*epsilon - *log( 1./(epsilon*epsilon+SMALL) ); + //McLaughlin (1991), Eq (3.27): J = 32 * pi^2 * epsilon^5 * ln(1 / epsilon^2) + J_star = -140.0 * epsilonSqr * epsilonSqr * epsilon * log(1. / (epsilonSqr+SMALL)); } - else if(epsilon > 20) + else if (epsilon > 20.0) //epsilon >> 1 { - J_star = 1.0-0.287/(epsilon*epsilon+SMALL); + //McLaughlin (1991), Eq (3.26): J = 2.255 - 0.6463 / epsilon^2 + J_star = 1.0 - 0.287 / epsilonSqr; } else { + //Mei (1992), Eq (10) + //Loth and Dorgan (2009), Eq (32) J_star = 0.3 - *( 1.0 - +tanh( 2.5 * log10(epsilon+0.191) ) - ) - *( 0.667 - +tanh( 6.0 * (epsilon-0.32) ) - ); + * (1.0 + tanh(2.5 * (log10(epsilon) + 0.191))) + * (0.667 + tanh(6.0 * ( epsilon - 0.32 ))); } - Cl = J_star * 4.11 * epsilon; //multiply McLaughlin's correction to the basic Saffman model + //Loth and Dorgan (2009), Eq (31): Saffman lift-coefficient: ClSaff = 12.92 / pi * epsilon ~ 4.11 * epsilon + //Loth and Dorgan (2009), Eq (32) + Cl = J_star * 4.11 * epsilon; //multiply correction to the basic Saffman model - //Second order terms given by Loth and Dorgan 2009 - if(useSecondOrderTerms_) + //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; + scalar sqrtRep = sqrt(Rep); + //Loth and Dorgan (2009), Eq (34) + Cl_star = 1.0 - (0.675 + 0.15 * (1.0 + tanh(0.28 * (alphaStar - 2.0)))) * tanh(0.18 * sqrtRep); + //Loth and Dorgan (2009), Eq (38) + Omega_eq = alphaStar * (1.0 - 0.0075 * Rew) * (1.0 - 0.062 * sqrtRep - 0.001 * Rep); + //Loth and Dorgan (2009), Eq (39) + Cl += Omega_eq * Cl_star; } - lift = 0.125*M_PI - *rho - *Cl - *magUr*Ur^vorticity/magVorticity - *ds*ds; + //Loth and Dorgan (2009), Eq (27) + lift = 0.125 * constant::mathematical::pi + * rho + * Cl + * magUr * Ur ^ vorticity / magVorticity + * ds * ds; if (modelType_ == "B") { diff --git a/src/lagrangian/cfdemParticle/subModels/forceModel/MeiLift/MeiLift.H b/src/lagrangian/cfdemParticle/subModels/forceModel/MeiLift/MeiLift.H index 6a6016b9..81ecfaee 100644 --- a/src/lagrangian/cfdemParticle/subModels/forceModel/MeiLift/MeiLift.H +++ b/src/lagrangian/cfdemParticle/subModels/forceModel/MeiLift/MeiLift.H @@ -5,7 +5,8 @@ www.cfdem.com Christoph Goniva, christoph.goniva@cfdem.com Copyright 2009-2012 JKU Linz - Copyright 2012- DCS Computing GmbH, Linz + Copyright 2012-2015 DCS Computing GmbH, Linz + Copyright 2015- JKU Linz ------------------------------------------------------------------------------- License This file is part of CFDEMcoupling. @@ -28,18 +29,20 @@ 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). - This function is based on the derivation in R. Kurose, S. Komori, - Drag and lift forces on a rotating sphere in a linear shear flow, - Journal of Fluid Mechanics. 384 (1999) 183-206. + This function is based on the derivation in R. Mei, + An approximate expression for shear lift force on a spherical particle at a + finite Reynolds number, + Int. J. Multiph. Flow 18 (1992) 145–147 - The data for this functions is based on J.B. Mclaughlin, - Inertial migration of a small sphere in linear shear flows, - Journal of Fluid Mechanics. 224 (1991) 261-274. + The data for this functions is based on J.B. Mclaughlin, + Inertial migration of a small sphere in linear shear flows, + Journal of Fluid Mechanics. 224 (1991) 261-274. - The second order terms are based on: - Mei Lift force following Loth and Dorgan 2009, - and can be added to the lift coefficient if desired - (contribution from RQ) + The second order terms are based on E. Loth and A. J. Dorgan, + An equation of motion for particles of finite Reynolds number and size, + Environ. Fluid Mech. 9 (2009) 187–206 + and can be added to the lift coefficient if desired + (contribution from RQ) - including interpolation of the velocity to the particle position (optional) - including output to file for testing/data analysis (optional) diff --git a/src/lagrangian/cfdemParticle/subModels/locateModel/engineSearchIB/engineSearchIB.C b/src/lagrangian/cfdemParticle/subModels/locateModel/engineSearchIB/engineSearchIB.C index 3e79047e..a4b95a3f 100644 --- a/src/lagrangian/cfdemParticle/subModels/locateModel/engineSearchIB/engineSearchIB.C +++ b/src/lagrangian/cfdemParticle/subModels/locateModel/engineSearchIB/engineSearchIB.C @@ -30,14 +30,10 @@ Description \*---------------------------------------------------------------------------*/ #include "error.H" - - #include "engineSearchIB.H" #include "addToRunTimeSelectionTable.H" #include "mathematicalConstants.H" -#include - // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // namespace Foam @@ -64,15 +60,19 @@ engineSearchIB::engineSearchIB cfdemCloud& sm ) : - engineSearch(dict.subDict(typeName + "Props"),sm), + engineSearch(dict.subDict(typeName + "Props"), sm), propsDict_(dict.subDict(typeName + "Props")), zSplit_(readLabel(propsDict_.lookup("zSplit"))), xySplit_(readLabel(propsDict_.lookup("xySplit"))), - checkPeriodicCells_(false) + thetaSize_(180./zSplit_), + phiSize_(360./xySplit_), + deg2rad_(constant::mathematical::pi/180.), + numberOfSatellitePoints_((zSplit_-1)*xySplit_+2) { - - if(propsDict_.found("checkPeriodicCells")) checkPeriodicCells_=true; - + for (int countPoints = 0; countPoints < numberOfSatellitePoints_; ++countPoints) + { + satellitePoints_.push_back(generateSatellitePoint(countPoints)); + } } @@ -93,87 +93,91 @@ label engineSearchIB::findCell int size ) const { + bool checkPeriodicCells(particleCloud_.checkPeriodicCells()); const boundBox& globalBb = particleCloud_.mesh().bounds(); vector position; - for(int index = 0;index < size; ++index) + for (int index = 0; index < size; ++index) { - cellIDs[index][0]=-1; - double radius=particleCloud_.radius(index); - //if(mask[index][0] && radius > SMALL) - if(radius > SMALL) + cellIDs[index][0] = -1; + double radius = particleCloud_.radius(index); + + if (radius > SMALL) { // create pos vector - for(int i=0;i<3;i++) position[i] = positions[index][i]; + for (int i = 0; i < 3; i++) position[i] = positions[index][i]; // find cell label oldID = cellIDs[index][0]; - cellIDs[index][0] = findSingleCell(position,oldID); - //cellIDs[index][0] = particleCloud_.mesh().findCell(position); + cellIDs[index][0] = findSingleCell(position, oldID); - //mod by alice upon from here - if(cellIDs[index][0] < 0) + if (cellIDs[index][0] < 0) { - vector pos = position; label altStartPos = -1; - label numberOfPoints = (zSplit_-1)*xySplit_ + 2; // 1 point at bottom, 1 point at top - label thetaLevel = 0; - scalar theta, phi; - const scalar thetaSize = 180./zSplit_, phiSize = 360./xySplit_; - const scalar deg2rad = M_PI/180.; - for(int countPoints = 0; countPoints < numberOfPoints; ++countPoints) + for (int countPoints = 0; countPoints < numberOfSatellitePoints_; ++countPoints) { - pos = position; - if(countPoints == 0) - { - pos[2] += radius; - } - else if(countPoints == 1) - { - pos[2] -= radius; - } - else - { - thetaLevel = (countPoints - 2) / xySplit_; - theta = deg2rad * thetaSize * (thetaLevel+1); - phi = deg2rad * phiSize * (countPoints - 2 - thetaLevel*xySplit_); - pos[0] += radius * sin(theta) * cos(phi); - pos[1] += radius * sin(theta) * sin(phi); - pos[2] += radius * cos(theta); - } + vector pos = getSatellitePoint(index, countPoints); + + altStartPos = findSingleCell(pos,oldID); - altStartPos=findSingleCell(pos,oldID); //particleCloud_.mesh().findCell(pos);// //check for periodic domains - if(checkPeriodicCells_) + if (checkPeriodicCells) { - for(int iDir=0;iDir<3;iDir++) + for (int iDir = 0; iDir < 3; iDir++) { - if( pos[iDir] > globalBb.max()[iDir] ) + if (pos[iDir] > globalBb.max()[iDir]) { - pos[iDir]-=globalBb.max()[iDir]-globalBb.min()[iDir]; + pos[iDir] -= globalBb.max()[iDir] - globalBb.min()[iDir]; } - else if( pos[iDir] < globalBb.min()[iDir] ) + else if (pos[iDir] < globalBb.min()[iDir]) { - pos[iDir]+=globalBb.max()[iDir]-globalBb.min()[iDir]; + pos[iDir] += globalBb.max()[iDir] - globalBb.min()[iDir]; } } - altStartPos=findSingleCell(pos,oldID); //particleCloud_.mesh().findCell(pos);// + + altStartPos = findSingleCell(pos, oldID); } - if(altStartPos >= 0) // found position, we're done + if (altStartPos >= 0) // found position, we're done { cellIDs[index][0] = altStartPos; break; } - } + } } } } + return 1; } +vector engineSearchIB::generateSatellitePoint(int countPoints) const +{ + // 1 point at bottom, 1 point at top + if (countPoints == 0) + { + return vector(0., 0., 1.); + } + else if (countPoints == 1) + { + return vector(0., 0., -1.); + } + else + { + const scalar thetaLevel = (countPoints - 2) / xySplit_; + const scalar theta = deg2rad_ * thetaSize_ * (thetaLevel + 1); + const scalar phi = deg2rad_ * phiSize_ * (countPoints - 2 - thetaLevel * xySplit_); + return vector(sin(theta) * cos(phi), sin(theta) * sin(phi), cos(theta)); + } +} + +vector engineSearchIB::getSatellitePoint(int index, int countPoints) const +{ + return particleCloud_.position(index) + + particleCloud_.radius(index) * satellitePoints_[countPoints]; +} // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // diff --git a/src/lagrangian/cfdemParticle/subModels/locateModel/engineSearchIB/engineSearchIB.H b/src/lagrangian/cfdemParticle/subModels/locateModel/engineSearchIB/engineSearchIB.H index 5a7bda5f..8dbdb4e7 100644 --- a/src/lagrangian/cfdemParticle/subModels/locateModel/engineSearchIB/engineSearchIB.H +++ b/src/lagrangian/cfdemParticle/subModels/locateModel/engineSearchIB/engineSearchIB.H @@ -47,7 +47,7 @@ namespace Foam { /*---------------------------------------------------------------------------*\ - Class noDrag Declaration + Class engineSearchIB Declaration \*---------------------------------------------------------------------------*/ class engineSearchIB @@ -59,11 +59,17 @@ private: dictionary propsDict_; - const label zSplit_; + const label zSplit_; - const label xySplit_; - - bool checkPeriodicCells_; + const label xySplit_; + + const scalar thetaSize_; + const scalar phiSize_; + const scalar deg2rad_; + + const label numberOfSatellitePoints_; + + std::vector satellitePoints_; public: @@ -94,6 +100,17 @@ public: int size ) const; +private: + vector generateSatellitePoint + ( + int countPoints + ) const; + + vector getSatellitePoint + ( + int index, + int countPoints + ) const; }; diff --git a/src/lagrangian/cfdemParticle/subModels/voidFractionModel/GaussVoidFraction/GaussVoidFraction.C b/src/lagrangian/cfdemParticle/subModels/voidFractionModel/GaussVoidFraction/GaussVoidFraction.C index 231e7e97..3964b074 100644 --- a/src/lagrangian/cfdemParticle/subModels/voidFractionModel/GaussVoidFraction/GaussVoidFraction.C +++ b/src/lagrangian/cfdemParticle/subModels/voidFractionModel/GaussVoidFraction/GaussVoidFraction.C @@ -30,7 +30,7 @@ Description \*---------------------------------------------------------------------------*/ #include "error.H" - +#include "mathExtra.H" #include "GaussVoidFraction.H" #include "addToRunTimeSelectionTable.H" #include "locateModel.H" @@ -69,7 +69,7 @@ GaussVoidFraction::GaussVoidFraction alphaLimited_(0) { 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; + FatalError << "\n\n This model does not yet work properly! \n\n" << endl; //reading maxCellsPerParticle from dictionary maxCellsPerParticle_=readLabel(propsDict_.lookup("maxCellsPerParticle")); @@ -115,12 +115,12 @@ void GaussVoidFraction::setvoidFraction(double** const& mask,double**& voidfract label particleCenterCellID=particleCloud_.cellIDs()[index][0]; radius = particleCloud_.radius(index); - volume = 4.188790205*radius*radius*radius*scaleVol; + volume = constant::mathematical::fourPiByThree*radius*radius*radius*scaleVol; radius *= scaleRadius; vector positionCenter=particleCloud_.position(index); - scalar core; - scalar dist; + scalar core; + scalar dist; if (particleCenterCellID >= 0) { @@ -151,8 +151,8 @@ void GaussVoidFraction::setvoidFraction(double** const& mask,double**& voidfract //==========================// //setting the voidfractions - dist = mag(particleCloud_.mesh().C()[particleCenterCellID]-particleCloud_.position(index)); - core = pow(2.0/radius/radius/M_PI,1.5)*exp(-dist*dist/2.0/radius/radius)*particleCloud_.mesh().V()[particleCenterCellID]; + dist = mag(particleCloud_.mesh().C()[particleCenterCellID]-particleCloud_.position(index)); + core = pow(2.0/radius/radius/M_PI,1.5)*exp(-dist*dist/2.0/radius/radius)*particleCloud_.mesh().V()[particleCenterCellID]; // volume occupied in every covered cell scalar occupiedVolume = volume*core; @@ -164,8 +164,8 @@ void GaussVoidFraction::setvoidFraction(double** const& mask,double**& voidfract particleVolumes[index][0] += occupiedVolume; particleV[index][0] += occupiedVolume; - //Info << "Centre:set voidfraction in cellI=" << particleCenterCellID - // << ", voidfraction =" << voidfractionNext_[particleCenterCellID] << endl; + //Info << "Centre:set voidfraction in cellI=" << particleCenterCellID + // << ", voidfraction =" << voidfractionNext_[particleCenterCellID] << endl; // correct volumefraction of sub-cells for(label i=0;i 0.0) + coronaPoint += (cellCentrePosition - minPeriodicParticlePos) * (corona / centreDist); if (pointInParticle(index, minPeriodicParticlePos, coronaPoint) < 0.0) { @@ -249,7 +252,7 @@ void IBVoidFraction::setvoidFraction(double** const& mask,double**& voidfraction label partCellId = particleCloud_.mesh().findNearestCell(particlePosList[iPeriodicImage]); particleLabelList.append(partCellId); - buildLabelHashSet(radius, particlePosList[iPeriodicImage], particleLabelList[iPeriodicImage], hashSett, false); + buildLabelHashSet(index, particlePosList[iPeriodicImage], particleLabelList[iPeriodicImage], hashSett, false); } } //end checkPeriodicCells_ @@ -316,8 +319,10 @@ void IBVoidFraction::buildLabelHashSet scalar centreDist = mag(cellCentrePosition-position); 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); + scalar corona = 0.5 * sqrtThree_ * pow(particleCloud_.mesh().V()[neighbor], 1./3.); + vector coronaPoint = cellCentrePosition; + if (centreDist > 0.0) + coronaPoint += (cellCentrePosition - position) * (corona / centreDist); if (!hashSett.found(neighbor) && pointInParticle(index, position, coronaPoint) < 0.0) { diff --git a/src/lagrangian/cfdemParticle/subModels/voidFractionModel/IBVoidFraction/IBVoidFraction.H b/src/lagrangian/cfdemParticle/subModels/voidFractionModel/IBVoidFraction/IBVoidFraction.H index e2072025..0e430dc9 100755 --- a/src/lagrangian/cfdemParticle/subModels/voidFractionModel/IBVoidFraction/IBVoidFraction.H +++ b/src/lagrangian/cfdemParticle/subModels/voidFractionModel/IBVoidFraction/IBVoidFraction.H @@ -68,6 +68,8 @@ private: const scalar scaleUpVol_; //NP scaling radius, keeping volume of particle + const scalar sqrtThree_; + public: //- Runtime type information diff --git a/src/lagrangian/cfdemParticle/subModels/voidFractionModel/bigParticleVoidFraction/bigParticleVoidFraction.C b/src/lagrangian/cfdemParticle/subModels/voidFractionModel/bigParticleVoidFraction/bigParticleVoidFraction.C index 63405dca..d176668f 100644 --- a/src/lagrangian/cfdemParticle/subModels/voidFractionModel/bigParticleVoidFraction/bigParticleVoidFraction.C +++ b/src/lagrangian/cfdemParticle/subModels/voidFractionModel/bigParticleVoidFraction/bigParticleVoidFraction.C @@ -30,7 +30,7 @@ Description \*---------------------------------------------------------------------------*/ #include "error.H" - +#include "mathExtra.H" #include "bigParticleVoidFraction.H" #include "addToRunTimeSelectionTable.H" #include "locateModel.H" @@ -113,7 +113,7 @@ void bigParticleVoidFraction::setvoidFraction(double** const& mask,double**& voi //collecting data label particleCenterCellID=particleCloud_.cellIDs()[index][0]; radius = particleCloud_.radius(index); - volume = 4.188790205*radius*radius*radius*scaleVol; + volume = constant::mathematical::fourPiByThree*radius*radius*radius*scaleVol; radius *= scaleRadius; vector positionCenter=particleCloud_.position(index); diff --git a/src/lagrangian/cfdemParticle/subModels/voidFractionModel/centreVoidFraction/centreVoidFraction.C b/src/lagrangian/cfdemParticle/subModels/voidFractionModel/centreVoidFraction/centreVoidFraction.C index a9237811..0f1fb2bf 100644 --- a/src/lagrangian/cfdemParticle/subModels/voidFractionModel/centreVoidFraction/centreVoidFraction.C +++ b/src/lagrangian/cfdemParticle/subModels/voidFractionModel/centreVoidFraction/centreVoidFraction.C @@ -30,8 +30,8 @@ Description \*---------------------------------------------------------------------------*/ #include "error.H" - #include "centreVoidFraction.H" +#include "mathExtra.H" #include "addToRunTimeSelectionTable.H" // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // @@ -101,7 +101,7 @@ void centreVoidFraction::setvoidFraction(double** const& mask,double**& voidfrac { cellVol = voidfractionNext_.mesh().V()[cellI]; radius = particleCloud_.radius(index); - volume = 4.188790205*radius*radius*radius*scaleVol; + volume = constant::mathematical::fourPiByThree*radius*radius*radius*scaleVol; // store volume for each particle particleVolumes[index][0] = volume; diff --git a/src/lagrangian/cfdemParticle/subModels/voidFractionModel/dividedVoidFraction/dividedVoidFraction.H b/src/lagrangian/cfdemParticle/subModels/voidFractionModel/dividedVoidFraction/dividedVoidFraction.H index 9f4d91a4..e14062fb 100644 --- a/src/lagrangian/cfdemParticle/subModels/voidFractionModel/dividedVoidFraction/dividedVoidFraction.H +++ b/src/lagrangian/cfdemParticle/subModels/voidFractionModel/dividedVoidFraction/dividedVoidFraction.H @@ -39,6 +39,7 @@ SourceFiles #ifndef dividedVoidFraction_H #define dividedVoidFraction_H +#include "mathExtra.H" #include "voidFractionModel.H" #include "interpolationCellPoint.H" @@ -79,8 +80,8 @@ private: virtual inline scalar Vp(int index, scalar radius, scalar scaleVol) const { - return 4.188790205*radius*radius*radius*scaleVol; //4/3*pi=4.188790205 - }; + return constant::mathematical::fourPiByThree*radius*radius*radius*scaleVol; + } public: diff --git a/src/lagrangian/cfdemParticle/subModels/voidFractionModel/dividedVoidFraction/setWeightedSource.H b/src/lagrangian/cfdemParticle/subModels/voidFractionModel/dividedVoidFraction/setWeightedSource.H index 8a2ebaa5..d8a418db 100755 --- a/src/lagrangian/cfdemParticle/subModels/voidFractionModel/dividedVoidFraction/setWeightedSource.H +++ b/src/lagrangian/cfdemParticle/subModels/voidFractionModel/dividedVoidFraction/setWeightedSource.H @@ -67,7 +67,7 @@ } - particleWeights[index][storeInIndex] += 1/static_cast(nPoints); + particleWeights[index][storeInIndex] += 1.0/nPoints; particleVolumes[index][storeInIndex] += particleVolume; particleV[index][0] += particleVolume; //====================================================// diff --git a/src/lagrangian/cfdemParticle/subModels/voidFractionModel/trilinearVoidFraction/trilinearVoidFraction.C b/src/lagrangian/cfdemParticle/subModels/voidFractionModel/trilinearVoidFraction/trilinearVoidFraction.C index b58e6ad7..5e395e31 100644 --- a/src/lagrangian/cfdemParticle/subModels/voidFractionModel/trilinearVoidFraction/trilinearVoidFraction.C +++ b/src/lagrangian/cfdemParticle/subModels/voidFractionModel/trilinearVoidFraction/trilinearVoidFraction.C @@ -29,7 +29,7 @@ Description \*---------------------------------------------------------------------------*/ #include "error.H" - +#include "mathExtra.H" #include "trilinearVoidFraction.H" #include "addToRunTimeSelectionTable.H" @@ -95,7 +95,6 @@ void trilinearVoidFraction::setvoidFraction(double** const& mask,double**& voidf 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.); @@ -141,7 +140,7 @@ void trilinearVoidFraction::setvoidFraction(double** const& mask,double**& voidf if (cellI >= 0) // particel centre is in domain { radius = particleCloud_.radius(index); - volume = fourPiByThree * radius * radius * radius * scaleVol; + volume = constant::mathematical::fourPiByThree * radius * radius * radius * scaleVol; // store volume for each particle particleVolumes[index][0] = volume;