基于多材料混合的物体模拟仿真研究_毕设论文

学士学位论文 基于多材料混合的物体模 拟仿真研究 姓名 欧阳文清 院系 数学科学学院 学号 PB13001120 导师 刘利刚 教授 完成时间 二 一七年五月 University of Science and Technology of China A dissertation for bachelor s degree Multi Material based Simulation Optimization Author Wenqing Ouyang Department College of mathematical Science Student ID PB13001120 Supervisor Prof Ligang Liu Finished Time May 2017 中国科学技术大学学士学位论文 致谢 时光匆匆如流水 转眼便是大学毕业时节 春梦秋云 聚散真容易 离校日期 已日趋渐进 毕业论文的完成也随之进入了尾声 从开始进入课题到论文的顺利 完成 一直都离不开老师 同学 朋友给我热情的帮助 在这里请接受我诚挚的 谢意 首先感谢导师刘利刚教授指导和教诲 是他们把我带到了计算机图形学 的研究领域 刘老师严谨的研究态度及忘我的工作精神 认真细致的治学态度及 宽广的胸怀 都将使我受益终身 感谢张举勇副教授对我的指导 在我遇到瓶颈 的时候 张老师的指点总能让我有茅塞顿开的感觉 感谢柴双明和李建方两位师 兄的指点和帮助 是李建方师兄带我进入科研的大门 是柴双明师兄一直给我 实验和论文的帮助 他们在百忙之中抽出时间帮助我解答我的疑惑 使我受益 匪浅 感谢柴双明师兄之前的工作给予我的帮助和启发 同时还要感谢肖世伟 和董智超两位师兄对我的指点 在我对一些技术上的细节不太清楚的时候 是 两位师兄给了我指点让我能不受这些细节问题的困扰 感谢我的两位组员宿建 平和邹倩芳同学 和你们交流项目的问题我总是受益匪浅 也是你们对项目的 坚定信念推动了整个项目的进展 感谢汪子琦 李弢 曲昊男 周晓岑 周尚 几位同学对我的帮助和支持 每一次和你们的交流都让我感到茅塞顿开 同时 我还要感谢我的同学和其他的老师们在学习和生活上给予了我莫大的帮助 我们 这一路走来 也留下了美好的回忆 感谢科大 感谢一路走过来的兄弟姐妹们 在 最宝贵年华里 是你们伴随着我的成长 最后 感谢我家人一贯的鼓励和支持 你 们是我追求学业的坚强后盾 欧阳文清 2017 年 6 月 25 日 I 中国科学技术大学学士学位论文 目录 致谢 I 目录 III 摘要 V ABSTRACT VII 第一章 绪论 1 1 1 3D 打印简介 1 1 2 多材料的 3D 打印 2 1 3 相关工作 2 1 4 本文的内容和形式 3 第二章 算法原理 5 2 1 问题概述 5 2 2 问题的数学描述 5 2 3 弹力的定义 7 2 4 优化过程 8 2 5 仿真框架 9 2 5 1 仿真基础 10 2 5 2 欧拉隐式积分格式 10 2 5 3 牛顿 拉弗森方法 11 2 6 超弹性仿真 12 第三章 结果 13 3 1 超弹性仿真结果 13 3 2 算法的实验结果 14 第四章 总结和展望 17 4 1 总结 17 4 2 对未来的展望 17 参考文献 19 III 中国科学技术大学学士学位论文 摘要 近年来随着 3D 打印技术的发展 多材料打印逐渐变成了现实 这项技术使 得打印机可以在同一个模型上打印多种材料 从而体现一些特殊的性质 多材料 物体的可计算模型一直是图形学研究的热点 研究者们发明了很多种仿真模型去 仿真多材料物体的形变 本文提出了一个基于多材料物体的优化算法 在输入源模型和目标模型的 情况下 算法通过优化我们事先设定的一些能量 可以计算出合适的材料分布和 稀疏外力 使得在这样的材料分布和外力作用的情况下 源模型可以变形到尽量 接近目标模型 我们的能量包括逼近目标模型 保证外力稀疏 保证体积变动不 大 保证材料分布连续等 由于 co rotational 框架所定义的弹力的导函数没有解析 形式 本文在此处采用超弹性的 Neo Hookean 模型 同时由于 Neo Hookean 模型 的弹力的二阶导数计算起来非常费时 为了算法的速度 本文采用拟牛顿法中的 BFGS 方法去优化整个目标函数 由于目标函数的高度非线性 本文采取先固定 一个变量优化另一个变量 再固定另一个变量优化原来的变量的方法 为了验证 结果的正确性 本文还设计了一个超弹性的仿真框架 关键词 多材料 3D 打印 计算几何 模拟仿真 V 中国科学技术大学学士学位论文 ABSTRACT As the development of the 3D print technology in the recent years multi material 3D printing has became a reality This technology would let printer be able to print at the same object with several materials which can let the object show some special prop erty The computational model of multi material object is always hot issue in computer graphics research researchers invent lots of simulation model to simulate the deforma tion of multi materials object Thisarticleraisesaoptimizationalgorithmbasingonmulti materialsobject Aswe input the source mesh and the target mesh this algorithm can compute proper external force and material distribution by optimize some energy we set in advance in which case the object can deform to approximate the target mesh as well as possible Our ener gy include the approximation of the target mesh the sparsity of external force the little changeinvolume thecontinuityoftheyoung smodule etc Becausetheelasticforcede finedbyco rotationalframedoesnothaveanalyticform weuseNeo Hookeanmodelin hyper elastic model here At the same time because of the large time cost of the compu tation of the second derivation of the elastic force defined by Neo Hookean model we use BFGS of the quasi newton methods to optimize the objective function For the high no linearity of the objective function we fix one variable to optimize another one and then fix another one to optimize the formal one This article has designed a simulation frame using hyper elastic material to prove our result s correction in addition Keywords multi materials 3D printing computational design simulation VII 中国科学技术大学学士学位论文 第一章绪论 1 13D 打印简介 3D 打印 也称快速成型技术 或增材制造 Additive Manufacturing AM 它是一种以三维数字模型文件为输入 使用可粘合材料如粉末状金属或塑料 等 通过逐层堆叠累积的方式来构造物体的技术 传统工艺通常使用切 削 钻孔等方法去制造最终产品 可以称之为 减材 制造 与之不同 3D 打印 也称增材制造技术 是在系统的控制之下 逐层地打印材料 然后再将材料堆 积成型 3D 打印可以制造出一些使用传统方法难以制造的复杂形状的三维曲 面 如图 1 1 三维打印通常是采用数字技术材料打印机来实现的 现在存在 着许多不同的三维打印技术 它们的不同之处在于以可用的材料的方式 并以 不同层构建创建部件 以打印材料的性质和形态区分 可以分成挤压型 线型 粒状型 粉末层型 层压型 光聚合型 挤压型的打印技术包括熔融沉积式 Fused deposition modeling FDM 线型的打印技术包括电子束自由成型制造 Electron Beam Freeform Fabrication EBF 粒状型的打印技术包括直接金属 激光烧结 Direct metal laser sintering DMLS 电子束熔化成型 Electron beam melting EBM 选择性激光熔化成型 Selectivelasermelting SLM 选择性热 烧结 Selective heat sintering SHS 选择性激光烧结 Selective laser sintering SLS 粉末层型的打印技术包括石膏 3D 打印 Plaster based 3D printing PP 层 压型的打印技术包括分层实体制造 laminated object manufacturing LOM 光 聚合型的打印技术包括立体平板印刷 Stereolithography SLA 数字光处理 Digital Light Processing DLP 图 1 1使用传统方法难以制造的三维曲面 图片来源 1 1 中国科学技术大学学士学位论文 1 2多材料的 3D 打印 随着 3D 打印技术的发展 多材料混合的 3D 打印也逐渐成为了现实 多材 料的混合 3D 打印方式能够创造一个本身具有不同属性的产品而无需组装 其 目的是通过减少制造产品的步骤来提高效率 与单一材料 3D 打印相比 它可以 一次制造拥有有多种的功能或物理属性的产品 而不需要再把各种部件组装起 来 多材料混合 3D 打印技术加快了拥有日益复杂部件的产品推向市场的速度 并可以精确计算所需的原材料数量 减少了生产浪费 从美观角度来说 多材 料混合 3D 打印打印出来的模型也更加丰富多彩 图 1 2 和图 1 3 都是多材料混 合 3D 打印的一个例子 图 1 2多变扭曲六边盒子 图 片来源 2 图 1 3彩色的孔明锁 图片来源 2 1 3相关工作 本文研究如何在给定源模型和目标模型的情况下 优化出适合的材料分布 和外力分布 在给出材料分布和外力分布的同时 要保证如下几点 1 在我们计算出的外力和材料分布的作用下 源模型的实际变形结果要尽 量接近目标模型 2 我们所计算出的外力应保持稀疏 否则在施加外力时会比较困难 3 我们所计算出的材料分布应保持连续 否则打印起来十分费时 4 由于我们仅仅使用两种材料 因而每个体素的杨氏模量应为两种材料之 一 在这一节 我们将回顾一下最近的可编程材料的相关进展 我们的工作和最 近的关于多材料打印变形模型 可编程材料 模拟仿真的研究息息相关 2 中国科学技术大学学士学位论文 变形物体的仿真是计算机图形学中很成熟的一个问题 其中有限元方法 Finite Element Method FEM 是最流行的方法之一 3 和 4 给出了这方面 工作的一个概述 为了构建一个视觉效果好的模拟过程 我们必须选择一个合 适的材料模型并且调试好该模型的参数 为了规避调参的复杂性 一些作者选 择从实际物体得到材料参数 很多方法被发明出来用以测量材料的杨氏模量和 泊松比 5 和 6 以及 7 随着高维度的所有四面体的杨氏模量材料分布优化 8 通过力和位移的测量获得了真实物体的非均匀材料分布 研究者们发明了大 量用于构建预期变形行为的数据驱动材料设计方法 9 和 10 本文的工作包 含了最近在多材料 3D 打印方面的进展 11 和 12 和用以构造非均匀空间材料 分布的物体的交互式材料设计 13 本文的工作和 14 的可计算设计有些类似 但也有很多不同之处 例如 本文的的输入模型均以三维体网格作为输入 而 14 的输入是以二维面网格作为输入 除此之外 本文的优化方法也和 14 的 不尽相同 最近一年 15 通过改变非均匀空心模型的厚度设计了一种弯曲变 形方法 实现了期望的弯曲变形 16 介绍了一种允许多材料形变的可编程材 料 1 4本文的内容和形式 本文主要研究如何在给定源模型和目标模型的情况下 优化出适合的材料 分布和外力分布 本文的结构安排如下 第一章为绪论 简要介绍了 3D 打印和多材料混合 3D 打印 并且介绍了本 文的相关工作 第二章介绍本文的问题和相关约束 并且详细介绍了本文的算法和实现 细节 第三章展示了本文的算法制造的结果 第四章是对本文工作的总结以及对未来的展望 3 中国科学技术大学学士学位论文 第二章算法原理 2 1问题概述 本问题的输入是源模型的四面体网格和目标形状的四面体网格 以及两种 材料的杨氏模量 本文的算法也支持多种材料 为了简化起见 这里我们只讨 论两种材料的情况 输出有三个 第一个是源模型的材料分布 第二个是逼近 目标形状的实际变形形状 第三个是变形所需要的外力的位置和大小以及方向 图 2 1 基本说明了本文的工作流程 首先 我们接收两个体网格作为程序的 输入 然后我们利用 example based learning 方法 给出材料分布的初值 并且 设置好实际变形形状的初值 然后我们输出计算结果并将计算结果导入到仿真 框架中 验证我们的计算结果的有效性和合理性 我们需要考虑如下四个约束 1 在我们计算出的外力和材料分布的作用下 源模型的实际变形结果要尽 量接近目标模型 2 我们所计算出的外力应保持稀疏 否则在施加外力时会比较困难 3 我们所计算出的材料分布应保持连续 否则打印起来十分费时 4 由于我们仅仅使用两种材料 因而每个体素的杨氏模量应为两种材料之 一 我们利用在目标函数中加入对应的能量亦或是在优化问题加入等式约束来 使得本问题的输出能够满足我们所指定的四个约束 因为各个约束的能量是互 相制约的 所以调试参数也是本文很重要的工作之一 在实验的过程 我们发现优化过程中内部力应当保持为 0 否则优化的过 程中内部的四面体将会发生翻转 这是模拟仿真所不允许的情况 因为四面体 一旦翻转 就不符合物理现实了 这样有限元仿真就和实际情况相去甚远了 2 2问题的数学描述 首先我们给出一些记号 记输入的四面体网格的边界点集为B 输入的源 模型的所有点的坐标构成一个向量 P R3np 其中 np代表顶点的数目 对应的目 标形状所有点的坐标构成的向量记为 t R3np 我们把实际变形结果的的所有点 的坐标记为 x R3np 记整个模型的杨氏模量分布为 E Rnt nt代表四面体单元 的数目 记T代表整个模型的四面体单元集合 图 2 1我们的工作流程 5 中国科学技术大学学士学位论文 然后我们考虑优化问题需要加入的能量 首先要求杨氏模量 因此加入能 量 ETLE 其中 L 为定义在四面体网格上的标量场 在此处为 E 的拉普拉斯矩 阵 离散拉普拉斯矩阵定义如下 LE i m j 1 i j Ei Ej 2 1 这里 i j 1 如果四面体 i 和四面体 j 是面相邻的 其他情况 i j 0 显然拉 普拉斯矩阵 L 是半正定的而且是非满秩的 L 的零空间是一维的 而且包含均 匀分布 即所有的 Ee都相同的情况 通过计算可以得到 ETLE 1 i 2 来保证优化结果是稀疏的 ksparse为稀疏能量的系数 同时为了保证有限元仿真的精度 如果实际变形结果导致四面体单元的体积 变化过大 那么有限元仿真的实际精度便会下降 因此我们需要变形结果保体 积 因此加入保体积能量 i V i V i prim 2 为了节约成本 我们希望尽可能少地 使用柔性材料 因此加入能量 e Ee E1 2 最后内部点的弹力应当为 0 加入能 量 e e B fe int 2 6 中国科学技术大学学士学位论文 下面考虑约束 唯一的约束就是每个单元的杨氏模量应当是两种材料的杨氏 模量之一 在一开始 我曾经考虑过将内部点弹力为 0 的条件作为等式约束加进 来 后来在实际的计算中发现这样会导致收敛速度过慢 原因很简单 这样做 的话我们加了两个等式约束 注意到弹力同时和实际变形结果 x 以及材料分布 E 都相关 如果使用增广拉格朗日法去优化 想要两个等式约束同时成立 增 广系数必须相当大 而只加入杨氏模量这一个约束的话 我们可以只在优化 E 的时候增大增广系数 而且由于程序使用的模型特殊性 优化 E 的速度相比优 化 x 快很多 后来我发现 即便不把内力为 0 的条件作为等式约束加进去 单 纯增大内部力对应的能量项也能达到不错的效果 所以为了速度 我舍弃了这 个想法 综上所述 这个优化问题可以描述为 min E x a1ETLE a2 x t 2 a3Rsparse fe int x E e B a4 i V i V i prim a5 e e B fe int 2 a6 e Ee E1 2 2 6 s t Ee E1or Ee E2 e T 2 7 其中 ai i 1 2 3 4 5 均为要调试的参数 2 3弹力的定义 给定一个未变形模型和一个变形模型 我们应该如何计算变形模型上每一 点的弹力呢 如果不是大变形的情况的话 我们可以直接使用胡克模型 然而 一般而言 本问题中遇到的实际情况都是复杂的大变形情况 在计算机图形学 中 对于大变形计算弹力的问题 我们一般有两种模型可供计算 一是超弹性 模型 一是 co rotational 模型 简要说明的话 超弹性模型是假设物体有一个应 变势能 当物体发生应变时 物体产生的应力取决于当前能量在此处的负梯度 与之对应 co rotational 模型则是假设物体的每一小段形变都是线性的 然后逐 步变形到目标形状 实际应用中 co rotational 的计算速度比超弹性模型快很 多 因此在模拟仿真中 co rotational 模型有着广泛的应用 在我们的问题中 我们需要弹力对实际变形结果的坐标 x 和杨氏模量 E 的导函数用于优化 使用 co rotational 模型计算这个导函数事实上是不可能的 所以我们只能选择超弹性 模型 既然超弹性模型是假设物体有一个应变势能 那么根据应变势能 更准确 的说是应变势能密度 的定义不同 超弹性模型也有很多种模型 大体上可以 分成几何非线性和材料非线性两类 几何非线性指应变和位移的关系是非线性 关系 材料非线性指应力和应变的关系是非线性关系 下面我们介绍几种超弹 性模型 7 中国科学技术大学学士学位论文 首先还是给出一些记号 下面的应变势能密度定义都是针对单个四面 体单元的 设四面体变形前四个点的坐标为 Xe 0 X e 1 X e 2 X e 3 R3 对应的变 形后四个点的坐标为 xe 0 x e 1 x e 2 x e 3 R3 定义矩阵 d xe i xe 0 i 1 2 3 D Xe i Xe 0 i 1 2 3 那么变形梯度定义为 F e dD 1 柯西格林张量定义为 Ce Fe TFe 拉格朗日格林张量定义为 Ee Ce I 参数 为一阶拉梅常 数 为二阶拉梅常数 Neo Hookean模型 其应变势能密度的定义为 e Fe tr Ce 3 det Fe 2 3 2 det F e 1 2 不可压缩的 Neo Hookean 模型 其应变势能密度的定义为 e Fe tr Ce 3 2 det F e 1 2 SaintVenant Kirchhoff模型 其应变势能密度的定义为 e Fe tr Ee 2 2 tr E e 2 在本文中 我采用的是 Neo Hookean 模型 在给出了应变势能密度的之 后 我们就可以计算单个四面体元的应变势能 We Xe xe V e e Xe xe 然 后将每个四面体元的应变势能加起来就得到整个物体的应变势能 W X x e We Xe xe 这样就可以定义弹力 f x W X x x 从中我们也可以直接得 到弹力的导函数 2 4优化过程 首先我们考虑初值问题 从 3 1 给出的优化问题来看 我们所优化的变 量包括 E x x 的初值可以自然地设置成目标形状 t 但是给出一个合适的 E 的初值并不简单 注意到我们的输入只有源模型和目标形状的体网格 因 此我们要从两个体网格直接导出材料分布 一种自然的想法的是给源模型 中变形较大的四面体单元赋以柔性材料 给变形较小的四面体单元赋以刚 性材料 那么 如何衡量四面体单元的变形度呢 这里我们准备采用参考文 献 17 的 LSD 变形度定义 为了不再累述这些符号的定义 此处我直接引 用 3 2 中的符号 并且假设 E1代表刚性材料的杨氏模量 LSD 变形度的定义 为 L Fe L2 Fe 2det Fe L2 Fe 2det Fe 1 3 withL2 Fe 1 3 Fe 2 F 通过此处变形度的 定义我们可以得到材料分布 E 的初值 然后我们需要考虑整个优化问题的优化方法 这是一个约束优化问题 而 且约束为纯等式约束 在实践中 f 的二阶导数求起来实在太过费时 事实上 求 f 的一阶导数就已经不是一件简单的事情了 调试程序的过程中我发现大量 的时间都浪费在求梯度的计算中 原因在于 f 和 f 的导函数的每个分量计算起 来要上千次运算 即便是并行计算都非常的慢 因此为了速度我们准备采用拟 牛顿法中的 BFGS 方法去优化这个函数 我们准备采用增广拉格朗日法去处理 这个约束优化问题 由于此处我们的优化方法是拟牛顿法 如果使用较为简单 8 中国科学技术大学学士学位论文 的二次罚函数法 当系数较大时 海森矩阵的病态性质会极大地影响我们的优 化效果 此处目标函数的变量有两个 E 和 x 我们采取常用的方法 先固定 E 再优化 x 然后再固定 x 优化 E 算法 2 1 详细描述了我们的优化算法的主要 流程 1Initialize 利用 example based learning 方法给出 E 的初值 并且将目标形状 t 设为变形形状 x 的初值 并设置增广拉格朗日系数 和 的初值 2while d e Ee E1 2 Ee E2 2 1do 3while f1 f2 2do 4f1 f2 5固定 E 优化 x 这一步的优化可以描述为一个无约束优化问题 min x a2 x t 2 a3Rsparse fe int x E e B a4 i V i V i prim 2 a5 e e B fe int 2 6 我们采用拟牛顿法去优化这个能量 7固定 x 优化 E 这一步的优化可以描述为一个无约束优化问题 min E a1ETLE a3Rsparse fe int x E e B a5 e e B fe int 2 a6 e Ee E1 2 e Ee E1 Ee E2 2 e e E1 Ee Ee E2 我们采用拟牛顿法去优化这个能量 8计算 f2 9end 10 E1 Ee Ee E2 11 10 12end 算法 2 1 多材料混合问题的优化算法 2 5仿真框架 当我们计算出所期望的材料分布 实际变形结果 稀疏外力 此处外力即 为边界上弹力的方向相反 大小相等的力 之后 需要输出结果到我们的仿真 框架中用于验证结果的有效性和合理性 本节将详细叙述仿真的理论基础以及 超弹性模型的仿真结果和 co rotational 模型的结果对比 9 中国科学技术大学学士学位论文 2 5 1仿真基础 仿真算法的关键在于牛顿第二定律 f m x 变形得到 x f m 这里 x 代表 质点坐标对时间的二阶导数 此公式可以用于计算质点的加速度 为了便于实 际计算 我们将这个二阶的微分方程分成两个一阶的微分方程 v f x v m 2 8 x v 2 9 这个微分方程的积分形式为 v t v0 w t t0 f t mdt 2 10 x t x0 w t t0 v t dt 2 11 这里 v0 v t0 x0 x t0 仿真本质意义上也就是从 t0开始计算得到每个 t 对 应的 v t 和 x t 因此 仿真 和 对时间积分 这两个词经常可以混用 我 们用差分去逼近 2 10 和 2 11 两个式子 v vt 1 vt t O t2 2 12 x xt 1 xt t O t2 2 13 这里 t 是离散的时间步长 结合式 2 8 和式 2 9 我们得到 vt 1 vt tf xt vt m 2 14 xt 1 xt tvt 2 15 这便是常用的欧拉显式积分格式 由于下一个时间节点的信息可以完全由当前 时间节点的状态的一个显式公式 因此这是显式格式 一个常用的使得该格式 更加稳定的技巧是在式 2 13 的右边用 vt 1代替 vt 欧拉显式积分格式是最 简单的积分方法之一 但是它有一个非常不好的缺点 显式格式只对充分小的 时间步长稳定 我们定性地分析这个问题的原因 显式格式假定在时间步长之 内 力 f 是常量 因此如果时间步长取得过大 那么力 f 甚至可能改变至相反 的方向 这导致能量不守恒 质点在这个积分过程中获得了额外的能量 这样 的过程叠加起来 最终质点可能获得相当大的能量 从而导致格式不稳定 解 决这个问题的方法便是我们即将提到的欧拉隐式积分格式 2 5 2欧拉隐式积分格式 在实时仿真中 例如电脑游戏 一个无条件稳定的格式意味着在任何时间 步长下该格式都是稳定的 这非常重要 一种解决这个问题的方法便是使用隐 10 中国科学技术大学学士学位论文 式积分格式 而计算机图形学中最常用的隐式积分格式便是欧拉隐式积分格式 式 2 12 和式 2 13 的显式积分格式做如下修改即可变成隐式格式 vt 1 vt tf xt 1 m 2 16 xt 1 xt tvt 1 2 17 可以看出 最核心的改变在于在方程的右边使用了下一个时间节点的位置 和速度代替当前时间节点的位置和速度 现在不可能通过这两个方程直接 解出下一个时间节点的信息了 取而代之的是两个以下一个时间节点的位 置和速度为变量的非线性方程 我们定义质量矩阵 M R3n 3n为对角元为 m1 m1 m1 mn mn mn的对角矩阵 将质量矩阵代入式 2 14 和式 2 15 得 到 Mvt 1 Mvt tf xt 1 2 18 将式 2 17 代入式 2 18 得到 Mvt 1 Mvt tf xt vt 1 2 19 这个动力系统是非线性的 因为力和位置的关系一般而言是非线性的 2 5 3牛顿 拉弗森方法 通常解决这样的一个系统的方法是使用牛顿拉弗森方法 这个方法从一个 vt 1的猜测值开始 不断地迭代得到近似值 最后 线性化这个方程用以找到 一个更好的逼近 这个过程重复直到误差降到某个容忍值以下 对实时的应用而言 多次迭代找寻精确值代价太大了 所以一般而言只线 性化一次 并且使用 vt作为 vt 1的猜测值 牛顿第一定律指出在没有外力的作 用下 速度会保持不变 因此这说明 vt是 vt 1很好的一个猜测值 在 xt处线 性化这个方程得到 Mvt 1 Mvt t f xt xf x t tvt 1 2 20 Mvt tf xt t2Kvt 1 2 21 这里 K R3n 3n是力 f 的雅可比矩阵 K 包含 3n 个力的分量对于 3n 个坐标分 量的所有导数 K 是通过当前的质点的位置来计算的 整理上述方程得到一个 标准化的线性化方程 M t2K vt 1 Mvt tf xt 2 22 Avt 1 b 2 23 这里 A R3n 3n 这个方程可以通过共轭梯度法去解出来 右边的 b 则可以通 过已知的信息直接计算出来 11 中国科学技术大学学士学位论文 2 6超弹性仿真 根据 2 3 节所述 在超弹性模型中 弹力的定义为 f x W X x x 2 24 自然的可以定义 2 5 3 中的 K Ki j 2W X x xi xj 2 25 结合 2 5 节我们可以得到超弹性模型仿真的算法 1Initialize x0 v0 2for j 1 j maxtime j do 3compute f 4compute K 5compute b Mvt tf xt 6compute A M t2K 7solve vt 1by M t2K vt 1 Mvt tf xt update xt 1by xt 1 xt vt 1 t 8end 算法 2 2 超弹性模型仿真算法 12 中国科学技术大学学士学位论文 第三章结果 3 1超弹性仿真结果 首先我们展示三个使用超弹性仿真框架仿真得到的例子 用以证明我们的 超弹性仿真框架是正确的 值得注意的是 由于超弹性仿真的速度太慢 以下 的仿真结果中时间步长 t 0 01s 取的非常大 不过由于我们使用的隐式格 式 所以一般而言是没有什么影响的 图 4 1 的例子是固定了左端的顶点 然后施加重力 物体在重力的作用下 弯曲 其中物体为均质材料 由于我在整个系统中加了阻尼 经过足够长的时 间后 物体的弯曲会停 下来 图 4 2 的例子是固定了左端顶点 然后最右端的一层以及最左端的一层小 正方体的材料为刚性材料 其余均为柔性材料 然后在最右端的四个角顶点施 加和边方向相同的外力 由此来导致物体的扭曲 图 4 3 的例子没有固定任何顶点 最左端和最右端的一层小正方体的材料 为刚性材料 其余均为柔性材料 然后在最左端和最右端共计 8 个角顶点上施 图 3 1固定左边顶点 重力作用下弯曲的例子 图 3 2固定左边顶点 在最右边施加 4 个稀疏外力使得物体扭曲的例子 13 中国科学技术大学学士学位论文 图 3 3不固定顶点 施加 8 个稀疏外力使得物体扭曲的例子 加方向相反稀疏外力 由此来导致物体的扭曲 可以看出 仿真结果大体上还是得到了我们想要的结果 可以证明我们的 超弹性仿真框架是正确的 3 2算法的实验结果 在展示完超弹性的仿真结果之后 我们来证明我们的算法的有效性和合理 性 我们使用 co rotational 仿真框架构造了两组输入 这是为了得到的输入并不 是那么好 co rotational 模型和超弹性模型有一些区别 然后将两组输入导入 到我们的算法程序中 然后算法程序输出外力 材料分布 和实际变形结果 我们将外力 材料分布 源模型导入到超弹性的仿真框架中再得到仿真结果 然后对比仿真结果和实际变形结果 图 4 4 展示了我们的算法对扭曲情况的结果 其中的目标形状我是用 co rotational 仿真框架加 8 个外力输入得到 在设置稀疏外力的时候 取边界上最 大的 10 个外力 其中仿真结果那张图中 红色部分代表刚性材料 蓝色部分代 表柔性材料 坐标对比的那张图中 颜色是渐变的 偏蓝色说明坐标和目标形 状差距越大 偏红色说明坐标和目标形状差距越小 这里需要注意我们的算法 仅仅逼近边界点 从图中可以看出 大部分点都是偏红色的 因此我们的算法 还是取得了不错的结果 图 4 5 展示了我们的算法对弯曲情况的结果 其中的目标形状是由 co rotational 仿真框架施加 15 个外力得到的 在设置稀疏外力的时候 取边界上最 大的 20 个外力 与图 4 4 一样 其中仿真结果那张图中 红色部分代表刚性材 料 蓝色部分代表柔性材料 坐标对比的那张图中 颜色是渐变的 偏蓝色说 明坐标和目标形状差距越大 偏红色说明坐标和目标形状 差距越小 图中大部 分点都是偏红色的 故我们的算法的结果还是不错的 14 中国科学技术大学学士学位论文 图 3 4算法对扭曲情况的结果 图 3 5算法对弯曲情况的结果 15 中国科学技术大学学士学位论文 第四章总结和展望 4 1总结 本文聚焦于可编程的材料分布计算问题和多材料物体的模拟仿真问题 本 文的主要贡献有以下几点 1 针对给定的变形结果 本文可以计算出合理的外力分布和材料分布使得 变形结果逼近目标形状 2 本文解决了 14 的 3D 情形的问题 3 本文系统地叙述了仿真的理论基础 并且使用仿真框架来验证本文计算 结果的合理性 4 本文系统地介绍了超弹性材料的性质 并且针对 Neo Hookean 模型完成 了超弹性材料的仿真框架 4 2对未来的展望 可以看出本文给出的实验例子都非常简单 这并非是我不想使用复杂的例 子 而是计算速度所限 对一个非常简单的 376 个顶点的例子 本文的计算程 序需要 5 个小时的时间 而一般的复杂模型有超过 5000 个顶点 此时 BFGS 不 太适用 稠密矩阵的相乘运算会占用大量的时间 但即使使用 L BFGS 方法去 优化 按照本文的算法 时间复杂度是线性增加的 那么需要超过一周的时间 未来的工作应该聚焦在加速程序的运行 超弹性模型的仿真也非常之慢 主要原因在于应变势能的导函数形式实在 太过复杂 未来的工作可以考虑将超弹性模型线性化 用以加速其运算 17 中国科学技术大学学士学位论文 参考文献 1 制造业革命 3D 打印的主流工艺盘点 0130407312949 html 2 国内发布首台三色混合 3D 打印机 3 NealenA M llerM KeiserR etal PhysicallyBasedDeformableModelsinComputerGraph ics 2005 4 Sifakis E Barbic J FEM Simulation of 3D Deformable Solids A Practitioner s Guide to Theory DiscretizationandModelReduction ProceedingsofACMSIGGRAPH2012Courses New York NY USA ACM 2012 20 1 20 50 5 VerstraetenB SermeusJ SalenbienR etal Determinationofthermoelasticmaterialproperties bydifferentialheterodynedetectionofimpulsivestimulatedthermalscattering Photoacoustics 2015 3 2 64 77 6 Becker M Teschner M Robust and Efficient Estimation of Elasticity Parameters us