[Unity]RigidbodyのDragから終端速度を得る

Unity の Rigidbody に Dragと一定外力が与えられたときの、終端速度(terminal velocity)を出してみよう。

#実装

通常の運動方程式をプログラムで書くと、次のようなコードになる。

void update() {
    float a = F/m;   // 加速度は外力F割ることの質量m
    v = v + a*dt;    // 速度に加速度×Δtを足す
    x = x + v*dt;    // 位置に速度×Δtを足す
}

簡単。シンプル。美しい。今回は理論の話になるので += 演算子はあえて使わないでおく。
で、Unity では抵抗力として Drag というのがある。この記事は数式が多いので、わかりやすいようにカタカナでドラッグと呼ぼう。
このドラッグの値を k とおく。以下頻繁に出てくる k はドラッグだよ。覚えといてくださいね。

んでドラッグの影響について、Unityでは次のようなコードで実装されている。けっこう自信ある。

void update() {
    float a = F/m;
    v = v + a*dt;
    v = v * (1f - k*dt);  // これ。k は Drag
    x = x + v*dt;
}

増えた行がドラッグの適用だ。このようにドラッグは抵抗力という定義にもかかわらず外力として扱われておらず、速度を直接補正している。
なぜこんなことをしているかというと、ドラッグは逆方向に力を発生させるので、ドラッグが大きな値でなおかつΔtも大きい場合に発散してしまうケースがあるからだ。とくに強いバネを扱うと、簡単に発散する。抵抗の適用を遅らせることで、同フレームの速度を参照することができるので、発散を抑制することができる。ぼくもかつて自分のプロジェクトで発散に悩まされて、この方法に辿り着いた。
抵抗力が力として実装されていないおかげで解析はけっこう面倒なことになるのだけど、頑張って解いてみよう。

#解析

##抵抗式の外力換算

プログラムの
v = v * (1f - kdt);
に注目すると、ドラッグ k は毎フレーム -vk の加速度を発生させているのと同義だ。という理解で、以下ドラッグによる影響を外力に換算する。

v' = v(1-k\Delta t) 　・・・・・・(1)\\
(1)より、 \\
v' = v - vk\Delta t \\
v を移項して　\\
v'-v = -vk\Delta t \\
v' は計算後、v は計算前なので v'-v は \Delta v の意味になるから \\
\Delta v = -vk\Delta t \\
\Delta t で両辺を割って　\\
\frac{\Delta v}{\Delta t} = -vk \\
左辺は加速度の意味になるので、運動方程式 F = ma よりドラッグの影響力 F_{k} は\\
F_{k} = m\frac{\Delta v}{\Delta t} = -mvk ・・・・・・(2)\\

やったね。これでドラッグkが及ぼす外力は、 -mvk （質量と速度とドラッグを乗算して符号を反転させたもの）であることがわかった。

##運動方程式の解析

次はドラッグと一定の外力が加わったときの運動方程式を立て、終端速度を導く。

ma = F \\
いわゆる運動方程式。このFに一定の外力 F_{0} とドラッグの F_{k} を入れて\\
ma = F_{0} + F_{k} \\
F_{k}を先ほどの結果(2)で置き換えて、\\
ma = F_{0} - mkv \\
両辺をmで割って\\
a = \frac{F_{0}}{m} - kv \\
加速度を速度の微分で表すと \\
\frac{\Delta v}{\Delta t} = \frac{F_{0}}{m} - kv \\
右辺をvの係数でくくると \\
\frac{\Delta v}{\Delta t} = -k(v - \frac{F_{0}}{mk}) \\
vに関する項で両辺を割ると\\
\frac{1}{v-\frac{F_{0}}{mk}} \frac{\Delta v}{\Delta t} = -k \\
\Delta tを右辺に移して \\
\frac{\Delta v}{v-\frac{F_{0}}{mk}} = -k \Delta t \\
両辺を積分して \\
\int{\frac{\Delta v}{v-\frac{F_{0}}{mk}}} = \int{-k \Delta t} \\
積分を計算して \\
\ln{|v-\frac{F_{0}}{mk}|} = -kt + C \\
logを展開して \\
v-\frac{F_{0}}{mk} = \pm{e^{-kt+C}} \\
v について解くと \\
v = \frac{F_{0}}{mk} \pm{e^{-kt+C}} \\
\lim_{t \to \infty}{\pm{e^{-kt+C}}} = 0 より \\
\lim_{t \to \infty}{v} = \frac{F_{0}}{mk} ・・・・・・(3)\\

というわけで一定の外力が加わったときの終端速度を求めることができた。

##補正

が！まだ終わらない。最初に述べたドラッグの変則的な適用は、１フレーム未来の速度を計算していることになるため、１フレーム分の補正が必要となる。
ドラッグの式を再び適用する。

(1)を再び\\
v' = v(1-k\Delta t) \\
vに先ほど(3)で算出した終端速度を代入して \\
v' = \frac{F_{0}}{mk}(1-k\Delta t) \\

このv'が、正真正銘の終端速度になる。ばんざい。

#検証

逆に、目的の速度に収束させようと思った場合の外力は、上式を変形して

F_{0} = \frac{mkv'}{1-k\Delta t} \\

こうなる。この式を実際のプログラムで試してみよう。

void FixedUpdate()
{
　　　　　　　　var rb = GetComponent<Rigidbody>();
　　　　　　　　var target_velocity = new Vector3(10f, 0f, 0f); // 10m/s の速度が目標
　　　　　　　　rb.AddForce(target_velocity*rb.mass*rb.drag/(1f-rb.drag*Time.fixedDeltaTime));
　　　　　　　　Debug.Log(rb.velocity.magnitude);
}

これをRigidbodyを持つ適当なオブジェクトに適用する。再生しながら drag や mass を変更しても、速度は変わらないことが確認できる。AddForceに渡す値が正しく調整されているからだ。そしてその速度は狙った通りの値になる。理屈と現実が一致するこの嬉しさよ。めでたしめでたし。