自然物体有各自鲜明、千姿百态的力学特性.如何逼真、高效地实现物体的物理仿真, 一直是图形学的重要研究方向之一, 此课题在交互式的虚拟现实、影视动画、游戏等领域有着重要和广泛的应用.其中, 弹性形变体的形变行为的模拟是物理仿真领域中一个非常重要的分支^[1,2].形变体的形变效果很大程度上取决于材质的应力应变关系, 即, 本构模型的选择.给定一种材质以及合适的参数, 就能够模拟出相应的运动过程, 因此, 材质属性是形变体仿真方法中非常重要的参数.

材料学为图形学形变体仿真提供了许多基本材料的本构方程, 包括各种线性模型以及非线性的材质模型, 如StVK(saint venant-kirchhoff), NeoHookean, Mooney-Rivlin等.除此之外, 物质世界还存在这样一类物质, 它们是由混合材质构造出来的, 比如混凝土等.为了更逼真地模拟这类物质的力学特性, 图形学研究者们提出了非均质材质的构造方法, 以获取混合材料物体的仿真模型.

从形变行为上分析, 基本的本构模型虽然可以模拟出现实世界中简单形变体的运动特征, 但由于其应力应变关系是恒定的, 所能调节的物质参数有限, 从而导致仿真的形变体只具有特定的运动规律, 形变简单, 无法获取丰富的视觉效果, 不利于游戏、影视动画等实际领域的应用.由此, 为获得更加复杂的本构模型, 依赖于传统材质, 产生了一系列研究方法, 如直接调整参数曲线改变传统的材质模型受力后的形变趋势.这种调节参数曲线的方式一定程度上提高了传统材质的形变范围, 但形变还是遵循运动规律, 不具有艺术指导效果, 无法获得夸张的形变效果.为了获得更自由的形变以满足设计师的动画模拟需求, 研究人员提出直接从设计形变效果的角度出发, 通过给定关键帧改变传统材质模型的形变轨迹, 使仿真模型能够获得更生动的动画效果.

本文第1节阐述形变体仿真中普遍采用的均质材质的几种基本本构模型的基本特性.第2节针对非均质材质的仿真方法进行归类, 梳理目前已有的研究思路和方法.第3节针对可编辑材质分类讨论, 一类是对材质的属性和结构的编辑方法, 另外一类通过编辑形变体的形变效果获得新的材质.第4节回顾现有的研究工作, 并针对几类典型的本构模型作对比分析.第5节对全文总结, 探讨在形变体仿真中材质本构模型的应用与编辑领域的研究方向.

特别需要指出的是, 本文侧重于讨论弹性形变行为, 即三维形变体在虚拟仿真中材质本构模型的应用与设计, 而不讨论当固体形变超过某个阈值时物体所体现出的塑性行为或者碎裂等现象.

1 均质材质

虽然现实世界的模型多是非均质材质组成的, 但均质材质的形变体仿真更容易建模, 处理起来相对容易, 因此在图形学现有的研究中, 对均质材质模型仿真方法的研究占有很大比重.弹性材质模型在形变过程中应力所作的功仅与初始状态和最终状态有关, 而与形变过程无关.其本构关系描述的是固体在弹性形变时应力与应变之间的关系, 通常表示成形变体内单位体积存储的应变能或者弹性势能与应变关系, 这里列举了一些在图形学领域常用的均质材质的本构关系.

1.1 线性弹性模型

在所有材质模型中, 线性模型是最简单的本构模型, 应力张量和应变张量表示为线性关系σ=C:ε, σ, ε分别表示柯西应力和柯西应变张量, 在三维情况下为二阶张量, 将其分别表示成向量的形式为

$ \varepsilon = {[{\varepsilon _{11}}\;{\varepsilon _{22}}\;{\varepsilon _{33}}\;{\varepsilon _{12}}\;{\varepsilon _{23}}\;{\varepsilon _{31}}]^T} $

(1)

$ \sigma = {[{\sigma _{11}}\;{\sigma _{22}}\;{\sigma _{33}}\;{\sigma _{12}}\;{\sigma _{23}}\;{\sigma _{31}}]^T} $

(2)

C为弹性张量, 是四阶张量, 包含了物质的81个材质参数, 由于应力张量和应变张量的对称性, C可简化为具有36个独立参数, 对应公式(1)、公式(2) 表示成矩阵形式:

$ C = \left[{\begin{array}{*{20}{c}} {{C_{11}}}&{{C_{12}}}&{{C_{13}}}&{{C_{14}}}&{{C_{15}}}&{{C_{16}}}\\ {}&{{C_{22}}}&{{C_{23}}}&{{C_{24}}}&{{C_{25}}}&{{C_{26}}}\\ {}&{}&{{C_{33}}}&{{C_{34}}}&{{C_{35}}}&{{C_{36}}}\\ {}&{}&{}&{{C_{44}}}&{{C_{45}}}&{{C_{46}}}\\ {}&{}&{}&{}&{{C_{55}}}&{{C_{56}}}\\ {Sym}&{}&{}&{}&{}&{{C_{66}}} \end{array}} \right] $

(3)

在最简单的情况下, 当形变体由各向同性材质构成时, C只具有2个独立参数.通常用杨氏模量E和泊松比ν表示, 杨氏模量描述材质的刚性, 泊松比反映了形变体积变化, 取值范围为(-1, 0.5).根据对称性, 将应力张量和应变张量表示为向量, 则线性弹性模型的本构关系表示为

$ \left[{\begin{array}{*{20}{c}} {{\sigma _{11}}}\\ {{\sigma _{22}}}\\ {{\sigma _{33}}}\\ {{\sigma _{12}}}\\ {{\sigma _{23}}}\\ {{\sigma _{31}}} \end{array}} \right] = \frac{E}{{(1 + \nu )(1 - 2\nu )}}\left[{\begin{array}{*{20}{c}} {1-\nu }&\nu &\nu &0&0&0\\ {}&{1-\nu }&\nu &0&0&0\\ {}&{}&{1-\nu }&0&0&0\\ {}&{}&{}&{1 - 2\nu }&0&0\\ {}&{}&{}&{}&{1 - 2\nu }&0\\ {Sym}&{}&{}&{}&{}&{1 - 2\nu } \end{array}} \right]\left[{\begin{array}{*{20}{c}} {{\varepsilon _{11}}}\\ {{\varepsilon _{22}}}\\ {{\varepsilon _{33}}}\\ {{\varepsilon _{12}}}\\ {{\varepsilon _{23}}}\\ {{\varepsilon _{31}}} \end{array}} \right] $

(4)

另外一种形式的本构关系用应变能量密度表示, 线性本构关系的应变能量密度为

$ \mathit{\Psi }(F) = \mu {\varepsilon _C}:{\varepsilon _C} + \frac{\lambda }{2}t{r^2}({\varepsilon _C}) $

(5)

其中, F表示形变梯度; ${\varepsilon _C} = \frac{1}{2}(F + {F^T}) - I$表示柯西应变张量, I为单位矩阵; μ, λ表示莱姆系数, 用杨氏模量E与泊松比ν表示为$\mu = \frac{E}{{2(1 + \nu )}}, \lambda = \frac{{E\nu }}{{(1 + \nu )(1 - 2\nu )}}.$

1.2 共旋的线性弹性模型

因为计算效率高, 图形学中大量使用柯西应力表示的线性模型.但线性模型描述的应力应变关系是线性的, 所以它只能描述小形变, 无法表达较为复杂的运动, 且该模型不具有旋转不变性.

2002年, Müller等人^[3]提出了刚度扭曲(stiffness warping)方法, 其核心思想是:利用物体形变的旋转分量, 将形变后的物体旋转到物质空间, 得到的网格与初始形状计算位移, 从而获得线性力, 然后再将旋转分量作用到物体的运动中.线弹性模型中, 四面体单元的弹性力表示为刚度矩阵与位移之间的线性关系f=K(x-X), 其中:K∈R^12×12表示刚度矩阵, 是一个常量矩阵; f, x, X∈R¹²分别表示四面体单元4个顶点的弹性力、形变后和初始时刻位置向量.根据刚度扭曲方法计算的弹性力模型为f=RK(R^-1x-X).其中, R∈R^12×12表示形变体四面体单元上各个顶点构成的旋转矩阵.他们提出了两种方法提取每个仿真单元结点的旋转变换:第1种是根据每个仿真结点处形变前后的方向变化构造一个变换矩阵, 然后对该矩阵进行SVD分解, 得到变换中的旋转分量; 第2种是构造每个仿真结点处形变前后的局部标架, 根据标架的变化得到仿真结点的旋转变换.由于提取的是每个仿真结点的旋转变换, 因此每个单元的内力合力不为0, 会产生虚拟力(ghost force).随后的研究中, Müller等人^[4]通过对每个单元的形变梯度的极分解得到以单元为单位的旋转变化, 可以有效避免这个问题.但这种通过矩阵分解得到旋转变换的方式存在不稳定性因素:当物体的形变较大时, 待分解的矩阵可能条件不够好, 导致分解失败.因此, Georgii等人^[5]从形状匹配的角度为每个单元定义待求旋转变换并反作用到形变后状态, 利用形变后与初始状态的距离定义的一个二次能量函数, 最小化该能量函数得到每个单元的旋转变换.这样得到的旋转变换无需进行矩阵分解, 比之前的方法更加稳定, 并且也不会产生虚拟力.

共旋的线性弹性模型(corotational linear elasticity)^[6-9]源于这种刚度扭曲的思想, 该模型既保持了柯西应变和线性弹性模型速度快的优点, 还能近似模拟大形变.利用形变梯度的极值分解F=RS, 抽取出形变梯度F的旋转分量R, 创建一个新的线性形变度量ε_C=S-I, 从而解决了普通的线性模型不能模拟大形变的问题.替换掉公式(5) 的应变张量, 得到共旋的线性模型的能量密度表示为$\mathit{\Psi } (F) = \mu ||S - I||_F^2 + \frac{\lambda }{2}t{r^2}(S - I).$

1.3 非线性模型

共旋的线性弹性模型虽然可以一定程度上改善线性形变的不足, 但遇到大的拉伸形变或者压缩形变时, 还是受限于其线性的材质属性.若要更好地描述材质的非线性形变行为, 需要利用非线性的本构关系.

在非线性模型里, 有两类非线性:一种叫做几何非线性模型, 描述的是位移和应变张量的非线性关系; 另外一种是材质的非线性, 描述的是应变张量和应力张量的非线性关系.常用的应变张量有两种:一种是前文描述线性模型时提到的柯西应变张量, 另外一种是格林应变张量.根据格林应变张量的表达式, 其具有旋转不变性, 而柯西应变不具有这个特性.若用形变梯度表示, 格林应变表示为${\varepsilon _G} = \frac{1}{2}({F^T}F - I).$

经常用来模拟非线性的材质模型StVK材质模型^[10,11]就是一种几何非线性模型.把公式(1) 中的柯西应力和柯西应变分别换成第二Piola-Kirchhoff应力和具有旋转不变性的二次格林应变张量, 就变成了StVK材质模型, 其应变张量和位移之间是非线性关系, 具有几何非线性.但应力张量和应变张量是线性关系, 因此从材质上还是线性模型, 其能量密度函数表示为$\mathit{\Psi }(F) = \mu {\varepsilon _G}:{\varepsilon _G} + \frac{\lambda }{2}t{r^2}({\varepsilon _G}).$

还有另外一类常用的非线性模型, 如NeoHookean材质模型, 这种材质的应力张量和应变张量是非线性关系, 其能量密度表示为$\mathit{\Psi }(F) = \frac{\mu }{2}(trC - 3) - \mu \ln J + \frac{\lambda }{2}{(\ln J)^2}, $其中, C=F^TF, J=detF.

除上述模型外, 图形学中还有一些常用的非线性材质模型, 如Mooney-Rivlin, Ogden, Yeohmo模型等^[12].

基于四面体离散化方法的有限元仿真方法中, 当发生大形变或者仿真选取比较大的时间步长时, 四面体可能出现翻转的情况.即使针对不可压缩材料, 采用离散化的方式表示连续的物质, 由于离散化的误差也不能完全保证四面体单元是保体积的, 因此也会出现翻转的情况.

早期在质点弹簧系统中, 利用高度弹簧(altitude spring)解决翻转的问题^[13].这种简单的处理翻转的方法也同样适用于有限元方法, 但采用这种方式会改变有限元仿真的微分方程, 并且失去传统有限元方法的许多优势.研究人员也利用高度弹簧的方法做了相似的研究, 如:Teschner等人^[14]通过添加保体积项作为惩罚项避免翻转的情况; Müller等人^[4]移除旋转项, 利用线性模型在无旋转的空间中仿真, 可以很好地避免单元翻转, 但这种方法在材质属性上有很大的局限性; Irving等人^[15]提出了可翻转的有限元框架, 能够检测出有限元仿真过程中四面体单元是否出现了翻转, 该方法的最大优点是适用于任意的弹性材质的本构模型, 并在其随后的研究中^[16]将该方法推广到基于六面体单元离散化方法的有限元仿真中; Teran等人^[17]利用隐式积分, 在有限元仿真框架中提出一种改进的Newton-Raphson迭代算法, 更好地保证形变至翻转状态时的稳定性.

以上方法都是利用形变梯度奇异值分解算法, 对形变梯度进行对角化检测单元是否翻转, 并通过截断或者线性外插的方式计算翻转状态下第一Piola-Kirchhoff应力.而Stomakhin等人^[18]提出基于能量的方法, 定义了翻转状态的形变能量, 进一步提高方法处理翻转问题的鲁棒性, 图 1表达了被固定手脚模型的单元翻转与恢复.

1.4 各向同性材质模型

各向同性材质(isotropic material)表示材质在各个方向上抵抗外力产生的形变是一致的.若一个超弹性本构关系是各向同性的, 那么其应变能量密度函数应满足Ψ(FQ)=Ψ(F).其中, Q表示任意3×3的旋转矩阵.上述的共旋线性弹性模型、StVK模型以及NeoHookean模型都是各向同性材质模型^[2].

对于各向同性的材料, 为避免对形变梯度进行SVD分解的复杂操作, 利用形变梯度的3个各向同性的不变量(isotropic invariants)表示能量.3个不变量利用形变梯度表示为

$ {I_1}\left( F \right) = tr\left( {{F^T}F} \right), {I_2}\left( F \right) = tr\left[{{{\left( {{F^T}F} \right)}^2}} \right], {I_3}\left( F \right) = {\rm{det}}\left( {{F^T}F} \right) = {\left( {{\rm{det}}F} \right)^2} $

(6)

利用各向同性本构模型的不变量表示的NeoHookean本构模型的弹性能量表示为

$ \mathit{\Psi }({I_1}, {I_3}) = \frac{\mu }{2}({I_1} - \log ({I_3}) - 3) + \frac{\lambda }{8}{\log ^2}({I_3}) $

(7)

1.5 各向异性材质模型

与各向同性材料相比, 各向异性材料是更为复杂的材质模型, 对各向异性材料的设计和应用有很多研究工作^[19].如布料仿真领域, 由于布料纤维的构成, 模拟各向异性的行为对于布料运动的真实度具有非常重要的意义.Etzmuβ等人^[20]利用正交异性有限元方法模拟布料运动, Huber等人^[21]、Garg等人^[22]利用正交异性模型分别模拟了湿衣服和布料的弯曲效果.Thomaszewski等人^[23]利用约束的方式限制柯西应变张量每个方向上的独立分量模拟各向异性材料, Wang等人^[24]提出了分段线性的各向异性的布料模型, Allard等人^[25]利用二维各向异性材料模拟了轻薄软布料的撕裂效果.对各向异性材料的研究主要集中于横向各向同性材料(transversely isotropic)、正交材料(orthotropic material)以及一般化的各向异性材料设计, 目前, 在图形学领域侧重于前两种各向异性材料的研究.

1.5.1 正交材料

在工程领域, 正交材料的9个独立的正交参数都是通过实际的物质测量获得的.在图形学领域, Li和Barbič等人^[26]首次提出了一种比较直观的方式调整正交材料参数.

利用第1.1节的线性弹性模型, 如公式(1)~公式(3) 中, 将张量表示为向量形式, 分量11, 22, 33是正分量(或法向分量normal components), 12, 23, 31是剪切分量(shear components)(为前后文一致, 这里的弹性张量依然用C表示, 而没有用文献[26]中表示弹性张量符号).

与各向同性的材料中只有2个参数的杨氏模量和泊松比不同, 正交材料有3个不同的杨氏模量E₁, E₂, E₃, 每个量代表一个正交方向, 6个泊松比ν_ij, i≠j.杨氏模量E_i表示正交方向i的刚性; 泊松比ν_ij表示当第i个方向拉伸时, 第j个方向上的收缩系数.对于一般的各向异性材料, 在正交方向和剪切方向上的应变张量都受正交和剪切方向的应力影响.因此, 弹性张量是个稠密矩阵, 对于正交材料, 正方向和剪切方向解耦合, 正方向的应力只影响正方向的应变, 剪切方向的应力只影响剪切方向的应变, 并且剪切方向的12, 23, 31方向的应变通过一个剪切模量(shear modulus)μ_ij受相应方向的剪切应力影响.根据正交材料的假设, 弹性张量C有9个独立分量, 表示为分块对角形式.

$ {C_{ortho}} = \left[{\begin{array}{*{20}{c}} {\frac{1}{{{E_1}}}}&{-\frac{{{\nu _{21}}}}{{{E_2}}}}&{-\frac{{{\nu _{31}}}}{{{E_3}}}}&0&0&0\\ {-\frac{{{\nu _{12}}}}{{{E_1}}}}&{\frac{1}{{{E_2}}}}&{ - \frac{{{\nu _{32}}}}{{{E_3}}}}&0&0&0\\ { - \frac{{{\nu _{13}}}}{{{E_1}}}}&{ - \frac{{{\nu _{23}}}}{{{E_2}}}}&{\frac{1}{{{E_3}}}}&0&0&0\\ 0&0&0&{\frac{1}{{{\mu _{12}}}}}&0&0\\ 0&0&0&0&{\frac{1}{{{\mu _{23}}}}}&0\\ 0&0&0&0&0&{\frac{1}{{{\mu _{31}}}}} \end{array}} \right] $

(8)

Li等人通过一定的约束和假设条件, 提供了9个独立分量的计算方式.在随后的研究中, 他们把对正交材料调整参数的方法推广到更加一般化的各向异性材料^[27].

1.5.2 横向各向同性材料

横向各向同性材料实质上是一种特殊的正交材料, 当3个正交方向中的两个刚性系数相同时, 即为横向各向同性材料.对于这种材料, 材料的一个平面方向是各向同性的, 这种材料的弹性张量有5个独立分量.

横向各向同性材料最早由Bonet等人^[28]提出, 在线性和非线性有限元方法中应用比较广泛.在医学仿真领域, Picinbono等人^[29]将线性的横向同性材料用于虚拟手术仿真领域, 达到了交互级的时间效率; 之后, 又把材料扩展到非线性的横向同性材料^[30].Liao等人^[31]利用横向和正交材料从CT数据中获取骨骼材料.Sermesant等人^[32,33]和Talbot等人^[34]采用横向各向同性材料建立心脏的机电模型.Comas等人^[35]在GPU上构建了横向各向同性的超弹性粘性模型.Teran等人^[36]和Sifakis等人^[37]利用横向各向同性材料分别建立了人体肌肉组织模型.

在调整各向异性材料的仿真方法中, 有一类称为应变约束方法, 其最早应用于基于质点弹簧系统模型的布料仿真^[38,39].Desbrun等人^[40]改进了基于质点弹簧模型的布料交互仿真过程中过拉伸的问题.Müller^[41]利用一种多分辨率的方法, 在质点弹簧系统上施加了应变约束.之后, 由Thomaszewski等人^[42]推广到有限元方法中.由于三角单元的物质形变行为是由物质的刚度属性决定的, 可以利用简单的迭代方法投影避免过度形变, 投影的过程中, 利用每个方向上应变约束调节结点的速度和位置值, 计算正确的速度向量.

利用有限元方法并以四面体网格作为形变体的离散化方式的计算应变约束方法中, 针对线性的材质模型, 应力张量在每个四面体单元内是常量.根据Hernandez等人^[43]的方法, 沿着每个四面体的形变梯度的主方向限制应变值的大小, 首先将每个四面体的形变梯度进行SVD分解:

$ F = U\Sigma {V^T}, \Sigma = \left( {\begin{array}{*{20}{c}} {{\lambda _1}}&0&0\\ 0&{{\lambda _2}}&0\\ 0&0&{{\lambda _3}} \end{array}} \right) $

(9)

其中, U, V分别表示旋转正交矩阵, 奇异值{λ₁, λ₂, λ₃}表示沿着3个主轴方向的拉伸或压缩程度.当奇异值都为1时, 表示没有形变.应变约束的方法就是给定一个压缩约束λ_min, 一个拉伸约束λ_max, 使得形变梯度的每个奇异值都在约束范围内, 即λ_min≤λ_i≤λ_max.

针对各向同性的形变材料, Wang等人^[44]提出了多种分辨率下, 各向同性的应变施加约束的方法, 通过在每个形变梯度的主轴方向设置同样的应变约束值, 直观地获得控制材料应变约束的方式.

相对于各向同性的材料, 各向异性的材料模型的应变约束复杂得多.Picinbono等人^[45]通过添加一个能量项模拟超弹性形变, 解决了横向各向异性约束问题.该方法只限制横向拉伸方向的值, 其他方向自由形变, 采用这种方式实现了一个简单的各向异性行为, 由于只约束一个方向, 因此该方法只是针对应变张量沿着横向轴的一个投影, 并不是真正意义上的各向异性.若要模拟更一般化的各向异性材料, 各个方向上的约束需要取不同的值.此外, 由于拉伸和压缩约束都是定义在全局的参考坐标系下, 形变方向定义为沿着形变梯度SVD分解的主轴方向, 而各个单元上的标架是不同的.通常情况下, 两种坐标系是不匹配的, 因此需要计算拉伸和压缩约束沿着形变梯度的主轴方向的约束值, 这种投影操作非常复杂.针对这个问题, Hernandez等人^[43]提出了沿着全局坐标系定义形变方向的计算方法, 实现了应变约束在各向异性材料中的仿真模型, 成功地模拟了手指等物体的各向异性材料的形变行为.

通过上述材质的分析及现有成果, 可以总结出一些均质材质的特点和主要研究方向.

(1) 线性模型非常简单, 但不能模拟复杂形变.在此基础上的共旋的线性模型虽然一定程度上改善了旋转问题, 但对于复杂的非线性形变不能有效地模拟.以StVK材质为代表的几何非线性模型虽然不完全是非线性模型, 但在图形学的仿真模拟中仍然应用广泛.尤其是随着Barbič等人^[46]提出的降维仿真方法发现StVK材质的内力可以表现为位移的三次多项式, 并在预计算阶段求解其系数, 使得仿真效率获得大幅度的提高.2008年, An等人^[47]扩展了Barbič等人在降维仿真中针对StVK材质的内力表示形式, 利用优化的积分单元替代全局的单元计算内力, 提高仿真效率的同时, 使得降维仿真方法不局限于StVK材质, 具有非常重要的意义.从仿真趋势来看, 对于均质材质, 利用NeoHookean等人非线性材质的仿真方法更能说明方法的通用性;

(2) 单一的各向同性材质的模拟从效果上具有一定的局限性, 模拟各向异性材料, 并且是各个方向上形变不同的一般化的各向异性的材质具有一定的难度, 但同时也是未来的研究趋势.

2 非均质材质

前文中提到的各种材质是从微观描述材料属性, 但在现实生活中, 大多数材质表现的是宏观属性, 如由多种物质构成的复合材料, 即, 非均质材质.相对于均质材质, 非均质材质在CAD、CAE和CAM等领域具有更大的研究价值, 在几何建模方法以及复杂的材质属性分配策略方面已经取得了很多研究成果^[48,49], 模拟非均质材质的形变体仿真具有重要的意义.在基于物理的仿真领域, 非均质材质的各种仿真方法有各自的优势和难点, 但总结各项研究策略, 非均质材质的仿真大体遵循以下思路.

(1) 均质化方法.高效的模拟复杂材质的形变行为是基于物理的仿真方法追求的目标, 为提高仿真效率, 现有的研究包括利用内嵌网格的方式、降维的方法、基于几何的形状匹配的方法等.利用降维仿真或者基于几何的全局形状匹配的方式, 虽然能够提高效率, 但无法捕捉局部的细节信息, 因此在模拟非均质材质的仿真时, 比较重要的思路是:采用低分辨率的模型代替高分辨率模型提高仿真速度的同时, 尽可能保证高分辨率的非均质材质的形变特征;

(2) 材质细分方法.对于一个给定的非均质的仿真模型, 要表现模型宏观的材质特征, 模拟计算的思路是:尽可能将不同材质进行分离, 使得计算单元内部只包含单一的均质材质, 或者在计算过程中能够体现材质的属性等方法;

(3) 数据驱动方法.第3种研究方法基于数据驱动的方式, 通过给定的一组数据或目标形变, 利用优化求解的方式获得模型的材质属性.

2.1 均质化方法

在模型离散化过程中, 利用高分辨率的模型模拟非均质材质更能体现材质的形变细节, 但仿真效率低, 利用粗糙的模型代替高分辨率模型提高仿真效率并保持材质特征, 是研究的重点.针对有限元线弹性模型, Nesme等人^[50]将模型用六面体体素化形式表示, 利用模型不断地合并获得低分辨率的模型.通过对弹性张量求平均的方法, 近似表示计算单元内非均一的刚度分布, 获得了任意分辨率网格单元的材质属性; 同时, 鲁棒地处理了退化情况.但该方法采用简单的求平均策略获得粗糙网格内部的质量矩阵和刚度矩阵, 无法准确描述材料分布.

均质化理论^[51-53]能够利用宏观的材质表现微观的细节, 2009年, Kharevych等人^[54]首次在图形学领域借鉴均质化理论, 用非均质的粗糙模型替代精细仿真模型提高计算效率.其核心思想是:输入的精细分辨率的仿真物体的每个单元由均质材质构成, 每个单元之间可以具有不同的物质属性, 通过将粗糙单元与精细单元对应, 利用预先定义的平移和剪切方向的6组协调位移(如图 2所示)计算能量, 建立两个分辨率模型能量的等价关系, 从而计算出对应的粗糙分辨率的非均质材质属性.

具体地, 对于一个高分辨率网格D, 每个四面体T_p(p=1, …, |D|)单元上都具有不同的弹性张量${C_{{T_p}}}$, 对应的粗糙四面体单元个数|D| < < |D|, 需要计算出对应的每个粗糙单元的弹性张量${\mathbb{C}_{{\mathbb{T}_p}}}$, 建立两个分辨率之间的单元对应关系及能量等价关系, 分别表示为

$ (D = \{ {T_p}\}, C = \{ {C_{{T_p}}}\} )\mathop \to \limits^{粗糙化} (\mathbb{D} = \{ {\mathbb{T}_p}\}, \mathbb{C} = \{ {\mathbb{C}_{{\mathbb{T}_p}}}\} ) $

(10)

$ \int_{{T_q}} {\varepsilon (u):C:\varepsilon (u){\rm{d}}V} = \int_{{\mathbb{T}_q}} {\varepsilon (\mathbb{U}):\mathbb{C}:\varepsilon (\mathbb{U}){\rm{d}}V} $

(11)

求解的物理量是粗糙单元的弹性张量, 根据其对称性分析, 具有21项独立分量.这个方法提供了很好的构建低分辨率网格非均质材质的方式.但该方法利用的位移是比较简单的线性形变, 并且计算能量的应力应变关系是线性关系, 使得仿真方法具有一定的局限性.2014年, Torres等人^[55]扩展了Kharevych等人的方法, 在粗糙网格上模拟共旋的线性弹性模型, 提高了模拟形变的范围, 并将其应用于医学研究领域.

2.2 材质细分方法

在粗糙网格模拟非均质材质的形变体仿真方法中, 当同一网格内部存在多种材质时, 按照传统方式不易区分, 如何有效分离不同材质的求解过程, 是需要考虑的主要问题.2009年, Nesme等人^[56]针对粗糙的网格模拟内嵌在其中的复杂材质结构的仿真问题, 提出在预处理阶段, 利用高分辨率网格上特定边界条件的准静态分析方法, 得到粗糙网格内部的插值关系, 使粗糙网格单元内部的插值能够反映材质分布.该方法通过复制单元的方式, 有效区分了几何位置上处于同一网格但模型上未连接区域单元内部的不同分支结构; 并且, 他们还考虑了无材质单元的仿真问题.

借鉴于Kharevych等人^[54]的方法, Schumacher等人^[57]打破了3D打印依赖单一材质只能表示一种物质属性的限制.该方法实质上是均质化方法的逆向过程, 通过给定的粗糙的非均质材质的模型, 利用打印材料不同的微观结构单元, 近似表示给定粗糙模型的材质属性, 采用柯西应变表示的线性模型表示能量, 通过基础材质结构优化拟合目标材质, 获得了能够打印的非均质模型, 如图 3所示.

Zhao等人^[58]针对传统物质点法(material point method)模拟非均质材质时无法区分多种材质的问题, 提出了增加求解自由度细分材质的策略.传统物质点法利用固定不变的背景网格结点求解方程, 利用网格内部的物质点与网格结点传递物理量, 同一网格内部的粒子与背景计算网格插值仅考虑几何位置分布, 并不区分不同的物质属性, 因此很难模拟出同一网格内部的非均质材质的属性.他们提出引入粒子区域, 在预处理阶段检测构成物质的不同材质的边界区域, 仿真过程中, 在物质边界区域包含的粒子所在的背景网格中的所有粒子利用粒子区域求解, 增加了求解的自由度, 有效地利用物质点法模拟了非均质材质的弹性形变, 拓展了物质点法的适用范围.该方法的缺点是:在传统物质点法基础上引入粒子区域, 每次都需要更新粒子区域的位置关系, 并且需要从粒子到物质区域再到网格结点两层映射, 时间效率上有一定影响.

此外, 利用粗糙的网格模拟多种材质的仿真方法中, 还包括一类方法在粗糙网格内部区分材质, 思路是将材质属性作为传统的基于位置的形函数的一部分模拟非均质材质的仿真方法.典型的代表是Faure等人^[59]在无网格框架下, 提出与材质相关的形函数, 在构建插值形函数时考虑形函数作用域内材质的刚度属性, 使得形函数不仅与几何相关, 也与材质属性的分布相关.该方法可以模拟一些非均质材质的形变属性, 如图 4展示的T骨牛排的形变, 牛排包括一个比较硬的骨头以及比较软的肉, 其刚度属性分布如第2张图所示, 第3张图是其离散化表示, 第4张表示其提出的柔性(compliance)距离, 最后一张是其形变效果.从仿真结果来看:虽然该方法从一定程度上模拟了非均质材质的特征, 但由于新的形函数的材质属性只考虑了一个维度的刚度系数的标量, 因此对于全面模拟三维非均质材质的复杂形变有一定局限性.

2.3 数据驱动的方法

针对非线性材质模型, Bickel等人^[60]提出一种数据驱动的方式, 利用捕捉设备得到一组形变位移与力的关系数据, 表示成应力应变关系曲线.仿真过程中, 利用关系曲线之间的非线性插值, 得到仿真的外力对应的形变效果.该方法的缺点是外力构建形式比较简单, 得到的形变效果有限, 无法准确模拟复杂运动的形变.如图 5所示:首先, 利用力传感器设备从真实物体上获得几组样例形变的力的大小和方向; 然后, 利用这些数据得到应力应变关系曲线.仿真过程中, 利用非线性插值, 从采样得到的应变空间计算模型的材质.

从材质设计的角度来说, 要想设计一套符合运动规律及理想的非均质材质模型, 在材质设计过程中, 每次改变模型的材质属性, 都需要长时间的预处理过程才能获得对应的非均质粗糙单元的材质.如果所设计的形变不满足设计师要求, 这个处理过程需要不断重复, 非常耗时.针对这种设计缺陷, Chen等人^[61]提出利用规则六面体网格构建仿真物体, 针对每个小立方体的均质材质, 按照不同的材质组合采用8个立方体合成1个的方法, 分层级求解出不同分辨率模型的材质属性保存在数据库中.仿真过程中, 根据需求选取材质, 调节材质参数非常便利, 只需查询数据库中对应物质属性, 整个材质设计过程高效易行.采用规则立方体的方式也非常便利于精细分辨率到粗糙分辨率的对应与合成, 同时, 嵌入表面网格的方式可以获得非常精细的绘制结果.

综合现有的非均质材质的研究, 从方法和效果来看, 今后的研究思路可以总结为以下几点.

(1) 虽然提高网格分辨率, 利用更精细的模型仿真, 可以获得更加准确的非均质材质的仿真结果, 但仿真效率受到很大影响, 因此, 未来的趋势还是探讨如何利用粗糙的分辨率模型更准确地模拟非均质材质的形变效果;

(2) 在有限元仿真方法中, 大多利用四面体单元或者六面体单元进行离散化, 从模型的粗糙化角度来说, 规则的六面体单元相对于四面体单元具有天然的优势, 粗糙模型与精细模型之间可以获得更加精准的对应关系, 避免四面体模型在两种分辨率下边界无法对应的问题;

(3) 针对应力应变之间的关系, 虽然线性关系处理比较简单, 但其无法处理更加复杂的形变行为, 因此, 非线性的形变关系是发展的必然趋势.

3 材质编辑

随着3D打印技术、VR技术、计算机动画的发展, 传统的材质模型已经不满足需求, 在实际应用中, 设计师往往需要更复杂的材质模型描述形变.因此, 在图形学领域, 出现了关于编辑和设计形变体材质的研究.其中, 一类针对材质属性直接调整参数的方式设计新材质; 另外一类不调整材质参数, 而是通过给定一个仿真目标或嵌入式结构改变形变方式, 展现丰富的动画效果.

3.1 基于材质编辑的方法

最初的材质属性是通过测量真实物体的物质属性获得的, 主要参数包括杨氏模量或泊松比^[60,62,63]、粘弹性属性^[64,65]、塑性参数^[66]、布料的拉伸和弯曲系数^[67]等.之后的研究主要通过添加一些外部条件获得物质参数, 将仿真物体离散化表示, 在顶点上施加力, 通过顶点的形变位移计算材质属性, 或者在仿真过程中同时优化物体的拓扑变化和物质参数获得模型的材质信息.

2015年, Xu等人^[68]受工程上研究不可压缩材料拟合方法的启发^[69], 提出了利用主拉伸方向的非线性可压缩材料的设计方法.方法思路是:将内力和刚度矩阵表示成主拉伸方向的几个物理量的函数形式, 通过编辑应力应变曲线设计新的各向同性的材料; 对于各向异性材料的设计, 除了对各向同性参数曲线的修改, 还添加材料的3个正交的拉伸方向的曲线调整.该方法基于有限元仿真方法, 将能量表示为Ψ(λ₁, λ₂, λ₃)形式, 其中, λ₁, λ₂, λ₃为形变梯度经过SVD分解的3个主拉伸方向上的值.应变能量表示成不变量的形式为

$ {{\rm{I}}_C} = trace(C) = \lambda _1^2 + \lambda _2^2 + \lambda _3^2, {\rm{I}}{{\rm{I}}_C} = \det (C) = \lambda _1^2\lambda _2^2\lambda _3^2, {\rm{II}}{{\rm{I}}_C} = C:C = \lambda _1^4 + \lambda _2^4 + \lambda _3^4 $

(12)

能量密度函数表示为

$ {\mathit{\Psi }}({\lambda _1}, {\lambda _2}, {\lambda _3}) = f({\lambda _1}) + f({\lambda _2}) + f({\lambda _3}) + g({\lambda _1}{\lambda _2}) + g({\lambda _2}{\lambda _3}) + g({\lambda _3}{\lambda _1}) + h({\lambda _1}{\lambda _2}{\lambda _3}) $

(13)

其中, f, g, h分别为非线性应变能量密度函数的单轴(长度)、双轴(面积)、三轴(体积)形式.通过将能量密度函数对主拉伸方向的第一个分量求导, 得到如下表示形式:

$ \frac{{\partial \mathit{\Psi }}}{{\partial {\lambda _1}}} = f'({\lambda _1}) + g'({\lambda _1}{\lambda _2}){\lambda _2} + g'({\lambda _3}{\lambda _1}){\lambda _3} + h'({\lambda _1}{\lambda _2}{\lambda _3}){\lambda _2}{\lambda _3} $

(14)

其他分量的导数类似求得, 通过在一定的约束下调整f', g', h'的曲线, 可以获得能量密度的唯一计算方式.对于正交的非线性材料的设计, 定义的应变能量密度函数表示为各向同性和正交项两部分.

$ \mathit{\Psi } = {\mathit{\Psi }_{iso}}({\lambda _1}, {\lambda _2}, {\lambda _3}) + {\mathit{\Psi }_{ortho}}, {\Psi _{ortho}} = {\omega _1}({\bar \lambda _1}) + {\omega _2}({\bar \lambda _2}) + {\omega _3}({\bar \lambda _3}) $

(15)

其中, ${\omega _i}({\bar \lambda _i})$表示材料的3个正交的拉伸方向的能量项.对于这种非线性材质, 除了调整前文所述的各向同性的曲线之外, 还需调整3个正交项曲线的参数, 以达到各向异性的非线性模型的形变效果.他们的方法提供了一种非常直观灵活地调节材料属性的方式, 这是传统非线性材质不具备的优势, 但该方法并没有给出保证材质形变稳定性的参数调整范围.通过参数调整, 获得不同的仿真结果的例子如图 6所示, 球模型由于调整了参数曲线, 获得了不同的压缩和拉伸效果.

前文所述的通过施加力和位移的方式获得物质参数的方式, 忽略了一个问题, 就是所施加的力并不是守恒的.2016年, Miguel等人^[70]针对这个问题提出添加一个能量约束模拟真实世界的非线性的形变, 保证了弹性力在创建过程中是恒定的.

随着3D打印机的诞生以及打印材料类型的增多, 人们希望能够打印出更符合实际需求的具有复杂外观、表面光学特征及力学特性的物体.这一需求催生了3D打印中的一个重要问题:如何确定一个物体对象的材料组成, 使其满足给定的表面外观效果或变形的功能要求^[71].这个问题涉及到结构设计优化^[72]、结构的匹配连接以及整体材质的合成等.现有的研究针对3D打印中模型的设计已经取得一定的成果, 如针对指定形状的模型的平衡问题^[73]、关节结构的角色设计^[74]以及通过给定的形变形状设计初始形状^[75,76].

针对3D打印技术的材料设计问题, 利用被称为逆均质化(inverse homogenization)的方法将多个基础的微结构组合成符合模型整体形变特征的过程.针对基础微结构的设计已经取得一些进展, 包括弯曲的板型结构、球状的薄壳体、刚体单元等.利用这些基础微结构的拓扑优化问题, 为材质的设计提供了更加多样化的选择空间, 但同时也引入了如计算开销过大等问题.Rodrigues等人^[77]和Coelho等人^[78]提出利用分层的拓扑优化方法计算粗糙的材质对应的精细微结构的材质属性, 虽然每个微观的结构单元可以独立优化, 但其优化方法非常耗时, 并且无法保证邻域结构之间的连接关系.针对这个问题, Zhou等人^[79]提出了保证微结构边界连接关系的方法, 提高了匹配效率.

单一打印材料的形变行为是有限的, 如何从基础材料设计出更加符合实际需求的复合材料, 是非常有意义的工作, Bickel等人^[80]提供了一种设计和组装新材料的方法, 首先获得几种基础材料的形变数据.该方法整个优化求解新材料的过程是寻找最符合用户指定样例的形变方式, 通过给定的基本的均质材料及其力学性能曲线, 将基本材料混合, 获得符合指定力反馈的混合材料, 其流程如图 7所示.采用的按层方式混合不同的材料, 并引入一个优化过程来获得最接近指定材料性能的混合结果.

从均质化的角度讲, 他们的方法是对Kharevych等人^[54]非均质材质模拟工作的线性应力应变模型到非线性材质的扩展, 不同的是, Bickel等人工作的基础材料, 即, 高分辨率的精细模型的形变是已知的, 如何生成一个能够满足目标的粗糙模型的力和位移关系的非均质材质是需要解决的问题.他们方法的优势是生成的复合材料是可以真实打印的, 具有很好的实用性, 缺点是调整参数的方式比较复杂, 并且当基础材料超过5层的时候, 时间效率比较低.为了提高效率, Xu等人^[81]提出了利用降维空间, 采用交互式的方式设计物质材料, 提高仿真速度的同时, 提供了一种直观地调整物质参数的方式.其核心思想是:通过对部分网格顶点指定位移和内力, 利用有限元方法快速求解一个优化问题计算各向异性的物质属性的空间分布.

Schumacher等人^[57]利用均质化理论, 为3D打印技术提供了一种通过单一的打印材质将虚拟的非均质材质模型变成可打印的真实的物体的方法.他们借助数值优化的方法设计微观单元的结构, 模拟宏观上非均质的材质属性.其核心思想是:在预处理过程中建立一个超材料的数据库, 将材质参数与微观结构作一对一的映射.如图 8所示, 每个各向同性的物质参数(杨氏模量和泊松比)构建的坐标系上的点对应一种微观结构.通过插值优化得到的不同微观结构的物质属性可能重叠, 为了获得指定的小单元的物质属性并考虑邻域结构之间的连接关系, 采用全局优化的策略选择最合适的微观结构.

相对于Bickel等人的研究需要指定一系列的基础材料, 并逐层按照目标形状拟合的方法, Schumacher等人的方法无需指定基础材料, 自动计算精细单元的结构, 并且从效率上可以计算数千层(或单元), 只需要秒级别的时间, 计算效率远远高于前者.同年, 针对3D打印技术中物体的材质问题, Panetta等人^[82]也提出了一种在空间上设计微结构以满足不同材质属性的方法, 该方法通过调整基本材质单元边的厚度, 并为顶点位置添加补偿量的方式调整基础结构获得新的微结构.在设置打印阶段, 根据要打印模型的每个单元的杨氏模量和泊松比, 查找一个对应的填充形状, 并获得使邻域互相连接的微观结构.该方法不仅建立各向异性的打印材料, 还能够模拟指定的形变行为.该方法的缺陷是采用了线性的材质模型, 因此生成的材料对于某些非线性的形变无法满足.

针对多材料打印问题, Vidimče等人^[83]给出了一个可编程流水线OpenFab解决多种材料打印的合成问题. Chen等人^[84]提出了减速调谐器(reduced-tuner)优化调整调谐节点, 调谐节点构成了一个调谐网络(tuner network), 调谐节点利用优化管线与调谐网络优化其相关联的简化节点参数, 实现系统对材质和几何形状的优化目标.2013年, Skouras等人^[85]提出了一种利用软硬两种材质构建可打印材质的方法, 通过虚拟的形变物体以及一组形变姿势作为输入, 获得一个与输入形变和运动逼近的包含内部物质属性分配的可打印的模型.

3.2 基于形变编辑的方法

前文所述的传统仿真方法中, 一旦仿真物体的材料模型确定, 它的形变行为也随即确定.即便通过调整应力应变曲线能够获得更加丰富的材质属性, 但其形变方式符合应变曲线规律, 获得的动画效果有限.事实上, 除了通过直观地调节物质参数或优化微观结构改变形变体的运动方式外, 还有一类设计方法通过调整物质的形变行为展现更加丰富的动画效果, 如基于样例材料的仿真.从仿真物体角度考虑, 基于样例的材料仿真包括针对三维固体形变体、薄壳体^[86]、刚体, 从仿真现象包含形变体的弹性形变、塑性形变以及刚体的塑性^[87]、碎裂^[88]等多种效果, 这里, 我们主要讨论形变体的弹性形变.

2011年, Martin等人^[89]提出了一种基于样例材料的方法模拟复杂弹性材料的形变行为, 通过给定一组仿真物体的形变姿势, 利用形状插值的方式获得一个形变空间.仿真过程中, 样例流形(example manifold)作为一个额外的弹性力, 引导仿真物体向着给定样例姿势的方向形变, 从而达到设计各向异性和非均质材质形变体仿真的目的.方法的核心思想是:在物体的运动方程中添加一项额外的弹性势能, 将物体引导至由给定的一些样例组成的形变空间.基于样例材料的有限元仿真基本方程为

$ M\ddot x + {\nabla _x}{W_c}(x) + {\nabla _x}{W_p}({x_\omega }, w, x) = {f_{ext}} $

(16)

其中, M表示仿真物体的质量矩阵, $x, \ddot x$分别表示位移和加速度, $\nabla $ 表示微分算子, W_c, W_p表示传统的形变能量和基于样例的能量, w表示样例之间的插值权重, x_ω表示样例空间上获得的插值形状, f_ext表示系统的外力.

该方法利用格林应变作为样例描述符, 假设给定的样例由m个四面体组成, 则该样例的描述向量为6m维, 每个四面体的应变是定值, 表示为E=[E₁, E₂, …, E_m]^T∈R^6m.给定的样例组成的样例空间描述了目标的形变, 样例之间线性插值得到应变状态E_ω.但直接线性插值得到的应变状态E_ω可能是不符合真实的, 即:不存在物体的一个状态x_ω, 使得每个单元的应变都等于插值得到的值.因此, 作者利用非线性的优化问题, 找到距离插值得到的形变最接近且真实的形变状态, 优化公式为

$ \mathop {\min }\limits_{{x_\omega }} {W_I}({x_\omega }, \omega ) = \mathop {\min }\limits_{{x_\omega }} \frac{1}{2}|E({x_\omega }) - E(\omega )|_F^2 $

(17)

在仿真过程中, 通过定义弹性势能, 在样例空间中找到符合真实且距离物体当前状态弹性势能最小的一个状态, 将物体从当前状态投影到样例空间中.假设当前形变状态x, 目标是样例空间中找到的x_ω, 即:

$ \mathop {\min }\limits_{{x_\omega }, \omega } W({x_\omega }, x)\;\;\;{\rm{s}}{\rm{.t}}{\rm{.}}\;\;{\nabla _{{x_\omega }}}{W_I}({x_\omega }, \omega ) = 0 $

(18)

其中, 约束项是为了满足得到的x_ω是真实的.得到x_ω后, W(x_ω,x)被作为额外的一个弹性势能项加入到物体的运动方程中, 引导物体的形变.如图 9所示, 汽车模型在不同样例作用下获得的仿真结果.

这种基于样例的方法能够得到与给定的目标样例很接近的弹性形变, 而且物体在形变过程中运动比较真实自然.但由于在样例空间中寻找符合真实的插值形状时需要求解一个非线性的优化问题, 因此该方法的计算开销比较大.效率问题是基于样例材料仿真方法面临的一个主要问题, 为缓解这个问题, Schumacher等人^[90]提出了一种高效的基于样例的形变仿真方法, 能够得到类似的效果, 且速度加速数十倍.该方法对于给定的样例以及样例空间中插值得到的中间状态都使用不共用顶点的单元表示, 直接线性插值单元的顶点位置, 而不像Martin等人先得到通过应变线性插值的中间状态, 再求解非线性优化问题得到符合真实的物体形状.

为提高基于样例材料的仿真效率, Koyama等人^[91]从几何角度出发, 利用无网格的形状匹配方法, 使得整个仿真方法达到了实时.由于他们用的是形状匹配方式, 且将样例的线性插值直接作为样例空间, 因此能达到非常快的仿真效率, 但效果不如Martin等人的方法.2015年, Zhang等人^[92]利用有限元方法提出在降维子空间中求解基于样例材料的实时仿真方法, 效率提高的同时, 避免了几何方法在仿真过程中一些物理真实性不足的问题.他们通过定义一种基于样例的格林应变张量计算势能, 通过这个势能将形变引导到样例空间, 样例权重可以显式获得, 内力变为传统有限元计算内力部分和由样例材料引导部分.利用降维子空间, 将两部分内力都表示成位移多项式的形式, 达到了实时的仿真效率.但该方法的缺陷是计算内力的方式只适用于StVK材质, 不具有通用性.

传统的基于样例的仿真方法, 样例与形变物体通常采用同样的拓扑结构, 这使得形变材料的设计不够灵活.而突破拓扑结构的问题, 最大的挑战在于找到一种描述方式, 使不同的模型在同一个空间下表示.为解决这个问题, Zhu等人^[93]提出利用Laplace-Beltrami形状空间作为形状插值空间, 利用仿真物体与样例之间的等距不变性特征, 使仿真物体与样例具有不同的拓扑关系以及网格分辨率, 并且由于构建的形状空间是一个降维子空间, 求解插值目标形状时具有更高的效率.他们利用Laplace-Beltrami算子为不同形状定义了一个放缩无关、拓扑无关的形状空间, 输入的样例构成了该形状空间中的一个流形作为形状插值空间.通过高效地求解一个降维的最小二乘优化问题, 得到一个目标形状将仿真物体引导至输入样例.如图 10所示, 利用Laplace-Beltrami形状空间, 长条形状的物体构建形变体仿真算法运行时的示意图.仿真物体与样例1采用了相同的模型, 不同的样例之间具有不同的分辨率与拓扑结构, 样例在各自的特征函数上投影获得一个样例流形, 仿真物体经过形变之后得到一个当前构型在欧氏空间中的几何描述, 仿真中, 通过将当前形状投影到样例流形上, 在样例空间中寻找与当前形状最接近的形状, 得到目标形状之后, 在欧氏空间下重建得到新的形变状态.

对样例流形上的形状, 定义其相对于所有样例的形变能量的加权之和, 用于指导目标形状t是最小化该加权能量和的对应形状, 要求解的问题描述为

$ \mathop {\min }\limits_t \sum\limits_{i = 1}^k {{\omega _i}E(t, {e_i})} $

(19)

其中, E(·, ·)表示两个形状之间的非线性形变能量, ω_i(1≤i≤k)衡量每个样例的指导权重.通过上式定义的目标形状在以形变能量为距离度量的标准下, 距离所有输入样例最近.在仿真过程中, 根据仿真物体的当前形状c与每个样例的接近程度, 动态决定每个样例的控制强度.

$ \mathop {\min }\limits_{{\omega _i}} \frac{1}{2}\left\| {\sum\limits_{i = 1}^k {{\omega _i}{e_i} - c} } \right\|_F^2\;\;{\rm{s}}{\rm{.t}}{\rm{.}}\;\;0 \le {\omega _i} \le 1, \sum\nolimits_{i = 1}^k {{\omega _i}} = 1 $

(20)

如图 11所示, 猩猩模型在armadillo样例的引导下与armadillo模型具有相似的运动行为.该方法的缺点是:在效率上, 虽然是在降维的Laplace-Beltrami形状空间上求解, 但计算能量的方式是在欧氏空间中利用表面网格顶点计算的形变能量, 虽然迭代只需要几个时间步, 但达不到实时的效率.另外, 作为引导的样例与仿真物体需要具有等距不变性, 特征函数个数的选择及仿真物体与样例特征函数的配准结果对仿真效果有一定影响.

除了基于样例的仿真, 还有类似的方法也达到了设计形变体形变行为的目的.Tan等人^[94]在物体内部嵌入肌肉纤维, 根据运动控制理论优化出特定运动目标对应的纤维长度, 然后将纤维的收缩力作用到物体上, 模拟软体动物的运动.Coros等人^[95]提出在仿真过程中根据特定的运动目标, 通过修改物体的初始状态获得内能, 使得物体在不借助外力的作用下, 运动到目标状态.Liu等人^[96]提出了一种在固体内部嵌入纤维, 通过调整嵌入纤维的形状控制固体形变的方法, 得到了与Martin等人^[89]类似的形变效果.该方法在无网格框架下, 定义一组仿真物体的离散化的物质粒子以及由植入纤维结构离散化表示的结构粒子, 物质粒子表示均质的各向同性的物质属性, 植入纤维的结构粒子按照线性绳索的能量方式表示拉伸、弯曲等形变模式来影响结构粒子周围粒子的运动方式, 从而体现出整个材质的非线性形变, 达到设计材质的目的.

时空约束方法也是一种比较经典的控制物体运动的方法, 该方法将物体的整个运动过程转化为一个受约束的优化问题, 物体运动遵循的物理规律及设计者给定的目标作为与时间空间相关的约束, 在最小化特定的目标函数后, 获得物体在整个时间段的运动状态.该方法能够获得既符合物理真实又符合设计要求的材质运动效果, 但由于求解的是一个时间空间都相关的优化问题, 因此计算开销大.Barbič等人^[97]提出一种在降维空间求解时空约束的方法, 根据给定的关键帧, 得到符合设计的物理仿真过程, 且控制物体运动所施加的控制力尽量小.虽然他们的方法在降维空间求解时空约束的问题, 减小了计算量, 使得仿真时间从小时级别降到分钟级别, 但仍然无法做到实时的效率.Hildebrandt等人^[98]提出了一种用时空约束控制物体形变的方法, 达到交互级别的速率, 该方法同样也在降维空间求解问题, 但与Barbič不同的是, 作者通过模态分析得到降维的基底, 在降维空间中, 问题被转化为多个一维的时空优化问题, 获得了快速求解问题的方法.

综上所述, 无论是编辑材质还是编辑形变行为, 本质上都是解决单纯依靠给定模型材质参数的调整无法模拟的形变效果, 丰富了形变体仿真的思路, 为虚拟仿真及动画模拟提供了一种新的方向.从现有的方法来看:

(1) 所有的编辑工作大体思路都是以传统材质模型为基础的一定程度的扩展和修正, 研究方向都是从各向同性的材质编辑到更加复杂和一般化的各向异性的材质编辑.从仿真方法而言, 现有研究采用有限元方法较多, 也包含一些无网格方法、形状匹配方法等.因此, 扩展仿真方法所适用的材质范围是可以研究的方向;

(2) 针对编辑形变方法的基于样例形变材料仿真中, 大多数的研究都包含了针对不同局部区域的形变方式, 使得同一个模型可以分区域按照不同的样例引导形变, 既有利于算法并行效率的提升, 又可以丰富形变体形变效果, 也是对非均质材质形变效果的扩展;

(3) 对材质形变效率的研究也是编辑工作的研究方向之一, 无论是根据应力应变曲线还是形变行为的引导, 都希望能获得交互级别甚至是实时的仿真效率.目前的研究方法中, 也针对效率的提升做出了一些探索, 多采用降维的方式或者降低网格分辨率的方式, 两种方法各有优缺点:降维仿真方法的形变基底在仿真过程中是不变的, 不利于模拟拓扑结构发生变化的碎裂等仿真现象, 并且降维方法也会丢失细节; 而降低分辨率的方式模拟材质一定程度上需要材质的近似计算, 模拟结果不够准确.因此, 从效率出发, 还有很多研究工作值得探讨.

4 比较分析

针对前文所述不同的材质模型和仿真方法, 表 1中列出了目前在图形学中比较重要的研究工作的仿真参数的对比, 包括三维模型离散化表示方式、仿真方法、积分求解器选择、材质模型以及仿真效率和模拟的仿真物体的适用范围.从表中可以看出:形变体仿真中, 由于不涉及到模型的拓扑结构的变化, 因此常用的方法是有限元方法在利用有限元方法的离散几何表示上, 一般采用四面体单元或者六面体单元; 从稳定性的角度, 利用隐式积分更好; 从效率上来说, 部分仿真方法已经可以达到实时的效率; 针对仿真对象, 很多研究方法不局限于三维固体, 同时兼顾了二维薄壳体以及布料等.各种方法针对研究的侧重点各有优缺点.

针对表中所列出的各类研究方法, 我们总结出关于材质发展和研究领域的一个整体脉络.

(1) 针对表中未列出的基本材质模型, 最初采用最简单的线性模型, 但由于其只能处理小形变, 能模拟的形变效果有限.之后提出的共旋的线性模型^[6-9]通过形变梯度的分解, 提取旋转分量的方式, 巧妙地处理了关于旋转的问题, 提高了处理大形变的能力, 并且模型非常简单.至今, 在很多形变体仿真相关的研究中, 共旋的线性模型都作为基础的形变模型应用非常广泛.由于在模拟过程中, 当模型的形变非常剧烈或者时间步长非常大时, 模型单元不可避免地出现翻转的状态, 普遍采用修正形变梯度的方式处理翻转的单元, 这种方法不限制材质模型的选择^[15-18], 具有非常直观的效果.随着研究的深入, 简单的均质材质已经不能满足要求, 研究开始转向各向异性材质, 用于处理实际中更加复杂的虚拟物体的建模, 如手指、组织器官、布料纤维等, 目前研究比较多的是正交各向异性材料和横向各向同性材料, 更加一般化的各向异性材料材质还需进一步研究;

(2) 传统的均质材质虽然可以利用各向异性材质在不同方向上模拟不同的形变, 但无法表现复合材料的形变特征.目前, 从非均质材质的离散化方法, 如有限元方法^{[54-56,60,61]}、物质点法^[58]、无网格方法^[59]等, 针对其各自的优势和问题构建了不同的仿真策略.针对微观和宏观两种尺度模型, 提出了均质化方法及其逆过程的设计思路.其中, Kharevych等人^[54]于2009年提出的方法具有开创性的意义, 根据他们的思路, 陆续有文献在不同方面进行改进, 包括从线性模型到共旋线性模型^[55]再到非线性模型^[60]逐步提高处理形变范围; 借鉴其思路的逆均质化方法设计3D打印模型中的微观结构^[57]以及满足不同材质的层级优化结构^[80]的方法.另外, 给定初始数据, 基于数据拟合优化求解材质参数的方法^[60,81]以及通过建立预计算材质数据库的材质设计方法^[61], 材质数据库的构建节省了仿真过程中调试材质的时间, 加快了动画生成过程, 是非常值得推广以及深入研究的课题;

(3) Martin等人^[89]提出的基于样例的弹性材料的仿真方法为材质设计提供了新思路, 为获得复杂的具有艺术指导的形变效果引入了新的研究方向.跟随其仿真思路, 随后的研究者针对不同的形变模型(集中于薄壳体^[86]和三维形变体)、形变效果(弹性、塑性^[87]、碎裂^[88]等)、仿真效率的提升^[90-92]以及使样例与形变目标具有不同拓扑结构^[93]等多个方面展开研究, 扩展了传统材质模型在动画效果上依赖于应力应变曲线的形变路径的方法, 在动画、游戏及影视制作等方面具有很强的应用价值;

(4) 基于样例的弹性材料虽然可以获得非常多样化的形变效果, 但多数用于虚拟形变的动画, 这种调节材质的方式很难将样例材料与物理真实材料相对应.随着3D打印机的出现, 使得材质设计已经从虚拟仿真走向真实的应用场景, 并衍生出很多几何设计构造方面的研究问题, 包括几何优化、结构分析、表面材料效果定制、机构设计、自支撑结构设计、内部结构设计等^[71].多材料打印机的出现, 使生成的模型从单一材质的简单结构逐渐发展到多种材质混合的非均质材质^[80], 以及利用复杂微观结构组合得到宏观的多样化的材料模型^[57,75,81], 硬件设备的进步与软件设计上的不断优化, 使得模型制造越来越符合人们实际需求;

(5) 对于各向异性材质、非均质材质、编辑材质等复杂结构, 是从材质模型角度分析, 从宏观上来说具有复杂的结构特征, 但从微观上来说, 在离散化表示的单元内部, 其应力应变关系可以用线性表示, 也可以是非线性的.

5 总结

本文尝试从动画模拟领域中形变体的本构关系角度, 总结目前仿真中通常采用的均质材质、非均质材质以及更具灵活性的编辑材质.针对形变体材质模型在各个领域的应用, 归纳以下几点材质设计和研究方面需要进一步考虑的问题和未来可能的研究方向.

(1) 真实物体和虚拟物体设计.针对计算机动画、电影、游戏等领域, 需要比较夸张的形变行为, 对应的形变模型是虚拟的, 其运动与真实世界可能存在差异; 而对于3D打印产业来说, 则需要设计合理的具有复杂结构的真实的材质模型.因此, 从材质设计角度而言, 设计虚拟的材料以及真实的材料都具有重要的意义;

(2) 从扩展材质的角度而言, 材质的设计和编辑工作将是未来最主要的研究方向; 而对于编辑材质的稳定性研究、材质的属性编辑范围以及对效果的定性和定量分析还存在一定的难度, 如何获得更稳定的材质模型是未来研究的重点;

(3) 时间效率.基于物理的仿真方法中, 时间效率是不容忽视的因素, 在高逼真动画效果与时间效率二者之间追求平衡, 是各个应用领域需要考虑的重要问题.如何高效地获得满足设计师动画需求的材质, 是未来材质编辑工作面临的挑战.对于形变体仿真而言, 为获得满意的动画效果, 调整参数繁琐而耗时.目前, 已经有一些研究工作是利用数据驱动的方式构建数据库或者采用交互地调整参数方式直观地获得材质模型, 在效率上已经有了大幅度改进.随着机器学习领域的深入研究, 借鉴其研究思路也会给形变体的材质设计带来新的研究方向.