single-layer two-phase formulation

SPH で流体とソリッドをカップリングさせる方法は比較的古くから提案されています。

以前、PersianSPHを触っていた際に文献でおおむね理解していたのですが、実装目的であらためて詳細を確認しました。

基本は土質力学レベルです。それを丁寧に組み合わせると two-phase formulation の出来上がりです。SPH の場合はこれを2レイヤーで計算します。ここで倍の計算量。さらに相互作用を反映するため、土粒子の位置での水粒子の動き、水粒子位置での土粒子の動きを計算します。単純に4倍ではありませんが、計算量は飛躍的に増えます。

これを解決した、という文献がコチラ⇓
Two-phase fully-coupled smoothed particle hydrodynamics (SPH) model for unsaturated soils and its application to rainfall-induced slope collapse - ScienceDirect
A general smoothed particle hydrodynamics (SPH) formulation for coupled liquid flow and solid deformation in porous media - ScienceDirect
1レイヤーで2相の計算が可能に工夫されています。不飽和浸透も取り扱えるので一通りの計算ができそうです。が、境界条件が面倒。いえ、実装されたものを使用するのであればDtransu を利用するのと変わらないでしょうが、汎用性を持たせて実装しようとすると気が遠くなります。

著者には DualSPHysics の開発者が含まれています。GPU対応で開発されているのでしょう。実装されるまで待ちましょうか。

火山地帯の観測機器

先週は火山につけた観測機器のメンテナンスでした。

いくつか新しい機器を設置し、古い機器を入れ替えて帰社したのですが、あらゆるモノが錆びていました。火山ガスと雨でしょうね。防食、大事です。

PCは火山灰を吸い込み異音が発生。分解、清掃しても、治るかどうか。寿命を短くしたのは間違いなさそうです。

こういった環境では、安価な機器を使い捨て＆ばらまきにするか、お金をかけて防食し長期運用＆集中させるかのどちらかでしょう。費用が１～2桁違いますので、重要度に応じて使い分けるのが良いのでしょう。

SPHのベンチマーク

組み終わったSPHで静水圧のチェックは終わりました。残る弾塑性のチェックをどうしようか、と調べてみましたが、適当なベンチマークはなさそうでした。

流体のベンチマークは多く見かけます。が、弾塑性のものは少ないようです。FEMと同様に、室内３軸圧縮試験であったり、支持力公式との比較をいくつかの文献で見かけますが、ベンチマークとして共有されるには至っていないようです。

３軸は少し面倒ですので、支持力を試してみました。
が、文献に多くあるように、解析値のほうがやや大きく出る結果になります。うーん。ま、FEMと同程度に大きな値になっているので、あっているともいえるのですが。

新しい機能も組み込みつつ、もう少し探してみましょう。

EOS

equation of state

WCSPH で圧力を計算する際に用いられる式です。Monahan らが提案した Tait 型だそうですが、これ、簡単に見えて扱いが難しい。

係数を決めるのに最大流速を考慮しないといけないのですが、流速は計算しないとわかりません。ま、感覚で大体この程度、圧縮率は〇%程度をネラッテ音速ハコノクライ、みたいな決め方をします。DualSPHysics に使われている multi-phase （liquid-sediment）の文献では、音速はドメイン内の最大流速の10倍以上が利用されています。これは、マッハ数が0.1なので、概ね0.5%までの圧縮を許容するということでしょう。
https://research.manchester.ac.uk/en/publications/modelling-multi-phase-liquid-sediment-scour-and-resuspension-indu

では、個体が弾性から塑性になり、さらに非ニュートン流体として流動するような現象を追う場合はどうするか？

弾性、塑性では変形係数とポアソン比から音速を求められますので、それぞれ別の値を利用するという方法が考えられます。流体は水のように扱うことも可能です。
が、それぞれ別の音速を使うと計算が破綻しやすくなります。試算では塑性の計算で圧力が低くなり、塑性領域が広くなって流動化しやすくなるような結果を得ました。あと、2Dでは動くが3Dだと進まなくなるとか（2Dの計算って、問題やバグが内在しても露見しにくいことを、最近強く感じています）。

結論としては、弾性として求めた音速をそのまま利用するか、徐々に切り替えるか、なのでしょう。文献には出てこないのでノウハウの部類なのでしょうが、皆さんどうされているのか知りたいところです。