problems with ShirgaonkarIB forceModel with cfdemSolverIB solver?

Hello,
I am trying to characterize the cfdemSolverIB.
I adapted the twoSpheresGlowinskiMPI example, removing one sphere to have the simple one falling particle. All my units are consistent with SI system (in both DEM and CFD sides). I want to measure the z velocity of the sphere.
I started testing this using ArchimidesIB force as the only interaction. I used different densities and gravities values. All those cases work fine and they match to the theoretical prediction.
The problem is when I add the ShirgaonkarIB force (using a Kinematic viscosity nu [ 0 2 -1 0 0 0 0 ] 0.00002;). I run the same case and compare the velocity with the theoretical prediction given by

vel(time)=2./9.*(rho_particle-rho_fluid) * gravity * (radius*radius) /dynamic_viscosity * (1.-exp(-9/2.*dynamic_viscosity / (rho_particle * (r*radius) * time)))

(I am using the proper dynamic viscosity = rho_fluid*kinematic_viscosity).
these are the simulation parameters I am using g=9.81
gravity=9.81
rho_fluuid=100.
rho_sphere=130.
radius_1=0.0025
radius_2=0.002
radius_3=0.0045/2.
dynamic_viscosity=0.002
I am having half of the expected velocity. If I compare this theoretical equation with 2 times the particle velocity, it agrees nicely, even the exponential behaviour. As the gravity is the only variable that I can change without affect the the exponential behaviour the simulation looks like the value of gravity is taken wrongly (a half).

I am attaching a plot of this, I have the same behaviour for different particle radii. this just happens with cfdemSolverIB ( ShirgaonkarIB ).
Anyone that could help me with this?? I am really lost at this point about this problem.

Here it is my coupling properties file, just in case. thank you very much.

modelType none;
couplingInterval 100;
depth 0;
voidFractionModel IB;
locateModel engineIB;
meshMotionModel noMeshMotion;
regionModel allRegion;
dataExchangeModel twoWayMPI;
IOModel basicIO;
probeModel off;
averagingModel dilute;
clockModel off;
smoothingModel off;
forceModels
(
ShirgaonkarIB
ArchimedesIB
);
momCoupleModels
(
);
turbulenceModelType turbulenceProperties; //OFversion30x
//===========================================================================//
// sub-model properties
ShirgaonkarIBProps
{ velFieldName "U";
pressureFieldName "p";
}
ArchimedesIBProps
{
gravityFieldName "g";
voidfractionFieldName "voidfractionNext";
}
twoWayMPIProps
{
maxNumberOfParticles 10100;
liggghtsPath "../DEM/in.liggghts_run";
}
IBProps
{
maxCellsPerParticle 6000;
alphaMin 0.01;
scaleUpVol 1.0;
}
engineIBProps
{
treeSearch false;
zSplit 8;
xySplit 16;
}