Hello everybody,
I am currently writing a rigid bond contact model (as from "A bonded-particle model for rock" Potyondy and Cundall 2004) using the contact style cohesion. I have everything working, but I came across an odd behavior:
In the attached case I have a beam with one side fixed and a force acting on the other side. Depending on the direction of the force ( fix F_force pull addforce 0 0 -0.05 or fix F_force pull addforce 0 0 0.05) the trajectory is different! I have found that the point, where the behaviour differs is, when the reordering of the particle indices starts.
Case: fix F_force pull addforce 0 0 0.05
In the dump file, there is no reordering of the indices, the list remains steadily increasing from 1-16.
Case: fix F_force pull addforce 0 0 -0.05
From timestep 37000 to 38000 the list is reordered and that is when the trajectory starts to change
I could not find anything in other contact models, that refers to any kind of reordering and I could not find me using the atom->tag instead of the local index. But I may be wrong. If anyone has an idea what could cause this, I would be greatful!
Regards,
Valentin
P.S.: I am using LIGGGHTS 3.5.0
| Attachment | Size |
|---|---|
 cohesion_model_rigid_bond_2.h.tar_.gz | 40 KB |
| data.beam_.tar_.gz | 10 KB |
| in.liggghts.tar_.gz | 10 KB |
richti83 | Tue, 03/14/2017 - 15:39
unfortunately the file links
unfortunately the file links are broken, could you upload the header file once again ? (maybe compressed as tar.gz)
Baric | Wed, 03/15/2017 - 09:26
New files uploaded
The files in the first post are updated now.
When I was comparing my code to cohesion_model_easo_capillary_viscous.h I have found that sidata has sidata.cri/sidata.crj, which scdata does not provide. What is this value? In surfacesClose radi/radj is used instead.
richti83 | Wed, 03/15/2017 - 12:11
I have had a short look in
I have had a short look in the code and It compiled nicely in 3.6.0, thanks for shoring this ! I could not test your issue because the in.liggghts file contains only zeros after unpacking.
Baric | Wed, 03/15/2017 - 12:59
Weird
I have no idea what went wrong. I have uploaded it again and checked it, it looks good for now. Just in case, here is the content:
################################################################################
# Define the atomstyle to be a hybrid of granular and molecular to use the molecule flag
atom_style hybrid granular molecular
units si
boundary f f f
newton off
variable interv equal 1000
communicate single vel yes
neighbor 3.0e-03 bin
neigh_modify delay 0
# Read the inputdata from datafile
read_data data.beam
# Read the include file with the potentials and the interaction parameters
# Kn, Kt, gamman, gammat, etc.
fix m1 all property/global coefficientRollingFriction peratomtypepair 1 0.5
fix m2 all property/global coefficientFriction peratomtypepair 1 0.5
fix m3 all property/global youngsModulus peratomtype 5.68e9
fix m4 all property/global poissonsRatio peratomtype 0.5
fix m5 all property/global coefficientRestitution peratomtypepair 1 0.1
fix m8 all property/global bondTensileStrength peratomtypepair 1 456e6
fix m9 all property/global bondTensileStrengthShear peratomtypepair 1 456e6
fix m10 all property/global bondRadius peratomtypepair 1 1.25e-3
fix m11 all property/global bondMaxSeparationDistance peratomtypepair 1 0.0055
fix m12 all property/global bondSn peratomtypepair 1 3e9
fix m13 all property/global bondSt peratomtypepair 1 1e9
fix m14 all property/global bondRf peratomtypepair 1 0.5
pair_style gran model hertz tangential history cohesion rigid/bond2 rolling_friction epsd2 full_length on per_molecule on
pair_coeff * *
shell rm generated_bonds.txt
timestep 1e-7
run 2
group freeze id 1
group pull id 16
group rest subtract all freeze
fix F_nve rest nve/sphere
fix F_frozen1 freeze setforce 0.0 0.0 0.0
fix F_force pull addforce 0 0 0.05
shell mkdir Post
shell rm broken_bonds.txt
thermo ${interv}
dump trjdmp all custom ${interv} Post/smd.trj id mol type radius x y z vx vy vz fx fy fz omegax omegay omegaz
run 500000
Baric | Thu, 03/23/2017 - 14:28
Problem solved
I finally got it...
I was using the hsetup as used in other models. In most of them only a contactflag is stored but I was storing vectors there. Anyways, the bond was set for example particle 1 and 2. Therefore, all relative velocities and the deltax, deltay and so on, were calculated for particle 1 - particle 2. If all of the sudden the combination is switched to particle 2 - particle 1, these calculations are not correct anymore.
Therefore, I had to include a function that all those differences are always for the same combination (e.g. 1-2) and not the flipped one and furthermore, the forces hat to be transferred in the correct way at the end.
Thanks to everyone who was trying to help me. As soon as I have cleaned my code and made the final verification steps I will share the model with you.
Regards,
Valentin