自由表面を伴う流れ
流体の流れの問題の多くは、複雑な幾何形状において自由表面を伴うもので、その大半は極めて非定常な現象でもあります。流動現象の解析対象としては、洪水吐や河川での流れ、橋杭の周りの流れ、洪水などによる越流、さらには水路や水門をはじめとする各種の構造物を伴う環境での流れなどがあります。こうした流れの現象を数値解析用にモデル化できる機能は魅力的ですが、ある程度の計算リソースで正確に計算できることが前提になります。また、有用なシミュレーションを行うには、物理モデルを使用した場合よりも、短時間かつ低コストで実現できることが不可欠です。
コンピュータプログラムの多くは、流体のさまざまな現象を記述した偏微分方程式を解くことができます。しかし、自由表面を伴うシミュレーションを実行できるものは、それほど多くありません。この原因は、自由境界問題と呼ばれるよく知られた数学的な問題です。自由境界の取り扱いが難しいのは、表面が移動するに従って計算領域が変化する一方で、その表面移動自体が計算によって決定される点にあります。計算領域の変化は、そのサイズや形状の変動だけでなく、場合によっては領域の合体や分割(すなわち、自由表面の発生や消滅)を伴います。
本書では、任意の自由表面を伴う流れの現象を数値解析用にモデル化するための手法について説明します。この手法はVOF(Volume-of-Fluid)法に基づくもので、特に自由表面流れに適した、さまざまな特長を備えています。本書では、VOF法が自由表面とその発生や消滅を捉えるのに最も自然かつ極めて効率的な手法である根拠を示します。
VOF法の特長をよく示すために、シンプルですが非常に重要な、流動現象に関する問題を取り上げます。ここでは、階段工での落下流を例にとります。概念的に単純な流れであるとともに、結果の妥当性を確認するための良好な実験データも提供されています(N. Rajaratnam and M.R. Chamani、『Energy Loss at Drops』、J. Hydraulic Res. Vol. 33、p.373、1995を参照)。
自由表面を伴う流動現象のプロトタイプ
図1aは、定常状態に達した後の流れの問題を示しています。階段工上部からの越流(液体またはナップのシート)には、上下両方の自由表面があります。越流の下側には、越流と階段工の面の間にプールが形成されており、下流では、液体は平坦な定常表面で右に向かって流れています。厳密に言うと、プール領域の流れの状態は定常です。これは、衝突する液体によって、プールに乱流混合が発生しているからです。ただし、平均的な構成が存在し、それが実験で報告されます。
実用目的では、流れは常に2次元です。つまり、図1aの図に垂直な方向では大きな変動はありません。現実に、プールの上方に空間を作るには、大気への何らかの開口部が必要で、それがなければ閉じてしまいます。
階段工の上部の流速は重要です。つまり、これは表面波と同じかそれ以上の速度であるため、下流からの擾乱がこの領域を貫通して上流の流れ(階段工の左側)に影響することはできません。そのため、この領域での流れは例外的に滑らかであり、定常です。
この問題には、数値シミュレーションと比較できる幾何形状機能が多数あります。たとえば、階段工の前後の流れの高さ、越流が底部に衝突するときの角度、越流の下に形成されるプールの深さなどです。さらに、実用化のための重要な比較として、階段工を通過する落下流によって失われるエネルギーの量(運動エネルギーと位置エネルギーの合計)があります。
プロトタイプ問題のシミュレーション
図1aはシミュレーションの結果です。この例では、実験で使用したすべての幾何形状および物質の特性がシミュレーションで使用されました。実験室のテストで使用した階段工の高さは62cmで、液体は普通の水(密度=1.0gm/cc、動粘性=0.01dynes/cm)です。計算領域に入る水の深さは15.5cmで、速度は臨界に近い123.0cm/秒でした。もちろん、重力は垂直方向で、大きさはg=-980cm/s2です。
図1a. 階段工での落下流のシミュレーション(左)、1b. シミュレーションで使用した格子(右)
越流の左側にあるプールに乱流が発生することが予測されたため、シミュレーションでは乱流モデル(繰り込み群、つまりRNGモデル)を使用しました。その後、乱流モデルを使用せずに行ったシミュレーションでもよく似た結果が得られましたが、これはそれほど意外なことではありません。流れの重要な要素のほとんどは、滑らかな(つまり、乱流ではない)流入、流出、越流だからです。
図1bのシミュレーション領域は、幅170cm、高さ100cmで、横80個、縦60個、合計4800個のセルで構成される、同じサイズの長方形セルの格子に細分されています。この格子は、流体力学の支配微分方程式(ナビエストークス方程式)の有限差分近似の基底として使用されます。格子セルの数とサイズは、流れの中で予測される最小の特性を捉える目的で選択されました。結果を見て何らかの調整が必要と思われる場合は、数を簡単に増減できます。実際に、解像度を変えてシミュレーションを繰り返すことにより、計算がそうした変化に影響されすぎないことを確認することをお勧めします。
左側の境界は、指定した速度境界です(流体の高さも指定)。右側の境界は流出境界で、すべての流量が境界に対して垂直なゼロ勾配であり、均一な流出が促進されます。上下の境界は剛壁で、3番目の方向では、境界は対称の面(粘性抵抗ゼロの壁)として処理されました。階段工の表面も、自由すべり境界として処理されました。
初期条件は、予測される流れの配列を粗近似するように設定することもできましたが、流れの構成は計算したいものの1つであるため、流体がどのように分布されるかがわからない場合は特に、よりシンプルなアプローチが必要です。この例では非定常の流れシミュレータを使用したため、図1aの階段工の上に流体のブロックだけがあり、左側の境界に同じ水平速度と高さが割り当てられている、シンプルな初期条件を定義できました。シミュレーションはその後、定常流れに発展しますが、これはおよそ8.0秒後に発生します。シミュレーションは、定常状態に達したことを確実にするために、10.0秒の時間まで実行されました。図2は、中間時間を2つ示しています。図2bは0.2秒、図2cは0.5秒の時点で、図2dは最後の10.0秒の時点を示しています。
図2a～2d (左から順に)0.0秒、0.2秒、0.5秒、10.0秒経過時のシミュレーション
最初は単一の結合している自由表面だったものが、流体が底部に衝突した後で2つの独立した自由表面(上下のナップ表面)に変化することに注目してください。下境界の衝撃点の左右に流れが分離しても、問題はありません。これについては、次のセクションで詳細に説明します。
実験とシミュレーションの比較は、以下の表のとおりで、極めてよく一致しています。
これらの結果を考慮すると、このような精度を達成するには、かなりの計算時間が必要になると思われるかもしれません。しかし実際は、Pentium 4、3.20GHzのデスクトップコンピュータでの合計CPU時間は、たったの88秒でした。計算時間がこれほど短くなることには説明が必要です。それが、これ以降のセクションの目的です。
VOF法が適している理由
VOF法の仕組みと、それが非常に効率的な手法である理由を理解するために、さまざまな計算法、その中でも特にVOF法に関するいくつかの基本概念を示します。
基本理論
どのような数値法でも、流れの問題を単純な算術演算で計算できるように有限の数値セットに簡略化する必要があります。連続流体を離散化された数値セットに近似するために一般的に用いられるのが、流体が占める空間を格子に分割する方法です。この格子は、通常は、多数の小さな長方形のブロック(要素)で構成されます。これらの各要素に対して平均化処理を施すことで、その要素における流体の圧力、密度、速度、および温度の代表値を得ることができます。
単純な式を使い、ある時間にわたる各要素の値と隣接要素との相互作用を近似することができます。たとえば、ある要素の密度は、その要素と隣接要素との間で(質量保存による)正味の質量流量が交換されたときのみ変化します。要素間で質量が交換される物質の速度は運動量保存則によって計算され、通常はナビエストークス方程式で表されます。ナビエストークス方程式は、隣接する要素間で作用している圧力と粘性応力を用いて、要素内で変化する流体速度を近似します。
こうした要素と隣接要素の間での相互作用に基づく考え方は、偏微分方程式で近傍の量の変化によって生じる小さな変動の効果を評価することと本質的には同じです。工学系のテキストなどでは、小さなコントロールボリュームを使用して、そのサイズを無限小まで小さくした近似の極限として偏微分方程式が導かれます。数値シミュレーションでも同じ方法がとられますが、要素数が多すぎると追跡が困難になるため、コントロールボリュームのサイズを極限まで小さくすることはできません。実際のシミュレーションでは、現象を解くのに十分かつ計算時間を最小限まで抑えることができる要素数を設定することが目標となります。
要素に対して用いられる算術演算は、基本的には加算、減算、乗算、および除算のみを伴う単純なものです。たとえば、ある要素での質量の変化は、一定の時間間隔にわたって要素の面から流入および流出した質量の加算や減算で求めることが可能です。ただし、シミュレーションでは、こうした演算を数千、ときには数百万もの要素に対して、非常に短い時間間隔について繰り返して計算する必要があります。そのため、このような反復計算の高速処理にはコンピュータが適しています。
自由表面を伴う流体運動のシミュレーションでは、形状が変化する計算領域を扱わなければなりません。この複雑さに対応できる解析法が、次に説明するVOF法です。
VOF法の概念
VOF法は、各格子セルに、セル体積のうち液体が占める割合、すなわち体積占有率を記録するという考え方に基づきます。通常、この体積占有率は量Fで表されます。Fは体積占有率であるため、値がとりうる範囲は0.0～1.0です。
液体内部の領域では、Fの値は1.0になり、液体の外部、つまり(空気などの)気体領域では、Fの値は0になります。Fの値が0.0と1.0の間で変化する場所が、自由表面が存在する位置です。すなわち、0.0より大きく、かつ1.0未満のF値を持つ要素は、必ず表面を有します。
ここで留意すべきことは、VOF法では自由表面を直接的に定義するのではなく、バルク流体の位置を定義するという点です。これにより、計算上の難しさをもたらすことなく流体領域を合体または分割できます。自由表面は、単に、流体の体積占有率が1.0と0.0の間で変化する場所であると定義されます。これは、自由表面を伴うほぼすべての問題に適用できるVOF法の非常に優れた特長でもあります。
また、格子の各要素に単一の数値(F)を割り当てることで流体の位置を記録できる点も、VOF法の重要な特長です。これは、圧力成分や速度成分など、他のすべての流体物性の平均値を要素内に記録できる機能とも一貫しています。
VOF法の詳細
図3 1D要素列内の表面
図4 2D要素格子内の表面
1次元のケースのように、まず、近隣要素をテストすることによって、表面のおおよその方向を見つける必要があります。図4では、外向きの法線が上向きの方向に最も近くなります。これはその方向の近隣の値の差が、他の方向より大きいからです。次に、ほぼ垂直方向にある要素列で、表面の局所的な高さが計算されます。図4の2次元のケースでは、これらの高さを矢印で示しています。最後に、表面要素を含む列内の高さにより、その要素内の表面の場所が特定されます。他の2つの高さを使用すると、局所的な表面傾斜と表面曲率を計算できます。
3次元でも同じ手順を使用しますが、表面要素の周囲にある9個の列に関して列の高さを求める必要があります。必要な計算は少し多くなりますが、主な内容は、列内の単純な足し算と、傾斜や曲率を求めるための列の高さの加減算です。この説明に基づいて、流体の体積占有率を使用して自由表面の定義に必要なすべての情報を短時間で簡単に求める方法が、読者の皆さんにもお分かりいただけたでしょう。
扱わなければならない問題があと2つ残っています。1つは、図1および2のようなシミュレーションが、流体が存在する領域の流体力学のみに対応していることです。これは、VOF法の計算効率が高いことの、もう1つの理由です。階段工での落下流の問題で流体が占有する領域は、計算格子のオープン領域の半分以下です。液体を取り囲む気体の流れも計算する必要があるとしたら、必要な計算時間は大幅に増えます。ただし、液体のみで計算を行うには、自由表面で境界条件を指定する必要があります。この条件とは、接線応力の消滅と、気体の圧力に等しい標準圧力を表面に加えることです。
2番目の問題は、自由表面が流体と共に動くときの動きと変形を、流体占有率変数Fを求めることによって計算しなければならないことです。変数Fは不連続(主に0.0または1.0)であるため、計算格子を移動するときにこの不連続性が維持されるように注意する必要があります。VOF法では、この目的で特殊な移流アルゴリズムが使用されています。
VOF法による自由表面の追跡の図解
図5aは、これがいかに適しているかを示しています。流体の体積占有率は、格子要素ごとに均一に色分けされ、その要素における値を表しています。自由表面は、ほぼすべての場所で鮮明に定義されています。ナップの一番低い、一番狭い部分にのみ、鮮明な流体分布の損失が確認できます(図5b)。これは予測どおりです。この領域では、ナップの厚さは3つの要素より小さく、これにより、部分充填された表面要素に関連付けられている小さいF値がいくつか中心要素(値1.0)に混入するからです。計算目的では、このことはあまり問題になりません。このシミュレーション法では、液体内部の要素は純粋な液体要素と同様に処理されるからです。
図5bに示す領域では、実際の実験で乱流と空気混入が観察されたということも、指摘しておかなければなりません。したがって、流体占有率の値を1より少し小さく見せることがやや現実的です。このことは、全く予想外というわけではありません。液体の噴流がプールと交わることは、乱流や空気混入の原因となりますが、流体占有率の値を液体内部に「取り込む」原因ともなるからです。
図5a 要素内の流体占有率の値、表面定義の鮮明さを示す(左)、5b 越流が底部に衝突する地点の流体占有率値の拡大図(右)
まとめ
最初は、コンピュータで数値の配列の算術演算を繰り返し行うだけで、複雑な、時間に依存する流体力学問題の現実的なシミュレーションを生成できるなんて、なんだか魔法のようだと思われたかもしれません。これを比較的初歩的な手順で行うアプローチについて説明することが、ここでの目的でした。
ここでは、シンプルながらも重要な流動現象の例を使用して、コンピュータのシミュレーションによって、物理的測定と極めてよく一致する詳細な結果が得られるということを示してきました。また、VOF (Volume of Fluid)法に基づくシミュレーションでは、正確かつ効率的な単純近似法を使用することも説明してきました。
確かに、実世界の例では、水力発電所で使用されるような複雑な水力施設が関わってくるため、役に立つ結果を得るには、ここでの例で使用した数秒という計算時間よりも長い時間がかかることは間違いありません。それでも、その結果は(人にとってもコンピュータにとっても)適正な時間で生成でき、物理的実験ではほとんど不可能な詳細を豊富に含めることができます。例として、弊社の水への応用に関するページを参照してください。さらに、幾何形状、流れの状態、流体特性のどんな種類の変化による影響でも簡単にテストできるということは、シミュレーションを採用するもう1つの強力な理由です。流動シミュレーション用の最新のソフトウェアおよびハードウェアは、従来の物理的モデル化に比べてコスト面で大きなメリットがあります。
追記
VOF法は、1981年、C.W. Hirt氏とB.D. Nichols氏によって、J. Comp. Phys., 39, p.201で初めて詳細に説明されました。本書に出てくるすべてのシミュレーションは、米国Flow Science, Inc. (Santa Fe, NM, 87505)が開発した市販のソフトウェアパッケージFLOW-3Dを使用して実行しました。このプログラムでは、VOFの概念の強化版であるTruVOFを使用しています。