粒子法FLOSSに対する妥当性確認試験

水柱崩壊問題

Zhou et al.(1999)の実験の再現

妥当性確認の結果

まとめ

• DualSPHysicsでは
  静水圧問題が解けなかった
• OpenMpsのGC法が不自然な
  運動量生成をしていた
• OpenMpsにSPP法を導入することで
  圧力擾乱が抑制できた

数値流体力学の応用範囲

• 理学的な理論検証
• 工学的な現象解明
• 実務ものづくりの設計などへの応用（CAE）

本研究：CAEにおける数値流体力学を想定

数値流体でCAEしたい！となった時

どんな方法？

• 格子法
• 粒子法

どこから入手？

• 商用
• オープンCAE(FLOSS)

格子法

従来、よく用いられてきた

• 相互作用の相手が初期に固定
• 計算点が空間に固定されていることが多い（オイラー法）
• 研究事例が豊富
• 計算負荷は低い
• 数学的展開が容易

粒子法

近年、用いられるようになってきた

• 相互作用の相手が毎回変化
• 計算点が空間に固定されていない（ラグランジュ法）
• 移流項による誤差拡大がない
• 境界の大変形に追従しやすい
• 一般的に計算負荷は高い

粒子法をCAEで使いたい場合

自作ソルバー

一部の熟練者でないと困難なため避けたい

商用ソルバー

• Particleworks（プロメテック・ソフトウェア）
• MPS-RYUJIN（富士テクニカルリサーチ）
• Abaqus/Explicit (Dassault Systmes)

オープンCAE (FLOSS)

• DualSPHysics
• OpenMps

オープンCAEとは

オープンCAE学会などが推進

• Free/Libre Open Source Software (FLOSS)でCAEする
• ソフトウェア開発を
  利用者コミュニティーも一緒に進める

商用 vs オープンCAE

本研究の目標

• 粒子法のための妥当性確認試験を整備
• 粒子法FLOSSに対して妥当性確認を実施
• （結果を元にして、事業化の参考）

本研究の動機

顧客のソフトウェアを高速化することが主力事業

（組み込みSoCからスパコンまで、応用分野を問わず）

1. 高速化の知見を活かした自社事業ができないか？
2. 計算負荷が高い粒子法を題材にすると良さそう
3. 応用分野は海岸・海洋工学（土木工学）が良さそう
4. 元にする実装はどれがいいか？を知りたい

対象にしたいソフトウェアの特徴

• 海岸・海洋工学向きである
• ≒ニュートン流体の重力下における自由表面大変形を伴う非圧縮性流れが解ける
• FLOSSである
• 既往研究が極端に少なくない

→DualSPHysicsとOpenMPS

実施したいベンチマークの特徴

• 海岸・海洋工学で重要な物理量を計測する
• ≒水深・波圧など
• （まずは）あまり大規模・複雑ではない
• ≒基礎的な理論・実験など

→中心重力問題と水柱崩壊問題

対象ソフトウェア

DualSPHysics

• 既に多くの利用事例がある
• CUDA(GPU)対応されている（が、高速化が不十分）
• 弱圧縮SPH法なので、精度が心配

OpenMps

• 利用者が少ない（ほぼ皆無）
• 素直なC++なので分かりやすい
• 高精度MPS法が実装されているので、計算負荷が高そう

時間の都合上、詳細は割愛
（概要の第2章を参照してください）

中心重力問題とは？

• 重力が原点向き\[ \mathbf{g} \left( \mathbf{x} \right) = -g_z \frac{\mathbf{x}}{\left| \mathbf{x} \right|} \]
• 半径Rで定常になる
• 原点から半径方向に静水圧になる\[ p = \rho g_z \left( R - \left| \mathbf{x} \right| \right) \]

壁面境界が不要なので純粋な流体ソルバーの確認向き

DualSPHysicsの結果

時間が経過すると自由表面が膨張？し、圧力が低下

原因

駆動力（圧力）が依存している密度計算法の違いに起因

通常のSPH法（カーネル和型）

\[\rho_i = m \sum_{j \neq i} W_{ji}\]

DualSPHysics/WCSPH（時間発展型）

\[\frac{\mathrm{D} \rho }{\mathrm{D} t} = m \sum_{j \neq i} \left( \mathbf{u}_j - \mathbf{u}_i \right) \cdot \mathbf{\nabla} W_{ji}\]

僅かな数値誤差により拡大している

水深の時間発展

※DualSPHysicsではなく自作の簡易実装で試した結果

結論

• DualSPHysicsは静水圧問題が解けないように見える
• →浮体問題などで致命的
• 解決にはDualSPHysicsの修正が必要
• →高難易度
  （コードがあまりきれいではないので）

∴以降は、DualSPHysicsの妥当性確認は実施しない

OpenMPSでの結果

自由表面の振動が増大？

原因

勾配修正法（GC法、高精度粒子法の1つ）が原因

• \[\mathbf{\nabla} \phi_i = \frac{D}{n_0} \sum_{j \neq i} \frac{\phi_j - \phi_i}{\left| \mathbf{r}_{ji} \right|} \frac{{C_i}^{-1} \mathbf{r}_{ji}}{\left| \mathbf{r}_{ji} \right|} w_{ji} \]
• \[C_i = \frac{D}{n_0} \sum_{j \neq i} \frac{\mathbf{r}_{ji} \otimes \mathbf{r}_{ji}}{{\left| \mathbf{r}_{ji} \right|}^2} w_{ji} \]

$C_i$が異なるため、i-jの相互作用力が非対称

GC法なし

GC法なし（圧力）

結論

• GC法のせいで静水圧にならない
• 理由は、実装か理論かはまだ区別できていない

以降ではGC法を切って続行した

Oka&Koshizuka(1996)の実験

水柱崩壊（ダムブレーク）問題

原始MPS法での妥当性確認に使用されたもの

計算結果

崩壊波先端の位置

結論

• 波先端位置は実験に近い結果
• ≒自由表面の形状は良好に再現
• 原始MPS法に比べて良い結果
• →高精度粒子法による改善

Zhou et al.(1999)の実験

水柱崩壊（ダムブレーク）問題の別バージョン

定点の水深(H1, H2)・圧力(P2)を比較できる

計算結果

圧力表示

圧力の擾乱がひどく、まともに評価できていない！

原因と解決策

1. 粒子が移動する
2. わずかに空隙が発生する
3. 自由表面だと判定される
4. 圧力がゼロの境界条件が設定される

空隙を埋めるSPP法を導入すれば解決するはず

圧力表示（SPP法あり）

圧力(P2)時系列

水深(H2)時系列

水深(H1)時系列

結論

• SPP法を導入しないと、圧力擾乱が酷く使い物にならなかった
• 導入後も多少圧力擾乱があるため、引き続き改善を試みたい

まとめ

1. 粒子法FLOSSの妥当性確認試験を整備した

   • 中心重力問題
   • 水柱崩壊問題：Oka&Koshizuka(1996)
   • 水柱崩壊問題：Zhou et al.(1999)

2. DualSPHysicsとOpenMpsに対して妥当性確認をした

   • DualSPHysicsでは静水圧問題が解けなかった
   • OpenMpsのGC法には問題があった
   • OpenMpsにSPP法を導入することで圧力擾乱が抑制できた

今後の課題

• DualSPHysicsにカーネル和型密度計算の実装
• OpenMpsのGC法の問題改善
• より現実的なベンチマークの追加
  • 強制液面振動（スロッシング）
  • 造波
  • 消波
  • などなど・・・

return 0;