By Jalil Ouazzani, ArcoFluid ? Bordeaux, France, Zaim Ouazzani, Arcofluid Consulting LLC,
Orlando, USA ? John Ludwig CHAM Ltd, London, UK.
はじめに
ここでは、混ざりあわない3つの異なる流体による三相流をシミュレートするための、既存のPHOENICS VOF(Vomume of Fluid) オプションの拡張について報告します。例としては、油水分離器の水/油/空気界面や、ガス撹拌冶金取鍋の液体鉄鋼/液体スラグ/ガス界面に見られるような、密度の異なる液体と気体の組み合わせを含むシステムが挙げられます。3 つの流体への拡張には、第3相を表す追加のカラー関数C3 の保存則を解くことが含まれます。3つの異なる流体を処理するために、混合物の運動量方程式における表面張力を修正します。
PHOENICS 3相 VOF の実装は、気泡上昇について２つのケースでテストされました。１番目は単一液柱内で上昇する気泡の問題、２番目は2つの液体が層状に液体を含む塔内で上昇する気泡の場合です。最初のケースは、VOFの拡張で、単純な2相システムについて公表されている数値結果と一致する結果が得られることを検証するために重要です。2番目のテストケースは、完全な三相システムをシミュレートして、各相の物性の変化に応じた気泡の上昇を定性的に再現できるかか調査することです。
理論的考察
C3の輸送方程式は、2相VOFシミュレーションで使用される既存のカラー関数C1と同じ形式です。次のようになります。
DC3/Dt=0 (1)
PHOENICSには、物理的な問題に応じて、(１)式を保存形式もしくは非保存形式で解くオプションがあります。∑ Cn = 1とすることで体積の連続性を強制し、2つのカラー関数を連結させます。ここで Cn はn相のカラー関数です。物性はΦ = ∑ Cn Φn によって平均化されます。Φ は、密度、動粘度、比熱、熱伝導率、および体積膨張係数です。
PHOENICSの2相 VOF法では、Brackbill et al (1992) の標準的な連続表面力 (CSF)アプローチを使用して、表面張力を等価物体力fcap = σ ki δ niの形で運動量方程式に導入します。
ここで ni = -∇Ci/|∇Ci| は界面から外向きの法線ベクトル、Ciはi相のカラー関数、δ = |∇Ci| は界面を中心とするディラック デルタ関数、ki = -(∇・ni) は界面の曲率です。CSF アプローチの欠点は、隣接する相の密度が異なると、運動量方程式に導入された表面張力により、界面の位置に対して加速場の非対称な分布が生じることです。たとえば、加速度fcap /ρ (ρ はローカル VOF 相の密度) は、密度が低い位相でははるかに高く、その逆も同様です。CSF アプローチでは、ベクトル fcapの方向に応じて、相間の滑らかな移行領域が薄くなったり厚くなったりします。この問題は、Brackbill et al (1992) によって、次のように CSF (DS-CSF) の密度スケーリングを使用することにより、二相システムに対して解決されました。
fcap = - σ ki ∇Ci ρ/<ρ> (2)
ここで、<ρ> = (ρ1 + ρ2)/2 は、隣接する第１相と第2相の平均密度です。これにより、界面に関して対称的な加速度分布が得られます。この研究では、Tofighi and Yildiz (2013) に従って、結果として生じる表面張力を3つの構成要素 (各相に1つ) に分割することにより、DS-CSFを3相に拡張しました。これらの相固有の力はそれぞれ、上記の式 (2) によって与えられますが、界面表面張力を使用する代わりに、3つの相に固有の表面張力が使用されます。これらの力にはσn (n = 1,2,3) が使用されます。後で説明するように、このアプローチは三相システムにのみ有効です。特定の位相 n に注目する場合、密度スケーリングの考え方は、他の 2 つの位相を、空間的に密度が変化する単一の n 隣接位相として扱うことです。二相系と同様に、DS-CSF の n隣接相の密度を使用すると、三相系の表面張力は次のように計算できます。
fcap = ∑ fn,cup = - ∑ σn kn ∇Cn ρ/<ρ>a (3)
ここで、fn,cap は、<ρ>n = (ρn + ρn_adjacent)/2 の n相 の fcap と等価です。
この定式化は、対称的な加速度を生み出すように、界面全体に表面力を再分配します。相固有の表面張力の値を定義することはまだ残っています。このアイデアは、結果として生じる力ベクトルを 3 つの構成要素である相固有の力に分解することに基づいています (Tofighi and Yildiz (2013) を参照)。これらの相固有の力は、1 種類の界面のみが可能な 2 相システムの表面力と同じ方法で個別に処理されます。この目的のために、相 n と βの間の界面張力は、人工的に導入された相固有の表面張力によって表現され、 σnβ = σn + σβ
ここで、:
σ1 = 0.5(σ12 + σ13 - σ23)
σ2 = 0.5(σ12 + σ23 - σ13) (4)
σ3 = 0.5(σ13 + σ23 - σ12)
三相システムの問題点の 1 つは、すべての相が直接接触する可能性があることです。しかしながら、これらの状況は、前述の表面張力を相固有の表面張力の合計に分解することによって自動的に反映されます。
単一液体中を上昇する気泡の問題（二相系）
考慮される流れは、Hysing et al (2009) のテストケース 2 として定義されている、液柱内を上昇する2次元の気泡です。このケースは、密度が密度よりもはるかに低い気泡に関するものであるため、産業用途を代表するものです。解析領域を図 1 に示します。液体 (第1相) の特性は、ρ1=1000 kg/m3 および μ1=10 Ns/m2 と見なされます。ガス (第2相) の特性は、ρ2=1 kg/m3 および μ2=0.1 Ns/m2 に設定されます。表面張力と重力加速度はそれぞれ σ=1.96 N/m と g=0.98 m/s2 に設定されます。
図１ 二次元気泡のテストケース解析モデル図
この問題では、気泡の重心の位置、気泡の上昇速度、気泡の真円度/球形度、気泡面積、および表面周囲長を予測することです。これらのパラメータについて PHOENICSの結果と、他で行った結果と比較します。円形度は、二次元の面積相当円の外周長Paに対する気泡表面の外周長Pbの比です。これは計算の開始時に 1 の値を示し、バブルが変形するにつれて減少します。
PHOENICS VOF シミュレーションは、h =1/40、1/80、1/160、および1/320 で定義される各座標方向の均一なメッシュ サイズ h を使用して、シミュレーション時間を3 秒間として実行しました。比較のために、カラー関数の輸送方程式における非線形対流項の離散化について、CICSAM および THINC 界面解像度スキームを使用して計算しました。
結果を図 2 ～ 5 に示します。まず、図 2 は、グリッド解像度 320 でのバブルの時間変化のスナップショットを示しています。図 3 は、時間 t = 3 秒でグリッド解像度 320 の PHOENICS で得られた気泡の形状を示しています。Gamet et al (2018) が Interfoam と InterIsoFoam を使用して得た数値結果も併せて示します。図 4 は、PHOENICS VOF-CICSAM で得られた結果ですが、この図は、グリッド解像度 40、80、160、および 320 での質量中心、上昇速度、真円度を示しています。図 5 は、解像度 320 での重心、上昇速度、真円度について、PHOENICS VOF-CICSAM の結果と PHOENICS VOF-THINC で得られた結果を比較したものです。
図２ PHOENICS VOF-THINC 計算結果
図３ PHOENICS(1/h=1/320),InterFoam(1/h=1/160),InterisoFoam(1/h=160)のt=3sの気泡形状
図４ PHOENICS VOF-CICSAMの時間変化－メッシュサイズの比較 (a) 質量中心、(b)上昇速度、(c)真円度、(d)閉じた気泡の真円度
図５ PHOENICS VOF-CICSAMとTHINCの時間変化の比較 (a) 質量中心、(b)上昇速度、(c)真円度、(d)閉じた気泡の真円度
VOF 方式を使用においてCICSAM,THINCのどちらでも良好な動作が得られました。
さらに、この結果は他の公表された結果とよく比較されますが、より深い研究はまだ残っています。今後、スムージングレベル、ディラック関数のカットオフ、掃引回数など、PHOENICS VOF メソッドのいくつかのパラメータを変化させた場合の影響を調査する予定です。
2相液柱内の上昇気泡による3相系への応用
このセクションでは、 層状になっている密度の異なる２つの液体を通って上昇する気泡について考えます。表面張力、粘度、重力が考慮されます。
計算領域は正方形の形状 Ω = [0.0, 1.0] × [0.0, 1.0] をとり、時間ゼロでは、バブル (フェーズ 1) は (((?? ?0.5)/??)^2 + ((?? ? 0.325)/??)^2= 1で定義される楕円形状を持ちます。（長軸と短軸はそれぞれ a=0.15 と b=0.075）第２相と第３相はy=0.5の高さで上下に分離し、密度はそれぞれ 1000 kg/m3 と 1500 kg/m3 です。重力は、鉛直方向の負向きに 9.8 m/s2 の大きさで作用します。計算は３パターン行い、表 1 に示すように、3つの相の物性を指定してします。メッシュ分割は、200 × 200 セル にします。
表１ ケース番号と物性値
図６ ３相流気泡上昇問題における密度コンター図
4 つの異なる時間 (t = 0.02、0.12、0.28、0.5 秒) での数値結果を図 6 に示します。
気泡は浮力によって上昇し、上方に移動するにつれてその形状が崩れる傾向があることがわかります。気泡は浮力によって上昇し、上に行くほど形状が崩れやすいことがわかります。これは、表面力が円形を維持できるほど大きくないためです。粘度の影響により、気泡の上昇速度とその形状が減衰します。
ケース 1、2、3 のレイノルズ数は、それぞれ 8267.43、82.67、 82.67 です。レイノルズ数が大きいほど、気泡の歪みは大きくなります。2 番目のケースでは、他のケースよりも高い動粘度を使用するため、レイノルズ数が小さくなります。図からケース 2 では気泡の歪みが少ないことを示しています。ケース 1 (左列) とケース 3 (右列) の結果を比較すると、表面張力の影響がわかります。ケース 3 の結果は、ケース 1 に比べて表面張力が大きくなったため、気泡が分裂しないことがわかります。ケース 3 の場合、気泡が最終的にフェーズ 2 の流体に完全に浸漬されることを示しています。
結論
PHOENICS VOF オプションは、3 つの異なる混ざりあわない流体の流れをシミュレートするように拡張されました。この実装は、層状に分離された液体内を気泡が上昇する問題で検証されました。各流体の物性をパラメーターにして、その変化に応じて気泡の動きと歪みが定性的に確認できました。今後の研究は、公表されている数値結果との比較による三相システムの定量的な検証を目的としています。
|
