GAMES201: 高级物理引擎实战指南 2020
|
115
|
0
|
笔记-GAMES

6801 字
|
39 分钟

1. 关于 GAMES 201 课程

  • 课程简介：本课程将会介绍基于物理的动画（Physically based animation）的基础和前沿知识，从拉格朗日、欧拉、混合欧拉-拉格朗日三大视角，介绍刚体、布料、烟雾、液体、弹塑性体（雪、泥沙、果冻、橡皮泥等）的模拟。通过本课程的学习，配合Taichi编程语言的使用，同学们可以独立从零开始编写最先进的高性能影视级物理求解器，并且利用自己的渲染器生成自己的特效动画。
  • 课程内容：Taichi语言基础、刚体、液体、烟雾、弹塑性体、PIC/FLIP法、Krylov-子空间求解器、预条件、无矩阵法、多重网格、弱形式与有限元、隐式积分器、辛积分器、拓扑优化、带符号距离场、自由表面追踪、物质点法、大规模物理效果渲染、现代处理器微架构、内存层级、并行编程、GPU编程、稀疏数据结构、可微编程…

物理引擎（Physics Engine）：A physics engine is computer software that provides an approximate simulation of certain physical systems, such as rigid body dynamics (including collision detection), soft body dynamics, and fluid dynamics, of use in the domains of computer graphics, video games and film.

Keywords: Discretization(离散化), Efficient Solvers(高效率求解器), Productivity(提高生产力), Performance(性能), Hardware Architecture(硬件架构), Data Structures(数据结构), Parallelization(并行化), Differentiability(可微性).

1.1. Taichi 编程语言 简洁介绍

Taichi 是一个嵌入在python的高性能领域特定的语言。

  • 设计初衷自娱提高计算机图形学应用的生产力和可移植性
  • 数据导向、并行、megakernels（将多个小内核合并成一个更大的内核，以减少内核启动的开销并优化硬件资源利用率）
  • 将数据结构从计算中解耦（decouple）
  • 以密集方式访问空间稀疏张量
  • 可微分编程支持

1.2. 环境搭建

在 python 环境下按如下指令搭建即可：

pip install taichi

加载 demo gallery

ti gallery

1.3. 初始化

在使用 taichi 前，一定要先初始化 taichi

ti.init(arch=ti.gpu)

最频繁使用的参数：arch

  • ti.x64/arm/cuda/vulkan/opengl/metal：指定使用某个后端；
  • ti.cpu：（默认），自动检测 x64/arm CPUs；
  • ti.gpu：按顺序尝试使用 ti.cuda/vulkan/opengl/metal，如果没有可用的 GPU 架构，会使用 CPU。

1.4. 数据类型

Taichi 是强类型（Strongly Typed，显式类型声明、严格类型检查、限制隐式转换），支持如下基本类型：

  • Signed integers：ti.i8/i16/i32/i64
  • Unsigned integers：ti.u8/u16/u32/u64
  • Float-point numbers：ti.f32/f64

最常用：ti.i32 and ti.f32，布尔值用 ti.i32 表示。

在 Taichi 中，tensors（张量）是 first-class citizens（某种数据类型或对象能够像其他数据类型一样被平等地使用、传递、赋值、返回和操作）。

  • 张量本质上是多维数组
  • 张量的元素可以是标量（var），向量（ti.Vector），或矩阵（ti.Matrix）
  • 张量的元素总是通过 a[i,j,k] 语法来访问（没有指针！）。
  • 访问越界会导致未定义行为
  • （高级）张量可以是空间稀疏的
  • 尽管数学上矩阵被视为 2D 张量，但在 Taichi 中，张量和矩阵是两个完全不同的概念，矩阵可以作为张量元素使用。

tensor 使用示例：

import taichi as ti
ti.init()
a = ti.var(dt=ti.f32, shape=(42, 63))  # A tensor of 42x63 scalars
b = ti.Vector(3, dt=ti.f32, shape=4)  # A tensor of 4x3D vectors
c = ti.Matrix(2, 2, dt=ti.f32, shape=(3, 5))  # A tensor of 3x5 2x2 matrices
loss = ti.var(dt=ti.f32, shape=())  # A (0-D) tensor of a single scalar

a[3, 4] = 1
print('a[3, 4] =', a[3, 4])  # "a[3, 4] = 1.000000"

b[2] = [6, 7, 8]
print('b[0] =', b[0][0], b[0][1], b[0][2])  

loss[None] = 3
print(loss[None])  # 3

1.5. 计算

1.5.1. kernel 和 function

在 Taichi 中，计算位于（reside）内核（kernel），kernel 无法调用 kernel。

  • Taichi 内核和函数使用的语言类似于 Python
  • Taichi 内核语言是 compiled（编译型），statically-typed（静态类型的），lexically-scoped（词法作用域的），parallel（并行）和 differentiable（可微分）
  • 内核必须使用 @ti.kernel 修饰，这样内核才会被标记并编译
  • 内核的参数和返回值必修进行 type-hinted 类型提示（比如 ti.i32）

示例：

@ti.kernel
def hello(i: ti.i32):
    a = 40
    print('Hello world!', a + i)

hello(2)  # "Hello world! 42"

@ti.kernel
def calc() -> ti.i32:
    s = 0
    for i in range(10):
        s += i
    return s  # 45

Taichi 函数可以由 Taichi 内核和其他 Taichi 函数调用。它们必须使用 @ti.func 装饰器。（Taichi 中内核和函数的区别：Python 可以直接调用 kernel 但不能直接调用 function）

示例：

@ti.func
def triple(x):
    return x * 3

@ti.kernel
def triple_array():
    for i in range(128):
        a[i] = triple(a[i])
  • Taichi 函数会被强制内联，当前不支持递归（2020）
  • Taichi 函数最多只能包含一个返回语句

1.5.2. 数学计算

大部分 Python 的数学运算函数在 Taichi 中均支持，且支持链式比较（chaining comparisons, e.g. a < b <= c != d）。

  • ti.Matrix 仅适用于小矩阵（例如 3×3）。如果需要处理 64×64 的矩阵，应该考虑使用一个 2D 张量来存储标量。
  • ti.Vector 和 ti.Matrix 相同，只不过它只有一列。
  • 区分元素级（element-wise）乘积 * （会把矩阵 A 和 B 对应位置的元素取出来相乘得到新矩阵）和矩阵乘积 @。
  • 常见的矩阵操作：
A.transpose() # 转置矩阵
A.inverse() # 求矩阵的逆
A.trace() # 求矩阵的迹（对角线元素之和）
A.determinant(type) # 求矩阵的行列式
A.normalized() # 将矩阵归一化
A.cast(type) # 转换矩阵的数据类型
R, S = ti.polar_decompose(A, ti.f32) # 极分解
U, sigma, V = ti.svd(A, ti.f32) # 奇异值分解（SVD），sigma 是对角矩阵
ti.sin(A)/cos(A) # （元素级）正弦除以余弦

1.5.3. 循环

Taichi 中的 for 循环有两种形式：

  • 范围循环（Range-for loops）：
    语法上与 Python 的 for 循环没有区别，但在 Taichi 内核中最外层作用域 使用时，会被自动并行化。范围循环是可以嵌套使用的。
  • 结构循环（Struct-for loops）：
    用于遍历张量（尤其是稀疏 sparse 张量）中已经分配的元素。（稍后会详细介绍这一点。）
  • 在 Taichi 内核中，处于最外层作用域的 for 循环会自动并行化执行。

Range-for loops 示例：

@ti.kernel
def fill():
    for i in range(10):  # Parallelized
        x[i] += i
        s = 0
        for j in range(5):  # Serialized in each parallel thread
            s += j
        y[i] = s

@ti.kernel
def fill_3d():
    # Parallelized for all 3 <= i < 8, 1 <= j < 6, 0 <= k < 9
    for i, j, k in ti.ndrange((3, 8), (1, 6), 9):
        x[i, j, k] = i + j + k

只有处于最外层作用域的 for 循环才会被并行化，而不是代码中最外层的循环。

@ti.kernel
def foo():
    for i in range(10):  # ✅ 并行化执行！
        ...

@ti.kernel
def bar(k: ti.i32):
    if k > 42:
        for i in range(10):  # ❌ 串行执行（因为不在最外层作用域）
            ...

Struct-for 示例：

from taichi import ti

ti.init(arch=ti.gpu)

n = 320
pixels = ti.var(dt=ti.f32, shape=(n * 2, n))

@ti.kernel
def paint(t: ti.f32):
    for i, j in pixels:  # ✅ 并行遍历所有像素
        pixels[i, j] = i * 0.001 + j * 0.002 + t

paint(0.3)
# Struct-for 循环会遍历张量 pixels 中所有的坐标，也就是：(0, 0), (0, 1), (0, 2), ..., (0, 319),(1, 0), ..., (639, 319)

1.5.4. 原子操作

在 Taichi 中，增强赋值操作（例如 x[i] += 1）会自动执行为 原子操作，以确保并行环境下数据安全。

当你在并行环境中修改全局变量时，务必使用原子操作，否则结果将是错误的。例如，想要对张量 x 中所有元素求和：

@ti.kernel
def sum():
    for i in x:
        # 方法 1：正确，使用 +=，自动原子操作
        total[None] += x[i]

        # 方法 2：正确，显式使用 atomic_add
        ti.atomic_add(total[None], x[i])

        # 方法 3：❌ 错误，非原子操作，可能导致结果错误
        total[None] = total[None] + x[i]

1.5.5. Taichi-scope v.s. Python-scope

  • Taichi 作用域（Taichi-scope）：
    所有被 @ti.kernel 或 @ti.func 装饰的代码。
  • Python 作用域（Python-scope）：
    所有在 Taichi 作用域之外的普通 Python 代码。
  • Note：
    • Taichi 作用域 中的代码会被 Taichi 编译器编译，并在 并行设备上执行（如 GPU 或 CPU 多线程）。
    • Python 作用域 中的代码是普通的 Python 代码，由 Python 解释器 逐行执行。

1.5.6. Taichi 程序的执行阶段（Phases）

  1. 初始化阶段（Initialization）：ti.init(…)
  2. 张量分配阶段（Tensor Allocation）：ti.var, ti.Vector, ti.Matrix
  3. 计算阶段（Computation）（启动内核，在 Python 作用域中访问张量）：在初次加载一个 kernel 或在 Python-scope 访问一个 tensor 的时候会从张量分配阶段转向计算阶段，此时无法再分配新的张量（2020）
  4. （可选）重置阶段（Optional Reset）：重置 Taichi 状态（清除所有已分配的张量和内核，释放内存）：ti.reset()

1.5.7. 综合示例：fractal.py

import taichi as ti

ti.init(arch=ti.gpu)

n = 320
pixels = ti.var(dt=ti.f32, shape=(n * 2, n))

@ti.func
def complex_sqr(z):
    return ti.Vector([z[0] ** 2 - z[1] ** 2, z[1] * z[0] * 2])

@ti.kernel
def paint(t: ti.f32):
    for i, j in pixels:  # 并行绘制所有像素
        c = ti.Vector([-0.8, ti.cos(t) * 0.2])
        z = ti.Vector([i / n - 1, j / n - 0.5]) * 2
        iterations = 0
        while z.norm() < 20 and iterations < 50:
            z = complex_sqr(z) + c
            iterations += 1
        pixels[i, j] = 1 - iterations * 0.02

gui = ti.GUI("Julia Set", res=(n * 2, n))

for i in range(1000000):
    paint(i * 0.03)
    gui.set_image(pixels)
    gui.show()

1.6. Debug mode

以 debug mod 初始化 Taichi：ti.init(debug=True, arch=ti.cpu) ，启用越界检查（仅限 CPU 后端）。

示例：

import taichi as ti

ti.init(debug=True, arch=ti.cpu)

a = ti.var(ti.i32, shape=(10))
b = ti.var(ti.i32, shape=(10))

@ti.kernel
def shift():
    for i in range(10):
        a[i] = b[i + 1]  # 在 Runtime error in debug mode

shift()

Python 代码格式化可以使用 yapf

2. 拉格朗日视角（1）：Mass-Spring System and SPH

2.1. Lagrangian v.s. Eulerian: Two Views of Continuums（连续体）

  1. Lagrangian（拉格朗日视角）
    • 核心思想：跟踪每一个粒子（或质点）随时间的运动和状态变化。
    • 像是“坐在物体上观察世界”。
    • 常用于：固体力学、粒子方法（如 MPM、SPH）。
    • 例子：追踪一滴水的路径，从出发点到当前位置。
  2. Eulerian（欧拉视角）
    • 核心思想：固定在空间的某个位置，观察有多少粒子经过这里，它们的速度、密度、压力等如何变化。
    • 像是“站在原地看物体经过你”。
    • 常用于：传统流体模拟（如基于网格的 Navier-Stokes 解法）。
    • 例子：在河边观察水流经过你的位置。

2.2. Mass-Spring Systems（弹簧质点系统）

2.2.1. Force Formulation（力学建模）

Hooke’s Law: \[\mathbf{f}_{ij} = -k(\|\mathbf{x}_i – \mathbf{x}_j\|_2 – l_{ij})(\widehat{\mathbf{x}_i – \mathbf{x}_j})\]

\[\mathbf{f}_i = \sum\limits_j^{j \ne i}\mathbf{f}_{ij}\]

Newton’s second law of motion: \[\frac{\partial \mathbf{v}_i}{\partial t} = \frac{1}{m_i}\mathbf{f}_i\]

\[\frac{\partial \mathbf{x}_i}{\partial t} = \mathbf{v}_i\]

  • \(k\)：弹簧刚度（spring stiffness）
  • \(l_{ij}\)：粒子 \(i\) 和粒子 \(j\) 之间弹簧的原长（spring rest length）
  • \(m_i\)​：粒子 \(i\) 的质量（mass of particle \(i\)）
  • \(\widehat{\mathbf{x}_i – \mathbf{x}_j}\)：从粒子 \(j\) 指向粒子 \(i\) 的方向向量
  • \(\widehat{□}\)：表示单位化（归一化，normalization），即方向向量的归一版本（单位向量）

2.2.2. Time integration（时间积分）

  1. Forward Euler（explicit，显式欧拉）\[\mathbf{v}_{t+1} = \mathbf{v}_{t} + \Delta t \frac{\mathbf{f}_{t}}{m}\]\[\mathbf{x}_{t+1} = \mathbf{x}_{t} + \Delta t \mathbf{v}_{t}\]
  2. Semi-implicit Euler（aka. symplectic Euler, explicit，半隐式欧拉）\[\mathbf{v}_{t+1} = \mathbf{v}_{t} + \Delta t \frac{\mathbf{f}_{t}}{m}\]\[\mathbf{x}_{t+1} = \mathbf{x}_{t} + \Delta t \mathbf{v}_{t+1}\]
  3. Backward Euler（often with Newton’s method, implicit，隐式欧拉，需迭代求解）

简单实现示例（基于 Semi-implicit Euler）：计算新速度，与地面碰撞，计算新位置

@ti.kernel
def substep():
    n = num_particles[None]

    # 1. 计算弹簧力 + 阻尼力 + 重力，加速度 -> 更新速度
    for i in range(n):
        # 阻尼（指数衰减）
        v[i] *= ti.exp(-dt * damping[None])
        # 初始总力 = 重力
        total_force = ti.Vector(gravity) * particle_mass
        for j in range(n):
            if rest_length[i, j] != 0:
                x_ij = x[i] - x[j]
                total_force += -spring_stiffness[None] * (x_ij.norm() - rest_length[i, j]) * x_ij.normalized()
        # 更新速度
        v[i] += dt * total_force / particle_mass

    # 2. 与地面碰撞（y 轴下界）
    for i in range(n):
        if x[i].y < bottom_y:
            x[i].y = bottom_y
            v[i].y = 0

    # 3. 更新位置
    for i in range(n):
        x[i] += v[i] * dt

2.2.3. 显式 v.s. 隐式 时间积分器（time integrators）

  • 显式方法（Explicit）如：Forward Euler、Symplectic Euler、Runge-Kutta（RK） 等
    • 未来状态只依赖于过去（无需迭代）
    • 实现简单
    • 容易“爆炸”（数值不稳定）： \[\Delta t \leq c \sqrt{\frac{m}{k}} \quad (c \approx 1)\]即时间步长受限于质量和刚度（适用于弹簧质点等系统）
    • 不适合处理高刚度（stiff）材料
  • 隐式方法（Implicit）如：Backward Euler、中点法（Midpoint method）等
    • 未来状态依赖于过去与未来自身
    • 存在“先有鸡还是先有蛋”的问题 —— 需要求解一个线性或非线性方程系统
    • 实现相对复杂
    • 每一步计算开销更大，但允许使用更大的时间步长（提升整体效率）：有时能带来很大收益（如数值稳定性），但并不总是划算
    • 存在数值阻尼（Numerical Damping）和锁死现象（Locking）
  • 隐式时间积分：\[\mathbf{x}_{t+1} = \mathbf{x}_{t} + \Delta t \mathbf{v}_{t+1}\]\[\mathbf{v}_{t+1} = \mathbf{v}_{t} + \Delta t M^{-1} \mathbf{f}(\mathbf{x}_{t+1})\]
  • 消除 \(\mathbf{x}_{t+1}\)：\[\mathbf{v}_{t+1} = \mathbf{v}_{t} + \Delta t \mathbf{M}^{-1} \mathbf{f}(\mathbf{x}_{t} + \Delta t \mathbf{v}_{t+1})\]
  • 一阶牛顿线性化（First-order Newton Approximation），使用泰勒展开，对力项线性近似（仅用牛顿法的一步）：\[\mathbf{v}_{t+1} = \mathbf{v}_{t} + \Delta t \mathbf{M}^{-1} \big[\mathbf{f}(\mathbf{x}_t) + \frac{\partial \mathbf{f}}{\partial \mathbf{x}} (\mathbf{x}_t) \Delta t \mathbf{v}_{t+1}\big]\]
  • 整理（Clean up）：\[\big[ \mathbf{I} – \Delta t^2 \mathbf{M}^{-1} \frac{\partial \mathbf{f}}{\partial \mathbf{x}} (\mathbf{x}_t) \big] \mathbf{v}_{t+1} = \mathbf{v}_t + \Delta t \mathbf{M}^{-1} \mathbf{f}(\mathbf{x}_t)\]\[\mathbf{A} = \big[ \mathbf{I} – \Delta t^2 \mathbf{M}^{-1} \frac{\partial \mathbf{f}}{\partial \mathbf{x}} (\mathbf{x}_t) \big]\]\[\mathbf{b} = \mathbf{v}_t + \Delta t \mathbf{M}^{-1} \mathbf{f}(\mathbf{x}_t)\]\[\mathbf{A}\mathbf{v}_{t+1} = \mathbf{b}\]得到一个线性系统
  • 求解方式：
    • Jacobi（实现简单）：每一轮使用上一轮的解去更新所有变量
    • Gauss-Seidel 迭代法（实现简单）：使用最新计算出的值即使更新后续变量
    • Conjugate gradients 共轭梯度法

使用 Jacobi 迭代的示例（Demo）：

A = ti.var(dt=ti.f32, shape=(n, n))  # 系数矩阵 A
x = ti.var(dt=ti.f32, shape=n)       # 当前解向量 x
new_x = ti.var(dt=ti.f32, shape=n)   # 下一轮迭代解
b = ti.var(dt=ti.f32, shape=n)       # 右侧常数向量 b

@ti.kernel
def iterate():
    # 每个维度更新一次
    for i in range(n):
        r = b[i]
        for j in range(n):
            if i != j:
                r -= A[i, j] * x[j]  # 减去非对角项
        new_x[i] = r / A[i, i]       # 使用对角元素除以得到新值
    # 把 new_x 拷贝回 x，准备下一轮
    for i in range(n):
        x[i] = new_x[i]

2.2.4. 统一显式与隐式积分器

\[\big[ \mathbf{I} – \beta \Delta t^2 \mathbf{M}^{-1} \frac{\partial \mathbf{f}}{\partial \mathbf{x}} (\mathbf{x}_t) \big] \mathbf{v}_{t+1} = \mathbf{v}_t + \Delta t \mathbf{M}^{-1} \mathbf{f}(\mathbf{x}_t)\]

  1. \(\beta = 0\)：forward / semi-implicit Euler（explicit）
  2. \(\beta = 1 / 2\)：middle-point（implicit）
  3. \(\beta = 1\)：backward Euler（implicit）

2.2.5. 高效求解思路

如果有巨量的质点和弹簧（如上百万）：

2.3. 拉格朗日法流体模拟：平滑粒子流体动力学（Smoothed Particle Hydrodynamics, SPH）

高层思想：使用粒子（particles）来携带物理量的采样值，并借助一个核函数 \(W\) 来近似表示连续场（其中 \(A\) 可以是几乎任意随空间变化的物理属性：例如密度（density）、压力（pressure）等，但如果是微分项，则需要用不同的处理方式。）。

\[A(\mathbf{x}) = \sum\limits_{i}A_i \frac{m_i}{\mathbf{\rho}_i} W(\|\mathbf{x} – \mathbf{x}_j\|_2, h)\]

Figure: SPH particles and their kernel (Source: Wikipedia)

2.3.1. 使用状态方程（Equation of States, EOS）实现 SPH

也称为弱可压缩 SPH（Weakly Compressible SPH, WCSPH）。M. Becker and M. Teschner (2007). “Weakly compressible SPH for free surface flows”. In: Proceedings of the 2007 ACM SIGGRAPH/Eurographics symposium on Computer animation. Eurographics Association, pp. 209–217.

动量（Momentum）方程：（\(D\)：物质导数（粒子关于时间的导数） \(\rho\)：密度 \(B\)：体积模量（bulk modulus）\(\gamma\)：常数，通常 \(\sim 7\)）

\[\frac{D\mathbf{v}}{Dt} = -\frac{1}{\rho}\nabla p + \mathbf{g}, \quad p = B\big[\big( \frac{\rho}{\rho_0} \big)^{\gamma} – 1\big]\]\[A(\mathbf{x}) = \sum\limits_{i}A_i \frac{m_i}{\mathbf{\rho}_i} W(\|\mathbf{x} – \mathbf{x}_j\|_2, h), \quad \rho_i = \sum\limits_{j}m_jW(\|\mathbf{x}_i – \mathbf{x}_j\|_2, h)\]

Extras：表面张力（surface tension），粘性项（viscosity）；tutorial：D. Koschier et al. (2019). “Smoothed Particle Hydrodynamics Techniques for the Physics Based Simulation of Fluids and Solids”, website.

2.3.2. SPH 中的梯度

\[A(\mathbf{x}) = \sum\limits_{i}A_i \frac{m_i}{\mathbf{\rho}_i} W(\|\mathbf{x} – \mathbf{x}_j\|_2, h)\]\[\nabla A_i = \rho_i \sum\limits_{j}m_j \big(\frac{A_i}{\rho_i^2} + \frac{A_j}{\rho_j^2}\big) \nabla_{\mathbf{x}_i} W(\|\mathbf{x} – \mathbf{x}_j\|_2, h)\]

并不是很准确且无法推导得到，但是至少是对称且动量守恒，能够计算 \(\nabla p_i\)，如果不用如上公式，系统容易出现各种各样的不稳定性。

扩展：高阶梯度 Laplace operator (viscosity etc.)…

2.3.3. SPH 模拟循环

回顾动量方程和状态方程：\[\frac{D\mathbf{v}}{Dt} = -\frac{1}{\rho}\nabla p + \mathbf{g}, \quad p = B\big[\big( \frac{\rho}{\rho_0} \big)^{\gamma} – 1\big]\]

  1. 对于每一个粒子 \(i\)，计算 \(\rho_i = \sum\limits_{j} m_jW(\|\mathbf{x}_i – \mathbf{x}_j\|_2, h)\)（注意粒子本身不携带密度，而是通过质量与周围平滑估计）
  2. 对于每一个粒子 \(i\)，使用上述梯度算符计算 \(\nabla p_i\)
  3. Symplectic Euler step：\[\mathbf{v}_{t+1} = \mathbf{v}_t + \Delta t \frac{D\mathbf{v}}{Dt}, \quad \mathbf{x}_{t+1} = \mathbf{x}_t+\Delta t\mathbf{v}_{t+1}\]

2.3.4. SPH 的变种

2.3.5. Courant–Friedrichs–Lewy (CFL) 条件

CFL 条件是数值计算中用于控制时间步长（\(\Delta t\)）的稳定性约束，特别适用于显式时间积分方法（explicit integrators），用于解决 stiff 过大时的数值爆炸问题。

公式（用于估算最大稳定时间步长）：\[C = \frac{u \Delta t}{\Delta x} \le C_{\text{max}} \sim 1\]

  • \(C\)：CFL 数（Courant 数）
  • \(\Delta t\)：时间步长
  • \(\Delta x\)：空间步长（如粒子半径、网格尺寸等）
  • \(u\)：最大速度（流场或粒子）
  • \(C_{\text{max}}\)：稳定性所允许的最大 CFL 数，通常 \(\le 1\)
  • CFL 条件用于显式积分法中估计“允许的最大时间步长”，确保仿真在数值上稳定。若违反 CFL 条件，可能导致数值发散或模拟崩溃。常见的 \(C_{\text{max}}\)​ 推荐值：
    • SPH: \(\sim 0.4\)
    • MPM: \(0.3 \sim 1\)
    • FLIP fluid (smoke): \(1 \sim 5+\)

2.3.6. 加速 SPH：邻居搜索

  • 目前为止，每一步 SPH 子步骤的时间复杂度是 \(\mathcal{O}(n^2)\)，因为每个粒子都要与所有其他粒子计算交互。这种计算成本太高，在实际中不可接受。
  • 在实际应用中，通常会构建一种空间数据结构（如体素网格 voxel grid 或哈希表（Compact hashing）），用来加速邻域查找。这样可以显著降低时间复杂度：\(\mathcal{O}(n^2) \rightarrow \mathcal{O}(n)\)
Figure: Neighborhood search with hashing. (Source: Koschier et al. 2019.)

2.3.7. 其他基于粒子的模拟方法

  1. 离散元方法（Discrete Element Method，DEM）
  2. 移动粒子半隐式法（Moving Particle Semi-implicit，MPS）
  3. 功率粒子法（Power Particles）
  4. 基于周边动力学的弹簧-质点断裂模拟（A Peridynamic Perspective on Spring-Mass Fracture）

3. 拉格朗日视角（2）：形变、弹性与有限元的基础，Taichi 进阶

  • 模拟弹性材料：
    • 视觉效果
    • 实现难度不高
    • 是其他材料模拟的基础（如：黏弹性（viscoelastic）、弹塑性（elastoplastic）、黏塑性（viscoplastic））

3.1. 形变（Deformation）

  • 形变映射（Deformation map）\(\boldsymbol{\phi}\)：是一个 向量到向量的函数，用于描述从参考（静止）材料位置到变形后位置的映射关系：\[\mathbf{x}_{\text{deformed​}} = \boldsymbol{\phi} (\mathbf{x}_{\text{rest}}​)\]
  • 形变梯度（Deformation Gradient）\(\mathbf{F}\)：\[\mathbf{F} = \frac{\partial \mathbf{x}_{\text{deformed}}}{\partial \mathbf{x}_{\text{rest}}}\]
  • 形变梯度的平移不变性（Translational Invariance）：两种映射只相差一个常向量偏移，那么它们的形变梯度是相同的。即：形变梯度与整体平移无关，只与相对形变相关。
  • 体积比（deform/rest volume ratio）\(J = \text{det}(\mathbf{F})\)（Jacobi Determinante）

3.2. 弹性（Elasticity）

3.2.1. 超弹性（Hyperelasticity）

  • 超弹性材料是指其应力-应变关系可以通过一个应变能量密度函数（strain energy density function）来定义的材料，数学形式为：\[\boldsymbol{\psi} = \boldsymbol{\psi}(\mathbf{F})\]
    • 其中：\(\boldsymbol{\psi}\)：应变能量密度函数，\(\mathbf{F}\)：形变梯度
  • 直观理解：\(\boldsymbol{\psi}\) 可以被看作一个势能函数，它对材料的形变进行“惩罚”（即能量越高表示形变越大、不自然）。模拟时，系统会试图最小化总能量 \(\boldsymbol{\psi}\)，从而实现弹性恢复。
  • “应力”（Stress）：表示材料内部因变形而产生的弹性力。
  • “应变”（Strain）：在这里可以近似理解为形变梯度 \(\mathbf{F}\)。
  • 注意区分：
    • \(\boldsymbol{\psi}\)：应变能量密度函数
    • \(\boldsymbol{\phi}\)：形变映射函数，描述物体从参考状态到当前状态的位置变化

3.2.2. 应力张量（Stress Tensor）

  • 应力（stress）描述的是材料内部微小单元之间由于变形而产生的相互作用力，即一个微元对其邻近单元施加的内力密度。
  • 应力张量（stress tensor）是一个 \(3 \times 3\) 的矩阵或二阶张量，用于全面描述各方向上的应力分布。在一个给定方向的微小截面上，其单位法向量为 \(\mathbf{n}\)，通过应力张量 \(\boldsymbol{\sigma}\) 与 \(\mathbf{n}\) 点乘：\[\text{牵引力（traction）}~ \mathbf{t} = \boldsymbol{\sigma}^T \mathbf{n}\] 可得到该截面上单位面积所受的牵引力（traction vector），即邻近单元施加在该面的内力密度。
  • 根据不同的需求会使用不同类型的应力张量：
    • 第一类 Piola–Kirchhoff 应力张量（PK1）\[\mathbf{P}(\mathbf{F}) = \frac{\partial \boldsymbol{\psi} (\mathbf{F})}{\partial \mathbf{F}}​\]
      • 表示参考空间（rest space）中的应力
      • 优点：计算简单，广泛用于模拟中
      • 缺点：不对称，不守恒动量
    • Kirchhoff 应力张量：\(\boldsymbol{\tau}\)
      • 一种中间变量，应力张量乘以体积比 J=det⁡(F)J = \det(\mathbf{F})J=det(F) 得到
      • 与 Cauchy 应力相关
    • Cauchy 应力张量：\(\boldsymbol{\sigma}\)
      • 在当前变形空间下定义的真实应力
      • 对称张量，因为角动量守恒要求内力对称
  • 应力张量间关系：\[\boldsymbol{\tau} = J\boldsymbol{\sigma} = \mathbf{P} \mathbf{F}^T \quad \mathbf{P} = J \boldsymbol{\sigma} \mathbf{F}^{-T}\]
  • 直观理解（关于 \(\mathbf{P} = J \boldsymbol{\sigma} \mathbf{F}^{-T}\)）：
    • \(J\)：体积变化因子
    • \(\boldsymbol{\sigma}\)：真实空间中的应力
    • \(\mathbf{F}^{-T}\)：将法向量 n\mathbf{n}n 从当前空间变换回参考空间
    • 使用 \(\mathbf{F}^{-T}\) 而不是 \(\mathbf{F}^{-1}\)，是因为我们在变换的是单位法向量（转置是 \(^{-T}\)）而不是点的位置

3.2.3. 弹性模量（Elastic Moduli，针对各向同性材料（isotropic materials））

  • 杨氏模量（Young’s modulus）\(E\) \[E = \frac{\sigma}{\varepsilon}\]描述材料抵抗拉伸/压缩变形的能力（应力/应变）。
  • 体积模量（Bulk modulus）\(K\) \[K = −V \frac{dV}{dP}​\]描述材料抵抗拉伸/压缩变形的能力（应力/应变）。
  • 泊松比（Poisson’s ratio）\(\nu\)
    • 表示材料在一个方向受拉时，垂直方向的收缩程度
    • 取值范围： \(\nu \in [0, 0.5)\)
    • 某些“负泊松比材料”（Auxetics）具有 \(\nu \lt 0\)，会在拉伸时横向也膨胀。
  • Lamé 常数（Lamé parameters）
    • 第一 Lamé 常数：\(\lambda\)
    • 第二 Lamé 常数：\(\mu\)，又称剪切模量（Shear modulus），也常记作 \(G\)
  • 常用的转换公式：\[K = \frac{E}{3(1 − 2\nu)}​ \quad \lambda = \frac{E\nu}{(1 + \nu)(1 – 2\nu)} \quad \mu = \frac{E}{2(1 + \nu)}​\]
  • 刻画材料习惯于使用杨氏模量和泊松比，因为 Lamé 常数不太直观，一般指定 \(E\) 和 \(\nu\) 后计算得到 \(\lambda\) 和 \(\mu\)。

GAMES笔记
暂无评论

发送评论 编辑评论 


				

			

						

				

					

						

							

								
							

							
						

					

				

				

					

						

							

								
							

							
						

					

				

				

					

						

							

								
							

							
							
													

					

				

			

			
			

				

					
				

			

				

											

							
							
						

																				
					
											
						

	

		
|´・ω・)ノ
ヾ(≧∇≦*)ゝ
(☆ω☆)
（╯‵□′）╯︵┴─┴
￣﹃￣
(/ω＼)
∠( ᐛ 」∠)＿
(๑•̀ㅁ•́ฅ)
→_→
୧(๑•̀⌄•́๑)૭
٩(ˊᗜˋ*)و
(ノ°ο°)ノ
(´இ皿இ｀)
⌇●﹏●⌇
(ฅ´ω`ฅ)
(╯°A°)╯︵○○○
φ(￣∇￣o)
ヾ(´･ ･｀｡)ノ"
( ง ᵒ̌皿ᵒ̌)ง⁼³₌₃
(ó﹏ò｡)
Σ(っ °Д °;)っ
( ,,´･ω･)ﾉ"(´っω･｀｡)
╮(╯▽╰)╭
o(*////▽////*)q
＞﹏＜
( ๑´•ω•) "(ㆆᴗㆆ)
😂
😀
😅
😊
🙂
🙃
😌
😍
😘
😜
😝
😏
😒
🙄
😳
😡
😔
😫
😱
😭
💩
👻
🙌
🖕
👍
👫
👬
👭
🌚
🌝
🙈
💊
😶
🙏
🍦
🍉
😣
Source: github.com/k4yt3x/flowerhd
	

	

		
颜文字
Emoji
小恐龙
花!
	


									

			

			
			
		

	






上一篇FLAME: Learning a model of facial shape and expression from 4D scans（SIGGRAPH2017）环境配置 与 集成扩展示例
下一篇 DampsEngine：gpu加速刚体引擎项目