計算解析の対象となる連続体が複数の粒子により自由表面を伴って集合的にモデリングされるところの粒子法の計算処理の過程において、ある粒子が自由表面粒子に属するかどうかを判定する方法であって、計算処理の過程における以前のタイムステップにおけるメモリに記憶されている自由表面識別情報を参照して、ある粒子を基準にした所定の範囲内に、自由表面に属していた粒子が存在しているかどうかを判定するステップと、判定するステップにおいて、自由表面に属していた粒子が存在していると判定されたら近自由表面であるとし、存在していないと判定されたら十分内部であるとして、以後のタイムステップの参照のために近自由表面識別情報としてメモリに記憶するステップとを有する、方法。計算解析の対象となる連続体の複数の粒子の全てについて、請求項1に記載の全てのステップを繰り返すことを特徴とする、請求項1に記載の方法。さらに、メモリに記憶されている近自由表面識別情報を参照して、ある粒子が十分内部であると判定された場合には、内部であるとして、以後のタイムステップの参照のために自由表面識別情報としてメモリに記憶するステップと、を有する、請求項1に記載の方法。さらに、メモリに記憶されている近自由表面識別情報を参照して、ある粒子が近自由表面であると判定された場合には、以前のタイムステップにおける粒子数密度と所定の粒子数密度の閾値とを比較するステップと、比較するステップにおいて、所定の粒子数密度の閾値よりも低い場合には、自由表面であるとして、以後のタイムステップの参照のために自由表面識別情報としてメモリに記憶するステップと、比較するステップにおいて、所定の粒子数密度の閾値よりも低くない場合には、内部であるとして、以後のタイムステップの参照のために自由表面識別情報としてメモリに記憶するステップとを有する、請求項3に記載の方法。計算解析の対象となる連続体が複数の粒子により自由表面を伴って集合的にモデリングされるところの粒子法の計算処理の過程において、ある粒子が自由表面粒子に属するかどうかを判定する装置であって、計算処理の過程における以前のタイムステップにおけるメモリに記憶されている自由表面識別情報を参照して、ある粒子を基準にした所定の範囲内に、自由表面に属していた粒子が存在しているかどうかを判定する手段と、判定する手段において、自由表面に属していた粒子が存在していると判定されたら近自由表面であるとし、存在していないと判定されたら十分内部であるとして、以後のタイムステップの参照のために近自由表面識別情報としてメモリに記憶する手段とを有する、装置。計算解析の対象となる連続体の複数の粒子の全てについて、請求項5に記載の各手段の処理を繰り返すことを特徴とする、請求項5に記載の装置。さらに、メモリに記憶されている近自由表面識別情報を参照して、ある粒子が十分内部であると判定された場合には、内部であるとして、以後のタイムステップの参照のために自由表面識別情報としてメモリに記憶する手段と、を有する、請求項5に記載の装置。さらに、メモリに記憶されている近自由表面識別情報を参照して、ある粒子が近自由表面であると判定された場合には、以前のタイムステップにおける粒子数密度と所定の粒子数密度の閾値とを比較する手段と、比較する手段において、所定の粒子数密度の閾値よりも低い場合には、自由表面であるとして、以後のタイムステップの参照のために自由表面識別情報としてメモリに記憶する手段と、比較する手段において、所定の粒子数密度の閾値よりも低くない場合には、内部であるとして、以後のタイムステップの参照のために自由表面識別情報としてメモリに記憶する手段とを有する、請求項7に記載の装置。計算解析の対象となる連続体が複数の粒子により自由表面を伴って集合的にモデリングされるところの粒子法の計算処理の過程において、ある粒子が自由表面粒子に属するかどうかを判定するプログラムであって、計算処理の過程における以前のタイムステップにおけるメモリに記憶されている自由表面識別情報を参照して、ある粒子を基準にした所定の範囲内に、自由表面に属していた粒子が存在しているかどうかを判定するステップと、判定するステップにおいて、自由表面に属していた粒子が存在していると判定されたら近自由表面であるとし、存在していないと判定されたら十分内部であるとして、以後のタイムステップの参照のために近自由表面識別情報としてメモリに記憶するステップと有しており、これらのステップをコンピュータに実行させる、プログラム。本願発明の目的は、以上のような背景に鑑みてなされたものであり、流体の自由表面に位置する粒子を正確に判定することにある。本願発明の方法では、計算解析の対象となる連続体が複数の粒子により自由表面を伴って集合的にモデリングされるところの粒子法PMの計算処理の過程において、ある粒子が自由表面粒子に属するかどうかを判定する方法であって、計算処理の過程における以前のタイムステップにおけるメモリに記憶されている自由表面識別情報を参照して、ある粒子を基準にした所定の範囲内に、自由表面に属していた粒子が存在しているかどうかを判定するステップと、判定するステップにおいて、自由表面に属していた粒子が存在していると判定されたら近自由表面であるとし、存在していないと判定されたら十分内部であるとして、以後のタイムステップの参照のために近自由表面識別情報としてメモリに記憶するステップとを有するものが開示される。各ステップの処理を実行する手段、または、コンピュータにこれらのステップを実行させるプログラムという態様でも実現される。計算解析の対象となる連続体の複数の粒子の全てについて、ステップを繰り返す。さらに、メモリに記憶されている近自由表面識別情報を参照して、ある粒子が十分内部であると判定された場合には、内部であるとして、以後のタイムステップの参照のために自由表面識別情報としてメモリに記憶するステップとを有するものが開示される。さらに、メモリに記憶されている近自由表面識別情報を参照して、ある粒子が近自由表面であると判定された場合には、以前のタイムステップにおける粒子数密度が所定の粒子数密度の閾値と比較するステップと、比較するステップにおいて、所定の粒子密度の閾値よりも低い場合には、自由表面であるとして、以後のタイムステップの参照のために自由表面識別情報としてメモリに記憶するステップと、比較するステップにおいて、所定の粒子密度の閾値よりも低くない場合には、内部であるとして、以後のタイムステップの参照のために自由表面識別情報としてメモリに記憶するステップとを有するものが開示される。より正確に自由表面粒子の判定を行うことが可能になる。従来手法と本願発明との比較をした具体的なシミュレーション結果を、実施例において示すことにする。図1は、本願発明である自由表面粒子判定装置を含む形で構成された、粒子法を用いて流体の挙動を解析する流体解析装置1の概略図である。流体解析装置1は、RAM2、初期設定データベース3、CPU4、を主な構成要素として備える。RAM2は、CPU4が実行する流体解析装置1の起動プログラムや流体解析プログラム等のコンピュータプログラム及び各種データを一時的に格納するワーキングエリアを提供する。初期設定DB3は、粒子法の一種であるMPS法を例にとるとMPS法によって流体解析を実行するために必要な各種パラメータを格納する。具体的には、メモリとして、標準的なMPS法で必要となる、粒子の位置、粒子の速度、粒子の密度、動粘性係数等の物性に関するデータや、重力加速度等の外力に関するデータ、初期配置における粒子間距離に関するデータ、粒子数密度を求めるための基準半径等の計算データが格納される。更に、本願発明の実施形態ではこれに加えて、流体を構成する各粒子に対して、自由表面に配置されているかどうかを示す自由表面識別情報が記憶される。このように、ハードウエア資源としての構成、ソフトウエア資源としての構成、または、ハードウエア資源とソフトウエア資源とが協働する構成に成っている流体解析装置1は、流体解析処理を実行する。ソフトウエア資源は、外部からダウンロードするようにすることもでき、コンピュータに実行させるためのプログラムとしても提供することができる。図2は、図2Aが、粒子法の一例であるMPS法の流体解析方法の概要を説明する図である。非圧縮性流体の運動方程式を、第1段階と第2段階とに分けて、処理を進めていくものである。図2Bは、境界条件設定の概要を説明する図である。この図2Aおよび図2Bで説明している処理自体は、本願発明の出願時点で既に公知であるので、ごく簡単な説明に留めている。以下、図3に示すフローチャートを参照して、この流体解析装置1の動作について詳しく説明する。図3に示すフローチャートでは、まず、流体解析処理はステップS1の処理へ進む。ステップS1では、CPU4が、初期設定DB3からRAM2内に解析対象となる非圧縮性流体に関する初期の計算データを読み出して取得する。ステップS2では、近自由表面設定部11が、近自由表面設定の処理を行う。この処理は、ある粒子または全ての粒子に対して、この時点でその粒子が自由表面であるか、もしくはその粒子を基準にした所定の範囲内に自由表面粒子が存在すると判定されれば近自由表面という識別情報を設定し、自由表面粒子が所定の範囲内に存在しないと判定された粒子に対しては、十分内部という識別情報を設定する。これは内部という識別情報との区別においての命名にすぎず、必ずしも実際の計算解析の対象における空間において十分な内部に粒子があることを意味しているものではない。この識別情報は、近自由表面識別情報としてRAM2他のどこかのメモリでもよいに一時的に記憶しておくことができる。この識別情報は、ある粒子または全ての粒子に対して設定されるものであるため、流速の情報や位置の情報が通常取り扱われるのと同様、配列のデータ保持形式で、計算処理のタイムステップを進めていく過程において、これを流体の挙動をシミュレートするために必要な時点まで保持しておけばよい。このステップS2の処理は、以後の処理である自由表面判定を行う際に、その判定を適用する粒子を近自由表面に属する粒子に制限し、十分内部にある粒子が誤って自由表面であると判定されることを防ぐための処理である。流体解析のタイムステップの発展tdtは、計算が数値的に安定に実行できるように、1タイムステップdtの刻み幅を制限している。これによれば、流体全体の形状が大幅に変更されることがないので、自由表面から十分に離れ、流体の内部にある粒子が、高々1タイムステップdtの発展によって、自由表面粒子となることは考えにくい。従って、この段階で十分内部と判定された粒子に対しては、以後の境界条件設定ステップにおいて、自由表面であると判定することを避けることが望ましい。このステップS2のより詳細は、次の図4に示されたフローチャートによって示される処理として実行する。ステップS21では、ある粒子または全ての粒子について、位置xyと自由表面識別情報をメモリを参照して読み出すことで取得する。ステップS22では、ある粒子に対する自由表面識別情報を取得する。ステップS23では、ある粒子に対する自由表面識別情報が、自由表面であればステップS27へ、自由表面でなければステップS24へ処理を進める。ここで、粒子が自由表面であると識別される場合、以後の自由表面設定を行う際、再び自由表面粒子となる可能性が高い。ステップS24では、まず、粒子と、その他の粒子それぞれの位置xyから、それらの距離を求める。粒子と粒子との距離とは、次式2で定義されるものである。ここで、xiyiは粒子のxy座標、xjyjは粒子のxy座標である。その距離が所定範囲Rメモリに記憶されているところの所定の閾値内にあると判定されれば、即ち、次式3が成立すれば、その粒子の自由表面識別情報を参照する。そうでなければ、処理をステップS25に進める。更に、粒子の自由表面識別情報が自由表面であれば、処理をステップS27へ進める。そうでなければ、処理をステップS25に進める。ここで、粒子に対して所定範囲内にある粒子が自由表面であれば、粒子は、自由表面の近くに存在しており、以後の自由表面設定ステップにおいて、粒子自身が自由表面粒子である可能性が高い。ステップS25では、粒子に対して、全ての粒子に対するステップS24の処理が終了したかどうかを判定する。全ての粒子に対する処理が完了していなければ、処理をステップS24に戻す。全ての粒子に対する処理が完了していれば、処理をステップS26に進める。ステップS26では、粒子に対する近自由表面識別情報として十分内部という情報を設定する処理どこかのメモリに記憶しておくこともできるを行い、処理をステップS28へ進める。ステップS27では、粒子に対する近自由表面識別情報として近自由表面という情報を設定する処理を行い、処理をステップS28へ進める。ステップS28では、全ての粒子に対する近自由表面識別情報の設定処理が終了したかどうかを判定する。全ての粒子に対する処理が完了していなければ、処理をステップS22に戻す。全ての粒子に対する処理が完了していれば、近自由表面設定の全ての処理S2を終了し、処理をステップS3図3へ進める。ステップS3からS6では、各粒子に対して外力や粘性により及ぼされる力を計算し、それに基づいて流速仮速度及び位置仮位置が計算される。第1段階の計算結果は、第2段階の計算に使用するために、どこかのメモリに記憶しておくことができる。この段階で得られる流速や位置は仮のものであり、以降の処理である第2段階の計算によって修正される。なお、これらのステップでの処理については、本願発明の出願時点で既に公知であるため、図2で紹介しておいた程度の説明に留めている。ステップS7では、粒子数密度が計算される。粒子の位置における粒子数密度niを求めるには、粒子に対する他の粒子の重み関数wの総和とされる。なお、これらのステップでの処理については、本願発明の出願時点で既に公知であるため、図2で紹介しておいた程度の説明に留めている。ステップS8では、自由表面設定部16が、自由表面かどうかの判定を行って自由表面識別情報の設定を行う。この処理は、各粒子に対して、近自由表面識別情報と粒子数密度の情報を用いて、自由表面判定を行う。前述のように、MPS法では、第1段階の計算によって仮の位置に移動した粒子の分布は一般的にはある程度のばらつきを持って配置されることになってしまう。このため、この段階では、十分内部にあるはずの粒子についても、その粒子数密度は低下した状態になってしまうことも少なくない。全ての粒子に対して、従来の粒子数密度に基づく閾値判定を適用してしまっては、このような粒子に対しても自由表面であるという判定結果を与えることになってしまう可能性があり、それがPoisson方程式の解の擾乱を招いていた。これに対して、ここでは、近自由表面識別情報によって、第1段階の計算結果よりも以前のタイムステップにおいて、自由表面粒子であったか、もしくは、自由表面付近にあった粒子に対してだけ、粒子数密度に対する閾値評価を適用し、自由表面かどうかの判定を行う。以前のタイムステップが、直前のタイムステップに近いものであればあるほど流体の挙動の変化も少ないであろうから、より効果的であることが予測されるが、これに限られるものではない。これによって、前述の、十分内部にあるはずの粒子が自由表面であるとしてしまう誤判定を原理的に排除することができる。このステップS8のより詳細は、次の図5に示されたフローチャートによって示される処理を実行する。ステップS81では、ある粒子または全ての粒子について、位置xy、粒子数密度n及び近自由表面識別情報をメモリを参照して読み出すことで取得する。ステップS82では、粒子に対する近自由表面識別情報を取得する。ステップS83では、粒子に対する自由表面識別情報が、十分内部であればステップS87へ、近自由表面であればステップS84へ処理を進める。ステップS84では、粒子の粒子数密度nを取得する。ステップS85では、粒子の粒子数密度が、所定の閾値よりも低ければ処理をステップS87へ進める。粒子の粒子数密度が、所定の閾値よりも低くなければ処理をステップS86へ進める。ステップS86では、粒子に対する自由表面識別情報として自由表面という情報を設定する処理どこかのメモリに記憶しておくこともできるを行い、処理をステップS88へ進める。ステップS87では、粒子に対する自由表面識別情報として内部という情報を設定する処理を行い、処理をステップS88へ進める。ステップS88では、全ての粒子に対する自由表面識別情報の設定処理が終了したかどうかを判定する。全ての粒子に対する処理が完了していなければ、処理をステップS82に戻す。全ての粒子に対する処理が完了していれば、自由表面設定の全ての処理S8を終了し、処理をステップS9図3へ進める。非圧縮性流体の場合、各粒子の粒子数密度が一定となるようにしなければならない。しかし、上記第1段階の計算で得られた仮の位置にある粒子の分布において求められる粒子数密度には、一般にはばらつきが生じてしまう。ステップS9からS11では、このばらつきを均質な状態に修正する為に、Poisson方程式の解を用いて圧力項を決定し、この圧力勾配によって生じる力に基づいて流速と位置を修正する処理を行う。なお、これらのステップでの処理については、本願発明の出願時点で既に公知であるため、図2で紹介しておいた程度の説明に留めている。ステップS12では、解析を進めるべくタイムステップのインクリメントdtを行う。ステップS13では、解析時間が計算終了時刻であるか否かを判定する。判定の結果、解析終了時刻である場合には処理を終了し、そうでない場合には、ステップS2の処理に戻す。本願発明は、上記実施形態に限定されること無く適時設計変更可能である。例えば、フローチャートにおけるS2の実行のタイミングは、当業者であれば適時変更が可能である。また、本願発明の自由表面判定方法は、MPS法に限らず、流体の自由表面判定を要する他の手法に対しても適用することができる。更に、近自由表面識別情報と合わせて用いる自由表面判定方法についても、本願発明の技術的思想の範囲内において、実施形態で用いた粒子数密度に基づく判定部分に替えて、他の判定部分を採用することも可能である。図6は、2次元水柱崩壊という流体の挙動をシミュレーションする概要を示す図である。概要としては、左端に置かれた水柱に重力g外力fがかかっており、その力が圧力に変換されて、下部から右に流れる現象をシミュレーションしている。Time0Time1Time2Time3は、時間が経過している流体の挙動の各状態を示している。図7は、aが従来手法を用いた場合の自由表面粒子の分布と圧力の分布を示すシミュレーションの結果を示す図であり、bが本願発明を用いた場合の自由表面粒子の分布と圧力の分布を示すシミュレーションの結果を示す図である。図7aの従来手法においては、Time2と、Time3とにおいて、自由表面粒子であると誤判定された内部粒子や、誤った表面設定条件による圧力場の空間的振動が現れてしまっている。これと比較して、図7bの本願発明においては、それらが現れていないという効果を確認することができる。図8は、本願発明における、次のタイムステップでの自由表面粒子の候補をどのあたりの部分から選択するかを示す模式図である。従来技術においては、全ての粒子に対して、その周りにある粒子の情報から、粒子数密度を定義し、その値が大きければ流体の内部、小さければ自由表面、としていた。一方で、本願発明では、自由表面近くの粒子のインデックスをもっておいて、その中からしか、次のタイムステップでの自由表面粒子を選択しないようにしている、ことが特徴である。