メールファイルの移動、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並列にするか。
続く。
SPH でも弾塑性を扱いたい、ということで PersianSPH の登場です。
まずは、インストール。HPに従って進めます。
https://www.m2clab.com/persiansph/installation
$ wget "https://github.com/mghkorzani/persiansph/blob/master/persiansph_install.sh"
$ bash persiansph_install.sh
ここで止まります。2年前に更新が止まっているようですので、仕方なし。
GitHub のリポジトリにあるコードを zip で落として展開。直下に persiansph_install.sh があるので、これを使いましょう。
$ bash persiansph_install.sh
Installation is completed.
Close all terminal windows and open a new one to take effect the defined environment variables
Please refer to the tutorial to run a simulation.
persiansph フォルダが作成されます。その中の、CmakeLists.txt を書き換えて目的の cpp ファイルを指定（サンプルから一つを選択）。
その後、新しいフォルダを作り、そこからcmake。
$ mkdir (New Path)
$ cd (New Path)
$ cmake $SPH
$ make
実行ファイルが（New Path）に作成されましたので、これを動かすと計算開始。フォルダ内に output ファイルが書き込まれます。ファイル形式は hdf5 と xmf。表示は Visit or ParaView。hdf5 は読めませんでしたが、xmf は両者とも OK でした（ParaView では XDMF Reader or Xdmf3Reader）。2D だとVisitのほうが見やすいかな。
続きは後日。
**********************
20220603修正
ParaView では XDMF Reader が正解でした。Xdmf3フィルターをかけた場合に正常に表示されない場合が出てきました。
hdf5 と xmf はセットでした。xmfを参照し、hdf5からデータを取り出しているようです。
ダイヤさんの Dtransu の Version が 2.0になっていました。
大きな変更点は、並列化対応でしょうか。
OpenMP ですので「ようやく」感は否めません。手元にあった ver.1 改良ソースを見ると、OpenMP で並列化したのが11年前。その他いくつか改良を独自に加えていましたが、それらの8割ほどは今回の Ver.UP で追いつかれました。
Post 処理でも ParaView に対応です。これは個人的に嬉しい。
あとは pre だけですね。2D版 は Gmsh 用の変換プログラムが作られていましたが、3D版は「商用ソフトで」となっていました。残念。
１０年以上、改良成果を無償で公表される取り組みは称賛に値します。次の技術者も育っているのでしょうね。
今日、国会図書館に文献の複写を依頼しました。
まだ、紙での提供のようでした。ペーパレスの時代にメールで送ってもらえないのは、資源も時間も無駄です。
著作権が引っかかっているようなのですが、問題になる箇所を詳しく理解していません。
https://crd.ndl.go.jp/reference/modules/d3ndlcrdentry/index.php?page=ref_view&id=1000177352
海外の雑誌は PDF で購入できます。日本でも PDF に対応したサービスはあるようですが、コピーに配慮が必要なのは紙の時代と同じです。
https://jdream3.com/question/copy-etc.html
国会図書館をはじめ国内で PDF 配布されるようになるのは、まだまだかかるのでしょうか。
