Does LIGGGHTS v3.8 support breakable bonds as a function of stress/force?

Submitted by dwhite on Fri, 08/24/2018 - 02:40

Hi,

I see that bonds that break as a function of stress/force are included in the premium version, however I am wondering if there is anything that is applicable in the public version. The documentation mentions the bond/create and bond/break commands, however these depend solely on distance. Furthermore, when attempting to use these commands we had numerous errors pop up throughout our simulations. After debugging we were able to trace the error back to the fix bond/create command, which was not recognized in our current version of LIGGGHTS. We tried enabling the molecular packages and recompiling LIGGGHTS with the MC package from LAMMPS as suggested in the documentation, but were unsuccessful in fixing the problem.

We already tried the LIGGGHTS-WITH-BONDS repo (https://github.com/richti83/LIGGGHTS-WITH-BONDS based on v3.3) discussed in other articles on this topic however it won't work for our particular application, as we need a contact model that is only available in LIGGGHTS v3.8 . We are trying to use bonds with deformable particles (Luding model). We tried porting over the relevant source code from LIGGGHTS-WITH-BONDs in to LIGGGHTS, as well as moving the Ludig model source code files from LIGGGHTS v3.8 into the LIGGGHTS-WITH-BONDS v3.3 both of which resulted in several compilation errors.

Any advice that could be provided would be much appreciated. Up until this point LIGGGHTS has been very useful for our application.

Thanks,
Doug

richti83's picture

richti83 | Fri, 08/24/2018 - 13:09

I stopped working on the bond problem as there is no funding for it, but mschramm got it to work with 3.7.0 (should be easy to push it to 3.8.0), start here:
https://www.cfdem.com/forums/bonded-particle-damping-extending-richti83s...
and copy the needed contact_* models to this src or the bond-files to public 3.8.0.

I'm not an associate of DCS GmbH and not a core developer of LIGGGHTS®
ResearchGate | Contact

dwhite | Fri, 08/24/2018 - 17:06

Thanks for pointing me towards the updated version. I managed to transfer over the Luding model and get it to compile. Do you know if this version includes all of the original functionality from LIGGGHTS-WITH-BONDS? This flexible bond fiber setup requires hard particles to be used and I am afraid this will lead to instability when used in combination with the deformable particles models. I will do some additional testing to see what I can get to work using this setup.

dwhite | Sat, 09/15/2018 - 00:58

Thanks so far for all the help. I would like to create a visualization of bonds that break as well. How do the IDs posted in the command:

compute b1 all property/local batom1x batom1y batom1z batom2x batom2y batom2z batom1 batom2 btype

relate to the atom IDs exported in the dump custom command shown below:

dump test all custom 100 dump.*.ids id type x y z

They should be the same correct?

dwhite | Fri, 09/28/2018 - 04:22

After getting the bond visualization working we realized one of the issues with our simulations was a lack of overall bond formation. We were able to adjust the parameters to create more bonds in our simulations, but after doing this we noticed some weird behavior.

Here are the basics of what we are trying to do (for some background):
1) We start by compressing a set of particles in to a cylinder. Particle-Particle interactions are modeled using the luding model. Particle wall interactions are modeled using the hertz model.
2) Once the compression is finished we create a set of bonds between the particles using the functionality provided by liggghts_flexible_fibers
3) Then, we perform uniaxial force testing by pressing the side of the cylinder with a moving mesh wall and recording the force

Steps 1 and 2 work well, and we can get robust formation of bonds between particles.

The main issue we are currently dealing with is that the particle are moving up the simulation domain after the bonds are formed. We suspect this has to do with instability in the bonds, but are unsure of where this could be coming from. We tried including dampening, but that made the particles more unstable and exit the simulation domain faster. Reducing the timestep (we tried 1E-6 and 5E-7) did not appear to fix the issue.

We are using the following parameters:
particle_radius = 6.6E-4
ro = 0.15 #outer bond radius
ri = 0 #inner bond radius
lb = 0.6 #particles are deformable, so a reduced bond length is required
bond_damp = 0.0
Y_bond = 1E8
G_bond = 3.4E7
time_step = 5E-6
break_stress = 5E7

If anyone could provide some guidance on how to optimize the parameters to achieve bond stability that would be much appreciated.

FelipeL | Mon, 12/24/2018 - 06:13

Hi dwhite, can you tell me how did you get the bond visualization? . I'm using the information provided by:

compute b1 all property/local batom1x batom1y batom1z batom2x batom2y batom2z batom1 batom2 btype

then i dump it but i cant read that information in paraview, even if them are in vtk format. Can you give me some guidance please? i Also need to visualize bonds that break. THanks in advance :)

Felipe

mschramm | Fri, 09/28/2018 - 18:29

Hello,
How did you determine your time step?
Using what you provided and Yu Guo's work on determining stable time steps, I found that your particle density would need to be around 1600 kg/m3.
It should be noted that adding damping drastically reduces the stable time step.
This is a hurdle in my research and most of my time has been in working with the damping of the bond.

I built the package with my research in mind which involves fibers (like wheat straw or corn stalks). You may need to substitute some of my code with what
was originally there from richti83 to get what you are looking for.

Lines that may need changed are
314: sndt = Kn * (r-bondLength)*rinv; --> sndt = Kn*dt;
315: fn_bond[0] = - sndt*delx; --> fn_bond[0] = -vn1*sndt;
316: fn_bond[1] = - sndt*dely; --> fn_bond[1] = -vn2*sndt;
317: fn_bond[2] = - sndt*delz; --> fn_bond[2] = -vn3*sndt;
at 363, add the following lines
// rotate normal force
rot = bondhistlist[n][0]*delx + bondhistlist[n][1]*dely + bondhistlist[n][2]*delz;
rot *= rsqinv;
vel_temp[0] = rot*delx;
vel_temp[1] = rot*dely;
vel_temp[2] = rot*delz;
vel_norm = sqrt (vel_temp[0]*vel_temp[0]+vel_temp[1]*vel_temp[1]+vel_temp[2]*vel_temp[2]);
f_norm = bondhistlist[n][0]*bondhistlist[n][0] + bondhistlist[n][1]*bondhistlist[n][1] + bondhistlist[n][2]*bondhistlist[n][2];
if (vel_norm == 0) f_norm =0;
else f_norm = sqrt (f_norm) /vel_norm;

bondhistlist[n][0] = f_norm*vel_temp[0];
bondhistlist[n][1] = f_norm*vel_temp[1];
bondhistlist[n][2] = f_norm*vel_temp[2];

410: bondhistlist[n][ 0] + fn_bond[0] --> bondhistlist[n][ 0] = bondhistlist[n][ 0] + fn_bond[0];
411: bondhistlist[n][ 1] + fn_bond[1] --> bondhistlist[n][ 1] = bondhistlist[n][ 1] + fn_bond[1];
412: bondhistlist[n][ 2] + fn_bond[2] --> bondhistlist[n][ 2] = bondhistlist[n][ 2] + fn_bond[2];

Then any mention below line 412 to 499 of fn_bond[k] needs to be replaced with bondhistlist[n][k]

dwhite | Mon, 10/01/2018 - 19:19

Thanks for your comment. Our current particle density is set as 1550 kg/m3. I will try incorporating the changes you suggested to see how this affects the stability.

Would you recommend removing the bond dampening all together (i.e. setting this parameter to 0)?

dwhite | Tue, 10/02/2018 - 17:15

Excellent. I started the process of swapping out the lines. You mentioned swapping these lines:

410: bondhistlist[n][ 0] + fn_bond[0] --> bondhistlist[n][ 0] = bondhistlist[n][ 0] + fn_bond[0];
411: bondhistlist[n][ 1] + fn_bond[1] --> bondhistlist[n][ 1] = bondhistlist[n][ 1] + fn_bond[1];
412: bondhistlist[n][ 2] + fn_bond[2] --> bondhistlist[n][ 2] = bondhistlist[n][ 2] + fn_bond[2];

After I added in the other lines at 363, I calculated these should correspond to lines 422 or so. However, when I look there I saw the following lines instead:
bondhistlist[n][ 0] = fn_bond[0]
bondhistlist[n][ 1] = fn_bond[1]
bondhistlist[n][ 2] = fn_bond[2]

These were the lines you referring to right? The plus sign above is just a typo and should be an equals sign?

dwhite | Wed, 10/03/2018 - 01:14

Thank you for the help. I made the changes to the code, recompiled successfully and ran another test. This appears to have fixed the issues with the bonded particle ensemble migrating up the simulation, and with the particle ensemble appearing to contract upon initiation of the bonds. I will run a few more tests, but for now this appears to be working - I really appreciate all of your help!

dwhite | Tue, 10/16/2018 - 17:06

After some testing we managed to get the simulations working as expected with respect to bond formation and breaking, thank you for all your previous help.

As mentioned above our application is for hardness testing after particle compression. Currently, the simulations lead to bonds breaking throughout defined regions in the tablet, but don't lead to the fissures we would expect. This results in a material that is extremely brittle, which displays a break pattern more like glass. After examining the force across the bonds, it appears that the force is oscillating somewhat during the hardness testing. We think this could be due to the dampening parameter, and the way this is currently implemented in LIGGGHTS.

Any thoughts on how we could go about defining this parameter correctly? In the modified version of the source code we are using, is there a maximum value that the dampening parameter can take?

dwhite | Wed, 10/17/2018 - 19:09

Thanks for sharing, that's very interesting. Is there a relationship between the necessary timestep and the dampening parameter?

dwhite | Wed, 10/24/2018 - 18:53

OK - after some more testing we have noticed that increasing the dampening parameter does seem to increase simulation stability.

However, we are still not observing fracture patterns that are consistent with our expectations. Now, the bonds either break at the edge of the simulations. We suspect there are a number or reasons why this could be occurring, one of which is related to the uniformity of bond lengths in the simulations.

Do you know if there is a way to capture different bond lengths between particles in LIGGGHTS?

dwhite | Thu, 10/25/2018 - 19:57

Thanks for the quick reply.

Let me check my understanding of how multiple bonds types would work. Let's say I have a bond type "type 1" with a bond length of 0.5 that I apply across all particle pairs within a distance of two times the particle radius. If I create a second bond "type 2" with bond length of 0.6 across all atoms within a distance of 2.2 times the particle radius, would this create particles which have two sets of bonds (type 2 and type 1)? Or can only one set of bonds be formed between each set of particles?

Is there any particular reason why LIGGGHTS doesn't support creating bonds which have individual equilibrium distances according to the initial locations of the atoms when the bonds are formed? How difficult would it be to modify the source code to accomplish this?

mschramm | Fri, 10/26/2018 - 15:16

You are correct but remember that each particle has a maximum number of bonds that it can have.
As to setting up the code to automatically get equilibrium distances... this was needed for my research.
It is not the main focus of my research but it was necessary to complete it.
Let me contact them today and see if I can release how to do the initial bond length calculation early. (I have made many changes to the bond model since last year...)
But it wasn't that hard to do (though I do question if I did it in the "best way").

If I can release it I will post step by step how to do it. This would make the parameter lb useless but leaving it in the code will not cause any problems.

mschramm | Fri, 10/26/2018 - 18:14

Make sure line 68 says
n_granhistory(13);
This is where we are going to put the eq length...

Delete line 261 and put
if (bondhistlist[n][12] == 0.0) # When a bond is created, all history values are set to 0.0
bondhistlist[n][12] = r;
// Set Bond Length. fabs is needed because we could get a negative bondlength.
bondLength = fabs(bondhistlist[n][12]);

I will also be making this change to the github soon.

dwhite | Wed, 10/31/2018 - 15:01

Thank you for all of your assistance so far. We incorporated these changes in to our branch of the script, and noticed that you made these changes to the Github master branch. While this did address some of the issues with bond formation and bond breaking, we still notice two main issues which we can't seem to get around:

1) When we start the simulation with an appropriately bonded tablet, once compression starts, all of the bonds that are parallel to the direction of compression disappear without changing the resultant force. This leaves a set of perpendicular bonds only which results in the tablet being much more compressible than it should be.
2) All of the bonds are still breaking primarily at the edge of the tablet. We checked the force across the particles, and while we can see changes in the forces when bonds break, it is unclear that the force is being transmitted across all of the particles correctly.

Any thoughts you could provide on this would be very helpful. This is probably our last round of attempts before we conclude that public version of LIGGGHTS is not useful for our application.

mschramm | Thu, 11/01/2018 - 02:50

Hello,
Sorry for the trouble this is causing you.
I must admit that all of my testing has only been for flexible fibers.
Could you post your test script and I will see what I can do over the weekend.
The only thing that is different from my implementation and the original work by Potyondya and Cundall
is how they did the bond damping. They used a force based design but looked to the
atom's velocities to determine the direction of damping and not their relative velocities.

Bias's picture

Bias | Thu, 11/01/2018 - 08:41

Hi,
I have been following this discussion a bit. I'm also interested looking into running Liggghts with breakable bonds as a function of stress/force.
Could anyone link to a pedagogical formulation of the bond model you are using? I have not seen that so far so have not understood the mathematical model of the bonds.
Does anyone have a working implementation for version 3.8.0 that is shared on GitHub? I would be happy to try it :)
I saw the article "Granular Shear Flows of Flexible Rod-like Particles" by Y. Guo. Is this bonds only working for rods i.e. parallel bonds or could it by bonds in multiple directions?

I would like to model snow as ice agglomerates bonded together that then will impact against a solid wall. I would like that the bonds grow depending on contact time but as a start I think I will try to just have bonds between the particles and not with the wall.

Best regards
Tobias

dwhite | Wed, 11/07/2018 - 20:00

We may have located the issue that is causing some of the instabilities I talked about in my previous post.

In the meantime - I have a question about how the bond deletion is carried out in LIGGGHTS. In the LIGGGHTS-flexible-fibers code (https://github.com/richti83/LIGGGHTS-WITH-BONDS/blob/master/src/bond_gra...), it appears that bonds are deleted by setting the bondlist[n][3] to 1. This occurs in multiple places throughout the code, one of which I have included below:

//2nd check if bond overlap the box-borders
if (x[i1][0]<(domain->boxlo[0]+cutoff)) {
bondlist[n][3]=1;
continue;

However in the compute local/property gran command source code which we use to output information on the total number of bonds (https://github.com/richti83/LIGGGHTS-WITH-BONDS/blob/master/src/compute_...) the following comment can be found:
/* ----------------------------------------------------------------------
count bonds on this proc
only count bond once if newton_bond is off
all atoms in interaction must be in group
all atoms in interaction must be known to proc
if bond is deleted (type = 0), do not count
if bond is turned off (type < 0), still count
------------------------------------------------------------------------- */

This seems to indicate that bonds are still counted if the type is not equal to 0, however in the bond_gran code above, they are set to 1.

Can you provide any comments on this to help me understand what the discrepancy is? It's important for our application that LIGGGHTS correctly output the number of bonds.

mschramm | Thu, 11/08/2018 - 04:48

Hello,
As far as I can tell, the number of bonds is remapped during neighbor finding during the integration step...
I have not done much with bond tracking but I have been using the numbond thermo command to count the number of bonds...

If this does not work for you, you could do something similar to compute_bonds but use neighbors. Something like
int nbondlist = neighbor->nbondlist;
int **bondlist = neighbor->bondlist;
tot_bonds = 0
for (k = 0; k < nbondlist; ++k) {
if (bondlist[k][3] == 1)
continue;
++tot_bonds;
}

dwhite | Thu, 11/08/2018 - 19:20

Thanks for the quick response. Can you elaborate on the usage of the numbond thermo command? I checked your source code and it looks like you use this as follows:

thermo_style custom step atoms numbond v_sim_time cpu cpuremain v_px2 v_pz2 v_vx2 v_vz2
thermo ${thermo_out}
thermo_modify lost ignore norm no
compute_modify thermo_temp dynamic yes
....
dump dmp all custom ${filestep} post/dump*.liggghts id type x y z vx vy vz fx fy fz omegax omegay omegaz radius

If I understand this correctly, this gets you the total number of bonds, but doesn't output the positions of the bonds.

We are currently getting the bond information using:

compute b1 all property/local batom1x batom1y batom1z batom2x batom2y batom2z batom1 batom2 btype
dump bnd all local ${output_ts} postbond/bonds*.bond c_b1[1] c_b1[2] c_b1[3] c_b1[4] c_b1[5] c_b1[6] c_b1[7] c_b1[8] c_b1[9]

Is the "compute all property/local" command supported with the changes you have made?

mschramm | Thu, 11/08/2018 - 19:59

Hello,
Yes, numbond only outputs the total number of bonds.

I have made some changes to compute_property_local to also output torques and eq_dist and should be fine.
I also had to make changes to the paraview bond reader that richti83 provided (I have forked this and have made the changes there).

I have been looking into this the past few days and the bond list that my code uses is rebuilt every time neighbors are calculated (if they need to be...).
So the check
if (x[i1][0]<(domain->boxlo[0]+cutoff)) {
bondlist[n][3]=1;
continue;

Is only there if the list isn't rebuilt before bond forces are recalculated.

The code in fix_bond_propagate_gran is where bonds are deleted for use in compute_property_local...

jonathancuba | Sun, 05/26/2019 - 23:28

Hello, I am working on quantifying the fragmentation of material by breaking bonds. I use liggghts-with-bonds, I would like to quantify how many bonds were broken to obtain the particle size of the rock after being fragmented by compression.
Could you tell me commands to give me information about the status of the bonuses, please?

compute b1 all property / local batom1x batom1y batom1z batom2x batom2y batom2z batom1 batom2 btype
dump bnd all local 100 posts / bonds * .bond c_b1 [1] c_b1 [2] c_b1 [3] c_b1 [4] c_b1 [5] c_b1 [6] c_b1 [7] c_b1 [8] c_b1 [9]
and
thermo numbond

Are there other commands that give information on the status of the bonds?