A bonded particle model

Submitted by zamir on Thu, 07/31/2014 - 19:04

Hi All,

I have coded up a bonded particle model based on the paper:
Obermayr, M., Dressler, K., Vrettos, C., Eberhard, P., "A bonded-particle model for cemented sand". Computers and Geotechnics. 01/2013; 49:299 – 313. DOI: 10.1016/j.compgeo.2012.09.001

The attached in-file shows how to use it. The bonds are defined in a read_data file.

In order to use this, you must compile liggghts with the molecule and asphere packages, if you don't you will get errors. There are 2 new operators that must be compiled with liggghts:
1) An integrator which keeps track of orientation, fix_nve_sphere_orientation.cpp/h
2) A bond formulation, bond_beam.cpp/h

Pros:
Works in parallel as far as I can tell. No new particle property is created that needs to be passed in parallel, only existing LAMMPS infrastructure is used.
No memory leaks yet!
No gimbal lock or breakdown of trig functions
Detects neighbor list rebuilds and adjusts itself
There is existing infrastructure in the code to implement a damping model - currently local damping is 0 (however, global damping will still act on the particles normally).

Cons:
The shrink style bounding box does not always work. I get consistent results by using a fixed bounding box (ie boundary f f f). I will investigate this further.
You have to define quite a lot of things in the read data file because the quaternion property belongs to ellipsoidal particles in LIGGGHTS.
The c++ code does not follow RAII, there is a lot of manual memory management, which opens the possibility for memory leaks.
Bond breakage has not been implemented (but its not too difficult if someone wants to try).

Please, please, I beg people to try it out. I really need some good understanding of where it works and where it doesn't, and if there are memory leaks.
Notes: all files are renamed to .txt, you must rename them before compiling and running LIGGGHTS. There are two ogv video files also attached renamed to txt.

Thats all for now,
Zamir

ckloss's picture

ckloss | Wed, 09/17/2014 - 10:43

Hi Zamir,

I'll try to have a short look at it at some point and will then come back to you!

Best wishes
Christoph

ckloss's picture

ckloss | Thu, 10/30/2014 - 14:30

Hi Zamir,

Thanks for sharing this - looks very interesting! We might contact you at some point in the near future (next months) to work towards incorporating this or a similar model into LIGGGHTS

Very best wishes
Christoph

zhuang | Wed, 07/08/2015 - 22:37

Hi, Zamir,

Thanks for your sharing! I'm trying to couple your codes in my work now. I want to know if you have considered the bond breakage now, which is also an important part of the BPM. BTW, I'd appreciate it if you could give any contact info for further communication.

Zhuang

ckloss's picture

ckloss | Wed, 07/15/2015 - 22:37

Hi Zhuang,

thanks for posting and sharing the model. Honestly, I did not have time to look into it - if I do, I will contact you

Thanks for sharing - I know there are many people interested in this modelling.

Best wishes
Christoph

richti83's picture

richti83 | Mon, 12/07/2015 - 08:00

I don't know the difference between the bond model of Prof. Eberhard et all and the model proposed by Potyondy and Cundall [1] (this model is used in commercial DEM Software PFC³D), but the second one can be found here:
https://github.com/richti83/LIGGGHTS-WITH-BONDS
It includes the possebility to break bonds under stress / tension or by a defined cutoff distance. It has been initial released by C. Kloss and I try to keep the source up to date with latest Liggghts (3.3.0 at this moment).
There are still some limitations:
- no damping of the bond-force
- memory expensive
- not fully tested with all fixes

[1] http://www.sciencedirect.com/science/article/pii/S1365160904002874

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

xlibp | Wed, 01/20/2016 - 10:41

hi, richti83

Many thanks for your sharing.
One question: Does the limitation "no damping of the bond-force" means the information of bonds cannot be written to a restart file?

The reason why asking this:
Step1: I generate a particle system with no bonds formed and then write the particle sample information to a restart file.
Step2: Call the restart file in further simulation.
Everything is fine.

But if there are bonds formed in Step 1, then in Step2, it always collapse with the following error:

LIGGGHTS (Version LIGGGHTS-PUBLIC 3.3.0, compiled 2016-01-20-14:44:37 by root, git commit based on LAMMPS 23 Nov 2013)
Reading restart file ...
orthogonal box = (0 0 0) to (4 2 0.4)
1 by 1 by 1 MPI processor grid
WARNING: Bond granular: This is a beta version - be careful! (../bond_gran.cpp:83)
[lxy-OptiPlex-990:04315] *** Process received signal ***
[lxy-OptiPlex-990:04315] Signal: Segmentation fault (11)
[lxy-OptiPlex-990:04315] Signal code: Address not mapped (1)
[lxy-OptiPlex-990:04315] Failing at address: (nil)
[lxy-OptiPlex-990:04315] [ 0] /lib/x86_64-linux-gnu/libpthread.so.0(+0x10340) [0x7fcb5f231340]
[lxy-OptiPlex-990:04315] [ 1] /home/lxy/LIGGGHTS/LIGGGHTS-PUBLIC/src/lmp_fedora() [0x4282d9]
[lxy-OptiPlex-990:04315] [ 2] /home/lxy/LIGGGHTS/LIGGGHTS-PUBLIC/src/lmp_fedora() [0x43cd64]
[lxy-OptiPlex-990:04315] [ 3] /home/lxy/LIGGGHTS/LIGGGHTS-PUBLIC/src/lmp_fedora() [0x9d3ef5]
[lxy-OptiPlex-990:04315] [ 4] /home/lxy/LIGGGHTS/LIGGGHTS-PUBLIC/src/lmp_fedora() [0x91f7f6]
[lxy-OptiPlex-990:04315] [ 5] /home/lxy/LIGGGHTS/LIGGGHTS-PUBLIC/src/lmp_fedora() [0x91d47f]
[lxy-OptiPlex-990:04315] [ 6] /home/lxy/LIGGGHTS/LIGGGHTS-PUBLIC/src/lmp_fedora() [0x91dfb6]
[lxy-OptiPlex-990:04315] [ 7] /home/lxy/LIGGGHTS/LIGGGHTS-PUBLIC/src/lmp_fedora() [0x40e016]
[lxy-OptiPlex-990:04315] [ 8] /lib/x86_64-linux-gnu/libc.so.6(__libc_start_main+0xf5) [0x7fcb5ee7dec5]
[lxy-OptiPlex-990:04315] [ 9] /home/lxy/LIGGGHTS/LIGGGHTS-PUBLIC/src/lmp_fedora() [0x40ef3f]
[lxy-OptiPlex-990:04315] *** End of error message ***
Segmentation fault (core dumped)

Here is the code:
Step1: in.bonds1

atom_style hybrid granular bond/gran n_bondtypes 1 bonds_per_atom 20
atom_modify map array
boundary f f f
newton off
communicate single vel yes
units si
region reg block 0 4 0 2 0 0.4 units box
create_box 2 reg
neighbor 0.003 bin
neigh_modify delay 0
#New pair style
pair_style gran model hertz tangential history rolling_friction cdt #Hertzian without cohesion with rolling friction
pair_coeff * *
#bond
bond_style gran
variable simplebreak equal 0
variable stressbreak equal 1
#for simplebreak
#bond_coeff 1 0.0025 10000000000 10000000000 ${simplebreak} 0.002501
#for stressbreak
bond_coeff 1 1 1e8 1e8 ${stressbreak} 1e4 1e4
#added notes bond_type bond-radius_multiplier normal_stiffness shear_stiffness tensile_strength shear_strength
#walls
fix xwall1 all wall/gran model hertz tangential history rolling_friction cdt primitive type 2 xplane 0
fix xwall2 all wall/gran model hertz tangential history rolling_friction cdt primitive type 2 xplane 4
fix zwall1 all wall/gran model hertz tangential history rolling_friction cdt primitive type 2 zplane 0
fix zwall2 all wall/gran model hertz tangential history rolling_friction cdt primitive type 2 zplane 0.3
timestep 0.0000005
fix gravi all gravity 9.81 vector 0.0 1.0 0.0
# create particles IN A DEFINED REGION
region bc block 0 4 0 0.05 0 0.4 units box
lattice custom 0.05 a1 1.0 0.0 0.0 a2 0.0 1.0 0.0 a3 0.0 0.0 1.0 basis 0.5 0.5 0.5
create_atoms 1 region bc
group bb region bc
set group all density 2500 diameter 0.05
#set group all vz 0 vx 0
#Material properties required for new pair styles
fix m1 all property/global youngsModulus peratomtype 7e10 7e11
#E=70GPa units Pa
fix m2 all property/global poissonsRatio peratomtype 0.3 0.3
fix m3 all property/global coefficientRestitution peratomtypepair 2 0.3 0.3 0.3 0.3
fix m4 all property/global coefficientFriction peratomtypepair 2 0.7 0.7 0.7 0.7
fix m5 all property/global coefficientRollingFriction peratomtypepair 2 0 0 0 0
group nve_group region reg
#mass 1 1.0 #dummy
mass * 1.0
##############################################==============================================================
#fix bondcr all bond/create/gran 1 1 1 0.041 1 20 ###########NO bonds formed and NO Error occurs when running in.bonds2
fix bondcr all bond/create/gran 1 1 1 0.051 1 20 ############Bonds formed and Error occurs when running in.bonds2
##############################################==============================================================
#-import mesh from cad:
fix bottom all mesh/surface/stress file bottom.stl type 2
fix wall all wall/gran model hertz tangential history rolling_friction cdt mesh n_meshes 1 meshes bottom
#apply nve integration to all particles that are inserted as single particles
fix integr all nve/sphere
#output settings, include total thermal energyg
fix ts all check/timestep/gran 10000 0.9 0.9
compute 1 all erotate/sphere
variable etotal equal ke+c_1
thermo_style custom step atoms numbond ke c_1 vol
#thermo_style custom step atoms ke c_1 vol
thermo 10000
thermo_modify lost ignore norm no
compute_modify thermo_temp dynamic yes
#insert the first particles so that dump is not empty
run 1
#fix_modify bondcr every 0
dump dmp all custom 10000 post1/dump*.liggghts id type x y z vx vy vz fx fy fz omegax omegay omegaz radius
dump dumpstl all mesh/stl 10000 post1/dump*.stl
run 30000
write_restart liggghts.restart1

Step2: in.bonds2

#bond
atom_style hybrid granular bond/gran n_bondtypes 1 bonds_per_atom 20
#atom_style granular
atom_modify map array
boundary f f f
newton off
communicate single vel yes
units si
#read the restart file
read_restart liggghts.restart1
neighbor 0.003 bin
neigh_modify delay 0
#New pair style
pair_style gran model hertz tangential history rolling_friction cdt #Hertzian without cohesion with rolling friction
pair_coeff * *
#bond
bond_style gran
variable simplebreak equal 0
variable stressbreak equal 1
#for simplebreak
#bond_coeff 1 0.0025 10000000000 10000000000 ${simplebreak} 0.002501
#for stressbreak
bond_coeff 1 1 1e8 1e8 ${stressbreak} 1e4 1e4
#walls
fix xwall1 all wall/gran model hertz tangential history rolling_friction cdt primitive type 2 xplane 0
fix xwall2 all wall/gran model hertz tangential history rolling_friction cdt primitive type 2 xplane 4
fix zwall1 all wall/gran model hertz tangential history rolling_friction cdt primitive type 2 zplane 0
fix zwall2 all wall/gran model hertz tangential history rolling_friction cdt primitive type 2 zplane 0.4
reset_timestep 0
timestep 0.0000005
fix gravi all gravity 9.81 vector 0.0 1.0 0.0
#Material properties required for new pair styles
fix m1 all property/global youngsModulus peratomtype 7e10 7e11
#E=70GPa units Pa
fix m2 all property/global poissonsRatio peratomtype 0.3 0.3
fix m3 all property/global coefficientRestitution peratomtypepair 2 0.3 0.3 0.3 0.3
fix m4 all property/global coefficientFriction peratomtypepair 2 0.7 0.7 0.7 0.7
fix m5 all property/global coefficientRollingFriction peratomtypepair 2 0 0 0 0
#mass 1 1.0 #dummy
mass * 1.0
#-import mesh from cad:
fix bottom all mesh/surface/stress file bottom.stl type 2
fix wall all wall/gran model hertz tangential history rolling_friction cdt mesh n_meshes 1 meshes bottom
#apply nve integration to all particles that are inserted as single particles
fix integr all nve/sphere
#output settings, include total thermal energyg
fix ts all check/timestep/gran 100000 0.9 0.9
compute 1 all erotate/sphere
variable etotal equal ke+c_1
thermo_style custom step atoms numbond ke c_1 vol
thermo 10000
thermo_modify lost ignore norm no
compute_modify thermo_temp dynamic yes
#insert the first particles so that dump is not empty
run 1
dump dmp all custom 10000 post2/dump*.liggghts id type x y z vx vy vz fx fy fz omegax omegay omegaz radius
dump dumpstl all mesh/stl 10000 post2/dump*.stl
run 100000

Any comments are appreciated.

Best regards.

Maryam | Sat, 10/15/2016 - 13:28

Hi Richti,

Thanks a lot for sharing the bonded particle model that you've developed. I have tried using it, mostly with success, a few times. The only problem I get is that sometimes when running on parallel it randomly terminates the run, while the same case runs perfectly in serial.

So I'm wondering if this is because this model is not fully tested, or I might not be aware of specific tricks.

Your reply is greatly appreciated!

Maryam

zamir | Tue, 12/08/2015 - 22:53

Hi Richti83,

I love your bonded model. I chose the obermayr model because quaternions were readily available in liggghts and I could see a way to parallelize the model. Do you know of a good way to parallelize the Potyondy and Cundall model?

I think that a lot of researchers would very much benefit from a parallel bonded model.

brucef | Sun, 01/10/2016 - 20:34

Hi Richti83,

Your version of bonded particle model (BPM) works beautifully! Thank you for sharing this!
I just have two questions:

1- Is there currently an option in your implementation for NO bond breakage after initial bond formation?
2- Is it possible to form bonds between particles and granular wall?

Thank you again.

I am using this for modeling rock and fault friction simulations. I am also going to look at the influence of "neigh_modify exclude" of bonded particles for this simulation purposes and will report it back here.

Best regards,
Behrooz (Bruce) Ferdowsi, Dr. sc. ETH Zurich
Synthesis Postdoc Fellow and Postdoctoral researcher
National Center for Earth-Surface Dynamics (NCED) and
Earth and Environmental Science, University of Pennsylvania
58A Hayden Hall, 240 South 33rd Street
Philadelphia, PA 19104-6316, USA
behroozf.github.io

Behrooz (Bruce) Ferdowsi, Dr. sc. ETH Zurich
Hess Fellow and Postdoctoral Research Associate
Department of Geosciences, Princeton University
Room 314, Guyot Hall, Princeton, NJ 08544
behroozf.github.io

richti83's picture

richti83 | Mon, 01/11/2016 - 08:25


1- Is there currently an option in your implementation for NO bond breakage after initial bond formation?

Just set nstrengt / ssthrength or bond-break radius very very high. Or comment line 352 and 375 in bond_gran.cpp.

2- Is it possible to form bonds between particles and granular wall?

no, but you could make a wall of freezed particles

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

VBaric | Fri, 08/05/2016 - 15:21

Hi Richti83,
hi everyone,

I have just recently started working on DEM Simulations with ligggths using the BPM for rigid bonds between particles. It has worked quite nicely for me in the beginning, but after looking more into detail I have found an Issue that causes me headache. I have also found the same issue as xlibp (http://www.cfdem.com/comment/16298#comment-16298) and I have furthermore issues with the mpirun:

1. When I just simply pull 2 particles apart which are bond, I get a periodic breakup of the same bond. So the particles move away from each other, a force builds up until the bond breaks. In order to detect that I have uncommented line 334 in bond_gran.cpp. After that the force is 0. A few timesteps later the force builds up again! However, I do not get the message that a new bond is formed. This is periodically repeated. I have found out that the distance between the breakups equals have the distance defined for the neighborslist (after changing this, the breakup distance also changes). This goes on forever.

2. Reading a restart file leads to an error as stated by xlibp.

3. Starting simulations with mpirun -np X lmp_liggghts < in.liggghts <\em> leads to an error (inputfiles see below):

Step Atoms KinEng CPU
1 2 0 0
creating bond btw atoms 0 and 1 (i has now 1 bonds) at step 2
[hmipc17:13819] *** Process received signal ***
[hmipc17:13819] Signal: Segmentation fault (11)
[hmipc17:13819] Signal code: Address not mapped (1)
[hmipc17:13819] Failing at address: (nil)
[hmipc17:13819] [ 0] /lib/x86_64-linux-gnu/libpthread.so.0(+0x10330) [0x7fbe48aec330]
[hmipc17:13819] [ 1] /home/baric/bin/liggghts_3.4.1_dev(_ZN9LAMMPS_NS17FixBondCreateGran14post_integrateEv+0x961) [0x81c1c1]
[hmipc17:13819] [ 2] /home/baric/bin/liggghts_3.4.1_dev(_ZN9LAMMPS_NS6Modify14post_integrateEv+0xfc) [0x1319b6c]
[hmipc17:13819] [ 3] /home/baric/bin/liggghts_3.4.1_dev(_ZN9LAMMPS_NS6Verlet3runEi+0xf9) [0x1467829]
[hmipc17:13819] [ 4] /home/baric/bin/liggghts_3.4.1_dev(_ZN9LAMMPS_NS3Run7commandEiPPc+0xc6f) [0x13d29af]
[hmipc17:13819] [ 5] /home/baric/bin/liggghts_3.4.1_dev(_ZN9LAMMPS_NS5Input15command_creatorINS_3RunEEEvPNS_6LAMMPSEiPPc+0x29) [0x12e5f39]
[hmipc17:13819] [ 6] /home/baric/bin/liggghts_3.4.1_dev(_ZN9LAMMPS_NS5Input15execute_commandEv+0x243b) [0x12e0d3b]
[hmipc17:13819] [ 7] /home/baric/bin/liggghts_3.4.1_dev(_ZN9LAMMPS_NS5Input4fileEv+0x1dd) [0x12dcb3d]
[hmipc17:13819] [ 8] /home/baric/bin/liggghts_3.4.1_dev(main+0x94) [0x12fe854]
[hmipc17:13819] [ 9] /lib/x86_64-linux-gnu/libc.so.6(__libc_start_main+0xf5) [0x7fbe48738f45]
[hmipc17:13819] [10] /home/baric/bin/liggghts_3.4.1_dev() [0x71ed19]
[hmipc17:13819] *** End of error message ***
--------------------------------------------------------------------------
mpirun noticed that process rank 1 with PID 13819 on node hmipc17 exited on signal 11 (Segmentation fault).
--------------------------------------------------------------------------
<\code>

Here you find my inputfiles:
in.liggghts:

atom_style hybrid granular bond/gran n_bondtypes 1 bonds_per_atom 30
units si
boundary p p p
newton off
read_data data.2Particles
communicate single vel yes
# Neighborslist cutoff
neighbor 3.0e-9 bin
neigh_modify delay 0
# Interaction
fix m1 all property/global youngsModulus peratomtype 5.6e10
fix m2 all property/global poissonsRatio peratomtype 0.28
fix m3 all property/global coefficientRestitution peratomtypepair 1 0.94
fix m4 all property/global coefficientRollingFriction peratomtypepair 1 2.5e8
fix m5 all property/global coefficientFriction peratomtypepair 1 6.5e8
fix m6 all property/global cohesionEnergyDensity peratomtypepair 1 2e8
pair_style gran model hertz tangential history rolling_friction cdt
pair_coeff * *
# bonds
bond_style gran
bond_coeff 1 2.400000e-09 2.800000e+18 1.093750e+18 1 1.40e+09 1.40e+09
# Initialize
timestep 5e-13
run 1
# Create Bond
fix bond1 all bond/create/gran 1 1 1 1.010000e-08 1 9
run 1
# disable bond creations for further simulation
fix_modify bond1 every 0
# Movement
group pull id 1
group freeze id 2
fix F_frozen1 pull setforce 0.0 0.0 0.0
fix F_frozen2 freeze setforce 0.0 0.0 0.0
fix F_move pull move linear 0.0 -0.5 0.0 units box
fix F_force all ave/time 1 500 500 f_F_frozen1[1] f_F_frozen1[2] f_F_frozen1[3] f_F_frozen2[1] f_F_frozen2[2] f_F_frozen2[3] file forces.out
thermo 500
dump dmp all custom 500 trajectory.trj id mol type radius x y z ix iy iz vx vy vz fx fy fz omegax omegay omegaz
run 100000
<\code>

and data.2Particles:

LIGGGHTS data file from DLA Simulation
2 atoms
1 atom types
1 bond types
0 angle types
24 extra bond per atom
-5.0e-07 5.0e-07 xlo xhi
-5.0e-07 5.0e-07 ylo yhi
-5.0e-07 5.0e-07 zlo zhi
Masses
1 2.315992e-21
Atoms
1 1 0.000e+00 0.000e+00 5.000e-09 1.015e-08 4.230e+03 1
2 1 0.000e+00 1.000e-08 5.000e-09 1.015e-08 4.230e+03 1
<\code>

I hope you can help me face these issues!
Thank you
VBaric

VBaric | Fri, 08/05/2016 - 15:24

P.S.: I have used the newest versions of both LIGGGHTS and the bond model

VBaric | Thu, 08/11/2016 - 15:37

Hi,

I have found that there was indeed a version conflict responsible for the erroneous periodical bond breakup. So I have fixed that.
But still the mpirun does not work and also the restarts do not work. Did someone experience the same?

Regards
VBaric

richti83's picture

richti83 | Mon, 08/15/2016 - 15:00

The order of commands is important, you need first specify the bond_style gran, than read the datafile and than set the bond properties.
(Otherwise no bond history is created and the create_bond command fails).

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

NTT1508 | Fri, 08/12/2016 - 08:09

Dear Zami,
the bond model by Potyondy and Cundall is more simple than the one by Obermayr which might be able to capture the elasto-plastic behaviour and ductile breakage. I would like to compare in details by applying these models for particular cases.

Zami, Could you kindly share your codes about the bond model proposed by Obermayr et al again ? Your previous links have died.

Thank you,

Regards,

Ermek Asylbekov | Wed, 06/21/2017 - 16:15

Hi everyone,
I'm also using the bonded particles model that was proposed by richti83. It works great even in CFDEM! However if I use periodic boundary all particles that exceed the boundary and reappear on the other side of my geometry lose all their bonds and are calculated as single unbounded particles from there on.
Am I doing something wrong or is it just a limitation of the current model?

Thx and regards

Ermek Asylbekov | Fri, 06/23/2017 - 12:46

Hi,

I think I 've found the reason for it:
in bond_gran.cpp in lines 170ff. a bond is defined as broken if it overlaps the box borders.

I'll try to fix it for a periodic box. I'll post here if I happen to improve the code ;)

arnom's picture

arnom | Tue, 06/27/2017 - 15:33

Just a small hint, periodic boundaries have some issues in 3.6. In a few days we will release 3.7 and we hope to fix those issues with that. Obviously that does not solve any issues in the bond model by Richie.

DCS team member & LIGGGHTS(R) core developer