TR-797の実装

TR-797”非平面翼の最適設計-揚力と翼根曲げモーメントを与えた時の最小誘導抵抗-”のExcel VBAでの実装を行う

目的とフローチャート

論文の解説自体は前回(≫TR-797の解説)で行ったので,今回はExcel VBAで実装してみる

テンプレートファイルのダウンロード

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

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

前々回(≫重量推算②)の続きからやりたい人は,まずセルAW12:AW15に構造的制限の有無,翼に沿った長さ\(l_{e}\),パネルの反復\(\Delta S\),構造的制限\(\beta\)を入力,セルBB22:BB23に対気速度と高度を入力する

またAF列の45行目以下に,テンプレートシートから上反角分布をもってきておく

フローチャート

循環分布を最適化するまでのフローチャートを以下に示す

機体形状を定式化し,循環-吹きおろしの影響係数を求める.揚力のつり合い,翼根曲げモーメントの制限,誘導抗力が最小になるという条件を詰め込んだ連立一次方程式を解くことで循環分布が求まり,それをもとに誘導抗力などを計算する

どちらかというと渦格子法に近いプログラムである

1つずつ解説していく

変数の宣言

今回は翼分割数span_div,構造的制限beta,高度hEを引数として,循環分布,吹きおろし分布を戻り値として返すFunctionとする

【VBA入門】Functionの使い方(呼び出し、引数、戻り値)

使用している変数を以下に示す

これからのプログラムはシートをまたいで実行されることが多い.そういうときはワークシートプロジェクトを定義する

【VBA入門】WorksheetsからWorksheetオブジェクトを取得し操作する

さすがに変数が増えてきて,動的変数を再定義するだけで行数がばかにならないので可読性を上げるために配列の宣言と動的変数の再定義を同じ行で行っている

値の読み込み

必要な値を読み込む

シートをまたいで値を読み込むときにはWithステートメントを使うとよい

【VBA入門】Withの使い方、入れ子(ネスト)で使う方法

状態変数を読み込むときに使っている関数は次のようなものである

構造体については次の記事を参考にしてほしい

設計シートに使う構造体について

構造体は便利なのでこれからもどんどん使っていく

式(9)の計算

式(9)を計算する

\begin{eqnarray}
y’_{ij} &=& (y_{i}-y_{j}) cos(\phi_{j})+(z_{i}-z_{j}) sin(\phi_{j}) \\
z’_{ij} &=& -(y_{i}-y_{j}) sin(\phi_{j})+(z_{i}-z_{j}) cos(\phi_{j}) \\
y”_{ij} &=& (y_{i}+y_{j}) cos (\phi_{j}) -(z_{i}-z_{j}) sin (\phi_{j}) \\
z”_{ij} &=& (y_{i}+y_{j}) sin (\phi_{j}) +(z_{i}-z_{j}) cos (\phi_{j}) \\
{y’_{ij}}_{mirror} &=& (y_{i}-y_{j}) cos(-\phi_{j})+(z_{i}-(z_{j}-2h_{E})) sin(-\phi_{j}) \\
{z’_{ij}}_{mirror} &=& -(y_{i}-y_{j}) sin(-\phi_{j})+(z_{i}-(z_{j}-2h_{E})) cos(-\phi_{j}) \\
{y”_{ij}}_{mirror} &=& (y_{i}+y_{j}) cos (\phi_{j}) +(z_{i}-(z_{j}-2h_{E}) sin (\phi_{j}) \\
{z”_{ij}}_{mirror} &=& -(y_{i}+y_{j}) sin (\phi_{j}) +(z_{i}-(z_{j}-2h_{E})) cos (\phi_{j}) \\
\end{eqnarray}

式(10)の計算

式(10)を計算する

\begin{eqnarray}
{R_{+ij}}^2 &=& (y’_{ij}-\Delta S_{j})^2+{z’_{ij}}^2 \\ {R_{-ij}}^2 &=& (y’_{ij}+\Delta S_{j})^2+{z’_{ij}}^2 \\ {R’_{+ij}}^2 &=& (y”_{ij}+\Delta S_{j})^2+{z”_{ij}}^2 \\ {R’_{-ij}}^2 &=& (y”_{ij}-\Delta S_{j})^2+{z”_{ij}}^2 \\ {{R_{+ij}}_{mirror}}^2 &=& ({y’_{ij}}_{mirror}-\Delta S_{j})^2+{{z’_{ij}}_{mirror}}^2 \\
{{R_{-ij}}_{mirror}}^2 &=& ({y'{ij}}_{mirror}+\Delta S_{j})^2+{{z’_{ij}}_{mirror}}^2 \\{{R’_{+ij}}_{mirror}}^2 &=& ({y”_{ij}}_{mirror}+\Delta S_{j})^2+{{z”_{ij}}_{mirror}}^2 \\
{{R’_{-ij}}_{mirror}}^2 &=& ({y”_{ij}}_{mirrro}-\Delta S_{j})^2+{{z”_{ij}}_{mirror}}^2 \\
\end{eqnarray}

式(8)の計算

式(8)を計算する

\begin{eqnarray}
Q_{ij} &=& \frac{1}{2\pi}((-\frac{y’_{ij}-\Delta S_{j}}{{R_{+ij}}^2}+\frac{y’_{ij}+\Delta S_{j}}{{R_{-ij}}^2}) cos(\phi_{i}-\phi_{j}) \\
&& +(-\frac{z’_{ij}}{{R_{+ij}}^2}+\frac{z’_{ij}}{{R_{-ij}}^2}) sin(\phi_{i}-\phi_{j}) \\
&& +(-\frac{y”_{ij}-\Delta S_{j}}{{R’_{-ij}}^2}+\frac{y”_{ij}+\Delta S_{j}}{{R’_{+ij}}^2}) cos(\phi_{i}+\phi_{j}) \\
&& +(-\frac{z”_{ij}}{{R’_{-ij}}^2}+\frac{z”_{ij}}{{R’_{+ij}}^2}) sin(\phi_{i}+\phi_{j}) \\
&& -(-\frac{{y’_{ij}}_{mirror}-\Delta S_{j}}{{{R_{+ij}}_{mirror}}^2}+\frac{{y’_{ij}}_{mirror}+\Delta S_{j}}{{{R_{-ij}}_{mirror}}^2}) cos(\phi_{i}-(-\phi_{j})) \\
&& -(-\frac{{z’_{ij}}_{mirror}}{{{R_{+ij}}_{mirror}}^2}+\frac{{z’_{ij}}_{mirror}}{{{R_{-ij}}_{mirror}}^2}) sin(\phi_{i}-(-\phi_{j})) \\
&& -(-\frac{{y”_{ij}}_{mirror}-\Delta S_{j}}{{{R’_{-ij}}_{mirror}}^2}+\frac{{y”_{ij}}_{mirror}+\Delta S_{j}}{{{R’_{+ij}}_{mirror}}^2}) cos(\phi_{i}-\phi_{j}) \\
&& -(-\frac{{z”_{ij}}_{mirror}}{{{R’_{-ij}}_{mirror}}^2}+\frac{{z”_{ij}}_{mirror}}{{{R’_{+ij}}_{mrror}}^2}) sin(\phi_{i}-\phi_{j})) \\
\end{eqnarray}

頑張って入力しよう

がんばった

式(20),(21),(22)の計算

式(20),(21),(22)を計算する

\begin{eqnarray}
c_{i} &=& 2cos{\phi_{i}} \Delta \sigma
\tag{20} \\
b_{i} &=& \frac{3\pi}{2}(\eta_{i} cos{\phi_{i}}+\zeta_{i} sin{\phi_{i}})\Delta \sigma
\tag{21} \\
A_{ij} &=& \pi q_{ij} \Delta \sigma_{i}
\tag{22} \\
\end{eqnarray}

式(33)の計算

式(33)を計算する

\begin{eqnarray}
\left[
\begin{array}{cccccccc}
2A_{11} & A_{21}+A_{12} & \ldots & \ldots & \ldots & A_{N1}+A_{1N} & -c_{1} & -b_{1} \\
A_{12}+A_{21} & 2A_{22} & & & & \vdots & -c_{2} & -b_{2} \\
\vdots & \vdots & \ddots & & & \vdots & \vdots & \vdots \\
A_{1i}+A_{i1} & \vdots & \ldots & 2A_{ii} & \ldots & \vdots & \vdots & \vdots \\
\vdots & \vdots & & & \ddots & \vdots & \vdots & \vdots \\
A_{1N}+A_{N1} & \ldots & \ldots & \ldots & \ldots & 2A_{NN} & -c_{N} & -b_{N} \\
-c_{1} & -c_{2} & \ldots & \ldots & \ldots & -c_{N} & 0 & 0 \\
-b_{1} & -b_{2} & \ldots & \ldots & \ldots & -b_{N} & 0 & 0 \\
\end{array}
\right]
\cdot \left[
\begin{array}{c}
g_{1} \\
g_{2} \\
\vdots \\
\vdots \\
\vdots \\
g_{N} \\
\mu_1 \\
\mu_2 \\
\end{array}
\right]
= \left[
\begin{array}{c}
0 \\
0 \\
\vdots \\
\vdots \\
\vdots \\
0 \\
-1 \\
\beta \\
\end{array}
\right]
\end{eqnarray}

ここでは構造的制限の有無で場合分けし,構造的制限のない場合はN+1次元連立1次方程式を解いている

式(19)の計算

連立方程式を解いて\({\bf g}\)が求まったので,それをもとに循環\({\bf \Gamma}\)を計算する

\begin{eqnarray}
g_{i} &=& \frac{2 l_{e} \rho U \Gamma_{i}}{L}
\end{eqnarray}

式(7)の計算 式(4),(5),(6)の計算

循環\({\bf \Gamma}\)がわかったので,式(7)を用いて吹きおろしを計算する

\begin{eqnarray}
V_{ni} &=& \sum_{j=1}^{N} Q_{ij} \Gamma_{i} \\
\end{eqnarray}

ただしこの吹きおろしはTrefftz面での吹きおろしなので,翼面上での吹きおろしはこれの半分であることに注意する

最後に式(4),式(5),式(6)を計算する

\begin{eqnarray}
L &=& 4\rho U \sum_{i=1}^{N} \Gamma_{i} cos{\phi_{i}} \Delta S_{i} \\
B &=& 2\rho U \sum_{i=1}^{N} \Gamma_{i} (y_{i}cos{\phi_{i}}+z_{i}sin{\phi_{i}}) \Delta S_{i} \\
D_{i} &=& 2\rho U \sum_{i=1}^{N} \Gamma_{i} V_{ni} \Delta S_{i} \\
\end{eqnarray}

結果の出力

計算した結果を出力する

循環分布と吹きおろし分布を戻り値として返している

以上でプログラムの解説は終了である

プログラムの実行

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

実行には次のコードを用いる

実際に実行してみると次のようになる

楕円分布よりもさらに翼根側に揚力が集中した循環分布が計算できた

まとめ

地面効果を含めた循環分布の最適化ができた

あとは煮るなり焼くなり好きに最適化しよう


TR-797の解説

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

記事一覧に戻る

Sponsored Link

コメント

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