VSPAERO の .stab ファイルを入力として、航空機の定常旋回パラメータを計算する Python スクリプトについて説明する。
はじめに
この記事では、VSPAERO の .stab ファイルを入力として、航空機の定常旋回パラメータを計算する Python スクリプトについて説明する。
定常旋回では、速度、迎角、横滑り角、バンク角、ピッチ角、旋回率、舵角、推力が互いに連成する。
たとえば、速度と旋回率を指定しても、それだけで旋回状態が決まるわけではない。
高度を維持するのか、横滑りをゼロにするのか、エルロンを使わないのか、滑空機として推力をゼロにするのかによって、解くべき未知量が変わる。
定式化の詳細は、前の記事「航空機の定常旋回運動の定式化」に譲る。
本記事では、その定式化を Python でどのように解くかに集中する。
それではいってみよう。
プログラムの全体像
今回の中心となるファイルは TrimTurnSolver.py である。
このファイルでは、動力機の定常・高度維持旋回と、滑空機の定常降下旋回を扱う。
動力機の定常・高度維持旋回では、主要なパラメータを
\begin{aligned} V,\alpha,\beta,\phi,\theta,\Omega,\delta_e,\delta_a,\delta_r,T \end{aligned}
とする。
この 10 個のうち 3 個を固定し、残り 7 個を、力のつり合い 3 式、モーメントのつり合い 3 式、高度維持条件 1 式から求める。
滑空機の定常降下旋回では、推力を
\begin{aligned} T=0 \end{aligned}
に固定する。
したがって、解く対象は
\begin{aligned} V,\alpha,\beta,\phi,\theta,\Omega,\delta_e,\delta_a,\delta_r \end{aligned}
の 9 個である。
このうち 3 個を固定し、残り 6 個を、力のつり合い 3 式とモーメントのつり合い 3 式から求める。
高度維持条件は課さず、得られた姿勢と速度から降下率を後処理で計算する。
処理の流れは次のようになる。
.stabファイルを読み込む。- 参照面積、参照翼弦、参照翼幅、基準飛行条件を取得する。
- 安定微係数表を取得する。
- Control Surface Group の列を \(\delta_a\)、\(\delta_e\)、\(\delta_r\) に対応づける。
- 固定する 3 つのパラメータを受け取る。
- 未知量の初期値を作る。
- 高度維持条件、力のつり合い、モーメントのつり合いの残差を作る。
scipy.optimize.least_squares()で残差を小さくする。- 解、残差、派生量、収束情報を辞書として返す。
入力データ
.stab ファイル
入力として使う .stab ファイルは、OpenVSP / VSPAERO の STABILITY_DEFAULT 解析で出力される安定微係数ファイルである。
.stab には、代表的に次の情報が含まれる。
- 参照面積 \(S_{\mathrm{ref}}\)
- 参照翼弦 \(c_{\mathrm{ref}}\)
- 参照翼幅 \(b_{\mathrm{ref}}\)
- 重心位置
- Mach 数
- 迎角
- 横滑り角
- 密度
- 飛行速度
- 基準ケースの空力係数
- 迎角、横滑り角、角速度、舵角に対する微係数
本スクリプトでは、.stab の CL、CD、CS、CMl、CMm、CMn を用いる。
VSPAERO の .stab では、力係数とモーメント係数の座標系に注意が必要である。
本スクリプトでは、.stab の CL、CD、CS から、一般的な機体軸
\begin{aligned} +x_b &: \text{前方}, \\ +y_b &: \text{右方}, \\ +z_b &: \text{下方} \end{aligned}
に対応する \(C_X\)、\(C_Y\)、\(C_Z\) を計算して使う。
Control Surface Group
VSPAERO の舵効き微係数は、Control Surface Group に対して出力される。
個別舵面の幾何舵角に対する微係数とは扱いが異なる。
そのため、.stab の ConGrp_1、ConGrp_2、ConGrp_3 が、それぞれエルロン、エレベータ、ラダーのどれに対応するかを決める必要がある。
TrimTurnSolver.py では、resolve_control_columns_from_stab() がこの対応付けを行う。
対応の考え方は以下の通り。
control_mapを明示した場合は、その対応を使う。control_mapがない場合は、Control Group 名に AILERON、ELEVATOR、RUDDER が含まれるかを見て対応を決める。- それでも決まらない場合は、G103A 形式の慣例として
ConGrp_1、ConGrp_2、ConGrp_3を使う。
プログラムの解説
ファイル構成
TrimTurnSolver.py の主な関数は以下の通り。
| 関数 | 役割 |
|---|---|
read_vspaero_stab() | .stab ファイルから基準値、基準空力係数、安定微係数表を読む |
resolve_control_columns_from_stab() | .stab の Control Group 列を \(\delta_a\)、\(\delta_e\)、\(\delta_r\) に対応づける |
_steady_turn_derived() | 速度成分、無次元角速度、高度方向速度などの派生量を計算する |
_linear_aero_coefficients() | 線形化した空力係数とモーメント係数を計算する |
_steady_turn_residuals() | 力のつり合い、モーメントのつり合い、高度維持条件の残差を作る |
solve_steady_level_turn() | 3 個の量を固定して定常・高度維持旋回を解く |
solve_steady_gliding_turn() | 3 個の量を固定して無動力の定常降下旋回を解く |
solve_rudder_limit_turn() | 正負のラダー限界における定常旋回を解き、成立した解から最大の \( |
.stab ファイルの読み取り
.stab の読み取りは read_vspaero_stab() で行う。
この関数では、まず .stab のヘッダから、次の値を読む。
\begin{aligned} S_{\mathrm{ref}},c_{\mathrm{ref}},b_{\mathrm{ref}},x_{\mathrm{cg}},y_{\mathrm{cg}},z_{\mathrm{cg}} \end{aligned}
次に、基準飛行条件として、
\begin{aligned} M,\alpha_0,\beta_0,\rho,V \end{aligned}
を読む。
さらに、Base_Aero 行から基準空力係数を読み、Coef Base Aero ... で始まる表から安定微係数を読み取る。
読み取った内容は、次のような辞書と表に整理される。
stab.references
stab.base_condition
stab.base_aero
stab.derivatives
stab.control_groups
この段階では、.stab に書かれている値をそのまま読む。
単位変換は行わない。
空力係数の計算
定常旋回の残差を作るには、任意の状態量と舵角における空力係数が必要である。
スクリプトでは、.stab の基準ケースと安定微係数を使い、空力係数を一次近似で計算する。
横方向の力係数は
\begin{aligned} C_Y =C_{Y0} +C_{Y\beta}\Delta\beta +C_{Y\hat p}\hat p +C_{Y\hat r}\hat r +C_{Y\delta_a}\delta_a +C_{Y\delta_r}\delta_r \end{aligned}
である。
ピッチングモーメント係数は
\begin{aligned} C_m =C_{m0} +C_{m\alpha}\Delta\alpha +C_{m\hat q}\hat q +C_{m\delta_e}\delta_e \end{aligned}
である。
ローリングモーメント係数とヨーイングモーメント係数は
\begin{aligned} C_l &= C_{l0} +C_{l\beta}\Delta\beta +C_{l\hat p}\hat p +C_{l\hat r}\hat r +C_{l\delta_a}\delta_a +C_{l\delta_r}\delta_r, \\ C_n &= C_{n0} +C_{n\beta}\Delta\beta +C_{n\hat p}\hat p +C_{n\hat r}\hat r +C_{n\delta_a}\delta_a +C_{n\delta_r}\delta_r \end{aligned}
である。
ここで、\(\Delta\alpha\) と \(\Delta\beta\) は、.stab の基準迎角・基準横滑り角からの差である。
派生量の計算
旋回パラメータから、速度成分と角速度を計算する。
速度成分は
\begin{aligned} u &= V\cos\alpha\cos\beta, \\ v &= V\sin\beta, \\ w &= V\sin\alpha\cos\beta \end{aligned}
である。
定常旋回では、方位角速度を \(\Omega\) とすると、機体軸角速度は
\begin{aligned} p &= -\Omega\sin\theta, \\ q &= \Omega\sin\phi\cos\theta, \\ r &= \Omega\cos\phi\cos\theta \end{aligned}
である。
VSPAERO の安定微係数は無次元角速度に対する微係数なので、
\begin{aligned} \hat p &= \frac{p b_{\mathrm{ref}}}{2V}, \\ \hat q &= \frac{q c_{\mathrm{ref}}}{2V}, \\ \hat r &= \frac{r b_{\mathrm{ref}}}{2V} \end{aligned}
に変換して使う。
高度方向の速度は、NED 座標系の下向き速度として
\begin{aligned} \dot z_e =-u\sin\theta +v\sin\phi\cos\theta +w\cos\phi\cos\theta \end{aligned}
である。
高度を \(h=-z_e\) とすれば、
\begin{aligned} \dot h=-\dot z_e \end{aligned}
である。
滑空機の定常降下旋回では、この値から降下率を確認する。
水平飛行経路速度は、機首方位方向と、その右向き方向の水平速度成分を
\begin{aligned} V_{H\parallel} &= u\cos\theta +v\sin\phi\sin\theta +w\cos\phi\sin\theta, \\ V_{H\perp} &= v\cos\phi-w\sin\phi \end{aligned}
として、
\begin{aligned} V_H =\sqrt{V_{H\parallel}^2+V_{H\perp}^2} \end{aligned}
で求める。
定常旋回では、水平飛行経路方位角も \(\Omega\) で回転する。
したがって、実際の水平飛行経路半径は、
\begin{aligned} R_H= \begin{cases} \dfrac{V_H}{|\Omega|}, & \Omega\neq0, \\ \infty, & \Omega=0 \end{cases} \end{aligned}
である。
残差の作り方
_steady_turn_residuals() では、力とモーメントの残差を作る。
推力が機体軸 \(x_b\) 正方向に働くとして、空気力を
\begin{aligned} X &= T+q_\infty S C_X, \\ Y &= q_\infty S C_Y, \\ Z &= q_\infty S C_Z \end{aligned}
とする。
動圧は
\begin{aligned} q_\infty=\frac{1}{2}\rho V^2 \end{aligned}
である。
定常旋回では、機体軸方向の力のつり合いは
\begin{aligned} X &= m(qw-rv+g\sin\theta), \\ Y &= m(ru-pw-g\sin\phi\cos\theta), \\ Z &= m(pv-qu-g\cos\phi\cos\theta) \end{aligned}
である。
モーメントのつり合いは、角速度 2 次の慣性連成を無視した範囲では
\begin{aligned} C_l &= 0, \\ C_m &= 0, \\ C_n &= 0 \end{aligned}
である。
定常・高度維持旋回では、さらに高度維持条件
\begin{aligned} -u\sin\theta+v\sin\phi\cos\theta+w\cos\phi\cos\theta=0 \end{aligned}
を加える。
したがって、動力機の定常・高度維持旋回では、残差は 7 成分である。
実装では、力残差を基準空気力 \(q_\infty S\) で、鉛直速度残差を飛行速度 \(V\) で正規化する。
\begin{aligned} \mathbf{r} =\begin{bmatrix} \dfrac{X-X_{\mathrm{req}}}{q_\infty S} \\ \dfrac{Y-Y_{\mathrm{req}}}{q_\infty S} \\ \dfrac{Z-Z_{\mathrm{req}}}{q_\infty S} \\ C_l \\ C_m \\ C_n \\ \dfrac{\dot z_e}{V} \end{bmatrix} \end{aligned}
ここで、
\begin{aligned} X_{\mathrm{req}}&=m(qw-rv+g\sin\theta), \\ Y_{\mathrm{req}}&=m(ru-pw-g\sin\phi\cos\theta), \\ Z_{\mathrm{req}}&=m(pv-qu-g\cos\phi\cos\theta) \end{aligned}
である。
無動力の定常降下旋回では、\(T=0\) を固定し、高度維持条件を残差に含めない。
したがって、残差は 6 成分である。
\begin{aligned}
\mathbf{r}_{\mathrm{glide}}
=\begin{bmatrix}
\dfrac{X-X_{\mathrm{req}}}{q_\infty S} \\
\dfrac{Y-Y_{\mathrm{req}}}{q_\infty S} \\
\dfrac{Z-Z_{\mathrm{req}}}{q_\infty S} \\ C_l \\ C_m \\ C_n \end{bmatrix}
\end{aligned}
したがって、max_abs_residual と residual_tol は、N や m/s ではなく、上記の無次元残差に対する値である。
solve_steady_level_turn()
solve_steady_level_turn() は、動力機の定常・高度維持旋回を解く関数である。
呼び出し時には、10 個の旋回パラメータのうち、ちょうど 3 個を fixed に指定する。
result = solve_steady_level_turn(
fixed={
"V": 30.0,
"Omega": 0.08,
"beta": 0.0,
},
stab_path="G103A.stab",
mass=580.0,
rho=1.225,
)
この例では、速度 \(V\)、旋回率 \(\Omega\)、横滑り角 \(\beta\) を固定している。
これは、指定速度・指定旋回率で横滑りなしの旋回を求める使い方である。
ラダーのみ旋回を解きたい場合は、エルロン舵角を固定する。
result = solve_steady_level_turn(
fixed={
"V": 30.0,
"Omega": 0.08,
"delta_a": 0.0,
},
stab_path="G103A.stab",
mass=580.0,
rho=1.225,
)
この場合は、速度 \(V\)、旋回率 \(\Omega\)、エルロン舵角 \(\delta_a\) を固定し、横滑り角 \(\beta\) やラダー舵角 \(\delta_r\) は結果として求まる。
solve_steady_gliding_turn()
solve_steady_gliding_turn() は、無動力の定常降下旋回を解く関数である。
滑空機として扱うため、推力は内部で
\begin{aligned} T=0 \end{aligned}
に固定される。
したがって、fixed に \(T\) を指定しない。
result = solve_steady_gliding_turn(
fixed={
"V": 30.0,
"Omega": 0.08,
"beta": 0.0,
},
stab_path="G103A.stab",
mass=580.0,
rho=1.225,
)
この例では、速度 \(V\)、旋回率 \(\Omega\)、横滑り角 \(\beta\) を固定し、無動力で釣り合う姿勢と舵角を求める。
滑空機では高度維持条件を解かない。
そのため、結果には derived["z_dot"]、derived["h_dot"]、derived["sink_rate"] が含まれる。
solve_rudder_limit_turn()
solve_rudder_limit_turn() は、正負の最大ラダー舵角における定常旋回トリムをそれぞれ解き、成立した解のうちバンク角の絶対値が大きい方を選ぶ関数である。
評価する 2 ケースは、
\begin{aligned} V &= V_0, & \delta_a &= \text{一定}, & \delta_r &= -\delta_{r,\max}, \\ V &= V_0, & \delta_a &= \text{一定}, & \delta_r &= +\delta_{r,\max} \end{aligned}
である。
mode="level" では solve_steady_level_turn()、mode="gliding" では solve_steady_gliding_turn() を内部で使う。
正負の両方でトリムが成立した場合は、
\begin{aligned}
|\phi|_{\max,\mathrm{trim}}
=\max\left( \left|\phi\left(-\delta_{r,\max}\right)\right|, \left|\phi\left(+\delta_{r,\max}\right)\right| \right)
\end{aligned}
を返す。
片側だけが成立した場合は、その解を返し、rudder_limit_complete=False とする。
この関数が直接評価するのは正負のラダー限界である。
中間舵角は掃引しないため、舵角範囲全体における単調性と大域的な最大値を確認する場合は、複数の \(\delta_r\) を指定して追加のトリム計算を行う。
使い方
ディレクトリ構成
最小構成は次のようになる。
project/
├── TrimTurnSolver.py
├── AnalysisVSPAERO.py
├── G103A.stab
└── run_trim_turn.py
.stab ファイルは、OpenVSP / VSPAERO の STABILITY_DEFAULT 解析で作成しておく。
AnalysisVSPAERO.py を使って .stab を作る場合は、先にモデルの事前チェックと VSPAERO の安定微係数計算を実行しておく。
入力する値
solve_steady_level_turn() と solve_steady_gliding_turn() で共通して必要な主な入力は以下の通り。
| 引数 | 意味 |
|---|---|
fixed | 固定する 3 個の旋回パラメータ |
stab_path | VSPAERO の .stab ファイル |
mass | 機体質量 |
rho | 密度。省略時は .stab の Rho を使う |
initial_guess | 未知量の初期値 |
bounds | 未知量の下限・上限 |
control_map | Control Group と舵角の対応 |
residual_tol | 収束判定に使う残差許容値 |
fixed に指定できる値は、定常・高度維持旋回では次の 10 個である。
\begin{aligned} V,\alpha,\beta,\phi,\theta,\Omega,\delta_e,\delta_a,\delta_r,T \end{aligned}
定常降下旋回では \(T=0\) が内部で固定されるため、fixed に指定できる値は次の 9 個である。
\begin{aligned} V,\alpha,\beta,\phi,\theta,\Omega,\delta_e,\delta_a,\delta_r \end{aligned}
単位系
角度は rad で指定する。
import math
beta = math.radians(0.0)
phi = math.radians(20.0)
delta_a = math.radians(0.0)
速度、密度、質量、面積、重力加速度は、一貫した単位系で指定する。
たとえば、長さを m、時間を s、質量を kg として扱うなら、
\begin{aligned} V &: \mathrm{m/s}, \\ \rho &: \mathrm{kg/m^3}, \\ m &: \mathrm{kg}, \\ S &: \mathrm{m^2}, \\ g &: \mathrm{m/s^2} \end{aligned}
でそろえる。
rho=None とした場合は、.stab の Rho がそのまま使われる。
G103A の以下の例では、質量を 580.0 kg、密度を 1.225 kg/m^3 として明示し、.stab に保存された密度値を暗黙には使わない。
実行例:協調旋回
協調旋回を計算する場合は、横滑り角をゼロに固定する。
import math
from TrimTurnSolver import solve_steady_level_turn
result = solve_steady_level_turn(
fixed={
"V": 30.0,
"Omega": 0.06,
"beta": 0.0,
},
stab_path="G103A.stab",
mass=580.0,
rho=1.225,
initial_guess={
"alpha": math.radians(3.0),
"phi": math.radians(10.0),
"theta": math.radians(3.0),
},
)
print(result["passed"])
print(result["solution"])
print(result["derived"])
結果の solution には、固定値と求解された値を合わせた旋回パラメータが入る。
実行例:ラダーのみ旋回
ラダーのみ旋回を計算する場合は、エルロン舵角をゼロに固定する。
import math
from TrimTurnSolver import solve_steady_level_turn
result = solve_steady_level_turn(
fixed={
"V": 30.0,
"Omega": 0.06,
"delta_a": 0.0,
},
stab_path="G103A.stab",
mass=580.0,
rho=1.225,
initial_guess={
"alpha": math.radians(3.0),
"beta": math.radians(3.0),
"phi": math.radians(10.0),
"theta": math.radians(3.0),
"delta_r": math.radians(5.0),
},
)
print(result["passed"])
print(result["solution"])
print(result["residuals"])
この例では、\(\beta\) や \(\delta_r\) は、力とモーメントのつり合いから決まる結果量である。
実行例:滑空機の定常降下旋回
滑空機として、推力ゼロの定常降下旋回を計算する場合は、solve_steady_gliding_turn() を使う。
import math
from TrimTurnSolver import solve_steady_gliding_turn
result = solve_steady_gliding_turn(
fixed={
"V": 30.0,
"Omega": 0.06,
"beta": 0.0,
},
stab_path="G103A.stab",
mass=580.0,
rho=1.225,
initial_guess={
"alpha": math.radians(3.0),
"phi": math.radians(10.0),
"theta": math.radians(-3.0),
},
)
print(result["passed"])
print(result["solution"])
print(result["derived"]["sink_rate"])
この関数では高度維持条件を課さないため、得られた sink_rate を見て、滑空状態として妥当かを確認する。
実行例:ラダー限界における最大バンク角
正負の最大ラダー舵角における無動力の定常降下旋回を解き、成立した解から最大の \(|\phi|\) を選ぶ場合は、solve_rudder_limit_turn() を使う。
import math
from TrimTurnSolver import solve_rudder_limit_turn
result = solve_rudder_limit_turn(
stab_path="G103A.stab",
mass=580.0,
rho=1.225,
delta_r_max=math.radians(10.0),
mode="gliding",
V=30.0,
delta_a=0.0,
initial_guess={
"alpha": math.radians(3.0),
"beta": math.radians(3.0),
"phi": math.radians(5.0),
"theta": math.radians(-3.0),
"Omega": math.radians(3.0),
"delta_e": 0.0,
},
)
print(result["passed"])
print(result["rudder_limit_complete"])
print(result["selected_side"])
print(math.degrees(result["max_abs_phi"]))
print(result["negative_trim"])
print(result["positive_trim"])
passed=True は、正負のラダー限界の少なくとも一方でトリム解が成立したことを表す。
rudder_limit_complete=True は、正負の両方でトリム解が成立したことを表す。
selected_side には、採用した側として negative または positive が入る。
Control Group の対応を明示する場合
.stab の Control Group 名が標準的な名前でない場合は、control_map を使って対応を明示する。
result = solve_steady_level_turn(
fixed={
"V": 30.0,
"Omega": 0.06,
"delta_a": 0.0,
},
stab_path="G103A.stab",
mass=580.0,
rho=1.225,
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",
}
出力データ
返り値は辞書である。
主な要素は以下の通り。
| キー | 内容 |
|---|---|
passed | least_squares() が正常終了し、最大絶対残差が許容値以下なら True |
message | least_squares() の終了メッセージ |
fixed | 呼び出し側で固定した 3 個のパラメータ |
unknown_names | 数値的に解いたパラメータ名 |
solution | 固定値と求解値を合わせた旋回パラメータ |
derived | 速度成分、角速度、無次元角速度、水平飛行経路半径、降下率など |
coefficients | 解における \(C_L,C_D,C_S,C_l,C_m,C_n,C_X,C_Y,C_Z\) |
residuals | 正規化された力残差とモーメント残差。高度維持旋回では高度残差も含む |
max_abs_residual | residuals の最大絶対値 |
residual_tol | passed 判定に用いる無次元残差許容値 |
cost | least_squares() のコスト |
optimality | least_squares() の一次最適性指標 |
nfev | 残差関数の評価回数 |
rho | 実際に使用した密度 |
rho_source | 密度が引数由来か .stab 由来か |
control_columns | \(\delta_a\)、\(\delta_e\)、\(\delta_r\) と .stab 列の対応 |
stab | 使用した .stab の参照値、基準条件、Control Group 情報 |
まず見るべき値は passed、solution、residuals である。
solve_rudder_limit_turn() では、上の要素に加えて次の値を確認する。
| キー | 内容 |
|---|---|
rudder_limit_complete | 正負のラダー限界の両方でトリムが成立した場合は True |
delta_r_max | 入力した最大ラダー舵角 |
limiting_delta_r | 選択されたトリム解のラダー舵角 |
max_abs_phi | 選択された最大バンク角の絶対値 |
negative_trim | \(\delta_r=-\delta_{r,\max}\) における完全なトリム結果 |
positive_trim | \(\delta_r=+\delta_{r,\max}\) における完全なトリム結果 |
passed が False の場合は、次の点を確認する。
- 固定した 3 個のパラメータが物理的に無理な組合せになっていないか。
- 初期値が解から離れすぎていないか。
boundsが狭すぎないか。.stabの単位系とmass、rho、gの単位系が一致しているか。- Control Group の対応が正しいか。
おわりに
VSPAERO の .stab ファイルを入力として、航空機の定常旋回パラメータを計算する Python スクリプトについて説明した。


コメント