揚力 Lift

揚力の量は、次の計算式で計算できます。
(1/2) * p * ( V * V ) * S * Cl 

• p 空気密度 Pressure / Air density [kg/m^3]
• V 速度 Velocity [m/s]
• S 翼面積 Surface (size) [m^2]
• C 揚力係数 Coefficient (of Lift) [-]
  
  

つまり「0.5 × 空気密度 × (速度 × 速度) × 翼面積 × 揚力係数 」で計算できるということですね。

public float CalculateLift(float coefficient, float surfaceSize, float velocity, float airDensity)

{

    return 0.5f * airDensity * (velocity * velocity) * surfaceSize * coefficient;

}

どのパラメーターも意外と簡単です。順番に見ていきましょう。

上の揚力曲線カーブのX軸は迎角になっていますが、この画像の通り『気流』と『翼弦線』との間にできる角度が迎角です。空間が無風の場合には、飛行機の移動ベクトルを逆にすると、それがそのまま飛行機に向かってくる気流になります。

ここでは話を簡単にするために無風を想定しているので、(-rigidbody.velocity)ベクトルに対する翼弦線のピッチ角を求めればよいです。

下の計算では、rigidBody.velocityを翼自身のローカル空間に取り込んでAtan2で角度を求めています。Atan2の結果の単位がラジアン(Radian)で返ってきますが、今回の実装では迎角は度(Degree)で扱っているのでMathf.Rad2Degをかけて度に変換しておきます。

// Convert to local vector ローカルなベクトルに変換

Vector3 localVelocity = this.transform.InverseTransformVector(rigidBody.velocity);

// Get Angle-of-attack(pitch) from local vector ローカルなベクトルから迎角(ピッチ角)を計算

float aoa = -Mathf.Atan2(localVelocity.y, localVelocity.z) * Mathf.Rad2Deg;

抗力 Drag

抗力は次の計算式で計算できます。

(1/2) * p * ( V * V ) * S * Cd 

これは、上で紹介した揚力の計算式と全く同じですね。
抗力係数が違うだけで、揚力のメソッドを再利用すればよいです。

揚力係数と比べると1桁小さな値です。
揚力が直線的な増加カーブだったのに対し、こちらはなんか・・・2乗っぽいカーブ?です。
特徴として、失速迎角あたりから急激に値が増えていきます。

コードまとめ

using UnityEngine;

[RequireComponent(typeof(Rigidbody))]

public class SimpleWing : MonoBehaviour

{

    // 翼面積 [m^2] Wing area 

    public float wingArea;

    void FixedUpdate()

    {

        Rigidbody rigidBody = this.GetComponent<Rigidbody>();

        // ローカルなベクトルに変換 Convert to local vector 

        Vector3 localVelocity = this.transform.InverseTransformVector(rigidBody.velocity);

        // 速度 Speed 

        float v = rigidBody.velocity.magnitude;

        // ローカルなベクトルから迎角(ピッチ角)を計算 Get Angle-of-attack(pitch) from local vector 

        float aoa = -Mathf.Atan2(localVelocity.y, localVelocity.z) * Mathf.Rad2Deg;

//------------------------

//　▼ 揚力

//------------------------

        // 揚力係数（簡易版） Lift coefficient 

        float cl = aoa * 0.1f;

        // 揚力を計算 Calculate lift 

        float lift = CalculateLiftOrDrag(cl, wingArea, v);

        // 揚力ベクトルを計算 Calculate lift direction (normalized) vector 

        Vector3 liftVector = Vector3.Cross(Vector3.Cross(rigidBody.velocity, this.transform.up), rigidBody.velocity).normalized;

//------------------------

//　▼ 抗力

//------------------------

        // 抗力係数（簡易版） Drag coefficient 

        float cd = Mathf.Clamp01(0.005f + Mathf.Pow(Mathf.Abs(aoa)*0.0315f,5));

        // 抗力を計算 Calculate drag

        float drag = CalculateLiftOrDrag(cd, wingArea, v);

        // 抗力ベクトルを計算 Calculate drag direction (normalized) vector

        Vector3 dragVector = -rigidBody.velocity.normalized; //ベロシティの逆ベクトルを正規化するだけ

//------------------------

//　▼ 力を適用

//------------------------

        // 揚力を適用 Apply lift

        rigidBody.AddForce(liftVector * lift);

        // 抗力を適用 Apply drag

        rigidBody.AddForce(dragVector * drag);

    }

    // 揚力・抗力を計算する共通メソッド

    public static float CalculateLiftOrDrag(float coefficient, float surfaceSize, float velocity, float airDensity = 1.293f)

    {

        float q = 0.5f * (velocity * velocity) * airDensity; //動圧

        return q * surfaceSize * coefficient;

    }

}