もう1ヶ月以上前になるが、Dr.Jasak氏によるOpenFOAMのセミナーを受講してきた。
前半のほとんどのセッションはすでにワークショップなどで公開されている資料をベースのお話であったが、後半の「市販ソルバーからOpenFOAMへの移行操作」と、「OpenFOAMにおけるプログラミング」のセッションは多分初めてお目にかかったもので、とくに離散化方法（Discretisation Settings）について具体的なガイドラインまであって、帰ったら是非とも試してみたい・・・と思っていた項目について、ようやくの検証レポートです。
ちなみに、自分がこれまで問題に応じて変更したこのとあるスキームは、divSchemeくらいで、Laplacian schemes については、これまで馬鹿の一つ覚え的に、「linear」で良いんだろう・・・と思っていたのだが、実はそうでないということを教えられた。
メッシュの品質に応じて設定を変えたほうが良いということで、ヘキサメッシュ（snappyやポリヘドラルも含む）の場合、checkMeshを実行した結果のmax-non-orthogonality の値が60以下であれば、ほとんどのチュートリアルケースのデフォルト設定
Gauss linear corrected
を使ってよいが、70以上の場合は、
Gauss linear limited 0.5
を使いなさいというものであった。
現実の製品モデルなど複雑形状を対象としたメッシュでは、ほとんどの場合、60以下に収まることはない。そこで実際に適用してみたんだが、何と、これまで直ちに発散してしまっていたケースが安定して計算できるようになったという事例がたくさん出てきました。会社で扱っているケースの例はここで公開できないのですが、特に非定常問題で効果があったということくらいの言及はしておきます。 また、むやみに 0.5の値を小さくすれば、すべからく安定化するかというと、そういうことでもないので話は少々厄介なんですが・・・まぁ、計算が安定しない場合にはここを変更したスタディは大いにやってみる価値はあるということです。
ここでは、定常問題ですがDEXCSの標準テンプレートケースを使って検証した結果を紹介します。
DEXCSの標準テンプレートケースでは、DEXCSの３次元フォント周りの流れ計算をやっていますが、通常、ポテンシャル流れ計算によって初期流れ場状態を作成してからsimpleFoamによる流れ計算（非圧縮性乱流解析）をやらないと、安定した収束解が得られませんでした。このやり方をした場合の残差グラフと、流れ場の表示例を下に示しておきます。
これをポテンシャル計算なしで、いきなりsimpleFoam計算をすると、10数回のステップを経て発散してしまうか、下図の残差グラフで示すような一見安定して計算しているかのように見えて、とんでもない流れ場になってしまっています。
ちなみにcheckMeshの結果は以下の通り
Checking geometry...
Overall domain bounding box (-7 -5.5 -1.8) (1 -1.5 2.2)
Mesh (non-empty, non-wedge) directions (1 1 1)
Mesh (non-empty) directions (1 1 1)
Boundary openness (9.394e-18 1.40393e-16 -1.2532e-16) OK.
Max cell openness = 8.99745e-16 OK.
Max aspect ratio = 21.6098 OK.
Minumum face area = 1.53498e-05. Maximum face area = 0.160434. Face area magnitudes OK.
Min volume = 1.13903e-07. Max volume = 0.0641652. Total volume = 127.892. Cell volumes OK.
Mesh non-orthogonality Max: 64.8335 average: 8.61539
Non-orthogonality check OK.
Face pyramids OK.
Max skewness = 3.45097 OK.
Mesh OK.
問題のmax-non-orthogonalityは、64.8335 ということで、上に述べた基準値の中間的な値。
そこで、上述のGauss linear limited 0.5を使ったらどうなるかです。
イタレーションの回数は大幅に増えてはいるものの、ちゃんと収束して、流れ場もポテンシャル流れ場から始めた計算結果と比べて見た目での違いは区別できません。
なお、以上の結果は、まだ製作途上にあるDEXCS2011（OpenFOAM-2.0.ｘ）で行ったもので、モデルは基本的にDEXCS2010と同じものを使っていますが、フィーチャーエッジの影響などあって、上の結果はDEXCS2010でそのまま当てはまりません。また、DEXCS2010も、32bit版と64bit版、単体計算か並列計算かによって微妙に結果が異なってくるようです。多分、メッシュの品質が違うんだと思いますが正確なところはわかりません。
たとえば、DEXCS2010（OpenFOAM-1.7.x）の64bit版の場合、
Gauss linear limited 0.5 では発散
Gauss linear limited 0.33 とすると単体計算では収束方向の結果が得られたが２並列の計算では発散。良好基準値（60）はオーバーしているので、効果があったということでしょう。
Gauss linear uncorrecter ( Gauss linear limited 0 と同じ）
では、２並列の計算でも収束し、以下の結果が得られました。（20.png）
また、500から1000ステップあたりで残差がほぼフラットに推移する領域がありますが、十分収束しているとは言い難い状態で、
あれこれ解の方向を探して1300ステップ以降は一気呵成でした
ちなみに、checkMeshの結果は、
Checking geometry...
Overall domain bounding box (-7 -5.5 -1.8) (1 -1.5 2.2)
Mesh (non-empty, non-wedge) directions (1 1 1)
Mesh (non-empty) directions (1 1 1)
Boundary openness (-3.30528e-17 2.01901e-17 1.698e-17) OK.
Max cell openness = 6.75923e-16 OK.
Max aspect ratio = 18.0608 OK.
Minumum face area = 1.56449e-05. Maximum face area = 0.160451. Face area magnitudes OK.
Min volume = 1.52447e-07. Max volume = 0.0641757. Total volume = 127.893. Cell volumes OK.
Mesh non-orthogonality Max: 61.5977 average: 8.12569
Non-orthogonality check OK.
Face pyramids OK.
Max skewness = 2.023 OK.
Mesh OK.
となっており、non-orthogonality Max: 61.5977 で、DEXCS2011の場合と比べて、さほど品質が悪くない、ということもあって効果が今ひとつだったんでしょう。