(57)【要約】      （修正有）【課題】回転による幾何学的非線形に対応し且つ分子動力学ソルバを用いた構造物のＦＥＭ解析方法を提供する。【解決手段】構造物を複数の四面体要素ｅで表現したＦＥＭモデルデータに基づき集中質量行列と要素剛性行列Ｋｅを要素毎に算出し、第１節点を基準とした変形行列Ａを算出し、右極分解により歪変形行列Ｑと回転変形行列Ｒを算出し、歪変形行列Ｑを用いて回転を伴わない状態での変位ベクトルを算出し、要素剛性行列と変位ベクトルとに基づき歪のみに起因する力を算出し、回転変形行列Ｒを用いて、回転及び歪を伴った各節点に作用する力を算出し、時点（ｔ）における各節点の位置と速度と回転及び歪を伴った各節点に作用する力と集中質量行列とに基づき分子動力学ソルバに次の時点における各節点の位置及び速度を算出させる。初期時点から目標の時点に到達するまで、各処理を繰り返す。【選択図】図７
2020177388.tif 000002
【特許請求の範囲】
【請求項１】
１又は複数のプロセッサが実行する方法であって、
（ａ）構造物を複数の四面体要素で表現したＦＥＭモデルデータであって、各要素を区画する４つの節点および物性値が設定されたＦＥＭモデルデータを取得すること、
（ｂ）前記ＦＥＭモデルデータに基づき、前記構造物の質量を各節点に配分して表す集中質量行列と、各節点の初期位置及び初期速度と、を取得すること、
（ｃ）前記各節点の初期位置及び前記物性値に基づき要素剛性行列を要素毎に算出すること、
（ｄ１）初期位置における第１節点を基準とした第２〜４節点の相対位置ベクトルと、時点（ｔ）における第１節点を基準とした第２〜４節点の相対位置ベクトルと、に基づき各節点の変形を示す変形行列Ａを算出すること、
（ｄ２）右極分解により前記変形行列Ａから、歪による変形を表す歪変形行列Ｑと、回転による変形を示す回転変形行列Ｒと、を算出すること、
（ｄ３）前記歪変形行列Ｑを用いて、回転を伴わない状態での歪に起因する、各節点の初期位置に対する変位を表す変位ベクトルを要素毎に算出すること、
（ｅ１）前記要素剛性行列と前記変位ベクトルとに基づき、歪のみに起因して各節点に作用する力を要素毎に算出すること、
（ｅ２）前記回転変形行列Ｒを用いて前記歪のみに起因して各節点に作用する力の向きを変換することで、回転及び歪を伴った各節点に作用する力を算出すること、
（ｆ）時点（ｔ）における前記各節点の位置と速度と前記回転及び歪を伴った各節点に作用する力と前記集中質量行列とに基づき分子動力学ソルバに次の時点（ｔ＋単位時間）における各節点の位置及び速度を算出させること、
を含み、
初期時点から目標の時点に到達するまで、前記（ｄ１）（ｄ２）（ｄ３）（ｅ１）（ｅ２）及び（ｆ）の処理を、前記時点（ｔ）を単位時間経過させつつ繰り返し実行する、構造物のＦＥＭ解析方法。
【請求項２】
前記（ｃ）における前記要素剛性行列の算出は、前記単位時間経過する毎に算出される、請求項１に記載の方法。
【請求項３】
要素ｅの要素剛性行列Ｋｅは、式（１）〜（３）で表現される、請求項１又は２に記載の方法。
【数１】
2020177388.tif 000012
Ｖｅは要素ｅの体積を示し、Ｌｅは要素ｅにおける節点数を示し、Ｎｅｉ（ｒ）は要素ｅの節点ｉ（ｉ＝１，…，Ｌｅ）の形状関数を示し、Ｅは構造物のヤング率を示し、ｖは構造物のポアソン比を示す。
【請求項４】
前記変形行列Ａは、式（９）で表現される、請求項１〜３のいずれかに記載の方法。
【数２】
2020177388.tif 000013
ただし、初期位置における第１〜４節点のｘｙｚ座標が、それぞれ（ｘ１、ｙ１、ｚ１）、（ｘ２、ｙ２、ｚ２）、（ｘ３、ｙ３、ｚ３）及び（ｘ４、ｙ４、ｚ４）であり、時刻（ｔ）における第１〜４節点のｘｙｚ座標が、それぞれ（ｘ’１、ｙ’１、ｚ’１）、（ｘ’２、ｙ’２、ｚ’２）、（ｘ’３、ｙ’３、ｚ’３）及び（ｘ’４、ｙ’４、ｚ’４）である。
【請求項５】
前記回転を伴わない状態での歪に起因する、各節点の初期位置に対する変位を表す変位ベクトル（ｓｅ、ｕｅ、ｗｅ）は、式（１１）で表現される、請求項１〜４のいずれかに記載の方法。
【数３】
2020177388.tif 000014
ただし、変位ベクトルｓｅは要素ｅのｘ成分の変位を示し、変位ベクトルｕｅは要素ｅのｙ成分の変位を示し、変位ベクトルｗｅは要素ｅのｚ成分の変位を示す。Ｓｅｉは要素ｅの節点ｉにおけるｘ成分の変位を示し、Ｕｅｉは要素ｅの節点ｉにおけるｙ成分の変位を示し、Ｗｅｉは要素ｅの節点ｉにおけるｚ成分の変位を示し、初期位置における第１〜４節点のｘｙｚ座標が、それぞれ（ｘ１、ｙ１、ｚ１）、（ｘ２、ｙ２、ｚ２）、（ｘ３、ｙ３、ｚ３）及び（ｘ４、ｙ４、ｚ４）である。
【請求項６】
前記歪のみに起因して各節点に作用する力ｆｅは、式（１２）で表現される、請求項５に記載の方法。
【数４】
2020177388.tif 000015
ただし、ｆｅｉ，ｘは要素ｅの節点ｉに作用する力のｘ成分を示し、ｆｅｉ，ｙは要素ｅの節点ｉに作用する力のｙ成分を示し、ｆｅｉ，ｚは要素ｅの節点ｉに作用する力のｚ成分を示す。
【請求項７】
前記回転及び歪を伴った各節点に作用する力Ｆｅは、式（１３）で表現される、請求項６に記載の方法。
【数５】
2020177388.tif 000016
ただし、Ｆｅｉ，ｘは要素ｅの節点ｉに作用する回転及び歪を伴った力のｘ成分を示し、Ｆｅｉ，ｙは要素ｅの節点ｉに作用する回転及び歪を伴った力のｙ成分を示し、Ｆｅｉ，ｚは要素ｅの節点ｉに作用する回転及び歪を伴った力のｚ成分を示す。
【請求項８】
構造物を複数の四面体要素で表現したＦＥＭモデルデータであって、各要素を区画する４つの節点および物性値が設定されたＦＥＭモデルデータを取得するＦＥＭモデル取得部と、
前記ＦＥＭモデルデータに基づき、前記構造物の質量を各節点に配分して表す集中質量行列を取得する集中質量行列取得部と、
前記ＦＥＭモデルデータに基づき、各節点の初期位置及び初期速度を取得する初期位置速度取得部と、
前記各節点の初期位置及び前記物性値に基づき要素剛性行列を要素毎に算出する要素剛性行列算出部と、
初期位置における第１節点を基準とした第２〜４節点の相対位置ベクトルと、時点（ｔ）における第１節点を基準とした第２〜４節点の相対位置ベクトルと、に基づき各節点の変形を示す変形行列Ａを算出する変形行列算出部と、
右極分解により前記変形行列Ａから、歪による変形を表す歪変形行列Ｑと、回転による変形を示す回転変形行列Ｒと、を算出する右極分解部と、
前記歪変形行列Ｑを用いて、回転を伴わない状態での歪に起因する、各節点の初期位置に対する変位を表す変位ベクトルを要素毎に算出する変位ベクトル算出部と、
前記要素剛性行列と前記変位ベクトルとに基づき、歪のみに起因して各節点に作用する力を要素毎に算出する力算出部と、
前記回転変形行列Ｒを用いて前記歪のみに起因して各節点に作用する力の向きを変換することで、回転及び歪を伴った各節点に作用する力を算出する力方向変換部と、
時点（ｔ）における前記各節点の位置と速度と前記回転及び歪を伴った各節点に作用する力と前記集中質量行列とに基づき次の時点（ｔ＋単位時間）における各節点の位置及び速度を算出する分子動力学ソルバと、
を備え、
初期時点から目標の時点に到達するまで、前記変形行列算出部、前記右極分解部、前記変位ベクトル算出部、前記力算出部、前記力方向変換部及び前記分子動力学ソルバによる処理を、前記時点（ｔ）を単位時間経過させつつ繰り返し実行するように構成されている、構造物のＦＥＭ解析システム。
【請求項９】
前記要素剛性行列算出部による前記要素剛性行列の算出は、前記単位時間経過する毎に算出される、請求項８に記載のシステム。
【請求項１０】
要素ｅの要素剛性行列Ｋｅは、式（１）〜（３）で表現される、請求項８又は９に記載のシステム。
【数６】
2020177388.tif 000017
Ｖｅは要素ｅの体積を示し、Ｌｅは要素ｅにおける節点数を示し、Ｎｅｉ（ｒ）は要素ｅの節点ｉ（ｉ＝１，…，Ｌｅ）の形状関数を示し、Ｅは構造物のヤング率を示し、ｖは構造物のポアソン比を示す。
【請求項１１】
前記変形行列Ａは、式（９）で表現される、請求項８〜１０のいずれかに記載の方法。
【数７】
2020177388.tif 000018
ただし、初期位置における第１〜４節点のｘｙｚ座標が、それぞれ（ｘ１、ｙ１、ｚ１）、（ｘ２、ｙ２、ｚ２）、（ｘ３、ｙ３、ｚ３）及び（ｘ４、ｙ４、ｚ４）であり、時刻（ｔ）における第１〜４節点のｘｙｚ座標が、それぞれ（ｘ’１、ｙ’１、ｚ’１）、（ｘ’２、ｙ’２、ｚ’２）、（ｘ’３、ｙ’３、ｚ’３）及び（ｘ’４、ｙ’４、ｚ’４）である。
【請求項１２】
前記回転を伴わない状態での歪に起因する、各節点の初期位置に対する変位を表す変位ベクトル（ｓｅ、ｕｅ、ｗｅ）は、式（１１）で表現される、請求項８〜１１のいずれかに記載の方法。
【数８】
2020177388.tif 000019
ただし、変位ベクトルｓｅは要素ｅのｘ成分の変位を示し、変位ベクトルｕｅは要素ｅのｙ成分の変位を示し、変位ベクトルｗｅは要素ｅのｚ成分の変位を示す。Ｓｅｉは要素ｅの節点ｉにおけるｘ成分の変位を示し、Ｕｅｉは要素ｅの節点ｉにおけるｙ成分の変位を示し、Ｗｅｉは要素ｅの節点ｉにおけるｚ成分の変位を示し、初期位置における第１〜４節点のｘｙｚ座標が、それぞれ（ｘ１、ｙ１、ｚ１）、（ｘ２、ｙ２、ｚ２）、（ｘ３、ｙ３、ｚ３）及び（ｘ４、ｙ４、ｚ４）である。
【請求項１３】
前記歪のみに起因して各節点に作用する力ｆｅは、式（１２）で表現される、請求項１２に記載の方法。
【数９】
2020177388.tif 000020
ただし、ｆｅｉ，ｘは要素ｅの節点ｉに作用する力のｘ成分を示し、ｆｅｉ，ｙは要素ｅの節点ｉに作用する力のｙ成分を示し、ｆｅｉ，ｚは要素ｅの節点ｉに作用する力のｚ成分を示す。
【請求項１４】
前記回転及び歪を伴った各節点に作用する力Ｆｅは、式（１３）で表現される、請求項１３に記載の方法。
【数１０】
2020177388.tif 000021
ただし、Ｆｅｉ，ｘは要素ｅの節点ｉに作用する回転及び歪を伴った力のｘ成分を示し、Ｆｅｉ，ｙは要素ｅの節点ｉに作用する回転及び歪を伴った力のｙ成分を示し、Ｆｅｉ，ｚは要素ｅの節点ｉに作用する回転及び歪を伴った力のｚ成分を示す。
【請求項１５】
請求項１〜７のいずれかに記載の方法を１又は複数のプロセッサに実行させるプログラム。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、構造物のＦＥＭ解析方法、システム及びプログラムに関する。
【背景技術】
【０００２】
構造物の動解析には有限要素法（ＦＥＭ；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）やｇｒｏｍａｃｓ等があり、本実施形態ではＬＡＭＭＰＳを用いている。分子動力学ソルバに渡すデータは、材料を構成する複数の粒子の初期位置及び初期速度、各粒子に作用する相互作用（力、ポテンシャル）がある。相互作用は粒子間距離に応じて定まる。分子動力学ソルバは、ポテンシャル算出部と、粒子位置更新部とを有する。一般的な分子動力学ソルバの動作の流れの概要として、ポテンシャル算出部が、時点（ｔ）における各粒子の位置及び速度に基づき時点（ｔ）における各粒子に作用する力（ポテンシャル）を粒子間距離に応じて算出する。次に、粒子位置更新部が、時点（ｔ）における各粒子の位置及び速度、各粒子に作用する力（ポテンシャル）に基づき、次の時点（ｔ＋１）における各粒子の位置及び速度を算出する。このように、粒子に作用する力（ポテンシャル）の算出と、時点（ｔ）における各粒子の位置及び速度、各粒子に作用する力に基づく次の時点（ｔ＋１）における各粒子の位置及び速度の更新と、が繰り返し実行される。この繰り返し処理は、目標となる時点まで実行される。本実施形態において、分子動力学ソルバ１６のポテンシャル算出部が、ＦＥＭ解析ではそのまま利用できないため、代わりに力算出部１５を設けている。分子動力学ソルバ１６の粒子位置更新部はそのまま利用している。
以下、具体的に説明する。
【００１７】
図１に示すＦＥＭモデル取得部１０は、構造物を複数の要素で表現したＦＥＭモデルデータＤ１を取得する。ＦＥＭモデルデータＤ１は、各要素を区画する複数の節点により構造物の立体形状を示し、要素に物性値（ヤング率Ｅ、ポアソン比ｖ、質量密度ρ）が設定されている。図３の例示は、１要素が４つの節点で表現される四面体であり、要素ｅ＝１、２で示される２つの要素を例示している。要素ｅ＝１を構成する節点ｉ＝１〜４であり、要素ｅ＝２を構成する節点ｉ＝２〜５である。節点ｉはそれぞれｘｙｚ座標の値を有し且つどの節点と接続されているかとどの要素に属するかを示す情報を有する。ＦＥＭモデルデータＤ１は、一般の有限要素法シミュレーションにて多用されている。なお、図３に示す本実施形態では、ＦＥＭモデルデータＤ１の要素は四面体であるが、これに限定されない。原理的には任意の多面体要素が利用可能である。
【００１８】
図１に示す集中質量行列取得部１１は、構造物の質量を各節点に配分して示す集中質量行列Ｍｅを取得する。集中質量行列取得部１１は、ＦＥＭモデルデータＤ１に集中質量行列Ｍｅが定義されている場合は集中質量行列Ｍｅをそのまま取得し、集中質量行列Ｍｅが定義されていない場合には算出する。要素ｅの集中質量行列Ｍｅは次の式（６）で表現される。集中質量行列取得部１１は、全ての要素について集中質量行列Ｍｅを取得又は算出する。図３に示す実施例では、要素数が５０７０であるので、５０７０の集中質量行列Ｍｅを取得する。集中質量行列Ｍｅは、対角上にのみ質量を有する。標準的なＦＥＭ解析では非対角要素をもつ整合質量行列を使用することもあるが、分子動力学ソルバの大幅な変更を避けるために集中質量行列を用いている。
【数１】
2020177388.tif 000003
ただし、ρは構造物の質量密度を示し、Ｖｅは要素ｅの体積を示し、Ｌｅは要素ｅにおける節点数を示す。
【００１９】
図１に示す初期位置速度取得部１２は、ＦＥＭモデルデータＤ１に基づき各節点の初期位置及び初期速度を取得し、これらの値をワーキングメモリＤ２に記憶する。各節点の初期位置は、ＦＥＭモデルデータＤ１における各節点の座標（ｘｙｚ座標）である。各節点の初期速度は０である。図３に示す実施例では、要素数が５０７０で全節点数が１２０４であるので、１２０４個の節点の座標を取得する。
【００２０】
図１に示す要素剛性行列算出部１３は、各節点の初期位置及び物性値に基づき要素剛性行列Ｋｅを要素毎に算出する。要素ｅの要素剛性行列Ｋｅは、式（１）〜（３）で表現される。要素剛性行列算出部１３は、全ての要素について要素剛性行列Ｋｅを算出する。図３に示す実施例では、要素数が５０７０であるので、５０７０の要素剛性行列Ｋｅを算出する。
【数２】
2020177388.tif 000004
Ｖｅは要素ｅの体積を示し、Ｌｅは要素ｅにおける節点数を示し、Ｎｅｉ（ｒ）は要素ｅの節点ｉ（ｉ＝１，…，Ｌｅ）の形状関数を示し、Ｅは構造物のヤング率を示し、ｖは構造物のポアソン比を示す。
【００２１】
図１に示す変位行列算出部１４は、各節点の初期位置に対する時点（ｔ）の位置の変位を示す変位行列ｄｅを要素毎に算出する。変位行列ｄｅは式（５）で表現される。変位行列算出部１４は、全ての要素について当該要素を構成する全節点の変位を示す変位行列ｄｅを算出する。図３に示す実施例では、要素数が５０７０であるので、５０７０の変位行列ｄｅを算出する。
【数３】
2020177388.tif 000005
ただし、Ｓｅｉは要素ｅの節点ｉにおけるｘ成分の変位を示し、Ｕｅｉは要素ｅの節点ｉにおけるｙ成分の変位を示し、Ｗｅｉは要素ｅの節点ｉにおけるｚ成分の変位を示す。
【００２２】
図１に示す力算出部１５は、要素剛性行列算出部１３が算出した要素剛性行列Ｋｅと変位行列算出部１４が算出した変位行列ｄｅとに基づき、要素ｅを構成する節点ｉ（ｉ＝１，…，Ｌｅ）に作用する力（Ｆｅｉ，ｘ、Ｆｅｉ，ｙ、Ｆｅｉ，ｚ）を要素毎に算出する。要素ｅを構成する各節点ｉに作用する力Ｆｅは式（４）で表される。
【数４】
2020177388.tif 000006
ただし、Ｆｅｉ，ｘは要素ｅの節点ｉに作用する力のｘ成分を示し、Ｆｅｉ，ｙは要素ｅの節点ｉに作用する力のｙ成分を示し、Ｆｅｉ，ｚは要素ｅの節点ｉに作用する力のｚ成分を示す。
【００２３】
図１に示す分子動力学ソルバ１６は、ワーキングメモリＤ２に記憶されている、時点（ｔ）における全ての要素を構成する各節点の位置及び速度と、各節点に作用する力と、集中質量行列Ｍｅとに基づいて、次の時点（ｔ＋単位時間）における各節点の位置及び速度を算出する。算出された各節点の位置及び速度はワーキングメモリＤ２に記憶される。
【００２４】
［構造物のＦＥＭ解析方法］
図１に示すシステム１における１又は複数のプロセッサが実行する、構造物のＦＥＭ解析方法について、図２を用いて説明する。
【００２５】
まず、ステップＳＴ１０１において、ＦＥＭモデル取得部１０は、構造物を複数の要素で表現したＦＥＭモデルデータＤ１であって、各要素を区画する複数の節点および物性値（ヤング率Ｅ、ポアソン比ｖ、質量密度ρ等）が設定されたＦＥＭモデルデータＤ１を取得する。
【００２６】
次のステップＳＴ１０２において、集中質量行列取得部１１がＦＥＭモデルデータＤ１に基づき、構造物の質量を各節点に配分して表す集中質量行列Ｍｅを要素毎に取得し、初期位置速度取得部１２が、ＦＥＭモデルデータＤ１に基づき各節点の初期位置及び初期速度を取得する。
【００２７】
ステップＳＴ１０３、ＳＴ１０４、ＳＴ１０５は、全ての要素（図３の例では、５０７０個の要素）について要素毎（ｅ＝１〜５０７０）に実行する。ステップＳＴ１０３、ＳＴ１０４、ＳＴ１０５の処理を全ての要素について実行すれば、ステップＳＴ１０６の処理へ移行する。
【００２８】
ステップＳＴ１０３において、要素剛性行列算出部１３は、要素ｅの各節点の初期位置及び物性値に基づき要素剛性行列を算出する。要素ｅの要素剛性行列Ｋｅは、式（１）〜（３）で表現される。
【００２９】
ステップＳＴ１０４において、変位行列算出部１４は、時点（ｔ）における要素ｅの変位行列ｄｅを算出する。変位行列ｄｅは、要素ｅを構成する各節点ｉ（ｉ＝１〜Ｌｅ）の初期位置に対する現在位置の変位を示す。変位行列ｄｅは式（５）で表現される。
【００３０】
ステップＳＴ１０５において、力算出部１５は、要素ｅの要素剛性行列Ｋｅと変位行列ｄｅとに基づき、要素ｅの各節点に作用する力を算出する。要素ｅを構成する各節点ｉ（ｉ＝１，…，Ｌｅ）に作用する力（Ｆｅｉ，ｘ、Ｆｅｉ，ｙ、Ｆｅｉ，ｚ）は、式（４）〜（５）で表現される
【００３１】
全ての要素についてステップＳＴ１０３、ＳＴ１０４、ＳＴ１０５の処理を実行した後に、次にステップＳＴ１０６において、分子動力学ソルバ１６は、時点（ｔ）における各節点の位置と速度と各節点に作用する力と集中質量行列Ｍｅとに基づき次の時点（ｔ＋単位時間）における各節点の位置及び速度を算出する。
【００３２】
次のステップＳＴ１０７において、ＳＴ１０３〜ＳＴ１０６の処理を、初期時点から目標時点に到達するまで、時点（ｔ）を単位時間経過させつつ繰り返し実行したか否かを判定する。目標時点に到達していない場合（ＳＴ１０７：ＮＯ）には、時点（ｔ）を単位時間経過させて、ＳＴ１０３の処理に戻る。目標時点に到達した場合（ＳＴ１０７：ＹＥＳ）には、ＦＥＭ解析処理を終了する。
【００３３】
なお、本実施形態では、ステップＳＴ１０３における要素剛性行列Ｋｅの算出は、単位時間経過する毎に実行されている。これは、要素剛性行列Ｋｅは要素数５０７０あるため、メモリに記憶するとメモリ消費量が無視できないからであり、単位時間経過するたびに計算しなおしても、計算コストが低いためである。勿論、全ての要素についての要素剛性行列Ｋｅを算出した時点でメモリに記憶しておき、単位時間経過するたびに算出しないようにしてもよい。
【００３４】
第１実施形態のシミュレーション方法の結果を説明する。
図４及び図５は、半径１０ｍｍの弾性球の落下、反発運動をシミュレーションした結果である。図４は、初期位置のＦＥＭモデルの各節点位置を示し、図５は、壁との衝突時におけるＦＥＭモデルの各節点の位置を示す。物性値は、ヤング率Ｅが１ＭＰａ、ポアソン比ｖが０．４９、質量密度ρは９７０ｋｇ／ｍ３であり、ゴムボールに相当する構造物である。要素タイプは、１要素が４つの節点を有する四面体１次要素であり、全体として５０７０要素、１２０４節点ある。弾性球の中心の高さを２０ｍｍとして重力により自然落下させ、高さ０ｍｍに剛体壁を設定している。時間ステップ数は３００００（約１．０秒に相当）である。
【００３５】
並列計算機能が実現できたことを示すために、プロセス数を１、２、４、８、１６と設定し、計算時間を次に示す。下記の通り、プロセス数に応じて計算時間が短縮されているので、ＦＥＭ解析に並列計算機能をリーズナブルに実装できていることがわかる。
プロセス数：１      計算時間（秒）：１１０６．８
プロセス数：２      計算時間（秒）：５６０．４
プロセス数：４      計算時間（秒）：２８２．６
プロセス数：８      計算時間（秒）：１５０．２
プロセス数：１６    計算時間（秒）：８３．０
【００３６】
以上のように、第１実施形態の構造物のＦＥＭ解析方法は、
１又は複数のプロセッサが実行する方法であって、
（ａ）構造物を複数の要素ｅで表現したＦＥＭモデルデータＤ１であって、各要素ｅを区画する複数の節点ｉおよび物性値（ヤング率Ｅ、ポアソン比ｖ、質量密度ρ）が設定されたＦＥＭモデルデータＤ１を取得すること（ＳＴ１０１）、
（ｂ）ＦＥＭモデルデータＤ１に基づき、構造物の質量を各節点に配分して表す集中質量行列Ｍｅと、各節点の初期位置及び初期速度と、を取得すること（ＳＴ１０２）、
（ｃ）各節点の初期位置及び物性値に基づき要素剛性行列Ｋｅを要素毎に算出すること（ＳＴ１０３）、
（ｄ）各節点の初期位置に対する時点（ｔ）の位置の変位を表す変位行列ｄｅを要素毎に算出すること（ＳＴ１０４）、
（ｅ）要素剛性行列Ｋｅと変位行列ｄｅとに基づき、各節点に作用する力Ｆｅを要素毎に算出すること（ＳＴ１０５）、
（ｆ）時点（ｔ）における前記各節点の位置と速度と前記各節点に作用する力と前記集中質量行列とに基づき分子動力学ソルバに次の時点（ｔ＋単位時間）における各節点の位置及び速度を算出させること（ＳＴ１０６）、
を含み、
初期時点から目標の時点に到達するまで、（ｄ）（ｅ）及び（ｆ）の処理を、時点（ｔ）を単位時間経過させつつ繰り返し実行する。
【００３７】
第１実施形態の構造物のＦＥＭ解析システムは、
構造物を複数の要素ｅで表現したＦＥＭモデルデータＤ１であって、各要素ｅを区画する複数の節点ｉおよび物性値（ヤング率Ｅ、ポアソン比ｖ、質量密度ρ）が設定されたＦＥＭモデルデータＤ１を取得するＦＥＭモデル取得部１０と、
ＦＥＭモデルデータＤ１に基づき、構造物の質量を各節点に配分して表す集中質量行列Ｍｅを取得する集中質量行列取得部１１と、
ＦＥＭモデルデータＤ１に基づき、各節点の初期位置及び初期速度を取得する初期位置速度取得部１２と、
各節点の初期位置及び物性値に基づき要素剛性行列Ｋｅを要素毎に算出する要素剛性行列算出部１３と、
各節点の初期位置に対する時点（ｔ）の位置の変位を表す変位行列ｄｅを要素毎に算出する変位行列算出部１４と、
要素剛性行列Ｋｅと変位行列ｄｅとに基づき、各節点に作用する力Ｆｅを要素毎に算出する力算出部１５と、
時点（ｔ）における前記各節点の位置と速度と前記各節点に作用する力と前記集中質量行列とに基づき次の時点（ｔ＋単位時間）における各節点の位置及び速度を算出する分子動力学ソルバ１６と、
を備え、
初期時点から目標の時点に到達するまで、変位行列算出部１４、力算出部１５及び分子動力学ソルバ１６による処理を、時点（ｔ）を単位時間経過させつつ繰り返し実行するように構成されている。
【００３８】
このように、ＦＥＭモデルにおける節点を１粒子として取り扱い、ＦＥＭモデルに基づく要素剛性行列と初期位置に対する位置変位とに基づき節点に作用する力を算出するので、分子動力学ソルバを用いて時点（ｔ）における力と節点位置から次の時点（ｔ＋単位時間）の各節点の位置が算出でき、ＦＥＭの計算を、分子動力学ソルバを用いて実行可能となる。それでいて、分子動力学ソルバは無償利用可能なソフトを含めて並列計算機能を有するので、構造物のＦＥＭ解析をリーズナブルに並列計算で実現可能となる。
【００３９】
第１実施形態のように、前記（ｃ）における要素剛性行列算出部１３による要素剛性行列Ｋｅの算出は、単位時間経過する毎に算出されることが好ましい。
【００４０】
このようにすれば、要素剛性行列Ｋｅを単位時間経過する毎に毎回算出しても計算コストが大きくなく、全ての要素の要素剛性行列Ｋｅをメモリに記憶することによって生じるメモリの消費を抑制することが可能となる。
【００４１】
第１実施形態のように、要素ｅの要素剛性行列Ｋｅは、式（１）〜（３）で表現されることが好ましい。好適な実施形態である。
【００４２】
第１実施形態のように、要素ｅを構成する各節点ｉ（ｉ＝１，…，Ｌｅ）に作用する力（Ｆｅｉ，ｘ、Ｆｅｉ，ｙ、Ｆｅｉ，ｚ）は、式（４）〜（５）で表現されることが好ましい。好適な実施形態である。
【００４３】
第１実施形態に係るプログラムは、上記方法をコンピュータに実行させるプログラムである。このプログラムを実行することによっても、上記方法の奏する作用効果を得ることが可能となる。
【００４４】
以上、本発明の実施形態について図面に基づいて説明したが、具体的な構成は、これらの実施形態に限定されるものでないと考えられるべきである。本発明の範囲は、上記した実施形態の説明だけではなく特許請求の範囲によって示され、さらに特許請求の範囲と均等の意味および範囲内でのすべての変更が含まれる。
【００４５】
例えば、図１に示す各部１０〜１５は、所定プログラムをコンピュータのプロセッサで実行することで実現しているが、各部を専用回路で構成してもよい。また、本実施形態では１つのコンピュータにおけるプロセッサが各部１０〜１５を実装しているが、少なくとも１又は複数のプロセッサに分散して実装してもよい。
【００４６】
＜第２実施形態＞
第２実施形態について、図６〜図８を用いて説明する。第１実施形態と同じ箇所は同じ符号を付けて説明を省略する。
【００４７】
［第２実施形態：構造物のＦＥＭ解析システム］
図６に示すように、第２実施形態のシステム１０１は、第１実施形態のシステム１と同様に、ＦＥＭモデル取得部１０、集中質量行列取得部１１、初期位置速度取得部１２、要素剛性行列算出部１３、及び、分子動力学ソルバ１６を有する。第２実施形態のシステム１０１は、第１実施形態のシステム１が有する変位行列算出部１４及び力算出部１５の代わりに、変形行列算出部１１４ａ、右極分解部１１４ｂ、変位ベクトル算出部１１４ｃ、力算出部１１５ａ及び力方向変換部１１５ｂを有する。
【００４８】
第１実施形態のＦＥＭモデル取得部１０が取得するＦＥＭモデルの要素の形状は、四面体要素に限定されず、六面体などの任意の形状でよい。一方、第２実施形態のＦＥＭモデル取得部１０が取得するＦＥＭモデルの要素の形状は、四面体要素に限定される。
【００４９】
図６に示す変形行列算出部１１４ａは、各節点の初期位置に対する時点（ｔ）の位置への変形を示す変形行列Ａを算出する。変形行列Ａは、初期位置における第１節点を基準とした第２〜４節点の相対位置ベクトル（ｒ１２、ｒ１３、ｒ１４）と、時点（ｔ）における第１節点を基準とした第２〜４節点の相対位置ベクトル（ｒ’１２、ｒ’１３、ｒ’１４）とに基づき算出される。
【００５０】
例えば、図８に示すように、初期位置にあった要素ｅ＝１が、時刻（ｔ）における位置へ変形したとする。初期位置における要素ｅ＝１の第１節点（ｉ＝１）、第２節点（ｉ＝２）、第３節点（ｉ＝３）及び第４節点（ｉ＝４）のｘｙｚ座標が、それぞれ（ｘ１、ｙ１、ｚ１）、（ｘ２、ｙ２、ｚ２）、（ｘ３、ｙ３、ｚ３）及び（ｘ４、ｙ４、ｚ４）とする。同様に、時刻（ｔ）における第１〜４節点のｘｙｚ座標が、（ｘ’１、ｙ’１、ｚ’１）、（ｘ’２、ｙ’２、ｚ’２）、（ｘ’３、ｙ’３、ｚ’３）及び（ｘ’４、ｙ’４、ｚ’４）とする。この場合、初期位置における第１節点を基準とした第２〜４節点の相対位置ベクトル（ｒ１２、ｒ１３、ｒ１４）は式（７）で表現される。時点（ｔ）における第１節点を基準とした第２〜４節点の相対位置ベクトル（ｒ’１２、ｒ’１３、ｒ’１４）は式（８）で表現される。
【数５】
2020177388.tif 000007
【００５１】
下式（９）に示すように、初期位置における相対位置ベクトル（ｒ１２、ｒ１３、ｒ１４）に対応する行列に変形行列Ａをかければ、変形後の相対位置ベクトル（ｒ’１２、ｒ’１３、ｒ’１４）に対応する行列となる。よって、変形行列Ａは、式（９）で表現される。
【数６】
2020177388.tif 000008
【００５２】
上記では、ＦＥＭモデルの要素が四面体に限定されると述べた理由は、４面体であれば、相対位置ベクトルが３つであり、変形行列Ａを一意に求めることができる。しかし、４面体以外の例えば六面体であれば節点が８個となり、相対位置ベクトルが７つとなり、変形行列Ａを一意に求めることができないからである。
【００５３】
図６に示す右極分解部１１４ｂは、右極分解により変形行列Ａから、歪による変形を表す歪変形行列Ｑと、回転による変形を示す回転変形行列Ｒと、を算出する。変形行列Ａと、歪変形行列Ｑと、回転変形行列Ｒの関係は、式（１０）で表現される。式（１０）におけるＩは単位行列を示す。
Ａ＝ＲＱ、  ＲＲＴ＝Ｉ、  Ｑ＝ＱＴ      …（１０）
【００５４】
図６に示す変位ベクトル算出部１１４ｃは、歪変形行列Ｑを用いて、回転を伴わない状態での歪に起因する、各節点の初期位置に対する変位を表す変位ベクトル（ｓｅ、ｕｅ、ｗｅ）を要素毎に算出する。変位ベクトル（ｓｅ、ｕｅ、ｗｅ）は、式（１１）で表現される。
【数７】
2020177388.tif 000009
ただし、変位ベクトルｓｅは要素ｅのｘ成分の変位を示し、変位ベクトルｕｅは要素ｅのｙ成分の変位を示し、変位ベクトルｗｅは要素ｅのｚ成分の変位を示す。Ｓｅｉは要素ｅの節点ｉにおけるｘ成分の変位を示し、Ｕｅｉは要素ｅの節点ｉにおけるｙ成分の変位を示し、Ｗｅｉは要素ｅの節点ｉにおけるｚ成分の変位を示す。
【００５５】
図６に示す力算出部１１５ａは、要素剛性行列算出部１３が算出した要素剛性行列Ｋｅと、変位ベクトル算出部１１４ｃが算出した変位ベクトル（ｓｅ、ｕｅ、ｗｅ）と、に基づき、歪のみに起因して要素ｅの各節点ｉ（ｉ＝１〜４）に作用する力ｆｅ（ｆｅｉ，ｘ、ｆｅｉ，ｙ、ｆｅｉ，ｚ）を要素毎に算出する。歪のみに起因して要素ｅを構成する各節点ｉに作用する力ｆｅは、式（１２）で表される。
【数８】
2020177388.tif 000010
ただし、ｆｅｉ，ｘは要素ｅの節点ｉに作用する力のｘ成分を示し、ｆｅｉ，ｙは要素ｅの節点ｉに作用する力のｙ成分を示し、ｆｅｉ，ｚは要素ｅの節点ｉに作用する力のｚ成分を示す。
【００５６】
図６に示す力方向変換部１１５ｂは、回転変形行列Ｒを用いて歪のみに起因して各節点に作用する力（ｆｅｉ，ｘ、ｆｅｉ，ｙ、ｆｅｉ，ｚ）の向きを変換することで、回転及び歪を伴った各節点に作用する力Ｆｅを算出する。力Ｆｅは、式（１３）で表現される。この処理は、力算出部１１５ａが算出した力の向きを、回転変形行列Ｒを用いて変更する処理である。
【数９】
2020177388.tif 000011
ただし、Ｆｅｉ，ｘは要素ｅの節点ｉに作用する回転及び歪を伴った力のｘ成分を示し、Ｆｅｉ，ｙは要素ｅの節点ｉに作用する回転及び歪を伴った力のｙ成分を示し、Ｆｅｉ，ｚは要素ｅの節点ｉに作用する回転及び歪を伴った力のｚ成分を示す。
【００５７】
［第２実施形態：構造物のＦＥＭ解析方法］
図６に示すシステム１０１における１又は複数のプロセッサが実行する、構造物のＦＥＭ解析方法について、図７を用いて説明する。
【００５８】
図７に示すように、第１実施形態のステップＳＴ１０１、ＳＴ１０２、ＳＴ１０６、ＳＴ１０７は第２実施形態でも同じである。第２実施形態では、第１実施形態におけるステップＳＴ１０３、ＳＴ１０４、ＳＴ１０５の代わりに、ＳＴ２０３、ＳＴ２０４、ＳＴ２０５、ＳＴ２０６、ＳＴ２０７を実行する。
【００５９】
ステップＳＴ２０３において、変形行列算出部１１４ａが、初期位置における第１節点を基準とした第２〜４節点の相対位置ベクトル（ｒ１２、ｒ１３、ｒ１４）と、時点（ｔ）における第１節点を基準とした第２〜４節点の相対位置ベクトル（ｒ’１２、ｒ’１３、ｒ’１４）と、に基づき各節点の変形を示す変形行列Ａを算出する。
【００６０】
次のステップＳＴ２０４において、右極分解部１１４ｂが、右極分解により変形行列Ａから、歪による変形を表す歪変形行列Ｑと、回転による変形を示す回転変形行列Ｒと、を算出する。
【００６１】
次のステップＳＴ２０５において、変位ベクトル算出部１１４ｃが、歪変形行列Ｑを用いて、回転を伴わない状態での歪に起因する、各節点の初期位置に対する変位を表す変位ベクトル（ｓｅ、ｕｅ、ｗｅ）を要素毎に算出する。
【００６２】
次のステップＳＴ２０６において、力算出部１１５ａが、要素剛性行列Ｋｅと変位ベクトル（ｓｅ、ｕｅ、ｗｅ）とに基づき、歪のみに起因して各節点に作用する力ｆｅ（ｆｅｉ，ｘ、ｆｅｉ，ｙ、ｆｅｉ，ｚ）を要素毎に算出する。
【００６３】
次のステップＳＴ２０７において、力方向変換部１１５ｂが、回転変形行列Ｒを用いて歪のみに起因して各節点に作用する力ｆｅ（ｆｅｉ，ｘ、ｆｅｉ，ｙ、ｆｅｉ，ｚ）の向きを変換することで、回転及び歪を伴った各節点に作用する力Ｆｅを算出する。
【００６４】
第２実施形態のシミュレーション方法の結果を説明する。
図９は、立方体の初期位置を示す側面図である。図１０は、ねじり変形させている様子を示す側面図である。図９及び図１０に示すように、３０ｍｍ×１０ｍｍ×１０ｍｍの直方体のねじり変形を模擬した。直方体の物性値は、ヤング率Ｅが１ＭＰａ、ポアソン比ｖが０．４９、質量密度ρは９７０ｋｇ／ｍ３であり、ゴムに相当する構造物である。要素タイプは、１要素が４つの節点を有する四面体１次要素である。ｙｚ底面（図９〜１０の左端面）を固定し、ｙｚ上面（図９〜１０の右端面）にｘ軸回りに１０００Ｎｍのトルクを作用させた。
【００６５】
図１１は、第１実施形態のシステム１を用いた結果であり、直方体のｙｚ上面を示す平面図である。図１１に示すように、ｙｚ上面の面積が１５２．８ｍｍ２であり、回転運動が圧縮変形と誤認され、ｙｚ上面が拡張しており、回転による幾何学的非線形に対応していない。
【００６６】
これに対し、図１２は、第２実施形態のシステム１０１を用いた結果であり、直方体のｙｚ上面を示す平面図である。図１２に示すように、ｙｚ上面の面積が１００．２ｍｍ２であり、回転による幾何学的非線形に対応させることができていることがわかる。
【００６７】
以上のように、第２実施形態の構造物のＦＥＭ解析方法は、
１又は複数のプロセッサが実行する方法であって、
（ａ）構造物を複数の四面体要素ｅで表現したＦＥＭモデルデータＤ１であって、各要素ｅを区画する４つの節点ｉ（ｉ＝１〜４）および物性値（ヤング率Ｅ、ポアソン比ｖ、質量密度ρ）が設定されたＦＥＭモデルデータＤ１を取得すること（ＳＴ１０１）、
（ｂ）ＦＥＭモデルデータＤ１に基づき、構造物の質量を各節点に配分して表す集中質量行列Ｍｅと、各節点の初期位置及び初期速度と、を取得すること（ＳＴ１０２）、
（ｃ）各節点の初期位置及び物性値に基づき要素剛性行列Ｋｅを要素毎に算出すること（ＳＴ１０３）、
（ｄ１）初期位置における第１節点を基準とした第２〜４節点の相対位置ベクトルと、時点（ｔ）における第１節点を基準とした第２〜４節点の相対位置ベクトルと、に基づき各節点の変形を示す変形行列Ａを算出すること（ＳＴ２０３）、
（ｄ２）右極分解により前記変形行列Ａから、歪による変形を表す歪変形行列Ｑと、回転による変形を示す回転変形行列Ｒと、を算出すること（ＳＴ２０４）、
（ｄ３）前記歪変形行列Ｑを用いて、回転を伴わない状態での歪に起因する、各節点の初期位置に対する変位を表す変位ベクトルを要素毎に算出すること（ＳＴ２０５）、
（ｅ１）前記要素剛性行列と前記変位ベクトルとに基づき、歪のみに起因して各節点に作用する力を要素毎に算出すること（ＳＴ２０６）、
（ｅ２）前記回転変形行列Ｒを用いて前記歪のみに起因して各節点に作用する力の向きを変換することで、回転及び歪を伴った各節点に作用する力を算出すること（ＳＴ２０７）、
（ｆ）時点（ｔ）における前記各節点の位置と速度と前記回転及び歪を伴った各節点に作用する力と前記集中質量行列とに基づき分子動力学ソルバに次の時点（ｔ＋単位時間）における各節点の位置及び速度を算出させること（ＳＴ１０６）、
を含み、
初期時点から目標の時点に到達するまで、（ｄ１）（ｄ２）（ｄ３）（ｅ１）（ｅ２）及び（ｆ）の処理を、時点（ｔ）を単位時間経過させつつ繰り返し実行する。
【００６８】
第２実施形態の構造物のＦＥＭ解析システムは、
構造物を複数の四面体要素ｅで表現したＦＥＭモデルデータＤ１であって、各要素ｅを区画する４つの節点ｉ（ｉ＝１〜４）および物性値（ヤング率Ｅ、ポアソン比ｖ、質量密度ρ）が設定されたＦＥＭモデルデータＤ１を取得するＦＥＭモデル取得部１０と、
ＦＥＭモデルデータＤ１に基づき、構造物の質量を各節点に配分して表す集中質量行列Ｍｅを取得する集中質量行列取得部１１と、
ＦＥＭモデルデータＤ１に基づき、各節点の初期位置及び初期速度を取得する初期位置速度取得部１２と、
各節点の初期位置及び物性値に基づき要素剛性行列Ｋｅを要素毎に算出する要素剛性行列算出部１３と、
初期位置における第１節点を基準とした第２〜４節点の相対位置ベクトルと、時点（ｔ）における第１節点を基準とした第２〜４節点の相対位置ベクトルと、に基づき各節点の変形を示す変形行列Ａを算出する変形行列算出部１１４ａと、
右極分解により前記変形行列Ａから、歪による変形を表す歪変形行列Ｑと、回転による変を示す回転変形行列Ｒと、を算出する右極分解部１１４ｂと、
前記歪変形行列Ｑを用いて、回転を伴わない状態での歪に起因する、各節点の初期位置に対する変位を表す変位ベクトルを要素毎に算出する変位ベクトル算出部１１４ｃと、
前記要素剛性行列と前記変位ベクトルとに基づき、歪のみに起因して各節点に作用する力を要素毎に算出する力算出部１１５ａと、
前記回転変形行列Ｒを用いて前記歪のみに起因して各節点に作用する力の向きを変換することで、回転及び歪を伴った各節点に作用する力を算出する力方向変換部１１５ｂと、
時点（ｔ）における前記各節点の位置と速度と前記回転及び歪を伴った各節点に作用する力と前記集中質量行列とに基づき分子動力学ソルバに次の時点（ｔ＋単位時間）における各節点の位置及び速度を算出する分子動力学ソルバ１６と、
を備え、
初期時点から目標の時点に到達するまで、変形行列算出部１１４ａ、右極分解部１１４ｂ、変位ベクトル算出部１１４ｃ、力算出部１１５ａ、力方向変換部１１５ｂ及び前記分子動力学ソルバによる処理を、時点（ｔ）を単位時間経過させつつ繰り返し実行するように構成されている。
【００６９】
このように、ＦＥＭモデルにおける節点を１粒子として取り扱い、ＦＥＭモデルに基づく要素剛性行列と各節点に作用する力を算出するので、分子動力学ソルバを用いて時点（ｔ）における力と節点位置から次の時点（ｔ＋単位時間）の各節点の位置が算出でき、ＦＥＭの計算を、分子動力学ソルバを用いて実行可能となる。それでいて、分子動力学ソルバは無償利用可能なソフトを含めて並列計算機能を有するので、構造物のＦＥＭ解析をリーズナブルに並列計算で実現可能となる。
さらに、各節点の初期位置から時点（ｔ）の位置に基づく変形行列Ａから、歪変形行列Ｑと回転変形行列Ｒを求め、歪変形行列Ｑを用いて歪のみに起因する力を算出し、その後、力の向きを回転変形行列Ｒで変換しているので、幾何学的非線形に対応可能となる。
【００７０】
第２実施形態のように、前記（ｃ）における要素剛性行列算出部１３による要素剛性行列Ｋｅの算出は、単位時間経過する毎に算出されることが好ましい。
【００７１】
このようにすれば、要素剛性行列Ｋｅを単位時間経過する毎に毎回算出しても計算コストが大きくなく、全ての要素の要素剛性行列Ｋｅをメモリに記憶することによって生じるメモリの消費を抑制することが可能となる。
【００７２】
第２実施形態のように、要素ｅの要素剛性行列Ｋｅは、式（１）〜（３）で表現されることが好ましい。好適な実施形態である。
【００７３】
第２実施形態のように、変形行列Ａは、式（９）で表現されることが好ましい。好適な実施形態である。
【００７４】
第２実施形態のように、回転を伴わない状態での歪に起因する、各節点の初期位置に対する変位を表す変位ベクトル（ｓｅ、ｕｅ、ｗｅ）は、式（１１）で表現されることが好ましい。好適な実施形態である。
【００７５】
第２実施形態のように、歪のみに起因して各節点に作用する力ｆｅは、式（１２）で表現されることが好ましい。好適な実施形態である。
【００７６】
第２実施形態のように、回転及び歪を伴った各節点に作用する力Ｆｅは、式（１３）で表現されることが好ましい。好適な実施形態である。
【００７７】
第２実施形態に係るプログラムは、上記方法をコンピュータに実行させるプログラムである。このプログラムを実行することによっても、上記方法の奏する作用効果を得ることが可能となる。
【００７８】
上記の各実施形態で採用している構造を他の任意の実施形態に採用することは可能である。各部の具体的な構成は、上述した実施形態のみに限定されるものではなく、本発明の趣旨を逸脱しない範囲で種々変形が可能である。
【符号の説明】
【００７９】
１０      ＦＥＭモデル取得部
１１      集中質量行列取得部
１２      初期位置速度取得部
１３      要素剛性行列算出部
１１４ａ  変形行列算出部
１１４ｂ  右極分解部
１１４ｃ  変位ベクトル算出部
１１５ａ  力算出部
１１５ｂ  力方向変換部
１６      分子動力学ソルバ
【図１】
2020177388.tif 000022
【図２】
2020177388.tif 000023
【図３】
2020177388.tif 000024
【図４】
2020177388.tif 000025
【図５】
2020177388.tif 000026
【図６】
2020177388.tif 000027
【図７】
2020177388.tif 000028
【図８】
2020177388.tif 000029
【図９】
2020177388.tif 000030
【図１０】
2020177388.tif 000031
【図１１】
2020177388.tif 000032
【図１２】
2020177388.tif 000033