tmp

ラダーのみを用いて旋回する航空機について、垂直尾翼容積比 \(V_v\) と等価上反角 \(\Gamma_{\mathrm{eff}}\) を軸にした設計チャートの作り方を整理する。

スポンサーリンク

はじめに

ラダーのみを用いて旋回する航空機では、ラダー入力によってヨーイング角速度と横滑りを作り、ヨーとロールの連成項や上反角効果を通じてロールを作る。

したがって、垂直尾翼・上反角の設計上は、スパイラル安定性や Dutch roll 特性だけでなく、旋回性能、横転性能、横風突風応答性能、抵抗性能などをまとめて見つつ、バランスの良い垂直尾翼と上反角を設定する必要がある。

そこで、本記事では、ラダーのみを用いて旋回する航空機の垂直尾翼・上反角設計について、横軸に垂直尾翼容積比 \(V_v\)、縦軸に等価上反角 \(\Gamma_{\mathrm{eff}}\) を取った \(V_v-\Gamma_{\mathrm{eff}}\) チャートを用いて検討してみる。

↓対象とする機体はこちら

鳥コン滑空機を想定したSampleGliderのOpenVSPモデル
鳥コン滑空機を想定したSampleGliderのOpenVSPモデルおよび数値例を示す。

↓最初から理論を追いたい人はこちらから

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

それではいってみよう。

\(V_v-\Gamma_{\mathrm{eff}}\) チャートの考え方

\(V_v-\Gamma_{\mathrm{eff}}\) チャートでは、垂直尾翼と上反角の組合せを

横軸に垂直尾翼容積比:

\begin{align} V_v=\frac{S_vl_v}{Sb} \end{align}

縦軸に等価上反角:

\begin{align} \Gamma_{\mathrm{eff}} =\frac{ \displaystyle\int_0^1 x\sqrt{1-x^2}\,\theta_b(x)\,dx }{ \displaystyle\int_0^1 x\sqrt{1-x^2}\,dx } \end{align}

を取った 2 次元平面として考える。

ここで、等価上反角 \(\Gamma_{\mathrm{eff}}\) は、スパン方向に分布するたわみ角を、ローリングモーメントアームで重みづけした平均値として定義する。

\(V_v-\Gamma_{\mathrm{eff}}\) チャート上で、スパイラル安定性、旋回性能、横転性能、横風突風応答性能、抵抗性能などを重ねて見ることで、ラダーのみ旋回に必要な複数の条件を同時に満たす垂直尾翼と上反角の存在領域を探すことができる。

OpenVSP による設計空間の生成

\(V_v-\Gamma_{\mathrm{eff}}\) チャートを作成する流れは以下の通り。

  • OpenVSP の基準形状を準備する。
  • 垂直尾翼容積比 \(V_v\) を変化させた形状を作る。
  • 翼端たわみ量 \(w_{\mathrm{tip}}\) から主翼のたわみ角分布を作る。
  • たわみ角分布を OpenVSP の Dihedral へ反映する。
  • VSPAERO で安定微係数を取得する。
  • .stab ファイルからラダーのみ旋回の評価指標を計算する。
  • \(V_v-\Gamma_{\mathrm{eff}}\) 平面に等値線を描く。

1 つずつ順番に説明していく。

基準形状に必要な要件

\(V_v-\Gamma_{\mathrm{eff}}\) チャートを作るには、まず OpenVSP の基準形状が必要である。

本記事では、次のような基準形状を想定する。

  • 胴体、主翼、水平尾翼、垂直尾翼を持つ。
  • 主翼、水平尾翼、垂直尾翼は VSPAERO の薄翼解析対象として扱う。
  • ラダー、エレベータ、エルロンの Control Surface Group が定義されている。
  • VSPAERO の Reference Geom が主翼に設定されている。
  • 安定微係数解析で使う Set が整理されている。

AnalysisVSPAERO.py では、上記のモデルを検証するため、Geom 名、Geom Type、SubSurface、Set、Control Surface Group、Reference Geom について事前確認する実装になっている。

具体的には、次のモデルを想定する。

要素
胴体FuselageGeom
主翼WingGeom
水平尾翼HTailGeom
垂直尾翼VTailGeom
エルロンAILERON_GROUP
エレベータELEVATOR_GROUP
ラダーRUDDER_GROUP
Reference GeomWingGeom

VSPAERO の安定微係数解析では、Control Surface に対する微係数も計算するので、Control Surface Group の対応が崩れるとチャートの評価指標も崩れてしまう。

この記事の内容を試してみたい人は、とりあえずこの記事の例の vsp3 モデルをベースに、Geom の名前などは変えずに寸法だけ変えていくのが良いと思われる。

基準形状で確認する量

基準形状から、次の量を取得する。

\begin{align} S &= S_{\mathrm{ref}}, \\ b &= b_{\mathrm{ref}}, \\ l_v &= x_{\mathrm{ac},v}-x_{\mathrm{cg}}, \\ S_v &= \text{垂直尾翼面積} \end{align}

ここで、\(S\) は基準主翼面積、\(b\) は基準翼幅、\(l_v\) は重心から垂直尾翼空力中心までの距離である。

垂直尾翼容積比は

\begin{align} V_v=\frac{S_vl_v}{Sb} \end{align}

で定義する。

この式を \(S_v\) について解くと、目標の垂直尾翼容積比 \(V_v\) に対して必要な垂直尾翼面積は

\begin{align} S_v=\frac{V_vSb}{l_v} \end{align}

である。

垂直尾翼容積比 \(V_v\) を変化させる方法

VvGammaChart.py の update_vtail_volume() は、基準形状の主翼面積 \(S\)、主翼翼幅 \(b\)、重心位置 \(x_{\mathrm{cg}}\)、垂直尾翼位置から \(l_v\) を求め、

  • アスペクト比固定
  • スパン固定

の 2 通りの方法で、\(V_v\) に対応する垂直尾翼面積へ形状を決定する。

目標面積は

\begin{align} S_{v,\mathrm{target}}=\frac{V_vSb}{l_v} \end{align}

である。

垂直尾翼の現状面積を \(S_{v,0}\) とすれば、面積倍率 \(\lambda_S\) は

\begin{align} \lambda_S=\frac{S_{v,\mathrm{target}}}{S_{v,0}} \end{align}

である。

アスペクト比固定の場合は、スパン方向倍率 \(\lambda_b\) とコード方向倍率 \(\lambda_c\) を同じにして拡大縮小する。

垂直尾翼面積は、代表的にスパンとコードの積に比例するので、

\begin{align} \lambda_S=\lambda_b\lambda_c \end{align}

である。

アスペクト比を保つために

\begin{align} \lambda_b=\lambda_c \end{align}

とおくと、

\begin{align} \lambda_S &= \lambda_b^2, \\ \lambda_b &= \sqrt{\lambda_S}, \\ \lambda_c &= \sqrt{\lambda_S} \end{align}

となる。

この場合、アスペクト比は

\begin{align} AR_{v,\mathrm{new}} &=\frac{(\lambda_b b_v)^2}{\lambda_S S_v} \\ &=\frac{\lambda_S b_v^2}{\lambda_S S_v} \\ &=AR_{v,0} \end{align}

となり、基準形状の垂直尾翼アスペクト比を保ったまま垂直尾翼面積を変更できる。

一方、スパン固定の場合は、垂直尾翼のスパンを変えずにコードだけを変える。

したがって、

\begin{align} \lambda_b=1 \end{align}

とすると、面積倍率 \(\lambda_S\) は

\begin{align} \lambda_S=\lambda_b\lambda_c=\lambda_c \end{align}

となる。

この切り替えは、update_vtail_volume() の vtail_area_scale_mode で指定でき、vtail_area_scale_mode="fixed_aspect_ratio" では

\begin{align} \lambda_b &= \sqrt{\lambda_S}, \\ \lambda_c &= \sqrt{\lambda_S} \end{align}

vtail_area_scale_mode="fixed_span" では

\begin{align} \lambda_b &= 1, \\ \lambda_c &= \lambda_S \end{align}

となる。

このようにして、チャートの横軸である \(V_v\) を直接指定し、その値に対応する OpenVSP 形状を生成する。

等価上反角

長スパン・弾性翼では、荷重によって翼が上方にたわむ。

たわみ角を \(\theta_b(y)\) とすれば、小角近似では、弾性たわみによる上反角は次で表される。

\begin{align} \Gamma_{\mathrm{elastic}}(y) \simeq \frac{dw}{dy} =\theta_b(y) \end{align}

元の幾何上反角を \(\Gamma_{\mathrm{rigid}}(y)\) とすれば、たわみ後の局所実効上反角は

\begin{align} \Gamma_{\mathrm{eff}}(y) =\Gamma_{\mathrm{rigid}}(y) +\theta_b(y) \end{align}

である。

長スパン・弾性翼では、主翼のたわみにより上反角がスパン方向に分布するため、チャートの縦軸として使える代表的な上反角を単一の値として表す必要がある。

この記事では、それを等価上反角 \(\Gamma_{\mathrm{eff}}\) とし、翼端たわみ量 \(w_{\mathrm{tip}}\) から作ったたわみ角分布 \(\theta_b(x)\) を、楕円平面形に対応する重みで平均した値として定義する。

横滑り角 \(\beta\) と上反角 \(\Gamma(y)\) による局所揚力変化を

\begin{align} d(\Delta L) \simeq q_\infty c(y)a(y)\beta\Gamma(y)\,dy \end{align}

と見る。

この揚力差がモーメントアーム \(y\) を持つため、ローリングモーメント係数への寄与は

\begin{align} C_{l\beta,\Gamma} \simeq -\frac{2}{Sb} \int_0^{b/2} y c(y)a(y)\Gamma(y)\,dy \end{align}

である。

同じ \(C_{l\beta,\Gamma}\) を与える一様上反角 \(\Gamma_0\) は

\begin{align} \Gamma_0 =\frac{ \displaystyle\int_0^{b/2}y c(y)a(y)\Gamma(y)\,dy }{ \displaystyle\int_0^{b/2}y c(y)a(y)\,dy } \end{align}

と定義できる。

つまり、等価上反角は、局所上反角分布 \(\Gamma(y)\) をロールモーメントアーム \(y\) と局所空力寄与 \(c(y)a(y)\) で重み付けした平均である。

今回は、楕円コード長分布・桁径一定を仮定して、翼端たわみ量 \(w_{\mathrm{tip}}\) からたわみ形状を作成し、等価上反角を計算した。

詳細な計算式は、Appendix を参照してほしい。

安定微係数の計算

安定微係数の詳しい計算方法は、別記事「OpenVSP の Python API で動安定微係数を計算する」に譲る。

本記事では、OpenVSP / VSPAERO により .stab ファイルを生成し、その .stab ファイルから必要な微係数を読み取るものとして扱う。

主に使う微係数は次である。

\begin{align} C_{l\beta}, C_{n\beta}, C_{l\hat p}, C_{l\hat r}, C_{n\hat p}, C_{n\hat r}, C_{l\delta_r}, C_{n\delta_r} \end{align}

TrimTurnSolver.py では、.stab ファイルを読み取り、Base Aero、安定微係数表、Control Surface Group を取り出す。

プログラムの説明

\(V_v-\Gamma_{\mathrm{eff}}\) チャート作成プログラムは、主に VvGammaChart.py に実装されている。

処理は、次の段階に分けられる。

基準形状の読み取り

read_reference_wing_summary() と read_wing_section_table() により、基準主翼の面積、翼幅、断面情報を取得する。

主翼の面積と翼幅は、垂直尾翼容積比

\begin{align} V_v=\frac{S_vl_v}{Sb} \end{align}

の基準量として使う。

垂直尾翼容積比の変更

update_vtail_volume() により、指定した \(V_v\) から目標垂直尾翼面積を求め、垂直尾翼形状を更新する。

\begin{align} S_{v,\mathrm{target}} =\frac{V_vSb}{l_v} \end{align}

翼端たわみ量から上反角を作る

elastic_wing_deflection_distribution() により、翼端たわみ量 \(w_{\mathrm{tip}}\) からたわみ角分布 \(\theta_b(x)\) を作る。

apply_wing_deflection_as_dihedral() により、得られたたわみ角分布を OpenVSP の Dihedral へ近似的に反映する。

VSPAERO 安定微係数解析

run_vspaero_stability_case() により、生成した .vsp3 に対して VSPAERO の安定微係数解析を実行する。

基準形状の検証と安定微係数解析については、AnalysisVSPAERO.py の validate_vsp3_for_stability_derivatives() と vsp_stability_derivatives() の考え方に従う。

後処理

postprocess_vv_gamma_cases() により、各ケースの .stab から次の量を計算する。

  • \(w_{\mathrm{tip}}\) から計算した \(\Gamma_{\mathrm{eff}}\)
  • \(C_{l\beta}\)、\(C_{n\beta}\) などの安定微係数
  • spiral_margin
  • 旋回性能指標
  • 横転性能指標
  • 横風突風応答指標

可視化

plot_vv_gamma_contour_panel() により、\(V_v-\Gamma_{\mathrm{eff}}\) 平面上に各指標の等値線を描く。

使い方

test_vv_gamma_chart.ipynb.py と plot_vv_gamma_chart.ipynb.py では、notebook からチャート作成処理を呼び出す例が示されている。

基本的な使い方は次の流れである。

入力条件を設定する

まず、基準 OpenVSP ファイル、出力ディレクトリ、\(V_v\) の候補、翼端たわみ量の候補を設定する。

base_vsp3_path = Path('../../models/BRGlider/BRGlider.0G.vsp3')
output_dir = Path('')

vv_values = np.linspace(0.001, 0.006, 11)
tip_deflections = np.linspace(0.0, 2.0, 21)

次に、VSPAERO 解析条件を設定する。

flight_condition = {
    'alpha_deg': 0.0,
    'mach': 0.0,
    'reynolds': 1.0e6,
}

さらに、形状変更に必要な設定を与える。

geometry_config = {
    'lv': 4.5,
    'wing_area': 18,
    'wing_span': 27,
    'xcg': 1.4,
    'wing_name': 'WingGeom',
    'vtail_name': 'VTailGeom',
    'n_span': 101,
    'vtail_area_scale_mode': 'fixed_aspect_ratio',
}

VSPAERO 安定微係数スイープを実行する

sweep = run_vv_wtip_stability_sweep(
    base_vsp3_path,
    vv_values,
    tip_deflections,
    flight_condition,
    geometry_config,
    output_dir,
    validate_base_model=True,
    fixed_wake_flag=False,
    wake_num_iter=10,
    ncpu=8,
    verbose=1,
)

この処理により、\(V_v\) と翼端たわみ量の組合せごとに .vsp3 と .stab が作成される。

1 形態 4~5 分で、200 形態ほど計算すると 12 時間くらいかかる。

後処理を行う

既存の .stab 群から、等価上反角、安定微係数、旋回トリム、6DOF ロール応答指標を計算する。

results = postprocess_vv_gamma_cases(
    sweep,
    mass=mass,
    inertia=inertia,
    delta_r=math.radians(10.0),
    target_delta_phi=math.radians(2.0),
    turn_trim_mode='gliding',
    turn_trim_phi=math.radians(2.0),
    output_csv_path=output_dir / 'vv_gamma_metrics.csv',
    history_output_dir=output_dir / '6dof_history',
    write_6dof_history=False,
    plot_6dof_history=False,
    verbose=1,
)

だいたい 1 形態 10~15 秒で、200 形態計算するのに 30 分ほどかかる。

チャートを描く

後処理結果 vv_gamma_metrics.csv を読み込み、必要な列を等値線として描く。

fig, ax = plt.subplots()
plot_vv_gamma_contour(results, 'spiral_margin', ax=ax)
plt.show()

複数のパネルへ 1 指標ずつ描く場合は、plot_vv_gamma_contour_panel() を使う。

複数指標の等値線を 1 つの \(V_v-\Gamma_{\mathrm{eff}}\) 平面へ重ねる場合は、plot_vv_gamma_contour_lines() を使う。

考察

ここからは、実際に \(V_v-\Gamma_{\mathrm{eff}}\) チャートを使って、ラダーのみで旋回する航空機の垂直尾翼・上反角の設計について考えてみる。

この章で確認する指標は以下の通り。

見る項目主に見る列何を判断するか
安定微係数・ラダー効き\(C_{Y\beta}\)、\(C_{l\beta}\)、\(C_{n\beta}\)、\(C_{l\hat r}\)、\(C_{n\hat r}\)、\(C_{l\delta_r}\)、\(C_{n\delta_r}\)ラダー入力、横滑り角、ヨーレート、ロールレートがどの経路でつながるか
安定性\(C_{l\beta}C_{n\hat r}-C_{n\beta}C_{l\hat r}\)、\(\delta_{r,\mathrm{trim}}\)スパイラル安定性と、小バンクのラダーのみ定常旋回が成立するか
旋回性能\(\phi
横転性能\(\dfrac{\Delta\phi/\Delta t}{\delta_r}\)、\(K_{\mathrm{rudder,roll}}=(C_{l\hat r}-C_{l\beta})C_{n\delta_r}\)ラダー入力後に、有限時間でどれだけバンク角を作れるか
横風突風応答\(\dfrac{\Delta\phi(\tau_e)}{\beta_{g,\max}}\)、\(K_{\mathrm{gust,roll}}=\dfrac{8I_z}{\rho Sb^3}C_{l\beta}+C_{l\hat r}C_{n\beta}\)横風外乱に対してロールやヨーがどれだけ出やすいか

まず安定微係数で各経路の強さを確認し、その後に安定性、旋回性能、横転性能、横風突風応答を順に見ていく。

安定微係数・ラダー効き

スパイラル安定性、旋回性能、横転性能、横風突風応答、抵抗性能などの複合指標を見る前に、各安定微係数とラダー効きを \(V_v-\Gamma_{\mathrm{eff}}\) チャート上で確認する。

横・方向の線形空力係数は、入力を \(\beta\)、\(\hat p\)、\(\hat r\)、\(\delta_r\)、出力を \(C_Y\)、\(C_l\)、\(C_n\) として整理する。

入力横力ローリングモーメントヨーイングモーメント
横滑り角 \(\beta\)\(C_{Y\beta}\)\(C_{l\beta}\)\(C_{n\beta}\)
ロールレート \(\hat p\)\(C_{Y\hat p}\)\(C_{l\hat p}\)\(C_{n\hat p}\)
ヨーレート \(\hat r\)\(C_{Y\hat r}\)\(C_{l\hat r}\)\(C_{n\hat r}\)
ラダー舵角 \(\delta_r\)\(C_{Y\delta_r}\)\(C_{l\delta_r}\)\(C_{n\delta_r}\)

それぞれの設計上の役割は、次のように整理できる。

主に効く設計変数設計上の意味
\(C_{Y\beta}\)\(V_v\)、\(\Gamma_{\mathrm{eff}}\)横滑り角に対する横力。旋回時の横力のつり合いと突風応答に効く
\(C_{l\beta}\)\(\Gamma_{\mathrm{eff}}\)横滑り角をロールへ変換する上反角効果。横転性能と横風突風応答に効く
\(C_{n\beta}\)\(V_v\)方向静安定。スパイラル安定性、旋回性能、横風突風応答に効く
\(C_{Y\hat p}\)\(\Gamma_{\mathrm{eff}}\)、\(V_v\)ロールレートに伴う横力。主翼上反角と垂直尾翼の上下位置に依存する
\(C_{l\hat p}\)主翼スパン、揚力傾斜、上反角分布ロール減衰。有限時間のロール応答を抑える
\(C_{n\hat p}\)主翼空気力、上反角分布、垂直尾翼位置ロールレートに伴うヨーイングモーメント
\(C_{Y\hat r}\)\(V_v\)、\(\Gamma_{\mathrm{eff}}\)ヨーレートに伴う横力。横滑り運動とヨー運動を結ぶ
\(C_{l\hat r}\)主翼揚力分布、スパン、上反角分布ヨーレートをロールへ変換する経路。ラダーのみ旋回で重要
\(C_{n\hat r}\)\(V_v\)ヨー減衰。スパイラル安定性と突風応答に効く
\(C_{Y\delta_r}\)\(V_v\)、ラダー効きラダー入力による横力
\(C_{l\delta_r}\)垂直尾翼の上下位置、ラダー効きラダー入力による直接ロール
\(C_{n\delta_r}\)\(V_v\)、ラダー効きラダー入力によるヨーイングモーメント

定性的な妥当性を確認するための推算式は、以下の記事を参考にする。

鳥コン滑空機における横・方向の微小擾乱理論
鳥コン滑空機のような大きな上反角を持つ機体に対する横・方向の微小擾乱理論について説明する

参照記事では \(c_y\)、\(p\)、\(r\) を用いるが、本記事と csv では、\(c_y\) を \(C_Y\)、\(p\) と \(r\) に関する微係数を、

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

に対する微係数として読む。

また、参照記事の \(S_f\)、\(l_f\)、\(z_f\)、\(a_f\) は、それぞれ垂直尾翼面積、重心から垂直尾翼までの前後方向距離、上下方向距離、垂直尾翼揚力傾斜を表す。\(a_1\)、\(c=c(y)\)、\(\Gamma=\Gamma(y)\)、\(z=z(y)\) は、主翼の局所揚力傾斜、コード長、上反角、上下位置である。

横滑り角による作用

今回のチャートでは、

\begin{align} -0.213923 &\le C_{Y\beta} \le -0.015450, \\ -0.222588 &\le C_{l\beta} \le 0.001263, \\ -0.023719 &\le C_{n\beta} \le 0.015016 \end{align}

である。

参照記事より、横滑り角に関する 3 係数は、

\begin{align} C_{Y\beta} &=-\frac{2}{S} \int_0^{b/2} a_1\sin^2\Gamma\,c\,dy -\frac{S_f}{S}a_f, \\ C_{l\beta} &=-\frac{2}{Sb} \int_0^{b/2} a_1\sin\Gamma \left( y\cos\Gamma+z\sin\Gamma \right)c\,dy -\frac{S_f}{S}a_f\frac{z_f}{b}, \\ C_{n\beta} &=-\frac{2}{Sb} \int_0^{b/2} \left( a_1\sin\alpha +C_L\cos\alpha +C_D\sin\alpha \right) \sin\Gamma\,y\,c\,dy +\frac{S_f}{S}a_f\frac{l_f}{b} \end{align}

である。

サイドウォッシュ微分を残す場合、\(C_{Y\beta}\) の垂直尾翼寄与は、

\begin{align} -\frac{S_f}{S}a_f \left( 1-\frac{\partial\sigma}{\partial\beta} \right) \end{align}

となる。

\(C_{Y\beta}\) の主翼寄与には \(\sin^2\Gamma\)、垂直尾翼寄与には \(S_f/S\) が入る。

このため、今回のチャートでは \(C_{Y\beta}\) が \(\Gamma_{\mathrm{eff}}\) と \(V_v\) の両方に感度を持つ。

\(\Gamma_{\mathrm{eff}}\) を増やすと主翼横力の寄与が増え、\(V_v\) を増やすと垂直尾翼横力の寄与が増える。

\(C_{l\beta}\) の主翼寄与には \(\sin\Gamma\) とローリングモーメントアーム \(y\cos\Gamma+z\sin\Gamma\) が入るため、今回のチャートでは \(\Gamma_{\mathrm{eff}}\) 感度が強い。

\(C_{n\beta}\) の垂直尾翼寄与には \(S_fa_fl_f/(Sb)\) が入り、今回のチャートでは \(V_v\) 感度が強い。

ロールレートによる作用

今回のチャートでは、

\begin{align} -0.451055 &\le C_{Y\hat p} \le -0.019076, \\ -0.776207 &\le C_{l\hat p} \le -0.771120, \\ -0.0878647 &\le C_{n\hat p} \le -0.0870005 \end{align}

である。

参照記事より、ロールレートに関する 3 係数は、

\begin{align} C_{Y\hat p} &=-\frac{4}{Sb} \int_0^{b/2} a_1y\sin\Gamma\cos\Gamma\,c\,dy -\frac{S_f}{S}a_f\frac{2z_f}{b}, \\ C_{l\hat p} &=-\frac{4}{Sb^2} \int_0^{b/2} a_1y\cos\Gamma \left( y\cos\Gamma+z\sin\Gamma \right)c\,dy -\frac{S_f}{S}a_f\frac{2z_f^2}{b^2}, \\ C_{n\hat p} &=-\frac{4}{Sb^2} \int_0^{b/2} \left( a_1\sin\alpha +C_L\cos\alpha +C_D\sin\alpha \right) \cos\Gamma\,y^2c\,dy +\frac{S_f}{S}a_f\frac{2z_fl_f}{b^2} \end{align}

である。

\(C_{Y\hat p}\) の主翼寄与には \(\sin\Gamma\cos\Gamma\) が入るため、今回のチャートでは \(\Gamma_{\mathrm{eff}}\) 感度が強い。

垂直尾翼寄与には \(S_fz_f/(Sb)\) が入るため、\(V_v\) の影響も加わる。

\(C_{l\hat p}\) の主翼寄与は、\(a_1\)、\(y\cos\Gamma\)、ローリングモーメントアーム \(y\cos\Gamma+z\sin\Gamma\) で決まる。

\(C_{n\hat p}\) の主翼寄与は、\(a_1\sin\alpha+C_L\cos\alpha+C_D\sin\alpha\)、\(\cos\Gamma\)、\(y^2\) で決まる。

今回の形状変更では、これらを支配する主翼平面形と飛行条件をほぼ維持しているため、両係数の変化幅は小さい。

ヨーレートによる作用

今回のチャートでは、

\begin{align} 0.008091 &\le C_{Y\hat r} \le 0.158519, \\ 0.220813 &\le C_{l\hat r} \le 0.224332, \\ -0.011632 &\le C_{n\hat r} \le -0.007140 \end{align}

である。

参照記事より、ヨーレートに関する 3 係数は、

\begin{align} C_{Y\hat r} &=\frac{8}{Sb} \int_0^{b/2} C_Ly\sin\Gamma\,c\,dy +\frac{S_f}{S}a_f\frac{2l_f}{b}, \\ C_{l\hat r} &=\frac{8}{Sb^2} \int_0^{b/2} C_Ly \left( y\cos\Gamma+z\sin\Gamma \right)c\,dy +\frac{S_f}{S}a_f\frac{2l_fz_f}{b^2}, \\ C_{n\hat r} &=\frac{8}{Sb^2} \int_0^{b/2} \left( C_L\sin\alpha-C_D\cos\alpha \right)y^2c\,dy -\frac{S_f}{S}a_f\frac{2l_f^2}{b^2} \end{align}

である。

\(C_{Y\hat r}\) の主翼寄与には \(\sin\Gamma\)、垂直尾翼寄与には \(l_f\) が入る。

\(C_{l\hat r}\) の主翼寄与は、主翼揚力分布とローリングモーメントアームで決まり、垂直尾翼寄与には \(l_fz_f\) が入る。

\(C_{n\hat r}\) の垂直尾翼寄与には \(-l_f^2\) が入り、ヨー減衰を作る。

このため、今回のチャートでは \(C_{Y\hat r}\) と \(C_{n\hat r}\) が \(V_v\) に感度を持つ。

\(C_{l\hat r}\) は主翼揚力分布と平面形の影響が大きく、今回のチャートでは変化幅が小さい。

ラダー舵角による作用

今回のチャートでは、

\begin{align} 0.010920 &\le C_{Y\delta_r} \le 0.060707, \\ -0.000612 &\le C_{l\delta_r} \le 0.002702, \\ -0.010628 &\le C_{n\delta_r} \le -0.001633 \end{align}

である。

参照記事より、ラダー舵角に関する 3 係数は、

\begin{align} C_{Y\delta_r} &=\frac{S_f}{S}a_f\tau, \\ C_{l\delta_r} &=\frac{S_f}{S}a_f\tau\frac{z_f}{b}, \\ C_{n\delta_r} &=-\frac{S_f}{S}a_f\tau\frac{l_f}{b} \end{align}

である。

\(C_{Y\delta_r}\) はラダーによる垂直尾翼横力、\(C_{l\delta_r}\) は横力と上下方向モーメントアーム \(z_f\)、\(C_{n\delta_r}\) は横力と前後方向モーメントアーム \(l_f\) で決まる。

したがって、今回のチャートでは \(C_{Y\delta_r}\) と \(C_{n\delta_r}\) が主として \(V_v\) に支配される。

\(C_{l\delta_r}\) は垂直尾翼とラダーの上下位置に依存し、今回の範囲では絶対値が小さい。

安定性

安定性では、スパイラルモードの安定性と、小バンクのラダーのみ定常旋回に必要なラダー舵角を確認する。

スパイラルモードの判別式は、

\begin{align} D_{\mathrm{sp}} =C_{l\beta}C_{n\hat r} -C_{n\beta}C_{l\hat r} \end{align}

である。

今回のチャートでは、

\begin{align} -0.003321 \le D_{\mathrm{sp}} \le 0.006916 \end{align}

となった。

\(C_{l\beta}\) は主として \(\Gamma_{\mathrm{eff}}\)、\(C_{n\beta}\) と \(C_{n\hat r}\) は主として \(V_v\) で変化するため、中立線 \(D_{\mathrm{sp}}=0\) は両設計変数の組合せで決まる。

小バンクのラダーのみ定常旋回に必要なラダー舵角は、6 自由度定常トリムから得た \(\delta_{r,\mathrm{trim}}\) で確認する。

今回、トリムが成立した 157 ケースでは、

\begin{align} -26.35\ \mathrm{deg} \le \delta_{r,\mathrm{trim}} \le 21.08\ \mathrm{deg} \end{align}

であった。

ラダーのみ定常旋回に必要なラダー舵角の簡易式は、

\begin{align} \delta_r \simeq -\frac{\Omega_0b}{2V_0} \frac{ C_{l\beta}C_{n\hat r} -C_{n\beta}C_{l\hat r} }{ C_{l\beta}C_{n\delta_r} } \end{align}

である。

必要な値がそろう 157 ケースでは、6 自由度定常トリムとの Pearson 相関係数は、

\begin{align} r\left( \delta_{r,\mathrm{trim}}, \delta_{r,\mathrm{simple}} \right) &=0.9988, \\ r\left( |\delta_{r,\mathrm{trim}}|, |\delta_{r,\mathrm{simple}}| \right) &=0.9962 \end{align}

であった。

簡易式による舵角の範囲は、

\begin{align} -28.30\ \mathrm{deg} \le \delta_{r,\mathrm{simple}} \le 19.84\ \mathrm{deg} \end{align}

であり、今回の設計平面では 6 自由度定常トリムの分布とよく対応している。

\(V_v\) を増やすと、主に \(C_{n\beta}\)、\(C_{n\hat r}\)、\(C_{n\delta_r}\) が変化する。

\(|C_{n\delta_r}|\) の増加は必要ラダー舵角を小さくする方向に働くが、方向安定とヨー減衰も同時に変わるため、必要舵角はラダー効きだけでは決まらない。

\(\Gamma_{\mathrm{eff}}\) を増やすと、主に \(C_{l\beta}\) の絶対値が増える。

\(C_{l\beta}\) は簡易式の分子と分母の両方に入るため、\(\delta_{r,\mathrm{trim}}\) の \(\Gamma_{\mathrm{eff}}\) 感度は単純な反比例にはならない。

したがって、必要ラダー舵角は、\(V_v\) による方向安定・ヨー減衰・ラダー効きの変化と、\(\Gamma_{\mathrm{eff}}\) による上反角効果の変化を合わせて判断する。

符号は旋回方向と舵角規約に依存するため、設計上は必要舵角の絶対値と、設定した最大ラダー舵角に収まるかを見る。

↓詳細な導出はこちら。

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

旋回性能

旋回性能では、最大ラダー舵角の範囲内で定常的につり合える最大バンク角を確認する。

数値計算指標は、ラダー舵角を正負の限界値に固定して定常トリムを解き、成立した解の最大絶対バンク角を選んだ

\begin{align} |\phi|_{\max,\delta_r} =\max \left( |\phi(-\delta_{r,\max})|, |\phi(+\delta_{r,\max})| \right) \end{align}

である。

片側でトリムが成立しない場合は、成立した側だけを候補とする。

今回、ラダー限界定常旋回が成立した 182 ケースでは、

\begin{align} 0.0241\ \mathrm{deg} \le |\phi|_{\max,\delta_r} \le 7.495\ \mathrm{deg} \end{align}

であった。

同じ状態における旋回率と横滑り角は、

\begin{align} -4.433\ \mathrm{deg/s} &\le \Omega_{\delta_r,\lim} \le 5.886\ \mathrm{deg/s}, \\ -19.00\ \mathrm{deg} &\le \beta_{\delta_r,\lim} \le 18.52\ \mathrm{deg} \end{align}

である。

ラダー舵角に対する定常バンク角のゲインは、

\begin{align} K_\phi \simeq -\frac{ C_{n\delta_r} }{ C_{l\beta}C_{n\hat r} -C_{n\beta}C_{l\hat r} } \left( \frac{4m}{\rho Sb}C_{l\beta} +C_{Y\beta}C_{l\hat r} \right) \end{align}

である。

このとき、指定ラダー舵角に対するバンク角は、

\begin{align} \phi \simeq \frac{q_\infty S}{mg} K_\phi\delta_r \end{align}

となる。

したがって、最大ラダー舵角に対する簡易最大バンク角を完全に書き下すと、

\begin{align} |\phi|_{\max,\mathrm{simple}} \simeq \frac{q_\infty S}{mg} \left| -\frac{ C_{n\delta_r} }{ C_{l\beta}C_{n\hat r} -C_{n\beta}C_{l\hat r} } \left( \frac{4m}{\rho Sb}C_{l\beta} +C_{Y\beta}C_{l\hat r} \right) \right| \delta_{r,\max} \end{align}

である。

今回、数値計算指標と簡易指標を同時に得られた 182 ケースの Pearson 相関係数は、

\begin{align} r=0.989 \end{align}

であり、設計平面上の傾向はよく対応した。

\(V_v\) を増やすと \(C_{n\delta_r}\) と \(C_{n\beta}\) の絶対値が増え、ラダー効きと方向安定が同時に変化する。

\(\Gamma_{\mathrm{eff}}\) を増やすと \(C_{l\beta}\) の絶対値が増え、横滑りからバンク角を作る経路が強くなる。

最大バンク角は、この両者とスパイラル安定性の組合せで決まる。

\(K_\phi\) はスパイラル判別式を分母に含むため、スパイラル中立付近では特異性を持つ。

この領域では定常トリム結果を優先する。

↓詳細はこちら。

航空機のラダーのみ旋回における旋回性能
ラダーのみを用いて旋回する航空機の旋回性能について説明する。

横転性能

横転性能では、ラダー入力後の有限時間内にどれだけバンク角を作れるかを確認する。

6 自由度有限時間応答から、平均バンク角速度をラダー舵角で正規化した指標を、

\begin{align} K_{\mathrm{6DOF,roll}} =\frac{\Delta\phi/\Delta t}{\delta_r} \end{align}

とする。

今回、応答計算が成立した 174 ケースでは、

\begin{align} -0.07277 \le K_{\mathrm{6DOF,roll}} \le -0.009356 \end{align}

であった。

簡易ラダー横転指数は、主要な二つのロール生成経路をまとめて、

\begin{align} K_{\mathrm{rudder,roll}} =\left( C_{l\hat r} -C_{l\beta} \right) C_{n\delta_r} \end{align}

である。

今回の範囲は、

\begin{align} -0.004677 \le K_{\mathrm{rudder,roll}} \le -0.000371 \end{align}

である。

両指標を同時に得られた 174 ケースの Pearson 相関係数は

\begin{align} 0.973 \end{align}

であり、設計平面上の傾向はよく対応した。

\(V_v\) を増やすと \(C_{n\delta_r}\) の絶対値が増え、ラダー入力からヨーレートを作る経路が強くなる。

\(\Gamma_{\mathrm{eff}}\) を増やすと \(C_{l\beta}\) の絶対値が増え、横滑りからロールへ戻る経路が強くなる。

したがって、横転性能を高めるには、両設計変数の組合せを確認する必要がある。

↓詳細な導出はこちら。

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

横風突風応答

横風突風応答では、1-cosine 横風突風を受けたときのバンク応答と、安定微係数だけから作る簡易指標を確認する。

数値計算指標には、突風通過終了時のバンク角変化をピーク突風横滑り角で正規化した、

\begin{align} K_{\phi,g} =\frac{\Delta\phi(\tau_e)}{\beta_{g,\max}} \end{align}

を用いる。

今回の範囲は、

\begin{align} -0.17672 \le K_{\phi,g} \le 0.05226 \end{align}

である。

簡易横風突風ロール指数は、横滑り角から直接ロールする経路と、方向安定からヨーレートを介してロールする経路をまとめて、

\begin{align} K_{\mathrm{gust,roll}} =\frac{8I_z}{\rho Sb^3}C_{l\beta} +C_{l\hat r}C_{n\beta} \end{align}

とする。

今回の範囲は、

\begin{align} -0.009390 \le K_{\mathrm{gust,roll}} \le 0.003321 \end{align}

である。

両指標の Pearson 相関係数は、

\begin{align} r=0.997 \end{align}

であり、今回の設計平面では非常によく対応した。

\(\Gamma_{\mathrm{eff}}\) を増やすと \(C_{l\beta}\) の絶対値が増え、横風から直接ロールする経路が強くなる。

\(V_v\) を増やすと \(C_{n\beta}\) と \(C_{n\hat r}\) が変化し、ヨー経由のロール応答とヨー減衰が同時に変わる。

↓詳細な導出はこちら。

航空機の横風突風応答
航空機の横風突風応答について説明する。

スパイラル安定性との対応

スパイラルモードの判別式を、

\begin{align} D_{\mathrm{sp}} =C_{l\beta}C_{n\hat r} -C_{n\beta}C_{l\hat r} \end{align}

とする。

今回のチャートの 228 ケースでは、\(D_{\mathrm{sp}}\) と突風終了時のバンク応答 \(K_{\phi,g}\) の Pearson 相関係数は、

\begin{align} r\left( D_{\mathrm{sp}}, K_{\phi,g} \right) &=-0.9953, \\ r\left( |D_{\mathrm{sp}}|, |K_{\phi,g}| \right) &=0.9761 \end{align}

であった。

符号付き相関が負になることは、今回の符号規約と、次の恒等変形に現れる \(-D_{\mathrm{sp}}\) に対応する。

\begin{align} K_{\mathrm{gust,roll}} =\left( \frac{8I_z}{\rho Sb^3} +C_{n\hat r} \right)C_{l\beta} -\left( C_{l\beta}C_{n\hat r} -C_{n\beta}C_{l\hat r} \right) \end{align}

今回の機体諸元では、

\begin{align} \frac{8I_z}{\rho Sb^3} =0.01843 \end{align}

である。

一方、\(C_{n\hat r}\) は、

\begin{align} -0.01163 \le C_{n\hat r} \le -0.00714 \end{align}

で、中央値は、

\begin{align} \operatorname{median} \left( C_{n\hat r} \right) =-0.00975 \end{align}

である。

したがって、

\begin{align} \frac{8I_z}{\rho Sb^3} \simeq 1.89 \left| \operatorname{median} \left( C_{n\hat r} \right) \right| \end{align}

となり、両者は同じオーダーである。

さらに、

\begin{align} 0.00680 \le \frac{8I_z}{\rho Sb^3} +C_{n\hat r} \le 0.01129 \end{align}

で、中央値は約 0.00869 となる。

正の慣性項と負の \(C_{n\hat r}\) が部分的に相殺し、\(C_{l\beta}\) に掛かる残差係数が小さくなったため、今回のチャートでは \(-D_{\mathrm{sp}}\) の分布が目立ち、横風突風応答とスパイラル安定性の対応が強くなった。

機体諸元や飛行条件が変われば、この対応の強さも変わる。

抵抗性能

抵抗性能では、基準状態の抵抗係数と、ラダーのみ定常旋回における横滑り角を見る。

今回の基準抵抗係数は、

\begin{align} 0.021355 \le C_{D,\mathrm{base}} \le 0.021748 \end{align}

であり、今回の \(V_v-\Gamma_{\mathrm{eff}}\) 範囲では変化幅は小さい。

ただし、VSPAERO の基準抵抗係数だけでは、実機の表面摩擦、支持部、隙間、ラダー偏向に伴う抵抗増加、構造重量増加をすべて評価できない。

垂直尾翼を大きくした設計では、別途、寄生抵抗と重量を確認する必要がある。

定常旋回中の横滑り角は、横力と誘導抵抗、ラダー偏向に伴う抵抗増加の目安になる。

今回のラダー限界定常旋回では、

\begin{align} -19.00\ \mathrm{deg} \le \beta_{\delta_r,\lim} \le 18.52\ \mathrm{deg} \end{align}

であった。

最大バンク角だけが大きくても、大きな横滑り角を必要とする設計は、抵抗と飛行姿勢の観点から採用しにくい。

したがって、\(|\phi|_{\max,\delta_r}\) と \(|\beta_{\delta_r,\lim}|\) を同時に確認する。

垂直尾翼と上反角設計のまとめ

ここまでの結果をまとめると、ラダーのみで旋回する航空機の垂直尾翼と上反角設計は、次の 5 つの性能のバランスで決まる。

  • スパイラル安定性
  • 旋回性能
  • 横転性能
  • 横風突風応答
  • 抵抗性能

設計変数ごとの効き方をまとめると、次のようになる。

設計指標としては、まず \(V_v\) と \(\Gamma_{\mathrm{eff}}\) の組合せで、スパイラルモードが中立からやや安定側に入る領域を探すのがよい。

上反角の大きさは、横転性能と横風突風応答のトレードオフで決める。

横転性能を高めるには大きな上反角が有利だが、横風突風に対する耐性を高めるには小さな上反角が有利である。

したがって、上反角は、垂直尾翼とのバランスでスパイラルモードを中立からやや安定側に置きつつ、横転性能が必要量を満たす範囲に設定する。

垂直尾翼の大きさは、操縦性と抵抗特性のトレードオフで決める。

ラダーのみ旋回を成立させるには、垂直尾翼をある程度大きくして、方向安定、ラダー効き、ヨー減衰を確保する必要がある。

一方で、過大な垂直尾翼は抵抗と重量を増やす。

したがって、フライトプランをよく練り、旋回性能、横転性能、横風突風応答、抵抗性能のどれを優先するかを明確にした上で、\(V_v\) を決める。

このように、\(V_v-\Gamma_{\mathrm{eff}}\) チャートは、ラダーのみで旋回する航空機に必要な操縦性、安定性、突風耐性を同じ平面上で比較する設計ツールとして活用できる。

おわりに

この記事では、ラダーのみを用いて旋回する航空機について、垂直尾翼容積比 \(V_v\) と等価上反角 \(\Gamma_{\mathrm{eff}}\) を軸にした設計チャートの作り方を整理した。

\(V_v-\Gamma_{\mathrm{eff}}\) チャートを使うと、旋回性能、横転性能、突風応答性能、スパイラル安定性を同じ平面上で比較できる。









Appendix

ここでは、翼端たわみ量 \(w_{\mathrm{tip}}\) から代表的なたわみ形状を作るため、楕円揚力分布を受ける桁径一定の片持ち翼を考える。

楕円揚力分布を受ける片持ち翼

半スパンを \(s\)、無次元スパン座標を

\begin{align} x=\frac{y}{s} \end{align}

とする。

単位スパンあたり揚力を次でおく。

\begin{align} \ell(y)=\ell_0\sqrt{1-\left(\frac{y}{s}\right)^2} \end{align}

全翼揚力を \(L_w\) とすると、

\begin{align} L_w &=2\int_0^s\ell(y)\,dy \\ &=2\ell_0s\int_0^1\sqrt{1-x^2}\,dx \\ &=2\ell_0s\frac{\pi}{4} \\ &=\frac{\pi}{2}\ell_0s \end{align}

したがって、

\begin{align} \ell_0=\frac{2L_w}{\pi s} \end{align}

である。

片持ち翼の断面 \(y\) における曲げモーメントは、翼端側の分布荷重によって

\begin{align} M(y)=\int_y^s\ell(\eta)(\eta-y)\,d\eta \end{align}

である。

無次元座標で書くと、

\begin{align} M(x)=\ell_0s^2\int_x^1\sqrt{1-\xi^2}(\xi-x)\,d\xi \end{align}

である。

積分を実行すると、

\begin{align} M(x)=\ell_0s^2\left[ \frac{x^2\sqrt{1-x^2}}{6} +\frac{x\sin^{-1}x}{2} -\frac{\pi x}{4} +\frac{\sqrt{1-x^2}}{3} \right] \end{align}

となる。

桁径一定の場合のたわみ角分布

桁径一定、すなわち断面二次モーメント \(I\) がスパン方向に一定の場合を考える。

Euler-Bernoulli はりとして、

\begin{align} EI\frac{d^2w}{dy^2}=M(y) \end{align}

を用いる。

たわみ角は

\begin{align} \theta_b(y)=\frac{dw}{dy}=\int_0^y\frac{M(\eta)}{EI}\,d\eta \end{align}

である。

無次元座標で書くと、

\begin{align} \theta_b(x) =\frac{\ell_0s^3}{EI} \int_0^x \left[ \frac{\xi^2\sqrt{1-\xi^2}}{6} +\frac{\xi\sin^{-1}\xi}{2} -\frac{\pi\xi}{4} +\frac{\sqrt{1-\xi^2}}{3} \right]d\xi \end{align}

である。

この積分を実行すると、

\begin{align} \theta_b(x) =\frac{\ell_0s^3}{EI} \left[ \frac{4x^2+1}{16}\sin^{-1}x -\frac{\pi x^2}{8} +\frac{x\sqrt{1-x^2}(13+2x^2)}{48} \right] \end{align}

である。

\(\ell_0=2L_w/(\pi s)\) を代入すると、

\begin{align} \theta_b(x) =\frac{2L_ws^2}{\pi EI} \left[ \frac{4x^2+1}{16}\sin^{-1}x -\frac{\pi x^2}{8} +\frac{x\sqrt{1-x^2}(13+2x^2)}{48} \right] \end{align}

である。

翼端では \(x=1\) であるから、

\begin{align} \theta_{b,\mathrm{tip}} &=\frac{\ell_0s^3}{EI}\left(\frac{5\pi}{32}-\frac{4\pi}{32}\right) \\ &=\frac{\ell_0s^3}{EI}\frac{\pi}{32} \end{align}

である。

\(\ell_0=2L_w/(\pi s)\) を代入すると、

\begin{align} \theta_{b,\mathrm{tip}}=\frac{L_ws^2}{16EI} \end{align}

である。

さらに \(s=b/2\) を代入すれば、

\begin{align} \theta_{b,\mathrm{tip}}=\frac{L_wb^2}{64EI} \end{align}

となる。

桁径一定の場合のたわみ分布

たわみは、

\begin{align} w(y)=\int_0^y(y-\eta)\frac{M(\eta)}{EI}\,d\eta \end{align}

である。

無次元座標で書くと、

\begin{align} w(x) =\frac{\ell_0s^4}{EI} \int_0^x (x-\xi) \left[ \frac{\xi^2\sqrt{1-\xi^2}}{6} +\frac{\xi\sin^{-1}\xi}{2} -\frac{\pi\xi}{4} +\frac{\sqrt{1-\xi^2}}{3} \right]d\xi \end{align}

である。

この積分を実行すると、

\begin{align} w(x) =\frac{\ell_0s^4}{EI} \left[ \frac{x(4x^2+3)}{48}\sin^{-1}x -\frac{\pi x^3}{24} +\frac{\sqrt{1-x^2}(6x^4+83x^2+16)}{720} -\frac{1}{45} \right] \end{align}

である。

\(\ell_0=2L_w/(\pi s)\) を代入すると、

\begin{align} w(x) =\frac{2L_ws^3}{\pi EI} \left[ \frac{x(4x^2+3)}{48}\sin^{-1}x -\frac{\pi x^3}{24} +\frac{\sqrt{1-x^2}(6x^4+83x^2+16)}{720} -\frac{1}{45} \right] \end{align}

である。

翼端たわみは \(x=1\) を代入して、

\begin{align} w_{\mathrm{tip}} =\frac{\ell_0s^4}{EI} \left(\frac{\pi}{32}-\frac{1}{45}\right) \end{align}

である。

\(\ell_0=2L_w/(\pi s)\) を代入すると、

\begin{align} w_{\mathrm{tip}} =\frac{2L_ws^3}{\pi EI} \left(\frac{\pi}{32}-\frac{1}{45}\right) \end{align}

である。

さらに \(s=b/2\) を代入すれば、

\begin{align} w_{\mathrm{tip}} =\frac{L_wb^3}{4\pi EI} \left(\frac{\pi}{32}-\frac{1}{45}\right) \end{align}

である。

翼端たわみ量からたわみ角分布を作る

チャート作成では、構造荷重や剛性の詳細指定よりも、翼端たわみ量 \(w_{\mathrm{tip}}\) を直接指定する方が扱いやすいため、前節までに含まれていた \(L_w/EI\) を、指定した翼端たわみ量 \(w_{\mathrm{tip}}\) で消去する。

翼端たわみの完全書き下し式より、

\begin{align} w_{\mathrm{tip}} =\frac{2L_ws^3}{\pi EI} \left(\frac{\pi}{32}-\frac{1}{45}\right) \end{align}

である。

これを \(L_w/EI\) について解くと、

\begin{align} \frac{L_w}{EI} =\frac{\pi w_{\mathrm{tip}}} {2s^3\left(\frac{\pi}{32}-\frac{1}{45}\right)} \end{align}

である。

この関係をたわみ角分布の \(L_w/EI\) に代入すると、

\begin{align} \theta_b(x) =\frac{w_{\mathrm{tip}}} {s\left(\frac{\pi}{32}-\frac{1}{45}\right)} \left[ \frac{4x^2+1}{16}\sin^{-1}x -\frac{\pi x^2}{8} +\frac{x\sqrt{1-x^2}(13+2x^2)}{48} \right] \end{align}

である。

さらに \(s=b/2\) を代入すれば、

\begin{align} \theta_b(x) =\frac{2w_{\mathrm{tip}}} {b\left(\frac{\pi}{32}-\frac{1}{45}\right)} \left[ \frac{4x^2+1}{16}\sin^{-1}x -\frac{\pi x^2}{8} +\frac{x\sqrt{1-x^2}(13+2x^2)}{48} \right] \end{align}

となる。

これにより、翼端たわみ量 \(w_{\mathrm{tip}}\) を直接入力として、弾性上反角分布を作ることができる。

VvGammaChart.py では、この考え方に対応して、elastic_wing_deflection_distribution() が翼端たわみ量からたわみ角分布を作り、apply_wing_deflection_as_dihedral() がそのたわみ角を OpenVSP の WingGeom の Dihedral へ近似的に加える。

楕円コード長分布・桁径一定を仮定したときの等価上反角

ここでは、前節で定義した重み付き平均に、楕円コード長分布と桁径一定のたわみ角分布を代入する。

半スパンを \(s\)、無次元スパン座標を

\begin{align} x=\frac{y}{s} \end{align}

とする。

楕円コード長分布を

\begin{align} c(x)=c_0\sqrt{1-x^2} \end{align}

とおく。

等価上反角の重みは、ロールモーメントアーム \(y\) と局所空力寄与 \(c(y)a(y)\) の積である。

局所揚力傾斜 \(a(y)\) を一定とすると、重みは

\begin{align} yc(y)a(y)\,dy \propto (sx)c_0\sqrt{1-x^2}\,s\,dx \propto x\sqrt{1-x^2}\,dx \end{align}

となる。

したがって、楕円コード長分布・局所揚力傾斜一定のもとで、等価上反角は

\begin{align} \Gamma_{\mathrm{eff}} =\frac{ \displaystyle\int_0^1x\sqrt{1-x^2}\,\theta_b(x)\,dx }{ \displaystyle\int_0^1x\sqrt{1-x^2}\,dx } \end{align}

である。

分母は

\begin{align} \int_0^1x\sqrt{1-x^2}\,dx=\frac{1}{3} \end{align}

である。

一方、桁径一定翼のたわみ角分布は、前節で導出したように

\begin{align} \theta_b(x) =\frac{2w_{\mathrm{tip}}} {b\left(\frac{\pi}{32}-\frac{1}{45}\right)} \left[ \frac{4x^2+1}{16}\sin^{-1}x -\frac{\pi x^2}{8} +\frac{x\sqrt{1-x^2}(13+2x^2)}{48} \right] \end{align}

である。

これを等価上反角の定義に代入すると、

\begin{align} \Gamma_{\mathrm{eff}} =\frac{ \displaystyle \int_0^1 x\sqrt{1-x^2} \frac{2w_{\mathrm{tip}}} {b\left(\frac{\pi}{32}-\frac{1}{45}\right)} \left[ \frac{4x^2+1}{16}\sin^{-1}x -\frac{\pi x^2}{8} +\frac{x\sqrt{1-x^2}(13+2x^2)}{48} \right]dx }{ \displaystyle\int_0^1x\sqrt{1-x^2}\,dx } \end{align}

である。

さらに分母の積分値を代入すると、

\begin{align} \Gamma_{\mathrm{eff}} =3\int_0^1 x\sqrt{1-x^2} \frac{2w_{\mathrm{tip}}} {b\left(\frac{\pi}{32}-\frac{1}{45}\right)} \left[ \frac{4x^2+1}{16}\sin^{-1}x -\frac{\pi x^2}{8} +\frac{x\sqrt{1-x^2}(13+2x^2)}{48} \right]dx \end{align}

となる。

この積分を実行すると、

\begin{align} \int_0^1 x\sqrt{1-x^2} \left[ \frac{4x^2+1}{16}\sin^{-1}x -\frac{\pi x^2}{8} +\frac{x\sqrt{1-x^2}(13+2x^2)}{48} \right]dx =\frac{128}{1575}-\frac{\pi}{60} \end{align}

である。

したがって、完全に整理すると、

\begin{align} \Gamma_{\mathrm{eff}} &= \frac{ 6\left(\frac{128}{1575}-\frac{\pi}{60}\right) }{ \left(\frac{\pi}{32}-\frac{1}{45}\right) } \frac{w_{\mathrm{tip}}}{b} \\ &= \frac{48(512-105\pi)}{35(45\pi-32)} \frac{w_{\mathrm{tip}}}{b} \end{align}

である。

数値的には

\begin{align} \Gamma_{\mathrm{eff}}\simeq2.284\frac{w_{\mathrm{tip}}}{b} \end{align}

である。

VvGammaChart.py の calculate_gamma_eff_from_tip_deflection() は、この定義に対応している。

まず elastic_wing_deflection_distribution() により、翼端たわみ量 \(w_{\mathrm{tip}}\) と半スパンから、楕円揚力分布・桁径一定モデルのたわみ角分布 \(\theta_b(x)\) を作る。

次に、

\begin{align} \Gamma_{\mathrm{eff}} =\frac{ \displaystyle\int_0^1x\sqrt{1-x^2}\,\theta_b(x)\,dx }{ \displaystyle\int_0^1x\sqrt{1-x^2}\,dx } \end{align}

を数値積分で評価する。

質問・感想・意見などあれば気軽にTwitterのDMかコメントお願いします!
スポンサーリンク