ドローン航空力学：物理基礎と Python シミュレーション

親しみやすく、卓越した航空力学の性能を備えたドローンの開発・理解に不可欠な物理基礎を解説します。本稿は、マルチローター向けの Python シミュレーションフレームワークに関する研究および 2023 年の AIAA DASC での発表内容に基づいて構成されています。また、表記法はチャールズ・タイトラー氏の研究成果を参考にしています。

読者は線形代数、微積分、古典力学の基礎知識をご確認ください。X、Reddit、HackerNews で議論に参加しましょう。

---

1. 飛行器の記述

マルチローター UAV（無人航空機）は、6 つの自由度を持つものとしてモデル化されます：

• 線形運動: 3 つの直交軸 ($x, y, z$)
• 回転運動: 3 つの角速度軸（ロール $\phi$、ピッチ $\theta$、ヨー $\psi$）

1.1 座標系 (NED)

共通合意された North-East-Down (NED) 系が採用されます。正方向は以下の通りです。

• $x$: 前方（北の方位）
• $y$: 右方（東の方位）
• $z$: 下方

1.2 参照系 (フレーム)

状態を表現するために、以下の 2 つの参照系を使用します。

• 基準的・慣性系 ($n$)
  • 静的な参照系で、恣意的な全般的な全球方向と整列しています。
  • ベクトル表記：$\hat{n} = [\hat{x}^n, \hat{y}^n, \hat{z}^n]^T$
• 体固定・非慣性系 ($b$)
  • 剛体の重心を中心として整列された参照系です。
  • 飛行器とともに移動し回転します（原点をドローンの重心に接続）。
  • ベクトル表記：$\hat{b} = [\hat{x}^b, \hat{y}^b, \hat{z}^b]^T$

1.3 線形表現と角表現

位置と速度

• 飛行器の位置: 慣性系に対する体固定系の離位 $\hat{r}^n$ として表現されます。
• 飛行器の速度: 体固定系に対する慣性系への相対速度を、体固定系の座標で表現します ($\hat{v}^b$)。

姿勢角と角速度 (Tait-Bryan angles)

慣性参照系における姿勢は、以下の順序で記述されます：

1. ヨー ($\psi$): 慣性 $z$ 軸周りの回転 $R(\psi)$
2. ピッチ ($\theta$): 新しい $y$ 軸周りの回転 $R(\theta)$
3. ロール ($\phi$): さらに新しい $x$ 軸周りの回転 $R(\phi)$

最終的な体参照系は $\hat{n}_{\psi,\theta,\phi}$ で表されます。

角速度の行列方程式: 慣性系での各軸周りの角速度（$\dot{\phi}, \dot{\theta}, \dot{\psi}$）と、体固定系で観測される角速度 ($\omega_x, \omega_y, \omega_z$) の間には以下の関係があります。

$$ \begin{bmatrix} \omega_x \ \omega_y \ \omega_z \end{bmatrix} = R(\phi)\cdot R(\theta) \begin{bmatrix} 0 \ 0 \ \dot{\psi} \end{bmatrix} + R(\phi) \begin{bmatrix} 0 \ \dot{\theta} \ 0 \end{bmatrix} + \begin{bmatrix} \dot{\phi} \ 0 \ 0 \end{bmatrix}

\begin{bmatrix} \dot{\phi} + \dot{\psi} \sin{\theta}\

• \dot{\psi} \sin{\phi} \cos{\theta} + \dot{\theta}\ \dot{\psi} \cos{\phi} \cos{\theta} \end{bmatrix} $$

回転行列 (Rotation Matrix)

体固定系から慣性系へのベクトル変換は、以下の回転行列 $R_b^n$ を使用します（ただし記号の定義により注意が必要です。ここでは各軸周りの基本回転行列を定義します）。

余弦 ($c$) と正弦 ($s$) はそれぞれ $\cos, \sin$ を表します。

$$ \begin{align} R(\phi) &= \begin{bmatrix} 1 & 0 & 0 \ 0 & c\phi & -s\phi \ 0 & s\phi & c\phi \end{bmatrix} \ R(\theta) &= \begin{bmatrix} c\theta & 0 & s\theta \ 0 & 1 & 0 \ -s\theta & 0 & c\theta \end{bmatrix} \ R(\psi) &= \begin{bmatrix} c\psi & -s\psi & 0 \ s\psi & c\psi & 0 \ 0 & 0 & 1 \end{bmatrix} \end{align} $$

輸送定理 (Transport Theorem)

参照系が回転している場合、ベクトルの変化率は単純ではありません。瞬間的な角速度 $\hat{\omega}$ を含む以下の式（輸送定理）を用います。

$$ \frac{d \hat{\mathcal{V}}^b}{d t} = \hat{b} \cdot \frac{d \hat{\mathcal{V}}}{d t} + \hat{\omega} \times \hat{\mathcal{V}}^b $$

これは、体固定系での時間微分と**回転による擬似項（外積）**の和として表現されます。これにより、非慣性系でもニュートンの法則を適用できるようになります。

---

2. 飛行器の状態変数

飛行器のダイナミクスは以下の 12 つの状態変数によって記述されます。

• $\hat{r}^n=[x,y,z]^T$: 慣性系のナビゲーション座標（位置）
• $\hat{v^b}=[v_x, v_y, v_z]^T$: 体固定系軸方向の速度
• $\hat{\Phi}=[\phi, \theta, \psi]^T$: オイラー角（姿勢：ロール、ピッチ、ヨー）
• $\hat{\omega}=[\omega_x, \omega_y, \omega_z]^T$: 体固定系の瞬間的な角速度

---

3. マルチローター UAV のダイナミクス

状態変数の時間変化は、線形運動方程式と角運動方程式に分割されます。ここでは力およびモーメントを体固定系の座標 ($b$) で表現します。

3.1 力のモデル

飛行器に加えられている主な力は以下の通りです：

• 重力: 慣性系 $\hat{F_g}^n=[0,0,-mg]^T$（注意：原文の記号定義に合わせて整合させる必要がありますが、物理的には下方）。体固定系では回転行列で変換されます。
• 推力 ($T$): プロペラにより生じ、体固定系の $z$ 軸方向にのみ働きます。

$$ \hat{F}^b = \hat{T} + R_n^b \hat{F_g}^n $$

3.2 線形運動方程式 (Newton's 2nd Law for Translation)

体固定系での運動量は時間変化するため、輸送定理を適用します。

$$ \hat{F}^b = m(\hat{\dot{v}}^b + \hat{\omega} \times \hat{v}^b) $$

これより、加速度 $\hat{\dot{v}}^b$ を導出すると：

$$ \hat{\dot{v}}^b = \frac{\hat{T} + R_n^b \hat{F_g}^n}{m} - \hat{\omega} \times \hat{v}^b $$

各成分展開形（簡易化された重力成分を考慮）：

$$ \begin{bmatrix} \dot{v}_x \ \dot{v}_y \ \dot{v}_z \end{bmatrix} = \begin{bmatrix}

• \omega_{y} v_{z} + \omega_{z} v_{y} + g \sin{\theta}\ \omega_{x} v_{z} - \omega_{z} v_{x} - g \sin{\phi} \cos{\theta}\ -\frac{T}{m} - \omega_{x} v_{y} + \omega_{y} v_{x} + g\cos{\phi} \cos{\theta} \end{bmatrix} $$

位置の変化率は、速度ベクトルを回転行列で変換することで得られます： $$ \begin{bmatrix} \dot{x} \ \dot{y} \ \dot{z} \end{bmatrix} = R_b^n \begin{bmatrix} v_x & v_y & v_z \end{bmatrix}^T $$

3.3 角運動方程式 (Newton's 2nd Law for Rotation)

トルク ($\tau$) は角運動量の時間変化率です。

$$ \hat{\tau} = \frac{d}{dt}(I\hat{\omega}) $$

慣性行列 ($I$): ドローンを剛体とみなす場合、各軸周りの対角要素 $I_{xx}, I_{yy}, I_{zz}$ (慣性モーメント) と非対角要素 (慣性積) で構成されます。対称性を仮定すると：

$$ I = \begin{bmatrix} I_{xx} & 0 & 0 \ 0 & I_{yy} & 0 \ 0 & 0 & I_{zz} \end{bmatrix} $$

慣性モーメントの計算:

• 実体球 (重心): $\frac{2}{5} m R^2$
• 実体円柱 (対称軸): $\frac{1}{2} m R^2$
• 実体円柱 (垂直軸): $\frac{1}{4} m R^2 + \frac{1}{12} m L^2$
• 点質量: $m r^2$

平行軸の定理を使用して、各部品の慣性モーメントを合計します。

3.4 輸送定理による角運動量の展開

体固定系では $\hat{\tau}^b$ が定義されるため、以下の関係式が成立します。

$$ \hat{\tau}^b = I\frac{d}{dt}\hat{\omega} + \hat{\omega} \times I\hat{\omega} $$

ここで、$\frac{d}{dt}$ は体固定系内での変化率です。これより角加速度 $\dot{\hat{\omega}}$ を求めます：

$$ \dot{\hat{\omega}} = I^{-1} (\hat{\tau}^b - \hat{\omega} \times I\hat{\omega}) $$

トルクの構成

総トルク $\hat{\tau}^b$ は 2 つの成分で構成されます。

1. 推力によるモーメント ($\hat{\tau}_T^b$): プロペラ位置ベクトル $r_i$ と推力 $T_i$ の外積 ($r_i \times T_i$)。これによりロール・ピッチが生じます。
2. ヨーモーメント ($\hat{\tau}_\tau^b$): モーターのトルク $\tau_i$ の和（反作用トルク）。これによりヨーが生じます。

$$ \begin{align} \dot{\phi} &= \omega_x - \frac{\omega_z \sin{\theta}}{\cos{\phi} \cos{\theta}} \ \dot{\theta} &= \omega_y + \frac{\omega_z \sin{\phi}}{\cos{\phi}} \ \dot{\psi} &= \frac{\omega_z}{\cos{\phi} \cos{\theta}} \end{align} $$

これら 12 つの状態変数の微分方程式は、シミュレーションの積分ステップにおいて連続的に更新されます。

---

4. 力およびトルクの発生メカニズム

4.1 プロペラ (Blade Elements)

プロペラが回転すると、空気へ力が及び、その反作用として推力と抵抗トルクが発生します。

• 推力 ($T$): プロペラ速度（回転数）$\Omega$ の 2 乗に比例します。 $$ T = k_{T} \cdot \Omega^2 $$
• 誘起トルク ($\tau$): 同様に $\Omega^2$ に比例します。 $$ \tau = k_{d} \cdot \Omega^2 $$

$k_T, k_d$ はプロペラ定数です（形状とピッチに依存）。

4.2 モーター (BLDC)

ブラシレス直流モーターは以下のモデルで記述されます。

• 電流 ($i_{BLDC}$): $$ i_{BLDC} = \frac{v_{BLDC} - k_e \cdot \Omega}{R_{BLDC}} $$ ($k_e$: 逆起電力定数, $R_{BLDC}$: 内部抵抗)
• 発生トルク ($\tau_{BLDC}$): $$ \tau_{BLDC} = k_\tau \cdot i_{BLDC} $$ (理想では $k_\tau \approx k_e$)
• 総トルク ($\tau$): 発生トルクから消耗（摩擦、空気抵抗）を引きます。 $$ \tau = \tau_{BLDC} - k_{DF}\cdot\Omega - k_d\cdot\Omega^2 $$
• 角速度変化: $$ \frac{d \Omega}{dt} = \frac{\tau}{J} $$

---

5. コントロールシステム (Control System)

5.1 コントローラー配分 (Control Allocation)

望ましいダイナミクス（力・トルク）を実現するために、各モーターの回転速度 $\hat{\Omega}$ を決定します。

プロペラ $p$ の個数に基づき、以下のような行列関係 $A \cdot \hat{\Omega}^2 = D$ が成り立ちます。

• $D$: 必要な力・トルクベクトル (6x1)
• $A$: 配分行列 (6xp)，プロペラの位置と推力定数に基づき生成されます。

対称的な配置において、$A$ は正方でない場合が多いです。この場合は**疑似逆行列 (Pseudoinverse)**を使用します。

$$ \hat{\Omega} = \sqrt{A^+ D}, \quad A^+ = A^T (A A^T)^{-1} $$

数値的安定性の観点から、特異値分解 (SVD) を使用して疑似逆行列を計算する方法が推奨されます。

5.2 参照追跡 (Reference Tracking)

ウェイポイント $\hat{r}_{wp}$ に従うための制御を行います。PID コントローラーが標準的に使用されます。

• 位置誤差 $\rightarrow$ 目的速度指令
• 速度誤差 $\rightarrow$ 目的姿勢（ロール・ピッチ）指令
• 姿勢誤差 $\rightarrow$ 目的角速度指令

カスケード PID コントローラー構成

1. 外ループ: 位置追跡を行い、必要な推力と姿勢角度を計算します。
2. 中ループ: 姿勢コントローラー（ロール・ピッチ・ヨー）が入力を受け、必要なトルクを生成します。
3. 内ループ: モーター制御がトルク指令に対して角速度 $\Omega$ を調整します。

生成されたトルク指令は、前述のコントローラー配分を通じてモーター速度にマッピングされます。

5.3 パス計画および監視制御

より上位層の機能として、以下のタスクがあります。

• ウェイポイント選択: 飛行エンベロープ（最大速度、加速度制約など）を考慮した軌道設計。
• 軌道生成: S 曲線や B スプラインを使用。
• 最適化: モデル予測制御 (MPC) や強化学習を用いた監視制御。

これらは本稿の範囲外ですが、ArduPilot などのフレームワークで実装可能です。