JP7239309B2

JP7239309B2 - シミュレーション装置及びプログラム

Info

Publication number: JP7239309B2
Application number: JP2018231799A
Authority: JP
Inventors: 公則坂井
Original assignee: Sumitomo Heavy Industries Ltd
Current assignee: Sumitomo Heavy Industries Ltd
Priority date: 2018-12-11
Filing date: 2018-12-11
Publication date: 2023-03-14
Anticipated expiration: 2038-12-11
Also published as: JP2020095400A; US20200184029A1

Description

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

分子動力学法、繰り込み群分子動力学法等の動的陽解法を用いる機構構造解析では、計算を発散させないためにメッシュサイズ（粒子間距離）に応じた時間刻み幅以下で時間発展を行うことが好ましい。メッシュサイズが小さいほど（粒子間距離が短いほど）、好ましい時間刻み幅は小さくなる。解析対象の系のメッシュサイズが場所によって異なっている場合、系内の最も小さなメッシュサイズに合わせて時間刻み幅を小さくしなければならない。このため、演算量が増え、解析時間が長くなってしまう。

下記の特許文献１に、サルコメア内の分子の運動と、筋肉の連続体モデルを連成させて解析を行う生体シミュレーション装置が開示されている。サルコメアモデルのモンテカルロシミュレーション処理の時間刻み幅と、連続体モデルの有限要素解析の時間刻み幅とは異なっている。

特開２０１４－１４９６４９号公報

従来例において、動的陽解法を用いた連続体モデルの有限要素解析の時間刻み幅は連続体モデルのメッシュサイズに依らず一定である。このため、最も小さなメッシュサイズに対応して時間刻み幅を短くしなければならず、解析時間を短縮することが困難である。

本発明の目的は、解析時間を短縮することが可能なシミュレーション装置及びプログラムを提供することである。

本発明の一観点によると、
シミュレーション対象物の形状を定義する形状定義データ、及びシミュレーション対象物をメッシュ分割するときのメッシュサイズに関するメッシュ分割条件が入力される入力部と、
前記入力部に入力された前記形状定義データ及び前記メッシュ分割条件に基づいて、シミュレーション対象物をメッシュ分割し、生成されたメッシュの節点に粒子を配置し、前記粒子の間の相互作用ポテンシャルに基づいて前記粒子の状態を時間発展させる処理部と
を有し、
前記処理部は、
前記粒子の状態を時間発展させる時間刻み幅として、長さの異なる複数の時間刻み幅を定義し、かつ複数の前記時間刻み幅を昇順に並べたとき、前記時間刻み幅のそれぞれが、１つ前の時間刻み幅の整数倍になるように前記時間刻み幅を定義し、
前記シミュレーション対象物の弾性の程度を表す物性値、前記シミュレーション対象物の密度、及び生成されたメッシュの形状に基づいて前記粒子の間の相互作用ポテンシャル及び前記粒子の質量を決定し、
前記粒子の間の相互作用ポテンシャル及び前記粒子の質量から求まる固有振動の周期を求め、
前記粒子のそれぞれについて、当該粒子が含まれる粒子ペアの固有振動の周期に基づいて、複数の前記時間刻み幅から１つの時間刻み幅を決定し、前記粒子の状態を時間発展させる際に、決定された時間刻み幅を用いて時間発展させるシミュレーション装置が提供される。

本発明の他の観点によると、
シミュレーション対象物の形状を定義する形状定義データ、及びシミュレーション対象物をメッシュ分割するときのメッシュサイズに関するメッシュ分割条件を取得する機能と、
前記形状定義データ及び前記メッシュ分割条件に基づいて、シミュレーション対象物をメッシュ分割し、生成されたメッシュの節点に粒子を配置し、前記粒子の間の相互作用ポテンシャルに基づいて前記粒子の状態を時間発展させる機能と、
前記粒子の状態を時間発展させる時間刻み幅として、長さの異なる複数の時間刻み幅を定義し、かつ複数の前記時間刻み幅を昇順に並べたとき、前記時間刻み幅のそれぞれが、１つ前の時間刻み幅の整数倍になるように前記時間刻み幅を定義する機能と、
前記シミュレーション対象物の弾性の程度を表す物性値、前記シミュレーション対象物の密度、及び生成されたメッシュの形状に基づいて前記粒子の間の相互作用ポテンシャル及び前記粒子の質量を決定する機能と、
前記粒子の間の相互作用ポテンシャル及び前記粒子の質量から求まる固有振動の周期を求める機能と、
前記粒子のそれぞれについて、当該粒子が含まれる粒子ペアの固有振動の周期に基づいて、複数の前記時間刻み幅から１つの時間刻み幅を決定し、前記粒子の状態を時間発展させる際に、決定された時間刻み幅を用いて時間発展させる機能と
をコンピュータに実現させるプログラムが提供される。

粒子の状態を時間発展させる時間刻み幅を、前記粒子の間の距離に応じて異ならせることにより、解析時間を短縮することが可能になる。

図１は、実施例によるシミュレーション装置のブロック図である。図２は、本実施例によるシミュレーション方法のフローチャートである。図３は、固有周期Ｐと、時間刻み幅Δｔとの関係の一例を示すグラフである。図４は、第１ループＬ１～第４ループＬ４を実行するタイミングチャートである。図５は、シミュレーション対象物の形状から生成されたメッシュの複数の節点に配置した複数の粒子の一例を示す図である。図６は、ステップＳ６（図２）の詳細な手順を示すフローチャートである。図７は、第１ループＬ１の処理（図６のステップＳ６７）の手順を示すフローチャートである。図８Ａは、ステップＳ６６（図６）の処理の手順を示すフローチャートであり、図８Ｂは、ステップＳ６４（図６）の処理の手順を示すフローチャートである。図９は、ステップＳ６２（図６）の処理の手順を示すフローチャートである。図１０Ａは、シミュレーション対象の２つの物体の形状モデルを示す図であり、図１０Ｂは、２つの物体が接触した状態の形状モデルを示す図である。図１１Ａは、検証モデルの概略斜視図であり、図１１Ｂは、検証モデルの固定されていない方の端部の位置の時間変化のシミュレーション結果を示すグラフである。図１２は、物体内相互作用計算、及び位置、速度計算に要した１タイムステップあたりの計算時間を、実施例と比較例とで対比させて示すグラフである。

図１～図１２を参照して実施例によるシミュレーション装置、シミュレーション方法について説明する。

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

処理部２１は、入力されたシミュレーション条件及び指令に基づいてシミュレーション対象物の形状をメッシュ分割し、メッシュの節点に仮想的に粒子を配置して分子動力学法によるシミュレーションを行う。さらに、シミュレーション結果を出力部２２に出力する。シミュレーション結果には、シミュレーション対象物の形状の時間的変化を表す情報が含まれる。処理部２１は、例えばコンピュータを含み、くりこみ群分子動力学法によるシミュレーションをコンピュータに実行させるためのプログラムが記憶部２３に記憶されている。出力部２２は、通信装置、リムーバブルメディア書込み装置、ディスプレイ等を含む。

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

まず、処理部２１は、シミュレーション対象物の形状を定義する形状定義データ、シミュレーション対象物の物性値（密度、弾性の程度を表す物性値等）、シミュレーションの初期条件、境界条件、メッシュ分割条件等を取得する（ステップＳ１）。これらのデータは入力部２０に入力され、処理部２１が入力部２０から取得する。形状定義データを取得すると、処理部２１は形状定義データで定義される形状をメッシュ分割条件に基づいてメッシュ分割することにより、形状モデルを生成する（ステップＳ２）。メッシュ分割処理には、例えば公知のメッシュ分割アルゴリズムを使用することができる。本実施例では、シミュレーション対象物をテトラメッシュに分割する。

メッシュが生成されると、処理部２１は、メッシュの節点に仮想的に粒子を配置し、各粒子に質量を付与し、粒子間の相互作用ポテンシャルを決定する（ステップＳ３）。粒子の質量は、例えば、シミュレーション対象物の密度と、粒子の分布とに基づいて、形状モデルの密度がシミュレーション対象物の密度と同一になるように決定する。例えば、シミュレーション対象物の密度が均一である場合、粒子の分布密度が高い領域では粒子の質量が相対的に小さくなり、粒子の分布密度が低い領域では粒子の質量が相対的に大きくなる。

粒子間の相互作用ポテンシャルとして、例えばバネマスモデルのポテンシャルを用いる。バネ定数の決定方法は、例えば特開２００９－３７３３４号公報に詳しく説明されている。ここでは、バネ定数の決定方法について簡単に説明する。

まず、テトラメッシュの節点に配置された粒子のボロノイ多面体解析を行う。ボロノイ多面体を構成する複数の界面のうち粒子ｉと粒子ｊとを両端とする線分を横切る界面の面積をＳ_ｉｊで表す。シミュレーション対象物のヤング率と面積Ｓ_ｉｊとからバネ定数ｋ_ｉｊを決定することができる。

次に、バネで相互に接続されている２つの粒子（粒子ペア）ごとに、固有振動の周期（固有周期）を算出する（ステップＳ４）。一般的に、粒子ペアを構成する２つの粒子の質量は同一ではないため、一方の粒子を固定した状態における固有周期と、他方の粒子を固定した状態における固有周期とは異なる。粒子のペアの固有周期として、例えば、２つの固有周期のうち短い方の固有周期の値を採用する。

処理部２１は、ステップＳ４で算出された固有周期に基づいて、粒子ペアごとに、粒子の状態の時間発展の演算を行う際の時間刻み幅を決定する（ステップＳ５）。固有周期が短いほど、一定時間に粒子の状態が大きく変化する。粒子の状態の変化を精度よく解析するために、固有周期が短いほど時間刻み幅を短くするとよい。例えば、時間刻み幅は、固有周期の１／２０以上１／１０以下の範囲から選択するとよい。１つの形状モデルに対して、複数の時間刻み幅が定義される。

時間刻み幅が決定されると、処理部２１は、決定された時間刻み幅で粒子の状態の変化を解析する（ステップＳ６）。具体的には、運動方程式を解くことにより粒子の速度と位置の時間変化を求める。解析が終了すると、解析結果を出力部２２（図１）に出力する（ステップＳ７）。

図３は、固有周期Ｐと、時間刻み幅Δｔとの関係の一例を示すグラフである。横軸は固有周期Ｐを表し、縦軸は時間刻み幅Δｔを表す。例えば、固有周期ＰがＰ１未満、Ｐ１以上Ｐ２未満、Ｐ２以上Ｐ３未満、及びＰ３以上のとき、それぞれ時間刻み幅Δｔを、Δｔ_１、Δｔ_２、Δｔ_３、及びΔｔ_４とする。ここで、Δｔ_１＜Δｔ_２＜Δｔ_３＜Δｔ_４の大小関係が成立する。さらに、Δｔ_４はΔｔ_３の整数倍であり、Δｔ_３はΔｔ_２の整数倍であり、Δｔ_２はΔｔ_１の整数倍である。

固有周期ＰがＰ３以上の粒子ペアについては、時間刻み幅Δｔ_４で第４ループＬ４の処理を実行することにより、粒子間に働く力の計算、及び速度の時間発展を行う。固有周期ＰがＰ２以上Ｐ３未満の粒子ペアについては、時間刻み幅Δｔ_３で第３ループＬ３の処理を実行する。固有周期ＰがＰ１以上Ｐ２未満の粒子ペアについては、時間刻み幅Δｔ_２で第２ループＬ２の処理を実行する。固有周期ＰがＰ１未満の粒子ペアについては、時間刻み幅Δｔ_１で第１ループＬ１の処理を実行する。

図４は、第１ループＬ１～第４ループＬ４を実行するタイミングチャートである。第１ループＬ１の処理は、時間刻み幅Δｔ_１ごとに実行される。言い換えると、第１ループＬ１の処理を１回実行すると、時間刻み幅Δｔ_１だけ時間発展する。この時間発展の累積が時間刻み幅Δｔ_２だけ増加するごとに、第２ループＬ２の処理が実行される。時間発展の累積が時間刻み幅Δｔ_３だけ増加するごとに、第３ループＬ３の処理が実行される。時間発展の累積が時間刻み幅Δｔ_４だけ増加するごとに、第４ループＬ４の処理が実行される。時間刻み幅Δｔ_１だけ時間発展させる１回の処理を１つのタイムステップということとする。

図５は、シミュレーション対象物の形状から生成されたメッシュの複数の節点に配置した複数の粒子の一例を示す図である。図５に示した複数の粒子から複数の粒子ペアが定義され、複数の粒子ペアにそれぞれ固有周期に基づいて第１ループＬ１～第４ループＬ４のいずれかが対応付けられる。例えば、粒子ペアＡＢ、ＢＣに第１ループＬ１が対応付けられている。粒子ペアＢＤ、ＢＥ、ＣＥに第２ループＬ２が対応付けられている。粒子ペアＤＥ、ＤＦに第３ループＬ３が対応付けられている。粒子ペアＥＦ、ＥＧに第４ループＬ４が対応付けられている。

第１ループＬ１（図４）のみを実行するタイムステップでは、粒子ペアＡＢ及びＢＣについて、粒子の速度の時間発展、及び力の計算を行う。第１ループＬ１及び第２ループＬ２（図４）を実行するタイムステップでは、さらに、粒子ペアＢＤ、ＢＥ、及びＣＥについても、粒子の速度の時間発展、及び力の計算を行う。第１ループＬ１、第２ループＬ２、及び第３ループＬ３（図４）を実行するタイムステップでは、さらに、粒子ペアＤＥ、及びＤＦについても、粒子の速度の時間発展、及び力の計算を行う。第１ループＬ１～第４ループＬ４の全てを実行するタイムステップでは、さらに粒子ペアＥＦ、及びＥＧについても、粒子の速度の時間発展、及び力の計算を行う。

図６は、粒子の状態の変化を解析する処理（図２のステップＳ６）の詳細な手順を示すフローチャートである。以下の説明において、適宜、図５を参照する。

シミュレーションの演算終了条件が満たされるまで、第４ループＬ４の処理を実行する。演算終了条件として、例えばタイムステップ数、演算時間等が与えられる。演算終了条件が満たされない場合には、まず、第４ループＬ４の処理対象となる粒子ペアＥＦ、及びＥＧ（図５）に含まれる粒子Ｅ、Ｆ、Ｇについて運動方程式を解き、時間刻み幅Δｔ_４の１／２だけ各粒子の速度を時間発展させる（ステップＳ６１）。

ステップＳ６１を実行した後、第３ループＬ３の処理対象となる粒子ペアＤＥ、及びＤＦ（図５）に含まれる粒子Ｄ、Ｅ、Ｆについて運動方程式を解き、時間刻み幅Δｔ_３の１／２だけ各粒子の速度を時間発展させる（ステップＳ６３）。ここで、ステップＳ６１で速度を時間発展させた粒子、例えば粒子Ｅ、Ｆについては、時間発展後の最新の速度を、さらに時間発展させる。以降の速度の時間発展の処理においても同様に、最新の速度をさらに時間発展させる。

ステップＳ６３を実行した後、第２ループＬ２の処理対象となる粒子ペアＢＤ、ＢＥ、及びＣＥ（図５）に含まれる粒子Ｂ、Ｃ、Ｄ、Ｅについて運動方程式を解き、時間刻み幅Δｔ_２の１／２だけ各粒子の速度を時間発展させる（ステップＳ６５）。

ステップＳ６５を実行した後、全ての粒子について第１ループＬ１の処理を実行する（ステップＳ６７）。第１ループＬ１の処理の詳細については、後に図７を参照して説明する。

第１ループＬ１の処理（タイムステップ）を複数回繰り返して、ステップＳ６５以降の時間発展の累積値がΔｔ_２だけ増加すると、第２ループＬ２の後処理を実行する（ステップＳ６６）。第２ループＬ２の後処理（ステップＳ６６）では、時間刻み幅Δｔ_２の１／２だけ速度をさらに時間発展させる。第２ループＬ２の後処理（ステップＳ６６）については、後に図８Ａを参照して詳細に説明する。このように、速度の時間発展を時間刻み幅の１／２ずつ実行する方法は、速度ベルレ法と呼ばれる。

第２ループＬ２の処理を複数回繰り返して、ステップＳ６３以降の時間発展の累積値がΔｔ_３だけ増加すると、第３ループＬ３の後処理を実行する（ステップＳ６４）。第３ループＬ３の後処理（ステップＳ６４）では、時間刻み幅Δｔ_３の１／２だけ速度をさらに時間発展させる。第３ループＬ３の後処理（ステップＳ６４）については、後に図８Ｂを参照して詳細に説明する。

第３ループＬ３の処理を複数回繰り返した後、ステップＳ６１以降の時間発展の累積値がΔｔ_４だけ増加すると、第４ループＬ４の後処理を実行する（ステップＳ６２）。第４ループＬ４の後処理（ステップＳ６２）では、時間刻み幅Δｔ_４の１／２だけ速度をさらに時間発展させる。第４ループＬ４の後処理（ステップＳ６２）については、後に図９を参照して詳細に説明する。

図７は、第１ループＬ１の処理（図６のステップＳ６７）の手順を示すフローチャートである。まず、第１ループＬ１の処理対象となる粒子ペアＡＢ、及びＢＣ（図５）に含まれる粒子Ａ、Ｂ、Ｃについて運動方程式を解き、時間刻み幅Δｔ_１の１／２だけ速度を時間発展させる（ステップＳ６７１）。次に、全ての粒子の位置を、各粒子の最新の速度に基づいて時間刻み幅Δｔ_１だけ時間発展させる（ステップＳ６７２）。

全ての粒子の位置を時間発展させた後、物体間の接触の有無を判定する（ステップＳ６７３）。例えば、シミュレーション対象物が複数の物体で構成されている場合、同一の物体内の粒子の間には、ステップＳ３（図２）で決定した相互作用ポテンシャルに基づく力が作用する。異なる物体の粒子の間には、物体が接触していない状態では力は作用しない。ところが、２つの物体が接触すると、物体の接触箇所に位置する粒子の間で力を及ぼし合う。ステップＳ６７３で物体間の接触があったと判定された場合には、次のタイムステップにおいて、物体の接触箇所に位置する粒子に対して、異なる物体の粒子間に作用する相互作用ポテンシャルを導入して速度及び位置を時間発展させる。物体間の接触の判定、及び異なる物体の粒子間の相互作用ポテンシャルについては、後に図１１Ａ及び図１１Ｂを参照して説明する。

次に、粒子の時間発展後の位置に基づいて、粒子に作用する力を計算する（ステップＳ６７４）。このとき、粒子間の相互作用ポテンシャルに基づいて生じる力については、第１ループＬ１の処理対象の粒子ペアＡＢ及びＢＣ（図５）についてのみ計算する。また、重力や境界条件として与えられた荷重等の外力が作用する粒子については、これらの外力を計算する。第２ループＬ２～第４ループＬ４の処理対象の粒子ペアについては、相互作用ポテンシャルに基づく力の計算は行わない。各粒子について運動方程式を解く際には、当該粒子に作用する力を合成した合成力に基づいて計算を行う。

粒子に作用する力が求まると、第１ループＬ１の処理対象となる粒子ペアＡＢ及びＢＣを構成する粒子Ａ、Ｂ、Ｃの各々について運動方程式を解き、時間刻み幅Δｔ_１の１／２だけ粒子の速度を時間発展させる（ステップＳ６７５）。ステップＳ６７１とステップＳ６７５とで、それぞれ粒子の速度が時間刻み幅Δｔ_１の１／２ずつ時間発展するため、粒子の速度は時間刻み幅Δｔ_１だけ時間発展することになる。

図８Ａは、ステップＳ６６（図６）の処理の手順を示すフローチャートである。まず、粒子に作用する力を計算する（ステップＳ６６１）。このとき、粒子間の相互作用ポテンシャルに基づいて生じる力については、第２ループＬ２の処理対象の粒子ペアＢＤ、ＢＥ、及びＣＥ（図５）についてのみ計算する。粒子に作用する力が求まると、第２ループＬ２の処理対象となる粒子ペアＢＤ、ＢＥ、及びＣＥ（図５）を構成する粒子Ｂ、Ｃ、Ｄ、Ｅの各々について運動方程式を解き、時間刻み幅Δｔ_２の１／２だけ粒子の速度を時間発展させる（ステップＳ６６２）。ステップＳ６５（図６）とステップＳ６６２との２つのステップで、粒子の速度が時間刻み幅Δｔ_２だけ時間発展することになる。

図８Ｂは、ステップＳ６４（図６）の処理の手順を示すフローチャートである。ステップＳ６６の場合と同様に、第３ループＬ３の処理対象の粒子ペアＤＥ及びＤＦ（図５）について力を計算する（ステップＳ６４１）。その後、第３ループＬ３の処理対象となる粒子ペアＤＥ及びＤＦ（図５）を構成する粒子Ｄ、Ｅ、Ｆの各々について運動方程式を解き、時間刻み幅Δｔ_３の１／２だけ粒子の速度を時間発展させる（ステップＳ６４２）。ステップＳ６３（図６）とステップＳ６４２との２つのステップで、粒子の速度が時間刻み幅Δｔ_３だけ時間発展することになる。

図９は、ステップＳ６２（図６）の処理の手順を示すフローチャートである。ステップＳ６６の場合と同様に、第４ループＬ４の処理対象の粒子ペアＥＦ及びＥＧ（図５）について力を計算する（ステップＳ６２１）。その後、第４ループＬ４の処理対象となる粒子ペアＥＦ及びＥＧ（図５）を構成する粒子Ｅ、Ｆ、Ｇの各々について運動方程式を解き、時間刻み幅Δｔ_４の１／２だけ粒子の速度を時間発展させる（ステップＳ６２２）。ステップＳ６１（図６）とステップＳ６２２との２つのステップで、粒子の速度が時間刻み幅Δｔ_４だけ時間発展することになる。

ステップＳ６２２の後、シミュレーション対象の系の持つエネルギを計算する（ステップＳ６２３）。エネルギには、粒子の運動エネルギ、粒子ペアの弾性エネルギ、物体間の接触エネルギ、及び外部と系との間で入出力された入出力エネルギが含まれる。エネルギ保存の法則により、これらのエネルギの合計はほぼ一定になる。逆に、これらのエネルギの合計がほぼ一定であるならば、シミュレーション結果が正常であると推認される。

ステップＳ６２３の後、シミュレーション結果を出力部２２（図１）に出力する。出力部２２に出力される情報には、例えばシミュレーション対象物の形状の時間変化、代表粒子の位置の時間変化、ステップＳ６２３で求められたエネルギの時間変化等を含めるとよい。

次に、図１０Ａ及び図１０Ｂを参照して物体間の接触の有無を判定する方法（図７のステップＳ６７３）について説明する。

図１０Ａは、シミュレーション対象の２つの物体の形状モデルを示す図である。一方の物体３０の形状モデルが複数の粒子３１の集合体で表され、他方の物体４０の形状モデルが複数の粒子４１の集合体で表されている。物体３０と物体４０とが接触していないときは、一方の物体３０を構成する粒子３１と、他方の物体４０を構成する粒子４１とは、相互作用を及ぼさない。

２つの物体３０、４０が接触したか否かは、例えば、一方の物体３０の表面と他方の物体４０の表面との距離に関する情報に基づいて判定することができる。

図１０Ｂは、物体３０と物体４０とが接触した状態の形状モデルを示す図である。２つの物体３０、４０の表面間の距離が、図１０Ａの場合よりも短い。

次に、２つの物体３０、４０が相互に接触した場合に、一方の物体３０の形状モデルを構成する粒子３１と、他方の物体４０の形状モデルを構成する粒子４１とが及ぼしあう相互作用について説明する。

物体３０と４０とが接触している場合、接触箇所の近傍の粒子３１と４１との間に斥力が作用するように、相互作用ポテンシャルとしてバネマスポテンシャルを適用する。バネ定数は、物体間の食い込み量（計算誤差）をどこまで許容できるかによって決まる。計算誤差を許容できる場合には、バネ定数を小さくすることができる。

各粒子の速度を時間発展させるときに、接触箇所の近傍の粒子３１と４１とについて、相互作用ポテンシャルを考慮する。時間発展させる時間刻み幅Δｔは、同一物体内の粒子ペアと同様に固有周期に基づいて決定する。粒子３１と４１とのペアについて速度の時間発展の計算、及び力の計算を、固有周期に基づいて、第１ループＬ１～第４ループＬ４のいずれのループで行うかが決定される。

次に、上記実施例の優れた効果について説明する。
上記実施例では、シミュレーション対象の形状モデルを構成する複数の粒子に基づく粒子ペアごとに、時間発展の時間刻み幅を固有周期に応じて設定している。粒子間の距離が長い領域では、粒子の分布密度が小さくなるため、粒子の各々に割り振られる質量が大きくなる。このため、固有周期は長くなる。すなわち、粒子ペアの固有周期は、粒子間の距離に依存する。従って、上記実施例において時間刻み幅を固有周期に応じて設定することは、粒子間の距離に応じて時間刻み幅を設定しているともいえる。

このように、１つの形状モデル内で粒子間の距離が相対的に長い粒子ペアについては時間刻み幅を大きくしているため、全体としての計算負荷を軽減することができる。さらに、時間刻み幅を粒子ペアの固有周期に基づいて決定しているため、試行錯誤的に時間刻み幅を決定する場合に比べて、好ましい時間刻み幅を容易に決定することができる。

また、上記実施例では、２つの物体の接触を考慮し、接触箇所に適用する時間発展の時間刻み幅を、物体内の粒子ペアに適用する時間刻み幅と同様の方法で決定している。このため、機構構造系の多体接触に関するシミュレーションに上記実施例を適用する場合にも、計算負荷を軽減させることができる。

次に、図１１Ａ～図１２を参照して、上記実施例の効果を検証するために行ったシミュレーションについて説明する。

図１１Ａは、検証モデルの概略斜視図である。一方向に長い板状部材の長手方向の一方の端部５０を固定し、他方の端部５１に厚さ方向の荷重５２を印加し、板状部材の形状の時間変化をシミュレーションした。すなわち、検証モデルとして、片持ち梁構造の先端に荷重５２を印加するというモデルを採用した。固定された方の端部５０の近傍においてメッシュサイズを相対的に小さくした。メッシュの節点（粒子）の個数は１５，２８３個であり、粒子ペアの個数は９３，７４０個であった。

粒子の状態を時間発展させる３つの異なる時間刻み幅を設定した。時間刻み幅を、短い順にΔｔ_１、Δｔ_２、Δｔ_３と表記する。時間刻み幅Δｔ_２、Δｔ_３を、それぞれΔｔ_１の２倍、４倍とした。粒子ペアの約６０％に時間刻み幅Δｔ_３を適用し、約２５％に時間刻み幅Δｔ_２を適用することが可能であった。残りの約１５％の粒子ペアに時間刻み幅Δｔ_１が適用される。比較のために、全ての粒子ペアに時間刻み幅Δｔ_１を適用した比較例によるシミュレーションも行った。

図１１Ｂは、検証モデルの固定されていない方の端部５１（図１１Ａ）の位置の時間変化のシミュレーション結果を示すグラフである。横軸は時間を単位「ｍｓ」で表し、縦軸は位置を単位「ｍｍ」で表す。図１１Ｂのグラフにおいて、中実の丸記号及び白抜きの丸記号は、それぞれ実施例及び比較例によるシミュレーション結果を示す。両者の丸記号はほぼ重なっており、実施例においても比較例と同等の精度が得られていることが確認された。

図１２は、物体内相互作用の計算、及び位置、速度の計算に要した１タイムステップあたりの計算時間を、実施例と比較例とで対比させて示すグラフである。物体内相互作用の計算は、粒子ペアごとに行う力の計算（図７のステップＳ６７４、図８ＡのステップＳ６６１、図８ＢのステップＳ６４１、図９のステップＳ６２１）に相当する。位置、速度計算は、速度の時間発展、及び位置の時間発展を行う計算に相当する。

一例として図５に示したモデルの場合、粒子Ｅは、第２ループＬ２に対応する粒子ペアＣＥ、ＢＥ、第３ループＬ３に対応する粒子ペアＤＥ、及び第４ループＬ４に対応する粒子ペアＥＦ、ＥＧに含まれている。このため、粒子Ｅの速度の時間発展については、第２ループＬ２、第３ループＬ３、及び第４ループＬ４の全てで実行される。従って、位置、速度の時間発展を行う処理では、実施例の計算時間が比較例の計算時間より長くなっている。

物体内相互作用の計算では、多くの粒子ペアについて時間刻み幅が長くなるため、実施例の計算時間が比較例の計算時間の約４０％まで短くなっている。位置、速度の更新まで含めた全体の計算時間については、実施例の計算時間が比較例の計算時間の約４６％であった。実施例によるシミュレーション方法を適用することにより、比較例と比べて計算時間を短くできることが確認された。

次に、上記実施例の変形例について説明する。
上記実施例では固有周期に基づく時間刻み幅Δｔとして、４種類の時間刻み幅Δｔ_１、Δｔ_２、Δｔ_３、及びΔｔ_４を用いたが、時間刻み幅は２種類以上にすればよい。

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

２０入力部
２１処理部
２２出力部
２３記憶部
３０、４０シミュレーション対象の物体
３１、４１粒子
５０固定された端部
５１荷重が印加される端部
５２荷重

Claims

シミュレーション対象物の形状を定義する形状定義データ、及びシミュレーション対象物をメッシュ分割するときのメッシュサイズに関するメッシュ分割条件が入力される入力部と、
前記入力部に入力された前記形状定義データ及び前記メッシュ分割条件に基づいて、シミュレーション対象物をメッシュ分割し、生成されたメッシュの節点に粒子を配置し、前記粒子の間の相互作用ポテンシャルに基づいて前記粒子の状態を時間発展させる処理部と
を有し、
前記処理部は、
前記粒子の状態を時間発展させる時間刻み幅として、長さの異なる複数の時間刻み幅を定義し、かつ複数の前記時間刻み幅を昇順に並べたとき、前記時間刻み幅のそれぞれが、１つ前の時間刻み幅の整数倍になるように前記時間刻み幅を定義し、
前記シミュレーション対象物の弾性の程度を表す物性値、前記シミュレーション対象物の密度、及び生成されたメッシュの形状に基づいて前記粒子の間の相互作用ポテンシャル及び前記粒子の質量を決定し、
前記粒子の間の相互作用ポテンシャル及び前記粒子の質量から求まる固有振動の周期を求め、
前記粒子のそれぞれについて、当該粒子が含まれる粒子ペアの固有振動の周期に基づいて、複数の前記時間刻み幅から１つの時間刻み幅を決定し、前記粒子の状態を時間発展させる際に、決定された時間刻み幅を用いて時間発展させるシミュレーション装置。
前記シミュレーション対象物が複数の物体を含み、
前記処理部は、
時間発展によって前記物体の間で接触があったか否かを判定し、
前記物体の間で接触があったと判定された場合、接触箇所に配置されている前記粒子の質量及び異なる物体間に作用する相互作用ポテンシャルを導入して、接触箇所に位置する前記粒子の速度及び位置を時間発展させる請求項１に記載のシミュレーション装置。
シミュレーション対象物の形状を定義する形状定義データ、及びシミュレーション対象物をメッシュ分割するときのメッシュサイズに関するメッシュ分割条件を取得する機能と、
前記形状定義データ及び前記メッシュ分割条件に基づいて、シミュレーション対象物をメッシュ分割し、生成されたメッシュの節点に粒子を配置し、前記粒子の間の相互作用ポテンシャルに基づいて前記粒子の状態を時間発展させる機能と、
前記粒子の状態を時間発展させる時間刻み幅として、長さの異なる複数の時間刻み幅を定義し、かつ複数の前記時間刻み幅を昇順に並べたとき、前記時間刻み幅のそれぞれが、１つ前の時間刻み幅の整数倍になるように前記時間刻み幅を定義する機能と、
前記シミュレーション対象物の弾性の程度を表す物性値、前記シミュレーション対象物の密度、及び生成されたメッシュの形状に基づいて前記粒子の間の相互作用ポテンシャル及び前記粒子の質量を決定する機能と、
前記粒子の間の相互作用ポテンシャル及び前記粒子の質量から求まる固有振動の周期を求める機能と、
前記粒子のそれぞれについて、当該粒子が含まれる粒子ペアの固有振動の周期に基づいて、複数の前記時間刻み幅から１つの時間刻み幅を決定し、前記粒子の状態を時間発展させる際に、決定された時間刻み幅を用いて時間発展させる機能と
をコンピュータに実現させるプログラム。