(57)【要約】【課題】解析時間を短縮することが可能なシミュレーション装置を提供する。【解決手段】シミュレーション対象物の形状を定義する形状定義データ、及びシミュレーション対象物をメッシュ分割するときのメッシュサイズに関するメッシュ分割条件が入力部に入力される。処理部が、入力部に入力された形状定義データ及びメッシュ分割条件に基づいて、シミュレーション対象物をメッシュ分割し、生成されたメッシュの節点に粒子を配置し、粒子の間の相互作用ポテンシャルに基づいて粒子の状態を時間発展させる。処理部は、粒子の状態を時間発展させる時間刻み幅を、粒子の間の距離に応じて異ならせる。【選択図】図２
2020095400.tif 000002
【特許請求の範囲】
【請求項１】
シミュレーション対象物の形状を定義する形状定義データ、及びシミュレーション対象物をメッシュ分割するときのメッシュサイズに関するメッシュ分割条件が入力される入力部と、
前記入力部に入力された前記形状定義データ及び前記メッシュ分割条件に基づいて、シミュレーション対象物をメッシュ分割し、生成されたメッシュの節点に粒子を配置し、前記粒子の間の相互作用ポテンシャルに基づいて前記粒子の状態を時間発展させる処理部と
を有し、
前記処理部は、前記粒子の状態を時間発展させる時間刻み幅を、前記粒子の間の距離に応じて異ならせるシミュレーション装置。
【請求項２】
前記処理部は、
前記シミュレーション対象物の弾性の程度を表す物性値、前記シミュレーション対象物の密度、及び生成されたメッシュの形状に基づいて前記粒子の間の相互作用ポテンシャル及び前記粒子の質量を決定し、
前記粒子の間の相互作用ポテンシャル及び前記粒子の質量から求まる固有振動の周期に基づいて粒子ペアごとに前記時間刻み幅を決定する請求項１に記載のシミュレーション装置。
【請求項３】
前記シミュレーション対象物が複数の物体を含み、
前記処理部は、
時間発展によって前記物体の間で接触があったか否かを判定し、
前記物体の間で接触があったと判定された場合、接触箇所に配置されている前記粒子の質量及び相互作用ポテンシャルから求まる固有振動の周期に基づいて、接触箇所に位置する前記粒子について前記時間刻み幅を決定する請求項１または２に記載のシミュレーション装置。
【請求項４】
シミュレーション対象物の形状をメッシュ分割し、生成されたメッシュの節点に粒子を配置し、
前記粒子の状態を時間発展させる時間刻み幅を、前記粒子の間の距離に応じて異ならせて、前記粒子の間の相互作用ポテンシャルに基づいて前記粒子の状態を時間発展させるシミュレーション方法。
【請求項５】
シミュレーション対象物の形状を定義する形状定義データ、及びシミュレーション対象物をメッシュ分割するときのメッシュサイズに関するメッシュ分割条件を取得する機能と、
前記形状定義データ及び前記メッシュ分割条件に基づいて、シミュレーション対象物をメッシュ分割し、生成されたメッシュの節点に粒子を配置し、前記粒子の間の相互作用ポテンシャルに基づいて前記粒子の状態を時間発展させる機能と、
前記粒子の状態を時間発展させる時間刻み幅を、前記粒子の間の距離に応じて異ならせる機能と
コンピュータに実行させるプログラム。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、シミュレーション装置、シミュレーション方法、及びプログラムに関する。
【背景技術】
【０００２】
分子動力学法、繰り込み群分子動力学法等の動的陽解法を用いる機構構造解析では、計算を発散させないためにメッシュサイズ（粒子間距離）に応じた時間刻み幅以下で時間発展を行うことが好ましい。メッシュサイズが小さいほど（粒子間距離が短いほど）、好ましい時間刻み幅は小さくなる。解析対象の系のメッシュサイズが場所によって異なっている場合、系内の最も小さなメッシュサイズに合わせて時間刻み幅を小さくしなければならない。このため、演算量が増え、解析時間が長くなってしまう。
【０００３】
下記の特許文献１に、サルコメア内の分子の運動と、筋肉の連続体モデルを連成させて解析を行う生体シミュレーション装置が開示されている。サルコメアモデルのモンテカルロシミュレーション処理の時間刻み幅と、連続体モデルの有限要素解析の時間刻み幅とは異なっている。
【先行技術文献】
【特許文献】
【０００４】
【特許文献１】特開２０１４−１４９６４９号公報
【発明の概要】
【発明が解決しようとする課題】
【０００５】
従来例において、動的陽解法を用いた連続体モデルの有限要素解析の時間刻み幅は連続体モデルのメッシュサイズに依らず一定である。このため、最も小さなメッシュサイズに対応して時間刻み幅を短くしなければならず、解析時間を短縮することが困難である。
【０００６】
本発明の目的は、解析時間を短縮することが可能なシミュレーション装置、シミュレーション方法、及びプログラムを提供することである。
【課題を解決するための手段】
【０００７】
本発明の一観点によると、
シミュレーション対象物の形状を定義する形状定義データ、及びシミュレーション対象物をメッシュ分割するときのメッシュサイズに関するメッシュ分割条件が入力される入力部と、
前記入力部に入力された前記形状定義データ及び前記メッシュ分割条件に基づいて、シミュレーション対象物をメッシュ分割し、生成されたメッシュの節点に粒子を配置し、前記粒子の間の相互作用ポテンシャルに基づいて前記粒子の状態を時間発展させる処理部と
を有し、
前記処理部は、前記粒子の状態を時間発展させる時間刻み幅を、前記粒子の間の距離に応じて異ならせるシミュレーション装置が提供される。
【０００８】
本発明の他の観点によると、
シミュレーション対象物の形状をメッシュ分割し、生成されたメッシュの節点に粒子を配置し、
前記粒子の状態を時間発展させる時間刻み幅を、前記粒子の間の距離に応じて異ならせて、前記粒子の間の相互作用ポテンシャルに基づいて前記粒子の状態を時間発展させるシミュレーション方法が提供される。
【０００９】
本発明のさらに他の観点によると、
シミュレーション対象物の形状を定義する形状定義データ、及びシミュレーション対象物をメッシュ分割するときのメッシュサイズに関するメッシュ分割条件を取得する機能と、
前記形状定義データ及び前記メッシュ分割条件に基づいて、シミュレーション対象物をメッシュ分割し、生成されたメッシュの節点に粒子を配置し、前記粒子の間の相互作用ポテンシャルに基づいて前記粒子の状態を時間発展させる機能と、
前記粒子の状態を時間発展させる時間刻み幅を、前記粒子の間の距離に応じて異ならせる機能と
コンピュータに実行させるプログラムが提供される。
【発明の効果】
【００１０】
粒子の状態を時間発展させる時間刻み幅を、前記粒子の間の距離に応じて異ならせることにより、解析時間を短縮することが可能になる。
【図面の簡単な説明】
【００１１】
【図１】図１は、実施例によるシミュレーション装置のブロック図である。
【図２】図２は、本実施例によるシミュレーション方法のフローチャートである。
【図３】図３は、固有周期Ｐと、時間刻み幅Δｔとの関係の一例を示すグラフである。
【図４】図４は、第１ループＬ１〜第４ループＬ４を実行するタイミングチャートである。
【図５】図５は、シミュレーション対象物の形状から生成されたメッシュの複数の節点に配置した複数の粒子の一例を示す図である。
【図６】図６は、ステップＳ６（図２）の詳細な手順を示すフローチャートである。
【図７】図７は、第１ループＬ１の処理（図６のステップＳ６７）の手順を示すフローチャートである。
【図８】図８Ａは、ステップＳ６６（図６）の処理の手順を示すフローチャートであり、図８Ｂは、ステップＳ６４（図６）の処理の手順を示すフローチャートである。
【図９】図９は、ステップＳ６２（図６）の処理の手順を示すフローチャートである。
【図１０】図１０Ａは、シミュレーション対象の２つの物体の形状モデルを示す図であり、図１０Ｂは、２つの物体が接触した状態の形状モデルを示す図である。
【図１１】図１１Ａは、検証モデルの概略斜視図であり、図１１Ｂは、検証モデルの固定されていない方の端部の位置の時間変化のシミュレーション結果を示すグラフである。
【図１２】図１２は、物体内相互作用計算、及び位置、速度計算に要した１タイムステップあたりの計算時間を、実施例と比較例とで対比させて示すグラフである。
【発明を実施するための形態】
【００１２】
図１〜図１２を参照して実施例によるシミュレーション装置、シミュレーション方法について説明する。
【００１３】
図１は、実施例によるシミュレーション装置のブロック図である。実施例によるシミュレーション装置は、入力部２０、処理部２１、出力部２２、及び記憶部２３を含む。入力部２０から処理部２１にシミュレーション条件等が入力される。さらに、オペレータから入力部２０に各種指令（コマンド）等が入力される。入力部２０は、例えば通信装置、リムーバブルメディア読取装置、キーボード等で構成される。
【００１４】
処理部２１は、入力されたシミュレーション条件及び指令に基づいてシミュレーション対象物の形状をメッシュ分割し、メッシュの節点に仮想的に粒子を配置して分子動力学法によるシミュレーションを行う。さらに、シミュレーション結果を出力部２２に出力する。シミュレーション結果には、シミュレーション対象物の形状の時間的変化を表す情報が含まれる。処理部２１は、例えばコンピュータを含み、くりこみ群分子動力学法によるシミュレーションをコンピュータに実行させるためのプログラムが記憶部２３に記憶されている。出力部２２は、通信装置、リムーバブルメディア書込み装置、ディスプレイ等を含む。
【００１５】
図２は、本実施例によるシミュレーション方法のフローチャートである。図２に示した各ステップの処理は、処理部２１（図１）が記憶部２３に記憶されたプログラムを実行することにより実現される。
【００１６】
まず、処理部２１は、シミュレーション対象物の形状を定義する形状定義データ、シミュレーション対象物の物性値（密度、弾性の程度を表す物性値等）、シミュレーションの初期条件、境界条件、メッシュ分割条件等を取得する（ステップＳ１）。これらのデータは入力部２０に入力され、処理部２１が入力部２０から取得する。形状定義データを取得すると、処理部２１は形状定義データで定義される形状をメッシュ分割条件に基づいてメッシュ分割することにより、形状モデルを生成する（ステップＳ２）。メッシュ分割処理には、例えば公知のメッシュ分割アルゴリズムを使用することができる。本実施例では、シミュレーション対象物をテトラメッシュに分割する。
【００１７】
メッシュが生成されると、処理部２１は、メッシュの節点に仮想的に粒子を配置し、各粒子に質量を付与し、粒子間の相互作用ポテンシャルを決定する（ステップＳ３）。粒子の質量は、例えば、シミュレーション対象物の密度と、粒子の分布とに基づいて、形状モデルの密度がシミュレーション対象物の密度と同一になるように決定する。例えば、シミュレーション対象物の密度が均一である場合、粒子の分布密度が高い領域では粒子の質量が相対的に小さくなり、粒子の分布密度が低い領域では粒子の質量が相対的に大きくなる。
【００１８】
粒子間の相互作用ポテンシャルとして、例えばバネマスモデルのポテンシャルを用いる。バネ定数の決定方法は、例えば特開２００９−３７３３４号公報に詳しく説明されている。ここでは、バネ定数の決定方法について簡単に説明する。
【００１９】
まず、テトラメッシュの節点に配置された粒子のボロノイ多面体解析を行う。ボロノイ多面体を構成する複数の界面のうち粒子ｉと粒子ｊとを両端とする線分を横切る界面の面積をＳｉｊで表す。シミュレーション対象物のヤング率と面積Ｓｉｊとからバネ定数ｋｉｊを決定することができる。
【００２０】
次に、バネで相互に接続されている２つの粒子（粒子ペア）ごとに、固有振動の周期（固有周期）を算出する（ステップＳ４）。一般的に、粒子ペアを構成する２つの粒子の質量は同一ではないため、一方の粒子を固定した状態における固有周期と、他方の粒子を固定した状態における固有周期とは異なる。粒子のペアの固有周期として、例えば、２つの固有周期のうち短い方の固有周期の値を採用する。
【００２１】
処理部２１は、ステップＳ４で算出された固有周期に基づいて、粒子ペアごとに、粒子の状態の時間発展の演算を行う際の時間刻み幅を決定する（ステップＳ５）。固有周期が短いほど、一定時間に粒子の状態が大きく変化する。粒子の状態の変化を精度よく解析するために、固有周期が短いほど時間刻み幅を短くするとよい。例えば、時間刻み幅は、固有周期の１／２０以上１／１０以下の範囲から選択するとよい。１つの形状モデルに対して、複数の時間刻み幅が定義される。
【００２２】
時間刻み幅が決定されると、処理部２１は、決定された時間刻み幅で粒子の状態の変化を解析する（ステップＳ６）。具体的には、運動方程式を解くことにより粒子の速度と位置の時間変化を求める。解析が終了すると、解析結果を出力部２２（図１）に出力する（ステップＳ７）。
【００２３】
図３は、固有周期Ｐと、時間刻み幅Δｔとの関係の一例を示すグラフである。横軸は固有周期Ｐを表し、縦軸は時間刻み幅Δｔを表す。例えば、固有周期ＰがＰ１未満、Ｐ１以上Ｐ２未満、Ｐ２以上Ｐ３未満、及びＰ３以上のとき、それぞれ時間刻み幅Δｔを、Δｔ１、Δｔ２、Δｔ３、及びΔｔ４とする。ここで、Δｔ１＜Δｔ２＜Δｔ３＜Δｔ４の大小関係が成立する。さらに、Δｔ４はΔｔ３の整数倍であり、Δｔ３はΔｔ２の整数倍であり、Δｔ２はΔｔ１の整数倍である。
【００２４】
固有周期ＰがＰ３以上の粒子ペアについては、時間刻み幅Δｔ４で第４ループＬ４の処理を実行することにより、粒子間に働く力の計算、及び速度の時間発展を行う。固有周期ＰがＰ２以上Ｐ３未満の粒子ペアについては、時間刻み幅Δｔ３で第３ループＬ３の処理を実行する。固有周期ＰがＰ１以上Ｐ２未満の粒子ペアについては、時間刻み幅Δｔ２で第２ループＬ２の処理を実行する。固有周期ＰがＰ１未満の粒子ペアについては、時間刻み幅Δｔ１で第１ループＬ１の処理を実行する。
【００２５】
図４は、第１ループＬ１〜第４ループＬ４を実行するタイミングチャートである。第１ループＬ１の処理は、時間刻み幅Δｔ１ごとに実行される。言い換えると、第１ループＬ１の処理を１回実行すると、時間刻み幅Δｔ１だけ時間発展する。この時間発展の累積が時間刻み幅Δｔ２だけ増加するごとに、第２ループＬ２の処理が実行される。時間発展の累積が時間刻み幅Δｔ３だけ増加するごとに、第３ループＬ３の処理が実行される。時間発展の累積が時間刻み幅Δｔ４だけ増加するごとに、第４ループＬ４の処理が実行される。時間刻み幅Δｔ１だけ時間発展させる１回の処理を１つのタイムステップということとする。
【００２６】
図５は、シミュレーション対象物の形状から生成されたメッシュの複数の節点に配置した複数の粒子の一例を示す図である。図５に示した複数の粒子から複数の粒子ペアが定義され、複数の粒子ペアにそれぞれ固有周期に基づいて第１ループＬ１〜第４ループＬ４のいずれかが対応付けられる。例えば、粒子ペアＡＢ、ＢＣに第１ループＬ１が対応付けられている。粒子ペアＢＤ、ＢＥ、ＣＥに第２ループＬ２が対応付けられている。粒子ペアＤＥ、ＤＦに第３ループＬ３が対応付けられている。粒子ペアＥＦ、ＥＧに第４ループＬ４が対応付けられている。
【００２７】
第１ループＬ１（図４）のみを実行するタイムステップでは、粒子ペアＡＢ及びＢＣについて、粒子の速度の時間発展、及び力の計算を行う。第１ループＬ１及び第２ループＬ２（図４）を実行するタイムステップでは、さらに、粒子ペアＢＤ、ＢＥ、及びＣＥについても、粒子の速度の時間発展、及び力の計算を行う。第１ループＬ１、第２ループＬ２、及び第３ループＬ３（図４）を実行するタイムステップでは、さらに、粒子ペアＤＥ、及びＤＦについても、粒子の速度の時間発展、及び力の計算を行う。第１ループＬ１〜第４ループＬ４の全てを実行するタイムステップでは、さらに粒子ペアＥＦ、及びＥＧについても、粒子の速度の時間発展、及び力の計算を行う。
【００２８】
図６は、粒子の状態の変化を解析する処理（図２のステップＳ６）の詳細な手順を示すフローチャートである。以下の説明において、適宜、図５を参照する。
【００２９】
シミュレーションの演算終了条件が満たされるまで、第４ループＬ４の処理を実行する。演算終了条件として、例えばタイムステップ数、演算時間等が与えられる。演算終了条件が満たされない場合には、まず、第４ループＬ４の処理対象となる粒子ペアＥＦ、及びＥＧ（図５）に含まれる粒子Ｅ、Ｆ、Ｇについて運動方程式を解き、時間刻み幅Δｔ４の１／２だけ速度を時間発展させる（ステップＳ６１）。
【００３０】
ステップＳ６１を実行した後、第３ループＬ３の処理対象となる粒子ペアＤＥ、及びＤＦ（図５）に含まれる粒子Ｄ、Ｅ、Ｆについて運動方程式を解き、時間刻み幅Δｔ３の１／２だけ各粒子の速度を時間発展させる（ステップＳ６３）。ここで、ステップＳ６１で速度を時間発展させた粒子、例えば粒子Ｅ、Ｆについては、時間発展後の最新の速度を、さらに時間発展させる。以降の速度の時間発展の処理においても同様に、最新の速度をさらに時間発展させる。
【００３１】
ステップＳ６３を実行した後、第２ループＬ２の処理対象となる粒子ペアＢＤ、ＢＥ、及びＣＥ（図５）に含まれる粒子Ｂ、Ｃ、Ｄ、Ｅについて運動方程式を解き、時間刻み幅Δｔ２の１／２だけ各粒子の速度を時間発展させる（ステップＳ６３）。
【００３２】
ステップＳ６５を実行した後、全ての粒子について第１ループＬ１の処理を実行する（ステップＳ６７）。第１ループＬ１の処理の詳細については、後に図７を参照して説明する。
【００３３】
第１ループＬ１の処理（タイムステップ）を複数回繰り返して、ステップＳ６５以降の時間発展の累積値がΔｔ２だけ増加すると、第２ループＬ２の後処理を実行する（ステップＳ６６）。第２ループＬ２の後処理（ステップＳ６６）では、時間刻み幅Δｔ２の１／２だけ速度をさらに時間発展させる。第２ループＬ２の後処理（ステップＳ６６）については、後に図８Ａを参照して詳細に説明する。このように、速度の時間発展を時間刻み幅の１／２ずつ実行する方法は、速度ベルレ法と呼ばれる。
【００３４】
第２ループＬ２の処理を複数回繰り返して、ステップＳ６３以降の時間発展の累積値がΔｔ３だけ増加すると、第３ループＬ３の後処理を実行する（ステップＳ６４）。第３ループＬ３の後処理（ステップＳ６４）では、時間刻み幅Δｔ３の１／２だけ速度をさらに時間発展させる。第３ループＬ３の後処理（ステップＳ６４）については、後に図８Ｂを参照して詳細に説明する。
【００３５】
第３ループＬ３の処理を複数回繰り返した後、ステップＳ６１以降の時間発展の累積値がΔｔ４だけ増加すると、第４ループＬ４の後処理を実行する（ステップＳ６２）。第４ループＬ４の後処理（ステップＳ６２）では、時間刻み幅Δｔ４の１／２だけ速度をさらに時間発展させる。第４ループＬ４の後処理（ステップＳ６２）については、後に図９を参照して詳細に説明する。
【００３６】
図７は、第１ループＬ１の処理（図６のステップＳ６７）の手順を示すフローチャートである。まず、第１ループＬ１の処理対象となる粒子ペアＡＢ、及びＢＣ（図５）に含まれる粒子Ａ、Ｂ、Ｃについて運動方程式を解き、時間刻み幅Δｔ１の１／２だけ速度を時間発展させる（ステップＳ６７１）。次に、全ての粒子の位置を、各粒子の最新の速度に基づいて時間刻み幅Δｔ１だけ時間発展させる（ステップＳ６７２）。
【００３７】
全ての粒子の位置を時間発展させた後、物体間の接触の有無を判定する（ステップＳ６７３）。例えば、シミュレーション対象物が複数の物体で構成されている場合、同一の物体内の粒子の間には、ステップＳ３（図２）で決定した相互作用ポテンシャルに基づく力が作用する。異なる物体の粒子の間には、物体が接触していない状態では力は作用しない。ところが、２つの物体が接触すると、物体の接触箇所に位置する粒子の間で力を及ぼし合う。ステップＳ６７３で物体間の接触があったと判定された場合には、次のタイムステップにおいて、物体の接触箇所に位置する粒子に対して、異なる物体の粒子間に作用する相互作用ポテンシャルを導入して速度及び位置を時間発展させる。物体間の接触の判定、及び異なる物体の粒子間の相互作用ポテンシャルについては、後に図１１Ａ及び図１１Ｂを参照して説明する。
【００３８】
次に、粒子の時間発展後の位置に基づいて、粒子に作用する力を計算する（ステップＳ６７４）。このとき、粒子間の相互作用ポテンシャルに基づいて生じる力については、第１ループＬ１の処理対象の粒子ペアＡＢ及びＢＣ（図５）についてのみ計算する。また、重力や境界条件として与えられた荷重等の外力が作用する粒子については、これらの外力を計算する。第２ループＬ２〜第４ループＬ４の処理対象の粒子ペアについては、相互作用ポテンシャルに基づく力の計算は行わない。各粒子について運動方程式を解く際には、当該粒子に作用する力を合成した合成力に基づいて計算を行う。
【００３９】
粒子に作用する力が求まると、第１ループＬ１の処理対象となる粒子ペアＡＢ及びＢＣを構成する粒子Ａ、Ｂ、Ｃの各々について運動方程式を解き、時間刻み幅Δｔ１の１／２だけ粒子の速度を時間発展させる（ステップＳ６７５）。ステップＳ６７１とステップＳ６７５とで、それぞれ粒子の速度が時間刻み幅Δｔ１の１／２ずつ時間発展するため、粒子の速度は時間刻み幅Δｔ１だけ時間発展することになる。
【００４０】
図８Ａは、ステップＳ６６（図６）の処理の手順を示すフローチャートである。まず、粒子に作用する力を計算する（ステップＳ６６１）。このとき、粒子間の相互作用ポテンシャルに基づいて生じる力については、第２ループＬ２の処理対象の粒子ペアＢＤ、ＢＥ、及びＣＥ（図５）についてのみ計算する。粒子に作用する力が求まると、第２ループＬ２の処理対象となる粒子ペアＢＤ、ＢＥ、及びＣＥ（図５）を構成する粒子Ｂ、Ｃ、Ｄ、Ｅの各々について運動方程式を解き、時間刻み幅Δｔ２の１／２だけ粒子の速度を時間発展させる（ステップＳ６６２）。ステップＳ６５（図６）とステップＳ６６２との２つのステップで、粒子の速度が時間刻み幅Δｔ２だけ時間発展することになる。
【００４１】
図８Ｂは、ステップＳ６４（図６）の処理の手順を示すフローチャートである。ステップＳ６６の場合と同様に、第３ループＬ３の処理対象の粒子ペアＤＥ及びＤＦ（図５）について力を計算する（ステップＳ６４１）。その後、第３ループＬ３の処理対象となる粒子ペアＤＥ及びＤＦ（図５）を構成する粒子Ｄ、Ｅ、Ｆの各々について運動方程式を解き、時間刻み幅Δｔ３の１／２だけ粒子の速度を時間発展させる（ステップＳ６４２）。ステップＳ６３（図６）とステップＳ６４２との２つのステップで、粒子の速度が時間刻み幅Δｔ３だけ時間発展することになる。
【００４２】
図９は、ステップＳ６２（図６）の処理の手順を示すフローチャートである。ステップＳ６６の場合と同様に、第４ループＬ４の処理対象の粒子ペアＥＦ及びＥＧ（図５）について力を計算する（ステップＳ６２１）。その後、第４ループＬ４の処理対象となる粒子ペアＥＦ及びＥＧ（図５）を構成する粒子Ｅ、Ｆ、Ｇの各々について運動方程式を解き、時間刻み幅Δｔ４の１／２だけ粒子の速度を時間発展させる（ステップＳ６２２）。ステップＳ６１（図６）とステップＳ６２２との２つのステップで、粒子の速度が時間刻み幅Δｔ４だけ時間発展することになる。
【００４３】
ステップＳ６２２の後、シミュレーション対象の系の持つエネルギを計算する（ステップＳ６２３）。エネルギには、粒子の運動エネルギ、粒子ペアの弾性エネルギ、物体間の接触エネルギ、及び外部と系との間で入出力された入出力エネルギが含まれる。エネルギ保存の法則により、これらのエネルギの合計はほぼ一定になる。逆に、これらのエネルギの合計がほぼ一定であるならば、シミュレーション結果が正常であると推認される。
【００４４】
ステップＳ６２３の後、シミュレーション結果を出力部２２（図１）に出力する。出力部２２に出力される情報には、例えばシミュレーション対象物の形状の時間変化、代表粒子の位置の時間変化、ステップＳ６２３で求められたエネルギの時間変化等を含めるとよい。
【００４５】
次に、図１０Ａ及び図１０Ｂを参照して物体間の接触の有無を判定する方法（図７のステップＳ６７３）について説明する。
【００４６】
図１０Ａは、シミュレーション対象の２つの物体の形状モデルを示す図である。一方の物体３０の形状モデルが複数の粒子３１の集合体で表され、他方の物体４０の形状モデルが複数の粒子４１の集合体で表されている。物体３０と物体４０とが接触していないときは、一方の物体３０を構成する粒子３１と、他方の物体４０を構成する粒子４１とは、相互作用を及ぼさない。
【００４７】
２つの物体３０、４０が接触したか否かは、例えば、一方の物体３０の表面と他方の物体４０の表面との距離に関する情報に基づいて判定することができる。
【００４８】
図１０Ｂは、物体３０と物体４０とが接触した状態の形状モデルを示す図である。２つの物体３０、４０の表面間の距離が、図１０Ａの場合よりも短い。
【００４９】
次に、２つの物体３０、４０が相互に接触した場合に、一方の物体３０の形状モデルを構成する粒子３１と、他方の物体４０の形状モデルを構成する粒子４１とが及ぼしあう相互作用について説明する。
【００５０】
物体３０と４０とが接触している場合、接触箇所の近傍の粒子３１と４１との間に斥力が作用するように、相互作用ポテンシャルとしてバネマスポテンシャルを適用する。バネ定数は、物体間の食い込み量（計算誤差）をどこまで許容できるかによって決まる。計算誤差を許容できる場合には、バネ定数を小さくすることができる。
【００５１】
各粒子の速度を時間発展させるときに、接触箇所の近傍の粒子３１と４１とについて、相互作用ポテンシャルを考慮する。時間発展させる時間刻み幅Δｔは、同一物体内の粒子ペアと同様に固有周期に基づいて決定する。粒子３１と４１とのペアについて速度の時間発展の計算、及び力の計算を、固有周期に基づいて、第１ループＬ１〜第４ループＬ４のいずれのループで行うかが決定される。
【００５２】
次に、上記実施例の優れた効果について説明する。
上記実施例では、シミュレーション対象の形状モデルを構成する複数の粒子に基づく粒子ペアごとに、時間発展の時間刻み幅を固有周期に応じて設定している。粒子間の距離が長い領域では、粒子の分布密度が小さくなるため、粒子の各々に割り振られる質量が大きくなる。このため、固有周期は長くなる。すなわち、粒子ペアの固有周期は、粒子間の距離に依存する。従って、上記実施例において時間刻み幅を固有周期に応じて設定することは、粒子間の距離に応じて時間刻み幅を設定しているともいえる。
【００５３】
このように、１つの形状モデル内で粒子間の距離が相対的に長い粒子ペアについては時間刻み幅を大きくしているため、全体としての計算負荷を軽減することができる。さらに、時間刻み幅を粒子ペアの固有周期に基づいて決定しているため、試行錯誤的に時間刻み幅を決定する場合に比べて、好ましい時間刻み幅を容易に決定することができる。
【００５４】
また、上記実施例では、２つの物体の接触を考慮し、接触箇所に適用する時間発展の時間刻み幅を、物体内の粒子ペアに適用する時間刻み幅と同様の方法で決定している。このため、機構構造系の多体接触に関するシミュレーションに上記実施例を適用する場合にも、計算負荷を軽減させることができる。
【００５５】
次に、図１１Ａ〜図１２を参照して、上記実施例の効果を検証するために行ったシミュレーションについて説明する。
【００５６】
図１１Ａは、検証モデルの概略斜視図である。一方向に長い板状部材の長手方向の一方の端部５０を固定し、他方の端部５１に厚さ方向の荷重５２を印加し、板状部材の形状の時間変化をシミュレーションした。すなわち、検証モデルとして、片持ち梁構造の先端に荷重５２を印加するというモデルを採用した。固定された方の端部５０の近傍においてメッシュサイズを相対的に小さくした。メッシュの節点（粒子）の個数は１５，２８３個であり、粒子ペアの個数は９３，７４０個であった。
【００５７】
粒子の状態を時間発展させる３つの異なる時間刻み幅を設定した。時間刻み幅を、短い順にΔｔ１、Δｔ２、Δｔ３と表記する。時間刻み幅Δｔ２、Δｔ３を、それぞれΔｔ１の２倍、４倍とした。粒子ペアの約６０％に時間刻み幅Δｔ３を適用し、約２５％に時間刻み幅Δｔ２を適用することが可能であった。残りの約１５％の粒子ペアに時間刻み幅Δｔ１が適用される。比較のために、全ての粒子ペアに時間刻み幅Δｔ１を適用した比較例によるシミュレーションも行った。
【００５８】
図１１Ｂは、検証モデルの固定されていない方の端部５１（図１１Ａ）の位置の時間変化のシミュレーション結果を示すグラフである。横軸は時間を単位「ｍｓ」で表し、縦軸は位置を単位「ｍｍ」で表す。図１１Ｂのグラフにおいて、中実の丸記号及び白抜きの丸記号は、それぞれ実施例及び比較例によるシミュレーション結果を示す。両者の丸記号はほぼ重なっており、実施例においても比較例と同等の精度が得られていることが確認された。
【００５９】
図１２は、物体内相互作用の計算、及び位置、速度の計算に要した１タイムステップあたりの計算時間を、実施例と比較例とで対比させて示すグラフである。物体内相互作用の計算は、粒子ペアごとに行う力の計算（図７のステップＳ６７４、図８ＡのステップＳ６６１、図８ＢのステップＳ６４１、図９のステップＳ６２１）に相当する。位置、速度計算は、速度の時間発展、及び位置の時間発展を行う計算に相当する。
【００６０】
一例として図５に示したモデルの場合、粒子Ｅは、第２ループＬ２に対応する粒子ペアＣＥ、ＢＥ、第３ループＬ３に対応する粒子ペアＤＥ、及び第４ループＬ４に対応する粒子ペアＥＦ、ＥＧに含まれている。このため、粒子Ｅの速度の時間発展については、第２ループＬ２、第３ループＬ３、及び第４ループＬ４の全てで実行される。従って、位置、速度の時間発展を行う処理では、実施例の計算時間が比較例の計算時間より長くなっている。
【００６１】
物体内相互作用の計算では、多くの粒子ペアについて時間刻み幅が長くなるため、実施例の計算時間が比較例の計算時間の約４０％まで短くなっている。位置、速度の更新まで含めた全体の計算時間については、実施例の計算時間が比較例の計算時間の約４６％であった。実施例によるシミュレーション方法を適用することにより、比較例と比べて計算時間を短くできることが確認された。
【００６２】
次に、上記実施例の変形例について説明する。
上記実施例では固有周期に基づく時間刻み幅Δｔとして、４種類の時間刻み幅Δｔ１、Δｔ２、Δｔ３、及びΔｔ４を用いたが、時間刻み幅は２種類以上にすればよい。
【００６３】
上述の実施例は例示であり、本発明は上述の実施例に制限されるものではない。例えば、種々の変更、改良、組み合わせ等が可能なことは当業者に自明であろう。
【符号の説明】
【００６４】
２０  入力部
２１  処理部
２２  出力部
２３  記憶部
３０、４０  シミュレーション対象の物体
３１、４１  粒子
５０  固定された端部
５１  荷重が印加される端部
５２  荷重
【図１】
2020095400.tif 000003
【図２】
2020095400.tif 000004
【図３】
2020095400.tif 000005
【図４】
2020095400.tif 000006
【図５】
2020095400.tif 000007
【図６】
2020095400.tif 000008
【図７】
2020095400.tif 000009
【図８】
2020095400.tif 000010
【図９】
2020095400.tif 000011
【図１０】
2020095400.tif 000012
【図１１】
2020095400.tif 000013
【図１２】
2020095400.tif 000014
【手続補正書】
【提出日】令和1年9月13日(2019.9.13)
【手続補正１】
【補正対象書類名】明細書
【補正対象項目名】００２９
【補正方法】変更
【補正の内容】
【００２９】
シミュレーションの演算終了条件が満たされるまで、第４ループＬ４の処理を実行する。演算終了条件として、例えばタイムステップ数、演算時間等が与えられる。演算終了条件が満たされない場合には、まず、第４ループＬ４の処理対象となる粒子ペアＥＦ、及びＥＧ（図５）に含まれる粒子Ｅ、Ｆ、Ｇについて運動方程式を解き、時間刻み幅Δｔ４の１／２だけ各粒子の速度を時間発展させる（ステップＳ６１）。
【手続補正２】
【補正対象書類名】明細書
【補正対象項目名】００３１
【補正方法】変更
【補正の内容】
【００３１】
ステップＳ６３を実行した後、第２ループＬ２の処理対象となる粒子ペアＢＤ、ＢＥ、及びＣＥ（図５）に含まれる粒子Ｂ、Ｃ、Ｄ、Ｅについて運動方程式を解き、時間刻み幅Δｔ２の１／２だけ各粒子の速度を時間発展させる（ステップＳ６５）。
【手続補正３】
【補正対象書類名】図面
【補正対象項目名】図６
【補正方法】変更
【補正の内容】
【図６】
2020095400.tif 000015