PR

航空機の定常旋回パラメータを計算する Python スクリプト

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 式から求める。

高度維持条件は課さず、得られた姿勢と速度から降下率を後処理で計算する。

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

  1. .stab ファイルを読み込む。
  2. 参照面積、参照翼弦、参照翼幅、基準飛行条件を取得する。
  3. 安定微係数表を取得する。
  4. Control Surface Group の列を \(\delta_a\)、\(\delta_e\)、\(\delta_r\) に対応づける。
  5. 固定する 3 つのパラメータを受け取る。
  6. 未知量の初期値を作る。
  7. 高度維持条件、力のつり合い、モーメントのつり合いの残差を作る。
  8. scipy.optimize.least_squares() で残差を小さくする。
  9. 解、残差、派生量、収束情報を辞書として返す。

入力データ

.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_1ConGrp_2ConGrp_3 が、それぞれエルロン、エレベータ、ラダーのどれに対応するかを決める必要がある。

TrimTurnSolver.py では、resolve_control_columns_from_stab() がこの対応付けを行う。

対応の考え方は以下の通り。

  • control_map を明示した場合は、その対応を使う。
  • control_map がない場合は、Control Group 名に AILERON、ELEVATOR、RUDDER が含まれるかを見て対応を決める。
  • それでも決まらない場合は、G103A 形式の慣例として ConGrp_1ConGrp_2ConGrp_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_pathVSPAERO の .stab ファイル
mass機体質量
rho密度。省略時は .stab の Rho を使う
initial_guess未知量の初期値
bounds未知量の下限・上限
control_mapControl 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",
}

出力データ

返り値は辞書である。

主な要素は以下の通り。

キー内容
passedleast_squares() が正常終了し、最大絶対残差が許容値以下なら True
messageleast_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_residualresiduals の最大絶対値
residual_tolpassed 判定に用いる無次元残差許容値
costleast_squares() のコスト
optimalityleast_squares() の一次最適性指標
nfev残差関数の評価回数
rho実際に使用した密度
rho_source密度が引数由来か .stab 由来か
control_columns\(\delta_a\)、\(\delta_e\)、\(\delta_r\) と .stab 列の対応
stab使用した .stab の参照値、基準条件、Control Group 情報

まず見るべき値は passedsolutionresiduals である。

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 の単位系と massrhog の単位系が一致しているか。
  • Control Group の対応が正しいか。

おわりに

VSPAERO の .stab ファイルを入力として、航空機の定常旋回パラメータを計算する Python スクリプトについて説明した。

コメント