用Unity试着做翅膀的故事
Unityで翼を作ってみた話
はじめに
学校の制作でついにフライトシミュレータを作る機会が到来したので、同じフライトシムの制作に挑戦している方などの参考になれると嬉しいです。
本記事の内容は「無限翼」の簡単なモデルとなっており、「有限翼」の「翼端渦」による誘導抗力までは実装できていません。
また、コードは完全なものではなく、学習用の簡易実装です。
実装の参考書として「ジョン・D・アンダーソン」の「Fundamentals of Aerodynamics」第六版を利用しております。
自分は航空宇宙工学を正式に勉強したことがないので、もし間違っているところがありましたらご指摘いただけると泣いて喜びます。
それではまず問題へのアプローチについて考察していきます。
目的
Unity Engineには物理エンジンが実装されていますが、さすがに空気力学の物理演算や便利な関数なんて用意されていません😭
とはいっても、剛体の運動はその剛体（Rigidbody）に働く合力とその合力のモーメント（トルク）ですべて決まります。最終的には翼に働く力とその力のモーメントを算出し、以下の関数でUnityの物理エンジンに渡せばOKです。これが最終目標となります。
//力
rigidbody.AddRelativeForce(力);
//モーメント（トルクとモーメントは同じことだと思っていい）
rigidbody.AddRelativeTorque(モーメント);
さて、「どれくらいの力をどこに加えればいいの？」とまだいろいろ疑問があると思いますがひとまずそれを置いといて、実際の翼型について考えてみましょう。
ちょっとした空気力学
ざっくりとした説明になるので、詳しい数式や定義は「Fundamentals of Aerodynamics」の1章4節から参照してください。
いきなり三次元の翼ではなく、二次元の翼型（球形の牛あるある）にどういう力が働いているか分析してみましょう。
まずはこれから使う用語や座標系を紹介したいと思います。
用語と座標系
翼型のFBD
翼型の前縁と後縁を結んだ直線を翼弦線（Chord）と呼び、その長さを翼弦長（Chord length）と呼びます。
上記の図の青い線は翼弦線を延長した線で、気流の方向と翼弦線のなす角度を迎角（Angle of Attack, AoA）と呼びます。
翼に働く合力の気流に対して垂直な成分を揚力(Lift)、平行な成分を抗力(Drag)と定義します。
次に、気流が翼型に与える力について考えてみます。
翼型に作用する空気力
固定された翼型に対して左から右の方向に一様な空気を送り込むと下記の画像のように空気の分子が翼に衝突し、翼型に圧力と摩擦応力をかけます。
圧力と摩擦応力は面積単位あたりの力（Nm^{-2} = {Pa}）で力（N）ではないのでお気をつけください
翼型の一点一点を見ると圧力による力は曲線に対して垂直、摩擦応力による力は曲線に対して平行になっていることがわかります。
それ以外の力はないのでこれらを足せば（積分）合力と合力モーメントが計算できるのでしょう。
ならば早速実装してみましょう！
…ん？結局何を足せばいいの？…😅
そうです。残念ながら我々には翼型の一点一点にかかっている圧力の情報がありません。
なので積分する対象も💩もないです。
ここで「空気の分子をパーティクルとしてシミュレーションしたり、ナビエ–ストークス方程式を解いて流体シミュレーションを行えば圧力と摩擦応力の分布を算出できるのでは？」と思っている方もいるかもしれませんが、NASAのスーパーコンピュータがお持ちでしたらどうぞやってみてください。
余談ですが、実はMSFS2024にはCFD（Computational Fluid Dynamics/数値流体力学）によるリアルタイムシミュレーションが実装されていますし、CFDはGPUに向いていますが、この記事の範囲外です。
今からCFDを勉強する時間もないし使える風洞もあるわけがない…
では、我々にはもう手段はないでしょうか🥺…
…と思うかもしれませんがまだ諦めないでください。なぜなら世の中の研究者が風洞実験やXFOILによる翼型データを公開してくれています！
次のセクションで翼型データの読み方、使い方について解説していきます。
翼型データと空力係数
それでは翼型のNACA4412をググってみましょう。
NACA4412の検索結果
検索結果の上位には、おそらく「Airfoil Tools」というサイトが表示されているかと思います。
サイトを開いて、上から下の順番に何が書かれているか解説していきます。
閲覧中の翼型の図と型が表示されています
ページの一番上に検索した翼型の図と型が表示されています。
肝心な翼型のデータがここにあります
次に「Similar airfoils」という今閲覧している翼型に近い翼型のリストが表示されています。似た翼型のデータを比較したいときに便利ですが、今回は使いません。その下が我々の求めている翼型のデータです。
いろいろな係数と迎角の関係を表すグラフ集
サイトの一番下には上にあるチェックマークの入れたデータのグラフが集まっています。
…c_l、c_d、c_m、α？なんの変数？今まで話してきた圧力分布は🤨？
と思いますよね。これから説明していきますのでご安心ください。
空力係数について
注意：ここからすこし数学が重くなります。
まず、これからの数式によくでてくる空気の動圧 q_{\infty}を定義します。
\begin{equation} q_{\infty} ≡ \dfrac{1}{2}\rho_{\infty}V_{\infty}^{2} \end{equation}
ここでは\rho_{\infty}が空気の密度でV_{\infty}が気流の速度です。要するに空気の動圧は単位体積あたりの空気の運動エネルギーです。もちろん気流は場所によって速度や密度が変わったりしますのであくまで翼型の遥か前方にある空気の話です。またここの速度は翼型に対する相対速度になります。
次に、いよいよ空力係数について解説していきます。
\begin{equation} c_l ≡ \dfrac{L'}{q_{\infty}c} \end{equation}
c_lは翼型の揚力係数といいます。分母のcは翼弦長で、L'は翼型に働く単位翼幅あたりの揚力です。この式を用いてL'を求め、翼幅に掛けると実際の揚力を計算できます。
\begin{equation} c_d ≡ \dfrac{D'}{q_{\infty}c} \end{equation}
c_dは翼型の抗力係数で、D'は翼型に働く単位翼幅あたりの抗力です。揚力と同様、D'をこの式で求め、翼幅にかけると翼に働く抗力が得られます。
\begin{equation} c_m ≡ \dfrac{M'}{q_{\infty}c^{2}} \end{equation}
c_mはピッチングモーメント係数といいます。ここでM'が翼型に働く単位翼幅あたりのモーメントです。しかし力のモーメントは基準点に依存する値なので、「参照しているデータはどこを基準点にしているか」を判明しなければ使えません。
ほとんどの翼型は25%コード長あたりに空力中心がある
実は風洞試験データや本記事で使っている「Xfoil」のデータは25%コード長を基準点として整理されていることが主流です。理由は、ほとんどの翼型で空力中心が25%コード長付近にあるからです。空力中心まわりのモーメントは、理論上は迎角によらず一定です。しかし、我々が使うデータはあくまで25%コード長の点を基準にしていますので、若干変化します。空力中心について詳しく知りたい方は、「Fundamentals of Aerodynamics」の4章3節から参照してください。
最後に、上の２つの空力係数と同じく、(4)の数式でM'を求め、それを翼幅にかけると25%コード長まわりのモーメントが得られます。
また、空力係数の数式はバッキンガムのパイ定理による次元解析で導出できます。興味ある方は「Fundamentals of Aerodynamics」の1章7節から読んでみてください。
これらの空力係数はすべて無次元（単位がない）ため、翼型のサイズや気流の速度に依存しません。同じ形状・同じ迎角・同じレイノルズ数なら、どんなサイズ・速度でも同じ値になります。翼型データが、生の力のデータではなく係数として公開されている理由もこれです。
では、これらの係数はどこから来るのでしょうか？ 係数は、翼に働く圧力分布を積分して得られる数値です。つまり、誰かが一度風洞試験したり、CFDを使ったりして圧力分布を計算し、それを無次元化した結果が係数として公開されています。我々はその計算を自分で行う必要はなく、結果だけを受け取ればいいのです。
空力係数と迎角の関係
グラフをより見やすいようにレイノルズ数が100万、マッハ数が0、NCritが9になっているデータ以外からチェックマークを外し、テーブルの左下の「Update plots」をクリックします。
グラフに表示するデータを1セットに絞って見やすくします
ここでレイノルズ数、マッハ数、NCritについて軽く説明します。
レイノルズ数は流体において慣性力と粘性力の割合を表した数値です。つまりレイノルズ数の高い流れは乱流になりやすく、低い流れは層流になりやすいです。
マッハ数は気流の速度を音速で割った比（無次元）です。圧縮性の影響を表し、音速に近づくほど重要になりますが、本記事は低マッハの飛行を想定しているので無視して大丈夫です。
NCritは気流の乱れを表した数値です。低いと乱流になりやすく、高いと層流が維持しやすいという感覚で捉えていただければ大丈夫です。
では、揚力係数と迎角のグラフを見てみましょう。
揚力係数と迎角の関係を表すグラフ
まず、迎角が0度のとき、揚力係数が0.50くらいになっていますね。NACA4412はキャンバー
（反り）のある翼型なので、迎角が0のときでも揚力が発生しています。
迎角が16度付近を超えると急激に揚力係数が落ち、翼型が失速します。この角度を臨界迎角といいます。また-15度付近より低下すると同じく失速します。
抗力係数と迎角の関係を表すグラフ
抗力係数は迎角の絶対値が大きくなるにつれて急激に増加します。これは翼が気流に対して垂直に近づくほど、空気を押しのける面積が増えるからです。
ピッチングモーメント係数と迎角の関係を表すグラフ
空力中心まわりのモーメントは理論上一定ですが、XFOILが報告する値は25%コード長まわりであり、これは空力中心の近似位置です。そのため、実際のデータでは迎角に対して若干変化します。
これらのグラフの形を覚えておくと、実装後に翼の挙動がおかしいときのデバッグに役立ちます。次のセクションでは、これらのデータを実際にLUT（Look Up Table）として読み込み、実行時にどう使うかを見ていきます。
実装編
データのダウンロード
スクリプトを作っていく前にデータをダウンロードしましょう。
データの右端にある「Details」をクリックしてダウンロードページに移動します
データの右端にある「Details」をクリックすると…
このページの左の「Details」パネルの「」Daunload as CSV file」でデータをダウンロードできます
このようなページが開かれます。
左の「Details」のパネルにある「Download as CSV file」でCSVファイルとしてダウンロードします。
ファイルを開くと、このようになっているはずです。
ダウンロードしたCSVファイルの中身
ファイルの上にメタ情報が並んでおり、下に肝心なデータテーブルが入っていますね。
そして我々に必要なのは\alpha、c_l、c_d、c_mのデータだけです。
カンマで文字列を分けて構文解析していけばいいのでデータ読み込み用のスクリプトは省略します。
翼の物理演算
ではいよいよ翼にアタッチするスクリプトを作成していきます。
using UnityEngine;
//空力係数のラッパー
public struct AerodynamicCoefficients
{
//コンストラクター
public AerodynamicCoefficients(float cl, float cd, float cm)
{
this.cl = cl;
this.cd = cd;
this.cm = cm;
}
public float cl;
public float cd;
public float cm;
}
public class SimpleWing : MonoBehaviour
{
}
まずは「SimpleWing」というUnityのMonoBehaviourから継承するクラスを作ります。
ついでに、データテーブルを参照して空力係数を調べる度に、調べた係数を戻り値として返せるよう、「AerodynamicCoefficients」という構造体を定義します。
次に、翼の計算に必要な基本情報を追加します。
public class SimpleWing : MonoBehaviour
{
//翼弦長(m)
public float chord;
//翼幅(m)
public float span;
//25%コード長の相対位置（m）
private Vector3 quarter_chord_pos;
//リジッドボディコンポーネント
[SerializeField]
private Rigidbody rb;
//重心の相対位置（m）
private Vector3 com;
//海面レベルの空気密度(kg・m^-3)
private const float AIR_DENSITY_SEA_LEVEL = 1.225f;
//ダウンロードしたデータに指すパス、Unityのプロジェクトパスに相対する("./Assets/.../xf-naca4412-il-1000000.csv")
[SerializeField]
private string xfoil_file_path;
//読み込んだデータテーブル
//ダウンロードしたデータを格納するだけのクラスなので実装は省略します
private AirfoilDataTable data_table;
}
Start()の関数でCSVファイルからデータを読み込み、25%コード長の位置と重心の位置を取得します。
private void Start()
{
//data_tableの初期化（データをdata_tableに格納する）
FillDataTable(xfoil_file_path);
//アタッチするオブジェクトのピボット位置を翼の前縁とします
quarter_chord_pos = new Vector3(0, 0, -chord / 4);
com = rb.worldCenterOfMass;
}
UnityのFixedUpdate()は物理演算システムと同期して固定された時間間隔で呼び出されるので、翼まわりの力の計算はこれを使って行います。
private void FixedUpdate()
{
//翼に働く力を計算し、Rigidbodyに加える関数をここで呼び出す
ComputeForcesAndMoments();
}
そうしたらこの関数の中身を作っていきます。
private void ComputeForcesAndMoments()
{
//25%コード長の位置の速度を求めます
//rigidbodyの回転による角速度も考慮する必要があります
Vector3 v_quarter_chord_pos = rb.linearVelocity + Vector3.Cross(rb.angularVelocity, transform.TransformPoint(quarter_chord_pos) - com);
//そうしたら、計算した速度を使って、空気の動圧を算出します
float q_fs = ComputeDynamicPressure(v_quarter_chord_pos);
}
このComputeDynamicPressureという関数は(1)の数式を使って計算しているだけです。とりあえず高度による密度変化を考えず、海面レベルの空気密度を使っていきます。もし翼の高度から空気密度を計算しているなら、ここのAIR_DENSITY_SEA_LEVELを計算した結果に置き換えるいいとでしょう。
private float ComputeDynamicPressure(Vector3 v_fs)
{
return 0.5f * AIR_DENSITY_SEA_LEVEL * v_fs.magnitude * v_fs.magnitude;
}
次に、迎角である\alphaを計算し、データテーブルをサンプリングして空力係数を出します。
void ComputeForcesAndMoments(){
//...
//速度がゼロだったらなにもしない
if (v_quarter_chord_pos.sqrMagnitude < 0.0001f) return;
Vector3 velocity_norm = v_quarter_chord_pos.normalized;
//-1をかけると上に傾けたときの迎角がプラス、下に傾けてるときの迎角がマイナスになり、データテーブルの座標系に合わせます
float up_component = -1f * Vector3.Dot(transform.up, velocity_norm);
float forward_component = Vector3.Dot(transform.forward, velocity_norm);
float alpha = Mathf.Atan2(up_component, forward_component) * 180f / Mathf.PI;
//データテーブルをサンプリング
AerodynamicCoefficients coefficients = LookUpTable(alpha);
}
もし\alphaがデーターテーブルの範囲外だった場合、c_lとc_mを0にして、c_dを最後に参照した値のままにすればとりあえず失速させることができます。今後はパネル法などを利用して、より現実的な失速を演出させたいと考えていますがそれまではこれでしのぎたいと思います。
これで翼に働く力を計算するための情報が揃ったので(2)、(3)、(4)の式を使って実際の揚力、抗力、25%コード長まわりのモーメントを計算します。
void ComputeForcesAndMoments(){
//...
float lift_magnitude = q_fs * coefficients.cl * chord * span;
float drag_magnitude = q_fs * coefficients.cd * chord * span;
float moment_magnitude = q_fs * coefficients.cm * chord * chord * span;
Vector3 lift = lift_magnitude * (transform.up - velocity_norm * Vector3.Dot(velocity_norm, transform.up)).normalized;
Vector3 drag = -1 * drag_magnitude * velocity_norm;
Vector3 net_force = lift + drag;
}
UnityのAddRelativeForce()、AddRelativeTorque()を使って、力とモーメントを与えていくので、計算した力とモーメントをRigidbodyコンポーネントのローカル座標に変換します。
void ComputeForcesAndMoments(){
//...
//合力の座標変換
Vector3 net_force_rb =
Vector3.Dot(net_force, rb.transform.forward) * Vector3.forward +
Vector3.Dot(net_force, rb.transform.up) * Vector3.up +
Vector3.Dot(net_force, rb.transform.right) * Vector3.right;
//25%コード長の位置から重心までの変位
Vector3 quarter_chord_pos_world = transform.position - transform.forward * (0.25f * chord);
Vector3 quarter_chord_pos_to_com = com - quarter_chord_pos_world;
//その変位をRigidbodyのローカル空間に座標変換します
Vector3 quarter_chord_pos_to_com_rb =
Vector3.Dot(quarter_chord_pos_to_com, rb.transform.forward) * Vector3.forward +
Vector3.Dot(quarter_chord_pos_to_com, rb.transform.up) * Vector3.up +
Vector3.Dot(quarter_chord_pos_to_com, rb.transform.right) * Vector3.right;
//重心まわりのモーメントを算出
//空力の符号規約では正のcm = 機首上げ。
//Unityでは+X方向のトルクが機首下げなので、Vector3.left（-X）で符号を合わせる。
Vector3 moment_com = (moment_magnitude * Vector3.left) - Vector3.Cross(quarter_chord_pos_to_com_rb, net_force_rb);
}
ここまで読んでくれた方、お疲れ様です。計算した結果をAddRelativeForce、AddRelativeTorqueに渡せば完成です。
void ComputeForcesAndMoments(){
//...
rb.AddRelativeForce(net_force_rb);
rb.AddRelativeTorque(moment_com);
}
挙動の確認
いよいよ翼を飛ばすところまで来たので、翼のゲームオブジェクトを作っていきます。
作成したスクリプトとRigidbodyをゲームオブジェクトにアタッチし、Rigidbodyの質量を8kg、翼の翼幅、翼弦長をそれぞれ2mと1mに設定します。そして、見た目として翼の大きさにスケーリングされた赤いキューブのメッシュを子どもとしてつけると完成です。
では適当に高いところから落としてみましょう…
滑空する翼
なんと美しい🥹…
NACA4412のデータによりますとピッチングモーメント係数が常にマイナスで、予想通り、機首下げの方向に曲がりながら飛んでいきましたね。
次に、翼をプレハブ化して、飛行機作ってみましょう。
翼のプレハブとカプセルでできた飛行機
どっかで見たことある爆弾飛行機ですね。尾翼はキャンバーのない翼型が望ましいのでNACA0006のデータを使っております。また水平尾翼を若干下向きに傾けさせています。
では適当に高い場所から落としてみましょう…
空を飛ぶ飛行機
素晴らしい😭…
また、翼のスクリプトはそのついているオブジェクトの向きから迎角を計算していますので、そのオブジェクトを何らかの方法で制御できるようにすれば飛行機の操縦ができます。
結論
ここまで読んでいただきありがとうございます。
本記事では、XFOILの空力データを使って翼に働く揚力・抗力・モーメントを計算し、Unityの物理エンジンに渡すことで、滑空する翼と簡単な飛行機を実装しました。生の圧力分布を積分する代わりに、無次元の空力係数をLUTとして参照するこの方法は、リアルタイムアプリケーションでは標準的なアプローチです。
ただし、本実装は「無限翼」の2D翼型モデルであり、有限翼の翼端渦による誘導抗力は含まれていません。また、失速後の挙動も簡易的なもので、実際の翼の複雑なはく離現象は再現できていません。パラメーターの設定をミスすると、とんでもない挙動になります。
今後は、BET（Blade Element Theory）によるスパン方向の積分、誘導抗力の計算、制御翼（エルロン、エレベーター、ラダー）の実装を進める予定です。これらを実装すれば、より現実的なフライトシミュレータに近づくと思います。
僕もまだ空力学を勉強中の身ですので、間違いがあれば指摘いただけると嬉しいです。本記事が、自作フライトシミュレータを目指す方の一助になれば幸いです。
来源：zenn · 游戏开发社区 · zenn.dev