計算工学講演会論文集

Vol.9 (2004年5月)

計算工学会

有限要素法によるフレキシブルマニピュレータの軌道計算
Calculating Trajectories for Flexible Manipulators by Using Finite Element Method○ 加藤 昭博1)， 磯部大吾郎2)
Akihiro KATO and Daigoro ISOBE
1)筑波大学大学院 （〒305-8573 茨城県つくば市天王台1-1-1, akato@kz.tsukuba.ac.jp）
2)工博

筑波大学助教授 機能工学系（〒305-8573 茨城県つくば市天王台1-1-1, isobe@kz.tsukuba.ac.jp）

In this paper, a scheme using the Finite Element Method (FEM) for calculating trajectories of flexible
manipulators is proposed. The trajectories are required in inverse dynamics calculation using the formerly
developed parallel solution scheme, where the inverse dynamics can be easily calculated regardless of the rigidity or
the connectivity of the constituting link members. The proposed scheme is formulated based upon the principle ofvirtual work in consideration of motion of the entire system, by dividing nodal displacements into values ofoperation distance and deformation. It enables us to discuss on forces generated by rotational motions. The cubic
beam element based on the Bernoulli-Euler hypothesis is used for modeling link members in order to improveaccuracy against elastic bending deformation. In this paper, some numerical tests and experiments are carried out ona simple flexible link system, and the results are compared to show the validity of the proposed scheme.
Key Words: Flexible Manipulators, Trajectory, Finite Element Method

1．
．はじめに連続体力学に基づく数値解析手法として広く利用される有限要素法（FEM）は，全体直交座標系に基づいて系全体を微小要素に離散化し，全体方程式にまとめて解を求める汎用手法である．そのため，各要素の情報が並列的に取り扱われ，系の構造の複雑さ如何に関わらず，節点力や変位・ひずみ等の算出が可能となる．この特長を利用し，リンクの形態に依存せずに逆動力学が計算可能な並列的解法が開発された[1]．従来の手法の一つ，ニュートン・オイラー法が相対回転座標系に基づいて再帰的な処理により動力学方程式を求める直列的なアプローチをとるのに対し，離散化方程式から全関節トルクを一気に求める本解法は，並列的なアプローチをとると言える．
本解法では，動的な効果を含めた各節点力から，力学的な関係に基づいて各関節トルクが算出される．また，従来の動力学方程式が全ての成分が混在した形をとるのに対し，本解法のトルク算出式は，力の次元に関する成分，座標変換に関する成分，長さの次元に関する成分ごとのマトリクスに分離される．そのため，柔軟性・拡張性に富み，機構構造が複雑な場合や制御の途中で構造が変化する場合にも，入力データの変更のみで対応できる．これより，従来の逆動力学計算法では困難であった閉ループ系や，連続的に開から閉ループ系などへと形態が移行する機構に対する逆動力学計算[2]，さらにこれらの機構のフィードフォワード制御[3]などが，動力学方程式や計算アルゴリズムを変更することなく実行可能となった．
ロボットの動作のさらなる高速化が求められている現状
において，これは非常に重要な意義を持つものである．
他方，ロボット機構がタスクを行う際のエネルギ消費量の削減やロボット自体のコスト削減を目的として，機構のさらなる軽量化の必要性が高まっている．しかし，部材の軽量化はリンク剛性の低下を招き，機構内に発生する弾性たわみや曲げ振動が無視できなくなるという問題が生じる．この問題は，機構を完全剛体と仮定して設計することが多い現状において，その制御を非常に困難なものとする要因の一つとなっている．そのため，弾性たわみを考慮したロボットアームのモデル化[4][5]や逆動力学計算[6]などの研究も盛んに行われている．しかし，従来の動力学方程式ではトルクの次元で求解するため，変数分離された形をとらず，さらにそれに加えて弾性たわみの影響を考慮する必要がある．したがって，その逆動力学計算には大変複雑な処理を要する．一方，前述の並列的解法では，動的な効果を含めた節点力を求めれば全関節トルクが算出されるため，剛体リンク系と柔軟リンク系との間では，その計算アルゴリズム自体に相違が無い．この特徴を利用して，並列的逆動力学計算法に対し，有限要素によって弾性たわみを考慮した軌道計算アルゴリズムを付加することで，柔軟リンク系に対しても適用可能な解法が開発された[7]．また，軌道計算と逆動力学計算を同一モデルで行うことにより，リンク機構に任意のタスクを与えた場合に，その軌道の算出から必要な関節トルクの算出までの過程を包括的に一つのアルゴリズムとして扱うことが可能になった．
本報告では，文献[7]の軌道計算アルゴリズムに対し，Bernoulli-Euler の理論に基づく 3 次はり要素を導入して曲げ変形精度を高め，計算のさらなる高速化を図った．さらに，実機による検証を行い，アルゴリズムの有効性を検証した．

２．軌道計算アルゴリズム
２．軌道計算アルゴリズム
本節では，フレキシブルマニピュレータに任意のタスクを与えた際に，その動作形態を求める軌道計算アルゴリズムについて解説する．本アルゴリズムは，骨組構造解析に広く用いられる解析手法を，能動的に動作する機

0.2
0.1
0-0.1
0(4)

0.1

0.2
time(sec)

0.3

0.4

(a) Lumped mass matrix

(1)
(2)
(3)

ここで， [M ] は全体質量マトリクス， [C ] は全体減衰マト
リクス，[K ] は全体剛性マトリクス，{F } は外力ベクトル，{R} は内力ベクトル， {um } は系の移動量のみのベクトル，{u d } は系の変形量のみのベクトルである．また，全体減衰マトリクス [C ] は比例減衰を用いることとし，以下の式で定義した．
[C ] = α [ M ] + β [ K ]0.3

1element
2elements
4elements
8elements
16elements

0.4
0.3

X(m)

&&
[M]{u&u&u&∆ud}={F}t+∆t −{R}tm}t+∆t +[M]{d}t+∆t +[C]{d}t+∆t +[K]{{u m }t + ∆t = {um }t + {∆u m }{u d }t + ∆t = {u d }t + {∆u d }1element
2elements
4elements
8elements
16elements

0.4x(m)構に適用したものである．フレキシブルマニピュレータの挙動を動的非線形問題として定式化し，その増分理論と し て は 非 線 形 問 題 で 一 般 的 に 使 用 さ れ る Updated
Lagrangian Formulation (U.L.F)を用いた．また，時間積分法には，代表的な陰解法の一つである Newmark のβ法を
δ = 1 / 2 ， β = 1 / 4 として使用した．リンク部材のモデル化には前述の 3 次はり要素を用い，曲げ変形に対する精度の向上を図った．
系の移動によって生じる慣性力を考慮すると，仮想仕事の原理より，時刻 t + ∆t における運動方程式は以下のように定式化される．

0.20.1
0-0.1
00.1

0.2
time(sec)

0.3

0.4(b) Distributed mass matrix
Fig.2 Numerical estimations

ここで用いる比例定数α，βは，後に実験からフィードバックして決定する．式(1)に，Newmark のβ法を適用すると次式が得られる．

1 − 4β1&d }t ∆t )){u&− 1){u&d }t + (2β4β

0.4

0.2

0.2

00-0.2

-0.2

-0.4

X(m)
00.2

(a) Rigid body model

0.4

Y(m)

-0.4

X(m)
0(6)

３．数値的評価ることにより，変形量増分 {∆ud } が得られる．最終的にY(m)m(5)

ここで，入力値となる {∆u m } ，Newmark のβ法から求ま&る速度・加速度成分 {u&}t , {u&d }t , {u&}t , {u&&d }t を式(5)に代入す

0.4dより，機構全体の移動量と変形量を含んだ変位量増分{∆u} が得られる．このように，機構の変位量を移動量と変形量に分離して定式化することにより，遠心力などの系の動作に伴う内力変化が考慮可能となる．

11[M ] +[C ]){∆u d } = {F }t + ∆t − {R}tβ∆t 22 β∆t11111&&{∆um}+ {u&+[M]( {u&−1){u&−1){u&d}t +(d}t −m}t +(m}t )2β2ββ∆tβ∆t2β∆t([ K ] ++ [C ](({∆u} = {∆u } + {∆u }0.2(b) Flexible modelFig.1 Target trajectory obtained by FEM

0.4

十分な精度を得るために必要な要素分割数，および質量分布に関して考察を行った．本報告では，要素両端に要素質量の半分をそれぞれ集中させる集中質量マトリク
スと，分布質量マトリクスの 2 つを用いて比較した．前節の軌道計算アルゴリズムにより算出される軌道の精度を評価するために，図 1 に示すような１リンクフレキシブルマニピュレータ（リンク長 0.4 [m]，リンク質量 0.215[kg]，曲げ剛性 EI=0.46 [Nm2]）に対し，水平面内でリンク根元を 0.2 [s]で π [rad]回転させ，その後 0.2 [s]静止させた場合の解析を実施した．時間増分は，この部材の曲げ振動固有周期の理論値が 99.8 [ms]であることを考慮し，十分に細かい 0.4 [ms]とした．また，ここでは減衰と重力は考慮しないものとする．図 1(a)には，初期値として入力した剛体リンク系の軌道を示す．また，図 1(b)には軌
道計算アルゴリズムによって出力された軌道を示すが，リンクの部材剛性が考慮された結果となっていることがわかる．図 2(a),(b)にはそれぞれ，集中質量マトリクスと分布質量マトリクスを用いて得られたリンク先端の X 座標値を示す．両図を比較すると，どちらも要素分割数が細かくなるにつれ一つの曲線に収束するが，分布質量マトリクスを用いた計算結果は，1 要素分割でも精度良く計Fig.3 1-link flexible manipulator used in
the experiments

水平面内を2.0 [s]で π [rad]回転させ，その後2.0 [s]静止させるものを与えた．
表１(a)に示す部材Aに対し，部材先端に30 [g]の質点を付加した上で, 時間増分を様々に変化させ軌道計算を行った結果を図5に示す．ここで，減衰は考慮しない．図から, 10[ms]の時間増分で振幅，周期ともに十分に収束解に達していると判断できる．実験結果との比較では，時間増分としてこの値を使用したものを用いる．
次に，表1に示す４種類の部材をリンク部材として用い，これらに実際に先のタスクを与えた．実機が描く軌道は，実機上方に設置したカメラによって撮影し，得られた画像から軌道を計測した．また，観測角度による誤差を小さくするため，最終姿勢における自由振動時の先端節点の座標値を，解析値との比較対象とした．図6(a),(b),(c),(d)に，それぞれのリンク部材の場合における先端節点X座標の解析値および実測値を示す．それぞれの図には，減衰を考慮しない場合の解析値，実験から求めた減衰定数を導入した場合の解析値，および実験値をプロットした．
全てのリンク部材に対し，減衰定数を導入した解析結果は実験結果と良く一致していることがわかる．特に，図6(b)での減衰を考慮しない解析結果は，自由振動時にほとんど振動を生じない．これは，部材Bに30 [g]の質点を付加した場合の曲げ振動固有周期が，運動の終了する2 [s]Table 1(a) Parameter of Link member A (Acrylicplastic)

Fig.4 Connected joint
算されることがわかる．さらに，収束解の曲げ振動固有周期は 99.2 [ms]となり，理論値（99.8 [ms]）と良好に一致する．以上から，今後の計算では，精度維持のためには分布質量マトリクスを用いた 1 要素分割で十分と判断し，計算コスト削減のためにこれを用いることにした．

４．実機による検証
本節では，２節で示した軌道計算アルゴリズムの有効性を示すため，実機による検証を行う．実験に用いた１
リンクフレキシブルマニピュレータの概観を図3，リンク部材とモータとの結合部を図4に示す．リンク部材の拘束点がモータの軸から大きく離れるため，数値解析では軸と拘束点との間に完全剛体に近い物性値を持つ要素を一つ設け，全体を計2要素でモデル化した．タスクとしては，
0.1ms
1ms
10ms
20ms
50ms

X(m)

0.4

0.2

Parameter
Length: L [m]Sectional area: A [m2]Mass: m [kg]Young’s modulus: E [GPa]Moment of inertia: Iy [m4]

Value
0.4
0.2×10-4
0.0952×10-1
3.2
6 667×10-12
Table 1(b) Parameter of Link member B (Acrylic plastic)
Parameter
Length: L [m]Sectional area: A [m2]Mass: m [kg]Young’s modulus: E [GPa]Moment of inertia: Iy [m4]

Value
0.4
0.4×10-4
0.1904×10-1
3.2
1.333×10-11

Table 1(c) Parameter of Link member C (Aluminum)
Parameter
Length: L [m]Sectional area: A [m2]Mass: m [kg]Young’s modulus: E [GPa]Moment of inertia: Iy [m4]

Value
0.4
0.1×10-4
0.112×10-1
69.08.333×10-13

Table 1(d) Parameter of Link member D (Aluminum)
0012time(sec)34

Fig.5 Numerical estimation on time increments

Parameter
Length: L [m]Sectional area: A [m2]Mass: m [kg]Young’s modulus: E [GPa]Moment of inertia: Iy [m4]

Value
0.4
0.2×10-4
0.224×10-1
69.0
1.333×10-11

0.5
0.4

0.4
0.3X(m)

X(m)

0.2

0

0.2
0.1
0Analysis without damping
Analysis with damping
α=0.14 β=1.0∗10 −4Experiment-0.2
012time(sec)Analysis without dampingAnalysis with damping
α=1.0 β=1.7∗10 −4
Experiment-0.1
34-0.2

00.5

0.5

0.4

0.4

0.3

0.3

0.2

0.2

0.1

3

4

0.1

00
Analysis without damping
Analysis with damping
α=0.14 β=1.0∗10 −4
Experiment

-0.1
-0.2

2time(sec)(b) Link member BX(m)X(m)(a) Link member

10

12time(sec)Analysis without damping
Analysis with damping
α=1.4 β=1.0∗10 −4
Experiment-0.1
34(c) Link member C-0.2

012time(sec)

34(d) Link member DFig.6 Comparison of trajectories for 1-link manipulator
でちょうど部材変形が無くなるように作用したからと思われるが，減衰を考慮すると振幅が大きくなり，実験結果と良好に一致するようになった．以上より，様々な剛性値を持つリンク部材による実験結果と比較しても，本軌道計算アルゴリズムはそれぞれ十分に妥当な解を算出することが確認できた．

５．結論
本研究では，様々なリンク系に対する統一的な逆動力学計算法の開発の一環として，フレキシブルマニピュレータの高精度な軌道計算法を構築した．また，本解法の有効性を実験により確認することができた．
今後は，本解法をフレキシブルマニピュレータの逆動力学計算に導入し，さらにそれを用いたフィードフォワード制御系を構築する予定である．

参考文献[1] 磯部大吾郎, 有限要素法を用いたリンク機構の逆動力学計算, 日本ロボット学会誌, 20(6), 2002, 647-653.
[2] D.Isobe, D.Imaizumi and A.Yagi, A Parallel Solution
Scheme for Inverse Dynamics of Link Mechanisms, JSMEInternational Journal Series C, 46(2), 2003, 766-771.

[3] D. Isobe, Y. Chikugo, D. Imaizumi, S. Sato and A. Yagi,
Feed-forward control of link mechanisms under various
boundary conditions by using a parallel solution scheme, Proc.
2003 IEEE/RSJ International Conference on Intelligent Robotsand Systems (IROS 2003), 2003, 2773-2778, Las Vegas U.S.A.
[4] E.Boya, A Finite-Element Approach to Control theend-Point Motion of a Single-Link Flexible Robot, J. ofRobotic Systems, 4(1), 1987, 63-75.
[5] A.A.Shabana, Dynamics of Flexible Bodies Using
Generalized Newton-Euler Equations, Trans. of ASME, J. ofDynamic Systems, Measurement, and Control, 112, 1990,496-503.
[6] H.Asada, Z. -D.Ma and H.Tokumaru, Inverse Dynamics ofFlexible Robot Arms: Modeling and Computation for
Trajectory Control, Trans. of ASME, J. of Dynamics Systems,Measurement, and Control, 112, 1990, 177-185.
[7] 磯部大吾郎，今泉大作，リンク系の部材剛性に依存しない統一的逆動力学計算法，
日本機械学会論文集（C 編），70(691), 2004, 728-735.