書店で面白そうな図書を見かけました。
中川(監訳)「環境保全のための地下水水質化学-地球化学，地下水および汚染-（上）·（下）」
国内の図書にしては珍しく PHREEQC の input data が載っている!と思いきや、ほぼ記憶にある内容。あれ？と思って表紙を見てみると、C.A.J. アペロ「Geochemistry, Groundwater and Pollution, 2nd edition」の翻訳書でした。日本語版が出版されたのですね。昨年発売ですので、遅れること16年。国内にも細々と使われている方がいらっしゃったようです。
この分野、国内には良い図書や文献がほとんどない状況。理解したいなら海外の知見に頼らざるを得ません。そうすると、最初から英語で学んだ方が良いように感じます。が、この洋書1冊を理解するのにかかった時間を考えると、最初は訳本を併用した方が効率的かもしれません。
ぜひ、hydrogeochemistry の裾野を広げる契機となってほしいものです。
OpenFoam で河川の計算に再トライ。
その中で、地下水を連成する際に必要な水位固定の境界条件を探していました。が、どうも、デフォでは実装されていないようです。古い情報ですが、groovyBC を使えば inlet 適用できそうでした。が、側方に張る場合はどうすべきでしょう？
https://www.cfd-online.com/Forums/openfoam-pre-processing/122819-inlet-boundary-condition-based-specified-water-level.html
2週間ほど調べましたが、良い設定が見つからず。
ひとまず、側方は不透水にしてin-out のみの境界条件を張って計算してみると、川と地下水の流れは定性的に計算できそうでした(中流で池のように溜まっていましたが)。山岳地の河川で上流端と下流端の流量の違いを表現したいような場合には使えるでしょう。
はい、眠らせましょう。
次は、河川・地下水間の移流拡散。これはinterFOAM をベースに作ろうかと思ったのですが、既に実装されていたようでした。
interFoam with scalar transport in single phase
https://www.cfd-online.com/Forums/openfoam-solving/234344-interfoam-scalar-transport-single-phase.html
相間のスカラー輸送を計算する場合と、単一相内での輸送は切り分けないといけないですね。当然なのですが、これを読むまで気づきませんでした。
早速、リンク先の tutorial を動かそうとしたのですが、コンパイルが必要でした。手元にあるのは native Win ＋ ESI 版2012。Ubuntu の入った laptop が故障中のため、新しく環境をつくらないといけない？面倒だなと思いつつ一旦休憩へ。
シャワーを浴びながら、この方法では問題があることに気づきました。流速の問題です。OpenFOAM では Darcy-Forchheimer Equation を用い、porous media に対して抵抗としての圧力（勾配）や方向を補正し流速を導く計算でした。間隙比を陽に扱っていないため、実流速を計算できせん。汚染の計算はこれではダメ。地表流のみならOKですが、地下水部分の汚染の拡大が遅くなります。そもそも、地下水の量もダメでしょうね。目の詰まった遅い流れは対象外なのでしょう。
そうなると、汚染の計算は河川のみ。それなら iRIC を選択すべきとなります。はい、眠らせましょう。
見過ごした情報もあるかもしれませんが、3週間でたどり着いたのが「使えなさそう」。なかなか良いツールはありませんね。
地質の3次元化の問い合わせが寄せられるようになりました。
設計に使いたいと仰る積極的な方や、提案書に書いてしまったので仕方なくという消極的な方がいらっしゃいます。お客様の中には、明日までに事例が欲しいと仰る方も。皆さん、タイムリミットが迫って腰を上げ始めたというのは共通のようです。
現時点で、地質のモデリングは Civil3D＋GEORAMA を選択せざるを得ないのでは、と考えていたのですが、設計者の視点では V-nas Clair ＋ GEO_Kit で作成した支持層サーフェスで充分、とおっしゃる方もいらっしゃいます。後者のメッシュは粗く、ユーザー側で制御できないため意図した面ができないのですが、意外と受け入れられているようです。
面の推定方法は GEORAMA と GEO-Kit で同様のアルゴリズムが使われています。
面の優先順位指定方法は GEORAMA に実装されていますが、GEO-Kit にはありません。
ボーリング読み込みは両者とも xml 読み込みに対応しています。
属性は Civil3D → Navis で概ね作成できますが、IFCに出力できません。Civil3D 上で付与可能なソフトも販売されているようですので、今後はそちらが主流になるのかな。
V-nas Clair では i-ConCIM_Kit で 属性付与可能です。といっても、外部参照の属性ファイルはイチから作ることになります。外部参照しなくても、IFC 書き出し時に柱状図の属性が少しだけ付与されます。
地質のモデリングが必要であれば GEORAMA、ボーリングモデルの納品のみであれば GEO_Kit でしょうか。設計者の意向も確認したいですね。
隣の地質を担当している部署では今でも2次元の世界から出てこないオジサマが多く、CIM に迫られても対応しようと手を動かしている方はいません。CIM を振られたら、他支店に投げるのでしょう。一番のネックはソフトでなくココかもしれません。
GEORAMA v3.3 から吐き出した Nastran のメッシュを、FLAC に渡そうとして思案。
FLAC 7.0 で対応している外部メッシュは Ansys の lis とAbaqus の inp。どちらも手元にありません。
何かないか？と探してみると、ありました。
meshio
https://github.com/nschloe/meshio
conda install して使ってみましたが、nas の読み込み時にエラーが発生します。nas の Ver. があってないのかな。FLAC の f3grid に対応していましたし、ParaView のプラグインもあったのに惜しい！
次。
こちらはOpenFOAM で利用していました。フリーのメッシャーとしか認識していませんでしたが、色々なフォーマットの入出力に対応していました。メッシュ変換ソフトとしても利用できそうです。
なぜか1台のPCでは export できなかったのですが、別のPCでは問題なく変換できました。感謝です。
どこかの商用メッシャーよりも優れています。使用しているメッシャーの種類を聞けば、その方々の技術力もある程度ふるい分けできそうな予感。
地層研さんの Geo-Graphia Ver.2.3 は、どこまで CIM に対応しているかを調べてみました。
別部署の先輩が購入したのですが、その部下がモデルの画像キャプチャーに使う程度で、ほとんど稼働していないようでした。
個人的には、地層研さんの地下水プリポストは問題が多く、サポートの対応も良くなかったので、現在使用していません（ブログ内をG-TRANで検索して、詳細を思い出しました）。が、Geo-Graphia は触ったことがなく、空いているなら CIM 機能だけでも把握しておこうと思い触ってみることにしました。
地層研さんのこだわりなのでしょうか、undo はいまだに実装されていません。操作性もイマイチ。後輩さん曰く、拡大すると落ちる現象は、引き継がれているようです。
そのあたりの作りは大きく変わらなかったのですが、export できる内容は増えていました。が、これもマニュアルを読むだけでは具体内容を理解できず。実際に試してみました。
export:
サーフェス（地形）
dwg:連結されたポリメッシュで出力
→これはCivil3DでTINサーフェス化可能ですが、地形はどのソフトでも作成できるので出番はないか。
地層モデル（地質）
dwg:地層境界面が連結されていないポリメッシュor3D面で出力
inp:ファイルは作成されるがParaViewで読めず。
stl:ファイルが作成されない。
地層モデル（ボーリング）
dwg：ファイルは作成されるが何も書き出されず。
stl:ファイルが作成されない。
ソリッドモデル
↑サーフェスモデラーのようで、ソリッドと言いながらもstl（サーフェス）を吐き出します。
dwg:ファイルは作成されるが何も書き出されず。
inp:ファイルは作成されるが何も書き出されず。
stl:ジオメトリをパーツ単位でEXPORT。
→相対位置関係を保ったまま blender でFBXへ変換可能。
本来はマニュアルに記載しておくべき内容ですが、多くの制限がありそうです。今回の試行で使えそうだったのは、、ナシかな。一般的なソリッドではなく、面なので。妥協して外側の形を見せるだけとしたら、stl でしょうか。
メッシング機能も以前と大きく変わらず、SoilPlus や FEA NX に比べると機能不足。
残念ながら、CIM にも対応しているとは言えません。柱状図xmlを読めませんし、N値も表示できません。ソリッドモデルをIFCに書き出すことも、属性を与えることもできません。あくまで3次元可視化ソフトの位置付けなのでしょう。
やはり私は使用することはないでしょうし、お勧めもしません。が、それを再確認出来た点はヨカッタ。
先輩さんの部署では元を取るまで使っていかざるを得ないでしょうね。オプションも含めると結構な金額でしたので。騙されたと思われるか、自分の目利きを後悔されるかはわかりませんが、いずれにしてもお気づきになるまで相当の時間がかかるでしょう。
次は Nays2D+ Ver1.0.8。
こちらは iRIC Ver.3 から登場とのこと。平面２次元の計算にて疑似的な３次元（鉛直方向）の流速まで表現できるそうです。が、興味があったのは濃度の移流拡散機能。河川の流れに物質を乗せて運んでくれます。
拡散係数は渦動粘性係数に乗じる値となっています。Nays2D＋に説明はありませんが、UTT で以下の解説があります。
K=aν+b （ν：渦動粘性係数）
https://i-ric.org/yasu/UTT_JP/02_Overview.html
ガウス分布ベクトル L に乗じて適度に散らす効果をあたえているようですが、この方法は知りませんでした。河川分野ではよく使われているのでしょうか？負荷の小さい計算で濃度をあわせ込めるのであれば、お手軽です。
で、計算。
適当な河川を選んで 5mLP をインポート。今回は計算条件を設定して不等流計算を行い、定常になってから物質を流します。流し始める位置は土石流同様にポリゴンで指定。簡単です。
最上流から500秒間、1%濃度の溶質を流した結果です。
背景は国土地理院の色別標高図を表示しています。
https://maps.gsi.go.jp/development/ichiran.html#relief
お手軽ですね。
調べてみると、TELEMAC の ポスト処理としてＢlueKenue にも似たような機能が搭載されていました（こちらは sediment transport)。案外、広く使われているのかもしれません。
現段階では境界条件として時系列の濃度変化までは設定できないようですし、吸着も扱えません。が、そのうち実装されるでしょう。コチラも気長に待ちましょう。
iRIC Ver.3 に含まれている土石流解析ソルバーの Morpho2DH を使用。Ver.1.0です。
ソルバーマニュアルには、「Morpho2DH は平面二次元の土石流・泥流モデルを主体とした解析ソルバであり，斜面崩壊を初期条件とした土石流・泥流の流動・堆積過程を表現可能なモデルです」とあります。
まず、tutorial。前半は問題なく進みましたので、動くようです。
プレゼンで示されていたように、非常にお手軽に計算ができます。tutorial では xyz の csv データが用意されていましたが、iRIC のプレは基盤地図情報からのダウンロード機能を備えています。座標系を指定しておけば背景に地理院地図等を表示できますので、それを見ながら DL する範囲の指定が可能。範囲内に ５mLP がなければ、ある部分だけで点群を作ってくれます。この一連機能は秀逸。LandXML への吐き出しも実装されているようですので、土石流や河床変動以外の解析でも、この部分だけ利用する場面が出てきそうです。
次は、後半の固定床を飛ばして、実例をトレースすることに。
5mLP をインポートし、流下方向を考慮しながらメッシュを切る。空中写真を背景に頭部崩壊領域を作成し、メッシュに属性を付加、計算条件を入力したら設定完了です。
今のところ GUI は、こけていません。たまに属性のマッピングが反映されていないケースに遭遇しましたが、どういうタイミングで発生するのか掴めていません。念のため保存して終了し、再起動後に読み込んで入力値が入っていたらOKです。同じ類の Hyper-Kanako に比べ、操作性・安定性はマシマシです。
で、計算！
はい、結果が出ません。すぐに発散しています。ソルバーの中にタイムステップを自動調整してくれる機能がないようです。しかも並列化されていないので、計算が遅い。
タイムステップを 0.01 秒から 0.001 に切り替えてもダメ。
メッシュを 10m から 5m に切りなおすと、動きました。
ひとまず一つの流路では動いたので、同時多発を試行。地理院の空中写真を背景に切り替えて10か所超の頭部崩壊領域を設けました。発生時間をずらすように設定し、計算スタート!
結果はそこそこ。
建物の影響を考慮していない地形ですので実際の堆積域とは一致しません。が、計算は走りました。力技ですが、メッシュを細かく切り「障害物」として建物を指定すると、その影響は取り込めるようです（さらに計算は遅くなりますが）。
iRIC のポスト処理として、地理院の空中写真を背景に層厚や速度を重ねることができます。動画も書き出せます。操作が簡単で安定していることを含め、フリーのプリポストとしては作りこまれている感じです。
Morpho2DH としては、同時多発の土石流を一度に計算できることが特徴でしょう。改善要望としては、並列化と建物の崩壊を含めた影響検討ですかね。いつか実装されるかもしれませんので、気長に待ちましょう。
一週間ほど、OpenFOAM の interFOAM (VOF 法)で、河川の自由表面流れを計算できないか探っていました。
結論としては達成できず。力不足です。
tutorial がなく、計算例も web 上に見当たらず。どうも境界条件の与え方が良くないのだろうとあたりを付けたものの、どのように修正して試行すればよいのかわかりません。こうなるとサポートのない Open Source では、初心者は行き詰ってしまいます。
で、寝かすことに。
地下水と連成しないのであれば iRIC でよいか、とそちらへ移行。
これ、実務では使ったことがありません。これまで実務で扱った土石流や河床変動計算は、現場によって独自の条件が加わっており、パッケージソフトでは対応できなかったのが実情でした（河川のプロはメッシュ作成に使っていると聞いています）。
以前、iRIC では簡単な土石流計算が手軽にできるといったプレゼンを見せていただいてから、いつか私も使ってみようと考えていました。
で、購入した図書がコチラ↓
木村「iRICによる河川シミュレーション」
ソフトは Ver.3.0.19 を DL。
早速、図書の順番通りに 1D から動かしてみようと思いきや、Ver.3 には 1D が入っていない模様（別途入れるのでしょうか？）。
2D で 1D 計算をすればよいかと動かしてみたら、今度は fortran の error。
作成したファイルをもう一台のPCで動かすと、計算してくれました。原因がわかりませんが、気にしないことに。
2章の内容を Nays2DH を使用してトレースしただけですが、概ね手順を理解できました。なかなかお手軽です。ポスト処理が楽ですね。
そのまま4章まで読みましたが、難しいことは書かれていません。CIP：空間3次、風上差分：１次を実装しているようですが、前者で振動するなら後者で、だそうです。
実際の地形を考慮するのは５章以降。引き続き読んでみましょう。
*******************************
20220213
error はバグとのこと。以下で動きました。
https://i-ric.org/forum/tmp%e3%83%95%e3%82%a1%e3%82%a4%e3%83%ab%e3%81%ab%e3%81%a4%e3%81%84%e3%81%a6/1D
1Dソルバーは 手作業で追加でした。
iRICインストールフォルダ内の「solvers」へ。
https://i-ric.org/download/solver-install-manual/
OpenFOAM のDarcy-Forchheimer model で設定する D,F を整理。ついでに透水係数との関係も整理。
まずは、公式wikiより。なぜか太字が出ないのは御愛嬌。
https://openfoamwiki.net/index.php/DarcyForchheimer
The Darcy Forchheimer acts in the momentum equation as a sink term Sm. Considering the momentum equation, it follows:
\begin{align*}\frac{\partial\rho U}{\partial t}+\mathrm{\nabla}\left(\rho UU\right)=\mathrm{\nabla\sigma}+S_m
\end{align*}
The main important term is the source term Sm which is given as:
\begin{align*}S_m=-\left(\mu D+\frac{1}{2}\rho tr\left(U\cdot I\right)F\right)U\end{align*}
As the source term Sm can be expressed as a pressure gradient it follows Sm=∇p. Therefore, we can write:
\begin{align*}\nabla p=-\left(\mu D+\frac{1}{2}\rho tr\left(U\cdot I\right)F\right)U\end{align*}
While switching to the Cartesian coordinate system, we can achieve something similar to:
\begin{align*}\nabla p=-\left(\mu D_i-\frac{1}{2}\rho F_i\left|u_{kk}\right|\right)u_i\end{align*}
If we know the velocity dependent pressure we can calculate a polynomial function such as:
\begin{align*}-\mathrm{\nabla p}=Au+Bu^2\end{align*}
\begin{align*}A=\mu D\end{align*}
\begin{align*}B=\frac{1}{2}\rho F\end{align*}
ｐ：圧力損失P[N/m2], u：流速［m/s］, μ：動粘性係数[N/m2・s]=[kg*m/s2/m2*s]=[kg/m/s], ρ：流体密度[kg/m3]
CfdOF では、D、F を constant/fvOptions で指定しています。
F＝0 でダルシー則。単位は以下。
∇P [N/m3] ＝Au[N/m4*s*m/s]＋Bu２[N/m5*s2*m/s*m/s]
係数A[N/m4*s]=[kg*m/s2/m4*s]=[kg/m3/s]
係数B[N/m5*s2]=[ kg*m/s2/m5*s2]=[kg/m4]
D=A/μ[ N/m4*s*/N*m2/s]=[1/m2]
F=2B/ρ [kg/m4/kg*m3] =[1/m]
吉岡ほか（2010）「高透水性多孔質体中の非ダルシー流れに関する考察」では、ρg [kg/m3*m/s2=N/m3] で除し無次元化した式が紹介されています。係数 a、b に関して Ergun 等の近似式を検証しています。こちらの式の方が土木では‘なじみ’があります。
\begin{align*}-\mathrm{\nabla h}=au+bu^2\end{align*}
\begin{align*}a=\frac{150v\left(1-\emptyset\right)^2}{g\emptyset^3d^2}\end{align*}
\begin{align*}b=\frac{1.75\left(1-\emptyset\right)}{g\emptyset^3d}\end{align*}
h：水頭［m］, a：係数［s/m］, b：係数［s２/m２］, ｄ：粒径[m], φ：間隙率［－］, ν：動粘性係数［m２/s］, g：重力加速度［m/s２］
u=-1/a∇hより、透水係数 k＝1/a[m/s]
換算は以下の通り。
\begin{align*}D=A/\mu＝aρg/μ=\frac{150v\left(1-\emptyset\right)^2}{g\emptyset^3d^2}ρg/μ\end{align*} [s/m*N/m3/N*m2/s]=[ 1/m/m3*m2]=[1/m2]
\begin{align*}F=2B/\rho=2b\rho g/\rho=2\frac{1.75\left(1-\emptyset\right)}{g\emptyset^3d}\rho g/\rho\end{align*} [1/m]
あるいは、
\begin{align*}D=\frac{\rho g}{k\mu}\end{align*}
***********************
20220611追記
見通しを良くしてみました。
https://phreeqc.blogspot.com/2022/06/blog-post_11.html
集水地形に雨を降らせて、地下水と表流水の変化を見てみることにしました。
といっても、まだ実務レベルではありません。先日の CfdOF (OpenFOAM) で、どこまで手間をかけずに計算できるか？という確認レベルです。
まずは、モデル作成。
あとは実行。
ですが、途中で止まります。やや乱暴な境界条件でしたので「それが原因かな？」と思いつつ、それを修正するにはジオメトリ修正からなので後回し。簡単に修正可能な設定から見直すことにしました。
「OpenFOAMによる熱移動と流れの数値解析 第2版」と著者のサイトを参考に、3日ほどあれこれ試していたら、走るようになりました。効果的だったのは、以下の3点。
ひとまずこれでOK。出てきた結果がコチラ↓
うーん、コレジャナイ。いえ、最初だから砂利山のような高透水を設定したので、これで正解なのですが。
透水性を下げてみると、
それらしくなりました。
ひとまず、エラーは回避できました。
