Tsubu2 （粒々）Simulator (Ver.1.0.0)
開発成果資料りゅうりゅうしんく『粒粒辛苦』 by 李伸
一つのことを成し遂げるためにこつこつと小さな努力を積み重ねること

2019/11/14
株式会社HASL
Copyright© 2010 Hyper Advanced Simulation Laboratory Co., Ltd. All Rights Reserved

りゅうりゅうかんらく『粒粒甘楽』 by
粒々のシミュレーションを楽しむこと

技術開発の狙い◆ 溶融樹脂内 Filler の配向、破砕解析（マクロな立場で）◆ スクリュ内に投入されるペレットの固体輸送/溶融可塑化挙動の新たな定量化モデルの構築を目指して◆ ポストプロセッサ（特にトレーサ粒子運動解析機能と流線表示機能）の性能向上：全ソフト共通

Copyright© 2010 Hyper Advanced Simulation Laboratory Co., Ltd. All Rights Reserved

1

技術開発の背景
FEM 1956~ 起源：航空工学
Clough, R. W., Martin, H. C., Topp, L. J., and Turner, M.J.
‘Stiffness and deflection analysis of complex structures’,
Journal of Aeronautical Sciences, 23(9),1956.

DEM 1971~ 起源：岩盤工学
Peter A. Cundall ,‘The measurement and analysis of accelerations in rock slope’,PhD thesis, Imperial College London (1971).
Computer Graphics, Games, Vessel Surgery,Cloth, Hair, Multi body dynamics, etc.
最近では、様々な分野へ適用範囲拡大中

Copyright© 2010 Hyper Advanced Simulation Laboratory Co., Ltd. All Rights Reserved2

Tsubu2 Simulator
離散要素モデル◆ Tetrahedral mass spring model*)◆ Standard DEM modelTetraTorusSphere◆ Jeffery’s ellipsoidal modelCubeDisc

Elipsoid

Fiber
SphereCopyright© 2010 Hyper Advanced Simulation Laboratory Co., Ltd. All Rights Reserved

3Cylinder
・ Arbitrary shape(defined by using 3D CAD)
・Deformable(visco-elastic flexible body)*) Mollemans, W. ,et al. ,’Tetrahedral mass spring model
for fast soft tissue deformation’, Surgery simulation andsoft tissue modeling, procededings,145(2003)

TMS model simulation samples

Simulation of a flexible ball reboundingon the floor
柔軟性（変形）の表現が可能。

Copyright© 2010 Hyper Advanced Simulation Laboratory Co., Ltd. All Rights Reserved

Dynamic simulation of multi-rigid
bodies falling on the slope floor
剛直性、モデル同士の接触の表現が可能。
剛柔併せ持つモデル。
4

DEM理論背景
運動方程式:

Newton 運動方程式1665~

dv i
mi=  (Fijn + Fijt + Fijd ) + mi gdt
j ( j ≠i )Ii

jidωi=  rij × Fijt
dt
j ( j ≠i )
粒子間相互作用力

mi : 質量v：
速度ベクトルig:

重力加速度

Ι i : 慣性モーメント

Fijn : 法線方向力ω：i 角速度ベクトル

Fijt : 接線方向力

下付き添え字は粒子番号

Copyright© 2010 Hyper Advanced Simulation Laboratory Co., Ltd. All Rights Reserved

5Fijd : 減衰力

衝突判定アルゴリズム：

Copyright© 2010 Hyper Advanced Simulation Laboratory Co., Ltd. All Rights Reserved6高分子材料の衝突運動をばねモデルで記述する場合の代表的なタイムファクター
高密度ポリエチレンヤング率E：1.5GPa
高密度ポリエチレン密度ρ：950 kg/m3
半径r(1mm)の等価ばね定数＊）：
半径r(1mm)の球の質量：

ばね振動周期:k=πrE24m = ρ π r33T=

12πm∝rkk = 2.4×106 N / mm = 4.0 × 10−6 kg

T=2.0×10-7 sec

数値シミュレーションで衝突運動状態を評価する場合、計算時間刻みは振動周期未満に設定する必要がある。100万回の計算ステップで0.2sec 未満!!!
（kを4桁下げれば、あるいはmを4桁上げれば20sec未満）*)Yamamoto, S.; Matsuoka, T. ‘A method for dynamic
simulation of rigid and flexible fibers in a flow field.’ J.
Chem. Phys, 98, 644-650(1993).
Copyright© 2010 Hyper Advanced Simulation Laboratory Co., Ltd. All Rights Reserved7ばね定数とヤング率の関係
F = kx力＝ばね定数×変位

σ = Eε応力＝ヤング率×ひずみ

Fx
=ESL

長さL π r2 S ∴F =  E  x = Ex
L  L 

半径r
力F

力F

長さ2r
半径r

半径r=πr F =Ex
 2 ???

円柱Copyright© 2010 Hyper Advanced Simulation Laboratory Co., Ltd. All Rights Reserved

k球8

力の作用面積をπr2とするのは？？？

計算高速化対策
1) 壁面情報のブロック化（ソフト的対策）
非定常計算サイクル数：Ncycle (例えば100000:105)
粒子数：Npart(例えば5000:5×103)
壁面要素数：Nwall(例えば20000:2×104)

粒子 vs. 壁面接触判定の概算演算回数：
Ncyle(105)×Npart(5×103)×Nwall(2×104)=1013回!!!
スーパーコンピュータ京（ケイ）：10ペタフロップス＝1016フロップスなら0.001秒→1エクサプロップス＝1018富岳（ふがく） 100倍高速化2021年稼働予定（理化学研究所）
スーパーコンピュータSummit(米国）：150ペタフロップス＝1.5×1017フロップスなら0.000067秒
Core i7(4コア）=50Gフロップス=5×1010フロップスで200秒

Wikipedia 情報

1ギガ＝109,1テラ=1012,1ペタ＝1015,1億＝108,1兆=1012, 1京＝1016
Copyright© 2010 Hyper Advanced Simulation Laboratory Co., Ltd. All Rights Reserved

9愚直な計算

壁面との接触計算回数：1粒子当たり
壁面要素数Nwall(20000)回

壁面情報をブロック化した計算

壁面との接触演算回数：1粒子当たりブロックの検出計算回数100回＋ブロック内
壁面要素接触判定計算回数(200)回＝300回

演算回数を300/20000=3/200=0.015倍に低減
ブロック 12100
3
各ブロック毎に壁面要素情報(200要素分)を格納

Copyright© 2010 Hyper Advanced Simulation Laboratory Co., Ltd. All Rights Reserved10Z軸周りに回転するスクリュモデルをZ方向に対してブロック化するのは特に有効（非定常サイクル毎にブロック情報を更新する必要が無い）

2) GPU利用（ハードウェア的対策）
理論性能上は約100倍速い（人気ゲームのスパーチャージャー）
2019/2/22販売開始 価格4万円前後のミドルクラス

NVIDIA® GeForce GTX® 1660Ti

5.437 TFLOPS(単精度演算理論性能)
Core i7(4コア）53.28 GFLOPS(理論性能) ：Wikipedia 情報
GPU 対応Fortranコンパイラー導入検討中

株式会社ソフテック様資料より引用
https://www.softek.co.jp/SPG/Pgi/product.htmlCopyright© 2010 Hyper Advanced Simulation Laboratory Co., Ltd. All Rights Reserved

11Jefferyモデル*)理論背景z

せん断流動場内の楕円体の回転運動を記述する方程式：

 ∂w  ∂v (b2 + c2 )ω1 = b2   − c2   , ∂z  ∂y  ∂u  ∂w (c2 + a 2 )ω2 = c2   − a 2   , ∂z  ∂x cabx ∂v  2  ∂u (a + b )ω3 = a   − b   ∂x  ∂y 2y:粒子座標系( u, v , w ):粒子座標系流速ベクトル(ω1 , ω 2 , ω 3 ) :粒子座標系角速度ベクトル( a , b, c ):楕円体の三軸半径22( x, y, z )*)G.B.Jeffery, ‘The motion of ellipsoidal particles immersed ina viscous fluid’, Proc. Royal Soc. London. Ser. A 102,161(1922)

https://royalsocietypublishing.org/doi/pdf/10.1098/rspa.1922.0078にて閲覧可能
Copyright© 2010 Hyper Advanced Simulation Laboratory Co., Ltd. All Rights Reserved12

Time derivative formula for quaternion：

ウイリアム
Quaternion(四元数) 1843年アイルランドの数学/物理学者・ローワン・ハミルトンが偶然発見

Jeffery model•  q0  − q1ω 1 − q 2ω 2 − q 3ω 3 • •qqq
ω
−ω+ω11q

0 1
3 2
2 3 q =  •1 = ω q =   22 q0ω 2 − q1ω 3 + q 3ω 1 q
 2 q ω −q ω +q ω • 2 11 2  0 3
 q3 

Leap-frog time integration

Particle (material) coordinate ⇔ Physical (Cartesian) coordinate transform formula：
 ∂u
 ∂x
 ∂v
 ∂x
 ∂w
 ∂x


∂u
∂y∂v
∂y∂w
∂y∂u  ∂U∂z  ∂X∂v ∂V= R ∂X∂z  ∂W∂w  ∂X∂z ∂U∂Y∂V∂Y∂W
∂Y

∂U 
∂Z 
 q02 + q12 − q22 − q32∂V  TR , R =  2(q1q2 + q0 q3 )∂Z ∂W  2(q1q3 − q0 q2 )∂Z 

( X ,Y , Z )(U , V , W )Copyright© 2010 Hyper Advanced Simulation Laboratory Co., Ltd. All Rights Reserved132(q1q2 − q0 q3 )q02 − q12 + q22 − q322(q2 q3 + q0 q1 )2(q1q3 + q0 q2 ) 2(q2 q3 − q0 q1 ) 
q02 − q12 − q22 + q32 :物理座標系:物理座標系流速ベクトルTsubu2 Simulator 運用方法
プリプロセッシング
1）離散要素生成機能
1-1）中心座標、半径（楕円球の場合は、3軸半径）指定による生成

1-2）矩形領域内、円形領域内生成

矩形、円形の選択
要素間距離と領域径の指定

中心座標
半径Copyright© 2010 Hyper Advanced Simulation Laboratory Co., Ltd. All Rights Reserved

141-3）テトラソリッド要素(3D FEM解析モデル）の節点を離散要素粒子に変換(TMSモデル）

3D-CAD

AutoCAD,Solidworks,
FreeCAD,etc.

Femap,Netgen neutral file インポート機能

IGES file

3D Tetra Mesher

Femap,Netgen,etc.

Mesh file

Tsubu2 simulator

Copyright© 2010 Hyper Advanced Simulation Laboratory Co., Ltd. All Rights Reserved

15Copyright© 2010 Hyper Advanced Simulation Laboratory Co., Ltd. All Rights Reserved

162）壁面情報生成/インポート機能
2-1） CAD 作成STLファイル
インポート機能
STL file インポート機能

3D-CAD

AutoCAD,Solidworks,
FreeCAD,etc.

STL file

Tsubu2 simulator

Copyright© 2010 Hyper Advanced Simulation Laboratory Co., Ltd. All Rights Reserved

172-2） Single Screw Simulator, Twin Screw Simulator テンプレートで作成した2.5D FEM解析用モデルの転用
Single Screw Simulator

Twin Screw Simulator

Copyright© 2010 Hyper Advanced Simulation Laboratory Co., Ltd. All Rights Reserved

18アナリーシス

解析種別

解析内容

関連ソフトウェア

DEM

標準の球形離散要素解析

SSS,TSS,FS3D,
Multi purpose

Tetrahedral mass spring

非球形離散要素解析
Filler の配向/破砕解析

SSS,TSS,FS3D,
Multi purpose

Jeffery(Lagrangian)

Jeffery回転楕円体の配向解析（Lagrangian 描像）

Jeffery(Eulerian)

Jeffery回転楕円体の配向解析（Eulerian 描像）

Flow tracer

従来の粒子運動解析（流線描画）

SSS,TSS,FS3D,
Flat simulator,
Spiral simulator

SSS : Single Screw Simulator
TSS : Twin Screw Simulator
FS3D : Flow Simulator 3D
Copyright© 2010 Hyper Advanced Simulation Laboratory Co., Ltd. All Rights Reserved

19

Copyright© 2010 Hyper Advanced Simulation Laboratory Co., Ltd. All Rights Reserved20ポストプロセッシング
Bmp file output チェックボックスをチェック状態とすると指定フォルダー内に描画
情報をBmpファイル形式で保存

・壁面（移動境界）アニメーション表示
・離散要素（色付き）アニメーション表示
・各種履歴情報グラフ図表示
・表示結果Bmpファイル出力・表示結果Bmpファイル再生表示

Copyright© 2010 Hyper Advanced Simulation Laboratory Co., Ltd. All Rights Reserved

21フォルダー内に複数のBmpファイルが自動生成①Tsubu2 simulator で再生②Power point に張り付け
プレゼン資料作成

Copyright© 2010 Hyper Advanced Simulation Laboratory Co., Ltd. All Rights Reserved

22③Title, Caption 等を付けて動画ファイル作成

Test samples
Feed simulation for a single screw extruder

粒子直径3mm・単軸スクリュ形状は充満搬送され易い薄肉流路形状。
・Darnell & Mol固体輸送モデルでは、固体の挙動を一定速の栓流扱いするが、実際は離散要素の移動速度に分布が生じる。
・ひずみ速度（肉厚方向の速度勾配）を正確に表現するのは至難の業。
Copyright© 2010 Hyper Advanced Simulation Laboratory Co., Ltd. All Rights Reserved

23

Feed simulation for a corotating self wipe twin screw extruder
Transport direction

Upperside view

Lowerside view

・2軸スクリュは未充満搬送され易い厚肉流路形状。
・未充満状態では、重力の影響を受け易く、底面側の充満状態が高い。
・天面側では、噛合領域にペレットの溜が発生。
・予想通り、ニーディングや逆ねじ領域は、充満率が高く、滞留時間も長くなり、混練が促進される様子が視覚化されている。

Copyright© 2010 Hyper Advanced Simulation Laboratory Co., Ltd. All Rights Reserved

24Breakage simulation for a spherical shape filler in viscous simple shear flowyu f = γɺ y : 周囲流体のx方向流速
vf = 0 :周囲流体のy方向流速
Filler diameter : 6 mm
Strain rate γɺ : 50 /sec
Particle diameter : 0.5 mm
Fluid viscosity η : 100 Pa・sx
Breakage criteria: 1.5Copyright© 2010 Hyper Advanced Simulation Laboratory Co., Ltd. All Rights Reserved25単純せん断流動内の半径 r の剛体球の回転速度、回転数は？

y球の両極で回転速度が流体速度と一致していたとすると球の1秒当たりの回転数Nは、剛体球の半径に依存せず、

us = γɺr = u fγɺ :N=γɺrγɺ=2π r 2πx50 /sec の場合、約8回転

us = −γɺr = u f

Copyright© 2010 Hyper Advanced Simulation Laboratory Co., Ltd. All Rights Reserved26

y

この点の球表面上の速度：

us = γɺr = u f

us = γɺr sin θ = u f , 不思議なことにx方向の分速度は流体速と一致。
vs = −γɺr cos θ ≠ v f = 0 が生じるが、周囲流体のv は0．この剛体球表面では下向きの分速度成分f

速度差は、粘性抵抗力となり剛体球の回転運動を阻害する。

θγɺ
∴N <
2π

θ

x
回転を阻害する
粘性抵抗力

Jeffery model

us = −γɺr = u fb aN=γɺ≤γɺ
,4π1 
2π  AR +AR bAR = : Aspect ratioaNAR
アスペクト比を増加させると回転数は減少する（配向し易くなる）。
Copyright© 2010 Hyper Advanced Simulation Laboratory Co., Ltd. All Rights Reserved27

Mass Spring Constant : 100 N/m

・TMSモデルで剛体球の回転運動を表現可能。
・0.28秒で約1回転強、1秒で約4回転：試算値（約8回転以内）。
Jeffery model から評価される回転数（約4回）と一致。
Copyright© 2010 Hyper Advanced Simulation Laboratory Co., Ltd. All Rights Reserved28

Mass Spring Constant : 10 N/m

・TMSモデルで柔軟体の変形表現可能。
・斜めに配向した変形回転運動。

Copyright© 2010 Hyper Advanced Simulation Laboratory Co., Ltd. All Rights Reserved29Mass Spring Constant : 2 N/m

・TMSモデルで柔軟体の破砕表現可能。
・回転しながらの連鎖的破砕挙動。

Copyright© 2010 Hyper Advanced Simulation Laboratory Co., Ltd. All Rights Reserved30Orientation simulation for Jeffery’s ellipsoidal modelin viscous shear & elongational flow
MPa1Pressure0mm/s
8000

FS3D流動解析1秒以内
Velocity
Copyright© 2010 Hyper Advanced Simulation Laboratory Co., Ltd. All Rights Reserved310Aspect ratio :a/b/c:2.5/3/2.5
cba
Aspect ratio :a/b/c:1/3/1c・Jeffery model 解析時間約20秒(2500サイクル)b・解析結果は上下鏡面対称（当然）。
上部では反時計回り、下部では時計回りの回転運動。流動中央部では並進運動。

aAspect ratio :a/b/c:0.5/3/0.5c

・アスペクト比:b/aの増加に伴って回転周期が増加（回転し難く、配向し易くなる）。

ba
Copyright© 2010 Hyper Advanced Simulation Laboratory Co., Ltd. All Rights Reserved

32

Deformation & breakage simulation for a spherical shape filler
in viscous shear & elongational flow

Copyright© 2010 Hyper Advanced Simulation Laboratory Co., Ltd. All Rights Reserved33伸長流動による変形並進運動
伸長流動による（流動方向に延伸） 延伸促進
破砕開始

破砕後の並進運動

流速差を反映した破砕要素の凝縮

分散性能〇
分配性能×Copyright© 2010 Hyper Advanced Simulation Laboratory Co., Ltd. All Rights Reserved34

今後の取り組み
Single simulator, Twin screw simulator
・改良成果は、随時プレリリース
・正規版リリース2020/4月予定
Tsubu2 simulator
・高速化対策、ハード/GPUコンパラー導入後、出荷準備スタート
・正規版リリース2020/4月予定
Flat simulator, Spiral simulator
・ 2.5D FEM 多層流動解析機能のブラシュアップ
・検証解析継続、成果の学会誌投稿予定(2020/3月)
・正規版リリース2020/6月予定
Hyper blow
・ 高速化対策(金型接触判定計算）
・ 要素ライブラリー追加（軸対称ソリッド要素）
・ ポストプロセッサ機能改良
・ 正規版リリース2020/6月予定

Copyright© 2010 Hyper Advanced Simulation Laboratory Co., Ltd. All Rights Reserved35