(57)【要約】【課題】本発明の目的は、液滴や気泡の運動、変形、分裂、結合、固体表面の表層流れ、トラップ・配管内流れおよび波のシミュレーションにおいて必要となる気液自由界面の界面計算を高速、高精度かつ低価格に実現するシミュレーション装置を提供する。【解決手段】気液自由界面からの距離を各計算格子において計算し、指定値よりも距離が大きくかつ保存型解法により移流方程式を計算することで得られた該格子の流体率もしくは密度が数値拡散により誤差を生じた場合、正しい値との差値を気液自由界面方向に移動させることで数値拡散を補正する。【選択図】　図１
2004118394.tif 000002
【特許請求の範囲】
【請求項１】
流速が音速に比べて十分に遅い非圧縮流れの場において、気体と液体との二相からなりかつ気液自由界面を有する一次元もしくは多次元系の速度分布、圧力分布、密度分布および気液自由界面形状をナビエ・ストークス運動方程式および固定された直交もしくは非直交格子を用いてシミュレーションする気液自由界面シミュレーション装置において、気液自由界面からの距離を各計算格子において計算し、気液自由界面近傍であるか否かを判別するあらかじめ指定された距離値よりも大きい場合に、当計算格子は気液自由界面に含まれないと判別される気液自由界面外格子検出手段と、該計算格子における該距離値の符号により該計算格子が気相もしくは液相であると判別される気液判別手段と、移流方程式を計算することで得られた該格子の流体率もしくは密度が、判別された気相もしくは液相での正しい流体率値もしくは密度値と異なる場合、正しい値との差値を気液自由界面方向に移動させることで、移流方程式計算で生じた密度の数値拡散を補正する密度補正手段とを具えたことを特徴とする気液自由界面シミュレーション装置。
【請求項２】
請求項１において、各計算格子における気液自由界面からの距離の計算において、気相の流体率と液相の流体率の中間値を気液自由界面での値とし、移流方程式を計算することで得られた各計算格子の流体率から該中間値を引き、気液自由界面近傍であるか否かを判別する該距離値を乗算した値を距離計算の初期値として格納することを特徴とする請求項１に記載の気液自由界面シミュレーション装置。
【請求項３】
請求項１の多次元系において、正しい値との差値を気液自由界面方向に移動させる補正において、気液自由界面からの距離分布に対し各軸方向の微係数をとることで得られる法線ベクトルを気液自由界面方向とみなすことを特徴とする請求項１に記載の気液自由界面シミュレーション装置。
【請求項４】
請求項１において、各計算格子における気液自由界面からの距離の計算において、繰り返し計算を行い距離値を収束させる過程で、請求項３に記載の法線ベクトル方向の近傍点を参照し、該近傍点での距離値を該近傍点を囲む全ての格子点の距離値より各軸方向一次関数を構成して補間計算し、該近傍点と該計算格子との距離を加えることで該計算格子での距離を計算することを特徴とする請求項１に記載の気液自由界面シミュレーション装置。
【請求項５】
請求項１の多次元系において、正しい値との差値を気液自由界面方向に移動させる補正において、該格子の格子形状を気液自由界面方向にずらし、ずらされた格子形状と一部分もしくは全部分が重なる格子において、ずらされた該格子の全領域に対する重なる格子領域の割合に該差値を乗算した値を、重ねられた格子における流入補正量とすることを特徴とする請求項１に記載の気液自由界面シミュレーション装置。
【請求項６】
請求項１において、正しい値との差値を気液自由界面方向に移動させる補正において、該格子の距離値から決められる緩和係数を該差値に乗算することで計算を安定化させることを特徴とする請求項１に記載の気液自由界面シミュレーション装置。
【請求項７】
請求項１において、正しい値との差値を気液自由界面方向に移動させる補正において、気液自由界面近傍であるか否かを判別するあらかじめ指定する距離値を格子間隔の１〜１．５倍とすることを特徴とする請求項１に記載の気液自由界面シミュレーション装置。
【請求項８】
請求項１において、保存型移流方程式解法を移流方程式計算に用いることを特徴とする請求項１に記載の気液自由界面シミュレーション装置。
【請求項９】
請求項１において、該気液自由界面外格子検出、気液判別および密度補正の一連の手順を、各時刻ステップにおける計算サイクルでの移流方程式計算後に毎サイクル実施することを特徴とする請求項１に記載の気液自由界面シミュレーション装置。
【発明の詳細な説明】
【０００１】
【発明の属する技術分野】
本発明は、流速が音速に比べて十分に遅い非圧縮流れの場において、気体と液体との二相からなりかつ気液自由界面を有する一次元もしくは多次元系の速度分布、圧力分布、密度分布および気液自由界面形状を求める気液自由界面シミュレーション装置および方法に関する。
【０００２】
また、本発明は、液滴や気泡の運動、変形、分裂、結合、固体表面の表層流れ、トラップ・配管内流れおよび波のシミュレーションにおいて必要となる気液自由界面の界面計算を高速、高精度かつ安価に実現する演算装置および方法に関する。
【０００３】
【従来の技術】
従来、気液自由界面流体シミュレーションでは、最も基本的な手法として、例えば各格子の流体率を気体を０、液体を１としてあらわし、流体率分布を質量保存式から得られる非保存型移流方程式で移流させる方法がある。本手法では、移流方程式の解法になるべく精度の高い方法を用いることで、気液自由界面シミュレーションを試みる。
【０００４】
しかしながら前記の従来手法では、移流方程式解法の精度に限界があるため、実際的には計算サイクルが増えると気液界面に０と１の中間的な値を持つ数値拡散層が大幅に増大し、気液界面がほやける課題を有し、且つ、気液界面形状から計算される表面張力の計算において正しい曲率計算ができなくなるという問題点を有する。
【０００５】
一方、最も一般的な気液自由界面シミュレーション手法であるＶＯＦ法では、流体率の移流計算後、ドナーアクセプタ法などで気液界面の探索や再構築を行う。本手法では、数値拡散は補正され、質量や体積の保存性も保証される（例えば、非特許文献１，２参照）。
【０００６】
しかしながら、前記の従来手法においては、気液界面の探索や再構築に複雑な条件判定が必要となり、計算プログラムが複雑化する。また気液界面の再構築アルゴリズムは物理的現象に基づくものとはいえず、計算プログラムの構成方法次第で計算結果が大きく変わるため、計算精度においても課題を有する。
【０００７】
他方、ＣＩＰ法とよぶ高次精度移流計算解法およびデジタイザ法を組み合わせた気液自由界面シミュレーション方法がある。ＣＩＰ法では移流方程式の補間関数に三次式を用い、次の時刻ステップでの参照点における値をセミラグランジェ的に持ってくる。またデジタイザ法ではタンジェント関数を用いて流体率を変数変換する。本手法においては、計算プログラムは極めて簡易でありかつ数値拡散もほとんどない結果が得られる（例えば、非特許文献３参照）。
【０００８】
しかし、前記気液自由界面シミュレーション方法においては、質量や体積の保存性が保証されず、例えば微少の液滴や薄膜の液体が、計算サイクル数が多くなるに従い消失していくといった課題を有する。
【０００９】
このため、上記のような従来方法に対して、保存型ＣＩＰ法やＰＰＭ法などの保存型移流方程式解法を用いて移流方程式を解くことで、保存性を満足しつつＣＩＰ法やＰＰＭ法の持つ高い精度を利用する方法がある（例えば、非特許文献４，５参照）。
【００１０】
【非特許文献１】
Ｃ．Ｗ．Ｈｉｒｔ　ａｎｄ　Ｂ．Ｄ．Ｎｉｃｈｏｌｓ，　Ｖｏｌｕｍｅ　ｏｆ　ｆｌｕｉｄ　（ＶＯＦ）　ｍｅｔｈｏｄ　ｆｏｒ　ｔｈｅ　ｄｙｎａｍｉｃｓ　ｏｆ　ｆｒｅｅ　ｂｏｕｎｄａｒｉｅｓ，　　Ｊ．　Ｃｏｍｐ．　Ｐｈｙｓ．，　　３９，　２０１（１９８１）．
【非特許文献２】
Ｄ．Ｌ．Ｙｏｕｎｇｓ，　Ｔｉｍｅ−ｄｅｐｅｎｄｅｎｔ　ｍｕｌｔｉ−ｍａｔｅｒｉａｌ　ｆｌｏｗ　ｗｉｔｈ　ｌａｒｇｅ　ｆｌｕｉｄ　ｄｉｓｔｏｒｔｉｏｎ，　｛￥ｅｍ　Ｎｕｍｅｒ．　Ｍｅｔｈｏｄｓ　ｆｏｒ　Ｆｌｕｉｄｓ　Ｄｙｎａｍｉｃｓ｝，　ｅｄｉｔｅｄ　ｂｙ　Ｋ．Ｗ．Ｍｏｒｔｏｎ　ａｎｄ　Ｍ．Ｊ．Ｂａｉｎｅｓ，　ｐ２７３
（１９８２）．
【非特許文献３】
Ｔ．Ｙａｂｅ　ａｎｄ　Ｆ．Ｘｉａｏ，　Ｄｅｓｃｒｉｐｔｉｏｎ　ｏｆ　Ｃｏｍｐｌｅｘ　ａｎｄ　Ｓｈａｒｐ　Ｉｎｔｅｒｆａｃｅ　ｄｕｒｉｎｇ　Ｓｈｏｃｋ　Ｗａｖｅ　Ｉｎｔｅｒａｃｔｉｏｎ　ｗｉｔｈ　Ｌｉｑｕｉｄ　Ｄｒｏｐ，　Ｊ．　Ｐｈｙｓ．　Ｓｏｃ．　Ｊａｐａｎ，　　６２，　２５３７（１９９３）．
【非特許文献４】
Ｔ．Ｙａｂｅ，　Ｒ．Ｔａｎａｋａ，　Ｔ．Ｎａｋａｍｕｒａ　ａｎｄ　Ｆ．Ｘｉａｏ，　Ｅｘａｃｔｌｙ　ｃｏｎｓｅｒｖａｔｉｖｅ　ｓｅｍｉ−Ｌａｇｒａｎｇｉａｎ　ｓｃｈｅｍｅ　（ＣＩＰ−ＣＳＬ）　ｉｎ　ｏｎｅ　ｄｉｍｅｎｓｉｏｎ，　Ｍｏｎ．Ｗｅａ．Ｒｅｖ．，１２９，　３３２−３４４．（２００１）
【非特許文献５】
Ｐ．Ｃｏｌｅｌｌａ　ａｎｄ　Ｐ．Ｒ．Ｗｏｏｄｗａｒｄ，　Ｔｈｅ　ｐｉｅｃｅｗｉｓｅ　ｐａｒａｂｏｌｉｃ　ｍｅｔｈｏｄ　（ＰＰＭ）　ｆｏｒ　ｇａｓ−ｄｙｎａｍｉｃａｌ　ｓｉｍｕｌａｔｉｏｎｓ，　Ｊ．　Ｃｏｍｐ．　Ｐｈｙｓ．｝，　　５４，　１７４（１９８４）．
【００１１】
【発明が解決しようとする課題】
しかしながら、上記のような従来技術の構成をとった場合、以下のような課題がある。
ア）二次もしくは三次の補間関数をもつＣＩＰ法やＰＰＭ法では、数値拡散が防止できず、計算サイクル数が増大すると、界面がぼやける。
イ）四次の高次補間関数をＣＩＰ法に適用することで大幅に数値拡散を低減できるものの、数値拡散を０にすることはできず、計算サイクル数が非常に大きい場合、やはり界面がぼやける。
ウ）保存型移流方程式解法を適用する場合は移流方程式も保存型となるため、非保存型ＣＩＰ法で適用するデジタイザ法を用いることはできず、界面の数値拡散を補正することができない。
エ）ぼやけた界面に対し、例えば流体率が０．５の位置を気液界面とみなしても、界面の数値拡散は気相・液相両側に対し均等に数値拡散することは保証されていないため、気液界面で囲まれた液体もしくは気体体積は保存されない。
【００１２】
ここに本発明は、上記保存型移流方程式解法を用いた気液自由界面シミュレーション装置を改良し、保存性を満足しつつ数値拡散を防止する簡易かつ高精度なシミュレーション装置を提供することにある。
【００１３】
【課題を解決するための手段】
上記課題を解決するため本発明の請求項１は、流速が音速に比べて十分に遅い非圧縮流れの場において、気体と液体との二相からなりかつ気液自由界面を有する一次元もしくは多次元系の速度分布、圧力分布、密度分布および気液自由界面形状をナビエ・ストークス運動方程式および固定された直交もしくは非直交格子を用いてシミュレーションする気液自由界面シミュレーション装置において、気液自由界面からの距離を各計算格子において計算し、気液自由界面近傍であるか否かを判別するあらかじめ指定された距離値よりも大きい場合に、当計算格子は気液自由界面に含まれないと判別される気液自由界面外格子検出手段と、該計算格子における該距離値の符号により該計算格子が気相もしくは液相であると判別される気液判別手段と、移流方程式を計算することで得られた該格子の流体率もしくは密度が、判別された気相もしくは液相での正しい流体率値もしくは密度値と異なる場合、正しい値との差値を気液自由界面方向に移動させることで、移流方程式計算で生じた密度の数値拡散を補正する密度補正手段とを具えるようにした。
これにより、当計算格子は正しい値に修正され、また気液自由界面方向の隣接格子の値が該差値の分上昇するため、移流方程式計算で生じた密度の数値拡散に対し、あらかじめ指定した気液自由界面外の領域での数値拡散を除去することが可能となる。
【００１４】
また、本発明の請求項２では、各計算格子における気液自由界面からの距離の計算において、気相の流体率と液相の流体率の中間値を気液自由界面での値とし、移流方程式を計算することで得られた各計算格子の流体率から該中間値を引き、気液自由界面近傍であるか否かを判別する該距離値を乗算した値を距離計算の初期値として格納するようにした。
これにより、該中間値を表わす等値面上では気液自由界面からの距離が０、すなわち気液自由界面そのものとなるため、気液界面を動かさずに自由界面との距離を計算する初期値設定方法を構成することが可能となる。
【００１５】
また、本発明の請求項３では、多次元系において、気液自由界面近傍外と判別された領域における正しい流体率値もしくは密度値との差値を気液自由界面方向に移動させる補正において、気液自由界面からの距離分布に対し各軸方向の微係数をとることで得られる法線ベクトルを気液自由界面方向とみなすようにした。これにより、格子数を増やしていくと正しい気液界面方向に収束することが保証されるため、簡易な計算ながら理論的に誤りのない気液自由界面方向検出方法を構成することが可能となる。
【００１６】
また、本発明の請求項４では、各計算格子における気液自由界面からの距離の計算において、繰り返し計算を行い距離値を収束させる過程で、気液自由界面方向の近傍点を参照し、該近傍点での距離値を該近傍点を囲む全ての格子点の距離値より各軸方向一次関数を構成して補間計算し、該近傍点と該計算格子との距離を加えることで該計算格子での距離を計算するようにした。
これにより、格子間隔および該距離が十分小さい場合、該計算格子においてほぼ正確な距離となることを保証するため、簡易な各軸方向一次関数による補間計算を用いて、複雑な気液自由界面形状を崩さず気液自由界面からの距離分布を計算する方法を構成することが可能となる。
【００１７】
また、本発明の請求項５では、多次元系において、気液自由界面近傍外と判別された領域における正しい流体率値もしくは密度値との差値を気液自由界面方向に移動させる補正において、該格子の格子形状を気液自由界面方向にずらし、ずらされた格子形状と一部分もしくは全部分が重なる格子において、ずらされた該格子の全領域に対する重なる格子領域の割合に該差値を乗算した値を、重ねられた格子における流入補正量とするようにした。
これにより、ずらした方向、すなわち数値拡散が生じた出発点方向に最も近い格子での流入補正量が最も大きくなるため、流体の質量や体積を保存しつつ、数値拡散が生じた出発点に戻すという物理的イメージに沿った方法を構成することが可能となる。
【００１８】
また、本発明の請求項６では、気液自由界面近傍外と判別された領域における正しい流体率値もしくは密度値との差値を気液自由界面方向に移動させる補正において、該格子の距離値から決められる緩和係数を該差値に乗算することで計算を安定化させるようにした。
これにより、気液界面から遠ざかるにつれて緩やかに補正量が上昇するため、移流計算による数値拡散が大きく生じかつ気液自由界面の形状が複雑な場合、急激な補正による計算の不安定化やオーバーシュート、アンダーシュートを抑え、計算の発散を防止する方法を構成することが可能となる。
【００１９】
また、本発明の請求項７では、気液自由界面近傍外と判別された領域における正しい流体率値もしくは密度値との差値を気液自由界面方向に移動させる補正において、気液自由界面近傍であるか否かを判別するあらかじめ指定する距離値を格子間隔の１〜１．５倍とすることで、気液自由界面近傍であると判別される格子の気液自由界面法線方向の格子数を２以内とするようにした。
これにより、気液自由界面近傍における気液界面法線方向の流体率分布の偏りがなくなるため、該補正による該格子の該法線方向の流体率もしくは密度分布のばらつきを防止する方法を構成することが可能となるしたことにある。
【００２０】
また、本発明の請求項８では、保存型移流方程式解法を移流方程式計算に用いるようにした。
これにより、移流方程式の計算において保存性を満足し、且つ数値拡散補正における保存性も満足されるため、いかなる計算サイクルにおいても保存性を満足させる方法を構成することが可能となる。
【００２１】
また、本発明の請求項９では、気液自由界面外格子検出、気液判別および密度補正の一連の手順を、各時刻ステップにおける計算サイクルでの移流方程式計算後に毎サイクル実施するようにした。
これにより、毎サイクル計算後には数値拡散は十分小さくなっていることを保証されるため、いかなる計算サイクルにおいても数値拡散をほとんど生じさせない方法を構成することが可能となる。
【００２２】
【発明の実施の形態】
以下、本発明の詳細実施内容を図面を参照して説明する。
第１図は本発明のアルゴリズムを示す。まず、数１式に示す流体率Ｆに関する次の保存型移流方程式を、保存型ＣＩＰ法やＰＰＭ法を用いて解く。
【００２３】
【数１】
2004118394.tif 000003
【００２４】
次いで、得られた各格子ごとの流体率を用い、求める気液自由界面からの距離分布Ｇの初期値を、数２式により計算する。
【００２５】
【数２】
2004118394.tif 000004
【００２６】
ここにαはあらかじめ指定されるパラメータであり、αよりも気液自由界面からの距離が小さい領域では、気液自由界面近傍とみなされる。本計算より、流体率０．５の気液界面の位置では距離Ｇは０であり、気液界面の位置を正しく捉えることができる。
【００２７】
気液自由界面からの距離分布の計算は、繰り返し計算により行う。Ｎ＋１サイクル目のある位置Ｘの格子におけるＧは、Ｎサイクル目のＧを用いて、数３式のように計算する。数３式でｓｇｎ（Ｇ）はＧの符号関数であり、Ｇが正であれば１、Ｇが負であれば−１である。
【００２８】
【数３】
2004118394.tif 000005
【００２９】
ここにδは微小ベクトルであり、Ｘ＋δはＸから気液自由界面方向に向かって距離｜δ｜離れた参照点を表す。本計算サイクルを繰り返し行うことにより、計算領域全体における気液自由界面からの距離分布を求めることができる。
【００３０】
気液自由界面方向は、距離分布において各軸方向の微係数をとることで決定する。すなわち以下のように数３式のδを決定する。
【００３１】
【数４】
2004118394.tif 000006
【００３２】
参照点までの距離｜δ｜は、以下の方法で決定する。まず仮の距離｜δ｜を格子間隔Δを越えない値を設定し、数３式のＧを計算する。次に数５式の符号を調べる。
【００３３】
【数５】
2004118394.tif 000007
【００３４】
符号が負の場合、参照点は気液自由界面を越えて参照していることが分かるため、既に計算したＮ＋１サイクルのＧを破棄し、数６式により｜δ｜を小さくし、参照点をさらに近傍に変更する。
【００３５】
【数６】
2004118394.tif 000008
【００３６】
ここでγは１より小さい正数である。｜δ｜変更後、再度数３式を計算し、数５式により符号判定を行い、Ｓが正となるまで以上のステップを繰り返す。
【００３７】
数３式右辺におけるＮサイクル目の参照点Ｘ＋δでのＧは、Ｘ＋δは格子点位置ではないので、Ｘ＋δを囲む周囲の格子点のＮサイクル目のＧより補間計算を行う。まず一次元系の場合は、格子ｉと格子ｉ＋１との間に参照点Ｘ＋δが位置すると仮定すると、格子ｉの位置をＸｉ、格子ｉ＋１の位置をＸｉ＋１とおくと、数７式によりｔを計算する。
【００３８】
【数７】
2004118394.tif 000009
【００３９】
したがって参照点におけるＧは、一次補間式を用いて数８式により計算することができる。
【００４０】
【数８】
2004118394.tif 000010
【００４１】
二次元系の場合も同様に、格子（ｉ、ｊ）、格子（ｉ＋１、ｊ）、格子（ｉ、ｊ＋１）、格子（ｉ＋１、ｊ＋１）の四点に参照点が囲まれているとすると、数９式、数１０式でｔを計算する。
【００４２】
【数９】
2004118394.tif 000011
【００４３】
【数１０】
2004118394.tif 000012
【００４４】
したがって参照点におけるＧは、各軸方向一次補間式を用いて、数１１式のように表される。
【００４５】
【数１１】
2004118394.tif 000013
【００４６】
三次元系の場合の各軸方向一次補間式も、同様に数１２式で表すことができる。
【００４７】
【数１２】
2004118394.tif 000014
【００４８】
数３式による気液自由界面との距離分布計算が完了したら、次に各格子ごとに数３式で得られた気液自由界面との距離Ｇおよび数２式のαより、気液自由界面近傍か否かの判別を行う。すなわち該格子におけるＧの絶対値がαよりも大きい場合、該格子は気液自由界面近傍外とみなし、移流方程式計算で生じた数値拡散の補正対象格子とする。
【００４９】
補正対象格子について、Ｇが正ならば液相、負ならば気相であり、前者の場合１−Ｆ、後者の場合０−Ｆ＝−Ｆが補正量ΔＦとなる。すなわち、Ｆに補正量ΔＦを加えた補正後流体率が、液相ならば１、気相ならば０となる。ただし数値拡散が大きい場合補正量が多すぎ、計算が発散したり不安定になる可能性があるため、ΔＦに緩和係数を乗じる。緩和係数は以下の数１３式のように定める。
【００５０】
【数１３】
2004118394.tif 000015
【００５１】
ｎは正数である。すなわち、｜Ｇ｜＝αの格子点では補正を行わず、αが大きくなるに従い補正量を多くすることで、計算の安定化を図る。
【００５２】
補正対象格子に補正をするのみでは、全体の質量や体積が変化してしまうため、補正による変化分を本来数値拡散が生じた出発点である気液自由界面方向に移動させる。気液自由界面方向を表す単位ベクトルは、数１４式で計算される。
【００５３】
【数１４】
2004118394.tif 000016
【００５４】
ある格子において、ある距離ｒの分だけ、格子体積が気液自由界面方向に移動したとする。移動後の様子を第２図に示す。もとの格子位置において移動後の格子体積が残っている体積は、三次元系では数１５式で計算される。
【００５５】
【数１５】
2004118394.tif 000017
【００５６】
ここでΔｘ、Δｙ、ΔｚはＸ、Ｙ、Ｚ軸方向の格子間隔である。したがって、ΔＦが負の場合、補正後における該格子での流量を考慮すると下式が成立し、数１５式のｒを計算することができる。本式はｒに関する三次方程式となるが、ニュートン法を用いることにより、数１６式に示すように、比較的容易に計算できる。
【００５７】
【数１６】
2004118394.tif 000018
【００５８】
またΔＦが正の場合、補正前における該格子での流量を考慮すると数１７式が成立し、数１５式のｒを計算することができる。
【００５９】
【数１７】
2004118394.tif 000019
【００６０】
計算したｒより、該格子に隣接した格子において、第２図の移動格子との交差体積を計算でき、移動前後の流量収支を計算することができる。また該格子における流体率はΔＦ加算する。以上より、数値拡散補正前後の総流体体積を変化させることなく、補正を行うことができる。
【００６１】
用いる移流計算スキームの特性によっては、Ｆが微小の負の値、すなわちアンダーシュートを持つことがありうる。その場合、数１６式の計算において、特異な値を持つことになる。その対策として、本数値拡散補正ステップの前段階に、Ｆのアンダーシュートを消すような微小流量を全格子に均一に加え、数値拡散補正ステップ終了後に加えた微小流量を減算する。その場合、補正後流体率は０もしくは１ではなく、それらの値に微小流量加えた分の流体率を加算する。
【００６２】
また、数１５式におけるより簡易な計算方法として、移動体積における流体率をΔＦと考え、数１８式のようにおくことも可能である。この場合ｒを求めるのに、ニュートン法を使用する必要はない。
【００６３】
【数１８】
2004118394.tif 000020
【００６４】
以上の手順で流体率の数値拡散を補正後、ナビエ・ストークス運動方程式を解き、流速分布を計算する。さらに流体率から密度分布を計算し、圧力ポアソン方程式を解いて、流速分布を修正する。密度分布は、数３式で得られた気液自由界面からの距離Ｇをもとに重み関数を計算して密度を計算する。
【００６５】
以上の手順により、１時刻ステップにおける流速分布、圧力分布および密度分布が求められる。求める時刻に到達するまで、本手順を繰り返し計算する。
【００６６】
【発明の効果】
本発明により、保存型移流方程式解法と数値拡散補正手法との組合せにより、保存性を保ちつつ気液界面の数値拡散を防止することができる。
【００６７】
本発明により、数値拡散補正において複雑な条件判定処理を行うことなく簡易かつ規則的な計算手順により実現できるため、計算速度の向上が実現できる。
【００６８】
本発明により、数値拡散補正において複雑な条件判定処理を行うことなく規則的な計算手順により実現できるため、計算領域を複数に分割して複数のパーソナルコンピュータで並列に計算でき、大幅な計算時間短縮を実現できる。
【００６９】
本発明では、保存型移流方程式解法と数値拡散補正手法とは独立しており、保存型移流方程式解法は任意に選択できるため、ＣＩＰ法やＰＰＭ法などの高精度な移流方程式解法を選択でき、気液二相流体解析を高精度に実現できる。高次補間関数を用いたＣＩＰ法では、数値拡散が非常に少ないため、数値拡散補正による補正量も低減でき、気液自由界面の精度を高くすることが可能である。
【図面の簡単な説明】
【図１】本発明による気液自由界面シミュレーション方法を示すフローチャート図
【図２】本発明による数値拡散補正方法における、数値拡散の移動体積を示す概念図。
【符号の説明】
ΔＸ…Ｘ軸方向格子間隔
ΔＹ…Ｙ軸方向格子間隔
Ｖ’…もとの格子位置において、移動格子が残っている体積
ｎ…気液界面方向法線ベクトル
【図１】
2004118394.tif 000021
【図２】
2004118394.tif 000022