(57)【特許請求の範囲】
【請求項１】
１つ以上の変形体の運動の解析を行う解析装置であって、
その中心が前記変形体を分割する四面体の頂点を構成する複数の仮想的な粒子が設けられた当該変形体のモデルを生成するモデル生成手段と、
前記モデル生成手段によって生成されたモデルにおいて、前記粒子のうち前記変形体の表面に設けられると共に互いに隣接しない粒子の間に生じる第１の力を算出する第１の力算出手段と、
前記モデル生成手段によって生成されたモデルにおいて、前記変形体を分割する四面体に生じる応力を算出し、当該四面体を分割する面を求めることによって、算出した応力から当該四面体を構成する各粒子に生じる第２の力を算出する第２の力算出手段と、
各粒子について、前記第１の力算出手段によって算出された第１の力と、前記第２の力算出手段によって算出された第２の力との和を取ることで各粒子に生じる力を算出して、算出した力に基づいて運動方程式を解くことで各粒子の位置を更新する粒子位置更新手段と、
を備える解析装置。
【請求項２】
前記モデル生成手段は、変形体の数に応じて配置された前記粒子と予め設定した位置関係にあると共に前記四面体の頂点を構成する粒子を設定された終了条件を満たすまで更に配置することで、前記変形体のモデルを生成することを特徴とする請求項１に記載の解析装置。
【請求項３】
前記第１の力算出手段は、第１の力を、各粒子の位置、並びに予め設定される粒子のパラメータに基づいて算出することを特徴とする請求項１又は２に記載の解析装置。
【請求項４】
予め設定される粒子のパラメータは、粒子の半径、質量及びバネ係数であることを特徴とする請求項３に記載の解析装置。
【請求項５】
前記第２の力算出手段は、前記変形体を分割する四面体に生じる応力を、当該四面体を構成する各粒子の位置の時間変化又は速度に基づいて算出することを特徴とする請求項１〜４の何れか一項に記載の解析装置。
【請求項６】
１つ以上の変形体の運動の解析を行う解析装置による解析方法であって、
その中心が前記変形体を分割する四面体の頂点を構成する複数の仮想的な粒子が設けられた当該変形体のモデルを生成するモデル生成ステップと、
前記モデル生成ステップにおいて生成されたモデルにおいて、前記粒子のうち前記変形体の表面に設けられると共に互いに隣接しない粒子の間に生じる第１の力を算出する第１の力算出ステップと、
前記モデル生成ステップにおいて生成されたモデルにおいて、前記変形体を分割する四面体に生じる応力を算出し、当該四面体を分割する面を求めることによって、算出した応力から当該四面体を構成する各粒子に生じる第２の力を算出する第２の力算出ステップと、
各粒子について、前記第１の力算出ステップにおいて算出された第１の力と、前記第２の力算出ステップにおいて算出された第２の力との和を取ることで各粒子に生じる力を算出して、算出した力に基づいて運動方程式を解くことで各粒子の位置を更新する粒子位置更新ステップと、
を含む解析方法。
【請求項７】
コンピュータに、１つ以上の変形体の運動の解析を行わせる解析プログラムであって、
前記コンピュータを、
その中心が前記変形体を分割する四面体の頂点を構成する複数の仮想的な粒子が設けられた当該変形体のモデルを生成するモデル生成手段と、
前記モデル生成手段によって生成されたモデルにおいて、前記粒子のうち前記変形体の表面に設けられると共に互いに隣接しない粒子の間に生じる第１の力を算出する第１の力算出手段と、
前記モデル生成手段によって生成されたモデルにおいて、前記変形体を分割する四面体に生じる応力を算出し、当該四面体を分割する面を求めることによって、算出した応力から当該四面体を構成する各粒子に生じる第２の力を算出する第２の力算出手段と、
各粒子について、前記第１の力算出手段によって算出された第１の力と、前記第２の力算出手段によって算出された第２の力との和を取ることで各粒子に生じる力を算出して、算出した力に基づいて運動方程式を解くことで各粒子の位置を更新する粒子位置更新手段と、
して機能させる解析プログラム。
【請求項８】
コンピュータに、１つ以上の変形体の運動の解析を行わせる解析プログラムを記録したコンピュータ読み取り可能な記録媒体であって、
前記解析プログラムは、
前記コンピュータを、
その中心が前記変形体を分割する四面体の頂点を構成する複数の仮想的な粒子が設けられた当該変形体のモデルを生成するモデル生成手段と、
前記モデル生成手段によって生成されたモデルにおいて、前記粒子のうち前記変形体の表面に設けられると共に互いに隣接しない粒子の間に生じる第１の力を算出する第１の力算出手段と、
前記モデル生成手段によって生成されたモデルにおいて、前記変形体を分割する四面体に生じる応力を算出し、当該四面体を分割する面を求めることによって、算出した応力から当該四面体を構成する各粒子に生じる第２の力を算出する第２の力算出手段と、
各粒子について、前記第１の力算出手段によって算出された第１の力と、前記第２の力算出手段によって算出された第２の力との和を取ることで各粒子に生じる力を算出して、算出した力に基づいて運動方程式を解くことで各粒子の位置を更新する粒子位置更新手段と、
して機能させる記録媒体。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、１つ以上の変形体の運動の解析を行う解析装置、解析方法、解析プログラム及び記録媒体に関する。
【背景技術】
【０００２】
従来から、衝突や接触がある物体（要素）の運動解析を行う方法として個別要素法（ＤＥＭ）等の手法が知られている。更に、解析の対象となる要素の変形を考慮した不連続変形法（Discontinuous Deformation Analysis：ＤＤＡ）と呼ばれる方法（例えば、特許文献１参照）や、ＱＤＥＭ（Quadraple Discrete Element Method）と呼ばれる方法が知られている（例えば、非特許文献１参照）。
【先行技術文献】
【特許文献】
【０００３】
【特許文献１】特開平７−１８１８７９号公報
【非特許文献】
【０００４】
【非特許文献１】Hide Sakaguchi、“New computational scheme of discrete element approach for the solidearth multi-materials simulation using three-dimensional four particleinteraction”、FRONTIER RESEARCH ON EARTH EVOLUTION、VOL.2、(2004)
【発明の概要】
【発明が解決しようとする課題】
【０００５】
ところで、砂利やワイヤーやケーブルは、個々の石や鋼線が多数集まったものであるが、個々の石や鋼線の変形に関する性質が全体の挙動に大きく影響する。しかし、その集合体としての性質や、個々の石や鋼線の性質だけからは決まらない。例えば、砂利ならばその敷きならし方や盛り上げた形状に、ワイヤーやケーブルならば鋼線の撚り方によって全体の性質は大きく変化する。
【０００６】
このような物体に対しては、個々の石や鋼線等の要素間の衝突や接触と、個々の要素自体の変形との両方を考慮しなければ、適切な解析、例えば要素内の弾性波伝播や破壊に伴う運動を正確にシミュレーションすることができない。石や鋼線等の要素は接触圧が大きく摩擦力が大きい。即ち、上記の接触や衝突としては、具体的には、摩擦力やすべり特性を考慮する必要がある。
【０００７】
しかしながら、上述した不連続変形法では、変形を考慮しているものの個々の要素についての変形に係る計算を行った上で衝突や接触の計算を行うため、個々の石や鋼線等の要素が多くなると、計算負荷が非常に大きくなるため計算時間が大きくなってしまうという問題点があった。そのため、このような物体に対しては経験的な評価しか行えないのが原状である。また、ＤＥＭでは、複数個の球形粒子を連結して一つの変形要素をモデル化する方法も有るが、正しい物性を再現する構成則が存在しないために要素自体の変形や要素内部の応力解析を正しく行うことができない。更に、ＱＤＥＭでは、単一要素の変形や要素内の応力を構成則に基づいて解析することは可能であるが、複数の要素間の接触や衝突判定方法及び相互作用力の算出機能を有していないため、ここで目的とする多数の変形する要素の運動を解析することは不可能である。
【０００８】
本発明は、上記の問題点を鑑みてなされたものであり、対象とする要素である変形体（物体）が多い場合でも、変形体間の衝突や接触と変形体自体の変形との両方を考慮して高速に解析を行うことができる解析装置、解析方法、解析プログラム及び記録媒体を提供することを目的とする。
【課題を解決するための手段】
【０００９】
上記の目的を達成するために、本発明に係る解析装置は、１つ以上の変形体の運動の解析を行う解析装置であって、その中心が変形体を分割する四面体の頂点を構成する複数の仮想的な粒子が設けられた当該変形体のモデルを生成するモデル生成手段と、モデル生成手段によって生成されたモデルにおいて、粒子のうち変形体の表面に設けられると共に互いに隣接しない粒子の間に生じる第１の力を算出する第１の力算出手段と、モデル生成手段によって生成されたモデルにおいて、変形体を分割する四面体に生じる応力を算出し、当該四面体を分割する面を求めることによって、算出した応力から当該四面体を構成する各粒子に生じる第２の力を算出する第２の力算出手段と、各粒子について、第１の力算出手段によって算出された第１の力と、第２の力算出手段によって算出された第２の力との和を取ることで各粒子に生じる力を算出して、算出した力に基づいて運動方程式を解くことで各粒子の位置を更新する粒子位置更新手段と、を備える。
【００１０】
本発明に係る解析装置では、変形体間の衝突や接触に係る第１の力と変形体の変形に係る第２の力とがそれぞれ算出されて、それらの和の力に基づいて各粒子の位置が更新される。また、第１の力の算出と第２の力の算出との何れにおいても、変形体に設けられた仮想的な粒子が用いられる。従って、本発明に係る解析装置では、第１の力の算出と第２の力の算出との両方が、上記の仮想的な粒子という共通の項目に対する計算にて行うことができる。これにより、本発明に係る解析装置によれば、変形体が多い場合でも、変形体間の衝突や接触と変形体自体の変形との両方を考慮して高速に解析を行うことができる。
【００１１】
モデル生成手段は、変形体の数に応じて配置された粒子と予め設定した位置関係にあると共に四面体の頂点を構成する粒子を設定された終了条件を満たすまで更に配置することで、変形体のモデルを生成することとしてもよい。この構成によれば、多数の変形体に対して規則的で品質のよい四面体が設けられたモデルを効率的に生成することができる。
【００１２】
第１の力算出手段は、第１の力を、各粒子の位置、並びに予め設定される粒子のパラメータに基づいて算出することとしてもよい。予め設定される粒子のパラメータは、粒子の半径、質量及びバネ係数であることとしてもよい。これらの構成によれば、適切かつ確実に第１の力を算出することができる。
【００１３】
前記第２の力算出手段は、前記変形体を分割する四面体に生じる応力を、当該四面体を構成する各粒子の位置の時間変化又は速度に基づいて算出することとしてもよい。この構成によれば、適切かつ確実に第２の力を算出することができる。
【００１４】
ところで、本発明は、上記のように解析装置の発明として記述できる他に、以下のように解析方法及び解析プログラムの発明としても記述することができる。これはカテゴリ等が異なるだけで、実質的に同一の発明であり、同様の作用及び効果を奏する。
【００１５】
本発明に係る解析方法は、１つ以上の変形体の運動の解析を行う解析装置による解析方法であって、その中心が変形体を分割する四面体の頂点を構成する複数の仮想的な粒子が設けられた当該変形体のモデルを生成するモデル生成ステップと、モデル生成ステップにおいて生成されたモデルにおいて、粒子のうち変形体の表面に設けられると共に互いに隣接しない粒子の間に生じる第１の力を算出する第１の力算出ステップと、モデル生成ステップにおいて生成されたモデルにおいて、変形体を分割する四面体に生じる応力を算出し、当該四面体を分割する面を求めることによって、算出した応力から当該四面体を構成する各粒子に生じる第２の力を算出する第２の力算出ステップと、各粒子について、第１の力算出ステップにおいて算出された第１の力と、第２の力算出ステップにおいて算出された第２の力との和を取ることで各粒子に生じる力を算出して、算出した力に基づいて運動方程式を解くことで各粒子の位置を更新する粒子位置更新ステップと、を含む。
【００１６】
本発明に係る解析プログラムは、コンピュータに、１つ以上の変形体の運動の解析を行わせる解析プログラムであって、コンピュータを、その中心が変形体を分割する四面体の頂点を構成する複数の仮想的な粒子が設けられた当該変形体のモデルを生成するモデル生成手段と、モデル生成手段によって生成されたモデルにおいて、粒子のうち変形体の表面に設けられると共に互いに隣接しない粒子の間に生じる第１の力を算出する第１の力算出手段と、モデル生成手段によって生成されたモデルにおいて、変形体を分割する四面体に生じる応力を算出し、当該四面体を分割する面を求めることによって、算出した応力から当該四面体を構成する各粒子に生じる第２の力を算出する第２の力算出手段と、各粒子について、第１の力算出手段によって算出された第１の力と、第２の力算出手段によって算出された第２の力との和を取ることで各粒子に生じる力を算出して、算出した力に基づいて運動方程式を解くことで各粒子の位置を更新する粒子位置更新手段と、して機能させる。
【００１７】
本発明に係る記録媒体は、コンピュータに、１つ以上の変形体の運動の解析を行わせる解析プログラムを記録したコンピュータ読み取り可能な記録媒体であって、解析プログラムは、コンピュータを、その中心が変形体を分割する四面体の頂点を構成する複数の仮想的な粒子が設けられた当該変形体のモデルを生成するモデル生成手段と、モデル生成手段によって生成されたモデルにおいて、粒子のうち変形体の表面に設けられると共に互いに隣接しない粒子の間に生じる第１の力を算出する第１の力算出手段と、モデル生成手段によって生成されたモデルにおいて、変形体を分割する四面体に生じる応力を算出し、当該四面体を分割する面を求めることによって、算出した応力から当該四面体を構成する各粒子に生じる第２の力を算出する第２の力算出手段と、各粒子について、第１の力算出手段によって算出された第１の力と、第２の力算出手段によって算出された第２の力との和を取ることで各粒子に生じる力を算出して、算出した力に基づいて運動方程式を解くことで各粒子の位置を更新する粒子位置更新手段と、して機能させる。
【発明の効果】
【００１８】
本発明では、変形体間の衝突や接触に係る第１の力の算出と変形体の変形に係る第２の力の算出との両方が、仮想的な粒子という共通の項目に対する計算にて行うことができる。これにより、本発明によれば、変形体が多い場合でも、変形体間の衝突や接触と変形体自体の変形との両方を考慮して高速に解析を行うことができる。
【図面の簡単な説明】
【００１９】
【図１】本発明の実施形態に係る解析装置の機能構成を示す図である。
【図２】解析で用いるモデルの生成を説明するための図であり、面心立方格子の基本ユニットを示す図である。
【図３】基本ユニットにおいて生成される正四面体を示す図である。
【図４】基本ユニットにおいて生成される四角錐（の空洞）を示す図である。
【図５】隣接する基本ユニット間において生成される八面体を示す図である。
【図６】基本ユニットにおいて生成される四面体、及び隣接する基本ユニット間において当該四面体から構成される八面体を示す図である。
【図７】八面体の分割を示す図である。
【図８】粒子の成長方向の例を示す図である。
【図９】モデルの生成の具体例を示す図である。
【図１０】四面体に係る弾性応力テンソルを説明するための図である。
【図１１】四面体に係る粘性応力テンソルを説明するための図である。
【図１２】四面体に生じるテンソルから各粒子に生じる力の計算を説明するための図である。
【図１３】本発明の実施形態に係る解析装置で実行される処理（解析方法）を示すフローチャートである。
【図１４】本発明の実施形態に係る解析装置による解析結果の例を示す図である。
【図１５】本発明の実施形態に係る解析プログラムの構成を、記録媒体と共に示す図である。
【発明を実施するための形態】
【００２０】
以下、図面と共に本発明に係る解析装置、解析方法及び解析プログラムの好適な実施形態について詳細に説明する。なお、図面の説明においては同一要素には同一符号を付し、重複する説明を省略する。
【００２１】
図１に本実施形態に係る解析装置１０を示す。解析装置１０は、１つ以上の変形体の運動の解析（シミュレーション）を行う。変形体は、例えば、石や鋼線等の物体である。変形体は、固体（弾性体あるいは粘弾性体）でもよいし、あるいは粘性流体や粘弾性流体でもよい。本実施形態に係る解析装置１０は、変形体をモデル化して、解析により時間（時刻）毎の変形体の位置（例えば、三次元座標）を算出する。その際、変形体間、あるいは１つの変形体の異なる部分の間の接触（衝突）を考慮する。また、変形体に加わる力によって、変形体が変形することも考慮する。
【００２２】
解析装置１０は、例えば、ＣＰＵ（Central Processing Unit）、ＧＰＵ（Graphics Processing Unit）、メモリ、ハードディスク、ディスプレイ等のハードウェアを備えるコンピュータとして構成される。これらの構成要素がプログラム等により動作することによって、後述する解析装置１０としての機能が発揮される。
【００２３】
図１に示すように解析装置１０は、記憶部１１と、モデル生成部１２と、第１の力算出部１３と、第２の力算出部１４と、粒子位置更新部１５と、出力部１６とを備える。
【００２４】
記憶部１１は、解析装置１０による解析に必要な情報を記憶する手段である。記憶部１１は、解析に予めパラメータとして必要である情報については予め記憶している。なお、この情報については、例えば、コントロールファイル等として解析装置１０の利用者によって予め解析装置１０に入力されている。また、解析の過程で、算出される情報については必要に応じて後述する各構成要素から、記憶部１１に入力される（書き込まれる）。記憶部１１が記憶している情報は、必要に応じて各構成要素（モデル生成部１２、第１の力算出部１３、第２の力算出部１４、粒子位置更新部１５及び出力部１６）から、参照される。具体的に、記憶部１１がどのような情報が記憶するかは、その情報の利用、算出の時点で説明する。
【００２５】
モデル生成部１２は、解析の対象となる変形体のモデルを生成するモデル生成手段である。モデル生成部１２は、変形体のモデルを、三次元形状の変形体（が存在する領域）に複数の仮想的な粒子を設けることによって生成する。当該粒子は、その中心が変形体を分割する四面体の頂点を構成する。粒子が位置する領域に上記の変形体があるものとして解析が行われる。
【００２６】
モデル生成部１２は、変形体の形状及び初期配置を示す情報を記憶部１１から読み出して、当該情報に基づいて、有限要素法（ＦＥＭ：Finite Element Method）等で利用される既存の四面体メッシュの自動生成方法によって、四面体に粒子を設けることとしてもよい。即ち、１）変形体の表面に三角形ポリゴン網を生成する点を張り、２）変形体内部空間に三角形ポリゴンと同程度の数密度を保った多数点を発生させ、３）隣接４点間の最小外接球を数値的に見出すことで四面体を逐次決定する。
【００２７】
このとき、各四面体の品質については、四面体の稜線（辺）の長さをより均一に近づけることで管理される。なお、ここで四面体の品質とは、例えば、四面体を構成する辺の長さや四面体間の大きさのばらつきを指す（ばらつきが小さい方が品質がよい）。また手順として、より大まかで比較的均質な四面体メッシュを予め見つけた後に各四面体のメッシュを段階的に細分化することで、四面体品質を保ちながら微小四面体分割を行う手法を用いてもよい。上記の方法で生成した四面体の頂点を中心とした粒子を設けて、モデルを生成する。
【００２８】
上記の一般的な四面体メッシュ生成方法を用いた方法以外でも、モデル生成部１２は、以下の方法によって変形体のモデルを生成することとしてもよい。即ち、モデル生成部１２は、変形体の数に応じて仮想的な粒子を配置し、当該粒子と予め設定した位置関係にあると共に四面体の頂点を構成する粒子を設定された終了条件を満たすまで更に配置する（新たに配置していく）ことで、変形体のモデルを生成することとしてもよい。例えば、以下のようにモデルを生成する。
【００２９】
上記の一般的な四面体メッシュ生成方法は、本来１つの物体や領域の四面体分割を目的として作られたもので、物体形状や領域形状を予め入力として与えなければならない。従って、多数の物体を対象としたり、更に個々の物体が不規則形状を有したりする場合には、物体群を自動生成しつつ個々の物体を多数の四面体に自動分割することには、本質的に不向きである。例えば、解析の対象を多数の砂利（によって構成される石）のような多数の不規則形状の変形体とする場合には、以下のような、変形体を自動発生させながら内部を四面体分割する方法を用いることが望ましい。
【００３０】
そこで本手法では、変形体内部を四面体分割するのではなく、変形体内部の点から面心立方格子とその内部に包含される２種類の四面体を用いて、解析に必要とされる変形体の形状まで成長させることで、多数の物体に対して規則的で品質のよい四面体メッシュを自動発生させる。
【００３１】
図２に面心立方格子の基本ユニットを示す。図中に示される（０，０，０）を原点とする直交座標座標系ｘ−ｙ−ｚを定義する。そこに原点を重心として一辺の長さが２である３枚の正方形をそれぞれｘｙ平面（ｚ＝０）、ｙｚ平面（ｘ＝０）、ｚｘ平面上（ｙ＝０）に考える。すると３枚の正方形の１２個の頂点（１，１，０），（−１，１，０），（−１，−１，０），（１，−１，０），（０，１，１），（０，−１，１），（０，−１，−１），（０，１，−１），（１，０，１），（−１，０，１），（−１，０，−１），（１，０，−１）は、全て原点から√２の等距離に位置することが分かる。この１２点と原点を合わせた合計１３点が面心立方格子の基本ユニットの格子点である。
【００３２】
面心立方格子の基本ユニットである１３点の格子点の中から隣接格子点距離が全て√２である格子点４点を選んで四面体を作ると、図３に示すように一辺の長さが√２の正三角形４面からなる正四面体となる。ｙｚ平面の正方形の表側（ｘ＞０）を考えると、図３に示す四面体と合同な正四面体が４つできる。同様にｙｚ平面の正方形の裏側（ｘ＜０）にも同じ四面体が４つ考えることができ、この基本ユニットの中には合計８つの正四面体が存在することが分かる。
【００３３】
上述の方法でできるｙｚ平面の表側（ｘ＞０）の領域の４つの正四面体の間には、図４（ａ）に示すように底面が１辺の長さ√２の正方形で高さが１の四角錐の空洞ができる。同様に、この四角錐の空洞はｘ＜０側にもでき、更にｚｘ平面のｙ＞０側とｙ＜０側、ｘｙ平面のｚ＞０側とｚ＜０側の合計６方向に生じる。図４（ａ）に見られる四角錐の空洞は、（０，０，０）を中心とする面心立方格子の全ての隣のユニットにも同様に生じる。例えば、（２，０，０）を中心とするユニットのｙｚ平面のｘ＜０側の４つの正四面体の間には、図４（ｂ）に示す四角錐の空洞が生じる。
【００３４】
図４（ａ）に示す空洞と、図４（ｂ）に示す空洞とは、底面で同一正方形を共有し、かつ２つの四角錐の頂点を結ぶｘ軸上にある線分（０，０，０）−（２，０，０）は、この正方形と重心で垂直に交わる。すると、図４（ａ）に示す四角錐と、図４（ｂ）に示す四角錐とを合体させてできる立体は、図５に示すように一辺の長さ√２の正三角形８枚からなる正八面体となる。
【００３５】
一方、図３に示す正四面体について、頂点を原点（０，０，０）とし底面三角形を（１，１，０）−（０，１，１）−（１，０，１）と考えるとき、底面三角形において頂点と反対側の空間にも（０，０，０）のユニットには空洞が生じている。この空洞を埋め尽くすために、図３に示したｘｙ平面（ｚ＝０）上の正方形に対してｚ方向の隣接ユニットである（０，０，２）を中心とするｘｙ平面の正方形の頂点のうち点（１，０，１）と（０，１，１）に最も隣接している（１，１，２）との位置関係を考える。するとそこには、図６（ａ）に示すような４つの頂点（１，１，２），（０，１，１），（１，０，１），（１，１，０）から構成される四面体領域ができる。
【００３６】
同様に四面体領域は、線分（１，１，０）−（１，１，２）を囲んで中心が（２，０，０），（２，２，０），（０，２，０）の３つの隣接ユニットにも生じる。これら４つの隣接ユニットにおける線分（１，１，０）−（１，１，２）を共有する四面体領域を全て足し合わせると図６（ｂ）に示すような八面体となることが分かる。この八面体表面を構成する三角形の一辺の長さは全て√２であることから、図６（ｂ）に示す八面体は正八面体となり、同時に図５に示す正八面体と合同である。図６（ｂ）に示す八面体と同様の正八面体は、図３で考えた８つの正四面体の全ての底面三角形について考えることができる。
【００３７】
ここまで示されたように、図２に示した１つの面心立方格子基本ユニットの空間は、８つの正四面体と、隣接ユニットと共有される図５に示す正八面体６つと、図６（ｂ）に示す正八面体８つとによって埋め尽くされていることが分かる。図５に示す正八面体と、図６（ｂ）に示す正八面体とは前述の通り合同であったので、この１４個の正八面体をそれぞれ四面体に分割すると、面心立方と四面体との関係が決定することになる。
【００３８】
正八面体を四面体に分割する方法は一意ではないが、例えば、元の正八面体の頂点のみを用いて新たな頂点を増やさずに４つの合同な四面体に分割することとすればよい。その方法は、分割によって新たに追加される四面体の稜（四面体を構成する三角形の辺）の方向の違いによって３通り存在する。
【００３９】
便宜上、図７（ａ）に示すように正八面体の重心を原点に置き、３組の向かい合う２つの頂点をそれぞれｘ軸、ｙ軸、ｚ軸上に乗るように配置する。すると、図７（ｂ）〜図７（ｄ）に示されるように、正八面体を分割するために追加される辺をｘ軸にとる場合、ｙ軸にとる場合、ｚ軸にとる場合に対応して、３通りの分割方法が存在することが分かる。その結果、図７（ｂ）〜図７（ｄ）に示すように何れも１つの正八面体を２面の直角二等辺三角形と２面の正三角形からなる合同な四面体４つに分割することができる。
【００４０】
図７を用いて説明した３通りの分割方法は、正八面体の頂点の位置さえ決まれば、他に何も依存することなく決めることができる。そこで、例えば、乱数によって図７を用いて説明した３通りの分割方法から１つを等確率で与えて分割方法を決めることすればよい。これにより、大域的なメッシュ全体で３種類の分割方法が均等に出現することになり、四面体分割の等方性を確保することができる。
【００４１】
１つの面心立方格子基本ユニットを囲む正八面体のうち、図５に示す１つの正八面体は２つのユニットに共有されている。このパターンが面心立方の立方体６方向の面で隣接するユニットに対して存在し、１つの正八面体は図７の何れかの方法で４つの四面体に分割される。従って、１つのユニットには、図５のパターンの四面体が４×６÷２＝１２個分が生成される。また、一方、図６（ｂ）に示す１つの正八面体は８つの隣接ユニットに共有されている。このパターンが面心立方の立方体８方向の頂点で隣接するユニットに対して存在し、１つの正八面体は図７の何れかの方法で４つの四面体に分割される。従って、１つのユニットには、図６（ｂ）のパターンの四面体が４×８÷８＝４個分が生成される。
【００４２】
以上までの説明で分かるように、図２で示した１つの面心立方格子基本ユニットの空間は、８つの正四面体と、直角二等辺三角形２枚及び正三角形２枚からなる四面体１２＋４＝１６個とに分割される。
【００４３】
上記の面心立方格子を成長させることによって、変形体のモデルを生成する。まず、ｎ個の変形体（ここでは、マクロパーティクル、ＭＰと呼ぶ）に対して、ＭＰの種（ここでは、シードパーティクル、ＳＰと呼ぶ）を所定の三次元空間領域に発生させる。ここで、ｎは予め設定される変形体の数である。ＭＰがランダムに配置されることを仮定するときは、三次元座標（ｘ，ｙ，ｚ）を与えるｎ組の一様乱数を用い、ＭＰが予め設定された場所に位置することを仮定するときは、当該位置に応じたｎ個の座標値（例えば、ＭＰの端点やＭＰの重心の位置）を与える。
【００４４】
次に、ｎ個のＳＰに対して、面心立方格子の方向（＝ＭＰの成長方向）を与える。このとき、ＭＰの成長方向にランダム性を仮定するときは、乱数によって与えられる３つの座標軸の回転角（θｘ，θｙ，θｚ）から決まるｎ組の座標変換マトリックスを用いて、予め設定されたＭＰの成長方向を仮定するときは、当該成長方向に応じて予め設定されたｎ個の所定の座標変換マトリックスを用いる。図８（ａ）にＭＰの成長方向をランダムな成長方向とした場合の例を示し、図８（ｂ）にＭＰの成長方向を規則的な成長方向とした場合の例を示す。
【００４５】
ＳＰの配置とＭＰの成長方法とが決まったら、ＭＰを次の手順で成長させる（変形体を構成する粒子を新たに配置する）。
（１）成長させるＭＰの選択：整数乱数又は予め設定される所定のルールに基づいて、成長させるＭＰを決める。例えば、ＭＰに一意に特定可能な番号（情報）を設定しておき、成長させるＭＰの番号を決める。
（２）成長させる格子点の選択：１つのＭＰにＳＰを原点として、ＭＰの成長方向に座標軸が定められた面心立方格子を仮定しＳＰから順に隣接する面心立方格子点に成長粒子（ＧＰと呼ぶ）を配置させる。１つのＭＰにＳＰを含む複数のＧＰが存在するとき、成長させるべきＧＰを整数乱数又は予め設定される所定のルールによって決める。例えば、ＳＰ及びＧＰに一意に特定可能な番号（情報）を設定しておき、成長させるＳＰ及びＧＰの番号を決める。
【００４６】
（３）成長方向の選択：１つの格子点は１２個の隣接格子点を有するため、１２個の中から１つを乱数又は予め設定される所定のルールによって成長方向を決める。このとき、１２個の隣接格子点のうち既にＧＰが配置されている隣接格子点は選択肢から排除する。また、ある格子点の１２個の隣接格子点が全てＧＰで占有されている場合は、その格子点は飽和しているものとして飽和フラグＳＦ＝１として、その格子点の成長を終了させる。
（４）ＭＰ間の成長競争：もし、あるＭＰのあるＧＰのある成長方向が、上述した（１）〜（３）の手順で選択されたときに、その成長方向には既に別のＭＰのＧＰがあり、そのＧＰとの重なりが生じる場合は、その方向の成長はできないものとする。なお、記憶部１１には、変形体に含まれる仮想的な粒子の大きさの情報（例えば、粒子の半径の情報）が予め記憶されており、モデル生成部１２は、当該情報を参照してＧＰ（粒子）間の重なり合いの判断を行う。
【００４７】
（５）ＭＰの成長限界：もし、あるＭＰに特定の形状が与えられているときには、上述した（１）〜（４）の手順で成長させるＧＰがその形状外に出てしまうならば、その方向は成長できないものとする。なお、記憶部１１には、ＭＰの形状を示す情報が予め記憶されており、モデル生成部１２は、当該情報を参照して上記の判断を行う。
（６）ＭＰの成長終了：ＭＰの成長回数に制限がある場合、全てのＭＰに対してその成長回数に達した場合、あるいは全てのＭＰに対するＳＦが１になった場合、設定された終了条件を満たしたものとしてＭＰの成長を終了させる。なお、記憶部１１には、設定された終了条件を示す情報（例えば、成長回数）が予め記憶されており、モデル生成部１２は、当該情報を参照して上記の判断を行う。
【００４８】
図９に、所定領域（直方体の領域）にＳＰの座標、ＳＰの成長方向、成長ＭＰの選択、成長格子の選択、及び成長方向の選択を全て乱数で与えて、所定の成長回数で止めた場合のＭＰの成長の様子を示す。図９（ａ）に示す状態から図９（ｃ）に示す状態まで、徐々にＭＰが成長している。以上の処理で決められた各ＭＰについて、面心立方格子上に配置されたＧＰの情報から、上述した四面体を構成する処理を実行することで、多数のＭＰに対して自動的に四面体で分割されたモデルを生成することができる。なお、１つのユニットの大きさについては、変形体の大きさや解析の目的等に応じて予め設定されている。なお、図９に示すように生成したモデルの情報を解析装置１０が備えるディスプレイに表示させてもよい。
【００４９】
モデル生成部１２は、生成したモデルの情報を記憶部１１に格納して、後の解析に利用できるようにする。モデルの情報としては、具体的には例えば、解析の対象である変形体に対応付けられた仮想的な粒子の三次元座標と、当該粒子が構成する四面体の辺の情報（どの粒子とどの粒子とが辺を構成しているかの情報）とである。また、仮想的な粒子のうち、変形体の表面に位置する粒子については、フラグが付与される等、表面に位置するものと判断できるようになっている。なお、モデル生成部１２を主要な構成として実現されるモデル生成の機能は、それ自体が、本発明の実施形態に係る解析等とは独立した、モデル生成装置として構成されることとしてもよい。
【００５０】
また、解析に用いられる情報として、以下の情報も記憶部１１に予め記憶されている。ＭＰ四面体メッシュ（変形体）に対して、以下の物性が与えられている。なお、ここでは、各ＭＰには弾性体が仮定されているケースについてのみ示すが、弾性体以外の構成関係（粘弾性、弾塑性等）に拡張したモデルに拡張することも可能である。
（変形体の）密度：ρ
弾性定数：（ヤング率Ｅ，ポアソン比ν）
減衰定数：ηｔ
【００５１】
ＭＰ表面間の相互作用に関するパラメータとして、以下のものが与えられている。後述するようにＭＰ表面間の相互作用には通常の個別要素法と同様の扱いを施すものとし、ＭＰ表面を構成する粒子同士が衝突する際の法線方向のバネ定数及び反発係数、接線方向のバネ定数及び反発係数、摩擦係数を入力値とする。
粒子間法線方向バネ定数：Ｋｎ
粒子間法線方向反発係数：ｅｎ
粒子間接線方向バネ係数：Ｋｔ
粒子間接線方向反発係数：ｅｔ
粒子間摩擦係数：μ
【００５２】
ＭＰ表面と壁面との境界間の相互作用に関するパラメータとして、以下のものが与えられている。なお、壁面とは、解析において変形物が衝突しえる（変形物以外の）物体である。ＭＰ表面と壁面との境界間の相互作用には通常の個別要素法と同様の扱いを施すものとし、ＭＰ表面を構成する粒子が壁面に衝突する際の法線方向のバネ定数及び反発係数、接線方向のバネ定数及び反発係数、摩擦係数を入力値とする。
粒子−壁面境界間法線方向バネ定数：Ｋｎ＿ｗａｌｌ
粒子−壁面境界間法線方向反発係数：ｅｎ＿ｗａｌｌ
粒子−壁面境界間接線方向バネ係数：Ｋｔ＿ｗａｌｌ
粒子−壁面境界間接線方向反発係数：ｅｔ＿ｗａｌｌ
粒子−壁面境界間摩擦係数：μ＿ｗａｌｌ
【００５３】
ＭＰ表面間及びＭＰ表面と壁面との境界間の接触判定には、異なるＭＰ表面に配置される仮想粒子の中心間距離と粒子径との比較によって行われる。そこで、この接触判定に必要なパラメータとしてＭＰ表面の仮想粒子の半径を与える。なお、仮想粒子の半径は、上述したＭＰの成長の際に判断に用いたものと同様のものを用いることが望ましい。
仮想粒子半径：ｒ＿ｖ
【００５４】
数値積分の時間増分に関するパラメータとして、以下のものが与えられている。後述するように変形体を分割する四面体の頂点（を中心とする仮想的な粒子）が受け持つ力が算出され、頂点に分配された質量に基づいて数値積分し、各頂点の位置を更新するにあたり必要なパラメータとして数値積分の時間増分Δｔを入力値として与える。なお、Δｔに関しては、計算の開始から終了まで固定値とする場合の他、衝撃力を伴う現象の解析を扱う場合のように短時間に劇的に大きな力を四面体頂点に生じることが予想される場合は、可変値とすることも可能である。
数値積分時間増分：Δｔ
【００５５】
計算対象領域の指摘に関するパラメータとして、以下のものが与えられている。後述するように、ＭＰ表面間の接触判定を行うために表面仮想粒子に対する接触候補探索が行われる。この探索には空間のセル分割が必要なのでその範囲を予め指定する必要がある。
計算領域ｘ方向最小値、最大値：ｘｍｉｎ，ｘｍａｘ
計算領域ｙ方向最小値、最大値：ｙｍｉｎ，ｙｍａｘ
計算領域ｚ方向最小値、最大値：ｚｍｉｎ，ｚｍａｘ
重力や外部から加わる力、その他固定境界、移動境界等の条件を入力として与える。
【００５６】
第１の力算出部１３は、モデル生成部１２によって生成されたモデルにおいて、仮想的な粒子のうち変形体の表面に設けられると共に互いに隣接しない粒子の間に生じる第１の力を算出する第１の力算出手段である。具体的には、第１の力算出部１３は、各粒子の位置、並びに予め設定される粒子のパラメータに基づいて第１の力を算出する。予め設定される粒子のパラメータは、上述したような粒子の半径、質量及びばね係数等である。
【００５７】
第１の力算出部１３によって算出される第１の力は、変形体間の衝突や接触に係る力である。第１の力算出部１３は、変形体の表面に設けられた仮想的な粒子を対象として、粒子間相互作用力を算出する。具体的には、例えば、個別要素法や分子動力学法に用いられる方法により、以下のように第１の力を算出する。
【００５８】
まず、第１の力算出部１３は、計算領域である三次元領域を複数のセルに分割する。計算領域やセルのサイズを示す情報については、予め記憶部１１に記憶されており、その情報が参照される。例えば、ｘ、ｙ、ｚ方向の領域幅とｘ、ｙ、ｚ方向のセル幅がそれぞれ記憶部１１に記憶されている。続いて、第１の力算出部１３は、仮想的な粒子の三次元座標に基づき、仮想的な粒子がどのセルに含まれるか判断する。そして、着目粒子（接触候補相手の粒子を探す粒子）が含まれるセルと同じセル及び隣接セルに含まれる他の粒子のうち所定の粒子間距離以内（例えば粒子径×１．１倍）にある粒子を接触候補として、接触候補リストを作成する。ここで、接触候補リストとは、接触候補とされた２つの粒子番号を一つのペアとした場合、全てのペアの粒子番号を並べて表記したリストである。そして接触候補リストにおいて、接触候補とされた粒子同士について、仮想的な粒子の三次元座標と粒子の半径とから、粒子同士の接触を判定する。ここで、同一の変形体において隣接する粒子（四面体の同一の辺の両端に設けられた粒子）同士の接触は判定しない（接触候補リストに含めない）。なお、同一の変形体であっても、隣接していない粒子（四面体の同一の辺の両端に設けられたものではない粒子）同士の接触は判定する。
【００５９】
そして、接触したと判断された粒子について、上述したバネ定数、反発係数及び摩擦係数等のパラメータに基づいて、接触している粒子間に生じる接触力が計算されて、当該接触力に基づき個々の粒子に生じる第１の力を算出する。接触力の計算は以下の式によって法線方向と接線方向の接触力（ＦｎとＦｔ）がそれぞれ計算される。
【数１】
0005800145.tif 000002
ここでｘとｖは相対変位および相対速度ベクトル、ｍｉｎは絶対値が小さい方の値を表す関数、ｎｔは接線方向単位ベクトル、ηはダッシュポットの粘性減衰係数であり、反発係数ｅと次式の関連がある。ここでｍｐは粒子質量である。
【数２】
0005800145.tif 000003
接触力には、この他に付着力や静電気力、ファンデルワールス力などを含めても良い。第１の力は、例えばベクトルによって示される。第１の力算出部１３は、解析の時間（時刻）毎（解析の時間（時刻）が数値積分時間増分Δｔ進められる毎）に第１の力を算出する。第１の力算出部１３は、変形体の表面に設けられる各粒子について算出した第１の力を示す情報を粒子位置更新部１５に出力する。
【００６０】
なお、接触候補リストの作成は、全ての解析上の時刻において行われるのではなく、予め設定された接触候補更新条件を満たした場合のみに行われればよい。接触候補更新条件としては、例えば、設定された回数の時刻が経過したことや、接触候補リストを作成してからの粒子の移動距離が設定された閾値を超えたことを条件とする。
【００６１】
また、第１の力算出部１３は、粒子間の接触だけでなく、変形体の表面に設けられた仮想的な粒子と壁面との間に生じる力や、流体から受ける抵抗力、電磁場から受ける力等についても算出してもよい。また、変形体が流体である場合には、表面張力を算出してもよい。また、上述したように第１の力は、変形体の表面に設けられる粒子を対象として個別要素法によって算出しえるものであるので、第１の力の算出には個別要素法に係る手法が適宜適用されてもよい。
【００６２】
第２の力算出部１４は、モデル生成部１２によって生成されたモデルにおいて、変形体を分割する四面体に生じる応力を算出し、当該四面体を分割する面を求めることによって、算出した応力から当該四面体を構成する各粒子に生じる第２の力を算出する第２の力算出手段である。具体的には、第２の力算出部１４は、四面体に生じる応力を、該四面体を構成する各粒子の位置の時間変化又は速度に基づいて、算出する。
【００６３】
第２の力算出部１４によって算出される第２の力は、変形体の変形に係る力である。上述したように第２の力は、変形体を分割する四面体を単位として計算される。なお、変形体の変形には、破壊が含まれていてもよい。例えば、四面体を構成する粒子間の予め定められた距離以上となった場合（粒子間の近傍関係が失われた場合）は、当該四面体は破壊されたものとしてもよい。破壊が生じた場合、変形体の表面に設けられたものでない（変形体の内部に位置していた）仮想的な粒子を、変形体の表面に設けられた粒子として扱うこととしてもよいし、新たに仮想的な粒子を配置して変形体の表面に新たな四面体を作製してもよい（破壊により変形体の内部の部分が露出したものとみなす）。
【００６４】
本実施形態では、第２の力を算出するための方法をＱＤＥＭ（QuadrapleDiscrete Element Method）と呼ぶ。ＱＤＥＭは、４体（４つの仮想的な粒子）間の相互作用のモデルを用いた手法である。構成則及び破壊基準は三次元である。現実の三次元物質のマクロな構成則との整合性がある。
【００６５】
上述したようにＱＤＥＭで用いられるモデルは、４つの仮想的な粒子の中心が十分に近い、当該中心を頂点とする四面体から構成されている。ＱＤＥＭでは、４粒子の相互作用から四面体に生じる１つの応力を決める。変形体が個体である場合、変形体に生じる応力とは、変形した四面体が元の形状に戻ろうとするために生じる応力である。これは弾性である。本実施形態においては、固体である場合の変形体は弾性体であるものとする。なお、四面体に生じる応力を算出することが可能であれば、変形体が弾性体以外のものであることを仮定してもよい。弾性体は、応力テンソルＴ（ｔ）が、現時刻での変形勾配テンソルＦ（ｔ）だけで定まる物質である。この場合、四面体に生じる応力は、弾性応力テンソルＴで表すことができる。
【００６６】
四面体に係る弾性応力テンソルＴは、以下のように算出することができる。まず、変形勾配テンソルＦが以下のように定義される。
ｄｘ＝Ｆ・ｄＸ
ここで、図１０に示すようにｄｘは現配置における、四面体の特定の頂点の座標から別の３つの頂点の座標へのベクトルから構成される行列である。ｄＸは、現配置よりも前の時刻の基準配置における、四面体の特定の頂点の座標から別の３つの頂点の座標へのベクトルから構成される行列である。四面体が微小ならば、ｄＸからｄｘへの変換は、線形変換である。Ｆの９つの成分を決めるためには３組の一次独立なベクトルの変換を調べればよい。３組の一次独立なベクトルは、上記のように四面体の４つの頂点の座標から算出できるものである。なお、基準配置と現配置との時間間隔は、上述した数値積分時間増分Δｔである。基準配置からΔｔ後の原配置は、後述するように基準配置に基づいて粒子位置更新部１５によって算出される。
【００６７】
第２の力算出部１４は、基準配置と現配置とから上記の式に基づいて、各四面体について変形勾配テンソルＦを算出する。続いて、第２の力算出部１４は、以下のＬｅｆｔ  Ｃａｕｃｈｙ−Ｇｒｅｅｎ変形テンソルＢを求める。
Ｂ＝Ｆ・ＦＴ
これは、Ｆを左極分解すると、Ｆ＝Ｖ・Ｒとなる。ここで、Ｖは正味の変形、Ｒは剛体回転を表している。従って、Ｂ＝Ｖ・Ｒ・ＲＴ・Ｖ＝Ｖ２となり、剛体回転の影響を受けない便利な量である。
【００６８】
ＢからＴへの関係を示す等方弾性体の構成式は、表示定理より以下の式により表される。
Ｔ＝φ０（ＩＢ，ＩＩＢ，ＩＩＩＢ）Ｉ＋φ１（ＩＢ，ＩＩＢ，ＩＩＩＢ）Ｂ＋φ２（ＩＢ，ＩＩＢ，ＩＩＩＢ）Ｂ２
ここで、ＩＢ，ＩＩＢ，ＩＩＩＢは、それぞれテンソルＢの３つの不変量である。また、Ｉは単位行列である。Ｃａｙｌｅｙ−Ｈａｍｉｌｔｏｎの定理を使って定式を変形して、Ｔが次式で与えられる。ここで、φ０、φ１、φ２を予め設定される物質定数の入力パラメータとする。
Ｔ＝（φ０−φ２ＩＩＢ）Ｉ＋（φ１＋φ２ＩＢ）Ｂ＋（φ２）ＩＩＩＢＢ−１
微小変形、微小歪みの仮定を使わないので、破壊の寸前まで扱うことができる。第２の力算出部１４は、上記の式を用いて四面体に係る弾性応力テンソルＴを算出する。
【００６９】
一方で変形体が粘弾性固体や流体である場合には、変形体に生じる応力とは、四面体内で速度勾配があるときに運動量の交換によって生じる応力である。これは粘性である。この場合、四面体に生じる応力は、粘性応力テンソルＴで表すことができる。四面体に係る粘性応力テンソルＴは、以下のように算出することができる。まず、速度勾配テンソルＬが以下のように定義される。
ｄｖ＝Ｌ・ｄｘ
ここで、図１１に示すようにｄｖは、現配置における四面体の特定の頂点における速度ベクトルｖ１を、別の３つの頂点における速度ベクトルｖ２，ｖ３，ｖ４から引いた３つの相対速度ベクトル（ｄｖａ＝ｖ２−ｖ１、ｄｖｂ＝ｖ３−ｖ１、ｄｖｃ＝ｖ４−ｖ１）から構成される行列である。四面体が微小ならば、ｄｘからｄｖへの変換は、線形変換である。Ｌの９つの成分を決めるためには３組の一次独立なベクトルの変換を調べればよい。３組の一次独立なベクトルは、上記のように四面体の４つの頂点の座標と相対速度とから算出できるものである。なお、現配置における四面体の特定の頂点における速度ベクトルは、後述するように粒子位置更新部１５によって算出される。
【００７０】
第２の力算出部１４は、基準配置と現配置とから上記の式に基づいて、各四面体について速度勾配テンソルＬを算出する。速度勾配テンソルＬは、以下のように加算分解することができる。
Ｌ＝Ｄ＋Ｗ
Ｄ＝１／２（Ｌ＋ＬＴ）  歪速度テンソル
Ｗ＝１／２（Ｌ−ＬＴ）  回転速度テンソル
【００７１】
ＤからＴへの関係を示す構成式は、表示定理より以下の式により表される。
Ｔ＝（−ｐ＋ψ０）Ｉ＋ψ１Ｄ＋ψ２Ｄ２
ここで、ｐは圧力項であり、予め設定された体積圧縮係数に四面体の体積変化率Ｊ＝ｄｅｔＦを乗じることによって得られる圧力増分δｐと、予め設定される初期圧力ｐ０から求める。ここで、ψ０、ψ１、ψ２を予め設定される物質定数（粘性係数）の入力パラメータとする。第２の力算出部１４は、上記の式を用いて四面体に係る弾性応力テンソルＴを算出する。
【００７２】
第２の力算出部１４は、各変形体の各四面体について、上述のようにして求めた応力（テンソル）を、当該四面体の４つの頂点を中心とする粒子に分配し、各粒子に生じる第２の力を算出する。具体的には、以下のように四面体を４つの等体積のブロックに分割して、分割面の法線ベクトルと面積とを算出する。まず、四面体の各辺の重心（２点の重心だから２重心と呼ぶ）の座標を算出する。四面体には辺が６本あるので６点ある。続いて、四面体の各面の重心（３点の重心だから３重心と呼ぶ）の座標を算出する。四面体には面が４個あるので４点ある。続いて、四面体の重心（４点の重心だから４重心と呼ぶ）の座標を算出する。
【００７３】
ここで、４重心と４つの３重心の中から選ばれた２つの３重心との３点からなる三角形は６個ある（４Ｃ２）。２重心とその２重心を含む２面の３重心とからなる三角形も６個ある（６Ｃ１）。これら１２個の三角形で四面体に切れ目を入れると、図１２に示すように４面体は等体積の４つのブロックに分割できる。
【００７４】
１つの粒子ｉに加わる力ｆｉは、四面体の応力テンソルＴ（弾性応力テンソルＴ、粘性応力テンソルＴ）と、四面体を４つのブロックに分割する三角形とから求められる。第２の力算出部１４は、四面体を４つのブロックに分割する全ての三角形の単位法線ベクトルｎと、当該三角形の面積ｄＳとを算出する。第２の力算出部１４は、四面体の応力テンソルＴと単位法線ベクトルｎとから、以下の式によって四面体内のある微小断面（四面体を４つのブロックに分割する三角形）上での応力ベクトルｔｎを算出する。
ｔｎ＝ＴＴ・ｎ
その微小断面の面積がｄＳであると、第２の力算出部１４は、面に働く力ｄｆを以下の式により算出する。
ｄｆ＝ｔｎｄＳ
１つの粒子ｉに加わる力ｆｉは、
ｆｉ＝Σｄｆ
である。ここで、Σは粒子ｉが頂点となっている全ての四面体から受ける力を足し合わせることを意味する。第２の力算出部１４は、上記の式から個々の粒子に加わる力を算出する。
【００７５】
上述したように、第２の力算出部１４は、四面体の各頂点の物理量と位置とから、その勾配テンソルを算出する。ここで、勾配テンソルは、弾性応力に対しては変形勾配テンソルであり、粘性応力に対しては速度勾配テンソルである。第２の力算出部１４は、算出した勾配テンソルを独立変数とするテンソル値関数である構成式を使って四面体に生じる応力を算出する。第２の力算出部１４は、応力テンソルと四面体内を分割する三角形の面の法線ベクトルから、当該三角形の面に生じる応力ベクトルを算出する。第２の力算出部１４は、算出した応力ベクトルと当該三角形の面の面積とから、四面体を構成する各粒子に生じる第２の力を算出する。
【００７６】
第２の力は、例えばベクトルによって示される。第２の力算出部１４は、解析の時間（時刻）毎（解析の時間（時刻）が数値積分時間増分Δｔ進められる毎）に第２の力を算出する。第２の力算出部１４は、変形体の表面に設けられる各粒子について算出した第２の力を示す情報を粒子位置更新部１５に出力する。
【００７７】
粒子位置更新部１５は、変形体の各粒子について、第１の力算出部１３によって算出された第１の力と、第２の力算出部１４によって算出された第２の力との和を取ることで各粒子に生じる力を算出して、算出した力に基づいて運動方程式を解くことで各粒子の位置を更新する粒子位置更新手段である。具体的には、第１の力と第２の力との和を取って算出された各粒子に生じる力をもとに、粒子の運動方程式（ｆ＝ｍａ、ｆ：各粒子に加わる力、ｍ：粒子の質量、ａ：粒子の速度）を数値積分して、数値積分時間増分Δｔだけ進められた時刻の粒子の位置座標及び速度を更新する。この際に、予め記憶部１１に記憶されている粒子の質量（密度から計算されてもよい）等の物理量が参照される。なお、第１の力は変形体の表面の粒子のみについて算出されているため、和が取られるのは変形体の表面の粒子についてのみである。変形体の表面以外の粒子については、第２の力のみが考慮される。
【００７８】
また、各粒子に加わる力として、第１の力及び第２の力だけでなく、予め設定された重力や外部からの力が粒子に加わることとして、粒子位置を計算することとしてもよい。計算された粒子の位置座標及び速度は、次の時刻（数値積分時間増分Δｔだけ進められた時刻）の粒子の位置座標及び速度として、記憶部１１に記憶され、次の時刻の計算に用いられる。
【００７９】
粒子位置更新部１５は、また、解析の終了条件を満たしているか否かを判断して、満たしていると判断した場合は、解析を終了させる。終了条件は、例えば、計算を行う時刻の回数が予め設定された回数実行されたことや、粒子の位置座標が収束している等の条件が用いられる。終了条件については、予め記憶部１１に記憶され、粒子位置更新部１５に参照される。
【００８０】
出力部１６は、解析装置１０による解析の結果を出力する手段である。例えば、出力部１６は、粒子位置更新部１５によって計算されて記憶部１１に記憶されている各時刻における粒子の位置及び速度を示す情報、即ち、変形体の状態を示す情報を出力する。出力は、例えば、解析装置１０の利用者が、解析の結果を確認できるように解析装置１０が備えるディスプレイに表示することによって行われる。あるいは、他の装置に解析結果の情報を出力することとしてもよい。以上が、解析装置１０の構成である。
【００８１】
引き続いて、図１３のフローチャートを用いて、本実施形態に係る解析装置１０で実行される処理（解析方法）を説明する。この処理は、解析装置１０の利用者が、解析装置１０に対して解析を開始する操作を行うことで開始される。
【００８２】
解析装置１０では、モデル生成部１２によって、上述したように、解析の対象となる変形体のモデルが生成される（Ｓ０１、モデル生成ステップ）。当該モデルは、三次元形状の変形体に、その中心が変形体を分割する四面体の頂点を構成する複数の仮想的な粒子が設けられたものである。生成されたモデルの情報は、記憶部１１に記憶され、必要に応じて各構成要素（第１の力算出部１３、第２の力算出部１４、粒子位置更新部１５及び出力部１６）から参照される。
【００８３】
続いて、第１の力算出部１３によって、モデル生成部１２によって生成されたモデルにおいて、第１の力の算出対象となる変形体の表面に設けられる粒子を対象として、接触候補リストが生成される（Ｓ０２、第１の力算出ステップ）。なお、隣接する粒子（四面体の辺の両端となっている粒子）同士は、接触候補とはならない。続いて、第１の力算出部１３によって、接触候補リストにおいて接触候補となっている粒子同士について接触しているかが判断され、接触している粒子同士についての接触力が計算される（Ｓ０３、第１の力算出ステップ）。続いて、第１の力算出部１３によって、当該接触力に基づいて個々の粒子に生じる第１の力が算出される（Ｓ０４、第１の力算出ステップ）。算出された第１の力を示す情報は、第１の力算出部１３から粒子位置更新部１５に出力される。
【００８４】
Ｓ０２の処理の後、一方では第２の力算出部１４によって、ＱＤＥＭにより、モデル生成部１２によって生成されたモデルにおいて変形体を分割する四面体に生じる応力が、当該四面体を構成する各粒子の位置の時間変化又は速度に基づいて算出される（Ｓ０５、第２の力算出ステップ）。続いて、第２の力算出部１４によって、当該四面体を分割する面が求められて、当該面が用いられて四面体に生じる応力から、個々の粒子に生じる第２の力が算出される（Ｓ０６、第２の力算出ステップ）。算出された第２の力を示す情報は、第２の力算出部１４から粒子位置更新部１５に出力される。なお、上述したＳ０３及びＳ０４の処理と、Ｓ０５及びＳ０６の処理とを並列に行うこととしていたが、必ずしも並列で行う必要はなく、何れか一方の処理を行う等順番に処理が行われてもよい。
【００８５】
第１の力及び第２の力が算出されると、続いて、粒子位置更新部１５によって、各粒子について第１の力と第２の力とが足し合わせられる（Ｓ０７、粒子位置更新ステップ）。続いて、粒子位置更新部１５によって、足し合わされた力に基づいて、粒子の運動方程式の数値積分が行われて、次の時刻（数値積分時間増分Δｔだけ進められた時刻）の各粒子の位置座標及び速度が計算されて更新される（Ｓ０８、粒子位置更新ステップ）。
【００８６】
続いて、粒子位置更新部１５によって、解析の終了条件を満たしているか否かが判断される（Ｓ０９）。ここで、終了条件とは、例えば、解析時間（シミュレーション時間）の上限値であったり粒子の座標変化に対する下限値であったりしてもよい。解析の終了条件を満たしていないと判断される（Ｓ０９のＮＯ）と、以下のように解析が継続される。まず、第１の力算出部１３によって、接触候補更新条件が満たされているか否かが判断される（Ｓ１０）。接触候補更新条件が満たされていると判断される（Ｓ１０のＹＥＳ）と、続いて、第１の力算出部１３によって、変形体の表面に設けられる粒子を対象として、再度、接触候補リストが生成される（Ｓ０２）。接触候補リストは、Ｓ０８によって更新された粒子の（最新の）位置座標に基づいて行われる。
【００８７】
Ｓ１０において接触候補更新条件が満たされていないと判断された場合（Ｓ１０のＮＯ）、あるいは、Ｓ０２において再度、接触候補リストが生成された場合、続いて、上述したＳ０３及びＳ０５以降の処理が繰り返される。なお、第１の力及び第２の力の算出は、Ｓ０８によって更新された粒子の（最新の）位置座標及び速度に基づいて行われる。
【００８８】
Ｓ０９において、解析の終了条件が満たされていると判断される（Ｓ０９のＹＥＳ）と、出力部１６によって解析の結果が出力されて、処理が終了する（Ｓ１１、出力ステップ）。なお、解析の結果の出力は、処理の最後のみに行われるのではなく、Ｓ０８の処理の後毎に各粒子の位置座標及び速度を示す情報が出力される等によって行われてもよい。また、Ｓ１１の出力と同時に記憶部１１に各粒子の位置座標及び速度を示す情報等が記憶されてもよい。以上が、本実施形態に係る解析装置１０で実行される処理（解析方法）である。
【００８９】
上述したように本実施形態では、個別要素法等の方法により、変形体間の衝突や接触に係る第１の力が算出され、また、ＱＤＥＭにより、変形体の変形に係る第２の力が算出され、それらの和の力に基づいて各粒子の位置が更新される。また、第１の力の算出と第２の力の算出との何れにおいても、変形体に設けられた仮想的な粒子が用いられる。従って、本実施形態では、第１の力の算出と第２の力の算出との両方が、上記の仮想的な粒子という共通の項目に対する計算にて行うことができる。これにより、本実施形態によれば、変形体の変形と変形体間の衝突とを別々の処理で計算する従来の方法と比べて、変形体間の複雑な接触判定や接触力の計算を変形体表面の四面体メッシュ頂点に配置された仮想的な表面粒子間の接触問題として扱うことで個別要素法等によって高速に行うことができ、更に、変形体の変形量や変形により生じる力は四面体メッシュの頂点である仮想的な粒子の座標や速度といった仮想的な粒子が持つ情報等から計算され、その力は変形体間の接触力を計算するときに使用した共通の仮想的な粒子に加わる力として処理される。従って、変形体間の衝突や接触と変形体自体の変形との両方を共通の仮想的な粒子の運動のみを計算することで処理することができるため、従来のＦＥＭとＤＥＭを組み合わせた手法等で見られる様な変形体の変形による力の計算と変形体間の接触力の計算との間で必要な情報のやり取りが少なく、変形体間の接触力を変形体の頂点に振り分ける処理も必要で無いため、変形体が多い場合でも、変形体間の衝突や接触と変形体自体の変形との両方を考慮して計算負荷が低く高速に解析を行うことができる。また、本手法を並列計算機を用いて並列計算を行う場合、従来のＦＥＭとＤＥＭを組み合わせた手法と比較して、ＦＥＭの計算に用いられる様な変形体の全体剛性マトリックス等の対象領域全てに関する巨大なマトリックスを作る必要が全くないので、並列計算で全ノード間通信の必要が全くない点で高速計算するにあたって有利である。
【００９０】
また、上述した粒子を成長させて変形体のモデルを生成することとすれば、多数の変形体に対して規則的で品質のよい四面体が設けられたモデルを効率的に生成することができる。具体的には変形体が多い場合、上述したポリゴンを用いる方法では、変形体の形状に合うように変形体に四面体メッシュを設けるのは計算時間がかかるが、本実施形態に係る粒子を成長させる方法によれば、ポリゴンを用いる方法に比べて変形体を分割する四面体を効率的に生成することができる。但し、変形体の形状が、ポリゴンを用いた場合でも効率的に四面体を生成できる形状の場合等には、必ずしも粒子を成長させて変形体のモデルを生成する必要はない。
【００９１】
また、解析に用いる粒子のパラメータとして、上述したように粒子の半径、質量及びバネ係数を用いることができる。この構成によれば、適切かつ確実に第１の力を算出することができる。
【００９２】
引き続いて、本実施形態に係る解析装置１０（解析方法）で実際に解析を行った例を示す。この例は、固定床面から０．５ｍの高さの空中に置かれたマットレスと同じ材料の球をマットレスと同時に落下させる解析である。この解析では、マットレスの床との衝突及び玉がマットレスに衝突することによる、マットレスの変形及び反発特性を同時に調べることができる。
【００９３】
この例では、以下のようにパラメータの値を与えた。密度：ρ＝２０００．０、弾性定数：（ヤング率Ｅ，ポアソン比ν）＝１．０×１０７，０．２５、減衰定数：ηｔ＝０．００１、粒子間法線方向バネ定数：Ｋｎ＝１．０×１０６、粒子間法線方向反発係数：ｅｎ＝０．８、粒子間接線方向バネ係数：Ｋｔ＝１．６×１０６、粒子間接線方向反発係数：ｅｔ＝０．８、粒子間摩擦係数：μ＝０．６、粒子−壁面境界間法線方向バネ定数：Ｋｎ＿ｗａｌｌ＝１．０×１０６、粒子−壁面境界間法線方向反発係数：ｅｎ＿ｗａｌｌ＝０．８、粒子−壁面境界間接線方向バネ係数：Ｋｔ＿ｗａｌｌ＝１．０×１０６、粒子−壁面境界間接線方向反発係数：ｅｔ＿ｗａｌｌ＝０．８、粒子−壁面境界間摩擦係数：μ＿ｗａｌｌ＝０．６、仮想粒子半径：ｒ＿ｖ＝０．０１、数値積分時間増分：Δｔ＝１．０×１０−５、計算領域ｘ方向最小値、最大値：ｘｍｉｎ，ｘｍａｘ＝０．０，５．０、計算領域ｙ方向最小値、最大値：ｙｍｉｎ，ｙｍａｘ＝０．０，５．０、計算領域ｚ方向最小値、最大値：ｚｍｉｎ，ｚｍａｘ＝０．０，５．０、重力加速度＝９．８。
【００９４】
図１４（ａ）に初期状態を示す。図の濃淡は、四面体に生じる歪みの絶対値を示す。図１４（ｂ）にマットレスと床面とが衝突した状態を示す。図１４（ｃ）に上昇するマットレスと球とが１度目に衝突した状態を示す。図１４（ｄ）にマットレスと球とが２度目に衝突した状態を示す。図１４（ｅ）にマットレスと球との３度目以降の衝突があり球が転がっている状態を示す。図１４（ｆ）に最終状態を示す。この例に示すように、本実施形態に係る解析装置１０（解析方法）では、変形体間の衝突と変形体の変形とが考慮された解析が実現できている。なお、図１４に示すように生成した解析の結果を示す情報を解析装置１０が備えるディスプレイに表示させてもよい。
【００９５】
引き続いて、上述した一連の解析装置１０の、１つ以上の変形体の運動の解析を行う処理をコンピュータに実行させるための解析プログラムを説明する。図１５に示すように、解析プログラム２１は、コンピュータに挿入されてアクセスされる、あるいはコンピュータが備える記録媒体２０に形成されたプログラム格納領域２０ａ内に格納される。
【００９６】
解析プログラム２１は、解析処理を統括的に制御するメインモジュール２１ａと、記憶モジュール２１ｂと、モデル生成モジュール２１ｃと、第１の力算出モジュール２１ｄと、第２の力算出モジュール２１ｅと、粒子位置更新モジュール２１ｆと、出力モジュール２１ｇとを備えて構成される。記憶モジュール２１ｂと、モデル生成モジュール２１ｃと、第１の力算出モジュール２１ｄと、第２の力算出モジュール２１ｅと、粒子位置更新モジュール２１ｆと、出力モジュール２１ｇとを実行させることにより実現される機能は、上述した解析装置１０の記憶部１１と、モデル生成部１２と、第１の力算出部１３と、第２の力算出部１４と、粒子位置更新部１５と、出力部１６との機能とそれぞれ同様である。
【００９７】
なお、解析プログラム２１は、その一部若しくは全部が、通信回線等の伝送媒体を介して伝送され、他の機器により受信されて記録（インストールを含む）される構成としてもよい。また、解析プログラム２１の各モジュールは、１つのコンピュータでなく、複数のコンピュータのいずれかにインストールされてもよい。その場合、当該複数のコンピュータによるコンピュータシステムよって上述した一連の解析プログラム２１の１つ以上の変形体の運動の解析を行う処理が行われる。また、解析プログラム２１をネットワークに接続されたサーバに格納しておき、所定のハードウェア装置にダウンロードしてもよい。
【産業上の利用可能性】
【００９８】
本発明の解析（シミュレーション）方法は、土木工学分野における地盤崩壊、建築構造物の破壊等、防災関連シミュレーション、化学工学分野における鉄鉱石コークスの輸送、荷下がり挙動解析、機械工学分野における建設機械の設計、ロボティクス、宇宙インフラの開発、海洋工学分野における海底ケーブル、通信ケーブルの耐久性評価、製薬分野における充填時のカプセルの挙動、耐久性評価、食品加工分野における粒状体（米、菓子）の輸送、充填等、変形できる多面体要素の運動解析に応用が期待される。
【符号の説明】
【００９９】
１０…解析装置、１１…記憶部、１２…モデル生成部、１３…第１の力算出部、１４…第２の力算出部、１５…粒子位置更新部、１６…出力部、２０…記録媒体、２０ａ…プログラム格納領域、２１…解析プログラム、２１ａ…メインモジュール、２１ｂ…記憶モジュール、２１ｃ…モデル生成モジュール、２１ｄ…第１の力算出モジュール、２１ｅ…第２の力算出モジュール、２１ｆ…粒子位置更新モジュール、２１ｇ…出力モジュール。
【図１】
0005800145.tif 000004
【図２】
0005800145.tif 000005
【図３】
0005800145.tif 000006
【図４】
0005800145.tif 000007
【図５】
0005800145.tif 000008
【図６】
0005800145.tif 000009
【図７】
0005800145.tif 000010
【図８】
0005800145.tif 000011
【図１０】
0005800145.tif 000012
【図１１】
0005800145.tif 000013
【図１２】
0005800145.tif 000014
【図１３】
0005800145.tif 000015
【図１５】
0005800145.tif 000016
【図９】
0005800145.tif 000017
【図１４】
0005800145.tif 000018