디스크 골프 [#1] - 프리스비 운동 R&D

디스크 골프는 플라스틱 원반을 던져서 바스켓에 넣는 게임이다. 규칙은 골프와 동일하며 구현할 것 자체는 많이 없지만 핵심이 되는 프리스비 운동이 곧 디스크 골프 게임의 퀄리티에 직결하므로 프리스비 운동을 신경써서 R&D하기로 했다.

 

가장 먼저 해본일은 프리스비가 어떻게 움직이는지 알아보는 것이었다. 유튜브에서 디스크 골프 대회 영상을 참고하였다.

https://youtu.be/Ldgt7ZYxTcU?si=AK4Oq9DkA-AZsfSH

 

 

프리스비는 수평에서 살짝 각을 줘서 날리면 다음과 같이 곡선을 그리며 상승하다가 낙하한다. 

이걸 보고 든 생각은 이걸 물리학을 적용해서 구현해야 할지 아니면 곡선을 그리는 방정식을 활용해서 구현해야 할지 고민이 되었다. 아무래도 물리학을 사용해서 구현하는게 좀 더 보람이 있고 자연스럽게 작동할 것 같아서 가장 먼저 물리를 사용하는 방식을 R&D해보기로 하였다. 

 

관련 자료를 찾아봤는데 게임커뮤니티엔 프리스비의 운동에 관한 정보가 많이 없었다. 그래서 일단은 프리스비 비행의 기본이 되는 공기역학에 대해서 공부하기로 하였다. 이 부분에 관해선 유니티는 프리스비의 운동에 대해선 자료가 많이 없지만 항공기 시뮬레이터에 대해선 자료가 많았으므로 참고가 많이 되었다.

 

https://youtu.be/p3jDJ9FtTyM?si=To0o9KwxXr9IOGkr 

 

비행하는 물체에는 4가지 힘이 작용한다 

 

 

Lift (양력)

Drag (항력)
Gravity (중력)
Thrust (추력)

 

의 4가지로

 

Thrust는 기계적인 힘이므로 제외하고 Gravity는 고정된 값이기 때문에 프리스비의 운동에서 눈여겨 봐야할 건 Drag와 Lift힘이다.

 

공식 유도를 제외하고 간단하게 설명하면

 

Drag의 공식은 1/2 * 유체밀도 * 속도의 제곱 * 반경(프리스비에선 프리스비 면적) * Drag계수(a각에서의)

Lift의 공식은 1/2 * 유체밀도 * 속도의 제곱 * 반경(프리스비에선 프리스비 면적) * Lift계수(a각에서의)

 

유체밀도가 높아질수록, 속도가 빨라질수록, 물체의 반경이 커질수록, 수평과의 각도가 커질수록 Drag와 Lift의 값은 커진다.

 

https://web.mit.edu/womens-ult/www/smite/frisbee_physics.pdf

 

다음 자료의 자바코드를 참고해서 유니티의 C#코드로 바꾼다음 한번 적용해보았다.

 

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

public class FrisbeeThesis : MonoBehaviour
{
    private Rigidbody rb;
    // 중력가속도(m/s^2)
    private const double gravity = -9.81f;
    // 프리스비의 질량(kg)
    private const double mass = 0.175f;
    // 공기밀도 (kg/m^3)
    private const double rho = 1.23f;
    // 프리스비의 표면적 (πr^2)
    private const double area = 0.0568f;
    // 알파에서의 lift 계수
    private const double clo = 0.1f;
    // 알파와의 상대적인 lift 계수
    private const double cla = 1.4f;
    // 알파에서의 drag계수
    private const double cdo = 0.08f;
    // 알파와의 상대적인 drag계수
    private const double cda = 2.72f;
    private const double alpha0 = -4f;

    private float deltaT = 0.02f;
    void Start()
    {
        this.rb = GetComponent<Rigidbody>();
        this.rb.AddForce((Vector3.forward + Vector3.up).normalized * 1000f);
    }

    // Update is called once per frame
    private void FixedUpdate()
    {
        float alpha = Vector3.Angle(this.transform.forward,Vector3.forward);
        double cl = clo + cla * alpha * Math.PI / 180;
        double cd = cdo + cda * Math.Pow((alpha - alpha0) * Math.PI / 180, 2);
        double deltavy = (rho * Math.Pow(rb.velocity.z, 2) * area * cl / 2 / mass) * deltaT;
        double deltavx = -rho * Math.Pow(rb.velocity.z, 2) * area * cd * deltaT;
        rb.velocity = rb.velocity + new Vector3(0f, (float)deltavy, (float)deltavx);
    }
}

 

프리스비의 표면적은 유니티의 게임오브젝트를 사용해도 되지만 더 정확한 결과를 얻기 위해 자료의 표면적을 그대로 사용하였고, 나머지 계수들은 마찬가지로 사이트의 시뮬레이션 해놓은 값을 그대로 사용하였다. 월드 좌표의 forward와 프리스비의 로컬좌표를 각을 계산해서 alpha값으로 사용했으며, 프리스비에는 리지드바디를 부착해서 중력의 영향을 받게 해서 Lift는 상승힘만 작용하도록 바꾸었다.

 

 

 

코드를 적용하고 실험해보았는데 원하는 결과가 나오지 않았다. 이유를 알기 위해 계속 찾아보았는데 프리스비의 운동에는 2가지 요소가 작용한다는걸 알게 되었다. 첫번째는 공기역학적 힘인데 위 코드로 구현한 것이 공기역학적 힘이다.

 

두번째는 자이로스코프로이다. 위 코드에서 구현하지 못한 요소로 프리스비의 회전이 운동에 반영되지 않았기 때문에 다음과 같은 결과가 나온 것이다.

 

프리스비의 자이로스코프적인 힘을 구현하기 위해서 프리스비 시뮬레이터의 코드를 활용하기로 했다. 

 

https://github.com/tmcclintock/FrisPy/blob/main/frispy/disc.py

 

위 깃코드는 파이썬으로 되있어서 일단은 chat gpt를 사용하여 파이썬에서 C#으로 바꾼다음 사용하였다.

 

using System;

public class Model
{
    public float PL0 { get; set; } = 0.33f;
    public float PLa { get; set; } = 1.9f;
    public float PD0 { get; set; } = 0.18f;
    public float PDa { get; set; } = 0.69f;
    public float PTxwx { get; set; } = -0.013f;
    public float PTxwz { get; set; } = -0.0017f;
    public float PTy0 { get; set; } = -0.082f;
    public float PTya { get; set; } = 0.43f;
    public float PTywy { get; set; } = -0.014f;
    public float PTzwz { get; set; } = -0.000034f;
    public float Alpha0 { get; set; } = 4 * (float)Math.PI / 180;

    public float C_lift(float alpha)
    {
        return PL0 + PLa * alpha;
    }

    public float C_drag(float alpha)
    {
        return PD0 + PDa * (float)Math.Pow(alpha - Alpha0, 2);
    }

    public float C_x(float wx, float wz)
    {
        return PTxwx * wx + PTxwz * wz;
    }

    public float C_y(float alpha, float wy)
    {
        return PTy0 + PTywy * wy + PTya * alpha;
    }

    public float C_z(float wz)
    {
        return PTzwz * wz;
    }
}

 

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

public class EOM
{
    private float area;
    private float diameter;
    private float I_xx;
    private float I_zz;
    private float mass;
    private Model model;
    private float air_density;
    private float g;
    private float force_per_v2;
    private float torque_per_v2;
    private Vector3 F_grav;
    private Vector3 z_hat;

    public EOM(float area, float I_xx, float I_zz, float mass, float air_density = 1.225f, float g = 9.81f, Model model = null)
    {
        this.area = area;
        this.diameter = 2 * Mathf.Sqrt(area / Mathf.PI);
        this.I_xx = I_xx;
        this.I_zz = I_zz;
        this.mass = mass;
        this.model = model ?? new Model();
        this.air_density = air_density;
        this.g = g;

        this.force_per_v2 = 0.5f * air_density * area;
        this.torque_per_v2 = force_per_v2 * diameter;
        this.F_grav = new Vector3(0, 0, -mass * g);
        this.z_hat = new Vector3(0, 0, 1);
    }

    public static Matrix4x4 RotationMatrixFromPhiTheta(float phi, float theta)
    {
        float sp = Mathf.Sin(phi);
        float cp = Mathf.Cos(phi);
        float st = Mathf.Sin(theta);
        float ct = Mathf.Cos(theta);
        return RotationMatrix(sp, cp, st, ct);
    }

    private static Matrix4x4 RotationMatrix(float sp, float cp, float st, float ct)
    {
        Matrix4x4 matrix = new Matrix4x4();
        matrix.SetColumn(0, new Vector4(ct, 0, st, 0));
        matrix.SetColumn(1, new Vector4(sp * st, cp, -st * cp, 0));
        matrix.SetColumn(2, new Vector4(-st * cp, -sp * st, cp * ct, 0));
        matrix.SetColumn(3, new Vector4(0, 0, 0, 1));
        return matrix;
    }

    public static (float, float, float, float, float, Matrix4x4, float, Vector3) ComputeAngleOfAttack(float phi, float theta, Vector3 velocity, bool returnAllVariables = false)
    {
        float sp = Mathf.Sin(phi);
        float cp = Mathf.Cos(phi);
        float st = Mathf.Sin(theta);
        float ct = Mathf.Cos(theta);

        Matrix4x4 rotationMatrix = RotationMatrix(sp, cp, st, ct);

        Vector3 zhat = rotationMatrix.GetColumn(2);
        float vDotZhat = Vector3.Dot(velocity, zhat);
        Vector3 vInPlane = velocity - zhat * vDotZhat;
        float angleOfAttack = -Mathf.Atan2(vDotZhat, vInPlane.magnitude);

        if (returnAllVariables)
        {
            return (angleOfAttack, sp, cp, st, ct, rotationMatrix, vDotZhat, vInPlane);
        }
        else
        {
            return (angleOfAttack, 0, 0, 0, 0, Matrix4x4.identity, 0, Vector3.zero);
        }
    }

    public Dictionary<string, object> ComputeForcesAndTorques(float phi, float theta, Vector3 velocity, Vector3 angularVelocity)
    {
        Dictionary<string, object> result = new Dictionary<string, object>();

        var angleOfAttackResult = ComputeAngleOfAttack(phi, theta, velocity, true);
        float angleOfAttack = angleOfAttackResult.Item1;
        float sp = angleOfAttackResult.Item2;
        float cp = angleOfAttackResult.Item3;
        float st = angleOfAttackResult.Item4;
        float ct = angleOfAttackResult.Item5;
        Matrix4x4 rotationMatrix = angleOfAttackResult.Item6;
        float vDotZhat = angleOfAttackResult.Item7;
        Vector3 vInPlane = angleOfAttackResult.Item8;

        Vector3 zhat = rotationMatrix.GetColumn(2);
        Vector3 xhat = vInPlane.normalized;
        Vector3 yhat = Vector3.Cross(zhat, xhat);

        float forceAmplitude = force_per_v2 * velocity.sqrMagnitude;

        // Lift force
        Vector3 F_lift = model.C_lift(angleOfAttack) * forceAmplitude * Vector3.Cross(velocity.normalized, yhat);

        // Drag force
        Vector3 F_drag = model.C_drag(angleOfAttack) * forceAmplitude * (-velocity.normalized);

        // Gravitational force
        Vector3 F_grav = this.F_grav;

        // Total force
        Vector3 F_total = F_lift + F_drag + F_grav;

        // Acceleration
        Vector3 acceleration = F_total / mass;

        result.Add("F_lift", F_lift);
        result.Add("F_drag", F_drag);
        result.Add("F_grav", F_grav);
        result.Add("F_total", F_total);
        result.Add("Acceleration", acceleration);

        // Torque
        float torqueAmplitude = torque_per_v2 * velocity.sqrMagnitude;
        Vector3 w = new Vector3(angularVelocity.x * ct, angularVelocity.y, angularVelocity.x * st + angularVelocity.z);
        Vector3 T_x_lab = model.C_x(w.x, w.z) * torqueAmplitude * xhat;
        Vector3 T_y_lab = model.C_y(angleOfAttack, w.y) * torqueAmplitude * yhat;
        Vector3 T_z = model.C_z(w.z) * torqueAmplitude * z_hat;

        // Rotate into the disc frame
        Vector3 T_x = rotationMatrix * T_x_lab;
        Vector3 T_y = rotationMatrix * T_y_lab;

        // Total torque
        Vector3 T_total = T_x + T_y + T_z;

        result.Add("T_x_lab", T_x_lab);
        result.Add("T_y_lab", T_y_lab);
        result.Add("T_z", T_z);
        result.Add("T_x", T_x);
        result.Add("T_y", T_y);
        result.Add("T_total", T_total);

        return result;
    }

    public float[] ComputeDerivatives(float time, float[] coordinates)
    {
        float[] derivatives = new float[coordinates.Length];

        // If the disk hit the ground, then stop calculations
        if (coordinates[2] <= 0)
        {
            Array.Fill(derivatives, 0);
            return derivatives;
        }

        Vector3 velocity = new Vector3(coordinates[3], coordinates[4], coordinates[5]);
        Vector3 angularVelocity = new Vector3(coordinates[9], coordinates[10], coordinates[11]);

        var forcesAndTorques = ComputeForcesAndTorques(coordinates[6], coordinates[7], velocity, angularVelocity);

        // Copy velocity derivatives
        Array.Copy(velocity.ToArray(), 0, derivatives, 0, 3);

        // Copy acceleration derivatives
        Array.Copy(forcesAndTorques["Acceleration"] as Vector3, 0, derivatives, 3, 3);

        // Copy angular velocity derivatives
        Array.Copy(angularVelocity.ToArray(), 0, derivatives, 6, 3);

        // Copy torque derivatives
        Array.Copy((forcesAndTorques["T_total"] as Vector3).ToArray(), 0, derivatives, 9, 3);

        return derivatives;
    }
}

 

Model은 프리스비의 계수와 토크를 계산해서 리턴해주는 메서드가 있는 클래스이고

 

EOM은 Equation of Motion의 약자로 시뮬레이터를 돌리는데 필요한 힘을 직접 계산해주는 메서드가 있는 클래스이다.

 

Disc도 있지만 시뮬레이터를 돌리는 부분이라 일단은 생략하고 EOM의 ComputeForcesAndTorques부분부터 R&D해보기로 했다.