LAMMPS WWW Site - LAMMPS Documentation - LAMMPS Mailing List Archives
Re: [lammps-users] Errors in lubricate/poly ?
[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]

Re: [lammps-users] Errors in lubricate/poly ?

From: Ranga Radhakrishnan <r.radhakrishnan@...652...>
Date: Thu, 13 Jul 2017 12:14:27 +0100


Sorry for another email addressing the same issue. Maybe I am beating a dead horse here, but another way to highlight the same error in omega calculation in both the lubricate and lubricate/poly pair styles would be to do a two particle shearing example as described earlier and in the attached input script (with one particle at the bottom edge and the other in contact) and increase the simulation box size. Even if one assumes that there should be a force between the particles, it should surely not depend on the box size. However, there is a box size dependency with these lubricate pair styles.

I also want to make it clear that the problem with omega is not the only problem with the current implementation of lubricate pair styles. Another interesting situation is when particles overlap when shearing (it can happen in simulations of dense suspensions paired with granular forces), where one has to be careful about how the lubrication force is calculated. There's also an issue with the non-symmetricity of the stress tensor due to lubrication force as currently calculated in LAMMPS. Specifically for lubricate/poly, there are other issues that I pointed out earlier.

I can discuss all these issues in more detail if the LAMMPS strategy is to modify the current implementation.

One clarification about one of my earlier comments. I quote from my earlier email.

"4) 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."

The first term between Eq. 9. 33 and Sec. 11.2.2 of Kim and Karrila is consistent, since they scaled the gap distance by particle radius in the first case, however the second term is still inconsistent between the two sections of the book. I think that the Chapter 11 is the correct version in the sense that the terms are consistent with Eq. (3.19a, b) of Jefferey and Onishi, J. Fluid Mech. (1984), wol. 139, pp. 261-290, which is the original research article.

I've also cc'ed the LAMMPS mailing list this time, apologies for not doing in my last emails.


On 11/07/17 18:56, Ranga Radhakrishnan wrote:

Hi David,

I was trying to illustrate a problem with the current lubricate implementation by using a specific example where we know what to expect. That is the reason that I chose to use fix nve/sphere and fix deform. According to LAMMPS documentation, the command "fix fid all deform 1 xy erate ${arg} remap v" will add a relative velocity (V=${arg}*box_length) only if the particles are across the periodic boundary. The way to get the particles that don't cross the boundary have a relative velocity with each other is to use a fix/nvt sslod, or another NEMD technique.

So for the case that I was describing (using fix nve/sphere), both the top particle at (0,1.53,0) and the bottom particle (0,0,0) will have the same velocity (zero relative velocity) even with fix deform, and hence my problem with the lubricate force in LAMMPS. As I was describing earlier, this problem seems to be less important when  "omega[i][0] += 0.5*h_rate[3]/domain->zprd;" etc. is corrected in pair_lubricate.cpp.

On the other hand, if the particles are on opposite sides of the cubic box (length=10) at say (0,0.5,0) and (0,9.53,0), they have a relative velocity V between them and again we should be able to match this to expected result. In the dump file we may not be able to see that these particle have a relative velocity between them (at least on time step=1).

Sorry for not ccing this discussion on the LAMMPS forum. Please let me know if I should.



On 11/07/17 16:39, Heine, David R wrote:



Can you clarify a point for me?  With particles separated in the y direction and shear in the xy direction, why do you say there is no relative velocity between the particles?  The top particle will have a larger x velocity than the bottom particle in proportion to the shear rate.  If you instead placed the second particle at (0,0,1.03), then they would have the same relative velocity in the presence of xy shear.




From: Ranga Radhakrishnan [mailto:r.radhakrishnan@...6990.....]
Sent: Thursday, July 06, 2017 12:54 PM
To: Heine, David R; Steve Plimpton; Bolintineanu, Dan Stefan (-EXP)
Subject: Re: [lammps-users] Errors in lubricate/poly ?


Hi David,

I think that the below example is sufficient to illustrate the importance of updating lubricate and lubricate/poly pair styles.

Consider two particles of diameter 1.0 placed at (0,0,0) and (0,1.03,0) in a cubic box of length=10 with velocity=0, and only interacting with lubricate (/poly) pair style whose outer cutoff is 1.05, and inner cutoff is 1.001 (flagfld=0). If the box is not sheared, then the force acting between the two particles is 0. If we used fix nve/sphere for integration and the box is deformed in the xy direction, the force should still be equal to 0 since there is 0 relative velocity between the particles. However, this is not the result that one obtains from lubricate, instead, one gets a force that depends on the box deformation rate. Please see the attached collision.lmp input file for other details of my simulation setup.

As proposed in my previous emails, I have developed a pair style that is based on the simplification of the Grand Resistance Matrix (GRM) given by Kim and Karilla's textbook. My latest document and code are on bitbucket ( From my last email, I have made a few changes to the calculation of stress in my force class, and to the Stokes drag force calculation. I have verified my code against Tim Najuch's "lubricate/GRM" formulation as well as some analytical results.

I hope this is useful.






On 01/07/17 17:04, Ranga Radhakrishnan wrote:

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 non-Brownian 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.





On 30/06/17 15:18, Heine, David R wrote:



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 near-contact spheres with the three versions of pair_lubricate we now have.




From: Ranga Radhakrishnan [mailto:r.radhakrishnan@...12...652...]
Sent: Wednesday, June 28, 2017 12:53 PM
To: Heine, David R; Steve Plimpton; Bolintineanu, Dan Stefan (-EXP)
Subject: Re: [lammps-users] Errors in lubricate/poly ?


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.

1) \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.

2) 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 non-dimensionalised by (radi+radj)/2.

3) 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).

4) 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.

4) 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.

5) 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).





On 28/06/17 14:51, Heine, David R wrote:



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 semi-dilute 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,




From: Steve Plimpton [mailto:sjplimp@...24...]
Sent: Wednesday, June 21, 2017 10:23 AM
To: Ranga Radhakrishnan; Heine, David R; Bolintineanu, Dan Stefan (-EXP)
Subject: Re: [lammps-users] Errors in lubricate/poly ?


I'm CCing Dan Bolintineanu and Dave Heine who can likely

answer these Qs.



On Wed, Jun 21, 2017 at 4:54 AM, Ranga Radhakrishnan <r.radhakrishnan@...652...> wrote:


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 gap-distance (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.



Check out the vibrant tech community on one of the world's most
engaging tech sites,!
lammps-users mailing list





variable DIA1 equal 1.0
variable MASS1 equal (PI*(${DIA1}^3.0)/6.0)
variable DIA2 equal 1.0
variable MASS2 equal (PI*(${DIA2}^3.0)/6.0)
variable FLD_VISC equal 1.0
variable T_STEP equal 0.00001
variable LEN equal 10

atom_style	sphere
boundary	p p p
newton		off
comm_modify	mode single vel yes

region		reg prism 0 ${LEN} 0 ${LEN} 0 ${LEN} 0 0 0 units box
create_box	2 reg

neigh_modify	delay 0

print "Vol=$(vol)"
create_atoms 1 single  0.0 1.43 0.0    
group right type 1
velocity right set 0 0 0 
set type 1 omega 0 0 0 

create_atoms 2 single 0.0 0.5 0.0    
group left type 2
velocity left set 0 0 0 
set type 2 omega 0 0 0

set type 1 diameter ${DIA1}
set type 1 mass     ${MASS1}
set type 2 diameter ${DIA2}
set type 2 mass     ${MASS2}

timestep ${T_STEP}
fix		1 all nve/sphere
pair_style      lubricate ${FLD_VISC} 1 0 1.001 1.05 1 0
pair_coeff * *
#pair_coeff      1 1 1.001 1.05
#pair_coeff      2 2 1.4014 1.47
#pair_coeff      1 2 1.2012 1.26

compute mytemp all temp
compute mypres all pressure NULL pair

fix 2 all deform 1 xy erate 1 remap v

thermo_style    custom step c_mytemp c_mypres[*]
thermo		1
thermo_modify	lost ignore norm no
compute_modify	thermo_temp dynamic yes

dump	id all custom 1 poly.lammpstrj id type radius mass x y z vx vy vz fx fy fz tqx tqy tqz omegax omegay omegaz
run 50
The University of Edinburgh is a charitable body, registered in
Scotland, with registration number SC005336.