OpenVSP の Python API で安定微係数を計算する

OpenVSP の Python API から VSPAERO の安定微係数解析を実行し、基準空力係数、静安定微係数、動安定微係数、Control Surface Group 微係数を取得する。

GitHub - mtkbirdman/OpenVSP

Contribute to mtkbirdman/OpenVSP development by creating an account on GitHub.

github.com

はじめに

VSPAERO には、基準状態へ微小な迎角、横滑り角、角速度、Mach 数、Control Surface Group 舵角を加え、安定微係数を計算する機能がある。

本記事では、Python API の VSPAEROSweep から UnsteadyType = STABILITY_DEFAULT を指定し、基準状態まわりの安定微係数を計算する。

解析には複座グライダー G103A の OpenVSP モデルを使用する。

GROB G 103 Twin II の OpenVSPモデル

実機グライダーG103A（GROB G 103 Twin II）のOpenVSPモデルおよび数値例を示す。

mtkbirdman.com

2026.07.16

本記事で扱う主なファイルは次の通り。

G103A.vsp3
AnalysisVSPAERO.py
G103A.vspaero
G103A.stab
G103A.flt
run_wake_iter_sensitivity.ipynb
plot_wake_iter_sensitivity.ipynb
wake_iter_sensitivity.csv

↓OpenVSP の Python API のインストール方法はこちら

【Windows】OpenVSP Python APIのインストール

Windows11にOpenVSPのPython APIをインストールする

mtkbirdman.com

2024.11.30

↓OpenVSP の Python API でポーラーカーブを計算する方法はこちら

OpenVSPのPythonAPIでポーラーカーブを計算する

OpenVSPのPythonAPIで作成済みのモデルのポーラーカーブの計算（αスイープの計算）を実行する

mtkbirdman.com

2024.12.22

↓OpenVSP の Python API で舵角を設定する方法はこちら

OpenVSPのPythonAPIで舵角を設定して計算する

OpenVSPのPythonAPIで作成済みのモデルの舵角を設定してαスイープの計算を行う

mtkbirdman.com

2024.12.26

それではいってみよう。

OpenVSP/VSPAERO で計算する安定微係数

安定微係数とは

基準状態まわりの空力係数 \(C\) を線形化すると、次のように書ける。

\begin{align} C = C_0 + C_\alpha(\alpha-\alpha_0) + C_\beta(\beta-\beta_0) + C_{\hat{p}}\hat{p} + C_{\hat{q}}\hat{q} + C_{\hat{r}}\hat{r} + C_u\frac{V-V_0}{V_0} + C_{\delta_a}\delta_a + C_{\delta_e}\delta_e + C_{\delta_r}\delta_r. \end{align}

無次元角速度は次の通り。

\begin{align} \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{align}

\(C_\alpha\)、\(C_\beta\) は静安定微係数、\(C_{\hat{p}}\)、\(C_{\hat{q}}\)、\(C_{\hat{r}}\) は動安定微係数、\(C_\delta\) は操舵微係数である。

VSPAERO で利用できる Stability Type

Python API で指定できる Stability Type は次の通り。

Python enum	内容
STABILITY_OFF	安定解析を実行しない
STABILITY_DEFAULT	定常 6 自由度の有限差分安定解析
STABILITY_P_ANALYSIS	time-accurate roll analysis
STABILITY_Q_ANALYSIS	time-accurate pitch analysis
STABILITY_R_ANALYSIS	time-accurate yaw analysis
STABILITY_PITCH	simplified pitch analysis
STABILITY_ADJOINT	adjoint を用いた定常安定解析

STABILITY_DEFAULT と time-accurate P/Q/R analysis の使い分け

STABILITY_DEFAULT は、基準ケースと一定角速度を加えた摂動ケースの差から微係数を求める。

P/Q/R analysis は、機体へ微小な周期運動を与え、後半 1 周期の空力係数履歴から角速度と同相の成分を抽出する。

項目	STABILITY_DEFAULT	P/Q/R analysis
運動	一定角速度の微小摂動	微小な周期運動
数値処理	基準ケースとの前進差分	履歴の調和成分を積分
roll 出力	\(C_{\hat{p}}\)	\(C_{\hat{p}}\)
pitch 出力	\(C_{\hat{q}}\)	\(C_{\hat{q}+\dot{\alpha}}\)
yaw 出力	\(C_{\hat{r}}\)	\(C_{\hat{r}-\dot{\beta}}\)

Q analysis は \(C_{\hat{q}+\dot{\alpha}}\)、R analysis は \(C_{\hat{r}-\dot{\beta}}\) として出力される。

これらは周期運動で同時に生じる角速度と流入角変化の複合微係数である。

定常角速度に対する局所的な微係数が必要な場合は STABILITY_DEFAULT を使用する。

周期運動と時間発展する wake を含む減衰特性を評価する場合は P/Q/R analysis を使用する。

本記事の 6 自由度線形モデルでは、STABILITY_DEFAULT の \(C_{\hat{p}}\)、\(C_{\hat{q}}\)、\(C_{\hat{r}}\) を使用する。

STABILITY_DEFAULT の計算方法

基準ケースと摂動ケース

各 \((M,\alpha,\beta)\) 条件について、次のケースを計算する。

ケース	摂動量
Base Aero	基準状態
Alpha	\(\Delta\alpha = 0.01\ \mathrm{deg}\)
Beta	\(\Delta\beta = 0.01\ \mathrm{deg}\)
Roll Rate	\(\Delta p = 0.01\ \mathrm{rad/Tunit}\)
Pitch Rate	\(\Delta q = 0.01\ \mathrm{rad/Tunit}\)
Yaw Rate	\(\Delta r = 0.01\ \mathrm{rad/Tunit}\)
Mach	\(\Delta M = 0.1\)
Control Surface Group	\(\Delta\delta_g = 0.1\ \mathrm{deg}\)

微係数は基準ケースと正側摂動ケースの片側前進差分から計算する。

VSPAERO は CG からの相対位置に応じた回転流速を一様流へ加える。

\begin{align} x_r &= x-x_{\mathrm{cg}}, \\ y_r &= y-y_{\mathrm{cg}}, \\ z_r &= z-z_{\mathrm{cg}}. \end{align}

VSPAERO 内部座標における局所流速増分は、

\begin{align} \Delta u &= -z_r q-y_r r, \\ \Delta v &= -z_r p+x_r r, \\ \Delta w &= y_r p+x_r q. \end{align}

この局所流速を含む条件で、循環、wake、全機の力・モーメントを再計算する。

迎角・横滑り角微係数

迎角微係数は、degree で与えた摂動量を rad へ変換して計算する。

\begin{align} C_\alpha = \frac{C(\alpha+\Delta\alpha)-C(\alpha)}{\Delta\alpha\pi/180}. \end{align}

横滑り角微係数は、

\begin{align} C_\beta = \frac{C(\beta+\Delta\beta)-C(\beta)}{\Delta\beta\pi/180}. \end{align}

無次元角速度と動安定微係数

ロール角速度の無次元増分は、

\begin{align} \Delta\hat{p} = \frac{\Delta p\,b_{\mathrm{ref}}}{2V_\infty}. \end{align}

したがって、

\begin{align} C_{\hat{p}} = \frac{C(p=\Delta p)-C(p=0)}{\Delta p\,b_{\mathrm{ref}}/(2V_\infty)}. \end{align}

ピッチ角速度の無次元増分は、

\begin{align} \Delta\hat{q} = \frac{\Delta q\,c_{\mathrm{ref}}}{2V_\infty}. \end{align}

したがって、

\begin{align} C_{\hat{q}} = \frac{C(q=\Delta q)-C(q=0)}{\Delta q\,c_{\mathrm{ref}}/(2V_\infty)}. \end{align}

ヨー角速度の無次元増分は、

\begin{align} \Delta\hat{r} = \frac{\Delta r\,b_{\mathrm{ref}}}{2V_\infty}. \end{align}

したがって、

\begin{align} C_{\hat{r}} = \frac{C(r=\Delta r)-C(r=0)}{\Delta r\,b_{\mathrm{ref}}/(2V_\infty)}. \end{align}

Mach 微係数と U 微係数

Mach 数微係数は、

\begin{align} C_M = \frac{C(M+\Delta M)-C(M)}{\Delta M}. \end{align}

Control Surface Group 微係数

Control Surface Group \(g\) の command deflection を \(\delta_g\)、group 内の surface \(j\) に設定された Gain を \(G_{g,j}\) とする。

各 surface へ適用される舵角は、

\begin{align} \delta_{g,j} = G_{g,j}\delta_g. \end{align}

基準ケースでは、

\begin{align} \delta_{g,j}^{(0)} = G_{g,j}\delta_g^{(0)}. \end{align}

摂動ケースでは、

\begin{align} \delta_{g,j}^{(+)} = G_{g,j}\left(\delta_g^{(0)}+\Delta\delta_g\right). \end{align}

したがって、Control Surface Group 微係数は、

\begin{align} C_{\delta_g} = \frac{C\left(\left\{G_{g,j}(\delta_g^{(0)}+\Delta\delta_g)\right\}\right)-C\left(\left\{G_{g,j}\delta_g^{(0)}\right\}\right)}{\Delta\delta_g\pi/180}. \end{align}

局所線形性を仮定すると、

\begin{align} C_{\delta_g} = \sum_j G_{g,j}C_{\delta_{g,j}}. \end{align}

単位

.stab の列	単位・正規化
Alpha	per rad
Beta	per rad
p	per \(\hat{p}\)
q	per \(\hat{q}\)
r	per \(\hat{r}\)
.stab の列	単位・正規化
Mach	per Mach
U	per \(\Delta V/V_0\)
ConGrp_N	group command の per rad

derivative table は rad あたりの値である。

座標系と係数の符号

VSPAERO 内部座標と機体軸

6 自由度運動シミュレーションで一般に使用する機体軸は、\(x\) 前方、\(y\) 右方、\(z\) 下方である。

OpenVSP の内部座標は、\(X\) 後方、\(Y\) 右方、\(Z\) 上方である。

roll と yaw の符号を読むときは、この座標方向の違いを含めて確認する。

CL、CD、CS

.stab には、VSPAERO 内部軸の CFx、CFy、CFz と、aerodynamic coefficients の CL、CD、CS が出力される。

CFx、CFy、CFz は、それぞれ OpenVSP 内部の \(X\) 軸、\(Y\) 軸、\(Z\) 軸方向の力係数である。

OpenVSP 内部座標では \(X\) 軸が機体後方、\(Y\) 軸が機体右方、\(Z\) 軸が機体上方を向くため、各係数は次の向きを正とする。

係数	正方向
CFx	機体後方
CFy	機体右方
CFz	機体上方

したがって、機体軸を \(x\) 前方、\(y\) 右方、\(z\) 下方と定義すると、軸方向の力係数には次の関係がある。

\begin{align} C_X &= -C_{Fx}, \\ C_Y &= C_{Fy}, \\ C_Z &= -C_{Fz}. \end{align}

一方、CL、CD、CS は一様流に対する揚力、抗力、横力の係数である。

6 自由度運動方程式へ使用するときは、迎角と横滑り角を用いて \(C_X\)、\(C_Y\)、\(C_Z\) へ変換する。

CL は揚力係数、CD は抗力係数、CS は side-force coefficient である。

CMl、CMm、CMn

CMx、CMy、CMz は、それぞれ OpenVSP 内部の \(X\) 軸、\(Y\) 軸、\(Z\) 軸まわりのモーメント係数である。

これらは .stab に VSPAERO 内部座標の係数として出力される。

OpenVSP 内部の \(X\) 軸と \(Z\) 軸は、飛行力学で一般に用いる機体軸の \(x\) 軸と \(z\) 軸に対して反対向きである。

このため、roll、pitch、yaw moment coefficient として使用する CMl、CMm、CMn とは次の関係がある。

\begin{align} C_l &= -C_{Mx}, \\ C_m &= C_{My}, \\ C_n &= -C_{Mz}. \end{align}

.stab の CMl、CMm、CMn は、飛行力学モデルの \(C_l\)、\(C_m\)、\(C_n\) として使用する。

CL、CD、CS から機体軸力係数への変換

6 自由度運動シミュレーションなどを行う場合は、\(C_L\)、\(C_D\)、\(C_S\) を \(C_X\)、\(C_Y\)、\(C_Z\) へ変換する。

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

特に、横滑り角がゼロの場合は、

\begin{align} C_Y &= -C_D\sin 0+C_S\cos 0 \\ &= C_S. \end{align}

Thin と Thick

VSPAERO で解析対象を設定するときは、各コンポーネントを Thin または Thick として扱う。

ここでは両者の違いを整理したうえで、G103A モデルに用いた Set 構成を確認する。

Thin geometry

Thin geometry は、翼を厚みのない代表面として扱う。

主翼、水平尾翼、垂直尾翼、カナード、winglet など、循環と wake が主要な揚力面に適している。

Thick geometry

Thick geometry は、コンポーネントの外形面を厚みのある物体として扱う。

胴体、ナセル、ポッド、外部タンクなど、体積排除と body-induced flow を表現したいコンポーネントに適している。

公式の基本方針

OpenVSP 3.45 以降では、コンポーネントごとに Thin または Thick を指定する。

基本構成は次の通り。

• 揚力面：Thin
• 非揚力 body：Thick

Thin 同士、Thin と Thick、Thick 同士の intersection と trimming は OpenVSP 側で処理する。

本プログラムでの Set 構成

この方針を適用した場合の G103A モデルの Set 構成を確認する。

ThickGeom
└─ FuselageGeom

ThinGeom
├─ WingGeom
├─ HTailGeom
└─ VTailGeom

All や Shown へ依存すると、モデル編集や表示状態の変更により解析対象が変化する可能性があるため、G103A では解析専用の named Set を使用する。

Python API の解析入力

VSPAEROSweep には、解析対象、飛行条件、reference quantities、CG、solver、wake、Stability Type に関する入力がある。

vsp_stability_derivatives() で引数として明示した値は Python 側から設定し、指定しなかった optional input は .vsp3 の VSPAEROSettings を維持する。

解析対象と Set

主な入力は次の通り。

入力	内容
GeomSet	Thick geometry の Set
ThinGeomSet	Thin geometry の Set
UseModeFlag	Set の代わりに Mode を使用するフラグ
入力	内容
ModeID	使用する Mode ID

飛行条件

入力	内容
AlphaStart, AlphaEnd, AlphaNpts	迎角条件
BetaStart, BetaEnd, BetaNpts	横滑り角条件
MachStart, MachEnd, MachNpts	Mach 数条件
ReCref, ReCrefEnd, ReCrefNpts	Reynolds 数条件
Vinf	一様流速
Rho	密度

単一条件では、各 Npts を 1 とし、Start と End を同じ値へ設定する。

reference quantities

入力	内容
RefFlag	manual または component reference
WingID	reference wing の Geom ID
Sref	基準面積
bref	基準スパン
cref	基準コード長
MACFlag	MAC を使用するフラグ
ScurveFlag	curve area を使用するフラグ

動安定微係数の正規化に \(b_{\mathrm{ref}}\)、\(c_{\mathrm{ref}}\)、\(V_\infty\) を使用するため、

\begin{align} b_{\mathrm{ref}}>0,\quad c_{\mathrm{ref}}>0,\quad V_\infty>0 \end{align}

を満たす必要がある。

CG と moment reference

入力	内容
Xcg	moment reference の \(x\) 座標
Ycg	moment reference の \(y\) 座標
Zcg	moment reference の \(z\) 座標
CGGeomSet	CG 計算に使用する Set
UseCGModeFlag	CG 計算に Mode を使用するフラグ
CGModeID	CG 計算用 Mode ID

CG は moment coefficient の基準点であり、P/Q/R analysis では周期回転の中心でもある。

solver・wake

入力	内容
NCPU	CPU 数
FixedWakeFlag	fixed wake の使用
WakeNumIter	wake 反復回数
NumWakeNodes	wake node 数
ImplicitWake	implicit wake の使用
WakeRelax	wake relaxation factor
ForwardGMRESConvergenceFactor	forward GMRES の収束係数
AdjointGMRESConvergenceFactor	adjoint GMRES の収束係数
NonLinearConvergenceFactor	nonlinear solver の収束係数

Stability Type

UnsteadyType へ Stability Type を設定する。

vsp.SetIntAnalysisInput(
    analysis_name,
    "UnsteadyType",
    [vsp.STABILITY_DEFAULT],
    0,
)

プログラムの全体像

ここまでで、安定微係数の計算方法、座標系、Thin/Thick、Python API の解析入力を整理した。

ここからは、これらの設定を AnalysisVSPAERO.py でどのように実行するかを確認する。

対象ファイル

解析に使用する入力ファイル、実行用スクリプト、VSPAERO が生成する出力ファイル、収束性確認用 Notebook の役割は次の通り。

ファイル	役割
AnalysisVSPAERO.py	VSPAERO 解析の実行
G103A.vsp3	OpenVSP モデル
G103A.vspaero	VSPAERO へ渡す最終設定
G103A.stab	基準ケース、摂動ケース、微係数
G103A.flt	飛行力学形式の微係数
run_wake_iter_sensitivity.ipynb	WakeNumIter 感度解析
plot_wake_iter_sensitivity.ipynb	感度解析結果の可視化

処理フロー

プログラムは、モデルの読み込みから結果の取得までを次の順番で実行する。

事前検証を geometry 生成より前に行うことで、Geom 名や Set 構成の不備を solver 実行前に検出する。

1. .vsp3 を読み込む
2. Geom、SubSurface、Set、Control Surface Group、reference wing を検証する
3. VSPAEROComputeGeometry を実行する
4. VSPAEROSweep の入力を設定する
5. UnsteadyType = STABILITY_DEFAULT を設定する
6. VSPAERO を実行する
7. Results Manager から結果を取得する
8. .vspaero、.stab、.flt を保存する
9. 必要に応じて .stab を後続の飛行力学解析へ渡す

関数構成

主な関数は次の通り。

validate_vsp3_for_stability_derivatives()
vsp_stability_derivatives()

validate_vsp3_for_stability_derivatives() は、G103A モデルの命名規約と解析設定を確認する。

vsp_stability_derivatives() は、geometry 生成、解析入力設定、solver 実行、結果取得を上から順に行う。

モデルの事前検証

安定微係数解析を実行する前に、モデルの Geom、SubSurface、Thin/Thick Set、Control Surface Group、reference wing が想定通りに設定されていることを確認する。

Geom

G103A では次の Geom を要求する。

FuselageGeom
WingGeom
HTailGeom
VTailGeom

SubSurface

WingGeom : AILERON
HTailGeom : ELEVATOR
VTailGeom : RUDDER

Thin/Thick Set

ThickGeom : FuselageGeom
ThinGeom  : WingGeom, HTailGeom, VTailGeom

Control Surface Group

AILERON_GROUP
ELEVATOR_GROUP
RUDDER_GROUP

reference wing

reference wing は WingGeom である。

G103A の .vsp3 には WingGeom の Geom ID が RefGeomID として保存されている。

安定微係数解析の実行

事前検証に合格したモデルに対して、VSPAERO 用 geometry を生成し、単一の飛行条件を VSPAEROSweep へ設定して安定微係数解析を実行する。

処理は VSPAEROComputeGeometry、VSPAEROSweep の初期化、解析条件の設定、solver 実行、結果取得の順に進む。

VSPAEROComputeGeometry

最初に VSPAERO 用 geometry を生成する。

compgeom_name = "VSPAEROComputeGeometry"

vsp.SetAnalysisInputDefaults(compgeom_name)
vsp.SetIntAnalysisInput(compgeom_name, "GeomSet", [thick_set_index], 0)
vsp.SetIntAnalysisInput(compgeom_name, "ThinGeomSet", [thin_set_index], 0)

compgeom_result_id = vsp.ExecAnalysis(compgeom_name)

VSPAEROSweep

単一条件の安定微係数解析も VSPAEROSweep から実行する。

analysis_name = "VSPAEROSweep"
vsp.SetAnalysisInputDefaults(analysis_name)

解析条件の設定

VSPAEROSweep を初期化したら、Thin/Thick Set、迎角、Mach 数、Reynolds 数、Stability Type を設定する。

本記事では 1 つの基準状態まわりの微係数を計算するため、迎角、Mach 数、Reynolds 数の各点数を 1 とする。

vsp.SetIntAnalysisInput(analysis_name, "GeomSet", [thick_set_index], 0)
vsp.SetIntAnalysisInput(analysis_name, "ThinGeomSet", [thin_set_index], 0)

vsp.SetDoubleAnalysisInput(analysis_name, "AlphaStart", [alpha], 0)
vsp.SetDoubleAnalysisInput(analysis_name, "AlphaEnd", [alpha], 0)
vsp.SetIntAnalysisInput(analysis_name, "AlphaNpts", [1], 0)

vsp.SetDoubleAnalysisInput(analysis_name, "MachStart", [mach], 0)
vsp.SetDoubleAnalysisInput(analysis_name, "MachEnd", [mach], 0)
vsp.SetIntAnalysisInput(analysis_name, "MachNpts", [1], 0)

vsp.SetDoubleAnalysisInput(analysis_name, "ReCref", [reynolds], 0)
vsp.SetDoubleAnalysisInput(analysis_name, "ReCrefEnd", [reynolds], 0)
vsp.SetIntAnalysisInput(analysis_name, "ReCrefNpts", [1], 0)

vsp.SetIntAnalysisInput(
    analysis_name,
    "UnsteadyType",
    [vsp.STABILITY_DEFAULT],
    0,
)

solver 実行

解析条件の設定が完了したら、ExecAnalysis() を呼び出して VSPAERO を実行する。

result_id = vsp.ExecAnalysis(analysis_name)

返された result_id は、Results Manager から解析結果を取得するときに使用する。

RedirectFile で solver 出力先を指定できる。

StopBeforeRun = True では入力ファイル生成後、solver 実行前で停止する。

結果の取得

solver 実行後は、result_id を用いて Results Manager から data name、double result、string result を取得する。

解析結果を Python 上で処理するとともに、VSPAERO が生成したファイルも保存する。

ファイルとしては、主に次を使用する。

G103A.vspaero
G103A.stab
G103A.flt

.vspaero は solver へ渡した最終設定、.stab は基準・摂動ケースと微係数、.flt は飛行力学向けの係数名を保持する。

Control Surface Group の扱い

Control Surface Group 微係数の符号は、group command と各舵面の Gain で決まる。

ここでは G103A に設定したエルロン、エレベータ、ラダーの group 構成を確認し、.stab の ConGrp_N との対応を整理する。

G103A の Gain は次の通り。

AILERON_GROUP
├─ WingGeom_Surf0_AILERON : Gain =  1
└─ WingGeom_Surf1_AILERON : Gain =  1

ELEVATOR_GROUP
├─ HTailGeom_Surf0_ELEVATOR : Gain =  1
└─ HTailGeom_Surf1_ELEVATOR : Gain = -1

RUDDER_GROUP
└─ VTailGeom_Surf0_RUDDER : Gain = -1

左右 surface の幾何学的な正方向には reflection が含まれる。

group command の正方向は、各 surface の Gain を含めて確認する。

各 group の基準 command deflection は 0 deg である。

.stab の ConGrp_N との対応は以下の通り

delta_a -> ConGrp_1 -> AILERON_GROUP
delta_e -> ConGrp_2 -> ELEVATOR_GROUP
delta_r -> ConGrp_3 -> RUDDER_GROUP

使い方

環境

OpenVSP 本体と、同じリリースに対応する Python API を使用する。

環境構築はこちら。

【Windows】OpenVSP Python APIのインストール

Windows11にOpenVSPのPython APIをインストールする

mtkbirdman.com

2024.11.30

ディレクトリ構成

AnalysisVSPAERO.py と共通モジュールを bin へ置き、解析対象の .vsp3 と実行スクリプトを機体ごとのディレクトリへ置く。

この記事で想定するディレクトリ構成は次の通り。

./
├─ bin/
│  ├─ AnalysisVSPAERO.py
│  ├─ util.py
│  └─ ISAspecification.py
├─ G103A/
│  ├─ G103A.vsp3
│  └─ test_stability.py
└─ results/

モデルの準備

Python スクリプトを実行する前に、OpenVSP の GUI で解析対象、Control Surface Group、reference quantities を設定する。

G103A モデルでは次の項目を準備する。

1. Geom 名を確認する
2. AILERON、ELEVATOR、RUDDER の SubSurface を作成する
3. ThickGeom と ThinGeom を作成する
4. Control Surface Group を作成する
5. Gain と基準舵角を設定する
6. WingGeom を reference wing に設定する
7. \(S_{\mathrm{ref}}\)、\(b_{\mathrm{ref}}\)、\(c_{\mathrm{ref}}\)、CG を確認する

関数の引数

vsp_stability_derivatives() には解析対象の .vsp3 を必須引数として与え、飛行条件と solver 設定をキーワード引数で指定する。

関数定義は次の通り。

def vsp_stability_derivatives(
    vsp3_path,
    *,
    alpha=2.0,
    mach=0.1,
    reynolds=4.4e6,
    verbose=1,
    vspaero_verbose=0,
    ncpu=None,
    wake_num_iter=None,
    fixed_wake_flag=None,
    redirect_file="",
    stop_before_run=False,
):

引数	内容
vsp3_path	解析対象の .vsp3
alpha	迎角、deg
mach	Mach 数
reynolds	\(Re_{c_{\mathrm{ref}}}\)
verbose	Python 側の表示レベル
vspaero_verbose	VSPAERO 側の表示レベル
ncpu	CPU 数
wake_num_iter	wake 反復回数
fixed_wake_flag	fixed wake の指定
redirect_file	solver 出力先
stop_before_run	solver 実行前で停止するフラグ

None を指定した optional input は、.vsp3 の VSPAEROSettings を維持する。

返り値

返り値には、解析の成否、Result ID、取得した結果、警告、エラー、出力ファイル、実行時間を含む。

正確な key は使用中の AnalysisVSPAERO.py の return statement を確認する。

実行例

G103A を迎角 2 deg、Mach 数 0.1、基準コード Reynolds 数 \(4.4\times 10^6\) で解析する例は次の通り。

from AnalysisVSPAERO import vsp_stability_derivatives

result = vsp_stability_derivatives(
    "G103A.vsp3",
    alpha=2.0,
    mach=0.1,
    reynolds=4.4e6,
    ncpu=4,
    wake_num_iter=10,
    fixed_wake_flag=False,
    redirect_file="",
    stop_before_run=False,
)

出力ファイル

解析を実行すると、入力した .vsp3 と同じ basename を持つ VSPAERO 出力ファイルが生成される。

安定微係数の確認には主に .stab を使用し、solver 条件や収束履歴を確認するときは .vspaero と .history も参照する。

ファイル	内容
.vspaero	solver 入力条件
.stab	基準・摂動ケース、微係数
.flt	飛行力学向け微係数
.history	solver 履歴
.polar	polar result

WakeNumIter の収束性確認

WakeNumIter は、steady wake を更新して流れ場を収束させる反復回数である。

wake 形状と循環分布は相互に影響するため、反復回数が少ない場合は両者が十分に整合しない可能性がある。

安定微係数は、基準ケースと摂動ケースの小さな係数差を微小な摂動量で除して求める。

このため、基準空力係数では目立たない未収束成分が、安定微係数や舵効きの微係数では増幅される場合がある。

そこで、run_wake_iter_sensitivity.ipynb から WakeNumIter を変更しながら vsp_stability_derivatives() を実行し、基準空力係数、安定微係数、舵効きの微係数、実行時間を比較する。

確認範囲は、

\begin{align} 4\leq N_{\mathrm{wake}}\leq 40 \end{align}

とした。

解析した 23 ケースはすべて正常に終了している。

基準空力係数

はじめに、基準ケースの \(C_L\)、\(C_D\)、\(C_Y\)、\(C_l\)、\(C_m\)、\(C_n\) を確認する。

\(N_{\mathrm{wake}}=4\) から 40 までの範囲では、いずれの基準空力係数も変化は小さい。

特に \(C_L\) と \(C_D\) は、最小の反復回数からおおむね一定である。

\(C_m\) は基準空力係数の中では比較的大きく変化しているものの、縦軸の範囲に対しては小さい。

全 23 ケースにおける最小値、最大値、変動幅は次の通り。

係数	最小値	最大値	変動幅
\(C_L\)	0.5344976	0.5349229	0.0004253
\(C_D\)	0.0149824	0.0150022	0.0000198
\(C_Y\)	0.0001798	0.0001996	0.0000198
\(C_l\)	0.0002635	0.0002763	0.0000128
\(C_m\)	-0.0586810	-0.0579700	0.0007110
\(C_n\)	-0.0000016	0.0000038	0.0000054

\(N_{\mathrm{wake}}\geq 10\) に限定すると、変動幅はさらに小さくなる。

係数	\(N_{\mathrm{wake}}\geq 10\) における変動幅
\(C_L\)	0.0000946
\(C_D\)	0.0000077
\(C_Y\)	0.0000103
\(C_l\)	0.0000010
\(C_m\)	0.0000769
\(C_n\)	0.0000024

基準空力係数だけを見ると、比較的小さな WakeNumIter でも収束しているように見える。

ただし、この結果だけでは安定微係数解析の収束性を判断できない。

基準ケースと摂動ケースの差に残る未収束成分は、次に示す微係数で明確になる。

安定微係数・舵効きの微係数

次に、縦・横・方向の安定微係数と舵効きの微係数を確認する。

\(N_{\mathrm{wake}}<10\) では、複数の安定微係数と舵効きの微係数に大きな変動が残っている。

特に、\(C_{m\hat{q}}\)、\(C_{m\beta}\)、\(C_{m\alpha}\)、\(C_{L\alpha}\)、\(C_{L\hat{q}}\)、\(C_{L\beta}\) では、反復回数の小さいケースが後続の値から大きく外れる。

横滑り角、角速度、舵角に対する絶対値の小さい交差微係数では、符号の反転や相対的に大きな振動も見られる。

したがって、\(N_{\mathrm{wake}}<10\) の結果は、安定微係数や舵効きの微係数を用いる定量評価には不十分である。

\(N_{\mathrm{wake}}\geq 10\) では、反復回数の小さい領域に見られた大きな変動は減少し、主要な微係数の符号もおおむね安定する。

代表的な微係数について、\(N_{\mathrm{wake}}\geq 10\) における範囲を次に示す。

微係数	最小値	最大値	変動幅
\(C_{L\alpha}\)	5.7045420	5.7432217	0.0386797
\(C_{D\alpha}\)	0.1489467	0.1523004	0.0033537
\(C_{m\alpha}\)	-1.0907280	-0.9502114	0.1405166
\(C_{Y\beta}\)	-0.2074396	-0.1708539	0.0365857
\(C_{l\beta}\)	-0.0666763	-0.0604742	0.0062021
\(C_{n\beta}\)	0.0266749	0.0355092	0.0088343
\(C_{l\hat{p}}\)	-0.6672625	-0.6650436	0.0022189
\(C_{m\hat{q}}\)	-24.5824014	-24.1637981	0.4186033
\(C_{n\hat{r}}\)	-0.0417649	-0.0340175	0.0077474
\(C_{l\delta_a}\)	-0.3743014	-0.3732662	0.0010352
\(C_{m\delta_e}\)	-1.5175926	-1.5034014	0.0141912
\(C_{n\delta_r}\)	-0.0399974	-0.0375290	0.0024684

\(C_{l\hat{p}}\)、\(C_{l\delta_a}\)、\(C_{m\delta_e}\) は比較的狭い範囲に収まる。

\(C_{m\alpha}\)、\(C_{m\hat{q}}\)、\(C_{n\hat{r}}\) には、\(N_{\mathrm{wake}}\geq 10\) でも確認できる変動が残っている。

絶対値が小さい微係数では、変動幅そのものが小さくても相対変動が大きくなるため、後続計算への影響を個別に確認する必要がある。

実行時間

実行時間は概ね WakeNumIter とともに増加する。

WakeNumIter	VSPAERO 実行時間の目安
4	約 145 s
10	約 339 s
20	約 638 s
40	約 1,258 s

\(N_{\mathrm{wake}}=40\) の実行時間は、\(N_{\mathrm{wake}}=10\) の約 3.7 倍である。

反復回数を大きくすればよいとは限らないため、必要な精度と実行時間の両方を考慮する。

今回の G103A モデルでは、\(N_{\mathrm{wake}}=10\) を安定微係数解析に用いる最低限の目安とした。

おわりに

OpenVSP の Python API を用いて、VSPAERO の STABILITY_DEFAULT による安定微係数解析を行う方法について説明した。

↓関連記事

【まとめ】OpenVSP入門

NASAが開発したオープンソースの航空機概念設計ソフト「OpenVSP」についての記事の一覧

mtkbirdman.com

2025.01.22