lubricateU in oscillatory simple shear

Dear LAMMPS users,

I am trying to set up a simulation of a particle fluid suspension under oscillatory shear using the Fast Lubrication Dynamics model from the COLLOID package. The documentation has been pretty helpful thus far and I’ve been able to simulate the system using the lubricate pair style. I would like to add additional interaction models and would like to use the lubricateU pair style in order to make the simulations computationally more feasible.

I want to investigate the suspensions under oscillatory shear using the fix deform command. The lubricateU pair style has the argument gdot which is the strain rate. The documentation says that lubricateU adjusts the strain rate automatically using the fix deform (I use remap v) command. It also says that lubricateU needs to be the last pair style to be defined which is what I did.

However, for some reason my suspension only seems to be sheared in steady shear. Am I missing some definitions or using the wrong integrator? I copied part of my input script in the email below.

Thank you for your help!

It sounds like the problem is with the fix deform command. I suggest you first
simplify your model (simple pair style) and focus on getting the oscillatory deformation
to work as you want. You can try the fix deform xy wiggle option. You can also
the provided variable functions for swiggle and cwiggle to see if they work. Then try
your custom variables - if they don’t work you will have to debug what they are doing.


Hi everyone,

I had a look at the source code and it looks a lot like for time-dependent shear-rates the argument gdot for the pair style lubricateU is not recalculated every time-step. The strain rate is only calculated once upon initialization of the pair style. Consequently, for non-steady shear (such as sinusoidal excitations) modifications would have to be made? I am relatively new to LAMMPS, is there anyhow or anyone I could check back with who is more experienced using this particular pair style?

I added 3 lines to the source code that re-calculate gdot every time step based on the deformation of the domain. Now the results look like the particles are being sheared back and forth periodically, however, when I plot the shear stress as a function of time I get that the effective shear stress oscillates about a non-zero mean. This should not be the case, the shear stress should oscillate about zero mean for sinusoidal excitations.

I also noticed that the fix adapt doesn’t seem to work for the lubricate pair style. According to the documentation for fix adapt this fix should work for the pair style lubricate when called as

fix 1 all adapt 1 pair lubricate mu 1 1 v_VISC

I had a look at the source code and for the pair style lubricate there exists no extract(…) method. I guess this is why the fix adapt does not work? I get the error msg Fix adapt pair style param not supported.

I attached a bare-bones script to this email to that demonstrates some of the issues.

Is there anyone that has more experience with this pair style. Any help would be greatly appreciated.

Thank you,


in.2Dshear_oscill_hybrid_U (2.66 KB)