【０００１】
【発明の属する技術分野】
本発明は、流速が音速に比べて十分に遅い非圧縮流れの場において、気体と液体との二相からなりかつ気液自由界面を有する一次元もしくは多次元系の速度分布、圧力分布、密度分布および気液自由界面形状を求める気液自由界面シミュレーション装置および方法に関する。
【０００２】
また、本発明は、液滴や気泡の運動、変形、分裂、結合、固体表面の表層流れ、トラップ・配管内流れおよび波のシミュレーションにおいて必要となる気液自由界面の界面計算を高速、高精度かつ安価に実現する演算装置および方法に関する。
【０００３】
【従来の技術】
従来、気液自由界面流体シミュレーションでは、最も基本的な手法として、例えば各格子の流体率を気体を０、液体を１としてあらわし、流体率分布を質量保存式から得られる非保存型移流方程式で移流させる方法がある。本手法では、移流方程式の解法になるべく精度の高い方法を用いることで、気液自由界面シミュレーションを試みる。
【０００４】
しかしながら前記の従来手法では、移流方程式解法の精度に限界があるため、実際的には計算サイクルが増えると気液界面に０と１の中間的な値を持つ数値拡散層が大幅に増大し、気液界面がほやける課題を有し、且つ、気液界面形状から計算される表面張力の計算において正しい曲率計算ができなくなるという問題点を有する。
【０００５】
一方、最も一般的な気液自由界面シミュレーション手法であるＶＯＦ法では、流体率の移流計算後、ドナーアクセプタ法などで気液界面の探索や再構築を行う。本手法では、数値拡散は補正され、質量や体積の保存性も保証される（例えば、非特許文献１，２参照）。
【０００６】
しかしながら、前記の従来手法においては、気液界面の探索や再構築に複雑な条件判定が必要となり、計算プログラムが複雑化する。また気液界面の再構築アルゴリズムは物理的現象に基づくものとはいえず、計算プログラムの構成方法次第で計算結果が大きく変わるため、計算精度においても課題を有する。
【０００７】
他方、ＣＩＰ法とよぶ高次精度移流計算解法およびデジタイザ法を組み合わせた気液自由界面シミュレーション方法がある。ＣＩＰ法では移流方程式の補間関数に三次式を用い、次の時刻ステップでの参照点における値をセミラグランジェ的に持ってくる。またデジタイザ法ではタンジェント関数を用いて流体率を変数変換する。本手法においては、計算プログラムは極めて簡易でありかつ数値拡散もほとんどない結果が得られる（例えば、非特許文献３参照）。
【０００８】
しかし、前記気液自由界面シミュレーション方法においては、質量や体積の保存性が保証されず、例えば微少の液滴や薄膜の液体が、計算サイクル数が多くなるに従い消失していくといった課題を有する。
【０００９】
このため、上記のような従来方法に対して、保存型ＣＩＰ法やＰＰＭ法などの保存型移流方程式解法を用いて移流方程式を解くことで、保存性を満足しつつＣＩＰ法やＰＰＭ法の持つ高い精度を利用する方法がある（例えば、非特許文献４，５参照）。
【００１０】
【非特許文献１】
Ｃ．Ｗ．Ｈｉｒｔ ａｎｄ Ｂ．Ｄ．Ｎｉｃｈｏｌｓ， Ｖｏｌｕｍｅ ｏｆ ｆｌｕｉｄ （ＶＯＦ） ｍｅｔｈｏｄ ｆｏｒ ｔｈｅ ｄｙｎａｍｉｃｓ ｏｆ ｆｒｅｅ ｂｏｕｎｄａｒｉｅｓ， Ｊ． Ｃｏｍｐ． Ｐｈｙｓ．， ３９， ２０１（１９８１）．
【非特許文献２】
Ｄ．Ｌ．Ｙｏｕｎｇｓ， Ｔｉｍｅ−ｄｅｐｅｎｄｅｎｔ ｍｕｌｔｉ−ｍａｔｅｒｉａｌ ｆｌｏｗ ｗｉｔｈ ｌａｒｇｅ ｆｌｕｉｄ ｄｉｓｔｏｒｔｉｏｎ， ｛￥ｅｍ Ｎｕｍｅｒ． Ｍｅｔｈｏｄｓ ｆｏｒ Ｆｌｕｉｄｓ Ｄｙｎａｍｉｃｓ｝， ｅｄｉｔｅｄ ｂｙ Ｋ．Ｗ．Ｍｏｒｔｏｎ ａｎｄ Ｍ．Ｊ．Ｂａｉｎｅｓ， ｐ２７３
（１９８２）．
【非特許文献３】
Ｔ．Ｙａｂｅ ａｎｄ Ｆ．Ｘｉａｏ， Ｄｅｓｃｒｉｐｔｉｏｎ ｏｆ Ｃｏｍｐｌｅｘ ａｎｄ Ｓｈａｒｐ Ｉｎｔｅｒｆａｃｅ ｄｕｒｉｎｇ Ｓｈｏｃｋ Ｗａｖｅ Ｉｎｔｅｒａｃｔｉｏｎ ｗｉｔｈ Ｌｉｑｕｉｄ Ｄｒｏｐ， Ｊ． Ｐｈｙｓ． Ｓｏｃ． Ｊａｐａｎ， ６２， ２５３７（１９９３）．
【非特許文献４】
Ｔ．Ｙａｂｅ， Ｒ．Ｔａｎａｋａ， Ｔ．Ｎａｋａｍｕｒａ ａｎｄ Ｆ．Ｘｉａｏ， Ｅｘａｃｔｌｙ ｃｏｎｓｅｒｖａｔｉｖｅ ｓｅｍｉ−Ｌａｇｒａｎｇｉａｎ ｓｃｈｅｍｅ （ＣＩＰ−ＣＳＬ） ｉｎ ｏｎｅ ｄｉｍｅｎｓｉｏｎ， Ｍｏｎ．Ｗｅａ．Ｒｅｖ．，１２９， ３３２−３４４．（２００１）
【非特許文献５】
Ｐ．Ｃｏｌｅｌｌａ ａｎｄ Ｐ．Ｒ．Ｗｏｏｄｗａｒｄ， Ｔｈｅ ｐｉｅｃｅｗｉｓｅ ｐａｒａｂｏｌｉｃ ｍｅｔｈｏｄ （ＰＰＭ） ｆｏｒ ｇａｓ−ｄｙｎａｍｉｃａｌ ｓｉｍｕｌａｔｉｏｎｓ， Ｊ． Ｃｏｍｐ． Ｐｈｙｓ．｝， ５４， １７４（１９８４）．
【００１１】
【発明が解決しようとする課題】
しかしながら、上記のような従来技術の構成をとった場合、以下のような課題がある。
ア）二次もしくは三次の補間関数をもつＣＩＰ法やＰＰＭ法では、数値拡散が防止できず、計算サイクル数が増大すると、界面がぼやける。
イ）四次の高次補間関数をＣＩＰ法に適用することで大幅に数値拡散を低減できるものの、数値拡散を０にすることはできず、計算サイクル数が非常に大きい場合、やはり界面がぼやける。
ウ）保存型移流方程式解法を適用する場合は移流方程式も保存型となるため、非保存型ＣＩＰ法で適用するデジタイザ法を用いることはできず、界面の数値拡散を補正することができない。
エ）ぼやけた界面に対し、例えば流体率が０．５の位置を気液界面とみなしても、界面の数値拡散は気相・液相両側に対し均等に数値拡散することは保証されていないため、気液界面で囲まれた液体もしくは気体体積は保存されない。
【００１２】
ここに本発明は、上記保存型移流方程式解法を用いた気液自由界面シミュレーション装置を改良し、保存性を満足しつつ数値拡散を防止する簡易かつ高精度なシミュレーション装置を提供することにある。
【００１３】
【課題を解決するための手段】
上記課題を解決するため本発明の請求項１は、流速が音速に比べて十分に遅い非圧縮流れの場において、気体と液体との二相からなりかつ気液自由界面を有する一次元もしくは多次元系の速度分布、圧力分布、密度分布および気液自由界面形状をナビエ・ストークス運動方程式および固定された直交もしくは非直交格子を用いてシミュレーションする気液自由界面シミュレーション装置において、気液自由界面からの距離を各計算格子において計算し、気液自由界面近傍であるか否かを判別するあらかじめ指定された距離値よりも大きい場合に、当計算格子は気液自由界面に含まれないと判別される気液自由界面外格子検出手段と、該計算格子における該距離値の符号により該計算格子が気相もしくは液相であると判別される気液判別手段と、移流方程式を計算することで得られた該格子の流体率もしくは密度が、判別された気相もしくは液相での正しい流体率値もしくは密度値と異なる場合、正しい値との差値を気液自由界面方向に移動させることで、移流方程式計算で生じた密度の数値拡散を補正する密度補正手段とを具えるようにした。
これにより、当計算格子は正しい値に修正され、また気液自由界面方向の隣接格子の値が該差値の分上昇するため、移流方程式計算で生じた密度の数値拡散に対し、あらかじめ指定した気液自由界面外の領域での数値拡散を除去することが可能となる。
【００１４】
また、本発明の請求項２では、各計算格子における気液自由界面からの距離の計算において、気相の流体率と液相の流体率の中間値を気液自由界面での値とし、移流方程式を計算することで得られた各計算格子の流体率から該中間値を引き、気液自由界面近傍であるか否かを判別する該距離値を乗算した値を距離計算の初期値として格納するようにした。
これにより、該中間値を表わす等値面上では気液自由界面からの距離が０、すなわち気液自由界面そのものとなるため、気液界面を動かさずに自由界面との距離を計算する初期値設定方法を構成することが可能となる。
【００１５】
また、本発明の請求項３では、多次元系において、気液自由界面近傍外と判別された領域における正しい流体率値もしくは密度値との差値を気液自由界面方向に移動させる補正において、気液自由界面からの距離分布に対し各軸方向の微係数をとることで得られる法線ベクトルを気液自由界面方向とみなすようにした。これにより、格子数を増やしていくと正しい気液界面方向に収束することが保証されるため、簡易な計算ながら理論的に誤りのない気液自由界面方向検出方法を構成することが可能となる。
【００１６】
また、本発明の請求項４では、各計算格子における気液自由界面からの距離の計算において、繰り返し計算を行い距離値を収束させる過程で、気液自由界面方向の近傍点を参照し、該近傍点での距離値を該近傍点を囲む全ての格子点の距離値より各軸方向一次関数を構成して補間計算し、該近傍点と該計算格子との距離を加えることで該計算格子での距離を計算するようにした。
これにより、格子間隔および該距離が十分小さい場合、該計算格子においてほぼ正確な距離となることを保証するため、簡易な各軸方向一次関数による補間計算を用いて、複雑な気液自由界面形状を崩さず気液自由界面からの距離分布を計算する方法を構成することが可能となる。
【００１７】
また、本発明の請求項５では、多次元系において、気液自由界面近傍外と判別された領域における正しい流体率値もしくは密度値との差値を気液自由界面方向に移動させる補正において、該格子の格子形状を気液自由界面方向にずらし、ずらされた格子形状と一部分もしくは全部分が重なる格子において、ずらされた該格子の全領域に対する重なる格子領域の割合に該差値を乗算した値を、重ねられた格子における流入補正量とするようにした。
これにより、ずらした方向、すなわち数値拡散が生じた出発点方向に最も近い格子での流入補正量が最も大きくなるため、流体の質量や体積を保存しつつ、数値拡散が生じた出発点に戻すという物理的イメージに沿った方法を構成することが可能となる。
【００１８】
また、本発明の請求項６では、気液自由界面近傍外と判別された領域における正しい流体率値もしくは密度値との差値を気液自由界面方向に移動させる補正において、該格子の距離値から決められる緩和係数を該差値に乗算することで計算を安定化させるようにした。
これにより、気液界面から遠ざかるにつれて緩やかに補正量が上昇するため、移流計算による数値拡散が大きく生じかつ気液自由界面の形状が複雑な場合、急激な補正による計算の不安定化やオーバーシュート、アンダーシュートを抑え、計算の発散を防止する方法を構成することが可能となる。
【００１９】
また、本発明の請求項７では、気液自由界面近傍外と判別された領域における正しい流体率値もしくは密度値との差値を気液自由界面方向に移動させる補正において、気液自由界面近傍であるか否かを判別するあらかじめ指定する距離値を格子間隔の１〜１．５倍とすることで、気液自由界面近傍であると判別される格子の気液自由界面法線方向の格子数を２以内とするようにした。
これにより、気液自由界面近傍における気液界面法線方向の流体率分布の偏りがなくなるため、該補正による該格子の該法線方向の流体率もしくは密度分布のばらつきを防止する方法を構成することが可能となるしたことにある。
【００２０】
また、本発明の請求項８では、保存型移流方程式解法を移流方程式計算に用いるようにした。
これにより、移流方程式の計算において保存性を満足し、且つ数値拡散補正における保存性も満足されるため、いかなる計算サイクルにおいても保存性を満足させる方法を構成することが可能となる。
【００２１】
また、本発明の請求項９では、気液自由界面外格子検出、気液判別および密度補正の一連の手順を、各時刻ステップにおける計算サイクルでの移流方程式計算後に毎サイクル実施するようにした。
これにより、毎サイクル計算後には数値拡散は十分小さくなっていることを保証されるため、いかなる計算サイクルにおいても数値拡散をほとんど生じさせない方法を構成することが可能となる。
【００２２】
【発明の実施の形態】
以下、本発明の詳細実施内容を図面を参照して説明する。
第１図は本発明のアルゴリズムを示す。まず、数１式に示す流体率Ｆに関する次の保存型移流方程式を、保存型ＣＩＰ法やＰＰＭ法を用いて解く。
【００２３】
【数１】
【００２４】
次いで、得られた各格子ごとの流体率を用い、求める気液自由界面からの距離分布Ｇの初期値を、数２式により計算する。
【００２５】
【数２】
【００２６】
ここにαはあらかじめ指定されるパラメータであり、αよりも気液自由界面からの距離が小さい領域では、気液自由界面近傍とみなされる。本計算より、流体率０．５の気液界面の位置では距離Ｇは０であり、気液界面の位置を正しく捉えることができる。
【００２７】
気液自由界面からの距離分布の計算は、繰り返し計算により行う。Ｎ＋１サイクル目のある位置Ｘの格子におけるＧは、Ｎサイクル目のＧを用いて、数３式のように計算する。数３式でｓｇｎ（Ｇ）はＧの符号関数であり、Ｇが正であれば１、Ｇが負であれば−１である。
【００２８】
【数３】
【００２９】
ここにδは微小ベクトルであり、Ｘ＋δはＸから気液自由界面方向に向かって距離｜δ｜離れた参照点を表す。本計算サイクルを繰り返し行うことにより、計算領域全体における気液自由界面からの距離分布を求めることができる。
【００３０】
気液自由界面方向は、距離分布において各軸方向の微係数をとることで決定する。すなわち以下のように数３式のδを決定する。
【００３１】
【数４】
【００３２】
参照点までの距離｜δ｜は、以下の方法で決定する。まず仮の距離｜δ｜を格子間隔Δを越えない値を設定し、数３式のＧを計算する。次に数５式の符号を調べる。
【００３３】
【数５】
【００３４】
符号が負の場合、参照点は気液自由界面を越えて参照していることが分かるため、既に計算したＮ＋１サイクルのＧを破棄し、数６式により｜δ｜を小さくし、参照点をさらに近傍に変更する。
【００３５】
【数６】
【００３６】
ここでγは１より小さい正数である。｜δ｜変更後、再度数３式を計算し、数５式により符号判定を行い、Ｓが正となるまで以上のステップを繰り返す。
【００３７】
数３式右辺におけるＮサイクル目の参照点Ｘ＋δでのＧは、Ｘ＋δは格子点位置ではないので、Ｘ＋δを囲む周囲の格子点のＮサイクル目のＧより補間計算を行う。まず一次元系の場合は、格子ｉと格子ｉ＋１との間に参照点Ｘ＋δが位置すると仮定すると、格子ｉの位置をＸｉ、格子ｉ＋１の位置をＸｉ＋１とおくと、数７式によりｔを計算する。
【００３８】
【数７】
【００３９】
したがって参照点におけるＧは、一次補間式を用いて数８式により計算することができる。
【００４０】
【数８】
【００４１】
二次元系の場合も同様に、格子（ｉ、ｊ）、格子（ｉ＋１、ｊ）、格子（ｉ、ｊ＋１）、格子（ｉ＋１、ｊ＋１）の四点に参照点が囲まれているとすると、数９式、数１０式でｔを計算する。
【００４２】
【数９】
【００４３】
【数１０】
【００４４】
したがって参照点におけるＧは、各軸方向一次補間式を用いて、数１１式のように表される。
【００４５】
【数１１】
【００４６】
三次元系の場合の各軸方向一次補間式も、同様に数１２式で表すことができる。
【００４７】
【数１２】
【００４８】
数３式による気液自由界面との距離分布計算が完了したら、次に各格子ごとに数３式で得られた気液自由界面との距離Ｇおよび数２式のαより、気液自由界面近傍か否かの判別を行う。すなわち該格子におけるＧの絶対値がαよりも大きい場合、該格子は気液自由界面近傍外とみなし、移流方程式計算で生じた数値拡散の補正対象格子とする。
【００４９】
補正対象格子について、Ｇが正ならば液相、負ならば気相であり、前者の場合１−Ｆ、後者の場合０−Ｆ＝−Ｆが補正量ΔＦとなる。すなわち、Ｆに補正量ΔＦを加えた補正後流体率が、液相ならば１、気相ならば０となる。ただし数値拡散が大きい場合補正量が多すぎ、計算が発散したり不安定になる可能性があるため、ΔＦに緩和係数を乗じる。緩和係数は以下の数１３式のように定める。
【００５０】
【数１３】
【００５１】
ｎは正数である。すなわち、｜Ｇ｜＝αの格子点では補正を行わず、αが大きくなるに従い補正量を多くすることで、計算の安定化を図る。
【００５２】
補正対象格子に補正をするのみでは、全体の質量や体積が変化してしまうため、補正による変化分を本来数値拡散が生じた出発点である気液自由界面方向に移動させる。気液自由界面方向を表す単位ベクトルは、数１４式で計算される。
【００５３】
【数１４】
【００５４】
ある格子において、ある距離ｒの分だけ、格子体積が気液自由界面方向に移動したとする。移動後の様子を第２図に示す。もとの格子位置において移動後の格子体積が残っている体積は、三次元系では数１５式で計算される。
【００５５】
【数１５】
【００５６】
ここでΔｘ、Δｙ、ΔｚはＸ、Ｙ、Ｚ軸方向の格子間隔である。したがって、ΔＦが負の場合、補正後における該格子での流量を考慮すると下式が成立し、数１５式のｒを計算することができる。本式はｒに関する三次方程式となるが、ニュートン法を用いることにより、数１６式に示すように、比較的容易に計算できる。
【００５７】
【数１６】
【００５８】
またΔＦが正の場合、補正前における該格子での流量を考慮すると数１７式が成立し、数１５式のｒを計算することができる。
【００５９】
【数１７】
【００６０】
計算したｒより、該格子に隣接した格子において、第２図の移動格子との交差体積を計算でき、移動前後の流量収支を計算することができる。また該格子における流体率はΔＦ加算する。以上より、数値拡散補正前後の総流体体積を変化させることなく、補正を行うことができる。
【００６１】
用いる移流計算スキームの特性によっては、Ｆが微小の負の値、すなわちアンダーシュートを持つことがありうる。その場合、数１６式の計算において、特異な値を持つことになる。その対策として、本数値拡散補正ステップの前段階に、Ｆのアンダーシュートを消すような微小流量を全格子に均一に加え、数値拡散補正ステップ終了後に加えた微小流量を減算する。その場合、補正後流体率は０もしくは１ではなく、それらの値に微小流量加えた分の流体率を加算する。
【００６２】
また、数１５式におけるより簡易な計算方法として、移動体積における流体率をΔＦと考え、数１８式のようにおくことも可能である。この場合ｒを求めるのに、ニュートン法を使用する必要はない。
【００６３】
【数１８】
【００６４】
以上の手順で流体率の数値拡散を補正後、ナビエ・ストークス運動方程式を解き、流速分布を計算する。さらに流体率から密度分布を計算し、圧力ポアソン方程式を解いて、流速分布を修正する。密度分布は、数３式で得られた気液自由界面からの距離Ｇをもとに重み関数を計算して密度を計算する。
【００６５】
以上の手順により、１時刻ステップにおける流速分布、圧力分布および密度分布が求められる。求める時刻に到達するまで、本手順を繰り返し計算する。
【００６６】
【発明の効果】
本発明により、保存型移流方程式解法と数値拡散補正手法との組合せにより、保存性を保ちつつ気液界面の数値拡散を防止することができる。
【００６７】
本発明により、数値拡散補正において複雑な条件判定処理を行うことなく簡易かつ規則的な計算手順により実現できるため、計算速度の向上が実現できる。
【００６８】
本発明により、数値拡散補正において複雑な条件判定処理を行うことなく規則的な計算手順により実現できるため、計算領域を複数に分割して複数のパーソナルコンピュータで並列に計算でき、大幅な計算時間短縮を実現できる。
【００６９】
本発明では、保存型移流方程式解法と数値拡散補正手法とは独立しており、保存型移流方程式解法は任意に選択できるため、ＣＩＰ法やＰＰＭ法などの高精度な移流方程式解法を選択でき、気液二相流体解析を高精度に実現できる。高次補間関数を用いたＣＩＰ法では、数値拡散が非常に少ないため、数値拡散補正による補正量も低減でき、気液自由界面の精度を高くすることが可能である。
【図面の簡単な説明】
【図１】本発明による気液自由界面シミュレーション方法を示すフローチャート図
【図２】本発明による数値拡散補正方法における、数値拡散の移動体積を示す概念図。
【符号の説明】
ΔＸ…Ｘ軸方向格子間隔
ΔＹ…Ｙ軸方向格子間隔
Ｖ’…もとの格子位置において、移動格子が残っている体積
ｎ…気液界面方向法線ベクトル[0001]
TECHNICAL FIELD OF THE INVENTION
The present invention relates to a velocity distribution, a pressure distribution, and a density distribution of a one-dimensional or multidimensional system having two phases of a gas and a liquid and having a gas-liquid free interface in an incompressible flow field where the flow velocity is sufficiently slower than the sound velocity. The present invention relates to a gas-liquid free interface simulation apparatus and method for obtaining a distribution and a gas-liquid free interface shape.
[0002]
In addition, the present invention provides high-speed, high-precision calculation of the interface of the gas-liquid free interface, which is necessary for simulation of the motion, deformation, splitting, and binding of droplets and bubbles, surface flow on solid surfaces, flows in traps and pipes, and waves. The present invention relates to an arithmetic device and a method that can be realized at low cost.
[0003]
[Prior art]
Conventionally, in the gas-liquid free interface fluid simulation, as the most basic method, for example, the fluid rate of each lattice is represented by gas as 0 and the liquid as 1, and the fluid rate distribution is expressed by a non-conservative advection equation obtained from a mass conservation equation. There is a method of advection. In this method, a gas-liquid free interface simulation is attempted by using a method that is as accurate as possible to solve the advection equation.
[0004]
However, in the above-mentioned conventional method, since the accuracy of the advection equation solution is limited, the number of numerical diffusion layers having an intermediate value between 0 and 1 at the gas-liquid interface increases substantially when the number of calculation cycles increases, There is a problem that the gas-liquid interface is blurred, and that a correct curvature cannot be calculated in the calculation of the surface tension calculated from the gas-liquid interface shape.
[0005]
On the other hand, in the VOF method, which is the most general gas-liquid free interface simulation method, after calculating the advection of the fluid ratio, search and reconstruction of the gas-liquid interface are performed by a donor acceptor method or the like. In this method, the numerical diffusion is corrected, and the conservation of mass and volume is also guaranteed (for example, see Non-Patent Documents 1 and 2).
[0006]
However, in the above-mentioned conventional method, complicated condition determination is required for searching and reconstructing the gas-liquid interface, and the calculation program becomes complicated. Further, the algorithm for reconstructing the gas-liquid interface is not based on physical phenomena, and the calculation results vary greatly depending on the configuration method of the calculation program.
[0007]
On the other hand, there is a gas-liquid free interface simulation method combining a higher-order advection calculation solution called the CIP method and a digitizer method. In the CIP method, a cubic expression is used as an interpolation function of an advection equation, and a value at a reference point at the next time step is obtained in a semi-Lagrange manner. In the digitizer method, the fluid ratio is converted into a variable using a tangent function. In this method, a calculation program is very simple and a result with almost no numerical diffusion is obtained (for example, see Non-Patent Document 3).
[0008]
However, the gas-liquid free interface simulation method has a problem that the preservability of mass and volume is not guaranteed, and, for example, there is a problem that minute droplets or thin film liquids disappear as the number of calculation cycles increases.
[0009]
For this reason, by solving the advection equation using a conservative advection equation solution method such as the conservative CIP method or the PPM method with respect to the conventional method described above, the CIP method or the PPM method has There is a method that utilizes high accuracy (for example, see Non-Patent Documents 4 and 5).
[0010]
[Non-Patent Document 1]
C. W. Hirt and B.A. D. Nichols, Volume of fluid (VOF) method for the dynamics of free boundaries, J. Mol. Comp. Phys. , 39, 201 (1981).
[Non-patent document 2]
D. L. Youngs, Time-dependent multi-material flow with large fluid distortion, ｛￥ em Numer. Methods for Fluids Dynamics｝, edited by K.M. W. Morton and M.S. J. Baines, p. 273
(1982).
[Non-Patent Document 3]
T. Yabe and F.S. Xiao, Description of Complex and Sharp Interface during Shock Wave Interaction with Liquid Drop, J. J .; Phys. Soc. Japan, 62, 2537 (1993).
[Non-patent document 4]
T. Yabe, R.A. Tanaka, T .; Nakamura and F.S. Xiao, Exactly conservative semi-Lagrangian scheme (CIP-CSL) in one dimension, Mon. Wea. Rev .. , 129, 332-344. (2001)
[Non-Patent Document 5]
P. Collella and P.M. R. Woodward, The Piecewise Parabolic Method (PPM) for Gas-Dynamic Simulations, J. Amer. Comp. Phys. ｝, 54, 174 (1984).
[0011]
[Problems to be solved by the invention]
However, when the above-described conventional technology configuration is employed, there are the following problems.
A) In the CIP method or the PPM method having a quadratic or cubic interpolation function, numerical diffusion cannot be prevented, and the interface becomes blurred when the number of calculation cycles increases.
B) Numerical diffusion can be greatly reduced by applying the fourth-order higher-order interpolation function to the CIP method, but the numerical diffusion cannot be reduced to 0, and the interface is also blurred when the number of calculation cycles is very large. .
C) When the conservative advection equation solution method is applied, the advection equation is also of the conservative type, so the digitizer method applied in the non-conservative CIP method cannot be used, and the numerical diffusion at the interface cannot be corrected.
D) Regarding the blurred interface, even if the position where the fluid ratio is 0.5 is regarded as the gas-liquid interface, it is not guaranteed that the numerical diffusion at the interface is evenly distributed to both the gas and liquid phases. Therefore, the liquid or gas volume surrounded by the gas-liquid interface is not preserved.
[0012]
It is an object of the present invention to improve a gas-liquid free interface simulation device using the above conservative advection equation solution method, and to provide a simple and highly accurate simulation device that satisfies preservation and prevents numerical diffusion.
[0013]
[Means for Solving the Problems]
In order to solve the above-mentioned problem, the present invention relates to an incompressible flow field whose flow velocity is sufficiently slower than the speed of sound. In a gas-liquid free interface simulation system that simulates the velocity distribution, pressure distribution, density distribution and gas-liquid free interface shape of a three-dimensional system using the Navier-Stokes equation of motion and a fixed orthogonal or non-orthogonal grid, Is calculated at each calculation grid, and it is determined that the calculation grid is not included in the gas-liquid free interface when the distance is larger than a predetermined distance value for determining whether the calculation grid is near the gas-liquid free interface. Gas-liquid free interface outside grid detection means, and gas-liquid discrimination means for determining that the calculation grid is a gas phase or a liquid phase based on the sign of the distance value in the calculation grid. If the fluid rate or density of the lattice obtained by calculating the advection equation is different from the correct fluid rate or density value in the determined gas phase or liquid phase, the difference between the correct value and the gas-liquid free value is calculated. It is provided with a density correcting means for correcting the numerical diffusion of the density caused by the advection equation calculation by moving in the interface direction.
As a result, the calculation grid is corrected to the correct value, and the value of the adjacent grid in the direction of the gas-liquid free interface increases by the difference value, so that the numerical diffusion of the density generated in the advection equation calculation is specified in advance. Numerical diffusion in a region outside the gas-liquid free interface can be removed.
[0014]
Further, in claim 2 of the present invention, in calculating the distance from the gas-liquid free interface in each calculation grid, the intermediate value of the gas phase fluid rate and the liquid phase fluid rate is taken as the value at the gas-liquid free interface, The intermediate value is subtracted from the fluid rate of each calculation grid obtained by calculating the equation, and the value obtained by multiplying by the distance value for determining whether or not it is near the gas-liquid free interface is stored as the initial value of the distance calculation. I did it.
As a result, the distance from the gas-liquid free interface on the iso-surface representing the intermediate value is 0, that is, the gas-liquid free interface itself. Therefore, the initial value for calculating the distance to the free interface without moving the gas-liquid interface. It is possible to configure the setting method.
[0015]
According to a third aspect of the present invention, in a multidimensional system, in a correction for moving a difference value between a correct fluid rate value or a density value in a region determined to be outside the vicinity of a gas-liquid free interface in a gas-liquid free interface direction, The normal vector obtained by taking the derivative in each axial direction with respect to the distance distribution from the gas-liquid free interface is regarded as the gas-liquid free interface direction. As a result, it is guaranteed that convergence in the correct gas-liquid interface direction is achieved as the number of lattices is increased, so that a gas-liquid free interface direction detection method that is theoretically error-free with simple calculations can be configured. .
[0016]
Further, in claim 4 of the present invention, in the calculation of the distance from the gas-liquid free interface in each calculation grid, in the process of repeatedly calculating and converging the distance value, a point near the gas-liquid free interface direction is referred to. The distance value at the neighboring point is interpolated by constructing a linear function in each axis from the distance values of all the grid points surrounding the neighboring point, and the distance between the neighboring point and the computational grid is added. Calculate distance at
In this way, when the grid spacing and the distance are sufficiently small, a complicated gas-liquid free interface shape can be obtained using interpolation calculations by simple axial linear functions in order to ensure that the calculation grid has a substantially accurate distance. It is possible to configure a method of calculating the distance distribution from the gas-liquid free interface without breaking the equation.
[0017]
According to a fifth aspect of the present invention, in the correction for moving a difference value between a correct fluid rate value or a density value in a region determined to be outside the vicinity of the gas-liquid free interface in the gas-liquid free interface direction in a multidimensional system, The lattice shape of the lattice is shifted in the direction of the gas-liquid free interface, and the ratio of the overlapping lattice region to the entire region of the shifted lattice is multiplied by the difference value in the lattice in which the shifted lattice shape partially or entirely overlaps. The value was set as the inflow correction amount in the superimposed grid.
As a result, the amount of inflow correction in the shifted direction, that is, the lattice closest to the starting point direction in which the numerical diffusion has occurred, is the largest, so that the fluid is returned to the starting point in which the numerical diffusion has occurred while preserving the mass and volume of the fluid. It is possible to configure a method according to the physical image.
[0018]
According to a sixth aspect of the present invention, in the correction for moving the difference value between the correct fluid rate value or the density value in the region determined to be outside the vicinity of the gas-liquid free interface in the direction of the gas-liquid free interface, the distance value of the grid is used. The calculation is stabilized by multiplying the difference value by a relaxation coefficient determined from.
As a result, the correction amount gradually increases as the distance from the gas-liquid interface increases, and when numerical diffusion by advection calculation is large and the shape of the gas-liquid free interface is complicated, the calculation becomes unstable due to abrupt correction and overshoot. , A method of suppressing undershoot and preventing divergence of calculation can be configured.
[0019]
According to a seventh aspect of the present invention, in the correction for moving the difference value between the correct fluid rate value or the density value in the region determined to be outside the vicinity of the gas-liquid free interface in the direction of the gas-liquid free interface, By setting the distance value specified in advance to be 1 to 1.5 times the grid interval, the grid in the normal direction of the gas-liquid free interface of the grid determined to be near the gas-liquid free interface is determined. The number was limited to 2 or less.
This eliminates the bias of the fluid ratio distribution in the normal direction of the gas-liquid interface near the gas-liquid free interface, and thus constitutes a method for preventing the fluid rate or the density distribution of the lattice in the normal direction from being varied due to the correction. It has become possible.
[0020]
Further, in claim 8 of the present invention, a conservative advection equation solution method is used for advection equation calculation.
This satisfies the preservation in the calculation of the advection equation and also satisfies the preservation in the numerical diffusion correction, so that it is possible to configure a method that satisfies the preservation in any calculation cycle.
[0021]
According to the ninth aspect of the present invention, a series of procedures of detection of a gas-liquid free interface outside grid, gas-liquid discrimination, and density correction are executed in each cycle after calculating an advection equation in a calculation cycle at each time step.
As a result, since it is guaranteed that the numerical diffusion is sufficiently small after each cycle calculation, it is possible to configure a method in which the numerical diffusion hardly occurs in any calculation cycle.
[0022]
BEST MODE FOR CARRYING OUT THE INVENTION
Hereinafter, details of the present invention will be described with reference to the drawings.
FIG. 1 shows the algorithm of the present invention. First, the following conservative advection equation for the fluid rate F shown in Equation 1 is solved by using the conservative CIP method or the PPM method.
[0023]
(Equation 1)
Next, using the obtained fluid ratio for each lattice, the initial value of the distance distribution G from the gas-liquid free interface to be calculated is calculated by Expression 2.
[0025]
(Equation 2)
Here, α is a parameter specified in advance, and in a region where the distance from the gas-liquid free interface is smaller than α, it is regarded as the vicinity of the gas-liquid free interface. From this calculation, the distance G is 0 at the position of the gas-liquid interface with the fluid rate of 0.5, and the position of the gas-liquid interface can be correctly grasped.
[0027]
Calculation of the distance distribution from the gas-liquid free interface is performed by iterative calculation. G in the lattice at a certain position X in the (N + 1) th cycle is calculated as shown in Expression 3 using G in the Nth cycle. In Equation 3, sgn (G) is the sign function of G, and is 1 if G is positive and -1 if G is negative.
[0028]
(Equation 3)
Here, δ is a minute vector, and X + δ represents a reference point at a distance | δ | away from X in the direction of the gas-liquid free interface. By repeating this calculation cycle, the distance distribution from the gas-liquid free interface in the entire calculation region can be obtained.
[0030]
The gas-liquid free interface direction is determined by taking a differential coefficient in each axial direction in the distance distribution. That is, δ of Expression 3 is determined as follows.
[0031]
(Equation 4)
The distance | δ | to the reference point is determined by the following method. First, the temporary distance | δ | is set to a value that does not exceed the lattice spacing Δ, and G in Equation 3 is calculated. Next, the sign of Equation 5 is examined.
[0033]
(Equation 5)
If the sign is negative, it can be seen that the reference point refers beyond the gas-liquid free interface, so the previously calculated N + 1 cycle G is discarded, and | δ | Further, it is changed to the vicinity.
[0035]
(Equation 6)
Here, γ is a positive number smaller than 1. After changing | δ |, the equation (3) is calculated again, the sign is determined by the equation (5), and the above steps are repeated until S becomes positive.
[0037]
Regarding G at the reference point X + δ at the Nth cycle on the right side of Equation 3, since X + δ is not a grid point position, interpolation calculation is performed from G at the Nth cycle of surrounding grid points surrounding X + δ. First, in the case of a one-dimensional system, assuming that the reference point X + δ is located between the grid i and the grid i + 1, if the position of the grid i is Xi and the position of the grid i + 1 is Xi + 1, t is calculated by Expression 7. I do.
[0038]
(Equation 7)
Therefore, G at the reference point can be calculated by equation 8 using a linear interpolation equation.
[0040]
(Equation 8)
Similarly, in the case of the two-dimensional system, assuming that the reference point is surrounded by four points of the grid (i, j), the grid (i + 1, j), the grid (i, j + 1), and the grid (i + 1, j + 1). T is calculated by Expressions 9 and 10.
[0042]
(Equation 9)
(Equation 10)
Therefore, G at the reference point is expressed by Equation 11 using the primary interpolation formula in each axis direction.
[0045]
(Equation 11)
In the case of a three-dimensional system, the primary interpolation formula in each axis direction can be similarly expressed by Expression 12.
[0047]
(Equation 12)
After the calculation of the distance distribution to the gas-liquid free interface by Equation 3 is completed, the gas-liquid free interface is then calculated from the distance G to the gas-liquid free interface obtained by Equation 3 for each lattice and α in Equation 2 It is determined whether it is near or not. That is, when the absolute value of G in the lattice is larger than α, the lattice is regarded as outside the vicinity of the gas-liquid free interface, and is set as a correction target lattice for numerical diffusion generated by advection equation calculation.
[0049]
Regarding the correction target grid, if G is positive, it is a liquid phase, and if it is negative, it is a gas phase. In the former case, 1-F, and in the latter case, 0-F = -F is the correction amount ΔF. That is, the corrected fluid ratio obtained by adding the correction amount ΔF to F is 1 for a liquid phase and 0 for a gas phase. However, when the numerical diffusion is large, the correction amount is too large, and the calculation may diverge or become unstable. Therefore, ΔF is multiplied by the relaxation coefficient. The relaxation coefficient is determined as in the following equation (13).
[0050]
(Equation 13)
n is a positive number. That is, the correction is not performed at the lattice point of | G | = α, and the correction amount is increased as α increases, thereby stabilizing the calculation.
[0052]
Simply correcting the grid to be corrected changes the overall mass and volume, so the amount of change due to the correction is moved toward the gas-liquid free interface, which is the starting point where the numerical diffusion originally occurred. The unit vector representing the direction of the gas-liquid free interface is calculated by Expression 14.
[0053]
[Equation 14]
In a certain grid, it is assumed that the grid volume has moved toward the gas-liquid free interface by a certain distance r. The state after the movement is shown in FIG. The volume in which the moved lattice volume remains at the original lattice position is calculated by Equation 15 in a three-dimensional system.
[0055]
(Equation 15)
Here, Δx, Δy, and Δz are lattice intervals in the X, Y, and Z axis directions. Therefore, when ΔF is negative, the following expression is satisfied when the flow rate in the grid after the correction is considered, and r of Expression 15 can be calculated. This equation is a cubic equation for r, but can be calculated relatively easily by using Newton's method, as shown in equation (16).
[0057]
(Equation 16)
When ΔF is positive, Expression 17 is satisfied when the flow rate in the grid before correction is considered, and r of Expression 15 can be calculated.
[0059]
[Equation 17]
From the calculated r, the volume of intersection with the moving grid in FIG. 2 can be calculated in the grid adjacent to the grid, and the flow balance before and after the movement can be calculated. Further, the fluid rate in the lattice is added by ΔF. As described above, the correction can be performed without changing the total fluid volume before and after the numerical diffusion correction.
[0061]
Depending on the characteristics of the advection calculation scheme used, F may have a small negative value, ie, undershoot. In that case, in the calculation of the expression 16, it has a unique value. As a countermeasure, a minute flow that eliminates the undershoot of F is uniformly applied to all the lattices before the numerical diffusion correction step, and the minute flow added after the end of the numerical diffusion correction step is subtracted. In this case, the fluid rate after correction is not 0 or 1, and the fluid rate corresponding to the minute flow rate is added to those values.
[0062]
Further, as a simpler calculation method in Expression 15, it is possible to consider the fluid rate in the moving volume as ΔF and to set it as in Expression 18. In this case, it is not necessary to use Newton's method to determine r.
[0063]
(Equation 18)
After correcting the numerical diffusion of the fluid rate by the above procedure, solve the Navier-Stokes equation of motion and calculate the flow velocity distribution. Further, the density distribution is calculated from the fluid fraction, and the pressure Poisson equation is solved to correct the flow velocity distribution. The density distribution is calculated by calculating a weighting function based on the distance G from the gas-liquid free interface obtained by Expression 3 to calculate the density.
[0065]
According to the above procedure, the flow velocity distribution, pressure distribution, and density distribution at one time step are obtained. Repeat this procedure until the required time is reached.
[0066]
【The invention's effect】
According to the present invention, it is possible to prevent the numerical diffusion at the gas-liquid interface while maintaining the preservation by combining the conservation type advection equation solving method and the numerical diffusion correction method.
[0067]
According to the present invention, the numerical diffusion correction can be realized by a simple and regular calculation procedure without performing a complicated condition determination process, so that the calculation speed can be improved.
[0068]
According to the present invention, since the numerical diffusion correction can be realized by a regular calculation procedure without performing complicated condition determination processing, the calculation area can be divided into a plurality of sections and calculations can be performed in parallel by a plurality of personal computers, thereby greatly reducing the calculation time. Can be realized.
[0069]
In the present invention, the conservation type advection equation solution and the numerical diffusion correction method are independent, and the conservation type advection equation solution can be arbitrarily selected. Gas-liquid two-phase fluid analysis can be realized with high accuracy. In the CIP method using the higher-order interpolation function, since the numerical diffusion is very small, the correction amount by the numerical diffusion correction can be reduced, and the accuracy of the gas-liquid free interface can be increased.
[Brief description of the drawings]
FIG. 1 is a flowchart illustrating a gas-liquid free interface simulation method according to the present invention. FIG. 2 is a conceptual diagram illustrating a moving volume of numerical diffusion in a numerical diffusion correction method according to the present invention.
[Explanation of symbols]
ΔX: X-axis direction lattice spacing ΔY: Y-axis direction lattice spacing V ′: Volume at which the moving lattice remains at the original lattice position n: Gas-liquid interface direction normal vector
