CN101944144A - 一种基于无网格的布类仿真方法 - Google Patents

Abstract

一种基于无网格的布类仿真方法，包括：(a)将给定的布类模型Ω参数化到中面Λ上，确定变形前的构型和开始计算时的参考构型；(b)在参数化的中面Λ中选择有限数目的节点，并确定每个节点的局部子域Ω_s；(c)在每个节点的定义域上计算形函数Φ_I，并计算形函数的导数；(d)确定所有节点的高斯积分点，计算数值积分并组装刚度矩阵；(e)组装节点质量矩阵和粘性矩阵；(f)根据仿真时间t和设置的时间步间隔δt确定循环次数，将下述步骤循环t/δt次：计算旋转后的全局刚度矩阵，在全局刚度矩阵、节点质量矩阵和粘性矩阵上实施约束，通过时间积分获得系统线性方程组并求解，得到节点位移，实施无网格碰撞检测和响应算法修改节点位移，重新采样后，生成每一时间步所有采样点的位移。

Description

一种基于无网格的布类仿真方法

技术领域

本发明涉及布类模拟仿真技术领域，尤其涉及一种基于无网格的布类离散方法和无网格布类的碰撞检测方法。

背景技术

基于物理的布类仿真已经广泛的应用于虚拟现实、计算机动画以及纺织工业计算机辅助设计中。布类动画模型经历了从纯几何阶段到基于力学模型的阶段。其主要应用的模型包括离散的质点弹簧模型和连续的有限元以及有限差分模型。

质点弹簧模型由于其具有模型简单、易于实现已经广泛应用到了布类模拟中。质点弹簧是一种简单的可变形模型，属于离散模型。正如其名称所示，质点弹簧模型只是简单的包括了质点，以及与其相连接的弹簧，他们一起组成相应网络。也能够取得比较逼真的动画效果，这种方法由于缺少连续的构型，具有其固有的缺陷，如材料不能够得到连续的模拟，模拟的结果很大程度上依赖于弹簧网络，很难获得指定材料的力学特性；在模拟中缺乏相应的弹簧参数。这给应用带来本质上的困难，例如不能真实的体现不同材料的服装获得的动画效果的细微差别，这些材料力学行为上的细微差别，在计算机动画中的视觉上却存在很明显的不真实的体验。

然而另一方面，在纺织工业界，对于模拟结果的真实性要求很高，因此有技术人员转向了连续模型，如有限元方法，来解决这个问题。通过连续模型，可以精确的建模材料属性，而不依赖于离散化的精度。因此，连续模型是发展的必然。

近年来，无网格方法作为传统的基于网格方法的替代模型被引入到计算机图形领域，其具有优于基于网格的方法如有限元方法所不具备的优点。与此同时，点采样模型，这种不需要存储和保持拓扑结构信息的模型也被应用的无网格方法中。目前这一技术题目成为了当前前沿的方向，受到技术人员的广泛关注。

尽管无网格法在计算机动画中的技术已经取得了很多的成果，此项技术领域仍然具有很多有待发展的方面，无网格布模型就是其中之一。

目前已经使用的连续模型主要是简化的有限元模型，但由于有限元模型中网格本身的限制，以及布类大变形的固有特征，使得有限元模型得到的动画与网格划分有很大影响，结果差强人意。而无网格法不需重新划分网格，具有仿真布类的天然优势。但是由于无网格模型本身的技术时间不长，应用方面的技术更是少之又少，目前还没有相关文献的无网格布类应用。这使得无网格布类技术变得十分有价值。

发明内容

为了克服上述缺陷，本发明提出了一种基于无网格的布类仿真方法，包括：

(a)将给定的布类模型Ω参数化到中面Λ上，确定初始构型和参考构型；

(b)在Λ中确定有限数目的节点及其局部子域Ω_s；

(c)在每个节点的定义域上计算形函数Φ_I，计算形函数的导数；

(d)对所有全局区域内的节点循环步骤(d1)和(d2)：

(d1)确定节点的子域Ω_s上的所有高斯积分点X_Q；

(d2)对所述所有高斯积分点X_q计算数值积分，组装全局刚度矩阵K；

(e)组装节点质量矩阵M和粘性矩阵D；

(f)根据仿真时间长度t和所设置的时间步间隔δt确定循环次数，循环步骤(f1)至(f8)t/δt次：

(f1)计算旋转后的全局刚度矩阵；

(f2)计算强加约束后的全局刚度矩阵、节点质量矩阵和粘性矩阵；

(f3)根据不同的时间积分方法合成线性方程组；

(f4)求解生成的线性方程组，并获得节点位移；

(f5)更新采样点的未知量和它们的导数；

(f6)碰撞检测和响应，更新节点位移；

(f7)提取新的旋转场R；

(f8)重新采样得到采样点的位移。

其中，所述将给定的布类模型Ω参数化到中面Λ上具体包括：将布类模型分成多个单独的布片，确定所述多个布片的相对位置关系，按照相对位置关系将所述多个单独的布片组合在一起，从而保证所述多个单独的布片具有位移上的一致性，并沿经线和纬线两个方向对上述单独的布片进行参数化。

其中形函数Φ_I具体为：

Φ_I＝[R^T(X)，P^T(X)]G。

其中R^T(x)＝[R₁(x)，R₂(x)，...，R_n(x)]是以n个离散点为中心的一系列径向基函数。P^T(x)＝[p₁(x)，p₂(x)，...，p_m(x)]是m阶的单项式基函数。矩阵G是由散布点上的R和P组成的矩阵，

$G={\left[\begin{array}{cc}{R}_0& P_0\\ P_0^T& 0\end{array}\right]}^{-1},$

R₀＝[R^T(x_i)]，P₀＝[P^T(x_i)]，i＝1，2，...，n

其中，所述计算旋转后的全局刚度矩阵K具体为：对拉伸和弯曲的应变进行线性化，通过旋转场使用共旋的应变形式，获得旋转后的线性应变张量，通过极分解将变形梯度张量分解成一个旋转张量和一个纯变形张量，从变形梯度中提取旋转场R，从而获得旋转后的全局刚度矩阵RK。

其中，所述碰撞检测具体为：假设一个广义的参数曲面

所述参数曲面

是由模型的所有离散点通过近似函数定义的曲面，且

可以通过Λ上的离散点近似出来，广义的参数曲面可以理解为模型参数曲面在边界上向四周的无限扩展；

确定碰撞点W是否存在，计算W点的坐标；

在广义的参数曲面

上，判断W点是否在由Λ定义的布的表面上。

其中，将采样点的初始位置和各个时间步的位移组合成布类仿真数据。

附图说明

图1为本发明的基于无网格的布类仿真方法的步骤框图；

图2是KL薄壳模型以及无网格方法的局部子域的示意图；

图3示出了寻找M点的方法的示意图；

图4示出了最典型的在圆筒上构建布类的示意图。

具体实施方式

布类动画模型经历了从纯几何阶段到基于力学模型的阶段。其主要应用的模型包括离散的质点弹簧模型和连续的有限元以及有限差分模型。离散的质点弹簧模型应用的技术已经开始十年左右，也能够取得比较逼真的动画效果，但是由于其对弹簧参数的依赖，很难获得指定材料的力学特性，这给应用带来本质上的困难，例如不同材料的服装获得的动画效果不应相同。

因此，连续模型是发展的必然。目前已经使用的连续模型主要是简化的有限元模型，但由于有限元模型中网格本身的限制，以及布类大变形的固有特征，使得有限元模型得到的动画与网格划分有很大影响，结果差强人意。而无网格法不需重新划分网格，具有仿真布类的天然优势。但是由于无网格模型本身的技术时间不长，应用方面的技术更是少之又少，目前还没有相关技术的无网格布类应用。这使得无网格布类的模拟技术变得十分有价值。

质点弹簧模型是现有技术中使用较多的一种离散模型。质点弹簧模型可细分为基于力的和基于能量的两种。基于力的质点弹簧模型是从每一质点的受力分析入手，根据物体上每一粒子的运动要满足牛顿第二定律。基于能量的方法是从整体上最小化系统能量从而获得下一个时间步系统的状态。

质点弹簧模型最早是用来求解静态问题，之后有学者将这种模型应用到皮肤和肌肉的建模上。Breen提出了粒子系统来建模布类，它将布类看成了由经线和纬线组成的纤维，他使用了Kawabata Evaluation System来获得测量数据，使用能量函数来建模弯曲、拉伸和剪切力，可以获得较为真实的布类运动。

Xavier提出了质点弹簧模型，他的思想是将模型视为由m×n个质点组成的网格，质点之间用无质量，自然长度不为零的弹簧连接。每一个质点和每一条弹簧分别处理，最后获得系统的整体运动。他将质子间的弹簧的作用分成几种类型，包括

●结构弹簧：连接两个质点，处理张力。

●剪切弹簧：连接对角线的质点，处理剪切力。

●弯曲弹簧：连接相隔的两个质点，处理弯曲力。

在Xavier的基础上，Choi和Ko将弯曲力和压缩力用一个单独的非线性模型代替，并且增加了对角线方向跨度为二的弹簧，以模拟大的弯曲，同时采取了两种不同的刚度，较大的用于拉伸和剪切，较小的用于弯曲和压缩，这样增加了系统的稳定性。

质点弹簧的不足之处在于网格的拓扑结构很大程度上决定了系统的行为，因而若拓扑结构设计不合理就会导致失真。如果使用细密的网格，得到的结果更精确，却导致了计算量的增加。如果选择了不正确的参数值就会导致系统不稳定，并且有一些材料的属性也没有在模型中得到很自然的表达。此模型也不适合表达刚度大的材料，会产生刚化问题。在应用中这种模型的主要瓶颈还在于细密的网格计算代价过高，因此很多模拟工作也集中到如何提高计算效率，改进计算方法上面。如Grzeszczuk将神经网络方法应用到基于物理模型的动画中，通过预先学习动力学模型，然后利用已经学习的模型进行模拟。该方法能生成接近物理真实的动画，但计算速度却快一至二个数量级。

针对布类仿真动画，比较突出的几项工作包括Baraff和Witkin提出将隐式积分方法引入布类的模拟，提高了系统的稳定性。它采用的是一种混合的方法，即基于能量函数构建系统方程，同时采用三角形网格模型，并将边和顶点作为物理模型中的弹簧和质子对待。在其能量函数计算中使用了稀疏的Jacobian矩阵，以加速计算。能量函数定义在有限三角形上，每个三角形平面内的变形能量从可以从力学方程得出，弯曲变形能量基于相邻三角形的角度。另外，在他们的方法中也使用了新的隐式时间步长求解数值计算方法，使他们的方法速度得到提高，他们的方法被应用到Maya软件的布类制作中。

Bridson提出了混合的显式隐式(implicit-explicit IMEX)的数值计算方法，同时采用了碰撞检测预测的技术，通过预处理的方法很好的得到了大型褶皱的布类动画。Boxerman和Ascher使用Choi的物理模型，考虑了局部性，将质点网格分解，并单独进行处理的改进方法。他们采用了自适应的显式隐式数值计算格式，提高了在每一个时间步内的效率。此模型使用了网格的有限差分方法进行离散化，在每一个时间步产生了一个稀疏的系统，此系统使用了Gauss-Seidel方法求解。MIRALab也开发了相应的系统，使用了弹性模型和碰撞响应。该方法主要特点是使用了连续介质理论和有限差分方法，将其应用在了旗子和下垂的实验中。

在此基础上，Aono考虑了约束的使用，将布类按照杨氏模量(Young Modulus)线性弹性模型进行建模，同时融入了其在平衡状态使用D′Alembert准则的Poisson系数。这种布类模型的假设是布类永远不会在其表面的法向变形。这种模型的优点是可以精确的表达褶皱在一个简单布类表面的变形传播，但是它不能够产生复杂的下垂运动，不能够处理碰撞。

有限元模型也同样被应用到了柔软布类的建模中。有限元作为一种在工程领域应用非常广泛的强大的连续模型，也被部分引用到布类仿真中。Eischen技术如何使用三角化布类表面来精确建立一个非线性的弹性模型，他使用了在纺织界应用的KES实验数据和曲线。此项技术重点在精确的建模布类的弯曲属性。使用有限元方法具有可以得到精确模型的优点，但是同时也具有求解速度慢的缺点。同时他的技术也只是限于几种有限的规律网格和材料，无法得到更多的推广，而三角化也是关键的瓶颈。

目前已有的布类仿真技术使用角度表达式来建模弯曲能量或弯曲力。Breen等使用线性梁理论来刻画弯曲时的曲率。曲率通过与之相应的两条特定的边和三个点来定义。Courshesnes等使用了相似的方法，但是使用了二面角来建模两个相邻三角形。Choi提出了同时计算材料平面压缩力和褶皱的弯曲模型。该方法假定，产生褶皱后的状态和相应的能量导致了弯曲力的产生。Bridson等确定了一个独立的弯曲模型，他的模型不会影响物体刚体的运动和in-plane变形。因此，他们获得了很多的与弯曲力相关的包含了2个邻居三角形的弯曲单元。Grinspun使用了离散曲率的方法来估计柔软的薄壳物体的弯曲能量。因为斜率的计算非常复杂，技术人员建议使用自动差分的方法。Wicke提出了使用点采样表面和薄壳动力学来进行快速动画。这种方法在发生切断和撕扯等拓扑结构需要较大变化的场景中非常适合。然而，这种方法由于使用了显式积分求解使其应用受到了限制，同时他的计算框架也不是十分明确。

以上所述传统的布类动画方法虽然能够产生动画，但是也存在着不可避免的问题，如碰撞问题，屈曲问题。同时具有本发明前述部分的一些固有缺陷。因此，提出本发明的无网格布类仿真方法。

首先介绍运动学描述：

Kirchhoff-Love(KL)模型假设需要求解的壳体是很薄的，这符合本发明模拟布类的需求。用全局笛卡尔坐标系E^I来表示三维空间中的壳。

定义了壳体上的任意一点的位置。同时，

定义为薄壳正中间位置的横截面上任意一点的位置。a₃是一个单位向量，作为薄壳表面的法向量。图3是KL薄壳模型的表示以及无网格方法的局部子域的示意图。Ω构型可以表示为：

其中，ξ¹，ξ²∈Λ，并且ξ∈<h^-，h⁺>。这里Λ表示参数空间，<h^-，h⁺>是参考构型中薄壳上、下表面的距离。

定义一组基向量g_I：

其中大写拉丁字母下标I表示从1到3。

在参数空间Λ中，本发明定义

$a_3=\frac{a_1\&times;a_2}{|a_1\&times;a_2|},---\left(3.2\right)$

其中希腊字母α的范围是从1到2。薄壳的a₃方向和中面的法向相同，并具有如下特性：

a_α·a₃＝0，|a₃|＝1.

在Ω构型中，本发明可以得到协变基向量为

$g_{\mathrm{\&alpha;}}=\frac{\&PartialD;x({\mathrm{\&xi;}}^1,{\mathrm{\&xi;}}^2)}{{\&PartialD;\mathrm{\&xi;}}^I}=a_{\mathrm{\&alpha;}}+\mathrm{\&xi;}{a}_{3,\mathrm{\&alpha;}},$

$g_3=\frac{\&PartialD;x({\mathrm{\&xi;}}^1,{\mathrm{\&xi;}}^2)}{\&PartialD;\mathrm{\&xi;}}=a_3.---\left(3.3\right)$

本发明的以下内容将使用上划线表示参考构型的量，例如

是指参考构型中表面上的一点；而没有上划线的量是指在当前构型中。

壳的Green-Lagrange应变张量定义为：

$E_{\mathrm{IJ}}=\frac{1}{2}(g_{\mathrm{IJ}}-{\stackrel{\&OverBar;}{g}}_{\mathrm{IJ}})$

$={\mathrm{\&alpha;}}_{\mathrm{IJ}}-{\mathrm{\&xi;\&beta;}}_{\mathrm{IJ}},---\left(3.4\right)$

张量α_IJ和β_IJ的非零项和壳的变形相关。

本发明改进了无网格法应变的共旋形式，这部分将在第(3.5)节中讨论。

所以本发明只需要导出线性的运动学方程：

$x({\mathrm{\&xi;}}^1,{\mathrm{\&xi;}}^2)=\stackrel{\&OverBar;}{x}({\mathrm{\&xi;}}^1,{\mathrm{\&xi;}}^2)+u({\mathrm{\&xi;}}^1,{\mathrm{\&xi;}}^2),$

其中u(ξ¹，ξ²)是薄壳中面的位移场。本发明可以将拉伸和弯曲应变表示为： ${\mathrm{\&alpha;}}_{\mathrm{\&alpha;\&beta;}}=\frac{1}{2}({\stackrel{\&OverBar;}{a}}_{\mathrm{\&alpha;}}\&CenterDot;u_{,\mathrm{\&beta;}}+u_{,\mathrm{\&alpha;}}\&CenterDot;{\stackrel{\&OverBar;}{a}}_{\mathrm{\&beta;}}),---\left(3.5\right)$

${\mathrm{\&beta;}}_{\mathrm{\&alpha;\&beta;}}=-u_{,\mathrm{\&alpha;\&beta;}}\&CenterDot;{\stackrel{\&OverBar;}{a}}_3$

$+\frac{1}{|{\stackrel{\&OverBar;}{a}}_1\&times;{\stackrel{\&OverBar;}{a}}_2|}[u_{,1}\&CenterDot;({\stackrel{\&OverBar;}{a}}_{\mathrm{\&alpha;},\mathrm{\&beta;}}\&times;{\stackrel{\&OverBar;}{a}}_2)+u_{,2}\&CenterDot;({\stackrel{\&OverBar;}{a}}_1\&times;{\stackrel{\&OverBar;}{a}}_{\mathrm{\&alpha;},\mathrm{\&beta;}})]$

$+\frac{{\stackrel{\&OverBar;}{a}}_3\&CenterDot;{\stackrel{\&OverBar;}{a}}_{\mathrm{\&alpha;},\mathrm{\&beta;}}}{|{\stackrel{\&OverBar;}{a}}_1\&times;{\stackrel{\&OverBar;}{a}}_2|}[u_{,1}\&CenterDot;({\stackrel{\&OverBar;}{a}}_2\&times;{\stackrel{\&OverBar;}{a}}_3)+u_{,2}\&CenterDot;({\stackrel{\&OverBar;}{a}}_3\&times;{\stackrel{\&OverBar;}{a}}_1)].---\left(3.6\right)$

从以上表达式本发明可以得出，薄壳中面的位移场u具备对薄壳变形的完全表达。所以根据以下关系，本发明将u作为分析中的主未知量：得益于Kirchhoff-Love的假设，薄壳上所有点的应变度量都可以从薄壳中面的变形得到。

本发明的无网格数值离散包括一下部分：近似方法、表面和位移近似、动力学系统方程、数值积分方案和时间积分方法。

其中，本发明的近似方法具体为：

在实现无网格方法时，需要使用无网格的近似，用来近似在求解域上的试函数。事实上，也通过各种不同的无网格近似函数来划分各种无网格方法。

在本发明中，使用了扩充的径向基函数(Augmented Radial basis functions，ARBF)近似来构成形函数。与目前应用较多的移动最小二乘(moving least squares，MLS)近似相比，径向基函数近似具有非常明显的优势。优势体现在以下方面，首先形函数拥有delta函数的特性，而且形函数的二阶导数更光滑；其次，径向基函数近似的计算量也比移动最小二乘近似减小很多。它的delta函数特性使得该方法可以直接的施加基本的边界条件。

图3是KL薄壳模型的表示以及无网格方法的局部子域的示意图。考虑一个点x的局部求解子域Ω_s，为了在域Ω_s上近似函数u的分布，需要通过一系列散布的点{x_i}，(i＝1，2，...，n)进行计算。无论是径向基函数近似还是移动最小二乘近似，其局部插值函数都可以表达为以下形式：

$u\left(x\right)={\mathrm{\&Phi;}}^T\left(x\right)u,\&ForAll;x\&Element;{\mathrm{\&Omega;}}_s,---\left(3.7\right)$

其中Φ(x)为形函数。

u(x)的导数可以使用形函数Φ_I(x)的导数来表示，

u_，j(x)＝{Φ_I，j(x)}^T{u_I}.    (3.8)

本发明以径向基函数插值作为例子：

Φ^T(x)＝[R^T(x)，P^T(x)]G，    (3.9)

其中R^T(x)＝[R₁(x)，R₂(x)，...，R_n(x)]是以n个离散点为中心的一系列径向基函数。P^T(x)＝[p₁(x)，p₂(x)，...，p_m(x)]是m阶的单项式基函数。G是由在散布点上的R和P组成的矩阵：

R₀＝[R^T(x_i)]，P₀＝[P^T(x_i)]，i＝1，2，...，n

由于近似函数构造的好坏直接影响到算法的精度和收敛性，在构造近似函数的时候，节点的分布，以及基函数(或核函数、权函数)的选择等问题也是需要注意的。关于节点的分布和数目，虽然无网格法不依赖于网格，但是节点分布会影响求解逆矩阵的条件数。所以必须产生足够的节点来满足单元矩阵(比如上述矩阵G)的可逆性。产生节点的方法有Biting Method，基于八叉树的节点生成技术，以及Voronoi图和加权球填充法等。Li和Liu给出了可行粒子分布(Admissible Particle Distribution)需要满足的条件。而在已知节点分布的情况下，可以使用适当的节点选取方法来防止产生奇异矩阵。

薄壳上的形函数Φ_I(x)至少需要C¹连续性，因此在对具有几何不连续表面建模时需要使用一些特殊技术来表达这种不连续性，如在布类表面产生折痕、缝纫线造成的折角等，这种表面需要在不连续处强加使用适当的边界条件。

本发明的表面和位移近似具体为：

在基于点的图形学中，类似于布类的物体可以通过曲面上散布的点的坐标进行渲染。本发明当前使用的方法是将布作为一种壳对待，这意味着本发明需要处理表面形状近似的问题。计算形函数的方法可以很自然地用来近似由一系列散布的点构成的表面。

径向基函数和移动最小二乘函数可以用来进行曲面估计。假设给出空间中m个散布的点，这些点在要估计的曲面上的位置为x_I。曲面的估计和位移的估计如下：

其中ξ¹，ξ²是曲面的参数坐标。注意为了估计位移向量u(ξ¹，ξ²)的导数，本发明只需要得到形函数Φ_I(ξ¹，ξ²)的导数。这意味着位移函数的光滑程度依赖于形函数。为了计算应变，本发明只需要获得对于参数ξ¹，ξ²的二阶导数。本发明的动力学系统方程具体为：

根据求解微分方程的离散形式的不同，无网格法可以分为强形式和弱形式。本发明可以通过加权余量法直接从平衡方程与边界条件得到无网格法的弱形式；也可以通过虚功原理或虚位移原理直接导出弱形式。

通过材料力学法则

S＝C(E).

将功共轭的一对应变E和应力S联系起来。

凡是物体几何约束所允许的位移就称为可能位移，取其任意微小的变化量就是虚位移δu，也就是几何上可能位移的变分。根据虚功原理或者虚位移原理：

${\&Integral;}_{\mathrm{\&Omega;}}\mathrm{\&delta;E}:\mathrm{Sd\&Omega;}-{\&Integral;}_{\mathrm{\&Omega;}}\mathrm{\&delta;u}\&CenterDot;\mathrm{qd\&Omega;}+{\&Integral;}_{\mathrm{\&Omega;}}\mathrm{\&nu;\&delta;u}\&CenterDot;\stackrel{\&CenterDot;}{u}\mathrm{d\&Omega;}+{\&Integral;}_{\mathrm{\&Omega;}}\mathrm{\&rho;\&delta;u}\&CenterDot;\stackrel{\&CenterDot;\&CenterDot;}{u}\mathrm{d\&Omega;}=0---\left(3.12\right)$

方程中的三项分别表示应力在虚应变上所做的功，以及由体积力、粘性力和惯性力在虚位移上所做的功。通过定义拉伸和弯曲应变和应力，应力在虚应变上所做的功可以重写为：

∫_ΩδE：SdΩ＝∫_Ω(δα^TH_mα+δβ^TH_bβ)dΩ，    (3.13)

其中H_m和H_b是和拉伸以及弯曲相关的材料参数的矩阵，这和有限元方法一致，更多关于H的细节可以参考(Cirak et al.，2000)。

因为公式(3.5)到(3.6)中的α，β和u具有线性关系，以及在公式(3.7)中，位移u的近似是由位移u_I计算获得的，本发明可以得到：

$\begin{array}{cc}\mathrm{\&alpha;}({\mathrm{\&xi;}}^1,{\mathrm{\&xi;}}^2)=\underset{I=1}{\overset{m}{\mathrm{\&Sigma;}}}{R}_m^I({\mathrm{\&xi;}}^1,{\mathrm{\&xi;}}^2)u_I,& \left(3.14\right)\\ \mathrm{\&beta;}({\mathrm{\&xi;}}^1,{\mathrm{\&xi;}}^2)=\underset{I=1}{\overset{m}{\mathrm{\&Sigma;}}}{R}_b^I({\mathrm{\&xi;}}^1,{\mathrm{\&xi;}}^2)u_I,& \left(3.15\right)\end{array}$

其中

和分别对应于拉伸和弯曲应变的节点位移。

以时间的二阶常微分方程表示的动力学方程如下：

$M\stackrel{\&CenterDot;\&CenterDot;}{u}+D\stackrel{\&CenterDot;}{u}+\mathrm{Ku}=f,---\left(3.16\right)$

其中M是对角节点质量矩阵，D是粘性矩阵，K是刚度矩阵，节点位移向量是u，力向量是f。

已上偏微分方程的离散形式为全局弱形式，而在局部Petrov-Galerkin法中，本发明可以写出一个点x_k的局部子域Ω_s上的弱形式，局部子域可以是任意的形状。应用公式(3.12)到(3.15)，对应于动力学公式(3.16)的刚度矩阵和节点力向量的弱形式可以表示为一种一般形式：

$K_k^{\mathrm{IJ}}={\&Integral;}_{{\mathrm{\&Omega;}}_s}[{\left(R_m^I\right)}^T{H}_m{R}_m^J+{\left(R_b^I\right)}^T{H}_b{R}_b^J]\mathrm{d\&Omega;},---\left(3.17\right)$

$f_k^I={\&Integral;}_{{\mathrm{\&Omega;}}_s}{\mathrm{\&Phi;}}_I\mathrm{qd\&Omega;}.---\left(3.18\right)$

以上推导结果是局部Petrov-Galerkin形式，它是由一般的Galerkin形式推广出来的，在局部子域满足平衡方程和边界条件。可见它的格式与有限元法的格式非常相似。

本发明的数值积分方案具体为：

一般采用全局积分方程形式的算法产生的整体矩阵为满阵，整体方程的求解将非常困难；而本发明采用的局部积分方程形成的矩阵为类似于有限元的稀疏矩阵，在计算效率上有明显的优势。局部刚度矩阵可以通过局部子域Ω_s的数值积分得到。由于局部子域可以是任意的形状，也就是说局部积分区域比较确定，往往为比较简单的规则形状，如圆、球形或这些图形的交叉部分，积分的收敛性比基于背景网格的积分方法好，因此局部弱形式的积分不需要任何网格或背景网格。一般情况下可以使用De和Bathe等提出的适用于圆域或球域的特殊高斯积分方案或分片中点积分方案。基于MLPG法导出的局部积分其在被积分区域为圆域或扇形，并且被积函数为非多项式。在这种情况下，分片中点积分方案在数值稳定性和精度上都优于高斯积分方案。

在计算机动画的解决方案中，数值积分方案是一个很灵活的步骤。根据不同的需求，可以采用不同的数值积分方案。比如，在实时性要求较高的交互式动画中，可以使用配点型积分方案。Liu提出了MLPG混合的配点法，这种方法具有稳定的收敛速度，而且比MLPG混合的有限体积法具有更高的求解效率。而在精度要求相对较高的影视动画等应用中，可以在分片中点积分方案或高斯积分方案，并使用较多的积分点。当然，使用各种积分方案的前提是保证计算的稳定性。

本发明的时间积分方法具体为：

对于以上动力学方程的数值解法，可以使用任何一种合适的时间积分方法，将公式(3.16)转换成一阶的常微分方程来求解。时间积分方法在布类模拟中已经获得了技术人员的重点关注。布类是一个高度的各向异性的材料，这是由于其材料属性决定：它对弯曲产生阻碍比较微弱，但是对于拉伸会产生很大的阻力。当模拟一个高硬度的体材料的时候，本发明可以使用一个精简的变形模型，例如使用模态分析的方法，将其视为一个坚硬的物体。显然，这对于布类不适用。本发明对于大变形的具有平面外的弯曲和褶皱的布类行为建模和可视化感兴趣，因此就要在本发明的模型中处理这些平面内能量。从数值角度，微分方程是有刚性的，也就是说它拥有很广范围的特征值。在使用显式方法时这种力一定要使用小步长的时间步进行。

Xavier使用了减小拉伸弹簧的能量，之后在每一个时间步的结尾处理布类的网格的方法来解决这个问题。约束弹簧的拉伸不能超过10％，使得邻居节点的拉伸被释放，因此可以获得收敛。事实上，这反映了布类运动的二阶段性，首先是很微弱的抑制力，之后达到一个阈值之后，刚度就会明显的增加。在应用中，这种方法给出了很高效，很真实的结果，但是这种问题的收敛性还没有被证实。

本发明使用一种线性的隐式格式以保证大的时间步的收敛特性。这种方法已经被证明是一种鲁棒而且高效的解决刚性问题的方法，是一种应用十分广泛的技术。特别的，它应用了后向欧拉方法，给出公式：

$\left[\begin{array}{c}\mathrm{\&Delta;x}\\ \mathrm{\&Delta;v}\end{array}\right]=\mathrm{\&Delta;t}\left[\begin{array}{c}{v}_0+\mathrm{\&Delta;v}\\ M^{-1}F(x_0+\mathrm{\&Delta;x},v_0+\mathrm{\&Delta;v})\end{array}\right]---\left(3.19\right)$

其中F_i＝f_i-K_iu_i-D_iu_i。这是一个对于Δx和Δv的非线性方程。通过一阶泰勒展开，可以线性化方程(3.19)：

$F(x_0+\mathrm{\&Delta;x},v_0+\mathrm{\&Delta;v})=F_0+\frac{\&PartialD;F}{\&PartialD;x}\mathrm{\&Delta;x}+\frac{\&PartialD;F}{\&PartialD;v}\mathrm{\&Delta;v}---\left(3.20\right)$

在此方程

中，是为了估计状态(x₀，v₀)，同

理可得。将此泰勒展开带入方程3.19，可以获得线性方程组

$\left(\begin{array}{c}\mathrm{\&Delta;x}\\ \mathrm{\&Delta;v}\end{array}\right)=h\left(\begin{array}{c}{v}_0+\mathrm{\&Delta;v}\\ M^{-1}(F_0+\frac{\&PartialD;F}{\&PartialD;x}\mathrm{\&Delta;x}+\frac{\&PartialD;F}{\&PartialD;v}\mathrm{\&Delta;v})\end{array}\right)---\left(3.21\right)$

由Δx＝h(v₀+Δv)，可以获得

$\mathrm{\&Delta;v}=h{M}^{-1}(F_0+\frac{\&PartialD;F}{\&PartialD;x}h(v_0+\mathrm{\&Delta;v})+\frac{\&PartialD;F}{\&PartialD;v}\mathrm{\&Delta;v})---\left(3.22\right)$

使用I来表示单位矩阵，方程重新写成：

$(I-h{M}^{-1}\frac{\&PartialD;F}{\&PartialD;v}-h^2{M}^{-1}\frac{\&PartialD;F}{\&PartialD;x})\mathrm{\&Delta;v}=h{M}^{-1}(F_0+h\frac{\&PartialD;F}{\&PartialD;x}{v}_0)---\left(3.23\right)$

这样我们获得了Δv的解，给定Δv，可以很方便的计算出Δx＝h(v₀+Δv)。

尽管求解此方程增加了每一个时间步的计算代价，但是由于可以使用相对较大步长，此方法还是具有一定的效果。而这种方法也被广泛的被应用到布类模拟器的求解中。此外，近年来随着计算技术的发展，结合显式隐式方法的计算也同样被应用到布类的模拟中。

以节点速度为主要未知量生成的线性方程组一般规模很大，而且通常是稀疏的方程组，可以直接通过直接分解的方法或者迭代分解来求解。

以上方法允许本发明使用灵活的时间步长，并且同时保证稳定性。然而，有些时候本发明仍然需要减小时间步长来避免发散。大多数改变时间步长的方法，不管是针对显式时间积分还是隐式时间积分的，都是以提高模拟的准确性为目的。而本发明的研究目标是计算机动画，相比于工程计算，本发明更希望得到视觉上更好的结果，也就是说，数值稳定性比准确性更重要。数值稳定性是本发明的时间积分方法的决定性因素，所以需要在不稳定因素出现之前将它识别出来。

计算中任何潜在的不稳定因素都来源于刚度，而它几乎完全是由布类的强拉伸力产生的。在一个隐式时间步结束后，本发明会得到一个位移差Δx，作为这个时间步系统状态的改变参量。对每一个节点，本发明检查更新后的系统状态中的拉伸项，如果任何一个节点的拉伸过大，本发明则放弃该时间步，退回到上个时间步的状态，然后缩小时间步长，并重新计算。用这种方法，很容易提前发现不稳定的状态并避免它。

系统中有一个代表可容忍的最大时间步长的变量，这个变量可以由用户设定，而且应该总是小于或等于产生的计算机动画的帧步长(一般来说，帧步长设为0.02秒)。当计算中时间步长被缩小之后，如果成功计算超过两步，则系统尝试增加时间步长。如果大步长的模拟失败，则系统继续维持小步长，并且提高维持小步长的步数。这种方法在实践中证明是可行的。

本发明的无网格大变形共旋形式具体为：

当考虑有限应变时，应变的度量为非线性，即几何非线性，这时，模拟求解的方程变化为非线性系统。不幸的是，非线性系统的求解代价很高。应变的几何非线性和有限应变的旋转不变性相关。线性柯西应变不具备旋转不变性，而在本发明的布类模拟中，薄壳模型使用Green-Lagrange应变，它对位移是非线性的。

Green-Lagrange应变会产生一个非线性的偏微分方程。而共旋形式(corotational formulation)的目的就是削除几何非线性。共旋形式在每个节点上都具有局部坐标系，共旋形式的思路是跟踪物体上的每个节点的局部坐标的旋转分量。在计算机图形学领域，Muller使用方法导出和求解共旋形式。而本发明通过极分解的形式，虽然计算开销稍大一点，但是可以获得更好的稳定性和准确性，而稳定性和准确性在计算机动画的多时间步迭代的动力学计算中是非常重要的。而且，本发明的共旋形式是针对无网格方法的，这与其他文献中基于有限元方法的共旋形式有所不同。

本发明对方程(3.5)和(3.6)中的应变进行线性化，由此使得在前面的运动学描述的应变对位移是线性的，但是不再具备旋转不变性。然而，如果旋转场R已知，本发明就可以使用共旋的应变形式，并获得在旋转后的当前构型上的线性应变张量.

由方程(3.1)和方程(3.2)，变形的梯度可以写为：

$F=\frac{\&PartialD;x}{\&PartialD;\stackrel{\&OverBar;}{x}}=\frac{\&PartialD;x}{\&PartialD;{\mathrm{\&xi;}}^I}\&CircleTimes;{\stackrel{\&OverBar;}{g}}^I=[a_{\mathrm{\&alpha;}}+\mathrm{\&xi;}{a}_{3,\mathrm{\&alpha;}}]\&CircleTimes;{\stackrel{\&OverBar;}{g}}^{\mathrm{\&alpha;}}+a_3\&CircleTimes;{\stackrel{\&OverBar;}{g}}^3.---\left(3.24\right)$

通过极分解，变形梯度张量F可以分解成一个旋转张量R和一个纯变形张量U：

F＝RU

本发明使用一种高效的二次收敛的迭代方法，可以从变形梯度中提取旋转场：

R⁰：＝F    (325)

$R^{n+1}:=\frac{1}{2}(R^n+{\left(R^n\right)}^{-T}).---\left(3.26\right)$

对于计算R，这是一种很高效而且准确的方法，同时可以通过迭代的次数来控制精度。

获得了旋转场R，刚度矩阵Ku变成了RKu^R。从节点I的角度考虑，对节点J的位移影响是

因此，动力学系统方程(3.16)可以重新写为

$M\stackrel{\&CenterDot;\&CenterDot;}{u}+D\stackrel{\&CenterDot;}{u}+{\mathrm{RKR}}^{T}u=f+\mathrm{RK}(R^T-I)(x_{\mathrm{node}}-x).---\left(3.27\right)$

这样，将全局非线性问题转化为了局部线性问题。每个时间步需要求解的系统方程可以保持线性。和传统的线性无网格方法相比，线性系统不是固定的，当参考构型更新时，线性系统随之改变。而且，可以使用原始构型作为参考构型(区别于Updated Lagrange格式)。共旋形式的动力学方程相当于将无网格计算节点与其相关的周围节点所构成的假想“单元”，按照求得的局部旋转坐标反向旋转。对于每个“单元”，应变是线性的；而对于全局区域，应变是非线性的。

本发明的无网格表面的碰撞检测具体为：

对于碰撞检测，由于纯的无网格方法不保留表面的网格，因此传统的基于网格的碰撞检测方法就会失去作用。比如，广泛应用的点面检测和面面检测，在无网格法中无法使用。而且，由无网格节点构成的参数曲面的碰撞检测，还不同于普通的三维无网格固体的碰撞检测技术。尽管后者在无网格接触问题的技术中已经存在几种解决方法，但这些方法只适用于普通的三维无网格固体，在本发明的无网格参数表面的问题中并不适用。

本发明提出了一种新的碰撞检测方法，它适用于由点构成的无网格表面，是一种基于形函数的碰撞检测方法。这种方法用来判断在一个无网格离散的表面(如无网格布类)上，一个节点是否会与的表面的其它部分发生碰撞。

首先，其中所述碰撞检测具体为：假设一个广义的参数曲面

所述参数曲面

是由模型的所有离散点通过近似函数定义的曲面，且

可以通过A上的离散点近似出来，广义的参数曲面可以理解为模型参数曲面在边界上向四周的无限扩展。举例说明，假设在一块布片中间裁剪出一个圆孔，

那么全局参数曲面是包括这个圆孔上的中面的近似曲面的。

以在位置为P的节点为例，设在下一个时间步，它的位置是Q。如图3所示，如果P和广义的参数曲面

具有一个碰撞，那么PQ就会穿过

假设其交点为M。

这样问题就可以分成两个步骤解决，首先本发明决定碰撞点M是否存在，以及计算M点的参数坐标；第二步，在广义参数曲面

上，判断M点是否在由Λ定义的布的表面上。注意到，并不能从广义参数曲面上M点的参数坐标来直接判断M点是否也在布片模型的参数曲面上。

为了获得M，本发明需要将P点映射到曲面上。本发明假设P′点是P点的映射，满足方程(3.1)，

$P=P^{\&prime;}({\mathrm{\&xi;}}^1,{\mathrm{\&xi;}}^2)+\mathrm{\&xi;}{a}_3^P({\mathrm{\&xi;}}^1,{\mathrm{\&xi;}}^2),---\left(3.28\right)$

P′不能够直接映射，首先本发明需要找到距离最近的节点X，通过以下的评价标准：

$-\mathrm{\&Delta;}<\frac{\mathrm{XP}\&CenterDot;a_3}{\left|\mathrm{XP}\right|}<+\mathrm{\&Delta;},---\left(3.29\right)$

其中Δ是一个很小的值，a是在方程(3.2)中定义的X的局部坐标，

这个评价标准表示，PX垂直于X的局部子域(对于布类的例子，局部子域是一个曲面)。

因此P′和X邻近，而且P′点的参数坐标可以写为

$({\mathrm{\&xi;}}_X^1+\mathrm{XP}\&CenterDot;a_1,{\mathrm{\&xi;}}_X^2+\mathrm{XP}\&CenterDot;a_2)---\left(3.30\right)$

之后本发明可以很容易的获得以下结论：

如果XP·a₃与QP·a₃同号，同时|XP·a₃|≥|QP·a₃|(注意到|XP.a₃|＝|MP.a₃|)，那么PQ就与全局参数曲面有一个碰撞点M，同时

$M=P+\mathrm{PQ}\frac{\left|\mathrm{PM}\right|}{\left|\mathrm{PQ}\right|}=P+\mathrm{PQ}\frac{\mathrm{XP}\&CenterDot;a_3}{\mathrm{PQ}\&CenterDot;a_3}---\left(3.31\right)$

因此M的参数坐标就可以表达为

图3给出了上述的寻找M点的算法的示意图。

一旦本发明获得了M的参数坐标，第二步就简化为，在二维空间中无网格曲面的碰撞问题。

在本发明的无网格插值方法中，形函数(3.9)的构建需要计算矩阵

$N=\left[\begin{array}{cc}{R}_0& P_0\\ P_0^T& 0\end{array}\right]---\left(3.32\right)$

技术人员已经发现这个矩阵具有跟踪连续物体表面的指示作用。基于从形函数得到的矩阵N，二维参数空间碰撞检测可以准确而容易的获得。Li证明一个连续区域的内部和外部能够通过计算矩阵N的行列式来区分。在任何形状的连续区域的内部，det{N(x)}具有正值，而在域外部，det{N(x)}→0或为负值。上述这种非常有用的属性让本发明能够精确的跟踪任何形状连续体内的任何位置，而不需要提取它的边界。这在普通的方法中是基本不可能获得的。因此本发明可以判断碰撞点M是在布的区域内，或者是在剪裁掉的区域内。所以，尽管本发明不能直接从全局参数曲面上M点的参数坐标来判断M点是否也在布片模型的参数曲面上，但是通过计算矩阵N的行列式，本发明达到了这个目的。

碰撞相应具有几种不同的方法。物理上，它可以通过添加惩罚力或在接触点间添加约束来实现。幸运的是，这些碰撞响应方法都是和基于网格的方法相同的。碰撞点M的a₃(通过点P的位置在a₃和-a₃之间选择)可以用来指示碰撞相应的方向。

PQ穿透的比例，可以用来指示惩罚力的大小或约束的位置。

本发明的约束处理具体为：

本发明的一些实施方式提供给节点强加约束的方法。这里本发明所说的约束是指：用户指定的约束(如，固定一个节点，或强制一个节点运动的轨迹)，以及由模拟中的背景物体(在模拟中作为刚体，如人体模型)和节点之间的接触产生的约束。

在一个计算时间步中，一个节点的状态可以是完全不受约束的，或者受到一个方向的约束(只能在一个平面上运动)，或者受到两个正交方向上的约束(只能在一条直线上运动)，或者受到三个正交方向上的约束(固定不动)。比如，一个节点和一个刚体表面接触的时候，它只能在这个刚体表面上滑动，这时该节点受到刚体表面法向上的一个方向的约束。

传统的约束方法可以总结为三种：缩减坐标系法、罚函数法以及拉格朗日乘子法。但是它们并不符合本发明的要求。

缩减坐标系法是一种非常直接和准确的方法，它直接减少受约束节点的坐标数目来达到约束的目的。完全约束的节点是没有坐标的，而受一维约束的节点具有两维坐标。这种方法只能在局部坐标系的方向上施加约束，对布类和刚体接触的约束来说，这是不够的。而且坐标系的缩减会让系统变得混乱，如果本发明更改所有节点的坐标数量，那么必须改变方程中导数矩阵，以及形成的系统刚度矩阵的结构。考虑到这些，缩减坐标系法是不适用的。

罚函数法通过将基本的边界条件约束用罚函数引入偏微分方程弱形式。罚函数法在理论上是可行的，在实际计算中的缺点是罚因子的取值难于把握，太小起不到惩罚作用，太大则由于误差的影响会导致错误。罚函数法施加约束的结果并不准确，在计算机动画中会产生“穿透”现象等非常敏感的视觉错误。而且，罚函数法给系统带来了附加刚度，使系统稳定性变差。

拉格朗日乘子法利用拉格朗日乘子将约束条件加入弱形式中。该方法增加了拉格朗日系数，从而增加了整体方程的数目。对于动力学问题或非线性问题，生成增大的方程组会耗费更多的计算复杂度，而且生成的新矩阵还可能非正定，可能引起算法稳定性问题。

本发明改进了Baraff提出的一种通过修改质量矩阵来达到控制约束目的的方法。虽然该方法的提出是针对质点弹簧法产生的动力学方程的，本发明将其改进并应用到无网格法的弱形式上。这种方法应用了质量矩阵的逆M^-1。考虑节点i，由方程3.16得

${\stackrel{\&CenterDot;\&CenterDot;}{u}}_i=\frac{1}{m_i}{F}_i$

其中 $F_i=f_i-K_i{u}_i-D_i{\stackrel{\&CenterDot;}{u}}_i$

假设本发明需要约束节点i的速度不变，那么本发明可以设

为零，即给节点i一个无穷大的质量值，使得方程忽略所有施加在节点i上的力。对节点加速度的完全约束可以通过设置该节点的质量的逆为零，但是本发明希望可以只在一个方向或两个方向上约束节点。当然，约束方向可以是任意的，并不局限于局部坐标系的方向。将质量写成矩阵的形式，动力学方程可写为，

${\stackrel{\&CenterDot;\&CenterDot;}{u}}_i=\left(\begin{array}{ccc}\frac{1}{m_i}& 0& 0\\ 0& \frac{1}{m_i}& 0\\ 0& 0& \frac{1}{m_i}\end{array}\right)F_i$

并给定一个单位向量p_i∈R³，因为

所以可以用

$\frac{1}{m_i}(I-p_i{p}_i^T)$

代替方程中质量矩阵的逆，来实现对节点i在p_i方向加速度的约束。同样，给定两个单位向量p_i和q_i，可以用

$\frac{1}{m_i}(I-p_i{p}_i^T-q_i{q}_i^T)$

代替方程中质量矩阵的逆，来实现对p和q两个方向加速度的约束。这种方法相当于除去约束方向在质量矩阵的逆矩阵上的投影。设ndof(i)为节点i的未受约束的自由度，修正的质量矩阵的逆设为W，在对角元上，

本发明有

$S_i=\left\{\begin{array}{cc}I& \mathrm{ndof}\left(i\right)=3\\ (I-p_i{p}_i^T)& \mathrm{ndof}\left(i\right)=2\\ (I-p_i{p}_i^T-q_i{q}_i^T)& \mathrm{ndof}\left(i\right)=1\\ 0& \mathrm{ndof}\left(i\right)=0\end{array}\right.$

对于节点i，设z_i为本发明期望的施加约束后在约束方向上速度的改变量，将修正的质量矩阵的逆应用到由时间积分方案得出的线性方程组3.23中，得

通过z本发明可以量化的约束一个节点。

上式得出约束后的速度改变量□v，对于完全约束的节点，□v_i＝z_i；对于部分约束的节点，约束部分的加速度分量等于z_i。

本发明的表面参数化和布类边界缝纫具体为

必须定义方程(3.1)中的参数空间Λ，以直接使用建立在参数($ξ¹，ξ²$)的上的形函数。对于由分块的布缝纫成的服装的参数化，从CAD数据中获得每一张单独的布表面上沿着经线和纬线两个方向进行参数化，是最直接而且简单的方法。所有的从CAD数据中获得的布片都需要组合在一起而形成服装，而这意味着在原始的布片是简单的形状，并且不存在自身交叉。

这样，每一张布片的参数化都很容易。

另外一个使用经线和纬线作为参数坐标的优势就是，由于通过这两个方向进行参数化可以很方便的建立在经线和纬线上各向异性的布。通过这种参数方法，参数坐标就可以很容易的计算出来。

单一布片的展平模式代表了布类的自由静止状态。一件织物可能由多块这样的单一布片组成，所以构造服装模拟的第一步是确定这些布片“穿”在3D人体模型上的相对位置关系，

并将它们按照这个位置关系缝纫到一起。图4显示了最典型的在圆筒上构建布类的例子，这需要将四边形的两个边缝纫在一起。

在将服装缝纫起来之后，本发明的模拟器需要连接各块布的缝纫线，以保持缝纫在一起的节点在位移上连续性。拉格朗日乘子(Lagrange Multiplier)可以用来约束缝纫点具有相同的位移。但是正如前面所提到的，尽管已经约束到一起，每一个布片仍然需要单独进行参数化。

同样，每一块布片上的形函数也只在当前布片上计算，而不会跨越到链接在一起的邻居布片上，即使是在缝纫边界处也是如此。全局的刚度矩阵可以通过整体服装的各个部分进行组合。

因为拉格朗日乘子法扩大了全局刚度矩阵，本发明只利用拉格朗日乘子法来约束缝纫边界，因为缝纫边界的约束条件是一种隐式的约束处理，需要在方程求解之后才能得到约束的显式条件。而对于另一种类型的约束，即由基于约束处理碰撞响应所产生的约束条件，本发明采用了Baraff中提出的约束方法。此方法对刚度矩阵K和速度矩阵D进行了滤波，是一种显式的约束处理，它不会增加刚度矩阵的大小。

相比于质点弹簧法和有限元法等基于网格的方法，本发明的无网格法在裁剪面的不规则边界的处理上具有很大优势，它不用破坏和重建网格。而对于质点弹簧法，边界上的不规则网格给计算带来很大问题；对于有限元法，网格的重新划分是非常困难的问题。

Claims (6)

1.一种基于无网格的布类仿真方法，包括：

(a)将给定的布类模型Ω参数化到中面Λ上，确定初始构型和参考构型；

(b)在Λ中确定有限数目的节点及其局部子域Ω_s；

(c)在每个节点的定义域上计算形函数Φ_I，计算形函数的导数；

(d)对所有全局区域内的节点循环步骤(d1)和(d2)：

(d1)确定节点的子域Ω_s上的所有高斯积分点X_Q；

(d2)对所述所有高斯积分点x_q计算数值积分，组装全局刚度矩阵K；

(e)组装节点质量矩阵M和粘性矩阵D；

(f)根据仿真时间长度t和所设置的时间步间隔δt确定循环次数，将步骤(f1)至(f8)循环t/δt次：

(f1)计算旋转后的全局刚度矩阵；

(f2)计算强加约束后的全局刚度矩阵、节点质量矩阵和粘性矩阵；

(f3)根据不同的时间积分方法合成线性方程组；

(f4)求解生成的线性方程组，并获得节点位移；

(f5)更新采样点的未知量和它们的导数；

(f6)碰撞检测和响应，更新节点位移；

(f7)提取新的旋转场R；

(f8)重新采样得到采样点的位移。

2.根据权利要求1所述的方法，其特征在于，所述将给定的布类模型Ω参数化到中面Λ上具体包括：将布类模型分成多个单独的布片，确定所述多个布片的相对位置关系，按照相对位置关系将所述多个单独的布片组合在一起，从而保证所述多个单独的布片具有位移上的一致性，并沿经线和纬线两个方向对上述单独的布片进行参数化。

3.根据权利要求1所述的方法，其特征在于，其中形函数Φ_I具体为：

Φ_I＝[R^T(X)，P^T(X)]G。

其中R^T(x)＝[R₁(x)，R₂(x)，...，R_n(x)]是以n个离散点为中心的一系列径向基函数。P^T(x)＝[p₁(x)，p₂(x)，...，p_m(x)]是m阶的单项式基函数。矩阵G是由离散点上的R和P组成的矩阵，

$G={\left[\begin{array}{cc}{R}_0& P_0\\ P_0^T& 0\end{array}\right]}^{-1},$

R₀＝[R^T(x_i)]，P₀＝[P^T(x_i)]，i＝1，2，...，n

4.根据权利要求1所述的方法，其特征在于，其中所述计算旋转后的全局刚度矩阵K具体为：通过旋转场使用共旋的应变形式，获得旋转后的线性应变张量，通过极分解将变形梯度张量分解成一个旋转张量和一个纯变形张量，从变形梯度中提取旋转场R，从而获得旋转后的全局刚度矩阵RK。

5.根据权利要求1所述的方法，其特征在于，其中所述碰撞检测具体为：假设一个广义的参数曲面所述参数曲面

是由模型的所有离散点通过近似函数定义的曲面，且

可以通过Λ上的离散点近似出来，广义的参数曲面可以理解为模型参数曲面在边界上向四周的无限扩展；

确定碰撞点W是否存在，计算W点的坐标；在广义的参数曲面上，判断W点是否在由Λ定义的布的表面上。

6.根据权利要求1所述的方法，其特征在于，将采样点的初始位置和各个时间步的位移组合成布类仿真数据。

Images