自由表面モデル化の方法
気体と液体の間の界面は、多くの場合、自由表面と呼ばれます。「自由」という呼び方になったのは、気体と液体の密度が大きく異なるからです(たとえば、水の空気に対する密度比は1000です)。気体の密度が低いということは、液体の慣性と比べて、気体の慣性は通常は無視できるということを意味します。この意味で、液体は、気体に対して独立して、つまり自由に動きます。気体の唯一の影響は、液体の表面にかかる圧力です。つまり、気体と液体の表面は、制約されているのではなく、自由であるということです。
熱伝達に関する文書では、自由境界の問題を表す際に、スティーブンの問題という用語がよく使用されます。ただし、ここでは、境界は相境界、つまり、対流的な流体の流れによって供給された熱に反応して変化する氷と水の間の境界などを指します。
名前が何であろうと、自由または可動の境界が存在すると、どんなタイプの解析でもかなり複雑な問題になることは明らかです。最も簡単な問題を除けば、数値解法という手段に出る必要があります。その場合でも、自由平面では、特殊な手法を導入して、場所や動き、流れへの影響を定義する必要があります。
以下の説明では、自由表面のモデル化に使用されてきたさまざまなタイプの数値的アプローチについて簡単に見直して、各手法の長所と短所を説明します。どの手法を採用するかにかかわらず、自由表面を適切にモデル化するには、以下の3つの機能が必須です。
- 表面の形状と場所を説明するスキーム
- 時間の経過と共に形状と場所を発展させるアルゴリズム
- 自由表面境界条件を表面に適用すること
ラグランジュ格子法
概念上、自由表面を定義および追跡する最も簡単な方法は、流体に埋め込まれて共に移動するラグランジュ格子を構築することです。多くの有限要素法で、このアプローチが使用されています。格子と流体は共に動くため、格子は自動的に自由表面を追跡します。
表面では、適切な境界条件を含め、境界の片側のみに流体が存在することを説明するように、近似方程式を修正する必要があります。これを行わないと、非対称が進展し、最終的にはシミュレーションの精度が失われます。
ラグランジュ法の第一の制限事項は、分散している表面や交差している表面を追跡できないことです。表面の動きの振幅が大きい場合であっても、ALE (Arbitrary-Lagrangian-Eulerian)法のような格子再生成手法を導入せずに追跡するのは困難です。これらのアプローチの初期の例については、参考文献の1970および1974を参照してください。
ここで説明する残りの自由表面手法では、より複雑な表面の動きを処理できるように、固定のオイラー格子を計算の基礎として使用します。
表面高さ法
振幅の小さいスロッシングや浅水波、およびその他の、表面が水平からそれほど逸脱していない自由表面の動きは、何らかの基準高度に対して相対的な表面の高さHで表すことができます。高さの時間発展は、(u,v,w)を(x,y,z)方向の流体速度とする運動学方程式によって支配されます。この方程式は、表面が流体と共に移動しなければならないという事実を数学的に表現したものです。
この方程式の有限差分近似は、簡単に実行できます。さらに、3次元数値解法のメモリ要件が極度に小さくなるように、同じ高さの場所の値のみを記録する必要があります。最後に、自由表面の境界条件の適用も、ほぼ水平なままという表面の条件によって簡潔化されます。この手法の例については、参考文献の1971および1975を参照してください。
MAC (Marker-and-Cell)法
時間依存性を持つ自由表面の流れの問題に関して最初に考案された数値法が、MAC (Marker-and-Cell)法です(参考文献1965を参照)。このスキームは、コントロールボリュームの固定のオイラー格子に基づいています。格子内の流体の場所は、流体と共に動き、それ以外は体積、質量、その他の特性を持たない、一連のマーカ粒子によって決定されます。
マーカを含む格子セルは流体で占有されていると見なされ、マーカがない格子セルは空(ボイド)です。粒子を含み、少なくとも1つの隣接格子セルがボイドである格子に、自由表面は存在すると定義されます。セル内の表面の場所と向きは、当初のMAC法には含まれていませんでした。
表面の発展は、局所的に補間された流体速度でマーカを動かすことによって計算されました。新しく充填された格子セルの流体特性を定義したり、空になったセルの値をキャンセルしたりするには、特殊な処理が必要でした。
自由表面の境界条件を適用するには、すべての表面セルに気体の圧力を割り当てる必要がありました。また、速度成分は、非圧縮性とゼロ表面せん断応力の条件を近似する方法で、表面上、または表面のすぐ外側のすべての場所に割り当てられました。
幅広く複雑な自由表面の流れの問題を解く上でMAC法が素晴らしい成功を収めていることは、数々の文献で十分に立証されています。この成功の理由の1つは、マーカが表面を直接追跡するのではなく、流体の体積を追跡することです。表面は体積の境界にすぎず、そういう意味では、表面は、分割または合体された体積として出現、マージ、消滅する可能性があります。
さまざまな改善により、当初のMAC法の精度と適用性は向上してきました。たとえば、セル内の補間された表面の各所に気体の圧力を加えると、静水圧力の影響を受ける問題の精度が向上し、一方で、表面張力を含めると、MAC法の適用範囲がより幅広い問題に拡大されます(参考文献1969、1975を参照)。
数々の成功にもかかわらず、MAC法は主に2次元のシミュレーションに使用されてきました。これは必要な数のマーカ粒子に対応するには、かなりのメモリとCPU時間が必要だからです。通常は、大きな変形を伴う表面を正確に追跡するためには、格子セルごとに平均で約16個のマーカが必要です。
マーカ粒子のもう1つの制限は、流れが収束/発散する領域で流動過程を追跡する能力が、あまり十分ではないことです。マーカは通常、小さい流体要素の中心を追跡するものと解釈されます。しかし、その流体要素が細長い渦巻状になると、マーカは流体の構成を十分に示すことができなくなる場合があります。この現象は、たとえば、流れのよどみ部分で、マーカが一方向に集積するものの、垂直方向に引き離される場合などに確認できます。十分に引き離された場合(格子セル1個の幅以上)は、流れの中に非物理的ボイドが発生することがあります。
表面マーカ法
マーカのメモリおよびCPU時間の消費を制限する方法の1つは、マーカ粒子を、流体領域の内部ではなく、表面のみに保持することです。もちろん、これによってMAC法の体積追跡特性が排除されるため、表面が分割または合体するタイミングや方法を特定するための論理を追加する必要があります。
2次元の場合は、表面上のマーカ粒子は表面に沿って線形に配置できます。この配列には、粒子の間隔を均一に保つことができる、別々の表面が交わる部分の計算が簡単になる、など、いくつかのメリットがあります。また、表面マーカを使用すると、境界条件を適用する際に、格子セル内で表面を簡単な方法で見つけることができます。
あいにく、3次元では、表面上に粒子を並べる簡単な方法はありません。このことが、表面マーカ手法の大きなデメリットとなっています。表面が拡大していて、空間がマーカで充填されていない領域が存在する場合があります。マーカがないと、表面の構成は不明なため、結果として、マーカを追加する方法がなくなります。
参考文献1975には、この方法のメリットと制限を示す例が含まれています。
VOF (Volume-of-Fluid)法
最後に説明する方法は、流体の体積占有率の概念に基づいています。このアプローチは当初、MAC法の強力な体積追跡機能を、メモリやCPUのコストを抑えて使用するための方法として考えられました。
各格子セル(コントロールボリューム)内では、流量ごとに1つの値(圧力、速度、温度など)だけを保持することが慣習となっています。この理由により、自由表面を見つけるためにこれより多くの情報を保持しても、ほとんど意味がありません。この推論に従うと、単一の量、各格子セルの流体体積占有率を使用することは、他の流量を解決することと整合性がとれています。
各セル内の流体の量がわかっている場合は、表面を見つけることも、表面傾斜や表面曲率を特定することも可能です。表面は、流体が部分充填されているセル内、または流体が全体に充填されているセルと流体が全くないセルの間に存在するため、簡単に見つけることができます。
傾斜と曲率は、隣接セルの流体体積占有率を使用して計算されます。体積占有率は階段関数でなければならないこと、つまり、値が1か0であることを覚えておくことは重要です。このことを知っていれば、隣接セルの体積占有率を使用して、特定のセル内の流体の位置(およびその傾斜と曲率)を見つけることができます。
自由表面の境界条件を、MAC法と同様に適用する必要があります。つまり、適切な気体圧力(および対応する表面張力)を割り当て、また、表面におけるゼロせん断応力条件を満たすには表面外側のどの速度成分を使用する必要があるかを特定します。実際には、表面における速度成分の代わりに速度勾配を割り当てる方が簡単な場合があります。
最後に、表面の時間発展を計算するには、分布の階段関数の性質が保持されるような方法で、格子を通過して体積占有率を移動する手法が必要です。流体占有率の基本的な運動学方程式は、高さ関数法と同様です。Fは流体占有率関数です。
この方程式をモデル化する際は、簡単な数値近似は使用できません。数値の拡散や分散エラーによって、F分布の明確な階段関数の性質が損なわれるからです。
F分布が0または1の値を保持するような1次元では、この方程式の解を正確にモデル化することは簡単です。1列のセルに上から下まで流体が充填されている場合を想像してください。ある瞬間に、流体の界面は、あるセルの中間領域にあり、その下側の隣接セルは充填されていて、上側の隣接セルは空になっています。隣接セル内の流体の向きは、界面とセルの下端との距離が、セル内の流体占有率と等しくなければならないことを意味します。次に、まず表面を含むセルの空の領域を充填してから次のセルに流体を送るように、上側の空のセルに移動する流体の量の計算を変更できます。
2次元や3次元でも、隣接セルからの情報を使用する同様の手順を使用できますが、1次元の場合ほど正確にすることは不可能です。2次元以上の場合の問題は、表面の形状と場所を正確に特定できないことです。それでも、VOF法を使用して達成された多数の成功例からもわかるように、この手法をうまく機能させることはできます。この手法に関する当初の研究については、参考文献の1975、1980、1981を参照してください。
VOF法は、MAC法と同じくらい強力な手法をオーバーヘッドなしで提供するという目標を達成してきました。表面追跡ではなく体積追跡機能を使用することは、流体質量の分割と合体を処理するのに十分な堅牢性を持っていることを意味します。さらに、連続関数を使用するため、離散化された粒子で発生する、数値が割り切れないという問題に見舞われることもありません。
VOF法の可変密度近似
VOF法の特殊な処理が必要な機能の1つは、境界条件の適用です。表面が格子を通過して移動するとき、流体を含むセルは絶えず変化します。つまり、計算領域も変化しているということです。この変化している領域の自由境界には、適切な自由表面応力条件も適用する必要があります。
流体領域の更新や境界条件の適用は、重要な作業です。このため、液体と気体の両方の領域で流れが計算される、VOF法への近似が使用されてきました。通常、これは、可変密度を持つ単一の流体として流れを処理することによって行われます。密度を定義するには、F関数を使用します。それから、流れ方程式は液体と気体の両方の領域で計算されるため、界面の境界条件を設定する必要はないという論証が行われます。
あいにく、このアプローチは、2つの理由によって、実際にはあまりうまく機能しません。1つは、圧力の変化に対する気体領域の感度が、通常は液体領域より大幅に大きいことです。このため、圧力-速度結合解法で収束を達成することは困難です。この手法では、必要なCPU時間が非常に大きくなる場合があります。
2番目の、より重要な理由は、界面において接線速度が不連続になる可能性に関連しています。圧力に対する反応がそれぞれ異なるため、界面における気体と液体の速度は、通常は大きく異なります。可変密度モデルでは、界面は平均速度で動きますが、これによって界面の動きが非現実的になる場合が多くあります。
可変密度法は、VOF法と呼ばれることもありますが、流体占有率関数を使用するため、この呼び方は正しくありません。液体と気体の明確な界面を正確に追跡するためには、実際は界面を不連続として処理する必要があります。これはつまり、界面の不連続性を定義する手法と、その界面に適切な境界条件を課す方法が必要ということです。また、特殊な数値法を使用して、不連続という特徴を損なうことなく、格子を通過する界面の動きを追跡することも必要です。
まとめ
ここでは、自由表面を数値的にモデル化する際に使用するさまざまな手法について、相対的な長所と短所に関する解説も含めて、簡単に説明してきました。長年に渡って、こうした基本的な手法が数多く提案されてきたことを知っても、読者の皆さんは驚かないでしょう。おそらく、最も成果を上げた方法は、簡潔かつ堅牢なVOF法でしょう。この方法に、一部改良を加えたものが、FLOW-3Dプログラムで使用されています。
VOF法の改良では、より良い、より正確な方法で流体占有率を格子を通過して移動することに重点が置かれてきました。その他の開発では、物体適合格子に関連して手法を適用したり、複数の流体成分をモデル化するために複数の流体占有率関数を採用したりしてきました。こうした開発についての議論は、ここでの説明の範囲外です。
参考文献
1965 F.H. Harlow、J.E. Welch、『Numerical Calculation of Time-Dependent Viscous Incompressible Flow』、Phys. Fluids 8、2182
1969 B.J. Daly、『Numerical Study of the Effect of Surface Tension on Interface Instability』、Phys. Fluids 12、1340
1970 C.W. Hirt、J.L. Cook、T.D. Butler、『A Lagrangian Method for Calculating the Dynamics of an Incompressible Fluid with Free Surface』、J. Comp. Phys. 5、103
1971 B.D. Nichols、C.W. Hirt、『Calculating Three-Dimensional Free Surface Flows in the Vicinity of Submerged and Exposed Structures』、J. Comp. Phys. 12、234
1974 C.W. Hirt、A.A. Amsden、J.L. Cook、『An Arbitrary Lagrangian-Eulerian Computing Method for all Flow Speeds』、J. Comp. Phys.、14、227
1975 B.D. Nichols、C.W. Hirt、『Methods for Calculating Multidimensional, Transient Free Surface Flows Past Bodies』、Proc. of the First International Conf. On Num. Ship Hydrodynamics、Gaithersburg、ML、Oct. 20-23
1980 B.D. Nichols、C.W. Hirt、『Numerical Simulation of BWR Vent-Clearing Hydrodynamics』、Nucl. Sci. Eng. 73、196
1981 C.W. Hirt、B.D. Nichols、『Volume of Fluid (VOF) Method for the Dynamics of Free Boundaries』、J. Comp. Phys. 39、201