特許5945183 | 知財ポータル「IP Force」

知財求人 - 知財ポータルサイト「IP Force」

▶ 住友化学株式会社の特許一覧

特許5945183シミュレーション装置、シミュレーション方法、およびシミュレーションプログラム

書誌要約請求の範囲詳細な説明課題実施例実施するための形態図面の説明

目に優しい文字サイズ小中大 PDF Top

< >

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30

(19)【発行国】日本国特許庁(JP)

(12)【公報種別】特許公報(B2)

(11)【特許番号】5945183

(24)【登録日】2016年6月3日

(45)【発行日】2016年7月5日

(54)【発明の名称】シミュレーション装置、シミュレーション方法、およびシミュレーションプログラム

(51)【国際特許分類】

G06F 19/00 20110101AFI20160621BHJP

【ＦＩ】

G06F19/00 110

【請求項の数】6

【全頁数】33

(21)【出願番号】特願2012-164698(P2012-164698)

(22)【出願日】2012年7月25日

(65)【公開番号】特開2014-26379(P2014-26379A)

(43)【公開日】2014年2月6日

【審査請求日】2015年5月22日

(73)【特許権者】

【識別番号】000002093

【氏名又は名称】住友化学株式会社

(74)【代理人】

【識別番号】110001195

【氏名又は名称】特許業務法人深見特許事務所

(72)【発明者】

【氏名】ダールアビナーブ

(72)【発明者】

【氏名】島田直樹

【審査官】宮地匡人

(56)【参考文献】

【文献】 Xiao F，A simple algebraic interface capturing scheme using hyperbolic tangent function，Int. J. Numer. Meth. Fluids，２００５年４月７日，No.48，pp.1023-1040

【文献】伊井仁志，THINC法を発展させた多次元連続関数を用いる界面捕獲手法の提案，第23回計算力学会講演会CD-ROM論文集，日本機械学会，２０１０年９月２５日，No.10-2，pp.29-30

【文献】李光浩，新たな代数的な気液混相流の境界追跡法の提案，土木学会中部支部研究発表会講演概要集，２０１２年３月，pp.137-138

(58)【調査した分野】（Int.Cl.，ＤＢ名）

Ｇ０６Ｆ１９／００

ＣｉＮｉｉ

ＪＳＴＰｌｕｓ（ＪＤｒｅａｍＩＩＩ）

(57)【特許請求の範囲】

【請求項1】

２種類以上の流体からなる混相流における自由界面をシミュレーションするためのシミュレーション装置であって、
シミュレーション対象が分割された複数のセルのそれぞれにおける前記２種類以上の流体のうちの１の流体の占有率αを算出するための算出手段と、
変数βおよびγを含む式（Ａ１）および式（Ａ２）を用いるＴＨＩＮＣ（tangent of hyperbola for interface capturing）法に従って、各前記セルの第１の軸に沿った界面χ（ｘ）の再構築を行なうための再構築手段とを備え、

【数1】

前記再構築手段は、各前記セルの前記第１の軸に沿った方向についての前記変数βを、式（Ａ３）を用いて算出する、シミュレーション装置。

【数2】

（式（Ａ３）において、文字ｋは前記第１の軸を表し、文字ｎ_ｋは前記第１の軸方向の法線ベクトルを表し、文字ｊはシミュレーションにおいて定義される直行座標軸の前記第１の軸以外の軸を表し、文字ｎ_ｊは前記第１の軸以外の軸の法線ベクトルを表す。）

【請求項2】

【数3】

前記再構築手段は、各前記セルの前記第１の軸に沿った方向についての前記変数βを、式（Ａ４）または式（Ａ５）を用いて算出する、シミュレーション装置。

【数4】

（式（Ａ４）および式（Ａ５）において、文字ｎ_ｘはｘ軸方向の法線ベクトルを表し、文字ｎ_ｙはｙ軸方向の法線ベクトルを表し、文字ｎ_ｚはｚ軸方向の法線ベクトルを表す。）

【請求項3】

２種類以上の流体からなる混相流における自由界面をシミュレーションするためのコンピュータによって実行されるシミュレーション方法であって、
シミュレーション対象が分割された複数のセルのそれぞれにおける前記２種類以上の流体のうちの１の流体の占有率αを算出するステップと、
変数βおよびγを含む式（Ａ１）および式（Ａ２）を用いるＴＨＩＮＣ（tangent of hyperbola for interface capturing）法に従って、各前記セルの第１の軸に沿った界面の再構築を行なうステップと、

【数5】

前記式（Ａ１）における前記変数βを求めるための定数であって、前記シミュレーションにおいて定義される前記第１の軸を含む２以上の軸ごとに設定された定数を記憶するステップとを備え、
前記再構築を行なうステップは、各前記セルの前記第１の軸に沿った方向についての前記変数βを、式（Ａ３）を用いて算出する、シミュレーション方法。

【数6】

【請求項4】

【数7】

前記式（Ａ１）における前記変数βを求めるための定数であって、前記シミュレーションにおいて定義される前記第１の軸を含む２以上の軸ごとに設定された定数を記憶するステップとを備え、
前記再構築を行なうステップは、各前記セルの前記第１の軸に沿った方向についての前記変数βを、式（Ａ４）または式（Ａ５）を用いて算出する、シミュレーション方法。

【数8】

【請求項5】

２種類以上の流体からなる混相流における自由界面をシミュレーションするためのコンピュータによって実行されるシミュレーションプログラムであって、
前記シミュレーションプログラムは、前記コンピュータに、
シミュレーション対象が分割された複数のセルのそれぞれにおける前記２種類以上の流体のうちの１の流体の占有率αを算出するステップと、
変数βおよびγを含む式（Ａ１）および式（Ａ２）を用いるＴＨＩＮＣ（tangent of hyperbola for interface capturing）法に従って、各前記セルの第１の軸に沿った界面の再構築を行なうステップと、

【数9】

前記式（Ａ１）における前記変数βを求めるための定数であって、前記シミュレーションにおいて定義される前記第１の軸を含む２以上の軸ごとに設定された定数を記憶するステップとを実行させ、
前記再構築を行なうステップは、各前記セルの前記第１の軸に沿った方向についての前記変数βを、式（Ａ３）を用いて算出する、シミュレーションプログラム。

【数10】

【請求項6】

【数11】

前記式（Ａ１）における前記変数βを求めるための定数であって、前記シミュレーションにおいて定義される前記第１の軸を含む２以上の軸ごとに設定された定数を記憶するステップとを実行させ、
前記再構築を行なうステップは、各前記セルの前記第１の軸に沿った方向についての前記変数βを、式（Ａ４）または式（Ａ５）を用いて算出する、シミュレーションプログラム。

【数12】

【発明の詳細な説明】

【技術分野】

【0001】

本発明は、シミュレーション装置に関し、特に、混相流における自由界面を体積追跡法に従ってシミュレーションするための、シミュレーション装置、シミュレーション方法、および、シミュレーションプログラムに関する。

【背景技術】

【0002】

従来から、混相流における自由界面の解析法について、種々の提案がなされている。その中の一つである体積追跡法は、体積占有率を用いて格子上に界面を表現する方法である。

【0003】

体積占有率は、次の式（１）で示される。二相流系において、ある格子が完全に一方の流体で占められている場合は当該関数の値は１であり、完全に他方の流体で占められている場合は当該関数の値は０である。そして、この関数の値は、一方の流体の占める割合に応じて、０から１の値をとる。

【0004】

【数1】

【0005】

式（１）において、ｔは時間を示し、αは各セルの一方の流体の占有率を示す。また、ｖベクトル（式（１）中のｖの上に矢印が付されて表記されたもの）は、上記一方の流体の、各セルの図２８に示されるような流速ｕの、再構築の方向についてのベクトル量である。

【0006】

各セルの占有率について、図２９および図３０を参照して説明する。
図２９は、２つの流体が存在する状態を模式的に示す図である。図２９では、ハッチングを付されていない部分が一方の流体の存在領域を示し、ハッチングを付された部分が他方の流体の存在領域を示す。なお、図２９では、６×６の３６個のセルが定義された状態が示されている。そして、図２９の３６個のセルのそれぞれにおける、上記一方の流体の占有率を、図３０に示す。

【0007】

図３０に示されるような各セルの或る流体の占有率から、図２９に示されるような界面を再構築するための方法も、従来から種々提案されている。

【0008】

たとえば、非特許文献１には、ＴＨＩＮＣ（tangent of hyperbola for interface capturing）法が開示されている。この方法は、１）移流量の正確な保存、２）数値誤差の効率的な除去、および、３）体積占有率における階段的な変化の近傍でのスプリアスな発振（spurious oscillation）を発生の回避を目的として、特徴的な補完関数を採用している。この方法では、補完関数として、双曲線正接関数（the hyperbolic tangent functions）が利用されている。

【0009】

さらに、非特許文献２には、ｘ軸とｙ軸の両方の界面の単位法線ベクトルを重み付けして再構築に用いる、ＴＨＩＮＣ／ＷＬＩＣ（weighted line interface calculation）法が開示されている。

【先行技術文献】

【非特許文献】

【0010】

【非特許文献1】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

【非特許文献2】K. Yokoi, "Efficient implementation of THINC scheme: A simple and practical smoothed VOF algorithm", Journal of Computational Physics 226 (2007) 1985-2002

【発明の概要】

【発明が解決しようとする課題】

【0011】

しかしながら、複雑な自由界面の流れのシミュレーションは、まだ現実の流れを正確に再現できているとは言えず、検討の余地が残されている。

【0012】

本発明は、かかる実情に鑑み考え出されたものであり、その目的は、より正確に自由界面を再構築できるシミュレーション装置、シミュレーション方法、および、シミュレーションプログラムを提供することである。

【課題を解決するための手段】

【0013】

本発明に従ったシミュレーション装置は、２種類以上の流体からなる混相流における自由界面をシミュレーションするためのシミュレーション装置であって、シミュレーション対象が分割された複数のセルのそれぞれにおける２種類以上の流体のうちの１の流体の占有率αを算出するための算出手段と、変数βおよびγを含む式（Ａ１）および式（Ａ２）を用いるＴＨＩＮＣ法に従って、各セルの第１の軸に沿った界面χ（ｘ）の再構築を行なうための再構築手段とを備え、

【0014】

【数2】

【0015】

再構築手段は、各セルの第１の軸に沿った方向についての変数βを、各セルの法線ベクトルの、当該第１の軸に対する勾配として算出する。

【0016】

本発明に従ったシミュレーション方法は、２種類以上の流体からなる混相流における自由界面をシミュレーションするためのコンピュータによって実行されるシミュレーション方法であって、シミュレーション対象が分割された複数のセルのそれぞれにおける２種類以上の流体のうちの１の流体の占有率αを算出するステップと、変数βおよびγを含む式（Ａ１）および式（Ａ２）を用いるＴＨＩＮＣ法に従って、各セルの第１の軸に沿った界面の再構築を行なうステップと、

【0017】

【数3】

【0018】

式（Ａ１）における変数βを求めるための定数であって、シミュレーションにおいて定義される第１の軸を含む２以上の軸ごとに設定された定数を記憶するステップとを備え、再構築を行なうステップは、各セルの第１の軸に沿った方向についての変数βを、各セルの法線ベクトルの、当該第１の軸に対する勾配として算出する。

【0019】

本発明に従ったシミュレーションプログラムは、２種類以上の流体からなる混相流における自由界面をシミュレーションするためのコンピュータによって実行されるシミュレーションプログラムであって、シミュレーションプログラムは、コンピュータに、シミュレーション対象が分割された複数のセルのそれぞれにおける２種類以上の流体のうちの１の流体の占有率αを算出するステップと、変数βおよびγを含む式（Ａ１）および式（Ａ２）を用いるＴＨＩＮＣ法に従って、各セルの第１の軸に沿った界面の再構築を行なうステップと、

【0020】

【数4】

【0021】

式（Ａ１）における変数βを求めるための定数であって、シミュレーションにおいて定義される第１の軸を含む２以上の軸ごとに設定された定数を記憶するステップとを実行させ、再構築を行なうステップは、各セルの第１の軸に沿った方向についての変数βを、各セルの法線ベクトルの、当該第１の軸に対する勾配として算出する。

【発明の効果】

【0022】

本発明によれば、シミュレーションにおける自由界面の再構築の精度を向上させることができる。

【図面の簡単な説明】

【0023】

【図1】この発明の実施の形態に従うシミュレーション装置を実現するための代表的なハードウェア構成であるコンピュータの概略構成図である。

【図2】ＴＨＩＮＣ法による、一次元の占有率の変化を示す図である。

【図3】実施の形態における変数の算出方法を説明するための図である。

【図4】ＴＨＩＮＣ／ＷＬＩＣ法の概略を示す図である。

【図5】ＴＨＡＩＮＣ法の概略を説明するための図である。

【図6】本実施の形態の計算方法の概略を説明するための図である。

【図7】本実施の形態のシミュレーションにおける処理のフローチャートである。

【図8】図７の処理内容を説明するための図である。

【図9】図７の処理内容を説明するための図である。

【図10】図７の処理内容を説明するための図である。

【図11】図７の処理内容を説明するための図である。

【図12】図７の処理内容を説明するための図である。

【図13】図１１に対応した各セルの密度を示す図である。

【図14】図１２に対応した各セルの密度を示す図である。

【図15】（x,y,z）座標系において中心座標を（i,j,k）とするセルＣi,j,kの、ｘ軸方向の法線ベクトルｎ_ｘの算出を説明するための図である。

【図16】（x,y,z）座標系において中心座標を（i,j,k）とするセルＣi,j,kの、ｘ軸方向の法線ベクトルｎ_ｘの算出を説明するための図である。

【図17】（x,y,z）座標系において中心座標を（i,j,k）とするセルＣi,j,kの、ｘ軸方向の法線ベクトルｎ_ｘの算出を説明するための図である。

【図18】図１７の線Ｌ１の断面を模式的に示す。

【図19】図１７の線Ｌ２の断面を模式的に示す。

【図20】３次元シミュレーションにおいて定義されるセルの格子の一例を模式的に示す図である。

【図21】３次元シミュレーションにおける各セルの占有率の計算結果の一例を模式的に示す図である。

【図22】３次元シミュレーションにおける各セルの密度の計算結果の一例を模式的に示す図である。

【図23】液相中の気泡を５回回転移動させるシミュレーションをした際の界面の再構築の結果を示す図である。

【図24】２次元シミュレーションにおいて、液相中の気泡を移動させた際の界面の再構築の結果を示す図である。

【図25】３次元シミュレーションにおいて、液相中の気泡を移動させた際の界面の再構築の結果を示す図である。

【図26】本実施の形態の計算方法におけるβの内容を説明するための図である。

【図27】本実施の形態の計算方法におけるβの内容を説明するための図である。

【図28】従来の技術を説明するための図である。

【図29】従来の技術を説明するための図である。

【図30】従来の技術を説明するための図である。

【発明を実施するための形態】

【0024】

本発明の実施の形態について、図面を参照しながら詳細に説明する。なお、図中の同一または相当部分については、同一符号を付してその説明は繰返さない。

【0025】

＜ハードウェア構成＞
図１は、この発明の実施の形態に従うシミュレーション装置を実現するための代表的なハードウェア構成であるコンピュータ１００の概略構成図である。

【0026】

図１を参照して、コンピュータ１００は、ＦＤ（Flexible Disk）駆動装置１１１およびＣＤ−ＲＯＭ（Compact Disk-Read Only Memory)駆動装置１１３を搭載したコンピュータ本体１０１と、モニタ１０２と、キーボード１０３と、マウス１０４とを含む。

【0027】

コンピュータ本体１０１は、ＦＤ駆動装置１１１およびＣＤ−ＲＯＭ駆動装置１１３に加えて、相互にバスで接続された、演算装置であるＣＰＵ（Central Processing Unit）１０５と、メモリ１０６と、記憶装置である固定ディスク１０７と、通信インターフェース１０９とを含む。

【0028】

本実施の形態に従うシミュレーションは、ＣＰＵ１０５がメモリ１０６などのコンピュータハードウェアを用いて、プログラムを実行することで実現される。一般的に、このようなプログラムは、ＦＤ１１２やＣＤ−ＲＯＭ１１４などの記録媒体に格納されて、またはネットワークなどを介して流通する。そして、このようなプログラムは、ＦＤ駆動装置１１１やＣＤ−ＲＯＭ駆動装置１１３などにより記録媒体から読取られて、または通信インターフェース１０９にて受信されて、固定ディスク１０７に格納される。さらに、このようなプログラムは、固定ディスク１０７からメモリ１０６に読出されて、ＣＰＵ１０５により実行される。

【0029】

ＣＰＵ１０５は、様々な数値論理演算を行なう演算処理部であり、プログラムされた命令を順次実行することで、本実施の形態に従う制御機能を提供する。メモリ１０６は、ＣＰＵ１０５のプログラム実行に応じて各種の情報を記憶する。

【0030】

モニタ１０２は、ＣＰＵ１０５が出力する情報を表示するための表示部であって、一例としてＬＣＤ（Liquid Crystal Display）やＣＲＴ（Cathode Ray Tube）などから構成される。すなわち、モニタ１０２には、シミュレーション結果などが表示される。

【0031】

マウス１０４は、クリックやスライドなどの動作に応じたユーザから指令を受付ける。キーボード１０３は、入力されるキーに応じたユーザから指令を受付ける。

【0032】

通信インターフェース１０９は、コンピュータ１００と他の装置との間の通信を確立するための装置であり、各種データを外部から受付可能である。

【0033】

なお、上述したようなコンピュータに代えて、その一部または全部を専用のハードウェアによって実現してもよい。

【0034】

また、コンピュータ１００を本実施の形態に従うように機能させるプログラムが記録される記録媒体は、上記したＦＤやＣＤ−ＲＯＭに限定されない。記録媒体の他の例としては、ＤＶＤ−ＲＯＭ（Digital Versatile Disk - Read Only Memory）、ＵＳＢ（Universal Serial Bus）メモリ、メモリカード、ハードディスク、磁気テープ、カセットテープ、ＭＯ（Magnetic Optical Disc）、ＭＤ（Mini Disc）、ＩＣ（Integrated Circuit）カード（メモリカードを除く）、光カード、マスクＲＯＭ、ＥＰＲＯＭ、ＥＥＰＲＯＭ（Electronically Erasable Programmable Read-Only Memory）などの、不揮発的にプログラムを格納する媒体が挙げられる。

【0035】

＜本実施の形態のシミュレーション方法の概要＞
上記したように、自由界面の再構築に関し、従来、ＴＨＩＮＣ法およびＴＨＩＮＣ／ＷＬＩＣ法が提案されている。一方、本実施の形態のシミュレーション方法は、ＴＨＩＮＣ法において値を固定されていた定数（β）を、法線ベクトルから各輸送方向に対する勾配として設定する。

【0036】

このことから、本明細書では、まず、本実施の形態のシミュレーション方法の前提となるＴＨＩＮＣ法およびＴＨＩＮＣ／ＷＬＩＣ法について説明した上で、本実施の形態のシミュレーション方法を説明する。また、本明細書では、参考として、上記定数（β）を本実施の形態とは異なる態様で設定するＴＨＡＩＮＣ（tangent of hyperbola with adaptive function for interface capturing）法についても説明を行なう。

【0037】

なお、本明細書では、本実施の形態のシミュレーション方法を、適宜「ＴＨＡＩＮＣ＿Ｍ（tangent of hyperbola with adaptive function for interface capturing_modified）法」と呼ぶ。

【0038】

また、本明細書では、シミュレーションにおいて定義される軸として、主にｘ軸とｙ軸の２軸を想定して説明がなされる。つまり、２次元でのシミュレーションが主に説明される。ただし、本実施の形態において開示されるシミュレーション方法は、公知の技術に従って３次元に拡張可能である。なお、３次元への拡張の方法については、公知の技術を採用可能であるため、ここでは説明は繰り返さない。

【0039】

＜ＴＨＩＮＣ法＞
図２は、ＴＨＩＮＣ法による、一次元の占有率の変化を示す図である。ＴＨＩＮＣ法は、図２に示されるように、平滑化された双曲線正接関数を利用している。なお、図２の破線は、体積追跡法による変化を示している。

【0040】

ＴＨＩＮＣ法では、特殊関数（characteristic function）として、次の式（２）が使用される。

【0041】

【数5】

【0042】

式（２）において、γは次の式（３）で表されるように、ｎ_ｘ，ｉの値に応じて「１」または「−１」の値をとる。

【0043】

【数6】

【0044】

式（３）のｎ_ｘ，ｉは、法線ベクトルのｘ成分であり、計算対象とするセルを３×３の格子状の中心に置いた場合の、残りの９個のセルの占有率αに基づいて算出される。具体的には、図３に示すように３×３の格子状の９個のセルを想定する。各セルの占有率は、上段では、左から、順に、α_{ｉ−１，ｊ＋１}、α_{ｉ，ｊ＋１}、α_{ｉ＋１，ｊ＋１}である。中段では、左から、順に、α_{ｉ−１，１}、α_ｉ，ｊ、α_{ｉ＋１，ｊ}である。下段では、左から、順に、α_{ｉ−１，ｊ−１}、α_{ｉ，ｊ−１}、α_{ｉ＋１，ｊ−１}である。そして、これらの占有率に基づき、次の式（４−１）〜式（４−４）および式（５）から、ｎ_ｘ，ｉが算出される。なお、これらの式における、「ＮＸＲＴ」「ＮＸＬＴ」「ＮＸＲＢ」「ＮＸＬＢ」は、それぞれ、図３に示される、中心のセルの４個の頂点における、ｘ軸方向に微小距離移動したときの占有率の変化量を示している。式（４−１）において、「ＮＸＲＴ」は、中心セルの右上部分の占有率の変化量を表している。式（４−２）において、「ＮＸＬＴ」は、中心セルの左上部分の占有率の変化量を表している。式（４−３）において、「ＮＸＲＢ」は、中心セルの右下部分の占有率の変化量を表している。式（４−４）において、「ＮＸＬＢ」は、中心セルの左下部分の占有率の変化量を表している。

【0045】

【数7】

【0046】

【数8】

【0047】

式（２）に戻って、βは、図２に示された双曲線正接関数の傾きを定めるパラメータである。一般的には、βとして、３．５という固定された値が採用される。

【0048】

式（２）中の「ｘ_ｉ〜」は、式（２）のχ（ｘ）の平均を用い、次の式（６）に従って算出される。なお、本明細書では、チルダ記号が「〜」で表現される。

【0049】

【数9】

【0050】

つまり、上記のように設定されたγ，βを用いて式（２）のχ（ｘ）が構築され、そして、式（６）に基づいて、「ｘ_ｉ〜」が算出される。

【0051】

＜ＴＨＩＮＣ／ＷＬＩＣ法＞
ＴＨＩＮＣ／ＷＬＩＣ法では、ｘ方向の流束は、式（７）に従って算出される。なお、式（７）中の「ＲＦ」はRight Faceを意味する。つまり、右端の面（後述する図５のｕで示された黒丸が位置する縦線が存在する面）を意味し、より具体的にはセル番号「ｉ＋１／２」を意味する。

【0052】

【数10】

【0053】

式（７）の右辺において、第１項はＴＨＩＮＣ法に従った流束であり、第２項はＵＰＷＩＮＤ法（一次風上法）に従った流束である。

【0054】

式（７）において、ω_ｘ，ｉは、次の式（８）に従って算出される。

【0055】

【数11】

【0056】

式（８）において、ｎ_ｘ，ｉは、上記のように、式（４−１）〜式（４−４）および式（５）から算出される。また、ｎ_ｙ，ｉは、図３のような９個のセルの占有率に基づき、次の式（９−１）〜式（９−４）および式（１０）から算出される。

【0057】

【数12】

【0058】

【数13】

【0059】

つまり、ＴＨＩＮＣ／ＷＬＩＣ法では、図４に示されるように、ｎ_ｙ，ｉが０の場合には、式（８）に従ってω_ｘ，ｉが１になり、これにより、ｘ方向の流束は式（７）の右辺の第１項のみに基づいて算出される。また、ｎ_ｘ，ｉが０の場合には、式（８）に従ってω_ｘ，ｉが０になり、これにより、ｘ方向の流束は式（７）の右辺の第２項のみに基づいて算出される。

【0060】

そして、これらの場合の間の状態では、式（８）において、速度ベクトルのｘ成分とｙ成分の割合に応じてω_ｘ，ｉが算出される。これにより、ｘ方向の流束は、算出されたω_ｘ，ｉの値に基づき、式（７）の右辺の第１項と第２項が式（７）に従って組み合わされて、算出される。なお、式（７）の右辺の第１項は、ＴＨＩＮＣ法による流束の算出に対応し、第２項は、ＵＰＷＩＮＤ法による流束の算出に対応する。

【0061】

＜ＴＨＡＩＮＣ法＞
本実施の形態の比較例のシミュレーション方法であるＴＨＡＩＮＣ法について説明する。ＴＨＡＩＮＣ法は、ＴＨＩＮＣ法において値を固定されていた定数（β）を、軸ごとに設定された定数（後述するβｕ，βｌ）を用い、さらに、占有率を考慮する流体の速度ベクトルにおける軸ごとの寄与度（後述するω）に基づいて設定する方法を挙げる。

【0062】

図５は、ＴＨＡＩＮＣ法の概略を説明するための図である。
図５を参照して、ＴＨＡＩＮＣ法は、式（２）で示されたＴＨＩＮＣ法の特殊関数において、変数βの設定方法に特徴を有する。

【0063】

具体的には、変数βを、図５中の式（１１）として示されるように、β_ｕとβ_ｌとω_ｘ，ｉとを用いて算出して、設定する。β_ｕとβ_ｌの値は、たとえば固定ディスク１０７に格納されている。

【0064】

β_ｕは、従来βとして利用されてきた値３．５よりも界面付近を急激に変化させる値であり、たとえば１０とされる。なお、図５のβ_ｕは、式（２）においてγ≧０の場合の結果が示されている。β_ｌは、従来βとして利用されてきた値３．５よりも界面付近の変化を緩くさせる値であり、たとえば０．００１とされる。図５には、β_ｕとβ_ｌのそれぞれが式（２）中のβとして利用された場合の、界面付近の占有率の変化が示されている。

【0065】

ＴＨＡＩＮＣ法によれば、βは、法線ベクトルの大小関係（２次元シミュレーションであれば、ｎ_ｘ，ｉとｎ_ｙ，ｉの大小関係）に応じて設定される。

【0066】

具体的には、ｎ_ｘ，ｉがｎ_ｙ，ｉに比べて大きいときには、つまり、セルの速度ベクトルのｘ成分がｙ成分よりも大きい場合には、当該セルの再構築において、従来よりも、ｘ方向の界面の流体の占有率の変化をより急激にするように、βが設定される。式（１１）に従ったβの設定において、β_ｕの方がβ_ｌよりも大きい割合で反映されるからである。

【0067】

一方、ｎ_ｘ，ｉがｎ_ｙ，ｉに比べて小さいときには、つまり、セルの速度ベクトルのｘ成分がｙ成分よりも小さい場合には、当該セルの再構築において、従来よりも、ｘ方向の界面の流体の占有率の変化がより緩やかになるように、βが設定される。式（１１）に従ったβの設定において、β_ｌの方がβ_ｕよりも大きい割合で反映されるからである。

【0068】

つまり、本実施の形態に従った界面の再構築では、複数の軸を対象とする場での速度ベクトルにおける、再構築の軸方向の成分の寄与度（本実施の形態では、ｎ_ｘ，ｉとｎ_ｙ，ｉの和に対するｎ_ｘ，ｉの寄与度）が、より大きく反映されることになる。

【0069】

ＴＨＡＩＮＣ法における計算フローは、図７を参照して説明された本実施の形態のシミュレーションにおける計算フローと同様とすることができる。

【0070】

つまり、まずステップＳ１０において、各セルの占有率αと速度ｖの初期値が読み込まれ、ステップＳ２０で、各セルの速度ｖならびに当該速度のｘ方向およびｙ方向の法線ベクトルが算出され、ステップＳ３０で、最新の占有率αについての「∇α」が算出され、ステップＳ４０で、各セルの占有率αと速度ｖの積についての「∇αｖ」が算出され、ステップＳ５０で、各セルの占有率αに対してΔｔ後の占有率αが算出され、そして、ステップＳ６０で、最新の占有率αに基づいて、式（３）に従ってγを算出し、式（１１）に従ってβが算出され、式（２）および式（６）に従って「ｘｉ〜」が算出されることにより、各セルの界面が再構築される。

【0071】

そして、ステップＳ７０で、シミュレーションの期間（Δｔの総和）Ｔ_sumが指定された期間Ｔを超えると、処理が終了する。

【0072】

ＴＨＡＩＮＣ法についても、３次元シミュレーションへの拡張は、周知の技術を採用することによって実現可能であるため、ここでは詳細な説明は繰り返さない。たとえば、ω_ｘ，ｉは、３次元シミュレーションでは、次の式（１２）に示されるように、ｘ軸方向の法線ベクトルｎ_ｘ、ｙ軸方向の法線ベクトルｎ_ｙ、および、ｚ軸方向の法線ベクトルｎ_ｚを用いて表される。

【0073】

【数14】

【0074】

＜本実施の形態の計算方法＞
本実施の形態の計算方法では、基本的には、ＴＨＡＩＮＣ法に従った計算方法で、各セルの再構築が行なわれる。

【0075】

図６は、本実施の形態の計算方法の概略を説明するための図である。
図６を参照して、本実施の形態の計算方法は、式（２）で示されたＴＨＩＮＣ法の特殊関数において、変数βの設定方法に特徴を有する。具体的には、変数βを、図６中に式（１３）として示されるように算出して、設定する。

【0076】

式（１３）には、ｋ軸（シミュレーションにおいて定義される軸の１つ）方向の再構築に利用されるβが示されている。ｊは、式（１３）でのみ、シミュレーションにおいて定義される直行座標軸のｋ軸以外の軸を意味する。「ｋ」は、式（１３）でのみ、軸を表す文字として使用される。

【0077】

次の式（１４）〜式（１６）に、シミュレーションモデルが、１次元モデルの場合（１Ｄ）、２次元モデルの場合（２Ｄ）、３次元モデルの場合（３Ｄ）のそれぞれについてのβの算出式を示す。

【0078】

【数15】

【0079】

まず、シミュレーションモデルが１次元モデルである場合には、ｊ軸が存在しないため、βは、式（１４）に示されるように定数（３．５）とされる。

【0080】

シミュレーションモデルが２次元モデルである場合には、式（１５）に示されるように、βは、ｘ軸方向の法線ベクトルｎ_ｘ、および、ｙ軸方向の法線ベクトルｎ_ｙを用いて表される。具体的には、２次元シミュレーションにおいてｘ軸とｙ軸が定義されている場合、ｋ軸がｘ軸であればｊ軸はｙ軸であり、ｋ軸がｙ軸であればｊ軸はｘ軸である。そして、ｘ軸方向の再構築に利用されるβ（つまり、β_ｘ）は、ｎ_ｘをｎ_ｙで除した値として算出される。また、ｙ軸方向の再構築に利用されるβ（つまり、β_ｙ）は、ｎ_ｙをｎ_ｘで除した値として算出される。

【0081】

シミュレーションモデルが３次元モデルである場合には、式（１６）に示されるように、βは、ｘ軸方向の法線ベクトルｎ_ｘ、ｙ軸方向の法線ベクトルｎ_ｙ、および、ｚ軸方向の法線ベクトルｎ_ｚを用いて表される。具体的には、３次元シミュレーションにおいてｘ軸とｙ軸とｚ軸が定義されている場合、ｋ軸がｘ軸であればｊ軸はｙ軸とｚ軸であり、ｋ軸がｙ軸であればｊ軸はｚ軸とｘ軸であり、また、ｋ軸がｚ軸であればｊ軸はｘ軸とｙ軸である。そして、ｘ軸方向の再構築に利用されるβ（つまり、β_ｘ）は、ｎ_ｘをｎ_ｙとｎ_ｚの二乗の和の平方根で除した値として算出される。ｙ軸方向の再構築に利用されるβ（つまり、β_ｙ）は、ｎ_ｙをｎ_ｚとｎ_ｘの二乗の和の平方根で除した値として算出される。そして、ｚ軸方向の再構築に利用されるβ（つまり、β_ｚ）は、ｎ_ｚをｎ_ｘとｎ_ｙの二乗の和の平方根で除した値として算出される。

【0082】

＜本実施の形態における計算フロー＞
図７は、本実施の形態のシミュレーションにおける処理のフローチャートである。

【0083】

以下の説明では、シミュレーションの対象として、液体中の気体の移動を例示する。これにより、占有率αとして、各セルにおける気体の占有率が例示される。なお、本実施の形態のシミュレーション対象となる混相流における相の組合せはこれに限定されるものではない。異なる２以上の流体が存在する系であれば、異なる２以上の液体が存在する系であっても良い。

【0084】

図７を参照して、コンピュータ１００に対して計算開始の指示が入力されると、ＣＰＵ１０５は、まずステップＳ１０において、各セルの占有率αと速度ｖの初期値を読み込んで、ステップＳ２０へ処理を進める。なお、これらの値は、たとえば固定ディスク１０７に格納されている。

【0085】

ステップＳ２０では、ＣＰＵ１０５は、各セルの速度ｖならびに当該速度のｘ方向およびｙ方向の法線ベクトルを算出して、ステップＳ３０へ処理を進める。

【0086】

ステップＳ３０では、ＣＰＵ１０５は、最新の占有率αについての「∇α」を算出して、ステップＳ４０へ処理を進める。

【0087】

ステップＳ４０では、各セルの占有率αと速度ｖの積についての「∇αｖ」を算出して、ステップＳ５０へ処理を進める。

【0088】

ステップＳ５０では、各セルの占有率αに対してΔｔ後の占有率αを算出して、ステップＳ６０へ処理を進める。

【0089】

ステップＳ６０では、最新の占有率αに基づいて、式（３）に従ってγを算出し、式（１３）に従ってβを算出し、式（２）および式（６）に従って「ｘ_ｉ〜」を算出することにより、各セルの界面を再構築して、ステップＳ７０へ処理を進める。

【0090】

ステップＳ７０では、ＣＰＵ１０５は、シミュレーションの期間（Δｔの総和）Ｔsumが指定された期間Ｔを超えたか否かを判断し、まだ超えていなければステップＳ２０へ処理が戻され、超えていれば処理を終了させる。

【0091】

以上、図７を参照して説明した処理により、微小時間Δｔごとに、各セルの占有率が更新され、そして、更新後の占有率に従って各セルの界面が、図６等を参照して説明された本実施の形態の計算方法に従って、再構築される。

【0092】

なお、本実施の形態のシミュレーションでは、シミュレーションの対象の領域として、図８に示されるようなセルの格子（図８では１００×１００のセルの格子）が定義される。そして、ステップＳ２０からステップＳ６０の処理が一回行なわれるごとに、各セルの占有率αの計算結果が更新され、そして、各セルにおける界面の再構築が行なわれる。

【0093】

図９および図１０は、図８に示された格子状のセルの一部（１０×１０のセルの格子）についての、各セルの占有率αの計算結果を示す。なお、図１０は、図９の計算結果に対する、時間ΔＴ後の計算結果を示すものとする。

【0094】

また、図１１，図１２は、それぞれ、図９，図１０の計算結果に対応した、界面の再構築の状態を示す図である。つまり、図１２は、図１１に示された結果に対する、時間ΔＴ後の結果を示す。

【0095】

図９では、３〜６列目および３〜６行目にある１６個のセル以外のセルでは、占有率αは「０．０」とされている。一方、これらの１６個のセルでは、占有率αとして、０．４，０．８，１．０のいずれかの値が算出されている。そして、各占有率αに対応して、図１１に示されるように界面が再構築されている。

【0096】

図１０では、図９に対して、占有率αが０．０以外の値を有する１６個のセルが、それぞれ１列ずつ右方向に移動している。図１２では、これに対応して、再構築される界面が１列ずつ、図１１に示されたものに対して右に移動している。

【0097】

また、図１１および図１２に示されたような各セルの占有率αに基づき、次の式（１４）に従って、各セルの密度が算出される。図１１に対応した密度の算出結果を図１３に、図１２に対応した密度の算出結果を図１４に、それぞれ示す。なお、式（１４）において、ρ_Gは気体の密度であり、ρ_ｌは液体の密度である。また、図１３および図１４では、ρ_Gの値の一例として１kg/m³が、ρ_ｌの値の一例として１０００kg/m³が、それぞれ採用されている。

【0098】

【数16】

【0099】

＜３次元への拡張＞
本実施の形態のシミュレーション方法（ＴＨＡＩＮＣ）では、図６の式（１３）を参照して説明されたように、２次元以上の次元でシミュレーションが行なわれる場合、再構築対象のセルの各軸についての法線ベクトルが利用された。３次元シミュレーションにおける、各セルの法線ベクトルの算出について、説明する。

【0100】

図１５は、（x,y,z）座標系において中心座標を（i,j,k）とするセルＣ_i,j,kの、ｘ軸方向の法線ベクトルｎ_ｘの算出を説明するための図である。

【0101】

図１５では、Ｃ_i,j,kが立方体で示され、そして、当該立方体の８個の頂点近傍における、ｘ軸方向に微小距離移動したときの占有率の変化量が、それぞれ「ＮＸＲＴＦ」「ＮＸＬＴＦ」「ＮＸＲＢＦ」「ＮＸＬＢＦ」「ＮＸＲＴＢ」「ＮＸＬＴＢ」「ＮＸＲＢＢ」「ＮＸＬＢＢ」で示されている。

【0102】

これらの変化量は、次の式（１６−１）〜式（１６−８）に従って算出される。

【0103】

【数17】

【0104】

式（１６−１）〜式（１６−８）から理解されるように、上記８個の変化量の算出には、計算対象であるセルと当該セルに隣接するセルの占有率αが利用される。

【0105】

図１６に、上記８個の変化量の算出に利用されるセルを示す。
図１６では、計算対象のセルＣ_i,j,kを中心として３軸方向に３個ずつ配列された２７個のセルが示されている。２７個のセルについて、図１６では、ｘ軸方向に並ぶ３個のセルのそれぞれのｘ座標が「ｉ−１」「ｉ」「ｉ＋１」で示されている。また、ｙ軸方向に並ぶ３個のセルのｙ座標が「ｊ−１」「ｊ」「ｊ＋１」で示されている。また、ｚ軸方向に並ぶ３個のセルのｚ座標が「ｋ−１」「ｋ」「ｋ＋１」で示されている。つまり、２７個のセルの３次元座標は、（i-1,j-1,k-1）（i-1,j,k-1）（i-1,j+1,k-1）（i-1,j-1,k）（i-1,j,k）（i-1,j+1,k）（i-1,j-1,k+1）（i-1,j,k+1）（i-1,j+1,k+1）（i,j-1,k-1）（i,j,k-1）（i,j+1,k-1）（i,j-1,k）（i,j,k）（i,j+1,k）（i,j-1,k+1）（i,j,k+1）（i,j+1,k+1）（i+1,j-1,k-1）（i+1,j,k-1）（i+1,j+1,k-1）（i+1,j-1,k）（i+1,j,k）（i+1,j+1,k）（i+1,j-1,k+1）（i+1,j,k+1）（i+1,j+1,k+1）となる。

【0106】

そして、法線ベクトルｎ_ｘは、これらの８個の変化量を利用して、次の式（１９）に従って算出される。

【0107】

【数18】

【0108】

つまり、法線ベクトルｎ_ｘは、上記８個の変化量の平均値として算出される。
上記８個の変化量の算出に話を戻して、各変化量は、当該頂点近傍に位置する８個のセルの占有率を利用して算出される。図１６には、セルＣ_i,j,kの８個の頂点の１つである「ＮＸＬＴＦ」が示されている。

【0109】

たとえば、図１７に示されるように、「ＮＸＬＴＦ」で示される頂点近傍の、ｘ軸方向の微小領域における占有率の変化量は、計算対象のセルＣ_i,j,kと、ＮＸＬＴＦ近傍に位置する７個のセル（Ｃ_i,j+1,k，Ｃ_i,j+1,k-1，Ｃ_i,j,k-1，Ｃ_i-1,j+1,k，Ｃ_i-1,j,k-1）の占有率が利用される。Ｃに対する添え字は、各セルの座標を表す。これらの各セルの占有率は、α_i,j,k，α_i,j+1,k，α_i,j+1,k-1，α_i,j,k-1，α_i-1,j+1,k，α_i-1,j+1,k-1，α_i-1,j,k，α_i-1,j,k-1で表される。

【0110】

図１８は、図１７の線Ｌ１の断面を模式的に示す。当該断面は、図１７の８個のセルのｚ＝ｋにおける断面に相当する。また、図１９は、図１７の線Ｌ２の断面を模式的に示す。図１８と図１９には、それぞれのセルの占有率が示されている。当該断面は、図１７の８個のセルのｚ＝ｋ−１における断面に相当する。

【0111】

そして、「ＮＸＬＴＦ」は、式（１６−２）に示されるように、図１７の８個のセルについて、法線ベクトルを求める軸（ｘ軸）の方向に隣接する２つのセルの占有率の差をｘ軸方向の微小距離（△ｘ）で除した値の平均として算出される。つまり、これらの８個のセルにおいて、ｘ軸方向に隣接するセルの組合せは４通り（Ｃ_i,j+1,k-1とＣ_i-1,j+1,k-1、Ｃ_i,j,k-1とＣ_i-1,j,k-1、Ｃ_i,j+1,kとＣ_i-1,j+1,k、Ｃ_i,j,kとＣ_i-1,j,k）ある。そこで、式（１６−２）の分子では、これらの４組の各組における占有率の差がたしあわされているので、式（１６−２）の分母では、これらの平均を取るために、微小距離（△ｘ）とともに、「４」が記載されている。

【0112】

式（１６−１）〜式（１６−８）の他の式でも、同様に、セルＣ_i,j,kの８個の頂点付近の、ｘ軸の方向に隣接する２つのセルの占有率の差をｘ軸方向の微小距離（△ｘ）で除した値の平均が、算出される。そして、これらの式によって算出された結果が利用されて、式（１９）に従って、ｘ方向の法線ベクトルｎ_ｘが算出される。

【0113】

また、ｙ方向の法線ベクトルｎ_yとｚ方向の法線ベクトルｎ_zは、それぞれ、次の式（２０）と式（２１）に従って算出される。

【0114】

【数19】

【0115】

式（２０）中の「ＮＹＲＴＦ」「ＮＹＬＴＦ」「ＮＹＲＢＦ」「ＮＹＬＢＦ」「ＮＹＲＴＢ」「ＮＹＬＴＢ」「ＮＹＲＢＢ」「ＮＹＬＢＢ」の８個の数値、および、式（２１）中の「ＮＺＲＴＦ」「ＮＺＬＴＦ」「ＮＺＲＢＦ」「ＮＺＬＢＦ」「ＮＺＲＴＢ」「ＮＺＬＴＢ」「ＮＺＲＢＢ」「ＮＺＬＢＢ」の８個の数値は、それぞれ、図１５に示したような再構築の対象のセルＣ_i,j,kの８個の頂点付近における微小領域での占有率の変化率である。なお、式（２０）中の８個の数値は、ｙ方向の法線ベクトルｎ_yを求めるための数値であるため、ｙ軸方向に隣接する４組の２セルの占有率の差をｙ軸方向の微小距離（△ｙ）で除した値の平均として、算出される。また、式（２１）中の８個の数値は、ｚ方向の法線ベクトルｎ_zを求めるための数値であるため、ｚ軸方向に隣接する４組の２セルの占有率の差をｚ軸方向の微小距離（△ｚ）で除した値の平均として、算出される。

【0116】

以上説明したように、３次元シミュレーションでは、ｘ方向の法線ベクトルｎ_xとｙ方向の法線ベクトルｎ_yとｚ方向の法線ベクトルｎ_zとが算出される。そして、これらが式（１６）に利用されて、ｘ方向、ｙ方向、ｚ方向の各方向についてのβが算出されることにより、３次元シミュレーションに基づく各セルの再構築が実現される。

【0117】

３次元シミュレーションでは、シミュレーションの対象の領域として、３次元方向に配列されたセルの格子（図２０では１００×１００×１００のセルの格子）が定義される。そして、図７のステップＳ２０からステップＳ６０の処理が一回行なわれるごとに、図２１に示されるように各セルの占有率αの計算結果が更新され、更新後の各セルの占有率αが利用されて、各セルにおける界面の再構築が行なわれる。

【0118】

また、図２１に示されたような各セルの占有率αに基づき、上記の式（１７）に従って各セルの密度が算出される。図２２は、図２１に対応した密度の算出結果を示す図である。図２２では、ρ_Gの値の一例として１kg/m³が、ρ_ｌの値の一例として１０００kg/m³が、それぞれ採用されている。

【0119】

＜計算結果の比較＞
次に、本実施の形態による界面の再構築の結果を、これまでの種々の方法による再構築の結果と比較する。

【0120】

図２３〜図２５は、異なったシミュレーション内容についての、ＴＨＩＮＣ／ＷＬＩＣ法、ＴＨＡＩＮＣ法、および、ＴＨＡＩＮＣ＿Ｍ法のそれぞれによるシミュレーションの結果を示す図である。各図では、シミュレーションの結果がプロットされたグラフ、当該結果の数値を示すテーブル、および、シミュレーションの内容の説明図が、上から順に示されている。

【0121】

グラフでは、各方法の、気泡についての、初期形状に対する誤差（ＥＲＲ）が、クーラン数（ＣＯ）に対してプロットされている。誤差は、百分率で示されている。ＴＨＩＮＣ／ＷＬＩＣ法の結果はひし形で、ＴＨＡＩＮＣ法の結果は四角で、そして、ＴＨＡＩＮＣ＿Ｍ法の結果は三角で、それぞれプロットされている。この誤差は、次の式（２２）に従って、誤差（Ｅ_ＳＰ）として算出される。

【0122】

【数20】

【0123】

テーブルでは、クーラン数（ＣＯ）と、ＣＯに対応する誤差の百分率の値とが、示されている。

【0124】

（５回回転）
図２３は、液相中の気泡を５回回転移動させるシミュレーションをした際の界面の再構築の結果を示す図である。このシミュレーションは、説明図に示されるように、気泡ＢＵで示される初期状態の気泡を液体ＬＱ内を矢印Ｒ１で示される回転方向に５回転させて元の位置に戻すというものである。

【0125】

図２３のグラフおよびテーブルから理解されるように、ＴＨＡＩＮＣ＿Ｍ法は、示されたクーラン数のほぼ全域に渡って、ＴＨＩＮＣ／ＷＬＩＣ法に対して誤差が１％程度低く、また、ＴＨＡＩＮＣ法に対しても誤差が低くなっている。

【0126】

（２次元シミュレーションでの移動）
図２４は、２次元シミュレーションにおいて、液相中の気泡を移動させた際の界面の再構築の結果を示す図である。このシミュレーションは、説明図に示されるように、気泡ＢＵで示される初期状態の気泡を液体ＬＱ内を矢印Ｒ２に従って移動させるものである。

【0127】

図２４のグラフおよびテーブルから理解されるように、ＴＨＡＩＮＣ＿Ｍ法は、示されたクーラン数のほぼ全域に渡って、ＴＨＩＮＣ／ＷＬＩＣ法およびＴＨＡＩＮＣ法に対して誤差が低くなっている。特に、ＴＨＡＩＮＣ法に対し、ＴＨＡＩＮＣ＿Ｍ法は、ＣＯが０．５の場合には１．５％程度、ＣＯが０．２５の場合には３％近く、そして、ＣＯが０．１の場合には３％以上、誤差が低くなっている。

【0128】

（３次元シミュレーションでの移動）
図２５は、３次元シミュレーションにおいて、液相中の気泡を移動させた際の界面の再構築の結果を示す図である。このシミュレーションは、説明図に示されるように、気泡ＢＵで示される初期状態の気泡を液体ＬＱ内を矢印Ｒ３に従って３次元空間内を移動させるものである。図２５では、移動後の気泡が、気泡ＢＵ０１で示されている。

【0129】

図２５のグラフおよびテーブルから理解されるように、ＴＨＩＮＣ／ＷＬＩＣ法とＴＨＡＩＮＣ法の間では、ＣＯの値によって誤差の大小関係が入れ替わったりしているが、ＴＨＡＩＮＣ＿Ｍ法は、示されたクーラン数のほぼ全域に渡って、３％程度も、ＴＨＩＮＣ／ＷＬＩＣ法およびＴＨＡＩＮＣ法に対して誤差が低くなっている。

【0130】

＜実施の形態のまとめ＞
以上説明した本実施の形態のＴＨＡＩＮＣ＿Ｍ法は、２以上の軸が定義された場合に、式（２）（図６参照）に示されたＴＨＩＮＣ法の特殊関数のβの算出方法に特徴がある。具体的には、各セルの、各軸についての法線ベクトルを用いて、βが算出される。より具体的には、式（１５）および式（１６）等に示されるように、第１の軸方向についての再構築に用いるβは、当該第１の軸の法線ベクトルの、それ以外の軸の法線ベクトルの合成されたベクトルに対する比として、算出される。

【0131】

図２６には、２次元シミュレーションにおける或るセルについての法線ベクトルのｘ成分（ベクトルｎ_ｘ）とｙ成分（ベクトルｎ_ｙ）が、界面ＳＢおよび法線ベクトル（ベクトルｎ）とともに、模式的に示されている。図２６に示された例では、界面ＳＢの左右で注目している流体の占有率αの値が変化している。界面ＳＢより左方ではα＝０であり、右方ではα＝１である。そして、この場合には、式（１５）を参照して説明したように、ｘ軸方向の再構築に利用されるβ（つまり、β_ｘ）は、法線ベクトルｎのｘ軸方向に対する勾配として算出される。また、ｙ軸方向の再構築に利用されるβ（つまり、β_ｙ）は、法線ベクトルｎのｙ軸方向に対する勾配として算出される。

【0132】

図２７には、３次元シミュレーションにおける或るセルについての法線ベクトルのｘ成分（ベクトルｎ_ｘ）と、法線ベクトル（ベクトルｎ）と、そのｙ成分とｚ成分の合成ベクトル（ベクトルｎ_ｙｚ）とが、模式的に示されている。なお、図２７では、法線ベクトルｎのｙ成分がｙ軸上に「ｎ_ｙ」で示され、また、法線ベクトルｎのｚ成分がｚ軸上に「ｎ_ｚ」で示されている。また、当該合成ベクトル（ベクトルｎ_ｙｚ）のスカラー量は、次の式（２３）で表される。

【0133】

【数21】