Unity 3D中实现XPBD绳索动力学模拟遇问题求助

我尝试在Unity中基于XPBD（扩展基于位置的动力学）实现软体模拟，参考了Ten Minute Physics频道Matthias Muller的教程，将原代码转换为Unity风格（比如使用Vector3对象），但无论如何调整参数，模拟要么出现爆炸效果，要么无法对抗重力维持形状。即便简化为仅使用距离约束的单条绳索，仍未达到预期效果。

以下是仅使用距离约束的脚本代码（挂载在世界原点的空GameObject上）：

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

public class RopeDynamics : MonoBehaviour
{
    public MeshFilter meshFilter;

    public bool showIndexes = true;
    public int indexLabelSize = 16;

    public int subSteps = 10;
    public float edgeCompliance = 100f;
    public float mass = 1f;


    public Vector3[] nodes;
    public int[] edgeNodeIds;

    [SerializeField]
    private float[] _restEdgeLengths;    
    [SerializeField]
    private float[] _inverseMasses;
    [SerializeField]
    private Vector3[] _velocities;
    [SerializeField]
    private Vector3[] _prevPositions;

    private Vector3 _gravity;
    private Vector3[] _gradients = new Vector3[4];

    public int NodeCount
    {
        get
        {
            return nodes.Length;
        }
    }

    public int EdgeCount
    {
        get
        {
            return edgeNodeIds.Length / 2;
        }
    }

    private void Awake()
    {
        //Create the Rope Nodes/Vertices.
        nodes = new Vector3[]
        {
            new Vector3(0f,2f,0f),
            new Vector3(0f,1.75f,0f),
            new Vector3(0f,1.5f,0f),
            new Vector3(0f,1.25f,0f),
            new Vector3(0f,1f,0f),
            new Vector3(0f,0.75f,0f),
            new Vector3(0f,0.5f,0f)
        };

        edgeNodeIds = new int[]
        {
            0,1,
            1,2,
            2,3,
            3,4,
            4,5,
            5,6
        };

        var mesh = new Mesh();
        mesh.vertices = nodes;
        mesh.SetIndices(edgeNodeIds, MeshTopology.Lines, 0);
        mesh.RecalculateBounds();
        meshFilter.mesh = mesh;

        //Get Local Vector of Gravity.
        _gravity = transform.InverseTransformDirection(Physics.gravity);

        InitializeSimulation();
    }


    private void FixedUpdate()
    {
        //Get the Substep.
        var sdt = Time.fixedDeltaTime / subSteps;

        for (int i = 0; i < subSteps; i++)
        {

            PreSolve(sdt, _gravity);
            Solve(sdt);
            PostSolve(sdt);
        }
    }


    private void InitializeSimulation()
    {
        //Initialize the Arrays.
        _restEdgeLengths = new float[EdgeCount];
        //_restVolumes = new float[tetMesh.TetrahedronCount];
        _inverseMasses = new float[NodeCount];
        _velocities = new Vector3[NodeCount];
        _prevPositions = new Vector3[NodeCount];

        //Get all of the Rest Lengths for each Edge.
        for (int i = 0; i < EdgeCount; i++)
        {
            //_restEdgeLengths[i] = tetMesh.GetSqrEdgeLength(i);
            _restEdgeLengths[i] = GetMagEdgeLength(i);
        }

        var w = 1f / (mass / NodeCount);

        //Get all of the Inverse Masses.
        for (int i = 0; i < NodeCount; i++)
        {
            _inverseMasses[i] = w;            
        }

        _inverseMasses[0] = 0f;
    }

    private void PreSolve(float dTime, Vector3 gravity)
    {
        //For each node.
        for (int i = 0; i < NodeCount; i++)
        {
            //Skip if the Inverse Mass is Zero/Infinite.
            if (_inverseMasses[i] == 0f)
                continue;

            //Get selected Velocity and add the Gravity vector.
            _velocities[i] += gravity * dTime;

            //Cache Previous Position.
            _prevPositions[i] = nodes[i];

            //Add Velocity vector to the Nodes Position vector.
            nodes[i] += _velocities[i] * dTime;
        }

    }

    private void Solve(float dTime)
    {
        SolveEdges(dTime, edgeCompliance);        
    }

    private void PostSolve(float dTime)
    {
        for (int i = 0; i < NodeCount; i++)
        {
            //Skip if the Inverse Mass is Zero/Infinite.
            if (_inverseMasses[i] == 0f)
                continue;

            //Update the selected Velocity.
            _velocities[i] = (nodes[i] - _prevPositions[i]) / dTime;
        }

        //Update the Mesh.
        UpdateMesh();
    }

    private void SolveEdges(float dTime, float compliance)
    {
        //Calculate the alpha for stiffness.
        var alpha = compliance / (dTime * dTime);

        //For each Edge.
        for (int i = 0; i < EdgeCount; i++)
        {
            //Get the Node Id's for the selected Edge.
            var ids = GetEdgeNodeIds(i);
            var eNodes = GetEdgeNodes(i);

            //Get the Inverse Masses for each Node.
            var w0 = _inverseMasses[ids[0]];
            var w1 = _inverseMasses[ids[1]];

            var w = w0 + w1;

            if (w == 0f)
                continue;

            //Get the current length of the edge.
            //var len = tetMesh.GetSqrEdgeLength(i);
            var len = GetMagEdgeLength(i);

            //Check if current length is Zero.
            if (len == 0f)
                continue;

            //Get the error between the current length and the rest length.
            var c = len - _restEdgeLengths[i];

            //Get the Vector difference between the two Nodes.
            var diff = eNodes[1] - eNodes[0];

            //Calculate the Gradient.
            _gradients[0] = diff / len;
            _gradients[1] = -_gradients[0];

            var gradMag0 = Vector3.Magnitude(_gradients[1]);
            var gradMag1 = Vector3.Magnitude(_gradients[0]);

            var gradSqr0 = gradMag0 * gradMag0;
            var gradSqr1 = gradMag1 * gradMag1;

            var s = -c / ((w0 * gradSqr0) + (w1 * gradSqr1) + alpha);

            var dX0 = (s * w0) * _gradients[1];
            var dX1 = (s * w1) * _gradients[0];
            

            nodes[ids[0]] += dX0;
            nodes[ids[1]] += dX1;
        }
    }

    private void UpdateMesh()
    {
        meshFilter.mesh.vertices = nodes;
    }

    public float GetMagEdgeLength(int index)
    {
        //Get the Nodes of the selected Edge.
        var node1 = nodes[edgeNodeIds[2 * index]];        
        var node2 = nodes[edgeNodeIds[2 * index + 1]];

        //Return the Squared Manitude of the 2 Points.
        return Vector3.Magnitude(node2 - node1);
    }

    public Vector3[] GetEdgeNodes(int index)
    {
        return new Vector3[2]
        {
            nodes[edgeNodeIds[2 * index]],
            nodes[edgeNodeIds[2 * index + 1]],
        };
    }

    public int[] GetEdgeNodeIds(int index)
    {
        return new int[2]
        {
            edgeNodeIds[2 * index],
            edgeNodeIds[2 * index + 1]
        };
    }

    //Used to draw spheres at the Nodes/Vertices.
    private void OnDrawGizmos()
    {
        Gizmos.color = Color.blue;
        GUI.contentColor = Color.blue;

        if (showIndexes)
        {

            for (int i = 0; i < NodeCount; i++)
            {
                var point = transform.TransformPoint(meshFilter.sharedMesh.vertices[i]);

                Gizmos.DrawSphere(point, 0.008f);

                Handles.Label(point, i.ToString(), new GUIStyle(GUI.skin.label) { fontSize = indexLabelSize });
            }
        }

    }

}

我始终无法定位问题所在，可能是某个细节错误被忽略，恳请提供帮助，非常感谢！

内容的提问来源于stack exchange，提问作者Xigtim