2D simulation

>>Could this affect something?
yes - but it will make it more realistic

I would use a 3D (slot) model with periodic boundary conditions - simply because granular physics is not 2D - it's 3D. That's also the reason why I personally won't implement 2D for new features. Imho 2D DEM is not state of the art any longer

Christoph

Alright, that was a good hint.

Right now the simulation works fine, but paraview crashes at certain dump-files, and I don't have a clue what the problem is...

All my particles are falling after been inserted by insert/pack, touching the ground, touching each other, time is passing and suddenly, the whole simulation box converts into a big bubble (at dump-file 136)...

Can anyone please, check my input-script out and tell me, what the problem is?

Thanks very much in advance,
Sebastian

####################################################
# LIGGGHTS version: 1-5-2
# ##################################################

variable NAME index Test2D_mesh_nein_10

atom_style granular
atom_modify map array
boundary f f f
newton off

communicate single vel yes

units si

variable scale equal 1
variable scale equal 0.4 # just for first glance
variable scaleCAD equal 0.001*${scale}

variable Sp_Hoehe_ equal 2
variable Sp_Radius_ equal 0.7
variable Dauer equal 21600 #in [s] = 6h
variable Partikelradius equal 0.0315
variable Partikeldichte equal 3000

variable Sp_Hoehe_skal equal ${Sp_Hoehe_}*${scale}
variable Sp_Radius_skal equal ${Sp_Radius_}*${scale}

variable Einfuellhoehe equal 2*${Sp_Hoehe_skal}-15*${Partikelradius}
variable Einfuellmitte_plus equal ${Sp_Hoehe_skal}+${Partikelradius}
variable Einfuellmitte_minus equal ${Sp_Hoehe_skal}-${Partikelradius}
variable Einfuellradius equal ${Sp_Radius_skal}-2*${Partikelradius}

variable Einheitsdicke equal 0.032

region reg block -${Sp_Radius_} ${Sp_Radius_} -${Einheitsdicke} ${Einheitsdicke} 0 ${Einfuellhoehe} units box
region reg_1 block -${Einfuellradius} ${Einfuellradius} -${Einheitsdicke} ${Einheitsdicke} 0 ${Einfuellmitte_minus} units box
region reg_2 block -${Einfuellradius} ${Einfuellradius} -${Einheitsdicke} ${Einheitsdicke} ${Einfuellmitte_plus} ${Einfuellhoehe} units box

create_box 1 reg

neighbor 0.006 bin # 0.1*min(particle diameter)
neigh_modify delay 0

variable YM equal 30.e9 #YoungsModulus
timestep 0.000001

# für schnelles Modell # just for first glance
variable YM equal 30.e7 #YoungsModulus
timestep 0.00001

#Material properties required for new pair styles

fix m1 all property/global youngsModulus peratomtype ${YM}
fix m2 all property/global poissonsRatio peratomtype 0.250
fix m3 all property/global coefficientRestitution peratomtypepair 1 0.1 #je kleiner der Wert desto schneller ist die Schüttung gesetzt. Volker hat Parameterstudie dazu gemacht, auf Thermomechanik hat dieser Wert keinen Einfluss
fix m4 all property/global coefficientFriction peratomtypepair 1 0.5
fix m5 all property/global characteristicVelocity scalar 2.

#New pair style
pair_style gran/hooke/history 1 0 #Hertzian without cohesion
pair_coeff * *

fix gravi all gravity 9.81 vector 0.0 0.0 -1.0

#walls
fix force1 all store/force
fix xwalls all wall/gran/hooke/history 1 0 xplane -${Sp_Radius_skal} ${Sp_Radius_skal} 1
fix zwalls all wall/gran/hooke/history 1 0 zplane 0 NULL 1
fix force2 all store/force

#particle distributions
fix pts1 all particletemplate/sphere 1 atom_type 1 density constant ${Partikeldichte} radius constant ${Partikelradius}
fix pdd1 all particledistribution/discrete 1. 1 pts1 1

#region and insertion
fix ins all insert/pack seed 100001 distributiontemplate pdd1 &
insert_every once overlapcheck yes all_in yes vel constant 0. 0. -1. &
volumefraction_region 0.63 region reg_1 ntry_mc 1000

#apply nve integration
fix integr all nve/sphere

variable dump equal 100
variable Nfreq equal v_dump
variable Nevery equal v_Nfreq/10
variable Nrepeat equal v_Nfreq/v_Nevery

#output settings, include total thermal energy
compute 1 all erotate/sphere
thermo_style custom step atoms ke c_1 vol
thermo ${dump}
thermo_modify lost ignore norm no
compute_modify thermo_temp dynamic yes

# letztes fix
fix setforce all setforce NULL 0 NULL

#insert the first particles so that dump is not empty
run 1

variable xwallforce atom f_force2[1]-f_force1[1]
variable ywallforce atom f_force2[2]-f_force1[2]
variable zwallforce atom f_force2[3]-f_force1[3]

variable betrag_kwx atom sqrt(xwallforce*xwallforce)
variable betrag_kwy atom sqrt(ywallforce*ywallforce)
variable betrag_kwz atom sqrt(zwallforce*zwallforce)

# => Diese Mittelung hier in LIGGGHTS bringts nicht, da ansonsten die Partikel, die keine Wandkraft erfahren in die Mittelung mit einbezogen werden
# => Auswerten der variablen betrag_kw in Matlab
#fix ave_F_Wand all ave/spatial ${Nevery} ${Nrepeat} ${Nfreq} z lower ${Hoehenmittelung} betrag_kwx betrag_kwy betrag_kwz ave one file ave_F_Wand_betrag units box
#fix ID group-ID ave/spatial Nevery Nrepeat Nfreq dim origin delta value1 value2 ... keyword args ...

# Mittelung der Schüttbettbewegung über die Höhe
variable Hoehenmittelung equal ${Sp_Hoehe_skal}/10
fix ave_v all ave/spatial ${Nevery} ${Nrepeat} ${Nfreq} z lower ${Hoehenmittelung} vx vy vy ave one file post_${NAME}/00_ave_v units box

dump dmp_setzen all custom ${dump} post_${NAME}/dump*.setzen id type type x y z ix iy iz vx vy vz fx fy fz omegax omegay omegaz radius v_xwallforce v_ywallforce v_zwallforce f_force1[1] f_force1[2] f_force1[3]
dump_modify dmp_setzen first yes

restart ${dump} post_${NAME}/restart.${NAME}.

#insert particles
run 1500000
undump dmp_setzen