Lubrication force implementation

Dear LAMMPS community,

I had a look at the code of pair_lubricate.cpp / pair_lubricate_poly.cpp
and would like to raise some points with regard to the lubrication force

Brief overview:

1) One additional "beta0" seems to be missing in 2nd squeeze term
2) Dimensionless gap distance definition different to literature
3) Inconsistency between "pair_lubricate_poly" and "pair_lubricate" for
"a_pu" term

In detail:

1) I think that there is a "beta" missing in the second term of the
squeeze term:

Current implementation (line 298 in pair_lubricate_poly.cpp):

"a_sq = ... +

Comparing it to Kim & Karilla "Microhydrodynamics: Principles and
Selected Applications", it should be multiplied by an additional "beta0":

"a_sq = ... +

2) Dimensionless gap distance "h_sep" in pair_lubricate_poly.cpp is:

"h_sep = (r - radi - radj)/radi"

In Kim & Karilla it is however:

"h_sep = 2*(r- radi - radj)/(radi+radj)"

This will affect the terms "log(1.0/h_sep)" and make them different to
the literature.

Was the former equation implemented deliberately like this and if so,
what is the reasoning behind it?

3) It seems that the leading order term for "a_pu" in pair_lubricate.cpp
is not consistent to "a_pu" in pair_lubricate_poly.cpp when monosized
spheres are chosen.
If we set "beta0 = 1" in pair_lubricate_poly, we obtain for the first
term in "a_pu = beta0*(4.0+beta0)/10.0/beta1/beta1*log(1.0/h_sep);"
(line 309):

"a_pu = 1/8*log(1.0/h_sep);"

Comparison to pair_lubricate.cpp (line 311):

"a_pu = ... *(3.0/160.0*log(1.0/h_sep));

I would have expected it to be the same.

Currently, I am having a look at all the resistance functions in Kim &
Karilla and it looks as if the implemented "a_pu" formula in
pair_lubricate_poly consists of the resistance functions "Y_11^B" and
"Y_12^B" instead of "Y_11^C" and "Y_12^C" (what I'd have expected).

I would be grateful if someone had any ideas to points 2) and 3).

Thank you in advance,


I’m CCing two people who worked on lubrication forces in

LAMMPS. Hopefully then can answer your Qs.


To add to this.

There was a post last summer that ran into similar issues with pair_lubricate_poly in combination with pair_style Brownian, that I fully embarrassed myself with trying to answer. (search - Mysterious attraction between hard spheres with pair lubricate) There he was unable to find the correct equilibrium behavior. From there, the conversation went offline on trying to correct it, and I never saw if they found the issue and resolved it.

Thanks for hunting these down Tim. I’m sure Jeremy and gang will be able to point you in the right direction.

Hi Tim,

Thanks for pointing this out. I’ve been looking into it; trying to recall details concerning the creation of lub_poly. May be we can work this out together. Here are some initial responses to your queries:

  1. The reference to Kim & Karrila you cite, are you referring to the equation in section 11.2.2 on page 273 (2005 version)? This should be compared and contrasted with equation 9.33 on page 231. I believe it’s equation 9.33 that was implemented in lub_poly. You’ll notice that the factor of beta in question does not show up in equation 9.33. I’m wondering if this is a typo in the book – I am currently looking into the differences. Can you confirm there is a problem with equation 9.33, or that it is not the correct form to use here?

  2. This also relates to the difference between of “epsilon” (Chapter 9) versus “xi” (Chapter 11). Our h_sep = epsilon. Notice, surficial, the definition of beta in chaper 11 mean 0 <= beta <= 1. This does not seem to be an issue in Chapter 9, but not clear why? Also, I think this form was chosen due to issues of the lub_poly not conserving momentum – but I’m not certain. If you interchange particle a and b, I don’t think the forces are not equal (I am trying to re-confirm/refresh my memory on this), some momentum should be imparted to the fluid… So, yes, I think this form was chosen for a reason.

  3. The choice, of course, depends on the type of relative motion that is meant to be captured. I believe the choice here is of Kim & Karrila Figure 9.2 and equation 9.26. Is that what you were expecting?


Thank you for your reply Jeremy.
Also, thank you Eric for chiming in with your experience using
lubricate_poly. I'll check also the old thread you mentioned.

To points 1-3:

1) Yes, I am referring to chapter 11 (eq. 11.2.2 or more generally to
the resistance functions given in chapter 11 (X_11^A and X_12^A for the
squeezing motion)). Thanks for pointing out that there are differences
to chapter 9, which I only skimmed because it looks as if chapter 9 is
only applicable when one sphere is stationary, whereas chapter 11 gives
all the resistance functions for the resistance matrix, i.e. both
spheres can move.
To be honest, it confuses me that the formula between chapter 9 and 11
is different.
I checked the paper by Jeffrey and Onishi (1984) "Calculation of the
resistance and mobility functions for two unequal rigid spheres in
low-Reynolds-number flow" which is more or less Kim&Karilla's chapter
11. It seems that there's also an additional "beta" in the paper. So,
chapter 11 gives apparently at least the same formulas as Jeffrey and

Personally, I'd tend to use the formulas found in chapter 11 because we
find this formulas also in the original paper by Jeffrey and Onishi + it
seems to be applicable to a more general case (both spheres can move and
not one stationary sphere).

2) I'll need additional time to think about the implications of having a
differently defined dimensionless gap distance.

3) Yes, it depends on the relative motion. In fig. 9.2, we have two
spheres without any translational motion and one sphere rotates.
I'd have expected eq. (9.27) for the torque calculation instead of eq
(9.26) which calculates the force (and not the torque).
If I stick to chapter 11, I can write for the torque acting on particle 1:

T_1 = mu*(Y_11^B * (U_1 - U) + Y_12^B * (U_2 - U) + Y_11^C * (Omega_1 -
Omega) + Y_12^C * (Omega_2 - Omega))

So, assuming that "a_pu" considers only the particle rotation (and not
translational motion) for the torque calculation, it would be the Y^C
terms which should contribute to the torque.
Currently, implemented is (9.26) (Y^B terms).
But, perhaps I am missing here something and there was a specific reason
for choosing the eq. (9.26) / Y^B terms.

Thank you again for taking your time!

Maybe it would help to contact the contributing author Dave Heine. I think he was responsible to taking the pair_lubricate that Amit Kumar from John Higdon's group had implemented in lammps and generalizing for polydisperse systems.


CCing Dave on this thread …


Thought I could add something substantive to this thread.

I altered some of my old input scripts for shearing cementitious particle slurries so that only granular and lubrication forces were present. The examples attached are stiffer than they need to be. I sheared systems of 45% and 55% volume fraction, using lubricate and lubricate/poly but both with monodispersed spheres. Latest stable release of LAMMPS was used Feb 16.

In the attached plots there are some differences in the pressure and shear stress response. The largest differences can be seen in the shear stress for the 55% vol. fraction - around 10% difference in the . Also there is some strange large fluctuations in the lubricate/poly with 45% volume fraction early on. Pressures show less difference, and the so called Shears stress ratio (shear stress/pressure) used in the granular community, probably shows the best agreement between lubricate/poly and lubricate.

In essence its hard to say if the two truly give different answers from what I’ve done, though the fluctuations in lubricate/poly for the same preparation procedure are strange.

Of course these were scripts I just had lying around. I’d suggest to Tim to look at two particle systems, with prescribed kinematics. That would be the ideal way to go for comparison.




Lubricatepolycompare.tar.gzip (928 KB)

Hi Eric,

Thanks for your input.

So far, without having had a look at your files: I am not surprised that
there's a difference because the implemented formulas in lubricate and
lubricate_poly are slightly different.
The lubricate class does not contain all of the terms found in
lubricate_poly. However, the dimensionless gap distance and the
coefficients for the implemented terms are the same in case of
monodisperse particles apart for the a_pu coefficient.

I'll have a look at your scripts and run them with modified lubricate /
lubricate_poly styles.

Best regards,



I pulled out my notes from 2010 and can at least settle the questions about the equations used in lubricate_poly. The equations did come from Chapter 9 instead of Chapter 11. Neither are perfect as in one case you assume one particle is at rest and in the other you assume the two particles are moving with the same velocity magnitude relative to the fluid. At the time, I probably looked at the two and said the Chapter 9 equations would be more computationally efficient even though I meant they would be easier for me to implement.

Regarding the pumping term, I believe there was some confusion between equation 9.27 and 9.25, not 9.26. 9.25 is also a torque, but doesn't apply to the rotation about the y axis. I have equations 9.24 and 9.25 highlighted in my book, so those must be the ones that were implemented. Regardless, it has to be fixed to give the proper pumping term.

I spent a little time plotting the forces in Chapter 9 and Chapter 11 and found that neither of them match the pair_lubricate force when the spheres are equal sized, although one is close at short separations and the other is close at long separations, so I'm not sure if pair_lubricate assumes one particle is stationary or they have the same magnitude or what.

Looking through what you posted and the Jeffrey paper, maybe following the Chapter 11 approach is generally closer to the typical particle interaction. I can work on implementing it when I get a chance, although I ran into a problem you may be able to explain. In both Chapter 11 and the Jeffrey paper, they include only the W^X term in the tabulated values that they say is a combination of L11^X(beta)+(1+beta)L12^X(beta)+betaL22^x(beta). The individual L11 and L12 terms are not given because "only this combination is needed when calculating the asymptotic behavior of the mobility functions". It looks to me like we do need the individual terms to calculate the resistance functions for the squeeze term. Is there a way around that, or would we have to re-derive the full expression for the squeeze term to determine L11 and L12?


Hi David,

Thank you for having a look at your old notes and clarifying which
equations were implemented.

You might want to wait a bit with implementing the lubrication forces
found in chapter 11 because I have already started to implement the
resistance matrix due to all the confusion. As soon as I am confident to
have it implemented correctly, I can send it to you for revision and
perhaps submission to the LAMMPS code itself if there is any interest.

Yes, I think also that we should ideally implement the L11 and L12 terms.
The thing is that I am not able to find a formula to compute the L11 and
L12 terms. It seems that those values have to be actually matched to
numerical data, as it is done by e.g. Jeffrey (1982)
"Low-Reynolds-number flow between converging spheres".
Also, it seems that the L11 and L12 terms have only a "stronger"
influence for large gap distances (see attached pdf in which I plot the
contributions of the different terms to the lubrication force).
I'd say that we can perhaps neglect those terms (depending on the
cut-off distance).

For the sake of completeness, I'd like to mention that the A11 and A12
terms can be apparently computed with formulas given by Jeffrey and
Onishi. However, it seems to be computationally relatively expensive so
that I would stick here to implementing the tabulated values. Also, I
think that the A11 and A12 terms are actually important to implement
(see plot in attached pdf).

Best regards,


LAMMPS_sq_lubricate_tim.pdf (179 KB)