I have been running simple shear stress simulations in LIGGGHTS using moving layers of particles along the top and bottom of the domain to mimic a Lees-Edwards boundary condition. I am able to get reasonable agreement with kinetic theory for normal and shear stresses for spherical particles using the compute stress/atom command. However, when I create spherical superquadric particles to run the same cases, I would expect very close to the same, exact result. I consistently get normal and shear stresses about a factor of 2 below what the were for granular particles (which also matched kinetic theory).
Has anyone noticed issues when using compute stress/atom with superquadric particles? I will include some of my output command below so you can see how I am computing normal and shear stress.
fix pts1 all particletemplate/superquadric 123457 atom_type 1 density constant 2500 shape constant 0.015 0.015 0.015 blockiness constant 2.0 2.0
compute peratom nve_group stress/atom
compute p nve_group reduce/region main sum c_peratom[1] c_peratom[2] c_peratom[3] c_peratom[4]
variable pyy equal -(c_p[2])
variable pxy equal (c_p[4])
Thank you.
| Attachment | Size |
|---|---|
 in_granular.shearflow.txt | 2.28 KB |
| in_superquadric.shearflow.txt | 2.35 KB |
deepakpawar.2310 | Fri, 01/08/2021 - 19:21
Hi
Hi
I tried with the following snippt, and its working for the superquadric particles system.
#(1.1) Atom stresses (#stresses along xx, yy, zz, xy, xz and yz)---------------------------------------------------------
compute sig all stress/atom ke pair
compute pr all reduce sum c_sig[1] c_sig[2] c_sig[3] c_sig[4] c_sig[5] c_sig[6]
variable sigma11 equal -c_pr[1]
variable sigma22 equal -c_pr[2]
variable sigma33 equal -c_pr[3]
variable sigma12 equal c_pr[4]
variable sigma13 equal c_pr[5]
variable sigma23 equal c_pr[6]
fix data11 all print ${writeOutEvery} "${t}, $(v_sigma11), $(v_sigma22), $(v_sigma33), $(v_sigma12), $(v_sigma13), $(v_sigma23)" &
file stressTensor.csv screen no title "time, s11(KPa), s22, s33, s12, s13, s23"
Thanks
Happy learning !!
RonMexico | Fri, 01/08/2021 - 21:10
Thank you, but my issue is with the output itself
I believe there is a bug in the code (or I am not understanding something). I have attached two input scripts that seed five 0.15 m radius particles in a 1 m^3 domain. One uses the granular formulation and the other uses superquadric. You will see that all outputs are identical except pyy and pxy when there is a contact.
Here is the output from granular simulation:
Step Time Atoms KinEng trans rotat pyy pxy num_cont
1 1e-07 5 8060.6948 21.381799 -5.7326602 4795.3895 -2490.4302 0
1000 0.0001 5 8060.6948 21.381799 -5.7453222 4795.3895 -2490.4302 0
2000 0.0002 5 8060.6948 21.381799 -5.7578616 4795.3895 -2490.4302 0
3000 0.0003 5 8060.6948 21.381799 -5.7702636 4795.3895 -2490.4302 0
4000 0.0004 5 8060.6948 21.381799 -5.782526 4795.3895 -2490.4302 0
5000 0.0005 5 8060.6948 21.381799 -5.7946469 4795.3895 -2490.4302 0
6000 0.0006 5 8060.6948 21.381799 -5.806624 4795.3895 -2490.4302 0
7000 0.0007 5 8060.4635 21.380927 -5.8184554 5427.9698 -3117.822 1
8000 0.0008 5 8008.845 21.186193 -5.8300099 13223.901 -10823.086 1
9000 0.0009 5 7838.3753 20.543085 -5.840606 25215.75 -22615.898 1
10000 0.001 5 7552.9118 19.466156 -5.8491795 38608.575 -35709.411 1
11000 0.0011 5 7221.2099 18.214791 -5.8545704 50877.919 -47613.141 1
12000 0.0012 5 6950.8547 17.194859 -5.8557643 59737.879 -56098.184 1
13000 0.0013 5 6839.4475 16.774568 -5.8520987 63509.463 -59553.914 1
14000 0.0014 5 6925.7951 17.10032 -5.8434002 61476.82 -57326.944 1
15000 0.0015 5 7170.3327 18.022853 -5.8300205 54051.095 -49867.461 1
16000 0.0016 5 7476.9217 19.179479 -5.8127641 42676.011 -38622.719 1
17000 0.0017 5 7741.4355 20.177374 -5.7927265 29521.102 -25726.917 1
18000 0.0018 5 7897.8223 20.767353 -5.7710867 17130.823 -13653.745 1
19000 0.0019 5 7940.8137 20.92954 -5.7488873 8563.5623 -5347.7268 1
20000 0.002 5 7939.5877 20.924915 -5.7267706 8451.4866 -5239.6112 0
21000 0.0021 5 7939.5877 20.924915 -5.7048289 8451.4866 -5239.6112 0
22000 0.0022 5 7939.5877 20.924915 -5.6830609 8451.4866 -5239.6112 0
23000 0.0023 5 7939.5877 20.924915 -5.6614667 8451.4866 -5239.6112 0
24000 0.0024 5 7939.5877 20.924915 -5.6400465 815.76588 -516.16203 0
25000 0.0025 5 7939.5877 20.924915 -5.6188003 815.76588 -516.16203 0
For superquadric simulation:
Step Time Atoms KinEng trans rotat pyy pxy num_cont
1 1e-07 5 8060.6948 21.381799 -5.7326602 4795.3895 -2490.4302 0
1000 0.0001 5 8060.6948 21.381799 -5.7453222 4795.3895 -2490.4302 0
2000 0.0002 5 8060.6948 21.381799 -5.7578616 4795.3895 -2490.4302 0
3000 0.0003 5 8060.6948 21.381799 -5.7702636 4795.3895 -2490.4302 0
4000 0.0004 5 8060.6948 21.381799 -5.782526 4795.3895 -2490.4302 0
5000 0.0005 5 8060.6948 21.381799 -5.7946469 4795.3895 -2490.4302 0
6000 0.0006 5 8060.6948 21.381799 -5.806624 4795.3895 -2490.4302 0
7000 0.0007 5 8060.4635 21.380927 -5.8184554 4795.5146 -2490.5154 1
8000 0.0008 5 8008.845 21.186193 -5.8300099 4837.9386 -2523.6799 1
9000 0.0009 5 7838.3753 20.543085 -5.840606 4997.349 -2650.1928 1
10000 0.001 5 7552.9118 19.466156 -5.8491795 5301.4745 -2890.9035 1
11000 0.0011 5 7221.2099 18.214791 -5.8545704 5740.0814 -3235.0816 1
12000 0.0012 5 6950.8547 17.194859 -5.8557643 6270.8316 -3646.1632 1
13000 0.0013 5 6839.4475 16.774568 -5.8520987 6831.5838 -4073.0929 1
14000 0.0014 5 6925.7951 17.10032 -5.8434002 7357.8044 -4465.467 1
15000 0.0015 5 7170.3327 18.022853 -5.8300205 7798.5336 -4786.4002 1
16000 0.0016 5 7476.9217 19.179479 -5.8127641 8125.2057 -5018.3326 1
17000 0.0017 5 7741.4355 20.177374 -5.7927265 8332.1695 -5161.4723 1
18000 0.0018 5 7897.8223 20.767353 -5.7710867 8432.2871 -5228.5711 1
19000 0.0019 5 7940.8137 20.92954 -5.7488873 8453.4345 -5241.2767 1
20000 0.002 5 7939.5877 20.924915 -5.7267706 8451.4866 -5239.6112 0
21000 0.0021 5 7939.5877 20.924915 -5.7048289 8451.4866 -5239.6112 0
22000 0.0022 5 7939.5877 20.924915 -5.6830609 8451.4866 -5239.6112 0
23000 0.0023 5 7939.5877 20.924915 -5.6614667 8451.4866 -5239.6112 0
24000 0.0024 5 7939.5877 20.924915 -5.6400465 815.76588 -516.16203 0
25000 0.0025 5 7939.5877 20.924915 -5.6188003 815.76588 -516.16203 0
Any suggestions is appreciated. Thank you.
RonMexico | Thu, 01/14/2021 - 21:25
I have dug deeper into this, and the disagreement
appears to be in the pairwise contribution term output from compute pressure and compute stress/atom. I found this using the "ke" and "pair" keyword options within the compute command. The forces on the particles do match throughout the simulation between granular and superquadric models. To me, this points to a problem in the position vector r in the virial formula presented in the LIGGGHTS documentation. Has anyone noticed similar problems?
Thanks.
RonMexico | Tue, 01/19/2021 - 18:34
Going into this further, the position vectors also match
Using the atom vector variables, I have verified that both position and force values are identical between granular and superquadric simulations. This means the pairwise virial term should also be identical, but it is not. I suspect that the pairwise virial computation is not being called, or its result is getting altered somewhere in the superquadric code. The other possibility is that I am missing a flag that is required for superquadric but not granular formulation in the input files. Both input files are attached in the first post.