Hello people. I made another thread in developer forum, asking how to get the force in a simulation using granular bonds. The problem is that the values shown in dump files, aswell as the values calculated with compute property/local are way low than they should be (for the example i uploaded it should be around 31N, and its showing in both cases something close to 1e-6). What i need to validate y simulation, a tensile test, i need to evaluate the forces at bonds, as they are calculated in bond_gran.cpp, specifically the values defined from line 358 to 361: nforce_mag, tforce_mag, ntorque_mag and ttorque_mag.'?
Is there a way i can do this using commands inside liggghts? or should i code it inside the program?
Im uploading the example code here too, a single bond between two particles.
Thanks in advance
| Attachment | Size |
|---|---|
 simple.txt | 2.76 KB |
mschramm | Wed, 06/07/2017 - 15:45
mag values
lines 358 to 361 give you the bondhistlist[n][i] values for the nth bond
you can get these values using
compute b1 all property/local batom1x batom1y batom1z batom2x batom2y batom2z batom1 batom2 btype bforceX bforceY bforceZ
For example bforceX will contain bondhistlist[n][0] + bondhistlist[n][3]
I do not believe there is a way to get torque values.
You could extend the property local command to include btorqueX through btorqueZ
mschramm | Wed, 06/07/2017 - 23:24
Reduce your time step
I downloaded you code and using r == 1 (in the code it is lambda and not your bond radius), I got a bond force value of 31.419047 when it broke.
I am using a time step of
variable tstep equal 0.0000001*${raytime}
and the bond broke at step 9936
Also, since you only have two spheres and one is frozen while the other is being moved, you do not need the command
fix integr nve_group nve/sphere
Setting up run at Wed Jun 7 16:22:08 2017
Memory usage per processor = 10.9194 Mbytes
Step Atoms numbond KinEng bfx bfy bfz
0 2 0 0 -1e+20 -1e+20 -1e+20
creating bond btw atoms 0 and 1 (i has now 1 bonds) at step 1
creating bond btw atoms 1 and 0 (i has now 1 bonds) at step 1
Created 1 bonds at timestep 1
1000 2 1 6.5449847e-11 0 0 3.1621423
2000 2 1 6.5449847e-11 0 0 6.3242847
3000 2 1 6.5449847e-11 0 0 9.4864271
4000 2 1 6.5449847e-11 0 0 12.648569
5000 2 1 6.5449847e-11 0 0 15.810712
6000 2 1 6.5449847e-11 0 0 18.972854
7000 2 1 6.5449847e-11 0 0 22.134997
8000 2 1 6.5449847e-11 0 0 25.297139
9000 2 1 6.5449847e-11 0 0 28.459281
broken bond 0 at step 9936
it was nstress
sigman_break == 4.000000e+07
mag_force == 4.000397e+07
tau_break == 2.309401e+07
mag_force == 2.793713e-13
10000 2 1 6.5449847e-11 0 0 31.419047
11000 2 1 6.5449847e-11 0 0 31.419047
12000 2 1 6.5449847e-11 0 0 31.419047
13000 2 1 6.5449847e-11 0 0 31.419047
14000 2 1 6.5449847e-11 0 0 31.419047
15000 2 1 6.5449847e-11 0 0 31.419047
16000 2 1 6.5449847e-11 0 0 31.419047
felipe.merino | Thu, 06/08/2017 - 17:37
WOW
Many thanks for your answer, as you know from the other post, i tried this all, knowing that the first coefficient represent a constant and is not the radius itself. It resulted to be the fix nve/sphere. For some reason, that command is messing up the coefficient. Only using low values for lambda (as the radius) it broke at the right timestep, but with a wrong force value. Using the fix with lambda = 1.0 doesnot break the bond.
I really do not fully understand that fix, or why it could be messing the force values.
Again thanks for your response and patience
mschramm | Thu, 06/08/2017 - 18:35
double counting
If you leave the fix in, you should see a warning message stating that "some particles are integrated twice". nve/sphere updates position, velocity, and angular values, but freeze and move alter those updates. This might have been the reason.