こんにちは(@t_kun_kamakiri)
皆さん、流体解析をやっていて「発散した!」って慌てることがよくありますよね。
むしろ、よくありすぎてもはや慌てませんかね。
逆に発散したときに「また、おれのノウハウが蓄積される」とワクワクしてきますよね(´;ω;`)ウッ…
本記事ではまだまだ未熟ながら流体解析に携わって、計算が発散した場合の原因と対処方法をいくつか例を出しておきたいと思います。
流体解析で計算が発散してしまった場合の原因と対処方法の一般的な話をまとめます
普段はシミュレーションとしてOpenFOAMを使用しているため、OpenFOAMに沿った内容になりますが商用の解析ソフトでも同じ考えが適用できることでしょう。
また、気付いた点は随時アップデートしていくため、計算が発散した際の対処方法はこれがすべてではありません。
他の記事でもいくつか流体解析において計算を安定させるための指針を示してくださっている方がいます。
- 計算の安定化のための指針
- FOAM FATAL ERROR: Maximum number of iterations exceeded
- 【CATCFDzero Tips】計算が発散する・収束しない場合の対処法
過去に自分が熱流体で計算が発散した場合に試したことをこちらにまとめておきました。
【OpenFOAM】熱流体のチュートリアルで発散した。改善の着眼点。
GoogleグループでもOpenFOAMに関して自由に質問することができ、有志で質問に答えてくれる方がいらっしゃいます。そのおかげで過去に起こったエラーはこちらで検索を掛けるととてもためになるやり取りをされている場合があります。
以下参考になりそうなものを添付しておきます。
単純な設定ミスや、設定はあっているが物理的に成り立たな設定になってしまっているものまで原因は様々です。
こういった”体系的な数値解析における誤り”をまとめるようと思うと少々力んでしまうので、本記事では思いつくままにどんどん追加していきます。
原因と対処法を言語化しておくことで指針が固まってくるものです。
人に質問したり調べたりする際には、よく観察して言葉にするという工程が必要です。
知識不足であるがために、原因まで深堀できないということもありますが、単に質問が壊滅的に下手であるがゆえに原因と対処法を深堀できていないという場面もしばしば見受けられます。
このあたりを一読しておくとよいでしょう。
設定のミス
流体解析で計算を実行した初期からエラーで計算が止まってしまう場合があります。
その場合はエラー内容を見て落ち着いて対処すれば良いです。
そもそも、初めての計算で一発で計算がうまくいくなんてことの方が稀なので、エラーが出たら「ラッキー🤗」って思って良いです。
エラー文にエラー内容が書いているので、修正はたやすいです。
この時点でのエラーは以下のものがあります。
- 境界条件とメッシュのタイプがあっていない
壁条件で設定する場合は「constant/boundary」のパッチタイプはwallでなければいけませんが、タイプが間違っている。 - メッシュ数と内部領域のリストの数があっていない(interFoamでありがち)
interFoamでsetFieldsを行うとalpha.waterに内部値が割り当てられ、値がリスト化されます。その際に内部状態のリストの数は決まっていますが、仮に後でメッシュを細かくしたいと思ってblockMeshでメッシュを切ってしまった場合、メッシュの数と内部値の情報の数が合わずにエラーが生じます。 - setFieldsをしているが0ディレクトリがない
OpenFOAMの場合、すぐに初期設定に戻せるように「0.orig」というディレクトリで境界条件を設定し、「cp -r 0.0orig 0」で0ディレクトリにコピーして進めますが、そのことを忘れてしまっている場合があります。
setFieldsコマンドでもエラーが生じるので注意が必要です。 - 解析設定が間違っている
ほとんどOpenFOAMの内容になってしまいましたが、とにかく初歩的ミスでエラーになっていることもよくあるということです。
エラー文をよく見て対処しましょう。
計算途中で発散する
計算が途中で発散してしまう場合は一番厄介です。
発散の原因
計算が発散してしまう原因を思いつくままに書いてみました。
ここから対処方法も見えてくるかもしれません。
- メッシュ品質が悪い(checkMeshで事前に確認する)
メッシュ品質の基準は、こちらの記事で「メッシュ品質」と検索して確認してみてください。 - 温度が負になる
メッシュ品質が悪い、境界条件がおかしい(全て流速固定するなど)、緩和係数や時間刻みが大きすぎる。
まずはメッシュ品質を疑う。 - 流速が速すぎる(断面積が小さいと流速が速くなる)
- 外部流出の場所で逆流が起こっている
外部流れを計算する場合は、解析領域を十分広くとる必要があります。
こちらに解説があるように物体の5~10倍くらいの広さを目安に流出領域をとる。
チュートリアル:簡易車体モデルでの抗力係数の測定 – XSim - 乱流モデルの境界条件がおかしい
- AMI系の境界条件でペア面のメッシュサイズが違い過ぎる
boundary設定のマッチングトレランスが合わなくて(たぶん連続式が満たされなくて)発散する - DyM系のソルバーだと計算途中のメッシュの変形で歪みが大きくなって発散
(Overset系ソルバーを使うなど別で対処)
エラーに対して対策は一対一に対応しておらず原因は複合的なので、ここでは「これに対してはこれ」といった明確な対処方法を示すことはしません。
対処方法
上で挙げた原因に対する具体的な対策というよりは、一般的に計算が発散した際にすることを以下に列挙します。
- メッシュ品質を見直し、メッシュを作り直す
- 緩和係数を小さくする(定常解析の場合)
- 時間刻みを小さくする or クーラン数を小さくする(非定常解析の場合)
(主流流速)✕(時間刻み)<(メッシュサイズ)を基準にしておく - 対流項スキームを変える(1次風上差分のupwind、TVDスキームにする)
スキームを変えるだけでこれだけ結果が変わりますよという例をこちらの記事に書いています。 - 安定的な行列ソルバーを使う
- 内部初期状態を確認する
定常解析の場合、内部の初期値は適当でも良いが境界との差がありすぎると問題になる - 非定常解析に変えてみる
そもそも非定常の現象なのにむりやり定常解析で行おうとして安定していない場合もある - 乱流モデルを使う
乱流モデルを使っていない場合、極端に粘性係数を小さくしてみると計算がうまくいく場合がある。粘性を小さくした際に計算ができるなら、乱流モデルを入れることで渦粘性が効いてくいているので乱流モデルを使うと安定する可能性がある。 - そもそも乱流モデルを使った方が良いのか検討する
層流なのに乱流モデルを入れると渦粘性は計算に入るので、こいつが悪さをする場合がある。レイノルズ数、グラフホフ数で乱流状態かどうかを見積もるのも大事 - 流速と圧力の両方を固定値にしていないか
流速が固定値なら圧力は0勾配
圧力が固定値なら流速は0勾配
とするのが良い。一つの境界でどちらも固定値にした場合に、拘束が厳しすぎて途中で計算が破綻する場合がある。
境界条件については「流体解析の境界条件について」が参考になるでしょう。 - メモリ不足
メッシュ生成でのkilledを起こす場合がある。
よくわからずエラーになるのでびっくりするが、単なるPCのメモリ不足。
Linuxであれば「free -h」で調べることができます。
1 2 3 4 | $ free -h total used free shared buff/cache available Mem: 31Gi 4.3Gi 26Gi 0.0Ki 279Mi 26Gi Swap: 100Gi 0B 100Gi |
メモリ不足→実メモリを増やす、メッシュ数を減らす、swapして仮想メモリを使うなどで対処する。
物理的におかしい
計算はエラーを起こさずに物理的に不自然な現象を起こすこともあります。
ここからは個別に問題と解決策を順に紹介します。
- 自然対流で上側へ自然に抜けない問題
圧力境界条件をprghPressureでは改善(参考)
重力を考慮した気液二相流の出口境界条件について
浮力による流出条件
出口境界近すぎ問題もありそう。
変なエラーを起こさないために
変なエラーを起こさないために最低限確認しておく項目を書いておきます。
こちらに書いている内容は、エラーが起こらずたまたま計算が進んでしまった場合にも必ず確認しておくと良い項目です。
【計算前】checkMeshでメッシュ品質を確認する確認
以下のコマンドでメッシュ品質を確認します。
1 | checkMesh |
必要に応じてオプションをつけて実行します。
出てきた結果に対してメッシュ品質やメッシュ数を確認します。
ちなみにメッシュ品質の基準はOpenFOAMのコードの中に書いてあります。
例えば以下のコマンドでファイルのありかを調べます。
1 | find $FOAM_SRC -name "primitiveMeshCheck" |
以下にあることがわかります。
1 | /opt/OpenFOAM/OpenFOAM-v2012/src/OpenF/OpenFOAM/meshes/primitiveMesh/primitiveMeshCheck |
パスはお使いのバージョンによります。
viコマンドで中身を確認。
1 | vi $FOAM_SRC/OpenFOAM/meshes/primitiveMesh/primitiveMeshCheck/primitiveMeshCheck.C |
※$FOAM_SRC=/opt/OpenFOAM/OpenFOAM-v2012/src/
品質基準は以下を参照。
1 2 3 4 5 | Foam::scalar Foam::primitiveMesh::closedThreshold_ = 1.0e-6; Foam::scalar Foam::primitiveMesh::aspectThreshold_ = 1000; Foam::scalar Foam::primitiveMesh::nonOrthThreshold_ = 70; // deg Foam::scalar Foam::primitiveMesh::skewThreshold_ = 4; Foam::scalar Foam::primitiveMesh::planarCosAngle_ = 1.0e-6; |
メッシュ品質が悪ければ以下のコマンで、vtkファイルを作成することができます。
1 2 | foamToVTK -faceSet skewFaces -time 0 foamToVTK -faceSet nonOrthoFaces -time 0 |
ls constant/polyMesh/setsで品質の悪い項目のファイルを確認。
constant/polyMesh/sets/ファイル名
ファイル名 = skewFaces
paraviewでvtkを読み込む
- 高アスペクト比(Aspect Ratio)
非常に微細な境界層で現れる。
ソルバーの安定性にとって致命的なものではありませんが、収束速度を著しく低下させる可能性がある - 非直交性(Orthogonality)
3つの値の範囲を定義する
・n0 < 70 – 安全な値
・70 < n0 < 90 – fvSolutionのnonOrthoCorrectorやfvSchemesの数値スキームに対して特別な処理を行う必要あり
・n0 > 90 – シミュレーションに使用できない悪いメッシュ。メッシュの作り直し。 - 歪度(skewness)
値が大きいと結果の品質(正確さ)が損なわれる可能性がありますが、適度な大きさの歪度であれば、シミュレーションに使用することができる - Smoothness.
メッシュ間隔の最大変化は20%未満が理想
【計算中】確認する事項一覧
【計算中】残差を確認
残差も反復計算ののち収束していく値であるため、これが極端に大きな値になってしまう場合は途中でエラーになって計算がストップしてしまいます。
自分は連続式の誤差と残差は必ず確認するようにしています。
例えば、system/controlDictには残差や連続の式の誤差を出力する設定を記述しておけば、計算実行時に「postProcessing」に出力してくれます。
system/controlDict
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 | application simpleFoam; startFrom startTime; startTime 0; stopAt endTime; endTime 500; deltaT 1.0; adjustTimeStep no; maxCo 0.9; writeControl timeStep; writeInterval 50; purgeWrite 0; writeFormat ascii; writePrecision 6; writeCompression off; timeFormat general; timePrecision 6; runTimeModifiable yes; functions { continuityError1 { type continuityError; libs ( fieldFunctionObjects ); // Mandatory entries (unmodifiable)// Optional entries (runtime modifiable) phi phi; // Optional (inherited) entries writePrecision 8; writeToFile yes; useUserTime yes; region region0; enabled yes; log yes; timeStart 0; timeEnd 6000; executeControl timeStep; executeInterval 1; writeControl timeStep; writeInterval 50; } residuals { type solverInfo; libs ("libutilityFunctionObjects.so"); fields (U p k omega); } #includeFunc yPlus; #includeFunc CourantNo; } |
残差は必ず確認するようにしましょう。
【計算中】連続式の誤差を確認
OpenFOAMの計算実行時のログで以下を見かけると思います。
1 | time step continuity errors : sum local = 1.92864822495e-12, global = -6.75896580941e-13, cumulative = -3.1325829691e-08 |
連続の式の誤差
- sum local : 誤差絶対値の格子体積重み付け平均
- global : 誤差(符号あり)の格子体積重み付け平均
- cumulative : globalの累積
こちらは連続式の誤差です。
本来連続式は質量保存の観点から0になるのですが、これが数値的には誤差があり完全に0になることはなく、上記のログのように小さい値ですが0より大きい値を持ちます。
しかし、計算が発散する場合はこちらの誤差が溜まりにたまって発散している場合があります。もし、計算が発散している場合はこちらの値も注意深く見てみてください。
すぐに思いつく原因はとしては、
- 境界条件がおかしい
- メッシュ品質が悪すぎる
- メッシュの粗密が大きすぎ
- 緩和係数が大きすぎる(時間刻みが大きすぎる)
です。
とりあえず緩和係数(時間刻み or クーラン数)を小さくしてどうかはまず確認しています。
残差と連続式の誤差もデータとして出力されるので、可視化しておくと便利ですよね。
例えば以下のようにPythonで書いておくとデータをすぐに可視化することができます。
グラフを出力するPythonスクリプトはこちらです。
連続の式の誤差
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 | import numpy as np import matplotlib.pyplot as plt import pandas as pd import os def graph_func(df_continuityError, imin, imax, iposi,icol): plt.subplot(imin,imax,iposi) plt.plot(df_continuityError['time'],df_continuityError[icol],linestyle='solid',color='b',marker="o",label=icol) plt.xlabel('time') plt.ylabel(icol) plt.legend() plt.grid() def df_continuityError_func(dir_): data_continuityError = np.loadtxt(f'./postProcessing/continuityError1/{dir_}/continuityError.dat') df_continuityError = pd.DataFrame(data_continuityError, columns=['time','Local','Global','Cumulative']) return df_continuityError def df_concat(dir_List): df = pd.DataFrame() for dir_ in dir_List: try: df_ = df_continuityError_func(dir_) df = pd.concat([df, df_]) except Exception: print("Empty") return df dir_List = os.listdir('./postProcessing/continuityError1') df_continuityError = df_concat(dir_List) plt.figure(figsize=(20,3)) graph_func(df_continuityError, 1,3,1,'Local') graph_func(df_continuityError, 1,3,2,'Global') graph_func(df_continuityError, 1,3,3,'Cumulative') print('time : ',df_continuityError['time'].iloc[-1], 'continuityError : ', df_continuityError['Local'].iloc[-1]) df_continuityError.head() |
残差
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 | import pandas as pd import matplotlib.pyplot as plt import os def graph_layout(): plt.grid() plt.legend(loc='best',bbox_to_anchor=(1, 1)) plt.yscale('log') plt.xlabel('time', fontsize=16) plt.ylabel('residual', fontsize=16) def df_residual_func(dir_): df_residual = pd.read_table(f'./postProcessing/residuals/{dir_}/solverInfo.dat',skiprows=1) df_residual = pd.DataFrame(df_residual) return df_residual def df_concat(dir_List): df = pd.DataFrame() for dir_ in dir_List: try: df_ = df_residual_func(dir_) df = pd.concat([df, df_]) except Exception: print("Empty") return df dir_List = os.listdir('./postProcessing/residuals') df_residual = df_concat(dir_List) initial_residial = [data for data in df_residual.columns if "initial" in data] final_residial = [data for data in df_residual.columns if "final" in data] df_residual.plot(x="# Time ",y=initial_residial, figsize=(6, 4)) graph_layout() #df_residual.plot(x="# Time ",y=final_residial) |
残差の分布を出力す量にカスタマイズした有益な議論があります。
グラフ化の方法は他にはgnuplotを使う方法もあります。自身の環境と好みに合わせて好きなものを使えば良いかなと思います。
【計算中】y+
とある資料にメモした pic.twitter.com/nSUijF5Na2
— カマキリ🐲CAE頑張る (@t_kun_kamakiri) June 14, 2023
まとめると・・・
- 高Re型の標準k-εの場合は30<y+<300の対数則域を推奨
- k-ωSSTのようにy+広い範囲で適用できるモデルもある(1<y+<300 程度)
- y+<1の場合は低Re型モデル、壁関数を使用せず壁面のレイヤーを10~20層入れる
1 2 3 4 5 | yPlus { type yPlus; libs (fieldFunctionObjects); } |
「#includeFunc yPlus;」で指定すれば、sytemディレクトリにyPlusがあればそれを読み込みますが、もしなければ「$FOAM_ETC/caseDicts/postProcessing/fields/」からファイルを探してくるようになっています。
※$FOAM_ETC = /usr/lib/openfoam/openfoam2212/etc
yPlusの他にも出力設定ファイルは用意されているので、必要に応じて使えば良いです。
1 2 3 4 5 6 | ls $FOAM_ETC/caseDicts/postProcessing/fields/ AMIWeights MachNo R ddt fieldMinMax log randomise turbulenceFields writeCellCentres CourantNo ObukhovLength XiReactionRate div flowType mag randomise.cfg vorticity writeCellVolumes LambVector PecletNo add energySpectrum grad magSqr streamFunction wallHeatFlux writeObjects Lambda2 Q components enstrophy heatTransferCoeff processorField subtract wallShearStress yPlus |
【計算中】流量の整合性(収支)
流入量と流出量が釣り合っていないとおかしい。
全ての境界面を流入量と流出量で固定値とすると、境界条件に強い制限がかかるので計算が成り立たなくなる。どこかを自然流出にして、流出量を自然に決まるものとして計算するなどする。その際、流入量と流出量がつりあっていることを確認すると設定ミスに気付くことができる
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 | surfaceFieldValue1 //任意の名前 { type surfaceFieldValue; libs (fieldFunctionObjects); fields (phi); operation sum; regionType patch; name surfaceFieldValue1;//patch名 } //以下のようにすると設定を使いまわしできる surfaceFieldValue2//任意の名前 { $surfaceFieldValue1; name surfaceFieldValue2; //patch名 } |
バグによるエラー
オープンソースであろうと商用ソフトであろうとプログラムのバグは多少なりともあります。
特にバージョンアップした際のリリース直後にはバグが潜んでいることもあります。
なので、例えばver9.0が出た場合はver9.1くらいになってから使うなど、様子見が必要な場合もあります。
バグではないにしてもバージョンアップの際に、しれっとデフォルトを変えられている場合もあるので注意が必要です。
- リリースノートをよく読む
- 前のバージョンと同じ動作をするか確認
これくらいは最低限確認しておきます。
バグを発見したらベンダーか開発元に教えてあげると親切ですね。
計算が正常であることと実現象を再現しているかは別
計算がうまくいっても実現象を再現できているかどうかはまた別の話なのがシミュレーションの難しいところです。
CAEの目的ってなんだろうかというのを前に考えたことがあります。
シミュレーションは、
- 実験の代替えなのか?
- 実現象を完全に再現した便利なツールなのか?
- はたまた趣味なのか?
自分は趣味でシミュレーションしている面もありますが、シミュレーションの本来の目的は「製品開発の設計に活かすこと」であります。
製品開発には製品として機能するのかどうかの実験による検証があり、上手くいかない場合には「原因解析→対策→再試験」を繰り返します。
そもそも原因解析の場合に実験データだけでは限界があったりしますので、シミュレーションの出番といったところになります。しかし、シミュレーションが実験を全く再現していないトンチンカンなものだったら使うことができませんし、使いたくもありません。
そこで、シミュレーションが実験結果を再現しているか「CAEと実験とのコリレーション」という活動が必要になります。
実験がない場合は上記のようなベンチマーク試験を題材にして、自身の解析テーマに沿ってシミュレーションの妥当性を検証する必要があります。
シミュレーションもあくまで製品開発・設計のためのツールであり、使いこなせるかどうかは使っている人の意識の問題ということになるでしょう。
CAEと実験データとの比較には地道な検証が必要です。
以下に精度アップの検討材料を列挙しておきます。
精度アップの検討材料
- 定常解析と非定常解析
- 時間刻み or クーラン数
- 不足緩和係数(定常解析で発散する場合や残差が振動しすぎる場合など)
- メッシュの変更(ヘキサ、ポリヘドラルメッシュ)
- スキーム(特に対流項)
- 乱流モデル(他のRANS、LES)
- 流出位置
- 軸対称、2次元計算ではなくフルモデル
- ソルバの変更
- その他条件(壁関数の種類、VOF法の界面圧縮項など)
おすすめ参考図書
☟こちらは、OpenFOAMの日本語書籍が無い中唯一わかりやすい参考書だと思います。
☟こちらもOpenFOMの古いバージョンでの和訳になります。さすがにこちらはバージョンに対する日本語でのケアはしていないので、OpenFOAMに慣れている方は購入しても良いかと思います。僕は「日本語でまとまっている内容」なので少し重宝しています。
☟以下に、もっと初心者向けの同人誌を紹介しておきます。
初心者は「ってか、まずどうやってOpenFOAMをインストールするの?」というところからつまずきがちです。
そんな時は、以下の書籍をおすすめします。
改訂新版 OpenFOAMの歩き方 (技術の泉シリーズ(NextPublishing))
インストール方法とチュートリアルで流体解析を体験・・・ちょっと高度な解析まで解説があります。著者曰くOpenFOAMのバージョンの追跡を行いながら、書籍をアップデートするようなので安心ですね。
こちらはFoundation版ですが非常に細かくコードライブラリ書かれていて勉強になります。