Merge pull request #28 from ParticulateFlow/feature/cherry_pick_PUBLIC

Feature/cherry pick public

src/lagrangian/cfdemParticle/cfdTools/debugInfo.H

@@ -1,7 +1,7 @@
-
+{
   //========================================================================//
-  scalar countCell=0;              // number of cells touched by particles
+  label countCell=0;               // number of cells touched by particles
-  int points=0;                    // number of particles and sub-points
+  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;
+}

src/lagrangian/cfdemParticle/cfdTools/mathExtra.H

@@ -36,14 +36,36 @@ Description
 #define CFDEM_MATH_EXTRA_H
 
 #include <vector>
-//#include "math.h"
 #include <stdio.h>
 #include <string.h>
 #include <error.h>
 #include <ctype.h>
+#include "mathematicalConstants.H"
 
 #define TOLERANCE_ORTHO 1e-10
 
+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
 {
 

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"
@@ -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);
@@ -148,14 +164,12 @@ void MeiLift::setForce() const
                 if (forceSubM(0).interpolation())
                 {
                     position  = particleCloud_.position(index);
-                    Ur        = UInterpolator_.interpolate(position,cellI)
+                    Ur        = UInterpolator_.interpolate(position,cellI) - Us;
-                                - Us;
                     vorticity = VorticityInterpolator_.interpolate(position,cellI);
                 }
                 else
                 {
-                    Ur =  U_[cellI]
+                    Ur        = U_[cellI] - Us;
-                          - Us;
                     vorticity = vorticityField[cellI];
                 }
 
@@ -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,42 +186,50 @@ void MeiLift::setForce() const
                     Rep = ds*magUr/nuf;
                     Rew = magVorticity*ds*ds/nuf;
 
-                    alphaStar  = magVorticity*ds/magUr/2.0;
+                    omega_star = magVorticity * ds / magUr;
-                    epsilon    = sqrt(2.0*alphaStar /Rep );
+                    alphaStar = 0.5 * omega_star;
-                    omega_star = 2.0*alphaStar;
+                    epsilonSqr = omega_star / Rep;
+                    epsilon = sqrt(epsilonSqr);
 
                     //Basic model for the correction to the Saffman lift
-                    //Based on McLaughlin (1991)
+                    //McLaughlin (1991), Mei (1992), Loth and Dorgan (2009)
-                    if(epsilon < 0.1)
+                    //J_star = 0.443 * J
+                    if (epsilon < 0.1) //epsilon << 1
                     {
-                        J_star = -140 *epsilon*epsilon*epsilon*epsilon*epsilon
+                        //McLaughlin (1991), Eq (3.27): J = 32 * pi^2 * epsilon^5 * ln(1 / epsilon^2)
-                                             *log( 1./(epsilon*epsilon+SMALL) );
+                        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
+                                 * (1.0   + tanh(2.5 * (log10(epsilon) + 0.191)))
-                                      +tanh(  2.5 * log10(epsilon+0.191)  )
+                                 * (0.667 + tanh(6.0 * (      epsilon  - 0.32 )));
-                                 )
-                                *(    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
+                    //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);
+                        scalar sqrtRep = sqrt(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));
+                        //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
+                    //Loth and Dorgan (2009), Eq (27)
+                    lift = 0.125 * constant::mathematical::pi
                            * rho
                            * Cl
                            * magUr * Ur ^ vorticity / magVorticity

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,16 +29,18 @@ 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, 
+    This function is based on the derivation in R. Mei,
-     Drag and lift forces on a rotating sphere in a linear shear flow, 
+    An approximate expression for shear lift force on a spherical  particle at a
-     Journal of Fluid Mechanics. 384 (1999) 183-206.
+    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: 
+    The second order terms are based on E. Loth and A. J. Dorgan,
-     Mei Lift force following Loth and Dorgan 2009,
+    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)
 

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 <mpi.h>
-
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 namespace Foam
@@ -68,11 +64,15 @@ engineSearchIB::engineSearchIB
     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)
 {
-
+    for (int countPoints = 0; countPoints < numberOfSatellitePoints_; ++countPoints)
-    if(propsDict_.found("checkPeriodicCells")) checkPeriodicCells_=true;
+    {
-
+        satellitePoints_.push_back(generateSatellitePoint(countPoints));
+    }
 }
 
 
@@ -93,6 +93,7 @@ label engineSearchIB::findCell
     int size
 ) const
 {
+    bool checkPeriodicCells(particleCloud_.checkPeriodicCells());
     const boundBox& globalBb = particleCloud_.mesh().bounds();
 
     vector position;
@@ -100,7 +101,7 @@ label engineSearchIB::findCell
     {
         cellIDs[index][0] = -1;
         double radius = particleCloud_.radius(index);
-        //if(mask[index][0] && radius > SMALL)
+
         if (radius > SMALL)
         {
             // create pos vector
@@ -109,43 +110,19 @@ label engineSearchIB::findCell
             // find cell
             label oldID = cellIDs[index][0];
             cellIDs[index][0] = findSingleCell(position, oldID);
-            //cellIDs[index][0] = particleCloud_.mesh().findCell(position);
 
-            //mod by alice upon from here
             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;
+                    vector pos = getSatellitePoint(index, countPoints);
-                    if(countPoints == 0)
+
-                    {
+                    altStartPos = findSingleCell(pos,oldID);
-                        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);
-                    }
 
-            		altStartPos=findSingleCell(pos,oldID); //particleCloud_.mesh().findCell(pos);//
                     //check for periodic domains
-                    if(checkPeriodicCells_)
+                    if (checkPeriodicCells)
                     {
                         for (int iDir = 0; iDir < 3; iDir++)
                         {
@@ -158,7 +135,8 @@ label engineSearchIB::findCell
                                 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
@@ -171,9 +149,35 @@ label engineSearchIB::findCell
             }
         }
     }
+
     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];
+}
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 

src/lagrangian/cfdemParticle/subModels/locateModel/engineSearchIB/engineSearchIB.H

@@ -47,7 +47,7 @@ namespace Foam
 {
 
 /*---------------------------------------------------------------------------*\
-                           Class noDrag Declaration
+                           Class engineSearchIB Declaration
 \*---------------------------------------------------------------------------*/
 
 class engineSearchIB
@@ -63,7 +63,13 @@ private:
 
         const label xySplit_;
 
-    	bool checkPeriodicCells_;
+        const scalar thetaSize_;
+        const scalar phiSize_;
+        const scalar deg2rad_;
+
+        const label numberOfSatellitePoints_;
+
+        std::vector<vector> satellitePoints_;
 
 public:
 
@@ -94,6 +100,17 @@ public:
             int size
         ) const;
 
+private:
+        vector generateSatellitePoint
+        (
+            int countPoints
+        ) const;
+
+        vector getSatellitePoint
+        (
+            int index,
+            int countPoints
+        ) const;
 };
 
 

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,7 +115,7 @@ 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);
@@ -183,7 +183,6 @@ void GaussVoidFraction::setvoidFraction(double** const& mask,double**& voidfract
 
                         //Info << "AFTER:set voidfraction in cellI=" << cellI
                         //     << ", voidfraction =" << voidfractionNext_[cellI] << endl;
-
                     }
 
                     // debug
@@ -235,12 +234,13 @@ void GaussVoidFraction::buildLabelHashSet
     forAll(nc,i)
     {
         label neighbor=nc[i];
-        if(!hashSett.found(neighbor) && mag(position-particleCloud_.mesh().C()[neighbor])<radius){
+        if(!hashSett.found(neighbor) && mag(position-particleCloud_.mesh().C()[neighbor])<radius)
+        {
             buildLabelHashSet(radius,position,neighbor,hashSett);
         }
     }
-
 }
+
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 } // End namespace Foam

src/lagrangian/cfdemParticle/subModels/voidFractionModel/IBVoidFraction/IBVoidFraction.C

@@ -66,7 +66,8 @@ IBVoidFraction::IBVoidFraction
     propsDict_(dict.subDict(typeName + "Props")),
     alphaMin_(readScalar(propsDict_.lookup("alphaMin"))),
     alphaLimited_(0),
-    scaleUpVol_(readScalar(propsDict_.lookup("scaleUpVol")))
+    scaleUpVol_(readScalar(propsDict_.lookup("scaleUpVol"))),
+    sqrtThree_(sqrt(3.0))
 {
     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"));
@@ -128,8 +129,10 @@ void IBVoidFraction::setvoidFraction(double** const& mask,double**& voidfraction
                 }
 
                 scalar centreDist = mag(cellCentrePosition - minPeriodicParticlePos);
-                scalar corona = 0.5 * sqrt(3.0) * pow(particleCloud_.mesh().V()[particleCenterCellID], 1./3.);
+                scalar corona = 0.5 * sqrtThree_ * pow(particleCloud_.mesh().V()[particleCenterCellID], 1./3.);
-                vector coronaPoint = cellCentrePosition + (cellCentrePosition - minPeriodicParticlePos) * (corona / centreDist);
+                vector coronaPoint = cellCentrePosition;
+                if(centreDist > 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.);
+        scalar corona = 0.5 * sqrtThree_ * pow(particleCloud_.mesh().V()[neighbor], 1./3.);
-        vector coronaPoint = cellCentrePosition + (cellCentrePosition - position) * (corona / centreDist);
+        vector coronaPoint = cellCentrePosition;
+        if (centreDist > 0.0)
+            coronaPoint += (cellCentrePosition - position) * (corona / centreDist);
 
         if (!hashSett.found(neighbor) && pointInParticle(index, position, coronaPoint) < 0.0)
         {

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

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);
 

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;

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:
 

src/lagrangian/cfdemParticle/subModels/voidFractionModel/dividedVoidFraction/setWeightedSource.H

@@ -67,7 +67,7 @@
 
         }
 
-        particleWeights[index][storeInIndex] += 1/static_cast<scalar>(nPoints);
+        particleWeights[index][storeInIndex] += 1.0/nPoints;
         particleVolumes[index][storeInIndex] += particleVolume;
         particleV[index][0] += particleVolume;
         //====================================================//

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;