中で何を計算しているか
本プログラムは目的に応じて 6 種類の計算エンジンを切り替えられます。 「シミュレータ(F) → シミュレータ設定」で選びます。
| エンジン | 自由度 | 内容 |
|---|---|---|
| 3DOF 質点 | 3 | 重心位置 X・Y とピッチ角の 3 自由度。もっとも軽く、素早く傾向を見たいとき向きです。 |
| 3DOF | 3 | 縦 3 自由度に尾翼操作を加えたもの。ダウンウォッシュの遅れも扱います。 |
| 6DOF | 6 | 並進 3 + 回転 3 の完全な 6 自由度。補助翼・方向舵によるロールとヨーが効きます。 |
| 非線形 6DOF | 6 | 翼を細かく分割し、断面ポーラを使って各要素の空気力を積み上げる非線形モデル |
| 6DOFNLFS | 6 | 非線形飛行シミュレーション系のモデル |
| 6DOFGame | 6 | クォータニオンを使ったゲーム向け 6DOF（Physics for Game Developers 参照） |
数値積分には 4 次のルンゲ＝クッタ法を使い、毎フレームの位置・速度・加速度を求めています。
主翼は翼幅方向に細かく翼要素へ分割し、各要素ごとに、 そのときのねじれ・上反角・地面効果を考慮した揚力係数を求めて計算しています。
翼型の断面特性は、翼型を割り当てていれば XFOIL などで実際に計算したポーラをそのまま使います。 迎角と Reynolds 数の両方で補間するので、飛行速度が変われば断面特性も変わります。
各要素の翼弦長から局所 Reynolds 数を求め、上反角 Γ については cos Γ の成分と投影翼幅で誘導角を補正しています。
揚力線理論群では、求めた循環 Γ から Kutta–Joukowski の定理で翼幅荷重を出しています。
地面効果はアスペクト比の補正という形で取り込んでいます。 高さ h、翼幅 b、もとのアスペクト比 A に対して、有効アスペクト比 gA は次のようになります。
誘導角にはこの逆比 A/gA を掛けます。地面すれすれ（h → 0）では誘導角が 0 に近づき、 十分高い（h ≫ b）ところでは 1 に戻って地面効果が消えます。
参考文献：『航空工学（下）』日本航空技術協会
たわみは、たわみによって生じる各翼要素の揚力方向と位置が毎フレーム変化していく 幾何学的非線形問題として解いています。
具体的には、各翼要素のモーメントを求め、モールの定理によりこのモーメントを 曲げ剛性 EI で割った値 M/EI を仮想荷重とし、再度モーメントを解くことで得た値を たわみとしています。
また、このときに生じる各翼のアスペクト比も毎フレームごとに求め直しています。
翼型解析には 9 種類のエンジンを用意しています。大きく 3 系統に分かれます。
| XFOIL 6.99 | M. Drela によるパネル法＋境界層法。低速翼型の標準的な解析ツールです。 抵抗・失速まで求まりますが、Reynolds 数の指定が必要で、計算時間もかかります。 |
| Eppler PROFIL 1.1 | Eppler のパネル法と境界層法による粘性・圧縮性解析 |
| Theodorsen NACA TR-452 |
Theodorsen・Garrick 理論。Reynolds 数に依存しません。 |
| 今井 1942 |
今井功『任意翼型の理論』の逐次近似法 |
| 守屋 1937/1938/1941 |
守屋富次郎の任意翼型理論。 輪郭上に渦を配列した積分方程式(1937)、翼型を Fourier 級数で表す第一次近似(1938)、 円への逐次写像による正解(1941)の 3 法を選べます。 圧力分布は補遺(1942)の訂正に従います。 |
| 森口 1938 |
森口繁一の共役 Fourier 積分・逐次写像による解析 |
| 水野 BEM SNG_B |
線形分布渦による非粘性境界要素法。最大 49 要素で解析します。 |
| Eppler PROFIL TRAPRO |
Eppler の等角写像法。設計迎角 α* の分布から翼型を作ります。 |
| 佐藤 1966 |
厳密二次元ポテンシャル理論。目標流速分布 qα(θ) や 写像関数 g(θ) から翼型を直接作ります。 |
| Lighthill A.R.C. R&M No.2112 |
log q0(θ) の Fourier 級数から、式(7) → 式(10) → 式(8) の順に逆設計します。 |
Qspec逆設計には 2 つの解き方があります。どちらも、描いた速度分布にいちばん近い翼型を Levenberg–Marquardt 法（減衰つきガウス・ニュートン法）で探します。 指定が実現不可能でも、最小二乗がいちばん近い形を選ぶので計算は破綻しません。
| 解法 | 中で解いているもの |
|---|---|
| 混合逆解法 XFOIL QDES 相当 |
流れは二次元パネル法（線形渦分布＋Kutta 条件）で解きます。表面速度は節点の渦の強さそのもので、 揚力は表面圧力の積分で出します。形は、書き換えた区間の節点を法線方向に δ(t) = sin²(πt)·Σ ck cos(kπt) だけ動かします。sin²(πt) が区間の両端で値も傾きも 0 にするので、 動かさない部分との継ぎ目に折れ目ができません。動かせる範囲は合わせる範囲より前後に 6 割広く取り、 速度分布の山の両脇に出るへこみを作る余地を残しています。 |
| 完全逆解法 XFOIL MDES 相当 |
Theodorsen の写像で、翼型を円周上の関数 ψ(φ) のフーリエ係数で表します。 表面速度は Q = N(φ)·s(φ) に分かれ、s は迎角にも零揚力角にもよらないので、 u = ln s を設計変数にして係数を解きます。後縁は元の翼型の値に拘束します。 薄い翼型のために、式(13) は逐次代入ではなく Newton 法で解きます。 |
UIUC/Selig 形式の翼型座標 1,665 点を SQLite3 データベースに収めています。 座標だけでなく、翼厚比・キャンバー比・点数といった幾何量も計算して格納しているので、 「翼厚 12 % 前後の翼型」といった絞り込みができます。
解析結果もこのデータベースに記録されます。同じ条件で再度実行すると、 計算せずに前回の結果を返します。
全翼型バッチ解析・GA 探索・GA最適化は、複数の CPU コアを使います。 この 3 つは同時には 1 つしか動かしません。
| XFOIL のバッチ | 1 ケース = xfoil.exe 1 プロセスなので、スレッドで並列にします |
| 理論ソルバーのバッチ | Python で計算するため、プロセスプールへ投げて GUI が止まらないようにしています |
| GA 探索 | 個体の評価をプロセスプールで並列に行います |
| GA最適化 | 1 個体（1 回の飛行）をプロセスプールで並列に飛ばします。エリートは前の世代の成績を使い回し、飛ばし直しません。 停止を押すと、飛行中の個体も 50 ステップごとに打ち切りを確かめて、すぐに止まります |
ワーカープロセスは親プロセスを監視しており、 親プロセスが終了すれば自動的に終了します。アプリを閉じた後にプロセスが残ることはありません。
3D モデリング・シミュレータ・散策モードでは、画面表示と入出力に共通のモデリング座標系を使っています。 既存の空力・構造・飛行計算は従来の解析内部座標を維持し、境界で安全に座標変換します。
| +X | 右方向（翼幅方向）— 軸コンパスでは赤 |
| −Y | 湖の方向（飛行方向）— 軸コンパスでは青が +Y |
| +Z | 上方向（高さ）— 軸コンパスでは緑 |
散策モードで画面右下に出る座標とコピーされる modeling 座標も同じ座標系です。