(57)【要約】【課題】連続体の運動を解く数値計算シミュレーションを精度良く行うこと。【解決手段】シミュレーションプログラム３３０ａは、コンピュータ３００に、変形可能な連続体を複数の粒子で表した場合における複数の粒子のそれぞれの球形状のポテンシャル領域を、複数の粒子のそれぞれの変形量に応じて、楕円球形状に変形させる。シミュレーションプログラム３３０ａは、コンピュータ３００に、楕円球形状に変形された、複数の粒子のそれぞれのポテンシャル領域に基づいて、粒子間の反力を計算させる。【選択図】図７
2014146299.tif 000002
【特許請求の範囲】
【請求項１】
コンピュータに、
変形可能な連続体を複数の粒子で表した場合における前記複数の粒子のそれぞれの球形状のポテンシャル領域を、前記複数の粒子のそれぞれの変形量に応じて、楕円球形状に変形させ、
楕円球形状に変形された、前記複数の粒子のそれぞれのポテンシャル領域に基づいて、粒子間の反力を計算させることを特徴とするシミュレーションプログラム。
【請求項２】
前記複数の粒子のそれぞれの球形状のポテンシャル領域の楕円球形状への変形は、前記複数の粒子のそれぞれの球形状のポテンシャル領域を、前記複数の粒子のそれぞれの変形勾配テンソルを用いて、楕円球形状に変形させることを特徴とする請求項１記載のシミュレーションプログラム。
【請求項３】
コンピュータが、
変形可能な連続体を複数の粒子で表した場合における前記複数の粒子のそれぞれの球形状のポテンシャル領域を、前記複数の粒子のそれぞれの変形量に応じて、楕円球形状に変形し、
楕円球形状に変形された、前記複数の粒子のそれぞれのポテンシャル領域に基づいて、粒子間の反力を計算することを特徴とするシミュレーション方法。
【請求項４】
変形可能な連続体を複数の粒子で表した場合における前記複数の粒子のそれぞれの球形状のポテンシャル領域を、前記複数の粒子のそれぞれの変形量に応じて、楕円球形状に変形させる変形部と、
楕円球形状に変形された、前記複数の粒子のそれぞれのポテンシャル領域に基づいて、粒子間の反力を計算する計算部と
を有することを特徴とするシミュレーション装置。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、シミュレーションプログラム、シミュレーション方法及びシミュレーション装置に関する。
【背景技術】
【０００２】
従来、流体や弾性体等の連続体の運動を解く数値計算の手法として、格子をベースにして微分方程式の近似解を求解する有限差分法や有限要素法、有限体積法などが用いられている。また、近年では、数値計算をＣＡＥ（Computer  Aided  Engineering）等の応用分野で活用するため、連続体の状態を解く数値計算の手法も発展し、流体と構造物とが相互作用する問題が解かれる。しかしながら、格子を用いる数値計算の手法では、自由表面等の界面の存在する問題や、流体と構造とを連成して解く流体構造連成解析問題などの移動境界が発生する場合には、連続体の取り扱いが複雑になるため、プログラム作成が困難になる場合がある。
【０００３】
ここで、格子を用いない数値計算の手法として、粒子法がある。粒子法とは、連続体の運動を有限の数の粒子の運動として解析する手法である。現在提案されている代表的な粒子法としては、ＳＰＨ（Smoothed  Particles  Hydrodynamics）法やＭＰＳ（Moving  Particles  Semi-implicit）法といったものがある。粒子法では、移動境界の取り扱いに特別な処置をせずに連続体の運動を解析することができる。そのため、近年、連続体の運動を解く数値計算の手法として広く用いられるようになってきている。
【０００４】
構造解析の分野では、物体同士の衝突などの接触問題の計算を行う場合がある。有限要素法をベースとしたＬＳ−ＤＹＮＡ（登録商標）等のソフトウェアは、計算格子の幾何形状から接触問題を扱う。しかしながら、粒子法は、計算格子を生成しない。このため、粒子法は、接触問題を計算する場合に、例えば、粒子に予め指定した半径ｈｓ内に別の粒子が入ってきた際に、互いの粒子に反力が働くように計算する方法を用いる。図８は、従来の粒子法により計算される反力の一例を示す図である。図８の例は、粒子９０の半径ｈｓ内に粒子９１が存在する場合を示す。また、図８の例は、粒子９１の半径ｈｓ内に粒子９０が存在する場合を示す。図８の例が示す場合において、粒子法は、粒子９０からの粒子９１に働く反力９０ａを計算する。また、粒子法は、粒子９１からの粒子９０に働く反力９１ａを計算する。
【０００５】
また、粒子法は、例えば、粒子ｉ，ｊ間に働くポテンシャル関数として以下の式（１）を用いる。
【０００６】
【数１】
2014146299.tif 000003
【０００７】
ここで、ｘｉ，ｘｊは、それぞれ、粒子ｉ，粒子ｊの位置ベクトルである。また、ｃは、定数である。
【０００８】
式（１）が示すポテンシャルにより粒子ｉに働く粒子ｊからの反力が、
【数２】
2014146299.tif 000004
として得られ、粒子ｊに働く粒子ｉからの反力が、
【数３】
2014146299.tif 000005
として得られる。
【０００９】
また、粒子法で接触問題を計算する方法の他の例としては、粒子の接触点を求める際に球形以外の閉曲面の接触点を求める方法もある。
【先行技術文献】
【特許文献】
【００１０】
【特許文献１】特開２００９−２６２７９号公報
【発明の概要】
【発明が解決しようとする課題】
【００１１】
しかしながら、従来の粒子法は、ポテンシャル領域を球形で表現する。このため、従来の粒子法では、弾性体などの変形体同士の接触問題を扱う際に、変形前は、粒子が侵入できなかった領域にも、変形後は、粒子が侵入できるようになってしまう場合がある。かかる場合、ある粒子の計算領域上に他の粒子が存在するため、従来の粒子法では、反力などの計算の精度が悪くなるという問題がある。図９及び図１０は、従来の粒子法の問題点の一例を説明するための図である。図９は、変形前、例えば、従来の粒子法を用いたシミュレーションにおける初期の時点での弾性体の粒子８０の配置の一例を示す。図９の例に示すように、２つの粒子８０は、互いに、半径ｈｓ内の領域が重なるように、配置されている。このため、粒子８０を含む弾性体とは異なる弾性体の粒子８１は、２つの粒子８０の間に侵入できない。
【００１２】
図１０は、変形後、例えば、従来の粒子法を用いたシミュレーションにおいてある程度シミュレーションが進んだ時点での弾性体の粒子８０の配置の一例を示す。図１０の例に示すように、２つの粒子８０は、互いに、離れる方向に変形され、互いの半径ｈｓ内の領域が重ならないように変形されている。このため、粒子８０を含む弾性体とは異なる弾性体の粒子８１は、２つの粒子８０の間に侵入することが可能となる。粒子８１が、２つの粒子８０の間に侵入するような場合には、粒子８１の計算領域上に粒子８０が存在するとともに、粒子８０の計算領域上に粒子８１が存在する。このため、図１０の例が示す場合では、反力などの計算の精度が悪くなる。
【００１３】
１つの側面では、連続体の運動を解く数値計算シミュレーションを精度良く行うことを目的とする。
【課題を解決するための手段】
【００１４】
本願の開示するシミュレーションプログラムは、１つの態様において、コンピュータに次の処理を実行させる。すなわち、シミュレーションプログラムは、コンピュータに、変形可能な連続体を複数の粒子で表した場合における複数の粒子のそれぞれの球形状のポテンシャル領域を、複数の粒子のそれぞれの変形量に応じて、楕円球形状に変形させる。シミュレーションプログラムは、コンピュータに、楕円球形状に変形された、複数の粒子のそれぞれのポテンシャル領域に基づいて、粒子間の反力を計算させる。
【発明の効果】
【００１５】
連続体の運動を解く数値計算シミュレーションを精度良く行うことができる。
【図面の簡単な説明】
【００１６】
【図１】図１は、実施例に係るシミュレーション装置の機能構成の一例を示す図である。
【図２】図２は、実施例に係るシミュレーション装置が実行する処理の一例を説明するための図である。
【図３】図３は、実施例に係るシミュレーション装置が実行する処理の一例を説明するための図である。
【図４】図４は、従来の粒子法を用いた場合の処理の一例を説明するための図である。
【図５】図５は、実施例に係るシミュレーション処理の手順を示すフローチャートである。
【図６】図６は、実施例に係る反力計算処理の手順を示すフローチャートである。
【図７】図７は、シミュレーションプログラムを実行するコンピュータを示す図である。
【図８】図８は、従来の粒子法により計算される反力の一例を示す図である。
【図９】図９は、従来の粒子法の問題点の一例を説明するための図である。
【図１０】図１０は、従来の粒子法の問題点の一例を説明するための図である。
【発明を実施するための形態】
【００１７】
以下に、本願の開示するシミュレーションプログラム、シミュレーション方法及びシミュレーション装置の実施例を図面に基づいて詳細に説明する。なお、実施例は開示の技術を限定するものではない。
【００１８】
［シミュレーション装置の構成］
実施例に係るシミュレーション装置について説明する。本実施例に係るシミュレーション装置は、弾性体、弾塑性体や粘弾性体などの変形可能な物体が、他の変形可能な物体との接触問題を計算する。本実施例に係るシミュレーション装置は、各粒子の位置、速度、加速度、変形勾配テンソル、温度、密度をタイプステップｔｔｓ毎に計算する。図１は、実施例に係るシミュレーション装置の機能構成の一例を示す図である。図１に示すように、シミュレーション装置１０は、入力部１１と、表示部１２と、記憶部１３と、制御部１４とを有する。
【００１９】
入力部１１は、制御部１４に情報を入力する。例えば、入力部１１は、ユーザから後述のシミュレーション処理を実行する指示であるシミュレーション実行指示を受け付けて、受け付けたシミュレーション実行指示を制御部１４に入力する。また、入力部１１は、ユーザから、初期状態における各粒子の初期値を受け付けて、受け付けた各粒子の初期値を制御部１４に入力する。ここで、初期状態における各粒子の初期値には、各粒子の位置、速度、加速度、変形勾配テンソル、温度、密度が含まれる。入力部１１のデバイスの一例としては、キーボードやマウスなどが挙げられる。
【００２０】
表示部１２は、各種の情報を表示する。例えば、表示部１２は、後述の表示制御部１４ｃの制御によりシミュレーション結果を表示する。表示部１２のデバイスの一例としては、液晶ディスプレイなどが挙げられる。
【００２１】
記憶部１３は、制御部１４で実行される各種プログラムを記憶する。また、記憶部１３は、変形モデルデータ１３ａを記憶する。変形モデルデータ１３ａは、弾性体、弾塑性体や粘弾性体などの変形可能な連続体を複数の粒子として表したモデルを示す。
【００２２】
図１に戻り、記憶部１３は、例えば、フラッシュメモリなどの半導体メモリ素子、または、ハードディスク、光ディスクなどの記憶装置である。なお、記憶部１３は、上記の種類の記憶装置に限定されるものではなく、ＲＡＭ（Random  Access  Memory）、ＲＯＭ（Read  Only  Memory）であってもよい。
【００２３】
制御部１４は、各種の処理手順を規定したプログラムや制御データを格納するための内部メモリを有し、これらによって種々の処理を実行する。図１に示すように、制御部１４は、計算部１４ａと、変形部１４ｂと、表示制御部１４ｃとを有する。
【００２４】
計算部１４ａは、各種の情報を計算する。例えば、計算部１４ａは、各粒子の位置、速度、加速度、変形勾配テンソル、温度、密度をタイプステップｔｔｓ毎に計算する。また、計算部１４ａは、後述の変形部１４ｂにより楕円球形状に変形された、複数の粒子のそれぞれのポテンシャル領域に基づいて、粒子間の反力を計算する。
【００２５】
計算部１４ａの一態様について説明する。例えば、計算部１４ａは、入力部１１からシミュレーション実行指示が入力された場合には、まず、タイムステップｔｔｓの値を０に設定する。そして、計算部１４ａは、入力部１１から各粒子の初期値が入力されたか否かを判定する。初期値が入力された場合には、計算部１４ａは、タイムステップｔｔｓの値を１だけインクリメントする。また、計算部１４ａは、タイムステップｔｔｓの値がシミュレーションの最後のタイムステップＮＬ以下であると表示制御部１４ｃにより判定された場合にも、タイムステップｔｔｓの値を１だけインクリメントする。
【００２６】
続いて、計算部１４ａは、タイムステップｔｔｓにおける、各粒子にかかる力、例えば、応力や外力などを計算する。
【００２７】
続いて、計算部１４ａは、全ての粒子のうち、粒子ｊとして未選択の粒子があるか否かを判定する。粒子ｊとして未選択の粒子がある場合には、計算部１４ａは、未選択の粒子を粒子ｊとして１つ選択する。そして、計算部１４ａは、全ての粒子から、新たに粒子ｊとして選択された粒子を除いた粒子のうち、粒子ｉとして未選択の粒子があるか否かを判定する。粒子ｉとして未選択の粒子がある場合には、計算部１４ａは、未選択の粒子を粒子ｉとして１つ選択する。ここで、変形後の粒子ｉの位置ベクトルをｘｉ、変形前の粒子ｉの位置ベクトルをｘｉ＊とする。本実施例では、各粒子は、変形量を表す変形勾配テンソルを持つ。例えば、粒子ｉが持つ変形勾配テンソルをＦｉとする。そして、本実施例では、各粒子にかかる力が計算された後に、各粒子にかかる反力が計算部１４ａ及び変形部１４ｂによって計算される。
【００２８】
変形部１４ｂについて説明する。例えば、変形部１４ｂは、計算部１４ａにより未選択の粒子が粒子ｉとして選択された場合に、下記の式（２）に従って、粒子ｊの持つ変形勾配テンソルＦｊを用いて、変形前の座標系での相対位置ベクトルｘｉｊ，ｊ＊を計算する。
【００２９】
【数４】
2014146299.tif 000006
【００３０】
ただし、粒子ｉと粒子ｊとの相対位置ベクトルをｘｉｊ＝ｘｉ−ｘｊとする。粒子ｊが変形勾配ベクトルＦｊに応じて変形を受けているので、式（２）が示す逆変換により、粒子ｊが変化していない状態での粒子ｉと粒子ｊとの相対距離を示す相対位置ベクトルが算出される。
【００３１】
そして、変形部１４ｂは、下記の式（３）を用いて、粒子ｊからの反力ポテンシャルを定義する。
【００３２】
【数５】
2014146299.tif 000007
【００３３】
ここで、ｃｉｊ，ｊは、定数である。
【００３４】
上述したように、変形部１４ｂは、式（２）を用いて、変形前の座標系での相対位置ベクトルｘｉｊ，ｊ＊を計算し、式（３）を用いて、粒子ｊからの反力ポテンシャルを定義する。これにより、変形部１４ｂは、粒子ｊの球形状のポテンシャル領域を粒子ｊの変形量に応じて楕円球形状に変形させる。例えば、変形部１４ｂは、粒子ｊの持つ変形勾配テンソルＦｊを用いて、粒子ｊの球形状のポテンシャル領域を楕円球形状に変形させる。図２は、実施例に係るシミュレーション装置が実行する処理の一例について説明するための図である。図２の例では、粒子７０が粒子ｊに対応する。図２の例に示すように、本実施例によれば、粒子７０のポテンシャル領域の形状が楕円球形状に変形される。また、図２の例に示すように、粒子７０のポテンシャル領域の形状が楕円球形状に変形されると、粒子７１が侵入できないようになる。そのため、本実施例によれば、粒子７０の計算領域上に粒子７１が存在してしまうような状況が発生することが抑制される。そのため、本実施例によれば、精度良く反力などの計算を行うことができる。
【００３５】
そして、計算部１４ａは、式（３）のポテンシャルの勾配計算から、下記の式（４）にしたがって、粒子ｉにかかる力ｆｉｊ，ｊを計算する。
【００３６】
【数６】
2014146299.tif 000008
【００３７】
ここで、ｘｉ，ｊ＊＝Ｆｊ−１ｘｉである。さらに、粒子ｊには、式（４）が示す力ｆｉｊ，ｊと同じ大きさで逆向きの力ｆｊｉ，ｊ＝−ｆｉｊ，ｊを反作用として受けるものとする。計算部１４ａは、下記の式（５）にしたがって、粒子ｉが粒子ｊに与える反力の反作用を加味して、粒子ｉが粒子ｊから受ける反力ｆｃ，ｉｊを計算する。
【００３８】
【数７】
2014146299.tif 000009
【００３９】
そして、計算部１４ａは、全ての粒子から、新たに粒子ｊとして選択された粒子を除いた粒子のうち、粒子ｉとして未選択の粒子があるか否かを判定する上述した処理以降の処理を再び行う。これにより、計算部１４ａは、粒子ｊとして選択された粒子から、粒子ｊ以外の全ての粒子のそれぞれが受ける反力を計算することができる。
【００４０】
そして、全ての粒子から、新たに粒子ｊとして選択された粒子を除いた粒子のうち、粒子ｉとして未選択の粒子があるか否かを判定する処理において、粒子ｉとして未選択の粒子がないと判定された場合には、計算部１４ａは、次の処理を行う。すなわち、計算部１４ａは、全ての粒子を粒子ｉとして未選択であると定める。そして、計算部１４ａは、全ての粒子のうち、粒子ｊとして未選択の粒子があるか否かを判定する上述した処理以降の処理を再び行う。これにより、全ての粒子間の反力が計算される。
【００４１】
そして、全ての粒子のうち、粒子ｊとして未選択の粒子がないと判定した場合には、計算部１４ａは、タイムステップｔｔｓにおける、各粒子に加わる力から各粒子の加速度を計算する。
【００４２】
そして、計算部１４ａは、各粒子の速度を、各粒子の加速度から計算した速度を用いて更新する。続いて、計算部１４ａは、各粒子の位置を、各粒子の速度から計算した位置を用いて更新する。その後、計算部１４ａは、各粒子の変形勾配テンソルを、速度から求まる各粒子の変形量を用いて更新する。続いて、計算部１４ａは、全ての粒子の更新結果と、タイムステップｔｔｓとを対応付けて、記憶部１３に格納する。
【００４３】
表示制御部１４ｃは、各種の情報の表示を制御する。表示制御部１４ｃの一態様について説明する。例えば、表示制御部１４ｃは、全ての粒子の更新結果とタイムステップｔｔｓとが対応付けられて計算部１４ａにより記憶部１３に格納された場合に、タイムステップｔｔｓの値が、シミュレーションの最後のタイムステップＮＬ以下であるか否かを判定する。タイムステップｔｔｓの値が、シミュレーションの最後のタイムステップＮＬ以下でない場合には、表示制御部１４ｃは、次の処理を行う。すなわち、表示制御部１４ｃは、記憶部１３に格納されたタイムステップごとの全ての粒子の更新結果を全てのタイムステップについて取得する。そして、表示制御部１４ｃは、シミュレーション結果（全てのタイムステップにおける全ての粒子の更新結果）を表示するように表示部１２の表示を制御する。
【００４４】
制御部１４は、ＡＳＩＣ（Application  Specific  Integrated  Circuit）やＦＰＧＡ（Field  Programmable  Gate  Array）などのハードワイヤードロジックである。または、制御部１４は、ＣＰＵ（Central  Processing  Unit）やＭＰＵ（Micro  Processing  Unit）などにプログラムを実行させることにより実現される。
【００４５】
ここで、従来の粒子法を用いた場合の処理の一例について説明する。図４は、従来の粒子法を用いた場合の処理の一例を説明するための図である。図４の例では、弾性体５０を離散化（粒子化）して、複数の粒子５０ａを含むモデルを生成し、弾性体５１を離散化して、複数の粒子５１ａを含むモデルを生成した場合を示す。従来の粒子法では、複数の粒子５０ａを含む弾性体５０のモデルと、複数の粒子５１ａを含む弾性体５１のモデルとを衝突させて変形計算を行うと、図４に示すように、例えば、２つの粒子５０ａが互いに遠ざかるように変形する場合がある。この場合、図４に示すように、２つの粒子５０ａの間に、粒子５１ａが侵入することができる。すなわち、粒子５１ａは、粒子５０ａからの反力を受けない。このとき、粒子５０ａは、粒子５１ａの計算領域上に存在し、粒子５１ａは、粒子５０ａの計算領域上に存在することになる。それゆえ、従来の粒子法では、反力などの計算の精度が悪くなるという問題がある。
【００４６】
一方、本実施例によれば、先の図２の例に示すように、変形して粒子間の距離が離れても、粒子のそれぞれの球形状のポテンシャル領域が楕円球形状に変形されるので、計算領域に他の粒子が存在するような事態が発生することが抑制される。それゆえ、本実施例では、精度良く、反力などを計算することができる。
【００４７】
［処理の流れ］
次に、本実施例に係るシミュレーション装置１０の処理の流れを説明する。図５は、実施例に係るシミュレーション処理の手順を示すフローチャートである。このシミュレーション処理の実行タイミングとしては様々なタイミングが考えられる。例えば、シミュレーション処理は、シミュレーション処理を実行するシミュレーション実行指示が入力部１１から入力された場合に、制御部１４により実行される。
【００４８】
図５に示すように、計算部１４ａは、タイムステップｔｔｓの値を０に設定する（Ｓ１０１）。そして、計算部１４ａは、入力部１１から各粒子の初期値が入力されたか否かを判定する（Ｓ１０２）。初期値が入力されていない場合（Ｓ１０２；Ｎｏ）には、計算部１４ａは、再び、Ｓ１０２での判定を行う。一方、初期値が入力された場合（Ｓ１０２；Ｙｅｓ）には、計算部１４ａは、タイムステップｔｔｓの値を１だけインクリメントする（Ｓ１０３）。
【００４９】
続いて、計算部１４ａは、タイムステップｔｔｓにおける、各粒子にかかる力、例えば、応力や外力などを計算する（Ｓ１０４）。続いて、計算部１４ａは、反力計算処理を実行する（Ｓ１０５）。そして、計算部１４ａは、タイムステップｔｔｓにおける、各粒子に加わる力から各粒子の加速度を計算する（Ｓ１０６）。
【００５０】
そして、計算部１４ａは、各粒子の速度を、各粒子の加速度から計算した速度を用いて更新する（Ｓ１０７）。続いて、計算部１４ａは、各粒子の位置を、各粒子の速度から計算した位置を用いて更新する（Ｓ１０８）。その後、計算部１４ａは、各粒子の変形勾配テンソルを、速度から求まる各粒子の変形量を用いて更新する（Ｓ１０９）。続いて、計算部１４ａは、全ての粒子の更新結果と、タイムステップｔｔｓとを対応付けて、記憶部１３に格納する（Ｓ１１０）。
【００５１】
そして、表示制御部１４ｃは、タイムステップｔｔｓの値が、シミュレーションの最後のタイムステップＮＬ以下であるか否かを判定する（Ｓ１１１）。タイムステップｔｔｓの値が、シミュレーションの最後のタイムステップＮＬ以下である場合（Ｓ１１１；Ｙｅｓ）には、Ｓ１０３に戻る。一方、タイムステップｔｔｓの値が、シミュレーションの最後のタイムステップＮＬ以下でない場合（Ｓ１１１；Ｎｏ）には、表示制御部１４ｃは、次の処理を行う。すなわち、表示制御部１４ｃは、記憶部１３に格納されたタイムステップごとの全ての粒子の更新結果を全てのタイムステップについて取得する。そして、表示制御部１４ｃは、シミュレーション結果（全てのタイムステップにおける全ての粒子の更新結果）を表示するように表示部１２の表示を制御し（Ｓ１１２）、処理を終了する。
【００５２】
図６は、実施例に係る反力計算処理の手順を示すフローチャートである。図６に示すように、計算部１４ａは、全ての粒子のうち、粒子ｊとして未選択の粒子があるか否かを判定する（Ｓ２０１）。粒子ｊとして未選択の粒子がある場合（Ｓ２０１；Ｙｅｓ）には、計算部１４ａは、未選択の粒子を粒子ｊとして１つ選択する（Ｓ２０２）。そして、計算部１４ａは、全ての粒子から、新たに粒子ｊとして選択された粒子を除いた粒子のうち、粒子ｉとして未選択の粒子があるか否かを判定する（Ｓ２０３）。粒子ｉとして未選択の粒子がある場合（Ｓ２０３；Ｙｅｓ）には、計算部１４ａは、未選択の粒子を粒子ｉとして１つ選択する（Ｓ２０４）。
【００５３】
変形部１４ｂは、式（２）に従って、粒子ｊの持つ変形勾配テンソルＦｊを用いて、変形前の座標系での相対位置ベクトルｘｉｊ，ｊ＊を計算する（Ｓ２０５）。そして、変形部１４ｂは、式（３）を用いて、粒子ｊからの反力ポテンシャルを定義する（Ｓ２０６）。
【００５４】
そして、計算部１４ａは、式（３）のポテンシャルの勾配計算から、式（４）にしたがって、粒子ｉにかかる力ｆｉｊ，ｊを計算する（Ｓ２０７）。そして、計算部１４ａは、式（５）にしたがって、粒子ｉが粒子ｊに与える反力の反作用を加味して、粒子ｉが粒子ｊから受ける反力ｆｃ，ｉｊを計算し（Ｓ２０８）、Ｓ２０３に戻る。
【００５５】
一方、粒子ｉとして未選択の粒子がない場合（Ｓ２０３；Ｎｏ）には、計算部１４ａは、全ての粒子を粒子ｉとして未選択であると定め（Ｓ２０９）、Ｓ２０１に戻る。
【００５６】
また、粒子ｊとして未選択の粒子がない場合（Ｓ２０１；Ｎｏ）には、計算部１４ａは、処理結果を内部メモリに格納し、リターンする。
【００５７】
上述してきたように、本実施例に係るシミュレーション装置１０は、変形可能な連続体を複数の粒子として表した場合の複数の粒子のそれぞれの球形状のポテンシャル領域を、複数の粒子のそれぞれの変形量に応じて、楕円球形状に変形させる。そして、シミュレーション装置１０は、楕円球形状に変形された、複数の粒子のそれぞれのポテンシャル領域に基づいて、粒子間の反力を計算する。それゆえ、本実施例に係るシミュレーション装置１０によれば、粒子間の距離が離れるような変形が発生しても、互いの粒子のポテンシャル領域が重なっているような状況が発生する。そのため、本実施例に係るシミュレーション装置１０によれば、粒子間に他の粒子が侵入するような事象の発生が抑制される。したがって、本実施例に係るシミュレーション装置１０によれば、計算領域に他の粒子が存在するような事態が発生することが抑制される。それゆえ、本実施例に係るシミュレーション装置１０によれば、連続体の運動を解く数値計算シミュレーションを精度良く行うことができる。
【００５８】
また、本実施例に係るシミュレーション装置１０は、複数の粒子のそれぞれの球形状のポテンシャル領域を、複数の粒子のそれぞれの変形勾配テンソルを用いて、楕円球形状に変形させることもできる。
【００５９】
また、シミュレーション装置１０によれば、式（２）に示すように逆変換を用いて、楕円球形状のポテンシャル領域を球形状に戻すことにより、接触点を２点間の相対位置ベクトルから求めるという簡易な方法で接触点を計算することができる。したがって、シミュレーション装置１０によれば、接触点を計算する際の計算時間を短縮することができる。
【００６０】
さて、これまで開示の装置に関する実施例について説明したが、本発明は上述した実施例以外にも、種々の異なる形態にて実施されてよい。例えば、実施例などにおいて説明した各処理のうち、自動的に行われるものとして説明した処理の全部または一部を手動的に行うこともできる。
【００６１】
また、各種の負荷や使用状況などに応じて、実施例などにおいて説明した各処理の各ステップでの処理を任意に細かくわけたり、あるいはまとめたりすることができる。また、ステップを省略することもできる。
【００６２】
また、各種の負荷や使用状況などに応じて、実施例などにおいて説明した各処理の各ステップでの処理の順番を変更できる。
【００６３】
また、図示した装置の各構成要素は機能概念的なものであり、必ずしも物理的に図示の如く構成されていることを要しない。すなわち、装置の分散・統合の具体的状態は図示のものに限られず、その全部または一部を、各種の負荷や使用状況などに応じて、任意の単位で機能的または物理的に分散・統合して構成することができる。
【００６４】
［シミュレーションプログラム］
また、上記のシミュレーション装置１０のシミュレーション処理は、あらかじめ用意されたプログラムをパーソナルコンピュータやワークステーションなどのコンピュータシステムで実行することによって実現することもできる。そこで、以下では、図７を用いて、上記のシミュレーション装置１０と同様の機能を有するシミュレーションプログラムを実行するコンピュータの一例を説明する。
【００６５】
図７は、シミュレーションプログラムを実行するコンピュータを示す図である。図７に示すように、コンピュータ３００は、ＣＰＵ（Central  Processing  Unit）３１０、ＲＯＭ（Read  Only  Memory）３２０、ＨＤＤ（Hard  Disk  Drive）３３０、ＲＡＭ（Random  Access  Memory）３４０を有する。これら３００〜３４０の各部は、バス３５０を介して接続される。
【００６６】
ＨＤＤ３３０には、上記の実施例で示す計算部１４ａ、変形部１４ｂ、表示制御部１４ｃと同様の機能を発揮するシミュレーションプログラム３３０ａが予め記憶される。なお、シミュレーションプログラム３３０ａについては、適宜分離しても良い。
【００６７】
そして、ＣＰＵ３１０が、シミュレーションプログラム３３０ａをＨＤＤ３３０から読み出して実行する。
【００６８】
そして、ＨＤＤ３３０には、図１の例に示す記憶部１３に記憶された変形モデルデータが設けられる。
【００６９】
そして、ＣＰＵ３１０は、ＨＤＤ３３０からデータを読み出してＲＡＭ３４０に格納する。さらに、ＣＰＵ３１０は、ＲＡＭ３４０に格納された各種のデータを用いて、シミュレーションプログラム３３０ａを実行する。なお、ＲＡＭ３４０に格納される各データは、常に全てのデータがＲＡＭ３４０に格納されなくともよく、全てのデータのうち処理に用いられるデータのみがＲＡＭ３４０に格納されれば良い。
【００７０】
なお、上記したシミュレーションプログラム３３０ａについては、必ずしも最初からＨＤＤ３３０に記憶させなくともよい。
【００７１】
例えば、コンピュータ３００に挿入されるフレキシブルディスク（ＦＤ）、ＣＤ−ＲＯＭ、ＤＶＤディスク、光磁気ディスク、ＩＣカードなどの「可搬用の物理媒体」にプログラムを記憶させておく。そして、コンピュータ３００がこれらからプログラムを読み出して実行するようにしてもよい。  
【００７２】
さらには、公衆回線、インターネット、ＬＡＮ、ＷＡＮなどを介してコンピュータ３００に接続される「他のコンピュータ（またはサーバ）」などにプログラムを記憶させておく。そして、コンピュータ３００がこれらからプログラムを読み出して実行するようにしてもよい。
【符号の説明】
【００７３】
１０      シミュレーション装置
１３      記憶部
１３ａ    変形モデルデータ
１４      制御部
１４ａ    計算部
１４ｂ    変形部
１４ｃ    表示制御部
【図１】
2014146299.tif 000010
【図２】
2014146299.tif 000011
【図３】
2014146299.tif 000012
【図４】
2014146299.tif 000013
【図５】
2014146299.tif 000014
【図６】
2014146299.tif 000015
【図７】
2014146299.tif 000016
【図８】
2014146299.tif 000017
【図９】
2014146299.tif 000018
【図１０】
2014146299.tif 000019