有限要素法によるロボット機構の歩行動作における衝撃力の予測
筑波大学○磯部大吾郎

小澤一裕

Prediction of Impulsive Force in Walking Operation of Robotic
Architecture by Using Finite Element Method
Univ. of Tsukuba ○Daigoro ISOBE

Kazuhiro KOZAWA

Abstract: In general, impulsive force produced in constituting members of a robotic architecture was
not a main issue. But as the robotic architecture becomes maneuverable in high-speed motion, thestress wave produced in the structure is due to increase, as well as the danger of mechanicaldestruction. In this research, an algorithm calculating inertial force distributed in multi-linkmechanisms is constructed by using Finite Element Method. It is verified by carrying out some simple
numerical tests, along with a simulation of a biped. robotic model. By using the proposed scheme, theprocedure of shifting mass, for example, can be designed suitably to reduce the impulsive force
occurring in the mechanisms.
Key Words: Impulsive Force, Stress Propagation, Robotic Architecture, Finite Element Method,Biped. Robot

1.はじめに1.はじめに構造力学的見地からロボット機構を観察すると，これらはやはり，モータや部材を介して力を伝達すわず，機構内の発生力をある程度予測し，フィードフォワード的に制御する必要性が増すことが考えられる．

る 構造物 の一種である．介護，あるいはエンタ

そこで本研究では，リンク機構が目標軌道上を動ーテイメントロボットなどの人間共存型ロボット作する際に，機構内に発生する衝撃力を短時間で計に注目が集まる中，未だにロボット機構は完全な剛算する手法を開発することにした．そのアルゴリズ体として設計・製作され，剛体の固まりとして制御ムには，構造解析などのツールとして用いられる有されている感が否めない．人間と接触した場合に危限要素法(FEM)を適用した．また，リンク機構のモ険を伴うだけではなく，機構内の応力集中や衝撃波デル化には，骨組構造の有限要素解析でその有効性の伝播によって機構自体が破損する可能性が存在が立証されている Shifted Integration 法（SI 法）する（Fig.1 参照）．また，特に 2 足歩行ロボットなどでは，足の裏に作用する反力をセンサで検知し，これから機構全体の ZMP を求め，転倒しないように制御することが重要な課題の一つとなっている[1]．
このように，機構内に発生する力を知ることは，無駄のない機構設計や無理のない動作計画を行う上で，必要不可欠な要素となることが多い．
他方，多くのロボットは近い将来，高速で移動することを目標としており，その機構内に生じる衝撃力が増し，構成部材や部品の破壊の危険性が高くなると考えられる．また，動作が高速になるとセンサからの情報をフィードバックするだけでは間に合Fig.1 Impulsive force causing damage in robots

Joint (Servo motor)+Center ofgravity

Nodal point

Numerical integration
point

Link
Finiteelement+(a)Link mechanism

Fig.2 Linear Timoshenko beam element and(b)Finite elements

Fig.3
Fig.3 Modeling of link mechanism byShifted Integration technique

its physical equivalentを用い，要素数を削減して計算時間の短縮化を図っ

Fig.3 に，SI 法によるリンク機構のモデル化概念

た．FEM を用いたリンク機構の内力分布解析アル図を示す．図に示すように，モータとリンク部材にゴリズムを構築し，その動作確認を行った上で，歩よって構成されたリンク機構は，リンクの重心位置行動作を行う 2 足歩行ロボットモデルの衝撃力予に節点を設け，関節位置の反対側に数値積分点をシ測解析を実施し，重心移動が機構内に発生する衝撃

フトした２つの線形チモシェンコはり要素によっ力に与える影響について有益な情報を得た．て表現される[4]．

2. SI 法によるリンク機構のモデル化
本研究でリンク機構のモデル化に導入する SI 法
は，本来は骨組構造の有限要素解析に用いられる手Start
Input material constants andmodel data法である．Fig.2 に示す線形チモシェンコはり要素

istage <= nstageと，ヒンジが回転ばねおよびせん断ばねによって陽Read model
coodinates

noyesに表現される剛体ばねモデル(RBSM)のひずみエ

Make transformation matrix [T]ネルギ近似式の考察により，次式のような数値積分istage = 1

点とばねの位置関係が導出されている[2]．

Endnoyes

Make mass matrix [M]s1 = Är1 または r1 = Äs1(1)istage >= 2noyes

ここに，s1 および r1 はそれぞれ，数値積分点位置おJudge ground contact

よび形成されるべき塑性ヒンジ（またはばね）の位

Change bindingcondition

置である．
SI 法では，上式にしたがってヒンジが生じうる点と対称の位置に数値積分点をあらかじめ配置することにより，精度良く骨組構造物の崩壊荷重解を

yes

Remake mass matrix [M]Make stiffness matrix [K]Calculate simultaneous equations得ることができる．また，本手法を順応型に発展させることにより，様々な条件下において最小限の計

Output normalforce

算コストで精度の高い解が得られることが立証されている[3]．

Fig.4 Flow chartno3. FEM を用いた内力分布解析アルゴリズム

アルゴリズムでは，最後に，式(3)にしたがって算本研究で構築したアルゴリズムのフローチャー

出されたひずみ情報から部材内の軸力を算出する．

トを Fig.4 に示す．本アルゴリズムでは，動的変形を伴うモデルデータ（すなわち，目標軌道）を逐次4. 数値例入力し，ステップ毎には静的な解析を行い，自重に本章では，有限要素法による内力分布解析アルゴ伴う変形量を算出する．また，時間積分法の一種でリズムを用いていくつかの解析を行い，本手法の有ある Newmark のβ法（β=1/4）を用い，次式で表される速度・加速度情報をステップ毎に算出する．

Äf°ugt+Å t =311fÅug Ä (Ä 1)fug_ t2åÅt
2å
1 Ä 4å
f°ugt Åt4å

(2a)111fÅugÄfug
_ t Ä( Ä1)f°ugt
åÅt2åÅt2å(2b)Normal Force [N]fug_ t+Å t =[s]
T=1.0
T=3.0
T=5.0
T=10.0210

これらの値を逐次的に加算し，次式の運動方程式に代入することにより，機構が能動的に変形する際の-1

内力分布が算出される．

011([K] +[M ])fÅug =[C] +2åÅtåÅt2Fig.5(a) Normal force calculated from
orbit
orbit data11Ä 1)fug_ tÄ(Ä 1)f°ugt Åt)2å
4å

+[M ](11fug_ t+(Ä 1)f°ugt )åÅt2å

10Time [s]fF gt+Å t Ä fRgt+[C]((5(3)

外力ベクトル，{R}：内力ベクトル，Δt：時間増分である．今後，衝撃力予測システムを実機に搭載する際には減衰マトリックスの同定が必要となるが，現段階では，簡単のためにそれを無視した．

Normal Force [N]

ここで，[M]：全体質量マトリックス，[C]：全体減衰マトリックス，[K]：全体剛性マトリックス，{F}：

T=1.0
T=3.0
T=5.0
T=10.0[s]32
10

接地・離地を伴う運動を行う機構の場合には，各節点の鉛直方向座標を随時監視し，その接地・離地-1の判定後に，鉛直方向自由度を拘束または自由にする．またこの時，自由度番号の更新を行って質量マトリックスなどの情報を組み替える操作を行う．本

0510Time [s]Fig.5(b) Normal force calculated by FEM

1[m]

①②④

③⑤Fig.6 Trajectory for a 33-link mechanism
効性について検証した．

2.98[m]となりモデルの両端間距離とほぼ一致した．

4.1 軸力の検証

すなわち，機構内の応力波伝播の状態が正確に計算部材内に発生する軸力の計算結果を検証するたされていることが判明した．

め，加速・減速によって周回する円運動軌道を 1 リンク機構に与え，目標軌道データから求めた結果と本手法による結果とを比較した．Fig.5(a)に目標軌和の時刻歴，Fig.5(b)に FEM により算出した軸力の時刻歴を示す．運動周期が長い場合には両者の結果は一致しているが，運動周期が短くなるにつれて相違が生じる．これは，前者の計算にはコリオリ力のNormal Force [N]

道データから算出した遠心力と自重の軸方向成分0①-10影響が含まれていないためである．一方，FEM によに若干の振れが見られ，コリオリ力の影響が生じている様子が確認できる．

4.2 衝撃力の検証②-20る算出結果では，T=1.0[s]の場合に回転運動の最後⑤④③00.5

1Time [s]Fig.7 Impulsive force propagated in thelink member

次に，接地・離地する際に機構内に生じる衝撃力の計算結果を検証するため，Fig.6 のような 3 リン

-15.8

ク機構に接地・離地を繰り返す運動軌道を与えた．
破線の○で囲ったリンク部材に生じる軸力の時刻ず，他方の脚が離地する際にも相当量の衝撃力が発生していることが分かる．これは，機構の全重量が瞬時にこのリンク部材に載荷されることに起因しNormal Force [N]

歴を Fig.7 に示す．その部材が接地する時のみなら-16-16.2-16.4ている．さらに，Fig.7 の破線の○で囲った部分を拡大した Fig.8 を見ると，振動の半周期が求められる．この値と，部材のヤング率と密度から算定した応力伝播速度を用いて応力伝播距離を求めると，約-16.60.004

0.006

0.008

Tim e [s ]Fig.8 Cycle of stress propagation

ZZZ生じる振動であると考えられる．また，最後に着地する際に左脚に大きな振動が生じているが，重心をYXYXYX

Z適切に移動させることにより，この時に生じる衝撃力を最大約 40％減少できることが Fig.10, Fig.11の比較から分かった．また，今度は右脚に着目する

ZZと，自らが接地する場合と同様に，右脚のみが機構を支える場合にも大きな振動が続けて生じている
YYX

X

YX

ことが確認できる．これは，全重量が瞬時に右脚に掛かったことによる振動と，右脚の着地時に生じたZ衝撃が減衰する前に次の運動に移ったために生じた振動によるものであり，右脚を構成する部材にはYX相当な負荷が掛かっているものと思われる．またこの際の衝撃力も，Fig.9 に示すような重心移動を行

Fig.9 Trajectory for a biped. mechanism

うことにより，重心移動しない場合に比べ約 40％の緩和が可能であることを確認できた．

5. 歩行ロボットにおける衝撃力の予測

本解析は CPU: Intel Pentium Ⅲ600[MHz]，RAM:

本章では，第３章に示したアルゴリズムを用い，384[MB]を有する PC 上で実行した．１ステップ当り歩行動作を行う 5 リンク機構の衝撃力予測解析をの時間増分を 0.005[s]，合計 200 ステップで解析行った．また，衝撃を緩和する方法として重心移動を行った場合に要した計算時間は 0.35[s]であり，によるものを取り上げ，重心を移動する場合としな現時点においても十分リアルタイム処理に対応でい場合とを比較した．さらに，計算コストの面からきる結果である．

本手法の有効性を検証した．
Fig.9 に，2 足歩行ロボットを模擬した 5 リンク

6.最後に機構に与えた目標軌道を示す．紙面に向かって右側

本研究では，能動的に変形する機構を Shiftedを左脚，左側を右脚とし，両脚を連結する部分をボIntegration 法を用いてモデル化し，FEM による内ディ部と称することにする．脚は，長さ 0.3[ｍ]の力分布解析アルゴリズムを構築した．歩行動作を行部材が 2 本結合しており，ボディ部は，0.4[ｍ]のうモデルに対して解析を実施した結果，機構内に生部材 1 本から構成されている．また，ボディ部中央じる軸力の時刻歴を短い計算時間で得ることがでに 2.0[kg]の質量を配置し，この質量を図のようにきた．接地・離地の際に生じる衝撃力が十分な精度移動させることにより衝撃緩和を図った解析を行で求められ，さらに，適切なタイミングによる重心う．本例題では，ロボットが 1.5[s]で 2 歩前進す移動によって衝撃力が大幅に緩和できることなど，るように設定してあり，両脚が接地し次の動作に移有益な情報が得られることを確認した．長距離を速る前には，0.1[s]弱の静止時間を設けてある．

く歩くような場合，例えば競歩のフォームなどのよFig.10(a)，(b)には，重心を移動しない場合のそうに，上半身をひねって着地する脚と反対側の脚のれぞれの脚に作用する軸力の時刻歴を示す．また，方へ重心を交互に移動する方法は，衝撃力を緩和す

Fig.11(a)，(b)には，重心を移動した場合の時刻歴る上でも非常に合理的な歩行法である．

を示す．左脚に着目すると，反対の右脚を上げた瞬提案した手法は，リアルタイム処理も可能である間に衝撃力が生じていることが分かる．これは，両ため，今後はロボット機構の設計時およびオフライ脚で支えていた自重が瞬時に左脚に掛かるためにンによる事前の衝撃予測のみならず，実機に搭載し0-10-20

element 1-30

Normal Force [N]Normal Force [N]0-10-20-30
eleement 1000.501Time [s]0.51Time [s](a) Left leg(b) Right leg

Fig.10 Impulsive force propagated in each leg

0Normal Force [N]Normal Force [N]0-10-20-30-10-20-30element 10element 100.5

100.5Time [s]Time [s](a) Left leg(b) Right leg

1Fig.11 Influence of shifting mass on propagation of impulsive forceてオンライン制御に適用していく予定である．そのSystems (IROS’01), pp1093-1099, 2001.

際には，センサからの情報に頼る従来の制御とは違[2] 都井 裕, 骨組構造および回転対称シェル構造の有限い，予測情報に基づいた速い制御が可能となること

要素解析におけるShifted Integration 法について, 日本造を期待している．今後，ロボットの動作が高速化す

船学会論文集, 168, pp.357-369, 1990.

るにつれ，ロボットおよびその周辺の安全性を考え，[3] Isobe, D. and Toi, Y., Analysis of Structurally

構造力学的な見地からの設計および動作計画が練Discontinuous Reinforced Concrete Building Frames Using theられる必要性があると考える．

ASI Technique, Computers and Structures, 76‑4,pp.471‑481,762000.

参考文献[1] Kanehiro, F. et al., Virtual Humanoid Robot Platform toDevelop Controllers of Real Humanoid Robots without Porting,Proc. of 2001 IEEE/RSJ Int. Conf. on Intelligent Robots and[4] Isobe, D., A Unified Numerical Scheme for Calculating
Inverse Dynamics of Open/Closed Link Mechanisms , Proc. ofthe 27th Annual Conf. of the IEEE Industrial Electronics
Society (IECON'01), pp.341-344, 2001.