The Polynomial That Fixed 30 Years of Cloth Simulation

改变30年布料模拟的多项式 | Towards Data Science

物理

改变30年布料模拟的多项式

一家日本电子商务公司如何发表近年来最严谨的布料模拟论文，为什么对数障碍在紧绷的应变限制下会失败，以及如何在Python中复现其核心见解。

Ferran Alia

2026年6月8日

15分钟阅读

分享

Zozo的求解器示例 © ZOZO Technologies / Ryoichi Ando（Apache 2.0许可证）

几乎所有构建过的3D动画流水线代码中都隐藏着一个错误。当角色的袖子穿过自己的躯干时，当裙子在行走过程中穿过腿部时，当模拟的桌布像光一样穿过桌子边缘时，这个错误就会出现。高端工作室在电影上映前，会花费数千小时的艺术家时间逐帧查找这个错误。而在视频游戏和实时图形中，这个错误通常直接被发布。

剪切错误。布料不知道自己是固体。

2024年底，一家名为ZOZO（日本最大的时尚电子商务平台）的公司，在ACM Transactions on Graphics（该领域的顶级期刊）上发表了一篇名为《具有弹性包容动态刚度的立方障碍》的物理求解器，标题非常平淡。当时几乎没有引起任何新闻报道。随后，社区开始用他们拥有的所有其他求解器进行并排比较，而ZOZO的演示一直胜出。

来自算法简化版本的布料模拟，作者提供

这些演示展示了像布料一样行为的布料。不是橡胶，不是果冻，也不是幽灵：是真正的织物。五层网格覆盖在球体上，每一层都与下一层接触，没有任何一个三角形侵犯了另一个的空间。在他们的大规模压力测试中，求解器处理了超过1.84亿个同时接触对，而没有一次相互穿透。

要理解为什么多项式替换是一个如此重大的突破，你必须回顾之前几十年的数学妥协，以及如何用几百行Python代码复现ZOZO的核心逻辑。

没有人能完全解决的问题

在物理工程师关心的层面上，织物是由顶点通过三角形连接的网格。每个顶点都有一个质量和速度，每一步时间，引擎都会同时对所有顶点积分牛顿第二定律（F = ma）。力来自重力、织物的弹性阻力，以及布料接触的任何表面。弹性部分相对较容易解决：将每个三角形建模为一小块弹性材料，定义其如何抵抗拉伸和剪切，从变形梯度中推导出力。这是标准的有限元方法，效果足够好，不是瓶颈。

碰撞部分是所有聪明想法最终都会失败的地方。

不同碰撞处理的示意图，作者提供

天真的方法是惩罚力：当两个表面重叠时，用与穿透深度成比例的力将它们推开。简单且快速，但本质上是被动的：排斥力只有在穿透发生后才会激活。在大时间步长下，表面可能在力生效之前直接穿过彼此，如果你让惩罚力足够强以捕捉快速移动的表面，就会产生数值不稳定性，使模拟崩溃。

另一种经典方法是约束投影：在每一步结束时，检测交集并明确地将顶点分开以解决它们。这种方法在实践中更稳定，但它不能正确地保持能量守恒，长时间模拟会出现漂移现象，并且在有数千个接触点的场景中，要使其稳健变得真正困难，解决一个交集的同时，又会创造出两个新的交集。

几十年来，这些就是选择：接受裁剪、接受类似橡胶的拉伸，或者接受你的模拟偶尔会爆炸，通常是这三者混合在一起。

逐步潜在接触：学术界的黄金标准

2020年，一支来自美国的学术和工业研究团队发表了逐步潜在接触（Incremental Potential Contact，简称IPC），它完全重新定义了碰撞问题。其核心思想是：不要将接触建模为一种力或约束，而是将其建模为一种能量。

具体来说，向总能量中添加一个项，当两个表面之间的距离趋近于零时，该项会变得无限大。由于模拟本质上是能量最小化，求解器会自然避免表面接近的配置——不是因为你告诉它要这么做，而是因为进入这些配置会导致能量无限增加。穿透在数学上变得不可达，而不仅仅是受到惩罚。

在每个时间步中，你解决的是：

$$ x_{t+1} = \text{argmin}(E_{kinetic}(x) + E_{elastic}(x) + B(x)) $$

其中 $ B(x) $ 是障碍能量：当任何两个表面之间的距离小于阈值 $ \hat{d} $ 时，该函数会迅速增大。IPC的障碍函数是基于对数的：

$$ B_{\log}(d)= \begin{cases} -(d-\hat{d})^{2}\log\!\left(\frac{d}{\hat{d}}\right), & d<\hat{d},\\ 0, & d\ge\hat{d}. \end{cases} $$

当 $ d \to 0 $ 时，$ \log(d/\hat{d}) \to -\infty $，因此 $ B_{\log}(d) \to +\infty $。表面之间的分离可以渐近地接近零，但永远无法达到零。这是一个真正优雅的想法，结果也证明了这一点：薄布片堆叠在一起，每一层都正确地放在下一层上，没有相互穿透。

障碍函数的形状，作者提供

为什么对数形状是错误的

这就是ZOZO论文所利用的微妙之处。

对数障碍函数的二阶导数（其曲率，直接影响接触附近牛顿系统的刚度）随着表面接近而无限增长：

$$ \lim_{d\to 0^+} B_{\log}''(d)=\infty $$

这在数学上是必要的：你需要障碍函数在零点变得无限陡峭，以确保该点的能量无限大。但这对牛顿方法造成了严重的问题。每当两个表面接近时，每次迭代中需要求解的线性系统会变得病态：小的数值误差会被放大，求解器需要更多迭代才能收敛，在紧密配置中甚至可能完全失败。

当你考虑到实际织物的需求时，情况会变得更糟。棉布在拉伸超过百分之一或两之前就会产生强烈的抵抗。在模拟中，这意味着一个严格的应变限制：三角形不能从其静止面积扩展超过1-5%。要使用对数障碍函数来强制执行这一点，你需要 $ \hat{d} $ 非常小，这将障碍函数推入一个非常狭窄的间隙，迫使表面在障碍激活之前变得非常接近，这使得已经爆炸的曲率几乎没有空间来工作。你的牛顿系统的条件数会急剧上升。

除此之外，还有时间影响（TOI）锁定机制。IPC 使用加法连续碰撞检测（ACCD）来预测在给定时间步长内表面何时会发生碰撞，并通过减小步长来防止碰撞发生。在简单的场景中，这种方法是可行的。然而，在布料与布料之间的接触（数千个三角形相互滑动和摩擦）中，检测会级联发生，求解器被迫使用越来越小的时间步长，导致进展停滞。

立方体屏障：ZOZO 实际上做了哪些改变

ZOZO 的论文针对这两种失败模式进行了两个有针对性的修改。

首先：将对数函数替换为立方函数。

他们不再使用对数屏障，而是使用归一化距离的立方多项式。一个典型的表达形式如下：

$$ B_{\mathrm{cubic}}(d)= \begin{cases} \kappa\left(1-\frac{d}{\hat d}\right)^3, & d<\hat d,\\[6pt] 0, & d\ge \hat d. \end{cases} $$

其中 κ 是一个刚度参数。在整个区间 [0, d̂] 内，该函数值为正，接触时等于 κ，阈值时为零。其关键特性是它的二阶导数：

$$ B_{\mathrm{cubic}}”(d) = \frac{6\kappa}{\hat d^2} \left(1-\frac{d}{\hat d}\right) $$

在阈值 d = d̂ 处，曲率正好为零。随着 d 向接触点减少，曲率逐渐增加，直到 d = 0 时达到 6κ/d̂²。与对数屏障相比：

$$ \begin{aligned} B_{\log}”(d) &\to \infty && \text{as } d\to 0^+ && (\text{条件数爆炸})\\[4pt] B_{\mathrm{cubic}}”(d) &\to \frac{6\kappa}{\hat d^2} && \text{as } d\to 0^+ && (\text{始终有界}) \end{aligned} $$

这种逐渐增加的曲率就是关键所在。可以把它想象成一个配置良好的冷却风扇：当 CPU 变热时，风扇转得更快，而不是在电脑一启动时就全速运转。对数屏障（甚至二次函数的替代方案）在两个表面进入阈值区域的瞬间就达到最大刚度。而立方函数则温和地启动，只有当表面实际接近时才逐渐变硬。结果是，无论应变限制有多紧，牛顿系统在整个过程中都保持良好的条件。

曲率变化，作者提供

立方函数在零点时永远不会达到无限能量，这引发了一个显而易见的问题：没有这个无限尖峰，到底是什么阻止了穿透？我们需要诚实地面对这个问题，因为这是论文有时被误解的地方。

一个有限的屏障本身无法保证零穿透。如果一个顶点具有足够的动能，它可以在一个离散的时间步长内越过有限能量的障碍，直接穿过屏障，而求解器尚未有机会做出反应。ZOZO 的求解器仍然依赖 ACCD 来处理这种情况：它监视即将发生的碰撞，并在穿透发生之前减小时间步长，这与标准 IPC 的做法完全相同。立方屏障改变的并不是是否需要 ACCD，而是当表面接近时，牛顿系统退化程度的严重性。ACCD 负责处理严格的几何保证；立方函数使优化过程在数值上足够稳定，从而可靠地达到目标。这两种机制是互补的，而不是相互竞争的。

第二：将屏障刚度与局部弹性相结合。

在标准的 IPC 中，κ 是一个固定的全局常数，仅在初始化时进行调整，之后保持不变。而 ZOZO 则动态地（每个元素、每个牛顿步骤）根据材料在该点的局部弹性属性来计算 κ。这个思路很直观：一个已经接近其弹性极限的布料三角形，已经对变形产生了较强的抵抗，材料本身对整体刚度的贡献已经存在。此时，障碍项不需要独自去补偿。通过将 κ 作为局部弹性 Hessian 的函数进行计算（详细的推导见论文，值得一读），求解器确保了障碍项和弹性项在刚度上不会相互冲突。即使在紧致的应变限制下，整体系统依然保持良好的条件性，而立方多项式在接触点的有限峰值在实践中足以防止穿透。

结果是牛顿迭代次数减少，在应变极限下行为更加稳定，且能够扩展到足以使对数障碍求解器崩溃的接触数量。

在悬垂形状上的微妙差异，作者提供

在 Python 中实现核心思想

完整的求解器是用 CUDA、Rust 和 Python 包装 C++ 实现的（代码库中分别约占 49%、13% 和 22%）。但两种障碍项之间在条件性上的核心差异完全可以在一个简单的二维布料网格上通过一个笔记本实现，并且这也是最清晰地展示为什么多项式选择重要的方式。

设置是一个悬垂测试：一个平面三角化布料网格在重力作用下落到一个球体上，顶点之间可以相互接触，也可以与球体表面接触。我们在 2% 的应变限制下运行两次，一次使用对数障碍项，一次使用立方障碍项，并记录每一步牛顿系统的条件数。

一个微妙的陷阱是：解析障碍 Hessian 块中包含一个切向项，其特征值为 B′(d)/dist，对于任何排斥障碍项，这个值都是负的。不定的 Hessian 会产生非下降的牛顿方向；此时 Armijo 线搜索会拒绝每一步，模拟完全冻结。解决方法是仅保留法向分量，将 Hessian 块投影到其半正定部分。这在 IPC 文献中是标准做法，但在大多数教程实现中却缺失。

import numpy as np

def create_cloth_mesh(nx, ny, width, y_bottom):
    """三角化布料网格"""
    dx    = width / (nx - 1)
    xs    = np.linspace(-width / 2, width / 2, nx)
    ys    = [y_bottom + j * dx for j in range(ny)]
    verts = np.array([[x, y] for y in ys for x in xs], dtype=float)

    faces = []
    for j in range(ny - 1):
        for i in range(nx - 1):
            a = j * nx + i;  b = a + 1
            c = (j + 1) * nx + i;  d = c + 1
            faces += [[a, b, d], [a, d, c]]          # 每个四边形生成两个三角形
    faces = np.array(faces, dtype=int)

    edge_set = set()
    for f in faces:
        for k in range(3):
            edge_set.add(tuple(sorted([f[k], f[(k + 1) % 3]])))
    edges       = np.array(sorted(edge_set), dtype=int)
    rest_lengths = np.linalg.norm(verts[edges[:, 0]] - verts[edges[:, 1]], axis=1)

    pinned = np.arange((ny - 1) * nx, ny * nx)      # 顶部行固定不动
    return verts, edges, rest_lengths, faces, pinned

def sphere_sdf(x, cx, cy, r):
    return np.linalg.norm(x - np.array([cx, cy]), axis=1) - r

弹性能量是每个三角形的标准线性共旋转应变。屏障能量位于其上：对于每对距离小于d̂的非相邻顶点，计算B(d)，并将它的梯度和海森矩阵累加到全局系统中。两种屏障实现并排存在，这样结构差异在代码中立即可见——一种计算对数并观察其海森矩阵在接触附近发散，另一种计算多项式并保持稳定。

# ── IPC 对数屏障 ────────────────────────────────────────────────────────

def log_barrier(d, d_hat, kappa):
    if d >= d_hat: return 0.0
    d = max(d, 1e-10)           # 截断：log(0) = -∞, 能量 = +∞
    return -kappa * (d - d_hat)**2 * np.log(d / d_hat)

def log_barrier_d2(d, d_hat, kappa):
    """曲率随着d趋近于0而无限制增长。
       导致条件数趋向无穷大。"""
    if d >= d_hat: return 0.0
    d = max(d, 1e-10);  t = d - d_hat
    return -kappa * (2*np.log(d / d_hat) + 4*t / d - t**2 / d**2)

# ── ZOZO 三次屏障 ───────────────────────────────────────────────────────

def cubic_barrier(d, d_hat, kappa):
    """B(d) = κ·(1 - d/d̂)³。接触时有限；阈值时为零。"""
    if d >= d_hat: return 0.0
    return kappa * (1.0 - d / d_hat)**3

def cubic_barrier_d2(d, d_hat, kappa):
    """曲率——无论d如何接近零，始终有界于6κ/d̂²。"""
    if d >= d_hat or d <= 0.0: return 0.0
    return 6.0 * kappa / d_hat**2 * (1.0 - d / d_hat)

时间积分器是隐式欧拉法，故意设计得较慢，因为目标是透明性。每次矩阵装配、每次条件数测量、每次线搜索步骤都暴露出来，而不是隐藏在库调用中。当对数屏障的牛顿系统在接触附近开始失稳时，你可以在实时记录的条件数中看到这一过程。

def run_newton_step(x, x_tilde, mass, dt, edges, rest_lengths,
                    k_stretch, sphere_params, d_hat, kappa, barrier_type):

    E0, g, H = assemble(x, x_tilde, mass, dt, edges, rest_lengths,
                         k_stretch, sphere_params, d_hat, kappa, barrier_type)

    cond = float(np.linalg.cond(H))

    # Tikhonov 正则化
    eps     = max(1e-6 * np.trace(H) / H.shape[0], 1e-8)
    H_reg   = H + eps * np.eye(H.shape[0])
    delta   = np.linalg.solve(H_reg, -g)

    # Armijo 回溯
    alpha = 1.0
    for _ in range(max_ls):
        x_try = x + alpha * delta.reshape(-1, 2)
        E_try = total_energy(x_try, x_tilde, ...)
        if E_try <= E0 + ls_c * alpha * float(g.dot(delta)):
            break
        alpha *= ls_alpha

    return x + alpha * delta.reshape(-1, 2), cond

运行两种变体并在时间上绘制条件数。当2%的应变限制生效时，对数屏障的条件数在布料接触时迅速上升几个数量级，而三次屏障的条件数几乎不变。这是论文的核心条件数洞察，无需离开笔记本即可直观看到。

def run_simulation(barrier_type, cfg):
    x   = initial_mesh_positions.copy()
    v   = np.zeros_like(x)
    cond_per_frame = []

    for frame in range(cfg["n_frames"]):
        # 隐式欧拉预测器
        v[free_vertices] += cfg["dt"] * np.array([0, cfg["gravity_y"]])
        x_tilde = x + cfg["dt"] * v

牛顿迭代

for _ in range(cfg["max_newton"]): x, cond = run_newton_step(x, x_tilde, ..., barrier_type) cond_per_frame.append(cond) if gradient_norm(x, x_tilde, ...) < cfg["newton_tol"]: break

恢复速度

v = (x - x_tilde) / cfg["dt"] + v

return x, cond_per_frame

运行并绘图

log_drape, log_conds = run_simulation("log", CFG) cubic_drape, cubic_conds = run_simulation("cubic", CFG)

fig, ax = plt.subplots(figsize=(10, 4)) ax.semilogy(log_conds, color="#e63946", lw=2, label="Log barrier (IPC)") ax.semilogy(cubic_conds, color="#2a9d8f", lw=2, label="Cubic barrier (ZOZO)") ax.set_xlabel("Simulation frame"); ax.set_ylabel("Condition number (log scale)") ax.legend(); ax.grid(alpha=0.25)

条件数，数值越低越好，注意前几帧通常是最有问题的。作者

## 为什么一家时尚公司解决了这个问题

ZOZO 的动机指出了一个更广泛的转变，即如今严肃的图形研究正在哪里进行。ZOZO 并不是一个图形实验室，他们在线销售服装。他们对无穿透布料模拟的兴趣是商业性的：他们希望顾客在购买之前看到服装如何在身体上垂坠。在虚拟试穿流程中，一个发生穿透的模拟是不可信赖的，这直接关系到商业问题。

值得注意的是，求解器不仅仅用于布料。它能够处理壳体、固体和杆状物，即在同一个框架下处理布料、体积物体和类似纱线的结构。时尚的动机推动了这项工作，但最终的求解器是通用的。

这种商业压力塑造了每一个工程决策。求解器部署在 Docker（镜像大小约为 1 GB）中，运行在云 GPU 上，并且 README 中包含一个成本表，详细列出每个演示场景在 AWS L4 实例上每小时的费用（大约一美元）。Python API 有文档字符串，并且可以进行代码检查。GitHub Actions 连续十次运行每一个示例，使用抖动的初始位置来验证鲁棒性，将十次中的一次失败视为整个测试套件的失败。最能说明问题的是“我们是如何构建这个”的部分：早期的代码库是使用 GitHub Copilot 编写的，之后几乎所有代码都是通过与 Claude Code 和 Codex 的 vibe 编码完成的，所有代码在公开之前都经过人工审查。一个小型团队，使用现代的 AI 辅助工作流程，生产出了参考级的模拟软件，而以前这需要一个专门的学术实验室和数年的博士后研究。

在这里，“生产级”需要明确其含义，因为开发者对权衡是诚实的。这个求解器是为离线渲染和电子商务可视化构建的：在这些场景中，模拟可以运行数小时，正确性比帧率更重要，而穿透是真正不可接受的，而不是被认可的妥协。这与在 PlayStation 上发布是不同的生产定义，将两者混淆会低估求解器在其目标领域所实现的成就。（README 中包含一个明确的警告：“专为离线使用而构建；不适用于实时。”）

社区开发的 Blender 插件（如 AndoSim 和 ArzteZ-PPF-solver）让那些从未阅读过 ACM 论文的艺术家也能使用求解器，尽管官方的 ZOZO 插件尚未发布。短期内这些工具是否会取代现有工具，是另一个问题。作为这些理念的参考实现，并为后续基于这些理念的求解器奠定基础，它已经尽可能地让严肃的模拟研究变得易于接近。

立方体屏障是一个狭窄而精确的概念：将对数函数替换为一个曲率从零开始并逐渐上升的多项式函数，将刚度与局部弹性耦合，使两者按比例工作而非相互冲突，从而使布料停止穿透。数值方法中的难题往往并不需要更多的计算，而是需要一个条件更好的公式，以及一个足够具体的原因去寻找它。

作者：Ferran Alia

查看 Ferran Alia 的所有文章

计算机图形学

数据科学

深入解析

Python

分享本文

- 在 Facebook 上分享

- 在 LinkedIn 上分享

- 在 X 上分享

Towards Data Science 是一个社区出版物。提交你的见解，以触达全球读者，并通过 TDS 作者支付计划获得报酬。

请将 href 更新为你的实际投稿链接

为 TDS 撰写文章

✦ 结束 CTA ✦