# Can a non-inertial Fast Lubrication Dynamics scheme be used in combination with pair_style granular ?

Dear all,

I am trying to build a simulation aiming to describe the behavior of a dense colloidal suspensions made of silica particles.

I defined a DLVO potential (pair_style yukawa/colloid+colloid) for what concerns the conservative inter-particle interactions and I chose a non-inertial Fast Lubrication Dynamics (nFLD) scheme (pair_style brownian and lubricateU) in order to describe the hydrodynamics interactions. In such configuration, since the net force on the particles is zero (we assume no inertia), as reported in LAMMPS manual (pair_style lubricateU), we solve the following linear system in order to get velocities and angular velocities (vector U) of the particles (see Bolintineanu, Dan S., et al. , 2014 for further details):

-R_{FU}(U-U^{infty})= -R_{FE}(E^{infty})-F^{rest},

where (in my case), F^{rest} accounts for forces and torques relative to the DLVO potential and the Brownian dynamics.

What I would like to do now, is to add also forces and torques deriving from effective particle-particle contact.
Thus, I thought to define also the pair_style granular in order to have Hertz contact forces +adhesion+rolling resistance etc. But here it comes my problem:

First of all, I saw that LAMMPS automatically generates an "overlap error" (defined in pair_style colloid) if the center-to-center distance b/w two particles (r)
becomes smaller of the contact distance (d =2r_p, where r_p is the particle radius. We have particle of the same type here).
But r<d is the precise condition to have contact and no contact forces can be taken into account if r<d is not allowed !

Second problem is relative to the damping terms of the contact forces.
If for example we consider the normal damping force, i.e. F_{n,damp}=-C_damp U_{n,rel}, we have a force that is function of particles velocity
since U_{n,rel} is the relative component of the normal velocity b/w the particles in contact. It follows that we need to modify the coefficient matrix
(R_FU written in the balance of forces at the beginning) in order to account for this term. This procedure needs then to be repeated for all the damping
terms included in the contact forces (i.e. the damping term of the tangential contact force, of the rolling resistance torque etc.)

Is this done in LAMMPS? Can I simply define the pair_style granular and the coefficient matrix of the hydrodynamics interactions is automatically
corrected ?

Thank you all in advance,
all the best

Gianluca

Gianluca,

A couple of points to your question.

“First of all, I saw that LAMMPS automatically generates an “overlap error” (defined in pair_style colloid) if the center-to-center distance b/w two particles ®
becomes smaller of the contact distance (d =2r_p, where r_p is the particle radius. We have particle of the same type here).
But r<d is the precise condition to have contact and no contact forces can be taken into account if r<d is not allowed !”

This is true, but it’s not too hard to alter the pair style to include the DLVO potential. Pair style colloid works best if you have not contacts, but they of course happen. So don’t feel too strongly attached to this pair style as a ground truth. You also have to be careful/mindful of how you will treat ‘sigma’ (interatomic distance), the inner cut-off length-scale for the lubrication interaction, and the length-scale of the granules (your d). Their interpretation should be consistent. I’ve had to do that in the past when studying cement paste rheology - used the inertial version.

“Second problem is relative to the damping terms of the contact forces.
If for example we consider the normal damping force, i.e. F_{n,damp}=-C_damp U_{n,rel}, we have a force that is function of particles velocity
since U_{n,rel} is the relative component of the normal velocity b/w the particles in contact. It follows that we need to modify the coefficient matrix
(R_FU written in the balance of forces at the beginning) in order to account for this term. This procedure needs then to be repeated for all the damping
terms included in the contact forces (i.e. the damping term of the tangential contact force, of the rolling resistance torque etc.)”

I don’t know what leads you to that conclusion. The balance “-R_{FU}(U-U^{infty})= -R_{FE}(E^{infty})-F^{rest}” still holds whether or not you write F^{rest} in terms of a Resistance matrix. Not only that but not writing as one also allows for cleaner pair styles in my opinion.

To your ultimate question though, on the lubricateU page under restrictions
“Currently, these pair styles assume that all other types of forces/torques on the particles have been already been computed when it is invoked. This requires this style to be defined as the last of the pair styles, and that no fixes apply additional constraint forces. One exception is the fix wall/colloid commands, which has an “fld” option to apply their wall forces correctly.”

Seems straightforward enough.

Eric,

thanks for the reply.

I didn’t get some of the points you made in your message:

-“This is true, but it’s not too hard to alter the pair style to include the DLVO potential. Pair style colloid works best if you have not contacts, but they of course happen. So don’t feel too strongly attached to this pair style as a ground truth. You also have to be careful/mindful of how you will treat ‘sigma’ (interatomic distance), the inner cut-off length-scale for the lubrication interaction, and the length-scale of the granules (your d). Their interpretation should be consistent. I’ve had to do that in the past when studying cement paste rheology - used the inertial version.”

My point is not about including the DLVO potential. This is already there.
What I want is to actually have contacts b/w the particles in order to account for the contact forces (Hertz, tangential, rolling friction…).

I was just wondering if this was already taken into account in LAMMPS, meaning then that the error one gets regarding particle overlapping when using the pair_style colloid should

not appear when the granular model is also in there. That’s because is the condition of particle overlapping that generates the contact forces computed in the granular style !

Sure, I can modify by myself the pair_style colloid and get rid of the error but then it kind of make me think that such dynamics (long-range and short-range interactions plus contact interactions, all together)
has not been taken into account. My point was ultimately related to this. If LAMMPS has been thought to consistently do this.

-"I don’t know what leads you to that conclusion. The balance “-R_{FU}(U-U^{infty})= -R_{FE}(E^{infty})-F^{rest}” still holds whether or not you write F^{rest} in terms of a Resistance matrix. Not only that but not writing as one also allows for cleaner pair styles in my opinion.

To your ultimate question though, on the lubricateU page under restrictions
“Currently, these pair styles assume that all other types of forces/torques on the particles have been already been computed when it is invoked. This requires this style to be defined as the last of the pair styles, and that no fixes apply additional constraint forces. One exception is the fix wall/colloid commands, which has an “fld” option to apply their wall forces correctly.”

Seems straightforward enough."

What lead me to that conclusion is that we have a balance of this kind:

R(U_n)=F1+F2(U_n),

where U_n is the vector velocity at timestep n, R the resistance matrix, F1 all the forces that do not depend on particle velocity and F2(U_n) the forces that depend on particle velocity at the same timestep.

If we do not correct the resistance matrix, then how is solved this system ?

Best,

Gianluca

Well, you can’t get it all for free, if no one before you has found a use for it - or at the very least submitted the code to LAMMPS for the pair style to do what you want it to do. That’s just the reality of open source.

I get what you’re saying in the R_lub*(U_n) = F1 + R_damp*(U_n) conundrum. Yes. There would be some error there, since the U_n is a result of and not an input to the integration. Again, I doubt that anyone has taken care of this. You could easily do it by modifying lubricateU - but it wouldn’t be general if you still have to use a conservative granular model. A combined granular/lubricateU might help you, but it would of course just be messier.

However, there are some ways out of this - which I think are very reasonable.

1. In general, the R_lub dominates by orders of magnitude - you should be able to verify this for yourself by looking at R_lub at inner cut-off vs. b_damp in the granular model. In this case, just turn the damping off as its not likely to change the results. I also cannot stress enough that LubricateU already assumes this model is heavily over-damped (no bouncing here), while granular models usually have little dissipation relatively speaking.

On top of this - the damping in the granular model is already typically unphysically small, because in general you use much softer particles than you would otherwise find in nature to speed up computations. In most (but not all*), computations this affects neither kinematics nor stress. Spherical particle approximations are usually much more restrictive when it comes to stress prediction in athermal granular slurries.

1. Do a study comparing the dynamic model w/o contact damping, the dynamic model w/ contact damping, and the overdamped lubricateU model w/o contact damping and confirm for yourself, whether or not the results make sense - and if contact damping is necessary. Perhaps you would have to do this anyways to please a reviewer.

If you do care about this discrepancy, its likely up to you to fix it for yourself.

Sure,
I am not looking for an open source perfectly suited for my needs !
I was just asking, in order to understand if someone else had already face this problem.

In any case, thanks for all your tips and suggestions!

Best,
Gianluca