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