(57)【要約】  流体の表面張力をシミュレーションするシミュレーションプログラム、シミュレーション方法又はシミュレーション装置において、コンピュータに、入力された境界条件及び初期条件に基づいて、前記流体を粒子の集まりで表現した流体モデルの界面を算出させ、前記算出された界面の表面エネルギーを算出させ、前記算出された表面エネルギーに基づいて、前記界面の表面張力を算出させ、前記算出された表面張力に基づいて、所定時間間隔毎の前記流体の状態を出力させることを特徴とするシミュレーションプログラム、シミュレーション方法又はシミュレーション装置が提供される。
2014045416.tif 000002
【特許請求の範囲】
【請求項１】
流体の表面張力をシミュレーションするシミュレーションプログラムにおいて、
コンピュータに、
入力された境界条件及び初期条件に基づいて、前記流体を粒子の集まりで表現した流体モデルの界面を算出させ、
前記算出された界面の表面エネルギーを算出させ、
前記算出された表面エネルギーに基づいて、前記界面の表面張力を算出させ、
前記算出された表面張力に基づいて、所定時間間隔毎の前記流体の状態を出力させることを特徴とするシミュレーションプログラム。
【請求項２】
前記界面の表面エネルギーの算出は、
空気との界面の表面エネルギーを算出する第１項と前記空気以外との界面の表面エネルギーを算出する第２項との和で表現した計算式を用いることを特徴とする請求項１記載のシミュレーションプログラム。
【請求項３】
前記空気との界面の表面エネルギーの算出は、
前記第１項の表面張力係数と前記第２項の表面張力係数とを変えることを特徴とする請求項２記載のシミュレーションプログラム。
【請求項４】
前記空気との界面の表面エネルギーの算出は、前記粒子の表面分布の不均一性を補正するためのエネルギーを加算することを特徴とする請求項３記載のシミュレーションプログラム。
【請求項５】
前記空気との界面の表面エネルギーの算出は、前記流体の体積を補正するためのエネルギーを加算することを特徴とする請求項３記載のシミュレーションプログラム。
【請求項６】
流体の表面張力をシミュレーションするシミュレーション方法において、
コンピュータが、
入力された境界条件及び初期条件に基づいて、前記流体を粒子の集まりで表現した流体モデルの界面を算出し、
前記算出された界面の表面エネルギーを算出し、
前記算出された表面エネルギーに基づいて、前記界面の表面張力を算出し、
前記算出された表面張力に基づいて、所定時間間隔毎の前記流体の状態を出力することを特徴とするシミュレーション方法。
【請求項７】
流体の表面張力をシミュレーションするシミュレーション装置において、
境界条件及び初期条件を含む諸条件を入力する入力部と、
前記入力された境界条件及び初期条件に基づいて、前記流体を粒子の集まりで表現した流体モデルの界面を算出する界面算出部と、
前記算出された界面の表面エネルギーを算出する表面エネルギー算出部と、
前記算出された表面エネルギーに基づいて、前記界面の表面張力を算出する表面張力算出部と、
前記算出された表面張力に基づいて、所定時間間隔毎の前記流体の状態を出力する出力部を備えることを特徴とするシミュレーション装置。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、シミュレーションプログラム、シミュレーション方法及びシミュレーション装置に関する。
【背景技術】
【０００２】
近年の計算機パワーの向上に伴い、シミュレーションの手法も徐々に発展してきた。その結果として、様々な応用分野にシミュレーションが用いられるようになってきた。【０００３】
流体や弾性体等の連続体の問題を解く数値計算手法としては、格子をベースにして微分方程式の近似解を求解する有限差分法、有限要素法、又は有限体積法等が多く用いられてきた。近年では数値計算をＣＡＥ（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)
【発明の開示】
【００１１】
しかしながら、上述のような従来の技術においては、下記のような問題点がある。
【００１２】
すなわち、表面のみに影響する表面張力の効果を、粒子法で表すためのモデルが、表面以外にある粒子も影響を受けるような形で実装されており、非物理的な挙動が発生する。特に濡れ角の動的な運動を制御することが難しくなってしまう。
【００１３】
例えば、特許文献１、非特許文献１、非特許文献２、非特許文献３の手法による表面張力の計算では、表面だけではなく、内部の粒子にも影響が及んでしまう。
【００１４】
特に、特許文献１の手法では、全粒子に対して、引力を働かせているため、内部の粒子にも表面張力による力がかかることになる。
【００１５】
また、非特許文献１、非特許文献２、非特許文献３の手法は、周りの粒子の情報から、表面にある粒子を予想し、表面張力を導入する方法であるが、計算の中で、流体内部の粒子も、表面の粒子であると判定されてしまい、表面張力が加えられてしまうことが起こる。そのため、接触面を捕らえることが出来なくなり、濡れ角を持ちながら運動する場合などの扱いが困難となる。
【００１６】
例えば、粒子の分布を
【数１】
2014045416.tif 000003
（ここで、ｉは粒子の番号を表す）とすると、非特許文献１、非特許文献２では粒子ｉが表面の粒子であると、判定して表面張力の計算を、以下の手順１乃至３で実行する。
【００１７】
（手順１）
以下の量を計算する。
【００１８】
【数２】
2014045416.tif 000004
【００１９】
ここで、右辺の
2014045416.tif 000005
はカーネル関数（重み関数と呼ばれることもある）で、非特許文献１、非特許文献２では以下のスプライン関数を用いている。
【００２０】
【数３】
2014045416.tif 000006
【００２１】
ここで、
2014045416.tif 000007
は粒子間の影響半径で、初期状態の平均粒子間隔の２倍から３倍程度が良く用いられる。
2014045416.tif 000008
はカーネル関数の全空間積分量が１になるように調整された値で、２次元の場合は
【数４】
2014045416.tif 000009
、３次元の場合は
【数５】
2014045416.tif 000010
と決められる。
【００２２】
上記（式１）は、図３に示すように、カーネル関数による重みを付けた、粒子iを中心とした半径ｈ内の粒子の重心である。
【００２３】
（手順２）
（式１）の重心と粒子iの距離を計算する。
【数６】
2014045416.tif 000011
【００２４】
（手順３）
2014045416.tif 000012
の値が平均粒子間隔の所定倍数、例えば０．１５倍以上であれば粒子ｉを表面の粒子と見なして、表面張力を加える。
【００２５】
これらの手順１乃至３では、図４に示すように、周りの粒子分布が不均一になった粒子を境界の粒子と見なして表面張力を加えているが、内部の粒子が計算中にたまたま閾値“平均粒子間隔の０．１５倍”を超える
2014045416.tif 000013
を持ってしまった場合にも表面張力を適用してしまうことになる。
【００２６】
１つの側面では、本発明は、シミュレーション結果が非物理的な挙動を示すことを抑制することを目的とする。
【００２７】
本発明は、上記課題を解決するため、下記のような構成を採用した。
すなわち、本発明の一態様によれば、流体の表面張力をシミュレーションするシミュレーションプログラム、シミュレーション方法又はシミュレーション装置において、コンピュータに、入力された境界条件及び初期条件に基づいて、前記流体を粒子の集まりで表現した流体モデルの界面を算出させ、前記算出された界面の表面エネルギーを算出させ、前記算出された表面エネルギーに基づいて、前記界面の表面張力を算出させ、前記算出された表面張力に基づいて、所定時間間隔毎の前記流体の状態を出力させることを特徴とするシミュレーションプログラム、シミュレーション方法又はシミュレーション装置が提供される。
【００２８】
また、本発明は、前記界面の表面エネルギーの算出が、空気との界面の表面エネルギーを算出する第１項と前記空気以外との界面の表面エネルギーを算出する第２項との和で表現した計算式を用いることが望ましい。
【００２９】
また、本発明は、前記空気との界面の表面エネルギーの算出が、前記第１項の表面張力係数と前記第２項の表面張力係数とを変えることが望ましい。
【００３０】
また、本発明は、前記空気との界面の表面エネルギーの算出が、前記粒子の表面分布の不均一性を補正するためのエネルギーを加算することが望ましい。
【００３１】
また、本発明は、前記空気との界面の表面エネルギーの算出が、前記流体の体積を補正するためのエネルギーを加算することが望ましい。
【００３２】
開示されたシミュレーションプログラム、シミュレーション方法及びシミュレーション装置によれば、非物理的な挙動を示すことなく、適切なシミュレーション結果を出力することができる。
【図面の簡単な説明】
【００３３】
【図１】粒子法の問題点を説明するための図である。
【図２】粒子間の引力を説明するための図である。
【図３】重みを付けた粒子の重心を示す図である。
【図４】境界粒子の重心を示す図である。
【図５】本発明を適用したシミュレーション装置の構成例を示す図である。
【図６】点群を包括する閉曲線を求める方法を説明するための図である。
【図７】凸包となる閉曲線を示す図である。
【図８】点群を包括する閉曲線を求める処理の流れを示すフローチャートである。
【図９】凹部を求める方法を説明するための図である。
【図１０】凹部を求める方法により修正された閉曲線を示す図である。
【図１１】凹部を求める処理の流れを示すフローチャートである。
【図１２】流体の表面粒子の点列を示す図である。
【図１３】表面張力係数を説明するための図である。
【図１４】表面エネルギーのΧ（カイ）を幅εで滑らかに近似した関数を説明するための図である。
【図１５】線分上での積分を説明するための図である。
【図１６】接触核の表現を説明するための図である。
【図１７】固体からの距離とエネルギーの関係を示す図である。
【図１８】表面張力項の計算の流れを示すフローチャートである。
【図１９】体積を戻す効果が示された図である。
【図２０】表面粒子の負均一性を抑える効果が示された図である。
【図２１】２次元流体の例を示す図である。
【図２２】表面張力を加えて固体面との濡れ性を表現した図（その１）である。
【図２３】表面張力を加えて固体面との濡れ性を表現した図（その２）である。
【図２４】補正項を入れなかった場合の例を示す図である。
【図２５】補正項を入れた場合の例を示す図である。
【図２６】情報処理装置の構成図である。
【発明を実施するための形態】
【００３４】
以下、本発明の実施の形態について、図面を参照しながら詳細に説明する。
本発明は、コンピュータに実行させるシミュレーションプログラム、シミュレーション方法及びシミュレーション装置であって、シミュレーションの対象である流体を粒子の集まりとして捕らえ、その粒子の分布（点の分布）から、凸包構成法を基にした計算幾何手法により、当該シミュレーションの対象である流体と当該流体以外の気体又は固体との境界を構成する界面にある粒子から界面を求める。そして、界面を構成する粒子を用いて表面エネルギーを表現し、表面エネルギーの第１変分を計算することにより、表面張力項の計算を行う。
【００３５】
ここで、粒子法シミュレーションの対象である流体と界面を介して接する、液体状態にある相である液相、気体状態にある相である気相、固体状態にある相である固相の各々を「他相」という。
【００３６】
また、表面エネルギーを構成する際に、表面張力の大きさを、各面毎に変えることで、流体の濡れ性や、付着現象も表現できるモデルを採用し、差分化する。
【００３７】
図５は、本発明を適用したシミュレーション装置の構成例を示す図である。
図５において、シミュレーション装置５００は、処理部５０１、格納部５０２及び出力部５０３を備え、初期条件５１１に基づいて粒子法による数値計算を行って、シミュレーション結果５１２を出力する。
【００３８】
格納部５０２は、本発明を適用したシミュレーションプログラムを実行するための各計算式の情報を格納する。
【００３９】
処理部５０１は、第１の実施の形態乃至第３の実施の形態として後述する、本発明を適用したシミュレーション処理を実行する。
出力部５０３は、処理部５０１が実行したシミュレーション結果５１２を出力する。
【００４０】
（第１の実施の形態）
第１の実施の形態として、２次元の界面抽出法について説明する。
本第１の実施の形態における２次元の界面抽出は、まず、２次元上の粒子群を包括する閉曲線を求め、次に、閉曲線上の凹部を求めることで実現する。
【００４１】
まず、閉曲線の求め方を説明する。
２次元上の粒子群（点群）
2014045416.tif 000014
を考える。ここで、ｉは粒子の番号を表す。
【００４２】
このような点群を包括する閉曲線を求める方法を説明する。
図６は、点群を包括する閉曲線を求める方法を説明するための図である。
本第１の実施の形態では、凸包構成法の１つであるGift Wrapping法を用いて、第１乃至第４の手順で閉曲線を求める方法を説明する。
【００４３】
（第１の手順）
粒子群
2014045416.tif 000015
の中から最小のｘ座標の値をもつ粒子を
2014045416.tif 000016
とする（図６（Ａ））。
【００４４】
（第２の手順）
参考となる線分
【００４５】
【数７】
2014045416.tif 000017
と境界粒子以外の全ての粒子
2014045416.tif 000018
に対して、線分
2014045416.tif 000019
と
【数８】
2014045416.tif 000020
の成す角を測り、最小となるインデックスj2を探索し、
【数９】
2014045416.tif 000021
とする（図６（Ｂ））。
【００４６】
（第３の手順）
ｋ＝２から参考となる線分を
【数１０】
2014045416.tif 000022
として、
【数１１】
2014045416.tif 000023
との成す角を計算し、最小となるインデックスj2を探索し、
【数１２】
2014045416.tif 000024
とする。
【００４７】
（第４の手順）
ｋを１インクリメントし、上記第３の手順を
【数１３】
2014045416.tif 000025
となるまで繰り返す。
【００４８】
このような手順により得られた点列
2014045416.tif 000026
を繋ぐと、図７に示すように、凸包となる閉曲線を求めることができる。
【００４９】
図８は、点群を包括する閉曲線を求める処理の流れを示すフローチャートである。
【００５０】
ステップＳ８１において、粒子の入力データを取得する。このステップＳ８１が上記第１の手順に相当する。
【００５１】
ステップＳ８２において、上記第２の手順である初期粒子、参照線分を決定する。そして、ステップＳ８３において、全ての粒子に対して、参照線分と粒子間との成す角を計算する。これらステップＳ８２及びステップＳ８３が上記第２の手順に相当する。
【００５２】
そして、全ての境界粒子に対して、ステップＳ８４において、成す角が最小となる粒子を境界粒子とし、ステップＳ８５において、参照成分を更新する。これらステップＳ８４及びステップＳ８５が上記第３の手順及び第４の手順に相当する。
【００５３】
次に、上述のようにして求めた閉曲線（凸包）上の点列
2014045416.tif 000027
から、第５乃至第９の手順で凹部を求める求め方を説明する。
【００５４】
（第５の手順）
粒子ｉと粒子ｉ＋１間の距離
【数１４】
2014045416.tif 000028
を計算する。
【００５５】
（第６の手順）
上記粒子間の距離
2014045416.tif 000029
が所定の閾値（影響半径ｈの数倍程度）よりも大きい場合（長すぎる場合）、図９に示すように、粒子群の中から線分
【数１５】
2014045416.tif 000030
に射影を落とした距離（射影が線分上に無い場合は除く）が最小となる粒子
2014045416.tif 000031
を新たな界面粒子の候補とする。
【００５６】
（第７の手順）
【数１６】
2014045416.tif 000032
及び
【数１７】
2014045416.tif 000033
が両方とも
【数１８】
2014045416.tif 000034
より小さければ、
2014045416.tif 000035
をｉとｉ＋１番目の界面粒子の間に挿入する。
上述の大小関係が成立しない場合、後述の第９の手順に進む。
【００５７】
（第８の手順）
線分
【数１９】
2014045416.tif 000036
及び
【数２０】
2014045416.tif 000037
に対して、上記第５乃至第７の手順を繰り返し、距離が閾値より小さくなるか、界面候補の粒子が見つからなくなった場合には、処理を終了する。
【００５８】
（第９の手順）
上記第５乃至第８の手順を全てのiに対して実行する。
このような手順により、図１０に示すような凹部を求めることができる。
【００５９】
図１１は、凹部を求める処理の流れを示すフローチャートである。
ステップＳ１１１において、図６乃至図８を用いて説明した凸包構成法による表面粒子データを取得し、ステップＳ１１２において、表面粒子ｉと次の表面粒子ｉ＋１のデータを取得する。
【００６０】
ステップＳ１１３において、粒子ｉと粒子ｉ＋１の距離を計算する。このステップＳ１１３が上記第５の手順に相当する。
【００６１】
ステップＳ１１４において、表面粒子間の距離が所定の閾値より長いか否か判断する。そして、長い場合（ステップＳ１１４：Ｙｅｓ）、ステップＳ１１５において、他の粒子の中から候補ｋを選択し、ｋを上記ｉ＋１として上記ステップＳ１１３以降の処理を繰り返す。これらステップＳ１１４（Ｙｅｓ）及びステップＳ１１５からステップＳ１１３へ戻る流れが上記第６の手順乃至第８の手順に相当する。
【００６２】
そして、ステップＳ１１２乃至ステップＳ１１４を全ての境界粒子（ｉ）に対して実行する。これが上記第９の手順に相当する。
【００６３】
（第２の実施の形態）
第２の実施の形態として、表面張力項の計算法について説明する。
本第２の実施の形態における表面張力項の計算では、上述した第１の実施の形態により得られた閉曲面から表面張力を計算する。
【００６４】
上述の第１の実施の形態により、表面の粒子の位置として点列
2014045416.tif 000038
が得られており、２次元の場合、図１２に示すように、表面はｉ、ｉ＋１の粒子間で構成されているとする。３次元の場合は３点間の三角形要素などで表面は構成されているものとする。
【００６５】
流体の表面エネルギーは一般に以下のように定義できる。
【数２１】
2014045416.tif 000039
ここで、上記（式４）の右辺の積分領域は、流体の境界面全域であり、
2014045416.tif 000040
は流体が空気に触れている部分では１、それ以外の部分では０の値をとる。
2014045416.tif 000041
はそれぞれ、図１３に示すように、空気に触れている部分表面張力係数、固体に触れている部分の表面張力係数である。ここで、上記（式４）をどのように粒子法で表現するかが本発明のポイントとなる。
【００６６】
本発明では、界面抽出法により得られた点列
2014045416.tif 000042
を用いて、以下のように表面エネルギーを求める。
【数２２】
2014045416.tif 000043
【００６７】
【数２３】
2014045416.tif 000044
ここで、
2014045416.tif 000045
は、図１４に示すように、（式４）の
2014045416.tif 000046
を幅
2014045416.tif 000047
の間で、０から１までの曲線で滑らかに近似した関数であり、
2014045416.tif 000048
は、３次元上の点
2014045416.tif 000049
と固体との距離を表す関数である。例えば、ｙ＝０の平面に固体壁がある場合は
【数２４】
2014045416.tif 000050
である。
【００６８】
2014045416.tif 000051
は
【数２５】
2014045416.tif 000052
のときは０、
【数２６】
2014045416.tif 000053
のときは１の値を取り、
【数２７】
2014045416.tif 000054
の範囲は何らかの補間（線形補間や３次関数による補間等）を行った関数形を用いる。
【００６９】
積分範囲
2014045416.tif 000055
は、図１５に示すように、線分
【数２８】
2014045416.tif 000056
の上で取られる。積分計算は、ガウス積分法などを用いて数値積分することができる。
（式５）の右辺の第１項の和（ｊ）は、粒子間（ｊ，ｊ＋１の粒子間）で取られ、第２項の和（ｋ）は、全ての粒子で取るものである。右辺第１項は、（式４）の表面エネルギーを離散化したものに対して、右辺第２項は、粒子が壁を通り抜けないという境界条件を表したポテンシャルエネルギーであり、（式６）のような形式を用いると、粒子は固体との距離が０のときポテンシャルは無限大に発散するため、固体の壁を乗り越えて行く非物理的な挙動を抑えることができる。
【００７０】
（式５）のエネルギーを用いて、表面の粒子
2014045416.tif 000057
にかかる力は以下のように計算できる。
【数２９】
2014045416.tif 000058
なお、（式７）は差分近似して求めても良い。
【００７１】
上記（式５）は接触している面が、固体か気体かによって、表面張力係数が変わっている。図１６に示すように、この差によって、接触面が変わる際の接触角が表現できる形式となっている。
【００７２】
【数３０】
2014045416.tif 000059
の場合は固体と接触している場合の方がエネルギー的に低くなるため、固体に付着しようとする力が発生することになる。粒子法では、固体を通り抜けないという条件を壁近傍のポテンシャル（式６)として表すため、粒子が固体に完全に接触することはない。そのため、本第２の実施の形態では、図１７に示すように、（式５）のようにエネルギーを設定して、固体からの距離が
2014045416.tif 000060
以下の場合は反発力を受け、固体からの距離が
2014045416.tif 000061
以上
【数３１】
2014045416.tif 000062
未満の場合は付着力が働くようにした。
【００７３】
図１８は、表面張力項の計算の流れを示すフローチャートである。
まず、ステップＳ１８１において、シミュレーションの対象とする流体のデータを取得し、ステップＳ１８２において、現在の速度を用いて流体粒子の位置を時間ｄｔ／２分だけ更新する。そして、ステップＳ１８３において、物理モデルにより流体粒子にかかる力を算出する。
【００７４】
次に、ステップＳ１８４において、上述の第１の実施の形態により流体の表面粒子のデータを取得する。
【００７５】
そして、ステップＳ１８５において、全ての境界粒子に対して、上述のようにして表面エネルギーを求め、表面粒子にかかる力を計算する。
【００７６】
全ての境界粒子に対する計算が終了したら、ステップＳ１８６において、流体粒子の速度を更新し、ステップＳ１８７において、更新した速度を用いて流体粒子の位置を時間ｄｔ／２分だけ更新する。
【００７７】
そして、上記ステップＳ１８２乃至ステップＳ１８７を繰り返し実行する。
さらに、上記エネルギーに代え、表面粒子の不均一性と流体の体積変化を抑えるため以下の体積エネルギー及び接点エネルギーを導入することも可能である。
【００７８】
体積エネルギー：
【数３２】
2014045416.tif 000063
【００７９】
接点エネルギー：
【数３３】
2014045416.tif 000064
【００８０】
ここで、
2014045416.tif 000065
は流体の体積であり、
【数３４】
2014045416.tif 000066
と計算できる。ここで
2014045416.tif 000067
は外向き法線単位ベクトルであり、
2014045416.tif 000068
は定数であり、
2014045416.tif 000069
は表面粒子ｓ，ｉの質量である。
【００８１】
（式８）、（式９）のエネルギーによる補正項は以下のように計算できる。
【数３５】
2014045416.tif 000070
2014045416.tif 000071
は初期体積で、定数である。
体積エネルギーによる補正項は、図１９のように体積を
2014045416.tif 000072
に戻す効果があり、接点エネルギーによる補正項は、図２０のように表面粒子の不均一性を抑える効果がある。
【００８２】
（第３の実施の形態）
第３の実施の形態として、ＳＰＨ法による非圧縮性粘性流体の運動を計算する際の表面張力導入方法について説明する。
【００８３】
図２１は、２次元流体の例を示す図である。
図２１に示すような状況下（ｙ＝０の平面の上におかれた流体、ｙ方向は鉛直上向きとする）における、非圧縮性粘性流体の運動方程式を考える。
【００８４】
【数３６】
2014045416.tif 000073
【００８５】
【数３７】
2014045416.tif 000074
【００８６】
【数３８】
2014045416.tif 000075
2014045416.tif 000076
は、それぞれ、流体の密度場、速度場、圧力場、音速である。
【００８７】
（式１１）は質量保存則、（式１２）は運動量保存則、（式１３）は状態方程式である。
2014045416.tif 000077
は粘性応力テンソルで、流体の粘性係数を
2014045416.tif 000078
（定数）とすると、
【数３９】
2014045416.tif 000079
である。
【００８８】
（式１１）は流体が集まってくるような速度場がある場合には密度は上昇し、逆に流体が離れていくような速度場の場合には密度は低下する効果を表す。（式１２）の右辺第１項は圧力勾配項で、流体が圧力の大きい部分から圧力の小さい部分へ向かって力が発生する効果を表す。右辺第２項は粘性応力項で、流れにブレーキがかかるような効果をあらわす。右辺第３項は重力項である。
2014045416.tif 000080
は重力加速度、
2014045416.tif 000081
はｙ方向の単位ベクトルである。
【００８９】
本第３の実施の形態では、上記（式１３）のような密度と圧力のみの関係を用いているが、一般の温度や内部エネルギー、エントロピー等を用いた状態方程式を用いても良い。
【００９０】
また、空気との接触部には
2014045416.tif 000082
の表面張力係数（定数）、固体との接触部には
2014045416.tif 000083
の表面張力係数（定数）がかかっているものとする。
【００９１】
（式１１）乃至（式１３）に示したＳＰＨ法を用いて離散化し、上述の第２の実施の形態による表面張力項を導入すると、以下のようになる。
【００９２】
【数４０】
2014045416.tif 000084
【００９３】
【数４１】
2014045416.tif 000085
【００９４】
【数４２】
2014045416.tif 000086
【００９５】
【数４３】
2014045416.tif 000087
【００９６】
ここで、添え字は粒子の番号を表す。すなわちa番目の粒子の位置ベクトル、速度ベクトル、密度、圧力はそれぞれ
2014045416.tif 000088
である。
2014045416.tif 000089
はｂ番目の粒子の質量である。
2014045416.tif 000090
は粘性項を計算するために導入したパラメータで定数である。
【００９７】
（式１５）は、上記（式１２）の運動方程式を粒子法により離散化したもので、右辺第２項は圧力勾配項、第３項は粘性応力項を表したものである。ここで、
2014045416.tif 000091
は
2014045416.tif 000092
のy方向成分、
2014045416.tif 000093
は表面エネルギーであり、上記（式５）と同様の形式となっており、具体的に離散化した以下の形式を採用する。
【００９８】
【数４４】
2014045416.tif 000094
【００９９】
ここで
2014045416.tif 000095
はｙ＝０にある壁からの反発力を表すポテンシャルエネルギーであり、
【数４５】
2014045416.tif 000096
を用いる。
【０１００】
2014045416.tif 000097
は定数である。
2014045416.tif 000098
は上述の第１の実施の形態により抽出された表面粒子であり、粒子群
2014045416.tif 000099
の何れか１つと対応している。
【０１０１】
上記（式１８）の微分
【数４６】
2014045416.tif 000100
は以下のように計算する。
仮に、
2014045416.tif 000101
が表面粒子
2014045416.tif 000102
と対応していない場合は、
【数４７】
2014045416.tif 000103
【数４８】
2014045416.tif 000104
となり、仮に、
2014045416.tif 000105
が表面粒子
2014045416.tif 000106
と対応している場合は、
【０１０２】
【数４９】
2014045416.tif 000107
【０１０３】
【数５０】
2014045416.tif 000108
を計算し、
【数５１】
2014045416.tif 000109
が求められる。
【０１０４】
接点エネルギーは、２次元の場合以下のようになる。
【数５２】
2014045416.tif 000110
【０１０５】
【数５３】
2014045416.tif 000111
ここで
2014045416.tif 000112
，
2014045416.tif 000113
はそれぞれｘ方向単位ベクトル、ｙ方向単位ベクトルである。
【０１０６】
また、流体の体積
2014045416.tif 000114
は
【数５４】
2014045416.tif 000115
と計算できる。ここで
2014045416.tif 000116
はベクトルを
2014045416.tif 000117
回転させる行列である。
【０１０７】
また、体積エネルギー
2014045416.tif 000118
は（式８）のように計算でき、
【数５５】
2014045416.tif 000119
となる。
【０１０８】
さらに、
2014045416.tif 000120
は２次元の再標準化行列であり、
【数５６】
2014045416.tif 000121
【数５７】
2014045416.tif 000122
である。
【０１０９】
ここで、
【数５８】
2014045416.tif 000123
、
【数５９】
2014045416.tif 000124
、
【数６０】
2014045416.tif 000125
、
【数６１】
2014045416.tif 000126
である。
【０１１０】
2014045416.tif 000127
は粒子ａ、ｂ間での１次元リーマン問題を解いて求めた時空間の中間値である。具体的には以下のように決定する。
【０１１１】
粒子ａ、ｂに関して、以下の特性量を定める。
【数６２】
2014045416.tif 000128
【数６３】
2014045416.tif 000129
【数６４】
2014045416.tif 000130
【数６５】
2014045416.tif 000131
【０１１２】
さらに、各勾配を以下のよう計算する。
【数６６】
2014045416.tif 000132
【０１１３】
【数６７】
2014045416.tif 000133
【０１１４】
【数６８】
2014045416.tif 000134
【０１１５】
【数６９】
2014045416.tif 000135
【０１１６】
【数７０】
2014045416.tif 000136
【０１１７】
【数７１】
2014045416.tif 000137
【０１１８】
ここで、上記（式２７）の速度（ｖ）の上付き添え字は、各成分を表す。
これらを用いて、
2014045416.tif 000138
は以下のように決定する。
【０１１９】
【数７２】
2014045416.tif 000139
【０１２０】
【数７３】
2014045416.tif 000140
【０１２１】
【数７４】
2014045416.tif 000141
【０１２２】
【数７５】
2014045416.tif 000142
【０１２３】
【数７６】
2014045416.tif 000143
【０１２４】
図２２及び図２３は、上記の手法により表面張力を加えて固体面との濡れ性を表現したものである。図２２は、
【数７７】
2014045416.tif 000144
で、
【数７８】
2014045416.tif 000145
のときを示し、
【数７９】
2014045416.tif 000146
で、
【数８０】
2014045416.tif 000147
のときを示す。
【０１２５】
これらの図２２及び図２３は、静止しときの流体の形状である。それぞれ、表面張力係数を変化させることで、接触角の違いを表現することができている。
【０１２６】
また、図２４は、（式１９）及び（式２０）により計算される補正項を入れなかった場合を示し、図２５は、補正項を入れた場合を示す。
【０１２７】
図２４では不自然な穴ができ、計算が不安定になるのに対して、図２５では安定に計算できており、補正項の効果が現れている。
【０１２８】
シミュレーションのパラメータは以下の通りである。
なお、鉛直下向き方向をｙ軸とし、ｙ＝０を壁面とした。重力加速度は９．８
2014045416.tif 000148
、粒子数は９００個、初期形状は長さ０．００１６７［ｍ］の格子上に粒子を置き正方形にした。
【０１２９】
ここで、粘性係数
2014045416.tif 000149
は０．０１
2014045416.tif 000150
、
【数８１】
2014045416.tif 000151
、影響半径ｈは０．０００３３３［ｍ］、
【数８２】
2014045416.tif 000152
、密度は１０００
2014045416.tif 000153
、表面張力は
【数８３】
2014045416.tif 000154
、音速ｃは２３．７８
2014045416.tif 000155
、補正項のパラメータは
【数８４】
2014045416.tif 000156
とした。
【０１３０】
上述したように、本発明は粒子法の計算で、表面張力が効果的な流体運動の計算に対して適用可能である。特に溶けた金属の流し込みや、樹脂の流し込みのシミュレーション等には効果的である。
【０１３１】
図５のシミュレーション装置は、例えば、図２６に示すような情報処理装置（コンピュータ）を用いて実現することが可能である。図２６の情報処理装置は、ＣＰＵ（Central Processing Unit）２６０１、メモリ２６０２、入力装置２６０３、出力装置２６０４、外部記憶装置２６０５、媒体駆動装置２６０６及びネットワーク接続装置２６０７を備える。これらはバス２６０８により互いに接続されている。
【０１３２】
メモリ２６０２は、例えば、ＲＯＭ（Read Only Memory）、ＲＡＭ（Random Access Memory）、フラッシュメモリ等の半導体メモリであり、シミュレーション処理に用いられるプログラム及びデータを格納する。例えば、ＣＰＵ２６０１は、メモリ２６０２を利用してプログラムを実行することにより、上述のシミュレーション処理を行う。メモリ２６０２は、図５の格納部５０２としても使用できる。
【０１３３】
入力装置２６０３は、例えば、キーボード、ポインティングデバイス等であり、オペレータからの指示や情報の入力に用いられる。出力装置２６０４は、例えば、表示装置、プリンタ、スピーカ等であり、オペレータへの問い合わせや処理結果の出力に用いられる。出力装置２６０４は、図５の出力部５０３としても使用できる。
【０１３４】
外部記憶装置２６０５は、例えば、磁気ディスク装置、光ディスク装置、光磁気ディスク装置、テープ装置等である。この外部記憶装置２６０５には、ハードディスクドライブも含まれる。情報処理装置は、この外部記憶装置２６０５にプログラム及びデータを格納しておき、それらをメモリ２６０２にロードして使用することができる。
【０１３５】
媒体駆動装置２６０６は、可搬型記録媒体２６０９を駆動し、その記録内容にアクセスする。可搬型記録媒体２６０９は、メモリデバイス、フレキシブルディスク、光ディスク、光磁気ディスク等である。この可搬型記録媒体２６０９には、Compact Disk Read Only
Memory （ＣＤ−ＲＯＭ）、Digital Versatile Disk（ＤＶＤ）、Universal Serial Bus（ＵＳＢ）メモリ等も含まれる。オペレータは、この可搬型記録媒体２６０９にプログラム及びデータを格納しておき、それらをメモリ２６０２にロードして使用することができる。
【０１３６】
このように、シミュレーション処理に用いられるプログラム及びデータを格納するコンピュータ読み取り可能な記録媒体には、メモリ２６０２、外部記憶装置２６０５、及び可搬型記録媒体２６０９のような、物理的な（非一時的な）記録媒体が含まれる。
【０１３７】
ネットワーク接続装置２６０７は、通信ネットワーク２６１０に接続され、通信に伴うデータ変換を行う通信インタフェースである。情報処理装置は、プログラム及びデータを外部の装置からネットワーク接続装置２６０７を介して受け取り、それらをメモリ２６０２にロードして使用することができる。ネットワーク接続装置２６０７は、図５の出力部５０３としても使用できる。
【０１３８】
開示の実施形態とその利点について詳しく説明したが、当業者は、特許請求の範囲に明確に記載した本発明の範囲から逸脱することなく、様々な変更、追加、省略をすることができるであろう。
【図１】
2014045416.tif 000157
【図２】
2014045416.tif 000158
【図３】
2014045416.tif 000159
【図４】
2014045416.tif 000160
【図５】
2014045416.tif 000161
【図６】
2014045416.tif 000162
【図７】
2014045416.tif 000163
【図８】
2014045416.tif 000164
【図９】
2014045416.tif 000165
【図１０】
2014045416.tif 000166
【図１１】
2014045416.tif 000167
【図１２】
2014045416.tif 000168
【図１３】
2014045416.tif 000169
【図１４】
2014045416.tif 000170
【図１５】
2014045416.tif 000171
【図１６】
2014045416.tif 000172
【図１７】
2014045416.tif 000173
【図１８】
2014045416.tif 000174
【図１９】
2014045416.tif 000175
【図２０】
2014045416.tif 000176
【図２１】
2014045416.tif 000177
【図２２】
2014045416.tif 000178
【図２３】
2014045416.tif 000179
【図２４】
2014045416.tif 000180
【図２５】
2014045416.tif 000181
【図２６】
2014045416.tif 000182
【手続補正書】
【提出日】平成27年2月18日(2015.2.18)
【手続補正１】
【補正対象書類名】特許請求の範囲
【補正対象項目名】全文
【補正方法】変更
【補正の内容】
【特許請求の範囲】
【請求項１】
流体の表面張力をシミュレーションするシミュレーションプログラムにおいて、
コンピュータに、
入力された境界条件及び初期条件に基づいて、前記流体を粒子の集まりで表現した流体モデルの界面を、凸包構成アルゴリズムを基に算出させ、
前記算出された界面の表面エネルギーを算出させ、
前記算出された表面エネルギーに基づいて、前記界面の表面張力を算出させ、
前記算出された表面張力に基づいて、所定時間間隔毎の前記流体の状態を出力させることを特徴とするシミュレーションプログラム。
【請求項２】
前記界面の表面エネルギーの算出は、
空気との界面の表面エネルギーを算出する第１項と前記空気以外との界面の表面エネルギーを算出する第２項との和で表現した計算式を用いることを特徴とする請求項１記載のシミュレーションプログラム。
【請求項３】
前記空気との界面の表面エネルギーの算出は、
前記第１項の表面張力係数と前記第２項の表面張力係数とを変えることを特徴とする請求項２記載のシミュレーションプログラム。
【請求項４】
前記空気との界面の表面エネルギーの算出は、前記粒子の表面分布の不均一性を補正するためのエネルギーを加算することを特徴とする請求項３記載のシミュレーションプログラム。
【請求項５】
前記空気との界面の表面エネルギーの算出は、前記流体の体積を補正するためのエネルギーを加算することを特徴とする請求項３記載のシミュレーションプログラム。
【請求項６】
流体の表面張力をシミュレーションするシミュレーション方法において、
コンピュータが、
入力された境界条件及び初期条件に基づいて、前記流体を粒子の集まりで表現した流体モデルの界面を、凸包構成アルゴリズムを基に算出し、
前記算出された界面の表面エネルギーを算出し、
前記算出された表面エネルギーに基づいて、前記界面の表面張力を算出し、
前記算出された表面張力に基づいて、所定時間間隔毎の前記流体の状態を出力することを特徴とするシミュレーション方法。
【請求項７】
流体の表面張力をシミュレーションするシミュレーション装置において、
境界条件及び初期条件を含む諸条件を入力する入力部と、
前記入力された境界条件及び初期条件に基づいて、前記流体を粒子の集まりで表現した流体モデルの界面を算出する凸包構成アルゴリズムを基に界面算出部と、
前記算出された界面の表面エネルギーを算出する表面エネルギー算出部と、
前記算出された表面エネルギーに基づいて、前記界面の表面張力を算出する表面張力算出部と、
前記算出された表面張力に基づいて、所定時間間隔毎の前記流体の状態を出力する出力部を備えることを特徴とするシミュレーション装置。
【国際調査報告】
2014045416.tif 000183
2014045416.tif 000184
2014045416.tif 000185
2014045416.tif 000186