以前の記事で、DEMにおいて非球形粒子を用いた解析を行う方法を2つ紹介しました。

球形粒子を重ね合わせることで再現
超二次関数を用いて形状を表現

その記事では、超二次関数を用いる方法を紹介したので今回は球形粒子を重ね合わせる方法を紹介します。

粒子を重ね合わせる方法は、Liggghts-PublicではなくLiggghts-PFMというバージョンを使用した方が簡単なので、今回はLiggghts-PFMを使用した方法を紹介します。Liggghts-PFMについては、以下を参照ください。

全体の流れ

以下のように2stepに分けて非球形粒子を生成します。Liggghtsでは複数の粒子を重ね合わせて取り扱うことをclumpといいます。

Step1 解析したい粒子形状の形状ファイル作成(stl.vtkなど)
create_multisphere_clumpコマンドで対象形状に合わせ球形粒子をclump。その結果を出力(座標データ)
Step2 作成したデータを用いてmultisphereとして定義し、非球形粒子を生成。
単一粒子の場合と同じように取り扱い、任意個数非球形粒子を生成し、解析。

Liggghts-Public版は、Step2のみであれば実行可能です。multispereを構成する各粒子の座標位置を入力したファイルを手作業で作成すれば解析可能ですが、ちょっと現実的ではありませんね。。
以降、チュートリアルを用いて実行方法を確認してみたいと思います。
チュートリアルは、インストール先 > Liggghts-PFM > example > LIGGGHTS > PFM > multishare > Create_clump に保存してあるものを使用します。

Step1

チュートリアルフォルダは以下の構成になっています。
このチュートリアルでは２種類の非球形粒子をclumpで生成し、解析する事例になっています。

create_clump/
        ┣  data/                     #clump後の結果データ  
        ┃
        ┣  meshes/                   #対象の形状データ
        ┃     ┝ infile1.vtk               
        ┃     └ infile2.vtk                
        ┃ 
        ┣  post/                     #解析結果
        ┃ 
        ┣  in.create_clump　　　　　　#step1用inputファイル
        ┗  in.multisphere            #step2用inputファイル

Step1用のinputfileの中身を確認してみましょう。

in.create_clump

# Multisphere clump generation from surface mesh

create_multisphere_clump dmin absolute 0.0002 rmin 0.002 pmax 1.0 seed 2133421 &
                         surfacefile meshes/infile1.vtk &
                         subdivide linear 3 subdivisionfile meshes/subdiv1.vtk &
                         clumpfile data/clump1.dat post/clump1.vtk

create_multisphere_clump dmin absolute 0.0002 rmin 0.002 pmax 1.0 seed 4524661 &
                         surfacefile meshes/infile2.vtk &
                         subdivide linear 3 subdivisionfile meshes/subdiv2.vtk &
                         clumpfile data/clump2.dat post/clump2.vtk

clumpを実行するだけの簡単なファイルになっています。create_multisphere_clumpを2回実行しており、２種類の非球形粒子を生成していることが分かります。
create_multisphere_clumpの引数は次のような構文になっています。

create_multisphere_clump dmin dstyle dvalue rmin rvalue pmax pvalue seed svalue surfacefile file keyword values filetype fileargs

引数はそれぞれ以下のような意味になります。
2つのcreate_multisphere_clumpコマンドで異なっているのは、形状ファイルとsvalueの値となります。
２つの形状がどのようになっているか確認してみます。

dmin = 必須ワード
dstyle = absolute / radius_ratio
dvalue = 粒子と粒子生成に使用される次のサーフェスポイントとの間の最小距離
rmin =必須ワード
rvalue = 最小粒子半径
pmax = 必須ワード
pvalue = 球体生成に使用されるサーフェスポイントの最大割合
seed = 必須ワード
svalue = 粒子生成のシードポイントの数
surfacefile = 必須ワード
file = 形状ファイルの名前
keyword = invert_normals / subdivide #サーフェス再分割の設定
invert_normals value = yes / no
subdivide values = style N #再分割法に関する設定
style = linear or loop or butterfly
N = number of subdivision iterations
filetype = clumpfile / datafile #再分割結果の出力ファイル形式
clumpfile args = 再分割結果のファイル名
datafile args = filename / atom_type / density
filename = ファイル名　#クランプ結果
atom_type = 原子タイプ (オプション)
density = 粒子密度 (オプション)

形状ファイルをParaviewで開くと以下のような形状となっています。
ファイルは異なりますが同じ形状になっていることが分かりました。つまり2つのコマンドで異なっているのはsvalueのみということになります。では、実際に実行してみてどのようなclump結果になるか確認します。

以下のスクリプトで実行します。

.sh

liggghts -in in.multisphare

するとdataとpostフォルダに出力結果が保存されます。

create_clump/
        ┣  data/                     
        ┃     ┝ clump1.dat           #clump後の座標データ    
        ┃     └ clump2.dat           #clump後の結果データ
        ┃
        ┣  meshes/                   
        ┃     ┝ infile1.vtk 
        ┃     ┝ infile2.vtk  
        ┃     ┝ subdiv1.vtk          #細分化後の形状ファイル
        ┃     └ subdiv2.vtk          #細分化後の形状ファイル              
        ┃ 
        ┣  post/                     
        ┃     ┝ clump1.vtk           #可視化用vtk    
        ┃     └ clump2.vtk           #可視化用vtk
        ┃ 
        ┣  in.create_clump　　　　　　#step1用inputファイル
        ┗  in.multisphere            #step用inputファイル

それぞれの可視化用ファイルをParaviewで読み込んで、clumpによって生成された非球形粒子を確認してみます。
元の形状ファイルは同じですがseedの設定によって、結果が変わっていることが分かります。シード数が多かったclump2の方がより再現性が高くなっています。

Step2

step1で生成したclumpデータを使用して、非球形粒子を定義し、解析を実行します。
sep2用のin.multisphereファイルを確認します。

in.multishere

# Multisphere example using 2 different clumps

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.35 0.35 -0.35 0.35 0.0 0.5 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 104 ntry 1000000 spheres file data/clump1.dat scale 1.0 type 1

fix pts2 all particletemplate/multisphere 4573687 atom_type 1 density constant 2500 &
    nspheres 150 ntry 1000000 spheres file data/clump2.dat scale 1.0 type 2

fix pdd1 all particledistribution/discrete 521451 2 pts1 0.5 pts2 0.5

# region and insertion
region  bc cylinder z 0.0 0.0 0.15 0.05 0.45 units box

fix ins all insert/pack seed 100001 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 1000 volumefraction_region 0.15

# 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/msclump*.vtk id type mol x y z vx vy vz fx fy fz omegax omegay omegaz radius

run     50000 upto

先ほど作成したclumpデータから非球形粒子を定義し、その粒子をある領域内に発生させ落下させるシミュレーションを実行しています。
非球形粒子の定義は以下のコマンドで実行されています。

# templates and distribution for insertion
fix pts1 all particletemplate/multisphere 1272435 atom_type 1 density constant 2500 &
    nspheres 104 ntry 1000000 spheres file data/clump1.dat scale 1.0 type 1

このファイルを以下コマンドで実行します。

.sh

liggghts -in in.multisphare

実行するとpostフォルダー内にvtkファイルが出力されますので、Paraveiwで読み込んで可視化すると以下のような動画になります。
粒子を重ね合わせる方法で非球形粒子を取り扱うことができています。

次回は、超二次関数を使用した場合との比較結果を共有しようと思います。