主翼の空力および構造計算

主翼の空力および構造計算を行うプログラムを作成する

目的

今回の目的は,主翼の空力計算を行うことである.

XFLR5で行う空力計算とは異なり,たわみやねじれも同時に計算する

循環分布の計算についてはTR-797と同じ方法をとるので,詳しくは次の記事(≫TR-797の実装)を参考にしてほしい

構造設計については次の本を参考にした

グライダーの製作(Ⅰ) 『木製ボックス桁』・・・(文献1)

また,非線形の空力計算については次の文献を参考にした

Non-Linear Lifting Line Implementation and Validation for Aerodynamic and Stability Analysis・・・(文献2)

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

航空機力学入門|加藤寛一郎・・・(文献3)

ちなみに本家XFLR5の空力計算は次の文献を参考にしている

Method for calculating wing characteristics by lifting-line theory using nonlinear section lift data ・・・(文献4)

今回の記事では,文献1,文献2,文献3をもとに空気力を計算→翼のたわみとねじれを計算→空気力を計算→翼のたわみとねじれを計算→・・・を収束するまで行っていくことになる

また,参考文献を参照する際に助けとなるよう,式番号などはそれぞれの文献のものにそろえる

フローチャート

今回のプログラムのフローチャートを以下に示す

長い.それに尽きる

プログラムの解説

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

テンプレートファイル

エクセルシートの配置などの参考になるようにテンプレートファイルを添付しておく

※数式やマクロはすべて削除してあり,値も適当である

プログラムの構造

主翼の空力計算のプログラムは,引数を主翼の構造体Wing,状態変数の構造体state,設計シートへの出力の有無pとして,主翼の構造体に計算結果を返すFunctionとする

設計シートのプログラム全般で使う構造体については次の記事(≫設計シートに使う構造体について)を参考にしてほしい

プログラムの全体の流れ

今回のプログラムはかなり長いので,可読性を上げるためにもそれぞれの処理はプロシージャに分けて行った

プログラムの全体は次のようになる

引数が多すぎて逆に見づらくなっている気がしないでもない

変数の宣言

今回のプログラムで使う変数を以下に示す

行数の節約のために動的配列の大きさの宣言や値の計算は同じ行で行っている

Re_max,Re_min,alpha_max,apha_minは,繰り返し計算の途中で迎角やRe数がXFLR5であらかじめ計算した範囲を出てしまったときに起こる発散を防ぐための値である(≫XFLR5で二次元翼を解析して結果を出力する

値の読み込み

値の読み込みを行う前に,座標系の定義と番号の振り方を説明しておく

座標系については,機体後方をx軸正,右翼側をy軸正,機体上方をz軸正とする

翼を有限このパネルに分割する計算法において,空気力や循環が働くとするパネルの中心をcontrol point (cp),パネルの境界をsurface point (sp)と呼ぶ

リブ枚数をn(翼中心のリブ番号が0,翼端のリブ番号がn)とすると,cp,spの番号は次のようになる

わかりにくいときは下のように小さなnで具体的に考えてみるといい

この番号をしっかり意識しながらシートから必要な値を読み込んでいく

設計シート上で元々リブ位置(sp)で定義してある量は,cpの両側のspの値を平均している

このプロシージャについては以上である

繰り返し計算のループ

今回のプログラムでは循環分布を計算するために繰り返し計算を行う

ループ内の最後の2行はエクセルシートの下のステータスバーに進捗状況を表示するコードと,Excel VBAあるあるの”応答なし”を防ぐためのコードである

y,zの計算

まず,Qijを計算するのに必要な値を順次計算していく

翼のたわみやねじれをもとに翼の各spの座標や取り付け角を更新する

翼の上反角はTR-797の実装のときと同じように,左翼では負の値になっていることに注意する

このプロシージャについては以上である

cpの計算

cpの座標を計算する

ただ平均をとるだけである

このプロシージャについては以上である

yp,zpの計算

yp,zpを計算する

基本的にはTR-797の実装のときと同じだが,今回のプログラムでは右翼と左翼のたわみなどは別々に計算するので,左翼のyp,zpを計算するときには左翼のy,zの値を使う

このプロシージャについては以上である

Rijの計算,Qijの計算

Rij,Qijを計算する

TR-797の実装のときと同じである

Qijが計算できたので,循環分布がわかれば吹きおろし分布を計算できる

このプロシージャについては以上である

循環の更新

文献2の式(8)に従って循環分を更新する

\begin{eqnarray}
\Gamma_{input}=\Gamma_{old}+D(\Gamma_{new}-\Gamma_{old})
\tag{8}
\end{eqnarray}

ここで\(\Gamma_{input}\)は吹きおろしを計算するための循環分布,\(Gamma_{new}\)は新しく計算した循環分布,\(\Gamma_{old}\)は1つ前の循環分布,\(D\)はDampin factorで発散を防ぐための係数である

繰り返し計算の最初(iteration=0)のときはまだ循環分布を計算していないためcirculationには何の値も入っていない

iteration=1のときはじめてcirculationに1回目の循環分布の値が入っており,circulation_oldにその値を格納できる

よって式(8)の方法で循環分布を更新できるのはcirculation_oldに1回目の循環分布が格納されているiteration=2のときであり,そのため場合分けを行っている

このプロシージャについては以上である

吹きおろしの計算

吹きおろし,誘導迎角,平均誘導迎角の計算を行う

このプロシージャについては以上である

有効迎角の計算

ここからは循環分布を計算するために必要な値を順次計算していく

まず,有効迎角の計算を行う

特に意味はないが機体がロールしたときや横滑りを起こしたときの解析もできるようにしておく

各cpにおける有効迎角は次のように計算できる

\begin{eqnarray}
(有効迎角) &= & (全機迎角)+(取り付け角)-(誘導迎角)+(ロールによる影響)+(横滑りによる影響)
\end{eqnarray}

文献3によると,対気速度\(V\),横滑り角\(\beta\)で運動しているとき,上反角\(\Gamma_{dihedral_{i}}\)のcpに生じている吹き上げは次の式であらわされる

\begin{eqnarray}
V\sin{\beta}\sin({\Gamma_{dihedral_{i}}})
\tag{4.64′}
\end{eqnarray}

同じく,ロール角速度\(p\)で運動しているとき,スパン方向位置\(y_{i}\)にあるcpに生じている吹き上げは次の式であらわされる

\begin{eqnarray}
p y_{i}
\tag{4.74′}
\end{eqnarray}

また,ヨー角速度\(r\)で運動しているとき,スパン方向位置\(y_{i}\)にあるcpにおける対気速度は次の式であらわされる

\begin{eqnarray}
V-r y_{i}
\tag{4.92′}
\end{eqnarray}

よって,i番目のcpの有効迎角は次の式で計算できる

\begin{eqnarray}
\alpha_{effective}=\alpha+\alpha_{setting}-\alpha_{induced}
+\frac{V\sin{\beta}\sin({\Gamma_{dihedral_{i}}})}{V-r y_{i}}
+\frac{p y_{i}}{V-r y_{i}}
\end{eqnarray}

すべての有効迎角を計算した後に,XFLR5であらかじめ計算した範囲をはみ出た迎角にたいして処理を行っている

このプロシージャについては以上である

Re数の計算

Re数の計算を行う

このプロシージャについては以上である

空力係数の計算を行う

局所揚力係数をはじめとした空力係数,および局所揚力などの空気力を計算する

すでに翼型解析②:解析データを最小二乗法で多項式近似するで,二次元翼の空力係数を迎角,Re数の任意の次数の多項式であらわしたときの係数を求めてあるので,局所揚力係数\(C_{l}\),局所抗力係数\(C_{d_{p}}\),局所ピッチングモーメント係数\(C_{m_{ac}}\)は容易に計算できる

空力中心周り\(h_{ac}\)のピッチングモーメント係数\(C_{m_{ac}}\)と桁位置\(h_{spar}\)周りのピッチングモーメント係数\(C_{m_{cg}}\)の関係は次の式であらわされる

\begin{eqnarray}
Cm_{cg}=Cm_{ac}+C_{l}(h_{spar}-h_{ac})
\end{eqnarray}

さらにそこから各cpにはたらく空気力も次のように計算できる

\begin{eqnarray}
\Delta L_{i} &=& \frac{1}{2}\rho V^2 \Delta S_{i} C_{l_{i}} \\
\Delta D_{p_{i}} &=& \frac{1}{2}\rho V^2 \Delta S_{i} C_{d_{p_{i}}} \\
\Delta M_{ac_{i}} &=& \frac{1}{2}\rho V^2 \Delta S_{i} c_{i} C_{m_{ac_{i}}} \\
\Delta M_{cg_{i}} &=& \frac{1}{2}\rho V^2 \Delta S_{i} c_{i} C_{m_{cg_{i}}} \\
\Delta S_{i} &=& c_{cp_{i}} dy\cos({\Gamma_{dihedral_{i}}})
\end{eqnarray}

また,文献1のp.40にならって,空気力を構造軸(n-t軸)へ投影する

図より,空気力のN軸方向成分,T軸方向成分はそれぞれ次の式であらわされる

\begin{eqnarray}
N &=& L\cos{(\alpha_{effective}-\alpha_{setting})}+D_{p}\sin{(\alpha_{effective}-\alpha_{setting})} \\
T &=& -L\sin{(\alpha_{effective}-\alpha_{setting})}+D_{p}\cos{(\alpha_{effective}-\alpha_{setting})} \\
\tag{文献1,p.40}
\end{eqnarray}

このプロシージャについては以上である

空気力を計算する

循環分布,揚力,誘導抗力を計算する

循環\(\Gamma\),揚力\(L\),誘導抗力\(D_{i}\)はそれぞれ次のように計算できる

\begin{eqnarray}
\Gamma_{i} &=& \frac{1}{2}\rho c_{i} C_{l_{i}} \\
L &=& \sum \rho V_{i}\Gamma_{i}\Delta y\cos{(\Gamma_{dihedral_{i}})} \\
D_{i} &=& \sum \rho w_{i_{i}} \Gamma_{i} \Delta y \\
\end{eqnarray}

これで循環分布を計算できた

このプロシージャについては以上である

曲げモーメントを計算する

最後に翼の変形量を求めるのに必要な値を順次計算していく

まず,翼にはたらく曲げモーメント,剪断力,トルクを計算する

これらの値がパネルの両端(sp)で定義される量であることに注意すると,i番目(左翼)のspにおける曲げモーメント\(M_{bending_{i}}\),剪断力\(F_{i}\),トルク\(M_{torque_{i}}\)は次のように計算できる

\begin{eqnarray}
M_{bending_{i}} &=& \sum_{j=1}^i ((\Delta N_{j-1}\cos{(\Gamma_{dihedral_{j-1}})}-\Delta W_{i})\times|y_{cp_{j-1}}-y_{i}|+\Delta N_{j-1}\sin{(\Gamma_{dihedral_{j-1}})}\times|z_{cp_{j-1}}-z_{i}|) \\
F_{i} &=& \sum_{j=1}^i (\Delta N_{j-1}\cos{(\Gamma_{dihedral_{j-1}})}-\Delta W_{i}) \\
M_{torque_{i}} &=& \sum_{j=1}^i (\Delta M_{torque_{j-1}}+\Delta T \times|z_{cp_{j-1}}-z_{i}|) \\
\end{eqnarray}

ここで,トルク計算の右辺第二項の\(\Delta T \times|z_{cp_{j-1}}-z_{i}|\)は,上反角が大きい機体における空気力のx軸方向成分とパネル同士のZ方向距離によるトルクである

XFLR5 Guidelinesのp.26にある次の記述

The only use for the sweep and thedihedral in this implementation of the LLT is for the calculation of the pitching moment coefficient Cm.

XFLR5 v6.02 Guidelines page26

(後退角および上反角はCmの計算にしか使わない) とはこのことである

この計算方法だと翼中心の値が2倍になってしまうので,ループの外でその値を半分にする処理を行っている

このプロシージャについては以上である

たわみ,ねじれの計算

たわみ角,たわみ,ねじれ角の計算を行う

たわみ角\(\theta\),たわみ\(w\),ねじれ角\(\phi\)は次の式で計算できる

\begin{eqnarray}
\theta &=& \int \frac{M}{EI} dy \\
w &=& \int \theta dy \\
\phi &=& \int \frac{M_{torque}}{GJ} dy \\
\end{eqnarray}

これで翼の変形量を計算できた

このプロシージャについては以上である

収束判定

ループ内の処理の最後に収束判定を行う

繰り返し計算の揚力の変化量があらかじめ指定した閾値を下回ったとき判定とみなす

これでループ内の処理は終了である

構造体へ格納

もろもろの値を計算して構造体へ格納する

\begin{eqnarray}
\Delta S_{i} &=& c_{cp_{i}}dy\cos{(\Gamma_{dihedral_{i}})} \\
S &=& \sum_{i=0}^{2n-1} \Delta S_{i} \\
\overline{c} &=& \frac{2}{S}\sum_{i=0}^{n-1} c_{i}^2 \Delta y \cos{(\Gamma_{dihedral_{i}})}
\tag{文献3(2.52)} \\
\overline{y} &=& \frac{2}{S}\sum_{i=0}^{n-1} c_{i}y_{i} \Delta y \cos{(\Gamma_{dihedral_{i}})}
\tag{文献3(2.54)} \\
AR &=& \frac{b^2}{S}
\tag{文献3(2.40)} \\
C_{l_{\alpha}} &=& \frac{\partial}{\partial\alpha}(\frac{\sum_{i=0}^{2n-1}\frac{1}{2}\rho V^2 \Delta S C_{l_{i}}}{\frac{1}{2}\rho V^2 S}) \\
&=& \frac{\sum_{i=0}^{2n-1} C_{l_{\alpha_{i}}} \Delta S_{i}}{S} \\
D_{p} &=& \sum_{i=0}^{2n-1} \Delta D_{p_{i}} \\
L &=& -\sum_{i=0}^{2n-1} dN_{i}\cos{(\Gamma_{dihedral_{i}})}\times y_{cp_{i}} \\
M_{ac} &=& \sum_{i=0}^{2n-1} \Delta M_{ac_{i}}+dT_{i}\times z_{cp_{i}} \\
M_{cg} &=& M_{ac}+L_{wing}\times(h_{spar}-h_{ac}-\Delta h) \\
N &=& \sum_{i=0}^{2n-1} dT_{i}\times y_{cp_{i}} \\
D &=& D_{i}+D_{p} \\
C_{L} &=& \frac{L}{\frac{1}{2}\rho V^2 S} \\
C_{D_{p}} &=& \frac{D_{p}}{\frac{1}{2}\rho V^2 S} \\
C_{D_{i}} &=& \frac{D_{i}}{\frac{1}{2}\rho V^2 S} \\
C_{D} &=& \frac{D}{\frac{1}{2}\rho V^2 S} \\
C_{m_{ac}} &=& \frac{M_{ac}}{\frac{1}{2}\rho V^2 S} \\
C_{m_{cg}} &=& \frac{M_{cg}}{\frac{1}{2}\rho V^2 S} \\
e &=& \frac{\frac{C_{L}^2}{\pi AR}}{C_{D_{i}}} \\
a_{w} &=& \frac{C_{l_{\alpha_{i}}}}{1+\frac{C_{l_{\alpha}}}{\pi AR}} \\
C_{y_{\beta}} &=& -\frac{2}{S} \sum_{i=0}^{n-1} C_{l_{\alpha_{i}}}\sin{(\Gamma_{dihedral_{i}})}^2 \Delta S_{i} \\
C_{l_{\beta}} &=& -\frac{2}{Sb}\sum_{i=0}^{n-1}C_{l_{\alpha_{i}}}\sin{(\Gamma_{dihedral_{i*}})} y_{cp_{i}} \Delta S_{i}
\tag{文献3(4.65′)} \\
C_{l_{p}} &=& -\frac{4}{Sb^2}\sum_{i=0}^{n-1} C_{l_{\alpha_{i}}}y_{cp_{i}}^2 \Delta S_{i}
\tag{文献3(4.76′)} \\
C_{n_{p}} &=& -\frac{4}{Sb^2}\sum_{i=0}^{n-1} (C_{l_{i}}-C_{d_{\alpha_{i}}})y_{cp_{i}}^2 \Delta S_{i}
\tag{文献3(4.86′)} \\
C_{l_{r}} &=& -\frac{8}{Sb^2}\sum_{i=0}^{n-1} C_{l_{i}}y_{cp_{i}}^2 \Delta S_{i}
\tag{文献3(4.94′)} \\
C_{n_{r}} &=& -\frac{8}{Sb^2}\sum_{i=0}^{n-1} C_{d_{i}}y_{cp_{i}}^2 \Delta S_{i}
\tag{文献3(4.102′)} \\
\end{eqnarray}

ここで,\(C_{y_{\beta}}<0\)とは,機体が対気速度\(V\),正の横滑り角\(\beta\)で飛行しているときにy軸負の方向に行う並進運動を表した微係数であり,次の式で計算できる

\begin{eqnarray}
Y &=& \int_{\frac{-2}{b}}^{\frac{2}{b}} ((\frac{1}{2}\rho V^2 \Delta S C_{l_{\alpha}} \sin{(\Gamma_{dihedral})}\beta)\times\sin{(\Gamma_{dihedral})}) \\
C_{y} &=& \frac{Y}{\frac{1}{2}\rho V^2 S} \\
C_{y_{\beta}} &=& \frac{\partial C_{y}}{\partial\beta} \\
\end{eqnarray}

これは機体が横風に押し流される効果を表すものであり,主翼の上反角が大きい機体特有のもので文献3には載っていない

このプロシージャについては以上である

値を返す

関数に値を返す

結果の表示

結果をシートに表示する

できる限り配列に格納してまとめて貼り付ける

これでプログラムについての解説は終了である

プログラムの実行

それでは実際にプログラムを実行してみる

テンプレートファイルの下準備

前回(≫TR-797の実装)の続きからやりたい人は,まず主翼のシートのBB20:BB27に,解析の条件を入力する

たわみ・ねじりのシートのN4:N16にもろもろの値を各種パラメータのシートから持ってくる

列L~列Nの45行目以下に,曲げ剛性,ねじれ剛性,翼重量をそれぞれスパー寸法のシートと主翼のシートから持ってくる

翼重量の単位は[N]なので,主翼重量に先ほど持ってきた重力加速度をかけておく

これで準備は完了である

実行結果

実行用のプログラムは例えば次のようなものになる

これを実行すれば,数秒で計算が終了する

たわみ・ねじれのシートにさらにいくつかの値を主翼シートから持ってくる

たわみやねじれ,SFDやBMDも表示できる

繰り返し回数

繰り返し回数による結果の変化を見てみる

上図より,数十回の繰り返し計算で値が収束していることがわかる

下図を参考に変化量の閾値を決めるといい

XFLR5との比較

次の記事を参考にしてほしい

XFLR5で三次元翼を解析して結果を出力する

まとめ

Excel VBAで主翼の空力・構造計算のプログラムを作成した

XFLR5を使うよりもかなり時間の節約になるのではないかと思う


TR-797の実装に戻る

渦格子法を用いた尾翼の空力計算に進む

記事一覧に戻る

【更新情報】
20191219
・重心周りのピッチングモーメントの式の間違いを修正

Sponsored Link

コメント

タイトルとURLをコピーしました