メールファイルの移動、HPのボタン押しなどを Power Automate で実施。
使ってみて気づいたのは、以下の点。
後輩君に聞いてみたら、同じところ（自動実行）で困ったとのこと。どうして解決したのか聞いたら、VBAと組み合わせたらしい。微妙です。
・層流 : i=au
・層流から乱流への遷移状態、乱流 : i=au+bu2
・乱流 : i=bu2
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 入力で表現。
最後に、土‐水連成を作ってみましょう。
と古いサンプルに手を付けましたが、使えないコマンドがありました。
・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 の試算は終了。弾塑性と浸透＋越流の確認は十分とは言えませんので、今後、もう少し試してみましょう。
次は弾塑性。
見たところ、スタンダードな 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並列にするか。
続く。
