在Unity中基于Verlet算法模拟绳索
在Unity中运用Verlet算法加上胡克定律模拟绳索...
在Unity中基于Verlet算法模拟绳索
在Unity中基于Verlet算法模拟绳索
基本原理
Verlet 算法模拟运动
Verlet 算法推导
设某一时刻的点的位置函数为
\[y=r\left ( t \right )\]根据泰勒公式
\[f\left ( x \right )=\sum_{n=0}^{\infty}\frac{f\left ( x_{0} \right )}{n!}\cdot \left ( x-x_{0}\right )^{n}\]将位置函数二阶泰勒展开为
\[r\left ( t \right )=r\left ( t_{0} \right )+{r}'\left ( t_{0} \right )\cdot \left ( t-t_{0} \right )+\frac{\\{r}''\left ( t_{0} \right)\cdot \left ( t-t_{0} \right )^{2}}{2}+O\left ( t^{3} \right )\]令
\[t=t_{0}+\Delta t\]则①式：
\[\textcolor{blue}{r\left ( t_{0}+\Delta t \right )}=r\left ( t_{0} \right )+{r}'\left ( t_{0} \right )\cdot \Delta t+\frac{\\{r}''\left ( t_{0} \right)\cdot \left ( \Delta t \right )^{2}}{2}+O\left ( \Delta t^{3} \right )\]令
\[t=t_{0}-Delta t\]则②式：
\[\textcolor{red}{r\left ( t_{0}-\Delta t \right )}=r\left ( t_{0} \right )-{r}'\left ( t_{0} \right )\cdot \Delta t+\frac{\\{r}''\left ( t_{0} \right)\cdot \left ( \Delta t \right )^{2}}{2}+O\left ( \Delta t^{3} \right )\]①+②式得：
\[\textcolor{blue}{r\left ( t_{0}+\Delta t \right )}=2r\left ( t_{0} \right )-\textcolor{red}{r\left ( t_{0}-\Delta t \right )} +{r}''\left ( t_{0} \right)\cdot \left ( \Delta t \right )^{2}+O\left ( \Delta t^{3} \right )\]去掉泰勒余项进一步简化公式
\[\textcolor{blue}{r\left ( t_{0}+\Delta t \right )}=2r\left ( t_{0} \right )-\textcolor{red}{r\left ( t_{0}-\Delta t \right )} +{r}''\left ( t_{0} \right)\cdot \left ( \Delta t \right )^{2}\]处理$\textcolor{red}{{r}’'\left ( t_{0} \right)\cdot \left ( \Delta t \right )^{2}} $，一点点物理学小知识，位移对时间的一阶导数为此时刻的速度，速度对时间的一阶导数为此时刻的加速度
\[{r}'\left ( t_{0} \right)=v\left(t_{0}\right)\] \[v'\left(t_{0}\right)=a\left(t_{0}\right)\]即最后得到公式
\[\textcolor{blue}{r\left ( t_{0}+\Delta t \right )}=2r\left ( t_{0} \right )-\textcolor{red}{r\left ( t_{0}-\Delta t \right )} +{a}\left ( t_{0} \right)\cdot \left ( \Delta t \right )^{2}\]公式说明
由此可得到物体在$ t_{0}-\Delta t$、$ t_{0}$ 、$ t_{0}+\Delta t$ 这三个时刻的位置关系，将$ \Delta t$看作帧间隔,即得出结论：在已知物体的上一帧位置、此帧位置、物体的加速度 ，可以预测物体的下一帧位置。
重力模拟（$\textcolor{red}{a\left(t_{0}\right)}$项的处理）
在通常情况下，物体运动过程中所受合力为重力，即：我们忽略其他力，在只考虑重力的情况下。可将上式中的$\textcolor{red}{a\left(t_{0}\right)}$视为重力加速度。即可在物体运动模拟中添加重力的影响。
约束模拟
下面我们将点升维成线，即在点与点中考虑添加约束（应力）。 在上面的所有步骤中，我们完成对点运动和重力的模拟，因在这里只要考虑点与点之间应力对绳子内部各个点的影响即可。 笔者在这里初步设想： 将两个点之间看作一根弹簧，同时设定好约束长度。
当两点之间的距离超过约束长度时，同时为两点添加一个沿绳子方向向内的力
当两点之间的距离小于约束长度时，同时为两点添加一个沿绳子方向向外的力 基于此，运用胡克定律 $\textcolor{red}{F=k\cdot\Delta x}$,计算加速度$\textcolor{red}{a’}$
(注意：这里的$\textcolor{red}{\Delta x}$是指在这个时刻，当前点与连接点的距离，与Verlet算法中同一点上一时刻与下一时刻的位置变化量不同)
下一步单独计算应力的位移，方向为沿绳方向
\[x=a\cdot\left ( \Delta t \right )^{2}=\textcolor{red}{\frac{k\cdot \Delta x\cdot\left ( \Delta t \right )^{2}}{m}}\]除此之外，由于绳索上的点（除开端和结尾外）所受的应力有两部分，我们还需考虑上一个点对此点的应力和下一个点对此点的应力。
关于这一点，在编写程序中我们只需要创建一个类将他们储存起来即可。
由此，我们几乎完成了模拟绳索的所有步骤。
代码实现
创建类
首先创建一个Particle类，和一个Stick类，模拟点和两点之间绳子的约束
注意：通过设置Bool变量isLock来标记点是否被模拟算法影响
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
//创建粒子类
public class Particle
{
public Vector3 position;//当前位置
public Vector3 oldPosition;//上一帧位置
public bool isLock = false;//是否锁定
public Particle(Vector3 a)
{
position = a;
oldPosition = a;
}
}
//创建约束
public class Stick
{
public Particle particle_1;
public Particle particle_2;
public Stick(Particle a, Particle b)
{
particle_1 = a;
particle_2 = b;
}
}
初始化系统
因为归根结底要将算法应用于Unity，所以我将绳子的起始和结束分别关联到一个物体当中，设置好分段数目，然后通过计算自动分配点的位置
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
public GameObject startParticle;//起始位置
public GameObject endParticle;//结束位置
[SerializeField] private int particleCount = 10;//创建的粒子数目
private float stickLength;//约束的长度
private List<Particle> particles = new List<Particle>();
private List<Stick> sticks = new List<Stick>();
//初始化
private void InitSystem()
{
stickLength = (startParticle.transform.position - endParticle.transform.position).magnitude / (particleCount - 1);
//计算约束长度
Vector3 dir = (startParticle.transform.position - endParticle.transform.position).normalized;//开始点与结尾点的方向向量
for (int i = 0; i < particleCount; i++)
{
particles.Add(new Particle(startParticle.transform.position - dir * stickLength * i));
}
for (int i = 0; i < particles.Count - 1; i++)
{
sticks.Add(new Stick(particles[i], particles[i + 1]));
}
}
编写模拟算法
根据上边的结论：
Verlet 算法模拟
1
2
3
4
5
6
7
8
9
10
11
12
13
//Verlet 算法模拟
var verletDel = new Vector3();
for (int i = 0; i < particles.Count; i++)
{
if (particles[i].isLock == false)
{
Vector3 temp = particles[i].position;
verletDel = particles[i].position - particles[i].oldPosition + new Vector3(0, -gravity, 0) * Time.fixedDeltaTime * Time.fixedDeltaTime;//Verlet 算法算法下的位置增量
particles[i].position += verletDel;
particles[i].oldPosition = temp;//更新oldPosition
}
}
约束模拟
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
//约束
foreach (var k in sticks)
{
Vector3 dir = (k.particle_2.position - k.particle_1.position).normalized;//绳子方向的单位向量
float curLength = (k.particle_2.position - k.particle_1.position).magnitude;//当前节点长度
float deltaLength = curLength - stickLength;
var xDelta = ki * deltaLength * Time.fixedDeltaTime * Time.fixedDeltaTime;
if (!k.particle_1.isLock)
{
k.particle_1.position += xDelta * dir;
}
if (!k.particle_2.isLock)
{
k.particle_2.position -= xDelta * dir;
}
}
但是此时的约束模拟算法经实际验证还存在一些问题
我们遍历整个List<Sticks>时，是从顶端到末尾依次遍历，然后计算得出点的位置，一个一个的移动。
而实际上，绳子上每个点的移动，都会影响除这个点外的所有点的位置。
具体一点来讲，就像一串珠子，单独对某一个特定的珠子来看，影响其位置的仅仅只有与其相邻的两个珠子
但是整体来看，当我们拿起某个珠子时，因其连锁反应，这一串上面所有珠子都被拿起了。
因此，假设4个相邻位置的点A、B、C、D，
第一步：我们通过点A、C的位置确定B点的位置
第二步：我们通过点B、D的位置确定C点的位置
此时，C点的位置是一定会移动的，这C点的位置改变，又会导致第一步中使用未移动之前C点位置信息计算得来B点的位置错误。
这样来看我们永远无法正确模拟。
其实，可以使用多次计算来近似估计点的位置更新信息，在迭代多次后，每个点的移动距离几乎可以忽略不记，此时可以近似当作绳子的约束模拟完成。
加入迭代后代码如下：
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
public float ki = 10;//劲度系数
for (int i = 0; i < lterationStep; i++)
{
foreach (var k in sticks)
{
Vector3 dir = (k.particle_2.position - k.particle_1.position).normalized;//绳子方向的单位向量
float curLength = (k.particle_2.position - k.particle_1.position).magnitude;//当前节点长度
float deltaLength = curLength - stickLength;
var xDelta = ki * deltaLength * Time.fixedDeltaTime * Time.fixedDeltaTime;
if (!k.particle_1.isLock)
{
k.particle_1.position += xDelta * dir;
}
if (!k.particle_2.isLock)
{
k.particle_2.position -= xDelta * dir;
}
}
}
此时分析代码，很容易得出，迭代次数越大，每个点与点的距离越接近我们设定好的约束长度(stickLength)
由此可以借助调整迭代次数，来模拟弹性绳和非弹性绳
借助Unity自带组件lineRenderer渲染
1
2
3
4
5
6
7
8
private void Rendering()
{
for (int i = 0; i < particleCount; i++)
{
lineRenderer.SetPosition(i, particles[i].position);
}
endParticle.transform.position = particles[particles.Count - 1].position;//更新结束点位置
}
至此，大功告成。
完整代码
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
using System.Collections.Generic;
using UnityEngine;
namespace CordSimulation
{
[RequireComponent(typeof(LineRenderer))]//添加脚本时自动添加组件‘LineRenderer’
public class Cord : MonoBehaviour
{
public GameObject startParticle;//起始位置
public GameObject endParticle;//结束位置
public bool isStartLock;
public bool isEndLock;
public int lineNumCornerVertices = 1;
[SerializeField] private float gravity = 10;
[SerializeField] private int particleCount = 10;//创建的粒子数目
[SerializeField] private int lterationStep = 5;//迭代步数，越小越接近弹性绳
public float ki = 2000;//劲度系数
private float stickLength;//约束的长度
private LineRenderer lineRenderer;
private List<Particle> particles = new List<Particle>();
private List<Stick> sticks = new List<Stick>();
private void Awake()
{
lineRenderer = this.GetComponent<LineRenderer>();
lineRenderer.positionCount = particleCount;//设置linerender分段
lineRenderer.numCornerVertices = lineNumCornerVertices;//设置lineRender转角
}
void Start()
{
InitSystem();
}
//初始化系统
private void InitSystem()
{
stickLength = (startParticle.transform.position - endParticle.transform.position).magnitude / (particleCount - 1);//计算约束长度
Vector3 dir = (startParticle.transform.position - endParticle.transform.position).normalized;//开始点与结尾点的方向向量
for (int i = 0; i < particleCount; i++)
{
particles.Add(new Particle(startParticle.transform.position - dir * stickLength * i));
}
for (int i = 0; i < particles.Count - 1; i++)
{
sticks.Add(new Stick(particles[i], particles[i + 1]));
}
}
private void Simulation()
{
//实时更新绳子两端是否锁定
particles[particles.Count - 1].isLock = isEndLock;
particles[0].isLock = isStartLock;
var verletDel = new Vector3();
for (int i = 0; i < particles.Count; i++)
{
if (particles[i].isLock == false)
{
Vector3 temp = particles[i].position;
verletDel = particles[i].position - particles[i].oldPosition + new Vector3(0, -gravity, 0) * Time.fixedDeltaTime * Time.fixedDeltaTime;//Verlet 算法算法下的位置增量
particles[i].position += verletDel;
particles[i].oldPosition = temp;//更新oldPosition
}
}
//约束
for (int i = 0; i < lterationStep; i++)
{
// //约束方式1：位置偏移
// int j = 0;
// foreach (var k in sticks)
// {
// Vector3 dir = (k.particle_2.position - k.particle_1.position).normalized;//绳子方向的单位向量
// float curLength = (k.particle_2.position - k.particle_1.position).magnitude;//当前节点长度
// float deltaLength = curLength - stickLength;
// if (isStartLock)
// {
// sticks[0].particle_1 = new Particle(startParticle.transform.position);
// particles[0] = sticks[0].particle_1;
// }
// if (j == 0)
// {
// if (!k.particle_2.isLock)
// {
// k.particle_2.position -= deltaLength * dir;
// }
// }
// else if (j == sticks.Count - 1 && isEndLock)
// {
// if (!k.particle_1.isLock)
// {
// k.particle_1.position += deltaLength * dir;
// }
// }
// else if (j != 0)
// {
// if (!k.particle_1.isLock)
// {
// k.particle_1.position += 0.5f * deltaLength * dir;
// }
// if (!k.particle_2.isLock)
// {
// k.particle_2.position -= 0.5f * deltaLength * dir;
// }
// }
// j++;
// Debug.Log(j);
// }
//约束方式2
foreach (var k in sticks)
{
if (isStartLock)
{
sticks[0].particle_1 = new Particle(startParticle.transform.position);
particles[0] = sticks[0].particle_1;
}
Vector3 dir = (k.particle_2.position - k.particle_1.position).normalized;//绳子方向的单位向量
float curLength = (k.particle_2.position - k.particle_1.position).magnitude;//当前节点长度
float deltaLength = curLength - stickLength;
var xDelta = ki * deltaLength * Time.fixedDeltaTime * Time.fixedDeltaTime;
if (!k.particle_1.isLock)
{
k.particle_1.position += xDelta * dir;
}
if (!k.particle_2.isLock)
{
k.particle_2.position -= xDelta * dir;
}
}
}
}
private void FixedUpdate()
{
Simulation();
}
private void LateUpdate()
{
Rendering();
}
private void Rendering()
{
for (int i = 0; i < particleCount; i++)
{
lineRenderer.SetPosition(i, particles[i].position);
}
endParticle.transform.position = particles[particles.Count - 1].position;//更新结束点位置
}
}
//创建粒子类
public class Particle
{
public Vector3 position;//当前位置
public Vector3 oldPosition;//上一帧位置
public bool isLock = false;//是否锁定
public Particle(Vector3 a)
{
position = a;
oldPosition = a;
}
}
//创建约束
public class Stick
{
public Particle particle_1;
public Particle particle_2;
public Stick(Particle a, Particle b)
{
particle_1 = a;
particle_2 = b;
}
}
}
注意：
1
2
private int lterationStep = 5;//迭代步数，越小越接近弹性绳
public float ki = 2000;//劲度系数
这两项不宜设置过大
lterationStep过大会导致性能开销增加
ki过大会导致 $\Delta t$内点的位置增量变大，超过一定限度会使点永动，并无限拉伸撕裂。
注意：
代码中包含两种约束方式：
约束1为直接使用位置约束，其相关约束条件只有lterationStep
约束2为使用胡克定律位置约束
在约束2中调整的ki的值，能在更小lterationStep下，找到符合预期的绳索。
在约束1中因其直接使用位置约束，故不会出现约束2中ki值过大的无限拉伸撕裂情况。
参考文章
本文由作者按照 CC BY 4.0 进行授权