空力実装編 — `02_空力理論編.md` のVBA実装¶

02_空力理論編.mdで説明した揚力線理論の各要素が、実際のVBAコードのどこに対応するかをまとめます。各節の見出しは、コード側の' ●学習資料:アンカーコメントと対応しています。

1. 翼形状の展開(コサイン分布グリッドの生成)¶

💻 対応コード: clsWingGeometry.cls の Build

(コサイン分布グリッドがどう見えるかの図は02_空力理論編.md§2を参照)

「パラメータ_翼」シートの4つの点列(平面系・ねじり下げ分布・桁位置・翼型分布)は、Y座標が不揃いの疎な点列です。Buildは、これを理論編§2で説明したコサイン変数変換

m_Th(i) = pi * (i - 1) / n
m_Y(i)  = -BSpan / 2 * Cos(m_Th(i))

でn+1点(nは「パラメータ_翼」シートの「解析分割数」で指定、Get01_翼解析分割数経由、既定n=1000なので1001点)のコサイン分布グリッドへ展開し、各点の翼弦長・迎角・桁位置・翼型をMod03_分布補間の線形補間で求めます。

左右対称な翼という前提のもと、半分側(翼根θ=π/2〜翼端θ=π)だけ点列から補間し、もう半分(θ=0〜π/2)はそれを鏡映(ミラー) して埋めています(For i = 1 To K のループ)。これは理論編§2の「左右対称な翼は片側の情報だけで決まる」という性質をそのままデータ構造に反映したものです。

翼根位置(θ=π/2)は理論上Y=0になるはずですが、浮動小数点誤差で厳密な0からわずかにズレることがあるため、m_Y(K + 1) = 0で明示的に0へスナップしています。

2. フーリエ係数の連立方程式¶

💻 対応コード: Mod05_空力ソルバー.bas の SolveAero(A・Cvecの組み立てとMInverse/MMultの箇所)

理論編§3のモノプレーン方程式は、コード中で次のように組み立てられます(p, qは制御点・フーリエ項のインデックス、K個)。

Cvec(p, 1) = minfP * cordP / (4 * geo.BSpan) * Sin(thP) * abaP
A(p, q) = (XX(q) * minfP * cordP / (4 * geo.BSpan) + Sin(thP)) * Sin(XX(q) * thP)

minfPが局所揚力傾斜a₀、cordPが局所翼弦、abaPが(α−α_L0)に相当する「絶対迎角」です。XX(q) = 2q-1が奇数次のフーリエ次数nそのもの(理論編§2の「奇数次だけで足りる」に対応)です。

行列Aの逆行列とベクトルCの積で係数ベクトルBを求める部分は、このプロジェクト全体の方針([[MInverse/MMultの利用方針]])どおり、VBAで連立一次方程式を解く際は自前のガウス消去法などを書かず、Excelワークシート関数Application.WorksheetFunction.MInverse/MMultを呼び出しています。

Dim AInv As Variant: AInv = Application.WorksheetFunction.MInverse(A)
Bvec = Application.WorksheetFunction.MMult(AInv, Cvec)

3. 誘導迎角の分布復元とCL・スパン効率¶

💻 対応コード: Mod05_空力ソルバー.bas の SolveAero(CL/del_の計算、bn(j)ループ)

CL = pi * geo.AR * Bvec(1, 1)
del_ = delSum / (Bvec(1, 1) ^ 2)     ' delSum = Σ XX(q)*Bvec(q,1)^2  (q=2..K)

が理論編§4のCL = π・AR・A₁、δの式にそのまま対応します。

制御点だけでなく全n+1点での吹き下ろし・誘導迎角・局所揚力係数の復元(理論編§5)は、bn(j)ループで行われます。

bn(j) = Σ XX(q) * Bvec(q,1) * Sin(XX(q) * geo.Th(j))  / Sin(geo.Th(j))
W(j)  = V1 * bn(j)
Cl_dist(j) = minfJohen(j) * (aba(j) - W(j) / V1)

翼根(j=1)・翼端(j=n+1)はSin(Th)=0で数式が特異点になるため、bn(1) = bn(n+1) = 0.2という固定値で近似しています(元MATLABの実装をそのまま踏襲した経験則的な処理)。

4. プロファイル抗力・モーメント係数(翼型データベース参照)¶

💻 対応コード: Mod04_翼型データベース.bas、Mod05_空力ソルバー.bas の GetClParams/GetCd/GetCm 呼び出し箇所

理論編§6の「2次元翼型データのテーブル参照」は、Mod04_翼型データベースが担当します。翼型ごとに{翼型名}_Cd.csv/_Cm.csv/_Cl.csvをルックアップ用のクラス(clsGridTable/clsClTable)へ読み込み、局所Reynolds数・局所有効迎角をキーに補間参照します。CSVは初回参照時に遅延読み込みし、Scripting.Dictionaryでキャッシュします(同じ翼型を1001点全部で毎回ファイルから読み直すと遅いため)。

2つの翼型をブレンドする場合(ratio1)、SolveAero側で

CdDisJohen(i) = cd1 * geo.ratio1(i) + cd2 * (1 - geo.ratio1(i))

のように加重平均しています。

5. 圧力中心・ピッチングモーメントの積分¶

💻 対応コード: Mod05_空力ソルバー.bas の SolveAero(Cp(i)、MG(i)/Cmgの計算)

Cp(i) = 0.25 - CmDisJohen(i) / Cl_dist(i)
MG(i) = Lift(i) * (geo.KetaIti(i) - Cp(i)) * geo.Cord(i)

が理論編§7の圧力中心・ピッチングモーメントの式にそのまま対応します。geo.KetaIti(i)(桁位置、コード上は0〜1の割合)とCp(i)(圧力中心、翼弦割合)の差に翼弦長を掛けることで、桁位置を基準とした局所揚力のモーメントアーム[m]を求めています。

6. トリム速度の反復ループ¶

💻 対応コード: Mod05_空力ソルバー.bas の SolveAero(For iter = 1 To MAX_ITER ループ全体)

理論編§8のトリムループは、SolveAero関数の外側のFor iterループそのものです。

If Abs(Wtot - Weight) <= 0.001 Then Exit For
V1 = Sqr(2 * Weight / (RO * geo.s * CL_sum))

総揚力Wtot[kgf]が機体重量Weight[kgf]に0.001kgf以内で一致すれば収束とみなし、そうでなければV1 ∝ 1/√CL_sumの関係で飛行速度を更新して再計算します(最大MAX_ITER=200回)。ループの中で§2〜§5すべての計算(フーリエ係数の求解、CL、揚力分布)を毎回丸ごとやり直している点に注意してください — 効率的ではありませんが、元のMATLABコードの構成をそのまま踏襲しています。

7. 補間ユーティリティ¶

💻 対応コード: Mod03_分布補間.bas

Build(§1)が使う補間処理をまとめたモジュールです。InterpLinearは数値の点列の単純な線形補間、InterpAirfoilは翼型名(文字列)の「区間内で2つの翼型をブレンドする比率」を求める処理(理論編§6のブレンドの元ネタ)です。範囲外のクエリは境界値でクランプします。

SplitTwoColumns/SplitTwoColumnsTextは、シートのRange.Valueから得られる2次元Variant配列を、Y列・値列の1次元配列に分割するだけの補助関数です。

8. 結果の保持と出力¶

💻 対応コード: clsAeroResult.cls、Mod07_空力結果出力.bas

SolveAeroの計算結果(スカラー値のCL/CD/V1等、分布データのCl_dist/Lift/Drag等)はclsAeroResultインスタンスに詰められて返ります。Mod07_空力結果出力.WriteAeroResultToSheetは、この内容を目視確認用に「結果_空力」シートへ書き出すための関数で、揚力線理論の計算そのものには関与しません(01_全体アーキテクチャ.mdで説明したとおり、本番の構造解析パイプラインはこのシート出力を経由せずメモリ上のオブジェクトを直接使います)。

9. 実行エントリポイントとテスト¶

💻 対応コード: Mod08_空力実行.bas、Mod09_空力迎角スイープ.bas、Mod06_空力ソルバーテスト.bas

Mod08_空力実行.SolveAeroFromSheet — 「パラメータ_翼」系シートを読み、clsWingGeometry.Build→SolveAeroを1回実行して結果を返す共通処理。単一迎角実行(RunAero)・構造解析(Mod18)の両方から呼ばれる、この空力ソルバーへの唯一の入口です。
Mod09_空力迎角スイープ.RunAeroSweep — 「パラメータ設定」シートの迎角範囲を1本ずつSolveAeroFromSheetへ渡し、結果を2次元配列にまとめてイミディエイトウィンドウへ出力する(迎角ごとの特性を横に並べて比較したいときに使う)。ループはFor...StepではなくDo...Loop(刻み幅に対する相対許容誤差つきのExit Do)で書かれている。理由は§10を参照。
Mod06_空力ソルバーテスト.Test_AeroSolver — シートを介さず、コード内に直接埋め込んだ既知の入力データでSolveAeroを実行し、Python参照実装で計算しておいた期待値と突き合わせるテスト。揚力線理論の実装自体が壊れていないかを、シートの状態に左右されずに確認できます。

10. 結果のCSV出力¶

💻 対応コード: Mod20_空力結果CSV出力.bas、Mod19_出力共通.bas

元のMATLABプログラム(WING_IKI.m/geikaku_henka.m)は解析結果をCSVファイルへ書き出していました。このVBA移植版でも、ユーザー提供の共通ユーティリティOutputText(2次元配列→テキストファイル、文字コード・区切り文字を指定可能)を使って同等のCSVを再現しています。

元MATLABのファイル	VBA側の対応	呼び出し元
`output_wing_analyze_data.csv`	`WriteAeroSingleCsv`	`Mod08_空力実行.RunAero`
`heimenkei,aoa=X.csv`	`WriteAeroSweepDetailCsv`	`Mod09_空力迎角スイープ.RunAeroSweep`(迎角ごと)
`output_geikaku_henka.csv`	`WriteAeroSweepSummaryCsv`	`Mod09_空力迎角スイープ.RunAeroSweep`(ループ後に1回)

出力先はThisWorkbook.Path & "\output(機体名)\"(Mod19_出力共通.EnsureOutputFolder、機体名はGet01_機体名()、無ければ自動作成)、文字コードはShift-JIS固定(InputTextを使う既存の入力CSVと統一)です。列ごとに幅が異なる「分布本体+空行+key,valueフッター」という構成のCSVは、行をCollectionにいったん積んでから最大列数へ空文字パディングするMod19_出力共通.PadRowsToArray2Dで1回のOutputText呼び出しにまとめています。

元コードには2件の既知の不具合がありましたが、ここでは再現せず正しい値を出力しています。

全機誘導抗力係数: 元コードは局所値(翼端1点のCd_ind)を全機値として誤って使っていましたが、ここではclsAeroResult.Cdi_overall(=CL²(1+δ)/π/AR、Mod05_空力ソルバーが別途正しく計算・保持している全機スカラー値)を使います。
翼弦長・区間長の内訳: 元コードは区分点が3個であることを前提にハードコードしていたため、区分点数が3以外だと不正確でした。ここでは実際の区分点数に応じた可変長リストとして出力します(AppendBreakpointRows)。

「吹き降ろし分布」「誘導抗力係数分布」等、clsAeroResultに直接保持されていない量は、保持済みの値(Ae/ew/Cl_dist等)からMod05_空力ソルバー内部と同じ式で再構成しています(BuildAeroDistRow)。ヘッダーの文言(「縦ゆれモーメント係数」等)は元MATLABの表記に合わせてある(構造側05_構造実装編.md§11も同様)。

出力するファイルの取捨選択は、「設定」シートの実データチェック機構(Sh01_設定の「出力チェック」列、Enum_O出力データ一式という専用Enumで11種類のCSVをO01〜O11として管理)で行います。Mod19_出力共通.GetOutputCSVSetting()が1回の呼び出しで11種類ぶんのBoolean配列を返し、各Write*CsvはoutputSetting(Enum_O出力データ一式.O0X_...)で判定します(構造側も含め、以前使っていた文字列キー方式のShouldOutputはこちらに一本化して削除済み)。

RunAero/RunAeroSweepは実行完了後、ユーザー作成のMsgYesNoで完了と出力先フォルダを表示し、「はい」ならユーザー作成のOpenFolderで出力先フォルダをエクスプローラーで開きます。

空力実装編 — 02_空力理論編.md のVBA実装¶