テストデータ作成ツール

投稿日: 2014年4月28日 作成者: kiidax

マウスで書いた線をOctave用のデータにできます。IE8は未対応です。

20140428-01

カテゴリー: Blender | コメントをどうぞ

ピョンピョンはねないで着地したい数値解析

投稿日: 2014年4月13日 作成者: kiidax

GNU Octaveを使って数値解析をすることにした。

6つのグラフはばね定数kを変えてプロットしている。それぞれのグラフで、赤、緑、青、マゼンタ、シアンの順番にダンパー抵抗を１から200まで指数的に増やしている。横軸は時間で縦軸は車体の高さだ。初期位置は50cmの高さで、この時に10m/sで下に向かって落ちている。5m程度から落ちたケースとしている。これを見る限りだと調べている範囲では深く沈みこまないかどうかはほぼ、ダンパー抵抗で決まるようだ。5mから落ちているので若干沈み込むぐらいがいいだろう。k=50で、hは青の真ん中、h=14を選ぶことにした。

20140413-01

次に初速を0から20m/sまで変えてテストしてみる。5m/sを超えるのは厳しそうだ。

20140413-02

ここで、当初はk=100くらいを想定していたのだが、大分kが小さくなった。kが小さくなったことで、定常状態の位置がばねの自然長より短くなる。フックの法則より、mg/k ~= 0.20となり、ばねの自然長を余計に伸ばしておく必要がある。もともと自然長は1.2mだったため1.4mに変更。

suspension_test.m

function suspension_test()

  % Simulate spring and damper.

  spring_value = logspace(log10(10),log10(200),6);
  damper_value = logspace(log10(1),log10(200),15);

  clf;
  for i = 1:length(spring_value)
    subplot(2,3,i);
    hold on;
    axis([ 0, 5.0, -1, 1 ]);
    xlabel('Time');
    ylabel('Position');
    title(sprintf('k=%dm, h=1m,...,200m', spring_value(i)));
    for j = 1:length(damper_value)
      [t, x] = spring(spring_value(i), damper_value(j));
      plot(t, x, num2str(floor((j+2)/3)));
    end
    hold off;
  end
end

function [t, x] = spring(spring_value, damper_value)

  m = 1.0;
  t = linspace(0,5,300);
  x0 = [0.5; -10];

  A = [ 0 1 ; -spring_value / m -damper_value / m ];

  [temp, istate, msg] = lsode(@(x)A*x, x0, t);
  x = temp(:,1)';
end

suspension_test2.m

function suspension_test2()

  % Simulate spring and damper.

  spring_value = logspace(log10(10),log10(200),6);
  damper_value = 14;
  initial_speed = linspace(0,20,15);

  clf;
  for i = 1:length(spring_value)
    subplot(2,3,i);
    hold on;
    axis([ 0, 5.0, -1, 1 ]);
    xlabel('Time');
    ylabel('Position');
    title(sprintf('k=%dm, v_0=0,...,20', spring_value(i)));
    for j = 1:length(initial_speed)
      [t, x] = spring(spring_value(i), damper_value, initial_speed(j));
      plot(t, x, num2str(floor((j+2)/3)));
    end
    hold off;
  end
end

function [t, x] = spring(spring_value, damper_value, initial_speed)

  m = 1.0;
  t = linspace(0,5,300);
  x0 = [0.5; -initial_speed];

  A = [ 0 1 ; -spring_value / m -damper_value / m ];

  [temp, istate, msg] = lsode(@(x)A*x, x0, t);
  x = temp(:,1)';
end
カテゴリー: Blender | コメントをどうぞ

ピョンピョンはねないで着地したいが失敗

投稿日: 2014年4月12日 作成者: kiidax

20140412-01

レイキャストを使った２輪車は、車体から地面への距離を測って、距離と速度に応じて反発力をPythonで計算する。ソースコードはこんな感じだ。

import bge
import math
import mathutils

suspensionLength = 1.2
downVect = mathutils.Vector((0, 0, -1))
springValue = 90.0
damperValue = 5.0

def update(cont):
    own = cont.owner
    # Check if the foot touch on the ground.
    obj, point, normal = own.rayCast(own.worldPosition + downVect, None, suspensionLength * 1.2, "", 0, 0, 0)
    if obj:
        velocityVect = own.getLinearVelocity()
        dist = own.getDistanceTo(point)
        position = dist - suspensionLength
        force = - position * springValue - velocityVect.z * damperValue
        own.applyForce([0, 0, force])

springValueとdamperValueはそれぞればねとダンパーに関するパラメーターだが、この値によって、車体がピョンピョン跳ねるか滑らかに安定するかが変わってくる。

上記は大まかには、重力とばねの力が安定する場所（地面ではない）を原点、springValueをk、damperValueをhとすれば、以下のような微分方程式でオブジェクトを動かしていることになる。

\displaystyle m \ddot x = - k x - h \dot x

２階線形微分方程式 （ 同次 ）の解法より、特性方程式の解によってxの解が変わって車体の動き方が変わる。特性方程式は

\displaystyle m \lambda ^2 + h \lambda + k = 0

で

\displaystyle \lambda = \frac {-h \pm \sqrt { h^2 - 4mk } } {2m}

まず、h^2 - 4mk 0の場合は、判定式の解を\lambda_1, \lambda_2と置くと、

\displaystyle x = C_1 e^{\lambda_1 t} + C_2 e^{\lambda_2 t}

と書けるが、hが正であることから、{\lambda}_1, {\lambda}_2は共に負で、ピョンピョンはねることなく次第に止まっていく。

1mの高さから落ちたとして3秒後には振幅が3cm程度に軽減したいので、hを最初のケースと同じ程度に設定しておけばいいだろう。2.3にしておけばいいと思ったが、勘違いをしていることに気付いた。微分方程式では、車体が地面からどんなに離れていても力が働いていることになるが、現実にはそうではない。車輪が接地した瞬間にすでに速度がある程度ある状況から始めなければならない。

ちょっと、代数的に解析するのがめんどくさそうなので、いったん数値的に解析してみようと思う。

ところで、微分方程式は

\displaystyle \binom{\dot x}{\ddot x} = A \binom{x}{\dot x} \\ A = \begin{pmatrix} 0 & 1 \\ -k/m & -h/m \\ \end{pmatrix}

と状態方程式であらわされる。はじめての現代制御理論 講義09 状態フィードバックと極配置によるとシステムが安定するにはAの固有値の実数部が負になればよい。

\displaystyle \det (\lambda E - A) = \begin{vmatrix} \lambda & -1 \\ k/m & \lambda + h/m \\ \end{vmatrix} = 0

これは微分方程式の特性方程式と同じになる。この場合常に解の実数部は負になるので常に安定する。

カテゴリー: Blender | コメントをどうぞ

Blenderを触っていないような気がする

投稿日: 2014年3月15日 作成者: kiidax

20140315-01

Shape Keyで目を閉じれるようにしているけど、単に目を閉じても、よくあるマンガ的に目を閉じて笑う表情と違うな。

カテゴリー: Blender | コメントをどうぞ

並進運動する剛体の角運動方程式と非慣性系

投稿日: 2014年3月14日 作成者: kiidax

さて、角運動方程式を学んだところで倒立振子（とうりつしんし）の制御の説明を探してみて、気づいたのは多くのサイトで角運動の方程式を建てる場合に暗黙の内に非慣性系での運動方程式を立てていることだ。なんとなく、そういうものかとスルーしたくなるが、ちょっと調べてみた。

倒立振子に関して参考にしたのは倒立振子の制御で、下の図のような倒立振子を安定させる問題を扱っている。

20140313-01

並進運動の方程式は、台車と棒を別の剛体として立てているのだが、力のモーメントについて勉強のように重心の位置の運動方程式は力のモーメントのように作用点によって方程式が変わることはないから、台車の支点の水平位置をxとすると棒に関する縦方向の方程式は

\displaystyle m \frac{d^2}{dt^2} (x + l \sin \theta) = H \\ m \frac{d^2}{dt^2} (l \cos \theta) = V - mg

となる。一般には角運動量は中心をどこに取るかによって変わってくるので中心が動いてしまうと方程式が立てられないはず。まずは、固定された座標系の原点をもとに角運動量を求めてみる。Cは力のモーメントとしてどこに作用するかわからないのでとりあえずは0にしておく。角運動量を棒の微小な長さに対して積分して

\displaystyle xV - (x + l\sin\theta) mg = \frac{m}{2l} \frac{d}{dt} \int_0^{2l} \left\{ (x+p\sin\theta) \frac {d}{dt} (p\cos\theta) - p\cos\theta \frac {d}{dt} (x + p\sin\theta) \right\} dp \\ = \frac{m}{2l} \frac{d}{dt} \int_0^{2l} \left\{ (x+p\sin\theta) (-p\dot\theta\sin\theta) - p\cos\theta (\dot x + p\dot\theta\cos\theta) \right\} dp \\ = \frac{m}{2l} \frac{d}{dt} \int_0^{2l} \left\{ -px\dot\theta\sin\theta - p^2\dot\theta - p\dot x\cos\theta \right\} dp \\ = \frac{m}{2l} \frac{d}{dt} \big( -2l^2 x\dot\theta\sin\theta - \frac 8 3 l^3\dot\theta - 2l^2 \dot x\cos\theta \big) \\ = \frac{m}{2l} \big( -2l^2 \dot x\dot\theta\sin\theta -2l^2 x\ddot\theta\sin\theta -2l^2 x\dot\theta^2\cos\theta - \frac 8 3 l^3\ddot\theta - 2l^2 \ddot x\cos\theta + 2l^2 \dot x \dot\theta \sin\theta \big) \\ = \frac{m}{2l} \big( -2l^2 x\ddot\theta\sin\theta -2l^2 x\dot\theta^2\cos\theta - \frac 8 3 l^3\ddot\theta - 2l^2 \ddot x\cos\theta \big) \\ = m \big( -l x\ddot\theta\sin\theta -l x\dot\theta^2\cos\theta - \frac 4 3 l^2\ddot\theta - l \ddot x\cos\theta \big) \\
ここで
\displaystyle V = m \frac{d^2}{dt^2} (l \cos \theta) + mg \\ = m \frac{d}{dt} (-l\dot\theta\sin\theta) + mg \\ = m (-l\ddot\theta\sin\theta - l\dot\theta^2\cos\theta) + mg \\ = -ml\ddot\theta\sin\theta - ml\dot\theta^2\cos\theta + mg \\
としてVを代入。

\displaystyle -mlx\ddot\theta\sin\theta - mlx\dot\theta^2\cos\theta + mgx - (x + l\sin\theta) mg = m \big(-lx\ddot\theta\sin\theta - lx\dot\theta^2\cos\theta - \frac 4 3 l^2\ddot\theta - l\ddot x\cos\theta \big) \\ mgx - (x + l\sin\theta) mg = m \big( -\frac 4 3 l^2\ddot\theta - l\ddot x\cos\theta \big) \\ gx - (x + l\sin\theta) g = -\frac 4 3 l^2\ddot\theta - l\ddot x\cos\theta \\ -gl\sin\theta = -\frac 4 3 l^2\ddot\theta - l\ddot x\cos\theta \\ -g\sin\theta = -\frac 4 3 l\ddot\theta - \ddot x\cos\theta g\sin\theta = \frac 4 3 l\ddot\theta + \ddot x\cos\theta

棒の慣性モーメントJはJ=\frac 1 3 ml^2だから、これは倒立振子の制御の

\displaystyle ml\ddot x\cos\theta + (J + ml^2) \ddot\theta - mlg\sin\theta + C \dot \theta = 0

にC=0を代入したものになる。

倒立振子の制御の方程式は重心を中心とする角運動量に関するものだが、角運動量の中心が加速度運動をする場合は、非慣性系の運動方程式（慣性系 （ ガリレイ変換、遠心力、コリオリの力、フーコーの振り子 ）参照）が必要になる。一般に、質点P_1,...,P_nがx'_1,...x'_mの位置にあり、x_i + x_0 = x'_iとする。それぞれの質点P_iに力f_iがかかっている。質点P_1,...,P_nの質量の合計をM、重心の位置をx'_Mとする。x_M + x_0 = x'_M

\displaystyle \frac {d}{dt} \sum_i^n m_i x_i \times \dot x_i = \sum_i^n m_i x_i \times \ddot x_i \\ = \sum_i^n m_i x_i \times (\ddot x'_i - \ddot x_0) \\ = \sum_i^n m_i x_i \times \ddot x'_i - \sum_i^n m_i x_i \times \ddot x_0 \\ = \sum_i^n x_i \times f_i - M x_M \times \ddot x_0 \\

非慣性系では原点の加速度に応じて、「慣性力」が働くように見える。この慣性力による力のモーメントは-M x_M \times \ddot x_0で、重心に原点の加速度とは逆方向に力が働いているように見える。角運動の中心を重心に取ると、x_M=0になるので、-M x_M \times \ddot x_0 = 0となるので、たまたま、原点の加速度は角運動量には関係なくなる。

Cに関してだが、支点を中心とする角運動量を考えてみる。この場合、支点の抵抗による力のモーメントがC \dot \thetaになるのは、その他の場所を中心の角運動量を考えるより直観的にわかりやすい。支点を中心とする場合、作用点が回転の中心にあるのでVもHも角運動量に影響しなくなり、慣性力と重力だけが力のモーメントに寄与する。

\displaystyle -\ddot\theta \int_0^{2l} \frac {m}{2l} p^2 dp = -mgl \sin\theta + ml \ddot x \cos\theta + C \ddot\theta \\ -\ddot\theta \left[ \frac {m}{2l} \frac {1}{3} p^3 \right]_0^{2l} = -mgl \sin\theta + ml \ddot x \cos\theta + C \ddot\theta \\ -\ddot\theta \frac {m}{2l} \frac {8}{3} l^3 = -mgl \sin\theta + ml \ddot x \cos\theta + C \ddot\theta \\ ml \ddot x \cos\theta + m \frac {4}{3} l^2 \ddot\theta - mgl \sin\theta + C \ddot\theta = 0\\

結局、どこを角運動量の中心にとっても慣性力を適切にとれば、結果は同じになる。

カテゴリー: 物理 | コメントをどうぞ

行列の積の微分

投稿日: 2014年3月9日 作成者: kiidax

基礎知識なのかもしれないけど知らなかったのでメモ。

\displaystyle A = ( a_{ij} ) (i=1,...,l, j=1,...m)
\displaystyle B = ( b_{ij} ) (i=1,...,m, j=1,...n)
\displaystyle C = AB

とした場合、

\displaystyle \frac{dC}{dx} = \frac{dA}{dx} B + A \frac{dB}{dx}

のようにスカラー積の微分のように計算できる。

\displaystyle \frac{dC}{dx} = \frac{d}{dx} ( \sum_k^m a_{ik} b_{kj} ) (i=1,...l, j=1,...n)
\displaystyle = ( \sum_k^m (\frac{da_{ik}}{dx} b_{kj} + a_{ik} \frac{db_{kj}}{dx}) )
\displaystyle = ( \sum_k^m \frac{da_{ik}}{dx} b_{kj} ) + ( \sum_k^m a_{ik} \frac{db_{kj}}{dx} )
\displaystyle = \frac{dA}{dx} B + A \frac{dB}{dx}

ついでに回転行列の微分。

\displaystyle R(\theta) = \begin{pmatrix} \cos \theta & -\sin \theta \\ \sin \theta & \sin \theta \\ \end{pmatrix}

とすると、

\displaystyle \frac{dR(\theta)}{d\theta} = R(\theta + \frac{\pi}{2}) = R(\frac{\pi}{2}) R(\theta)

\displaystyle \frac{dR(\theta)}{d\theta} = \begin{pmatrix} -\sin \theta & -\cos \theta \\ \cos \theta & -\sin \theta \\ \end{pmatrix} = \begin{pmatrix} -\sin \theta & -\cos \theta \\ \cos \theta & -\sin \theta \\ \end{pmatrix}
\displaystyle = \begin{pmatrix} \cos (\theta + \frac{\pi}{2}) & -\sin (\theta + \frac{\pi}{2}) \\ \sin (\theta + \frac{\pi}{2}) & \cos (\theta + \frac{\pi}{2}) \\ \end{pmatrix} = R(\theta + \frac{\pi}{2}) = R(\frac{\pi}{2}) R(\theta)

極座標表現でも同様で、

\displaystyle \frac{d}{d\theta} e^{i\theta} = ie^{i\theta} = e^{\frac{\pi}{2}i}e^{i\theta} = e^{(\theta + \frac{\pi}{2})i}

カテゴリー: 物理 | コメントをどうぞ

角運動量の疑問

投稿日: 2014年3月4日 作成者: kiidax

角運動量について、説明サイトを見ていて理解に困難だった点は、角運動量について語る際に角運動量を自転として見る場合と公転として見る場合があることだ。Wikipediaの角運動量保存の法則を見ると、フィギュアスケートの高速スピンや「こま」の回転の話、後でケプラーの法則の話がでてくるが、前者は自転の話で、後者は公転の話になる。角運動量を理解するために運動量と角運動量の関係に関して検証してみた。

Wikipediaでは角運動量保存の法則の証明をしている。証明というからには、この法則が別の法則から導き出せる副次的なものということで、この別の法則は運動量保存の法則のことだ。質点に関して角運動量\mathbf L、力のモーメントを\mathbf Nとした場合、\mathbf L=m \mathbf x \times \mathbf {\dot x}を１階微分して

\displaystyle \mathbf {\dot L} = m \mathbf {\dot x} \times \mathbf {\dot x} + m \mathbf x \times \mathbf {\ddot x} = 0 + m \mathbf x \times \mathbf {\ddot x} = m \mathbf x \times \mathbf {\ddot x} = \mathbf N

から力のモーメントにより角運動量が変化することがわかる。

ところで、質点をたくさん集めてみたときに、その重心は運動の法則に従っているように見える。n個の質量m_1,..., m_nの質点\mathbf x_1,..., \mathbf x_nがあり、運動量の総和を\mathbf Pとすると、

\displaystyle \mathbf P = \sum_i m_i \mathbf x_i
\displaystyle \mathbf {\dot P} = \sum_i m_i \mathbf {\ddot x}_i

\mathbf x_iが外力\mathbf f_i、\mathbf x_jから内力\mathbf g_{ij}を受け取るとすると、

\displaystyle \mathbf {\dot P} = \sum_i m_i \mathbf {\ddot x}_i = \sum_i (\mathbf f_i + \sum_j \mathbf g_{ij}) = \sum_i \mathbf f_i + \sum_i \sum_j \mathbf g_{ij}

ここで、作用、反作用の法則から\mathbf g_{ij} + \mathbf g_{ji} = \mathbf 0だから、\displaystyle \sum_i \sum_j \mathbf g_{ij} = \mathbf 0。よって、

\displaystyle \mathbf {\dot P} = \sum_i \mathbf f_i

となり、ここで、重心を\mathbf xとすると、

\displaystyle (\sum_i m_i) \mathbf {\ddot x} = \sum_i m_i \mathbf {\ddot x}_i = \mathbf {\dot P} = \sum_i \mathbf f_i

より、重心は質点が受け取ったすべての外力の総和により、運動の法則に従っているようにふるまう。小さいものをまとめても同じ法則に従うのだ。

一方、角運動量の総和を\mathbf Lとすると、

\displaystyle \mathbf L = \sum_i m_i \mathbf x_i \times \mathbf {\dot x}_i
\displaystyle \mathbf {\dot L} = \sum_i m_i \mathbf x_i \times \mathbf {\ddot x}_i = \sum_i \mathbf x_i \times \mathbf f_i + \sum_i \sum_j \mathbf x_i \times \mathbf g_{ij}

となる。\displaystyle \sum_i \mathbf x_i \times \mathbf f_iは外力による力のモーメントの総和になるが、\displaystyle \sum_i \sum_j \mathbf x_i \times \mathbf g_{ij} = \mathbf 0とならなければ、質点を合成した場合に角運動量の保存則が成り立たない。質点を合成した場合の角運動量の保存則は、Wikipediaの作用、反作用の法則だけからは導けないらしい。運動量保存則だけでは不完全？でも同じ疑問を持った人がいるようだ。

作用、反作用の法則より、\mathbf g_{ij} + \mathbf g_{ji} = \mathbf 0だから、

\displaystyle \sum_i \sum_j \mathbf x_i \times \mathbf g_{ij} = \frac 1 2 (\sum_i \sum_j \mathbf x_i \times \mathbf g_{ij} - \sum_i \sum_j \mathbf x_j \times \mathbf g_{ij}) = \frac 1 2 \sum_i \sum_j (\mathbf x_i - \mathbf x_j) \times \mathbf g_{ij}

角運動量の保存則を成り立たせるための単純な法則を足すとすれば、\mathbf g_{ij}と\mathbf g_{ji}は\mathbf x_i - \mathbf x_jに平行とすればよい。これは重力でも成り立つので納得できる。

式を見る限りでは角運動量に、自転、公転の区別はないようだ。

外積の微分

\displaystyle \frac {\mathrm d} {\mathrm d x} (\mathbf a \times \mathbf b) = \frac {\mathrm d} {\mathrm d x} (a_2 b_3 - a_3 b_2, a_3 b_1 - a_1 b_3, a_1 b_2 - a_2 b_1) \\ = (\frac {\mathrm d} {\mathrm d x} a_2 b_3 - \frac {\mathrm d} {\mathrm d x} a_3 b_2 + a_2 \frac {\mathrm d} {\mathrm d x} b_3 - a_3 \frac {\mathrm d} {\mathrm d x} b_2, \frac {\mathrm d} {\mathrm d x} a_3 b_1 - \frac {\mathrm d} {\mathrm d x} a_1 b_3 + a_3 \frac {\mathrm d} {\mathrm d x} b_1 - a_1 \frac {\mathrm d} {\mathrm d x} b_3, \frac {\mathrm d} {\mathrm d x} a_1 b_2 - \frac {\mathrm d} {\mathrm d x} a_2 b_1 + a_1 \frac {\mathrm d} {\mathrm d x} b_2 - a_2 \frac {\mathrm d} {\mathrm d x} b_1) \\ = (\frac {\mathrm d} {\mathrm d x} a_2 b_3 - \frac {\mathrm d} {\mathrm d x} a_3 b_2, \frac {\mathrm d} {\mathrm d x} a_3 b_1 - \frac {\mathrm d} {\mathrm d x} a_1 b_3, \frac {\mathrm d} {\mathrm d x} a_1 b_2 - \frac {\mathrm d} {\mathrm d x} a_2 b_1) + (a_2 \frac {\mathrm d} {\mathrm d x} b_3 - a_3 \frac {\mathrm d} {\mathrm d x} b_2, a_3 \frac {\mathrm d} {\mathrm d x} b_1 - a_1 \frac {\mathrm d} {\mathrm d x} b_3, a_1 \frac {\mathrm d} {\mathrm d x} b_2 - a_2 \frac {\mathrm d} {\mathrm d x} b_1) \\ = \frac {\mathrm d} {\mathrm d x} \mathbf a \times \mathbf b + \mathbf a \times \frac {\mathrm d} {\mathrm d x} \mathbf b

カテゴリー: 物理 | コメントをどうぞ

力のモーメントについて勉強(2)

投稿日: 2014年3月2日 作成者: kiidax

力のモーメントについて勉強の棒による作用反作用の力と、その時の角加速度を求めてみる。本当は３次元で考えなければいけないと思うのだが、やり方がわからないので２次元で考える。まず、\mathbf x_1, \mathbf x_2を\mathbf x_0, \thetaで次のように表す。

\displaystyle \mathbf x_1 = \mathbf x_0 + L_1 (\cos \theta, \sin \theta)
\displaystyle \mathbf x_2 = \mathbf x_0 - L_2 (\cos \theta, \sin \theta)

１階微分すると

\displaystyle \dot{\mathbf x}_1 = \dot{\mathbf x}_0 + L_1 (-\sin \theta, \cos \theta) \dot \theta
\displaystyle \dot{\mathbf x}_2 = \dot{\mathbf x}_0 - L_2 (-\sin \theta, \cos \theta) \dot \theta

２階微分すると

\displaystyle \ddot{\mathbf x}_1 = \ddot{\mathbf x}_0 + L_1 (-\cos \theta, -\sin \theta) {\dot \theta}^2 + L_1 (-\sin \theta, \cos \theta) \ddot \theta
\displaystyle \ddot{\mathbf x}_2 = \ddot{\mathbf x}_0 - L_2 (-\cos \theta, -\sin \theta) {\dot \theta}^2 - L_2 (-\sin \theta, \cos \theta) \ddot \theta

となる。

\displaystyle m_1 \ddot{\mathbf x}_1 = \mathbf u - \mathbf f
\displaystyle m_2 \ddot{\mathbf x}_2 = \mathbf f

より、

\displaystyle \mathbf u - \mathbf f = m_1 \ddot{\mathbf x}_0 + m_1 L_1 (-\cos \theta, -\sin \theta) {\dot \theta}^2 + m_1 L_1 (-\sin \theta, \cos \theta) \ddot \theta
\displaystyle \mathbf f = m_2 \ddot{\mathbf x}_0 - m_2 L_2 (-\cos \theta, -\sin \theta) {\dot \theta}^2 - m_2 L_2 (-\sin \theta, \cos \theta) \ddot \theta

これを\mathbf uについて解くと、

\displaystyle \mathbf u = (m_1 + m_2) \ddot{\mathbf x}_0

これを使って、\mathbf fを表すと、

\displaystyle \mathbf f = \frac {m_2} {m_1+m_2} \mathbf u - m_2 L_2 (-\cos \theta, -\sin \theta) {\dot \theta}^2 - m_2 L_2 (-\sin \theta, \cos \theta) \ddot \theta

ここで棒は引っ張るか押すかしかできない（多分そうだと思う）ため、\mathbf fは(\cos \theta, \sin \theta)に平行になる。そのため、

\displaystyle \mathbf f \cdot (-\sin \theta, \cos \theta) = \frac {m_2} {m_1+m_2} \mathbf u \cdot (-\sin \theta, \cos \theta) - m_2 L_2 \ddot \theta = 0

\displaystyle \frac {m_2} {m_1+m_2} \mathbf u \cdot (-\sin \theta, \cos \theta) = m_2 L_2 \ddot \theta

両辺にLをかけて、

\displaystyle \frac {m_2} {m_1+m_2} L \mathbf u \cdot (-\sin \theta, \cos \theta) = m_2 L_2 L \ddot \theta
\displaystyle L_1 \mathbf u \cdot (-\sin \theta, \cos \theta) = \frac {m_1 m_2} {m_1+m_2} L^2 \ddot \theta

左辺は力のモーメントだから、Nとおき、また、I = \frac {m_1 m_2} {m_1+m_2} L^2より、

\displaystyle N = I \ddot \theta

よって、力のモーメントが角加速度に寄与していることがわかる。

カテゴリー: 物理 | コメントをどうぞ

力のモーメントについて勉強

投稿日: 2014年3月1日 作成者: kiidax

ネットの情報を漁っていると一輪車を倒れないようにすることは倒立振り子の安定化問題に似ていることが分かってきた。倒立振り子の安定化問題は制御工学の定番の課題らしく情報量は多い。運動方程式からモデルを作り安定するフィードバックを決定するのだが、高校物理の知識のまま安定化してしまった脳みそには回転を考慮した運動方程式が理解できないので剛体の力学を参考に勉強してみた。

回転を考えない場合、剛体にかかる力を求めて運動方程式に当てはめればいい。この場合、剛体のどの部分に力をかけるのかは考えない。一方、回転の場合は剛体にかかるトルクを求めて運動方程式に当てはめる。トルクは重心から力の作用点へのベクトルと力のベクトルの外積で定義される。外積なので重心方向に直角な分を取り出したものと考えることができる。ボールをボールの中心方向に押せば、回転せずに移動するが、表面と平行に押せば回転することから直観的に理解できる。

そして、剛体の運動は、高校物理でおなじみの並進運動の運動方程式と回転運動の運動方程式で２つの方程式で表すことができる。

ここでつまづいたのは、重心方向を向いている力は回転運動に寄与しないのであれば、重心方向を向いていない力は並進運動に寄与しないのではないかという勘違い。高校では力の合成と分解を勉強するが力の分解に似たことが並進運動と回転運動にも起こるのではないかと思ったが、そうでないことを実際検証してみた。下のように２つの物体が質量が無視できる変形しない棒でつながっていて、２つの物体に作用、反作用に法則により力\mathbf Fがかかっているとする。

20140301-01

２つの物体の重心を \mathbf x_0 とすると

\displaystyle m_0 = m_1 + m_2

\displaystyle m_0 \mathbf x_0 = m_1 \mathbf x_1 + m_2 \mathbf x_2

さらに、外部から片方の物体に力\mathbf uをかける。\mathbf {x_0,x_1,x_2,f,u}は時間に依存するとして、ニュートン記法を使うと、

\displaystyle m_1 \ddot{\mathbf x}_1 = \mathbf u - \mathbf f

\displaystyle m_2 \ddot{\mathbf x}_2 = \mathbf f

よって、

\displaystyle m_0 \ddot {\mathbf x}_0 = m_1 \ddot{\mathbf x}_1 + m_2 \ddot{\mathbf x}_2 = \mathbf u - \mathbf f + \mathbf f = \mathbf u

となり\mathbf fは相殺されて、２つの物体を合成した剛体はあたかも、質量m_0の物体に力\mathbf uをかけたかのように並進運動する。これから\mathbf uは重心を向いているかどうかに関係なく、剛体の並進運動に寄与することがわかる。

ところで、２つの物体の運動エネルギーの和は、

\displaystyle K_1 + K_2 = \frac 1 2 m_1 \dot {\mathbf x}_1^2 + \frac 1 2 m_2 \dot {\mathbf x}_2^2

２つの物体を合成した剛体の運動エネルギーは

\displaystyle K_0 = \frac 1 2 m_0 \dot {\mathbf x}_0^2

その差は

\displaystyle K' = K_1 + K_2 - K_0 = \frac 1 2 m_1 \dot {\mathbf x}_1^2 + \frac 1 2 m_2 \dot {\mathbf x}_2^2 - \frac 1 2 m_0 \dot {\mathbf x}_0^2

\displaystyle = \frac 1 2 ( m_1 \dot {\mathbf x}_1^2 + m_2 \dot {\mathbf x}_2^2 - \frac 1 {m_1 + m_2} (m_1 \dot {\mathbf x}_1 + m_2 \dot {\mathbf x}_2)^2 )

\displaystyle = \frac 1 {2(m_1+m_2)} ( m_1^2 \dot {\mathbf x}_1^2 + m_1 m_2 \dot {\mathbf x}_1^2 + m_1 m_2 \dot {\mathbf x}_2^2 + m_2^2 \dot {\mathbf x}_2^2 - m_1^2 \dot {\mathbf x}_1^2 - 2 m_1 m_2 \dot {\mathbf x}_1 \cdot \dot {\mathbf x}_2 - m_2^2 \dot {\mathbf x}_2^2 )

\displaystyle = \frac 1 {2(m_1+m_2)} ( m_1 m_2 \dot {\mathbf x}_1^2 + m_1 m_2 \dot {\mathbf x}_2^2 - 2 m_1 m_2 \dot {\mathbf x}_1 \cdot \dot {\mathbf x}_2 )

\displaystyle = \frac {m_1 m_2} {2(m_1+m_2)} ( \dot {\mathbf x}_1^2 + \dot {\mathbf x}_2^2 - 2 \dot {\mathbf x}_1 \cdot \dot {\mathbf x}_2 )

\displaystyle = \frac {m_1 m_2} {2(m_1+m_2)} ( \dot {\mathbf x}_1 - \dot {\mathbf x}_2 )^2

なので常に正になる。つまり、剛体を動かすのに使われたエネルギーは見かけの運動エネルギーとは別に何らかのエネルギーとして使われていることがわかる。

剛体が回転運動をしていないときには\dot {\mathbf x}_1^2 = \dot {\mathbf x}_2^2になるので、K' = 0になる。

棒の長さがLで剛体が角速度\omega回転している場合、それぞれの物体の重心からの距離は、

\displaystyle L_1 = \frac {m_2} {m_1+m_2} L

\displaystyle L_2 = \frac {m_1} {m_1+m_2} L

として、

\displaystyle | \dot {\mathbf x}_1 - \dot {\mathbf x}_0 | = L_1 \omega

\displaystyle | \dot {\mathbf x}_2 - \dot {\mathbf x}_0 | = L_2 \omega

でお互いに反対方向を向いてるので、

\displaystyle | \dot {\mathbf x}_1 - \dot {\mathbf x}_2 | = (L_1 + L_2) \omega = L \omega

よって、

\displaystyle K' = \frac {m_1 m_2} {2(m_1+m_2)} L^2 \omega^2

となる、一方、剛体の重心を中心とした場合の慣性モーメントは

\displaystyle I = m_1 L_1^2 + m_2 L_2^2 = \frac{m_1 m_2^2}{(m_1+m_2)^2}L^2 + \frac{m_1^2 m_2}{(m_1+m_2)^2}L^2 = \frac{m_1 m_2}{m_1+m_2}L^2

となるので、

K' = \frac 1 2 I \omega^2 となり、外部からの力で与えられたエネルギーとの差は剛体の見かけの運動エネルギーの差は回転エネルギーの定義と一致する。

試しにJavaScriptでシミュレーションをしてみる。剛体にかかる力を計算するのが面倒なので、赤と青の物体はばねでつながっている。緑の物体は赤と青を合成した剛体になっている。赤の物体と緑の物体に同じ力をかけて、その軌道を見ると、赤と青の物体の重心は緑の物体に見事一致している。

20140301-02

moment.html

<!DOCTYPE html>
<html lang="en">
    <head>
        <meta charset="utf-8" />
        <meta http-equiv="X-UA-Compatible" content="ie=edge" />
        <title>Moment Test</title>
        <script src="moment.js"></script>
    </head>
    <body>
        <h1>Moment Test</h1>
        <svg id="space" width="800" height="600" style="border: 5px solid blue; overflow: hidden;">
            <circle id="x1" cx="10" cy="10" r="2" fill="red"/>
            <circle id="x2" cx="10" cy="10" r="4" fill="blue"/>
            <circle id="x3" cx="10" cy="10" r="3" fill="green"/>
        </svg>
    </body>
</html>

moment.js

(function () {
    "use strict";

    window.addEventListener("load", function () {
        var frameIndex = 0;
        var startTime = new Date().getTime();
        var deltaTime = 1 / 24;
        var endFrameIndex = 10 / deltaTime;
        var x1Elem = document.getElementById("x1");
        var x2Elem = document.getElementById("x2");
        var x3Elem = document.getElementById("x3");
        var m1 = 1.0;
        var m2 = 2.0;
        var m3 = m1 + m2;
        var x1 = [1, 50];
        var x2 = [1, 58];
        var x3 = [(m1 * x1[0] + m2 * x2[0]) / m3, (m1 * x1[1] + m2 * x2[1]) / m3];
        var v1 = [0, -20];
        var v2 = [10, -10];
        var v3 = [(m1 * v1[0] + m2 * v2[0]) / m3, (m1 * v1[1] + m2 * v2[1]) / m3];
        var G = 9.8;
        var restDist = 6;
        var springValue = 30.0;

        function next() {
            var startTime = new Date().getTime();
            x1[0] += v1[0] * deltaTime;
            x1[1] += v1[1] * deltaTime;
            x2[0] += v2[0] * deltaTime;
            x2[1] += v2[1] * deltaTime;
            x3[0] += v3[0] * deltaTime;
            x3[1] += v3[1] * deltaTime;
            var dx = x1[0] - x2[0];
            var dy = x1[1] - x2[1];
            var dist = Math.sqrt(dx * dx + dy * dy);
            if (dist > 0) {
                var force = (dist - restDist) * springValue;
                v1[0] -= dx * force * deltaTime / (m1 * dist);
                v1[1] -= dy * force * deltaTime / (m1 * dist);
                v2[0] += dx * force * deltaTime / (m2 * dist);
                v2[1] += dy * force * deltaTime / (m2 * dist);
            }
            v2[1] += G * deltaTime / m2;
            v3[1] += G * deltaTime / m3;
        }

        function render() {
            x1Elem.setAttribute("cx", x1[0] * 10);
            x1Elem.setAttribute("cy", x1[1] * 10);
            x2Elem.setAttribute("cx", x2[0] * 10);
            x2Elem.setAttribute("cy", x2[1] * 10);
            x3Elem.setAttribute("cx", x3[0] * 10);
            x3Elem.setAttribute("cy", x3[1] * 10);
        }

        var timer = window.setInterval(function () {
            try {
                var nextFrameIndex = (new Date().getTime() - startTime) / (1000 * deltaTime);
                while (frameIndex <= nextFrameIndex) {
                    next();
                    frameIndex++;
                }
                if (frameIndex > endFrameIndex) {
                    window.clearInterval(timer);
                    timer = null;
                }
                render();
            } catch (e) {
                window.clearInterval(timer);
                timer = null;
            }
        }, deltaTime * 1000 * 0.9);
    });
})();
カテゴリー: 物理 | 4件のコメント

動作とアニメーションのスクリプトを分離

投稿日: 2014年2月11日 作成者: kiidax

同じ動きをするキャラクターで動作のためのスクリプトを共有したいが見た目が違う場合を考えて、
２輪車用のスクリプトはアニメーションのスクリプトを分離できるようにして、スケーターにも使えるようにしている。

20140211-02

具体的には、動作のためのスクリプトの

            self.leanAngle = math.asin(-normalVect.dot(worldLeftVect))

の部分でGameObjectに現在の車体の傾きを設定しておく。アニメーションのためのスクリプトではこの値を見てアニメーションを実行する。

20140211-01

動作のスクリプトはGirl Empty Objectに設定してあって、アニメーションのスクリプトはRigGirlに設定してある。GirlSkate_Leanは傾きごとの姿勢のためのアニメーションでフレーム1では右に傾き、61で正立、121で左に傾いている。アニメーションのactuatorのframeプロパティを設定して、傾きに従って姿勢が変わるようになっている。

20140211-03

import math

def update(cont):
    owner = cont.owner
    parent = owner.parent
    acc = cont.actuators["GirlSkate_Lean"]
    cont.activate(acc)
    frame = parent.leanAngle / math.pi * 180 + 61
    acc.frame = min(max(1, frame), 121)

20140211-05

20140211-04

カテゴリー: Blender | コメントをどうぞ