CN102393968A - 一种几何精确的不可拉伸条带模拟方法 - Google Patents

Abstract

本发明公开了一种几何精确的不可拉伸条带模拟方法。使用一类特殊的直纹面来描述不可拉伸条带对应的二维曲面，结合薄壳模型得出了一套用条带中心线描述的能量模型；与传统的将薄壳模型直接应用于二维曲面并进行三角或四边形离散的方法相比，算法复杂度大大降低，提高了求解效率。通过使用“四元数-中心线节点坐标-饶率曲率比”这一广义坐标来描述条带中心线，得到的运动方程能够稳定和高效地求解，外部力和力矩的施加也非常方便。最后通过利用条带扭转部分形变的频率远大于弯曲部分的频率，简化运动方程的求解为准静态地更新扭转部分的形变和动态更新弯曲部分的形变两步，这一处理在几乎不降低动态轨迹精度的同时进一步减少了计算开销。

Description

一种几何精确的不可拉伸条带模拟方法

技术领域

本发明涉及动画制作、工业生产中基于物理的弹性体模拟方法，尤其涉及一种几何精确的不可拉伸条带模拟方法。

技术背景

不可拉伸条带状物体特指在几何形状上介于线缆和布料之间的物体，即厚度可忽略，宽度远小于长度的物体；而在应变方面变现为易于弯曲和扭转，但极难拉伸或压缩。在宏观世界中不可拉伸条带状物体可以是装饰带、纸条、皮带甚至电影胶片，他们时常出现于电影，游戏或动画之中；而在微观世界，条带状物体也比比皆是，如石墨烯形成的条带，晶体生产形成的条带，DNA或者一些生化聚合物。一个精确并且高效的模拟方法，是艺术特效，游戏开发，工业生产，科学研究相关人员所希望拥有的工具。

目前可以被用于模拟不可拉伸条带的潜在方法主要是两种，一种是线缆类物体的模拟方法，另一类则是布料类物体的模拟方法，但这两种方法各有其弊端。

在这里线缆类模拟方法特指基于Cosserat或Kirchhoff-Love线缆模型的模拟方法。该模型仅使用线缆的中心线及中心线的材料标架表示线缆，并利用中心线和材料标架表示线缆对应的沿法向方向的弯曲势能，沿副法向方向的弯曲势能以及扭转势能，与体积或表面表示线缆的方法相比计算维度大大降低。代表作有(a)SPILLMANN，J.，AND TESCHNER，M.2007.CoRd E：Cosserat rodelements for the dynamic simulation of one dimensional elastic objects.InProceedings of the 2007 ACM SIGGRAPH/Eurographics symposium on Computeranimation，Eurographics Association，72.(b)BERGOU，M.，WARDETZKY，M.，ROBINSON，S.，AUDOLY，B.，AND GRINSPUN，E.2008.Discrete elastic rods.ACM Transactions on Graphics 27，3.(a)使用中心线加四元数的表示方法表示线缆，(b)则使用“曲线-角度”的表示方式。通过调节法向与副法向的弯曲势能比，该模型可近似条带状物体的弯曲弹性响应，然而不可拉伸条带曲面“可展”这一最主要的几何性质却无法被保留，使用该模型模拟经过扭转的条带，条带两侧被拉长。

另一类基于Kirchhoff-Love薄壳模型的通常用于布料类物体的模拟，代表作为GRINSPUN，E.，HIRANI，A.，DESBRUN，M.，AND SCHRODER，P.2003.Discrete shells.In Proceedings of the 2003 ACM SIGGRAPH/Eurographicssymposium on Computer animation，Citeseer.该方法将布料视为一个二维曲面，利用曲面曲率信息和弹性势能的关系，模拟布料类物体。尽管结合合适的长度约束处理策略该类方法能够保证不可拉伸条带曲面“可展”，然而该类方法在计算时需要对二维曲面进行三角或四边形离散，模拟条带将不必要地引入大量离散单元，使得计算效率低下。

发明内容

针对背景技术的不足，本发明的目的在于提供一种几何精确的不可拉伸条带模拟方法，该方法不仅具有Cosserat或Kirchhoff-Love线缆模型计算维度低的特点，同时又具备Kirchhoff-Love薄壳模型保证不可拉伸条带可展这一几何性质。

为实现上述的目的，本发明采用的技术方案它包括以下八个步骤：

(1)使用特殊的直纹面解析地描述不可拉伸条带对应的二维曲面，即不可拉伸条带形状由其中心线描述；

(2)根据步骤(1)建立不可拉伸条带上各点处的平均曲率与其中心线上各点的曲率、饶率以及饶率曲率比等的关系，结合薄壳模型，导出由中心线描述的弹性势能；

(3)根据步骤(1)建立不可拉伸条带上各点处的速度与其中心线上各点的速度、直纹变化速度的关系，导出由中心线描述的动能；

(4)“四元数-中心线节点坐标-饶率曲率比”这一广义坐标来描述根据步骤(1)中条带中心线，由该广义坐标导出的步骤(2)中描述的弹性势能、步骤(3)中描述的动能以及保证直纹面可展、广义坐标冗余自由度一致的约束数值无奇异；且易于用户操作相关约束的施加；

(5)借助有限元并利用分段线性基函数在空间上对步骤(4)中描述中心线的广义坐标进行离散化，导出离散可计算的弹性势能，动能以及约束；

(6)基于条带扭转部分形变的频率远大于弯曲部分的频率，分步求解由步骤(5)中离散能量构成的运动方程：首先准静态地更新扭转部分的形变，然后动态更新弯曲部分的形变。

(7)在步骤(6)求解过程中，利用运动方程分步求解这一特点同时借助时间域平行传输，将控制条带扭转的朝向约束转化为求解更为简单标架约束；

(8)在步骤(6)求解过程中，通过定义合适过渡函数，本条带模拟方法能够拓展到可定向或不可定向环状条带的模拟。

步骤(1)中所述的特殊的直纹面是一种保证曲面可展的直纹面，是对rectifying developable的拓展，提供了条带形变过程几何上精确的前提；rectifyingdevelopable由三维空间中的正则曲线描述，而平直或者波浪状的条带存在曲率为零的区域，即其中心线可能为非正则，因此rectifying developable被拓展以描述平直或者波浪状的条带，具体指对直纹的重新定义，使得在曲率为零的区域描述直纹所需的中心线法向量和饶率曲率比有定义；对于中心线法向量，始终使用条带一侧的在中心线处的表面法向；此外要求曲率为零区域的饶率曲率比平方和最小。

步骤(2)中所述的由薄壳模型推导所得的弹性势能中包含了一个指数项，当直纹相交于条带曲面内时该指数项将弹性势能惩罚至无穷，即弹性势能蕴含了“直纹不相交于条带曲面内”这一硬约束；因而在直纹相交状态下数值上奇异或无定义；为使得数值能够计算，指数项被泰勒展开，而直纹相交于条带曲面内这一原本蕴含于弹性势能的硬约束通过后期处理来满足；弹性势能形为：

$V=\underset{a=0}{\overset{\∞}{\mathrm{\Σ}}}\frac{w^{2a+1}D}{2^{2a+1}(2a+1)}\∫{\left({\mathrm{\η}}^{\′}\right)}^{2a}({\mathrm{\κ}}^2+{2\mathrm{\τ}}^2+{\mathrm{\τ}}^2{\mathrm{\η}}^2)\mathrm{du}$

其中a为泰勒展开的阶数，w为条带宽度，D为弯曲刚度，κ，τ，η分别为条带中心线的曲率，饶率以及饶率曲率比，′表示关于弧长u的导数。

步骤(3)中所述的由中心线描述的动能形为：

$T=\frac{1}{2}\mathrm{\ρ}({w\stackrel{\·}{r}}^2+\frac{w^3}{12}{\stackrel{\·}{\mathrm{\ω}}}^2+\frac{w^3}{6}{\mathrm{\η}}^{\′}\stackrel{\·}{r}\stackrel{\·}{\mathrm{\ω}})\mathrm{du}$

其中ρ为条带质量密度，w为条带宽度，

为条带中心线的线速度，

为条带直纹的变化速度。

步骤(4)中所述的“四元数-中心线节点坐标-饶率曲率比”这一广义坐标为g＝{q，r，η}；特点为条带中心线曲率和饶率可用四元数表示，即κ(q)，τ(q)；额外的广义坐标饶率曲率比η解决了曲率为零的区域计算饶率曲率比的数值奇异；引入的“相同饶率”一致性约束，c_sτ＝κ(q)η-τ(q)＝0则在后期处理时满足；显式四元数表示中心线标架和中心线节点坐标的使用易于用户交互，即中心线位置约束以及条带扭转约束的施加。

步骤(5)中所述的借助有限元并利用分段线性基函数在空间上对描述中心线的广义坐标进行离散化即可得到离散化的能量以及约束，得到能够计算的能量模型可逼近原解析表达式，而基函数的选取为更高阶。

步骤(6)中所述的分步求解运动方程：首先准静态地更新扭转部分的形变，然后动态更新弯曲部分的形变。具体特征为，动能包含中心线平动动能Tr和条带直纹的变化动能T_ω两部分：

$T_r=\frac{1}{2}\mathrm{\ρ}{\∫w\stackrel{\·}{r}}^2\mathrm{du},$ $T_{\mathrm{\ω}}=\frac{1}{2}\mathrm{\ρ}\∫(\frac{w^3}{12}{\stackrel{\·}{\mathrm{\ω}}}^2+\frac{w^3}{6}{\mathrm{\η}}^{\′}\stackrel{\·}{r}\stackrel{\·}{\mathrm{\ω}})\mathrm{du}$

当两者大小相当时T_ω对应形变在刚度上远大于T_r对应形变，即频高；而在刚度相当时，T_ω在大小上远小于T_r，即可忽略。所以T_ω对应的扭转部分形变可以被准静态地处理，而T_r对应的弯曲部分形变被动态处理。得运动方程求解流程：

$(r_{t_0},{\stackrel{\·}{r}}_{t_0},q_{t_0},{\mathrm{\η}}_{t_0})\stackrel{\mathrm{opt}}{\→}(r_{t_0},{\stackrel{\·}{r}}_{t_0},q_{t_1},{\mathrm{\η}}_{t_1})\stackrel{\mathrm{ode}}{\→}(r_{t_1},{\stackrel{\·}{r}}_{t_1},q_{t_1},{\mathrm{\η}}_{t_1})$

其中t₀，t₁为两个离散时刻，opt指通过优化准静态地处理扭转部分的形变而ode指求解微分方程动态处理弯曲部分的形变。

步骤(7)中所述的利用运动方程分步求解和时间域平行传输这两点来简化控制条带扭转的朝向约束为标架约束是指：在分步求解和中心线节点坐标的显式表示的前提下，对中心线法向d₁和副法向d₂做沿切向d₃时间积分所得曲线平行传输，朝向约束

即被转化为简单的标架约束

其中M＝(d₁，d₂，d₃)为中心线的标架，使用四元数q来表示。

步骤(8)中所述的过渡函数在条带接合点两侧定义了朝向不同的额外的能量、约束项，统一了非环状条带与可定向和不可定向环状条带的表达，且接合处在数值上与其他区域无差异。

与背景技术相比本发明具有的有益效果是：

本发明创造性的对rectifying developable进行拓展作为不可拉伸条带的几何描述并与Kirchhoff-Love薄壳模型相结合，进而利用Cosserat或Kirchhoff-Love线缆模型的广义坐标表示方法，提出了一套精确，高效且易于用户交互的不可拉伸条带模拟方法。其次，首先准静态地更新扭转部分的形变，然后动态更新弯曲部分的形变分步求解运动方程的策略将节约了将更多的计算资源分配给低频形变或对应动能较大的形变；这一处理又恰为控制扭转的约束通过时间域平行传输转化为更为简单的标架约束提供了可能。最后，结合提出的过渡函数，使得方法同时适用于非环状条带与可定向和不可定向环状条带的模拟。从结果可以看到，本发明方法能够精确展现条带各种有趣状态，而这是先前方法难以做到的。如紧绷下扭转条带所形成的三角形小块，放松后在拉紧条带开始卷起，扭转π所形成的条带中心线上只有一个零曲率点的Mobius环，以及扭转2π所形成的条带中心线上只有两个零曲率点的零字型环。

附图说明

图1是本发明的流程图，包含了关键部分和操作流程。

图2是非环状条带以及π和2π环状条带稳定状态与实验结果比较，模拟得到结果与真实结果几乎一致。

图3是2π环状带稳定状态的三视图，左为模拟结果右为实验结果，模拟得到结果与真实结果几乎一致。

图4是π和2π环状条带稳定状态中心线曲率、饶率分布图，得到的结果准确反映了条带中心线上的孤立零曲率点。

图5是图2左上与左下模拟过程条带总面积误差，模拟过程保持总面积不变，做到了不可拉伸。

图6是释放卷起条带后的运动轨迹。

图7是释放扭转2π并处于紧绷状态条带后的运动轨迹。

具体实施方式

下面结合附图和实施例对本发明作进一步说明。

如图1所示，本发明提出的一种几何精确的不可拉伸条带模拟方法包含8个关键部分：条带几何描述、弹性势能的导出、动能的导出、广义坐标的选取得到能量与约束、离散化能量与约束、分步求解运动方程策略的导出、朝向约束的简化以及过渡函数的导出。有了这8个关键部分之后，计算机的计算流程就是对运动方程的求解，具体分为6步，输入t₀时刻条带状态，准静态地更新扭转相关形变，后期处理A，动态地更新弯曲相关形变，后期处理B以及生成t1时刻条带形状。

本发明一种几何精确的不可拉伸条带模拟方法，包括以下关键部分：

(1)条带几何描述：

对于长度l宽度w的长方形条带，其曲面S可以被参数化表示为：

S(u，v)＝r(u)+vω(u)，ω(u)＝d₂(u)+η(u)d₃(u)，

η(u)＝τ(u)/κ(u)，u∈[0，l]，v∈[-w/2，w/2]

其中κ和τ分别为条带中心线r的曲率和饶率，η为饶率曲率比，而ω为直纹。三个单位向量k＝1，2，3组成了中心线上的材料标价

其中d₃＝r′/||r′||为r的单位切向量，d₁＝r″/||r″||为单位法向量，d₂＝d₃×d₁为单位副法向量。当r为正则曲线时S为传统的rectifying developable，而M为Frenet标架F。而非正则时，在曲率为零区域d₁无定义，在该区域前后d₁方向以及κ符号可能相反，导致直纹无定义或反转。因此统一使用S在中心线上一侧的法向作为d₁，并在曲率为零区域u∈[u_a，u_b]使用由以下极小化定义的η，

$\min {\∫}_{u_a}^{u_b}{\mathrm{\η}}^2\left(u\right)\mathrm{du},\mathrm{st}.c_{{\mathrm{\η}}^{\′}}\left(u\right),u\∈[u_a,u_b]$

为S所表示的曲可展，以下直纹不相交于条带曲面内约束必须被满足，

${\mathrm{\&eta;}}^{\&prime;}<\mathrm{sgn}\left({\mathrm{\&eta;}}^{\&prime;}\right)\frac{2}{w}$

由此可知，直纹面完全由其中心线r确定。图4是π和2π环状条带稳定状态中心线曲率、饶率分布图，使用本发明的几何描述能准确反映了条带中心线上的孤立零曲率点。

(2)弹性势能的导出：

自然平直且经受保尺度形变的条带其弹性势能可利用薄壳模型

来描述，其中H为平均曲率。利用条带可展的性质，即曲面S的主曲率之一为零，而另一主曲率为

因此可得，

$V=\frac{1}{2}D\&Integral;\frac{{\mathrm{\&kappa;}}^2{(1+{\mathrm{\&eta;}}^2)}^2}{{\mathrm{\&eta;}}^{\&prime;}}\ln \left(\frac{1+\frac{w}{2}{\mathrm{\&eta;}}^{\&prime;}}{1-\frac{w}{2}{\mathrm{\&eta;}}^{\&prime;}}\right)\mathrm{du}$

当

时，即在直纹不相交于条带曲面内约束违反时，能量函数奇异或无定义。因此首先忽略该约束，对以上能量函数中指数项做Taylor展开，

$V=\underset{a=0}{\overset{\&infin;}{\mathrm{\&Sigma;}}}\frac{w^{2a+1}D}{2^{2a+1}(2a+1)}{\&Integral;\left({\mathrm{\&eta;}}^{\&prime;}\right)}^{2a}{\mathrm{\&kappa;}}^2{(1+{\mathrm{\&eta;}}^2)}^2\mathrm{ds}$

以此作为弹性势能，并使用后期处理的策略来使直纹不相交于条带曲面内约束满足。尽可能消除可能导致不定项的η，得最终的弹性势能，

$V=\underset{a=0}{\overset{\&infin;}{\mathrm{\&Sigma;}}}\frac{w^{2a+1}D}{2^{2a+1}(2a+1)}\&Integral;{\left({\mathrm{\&eta;}}^{\&prime;}\right)}^{2a}({\mathrm{\&kappa;}}^2+{2\mathrm{\&tau;}}^2+{\mathrm{\&tau;}}^2{\mathrm{\&eta;}}^2)\mathrm{ds}$

图6是释放卷起条带后的运动轨迹，反映了弹性势能逐渐耗散的过程

(3)动能的导出：

通过累加S上所有质点的线速度之和，可得线动能，

$T=\frac{1}{2}\mathrm{\&rho;}\&Integral;({w\stackrel{\&CenterDot;}{r}}^2+\frac{w^3}{12}{\stackrel{\&CenterDot;}{\mathrm{\&omega;}}}^2+\frac{w^3}{6}{\mathrm{\&eta;}}^{\&prime;}\stackrel{\&CenterDot;}{r}\stackrel{\&CenterDot;}{\mathrm{\&omega;}})\mathrm{du}$

动能包含中心线平动动能T_r和条带直纹的变化动能T_ω两部分：

$T_r=\frac{1}{2}\mathrm{\&rho;}\&Integral;{w\stackrel{\&CenterDot;}{r}}^2\mathrm{du},$ $T_{\mathrm{\&omega;}}=\frac{1}{2}\mathrm{\&rho;}\&Integral;(\frac{w^3}{12}{\stackrel{\&CenterDot;}{\mathrm{\&omega;}}}^2+\frac{w^3}{6}{\mathrm{\&eta;}}^{\&prime;}\stackrel{\&CenterDot;}{r}\stackrel{\&CenterDot;}{\mathrm{\&omega;}})\mathrm{du}$

当两者大小相当时T_ω对应形变在刚度上远大于T_r对应形变，即频高；而在刚度相当时，T_ω在大小上远小于T_r，即可忽略。所以T_ω对应的扭转部分形变可以被准静态地处理，而T_r对应的弯曲部分形变被动态处理。图7是释放扭转2π并处于紧绷状态条带后的运动轨迹，简化的动能在降低计算开销的同时依然保证的扭转形变所致运动的准确性。

(4)广义坐标的选取得到能量与约束：

选用广义坐标“四元数-中心线节点坐标-饶率曲率比”，g＝{r，q，η}来表示条带中心线。显式用η表示饶率曲率比解决了平直区域η＝τ/κ带来的数值奇异，由此引入一致性约束：相同饶率约束，c_sτ＝κ(q)η-τ(q)＝0。其他由此广义坐标带来的内部约束包括：单位四元数约束，c_u＝||q||＝1；中心线不可拉伸约束，c_l＝||r′||＝1；切向平行约束，

副法相曲率为零约束c_γ＝γ＝-d₂·d′₃＝0；以及之前提及的直纹不相交于条带曲面内约束，

用户操作相关的约束则包括，位置约束c_p＝r(u_k)-p_k＝0；朝向约束见图2上部，是非环状条带扭转后的结果，模拟得到结果能在较少广义坐标的离散节点下模拟出于真实结果几乎一致的结果。图3是2π环状带稳定状态的三视图，左为模拟结果右为实验结果，模拟得到结果与真实结果几乎一致。

(5)离散化能量与约束：

借助有限元并利用分段线性基函数在空间上对描述中心线的广义坐标进行离散化话可得n个中心线节点r₀，...，r_n-1(形成n-1条边e⁰，...，e^n-2，其中eⁱ＝r_i+1-r_i)，中间n-2个节点上的n-2个标量η₀，...，η_n-3，以及每条边上的共n-1个四元数q⁰，...，q^n-2，由此可得基本的离散量，

${\left(r^{\&prime;}\right)}^i=\frac{r_{i+1}-r_i}{{\stackrel{\&OverBar;}{l}}^i},$ ${\left(q^{\&prime;}\right)}_i=\frac{2(q^{i+1}-q^i)}{{\stackrel{\&OverBar;}{l}}^{i+1}+{\stackrel{\&OverBar;}{l}}^i},$ ${\left({\mathrm{\&eta;}}^{\&prime;}\right)}^i=\frac{{\mathrm{\&eta;}}_{i+1}-{\mathrm{\&eta;}}_i}{{\stackrel{\&OverBar;}{l}}^{i+1}}$

${\stackrel{\widetilde{.}}{r}}^i=\frac{{\stackrel{.}{r}}_i+{\stackrel{.}{r}}_{i+1}}{2},$ ${\tilde{q}}_i=\frac{{\stackrel{\&OverBar;}{l}}^i{q}^i+{\stackrel{\&OverBar;}{l}}^{i+1}{q}^{i+1}}{{\stackrel{\&OverBar;}{l}}^{i+1}+{\stackrel{\&OverBar;}{l}}^i},$ ${\left({\widetilde{\mathrm{\&eta;}}}^{\&prime;}\right)}_i=\frac{{\mathrm{\&eta;}}_{i+1}-{\mathrm{\&eta;}}_{i-1}}{({\stackrel{\&OverBar;}{l}}^{i+1}+{\stackrel{\&OverBar;}{l}}^i)}$

其中带‘～’的量表示由基函数构造得到，带‘-’的量表示无弹性势能时的测得的量。离散副法向曲率，曲率以及饶率，

${\mathrm{\&gamma;}}_i\left(q\right)={2B}_1{\tilde{q}}_i\&CenterDot;{\left(q^{\&prime;}\right)}_i,$ ${\mathrm{\&kappa;}}_i\left(q\right)={2B}_2{\tilde{q}}_i\&CenterDot;{\left(q^{\&prime;}\right)}_i,$ ${\mathrm{\&tau;}}_i\left(q\right)={2B}_3{\tilde{q}}_i\&CenterDot;{\left(q^{\&prime;}\right)}_i$

其中

k＝1，2，3为反对称阵，

$B_{k}q=\left\{\begin{array}{cc}{(q_w,q_z,-q_y,-q_x)}^T& k=1\\ {({-q}_z,q_w,q_x,-q_y)}^T& k=2\\ {(q_y,-q_x,q_w,-q_z)}^T& k=3\end{array}\right.$

离散弹性势能V＝V_κ+V_τ+V_η，

$V_{\mathrm{\&kappa;}}=\underset{a=0}{\overset{\&infin;}{\mathrm{\&Sigma;}}}\frac{w^{2a+1}D}{2^{2a+2}(2a+1)}\underset{i=0}{\overset{n-3}{\mathrm{\&Sigma;}}}({\stackrel{\&OverBar;}{l}}^i+{\stackrel{\&OverBar;}{l}}^{i+1}){\left({\widetilde{\mathrm{\&eta;}}}^{\&prime;}\right)}_i^{2a}{\mathrm{\&kappa;}}_i^2\left(q\right)$

$V_{\mathrm{\&tau;}}=\underset{a=0}{\overset{\&infin;}{\mathrm{\&Sigma;}}}\frac{2w^{2a+1}D}{2^{2a+2}(2a+1)}\underset{i=0}{\overset{n-3}{\mathrm{\&Sigma;}}}({\stackrel{\&OverBar;}{l}}^i+{\stackrel{\&OverBar;}{l}}^{i+1}){\left({\widetilde{\mathrm{\&eta;}}}^{\&prime;}\right)}_i^{2a}{\mathrm{\&tau;}}_i^2\left(q\right)$

$V_{\mathrm{\&eta;}}=\underset{a=0}{\overset{\&infin;}{\mathrm{\&Sigma;}}}\frac{w^{2a+1}D}{2^{2a+2}(2a+1)}\underset{i=0}{\overset{n-3}{\mathrm{\&Sigma;}}}({\stackrel{\&OverBar;}{l}}^i+{\stackrel{\&OverBar;}{l}}^{i+1}){\left({\widetilde{\mathrm{\&eta;}}}^{\&prime;}\right)}_i^{2a}{\mathrm{\&tau;}}_i^2\left(q\right){\mathrm{\&eta;}}_i^2$

离散动能T，由于直纹的变化动能T_ω被准静态地处理，只需对平动动能T_r离散，

$T=\frac{1}{8}\mathrm{\&rho;w}\underset{i=0}{\overset{n-2}{\mathrm{\&Sigma;}}}{\stackrel{\&OverBar;}{l}}^i({\stackrel{.}{r}}_i+{\stackrel{.}{r}}_{i+1})\&CenterDot;({\stackrel{.}{r}}_i+{\stackrel{.}{r}}_{i+1}),$

离散的内部约束，

c_u，i＝||qⁱ||＝1，i∈{0，..，n-2}

$c_{l,i}=l^i\&CenterDot;l^i-{\stackrel{\&OverBar;}{l}}^i\&CenterDot;{\stackrel{\&OverBar;}{l}}^i=0,i\&Element;\{0,...,n-2\}$

$c_{\mathrm{para},i}=\frac{r_{i+1}-r_i}{{\stackrel{\&OverBar;}{l}}^i}-d_3^i=0,i\&Element;\{0,...,n-2\}$

$c_{\mathrm{\&gamma;},i}=-2B_1{\tilde{q}}_i\&CenterDot;{\left(q^{\&prime;}\right)}_i=0,i\&Element;\{0,...,n-3\}$

c_sτ，i＝κ_i(q)η_i-τ_i(q)＝0，i∈{0，..，n-3}

$c_{{\mathrm{\&eta;}}^{\&prime;},i}:\frac{{\mathrm{\&eta;}}_{i+1}-{\mathrm{\&eta;}}_i}{{\stackrel{\&OverBar;}{l}}_{i+1}}\&le;\mathrm{sgn}\left(\frac{{\mathrm{\&eta;}}_{i+1}-{\mathrm{\&eta;}}_i}{{\stackrel{\&OverBar;}{l}}_{i+1}}\right)\frac{2}{w}$

其中对于不等约束c_η′，i，使用积极集法将其转化为等式约束，具体为，在不带约束的情况下进行求解，判断是否违反不等约束c_η′，i，存在违反则收集为等式约束进行后处理使其满足，收集策略为，

$c_{{\mathrm{\&eta;}}^{\&prime;},\mathrm{eq},i}=\left\{\begin{array}{cc}\frac{{\mathrm{\&eta;}}_{i+1}-{\mathrm{\&eta;}}_i}{{\stackrel{\&OverBar;}{l}}_{i+1}}=-\frac{2}{w}& \frac{{\mathrm{\&eta;}}_{i+1}-{\mathrm{\&eta;}}_i}{{\stackrel{\&OverBar;}{l}}_{i+1}}<-\frac{2}{w}\\ \frac{{\mathrm{\&eta;}}_{i+1}-{\mathrm{\&eta;}}_i}{{\stackrel{\&OverBar;}{l}}_{i+1}}=\frac{2}{w}& \frac{{\mathrm{\&eta;}}_{i+1}-{\mathrm{\&eta;}}_i}{{\stackrel{\&OverBar;}{l}}_{i+1}}>\frac{2}{w}\\ \mathrm{none}& \mathrm{others}\end{array}\right.$

对于约束的处理，使用快速投影法(GOLDENTHAL，R.，HARMON，D.，FATTAL，R.，BERCOVIER，M.，AND GRINSPUN，E.2007.Efficient simulation ofinextensible cloth.ACM Transactions on Graphics(Proceedings of SIGGRAPH2007)26，3)这一后期处理策略。图5是图2左上与左下模拟过程条带总面积误差，通过快速投影法处理约束使得总面积保持不变，做到了不可拉伸。

分步求解运动方程策略的导出

运动方程求解流程为：

$(r_{t_0},{\stackrel{\&CenterDot;}{r}}_{t_0},q_{t_0},{\mathrm{\&eta;}}_{t_0})\stackrel{\mathrm{opt}}{\&RightArrow;}(r_{t_0},{\stackrel{\&CenterDot;}{r}}_{t_0},q_{t_1},{\mathrm{\&eta;}}_{t_1})\stackrel{\mathrm{ode}}{\&RightArrow;}(r_{t_1},{\stackrel{\&CenterDot;}{r}}_{t_1},q_{t_1},{\mathrm{\&eta;}}_{t_1})$

其中t₀，t₁为两个离散时刻，opt指通过优化准静态地处理扭转部分的形变，即优化函数，而ode指求解微分方程动态处理弯曲部分的形变，即求解方程，

(6)朝向约束的简化：

在分步求解和中心线节点坐标的显式表示的前提下，对于t₀，t₁两个离散时刻，t₁时刻的中心线法向d₁和副法向d₂可通过做沿切向d₃时间积分所得曲线平行传输得到，

$d_{j,t_1}^k=P_{\mathrm{dis}}(d_{3,t_0}^k,d_{3,t_1}^k)d_{j,t_0}^k,j\&Element;\left\{\mathrm{1,2}\right\}$

其中P_dis为平行传输操作，它可以是构成的最小旋转，结合由

可知的

，t₁时刻的朝向约束所对应的标架可以在ode完成后显式更新。因此，朝向约束被简化为标架约束，

(7)过渡函数的导出：

对于环状的条带，则需要接合点处额外的能量和约束。在使用的分段线性基函数下，具体包括，离散的环状约束，c_loop：r₀-r_n-1＝0；接合点处额外的弹性势能Vc；额外的副法相曲率为零约束c_γ，c；额外的直纹不相交于条带曲面内约束c_η′，c。在不可定向条带(顺时针旋转

n＝1，2，...接合所得)的例子中需要对计算这些额外量的基本量做方向(符号)的统一，如曲率的计算，对于0⁺侧的梯度则使用0⁺的d₁方向和0^-侧-d₁方向的量计算得到的曲率，而对于0^-侧的梯度则以0^-处的d₁和0⁺的-d₁方向的量计算得到的曲率。因此对于0^-侧的标架q^n-2，转化为与0⁺的q⁰统一标架

其中‘*’为四元数乘法，p为另一四元数，对应于绕轴e^n-2做

的旋转；而其他量η，κ(q)，τ(q)通过取负即可完成方向转换。图2下部和图3为使用过渡函数实现的环状结果，模拟得到的结果与实验结果几乎一致。

使用以上的关键部分，计算机的计算流程就是对运动方程的求解，具体分为5步：

(1)输入t₀时刻状态条带状态

(2)使用最速下降法准静态地更新扭转相关形变，得

(3)后期处理A

a.在曲率为零区域u∈[u_a，u_b]使用使用最速下降法求出由极小化对应函数所定义的η

b.使用快速投影法约束标架约束，并更新

c.使用快速投影法约束相同饶率约束和副法相曲率为零约束，并更新

d.收集并转化直纹不相交于条带曲面内约束为等式约束

e.使用快速投影法约束直纹不相交于条带曲面内约束，并更新

(4)使用半隐式Euler发动态地更新弯曲相关形变，得

(5)后期处理B

a.使用快速投影法约束更新中心线不可拉伸约束，并更新

b.使用快速投影法约束更新切向平行约束，并更新

c.使用平行传输更新朝向约束所对应的标架

(6)生成t₁时刻条带形状，

Claims (9)

1.一种几何精确的不可拉伸条带模拟方法，其特征在于包括以下八个步骤：

(1)使用特殊的直纹面解析地描述不可拉伸条带对应的二维曲面，即不可拉伸条带形状由其中心线描述；

(2)根据步骤(1)建立不可拉伸条带上各点处的平均曲率与其中心线上各点的曲率、饶率以及饶率曲率比等的关系，结合薄壳模型，导出由中心线描述的弹性势能；

(3)根据步骤(1)建立不可拉伸条带上各点处的速度与其中心线上各点的速度、直纹变化速度的关系，导出由中心线描述的动能；

(4)“四元数-中心线节点坐标-饶率曲率比”这一广义坐标来描述根据步骤(1)中条带中心线，由该广义坐标导出的步骤(2)中描述的弹性势能、步骤(3)中描述的动能以及保证直纹面可展、广义坐标冗余自由度一致的约束数值无奇异；且易于用户操作相关约束的施加；

(5)借助有限元并利用分段线性基函数在空间上对步骤(4)中描述中心线的广义坐标进行离散化，导出离散可计算的弹性势能，动能以及约束；

(6)基于条带扭转部分形变的频率远大于弯曲部分的频率，分步求解由步骤(5)中离散能量构成的运动方程：首先准静态地更新扭转部分的形变，然后动态更新弯曲部分的形变；

(7)在步骤(6)求解过程中，利用运动方程分步求解这一特点同时借助时间域平行传输，将控制条带扭转的朝向约束转化为求解更为简单标架约束；

(8)在步骤(6)求解过程中，通过定义合适过渡函数，本条带模拟方法能够拓展到可定向或不可定向环状条带的模拟。

2.根据权利要求1所述的一种几何精确的不可拉伸条带模拟方法，其特征在于：步骤(1)中所述的特殊的直纹面是一种保证曲面可展的直纹面，是对rectifying developable的拓展，提供了条带形变过程几何上精确的前提；rectifyingdevelopable由三维空间中的正则曲线描述，而平直或者波浪状的条带存在曲率为零的区域，即其中心线可能为非正则，因此rectifying developable被拓展以描述平直或者波浪状的条带，具体指对直纹的重新定义，使得在曲率为零的区域描述直纹所需的中心线法向量和饶率曲率比有定义；对于中心线法向量，始终使用条带一侧的在中心线处的表面法向；此外要求曲率为零区域的饶率曲率比平方和最小。

3.根据权利要求1所述的一种几何精确的不可拉伸条带模拟方法，其特征在于：步骤(2)中所述的由薄壳模型推导所得的弹性势能中包含了一个指数项，当直纹相交于条带曲面内时该指数项将弹性势能惩罚至无穷，即弹性势能蕴含了“直纹不相交于条带曲面内”这一硬约束；因而在直纹相交状态下数值上奇异或无定义；为使得数值能够计算，指数项被泰勒展开，而直纹相交于条带曲面内这一原本蕴含于弹性势能的硬约束通过后期处理来满足；弹性势能形为：

$V=\underset{a=0}{\overset{\&infin;}{\mathrm{\&Sigma;}}}\frac{w^{2a+1}D}{2^{2a+1}(2a+1)}\&Integral;{\left({\mathrm{\&eta;}}^{\&prime;}\right)}^{2a}({\mathrm{\&kappa;}}^2+{2\mathrm{\&tau;}}^2+{\mathrm{\&tau;}}^2{\mathrm{\&eta;}}^2)\mathrm{du}$

其中a为泰勒展开的阶数，w为条带宽度，D为弯曲刚度，κ，τ，η分别为条带中心线的曲率，饶率以及饶率曲率比，′表示关于弧长u的导数。

4.根据权利要求1所述的一种几何精确的不可拉伸条带模拟方法，其特征在于：步骤(3)中所述的由中心线描述的动能形为：

$T=\frac{1}{2}\mathrm{\&rho;}({w\stackrel{\&CenterDot;}{r}}^2+\frac{w^3}{12}{\stackrel{\&CenterDot;}{\mathrm{\&omega;}}}^2+\frac{w^3}{6}{\mathrm{\&eta;}}^{\&prime;}\stackrel{\&CenterDot;}{r}\stackrel{\&CenterDot;}{\mathrm{\&omega;}})\mathrm{du}$

其中ρ为条带质量密度，w为条带宽度，

为条带中心线的线速度，

为条带直纹的变化速度。

5.根据权利要求1所述的一种几何精确的不可拉伸条带模拟方法，其特征在于：步骤(4)中所述的“四元数-中心线节点坐标-饶率曲率比”这一广义坐标为g＝{q，r，η}；特点为条带中心线曲率和饶率可用四元数表示，即κ(q)，τ(q)；额外的广义坐标饶率曲率比η解决了曲率为零的区域计算饶率曲率比的数值奇异；引入的“相同饶率”一致性约束，c_sτ＝κ(q)η-τ(q)＝0则在后期处理时满足；显式四元数表示中心线标架和中心线节点坐标的使用易于用户交互，即中心线位置约束以及条带扭转约束的施加。

6.根据权利要求1所述的一种几何精确的不可拉伸条带模拟方法，其特征在于：步骤(5)中所述的借助有限元并利用分段线性基函数在空间上对描述中心线的广义坐标进行离散化即可得到离散化的能量以及约束，得到能够计算的能量模型可逼近原解析表达式，而基函数的选取为更高阶。

7.根据权利要求1所述的一种几何精确的不可拉伸条带模拟方法，其特征在于：步骤(6)中所述的分步求解运动方程：首先准静态地更新扭转部分的形变，然后动态更新弯曲部分的形变；具体特征为，动能包含中心线平动动能T_r和条带直纹的变化动能T_ω两部分：

$T_r=\frac{1}{2}\mathrm{\&rho;}{\&Integral;w\stackrel{\&CenterDot;}{r}}^2\mathrm{du},$ $T_{\mathrm{\&omega;}}=\frac{1}{2}\mathrm{\&rho;}\&Integral;(\frac{w^3}{12}{\stackrel{\&CenterDot;}{\mathrm{\&omega;}}}^2+\frac{w^3}{6}{\mathrm{\&eta;}}^{\&prime;}\stackrel{\&CenterDot;}{r}\stackrel{\&CenterDot;}{\mathrm{\&omega;}})\mathrm{du}$

当两者大小相当时T_ω对应形变在刚度上远大于T_r对应形变，即频高；而在刚度相当时，T_ω在大小上远小于T_r，即可忽略；所以T_ω对应的扭转部分形变可以被准静态地处理，而T_r对应的弯曲部分形变被动态处理，得运动方程求解流程：

$(r_{t_0},{\stackrel{\&CenterDot;}{r}}_{t_0},q_{t_0},{\mathrm{\&eta;}}_{t_0})\stackrel{\mathrm{opt}}{\&RightArrow;}(r_{t_0},{\stackrel{\&CenterDot;}{r}}_{t_0},q_{t_1},{\mathrm{\&eta;}}_{t_1})\stackrel{\mathrm{ode}}{\&RightArrow;}(r_{t_1},{\stackrel{\&CenterDot;}{r}}_{t_1},q_{t_1},{\mathrm{\&eta;}}_{t_1})$

其中t₀，t₁为两个离散时刻，opt指通过优化准静态地处理扭转部分的形变而ode指求解微分方程动态处理弯曲部分的形变。

8.根据权利要求1所述的一种几何精确的不可拉伸条带模拟方法，其特征在于：步骤(7)中所述的利用运动方程分步求解和时间域平行传输这两点来简化控制条带扭转的朝向约束为标架约束是指：在分步求解和中心线节点坐标的显式表示的前提下，对中心线法向d₁和副法向d₂做沿切向d₃时间积分所得曲线平行传输，朝向约束

即被转化为简单的标架约束

其中M＝(d₁，d₂，d₃)为中心线的标架，使用四元数q来表示。

9.根据权利要求1所述的一种几何精确的不可拉伸条带模拟方法，其特征在于：步骤(8)中所述的过渡函数在条带接合点两侧定义了朝向不同的额外的能量、约束项，统一了非环状条带与可定向和不可定向环状条带的表达，且接合处在数值上与其他区域无差异。

Images