Hello everybody,
I hope you are well. I want to report a bug within cfdemSolverIB solver and its imposing particle velocity step. As you may know, the particle velocity should be applied in the fluid domain in cfdemSolverIB according to this relation:
u = Ui + ωi × r
which u is the modified field velocity based on particles presence, Ui is the ith particle velocity from liggghts, ωi is the ith particle angular velocity from liggghts, and r which is the difference between cell centre occupied by particle and particle centre. The problem is with this vector. I should mention that this bug has been previously fixed by Kenneth Aycock (https://github.com/kenaycock/CFDEMcoupling-PUBLIC). When we have periodic boundary condition, the part of the rotating particle that is on the other side of the periodic box should take care of the vector. Indeed, in this case, the distance vector should pass the periodic side other than a straight line. So instead of this
https://github.com/CFDEMproject/CFDEMcoupling-PUBLIC/blob/b23726463721d5...
for(int subCell=0;subCell= 0)
{
// calc particle velocity
for(int i=0;i<3;i++) rVec[i]=U.mesh().C()[cellI][i]-position(index)[i];
for(int i=0;i<3;i++) angVel[i]=angularVelocities()[index][i];
velRot=angVel^rVec;
for(int i=0;i<3;i++) uParticle[i] = velocities()[index][i]+velRot[i];
// impose field velocity
U[cellI]=(1-voidfractions_[index][subCell])*uParticle+voidfractions_[index][subCell]*U[cellI];
}
}
We should check the periodicity of the box by:
for(int subCell=0;subCell= 0)
{
// calc particle velocity
vector positionI=position(index); //insert by Mo
if(checkPeriodicCells_) //insert by Mo
{ //insert by Mo
scalar r = radius(index);
vector cellPosition = U.mesh().C()[cellI];
// Expand r to make sure that we capture the boundary cells
r *= 1.5;
if(mag(cellPosition - positionI) > r) {
const boundBox& globalBb = U.mesh().bounds();
vector xMax = globalBb.max();
vector xMin = globalBb.min();
vector xRange = xMax - xMin;
for(int i=0; i<3; i++) {
if(positionI[i] > (xMax[i] - r) && cellPosition[i] < (xMin[i] + r) ) positionI[i] -= xRange[i];
else if (positionI[i] < (xMin[i] + r) && cellPosition[i] > (xMax[i] - r) ) positionI[i] += xRange[i];
}
}
} //insert by Mo
for(int i=0;i<3;i++) rVec[i]=U.mesh().C()[cellI][i]-positionI[i];
for(int i=0;i<3;i++) angVel[i]=angularVelocities()[index][i];
velRot=angVel^rVec;
for(int i=0;i<3;i++) uParticle[i] = velocities()[index][i]+velRot[i];
// impose field velocity
U[cellI]=(1-voidfractions_[index][subCell])*uParticle+voidfractions_[index][subCell]*U[cellI];
}
}
So, this can prevent the wrong velocity calculation for a "rotating" particle that is passing the periodic side.
I hope this bug can be fixed in the CFDEM-Coupling next versions.
Cheers,
Mo