(57)【要約】      （修正有）【課題】第１構造体と第２構造体とが接触して第１構造体が変形し、第１構造体と第２構造体との接触面が変化する際の、第１構造体の変形挙動を、比較的短時間で、かつ高精度に求める。【解決手段】第１構造体を複数の粒子の集合体として再現した第１モデル、および、前記第２構造体を再現した第２モデルを生成し、第１モデルに対して第２モデルを移動させることで、第２モデルと接触した第１モデルを時間をかけて変形させたときの、第１モデルの変形挙動を計算する際、第２モデルの物理量は第１モデルから何らの作用も受けないものとし、第１モデルに対して第２モデルを所定の時間ステップ毎に移動させて、移動の度に、第１モデルを構成する複数の粒子のうち、第２モデルとの距離が予め定められた所定値以下になった粒子の速度情報を、所定値以下の範囲にある第２モデルの速度情報と同一として、各時間ステップ毎の第１モデルの変形を計算する。【選択図】図４
2008152423.tif 000002
【特許請求の範囲】
【請求項１】
第１構造体と第２構造体とが接触して前記第１構造体が時間をかけて変形し、前記第１構造体と前記第２構造体との接触面が時間をかけて変化する際の、前記第１構造体の変形挙動を求める方法であって、
前記第１構造体を複数の粒子の集合体として再現した第１モデル、および、前記第２構造体を再現した第２モデルを生成するモデル生成ステップと、
前記第１モデルに対して前記第２モデルを移動させることで、前記第２モデルと接触した前記第１モデルを時間をかけて変形させたときの、前記第１モデルの変形挙動を計算するステップと、
前記計算結果を出力するステップと、を有し
前記変形挙動を計算するステップでは、前記第２モデルの物理量は前記第１モデルから何らの作用も受けないものとし、前記第１モデルに対して前記第２モデルを所定の時間ステップ毎に移動させて、移動の度に、前記第１モデルを構成する複数の粒子のうち、前記第２モデルとの距離が予め定められた所定値以下になった粒子の速度情報を、前記所定値以下の範囲にある前記第２モデルの速度情報と同一として、各時間ステップ毎の前記第１モデルの変形を計算することを特徴とする、粒子モデルを用いた変形挙動シミュレーション方法。
【請求項２】
前記第２モデルは、前記第２構造体を複数の粒子の集合体として再現したモデルであることを特徴とする請求項１記載の粒子モデルを用いた変形挙動シミュレーション方法。
【請求項３】
前記第２モデルは、前記第２構造体のうち、前記第１構造体の任意の粒子との距離が前記所定値以下になる可能性のある領域のみを再現したモデルであることを特徴とする請求項１または２記載の粒子モデルを用いた変形挙動シミュレーション方法。
【請求項４】
前記変形挙動を計算するステップでは、前記粒子のそれぞれについて設定されている、前記粒子の移動を規定する粒子移動規定条件に基づいて前記変形挙動を計算し、
前記粒子移動規定条件は、前記粒子のそれぞれに対応して近傍領域が定められるとともに、この近傍領域内の注目物理量の分布をこの近傍領域に含まれる粒子の持つ注目物理量を用いて近似することによって定められる条件であることを特徴とする請求項１〜３のいずれかに記載の粒子モデルを用いた変形挙動シミュレーション方法。
【請求項５】
前記粒子移動規定条件は、前記粒子によって構成されたモデルにおける近似したい物理量をｘ、前記近傍領域内の注目物理量の分布をｆ（ｘ）として、この注目物理量の分布ｆ（ｘ）を前記近傍領域に含まれる粒子の持つ注目物理量と所定のカーネル関数とを用いて下記式（１）で定まる近似物理量の分布＜ｆ（ｘ）＞で近似することによって定式化される方程式である請求項４に記載の粒子モデルを用いた変形挙動シミュレーション方法。
【数１】
2008152423.tif 000019
ただし、Ｗ（ｘ−ｘ’，ｈ）はカーネル関数、ｈは影響半径。
【請求項６】
第１構造体と第２構造体とが接触して前記第１構造体が時間をかけて変形し、前記第１構造体と前記第２構造体との接触面が時間をかけて変化する際の、前記第１構造体の変形挙動を求める方法をコンピュータに実行させるプログラムであって、
前記第１構造体を複数の粒子の集合体として再現した第１モデル、および、前記第２構造体を再現した第２モデルを生成させるモデル生成ステップと、
前記第１モデルに対して前記第２モデルを移動させることで、前記第２モデルと接触した前記第１モデルを時間をかけて変形させたときの、前記第１モデルの変形挙動を計算させるステップと、
前記計算結果を出力させるステップと、を有し、
前記変形挙動を計算させるステップでは、前記第２モデルの物理量は前記第１モデルから何らの作用も受けないものとし、前記第１モデルに対して前記第２モデルを所定の時間ステップ毎に移動させて、移動の度に、前記第１モデルを構成する複数の粒子のうち、前記第２モデルとの距離が予め定められた所定値以下になった粒子の速度情報を、前記所定値以下の範囲にある前記第２モデルの速度情報と同一として、各時間ステップ毎の前記第１モデルの変形を計算させることを特徴とする粒子モデルを用いた変形挙動シミュレーションプログラム。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、第１構造体と第２構造体とが接触して第１構造体が時間をかけて変形し、第１構造体と第２構造体との接触面が時間をかけて変化する際の、第１構造体の変形挙動を求める方法であって、特に、ゴムなどの超弾性体と剛体とが接触して超弾性体が変形し、超弾性体と剛体との接触面積が時間をかけて５％以上変化するような、いわゆる大変形問題における第１構造体の変形挙動を再現する方法に関する。
【背景技術】
【０００２】
ＳＰＨ（ｓｍｏｏｔｈｅｄ  Ｐａｒｔｉｃｌｅ  Ｈｙｄｒｏｄｙｎａｍｉｃｓ）法は、有限要素法とは異なり、点の集合のみで物質を表すメッシュ非依存の解析法である。ＳＰＨ法は、従来から広く利用されてきたＦＥＭ法と比較して、・メッシュレス法なので、要素作成が困難な複雑な３次元構造物でも容易にモデル化できる。・解析に応じた要素タイプの選別が不要。・精度を意識した要素分割の経験・技術が不要。・要素のつぶれによる解析不能が生じる大歪、大変形問題にも、柔軟な対応が期待できる。・構造解析、熱解析、衝突解析、流体解析、電磁場解析、等々で共通の粒子モデルが使用できる。等の長所がある。このため、ＦＥＭが主流となっている弾塑性解析や超弾性解析等について、ＳＰＨ法の応用が期待されている。
【０００３】
例えば、下記非特許文献１には、ＳＰＨ解析手法を用いた、固体の大変形弾塑性解析手法が提案されている。このような固体の大変形弾塑性解析手法に対し、本願発明者は、下記非特許文献２において、ＳＰＨ法を用いた、ゴム材料等の超弾性大変形問題の解析手法を提案している。非特許文献２では、複数の粒子の集合として表されたゴム構造体モデルに、同じく複数の粒子の集合として表された剛体プレートモデルを押し付けて、ゴム構造体モデルを準静的に変形させた際の、このゴム構造体モデル（の再現するゴム構造体）の変形挙動を計算している。
【０００４】
【非特許文献１】宋武、酒井譲、山下彰彦“ＳＰＨ法による弾塑性解析手法の検討（第１報  ２次元小変形問題）”、日本機械学会論文集（Ａ編）６８巻６６９号、ｐ７７２−７７８、２００２
【非特許文献２】河島庸一、酒井譲“粒子法によるゴムの超弾性解析”計算工学講演会論文集  第１１巻  第１号  Ｐ５３〜  ２００６年６月
【発明の開示】
【発明が解決しようとする課題】
【０００５】
非特許文献１や非特許文献２などに記載されている従来の方法では、ＳＰＨ法など点情報で表される２つの構造体（流体などの連続体を含む）の衝突など、２つの構造体の接触に基づく弾性体の変形を解析する場合、各構造体を同様に、粒子の集合体モデルで表して、各構造体のモデルの粒子同士の相互作用を計算していた。
【０００６】
例えば上記非特許文献２に記載されている例においても、ゴム構造体モデルと剛体プレートモデルとの特に接触面部分において、各構造体モデルを構成する複数の粒子同士の相互作用を計算している。しかし非特許文献２に記載されているような、ゴム構造体の準静的な変形では、剛体プレートが受ける相互作用については微々たるものであり、接触部分で剛体プレートは変形していないとみなせることは周知のとおり（前提）である。しかし、従来は、このように、ゴム構造体モデルの準静的な変形のみを計算したい場合などでも、ゴム構成部材の変形挙動とともに、実際に有意な解を得られないとわかっている、このような剛体プレートモデルの接触部分が受ける相互作用まで計算している。このため、粒子の集合体モデルを用いた従来の変形挙動シミュレーション方法では、計算に余分な時間やコストがかかってしまうといった問題点があった。
【０００７】
そこで、本発明は、第１構造体と第２構造体とが接触して第１構造体が時間をかけて変形し、第１構造体と第２構造体との接触面が時間をかけて変化する際の、第１構造体の変形挙動を、比較的短時間で、かつ高精度に求めることができる、粒子モデルを用いた変形挙動シミュレーション方法を提供することを目的とする。
【課題を解決するための手段】
【０００８】
上記目的を達成するために、本発明は、第１構造体と第２構造体とが接触して前記第１構造体が時間をかけて変形し、前記第１構造体と前記第２構造体との接触面が時間をかけて変化する際の、前記第１構造体の変形挙動を求める方法であって、前記第１構造体を複数の粒子の集合体として再現した第１モデル、および、前記第２構造体を再現した第２モデルを生成するモデル生成ステップと、前記第１モデルに対して前記第２モデルを移動させることで、前記第２モデルと接触した前記第１モデルを時間をかけて変形させたときの、前記第１モデルの変形挙動を計算するステップと、前記計算結果を出力するステップと、を有し、前記変形挙動を計算するステップでは、前記第２モデルの物理量は前記第１モデルから何らの作用も受けないものとし、前記第１モデルに対して前記第２モデルを所定の時間ステップ毎に移動させて、移動の度に、前記第１モデルを構成する複数の粒子のうち、前記第２モデルとの距離が予め定められた所定値以下になった粒子の速度情報を、前記所定値以下の範囲にある前記第２モデルの速度情報と同一として、各時間ステップ毎の前記第１モデルの変形を計算することを特徴とする、粒子モデルを用いた変形挙動シミュレーション方法を提案する。
【０００９】
なお、前記第２モデルは、前記第２構造体を複数の粒子の集合体として再現したモデルであることが好ましい。また、前記第２モデルは、前記第２構造体のうち、前記第１構造体の任意の粒子との距離が前記所定値以下になる可能性のある領域のみを再現したモデルであることが好ましい。
【００１０】
また、前記変形挙動を計算するステップでは、前記粒子のそれぞれについて設定されている、前記粒子の移動を規定する粒子移動規定条件に基づいて前記変形挙動を計算し、前記粒子移動規定条件は、前記粒子のそれぞれに対応して近傍領域が定められるとともに、この近傍領域内の注目物理量の分布をこの近傍領域に含まれる粒子の持つ注目物理量を用いて近似することによって定められる条件であることが好ましい。
【００１１】
また、前記粒子移動規定条件は、前記粒子によって構成されたモデルにおける近似したい物理量をｘ、前記近傍領域内の注目物理量の分布をｆ（ｘ）として、この注目物理量の分布ｆ（ｘ）を前記近傍領域に含まれる粒子の持つ注目物理量と所定のカーネル関数とを用いて下記式（１）で定まる近似物理量の分布＜ｆ（ｘ）＞で近似することによって定式化される方程式であることが好ましい。
【数２】
2008152423.tif 000003
ただし、Ｗ（ｘ−ｘ’，ｈ）はカーネル関数、ｈは影響半径。
【００１２】
本発明は、また、第１構造体と第２構造体とが接触して前記第１構造体が時間をかけて変形し、前記第１構造体と前記第２構造体との接触面が時間をかけて変化する際の、前記第１構造体の変形挙動を求める方法をコンピュータに実行させるプログラムであって、前記第１構造体を複数の粒子の集合体として再現した第１モデル、および、前記第２構造体を再現した第２モデルを生成させるモデル生成ステップと、前記第１モデルに対して前記第２モデルを移動させることで、前記第２モデルと接触した前記第１モデルを時間をかけて変形させたときの、前記第１モデルの変形挙動を計算させるステップと、前記計算結果を出力させるステップと、を有し、前記変形挙動を計算させるステップでは、前記第２モデルの物理量は前記第１モデルから何らの作用も受けないものとし、前記第１モデルに対して前記第２モデルを所定の時間ステップ毎に移動させて、移動の度に、前記第１モデルを構成する複数の粒子のうち、前記第２モデルとの距離が予め定められた所定値以下になった粒子の速度情報を、前記所定値以下の範囲にある前記第２モデルの速度情報と同一として、各時間ステップ毎の前記第１モデルの変形を計算させることを特徴とする粒子モデルを用いた変形挙動シミュレーションプログラムも、併せて提供する。
【発明の効果】
【００１３】
本発明の粒子モデルを用いた変形挙動シミュレーション方法、およびプログラムによれば、第１構造体と第２構造体とが接触して第１構造体が時間をかけて変形し、第１構造体と第２構造体との接触面が時間をかけて変化する際の、第１構造体の変形挙動を、比較的短時間で、かつ高精度に求めることができる。また、第２モデルを複数の粒子の集合体として再現することで、第１モデルと同様のモデルとして取り扱うことができ、シミュレーションにおける各種の数値計算処理が、同一のアプリケーションプログラムによってスムーズに実行されることができる。また、第２モデルの接触面部分の粒子数を少なくすることができ、さらに余分な相互作用を考慮する必要がないため、計算処理に係る時間を極端に少なくすることができる。
【発明を実施するための最良の形態】
【００１４】
以下、本発明の粒子モデルを用いた変形挙動シミュレーション方法およびプログラムについて、添付の図面に示される好適実施例を基に詳細に説明する。図１は、本発明の粒子モデルを用いた変形挙動シミュレーション方法を実施する変形挙動シミュレーション装置１０（装置１０）の概略の構成を示すブロック図である。装置１０では、例えばゴムなどの超弾性体である第１構造体と、例えば剛体である第２構造体とが接触して第１構造体が時間をかけて変形し、第１構造体と前記第２構造体との接触面が時間をかけて変化する際の、第１構造体の変形挙動を求める。装置１０は、モデル作成部１２、移動・圧縮演算部１４、変形演算部１６、メモリ１８、およびＣＰＵ（制御部）２０を有して構成される本体装置２２と、本体装置２２に接続されるマウスおよびキーボードからなる入力操作系２４と、本体装置２２に接続されるプリンタやディスプレイなどの出力部２６とを備えて構成されている。装置１０は、メモリ１８に記憶されたプログラムをＣＰＵ２０が実行することで上記各部位が機能するコンピュータである。
【００１５】
モデル作成部１２は、入力操作系２４を用いて入力されたモデル作成条件に応じて、第１構造体を複数の粒子の集合体として再現した第１モデル、および、第２構造体を再現した第２モデルを作成する。図２は、本発明の変形挙動シミュレーション方法によって作成されるモデルの一例を示している。図２に示す第１モデル３２は、超弾性体であるゴム体を、複数の粒子の集合で再現したモデルである。また、図２に示す第２モデル３４は、ゴム体と比較してはるかに硬度が高く、剛体として取り扱うことのできる金属体を、同じく複数の粒子の集合で再現したモデルである。以下、第２モデル３４が表す金属体平板によって、第１モデル３２が表すゴム体を十分に遅い速度で圧縮する際の、このゴム体の変形挙動を再現する実施形態について説明する。このような本実施形態における変形挙動は準静的な変形挙動であり、ゴム体と金属体との接触面積が、時間の経過とともに、初期状態と比べて５％以上も増減するような大変形問題であるといえる。
【００１６】
本実施形態では、このように、少なくとも第１の構造体を、複数の粒子で再現し、変形中の超弾性体のモデルを表現する、上記ＳＰＨ法を用いたシミュレーションを実施する。ＳＰＨ法では、従来の方法のように、空間格子や有限要素が設けられない。しかも、第１モデルの粒子は、粒子移動規定条件の下に自由に移動するので、大きな変形を超弾性体などの第１構造体に固定した座標系で記述するLagrangianに基づいて表現しても、数値計算上の適切な解を得ることができる。
【００１７】
粒子を用いて表す連続体の挙動は、公知のように、いわゆる連続の式、運動方程式およびエネルギー方程式の各式を支配方程式として用いて表される。超弾性体における各部分の速度（ゴムの流速）等の注目物理量の分布ｆ（ｘ）（ｘは近似したい物理量であり、本実施形態では、空間の３次元位置座標を表す）は、下記式（１）に示すように、所定の位置座標ｘ’における注目物理量ｆ（ｘ’）とカーネル関数Ｗ（ｘ−ｘ’，ｈ）（ｈは影響半径パラメータ）とを用いて求められる近似物理量の分布＜ｆ（ｘ）＞によって近似することができる。
【数３】
2008152423.tif 000004
【００１８】
粒子の持つ注目物理量ｆ（ｘ）は、粒子Ｊの密度ρＪと、粒子Ｊの質量ｍＪとを用いて（Ｊは粒子番号を表す）、下記式（２）のように離散化して表すことができる。
【数４】
2008152423.tif 000005
【００１９】
式（２）の微分系は、下記式（３）のように表され、粒子の運動方程式は、応力ダイバージェンスを用いて下記式（４）のように表される。式（４）におけるσは、Ｃａｕｃｈｙ応力テンソルである。
【数５】
2008152423.tif 000006
【数６】
2008152423.tif 000007
【００２０】
上記式（３），（４）より、Ｋｅｒｎｅｌ近似を用いた粒子の運動方程式は、式(５)で表され、注目する粒子Ｉの粒子加速度ａＩが求まる。Ｉは、複数の粒子のうち注目する粒子の番号を指す。以降注目する粒子をＩとする。
【００２１】
このように、速度等の注目物理量の分布ｆ（ｘ）を粒子Ｊの持つ注目物理量を用いて近似することによって定められる式（５）が作成される。粒子Ｉの位置と速度は、下記式（６）、（７）で求まる。本実施形態では、式（５）、（６）、（７）を、粒子Ｉの移動を規定する粒子移動規定条件として備えた第１モデルを作成する。
【数７】
2008152423.tif 000008
【数８】
2008152423.tif 000009
【数９】
2008152423.tif 000010
【００２２】
ここで、カーネル関数Ｗ（ｘ−ｘ’，ｈ）は、特に限定されないが、｜ｘ−ｘ’｜／ｈ＝νとして、例えば、下記式（８）に示すように、３次のＢ−スプライン関数を用いるとよい。
【数１０】
2008152423.tif 000011
【００２３】
式（８）に示す３次のＢ−スプライン関数をカーネル関数Ｗ（ｘ−ｘ’，ｈ）として用いると、νが２以上の場合、カーネル関数Ｗ（ｘ−ｘ’，ｈ）は０となり、νが２未満の場合正の値となるため、注目する粒子Ｉの挙動は、νが２未満の範囲内に位置する粒子Ｊに依存することになる。例えば、図２に示すように、粒子Ｍを中心とする半径２ｈの球の範囲（この範囲を粒子Ｍに対応した近傍領域ＲＭという）内に含まれる粒子Ｌによって、粒子Ｍの移動が規定される。この半径２ｈはスムージングレングスともいう。本発明では、このスムージングレングスとともに、図２に示すように、第１モデル３２の各粒子に対して、接触条件範囲ＣＮ（粒子Ｍについては、接触条件範囲ＣＭ）が設定される。接触条件範囲ＣＮとは、この範囲内に第２モデルが入った場合は、第２モデルと第１モデルとが接触したとみなす範囲である。本実施形態では、接触条件範囲Ｃを影響半径ｈの大きさに設定する。これら、スムージングレングスや接触条件範囲は、本実施形態の例に限定されず、任意の数値範囲を設定して構わない。本発明では、複数の粒子で構成された、このような第１モデルを作成する。
【００２４】
本実施形態では、図２に示すように、第１モデル３２と同様、第２モデル３４も、粒子の集合体として作成している。第２モデル３４を構成する粒子については、第２モデルから何らの作用も受けないよう設定されている。また、第１モデル３２の任意の粒子Ｋ（番号「Ｋ」については特に図示せず）の近傍領域ＲＫ内に、第２モデル３４の粒子が入ったとしても、第１モデルの任意の粒子Ｋは、この第２モデル３４の粒子によっては、その移動は規定されない。一方、第１モデル３２の任意の粒子Ｋの接触条件範囲ＣＫ内に、第２モデル３４の粒子が入った場合は、第１モデルの任意の粒子Ｋの速度ベクトルが、この第２モデル３４の粒子の速度ベクトルと同一のベクトルとなるよう、設定されている。この場合も、接触条件範囲ＣＫ内に入った第２モデル３４の粒子は、第１モデルの粒子から何らの作用を受けず、速度ベクトルも一切変らないよう設定されている。
【００２５】
本発明では、第２モデルについては特に限定されず、例えば、有限要素によってモデル化したり、接触面を表す関数によってモデル化しても構わない。ただし、第２モデルを複数の粒子の集合体として再現することで、第１モデルと同様のモデルとして取り扱うことができ、後述する各種の数値計算処理が、同一のアプリケーションプログラムによってスムーズに実行されることができる。
【００２６】
なお、図２に示す例では、第２モデルについて、金属体のうち、ゴム体と接触する可能性のある表面の一部分について再現しており、表面部分を構成する粒子は３層となっている。本発明において、第２モデルは、第２構造体のうち、第１構造体と接触する可能性のある領域のみを再現するモデルであればよい。図３（ａ)は、従来の変形挙動シミュレーション方法で作成される第１モデルおよび第２モデルの一例であり、図３（ｂ）は、本発明の変形挙動シミュレーション方法で作成される第１モデルおよび第２モデルの他の例である。図３（ａ）に示すモデルは、第２モデル３４ａを構成する粒子についても、第１モデル３２ａと同様に粒子移動規定条件を設定しておき（金属体と見なせるような条件）、第１モデルと第２モデルとの相互作用も計算しながら第２モデルの変形を計算するためのモデルである。このように従来では、第２モデルの粒子についても、移動規定条件や相互作用を考慮して計算する必要があり、第２モデルには、第１モデルと接触する領域以外にも、複数の粒子が必要であり（最低２個必要）、各粒子についての計算にかかる時間は多大なものとなっていた。本願発明では、第２モデルとしては、図３（ｂ）のように、第１構造体と接触する可能性のある領域のみを再現するモデル（第２モデル３４ｂ）であればよい。本願発明では、接触条件範囲ＣＮ内に入った第２モデルの粒子は、第１モデルの粒子から何らの作用を受けず、速度ベクトルも一切変らないよう設定されており、接触条件範囲ＣＮ内に入る可能性のある第２モデルの粒子のみをモデル化しておけば、第１モデルと接触することで第２モデルの表面形状は一切変らず、第２モデルを剛体モデルとして取り扱うことができる。図３（ａ）に示す第２モデル３４ａの粒子数は２３０８個であり、図３（ｂ）に示す第２モデル３４ｂの粒子数は２０２個である。図３（ａ）に示す場合に比べ、図３（ｂ）に示す場合の方が粒子数も単純に少なく、かつ余分な相互作用を考慮する必要がないため、計算処理に係る時間は極端に少ない。
【００２７】
移動・圧縮演算部１４は、入力操作系２４を用いて入力された圧縮条件に基づいて、第１モデル３２に対して第２モデルを移動させるシミュレーション演算を実施する。本実施形態では、第１モデル３２に対して第２モデル３４を一定速度で近づけていくシミュレーション演算を行う。移動・圧縮演算部１４では、単位時間間隔（時間ステップ）毎に、第１モデル３２に第２モデル３４を近づけていき、各時間間隔毎に、第１モデル３２の全粒子のうち、第２モデルと接触したとみなす粒子、すなわち接触条件範囲内に第２モデルが入った粒子を抽出する。そして、第２モデルと接触したとみなす粒子については、接触した段階における第２モデルの速度ベクトルを、この粒子の速度ベクトルとして付与する。この速度ベクトル（速度情報）の付与は、接触以前の粒子の位置や速度などの物理量には全く依存させることなく、無条件で強制的に行う。
【００２８】
具体的には、第１モデル３２の各粒子のうち、上記接触条件範囲Ｃ以内に第２モデルが入った粒子の速度ベクトルを、第２モデルの各粒子の速度ベクトルに設定する。この際、第２モデルの粒子について、第１モデルと衝突することによる一切の相互作用は無いものとみなし、第２モデルの速度ベクトルについては、接触・圧縮に起因した何らの変更も加えない。本実施形態のように準静的な変形では、接触・圧縮における第２構造体の移動速度は変化しないとみなすことができるので、このように速度ベクトルの強制付与によって、第１構造体と第２構造体の接触・圧縮状態を充分な精度で再現できる。本願発明では、第２構造体にかかる余分な相互作用を考慮する必要がないため、計算処理に係る時間は極端に少ない。また、同時に、第２モデルの起因する計算ノイズを低減することができる。なお、移動・圧縮演算部１４で用いる第１モデル３２は、１つ前の時間ステップにおいて、後述する変形演算部１６において計算された、１つ前の時間ステップ経過時点の圧縮変形状態にあるモデルである。なお、第１の時間ステップにおいて用いる第１モデル３２は、初期状態、例えば複数の粒子が等間隔に配列されたモデルを用いる。
【００２９】
変形演算部１６は、移動・圧縮演算部１４において第１モデル３２の接触部分の粒子に付与された速度ベクトルを外力条件として、現在の時間ステップｎにおける、第１モデルの圧縮変形状態を計算する。本実施形態では、非圧縮条件を体積一定の条件で解くアルゴリズムを適用した、非圧縮粘性流体解析法を用いる。非圧縮粘性流体解析法を用いるのは、ゴム等の高分子材料は、圧力によっても体積がほとんど変化しない非圧縮あるいは微圧縮性を示すが、ＳＰＨ法は基本的には圧縮性の取り扱いであり、そのままでは厳密には非圧縮状態を計算できないからである。なお、本発明の粒子モデルを用いた変形挙動シミュレーションは、圧縮条件および非圧縮条件のいずれの条件でも適用可能であり、変形挙動シミュレーションにおけるアルゴリズムの詳細については、特に限定されない。
【００３０】
まず、１つ前の時間ステップ（ｎ−１）経過時点の粒子速度と圧力場とを用いて、現時間ステップｎにおける仮の速度ｕｉ＊（下記式（９））を求める。
【数１１】
2008152423.tif 000012
ここに、ｕは粒子速度、Ｐは圧力、Ｖは動粘性係数、ｆは体積力である。このままでは、仮の粒子速度は連続の式を満たさないので、式（１０）に示すように、非回転速度場ｕ’を用いて、ｎステップ目の速度を表す。また、ｎステップ目の圧力場は、式（１１）のようになる。式（９）（１０）（１１）を用いて、ｎステップにおける粒子速度ｕｎを求めると、下記式（１２）のようになる。
【数１２】
2008152423.tif 000013
【数１３】
2008152423.tif 000014
【数１４】
2008152423.tif 000015
変形演算部１６では、このような式（１２）によって、仮の粒子速度を陽解法的に求める(外力による変形サブステップ〈ｎａ〉）。
【００３１】
ここで、ｕ＊は連続の式を満たさなくてはならないので、式（１２）の発散をとると、式（１３）の圧力に関するポアソン方程式が得られる。
【数１５】
2008152423.tif 000016
【００３２】
このポアソン方程式の左辺を離散化すると式（１４）となる。ここで、Ｐａｂ＝Ｐａ−Ｐｂ、ｒａｂ＝ｒａ−ｒｂである。ポアソン方程式の右辺は、式（１５）のように離散化される。
【数１６】
2008152423.tif 000017
【数１７】
2008152423.tif 000018
【００３３】
変形演算部１６では、外力による変形サブステップ〈ｎａ〉で求めた仮の粒子速度によって粒子を移動させた後、非圧縮性条件を満たすように式（１３）の圧力ポアソン方程式を陰解法的に解いて、修正された圧力と速度を求め、粒子を移動させる（位置修正サブステップ〈ｎｂ〉）。本実施形態では、このように、１つの時間ステップにおいて、外力による変形サブステップ〈ｎａ〉と位置修正サブステップ〈ｎｂ〉とを実施して、時間ステップｎにおける粒子の位置を求める。
【００３４】
また、本実施形態では、変形演算部１６において、変形サブステップ〈ｎａ〉、および位置修正サブステップ〈ｎｂ〉の各サブステップ毎に変形勾配を求め、時間ステップｎにおける各粒子の応力を求める。本実施形態では、Ｕｐｄａｔｅ  Ｌａｇｒａｎｇｉａｎ法を用いて、時間ステップｎにおける変形勾配テンソルＦｉｊを求める。Ｕｐｄａｔｅ  Ｌａｇｒａｎｇｉａｎ法に基く変形勾配テンソルＦは、１つ前のステップ＜ｎ−１＞までにアップデートされた変形勾配テンソル、すなわち、初期状態からステップ＜ｎ−１＞が終了するまでの変形勾配をＦ＜ｎ−１＞とすると、現在のステップでの変形勾配テンソルをｆ＜ｎ＞として、現ステップ＜ｎ＞終了時の変形勾配テンソルＦ＜ｎ＞は以下の式（１６）で表される
Ｆ＜ｎ＞＝ｆ＜ｎ＞Ｆ＜ｎ−１＞    ・・・（１６）
【００３５】
本実施形態では１つのステップは、上述のように、ＳＰＨ法に基く外力による変形サブステップ＜ｎａ＞と、粒子密度一定のための位置修正サブステップ＜ｎｂ＞の、２つのサブステップに分かれるので、さらに、各々のサブステップ毎に変形勾配を求め、応力を計算する。具体的には、ステップ＜ｎ−１＞終了時の変形勾配テンソルをＦ＜ｎ−１＞として、変形サブステップ＜ｎａ＞における変形勾配テンソルをｆａ＜ｎ＞として、このサブステップ終了時までの変形勾配テンソルＦａ＜ｎ＞を以下の式（１７）で求める。
Ｆａ＜ｎ＞＝ｆａ＜ｎ＞Ｆ＜ｎ−１＞    ・・・（１７）
【００３６】
そして、次に、粒子密度一定のための位置修正サブステップ＜ｎｂ＞における変形勾配テンソルをｆｂ＜ｎ＞として、このサブステップ終了時までの変形勾配テンソルＦｂ＜ｎ＞を、現ステップ＜ｎ＞終了時の変形勾配テンソルＦ＜ｎ＞として、以下の式（１８）により求める。
Ｆｂ＜ｎ＞＝ｆｂ＜ｎ＞Ｆａ＜ｎ＞＝Ｆ＜ｎ＞・・・（１８）
【００３７】
このようにして求められた変形勾配テンソルを元に、各ステップ終了時の各粒子の応力を求める。例えば、弾性ポテンシャル関数としてＭｏｏｎｅｙ−Ｒｉｖｌｉｎ体を用いて表した、Ｃａｕｃｈｙ応力テンソルを求めればよい。変形演算部１６では、このようにして、現在の時間ステップｎにおける第１モデルの圧縮変形状態、例えば、第１モデルの各粒子の位置および応力を計算する。Ｍｏｏｎｅｙ−Ｒｉｖｌｉｎ体を用いて表した、Ｃａｕｃｈｙ応力テンソルの導出や、上記各処理の詳細については、例えば上記非特許文献２に記載されている。
【００３８】
出力部２６は、公知のプリンタまたはディスプレイであって、変形演算部１６の計算結果を出力する。例えば、各時間ステップｎにおける、第１モデルの圧縮変形状態および応力分布を、画像として表示出力する。また、各粒子の位置や速度応力などの物理量の数値情報を表示出力してもよい。また、モデル作成部１２におけるモデル作成処理や、移動・圧縮演算部１４における移動・圧縮シミュレーションなど、各部の処理の経過についても、順次出力することができる。
【００３９】
図４は、本発明の粒子モデルを用いた変形挙動シミュレーション方法の一例のフローチャート図である。本発明の粒子モデルを用いた変形挙動シミュレーション方法では、まず、シミュレーション条件を設定する（ステップＳ１００）。例えば、第１構造体および第２構造体の大きさや形状、第１構造体に対する第２構造体の移動速度などの圧縮条件（どのくらいの速度で、どのくらいの時間ステップをかけて、どの程度圧縮するかなど）や、粒子の条件など、シミュレーションに必要な種々の条件を設定しておく。次に、モデル作成部１２において、上述のように、入力操作系２４を用いて入力されたモデル作成条件に応じて、第１構造体を複数の粒子の集合体として再現した第１モデル、および、第２構造体を再現した第２モデルを作成する（ステップＳ１０２）。
【００４０】
次に、移動・圧縮演算部１４が、入力操作系２４を用いて入力された圧縮条件に基づいて、第１モデル３２に対して第２モデルを移動させるシミュレーション演算を実施する（ステップＳ１０４）。本実施形態では、第１モデル３２に対して第２モデル３４を一定速度（例えば１０ｍｍ／ｓｅｃで近づけ）、接触した後も同一速度で第２モデルを移動させて、第１モデルを圧縮していくシミュレーション演算を行う（ステップＳ１０４）。このステップでは、所定の単位時間間隔（時間ステップ）分だけ、第１モデル３２に第２モデル３４を近づけていく。そして、移動・圧縮演算部１４が、この時間間隔分の移動・圧縮処理が終了した時点での、第１モデル３２の全粒子のうち、第２モデルと接触したとみなす粒子、すなわち接触条件範囲内に第２モデルが入った粒子を抽出する。そして、上述のように、第２モデルと接触したとみなす粒子については、接触した段階における第２モデルの速度ベクトルを、この粒子の速度ベクトルとして付与する（ステップＳ１０６）。この速度ベクトル（速度情報）の付与は、接触以前の粒子の位置や速度などの物理量には全く依存させることなく、無条件で強制的に行う。
【００４１】
そして、変形演算部１６は、上述のように、移動・圧縮演算部１４において第１モデル３２の接触部分の粒子に付与された速度ベクトルを外力条件として、現在の時間ステップが終了した時点での、第１モデルの圧縮変形状態（各粒子の位置や応力など）を計算する（ステップＳ１０８）。変形演算部１６における圧縮変形状態の演算が終了すると、予め定められた圧縮条件で規定された全ての時間ステップについて、圧縮変形状態が算出されたか否かが判定される（ステップＳ１１０）。このような判定は、例えばＣＰＵ２０が行なえばよい。すべての時間ステップについての圧縮変形状態の演算が終了していない場合は、次の時間ステップに変更されて（ステップ＜ｎ＞であった場合、ステップ＜ｎ＋１＞に変更）（ステップＳ１１２）、ステップＳ１０４〜Ｓ１１０の処理が繰り返し実施される。このような処理の繰り返しは、全ての時間ステップについて圧縮状態が演算されるまで、繰り返し実施される。
【００４２】
全ての時間ステップについて圧縮状態が演算されると、出力部２６は、シミュレーション結果を出力する（ステップＳ１１４）。例えば、各時間ステップ終了時点での第１モデルの変形状態を、各時間ステップ毎に連続表示することで、第１モデルと第２モデルが接触してから第１モデルが圧縮されてしまうまでの、第１モデルの変形挙動をアニメーション表示する。本発明の粒子モデルを用いた変形挙動シミュレーション方法は、このように実施される。
【００４３】
図５（ａ）〜（ｄ）は、本願発明の粒子モデルを用いた変形挙動シミュレーション方法の実施結果の一例であり、ゴム体である円柱形状の第１構造体を、略剛体の第２構造体によって圧縮した際の、第１構造体の変形挙動を再現シミュレーションを実施した結果である。図５（ａ）〜（ｄ）に示す例では、第１構造体をφ１００ｍｍのゴム体（超弾性材料）として、この第１構造体を再現する、複数の粒子からなる第１モデル４２を作成し、平板状の略剛体を再現する第２モデル４４を作成した。そして、第１モデル４２に対し、図中上側から下方向に向かって、第２モデル４４（の一方）を速度１０ｍｍ／ｓｅｃで移動させて、第１モデル４２の圧縮変形状態を演算した。なお、圧縮変形状態の演算では、各圧縮状態における（すなわち各時間ステップにおける）各粒子の位置および応力を求めた。なお、応力を求めるにあたり、第１モデル４２の弾性ポテンシャル関数としてＭｏｏｎｅｙ−Ｒｉｖｌｉｎ体を用いた。なお、Ｍｏｏｎｅｙ−Ｒｉｖｌｉｎ定数ｃ１＝１．０、ｃ２＝０．０４、微圧縮材料の体積変化の係数α＝７．５とした。
【００４４】
図５（ａ）は、第１モデル４２と第２モデル４４とが一部分が接触して、双方のモデルが停止している初期状態であり、図５（ｂ）は、第１モデル４２が図中上下方向に１０％圧縮された状態の結果であり、図５（ｃ）は、第１モデル４２が図中上下方向に２５％圧縮された状態における圧縮変形状態の結果であり、図５（ｄ）は、第１モデル４２が図中上下方向に４０％圧縮された状態における圧縮変形状態の結果である。図５（ａ）〜（ｄ）では、色が薄いほど圧縮される方向の応力が強く、色が濃いほど引っ張られる方向の応力が強い。図５（ａ）〜（ｄ）に示すように、本願発明の粒子モデルを用いた変形挙動シミュレーション方法によって、超弾性体である第１構造体の変形挙動を各圧縮状態毎に詳細に求めることができる。なお、図６（ａ）は、図５（ｄ）に示す状態の第１モデルについて、一部を拡大して示す図である。また、図６（ｂ）は、図５（ａ）〜（ｄ）に示す場合と同様の第１構造体および第２構造体をメッシュモデルによって再現したモデルを用い、ＦＥＭ（有限要素法）によって、図５（ａ）〜（ｄ）に示す場合と同様の圧縮変形解析を行なった結果である。図６（ｂ）は、図５（ｄ）と同様、第１モデルが図中上下方向に４０％圧縮された状態における、圧縮変形状態の結果である。図６（ａ）（図５（ｄ））に示す本発明のシミュレーション結果は、図６（ｂ）に示す有限要素法による解析結果とよく一致している。本発明の粒子モデルを用いた変形挙動シミュレーション方法は、ゴム体などの大変形問題を、十分な精度で解析可能である。
【００４５】
以上、本発明の粒子モデルを用いた変形挙動シミュレーション方法およびプログラムについて詳細に説明したが、本発明は上記実施形態に限定されず、本発明の主旨を逸脱しない範囲において、種々の改良や変更をしてもよいのはもちろんである。
【図面の簡単な説明】
【００４６】
【図１】本発明の粒子モデルを用いた変形挙動シミュレーション方法を実施する変形挙動シミュレーション装置の概略の構成を示すブロック図である。
【図２】本発明の変形挙動シミュレーション方法によって作成されるモデルの一例である。
【図３】（ａ)は、従来の変形挙動シミュレーション方法で作成されるモデルの一例であり、（ｂ）は、本発明の変形挙動シミュレーション方法で作成されるモデルの他の例である。
【図４】本発明の粒子モデルを用いた変形挙動シミュレーション方法の一例のフローチャート図である。
【図５】（ａ）〜（ｄ）は、本発明の粒子モデルを用いた変形挙動シミュレーション方法の実施結果の一例であり、第１構造体を第２構造体によって圧縮した際の第１構造体の変形挙動を再現するシミュレーションを実施した結果である。
【図６】（ａ）は、図５（ｄ）に示す状態の第１モデルについて、一部を拡大して示す図である。また（ｂ）は、図５（ａ）〜（ｄ）に示す場合と同様の第１構造体および第２構造体をメッシュモデルによって再現したモデルを用い、ＦＥＭ（有限要素法）によって、図５（ａ）〜（ｄ）に示す場合と同様の変形挙動を再現するシミュレーションを実施した結果である。
【符号の説明】
【００４７】
１０  シミュレーション装置
１２  モデル作成部
１４  移動・圧縮演算部
１６  変形演算部
１８  メモリ
２０  ＣＰＵ
２２  本体装置
２４  入力操作系
２６  出力部
３２、４２  第１モデル
３４、４４  第２モデル
【図１】
2008152423.tif 000020
【図２】
2008152423.tif 000021
【図４】
2008152423.tif 000022
【図５】
2008152423.tif 000023
【図３】
2008152423.tif 000024
【図６】
2008152423.tif 000025