用Unity实现简单的绳子模拟（一）

用Unity实现简单的绳子模拟（一）

说到Unity物理，一般都会想到内置的physX物理引擎。其实我们也可以用Unity的API去实现一些简单的物理算法。
本文会介绍如何从头实现一个简单的绳子模拟的小Demo。

物理模型

物理模拟的第一步，就是要对真实世界的事物做简化，转化成可计算的模型。对于绳子，我们最直观的简化方法就是把它想象成一堆由弹簧连接的小球（一般也叫粒子）。这就是Mass-Spring模型。

用胡克定律我们就可以简单的求出弹簧力。
F = k * (L - L₀)
其中k是弹簧系数，L是当前弹簧的长度，L₀是弹簧不受力时的长度（初始长度）。

那么每个粒子的受力就很容易算出来了。

时间积分

有了粒子受到的外力之后，就需要开始考虑粒子如何跟着受力在时间的维度里运动了。
物理模拟的时间间隔通常是恒定的，我们这里就把时间间隔记为dt。
通常，只有初始时刻的参数是已知的，我们会根据初始时刻的参数（位置，速度，加速度）来推断下一时刻的位置。

图中的p代表粒子位置，v代表速度，a代表加速度，F代表受到的合力。
结合高中物理知识，上面的图应该比较容易理解。实际上，我们是在时间维度上做积分。这里我们简化地认为在dt足够小的时候，这段时间的运动可以认为是匀速运动。上图中的积分模式叫做Explicit Euler，这是一种误差大且不太稳定的时间积分方法，有空会仔细讨论不同的时间积分方法以及他们的误差分析。

实现

OK，现在就开始实现这样一个简单的绳子模拟器吧。

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

namespace PhysicsLab
{
    public class RopeSpringSolver : MonoBehaviour
    {
        public Transform ParticlePrefab;
        public int SubStepCount = 10;
        public int Count = 3;
        public int Space = 1;
        public float SpringK = 1.0f;
        public float AirResistanceRatio = 0.1f;
        [Range(0, 1)]
        public float Damping = 0.1f;

        public Vector3 ExternForce = Vector3.zero;

        private List<Transform> chain = new List<Transform>();
        private List<SpringParticle> particleList = new List<SpringParticle>();

        void Start()
        {
            for (int i = 0; i < Count; i++)
            {
                var obj = Instantiate(ParticlePrefab, transform, true);
                obj.Translate(0, -i * Space, 0);
                chain.Add(obj);

                // Construct Particles
                var particle = new SpringParticle();
                particle.invMass = 1;
                particle.radius = 0.5f * Space;
                particle.pos = new Vector3(0, -i * Space, 0);
                particle.velocity = Vector3.zero;
                particleList.Add(particle);
            }
        }

        void FixedUpdate()
        {
            float dt = Time.fixedDeltaTime / SubStepCount;

            // Update Particle Position
            // Root Particle follow Transform
            particleList[0].pos = transform.position;
            for (int n = 0; n<SubStepCount; n++)
            {
                for (int i = 1; i < Count; i++)
                {
                    var particle = particleList[i];

                    // Calculate Spring Force
                    // Last Particle
                    Vector3 forceDir = particleList[i - 1].pos - particle.pos;
                    Vector3 springForce = SpringK * forceDir.normalized * (forceDir.magnitude - Space);

                    // Next Particle
                    if (i < Count - 1)
                    {
                        forceDir = particleList[i + 1].pos - particle.pos;
                        springForce += SpringK * forceDir.normalized * (forceDir.magnitude - Space);
                    }

                    // Update Particle Position according to Newton's 2nd Law
                    particle.pos += (1 - Damping) * particle.velocity * dt;

                    // Update velocity
                    Vector3 acceleration = (springForce + ExternForce - AirResistanceRatio * particle.velocity.magnitude * particle.velocity) * particle.invMass;
                    particle.velocity += acceleration * dt;
                }
            }

            // Apply Particle Position to Transform
            for (int i=0; i<Count; i++)
            {
                chain[i].position = particleList[i].pos;
            }
        }

        void OnDrawGizmos()
        {
            if (particleList == null || particleList.Count != Count) return;

            Gizmos.color = Color.blue;
            for (int i = 1; i < Count; i++)
            {
                var particleParent = particleList[i - 1];
                var particle = particleList[i];
                Debug.DrawLine(particleParent.pos, particle.pos);
            }
        }
    }

    public class SpringParticle
    {
        public float invMass; // 1 / mass
        public float radius;
        public Vector3 pos;
        public Vector3 velocity;
    }
}

Github地址：https://github.com/ossupero/UnityStrandSimulator/tree/master/Assets/Rope/SpringRope

讨论

如果亲自试验一把的话，相信很容易发现问题——绳子的移动好像非常缓慢。
这是因为默认设置把Damping和Air Resistance调得比较高。这两个参数的主要用途是在增加阻力，让速度不要变得太大。
那么，如果把这两个值调到0会怎么样？

绳子开始抽风，然后数值爆炸，出现NaN的错误提示。这就是Explicit Euler的不稳定性导致的结果。
所以弹簧模型在实际应用的时候，通常会用更稳定的Implicit Euler来做时间积分，但是这种积分方法实现更复杂，效率也更差。
另一种目前比较流行的物理建模方式是Position Based Dynamics，实现简单，效率高，会比这里的弹簧稳定很多：）