JP5892257B2

JP5892257B2 - シミュレーションプログラム、シミュレーション方法及びシミュレーション装置

Info

Publication number: JP5892257B2
Application number: JP2014536514A
Authority: JP
Inventors: 正喜風間
Original assignee: Fujitsu Ltd
Current assignee: Fujitsu Ltd
Priority date: 2012-09-21
Filing date: 2012-09-21
Publication date: 2016-03-23
Anticipated expiration: 2032-09-21
Also published as: US10031984B2; JPWO2014045416A1; WO2014045416A1; US20150161305A1

Description

本発明は、シミュレーションプログラム、シミュレーション方法及びシミュレーション装置に関する。

近年の計算機パワーの向上に伴い、シミュレーションの手法も徐々に発展してきた。その結果として、様々な応用分野にシミュレーションが用いられるようになってきた。

流体や弾性体等の連続体の問題を解く数値計算手法としては、格子をベースにして微分方程式の近似解を求解する有限差分法、有限要素法、又は有限体積法等が多く用いられてきた。近年では数値計算をＣＡＥ（Computer Aided Engineering）などの応用分野で活用するため、これらの数値計算手法も発展し、流体と構造物が相互作用する問題が解かれるようになってきた。

しかしながら、有限要素法又は有限体積法等の格子を用いる手法では、自由表面などの界面の存在する問題や、流体・構造連成問題などの移動境界が発生する場合には、取り扱いが複雑なため、プログラム作成が困難である場合が多い。

これに対して格子を用いないＭＰＳ（Moving Particle Semi-implicit）法、ＳＰＨ（Smoothed Particle Hydrodynamics）法等の粒子法では、移動境界の取り扱いに特別な処置を必要としない。それゆえ、近年粒子法は広く用いられるようになっている。

自由表面などの、変形し移動する境界を容易に扱いたいがために開発された粒子法であるが、離散化方法から、連続体が単なる粒子群となってしまう場合には、図１に示すように、連続体の境界面がどこかはっきりしなくなる。そのため、粒子法では表面張力などの境界を陽に扱う必要がある問題に対しては、統一的手法は出来ておらず、個別に独自開発された技術が、いつくかあるという状態である（例えば、特許文献１、非特許文献１、非特許文献２、非特許文献３参照。）。

例えば、特許文献１は、図２に示すように、粒子間距離からポテンシャルを構成し、粒子間に引力を働かせる手法である。表面張力は表面での分子間力の不均一性から生じるものである、という描像を元に作られた手法である。

また、非特許文献１、非特許文献２、非特許文献３の手法は、近傍粒子の数が減ったものを境界にある粒子と見なして、表面張力を境界の粒子に加える手法である。
また、非特許文献２、非特許文献３では、濡れ性を表現するモデルも導入されている。

特開２００８−１１１６７５号公報

T. Hongo, M. Shigeta, S. Izawa, and Y. Fukunishi, "３次元非圧縮ＳＰＨ法における気液界面に作用する表面張力モデル", 第２３回数値流体力学シンポジウム, A8-5 (2009) M. Agawa, M. Shigeta, S. Izawa, and Y. Fukunishi, "斜面上を流下する液体の非圧縮ＳＰＨシミュレーション", 第２３回数値流体力学シンポジウム, A9-4 (2009) K. Nomura, S. Koshizuka, Y. Oka and H. Obata, "Numerical Analysis of Droplet Breakup Behavior using Particle Method", Journal of Nuclear Science and Technology, Vol. 38, No. 12, pp. 1057-1064(2001)

しかしながら、上述のような従来の技術においては、下記のような問題点がある。

すなわち、表面のみに影響する表面張力の効果を、粒子法で表すためのモデルが、表面以外にある粒子も影響を受けるような形で実装されており、非物理的な挙動が発生する。特に濡れ角の動的な運動を制御することが難しくなってしまう。

例えば、特許文献１、非特許文献１、非特許文献２、非特許文献３の手法による表面張力の計算では、表面だけではなく、内部の粒子にも影響が及んでしまう。

特に、特許文献１の手法では、全粒子に対して、引力を働かせているため、内部の粒子にも表面張力による力がかかることになる。

また、非特許文献１、非特許文献２、非特許文献３の手法は、周りの粒子の情報から、表面にある粒子を予想し、表面張力を導入する方法であるが、計算の中で、流体内部の粒子も、表面の粒子であると判定されてしまい、表面張力が加えられてしまうことが起こる。そのため、接触面を捕らえることが出来なくなり、濡れ角を持ちながら運動する場合などの扱いが困難となる。

例えば、粒子の分布を

（ここで、ｉは粒子の番号を表す）とすると、非特許文献１、非特許文献２では粒子ｉが表面の粒子であると、判定して表面張力の計算を、以下の手順１乃至３で実行する。

（手順１）
以下の量を計算する。

ここで、右辺の

はカーネル関数（重み関数と呼ばれることもある）で、非特許文献１、非特許文献２では以下のスプライン関数を用いている。

ここで、

は粒子間の影響半径で、初期状態の平均粒子間隔の２倍から３倍程度が良く用いられる。

はカーネル関数の全空間積分量が１になるように調整された値で、２次元の場合は

、３次元の場合は

と決められる。

上記（式１）は、図３に示すように、カーネル関数による重みを付けた、粒子iを中心とした半径ｈ内の粒子の重心である。

（手順２）
（式１）の重心と粒子iの距離を計算する。

（手順３）

の値が平均粒子間隔の所定倍数、例えば０．１５倍以上であれば粒子ｉを表面の粒子と見なして、表面張力を加える。

これらの手順１乃至３では、図４に示すように、周りの粒子分布が不均一になった粒子を境界の粒子と見なして表面張力を加えているが、内部の粒子が計算中にたまたま閾値“平均粒子間隔の０．１５倍”を超える

を持ってしまった場合にも表面張力を適用してしまうことになる。

１つの側面では、本発明は、シミュレーション結果が非物理的な挙動を示すことを抑制することを目的とする。

本発明は、上記課題を解決するため、下記のような構成を採用した。
すなわち、本発明の一態様によれば、流体の表面張力をシミュレーションするシミュレーションプログラム、シミュレーション方法又はシミュレーション装置において、コンピュータに、入力された境界条件及び初期条件に基づいて、前記流体を粒子の集まりで表現した流体モデルの界面を算出させ、前記算出された界面の表面エネルギーを算出させ、前記算出された表面エネルギーに基づいて、前記界面の表面張力を算出させ、前記算出された表面張力に基づいて、所定時間間隔毎の前記流体の状態を出力させることを特徴とするシミュレーションプログラム、シミュレーション方法又はシミュレーション装置が提供される。

また、本発明は、前記界面の表面エネルギーの算出が、空気との界面の表面エネルギーを算出する第１項と前記空気以外との界面の表面エネルギーを算出する第２項との和で表現した計算式を用いることが望ましい。

また、本発明は、前記空気との界面の表面エネルギーの算出が、前記第１項の表面張力係数と前記第２項の表面張力係数とを変えることが望ましい。

また、本発明は、前記空気との界面の表面エネルギーの算出が、前記粒子の表面分布の不均一性を補正するためのエネルギーを加算することが望ましい。

また、本発明は、前記空気との界面の表面エネルギーの算出が、前記流体の体積を補正するためのエネルギーを加算することが望ましい。

開示されたシミュレーションプログラム、シミュレーション方法及びシミュレーション装置によれば、非物理的な挙動を示すことなく、適切なシミュレーション結果を出力することができる。

粒子法の問題点を説明するための図である。粒子間の引力を説明するための図である。重みを付けた粒子の重心を示す図である。境界粒子の重心を示す図である。本発明を適用したシミュレーション装置の構成例を示す図である。点群を包括する閉曲線を求める方法を説明するための図である。凸包となる閉曲線を示す図である。点群を包括する閉曲線を求める処理の流れを示すフローチャートである。凹部を求める方法を説明するための図である。凹部を求める方法により修正された閉曲線を示す図である。凹部を求める処理の流れを示すフローチャートである。流体の表面粒子の点列を示す図である。表面張力係数を説明するための図である。表面エネルギーのΧ（カイ）を幅εで滑らかに近似した関数を説明するための図である。線分上での積分を説明するための図である。接触核の表現を説明するための図である。固体からの距離とエネルギーの関係を示す図である。表面張力項の計算の流れを示すフローチャートである。体積を戻す効果が示された図である。表面粒子の負均一性を抑える効果が示された図である。２次元流体の例を示す図である。表面張力を加えて固体面との濡れ性を表現した図（その１）である。表面張力を加えて固体面との濡れ性を表現した図（その２）である。補正項を入れなかった場合の例を示す図である。補正項を入れた場合の例を示す図である。情報処理装置の構成図である。

以下、本発明の実施の形態について、図面を参照しながら詳細に説明する。
本発明は、コンピュータに実行させるシミュレーションプログラム、シミュレーション方法及びシミュレーション装置であって、シミュレーションの対象である流体を粒子の集まりとして捕らえ、その粒子の分布（点の分布）から、凸包構成法を基にした計算幾何手法により、当該シミュレーションの対象である流体と当該流体以外の気体又は固体との境界を構成する界面にある粒子から界面を求める。そして、界面を構成する粒子を用いて表面エネルギーを表現し、表面エネルギーの第１変分を計算することにより、表面張力項の計算を行う。

ここで、粒子法シミュレーションの対象である流体と界面を介して接する、液体状態にある相である液相、気体状態にある相である気相、固体状態にある相である固相の各々を「他相」という。

また、表面エネルギーを構成する際に、表面張力の大きさを、各面毎に変えることで、流体の濡れ性や、付着現象も表現できるモデルを採用し、差分化する。

図５は、本発明を適用したシミュレーション装置の構成例を示す図である。
図５において、シミュレーション装置５００は、処理部５０１、格納部５０２及び出力部５０３を備え、初期条件５１１に基づいて粒子法による数値計算を行って、シミュレーション結果５１２を出力する。

格納部５０２は、本発明を適用したシミュレーションプログラムを実行するための各計算式の情報を格納する。

処理部５０１は、第１の実施の形態乃至第３の実施の形態として後述する、本発明を適用したシミュレーション処理を実行する。
出力部５０３は、処理部５０１が実行したシミュレーション結果５１２を出力する。

（第１の実施の形態）
第１の実施の形態として、２次元の界面抽出法について説明する。
本第１の実施の形態における２次元の界面抽出は、まず、２次元上の粒子群を包括する閉曲線を求め、次に、閉曲線上の凹部を求めることで実現する。

まず、閉曲線の求め方を説明する。
２次元上の粒子群（点群）

を考える。ここで、ｉは粒子の番号を表す。

このような点群を包括する閉曲線を求める方法を説明する。
図６は、点群を包括する閉曲線を求める方法を説明するための図である。
本第１の実施の形態では、凸包構成法の１つであるGift Wrapping法を用いて、第１乃至第４の手順で閉曲線を求める方法を説明する。

（第１の手順）
粒子群

の中から最小のｘ座標の値をもつ粒子を

とする（図６（Ａ））。

（第２の手順）
参考となる線分

と境界粒子以外の全ての粒子

に対して、線分

と

の成す角を測り、最小となるインデックスj2を探索し、

とする（図６（Ｂ））。

（第３の手順）
ｋ＝２から参考となる線分を

として、

との成す角を計算し、最小となるインデックスj2を探索し、

とする。

（第４の手順）
ｋを１インクリメントし、上記第３の手順を

となるまで繰り返す。

このような手順により得られた点列

を繋ぐと、図７に示すように、凸包となる閉曲線を求めることができる。

図８は、点群を包括する閉曲線を求める処理の流れを示すフローチャートである。

ステップＳ８１において、粒子の入力データを取得する。このステップＳ８１が上記第１の手順に相当する。

ステップＳ８２において、上記第２の手順である初期粒子、参照線分を決定する。そして、ステップＳ８３において、全ての粒子に対して、参照線分と粒子間との成す角を計算する。これらステップＳ８２及びステップＳ８３が上記第２の手順に相当する。

そして、全ての境界粒子に対して、ステップＳ８４において、成す角が最小となる粒子を境界粒子とし、ステップＳ８５において、参照成分を更新する。これらステップＳ８４及びステップＳ８５が上記第３の手順及び第４の手順に相当する。

次に、上述のようにして求めた閉曲線（凸包）上の点列

から、第５乃至第９の手順で凹部を求める求め方を説明する。

（第５の手順）
粒子ｉと粒子ｉ＋１間の距離

を計算する。

（第６の手順）
上記粒子間の距離

が所定の閾値（影響半径ｈの数倍程度）よりも大きい場合（長すぎる場合）、図９に示すように、粒子群の中から線分

に射影を落とした距離（射影が線分上に無い場合は除く）が最小となる粒子

を新たな界面粒子の候補とする。

（第７の手順）

及び

が両方とも

より小さければ、

をｉとｉ＋１番目の界面粒子の間に挿入する。
上述の大小関係が成立しない場合、後述の第９の手順に進む。

（第８の手順）
線分

及び

に対して、上記第５乃至第７の手順を繰り返し、距離が閾値より小さくなるか、界面候補の粒子が見つからなくなった場合には、処理を終了する。

（第９の手順）
上記第５乃至第８の手順を全てのiに対して実行する。
このような手順により、図１０に示すような凹部を求めることができる。

図１１は、凹部を求める処理の流れを示すフローチャートである。
ステップＳ１１１において、図６乃至図８を用いて説明した凸包構成法による表面粒子データを取得し、ステップＳ１１２において、表面粒子ｉと次の表面粒子ｉ＋１のデータを取得する。

ステップＳ１１３において、粒子ｉと粒子ｉ＋１の距離を計算する。このステップＳ１１３が上記第５の手順に相当する。

ステップＳ１１４において、表面粒子間の距離が所定の閾値より長いか否か判断する。そして、長い場合（ステップＳ１１４：Ｙｅｓ）、ステップＳ１１５において、他の粒子の中から候補ｋを選択し、ｋを上記ｉ＋１として上記ステップＳ１１３以降の処理を繰り返す。これらステップＳ１１４（Ｙｅｓ）及びステップＳ１１５からステップＳ１１３へ戻る流れが上記第６の手順乃至第８の手順に相当する。

そして、ステップＳ１１２乃至ステップＳ１１４を全ての境界粒子（ｉ）に対して実行する。これが上記第９の手順に相当する。

（第２の実施の形態）
第２の実施の形態として、表面張力項の計算法について説明する。
本第２の実施の形態における表面張力項の計算では、上述した第１の実施の形態により得られた閉曲面から表面張力を計算する。

上述の第１の実施の形態により、表面の粒子の位置として点列

が得られており、２次元の場合、図１２に示すように、表面はｉ、ｉ＋１の粒子間で構成されているとする。３次元の場合は３点間の三角形要素などで表面は構成されているものとする。

流体の表面エネルギーは一般に以下のように定義できる。

ここで、上記（式４）の右辺の積分領域は、流体の境界面全域であり、

は流体が空気に触れている部分では１、それ以外の部分では０の値をとる。

はそれぞれ、図１３に示すように、空気に触れている部分表面張力係数、固体に触れている部分の表面張力係数である。ここで、上記（式４）をどのように粒子法で表現するかが本発明のポイントとなる。

本発明では、界面抽出法により得られた点列

を用いて、以下のように表面エネルギーを求める。

ここで、

は、図１４に示すように、（式４）の

を幅

の間で、０から１までの曲線で滑らかに近似した関数であり、

は、３次元上の点

と固体との距離を表す関数である。例えば、ｙ＝０の平面に固体壁がある場合は

である。

は

のときは０、

のときは１の値を取り、

の範囲は何らかの補間（線形補間や３次関数による補間等）を行った関数形を用いる。

積分範囲

は、図１５に示すように、線分

の上で取られる。積分計算は、ガウス積分法などを用いて数値積分することができる。
（式５）の右辺の第１項の和（ｊ）は、粒子間（ｊ，ｊ＋１の粒子間）で取られ、第２項の和（ｋ）は、全ての粒子で取るものである。右辺第１項は、（式４）の表面エネルギーを離散化したものに対して、右辺第２項は、粒子が壁を通り抜けないという境界条件を表したポテンシャルエネルギーであり、（式６）のような形式を用いると、粒子は固体との距離が０のときポテンシャルは無限大に発散するため、固体の壁を乗り越えて行く非物理的な挙動を抑えることができる。

（式５）のエネルギーを用いて、表面の粒子

にかかる力は以下のように計算できる。

なお、（式７）は差分近似して求めても良い。

上記（式５）は接触している面が、固体か気体かによって、表面張力係数が変わっている。図１６に示すように、この差によって、接触面が変わる際の接触角が表現できる形式となっている。

の場合は固体と接触している場合の方がエネルギー的に低くなるため、固体に付着しようとする力が発生することになる。粒子法では、固体を通り抜けないという条件を壁近傍のポテンシャル（式６)として表すため、粒子が固体に完全に接触することはない。そのため、本第２の実施の形態では、図１７に示すように、（式５）のようにエネルギーを設定して、固体からの距離が

以下の場合は反発力を受け、固体からの距離が

以上

未満の場合は付着力が働くようにした。

図１８は、表面張力項の計算の流れを示すフローチャートである。
まず、ステップＳ１８１において、シミュレーションの対象とする流体のデータを取得し、ステップＳ１８２において、現在の速度を用いて流体粒子の位置を時間ｄｔ／２分だけ更新する。そして、ステップＳ１８３において、物理モデルにより流体粒子にかかる力を算出する。

次に、ステップＳ１８４において、上述の第１の実施の形態により流体の表面粒子のデータを取得する。

そして、ステップＳ１８５において、全ての境界粒子に対して、上述のようにして表面エネルギーを求め、表面粒子にかかる力を計算する。

全ての境界粒子に対する計算が終了したら、ステップＳ１８６において、流体粒子の速度を更新し、ステップＳ１８７において、更新した速度を用いて流体粒子の位置を時間ｄｔ／２分だけ更新する。

そして、上記ステップＳ１８２乃至ステップＳ１８７を繰り返し実行する。
さらに、上記エネルギーに代え、表面粒子の不均一性と流体の体積変化を抑えるため以下の体積エネルギー及び接点エネルギーを導入することも可能である。

体積エネルギー：

接点エネルギー：

ここで、

は流体の体積であり、

と計算できる。ここで

は外向き法線単位ベクトルであり、

は定数であり、

は表面粒子ｓ，ｉの質量である。

（式８）、（式９）のエネルギーによる補正項は以下のように計算できる。

は初期体積で、定数である。
体積エネルギーによる補正項は、図１９のように体積を

に戻す効果があり、接点エネルギーによる補正項は、図２０のように表面粒子の不均一性を抑える効果がある。

（第３の実施の形態）
第３の実施の形態として、ＳＰＨ法による非圧縮性粘性流体の運動を計算する際の表面張力導入方法について説明する。

図２１は、２次元流体の例を示す図である。
図２１に示すような状況下（ｙ＝０の平面の上におかれた流体、ｙ方向は鉛直上向きとする）における、非圧縮性粘性流体の運動方程式を考える。

は、それぞれ、流体の密度場、速度場、圧力場、音速である。

（式１１）は質量保存則、（式１２）は運動量保存則、（式１３）は状態方程式である。

は粘性応力テンソルで、流体の粘性係数を

（定数）とすると、

である。

（式１１）は流体が集まってくるような速度場がある場合には密度は上昇し、逆に流体が離れていくような速度場の場合には密度は低下する効果を表す。（式１２）の右辺第１項は圧力勾配項で、流体が圧力の大きい部分から圧力の小さい部分へ向かって力が発生する効果を表す。右辺第２項は粘性応力項で、流れにブレーキがかかるような効果をあらわす。右辺第３項は重力項である。

は重力加速度、

はｙ方向の単位ベクトルである。

本第３の実施の形態では、上記（式１３）のような密度と圧力のみの関係を用いているが、一般の温度や内部エネルギー、エントロピー等を用いた状態方程式を用いても良い。

また、空気との接触部には

の表面張力係数（定数）、固体との接触部には

の表面張力係数（定数）がかかっているものとする。

（式１１）乃至（式１３）に示したＳＰＨ法を用いて離散化し、上述の第２の実施の形態による表面張力項を導入すると、以下のようになる。

ここで、添え字は粒子の番号を表す。すなわちa番目の粒子の位置ベクトル、速度ベクトル、密度、圧力はそれぞれ

である。

はｂ番目の粒子の質量である。

は粘性項を計算するために導入したパラメータで定数である。

（式１５）は、上記（式１２）の運動方程式を粒子法により離散化したもので、右辺第２項は圧力勾配項、第３項は粘性応力項を表したものである。ここで、

は

のy方向成分、

は表面エネルギーであり、上記（式５）と同様の形式となっており、具体的に離散化した以下の形式を採用する。

ここで

はｙ＝０にある壁からの反発力を表すポテンシャルエネルギーであり、

を用いる。

は定数である。

は上述の第１の実施の形態により抽出された表面粒子であり、粒子群

の何れか１つと対応している。

上記（式１８）の微分

は以下のように計算する。
仮に、

が表面粒子

と対応していない場合は、

となり、仮に、

が表面粒子

と対応している場合は、

を計算し、

が求められる。

接点エネルギーは、２次元の場合以下のようになる。

ここで

，

はそれぞれｘ方向単位ベクトル、ｙ方向単位ベクトルである。

また、流体の体積

は

と計算できる。ここで

はベクトルを

回転させる行列である。

また、体積エネルギー

は（式８）のように計算でき、

となる。

さらに、

は２次元の再標準化行列であり、

である。

ここで、

、

である。

は粒子ａ、ｂ間での１次元リーマン問題を解いて求めた時空間の中間値である。具体的には以下のように決定する。

粒子ａ、ｂに関して、以下の特性量を定める。

さらに、各勾配を以下のよう計算する。

ここで、上記（式２７）の速度（ｖ）の上付き添え字は、各成分を表す。
これらを用いて、

は以下のように決定する。

図２２及び図２３は、上記の手法により表面張力を加えて固体面との濡れ性を表現したものである。図２２は、

で、

のときを示し、

で、

のときを示す。

これらの図２２及び図２３は、静止しときの流体の形状である。それぞれ、表面張力係数を変化させることで、接触角の違いを表現することができている。

また、図２４は、（式１９）及び（式２０）により計算される補正項を入れなかった場合を示し、図２５は、補正項を入れた場合を示す。

図２４では不自然な穴ができ、計算が不安定になるのに対して、図２５では安定に計算できており、補正項の効果が現れている。

シミュレーションのパラメータは以下の通りである。
なお、鉛直下向き方向をｙ軸とし、ｙ＝０を壁面とした。重力加速度は９．８

、粒子数は９００個、初期形状は長さ０．００１６７［ｍ］の格子上に粒子を置き正方形にした。

ここで、粘性係数

は０．０１

、

、影響半径ｈは０．０００３３３［ｍ］、

、密度は１０００

、表面張力は

、音速ｃは２３．７８

、補正項のパラメータは

とした。

上述したように、本発明は粒子法の計算で、表面張力が効果的な流体運動の計算に対して適用可能である。特に溶けた金属の流し込みや、樹脂の流し込みのシミュレーション等には効果的である。

図５のシミュレーション装置は、例えば、図２６に示すような情報処理装置（コンピュータ）を用いて実現することが可能である。図２６の情報処理装置は、ＣＰＵ（Central Processing Unit）２６０１、メモリ２６０２、入力装置２６０３、出力装置２６０４、外部記憶装置２６０５、媒体駆動装置２６０６及びネットワーク接続装置２６０７を備える。これらはバス２６０８により互いに接続されている。

メモリ２６０２は、例えば、ＲＯＭ（Read Only Memory）、ＲＡＭ（Random Access Memory）、フラッシュメモリ等の半導体メモリであり、シミュレーション処理に用いられるプログラム及びデータを格納する。例えば、ＣＰＵ２６０１は、メモリ２６０２を利用してプログラムを実行することにより、上述のシミュレーション処理を行う。メモリ２６０２は、図５の格納部５０２としても使用できる。

入力装置２６０３は、例えば、キーボード、ポインティングデバイス等であり、オペレータからの指示や情報の入力に用いられる。出力装置２６０４は、例えば、表示装置、プリンタ、スピーカ等であり、オペレータへの問い合わせや処理結果の出力に用いられる。出力装置２６０４は、図５の出力部５０３としても使用できる。

外部記憶装置２６０５は、例えば、磁気ディスク装置、光ディスク装置、光磁気ディスク装置、テープ装置等である。この外部記憶装置２６０５には、ハードディスクドライブも含まれる。情報処理装置は、この外部記憶装置２６０５にプログラム及びデータを格納しておき、それらをメモリ２６０２にロードして使用することができる。

媒体駆動装置２６０６は、可搬型記録媒体２６０９を駆動し、その記録内容にアクセスする。可搬型記録媒体２６０９は、メモリデバイス、フレキシブルディスク、光ディスク、光磁気ディスク等である。この可搬型記録媒体２６０９には、Compact Disk Read Only
Memory （ＣＤ−ＲＯＭ）、Digital Versatile Disk（ＤＶＤ）、Universal Serial Bus（ＵＳＢ）メモリ等も含まれる。オペレータは、この可搬型記録媒体２６０９にプログラム及びデータを格納しておき、それらをメモリ２６０２にロードして使用することができる。

このように、シミュレーション処理に用いられるプログラム及びデータを格納するコンピュータ読み取り可能な記録媒体には、メモリ２６０２、外部記憶装置２６０５、及び可搬型記録媒体２６０９のような、物理的な（非一時的な）記録媒体が含まれる。

ネットワーク接続装置２６０７は、通信ネットワーク２６１０に接続され、通信に伴うデータ変換を行う通信インタフェースである。情報処理装置は、プログラム及びデータを外部の装置からネットワーク接続装置２６０７を介して受け取り、それらをメモリ２６０２にロードして使用することができる。ネットワーク接続装置２６０７は、図５の出力部５０３としても使用できる。

開示の実施形態とその利点について詳しく説明したが、当業者は、特許請求の範囲に明確に記載した本発明の範囲から逸脱することなく、様々な変更、追加、省略をすることができるであろう。

Claims

流体の表面張力をシミュレーションするシミュレーションプログラムにおいて、
コンピュータに、
入力された境界条件及び初期条件に基づいて、前記流体を粒子の集まりで表現した流体モデルの界面を、凸包構成アルゴリズムを基に算出させ、
前記算出された界面の表面エネルギーを算出させ、
前記算出された表面エネルギーに基づいて、前記界面の表面張力を算出させ、
前記算出された表面張力に基づいて、所定時間間隔毎の前記流体の状態を出力させることを特徴とするシミュレーションプログラム。
前記界面の表面エネルギーの算出は、
空気との界面の表面エネルギーを算出する第１項と前記空気以外との界面の表面エネルギーを算出する第２項との和で表現した計算式を用いることを特徴とする請求項１記載のシミュレーションプログラム。
前記空気との界面の表面エネルギーの算出は、
前記第１項の表面張力係数と前記第２項の表面張力係数とを変えることを特徴とする請求項２記載のシミュレーションプログラム。
前記空気との界面の表面エネルギーの算出は、前記粒子の表面分布の不均一性を補正するためのエネルギーを加算することを特徴とする請求項３記載のシミュレーションプログラム。
前記空気との界面の表面エネルギーの算出は、前記流体の体積を補正するためのエネルギーを加算することを特徴とする請求項３記載のシミュレーションプログラム。
流体の表面張力をシミュレーションするシミュレーション方法において、
コンピュータが、
入力された境界条件及び初期条件に基づいて、前記流体を粒子の集まりで表現した流体モデルの界面を、凸包構成アルゴリズムを基に算出し、
前記算出された界面の表面エネルギーを算出し、
前記算出された表面エネルギーに基づいて、前記界面の表面張力を算出し、
前記算出された表面張力に基づいて、所定時間間隔毎の前記流体の状態を出力することを特徴とするシミュレーション方法。
流体の表面張力をシミュレーションするシミュレーション装置において、
境界条件及び初期条件を含む諸条件を入力する入力部と、
前記入力された境界条件及び初期条件に基づいて、前記流体を粒子の集まりで表現した流体モデルの界面を算出する凸包構成アルゴリズムを基に界面算出部と、
前記算出された界面の表面エネルギーを算出する表面エネルギー算出部と、
前記算出された表面エネルギーに基づいて、前記界面の表面張力を算出する表面張力算出部と、
前記算出された表面張力に基づいて、所定時間間隔毎の前記流体の状態を出力する出力部を備えることを特徴とするシミュレーション装置。