OpenFOAMの計算安定化について

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

[Penguinitis](http://penguinitis.g1.xrea.com/study/OpenFOAM/stabilize_calculation.html)のページに既にまとめられているところですが、自分なりに得られた知見や経験も踏まえて、ここでもまとめておこうと思います。

# メッシュ

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

```c++
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`の組み合わせは、個人的な経験にはなりますが、リーズナブルな初期条件を吐き出してくれることが多いです。使用方法は[こちら](https://qiita.com/takaf05/items/04727e2cb6424812c81d)を参考ください。
- `potentialFoam`だけですと、速度を過大評価する傾向があります。

# 境界条件

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

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

```c++
// 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;
}
```


```c++
// 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`などを使用するのも良い。

```c++
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`

```c++: 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の設定を以下のようにする：

```c++
   nOuterCorrectors 7
   nCorrectors 2 // needed
   nonOrthogonalCorrectors 1 // enough
```


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