Theme NexT works best with JavaScript enabled

LemonSky's Blog
0%

物理模拟学习01-质点弹簧系统01

发表于 更新于 分类于

简介

质点弹簧系统（Mass-Spring Systems）是一个非常简单的模型，很适合入门物理模拟，它可以用于模拟形变物体，如：布料、绳子等等。质点弹簧系统假设物体由一系列质点和无质量的弹簧组成，弹簧将质点连接成一个整体。虽然模拟起来十分简单，但其行为依赖每个弹簧的参数，难以调整出所需的效果。

方法

物理定律

使用弹簧质点系统时，我们利用胡克定律计算质点所受的弹力，并通过牛顿第二定律计算加速度，进而得到速度及位置。

首先使用胡克定律计算质点 $j$ 对质点 $i$ 的弹力 $\mathbf{f}_{ij}$：

其中 $\mathbf{x}_i$ 表示质点 $i$ 的位置，$l_{ij}$ 表示弹簧的静止长度，$k$ 表示弹簧刚度，$\Vert\cdot\Vert$ 表示欧式距离，$\widehat{\mathbf{x}_i - \mathbf{x}_j}$ 表示质点 $j$ 到 $i$ 的单位方向向量。

然后考虑多个弹簧的作用力，并使用牛顿第二定律计算质点 $i$ 的加速度：

进而得到关于位置的微分方程：

时间积分

求解关于速度和位置的两个微分方程，即可得到 $t$ 时刻的速度和位置：

对时间进行离散化，并采用显式时间积分（半隐式欧拉法）计算 $t$ 时刻速度和位置：

显式时间积分的当前状态只依赖过去的状态，速度快且易实现，但是却容易导致积分结果发散。时间步长 $\Delta t$ 要满足 $\Delta t \leq c\sqrt{m/k} \; (c \approx 1)$ 才可保证积分结果收敛。

Demo

基于Taichi开发Demo，其特性如下：

  • 质点和弹簧构成一个二维网格
  • 固定左上方和右上方的两个质点
  • 鼠标点击窗口施加外力

Demo的完整代码如下：

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
import taichi as ti

ti.init(arch=ti.gpu)

# 时间步长
dt = 0.001
# 网格分辨率
resolution = 10

# 质点质量
mass = 0.2
# 弹簧刚度
stiffness = 1000.0
# 速度衰减
damping = 10

# 质点速度和位置
v = ti.Vector(2, dt=ti.f32, shape=(resolution * resolution))
x = ti.Vector(2, dt=ti.f32, shape=(resolution * resolution))
# 弹簧静止长度
l = ti.var(dt=ti.f32, shape=(resolution * resolution, resolution * resolution))

# 重力
gravity = ti.Vector(2, dt=ti.f32, shape=())
# 外力
ef = ti.Vector(2, dt=ti.f32, shape=())

@ti.kernel
def simulate():
    # 计算受力和速度
    for i in range(resolution * resolution):
        v[i] *= ti.exp(-dt * damping)
        f = gravity * mass + ef
        for j in range(resolution * resolution):
            if l[i, j] != 0:
                # 胡克定律
                dis = ti.sqrt((x[i] - x[j]).dot(x[i] - x[j]))
                f += -stiffness * (dis - l[i, j]) * (x[i] - x[j]) / dis
        v[i] += f / mass * dt
    # 固定住上边两个质点
    v[0] = [0, 0]
    v[resolution - 1] = [0, 0]
    # 更新位置
    for i in range(resolution * resolution):
        x[i] += v[i] * dt

# 初始化
s = 0.6 / (resolution - 1)
for i in range(resolution):
    for j in range(resolution - 1):
        l[i*resolution+j, i*resolution+j+1] = s * 0.6
        l[i*resolution+j+1, i*resolution+j] = s * 0.6
        l[j*resolution+i, (j+1)*resolution+i] = s * 0.6
        l[(j+1)*resolution+i, j*resolution+i] = s * 0.6
for i in range(resolution):
    for j in range(resolution):
        x[j*resolution+i] = [0.2+i*s, 0.8+j*-s]

gravity[None] = [0, -9.8]

gui = ti.GUI('Mass Spring System')
while True:
    # 处理事件
    for e in gui.get_events():
        if e.key == ti.GUI.ESCAPE:
            exit()
        elif e.key == ti.GUI.LMB and e.type == ti.GUI.PRESS:
            ef[None] = [(e.pos[0] - 0.5) * 20.0, (e.pos[1] - 0.5) * 20.0]
        elif e.key == ti.GUI.LMB and e.type == ti.GUI.RELEASE:
            ef[None] = [0, 0]
    # 模拟
    for step in range(10):
        simulate()
    # GUI显示
    pos = x.to_numpy()
    gui.circles(pos, radius=5)
    for i in range(resolution):
        for j in range(resolution-1):
            gui.line(pos[i*resolution+j], pos[i*resolution+j+1], radius=1)
            gui.line(pos[j*resolution+i], pos[(j+1)*resolution+i], radius=1)
    gui.show()

Demo演示效果：

将弹簧刚度调整为 $100000$，即导致积分结果发散，如下图：