Research

　潤滑や冷却を目的として気液二相流は工業機器・製品でよく用いられており、近年の環境意識の高まりからこれらの機器の高効率化が強く望まれています。気液二相流現象を解析することにより、機器の最適な設計を行うことによって損失を最小限に抑えることができます。

　歯車等の回転体に駆動される気液二相流は気液界面が大きく変形しながら分裂・合体を繰り返し、液相の分布は時々刻々複雑に変化します。このような流れ場を数値計算によって明らかにすることは、現象の複雑さゆえに非常に困難とされてきました[1]。我々は回転する円筒座標系上で記述された非圧縮Navier-Stokes方程式を解くことにします。回転運動に特有のコリオリ力と遠心力も現れます。気液二相流は気相・液相で物性値を切り替える一流体モデルを用いました。非圧縮条件を満足させるために圧力のPoisson方程式に対する疎行列を解く必要がありますが、気液二相流では気体と液体の密度が1000倍も異なるため、単相流の場合に比べて疎行列の反復計算が収束し難いことが知られています。本計算では、みずほ情報総研（株）製の行列ソルバライブラリ(MG-BiCGStab法)を用いることで、大規模な疎行列を高速かつ安定に計算しています。

　気液界面追跡法にはCLSVOF (Coupled Level Set and Volume of Fluid) 法を用いました。対象とするスケールでは気液界面は不連続と見なせるのですが、数値計算の安定性からVOF関数に遷移幅を持たせており、気液界面は有限な幅を持ったものとしてモデル化されています。また、高次精度の移流スキームを用いたとしても、VOF関数の移流を繰り返すうちに微小な誤差の積み重ねによって気液界面がぼやけてしまうのですが、VOF関数の再構築を行うことによって気液界面幅を常に一定に保つことができます。気液界面の単位法線ベクトルは急な勾配を持つVOF関数からではなく、滑らかな勾配を持つ符号付き距離関数（Level Set Function）から気液界面の単位法線ベクトルを算出しています。

　TSUBAME2.0上で128 CPU コアを用いて回転体に駆動される気液二相流を計算し、日産自動車（株）が実施した実験と比較を行いました。図1から分かるように、気液界面形状について計算は実験を良く再現できていることが確認できました[2]。

図1　実験と数値計算での気液界面形状の比較

また、撹拌抵抗値についても、図2から分かる通り実験と計算の良い一致が得られています[2]。

図2　攪拌抵抗値の実験と数値計算の比較

さらに、図3のように回転体にかかる圧力とせん断応力分布を解析することによって詳細な現象の理解も可能となり、図4に示すように圧力とせん断応力それぞれに起因する攪拌抵抗の時間変動を明らかにすることができました。図3と4のA〜Dは時刻が対応していて、回転体の歯面がオイルに突入する際に、歯面に大きな圧力がかかるため圧力起因の攪拌抵抗がピークを迎える一方、せん断応力起因の攪拌抵抗は比較的変動が小さいことが分かりました。全攪拌抵抗値は実験で計測できますが、せん断応力と圧力の攪拌抵抗への寄与度は数値計算により始めて明らかになったことであり、これらを機器の設計に反映できるなど産業への貢献が可能になります。

図3　回転体表面にかかる圧力とせん断応力の分布
(図をクリックすると拡大されます。)

図4　せん断応力と圧力に起因する攪拌抵抗の時間変動

参考文献