はじめに

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」の内容を書き換えて非球形粒子のシミュレーションを実行してみます。
オリジナルでは、球形粒子となっていますが立方体の粒子になるよう書き換えます。

input.liggghts

#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にて以下コマンドを実行し、解凍します。

.sh

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というフォークしたプロジェクトで開発が進めれれているようなので次回はこちらについて紹介しようと思います。