cfdemSolverIB for viscoplastic non-Newtonian fluid

Hi,
I am working on cfdemSolverIB and I have obtained good results for multiple Newtonian cases. However, I want to implement the Herschel-Bulkley non-Newtonian fluid model in my fluid field, which I faced some noteworthy issues.
We know that in Newtonian cases for cfdemSolverIB, based on the Patankar IB method, the particle velocity imposing in the fluid velocity field is equal to considering the below source term in NS equation:
rho*d(ubar-u)/dt
which ubar is the fluid field after considering the particle velocity, and u is the velocity value before that. BTW, when we have viscoplastic non-Newtonian viscoplastic fluid, the viscosity (mu) in the calculation should be updated as well, because of the change in velocity and strain rate. But this is an important issue to know which part of the code is the best place for this change. Before imposing the particle velocity or after that? If after that, why we didn't do the same thing for Newtonian cases? In those cases, we ignored the effect of the particle in viscosity and just consider it in velocity correction.
I hope one of the guys from cfdem-Coupling can help me to clarify this issue. Especially, Dr Alice Hager who is an expert in the cfdemSolverIB part of the package.
Thank you for your great help in advance.
Cheers,
Mo