用Unity试着做翅膀的故事

引言

在学校制作项目中，终于有了制作飞行模拟器的机会，如果能为同样在挑战制作飞行模拟器的人提供参考，我会很高兴。

本文内容是“无限翼”的简单模型，尚未实现“有限翼”中由“翼端涡”引起的诱导阻力。

此外，代码并不完整，是用于学习的简易实现。

作为实现时的参考书，使用了约翰·D·安德森的《Fundamentals of Aerodynamics》第六版。

我没有正式学习过航空航天工程，所以如果有错误之处，恳请指出，我将感激涕零。

那么首先，我们来探讨一下解决问题的方法。

目的

Unity Engine 中虽然实现了物理引擎，但当然并没有准备空气动力学的物理运算或便利的函数😭

话虽如此，刚体的运动完全由作用在该刚体（Rigidbody）上的合力及其合力矩（扭矩）决定。最终只要计算出作用在机翼上的力及其力矩，并通过以下函数传递给 Unity 的物理引擎即可。这就是最终目标。

//力
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」。

将图表中显示的数据缩减为一组，使其更易于查看

这里简单说明一下雷诺数、马赫数和NCrit。

雷诺数是表示流体中惯性力与黏性力之比的数值。也就是说，雷诺数高的流动容易成为湍流，雷诺数低的流动容易成为层流。

马赫数是气流速度除以音速得到的比值（无量纲）。它表示压缩性的影响，越接近音速就越重要，但本文设想的是低马赫飞行，因此可以忽略。

NCrit是表示气流扰动的数值。可以这样理解：它较低时容易成为湍流，较高时则容易维持层流。

那么，我们来看看升力系数与迎角的图表。

表示升力系数与迎角关系的图表

首先，当迎角为0度时，升力系数大约为0.50。NACA4412 是具有弯度
（camber）的翼型，因此即使迎角为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保持最后参照的值，就暂且能实现失速。今后我们想利用面板法（panel method）等方式来表现更真实的失速，不过在实现之前打算先用这个办法应付过去。

至此，计算作用在机翼上的力所需的信息已经齐备，所以使用(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 引用的方法，是实时应用中的标准做法。

不过，本实现是“无限翼”的二维翼型模型，并未包含有限翼因翼尖涡产生的诱导阻力。此外，失速后的行为也是简化处理，无法再现实际机翼复杂的分离现象。如果参数设置出错，就会出现极其离谱的行为。

今后计划推进基于 BET（Blade Element Theory，叶素理论）的展向积分、诱导阻力计算，以及操纵面（副翼、升降舵、方向舵）的实现。如果实现这些，我认为将更接近更真实的飞行模拟器。

我也还在学习空气动力学，如有错误，欢迎指出。希望本文能对以自制飞行模拟器为目标的人有所帮助。