JP2020173185A

JP2020173185A - 構造物のｆｅｍ解析方法、システム及びプログラム

Info

Publication number: JP2020173185A
Application number: JP2019075608A
Authority: JP
Inventors: 理日野; Tadashi Hino
Original assignee: Toyo Tire and Rubber Co Ltd; Toyo Tire Corp
Current assignee: Toyo Tire Corp
Priority date: 2019-04-11
Filing date: 2019-04-11
Publication date: 2020-10-22
Anticipated expiration: 2039-04-11
Also published as: JP7248487B2

Abstract

【課題】リーズナブルに並列計算を実行可能な構造物のＦＥＭ解析方法を提供する。【解決手段】構造物を複数の要素ｅで表現したＦＥＭモデルデータＤ１に基づき構造物の質量を各節点に配分して表す集中質量行列Ｍｅと各節点の初期位置及び初期速度とを取得し、各節点の初期位置及び物性値に基づき要素剛性行列Ｋｅを要素毎に算出し、（ｄ）各節点の初期位置に対する時点（ｔ）の位置の変位を表す変位行列ｄｅを要素毎に算出し、（ｅ）要素剛性行列Ｋｅと変位行列ｄｅとに基づき、各節点に作用する力Ｆｅを要素毎に算出し、（ｆ）時点（ｔ）における前記各節点の位置と速度と前記各節点に作用する力と前記集中質量行列とに基づき分子動力学ソルバに次の時点における各節点の位置及び速度を算出させる。初期時点から目標の時点に到達するまで、（ｄ）（ｅ）及び（ｆ）の処理を、時点（ｔ）を単位時間経過させつつ繰り返し実行する。【選択図】図２

Description

本発明は、構造物のＦＥＭ解析方法、システム及びプログラムに関する。

構造物の動解析には有限要素法（ＦＥＭ；Finite Element Method）が用いられる。ＦＥＭに用いるＦＥＭモデルデータは、解析対象となる構造物が複数の要素に分割されて表現され、複数の節点を結んだ閉領域を一つの要素として表現されている。ＦＥＭでは、節点毎に、構造物の物性値（ヤング率、ポアソン比、質量など）及び作用する力を用いた運動方程式を解くことにより、構造物の変形などを計算する。ＦＥＭを説明する文献として非特許文献１がある。

なお、構造物シミュレーションではないが、分子シミュレーションを実行するための分子動力学ソルバ（ソフトウェア）は、多粒子系の運動方程式を高効率で解くための並列計算機能を有し、ＬＡＭＭＰＳやｇｒｏｍａｃｓ等の無償利用可能なソフトウェアが知られている。分子動力学に関する文献として非特許文献２がある。

Hutton, D. V. (2017). Fundamentals of finite element analysis, Chapter 10, Structural Dynamics. McGraw-Hill. Frenkel, D., Smit, B., Tobochnik, J., McKay, S. R., & Christian, W. (1997). Understanding molecular simulation. Computers in Physics, 11(4), 351-354.

ＦＥＭ解析において解析精度を向上させるためには、要素のサイズを小さくしたモデルを用いればよいが、要素サイズを小さくすれば、要素数が増大し、計算時間が急速に増加してしまう。ＦＥＭ解析ソフトウェアの大半は、全ての要素を一体に取り扱う行列計算を行っているせいか、並列計算機能がないことが多く、並列計算機能を有するＦＥＭ解析ソフトウェアはライセンスコストが高価で実用的ではない。

本発明の目的は、リーズナブルに並列計算を実行可能な構造物のＦＥＭ解析方法、システム及びプログラムを提供することである。

本発明の構造物のＦＥＭ解析方法は、
１又は複数のプロセッサが実行する方法であって、
（ａ）構造物を複数の要素で表現したＦＥＭモデルデータであって、各要素を区画する複数の節点および物性値が設定されたＦＥＭモデルデータを取得すること、
（ｂ）前記ＦＥＭモデルデータに基づき、前記構造物の質量を各節点に配分して表す集中質量行列と、各節点の初期位置及び初期速度と、を取得すること、
（ｃ）前記各節点の初期位置及び前記物性値に基づき要素剛性行列を要素毎に算出すること、
（ｄ）各節点の初期位置に対する時点（ｔ）の位置の変位を表す変位行列を要素毎に算出すること、
（ｅ）前記要素剛性行列と前記変位行列とに基づき、各節点に作用する力を要素毎に算出すること、
（ｆ）時点（ｔ）における前記各節点の位置と速度と前記各節点に作用する力と前記集中質量行列とに基づき分子動力学ソルバに次の時点（ｔ＋単位時間）における各節点の位置及び速度を算出させること、
を含み、
初期時点から目標の時点に到達するまで、前記（ｄ）（ｅ）及び（ｆ）の処理を、前記時点（ｔ）を単位時間経過させつつ繰り返し実行する。

このように、ＦＥＭモデルにおける節点を１粒子として取り扱い、ＦＥＭモデルに基づく要素剛性行列と初期位置に対する位置変位とに基づき節点に作用する力を算出するので、分子動力学ソルバを用いて時点（ｔ）における力と節点位置から次の時点（ｔ＋単位時間）の各節点の位置が算出でき、ＦＥＭの計算を、分子動力学ソルバを用いて実行可能となる。それでいて、分子動力学ソルバは無償利用可能なソフトを含めて並列計算機能を有するので、構造物のＦＥＭ解析をリーズナブルに並列計算で実現可能となる。

本発明の構造物のＦＥＭ解析システムを示す図システムが実行する処理のフローチャートＦＥＭモデルを示す図実施例のＦＥＭモデルの各節点の初期位置を示す図実施例の壁に衝突したときのＦＥＭモデルの各節点の位置を示す図

以下、本発明の一実施形態を、図面を参照して説明する。

［構造物のＦＥＭ解析システム］
本実施形態のシステム１は、構造物を有限要素法で解析するシステムである。

図１に示すように、システム１は、ＦＥＭモデル取得部１０と、集中質量行列取得部１１と、初期位置速度取得部１２と、要素剛性行列算出部１３と、変位行列算出部１４と、力算出部１５と、分子動力学ソルバ１６と、を有する。これら各部１０〜１５は、プロセッサ、メモリ、各種インターフェイス等を備えたコンピュータにおいて予め記憶されている図２に示す処理ルーチンをプロセッサが実行することによりソフトウェア及びハードウェアが協働して実現される。本実施形態では、１つの装置におけるプロセッサが各部の処理を実行しているが、これに限定されない。例えば、ネットワークを用いて分散させ、複数のプロセッサが各部の処理を実行するように構成してもよい。すなわち、１又は複数のプロセッサが処理を実行する。

一般的にＦＥＭ解析は、ＦＥＭソルバ（ソフトウェア）を用い、ＦＥＭソルバに対してＦＥＭモデルデータ、ＦＥＭモデルに与える荷重、境界条件などの各種データを渡し、ＦＥＭソルバを動作させる。ＦＥＭは構造物を複数の要素に分割して解析するが、要素を構成する複数の節点を有し、節点間に分子動力学でいう相互作用（力、ポテンシャル）が作用する形で節点に作用する力を演算している、と発明者は認識した。すなわち、力の算出をＦＥＭ的に実装すれば、分子動力学ソルバが利用可能であると考えた。

そこで、本実施形態のＦＥＭ解析法ではＦＥＭソルバを用いずに、代わりに、ＦＥＭモデルの要素を構成する１つの節点を１つの粒子として取り扱い、複数の粒子で構成される材料の分子シミュレーションを実行する分子動力学ソルバを用いている。分子動力学ソルバには、無償利用可能なソフトウェアとして、ＬＡＭＭＰＳ（Large-scale Atomic/Molecular Massively Parallel Simulator）やｇｒｏｍａｃｓ等があり、本実施形態ではＬＡＭＭＰＳを用いている。分子動力学ソルバに渡すデータは、材料を構成する複数の粒子の初期位置及び初期速度、各粒子に作用する相互作用（力、ポテンシャル）がある。相互作用は粒子間距離に応じて定まる。分子動力学ソルバは、ポテンシャル算出部と、粒子位置更新部とを有する。一般的な分子動力学ソルバの動作の流れの概要として、ポテンシャル算出部が、時点（ｔ）における各粒子の位置及び速度に基づき時点（ｔ）における各粒子に作用する力（ポテンシャル）を粒子間距離に応じて算出する。次に、粒子位置更新部が、時点（ｔ）における各粒子の位置及び速度、各粒子に作用する力（ポテンシャル）に基づき、次の時点（ｔ＋１）における各粒子の位置及び速度を算出する。このように、粒子に作用する力（ポテンシャル）の算出と、時点（ｔ）における各粒子の位置及び速度、各粒子に作用する力に基づく次の時点（ｔ＋１）における各粒子の位置及び速度の更新と、が繰り返し実行される。この繰り返し処理は、目標となる時点まで実行される。本実施形態において、分子動力学ソルバ１６のポテンシャル算出部が、ＦＥＭ解析ではそのまま利用できないため、代わりに力算出部１５を設けている。分子動力学ソルバ１６の粒子位置更新部はそのまま利用している。
以下、具体的に説明する。

図１に示すＦＥＭモデル取得部１０は、構造物を複数の要素で表現したＦＥＭモデルデータＤ１を取得する。ＦＥＭモデルデータＤ１は、各要素を区画する複数の節点により構造物の立体形状を示し、要素に物性値（ヤング率Ｅ、ポアソン比ｖ、質量密度ρ）が設定されている。図３の例示は、１要素が４つの節点で表現される四面体であり、要素ｅ＝１、２で示される２つの要素を例示している。要素ｅ＝１を構成する節点ｉ＝１〜４であり、要素ｅ＝２を構成する節点ｉ＝２〜５である。節点ｉはそれぞれｘｙｚ座標の値を有し且つどの節点と接続されているかとどの要素に属するかを示す情報を有する。ＦＥＭモデルデータＤ１は、一般の有限要素法シミュレーションにて多用されている。なお、図３に示す本実施形態では、ＦＥＭモデルデータＤ１の要素は四面体であるが、これに限定されない。原理的には任意の多面体要素が利用可能である。

図１に示す集中質量行列取得部１１は、構造物の質量を各節点に配分して示す集中質量行列Ｍ_ｅを取得する。集中質量行列取得部１１は、ＦＥＭモデルデータＤ１に集中質量行列Ｍ_ｅが定義されている場合は集中質量行列Ｍ_ｅをそのまま取得し、集中質量行列Ｍ_ｅが定義されていない場合には算出する。要素ｅの集中質量行列Ｍ_ｅは次の式（６）で表現される。集中質量行列取得部１１は、全ての要素について集中質量行列Ｍ_ｅを取得又は算出する。図３に示す実施例では、要素数が５０７０であるので、５０７０の集中質量行列Ｍ_ｅを取得する。集中質量行列Ｍ_ｅは、対角上にのみ質量を有する。標準的なＦＥＭ解析では非対角要素をもつ整合質量行列を使用することもあるが、分子動力学ソルバの大幅な変更を避けるために集中質量行列を用いている。

ただし、ρは構造物の質量密度を示し、Ｖ_ｅは要素ｅの体積を示し、Ｌ_ｅは要素ｅにおける節点数を示す。

図１に示す初期位置速度取得部１２は、ＦＥＭモデルデータＤ１に基づき各節点の初期位置及び初期速度を取得し、これらの値をワーキングメモリＤ２に記憶する。各節点の初期位置は、ＦＥＭモデルデータＤ１における各節点の座標（ｘｙｚ座標）である。各節点の初期速度は０である。図３に示す実施例では、要素数が５０７０で全節点数が１２０４であるので、１２０４個の節点の座標を取得する。

図１に示す要素剛性行列算出部１３は、各節点の初期位置及び物性値に基づき要素剛性行列Ｋ_ｅを要素毎に算出する。要素ｅの要素剛性行列Ｋ_ｅは、式（１）〜（３）で表現される。要素剛性行列算出部１３は、全ての要素について要素剛性行列Ｋ_ｅを算出する。図３に示す実施例では、要素数が５０７０であるので、５０７０の要素剛性行列Ｋ_ｅを算出する。

Ｖ_ｅは要素ｅの体積を示し、Ｌ_ｅは要素ｅにおける節点数を示し、Ｎ^ｅ _ｉ（ｒ）は要素ｅの節点ｉ（ｉ＝１，…，Ｌ_ｅ）の形状関数を示し、Ｅは構造物のヤング率を示し、ｖは構造物のポアソン比を示す。

図１に示す変位行列算出部１４は、各節点の初期位置に対する時点（ｔ）の位置の変位を示す変位行列ｄ_ｅを要素毎に算出する。変位行列ｄ_ｅは式（５）で表現される。変位行列算出部１４は、全ての要素について当該要素を構成する全節点の変位を示す変位行列ｄ_ｅを算出する。図３に示す実施例では、要素数が５０７０であるので、５０７０の変位行列ｄ_ｅを算出する。

ただし、Ｓ^ｅ _ｉは要素ｅの節点ｉにおけるｘ成分の変位を示し、Ｕ^ｅ _ｉは要素ｅの節点ｉにおけるｙ成分の変位を示し、Ｗ^ｅ _ｉは要素ｅの節点ｉにおけるｚ成分の変位を示す。

図１に示す力算出部１５は、要素剛性行列算出部１３が算出した要素剛性行列Ｋ_ｅと変位行列算出部１４が算出した変位行列ｄ_ｅとに基づき、要素ｅを構成する節点ｉ（ｉ＝１，…，Ｌ_ｅ）に作用する力（Ｆ^ｅ _ｉ，ｘ、Ｆ^ｅ _ｉ，ｙ、Ｆ^ｅ _ｉ，ｚ）を要素毎に算出する。要素ｅを構成する各節点ｉに作用する力Ｆ_ｅは式（４）で表される。

ただし、Ｆ^ｅ _ｉ，ｘは要素ｅの節点ｉに作用する力のｘ成分を示し、Ｆ^ｅ _ｉ，ｙは要素ｅの節点ｉに作用する力のｙ成分を示し、Ｆ^ｅ _ｉ，ｚは要素ｅの節点ｉに作用する力のｚ成分を示す。

図１に示す分子動力学ソルバ１６は、ワーキングメモリＤ２に記憶されている、時点（ｔ）における全ての要素を構成する各節点の位置及び速度と、各節点に作用する力と、集中質量行列Ｍ_ｅとに基づいて、次の時点（ｔ＋単位時間）における各節点の位置及び速度を算出する。算出された各節点の位置及び速度はワーキングメモリＤ２に記憶される。

［構造物のＦＥＭ解析方法］
図１に示すシステム１における１又は複数のプロセッサが実行する、構造物のＦＥＭ解析方法について、図２を用いて説明する。

まず、ステップＳＴ１０１において、ＦＥＭモデル取得部１０は、構造物を複数の要素で表現したＦＥＭモデルデータＤ１であって、各要素を区画する複数の節点および物性値（ヤング率Ｅ、ポアソン比ｖ、質量密度ρ等）が設定されたＦＥＭモデルデータＤ１を取得する。

次のステップＳＴ１０２において、集中質量行列取得部１１がＦＥＭモデルデータＤ１に基づき、構造物の質量を各節点に配分して表す集中質量行列Ｍ_ｅを要素毎に取得し、初期位置速度取得部１２が、ＦＥＭモデルデータＤ１に基づき各節点の初期位置及び初期速度を取得する。

ステップＳＴ１０３、ＳＴ１０４、ＳＴ１０５は、全ての要素（図３の例では、５０７０個の要素）について要素毎（ｅ＝１〜５０７０）に実行する。ステップＳＴ１０３、ＳＴ１０４、ＳＴ１０５の処理を全ての要素について実行すれば、ステップＳＴ１０６の処理へ移行する。

ステップＳＴ１０３において、要素剛性行列算出部１３は、要素ｅの各節点の初期位置及び物性値に基づき要素剛性行列を算出する。要素ｅの要素剛性行列Ｋ_ｅは、式（１）〜（３）で表現される。

ステップＳＴ１０４において、変位行列算出部１４は、時点（ｔ）における要素ｅの変位行列ｄ_ｅを算出する。変位行列ｄ_ｅは、要素ｅを構成する各節点ｉ（ｉ＝１〜Ｌ_ｅ）の初期位置に対する現在位置の変位を示す。変位行列ｄ_ｅは式（５）で表現される。

ステップＳＴ１０５において、力算出部１５は、要素ｅの要素剛性行列Ｋ_ｅと変位行列ｄ_ｅとに基づき、要素ｅの各節点に作用する力を算出する。要素ｅを構成する各節点ｉ（ｉ＝１，…，Ｌ_ｅ）に作用する力（Ｆ^ｅ _ｉ，ｘ、Ｆ^ｅ _ｉ，ｙ、Ｆ^ｅ _ｉ，ｚ）は、式（４）〜（５）で表現される

全ての要素についてステップＳＴ１０３、ＳＴ１０４、ＳＴ１０５の処理を実行した後に、次にステップＳＴ１０６において、分子動力学ソルバ１６は、時点（ｔ）における各節点の位置と速度と各節点に作用する力と集中質量行列Ｍ_ｅとに基づき次の時点（ｔ＋単位時間）における各節点の位置及び速度を算出する。

次のステップＳＴ１０７において、ＳＴ１０３〜ＳＴ１０６の処理を、初期時点から目標時点に到達するまで、時点（ｔ）を単位時間経過させつつ繰り返し実行したか否かを判定する。目標時点に到達していない場合（ＳＴ１０７：ＮＯ）には、時点（ｔ）を単位時間経過させて、ＳＴ１０３の処理に戻る。目標時点に到達した場合（ＳＴ１０７：ＹＥＳ）には、ＦＥＭ解析処理を終了する。

なお、本実施形態では、ステップＳＴ１０３における要素剛性行列Ｋ_ｅの算出は、単位時間経過する毎に実行されている。これは、要素剛性行列Ｋ_ｅは要素数５０７０あるため、メモリに記憶するとメモリ消費量が無視できないからであり、単位時間経過するたびに計算しなおしても、計算コストが低いためである。勿論、全ての要素についての要素剛性行列Ｋ_ｅを算出した時点でメモリに記憶しておき、単位時間経過するたびに算出しないようにしてもよい。

本発明のシミュレーション方法の結果を説明する。
図４及び図５は、半径１０ｍｍの弾性球の落下、反発運動をシミュレーションした結果である。図４は、初期位置のＦＥＭモデルの各節点位置を示し、図５は、壁との衝突時におけるＦＥＭモデルの各節点の位置を示す。物性値は、ヤング率Ｅが１ＭＰａ、ポアソン比ｖが０．４９、質量密度ρは９７０ｋｇ／ｍ^３であり、ゴムボールに相当する構造物である。要素タイプは、１要素が４つの節点を有する四面体１次要素であり、全体として５０７０要素、１２０４節点ある。弾性球の中心の高さを２０ｍｍとして重力により自然落下させ、高さ０ｍｍに剛体壁を設定している。時間ステップ数は３００００（約１．０秒に相当）である。

並列計算機能が実現できたことを示すために、プロセス数を１、２、４、８、１６と設定し、計算時間を次に示す。下記の通り、プロセス数に応じて計算時間が短縮されているので、ＦＥＭ解析に並列計算機能をリーズナブルに実装できていることがわかる。
プロセス数：１計算時間（秒）：１１０６．８
プロセス数：２計算時間（秒）：５６０．４
プロセス数：４計算時間（秒）：２８２．６
プロセス数：８計算時間（秒）：１５０．２
プロセス数：１６計算時間（秒）：８３．０

以上のように、本実施形態の構造物のＦＥＭ解析方法は、
１又は複数のプロセッサが実行する方法であって、
（ａ）構造物を複数の要素ｅで表現したＦＥＭモデルデータＤ１であって、各要素ｅを区画する複数の節点ｉおよび物性値（ヤング率Ｅ、ポアソン比ｖ、質量密度ρ）が設定されたＦＥＭモデルデータＤ１を取得すること（ＳＴ１０１）、
（ｂ）ＦＥＭモデルデータＤ１に基づき、構造物の質量を各節点に配分して表す集中質量行列Ｍ_ｅと、各節点の初期位置及び初期速度と、を取得すること（ＳＴ１０２）、
（ｃ）各節点の初期位置及び物性値に基づき要素剛性行列Ｋ_ｅを要素毎に算出すること（ＳＴ１０３）、
（ｄ）各節点の初期位置に対する時点（ｔ）の位置の変位を表す変位行列ｄ_ｅを要素毎に算出すること（ＳＴ１０４）、
（ｅ）要素剛性行列Ｋ_ｅと変位行列ｄ_ｅとに基づき、各節点に作用する力Ｆ_ｅを要素毎に算出すること（ＳＴ１０５）、
（ｆ）時点（ｔ）における前記各節点の位置と速度と前記各節点に作用する力と前記集中質量行列とに基づき分子動力学ソルバに次の時点（ｔ＋単位時間）における各節点の位置及び速度を算出させること（ＳＴ１０６）、
を含み、
初期時点から目標の時点に到達するまで、（ｄ）（ｅ）及び（ｆ）の処理を、時点（ｔ）を単位時間経過させつつ繰り返し実行する。

本実施形態の構造物のＦＥＭ解析システムは、
構造物を複数の要素ｅで表現したＦＥＭモデルデータＤ１であって、各要素ｅを区画する複数の節点ｉおよび物性値（ヤング率Ｅ、ポアソン比ｖ、質量密度ρ）が設定されたＦＥＭモデルデータＤ１を取得するＦＥＭモデル取得部１０と、
ＦＥＭモデルデータＤ１に基づき、構造物の質量を各節点に配分して表す集中質量行列Ｍ_ｅを取得する集中質量行列取得部１１と、
ＦＥＭモデルデータＤ１に基づき、各節点の初期位置及び初期速度と、を取得する初期位置速度取得部１２と、
各節点の初期位置及び物性値に基づき要素剛性行列Ｋ_ｅを要素毎に算出する要素剛性行列算出部１３と、
各節点の初期位置に対する時点（ｔ）の位置の変位を表す変位行列ｄ_ｅを要素毎に算出する変位行列算出部１４と、
要素剛性行列Ｋ_ｅと変位行列ｄ_ｅとに基づき、各節点に作用する力Ｆ_ｅを要素毎に算出する力算出部１５と、
時点（ｔ）における前記各節点の位置と速度と前記各節点に作用する力と前記集中質量行列とに基づき次の時点（ｔ＋単位時間）における各節点の位置及び速度を算出する分子動力学ソルバ１６と、
を備え、
初期時点から目標の時点に到達するまで、変位行列算出部１４、力算出部１５及び分子動力学ソルバ１６による処理を、時点（ｔ）を単位時間経過させつつ繰り返し実行するように構成されている。

本実施形態のように、前記（ｃ）における要素剛性行列算出部１３による要素剛性行列Ｋ_ｅの算出は、単位時間経過する毎に算出されることが好ましい。

このようにすれば、要素剛性行列Ｋ_ｅを単位時間経過する毎に毎回算出しても計算コストが大きくなく、全ての要素の要素剛性行列Ｋ_ｅをメモリに記憶することによって生じるメモリの消費を抑制することが可能となる。

本実施形態のように、要素ｅの要素剛性行列Ｋ_ｅは、式（１）〜（３）で表現されることが好ましい。好適な実施形態である。

本実施形態のように、要素ｅを構成する各節点ｉ（ｉ＝１，…，Ｌ_ｅ）に作用する力（Ｆ^ｅ _ｉ，ｘ、Ｆ^ｅ _ｉ，ｙ、Ｆ^ｅ _ｉ，ｚ）は、式（４）〜（５）で表現されることが好ましい。好適な実施形態である。

本実施形態に係るプログラムは、上記方法をコンピュータに実行させるプログラムである。このプログラムを実行することによっても、上記方法の奏する作用効果を得ることが可能となる。

以上、本発明の実施形態について図面に基づいて説明したが、具体的な構成は、これらの実施形態に限定されるものでないと考えられるべきである。本発明の範囲は、上記した実施形態の説明だけではなく特許請求の範囲によって示され、さらに特許請求の範囲と均等の意味および範囲内でのすべての変更が含まれる。

例えば、図１に示す各部１０〜１５は、所定プログラムをコンピュータのプロセッサで実行することで実現しているが、各部を専用回路で構成してもよい。また、本実施形態では１つのコンピュータにおけるプロセッサが各部１０〜１５を実装しているが、少なくとも１又は複数のプロセッサに分散して実装してもよい。

上記の各実施形態で採用している構造を他の任意の実施形態に採用することは可能である。各部の具体的な構成は、上述した実施形態のみに限定されるものではなく、本発明の趣旨を逸脱しない範囲で種々変形が可能である。

１０ＦＥＭモデル取得部
１１集中質量行列取得部
１２初期位置速度取得部
１３要素剛性行列算出部
１４変位行列算出部
１５力算出部
１６分子動力学ソルバ

Claims

１又は複数のプロセッサが実行する方法であって、
（ａ）構造物を複数の要素で表現したＦＥＭモデルデータであって、各要素を区画する複数の節点および物性値が設定されたＦＥＭモデルデータを取得すること、
（ｂ）前記ＦＥＭモデルデータに基づき、前記構造物の質量を各節点に配分して表す集中質量行列と、各節点の初期位置及び初期速度と、を取得すること、
（ｃ）前記各節点の初期位置及び前記物性値に基づき要素剛性行列を要素毎に算出すること、
（ｄ）各節点の初期位置に対する時点（ｔ）の位置の変位を表す変位行列を要素毎に算出すること、
（ｅ）前記要素剛性行列と前記変位行列とに基づき、各節点に作用する力を要素毎に算出すること、
（ｆ）時点（ｔ）における前記各節点の位置と速度と前記各節点に作用する力と前記集中質量行列とに基づき分子動力学ソルバに次の時点（ｔ＋単位時間）における各節点の位置及び速度を算出させること、
を含み、
初期時点から目標の時点に到達するまで、前記（ｄ）（ｅ）及び（ｆ）の処理を、前記時点（ｔ）を単位時間経過させつつ繰り返し実行する、構造物のＦＥＭ解析方法。
前記（ｃ）における前記要素剛性行列の算出は、前記単位時間経過する毎に算出される、請求項１に記載の方法。
要素ｅの要素剛性行列Ｋ_ｅは、式（１）〜（３）で表現される、請求項１又は２に記載の方法。

Ｖ_ｅは要素ｅの体積を示し、Ｌ_ｅは要素ｅにおける節点数を示し、Ｎ^ｅ _ｉ（ｒ）は要素ｅの節点ｉ（ｉ＝１，…，Ｌ_ｅ）の形状関数を示し、Ｅは構造物のヤング率を示し、ｖは構造物のポアソン比を示す。
要素ｅを構成する各節点ｉ（ｉ＝１，…，Ｌ_ｅ）に作用する力（Ｆ^ｅ _ｉ，ｘ、Ｆ^ｅ _ｉ，ｙ、Ｆ^ｅ _ｉ，ｚ）は、式（４）〜（５）で表現される、請求項３に記載の方法。

ただし、Ｓ^ｅ _ｉは要素ｅの節点ｉにおけるｘ成分の変位を示し、Ｕ^ｅ _ｉは要素ｅの節点ｉにおけるｙ成分の変位を示し、Ｗ^ｅ _ｉは要素ｅの節点ｉにおけるｚ成分の変位を示し、Ｆ^ｅ _ｉ，ｘは要素ｅの節点ｉに作用する力のｘ成分を示し、Ｆ^ｅ _ｉ，ｙは要素ｅの節点ｉに作用する力のｙ成分を示し、Ｆ^ｅ _ｉ，ｚは要素ｅの節点ｉに作用する力のｚ成分を示す。
構造物を複数の要素で表現したＦＥＭモデルデータであって、各要素を区画する複数の節点および物性値が設定されたＦＥＭモデルデータを取得するＦＥＭモデル取得部と、
前記ＦＥＭモデルデータに基づき、前記構造物の質量を各節点に配分して表す集中質量行列を取得する集中質量行列取得部と、
前記ＦＥＭモデルデータに基づき、各節点の初期位置及び初期速度と、を取得する初期位置速度取得部と、
前記各節点の初期位置及び前記物性値に基づき要素剛性行列を要素毎に算出する要素剛性行列算出部と、
各節点の初期位置に対する時点（ｔ）の位置の変位を表す変位行列を要素毎に算出する変位行列算出部と、
前記要素剛性行列と前記変位行列とに基づき、各節点に作用する力を要素毎に算出する力算出部と、
時点（ｔ）における前記各節点の位置と速度と前記各節点に作用する力と前記集中質量行列とに基づき次の時点（ｔ＋単位時間）における各節点の位置及び速度を算出する分子動力学ソルバと、
を備え、
初期時点から目標の時点に到達するまで、前記変位行列算出部、前記力算出部及び前記分子動力学ソルバによる処理を、前記時点（ｔ）を単位時間経過させつつ繰り返し実行するように構成されている、構造物のＦＥＭ解析システム。
前記要素剛性行列算出部による前記要素剛性行列の算出は、前記単位時間経過する毎に算出される、請求項５に記載のシステム。
要素ｅの要素剛性行列Ｋ_ｅは、式（１）〜（３）で表現される、請求項５又は６に記載のシステム。

Ｖ_ｅは要素ｅの体積を示し、Ｌ_ｅは要素ｅにおける節点数を示し、Ｎ^ｅ _ｉ（ｒ）は要素ｅの節点ｉ（ｉ＝１，…，Ｌ_ｅ）の形状関数を示し、Ｅは構造物のヤング率を示し、ｖは構造物のポアソン比を示す。
要素ｅを構成する各節点ｉ（ｉ＝１，…，Ｌ_ｅ）に作用する力（Ｆ^ｅ _ｉ，ｘ、Ｆ^ｅ _ｉ，ｙ、Ｆ^ｅ _ｉ，ｚ）は、式（４）〜（５）で表現される請求項７に記載のシステム。

ただし、Ｓ^ｅ _ｉは要素ｅの節点ｉにおけるｘ成分の変位を示し、Ｕ^ｅ _ｉは要素ｅの節点ｉにおけるｙ成分の変位を示し、Ｗ^ｅ _ｉは要素ｅの節点ｉにおけるｚ成分の変位を示し、Ｆ^ｅ _ｉ，ｘは要素ｅの節点ｉに作用する力のｘ成分を示し、Ｆ^ｅ _ｉ，ｙは要素ｅの節点ｉに作用する力のｙ成分を示し、Ｆ^ｅ _ｉ，ｚは要素ｅの節点ｉに作用する力のｚ成分を示す。
請求項１〜４のいずれかに記載の方法を１又は複数のプロセッサに実行させるプログラム。