(57)【要約】【課題】より正確に自由界面を再構築できるシミュレーション装置、シミュレーション方法、および、シミュレーションプログラムを提供する。【解決手段】本発明に係るシミュレーション装置では、ＴＨＩＮＣ法に従って、各セルの界面の再構築が行なわれる。なお、当該装置では、第１の軸に沿った界面の再構築において、当該第１の軸（ｋ軸）方向の再構築に利用される変数βは、各セルの法線ベクトルの、当該第１の軸に対する勾配として算出される。【選択図】図５
2014026379.tif 000002
【特許請求の範囲】
【請求項１】
２種類以上の流体からなる混相流における自由界面をシミュレーションするためのシミュレーション装置であって、
シミュレーション対象が分割された複数のセルのそれぞれにおける前記２種類以上の流体のうちの１の流体の占有率αを算出するための算出手段と、
変数βおよびγを含む式（Ａ１）および式（Ａ２）を用いるＴＨＩＮＣ（tangent  of  hyperbola  for  interface  capturing）法に従って、各前記セルの第１の軸に沿った界面χ（ｘ）の再構築を行なうための再構築手段とを備え、
【数１】
2014026379.tif 000024
前記再構築手段は、各前記セルの前記第１の軸に沿った方向についての前記変数βを、各前記セルの法線ベクトルの、当該第１の軸に対する勾配として算出する、シミュレーション装置。
【請求項２】
２種類以上の流体からなる混相流における自由界面をシミュレーションするためのコンピュータによって実行されるシミュレーション方法であって、
シミュレーション対象が分割された複数のセルのそれぞれにおける前記２種類以上の流体のうちの１の流体の占有率αを算出するステップと、
変数βおよびγを含む式（Ａ１）および式（Ａ２）を用いるＴＨＩＮＣ（tangent  of  hyperbola  for  interface  capturing）法に従って、各前記セルの第１の軸に沿った界面の再構築を行なうステップと、
【数２】
2014026379.tif 000025
前記式（Ａ１）における前記変数βを求めるための定数であって、前記シミュレーションにおいて定義される前記第１の軸を含む２以上の軸ごとに設定された定数を記憶するステップとを備え、
前記再構築を行なうステップは、各前記セルの前記第１の軸に沿った方向についての前記変数βを、各前記セルの法線ベクトルの、当該第１の軸に対する勾配として算出する、シミュレーション方法。
【請求項３】
２種類以上の流体からなる混相流における自由界面をシミュレーションするためのコンピュータによって実行されるシミュレーションプログラムであって、
前記シミュレーションプログラムは、前記コンピュータに、
シミュレーション対象が分割された複数のセルのそれぞれにおける前記２種類以上の流体のうちの１の流体の占有率αを算出するステップと、
変数βおよびγを含む式（Ａ１）および式（Ａ２）を用いるＴＨＩＮＣ（tangent  of  hyperbola  for  interface  capturing）法に従って、各前記セルの第１の軸に沿った界面の再構築を行なうステップと、
【数３】
2014026379.tif 000026
前記式（Ａ１）における前記変数βを求めるための定数であって、前記シミュレーションにおいて定義される前記第１の軸を含む２以上の軸ごとに設定された定数を記憶するステップとを実行させ、
前記再構築を行なうステップは、各前記セルの前記第１の軸に沿った方向についての前記変数βを、各前記セルの法線ベクトルの、当該第１の軸に対する勾配として算出する、シミュレーションプログラム。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、シミュレーション装置に関し、特に、混相流における自由界面を体積追跡法に従ってシミュレーションするための、シミュレーション装置、シミュレーション方法、および、シミュレーションプログラムに関する。
【背景技術】
【０００２】
従来から、混相流における自由界面の解析法について、種々の提案がなされている。その中の一つである体積追跡法は、体積占有率を用いて格子上に界面を表現する方法である。
【０００３】
体積占有率は、次の式（１）で示される。二相流系において、ある格子が完全に一方の流体で占められている場合は当該関数の値は１であり、完全に他方の流体で占められている場合は当該関数の値は０である。そして、この関数の値は、一方の流体の占める割合に応じて、０から１の値をとる。
【０００４】
【数１】
2014026379.tif 000003
【０００５】
式（１）において、ｔは時間を示し、αは各セルの一方の流体の占有率を示す。また、ｖベクトル（式（１）中のｖの上に矢印が付されて表記されたもの）は、上記一方の流体の、各セルの図２８に示されるような流速ｕの、再構築の方向についてのベクトル量である。
【０００６】
各セルの占有率について、図２９および図３０を参照して説明する。
図２９は、２つの流体が存在する状態を模式的に示す図である。図２９では、ハッチングを付されていない部分が一方の流体の存在領域を示し、ハッチングを付された部分が他方の流体の存在領域を示す。なお、図２９では、６×６の３６個のセルが定義された状態が示されている。そして、図２９の３６個のセルのそれぞれにおける、上記一方の流体の占有率を、図３０に示す。
【０００７】
図３０に示されるような各セルの或る流体の占有率から、図２９に示されるような界面を再構築するための方法も、従来から種々提案されている。
【０００８】
たとえば、非特許文献１には、ＴＨＩＮＣ（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
【発明の概要】
【発明が解決しようとする課題】
【００１１】
しかしながら、複雑な自由界面の流れのシミュレーションは、まだ現実の流れを正確に再現できているとは言えず、検討の余地が残されている。
【００１２】
本発明は、かかる実情に鑑み考え出されたものであり、その目的は、より正確に自由界面を再構築できるシミュレーション装置、シミュレーション方法、および、シミュレーションプログラムを提供することである。
【課題を解決するための手段】
【００１３】
本発明に従ったシミュレーション装置は、２種類以上の流体からなる混相流における自由界面をシミュレーションするためのシミュレーション装置であって、シミュレーション対象が分割された複数のセルのそれぞれにおける２種類以上の流体のうちの１の流体の占有率αを算出するための算出手段と、変数βおよびγを含む式（Ａ１）および式（Ａ２）を用いるＴＨＩＮＣ法に従って、各セルの第１の軸に沿った界面χ（ｘ）の再構築を行なうための再構築手段とを備え、
【００１４】
【数２】
2014026379.tif 000004
【００１５】
再構築手段は、各セルの第１の軸に沿った方向についての変数βを、各セルの法線ベクトルの、当該第１の軸に対する勾配として算出する。
【００１６】
本発明に従ったシミュレーション方法は、２種類以上の流体からなる混相流における自由界面をシミュレーションするためのコンピュータによって実行されるシミュレーション方法であって、シミュレーション対象が分割された複数のセルのそれぞれにおける２種類以上の流体のうちの１の流体の占有率αを算出するステップと、変数βおよびγを含む式（Ａ１）および式（Ａ２）を用いるＴＨＩＮＣ法に従って、各セルの第１の軸に沿った界面の再構築を行なうステップと、
【００１７】
【数３】
2014026379.tif 000005
【００１８】
式（Ａ１）における変数βを求めるための定数であって、シミュレーションにおいて定義される第１の軸を含む２以上の軸ごとに設定された定数を記憶するステップとを備え、再構築を行なうステップは、各セルの第１の軸に沿った方向についての変数βを、各セルの法線ベクトルの、当該第１の軸に対する勾配として算出する。
【００１９】
本発明に従ったシミュレーションプログラムは、２種類以上の流体からなる混相流における自由界面をシミュレーションするためのコンピュータによって実行されるシミュレーションプログラムであって、シミュレーションプログラムは、コンピュータに、シミュレーション対象が分割された複数のセルのそれぞれにおける２種類以上の流体のうちの１の流体の占有率αを算出するステップと、変数βおよびγを含む式（Ａ１）および式（Ａ２）を用いるＴＨＩＮＣ法に従って、各セルの第１の軸に沿った界面の再構築を行なうステップと、
【００２０】
【数４】
2014026379.tif 000006
【００２１】
式（Ａ１）における変数βを求めるための定数であって、シミュレーションにおいて定義される第１の軸を含む２以上の軸ごとに設定された定数を記憶するステップとを実行させ、再構築を行なうステップは、各セルの第１の軸に沿った方向についての変数βを、各セルの法線ベクトルの、当該第１の軸に対する勾配として算出する。
【発明の効果】
【００２２】
本発明によれば、シミュレーションにおける自由界面の再構築の精度を向上させることができる。
【図面の簡単な説明】
【００２３】
【図１】この発明の実施の形態に従うシミュレーション装置を実現するための代表的なハードウェア構成であるコンピュータの概略構成図である。
【図２】ＴＨＩＮＣ法による、一次元の占有率の変化を示す図である。
【図３】実施の形態における変数の算出方法を説明するための図である。
【図４】ＴＨＩＮＣ／ＷＬＩＣ法の概略を示す図である。
【図５】ＴＨＡＩＮＣ法の概略を説明するための図である。
【図６】本実施の形態の計算方法の概略を説明するための図である。
【図７】本実施の形態のシミュレーションにおける処理のフローチャートである。
【図８】図７の処理内容を説明するための図である。
【図９】図７の処理内容を説明するための図である。
【図１０】図７の処理内容を説明するための図である。
【図１１】図７の処理内容を説明するための図である。
【図１２】図７の処理内容を説明するための図である。
【図１３】図１１に対応した各セルの密度を示す図である。
【図１４】図１２に対応した各セルの密度を示す図である。
【図１５】（x,y,z）座標系において中心座標を（i,j,k）とするセルＣi,j,kの、ｘ軸方向の法線ベクトルｎｘの算出を説明するための図である。
【図１６】（x,y,z）座標系において中心座標を（i,j,k）とするセルＣi,j,kの、ｘ軸方向の法線ベクトルｎｘの算出を説明するための図である。
【図１７】（x,y,z）座標系において中心座標を（i,j,k）とするセルＣi,j,kの、ｘ軸方向の法線ベクトルｎｘの算出を説明するための図である。
【図１８】図１７の線Ｌ１の断面を模式的に示す。
【図１９】図１７の線Ｌ２の断面を模式的に示す。
【図２０】３次元シミュレーションにおいて定義されるセルの格子の一例を模式的に示す図である。
【図２１】３次元シミュレーションにおける各セルの占有率の計算結果の一例を模式的に示す図である。
【図２２】３次元シミュレーションにおける各セルの密度の計算結果の一例を模式的に示す図である。
【図２３】液相中の気泡を５回回転移動させるシミュレーションをした際の界面の再構築の結果を示す図である。
【図２４】２次元シミュレーションにおいて、液相中の気泡を移動させた際の界面の再構築の結果を示す図である。
【図２５】３次元シミュレーションにおいて、液相中の気泡を移動させた際の界面の再構築の結果を示す図である。
【図２６】本実施の形態の計算方法におけるβの内容を説明するための図である。
【図２７】本実施の形態の計算方法におけるβの内容を説明するための図である。
【図２８】従来の技術を説明するための図である。
【図２９】従来の技術を説明するための図である。
【図３０】従来の技術を説明するための図である。
【発明を実施するための形態】
【００２４】
本発明の実施の形態について、図面を参照しながら詳細に説明する。なお、図中の同一または相当部分については、同一符号を付してその説明は繰返さない。
【００２５】
＜ハードウェア構成＞
図１は、この発明の実施の形態に従うシミュレーション装置を実現するための代表的なハードウェア構成であるコンピュータ１００の概略構成図である。
【００２６】
図１を参照して、コンピュータ１００は、ＦＤ（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）などの、不揮発的にプログラムを格納する媒体が挙げられる。
【００３５】
＜本実施の形態のシミュレーション方法の概要＞
上記したように、自由界面の再構築に関し、従来、ＴＨＩＮＣ法およびＴＨＩＮＣ／ＷＬＩＣ法が提案されている。一方、本実施の形態のシミュレーション方法は、ＴＨＩＮＣ法において値を固定されていた定数（β）を、法線ベクトルから各輸送方向に対する勾配として設定する。
【００３６】
このことから、本明細書では、まず、本実施の形態のシミュレーション方法の前提となるＴＨＩＮＣ法およびＴＨＩＮＣ／ＷＬＩＣ法について説明した上で、本実施の形態のシミュレーション方法を説明する。また、本明細書では、参考として、上記定数（β）を本実施の形態とは異なる態様で設定するＴＨＡＩＮＣ（tangent  of  hyperbola  with  adaptive  function  for  interface  capturing）法についても説明を行なう。
【００３７】
なお、本明細書では、本実施の形態のシミュレーション方法を、適宜「ＴＨＡＩＮＣ＿Ｍ（tangent  of  hyperbola  with  adaptive  function  for  interface  capturing_modified）法」と呼ぶ。
【００３８】
また、本明細書では、シミュレーションにおいて定義される軸として、主にｘ軸とｙ軸の２軸を想定して説明がなされる。つまり、２次元でのシミュレーションが主に説明される。ただし、本実施の形態において開示されるシミュレーション方法は、公知の技術に従って３次元に拡張可能である。なお、３次元への拡張の方法については、公知の技術を採用可能であるため、ここでは説明は繰り返さない。
【００３９】
＜ＴＨＩＮＣ法＞
図２は、ＴＨＩＮＣ法による、一次元の占有率の変化を示す図である。ＴＨＩＮＣ法は、図２に示されるように、平滑化された双曲線正接関数を利用している。なお、図２の破線は、体積追跡法による変化を示している。
【００４０】
ＴＨＩＮＣ法では、特殊関数（characteristic  function）として、次の式（２）が使用される。
【００４１】
【数５】
2014026379.tif 000007
【００４２】
式（２）において、γは次の式（３）で表されるように、ｎｘ，ｉの値に応じて「１」または「−１」の値をとる。
【００４３】
【数６】
2014026379.tif 000008
【００４４】
式（３）のｎｘ，ｉは、法線ベクトルのｘ成分であり、計算対象とするセルを３×３の格子状の中心に置いた場合の、残りの９個のセルの占有率αに基づいて算出される。具体的には、図３に示すように３×３の格子状の９個のセルを想定する。各セルの占有率は、上段では、左から、順に、αｉ−１，ｊ＋１、αｉ，ｊ＋１、αｉ＋１，ｊ＋１である。中段では、左から、順に、αｉ−１，１、αｉ，ｊ、αｉ＋１，ｊである。下段では、左から、順に、αｉ−１，ｊ−１、αｉ，ｊ−１、αｉ＋１，ｊ−１である。そして、これらの占有率に基づき、次の式（４−１）〜式（４−４）および式（５）から、ｎｘ，ｉが算出される。なお、これらの式における、「ＮＸＲＴ」「ＮＸＬＴ」「ＮＸＲＢ」「ＮＸＬＢ」は、それぞれ、図３に示される、中心のセルの４個の頂点における、ｘ軸方向に微小距離移動したときの占有率の変化量を示している。式（４−１）において、「ＮＸＲＴ」は、中心セルの右上部分の占有率の変化量を表している。式（４−２）において、「ＮＸＬＴ」は、中心セルの左上部分の占有率の変化量を表している。式（４−３）において、「ＮＸＲＢ」は、中心セルの右下部分の占有率の変化量を表している。式（４−４）において、「ＮＸＬＢ」は、中心セルの左下部分の占有率の変化量を表している。
【００４５】
【数７】
2014026379.tif 000009
【００４６】
【数８】
2014026379.tif 000010
【００４７】
式（２）に戻って、βは、図２に示された双曲線正接関数の傾きを定めるパラメータである。一般的には、βとして、３．５という固定された値が採用される。
【００４８】
式（２）中の「ｘｉ〜」は、式（２）のχ（ｘ）の平均を用い、次の式（６）に従って算出される。なお、本明細書では、チルダ記号が「〜」で表現される。
【００４９】
【数９】
2014026379.tif 000011
【００５０】
つまり、上記のように設定されたγ，βを用いて式（２）のχ（ｘ）が構築され、そして、式（６）に基づいて、「ｘｉ〜」が算出される。
【００５１】
＜ＴＨＩＮＣ／ＷＬＩＣ法＞
ＴＨＩＮＣ／ＷＬＩＣ法では、ｘ方向の流束は、式（７）に従って算出される。なお、式（７）中の「ＲＦ」はRight Faceを意味する。つまり、右端の面（後述する図５のｕで示された黒丸が位置する縦線が存在する面）を意味し、より具体的にはセル番号「ｉ＋１／２」を意味する。
【００５２】
【数１０】
2014026379.tif 000012
【００５３】
式（７）の右辺において、第１項はＴＨＩＮＣ法に従った流束であり、第２項はＵＰＷＩＮＤ法（一次風上法）に従った流束である。
【００５４】
式（７）において、ωｘ，ｉは、次の式（８）に従って算出される。
【００５５】
【数１１】
2014026379.tif 000013
【００５６】
式（８）において、ｎｘ，ｉは、上記のように、式（４−１）〜式（４−４）および式（５）から算出される。また、ｎｙ，ｉは、図３のような９個のセルの占有率に基づき、次の式（９−１）〜式（９−４）および式（１０）から算出される。
【００５７】
【数１２】
2014026379.tif 000014
【００５８】
【数１３】
2014026379.tif 000015
【００５９】
つまり、ＴＨＩＮＣ／ＷＬＩＣ法では、図４に示されるように、ｎｙ，ｉが０の場合には、式（８）に従ってωｘ，ｉが１になり、これにより、ｘ方向の流束は式（７）の右辺の第１項のみに基づいて算出される。また、ｎｘ，ｉが０の場合には、式（８）に従ってωｘ，ｉが０になり、これにより、ｘ方向の流束は式（７）の右辺の第２項のみに基づいて算出される。
【００６０】
そして、これらの場合の間の状態では、式（８）において、速度ベクトルのｘ成分とｙ成分の割合に応じてωｘ，ｉが算出される。これにより、ｘ方向の流束は、算出されたωｘ，ｉの値に基づき、式（７）の右辺の第１項と第２項が式（７）に従って組み合わされて、算出される。なお、式（７）の右辺の第１項は、ＴＨＩＮＣ法による流束の算出に対応し、第２項は、ＵＰＷＩＮＤ法による流束の算出に対応する。
【００６１】
＜ＴＨＡＩＮＣ法＞
本実施の形態の比較例のシミュレーション方法であるＴＨＡＩＮＣ法について説明する。ＴＨＡＩＮＣ法は、ＴＨＩＮＣ法において値を固定されていた定数（β）を、軸ごとに設定された定数（後述するβｕ，βｌ）を用い、さらに、占有率を考慮する流体の速度ベクトルにおける軸ごとの寄与度（後述するω）に基づいて設定する方法を挙げる。
【００６２】
図５は、ＴＨＡＩＮＣ法の概略を説明するための図である。
図５を参照して、ＴＨＡＩＮＣ法は、式（２）で示されたＴＨＩＮＣ法の特殊関数において、変数βの設定方法に特徴を有する。
【００６３】
具体的には、変数βを、図５中の式（１１）として示されるように、βｕとβｌとωｘ，ｉとを用いて算出して、設定する。βｕとβｌの値は、たとえば固定ディスク１０７に格納されている。
【００６４】
βｕは、従来βとして利用されてきた値３．５よりも界面付近を急激に変化させる値であり、たとえば１０とされる。なお、図５のβｕは、式（２）においてγ≧０の場合の結果が示されている。βｌは、従来βとして利用されてきた値３．５よりも界面付近の変化を緩くさせる値であり、たとえば０．００１とされる。図５には、βｕとβｌのそれぞれが式（２）中のβとして利用された場合の、界面付近の占有率の変化が示されている。
【００６５】
ＴＨＡＩＮＣ法によれば、βは、法線ベクトルの大小関係（２次元シミュレーションであれば、ｎｘ，ｉとｎｙ，ｉの大小関係）に応じて設定される。
【００６６】
具体的には、ｎｘ，ｉがｎｙ，ｉに比べて大きいときには、つまり、セルの速度ベクトルのｘ成分がｙ成分よりも大きい場合には、当該セルの再構築において、従来よりも、ｘ方向の界面の流体の占有率の変化をより急激にするように、βが設定される。式（１１）に従ったβの設定において、βｕの方がβｌよりも大きい割合で反映されるからである。
【００６７】
一方、ｎｘ，ｉがｎｙ，ｉに比べて小さいときには、つまり、セルの速度ベクトルのｘ成分がｙ成分よりも小さい場合には、当該セルの再構築において、従来よりも、ｘ方向の界面の流体の占有率の変化がより緩やかになるように、βが設定される。式（１１）に従ったβの設定において、βｌの方がβｕよりも大きい割合で反映されるからである。
【００６８】
つまり、本実施の形態に従った界面の再構築では、複数の軸を対象とする場での速度ベクトルにおける、再構築の軸方向の成分の寄与度（本実施の形態では、ｎｘ，ｉとｎｙ，ｉの和に対するｎｘ，ｉの寄与度）が、より大きく反映されることになる。
【００６９】
ＴＨＡＩＮＣ法における計算フローは、図７を参照して説明された本実施の形態のシミュレーションにおける計算フローと同様とすることができる。
【００７０】
つまり、まずステップＳ１０において、各セルの占有率αと速度ｖの初期値が読み込まれ、ステップＳ２０で、各セルの速度ｖならびに当該速度のｘ方向およびｙ方向の法線ベクトルが算出され、ステップＳ３０で、最新の占有率αについての「∇α」が算出され、ステップＳ４０で、各セルの占有率αと速度ｖの積についての「∇αｖ」が算出され、ステップＳ５０で、各セルの占有率αに対してΔｔ後の占有率αが算出され、そして、ステップＳ６０で、最新の占有率αに基づいて、式（３）に従ってγを算出し、式（１１）に従ってβが算出され、式（２）および式（６）に従って「ｘｉ〜」が算出されることにより、各セルの界面が再構築される。
【００７１】
そして、ステップＳ７０で、シミュレーションの期間（Δｔの総和）Ｔsumが指定された期間Ｔを超えると、処理が終了する。
【００７２】
ＴＨＡＩＮＣ法についても、３次元シミュレーションへの拡張は、周知の技術を採用することによって実現可能であるため、ここでは詳細な説明は繰り返さない。たとえば、ωｘ，ｉは、３次元シミュレーションでは、次の式（１２）に示されるように、ｘ軸方向の法線ベクトルｎｘ、ｙ軸方向の法線ベクトルｎｙ、および、ｚ軸方向の法線ベクトルｎｚを用いて表される。
【００７３】
【数１４】
2014026379.tif 000016
【００７４】
＜本実施の形態の計算方法＞
本実施の形態の計算方法では、基本的には、ＴＨＡＩＮＣ法に従った計算方法で、各セルの再構築が行なわれる。
【００７５】
図６は、本実施の形態の計算方法の概略を説明するための図である。
図６を参照して、本実施の形態の計算方法は、式（２）で示されたＴＨＩＮＣ法の特殊関数において、変数βの設定方法に特徴を有する。具体的には、変数βを、図６中に式（１３）として示されるように算出して、設定する。
【００７６】
式（１３）には、ｋ軸（シミュレーションにおいて定義される軸の１つ）方向の再構築に利用されるβが示されている。ｊは、式（１３）でのみ、シミュレーションにおいて定義される直行座標軸のｋ軸以外の軸を意味する。「ｋ」は、式（１３）でのみ、軸を表す文字として使用される。
【００７７】
次の式（１４）〜式（１６）に、シミュレーションモデルが、１次元モデルの場合（１Ｄ）、２次元モデルの場合（２Ｄ）、３次元モデルの場合（３Ｄ）のそれぞれについてのβの算出式を示す。
【００７８】
【数１５】
2014026379.tif 000017
【００７９】
まず、シミュレーションモデルが１次元モデルである場合には、ｊ軸が存在しないため、βは、式（１４）に示されるように定数（３．５）とされる。
【００８０】
シミュレーションモデルが２次元モデルである場合には、式（１５）に示されるように、βは、ｘ軸方向の法線ベクトルｎｘ、および、ｙ軸方向の法線ベクトルｎｙを用いて表される。具体的には、２次元シミュレーションにおいてｘ軸とｙ軸が定義されている場合、ｋ軸がｘ軸であればｊ軸はｙ軸であり、ｋ軸がｙ軸であればｊ軸はｘ軸である。そして、ｘ軸方向の再構築に利用されるβ（つまり、βｘ）は、ｎｘをｎｙで除した値として算出される。また、ｙ軸方向の再構築に利用されるβ（つまり、βｙ）は、ｎｙをｎｘで除した値として算出される。
【００８１】
シミュレーションモデルが３次元モデルである場合には、式（１６）に示されるように、βは、ｘ軸方向の法線ベクトルｎｘ、ｙ軸方向の法線ベクトルｎｙ、および、ｚ軸方向の法線ベクトルｎｚを用いて表される。具体的には、３次元シミュレーションにおいてｘ軸とｙ軸とｚ軸が定義されている場合、ｋ軸がｘ軸であればｊ軸はｙ軸とｚ軸であり、ｋ軸がｙ軸であればｊ軸はｚ軸とｘ軸であり、また、ｋ軸がｚ軸であればｊ軸はｘ軸とｙ軸である。そして、ｘ軸方向の再構築に利用されるβ（つまり、βｘ）は、ｎｘをｎｙとｎｚの二乗の和の平方根で除した値として算出される。ｙ軸方向の再構築に利用されるβ（つまり、βｙ）は、ｎｙをｎｚとｎｘの二乗の和の平方根で除した値として算出される。そして、ｚ軸方向の再構築に利用されるβ（つまり、βｚ）は、ｎｚをｎｘとｎｙの二乗の和の平方根で除した値として算出される。
【００８２】
＜本実施の形態における計算フロー＞
図７は、本実施の形態のシミュレーションにおける処理のフローチャートである。
【００８３】
以下の説明では、シミュレーションの対象として、液体中の気体の移動を例示する。これにより、占有率αとして、各セルにおける気体の占有率が例示される。なお、本実施の形態のシミュレーション対象となる混相流における相の組合せはこれに限定されるものではない。異なる２以上の流体が存在する系であれば、異なる２以上の液体が存在する系であっても良い。
【００８４】
図７を参照して、コンピュータ１００に対して計算開始の指示が入力されると、ＣＰＵ１０５は、まずステップＳ１０において、各セルの占有率αと速度ｖの初期値を読み込んで、ステップＳ２０へ処理を進める。なお、これらの値は、たとえば固定ディスク１０７に格納されている。
【００８５】
ステップＳ２０では、ＣＰＵ１０５は、各セルの速度ｖならびに当該速度のｘ方向およびｙ方向の法線ベクトルを算出して、ステップＳ３０へ処理を進める。
【００８６】
ステップＳ３０では、ＣＰＵ１０５は、最新の占有率αについての「∇α」を算出して、ステップＳ４０へ処理を進める。
【００８７】
ステップＳ４０では、各セルの占有率αと速度ｖの積についての「∇αｖ」を算出して、ステップＳ５０へ処理を進める。
【００８８】
ステップＳ５０では、各セルの占有率αに対してΔｔ後の占有率αを算出して、ステップＳ６０へ処理を進める。
【００８９】
ステップＳ６０では、最新の占有率αに基づいて、式（３）に従ってγを算出し、式（１３）に従ってβを算出し、式（２）および式（６）に従って「ｘｉ〜」を算出することにより、各セルの界面を再構築して、ステップＳ７０へ処理を進める。
【００９０】
ステップＳ７０では、ＣＰＵ１０５は、シミュレーションの期間（Δｔの総和）Ｔsumが指定された期間Ｔを超えたか否かを判断し、まだ超えていなければステップＳ２０へ処理が戻され、超えていれば処理を終了させる。
【００９１】
以上、図７を参照して説明した処理により、微小時間Δｔごとに、各セルの占有率が更新され、そして、更新後の占有率に従って各セルの界面が、図６等を参照して説明された本実施の形態の計算方法に従って、再構築される。
【００９２】
なお、本実施の形態のシミュレーションでは、シミュレーションの対象の領域として、図８に示されるようなセルの格子（図８では１００×１００のセルの格子）が定義される。そして、ステップＳ２０からステップＳ６０の処理が一回行なわれるごとに、各セルの占有率αの計算結果が更新され、そして、各セルにおける界面の再構築が行なわれる。
【００９３】
図９および図１０は、図８に示された格子状のセルの一部（１０×１０のセルの格子）についての、各セルの占有率αの計算結果を示す。なお、図１０は、図９の計算結果に対する、時間ΔＴ後の計算結果を示すものとする。
【００９４】
また、図１１，図１２は、それぞれ、図９，図１０の計算結果に対応した、界面の再構築の状態を示す図である。つまり、図１２は、図１１に示された結果に対する、時間ΔＴ後の結果を示す。
【００９５】
図９では、３〜６列目および３〜６行目にある１６個のセル以外のセルでは、占有率αは「０．０」とされている。一方、これらの１６個のセルでは、占有率αとして、０．４，０．８，１．０のいずれかの値が算出されている。そして、各占有率αに対応して、図１１に示されるように界面が再構築されている。
【００９６】
図１０では、図９に対して、占有率αが０．０以外の値を有する１６個のセルが、それぞれ１列ずつ右方向に移動している。図１２では、これに対応して、再構築される界面が１列ずつ、図１１に示されたものに対して右に移動している。
【００９７】
また、図１１および図１２に示されたような各セルの占有率αに基づき、次の式（１４）に従って、各セルの密度が算出される。図１１に対応した密度の算出結果を図１３に、図１２に対応した密度の算出結果を図１４に、それぞれ示す。なお、式（１４）において、ρGは気体の密度であり、ρｌは液体の密度である。また、図１３および図１４では、ρGの値の一例として１kg/m3が、ρｌの値の一例として１０００kg/m3が、それぞれ採用されている。
【００９８】
【数１６】
2014026379.tif 000018
【００９９】
＜３次元への拡張＞
本実施の形態のシミュレーション方法（ＴＨＡＩＮＣ）では、図６の式（１３）を参照して説明されたように、２次元以上の次元でシミュレーションが行なわれる場合、再構築対象のセルの各軸についての法線ベクトルが利用された。３次元シミュレーションにおける、各セルの法線ベクトルの算出について、説明する。
【０１００】
図１５は、（x,y,z）座標系において中心座標を（i,j,k）とするセルＣi,j,kの、ｘ軸方向の法線ベクトルｎｘの算出を説明するための図である。
【０１０１】
図１５では、Ｃi,j,kが立方体で示され、そして、当該立方体の８個の頂点近傍における、ｘ軸方向に微小距離移動したときの占有率の変化量が、それぞれ「ＮＸＲＴＦ」「ＮＸＬＴＦ」「ＮＸＲＢＦ」「ＮＸＬＢＦ」「ＮＸＲＴＢ」「ＮＸＬＴＢ」「ＮＸＲＢＢ」「ＮＸＬＢＢ」で示されている。
【０１０２】
これらの変化量は、次の式（１６−１）〜式（１６−８）に従って算出される。
【０１０３】
【数１７】
2014026379.tif 000019
【０１０４】
式（１６−１）〜式（１６−８）から理解されるように、上記８個の変化量の算出には、計算対象であるセルと当該セルに隣接するセルの占有率αが利用される。
【０１０５】
図１６に、上記８個の変化量の算出に利用されるセルを示す。
図１６では、計算対象のセルＣ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）となる。
【０１０６】
そして、法線ベクトルｎｘは、これらの８個の変化量を利用して、次の式（１９）に従って算出される。
【０１０７】
【数１８】
2014026379.tif 000020
【０１０８】
つまり、法線ベクトルｎｘは、上記８個の変化量の平均値として算出される。
上記８個の変化量の算出に話を戻して、各変化量は、当該頂点近傍に位置する８個のセルの占有率を利用して算出される。図１６には、セルＣi,j,kの８個の頂点の１つである「ＮＸＬＴＦ」が示されている。
【０１０９】
たとえば、図１７に示されるように、「ＮＸＬＴＦ」で示される頂点近傍の、ｘ軸方向の微小領域における占有率の変化量は、計算対象のセルＣ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で表される。
【０１１０】
図１８は、図１７の線Ｌ１の断面を模式的に示す。当該断面は、図１７の８個のセルのｚ＝ｋにおける断面に相当する。また、図１９は、図１７の線Ｌ２の断面を模式的に示す。図１８と図１９には、それぞれのセルの占有率が示されている。当該断面は、図１７の８個のセルのｚ＝ｋ−１における断面に相当する。
【０１１１】
そして、「ＮＸＬＴＦ」は、式（１６−２）に示されるように、図１７の８個のセルについて、法線ベクトルを求める軸（ｘ軸）の方向に隣接する２つのセルの占有率の差をｘ軸方向の微小距離（△ｘ）で除した値の平均として算出される。つまり、これらの８個のセルにおいて、ｘ軸方向に隣接するセルの組合せは４通り（Ｃ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）ある。そこで、式（１６−２）の分子では、これらの４組の各組における占有率の差がたしあわされているので、式（１６−２）の分母では、これらの平均を取るために、微小距離（△ｘ）とともに、「４」が記載されている。
【０１１２】
式（１６−１）〜式（１６−８）の他の式でも、同様に、セルＣi,j,kの８個の頂点付近の、ｘ軸の方向に隣接する２つのセルの占有率の差をｘ軸方向の微小距離（△ｘ）で除した値の平均が、算出される。そして、これらの式によって算出された結果が利用されて、式（１９）に従って、ｘ方向の法線ベクトルｎｘが算出される。
【０１１３】
また、ｙ方向の法線ベクトルｎyとｚ方向の法線ベクトルｎzは、それぞれ、次の式（２０）と式（２１）に従って算出される。
【０１１４】
【数１９】
2014026379.tif 000021
【０１１５】
式（２０）中の「ＮＹＲＴＦ」「ＮＹＬＴＦ」「ＮＹＲＢＦ」「ＮＹＬＢＦ」「ＮＹＲＴＢ」「ＮＹＬＴＢ」「ＮＹＲＢＢ」「ＮＹＬＢＢ」の８個の数値、および、式（２１）中の「ＮＺＲＴＦ」「ＮＺＬＴＦ」「ＮＺＲＢＦ」「ＮＺＬＢＦ」「ＮＺＲＴＢ」「ＮＺＬＴＢ」「ＮＺＲＢＢ」「ＮＺＬＢＢ」の８個の数値は、それぞれ、図１５に示したような再構築の対象のセルＣi,j,kの８個の頂点付近における微小領域での占有率の変化率である。なお、式（２０）中の８個の数値は、ｙ方向の法線ベクトルｎyを求めるための数値であるため、ｙ軸方向に隣接する４組の２セルの占有率の差をｙ軸方向の微小距離（△ｙ）で除した値の平均として、算出される。また、式（２１）中の８個の数値は、ｚ方向の法線ベクトルｎzを求めるための数値であるため、ｚ軸方向に隣接する４組の２セルの占有率の差をｚ軸方向の微小距離（△ｚ）で除した値の平均として、算出される。
【０１１６】
以上説明したように、３次元シミュレーションでは、ｘ方向の法線ベクトルｎxとｙ方向の法線ベクトルｎyとｚ方向の法線ベクトルｎzとが算出される。そして、これらが式（１６）に利用されて、ｘ方向、ｙ方向、ｚ方向の各方向についてのβが算出されることにより、３次元シミュレーションに基づく各セルの再構築が実現される。
【０１１７】
３次元シミュレーションでは、シミュレーションの対象の領域として、３次元方向に配列されたセルの格子（図２０では１００×１００×１００のセルの格子）が定義される。そして、図７のステップＳ２０からステップＳ６０の処理が一回行なわれるごとに、図２１に示されるように各セルの占有率αの計算結果が更新され、更新後の各セルの占有率αが利用されて、各セルにおける界面の再構築が行なわれる。
【０１１８】
また、図２１に示されたような各セルの占有率αに基づき、上記の式（１７）に従って各セルの密度が算出される。図２２は、図２１に対応した密度の算出結果を示す図である。図２２では、ρGの値の一例として１kg/m3が、ρｌの値の一例として１０００kg/m3が、それぞれ採用されている。
【０１１９】
＜計算結果の比較＞
次に、本実施の形態による界面の再構築の結果を、これまでの種々の方法による再構築の結果と比較する。
【０１２０】
図２３〜図２５は、異なったシミュレーション内容についての、ＴＨＩＮＣ／ＷＬＩＣ法、ＴＨＡＩＮＣ法、および、ＴＨＡＩＮＣ＿Ｍ法のそれぞれによるシミュレーションの結果を示す図である。各図では、シミュレーションの結果がプロットされたグラフ、当該結果の数値を示すテーブル、および、シミュレーションの内容の説明図が、上から順に示されている。
【０１２１】
グラフでは、各方法の、気泡についての、初期形状に対する誤差（ＥＲＲ）が、クーラン数（ＣＯ）に対してプロットされている。誤差は、百分率で示されている。ＴＨＩＮＣ／ＷＬＩＣ法の結果はひし形で、ＴＨＡＩＮＣ法の結果は四角で、そして、ＴＨＡＩＮＣ＿Ｍ法の結果は三角で、それぞれプロットされている。この誤差は、次の式（２２）に従って、誤差（ＥＳＰ）として算出される。
【０１２２】
【数２０】
2014026379.tif 000022
【０１２３】
テーブルでは、クーラン数（ＣＯ）と、ＣＯに対応する誤差の百分率の値とが、示されている。
【０１２４】
（５回回転）
図２３は、液相中の気泡を５回回転移動させるシミュレーションをした際の界面の再構築の結果を示す図である。このシミュレーションは、説明図に示されるように、気泡ＢＵで示される初期状態の気泡を液体ＬＱ内を矢印Ｒ１で示される回転方向に５回転させて元の位置に戻すというものである。
【０１２５】
図２３のグラフおよびテーブルから理解されるように、ＴＨＡＩＮＣ＿Ｍ法は、示されたクーラン数のほぼ全域に渡って、ＴＨＩＮＣ／ＷＬＩＣ法に対して誤差が１％程度低く、また、ＴＨＡＩＮＣ法に対しても誤差が低くなっている。
【０１２６】
（２次元シミュレーションでの移動）
図２４は、２次元シミュレーションにおいて、液相中の気泡を移動させた際の界面の再構築の結果を示す図である。このシミュレーションは、説明図に示されるように、気泡ＢＵで示される初期状態の気泡を液体ＬＱ内を矢印Ｒ２に従って移動させるものである。
【０１２７】
図２４のグラフおよびテーブルから理解されるように、ＴＨＡＩＮＣ＿Ｍ法は、示されたクーラン数のほぼ全域に渡って、ＴＨＩＮＣ／ＷＬＩＣ法およびＴＨＡＩＮＣ法に対して誤差が低くなっている。特に、ＴＨＡＩＮＣ法に対し、ＴＨＡＩＮＣ＿Ｍ法は、ＣＯが０．５の場合には１．５％程度、ＣＯが０．２５の場合には３％近く、そして、ＣＯが０．１の場合には３％以上、誤差が低くなっている。
【０１２８】
（３次元シミュレーションでの移動）
図２５は、３次元シミュレーションにおいて、液相中の気泡を移動させた際の界面の再構築の結果を示す図である。このシミュレーションは、説明図に示されるように、気泡ＢＵで示される初期状態の気泡を液体ＬＱ内を矢印Ｒ３に従って３次元空間内を移動させるものである。図２５では、移動後の気泡が、気泡ＢＵ０１で示されている。
【０１２９】
図２５のグラフおよびテーブルから理解されるように、ＴＨＩＮＣ／ＷＬＩＣ法とＴＨＡＩＮＣ法の間では、ＣＯの値によって誤差の大小関係が入れ替わったりしているが、ＴＨＡＩＮＣ＿Ｍ法は、示されたクーラン数のほぼ全域に渡って、３％程度も、ＴＨＩＮＣ／ＷＬＩＣ法およびＴＨＡＩＮＣ法に対して誤差が低くなっている。
【０１３０】
＜実施の形態のまとめ＞
以上説明した本実施の形態のＴＨＡＩＮＣ＿Ｍ法は、２以上の軸が定義された場合に、式（２）（図６参照）に示されたＴＨＩＮＣ法の特殊関数のβの算出方法に特徴がある。具体的には、各セルの、各軸についての法線ベクトルを用いて、βが算出される。より具体的には、式（１５）および式（１６）等に示されるように、第１の軸方向についての再構築に用いるβは、当該第１の軸の法線ベクトルの、それ以外の軸の法線ベクトルの合成されたベクトルに対する比として、算出される。
【０１３１】
図２６には、２次元シミュレーションにおける或るセルについての法線ベクトルのｘ成分（ベクトルｎｘ）とｙ成分（ベクトルｎｙ）が、界面ＳＢおよび法線ベクトル（ベクトルｎ）とともに、模式的に示されている。図２６に示された例では、界面ＳＢの左右で注目している流体の占有率αの値が変化している。界面ＳＢより左方ではα＝０であり、右方ではα＝１である。そして、この場合には、式（１５）を参照して説明したように、ｘ軸方向の再構築に利用されるβ（つまり、βｘ）は、法線ベクトルｎのｘ軸方向に対する勾配として算出される。また、ｙ軸方向の再構築に利用されるβ（つまり、βｙ）は、法線ベクトルｎのｙ軸方向に対する勾配として算出される。
【０１３２】
図２７には、３次元シミュレーションにおける或るセルについての法線ベクトルのｘ成分（ベクトルｎｘ）と、法線ベクトル（ベクトルｎ）と、そのｙ成分とｚ成分の合成ベクトル（ベクトルｎｙｚ）とが、模式的に示されている。なお、図２７では、法線ベクトルｎのｙ成分がｙ軸上に「ｎｙ」で示され、また、法線ベクトルｎのｚ成分がｚ軸上に「ｎｚ」で示されている。また、当該合成ベクトル（ベクトルｎｙｚ）のスカラー量は、次の式（２３）で表される。
【０１３３】
【数２１】
2014026379.tif 000023
【０１３４】
ｘ軸方向の再構築に利用されるβ（つまり、βｘ）は、法線ベクトルｎのｘ軸方向に対する勾配として、つまり、式（１６）に示されたように、算出される。
【０１３５】
以上説明した本実施の形態のＴＨＡＩＮＣ＿Ｍ法によれば、２種類以上の流体からなる混相流における自由界面のシミュレーションにおいて、シミュレーションの次元が変わっても、その誤差を、従来の計算方法よりも低く抑えることができる。
【０１３６】
＜その他の変形例等＞
以上説明した本実施の形態では、ステップＳ５０を実行するＣＰＵ１０５によって、各セルの占有率を算出するための算出手段が構成されている。そして、ステップＳ６０を実行するＣＰＵ１０５によって、各セルの界面の再構築を行なう再構築手段が構成されている。
【符号の説明】
【０１３７】
１００  コンピュータ、１０１  コンピュータ本体、１０２  モニタ、１０３  キーボード、１０４  マウス、１０６  メモリ、１０７  固定ディスク、１０９  通信インターフェース、１１１  ＦＤ駆動装置、１１３  ＣＤ−ＲＯＭ駆動装置。
【図１】
2014026379.tif 000027
【図２】
2014026379.tif 000028
【図３】
2014026379.tif 000029
【図４】
2014026379.tif 000030
【図５】
2014026379.tif 000031
【図６】
2014026379.tif 000032
【図７】
2014026379.tif 000033
【図８】
2014026379.tif 000034
【図９】
2014026379.tif 000035
【図１０】
2014026379.tif 000036
【図１１】
2014026379.tif 000037
【図１２】
2014026379.tif 000038
【図１３】
2014026379.tif 000039
【図１４】
2014026379.tif 000040
【図１５】
2014026379.tif 000041
【図１６】
2014026379.tif 000042
【図１７】
2014026379.tif 000043
【図１８】
2014026379.tif 000044
【図１９】
2014026379.tif 000045
【図２０】
2014026379.tif 000046
【図２１】
2014026379.tif 000047
【図２２】
2014026379.tif 000048
【図２３】
2014026379.tif 000049
【図２４】
2014026379.tif 000050
【図２５】
2014026379.tif 000051
【図２６】
2014026379.tif 000052
【図２７】
2014026379.tif 000053
【図２８】
2014026379.tif 000054
【図２９】
2014026379.tif 000055
【図３０】
2014026379.tif 000056