弹性体的仿真计算研究

学士学位论文 弹?的?计算?弹?的?计算? 作者姓名：倪宁 学科专业：信息与计算数学 导师姓名：刘利刚教授 完成时间：二一八年五月二十日 University of Science and Technology of China A dissertation for bachelors degree Simulation Computation of Elastic Objects Author: Ning Ni Speciality: Mathematics and Applied Mathematics Supervisor: Prof. Ligang Liu Finished time: May 20, 2018 中国科学技术大学本科毕业论文 致谢 时光荏苒，转眼间本科四年的学习生涯已接近尾声。在此，我衷心感谢那些 我学习过程中给予我帮助的同学和老师，使我能够不断进步。 首先，感谢刘利刚老师的指导与教诲，毕业设计过程中刘利刚老师为我提供 了平台和环境，同时也帮助我解决了毕业设计中的许多问题，刘利刚老师深厚的 学术功底，严谨的工作态度使我受益良多。衷心感谢刘利刚老师给予我的悉心教 导和热情帮助。 接下来，感谢欧阳文清学长在毕业设计过程中给予的指点与照顾，在毕业设 计过程中，学长始终不厌其烦，提供了许多的帮助，解决了许多学术上的相关的 疑惑。 同时，也要感谢在计算机图形学实验室一起学习的同学，他们给予了许多学 业上的帮助；在毕业设计过程中，他们与我一同讨论问题，我从中收获颇多。其 中，尤其要感谢孙宇澄同学，在这学期中，他为我提供了很多关心与照顾，孙宇 澄同学严谨认真的科研态度也在激励我不断进步。 最后，感谢我的家人，是你们在背后一直默默支持我，让我能够安心学习进 步。 中国科学技术大学本科毕业论文 目录 中文内容摘要 2 英文内容摘要 3 第一章绪论 4 第一节形变动画模拟 4 第二节相关工作 5 第三节本文的内容 6 第二章算法原理及实现 7 第一节问题的描述 7 第二节问题的数学描述 7 第三节弹力的求解 10 一、内力的计算 10 二、自动微分 10 第四节优化算法 12 一、优化问题的描述 12 二、收敛条件 13 三、一般二次优化问题 14 四、步长的搜索 15 五、优化框架 17 第三章结果与讨论 18 第四章总结和展望 21 第一节总结 21 第二节对未来的展望 21 参考文献 22 1 中国科学技术大学本科毕业论文 中文内容摘要 基于物理过程的动画制作一直是图形学研究的热点，用户通过一些简单的 控制，制作出一些接近真实物理过程的动画。 本文提出了一种控制弹性体变形的方法，用户通过一些目标能量控制弹性 体运动，算法通过优化适当的能量，计算出满足用户需求的动画，同时优化算法 保证优化结果尽可能满足运动方程，即满足真实物理过程。算法输入包括模型 的初始形状及材料分布和用户控制的目标能量，算法优化得到修正后的初始形 状和动画各帧。算法的能量包括目标能量，仿真约束能量和正则项，优化方法为 SQP 方法；对于模型的本构模型则采用 Neo-Hooken 模型，关于弹力的计算则使 用自动微分法。 关?关?：模拟仿真；动画控制；弹性变形 2 中国科学技术大学本科毕业论文 Abstract Physically-baseddeformableanimationisawell-studiedproblemincomputergraph- ics.Users generate physically-based deformable animation by simple control. This article presents a method for generating deformable body animations.Users control the motion of deformable body by objective energy.The algorithm generates the animation which satisfies the equation of motion as much as possible by optimiz- ing proper energy.The inputs of the algorithm includes the initial rest shape of the model,material distribution of the model and the objective energy which is used to con- trol the motion of deformable body.Then the algorithm generates all the frames of the animation and the correction of the rest shape.The energy of the algorithm includes objective,simulation constrain and regularizer.The method of the optimization is SQP method.ThealgorithmusesNeo-hookenmodelandcomputeselasticforcebyautomatic differentiation. Key Words: physically-based simulation;animation control;elastic deformation 3 中国科学技术大学本科毕业论文 第一章绪论 第一节形变动画模拟 在为游戏视频创造可信的虚拟世界中，物体的移动过程在其中扮演着至关 重要的角色，生成较为逼真的动画的问题已被广泛研究。 基于物理过程的模拟提供了一些方法生成形变体形变或运动的动画， 从 Ter- zopoulos1等人在有限元模拟 (FEM) 上的开创性工作开始，许多研究项目旨在 增加物理模拟的多样性，效率及准确性。目前已有许多针对不同物理模型所做 的本构模型建模及仿真工作：针对壳模型,Grinspun2等人提出了离散的壳模型， 用来对诸如帽子，树叶，纸张等较薄的柔性结构进行模拟仿真 (见图1.1, 图1.2)； 对于一维结构，Spillmann3及 Bergou4等人提出了一维弹性杆的仿真模型并推 导了相应的仿真模拟框架（见图1.3） ；对于固体体积单元，Irving5等人提出了 关于不可压缩固体的有限元模拟框架。 图 1.1利用离散壳模型模拟帽子的下落及碰撞（图片来源2） 4 中国科学技术大学本科毕业论文 图 1.2利用离散壳模型模拟海报（图片来源2） 图 1.3弹性杆模拟（图片来源3） 第二节相关工作 本文研究如何在用户给定所需的目标能量函数的情况下，生成尽可能满足 物理规律的运动动画，同时对物体的初始形状进行优化，对初始形状进行一定的 修正，同时在优化的能量中加入了正则项，保证优化的结果更加合理，并保证物 体变形更加平滑， 在这一节，我们回顾一下近年来形变模拟的建模及仿真的相关进展： 基于物理过程的的动画仿真模拟是一个在图形学领域被深入研究的问题。 其 中有限元方法（Finite Element Method, FEM）是使用最为广泛的方法之一。早期 的研究中，如 Terzopoulos1等人使用数值模拟生成物理过程的动画，这些方法 被广泛使用生成布料的合理仿真过程6， 在植物7、人体组织8的模拟中也使用 了这些技术，这些技术常常用来在游戏和电影中增加现实感。另一方面，在人体 或动物的模拟仿真中，已经能够在变形体内嵌入刚性骨架，同时 Hahn9等人提 5 中国科学技术大学本科毕业论文 出一些更加鲁棒的方法用来模拟骨骼与软组织间的相互作用。关于对形变体的 控制，Twigg et al.10提出了一种方法利用多个物体的模拟来生成不同场景下的 动画，并让用户通过交互来选择理想的结果；Popovi et al.11提出了一个系统来 自动处理各种模拟参数，如质量、惯性矩、弹性系数等，以便用户生成满足用户 需求的简单的动画；Martin et al.12建立了一个框架，借由实例，让用户引导物 体运动并与环境交互，从而生成满足需求的动画。 第三节本文的内容 本文研究在用户给定目标及模型初始形状的情况下，优化出模型合适的初 始形状并输出模拟动画各帧，本文结构安排如下： 第一章为绪论，简要介绍形变动画模拟，并介绍了模拟仿真的相关工作； 第二章介绍了本文的具体问题，并详细介绍了本文的算法和实现细节； 第三章对算法的结果进行了展示并作了相关讨论； 第四章是对本文工作的总结及对未来的展望。 6 中国科学技术大学本科毕业论文 第二章算法原理及实现 第一节问题的描述 问题的输入包括模型的参数及用户的目标，模型的参数包括模型的形状及 模型的材料分布。其中模型的形状通过四面体网格表示， 材料分布则通过给出每 个四面体内的材料参数表示，其中认为每个四面体只由一种材料填充，材料参数 由杨氏模量及泊松比表示；用户的目标则给出适当的目标函数，具体见下文。算 法的输出包括生成的逼近用户目标的动画，同时算法还对用户给出的初始形状 进行了适当修正，使得结果更加逼近用户的目标。如图2.1所示，用户希望让最 左侧的物体发生形变，经理中间的过程，得到最右侧的形状，且整个过程的运动 尽可能满足物理规律。 图 2.1问题目标形变过程（图片来源13） 第二节问题的数学描述 首先我们定义一些记号：记 npoint为初始模型的四面体网格中点的个数，记 X0为用户给出的模型的初始构型中所有点的坐标形成的 3 npoint维向量，记 X 为经过一定修正后的初始构型中所有点的坐标形成的 3 npoint维向量，记 动画模拟帧数为 nframe，记 xi，i 1,2,.,nframe，为第 i 帧中各点组成的 3npoint维向量，记 x(j) i 为第 i 帧中 j 个点的坐标，记 x 为所有帧的点的坐标组 成的 3 npoint nframe维向量。 7 中国科学技术大学本科毕业论文 然后我们考虑优化问题中的能量：首先是用户目标能量的给出，算法要求输 入的目标能量需满足如下的形式： Eobj= 1 2 nframe i=1 j ijx(j) i si2,(2.1) 其中，ij为第 i 帧中第 j 个点的权重，si则为第 i 帧所要逼近的目标。例如，用 户需要让物体作简单的知质心平移，可以取权重为 ij= mi i mi，取 si 为质心在 每一帧的目标位置；同时，目标函数也可以取多个，例如对于一些情况，比如用 户希望一部分点做一些运动，可以认为这些点每个点的运动都是一个目标；比如 对于标号为 k 的点的运动，可以在 j = k 时取 ij= 1，其他情况为 0，si则为该 点在每帧所要逼近的目标位置。对于多目标，例如模型的一部分点固定，一部分 点运动，可以将每个点的运动目标简单相加或加权相加，从而得到总共的目标能 量函数，表达式即为如下形式： Eobj= 1 2 nframe i=1 kU x(k) i s(k) i 2,(2.2) 其中，U 为控制点标号的集合，s(k) i 为第 k 个点在第 i 帧所要逼近的位置。 其次，我们希望模拟过程尽量满足物理规律，即满足如下运动方程: M x fert fint= 0,(2.3) 其中，fert和 fint分别为模型所受外力和内力，在本文中，内力为模型形变产生 的弹力；对时间离散得下式: M xi+1 2xi+ xi1 h2 fi ert f i int = 0,(2.4) 其中 h 为时间步长。优化过程中，我们希望式 (2.4) 尽可能成立，取仿真约束如 下形式： Esimulate= 1 2 nframe1 i=2 M(xi+1 2xi+ xi1) h2(fert+ fint(x,X)2,(2.5) 8 中国科学技术大学本科毕业论文 式中外力与内力均为3npoint维向量， 但随时间改变， 即每帧可以为不同的值， 在 本文中为讨论方便取外力恒为0。 在上述仿真约束过程中， 我们利用 xi12xi+2xi+1 h2 来当做对加速度的估计，这里实际上为显示的数值格式，可以将外力及内力取为 第 i + 1 帧的值，这样的计算结果为隐式法的数值格式，使方法在数值计算上更 加稳定。 最后，只有以上两个能量问题自由度仍过高，故加入正则项使优化能够更好 的收敛。首先，我们希望对用户的模型形状的修正不能过大，加入如下正则项: Eregular1= X X02,(2.6) 这里的正则项也可以取其他的函数，例如可以采用如下正则项： Eregular1= W(X0,X),(2.7) 其中 W 为以 X0为初始形状在 X 处的应变势能。 其次，本文使用 Neo-Hooken 作为本构模型来计算内力，优化过程中会出现 四面体网格翻转的情形，而一旦网格翻转，就导致错误的内力计算且优化很难收 敛，故希望网格中每个四面体体积变化不能过大，故加入如下正则项: Eregular2= i,j volj i volj 0i 2, (2.8) 其中 volj i 为第 i 帧第 j 个四面体的体积，volj 0i 则为相应的初始形状的体积，这 里的体积为带符号的体积，对于给定四面体（设其顶点为 ai,1 6 i 6 4） ，其体积 为向量 ai a4,1 6 i 6 3 的混合积的 1 6。 综上，得到我们优化的能量有如下形式： E(x,X) = k1Eobj+ k2Esimulate+ k3Eregular1+ k4Eregular2,(2.9) 而我们的问题则描述为下式: argmax x,X E(x,X),(2.10) 其中，ki，i 1,2,3,4 为能量的权重，需要根据具体问题调整。 9 中国科学技术大学本科毕业论文 第三节弹力的求解 一、内力的计算 本文物体内力计算假设物体为超弹性模型，采用 Neo-Hooken 模型。超弹性 模型通过定义物体内每点的应变能密度， 以此来表达应变与应力间的关系。本文 采用如下应变能密度如下所示： = 2 (I1 log(I3) 3) + 8 log2(I3),(2.11) 其中 I1= tr(FTF)，I3= (detF)2，F 为形变梯度， 与 与材料有关的参数，其 中 = E 2(1+)， = E (1+)(12)，这里 E 与 分别为杨氏模量与泊松比。 对于离散四面体网格，对于一个四面体单元，设其初始位置四个点坐标为 Xe 1,Xe2,Xe3,Xe4 R3，对应变形后点的坐标为 xe 1,xe2,xe3,xe4 R3，定义矩阵 D = Xe 1 Xe 4,Xe2 Xe 4,Xe3 Xe 4，d = xe1 xe 4,xe2 xe 4,xe3 xe 4，此处假设四面 体形变均为线性，故有 F = dD1，则对一个四面体单元可以定义出其应变势能 W(xe,Xe) = vole(xe,Xe)，对每个四面体相加即可得到整个模型的应变势能， 每个点弹力即可定义为：f(x) = W x 。 二、自动微分 除了对力的计算外，另一方面，优化过程中需要得到内力的梯度和 Hessian 阵，就需要对力求一阶导数和二阶导数，也等价于对所定义的应变势能求二阶导 和三阶导。求导过程可以手动求解或借助符号计算软件， 但仍过于繁琐且容易出 错，特别是 Hessian 阵的计算，故本文使用自动微分法来完成对导数的求解。 不同于数值微分强调直接数值近似求解，也不同于符号微分通过求解出导 函数具体表达式再带入数值计算，自动微分将符号微分应用于一些基本的算符， 如几个初等函数以及基本运算，在运算时保留中间结果，最后求得最终结果。自 动微分的原理为基于以下的观察，我们用到的大多函数为一些较为简单的函数 10 中国科学技术大学本科毕业论文 的复合或简单运算，其微分原理为求导链式法则。本文通过运算符重载实现自动 微分。 对于一阶导数的求解，设标量函数 F : x Rn7 y R 和 G : x Rn7 y R，我们可以定义二元组 u,du，其中 u 为函数值，du Rn为函数在所给 点 x Rn的梯度，我们定义二元组的运算： F,dF + G,dG = F + G,dF + dG,(2.12) F,dF G,dG = F G,G dF + F dG,(2.13) F,dF/G,dG = F G, G dF F dG F2 ,(2.14) 同时对于一些初等函数也可以类似定义，如 sinF,dF = sin(F),cos(F) dF。 对于复合函数 y = f(x1,x2,.,xn)，对其求导，通过链式法则有： f = f x1 x1+ f x2 x2+ . + f xn xn(2.15) 若 xj(1 6 j 6 n) 为自变量则有 xj= j i，其中 j i = 1,ifi = j 其他情形为 0。 至此对于一个复杂的函数，要求给定点函数值，可以通过对自变量做简单运算或 初等函数运算，再加上函数的复合得到函数值，而计算过程中，如果用以上二元 组计算，则除了能够得到函数值，也能够同时得到函数在给定点处梯度值，这里 运算均为数值运算。 对于二阶导数即 Hessian 阵的求解，也可以类似通过类似方法求解，通过链 式法则有： f = f(x1,x2,x3,.,xn),(2.16) f = fi xi,(2.17) f = fi xi+ fij xi xj,(2.18) . f = fi.xi+ 3fij xi xj+ fijk xi xj xk,(2.19) 11 中国科学技术大学本科毕业论文 这里采用 Einstein 求和约定，省略了求和记号，fijk= 3f xixjxk。由以上公式可 仿照一阶导数求法，得到函数的二阶导即 Hessian 阵，由 Hessian 阵对称，故求 解过程中只需求解上三角或下三角部分。 若要使用自动微分法求解更高阶的导数，将不再能使用如上方法，因为随着 导数阶数增加，链式法则求解过程的表达式会的复杂度将会大大增加，从而导致 计算效率大大降低，这时则需要利用截断的泰勒展式，实际上，函数的泰勒展式 第 n 项对应着其第 n 阶导数的信息，通过对复合函数的泰勒展开运算，能够求得 最终函数的泰勒展式。设 X(t) = p i=0 xiti，Y (t) = p i=0 yiti，Z(t) = p i=0 ziti，则有: Z(t) = X(t) + Y (t),zi= xi+ yi(2.20) Z(t) = X(t) Y (t),zi= xi yi(2.21) Z(t) = X(t) Y (t),zi= i j=0 xj yij(2.22) 由以上运算，及一些基本函数的泰勒展式，及函数复合即可得到最终函数的导 数。 第四节优化算法 一、优化问题的描述 对于一般非线性优化问题，可以表述为： min xRn f(x)(2.23) 满足如下约束： gi(x) 6 0,fori = 1,.,m1 gi(x) = 0,fori = m1+ 1,.,m (2.24) 其中 gi(x) 为约束项。 12 中国科学技术大学本科毕业论文 本文采用 SQP 方法优化上述非线性问题， 优化算法可以概述为如下过程， 对 给定初值点 x0，执行以下算法： 1. 检查 xk是否满足收敛条件，若满足则终止，否则转到 2。这里 xk为第 k 步计算得到的 x 的值。 2. 在点 xk处用一个二次的子问题来逼近原问题，用该二次子问题的解 dk 作为搜索的方向，转到 3。 3. 使用线性搜索的方法来评估函数值，并确定搜索步长 k，转到 4。 4. 更新点的位置，得到 xk+1= xk+ kdk，转到 1。 可简要概述为每一步通过一定方法，求解下一步的搜索方向，再确定搜索的 步长，从而得到下一步点的位置。 二、收敛条件 对于带约束的优化问题，可以通过拉格朗日函数将其转化为无约束的问题， 即： L(x,) := f(x) + Tg(x),(2.25) 其中， Rm为拉格朗日乘子。 对于问题的收敛条件采用 KKT 条件判别： (a)xL(x,) = 0, (b)gi(x) 6 0,i = 1,.,m1;gi(x) = 0,i = m1+ 1,.,m, (c) igi(x ) = 0,i = 1,.,m1, (d) i 0,i = 1,.,m1, (2.26) 其中 (x,) 为优化问题的一个局部最优解，称为 KKT 点；这里 KKT 条件为最 优化的必要条件，但不是充分条件。 我们称满足优化约束条件的点为可行解，可行解的集合为可行域。在 KKT 条件中，(a) 为拉格朗日函数取得局部极值的必要条件，(b) 为原问题的等式约束 13 中国科学技术大学本科毕业论文 和不等式约束； (c) 为松弛互补条件， 简单来说， 若没有不等式约束的问题极值点 已经满足不等式约束， 则该不等式约束相当于没有， 即 = 0， 否则， 加入约束后， 问题的极值点满足h(x) = 0， 这里设h(x)为约束， 即不等式约束转变为等式约束， 总之， (c) 一定成立； 而 (d) 与 (c) 相似， 若没有不等式约束的问题极值点满足不等 式约束则 = 0，否则，极值点满足 h(x) = 0，设区域 A = x Rn: h(x) 6 0， 故在极值点处，h(x) 梯度指向 A 的补集，而 f(x) 的梯度指向 A，又由极值点满足 等式： xf(x) + xh(x) = 0,(2.27) 可得 0。 三、一般二次优化问题 我们首先定义一般形式的二次优化问题，从而再将二次问题应用于 SQP 方 法。对于 Q Rnn，A Rm1n，C Rmm1n，c Rn，b Rm1，d Rmm1， 一般形式的二次规划问题为： min xRn 1 2x TQx + cTx, (2.28) 满足下式： Ax = b,(2.29) Cx 6 d.(2.30) 对于只有等式的二次规划问题，比较简单，可以使用 KKT 条件判别法，即 将该问题通过拉格朗日函数转化为无约束问题并进行求解，若没有不等式约束， 我们有： L(x,) = 1 2x TQx + cTx + T(Ax b), (2.31) 14 中国科学技术大学本科毕业论文 由等式 KKT 条件： QTx + c + AT Ax b = 0(2.32) 得到关于 x 的线性方程组，只需要求解该方程组即可得到原问题的解，求解方 程组可以直接求解，如高斯消元法或 LU 分解，对于规模较大的问题可以使用迭 代法求解。对于带不等式约束的二次规划问题，则更为复杂，一般通过有效集法 求解，大致思路是在每一步迭代中，确定哪些不等式约束是无效的，即不加上该 约束不影响这步迭代结果，这时可以不考虑该约束；而对于那些有效的约束，在 在该步迭代时将其作为等式约束考虑；由于本文只涉及到等式约束，在此不再赘 述。 对于我们的优化问题，在给定点 xk处，我们可以对其泰勒展开，并取到其 二次项, 原问题转化为： min p f(xk) + f(xk)Tp + 1 2p T2f(xk)p (2.33) 满足约束： gi(xk) + gi(xk)Tp 6 0,fori = 1,.,m1 gi(xk) + gi(xk)Tp = 0,fori = m1+ 1,.,m (2.34) 由此，我们可以将原来的非线性优化问题转化为一系列二次规划问题，以上问 题，可以由二次规划算法，计算得到下一步搜索的步长的方向。 四、步长的搜索 每一步中，由上述二次子问题，能够求解出迭代方向，而确定具体结果，则 还需要确定步长，为评价搜索得到的步长是否合适则需要引入容差评价函数，我 15 中国科学技术大学本科毕业论文 们将增广拉格朗日函数作为评价函数，有: La(x,;) = f(x) + m i=m1+1 igi(x) + 1 2 m i=m1+1 ig2 i(x) +1 2 m1 i=1 1 i (max0,i+ igi(x)2 2 i), (2.35) 也可以取： La(x,;) = f(x) + m1 i=1 imax(0,gi(x) + m i=m1+1 |gi(x)|.(2.36) 为确定每一步的搜索步长，我们定义： () := La(xk+ dk,k;k),(2.37) 如果 dk为评价函数下降方向，则必存在 0 使得 () 6 (0)，事实上，由 函数的连续性，我们可以知道只要 足够小，就能够使得上述不等式成立；但 同时，我们不希望搜索步长过小，否则评价函数下降也很小，且实际上迭代后的 点与原先的点相差不大，也会导致收敛速度过慢或难以收敛。这里使用 Armijo 准则来完成对步长的搜索，通过 Armijo 准则，能够得到一系列不断减小的步长， 从中我们选取第一个满足 Armijo 准则的步长，作为我们的结果。 对给定的实数 max (0,1， (0,1)， (0,1) 和自然数 K，可以定义数 列： i:= maxi,(2.38) 这里 i 0,1,.,K，通过对 i 从 0 开始判断，直到 i 满足如下不等式： (i) 6 (0) + i(0),(2.39) 由此可以得到迭代的步长。除了 Armijo 准则，也可以使用其他方法进行步长搜 索，例如滤子法，同时对于具体问题，可以先用 Armijo 准则判别，若搜索步长 失败，再使用滤子法搜索。 16 中国科学技术大学本科毕业论文 五、优化框架 综上，我们可以得到我们问题的大致框架： 1. 用户给出模型网格初始位置 X0及材料分布，同时给出 Eobj； 2. 设定优化初值 x，X，默认情况下每一帧均与初始位置 X0一致，X 也设 为 X0； 3. 计算能量各项的值、梯度和 Hessian 阵，这里关于仿真约束项的求导的计 算通过自动微分完成； 4. 将问题近似为二次规划问题，求解得到搜索方向 d； 5. 由 Armijo 准则进行线性搜索，得到合适的步长，从而得到迭代后的 x 和 X 值； 6. 判断结果是否满足 KKT 条件，若满足则结束迭代，否则重复 3、4、5 步。 17 中国科学技术大学本科毕业论文 第三章结果与讨论 首先这里展示几个通过对方块状弹性体作控制得到的例子。在弹性体的设 置中，取杨氏模量为 5000Pa，泊松比为 0.3，这里简单的取模型材料分布只有一 种材料；弹性体形状为长宽高分别为 1m、0.6m、0.6m，划分得到四面体网格，四 面体网格点数为 112 个，四面体个数为 270 个；取帧数为 40 帧，时间步长取为 0.08s，这里帧数不宜过大，否则会使得优化变量过多，从而使得优化过程极慢， 也不宜过少，否则在给定时间步长下，物体实际物理过程中理论上变形量应该很 小，达不到用户给的目标，可能也会在优化过程中出现问题；仿真约束中，可以 使用隐式法来求解，本文中使用显示法，也并无太大问题。 图 3.1弹性体的弯曲 首先为验证优化框架的正确性，对质心平移做了验证，即用户每一帧指定质 心的位置，算法能够 2 步迭代出结果，即得到模型的平行移动，因结果较为简单， 故不再展示。如图3.1所示，通过将模型一面固定，另一个面作一定方向的移动， 使得模型弯曲，以此为目标函数，从而得到动画；这里需要注意的是移动步长不 宜调整过大，否则优化过程中难以收敛，实际中发现，优化中会出现网格中部分 四面体翻转的情形，从而使得计算的力发生错误；若想对模型进行较大步长的平 移，可以将其划分为多个小位移的平移，以上个优化的结果作为下一个优化的处 置。 18 中国科学技术大学本科毕业论文 图 3.2弹性体的拉伸 图 3.3弹性体的弯曲 19 中国科学技术大学本科毕业论文 图3.2中展示的是模型的拉伸，取目标函数为下面的点固定，另一面点做一 定的平移运动；由图中可见，模型的拉伸能够做到不错的效果，基本能够与实际 物理效果相符。 图3.3中展示的是另一个模型的弯曲，取目标函数为下面的点固定，另一面 点做一定的平移运动；由图中可见，模型的基本弯曲能够做到不错的效果，但关 于初始形状的调整则使得模型看起来不是太光滑；这里有几种解决方法，一种是 将网格的剖分加细，用更好的网格来优化则能够得到更好的结果，例如能够得到 如图3.1的效果，一种是加大关于初始形状的正则项的系数，使得模型的形状变 化不宜过大，但该方法的效果有限，且会使得模型几乎不会修改形状；另一种方 法为通过 Cage 的结构，通过 Cage 上少部分点的坐标控制所有点的坐标，选取合 适的坐标可以使得原模型的形状修改更加平滑。 20 中国科学技术大学本科毕业论文 第四章总结和展望 第一节总结 本文聚焦于在用户给定控制下，生成尽可能满足物理规律的动画过程，本 文的主要内容有：借由用户的控制，可以生成基于物理过程的动画，并对初试 模型的形状作一定修正，使得结果逼近的更好；本文还系统的叙述了非线性优 化的 SQP 方法的理论基础，并介绍了自动微分法，从而使得我们能够在对复杂 函数优化时，方便的求得其梯度及 Hessian 阵；另外，本文也介绍了一些关于 Neo-Hooken 模型的内容。 第二节对未来的展望 本文中所给的例子非常简单，除了模型简单外，动画形变过程也较为简单， 同时模拟帧数也较少，这里很大的原因是计算速度所致；另外，用户控制的形变 不能过大，否则会导致优化不能收敛或收敛到错误的结果，因而制作大形变动 画需要将大形变进行分割，进行多次优化；这里优化失败的很大原因是所取的 Neo-Hooken 模型的弹性势能不保翻转，即模型的点翻转后，能量不能保证点能 够翻转回来，这可能需要寻找更好的弹性势能。另外，关于使得模型的初始形状 的优化也需要改进，可以加入 Cage 来利用部分点控制所有网格，一方面可以降 低问题维度，另一方面也可以使得优化的结果更加平滑，这里可以参考 Sanocki et al.14中利用调和坐标控制模型的工作。未来的主要工作还是主要集中在优化 方法的改进及优化能量的选取。 21 中国科学技术大学本科毕业论文 参 考 文 献 1Terzopoulos D. Elastically deformable models. Computer Graphics, 1987: 205-214. 2Grinspun E, Hirani A, Desbrun M, et al. Discrete shells/Proceedings of the 2003 ACM SIG- GRAPH/Eurographics Symposium on Computer Animation. 2003. 3Spillmann J, Teschner M. C o r d e: Cosserat rod ele