2C25

有限要素法を用いたマニピュレータの衝撃解析○守屋良昭，廣田直也(筑波大院)，磯部大吾郎(筑波大)Impact analysis of a manipulator using FEM* Yoshiaki MORIYA, Naoya HIROTA, Graduate School, Univ. of Tsukuba
Daigoro ISOBE, Univ. of Tsukuba
Abstract －In case of a manipulator operating at a high speed, a consideration of the performance of a drive
system such as an actuator is very important. However, the shock wave generated in the mechanism also becomes an important problem to prevent the destruction of itself. In this paper, a scheme for analyzing theshock wave generated in a manipulator is developed by using the Finite Element Method (FEM). In order tocheck the validity of the scheme, some impact analyses are carried out on a simple model. Then the shock
wave in a manipulator, which actively deforms its form, is analyzed.
Key Words: Impact Analysis, Finite Element Method, Manipulator

1. はじめに

2.1 Shifted Integration法

現在，マニピュレータはその生産性を向上させる
ため，今まで以上に安定し高速に稼動することが求められている．マニピュレータが高速で稼動する際，内部に発生する応力は増大する．また，高速移動時は接触によって思わぬ反力を衝撃的に受ける．このような応力増大や衝撃力が発生した場合，部品破損や振動による制御の困難を招く可能性がある．
これまでの研究では，マニピュレータやロボットといった能動的に変形する機構を完全剛体として扱うことが多く，それらが受ける衝撃に対し，機構の安全性を考慮した構造設計，姿勢決定，動作計画などについて議論されることは少なかった．今後，機構が動作することで発生する内力分布や，接触によって生じる衝撃力を設計段階で考慮することがますます重要になってくると考えられる．
そこで本研究では，マニピュレータのような能動的に動作する機構に対し，動作や接触を行う際に発生する，変形量や衝撃力を解析する手法を開発した．
本手法を用いて，マニピュレータをモデル化し，壁との接触を行う際の衝撃解析を行った．

2. FEM による解析アルゴリズム−1s10

本研究では，計算のツールとして有限要素法(FEM)を用いるが，解析を高速化するためにはモデル化の際に要素数の削減が必要とされる．そこで，機構のモデル化に際し，少ない要素数で高精度な解が得られ，構造工学分野の骨組構造解析でその有効性が確認されているShifted Integration法(SI法)1)を用いた．SI法では，有限要素と剛体・ばねモデルのそれぞれのひずみエネルギ近似式の等価条件を考察することにより，数値積分点の位置と回転・せん断ばねの位置の明確な関係が導出されている．あらかじめ関節が存在する位置が明確ならば，その位置に回転・せん断ばねを表現するように，入力データの段階で数値積分点をシフトする．Fig.1に，線形チモシェンコはり要素と，ヒンジまたは関節が回転ばね・せん断ばねによって陽に表現される剛体・ばねモデル(RBSM)との物理的等価関係を示す．両者の数値積分点位置と回転・せん断ばね位置との関係は以下の式のように表される．

s1 = − r1 または r1 = − s1(1)

s1 , r1 はそれぞれ数値積分点位置，回転・せん断ばね位置である．いま，1つのリンクを2つのチモシェ

1s

Link member
Center of gravity

Linear Timoshenko beam element−1r110r

Node No.1RBSM: NumericalTotal mass =M

Finite element model

integration point: Rotational and shear spring connecting rigid bars

Fig.1 Linear Timoshenko beam elementand its physical equivalent

第21回日本ロボット学会学術講演会（2003年9月20日〜22日）

Node No.3

Node No.2
Element No.1

Element No.2

M2

M

2
: Element Mass

Fig.2 Approximation of mass distribution
in a link memberンコはり要素でモデル化する場合，Fig.2のように，リンクの重心位置で節点1,2を持つ要素1と節点2,3を持つ要素2に分ける．このとき，2つの要素内で応力評価点をガウス積分点に相当する位置に配置することにより，曲げ変形を精度良く算出できる．そこで
本研究では，応力評価点に対応する数値積分点の値
を式(1)より s = m 1 − 2 3 3  とした．また，質量分布はFig.2に示すように，2つの要素によって構成される1本のリンクの全質量 M が重心位置に配置されるようにした．

2.2 運動方程式

仮想仕事の原理により，時刻 t = t + ∆t における運動方程式は，系の動作によって発生する慣性力を考慮すると以下のように定式化される．

{u&&}t + ∆t が求められたとき，変位量ベクトル {u }t + ∆t ，速度ベクトル {u& }t + ∆t を式(6)，式(7)のように求める．
{u }t + ∆ t = {u }t + {u }t ∆ t 1
+   − β  {u&&}t + β {u&&}t + ∆ t  ∆ t 221{u& }t + ∆ t = {u& }t + ({u&&}t + {u&&}t + ∆ t )∆ t2{u& m }t + ∆ t ={u&&m }t + ∆ t =(3){u& d }t + ∆ t =(4)

ここで， [M ] は全体質量マトリクス， {u m }は系の動作量ベクトル， {u d }は要素が変形することによって生じる変形量ベクトル， [C ] は全体減衰マトリクス，[K ] は全体剛性マトリクス， {F }は外力ベクトル，{R}は内力ベクトルである．また，全体減衰マトリク
ス [C ] は比例定数 a , b を用いて以下の式で表す2)．

[C ] = a[M ] + b[K ]

(5)

(7)

式(3)，式(4)，式(6)，式(7)より，動作量に対する速度，加速度ベクトル {u& m }t + ∆t , {u&&m }t + ∆t および変形量に対する速度，加速度ベクトル {u& d }t + ∆t , {u&&d }t + ∆t について次式が得られる．

[M ]{u&&m }t +∆t + [M ]{u&&d }t +∆t + [C]{u&d }t +∆t + [K ]{∆ud }(2)= {F}t +∆t − {R}t{u m }t + ∆t = {u m }t + {∆ u m }{u d }t + ∆t = {u d }t + {∆ u d }(6){u&&d }t + ∆ t =
− 1  {u&&m }t(8)1{∆ u m }−  1 − 1  {u& m }t −  1 − 4 β  {u&&m }t ∆ t2 β∆t 2β
 4β (9)1β∆t2{∆ u m }−1β∆t{u& m }t −  1 2β− 1  {u&&d }t
(10)1{∆ u d }−  1 − 1  {u& d }t −  1 − 4 β  {u&&d }t ∆ t2 β∆t 2β 4β (11)1β∆t2{∆ u d }−1β∆t{u& d }t −  1 2β

式(9)，式(10)，式(11)を式(2)に代入すると，運動方
程式は次式のようになる．
1
[K] +[M] + 1 [C]{∆ud } = {F}t +∆t −{R}t + [M] 1 {u&d }t +  1 −1{u&&d }t β∆t 22β∆t  2β  β∆t 1

1{u&m}t −  1 −1{u&&m}t  + [C] 1 −1{u&d }t + 1− 4β {u&&d }t ∆t − [M ] 2 {∆um}−∆∆
ββ2β2β4βtt




(12)

2.3 断面力の計算過程

FEMによる動的シミュレーションを行う場合には， 式(12)を {∆ u d }について解き，累積することにより変形量ベクトル {u d }t + ∆t を求める．また，入力データか解を逐次的に求める必要がある．その手順として，
ら動作量ベクトル {u m }t + ∆t は既知である．次に {u d }t + ∆ t ，マニピュレータの各節点の移動による位置データを{u m }t + ∆t を式(8)，式(9)，式(10)，式(11)に代入し，加毎ステップ入力し，そこから算出される動作・変形{u&&d }t + ∆t , {u&&m }t + ∆t および速度ベクトル
速度ベクトルによる速度・加速度データを次ステップの解析に逐&&{}{u}u,d t + ∆tm t + ∆ t を求める．さらに，式(12)で求まっ次的に加算する．これにより，マニピュレータの内た変形量ベクトル {u d }t + ∆t を要素座標系の情報へ変換
力分布の時刻歴解析が可能になる．本研究では，SIし，まとめて {δ }t + ∆t とおくと断面力 {σ }t + ∆t は法を用いてモデル化したため，少ない要素数で解析が可能であることから，陽解法と陰解法での各増分ステップにおける計算時間の差は大きくならない．
{σ }t + ∆t = [D ][B ]{δ }t + ∆t
(13)
また，陽解法では時間増分がCourant条件を満たす小さな値に設定する必要がある．したがって，全体とによって求められる．式(13)で [D ] はひずみ変位マトしては陰解法による計算時間の方が短くなると考え，リクス，[B ] は応力ひずみマトリクスである．以上の動的問題を計算するための時間積分スキームとして， ような手順より，マニピュレータ内に発生する衝撃代表的な陰解法の一つであるNewmarkのβ法を用い力の解析を行う．
ることにした．
Newmarkのβ法では，時刻 t における変位量ベクト
3. アルゴリズム検証
ル {u }t ，速度ベクトル {u& }t ，加速度ベクトル {u&&}t を既知とする．時刻 t + ∆ t において加速度ベクトル

3.1 準静的な外力作用下の内力解析

最初に本解析手法の有効性を示すために，Fig.3のように2リンク機構のマニピュレータを加速円運動させた時に生じる軸力を解析する．本解析に用いたマニピュレータの材質はアルミニウムを想定し，そのリンク部材の詳細をTable.1に示す．マニピュレー
タをモデル化する際，1リンク当り2つのはり要素で分割・モデル化し，リンクの重心位置に1リンク分の全質量 M を配置した．また，時間増分 ∆ t は0.01[ms]とし，動作による速度・加速度データを直接入力した．解析の結果，Fig.4(a), (b)に示すように，軌道より求まる遠心力とFEMで求めた軸力が完全に一致することがわかった．このことから準静的な外力に対し，内力分布の妥当な解を算出することが確認できた．

た．解析では時間増分 ∆ t を 0.01[ms]とした．Fig.6(a)は，このときに根元要素に発生する軸力の波形である．Fig.6(b)は載荷直後の波形を拡大したものである．
ここで，理論的な縦振動周期 T は次式となる．

3.2 撃力作用下の内力解析

4. マニピュレータの衝撃解析

次に，Fig.5のように2リンクマニピュレータを直線状に伸ばした状態で，先端部に軸に対して内向きに10.0 [N]のステップ荷重を与えた．減衰を考慮するため，式(5)の比例定数の値を a = 2.0 × 10 2，b = 0.0 とし

YXMMLL：Stress evaluation points
at Gaussian point

M : Mass
L : Length: Link joint

Fig.3 Finite element analysis of a rotating manipulator

T = 8Lρ(14)E

式(14)より，このマニピュレータの理論的な縦振動周期は T = 3.22 × 10 −4 [s]と計算される．ただし ρ は部材密度， E はヤング率を表す．Fig.6(b)より，解析された波形は周期 3.3 × 10 −4 [s]を持ち，理論値に近い値をとっていることが分かる．このことから，撃力に対しても妥当な固有周期をもつ振動が算出できることが確認できた．

4.1 解析モデル
本アルゴリズムを用い，Fig.7のような2リンクマニピュレータが水平面内で壁と接触をする際に発生する衝撃力の解析を行った．マニピュレータの目標軌
道は，Fig.7において，初期位置から①の動作をした後に一定時間静止し，その後②の動作をして再び初期位置に戻ってくるという軌道である．今回の解析
では，時間増分∆ t を0.01[ms]，全体の動作時間を2.5[s]とし，①の軌道で1.0[s]，静止時間が0.5[s]，②の軌道を1.0[s]という時間で目標軌道上を動作する．
部材の詳細はTable.1に示す．また，運動量の変化が全て力積に変わるものと仮定し，その力積が部材に加えられる時間を∆ t f とした．本解析では動作による速度・加速度データを直接入力した．

F

Table.1 Parameter of a link member
Parameter

Value

Length: L[ m ]

MLL

0.40

Sectional area: A[ m ]1.92 × 10 −4

Mass: M[ kg ]0.108

Young’s modulus: E[ N m 2 ]6.90 × 1010

Moment of inertia:I[ m 4 ]1.64 × 10 −82X

F ：10.0[N]M : Mass
L : Length：Stress evaluation points
at Gaussian point: Link joint

Fig.5 Finite element analysis of a manipulator understep load

1.5

0.5

0

1-100.5

0.5

1Time[s]1.5

2-10-20000
Normal Force[N]10Normal Force[N]Centrifugal Force[N]1.5Normal Force[N]MY00.5

1Time[s]1.5-202(a) Normal force by FEM(b) Centrifugal force
Fig.4 Time histories of normal force obtained by FEM
and centrifugal force

00.02

0.04 0.06
Time[s]0.08

0.1

00.001
Time[s]0.002(a) Entire process(b) Enlargement
Fig.6 Time histories of normal force0.6
Normal Force[N]②①
M0.6
Normal Force[N]

Wall

0.4

0.2

0-0.2

-0.20X：Stress evaluation points
at Gaussian point

0.20M

YL

0.4

1
Time[s]201
Time[s]2L(a) Base element(b) Tip element
Fig.8 Time histories of normal force (Vx≒0.0[m/s])

M : Mass
L : Length: Link joint0.2

初めに，壁との接触時に完全静止する目標軌道で，本研究で開発した手法により，動作時のマニピュレータ内に発生する衝撃力の解析が可能であることが確認できた．今後の課題として，高速化のため解析アルゴリズムのさらなる検証を行う必要がある．
また，実機に発生する振動を，AEセンサなどによって計測し，本解析アルゴリズムで求められた波形との比較・検討を行っていく予定である．

0.4

0.2

00-0.2∆ t f =0.01[ms]として解析を行った．そa =200， b =0，5. まとめ

Normal Force[N]0.4

接触時の衝撃解析-0.2
01
Time[s]201
Time[s]2(a) Base element(b) Tip element
Fig.9 Time histories of normal force (Vx=-0.015[m/s])

50
Normal Force[N]Normal Force[N]500-500-5011.02

Time[s]

1.04

1.06

11.02

1.04

Time[s]

1.06(a) Base element(b) Tip element
Fig.10 General outline of the shock wave
( a =200, ∆ t f =0.01[ms] ) (Vx=-0.015[m/s])
50
Normal Force[N]の結果をFig.8に示す．なお，Vxは接触直前のマニピュレータ先端のX軸方向における速度である．Fig.8
を見ると，①や②の動作の際には遠心力などによって部材に内力が発生していることが分かる．また，接触時にはほとんど速度を持たないため，大きな衝撃力は発生していないことが確認できる．
次に，壁との接触時に速度(Vx=-0.015[m/s])を持つ目標軌道で解析を行った．その結果をFig.9に示す．
Fig.8の場合と比較すると，小さな速度で接触しているにも関わらず大きな衝撃力が発生していることが分かる．衝撃波形全体図をFig.10に示す．Fig.10をみると，根元の要素，先端の要素ともに同様な衝撃力が発生していることが分かる．
最後に，減衰マトリクスの定数 a と力積を与える
時間長さ ∆ t f の条件をそれぞれ変えた場合の，先端要素に発生する衝撃力の解析結果を示す．∆ t f =0.1[ms]とした解析結果をFig.11(a)，∆ t f は変化させず減衰定数を a =2000とした解析結果をFig.11(b)に示す．∆ t fを大きくしたFig.11(a)の場合は，Fig.10(b)と比べて緩やかに力が与えられるために，発生する衝撃力のピーク値が小さくなっている．また，減衰定数を大きくしたFig.11(b)の場合には，Fig.10(b)と比べてピーク値の値は同じたが，より早く振動が減衰していることがわかる．今後，これらの値を実機での実験により決めていく必要があると考えられる．

0.650
Normal Force[N]4.20.6Normal Force[N]Fig.7 Finite element analysis of a manipulator in contact
with a wall

0-50

0-50
1

1.02

Time[s]

1.04

1.06

11.02

Time[s]

1.04

1.06

(b) a =2000, ∆ t f =0.01[ms](a) a =200, ∆ t f =0.1[ms]Fig.11 Time histories of normal force in tip element(Vx=-0.015[m/s])
参考文献
1) 都井 裕, 骨組構造および回転対称シェル構造の有限要素解析における Shifted Integration 法について, 日本
造船学会論文集, 第 168 号, (1991), p.369-377.

2) 鷲津 久一郎 他, 有限要素法ハンドブック Ⅰ基礎編,
培風館, p.133-134