航空機のラダーのみ旋回の 6 自由度剛体運動をシミュレーションする Python スクリプト

この記事では、ラダーのみ旋回に対する 6 自由度剛体運動をシミュレーションする Python スクリプトについて説明する。

はじめに

この記事では、ラダーのみ旋回に対する 6 自由度剛体運動をシミュレーションするPythonスクリプトについて説明する。

ラダーのみ旋回では、ラダー入力によってヨー運動を作り、その結果として横滑り角とヨーレートが生じる。

横滑り角とヨーレートは、上反角効果、ヨーレートロール、ロール減衰などを通じてロール運動につながる。

↓ラダーのみを用いて行う定常旋回についてはこちら

航空機のラダーのみ定常旋回

ラダーのみを用いた定常旋回について説明する。

mtkbirdman.com

2026.07.11

↓ラダーのみ旋回の横転性能についてはこちら

航空機のラダーのみで旋回における横転性能

ラダーのみで旋回する航空機の横転性能について説明する。

mtkbirdman.com

2026.07.12

本記事では、それらの定式化を使い、VSPAERO の .stab ファイルから得られる線形空力モデルを 6 自由度剛体運動方程式に入れて、ラダーステップ入力に対する時間応答を計算する。

最終的にはこんな感じになる

それではいってみよう。

プログラムの全体像

今回の中心となるファイルは RollRudderGain.py である。

このファイルでは、ラダー入力に対する横方向・方向運動を、次の 2 段階で評価する。

1. 線形横・方向モデルで有限時間応答を計算する。
2. 6 自由度剛体運動方程式で有限時間応答を計算する。

線形横・方向モデルでは、バンク角を含めた 4 状態を扱う。

\begin{align} \mathbf{x} =\begin{bmatrix} \beta & \hat{p} & \hat{r} & \phi \end{bmatrix}^{\mathrm{T}} \end{align}

6 自由度剛体運動では、機体軸速度、機体軸角速度、地球固定座標系の位置、姿勢角を時間積分する。

スクリプトが積分する状態量は、次の 12 個である。

\begin{align} \mathbf{x}_{6\mathrm{DOF}} =\begin{bmatrix} u & v & w & p & q & r & x_e & y_e & z_e & \phi & \theta & \psi \end{bmatrix}^{\mathrm T} \end{align}

ここで、\(u,v,w\) は機体軸速度、\(p,q,r\) は機体軸角速度、\(x_e,y_e,z_e\) は地球固定座標系の位置、\(\phi,\theta,\psi\) はオイラー角である。

処理の流れは次のようになる。

1. .stab ファイルを読み込む。
2. Control Surface Group からラダー舵角列を特定する。
3. 基準飛行条件、基準空力係数、安定微係数を取得する。
4. ラダーステップ入力 \(\delta_r\) を設定する。
5. 初期状態を作る。
6. 線形空力モデルで各時刻の空力係数を計算する。
7. 6 自由度剛体運動方程式の右辺を作る。
8. scipy.integrate.solve_ivp() で時間積分する。
9. バンク角変化、平均バンク角速度、有限時間ロール応答指標を後処理する。

それでは 1 つずつ順番に説明する。

入力データ

.stab ファイル

入力として使う .stab ファイルは、OpenVSP / VSPAERO の安定微係数解析で得られるファイルである。

.stab には、代表的に次の情報が含まれる。

• 参照面積 \(S_{\mathrm{ref}}\)
• 参照翼弦 \(c_{\mathrm{ref}}\)
• 参照翼幅 \(b_{\mathrm{ref}}\)
• 基準迎角
• 基準横滑り角
• 基準飛行速度
• 基準密度
• 基準空力係数
• 迎角、横滑り角、角速度、舵角に対する安定微係数

RollRudderGain.py では、read_vspaero_stab() と resolve_control_columns_from_stab() を使う。

これにより、.stab ファイルの読み取りと、Control Surface Group の対応付けを共通化している。

質量・慣性モーメント

6 自由度剛体運動では、質量と慣性モーメントが必要である。

主に使う値は次である。

\begin{align} m, I_{xx}, I_{yy}, I_{zz}, I_{xz} \end{align}

ここで、\(m\) は機体質量、\(I_{xx}\)、\(I_{yy}\)、\(I_{zz}\) は重心まわりの慣性モーメント、\(I_{xz}\) は慣性乗積である。

ラダー入力

ラダー入力は、基本的にはステップ入力として与える。

\begin{align} \delta_r(t) =\begin{cases} 0, & t<0, \\ \delta_{r,0}, & t\ge 0 \end{cases} \end{align}

ここで、\(\delta_{r,0}\) はステップ後のラダー舵角である。

スクリプトでは rad 単位で与える。

import math

delta_r = math.radians(-10.0)

プログラムの解説

関数構成

RollRudderGain.py の主な関数は次である。

関数	役割
calculate_linear_lateral_response_indices_from_stab()	線形横・方向モデルで有限時間応答を計算する
simulate_6dof_rudder_step_from_stab()	ラダーステップ入力に対する 6 自由度応答を計算する
write_6dof_history_csv()	6 自由度応答の時刻歴を CSV に保存する
plot_6dof_history()	6 自由度応答の時刻歴を描く
simulate_reduced_lateral_response_from_stab()	簡約横・方向モデルの時刻歴を計算する
estimate_roll_rate_gain_from_history()	応答時刻歴からロール率ゲインを推算する
calculate_roll_response_index_by_delta_phi()	指定したバンク角変化に到達するまでの有限時間指標を計算する

線形横・方向モデル

線形横・方向モデルでは、状態量を

\begin{align} \mathbf{x} =\begin{bmatrix} \beta & \hat{p} & \hat{r} & \phi \end{bmatrix}^{\mathrm{T}} \end{align}

とする。

無次元時間は

\begin{align} \tau = \frac{2V}{b}t \end{align}

である。

このとき、ラダーステップ入力に対する応答は、概念的に

\begin{align} \frac{d\mathbf{x}}{d\tau}=\mathbf{A}\mathbf{x}+\mathbf{B}\delta_r \end{align}

で表せる。

ここで、\(\mathbf{A}\) は横滑り角、ロール率、ヨーレート、バンク角の連成を表す行列であり、\(\mathbf{B}\) はラダー入力の入り方を表す列ベクトルである。

calculate_linear_lateral_response_indices_from_stab() は、この線形モデルを時間積分し、指定したバンク角変化 \(\Delta\phi\) に到達するまでの時間とロール応答指標を返す。

6 自由度剛体運動方程式

6 自由度剛体運動では、機体軸の並進運動と回転運動を同時に解く。

速度の式は、機体軸を \(+x_b\) 前方、\(+y_b\) 右方、\(+z_b\) 下方として、

\begin{align} \dot{u} &= rv-qw-g\sin\theta+\frac{X}{m}, \\ \dot{v} &= pw-ru+g\sin\phi\cos\theta+\frac{Y}{m}, \\ \dot{w} &= qu-pv+g\cos\phi\cos\theta+\frac{Z}{m} \end{align}

である。

角速度の式は、慣性乗積 \(I_{xz}\) を含めて、

\begin{align} \begin{bmatrix} I_{xx} & -I_{xz} \\ -I_{xz} & I_{zz} \end{bmatrix} \begin{bmatrix} \dot{p} \\ \dot{r} \end{bmatrix} &= \begin{bmatrix} L+(I_{yy}-I_{zz})qr+I_{xz}pq \\ N+(I_{xx}-I_{yy})pq-I_{xz}qr \end{bmatrix}, \\\\
I_{yy}\dot{q} &=M+(I_{zz}-I_{xx})pr-I_{xz}(p^2-r^2) \end{align}

である。

補助的に 2 行 2 列の逆行列を用いれば、\(\dot{p}\) と \(\dot{r}\) は次のように完全に書ける。

\begin{align} \dot{p} &= \frac{ I_{zz}\left[L+(I_{yy}-I_{zz})qr+I_{xz}pq\right] +I_{xz}\left[N+(I_{xx}-I_{yy})pq-I_{xz}qr\right] }{I_{xx}I_{zz}-I_{xz}^2}, \\\\
\dot{r} &= \frac{ I_{xz}\left[L+(I_{yy}-I_{zz})qr+I_{xz}pq\right] +I_{xx}\left[N+(I_{xx}-I_{yy})pq-I_{xz}qr\right] }{I_{xx}I_{zz}-I_{xz}^2} \end{align}

オイラー角の運動学は

\begin{align} \dot{\phi} &= p+q\sin\phi\tan\theta+r\cos\phi\tan\theta, \\ \dot{\theta} &= q\cos\phi-r\sin\phi, \\ \dot{\psi} &= \frac{q\sin\phi+r\cos\phi}{\cos\theta} \end{align}

である。

simulate_6dof_rudder_step_from_stab() は、この右辺を作り、solve_ivp() で時間積分する。

空力係数の作り方

6 自由度応答中の空力係数は、.stab の基準値と安定微係数を使って一次近似する。

迎角と横滑り角は、速度成分から

\begin{align} V &= \sqrt{u^2+v^2+w^2}, \\ \alpha &= \tan^{-1}\left(\frac{w}{u}\right), \\ \beta &= \sin^{-1}\left(\frac{v}{V}\right) \end{align}

として計算する。

無次元角速度は

\begin{align} \hat{p} &= \frac{pb}{2V}, \\ \hat{q} &= \frac{qc}{2V}, \\ \hat{r} &= \frac{rb}{2V} \end{align}

である。

まず、VSPAERO の .stab に含まれる風軸系の力係数 \(C_L,C_D,C_S\) と機体軸モーメント係数 \(C_l,C_m,C_n\) を、基準状態まわりで一次近似する。例えば、ローリングモーメント係数は

\begin{align} C_l =C_{l0} +C_{l\alpha}\Delta\alpha +C_{l\beta}\Delta\beta +C_{l\hat{p}}\hat{p} +C_{l\hat{q}}\hat{q} +C_{l\hat{r}}\hat{r} +C_{l\delta_a}\delta_a +C_{l\delta_e}\delta_e +C_{l\delta_r}\delta_r \end{align}

である。同じ方法で \(C_L,C_D,C_S,C_m,C_n\) を計算する。

6 自由度方程式へ入れる力係数には、.stab の CFx、CFy、CFz を直接使用しない。計算した \(C_L,C_D,C_S\) を、現在の \(\alpha,\beta\) を使って、ソルバーの機体軸 \(+x_b\) 前方、\(+y_b\) 右方、\(+z_b\) 下方へ変換する。

\begin{align}
C_X&=-\left(C_D\cos\beta+C_S\sin\beta\right)\cos\alpha+C_L\sin\alpha, \\
C_Y&=-C_D\sin\beta+C_S\cos\beta, \\
C_Z&=-\left(C_D\cos\beta+C_S\sin\beta\right)\sin\alpha-C_L\cos\alpha
\end{align}

この変換後の係数から、空気力と空力モーメントを

\begin{align} X &= q_\infty S C_X, \\ Y &= q_\infty S C_Y, \\ Z &= q_\infty S C_Z, \\ L &= q_\infty S b C_l, \\ M &= q_\infty S c C_m, \\ N &= q_\infty S b C_n \end{align}

として 6 自由度方程式に入れる。

ここで、動圧は

\begin{align} q_\infty = \frac{1}{2}\rho V^2 \end{align}

である。

ラダーステップ入力と初期エレベータ

現在の simulate_6dof_rudder_step_from_stab() は、.stab の基準速度、基準迎角、基準横滑り角から初期状態を作り、時刻 \(t=0\) 以降に指定されたラダー舵角を一定値として与える。

\begin{align} \delta_r(t) =\begin{cases} 0, & t<0, \\ \delta_{r,0}, & t\ge 0 \end{cases} \end{align}

delta_e=None かつ trim_elevator=True の場合、初期エレベータ舵角は、時刻 \(t=0\) 以降に使用するラダー舵角 \(\delta_{r,0}\) をすでに含めた空力係数から、初期ピッチングモーメント係数がゼロになるように計算される。

\begin{align} C_m\left(\alpha_0,\beta_0,p_0,q_0,r_0, \delta_a,\delta_e,\delta_{r,0}\right)=0 \end{align}

theta_hold=False の場合、この初期エレベータ舵角は時間積分中に固定される。

したがって、これは「ラダーだけを変化させ、他の舵角を入力前の値に保つ純粋な単独ステップ」とは異なる。

また、初期エレベータと初期推力の調整は初期時刻の釣り合いを改善するためのものであり、その後も定常飛行が維持されることを保証するものではない。

有限時間ロール応答指標

ラダー入力直後のロール能力は、ある時刻の瞬間的な機体軸ロール角速度 \(p\) だけでは判断しにくい。

そこで、指定したバンク角変化に到達するまでの平均バンク角速度として評価する。

初期バンク角を \(\phi_0\)、目標バンク角変化を \(\Delta\phi\) とすれば、到達時刻 \(t_{\Delta\phi}\) は

\begin{align} \phi(t_{\Delta\phi})-\phi_0=\Delta\phi \end{align}

で定義される。

このとき、目標到達までの平均バンク角速度は

\begin{align} \overline{\dot{\phi}}_{\Delta\phi} =\frac{\Delta\phi}{t_{\Delta\phi}} \end{align}

である。

平均バンク角速度を \(b/(2V)\) で無次元化し、ラダー舵角で割った有限時間ロール応答指標を

\begin{align} K_{\phi,\delta_r} =\frac{b}{2V} \frac{ \overline{\dot{\phi}}_{\Delta\phi} }{ \delta_r } =\frac{b}{2V} \frac{ \Delta\phi/t_{\Delta\phi} }{ \delta_r } \end{align}

と定義する。

calculate_roll_response_index_by_delta_phi() は、時刻歴から目標バンク角へ最初に到達する時刻を線形補間し、この指標を計算する。

使い方

ディレクトリ構成

最小構成は次のようになる。

project/
├── RollRudderGain.py
├── TrimTurnSolver.py
├── SampleGlider.stab
└── run_6dof_rudder.py

.stab ファイルは、OpenVSP / VSPAERO の安定微係数解析で作成しておく。

入力する値

6 自由度応答の計算で必要になる主な入力は次である。

引数	意味
stab_path	VSPAERO の .stab ファイル
mass	機体質量
Ixx	ロール軸まわりの慣性モーメント
Iyy	ピッチ軸まわりの慣性モーメント
Izz	ヨー軸まわりの慣性モーメント
Ixz	慣性乗積
delta_r	ラダーステップ入力
t_final	積分終了時刻
rho	密度。省略時は .stab の Rho を使う
control_map	Control Group と舵角の対応
target_delta_phi	有限時間ロール応答を評価するバンク角変化

角度と舵角は rad で指定する。

import math

delta_r = math.radians(-10.0)
target_delta_phi = math.radians(2.0)

速度、密度、質量、面積、慣性モーメントは、一貫した単位系で指定する。

G103A の SI 単位系の例では、質量と密度を次のように明示する。

mass = 580.0
rho = 1.225

.stab に保存された密度が別の単位系である場合は、rho を省略しない。

慣性モーメントには、同じ質量状態・同じ単位系で求めた値を使用する。

実行例：6 自由度ラダーステップ応答

example_rudder_step_response.ipynb と同じ SampleGlider.stab、質量、慣性モーメント、ラダー舵角を用いる例は次の通りである。

import math

from RollRudderGain import simulate_6dof_rudder_step_from_stab

history = simulate_6dof_rudder_step_from_stab(
    stab_path="SampleGlider.stab",
    mass=100.0,
    Ixx=1000.0,
    Iyy=75.0,
    Izz=1100.0,
    Ixz=0.0,
    delta_r=math.radians(-10.0),
    t_final=5.0,
    delta_a=0.0,
    delta_e=None,
    trim_elevator=True,
    trim_thrust=False,
    phi0=0.0,
    theta0=None,
    psi0=0.0,
    max_step=0.02,
)

print(history.columns)
print(history.tail())

この例では推力を使用せず、初期エレベータは \(\delta_r=-10\ \mathrm{deg}\) を含む初期ピッチングモーメントがゼロになるように計算される。

返り値は pandas.DataFrame であり、12 状態と空力・制御量の時刻歴が列として入る。

実行例：時刻歴を CSV へ保存する

6 自由度応答を CSV に保存する場合は、write_6dof_history_csv() を使う。

from RollRudderGain import write_6dof_history_csv

write_6dof_history_csv(
    history,
    "rudder_step_6dof_history.csv",
)

これにより、時刻、速度成分、角速度、姿勢角、迎角、横滑り角、舵角などを後で確認できる。

実行例：時刻歴を描く

時刻歴を確認する場合は、plot_6dof_history() を使う。

from RollRudderGain import plot_6dof_history

fig = plot_6dof_history(
    history,
    plot_path="rudder_step_6dof_history.png",
    show=True,
    degrees=True,
)

ラダー入力後に、横滑り角、ヨーレート、ロール率、バンク角がどのように立ち上がるかを見る。

実行例：有限時間ロール応答指標を計算する

指定したバンク角変化に到達するまでの有限時間ロール応答指標を計算する場合は、calculate_roll_response_index_by_delta_phi() を使う。

import math

from RollRudderGain import calculate_roll_response_index_by_delta_phi
from TrimTurnSolver import read_vspaero_stab

stab = read_vspaero_stab("SampleGlider.stab")

metric = calculate_roll_response_index_by_delta_phi(
    history,
    target_delta_phi=math.radians(2.0),
    delta_r=math.radians(-10.0),
    V=stab.V0,
    Bref=stab.Bref,
    phi0=0.0,
)

print(metric)

この結果から、目標バンク角変化に到達するまでの時間、平均バンク角速度、ラダー舵角あたりの平均バンク角速度、有限時間ロール応答指標を確認する。

Control Group の対応を明示する場合

.stab の Control Group 名が標準的でない場合は、control_map を明示する。

control_map = {
    "delta_a": "ConGrp_1",
    "delta_e": "ConGrp_2",
    "delta_r": "ConGrp_3",
}

ラダー応答だけを計算する場合でも、初期トリムや線形空力モデルでエルロン、エレベータの列を使うことがある。

Control Group の対応は、最初に確認しておくとよい。

出力データ

6 自由度時刻歴

6 自由度応答では、主に次の時刻歴を見る。

• 時刻 \(t\)
• 速度成分 \(u,v,w\)
• 角速度 \(p,q,r\)
• 位置 \(x_e,y_e,z_e\)
• 姿勢角 \(\phi,\theta,\psi\)
• 迎角 \(\alpha\)
• 横滑り角 \(\beta\)
• ラダー舵角 \(\delta_r\)
• バンク角変化 \(\Delta\phi\)

これらを見れば、ラダー入力からヨー運動が始まり、横滑り角とヨーレートを経由してロール運動が立ち上がる流れを確認できる。

応答指標

有限時間応答では、主に次の量を見る。

キー	意味
sixdof_roll_response_reached	指定した符号付きバンク角変化へ到達したか
sixdof_t_reach	目標バンク角変化へ到達した時刻
sixdof_dt_reach	初期時刻から目標到達までの所要時間
sixdof_roll_response_phi_rate	目標到達までの平均バンク角速度
sixdof_roll_response_phi_rate_per_delta_r	ラダー舵角あたりの平均バンク角速度
sixdof_finite_time_roll_index	平均バンク角速度を \(b/(2V)\) で無次元化した有限時間ロール応答指標
sixdof_roll_response_error	目標未到達などの後処理状態

sixdof_roll_response_reached=False の場合は、指定した t_final の範囲内で目標バンク角変化に到達していないことを意味する。

この場合は、t_final を伸ばすか、target_delta_phi を小さくして確認する。

おわりに

この記事では、ラダーのみ旋回に対する 6 自由度剛体運動をシミュレーションする Python スクリプトについて説明した。