過去の記事で、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のみとなりますので、注意して確認する必要がありますね。