はじめに

OpenFOAMは支配方程式を素直に解こうとするためか、適切なモデルの選定、初期条件・境界条件の設定をする必要があります。

Penguinitisのページに既にまとめられているところですが、自分なりに得られた知見や経験も踏まえて、ここでもまとめておこうと思います。

メッシュ

checkMeshでエラーが出ていないことが必須です。

Checking geometry...
    Overall domain bounding box (-10.028 -93.2695 -686.041) (45.8037 10.5564 -599.5)
    Mesh has 3 geometric (non-empty/wedge) directions (1 1 1)
    Mesh has 3 solution (non-empty) directions (1 1 1)
    Boundary openness (5.43168e-17 -2.89688e-16 -3.6735e-16) OK.
    Max cell openness = 3.36631e-16 OK.
    Max aspect ratio = 9.27753 OK.
    Minimum face area = 0.000164425. Maximum face area = 1.8177.  Face area magnitudes OK.
    Min volume = 7.14396e-05. Max volume = 0.612806.  Total volume = 92689.5.  Cell volumes OK.
    Mesh non-orthogonality Max: 82.3126 average: 7.17106
   *Number of severely non-orthogonal (> 70 degrees) faces: 27.
    Non-orthogonality check OK.
  <<Writing 27 non-orthogonal faces to set nonOrthoFaces
    Face pyramids OK.
 ***Max skewness = 4.83182, 3 highly skew faces detected which may impair the quality of the results
  <<Writing 3 skew faces to set skewFaces
    Coupled point location match (average 0) OK.

Failed 1 mesh checks.

上記の例ではskewnessでエラーが出ています。このようなときはメッシュの再作成をするのが良いです。

初期条件

適切な初期条件を用意します。
potentialFoam + applyBoundaryLayerの組み合わせは、個人的な経験にはなりますが、リーズナブルな初期条件を吐き出してくれることが多いです。使用方法はこちらを参考ください。
potentialFoamだけですと、速度を過大評価する傾向があります。

境界条件

スキーム

安定性の観点から、解析最初50回程度を1次スキームを用い、その後、2次スキームに上げていくの良いです。
非定常解析を開始する際は、2次スキームを用いた定常解析結果を用いるのが望ましいです。
定常解析において、解析最初20-50回の反復計算をupwindで行うと、安定した流れ場を得られやすいです。その後、2次以降のスキームに上げていくと良いです。
linearUpwindは2次精度であり、乱流パラメーターや他のスカラー量の解析に適しています。
linearUpwindを速度計算に用いることは可能であるが、RANS系のモデルでは予測精度が低下するので注意が必要です。
limitedLinearは2次精度のスキームである。速度勾配が大きい場に向いています。
LUSTは linear (75%) + linearUpwind (25%) のブレンドスキームです。RANSや境界層流れ、DES、LESなどで精度良く解析が可能(ESIのセミナーではこのスキームが推奨されていた)。

以下のような組み合わせで安定することが多いです。

// time 0-50
divSchemes
{
    div(phi, U)      bounded Gauss limitedLinear 1;
    div(phi,k)       bounded Gauss upwind 1;
    div(phi,epsilon) bounded Gauss upwind 1;
    div(phi,omega)   bounded Gauss upwind 1;
}

// time 50-
divSchemes
{
    div(phi,U)      Gauss LUST unlimited grad(U);
    div(phi,k)       Gauss linearUpwind grad(k);
    div(phi,epsilon) Gauss linearUpwind grad(epsilon);
    div(phi,omega)   Gauss linearUpwind grad(omega);
}

（epsilonなどの乱流パラメーターの計算で発散することが多い。その時はQUICKなどの使用も検討すると良いかも・・・）

計算安定性を高めたい場合は LUST のほか、CoBlendedやlocalBlendedなどを使用するのも良い。

divSchemes
{
    div(phi, U)   Gauss CoBlend <Co1> <scheme1> <Co2> <scheme2>
}

---
Co < Co1 : scheme1
Co > Co2 : scheme2
Co1 < Co < Co2 : linear blend between two

線形ソルバー

(ESIのセミナーでは）PBiCGStabを使用することを強く推奨。
(PBiCGは使用しないこと！)
pre-conditionerには DILU を使用するのが良い。
tolerance 設定は重要
輸送方程式には relTol 0.1
圧力方程式には relTol 0.01
timeActivatedFileUpdateを使用すれば、異なる値を設定することが可能。
p-U couplingの最終ループでは relTol 0

fvSolution example

p
{
    solver  GAMG;
    tolerance   0;
    relTol   0.01;
    smoother   DICGaussSeidel;
}

pFinal
{
    $p;
    tolerance   1e-6;
    relTol   0;
}

"(U|k|nuTilda)"
{
    solver   PBiCGStab;
    preconditioner   DILU;
    tolerance   1e-8;
    relTol   0.1;
}

"(U|k|nuTilda)Final"
{
    $U;
    tolerance   1e-6;
    relTol   0;
}

その他

定常解析で不安定な場合は、PIMPLEを使うのが良い。
smoothSolver よりは PBiCGStab の方がロバスト性や安定性が上。
より精度良く解析したい場合は、PIMPLEの設定を以下のようにする：

   nOuterCorrectors 7
   nCorrectors 2 // needed
   nonOrthogonalCorrectors 1 // enough

まとめ

順次、更新していきます。

OpenFOAMの計算安定化について