二重振り子の運動シミュレーション（粘性抵抗も考えるよ！）

卒論も終わり一息ついたので、ラグランジュの運動方程式の復習として二重振り子の運動方程式の導出と、導出した運動方程式を用いた運動シミュレーションをやってみた。

運動方程式の導出

上図のようにジョイント1、2の回転角をそれぞれ $\theta_1, \theta_2$ 、おもり1、2の質量をそれぞれ $m_1, m_2$ 、ロッド1、2の長さをそれぞれ $l_1, l_2$ とする。このとき、おもり1、2の位置ベクトルはそれぞれ次のようになる。

$\mathbf{r}_1 = \begin{bmatrix} l_1 \sin \theta_1 \\ - l_1 \cos \theta_1 \end{bmatrix}, \quad \mathbf{r}_2 = \begin{bmatrix} l_1 \sin \theta_1 + l_2 \sin \theta_2 \\ - l_1 \cos \theta_1 - l_2 \cos \theta_2 \end{bmatrix}$

これらを時間で微分して各おもりの速度を求める。

$\dot{\mathbf{r}}_1 = \begin{bmatrix} l_1 \dot{\theta}_1 \cos \theta_1 \\ l_1 \dot{\theta}_1 \sin \theta_1 \end{bmatrix}, \quad \dot{\mathbf{r}}_2 = \begin{bmatrix} l_1 \dot{\theta}_1 \cos \theta_1 + l_2 \dot{\theta}_2 \cos \theta_2 \\ l_1 \dot{\theta}_1 \sin \theta_1 + l_2 \dot{\theta}_2 \sin \theta_2 \end{bmatrix}$

二重振り子の運動エネルギー $K$ は次式のようになる。

$\begin{align} K &= \frac{1}{2} m_1 \|\dot{\mathbf{r}}_1\|^2 + \frac{1}{2} m_2 \|\dot{\mathbf{r}}_2\|^2 \\ &= \frac{1}{2} (m_1 + m_2) l_1^2 \dot{\theta}_1^2 + m_2 l_1 l_2 \dot{\theta}_1 \dot{\theta}_2 \cos (\theta_2 - \theta_1) + \frac{1}{2} m_2 l_2^2 \dot{\theta}_2^2 \\ \end{align}$

途中で加法定理を使った。重力加速度を $g$ とおくとポテンシャルエネルギー $U$ は次のようになる。

$U = -(m_1 + m_2) g l_1 \cos \theta_1 - m_2 g l_2 \cos \theta_2$

これらより、ラグランジアン $L$ は次のように求められる。

$\begin{align} L &= K - U \\ &= \frac{1}{2} (m_1 + m_2) l_1^2 \dot{\theta}_1^2 + m_2 l_1 l_2 \dot{\theta}_1 \dot{\theta}_2 \cos (\theta_2 - \theta_1) + \frac{1}{2} m_2 l_2^2 \dot{\theta}_2^2 \\ &\quad + (m_1 + m_2) g l_1 \cos \theta_1 + m_2 g l_2 \cos \theta_2 \end{align}$

（この辺りから計算が嫌になってくる）

一般化座標を $\theta_1, \theta_2$ とすると、ラグランジュの運動方程式は次の二式となる。

$\begin{align} \frac{d}{dt} \left( \frac{\partial L}{\partial \dot{\theta}_1} \right) - \frac{\partial L}{\partial \theta_1} &= u_{\theta_1} \\ \frac{d}{dt} \left( \frac{\partial L}{\partial \dot{\theta}_2} \right) - \frac{\partial L}{\partial \theta_2} &= u_{\theta_2} \end{align}$

ただし、 $u_{\theta_1}, u_{\theta_2}$ はそれぞれ $\theta_1, \theta_2$ に対する一般化力である。一般化力はあとで考えるから、とりあえずこのようにおいておく。ラグランジュの運動方程式に上で求めた $L$ を代入してゴリゴリ計算していく。

$\theta_1$ について：

$\begin{align} &\frac{\partial L}{\partial \dot{\theta}_1} = (m_1 + m_2) l_1^2 \dot{\theta}_1 + m_2 l_1 l_2 \dot{\theta}_2 \cos (\theta_2 - \theta_1) \\ \Rightarrow &\frac{d}{dt} \left( \frac{\partial L}{\partial \dot{\theta}_1} \right) = (m_1 + m_2 l_1^2) \ddot{\theta}_1 + m_2 l_1 l_2 \left[ \ddot{\theta}_2 \cos (\theta_2 - \theta_1) - \dot{\theta}_2 (\dot{\theta}_2 - \dot{\theta}_1) \sin (\theta_2 - \theta_1) \right] \\ &\frac{\partial L}{\partial \theta_1} = m_2 l_1 l_2 \dot{\theta}_1 \dot{\theta}_2 \sin (\theta_2 - \theta_1) - (m_1 + m_2) g l_1 \sin \theta_1 \end{align}$

$\theta_2$ について：

$\begin{align} &\frac{\partial L}{\partial \dot{\theta}_2} = m_2 l_1 l_2 \dot{\theta}_1 \cos (\theta_2 - \theta_1) + m_2 l_2^2 \dot{\theta}_2 \\ \Rightarrow &\frac{d}{dt} \left( \frac{\partial L}{\partial \dot{\theta}_2} \right) = m_2 l_1 l_2 \left[ \ddot{\theta}_1 \cos (\theta_2 - \theta_1) - \dot{\theta}_1 (\dot{\theta}_2 - \dot{\theta}_1) \sin (\theta_2 - \theta_1) \right] + m_2 l_2^2 \ddot{\theta}_2 \\ &\frac{\partial L}{\partial \theta_2} = - m_2 l_1 l_2 \dot{\theta}_1 \dot{\theta}_2 \sin (\theta_2 - \theta_1) - m_2 g l_2 \sin \theta_2 \end{align}$

（計算ミスがこわい）

以上から、運動方程式は次のように得られる。

二重振り子の運動方程式

$\begin{align} (m_1 + m_2) l_1^2 \ddot{\theta}_1 + m_2 l_1 l_2 \left[ \ddot{\theta}_2 \cos (\theta_2 - \theta_1) - \dot{\theta}_2^2 \sin (\theta_2 - \theta_1) \right] + (m_1 + m_2) g l_1 \sin \theta_1 &= u_{\theta_1} \\ m_2 l_1 l_2 \left[ \ddot{\theta}_1 \cos (\theta_2 - \theta_1) + \dot{\theta}_1^2 \sin (\theta_2 - \theta_1) \right] + m_2 l_2^2 \ddot{\theta}_2 + m_2 g l_2 \sin \theta_2 &= u_{\theta_2} \end{align}$

もう少し見やすくするために $\Delta \theta := \theta_2 - \theta_1$ として次のように行列形式に整理しておく。

二重振り子の運動方程式（行列形式）

$\begin{bmatrix} (m_1 + m_2) l_1^2 & m_2 l_1 l_2 \cos \Delta \theta \\ m_2 l_1 l_2 \cos \Delta \theta & m_2 l_2^2 \end{bmatrix} \begin{bmatrix} \ddot{\theta}_1 \\ \ddot{\theta}_2 \end{bmatrix} = \begin{bmatrix} u_{\theta_1} + m_2 l_1 l_2 \dot{\theta}_2^2 \sin \Delta \theta - (m_1 + m_2) g l_1 \sin \theta_1 \\ u_{\theta_2} - m_2 l_1 l_2 \dot{\theta}_1^2 \sin \Delta \theta - m_2 g l_2 \sin \theta_2 \end{bmatrix}$

左辺の行列を $\mathbf{M}$ 、右辺を $\mathbf{F}$ と書くことにすると、運動方程式は次のように角加速度について解ける( $\mathbf{M}$ は常に正則)。

$\begin{bmatrix} \ddot{\theta}_1 \\ \ddot{\theta}_2 \end{bmatrix} = \mathbf{M}^{-1} \mathbf{F}$

これを4次のルンゲ=クッタ法で時間発展させることで二重振り子の運動をシミュレーションする！

粘性抵抗について

ジョイント1およびジョイント2には粘性抵抗(ジョイントの摺動速度に比例する抵抗)がはたらくとする。粘性抵抗はこの場合トルクである。ジョイント1、2の粘性抵抗は比例係数をそれぞれ $c_1, c_2$ として次のように表せる。

$\tau_1 = -c_1 \dot{\theta}_1, \quad \tau_2 = -c_2 (\dot{\theta}_2 - \dot{\theta}_1)$

粘性抵抗はジョイントの回転速度に比例するのではなく、摺動速度(ジョイントの軸受けがすべって回転する速度)に比例するから、ジョイント2の抵抗を $\tau_2 = -c_2 \dot{\theta_2}$ とするのは不適切。

粘性抵抗は保存力ではないから、この抵抗は一般化力の項に入ってくる。問題はどのように入ってくるか。

これを考えるには二重振り子の状態が微小変化したときに抵抗がする微小仕事を考えればいい。抵抗による微小仕事は次のように与えられる。

$\begin{align} \delta^{\prime} W &= \tau_1 d\theta_1 + \tau_2 d(\theta_2 - \theta_1) \\ &= - c_1 \dot{\theta}_1 d\theta_1 - c_2 (\dot{\theta}_2 - \dot{\theta}_1) d(\theta_2 - \theta_1) \\ &= \left[ - c_1 \dot{\theta}_1 + c_2 (\dot{\theta}_2 - \dot{\theta}_1) \right] d\theta_1 - c_2 (\dot{\theta}_2 - \dot{\theta}_1) d\theta_2 \end{align}$

一方で、一般化力 $u_{\theta_1}, u_{\theta_2}$ は次で定義される。

$\delta^\prime W = u_{\theta_1} d\theta_1 + u_{\theta_2} d\theta_2$

よって、粘性抵抗による一般化力は次のようになる。

粘性抵抗による一般化力

$\begin{align} u_{\theta_1} &= - c_1 \dot{\theta}_1 + c_2 (\dot{\theta}_2 - \dot{\theta}_1) \\ u_{\theta_2} &= - c_2 (\dot{\theta}_2 - \dot{\theta}_1) \end{align}$

運動シミュレーション

二重振り子の各パラメータは下表のように設定した。

パラメータ	値	単位
(m₁, m₂)	(1, 1)	kg
(l₁, l₂)	(0.6, 0.6)	m

また、重力加速度はg = 9.80665 [m/s²]とし、シミュレーション刻みは0.001[s]とした。

計算はMATLABで行った。

粘性抵抗なし

まずは、粘性抵抗なしの場合(c₁ = c₂ = 0 [kgm²/s])についてシミュレーションをした。

初期条件を $\theta_1 = \theta_2 = \frac{2}{3}\pi \ \mathrm{[rad]}$ 、 $\dot{\theta}_1 = \dot{\theta}_2 = 0\ \mathrm{[rad/s]}$ とした場合：

リハビリとして2重振子の運動方程式の導出と，そのシミュレーションをやってみた pic.twitter.com/kLyZuAvuui
— スコル (@__scol__p) February 24, 2019

いわゆる「カオス運動」をしているのが見てとれる。振り子を二つ繋げただけなのにこんな複雑な運動をするのは相変わらずすごい。

粘性抵抗がない場合は保存系となるから力学的エネルギーは保存されるはず。そのため、力学的エネルギーの推移を見れば導いた運動方程式や計算結果の正しさをある程度判断することができる。さて、どうなっているか・・・

f:id:scol:20190227175848p:plain — 力学的エネルギーの推移

力学的エネルギーの変動は非常に小さいから正しいといえるだろう。図見ると、「大きく変動して少し下がる」というのを繰り返しているのが分かる。二重振り子が急激な運動をしているときに力学的エネルギーが大きく変動している。

今回はシミュレーション刻みを一定としたが、二重振り子の運動方程式ようにゆっくりとした変動と激しい変動がともに出てくる微分方程式(硬い方程式というらしい)の時間発展行う場合は、刻み幅を制御しながら計算する方法を用いるのがより適切らしい。埋め込み型ルンゲ=クッタ法というものを使えば誤差の評価と刻み幅の制御が行えるそうだ。

ルンゲ=クッタ法の説明と刻み幅制御 - シキノート

また、保存量をきちんと保存しながら数値積分を行える方法もあり、これをシンプレクティック数値積分というらしい。

シンプレクティック法をシミュレーションで理解する - ハローラスク

ところで、大角度ならいつでもカオスになると思っていたけど、カオスにならない場合もあった。それが次の場合である。

初期条件を $\theta_1 = \frac{2}{3}\pi \ \mathrm{[rad]}$ 、 $\theta_2 = \frac{1}{2}\pi \ \mathrm{[rad]}$ 、 $\dot{\theta}_1 = \dot{\theta}_2 = 0\ \mathrm{[rad/s]}$ とした場合：

カオス運動をすることで有名な2重振子だが，初期角度を120°，90°にしたら周期的な運動が出てきた．大角度でもカオスにならないことはあるんだな......
カオスにならないための条件って求められてたりするのかな？ pic.twitter.com/mtm9zGjDsz
— スコル (@__scol__p) February 24, 2019

$\theta_1$ と $\theta_2$ が同時に0になるのがカオスにならないポイントなのかなと思った。

粘性抵抗あり

ここからは、すべて初期条件を $\theta_1 = \theta_2 = \frac{2}{3}\pi \ \mathrm{[rad]}$ 、 $\dot{\theta}_1 = \dot{\theta}_2 = 0\ \mathrm{[rad/s]}$ とした。

まずは、2つのジョイントに弱い抵抗がある場合(c₁ = 0.04 [kgm²/s]、c₂ = 0.02 [kgm²/s])：

2つのジョイントに粘性抵抗(角速度に比例する抵抗)がある場合の2重振子
抵抗がある方が自然な感じがする pic.twitter.com/cjieQSPg1O
— スコル (@__scol__p) February 25, 2019

運動は次第に穏やかになり、非線形項の効果が小さくなって最終的には周期的な運動に落ち着いている。今回のシミュレーション条件では、微小かつ周期的な振動をしているとき、その振幅比は $\theta_1 : \theta_2 = 1 : \sqrt{2}$ か $\theta_1 : \theta_2 = 1 : -\sqrt{2}$ の2パターンであることが微小振動時の線形化運動方程式を調べることによりわかる(マイナスは逆向きの振動をするという意味)。

抵抗があるから、このときはもちろん力学的エネルギーは保存されない。

次は、ジョイント1には強い抵抗がありジョイント2には抵抗がない場合(c₁ = 5 [kgm²/s]、c₂ = 0 [kgm²/s])：

上のジョイント：強い抵抗
下のジョイント：抵抗なし

下の振子の振動が上の振子に吸収されている pic.twitter.com/Hvy2Tzs0u0
— スコル (@__scol__p) February 25, 2019

上の振り子が動吸振器の役割を果たし、振動は急速に減衰している。動吸振器というのは、一言でいえば「振動を吸収する装置」。地震の振動を抑えるために高いビルの屋上についている可動式のおもりは動吸振器の一種である。

最後は、ジョイント1には抵抗がなくジョイント2には強い抵抗がある場合(c₁ = 0 [kgm²/s]、c₂ = 5 [kgm²/s])：

上のジョイント：抵抗なし
下のジョイント：強い抵抗

そりゃそうだろうなという結果 pic.twitter.com/8HS3ccoRsw
— スコル (@__scol__p) February 25, 2019

そりゃそうだろなという結果。ジョイント2に強い抵抗があればそりゃ単振り子みたいな運動をするだろう。さっきと違って、この場合は振動はほとんど減衰されない。

まとめ

二重振り子の運動シミュレーションを通じてラグランジュの運動方程式を一般化力を含めて軽く復習できた。実際にシミュレーションを行って振り子がひょこひょこ動いているのを見ると何だか愛着が湧いてくるね。

シミュレーションプログラム(MATLAB)

github.com

間違いなどありましたらご指摘ください。

スコルの知恵袋

主にプログラミング関係の気になったこと等をまとめたブログ