CN103942377B - 一种面向弹性物体制造的逆向形状设计方法 - Google Patents

Abstract

本发明公开了一种面向弹性物体制造的逆向形状设计方法，包括以下步骤：由用户交互指定材料类型、目标形状和一组外力集合；基于实际制造物体的弹性变形特性，使用具有高度预测能力的超弹性材料模型作为物体变形的表征模型；基于渐近数值方法高效精准地求解高度非线性的逆向静力平衡方程，所得形状可在设定的外力条件下变形为用户设计的目标形状；基于核心算法的扩展，支持外力的交互式调整、多目标形状的逆向设计以及正向静力平衡问题的快速求解；使用快速成型等工艺手段制造出符合用户设计目标的实际物体。该方法无缝集成了形状设计和工艺制造，使得用户能够专注于目标形状的设计而不用担心弹性形变对最终产品外形的影响。

Description

一种面向弹性物体制造的逆向形状设计方法

技术领域

本发明主要涉及工业/艺术设计、产品制造/三维打印/快速成型、虚拟物体/角色创建等应用领域，尤其涉及一种面向弹性物体制造的逆向形状设计方法。

背景技术

本发明相关的技术背景简述如下：

一、面向制造的设计

近来研究领域涌现出不少用于创建几何形状的可计算设计工具，它们能在设计过程中施加制造或者物理属性方面的约束，比如结构稳定性（WHITING,E.,SHIN,H.,WANG,R.,OCHSENDORF,J.,AND DURAND,F.2012.Structuraloptimization of3d masonry buildings.ACM Trans.Graph.31,6(Nov.),159:1–159:11.;STAVA,O.,VANEK,J.,BENES,B.,CARR,N.,AND MˇE781CH,R.2012.Stress relief:Improving structural strength of3d printable objects.ACMTrans.Graph.31,4(July),48:1–48:11.）、变形（M.,BICKEL,B.,JAMES,D.L.,AND PFISTER,H.2012.Fabricating articulated characters from skinnedmeshes.ACM Transactions on Graphics(TOG)31,4,47.;CALì,J.,CALIAN,D.A.,AMATI,C.,KLEINBERGER,R.,STEED,A.,KAUTZ,J.,AND WEYRICH,T.2012.3D-printing of non-assembly,articulated models.ACM Trans.on Graphics(TOG)31,6.）、运动（ZHU,L.,XU,W.,SNYDER,J.,LIU,Y.,WANG,G.,AND GUO,B.2012.Motion-guided mechanical toy modeling.ACM Trans.Graph.31,6,127.;COROS,S.,THOMASZEWSKI,B.,NORIS,G.,SUEDA,S.,FORBERG,M.,SUMNER,R.W.,MATUSIK,W.,AND BICKEL,B.2013.Computational design of mechanicalcharacters.ACM Transactions on Graphics(TOG)32,4,83.;CEYLAN,D.,LI,W.,MITRA,N.J.,AGRAWALA,M.,AND PAULY,M.2013.Designing and fabricatingmechanical automata from mocap sequences.ACM Trans.Graph.32,6(Nov.).）和外观（DONG,Y.,WANG,J.,PELLACINI,F.,TONG,X.,AND GUO,B.2010.Fabricating spatially-varying subsurface scattering.ACM Trans.Graph.29,4(July),62:1–62:10.;LAN,Y.,DONG,Y.,PELLACINI,F.,AND TONG,X.2013.Bi-scale appearance fabrication.ACM Trans.Graph.32,4(July).;CHEN,D.,LEVIN,D.I.W.,DIDYK,P.,SITTHI-AMORN,P.,AND MATUSIK,W.2013.Spec2fab:Areducer-tuner model for translating specifications to3d prints.ACM Trans.Graph.32,4(July).）等。在建筑几何设计中也有类似工具被开发出来，但大都关注平面网格方面的设计（LIU,Y.,POTTMANN,H.,WALLNER,J.,YANG,Y.-L.,ANDWANG,W.2006.Geometric modeling with conical meshes and developable surfaces.ACM Trans.Graph.25,3(July),681–689.;LIU,Y.,XU,W.,WANG,J.,ZHU,L.,GUO,B.,CHEN,F.,AND WANG,G.2011.General planar quadrilateral mesh designusing conjugate direction field.ACM Trans.Graph.30,6(Dec.).;VOUGA,E.,M.,WALLNER,J.,AND POTTMANN,H.2012.Design ofself-supporting surfaces.ACM Trans.Graph.31,4(July),87:1–87:11.）。

为了使普通用户能够更便利地进行设计和制造，Mori等人提出了一种方法（MORI,Y.,AND IGARASHI,T.2007.Plushie:An interactive design system forplush toys.ACM Trans.Graph.26,3(July).），能够将基于草图的毛绒玩具形状设计与快速仿真方法结合在一起。Umetani等人最近的工作（UMETANI,N.,KAUFMAN,D.M.,IGARASHI,T.,AND GRINSPUN,E.2011.Sensitive couture forinteractive garment modeling and editing.ACM Trans.Graph.30,4(July),90:1–90:12.;UMETANI,N.,IGARASHI,T.,AND MITRA,N.J.2012.Guidedexploration of physically valid shapes for furniture design.ACM Trans.Graph.31,4(July),86:1–86:11.）在衣物和家具设计中使用敏感度分析技术为设计者提供快速的反馈，使其能获知形状变化对衣物形变和家具结构稳定性的影响。我们的逆向设计工具有相似的目的但是主要针对弹性物体的变形。我们需要处理的数值问题也与现有方法非常不同，主要涉及高度非线性的弹性制造。

与本发明最为相关的工作是一种橡皮气球的计算设计工具（SKOURAS,M.,THOMASZEWSKI,B.,BICKEL,B.,AND GROSS,M.2012.Computational designof rubber balloons.Comp.Graph.Forum31,2pt4(May),835–844.），其中会估算一个充气后能达到目标形状的气球的初始自然状态。他们在静力平衡问题下主要使用Kirchhoff模型来进行薄壳的仿真，并使用传统的牛顿法来最小化能量函数。与之不同，我们使用neo-Hookean模型处理弹性物体，并提出了新颖的数值求解器来快速求得精确解。

二、弹性制造

几何、材料和驱动参数是弹性制造算法中的主要设计变量。Bickel等（BICKEL,B.,M.,OTADUY,M.A.,LEE,H.R.,PFISTER,H.,GROSS,M.,AND MATUSIK,W.2010.Design and fabrication of materials with desireddeformation behavior.ACM Trans.Graph.29,4(July),63:1–63:10.）提出了一种分支界定算法来将基本材料组合成符合特定变形需求的复合材料。实体人脸克隆技术（BICKEL,B.,KAUFMANN,P.,SKOURAS,M.,THOMASZEWSKI,B.,BRADLEY,D.,BEELER,T.,JACKSON,P.,MARSCHNER,S.,MATUSIK,W.,AND GROSS,M.2012.Physical face cloning.ACM Transactions on Graphics(TOG)31,4,118.）在一个过程中同时优化几何、材料和驱动参数，来制造拥有指定形变特性的皮肤。Chen等（CHEN,D.,LEVIN,D.I.W.,DIDYK,P.,SITTHI-AMORN,P.,AND MATUSIK,W.2013.Spec2fab:A reducer-tuner model for translatingspecifications to3d prints.ACM Trans.Graph.32,4(July).）提出了一个“缩减-调整”模型来调制面向三维打印应用的材料结点。Skouras等（SKOURAS,M.,THOMASZEWSKI,B.,COROS,S.,BICKEL,B.,AND GROSS,M.2013.Computational design of actuated deformable characters.ACM Transactions onGraphics(TOG)32,4,82.）同时优化驱动位置和材料属性来实现可变形角色的直接操控。与以上部分工作类似，我们采用了neo-Hookean模型处理弹性形变，但更为重要的是，我们处理的设计问题中物体未变形前的自然状态是未知的，且本发明以一种高效、鲁棒的数值方法对此进行求解。

三、变形控制

控制变形行为是计算机动画领域中的一个重要课题。为了使动画制作更为便捷，许多算法可以生成符合关键帧变形需求的连续动画（HUANG,J.,TONG,Y.,ZHOU,K.,BAO,H.,AND DESBRUN,M.2011.Interactive shape interpolationthrough controllable dynamic deformation.Visualization and Computer Graphics,IEEE Transactions on17,7,983–992.;BARBIˇC,J.,DA SILVA,M.,ANDPOPOVI′C,J.2009.Deformable object animation using reduced optimal control.ACM Trans.Graph.28,3(July),53:1–53:9.;BARBIˇC688,J.,SIN,F.,ANDGRINSPUN,E.2012.Interactive Editing of Deformable Simulations.ACM Trans.Graph.31,4(July).;HILDEBRANDT,K.,SCHULZ,C.,VON TYCOWICZ,C.,AND POLTHIER,K.2012.Interactive spacetime control of deformable objects.ACM Trans.Graph.31,4(July),71:1–71:8.）。为了效率，这些方法往往使用简化模型，因为他们需要控制一系列的动态变形行为，寻找看起来可行但并非完全符合物理真实的结果。而我们需要的是能制造真实弹性物体的方法。

之前也有工作对物体的自然状态进行优化，使之符合特定的变形需求（DEROUET-JOURDAN,A.,BERTAILS-DESCOUBES,F.,AND THOLLOT,J.2010.Stable inverse dynamic curves.ACM Trans.Graph.29,6(Dec.),137:1–137:10.;TWIGG,C.D.,AND KAˇC I′C-ALESI′C786,Z.2011.Optimization for sag-freesimulations.In Proceedings of the2011ACM SIGGRAPH/Eurographics Symposiumon Computer Animation,ACM,New York,NY,USA,SCA’11,225–236.）。逆向头发动态建模（DEROUET-JOURDAN,A.,BERTAILS-DESCOUBES,F.,DAVIET,G.,AND THOLLOT,J.2013.Inverse dynamic hair modeling with frictional contact.ACM Trans.Graph.32,6(Nov.),159:1–159:10.）自动将一个未处理的头发几何转化成为一个可用于仿真的动态头发模型。Pathmanathan等（PATHMANATHAN,P.,CHAPMAN,S.J.,AND GAVAGHAN,D.J.2009.Inverse membrane problems inelasticity.QJMAM Mechanics&Applied Mathematics62,1,67–87.）提出了一种能针对重力进行变形补偿以获得所需仿真设置的方法。Li等（LI,H.,ALHASHIM,I.,ZHANG,H.,SHAMIR,A.,AND COHEN-OR,D.2012.Stackabilization.ACM Trans.Graph.31,6(Nov.).）提出了一种调整物体几何并将它们稳定地堆叠在一起的方法。与之不同，我们需要高保真的物理模型来对弹性制造进行可计算预测。

四、数值跟踪方法

为了求解逆向静力平衡问题，我们的渐近数值方法沿用了传统数值跟踪方法中的路径跟踪思想（ALLGOWER,E.L.,AND GEORG,K.1990.Numericalcontinuation methods,vol.13.Springer-Verlag Berlin.）。经典的方法包括“预测-纠正”方法和“分段线性”方法，其中的基本原则是对非线性解分支的逐步跟踪。虽然这些方法被广泛使用，但它们在迭代步长的确定上有困难，往往使用一个预定义值。太小的步长会造成缓慢的收敛，而太大的步长又会影响结果的精度。

发明内容

本发明的目的在于针对现有技术的不足，提供一种面向弹性物体制造的逆向形状设计方法。

本发明的目的是通过以下技术方案来实现的：一种面向弹性物体制造的逆向形状设计方法，包括如下步骤：

（1）初始设计信息的输入：由用户交互指定材料类型、目标形状和一组外力集合；

（2）超弹性材料模型的使用：基于实际制造物体的弹性变形特性，使用具有高度预测能力的超弹性材料模型作为物体变形的表征模型；

（3）逆向静力平衡问题的求解：基于渐近数值方法高效精准地求解高度非线性的逆向静力平衡方程，所得形状可在设定的外力条件下变形为用户设计的目标形状；

（4）逆向形状设计方法的扩展：基于核心算法的扩展，支持外力的交互式调整、多目标形状的逆向设计、以及正向静力平衡问题的快速求解；

（5）逆向求解结果的制造：基于上述步骤求解所得形状，包括步骤1中所选材料类型（参数），使用快速成型等工艺手段制造出符合用户设计目标的实际物体。

本发明的有益效果是：首次提出了一种面向弹性物体制造的逆向形状设计方法，使得设计者能够更为便捷地设计出面向制造的弹性物体形状。用户仅需指定制造材料属性、目标形状和外力设定，本方法即可进行逆向自动计算，得到可以直接用于制造的模型自然状态，而制造出的实际物体可以准确地在给定外力条件下达到用户所需的目标形状；本方法首次提出一种基于渐近数值方法和neo-Hookean超弹性模型来进行逆向形状设计求解的方法，相比传统的牛顿迭代一类的解法，求解速度提升了几个数量级，使得交互式的形状设计成为可能，并为用户提供流畅的设计体验，从而满足个性化弹性形状设计和制造的应用需求。

附图说明

图1是本发明中逆向形状设计方法的设计场景示意图，图中(a)为输入的目标形状，(b)为实现目标形状的一组外力设定，(c)为面向制造的计算得到的物体自然状态；

图2是本发明中渐近数值方法迭代步骤的示意图；

图3是本发明中路径参数a的收敛半径的示意图；

图4是本发明中基于逆向形状设计并制造的横条模型；

图5是本发明中基于逆向形状设计并制造的盆栽植物模型，图中，(a)为计算所得的自然状态，(b)为重力下的目标形状，(c)为本方法中的静力平衡求解器对其变形结果的准确预测，(d)为制造结果在重力下变形效果图，(e)为将目标形状送去制造后结果在重力下变形展示图；

图6是本发明中基于逆向形状设计并制造的手机支座模型，图中，(a)为计算所得的自然状态，(b)为工作状态下的目标形状以及相应的工作外力设定（两侧施加反向力以抓紧手机），(c)为制造结果，(d)为手机被紧握的效果；

图7是本发明中基于逆向形状设计并制造的衣架模型，图中，(a)为计算所得的自然状态，(b)为重力下的目标形状，(c)为工作状态下的目标形状以及工作外力设定（两侧悬挂力），(d)为制造结果在重力下的变形效果，(e)中展示了衣架悬挂衣服后的变形效果，(f)中展示了悬挂衣服的重量；

图8是本发明中基于逆向形状设计并制造的恐龙模型，图中，(a)为基于多个目标形状（(b)-(e)上排所示）计算所得的自然状态，而制造结果在相同力设定下展现出了与目标形状非常相似的形态（(b)-(e)下排所示）；

图9是两种求解方法错误曲线对比展示，图中，（a）为Levenberg-Marquardt优化方法错误曲线对比展示，（b）为本发明提出的渐近数值方法错误曲线对比展示。

具体实施方式

本发明的核心是基于用户给定材料、目标形状和外力设定高效地进行面向制造弹性物体的自然状态的求解。本发明的核心方法主要分为如下五个部分：初始设计信息的输入、超弹性材料模型的使用、逆向静力平衡问题的求解、逆向形状设计方法的扩展、逆向求解结果的制造。

1.初始设计信息的输入

本方法要求用户输入a）几何形状x，b）制造材料参数，c）一组施加于制造物体上的外力g。基于以上信息本方法会计算出一个形状X，依此制造的物体在外力g的作用下将会变形为用户指定的形状x。

在使用本方法的设计系统中，一个典型的设计场景会首先加载一个用户想达到的目标形状x。然后，用户交互式地指定用于最终制造的材料类型，这些材料都已预先校准好参数（当然，用户也可以直接指定材料参数）。接着，用户需要施加固定约束，并指定外力：重力由形状和材料密度直接决定，而工作外力的位置、大小和方向由用户指定。如附图1所示，用户可以固定鱼形手机支架的底座，并在鱼嘴两侧施加反向的工作外力来抓紧手机。设定好以上信息之后，系统会快速计算出自然状态X，并允许用户对外力g进行交互调整以得到附近的解。最后可导出形状X用于制造。

2.超弹性材料模型的使用

本方法使用有限元模型（四面体网格）来表示可变形物体，其中向量x保存了四面体网格所有顶点的三维坐标，因此对于有N结点的网格模型，x的长度为3N，结点为x_i。根据传统连续介质力学，x表示变形后的形状而X表示未变形前的形状。

由于制造出的弹性物体需要进行大尺度的变形，本方法使用neo-Hookean模型（BONET,J.,AND WOOD,R.D.1997.Nonlinear Continuum Mechanics forFinite Element Analysis.Cambridge University Press.;OGDEN,R.W.1997.Non-linear elastic deformations.Courier Dover Publications.）来精准地计算物体变形所产生的内力f。作为一种超弹性材料模型，它能够很好地预测材料在大尺度形变下的高度非线性弹性行为。

Neo-Hookean模型在几何形变和内部应力之间建立了一种高度非线性的关系，以下是应变的能量密度函数W（x,X）定义：

$W(x,X)=(\frac{\mathrm{\μ}}{2}(J^{-\frac{2}{3}}{I}_c-3)+\frac{\mathrm{\κ}}{2}{(J-1)}^2)$

其中材料参数μ和κ分别是切变模量和体变模量。J和I_c是形变相关量：若F表示变形梯度（），C=F^TF表示右柯西-格林变形张量，则J=det(F)为F的行列式值，而I_c=Tr(C)是C的迹，被称为柯西-格林张量的第一不变量。

从以上这个本构模型可推导出内力函数f的计算公式，具体的推导请参考教科书（BONET,J.,AND WOOD,R.D.1997.Nonlinear Continuum Mechanics forFinite Element Analysis.Cambridge University Press.），以下仅列举基本步骤。首先，第二Piola-Kirchhoff应力张量S可由下式计算：

$S=2\frac{\∂W}{\∂C}=\mathrm{\μ}{J}^{-\frac{2}{3}}I-\frac{\mathrm{\μ}}{3}{J}^{-\frac{2}{3}}{I}_c{C}^{-1}+\mathrm{\κ}(J-1){\mathrm{JC}}^{-1}$

其中I是3×3的单位矩阵。然后，基于S计算第一Piola-Kirchhoff应力张量P=FS。最后，用有限元的逼近方法，结点x_i处的内力f_i为t∈adj(x_i)表示与结点x_i相连的一个四面体，P_t是t中的第一Piola-Kirchhoff应力张量分段常量，而是结点x_i在四面体t处的有效法向。

3.逆向静力平衡问题的求解

3.1经典和逆向静力平衡问题

在经典的静力平衡问题中，已知一个弹性物体的自然状态X和外力g，求解能满足以下静力平衡方程的变形状态x：

f(x,X)+g=0

其中，g为作用在所有四面体网格结点上的外力向量，长度为3N，f是基于自然状态和变形后状态计算内力的函数，其内部采用neo-Hookean模型，具体形式在步骤2中已经给出。

与以上问题不同，本方法中求解的是逆向设计问题，即已知用户指定的目标形状x和外力设定g，求解能满足静力平衡方程的自然状态X。形状设计系统往往要求用户能够得到交互式的反馈和体验，但由于本方法中用于精确表征物理变形的非线性模型有比较复杂的形式，因此如何快速地对逆向静力平衡问题进行求解非常具有挑战性。

3.2渐近数值方法

相比传统的非线性求解器，渐近数值方法在效率、性能和鲁棒性上都要高出许多，因此被成功地应用于非线性工程结构分析等领域（ZAHROUNI,H.,COCHELIN,B.,AND POTIER-FERRY,M.1999.Computing finite rotations ofshells by an asymptotic-numerical method.Comput.Methods Appl.Mech.Eng.175,1,71–85.;LAZARUS,A.,MILLER,J.,AND REIS,P.2013.Continuation ofequilibria and stability of slender elastic rods using an asymptotic numerical method.J.Mech.Phys.Solids61,8.）。同样是迭代方法，但渐近数值方法的收敛速度要远远快于牛顿型方法（图9）。

本方法使用了渐近数值算法的基本框架，但针对逆向形状设计问题和neo-Hookean模型完整建立了底层的数学结构。我们的目标是求解静力平衡方程中的自然状态X，在渐近数值方法中，首先需要建立一个参数化的静力平衡方程（数学上称为同伦）：

f(x,X)+λg=0

其中λ∈[0,1]是负载参数。当λ为0时，上式的解为X=x（完全没有内力）；而当λ为1时，其解恰好为我们之前定义的逆向形状设计问题的解。渐近数值方法的基本思路继承自数值跟踪方法：以迭代的形式跟踪隐式曲线λ(a)（其中λ(0)=0）来改变参数λ。在每一个迭代中，对f(x,X(a))+λ(a)g=0求解，得到新的参数值a和当前解。由于每一个迭代中以一种最优的方式求解参数a，因此用很少的迭代次数就能达到λ(a)=1。

在渐近数值方法的迭代步骤中，定义a₀为当前参数值，则λ₀=λ(a₀)与其相应的解X₀满足先决条件f(x,X₀)+λ₀g=0。由于没有曲线λ(a)的显示表示，本方法使用局部渐近展开来表示λ(a)和其相应解：

$\left[\begin{array}{c}X\left(a\right)\\ \mathrm{\λ}\left(a\right)\end{array}\right]=\left[\begin{array}{c}{X}_0\\ {\mathrm{\λ}}_0\end{array}\right]+\underset{k=1}{\overset{n}{\mathrm{\Σ}}}{(a-a_0)}^k\left[\begin{array}{c}{X}_k\\ {\mathrm{\λ}}_k\end{array}\right]$

其中n是截断阶数，而集合{X_k,λ_k},k=1...n正是当前迭代需要计算的表达系数。在建立局部渐近展开之后，改变a的值直到满足λ(a)=1。当改变a使其远离起点a₀时，X(a)的渐近表示会逐渐远离真实解。当估算的残差超过一定的阈值之后，暂停跟踪并使用牛顿迭代法对当前解进行局部优化（由于离真实解很近，因此这个过程很快）。在得到足够精确的解之后，设a₀=a并继续进行下一个新迭代。重复这个过程直到λ(a)=1满足。在此终点，解X(a)恰好满足静力平衡方程。算法如下所示，具体实例可见图2，图中，x轴表示形变总位移X(a)的向量模大小，y轴表示负载参数λ(a)；根据图中内嵌的目标形状，渐近数值方法首先计算a=0处的渐近展开表示并用其局部地跟踪解分支，直到到达收敛半径为止，然后优化当前点的解并用其计算新的渐近展开表示，这个迭代过程一直持续到λ(a)=1为止。

算法1渐近数值追踪方法：

设X₀=x,λ₀=0,a₀=0；（初始点）

当λ<1

求解多项式系数{X_k,λ_k},k=1...n；

计算a的可信赖跟踪步长；

牛顿法优化X(a)；

设X₀=X(a),λ₀=λ(a),a₀=a；

重复以上循环。

根据以上算法，其中有两个关键步骤：a）渐近展开表示系数的高效求解；b）渐近表示X(a)与真实解的残差计算。

3.3计算渐近表示的系数

3.3.1基本数学原理

在给出系数计算的细节之前，先阐明本方法中高效系数计算的基本原理所在。首先假定内力函数f能够表示成X的一个二次型：

f(x,X)=L₀+L[X]+Q[X,X]

其中L[]和Q[,]分别是线性和双线性向量算子。将前述多项式局部渐近展开表示X(a)带入上式可得二次序列：

$\begin{array}{c}f(x,X\left(a\right))=L_0+L\left[X\right]+Q[X_0,X_0]+(a-a_0)\left(L\right[X_1]+2Q[X_0,X_1\left]\right)\\ +\underset{k=2}{\overset{n}{\mathrm{\&Sigma;}}}{(a-a_0)}^k\left(L\right[X_k]+2Q[X_0,X_k]+\underset{t=1}{\overset{k-1}{\mathrm{\&Sigma;}}}Q[X_t,X_{k-t}\left]\right)\end{array}$

将上式与λ(a)的展开形式一同带入f(x,X(a))+λ(a)g=0，然后通过将具有相同阶数(a-a₀)的项进行相互匹配来建立一组方程。其中，0阶项已经匹配，因为X₀是λ₀处的解，自然满足L₀+L[X₀]+Q[X₀,X₀]+λ₀g=0。而1阶系数需要满足以下方程：

L[X₁]+2Q[X₀,X₁]+λ₁g=0

由于Q中X₀是已知的，因此这是一个关于X₁的线性方程组。然而，X₁和λ₁都是未知量，这是一个拥有3N个方程和3N+1个未知量的线性系统。为了得到一个满秩系统，我们采用Cochelin等提出（COCHELIN,B.1994.A path-followingtechnique via an asymptotic numerical method. Computers & structures53,5,1181–1192.）的方法引入一个额外的约束：

(X(a)-X₀)^TX₁+(λ(a)-λ₀)λ₁=a

在此约束下，参数a衡量的是状态变量的增量(X-X₀,λ-λ₀)在局部切向量(X₁,λ₁)上的投影。类似地，将X(a)和λ(a)的展开形式带入上式并使相同(a-a₀)阶数的系数相等，可为每对(X_k,λ_k)建立一个额外的方程：

$X_k^T{X}_1+{\mathrm{\&lambda;}}_k{\mathrm{\&lambda;}}_1={\mathrm{\&delta;}}_{k1},k=1...n$

其中δ_k1是Kronecker delta：当k=1时δ_k1=1，否则δ_k1=0。此约束与前述1阶方程一起即构成了一个(X₁,λ₁)的满秩线性系统。类似地，如果已求得任意低阶j=1...k-1的系数(X_j,λ_j)，那么如下的k阶方程与以上的约束一起可构成(X_k,λ_k)的满秩线性系统：

$L\left[X_k\right]+2Q[X_0,X_k]+\underset{t=1}{\overset{k-1}{\mathrm{\&Sigma;}}}Q[X_t,X_{k-t}]+{\mathrm{\&lambda;}}_{k}g=0$

此线性系统可以非常快速地求解。事实上，如果线性系统有如下形式：

$A\left[\begin{array}{c}{X}_k\\ {\mathrm{\&lambda;}}_k\end{array}\right]=b$

则矩阵A正是非线性函数f(x,X)+λg在点(X₀,λ₀)处的雅可比矩阵，即：

$A=\frac{\&PartialD;}{\&PartialD;(X,\mathrm{\&lambda;})}(f(x,X)+\mathrm{\&lambda;g}){|}_{(X_0,{\mathrm{\&lambda;}}_0)}$

因此A的计算和每一步迭代结束时进行牛顿法优化所需计算完全一致，可重复利用上一步迭代中的优化计算来提升效率。

3.3.2二次型构造

前述建立(X_k,λ_k)的线性系统的方法有一个重要的前提，即f能够表示成X的一个二次型。然而，使用neo-Hookean模型的内力函数f无法直接写成X的二次型。本方法的一个主要贡献是将此f经过重新转换形成一个二次型表示，中间不存在近似，而是由原始f经过精确的数学推导所得。

首先用一个简单例子来表明基本原理。考虑一个简单函数f(x)=x^3/2。如果x可表示为a的多项式，即则我们的目标是将f(x)也表示为a的多项式形式。虽然f本身不是二次函数，但可通过引入一个新变量y将其写成二次型：

f(x,y)=xy,x=y²

现在，如果y也能表示为a的多项式，即则f(x,y)有如下的多项式形式：

$f(x,y)=\underset{k}{\overset{n}{\mathrm{\&Sigma;}}}{f}_i{a}^k,\mathrm{where}{f}_k=x_0{y}_k+x_k{y}_0+\underset{r=1}{\overset{k-1}{\mathrm{\&Sigma;}}}{x}_r{y}_{k-r}$

为计算f_k，需计算系数y_i,i=0...k。将x(a)和y(a)代入x=y²并使得a的同阶项相等，即得到计算y_k的方程：

$x_k=2y_0{y}_k+\underset{r=1}{\overset{k-1}{\mathrm{\&Sigma;}}}{y}_r{y}_{k-r}$

使用此方程可从低阶到高阶逐步计算y_k，从而得到f_k。以上例子表明对于任意非线性函数f，都可以通过引入辅助变量将之转化为二次型，并基于此二次型将f表示为a的多项式展开。

使用上述原理，我们现在可以计算(X_k,λ_k)。首先，由于逆向形状设计问题中变形后状态x是已知的，因此在变形空间中计算应力张量更为便利。本方法基于第二Piola-Kirchhoff应力张量来计算柯西应力张量σ：

$\mathrm{\&sigma;}=J^{-1}{\mathrm{FSF}}^T=\mathrm{\&mu;}{J}^{-\frac{5}{3}}b-\frac{\mathrm{\&mu;}}{3}{J}^{-\frac{5}{3}}{I}_{c}I+\mathrm{\&kappa;}(J-1)I$

相应地，在分段常量有限元表示中，结点x_i处的内力f_i可计算为： $f_i=\underset{t\&Element;\mathrm{adj}\left(x_i\right)}{\mathrm{\&Sigma;}}{\mathrm{\&sigma;}}^t{n}_i^t$

其中t∈adj(x_i)表示与结点x_i相连的一个四面体，σ^t是t中的柯西应力张量分段常量，而是结点x_i在变形后四面体t处的有效法向。

现在的目标是构造每个结点x_i处的内力f_i的渐近展开表示，即：

$\underset{t\&Element;\mathrm{adj}\left(x_i\right)}{\mathrm{\&Sigma;}}({\mathrm{\&sigma;}}_0^t+\underset{k=1}{\overset{n}{\mathrm{\&Sigma;}}}{a}^k{\mathrm{\&sigma;}}_k^t)n_i^t+({\mathrm{\&lambda;}}_0+\underset{k=1}{\overset{n}{\mathrm{\&Sigma;}}}{a}^k{\mathrm{\&lambda;}}_k)g_i=0$

其中是四面体t中柯西应力张量的渐近数值展开形式的系数。

如前所示，应力张量σ与X之间是一个高度非线性的关系。依照之前f(x)=x^3/2的例子，本方法引入一系列辅助变量，最终按照如下方式计算σ的k阶展开系数σ_k：

${\mathrm{\&sigma;}}_k=\mathrm{\&mu;}(s_k^{\left(5\right)}{b}_0+s_0^{\left(5\right)}{b}_k+\underset{r=1}{\overset{k-1}{\mathrm{\&Sigma;}}}{s}_r^{\left(5\right)}{b}_{k-r})-\frac{\mathrm{\&mu;}}{3}I(s_k^{\left(5\right)}{\left(I_c\right)}_k+\underset{r=1}{\overset{k-1}{\mathrm{\&Sigma;}}}{s}_r^{\left(5\right)}{\left(I_c\right)}_{k-r})+\mathrm{\&kappa;}{J}_{k}I$

其中J_k，(I_c)_k和b_k分别是J，I_c和b的k阶展开系数，而是辅助变量的k阶展开系数（这里为简洁忽略了t）。

至此，已得到f(x,X(a))的渐近表示形式。依照3.3.1中的方式，可构造关于(X_k,λ_k)的线性方程组：

${\mathrm{AX}}_k=-{\mathrm{\&lambda;}}_{k}g+f_{\left(k\right)}^{\mathrm{nl}}$

其中聚集了所有与X_k无关的项。在3.3.3中列举了计算A和的详细公式。如3.3.1中所述，结合以上线性方程组和(X_k,λ_k)上施加的额外线性约束，即可得到一个满秩的线性系统用于求解(X_k,λ_k)。同样，A的分解可以重复利用。

3.3.3完整二次型公式

引入一系列辅助变量，并定义σ和X之间的二次型关系如下：

其中 $s^{\left(5\right)}=J^{-\frac{5}{3}},s^{\left(2\right)}=J^{-\frac{2}{3}},s^{\left(1\right)}=J^{-\frac{1}{3}}.$

基于以上二次型关系，计算k阶多项式系数的完整迭代公式定义如下：

根据上式，σ_k显然与X_k呈线性关系。因此，矩阵A的计算仅仅依赖于辅助变量的初始0阶项，而的计算可由已求得的r阶项(1≤r<k)聚集所得。基于A和本方法使用以下策略计算最终的(X_k,λ_k)：

$\left\{\begin{array}{c}{\mathrm{AX}}_{\left(k\right)}^{\mathrm{nl}}=f_{\left(k\right)}^{\mathrm{nl}}\\ {\mathrm{\&lambda;}}_k=-{\mathrm{\&lambda;}}_1{X}_1^T{X}_{\left(k\right)}^{\mathrm{nl}}\\ X_k=\frac{{\mathrm{\&lambda;}}_k}{{\mathrm{\&lambda;}}_1}{X}_1+X_{\left(k\right)}^{\mathrm{nl}}\end{array}\right.$

3.4残差估计

最后，基于构造的X(a)局部展开解形式，还需要计算残差来定量地反映其和f(x,X)+λ(a)g=0精确解之间的差别精度。这个计算是在不知道精确解的前提下进行估算的。一个有益的观察是：当a-a₀在渐近展开的收敛半径内时，解曲线的两个连续阶（N和N-1）表示之间的差别也很小；而当达到收敛半径之后，两个连续阶曲线之间会很快分散开（见图3）。因此，残差r可计算如下：

$r=\frac{\left|\right|X{\left(a\right)}_{\mathrm{order\; n}}-X{\left(a\right)}_{\mathrm{order\; n}-1}\left|\right|}{\left|\right|X{\left(a\right)}_{\mathrm{order\; n}}-X\left(0\right)\left|\right|}=\frac{\left|\right|X_n{(a-a_0)}^n\left|\right|}{\left|\right|\underset{k=1}{\overset{n}{\mathrm{\&Sigma;}}}{X}_k{(a-a_0)}^k\left|\right|}$

如3.2中所述，r决定了参数a从a₀出发能跟踪的有效距离，并决定了何时需要进行一个新的迭代步骤。在实际试验中，本方法使用r≤1E-6。

4.逆向形状设计方法的扩展

基于逆向弹性形状设计求解的核心算法，我们还做了进一步的扩展以提供更为友好便捷的用户交互。

4.1交互式外力调整

渐近数值方法的一个极具吸引力的特点是能够交互式地调整外力g。这里的基本思路是采用渐近展开快速更新自然状态X（使用不同的外力大小λ）。在渐近数值方法的所有迭代步骤结束时，本方法会得到一个满足f(x,X)+g=0的自然状态X。进一步采用3.3中的算法得到{X,1}处的展开系数，而基于此处的渐近展开，本方法能让用户通过改变λ的值来调整外力的大小。假定更新的外力为，首先求解参数a：

$\stackrel{\&OverBar;}{\mathrm{\&lambda;}}=1+\underset{k=1}{\overset{n}{\mathrm{\&Sigma;}}}{(a-a_1)}^k{\mathrm{\&lambda;}}_k$

这里使用多项式求根方法求得最接近a₁（λ为1对应的参数a值）的根。最后，基于求解所得a值更新自然状态X。整个计算只涉及到低阶多项式的求根公式和级数的快速计算，非常地高效。

当用户调整使其偏离λ=1时，展开式的预测精度将会下降。使用3.4中的相同策略，本方法会在残差估算超过阈值r=1E-6时开始新的迭代步骤来更新X。

以上策略只能更新外力大小λ而不能改变外力方向，而对外力方向调整的方法需要更高一些的计算代价。当用户调整外力方向，从g₀变化到g₁时，使用新的渐近数值迭代步骤求解以下新的参数化逆向静力平衡方程：

f(x,X)+λ(g₁-g₀)+g₀=0

当λ=0时，解已知。使用渐近数值方法跟踪解曲线直到λ=1为止。如果更新的外力g₁与g₀较为接近，则本方法仍然能够提供交互式的反馈。否则，用户需要几秒钟来等待形状更新。

4.2多目标逆向形状设计

前述方法都使用一个目标形状作为输入。然而，一个制造物体往往会以多种方式工作，因此需要施加多种外力到一个物体上使其达到不同的目标形状（见图8）。下面形式化此问题：给定T个目标形状x_t,t=1...T，连同相应外力g_t,t=1...T，寻找一个可制造的自然状态使其在每个外力g_t作用下能够分别变形为目标形状x_t。这需要解一个过度约束的方程系统（有可能无解），也就是说很多情况下这样的自然状态事实上是不存在的。

因此，本方法寻找一个自然状态使得其在给定外力条件下能尽可能接近每一个目标形状x_t。一个直接的方法是将其作为一个优化问题：找到一个形状使其最小化以下力平衡残差之和：

$\stackrel{\&OverBar;}{X}=\arg {\min }_X\underset{t=1}{\overset{T}{\mathrm{\&Sigma;}}}\left|\right|f(x_t,X)+g_t\left|\right|$

然而，这个非线性优化问题会碰到各种性能和收敛性方面的问题（见表2）。为了实现一个响应式的设计工具，本方法采用另外一种快速估算的策略。渐近展开式中的系数X_i,i=0...n是独立向量，它们构成了一个表示的局部缩减子空间。在展开式中，这些生成向量的大小是单项式(a-a0)ⁱ，可以进一步将之放松为任何可能值来估算

具体来看，对于每一个目标形状x_t和相应的外力g_t，可以快速地独立求解一个自然状态在每个解处，计算展开系数并将其表示为X_t,i。将所有系数合并在一起得到一个基矩阵：

U=[X_1,0…X_1,nX_2,0…X_2,n…X_T,0…X_T,n]

然后基于这个基矩阵来找到一个形状使其与所有独立解之间最为接近（最小二乘）。即求解：

$q=\arg {\min }_q\underset{t=1}{\overset{T}{\mathrm{\&Sigma;}}}{\left|\right|\mathrm{Uq}-{\stackrel{\&OverBar;}{X}}_t\left|\right|}_2^2$

并求得以上计算只针对U求最小二乘解，比之前的方法要快捷许多。

由于使用近似，所求的变形结果很可能无法与每一个目标形状x_t完全一致。因此用户可使用以下的静力平衡求解方法来快速呈现X的变形结果与目标形状x_t之间的差别。

4.3经典静力平衡问题求解

虽然3中的数学方法是针对逆向形状设计问题，但渐近数值方法本身是个一般性的框架，能够有效地求解非线性系统，包括经典的静力平衡问题f(x,X)+g=0（其中自然状态X和外力g为已知，要求解变形后的形状x）。求解经典静力平衡问题不仅能够加速材料参数的校准过程，并且在多目标设计中提供了变形结果差别的快速反馈。更广泛来看，经典静力平衡问题出现在了许多可计算设计的研究工作（UMETANI,N.,KAUFMAN,D.M.,IGARASHI,T.,AND GRINSPUN,E.2011.Sensitive couture for interactive garment modeling and editing.ACM Trans.Graph.30,4(July),90:1–90:12.）和参数选择（MIGUEL,E.,BRADLEY,D.,THOMASZEWSKI,B.,BICKEL,B.,MATUSIK,W.,OTADUY,M.A.,ANDMARSCHNER, S. 2012. Data-driven estimation of cloth simulation models. Comp.Graph.Forum31,2(May),519–528.）中，因此一个快速求解方法能很好地提升这些方法的性能。

使用渐近数值方法求解静力平衡问题的方式与步骤3中一致，这里只说明不同的部分。首先，由于需要求解的是变形后状态x，其对应的展开形式为：

$\left[\begin{array}{c}x\left(a\right)\\ \mathrm{\&lambda;}\left(a\right)\end{array}\right]=\left[\begin{array}{c}{x}_0\\ {\mathrm{\&lambda;}}_0\end{array}\right]+\underset{k=1}{\overset{n}{\mathrm{\&Sigma;}}}{(a-a_0)}^k\left[\begin{array}{c}{x}_k\\ {\mathrm{\&lambda;}}_k\end{array}\right]$

另一个主要不同在于内力函数f的计算在这里使用的是第一Piola-Kirchhoff张量，P=FS。结点x_i处的内力f_i基于材料（未变形）空间中的有效法向进行计算：

$f_i=\underset{t\&Element;\mathrm{adj}\left(x_i\right)}{\mathrm{\&Sigma;}}{P}^t\stackrel{\&OverBar;}{n_i^t}$

其中t∈adj(x_i)表示与结点x_i相邻的四面体，P_t是第一Piola-Kirchhoff应力张量在t中的分段常量，而是结点x_i在未变形四面体t处的有效法向。最后，构造P的展开，并基于其在每一个迭代步骤中计算展开系数{x_k,λ_k},k=1...n。

渐近数值方法提供了快速的静力平衡求解性能。如表3所示，与之前使用运动阻尼的静力平衡求解方法（UMETANI,N.,KAUFMAN,D.M.,IGARASHI,T.,AND GRINSPUN,E.2011.Sensitive couture for interactive garment modeling andediting.ACM Trans.Graph.30,4(July),90:1–90:12.）相比，它提供了7倍的平均提速。同时，之前的方法显式地简化了一个二维仿真网格来达到交互级别的静力平衡求解速度。但是对于三维弹性物体的实际制造应用，必须使用高分辨率的网格来保证精度，因此直接应用Umetani等人的方法是不合适的。

4.3.1完整二次型公式

以下列举了求解经典静力平衡问题的完整二次型公式。首先引入一系列辅助变量来定义P和x之间的二次型关系：

基于以上二次型关系，计算k阶多项式系数的完整迭代公式定义如下：

其中，矩阵A和向量的计算，以及对x和λ的线性系统求解，都与3.3中的逆向形状设计方法类似，此处不再重复。

5.逆向求解结果的制造

导出3和4中计算所得形状（即物体自然状态），并将其转换为制造所需的模型数据格式，得到工艺制造所需的几何模型。然后基于用户指定的材料类型和参数，选择特定的工艺制造手段（三维打印、硅胶复模等）来制造出设计好的几何模型。最后，基于用户事先设定好的约束和外力，将制造好的物体在不同场景不同设定条件下变形为用户所需的目标形状。

实施实例

发明人设计并制造了7个不同的模型（统计数据见表1）来评价和验证本发明中提出的逆向形状设计方法。

表1：实验模型的统计数据

模型	顶点数目	四面体数目	目标数目	外力数目
					横条	4552	19552	1	1
植物	14842	47077	1	1
					手机支架	18753	72355	1	3
衣架	24323	98131	2	1/3
					老鹰	19307	71235	1	1
恐龙	10673	32953	4	1/3/2/2
					三杈	5093	24478	1	7

这些弹性物体的制造方式为：首先三维打印出固体模型用于制作硅胶模，然后用PU8400弹性材料进行铸造得到最终的弹性物体。以下是具体的实验分析。

1.逆向静力平衡求解的运行时间

发明人在一台配备了Intel i7-3770K中央处理器的桌面个人电脑上实现了本发明的实施实例。表2列举了实验中各实例的运行时间和迭代步骤。其中，Levmar求解器的时间统计包括两种不同的初始值设定。另外，对于老鹰和恐龙模型，Levmar求解器不能收敛到正确的解。

表2：逆向静力平衡求解的计算时间统计

对于简单的形状（如图4中的弹性条），渐近数值方法仅仅用了2.38秒进行求解。手机支架模型（图6）则需要更长的运行时间，因为施加的外力较大（40牛顿），所以需要更多的计算时间来使外力从0增长到40。对于多目标形状设计（包括图7的衣架和图8的恐龙），运行时间里包括了每一个目标单独的求解时间以及最终基于此进行自然状态估算的时间。

在表2中还比较了渐近数值方法和传统的牛顿型方法在求解逆向静力平衡问题上的性能差别。其中，牛顿型方法求解器基于Levenberg-Marquardt方法（LOURAKIS,M.I.2010.Sparse non-linear least squares optimization forgeometric vision.In European Conference on Computer Vision,vol.2,43–56.）实现。根据实验结果所示，对于所有测试例子，本方法要比牛顿型方法快2-3个数量级。另外，由于牛顿型方法的效果取决于不同的初始值选择，表2中也报告了使用两种不同初始值的Levenberg-Marquardt求解器运行时间。第一种策略使用目标形状作为初始值，而第二种策略用正向静力平衡求解器以如下方式计算一个初值：将目标形状x视为自然状态，并将外力g反向得到-g，然后求解一个静力平衡状态x'以满足f(x,x')-g=0，最后使用x'作为初值进行性能测试。这个策略往往能够提高Levenberg-Marquardt方法的求解效率，但需要额外的初值计算时间。更为重要的是，即便没有算上额外的计算时间，本方法仍然有1-2个数量级的速度优势。

2.经典静力平衡求解的运行时间

表3给出了本方法中静力平衡求解的运行时间和Umetani等人提出的运动阻尼仿真方法的比较。总的来说，本方法有3-10倍的速度优势。

表3：正向静力平衡求解的计算时间统计

3.物理验证

首先，发明人使用一个简单的条状模型进行验证。我们要求模型在重力下的变形目标形状是一个水平横条。图4中的照片展示了本方法计算出的用于制造的物体自然状态及其在重力下的形态，可以看见制造出的物体能很忠实地变形为水平形状。求解器受外力作用下的鲁棒性在图6中展示。在设计阶段，我们在鱼嘴的两侧施加了40牛顿的外力来使其张大以紧握手机（图6b）。最终制造的手机支架模型显示于图6c中，而外力作用下的真实变形在图6d中展示。

多目标优化算法的验证在图7中展示。重力和工作外力（衣架两侧各0.4牛顿）作用下的两个目标形状如图7b和图7c所示。求解所得自然状态如图7a所示。为验证算法的精度，发明人在制造出的弹性衣架上挂了一件90克的衣服，得到的变形效果（图7e）与设计目标（图7c）在视觉上完全一致。

另一个多目标优化的例子在图8中展示，其中的四个不同的设计目标，分别是重力、以及在手上、脖子上和头上施加外力作用下的变形结果。制造出的物体在这些外力条件下可以展现出与设计目标在视觉上非常相似的变形效果。

Claims (5)

1.一种面向弹性物体制造的逆向形状设计方法，其特征在于，包括如下步骤：

(1)初始设计信息的输入：由用户交互指定材料类型、目标形状和一组外力集合；具体包括以下子步骤：

(1.1)用户辅助交互，指定预先校准的材料类型或者直接给定材料参数，所述材料参数包括切变模量、体变模量、密度；

(1.2)用户提供所需目标形状模型，以实现特定功能或视觉效果展示；

(1.3)用户辅助交互，选择模型中特定部分，施加固定约束或者工作外力，调整外力的大小和方向；

(2)超弹性材料模型的使用：基于实际制造物体的弹性变形特性，使用具有高度预测能力的超弹性材料模型作为物体变形的表征模型；

(3)逆向静力平衡问题的求解：基于渐近数值方法高效精准地求解高度非线性的逆向静力平衡方程，所得形状可在设定的外力条件下变形为用户设计的目标形状；

(4)逆向形状设计方法的扩展：基于核心算法的扩展，支持外力的交互式调整、多目标形状的逆向设计、以及正向静力平衡问题的快速求解；

(5)逆向求解结果的制造：基于上述步骤求解所得形状，包括步骤(1)中所选材料类型或参数，使用快速成型工艺手段制造出符合用户设计目标的实际物体。

2.根据权利要求1所述面向弹性物体制造的逆向形状设计方法，其特征在于，所述步骤(2)包括以下子步骤：

(2.1)基于逆向形状设计问题所涉及的大尺度弹性变形特性，使用超弹性材料模型neo-Hookean模型来表征高度非线性的弹性形变；

(2.2)基于neo-Hookean模型，在几何形变与模型内部应力之间建立一种高度非线性的关系，并通过两个参数切变模量和体变模量来适应具有不同弹性属性的材料。

3.根据权利要求1所述面向弹性物体制造的逆向形状设计方法，其特征在于，所述步骤(3)包括以下子步骤：

(3.1)建立静力平衡方程，使得几何形变产生的模型内部应力与用户指定的外力之间达到静态平衡状态，即相互抵消；

(3.2)基于渐近数值方法和neo-Hookean模型，得到步骤(3.1)中静力平衡方程参数形式解的渐近展开表示形式，并计算其中所涉及到的各项系数；

(3.3)基于步骤(3.2)中获得的渐近表示解，使用数值跟踪方法最终得到平衡态的结果形状。

4.根据权利要求1所述面向弹性物体制造的逆向形状设计方法，其特征在于，所述步骤(4)包括以下子步骤：

(4.1)基于步骤(3)中的计算结果，变更渐近展开形式以支持最终结果附近的数值跟踪，使得用户对外力的交互式调整能得到快速及时的反馈；

(4.2)基于步骤(3)的计算结果，使用松弛子空间的投影法，得到多个目标形状在多组外力约束下的最优自然状态；

(4.3)求解正向静力平衡问题，得到自然状态施加给定外力之后的形变结果。

5.根据权利要求2所述面向弹性物体制造的逆向形状设计方法，其特征在于，所述步骤(5)包括以下子步骤：

(5.1)基于步骤(3)和步骤(4)中求解得到的结果，即物体自然状态，转换格式得到工艺制造所需几何模型；

(5.2)基于步骤(1)中提供的材料参数制造步骤(5.1)中所得几何模型；

(5.3)基于步骤(1)中设定的约束和外力，将步骤(5.2)中制造的物体变形为目标形状。