流体解析装置、流体解析方法及び流体解析プログラム
【課題】汎用性が高く、濡れ性モデルを統一的に取り扱うことができ、空間解像度が粗くても精度よく解析を行なうことができる流体解析装置、流体解析方法及び流体解析プログラムを提供する。
【解決手段】外力計算手段２１２は、所定範囲に含まれる第１の粒子から受ける外力は（３）式を用いて算出し、第２の粒子から受ける外力は（５）式を用いて算出する。この（５）式は、（３）式で算出した外力に対して、（１＋ｃｏｓθ）／２を乗算した算式である。粒子移動計算手段２１３は、算出したこれら外力の合計値を用いて、各粒子の仮速度を算出し、この仮速度から仮位置を算出する。粒子配置修正手段２１４は、仮速度及び仮位置及びポアソン方程式の計算を行ない、圧力勾配の計算を行なって各粒子の速度と位置を修正し、次のタイムステップの各粒子の位置、速度及び圧力を出力する。このようなタイムステップの処理を、反復数上限値まで繰り返す。
【解決手段】外力計算手段２１２は、所定範囲に含まれる第１の粒子から受ける外力は（３）式を用いて算出し、第２の粒子から受ける外力は（５）式を用いて算出する。この（５）式は、（３）式で算出した外力に対して、（１＋ｃｏｓθ）／２を乗算した算式である。粒子移動計算手段２１３は、算出したこれら外力の合計値を用いて、各粒子の仮速度を算出し、この仮速度から仮位置を算出する。粒子配置修正手段２１４は、仮速度及び仮位置及びポアソン方程式の計算を行ない、圧力勾配の計算を行なって各粒子の速度と位置を修正し、次のタイムステップの各粒子の位置、速度及び圧力を出力する。このようなタイムステップの処理を、反復数上限値まで繰り返す。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、粒子法を用いて連続体である流体の挙動を解析する流体解析装置、流体解析方法及び流体解析プログラムに関する。
【背景技術】
【０００２】
近年、流体の挙動を解析する方法の１つとして粒子法が知られている（例えば、非特許文献１及び非特許文献２参照。）。粒子法は、流体を複数の粒子の集まりとして表わし、これら粒子間の相互作用の計算によって、流体の挙動を解析する。この粒子法は、有限体積法や有限要素法の計算において、手間のかかる格子生成が不要という利点がある。
【０００３】
更に、この粒子法の１つとしてＭＰＳ（moving particle semi-implicit method）法がある。このＭＰＳ法は、連続体の離散的な計算を、粒子間相互作用モデルを用いて行なう。このＭＰＳ法においては、粒子に対して所定の基準半径以内にある近接の粒子が、その粒子に対して相互作用を与えると仮定する。なお、基準半径が小さい程、粒子の数が減少するので、計算時間が短くなるが、小さすぎると計算が不安定になる。そこで、計算時間と計算の安定性の観点から、一般に、基準半径は初期状態の粒子間距離の２〜４倍に設定される。この場合、２次元の座標系においては、基準半径以内にある粒子の数は１２〜４４個程度となる。
【０００４】
このＭＰＳ法を用いて非圧縮流体の挙動を解析する場合には、粒子数密度が一定であることが利用される。ここで、粒子数密度をｎiとし、粒子ｉ及びその近傍の粒子ｊの位置
ベクトルをそれぞれｒi，ｒjとすると、粒子数密度ｎiは、粒子ｉの位置における重み関
数ｗの和（次式）で示される。
ｎi＝Σｗ（│ｒj−ｒi│） ・・・（１）
ここで、各粒子の保持する質量が一定であると仮定すると、粒子数密度は流体の密度と比例する。従って、非圧縮流体の場合には粒子数密度が一定であり、ＭＰＳ法では、この条件を満足するように解析を行なう。
【０００５】
ところで、マイクロ流路における流体の挙動など、長さスケールが小さい流れにおいては、表面張力の影響が大きくなる。そこで、表面張力を考慮して解析を行なうことがある（非特許文献２のｐ．３５−３８参照）。非特許文献２に記載の技術では、表面張力係数σ、曲率κ、表面張力が表面のみに働くためのデルタ関数δ、表面と垂直な方向の単位ベクトルｖを用いて、運動量保存則の方程式に表面張力項を導入して解析を行なう。従って、表面張力を考慮する場合には、この表面張力項が示す部分を、粒子に働く力としてモデル化する必要がある。
【０００６】
従来、このモデル化においては、表面に位置する粒子を抽出し、この粒子に対する曲率計算を行なう。この曲率計算においては、上記（１)式に用いた重み関数とは異なる重み
関数（以下、曲率計算用重み関数という。）を用いた粒子数密度を算出する。この曲率計算用重み関数は、所定領域内では常に値が「１」で、領域外では常に値が「０」に設定された関数である。
【０００７】
そして、表面が平らな場合（すなわち曲率κが「０」の場合）の曲率計算用粒子数密度を、初期粒子配置において予め計算しておく。この曲率計算用粒子数密度は所定領域内にある単なる近傍粒子の数であるから、表面が凸の場合には、初期粒子配置における曲率計算用粒子数密度より小さくなり、表面が凹の場合には大きくなる。従って、曲率計算用粒子数密度の変化分から界面の曲率κを算出し、これを用いて表面張力を算出する。
【０００８】
また、表面張力の他の計算方法として、粒子間にポテンシャル力を働かせて計算する技術もある（非特許文献３を参照。）。この非特許文献３に記載の方法を用いることにより、比較的容易に表面張力の計算を行なうことができる。
【非特許文献１】越塚誠一著，「インテリジェント・エンジニアリング・シリーズ 数値流体力学」，初版，株式会社培風館，１９９７年，ｐ．１５１−１９４
【非特許文献２】越塚誠一著，「計算力学レクチャーシリーズ５ 粒子法」，丸善株式会社，２００５年，ｐ．１−７，ｐ．２１−３８
【非特許文献３】Ｎ．Shirakawa、他３名，「Journal of Nuclear Science and Technology」，２００１年，第３８巻，ｐ．３９２−４０２
【発明の開示】
【発明が解決しようとする課題】
【０００９】
上述した界面の曲率κを算出して表面張力を計算する場合、曲率計算用粒子数密度の計算精度を保つために、十分な粒子数（例えば直径に１０粒以上）が必要である。このため、空間解像度が粗い場合（すなわち計算に用いる粒子数が少ないように粗く設定する場合）には、曲率計算用粒子数密度の計算精度が悪くなり、流体の挙動を的確に解析することができなかった。
【００１０】
ところで、表面張力が影響するモデルとして、壁に対して液体の水滴が接触角θで静止した現象の「濡れ性」のモデルがある。この場合に用いる表面張力項の計算においても、上述したように、表面と垂直な方向の単位ベクトルｖを用いる。このため、「濡れ性」の挙動解析を行なう場合には、接触角θから、表面と垂直な方向のベクトル成分を幾何学的に計算し、その結果を用いていた。この場合、曲率計算用粒子数密度は、上述した表面張力の通常の計算よりも、誤差が大きくなっていた。更に、この「濡れ性」では接触角θをパラメータとして用いるため、解析対象の流体（液体）が流動する場合には、流動状態での接触角θを把握してパラメータとして設定する必要がある。しかし、このような接触角θを正確に把握することは困難である。そのため、ダイナミックな解析を行なうことは難しく、汎用性に乏しかった。
【００１１】
また、非特許文献３に記載の技術の場合、ポテンシャル力と表面張力係数σとの関係が不明であるため、パラメータ調整が必要となり、濡れ性の取り扱いについても考慮されていない。
【００１２】
本発明は、上述した問題に鑑みてなされ、その目的は、汎用性が高く、濡れ性モデルを統一的に取り扱うことができ、空間解像度が粗くても精度よく解析を行なうことができる流体解析装置、流体解析方法及び流体解析プログラムを提供することにある。
【課題を解決するための手段】
【００１３】
上記問題点を解決するために、請求項１に記載の発明は、流体を表した第１の粒子毎に、速度、位置、圧力に関するデータと、前記流体に接触する壁を表した第２の粒子毎に、速度、位置、圧力に関するデータとを記録する粒子データ記憶手段と、前記流体と前記壁との接触角θに関するデータと、各粒子について他の粒子からの作用力を受ける作用範囲に関するデータとを記録した計算条件データ記憶手段と、各粒子の挙動を計算する制御手段とから構成され、流体の挙動を解析する粒子法を用いた解析装置であって、前記制御手段が、所定のタイムステップにおいて、各第１の粒子の配置において前記作用範囲に含まれる他の粒子を特定する特定手段と、特定された他の粒子のうち第１の粒子については、第１の関数により第１の作用力を算出し、特定された他の粒子のうち第２の粒子については、前記第１の関数に対して（１＋ｃｏｓθ）／２を乗算した第２の関数により第２の作用力を算出する作用力算出手段と、前記各作用力により算出した各粒子の仮速度を用いて
、仮配置を決定する粒子移動計算手段と、前記仮配置を用いて、前記粒子の数密度を一定にする計算により、次回タイムステップの圧力を算出し、この圧力を用いて算出した圧力勾配から各粒子の修正後の速度を算出し、この速度を用いて各粒子の位置を修正し、前記粒子データ記憶手段に記録する粒子配置修正手段と、タイムステップを進めて、前記特定手段、前記粒子移動計算手段及び前記粒子配置修正手段を用いた各粒子の挙動の計算処理を繰り返す計算制御手段とを備えたことを要旨とする。
【００１４】
請求項２に記載の発明は、請求項１に記載の流体解析装置において、前記第１の関数は、作用力を及ぼす前記他の粒子との距離に応じたポテンシャルエネルギにより算出する関数であることを要旨とする。
【００１５】
請求項３に記載の発明は、流体を表した第１の粒子毎に、速度、位置、圧力に関するデータと、前記流体に接触する壁を表した第２の粒子毎に、速度、位置、圧力に関するデータとを記録する粒子データ記憶手段と、前記流体と前記壁との接触角θに関するデータと、各粒子について他の粒子からの作用力を受ける作用範囲に関するデータとを記録した計算条件データ記憶手段と、各粒子の挙動を計算する制御手段とから構成された解析装置を用いて、粒子法による流体の挙動を解析する方法であって、前記制御手段が、所定のタイムステップにおいて、各第１の粒子の配置において前記作用範囲に含まれる他の粒子を特定する特定段階と、特定された他の粒子のうち第１の粒子については、第１の関数により第１の作用力を算出し、特定された他の粒子のうち第２の粒子については、前記第１の関数に対して（１＋ｃｏｓθ）／２を乗算した第２の関数により第２の作用力を算出する作用力算出段階と、前記各作用力により算出した各粒子の仮速度を用いて、仮配置を決定する粒子移動計算段階と、前記仮配置を用いて、前記粒子の数密度を一定にする計算により、次回タイムステップの圧力を算出し、この圧力を用いて算出した圧力勾配から各粒子の修正後の速度を算出し、この速度を用いて各粒子の位置を修正し、前記粒子データ記憶手段に記録する粒子配置修正段階と、タイムステップを進めて、前記特定段階、前記粒子移動計算段階及び前記粒子配置修正段階による各粒子の挙動の計算処理を繰り返す計算制御段階とを実行することを要旨とする。
【００１６】
請求項４に記載の発明は、請求項３に記載の流体解析方法において、前記第１の関数は、作用力を及ぼす前記他の粒子との距離に応じたポテンシャルエネルギにより算出する関数であることを要旨とする。
【００１７】
請求項５に記載の発明は、流体を表した第１の粒子毎に、速度、位置、圧力に関するデータと、前記流体に接触する壁を表した第２の粒子毎に、速度、位置、圧力に関するデータとを記録する粒子データ記憶手段と、前記流体と前記壁との接触角θに関するデータと、各粒子について他の粒子からの作用力を受ける作用範囲に関するデータとを記録した計算条件データ記憶手段と、各粒子の挙動を計算する制御手段とから構成された解析装置を用いて、粒子法による流体の挙動を解析するプログラムであって、前記制御手段を、所定のタイムステップにおいて、各第１の粒子の配置において前記作用範囲に含まれる他の粒子を特定する特定手段、特定された他の粒子のうち第１の粒子については、第１の関数により第１の作用力を算出し、特定された他の粒子のうち第２の粒子については、前記第１の関数に対して（１＋ｃｏｓθ）／２を乗算した第２の関数により第２の作用力を算出する作用力算出手段、前記各作用力により算出した各粒子の仮速度を用いて、仮配置を決定する粒子移動計算手段、及び前記仮配置を用いて、前記粒子の数密度を一定にする計算により、次回タイムステップの圧力を算出し、この圧力を用いて算出した圧力勾配から各粒子の修正後の速度を算出し、この速度を用いて各粒子の位置を修正し、前記粒子データ記憶手段に記録する粒子配置修正手段と、タイムステップを進めて、前記特定手段、前記粒子移動計算手段及び前記粒子配置修正手段を用いた各粒子の挙動の計算処理を繰り返す計算制御手段として機能させることを要旨とする。
【００１８】
請求項６に記載の発明は、請求項５に記載の流体解析プログラムにおいて、前記第１の関数は、作用力を及ぼす前記他の粒子との距離に応じたポテンシャルエネルギにより算出する関数であることを要旨とする。
【００１９】
（作用）
本発明によれば、制御手段が、作用範囲に含まれる第１の粒子については、第１の関数により第１の作用力を算出し、作用範囲に含まれる第２の粒子については、第１の関数に対して（１＋ｃｏｓθ）／２を乗算した第２の関数により第２の作用力を算出する。制御手段は、各作用力により算出した各粒子の速度を用いて、仮配置を決定する。更に、制御手段は、仮配置を用いて、粒子の数密度を一定にする計算により、次回タイムステップの圧力を算出し、この圧力を用いて算出した圧力勾配から各粒子の修正速度を算出し、これから各粒子の位置を補正する。以上の処理を、制御手段は、タイムステップを進めて、繰り返して行なう。（１＋ｃｏｓθ）／２は、液相の界面エネルギと壁−液相の界面エネルギの関係式であるヤングの式から算出される係数である。このため、制御手段は、エネルギ関係式から第２の粒子の作用力を算出し、これら第１及び第２の粒子の作用力を用いて流体の挙動の解析を行なうので、濡れ性モデルを、表面張力モデルとして統一的に扱うことができる。
【００２０】
また、制御手段は、粒子数密度の減少分から算出される界面の曲率を用いずに、第２の粒子の作用力を算出しているので、空間解像度が粗くても、流体の挙動の解析を精度よく行なうことができる。
【００２１】
更に、エネルギ関係式は、流体や壁がどのような状態にあっても成立する。従って、制御手段は、第２の粒子の作用力をエネルギ関係式から算出しているので、例えば界面が動く場合など、種々の条件下でも用いることができ、汎用性が高い。
【００２２】
本発明によれば、第１の関数は、作用力を及ぼす前記他の粒子との距離に応じたポテンシャルエネルギにより算出する関数である。このため、ヤングの式から算出される係数と整合がよく、流体の挙動の解析を、より的確に良好に行なうことができる。
【発明の効果】
【００２３】
本発明によれば、粒子法において、汎用性が高く、濡れ性モデルを統一的に取り扱うことができ、空間解像度が粗くても精度よく解析を行なうことができる。
【発明を実施するための最良の形態】
【００２４】
以下、本発明を具体化した流体解析装置の一実施形態を図１〜図５に基づいて説明する。本実施形態では、壁が剛体で構成された収容器内の流体の挙動について解析する場合を想定する。ここでは、流体として、水のように圧縮率が小さい流体を用いる場合を想定する。また、本実施形態では、説明を簡単にするために、流体の流入及び流出がなく、２次元の解析を行なう場合を想定する。
【００２５】
図１に示すように、本実施形態の流体解析装置は、入力手段１０、出力手段１５及び流体解析手段２０を備えている。
入力手段１０は、流体解析手段２０に各種データを入力するための手段であり、例えばキーボードやポインティング・デバイスなどを用いる。そして、入力手段１０は、入力された計算条件データや初期条件における粒子データなどを流体解析手段２０に供給する。
【００２６】
出力手段１５は、流体解析手段２０による解析結果を出力するために用いられ、具体的には、ディスプレイやデータ記憶部などを用いる。
流体解析手段２０は、データの処理を行なうためのコンピュータシステムであって、制御手段２１、計算条件データ記憶部２２及び粒子データ記憶部２３を備える。
【００２７】
制御手段２１は、図示しないＣＰＵ、ＲＡＭ及びＲＯＭ等を有し、後述する処理（特定段階、作用力算出段階、粒子移動計算段階、粒子配置修正段階及び計算制御段階等を含む処理）を行なう。そして、このための流体解析プログラムを実行することにより、制御手段２１は、特定手段としても機能する計算制御手段２１１、作用力算出手段としての外力計算手段２１２、粒子移動計算手段２１３及び粒子配置修正手段２１４として機能する。
【００２８】
計算条件データ記憶部２２には、流体解析に用いる物性や計算制御などの計算条件データが記憶されている。この計算条件データは、流体解析を行なう前に、入力手段１０を介して登録される。本実施形態においては、解析対象の粒子としては、挙動を解析する流体を構成する第１の粒子Ａ１と、壁を構成する第２の粒子Ａ２がある。更に、壁の粒子Ａ２においては、第１の粒子Ａ１と接する壁粒子と、その外側の壁粒子とがある。このため、この計算条件データには、粒子識別子、内側壁粒子の種類識別子、壁粒子の種類識別子に関するデータが含まれる。更に、計算条件データには、密度、重力加速度、動粘性係数、接触角及び表面張力係数に関するデータが含まれる。
【００２９】
粒子識別子データ領域には、この流体解析の対象となる粒子の種類を特定するための識別子に関するデータが記録されている。このデータによって、制御手段２１は、計算に用いるデータのうち、どのデータが粒子の種類を特定するためのデータであるかを特定できる。
【００３０】
内側壁粒子の種類識別子データ領域には、第１の粒子Ａ１と接して圧力計算を行なう内側の壁の粒子の種類を特定するための識別子に関するデータが記録されている。これにより、内側の壁の粒子の種類識別子データと、内側壁粒子であるというフラグとが関連付けられることになる。この種類識別子によって特定される粒子Ａ２について、制御手段２１は、その位置は常に一定とし、この粒子Ａ２についての圧力を算出する。
【００３１】
壁粒子の種類識別子データ領域には、内側の壁の粒子Ａ２の外側に存在して壁を構成する粒子を特定するための識別子に関するデータが記録されている。これにより、制御手段２１は、この種類識別子によって特定される壁の粒子Ａ２について、その位置は常に一定とし、かつこの壁の粒子Ａ２の圧力は算出しない。
【００３２】
密度データ領域には、各粒子識別子に対応させて、各粒子によって構成される各物質の密度ρに関するデータが記録されている。
重力加速度データ領域及び動粘性係数データ領域には、それぞれ、計算に用いる重力加速度ｇ及び流体の動粘性係数νに関するデータが記録されている。
【００３３】
接触角データ領域には、接触角θに関するデータが記録されている。この接触角θは、図３に示すように、粒子Ａ２で表される壁（固体）の表面に、粒子Ａ１で表される流体（液体）が置かれた場合に、静的な状態での接触点において液体がなす角度である。
【００３４】
表面張力係数データ領域には、表面張力係数σに関するデータが記録されている。この表面張力係数σは、後述するポテンシャル関数のポテンシャルパラメータαを算出するために用いられる。
【００３５】
更に、この計算条件データには、反復数上限値、クーラン数上限値、圧力計算条件、粒子間距離、第１半径、第２半径、第３半径及び自由表面条件パラメータに関するデータが含まれる。
【００３６】
反復数上限値データ領域には、流体の挙動を計算する場合のタイムステップの上限値に関するデータが記録されている。制御手段２１は、この上限値に達した場合に計算処理を終了する。
【００３７】
クーラン数上限値データ領域には、クーラン数の上限値Ｃmaxに関するデータが記録さ
れている。クーラン数は、時間刻み幅Δｔを算出するために用いられる。
圧力計算条件データ領域には、圧力計算における収束条件や反復数の上限及び下限に関するデータが記録されている。
【００３８】
粒子間距離データ領域には、初期配置における各粒子と、この近傍の粒子との粒子間距離Ｌ0に関するデータが記録されている。
第１半径データ領域には、粒子数密度を計算するために用いる第１半径Ｒ1に関するデ
ータが記録されている。この第１半径Ｒ1は、粒子数密度n（0）を算出するために用いら
れ、具体的には、経験値から、粒子間距離Ｌ0の 2.1倍の値が用いられる。
【００３９】
第２半径データ領域には、後述する圧力のポアソン方程式の計算におけるラプラシアンモデルにおいて用いられる半径に関するデータが記録されている。具体的には、経験値から、この第２半径データとして、粒子間距離Ｌ0の４倍の値が用いられる。
【００４０】
第３半径データ領域には、後述するポテンシャル関数に用いる第３半径Ｒ3に関するデ
ータが記録されている。本実施形態では、この第３半径Ｒ3データとして、粒子間距離Ｌ0の３倍の値を用いる。この第３半径Ｒ3は、図４に示すように作用力を計算するために用
いる半径として用いられ、この第３半径Ｒ3内が特許請求の範囲の「作用範囲」に相当す
る。
【００４１】
自由表面条件パラメータ領域には、自由表面にある粒子を特定するためのパラメータβが記録されている。制御手段２１は、後述する陽的な計算が終了した時点における粒子数密度ｎ*が、β・n（0）より小さい粒子を自由表面上の粒子と特定する。自由表面の外側
には粒子が配置されないため、自由表面上の粒子の粒子数密度が低下するからである。そして、制御手段２１は、特定した自由表面上の粒子についての圧力を、圧力のポアソン方程式を解く際にＰ（k+1）＝０とする。なお、これにより、ディリクレ境界条件を設定し
たことになる。
【００４２】
粒子データ記憶部２３には、各タイムステップにおける粒子に関するデータが記録される。この粒子データ記憶部２３には、タイムステップ毎に、粒子毎の種類識別子、位置、速度及び圧力に関するデータが記録される。すなわち、粒子についての種類識別子、位置、速度及び圧力に関するデータは、各タイムステップにおいて粒子総数分だけ記憶されることになる。
【００４３】
タイムステップデータ領域には、各タイムステップを特定するためのデータが記録される。
種類識別子データ領域、位置データ領域、速度データ領域及び圧力データ領域には、このタイムステップにおける、それぞれの粒子の種類を特定するための識別子、位置、速度及び圧力に関するデータが記録される。本実施形態では、位置データとして粒子のｘ座標及びｙ座標を用い、速度データとして、ｘ方向速度成分及びｙ方向速度成分を用いる。
【００４４】
次に、流体解析手段２０の制御手段２１が備える各処理手段（２１１〜２１４）について詳述する。
計算制御手段２１１は、入力手段１０を介して、計算条件や粒子の初期配置に関するデ
ータを取得して、計算条件データ記憶部２２及び粒子データ記憶部２３にデータを記憶する。
【００４５】
更に、計算制御手段２１１は、時間刻み幅Δｔを算出する次式を記憶している。
Δｔ＝Ｌ0・Ｃmax／ｕmax ・・・（２）
ここで、時間刻み幅Δｔが大きい程、計算が効率的に行なえるので、時間刻み幅Δｔはできるだけ大きくすることが好ましい。しかし、数値安定性を保つためには、クーラン数の最大値に制限を設ける必要がある。上記（２）式において、時間刻み幅Δｔや粒子間距離Ｌ0は一定であるので、クーラン数が最大になる粒子は、速度の絶対値が最大値をもつ粒
子である。従って、計算制御手段２１１は、粒子Ａ１の最大値の速度と、クーラン数の上限値Ｃmaxと、粒子間距離Ｌ0とを上記（２）式に代入して、時間刻み幅Δｔを算出する。
【００４６】
また、計算制御手段２１１は、上記（１）式を記憶しており、外力計算手段２１２が算出した各粒子Ａ１，Ａ２の仮位置ｒ*に基づく粒子数密度ｎ*を、この（１）式を用いて算出する。
【００４７】
更に、計算制御手段２１１は、タイムステップｋにおいて算出した各粒子Ａ１，Ａ２の位置ｒ（k+1）、速度ｕ（k+1）及び圧力Ｐ（k+1）を、出力手段１５を介して出力する。
そして、計算制御手段２１１は、解析処理を終了するかを判断し、終了しない場合には次のタイムステップにおいて処理を繰り返して実行する。
【００４８】
一方、外力計算手段２１２は、流体の粒子Ａ１のそれぞれに対して、第３半径Ｒ3内に
ある粒子を特定する。
更に、この外力計算手段２１２は、各粒子の相互距離Ｌを算出する式を記憶している。このため、外力計算手段２１２は、この式と、各粒子の位置データとを用いて、各粒子の相互距離Ｌを算出する。
【００４９】
また、外力計算手段２１２は、流体の粒子Ａ１が、他の流体の粒子Ａ１に作用する第１の作用力としての外力ｆllを算出する次式（第１の関数）を記憶している。
ｆll＝α・（Ｌ−Ｌ０）・（Ｌ−Ｒ3） ・・・（３）
この（３）式は、界面付近の粒子のみに対して中心力が作用するとして分子間力を模擬したポテンシャル関数である。ただし、相互距離Ｌが第３半径Ｒ3以上の場合には、ｆll
＝０とする。
【００５０】
更に、外力計算手段２１２は、上記（３）式におけるポテンシャルパラメータαを算出する次式を記憶している。
α＝２・Ｌ０2・σ／Σｐe（Ｌi'j'） ・・・（４）
ここで、この（４）式の導出について説明する。表面張力係数σとは、単位面積の表面を形成するために必要なエネルギである。そこで、ポテンシャルエネルギをｐe〔＝（Ｌ
−Ｌ０）（Ｌ−Ｒ3）〕と表わした場合、ポテンシャル関数が〔α・ｐe〕のときに単位面積の表面を形成するために必要なエネルギが表面張力係数σとなるようにポテンシャルパラメータαを算出する。この算出においても、相互距離Ｌが第３半径Ｒ3以上の場合には
「０」とする。ここで、図５に示すように、２つの領域（第１領域Ｅ１、第２領域Ｅ２）から構成された物質を想定し、この物質から第１領域Ｅ１を削除して、この第１領域Ｅ１に接していた第２領域Ｅ２に界面を形成する場合を検討する。この場合、ポテンシャルパラメータαは、界面を形成していない場合（第１領域Ｅ１が第２領域Ｅ２に接している状態）と界面を形成した場合（第１領域Ｅ１を除去した状態）とのエネルギ差と、表面張力係数σとの関係式から算出される。具体的には、１つの粒子の面積ｄｓ＝Ｌ０２の表面のもつポテンシャルエネルギは、第１領域Ｅ１を構成する粒子ｊ'を、第２領域Ｅ２におい
て面積ｄｓ内に存在する粒子ｉ'（図５における領域ｅ２の粒子）から引き離すために必
要な仕事の半分である。粒子間相対位置ベクトルＬi'j'のポテンシャルエネルギをｐe（
Ｌi'j'）とした場合、以下のように表わされる。
σ＝（１／２）・〔Σｐe（Ｌi'j'）／Ｌ０2〕・α ・・・（４−１）
この（４−１）式を用いることにより、上記（４）を算出することができる。
【００５１】
外力計算手段２１２は、壁の粒子Ａ２が、流体の粒子Ａ１に作用する第２の作用力としての外力ｆlwを算出する次式（第２の関数）を記憶している。
ｆlw＝（１＋ｃｏｓθ）ｆll／２ ・・・（５）
ここで、この（５）式の導出について、液体粒子間のポテンシャル力Ｃfと液体固体間
のポテンシャル力Ｃfsの比とヤングの式との対応関係から説明する。
【００５２】
ヤングの式は、次式で示される。
σs−σfs−σf・ｃｏｓθ＝０ ・・・（５−１）
この（５−１）式において、σsは固体気体間の表面張力、σfsは固体液体間の界面張力
、σfは液体表面張力である。なお、ここでは、以下、表面張力σs＝０として計算する。
【００５３】
ここで、流体から単位面積の２つの流体表面を形成する場合、及び単位面積の流体固体表面から流体表面と固体表面を形成する場合を考える。それぞれの場合で、上述したポテンシャルエネルギｐe〔＝（Ｌ−Ｌ０）（Ｌ−Ｒ3）〕を用いて、外部より与えた仕事を計算する。具体的には、エネルギ保存則より、表面ポテンシャルエネルギの増分が与えた仕事と等しくなるので、
２σf＝〔Σｐe（Ｌi'j'）／Ｌ０2〕・Ｃf ・・・（５−２）
σf＋σs−σfs＝〔Σｐe（Ｌi'j'）／Ｌ０2〕・Ｃfs ・・・（５−３）
が成立する。従って、（５−１）式、（５−２）式、及び（５−３）式より、
Ｃfs−（Ｃf／２）・（１＋ｃｏｓθ）＝０ ・・・（５−４）
が導かれる。ここで、ポテンシャル力Ｃfs，Ｃfが、それぞれ作用力としての外力ｆlw，
ｆllに相当させることにより上記（５）式を導き出すことができる。
【００５４】
一方、粒子移動計算手段２１３は、粒子の移動を陽的に計算する手段であり、各粒子Ａ１，Ａ２についての仮速度ｕ*及び仮位置ｒ*を算出する。
このため、粒子移動計算手段２１３は、仮速度ｕ*を算出する次式を記憶している。
ｕ*＝ｕ（k）＋Δｔ[ν・∇２・ｕ＋ｆ＋ｇ]（k） ・・・（６）
ここで、ｆは上記ステップＳ１−３で算出した外力であり、ｇは重力加速度である。また、νは粒子Ａ１で表される流体の動粘性係数である。この（６）式の右辺は、タイムステップｋにおける値しか用いられていないため、その値を代入すれば、仮速度ｕ*を算出
できる。
【００５５】
更に、粒子移動計算手段２１３は、仮位置ｒ*を算出する次式を記憶している。
ｒ*＝ｒ（k）＋Δｔｕ* ・・・（７）
この（７）式の右辺も、タイムステップｋにおける値のみ用いられているため、その値を代入すれば、仮位置ｒ*を算出できる。
【００５６】
一方、粒子配置修正手段２１４は、陰的な計算を行なって、仮速度ｕ*及び仮位置ｒ*に対する粒子配置の修正を行なう。このため、粒子配置修正手段２１４は、次式の圧力のポアソン方程式を記憶している。
∇2・Ｐ（k+1）＝−[ρ／Δt2]・[(n*−n（0）)／n（0）] ・・・（８）
この（８）式は、流体の密度が粒子数密度に比例するために粒子数密度ｎ（0）が粒子
数密度ｎ*と修正項ｎ'との和で示せることと、流体の質量保存則の式と、陰的なナビア・ストークス方程式の圧力項の式とから算出されている。なお、この（８）式の左辺を、流体の粒子Ａ１の各位置におけるラプラシアンモデルで離散化すると、Ｐ（k+1）に対する
連立１次方程式が得られる。そして、粒子配置修正手段２１４は、得られた連立１次方程式を解くことにより、次のタイムステップ（ｋ＋１）の各粒子の圧力ｐ（k+1）を算出で
きる。この圧力を解く場合に、粒子配置修正手段２１４は、自由表面上の粒子に対してはＰ（k+1）＝０と設定するため、パラメータβと初期の粒子数密度n（0）との積を算出す
る。
【００５７】
また、粒子配置修正手段２１４は、速度の修正量ｕ'を算出する次式を記憶している。
ｕ'＝−（Δｔ／ρ0）∇ｐ（k+1） ・・・（９）
この（９）式には、右辺に勾配演算子があるため、公知のＭＰＳ法の勾配モデルを用いて計算を行なう。また、この場合、粒子配置修正手段２１４は、内側の壁粒子では圧力を計算し、この圧力を流体の粒子の圧力勾配の計算に用いる。これにより、壁に近づいた流体の粒子は壁から跳ね返されるようになる。
【００５８】
更に、粒子配置修正手段２１４は、流体の粒子Ａ１のそれぞれについて、修正した位置ｒ（k+1）を算出する次式を記憶している。
ｒ（k+1）＝ｒ*＋Δｔｕ' ・・・（１０）
【００５９】
次に、上述のように構成された流体解析装置において、流体の挙動についての解析を行なう場合の処理手順について、図２及び図４に基づいて説明する。
まず、制御手段２１は、計算条件の設定処理を実行する（ステップＳ１−１）。ここでは、制御手段２１は、入力手段１０を介して各種データを取得し、このデータを計算条件データ記憶部２２に記憶する。具体的には、計算制御手段２１１は、流体の粒子Ａ１と、壁の粒子Ａ２と、壁の粒子Ａ２のうち内側の壁を構成する粒子とを特定するための識別子が記憶される。更に、計算制御手段２１１は、各粒子Ａ１，Ａ２の密度ρ、重力加速度ｇ、動粘性係数ν、接触角θ及び表面張力係数σに関するデータを取得して、計算条件データ記憶部２２に記憶する。更に、計算制御手段２１１は、反復数上限値、クーラン数の上限値Ｃmax、圧力計算条件、粒子間距離Ｌ0、第１半径Ｒ1、第２半径、第３半径Ｒ3及び自由表面条件のパラメータβに関するデータを取得し、計算条件データ記憶部２２に記憶する。
【００６０】
次に、制御手段２１は、粒子の初期配置処理を実行する（ステップＳ１−２）。具体的には、計算制御手段２１１は、初期条件の流体及び壁に関するデータを、入力手段１０を介して取得し、このデータや粒子間距離Ｌ0のデータに基づいて、各粒子Ａ１，Ａ２の初
期の位置ｒ（0）、速度ｕ（0）、圧力Ｐ（0）に関するデータを設定して、粒子データ記
憶部２３に記憶する。この場合、計算制御手段２１１は、圧力を計算する内側の壁粒子における重み関数の広がりの範囲まで、壁の粒子Ａ２を配置するように、圧力を計算しない外側の壁の厚み（層）を決定する。具体的には、第１半径Ｒ1を粒子間距離Ｌ0の2.1倍と
しているので、本実施形態では、圧力を計算しない外側の壁粒子には２層が必要となり、これに対応するように壁の粒子Ａ２を配置する。
【００６１】
更に、計算制御手段２１１は、計算条件データ記憶部２２に記憶された粒子間距離Ｌ0
及び第１半径Ｒ1に基づいて、粒子数密度ｎ（0）を計算して保持する。また、計算制御手段２１１は、計算条件データ記憶部２２に記憶された粒子間距離Ｌ0、クーラン数の上限
値Ｃmax及び粒子速度から時間刻み幅Δｔを算出して保持する。
【００６２】
そして、制御手段２１は、次のステップＳ１−３以降の処理を、タイムステップ「０」から、順次、反復数上限値のタイムステップまで繰り返して行なう。ここでは、既にタイムステップｋにおいて算出された各粒子の位置ｒ（k）、速度ｕ（k）、圧力Ｐ（k）を用
いて、次のタイムステップ（ｋ＋１）における位置ｒ（k+1）、速度ｕ（k+1）、圧力Ｐ（k+1）を計算する。更に、各タイムステップにおける処理では、外力計算と、算出した外
力を用いて各粒子の仮速度及び仮位置とを陽的に計算した後、粒子毎の圧力を陰的に計算し、この圧力を用いて粒子の速度及び位置の修正量を陰的に計算する。
【００６３】
以下、所定のタイムステップｋにおける処理を説明する。ここでは、まず、制御手段２１は、各流体の粒子に作用する外力計算処理を実行する（ステップＳ１−３）。この外力計算処理では、各粒子Ａ１のそれぞれについて外力計算を行なう。この外力計算を行なう１の粒子Ａ１を着目粒子ａとし、この着目粒子ａに作用する外力の算出について、以下に説明する。
【００６４】
まず、制御手段２１の計算制御手段２１１が、着目粒子ａに対して、外力を作用させている粒子を特定する。具体的には、計算制御手段２１１は、着目粒子ａの位置ｒ（k）デ
ータと、計算条件データ記憶部２２に記憶されている第３半径データとから、着目粒子ａに対して外力計算を行なう粒子を特定するための境界領域を算出する。そして、計算制御手段２１１は、算出した境界領域と、他の粒子の位置ｒ（k）とを比較する。この比較結
果から、着目粒子ａの境界領域内にある粒子を、着目粒子ａに対して外力を作用させている粒子として特定する。例えば、図４は、着目粒子ａに対する、粒子Ａ１，Ａ２の影響を示した例である。この場合、計算制御手段２１１は、境界領域が点線で示される円により特定する。そして、この境界領域に含まれる粒子を、着目粒子ａに対して外力を作用させている粒子と特定する。
【００６５】
そして、計算制御手段２１１は、外力を作用させている粒子を特定するデータ（例えば各粒子の位置データ)を外力計算手段２１２に供給する。
外力計算手段２１２は、特定した粒子Ａ１，Ａ２の位置ｒ（k）と着目粒子ａの位置ｒ
（k）とから、着目粒子ａに対する各粒子Ａ１，Ａ２の相互距離Ｌを算出する。
【００６６】
そして、外力計算手段２１２は、第３半径Ｒ3内にあると特定した粒子Ａ１，Ａ２が着
目粒子ａに作用する外力ｆを算出する。
この場合、外力計算手段２１２は、流体の粒子Ａ１が、同じ種類の着目粒子ａに対して作用する外力ｆllを、上記（３）式を用いて算出する。具体的には、外力計算手段２１２は、計算条件データ記憶部２２から表面張力係数σを取得し、上記（４）式を用いて、ポテンシャルパラメータαを算出し、計算条件データ記憶部２２から粒子間距離Ｌ０を取得して第３半径Ｒ3を算出する。更に、外力計算手段２１２は、粒子間距離Ｌ０と、先に算
出した第３半径Ｒ3、相互距離Ｌ及びポテンシャルパラメータαとを（３）式に代入する
ことにより、各粒子Ａ１から着目粒子ａに作用する外力ｆllを算出する。
【００６７】
更に、外力計算手段２１２は、壁の粒子Ａ２が、異なる種類の着目粒子ａに対して作用する外力ｆlwを、上記（３）式及び（５）式を用いて算出する。具体的には、外力計算手段２１２は、計算条件データ記憶部２２から、接触角θ及びポテンシャルパラメータαを取得する。そして、外力計算手段２１２は、ポテンシャルパラメータαと、先に算出した相互距離Ｌとを（３）式に代入して外力ｆllを算出した後、この外力ｆllと接触角θとを（５）式に代入することにより、各粒子Ａ２から着目粒子ａに作用する外力ｆlwを算出する。
【００６８】
更に、外力計算手段２１２は、算出した各粒子Ａ１によって作用される外力ｆllと、各粒子Ａ２によって作用される外力ｆlwとを合計して、着目粒子ａの外力ｆを算出する。このとき、外力ｆll，ｆlwがベクトル量であるため、外力計算手段２１２は、ベクトルの計算により外力ｆを算出する。
【００６９】
以上の計算を、各流体の粒子Ａ１を順次、着目粒子ａとして、すべての粒子Ａ１について行ない、流体の粒子Ａ１に作用する各外力ｆを算出する。
次に、制御手段２１は、算出した外力ｆを用いて、各粒子Ａ１，Ａ２の仮速度ｕ*の算
出処理を実行する（ステップＳ１−４）。このとき、制御手段２１の粒子移動計算手段２１３は、上記（６）式を用いて、仮速度ｕ*を算出する。具体的には、粒子移動計算手段
２１３は、重力加速度ｇ及び動粘性係数νを計算条件データ記憶部２２から取得し、保持している時間刻みΔｔと、タイムステップｋにおける速度ｕ（k）のデータを粒子データ
記憶部２３から取得する。そして、粒子移動計算手段２１３は、取得したこれらのデータと、上記ステップＳ１−３において算出した外力ｆとを、上記（６）式に代入して、仮速度ｕ*を算出する。なお、この（６）式の算出においても、粒子移動計算手段２１３は、
ベクトルの計算を行なう。
【００７０】
次に、制御手段２１は、各粒子Ａ１，Ａ２の仮位置ｒ*の算出処理を実行する（ステッ
プＳ１−５）。具体的には、粒子移動計算手段２１３は、タイムステップｋにおける位置ｒ（k）を粒子データ記憶部２３から取得し、これと時間刻み幅Δｔと、ステップＳ１−
５で算出したｕ*とを、（７）式に代入して、各粒子Ａ１，Ａ２の仮位置ｒ*を算出する。
【００７１】
次に、制御手段２１は、陽的に計算した各粒子の仮速度ｕ*と仮位置ｒ*を用いて、圧力のポアソン方程式の計算処理を実行する（ステップＳ１−６）。具体的には、粒子配置修正手段２１４は、ステップ粒子移動計算手段２１３がステップＳ１−５において算出した仮位置ｒ*から、ステップＳ１−５が終了した時点での粒子数密度ｎ*を算出する。更に、粒子配置修正手段２１４は、算出した粒子数密度ｎ*と、保持している粒子数密度ｎ0とを用いて、上記（８）式の左辺を、流体の粒子Ａ１の各位置におけるラプラシアンモデルで離散化し、Ｐ（k+1）に対する連立１次方程式を算出する。
【００７２】
この離散化の際に、粒子配置修正手段２１４は、ラプラシアンモデルを適用して圧力変数ベクトルに対する係数行列を作成する。この場合、粒子配置修正手段２１４は、計算条件データ記憶部２２に記憶された第２半径データを取得し、このデータに応じたモデルの方程式を用いる。また、粒子配置修正手段２１４は、壁近傍の流体の粒子Ａ１や内側の壁の粒子Ａ２では、圧力を計算しない外側の壁の粒子Ａ２との係数をすべて「０」に設定する。ただし、係数行列の対角成分も修正し、非対角成分の和の絶対値と、対角成分の絶対値が等しくなるようにし、係数行列の対称性を維持する調整を実行する。なお、このようにすることで、粒子配置修正手段２１４は、壁面で圧力勾配「０」のノイマン境界条件を満たした計算を行なえる。
【００７３】
更に、粒子配置修正手段２１４は、計算条件データ記憶部２２に記憶されたパラメータβを取得し、このパラメータβと初期の粒子数密度n（0）との積と、算出した粒子数密度ｎ*とを比較する。そして、粒子配置修正手段２１４は、粒子数密度ｎ*が、β・n（0）より小さい粒子については、圧力のポアソン方程式においてＰ（k+1）＝０と設定する。
【００７４】
そして、上述の連立１次方程式を解くことにより、次のタイムステップ（ｋ＋１）の各粒子の圧力ｐ（k+1）が算出される。なお、この場合、粒子配置修正手段２１４は、計算
条件データ記憶部２２に記録されている圧力計算条件データが用いられる。
【００７５】
次に、制御手段２１は、圧力勾配の計算を行なうことにより、各粒子の速度と位置の修正処理を実行する（ステップＳ１−７）。具体的には、粒子配置修正手段２１４が、ステップＳ１−６で算出したタイムステップ（ｋ＋１）における各粒子の圧力ｐ（k+1）を、
上記（９）式に代入して、速度の修正量ｕ'を算出する。この場合、粒子配置修正手段２
１４は、内側の壁粒子では圧力を計算し、この圧力を流体の粒子の圧力勾配の計算に用いる。ただし、粒子配置修正手段２１４は、壁の粒子Ａ２自体の圧力勾配項は計算しない。
【００７６】
そして、粒子配置修正手段２１４は、算出した速度の修正量ｕ'を仮速度ｕ*に加算する
ことにより、修正した各粒子の速度ｕ（k+1）を算出する。
更に、粒子配置修正手段２１４は、算出した速度の修正量ｕ'と上記（１０）式とを用
いて、流体の粒子Ａ１のそれぞれについて、修正した位置ｒ（k+1）を算出する。壁の粒
子Ａ２の位置は常に固定されているため、粒子配置修正手段２１４は、壁の粒子Ａ２の位置ｒ（k+1）は計算しない。
【００７７】
次に、制御手段２１は、位置ｒ（k+1）、速度ｕ（k+1）及び圧力Ｐ（k+1）の出力処理
を実行する（ステップＳ１−８）。具体的には、計算制御手段２１１は、上記ステップＳ１−３〜Ｓ１−７の処理において算出したタイムステップ（ｋ＋１）における各粒子Ａ１，Ａ２の位置ｒ（k+1）、速度ｕ（k+1）及び圧力Ｐ（k+1）を、このタイムステップ（ｋ
＋１）と関連付けて、粒子データ記憶部２３に記憶する。更に、本実施形態では、計算制御手段２１１は、このタイムステップ（ｋ＋１）における各粒子Ａ１，Ａ２の配置を、出力手段１５において出力する。具体的には、各位置ｒ（k）に表示されていた各粒子Ａ１
，Ａ２のオブジェクトを、位置ｒ（k+1）に変更して表示し直す。この場合、各粒子Ａ１
，Ａ２の速度や圧力を、それぞれ配置された粒子Ａ１，Ａ２のオブジェクト上の矢印の大きさや色分けにより表示させるようにしてもよい。以上により、タイムステップｋにおける処理が終了する。
【００７８】
このようにしてタイムステップｋにおける処理（ステップＳ１−３〜Ｓ１−８）が終了すると、制御手段２１は、終了したタイムステップに「１」を加算した次のタイムステップ（ｋ＋１）について、上記ステップＳ１−３以降の処理を繰り返し行なう。そして、予め定めたタイムステップが終了するまで、上記処理（ステップＳ１−３〜Ｓ１−８）を繰り返し実行する。具体的には、計算制御手段２１１は、計算条件データ記憶部２２に記憶されている反復上限値データを取得し、この反復上限値と、処理が終了したタイムステップｋとを比較する。そして、計算制御手段２１１は、処理が終了したタイムステップｋに「１」を加算した値が反復上限値を超えた場合には、流体解析処理を終了する。
【００７９】
本実施形態によれば、以下のような効果を得ることができる。
・ 本実施形態では、外力計算手段２１２は、第３半径Ｒ3に含まれる流体の粒子Ａ１
から受ける外力ｆllは（３）式を用いて算出する。また、外力計算手段２１２は、同じ第３半径Ｒ3に含まれる壁の粒子Ａ２から受ける外力ｆlwは（５）式を用いて算出する。こ
の（５）式は、（３）式で算出した外力ｆllに対して、（１＋ｃｏｓθ）／２を乗算した算式である。この（１＋ｃｏｓθ）／２という係数は、液相の界面エネルギと壁−液相の界面エネルギ関係式であるヤングの式から算出される係数である。このため、外力計算手段２１２は、エネルギ関係式から壁の粒子Ａ２の作用力を算出している。そして、粒子移動計算手段２１３は、これら外力ｆlw，ｆllの合計値を用いて、各粒子の仮速度ｕ*を算
出し（ステップＳ１−４）、この仮速度ｕ*から仮位置ｒ*を算出する（ステップＳ１−５）。更に、粒子配置修正手段２１４は、仮速度ｕ*、仮位置ｒ*及びポアソン方程式の計算を行ない（ステップＳ１−６）、圧力勾配の計算を行なって各粒子の速度と位置を修正し（ステップＳ１−７）、次のタイムステップの各粒子の位置、速度及び圧力を出力する。そして、このようなタイムステップの処理を、反復数上限値まで繰り返す。
【００８０】
従って、本実施形態の制御手段２１は、表面張力モデルにおいて、エネルギ関係式を用いて壁の粒子Ａ２からの作用力を算出している。このため、表面張力モデルの１つである濡れ性モデルであってもエネルギ関係式が成り立つため、濡れ性モデルも表面張力モデルとして統一的に取り扱うことができる。この場合、外力計算手段２１２は、粒子数密度の減少分から算出される界面の曲率を用いずに第２の粒子の作用力を算出するので、空間解像度が粗くても、流体の挙動の解析を的確に行なうことができる。
【００８１】
更に、エネルギの関係式は、流体や壁がどのような状態にあっても成立する。外力計算
手段２１２は、このエネルギの関係式を用いて第２の粒子の作用力を算出しているため、界面が動く場合など種々の条件下で用いることができ、汎用性が高い。
【００８２】
・ 本実施形態では、上記（３）式として、他の粒子との距離に応じたポテンシャルエネルギの算出式を用いている。このため、エネルギの関係式である（５）式と整合がよく、流体の挙動の解析を、更に的確に行なうことができる。
【００８３】
また、上記実施形態は、以下のように変更してもよい。
○ 上記実施形態において、外力計算手段２１２は、外力ｆllを計算するポテンシャル関数として、上記（３）式を用いた。これに限らず、外力ｆllを算出する式は、他の関数を用いてもよい。この場合、（５）式との関係上、エネルギ関係式で外力ｆllを表現した関数を用いるのがよい。
【００８４】
○ 上記実施形態において、計算条件データ記憶部２２に記録された反復数上限値を用いて、制御手段２１は計算処理を終了させたが、計算処理の終了方法は、これに限定されるものではない。例えば、終了時間を記録して、これに対応させた反復数上限値により終了させてもよい。また、制御手段２１は、入力手段１０を介して終了信号が入力された段階で流体解析処理を終了するにしてもよい。
【００８５】
○ 上記実施形態においては、制御手段２１は、圧縮率がきわめて小さく圧縮性を考慮しない流体の挙動について解析を行なった。これに限らず、解析する対象によっては、流体の圧縮性を考慮してもよい。この場合、基準となる初期の密度及び圧力から、圧縮性を考慮した公知の修正を行なう。これにより、圧縮性が小さい流体において速い現象をも取り扱うことができる。
【図面の簡単な説明】
【００８６】
【図１】実施形態におけるシステムの概略図。
【図２】実施形態における処理手順を説明するための流れ図。
【図３】接触角θを説明するための説明図。
【図４】各粒子から受ける外力の算出を説明するための説明図。
【図５】（４)式の導出過程を説明するための説明図。
【符号の説明】
【００８７】
θ…接触角、Ａ１…第１の粒子、Ａ２…第２の粒子、Ｐ，ｐ…圧力、Ｒ3…所定範囲と
しての第３半径、ｆll，ｆlw…作用力としての外力、ｒ…位置、ｒ*…仮位置、ｕ…速度
、ｕ*…仮速度、２０…流体解析手段、２１…制御手段、２２…計算条件データ記憶手段
としての計算条件データ記憶部、２３…粒子データ記憶手段としての粒子データ記憶部、２１１…特定手段としての計算制御手段、２１２…外力計算手段、２１３…粒子移動計算手段、２１４…粒子配置修正手段。
【技術分野】
【０００１】
本発明は、粒子法を用いて連続体である流体の挙動を解析する流体解析装置、流体解析方法及び流体解析プログラムに関する。
【背景技術】
【０００２】
近年、流体の挙動を解析する方法の１つとして粒子法が知られている（例えば、非特許文献１及び非特許文献２参照。）。粒子法は、流体を複数の粒子の集まりとして表わし、これら粒子間の相互作用の計算によって、流体の挙動を解析する。この粒子法は、有限体積法や有限要素法の計算において、手間のかかる格子生成が不要という利点がある。
【０００３】
更に、この粒子法の１つとしてＭＰＳ（moving particle semi-implicit method）法がある。このＭＰＳ法は、連続体の離散的な計算を、粒子間相互作用モデルを用いて行なう。このＭＰＳ法においては、粒子に対して所定の基準半径以内にある近接の粒子が、その粒子に対して相互作用を与えると仮定する。なお、基準半径が小さい程、粒子の数が減少するので、計算時間が短くなるが、小さすぎると計算が不安定になる。そこで、計算時間と計算の安定性の観点から、一般に、基準半径は初期状態の粒子間距離の２〜４倍に設定される。この場合、２次元の座標系においては、基準半径以内にある粒子の数は１２〜４４個程度となる。
【０００４】
このＭＰＳ法を用いて非圧縮流体の挙動を解析する場合には、粒子数密度が一定であることが利用される。ここで、粒子数密度をｎiとし、粒子ｉ及びその近傍の粒子ｊの位置
ベクトルをそれぞれｒi，ｒjとすると、粒子数密度ｎiは、粒子ｉの位置における重み関
数ｗの和（次式）で示される。
ｎi＝Σｗ（│ｒj−ｒi│） ・・・（１）
ここで、各粒子の保持する質量が一定であると仮定すると、粒子数密度は流体の密度と比例する。従って、非圧縮流体の場合には粒子数密度が一定であり、ＭＰＳ法では、この条件を満足するように解析を行なう。
【０００５】
ところで、マイクロ流路における流体の挙動など、長さスケールが小さい流れにおいては、表面張力の影響が大きくなる。そこで、表面張力を考慮して解析を行なうことがある（非特許文献２のｐ．３５−３８参照）。非特許文献２に記載の技術では、表面張力係数σ、曲率κ、表面張力が表面のみに働くためのデルタ関数δ、表面と垂直な方向の単位ベクトルｖを用いて、運動量保存則の方程式に表面張力項を導入して解析を行なう。従って、表面張力を考慮する場合には、この表面張力項が示す部分を、粒子に働く力としてモデル化する必要がある。
【０００６】
従来、このモデル化においては、表面に位置する粒子を抽出し、この粒子に対する曲率計算を行なう。この曲率計算においては、上記（１)式に用いた重み関数とは異なる重み
関数（以下、曲率計算用重み関数という。）を用いた粒子数密度を算出する。この曲率計算用重み関数は、所定領域内では常に値が「１」で、領域外では常に値が「０」に設定された関数である。
【０００７】
そして、表面が平らな場合（すなわち曲率κが「０」の場合）の曲率計算用粒子数密度を、初期粒子配置において予め計算しておく。この曲率計算用粒子数密度は所定領域内にある単なる近傍粒子の数であるから、表面が凸の場合には、初期粒子配置における曲率計算用粒子数密度より小さくなり、表面が凹の場合には大きくなる。従って、曲率計算用粒子数密度の変化分から界面の曲率κを算出し、これを用いて表面張力を算出する。
【０００８】
また、表面張力の他の計算方法として、粒子間にポテンシャル力を働かせて計算する技術もある（非特許文献３を参照。）。この非特許文献３に記載の方法を用いることにより、比較的容易に表面張力の計算を行なうことができる。
【非特許文献１】越塚誠一著，「インテリジェント・エンジニアリング・シリーズ 数値流体力学」，初版，株式会社培風館，１９９７年，ｐ．１５１−１９４
【非特許文献２】越塚誠一著，「計算力学レクチャーシリーズ５ 粒子法」，丸善株式会社，２００５年，ｐ．１−７，ｐ．２１−３８
【非特許文献３】Ｎ．Shirakawa、他３名，「Journal of Nuclear Science and Technology」，２００１年，第３８巻，ｐ．３９２−４０２
【発明の開示】
【発明が解決しようとする課題】
【０００９】
上述した界面の曲率κを算出して表面張力を計算する場合、曲率計算用粒子数密度の計算精度を保つために、十分な粒子数（例えば直径に１０粒以上）が必要である。このため、空間解像度が粗い場合（すなわち計算に用いる粒子数が少ないように粗く設定する場合）には、曲率計算用粒子数密度の計算精度が悪くなり、流体の挙動を的確に解析することができなかった。
【００１０】
ところで、表面張力が影響するモデルとして、壁に対して液体の水滴が接触角θで静止した現象の「濡れ性」のモデルがある。この場合に用いる表面張力項の計算においても、上述したように、表面と垂直な方向の単位ベクトルｖを用いる。このため、「濡れ性」の挙動解析を行なう場合には、接触角θから、表面と垂直な方向のベクトル成分を幾何学的に計算し、その結果を用いていた。この場合、曲率計算用粒子数密度は、上述した表面張力の通常の計算よりも、誤差が大きくなっていた。更に、この「濡れ性」では接触角θをパラメータとして用いるため、解析対象の流体（液体）が流動する場合には、流動状態での接触角θを把握してパラメータとして設定する必要がある。しかし、このような接触角θを正確に把握することは困難である。そのため、ダイナミックな解析を行なうことは難しく、汎用性に乏しかった。
【００１１】
また、非特許文献３に記載の技術の場合、ポテンシャル力と表面張力係数σとの関係が不明であるため、パラメータ調整が必要となり、濡れ性の取り扱いについても考慮されていない。
【００１２】
本発明は、上述した問題に鑑みてなされ、その目的は、汎用性が高く、濡れ性モデルを統一的に取り扱うことができ、空間解像度が粗くても精度よく解析を行なうことができる流体解析装置、流体解析方法及び流体解析プログラムを提供することにある。
【課題を解決するための手段】
【００１３】
上記問題点を解決するために、請求項１に記載の発明は、流体を表した第１の粒子毎に、速度、位置、圧力に関するデータと、前記流体に接触する壁を表した第２の粒子毎に、速度、位置、圧力に関するデータとを記録する粒子データ記憶手段と、前記流体と前記壁との接触角θに関するデータと、各粒子について他の粒子からの作用力を受ける作用範囲に関するデータとを記録した計算条件データ記憶手段と、各粒子の挙動を計算する制御手段とから構成され、流体の挙動を解析する粒子法を用いた解析装置であって、前記制御手段が、所定のタイムステップにおいて、各第１の粒子の配置において前記作用範囲に含まれる他の粒子を特定する特定手段と、特定された他の粒子のうち第１の粒子については、第１の関数により第１の作用力を算出し、特定された他の粒子のうち第２の粒子については、前記第１の関数に対して（１＋ｃｏｓθ）／２を乗算した第２の関数により第２の作用力を算出する作用力算出手段と、前記各作用力により算出した各粒子の仮速度を用いて
、仮配置を決定する粒子移動計算手段と、前記仮配置を用いて、前記粒子の数密度を一定にする計算により、次回タイムステップの圧力を算出し、この圧力を用いて算出した圧力勾配から各粒子の修正後の速度を算出し、この速度を用いて各粒子の位置を修正し、前記粒子データ記憶手段に記録する粒子配置修正手段と、タイムステップを進めて、前記特定手段、前記粒子移動計算手段及び前記粒子配置修正手段を用いた各粒子の挙動の計算処理を繰り返す計算制御手段とを備えたことを要旨とする。
【００１４】
請求項２に記載の発明は、請求項１に記載の流体解析装置において、前記第１の関数は、作用力を及ぼす前記他の粒子との距離に応じたポテンシャルエネルギにより算出する関数であることを要旨とする。
【００１５】
請求項３に記載の発明は、流体を表した第１の粒子毎に、速度、位置、圧力に関するデータと、前記流体に接触する壁を表した第２の粒子毎に、速度、位置、圧力に関するデータとを記録する粒子データ記憶手段と、前記流体と前記壁との接触角θに関するデータと、各粒子について他の粒子からの作用力を受ける作用範囲に関するデータとを記録した計算条件データ記憶手段と、各粒子の挙動を計算する制御手段とから構成された解析装置を用いて、粒子法による流体の挙動を解析する方法であって、前記制御手段が、所定のタイムステップにおいて、各第１の粒子の配置において前記作用範囲に含まれる他の粒子を特定する特定段階と、特定された他の粒子のうち第１の粒子については、第１の関数により第１の作用力を算出し、特定された他の粒子のうち第２の粒子については、前記第１の関数に対して（１＋ｃｏｓθ）／２を乗算した第２の関数により第２の作用力を算出する作用力算出段階と、前記各作用力により算出した各粒子の仮速度を用いて、仮配置を決定する粒子移動計算段階と、前記仮配置を用いて、前記粒子の数密度を一定にする計算により、次回タイムステップの圧力を算出し、この圧力を用いて算出した圧力勾配から各粒子の修正後の速度を算出し、この速度を用いて各粒子の位置を修正し、前記粒子データ記憶手段に記録する粒子配置修正段階と、タイムステップを進めて、前記特定段階、前記粒子移動計算段階及び前記粒子配置修正段階による各粒子の挙動の計算処理を繰り返す計算制御段階とを実行することを要旨とする。
【００１６】
請求項４に記載の発明は、請求項３に記載の流体解析方法において、前記第１の関数は、作用力を及ぼす前記他の粒子との距離に応じたポテンシャルエネルギにより算出する関数であることを要旨とする。
【００１７】
請求項５に記載の発明は、流体を表した第１の粒子毎に、速度、位置、圧力に関するデータと、前記流体に接触する壁を表した第２の粒子毎に、速度、位置、圧力に関するデータとを記録する粒子データ記憶手段と、前記流体と前記壁との接触角θに関するデータと、各粒子について他の粒子からの作用力を受ける作用範囲に関するデータとを記録した計算条件データ記憶手段と、各粒子の挙動を計算する制御手段とから構成された解析装置を用いて、粒子法による流体の挙動を解析するプログラムであって、前記制御手段を、所定のタイムステップにおいて、各第１の粒子の配置において前記作用範囲に含まれる他の粒子を特定する特定手段、特定された他の粒子のうち第１の粒子については、第１の関数により第１の作用力を算出し、特定された他の粒子のうち第２の粒子については、前記第１の関数に対して（１＋ｃｏｓθ）／２を乗算した第２の関数により第２の作用力を算出する作用力算出手段、前記各作用力により算出した各粒子の仮速度を用いて、仮配置を決定する粒子移動計算手段、及び前記仮配置を用いて、前記粒子の数密度を一定にする計算により、次回タイムステップの圧力を算出し、この圧力を用いて算出した圧力勾配から各粒子の修正後の速度を算出し、この速度を用いて各粒子の位置を修正し、前記粒子データ記憶手段に記録する粒子配置修正手段と、タイムステップを進めて、前記特定手段、前記粒子移動計算手段及び前記粒子配置修正手段を用いた各粒子の挙動の計算処理を繰り返す計算制御手段として機能させることを要旨とする。
【００１８】
請求項６に記載の発明は、請求項５に記載の流体解析プログラムにおいて、前記第１の関数は、作用力を及ぼす前記他の粒子との距離に応じたポテンシャルエネルギにより算出する関数であることを要旨とする。
【００１９】
（作用）
本発明によれば、制御手段が、作用範囲に含まれる第１の粒子については、第１の関数により第１の作用力を算出し、作用範囲に含まれる第２の粒子については、第１の関数に対して（１＋ｃｏｓθ）／２を乗算した第２の関数により第２の作用力を算出する。制御手段は、各作用力により算出した各粒子の速度を用いて、仮配置を決定する。更に、制御手段は、仮配置を用いて、粒子の数密度を一定にする計算により、次回タイムステップの圧力を算出し、この圧力を用いて算出した圧力勾配から各粒子の修正速度を算出し、これから各粒子の位置を補正する。以上の処理を、制御手段は、タイムステップを進めて、繰り返して行なう。（１＋ｃｏｓθ）／２は、液相の界面エネルギと壁−液相の界面エネルギの関係式であるヤングの式から算出される係数である。このため、制御手段は、エネルギ関係式から第２の粒子の作用力を算出し、これら第１及び第２の粒子の作用力を用いて流体の挙動の解析を行なうので、濡れ性モデルを、表面張力モデルとして統一的に扱うことができる。
【００２０】
また、制御手段は、粒子数密度の減少分から算出される界面の曲率を用いずに、第２の粒子の作用力を算出しているので、空間解像度が粗くても、流体の挙動の解析を精度よく行なうことができる。
【００２１】
更に、エネルギ関係式は、流体や壁がどのような状態にあっても成立する。従って、制御手段は、第２の粒子の作用力をエネルギ関係式から算出しているので、例えば界面が動く場合など、種々の条件下でも用いることができ、汎用性が高い。
【００２２】
本発明によれば、第１の関数は、作用力を及ぼす前記他の粒子との距離に応じたポテンシャルエネルギにより算出する関数である。このため、ヤングの式から算出される係数と整合がよく、流体の挙動の解析を、より的確に良好に行なうことができる。
【発明の効果】
【００２３】
本発明によれば、粒子法において、汎用性が高く、濡れ性モデルを統一的に取り扱うことができ、空間解像度が粗くても精度よく解析を行なうことができる。
【発明を実施するための最良の形態】
【００２４】
以下、本発明を具体化した流体解析装置の一実施形態を図１〜図５に基づいて説明する。本実施形態では、壁が剛体で構成された収容器内の流体の挙動について解析する場合を想定する。ここでは、流体として、水のように圧縮率が小さい流体を用いる場合を想定する。また、本実施形態では、説明を簡単にするために、流体の流入及び流出がなく、２次元の解析を行なう場合を想定する。
【００２５】
図１に示すように、本実施形態の流体解析装置は、入力手段１０、出力手段１５及び流体解析手段２０を備えている。
入力手段１０は、流体解析手段２０に各種データを入力するための手段であり、例えばキーボードやポインティング・デバイスなどを用いる。そして、入力手段１０は、入力された計算条件データや初期条件における粒子データなどを流体解析手段２０に供給する。
【００２６】
出力手段１５は、流体解析手段２０による解析結果を出力するために用いられ、具体的には、ディスプレイやデータ記憶部などを用いる。
流体解析手段２０は、データの処理を行なうためのコンピュータシステムであって、制御手段２１、計算条件データ記憶部２２及び粒子データ記憶部２３を備える。
【００２７】
制御手段２１は、図示しないＣＰＵ、ＲＡＭ及びＲＯＭ等を有し、後述する処理（特定段階、作用力算出段階、粒子移動計算段階、粒子配置修正段階及び計算制御段階等を含む処理）を行なう。そして、このための流体解析プログラムを実行することにより、制御手段２１は、特定手段としても機能する計算制御手段２１１、作用力算出手段としての外力計算手段２１２、粒子移動計算手段２１３及び粒子配置修正手段２１４として機能する。
【００２８】
計算条件データ記憶部２２には、流体解析に用いる物性や計算制御などの計算条件データが記憶されている。この計算条件データは、流体解析を行なう前に、入力手段１０を介して登録される。本実施形態においては、解析対象の粒子としては、挙動を解析する流体を構成する第１の粒子Ａ１と、壁を構成する第２の粒子Ａ２がある。更に、壁の粒子Ａ２においては、第１の粒子Ａ１と接する壁粒子と、その外側の壁粒子とがある。このため、この計算条件データには、粒子識別子、内側壁粒子の種類識別子、壁粒子の種類識別子に関するデータが含まれる。更に、計算条件データには、密度、重力加速度、動粘性係数、接触角及び表面張力係数に関するデータが含まれる。
【００２９】
粒子識別子データ領域には、この流体解析の対象となる粒子の種類を特定するための識別子に関するデータが記録されている。このデータによって、制御手段２１は、計算に用いるデータのうち、どのデータが粒子の種類を特定するためのデータであるかを特定できる。
【００３０】
内側壁粒子の種類識別子データ領域には、第１の粒子Ａ１と接して圧力計算を行なう内側の壁の粒子の種類を特定するための識別子に関するデータが記録されている。これにより、内側の壁の粒子の種類識別子データと、内側壁粒子であるというフラグとが関連付けられることになる。この種類識別子によって特定される粒子Ａ２について、制御手段２１は、その位置は常に一定とし、この粒子Ａ２についての圧力を算出する。
【００３１】
壁粒子の種類識別子データ領域には、内側の壁の粒子Ａ２の外側に存在して壁を構成する粒子を特定するための識別子に関するデータが記録されている。これにより、制御手段２１は、この種類識別子によって特定される壁の粒子Ａ２について、その位置は常に一定とし、かつこの壁の粒子Ａ２の圧力は算出しない。
【００３２】
密度データ領域には、各粒子識別子に対応させて、各粒子によって構成される各物質の密度ρに関するデータが記録されている。
重力加速度データ領域及び動粘性係数データ領域には、それぞれ、計算に用いる重力加速度ｇ及び流体の動粘性係数νに関するデータが記録されている。
【００３３】
接触角データ領域には、接触角θに関するデータが記録されている。この接触角θは、図３に示すように、粒子Ａ２で表される壁（固体）の表面に、粒子Ａ１で表される流体（液体）が置かれた場合に、静的な状態での接触点において液体がなす角度である。
【００３４】
表面張力係数データ領域には、表面張力係数σに関するデータが記録されている。この表面張力係数σは、後述するポテンシャル関数のポテンシャルパラメータαを算出するために用いられる。
【００３５】
更に、この計算条件データには、反復数上限値、クーラン数上限値、圧力計算条件、粒子間距離、第１半径、第２半径、第３半径及び自由表面条件パラメータに関するデータが含まれる。
【００３６】
反復数上限値データ領域には、流体の挙動を計算する場合のタイムステップの上限値に関するデータが記録されている。制御手段２１は、この上限値に達した場合に計算処理を終了する。
【００３７】
クーラン数上限値データ領域には、クーラン数の上限値Ｃmaxに関するデータが記録さ
れている。クーラン数は、時間刻み幅Δｔを算出するために用いられる。
圧力計算条件データ領域には、圧力計算における収束条件や反復数の上限及び下限に関するデータが記録されている。
【００３８】
粒子間距離データ領域には、初期配置における各粒子と、この近傍の粒子との粒子間距離Ｌ0に関するデータが記録されている。
第１半径データ領域には、粒子数密度を計算するために用いる第１半径Ｒ1に関するデ
ータが記録されている。この第１半径Ｒ1は、粒子数密度n（0）を算出するために用いら
れ、具体的には、経験値から、粒子間距離Ｌ0の 2.1倍の値が用いられる。
【００３９】
第２半径データ領域には、後述する圧力のポアソン方程式の計算におけるラプラシアンモデルにおいて用いられる半径に関するデータが記録されている。具体的には、経験値から、この第２半径データとして、粒子間距離Ｌ0の４倍の値が用いられる。
【００４０】
第３半径データ領域には、後述するポテンシャル関数に用いる第３半径Ｒ3に関するデ
ータが記録されている。本実施形態では、この第３半径Ｒ3データとして、粒子間距離Ｌ0の３倍の値を用いる。この第３半径Ｒ3は、図４に示すように作用力を計算するために用
いる半径として用いられ、この第３半径Ｒ3内が特許請求の範囲の「作用範囲」に相当す
る。
【００４１】
自由表面条件パラメータ領域には、自由表面にある粒子を特定するためのパラメータβが記録されている。制御手段２１は、後述する陽的な計算が終了した時点における粒子数密度ｎ*が、β・n（0）より小さい粒子を自由表面上の粒子と特定する。自由表面の外側
には粒子が配置されないため、自由表面上の粒子の粒子数密度が低下するからである。そして、制御手段２１は、特定した自由表面上の粒子についての圧力を、圧力のポアソン方程式を解く際にＰ（k+1）＝０とする。なお、これにより、ディリクレ境界条件を設定し
たことになる。
【００４２】
粒子データ記憶部２３には、各タイムステップにおける粒子に関するデータが記録される。この粒子データ記憶部２３には、タイムステップ毎に、粒子毎の種類識別子、位置、速度及び圧力に関するデータが記録される。すなわち、粒子についての種類識別子、位置、速度及び圧力に関するデータは、各タイムステップにおいて粒子総数分だけ記憶されることになる。
【００４３】
タイムステップデータ領域には、各タイムステップを特定するためのデータが記録される。
種類識別子データ領域、位置データ領域、速度データ領域及び圧力データ領域には、このタイムステップにおける、それぞれの粒子の種類を特定するための識別子、位置、速度及び圧力に関するデータが記録される。本実施形態では、位置データとして粒子のｘ座標及びｙ座標を用い、速度データとして、ｘ方向速度成分及びｙ方向速度成分を用いる。
【００４４】
次に、流体解析手段２０の制御手段２１が備える各処理手段（２１１〜２１４）について詳述する。
計算制御手段２１１は、入力手段１０を介して、計算条件や粒子の初期配置に関するデ
ータを取得して、計算条件データ記憶部２２及び粒子データ記憶部２３にデータを記憶する。
【００４５】
更に、計算制御手段２１１は、時間刻み幅Δｔを算出する次式を記憶している。
Δｔ＝Ｌ0・Ｃmax／ｕmax ・・・（２）
ここで、時間刻み幅Δｔが大きい程、計算が効率的に行なえるので、時間刻み幅Δｔはできるだけ大きくすることが好ましい。しかし、数値安定性を保つためには、クーラン数の最大値に制限を設ける必要がある。上記（２）式において、時間刻み幅Δｔや粒子間距離Ｌ0は一定であるので、クーラン数が最大になる粒子は、速度の絶対値が最大値をもつ粒
子である。従って、計算制御手段２１１は、粒子Ａ１の最大値の速度と、クーラン数の上限値Ｃmaxと、粒子間距離Ｌ0とを上記（２）式に代入して、時間刻み幅Δｔを算出する。
【００４６】
また、計算制御手段２１１は、上記（１）式を記憶しており、外力計算手段２１２が算出した各粒子Ａ１，Ａ２の仮位置ｒ*に基づく粒子数密度ｎ*を、この（１）式を用いて算出する。
【００４７】
更に、計算制御手段２１１は、タイムステップｋにおいて算出した各粒子Ａ１，Ａ２の位置ｒ（k+1）、速度ｕ（k+1）及び圧力Ｐ（k+1）を、出力手段１５を介して出力する。
そして、計算制御手段２１１は、解析処理を終了するかを判断し、終了しない場合には次のタイムステップにおいて処理を繰り返して実行する。
【００４８】
一方、外力計算手段２１２は、流体の粒子Ａ１のそれぞれに対して、第３半径Ｒ3内に
ある粒子を特定する。
更に、この外力計算手段２１２は、各粒子の相互距離Ｌを算出する式を記憶している。このため、外力計算手段２１２は、この式と、各粒子の位置データとを用いて、各粒子の相互距離Ｌを算出する。
【００４９】
また、外力計算手段２１２は、流体の粒子Ａ１が、他の流体の粒子Ａ１に作用する第１の作用力としての外力ｆllを算出する次式（第１の関数）を記憶している。
ｆll＝α・（Ｌ−Ｌ０）・（Ｌ−Ｒ3） ・・・（３）
この（３）式は、界面付近の粒子のみに対して中心力が作用するとして分子間力を模擬したポテンシャル関数である。ただし、相互距離Ｌが第３半径Ｒ3以上の場合には、ｆll
＝０とする。
【００５０】
更に、外力計算手段２１２は、上記（３）式におけるポテンシャルパラメータαを算出する次式を記憶している。
α＝２・Ｌ０2・σ／Σｐe（Ｌi'j'） ・・・（４）
ここで、この（４）式の導出について説明する。表面張力係数σとは、単位面積の表面を形成するために必要なエネルギである。そこで、ポテンシャルエネルギをｐe〔＝（Ｌ
−Ｌ０）（Ｌ−Ｒ3）〕と表わした場合、ポテンシャル関数が〔α・ｐe〕のときに単位面積の表面を形成するために必要なエネルギが表面張力係数σとなるようにポテンシャルパラメータαを算出する。この算出においても、相互距離Ｌが第３半径Ｒ3以上の場合には
「０」とする。ここで、図５に示すように、２つの領域（第１領域Ｅ１、第２領域Ｅ２）から構成された物質を想定し、この物質から第１領域Ｅ１を削除して、この第１領域Ｅ１に接していた第２領域Ｅ２に界面を形成する場合を検討する。この場合、ポテンシャルパラメータαは、界面を形成していない場合（第１領域Ｅ１が第２領域Ｅ２に接している状態）と界面を形成した場合（第１領域Ｅ１を除去した状態）とのエネルギ差と、表面張力係数σとの関係式から算出される。具体的には、１つの粒子の面積ｄｓ＝Ｌ０２の表面のもつポテンシャルエネルギは、第１領域Ｅ１を構成する粒子ｊ'を、第２領域Ｅ２におい
て面積ｄｓ内に存在する粒子ｉ'（図５における領域ｅ２の粒子）から引き離すために必
要な仕事の半分である。粒子間相対位置ベクトルＬi'j'のポテンシャルエネルギをｐe（
Ｌi'j'）とした場合、以下のように表わされる。
σ＝（１／２）・〔Σｐe（Ｌi'j'）／Ｌ０2〕・α ・・・（４−１）
この（４−１）式を用いることにより、上記（４）を算出することができる。
【００５１】
外力計算手段２１２は、壁の粒子Ａ２が、流体の粒子Ａ１に作用する第２の作用力としての外力ｆlwを算出する次式（第２の関数）を記憶している。
ｆlw＝（１＋ｃｏｓθ）ｆll／２ ・・・（５）
ここで、この（５）式の導出について、液体粒子間のポテンシャル力Ｃfと液体固体間
のポテンシャル力Ｃfsの比とヤングの式との対応関係から説明する。
【００５２】
ヤングの式は、次式で示される。
σs−σfs−σf・ｃｏｓθ＝０ ・・・（５−１）
この（５−１）式において、σsは固体気体間の表面張力、σfsは固体液体間の界面張力
、σfは液体表面張力である。なお、ここでは、以下、表面張力σs＝０として計算する。
【００５３】
ここで、流体から単位面積の２つの流体表面を形成する場合、及び単位面積の流体固体表面から流体表面と固体表面を形成する場合を考える。それぞれの場合で、上述したポテンシャルエネルギｐe〔＝（Ｌ−Ｌ０）（Ｌ−Ｒ3）〕を用いて、外部より与えた仕事を計算する。具体的には、エネルギ保存則より、表面ポテンシャルエネルギの増分が与えた仕事と等しくなるので、
２σf＝〔Σｐe（Ｌi'j'）／Ｌ０2〕・Ｃf ・・・（５−２）
σf＋σs−σfs＝〔Σｐe（Ｌi'j'）／Ｌ０2〕・Ｃfs ・・・（５−３）
が成立する。従って、（５−１）式、（５−２）式、及び（５−３）式より、
Ｃfs−（Ｃf／２）・（１＋ｃｏｓθ）＝０ ・・・（５−４）
が導かれる。ここで、ポテンシャル力Ｃfs，Ｃfが、それぞれ作用力としての外力ｆlw，
ｆllに相当させることにより上記（５）式を導き出すことができる。
【００５４】
一方、粒子移動計算手段２１３は、粒子の移動を陽的に計算する手段であり、各粒子Ａ１，Ａ２についての仮速度ｕ*及び仮位置ｒ*を算出する。
このため、粒子移動計算手段２１３は、仮速度ｕ*を算出する次式を記憶している。
ｕ*＝ｕ（k）＋Δｔ[ν・∇２・ｕ＋ｆ＋ｇ]（k） ・・・（６）
ここで、ｆは上記ステップＳ１−３で算出した外力であり、ｇは重力加速度である。また、νは粒子Ａ１で表される流体の動粘性係数である。この（６）式の右辺は、タイムステップｋにおける値しか用いられていないため、その値を代入すれば、仮速度ｕ*を算出
できる。
【００５５】
更に、粒子移動計算手段２１３は、仮位置ｒ*を算出する次式を記憶している。
ｒ*＝ｒ（k）＋Δｔｕ* ・・・（７）
この（７）式の右辺も、タイムステップｋにおける値のみ用いられているため、その値を代入すれば、仮位置ｒ*を算出できる。
【００５６】
一方、粒子配置修正手段２１４は、陰的な計算を行なって、仮速度ｕ*及び仮位置ｒ*に対する粒子配置の修正を行なう。このため、粒子配置修正手段２１４は、次式の圧力のポアソン方程式を記憶している。
∇2・Ｐ（k+1）＝−[ρ／Δt2]・[(n*−n（0）)／n（0）] ・・・（８）
この（８）式は、流体の密度が粒子数密度に比例するために粒子数密度ｎ（0）が粒子
数密度ｎ*と修正項ｎ'との和で示せることと、流体の質量保存則の式と、陰的なナビア・ストークス方程式の圧力項の式とから算出されている。なお、この（８）式の左辺を、流体の粒子Ａ１の各位置におけるラプラシアンモデルで離散化すると、Ｐ（k+1）に対する
連立１次方程式が得られる。そして、粒子配置修正手段２１４は、得られた連立１次方程式を解くことにより、次のタイムステップ（ｋ＋１）の各粒子の圧力ｐ（k+1）を算出で
きる。この圧力を解く場合に、粒子配置修正手段２１４は、自由表面上の粒子に対してはＰ（k+1）＝０と設定するため、パラメータβと初期の粒子数密度n（0）との積を算出す
る。
【００５７】
また、粒子配置修正手段２１４は、速度の修正量ｕ'を算出する次式を記憶している。
ｕ'＝−（Δｔ／ρ0）∇ｐ（k+1） ・・・（９）
この（９）式には、右辺に勾配演算子があるため、公知のＭＰＳ法の勾配モデルを用いて計算を行なう。また、この場合、粒子配置修正手段２１４は、内側の壁粒子では圧力を計算し、この圧力を流体の粒子の圧力勾配の計算に用いる。これにより、壁に近づいた流体の粒子は壁から跳ね返されるようになる。
【００５８】
更に、粒子配置修正手段２１４は、流体の粒子Ａ１のそれぞれについて、修正した位置ｒ（k+1）を算出する次式を記憶している。
ｒ（k+1）＝ｒ*＋Δｔｕ' ・・・（１０）
【００５９】
次に、上述のように構成された流体解析装置において、流体の挙動についての解析を行なう場合の処理手順について、図２及び図４に基づいて説明する。
まず、制御手段２１は、計算条件の設定処理を実行する（ステップＳ１−１）。ここでは、制御手段２１は、入力手段１０を介して各種データを取得し、このデータを計算条件データ記憶部２２に記憶する。具体的には、計算制御手段２１１は、流体の粒子Ａ１と、壁の粒子Ａ２と、壁の粒子Ａ２のうち内側の壁を構成する粒子とを特定するための識別子が記憶される。更に、計算制御手段２１１は、各粒子Ａ１，Ａ２の密度ρ、重力加速度ｇ、動粘性係数ν、接触角θ及び表面張力係数σに関するデータを取得して、計算条件データ記憶部２２に記憶する。更に、計算制御手段２１１は、反復数上限値、クーラン数の上限値Ｃmax、圧力計算条件、粒子間距離Ｌ0、第１半径Ｒ1、第２半径、第３半径Ｒ3及び自由表面条件のパラメータβに関するデータを取得し、計算条件データ記憶部２２に記憶する。
【００６０】
次に、制御手段２１は、粒子の初期配置処理を実行する（ステップＳ１−２）。具体的には、計算制御手段２１１は、初期条件の流体及び壁に関するデータを、入力手段１０を介して取得し、このデータや粒子間距離Ｌ0のデータに基づいて、各粒子Ａ１，Ａ２の初
期の位置ｒ（0）、速度ｕ（0）、圧力Ｐ（0）に関するデータを設定して、粒子データ記
憶部２３に記憶する。この場合、計算制御手段２１１は、圧力を計算する内側の壁粒子における重み関数の広がりの範囲まで、壁の粒子Ａ２を配置するように、圧力を計算しない外側の壁の厚み（層）を決定する。具体的には、第１半径Ｒ1を粒子間距離Ｌ0の2.1倍と
しているので、本実施形態では、圧力を計算しない外側の壁粒子には２層が必要となり、これに対応するように壁の粒子Ａ２を配置する。
【００６１】
更に、計算制御手段２１１は、計算条件データ記憶部２２に記憶された粒子間距離Ｌ0
及び第１半径Ｒ1に基づいて、粒子数密度ｎ（0）を計算して保持する。また、計算制御手段２１１は、計算条件データ記憶部２２に記憶された粒子間距離Ｌ0、クーラン数の上限
値Ｃmax及び粒子速度から時間刻み幅Δｔを算出して保持する。
【００６２】
そして、制御手段２１は、次のステップＳ１−３以降の処理を、タイムステップ「０」から、順次、反復数上限値のタイムステップまで繰り返して行なう。ここでは、既にタイムステップｋにおいて算出された各粒子の位置ｒ（k）、速度ｕ（k）、圧力Ｐ（k）を用
いて、次のタイムステップ（ｋ＋１）における位置ｒ（k+1）、速度ｕ（k+1）、圧力Ｐ（k+1）を計算する。更に、各タイムステップにおける処理では、外力計算と、算出した外
力を用いて各粒子の仮速度及び仮位置とを陽的に計算した後、粒子毎の圧力を陰的に計算し、この圧力を用いて粒子の速度及び位置の修正量を陰的に計算する。
【００６３】
以下、所定のタイムステップｋにおける処理を説明する。ここでは、まず、制御手段２１は、各流体の粒子に作用する外力計算処理を実行する（ステップＳ１−３）。この外力計算処理では、各粒子Ａ１のそれぞれについて外力計算を行なう。この外力計算を行なう１の粒子Ａ１を着目粒子ａとし、この着目粒子ａに作用する外力の算出について、以下に説明する。
【００６４】
まず、制御手段２１の計算制御手段２１１が、着目粒子ａに対して、外力を作用させている粒子を特定する。具体的には、計算制御手段２１１は、着目粒子ａの位置ｒ（k）デ
ータと、計算条件データ記憶部２２に記憶されている第３半径データとから、着目粒子ａに対して外力計算を行なう粒子を特定するための境界領域を算出する。そして、計算制御手段２１１は、算出した境界領域と、他の粒子の位置ｒ（k）とを比較する。この比較結
果から、着目粒子ａの境界領域内にある粒子を、着目粒子ａに対して外力を作用させている粒子として特定する。例えば、図４は、着目粒子ａに対する、粒子Ａ１，Ａ２の影響を示した例である。この場合、計算制御手段２１１は、境界領域が点線で示される円により特定する。そして、この境界領域に含まれる粒子を、着目粒子ａに対して外力を作用させている粒子と特定する。
【００６５】
そして、計算制御手段２１１は、外力を作用させている粒子を特定するデータ（例えば各粒子の位置データ)を外力計算手段２１２に供給する。
外力計算手段２１２は、特定した粒子Ａ１，Ａ２の位置ｒ（k）と着目粒子ａの位置ｒ
（k）とから、着目粒子ａに対する各粒子Ａ１，Ａ２の相互距離Ｌを算出する。
【００６６】
そして、外力計算手段２１２は、第３半径Ｒ3内にあると特定した粒子Ａ１，Ａ２が着
目粒子ａに作用する外力ｆを算出する。
この場合、外力計算手段２１２は、流体の粒子Ａ１が、同じ種類の着目粒子ａに対して作用する外力ｆllを、上記（３）式を用いて算出する。具体的には、外力計算手段２１２は、計算条件データ記憶部２２から表面張力係数σを取得し、上記（４）式を用いて、ポテンシャルパラメータαを算出し、計算条件データ記憶部２２から粒子間距離Ｌ０を取得して第３半径Ｒ3を算出する。更に、外力計算手段２１２は、粒子間距離Ｌ０と、先に算
出した第３半径Ｒ3、相互距離Ｌ及びポテンシャルパラメータαとを（３）式に代入する
ことにより、各粒子Ａ１から着目粒子ａに作用する外力ｆllを算出する。
【００６７】
更に、外力計算手段２１２は、壁の粒子Ａ２が、異なる種類の着目粒子ａに対して作用する外力ｆlwを、上記（３）式及び（５）式を用いて算出する。具体的には、外力計算手段２１２は、計算条件データ記憶部２２から、接触角θ及びポテンシャルパラメータαを取得する。そして、外力計算手段２１２は、ポテンシャルパラメータαと、先に算出した相互距離Ｌとを（３）式に代入して外力ｆllを算出した後、この外力ｆllと接触角θとを（５）式に代入することにより、各粒子Ａ２から着目粒子ａに作用する外力ｆlwを算出する。
【００６８】
更に、外力計算手段２１２は、算出した各粒子Ａ１によって作用される外力ｆllと、各粒子Ａ２によって作用される外力ｆlwとを合計して、着目粒子ａの外力ｆを算出する。このとき、外力ｆll，ｆlwがベクトル量であるため、外力計算手段２１２は、ベクトルの計算により外力ｆを算出する。
【００６９】
以上の計算を、各流体の粒子Ａ１を順次、着目粒子ａとして、すべての粒子Ａ１について行ない、流体の粒子Ａ１に作用する各外力ｆを算出する。
次に、制御手段２１は、算出した外力ｆを用いて、各粒子Ａ１，Ａ２の仮速度ｕ*の算
出処理を実行する（ステップＳ１−４）。このとき、制御手段２１の粒子移動計算手段２１３は、上記（６）式を用いて、仮速度ｕ*を算出する。具体的には、粒子移動計算手段
２１３は、重力加速度ｇ及び動粘性係数νを計算条件データ記憶部２２から取得し、保持している時間刻みΔｔと、タイムステップｋにおける速度ｕ（k）のデータを粒子データ
記憶部２３から取得する。そして、粒子移動計算手段２１３は、取得したこれらのデータと、上記ステップＳ１−３において算出した外力ｆとを、上記（６）式に代入して、仮速度ｕ*を算出する。なお、この（６）式の算出においても、粒子移動計算手段２１３は、
ベクトルの計算を行なう。
【００７０】
次に、制御手段２１は、各粒子Ａ１，Ａ２の仮位置ｒ*の算出処理を実行する（ステッ
プＳ１−５）。具体的には、粒子移動計算手段２１３は、タイムステップｋにおける位置ｒ（k）を粒子データ記憶部２３から取得し、これと時間刻み幅Δｔと、ステップＳ１−
５で算出したｕ*とを、（７）式に代入して、各粒子Ａ１，Ａ２の仮位置ｒ*を算出する。
【００７１】
次に、制御手段２１は、陽的に計算した各粒子の仮速度ｕ*と仮位置ｒ*を用いて、圧力のポアソン方程式の計算処理を実行する（ステップＳ１−６）。具体的には、粒子配置修正手段２１４は、ステップ粒子移動計算手段２１３がステップＳ１−５において算出した仮位置ｒ*から、ステップＳ１−５が終了した時点での粒子数密度ｎ*を算出する。更に、粒子配置修正手段２１４は、算出した粒子数密度ｎ*と、保持している粒子数密度ｎ0とを用いて、上記（８）式の左辺を、流体の粒子Ａ１の各位置におけるラプラシアンモデルで離散化し、Ｐ（k+1）に対する連立１次方程式を算出する。
【００７２】
この離散化の際に、粒子配置修正手段２１４は、ラプラシアンモデルを適用して圧力変数ベクトルに対する係数行列を作成する。この場合、粒子配置修正手段２１４は、計算条件データ記憶部２２に記憶された第２半径データを取得し、このデータに応じたモデルの方程式を用いる。また、粒子配置修正手段２１４は、壁近傍の流体の粒子Ａ１や内側の壁の粒子Ａ２では、圧力を計算しない外側の壁の粒子Ａ２との係数をすべて「０」に設定する。ただし、係数行列の対角成分も修正し、非対角成分の和の絶対値と、対角成分の絶対値が等しくなるようにし、係数行列の対称性を維持する調整を実行する。なお、このようにすることで、粒子配置修正手段２１４は、壁面で圧力勾配「０」のノイマン境界条件を満たした計算を行なえる。
【００７３】
更に、粒子配置修正手段２１４は、計算条件データ記憶部２２に記憶されたパラメータβを取得し、このパラメータβと初期の粒子数密度n（0）との積と、算出した粒子数密度ｎ*とを比較する。そして、粒子配置修正手段２１４は、粒子数密度ｎ*が、β・n（0）より小さい粒子については、圧力のポアソン方程式においてＰ（k+1）＝０と設定する。
【００７４】
そして、上述の連立１次方程式を解くことにより、次のタイムステップ（ｋ＋１）の各粒子の圧力ｐ（k+1）が算出される。なお、この場合、粒子配置修正手段２１４は、計算
条件データ記憶部２２に記録されている圧力計算条件データが用いられる。
【００７５】
次に、制御手段２１は、圧力勾配の計算を行なうことにより、各粒子の速度と位置の修正処理を実行する（ステップＳ１−７）。具体的には、粒子配置修正手段２１４が、ステップＳ１−６で算出したタイムステップ（ｋ＋１）における各粒子の圧力ｐ（k+1）を、
上記（９）式に代入して、速度の修正量ｕ'を算出する。この場合、粒子配置修正手段２
１４は、内側の壁粒子では圧力を計算し、この圧力を流体の粒子の圧力勾配の計算に用いる。ただし、粒子配置修正手段２１４は、壁の粒子Ａ２自体の圧力勾配項は計算しない。
【００７６】
そして、粒子配置修正手段２１４は、算出した速度の修正量ｕ'を仮速度ｕ*に加算する
ことにより、修正した各粒子の速度ｕ（k+1）を算出する。
更に、粒子配置修正手段２１４は、算出した速度の修正量ｕ'と上記（１０）式とを用
いて、流体の粒子Ａ１のそれぞれについて、修正した位置ｒ（k+1）を算出する。壁の粒
子Ａ２の位置は常に固定されているため、粒子配置修正手段２１４は、壁の粒子Ａ２の位置ｒ（k+1）は計算しない。
【００７７】
次に、制御手段２１は、位置ｒ（k+1）、速度ｕ（k+1）及び圧力Ｐ（k+1）の出力処理
を実行する（ステップＳ１−８）。具体的には、計算制御手段２１１は、上記ステップＳ１−３〜Ｓ１−７の処理において算出したタイムステップ（ｋ＋１）における各粒子Ａ１，Ａ２の位置ｒ（k+1）、速度ｕ（k+1）及び圧力Ｐ（k+1）を、このタイムステップ（ｋ
＋１）と関連付けて、粒子データ記憶部２３に記憶する。更に、本実施形態では、計算制御手段２１１は、このタイムステップ（ｋ＋１）における各粒子Ａ１，Ａ２の配置を、出力手段１５において出力する。具体的には、各位置ｒ（k）に表示されていた各粒子Ａ１
，Ａ２のオブジェクトを、位置ｒ（k+1）に変更して表示し直す。この場合、各粒子Ａ１
，Ａ２の速度や圧力を、それぞれ配置された粒子Ａ１，Ａ２のオブジェクト上の矢印の大きさや色分けにより表示させるようにしてもよい。以上により、タイムステップｋにおける処理が終了する。
【００７８】
このようにしてタイムステップｋにおける処理（ステップＳ１−３〜Ｓ１−８）が終了すると、制御手段２１は、終了したタイムステップに「１」を加算した次のタイムステップ（ｋ＋１）について、上記ステップＳ１−３以降の処理を繰り返し行なう。そして、予め定めたタイムステップが終了するまで、上記処理（ステップＳ１−３〜Ｓ１−８）を繰り返し実行する。具体的には、計算制御手段２１１は、計算条件データ記憶部２２に記憶されている反復上限値データを取得し、この反復上限値と、処理が終了したタイムステップｋとを比較する。そして、計算制御手段２１１は、処理が終了したタイムステップｋに「１」を加算した値が反復上限値を超えた場合には、流体解析処理を終了する。
【００７９】
本実施形態によれば、以下のような効果を得ることができる。
・ 本実施形態では、外力計算手段２１２は、第３半径Ｒ3に含まれる流体の粒子Ａ１
から受ける外力ｆllは（３）式を用いて算出する。また、外力計算手段２１２は、同じ第３半径Ｒ3に含まれる壁の粒子Ａ２から受ける外力ｆlwは（５）式を用いて算出する。こ
の（５）式は、（３）式で算出した外力ｆllに対して、（１＋ｃｏｓθ）／２を乗算した算式である。この（１＋ｃｏｓθ）／２という係数は、液相の界面エネルギと壁−液相の界面エネルギ関係式であるヤングの式から算出される係数である。このため、外力計算手段２１２は、エネルギ関係式から壁の粒子Ａ２の作用力を算出している。そして、粒子移動計算手段２１３は、これら外力ｆlw，ｆllの合計値を用いて、各粒子の仮速度ｕ*を算
出し（ステップＳ１−４）、この仮速度ｕ*から仮位置ｒ*を算出する（ステップＳ１−５）。更に、粒子配置修正手段２１４は、仮速度ｕ*、仮位置ｒ*及びポアソン方程式の計算を行ない（ステップＳ１−６）、圧力勾配の計算を行なって各粒子の速度と位置を修正し（ステップＳ１−７）、次のタイムステップの各粒子の位置、速度及び圧力を出力する。そして、このようなタイムステップの処理を、反復数上限値まで繰り返す。
【００８０】
従って、本実施形態の制御手段２１は、表面張力モデルにおいて、エネルギ関係式を用いて壁の粒子Ａ２からの作用力を算出している。このため、表面張力モデルの１つである濡れ性モデルであってもエネルギ関係式が成り立つため、濡れ性モデルも表面張力モデルとして統一的に取り扱うことができる。この場合、外力計算手段２１２は、粒子数密度の減少分から算出される界面の曲率を用いずに第２の粒子の作用力を算出するので、空間解像度が粗くても、流体の挙動の解析を的確に行なうことができる。
【００８１】
更に、エネルギの関係式は、流体や壁がどのような状態にあっても成立する。外力計算
手段２１２は、このエネルギの関係式を用いて第２の粒子の作用力を算出しているため、界面が動く場合など種々の条件下で用いることができ、汎用性が高い。
【００８２】
・ 本実施形態では、上記（３）式として、他の粒子との距離に応じたポテンシャルエネルギの算出式を用いている。このため、エネルギの関係式である（５）式と整合がよく、流体の挙動の解析を、更に的確に行なうことができる。
【００８３】
また、上記実施形態は、以下のように変更してもよい。
○ 上記実施形態において、外力計算手段２１２は、外力ｆllを計算するポテンシャル関数として、上記（３）式を用いた。これに限らず、外力ｆllを算出する式は、他の関数を用いてもよい。この場合、（５）式との関係上、エネルギ関係式で外力ｆllを表現した関数を用いるのがよい。
【００８４】
○ 上記実施形態において、計算条件データ記憶部２２に記録された反復数上限値を用いて、制御手段２１は計算処理を終了させたが、計算処理の終了方法は、これに限定されるものではない。例えば、終了時間を記録して、これに対応させた反復数上限値により終了させてもよい。また、制御手段２１は、入力手段１０を介して終了信号が入力された段階で流体解析処理を終了するにしてもよい。
【００８５】
○ 上記実施形態においては、制御手段２１は、圧縮率がきわめて小さく圧縮性を考慮しない流体の挙動について解析を行なった。これに限らず、解析する対象によっては、流体の圧縮性を考慮してもよい。この場合、基準となる初期の密度及び圧力から、圧縮性を考慮した公知の修正を行なう。これにより、圧縮性が小さい流体において速い現象をも取り扱うことができる。
【図面の簡単な説明】
【００８６】
【図１】実施形態におけるシステムの概略図。
【図２】実施形態における処理手順を説明するための流れ図。
【図３】接触角θを説明するための説明図。
【図４】各粒子から受ける外力の算出を説明するための説明図。
【図５】（４)式の導出過程を説明するための説明図。
【符号の説明】
【００８７】
θ…接触角、Ａ１…第１の粒子、Ａ２…第２の粒子、Ｐ，ｐ…圧力、Ｒ3…所定範囲と
しての第３半径、ｆll，ｆlw…作用力としての外力、ｒ…位置、ｒ*…仮位置、ｕ…速度
、ｕ*…仮速度、２０…流体解析手段、２１…制御手段、２２…計算条件データ記憶手段
としての計算条件データ記憶部、２３…粒子データ記憶手段としての粒子データ記憶部、２１１…特定手段としての計算制御手段、２１２…外力計算手段、２１３…粒子移動計算手段、２１４…粒子配置修正手段。
【特許請求の範囲】
【請求項１】
流体を表した第１の粒子毎に、速度、位置、圧力に関するデータと、前記流体に接触する壁を表した第２の粒子毎に、速度、位置、圧力に関するデータとを記録する粒子データ記憶手段と、
前記流体と前記壁との接触角θに関するデータと、各粒子について他の粒子からの作用力を受ける作用範囲に関するデータとを記録した計算条件データ記憶手段と、
各粒子の挙動を計算する制御手段とから構成され、流体の挙動を解析する粒子法を用いた解析装置であって、
前記制御手段が、
所定のタイムステップにおいて、各第１の粒子の配置において前記作用範囲に含まれる他の粒子を特定する特定手段と、
特定された他の粒子のうち第１の粒子については、第１の関数により第１の作用力を算出し、特定された他の粒子のうち第２の粒子については、前記第１の関数に対して（１＋ｃｏｓθ）／２を乗算した第２の関数により第２の作用力を算出する作用力算出手段と、
前記各作用力により算出した各粒子の仮速度を用いて、仮配置を決定する粒子移動計算手段と、
前記仮配置を用いて、前記粒子の数密度を一定にする計算により、次回タイムステップの圧力を算出し、
この圧力を用いて算出した圧力勾配から各粒子の修正後の速度を算出し、この速度を用いて各粒子の位置を修正し、前記粒子データ記憶手段に記録する粒子配置修正手段と、
タイムステップを進めて、前記特定手段、前記粒子移動計算手段及び前記粒子配置修正手段を用いた各粒子の挙動の計算処理を繰り返す計算制御手段と
を備えたことを特徴とする流体解析装置。
【請求項２】
前記第１の関数は、作用力を及ぼす前記他の粒子との距離に応じたポテンシャルエネルギにより算出する関数であることを特徴とする請求項１に記載の流体解析装置。
【請求項３】
流体を表した第１の粒子毎に、速度、位置、圧力に関するデータと、前記流体に接触する壁を表した第２の粒子毎に、速度、位置、圧力に関するデータとを記録する粒子データ記憶手段と、
前記流体と前記壁との接触角θに関するデータと、各粒子について他の粒子からの作用力を受ける作用範囲に関するデータとを記録した計算条件データ記憶手段と、
各粒子の挙動を計算する制御手段とから構成された解析装置を用いて、粒子法による流体の挙動を解析する方法であって、
前記制御手段が、
所定のタイムステップにおいて、各第１の粒子の配置において前記作用範囲に含まれる他の粒子を特定する特定段階と、
特定された他の粒子のうち第１の粒子については、第１の関数により第１の作用力を算出し、特定された他の粒子のうち第２の粒子については、前記第１の関数に対して（１＋ｃｏｓθ）／２を乗算した第２の関数により第２の作用力を算出する作用力算出段階と、
前記各作用力により算出した各粒子の仮速度を用いて、仮配置を決定する粒子移動計算段階と、
前記仮配置を用いて、前記粒子の数密度を一定にする計算により、次回タイムステップの圧力を算出し、
この圧力を用いて算出した圧力勾配から各粒子の修正後の速度を算出し、この速度を用いて各粒子の位置を修正し、前記粒子データ記憶手段に記録する粒子配置修正段階と、
タイムステップを進めて、前記特定段階、前記粒子移動計算段階及び前記粒子配置修正段階による各粒子の挙動の計算処理を繰り返す計算制御段階と
を実行することを特徴とする流体解析方法。
【請求項４】
前記第１の関数は、作用力を及ぼす前記他の粒子との距離に応じたポテンシャルエネルギにより算出する関数であることを特徴とする請求項３に記載の流体解析方法。
【請求項５】
流体を表した第１の粒子毎に、速度、位置、圧力に関するデータと、前記流体に接触する壁を表した第２の粒子毎に、速度、位置、圧力に関するデータとを記録する粒子データ記憶手段と、
前記流体と前記壁との接触角θに関するデータと、各粒子について他の粒子からの作用力を受ける作用範囲に関するデータとを記録した計算条件データ記憶手段と、
各粒子の挙動を計算する制御手段とから構成された解析装置を用いて、粒子法による流体の挙動を解析するプログラムであって、
前記制御手段を、
所定のタイムステップにおいて、各第１の粒子の配置において前記作用範囲に含まれる他の粒子を特定する特定手段、
特定された他の粒子のうち第１の粒子については、第１の関数により第１の作用力を算出し、特定された他の粒子のうち第２の粒子については、前記第１の関数に対して（１＋ｃｏｓθ）／２を乗算した第２の関数により第２の作用力を算出する作用力算出手段、
前記各作用力により算出した各粒子の仮速度を用いて、仮配置を決定する粒子移動計算手段、及び
前記仮配置を用いて、前記粒子の数密度を一定にする計算により、次回タイムステップの圧力を算出し、
この圧力を用いて算出した圧力勾配から各粒子の修正後の速度を算出し、この速度を用いて各粒子の位置を修正し、前記粒子データ記憶手段に記録する粒子配置修正手段と、
タイムステップを進めて、前記特定手段、前記粒子移動計算手段及び前記粒子配置修正手段を用いた各粒子の挙動の計算処理を繰り返す計算制御手段
として機能させることを特徴とする流体解析プログラム。
【請求項６】
前記第１の関数は、作用力を及ぼす前記他の粒子との距離に応じたポテンシャルエネルギにより算出する関数であることを特徴とする請求項５に記載の流体解析プログラム。
【請求項１】
流体を表した第１の粒子毎に、速度、位置、圧力に関するデータと、前記流体に接触する壁を表した第２の粒子毎に、速度、位置、圧力に関するデータとを記録する粒子データ記憶手段と、
前記流体と前記壁との接触角θに関するデータと、各粒子について他の粒子からの作用力を受ける作用範囲に関するデータとを記録した計算条件データ記憶手段と、
各粒子の挙動を計算する制御手段とから構成され、流体の挙動を解析する粒子法を用いた解析装置であって、
前記制御手段が、
所定のタイムステップにおいて、各第１の粒子の配置において前記作用範囲に含まれる他の粒子を特定する特定手段と、
特定された他の粒子のうち第１の粒子については、第１の関数により第１の作用力を算出し、特定された他の粒子のうち第２の粒子については、前記第１の関数に対して（１＋ｃｏｓθ）／２を乗算した第２の関数により第２の作用力を算出する作用力算出手段と、
前記各作用力により算出した各粒子の仮速度を用いて、仮配置を決定する粒子移動計算手段と、
前記仮配置を用いて、前記粒子の数密度を一定にする計算により、次回タイムステップの圧力を算出し、
この圧力を用いて算出した圧力勾配から各粒子の修正後の速度を算出し、この速度を用いて各粒子の位置を修正し、前記粒子データ記憶手段に記録する粒子配置修正手段と、
タイムステップを進めて、前記特定手段、前記粒子移動計算手段及び前記粒子配置修正手段を用いた各粒子の挙動の計算処理を繰り返す計算制御手段と
を備えたことを特徴とする流体解析装置。
【請求項２】
前記第１の関数は、作用力を及ぼす前記他の粒子との距離に応じたポテンシャルエネルギにより算出する関数であることを特徴とする請求項１に記載の流体解析装置。
【請求項３】
流体を表した第１の粒子毎に、速度、位置、圧力に関するデータと、前記流体に接触する壁を表した第２の粒子毎に、速度、位置、圧力に関するデータとを記録する粒子データ記憶手段と、
前記流体と前記壁との接触角θに関するデータと、各粒子について他の粒子からの作用力を受ける作用範囲に関するデータとを記録した計算条件データ記憶手段と、
各粒子の挙動を計算する制御手段とから構成された解析装置を用いて、粒子法による流体の挙動を解析する方法であって、
前記制御手段が、
所定のタイムステップにおいて、各第１の粒子の配置において前記作用範囲に含まれる他の粒子を特定する特定段階と、
特定された他の粒子のうち第１の粒子については、第１の関数により第１の作用力を算出し、特定された他の粒子のうち第２の粒子については、前記第１の関数に対して（１＋ｃｏｓθ）／２を乗算した第２の関数により第２の作用力を算出する作用力算出段階と、
前記各作用力により算出した各粒子の仮速度を用いて、仮配置を決定する粒子移動計算段階と、
前記仮配置を用いて、前記粒子の数密度を一定にする計算により、次回タイムステップの圧力を算出し、
この圧力を用いて算出した圧力勾配から各粒子の修正後の速度を算出し、この速度を用いて各粒子の位置を修正し、前記粒子データ記憶手段に記録する粒子配置修正段階と、
タイムステップを進めて、前記特定段階、前記粒子移動計算段階及び前記粒子配置修正段階による各粒子の挙動の計算処理を繰り返す計算制御段階と
を実行することを特徴とする流体解析方法。
【請求項４】
前記第１の関数は、作用力を及ぼす前記他の粒子との距離に応じたポテンシャルエネルギにより算出する関数であることを特徴とする請求項３に記載の流体解析方法。
【請求項５】
流体を表した第１の粒子毎に、速度、位置、圧力に関するデータと、前記流体に接触する壁を表した第２の粒子毎に、速度、位置、圧力に関するデータとを記録する粒子データ記憶手段と、
前記流体と前記壁との接触角θに関するデータと、各粒子について他の粒子からの作用力を受ける作用範囲に関するデータとを記録した計算条件データ記憶手段と、
各粒子の挙動を計算する制御手段とから構成された解析装置を用いて、粒子法による流体の挙動を解析するプログラムであって、
前記制御手段を、
所定のタイムステップにおいて、各第１の粒子の配置において前記作用範囲に含まれる他の粒子を特定する特定手段、
特定された他の粒子のうち第１の粒子については、第１の関数により第１の作用力を算出し、特定された他の粒子のうち第２の粒子については、前記第１の関数に対して（１＋ｃｏｓθ）／２を乗算した第２の関数により第２の作用力を算出する作用力算出手段、
前記各作用力により算出した各粒子の仮速度を用いて、仮配置を決定する粒子移動計算手段、及び
前記仮配置を用いて、前記粒子の数密度を一定にする計算により、次回タイムステップの圧力を算出し、
この圧力を用いて算出した圧力勾配から各粒子の修正後の速度を算出し、この速度を用いて各粒子の位置を修正し、前記粒子データ記憶手段に記録する粒子配置修正手段と、
タイムステップを進めて、前記特定手段、前記粒子移動計算手段及び前記粒子配置修正手段を用いた各粒子の挙動の計算処理を繰り返す計算制御手段
として機能させることを特徴とする流体解析プログラム。
【請求項６】
前記第１の関数は、作用力を及ぼす前記他の粒子との距離に応じたポテンシャルエネルギにより算出する関数であることを特徴とする請求項５に記載の流体解析プログラム。
【図１】
【図２】
【図３】
【図４】
【図５】
【図２】
【図３】
【図４】
【図５】
【公開番号】特開２００８−１１１６７５（Ｐ２００８−１１１６７５Ａ）
【公開日】平成２０年５月１５日（２００８．５．１５）
【国際特許分類】
【出願番号】特願２００６−２９３２９８（Ｐ２００６−２９３２９８）
【出願日】平成１８年１０月２７日（２００６．１０．２７）
【国等の委託研究の成果に係る記載事項】（出願人による申告）平成１７〜１８年度独立行政法人新エネルギー・産業技術総合開発機構固体高分子型燃料電池実用化戦略的技術開発委託研究、産業活力再生特別措置法第３０条の適用を受ける特許出願
【出願人】（５９２１３１９０６）みずほ情報総研株式会社 (187)
【出願人】（５０４１３７９１２）国立大学法人 東京大学 (1,942)
【Ｆターム（参考）】
【公開日】平成２０年５月１５日（２００８．５．１５）
【国際特許分類】
【出願日】平成１８年１０月２７日（２００６．１０．２７）
【国等の委託研究の成果に係る記載事項】（出願人による申告）平成１７〜１８年度独立行政法人新エネルギー・産業技術総合開発機構固体高分子型燃料電池実用化戦略的技術開発委託研究、産業活力再生特別措置法第３０条の適用を受ける特許出願
【出願人】（５９２１３１９０６）みずほ情報総研株式会社 (187)
【出願人】（５０４１３７９１２）国立大学法人 東京大学 (1,942)
【Ｆターム（参考）】
[ Back to top ]