航空機が横方向の 1-cosine gust を受けたときの 6 自由度剛体運動シミュレーションを行う Python スクリプトについて説明する。
はじめに
この記事では、航空機が横方向の 1-cosine gust を受けたときの 6 自由度剛体運動を Python で時間積分するスクリプトについて説明する。
横風突風を受けると、機体から見た相対風の横方向成分が変化する。
微小擾乱として見れば、これは外乱横滑り角として空力係数へ入る。
横滑り角の外乱は、横力、ローリングモーメント、ヨーイングモーメントを通じて、横滑り角、ヨーレート、ロール率、バンク角を励起する。
横風突風応答の定式化、1-cosine gust の意味、線形横・方向モデル、Taylor 展開による短時間応答の詳細は、以下の記事を参照。
本記事では、その定式化を実装側から整理する。
具体的には、VSPAERO の .stab ファイルから読み取った基準空力係数と安定微係数を使い、横風突風速度を相対風へ入れながら 6 自由度剛体運動方程式を時間積分する。
最終的にはこんな感じになる
それではいってみよう。
プログラムの全体像
今回の中心となるファイルは RollRudderGain.py である。
このファイルはもともと、ラダー入力に対する線形横・方向応答と 6 自由度剛体運動応答を比較するための関数群を含んでいる。
その中に、横風突風応答を計算する関数もまとめて実装している。
横風突風応答に直接関係する主な関数は次である。
| 関数 | 役割 |
|---|---|
one_cosine_gust_velocity() | 指定時刻における 1-cosine gust の有次元速度を返す |
simulate_6dof_crosswind_gust_from_stab() | .stab 線形空力モデルを使い、横風突風に対する 6 自由度応答を積分する |
calculate_crosswind_gust_response_indices() | 横風突風応答の時刻歴からピーク値や評価指標を計算する |
write_6dof_history_csv() | 6 自由度応答の時刻歴を CSV に保存する |
plot_6dof_history() | 6 自由度応答の時刻歴を描く |
処理の流れは次のようになる。
- .stab ファイルを読み込む。
- Control Surface Group と delta_a、delta_e、delta_r の対応を決める。
- .stab の基準飛行条件から初期状態を作る。
- 必要に応じて初期エレベータと推力をトリムする。
- 時刻ごとに 1-cosine gust の横風速度を計算する。
- 横風速度を機体軸の突風速度として相対風へ入れる。
- 相対風にもとづいて線形空力係数を評価する。
- 空気力、空力モーメント、推力、重力を使って 6 自由度方程式の右辺を作る。
solve_ivp()で時間積分する。- 時刻歴を
pandas.DataFrameとして返す。 - 必要に応じて時刻歴 CSV、図、ピーク応答指標を作る。
図にすると、概念的には次の流れである。
.stab
↓
基準条件・安定微係数・Control Group 対応の読み取り
↓
初期状態・舵角・推力の設定
↓
1-cosine gust U_g(t) の生成
↓
相対風速度 u_a, v_a, w_a の計算
↓
alpha, beta, p_hat, q_hat, r_hat の計算
↓
.stab 線形空力モデルで C_X, C_Y, C_Z, C_l, C_m, C_n を評価
↓
6 自由度剛体運動方程式を時間積分
↓
時刻歴・ピーク値・応答指標の出力このように、本スクリプトは新しい空力モデルを作るものではなく、VSPAERO の .stab 線形空力モデルを 6 自由度運動方程式へ接続するためのものである。
入力データ
.stab ファイル
入力として使う .stab ファイルは、OpenVSP / VSPAERO の STABILITY_DEFAULT 解析で出力される安定微係数ファイルである。
.stab には、代表的に次の情報が含まれる。
- 参照面積 \(S_{\mathrm{ref}}\)
- 参照翼弦 \(c_{\mathrm{ref}}\)
- 参照翼幅 \(b_{\mathrm{ref}}\)
- 重心位置
- Mach 数
- 基準迎角
- 基準横滑り角
- 密度
- 飛行速度
- 基準空力係数
- 迎角、横滑り角、角速度、舵角に対する安定微係数
- Control Surface Group と ConGrp_* の対応
RollRudderGain.py では、.stab の読み取りに read_vspaero_stab() を使う。
また、Control Surface Group と舵角の対応付けには resolve_control_columns_from_stab() を使う。
質量・慣性モーメント
6 自由度剛体運動では、質量と慣性モーメントが必要である。
主に使う値は次である。
\begin{align} m, I_{xx}, I_{yy}, I_{zz}, I_{xz} \end{align}
ここで、\(m\) は機体質量、\(I_{xx}\)、\(I_{yy}\)、\(I_{zz}\) は重心まわりの慣性モーメント、\(I_{xz}\) は慣性乗積である。
.stab の長さ、速度、密度と、ここで与える質量・慣性モーメントは同じ単位系でそろえる。
OpenVSP / VSPAERO の .stab は Lunit、Tunit、Munit として値を出すため、Python 側で暗黙の単位変換は行わない。
横風突風条件
横風突風は、1-cosine gust として与える。
突風の有次元速度を \(U_g(t)\) とすると、
\begin{align} U_g(t) = \frac{U_{ds}}{2}\left[1-\cos\left(\frac{\pi Vt}{H}\right)\right],\quad 0 \le t \le \frac{2H}{V} \end{align}
である。
ここで、RollRudderGain.py に与える \(U_{ds}\) は、解析高度における符号付きの有次元ピーク突風速度であり、\(V\) と同じ長さ・時間単位で与える。
\(U_{ds}>0\) は右から左へ吹く横風突風を表す。
14 CFR 25.341 の設計突風速度を用いる場合、規定値は等価対気速度であるため、解析高度における真の突風速度へ変換する。
\begin{align} U_{ds,\mathrm{TAS}} =U_{ds,\mathrm{EAS}} \sqrt{ \frac{\rho_0}{\rho} } \end{align}
\(H\) は突風勾配距離である。
突風場を飛行速度 \(V\) で通過するとすれば、突風通過時間は、
\begin{align} t_e = \frac{2H}{V} \end{align}
である。
RollRudderGain.py の one_cosine_gust_velocity() は、この \(U_g(t)\) を返す。
突風通過後は、
\begin{align} U_g(t) = 0 \end{align}
として扱う。
舵角・初期状態・トリム条件
横風突風応答では、基本的には舵角を固定した状態で突風を入れる。
主な舵角は、
\begin{align} \delta_a, \delta_e, \delta_r \end{align}
である。
通常は、
\begin{align} \delta_a &= 0, \\ \delta_r &= 0 \end{align}
として、横風突風だけに対する機体応答を見る。
エレベータ舵角 \(\delta_e\) は、.stab の基準状態で初期ピッチングモーメントを小さくするため、必要に応じて自動トリムする。
推力も同様に、初期の機体軸 \(x\) 方向の力を釣り合わせるため、必要に応じて自動トリムする。
プログラムの解説
関数構成
RollRudderGain.py の横風突風応答に関係する関数は次である。
| 関数 | 説明 |
|---|---|
_initial_state_from_stab() | .stab の基準速度、迎角、横滑り角から 6 自由度初期状態を作る |
_linear_aero_from_stab() | 状態量、舵角、突風速度から線形空力係数を計算する |
one_cosine_gust_velocity() | 指定時刻における 1-cosine gust の速度を計算する |
_resolve_initial_controls() | 初期舵角を決める。必要に応じてエレベータトリムを行う |
_resolve_initial_thrust() | 初期推力を決める。必要に応じて機体軸 \\(x\\) 方向の力を釣り合わせる |
simulate_6dof_crosswind_gust_from_stab() | 横風突風に対する 6 自由度応答を時間積分する本体 |
calculate_crosswind_gust_response_indices() | 時刻歴からピーク横滑り角、ピークロール率、ピークバンク角などを計算する |
実際に使うときは、基本的には simulate_6dof_crosswind_gust_from_stab() を呼び出せばよい。
その後、必要に応じて calculate_crosswind_gust_response_indices()、write_6dof_history_csv()、plot_6dof_history() を使う。
1-cosine gust
1-cosine gust は、次の式で与える。
\begin{align} U_g(t) = \frac{U_{ds}}{2}\left[1-\cos\left(\frac{\pi V(t-t_0)}{H}\right)\right],\quad 0 \le t-t_0 \le \frac{2H}{V} \end{align}
ここで、\(t_0\) は突風開始時刻である。
突風開始前と突風通過後は、
\begin{align} U_g(t) = 0 \end{align}
である。
one_cosine_gust_velocity() は、概念的には次の処理を行う。
def one_cosine_gust_velocity(time, *, Uds, H, V, start_time=0.0):
s = V * (time - start_time)
if s < 0.0 or s > 2.0 * H:
return 0.0
return 0.5 * Uds * (1.0 - math.cos(math.pi * s / H))ここで、\(s=V(t-t_0)\) は突風場へ入ってからの進行距離である。
したがって、\(s=H\) で突風速度は最大値 \(U_{ds}\) になり、\(s=2H\) でゼロに戻る。
横風突風を相対風として入れる処理
機体軸速度を、
\begin{align} u, v, w \end{align}
とする。
ここで、機体軸は \(x_b\) 前方、\(y_b\) 右方、\(z_b\) 下方である。
本スクリプトでは、\(U_g(t)>0\) を右から左へ吹く横風突風として扱う。
機体軸は \(+y_b\) 方向を右方とするため、実際の大気速度ベクトルは、
\begin{align} \mathbf{v}_{\mathrm{airmass}} =\begin{bmatrix} 0 \\ -U_g(t) \\ 0 \end{bmatrix} \end{align}
である。
空力的に効くのは、機体速度そのものではなく、機体が感じる相対風である。
したがって、空力計算に使う相対速度成分は、
\begin{align} u_a &= u, \\ v_a &= v-\left[-U_g(t)\right]=v+U_g(t), \\ w_a &= w \end{align}
となる。
このとき、空力計算に使う対気速度、迎角、横滑り角は、
\begin{align} V_a &= \sqrt{u_a^2+v_a^2+w_a^2}, \\ \alpha_a &= \tan^{-1}\left(\frac{w_a}{u_a}\right), \\ \beta_a &= \sin^{-1}\left(\frac{v_a}{V_a}\right) \end{align}
である。
小擾乱として見れば、突風による外乱横滑り角は、
\begin{align} \beta_g \simeq \frac{U_g}{V} \end{align}
であり、有効横滑り角は、
\begin{align} \beta_a \simeq \beta+\beta_g \end{align}
である。
6 自由度実装では、この関係を線形化して直接使うのではなく、相対風速度から \(\alpha_a\) と \(\beta_a\) を計算する。
時刻歴では、横滑り角を次の列に分けて保存する。
| 列 | 意味 |
|---|---|
beta_body | 機体の剛体速度 \(u,v,w\) から求めた横滑り角。線形 4 状態モデルの状態量 \(\beta\) に対応する |
beta_g | 突風による外乱横滑り角 |
beta | 相対風 \(u_a,v_a,w_a\) から求めた有効対気横滑り角 \(\beta_a\) |
したがって、出力列 beta を線形モデルの状態量 \(\beta\) と読み替えてはならない。
小擾乱の範囲では、
\begin{align} \mathrm{beta} \simeq \mathrm{beta\_body} +\mathrm{beta\_g} \end{align}
である。
空力係数の計算
空力係数は、.stab の基準値と安定微係数を使って一次近似する。
基準迎角と基準横滑り角を、
\begin{align} \alpha_0, \beta_0 \end{align}
とする。
相対風から得た迎角・横滑り角を \(\alpha_a\)、\(\beta_a\) とすれば、摂動は、
\begin{align} \Delta\alpha &= \alpha_a-\alpha_0, \\ \Delta\beta &= \beta_a-\beta_0 \end{align}
である。
無次元角速度は、
\begin{align} \hat{p} &= \frac{p b_{\mathrm{ref}}}{2V_a}, \\ \hat{q} &= \frac{q c_{\mathrm{ref}}}{2V_a}, \\ \hat{r} &= \frac{r b_{\mathrm{ref}}}{2V_a} \end{align}
である。
小擾乱モデルを機体軸横力係数で表せば、横力係数は概念的には次の一次式になる。
\begin{align} C_Y ={}& C_{Y0}+C_{Y\alpha}\Delta\alpha+C_{Y\beta}\Delta\beta+C_{Y\hat{p}}\hat{p}+C_{Y\hat{q}}\hat{q}+C_{Y\hat{r}}\hat{r} \\ &+C_{Y\delta_a}\delta_a+C_{Y\delta_e}\delta_e+C_{Y\delta_r}\delta_r \end{align}
ただし、これは実装の入力データを直接示す式ではない。
実装では .stab の \(C_L\)、\(C_D\)、\(C_S\) をそれぞれ基準値と微係数から一次評価し、その後で \(C_X\)、\(C_Y\)、\(C_Z\) へ変換する。
ローリングモーメント係数は、
\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\) を .stab の基準値と微係数で計算する。
VSPAERO の .stab では、力係数として \(C_L\)、\(C_D\)、\(C_S\) が出力されるため、スクリプトではこれらを通常の機体軸、すなわち \(x_b\) 前方、\(y_b\) 右方、\(z_b\) 下方の係数へ戻して使う。
完全に書くと、
\begin{align} C_X &= -\left(C_D\cos\beta_a+C_S\sin\beta_a\right)\cos\alpha_a+C_L\sin\alpha_a, \\ C_Y &= -C_D\sin\beta_a+C_S\cos\beta_a, \\ C_Z &= -\left(C_D\cos\beta_a+C_S\sin\beta_a\right)\sin\alpha_a-C_L\cos\alpha_a \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}
である。
ここで、動圧は、
\begin{align} q_\infty = \frac{1}{2}\rho V_a^2 \end{align}
である。
6 自由度剛体運動方程式
6 自由度剛体運動では、速度、角速度、位置、姿勢角を時間積分する。
状態量は、概念的には次である。
\begin{align} \mathbf{x}=\begin{bmatrix}u & v & w & p & q & r & x_e & y_e & z_e & \phi & \theta & \psi\end{bmatrix}^{\mathrm{T}} \end{align}
速度の式は、
\begin{align} \dot{u} &= rv-qw-g\sin\theta+\frac{X+T}{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} \dot{p} &= \frac{I_{zz}\left[L+\left(I_{yy}-I_{zz}\right)qr+I_{xz}pq\right] +I_{xz}\left[N+\left(I_{xx}-I_{yy}\right)pq-I_{xz}qr\right]} {I_{xx}I_{zz}-I_{xz}^2}, \\ \dot{q} &= \frac{M+\left(I_{zz}-I_{xx}\right)pr-I_{xz}\left(p^2-r^2\right)}{I_{yy}}, \\ \dot{r} &= \frac{I_{xz}\left[L+\left(I_{yy}-I_{zz}\right)qr+I_{xz}pq\right] +I_{xx}\left[N+\left(I_{xx}-I_{yy}\right)pq-I_{xz}qr\right]} {I_{xx}I_{zz}-I_{xz}^2} \end{align}
である。
位置の運動学は、NED 座標系を使うと、
\begin{align} \dot{x}_e &= u\cos\theta\cos\psi +v\left(\sin\phi\sin\theta\cos\psi-\cos\phi\sin\psi\right) +w\left(\cos\phi\sin\theta\cos\psi+\sin\phi\sin\psi\right), \\ \dot{y}_e &= u\cos\theta\sin\psi +v\left(\sin\phi\sin\theta\sin\psi+\cos\phi\cos\psi\right) +w\left(\cos\phi\sin\theta\sin\psi-\sin\phi\cos\psi\right), \\ \dot{z}_e &= -u\sin\theta+v\sin\phi\cos\theta+w\cos\phi\cos\theta \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_crosswind_gust_from_stab() は、この右辺を作り、solve_ivp() で時間積分する。
応答指標の計算
横風突風応答では、時刻歴そのものに加えて、ピーク値を見ると設計比較がしやすい。
calculate_crosswind_gust_response_indices() は、主に次の量を計算する。
- 横滑り角 \(\beta\) のピーク値
- ロール率 \(p\) のピーク値
- ヨーレート \(r\) のピーク値
- 無次元ロール率 \(\hat{p}\) のピーク値
- 無次元ヨーレート \(\hat{r}\) のピーク値
- バンク角 \(\phi\) のピーク変化
- それぞれのピーク到達時刻
ピーク値は、符号付きの値と絶対値を分けて見る。
たとえば、ある量 \(x(t)\) のピーク絶対値は、
\begin{align} \lvert x\rvert_{\max}=\max_t\lvert x(t)\rvert \end{align}
である。
また、初期バンク角を \(\phi_0\) とすると、バンク角変化は、
\begin{align} \Delta\phi(t)=\phi(t)-\phi_0 \end{align}
である。
横風突風に対するロール応答を見るときは、\(\phi(t)\) そのものだけでなく \(\Delta\phi(t)\) も見る。
使い方
ディレクトリ構成
最小構成は次のようになる。
project/
├── RollRudderGain.py
├── TrimTurnSolver.py
├── SampleGlider.stab
└── run_crosswind_gust_6dof.py.stab ファイルは、OpenVSP / VSPAERO の安定微係数解析で作成しておく。
入力する値
simulate_6dof_crosswind_gust_from_stab() に与える主な入力は次である。
| 引数 | 意味 |
|---|---|
stab_path | VSPAERO の .stab ファイル |
mass | 機体質量 |
Ixx | ロール軸まわりの慣性モーメント |
Iyy | ピッチ軸まわりの慣性モーメント |
Izz | ヨー軸まわりの慣性モーメント |
Ixz | 慣性乗積 |
Uds | 右から吹く方向を正とする符号付きの有次元ピーク突風速度。 Part 25 の EAS 値を使う場合は解析高度の TAS 相当値へ変換する |
H | 突風勾配距離 |
t_final | 積分終了時刻 |
gust_start_time | 突風開始時刻 |
rho | 密度。省略時は .stab の Rho を使う |
control_map | Control Group と舵角の対応 |
角度と舵角は rad で指定する。
import math
phi0 = math.radians(0.0)
theta0 = None
psi0 = math.radians(0.0)長さ、速度、密度、質量、慣性モーメントは一貫した単位系で指定する。
実行例：横風突風 6 自由度応答
example_crosswind_gust_response.ipynb と同じ SampleGlider.stab、質量、慣性モーメント、突風条件を用いる例は次の通りである。
from RollRudderGain import simulate_6dof_crosswind_gust_from_stab
history = simulate_6dof_crosswind_gust_from_stab(
stab_path="SampleGlider.stab",
mass=100.0,
Ixx=1000.0,
Iyy=75.0,
Izz=1000.0,
Ixz=0.0,
Uds=3.0,
H=9.144,
t_final=3.0,
rho=1.225,
gust_start_time=0.0,
delta_a=0.0,
delta_e=None,
delta_r=0.0,
trim_elevator=True,
trim_thrust=False,
phi0=0.0,
theta0=None,
psi0=0.0,
max_step=0.02,
)
print(history.columns)
print(history.tail())ここでは、右から吹く突風として \(U_{ds}=3.0\ \mathrm{m/s}\)、\(H=9.144\ \mathrm{m}\) としている。
\(H=9.144\ \mathrm{m}\) は \(30\ \mathrm{ft}\) に対応する。
速度は .stab の基準飛行速度を使うため、突風通過時間は内部で、
\begin{align} t_e=\frac{2H}{V} \end{align}
として決まる。
trim_elevator=True により、エレベータ舵角は突風入力前の初期ピッチングモーメントが小さくなるように設定される。
trim_thrust=False なので推力はゼロである。
この初期化は完全な定常滑空トリムを保証しないため、横風突風による横・方向応答と同時に、速度、迎角、ピッチ姿勢も変化し得る。
実行例：時刻歴を CSV へ保存する
6 自由度応答を CSV に保存する場合は、write_6dof_history_csv() を使う。
from RollRudderGain import write_6dof_history_csv
write_6dof_history_csv(
history,
"crosswind_gust_6dof_history.csv",
)これにより、時刻、速度成分、角速度、姿勢角、迎角、横滑り角、突風速度、舵角などを後で確認できる。
実行例：時刻歴を描く
時刻歴を確認する場合は、plot_6dof_history() を使う。
from RollRudderGain import plot_6dof_history
fig = plot_6dof_history(
history,
plot_path="crosswind_gust_6dof_history.png",
show=False,
degrees=True,
)横風突風を受けたときに、横滑り角、ヨーレート、ロール率、バンク角がどのように立ち上がるかを見る。
実行例：応答指標を計算する
横風突風応答のピーク値を計算する場合は、calculate_crosswind_gust_response_indices() を使う。
from RollRudderGain import calculate_crosswind_gust_response_indices
metrics = calculate_crosswind_gust_response_indices(
history,
)
print(metrics)Vref と Bref を明示したい場合は、次のように与える。
metrics = calculate_crosswind_gust_response_indices(
history,
Vref=30.0,
Bref=17.5,
)これにより、無次元ロール率や無次元ヨーレートを評価する基準を明示できる。
線形モデルではピーク外乱横滑り角を \(U_{ds}/V\) と近似するが、この後処理では、
\begin{align} \beta_{g,\max}^{\mathrm{6DOF}} =\tan^{-1}\left(\frac{U_{ds}}{V_{\mathrm{ref}}}\right) \end{align}
を正規化に用いる。
crosswind_gust_beta_peak_input には、この幾何学的な値が入る。
Control Group の対応を明示する場合
.stab の Control Group 名が標準的でない場合は、control_map を明示する。
control_map = {
"delta_a": "ConGrp_1",
"delta_e": "ConGrp_2",
"delta_r": "ConGrp_3",
}Control Group 名で指定することもできる。
control_map = {
"delta_a": "AILERON_GROUP",
"delta_e": "ELEVATOR_GROUP",
"delta_r": "RUDDER_GROUP",
}横風突風応答だけを見る場合でも、初期エレベータトリムや舵角摂動の評価に Control Group の対応を使うため、最初に確認しておくとよい。
Vv-Gamma チャート後処理から使う場合
Vv-Gamma チャートの各ケースについて横風突風 6 自由度応答を計算する場合は、postprocess_vv_gamma_cases() から呼び出す。
from VvGammaChart import postprocess_vv_gamma_cases
results = postprocess_vv_gamma_cases(
"vv_wtip_stability_sweep.csv",
mass=100.0,
inertia={
"Ixx": 1000.0,
"Iyy": 75.0,
"Izz": 1000.0,
"Ixz": 0.0,
},
rho=1.225,
run_crosswind_gust_6dof=True,
crosswind_gust_Uds=3.0,
crosswind_gust_H=9.144,
crosswind_gust_t_final=3.0,
write_crosswind_gust_history=True,
plot_crosswind_gust_history=True,
history_output_dir="6dof_crosswind_history",
output_csv_path="vv_gamma_metrics.csv",
)この使い方では、各 .stab ケースに対して横風突風応答を計算し、ピーク値を vv_gamma_metrics.csv にまとめる。
これにより、垂直尾翼容積比 \(V_v\) と等価上反角 \(\bar{\Gamma}_{\mathrm{eff}}\) の平面上で、横風突風応答の大きさを比較できる。
出力データ
6 自由度時刻歴
simulate_6dof_crosswind_gust_from_stab() の返り値は pandas.DataFrame である。
主に見る列は次である。
- time
- u, v, w
- p, q, r
- x_e, y_e, z_e
- phi, theta, psi
alpha,beta：相対風から求めた有効迎角・有効対気横滑り角alpha_body,beta_body：機体の剛体速度から求めた迎角・横滑り角phat,qhat,rhat- delta_a, delta_e, delta_r
- gust_velocity_y
- airmass_velocity_y
- beta_g
- 空力係数や力・モーメントに関する列
列名は実装や保存設定によって増える場合がある。
まずは history.columns を確認するとよい。
gust_velocity_y と beta_g は右から吹く方向を正とする。
beta_body は機体運動側の横滑り角、beta は突風を含む相対風から求めた有効対気横滑り角である。
airmass_velocity_y は機体軸 \(+y_b\) 方向を正とする実際の大気速度成分であるため、
\begin{align} \mathrm{airmass\_velocity\_y} =-\mathrm{gust\_velocity\_y} \end{align}
の関係になる。
応答指標
calculate_crosswind_gust_response_indices() の返り値は辞書である。
主な要素は次である。
| キー | 内容 |
|---|---|
crosswind_gust_peak_beta | 有効対気横滑り角 beta、すなわち \(\beta_a\) の符号付きピーク値 |
crosswind_gust_max_abs_beta | 有効対気横滑り角 beta、すなわち \(\beta_a\) の最大絶対値 |
crosswind_gust_peak_p | ロール率の符号付きピーク値 |
crosswind_gust_max_abs_p | ロール率の最大絶対値 |
crosswind_gust_peak_r | ヨーレートの符号付きピーク値 |
crosswind_gust_max_abs_r | ヨーレートの最大絶対値 |
crosswind_gust_max_phi_delta | バンク角変化の符号付き最大応答 |
crosswind_gust_max_abs_phi_delta | バンク角変化の最大絶対値 |
crosswind_gust_max_abs_phi_delta_time | 最大バンク角変化が生じた時刻 |
crosswind_gust_bank_index_abs_per_beta_g | 最大バンク角変化をピーク外乱横滑り角で正規化した指標 |
crosswind_gust_phat_index_abs_per_beta_g | 最大無次元ロール率をピーク外乱横滑り角で正規化した指標 |
crosswind_gust_rhat_index_abs_per_beta_g | 最大無次元ヨーレートをピーク外乱横滑り角で正規化した指標 |
ピーク値だけでなく、ピーク到達時刻も一緒に見ると、応答が突風突入直後に出ているのか、突風通過後に遅れて出ているのかを確認できる。
おわりに
この記事では、横風突風応答に対する 6 自由度剛体運動をシミュレーションする Python スクリプトについて説明した。
コメント