(57)【特許請求の範囲】
【請求項１】
計算解析の対象となる連続体が複数の粒子により自由表面を伴って集合的にモデリングされるところの粒子法の計算処理の過程において、ある粒子(i)が自由表面粒子に属するかどうかを判定する方法であって、
計算処理の過程における以前のタイムステップにおけるメモリに記憶されている自由表面識別情報を参照して、ある粒子(i)を基準にした所定の範囲(R)内に、自由表面に属していた粒子が存在しているかどうかを判定するステップと、
判定するステップにおいて、自由表面に属していた粒子が存在していると判定されたら「近自由表面」であるとし、存在していないと判定されたら「十分内部」であるとして、以後のタイムステップの参照のために近自由表面識別情報としてメモリに記憶するステップとを有する、
方法。
【請求項２】
計算解析の対象となる連続体の複数の粒子の全てについて、請求項１に記載の全てのステップを繰り返すことを特徴とする、
請求項１に記載の方法。
【請求項３】
さらに、
メモリに記憶されている近自由表面識別情報を参照して、ある粒子(i)が十分内部であると判定された場合には、「内部」であるとして、以後のタイムステップの参照のために自由表面識別情報としてメモリに記憶するステップと、
を有する、請求項１に記載の方法。
【請求項４】
さらに、
メモリに記憶されている近自由表面識別情報を参照して、ある粒子(i)が近自由表面であると判定された場合には、以前のタイムステップにおける粒子数密度と所定の粒子数密度の閾値とを比較するステップと、
比較するステップにおいて、所定の粒子数密度の閾値よりも低い場合には、「自由表面」であるとして、以後のタイムステップの参照のために自由表面識別情報としてメモリに記憶するステップと、
比較するステップにおいて、所定の粒子数密度の閾値よりも低くない場合には、「内部」であるとして、以後のタイムステップの参照のために自由表面識別情報としてメモリに記憶するステップとを有する、請求項３に記載の方法。
【請求項５】
計算解析の対象となる連続体が複数の粒子により自由表面を伴って集合的にモデリングされるところの粒子法の計算処理の過程において、ある粒子(i)が自由表面粒子に属するかどうかを判定する装置であって、
計算処理の過程における以前のタイムステップにおけるメモリに記憶されている自由表面識別情報を参照して、ある粒子(i)を基準にした所定の範囲(R)内に、自由表面に属していた粒子が存在しているかどうかを判定する手段と、
判定する手段において、自由表面に属していた粒子が存在していると判定されたら「近自由表面」であるとし、存在していないと判定されたら「十分内部」であるとして、以後のタイムステップの参照のために近自由表面識別情報としてメモリに記憶する手段とを有する、
装置。
【請求項６】
計算解析の対象となる連続体の複数の粒子の全てについて、請求項５に記載の各手段の処理を繰り返すことを特徴とする、
請求項５に記載の装置。
【請求項７】
さらに、
メモリに記憶されている近自由表面識別情報を参照して、ある粒子(i)が十分内部であると判定された場合には、「内部」であるとして、以後のタイムステップの参照のために自由表面識別情報としてメモリに記憶する手段と、
を有する、請求項５に記載の装置。
【請求項８】
さらに、
メモリに記憶されている近自由表面識別情報を参照して、ある粒子(i)が近自由表面であると判定された場合には、以前のタイムステップにおける粒子数密度と所定の粒子数密度の閾値とを比較する手段と、
比較する手段において、所定の粒子数密度の閾値よりも低い場合には、「自由表面」であるとして、以後のタイムステップの参照のために自由表面識別情報としてメモリに記憶する手段と、
比較する手段において、所定の粒子数密度の閾値よりも低くない場合には、「内部」であるとして、以後のタイムステップの参照のために自由表面識別情報としてメモリに記憶する手段とを有する、請求項７に記載の装置。
【請求項９】
計算解析の対象となる連続体が複数の粒子により自由表面を伴って集合的にモデリングされるところの粒子法の計算処理の過程において、ある粒子(i)が自由表面粒子に属するかどうかを判定するプログラムであって、
計算処理の過程における以前のタイムステップにおけるメモリに記憶されている自由表面識別情報を参照して、ある粒子(i)を基準にした所定の範囲(R)内に、自由表面に属していた粒子が存在しているかどうかを判定するステップと、
判定するステップにおいて、自由表面に属していた粒子が存在していると判定されたら「近自由表面」であるとし、存在していないと判定されたら「十分内部」であるとして、以後のタイムステップの参照のために近自由表面識別情報としてメモリに記憶するステップと有しており、これらのステップをコンピュータに実行させる、
プログラム。
【発明の詳細な説明】
【技術分野】
【０００１】
本願発明は、粒子法を用いて流体等の連続体の挙動を数値的に解析するにあたり、流体等の連続体を近似的に表現する粒子のうち、自由表面に位置するものを正確に判定するための方法、装置及びプログラムに関する。
【背景技術】
【０００２】
粒子法は、連続体を複数の粒子の集まりで近似的に表現し、粒子の集まりが連続体の挙動を模倣するように、これらの粒子の間に代数的な関係式を定義し、数値計算によって解析を行う手法である。より一般化した言い方をするのであれば、計算解析の対象となる連続体が複数の粒子により自由表面を伴って集合的にモデリングされる。粒子法はメッシュ生成が不要であり、解析対象である流体の大変形にも容易に対応できるという利点がある。
【０００３】
流体に対する粒子法としては、圧縮性流体の挙動を陽解法で計算するSPH(Smoothed Particle Hydrodynamics)法や、非圧縮性流体の挙動を半陰解法で計算するMPS(Moving ParticleSemi-implicit)法やISPH法が知られている。非圧縮性流体に対する粒子法では、その計算処理の過程において、非圧縮性の条件を反映するように定義されたPoisson方程式の解として流体の圧力が計算される。しかし、実際に計算された圧力に不自然な擾乱が含まれてしまうことが、粒子法の実用化に向けた大きな問題となっている。
【０００４】
Poisson方程式の解として流体の圧力を計算するためには、流体を近似的に表現するある粒子が自由表面に位置するか否かを判定する必要が生じる。自由表面に位置すると判定された粒子に対しては、圧力の値を0に固定するという、Poisson方程式に境界条件を課すことに相当する処理が行われる。但し、Poisson方程式の解は、境界条件の違いによって大きく変化する。即ち、流体を近似的に表現するある粒子が自由表面に位置するか否かを正確に判定できなければ、粒子法で計算される圧力の分布に不自然な擾乱が引き起こされる。
【０００５】
ある粒子(i)が自由表面に位置するか否かを判定するための方法として最も代表的なものでは、ある粒子(i)の「粒子数密度」nｉ が、一定値n０に対して次式（１）の条件を満たすか否かで判定する。
【数１】
0005576309.tif 000002
粒子数密度 nｉは、ある粒子(i)を基準とした所定の範囲(r e)内に存在する他の粒子に対して定義される重みの総和として定義される。n０は基準とする粒子数密度であり、βはモデルによって定められる1未満の定数である。計算過程において非圧縮性が厳密に満たされていれば、流体の内部に存在する粒子の粒子数密度は、基準とする粒子数密度n０と一致する。この方法によって判定を行うのは、流体の外部には粒子が配置されないことから、自由表面に位置する粒子に対しては、重みの総和計算に影響する粒子の数が少なくなり、結果としてその粒子数密度が低下すると考えられるからである。他にも、本質的に同様の根拠によって定式化された閾値に基づく判定法がいくつか知られている。
【０００６】
しかしながら、このような閾値に基づく判定法では、自由表面に位置する粒子以外の粒子(解析対象の内部の粒子)についても自由表面粒子として誤判定される場合があった。これには、実際の粒子法の計算では、流体の内部に位置する粒子に対して、厳密には粒子数密度が一定に保たれず、ずれが生じてしまうことが大きく影響している。特に、本来は自由表面に位置すると判定されるべき粒子の粒子数密度よりも、流体の内部に位置する粒子の方が粒子数密度が小さくなってしまうことも少なくないため、このような方法を単独で用いることには、本質的な困難があるといえる。
【０００７】
粒子数密度に対する閾値判定を単独で適用するのではなく、粒子の分布の非対称性を同時に考慮する方法も知られている。この方法では、ある粒子を基準とした所定範囲内の他の粒子の分布の非対称性の尺度となる数値を定義し、その数値の閾値に基づく判定の結果が非対称であるとされ、かつ、粒子数密度が低いと判定されるものを自由表面に位置する粒子であると判定する。流体の外部には粒子が配置されないことから、自由表面に位置する粒子は単に粒子数密度が低下するだけでなく、その周辺には偏りを持って他の粒子が分布すると考えられることが、この判定法の根拠である。これは、特許文献1として公知である。
【０００８】
しかし、非対称性を考慮したこの判定方法を用いた場合でも、必ずしも正確に自由表面に位置する粒子が特定されるわけではない。例えば、全ての粒子の相対的な位置関係を維持し、かつ、座標系を固定したまま、ある粒子を中心にして適当な角度だけ回転させただけで、ある粒子における非対称性の尺度として用いられる数値は変化してしまう。これでは、本来は同じ判定結果であるべきものが、状況によって異なる判定結果となってしまうため、判定のための閾値を設定することは困難である。また、判定に用いるパラメータを変化させた場合に判定結果が受ける影響も非常に大きく、非対称性を考慮した判定法自体がほとんど有効に機能しない場合も少なくない。
【０００９】
一般に、流体の内部に位置するにもかかわらず自由表面に位置すると誤判定された粒子がある場合、その粒子が流体のより内部に位置するものであるほど、計算される圧力の分布に引き起こす擾乱はより大きくなる。ここまでの説明から、自由表面に位置する粒子を正しく判定しつつ、このような圧力の分布に大きな擾乱を引き起こす誤判定も原理的に防ぐということは、閾値などのパラメータを注意深く設定したとしても、従来手法を用いる限り不可能である。特に、多数の粒子を用いた大規模な計算による解析を行おうとした場合、状況はより深刻であり、大規模計算の本来の目的である高精度な解析を実現することは、極めて難しくなる。
【００１０】
ここまで述べてきたように、ある粒子が自由表面に位置するか否かの判定を正確に行うことは非常に困難であるが、粒子法を用いた解析を行う際には避けて通ることができない問題である。従って、従来手法のように簡便な指標で自由表面に位置する粒子を特徴付けつつ、圧力に特に大きな擾乱を引き起こす、十分内部の粒子を自由表面であるとしてしまう誤判定を原理的に防ぐこともできるような自由表面判定法を構成することができれば、粒子法の実用化に対して大きな進展を与えるものになると考えられる。
【先行技術文献】
【特許文献】
【００１１】
【特許文献１】国際公開ＷＯ／２０１０／０３２６５６号公報  −  特許協力条約（ＰＣＴ）における国際出願日から１年６月後の強制的な公開  −
【特許文献２】特開平７−３３４４８４号公報    −  日本国における特許出願日から１年６月後の強制的な公開  −
【発明の概要】
【発明が解決しようとする課題】
【００１２】
本願発明の目的は、以上のような背景に鑑みてなされたものであり、流体の自由表面に位置する粒子を正確に判定することにある。
【課題を解決するための手段】
【００１３】
本願発明の方法では、計算解析の対象となる連続体が複数の粒子により自由表面を伴って集合的にモデリングされるところの粒子法（ＰＭ）の計算処理の過程において、ある粒子(i)が自由表面粒子に属するかどうかを判定する方法であって、計算処理の過程における以前のタイムステップにおけるメモリに記憶されている自由表面識別情報を参照して、ある粒子(i)を基準にした所定の範囲(R)内に、自由表面に属していた粒子が存在しているかどうかを判定するステップと、判定するステップにおいて、自由表面に属していた粒子が存在していると判定されたら近自由表面であるとし、存在していないと判定されたら十分内部であるとして、以後のタイムステップの参照のために近自由表面識別情報としてメモリに記憶するステップとを有するものが開示される。各ステップの処理を実行する手段、または、コンピュータにこれらのステップを実行させるプログラムという態様でも実現される。
【００１４】
計算解析の対象となる連続体の複数の粒子の全てについて、ステップを繰り返す。
【００１５】
さらに、
メモリに記憶されている近自由表面識別情報を参照して、ある粒子(i)が十分内部であると判定された場合には、内部であるとして、以後のタイムステップの参照のために自由表面識別情報としてメモリに記憶するステップとを有するものが開示される。
【００１６】
さらに、
メモリに記憶されている近自由表面識別情報を参照して、ある粒子(i)が近自由表面であると判定された場合には、以前のタイムステップにおける粒子数密度が所定の粒子数密度の閾値と比較するステップと、比較するステップにおいて、所定の粒子密度の閾値よりも低い場合には、自由表面であるとして、以後のタイムステップの参照のために自由表面識別情報としてメモリに記憶するステップと、比較するステップにおいて、所定の粒子密度の閾値よりも低くない場合には、内部であるとして、以後のタイムステップの参照のために自由表面識別情報としてメモリに記憶するステップとを有するものが開示される。
【発明の効果】
【００１７】
より正確に自由表面粒子の判定を行うことが可能になる。従来手法と本願発明との比較をした具体的なシミュレーション結果を、実施例において示すことにする。
【図面の簡単な説明】
【００１８】
【図１】図１は、本願発明の実施形態に係る流体解析装置を示す概略図である。
【図２Ａ】図２Ａは、粒子法の一例であるMPS法の流体解析方法の概要を説明する図である。
【図２Ｂ】図２Ｂは、境界条件設定の概要を説明する図である。
【図３】図３は、粒子法（ＰＭ）の計算処理の過程において本願発明が部分的に適用される流体解析方法の手順を示すフローチャートである。
【図４】図４は、本願発明の実施形態に係る近自由表面判定（図３のS2の詳細）を行う手順を示すフローチャートである。
【図５】図５は、本願発明の実施形態に係る自由表面判定（図３のS8の詳細）を行う手順を示すフローチャートである。
【図６】図６は、２次元水柱崩壊という流体の挙動をシミュレーションする概要を示す図である。
【図７】図７は、（ａ）が従来手法を用いた場合の自由表面粒子の分布と圧力の分布を示すシミュレーションの結果の図であって、（ｂ）が本願発明を用いた場合の自由表面粒子の分布と圧力の分布を示すシミュレーションの結果の図である。
【図８】図８は、本願発明における、次のタイムステップでの自由表面粒子の候補をどのあたりの部分から選択するかを示す模式図である。
【発明を実施するための形態】
【００１９】
図１は、本願発明である自由表面粒子判定装置を含む形で構成された、粒子法を用いて流体の挙動を解析する流体解析装置１の概略図である。流体解析装置１は、RAM(Random Access Memory)２、初期設定データベース(DB)３、CPU(CentralProcessing Unit)４、を主な構成要素として備える。
【００２０】
RAM２は、CPU４が実行する流体解析装置１の起動プログラムや流体解析プログラム等のコンピュータプログラム及び各種データを一時的に格納するワーキングエリアを提供する。
【００２１】
初期設定DB３は、（粒子法の一種であるMPS法を例にとると）MPS法によって流体解析を実行するために必要な各種パラメータを格納する。具体的には、メモリとして、標準的なMPS法で必要となる、粒子の位置、粒子の速度、粒子の密度、動粘性係数等の物性に関するデータや、重力加速度等の外力に関するデータ、初期配置における粒子間距離に関するデータ、粒子数密度を求めるための基準半径等の計算データが格納される。更に、本願発明の実施形態ではこれに加えて、流体を構成する各粒子に対して、自由表面に配置されているかどうかを示す自由表面識別情報が記憶される。
【００２２】
このように、ハードウエア資源としての構成、ソフトウエア資源としての構成、または、ハードウエア資源とソフトウエア資源とが協働する構成に成っている流体解析装置１は、流体解析処理を実行する。ソフトウエア資源は、外部からダウンロードするようにすることもでき、コンピュータに実行させるためのプログラムとしても提供することができる。
【００２３】
図２は、図２Ａが、粒子法の一例であるMPS法の流体解析方法の概要を説明する図である。非圧縮性流体の運動方程式を、第１段階と第２段階とに分けて、処理を進めていくものである。図２Ｂは、境界条件設定の概要を説明する図である。この図２Ａおよび図２Ｂで説明している処理自体は、本願発明の出願時点で既に公知であるので、ごく簡単な説明に留めている。
【００２４】
以下、図３に示すフローチャートを参照して、この流体解析装置１の動作について詳しく説明する。図３に示すフローチャートでは、まず、流体解析処理はステップS1の処理へ進む。ステップS1では、CPU４が、初期設定DB３からRAM２内に解析対象となる非圧縮性流体に関する初期の計算データを読み出して取得する。
【００２５】
ステップS2では、近自由表面設定部１１が、近自由表面設定の処理を行う。この処理は、ある粒子(i)または全ての粒子に対して、この時点でその粒子が自由表面であるか、もしくはその粒子(i)を基準にした所定の範囲(R)内に自由表面粒子が存在すると判定されれば「近自由表面」という識別情報を設定し、自由表面粒子が所定の範囲(R)内に存在しないと判定された粒子に対しては、「十分内部」という識別情報を設定する。これは「内部」という識別情報との区別においての命名にすぎず、必ずしも実際の計算解析の対象における空間において「十分」な内部に粒子があることを意味しているものではない。
【００２６】
この識別情報は、近自由表面識別情報としてRAM２（他のどこかのメモリでもよい）に一時的に記憶しておくことができる。この識別情報は、ある粒子または全ての粒子に対して設定されるものであるため、流速の情報や位置の情報が通常取り扱われるのと同様、配列のデータ保持形式で、計算処理のタイムステップを進めていく過程において、これを流体の挙動をシミュレートするために必要な時点まで保持しておけばよい。
【００２７】
このステップS2の処理は、以後の処理である自由表面判定を行う際に、その判定を適用する粒子を近自由表面に属する粒子に制限し、十分内部にある粒子が誤って自由表面であると判定されることを防ぐための処理である。流体解析のタイムステップの発展（ｔ＋ｄｔ）は、計算が数値的に安定に実行できるように、１タイムステップ（ｄｔ）の刻み幅を制限している。
【００２８】
これによれば、流体全体の形状が大幅に変更されることがないので、自由表面から十分に離れ、流体の内部にある粒子が、高々１タイムステップ（ｄｔ）の発展によって、自由表面粒子となることは考えにくい。従って、この段階で十分内部と判定された粒子に対しては、以後の境界条件設定ステップにおいて、自由表面であると判定することを避けることが望ましい。このステップS2のより詳細は、次の図４に示されたフローチャートによって示される処理として実行する。
【００２９】
ステップS21では、ある粒子(i)または全ての粒子について、位置x, yと自由表面識別情報を（メモリを参照して読み出すことで）取得する。ステップS22では、ある粒子(i)に対する自由表面識別情報を取得する。ステップS23では、ある粒子(i)に対する自由表面識別情報が、自由表面であればステップS27へ、自由表面でなければステップS24へ処理を進める。ここで、粒子(i)が自由表面であると識別される場合、以後の自由表面設定を行う際、再び自由表面粒子となる可能性が高い。ステップS24では、まず、粒子(i)と、その他の粒子(j)それぞれの位置x, yから、それらの距離を求める。粒子(i)と粒子(j)との距離とは、次式（２）で定義されるものである。
【数２】
0005576309.tif 000003
【００３０】
ここで、x_i, y_iは粒子(i)のx, y座標、x_j, y_jは粒子(j)のx,y座標である。その距離が所定範囲Ｒ（メモリに記憶されているところの所定の閾値）内にあると判定されれば、即ち、次式（３）が成立すれば、その粒子(j)の自由表面識別情報を参照する。
【数３】
0005576309.tif 000004
【００３１】
そうでなければ、処理をステップS25に進める。更に、粒子(j)の自由表面識別情報が自由表面であれば、処理をステップS27へ進める。そうでなければ、処理をステップS25に進める。ここで、粒子(i)に対して所定範囲内にある粒子(j)が自由表面であれば、粒子(i)は、自由表面の近くに存在しており、以後の自由表面設定ステップにおいて、粒子(i)自身が自由表面粒子である可能性が高い。
【００３２】
ステップS25では、粒子(i)に対して、全ての粒子(j)に対するステップS24の処理が終了したかどうかを判定する。全ての粒子(j)に対する処理が完了していなければ、処理をステップS24に戻す。全ての粒子(j)に対する処理が完了していれば、処理をステップS26に進める。
【００３３】
ステップS26では、粒子(i)に対する近自由表面識別情報として「十分内部」という情報を設定する処理（どこかのメモリに記憶しておくこともできる）を行い、処理をステップS28へ進める。ステップS27では、粒子(i)に対する近自由表面識別情報として「近自由表面」という情報を設定する処理を行い、処理をステップS28へ進める。ステップS28では、全ての粒子(i)に対する近自由表面識別情報の設定処理が終了したかどうかを判定する。
【００３４】
全ての粒子(i)に対する処理が完了していなければ、処理をステップS22に戻す。全ての粒子(i)に対する処理が完了していれば、近自由表面設定の全ての処理S2を終了し、処理をステップS3（図３）へ進める。
【００３５】
ステップS3からS6では、各粒子に対して外力や粘性により及ぼされる力を計算し、それに基づいて流速（仮速度）及び位置（仮位置）が計算される(第１段階の計算)。第１段階の計算結果は、第２段階の計算に使用するために、どこかのメモリに記憶しておくことができる。この段階で得られる流速や位置は仮のものであり、以降の処理である第２段階の計算によって修正される。なお、これらのステップでの処理については、本願発明の出願時点で既に公知であるため、図２で紹介しておいた程度の説明に留めている。
【００３６】
ステップS7では、「粒子数密度」が計算される。粒子(i)の位置における粒子数密度ｎｉを求めるには、粒子(i)に対する他の粒子(j)の重み関数ｗの総和とされる。なお、これらのステップでの処理については、本願発明の出願時点で既に公知であるため、図２で紹介しておいた程度の説明に留めている。
【００３７】
ステップS8では、自由表面設定部１６が、自由表面かどうかの判定を行って自由表面識別情報の設定を行う。この処理は、各粒子に対して、近自由表面識別情報と粒子数密度の情報を用いて、自由表面判定を行う。前述のように、MPS法では、第１段階の計算によって仮の位置に移動した粒子の分布は一般的にはある程度のばらつきを持って配置されることになってしまう。このため、この段階では、十分内部にあるはずの粒子についても、その粒子数密度は低下した状態になってしまうことも少なくない。全ての粒子に対して、従来の粒子数密度に基づく閾値判定を適用してしまっては、このような粒子に対しても自由表面であるという判定結果を与えることになってしまう可能性があり、それがPoisson方程式の解の擾乱を招いていた。
【００３８】
これに対して、ここでは、近自由表面識別情報によって、第１段階の計算結果よりも以前のタイムステップにおいて、自由表面粒子であったか、もしくは、自由表面付近にあった粒子に対してだけ、粒子数密度に対する閾値評価を適用し、自由表面かどうかの判定を行う。以前のタイムステップが、直前のタイムステップに近いものであればあるほど流体の挙動の変化も少ないであろうから、より効果的であることが予測されるが、これに限られるものではない。これによって、前述の、十分内部にあるはずの粒子が自由表面であるとしてしまう誤判定を原理的に排除することができる。このステップS8のより詳細は、次の図５に示されたフローチャートによって示される処理を実行する。
【００３９】
ステップS81では、ある粒子(i)または全ての粒子について、位置x, y、粒子数密度n及び近自由表面識別情報を（メモリを参照して読み出すことで）取得する。ステップS82では、粒子(i)に対する近自由表面識別情報を取得する。ステップS83では、粒子(i)に対する自由表面識別情報が、十分内部であればステップS87へ、近自由表面であればステップS84へ処理を進める。ステップS84では、粒子(i)の粒子数密度nを取得する。ステップS85では、粒子(i)の粒子数密度が、所定の閾値よりも低ければ処理をステップS87へ進める。粒子(i)の粒子数密度が、所定の閾値よりも低くなければ処理をステップS86へ進める。
【００４０】
ステップS86では、粒子(i)に対する自由表面識別情報として「自由表面」という情報を設定する処理（どこかのメモリに記憶しておくこともできる）を行い、処理をステップS88へ進める。ステップS87では、粒子(i)に対する自由表面識別情報として「内部」という情報を設定する処理を行い、処理をステップS88へ進める。ステップS88では、全ての粒子(
i)に対する自由表面識別情報の設定処理が終了したかどうかを判定する。全ての粒子(i)に対する処理が完了していなければ、処理をステップS82に戻す。全ての粒子(i)に対する処理が完了していれば、自由表面設定の全ての処理S8を終了し、処理をステップS9（図３）へ進める。
【００４１】
非圧縮性流体の場合、各粒子の粒子数密度が一定となるようにしなければならない。しかし、上記第１段階の計算で得られた仮の位置にある粒子の分布において求められる粒子数密度には、一般にはばらつきが生じてしまう。ステップS9からS11では、このばらつきを均質な状態に修正する為に、Poisson方程式の解を用いて圧力項を決定し、この圧力勾配によって生じる力に基づいて流速と位置を修正する処理を行う。なお、これらのステップでの処理については、本願発明の出願時点で既に公知であるため、図２で紹介しておいた程度の説明に留めている。ステップS12では、解析を進めるべくタイムステップのインクリメント（＋ｄｔ）を行う。
【００４２】
ステップS13では、解析時間が計算終了時刻であるか否かを判定する。判定の結果、解析終了時刻である場合には処理を終了し、そうでない場合には、ステップS2の処理に戻す。
【００４３】
本願発明は、上記実施形態に限定されること無く適時設計変更可能である。例えば、フローチャートにおけるS2の実行のタイミングは、当業者であれば適時変更が可能である。また、本願発明の自由表面判定方法は、MPS法に限らず、流体の自由表面判定を要する他の手法(例えば、ISPH法等)に対しても適用することができる。
【００４４】
更に、近自由表面識別情報(S2)と合わせて用いる自由表面判定方法(S8)についても、本願発明の技術的思想の範囲内において、実施形態で用いた粒子数密度に基づく判定部分(S84S85)に替えて、他の判定部分を採用することも可能である。
【実施例】
【００４５】
図６は、２次元水柱崩壊という流体の挙動をシミュレーションする概要を示す図である。
【００４６】
概要としては、左端に置かれた水柱に重力ｇ（外力ｆ）がかかっており、その力が圧力に変換されて、下部から右に流れる現象をシミュレーションしている。Time 0, Time 1, Time 2, Time 3 は、時間が経過している流体の挙動の各状態を示している。
【００４７】
図７は、（ａ）が従来手法を用いた場合の自由表面粒子の分布と圧力の分布を示すシミュレーションの結果を示す図であり、（ｂ）が本願発明を用いた場合の自由表面粒子の分布と圧力の分布を示すシミュレーションの結果を示す図である。
【００４８】
図７（ａ）の従来手法においては、Time 2 と、Time 3とにおいて、自由表面粒子であると誤判定された内部粒子や、誤った表面設定条件による圧力場の空間的振動が現れてしまっている。これと比較して、図７（ｂ）の本願発明においては、それらが現れていないという効果を確認することができる。
【００４９】
図８は、本願発明における、次のタイムステップでの自由表面粒子の候補をどのあたりの部分から選択するかを示す模式図である。
【００５０】
従来技術においては、全ての粒子に対して、その周りにある粒子の情報から、粒子数密度を定義し、その値が大きければ流体の内部、小さければ自由表面、としていた。一方で、本願発明では、自由表面近くの粒子のインデックスをもっておいて、その中からしか、次のタイムステップでの自由表面粒子を選択しないようにしている、ことが特徴である。
【図１】
0005576309.tif 000005
【図２Ａ】
0005576309.tif 000006
【図２Ｂ】
0005576309.tif 000007
【図３】
0005576309.tif 000008
【図４】
0005576309.tif 000009
【図５】
0005576309.tif 000010
【図６】
0005576309.tif 000011
【図７】
0005576309.tif 000012
【図８】
0005576309.tif 000013