"1つ以上の変形体の運動の解析を行う解析装置であって、その中心が前記変形体を分割する四面体の頂点を構成する複数の仮想的な粒子が設けられた当該変形体のモデルを生成するモデル生成手段と、前記モデル生成手段によって生成されたモデルにおいて、前記粒子のうち前記変形体の表面に設けられると共に互いに隣接しない粒子の間に生じる第1の力を算出する第1の力算出手段と、前記モデル生成手段によって生成されたモデルにおいて、前記変形体を分割する四面体に生じる応力を算出し、当該四面体を分割する面を求めることによって、算出した応力から当該四面体を構成する各粒子に生じる第2の力を算出する第2の力算出手段と、各粒子について、前記第1の力算出手段によって算出された第1の力と、前記第2の力算出手段によって算出された第2の力との和を取ることで各粒子に生じる力を算出して、算出した力に基づいて運動方程式を解くことで各粒子の位置を更新する粒子位置更新手段と、を備える解析装置。前記モデル生成手段は、変形体の数に応じて配置された前記粒子と予め設定した位置関係にあると共に前記四面体の頂点を構成する粒子を設定された終了条件を満たすまで更に配置することで、前記変形体のモデルを生成することを特徴とする請求項1に記載の解析装置。前記第1の力算出手段は、第1の力を、各粒子の位置、並びに予め設定される粒子のパラメータに基づいて算出することを特徴とする請求項1又は2に記載の解析装置。予め設定される粒子のパラメータは、粒子の半径、質量及びバネ係数であることを特徴とする請求項3に記載の解析装置。前記第2の力算出手段は、前記変形体を分割する四面体に生じる応力を、当該四面体を構成する各粒子の位置の時間変化又は速度に基づいて算出することを特徴とする請求項1~4の何れか一項に記載の解析装置。1つ以上の変形体の運動の解析を行う解析装置による解析方法であって、その中心が前記変形体を分割する四面体の頂点を構成する複数の仮想的な粒子が設けられた当該変形体のモデルを生成するモデル生成ステップと、前記モデル生成ステップにおいて生成されたモデルにおいて、前記粒子のうち前記変形体の表面に設けられると共に互いに隣接しない粒子の間に生じる第1の力を算出する第1の力算出ステップと、前記モデル生成ステップにおいて生成されたモデルにおいて、前記変形体を分割する四面体に生じる応力を算出し、当該四面体を分割する面を求めることによって、算出した応力から当該四面体を構成する各粒子に生じる第2の力を算出する第2の力算出ステップと、各粒子について、前記第1の力算出ステップにおいて算出された第1の力と、前記第2の力算出ステップにおいて算出された第2の力との和を取ることで各粒子に生じる力を算出して、算出した力に基づいて運動方程式を解くことで各粒子の位置を更新する粒子位置更新ステップと、を含む解析方法。コンピュータに、1つ以上の変形体の運動の解析を行わせる解析プログラムであって、前記コンピュータを、その中心が前記変形体を分割する四面体の頂点を構成する複数の仮想的な粒子が設けられた当該変形体のモデルを生成するモデル生成手段と、前記モデル生成手段によって生成されたモデルにおいて、前記粒子のうち前記変形体の表面に設けられると共に互いに隣接しない粒子の間に生じる第1の力を算出する第1の力算出手段と、前記モデル生成手段によって生成されたモデルにおいて、前記変形体を分割する四面体に生じる応力を算出し、当該四面体を分割する面を求めることによって、算出した応力から当該四面体を構成する各粒子に生じる第2の力を算出する第2の力算出手段と、各粒子について、前記第1の力算出手段によって算出された第1の力と、前記第2の力算出手段によって算出された第2の力との和を取ることで各粒子に生じる力を算出して、算出した力に基づいて運動方程式を解くことで各粒子の位置を更新する粒子位置更新手段と、して機能させる解析プログラム。コンピュータに、1つ以上の変形体の運動の解析を行わせる解析プログラムを記録したコンピュータ読み取り可能な記録媒体であって、前記解析プログラムは、前記コンピュータを、その中心が前記変形体を分割する四面体の頂点を構成する複数の仮想的な粒子が設けられた当該変形体のモデルを生成するモデル生成手段と、前記モデル生成手段によって生成されたモデルにおいて、前記粒子のうち前記変形体の表面に設けられると共に互いに隣接しない粒子の間に生じる第1の力を算出する第1の力算出手段と、前記モデル生成手段によって生成されたモデルにおいて、前記変形体を分割する四面体に生じる応力を算出し、当該四面体を分割する面を求めることによって、算出した応力から当該四面体を構成する各粒子に生じる第2の力を算出する第2の力算出手段と、各粒子について、前記第1の力算出手段によって算出された第1の力と、前記第2の力算出手段によって算出された第2の力との和を取ることで各粒子に生じる力を算出して、算出した力に基づいて運動方程式を解くことで各粒子の位置を更新する粒子位置更新手段と、して機能させる記録媒体。ところで、砂利やワイヤーやケーブルは、個々の石や鋼線が多数集まったものであるが、個々の石や鋼線の変形に関する性質が全体の挙動に大きく影響する。しかし、その集合体としての性質や、個々の石や鋼線の性質だけからは決まらない。例えば、砂利ならばその敷きならし方や盛り上げた形状に、ワイヤーやケーブルならば鋼線の撚り方によって全体の性質は大きく変化する。このような物体に対しては、個々の石や鋼線等の要素間の衝突や接触と、個々の要素自体の変形との両方を考慮しなければ、適切な解析、例えば要素内の弾性波伝播や破壊に伴う運動を正確にシミュレーションすることができない。石や鋼線等の要素は接触圧が大きく摩擦力が大きい。即ち、上記の接触や衝突としては、具体的には、摩擦力やすべり特性を考慮する必要がある。しかしながら、上述した不連続変形法では、変形を考慮しているものの個々の要素についての変形に係る計算を行った上で衝突や接触の計算を行うため、個々の石や鋼線等の要素が多くなると、計算負荷が非常に大きくなるため計算時間が大きくなってしまうという問題点があった。そのため、このような物体に対しては経験的な評価しか行えないのが原状である。また、DEMでは、複数個の球形粒子を連結して一つの変形要素をモデル化する方法も有るが、正しい物性を再現する構成則が存在しないために要素自体の変形や要素内部の応力解析を正しく行うことができない。更に、QDEMでは、単一要素の変形や要素内の応力を構成則に基づいて解析することは可能であるが、複数の要素間の接触や衝突判定方法及び相互作用力の算出機能を有していないため、ここで目的とする多数の変形する要素の運動を解析することは不可能である。本発明は、上記の問題点を鑑みてなされたものであり、対象とする要素である変形体物体が多い場合でも、変形体間の衝突や接触と変形体自体の変形との両方を考慮して高速に解析を行うことができる解析装置、解析方法、解析プログラム及び記録媒体を提供することを目的とする。上記の目的を達成するために、本発明に係る解析装置は、1つ以上の変形体の運動の解析を行う解析装置であって、その中心が変形体を分割する四面体の頂点を構成する複数の仮想的な粒子が設けられた当該変形体のモデルを生成するモデル生成手段と、モデル生成手段によって生成されたモデルにおいて、粒子のうち変形体の表面に設けられると共に互いに隣接しない粒子の間に生じる第1の力を算出する第1の力算出手段と、モデル生成手段によって生成されたモデルにおいて、変形体を分割する四面体に生じる応力を算出し、当該四面体を分割する面を求めることによって、算出した応力から当該四面体を構成する各粒子に生じる第2の力を算出する第2の力算出手段と、各粒子について、第1の力算出手段によって算出された第1の力と、第2の力算出手段によって算出された第2の力との和を取ることで各粒子に生じる力を算出して、算出した力に基づいて運動方程式を解くことで各粒子の位置を更新する粒子位置更新手段と、を備える。本発明に係る解析装置では、変形体間の衝突や接触に係る第1の力と変形体の変形に係る第2の力とがそれぞれ算出されて、それらの和の力に基づいて各粒子の位置が更新される。また、第1の力の算出と第2の力の算出との何れにおいても、変形体に設けられた仮想的な粒子が用いられる。従って、本発明に係る解析装置では、第1の力の算出と第2の力の算出との両方が、上記の仮想的な粒子という共通の項目に対する計算にて行うことができる。これにより、本発明に係る解析装置によれば、変形体が多い場合でも、変形体間の衝突や接触と変形体自体の変形との両方を考慮して高速に解析を行うことができる。モデル生成手段は、変形体の数に応じて配置された粒子と予め設定した位置関係にあると共に四面体の頂点を構成する粒子を設定された終了条件を満たすまで更に配置することで、変形体のモデルを生成することとしてもよい。この構成によれば、多数の変形体に対して規則的で品質のよい四面体が設けられたモデルを効率的に生成することができる。第1の力算出手段は、第1の力を、各粒子の位置、並びに予め設定される粒子のパラメータに基づいて算出することとしてもよい。予め設定される粒子のパラメータは、粒子の半径、質量及びバネ係数であることとしてもよい。これらの構成によれば、適切かつ確実に第1の力を算出することができる。前記第2の力算出手段は、前記変形体を分割する四面体に生じる応力を、当該四面体を構成する各粒子の位置の時間変化又は速度に基づいて算出することとしてもよい。この構成によれば、適切かつ確実に第2の力を算出することができる。ところで、本発明は、上記のように解析装置の発明として記述できる他に、以下のように解析方法及び解析プログラムの発明としても記述することができる。これはカテゴリ等が異なるだけで、実質的に同一の発明であり、同様の作用及び効果を奏する。本発明に係る解析方法は、1つ以上の変形体の運動の解析を行う解析装置による解析方法であって、その中心が変形体を分割する四面体の頂点を構成する複数の仮想的な粒子が設けられた当該変形体のモデルを生成するモデル生成ステップと、モデル生成ステップにおいて生成されたモデルにおいて、粒子のうち変形体の表面に設けられると共に互いに隣接しない粒子の間に生じる第1の力を算出する第1の力算出ステップと、モデル生成ステップにおいて生成されたモデルにおいて、変形体を分割する四面体に生じる応力を算出し、当該四面体を分割する面を求めることによって、算出した応力から当該四面体を構成する各粒子に生じる第2の力を算出する第2の力算出ステップと、各粒子について、第1の力算出ステップにおいて算出された第1の力と、第2の力算出ステップにおいて算出された第2の力との和を取ることで各粒子に生じる力を算出して、算出した力に基づいて運動方程式を解くことで各粒子の位置を更新する粒子位置更新ステップと、を含む。本発明に係る解析プログラムは、コンピュータに、1つ以上の変形体の運動の解析を行わせる解析プログラムであって、コンピュータを、その中心が変形体を分割する四面体の頂点を構成する複数の仮想的な粒子が設けられた当該変形体のモデルを生成するモデル生成手段と、モデル生成手段によって生成されたモデルにおいて、粒子のうち変形体の表面に設けられると共に互いに隣接しない粒子の間に生じる第1の力を算出する第1の力算出手段と、モデル生成手段によって生成されたモデルにおいて、変形体を分割する四面体に生じる応力を算出し、当該四面体を分割する面を求めることによって、算出した応力から当該四面体を構成する各粒子に生じる第2の力を算出する第2の力算出手段と、各粒子について、第1の力算出手段によって算出された第1の力と、第2の力算出手段によって算出された第2の力との和を取ることで各粒子に生じる力を算出して、算出した力に基づいて運動方程式を解くことで各粒子の位置を更新する粒子位置更新手段と、して機能させる。本発明に係る記録媒体は、コンピュータに、1つ以上の変形体の運動の解析を行わせる解析プログラムを記録したコンピュータ読み取り可能な記録媒体であって、解析プログラムは、コンピュータを、その中心が変形体を分割する四面体の頂点を構成する複数の仮想的な粒子が設けられた当該変形体のモデルを生成するモデル生成手段と、モデル生成手段によって生成されたモデルにおいて、粒子のうち変形体の表面に設けられると共に互いに隣接しない粒子の間に生じる第1の力を算出する第1の力算出手段と、モデル生成手段によって生成されたモデルにおいて、変形体を分割する四面体に生じる応力を算出し、当該四面体を分割する面を求めることによって、算出した応力から当該四面体を構成する各粒子に生じる第2の力を算出する第2の力算出手段と、各粒子について、第1の力算出手段によって算出された第1の力と、第2の力算出手段によって算出された第2の力との和を取ることで各粒子に生じる力を算出して、算出した力に基づいて運動方程式を解くことで各粒子の位置を更新する粒子位置更新手段と、して機能させる。本発明では、変形体間の衝突や接触に係る第1の力の算出と変形体の変形に係る第2の力の算出との両方が、仮想的な粒子という共通の項目に対する計算にて行うことができる。これにより、本発明によれば、変形体が多い場合でも、変形体間の衝突や接触と変形体自体の変形との両方を考慮して高速に解析を行うことができる。以下、図面と共に本発明に係る解析装置、解析方法及び解析プログラムの好適な実施形態について詳細に説明する。なお、図面の説明においては同一要素には同一符号を付し、重複する説明を省略する。図1に本実施形態に係る解析装置10を示す。解析装置10は、1つ以上の変形体の運動の解析シミュレーションを行う。変形体は、例えば、石や鋼線等の物体である。変形体は、固体弾性体あるいは粘弾性体でもよいし、あるいは粘性流体や粘弾性流体でもよい。本実施形態に係る解析装置10は、変形体をモデル化して、解析により時間時刻毎の変形体の位置例えば、三次元座標を算出する。その際、変形体間、あるいは1つの変形体の異なる部分の間の接触衝突を考慮する。また、変形体に加わる力によって、変形体が変形することも考慮する。解析装置10は、例えば、CPUCentralProcessingUnit、GPUGraphicsProcessingUnit、メモリ、ハードディスク、ディスプレイ等のハードウェアを備えるコンピュータとして構成される。これらの構成要素がプログラム等により動作することによって、後述する解析装置10としての機能が発揮される。図1に示すように解析装置10は、記憶部11と、モデル生成部12と、第1の力算出部13と、第2の力算出部14と、粒子位置更新部15と、出力部16とを備える。記憶部11は、解析装置10による解析に必要な情報を記憶する手段である。記憶部11は、解析に予めパラメータとして必要である情報については予め記憶している。なお、この情報については、例えば、コントロールファイル等として解析装置10の利用者によって予め解析装置10に入力されている。また、解析の過程で、算出される情報については必要に応じて後述する各構成要素から、記憶部11に入力される書き込まれる。記憶部11が記憶している情報は、必要に応じて各構成要素モデル生成部12、第1の力算出部13、第2の力算出部14、粒子位置更新部15及び出力部16から、参照される。具体的に、記憶部11がどのような情報が記憶するかは、その情報の利用、算出の時点で説明する。モデル生成部12は、解析の対象となる変形体のモデルを生成するモデル生成手段である。モデル生成部12は、変形体のモデルを、三次元形状の変形体が存在する領域に複数の仮想的な粒子を設けることによって生成する。当該粒子は、その中心が変形体を分割する四面体の頂点を構成する。粒子が位置する領域に上記の変形体があるものとして解析が行われる。モデル生成部12は、変形体の形状及び初期配置を示す情報を記憶部11から読み出して、当該情報に基づいて、有限要素法FEM:FiniteElementMethod等で利用される既存の四面体メッシュの自動生成方法によって、四面体に粒子を設けることとしてもよい。即ち、1変形体の表面に三角形ポリゴン網を生成する点を張り、2変形体内部空間に三角形ポリゴンと同程度の数密度を保った多数点を発生させ、3隣接4点間の最小外接球を数値的に見出すことで四面体を逐次決定する。このとき、各四面体の品質については、四面体の稜線辺の長さをより均一に近づけることで管理される。なお、ここで四面体の品質とは、例えば、四面体を構成する辺の長さや四面体間の大きさのばらつきを指すばらつきが小さい方が品質がよい。また手順として、より大まかで比較的均質な四面体メッシュを予め見つけた後に各四面体のメッシュを段階的に細分化することで、四面体品質を保ちながら微小四面体分割を行う手法を用いてもよい。上記の方法で生成した四面体の頂点を中心とした粒子を設けて、モデルを生成する。上記の一般的な四面体メッシュ生成方法を用いた方法以外でも、モデル生成部12は、以下の方法によって変形体のモデルを生成することとしてもよい。即ち、モデル生成部12は、変形体の数に応じて仮想的な粒子を配置し、当該粒子と予め設定した位置関係にあると共に四面体の頂点を構成する粒子を設定された終了条件を満たすまで更に配置する新たに配置していくことで、変形体のモデルを生成することとしてもよい。例えば、以下のようにモデルを生成する。上記の一般的な四面体メッシュ生成方法は、本来1つの物体や領域の四面体分割を目的として作られたもので、物体形状や領域形状を予め入力として与えなければならない。従って、多数の物体を対象としたり、更に個々の物体が不規則形状を有したりする場合には、物体群を自動生成しつつ個々の物体を多数の四面体に自動分割することには、本質的に不向きである。例えば、解析の対象を多数の砂利によって構成される石のような多数の不規則形状の変形体とする場合には、以下のような、変形体を自動発生させながら内部を四面体分割する方法を用いることが望ましい。そこで本手法では、変形体内部を四面体分割するのではなく、変形体内部の点から面心立方格子とその内部に包含される2種類の四面体を用いて、解析に必要とされる変形体の形状まで成長させることで、多数の物体に対して規則的で品質のよい四面体メッシュを自動発生させる。図2に面心立方格子の基本ユニットを示す。図中に示される0,0,0を原点とする直交座標座標系x-y-zを定義する。そこに原点を重心として一辺の長さが2である3枚の正方形をそれぞれxy平面z=0、yz平面x=0、zx平面上y=0に考える。すると3枚の正方形の12個の頂点1,1,0,-1,1,0,-1,-1,0,1,-1,0,0,1,1,0,-1,1,0,-1,-1,0,1,-1,1,0,1,-1,0,1,-1,0,-1,1,0,-1は、全て原点から√2の等距離に位置することが分かる。この12点と原点を合わせた合計13点が面心立方格子の基本ユニットの格子点である。面心立方格子の基本ユニットである13点の格子点の中から隣接格子点距離が全て√2である格子点4点を選んで四面体を作ると、図3に示すように一辺の長さが√2の正三角形4面からなる正四面体となる。yz平面の正方形の表側x>0を考えると、図3に示す四面体と合同な正四面体が4つできる。同様にyz平面の正方形の裏側x<0にも同じ四面体が4つ考えることができ、この基本ユニットの中には合計8つの正四面体が存在することが分かる。上述の方法でできるyz平面の表側x>0の領域の4つの正四面体の間には、図4aに示すように底面が1辺の長さ√2の正方形で高さが1の四角錐の空洞ができる。同様に、この四角錐の空洞はx<0側にもでき、更にzx平面のy>0側とy<0側、xy平面のz>0側とz<0側の合計6方向に生じる。図4aに見られる四角錐の空洞は、0,0,0を中心とする面心立方格子の全ての隣のユニットにも同様に生じる。例えば、2,0,0を中心とするユニットのyz平面のx<0側の4つの正四面体の間には、図4bに示す四角錐の空洞が生じる。図4aに示す空洞と、図4bに示す空洞とは、底面で同一正方形を共有し、かつ2つの四角錐の頂点を結ぶx軸上にある線分0,0,0-2,0,0は、この正方形と重心で垂直に交わる。すると、図4aに示す四角錐と、図4bに示す四角錐とを合体させてできる立体は、図5に示すように一辺の長さ√2の正三角形8枚からなる正八面体となる。一方、図3に示す正四面体について、頂点を原点0,0,0とし底面三角形を1,1,0-0,1,1-1,0,1と考えるとき、底面三角形において頂点と反対側の空間にも0,0,0のユニットには空洞が生じている。この空洞を埋め尽くすために、図3に示したxy平面z=0上の正方形に対してz方向の隣接ユニットである0,0,2を中心とするxy平面の正方形の頂点のうち点1,0,1と0,1,1に最も隣接している1,1,2との位置関係を考える。するとそこには、図6aに示すような4つの頂点1,1,2,0,1,1,1,0,1,1,1,0から構成される四面体領域ができる。同様に四面体領域は、線分1,1,0-1,1,2を囲んで中心が2,0,0,2,2,0,0,2,0の3つの隣接ユニットにも生じる。これら4つの隣接ユニットにおける線分1,1,0-1,1,2を共有する四面体領域を全て足し合わせると図6bに示すような八面体となることが分かる。この八面体表面を構成する三角形の一辺の長さは全て√2であることから、図6bに示す八面体は正八面体となり、同時に図5に示す正八面体と合同である。図6bに示す八面体と同様の正八面体は、図3で考えた8つの正四面体の全ての底面三角形について考えることができる。ここまで示されたように、図2に示した1つの面心立方格子基本ユニットの空間は、8つの正四面体と、隣接ユニットと共有される図5に示す正八面体6つと、図6bに示す正八面体8つとによって埋め尽くされていることが分かる。図5に示す正八面体と、図6bに示す正八面体とは前述の通り合同であったので、この14個の正八面体をそれぞれ四面体に分割すると、面心立方と四面体との関係が決定することになる。正八面体を四面体に分割する方法は一意ではないが、例えば、元の正八面体の頂点のみを用いて新たな頂点を増やさずに4つの合同な四面体に分割することとすればよい。その方法は、分割によって新たに追加される四面体の稜四面体を構成する三角形の辺の方向の違いによって3通り存在する。便宜上、図7aに示すように正八面体の重心を原点に置き、3組の向かい合う2つの頂点をそれぞれx軸、y軸、z軸上に乗るように配置する。すると、図7b~図7dに示されるように、正八面体を分割するために追加される辺をx軸にとる場合、y軸にとる場合、z軸にとる場合に対応して、3通りの分割方法が存在することが分かる。その結果、図7b~図7dに示すように何れも1つの正八面体を2面の直角二等辺三角形と2面の正三角形からなる合同な四面体4つに分割することができる。図7を用いて説明した3通りの分割方法は、正八面体の頂点の位置さえ決まれば、他に何も依存することなく決めることができる。そこで、例えば、乱数によって図7を用いて説明した3通りの分割方法から1つを等確率で与えて分割方法を決めることすればよい。これにより、大域的なメッシュ全体で3種類の分割方法が均等に出現することになり、四面体分割の等方性を確保することができる。1つの面心立方格子基本ユニットを囲む正八面体のうち、図5に示す1つの正八面体は2つのユニットに共有されている。このパターンが面心立方の立方体6方向の面で隣接するユニットに対して存在し、1つの正八面体は図7の何れかの方法で4つの四面体に分割される。従って、1つのユニットには、図5のパターンの四面体が4×6÷2=12個分が生成される。また、一方、図6bに示す1つの正八面体は8つの隣接ユニットに共有されている。このパターンが面心立方の立方体8方向の頂点で隣接するユニットに対して存在し、1つの正八面体は図7の何れかの方法で4つの四面体に分割される。従って、1つのユニットには、図6bのパターンの四面体が4×8÷8=4個分が生成される。以上までの説明で分かるように、図2で示した1つの面心立方格子基本ユニットの空間は、8つの正四面体と、直角二等辺三角形2枚及び正三角形2枚からなる四面体124=16個とに分割される。上記の面心立方格子を成長させることによって、変形体のモデルを生成する。まず、n個の変形体ここでは、マクロパーティクル、MPと呼ぶに対して、MPの種ここでは、シードパーティクル、SPと呼ぶを所定の三次元空間領域に発生させる。ここで、nは予め設定される変形体の数である。MPがランダムに配置されることを仮定するときは、三次元座標x,y,zを与えるn組の一様乱数を用い、MPが予め設定された場所に位置することを仮定するときは、当該位置に応じたn個の座標値例えば、MPの端点やMPの重心の位置を与える。次に、n個のSPに対して、面心立方格子の方向=MPの成長方向を与える。このとき、MPの成長方向にランダム性を仮定するときは、乱数によって与えられる3つの座標軸の回転角θx,θy,θzから決まるn組の座標変換マトリックスを用いて、予め設定されたMPの成長方向を仮定するときは、当該成長方向に応じて予め設定されたn個の所定の座標変換マトリックスを用いる。図8aにMPの成長方向をランダムな成長方向とした場合の例を示し、図8bにMPの成長方向を規則的な成長方向とした場合の例を示す。SPの配置とMPの成長方法とが決まったら、MPを次の手順で成長させる変形体を構成する粒子を新たに配置する。1成長させるMPの選択:整数乱数又は予め設定される所定のルールに基づいて、成長させるMPを決める。例えば、MPに一意に特定可能な番号情報を設定しておき、成長させるMPの番号を決める。2成長させる格子点の選択:1つのMPにSPを原点として、MPの成長方向に座標軸が定められた面心立方格子を仮定しSPから順に隣接する面心立方格子点に成長粒子GPと呼ぶを配置させる。1つのMPにSPを含む複数のGPが存在するとき、成長させるべきGPを整数乱数又は予め設定される所定のルールによって決める。例えば、SP及びGPに一意に特定可能な番号情報を設定しておき、成長させるSP及びGPの番号を決める。3成長方向の選択:1つの格子点は12個の隣接格子点を有するため、12個の中から1つを乱数又は予め設定される所定のルールによって成長方向を決める。このとき、12個の隣接格子点のうち既にGPが配置されている隣接格子点は選択肢から排除する。また、ある格子点の12個の隣接格子点が全てGPで占有されている場合は、その格子点は飽和しているものとして飽和フラグSF=1として、その格子点の成長を終了させる。4MP間の成長競争:もし、あるMPのあるGPのある成長方向が、上述した1~3の手順で選択されたときに、その成長方向には既に別のMPのGPがあり、そのGPとの重なりが生じる場合は、その方向の成長はできないものとする。なお、記憶部11には、変形体に含まれる仮想的な粒子の大きさの情報例えば、粒子の半径の情報が予め記憶されており、モデル生成部12は、当該情報を参照してGP粒子間の重なり合いの判断を行う。5MPの成長限界:もし、あるMPに特定の形状が与えられているときには、上述した1~4の手順で成長させるGPがその形状外に出てしまうならば、その方向は成長できないものとする。なお、記憶部11には、MPの形状を示す情報が予め記憶されており、モデル生成部12は、当該情報を参照して上記の判断を行う。6MPの成長終了:MPの成長回数に制限がある場合、全てのMPに対してその成長回数に達した場合、あるいは全てのMPに対するSFが1になった場合、設定された終了条件を満たしたものとしてMPの成長を終了させる。なお、記憶部11には、設定された終了条件を示す情報例えば、成長回数が予め記憶されており、モデル生成部12は、当該情報を参照して上記の判断を行う。図9に、所定領域直方体の領域にSPの座標、SPの成長方向、成長MPの選択、成長格子の選択、及び成長方向の選択を全て乱数で与えて、所定の成長回数で止めた場合のMPの成長の様子を示す。図9aに示す状態から図9cに示す状態まで、徐々にMPが成長している。以上の処理で決められた各MPについて、面心立方格子上に配置されたGPの情報から、上述した四面体を構成する処理を実行することで、多数のMPに対して自動的に四面体で分割されたモデルを生成することができる。なお、1つのユニットの大きさについては、変形体の大きさや解析の目的等に応じて予め設定されている。なお、図9に示すように生成したモデルの情報を解析装置10が備えるディスプレイに表示させてもよい。モデル生成部12は、生成したモデルの情報を記憶部11に格納して、後の解析に利用できるようにする。モデルの情報としては、具体的には例えば、解析の対象である変形体に対応付けられた仮想的な粒子の三次元座標と、当該粒子が構成する四面体の辺の情報どの粒子とどの粒子とが辺を構成しているかの情報とである。また、仮想的な粒子のうち、変形体の表面に位置する粒子については、フラグが付与される等、表面に位置するものと判断できるようになっている。なお、モデル生成部12を主要な構成として実現されるモデル生成の機能は、それ自体が、本発明の実施形態に係る解析等とは独立した、モデル生成装置として構成されることとしてもよい。また、解析に用いられる情報として、以下の情報も記憶部11に予め記憶されている。MP四面体メッシュ変形体に対して、以下の物性が与えられている。なお、ここでは、各MPには弾性体が仮定されているケースについてのみ示すが、弾性体以外の構成関係粘弾性、弾塑性等に拡張したモデルに拡張することも可能である。変形体の密度:ρ弾性定数:ヤング率E,ポアソン比ν減衰定数:ηtMP表面間の相互作用に関するパラメータとして、以下のものが与えられている。後述するようにMP表面間の相互作用には通常の個別要素法と同様の扱いを施すものとし、MP表面を構成する粒子同士が衝突する際の法線方向のバネ定数及び反発係数、接線方向のバネ定数及び反発係数、摩擦係数を入力値とする。粒子間法線方向バネ定数:Kn粒子間法線方向反発係数:en粒子間接線方向バネ係数:Kt粒子間接線方向反発係数:et粒子間摩擦係数:μMP表面と壁面との境界間の相互作用に関するパラメータとして、以下のものが与えられている。なお、壁面とは、解析において変形物が衝突しえる変形物以外の物体である。MP表面と壁面との境界間の相互作用には通常の個別要素法と同様の扱いを施すものとし、MP表面を構成する粒子が壁面に衝突する際の法線方向のバネ定数及び反発係数、接線方向のバネ定数及び反発係数、摩擦係数を入力値とする。粒子-壁面境界間法線方向バネ定数:Kn_wall粒子-壁面境界間法線方向反発係数:en_wall粒子-壁面境界間接線方向バネ係数:Kt_wall粒子-壁面境界間接線方向反発係数:et_wall粒子-壁面境界間摩擦係数:μ_wallMP表面間及びMP表面と壁面との境界間の接触判定には、異なるMP表面に配置される仮想粒子の中心間距離と粒子径との比較によって行われる。そこで、この接触判定に必要なパラメータとしてMP表面の仮想粒子の半径を与える。なお、仮想粒子の半径は、上述したMPの成長の際に判断に用いたものと同様のものを用いることが望ましい。仮想粒子半径:r_v数値積分の時間増分に関するパラメータとして、以下のものが与えられている。後述するように変形体を分割する四面体の頂点を中心とする仮想的な粒子が受け持つ力が算出され、頂点に分配された質量に基づいて数値積分し、各頂点の位置を更新するにあたり必要なパラメータとして数値積分の時間増分Δtを入力値として与える。なお、Δtに関しては、計算の開始から終了まで固定値とする場合の他、衝撃力を伴う現象の解析を扱う場合のように短時間に劇的に大きな力を四面体頂点に生じることが予想される場合は、可変値とすることも可能である。数値積分時間増分:Δt計算対象領域の指摘に関するパラメータとして、以下のものが与えられている。後述するように、MP表面間の接触判定を行うために表面仮想粒子に対する接触候補探索が行われる。この探索には空間のセル分割が必要なのでその範囲を予め指定する必要がある。計算領域x方向最小値、最大値:xmin,xmax計算領域y方向最小値、最大値:ymin,ymax計算領域z方向最小値、最大値:zmin,zmax重力や外部から加わる力、その他固定境界、移動境界等の条件を入力として与える。第1の力算出部13は、モデル生成部12によって生成されたモデルにおいて、仮想的な粒子のうち変形体の表面に設けられると共に互いに隣接しない粒子の間に生じる第1の力を算出する第1の力算出手段である。具体的には、第1の力算出部13は、各粒子の位置、並びに予め設定される粒子のパラメータに基づいて第1の力を算出する。予め設定される粒子のパラメータは、上述したような粒子の半径、質量及びばね係数等である。第1の力算出部13によって算出される第1の力は、変形体間の衝突や接触に係る力である。第1の力算出部13は、変形体の表面に設けられた仮想的な粒子を対象として、粒子間相互作用力を算出する。具体的には、例えば、個別要素法や分子動力学法に用いられる方法により、以下のように第1の力を算出する。まず、第1の力算出部13は、計算領域である三次元領域を複数のセルに分割する。計算領域やセルのサイズを示す情報については、予め記憶部11に記憶されており、その情報が参照される。例えば、x、y、z方向の領域幅とx、y、z方向のセル幅がそれぞれ記憶部11に記憶されている。続いて、第1の力算出部13は、仮想的な粒子の三次元座標に基づき、仮想的な粒子がどのセルに含まれるか判断する。そして、着目粒子接触候補相手の粒子を探す粒子が含まれるセルと同じセル及び隣接セルに含まれる他の粒子のうち所定の粒子間距離以内例えば粒子径×1.1倍にある粒子を接触候補として、接触候補リストを作成する。ここで、接触候補リストとは、接触候補とされた2つの粒子番号を一つのペアとした場合、全てのペアの粒子番号を並べて表記したリストである。そして接触候補リストにおいて、接触候補とされた粒子同士について、仮想的な粒子の三次元座標と粒子の半径とから、粒子同士の接触を判定する。ここで、同一の変形体において隣接する粒子四面体の同一の辺の両端に設けられた粒子同士の接触は判定しない接触候補リストに含めない。なお、同一の変形体であっても、隣接していない粒子四面体の同一の辺の両端に設けられたものではない粒子同士の接触は判定する。そして、接触したと判断された粒子について、上述したバネ定数、反発係数及び摩擦係数等のパラメータに基づいて、接触している粒子間に生じる接触力が計算されて、当該接触力に基づき個々の粒子に生じる第1の力を算出する。接触力の計算は以下の式によって法線方向と接線方向の接触力FnとFtがそれぞれ計算される。ここでxとvは相対変位および相対速度ベクトル、minは絶対値が小さい方の値を表す関数、ntは接線方向単位ベクトル、ηはダッシュポットの粘性減衰係数であり、反発係数eと次式の関連がある。ここでmpは粒子質量である。接触力には、この他に付着力や静電気力、ファンデルワールス力などを含めても良い。第1の力は、例えばベクトルによって示される。第1の力算出部13は、解析の時間時刻毎解析の時間時刻が数値積分時間増分Δt進められる毎に第1の力を算出する。第1の力算出部13は、変形体の表面に設けられる各粒子について算出した第1の力を示す情報を粒子位置更新部15に出力する。なお、接触候補リストの作成は、全ての解析上の時刻において行われるのではなく、予め設定された接触候補更新条件を満たした場合のみに行われればよい。接触候補更新条件としては、例えば、設定された回数の時刻が経過したことや、接触候補リストを作成してからの粒子の移動距離が設定された閾値を超えたことを条件とする。また、第1の力算出部13は、粒子間の接触だけでなく、変形体の表面に設けられた仮想的な粒子と壁面との間に生じる力や、流体から受ける抵抗力、電磁場から受ける力等についても算出してもよい。また、変形体が流体である場合には、表面張力を算出してもよい。また、上述したように第1の力は、変形体の表面に設けられる粒子を対象として個別要素法によって算出しえるものであるので、第1の力の算出には個別要素法に係る手法が適宜適用されてもよい。第2の力算出部14は、モデル生成部12によって生成されたモデルにおいて、変形体を分割する四面体に生じる応力を算出し、当該四面体を分割する面を求めることによって、算出した応力から当該四面体を構成する各粒子に生じる第2の力を算出する第2の力算出手段である。具体的には、第2の力算出部14は、四面体に生じる応力を、該四面体を構成する各粒子の位置の時間変化又は速度に基づいて、算出する。第2の力算出部14によって算出される第2の力は、変形体の変形に係る力である。上述したように第2の力は、変形体を分割する四面体を単位として計算される。なお、変形体の変形には、破壊が含まれていてもよい。例えば、四面体を構成する粒子間の予め定められた距離以上となった場合粒子間の近傍関係が失われた場合は、当該四面体は破壊されたものとしてもよい。破壊が生じた場合、変形体の表面に設けられたものでない変形体の内部に位置していた仮想的な粒子を、変形体の表面に設けられた粒子として扱うこととしてもよいし、新たに仮想的な粒子を配置して変形体の表面に新たな四面体を作製してもよい破壊により変形体の内部の部分が露出したものとみなす。本実施形態では、第2の力を算出するための方法をQDEMQuadrapleDiscreteElementMethodと呼ぶ。QDEMは、4体4つの仮想的な粒子間の相互作用のモデルを用いた手法である。構成則及び破壊基準は三次元である。現実の三次元物質のマクロな構成則との整合性がある。上述したようにQDEMで用いられるモデルは、4つの仮想的な粒子の中心が十分に近い、当該中心を頂点とする四面体から構成されている。QDEMでは、4粒子の相互作用から四面体に生じる1つの応力を決める。変形体が個体である場合、変形体に生じる応力とは、変形した四面体が元の形状に戻ろうとするために生じる応力である。これは弾性である。本実施形態においては、固体である場合の変形体は弾性体であるものとする。なお、四面体に生じる応力を算出することが可能であれば、変形体が弾性体以外のものであることを仮定してもよい。弾性体は、応力テンソルTtが、現時刻での変形勾配テンソルFtだけで定まる物質である。この場合、四面体に生じる応力は、弾性応力テンソルTで表すことができる。四面体に係る弾性応力テンソルTは、以下のように算出することができる。まず、変形勾配テンソルFが以下のように定義される。dx=FdXここで、図10に示すようにdxは現配置における、四面体の特定の頂点の座標から別の3つの頂点の座標へのベクトルから構成される行列である。dXは、現配置よりも前の時刻の基準配置における、四面体の特定の頂点の座標から別の3つの頂点の座標へのベクトルから構成される行列である。四面体が微小ならば、dXからdxへの変換は、線形変換である。Fの9つの成分を決めるためには3組の一次独立なベクトルの変換を調べればよい。3組の一次独立なベクトルは、上記のように四面体の4つの頂点の座標から算出できるものである。なお、基準配置と現配置との時間間隔は、上述した数値積分時間増分Δtである。基準配置からΔt後の原配置は、後述するように基準配置に基づいて粒子位置更新部15によって算出される。第2の力算出部14は、基準配置と現配置とから上記の式に基づいて、各四面体について変形勾配テンソルFを算出する。続いて、第2の力算出部14は、以下のLeftCauchy-Green変形テンソルBを求める。B=FFTこれは、Fを左極分解すると、F=VRとなる。ここで、Vは正味の変形、Rは剛体回転を表している。従って、B=VRRTV=V2となり、剛体回転の影響を受けない便利な量である。BからTへの関係を示す等方弾性体の構成式は、表示定理より以下の式により表される。T=φ0IB,IIB,IIIBIφ1IB,IIB,IIIBBφ2IB,IIB,IIIBB2ここで、IB,IIB,IIIBは、それぞれテンソルBの3つの不変量である。また、Iは単位行列である。Cayley-Hamiltonの定理を使って定式を変形して、Tが次式で与えられる。ここで、φ0、φ1、φ2を予め設定される物質定数の入力パラメータとする。T=φ0-φ2IIBIφ1φ2IBBφ2IIIBB-1微小変形、微小歪みの仮定を使わないので、破壊の寸前まで扱うことができる。第2の力算出部14は、上記の式を用いて四面体に係る弾性応力テンソルTを算出する。一方で変形体が粘弾性固体や流体である場合には、変形体に生じる応力とは、四面体内で速度勾配があるときに運動量の交換によって生じる応力である。これは粘性である。この場合、四面体に生じる応力は、粘性応力テンソルTで表すことができる。四面体に係る粘性応力テンソルTは、以下のように算出することができる。まず、速度勾配テンソルLが以下のように定義される。dv=Ldxここで、図11に示すようにdvは、現配置における四面体の特定の頂点における速度ベクトルv1を、別の3つの頂点における速度ベクトルv2,v3,v4から引いた3つの相対速度ベクトルdva=v2-v1、dvb=v3-v1、dvc=v4-v1から構成される行列である。四面体が微小ならば、dxからdvへの変換は、線形変換である。Lの9つの成分を決めるためには3組の一次独立なベクトルの変換を調べればよい。3組の一次独立なベクトルは、上記のように四面体の4つの頂点の座標と相対速度とから算出できるものである。なお、現配置における四面体の特定の頂点における速度ベクトルは、後述するように粒子位置更新部15によって算出される。第2の力算出部14は、基準配置と現配置とから上記の式に基づいて、各四面体について速度勾配テンソルLを算出する。速度勾配テンソルLは、以下のように加算分解することができる。L=DWD=1/2LLT歪速度テンソルW=1/2L-LT回転速度テンソルDからTへの関係を示す構成式は、表示定理より以下の式により表される。T=-pψ0Iψ1Dψ2D2ここで、pは圧力項であり、予め設定された体積圧縮係数に四面体の体積変化率J=detFを乗じることによって得られる圧力増分δpと、予め設定される初期圧力p0から求める。ここで、ψ0、ψ1、ψ2を予め設定される物質定数粘性係数の入力パラメータとする。第2の力算出部14は、上記の式を用いて四面体に係る弾性応力テンソルTを算出する。第2の力算出部14は、各変形体の各四面体について、上述のようにして求めた応力テンソルを、当該四面体の4つの頂点を中心とする粒子に分配し、各粒子に生じる第2の力を算出する。具体的には、以下のように四面体を4つの等体積のブロックに分割して、分割面の法線ベクトルと面積とを算出する。まず、四面体の各辺の重心2点の重心だから2重心と呼ぶの座標を算出する。四面体には辺が6本あるので6点ある。続いて、四面体の各面の重心3点の重心だから3重心と呼ぶの座標を算出する。四面体には面が4個あるので4点ある。続いて、四面体の重心4点の重心だから4重心と呼ぶの座標を算出する。ここで、4重心と4つの3重心の中から選ばれた2つの3重心との3点からなる三角形は6個ある4C2。2重心とその2重心を含む2面の3重心とからなる三角形も6個ある6C1。これら12個の三角形で四面体に切れ目を入れると、図12に示すように4面体は等体積の4つのブロックに分割できる。1つの粒子iに加わる力fiは、四面体の応力テンソルT弾性応力テンソルT、粘性応力テンソルTと、四面体を4つのブロックに分割する三角形とから求められる。第2の力算出部14は、四面体を4つのブロックに分割する全ての三角形の単位法線ベクトルnと、当該三角形の面積dSとを算出する。第2の力算出部14は、四面体の応力テンソルTと単位法線ベクトルnとから、以下の式によって四面体内のある微小断面四面体を4つのブロックに分割する三角形上での応力ベクトルtnを算出する。tn=TTnその微小断面の面積がdSであると、第2の力算出部14は、面に働く力dfを以下の式により算出する。df=tndS1つの粒子iに加わる力fiは、fi=Σdfである。ここで、Σは粒子iが頂点となっている全ての四面体から受ける力を足し合わせることを意味する。第2の力算出部14は、上記の式から個々の粒子に加わる力を算出する。上述したように、第2の力算出部14は、四面体の各頂点の物理量と位置とから、その勾配テンソルを算出する。ここで、勾配テンソルは、弾性応力に対しては変形勾配テンソルであり、粘性応力に対しては速度勾配テンソルである。第2の力算出部14は、算出した勾配テンソルを独立変数とするテンソル値関数である構成式を使って四面体に生じる応力を算出する。第2の力算出部14は、応力テンソルと四面体内を分割する三角形の面の法線ベクトルから、当該三角形の面に生じる応力ベクトルを算出する。第2の力算出部14は、算出した応力ベクトルと当該三角形の面の面積とから、四面体を構成する各粒子に生じる第2の力を算出する。第2の力は、例えばベクトルによって示される。第2の力算出部14は、解析の時間時刻毎解析の時間時刻が数値積分時間増分Δt進められる毎に第2の力を算出する。第2の力算出部14は、変形体の表面に設けられる各粒子について算出した第2の力を示す情報を粒子位置更新部15に出力する。粒子位置更新部15は、変形体の各粒子について、第1の力算出部13によって算出された第1の力と、第2の力算出部14によって算出された第2の力との和を取ることで各粒子に生じる力を算出して、算出した力に基づいて運動方程式を解くことで各粒子の位置を更新する粒子位置更新手段である。具体的には、第1の力と第2の力との和を取って算出された各粒子に生じる力をもとに、粒子の運動方程式f=ma、f:各粒子に加わる力、m:粒子の質量、a:粒子の速度を数値積分して、数値積分時間増分Δtだけ進められた時刻の粒子の位置座標及び速度を更新する。この際に、予め記憶部11に記憶されている粒子の質量密度から計算されてもよい等の物理量が参照される。なお、第1の力は変形体の表面の粒子のみについて算出されているため、和が取られるのは変形体の表面の粒子についてのみである。変形体の表面以外の粒子については、第2の力のみが考慮される。また、各粒子に加わる力として、第1の力及び第2の力だけでなく、予め設定された重力や外部からの力が粒子に加わることとして、粒子位置を計算することとしてもよい。計算された粒子の位置座標及び速度は、次の時刻数値積分時間増分Δtだけ進められた時刻の粒子の位置座標及び速度として、記憶部11に記憶され、次の時刻の計算に用いられる。粒子位置更新部15は、また、解析の終了条件を満たしているか否かを判断して、満たしていると判断した場合は、解析を終了させる。終了条件は、例えば、計算を行う時刻の回数が予め設定された回数実行されたことや、粒子の位置座標が収束している等の条件が用いられる。終了条件については、予め記憶部11に記憶され、粒子位置更新部15に参照される。出力部16は、解析装置10による解析の結果を出力する手段である。例えば、出力部16は、粒子位置更新部15によって計算されて記憶部11に記憶されている各時刻における粒子の位置及び速度を示す情報、即ち、変形体の状態を示す情報を出力する。出力は、例えば、解析装置10の利用者が、解析の結果を確認できるように解析装置10が備えるディスプレイに表示することによって行われる。あるいは、他の装置に解析結果の情報を出力することとしてもよい。以上が、解析装置10の構成である。引き続いて、図13のフローチャートを用いて、本実施形態に係る解析装置10で実行される処理解析方法を説明する。この処理は、解析装置10の利用者が、解析装置10に対して解析を開始する操作を行うことで開始される。解析装置10では、モデル生成部12によって、上述したように、解析の対象となる変形体のモデルが生成されるS01、モデル生成ステップ。当該モデルは、三次元形状の変形体に、その中心が変形体を分割する四面体の頂点を構成する複数の仮想的な粒子が設けられたものである。生成されたモデルの情報は、記憶部11に記憶され、必要に応じて各構成要素第1の力算出部13、第2の力算出部14、粒子位置更新部15及び出力部16から参照される。続いて、第1の力算出部13によって、モデル生成部12によって生成されたモデルにおいて、第1の力の算出対象となる変形体の表面に設けられる粒子を対象として、接触候補リストが生成されるS02、第1の力算出ステップ。なお、隣接する粒子四面体の辺の両端となっている粒子同士は、接触候補とはならない。続いて、第1の力算出部13によって、接触候補リストにおいて接触候補となっている粒子同士について接触しているかが判断され、接触している粒子同士についての接触力が計算されるS03、第1の力算出ステップ。続いて、第1の力算出部13によって、当該接触力に基づいて個々の粒子に生じる第1の力が算出されるS04、第1の力算出ステップ。算出された第1の力を示す情報は、第1の力算出部13から粒子位置更新部15に出力される。S02の処理の後、一方では第2の力算出部14によって、QDEMにより、モデル生成部12によって生成されたモデルにおいて変形体を分割する四面体に生じる応力が、当該四面体を構成する各粒子の位置の時間変化又は速度に基づいて算出されるS05、第2の力算出ステップ。続いて、第2の力算出部14によって、当該四面体を分割する面が求められて、当該面が用いられて四面体に生じる応力から、個々の粒子に生じる第2の力が算出されるS06、第2の力算出ステップ。算出された第2の力を示す情報は、第2の力算出部14から粒子位置更新部15に出力される。なお、上述したS03及びS04の処理と、S05及びS06の処理とを並列に行うこととしていたが、必ずしも並列で行う必要はなく、何れか一方の処理を行う等順番に処理が行われてもよい。第1の力及び第2の力が算出されると、続いて、粒子位置更新部15によって、各粒子について第1の力と第2の力とが足し合わせられるS07、粒子位置更新ステップ。続いて、粒子位置更新部15によって、足し合わされた力に基づいて、粒子の運動方程式の数値積分が行われて、次の時刻数値積分時間増分Δtだけ進められた時刻の各粒子の位置座標及び速度が計算されて更新されるS08、粒子位置更新ステップ。続いて、粒子位置更新部15によって、解析の終了条件を満たしているか否かが判断されるS09。ここで、終了条件とは、例えば、解析時間シミュレーション時間の上限値であったり粒子の座標変化に対する下限値であったりしてもよい。解析の終了条件を満たしていないと判断されるS09のNOと、以下のように解析が継続される。まず、第1の力算出部13によって、接触候補更新条件が満たされているか否かが判断されるS10。接触候補更新条件が満たされていると判断されるS10のYESと、続いて、第1の力算出部13によって、変形体の表面に設けられる粒子を対象として、再度、接触候補リストが生成されるS02。接触候補リストは、S08によって更新された粒子の最新の位置座標に基づいて行われる。S10において接触候補更新条件が満たされていないと判断された場合S10のNO、あるいは、S02において再度、接触候補リストが生成された場合、続いて、上述したS03及びS05以降の処理が繰り返される。なお、第1の力及び第2の力の算出は、S08によって更新された粒子の最新の位置座標及び速度に基づいて行われる。S09において、解析の終了条件が満たされていると判断されるS09のYESと、出力部16によって解析の結果が出力されて、処理が終了するS11、出力ステップ。なお、解析の結果の出力は、処理の最後のみに行われるのではなく、S08の処理の後毎に各粒子の位置座標及び速度を示す情報が出力される等によって行われてもよい。また、S11の出力と同時に記憶部11に各粒子の位置座標及び速度を示す情報等が記憶されてもよい。以上が、本実施形態に係る解析装置10で実行される処理解析方法である。上述したように本実施形態では、個別要素法等の方法により、変形体間の衝突や接触に係る第1の力が算出され、また、QDEMにより、変形体の変形に係る第2の力が算出され、それらの和の力に基づいて各粒子の位置が更新される。また、第1の力の算出と第2の力の算出との何れにおいても、変形体に設けられた仮想的な粒子が用いられる。従って、本実施形態では、第1の力の算出と第2の力の算出との両方が、上記の仮想的な粒子という共通の項目に対する計算にて行うことができる。これにより、本実施形態によれば、変形体の変形と変形体間の衝突とを別々の処理で計算する従来の方法と比べて、変形体間の複雑な接触判定や接触力の計算を変形体表面の四面体メッシュ頂点に配置された仮想的な表面粒子間の接触問題として扱うことで個別要素法等によって高速に行うことができ、更に、変形体の変形量や変形により生じる力は四面体メッシュの頂点である仮想的な粒子の座標や速度といった仮想的な粒子が持つ情報等から計算され、その力は変形体間の接触力を計算するときに使用した共通の仮想的な粒子に加わる力として処理される。従って、変形体間の衝突や接触と変形体自体の変形との両方を共通の仮想的な粒子の運動のみを計算することで処理することができるため、従来のFEMとDEMを組み合わせた手法等で見られる様な変形体の変形による力の計算と変形体間の接触力の計算との間で必要な情報のやり取りが少なく、変形体間の接触力を変形体の頂点に振り分ける処理も必要で無いため、変形体が多い場合でも、変形体間の衝突や接触と変形体自体の変形との両方を考慮して計算負荷が低く高速に解析を行うことができる。また、本手法を並列計算機を用いて並列計算を行う場合、従来のFEMとDEMを組み合わせた手法と比較して、FEMの計算に用いられる様な変形体の全体剛性マトリックス等の対象領域全てに関する巨大なマトリックスを作る必要が全くないので、並列計算で全ノード間通信の必要が全くない点で高速計算するにあたって有利である。また、上述した粒子を成長させて変形体のモデルを生成することとすれば、多数の変形体に対して規則的で品質のよい四面体が設けられたモデルを効率的に生成することができる。具体的には変形体が多い場合、上述したポリゴンを用いる方法では、変形体の形状に合うように変形体に四面体メッシュを設けるのは計算時間がかかるが、本実施形態に係る粒子を成長させる方法によれば、ポリゴンを用いる方法に比べて変形体を分割する四面体を効率的に生成することができる。但し、変形体の形状が、ポリゴンを用いた場合でも効率的に四面体を生成できる形状の場合等には、必ずしも粒子を成長させて変形体のモデルを生成する必要はない。また、解析に用いる粒子のパラメータとして、上述したように粒子の半径、質量及びバネ係数を用いることができる。この構成によれば、適切かつ確実に第1の力を算出することができる。引き続いて、本実施形態に係る解析装置10解析方法で実際に解析を行った例を示す。この例は、固定床面から0.5mの高さの空中に置かれたマットレスと同じ材料の球をマットレスと同時に落下させる解析である。この解析では、マットレスの床との衝突及び玉がマットレスに衝突することによる、マットレスの変形及び反発特性を同時に調べることができる。この例では、以下のようにパラメータの値を与えた。密度:ρ=2000.0、弾性定数:ヤング率E,ポアソン比ν=1.0×107,0.25、減衰定数:ηt=0.001、粒子間法線方向バネ定数:Kn=1.0×106、粒子間法線方向反発係数:en=0.8、粒子間接線方向バネ係数:Kt=1.6×106、粒子間接線方向反発係数:et=0.8、粒子間摩擦係数:μ=0.6、粒子-壁面境界間法線方向バネ定数:Kn_wall=1.0×106、粒子-壁面境界間法線方向反発係数:en_wall=0.8、粒子-壁面境界間接線方向バネ係数:Kt_wall=1.0×106、粒子-壁面境界間接線方向反発係数:et_wall=0.8、粒子-壁面境界間摩擦係数:μ_wall=0.6、仮想粒子半径:r_v=0.01、数値積分時間増分:Δt=1.0×10-5、計算領域x方向最小値、最大値:xmin,xmax=0.0,5.0、計算領域y方向最小値、最大値:ymin,ymax=0.0,5.0、計算領域z方向最小値、最大値:zmin,zmax=0.0,5.0、重力加速度=9.8。図14aに初期状態を示す。図の濃淡は、四面体に生じる歪みの絶対値を示す。図14bにマットレスと床面とが衝突した状態を示す。図14cに上昇するマットレスと球とが1度目に衝突した状態を示す。図14dにマットレスと球とが2度目に衝突した状態を示す。図14eにマットレスと球との3度目以降の衝突があり球が転がっている状態を示す。図14fに最終状態を示す。この例に示すように、本実施形態に係る解析装置10解析方法では、変形体間の衝突と変形体の変形とが考慮された解析が実現できている。なお、図14に示すように生成した解析の結果を示す情報を解析装置10が備えるディスプレイに表示させてもよい。引き続いて、上述した一連の解析装置10の、1つ以上の変形体の運動の解析を行う処理をコンピュータに実行させるための解析プログラムを説明する。図15に示すように、解析プログラム21は、コンピュータに挿入されてアクセスされる、あるいはコンピュータが備える記録媒体20に形成されたプログラム格納領域20a内に格納される。解析プログラム21は、解析処理を統括的に制御するメインモジュール21aと、記憶モジュール21bと、モデル生成モジュール21cと、第1の力算出モジュール21dと、第2の力算出モジュール21eと、粒子位置更新モジュール21fと、出力モジュール21gとを備えて構成される。記憶モジュール21bと、モデル生成モジュール21cと、第1の力算出モジュール21dと、第2の力算出モジュール21eと、粒子位置更新モジュール21fと、出力モジュール21gとを実行させることにより実現される機能は、上述した解析装置10の記憶部11と、モデル生成部12と、第1の力算出部13と、第2の力算出部14と、粒子位置更新部15と、出力部16との機能とそれぞれ同様である。なお、解析プログラム21は、その一部若しくは全部が、通信回線等の伝送媒体を介して伝送され、他の機器により受信されて記録インストールを含むされる構成としてもよい。また、解析プログラム21の各モジュールは、1つのコンピュータでなく、複数のコンピュータのいずれかにインストールされてもよい。その場合、当該複数のコンピュータによるコンピュータシステムよって上述した一連の解析プログラム21の1つ以上の変形体の運動の解析を行う処理が行われる。また、解析プログラム21をネットワークに接続されたサーバに格納しておき、所定のハードウェア装置にダウンロードしてもよい。"