過去の記事で、Liggghtsで非球形粒子を扱う方法を２つ紹介しました。
- 超二次関数を用いた方法
- 粒子を重ね合わせて再現する方法(clump)
今回はこの２つの方法を比較してみようと思います。
比較
同じ形状を２つの方法でモデリングし、その結果を比較することで検証してみようと思います。
対象形状は以下とします。
超二次関数のパラメーターは以下の通りです。
超二次関数 \left(\left|\frac{x}{a}\right|^{n_2} + \left|\frac{y}{b}\right|^{n_2}\right)^\frac{n_1}{n_2} + \left|\frac{z}{c}\right|^{n_1} = 1
\text{パラメータ}
\begin{cases}
\text{$a$} = 0.05\\\
\text{$b$} = 0.05\\\
\text{c} = 0.05\\\
\text{$n_1$} = 10\\\
\text{$n_2$} = 10
\end{cases}
また、以下が複数の粒子を重ね合わせてモデル化するスクリプトです。任意のstlに書き換えて、粒径のパラメーターを変えるなどして使用できます。
clump.
# Multisphere clump generation from surface mesh
create_multisphere_clump dmin absolute 2.5 rmin 10 pmax 0.1 seed 10000 &
surfacefile meshes/super_ellipsoid.stl clumpfile data/clump_sup.dat post/clump_sup.vtk
検証解析
複数のDEM粒子を発生させ、自由落下時の挙動を解析します。
それぞれのインプットファイルは以下の通りです。
in.multishere
atom_style sphere
atom_modify map array sort 0 0
boundary f f f
newton off
communicate single vel yes
processors * * 1
units si
region reg block -0.8 0.8 -0.8 0.8 0.0 2.0 units box
create_box 1 reg
neighbor 0.004 bin
neigh_modify delay 0
# material properties required for granular pair styles
fix m1 all property/global youngsModulus peratomtype 1.e7
fix m2 all property/global poissonsRatio peratomtype 0.45
fix m3 all property/global coefficientRestitution peratomtypepair 1 0.3
fix m4 all property/global coefficientFriction peratomtypepair 1 0.5
# pair style
pair_style gran model hertz tangential history
pair_coeff * *
timestep 0.00001
fix gravi all gravity 9.81 vector 0.0 0.0 -1.0
fix zwalls all wall/gran model hertz tangential history primitive type 1 zplane 0.0
# templates and distribution for insertion
fix pts1 all particletemplate/multisphere 1272435 atom_type 1 density constant 2500 &
nspheres 159 ntry 1000000 spheres file data/clump_sup.dat scale 0.001 type 1
fix pdd1 all particledistribution/discrete 100 1 pts1 1.0
# region and insertion
region bc cylinder z 0.0 0.0 0.15 1.2 2.0 units box
fix bc all mesh/surface file msh/box.stl type 1
fix inface all mesh/surface file msh/plane.stl type 1
#fix ins all insert/stream seed 1001 d istributiontemplate pdd1 nparticles 1 vel constant 0. 0. -1.0 &
# orientation constant 10 10 10 10 omega constant 5.0 5.0 0.0 particlerate 100000 overlapcheck yes insertion_face inface extrude_length 0.4
fix ins all insert/pack seed 100000 distributiontemplate pdd1 vel constant 0. 0. -1.0 &
orientation random omega constant 5.0 5.0 0.0 &
insert_every once overlapcheck yes region bc ntry_mc 100000 particles_in_region 20
# integrator for multisphere rigid bodies
fix integr all multisphere
fix ts all check/timestep/gran 1000 0.1 0.1
# output settings
compute rke all erotate/multisphere
thermo_style custom step atoms ke c_rke
thermo 1000
thermo_modify lost ignore norm no
compute_modify thermo_temp dynamic yes
# insert the particles so that dump is not empty
run 1
dump dmp all custom/vtk 500 post/mlclump*.vtk id type mass x y z vx vy vz fx fy fz omegax omegay omegaz radius
run 200000 upto
in.superquadric
atom_style superquadric
atom_modify map array sort 0 0
boundary f f f
newton off
communicate single vel yes
processors * * 1
units si
region reg block -0.8 0.8 -0.8 0.8 0.0 2.0 units box
create_box 1 reg
neighbor 0.004 bin
neigh_modify delay 0
# material properties required for granular pair styles
fix m1 all property/global youngsModulus peratomtype 1.e7
fix m2 all property/global poissonsRatio peratomtype 0.45
fix m3 all property/global coefficientRestitution peratomtypepair 1 0.3
fix m4 all property/global coefficientFriction peratomtypepair 1 0.5
# pair style
pair_style gran model hertz tangential history surface superquadric
pair_coeff * *
timestep 0.00001
fix gravi all gravity 9.81 vector 0.0 0.0 -1.0
fix zwalls all wall/gran model hertz tangential history surface superquadric primitive type 1 zplane 0.0
# templates and distribution for insertion
fix pts1 all particletemplate/superquadric 10 atom_type 1 density constant 2500 shape constant 0.05 0.05 0.05 blockiness constant 10.0 10.0
fix pdd1 all particledistribution/discrete 100 1 pts1 1.0
# region and insertion
region bc cylinder z 0.0 0.0 0.15 1.2 2.0 units box
fix bc all mesh/surface file msh/box.stl type 1
fix inface all mesh/surface file msh/plane.stl type 1
#fix ins all insert/stream seed 1001 d istributiontemplate pdd1 nparticles 1 vel constant 0. 0. -1.0 &
# orientation constant 10 10 10 10 omega constant 5.0 5.0 0.0 particlerate 100000 overlapcheck yes insertion_face inface extrude_length 0.4
fix ins all insert/pack seed 10 distributiontemplate pdd1 vel constant 0. 0. -1.0 &
orientation random omega constant 5.0 5.0 0.0 &
insert_every once overlapcheck yes region bc ntry_mc 100000 particles_in_region 20
# integrator for multisphere rigid bodies
fix integr all nve/superquadric integration_scheme 1
fix ts all check/timestep/gran 1000 0.1 0.1
# output settings
compute rke all erotate/superquadric
thermo_style custom step atoms ke c_rke f_ts[1] f_ts[2] vol
thermo 1000
thermo_modify lost ignore norm no
compute_modify thermo_temp dynamic yes
# insert the particles so that dump is not empty
run 1
dump dmp all custom/vtk 500 post/supclump*.vtk id type mass x y z vx vy vz fx fy fz omegax omegay omegaz radius &
omegay omegaz radius shapex shapey shapez quat1 quat2 quat3 quat4 blockiness1 blockiness2 tqx tqy tqz angmomx angmomy angmomz
run 200000 upto
multispereの場合は、あらかじめ作成したclumpファイルを使用して粒子を生成しています。
2種類の粒子タイプでそれぞれ解析を実行可能です。
結果比較
２つの解析動画を示します。
全く同じ粒子配置で落下させているわけではないので、完全な比較は難しいですが、おおよそ同じような挙動になっているように思います。multisphereの方が粒子が回転に対して抵抗が働いているようにも見受けられますが、複数の粒子を重ね合わせて形状を近似している影響かもしれません。
今回は、正方形ということで超二次関数を用いた方法でも粒子を生成できましたが、もっと複雑な形状はmultisphereのみとなりますので、注意して確認する必要がありますね。