(19)【発行国】日本国特許庁(JP)
(12)【公報種別】特許公報(B2)
(11)【特許番号】
(24)【登録日】2023-02-27
(45)【発行日】2023-03-07
(54)【発明の名称】構造物のFEM解析方法、システム及びプログラム
(51)【国際特許分類】
G06F 30/23 20200101AFI20230228BHJP
【FI】
G06F30/23
(21)【出願番号】P 2019078393
(22)【出願日】2019-04-17
【審査請求日】2022-02-15
(73)【特許権者】
【識別番号】000003148
【氏名又は名称】TOYO TIRE株式会社
(74)【代理人】
【識別番号】110000729
【氏名又は名称】弁理士法人ユニアス国際特許事務所
(72)【発明者】
【氏名】日野 理
【審査官】松浦 功
(56)【参考文献】
【文献】特開2017-049971(JP,A)
【文献】特開2016-081183(JP,A)
【文献】田川和義 外2名,オンラインリメッシュ型共回転系変形モデルによる実時間非線形変形計算と力覚操作への応用,電子情報通信学会技術研究報告,一般社団法人電子情報通信学会,2014年03月10日,Vol. 113, No. 501 ,pp. 117-122
【文献】藤田智 外2名,分子動力学・有限要素弾性解析の動的結合技術の開発,日本金属学会講演概要 2006年秋期(第139回)大会,社団法人日本金属学会,2006年09月16日,p. 447
(58)【調査した分野】(Int.Cl.,DB名)
G06F 30/00 -G06F 30/398
G16B 5/00 -G16C 99/00
G16Z 99/00
Google Scholar
(57)【特許請求の範囲】
【請求項1】
1又は複数のプロセッサが実行する方法であって、
(a)構造物を複数の四面体要素で表現したFEMモデルデータであって、各要素を区画する4つの節点および物性値が設定されたFEMモデルデータを取得すること、
(b)前記FEMモデルデータに基づき、前記構造物の質量を各節点に配分して表す集中質量行列と、各節点の初期位置及び初期速度と、を取得すること、
(c)前記各節点の初期位置及び前記物性値に基づき要素剛性行列を要素毎に算出すること、
(d1)初期位置における第1節点を基準とした第2~4節点の相対位置ベクトルと、時点(t)における第1節点を基準とした第2~4節点の相対位置ベクトルと、に基づき各節点の変形を示す変形行列Aを算出すること、
(d2)右極分解により前記変形行列Aから、歪による変形を表す歪変形行列Qと、回転による変形を示す回転変形行列Rと、を算出すること、
(d3)前記歪変形行列Qを用いて、回転を伴わない状態での歪に起因する、各節点の初期位置に対する変位を表す変位ベクトルを要素毎に算出すること、
(e1)前記要素剛性行列と前記変位ベクトルとに基づき、歪のみに起因して各節点に作用する力を要素毎に算出すること、
(e2)前記回転変形行列Rを用いて前記歪のみに起因して各節点に作用する力の向きを変換することで、回転及び歪を伴った各節点に作用する力を算出すること、
(f)時点(t)における前記各節点の位置と速度と前記回転及び歪を伴った各節点に作用する力と前記集中質量行列とに基づき分子動力学ソルバに次の時点(t+単位時間)における各節点の位置及び速度を算出させること、
を含み、
初期時点から目標の時点に到達するまで、前記(d1)(d2)(d3)(e1)(e2)及び(f)の処理を、前記時点(t)を単位時間経過させつつ繰り返し実行する、構造物のFEM解析方法。
【請求項2】
前記(c)における前記要素剛性行列の算出は、前記単位時間経過する毎に算出される、請求項1に記載の方法。
【請求項3】
要素eの要素剛性行列K
eは、式(1)~(3)で表現される、請求項1又は2に記載の方法。
【数1】
Veは要素eの体積を示し、Leは要素eにおける節点数を示し、N
ei(r)は要素eの節点i(i=1,…,Le)の形状関数を示し、Eは構造物のヤング率を示し、vは構造物のポアソン比を示す。
【請求項4】
前記変形行列Aは、式(9)で表現される、請求項1~3のいずれかに記載の方法。
【数2】
ただし、初期位置における第1~4節点のxyz座標が、それぞれ(x
1、y
1、z
1)、(x
2、y
2、z
2)、(x
3、y
3、z
3)及び(x
4、y
4、z
4)であり、時刻(t)における第1~4節点のxyz座標が、それぞれ(x’
1、y’
1、z’
1)、(x’
2、y’
2、z’
2)、(x’
3、y’
3、z’
3)及び(x’
4、y’
4、z’
4)である。
【請求項5】
前記回転を伴わない状態での歪に起因する、各節点の初期位置に対する変位を表す変位ベクトル(s
e、u
e、w
e)は、式(11)で表現される、請求項1~4のいずれかに記載の方法。
【数3】
ただし、変位ベクトルs
eは要素eのx成分の変位を示し、変位ベクトルu
eは要素eのy成分の変位を示し、変位ベクトルw
eは要素eのz成分の変位を示す。S
e
iは要素eの節点iにおけるx成分の変位を示し、U
e
iは要素eの節点iにおけるy成分の変位を示し、W
e
iは要素eの節点iにおけるz成分の変位を示し、初期位置における第1~4節点のxyz座標が、それぞれ(x
1、y
1、z
1)、(x
2、y
2、z
2)、(x
3、y
3、z
3)及び(x
4、y
4、z
4)である。
【請求項6】
前記歪のみに起因して各節点に作用する力f
eは、式(12)で表現される、請求項5に記載の方法。
【数4】
ただし、f
e
i,xは要素eの節点iに作用する力のx成分を示し、f
e
i,yは要素eの節点iに作用する力のy成分を示し、f
e
i,zは要素eの節点iに作用する力のz成分を示す。
【請求項7】
前記回転及び歪を伴った各節点に作用する力F
eは、式(13)で表現される、請求項6に記載の方法。
【数5】
ただし、F
e
i,xは要素eの節点iに作用する回転及び歪を伴った力のx成分を示し、F
e
i,yは要素eの節点iに作用する回転及び歪を伴った力のy成分を示し、F
e
i,zは要素eの節点iに作用する回転及び歪を伴った力のz成分を示す。
【請求項8】
構造物を複数の四面体要素で表現したFEMモデルデータであって、各要素を区画する4つの節点および物性値が設定されたFEMモデルデータを取得するFEMモデル取得部と、
前記FEMモデルデータに基づき、前記構造物の質量を各節点に配分して表す集中質量行列を取得する集中質量行列取得部と、
前記FEMモデルデータに基づき、各節点の初期位置及び初期速度を取得する初期位置速度取得部と、
前記各節点の初期位置及び前記物性値に基づき要素剛性行列を要素毎に算出する要素剛性行列算出部と、
初期位置における第1節点を基準とした第2~4節点の相対位置ベクトルと、時点(t)における第1節点を基準とした第2~4節点の相対位置ベクトルと、に基づき各節点の変形を示す変形行列Aを算出する変形行列算出部と、
右極分解により前記変形行列Aから、歪による変形を表す歪変形行列Qと、回転による変形を示す回転変形行列Rと、を算出する右極分解部と、
前記歪変形行列Qを用いて、回転を伴わない状態での歪に起因する、各節点の初期位置に対する変位を表す変位ベクトルを要素毎に算出する変位ベクトル算出部と、
前記要素剛性行列と前記変位ベクトルとに基づき、歪のみに起因して各節点に作用する力を要素毎に算出する力算出部と、
前記回転変形行列Rを用いて前記歪のみに起因して各節点に作用する力の向きを変換することで、回転及び歪を伴った各節点に作用する力を算出する力方向変換部と、
時点(t)における前記各節点の位置と速度と前記回転及び歪を伴った各節点に作用する力と前記集中質量行列とに基づき次の時点(t+単位時間)における各節点の位置及び速度を算出する分子動力学ソルバと、
を備え、
初期時点から目標の時点に到達するまで、前記変形行列算出部、前記右極分解部、前記変位ベクトル算出部、前記力算出部、前記力方向変換部及び前記分子動力学ソルバによる処理を、前記時点(t)を単位時間経過させつつ繰り返し実行するように構成されている、構造物のFEM解析システム。
【請求項9】
前記要素剛性行列算出部による前記要素剛性行列の算出は、前記単位時間経過する毎に算出される、請求項8に記載のシステム。
【請求項10】
要素eの要素剛性行列K
eは、式(1)~(3)で表現される、請求項8又は9に記載のシステム。
【数6】
Veは要素eの体積を示し、Leは要素eにおける節点数を示し、N
ei(r)は要素eの節点i(i=1,…,Le)の形状関数を示し、Eは構造物のヤング率を示し、vは構造物のポアソン比を示す。
【請求項11】
前記変形行列Aは、式(9)で表現される、請求項8~10のいずれかに記載の方法。
【数7】
ただし、初期位置における第1~4節点のxyz座標が、それぞれ(x
1、y
1、z
1)、(x
2、y
2、z
2)、(x
3、y
3、z
3)及び(x
4、y
4、z
4)であり、時刻(t)における第1~4節点のxyz座標が、それぞれ(x’
1、y’
1、z’
1)、(x’
2、y’
2、z’
2)、(x’
3、y’
3、z’
3)及び(x’
4、y’
4、z’
4)である。
【請求項12】
前記回転を伴わない状態での歪に起因する、各節点の初期位置に対する変位を表す変位ベクトル(s
e、u
e、w
e)は、式(11)で表現される、請求項8~11のいずれかに記載の方法。
【数8】
ただし、変位ベクトルs
eは要素eのx成分の変位を示し、変位ベクトルu
eは要素eのy成分の変位を示し、変位ベクトルw
eは要素eのz成分の変位を示す。S
e
iは要素eの節点iにおけるx成分の変位を示し、U
e
iは要素eの節点iにおけるy成分の変位を示し、W
e
iは要素eの節点iにおけるz成分の変位を示し、初期位置における第1~4節点のxyz座標が、それぞれ(x
1、y
1、z
1)、(x
2、y
2、z
2)、(x
3、y
3、z
3)及び(x
4、y
4、z
4)である。
【請求項13】
前記歪のみに起因して各節点に作用する力f
eは、式(12)で表現される、請求項12に記載の方法。
【数9】
ただし、f
e
i,xは要素eの節点iに作用する力のx成分を示し、f
e
i,yは要素eの節点iに作用する力のy成分を示し、f
e
i,zは要素eの節点iに作用する力のz成分を示す。
【請求項14】
前記回転及び歪を伴った各節点に作用する力F
eは、式(13)で表現される、請求項13に記載の方法。
【数10】
ただし、F
e
i,xは要素eの節点iに作用する回転及び歪を伴った力のx成分を示し、F
e
i,yは要素eの節点iに作用する回転及び歪を伴った力のy成分を示し、F
e
i,zは要素eの節点iに作用する回転及び歪を伴った力のz成分を示す。
【請求項15】
請求項1~7のいずれかに記載の方法を1又は複数のプロセッサに実行させるプログラム。
【発明の詳細な説明】
【技術分野】
【0001】
本発明は、構造物のFEM解析方法、システム及びプログラムに関する。
【背景技術】
【0002】
構造物の動解析には有限要素法(FEM;Finite Element Method)が用いられる。FEMに用いるFEMモデルデータは、解析対象となる構造物が複数の要素に分割されて表現され、複数の節点を結んだ閉領域を一つの要素として表現されている。FEMでは、節点毎に、構造物の物性値(ヤング率、ポアソン比、質量など)及び作用する力を用いた運動方程式を解くことにより、構造物の変形などを計算する。FEMを説明する文献として非特許文献1がある。
【0003】
なお、構造物シミュレーションではないが、分子シミュレーションを実行するための分子動力学ソルバ(ソフトウェア)は、多粒子系の運動方程式を高効率で解くための並列計算機能を有し、LAMMPSやgromacs等の無償利用可能なソフトウェアが知られている。分子動力学に関する文献として非特許文献2がある。
【先行技術文献】
【非特許文献】
【0004】
【文献】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.
【発明の概要】
【発明が解決しようとする課題】
【0005】
FEM解析において解析精度を向上させるためには、要素のサイズを小さくしたモデルを用いればよいが、要素サイズを小さくすれば、要素数が増大し、計算時間が急速に増加してしまう。FEM解析ソフトウェアの大半は、全ての要素を一体に取り扱う行列計算を行っているせいか、並列計算機能がないことが多く、並列計算機能を有するFEM解析ソフトウェアはライセンスコストが高価で実用的ではない。
【0006】
そこで、分子動力学ソルバを用いて構造物のFEM解析を実現することが考えられる。
【0007】
構造物のFEM解析では、回転を含んだ変形及び運動を対象とする動解析が多い。通常の弾性論における変位の定義を使用すると、有限な回転運動は圧縮変形に誤認され、構造物の一部に非物理的な膨張が生じる場合がある。これは、回転による幾何学的非線形問題として知られている。一般のFEMソルバには、回転による幾何学的非線形を考慮した増分形式を用いたアルゴリズムが知られている。しかし、分子動力学ソルバには、幾何学的非線形に対処した機能は当然なく、上記増分形式を用いたアルゴリズムを分子動力学ソルバに適用することができない。
【0008】
本発明の目的は、回転による幾何学的非線形に対応し且つ分子動力学ソルバを用いた構造物のFEM解析方法、システム及びプログラムを提供することである。
【課題を解決するための手段】
【0009】
本発明の構造物のFEM解析方法は、
1又は複数のプロセッサが実行する方法であって、
(a)構造物を複数の四面体要素で表現したFEMモデルデータであって、各要素を区画する4つの節点および物性値が設定されたFEMモデルデータを取得すること、
(b)前記FEMモデルデータに基づき、前記構造物の質量を各節点に配分して表す集中質量行列と、各節点の初期位置及び初期速度と、を取得すること、
(c)前記各節点の初期位置及び前記物性値に基づき要素剛性行列を要素毎に算出すること、
(d1)初期位置における第1節点を基準とした第2~4節点の相対位置ベクトルと、時点(t)における第1節点を基準とした第2~4節点の相対位置ベクトルと、に基づき各節点の変形を示す変形行列Aを算出すること、
(d2)右極分解により前記変形行列Aから、歪による変形を表す歪変形行列Qと、回転による変形を示す回転変形行列Rと、を算出すること、
(d3)前記歪変形行列Qを用いて、回転を伴わない状態での歪に起因する、各節点の初期位置に対する変位を表す変位ベクトルを要素毎に算出すること、
(e1)前記要素剛性行列と前記変位ベクトルとに基づき、歪のみに起因して各節点に作用する力を要素毎に算出すること、
(e2)前記回転変形行列Rを用いて前記歪のみに起因して各節点に作用する力の向きを変換することで、回転及び歪を伴った各節点に作用する力を算出すること、
(f)時点(t)における前記各節点の位置と速度と前記回転及び歪を伴った各節点に作用する力と前記集中質量行列とに基づき分子動力学ソルバに次の時点(t+単位時間)における各節点の位置及び速度を算出させること、
を含み、
初期時点から目標の時点に到達するまで、前記(d1)(d2)(d3)(e1)(e2)及び(f)の処理を、前記時点(t)を単位時間経過させつつ繰り返し実行する。
【0010】
このように、FEMモデルにおける節点を1粒子として取り扱い、FEMモデルに基づく要素剛性行列と各節点に作用する力を算出するので、分子動力学ソルバを用いて時点(t)における力と節点位置から次の時点(t+単位時間)の各節点の位置が算出でき、FEMの計算を、分子動力学ソルバを用いて実行可能となる。それでいて、分子動力学ソルバは無償利用可能なソフトを含めて並列計算機能を有するので、構造物のFEM解析をリーズナブルに並列計算で実現可能となる。
さらに、各節点の初期位置から時点(t)の位置に基づく変形行列Aから、歪変形行列Qと回転変形行列Rを求め、歪変形行列Qを用いて歪のみに起因する力を算出し、その後、力の向きを回転変形行列Rで変換しているので、幾何学的非線形に対応可能となる。
【図面の簡単な説明】
【0011】
【
図1】本発明の第1実施形態の構造物のFEM解析システムを示す図
【
図2】第1実施形態のシステムが実行する処理のフローチャート
【
図4】第1実施形態の実施例のFEMモデルの各節点の初期位置を示す図
【
図5】第1実施形態の実施例の壁に衝突したときのFEMモデルの各節点の位置を示す図
【
図6】第2実施形態の構造物のFEM解析システムを示す図
【
図7】第2実施形態のシステムが実行する処理のフローチャート
【
図8】初期位置にある各節点の座標と、時間(t)における各節点の座標とを示す図
【
図9】第2実施形態の実施例の立方体FEMモデルの初期位置を示す側面図
【
図10】第2実施形態の実施例の立方体FEMモデルをねじり変形させている様子を示す側面図
【
図11】
図10に示すねじり変形を第1実施形態のシステムを用いた結果であり、直方体のyz上面を示す平面図
【
図12】
図10に示すねじり変形を第2実施形態のシステムを用いた結果であり、直方体のyz上面を示す平面図
【発明を実施するための形態】
【0012】
<第1実施形態>
以下、本発明の第1実施形態を、
図1~
図5を参照して説明する。第1実施形態は、幾何学的非線形を考慮せずに、分子動力学ソルバを用いた構造物のFEM解析を実現している。第2実施形態は、幾何学的非線形に対応した分子動力学ソルバを用いた構造物のFEM解析を実現している。
【0013】
[構造物のFEM解析システム]
第1実施形態のシステム1は、構造物を有限要素法で解析するシステムである。
【0014】
図1に示すように、システム1は、FEMモデル取得部10と、集中質量行列取得部11と、初期位置速度取得部12と、要素剛性行列算出部13と、変位行列算出部14と、力算出部15と、分子動力学ソルバ16と、を有する。これら各部10~15は、プロセッサ、メモリ、各種インターフェイス等を備えたコンピュータにおいて予め記憶されている
図2に示す処理ルーチンをプロセッサが実行することによりソフトウェア及びハードウェアが協働して実現される。本実施形態では、1つの装置におけるプロセッサが各部の処理を実行しているが、これに限定されない。例えば、ネットワークを用いて分散させ、複数のプロセッサが各部の処理を実行するように構成してもよい。すなわち、1又は複数のプロセッサが処理を実行する。
【0015】
一般的にFEM解析は、FEMソルバ(ソフトウェア)を用い、FEMソルバに対してFEMモデルデータ、FEMモデルに与える荷重、境界条件などの各種データを渡し、FEMソルバを動作させる。FEMは構造物を複数の要素に分割して解析するが、要素を構成する複数の節点を有し、節点間に分子動力学でいう相互作用(力、ポテンシャル)が作用する形で節点に作用する力を演算している、と発明者は認識した。すなわち、力の算出をFEM的に実装すれば、分子動力学ソルバが利用可能であると考えた。
【0016】
そこで、本実施形態のFEM解析法ではFEMソルバを用いずに、代わりに、FEMモデルの要素を構成する1つの節点を1つの粒子として取り扱い、複数の粒子で構成される材料の分子シミュレーションを実行する分子動力学ソルバを用いている。分子動力学ソルバには、無償利用可能なソフトウェアとして、LAMMPS(Large-scale Atomic/Molecular Massively Parallel Simulator)やgromacs等があり、本実施形態ではLAMMPSを用いている。分子動力学ソルバに渡すデータは、材料を構成する複数の粒子の初期位置及び初期速度、各粒子に作用する相互作用(力、ポテンシャル)がある。相互作用は粒子間距離に応じて定まる。分子動力学ソルバは、ポテンシャル算出部と、粒子位置更新部とを有する。一般的な分子動力学ソルバの動作の流れの概要として、ポテンシャル算出部が、時点(t)における各粒子の位置及び速度に基づき時点(t)における各粒子に作用する力(ポテンシャル)を粒子間距離に応じて算出する。次に、粒子位置更新部が、時点(t)における各粒子の位置及び速度、各粒子に作用する力(ポテンシャル)に基づき、次の時点(t+1)における各粒子の位置及び速度を算出する。このように、粒子に作用する力(ポテンシャル)の算出と、時点(t)における各粒子の位置及び速度、各粒子に作用する力に基づく次の時点(t+1)における各粒子の位置及び速度の更新と、が繰り返し実行される。この繰り返し処理は、目標となる時点まで実行される。本実施形態において、分子動力学ソルバ16のポテンシャル算出部が、FEM解析ではそのまま利用できないため、代わりに力算出部15を設けている。分子動力学ソルバ16の粒子位置更新部はそのまま利用している。
以下、具体的に説明する。
【0017】
図1に示すFEMモデル取得部10は、構造物を複数の要素で表現したFEMモデルデータD1を取得する。FEMモデルデータD1は、各要素を区画する複数の節点により構造物の立体形状を示し、要素に物性値(ヤング率E、ポアソン比v、質量密度ρ)が設定されている。
図3の例示は、1要素が4つの節点で表現される四面体であり、要素e=1、2で示される2つの要素を例示している。要素e=1を構成する節点i=1~4であり、要素e=2を構成する節点i=2~5である。節点iはそれぞれxyz座標の値を有し且つどの節点と接続されているかとどの要素に属するかを示す情報を有する。FEMモデルデータD1は、一般の有限要素法シミュレーションにて多用されている。なお、
図3に示す本実施形態では、FEMモデルデータD1の要素は四面体であるが、これに限定されない。原理的には任意の多面体要素が利用可能である。
【0018】
図1に示す集中質量行列取得部11は、構造物の質量を各節点に配分して示す集中質量行列M
eを取得する。集中質量行列取得部11は、FEMモデルデータD1に集中質量行列M
eが定義されている場合は集中質量行列M
eをそのまま取得し、集中質量行列M
eが定義されていない場合には算出する。要素eの集中質量行列M
eは次の式(6)で表現される。集中質量行列取得部11は、全ての要素について集中質量行列M
eを取得又は算出する。
図3に示す実施例では、要素数が5070であるので、5070の集中質量行列M
eを取得する。集中質量行列M
eは、対角上にのみ質量を有する。標準的なFEM解析では非対角要素をもつ整合質量行列を使用することもあるが、分子動力学ソルバの大幅な変更を避けるために集中質量行列を用いている。
【数1】
ただし、ρは構造物の質量密度を示し、Veは要素eの体積を示し、Leは要素eにおける節点数を示す。
【0019】
図1に示す初期位置速度取得部12は、FEMモデルデータD1に基づき各節点の初期位置及び初期速度を取得し、これらの値をワーキングメモリD2に記憶する。各節点の初期位置は、FEMモデルデータD1における各節点の座標(xyz座標)である。各節点の初期速度は0である。
図3に示す実施例では、要素数が5070で全節点数が1204であるので、1204個の節点の座標を取得する。
【0020】
図1に示す要素剛性行列算出部13は、各節点の初期位置及び物性値に基づき要素剛性行列Keを要素毎に算出する。要素eの要素剛性行列Keは、式(1)~(3)で表現される。要素剛性行列算出部13は、全ての要素について要素剛性行列Keを算出する。
図3に示す実施例では、要素数が5070であるので、5070の要素剛性行列Keを算出する。
【数2】
V
eは要素eの体積を示し、L
eは要素eにおける節点数を示し、N
e
i(r)は要素eの節点i(i=1,…,L
e)の形状関数を示し、Eは構造物のヤング率を示し、vは構造物のポアソン比を示す。
【0021】
図1に示す変位行列算出部14は、各節点の初期位置に対する時点(t)の位置の変位を示す変位行列d
eを要素毎に算出する。変位行列d
eは式(5)で表現される。変位行列算出部14は、全ての要素について当該要素を構成する全節点の変位を示す変位行列d
eを算出する。
図3に示す実施例では、要素数が5070であるので、5070の変位行列d
eを算出する。
【数3】
ただし、S
e
iは要素eの節点iにおけるx成分の変位を示し、U
e
iは要素eの節点iにおけるy成分の変位を示し、W
e
iは要素eの節点iにおけるz成分の変位を示す。
【0022】
図1に示す力算出部15は、要素剛性行列算出部13が算出した要素剛性行列K
eと変位行列算出部14が算出した変位行列d
eとに基づき、要素eを構成する節点i(i=1,…,L
e)に作用する力(F
e
i,x、F
e
i,y、F
e
i,z)を要素毎に算出する。要素eを構成する各節点iに作用する力F
eは式(4)で表される。
【数4】
ただし、F
e
i,xは要素eの節点iに作用する力のx成分を示し、F
e
i,yは要素eの節点iに作用する力のy成分を示し、F
e
i,zは要素eの節点iに作用する力のz成分を示す。
【0023】
図1に示す分子動力学ソルバ16は、ワーキングメモリD2に記憶されている、時点(t)における全ての要素を構成する各節点の位置及び速度と、各節点に作用する力と、集中質量行列M
eとに基づいて、次の時点(t+単位時間)における各節点の位置及び速度を算出する。算出された各節点の位置及び速度はワーキングメモリD2に記憶される。
【0024】
[構造物のFEM解析方法]
図1に示すシステム1における1又は複数のプロセッサが実行する、構造物のFEM解析方法について、
図2を用いて説明する。
【0025】
まず、ステップST101において、FEMモデル取得部10は、構造物を複数の要素で表現したFEMモデルデータD1であって、各要素を区画する複数の節点および物性値(ヤング率E、ポアソン比v、質量密度ρ等)が設定されたFEMモデルデータD1を取得する。
【0026】
次のステップST102において、集中質量行列取得部11がFEMモデルデータD1に基づき、構造物の質量を各節点に配分して表す集中質量行列Meを要素毎に取得し、初期位置速度取得部12が、FEMモデルデータD1に基づき各節点の初期位置及び初期速度を取得する。
【0027】
ステップST103、ST104、ST105は、全ての要素(
図3の例では、5070個の要素)について要素毎(e=1~5070)に実行する。ステップST103、ST104、ST105の処理を全ての要素について実行すれば、ステップST106の処理へ移行する。
【0028】
ステップST103において、要素剛性行列算出部13は、要素eの各節点の初期位置及び物性値に基づき要素剛性行列を算出する。要素eの要素剛性行列Keは、式(1)~(3)で表現される。
【0029】
ステップST104において、変位行列算出部14は、時点(t)における要素eの変位行列deを算出する。変位行列deは、要素eを構成する各節点i(i=1~Le)の初期位置に対する現在位置の変位を示す。変位行列deは式(5)で表現される。
【0030】
ステップST105において、力算出部15は、要素eの要素剛性行列Keと変位行列deとに基づき、要素eの各節点に作用する力を算出する。要素eを構成する各節点i(i=1,…,Le)に作用する力(Fe
i,x、Fe
i,y、Fe
i,z)は、式(4)~(5)で表現される
【0031】
全ての要素についてステップST103、ST104、ST105の処理を実行した後に、次にステップST106において、分子動力学ソルバ16は、時点(t)における各節点の位置と速度と各節点に作用する力と集中質量行列Meとに基づき次の時点(t+単位時間)における各節点の位置及び速度を算出する。
【0032】
次のステップST107において、ST103~ST106の処理を、初期時点から目標時点に到達するまで、時点(t)を単位時間経過させつつ繰り返し実行したか否かを判定する。目標時点に到達していない場合(ST107:NO)には、時点(t)を単位時間経過させて、ST103の処理に戻る。目標時点に到達した場合(ST107:YES)には、FEM解析処理を終了する。
【0033】
なお、本実施形態では、ステップST103における要素剛性行列Keの算出は、単位時間経過する毎に実行されている。これは、要素剛性行列Keは要素数5070あるため、メモリに記憶するとメモリ消費量が無視できないからであり、単位時間経過するたびに計算しなおしても、計算コストが低いためである。勿論、全ての要素についての要素剛性行列Keを算出した時点でメモリに記憶しておき、単位時間経過するたびに算出しないようにしてもよい。
【0034】
第1実施形態のシミュレーション方法の結果を説明する。
図4及び
図5は、半径10mmの弾性球の落下、反発運動をシミュレーションした結果である。
図4は、初期位置のFEMモデルの各節点位置を示し、
図5は、壁との衝突時におけるFEMモデルの各節点の位置を示す。物性値は、ヤング率Eが1MPa、ポアソン比vが0.49、質量密度ρは970kg/m
3であり、ゴムボールに相当する構造物である。要素タイプは、1要素が4つの節点を有する四面体1次要素であり、全体として5070要素、1204節点ある。弾性球の中心の高さを20mmとして重力により自然落下させ、高さ0mmに剛体壁を設定している。時間ステップ数は30000(約1.0秒に相当)である。
【0035】
並列計算機能が実現できたことを示すために、プロセス数を1、2、4、8、16と設定し、計算時間を次に示す。下記の通り、プロセス数に応じて計算時間が短縮されているので、FEM解析に並列計算機能をリーズナブルに実装できていることがわかる。
プロセス数:1 計算時間(秒):1106.8
プロセス数:2 計算時間(秒):560.4
プロセス数:4 計算時間(秒):282.6
プロセス数:8 計算時間(秒):150.2
プロセス数:16 計算時間(秒):83.0
【0036】
以上のように、第1実施形態の構造物のFEM解析方法は、
1又は複数のプロセッサが実行する方法であって、
(a)構造物を複数の要素eで表現したFEMモデルデータD1であって、各要素eを区画する複数の節点iおよび物性値(ヤング率E、ポアソン比v、質量密度ρ)が設定されたFEMモデルデータD1を取得すること(ST101)、
(b)FEMモデルデータD1に基づき、構造物の質量を各節点に配分して表す集中質量行列Meと、各節点の初期位置及び初期速度と、を取得すること(ST102)、
(c)各節点の初期位置及び物性値に基づき要素剛性行列Keを要素毎に算出すること(ST103)、
(d)各節点の初期位置に対する時点(t)の位置の変位を表す変位行列deを要素毎に算出すること(ST104)、
(e)要素剛性行列Keと変位行列deとに基づき、各節点に作用する力Feを要素毎に算出すること(ST105)、
(f)時点(t)における前記各節点の位置と速度と前記各節点に作用する力と前記集中質量行列とに基づき分子動力学ソルバに次の時点(t+単位時間)における各節点の位置及び速度を算出させること(ST106)、
を含み、
初期時点から目標の時点に到達するまで、(d)(e)及び(f)の処理を、時点(t)を単位時間経過させつつ繰り返し実行する。
【0037】
第1実施形態の構造物のFEM解析システムは、
構造物を複数の要素eで表現したFEMモデルデータD1であって、各要素eを区画する複数の節点iおよび物性値(ヤング率E、ポアソン比v、質量密度ρ)が設定されたFEMモデルデータD1を取得するFEMモデル取得部10と、
FEMモデルデータD1に基づき、構造物の質量を各節点に配分して表す集中質量行列Meを取得する集中質量行列取得部11と、
FEMモデルデータD1に基づき、各節点の初期位置及び初期速度を取得する初期位置速度取得部12と、
各節点の初期位置及び物性値に基づき要素剛性行列Keを要素毎に算出する要素剛性行列算出部13と、
各節点の初期位置に対する時点(t)の位置の変位を表す変位行列deを要素毎に算出する変位行列算出部14と、
要素剛性行列Keと変位行列deとに基づき、各節点に作用する力Feを要素毎に算出する力算出部15と、
時点(t)における前記各節点の位置と速度と前記各節点に作用する力と前記集中質量行列とに基づき次の時点(t+単位時間)における各節点の位置及び速度を算出する分子動力学ソルバ16と、
を備え、
初期時点から目標の時点に到達するまで、変位行列算出部14、力算出部15及び分子動力学ソルバ16による処理を、時点(t)を単位時間経過させつつ繰り返し実行するように構成されている。
【0038】
このように、FEMモデルにおける節点を1粒子として取り扱い、FEMモデルに基づく要素剛性行列と初期位置に対する位置変位とに基づき節点に作用する力を算出するので、分子動力学ソルバを用いて時点(t)における力と節点位置から次の時点(t+単位時間)の各節点の位置が算出でき、FEMの計算を、分子動力学ソルバを用いて実行可能となる。それでいて、分子動力学ソルバは無償利用可能なソフトを含めて並列計算機能を有するので、構造物のFEM解析をリーズナブルに並列計算で実現可能となる。
【0039】
第1実施形態のように、前記(c)における要素剛性行列算出部13による要素剛性行列Keの算出は、単位時間経過する毎に算出されることが好ましい。
【0040】
このようにすれば、要素剛性行列Keを単位時間経過する毎に毎回算出しても計算コストが大きくなく、全ての要素の要素剛性行列Keをメモリに記憶することによって生じるメモリの消費を抑制することが可能となる。
【0041】
第1実施形態のように、要素eの要素剛性行列Keは、式(1)~(3)で表現されることが好ましい。好適な実施形態である。
【0042】
第1実施形態のように、要素eを構成する各節点i(i=1,…,Le)に作用する力(Fe
i,x、Fe
i,y、Fe
i,z)は、式(4)~(5)で表現されることが好ましい。好適な実施形態である。
【0043】
第1実施形態に係るプログラムは、上記方法をコンピュータに実行させるプログラムである。このプログラムを実行することによっても、上記方法の奏する作用効果を得ることが可能となる。
【0044】
以上、本発明の実施形態について図面に基づいて説明したが、具体的な構成は、これらの実施形態に限定されるものでないと考えられるべきである。本発明の範囲は、上記した実施形態の説明だけではなく特許請求の範囲によって示され、さらに特許請求の範囲と均等の意味および範囲内でのすべての変更が含まれる。
【0045】
例えば、
図1に示す各部10~15は、所定プログラムをコンピュータのプロセッサで実行することで実現しているが、各部を専用回路で構成してもよい。また、本実施形態では1つのコンピュータにおけるプロセッサが各部10~15を実装しているが、少なくとも1又は複数のプロセッサに分散して実装してもよい。
【0046】
<第2実施形態>
第2実施形態について、
図6~
図8を用いて説明する。第1実施形態と同じ箇所は同じ符号を付けて説明を省略する。
【0047】
[第2実施形態:構造物のFEM解析システム]
図6に示すように、第2実施形態のシステム101は、第1実施形態のシステム1と同様に、FEMモデル取得部10、集中質量行列取得部11、初期位置速度取得部12、要素剛性行列算出部13、及び、分子動力学ソルバ16を有する。第2実施形態のシステム101は、第1実施形態のシステム1が有する変位行列算出部14及び力算出部15の代わりに、変形行列算出部114a、右極分解部114b、変位ベクトル算出部114c、力算出部115a及び力方向変換部115bを有する。
【0048】
第1実施形態のFEMモデル取得部10が取得するFEMモデルの要素の形状は、四面体要素に限定されず、六面体などの任意の形状でよい。一方、第2実施形態のFEMモデル取得部10が取得するFEMモデルの要素の形状は、四面体要素に限定される。
【0049】
図6に示す変形行列算出部114aは、各節点の初期位置に対する時点(t)の位置への変形を示す変形行列Aを算出する。変形行列Aは、初期位置における第1節点を基準とした第2~4節点の相対位置ベクトル(r
12、r
13、r
14)と、時点(t)における第1節点を基準とした第2~4節点の相対位置ベクトル(r’
12、r’
13、r’
14)とに基づき算出される。
【0050】
例えば、
図8に示すように、初期位置にあった要素e=1が、時刻(t)における位置へ変形したとする。初期位置における要素e=1の第1節点(i=1)、第2節点(i=2)、第3節点(i=3)及び第4節点(i=4)のxyz座標が、それぞれ(x
1、y
1、z
1)、(x
2、y
2、z
2)、(x
3、y
3、z
3)及び(x
4、y
4、z
4)とする。同様に、時刻(t)における第1~4節点のxyz座標が、(x’
1、y’
1、z’
1)、(x’
2、y’
2、z’
2)、(x’
3、y’
3、z’
3)及び(x’
4、y’
4、z’
4)とする。この場合、初期位置における第1節点を基準とした第2~4節点の相対位置ベクトル(r
12、r
13、r
14)は式(7)で表現される。時点(t)における第1節点を基準とした第2~4節点の相対位置ベクトル(r’
12、r’
13、r’
14)は式(8)で表現される。
【数5】
【0051】
下式(9)に示すように、初期位置における相対位置ベクトル(r
12、r
13、r
14)に対応する行列に変形行列Aをかければ、変形後の相対位置ベクトル(r’
12、r’
13、r’
14)に対応する行列となる。よって、変形行列Aは、式(9)で表現される。
【数6】
【0052】
上記では、FEMモデルの要素が四面体に限定されると述べた理由は、4面体であれば、相対位置ベクトルが3つであり、変形行列Aを一意に求めることができる。しかし、4面体以外の例えば六面体であれば節点が8個となり、相対位置ベクトルが7つとなり、変形行列Aを一意に求めることができないからである。
【0053】
図6に示す右極分解部114bは、右極分解により変形行列Aから、歪による変形を表す歪変形行列Qと、回転による変形を示す回転変形行列Rと、を算出する。変形行列Aと、歪変形行列Qと、回転変形行列Rの関係は、式(10)で表現される。式(10)におけるIは単位行列を示す。
A=RQ、 RR
T=I、 Q=Q
T …(10)
【0054】
図6に示す変位ベクトル算出部114cは、歪変形行列Qを用いて、回転を伴わない状態での歪に起因する、各節点の初期位置に対する変位を表す変位ベクトル(s
e、u
e、w
e)を要素毎に算出する。変位ベクトル(s
e、u
e、w
e)は、式(11)で表現される。
【数7】
ただし、変位ベクトルs
eは要素eのx成分の変位を示し、変位ベクトルu
eは要素eのy成分の変位を示し、変位ベクトルw
eは要素eのz成分の変位を示す。S
e
iは要素eの節点iにおけるx成分の変位を示し、U
e
iは要素eの節点iにおけるy成分の変位を示し、W
e
iは要素eの節点iにおけるz成分の変位を示す。
【0055】
図6に示す力算出部115aは、要素剛性行列算出部13が算出した要素剛性行列K
eと、変位ベクトル算出部114cが算出した変位ベクトル(s
e、u
e、w
e)と、に基づき、歪のみに起因して要素eの各節点i(i=1~4)に作用する力f
e(f
e
i,x、f
e
i,y、f
e
i,z)を要素毎に算出する。歪のみに起因して要素eを構成する各節点iに作用する力f
eは、式(12)で表される。
【数8】
ただし、f
e
i,xは要素eの節点iに作用する力のx成分を示し、f
e
i,yは要素eの節点iに作用する力のy成分を示し、f
e
i,zは要素eの節点iに作用する力のz成分を示す。
【0056】
図6に示す力方向変換部115bは、回転変形行列Rを用いて歪のみに起因して各節点に作用する力(f
e
i,x、f
e
i,y、f
e
i,z)の向きを変換することで、回転及び歪を伴った各節点に作用する力F
eを算出する。力F
eは、式(13)で表現される。この処理は、力算出部115aが算出した力の向きを、回転変形行列Rを用いて変更する処理である。
【数9】
ただし、F
e
i,xは要素eの節点iに作用する回転及び歪を伴った力のx成分を示し、F
e
i,yは要素eの節点iに作用する回転及び歪を伴った力のy成分を示し、F
e
i,zは要素eの節点iに作用する回転及び歪を伴った力のz成分を示す。
【0057】
[第2実施形態:構造物のFEM解析方法]
図6に示すシステム101における1又は複数のプロセッサが実行する、構造物のFEM解析方法について、
図7を用いて説明する。
【0058】
図7に示すように、第1実施形態のステップST101、ST102、ST106、ST107は第2実施形態でも同じである。第2実施形態では、第1実施形態におけるステップST103、ST104、ST105の代わりに、ST203、ST204、ST205、ST206、ST207を実行する。
【0059】
ステップST203において、変形行列算出部114aが、初期位置における第1節点を基準とした第2~4節点の相対位置ベクトル(r12、r13、r14)と、時点(t)における第1節点を基準とした第2~4節点の相対位置ベクトル(r’12、r’13、r’14)と、に基づき各節点の変形を示す変形行列Aを算出する。
【0060】
次のステップST204において、右極分解部114bが、右極分解により変形行列Aから、歪による変形を表す歪変形行列Qと、回転による変形を示す回転変形行列Rと、を算出する。
【0061】
次のステップST205において、変位ベクトル算出部114cが、歪変形行列Qを用いて、回転を伴わない状態での歪に起因する、各節点の初期位置に対する変位を表す変位ベクトル(se、ue、we)を要素毎に算出する。
【0062】
次のステップST206において、力算出部115aが、要素剛性行列Keと変位ベクトル(se、ue、we)とに基づき、歪のみに起因して各節点に作用する力fe(fe
i,x、fe
i,y、fe
i,z)を要素毎に算出する。
【0063】
次のステップST207において、力方向変換部115bが、回転変形行列Rを用いて歪のみに起因して各節点に作用する力fe(fe
i,x、fe
i,y、fe
i,z)の向きを変換することで、回転及び歪を伴った各節点に作用する力Feを算出する。
【0064】
第2実施形態のシミュレーション方法の結果を説明する。
図9は、立方体の初期位置を示す側面図である。
図10は、ねじり変形させている様子を示す側面図である。
図9及び
図10に示すように、30mm×10mm×10mmの直方体のねじり変形を模擬した。直方体の物性値は、ヤング率Eが1MPa、ポアソン比vが0.49、質量密度ρは970kg/m3であり、ゴムに相当する構造物である。要素タイプは、1要素が4つの節点を有する四面体1次要素である。yz底面(
図9~10の左端面)を固定し、yz上面(
図9~10の右端面)にx軸回りに1000Nmのトルクを作用させた。
【0065】
図11は、第1実施形態のシステム1を用いた結果であり、直方体のyz上面を示す平面図である。
図11に示すように、yz上面の面積が152.8mm
2であり、回転運動が圧縮変形と誤認され、yz上面が拡張しており、回転による幾何学的非線形に対応していない。
【0066】
これに対し、
図12は、第2実施形態のシステム101を用いた結果であり、直方体のyz上面を示す平面図である。
図12に示すように、yz上面の面積が100.2mm
2であり、回転による幾何学的非線形に対応させることができていることがわかる。
【0067】
以上のように、第2実施形態の構造物のFEM解析方法は、
1又は複数のプロセッサが実行する方法であって、
(a)構造物を複数の四面体要素eで表現したFEMモデルデータD1であって、各要素eを区画する4つの節点i(i=1~4)および物性値(ヤング率E、ポアソン比v、質量密度ρ)が設定されたFEMモデルデータD1を取得すること(ST101)、
(b)FEMモデルデータD1に基づき、構造物の質量を各節点に配分して表す集中質量行列Meと、各節点の初期位置及び初期速度と、を取得すること(ST102)、
(c)各節点の初期位置及び物性値に基づき要素剛性行列Keを要素毎に算出すること(ST103)、
(d1)初期位置における第1節点を基準とした第2~4節点の相対位置ベクトルと、時点(t)における第1節点を基準とした第2~4節点の相対位置ベクトルと、に基づき各節点の変形を示す変形行列Aを算出すること(ST203)、
(d2)右極分解により前記変形行列Aから、歪による変形を表す歪変形行列Qと、回転による変形を示す回転変形行列Rと、を算出すること(ST204)、
(d3)前記歪変形行列Qを用いて、回転を伴わない状態での歪に起因する、各節点の初期位置に対する変位を表す変位ベクトルを要素毎に算出すること(ST205)、
(e1)前記要素剛性行列と前記変位ベクトルとに基づき、歪のみに起因して各節点に作用する力を要素毎に算出すること(ST206)、
(e2)前記回転変形行列Rを用いて前記歪のみに起因して各節点に作用する力の向きを変換することで、回転及び歪を伴った各節点に作用する力を算出すること(ST207)、
(f)時点(t)における前記各節点の位置と速度と前記回転及び歪を伴った各節点に作用する力と前記集中質量行列とに基づき分子動力学ソルバに次の時点(t+単位時間)における各節点の位置及び速度を算出させること(ST106)、
を含み、
初期時点から目標の時点に到達するまで、(d1)(d2)(d3)(e1)(e2)及び(f)の処理を、時点(t)を単位時間経過させつつ繰り返し実行する。
【0068】
第2実施形態の構造物のFEM解析システムは、
構造物を複数の四面体要素eで表現したFEMモデルデータD1であって、各要素eを区画する4つの節点i(i=1~4)および物性値(ヤング率E、ポアソン比v、質量密度ρ)が設定されたFEMモデルデータD1を取得するFEMモデル取得部10と、
FEMモデルデータD1に基づき、構造物の質量を各節点に配分して表す集中質量行列Meを取得する集中質量行列取得部11と、
FEMモデルデータD1に基づき、各節点の初期位置及び初期速度を取得する初期位置速度取得部12と、
各節点の初期位置及び物性値に基づき要素剛性行列Keを要素毎に算出する要素剛性行列算出部13と、
初期位置における第1節点を基準とした第2~4節点の相対位置ベクトルと、時点(t)における第1節点を基準とした第2~4節点の相対位置ベクトルと、に基づき各節点の変形を示す変形行列Aを算出する変形行列算出部114aと、
右極分解により前記変形行列Aから、歪による変形を表す歪変形行列Qと、回転による変を示す回転変形行列Rと、を算出する右極分解部114bと、
前記歪変形行列Qを用いて、回転を伴わない状態での歪に起因する、各節点の初期位置に対する変位を表す変位ベクトルを要素毎に算出する変位ベクトル算出部114cと、
前記要素剛性行列と前記変位ベクトルとに基づき、歪のみに起因して各節点に作用する力を要素毎に算出する力算出部115aと、
前記回転変形行列Rを用いて前記歪のみに起因して各節点に作用する力の向きを変換することで、回転及び歪を伴った各節点に作用する力を算出する力方向変換部115bと、
時点(t)における前記各節点の位置と速度と前記回転及び歪を伴った各節点に作用する力と前記集中質量行列とに基づき分子動力学ソルバに次の時点(t+単位時間)における各節点の位置及び速度を算出する分子動力学ソルバ16と、
を備え、
初期時点から目標の時点に到達するまで、変形行列算出部114a、右極分解部114b、変位ベクトル算出部114c、力算出部115a、力方向変換部115b及び前記分子動力学ソルバによる処理を、時点(t)を単位時間経過させつつ繰り返し実行するように構成されている。
【0069】
このように、FEMモデルにおける節点を1粒子として取り扱い、FEMモデルに基づく要素剛性行列と各節点に作用する力を算出するので、分子動力学ソルバを用いて時点(t)における力と節点位置から次の時点(t+単位時間)の各節点の位置が算出でき、FEMの計算を、分子動力学ソルバを用いて実行可能となる。それでいて、分子動力学ソルバは無償利用可能なソフトを含めて並列計算機能を有するので、構造物のFEM解析をリーズナブルに並列計算で実現可能となる。
さらに、各節点の初期位置から時点(t)の位置に基づく変形行列Aから、歪変形行列Qと回転変形行列Rを求め、歪変形行列Qを用いて歪のみに起因する力を算出し、その後、力の向きを回転変形行列Rで変換しているので、幾何学的非線形に対応可能となる。
【0070】
第2実施形態のように、前記(c)における要素剛性行列算出部13による要素剛性行列Keの算出は、単位時間経過する毎に算出されることが好ましい。
【0071】
このようにすれば、要素剛性行列Keを単位時間経過する毎に毎回算出しても計算コストが大きくなく、全ての要素の要素剛性行列Keをメモリに記憶することによって生じるメモリの消費を抑制することが可能となる。
【0072】
第2実施形態のように、要素eの要素剛性行列Keは、式(1)~(3)で表現されることが好ましい。好適な実施形態である。
【0073】
第2実施形態のように、変形行列Aは、式(9)で表現されることが好ましい。好適な実施形態である。
【0074】
第2実施形態のように、回転を伴わない状態での歪に起因する、各節点の初期位置に対する変位を表す変位ベクトル(se、ue、we)は、式(11)で表現されることが好ましい。好適な実施形態である。
【0075】
第2実施形態のように、歪のみに起因して各節点に作用する力feは、式(12)で表現されることが好ましい。好適な実施形態である。
【0076】
第2実施形態のように、回転及び歪を伴った各節点に作用する力Feは、式(13)で表現されることが好ましい。好適な実施形態である。
【0077】
第2実施形態に係るプログラムは、上記方法をコンピュータに実行させるプログラムである。このプログラムを実行することによっても、上記方法の奏する作用効果を得ることが可能となる。
【0078】
上記の各実施形態で採用している構造を他の任意の実施形態に採用することは可能である。各部の具体的な構成は、上述した実施形態のみに限定されるものではなく、本発明の趣旨を逸脱しない範囲で種々変形が可能である。
【符号の説明】
【0079】
10 FEMモデル取得部
11 集中質量行列取得部
12 初期位置速度取得部
13 要素剛性行列算出部
114a 変形行列算出部
114b 右極分解部
114c 変位ベクトル算出部
115a 力算出部
115b 力方向変換部
16 分子動力学ソルバ