Drop CAD under gravity

Submitted by bchang5 on Mon, 02/24/2020 - 20:44

Hi all,

Is it possible to drop a CAD file under the action of gravity? For example, I would like to drop a sphere CAD file into a bed of particles. I would like to see the depth at which the sphere stops moving.

I looked at the "movecad" command, but it seems like it is actively driving the CAD file with set kinematic parameters. But I want my CAD to start from rest, fall into bed of particles, and stop naturally by the action of the particles resisting. Is this possible?

Best,
Brian

Forums:

mschramm | Wed, 02/26/2020 - 16:45

Hello,
You are going to want the rigid body solver
https://www.cfdem.com/rigid-body-dynamics-6-degrees-freedom

The problem... This is in the private release.
You can do this in a "dirty" way in your input file inside a loop or you can write your own fix class to do this.
I have done this once in the input file to get a pendulum to collide with my particles.

You are basically setting up an integration scheme in your input file where you set your STL file with initial conditions and every so many steps, update your stl file's velocity variables.

bchang5 | Tue, 03/03/2020 - 19:19

Do you have the input file on hand?

This makes sense for the free-fall case, but once the stl makes contact with the particles, how would I update the stl velocity?
Would I update the stl velocity variable based on the velocity of particles that are in contact with the stl?

tjleps | Tue, 03/03/2020 - 20:38

I think the idea is you sum the forces on the surface and then integrate F=ma. Basically make your own mini integrator inside of the LIGGGHTS framework.

mschramm | Tue, 03/03/2020 - 22:10

You are able to extract the forces and moments acting on the center of gravity of the file. So you can use this to generate a simple update scheme.
Since my system is constrained, I simply looked at the forces acting on the plunger to determine how to update it.
I am also updating the stl file at a different time step to speed up the computation a bit.
Ideally you would want to do this inside of a fix but I don't do much with MBD so I have not taken the time to do it.
(Input file is for use with my liggghts_flexible_fiber code)
Missing: All stl files

variable bondYoung equal 1.0e8
variable conYoung equal 1.0e8
variable poiRatio equal 0.3

variable bondShear equal ${bondYoung}/(2.0*(1+${poiRatio}))

atom_style hybrid granular bond/gran n_bondtypes 1 bonds_per_atom 6
hard_particles yes
atom_modify map array
boundary f f f
newton off
communicate single vel yes
shell mkdir post
units si
pair_style gran model hertz tangential history
bond_style gran
region sim_reg block -0.1 0.1 -0.25 0.10 0.0 0.25
create_box 1 sim_reg

neighbor 0.001 bin
neigh_modify delay 0

pair_coeff * *
bond_coeff 1 1.0 0.0 ${bondYoung} ${bondShear} 1 0.1 1 1.0e32 1.0e32

fix m1 all property/global youngsModulus peratomtype ${conYoung}
fix m2 all property/global poissonsRatio peratomtype ${poiRatio}
fix m3 all property/global coefficientRestitution peratomtypepair 1 0.3
fix m4 all property/global coefficientFriction peratomtypepair 1 0.3

fix pts1 all particletemplate/multiplespheres 15485863 atom_type 1 density constant 425 nspheres 61 ntry 150000 spheres &
0.000 0.0 0.0 0.0015 &
0.003 0.0 0.0 0.0015 &
0.006 0.0 0.0 0.0015 &
0.009 0.0 0.0 0.0015 &
0.012 0.0 0.0 0.0015 &
0.015 0.0 0.0 0.0015 &
0.018 0.0 0.0 0.0015 &
0.021 0.0 0.0 0.0015 &
0.024 0.0 0.0 0.0015 &
0.027 0.0 0.0 0.0015 &
0.030 0.0 0.0 0.0015 &
0.033 0.0 0.0 0.0015 &
0.036 0.0 0.0 0.0015 &
0.039 0.0 0.0 0.0015 &
0.042 0.0 0.0 0.0015 &
0.045 0.0 0.0 0.0015 &
0.048 0.0 0.0 0.0015 &
0.051 0.0 0.0 0.0015 &
0.054 0.0 0.0 0.0015 &
0.057 0.0 0.0 0.0015 &
0.060 0.0 0.0 0.0015 &
0.063 0.0 0.0 0.0015 &
0.066 0.0 0.0 0.0015 &
0.069 0.0 0.0 0.0015 &
0.072 0.0 0.0 0.0015 &
0.075 0.0 0.0 0.0015 &
0.078 0.0 0.0 0.0015 &
0.081 0.0 0.0 0.0015 &
0.084 0.0 0.0 0.0015 &
0.087 0.0 0.0 0.0015 &
0.090 0.0 0.0 0.0015 &
0.093 0.0 0.0 0.0015 &
0.096 0.0 0.0 0.0015 &
0.099 0.0 0.0 0.0015 &
0.102 0.0 0.0 0.0015 &
0.105 0.0 0.0 0.0015 &
0.108 0.0 0.0 0.0015 &
0.111 0.0 0.0 0.0015 &
0.114 0.0 0.0 0.0015 &
0.117 0.0 0.0 0.0015 &
0.120 0.0 0.0 0.0015 &
0.123 0.0 0.0 0.0015 &
0.126 0.0 0.0 0.0015 &
0.129 0.0 0.0 0.0015 &
0.132 0.0 0.0 0.0015 &
0.135 0.0 0.0 0.0015 &
0.138 0.0 0.0 0.0015 &
0.141 0.0 0.0 0.0015 &
0.144 0.0 0.0 0.0015 &
0.147 0.0 0.0 0.0015 &
0.150 0.0 0.0 0.0015 &
0.153 0.0 0.0 0.0015 &
0.156 0.0 0.0 0.0015 &
0.159 0.0 0.0 0.0015 &
0.162 0.0 0.0 0.0015 &
0.165 0.0 0.0 0.0015 &
0.168 0.0 0.0 0.0015 &
0.171 0.0 0.0 0.0015 &
0.174 0.0 0.0 0.0015 &
0.177 0.0 0.0 0.0015 &
0.180 0.0 0.0 0.0015 &
bonded yes

fix pdd1 all particledistribution/discrete 32452843 1 pts1 1.0

region ins_reg block -0.1 0.1 -0.1125 -0.06 0.05 0.25
fix ins all insert/rate/region seed 32452867 distributiontemplate pdd1 &
maxattempt 100000 nparticles 10 particlerate 5000000 insert_every 1000 &
orientation template overlapcheck yes region ins_reg ntry_mc 100000

#fix pen_arm1 all mesh/surface/stress file STL/Pendulum.stl type 1 scale 0.001 &
# move 0.0 0.0 0.5
fix pen_arm2 all mesh/surface/stress file STL/PendulumSmall.stl type 1 scale 0.001 &
move 0.0 0.0 0.5
fix u_hold1 all mesh/surface/stress file STL/U_Hold.stl type 1 scale 0.001 &
move 0.06 0.003 0.088
fix u_hold2 all mesh/surface/stress file STL/U_Hold.stl type 1 scale 0.001 &
move -0.06 0.003 0.088

fix mesh_wall all wall/gran model hertz tangential history mesh n_meshes 3 &
meshes pen_arm2 u_hold1 u_hold2

variable dt equal 1.0e-7
timestep ${dt}
fix integr all nve/sphere
variable deg equal 0.0
variable degPrint equal 180.0*${deg}/PI
variable simTime equal ${dt}*step
variable d_deg equal -5
thermo_style custom step atoms numbond v_simTime ke v_degPrint v_d_deg cpuremain cpu
thermo 5000

run 1

dump dmp all custom 5000 post/dump_*.liggghts id type x y z ix iy iz vx vy vz fx fy fz omegax omegay omegaz radius
#dump stl_pen1 all mesh/vtk 5000 post/stl_pen1_*.vtk stress stresscomponents pen_arm1
dump stl_pen2 all mesh/vtk 5000 post/stl_pen2_*.vtk stress stresscomponents pen_arm2
dump stl_hol all mesh/vtk 5000 post/stl_hol_*.vtk stress stresscomponents u_hold1 u_hold2
variable bondSkin equal 1.00000001*0.003

run 5000 upto
unfix ins
fix bondcr all bond/create/gran 1 1 1 ${bondSkin} 1 6
run 1
unfix bondcr
fix grav all gravity 9.81 vector 0.0 0.0 -1.0
velocity all set 0.0 0.0 -0.5
run 250000 upto

#fix move1 all move/mesh mesh pen_arm1 rotate/variable &
# origin 0. 0. 0.5 axis 1. 0. 0. omega v_d_deg
fix move2 all move/mesh mesh pen_arm2 rotate/variable &
origin 0. 0. 0.5 axis 1. 0. 0. omega v_d_deg

variable curStep equal step
variable simTime equal ${dt}*(step-${curStep})
variable totTime equal 0.5
variable totStep equal ceil(${totTime}/${dt})+${curStep}
variable stepCount equal 0
variable pengL equal 1.0
variable pengM equal 1.0
variable r2p equal 180.0/PI
variable penMu equal 5000
variable dtPen equal ${penMu}*${dt}
variable stepV equal ceil(${totStep}/${penMu})
variable dd_deg equal 0.0
label start_loop
run ${penMu} pre yes post no
variable stepCount equal ${stepCount}+1
variable d_deg equal ${d_deg}+0.5*${dtPen}*${dd_deg}
variable deg equal ${deg}+${dtPen}*${d_deg}
variable f1 equal -9.81*sin(${deg})/${pengL}
variable f2 equal f_pen_arm2[2]*cos(${deg})/(${pengL}*${pengM})
variable f3 equal f_pen_arm2[3]*sin(${deg})/(${pengL}*${pengM})
variable f4 equal -0.1*${d_deg}
variable dd_deg equal ${f1}+${f2}+${f3}+${f4}
variable d_deg equal ${d_deg}+0.5*${dtPen}*${dd_deg}
variable degPrint equal ${r2p}*${deg}
if "${stepCount}>${stepV}" then "jump in.liggghts end_loop"
# if "${deg}>0" then "jump in.liggghts end_loop"
jump in.liggghts start_loop
label end_loop

tjleps | Tue, 03/03/2020 - 11:53

If the ratio of the size of the sphere to the size of your grains isn't too large, you may be able to get away with using a large grain as your impactor, but your computational efficiency will take a bad hit as a result. If you don't care about lateral displacement you may be able to get away with a servo mesh with a maximum velocity equal to the analytically determined impact velocity, and a proper choice of ki and target_val for your desired mass.

bchang5 | Tue, 03/03/2020 - 19:13

Well, I am not actually using a sphere (it was just an example). I have some CAD models of various shapes (cubes, pyramids, cones) that I'd like to simulate dropping into granular media.

But taking your idea a little further, one possible work around that I will try is using the multisphere package. This would allow me to use individual spherical particles to define the shape of my CAD file.
My thought is to:
1. Generate a bed of spherical particles
2. Generate one multisphere particle (shaped as my CAD file) from a particular height and let it fall into bed of spherical particles.

Not sure if it will work though, as I'd like to obtain the forces acting on the CAD/multisphere particle. Do we know how forces and torques are calculated on a multisphere particle?