[CFD Q&A] simpleFoamで圧力残差が収束しない時の対処法

先週、simpleFoamが半日を持っていきました。圧力残差は 1e-3 の下へ降りませんでした。上がりました。最初に疑ったのはメッシュです。実際に触った場所は relaxation factor でした。

αp のスライダーを 0.6 より上へ押すと、圧力残差の曲線が 30 iter 付近で再び上へ折れます。

0.3 付近に置くと、同じ曲線は 200 iter まで降りていきます。下のログの Time=80、初期残差 0.891 が、その折れと同じ区間です。

問題の状況

OpenFOAM v2312。2D の長方形チャネルで、inlet から outlet。乱流モデルは standard k-ε。ソルバーは定常の simpleFoam です。

最初の 20 から 30 iter までは残差が降りていました。そのあと、圧力だけが跳ねました。

Time = 50
 
smoothSolver:  Solving for Ux, Initial residual = 0.0023, Final residual = 1.2e-5, No Iterations 3
smoothSolver:  Solving for Uy, Initial residual = 0.0019, Final residual = 9.8e-6, No Iterations 3
GAMG:  Solving for p, Initial residual = 0.412, Final residual = 0.038, No Iterations 8
time step continuity errors : sum local = 4.2e-5, global = 1.1e-5, cumulative = 3.4e-4
...
Time = 80
 
GAMG:  Solving for p, Initial residual = 0.891, Final residual = 0.21, No Iterations 20

速度は収まり、圧力だけが変です。圧力-速度カップリングが揺れているときの顔です。

原因分析のプロセス

ステップ1：メッシュ品質の確認#

まず checkMesh を回しました。

checkMesh -case .

最大 non-orthogonality は 32 度、max skewness は 0.8。メッシュを犯人にするには静かな数字です。

ステップ2：境界条件の点検#

inlet の圧力は zeroGradient、outlet は fixedValue 0。チャネルではよくある組です。ここまでは異常がありませんでした。

初期 p を開くと、inlet 側に大きな値が残っていました。消したつもりだった funkySetFields の試験の跡です。

ステップ3：fvSolution の確認#

// fvSolution (元の設定)
relaxationFactors
{
    fields
    {
        p               0.7;
    }
    equations
    {
        U               0.9;
        k               0.7;
        epsilon         0.7;
    }
}

U=0.9 は速い設定です。速度の緩和を上げると、見かけの収束は速くなります。圧力の更新が追いつかないと、ループが揺れます。

重なった原因は二つです。緩和が高かった。初期 p にゴミが残っていた。圧力補正が速度場とどう結び付くかは、圧力が格子縞に揺れるときのSIMPLE の側にほどいて書いてあります。

解決策

リラクゼーション・ファクタを下げる

// fvSolution (修正後)
relaxationFactors
{
    fields
    {
        p               0.3;
    }
    equations
    {
        U               0.7;
        k               0.5;
        epsilon         0.5;
    }
}

最初のランは低く始めます。残差が座ってから上げます。

初期 p フィールドの初期化#

# 0ディレクトリのpファイルをきれいに初期化
cp 0.orig/p 0/p

0/p を直接開いて internalField uniform 0; に合わせても同じです。

// 0/p
dimensions      [0 2 -2 0 0 0 0];
internalField   uniform 0;
 
boundaryField
{
    inlet
    {
        type            zeroGradient;
    }
    outlet
    {
        type            fixedValue;
        value           uniform 0;
    }
    walls
    {
        type            zeroGradient;
    }
}

検証

同じケースをもう一度回すと、ログはこう変わります。

Time = 50
 
smoothSolver:  Solving for Ux, Initial residual = 0.0018, Final residual = 8.3e-6, No Iterations 3
GAMG:  Solving for p, Initial residual = 0.089, Final residual = 0.0041, No Iterations 5
time step continuity errors : sum local = 2.1e-6, global = 4.3e-7, cumulative = 8.9e-6
 
Time = 200
 
GAMG:  Solving for p, Initial residual = 3.2e-4, Final residual = 1.1e-5, No Iterations 3

200 iter 付近で圧力残差は 1e-4 の水準まで降りました。内部反復も一桁です。

ログが示す倍率

直す前の Time=50 で、Ux は 0.0023 から 1.2e-5 へ落ちました。内部反復 3 回で約 192 倍です。Uy は 0.0019 から 9.8e-6 へ、約 194 倍です。同じ時刻の圧力は 0.412 から 0.038 まで、8 回で約 11 倍だけです。速度は約 200 倍近く降りたのに、圧力は約 11 倍で止まっています。

Time=80 は、そこから 30 ステップあとです。圧力の初期残差は 0.412 から 0.891 へ、約 2.2 倍になりました。内部反復は 8 から 20 へ 2.5 倍増えました。ところが減少倍率は 0.891 割る 0.21 で約 4.2 と、かえって減りました。反復を増やして、あまり降りていない記録です。

直したあとの Time=50 の圧力は 0.089 から 0.0041 へ、5 回で約 22 倍降りていきます。Time=200 では 3.2e-4 から 1.1e-5 へ、3 回で約 29 倍です。累積 continuity は 3.4e-4 から 8.9e-6 へ、約 38 倍小さくなりました。p の緩和 0.3 は、以前の 0.7 の約 43% です。U は 0.9 から 0.7 へ、k と epsilon は 0.7 から 0.5 へ下げました。圧力補正を半分より下へ削った代わりに、圧力の内部反復は 20 から一桁へ減りました。

Time=200 の最終残差 1.1e-5 を、直す前の Time=80 の最終値 0.21 と並べます。0.21 割る 1.1e-5 は約 19100 です。同じチャネル、同じメッシュで、変えたのは圧力の緩和と初期場だけです。メッシュを再び開ける前に見る数字は、この倍率です。

速度の内部反復は、直す前も後も 3 のままです。変わったのは圧力側の反復数と減少倍率だけです。Ux の最終残差も 1.2e-5 から 8.3e-6 へ、同じ 3 回の中で少し余計に降りただけです。ログが指す故障は速度ソルバーではありません。continuity の sum local も 4.2e-5 から 2.1e-6 へ減りました。4.2e-5 割る 2.1e-6 は 20 です。累積の 38 倍と同じ方向で、倍率だけが別です。

メッシュを開く前に下げる値

non-orthogonality 32 度と skewness 0.8 は、メッシュの作り直しから入る数字ではありません。先に下げる値は p の緩和です。0.7 は、圧力補正が速度 0.9 に追いつけなかった設定でした。最初のランは p 0.3、U 0.7 で見ます。Time=200 の残差が 1e-4 付近に座ってから上げます。

初期 p に残った funkySetFields の値も、同じランを壊しました。internalField uniform 0 で空にしないと、緩和を下げても入口のゴミ勾配が残ります。緩和を下げた解が以前と同じ物理かは、別に確認がいります。緩和を一つ変えただけで収束解が動いたRhie–Chow を見ると、0.7 を 0.3 にした選択が残差曲線の外まで残るかが、次の質問です。