Hello,
I am working on a rotating blade moving into a granular medium and would like to compute the force/torque exerted by the particles on the mesh. I use Liggghts version 2.1.1 and the fix "mesh/surface/stress".
My problem is that the output is invariably zero (see the helix.txt output file).
Can somebody help me finding the origin of the problem ?
The interesting part of my code is:
fix cad1 all mesh/surface file blade.stl type 1 scale 10 move 0. 0. 0.3
(insert particles and run some timesteps)
fix blade_wall all wall/gran/hooke/history mesh n_meshes 1 meshes cad1
fix MyStress all mesh/surface/stress file blade.stl type 1 stress on reference_point 0. 0. 0.3
variable ForceX equal f_MyStress [1]
variable ForceY equal f_MyStress [2]
variable ForceZ equal f_MyStress [3]
variable TorqueX equal f_MyStress [4]
variable TorqueY equal f_MyStress [5]
variable TorqueZ equal f_MyStress [6]
fix DumpText all print 200 "${ForceX} ${ForceY} ${ForceZ} ${TorqueX} ${TorqueY} ${TorqueZ}" screen no title "Fx Fy Fz Tx Ty Tz " file helix.txt
fix movecad1 all move/mesh mesh cad1 rotate origin 0. 0. 0. axis 0. 0. 1. period 0.2
fix movecad2 all move/mesh mesh cad1 linear -0. 0. -0.3
run 35000
The blade rotates around a vertical axis and, at t=0, the lower point of the axle is located at (0.0,0.0,0.3). Simultaneously, the blade goes down, deeper into the granular medium.
Any help/hint will be warmly appreciated.
Best regards, Olivier
ckloss | Thu, 09/20/2012 - 16:15
Hi Olivier, can you step-wise
Hi Olivier,
can you step-wise "debug" the simulation to see at which level of complexity the zero-output evolves?
Cheers, Christoph
bonnefoy | Thu, 09/20/2012 - 17:33
Hi Andreas, Hi
Hi Andreas, Hi Christoph,
thank you for your comments, guys.
Andreas, I tried to write one single fix as you suggest but it appears that "stress" and "reference_point" are not keywords for the mesh/surface fix (they are keywords only for the mesh/surface/stress fix).
Christoph, I reduced the complexity to the minimum I can but the output is still zero.zero. The visualization with Paraview shows nothing particular: some grains are lying on the blade, the other ones are below.
Maybe my input script can help:
atom_style granular
boundary m m m
newton off
communicate single vel yes
units si
#region reg block -.1 0.6 -.1 0.2 -0.3 0.25 units box
#region reg block -1 1 -.5 .5 -.5 .5 units box
region reg cylinder z 0 0 0.25 -.5 .5 units box
create_box 1 reg
neighbor 0.05 bin
neigh_modify delay 0
fix m1 all property/global youngsModulus peratomtype 5.e6
fix m2 all property/global poissonsRatio peratomtype 0.45
fix m3 all property/global coefficientRestitution peratomtypepair 1 0.9
fix m4 all property/global coefficientFriction peratomtypepair 1 0.05
fix m5 all property/global characteristicVelocity scalar 2.
pair_style gran/hooke/history #Hooke without cohesion
pair_coeff * *
timestep 0.00005
fix 1 all nve/sphere
fix 2 all gravity 9.81 vector 0.0 0.0 -1.0
#box walls
fix cylindre all wall/gran/hooke/history primitive type 1 zcylinder 0.25 0. 0.
fix top all wall/gran/hooke/history primitive type 1 zplane +0.5
fix bottom all wall/gran/hooke/history primitive type 1 zplane -0.5
#import mesh from cad:
fix cad1 all mesh/surface file blade.stl type 1 scale 10 move 0. 0. 0.
#thermo settings
compute 1 all erotate/sphere
thermo_style custom step atoms ke c_1 vol
thermo 1000
thermo_modify lost ignore norm no
compute_modify thermo_temp dynamic yes
#particle insertion
region bc cylinder z 0 0 0.24 .16 .18 units box
group nve_group region reg
fix ins nve_group pour/legacy 1500 1 1 vol 0.7 1000 diam 0.08 0.08 dens 2500 2500 vel 0. 0. 0. 0. -0.8 region bc
fix blade_wall all wall/gran/hooke/history mesh n_meshes 1 meshes cad1
run 16000
unfix ins
dump dmp all custom 200 post/dump*.movingMesh id type type x y z ix iy iz vx vy vz fx fy fz omegax omegay omegaz radius
dump dumpstl all mesh/stl 200 post/dump*.stl
fix monstress all mesh/surface/stress file blade.stl type 1 stress on reference_point 0. 0. 0.
variable ForceX equal f_monstress[1]
variable ForceY equal f_monstress[2]
variable ForceZ equal f_monstress[3]
variable CoupleX equal f_monstress[4]
variable CoupleY equal f_monstress[5]
variable CoupleZ equal f_monstress[6]
fix DumpTexte all print 200 "${ForceX} ${ForceY} ${ForceZ} ${CoupleX} ${CoupleY} ${CoupleZ}" screen no title "Fx Fy Fz Tx Ty Tz " file post/force_helice.txt
run 15000
aaigner | Thu, 09/20/2012 - 17:44
That's why you should use fix
Andreas, I tried to write one single fix as you suggest but it appears that "stress" and "reference_point" are not keywords for the mesh/surface fix (they are keywords only for the mesh/surface/stress fix).
That's why you should use fix mesh/surface/stress instead of mesh/surface.
mesh/surface/stress has all keywords from mesh/stress + stress/reference_point/...
Cheers,
Andreas
bonnefoy | Thu, 09/20/2012 - 18:04
Andreas, You were right with
Andreas,
You were right with your suggestion. It works absolutely fine.
I misread part of it and did not notice that you also changed the fix name into mesh/surface/stress.
Thank you again and have a wonderful day (US Time)
Olivier
aaigner | Thu, 09/20/2012 - 16:49
Double insert, but only one wall?!
Hi Olivier,
maybe your problem is that your mesh 'MyStress' isn't a wall.
I think your code should look like:
fix cad1 all mesh/surface/stress file blade.stl type 1 scale 10 move 0. 0. 0.3 stress on reference_point 0. 0. 0.3
(insert particles and run some timesteps)
fix blade_wall all wall/gran/hooke/history mesh n_meshes 1 meshes cad1
variable ForceX equal f_cad1[1]
variable ForceY equal f_cad1[2]
variable ForceZ equal f_cad1[3]
variable TorqueX equal f_cad1[4]
variable TorqueY equal f_cad1[5]
variable TorqueZ equal f_cad1[6]
fix DumpText all print 200 "${ForceX} ${ForceY} ${ForceZ} ${TorqueX} ${TorqueY} ${TorqueZ}" screen no title "Fx Fy Fz Tx Ty Tz " file helix.txt
fix movecad1 all move/mesh mesh cad1 rotate origin 0. 0. 0. axis 0. 0. 1. period 0.2
fix movecad2 all move/mesh mesh cad1 linear -0. 0. -0.3
run 35000
Cheers,
Andreas