CN102156785A - 一种动态布料仿真方法 - Google Patents

Abstract

一种动态布料仿真方法，在已知当前t时刻的布料模型中所有粒子的运动状态(x_t，v_t)，通过自适应混合积分方法求解获得一个时间步长h后布料模型的运动状态(x_t+h，v_t+h)；所述自适应混合积分方法通过一个自适应的时间步长计算方法：

来确定是通过显式积分求解方法还是隐式积分求解方法对一个时间步长h后布料模型的运动状态(x_t+h，v_t+h)的求解，求解后根据当前t时刻的布料模型中所有粒子的运动状态(x_t，v_t)以及求解得到的一个时间步长h后布料模型的运动状态(x_t+h，v_t+h)，布料模型的动态仿真完成。本发明的有意效果：仿真后的动态布料逼真，仿真不存在散乱失真现象。

Description

一种动态布料仿真方法

技术领域

本发明涉及三维织物动态仿真技术领域，特别是涉及到用一种动态布料仿真方法。

背景技术

由于柔性织物真实感仿真技术的广泛应用背景，因而吸引了众多计算机图形学家和纺织学家为此进行深入研究，他们分别从不同角度，采取多种研究方法进行研究，并已取得了一定的进展。从目前的研究结果来看，对织物仿真的建模有如下一些方法。

Jerry Weil采用余弦曲线及草几何变换模拟悬垂布料。将布料悬挂在一些约束点上，利用悬链线计算布料自由悬挂时的形状，根据实际需要，通过对悬链线和约束点构成的三角形不断细分产生织物皱褶，取得了较好模拟效果。但该方法只能模拟悬垂织物的外形，无法反映织物的真实属性。

B.K.Hind利用纯几何变换进行织物形态模拟。构造了基于等距面的交互服装设计系统，将服装作为一组围绕人体模型连接在一起的衣片进行生成。服装的褶皱通过参数上附加调和函数(正弦曲线等)产生。然而这种方法仅适用余合体服装的三维涉及，当涉及宽松服装时，直线与褶皱成为主要因素，此时该方法不能正确定义服装曲面，而且没有考虑服装织物的物理性质，不能真实展现服装质感和悬垂感。

针对几何模拟技术的缺陷，基于物理和能量特性的织物模拟仿真方法得到人们的广泛重视。尽管模拟计算复杂度很高，但该方法不仅能够反映柔性织物的外形特征，而且通过对柔性织物的内部结构的深入分析，可以得到较为逼真的织物空间状态。

Breen提出了织物的粒子模型(Particle Model)，该模型假设织物是一种不连续的物质，织物的经线与纬线每个交叉点被视为一个粒子，所有这些粒子集合的整体构成织物，织物的物理特性由粒子与邻近的四个粒子间的相互作用决定。采用能量最小化原理，求解织物最小能量状态，即最终平衡的状态。

Terzopoulos在薄板弹性变形方程基础上提出了对所有弹性体都适用的通用形变模型(Elastically Deformation Model)。把物体形状和运动的描述统一起来，该模型采用拉格朗日运动学方程描述变形体的机械运动。

Provot采用弹簧-质点模型来模拟柔性织物。织物首先被离散为规则的四边域网格，网格交点为质子，质子间以无质量弹簧相连。运用牛顿运动定律，给出运动方程。由于利用线性形变的弹簧去模拟非线性的织物形变，会出现失真情况，Provot采用了基于反演动力学的直接修正法，解决了该问题。

基于物理和能量特性的织物模拟仿真方法，常常要涉及到大量复杂微分方程组的求解，成为影响仿真效率的主要问题，研究人员提出了各种积分方法来解决这一任务。显式积分方法简单灵活，易于实现，但受稳定因素影响，无法实现具有刚性特征的织物动态模拟；隐式积分方法稳定性好，却忽略了非线性因素，而且计算复杂，直接影响到仿真的最终结果和实际效率。因此需要一种考虑系统受力形变的线性和非线性特征，最大限度地利用了时间域和空间域上的模拟局部参数，协调好计算效率和布料系统模型的稳定性方法。

发明内容

本发明要解决现有动态布料仿真方法存在仿真散乱失真现象、仿真效果差的问题，提供了一种仿真效率好、不存在失真现象的动态布料仿真方法。

本发明的技术方案：

一种动态布料仿真方法，其步骤如下：

1、构建布料的物理模型，布料模型是基于质点-弹簧模型的粒子系统，将布料抽象成一个由m×n个虚拟质点组成的网格，网格中每一个质点与其周围相邻的四个质点通过刚性的结构弹簧相连，与其对角线上的质点之间则通过刚性较小的剪切弹簧相连，质点与其空间上间隔一点相邻的质点间则通过非线性的弯曲弹簧相连；

布料的类型和运动取决于质点间的弹簧力和布料的拓扑结构，布料受到的作用力与其运动状态的关系可以由动力学方程来表示，如下式：

$\stackrel{\·\·}{x}=M^{-1}(-\frac{\∂E}{\∂x}+F)---\left(1\right)$

其中，M是质点-弹簧模型的粒子系统的质量矩阵，向量x包括了系统中所有质点的位置，F代表了所有的非保守力，非保守力指摩擦力、约束力等外力，表示了所有的保守力，保守力指重力、结构力、剪切力等；

将式(1)所述的具有位置二阶偏导的动力学方程分解为两个具有一阶偏导的运动方程，即：

$\left(\begin{array}{c}\stackrel{\·}{x}\\ \stackrel{\·}{v}\end{array}\right)=\left(\begin{array}{c}v\\ M^{-1}f\end{array}\right)---\left(2\right)$

其中，f代表合力；

2、已知当前t时刻的布料模型中所有粒子的运动状态(x_t，v_t)，通过自适应混合积分方法对式(2)进行求解获得一个时间步长h后布料模型的运动状态(x_t+h，y_t+h)；所述自适应混合积分方法如下：

设定一剪切或弯曲弹簧，定义其刚度为k_s、阻尼系数为k_d，原始长度为L，每一空间网格间距h是变化的，或者说局部网格参数m，L，k_s，k_d均是可变的，一个自适应的时间步长计算方法如下：

$\mathrm{\κ}=\frac{h}{m}(k_{s}h+{2k}_d)\≤0.2---\left(3\right)$

其中m是网格中质点个数；

在任一时间步长，利用当前空间网格间距h和m，k_s，k_d取值，可以计算某一类型弹簧连接的一对粒子作用在其上的力时，通过式(3)进行判断，如果判断为真，则通过显式积分求解方法对式(2)进行求解；否则通过隐式积分求解方法对式(2)进行求解；

其中，显式积分求解方法如下：标准的显式欧拉算法表达式如下：

$\left(\begin{array}{c}\stackrel{\·}{x}\\ \stackrel{\·}{v}\end{array}\right)=\left(\begin{array}{c}v\\ M^{-1}f(x,v)\end{array}\right)---\left(4\right)$

根据前后向欧拉法，可以推导得到：

$\left(\begin{array}{c}{x}^{n+1}-x^n\\ v^{n+1}-v^n\end{array}\right)=h\left(\begin{array}{c}{v}^{n+1}\\ M^{-1}{f}^{n+1}\end{array}\right)---\left(5\right)$

其中，xⁿ，vⁿ分别表示t＝tⁿ时刻的近似值，且tⁿ⁺¹＝tⁿ+h，速度v的更新使用前向欧拉法，而x的更新则使用后向欧拉法；

隐式积分求解方法如下：隐式后向欧拉积分表达式如下：

$\left[\begin{array}{c}\mathrm{\Δ}{x}_n\\ \mathrm{\Δ}{v}_n\end{array}\right]=\left[\begin{array}{c}{x}_{n+1}-x_n\\ v_{n+1}-v_n\end{array}\right]=h\left[\begin{array}{c}{v}_n+\mathrm{\Δ}{v}_n\\ M^{-1}f(x_n+\mathrm{\Δ}{x}_n,v_n+\mathrm{\Δ}{v}_n)\end{array}\right]---\left(6\right)$

将式(6)中的f(x_n+Δx_n，v_n+Δv_n)展开为一阶泰勒等式

$f(x_n+\mathrm{\Δ}{x}_n,v_n+\mathrm{\Δ}{v}_n)=f_n+\frac{\∂f}{\∂x}\mathrm{\Δ}{x}_n+\frac{\∂f}{\∂v}\mathrm{\Δ}{v}_n---\left(7\right)$

为力关于位置和速度的Jacobian矩阵，将其代入等式(6)，通过变换最终整理得到

$\mathrm{A\Δv}=(I-h{M}^{-1}\frac{\∂f}{\∂v}-h^2{M}^{-1}\frac{\∂f}{\∂x})\mathrm{\Δv}=h{M}^{-1}(f_n+h\frac{\∂f}{\∂x}{v}_n)---\left(8\right)$

其中质量矩阵M是以3×3对角矩阵为子矩阵的对角阵，而两个Jacobian矩阵和都是子矩阵为对称矩阵的稀疏矩阵；因此，总的系数矩阵

是一个稀疏矩阵，可以采用共轭梯度法迭代求解；

3、根据当前t时刻的布料模型中所有粒子的运动状态(x_t，v_t)以及求解得到的一个时间步长h后布料模型的运动状态(x_t+h，v_t+h)，布料模型的动态仿真完成。

本发明的技术构思，一般说来布料仿真时对弯曲弹簧显式积分求解，但是弯曲弹簧刚度变化范围幅度通常很大。当模拟布料对象具有较小的弯曲刚度时，这种显式求解方法的结果很不错。但是当模拟的布料对象具有大刚度特征时，这种处理方法就显得不合时宜了，此时隐式积分求解就显得更加合理。一般布料仿真中需要这两种求解方式混合运用，为了不再费力分析推测系统的哪一部分具体应用何种积分求解，建立一判断准则，依赖当前模拟参数和稳定性判断准则迅速自主灵活决定积分方法的运用。另外，在空间上参数局部变化的区域，则依赖弹簧的连接类型来确定。与混合积分方法相比，该法计算简单，降低了计算成本，改善了系统的稀疏性。采用自适应混合积分方法对布料进行动态仿真，使得仿真后的动态布料不存在散乱失真现象，有效而逼真地动态仿真了织物。

本发明的有意效果：仿真后的动态布料逼真，仿真不存在散乱失真现象。

附图说明

图1是本发明的仿真效果示意图。

具体实施方式

一种动态布料仿真方法，其步骤如下：

1、构建布料的物理模型，布料模型是基于质点-弹簧模型的粒子系统，将布料抽象成一个由m×n个虚拟质点组成的网格，网格中每一个质点与其周围相邻的四个质点通过刚性的结构弹簧相连，与其对角线上的质点之间则通过刚性较小的剪切弹簧相连，质点与其空间上间隔一点相邻的质点间则通过非线性的弯曲弹簧相连；

布料的类型和运动取决于质点间的弹簧力和布料的拓扑结构，布料受到的作用力与其运动状态的关系可以由动力学方程来表示，如下式：

$\stackrel{\·\·}{x}=M^{-1}(-\frac{\∂E}{\∂x}+F)---\left(1\right)$

其中，M是质点-弹簧模型的粒子系统的质量矩阵，向量x包括了系统中所有质点的位置，F代表了所有的非保守力，非保守力指摩擦力、约束力等外力，

表示了所有的保守力，保守力指重力、结构力、剪切力等；

将式(1)所述的具有位置二阶偏导的动力学方程分解为两个具有一阶偏导的运动方程，即：

$\left(\begin{array}{c}\stackrel{\·}{x}\\ \stackrel{\·}{v}\end{array}\right)=\left(\begin{array}{c}v\\ M^{-1}f\end{array}\right)---\left(2\right)$

其中，f代表合力；

2、已知当前t时刻的布料模型中所有粒子的运动状态(x_t，v_t)，通过自适应混合积分方法对式(2)进行求解获得一个时间步长h后布料模型的运动状态(x_t+h，v_t+h)；所述自适应混合积分方法如下：

设定一剪切或弯曲弹簧，定义其刚度为k_s、阻尼系数为k_d，原始长度为L，每一空间网格间距h是变化的，或者说局部网格参数m，L，k_s，k_d均是可变的，一个自适应的时间步长计算方法如下：

$\mathrm{\κ}=\frac{h}{m}(k_{s}h+{2k}_d)\≤0.2---\left(3\right)$

其中m是网格中质点个数；

在任一时间步长，利用当前空间网格间距h和m，k_s，k_d取值，可以计算某一类型弹簧连接的一对粒子作用在其上的力时，通过式(3)进行判断，如果判断为真，则通过显式积分求解方法对式(2)进行求解；否则通过隐式积分求解方法对式(2)进行求解；

其中，显式积分求解方法如下：标准的显式欧拉算法表达式如下：

$\left(\begin{array}{c}\stackrel{\·}{x}\\ \stackrel{\·}{v}\end{array}\right)=\left(\begin{array}{c}v\\ M^{-1}f(x,v)\end{array}\right)---\left(4\right)$

根据前后向欧拉法，可以推导得到：

$\left(\begin{array}{c}{x}^{n+1}-x^n\\ v^{n+1}-v^n\end{array}\right)=h\left(\begin{array}{c}{v}^{n+1}\\ M^{-1}{f}^{n+1}\end{array}\right)---\left(5\right)$

其中，xⁿ，vⁿ分别表示t＝tⁿ时刻的近似值，且tⁿ⁺¹＝tⁿ+h，速度v的更新使用前向欧拉法，而x的更新则使用后向欧拉法；

隐式积分求解方法如下：隐式后向欧拉积分表达式如下：

$\left[\begin{array}{c}\mathrm{\Δ}{x}_n\\ \mathrm{\Δ}{v}_n\end{array}\right]=\left[\begin{array}{c}{x}_{n+1}-x_n\\ v_{n+1}-v_n\end{array}\right]=h\left[\begin{array}{c}{v}_n+\mathrm{\Δ}{v}_n\\ M^{-1}f(x_n+\mathrm{\Δ}{x}_n,v_n+\mathrm{\Δ}{v}_n)\end{array}\right]---\left(6\right)$

将式(6)中的f(x_n+Δx_n，v_n+Δv_n)展开为一阶泰勒等式

$f(x_n+\mathrm{\Δ}{x}_n,v_n+\mathrm{\Δ}{v}_n)=f_n+\frac{\∂f}{\∂x}\mathrm{\Δ}{x}_n+\frac{\∂f}{\∂v}\mathrm{\Δ}{v}_n---\left(7\right)$

为力关于位置和速度的Jacobian矩阵，将其代入等式(6)，通过变换最终整理得到

$\mathrm{A\Δv}=(I-h{M}^{-1}\frac{\∂f}{\∂v}-h^2{M}^{-1}\frac{\∂f}{\∂x})\mathrm{\Δv}=h{M}^{-1}(f_n+h\frac{\∂f}{\∂x}{v}_n)---\left(8\right)$

其中质量矩阵M是以3×3对角矩阵为子矩阵的对角阵，而两个Jacobian矩阵和都是子矩阵为对称矩阵的稀疏矩阵；因此，总的系数矩阵

是一个稀疏矩阵，可以采用共轭梯度法迭代求解；

3、根据当前t时刻的布料模型中所有粒子的运动状态(x_t，v_t)以及求解得到的一个时间步长h后布料模型的运动状态(x_t+h，v_t+h)，布料模型的动态仿真完成。

本发明的技术构思，一般说来对弯曲弹簧显式积分求解，但是弯曲弹簧刚度变化范围幅度通常很大。当模拟布料对象具有较小的弯曲刚度时，这种显式求解方法的结果很不错。但是当模拟的布料对象具有大刚度特征时，这种处理方法就显得不合时宜了，此时隐式积分求解就显得更加合理。一般布料仿真中需要这两种求解方式混合运用，为了不再费力分析推测系统的哪一部分具体应用何种积分求解，建立一判断准则，依赖当前模拟参数和稳定性判断准则迅速自主灵活决定积分方法的运用。另外，在空间上参数局部变化的区域，则依赖弹簧的连接类型来确定。与混合积分方法相比，该法计算简单，降低了计算成本，改善了系统的稀疏性。采用自适应混合积分方法对布料进行动态仿真，使得仿真后的动态布料不存在散乱失真现象，有效而逼真地动态仿真了织物。

仿真后的效果图如图1，仿真效果表明本发明所述的方法保持了很好的数值稳定性，仿真系统始终保持稳定，无散乱失真现象。通过仿真过程中自主判定积分方法的使用，能够有效而逼真地动态仿真织物，不仅保证了精度，同时提高了计算效率。

就本发明方法与已有的积分方法在计算速度和计算精度上做了个对比，对比如下：

每种框架各自模拟100帧，记录总的计算时间，然后可以求得每帧动画的平均计算时间；此外，根据动画中总的步长数还可以求得平均单步计算时间，如下表所示。

不同积分方法下计算速度比较

选择布料动画的前100帧作为研究对象，各种计算框架的时间步长均取为0.002s，分别计算平均位置误差和平均角度(弧度)误差，如下表所示。这里显式方法迭代n次是指步长取为0.002/n。

不同积分方法的位置/角度精度比较

从对比中可以看出本发明所述方法在计算速度和计算精度即仿真真实性方面都取得了改善。

本说明书实施例所述的内容仅仅是对发明构思的实现形式的列举，本发明的保护范围的不应当被视为仅限于实施例所陈述的具体形式，本发明的保护范围也及于本领域技术人员根据本发明构思所能够想到的等同技术手段。

Claims (1)

1.一种动态布料仿真方法，其步骤如下：

1、构建布料的物理模型，布料模型是基于质点-弹簧模型的粒子系统，将布料抽象成一个由m×n个虚拟质点组成的网格，网格中每一个质点与其周围相邻的四个质点通过刚性的结构弹簧相连，与其对角线上的质点之间则通过刚性较小的剪切弹簧相连，质点与其空间上间隔一点相邻的质点间则通过非线性的弯曲弹簧相连；

布料的类型和运动取决于质点间的弹簧力和布料的拓扑结构，布料受到的作用力与其运动状态的关系可以由动力学方程来表示，如下式：

$\stackrel{\·\·}{x}=M^{-1}(-\frac{\∂E}{\∂x}+F)---\left(1\right)$

其中，M是质点-弹簧模型的粒子系统的质量矩阵，向量x包括了系统中所有质点的位置，F代表了所有的非保守力，非保守力指摩擦力、约束力等外力，表示了所有的保守力，保守力指重力、结构力、剪切力等；

将式(1)所述的具有位置二阶偏导的动力学方程分解为两个具有一阶偏导的运动方程，即：

$\left(\begin{array}{c}\stackrel{\·}{x}\\ \stackrel{\·}{v}\end{array}\right)=\left(\begin{array}{c}v\\ M^{-1}f\end{array}\right)---\left(2\right)$

其中，f代表合力；

2、已知当前t时刻的布料模型中所有粒子的运动状态(x_t，v_t)，通过自适应混合积分方法对式(2)进行求解获得一个时间步长h后布料模型的运动状态(x_t+h，v_t+h)；所述自适应混合积分方法如下：

设定一剪切或弯曲弹簧，定义其刚度为k_s、阻尼系数为k_d，原始长度为L，每一空间网格间距h是变化的，或者说局部网格参数m，L，k_s，k_d均是可变的，一个自适应的时间步长计算方法如下：

$\mathrm{\κ}=\frac{h}{m}(k_{s}h+{2k}_d)\≤0.2---\left(3\right)$

其中m是网格中质点个数；

在任一时间步长，利用当前空间网格间距h和m，k_s，k_d取值，可以计算某一类型弹簧连接的一对粒子作用在其上的力时，通过式(3)进行判断，如果判断为真，则通过显式积分求解方法对式(2)进行求解；否则通过隐式积分求解方法对式(2)进行求解；

其中，显式积分求解方法如下：标准的显式欧拉算法表达式如下：

$\left(\begin{array}{c}\stackrel{\·}{x}\\ \stackrel{\·}{v}\end{array}\right)=\left(\begin{array}{c}v\\ M^{-1}f(x,v)\end{array}\right)---\left(4\right)$

根据前后向欧拉法，可以推导得到：

$\left(\begin{array}{c}{x}^{n+1}-x^n\\ v^{n+1}-v^n\end{array}\right)=h\left(\begin{array}{c}{v}^{n+1}\\ M^{-1}{f}^{n+1}\end{array}\right)---\left(5\right)$

其中，xⁿ，vⁿ分别表示t＝tⁿ时刻的近似值，且tⁿ⁺¹＝tⁿ+h，速度v的更新使用前向欧拉法，而x的更新则使用后向欧拉法；

隐式积分求解方法如下：隐式后向欧拉积分表达式如下：

$\left[\begin{array}{c}\mathrm{\Δ}{x}_n\\ \mathrm{\Δ}{v}_n\end{array}\right]=\left[\begin{array}{c}{x}_{n+1}-x_n\\ v_{n+1}-v_n\end{array}\right]=h\left[\begin{array}{c}{v}_n+\mathrm{\Δ}{v}_n\\ M^{-1}f(x_n+\mathrm{\Δ}{x}_n,v_n+\mathrm{\Δ}{v}_n)\end{array}\right]---\left(6\right)$

将式(6)中的f(x_n+Δx_n，v_n+Δv_n)展开为一阶泰勒等式

$f(x_n+\mathrm{\Δ}{x}_n,v_n+\mathrm{\Δ}{v}_n)=f_n+\frac{\∂f}{\∂x}\mathrm{\Δ}{x}_n+\frac{\∂f}{\∂v}\mathrm{\Δ}{v}_n---\left(7\right)$

为力关于位置和速度的Jacobian矩阵，将其代入等式(6)，通过变换最终整理得到

$\mathrm{A\Δv}=(I-h{M}^{-1}\frac{\∂f}{\∂v}-h^2{M}^{-1}\frac{\∂f}{\∂x})\mathrm{\Δv}=h{M}^{-1}(f_n+h\frac{\∂f}{\∂x}{v}_n)---\left(8\right)$

其中质量矩阵M是以3×3对角矩阵为子矩阵的对角阵，而两个Jacobian矩阵和都是子矩阵为对称矩阵的稀疏矩阵；因此，总的系数矩阵

是一个稀疏矩阵，可以采用共轭梯度法迭代求解；

3、根据当前t时刻的布料模型中所有粒子的运动状态(x_t，v_t)以及求解得到的一个时间步长h后布料模型的运动状态(x_t+h，v_t+h)，布料模型的动态仿真完成。

Images