2022年2月19日土曜日

メッシュ変換

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 のプラグインもあったのに惜しい!

次。

Gmsh
https://gmsh.info/

こちらはOpenFOAM で利用していました。フリーのメッシャーとしか認識していませんでしたが、色々なフォーマットの入出力に対応していました。メッシュ変換ソフトとしても利用できそうです。
なぜか1台のPCでは export できなかったのですが、別のPCでは問題なく変換できました。感謝です。

どこかの商用メッシャーよりも優れています。使用しているメッシャーの種類を聞けば、その方々の技術力もある程度ふるい分けできそうな予感。

 

Geo-Graphia

地層研さんの 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次元可視化ソフトの位置付けなのでしょう。

やはり私は使用することはないでしょうし、お勧めもしません。が、それを再確認出来た点はヨカッタ。
先輩さんの部署では元を取るまで使っていかざるを得ないでしょうね。オプションも含めると結構な金額でしたので。騙されたと思われるか、自分の目利きを後悔されるかはわかりませんが、いずれにしてもお気づきになるまで相当の時間がかかるでしょう。

2022年2月15日火曜日

Nays2D+

次は Nays2D+ Ver1.0.8。

こちらは iRIC Ver.3 から登場とのこと。平面2次元の計算にて疑似的な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 の ポスト処理としてBlueKenue にも似たような機能が搭載されていました(こちらは sediment transport)。案外、広く使われているのかもしれません。

現段階では境界条件として時系列の濃度変化までは設定できないようですし、吸着も扱えません。が、そのうち実装されるでしょう。コチラも気長に待ちましょう。


2022年2月14日月曜日

Morpho2DH

iRIC Ver.3 に含まれている土石流解析ソルバーの Morpho2DH を使用。Ver.1.0です。
ソルバーマニュアルには、「Morpho2DH は平面二次元の土石流・泥流モデルを主体とした解析ソルバであり,斜面崩壊を初期条件とした土石流・泥流の流動・堆積過程を表現可能なモデルです」とあります。

まず、tutorial。前半は問題なく進みましたので、動くようです。
プレゼンで示されていたように、非常にお手軽に計算ができます。tutorial では xyz の csv データが用意されていましたが、iRIC のプレは基盤地図情報からのダウンロード機能を備えています。座標系を指定しておけば背景に地理院地図等を表示できますので、それを見ながら DL する範囲の指定が可能。範囲内に 5mLP がなければ、ある部分だけで点群を作ってくれます。この一連機能は秀逸。LandXML への吐き出しも実装されているようですので、土石流や河床変動以外の解析でも、この部分だけ利用する場面が出てきそうです。

次は、後半の固定床を飛ばして、実例をトレースすることに。
5mLP をインポートし、流下方向を考慮しながらメッシュを切る。空中写真を背景に頭部崩壊領域を作成し、メッシュに属性を付加、計算条件を入力したら設定完了です。
今のところ GUI は、こけていません。たまに属性のマッピングが反映されていないケースに遭遇しましたが、どういうタイミングで発生するのか掴めていません。念のため保存して終了し、再起動後に読み込んで入力値が入っていたらOKです。同じ類の Hyper-Kanako に比べ、操作性・安定性はマシマシです。

で、計算!

はい、結果が出ません。すぐに発散しています。ソルバーの中にタイムステップを自動調整してくれる機能がないようです。しかも並列化されていないので、計算が遅い。
タイムステップを 0.01 秒から 0.001 に切り替えてもダメ。
メッシュを 10m から 5m に切りなおすと、動きました。

ひとまず一つの流路では動いたので、同時多発を試行。地理院の空中写真を背景に切り替えて10か所超の頭部崩壊領域を設けました。発生時間をずらすように設定し、計算スタート!

結果はそこそこ。
建物の影響を考慮していない地形ですので実際の堆積域とは一致しません。が、計算は走りました。力技ですが、メッシュを細かく切り「障害物」として建物を指定すると、その影響は取り込めるようです(さらに計算は遅くなりますが)。

iRIC のポスト処理として、地理院の空中写真を背景に層厚や速度を重ねることができます。動画も書き出せます。操作が簡単で安定していることを含め、フリーのプリポストとしては作りこまれている感じです。

Morpho2DH としては、同時多発の土石流を一度に計算できることが特徴でしょう。改善要望としては、並列化と建物の崩壊を含めた影響検討ですかね。いつか実装されるかもしれませんので、気長に待ちましょう。

2022年2月11日金曜日

iRIC Ver.3

一週間ほど、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次、風上差分:1次を実装しているようですが、前者で振動するなら後者で、だそうです。

実際の地形を考慮するのは5章以降。引き続き読んでみましょう。

*******************************
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/

 

2022年2月5日土曜日

Darcy-Forchheimer model

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:圧力損失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]+Bu2[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:係数[s2/m2], d:粒径[m], φ:間隙率[-], ν:動粘性係数[m2/s], g:重力加速度[m/s2]

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

2022年2月4日金曜日

OpenFOAM の error

集水地形に雨を降らせて、地下水と表流水の変化を見てみることにしました。

といっても、まだ実務レベルではありません。先日の CfdOF (OpenFOAM) で、どこまで手間をかけずに計算できるか?という確認レベルです。

まずは、モデル作成。

  1. Infraworks のモデルビルダーで適当な集水地形を選択し、地表面サーフェスを作成
  2. Civil3D でモデルを取り込んで、サーフェスから下側のソリッド、それを含む全体のソリッドの2つを作成
  3. iges形式で書き出し
  4. ついでに下側だけ stl 書き出し(ParaView用)
  5. FreeCAD でiges をインポート
  6. CfdOF でメッシュ作成、条件設定
  7. 保存

あとは実行。
ですが、途中で止まります。やや乱暴な境界条件でしたので「それが原因かな?」と思いつつ、それを修正するにはジオメトリ修正からなので後回し。簡単に修正可能な設定から見直すことにしました。

「OpenFOAMによる熱移動と流れの数値解析 第2版」と著者のサイトを参考に、3日ほどあれこれ試していたら、走るようになりました。効果的だったのは、以下の3点。

  • 代数方程式のソルバーを GAMG/smoothSolver から PCG/PBiCGStab に変更 p92
  • PIMPLE法で nOuterCorrectors を1から2へ変更 p101
  • メッシャーを cfMesh に変更 

ひとまずこれでOK。出てきた結果がコチラ↓


うーん、コレジャナイ。いえ、最初だから砂利山のような高透水を設定したので、これで正解なのですが。 

透水性を下げてみると、

それらしくなりました。

ひとまず、エラーは回避できました。