(57)【特許請求の範囲】
【請求項１】
流体を、連続体ではなく、各微小領域を仮想的な粒子とし、各粒子に働く力の相互作用を表現する粒子間の距離をパラメータとした関数に基づいて、各粒子の位置の経時的変化を決定し、流体の形状の変化を再現するものであって、
上記関数は、各粒子間の距離がある臨界点距離となるとお互いに斥力を生じ、この臨界点距離を越えて一定の距離になるとお互いに引力を生じることを特徴とする曲線又は直線からなり、この曲線又は直線は、各粒子の温度変化、速度変化、又は時間変化に伴い、曲線又は直線自体が変化することにより、非ニュートン性粘度の変化を表現し、流体の形状の変化を再現するようにした粒子型流体シミュレーション方法。
【請求項２】
前記関数に基づいて、各粒子の圧力の経時的変化をも決定し、流体の形状の変化を再現する請求項１に記載の粒子型流体シミュレーション方法。
【請求項３】
流体を、連続体ではなく、各微小領域を仮想的な粒子とし、各粒子と固体表面との関係において親水性及び撥水性を再現するために、各粒子と固体表面との力の相互作用を固体表面と粒子間の距離をパラメータとした関数に基づいて、各粒子の位置の経時的変化を決定し、流体と固体表面との親水性及び撥水性等の現象を再現する粒子型流体シミュレーション方法。
【請求項４】
前記関数に基づいて、各粒子の圧力の経時的変化をも決定し、流体と固体表面との親水性及び撥水性等の現象を再現する請求項３に記載の粒子型流体シミュレーション方法。
【請求項５】
前記関数は、各粒子と固体表面との間の距離がある臨界点距離となるとお互いに斥力を生じ、この臨界点距離を越えて一定の距離になるとお互いに引力を生じることを特徴とする曲線又は直線からなり、この曲線又は直線は、各粒子の温度変化、速度変化、又は時間変化に伴い、曲線又は直線自体が変化することにより、非ニュートン性粘度の変化を表現し、流体と固体表面との親水性及び撥水性等の現象を再現するようにした請求項３又は４に記載の粒子型流体シミュレーション方法。
【発明の詳細な説明】
【０００１】
【発明の属する技術分野】
本発明は、流体（液体又は気体）の形状の変化、又は流体から固体までの表面形状（界面状態）を定量的に表現し、動的及び静的な流体挙動を予測する粒子型流体シミュレーション方法及び装置に関するものである。
【０００２】
【従来の技術】
従来の流体解析技術は、２種類の手法に分類でき、１つは、流体の基本方程式であるナビエ・ストークスの方程式（Ｎａｖｉｅｒ−Ｓｔｏｋｅｓの式）を解くものと、もう１つは、自由表面を扱うＶＯＦ法（Ａ  Ｖｏｌｕｍｅ  ｏｆ  Ｆｌｕｉｄ  Ａｌｇｏｒｉｔｈｍ）である。前者は、流体の表面形状には触れておらず、流体内部の流速や圧力や流体抵抗について、計算するものであり、後者は、流体の表面形状に関して計算可能だが、正確な計算ができていないといった問題点が存在する。
【０００３】
具体的な説明を図１９〜図２１で行う。まず、前者の手法では、液体そのものが経時的に変化しても、そのような変化については計算していない。図１９（Ａ）から（Ｂ）に示すように液体６２を入れた容器６０を斜めに傾けた際に容器６０の傾斜に応じて液体全体の形状が変化するにもかかわらず、この種の問題については、取り扱えない問題が存在する。図中、６３は容器６０内の液体６２の界面である。ただし、図２０のように、液体が管全体の中を流れる速さｖや流跡や圧力差ΔＰや流体抵抗を解く計算は行えるが、流体の表面について形状を予測したり、液体の分離及び合体のような表面形状そのものに関しては、取り扱うことができない。後者の手法についても、図２１を利用して説明すると、図２１（Ｂ）のように分割された領域の中で体積の占有率で液体７２の形状を把握するように計算する方法では、容器７０内の液体７２の境界（液体７２の界面７１）については、図２１の（Ｃ）及び（Ｄ）に示すように分割された領域において体積占有率が同一であれば、どのような形状であるか正確に判断できない問題点が存在する。従って、液体７２の表面形状について正確に計算することができない。同様に、液体の分離及び合体のような表面形状の予測も取り扱うことができない。
【０００４】
しかし、近年、様々な液体を利用して、液体を固体表面上に正確かつ均質に膜を形成する技術が、半導体のみならず様々な技術分野に必要不可欠となっている。このとき、液体の表面形状がどのようになっているかを予測することは、非常に重要度を増してきている。例えば、図２２（Ａ）のように圧力Ｐの作用により容器９１から吐出される液体９２の全体の流速ｖのみならず、液ヌレやハジキに伴う容器９１の固体表面上と液体表面との関係で液体９２が横飛び９４を起こす現象は、従来の解析手法では、前述した理由で予測困難であるという問題点が存在する。図２２（Ｂ）において、９３は液ヌレにより液滴が容器９１の開口端面に付着した状態を示している。
【０００５】
【発明が解決しようとする課題】
従来の流体解析技術やＶＯＦ法に代表される自由表面解析技術においては、前述の具体例で説明したように、ミクロな液滴形状そのものの算出や液体と固体表面形状とのヌレ及びハジキに伴う液の横飛びした現象は、計算不可能であった。そのため、このような現象を再現するためには、試行錯誤の実験を繰り返す他なかった。
従って、本発明の目的は、従来の流体解析技術やＶＯＦ法に代表される自由表面解析技術が有していた上記の問題点を解決することにあって、流体のミクロな表面形状そのものの算出やヌレ及びハジキに伴う現象などをより正確に計算することが可能な粒子型流体シミュレーション方法及び装置を提供することである。
【０００６】
【課題を解決するための手段及びその作用効果】
上記目的を達成するために、本発明は以下のように構成する。
すなわち、本発明の第１態様にかかる粒子型流体シミュレーション方法によれば、流体を、連続体ではなく、各微小領域を仮想的な粒子とし、各粒子に働く力の相互作用を表現する粒子間の距離をパラメータとした関数に基づいて、各粒子の位置の経時的変化を決定し、流体の形状の変化を再現するものであって、
上記関数は、各粒子間の距離がある臨界点距離となるとお互いに斥力を生じ、この臨界点距離を越えて一定の距離になるとお互いに引力を生じることを特徴とする曲線又は直線からなり、この曲線又は直線は、各粒子の温度変化、速度変化、又は時間変化に伴い、曲線又は直線自体が変化することにより、非ニュートン性粘度の変化を表現し、流体の形状の変化を再現するようにしている。
本発明の第２態様にかかる発明では、第１態様において、前記関数に基づいて、各粒子の圧力の経時的変化をも決定し、流体の形状の変化を再現することもできる。
【０００７】
本発明の第３態様にかかる粒子型流体シミュレーション方法によれば、流体を、連続体ではなく、各微小領域を仮想的な粒子とし、各粒子と固体表面との関係において親水性及び撥水性を再現するために、各粒子と固体表面との力の相互作用を固体表面と粒子間の距離をパラメータとした関数に基づいて、各粒子の位置の経時的変化を決定し、流体と固体表面との親水性及び撥水性等の現象を再現するようにしている。
本発明の第４態様によれば、第３態様において、前記関数に基づいて、各粒子の圧力の経時的変化をも決定し、流体と固体表面との親水性及び撥水性等の現象を再現することもできる。
本発明の第５態様によれば、第３又は４態様において、前記関数は、各粒子と固体表面との間の距離がある臨界点距離となるとお互いに斥力を生じ、この臨界点距離を越えて一定の距離になるとお互いに引力を生じることを特徴とする曲線又は直線からなり、この曲線又は直線は、各粒子の温度変化、速度変化、又は時間変化に伴い、曲線又は直線自体が変化することにより、非ニュートン性粘度の変化を表現し、流体と固体表面との親水性及び撥水性等の現象を再現するようにすることもできる。
【００１０】
従って、本発明の上記種々の態様によれば、上記流体を仮想的な粒子に置き換えて表現し、各粒子間に働く力の相互作用を粒子間の距離をパラメータとした関数に基づいて、各粒子の位置、速度、密度、圧力、又は温度の経時的変化を決定することができる。
また、本発明の上記種々の態様によれば、上記粒子と上記固体表面との力の相互作用を固体表面と粒子間の距離をパラメータとした関数に基づいて、各粒子の位置、速度、密度、圧力、温度の経時的変化を決定することにより、上記流体と上記固体表面とのヌレ（親水性）及びハジキ現象（撥水性）を再現することができる。
また、各粒子間に働く力の相互作用を温度変化、速度変化、又は時間的変化に伴い各粒子間の距離をパラメータとした関数を変更制御することにより、非ニュートン性粘度のチクソ性等による時間的変化に伴い、粘度が変化する物性に対応することができる。
また、各粒子の質量、体積、温度、粘度、及び表面張力等のパラメータを複数持つことにより、２種類以上の異なる物性からなる流体、又は、２相流（固体及び液体、固体及び気体、液体及び気体、固体及び液体及び気体）以上の状態からなる物体の状態を再現することができる。
よって、本発明によれば、流体（気体又は液体）から固体までの表面形状を粒子型モデルで粒子間力を計算することで、流体の粘性や表面張力による濡れ性やミクロな表面形状などをより正確にシミュレーションすることができ、流体（液体又は気体）の形状の変化、又は流体から固体までの表面形状（界面状態）を定量的に表現し、動的及び静的な流体挙動を予測することができる。よって、例えば、液体塗布装置など流体吐出時の流体収納容器又は吐出ノズルなどの最適な形状や吐出流体の最適粘度を決定することができる。
【００１１】
【発明の実施の形態】
以下に、本発明にかかる実施形態を図面に基づいて詳細に説明する。
本出願は、本出願人が既に出願した特許出願平成８−１７８８９４号（発明の名称「液体及び濡れの評価装置」）は、流体（気体又は液体）を、連続体ではなく、各微小領域を仮想的な粒子として扱う手法について記述したものであり、具体的な計算方法に関してまでは詳細に記載したしたものではなかった。そこで、本発明者は、鋭意研究の結果、流体を粒子で表現した場合の具体的な計算のアルゴリズムを明確にすることができる発明を創作するに至ったものである。その内容は大略２つの点からなり、１点目は粒子間に働く力の相互作用を粒子間の距離をパラメータとした関数に基づき決定することと、２点目は、粒子と固体表面にかかる相互作用においても粒子と固体表面の最短距離をパラメータとした関数に基づき決定することで液体と固体表面との親水性及び撥水性の性質を再現する手法について述べるものである。
【００１２】
言い換えれば、本実施形態では、流体を、連続体ではなく、各微小領域を仮想的な粒子として扱い、これらの粒子のうちの各粒子に働く力の相互作用を粒子間の距離をパラメータとした関数に基づいて、各粒子の位置（これ以外に、例えば、速度、密度、圧力、温度など）の経時的変化を決定し、流体の形状の変化を再現するものである。そして、各粒子に働く力の相互作用を粒子間の距離をパラメータとした関数に基づいて、上記各粒子の圧力の経時的変化をも決定し、流体の形状の変化を再現することもできる。さらに、上記各粒子間に働く力の相互作用を表現する上記距離のパラメータの関数は、各粒子間の距離がある臨界点距離となるとお互いに斥力を生じ、この臨界点距離を越えて一定の距離になるとお互いに引力を生じることを特徴とする曲線又は直線からなり、この曲線又は直線は、各粒子の温度変化、速度変化、又は時間変化に伴い、曲線又は直線自体が変化することにより、非ニュートン性粘度のチクソ性による時間的変化に伴う粘度の変化を表現することもできる。
【００１３】
また、本実施形態では、さらに、流体を、連続体ではなく、各微小領域を仮想的な粒子として扱い、各粒子と固体表面との関係において親水性及び撥水性を再現するために、各粒子と固体表面との力の相互作用を固体表面と粒子間の距離をパラメータとした関数に基づいて、各粒子の位置（これ以外に、例えば、速度、密度、圧力、温度など）の経時的変化を決定し、流体と固体表面との親水性及び撥水性等の現象を再現するものである。そして、各粒子と固体表面との力の相互作用を固体表面と粒子間の距離をパラメータとした関数に基づいて、上記各粒子の圧力の経時的変化をも決定し、流体と固体表面との親水性及び撥水性等の現象を再現することもできる。さらに、上記各粒子と固体表面との間に働く力の相互作用を表現する上記距離のパラメータの関数は、各粒子と固体表面との間の距離がある臨界点距離となるとお互いに斥力を生じ、この臨界点距離を越えて一定の距離になるとお互いに引力を生じることを特徴とする曲線又は直線からなり、この曲線又は直線は、各粒子の温度変化、速度変化、又は時間変化に伴い、曲線又は直線自体が変化することにより、非ニュートン性粘度のチクソ性による時間的変化に伴う粘度の変化を表現することもできる。
また、上記各粒子には、各粒子の質量、体積、粘度、表面張力の粒子パラメータを複数個持つようにして、種々のモデルを表現することができるようにすることもできる。
【００１４】
以下、本発明の一実施形態にかかる粒子型流体シミュレーション方法及び装置について図面を参照しながら説明する。
図１は、上記本発明の上記実施形態にかかる粒子型流体シミュレーション方法を実施するためのフローチャートである。図２は、上記本発明の上記実施形態にかかる粒子型流体シミュレーション方法を実施するための上記シミュレーション装置のブロック図である。
図２において、キーボード又はマウスなどの入力部１から、解析すべき流体及び境界条件などの解析条件などを入力して記憶部２に記憶させたのち、演算部３で記憶部２に記憶されたデータを基に上記シミュレーション方法を実施し、シミュレーション結果を表示画面又は印刷装置などの出力部４で出力する。
【００１５】
上記演算部３で行う上記シミュレーション方法について、図１のフローチャートに基づいて説明する。
まず、ステップＳ１では、初期設定として流体の物性値、又は、流体及び固体表面の物性値を入力部１から記憶部２に設定入力して記憶させるとともに、初期配置として流体の位置及び形状、又は、流体及び固体表面の位置及び形状を同様に入力部１から入力して記憶部２で記憶させる。
次いで、ステップＳ２では、演算部３で、注目する粒子の近傍粒子を検索し、近傍粒子が存在する場合には上記注目粒子とその近傍粒子との中心間の距離計算を行うか、又は、注目する粒子の近傍固体表面を検索し、近傍固体表面が存在すれば注目粒子と近傍固体表面との距離計算を行う。
【００１６】
次いで、ステップＳ３では、上記演算部３での演算結果に基づき、入力部１で、粒子間又は粒子固体表面間の距離に依存したパラメータの選択を行う。粒子間同士の距離に依存したパラメータの選択を行う場合には、ステップＳ４で、液体の種類などに応じて粒子間同士の距離に依存した相互作用のパラメータの参照を行う。例えば、液体Ａと液体Ａとの間、液体Ｂと液体Ｂとの間、又は、液体Ａと液体Ｂとの間など解析すべき対象のそれぞれに応じて、温度変化に伴うパラメータ、速度変化に伴うパラメータ、時間変化に伴うパラメータなどのうちから適当なものを選択する。その後、ステップＳ６に進む。一方、上記演算部３での演算結果に基づき、ステップＳ３で粒子と固体表面との間の距離に依存したパラメータを選択する場合には、ステップＳ５に進む。このステップＳ５では、固体表面に応じて粒子と固体表面との間の距離に依存した相互作用のパラメータの参照を行う。例えば、固体表面Ｃ、又は、固体表面Ｄなど解析すべき対象のそれぞれに応じて、温度変化に伴うパラメータ、速度変化に伴うパラメータ、時間変化に伴うパラメータなどのうちから適当なものを選択する。その後、ステップＳ６に進む。なお、上記液体の種類などに応じて粒子間同士の距離に依存した相互作用の各種パラメータ、及び、上記粒子と上記固体表面との間の各種パラメータは記憶部２に予め記憶されており、適宜、入力部１からの入力操作により記憶部２のこれらのパラメータのデータを参照して選択する。
【００１７】
次いで、ステップＳ６では、上記選択されたパラメータを使用して演算部３で各粒子の受ける力の計算を行う。
次いで、ステップＳ６での計算結果に基づき、ステップＳ７において、演算部３により各粒子の位置の変更を行う。このとき、各粒子の位置以外に、又は、位置は変更することなく、圧力、速度、密度なども必要に応じて変化させる。
次いで、ステップＳ８で、演算部３により解析が終了か否かを判断し、終了でないならば、ステップＳ２に戻る。終了ならば、ステップＳ９で結果を出力部４により適宜出力して、解析を終了する。
【００１８】
ここで、上記ステップＳ６の各粒子の受ける力の計算及びステップＳ７の各粒子の位置の変更について、図３，４に基づいて、さらに、詳細に説明する。
図３，４は、上記実施形態にかかるステップＳ６，Ｓ７の動作において行われるアルゴリズムを記述したものである。ここで、各粒子は、大きく分けて以下の４つの力を受けて、現在の時刻から比べて、次の時刻において各粒子の位置、速度、圧力、温度、及び密度を変化させていると考える。
【００１９】
まず、１つ目の力は、粒子間同士に働く力で、この力により、流体内部の粘性や延性を表すことが可能となる。その計算方法は、粒子ｉと粒子ｊとを考える場合、粒子ｉと粒子ｊの中心間の粒子間距離ｒによって粒子同士間に生じる力Ｆが変化すると考える。これは、図５に示すように、一般に、２つの分子間では分子間力が作用することに鑑みて、粒子同士間においても粒子同士の距離をパラメータとした相互作用が働くものとして近似的に考えるものである。具体的には、縦軸を上記粒子同士間に生じる力Ｆとし、横軸を上記粒子間距離ｒｘとした場合の両者の関係を図５にグラフで示す。このグラフが横軸を横切る点が以下に述べる臨海点距離ｒｃである。よって、このグラフにおいて、粒子ｉと粒子ｊとの中心間の距離である粒子間距離ｒが、図５（Ａ）に示すように、ある臨界点距離ｒｃより短くなると互いに斥力が生じ、図５（Ｂ）に示すように、粒子間距離ｒがある臨海点距離ｒｃに等しくなると粒子間には力が働かなくなり、図５（Ｃ）に示すように、粒子間距離ｒがある臨界点距離ｒｃより長くなると粒子間には互いに引き合う引力が生じる。さらに、図５には具体的には図示しないが粒子間距離がさらに大きくなると粒子間では互いに力を及ぼさなくなることがわかる。また、引力を表す面積が大きいほど基本的には粘度が大きくなり、非常に粘度が強くなると固体のように固まった物性になる。よって、この引力を表わす面積を変化させることにより、粘度を変化させることができる。また、臨界点の距離ｒｃを変えることで、粒子の密度が変化することにもなる。すなわち、臨界点距離ｒｃが大きくなれば粒子間の安定する距離が長くなるため粒子の密度が低くなる。これにより、粒子間の距離を変化させることにより、粒子の密度を変化させることが可能となる。
【００２０】
２つ目の力は、粒子と固体表面間に働く力で、この力により液体の親水性や撥水性を表現することが可能となる。その計算方法は、粒子間の力と同様に粒子の中心と固体表面との間の最短距離を計算し、その距離をパラメータとした関数として表現する。粒子同士の計算と基本的には、同様で、粒子と固体表面との相性によりその関数は変化する。
３つ目の力は、重力であり、全ての粒子に対して常に一定方向に重力が作用している。
４つ目の力は、その他の外力として、そのときにより、静電気力や慣性力など場面に応じて粒子に加えられる力である。
【００２１】
以上の４つの力を考慮して、上記ステップＳ６で、各粒子に加わるこれらの４つの力を全て計算した後に、次の時刻の各粒子の位置及び速度の計算を行い、計算結果に基づき、各粒子の位置及び速度などをステップＳ７で変更する。これを以下に具体的に説明する。
まず、１ｃｃ中に含まれる液体の粒子がＮ個であると仮定し、ｉを１〜Ｎまでの整数としたとき、対象とする粒子を粒子ｉで表す。図３において、ステップＳ２１で１番目の粒子について考えるため、ｉ＝１とする。
次いで、ステップＳ２２では、粒子ｉと該粒子ｉの近傍の粒子との中心間の距離ｒを計算する。すなわち、図４（Ａ）に示すように、粒子ｉと粒子ｊ、粒子ｉと粒子ｊ＋１、……、粒子ｉと粒子Ｎの各中心間距離を計算する。
【００２２】
次に、上記求められた中心間距離に基づき、ステップＳ２３において、図４（Ｂ）に示すように、粒子間における相互作用を考慮するため、各粒子から受ける力Ｆ１ｉの計算を「数１」により行う。ここで、粒子ｉが粒子ｊから受ける力をＦｉｊで表し、逆に、粒子ｊが粒子ｉから受ける力をＦｊｉで表す。
【数１】
0003566009.tif 000002
次に、ステップＳ２４では、各粒子ｉの中心から該粒子ｉの近傍の固体表面までの最短距離を計算する。例えば、図４（Ｃ）に示すように、粒子ｉの中心と固体表面Ｌとの最短距離ｒｉＬ、粒子ｉと固体表面Ｌ＋１との最短距離ｒｉＬ＋１を求める。
【００２３】
次に、上記求められた最短距離に基づき、ステップＳ２５において、各粒子ｉが固体表面から受ける力Ｆ２ｉを計算する。例えば、図４（Ｄ）に示すように、粒子ｉが固体表面Ｌから受ける力ＦｉＬ、粒子ｉが固体表面Ｌ＋１から受ける力ＦｉＬ＋１などを求める。
次に、ステップＳ２６では、各粒子ｉが重力により受ける力Ｆ３ｉ（＝ｍｇ）を計算する。ここで、１粒子の質量をｍｇとする。
さらに、ステップＳ２７では、各粒子ｉがその他の外力より受ける力Ｆ４ｉを計算する。
次に、ステップＳ２８では、各粒子ｉが受ける上記求められた４つの力の総和（Ｆｉ＝Ｆ１ｉ＋Ｆ２ｉ＋Ｆ３ｉ＋Ｆ４ｉ）を計算する。
次に、ステップＳ２９では、ｉをｉ＋１として、次の粒子を対象として考える。すなわち、ステップＳ２９の後、ステップＳ３０でｉ＜Ｎか否か、言い換えれば、総ての粒子についての計算が終了したか否か判断し、終了していないならば、ステップＳ２２に戻り、その粒子についてステップＳ２２〜Ｓ２８の計算を行う。終了しているならば、ステップＳ３１で、各粒子のΔｔ秒後の速度ｖｉを計算する。すなわち、現在の速度をｖｉ０とすると、ｖｉ＝ｖｉ０＋（Ｆｉ／ｍ）×Δｔを計算する。
次に、ステップＳ３２で各粒子のΔｔ秒後の位置を計算する。以上により、ステップＳ６の計算を終了し、その計算結果に基づき、図１のステップＳ７で各粒子の位置を変更する。
【００２４】
一方、図６（Ａ），（Ｂ）は、粒子の密度計算を行う際の考え方を示している。密度は、ある領域（体積Ｖ）内に存在する存在粒子数（Ｎ）により決定される。よって、ステップＳ５０において、ある領域（Ｖ）内の粒子２０の粒子数（Ｎ）を計算する。一方、圧力は容器中に平板を置き、その平板に加わる粒子の力の総和を計算し、平板面積で割られたものを圧力とする。よって、ステップＳ５１において、上記領域（Ｖ）内のＮ個の粒子２０の質量ｍの合計（ｍ×Ｎ）を計算する。この結果、ステップＳ５２において、密度ρは、ρ＝（ｍＮ／Ｖ）により求めることができる。
また、各粒子の温度は、内部エネルギーを温度エネルギーに変換し、外から受ける温度及び外へ逃げる温度の計算は、従来通りの熱伝導計算、又は、熱伝達計算と同様な計算方法で算出することができる。
【００２５】
図７，８は、液体を貯える同一容器から同一初期条件で、粘度の異なる別々の液体が上記容器から流れ出るときの違いを示したものである。図７が粘度の低い液体を流した図であり、図８が粘度の高い液体を流した図である。図７，８において、８１は液滴、８２は液体を貯える容器、８３は床面を示す。液体が吐出される際に、容器８２の出口部に液体が付着して液体が溜まる。その後、液体の表面張力が液体の自重に耐えられなくなり、下方の床面に液体が落下していく。このとき、液体の各粒子に加わる力Ｆは、「数２」で表現される。
【数２】
Ｆ＝Ｆ１＋Ｆ２＋Ｆ３
ここで、Ｆ１は重力、Ｆ２は粒子間力、Ｆ３は固体表面から受ける力である。
【００２６】
粒子ｉと粒子ｊ間の相互作用による力をｆｉｊとすると、「数３」のように表される。
【数３】
ｆｉｊ  ＝  ｆ（ｒｉｊ）
ここで、ｆｉｊは、粒子間距離ｒｉｊの関数として表現される。
例えば、図１０に示すように、ある粒子間に働く曲線が２次式からなる関数で表現されるとすると、「数４」のように表される。
【数４】
ｒｉｊ＜ｒ１  のとき          ｆｉｊ  ＝  −ｕ（ｒｉｊ−ｒ１）（ｒｉｊ−ｒ２）
ｒ１≦ｒｉｊ＜ｒ２  のとき    ｆｉｊ  ＝    ｕ（ｒｉｊ−ｒ１）（ｒｉｊ−ｒ２）
ｒ２≦ｒｉｊ  のとき          ｆｉｊ  ＝    ０
ここで、ｆは力の次元、ｒは距離の次元、ｕは距離の次元を掛けることにより力の次元となるような次元を持つパラメータ（係数）である。
【００２７】
このとき、上記したように、粒子間距離が非常に小さくなると（ｆｉｊ＜０のとき）、粒子間には、互いに反発する斥力が働き、一定の距離が保たれると（ｆｉｊ＞０）お互いに引き合う引力が働く。粒子間の距離ｒが、横軸ｒｉｊと交差する点の臨海点距離ｒ１と等しくなると粒子間に働く力の大きさが零となり、安定した状態となる。このとき、臨海点距離ｒ１の大きさを変更すると液体全体の密度や粘度が変更され、異なる物性を示す液体を表現することが可能となる。臨海点距離ｒ１と距離ｒ２の両方を変更したり、パラメータｕもこれに加えて変更することにより、種々のモデルに対応させることが可能である。
【００２８】
図９に、上記パラメータｕを変更したときの液体の状態の変化を示している。例えば、距離ｒ２を固定した状態（例えば距離ｒ２を３．６に固定した状態）で、粒子間の引力に関係するパラメータｕ又は臨海点距離ｒ１又はその両方を変更すると、２次式の関数で表される曲線の高さが大きくなったり、小さくなったりして、高粘度の液体から低粘度の液体の挙動を示すことが可能となり、さらには、固体から気体の性質までを表現することが可能となる。例えば、パラメータｕを大きくすれば、粒子間又は粒子と固体表面との間での引力が大きくなり、粘度が大きくなる。一方、パラメータｕを小さくすれば、粒子間又は粒子と固体表面との間での引力が小さくなり、粘度が小さくなる。具体的には、図９において、中央のＣの粘土状態において、臨海点距離ｒ１を大きくすれば、Ｃの粘土状態からＦの状態のように「ムース状」の粘土となる。逆に、臨海点距離ｒ１を小さくすれば、Ｃの状態からＧの状態のように「縮こまる」状態の粘土となる。また、パラメータｕを小さくすれば、Ｃの粘土状態からＢの液体状態に変化し、さらに、パラメータｕを小さくすれば、Ａの気体状態となる。一方、逆に、パラメータｕを大きくすれば、Ｃの粘土状態からさらに粘度が大きくなったＤの粘土状態となり、途中で固まって落下しないようになる。さらにパラメータｕを大きくすれば、最初に小さく縮まってしまい、小さな球になった後にまとまって落下するようなＥの固体状態となる。
【００２９】
同様に、粒子が固体表面から、受ける力Ｆ３においても、液と固体表面との相性や固体表面の粗さ及び加工精度の違いにおいて、粒子と固体表面との距離ｒｉｋで親水性や撥水性の表面形状を表現することが可能となる。
粒子ｉと固体表面ｋ間の相互作用による力をｆｉｋとすると、「数５」のように表される。
【数５】
ｆｉｋ  ＝  ｆ（ｒｉｋ）
すなわち、ｆｉｋは粒子と固体表面間距離ｒｉｋの関数として表現される。
例えば、図１１に示すように、ある粒子と固体表面間に働く曲線を２次式からなる関数で表現されるとすると、「数６」のように表される。
【数６】
ｒｉｋ＜ｒ１  のとき          ｆｉｋ＝−ｕ（ｒｉｋ−ｒ１）（ｒｉｋ−ｒ２）
ｒ１≦ｒｉｋ＜ｒ２  のとき    ｆｉｋ＝  ｕ（ｒｉｋ−ｒ１）（ｒｉｋ−ｒ２）
ｒ２≦ｒｉｋ  のとき          ｆｉｋ＝  ０
このとき、ｆｉｋ＞０となり粒子と固体表面との引き合う引力が強くなると粒子が表面に付着し易くなり、ｆｉｋ＜０となり弱くなると斥力が発生してハジキ易くなる。
【００３０】
このように、本実施形態によれば、粒子間に働く相互作用を粒子間の距離に関するパラメータｕで、液体の物性である粘度や表面張力や密度等の値を表現することが可能となる。
以下に、具体例を説明する。
図１２は、ある液滴を固体表面上に落とした際の液体の表面形状を表現している。図１２の左側と右側では、同一物性の液体を使用しているが、粒子と固体表面間に働く力の関数を変更することで、左側に向かうに従いヌレ易い状態を示し、右側に向かいに従いヌレにくい状態を表現している。これは、同一の液体においてもガラス表面に液滴を落としたときと、ポリエチレン系の素材に落としたときの違いを表現しており、前者はヌレ易く、後者はヌレにくいことになる。
一方、図１３においては、図１２の液体が経時的変化で、ｔ＝１の状態からｔ＝２０の最終状態になるまでの変化を示しており、時間的変化が液体の粘度と関係していることを示すものである。
【００３１】
さらに、図１４は、様々な液の流速変化等に伴う液の粘度の変化を表した表である。同一の液体であっても液体の内部温度や液体の流速によって、粘度が変化する性質は、一般に、よく知られている。例えば、図１４に示すように、粘度ηとせん断速度との関係において、せん断速度にかかわらず粘度が一定なものはニュートン性があり、せん断速度が大きくなるに従い粘度が小さくなるものは非ニュートン性でかつチクソトロピーを示し、粘度が大きくなるものは非ニュートン性でかつレオペキシーを示す。これと同様なことを再現するために、あらかじめ液体の種類により、流速変化や温度変化や時間変化に伴い異なる粒子間距離の関数を持ち、各粒子の流速変化や温度変化や時間的変化を観測制御することにより、通常と異なる粒子間の相互作用を行う。これにより、時間的変化に伴い粘度が変更し、実行したものが、図１５，１６である。
【００３２】
図１５において１４０、１４１、１４２、図１６の１４３、１４４、１４５の順に進むほど、時間経過に伴い液体の粘度が変化しており、図１５の１４０及び１４１の状態では、液体が高粘度状態であり、液体が液体容器の固体表面に付着していることを示す。１４２及び１４３の状態では、液体が比較的低粘度状態であり、液体が固体表面に付着しておらず、順次落下することを示す。さらに、図１６の１４４及び１４５の状態では、液体が再度高粘度状態となり、液体が固体表面に付着していることを示す。
さらに、図１７は、毛細管現象を表しており、液体を溜める容器内に細い管を挿入したとき、容器内の液体がその細い管を伝わって毛細管現象により上昇する状態を示している。
また、図１８は、あらかじめ液体を固体表面上に落として置き、その後、落とす位置をずらせて、さらに液体を固体表面上に落とした状態を示している。液体が、合体及び分割の状態を頻繁に起こす液体においては、従来の流体解析手法では計算不可能であった、界面の状態変化を計算することが容易となる。
【００３３】
図２３は、粒子数４００個が１ｃｃとしたときのある液体の高さ比（約１ｃｃの液体を数滴落下させたときの各液滴の高さのうちの最大高さで各液滴の高さを割った値）又は直径の比（各液滴の直径のうち最大直径で各液滴の直径を割った値）と体積比（各液滴の体積のうちの最大体積で各液滴の体積を割った値）との関係を示すグラフである。図２４はその結果を示す説明図である。図２３中、直径比と体積比について、一点鎖線Ｉは上記実施形態によるシミュレーションしたときの関係を示すグラフ、実線ＩＩは実験により測定して求められた関係を示すグラフである。また、高さ比と体積比にいて、太線ＩＩＩは上記実施形態によるシミュレーションしたときの関係を示すグラフ、実線ＩＩは実験により測定して求められた関係を示すグラフである。図２４において、（Ａ）は粒子数が４００個で液体約１ｃｃに相当する場合であり、実験により測定する実際の約１ｃｃの液体の状態を（Ｂ）に示す。同様に、（Ｃ）は粒子数が２００個で液体約０．５ｃｃに相当する場合であり、実験により測定する実際の約０．５ｃｃの液体の状態を（Ｄ）に示す。同様に、（Ｅ）は粒子数が４０個で液体約０．１ｃｃに相当する場合であり、実験により測定する実際の約０．１ｃｃの液体の状態を（Ｆ）に示す。図２３，２４より、実験により求められた関係とシミュレーションによる関係とが近似していることがわかる。
【図面の簡単な説明】
【図１】本発明の一実施形態にかかる粒子型流体シミュレーション方法を示すフローチャートである。
【図２】上記シミュレーション方法を実施するためのシミュレーション装置のブロック図である。
【図３】上記実施形態におけるより具体的な粒子の位置及び速度の計算アルゴリズムを示す図である。
【図４】（Ａ）〜（Ｄ）はそれぞれ図３のアルゴリズム中の距離計算動作の説明図である。
【図５】上記実施形態において、流体の不安定な分子間力を粒子間同士の相互作用で近似する考え方を示す図である。
【図６】（Ａ）及び（Ｂ）は上記実施形態における密度の計算方法を示すフローチャート及びその説明図である。
【図７】上記実施形態における液体のある物性のときの液体の落下状況の比較図である。
【図８】上記実施形態における液体の物性が図７と異なるときの液体の落下状況の比較図である。
【図９】上記実施形態における粒子間の距離に依存したパラメータ値又は距離ｒ１を変更させたときの液体の落下状況の比較図である。
【図１０】上記実施形態における粒子間の相互作用を決定する粒子間距離に依存した関数図である。
【図１１】上記実施形態における粒子と固体表面との相互作用を決定する粒子と固体表面の最短距離に依存した関数図である。
【図１２】上記実施形態における粒子と固体表面との距離に依存したパラメータを変更したときの比較図である。
【図１３】上記実施形態における液体形状の経時的変化図である。
【図１４】上記実施形態における様々な液体の速度変化、時間的変化に伴う液粘度の変化グラフである。
【図１５】上記実施形態における時間的変化に伴い粘度が変更した場合の液の落下図である。
【図１６】図１５に続く、上記実施形態における時間的変化に伴い粘度が変更した場合の液体の落下図である。
【図１７】上記実施形態における毛細管現象図である。
【図１８】上記実施形態における液体の合体図である。
【図１９】（Ａ），（Ｂ）は容器中に入れた液体の容器を傾ける前及び傾けたときの液体界面の経時的変化図である。
【図２０】管内の液体の流れ図である。
【図２１】（Ａ）〜（Ｄ）はＶＯＦ法による計算方法を示し、（Ａ）は液体が容器内に入った状態の斜視図、（Ｂ）は（Ａ）の一部拡大図、（Ｃ），（Ｄ）は同一占有率の比較図である。
【図２２】（Ａ），（Ｂ）は流体の流れと液のヌレによる横飛び現象が生じているときの説明図及びその拡大図である。
【図２３】ある液体を数滴落下させたときの液体の体積比と直径比又は高さ比との関係について、上記実施形態にかかるシミュレーションにより求めた関係と実験により測定した関係とを比較するグラフである。
【図２４】図２３においてシミュレーションにおける粒子数４００、２００、４０個の場合の状態とそれぞれの実際の液体の状態を示す比較図である。
【符号の説明】
１  入力部
２  記憶部
３  演算部
４  出力部
８１  液滴
８２  液を貯える容器
８３  床面
１４０、１４１  高粘度状態で液体が固体表面に付着している状態
１４２、１４３  比較的低粘度の液体に変更したときの付着状態
１４４、１４５  再度高粘度の液体に変更したときの付着状態
【図１】
0003566009.tif 000003
【図２】
0003566009.tif 000004
【図３】
0003566009.tif 000005
【図４】
0003566009.tif 000006
【図５】
0003566009.tif 000007
【図６】
0003566009.tif 000008
【図７】
0003566009.tif 000009
【図８】
0003566009.tif 000010
【図９】
0003566009.tif 000011
【図１０】
0003566009.tif 000012
【図１１】
0003566009.tif 000013
【図１２】
0003566009.tif 000014
【図１３】
0003566009.tif 000015
【図１４】
0003566009.tif 000016
【図１５】
0003566009.tif 000017
【図１６】
0003566009.tif 000018
【図１７】
0003566009.tif 000019
【図１８】
0003566009.tif 000020
【図１９】
0003566009.tif 000021
【図２０】
0003566009.tif 000022
【図２１】
0003566009.tif 000023
【図２２】
0003566009.tif 000024
【図２３】
0003566009.tif 000025
【図２４】
0003566009.tif 000026