文献：DualSPHysics with Chrono

 nvestigation of subaerial landslide–tsunamis generated by different mass movement types using smoothed particle hydrodynamics - ScienceDirect

AI要約

本研究は、陸上起源の落石・転倒（トップリング）など異なる質量移動タイプ（MMT）が引き起こすサブエリアル地すべり津波（SLT）の差異を、SPH法（DualSPHysics v5.0.4）とProject Chronoの連成により数値的に比較・評価したものである。従来は滑動（sliding）に偏った研究が多く、落下（falling）・転倒（overturning）の波特性は未解明部分が多かった。本研究は計26ケース（長方体・円柱ブロック、岩石想定密度ρs=2700 kg/m3）で、最大波高・振幅、減衰、方向分布などを整理し、無次元パラメータ（Froude数F、相対厚S=s/h、相対長L=l/h、相対幅B=b/h、相対質量M）の関数として新たな経験式を構築した。

土砂（水陸質量）と水の連成方法（数値手法・境界・結合）
連成の基本枠組み:
オープンソースのDualSPHysics v5.0.4（WCSPH）を用いて流体（自由表面水）をSPH粒子で離散化し、反転型の剛体運動はマルチフィジクスエンジンProject Chronoで剛体力学を解く弱結合（loose coupling）。SPH側で算出した流体力をChronoに受け渡し、Chronoが剛体の平動・回転加速度を更新し、これをSPHにフィードバックする時刻歴連成。ChronoはDVI（微分変分不等式）に基づく非平滑マルチボディタイムステッピングで接触・ヒンジ拘束を解く。両ソルバは独立タイムステップで、1 SPHステップ中にChronoが複数サブステップを持ちうる。データ連携はDSPHChronoLibインターフェース。落下型はChrono拘束を用いず重力自由落下。

流体方程式・数値安定化:
ナビエ–ストークスのSPH離散。弱圧縮性（WCSPH）で状態方程式を用い、音速比で圧縮性を制御（c0=71.2 m/s、密度変動<1%を確保）。数値粘性は人工粘性（α=0.01）を採用。

境界条件:
固体・壁はmDBC（modified Dynamic Boundary Condition）を採用。固定境界を静止流体粒子として扱い、ゴースト粒子の線形外挿で密度場を平滑化、壁近傍の圧力安定性と隙間抑制を図る。圧力場安定のため密度拡散項も併用。

計算設定例:
粒子間隔dp=12 mm、CFL=0.2、音速係数17、GPU（RTX A5000）で4秒計算に約40時間（約8300万粒子）。検証では氷密度相当、主試験では岩石を想定して剛体密度ρs=2700 kg/m3を使用。

検証・収束:
大規模実験（氷ブロックの落下・反転）に対し一次波峰a1は±10%、一次波高H1は±15%以内で整合。空中スプラッシュの影響が小さい測点ほど良好。空間収束は反転型で難しさが残るがdp≤12 mmで一次波は概ね良好。時間収束もCFL 0.05–0.20でaMの差は最大約1.3%。

結論
転倒は落下よりも大きいaM, HMを生じ、ただし沿軸減衰は速い。
新しい経験式（発生・伝播）はF, S, B, L, r/h, 方向γの関数で±50%程度の精度帯で整合。
実事例（カピトリオ）での適用は、適用範囲逸脱や地形・運動複合性により誤差が増えるが、初期評価に有効。
滑動との関係式により、MMT間の波規模変換の目安を提示。今後は粒状体、緩斜面、複合運動、適用範囲拡大が課題。

2025年。DualSPHysics で Chrono との連成は以前よりサポートされていますが、このようなことができるのですね。

文献：single-layer multi-phase formulation

 Modelling internal erosion using 2D smoothed particle hydrodynamics (SPH) - ScienceDirect

AI要約

1. 研究の背景と従来の課題
内部浸食現象のモデリングにおいて、従来の連続体モデルは微細粒子の進化や固相変形を十分に考慮できていなかった。また、有限差分法（FDM）や有限要素法（FEM）のようなメッシュベース手法は、大変形を伴う場合にメッシュの歪み問題に直面し、全破壊プロセスの解析が困難である。近年、メッシュフリー手法であるSPHが内部浸食の解析に適用され、ある程度の成功を収めているが、既存のモデルには詳細な説明の不足、検証の欠如、材料強度への浸食による影響の不考慮、数値安定化の課題などが残されていた。

2. 本研究の提案手法
本研究では、以下の点を改善した新しいモデルを提案。

• 飽和/不飽和条件および浸透誘発浸食、輸送、堆積を記述するためのu-w-p定式化に基づく3相5成分SPHモデル。計算は single-layer multi-phase。
• 微細粒子含有量と材料強度間の非線形関係を考慮するために修正された「Ma et al. (2022)」の微細粒子依存構成モデル。
• 浸透-浸食-変形連成プロセスに適応するように拡張された境界処理および拡散アルゴリズム。数値シミュレーションをさらに安定化させるための新しい粘性散逸項。

3. モデルの基礎と構成要素
浸食性多孔質材料を、固相（不浸食性固体骨格と浸食性固体微細粒子）、液相（液体化浸食微細粒子と間隙水）、気相（空気）の5成分からなる多相多成分連続体として扱う。各相/成分の体積分率が定義され、これにより各相の物理量を定量化する 。
支配方程式は、質量保存則と運動量保存則に基づいて定式化され、質量バランス方程式は固形骨格、微細粒子含有量、液体化浸食微細粒子、水密度について記述される。全応力はBishopの有効応力概念を用いて計算する。
構成則は以下のとおり。

• 浸食/堆積法則: 浸食と堆積のメカニズムを記述するため、「Cividini and Gioda (2004)」による浸食項を修正し、実験結果に合わせた経験的関係を利用。
• 土の構成モデル: 弾塑性理論とDrucker-Prager降伏基準が適用され、材料強度パラメータ（粘着力と摩擦角）は、ひずみ軟化、吸引圧による硬化、および微細粒子含有量への依存性が考慮される。特に、微細粒子含有量と材料強度の非線形関係を組み込むために「Ma et al. (2022)」のモデルが修正されており、浸食による強度低下が表現されている。
• 水の構成モデル: 間隙水圧は弱圧縮性流体の状態方程式から計算され、飽和度と吸引圧の関係にはVan Genuchtenモデルが用いられる。透水係数は飽和度とKozeny-Carman式に基づく空隙率の変化に応じて変化する。

2024年の発表です。実務で利用することはないですが、2Dとは言え他国はここまで進んでいます。キーワードに DualSPHysics が載っていますので、いずれ実装されるのでしょう。私は2相2成分で止まっていますが、何とかしないといけないと思わされます。

文献：SSRFEM + SPH

Post-Collapse Evolution of a Rapid Landslide from Sequential Analysis with FE and SPH-Based Models

AI要約

研究対象・地形・地質
Sant’Andrea地すべりは北東イタリア（アルプスのPerarolo di Cadore）に位置し、ボイテ川の左岸に隣接している。崩壊時に河川を堰き止めて下流集落に影響を及ぼす可能性が高い区域。
旧地すべり区域を含む72,000m²の範囲。主な地層は、厚さ30mのデブリ（A/B/C層）と、その下部に石膏・硬石膏を主とする岩盤（D層）。土壌・地質調査や地表変位のモニタリング（RTS等）が継続的に行われている。

FEM（有限要素法）解析
3D-FEM（Midas GTSNX）で現地地形・地質情報をもとにモデル化。地層・水理条件（降雨時の水位など）を詳細に考慮。
FEM解析では、剪断強度低減法（SRM: Strength Reduction Method）を使い、滑落面と不安定体積を評価。
乾燥状態・湿潤（雨後）状態でそれぞれ解析を実施し、「最悪ケース」の水理条件で土壌強度を低減した。

不安定体積（Unstable Volume）の詳細
地表変位分布から「iso-displacement surface（等変位面）」を用いて安定領域と不安定領域を識別。「地表変位の最大値の10%」を閾値として不安定体積を定義。
滑落面・分布は地質・水文条件で大きく変化。乾燥時は深く広範囲、湿潤時は浅く面積が狭い傾向。乾燥時の方が多いが、実際には湿潤時の方が危険度が高い（安全率が下がり引き金となるため）。

10%の根拠
・変位分布と実測データの適合性：モニタリングデータとモデル計算結果の分布が一致するように調整。

• 7.5%〜15%の範囲を候補とし、その中で安全側かつ現実的な値として10%を採用。
• 最大せん断ひずみゾーンと閾値により抽出された領域との整合性を確認し、閾値設定の妥当性を支持。
• 10%はより保守的（安全側）な設定であり、過小評価を防ぐ目的で選定。

SPH（Smoothed Particle Hydrodynamics）による伝播・堆積予測
上記FEMで得た不安定体積をSPHモデルに入力し、土砂流の伝播距離・堆積分布を解析。浅水方程式＋Völlmyモデル（粘性・摩擦抵抗の両効果）を採用。特に湿潤時には流動性が高まり、土砂がボイテ川谷でダム状堆積（厚さ15-20m）となった。摩擦角の値（μ=0.45, 0.58）などの設定も安全側・リスク評価に寄与。

感度解析の記載
「摩擦係数（μ）は0.36～0.62、摩擦角ϕとしては20度～32度の範囲で変化させて計算し、土砂の流下距離・堆積の広がりにどの程度影響するかを検証」。具体的には、摩擦角（ϕ）を20〜32度で7パターンに設定し、ξ＝500ms^-2で各パターンをシミュレーション。
結果として、「runout distance（流下距離）は摩擦角の値に敏感に反応したが、乱流係数（ξ）の影響は小さかった。

解析結果の詳細
乾燥時は不安定体積が多いが、流動性が低いため堆積量・範囲は限定的。湿潤時は不安定体積はやや少ないが、流動距離（runout）が長く、河川全体を埋める可能性がある。

結論・意義
複数の水理条件下で不安定体積と安全率（FS）を評価。
湿潤時（大雨後）が最も危険度が高く、比較的小さい体積でも流動性増加により下流域のリスクが顕在化。SPHモデルの摩擦係数（μまたは摩擦角ϕ）は、土砂流動・停止位置・堆積形状に対して極めて高い感度を示し、「runout distance」に大きく影響する。

SSRMを利用して不安定領域を求め、SPHで流す２段階の手法です。2021年です。私はSPHのみで解けるように組みましたが、発表するとしても2027年。「摩擦係数（μまたは摩擦角ϕ）は、土砂流動・停止位置・堆積形状に対して極めて高い感度を示す」という結果も同じ。のんびり組んでいたら時代遅れになっちゃいました。 

SPH on GPU

DualSPHysics が v5.4 になって数カ月が経過しました。

このVerでは Variable resolution model が実装されています。これはありがたい（と思いながら理屈を読んでいません）。
取り急ぎ解きたい問題がないので、まだ触っていません。が、定期的に新しく発表された文献の内容が取り入れられ保守されているのは、良いコード、チームだと感じます。

そういえば手元の 弾塑性 SPH コードが OpenMP 止まりだったので、そろそろ GPU を使おうと思い立ちました。FP64を潤沢に使える環境ではないので、単精度でも問題ない箇所のみを OpenACC で変更します。Cell-linked list  の作成なら単精度で問題ありません。ドメインが大きくなければ粒子同士の距離計算も問題にならないそうですが、これは完全に理解していないのでひとまず対象外としました。で、そこをOpenACC化。RTX8000ですので、CC75でコンパイル。

動きました。約3倍の速度。まだ効率が悪いのでしょう。が、10日の計算が3,4日で終わるなら上等です。
ひとまず様子を見ましょう。

Post-failure analysis

崩壊後の土砂流動を再現/予測するシミュレーション手法は、複数あります。

個人的には2次元平面‐浅水方程式（LSFLOW）、3次元DEM（PFC）、3次元SPH‐弾塑性＋ビンガム流体（自作）を利用したことがあります。

2次元平面の計算が軽く、3次元が重いのは仕方がありません。再現性はどれも同じでしょうか。

LSFLOWが流体っぽく対岸への乗り上げまで再現できていましたが、天然ダムの形成を再現するにはパラメーターを場所ごとに細かく設定する必要があり、予測には不向きです。

SPHは弾塑性ベース、pre-failure stage からスタートするので、土塊がまとまりながら落ちる（最初から広がらない）という実際の痕跡に近い挙動を示します。すべり面の発達過程や位置がわかる点で、他の2つよりも有利でしょうか。

DEMもSPHも、LSFLOWほどパラメータに敏感ではない印象があり、予測でも取り扱いやすいのが良い点でしょう。

以下では CEL: Coupled Eulerian-Lagrangian （ABAQUS）、MPM（ANURA3D ）、SPH（GEOXPM） を比較されています。意外と見かけない、ありがたい文献です。

Post-failure stage analysis of flow-type landslides using different numerical techniques - ScienceDirect

再現性優先でパラメータスタディするのではなく、パラメータ固定で再現性と速度、ファイルサイズを比較しています。MPIでも複数のノードを使わない条件での比較ですのでコードの能力を十分に発揮していないのですが、それでもSPHが速くて良いというのは予想と逆の結果でした。うーん、他の手法はそんなに重いのか？

もう少し情報が欲しいですね。

double-layer two-phase formulation その3

double-layer two-phase SPHの実装が概ね完了しました。

実装を通し、理論と実際を理解できました。
現状、2点の問題があるようです。

１つ目は seepage force (drag force) の計算式。

これまでの文献にも、DualSPHysicsでも同じ式が使われています。が、透水係数の大きな場合にしか実用的ではないようです。1e-4㎝/sオーダーだと、γ/kが大きくなりすぎて、少しの速度差が大きな拘束力を生みます。土塊がすべろうとしても、地下水がその動きに追随するまで拘束されるような挙動も見受けられました。

２つ目は密度変化。

土の間隙内の水が地表に浸出した場合、水相の密度は大きく変化しますが、この表現が困難。いえ、pressureの式で参照している基準密度を変更すれば容易に反映できるのですが、粒子が重なり合い、挙動がやや不安定でした。初期配置も面倒で、実用的ではないですね。水収支はともかく、土の動きだけを正しく反映したい、ということであれば計算の軽い single-layer two-phase formulation でよいと思います。

土相側でも層毎に密度を変えている場合、その境界で密度が訛るのは現状では仕方がないと考えています。が、この密度変化の問題は将来的に何とか解決されそうな気がします。たとえば、土層別にフラグを付けて、個別に計算するとか。それが実現象を正しく反映するようになれば、より実務に浸透していくのでしょう。

今年前半より取り組んできた SPH コードの修正と実装は、これで一旦終了です。
結果的には単相での計算はOK、2相は「定性的評価」まで、土木分野の実務で求められる諸問題には適用が難しい場合が多いと判断されます。が、もう少しだと思います。

将来に期待しましょう。

double-layer two-phase formulation その2

振り返ると、double-layer での実装に取り掛かってから１か月以上かかっています。

以下を参考に始めたのですが、試行錯誤の結果、最終的にはPersianSPHの文献と同じような形に落ち着きました。
A coupled fluid-solid SPH approach to modelling flow through deformable porous media - ScienceDirect

SPH approach for simulating hydro-mechanical processes with large deformations and variable permeabilities | Acta Geotechnica (springer.com)

２次元、３次元共に土、水の相互作用までは動くのですが、地下水が地上に出た場合の密度変化をうまく表現できません。相互作用の検証も終わっていませんので、まだ時間がかかりそうです。

SPHでの２層2相表現は定性的にとどまると伺ったことがあります。それほど新しいものでもないので、国内でも実装されている方はいらっしゃいます。その先生がおっしゃるのですから、そうなのかもしれません。

まだ１か月。もう少し追いかけてみましょう。

SPHの水圧計算

前回から約1か月かかりましたが、２相２層での水圧計算に目途が立ちました。

https://www.sciencedirect.com/science/article/pii/S002076831730286X

上記の78-82式を実装することで解けるようにはなったのですが、まだ土層を出た後の水の密度変化までは安定して解けません。最初の含水率のまま解くことはできるのですが、空中でそれは誤りですし、勝手に密度を上昇させることもできません。圧力が下がって自然に密度が上昇する、というように計算させるとすぐに発散します。そのあたり、この論文では明記されていません。土層の計算だけ追うのであれば、１層での計算が良いと思いますが、まだ何か見落としているようです。

ひとまず、次の drag force を安定して解けるか見てみましょう。

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では密度変化がないか小さい現象（弱圧縮）を扱うのが前提なので、単体の異なる土層とか飽和‐不飽和‐地表流などを扱うのが難しいのでしょうね。

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

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

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

弾塑性SPH その2

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

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

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

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

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

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

壁粒子の圧力

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”.

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

弾塑性 SPH

2020年の秋から始めた弾塑性SPHコードの修正が完了しました。

基本コードを頂いてから寝かせていたため実質1年程度ですが、SPHを知ってからは14年目です。長い。
何とかSPHの離散化や各種補正方法、境界条件について基本的なしくみを理解したと言える証になるでしょう。いや、よかった。 

いくつかのハマり処があったのですが、式の実装時の静水圧Pと静水圧応力σの符号の整理が効いたと思われます。多くの文献を引用する際に、これらの書き方の流儀を統一してから実装に落とし込む必要があったのでしょう。

粒子が凝集する現象を hourglass,  zero energy mode だと勘違いしたり、他にも似たような勘違いで多くの寄り道をしましたが、振り返るとそれらも良い勉強になったと思います。
https://phreeqc.blogspot.com/2022/05/blog-post.html

早くVerification を済ませたいですね。ベンチマークは何が良いのでしょうか？
それが済んだら高速化。倍精度を多く使っていますので、GPUではなくスパコン利用（MPI並列）のほうが現実的でしょうね。

まだまだ先は長い。

SPH の補正

An SPH study on viscoplastic surges overriding mobile beds: The many regimes of entrainment - ScienceDirect

DualSPHysics を用いた土砂移動の新たなモデルかなと期待して読み始めましたが、非ニュートン流体実装を用いた2Dの計算のようで途中でやめてしまいました。モデル自体はこちらがわかりやすいでしょうか。210629.pptx (jsce.or.jp)

不安定化に対する補正を加えられていますね。SPHではこれらの補正がないと色々と不具合が出るので仕方ないのですが、多くの補正が加わると本当に正しい計算なのか疑問に思えてくることがあります（いえ、私が作ったコードのことです）。

ま、それでも補正して動くようにしたものを公開される研究者には頭が下がります。

DEM + 連続体

面白そうです。

Using FDEM-SPH model coupling to model a dynamic high-level landslide blocking a river in a deep valley area

Author links open overlay panel

https://doi.org/10.1016/j.enggeo.2023.107108

今年も福井県で天然ダムができていたようですね。
ここ数年、天然ダムの形成に関する数値解析を行ってきましたが、結論として DEM だけではダメ、流体（連続体）だけでもダメ、です。こういったカップリングが（おそらく）正しい答えを出す。そこから逆に簡素化してDEMだけ、もしくは連続体だけでできることを切り分けてその時々の評価対象を絞って検討するべきなのでしょう。

読んでみましょう。

地下水～地表流まで

・層流 : i=au
・層流から乱流への遷移状態、乱流 : i=au+bu²
・乱流 : i=bu²

1. 層流（粒径小、速度低）

Darcy：\begin{align*}-\mathrm{\nabla h}=\frac{1}{k}\boldsymbol{u}=a\boldsymbol{u}\end{align*}　h=p/(ρg)より\begin{align*}-\mathrm{\nabla p}=\rho g\frac{1}{k}\boldsymbol{u}=\rho ga\boldsymbol{u}=\alpha\boldsymbol{u}\end{align*}
Darcy-Kozeny Carman:
\begin{align*}a=\frac{1}{k}=\frac{180\nu\left(1-n_e\right)^2}{g{\mathrm{\Phi}^2n}_e^3d^2}\end{align*}　ν=μ/ρより\begin{align*}\alpha=\frac{\rho g}{k}=\frac{180\mu\left(1-n_e\right)^2}{{\mathrm{\Phi}^2n}_e^3d^2}\end{align*}

2. 層流～乱流への遷移状態、乱流（粒径大、高速）

Forchheimer：
\begin{align*}-\mathrm{\nabla h}=a\boldsymbol{u}+b\left|\boldsymbol{u}\right|\boldsymbol{u}\end{align*}
Ergun(1952)
\begin{align*}a=\frac{1}{k}=\frac{150v\left(1-n_e\right)^2}{gn_e^3d^2}
　b=\frac{1.75\left(1-n_e\right)}{gn_e^3d}\end{align*}

Kadlec and Knight(1996):角ばった粒子に適用。
\begin{align*}a=\frac{1}{k}=\frac{255v\left(1-n_e\right)}{gn_e^{3.7}d^2}
　b=\frac{2\left(1-n_e\right)}{gn_e^3d}\end{align*}

3.乱流

Manning：
\begin{align*}u=\frac{1}{n}R^\frac{2}{3}I^\frac{1}{2}\end{align*}\begin{align*}I=\frac{n^2}{R^\frac{4}{3}}u^2=bu^2\end{align*}
h：水頭（m）
p：ポテンシャル（Pa）
ρ：密度（kg/m3）
k：透水係数（m/s)
u：速度（m/s）
α：係数（s/m）
β：係数（s2/m2）
Φ：球形度(0～1、球=1)
d：有効粒径（m）
μ：粘性係数（Pa・s）
ν：動粘性係数（m2/s）=μ/ρ
ne：間隙率
g：重力加速度(m/s2)
n：マニングの粗度係数(m^−1/3・s)
R：径深(m)
I：勾配

OpenFOAM
https://openfoamwiki.net/index.php/DarcyForchheimer
https://phreeqc.blogspot.com/2022/02/darcy-forchheimer-model.html?m=0
\begin{align*}-\mathrm{\nabla h}=a\boldsymbol{u}+b\left|\boldsymbol{u}\right|\boldsymbol{u}\end{align*}\begin{align*}-\mathrm{\nabla p}=\rho ga\boldsymbol{u}+\rho gb\left|\boldsymbol{u}\right|\boldsymbol{u}=\alpha\boldsymbol{u}+\beta\left|\boldsymbol{u}\right|\boldsymbol{u}=\mu D\boldsymbol{u}+\frac{1}{2}\rho F\left|\boldsymbol{u}\right|\boldsymbol{u}
\end{align*}
Darcy:
\begin{align*}\alpha=\frac{\rho g}{k}=\mu D\end{align*}\begin{align*}D=\frac{\rho g}{k\mu}\end{align*}\begin{align*}F＝0\end{align*}
Darcy-Kozeny Carman:
\begin{align*}\alpha=\frac{\rho g}{k}=\frac{180\mu\left(1-n_e\right)^2}{{\mathrm{\Phi}^2n}_e^3d^2}=\mu D\end{align*}\begin{align*}D=\frac{180\left(1-n_e\right)^2}{{\mathrm{\Phi}^2n}_e^3d^2}\end{align*}\begin{align*}F＝0\end{align*}
Forchheimer-Ergun:
\begin{align*}\alpha=\rho ga=\frac{\rho g}{k}=\frac{150\mu\left(1-n_e\right)^2}{n_e^3d^2}=\mu D\end{align*}\begin{align*}D=\frac{150\left(1-n_e\right)^2}{n_e^3d^2}\end{align*}\begin{align*}\beta=\rho gb=\frac{1.75\rho\left(1-n_e\right)}{n_e^3d}=\frac{1}{2}\rho F\end{align*}\begin{align*}F=\frac{3.5\left(1-n_e\right)}{n_e^3d}\end{align*}

PersianSPHでは、式と係数を選択。ソース内で係数変更可。
OpenFOAMでは、D,F 入力で表現。

PersianSPH その3

最後に、土‐水連成を作ってみましょう。
と古いサンプルに手を付けましたが、使えないコマンドがありました。

・dom.SeepageType
・dom.Time
・dom.KernelType
・dom.VisEq

Version 違いでしょうね。
Source を見てみると、SeepageType はParticle に対して指定するようです。無駄なので、後で修正されたのでしょう。Time は Private から出してやるとコンパイルが通りました。良いのかな？他に影響が出たらその時考えましょう。
後半２つは新しい？サンプルの書き方に修正することで通りました。

 土‐水連成のキーとなる変数は2つのようです。

・SWI
・SeepageType 

SWI:

0 => The seepage force + The bouyant unit weight of soil
1 => The seepage force + The surface erosion(Lift+Drag) + The bouyant unit weight of soil
2 => The seepage force + The pore water pressure from water particles
3 => Zero interaction force

1の文献が見当たりません。ソースは以下の通り。河床変動のようにせん断力を考慮しているのだろうと想像するのですが、何から引っ張ってきた式かは追えませんでした。

double Cd = 24.0*(P2->MuRef/P2->RefDensity)/(P1->d*norm(v)+0.01*h*h) + 2.0;
SFt = (3.0/(4.0*P1->d)*P2->RefDensity*(1.0-P1->n0)*Cd*norm(v)*v) *K;
SFt(1) += (P2->RefDensity*(1.0-P1->n0)*norm(v)*fabs(P2->S-P1->S)) *K;

SeepageType :

0 => Darcy's Law
1 => Darcy's Law & Kozeny–Carman Eq
2 => The Forchheimer Eq & Ergun Coeffs
3 => The Forchheimer Eq & Den Adel Coeffs

Kozeny–Carman だけかと思っていましたが、他にもありますね。https://phreeqc.blogspot.com/2021/01/kozeny-carman-equation.html
deとd15の片方しか指定する箇所がなさそうなので、こちらの文献に従って組まれたのだと思われます。
https://www.sciencedirect.com/science/article/abs/pii/S0266352X17302318
係数α, βはソース内指定です。透水係数に応じてソースを変更する感じでしょうか。

土‐水連成の場合、強度、透水性、SWIなど調整すべき項目が弾塑性に比べて増えます。計算時間は単純に倍。２次元でもパラスタが容易とは言えません。３次元ならなおさら。
ということで試算は２次元を選択。堤防の越流破堤を想定しモデルを作ってみましたが、SWI＝2の The surface erosion(Lift+Drag) が効きすぎ。ここを河床変動式に変更すれば、もっと現実っぽくなるのでしょう。
最終的には流入水の勢いで極端にえぐれるか、越流しても壊れないかの設定になり打ち切ってしまいましたが、現実的な実験値と比較すればそれなりの形状は得られるかもしれません。探してみましょうか。

以上で PersianSPH の試算は終了。弾塑性と浸透＋越流の確認は十分とは言えませんので、今後、もう少し試してみましょう。

PersianSPH その2

次は弾塑性。 

見たところ、スタンダードな SPH による解き方を実装しているようです。弾塑性のフレームワークも標準的。

弾塑性のサンプルがなかったっため、DamBreak の水を土に変更してみました。

2次元では問題なく動きます。

3次元ではバクハツ。
これは SPH のデフォなのでしょうか。ただし、zero-energy mode は出ていません。境界部分の過剰な正圧に起因するだけのようです。優秀ですね。
particle penetration を防ぐための Penalty Parameter が効きすぎていると踏んで、1→0.5に下げてみました。
バクハツは収まりましたが、これはダメ。底面で粒子の沈み込みが発生しています。
底面：0.95、側面：0.7程度が良さそうです。

粘着力を変えると、塊のまま移動する土塊とすべり層が表現されて良い感じ。
拡大すると、

せん断変形が孔内傾斜計の累積変位図のようにきれいに表現されています。引っ張りも安定して解けています。

問題は計算時間。
1秒の計算に8コアで4時間。時間がかかります。OpenMPなので、スパコン利用も不可。実務に持っていくには gpu に載せるか、MPI並列にするか。

続く。