【Unity道場スペシャル 2018仙台】物理シミュレーション完全マスター

෺ཧγϛϡϨʔγϣϯ
‫׬‬શϚελʔ 仙台版 フィールド・エンジニア 安原 祐二

積分 完全マスター

３個 ４列 ３×４＝12 ３個×４列＝12個

３m ４m ３×４＝12 2 ３m×４m＝12m

長さと長さを掛けると ３m 面積？ ４m ３×４＝12 2 ３m×４m＝12m

３m ３個 ４m ４列 12個 ぜんぜん違う話！ 12m2

長さと長さを掛けると ３m 面積？ ４m 積分の神秘は ３×４＝12 2 ３m×４m＝12m これと同じ

３ 極小の長さぶんを

３ 範囲ぶん足す

３m ４m 12m2 単位が変わる

y x x x 極小長さ： 方向の微小区間 dx（デルタ ）という

y x 面積を出す x で積分する

長さmを長さmで積分→面積m2

長さmを長さmで積分→面積m2 2 3 面積m を長さmで積分→体積m 積分すると必ず単位が変化する

Z 積分の計算には を使う （いんてぐらる） 計算は呪文のようなもの 理解そのものとは無関係

例：球の体積（部分） 呪文 V y = = S = r t = x = 理解 Z Z t Sdx r t ⇡(r2 x2 )dx r 3 x ⇡[r2 x ]t r 3 3 t ⇡{r2 t ( r2 r 3 t3 2 3 2 ( + r t + r )⇡ 3 3 円の面積を重ねる r3 )} 3

〜今回のおはなし〜 最小実装 外力 接触 力積 慣性テンソル

物理シミュレーションの 最小実装

2 加速度m/s を時間tで積分→速度m/s

2 加速度m/s を時間tで積分→速度m/s velocity += acceleration * dt;

2 加速度m/s を時間tで積分→速度m/s velocity += acceleration * dt; 速度m/sを時間tで積分→位置m

2 加速度m/s を時間tで積分→速度m/s velocity += acceleration * dt; 速度m/sを時間tで積分→位置m position += velocity * dt; これが物理シミュ

動画 最小実装はこれだ！

public class MyRigidbody : MonoBehaviour {
public Vector3 acceleration;
public Vector3 velocity;
public Vector3 position;
const float dt = 1f/60f;
public void AddForce(Vector3 force) {
acceleration += force;
}
void FixedUpdate() {
velocity += acceleration * dt;
position += velocity * dt;
if (position.y < 0.5f) {
velocity = -velocity;
}
transform.position = position;
acceleration = Vector3.zero;
}
}

Unityの物理更新と FixedUpdate

60fpsのゲームは dt １秒 60回 １秒に60回画面を更新している

60fpsのゲームは dt １秒 60回 １秒に60回画面を更新している

1/60秒 処理 余裕 余裕期間は待機

余裕 1/60秒

余裕 1/60秒

余裕 1/60秒

余裕 1/60秒

超過 1/60秒

待機 1/60秒

責任重大 1/60秒

1/60秒 通常更新 物理計算 FixedUpdate Update

1/60秒 通常更新 物理計算 FixedUpdate Update 物理更新

物理計算 FixedUpdate 物理計算 FixedUpdate 物理更新 責任重大 フレーム 物理更新 Update 物理更新は複数回実行される

責任重大 フレーム 物理計算 FixedUpdate 物理計算 FixedUpdate Time.fixedDeltaTime Update Time.deltaTime ※TimeManagerの設定で事情は変化 固定 Updateは 変動 dtが変動

物理計算 FixedUpdate 物理計算 FixedUpdate OnCollisionEnter OnCollisionExit OnTriggerEnter 物理更新 OnCollision*は から

外力

運動方程式 F = ma F :力 m:質量 a :加速度

運動方程式 F = ma F :力 m:質量 a :加速度 物理シミュレーションは 位置を知りたい F a= m 質量は固定なので 加速度は力

加速度があれば位置を出せる 2 加速度m/s を時間tで積分→速度m/s velocity += acceleration * dt; 速度m/sを時間tで積分→位置m position += velocity * dt; これが物理シミュ

力ベクトルは足せる

物理計算 FixedUpdate 物理更新

物理更新 物理計算 FixedUpdate

重力 物理更新 物理計算 様々な力を 反発力 摩擦力 FixedUpdate +) AddForce トータル外力 足す

AddForceは足すだけ void FixedUpdate() { Vector3 f = new Vector3(0f, -9.81f, 0f); // 重力 f += new Vector3(1f, 0f, 0f); // 風の影響 f += new Vector3(0f, 1f, 0f); // 浮力 rigidbody.AddForce(f); } こんなふうにあらかじめ足しても良い AddForceをそれぞれ呼ぶのと同じ

外力の応用 浮力

浮力とは 物体が押しのけた水の重さ 潜っている部分の体積 潜っている部分の重心 かなりたいへん

浮力の計算はたいへん 球で代用する 単純な式で出せる！ 向きを考慮しなくていい 波の対応も容易

動画 浮力っぽい

浮力の計算はたいへん y 球で代用する S r t x t3 2 V =( + r2 t + r3 )⇡ 3 3 単純な式で出せる！ 向きを考慮しなくていい 波の対応も容易

力積

運動方程式 F = ma F :力 m:質量 a :加速度 加速度は 微小期間の速度変化 dv a= dt

加速度は 微小期間の速度変化 F = ma dv a= dt dv F = m dt これも運動方程式

加速度は 微小期間の速度変化 力積 Fdt= mdv 運動量（速度）の変化 dv F = m dt

力積も運動量も F = ma の変形に過ぎない 加速度は 微小期間の速度変化 力積 Fdt= mdv 運動量（速度）の変化 dv F = m dt

「力積は運動量を変化させる」 Fdt= mdv dt期間に発生した力の合算を 速度変化に適用する術

「力積は運動量を変化させる」 Fdt= mdv dt期間に発生した力の合算を 速度変化に適用する術 瞬間的に発生する力に使う

爆発で吹っ飛ばしたい 爆発力が存在する期間 AddForce を呼べばいい が、爆発は一瞬 dtが大き過ぎて一瞬を作れない

ForceMode.Impulse rigidbody.AddForce(force, ForceMode.Impulse);

ForceMode.Impulse rigidbody.AddForce(force, ForceMode.Impulse); 内部的には速度の直接変更 （Fdt= mdv の右辺） 質量は考慮するが は考慮しない dt

ForceMode.Impulse rigidbody.AddForce(force, ForceMode.Impulse); 内部的には速度の直接変更 （Fdt= mdv の右辺） 質量は考慮するが は考慮しない dt dt に依存しない速度変化になる OnCollisionEnterなどの瞬間的な力

誤 dtを変えると動作が変わってしまう if (hit) { rb.AddForce(force); } 正 瞬間的な外力には力積を使用する if (hit) { rb.AddForce(force, ForceMode.Impulse); }

Drag

Dragとは 速度に比例した抵抗力 例：空気抵抗

v 指定終端速度 を得るための外力F v k m F = 1 kdt m:質量 k :Drag 参 考 [Unity]RigidbodyのDragから終端速度を得る https://qiita.com/yuji_yasuhara/items/1f438f0f27f5ef854a73

• https://qiita.com/yuji

考 参 〜実践的な物理入門〜 [Unity道場 札幌スペシャル] プロが教える脱初心者スクリプト術！ https://www.youtube.com/watch?v=FqjM9oujyNE

• https://www.youtube.com/watch?v

物理シミュレーションの ウソ

仮に物理シミュが完璧なら こんなに設定は必要ない

なにがウソなのか？ dt が大きすぎる あらゆる問題の根源

なにがウソなのか？ dt が大きすぎる あらゆる問題の根源 あらゆる工夫 物理シミュレーションの品質

実現困難な例 Wikipedia「Newton’s cradle」より https://en.wikipedia.org/wiki/Newton%27s_cradle

• https://en.wikipedia.org/wiki/Newton

接触

接触していない 移動

接触！（図形的に算出） 接触＝必ず潜り込みが発生

潜り込みを直す 位置と角度を調整 Solve(ソルブ) と呼ばれる処理

接触は 物理シミュレーションの 最大の難関 多くの進化や工夫は ここに集約される

下にも物体があったら 調整したせいで潜り込む！

何度も調整 負荷大 求める完全さは アプリの特徴次第

慣性テンソル

「質量m」 dvv 外力 F = m dt v：速度 質量は「動かしにくさ」

移動における「質量m」 dvv 外力 F = m dt v：速度 回転ににおける「慣性テンソル I 」 ! d! トルク T = I dt !：角速度 慣性テンソルは「回しにくさ」

動画 慣性テンソルを可視化

回しにくい 回しにくい 回しやすい

ここからは武装が必要 外積 テンソル 行列 トルク

外積 〜ベクトルの立体関係〜 A 外積 A×B B どちらのベクトルにも 垂直なベクトル

行列 〜ベクトルに掛けるもの〜 A = PB BにPを掛けたらA AもBもひとつの値（スカラー）なら大きさしか変えない が、ベクトルの場合は大きさも方向も変えるのでPは行列 0 1 Ax @ Ay A Az = 0 P00 @ P10 P20 P01 P11 P21 10 1 P02 Bx P12 A@ By A P22 Bz

テンソル 〜単なる呼びかた〜 スカラー ０階テンソル ベクトル １階テンソル 行列 ２階テンソル 今回扱うのは２階テンソル

トルク〜回転力の記述〜 回転には必ず回転軸がある 回転させる力をトルクと呼ぶ トルク（と角速度）は 軸ベクトルで表現

再掲 加速度があれば位置を出せる 2 加速度m/s を時間tで積分→速度m/s velocity += acceleration * dt; 速度m/sを時間tで積分→位置m position += velocity * dt; これが物理シミュ

加速度があれば位置を出せる 加速度があれば（速度を出せるから）位置を出せる

加速度があれば位置を出せる 加速度があれば（速度を出せるから）位置を出せる トルクがあれば姿勢を出せる

加速度があれば位置を出せる 加速度があれば（速度を出せるから）位置を出せる トルクがあれば姿勢を出せる トルクから角速度を出せば姿勢を出せる

加速度があれば位置を出せる 加速度があれば（速度を出せるから）位置を出せる トルクがあれば姿勢を出せる トルクから角速度を出せば姿勢を出せる トルクと角速度の関係は単純ではない

動画 衝撃的な事実

動画 期待した動き 実際の動き

角速度ω 点p r 重心 v ある点pにおける速度vは 重心からのベクトルrと 角速度ωと垂直なので v =!⇥r と外積で表せる

角速度ω 点p r 重心 v v =!⇥r f 速度vで移動している点pは 働いた力fの結果そうなっている dv m =f dt

トルクT f f f f 点p r 重心 f トルクTは すべての点で働いた力fを 距離を考慮して合計したもの 外積 r×f をぜんぶ足す T = X X r⇥f はぜんぶ足すの意味

角速度ω トルクT r 重心 f 質点が対称でない場合 トルクT（r×f の総計）と 角速度ω（運動の結果）は 軸が異なる

動画 慣性テンソルの性質

運動方程式 dv f =m dt v =!⇥r 角速度と速度の関係 d f = m (! ⇥ r) dt 力と角速度の関係！

力と角速度の関係 T = X T = r⇥f X d f = m (! ⇥ r) dt 力とトルクの関係 d r ⇥ m (! ⇥ r) dt トルクと角速度の関係！

T = X d r ⇥ m (! ⇥ r) dt d! 右辺はベクトル dt を線形変換している この変換は必ず行列 I とおけて（計算略） d! T =I dt トルクと角速度の関係

慣性テンソルとは d! T =I dt この I すなわち トルクベクトルと角速度ベクトルを関連づける 行列（２階テンソル）のこと 複雑な状況を行列ひとつで吸収できる驚異

参考 ~r ⇥ (~! ⇥ ~r) = (~r · ~r)~! = |~r|2 ! ~ (~r · ! ~ )~r (rx !x + ry !y + rz !z )~r0 1 ~ ⇥ (B ~ ⇥ C) ~ A ~ · C) ~ B ~ (A ~ · B) ~ C ~ = (A rx = |~r|2 ! ~ (rx !x + ry !y + rz !z ) @ ry A rz 0 2 1 密度と rx ! x + rx ry ! y + rz rx ! z = |~r|2 ! ~ @ rx ry !x + ry2 !y + ry rz !z A 物体形状で を 2 rz r x ! x + r y rz ! y + r z ! z 0 2 10 1 積分したものが rx rx ry rz rx !x ry2 ry r z A @ ! y A = (rx2 + ry2 + rz2 )~ ! @ rx r y !z rz r x r x r y rz2 0 =@ ry2 rz2 + rx r y rz rx rx r y rx2 + rz2 ry rz 1 rz rx ry rz A ! ~ rx2 + ry2 慣性テンソル I

例：立方体の慣性テンソル 0 @ ry2 rz2 + rx ry r z rx m abc Z rx ry rx2 + rz2 ry rz c 2 c 2 Z b 2 b 2 m 2 = (b + c2 ) 12 0 I=@ Z y 1 r z rx ry rz A rx2 + ry2 a 2 c b (y 2 + z 2 )dxdydz a 2 m 2 12 (b + c2 ) 0 0 x a z ９要素すべてを計算 1 0 0 m 2 2 A (c + a ) 0 12 m 2 2 0 (a + b ) 12

さて Rigidbody.inertiaTensor Rigidbody.inertiaTensorRotation これはなんだ

型が… Vector3 Rigidbody.inertiaTensor Rigidbody.inertiaTensorRotation Quaternion どういうこっちゃ

行列は二種類ある 対角化できるか、できないかだ！ 慣性テンソルは対角化できる行列

0 = a00 @ a10 a20 0 b00 @ b10 b20 a01 a11 a21 b01 b11 b21 1 対角化とは a02 a12 A a22 10 b02 A00 b12 A @ 0 b22 0 0 A11 0 10 0 b00 0 A @ b10 A22 b20 対角行列 対角行列に変換する術 b01 b11 b21 1 b02 b12 A b22 1

0 = I Ir a00 @ a10 a20 0 b00 @ b10 b20 a01 a11 a21 b01 b11 b21 1 a02 a12 A a22 10 b02 A00 b12 A @ 0 b22 0 0 0 b00 0 A @ b10 A22 b20 Id A11 0 10 b01 b11 b21 Ir 1 1 b02 b12 A b22 1

0 = I Ir a00 @ a10 a20 0 b00 @ b10 b20 a01 a11 a21 b01 b11 b21 1 a02 a12 A a22 10 b02 A00 b12 A @ 0 b22 0 0 0 b00 0 A @ b10 A22 b20 Id A11 0 10 b01 b11 b21 Ir I = IrIdIr 1 1 b02 b12 A b22 1 1

I = IrIdIr 1 Id Rigidbody.inertiaTensor 対角行列なのでVector3で表現可能 Ir Rigidbody.inertiaTensorRotation 回転行列なのでQuaternionで表現可能

動画 具体的に慣性テンソルを使う

public class SimpleRotator : MonoBehaviour {
void Update() {
if (Input.GetMouseButton(0)) {
var rb = GetComponent<Rigidbody>();
var omega = new Vector3(0f, 1f, 0f);
var R = transform.rotation;
var RI = Quaternion.Inverse(transform.rotation);
var Id = rb.inertiaTensor;
var Ir = rb.inertiaTensorRotation;
var IrI = Quaternion.Inverse(Ir);
var torque = R * Ir * Vector3.Scale(Id, IrI * RI * omega);
rb.AddTorque(torque, ForceMode.Impulse);
}
}
}

鉛直の角速度をもたらすトルクを得るサンプル
（現在の姿勢も考慮したもの）

おしまい