全機の計算

今まで作ったプログラムを用いて全機の計算を行う

目的とフローチャート
プログラムの解説
プログラムの実行
おまけ　動安定解析
まとめ

目的とフローチャート

今回の目的は，今まで作ったプログラムを用いて，全機の空力計算を行うことである．

今まで作ったプログラムについては次の記事を参考にしてほしい

≫設計シートに使う構造体について
≫主翼の空力および構造計算
≫渦格子法を用いた尾翼の空力計算

また，全機計算ではカウルにはたらく空気力も考慮する

カウルの空気力の計算については次の文献を参考にした

≫飛行機設計論　山名正夫著・・・(文献1)

その他の空力係数の計算については次の文献を参考にした

≫航空機力学入門　加藤寛一郎著・・・(文献2)
≫航空機の飛行力学と制御　片柳亮二著・・・(文献３)

フローチャート

フローチャートを以下に示す

ちなみにコードは次のようになる

Sub 全機計算()'全機の計算を行う
'シートの宣言
Set sht1 = Sheets("主翼(sht1)")
Set sht3 = Sheets("水平尾翼(sht3)")
Set sht4 = Sheets("垂直尾翼(sht4)")
Set sht8 = Sheets("全機計算(sht8)")
'構造体の宣言
Dim state As variables
Dim Plane As Specifications

state = ReadState(state, 8) '状態変数の読み込み

'それぞれのシートに状態変数の書き出し
sht1.Range("BB20") = state.alpha: sht3.Range("BB20") = state.alpha: sht4.Range("BB20") = state.alpha
sht1.Range("BB21") = state.beta: sht3.Range("BB21") = state.beta: sht4.Range("BB21") = state.beta
sht1.Range("BB22") = state.Vair: sht3.Range("BB22") = state.Vair: sht4.Range("BB22") = state.Vair
sht1.Range("BB23") = state.hE: sht3.Range("BB23") = state.hE: sht4.Range("BB23") = state.hE
sht1.Range("BB27") = state.dh: sht3.Range("BB27") = state.dh: sht4.Range("BB27") = state.dh
sht1.Range("BB28") = state.de: sht3.Range("BB28") = state.de: sht4.Range("BB28") = state.de
sht1.Range("BB29") = state.dr: sht3.Range("BB29") = state.dr: sht4.Range("BB29") = state.dr
'全機計算
Plane = Plane_calculation(1, Plane, state)
End Sub

Public Function Plane_calculation(ByVal p As Integer, ByRef Plane As Specifications, ByRef state As variables) As Specifications
'ワークシートの定義
Set sht8 = Sheets("全機計算(sht8)")
'構造体の宣言
Dim Wing As Specifications
Dim Tail As Specifications
Dim Fin As Specifications
Dim Cowl As Specifications

'状態変数の読み込み
Wing = Read_value(Wing, state, 1)
Tail = Read_value(Tail, state, 3)
Fin = Read_value(Fin, state, 4)
Wing = VLM_wing(p, Wing, state) '主翼計算
Plane = Wing 'PlaneにWingを代入
If sht8.Range("N19") = "Conventional" Then
    Tail.epsilon = 2 * Wing.epsilon '尾翼位置における吹きおろし角 [deg]
Else 'T字尾翼
    Tail.epsilon = 0
End If
Tail = VLM_tail(p, Tail, state) '水平尾翼の計算
Fin = VLM_Fin(p, Fin, state) '垂直尾翼の計算
Cowl = Cowl_Calculation(p, Plane, state, Cowl) 'カウルの計算
Call calculation_aero_coef(state, Plane, Wing, Tail, Fin, Cowl) '空力係数の計算
Plane_calculation = Plane '値を返す

End Function

全機計算のプログラム自体はそんなに難しくない

今までのプログラムを組み合わせるだけですむ

プログラムの解説

それでは実際に解説していく

前半部分は既に解説したものばかりなので，後半のカウルの計算と空力計算のみを解説する

シートの下準備

全機計算のシートに，それぞれのシートから値を持ってくる

カウル計算

カウルの抗力とピッチングモーメントを計算する

とはいいつつ正確な値を知るにはCFDで計算するか風洞試験を行うしかないので，ここで求めるのはあくまで概算値である

まず，計算に必要な次の値をシートに入力する

最小有害抗力係数
全長 [m]
取り付け角 [deg]　(カウルのコードラインと機体軸のなす角)
体積 [m^3]
前方投影面積 [m^2]

カウルの最小有害抗力係数を知るすべはないので，地面効果を含めない全機のL/Dが40になるように決定した

あとは二次元翼と同様に迎角とレイノルズ数で抗力を計算する

カウルのピッチングモーメントは，文献1のp.284を参考に，軸対称胴体が迎角\(\alpha\)，速度\(V\)の一様な非圧縮の流れの中に置かれたときのポテンシャル理論による次の理論式によって計算する

\begin{eqnarray}
C_{m}=2\alpha(\frac{\Delta}{S_{\pi}l}-\frac{S_{b}}{S_{\pi}})
\end{eqnarray}

ここで，\(S_{\pi}\)はカウルの前方投影面積，\(S_{b}\)はカウル後端の断面積，\(l\)はカウルの全長である

カウル後端がとがっていれば，\(S_{b}=0\)なので，カウルの\(C_{m}\)は次のようになる

\begin{eqnarray}
C_{m}=2\alpha\frac{\Delta}{S_{\pi}l}
\end{eqnarray}

今回のプログラムではこの式を使っている

また，胴体を含めた抗力係数は次の式で計算する

\begin{eqnarray}
C_{D_{cowl-fulelage}}=C_{interference}(C_{D_{cowl}}+\frac{S_{fuselage}}{S_{cowl}}C_{D_{fuselage}})
\end{eqnarray}

干渉抗力係数\(C_{interference}\)は適当に1.2とし，機体が左右対称なら\(S_{fuselage}=0\)とする

コードは次のようになる

Function Cowl_Calculation(ByVal p As Integer, ByRef Plane As Specifications, ByRef state As variables, ByRef Cowl As Specifications) As Specifications
'ワークシートの定義
Set sht14 = Sheets("カウル側面翼型解析(sht14)")
Set sht8 = Sheets("全機計算(sht8)")
Dim Fuselage As Specifications '構造体の宣言
'変数
Dim CD_interference As Double: CD_interference = sht8.Range("D28")
Dim CD0 As Double '迎角0のときの翼の抗力係数
Dim Cdp_coef(14) As Double '有害抗力係数を表す多項式の係数
'値の読み込み
Cowl.S = sht8.Range("N33")
Cowl.CD = sht8.Range("I34")
Cowl.Cm_cg = sht8.Range("I35")
Cowl.chord_mac = sht8.Range("I36")
Fuselage.S = sht8.Range("N34")
Fuselage.CD = sht8.Range("D27")
'カウルの抗力を計算
For i = 0 To 14
    Cdp_coef(i) = sht14.Cells(4, 7 + i)
Next i
alpha = state.alpha
Re = (state.Vair * Cowl.chord_mac) / state.rho
CD0 = Cdp_coef(0) '迎角0度のときの翼の抗力係数を計算
For i = 9 To 14
    CD0 = CD0 + Cdp_coef(i) * (Re ^ (i - 8))
Next i
CD0 = Cowl.CD / CD0 '翼の抗力係数をカウルの抗力係数と同じスケールに変換する係数
Cowl.CD = 0
For i = 0 To 8
    Cowl.CD = Cowl.CD + CD0 * Cdp_coef(i) * (alpha ^ i)
Next i
For i = 9 To 14
    Cowl.CD = Cowl.CD + CD0 * Cdp_coef(i) * (Re ^ (i - 8))
Next i
If p = 1 Then
    sht8.Range("D26") = Cowl.CD '結果の表示
End If
Cowl.CD = CD_interference * (Cowl.CD + (Fuselage.S / Cowl.S) * Fuselage.CD)
Cowl.Drag = Plane.dynamic_pressure * Cowl.S * Cowl.CD
Cowl.M_pitch = Plane.dynamic_pressure * Cowl.S * Cowl.Cm_cg * Cowl.chord_mac
Cowl_Calculation = Cowl
End Function

このプログラムについては以上である

空力係数の計算

機体の空力係数を計算する

~~いちいち解説するのが面倒くさい~~式は参考文献に書いてあるものをそのまま書き写すだけなので，各自でコードを確認してほしい

ただし，ピッチングモーメントの計算のみ，カウルのモーメントと主翼の抗力によるモーメントを加えている

\begin{eqnarray}
M&=&M_{wing}-D_{wing}\times Z_{cg}+M_{tail}+M_{cowl} \\
C_{m}&=&C_{m_{wing}}-\frac{Z_{cg}}{c_{mac}}C_{D_{wing}}+C_{m_{tail}}+C_{m_{cowl}}
\end{eqnarray}

ここで，\(Z_{cg}\)は機体重心と主翼-胴体接合部の距離，\(M_{tail}=L_{tail}\times l_{t}\)であり，低翼機なので主翼抗力によるモーメントは負としている

コードは次のようになる

Function calculation_aero_coef(ByRef state As variables, ByRef Plane As Specifications, ByRef Wing As Specifications, _
                               ByRef Tail As Specifications, ByRef Fin As Specifications, ByRef Cowl As Specifications)
Set sht8 = Sheets("全機計算(sht8)") 'ワークシートの定義
'変数
Dim eta As Double: eta = sht8.Range("N17") '尾翼効率
Dim Z_cg As Double: Z_cg = sht8.Range("I28") '機体重心のZ方向位置
Dim pi As Double: pi = 4 * Atn(1) '円周率

'全機の空力係数の計算
Plane.VH = (Tail.S * Tail.lt) / (Wing.S * Wing.chord_mac) '水平尾翼容積比
With Plane
    .CL = Wing.CL + (Tail.S / Wing.S) * Tail.CL '揚力係数
    .CD = Wing.CD + (Tail.S / Wing.S) * Tail.CD + (Fin.S / Wing.S) * Fin.CD + (Cowl.S / Wing.S) * Cowl.CD '抗力係数
    .CDi = Wing.CDi + (Tail.S / Wing.S) * Tail.CDi + (Fin.S / Wing.S) * Fin.CDi '誘導抗力係数
    .Cdp = Wing.Cdp + (Tail.S / Wing.S) * Tail.Cdp + (Fin.S / Wing.S) * Fin.Cdp + (Cowl.S / Wing.S) * Cowl.CD '有害抗力係数
    .Cm_cg = Wing.Cm_cg - Wing.CD * (Z_cg / Wing.chord_mac) - eta * Plane.VH * Tail.CL + ((Cowl.S * Cowl.chord_mac) / (Wing.chord_mac * Wing.S)) * Cowl.Cm_cg '重心まわりのピッチングモーメント係数
    .de_da = (2 * Wing.aw * (180 / pi)) / (pi * Wing.AR) '∂ε/∂α
    .aw = Wing.aw + eta * (Tail.S / Wing.S) * (1 - de_da) * Tail.aw '全機揚力傾斜
    .L_D = .CL / .CD '揚抗比
    .Lift = Wing.Lift + Tail.Lift '揚力
    .Drag = Wing.Drag + Tail.Drag + Fin.Drag + Cowl.Drag '抗力
    .Drag_induced = Wing.Drag_induced + Tail.Drag_induced + Tail.Drag_induced + Fin.Drag_induced '誘導抗力
    .Drag_parasite = Wing.Drag_parasite + Tail.Drag_parasite + Fin.Drag_parasite + Cowl.Drag '有害抗力
    .L_roll = Wing.L_roll + Tail.L_roll + Fin.L_roll 'ローリングモーメント
    .M_pitch = Wing.M_pitch - Wing.Drag * Z_cg + Tail.M_pitch + Cowl.M_pitch 'ピッチングモーメント
    .N_yaw = Wing.N_yaw + Fin.N_yaw 'ヨーイングモーメント
End With
End Function

構造体の便利さが遺憾なく発揮されたコードだといえる

このプログラムについては以上である