Hi,
I think that I have a few issues with the lubricate/poly implementation in LAMMPS based on my reading of Microhydrodynamics book by Kim and Karilla [1].
1) The gapdistance (h_sep) between the particles should be scaled by (radi+radj)/2 and not as radi, where radi, radj are the radii of the two particles.
2) The first term in the squeeze force seems to be missing a prefactor of 2.
3) \omega^\infty seems to have the wrong units. It is because "h_rate" has the units of length.
4) The pump term is also incorrect for particles of different sizes. Briefly, specific cases of calculation of torques in Ref. [1] cannot be used to write down a generalized version.
Please look at the attached pdf for a more detailed explanation on why I raised these concerns, and how to implement a "corrected" lubrication force if you agree with my concerns. Just to be clear, I have looked at previous messages in the mailing list before I send this message, and I don't think any of the previous messages have answered my concerns.
Regards,
Ranga
Forces_Torques.pdf (200 KB)
I’m CCing Dan Bolintineanu and Dave Heine who can likely
answer these Qs.
Steve
Ranga,
I didn’t get the attached pdf when the message was forwarded, so maybe sending it directly to me will help me understand the issue better.
In general, I followed the equations in Chapter 9 to incorporate polydispersity into pair_lubricate. As I was discussing with Tim Najuch, the text assumes you have particle A approaching another particle B, so being consistent with them, the separation distance is scaled by the radius of particle A. In the lubricate implementation, the forces on A and B are calculated separately, hence the requirement that “newton” is set to off. The grand resistance matrix approach in Chapter 11 that Tim was working on assumes the particles are approaching each other at the same speed, which may be a better approximation, but I don’t have a sense of how big the difference is when modeling things like highly filled systems as opposed to semidilute solutions. If you haven’t already talked to Tim about the grand resistance matrix implementation, maybe that will address some of your issues.
I don’t see the issues about specific terms you mention below, but again, maybe I need the pdf attachment to see your explanation. If you have a means of making this more generally applicable than what is provided in Kim and Karilla, then I am all in favor of it.
Best regards,
David
Hi David,
Thanks for your reply. I have attached the pdf as per your request. I will try to address your email in a enumerated list to make my points clear.

\omega^\infty seems to have the wrong units. It is because “h_rate” has the units of length/time (Please correct me if I am wrong). Specifically, I am referring to lines: “omega[i][0] += 0.5*h_rate[3]; …” which subtracts h_rate from omega. In case “h_rate” has the right units of 1/time, I am confused about the units of Ef (E^\infty) which should be the rate of strain tensor.

The results given in Chapter 9 of Ref. (1) are slightly misleading, because they cite Jeffrey and Onishi’s work (doi: 10.1017/S0022112084000355) before giving the final formulae for the forces, and torques. In the original work by Jefferey and Onishi, the gap distance is nondimensionalised by (radi+radj)/2.

Eqs. 9.26, and 9.27 in Ref. (1) are the solutions for force and torques for shearing of two surfaces only due to rotation. Why does the pump term account only for the torque? I don’t think the current formulation of the force considers the shearing of two surfaces only due to rotation correctly. (Please see Sec. IV of the attached document).

a.The squeeze term in lubricate/poly is taken from the force given in Eq. 9. 33 of Ref. (1). According to the resistance matrix formulation, the first term in Eq. 9. 33 should be multiplied by a prefactor of “2/(1+\beta)” , and the second term by “\beta” (apologies for the mistake in my previous email). One simple way to see that is that the magnitude of the leading order terms given in Sec. 11.2.2 should be twice of Eq. 9. 33, which is not the case in the textbook.

b. The squeeze or the shearing terms should be independent of the particle velocities and rotations, so Chapters 9 and 11 of Kim and Karilla should be consistent with each other. In case they are not, I have tried to refer to the original research articles and verify the same.

I don’t think we need to bring in volume fraction dependencies at the moment, because the issues that I have raised can be tracked down using just two particles of unequal sizes.
The general solution of the problem of two unequal spheres in a fluid is given in Chapter 11 of Kim and Karilla, or originally in Jeffrey’s research article (doi:10.1063/1.858494). In the attached pdf, I have mainly relied on Kim and Karilla as the reference. The results that you are referring to in Chapter 9 of Kim and Karilla can be derived as cases of the general result in Chapter 11 (as shown in Section IV of the attached pdf). As you rightly mention the grand (shear) resistance matrix formulation is slightly more involved to implement efficiently. However, one can get simplified expressions for forces and torques that are easier to implement in LAMMPS by considering only the first two leading order terms as shown in the attached document (Eqs. 12, or 13, and 23, or 24).
Regards,
Ranga
Forces_Torques.pdf (200 KB)
Ranga,
Yeah, it looks like you have done your homework. Tim submitted a GRM version of pair_lubricate based on Chapter 11, but maybe we are better off using the equations you highlighted here. Do you have this implemented in LAMMPS? It would be interesting to compare the behavior of two nearcontact spheres with the three versions of pair_lubricate we now have.
David
Hi David,
Thanks for your quick reply. Indeed, I have developed and tested a lubrication force class called lubricate/Simple, and it is currently a hack of the lubricate/poly. The link for my code is here:
I have took out most of the FLD parts of lubricate/poly since it is not needed for my purpose (dense nonBrownian suspensions). There are a few things that need to be included in my code as mentioned in the Readme. Hopefully, either you or some of the others in the LAMMPS community can help me with improving this code.
I have been in contact with Tim Najuch on the GRM version of the pair_lubricate as well. He has done a major revision of his code recently, and we expect to do a comparison of the results between our codes in a couple of weeks time. As for a quick comparison of the two particle results, I have done them for the analytical forms used in lubricate/poly in Sec IV D of my pdf (also in the bitbucket site).
Hope this helps.
Regards,
Ranga