特許 7547288 シミュレーション方法、シミュレーション装置、及びプログラム

(19)【発行国】日本国特許庁(JP)

(12)【公報種別】特許公報(B2)

(11)【特許番号】

(24)【登録日】2024-08-30

(45)【発行日】2024-09-09

(54)【発明の名称】シミュレーション方法、シミュレーション装置、及びプログラム

(51)【国際特許分類】

   G16Z 99/00 20190101AFI20240902BHJP

【ＦＩ】

G16Z99/00

【請求項の数】 5

(21)【出願番号】P 2021107625

(22)【出願日】2021-06-29

(65)【公開番号】P2023005604

(43)【公開日】2023-01-18

【審査請求日】2023-07-14

(73)【特許権者】

【識別番号】000002107

【氏名又は名称】住友重機械工業株式会社

(74)【代理人】

【識別番号】100105887

【弁理士】

【氏名又は名称】来山 幹雄

(72)【発明者】

【氏名】小林 義崇

【審査官】宮地 匡人

(56)【参考文献】

【文献】特開２０１７－０４９９７１（ＪＰ，Ａ）

【文献】SHIIHARA Yoshinori，Fast quasi-implicit NOSB peridynamic simulation based on FIRE algorithm，Mechanical Engineering Journal，2019年04月10日，Vol.6 No.3，Paper No.18-00363

(58)【調査した分野】(Int.Cl.，ＤＢ名)

Ｇ１６Ｚ ９９／００

ＪＳＴＰｌｕｓ（ＪＤｒｅａｍＩＩＩ）

ＪＳＴ７５８０（ＪＤｒｅａｍＩＩＩ）

(57)【特許請求の範囲】

【請求項1】

弾性を有する解析対象物を仮想的な複数の粒子の集合体として表し、
タイムステップごとに各粒子に作用する力を求め、
各粒子に作用する力に基づいて粒子ごとに運動方程式を解いて、各粒子の速度及び位置を更新するシミュレーション方法であって、
粒子ごとに運動方程式を解く際に、
前記解析対象物を剛体とみなし、タイムステップごとに求めた各粒子に作用する力、各粒子の速度及び位置に基づいて、前記解析対象物に作用する力及び前記解析対象物の速度を計算し、
前記解析対象物の運動に対してＦＩＲＥを適用して、前記解析対象物に作用させるべき力を計算し、
前記解析対象物に作用させるべき力を複数の粒子に分配することにより、各粒子に追加的に作用させる力を求め、
各粒子に作用する力に、前記追加的に作用させる力を追加して、運動方程式を解くシミュレーション方法。

【請求項2】

前記解析対象物に作用する力として、並進力及びトルクと計算し、前記解析対象物の速度として、並進速度及び角速度を計算し、
前記解析対象物の運動に対してＦＩＲＥを適用する際に、前記解析対象物の並進運動及び回転運動のそれぞれに対してＦＩＲＥを適用し、
前記解析対象物に作用させるべき力として、並進力及びトルクを計算する請求項１に記載のシミュレーション方法。

【請求項3】

弾性を有する解析対象物の挙動を、分子動力学法を用いてシミュレーションする装置であって、
シミュレーション条件が入力される入力部と、
前記入力部に入力されたシミュレーション条件に基づいて解析を行う処理部と
を備え、
前記処理部は、
前記解析対象物を仮想的な複数の粒子の集合体として表し、
タイムステップごとに各粒子に作用する力を求め、
各粒子に作用する力に基づいて粒子ごとに運動方程式を解いて、各粒子の速度及び位置を更新し、
粒子ごとに運動方程式を解く際に、
前記解析対象物を剛体とみなし、タイムステップごとに求めた各粒子に作用する力、各粒子の速度及び位置に基づいて、前記解析対象物に作用する力及び前記解析対象物の速度を計算し、
前記解析対象物の運動に対してＦＩＲＥを適用して、前記解析対象物に作用させるべき力を計算し、
前記解析対象物に作用させるべき力を複数の粒子に分配することにより、各粒子に追加的に作用させる力を求め、
各粒子に作用する力に、前記追加的に作用させる力を追加して、運動方程式を解くシミュレーション装置。

【請求項4】

前記処理部は、
前記解析対象物に作用する力として、並進力及びトルクと計算し、前記解析対象物の速度として、並進速度及び角速度を計算し、
前記解析対象物の運動に対してＦＩＲＥを適用する際に、前記解析対象物の並進運動及び回転運動のそれぞれに対してＦＩＲＥを適用し、
前記解析対象物に作用させるべき力として、並進力及びトルクを計算する請求項３に記載のシミュレーション装置。

【請求項5】

シミュレーション条件を取得する手順と、
弾性を有する解析対象物を仮想的な複数の粒子の集合体として表す手順と、
タイムステップごとに各粒子に作用する力を求める手順と、
各粒子に作用する力に基づいて粒子ごとに運動方程式を解いて、各粒子の速度及び位置を更新する手順と、
粒子ごとに運動方程式を解く際に、
前記解析対象物を剛体とみなし、タイムステップごとに求めた各粒子に作用する力、各粒子の速度及び位置に基づいて、前記解析対象物に作用する力及び前記解析対象物の速度を計算する手順と、
前記解析対象物の運動に対してＦＩＲＥを適用して、前記解析対象物に作用させるべき力を計算する手順と、
前記解析対象物に作用させるべき力を複数の粒子に分配することにより、各粒子に追加的に作用させる力を求める手順と、
各粒子に作用する力に、前記追加的に作用させる力を追加して、運動方程式を解く手順と
をコンピュータに実行させるプログラム。

【発明の詳細な説明】

【技術分野】

【0001】

本発明は、シミュレーション方法、シミュレーション装置、及びプログラムに関する。

【背景技術】

【0002】

複数の粒子の運動をシミュレーションする方法として、分子動力学法（ＭＤ法）、繰込み群分子動力学法（ＲＭＤ法）が知られている（特許文献１）。本明細書において、分子動力学法及び繰込み群分子動力学法を、単に「分子動力学法」ということとする。弾性を有する解析対象物を複数の粒子で表した解析モデルを生成し、粒子間に作用する力を定義し、これらの粒子の運動を、分子動力学法を用いて解析することにより、解析対象物の挙動をシミュレーションすることができる。

【0003】

分子動力学法を用いて粒子系の静解析を行う場合、求めるべき物理量の振動を短時間で収束させ、解を高速で得る手法として、ＦＩＲＥ（Fast Inertial Relaxation Engine）が知られている（非特許文献１）。

【先行技術文献】

【特許文献】

【0004】

【文献】特開２０２０－２０１８２８号公報

【非特許文献】

【0005】

【文献】Erik Bitzek rt. al., "Structural Relaxation Made Simple", Physical Review Letters 97, 170201 (2006)

【発明の概要】

【発明が解決しようとする課題】

【0006】

ＦＩＲＥは、単体の粒子に対して適用されるものであり、弾性を有する解析対象物を表す粒子の集合体に対して適用することはできない。本発明の目的は、弾性を有する解析対象物を複数の粒子の集合体として表し、粒子ごとの運動方程式を解いて解析対象物の挙動をシミュレーションする際に、解を高速で得ることが可能なシミュレーション方法、シミュレーション装置、及びプログラムを提供することである。

【課題を解決するための手段】

【0007】

本発明の一観点によると、
弾性を有する解析対象物を仮想的な複数の粒子の集合体として表し、
タイムステップごとに各粒子に作用する力を求め、
各粒子に作用する力に基づいて粒子ごとに運動方程式を解いて、各粒子の速度及び位置を更新するシミュレーション方法であって、
粒子ごとに運動方程式を解く際に、
前記解析対象物を剛体とみなし、タイムステップごとに求めた各粒子に作用する力、各粒子の速度及び位置に基づいて、前記解析対象物に作用する力及び前記解析対象物の速度を計算し、
前記解析対象物の運動に対してＦＩＲＥを適用して、前記解析対象物に作用させるべき力を計算し、
前記解析対象物に作用させるべき力を複数の粒子に分配することにより、各粒子に追加的に作用させる力を求め、
各粒子に作用する力に、前記追加的に作用させる力を追加して、運動方程式を解くシミュレーション方法が提供される。

【0008】

本発明の他の観点によると、
弾性を有する解析対象物の挙動を、分子動力学法を用いてシミュレーションする装置であって、
シミュレーション条件が入力される入力部と、
前記入力部に入力されたシミュレーション条件に基づいて解析を行う処理部と
を備え、
前記処理部は、
前記解析対象物を仮想的な複数の粒子の集合体として表し、
タイムステップごとに各粒子に作用する力を求め、
各粒子に作用する力に基づいて粒子ごとに運動方程式を解いて、各粒子の速度及び位置を更新し、
粒子ごとに運動方程式を解く際に、
前記解析対象物を剛体とみなし、タイムステップごとに求めた各粒子に作用する力、各粒子の速度及び位置に基づいて、前記解析対象物に作用する力及び前記解析対象物の速度を計算し、
前記解析対象物の運動に対してＦＩＲＥを適用して、前記解析対象物に作用させるべき力を計算し、
前記解析対象物に作用させるべき力を複数の粒子に分配することにより、各粒子に追加的に作用させる力を求め、
各粒子に作用する力に、前記追加的に作用させる力を追加して、運動方程式を解くシミュレーション装置が提供される。

【0009】

本発明のさらに他の観点によると、
シミュレーション条件を取得する手順と、
弾性を有する解析対象物を仮想的な複数の粒子の集合体として表す手順と、
タイムステップごとに各粒子に作用する力を求める手順と、
各粒子に作用する力に基づいて粒子ごとに運動方程式を解いて、各粒子の速度及び位置を更新する手順と、
粒子ごとに運動方程式を解く際に、
前記解析対象物を剛体とみなし、タイムステップごとに求めた各粒子に作用する力、各粒子の速度及び位置に基づいて、前記解析対象物に作用する力及び前記解析対象物の速度を計算する手順と、
前記解析対象物の運動に対してＦＩＲＥを適用して、前記解析対象物に作用させるべき力を計算する手順と、
前記解析対象物に作用させるべき力を複数の粒子に分配することにより、各粒子に追加的に作用させる力を求める手順と、
各粒子に作用する力に、前記追加的に作用させる力を追加して、運動方程式を解く手順と
をコンピュータに実行させるプログラムが提供される。

【発明の効果】

【0010】

弾性を有する解析対象物を複数の粒子の集合体として表し、粒子ごとの運動方程式を解いて解析対象物の挙動をシミュレーションする際に、解を高速で得ることができる。

【図面の簡単な説明】

【0011】

【図1】図１Ａは、弾性を有する解析対象物を複数の粒子の集合体で表した解析モデルを示す模式図であり、図１Ｂは、剛体とみなした解析対象物に定義される物理量を示す模式図である。

【図2】図２は、本実施例によるシミュレーション装置のブロック図である。

【図3】図３は、本実施例によるシミュレーション方法のフローチャートである。

【図4】図４は、回転体の一方の端面の一部の領域を固定した固定部材に加わるトルクの時間変化を示すグラフである。

【図5】図５は、回転体を剛体とみなしたときの角速度ωの時間変化を示すグラフである。

【図6】図６Ａ及び図６Ｂは、ＦＩＲＥを適用して解析する解析モデルの模式図である。

【発明を実施するための形態】

【0012】

［一般的なＦＩＲＥ］
実施例について説明する前に、図６を参照して、分子動力学法を用いたシミュレーションにおいて解を高速で得る一手法であるＦＩＲＥについて説明する。

【0013】

図６Ａ及び図６Ｂは、ＦＩＲＥを適用して解析する解析モデルの模式図である。解析空間に粒子２０が配置されている。粒子２０の速度をｖ、粒子２０に作用する力をｆと標記する。ＦＩＲＥを適用する際には、タイムステップごとに、以下の判定値Ｐ（ｔ）を計算する。

【数1】

【0014】

判定値Ｐ（ｔ）が負またはゼロのとき、すなわち、図６Ａに示すように速度ｖと力ｆとのなす角度が９０°以上のとき、次状態における粒子２０の速度ｖ（ｔ＋Δｔ）を、以下のように決定する。

【数2】

言い換えると、力が粒子２０を減速させる方向に作用する場合、その粒子２０の次状態の速度をゼロにする。

【0015】

判定値Ｐ（ｔ）が正のとき、すなわち、図６Ｂに示すように速度ｖと力ｆとのなす角度が９０°未満のとき、粒子２０について以下の運動方程式を解くことにより、次状態の速度ｖ（ｔ＋Δｔ）を求める。

【数3】

ここで、ｍは粒子２０の質量であり、γは減衰定数である。ベクトルに付されたハットは、単位ベクトルであることを意味する。式（３）は、粒子２０を、力ｆの方向に、より加速させることを意味する。

【0016】

［参考例によるＦＩＲＥ］
次に、弾性を有する解析対象物を粒子の集合体で表した解析モデルに適用する参考例によるＦＩＲＥについて説明する。

【0017】

参考例によるＦＩＲＥにおいては、判定値Ｐ（ｔ）を以下の式で定義する。

【数4】

ここで、ｆ_ｉはｉ番目の粒子に作用する力であり、ｖ_ｉは、ｉ番目の粒子の速度である。Ｎは粒子の個数であり、式（４）の右辺のΣは、解析対象物を構成するすべての粒子について合計することを意味する。

【0018】

判定値Ｐ（ｔ）が負またはゼロのとき、次状態におけるｉ番目の粒子２０の速度ｖ_ｉ（ｔ＋Δｔ）を、以下のように決定する。

【数5】

【0019】

判定値Ｐ（ｔ）が正のとき、ｉ番目の粒子２０について以下の運動方程式を解くことにより、次状態の速度ｖ_ｉ（ｔ＋Δｔ）を求める。

【数6】

ここで、ｍ_ｉは、ｉ番目の粒子の質量である。

【0020】

上記参考例によるＦＩＲＥは、以下の課題を有する。
参考例によるＦＩＲＥにおいては、式（４）に示すように、粒子のそれぞれに作用する力と粒子の速度との内積の合計値を判定値Ｐ（ｔ）として用いている。このため、解析対象物の重心速度がゼロであっても、内部振動によって判定値Ｐ（ｔ）がゼロではなく、ある値を持つ。解析対象物の温度が高い場合や、解析対象物の速度が収束に近づくと、判定値Ｐ（ｔ）が負という条件が度々満たされてしまう。この条件が満たされると、各粒子の速度がゼロになるため、収束が遅くなってしまう。

【0021】

さらに、判定値Ｐ（ｔ）が正のとき、式（６）に示すように、各粒子に対して、力の方向に向かって加速する力を追加的に付与した運動方程式を解くことになる。減衰定数γ（ｔ）は、粒子２０の各々に設定する必要がある。粒子２０の間で、質量や周囲との相互作用が異なるため、減衰定数γの設定が困難であり、計算が破綻する場合がある。

【0022】

［実施例］
次に、図１～図５を参照して、一実施例によるシミュレーション方法、シミュレーション装置、及びプログラムについて説明する。

【0023】

図１Ａは、弾性を有する解析対象物１０を複数の粒子１１の集合体で表した解析モデルを示す模式図である。なお、図１Ａでは、解析対象物１０を構成する粒子１１の一部のみを表している。解析空間にｘｙｚ直交座標系を定義し、ｉ番目の粒子１１の位置ベクトルをｒ_ｉと標記し、ｉ番目の粒子１１の速度をｖ_ｉと標記する。

【0024】

実施例においては、弾性を有する解析対象物１０を一つの剛体と仮定し、剛体とみなした解析対象物１０に対して物理量を定義する。

【0025】

図１Ｂは、剛体とみなした解析対象物１０に定義される物理量を示す模式図である。解析対象物１０の重心の位置ベクトルをＲと標記する。解析対象物１０の重心に作用する並進力をＦと標記し、重心を回転中心とした場合のトルクをＴと標記する。解析対象物１０の重心の並進速度をＶと標記し、角速度をωと標記する。例えば、解析対象物が、金属部材のように極わずかの変形しかしないような場合は、解析対象物を剛体とみなして、解析対象物に加わる並進力Ｆ、トルクＴ、並進速度Ｖ、及び角速度ωを定義することができる。

【0026】

図２は、本実施例によるシミュレーション装置のブロック図である。本実施例によるシミュレーション装置は、入力部４０、処理部４１、出力部４２、及び記憶部４３を含む。入力部４０から処理部４１にシミュレーション条件等が入力される。さらに、ユーザから入力部４０に各種指令（コマンド）等が入力される。入力部４０は、例えば通信装置、リムーバブルメディア読取装置、キーボード、ポインティングデバイス等で構成される。

【0027】

処理部４１は、入力されたシミュレーション条件及び指令に基づいて分子動力学法によるシミュレーションを実行する。さらに、シミュレーション結果を出力部４２に出力する。シミュレーション結果には、例えば、シミュレーション対象物を構成する粒子１１（図１Ａ）の位置の時間変化等を表す情報が含まれる。処理部４１は、例えばコンピュータの中央処理ユニット（ＣＰＵ）を含む。分子動力学法によるシミュレーションをコンピュータに実行させるためのプログラムが、記憶部４３に記憶されている。出力部４２は、通信装置、リムーバブルメディア書込み装置、ディスプレイ等を含む。

【0028】

図３は、本実施例によるシミュレーション方法のフローチャートである。処理部４１（図２）が記憶部４３（図２）に格納されているプログラムを実行することにより、シミュレーション方法の各ステップの処理が実行される。

【0029】

まず、処理部４１は、入力部４０から入力されたシミュレーション条件に基づいて、弾性を有する解析対象物１０（図１Ａ）を複数の粒子の集合体で表した解析モデルを生成する（ステップＳ１）。シミュレーション条件には、解析対象物１０の形状、解析対象物１０の物性値に依存する情報、シミュレーションの初期条件及び境界条件、時間刻み幅、解析終了条件等が含まれる。解析対象物１０の物性値に依存する情報には、解析モデルにおける粒子間相互作用ポテンシャルを定義する情報が含まれる。

【0030】

次に、粒子１１（図１Ａ）ごとに運動方程式を解いて、粒子１１のそれぞれの速度及び位置を１タイムステップ分更新する。

【0031】

次に、解析対象物１０（図１Ｂ）を剛体とみなして、粒子１１の位置及び速度に基づいて、解析対象物１０に作用する力、及び解析対象物１０の速度を計算する（ステップＳ３）。解析対象物１０の速度には、重心の並進速度と、重心を中心とした角速度とが含まれる。

【0032】

解析対象物１０に作用する並進力Ｆは、以下の式で計算する。

【数7】

すなわち、解析対象物１０に作用する並進力Ｆは、粒子１１のそれぞれに作用する力ｆ_ｉのベクトル和である。

【0033】

解析対象物１０の並進速度Ｖは、以下の式で計算する。

【数8】

ここで、Ｍは、解析対象物１０を構成するすべての粒子１１の質量の和である。解析対象物１０の並進速度Ｖは、解析対象物１０を構成する粒子１１の速度ｖ_ｉを、粒子１１の質量ｍ_ｉで重みづけして平均した値である。

【0034】

解析対象物１０に作用するトルクＴ＝（Ｔ_ｘ，Ｔ_ｙ，Ｔ_ｚ）は、以下の式で計算する。

【数9】

ここで、Ｔ_ｘ，Ｔ_ｙ，Ｔ_ｚは、それぞれ解析対象物１０に作用するトルクＴのｘ成分、ｙ成分、及びｚ成分である。ｘ_ｉ、ｙ_ｉ、ｚ_ｉは、それぞれ解析対象物１０の重心を原点としたときのｉ番目の粒子１１のｘ座標、ｙ座標、及びｚ座標である。ｆ_ｉ，ｘ、ｆ_ｉ，ｙ、ｆ_ｉ，ｚは、それぞれｉ番目の粒子１１に作用する力ｆ_ｉのｘ成分、ｙ成分、及びｚ成分である。

【0035】

解析対象物１０の角速度ω＝（ω_ｘ，ω_ｙ，ω_ｚ）は、以下の式で計算する。

【数10】

ここで、ω_ｘ、ω_ｙ、ω_ｚは、それぞれ角速度ωのｘ成分、ｙ成分、ｚ成分である。ｖ_ｉ，ｘ、ｖ_ｉ，ｙ、ｖ_ｉ，ｚは、それぞれｉ番目の粒子１１の速度のｘ成分、ｙ成分、及びｚ成分である。

【0036】

Ｉ_ｘ、Ｉ_ｙ、Ｉ_ｚは、それぞれｘ軸、ｙ軸、ｚ軸の周りの慣性モーメントであり、以下の式で定義される。

【数11】

【0037】

ステップＳ３（図３）の後、処理部４１（図２）は、剛体とみなした解析対象物１０の運動を規定する物理量の振動を緩和させるように、解析対象物１０の運動に対してＦＩＲＥを適用する。これにより、剛体とみなした解析対象物１０に作用させるべき力を計算する。より具体的には、解析対象物１０に作用させるべき並進力とトルクとを計算する（ステップＳ４）。

【0038】

以下、ステップＳ４の処理について詳細に説明する。
解析対象物１０の運動を、並進運動と回転運動とに分離し、それぞれに対してＦＩＲＥを適用する。まず、並進運動に対して適用するＦＩＲＥについて説明する。

【0039】

並進運動に対して適用するＦＩＲＥの判定値Ｐ_ｔ（ｔ）を、以下のように定義する。

【数12】

これは、個別の粒子に適用するＦＩＲＥの判定値（式（１））における粒子に作用する力ｆを解析対象物１０に作用する並進力Ｆに置き換え、粒子の速度ｖを解析対象物１０の並進速度Ｖに置き換えたものである。

【0040】

解析対象物１０の並進運動にＦＩＲＥを適用する場合には、判定値Ｐ_ｔ（ｔ）がゼロ、または負のとき、次状態における解析対象物１０の並進速度Ｖ（ｔ＋Δｔ）を、以下のように決定するとよい。

【数13】

式（１３）は、粒子の速度ｖをゼロにする式（２）に対応する。

【0041】

判定値Ｐ_ｔ（ｔ）が正のとき、解析対象物１０について以下の運動方程式を解くことにより、次状態の並進速度Ｖ（ｔ＋Δｔ）を求めるとよい。

【数14】

式（１４）は、粒子に対する運動方程式（式（３））に対応する。ここで、αは減衰定数である。

【0042】

回転運動に対して適用するＦＩＲＥの判定値Ｐ_ｒ（ｔ）を、以下のように定義する。

【数15】

式（１５）は、並進運動に対する式（１２）の並進力Ｆ及び並進速度Ｖを、それぞれトルクＴ及び角速度ωに置き換えたものである。

【0043】

解析対象物１０の回転運動にＦＩＲＥを適用する場合には、判定値Ｐ_ｒ（ｔ）がゼロ、または負のとき、次状態における解析対象物１０の角速度ω（ｔ＋Δｔ）を、以下のように決定するとよい。

【数16】

式（１６）は、解析対象物１０の並進速度Ｖをゼロにする式（１３）に対応する。

【0044】

判定値Ｐ_ｒ（ｔ）が正のとき、解析対象物１０について以下の運動方程式を解くことにより、次状態の角速度ω（ｔ＋Δｔ）を求めるとよい。

【数17】

式（１７）は、並進運動に対する運動方程式（式（１４））に対応する。ここで、βは減衰定数である。Ｉは、解析対象物１０の慣性モーメントであり、式（１１）で定義される。ここで、Ｔ（ｔ）／Ｉは、ベクトル（Ｔｘ／Ｉｘ、Ｔｙ／Ｉｙ，Ｔｚ／Ｉｚ）を意味する。

【0045】

式（１４）から、ＦＩＲＥによって解析対象物１０に追加的に作用させるべき並進力Ｆ_ａｄｄは、以下の式で記述される。

【数18】

【0046】

式（１７）から、ＦＩＲＥによって解析対象物１０に追加的に作用させるべきトルクＴ_ａｄｄのｘ成分は、以下の式で記述される。トルクＴ_ａｄｄのｙ成分及びｚ成分についても同様である。

【数19】

【0047】

ステップＳ４の後、処理部４１は、解析対象物１０に追加的に作用させる並進力Ｆ_ａｄｄ及びトルクＴ_ａｄｄを、各粒子１１に分配し、各粒子１１に追加的に作用させる力を計算する（ステップＳ５）。

【0048】

追加的に作用させる並進力Ｆ_ａｄｄは、以下の式に基づいて、各粒子１１に分配する。

【数20】

【0049】

追加的に作用させるトルクＴ_ａｄｄを、各粒子１１に分配したときのｉ番目の粒子１１に作用させるべき力のｘ成分ｔ_ｉ，ｘ、ｙ成分ｔ_ｉ，ｙ、ｚ成分ｔ_ｉ，ｚは、以下の式で記述される。

【0050】

【数21】

【0051】

解析対象物１０の並進運動に対する判定値Ｐ_ｔ（ｔ）がゼロまたは負のとき、粒子１１の各々の速度ｖ_ｉから現時点の解析対象物１０の並進速度Ｖを減じることにより、式（８）で定義される並進速度Ｖをゼロにする。回転運動に対する判定値Ｐ_ｒ（ｔ）がゼロまたは負のとき、粒子１１の各々の速度ｖ_ｉから、解析対象物１０の回転速度に相当する粒子１１の各々の速度ｒ_ｉωを減じることにより、式（１０）のω_ｘ、ω_ｙ、及びω_ｚをゼロにする。

【0052】

ステップＳ２からステップＳ５までの手順を、解析終了条件が満たされるまで繰り返す（ステップＳ６）。ステップＳ２で運動方程式を解くときには、ステップＳ５で求めた追加的に作用させる力を、粒子１１のそれぞれに追加する。

【0053】

次に、上記実施例の優れた効果について説明する。
上記実施例では、ステップＳ４（図３）において、剛体とみなした解析対象物１０の運動に対してＦＩＲＥを適用するため、解析対象物１０全体としてつり合い状態に達するまでの計算時間を短縮することができる。

【0054】

上記実施例の優れた効果を確認するために、上記実施例によるシミュレーション方法を用いて円柱状の回転体の静解析を行った。次に、図４及び図５を参照して、このシミュレーション結果について説明する。シミュレーションにおいて、回転体の一方の端面の一部の領域を固定部材に取り付けて固定し、回転体の外周面に一定のトルクを印加した。

【0055】

図４は、回転体の一方の端面の一部の領域を固定した固定部材に加わるトルクの時間変化を示すグラフである。横軸は、解析開始から一定時間経過したときからの経過時間を表し、縦軸は、印加したトルクで規格化したトルクを表す。図４の太い実線は、上記実施例によるシミュレーション方法を用いて求めたトルクを示し、細い実線は、式（４）～式（６）を用いた参考例によるシミュレーション方法を用いて求めたトルクを示し、破線は、ＦＩＲＥを適用しないでシミュレーションして求めたトルクを示す。つり合い状態では、固定部材に加わるトルクは、印加したトルクに等しくなる。すなわち、規格化したトルクが１になる。

【0056】

ＦＩＲＥを適用しない場合は、トルクの計算値が、つり合いの値を中心として振動する。参考例によるトルクの計算値は、オーバシュートした後、緩やかにつり合いの値に近づく。これに対して上記実施例によるトルクの計算値は、ほぼつり合いの値でトルクの計算値の上昇が緩やかになり、短時間でほぼつり合いの状態に達する。

【0057】

図５は、回転体を剛体とみなしたときの角速度ωの時間変化を示すグラフである。横軸は、解析開始から一定時間経過したときからの経過時間を表し、縦軸は、最大値を１として規格化した角速度を表す。図４の太い実線は、上記実施例によるシミュレーション方法を用いて求めた角速度を示し、細い実線は、式（４）～式（６）を用いた参考例によるシミュレーション方法を用いて求めた角速度を示し、破線は、ＦＩＲＥを適用しないでシミュレーションして求めた角速度を示す。つり合い状態では、角速度がゼロになる。

【0058】

ＦＩＲＥを適用しない場合は、角速度の計算値がゼロを中心として振動する。参考例による角速度の計算値は、最大になった後に低下し、つり合い状態の直前で不連続にゼロになっている。角速度の計算値が不連続にゼロになる時点は、式（４）の判定値Ｐ（ｔ）がゼロまたは負になったことを示している。これに対して上記実施例による角速度の計算値は、角速度の計算値がほぼ最大値になった直後に、角速度の計算値が不連続にゼロになっている。角速度の計算値が不連続にゼロになる時点は、式（１２）の判定値Ｐ_ｔ（ｔ）または式（１５）の判定値Ｐ_ｒ（ｔ）がゼロまたは負になったことを示している。なお、本シミュレーションにおいては、判定値Ｐ_ｒ（ｔ）がゼロまたは負になっている。

【0059】

図４及び図５に示したように、上記実施例によるシミュレーション方法を用いることにより、ＦＩＲＥを適用しない場合、及び参考例によるシミュレーション方法を用いた場合と比べて、つり合い状態に到達するまでの時間が短縮されていることがわかる。

【0060】

次に、上記実施例の変形例について説明する。
上記実施例では、ステップＳ３（図３）において、剛体とみなした解析対象物１０の速度として、並進速度と角速度とを求めたが、並進速度及び角速度の一方のみを求め、解析対象物１０の並進運動及び回転運動の一方のみに対してＦＩＲＥを適用してもよい。

【0061】

上述の実施例は例示であり、本発明は上述の実施例に制限されるものではない。例えば、種々の変更、改良、組み合わせ等が可能なことは当業者に自明であろう。

【符号の説明】

【0062】

１０ 解析対象物
１１、２０ 粒子
４０ 入力部
４１ 処理部
４２ 出力部
４３ 記憶部

【図1】

【図2】

【図3】

【図4】

【図5】

【図6】