揚力を計算して滑空する

概要

今モックで作成しているVRコンテンツに、滑空の要素を入れたかったので揚力の計算について調べてみたのでそのメモです。

ちなみに適当に飛行機っぽい形状を作って適用した動画はこんな感じです↓

Unityでグライダー風のものを作ってみた。 https://t.co/M2GwICz84g
— edom18@AR / MESON (@edo_m18) 2019年1月11日

また、実装にあたって参考にさせてもらったのは以下の記事です。

６４章：二次元翼の揚力と抗力

www.cfijapan.com

blog.goo.ne.jp

揚力の計算

揚力（Lift）は以下の数式から求めることができるようです。

ここで、 $C_L$ は無次元の係数で揚力係数、 $C_D$ は無次元の係数で抗力係数といいます。
また、数式の記号の意味は以下となります。

$ρ$ は空気密度、 $q$ は流速、 $A$ は翼の面積です。

それぞれの単位は以下のようになります。

流体の密度 $ρ$ ですが、場所によって異なるようで、海面高度の大気中は大体[tex: 1.2250kg/m³]となるようです。
また、揚力係数、抗力係数は実験的に求められているようで、こちらの記事から画像を引用させていただくと以下のようなグラフになるようです。

力のかかる方向

揚力、抗力については上で書いた通りです。
この力のかかる方向を図示すると以下のようになります。

f:id:edo_m18:20190113152919j:plain

$\theta$ が進行方向と流体（空気）との成す角度です。
そして $L$ が揚力（Lift）を、 $D$ が抗力（Drag）を、 $mg$ が重力方向を表しています。

この図を見てもらうと分かりますが、揚力は進行方向に対して垂直、抗力は進行方向と平行（ただし逆向き）となります。

揚力のかかる方向を求める

色々な記事を見ても、揚力のかかる方向が垂直であることは示されているものの、「じゃあ実際プログラムするときに垂直方向ってどっちさ？」となり、色々考えた結果、以下のようにして求めるようにしたところ、それっぽく動いているのでこれを採用しています。

実際は流体力学などから、流体、そして渦の生成、気圧など様々な条件から方向が定まり、実際に揚力のかかる方向を計算するのだと思いますが、今回はあくまで「それっぽく」動くことが目的だったので単純に垂直な方向を求めています。

考え方としては以下のような感じです。

進行方向との垂直方向、つまり外積の向く向きを採用する
しかし「飛行機が逆さま」になっていることも考慮すると、もうひとつの軸の取り方で結果が変わってしまう
平行に近い角度で進行している場合は、翼の「右」方向と進行方向との外積が「上」として都合が良さそう（進行方向との外積なので必ず進行方向に対して垂直となる）
反対の考慮は？
翼の「上」ベクトル（leftWing.transform.up）と進行方向との外積方向をまず求める
通常飛行（翼の前ベクトルと平行的な方向）の場合においては、翼の右ベクトル（leftWing.transform.rightと(5)のベクトルは概ね同じ方向を向く
翼に対して後ろ方向に進行している場合は右ベクトルとは反対方向を向くようになる
(5)のベクトルと右ベクトルとの「内積」を計算し、マイナスとなる場合は逆と判断
-leftWing.transform.rightを計算のベクトルに用いる

というようなフローで解決するようにしてみました。

実際に、翼の上ベクトルと進行方向との外積で求めたベクトルの動きを動画にしてみると以下のような感じになります。

進行方向の判別（揚力の適用向き用）

黄色のバーが進行方向を表し、紫色のバーが外積によって求めたベクトルです。
動画を見てもらうと分かりますが、進行方向が翼に対して前方に向いている場合は、翼の右ベクトルとの角度は鋭角になっており、進行方向が逆転した場合に鈍角になることが分かります。

このことから、翼の右ベクトルと求めたベクトルとの内積を取り、結果がマイナス（＝鈍角）なら進行方向が想定と逆、とみなすことができるわけです。

コードにすると以下のような感じです。

Vector3 dir = v.normalized;

// 進行方向に対しての「右」を識別する
Vector3 lcheckDir = Vector3.Cross(_leftWing.up, dir);

float lcheckd = Vector3.Dot(_leftWing.right, lcheckDir);

Vector3 lright = (lcheckd < 0) ? -_leftWing.right : _leftWing.right;

// ldirが揚力の働く方向（Lift Direction）
Vector3 ldir = Vector3.Cross(dir, lright);

コード全文

そんなに長くないのでコード全文を載せておこうと思います。

ちなみに「揚力係数」と「抗力係数」についてはアニメーションカーブを用いて、概ね以下のような感じで設定してそこから値を取得するようにしています。

揚力係数のアニメーションカーブ（グラフ）

f:id:edo_m18:20190113172042p:plain

抗力係数のアニメーションカーブ（グラフ）

f:id:edo_m18:20190113172125p:plain

using System.Collections;
using System.Collections.Generic;
using UnityEngine;

public class LiftTest : MonoBehaviour
{
    [SerializeField]
    private Transform _body;
    
    [SerializeField]
    private Transform _leftWing;

    [SerializeField]
    private Transform _rightWing;

    [SerializeField]
    private AnimationCurve _liftCoeff;

    [SerializeField]
    private AnimationCurve _dragCoeff;

    [SerializeField]
    private ParticleSystem _lparticle;

    [SerializeField]
    private ParticleSystem _rparticle;

    [SerializeField]
    private float _pitchSpeed = 0.5f;

    [SerializeField]
    private float _rollSpeed = 0.5f;

    [SerializeField]
    private float _acc = 1000f;

    [SerializeField]
    private float _rho = 1.225f;

    [SerializeField]
    private float _area = 12f;

    [SerializeField]
    private float _initVelocity = 30f;

    [SerializeField]
    private bool _useLift = true;

    private Rigidbody _rigid;

    private void Start()
    {
        _rigid = GetComponent();
        _rigid.velocity = transform.forward * _initVelocity;

        _lparticle.Stop();
        _rparticle.Stop();
    }

    private void Update()
    {
        Control();
    }

    private void FixedUpdate()
    {
        if (_useLift)
        {
            CalcLift();
        }
    }

    private void Control()
    {
        if (Input.GetKeyDown(KeyCode.Space))
        {
            _lparticle.Play();
            _rparticle.Play();
        }

        if (Input.GetKeyUp(KeyCode.Space))
        {
            _lparticle.Stop();
            _rparticle.Stop();
        }

        if (Input.GetKey(KeyCode.Space))
        {
            _rigid.AddForce(_body.forward * _acc);
        }

        if (Input.GetKey(KeyCode.DownArrow))
        {
            _body.Rotate(_body.worldToLocalMatrix.MultiplyVector(_body.right), -_pitchSpeed);
        }

        if (Input.GetKey(KeyCode.UpArrow))
        {
            _body.Rotate(_body.worldToLocalMatrix.MultiplyVector(_body.right), _pitchSpeed);
        }

        if (Input.GetKey(KeyCode.LeftArrow))
        {
            _body.Rotate(_body.worldToLocalMatrix.MultiplyVector(_body.forward), _rollSpeed);
        }

        if (Input.GetKey(KeyCode.RightArrow))
        {
            _body.Rotate(_body.worldToLocalMatrix.MultiplyVector(_body.forward), -_rollSpeed);
        }
    }

    private void CalcLift()
    {
        Vector3 lpos = _leftWing.position;
        Vector3 rpos = _rightWing.position;

        Vector3 lup = _leftWing.up;
        Vector3 rup = _rightWing.up;

        Vector3 v = _rigid.velocity;

        float m = v.magnitude;
        float velocitySqr = m * m;

        Vector3 dir = v.normalized;

        // 揚力、抵抗ともに使う係数の計算（ρ/2 * q^2 * A）
        // ρ ... 密度
        // q ... 速度
        // A ... 面積
        float k = _rho / 2f * _area * velocitySqr;

        Debug.DrawLine(_body.position, _body.position + dir, Color.black);

        float ldot = Vector3.Dot(lup, dir);
        float lrad = Mathf.Acos(ldot);

        float rdot = Vector3.Dot(rup, dir);
        float rrad = Mathf.Acos(rdot);

        float langle = (lrad * Mathf.Rad2Deg) - 90f;
        float rangle = (rrad * Mathf.Rad2Deg) - 90f;

        float lcl = _liftCoeff.Evaluate(langle);
        float rcl = _liftCoeff.Evaluate(rangle);

        // 単位: N = kg・m/s^2
        float ll = lcl * k;
        float rl = rcl * k;

        // 進行方向に対しての「右」を識別する
        Vector3 lcheckDir = Vector3.Cross(_leftWing.up, dir);
        Vector3 rcheckDir = Vector3.Cross(_rightWing.up, dir);

        float lcheckd = Vector3.Dot(_leftWing.right, lcheckDir);
        float rcheckd = Vector3.Dot(_rightWing.right, rcheckDir);

        Vector3 lright = (lcheckd < 0) ? -_leftWing.right : _leftWing.right;
        Vector3 rright = (rcheckd < 0) ? -_rightWing.right : _rightWing.right;

        Vector3 ldir = Vector3.Cross(dir, lright);
        Vector3 rdir = Vector3.Cross(dir, rright);

        Vector3 lv = ldir * ll;
        Vector3 rv = rdir * rl;

        Debug.DrawLine(_leftWing.position, _leftWing.position + lv, Color.cyan);
        Debug.DrawLine(_rightWing.position, _rightWing.position + rv, Color.cyan);

        float lcd = _dragCoeff.Evaluate(langle);
        float rcd = _dragCoeff.Evaluate(rangle);

        float ldrag = lcd * k;
        float rdrag = rcd * k;

        Vector3 drag = -dir * (ldrag + rdrag);

        Debug.DrawLine(_body.position, _body.position + drag);

        Vector3 force = (lv + rv + drag) * Time.deltaTime;;

        _rigid.AddForce(force);
    }
}

e.blog

主にUnity/UE周りのことについてまとめていきます

揚力を計算して滑空する

概要

揚力の計算

力のかかる方向

揚力のかかる方向を求める

コード全文