Turbidity current in Monterey Canyon

昨日の文献で示されていた、avalanche/tutorials/montereycanyon の Allrun を読んでみました。

Rauter, M. et al.(2018) Applications of the Finite Area Method on a Geographic Scale: From Dense Snow Avalanches to Turbidity Currents

The second model is based on the work of Parker et al. (1986) for turbidity currents. The solver models a turbulent mixture of sediment and water, including entrainment of sediment at the bottom and water at the top. It will be applied to a turbidity current in Monterey Canyon, off the coast of California.

参照先のソースを継ぎ接ぎし、要点のみ一読できる形にまとめ、見通し良くしました。
まず、この例では海底地形の STL ではなく、Esri ASCII ラスター形式読み込んでいます。ArcGIS で xyzポイント→ラスター変換→ASC保存 をお考えなのでしょう。文献では "GIS data is translated to OpenFOAM dictionaries using python scripts. "とありましたが、その script は見当たりませんでした。どこかにあるのでしょう。
解析範囲はshpを読み込んでいます。

gridToSTL
    system/gridToSTLDict
        //Output file
        stlname "constant/surface.stl";
        //Topography file
        gridname "constant/gisdata/bathymetry.asc";
            ncols        1194
            nrows        620
            xllcorner    574999.06
            yllcorner    457937.36
            dx           24.71
            dy           31.27
            NODATA_value  0
             -166.5 -164.13 -158.49 -154.69...
        .//Choose between boundary: fromShape or fromPoints:
        //Read Boundary from Shapefile
        boundary fromShape;
        shapeBoundary "constant/gisdata/aoi_UTM_N10";
        //Division of an edge of the boundary polygon
        divisions 200;
        //Height of the Domain
        domainHeight 1000.0;
        //Offset to reduce coordinate size
        offset (-574000.0 -4057000 0);

次に、polyhedral mesh と finite area mesh の作成。前者で meshDict を読んでメッシュを作成していますが、次の faMesh との連携がわかりません。文献ではcfMeshに関し"A mesh is generated that covers the volume above the terrain of interest; but only the bottom boundary mesh is utilized."とあり、faMesh が底のメッシュのみを取り出して2Dに変換する役割を担っているのかなあと想像しています。
release area は inlet で処理しているとのことでしたが、その詳細も見つけられませんでした。

pMesh: cfMesh-polyhedral mesh
    system/meshDict
        maxCellSize 200;
        surfaceFile "constant/surface.stl";
                additionalRefinementLevels 2;
                "wall0"
                    newName "inlet";
                "wall(9|10|11)"
                    newName "outlets";
                "top"
                    newName "top";
                "terrain"
                    newName "terrain";
makeFaMesh: finite area mesh
    constant/faMesh/faMeshDefinition

あとは通常運転。
０フォルダの作成（複製)。
restore0Dir: restore initial fields

ソルバーの設定、計算。
$(getApplication)
    system/controlDict
        faParkerFukushimaFoam
        #include "shapefileWrite":
            offset (-574000.0 -4057000 0);
        #include "gridfileWrite"
            ncols, nrows...offset
        #include "autoAreaToVolumeMapping"
            prefix  "fa_";

いくつか不明点は残りましたが、概ね見えてきました。

Gravitational mass flows + OpenFOAM

Rauter, M. et al.(2018) Applications of the Finite Area Method on a Geographic Scale: From Dense Snow Avalanches to Turbidity Currents
https://www.itn-slate.eu/wp-content/uploads/2019/10/ESR_14_MatthiasRauter-Abstract-7thESIOpenFOAM-Berlin2019.pdf

河川・砂防分野でよく使われる２次元浅水方程式にスカラー輸送、沈殿等を加味し OpenFOAM2006 に導入されています。$FOAM_TUTORIALS/incompressible/shallowWaterFoam/ のような簡単な例をベースに、雪崩や土石流などの Gravitational mass flows (GMF) に特化した形で実装されたのでしょう。

手元のWin11＋ OpenFOAM 2012 だと（当然ですが）コンパイルできませんでした。ダメもとでUbuntu20.04 ＋OpenFOAM９でも実施。こちらも Foundation 版には含まれていないソルバーを見に行ったみたいでダメでした。

正攻法の Ubuntu20.04 ＋ OpenFOAM2112 だとコンパイルが通りました。で、tutorial の Avalanche を Allrun。完走しました。
paraFoam で結果を見るとそこそこ。GIS上に持ち込めるようですが、中を読まないとわからないですね。説明が省略されていますので、地形の導入から結果の export まで、中を見たらわかるということなのでしょう。単純な GMF なら、手元の Fortran プログラムや iRIC の方が良いのかもしれません。

OpenFOAM ユーザーが多いのでしょう。多くのソルバーや解析例が公開されています。ありがたいのですが、初心者には敷居の高い open source です。

ObsPy その3

Obspy 1.2.2 で　import　できなくなりました。numpy で引っかかっています。

仮想環境を再構築してもダメ。
1.2.1 に douwngrade してもダメ。

Git で issue を検索するとヒット。
Import error for obspy v1.2.2 (windows, pip install) #2940
https://github.com/obspy/obspy/issues/2940

$ conda search numpy
$ conda install numpy=1.21.5

numpy を落とすと動きました。
不用意に upgrade してしまい、1時間ほど無駄にしました。

SciKit GStat

解析値と点群から作った地形グリッドの差分を可視化し、チェックする必要が出てきました。

数年前は Surfer で対応していましたし、Civil3Dでもよく実施していた内容です。
ただ、解析値が出るたびに手作業が発生するのは非効率。これ以外の計算や図化はプログラムに載せていましたので、この部分もその流れに含めたいところです。

さすがにクリギングぐらいは Python で容易にできる時代になっているでしょう、と思って検索したら、ありました。

SciKit GStat
https://scikit-gstat.readthedocs.io/en/latest/

地球統計学ライブラリは複数ありましたが、今回はこれを選択。必要な機能を最初に確認できたライブラリで、BESTということではありません。
解析値に対し Gaussian model で variogram を計算し、クリギングによる補間、グリッド値を算出します。これだけです。

あとは他のグリッド状の解析値と合成し、LPグリッドとの差分をとり、図化・出力。これで手作業から脱却できました。便利な時代です。ありがたい。

FloPy その2

自身も FloPy を動かしていなかったことに気づきました。

ひとまず tutorialと example をいくつか動かしてみました。
MODFLOW 関連ソフトは含まれているのかと思いきや、別途DLが必要でした（当然）。MT3DMS はUSGS版でない方が必要でした。

動かしてみて気づいたのは、MT3DMS にいくつかの吸着式や離散化スキームが実装されている点。昔、Visual-MODFLOWを触っていた時にも選択できたのでしょうが、気づきませんでした。
結果を見ると、風上差分が訛っています。これ、OpenFOAM の1次精度風上差分でも同様でした。2次精度では改善していましたので、MT3DMS の upstream は前者なのでしょう。MOC（とHMOC）だと訛らない、MMOCだと訛るというのは意外。やってみないとわかりません。

先の揚水量と透水係数の変化に応じた水面変化は、tutorialケースを修正することで容易に確認できました。無理な条件を設定すると収束しませんし、取水深度を考慮する必要があります。水位ではなく水頭コンターが表示されますので、同じ平面位置でも、どの深度に観測孔を設けるとどのような観測水位になるかも想像しやすい。EXCEL よりはるかに多くの情報を与えてくれます。

やはり、その時代に即したツールを使う方が効率的です。

技術者の階段 その2

数年前にお世話になった御客様から電話を頂きました。

どうも、御一緒した仕事の続きをされるみたいでした。
当時、地質統計学を用いた3次元確率分布の可視化を含め、好評だったと記憶しています。覚えていてくださったのはありがたい話です。で、今の担当部署にバトンタッチ。
結果、executive が「こんなことはできないよ」と一蹴されていました。担当部署では技術的に備えておらず、一般的な調査のみで仕事を片づけられたようです。私も当時使っていたソフトは部署移動時に手放してしまいましたのでサービスできず。残念です。

別部署の後輩さんとの雑談で、水位と透水係数の話が出てきました。

水位を変えるために透水係数をどう設定するかが分からないようでした。先輩さんの専門が異なるので、地下水に関する十分な指導を受けていないのでしょう。
昔 EXCEL で作っていた揚水時の2次元水位変化の book を渡しておきました。単純化しているので正確ではないのですが、透水係数や揚水量に応じてどのような水面が得られるかイメージを掴む程度のことは可能でしょう。
本来なら基礎理論を理解し、水理公式や FloPy などでイメージを掴んでから浸透流計算に移るべきです。が、後輩さんの部署は基礎理論はおろかプログラミングも計算もされないようですので EXCEL ぐらいしか選択肢がありません。残念です。

いずれも単純な技術力不足なのですが、それで通用するのがヒトの世。どこかで技術者の階段を必死に登ることをやめたヒトたちが集まって、バイアスがかかり「常識」が生まれ、技術力不足を感じなくなったのでしょう。怖いですね。

しかも、このような状況に陥らないようにするために技術士の総監があったはずです。専門技術のみならず、管理技術を勉強される方が増えると、お客様にとっても喜ばしいことです。そう思えば日々研鑽される方が出てくるのではないでしょか。

一人でも多く、そのような「技術者」と呼べる方が増えることを願います。

CfdOF その2

CfdOFに関しては、多孔質媒体を含む流体解析が目的。

OpenFOAM では Darcy Forchheimer（ダルシー・フォークハイマー）model として実装されています。
https://openfoamwiki.net/index.php/DarcyForchheimer

日本語の文献もありました。　
高田ほか（2019）防波堤基礎捨石より流入する津波や高潮に関する水理実験および数値解析
吉岡ほか（2010）高透水性多孔質体中の非ダルシー流れに関する考察

係数は Ergun の推定式が良いとされていますが、透水係数からでも求められます。今回は一つ目の文献と同条件で多孔質媒体の透水性を定め、テストしてみました。

VOF（Volume of Fluid）法での2相流です。越流からのエアのトラップも表現されていますが、メッシュをもっと細かくしないとダメですね。

粒径を半分にした場合がコチラ。

透気係数や空気の圧縮性が気になりますが、設定を見つけられませんでした。表面張力の設定は transportProperties の中にありました(sigma)。後からエディタで修正が必要でしょう。

OpenFOAM の基本は tutorial をコピーして利用するスタイルかと思いますが、こちらの方が便利。自由表面の流れと地下水を同時に扱う場合、これをベースに設定を詰めていくスタイルの方が楽そうです。

*************************
20220205追記

式を展開していくと、高田ほか（2019）の(12)式とは微妙に異なる形になりました。数値は一致したので、誤植でしょうね。