"流体の表面張力をシミュレーションするシミュレーションプログラムにおいて、コンピュータに、入力された境界条件及び初期条件に基づいて、前記流体を粒子の集まりで表現した流体モデルの界面を算出させ、前記算出された界面の表面エネルギーを算出させ、前記算出された表面エネルギーに基づいて、前記界面の表面張力を算出させ、前記算出された表面張力に基づいて、所定時間間隔毎の前記流体の状態を出力させることを特徴とするシミュレーションプログラム。前記界面の表面エネルギーの算出は、空気との界面の表面エネルギーを算出する第1項と前記空気以外との界面の表面エネルギーを算出する第2項との和で表現した計算式を用いることを特徴とする請求項1記載のシミュレーションプログラム。前記空気との界面の表面エネルギーの算出は、前記第1項の表面張力係数と前記第2項の表面張力係数とを変えることを特徴とする請求項2記載のシミュレーションプログラム。流体の表面張力をシミュレーションするシミュレーション方法において、コンピュータが、入力された境界条件及び初期条件に基づいて、前記流体を粒子の集まりで表現した流体モデルの界面を算出し、前記算出された界面の表面エネルギーを算出し、前記算出された表面エネルギーに基づいて、前記界面の表面張力を算出し、前記算出された表面張力に基づいて、所定時間間隔毎の前記流体の状態を出力することを特徴とするシミュレーション方法。流体の表面張力をシミュレーションするシミュレーション装置において、境界条件及び初期条件を含む諸条件を入力する入力部と、前記入力された境界条件及び初期条件に基づいて、前記流体を粒子の集まりで表現した流体モデルの界面を算出する界面算出部と、前記算出された界面の表面エネルギーを算出する表面エネルギー算出部と、前記算出された表面エネルギーに基づいて、前記界面の表面張力を算出する表面張力算出部と、前記算出された表面張力に基づいて、所定時間間隔毎の前記流体の状態を出力する出力部を備えることを特徴とするシミュレーション装置。しかしながら、上述のような従来の技術においては、下記のような問題点がある。すなわち、表面のみに影響する表面張力の効果を、粒子法で表すためのモデルが、表面以外にある粒子も影響を受けるような形で実装されており、非物理的な挙動が発生する。特に濡れ角の動的な運動を制御することが難しくなってしまう。例えば、特許文献1、非特許文献1、非特許文献2、非特許文献3の手法による表面張力の計算では、表面だけではなく、内部の粒子にも影響が及んでしまう。特に、特許文献1の手法では、全粒子に対して、引力を働かせているため、内部の粒子にも表面張力による力がかかることになる。また、非特許文献1、非特許文献2、非特許文献3の手法では、周りの粒子の情報から、表面にある粒子を予想し、表面張力を導入する方法であるが、計算の中で、流体内部の粒子も、表面の粒子であると判定されてしまい、表面張力が加えられてしまうことが起こる。そのため、接触面を捕らえることが出来なくなり、濡れ角を持ちながら運動する場合などの扱いが困難となる。例えば、粒子の分布をここで、iは粒子の番号を表すとすると、非特許文献1、非特許文献2では粒子iが表面の粒子であると、判定して表面張力の計算を、以下の手順1乃至3で実行する。手順1以下の量を計算する。ここで、右辺のはカーネル関数重み関数と呼ばれることもあるで、非特許文献1、非特許文献2では以下のスプライン関数を用いている。ここで、は粒子間の影響半径で、初期状態の平均粒子間隔の2倍から3倍程度が良く用いられる。はカーネル関数の全空間積分量が1になるように調整された値で、2次元の場合は、3次元の場合はと決められる。上記式1は、図3に示すように、カーネル関数による重みを付けた、粒子iを中心とした半径h内の粒子の重心である。手順2式1の重心と粒子iの距離を計算する。手順3の値が平均粒子間隔の0.15倍以上であれば粒子iを表面の粒子と見なして、表面張力を加える。これらの手順1乃至3では、図4に示すように、周りの粒子分布が不均一になった粒子を境界の粒子と見なして表面張力を加えているが、内部の粒子が計算中にたまたま閾値平均粒子間隔の0.15倍を超えるを持ってしまった場合にも表面張力を適用してしまうことになる。1つの側面では、本発明は、シミュレーション結果が非物理的な挙動を示すことを抑制することを目的とする。本発明は、上記課題を解決するため、下記のような構成を採用した。すなわち、本発明の一態様によれば、流体の表面張力をシミュレーションするシミュレーションプログラム、シミュレーション方法又はシミュレーション装置において、コンピュータに、入力された境界条件及び初期条件に基づいて、前記流体を粒子の集まりで表現した流体モデルの界面を算出させ、前記算出された界面の表面エネルギーを算出させ、前記算出された表面エネルギーに基づいて、前記界面の表面張力を算出させ、前記算出された表面張力に基づいて、所定時間間隔毎の前記流体の状態を出力させることを特徴とするシミュレーションプログラム、シミュレーション方法又はシミュレーション装置が提供される。また、本発明は、前記界面の表面エネルギーの算出が、空気との界面の表面エネルギーを算出する第1項と前記空気以外との界面の表面エネルギーを算出する第2項との和で表現した計算式を用いることが望ましい。また、本発明は、前記空気との界面の表面エネルギーの算出が、前記第1項の表面張力係数と前記第2項の表面張力係数とを変えることが望ましい。開示されたシミュレーションプログラム、シミュレーション方法及びシミュレーション装置によれば、非物理的な挙動を示すことなく、適切なシミュレーション結果を出力することができる。以下、本発明の実施の形態について、図面を参照しながら詳細に説明する。本発明は、コンピュータに実行させるシミュレーションプログラム、シミュレーション方法及びシミュレーション装置であって、シミュレーションの対象である流体を粒子の集まりとして捕らえ、その粒子の分布点の分布から、凸包構成法を基にした計算幾何手法により、当該シミュレーションの対象である流体と当該流体以外の気体又は固体との境界を構成する界面にある粒子から界面を求める。そして、界面を構成する粒子を用いて表面エネルギーを表現し、表面エネルギーの第1変分を計算することにより、表面張力項の計算を行う。ここで、粒子法シミュレーションの対象である流体と界面を介して接する、液体状態にある相である液相、気体状態にある相である気相、固体状態にある相である固相の各々を他相という。また、表面エネルギーを構成する際に、表面張力の大きさを、各面毎に変えることで、流体の濡れ性や、付着現象も表現できるモデルを採用し、差分化する。図5は、本発明を適用したシミュレーション装置の構成例を示す図である。図5において、シミュレーション装置500は、処理部501、格納部502及び出力部503を備え、初期条件511に基づいて粒子法による数値計算を行って、シミュレーション結果512を出力する。格納部502は、本発明を適用したシミュレーションプログラムを実行するための各計算式の情報を格納する。処理部501は、第1の実施の形態乃至第3の実施の形態として後述する、本発明を適用したシミュレーション処理を実行する。出力部503は、処理部501が実行したシミュレーション結果512を出力する。第1の実施の形態第1の実施の形態として、2次元の界面抽出法について説明する。本第1の実施の形態における2次元の界面抽出は、まず、2次元上の粒子群を包括する閉曲線を求め、次に、閉曲線上の凹部を求めることで実現する。まず、閉曲線の求め方を説明する。2次元上の粒子群点群を考える。ここで、iは粒子の番号を表す。このような点群を包括する閉曲線を求める方法を説明する。図6は、点群を包括する閉曲線を求める方法を説明するための図である。本第1の実施の形態では、凸包構成法の1つであるGiftWrapping法を用いて、第1乃至第4の手順で閉曲線を求める方法を説明する。第1の手順粒子群の中から最小のx座標の値をもつ粒子をとする図6A。第2の手順参考となる線分と境界粒子以外の全ての粒子に対して、線分との成す角を測り、最小となるインデックスj2を探索し、とする図6B。第3の手順k=2から参考となる線分をとして、との成す角を計算し、最小となるインデックスj2を探索し、とする。第4の手順kを1インクリメントし、上記第3の手順をとなるまで繰り返す。このような手順により得られた点列を繋ぐと、図7に示すように、凸包となる閉曲線を求めることができる。図8は、点群を包括する閉曲線を求める処理の流れを示すフローチャートである。ステップS81において、粒子の入力データを取得する。このステップS81が上記第1の手順に相当する。ステップS82において、上記第2の手順である初期粒子、参照線分を決定する。そして、ステップS83において、全ての粒子に対して、参照線分と粒子間との成す格を計算する。これらステップS82及びステップS83が上記第2の手順に相当する。そして、全ての境界粒子に対して、ステップS84において、成す格が最小となる粒子を境界粒子とし、ステップS85において、参照成分を更新する。これらステップS84及びステップS85が上記第3の手順及び第4の手順に相当する。次に、上述のようにして求めた閉曲線凸包上の点列から、第5乃至第9の手順で凹部を求める求め方を説明する。第5の手順粒子iと粒子i1間の距離を計算する。第6の手順上記粒子間の距離が所定の閾値影響半径hの数倍程度よりも大きい場合長すぎる場合、図9に示すように、粒子群の中から線分に射影を落とした距離射影が線分上に無い場合は除くが最小となる粒子を新たな界面粒子の候補とする。第7の手順及びが両方ともより小さければ、をiとi1番目の界面粒子の間に挿入する。上述の大小関係が成立しない場合、後述の第9の手順に進む。第8の手順線分及びに対して、上記第5乃至第7の手順を繰り返し、距離が閾値より小さくなるか、界面候補の粒子が見つからなくなった場合には、処理を終了する。第9の手順上記第5乃至第8の手順を全てのiに対して実行する。このような手順により、図10に示すような凹部を求めることができる。図11は、凹部を求める処理の流れを示すフローチャートである。ステップS111において、図6乃至図8を用いて説明した凸包構成法による表面粒子データを取得し、ステップS112において、表面粒子iと次の表面粒子i1のデータを取得する。ステップS113において、粒子iと粒子i1の距離を計算する。このステップS113が上記第5の手順に相当する。ステップS114において、表面粒子間の距離が所定の閾値より長いか否か判断する。そして、長い場合ステップS114:Yes、ステップS115において、他の粒子の中から候補kを選択し、kを上記i1として上記ステップS113以降の処理を繰り返す。これらステップS114Yes及びステップS115からステップS113へ戻る流れが上記第6の手順乃至第8の手順に相当する。そして、ステップS112乃至ステップS114を全ての境界粒子iに対して実行する。これが上記第9の手順に相当する。第2の実施の形態第2の実施の形態として、表面張力項の計算法について説明する。本第2の実施の形態における表面張力項の計算は、上述した第1の実施の形態により得られた閉曲面から表面張力を計算する。上述の第1の実施の形態により、表面の粒子の位置は点列が得られており、2次元の場合、図12に示すように、表面はi、i1の粒子間で構成されているとする。3次元の場合は3点間の三角形要素などで表面は構成されているものとする。流体の表面エネルギーは一般に以下のように定義できる。ここで、上記式4の右辺の積分領域は、流体の境界面全域であり、は流体が空気に触れている部分では1、それ以外の部分では0の値をとる。はそれぞれ、図13に示すように、空気に触れている部分表面張力係数、固体に触れている部分の表面張力係数である。ここで、上記式4をどのように粒子法で表現するかが本発明のポイントとなる。本発明では、界面抽出法により得られた点列を用いて、以下のように表面エネルギーを求める。ここで、は、図14に示すように、式4のを幅で滑らかに近似した関数であり、は、3次元上の点と固体との距離を表す関数である。例えば、y=0の平面に固体壁がある場合はである。はのときは0、のときは1の値を取り、の範囲は何らかの補間線形補間や3次関数による補間等を行った関数形を用いる。積分範囲は、図15に示すように、線分の上で取られる。積分計算は、ガウス積分法などを用いて数値積分することができる。式5の右辺の第1項の和jは、粒子間j,j1の粒子間で取られ、第2項の和kは、全ての粒子で取るものである。右辺第1項は、式4の表面エネルギーを離散化したものに対して、右辺第2項は、粒子が壁を通り抜けないという境界条件を表したポテンシャルエネルギーであり、式6のような形式を用いると、粒子は固体との距離が0のときポテンシャルは無限大に発散するため、固体の壁を乗り越えて行く非物理的な挙動を抑えることができる。式5のエネルギーを用いて、表面の粒子にかかる力は以下のように計算できる。なお、式7は差分近似して求めても良い。上記式5は接触している面が、固体か気体かによって、表面張力係数が変わっている。図16に示すように、この差によって、接触面が変わる際の接触角が表現できる形式となっている。の場合は固体と接触している場合の方がエネルギー的に低くなるため、固体に付着しようとする力が発生することになる。粒子法では、固体を通り抜けないという条件を壁近傍のポテンシャル式6として表すため、粒子が固体に完全に接触することはない。そのため、本第2の実施の形態では、図17に示すように、式5のようにエネルギーを設定して、固体からの距離が以下の場合は反発力を受け、固体からの距離が以上未満の場合は付着力が働くようにした。図18は、表面張力項の計算の流れを示すフローチャートである。まず、ステップS181において、シミュレーションの対象とする流体のデータを取得し、ステップS182において、現在の速度を用いて流体粒子の位置を時間dt/2更新する。そして、ステップS183において、物理モデルにより流体粒子にかかる力を算出する。次に、ステップS184において、上述の第1の実施の形態により流体の表面粒子のデータを取得する。そして、ステップS185において、全ての境界粒子に対して、上述のようにして表面エネルギーを求め、表面粒子にかかる力を計算する。全ての境界粒子に対する計算が終了したら、ステップS186において、流体粒子の速度を更新し、ステップS187において、更新した速度を用いて流体粒子の位置を時間dt/2更新する。そして、上記ステップS182乃至ステップS187を繰り返し実行する。第3の実施の形態第3の実施の形態として、SPH法による非圧縮性粘性流体の運動を計算する際の表面張力導入方法について説明する。図19は、2次元流体の例を示す図である。図19に示すような状況下y=0の平面の上におかれた流体、y方向は鉛直上向きとするにおける、非圧縮性粘性流体の運動方程式を考える。は、それぞれ、流体の密度場、速度場、圧力場、音速である。式8は質量保存則、式9は運動量保存則、式10は状態方程式である。は粘性応力テンソルで、流体の粘性係数を定数とすると、である。式8は流体が集まってくるような速度場がある場合には密度は上昇し、逆に流体が離れていくような速度場の場合には密度は低下する効果を表す。式9の右辺第1項は圧力勾配項で、流体が圧力の大きい部分から圧力の小さい部分へ向かって力が発生する効果を表す。右辺第2項は粘性応力項で、流れにブレーキがかかるような効果をあらわす。右辺第3項は重力項である。は重力加速度、はy方向の単位ベクトルである。本第3の実施の形態では、上記式10のような密度と圧力のみの関係を用いているが、一般の温度や内部エネルギー、エントロピー等を用いた状態方程式を用いても良い。また、空気と接触部にはの表面張力係数定数、固体と接触部にはの表面張力係数定数がかかっているものとする。式8乃至式10に示したSPH法を用いて離散化し、上述の第2の実施の形態による表面張力項を導入すると、以下のようになる。ここで、添え字は粒子の番号を表す。すなわちa番目の粒子の位置ベクトル、速度ベクトル、密度、圧力はそれぞれである。はb番目の質量である。は粘性項を計算するために導入したパラメータで定数である。式12は、上記式9の運動方程式を粒子法により離散化したもので、右辺第2項は圧力勾配項、第3項は粘性応力項を表したものである。ここで、はのy方向成分、は表面エネルギーであり、上記式5と同様の形式となっており、具体的に離散化した以下の形式を採用する。ここではy=0にある壁からの反発力を表すポテンシャルエネルギーであり、を用いる。は定数である。は上述の第1の実施の形態により抽出された表面粒子であり、粒子群の何れか1つと対応している。上記式15の微分は以下のように計算する。仮に、が表面粒子と対応していない場合は、となり、仮に、が表面粒子と対応している場合は、を計算し、が求められる。さらに、は2次元の再標準化行列であり、である。ここで、、、、である。は粒子a、b間での1次元リーマン問題を解いて求めた時空間の中間値である。具体的には以下のように決定する。粒子a、bに関して、以下の特性量を定める。さらに、各勾配を以下のよう計算する。ここで、上記式22の速度vの上付き添え字は、各成分を表す。これらを用いて、は以下のように決定する。図20及び図21は、上記の手法により表面張力を加えて固体面との濡れ性を表現したものである。図20は、で、のときを示し、で、のときを示す。これらの図20及び図21は、静止しときの流体の形状である。それぞれ、表面張力係数を変化させることで、接触角の違いを表現することができている。上述したように、本発明は粒子法の計算で、表面張力が効果的な流体運動の計算に対して適用可能である。特に溶けた金属の流し込みや、樹脂の流し込みのシミュレーション等には効果的である。図5のシミュレーション装置は、例えば、図22に示すような情報処理装置コンピュータを用いて実現することが可能である。図22の情報処理装置は、CPUCentralProcessingUnit2201、メモリ2202、入力装置2203、出力装置2204、外部記憶装置2205、媒体駆動装置2206及びネットワーク接続装置2207を備える。これらはバス2208により互いに接続されている。メモリ2202は、例えば、ROMReadOnlyMemory、RAMRandomAccessMemory、フラッシュメモリ等の半導体メモリであり、シミュレーション処理に用いられるプログラム及びデータを格納する。例えば、CPU2201は、メモリ2202を利用してプログラムを実行することにより、上述のシミュレーション処理を行う。メモリ2202は、図5の格納部502としても使用できる。入力装置2203は、例えば、キーボード、ポインティングデバイス等であり、オペレータからの指示や情報の入力に用いられる。出力装置2204は、例えば、表示装置、プリンタ、スピーカ等であり、オペレータへの問い合わせや処理結果の出力に用いられる。出力装置2204は、図5の出力部503としても使用できる。外部記憶装置2205は、例えば、磁気ディスク装置、光ディスク装置、光磁気ディスク装置、テープ装置等である。この外部記憶装置2205には、ハードディスクドライブも含まれる。情報処理装置は、この外部記憶装置2205にプログラム及びデータを格納しておき、それらをメモリ2202にロードして使用することができる。媒体駆動装置2206は、可搬型記録媒体2209を駆動し、その記録内容にアクセスする。可搬型記録媒体2209は、メモリデバイス、フレキシブルディスク、光ディスク、光磁気ディスク等である。この可搬型記録媒体2209には、CompactDiskReadOnlyMemoryCD-ROM、DigitalVersatileDiskDVD、UniversalSerialBusUSBメモリ等も含まれる。オペレータは、この可搬型記録媒体2209にプログラム及びデータを格納しておき、それらをメモリ2202にロードして使用することができる。このように、シミュレーション処理に用いられるプログラム及びデータを格納するコンピュータ読み取り可能な記録媒体には、メモリ2202、外部記憶装置2205、及び可搬型記録媒体2209のような、物理的な非一時的な記録媒体が含まれる。ネットワーク接続装置2207は、通信ネットワーク2210に接続され、通信に伴うデータ変換を行う通信インタフェースである。情報処理装置は、プログラム及びデータを外部の装置からネットワーク接続装置2207を介して受け取り、それらをメモリ2202にロードして使用することができる。ネットワーク接続装置2207は、図5の出力部503としても使用できる。開示の実施形態とその利点について詳しく説明したが、当業者は、特許請求の範囲に明確に記載した本発明の範囲から逸脱することなく、様々な変更、追加、省略をすることができるであろう。"