一种快速混合积分法的动态布料仿真

【摘要】针对布料动态模拟中快速稳定求解的瓶颈问题，提出了一种局部自适应的混合积分方法。在每一时间步长布料网格中质点利用其模拟参数，求解一局部稳定的判断准则，据此自适应判定质点相连弹簧不同弹性力部分引起的运动方程采用何种数值积分求解，有效提高了模拟效率且可以并行计算。此外，针对线性方程组的特点，用快速超松弛迭代法代替传统的共轭梯度法来求解，进一步提高了系统的性能表现。实验表明该方法具有近似线性的复杂度，便于并行计算，并有良好的稳定视觉效果。

【关键词】局部自适应；布料模拟；质点-弹簧模型

引言

长期以来，研究人员一直寻求布料仿真的高效、稳定、实时的算法[1,2,3]，目前已提出了多种模型来解决这一问题。这些研究成果广泛应用于3D电影、网络游戏虚拟现实等领域。人物角色一直是网络游戏和虚拟现实中最重要的对象之一，而服装又覆盖人体90%区域以上，因此视觉上真实稳定的动态布料仿真对网络游戏和虚拟现实的重要性不言而喻，但具有高度真实感的布料仿真往往计算成本很高，实时性欠缺。

本文采用典型的质点-弹簧模型，利用自适应混合积分方法来对微分方程求解，巧妙运用状态方程中的局部结构参数，计算得出自适应时间步长。仿真实验表明，该方法在保证系统稳定性同时加速了求解过程。

1.改进的质点-弹簧模型

目前研究人员提出了各种基于物理性能建模的图形表达。代表性模型分为两大类，即粒子系统模型和连续介质模型。本文中采用的布料模型与choi[4]的模型相类似，是一种粒子系统模型。这种模型的突出优点是对弯曲力和压缩力统一处理，因而得以避免了“后屈曲失稳”现象的发生，从而提供了很有吸引力的稳定结果。

本文将布料抽象为一个由m×n个虚拟质点组成的网格，质点之间用无质量的、自然长度为零的弹簧相连接。其连接关系有三种：“结构弹簧”，“剪切弹簧”，“弯曲弹簧”。

网格中每一个质点与其周围相邻的四个质点通过刚性的结构弹簧（Stretch spring）相连。与其对角线上的质点之间则通过刚性较小的剪切弹簧（Shear Spring）相连。最后，质点与其空间上间隔一点相邻的质点间则通过非线性的弯曲弹簧（bend spring）相连。

布料的物理模型是基于质点一弹簧模型的粒子系统，粒子间的连接是通过弹簧力。布料的类型和运动取决于粒子间的弹簧力和布料的拓扑结构。作用在粒子上的力可分为内力和外力，内力包括弹簧力，而外力则包括重力、风力等各种阻尼力。

1.1 内力分析

用fi,y表示质点i和质点j之间的弹簧作用在质点i上的力，这里弹簧遵循hooke定律，则有

（1）

这里xi,y=xi-xj，L是弹簧原长，ki,y是弹簧的弹性系数

该模型计算质点的受力非常简单直观，并且用三种弹簧模拟了三种力：拉伸力、剪切力和弯曲抗力。但该模型的缺陷在于弯曲弹簧的处理，使得布料很难生成褶皱。Choi采用了力矩平衡方程来求得弯曲抗力，成功解决了压杆稳定性问题，同时给布料带来的褶皱效果要比前者好的多，不过其计算相当复杂，这样在要求实时计算的条件下是不可取的，因此采用了更为简便的计算方式对弯曲抗力进行求解。弯曲抗力的计算公式如下：

（2）

注意，上述弯曲力计算公式表明弯曲弹簧仅仅在弹簧被压缩（即弯曲发生）的时候才起作用。这里的ki,y要比结构弹簧的刚度小的多，因此这个力对结构力的影响并不大，只用来产生布料的有效弯曲，即褶皱。

另外，考虑到结构弹簧和剪切弹簧的一致性，对其合并处理。最终采取的模型表示方法，其弹簧连接关系与Choi一致，结构力和剪切力计算方法类似方程（1）。

这个模型结合了provot和Choi两种代表性质点弹簧模型的优点，一方面它保持了provot模型的计算简单的特点，另一方面它采用了choi模型的方法，大幅提高了布料的稳定特性，可有效生成各种褶皱。

1.2 外力分析

以上给出了布料模型和内力的计算方法。但这是一个理想化的模型，它忽略了布料运动过程中能量的耗散。现实的布料运动会受到一些外力的影响，比如空气阻力、摩擦力等等。为了使布料仿真具有足够真实感的动态表现，布料物理模型必须考虑到这些外力的影响。为了让布料运动更具真实性，在布料的运动模拟中添加一个阻尼模型是有必要的。

阻尼力是用来抵抗粒子之间剪切弹簧引起的平面内运动，该力取决与相互作用粒子间的相对运动速度，我们采用的阻尼模型如下：

（3）

这里fd是连接质点i,y之间的阻尼力，Cd是阻尼系数，vi,vy是两质点的运动速度。

这个阻尼模型在布料的模拟中广为采用，在应用上也取得了良好的效果，有效地提高了算法的稳定性，增强了布料模拟的真实感。

2.数值积分

布料受到的作用力与其运动状态的关系可以由动力学方程来描述。一般来说，布料粒子系统的二阶动力学方程可以表示为：

（4）

这里，M是弹簧一质点系统的质量矩阵，向量x包括了系统中所有质点的位置。F代表了所有的非保守力，比如摩擦力、约束力等外力。表示了所有的保守力，比如重力、结构力、剪切力等，每种保守力都关联着一种势能，比如重力能、弹性能等。

基于物理的布料变形动画可以归结为沿时间轴求解布料动力学方程的初值问题。动力学方程中位置和速度都是时间的未知函数，然而根据位置、速度和力三者之间的物理关系，我们可以将上述具有位置二阶偏导的动力学方程分解为两个具有一阶偏导的运动方程，即:

（5）

这里f代表合力，为了获得布料粒子系统在每一帧的位置和速度，我们需要利用数值积分方法对上式进行求解。一般来说，只要已知当前t时刻的布料系统中所有粒子的运动状态(xt,vt)，并由方程构造显式、隐式或半隐式积分式，那么就可以计算出一个时间步长h后布料系统的运动状态(xt+h,vt+h)，在每个时间步长内，显式方法直接计算，而隐式方法需要迭代求解。由于布料物理模型的数值积分求解是影响布料仿真系统性能的一个重要因素，因此对数值积分求解器进行改进和优化始终是快速布料动画研究的一个重要内容。

2.1 混合积分

目前已有布料仿真的实践中，显式方法存在通常稳定性问题，所有提到的隐式求解方法只有一阶，因此精度较差，容易导致运动缓慢并可能丢失褶皱等细节。其他方法甚至不对系统求解，只是给出一些合理的近似。

当然，我们不是只有显式和隐式两种积分方法可以选择，我们也可以把两种方法结合到一起同时运用，这就是混合积分法。混合积分法的核心思想就是将常微分方程系统解耦为刚性部分和非刚性部分，对于刚性部分使用隐式的方法，对于非刚性部分则使用显式的方法，从而把隐式法的求解稳定性和显式法的计算简单性很好的结合起来。我们将方程（1）所示系统的力拆分为线性部分和非线性部分，考虑到布料在空间运动的刚性特征，线性力同布料的刚性特征密切相关，以隐式方法求解；非线性力部分则直接采用显式解法。

针对布料的质点—弹簧模型，我们对结构弹簧进行隐式求解，而剪切或弯曲弹簧则显式求解。不妨将系统合力f解耦为f=gI+gH，gI代表刚性部分，gH代表非刚性部分，则有：

（6）

由上述方程知系统的刚性部分运用后向欧拉法计算而非刚性部分则采用前后向欧拉法计算。

半隐式积分法中每一时间步长都需要进行系统的单步牛顿迭代计算，因而显著降低了每一时间步长的计算效率。混合积分法中不再需要显式的计算相邻质点间的系数矩阵Jacobian矩阵。更重要的是，Jacobian矩阵是巨型稀疏矩阵，导致矩阵的向量积计算相对于共轭梯度法的巨大计算成本降低了很多。

2.2 自适应混合积分法

针对目前布料仿真领域广泛使用的半隐式计算框架，在混合积分的基础上，我们对其做进一步改进，提出了一种新的数值求解方法，即自适应混合积分法。我们不再费力分析推测对系统的哪一部分应用何种积分求解，相反我们依赖当前模拟参数和一定的稳定性判断准则迅速自主灵活决定积分方法的运用。另外，在空间上参数局部变化的区域，则依赖弹簧的连接类型来确定。与其他方法不同，我们的方法没有牺牲准确性。与混合积分方法相比，该法计算简单，降低了计算成本，改善了系统的稀疏性。

一般说来我们对弯曲弹簧显式求解，但是弯曲弹簧刚度变化范围幅度通常很大。当模拟布料对象具有较小的弯曲刚度时，这种显式求解方法的结果很不错。但是当模拟的布料对象具有大刚度特征时，这种处理方法就显得不合时宜了，此时隐式求解就显得更加合理。

设定一剪切或弯曲弹簧，定义其刚度为ks、阻尼系数为kd，原始长度为L。理想状态下我们期待建立一个稳定性准则使我们能够在模拟过程中自动确定合适的时间步长如何对系统进行解耦。另外，针对自适应模拟技术，即每一空间网格间距h是变化的，或者说局部网格参数m,L,ks,kd均是可变的，我们需要一个在时间、空间上局部均可用的判断准则。提出的一个自适应的时间步长计算方法如下：

（7）

在任一时间步长，利用当前空间网格间距h和m,ks,kd取值，我们计算某一类型弹簧连接的一对粒子作用在其上的力时，我们利用上式进行判断，如果判断为真，则跳过相关Jacobian矩阵的计算，否则就如往常一样，计算相应的Jacobian矩阵。这样我们就可以主动控制系统的解耦计算而不是被动决定。

此外，该准则仅对线性化结构或剪切弹簧有效，对小刚度的弯曲弹簧，我们一般用显式方法处理，不会遇到稳定性问题。

实验结果表明，我们的自适应混合积分法是稳定的，这一点和标准的半隐式积分法几乎不分上下。比较这两种算法的计算效率，自适应混合积分法大约降低了20%左右的计算时间，且随着网格质点数增多计算成本节约更多。

3.实验结果和分析

基于上述方法，在Windows XP sp2操作系统和Intel pentium 4 2.4G处理器的软硬件环境下，在Visual C++ 6.0平台并结合Open GL图形库，对布料进行模拟。

为了测试本文局部自适应混合积分方法和其他常用数值积分方法对计算速度的提升差异，我们利用一10×10布料网格做单摆运动仿真实验，实验中分别应用显式龙格-库塔法（Runge-Kutta法）、半隐式欧拉法结合共轭梯度法（CG法）以及本文方法各自模拟100帧，记录总的计算时间，从而求出生成每帧的平均计算时间，利用这个参数来衡量不同数值方法的计算速度。另外根据仿真中总步长数还可以求解得到平均单步计算时间。实验结果如图（2）所示。由实验结果可知，尽管显式方法单步计算时间比本文方法要小，但考虑到本文方法在时间步长上要远大于它，因此本文方法的速度还是会达到超过显示方法，布料规模越大这种实时性能提升的效果越明显。

图3给出了16×12个质点，1014根弹簧构成的布料旗帜在风力作用下效果。

实验结果表明，自适应混合积分法保持了很好的数值稳定性，仿真系统始终保持稳定，无散乱失真现象。通过仿真过程中自适应混合积分方法的使用，能够有效而逼真地动态仿真织物，不仅保证了精度，同时提高了计算效率。

参考文献

[1]Ari Stern and Mathieu Desbrun.Discrete geometric mechanics for varia tional integrators (a primer for cs).ACM SIGGRAPH ’06 Course Notes on Discrete Di erential Geometry,2006.

[2]R.Bridson,S.Marino,and R.Fedkiw.Simulation of clothing with folds and wrinkles. In ACM SIGGRAPH/Eurographics Symposium Computer Animation,pages 28-36.ACM Press,2003.

[3]E.Boxerman and U. Ascher,“Decomposing cloth,”in Proc.ACM SIGGRAPH/Eurographics Symp.on Comput.Anim.,pp.153-161,2004.

[4]K.Choi and H.Ko.Stable but responsive cloth.In Proceedings of the 29th annual conference on Computer graphics and interactive techniques,pages 604-611.ACM Press,2002.

[5]D.Baraf,A.Witkin,and M.Kass.Untangling cloth.In ACM Trans.Graphics,pages 862{870.

ACM Press,2003.

[6]J.Teran,S.Salinas-Blemker,V.Ng,and R.Fedkiw.Finite volume methods for the simulation of muscle tissue.Eurographics/ACM Symp.on Comp.Animation,2003.

[7]P.Volino and N.Magnenat-Thalmann,“Implicit midpoint integration and adaptive damping for efficient cloth simulation,”Computer Animation and Virtual Worlds,vol.16,pp.163-175,2005.