今回は流体解析などで使われる緩和係数についてお話します。
目次
緩和係数とは
有限体積法の流体解析を行う場合、緩和係数を設定します。SIMPLE法などの反復解法において、ある計算ステップで計算された物理量をそのまま使って、次のステップの計算を行うと計算が不安定になることがあります。そこで、次の式のように、前のステップの値を加味することで安定に計算できるようにします。
$$x_C = x_C^* + \lambda (x'_C - x_C^*) \tag{1}$$
ここで、$x_C$は次のステップの値、$x_C^*$は前のステップの値、$x'_C$はそのステップで計算された値、$\lambda$は緩和係数。
式を見ると、緩和係数は前のステップとの変化率のようなものだということが分かります。
$0 < \lambda < 1$であれば、不足緩和（under-relaxation）と言います。$\lambda$が０に近いと、前のステップからの変化が少なく安定に計算できます。ただ、最終的に解が収束するまでの計算ステップが多くなり、計算時間がかかります。
$\lambda = 1$であれば、緩和はかかりません。そのステップで計算された値がそのまま次のステップで使われます。
$\lambda > 1$の場合は、加速緩和（over-relaxation）と呼ばれます。問題によっては収束が早まることがありますが、流体解析の一般的な問題ではほとんど使われません。流体解析では普通、不足緩和が使われます。
陰的な緩和方法
緩和は式上は(1)式のようになっていますが、解法上は陰的な方法が取られます[1]。
有限体積法の線形方程式は、$C$をそのセル、$F$を隣接セルとすると、
$$a_C x_C + \sum a_F x_F = b_C \tag{2}$$
となります。詳細は以下を参照。
(2)式より
$$x_C = \frac{b_C - \sum a_F x_F}{a_C} \tag{3}$$
を(1)式の$x'_C$に代入して整理すると、
となります。これは(2)式で、$a_C$ が $a_C/\lambda$ に、$b_C$ が $b_C + ((1-\lambda)/\lambda) a_C x_C^*$ に置き換わった形をしています。つまり、(2)式の代わりに(4)式を解くと、緩和された式を解いたことになるのです。
また、緩和係数$\lambda$は不足緩和（< 1）であることから、係数行列の対角成分 $a_C/\lambda$ は $a_C$ に比べ大きくなり、行列解法において安定に働きます。
緩和係数の値
だいたいよくある流体解析ソフトでは、緩和係数の値は、SIMPLE法による解法では、流速で0.7、圧力で0.3程度になっているのではないかと思います。実際の値は、アルゴリズムやソルバーによって違いがありますが、デフォルトではそのプログラムで最適な値が入っていると思われます。
※当サイトのオンライン流体解析ツールCATCFDzeroでもデフォルトで流速0.7、圧力0.3の緩和係数になっています。
ただし、問題によって緩和係数の値は変える必要があります。例えば、メッシュが歪んでいたり、境界条件や物性値などの条件が厳しい問題では、デフォルトの緩和係数では収束しないことがあります。また、収束しないだけでなく、発散して計算が破綻してしまうこともあります。
そのような場合は、緩和係数の値を小さくします（0.7から0.5などに変更）。こうすることで安定に計算できる場合が多いです。ただ、やみくもに小さくしすぎると、いつまで経っても収束せず計算時間がかかってしまいます。計算初期で発散傾向にある場合は、最初は緩和係数を小さくして、安定に計算できるようなら徐々にあげていくなどの対処をすることで、計算時間の増大を防いだりします。
まとめ
今回は、緩和係数についてお話しました。緩和係数は計算が不安定であったり、発散したりする場合に、まずはじめに調整するパラメータです。自分が使うソルバーで、緩和係数を変化させた時に計算の履歴がどのように変わるか確かめておくと、いざというとき調整しやすくなるので、いろいろ試してみてください。