Bonded Particle with Damping (Extending richti83's LIGGGHTS-With-Bonds)

Submitted by mschramm on Tue, 10/17/2017 - 20:56

I'm posting what I have done to incorporate damping into a flexible fiber model.
I am using Yu Guo's paper "Granular Shear Flows of Flexible Rod-like Particles"
This release currently works for LIGGGHTS 3.7.0
The only thing that i would like to explore further is this bond model interaction with the fix_template_multiplespheres command.
As this would allow the user to insert directly into the simulation rather than read in the particle and bond information from a text file.
I will also be uploading a FORTRAN code that allows the insertion of flexible fibers into a cylindrical region.

Here is a link to richti83's original bond code

Here is my link

felipe.merino | Wed, 10/18/2017 - 05:25

Thanks a lot for sharing this, i will try it on some tensile and flexure test i've been working on. Looks promising, as all my results before seem to be affected by the dissipation factor.
A question, what is the hard_particles line about? couldnt find anything related on documentation.
Again, thanks for sharing this.

mschramm | Wed, 10/18/2017 - 15:55

This is a new requirement for LIGGGHTS input scripts if you are using atoms with a young's modulus greater than a specific value (1e5 maybe).

Bond settings should be similar but I have added things for my specific research
bond_coeff 1 ${ro} ${ri} ${lb} ${Y_bond} ${G_bond} ${Bond_Damp} ${Break_Type} ...


ro is the outside radius of the "fiber bond" factor. This is a percentage of atom radius. So if your atom has radius 5 and your fiber bond needs the same radius, you would set ro to 1.

ri is the inside radius of the "fiber bond" factor. (see above)

lb is the bond length factor. This factor is based on the two atoms it is bonded with. Currently bond_length = lb*(radius(i1) + radius(2)).

Y_bond is the bond's young's modulus. The Kn factor is defined as Y_bond*A/bond_length

G_bond is the bond's shear modulus. The Kt factor is defined as G_bond*A/bond_length

Bond_Damp_Val is the damping value for the fiber. This has been fun... In testing, a value of 160 was needed to recreate lab tests. This was confirmed by Yu Guo.
Unfortunately, this required a very small timestep (1e-10) to run.

felipe.merino | Fri, 10/20/2017 - 04:39

I've used it a bit, and the first thing i noticed is a difference about the bond length. Did you redefined it? Im seeing that bonds created between particles that are not in contact, i.e. bond length is larger than the sum of particle radius, break instantly when they are created, a behavior that i never saw on the previous version. The force on the broken bond corresponds to a deformation of sqrt(2) times the bond length, which makes sense on my lattice distribution of particles. Maybe before the bond length was defined as the distance between centers of particles? I tried using length multiplier and it solves the problem, but i dont know if it will mess up the forces on regular bonds, making them to break at higher deformation than intended. Hope you understand my question, as english is not my native language.

mschramm | Fri, 10/20/2017 - 05:10

For what I have been using the bond model for. I defined a equilibrium bond length as is needed in Yu's paper to define Kn and Kt (Kn = Yb*Ab/Lb and Kt = Gb*Ab/Lb). Where Lb = lb*(radius(i1)+radius(i2)).

You may try using a smaller timestep during bond formation as this model can become VERY stiff.

Alberto | Tue, 10/24/2017 - 22:19

Hi, i have two question
1) In command Bond_coeff, the parameter ${Bond_Damp} have a value 10 in the example, but in the paper "Granular shear flow" it says this is in between 0 to 1

2) In the archive Bond_gran.cpp, when define the force of damping in the bond (line 338), there are a number two out of square, but in the paper it's inside the square

mschramm | Wed, 10/25/2017 - 15:33

I have been working with Yu to bring this out. He initially stated that beta needed to be between 0 and 1 but since then we have navigated away from that.
As for moving the 2 from inside to outside the sqrt, I do not know when he made that change, but that is what he is using now. A answer that I came up with
is that having the 2 where it is conforms to the damped harmonic oscilator
x'' + 2*beta*omega*x' + omega^2*x = 0 where omega^2 == kn/m and beta is the damping factor.
x'' + 2*beta*sqrt(kn/m)*x' + kn/m*x = 0
x'' + 2*beta*sqrt(kn*m)*x'/m + kn/m*x = 0
m*x'' = -2*beta*sqrt(kn*m)*x' - kn*x
F = F_bond_damp + F_bond
With this formulation, there exist no restraint on beta (as long as beta>0...). However, as beta increases the time step is driven further down.
This is where a lot of my research focus is...

Please let me know if you have any other questions

ceemrp | Thu, 11/02/2017 - 05:53

Hi mschramm,

Thank you so much for sharing this code :) I do have a simple question.

Does the bond fail exclusively under tensile stress? Will it fail under compression? I have taken a look at the source code (bond_gran.cpp), and found these lines [line 446 - 453],

nforce_mag = sqrt(fn_bond[0]*fn_bond[0] + fn_bond[1]*fn_bond[1] + fn_bond[2]*fn_bond[2]);
tforce_mag = sqrt(bondhistlist[n][3]*bondhistlist[n][3] + bondhistlist[n][ 4]*bondhistlist[n][ 4] + bondhistlist[n][ 5]*bondhistlist[n][ 5]);
ntorque_mag = sqrt(bondhistlist[n][6]*bondhistlist[n][6] + bondhistlist[n][ 7]*bondhistlist[n][ 7] + bondhistlist[n][ 8]*bondhistlist[n][ 8]);
ttorque_mag = sqrt(bondhistlist[n][9]*bondhistlist[n][9] + bondhistlist[n][10]*bondhistlist[n][10] + bondhistlist[n][11]*bondhistlist[n][11]);

nstress = sigma_break[type] < (nforce_mag/A + 2.*ttorque_mag/J*(rout-rin));
tstress = tau_break[type] < (tforce_mag/A + ntorque_mag/J*(rout-rin));
toohot = false;

which imply that the stress sign does not matter, and the bond fails if the stress magnitude is greater than the bond strength. I have not run your code to test this since I'm still struggling to compile it on Cygwin. In the meantime, I really appreciate if you can clarify this.


mschramm | Thu, 11/02/2017 - 15:07

The sign does not matter for this criteria. The code only checks if the magnitude of the normal force exceeds a specific value.

Are you using windows 10? If so, I would suggest using the Linux subshell option. This will give you an Ubuntu style experience.

achuth1992 | Tue, 05/08/2018 - 11:02

Hello mschramm,

Thank you for your excellent work. I have a doubt regarding the formulation of the normal force calculation.

sndt = Kn * (r-bondLength)*rinv;

Could you explain this part of the formulation? Or could you just point towards a source where I can clarify this part.

With regards

mschramm | Fri, 05/11/2018 - 07:51

Since we know what the bond length is, the true normal force can be calculated by
Fn = Kn*(r - bondlength)
You then need to get the x, y, and z components of Fn
Fn_x = Fn*dx/r
Fn_y = Fn*dy/r
Fn_z = Fn*dz/r

achuth1992 | Fri, 05/11/2018 - 10:23

Thank you for the fast reply. I am currently working with the bond formulation that you have provided and I am working on a non-linear normal stiffness formulation. I might have some questions later on and would like to know if it is alright if I contact you in the future to clarify them.

Best regards.

achuth1992 | Wed, 05/16/2018 - 10:40

I am trying to simulate two string of particles colliding with each other. I did two scenarios:
1) When damp = 0, the particles collide and keeps wobbling afterwards and does not arrive at an equilibrium state.
2) When damp = 50, the particles do not deform at all after the collision.

The function of the damp is to allow for the dissipation of the bond forces and moments. Still, in the case of damp = 50, there should be deformation right?

achuth1992 | Wed, 05/16/2018 - 14:57

Also, I would like to know if the terms 'dissipate' & 'damp' have some correlation. I have gone through the papers of Guo regarding the work performed for bond modelling. After analyzing, I came to the conclusion that the dissipation of the forces in the bonds (earlier meant to be performed with 'dissipate') was achieved using the velocity dependent damping.

After comparing your work and Richti's work, I saw that Richti has mentioned about a velocity-dependent damping that you were able to achieve using the damping forces ('damp').

Is my conclusion accurate?

Nathan | Wed, 06/27/2018 - 07:55

Hi mates,

I found these two papers below using the bond_model incorporated in LIGGGHTS to model natural fibres. However these papers address hydraulic conductivity, so they do not consider contact behaviour of fibres (no damping, no calibration for this aspect). Several validations are made with comparison to experimental tests.

I think the modified model made by Yu Guo will help improve those interested in modelling fibres in DEM. Contact and Collision of fibres are much different from hard particles.

Hi mschramm, I look at the paper by Yu Guo, but could not find any clear calibration. Do they compare the results with experimental measurement ? How could they estimate bond parameters: stiffness and damping, for example.


achuth1992 | Thu, 06/28/2018 - 13:54


I have played around with the bond model for quite a long time now and tried to simulate different scenarios with different conditions. In my opinion, the parameters of the bond can be scaled to the actual material parameters provided the bond radius and lengths are adjusted properly. Also, the damping will control the amount, rate and frequency of deformation.


Nathan | Tue, 07/10/2018 - 08:23

Hi Achuth,

Yes we can obtain parameters by trial and error scheme until the model matches the experimental data. I am confident with bond strength (e.g., tension), but for fibre contact it would be very tricky. Note that fibres are not hard particles and their damping ratio can be very high. I will try the damping bond model and compare with one without damping in my spare time.

Thanks for discussing.


achuth1992 | Wed, 07/11/2018 - 15:30

Hi Nathan,

I would suggest comparing the damping bond model with a spring model to find an analogy. I am not working with flexible fibres, but it involves deformation of the bonds. What I have observed is that the deformation with this bond model is plastic and will depend on the damping factor that we provide. For flexible fibres, this consideration works well.


achuth1992 | Sun, 11/18/2018 - 14:54


I have a question regarding periodic BC. When I tried a simulation with periodic BC, there was a loss of bonds across the boundaries. Has anyone come across a similar problem?

Setting a periodic BC will extensively reduce my simulation problem.


achuth1992 | Thu, 11/29/2018 - 16:25


I have solved the problem with the periodic boundaries. Include a check for the periodicity of the simulation box also. It is rather straight-forward.


mschramm | Fri, 11/30/2018 - 16:03

Thanks you,
could you post a quick line of how you did it or do a pull request?

I would be slightly careful about using this with periodic boundaries just because I do not know
how this model interacts with the rest of the code given a periodic boundary.
I simply did
x[i1][0]<(domain->boxlo[0]+cutoff) && domain->xperiodic == 0
x[i1][0]<(domain->boxlo[0]-cutoff) && domain->xperiodic == 0
and so on.

I will do some testing to see how this affects simulations.

achuth1992 | Mon, 12/03/2018 - 14:29


I basically did what you have mentioned. I also tested it with some DEM simulations with different combinations of boundary conditions and it seemed to work properly as far as I could see. I have included the modification below.

if ( x[i1][0] < (domain->boxlo[0]+cutoff) && !xperiodic) {
bondlist[n][3] = 1;
} else if (x[i1][0] > (domain->boxhi[0]-cutoff) && !xperiodic) {
bondlist[n][3] = 1;
} else if (x[i1][1] < (domain->boxlo[1]+cutoff) && !yperiodic) {
bondlist[n][3] = 1;
} else if (x[i1][1] > (domain->boxhi[1]-cutoff) && !yperiodic) {
bondlist[n][3] = 1;
} else if (x[i1][2] < (domain->boxlo[2]+cutoff) && !zperiodic) {
bondlist[n][3] = 1;
} else if (x[i1][2] > (domain->boxhi[2]-cutoff) && !zperiodic) {
bondlist[n][3] = 1;

if ( x[i2][0] < (domain->boxlo[0]+cutoff) && !xperiodic) {
bondlist[n][3] = 1;
} else if (x[i2][0] > (domain->boxhi[0]-cutoff) && !xperiodic) {
bondlist[n][3] = 1;
} else if (x[i2][1] < (domain->boxlo[1]+cutoff) && !yperiodic) {
bondlist[n][3] = 1;
} else if (x[i2][1] > (domain->boxhi[1]-cutoff) && !yperiodic) {
bondlist[n][3] = 1;
} else if (x[i2][2] < (domain->boxlo[2]+cutoff) && !zperiodic) {
bondlist[n][3] = 1;
} else if (x[i2][2] > (domain->boxhi[2]-cutoff) && !zperiodic) {
bondlist[n][3] = 1;


Min Zhang's picture

Min Zhang | Tue, 05/12/2020 - 22:04

Thanks for sharing!

I want to use the breakable bonds model to simulate a granular rock material. The bond has tensile and shear strengths.
Similar work is implemented in this paper through their in-house 2D DEM code:
Numerical modelling of backward front propagation in piping erosion by DEM-LBM coupling
D.K. Trana, N. Primeb, F. Froiioa∗, C. Callaric and E. Vincensa

Through some searching online, I found two bonded model codes so far:

Since your codes have been always updating and were initially based on richti83's original bond code, I want to try out yours. Could I ask some questions:
1. Can your code run in parallel?
2. Works for LIGGGHTS 3.7.0 only?
3. Does it have a requirement of the OpenFOAM version since I want to do CFD-DEM simulation?
4. Whether you could mention some pros and cons of your model? Or can you point out some doc I can refer to?
5. Do you have some tutorials for the modeling of a granular rock material?


mschramm | Mon, 05/25/2020 - 04:34

1. I have done some testing in parallel with the cantilever beam tutorial and found it to provide statistically similar results as to the single core results. This was repeated for the shear cell example and uniaxial compression example with similar results.

2. Code is currently working with 3.8 but has zero testing with superquadratics and the bond does not take into account this particle model.

3. The bond code itself, no but I recommend using 5.x. Simply because this is the latest version the CFDEM coupling software can use . My new research is into the drag model for a fiber particle so I myself am discovering any and all short comings. Currently, some users have reported reading in restart files can cause segmentation fault errors with mulisphere particles. I have done multiple tests with CFDEM with the fiber model and restarts and have not seen a similar result.

4. I have a doc of the bond model on the github. I should also have a powerpoint there as well. Biggest issue with the model is the damping model. Currently, it is recommended to use Yu Guo's velocity dependent damping, but this can cause very low time steps.

5. I do not.

Min Zhang's picture

Min Zhang | Mon, 05/25/2020 - 17:41

Question 2: So you are saying you don't consider the "bond model"? Or the bond model only works for LIGGGHTS 3.7.0?

I looked at Y. Guo's paper "Granular shear flows of flexible rod-like particles". I realized that there is a difference between this model and what I want.
(1) In my case, I have some singlesphese particles and there are bonds between these particles. What I want is like richti83's original bond code.
(2) In this paper, it has flexible rod-like particles/fibers and there is no connection between two fibers. There are bonds inside the fibers. I am wondering where these bonds can be broken or not.

Does your bond model have mechanical strengths?

So I would expect that I can define where the bonds are in the system. Then the above two cases should be consistent. Am I right?

Thanks again for all your kind help!!

mschramm | Mon, 05/25/2020 - 18:59

The code for bonds used in Yu Guo's paper and richti83 should be the same (I did some research with Yu Guo). The difference lies in how damping is used with Yu Guo (velocity based) and richti83 (no damping / artificial damping).

In the codes, the user states how bonds are formed by a distance multiplier. If two spheres are "close enough" then they form a bond. You can see this in my shear cell and uniaxial examples. I first insert a template of multiple spheres into the system then use a "bond_skin" small enough such that my templates form bonds while they do not form bonds with surrounding templates. If you make this "skin" larger, then templates would bond together.

Bonds can currently break by 1) being pulled to far apart, 2) exceeding a maximum tangential or normal stress, or 3) exceeding a max temperature.

Bonds are defined by inside and outside bond radii, bond young's modulus, bond shear modulus, damping coefficient, and break coefficients.

For bond formation, you can either use the bond_create fix to bond together spheres that are already in your system, or you can use the read_data command.

Min Zhang's picture

Min Zhang | Mon, 05/25/2020 - 19:22

Is that correct?

CFDEM®coupling version 3.8.0 released (released Dec 01, 2017)
by DCS computing GmbH

the release was developed for OF5.x and LIGGGHTS-3.8.0 (see versionInfo.H)
src and applications compile also with OF 3.0.x (compiling and execution tested!) - Note: syntax in tutorials is suited for OF5.x and must be changed if you run with OF3.0.x! Alice Hager (DCS Computing)
ability to compile with ParScale-PUBLIC (v1.2.1-beta) and LIGGGHTS 3.8.0 (minor changes in package-liggghts-list.txt) (DCS Computing)
check details in your (private) repository/doc directory, e.g. here.

Min Zhang's picture

Min Zhang | Wed, 05/27/2020 - 02:40


Except for the bond model part/feature, your liggghts should be the same as the official version, am I right? By "the same", I mean, when I am NOT using any feature related to the bond model, all the commands should be the same as the original official liggghts.

Thanks and best regards,

Min Zhang's picture

Min Zhang | Wed, 05/27/2020 - 20:11


If I find some errors/issues, where do you think is a good place to post? Here or your GitHub?

Thanks and best regards,

mschramm | Wed, 05/27/2020 - 21:43

If the question is related to how to use, I would say here. However, if it is an implementation issue or other code related issue, or if you find that the model should be altered, please create an issue on the github.

embastias | Tue, 06/30/2020 - 06:18

I hope you can help me I tried to install the LIGGGHTS Flexible Fibers version following the instructions in the pdf in the examples folder but I had some problems:

1. when I use the ln command to create the symbolic link I get this
ln: target '/usr/local/bin/liggghts' is not a directory

2 The command mpirun liggghts -in doesn't work but if I install the command liggghts I can use
Liggghts < in.chute_wear
to run this example and it works

3 But when I try the same with the example of the cantilever I get an error
ERROR: Invalid atom style (/build/liggghts-YO7u74/liggghts-3.8.0+repack1/src/atom.cpp:459)

This only happens with the bond/gran argument in the atom_style.

mschramm | Tue, 06/30/2020 - 06:49

Which distro are you using? If /usr/local/bin does not exist, you may use /usr/bin, just make sure you do not override an existing command.

What do you mean by "but if I install the command liggghts"? If you were not able to do (1) then (2) would also fail.

For 3) what is the compile date from the "liggghts" command that you used for the chute test from (2)? Does this compile date match the one in the src directory when running ./lmp_auto?

embastias | Tue, 06/30/2020 - 20:44

thank you for the quickly answer. I'm very new on linux system so your answer was very usefull.
i am using Ubuntu 18.04 LTS on wsl2, the both directories exist but I still getting this ln: target '/usr/local/bin/liggghts' is not a director

I achieve run the cantilever example with lmp_auto using this:

~/LIGGGHTS-Flexible-Fibers/src/lmp_auto < in.liggghts

on the cantilever example folder

The problem was the "liggght command" I was using, I followed a youtube tutorial that executes the command liggghts < in.chute_wear and the terminal suggested me to install liggghts (sudo apt install liggghts).

Do you have any idea how I can create the symbolic link?

embastias | Wed, 07/01/2020 - 00:36

I finally did it, I just needed to install the symlink. Thanks for the help.

mschramm | Sat, 07/04/2020 - 08:18

While working on a new project I have discovered an issue relating to cfd-dem simulations.
During the coupling process, the bond_hist variable in atom_vec_bond_gran gets reset.
This causes the equilibrium length to also be reset.

I am currently working on finding a fix for this but any help would be appreciated.