825
梁の動力学解析における節点座標のとり方とその解析結果の比較
（Absolute Nodal Coordinate 法と Large Rotation Vector 法の比較）
Comparison of formulations and calculated results between Absolute Nodal Coordinate formulation
and Large Rotation Vector formulation for the dynamic analysis of flexible beams
⃝正 高橋義考（いわき明星大） 正 清水信行（いわき明星大）
Yoshitaka Takahashi, Iwaki Meisei University
Nobuyuki Shimizu , Iwaki Meisei University
There are three basic finite element formulations which are used in multibody dynamics for solving the beam
problem undergoing large rotation and deformation. These are the floating frame reference approach, the incremental
method and the large rotation vector approach. And recently the absolute nodal coordinate formulation was proposed.
This formulation is said to have advantageous characteristics for computer calculations.
In this paper, the properties of the two methods, i. e. the large rotation vector formulation and the absolute
nodal coordinate formulation are compared and discussed from the theoretical aspect and the computational aspect for
a simple planar beam motion
K ey Wo rds : multibody dynamics, flexible beam, nodal coordinate, large rotation and deformation
1 . はじ めに
大変形をともなう弾性梁の動力学解析法として、線形有限要
素法（FEM）を応用した非線形有限要素法が研究されている。
これらの研究では、FEM を用いて弾性梁をモデル化するのに、
研究者により節点座標のとり方に工夫があり興味深い。例え
ば、
・floationg frame of reference（浮動基準枠）
・ convected coordinate system（対流座標系）
・ finite segment method（有限セグメント法）
・ large rotation vector method（大回転ベクトル法）
などが一般に扱われている (1) 。そして近年、コンピュータを
用いた解析において効率よく計算できる手法として、Shabana
らにより Absolute Nodal Coodinate method が発表された
( 2 ) 。従
来の手法と比べると、質量マトリクスが時間に対して不変とな
る点に大きな特徴がある。これは、運動方程式を数値計算で解
く場合、質量マトリクスの逆行列を１回だけ求めれば良く、計
算時間が有利であると言われているが従来の手法と具体例に対
して、その解析が比較検討されたことがないため、どのような
特徴を持っているかは明らかでない。
本研究では Absolute nodal coordinate 法（以下 A.N.C）と Large
rotation vector 法（以下 L.R.V.）の２つに注目し、弾性梁の平面
動力学シミュレーションを行うことで両手法の比較検討を行
い、その特徴を明らかにする。
の初期状態からのX,Y方向の全体座標からみた変位を u 1 , u 2 と
し、角度(rad)の変化を u 3 とする。同様に、節点 B の初期状態か
らの対応する変位と角度の変化を u 4 , u 5 , u 6 とする。
要素内の任意点 P の全体座標系から見た位置ベクトル r は
r = u 1 i + u 2 j + ξ i ′ + u ′p1 i ′ + u ′p 2 j ′ と表される。ここで、ξは要素座標系において要素が変形して
いないとしたときの P 点のX'軸の座標である。u ′p1 , u ′p 2 はこの点
から変形している実際の任意点 P までの要素座標系である。
u ′p1 , u ′p 2 は、ベルヌーイ・オイラー梁の形状関数
u ′p1 =
ξ
l
( 3 ) を用いて、
u4′ 2
3
2
3
ξ
ξ
ξ
ξ
− 2 u ′5 + l − + u ′6
l
l l
l
u ′p 2 = 3
(2)
(3)
と表すことができる。ここで、u4′ , u ′5 , u ′6 は、要素座標系に対し
て変形していない状態の要素を点 O' に左端節点 A を一致させ
X' 軸上に沿って置いた場合の右端節点 B から、実際の節点 B ま
での要素座標系の変位である。同様に節点Aについての要素座
Y′
B
u4
u6
Y
P
X′
A
u1
2 . 解析 理論
2 .1 L . R. V . 法 によ る定 式化 図 1 に 2 次元弾性梁の１要
素を示す。全体座標系を O-XY とし、要素の一端に固定した要
素座標系を O'-X' Y' とする。要素座標系の X' 軸は要素左端の節
点 A での接線方向としている。ここで扱う弾性梁はベルヌー
イ・オイラーの仮定が成り立つものとする。全体座標系のX、Y
(1)
u5
u2
Initial Condition
O
方向の単位ベクトルを i , j とし、要素座標系の単位ベクト
ルを i ′ , j ′ とする。要素が初期状態から移動したとき、節点 A
日本機械学会〔No.98-8〕機械力学・計測制御講演論文集〔'98.8.17∼20・札幌〕
ξ
u3
O′
P
X
ξ
l
Fig.1 Nodal coordinate of L.R.V
標 の 変 位 を u1′ u2′ u3′ と す る 。 こ の 要 素 の 節 点 の 要 素 座 標
u′
T
= {u1′ u2′ u3′ u4′ u5′ u6′ } と、全体座標 u には次のような関係が
d y / dr を e3 , e4 、梁の右端節点 B の X,Y 座標を e5 , e6 、勾配
d x / dr , d y / dr を e7 , e8 とする。
ある。
要素内の任意点 P の座標は次式で与えられる。
u1′ = u2′ = u3′ = 0 (4)
r = S e (
)
u4′ = (u4 − u1 + l ) cos u3 + u5 − u2 sin u3 − l (5)
(6)
u5′ = − (u4 − u1 + l ) sin u3 + (u5 − u2 ) cos u3 (7)
u6′ = u6 − u3 運動方程式を導出するにあたり、ハミルトンの原理を用い
る。まず、この要素の運動エネルギーは次式により得られる。
T =
1
2
∫ ρ r˙ r˙ dv T
(8)
V
ここで r˙ は r の時間による１階微分であり、梁上の任意点 P の
速度を表している。
次に時間 t 1 ∼ t 2 間において運動エネルギー
の変分をとり、
式(4),(5),(6),(7)を用い要素節点の変位を全体座標
u T = {u1 u2 u3 u4 u5 u6 } で表すと
δ T =
∫ δu
t2
T
t1
( Mu˙˙ + Qm ) dt (9)
となる。ここで M は質量マトリクス、Q m は遠心力、コリオリ
力成分である。
要素のポテンシャルエネルギーは節点の要素座標変位 u ′ を用
いて
1
u ′ T K u ′ (10)
2
と表される。、ここで K は要素の剛性マトリクスである。時間
t 1 ∼ t 2 間においてポテンシャルエネルギーの変分をとり、式(4)
,(5),(6),(7)を用い要素節点の全体座標変位 u を用いて書きなお
U =
すと
δ U =
∫ δu
t2
T
t1
Q k dt (11)
となる。ここで Q K は梁の復原力成分である。ハミルトンの原
理は、式(9)と式(11)を用いると次のようになる。
∫ δu
t2
t1
T
( Mu˙˙ + Qm − Q k ) dt = 0 (12)
よって、次式の運動方程式が導かれる。
˙˙ = Q (13)
Mu
Q = Q f − Qm + Qk
ここで S は形状関数であり次式で示される。
T
2
3
0
1 − 3ξ + 2ξ
2
3
1 − 3ξ + 2ξ
0
l(ξ − 2ξ 2 + ξ 3 )
0
l(ξ − 2ξ 2 + ξ 3 )
0
S=
0
3ξ 2 − 2ξ 3
(16)
3ξ 2 − 2ξ 3
0
3
2
0
l(ξ − ξ )
l(ξ 3 − ξ 2 )
0
ここで、l は梁の長さ、x は要素節点 A から任意点までの長さ、
ξ = x / l である。まず、慣性力を求めるために運動エネルギー
を求める。
T =
M =
r˙ dV =
1 T
e˙
2
∫ ρS
V
T
S dV e˙ (17)
∫ ρS
V
T
S dV (18)
d ξ (19)
0
により得られる。ここで、E は弾性係数、a は断面積、I は断面
U =
1
2
∫
l
2
2
∂u
∂ 2u
E a l + E I 2t
∂ξ
∂ξ
二次モーメントである。歪みエネルギを節点座標 e で偏微分す
ることにより弾性力が得られる。
QK = Ke = A11ei x2 + A22 ei y2 + ( A12 + A21 )ei x i y + B11ej x2
+ B22 ej y2 + ( B12 + B21 )ej x j y − A1i x − A2 i y
}
なわち、梁の左端節点 A の X,Y 座標を e1 , e2 、勾配 d x / dr ,
Y
B dr
dy
dx
dy
A dx
e1
T
量を ut とする。梁の歪みエネルギは
(14)
{
O
V
る要素までの、X'軸方向の軸変形量を ul 、Y'軸方向の曲げ変形
標による節点変位を e T = e1 e2 e3 e4 e5 e6 e7 e8 で表わす。す
dr
∫ ρ r˙
となる。ここで、S は時間に依存しないことから、質量マトリ
クスは時間に対して不変であることが分かる。
次に、弾性力を導出するために要素座標系を考える。要素両
端を結ぶ線を X' 軸とし、それに垂直方向を Y' とする。X' 軸に
変形していない要素を想定し、この要素から実際の変形してい
2 .2 A .N .C . 法 によ る定 式化 図 2に全体座標系 O-XY 内に
置かれた 2 次元弾性梁の１要素を示す。この要素節点の全体座
e2
1
2
式(17)より、質量マトリクスは
ここで、Q f は外力項である。
e6
(15)
e5
Fig.2 Nodal coordinate of A.N.C.
X
1
∂i
+ e T A11ei x + e T ( A12 + A21 )ei y − A1T e x
∂e
2
T
∂i y
1
+ e T A22 ei y + e T ( A12 + A21 )ei x − A2T e
2
∂e
T
1
∂j
+ e T B11ej x + e T ( B12 + B21 )ej y x
∂e
2
T
∂j y
1
+ e T B22 ej y + e T ( B12 + B21 )ej x
∂e
2
T
(20)
となる。ここで、A11 から B22 は、弾性係数、断面積、断面二次
モーメント、形状関数によって定義される定数マトリクスであ
る (2) 。
梁の運動方程式は式(18),(20)を用いて
M e˙˙ = Q (21)
Q = Q f − QK (22)
となる。
5
5
4
4
3
3
Y [m]
3. 1 要素内の任意 点の位置ベクトル 前章で述べた両手法
の定式化では、弾性梁はベルヌーイ・オイラーの仮定が成り立つ
とし、形状関数は同じとしている。L.R.V. 法では左端の節点での
接線が全体座標から傾く剛体回転角を定義し、
そこに固定した要
素座標系に対して節点の勾配は d y / dr をθと近似している。ま
た、軸方向の変位については、要素内において均一に軸方向に変
形していると仮定している。それに対して A.N.C. 法では勾配
Y [m]
3 . シミ ュレ ーシ ョン によ る比 較
2
1
0
1
0
1
2
3
X [m]
d x / dr 、d y / dr を近似せず、そのまま節点座標として用いてい
u T = { 1 , 1 , π / 4 , − 0.46 , 4.54 , π / 4 }
4
4
3
3
Y [m]
5
2
1
π
π
π
π
eT =
1 , 1 , cos , sin , 4.54 , 4.54 , cos , sin (24)
4
4
4
4
0
0
1
π
2
, sin
π
(26)
2
図 4 を比較すると、A.N.C. 法は梁両端が指定した角度になって
いるが、L.R.V. 法は梁右端の角度が一致していないことが分か
る。この原因は軸方向の変形量と、勾配を近似的に扱っているた
めである。
要素内の任意点にこのような誤差が生じることにより、運動
エネルギ、ポテンシャルエネルギにも誤差が生じるため、運動
方程式に誤差が出ると考えられる。そのため、L.R.V.法で変形の
大きい弾性梁の動力学解析を行う場合、梁の要素分割数を十分
に多くし、一つの要素の変形量を微少な範囲に押さえなければ
ならない。それに対して A.N.C. 法の場合は要素内の任意点の表
し方に近似式を用いていないので、
L.R.V.法ほど要素分割数を多
くしなくても解析できることが分かる。
3. 2 振子 の動力学解析 梁の動力学解析の例として振子を
取りあげ、両手法によりシミュレーションを行い比較検討する。
振子は曲げ剛さの大きい梁(case1)と、小さい梁(case2)の２種類を
想定し、材料特性は両梁とも長さ 1m、密度 5540 kg/m 3 、断面積
0.0018 m 2 、弾性係数 0.70 × 10 7 N/m 2 とし、断面二次モーメントだ
け c a s e に より 変 え た。 c a s e 1 は 1.215 × 10 − 6 m 4 、c a s e 2 は
1.215 × 10 − 8 m 4 であり、要素分割数を 5 とした。数値積分法には
Runge-Kutta-Gill 法を用い、タイムステップは 1 × 10 − 5 秒とした。
2
3
X [m]
(a) L.R.V.
2
3
X [m]
4
5
4
5
2
4
5
0
0
1
2
3
X [m]
(b) A.N.C.
Fig.4 Deformed configuration of beams
0.4
Y [m]
0
-0.4
XXXXXX X
X XXXX XX XXXXXXXXXXX
X
X
X XXX XXX X X
X X X XX X XX
X
X XX XX
0.0s
0.9s
0.3s
-0.8
0.6s
-1.2
-1.2
0
0.4
0.8
1.2
X [m]
Fig.5 Deformed configuration of the pendulum(case1)
calculated by L.R.V.
-0.8
-0.4
0.4
0
Y [m]
eT =
1 , 1 , cos 0 , sin 0 , 4.54 , 4.54 , cos
(25)
1
1
とした場合の梁要素を図 3 に示す。これは長さ５ m の梁を全体
座標に対してπ /4 [rad] 傾けた状態を表す。図 3 より、梁が変形
していない状態では、両手法の梁要素は一致することが確認で
きる。
本来 L.R.V. 法では微少変形を想定しているが、両手法を比較
するために、梁を大きく変形させた状態を考える。図 3の状態か
ら、要素左端の角度を 0 [rad]、右端の角度をπ /2 [rad]と変形さ
せた場合を考える。L.R.V. 法の節点座標、A.N.C. 法の節点座標
u T = { 1 , 1 , 0 , − 0.46 , 4.54 , π / 2 }
0
(b) A.N.C.
5
(23)
を次式に示し、式(1),(15)で計算した梁の状態を図 4 に示す。
0
5
Fig.3 Configuration of rigid beams
Y [m]
これらの違いにより、
要素内の任意点の位置ベクトルにどれだ
け差が出るかを式( 1) ,( 1 5) を用いて具体的に計算し確かめる。
L.R.V. 法と A.N.C. 法において同じ条件となる節点座標の値を代
入し、要素内の任意点を図に示す。最初は弾性梁が変形してい
ない状態として、L.R.V. 法の節点座標、A.N.C. 法の節点座標を
4
(a) L.R.V.
る。これは梁要素内の任意点 Pの位置を表すのに剛体回転を用い
ていないため、梁要素が全体座標に対して有限回転を有するとき
dr/dx= θと近似できないためである。また、任意点の X、Y 座標
を形状関数を用いているため、近似法は入っていない。
2
-0.4
-0.8
X
XXXX XX XXXXXXXXXXXXX
XXXXXXXX
X
X
X
X
X
X
X
X
X
X X X X XX
XX X XX
X
-1.2
-1.2
0.0s
0.9s
0.3s
0.6s
-0.8
-0.4
0
0.4
0.8
1.2
X [m]
Fig.6 Deformed configuration of the pendulum(case1)
calculated by A.N.C.
16
0.0s
0.3s
Y [m]
-1
0.6s
-2
-3
-4
-0.8
-0.4
0
X [m]
0.4
0.8
1.2
0.4
X XXXXXX
XXXXXXXXXXXXXXXXXXXX
XXXX XX X
X
X
X X XX XXXX X XX
XX X X
0.0s
Y [m]
-0.4
-0.8
-1.2
-1.2
0.9s
0.3s
0.6s
-0.8
-0.4
0
0.4
0.8
1.2
X [m]
Fig.8 Deformed configuration of the pendulum(case2)
calculated by A.N.C.
0.4
Y [m]
0
50 elements
-1.2
-1.2
10
8
6
4
2
1.0
10 15 20 25 30 35 40 45 50
Element number
Fig.10 The ratio of the calculation time
between L.R.V and A.N.C
5
3 .4 計算時 間の 比較 L.R.V.法に比べてA.N.C.法は質量
マトリクスが時間に対して一定になる特徴がある。そのため、
運動方程式を数値積分で解くときに質量マトリクスの逆行列
を計算するのは１回だけになり、計算時間は有利になる。
これを見るために両手法の計算時間の比較を行う。解析内
容は前節で行った case1 の振子とし、計算条件、材料特性も同
じものを用いた。A.N.C. 法での計算時間に対する L.R.V. 法の
計算時間の比と、要素数の関係を図 10 に示す。A.N.C. 法の計
算時間は要素数に比例するのに対して、L.R.V. 法は要素数の
二乗に比例して計算時間が増える。そのため、節点数の多い
計算の場合は A.N.C. 法が有利と思われる。しかし、A.N.C. 法
では梁の剛性が高い場合、数値積分の安定性が悪く、時間刻
みを十分に小さくしなければならないことが分かった。その
ため、固有振動数の高い梁を解析する場合には L.R.V. 法が有
利と思われる。
4 . 結 論
0.0s
0.9s
-0.4
-0.8
12
0
Fig.7 Deformed configuration of the pendulum(case2)
calculated by L.R.V.
0
14
0
0.9s
-5
-1.2
Ratio
XXXXXXX
XXXXX
X
XXXX XXXXXX
XX
X XX
X
XXX X X X
XX X X
XX X X
X
X
0
(L.R.V./A.N.C)
1
0.3s
0.6s
-0.8
-0.4
0
X [m]
0.4
0.8
1.2
Fig.9 Deformed configuration of the pendulum(case2)
calculated by L.R.V.
case1 のシミュレーション結果を図 5、6 に示す。また、case2
のシミュレーション結果を図 7、8 に示す。
case1 は梁の変形が少ないため、両手法のシミュレーション
結果は良く一致している。しかし、case2 は梁の変形が大きい
ため、L.R.V.法では解析がうまくできないことが分かる。これ
は要素分割数を増やすことで解決できる。図9に要素分割数を
50 にした場合の結果を示す。一方、A.N.C. 法の場合は要素分
割数が5でも解析できている。本研究の計算にはワークステー
ション HITACHI 3500 を用いた。case1 の計算時間は L.R.V. 法
で 6 分 40 秒、A.N.C. 法で 12 分 30 秒であった。
梁の動力学解析において、本論文で示した L.R.V. 法と、 A.N.C. 法の比較を行った。L.R.V. 法では要素座標系において
微少変形を仮定しているため、軸方向の変形量、勾配を近似
的に扱っている。この仮定により、要素座標系において従来
の有限要素法の剛性マトリクスを使うことができ、ポテン
シャルエネルギを簡単に求めることができる。しかし、要素
の変形が微小変形の仮定の範囲を超えると、正しく解析でき
ない。その解決策として梁の要素分割数を多くしなければな
らず計算時間が不利となる。
一方A.N.C.法は要素内の任意点を表すのに近似的手法を用
いていないため、梁の要素分割数が少なくても柔軟な梁の解
析ができる。しかし、ポテンシャルエネルギーを求める計算
が複雑となる。また、固有振動数の高い梁を解析する場合、解
を安定に解くためには数値積分の時間刻みを十分に小さくし
なければならないことが分かった。
参 考文 献
(1) A.A. Shabana, Flexible Multibody Dynamics:Review of Past and Recent
Developments, Multibody System Dynamics, 1997, vol.1, 189-222
(2) A.A. Shabana, H.A. Hussien, J.L. Escalona, Application of the Absolute
Nodal Coordinate Formulation to Large Rotation and Large Deformation
Problems, ASME, Journal of Mechanical Design, vol.120(1998),1-9
(3)鷲津久一郎 , エネルギ原理入門 , 培風館 ,1980, 56-71
© Copyright 2024 ExpyDoc
