(57)【特許請求の範囲】
【請求項１】
２種類以上の流体からなる混相流における自由界面をシミュレーションするためのシミュレーション装置であって、
シミュレーション対象が分割された複数のセルのそれぞれにおける前記２種類以上の流体のうちの１の流体の占有率αを算出するための算出手段と、
変数βおよびγを含む式（Ａ１）および式（Ａ２）を用いるＴＨＩＮＣ（tangent  of
hyperbola  for  interface  capturing）法に従って、各前記セルの第１の軸に沿った界面χ（ｘ）の再構築を行なうための再構築手段と、
【数１】
0005875409.tif 000017
前記式（Ａ１）における前記変数βを求めるための定数であって、前記シミュレーションにおいて定義される前記第１の軸を含む２以上の軸ごとに設定された定数βｕ，βｌを記憶するための記憶手段とを備え、
前記再構築手段は、前記軸ごとに設定された前記定数βｕ，βｌと、前記２種類以上の流体のうちの１の流体の速度ベクトルにおける前記軸ごとの寄与度ωとに基づいて、式（Ａ３）に従って前記変数βを設定し、
【数２】
0005875409.tif 000018
前記定数βｕは特定の値よりも大きい値であり、前記定数βｌは前記特定の値よりも小さい値である、シミュレーション装置。
【請求項２】
前記特定の値は、３．５である、請求項１に記載のシミュレーション装置。
【請求項３】
２種類以上の流体からなる混相流における自由界面をシミュレーションするためのコンピュータによって実行されるシミュレーション方法であって、
シミュレーション対象が分割された複数のセルのそれぞれにおける前記２種類以上の流体のうちの１の流体の占有率αを算出するステップと、
変数βおよびγを含む式（Ａ１）および式（Ａ２）を用いるＴＨＩＮＣ（tangent  of  hyperbola  for  interface  capturing）法に従って、各前記セルの第１の軸に沿った界面の再構築を行なうステップとを備え、
【数３】
0005875409.tif 000019
前記再構築を行なうステップは、前記式（Ａ１）における前記変数βを求めるための定数であって、前記シミュレーションにおいて定義される前記第１の軸を含む２以上の軸ごとに設定された定数βｕ，βｌと、前記２種類以上の流体のうちの１の流体の速度ベクトルにおける前記軸ごとの寄与度ωとに基づいて、式（Ａ３）に従って前記変数βを設定し、
【数４】
0005875409.tif 000020
前記定数βｕは特定の値よりも大きい値であり、前記定数βｌは前記特定の値よりも小さい値である、シミュレーション方法。
【請求項４】
２種類以上の流体からなる混相流における自由界面をシミュレーションするためのコンピュータによって実行されるシミュレーションプログラムであって、
前記シミュレーションプログラムは、前記コンピュータに、
シミュレーション対象が分割された複数のセルのそれぞれにおける前記２種類以上の流体のうちの１の流体の占有率αを算出するステップと、
変数βおよびγを含む式（Ａ１）および式（Ａ２）を用いるＴＨＩＮＣ（tangent  of  hyperbola  for  interface  capturing）法に従って、各前記セルの第１の軸に沿った界面の再構築を行なうステップとを備え、
【数５】
0005875409.tif 000021
前記再構築を行なうステップは、前記式（Ａ１）における前記変数βを求めるための定数であって、前記シミュレーションにおいて定義される前記第１の軸を含む２以上の軸ごとに設定された定数βｕ，βｌと、前記２種類以上の流体のうちの１の流体の速度ベクトルにおける前記軸ごとの寄与度ωとに基づいて、式（Ａ３）に従って前記変数βを設定し、
【数６】
0005875409.tif 000022
前記定数βｕは特定の値よりも大きい値であり、前記定数βｌは前記特定の値よりも小さい値である、シミュレーションプログラム。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、シミュレーション装置に関し、特に、混相流における自由界面を体積追跡法に従ってシミュレーションするための、シミュレーション装置、シミュレーション方法、および、シミュレーションプログラムに関する。
【背景技術】
【０００２】
従来から、混相流における自由界面の解析法について、種々の提案がなされている。その中の一つである体積追跡法は、体積占有率を用いて格子上に界面を表現する方法である。
【０００３】
体積占有率は、次の式（１）で示される。二相流系において、ある格子が完全に一方の流体で占められている場合は当該関数の値は１であり、完全に他方の流体で占められている場合は当該関数の値は０である。そして、この関数の値は、一方の流体の占める割合に応じて、０から１の値をとる。
【０００４】
【数１】
0005875409.tif 000002
【０００５】
式（１）において、ｔは時間を示し、αは各セルの一方の流体の占有率を示す。また、ｖベクトル（式（１）中のｖの上に矢印が付されて表記されたもの）は、上記一方の流体の、各セルの図２２に示されるような流速ｕの、再構築の方向についてのベクトル量である。
【０００６】
各セルの占有率について、図２３および図２４を参照して説明する。
図２３は、２つの流体が存在する状態を模式的に示す図である。図２３では、ハッチングを付されていない部分が一方の流体の存在領域を示し、ハッチングを付された部分が他方の流体の存在領域を示す。なお、図２３では、６×６の３６個のセルが定義された状態が示されている。そして、図２３の３６個のセルのそれぞれにおける、上記一方の流体の占有率を、図２４に示す。
【０００７】
図２４に示されるような各セルの或る流体の占有率から、図２３に示されるような界面を再構築するための方法も、従来から種々提案されている。
【０００８】
たとえば、非特許文献１には、ＴＨＩＮＣ（tangent  of  hyperbola  for  interface  capturing）法が開示されている。この方法は、１）移流量の正確な保存、２）数値誤差の効率的な除去、および、３）体積占有率における階段的な変化の近傍でのスプリアスな発振（spurious  oscillation）を発生の回避を目的として、特徴的な補完関数を採用している。この方法では、補完関数として、双曲線正接関数（the  hyperbolic  tangent  functions）が利用されている。
【０００９】
さらに、非特許文献２には、ｘ軸とｙ軸の両方の界面の単位法線ベクトルを重み付けして再構築に用いる、ＴＨＩＮＣ／ＷＬＩＣ（weighted  line  interface  calculation）法が開示されている。
【先行技術文献】
【非特許文献】
【００１０】
【非特許文献１】F. Xiao, Y. Honma, T. Kono, "A simple algebraic interface capturing scheme using hyperbolic tangent function", International Journal for Numerical Methods in Fluid 48 (2005) 1023-1040
【非特許文献２】K. Yokoi, "Efficient implementation of THINC scheme: A simple and practical smoothed VOF algorithm", Journal of Computational Physics 226 (2007) 1985-2002
【発明の概要】
【発明が解決しようとする課題】
【００１１】
しかしながら、複雑な自由界面の流れのシミュレーションは、まだ現実の流れを正確に再現できているとは言えず、検討の余地が残されている。
【００１２】
本発明は、かかる実情に鑑み考え出されたものであり、その目的は、より正確に自由界面を再構築できるシミュレーション装置、シミュレーション方法、および、シミュレーションプログラムを提供することである。
【課題を解決するための手段】
【００１３】
本発明に従ったシミュレーション装置は、２種類以上の流体からなる混相流における自由界面をシミュレーションするためのシミュレーション装置であって、シミュレーション対象が分割された複数のセルのそれぞれにおける２種類以上の流体のうちの１の流体の占有率αを算出するための算出手段と、変数βおよびγを含む式（Ａ１）および式（Ａ２）を用いるＴＨＩＮＣ（tangent  of  hyperbola  for  interface  capturing）法に従って、各セルの第１の軸に沿った界面χ（ｘ）の再構築を行なうための再構築手段と、
【００１４】
【数２】
0005875409.tif 000003
【００１５】
式（Ａ１）における変数βを求めるための定数であって、シミュレーションにおいて定義される第１の軸を含む２以上の軸ごとに設定された定数を記憶するための記憶手段とを備える。再構築手段は、軸ごとに設定された定数と、２種類以上の流体のうちの１の流体の速度ベクトルにおける軸ごとの寄与度とに基づいて変数βを設定し、軸ごとに設定された定数のうち、１の定数は特定の値よりも大きい値であり、他の定数は特定の値よりも小さい値である。
【００１６】
好ましくは、特定の値は、３．５である。
本発明に従ったシミュレーション方法は、２種類以上の流体からなる混相流における自由界面をシミュレーションするためのコンピュータによって実行されるシミュレーション方法であって、シミュレーション対象が分割された複数のセルのそれぞれにおける２種類以上の流体のうちの１の流体の占有率αを算出するステップと、変数βおよびγを含む式（Ａ１）および式（Ａ２）を用いるＴＨＩＮＣ（tangent  of  hyperbola  for  interface  capturing）法に従って、各セルの第１の軸に沿った界面の再構築を行なうステップと、
【００１７】
【数３】
0005875409.tif 000004
【００１８】
式（Ａ１）における変数βを求めるための定数であって、シミュレーションにおいて定義される第１の軸を含む２以上の軸ごとに設定された定数を記憶するステップとを備える。再構築を行なうステップは、軸ごとに設定された定数と、２種類以上の流体のうちの１の流体の速度ベクトルにおける軸ごとの寄与度とに基づいて変数βを設定し、軸ごとに設定された定数のうち、１の定数は特定の値よりも大きい値であり、他の定数は特定の値よりも小さい値である。
【００１９】
本発明に従ったシミュレーションプログラムは、２種類以上の流体からなる混相流における自由界面をシミュレーションするためのコンピュータによって実行されるシミュレーションプログラムであって、シミュレーションプログラムは、コンピュータに、シミュレーション対象が分割された複数のセルのそれぞれにおける２種類以上の流体のうちの１の流体の占有率αを算出するステップと、変数βおよびγを含む式（Ａ１）および式（Ａ２）を用いるＴＨＩＮＣ（tangent  of  hyperbola  for  interface  capturing）法に従って、各セルの第１の軸に沿った界面の再構築を行なうステップと、
【００２０】
【数４】
0005875409.tif 000005
【００２１】
式（Ａ１）における変数βを求めるための定数であって、シミュレーションにおいて定義される第１の軸を含む２以上の軸ごとに設定された定数を記憶するステップとを実行させる。再構築を行なうステップは、軸ごとに設定された定数と、２種類以上の流体のうちの１の流体の速度ベクトルにおける軸ごとの寄与度とに基づいて変数βを設定し、軸ごとに設定された定数のうち、１の定数は特定の値よりも大きい値であり、他の定数は特定の値よりも小さい値である。
【発明の効果】
【００２２】
本発明によれば、複数の軸を対象とする場での速度ベクトルにおける、再構築の軸方向の成分の寄与度が、より大きく反映されることになる。
【図面の簡単な説明】
【００２３】
【図１】この発明の実施の形態に従うシミュレーション装置を実現するための代表的なハードウェア構成であるコンピュータの概略構成図である。
【図２】ＴＨＩＮＣ法による、一次元の占有率の変化を示す図である。
【図３】実施の形態における変数の算出方法を説明するための図である。
【図４】ＴＨＩＮＣ／ＷＬＩＣ法の概略を示す図である。
【図５】本実施の形態の計算方法の概略を説明するための図である。
【図６】本実施の形態のシミュレーションにおける処理のフローチャートである。
【図７】図６の処理内容を説明するための図である。
【図８】図６の処理内容を説明するための図である。
【図９】図６の処理内容を説明するための図である。
【図１０】図６の処理内容を説明するための図である。
【図１１】図６の処理内容を説明するための図である。
【図１２】図１０に対応した各セルの密度を示す図である。
【図１３】図１１に対応した各セルの密度を示す図である。
【図１４】比較例の結果を示す図である。
【図１５】比較例の結果を示す図である。
【図１６】比較例の結果を示す図である。
【図１７】比較例の結果を示す図である。
【図１８】実施の形態の結果を示す図である。
【図１９】実施の形態の結果の初期状態を示す図である。
【図２０】比較例の結果を示す図である。
【図２１】実施の形態の結果を示す図である。
【図２２】従来の技術を説明するための図である。
【図２３】従来の技術を説明するための図である。
【図２４】従来の技術を説明するための図である。
【発明を実施するための形態】
【００２４】
本発明の実施の形態について、図面を参照しながら詳細に説明する。なお、図中の同一または相当部分については、同一符号を付してその説明は繰返さない。
【００２５】
＜ハードウェア構成＞
図１は、この発明の実施の形態に従うシミュレーション装置を実現するための代表的なハードウェア構成であるコンピュータ１００の概略構成図である。
【００２６】
図１を参照して、コンピュータ１００は、ＦＤ（Flexible  Disk）駆動装置１１１およびＣＤ−ＲＯＭ（Compact  Disk-Read  Only  Memory)駆動装置１１３を搭載したコンピュータ本体１０１と、モニタ１０２と、キーボード１０３と、マウス１０４とを含む。
【００２７】
コンピュータ本体１０１は、ＦＤ駆動装置１１１およびＣＤ−ＲＯＭ駆動装置１１３に加えて、相互にバスで接続された、演算装置であるＣＰＵ（Central  Processing  Unit
）１０５と、メモリ１０６と、記憶装置である固定ディスク１０７と、通信インターフェース１０９とを含む。
【００２８】
本実施の形態に従うシミュレーションは、ＣＰＵ１０５がメモリ１０６などのコンピュータハードウェアを用いて、プログラムを実行することで実現される。一般的に、このようなプログラムは、ＦＤ１１２やＣＤ−ＲＯＭ１１４などの記録媒体に格納されて、またはネットワークなどを介して流通する。そして、このようなプログラムは、ＦＤ駆動装置１１１やＣＤ−ＲＯＭ駆動装置１１３などにより記録媒体から読取られて、または通信インターフェース１０９にて受信されて、固定ディスク１０７に格納される。さらに、このようなプログラムは、固定ディスク１０７からメモリ１０６に読出されて、ＣＰＵ１０５により実行される。
【００２９】
ＣＰＵ１０５は、様々な数値論理演算を行なう演算処理部であり、プログラムされた命令を順次実行することで、本実施の形態に従う制御機能を提供する。メモリ１０６は、ＣＰＵ１０５のプログラム実行に応じて各種の情報を記憶する。
【００３０】
モニタ１０２は、ＣＰＵ１０５が出力する情報を表示するための表示部であって、一例としてＬＣＤ（Liquid  Crystal  Display）やＣＲＴ（Cathode  Ray  Tube）などから構成される。すなわち、モニタ１０２には、シミュレーション結果などが表示される。
【００３１】
マウス１０４は、クリックやスライドなどの動作に応じたユーザから指令を受付ける。キーボード１０３は、入力されるキーに応じたユーザから指令を受付ける。
【００３２】
通信インターフェース１０９は、コンピュータ１００と他の装置との間の通信を確立するための装置であり、各種データを外部から受付可能である。
【００３３】
なお、上述したようなコンピュータに代えて、その一部または全部を専用のハードウェアによって実現してもよい。
【００３４】
また、コンピュータ１００を本実施の形態に従うように機能させるプログラムが記録される記録媒体は、上記したＦＤやＣＤ−ＲＯＭに限定されない。記録媒体の他の例としては、ＤＶＤ−ＲＯＭ（Digital  Versatile  Disk  -  Read  Only  Memory）、ＵＳＢ（Universal  Serial  Bus）メモリ、メモリカード、ハードディスク、磁気テープ、カセットテープ、ＭＯ（Magnetic  Optical  Disc）、ＭＤ（Mini  Disc）、ＩＣ（Integrated  Circuit）カード（メモリカードを除く）、光カード、マスクＲＯＭ、ＥＰＲＯＭ、ＥＥＰＲＯＭ（Electronically  Erasable  Programmable  Read-Only  Memory）などの、不揮発的にプログラムを格納する媒体が挙げられる。
【００３５】
＜本実施の形態のシミュレーション方法の概要＞
上記したように、自由界面の再構築に関し、従来、ＴＨＩＮＣ法およびＴＨＩＮＣ／ＷＬＩＣ法が提案されている。一方、本実施の形態のシミュレーション方法は、ＴＨＩＮＣ／ＷＬＩＣ法において値を固定されていた定数（β）を、軸ごとに設定された定数（後述するβｕ，βｌ）を用い、さらに、占有率を考慮する流体の速度ベクトルにおける軸ごとの寄与度（後述するω）に基づいて設定する。
【００３６】
このことから、本明細書では、まず、本実施の形態のシミュレーション方法の前提となるＴＨＩＮＣ法およびＴＨＩＮＣ／ＷＬＩＣ法について説明した上で、本実施の形態のシミュレーション方法を説明する。
【００３７】
なお、本明細書では、シミュレーションにおいて定義される軸としてｘ軸とｙ軸の２軸を想定して記載がなされる部分も含む。つまり、２次元でのシミュレーションが説明される部分を含む。ただし、本実施の形態において開示されるシミュレーション方法は、公知の技術に従って３次元に拡張可能である。なお、３次元への拡張の方法については、公知の技術を採用可能であるため、ここでは説明は繰り返さない。
【００３８】
＜ＴＨＩＮＣ法＞
図２は、ＴＨＩＮＣ法による、一次元の占有率の変化を示す図である。ＴＨＩＮＣ法は、図２に示されるように、平滑化された双曲線正接関数を利用している。なお、図２の破線は、体積追跡法による変化を示している。
【００３９】
ＴＨＩＮＣ法では、特殊関数（characteristic  function）として、次の式（２）が使用される。
【００４０】
【数５】
0005875409.tif 000006
【００４１】
式（２）において、γは次の式（３）で表されるように、ｎｘ，ｉの値に応じて「１」または「−１」の値をとる。
【００４２】
【数６】
0005875409.tif 000007
【００４３】
式（３）のｎｘ，ｉは、法線ベクトルのｘ成分であり、計算対象とするセルを３×３の格子状の中心に置いた場合の、残りの９個のセルの占有率αに基づいて算出される。具体的には、図３に示すように３×３の格子状の９個のセルを想定する。各セルの占有率は、上段では、左から、順に、αｉ−１，ｊ＋１、αｉ，ｊ＋１、αｉ＋１，ｊ＋１である。中段では、左から、順に、αｉ−１，１、αｉ，ｊ、αｉ＋１，ｊである。下段では、左から、順に、αｉ−１，ｊ−１、αｉ，ｊ−１、αｉ＋１，ｊ−１である。そして、これらの占有率に基づき、次の式（４−１）〜式（４−４）および式（５）から、ｎｘ，ｉが算出される。
【００４４】
【数７】
0005875409.tif 000008
【００４５】
【数８】
0005875409.tif 000009
【００４６】
式（２）に戻って、βは、図２に示された双曲線正接関数の傾きを定めるパラメータである。一般的には、βとして、３．５という固定された値が採用される。
【００４７】
式（２）中の「ｘｉ〜」は、式（２）のχ（ｘ）の平均を用い、次の式（６）に従って算出される。なお、本明細書では、チルダ記号が「〜」で表現される。
【００４８】
【数９】
0005875409.tif 000010
【００４９】
つまり、上記のように設定されたγ，βを用いて式（２）のχ（ｘ）が構築され、そして、式（６）に基づいて、「ｘｉ〜」が算出される。
【００５０】
＜ＴＨＩＮＣ／ＷＬＩＣ法＞
ＴＨＩＮＣ／ＷＬＩＣ法では、ｘ方向の流束は、式（７）に従って算出される。なお、式（７）中の「ＲＦ」はRight Faceを意味する。つまり、右端の面（後述する図２２のｕで示された黒丸が位置する縦線が存在する面）を意味し、より具体的にはセル番号「ｉ＋１／２」を意味する。
【００５１】
【数１０】
0005875409.tif 000011
【００５２】
式（７）の右辺において、第１項はＴＨＩＮＣ法に従った流束であり、第２項はＵＰＷＩＮＤ法（一次風上法）に従った流束である。
【００５３】
式（７）において、ωｘ，ｉは、次の式（８）に従って算出される。
【００５４】
【数１１】
0005875409.tif 000012
【００５５】
式（８）において、ｎｘ，ｉは、上記のように、式（４−１）〜式（４−４）および式（５）から算出される。また、ｎｙ，ｉは、図３のような９個のセルの占有率に基づき、次の式（９−１）〜式（９−４）および式（１０）から算出される。
【００５６】
【数１２】
0005875409.tif 000013
【００５７】
【数１３】
0005875409.tif 000014
【００５８】
つまり、ＴＨＩＮＣ／ＷＬＩＣ法では、図４に示されるように、ｎｙ，ｉが０の場合には、式（８）に従ってωｘ，ｉが１になり、これにより、ｘ方向の流束は式（７）の右辺の第１項のみに基づいて算出される。また、ｎｘ，ｉが０の場合には、式（８）に従ってωｘ，ｉが０になり、これにより、ｘ方向の流束は式（７）の右辺の第２項のみに基づいて算出される。
【００５９】
そして、これらの場合の間の状態では、式（８）において、速度ベクトルのｘ成分とｙ成分の割合に応じてωｘ，ｉが算出される。これにより、ｘ方向の流束は、算出されたωｘ，ｉの値に基づき、式（７）の右辺の第１項と第２項が式（７）に従って組み合わされて、算出される。なお、式（７）の右辺の第１項は、ＴＨＩＮＣ法による流束の算出に対応し、第２項は、ＵＰＷＩＮＤ法による流束の算出に対応する。
【００６０】
＜本実施の形態の計算方法＞
本実施の形態の計算方法では、基本的には、図４等を参照して説明したＴＨＩＮＣ／ＷＬＩＣ法に従った計算方法で、各セルの再構築を行なう。
【００６１】
図５は、本実施の形態の計算方法の概略を説明するための図である。
図５を参照して、本実施の形態の計算方法は、式（２）で示されたＴＨＩＮＣ法の特殊関数において、変数βの設定方法に特徴を有する。
【００６２】
具体的には、変数βを、図５中に示された式（１１）に示されるように、βｕとβｌとωｘ，ｉとを用いて算出して、設定する。βｕとβｌの値は、たとえば固定ディスク１０７に格納されている。
【００６３】
βｕは、従来βとして利用されてきた値３．５よりも界面付近を急激に変化させる値であり、たとえば１０とされる。なお、図５のβｕは、式（２）においてγ≧０の場合の結果が示されている。βｌは、従来βとして利用されてきた値３．５よりも界面付近の変化を緩くさせる値であり、たとえば０．００１とされる。図５には、βｕとβｌのそれぞれが式（２）中のβとして利用された場合の、界面付近の占有率の変化が示されている。
【００６４】
本実施の形態によれば、βは、ｎｘ，ｉとｎｙ，ｉの大小関係に応じて設定される。
具体的には、ｎｘ，ｉがｎｙ，ｉに比べて大きいときには、つまり、セルの速度ベクトルのｘ成分がｙ成分よりも大きい場合には、当該セルの再構築において、従来よりも、ｘ方向の界面の流体の占有率の変化をより急激にするように、βが設定される。式（１１）に従ったβの設定において、βｕの方がβｌよりも大きい割合で反映されるからである。
【００６５】
一方、ｎｘ，ｉがｎｙ，ｉに比べて小さいときには、つまり、セルの速度ベクトルのｘ成分がｙ成分よりも小さい場合には、当該セルの再構築において、従来よりも、ｘ方向の界面の流体の占有率の変化がより緩やかになるように、βが設定される。式（１１）に従ったβの設定において、βｌの方がβｕよりも大きい割合で反映されるからである。
【００６６】
つまり、本実施の形態に従った界面の再構築では、複数の軸を対象とする場での速度ベクトルにおける、再構築の軸方向の成分の寄与度（本実施の形態では、ｎｘ，ｉとｎｙ，ｉの和に対するｎｘ，ｉの寄与度）が、より大きく反映されることになる。
【００６７】
＜本実施の形態における計算フロー＞
図６は、本実施の形態のシミュレーションにおける処理のフローチャートである。
【００６８】
以下の説明では、シミュレーションの対象として、液体中の気体の移動を例示する。これにより、占有率αとして、各セルにおける気体の占有率が例示される。なお、本実施の形態のシミュレーション対象となる混相流における相の組合せはこれに限定されるものではない。異なる２以上の流体が存在する系であれば、異なる２以上の流体または気体が存在する系であっても良い。
【００６９】
図６を参照して、コンピュータ１００に対して計算開始の指示が入力されると、ＣＰＵ１０５は、まずステップＳ１０において、各セルの占有率αと速度ｖの初期値を読み込んで、ステップＳ２０へ処理を進める。なお、これらの値は、たとえば固定ディスク１０７に格納されている。
【００７０】
ステップＳ２０では、ＣＰＵ１０５は、各セルの速度ｖならびに当該速度のｘ方向およびｙ方向の法線ベクトルを算出して、ステップＳ３０へ処理を進める。
【００７１】
ステップＳ３０では、ＣＰＵ１０５は、最新の占有率αについての「∇α」を算出して、ステップＳ４０へ処理を進める。
【００７２】
ステップＳ４０では、各セルの占有率αと速度ｖの積についての「∇αｖ」を算出して、ステップＳ５０へ処理を進める。
【００７３】
ステップＳ５０では、各セルの占有率αに対してΔｔ後の占有率αを算出して、ステップＳ６０へ処理を進める。
【００７４】
ステップＳ６０では、最新の占有率αに基づいて、式（３）に従ってγを算出し、式（１１）に従ってβを算出し、式（２）および式（６）に従って「ｘｉ〜」を算出することにより、各セルの界面を再構築して、ステップＳ７０へ処理を進める。
【００７５】
ステップＳ７０では、ＣＰＵ１０５は、シミュレーションの期間（Δｔの総和）Ｔsumが指定された期間Ｔを超えたか否かを判断し、まだ超えていなければステップＳ２０へ処理が戻され、超えていれば処理を終了させる。
【００７６】
以上、図６を参照して説明した処理により、微小時間Δｔごとに、各セルの占有率が更新され、そして、更新後の占有率に従って各セルの界面が、図５等を参照して説明された本実施の形態の計算方法に従って、再構築される。
【００７７】
なお、本実施の形態のシミュレーションでは、シミュレーションの対象の領域として、図７に示されるようなセルの格子（図７では１００×１００のセルの格子）が定義される。そして、ステップＳ２０からステップＳ６０の処理が一回行なわれるごとに、各セルの占有率αの計算結果が更新され、そして、各セルにおける界面の再構築が行なわれる。
【００７８】
図８および図９は、図７に示された格子状のセルの一部（１０×１０のセルの格子）についての、各セルの占有率αの計算結果を示す。なお、図９は、図８の計算結果に対する、時間ΔＴ後の計算結果を示すものとする。
【００７９】
また、図１０，図１１は、それぞれ、図８，図９の計算結果に対応した、界面の再構築の状態を示す図である。つまり、図１１は、図１０に示された結果に対する、時間ΔＴ後の結果を示す。
【００８０】
図８では、３〜６列目および３〜６行目にある１６個のセル以外のセルでは、占有率αは「０．０」とされている。一方、これらの１６個のセルでは、占有率αとして、０．４，０．８，１．０のいずれかの値が算出されている。そして、各占有率αに対応して、図１０に示されるように界面が再構築されている。
【００８１】
図９では、図８に対して、占有率αが０．０以外の値を有する１６個のセルが、それぞれ１列ずつ右方向に移動している。図１１では、これに対応して、再構築される界面が１列ずつ、図１０に示されたものに対して右に移動している。
【００８２】
また、図１０および図１１に示されたような各セルの占有率αに基づき、次の式（１２）に従って、各セルの密度が算出される。図１０に対応した密度の算出結果を図１２に、図１１に対応した密度の算出結果を図１３に、それぞれ示す。なお、式（１２）において、ρGは気体の密度であり、ρｌは液体の密度である。また、図１２および図１３では、ρGの値の一例として１kg/m3が、ρｌの値の一例として１０００kg/m3が、それぞれ採用されている。
【００８３】
【数１４】
0005875409.tif 000015
【００８４】
＜計算結果の比較＞
次に、本実施の形態による界面の再構築の結果を、これまでの種々の方法による再構築の結果と比較する。
【００８５】
（５回の回転）
図１４〜図１８は、液相中の気泡を、各方法によって５回回転移動させるシミュレーションをした際の界面の再構築の結果を示す図である。なお、図１４〜図１８の各図では、図１９に示されるように、気泡ＢＵで示される初期状態の気泡が矢印Ｒで示される回転方向に５回転された際の結果が示されている。各図のシミュレーション結果では、気泡の占有率が色の濃淡によって５段階（０．０，０．２５，０．５，０．７５，１．０）で示されている。
【００８６】
なお、図１４はＵＰＷＩＮＤ法（一次風上法）による結果を、図１５はＤＡ法（Ｄｏｎｏｒ−Ａｃｃｅｐｔｏｒ法）による結果を、図１６はＴＨＩＮＣ法による結果を、図１７はＴＨＩＮＣ／ＷＬＩＣ法による結果を、そして、図１８は本実施の形態による結果を、それぞれ示している。
【００８７】
なお、各図には、計算結果に関する数値が記載されている。上段は、界面によって規定される気泡の総面積の、初期形状に対する誤差（ＥＳＺ）の百分率を示す。下段は、界面によって規定される形状についての、初期形状に対する誤差（ＥＳＰ）の百分率を示す。各誤差は、次の式（１３）および式（１４）に従って、算出される。
【００８８】
【数１５】
0005875409.tif 000016
【００８９】
総面積については、図１６〜図１８に示された各結果は、図１４に示された結果に対して、１０９倍程度の精度を有し、また、図１５に示された結果に対しても、１０７倍程度の精度を有している。このことは、図１５の結果において界面で囲まれた領域の大きさが図１９の気泡ＢＵに対して大きくなっていたり、図１６の結果において回転軌道上に界面が残っていたりするのに対し、図１６〜図１８の各結果では、界面で囲まれた領域の大きさが図１９の気泡ＢＵとほぼ同様であることに対応する。
【００９０】
また、図１６〜図１８に示された各結果は、形状についても、図１４または図１５に示された各結果に対して、格段の精度を有している。
【００９１】
図１７のＴＨＩＮＣ／ＷＬＩＣ法は、図１６のＴＨＩＮＣ法に対して、総面積についての精度は多少劣るものの、形状については格段の精度を有している。
【００９２】
そして、図１８の本実施の形態は、ＴＨＩＮＣ／ＷＬＩＣ法は、総面積および形状のいずれについても、上回った精度を有し、また、ＴＨＩＮＣ法に対しても、総面積については多少精度は劣るものの、形状については格段の精度を有している。
【００９３】
（１８００回の回転）
図２０および図２１は、上記のような気泡の回転のシミュレーションとして、１８００回転された際の結果を示している。図２０は、ＴＨＩＮＣ／ＷＬＩＣ法の結果であり、図２１は、本実施の形態の結果である。また、図２０と図２１のそれぞれには、図１９に示されたような初期状態に対する形状の誤差の百分率が示されている。
【００９４】
図２０と図２１に示された界面の形状を比較すると、図２０では、初期状態の円形に対して横長の形状となているのに対し、図２１では、初期状態の円形が保持されている。このことは、形状について、本実施の形態の方がＴＨＩＮＣ／ＷＬＩＣ法よりも高い精度を達成できていることを示している。
【００９５】
＜その他の変形例等＞
以上説明した本実施の形態では、ステップＳ５０を実行するＣＰＵ１０５によって、各セルの占有率を算出するための算出手段が構成されている。そして、ステップＳ６０を実行するＣＰＵ１０５によって、各セルの界面の再構築を行なう再構築手段が構成されている。
【００９６】
以上説明した本実施の形態におけるシミュレーショ方法では、界面の再構築に利用される変数βの設定方法に特徴を有する。具体的には、図５を参照して説明したように、βの算出に、βｕとβｌが利用される。そして、本実施の形態では、βｕの一例として「１０」が採用され、また、βｌの一例として「０．００１」が採用された。なお、本実施の形態において採用されるβｕとβｌの値の組合せは、これに限定されない。βｕが、従来の技術においてβとして採用されていた値よりも大きく、そして、βｌが当該値よりも小さければ良い。
【００９７】
従来βとして利用されてきた値としては、「３．５」が挙げられる。そして、本実施の形態では、特に、βｕは、界面の再構築の精度向上の観点から、５．０以上であることが好ましく、１０．０以上であることがさらに好ましく、１５．０以上であることがさらに好ましい。
【００９８】
今回開示された実施の形態はすべての点で例示であって制限的なものではないと考えられるべきである。本発明の範囲は上記した説明ではなく、特許請求の範囲によって示され、特許請求の範囲と均等の意味および範囲内でのすべての変更が含まれることが意図される。また、実施の形態やその変形例に示された内容は、それぞれ独立して、または、可能な限り組み合わされて、実施されることが意図される。
【符号の説明】
【００９９】
１００  コンピュータ、１０１  コンピュータ本体、１０２  モニタ、１０３  キーボード、１０４  マウス、１０６  メモリ、１０７  固定ディスク、１０９  通信インターフェース、１１１  ＦＤ駆動装置、１１３  ＣＤ−ＲＯＭ駆動装置。
【図１】
0005875409.tif 000023
【図２】
0005875409.tif 000024
【図３】
0005875409.tif 000025
【図４】
0005875409.tif 000026
【図５】
0005875409.tif 000027
【図６】
0005875409.tif 000028
【図７】
0005875409.tif 000029
【図８】
0005875409.tif 000030
【図９】
0005875409.tif 000031
【図１０】
0005875409.tif 000032
【図１１】
0005875409.tif 000033
【図１２】
0005875409.tif 000034
【図１３】
0005875409.tif 000035
【図２２】
0005875409.tif 000036
【図２３】
0005875409.tif 000037
【図２４】
0005875409.tif 000038
【図１４】
0005875409.tif 000039
【図１５】
0005875409.tif 000040
【図１６】
0005875409.tif 000041
【図１７】
0005875409.tif 000042
【図１８】
0005875409.tif 000043
【図１９】
0005875409.tif 000044
【図２０】
0005875409.tif 000045
【図２１】
0005875409.tif 000046