不可伸展带状物体的几何性质描述模型及动力学模拟方法

摘要

本发明公开了一种不可伸展带状物体的几何描述模型及动力学模拟方法。对于不可伸展带状物体，采用正则中心曲线作为所述几何性质的描述模型，模型包含有中心曲线保长约束和母线不相交约束与一个边界条件。构造广义坐标来表示带状曲面，对矩形带状可展曲面建立弹性势能模型，该弹性势能模型用于计算该矩形带状可展曲面中的内应力；通过数值模拟计算方法计算矩形带状可展曲面的运动。本发明针对带状物体能够精确地描述其几何性质，能够精确地保证带状曲面的不可伸展性，在模拟时所需的自由度较少，不需要大量的保长约束，能够处理模拟中可能出现的退化情况，从而使得计算的效率与稳定性较之传统方法都有较大的提升。

说明书

技术领域

本发明涉及了一种模型和物理模拟方法，尤其涉及了一种不可伸展带状物 体的几何性质描述模型及动力学模拟方法，主要针对于矩形带状可展曲面。

技术背景

带状物体在日常生活中相当常见，比如胶卷、纸条、装饰彩带等等。该类 物体一般较窄，自然状态下平坦同时具有不可拉伸的性质。在几何上，带状物 体可被抽象为一张光滑的可展曲面，即直纹曲面。这种特殊的几何性质赋予带 状物体一些重要的静力学和动力学特性，使得现有方法在模拟这类物体时存在 很大的挑战。

虽然物理仿真领域已有不少针对杆状物和薄壳的模拟方法，但综合考虑几 何精确性、物理精确性和效率等方面时，这些方法并不适合用于模拟带状物体。 一种传统的方法就是在薄壳模拟的基础上施加硬约束来保证曲面的不可伸展 性。然而这种模式由于需要高分辨率的网格与大量的保长约束，会直接导致计 算性能的大幅下降。

发明内容

针对背景技术的不足，本发明的目的在于提出了一种不可伸展带状物体的 几何性质描述模型及动力学模拟方法，能够精确描述带状曲面，并基于模型实 现高效的物理模拟，在满足几何、物理精确性的同时能够兼顾计算效率。

为实现上述目的，本发明的采用的技术方案具有如下内容。

一、一种不可伸展带状物体的几何性质描述模型：

对于一张长为l，宽为w的矩形带状可展曲面，采用由以下正则中心曲线r 来描述的公式作为所述几何性质的描述模型，该模型在连续域上可由下述公式 表达：

S(u,v)＝r(u)+vω(u),

ω(u)＝d₂(u)+η(u)d₃(u),

η(u)＝τ(u)/κ(u)

u∈[0,l],v∈[-w/2,w/2]

式中，S(u,v)表示三维空间中的一张矩形带状可展曲面，u、v是曲面参数平 面上代表两个相互垂直方向的参数坐标，ω(u)代表母线；d₃(u),d₂(u)分别代表沿 着正则中心曲线r(u)的材质标架M中的单位切向量和单位副法向量，M表示沿 着正则中心曲线r(u)的材质标架，由d₁(u),d₂(u),d₃(u)三个分量构成，可表示为 M＝(d₁(u),d₂(u),d₃(u))∈R^3×3，其中R表示实数集，d₁(u)代表沿着r(u)的材质标架中 的单位法向量，由于r是S上的一条测地线，所以该标架M总是与Frenet标架 一致；v代表曲线上一点在母线ω(u)方向上的投影，η(u)表示母线ω(u)在单位切 向量上的投影，数值上等于正则中心曲线r(u)上一点挠率与曲率的比值；κ(u)和 τ(u)分别表示正则中心曲线r(u)上任意一点处的曲率和挠率，计算方法为 κ(u)＝d₃′(u)·d₁(u)，τ(u)＝d₂′(u)·d₃(u)，d3’(u)代表沿着r(u)的材质标架的单位法向 量对正则中心曲线参数坐标u的一阶导数，d2’(u)代表沿着r(u)的材质标架的单 位副法向量对正则中心曲线参数坐标u的一阶导数；上述字母含义关系如图1 所示。

所述几何性质描述模型包含有中心曲线保长约束和母线不相交约束的两个 约束与一个边界条件。

所述的中心曲线保长约束c_l为||r′||＝1，以保证曲面的不可拉伸性，其中r’ 表示正则中心曲线r(u)对中心曲线参数坐标u的一阶导数。

所述的母线不相交约束c_η′为|η′|＜2/w，以使得由母线交点所形成的曲线位 于S之外，η′表示η(u)对中心曲线参数坐标u的一阶导数。[可参考BOP.,WANG W.:Geodesic-controlleddevelopablesurfacesformodelingpaperbending.CGF(EG 07)26,3(2007),365–374]。

所述的边界条件为矩形带状可展曲面S两端的边与单位副法向量d₂平行， 即η(0)＝η(l)＝0，这样保证带状曲面上的所有点都能够被母线覆盖。

二、一种不可伸展带状物体的动力学模拟方法：

1)为了便于数值模拟，保证模拟的性能与稳定性，为模拟计算所述几何性 质描述模型构造广义坐标g来表示带状曲面，广义坐标g由以下三元组表示：

g＝{q,η,r}

其中，r表示中心曲线的位置坐标，q是代表中心曲线r上材质标架M的一 个四元数，η代表母线在单位切向量上的投影。

为了使上述广义坐标g所描述的带状物体具有不可伸展的几何性质，需要 包含如下约束：

A)单位四元数约束c_u:||q||＝1，q是表示中心曲线上材质标架的四元数；

B)平行约束r’是中心曲线r对参数坐标u的一阶导数， d₃(q)是材质标架中的单位切向量。

C)副法线零曲率约束c_γ:d₂(q)′·d₃(q)＝0，d₂(q)′是材质标架中单位副法向量对 中心曲线参数坐标u的一阶导数，d₃(q)是材质标架中的单位切向量。

D)保长约束c_l:||r′||-1＝0，r’是中心曲线位置坐标r对中心曲线参数坐标u的 一阶导数。

E)母线不相交约束c_η′：|η′|＜2/w，η′表示η对中心曲线参数坐标u的一阶 导数，w表示矩形带状可展曲面S的宽度；

F)等挠率约束c_sτ:κ(q)η-τ(q)＝0，τ(q),κ(q)分别表示挠率和曲率。

2)对矩形带状可展曲面S建立弹性势能模型，通过该弹性势能模型计算该 矩形带状可展曲面S中的内应力；

3)通过数值模拟计算方法计算矩形带状可展曲面S的运动。

广义坐标g中还可添加由用户施加的外部约束：

G)位置约束c_p:r(u_k)-p_k＝0，r(u_k)代表中心曲线的某一点，p_k代表空间中 固定的一点。

H)方向约束代表材质标架中单位法向量对时间的一 阶导数，d₂(q)代表材质标架中的单位副法向量。

所述步骤2)中的弹性势能V具体采用以下公式：

$V=\underset{a=0}{\overset{\infty }{\Sigma }}\frac{w^{2a+1}D}{(2a+1)}\int {\left({\eta }^{\prime }\right)}^{2a}({\kappa }^2+2{\tau }^2+{\tau }^2{\eta }^2)du---\left(1\right)$

式中，a代表泰勒展开式的阶数，D代表材质刚度，u表示曲面参数平面上 代表其中一个方向的参数坐标，w表示矩形带状可展曲面S的宽度，η代表母线 在单位切向量上的投影，κ和τ分别表示正则中心曲线r上任意一点处的曲率和 挠率。具体实施中可将上述弹性势能模型的泰勒展开式阶数a从5阶截断。

所述的弹性势能模型可采用以下方式推导得到：带状曲面在等距形变下的 弹性势能可通过下面的公式描述：

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

其中，V代表势能，D代表材质刚度。[可参考GIOMIL.,MAHADEVANL.: Statisticalmechanicsofdevelopableribbons.Phys.Rev.Lett.104,23(2010), 238104]。为了移除公式中可能出现的退化情况，对上述公知的势能函数中的对 数项进行泰勒展开进而获得弹性势能。

考虑到模拟的效率，本发明所述步骤3)的数值模拟计算方法只考虑带状物 体运动中的低频部分，把扭转运动按准静态方式处理。给定t₀时刻曲面的状态， 在一个时间步长内，方法包含了两个阶段，微分方程求解阶段(ode)和准静态 优化阶段(opt)，流程如下所示：

其中，分别代表q,η,r和r对时间的一阶导数在t₀时刻的取值。 分别代表q,η,r和r对时间的一阶导数在t₁时刻的取值。

所述步骤3)的数值模拟计算方法具体如下：在每一个时间步长内，构建运 动方程动态地更新中心曲线，更新中心曲线后通过后处理方式保持中心曲线保 长约束c_l，之后迭代进行下述过程直至结果收敛：通过准静态优化更新母线，更 新母线后通过后处理方式保持母线不相交约束c_η′、副法线零曲率约束c_γ、等挠 率约束c_st和单位四元数约束c_u。

所述构建运动方程动态地更新中心曲线，更新中心曲线后通过后处理方式 保持中心曲线保长约束c_l具体为：

3.1)构建以下公式的运动方程动态地更新中心曲线：

$M\stackrel{··}{r}=-{▿}_r{E}_{para}+F_e$

式中，M是代表中心曲线质量分布的对角矩阵，代表函数对正则中心曲 线r的一阶导数，代表正则中心曲线r对时间的二阶导数，E_para是正则中心曲 线r的势能函数，F_e是用户施加的力或重力；

在求解运动方程阶段，可使用辛欧拉法来积分运动方程，从而得到更新后 中心曲线上的节点位置。

所述的正则中心曲线r的势能函数E_para采用以下公式计算：

E_para＝w_para||c_para||²

其中，w_para是势能函数E_para的权重，c_para是平行约束，形式为： $\frac{r^{\prime }}{\left|\right|r^{\prime }\left|\right|}-d_3\left(q\right)=0.$

3.2)更新中心曲线后通过以下方式的后处理保持中心曲线保长约束c_l：

使用快速投影方法来计算以下拉格朗日函数的驻点：

$L_{\stackrel{·}{r}}=\frac{1}{2}{\stackrel{·}{r}}^{T}M\stackrel{·}{r}+c_l{\lambda }_l$

其中，λ_l是拉格朗日乘子，M是代表中心曲线质量分布的对角矩阵，是中 心曲线位置坐标对时间的一阶导数，是的转置。

[快速投影方法可参考GOLDENTHALR.,HARMOND.,FATTALR., BERCOVIERM.,GRINSPUNE.:Efficientsimulationofinextensiblecloth.ACM TOG(SIGGRAPH07)26,3(2007),49]。

所述通过准静态优化更新母线，更新母线后通过后处理方式保持所述四个 约束具体为：

3.3)在准静态优化阶段，为了保持母线与更新后中心曲线状态的一致性， 采用以下公式表示的无约束优化问题对母线进行优化更新，最优解记为q^*,η^*：

$\underset{(q,\eta )}{\min }V+E_{para}$

式中，q是表示中心曲线标架M的四元数；

3.4)通过后处理的方式保持母线不相交约束c_η′、副法线领取率约束c_γ、等 挠率约束c_sτ和单位四元数约束c_u等四个约束。方法为求解拉格朗日函数L_q,η的驻 点并将四元数归一化。拉格朗日函数的具体形式为：

$L_{q,\eta }=\frac{1}{2}{\mathrm{\Delta q}}^{T}I\Delta q+\frac{1}{2}{\mathrm{\Delta \eta }}^{T}I\Delta \eta +{\lambda }_{{\eta }^{\prime }}{c}_{{\eta }^{\prime }}+{\lambda }_{\eta }{c}_{\eta }+{\lambda }_{s\tau }{c}_{s\tau }$

其中，I为单位矩阵，Δq,Δη分别是q,η相对于最优解q^*,η^*的增量，即 Δq＝q-q^*,Δη＝η-η^*。Δq^T,Δη^T分别是Δq,Δη的转置。λ_η′,λ_η,λ_sτ分别是对应于约束 c_η′,c_η,c_sτ的拉格朗日乘子。

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

1)带状物体的几何描述。本发明使用基于中心曲线的表示方法来描述带状 物体的几何模型，这种描述方式相比于薄壳模型具有更少的自由度，并且能够 精确地保证曲面的不可伸展性，从而省去了数值求解时处理保距约束所需的计 算代价。同时，这种表示方法能够处理动力学模拟中出现的数值退化和中心曲 线的非正则退化。

2)能量模型。用户的输入、碰撞和数值方法有时会使势能趋向于无穷，进 而导致问题求解时的不稳定。本发明通过有限阶的泰勒展开和一系列约束来近 似带状物体的弹性势能，在解决了上述问题的同时并不损失原势能所具有的一 些关键特性。对于系统动能，本发明观察到带状曲面非中心曲线部分对系统的 动力学效应影响甚微，基于上述观察对运动方程做了相应的简化。

3)积分方法。本发明提出了一种高效的两阶段时间积分模式：根据运动方 程动态地演进中心线，然后通过准静态优化更新母线。这种模式避免了一些数 值求解上的困难，提高了求解的效率。

由此，本发明针对带状物体能够保持其精确地描述其几何性质，它能够精 确地保证带状曲面的不可伸展性，在模拟时所需的自由度较少，不需要大量的 保长约束，能够处理模拟中可能出现的退化情况，从而使得计算的效率与稳定 性较之传统的方法都有较大的提升。

附图说明

图1是本发明带状曲面描述模型参数化的图示。

图2是本发明数值模拟方法部分的流程图。

图3～图6是本发明模拟出的纸带效果图。

图7～图10是本发明方法与其他方法模拟效果的一个对比。模拟使用的方法 依次是：图7是本发明的方法；图8是采用相容有限元的薄壳模拟方法；图9 是采用非相容有限元的薄壳模拟方法；图10使用弹性杆模拟方法。

图11是两种情形下迭代次数与纸带总面积的相对误差的之间关系的统计 图。

具体实施方式

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

本发明文中的符号‘'’均代表变量对正则中心曲线参数坐标u的一阶导数， 符号‘‥’均代表变量对时间的二阶导数。

本发明的实施例如下：

1)构造出广义坐标g，广义坐标g由以下三元组表示：

g＝{q,η,r}

2)对矩形带状可展曲面S建立采用以下公式的弹性势能模型，通过该弹性 势能模型计算该矩形带状可展曲面S中的内应力：

$V=\underset{a=0}{\overset{\infty }{\Sigma }}\frac{w^{2a+1}D}{2^{2a+1}(2a+1)}\int {\left({\eta }^{\prime }\right)}^{2a}({\kappa }^2+2{\tau }^2+{\tau }^2{\eta }^2)du---\left(1\right)$

式中，a代表泰勒展开式的阶数，D代表材质刚度，u表示曲面参数平面上 代表其中一个方向的参数坐标，w表示矩形带状可展曲面S的宽度，η代表母线 在单位切向量d₃上的投影，κ和τ分别表示正则中心曲线r上任意一点处的曲率 和挠率。

3)通过数值模拟计算方法计算矩形带状可展曲面S的运动，具体实施方法 的流程图如图2所示。

a)在一个时间步长内包括如下步骤：

给定t₀时刻带状曲面的坐标作为初始状态然后利用辛欧拉法 对运动方程进行积分，得到更新之后中心曲线位置坐标。然后 通过后处理的形式来保持长度约束c_l，这一步使用快速投影方法来计算拉格朗日 函数的驻点，其中λ_l是拉格朗日乘子。

在求解完运动方程之后，通过对速度施加一个各项同性的耗散来为系统增 加阻尼。切向速度与法向速度可分别按照进行更新，其中μ_t是 切向速度的耗散系数，通常取做0.8-0.9；μ_nt是法向速度的耗散系数，通常取做 0.98。

b)接着准静态地更新母线的运动。这一步包分为无约束优化与约束投影两阶 段。具体地，在一个迭代步中，首先用L-BFGS方法解一个无约束优化问题：

$\underset{(q,\eta )}{\min }V+E_{para},$

记上述问题的最优解为q^*,η^*。接下来，约束的投影是通过求解拉格朗日 函数的驻点得到。该拉格朗日函数具有如下形式：

$L_{q,\eta }=\frac{1}{2}{\mathrm{\Delta q}}^{T}I\Delta q+\frac{1}{2}{\mathrm{\Delta \eta }}^{T}I\Delta \eta +{\lambda }_{{\eta }^{\prime }}{c}_{{\eta }^{\prime }}+{\lambda }_{\eta }{c}_{\eta }+{\lambda }_{s\tau }{c}_{s\tau }$

式中，I是单位矩阵，Δq＝q-q^*,Δη＝η-η^*。

约束投影过程中的一个迭代可通过解耦牛顿法来求解拉格朗日函数L_q,η的 一阶最优条件，即首先通过求解一个线性方程来算出拉格朗日乘子，然后再去 更新未知变量。使用局部投影高斯赛德尔迭代求解乘子。当连续两次迭代的梯 度变化小于1％、约束的范数(不含c_η′)小于10^-6并且c_η′≥0时终止迭代。另外， 当q发生变化时，单位四元数约束c_u可通过显式更新来满足，即有 q←q/||q||。

值得注意的是，在中心曲线上出现曲率退化的区间，即存在u∈[u_a,u_b]，有 κ(u)＝0时，坐标η需定义如下：

$\left(\begin{array}{cc}{\min }_{\eta }{\int }_{u_a}^{u_b}{\eta }^2\left(u\right)du,s.t& c_{{\eta }^{\prime }}\left(u\right),u\in [u_a,u_b]\end{array}\right)$

最后根据曲面的更新后的广义坐标生成网格并渲染结果。

本发明实施例对一条纸带施加不同的外部操作时，模拟后得到纸带的不同 形态如图3至图6所示，可见本发明能够得到较为真实的模拟结果。图3至图6 用户施加的外部操作依次为：a)拉伸并扭转纸带的两端；b)松开并固定纸带两端； c)将纸带一端扭转180°并与另一端相接。d)将纸带一端扭转360°并与另一端 相接。

表1

上表1给出了对应于这四种模拟情形的计算时耗分布表，可见本发明能够 快速地模拟出纸带的形变。如图11所示，显示当纸带一端被扭转180°和扭转 粘接形成一条莫比乌斯带的模拟中，纸带总面积的相对误差较小，说明了本发 明较好的保持住了纸带的不可伸展性。

如图7至图10所示，图7是本实施例的模拟效果，图8是本实施例采用相 容有限元的薄壳模拟方法的模拟效果，图9是本实施例采用非相容有限元的薄 壳模拟方法的模拟效果，图10是本实施例使用弹性杆模拟方法的模拟效果，由 此可显示出本发明的方法较之其他的方法，能更好地保持不可伸展带状物体的 一些特有的几何性质，这是其他方法难以做到的。

通过上述实施例，可以显示出本发明方法在性能和保持物体精确几何性质 等两方面上的优势。