2024年7月3日水曜日

double-layer two-phase formulation

SPHにおいて、two-phase を single-layer で実装するか、double-layerにするか?

前者は計算が速い!
でも、水の流れは表現できません。(工夫すればできるのかもしれませんが、single ではなくなります。)

後者は計算が遅い!しかも圧倒的に。
でも、水の流れは表現できます。

考えた結果、もともと計算してみたかった「崩土から水の抜ける表現」が可能な後者を選択しました。既存のプログラムを容易に変更できそうな点も good です。

2相 を double-layer で解く文献はたくさんあります。が、なぜか式が統一されていません。以下が最も詳しく書かれていましたので、これを参考に実装します。
A coupled fluid-solid SPH approach to modelling flow through deformable porous media - ScienceDirect

式を整理しながら不飽和の表現を single-layer の文献からいただきました。これで、飽和・不飽和浸透と力の関係を解けるはずです。

コツコツ、組んでみましょう。

**********************************
20240707追記
地下水の上下で飽和度(水の単体)を変更すると、境界で爆発してしまいます。計算過程を追えば当然なのですが、土層境界よりも挙動が激しいようでした。飽和度を低くするには文献のように水粒子をまばらに配置するしかなさそうです。
似た現象ですが、地表流がある直下で飽和度が過剰に高くなります。透水係数が実際よりも高く計算され、水が入りやすくなります。なかなか難しい。
SPHでは密度変化がないか小さい現象(弱圧縮)を扱うのが前提なので、単体の異なる土層とか飽和‐不飽和‐地表流などを扱うのが難しいのでしょうね。

2024年6月26日水曜日

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対応で開発されているのでしょう。実装されるまで待ちましょうか。

2024年6月23日日曜日

火山地帯の観測機器

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

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

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

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


2024年6月16日日曜日

SPHのベンチマーク

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

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

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

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

2024年6月8日土曜日

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の計算って、問題やバグが内在しても露見しにくいことを、最近強く感じています)。

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


2024年5月30日木曜日

弾塑性SPH その2

新たに頂いたコードを参考に、以下の組み合わせで計算ができるように改造しました。

・弾性
・弾塑性
・弾性‐流体
・弾塑性‐流体
・流体

これらのほとんどは日本語の図書でも紹介されています。個体と流体を同じように実装できるのはSPHならではでしょう。

弾塑性(一部流体)の挙動はコチラ。端から土圧により崩壊し始めます。

少しc・φをあげると、崩れにくくなります。その場合、まず角が楔のように崩れ、次に長辺の2/3が崩れ、最後に残りが崩れるような時間差が見られました。ちなみに、前回のアップしているイメージは「弾性‐流体」の計算結果で、全体がペシャっと潰れるような挙動を示しています。これだけ自由な表現ができると、多様な自然現象を再現することが可能でしょう(予測は相変わらず困難ですが)。

まだ、いくつか疑問点が残っています。資料を探して早く完成させましょう。

2024年5月26日日曜日

壁粒子の圧力

SPHで壁際の粒子が壁に張り付く現象を、後回しにしていました。
ある程度動くようになったので取り組んでみることに。

まずは mDBC で様子見。
効果はありますが、残念ながら完全ではありません。いくつかの粒子は壁に残り、計算時間も余計にかかりました。

負圧が効いているのは想像がつくのですが、他の修正方法を思い付きません。
そこで googling。
すると、「まさにコレ!」といった回答を見つけました。条件も全く同じです。
Trying my chances here; anyone knows why particles stick to the wall in my SPH simulation? - Specific Domains / Modelling & Simulations - Julia Programming Language (julialang.org)

Here the simplest fix is simply to put the boundary particle density to reference density, if it ever gets below that reference value.
In this way a wall can never produce suction, but only “push”.

当たり前すぎて一般的には解説されていないのでしょうね。初心者には思いつかないノウハウです。
密度でなく圧力を変更。数行の追加で計算時間もかからず、簡単に解決できました。