はじめに
DEMでは個々の粒子を球形として扱うことが多いです。これは粒子同士の接触判定・計算が単純であるためです。
では、非球形のものを解析したい場合はどのようにモデル化するのでしょうか？
よく用いられる方法として以下2つがあげられます。
-
球形粒子を重ね合わせることで再現
完全な表現は不可能ですが、多くの細かな粒子を重ね合わせれば、おおよそ再現が可能であり、複雑な形状でもモデル化できます。また、球形粒子をベースとしたモデル化になるため、接触判定・計算は単一粒子ものを拡張して実装可能です。 -
超二次関数を用いて形状を表現
超二次関するという方程式で目的の形状を再現します。ただし、この方程式で表現できるのは立方体や円柱など比較的単純なものに限られます。また、接触判定・計算のコストが高く大量の粒子を使用したシミュレーションには向いていません。
超二次関数 \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
超二次関数を用いる方法が実装されている商用ソフトはあまりないですが、Liggghtsでは実装されています。今回は、Liggghtsでの超二次関数を用いたシミュレーションを紹介したいと思います。
対象形状
超二次関数では、$a$・$b$・$c$・$n_1$・$n_2$が形状パラメータになっており、今回以下の値とすることで立方体を再現します。
\text{パラメータ}
\begin{cases}
\text{$a$} = 0.0085\\\
\text{$b$} = 0.0075\\\
\text{c} = 0.006\\\
\text{$n_1$} = 10\\\
\text{$n_2$} = 10
\end{cases}
この立方体形状でシミュレーションを行います。
解析
CFDEMのチュートリアル「conveyor」の内容を書き換えて非球形粒子のシミュレーションを実行してみます。
オリジナルでは、球形粒子となっていますが立方体の粒子になるよう書き換えます。
#Conveyor
atom_style superquadric
atom_modify map array
boundary m m m
soft_particles yes
newton off
communicate single vel yes
units si
region reg block -0.5 0.5 -0.2 0.2 -0.2 0.35 units box
create_box 1 reg
neighbor 0.01 bin
neigh_modify delay 0
#Material properties required for new pair styles
fix m1 all property/global youngsModulus peratomtype 1.e6
fix m2 all property/global poissonsRatio peratomtype 0.3
fix m3 all property/global coefficientRestitution peratomtypepair 1 0.2
fix m4 all property/global coefficientFriction peratomtypepair 1 0.6
#New 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
#granular walls
fix cv all mesh/surface file meshes/conveyor.stl type 1 surface_vel -4.5 0. 0.
fix bx all mesh/surface file meshes/box.stl type 1
fix inface all mesh/surface file meshes/insertion_face.stl type 1
fix wall all wall/gran model hertz tangential history surface superquadric mesh n_meshes 2 meshes bx cv
#distributions for insertion
fix pts1 all particletemplate/superquadric 15485863 atom_type 1 density constant 957 shape constant 0.0085 0.0075 0.006 blockiness constant 10.0 10.0
fix pdd1 all particledistribution/discrete 15485867 1 pts1 1.0
#region for insertion
group nve_group region reg
region bc block 0.3 0.5 -0.2 0.2 0.1 0.2 units box
#deprecated pour command
#fix ins nve_group pour/dev mass 30. 1 distributiontemplate pdd1 vol 0.25 200 massflowrate 30. vel uniform 0. 0. 0. 0. 0.0 region bc
#particle insertion
fix ins nve_group insert/stream seed 32452867 distributiontemplate pdd1 &
maxattempt 100 mass 1. massrate 1. overlapcheck yes vel constant 0. 0. -1.0&
insertion_face inface extrude_length 0.1
#apply nve integration to all particles that are inserted as single particles
fix integrator all nve/superquadric integration_scheme 1
#output settings, include total thermal energy
fix ts all check/timestep/gran 1000 0.1 0.1
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 flush yes norm no
#insert the first particles so that dump is not empty
run 1
dump dmpparticle all custom/vtk 1000 post_sugar_cube/particles_*.vtk id type mass x y z vx vy vz fx fy fz omegax &
omegay omegaz radius shapex shapey shapez quat1 quat2 quat3 quat4 blockiness1 blockiness2 tqx tqy tqz angmomx angmomy angmomz
#insert particles
run 140000 upto
unfix ins
以下のスクリプトで、非球形粒子を定義しています。constant以降の3つのパラメーターがそれぞれ$a$,$b$,$c$に該当します。blockiness constant以降の2つのパラメーターが$n_1$,$n_2$をしていしています。
fix pts1 all particletemplate/superquadric 15485863 atom_type 1 density constant 957 shape constant 0.0085 0.0075 0.006 blockiness constant 10.0 10.0
可視化
作成したinputファイルで解析を実行し、praviewで可視化を行います。
liggghtsで行った非球形粒子の解析を可視化するためにparaviewの専用ツールをインストールします。
以下リンク先からParaviewの専用プラグインをダウンロードします。大変便利なものを作っていただいているのでありがたく使わせていただきます。。。
ダウンロードした後、wslにて以下コマンドを実行し、解凍します。
tar -xzvf superquadric_visualiser_v2.0.tar_.gz
Paraviewを立ち上げてダウンロードしたプラグインを読み込みます。
Tools - Manage Plugins - Load New - LiggghtsSuperquadricVTK.xml
Tools - Manage Plugins - Load New - LiggghtsSuperquadricDUMP.xml
2つのプラグインのPropertyがLoadedになっていることを確認してください。
プラグインが読み込めたら、結果のvtkファイルを読み込みます。
読み込みが完了したら、ロードしたプラグインのフィルターを指定します。
Filters - Alphabetical - Liggghts-Superquadric visualizer(VTK)
パイプラインが追加されるので、Propertiesタブ - Aplyから非球形粒子を表示できます。
プラグインを使用して、今回実施したシミュレーションの動画を出力しました。
球形粒子を重ね合わせた場合の結果と比較したいと思いましたが、Public版のLiggghtsだと大変そうなので今回は断念しました。。。Public版は機能制限があり、開発も進んでいないので残念です。
現在は、Liggght-PFMというフォークしたプロジェクトで開発が進めれれているようなので次回はこちらについて紹介しようと思います。
