一种薄壳体3D打印优化方法

摘要

本发明公开了一种薄壳体3D打印优化方法，包括以下三个步骤，1、确定每个顶点位置能够达到的最大厚度；2、对模型进行分片操作，对区域边界相交处进行扩展得到过渡区域,对于每一个独立的片和每一个顶点均采用一个厚度参数；3、根据需求，对厚度参数进行优化，得到优化后的三维模型。通过本文提出的方法框架可以快速自动地计算在特定使用场景下，具有足够结构强度并且使用较少材料的可3D打印模型。方法明确，速度较快，结果鲁棒。可用于计算机辅助设计，生成3D打印模型等领域。

说明书

技术领域

本发明涉及计算机图形与辅助设计领域，3D打印技术领域，特别是 涉及薄壳体3D打印优化厚度参数得到可打印实体的方法。

背景技术

3D打印(三维打印)是增材制造技术(AM,additiveManufacturing) 的俗称，通过数字化设计模型，运用粉末状材料(塑料或金属)，通过热 熔粘合技术采用叠加成型的方式生成三维实体。智能数字化技术在3D打 印过程中扮演重要作用。借助计算机实现，能够只能的实现用户的意图并 生成可打印的实体。随着技术不断发展，3D打印机可以采用多种材料， 包括粉末，所料，金属等，以及在特殊应用场景下的巧克力，人体干细胞 等。通过分层加工，叠加成型的方式实现能够制造复杂模型。基于以上优 点，3D打印技术在工业制造，航空航天，机器人，生物医学等领域快速 发展。也由于3D打印机制造成本的降低，让这一技术在广大普通用户中 得到迅速推广。

在制造模型的过程中，用户借助3D建模软件(如Maya,3DSMax等) 建立网格模型。通过此种方法建立出来的模型通常是边界曲面模型，而这 样生成的实体通常是非常薄不能够直接用来打印。为了解决这个问题，大 多数3D打印软件采用以下两种方法：1、将曲面表面洞进行填充，整体曲 面封闭生成实体；2、对整体模型设定一个厚度，生成等距曲面封闭之后 形成实体。然而对于方法1，会显著改变模型表面结构，同时这样打印会 将内部空间全部打印成实体，在打印过程中会增加打印的时长，消耗更多 的材料。而对于方法2，厚度确定的过小会使得打印出的物体难以满足实 际需求，承受一定的外力；厚度过大会影响整体的外观。

近几年有一些研究者专注研究如何增强打印物体的结构强度。STAVA 和VANEK提出通过改善物体形状和增加支撑来增加物体的结构强度，同 时减少打印消耗的材料，参见STAVA,O.,VANEK,J.,BENES,B.,CARR, N.,ANDMˇECH,R.：Stressrelief:Improvingstructuralstrengthof3d printable，566，2012，ACMTrans.Graph.31,4,48:1–48:11。Wang等通过 构建框架支撑状结构实现坚固结构和节省材料的平衡。参见WANG,W., WANG,T.Y.,YANG,Z.,LIU,L.,TONG,X.,TONG,W.,DENG,J.,CHEN,F., ANDLIU,X.2013.Cost-effectiveprintingof3dobjectswithskin-frame structures.ACMTrans.Graph.32,6,177:1–177:10。然而以上方法对于薄壳 状结构的物体并不适用。由于薄壳状物体内外表面都可见，增加支撑的做 法会影响外观形状，同时难以满足用户设计时的使用需求。

发明内容

本发明提供了一种薄壳体3D打印优化方法，能够实现根据输入的模 型和需要达到的结构强度，自动计算每点厚度参数，并生成可打印的密闭 网格模型。

为了达到上述目的，本发明采用以下技术方案：一种薄壳体3D打印 优化方法，包括以下步骤：

(1)根据三维网格模型几何结构，确定每个顶点位置能够达到的最 大厚度和最小厚度；

(2)根据模型的几何结构，对模型进行分片操作，得到多个分片区 域，对分片区域边界相交处进行扩展得到过渡区域。根据步骤1得到的最 大厚度和最小厚度，对于分片区域中每一个独立的片采用一个厚度参数， 对于过渡区域中每一个顶点采用一个厚度参数；

(3)根据模型需要达到的使用中的应力强度要求，承受外力的情况 和打印所使用的材料，使用特定有限元模型模拟仿真和敏感度分析技术， 通过交替迭代对步骤2中的厚度参数进行优化，求得网格模型上每一个顶 点的最佳厚度，并生成优化后的三维模型。

进一步地，步骤(1)中确定每个顶点最大厚度方法，具体为：根据输入 三维网格模型，用M表示，求得其包围盒在x,y,z轴的尺寸为(a,b,c)，则 模型每个点最大厚度设定为：

T_Max＝0.05*(a+b+c-min{a,b,c}-max{a,b,c})

最小厚度根据3D打印机能够打印的最薄厚度确定：

T_min＝0.5mm

根据模型网格几何特征对M上每一点进行厚度修正，采用averaged distancefunctions函数计算网格模型M的距离场函数d(X,M)，表示空间中 点X到网格模型M的距离，点x是位于模型M上的点，F(x,t)表示通过点x 的积分曲线，表示点x指向F(x,t)的向量，N(x)表示模型M上x点法 向方向。

在曲面与法向相同方向一侧，最大厚度为：

同理在与法向相反一侧，最大厚度为：

对于M上点x，修正后的最大厚度T_G(x)表示为

$T_G\left(x\right)=2\min \{T_G^{-}\left(x\right),T_G^{+}\left(x\right),0.5T_{Max}\}$

进一步地，步骤(2)中根据模型几何特征，采用fuzzycut方法对模型 进行分片，具体为：

扩展分片之间的边界：输入模型的顶点数用n来表示，设定扩展边界 的迭代次数N_expand(至少进行一次)，在每次迭代过程中，将过渡区域中所 有顶点的单环邻域内的顶点添加到过渡区域。经过N_expand迭代后，得到的 过渡区域表示为B_T；

对于分片区域中每一个独立的片的厚度参数以及过渡区域中每一个 顶点的厚度参数采用如下办法获得：

过渡区域中顶点个数为n_T，每一个顶点赋予一个厚度参数，用向量β 表示。β中每个分量对应相应顶点的厚度参数，表示为

在进行边界扩展后，剩下的部分被分割成彼此不相邻的部分，称为分 块区域，表示为B_S，不相邻的分块数为n_S，每一个单独的分块区域采用同 样的分片厚度，用向量α表示。α每个分量对应相应分块区域的厚度，表 示为

表示α的厚度上限，表示α的厚度 下限。对于α_i的厚度参数变化范围由其分片区域包含的所有顶点集合决定，即

$T_{\min }<{\alpha }_i<\min \left\{T_G\left(x\right)\right|x\in {\mathrm{sx}}_{{\alpha }_i}\}$

表示β的厚度上限，表示β的厚度 下限。对于β_i的厚度参数变化范围由其对应的模型顶点决定，即

$T_{\min }<{\beta }_i<T_G\left(x_{{\beta }_i}\right)$

由此，根据分片厚度α及预计算的每个顶点厚度参数β可确定模型M 上每个顶点的厚度t。即对于给定的α,β可以求得t。

进一步地，步骤(3)中包括步骤：

(4.1)将模型M拆分为二次壳体单元和线性平板单元，加快计算速 度整体刚度矩阵，可以表示为厚度参数α,β的多项式。

(4.2)通过有限元方程建立厚度参数α,β变化对模型每个顶点von Mises力σ变化的影响

根据有限元方程，在外力F作用下，根据整体刚度矩阵K_sys可以求解 模型每个顶点的位移U：

K_sysU＝F

对于在模型M上顶点i处的厚度t_i变化对整体位移产生的影响可以表 示为：

$\frac{\partial F}{\partial t_i}=\frac{\partial K_{sys}}{\partial t_i}U+K_{sys}\frac{\partial U}{\partial t_i}$

由于外力保持不变，不随着厚度t_i变化，上式可简化为：

$\frac{\partial U}{\partial t_i}=-K_{sys}^{-1}\frac{\partial K_{sys}}{\partial t_i}U$

结合vonMises计算公式，可以根据位移对厚度参数的变化求得

(4.3)建立优化目标函数，通过交替迭代优化过程，分别优化α,β。

优化目标是在模型承载指定外力的条件下，保证材料不断裂，满足结 构强度，保持所需要的打印材料最小并且保持生成曲面光滑。即可用如下 形式表达为：

$\underset{{\mathrm{\Delta t}}_i^k}{\min }\underset{i=1}{\overset{n}{\Sigma }}(t_i^k+{\mathrm{\Delta t}}_i^k)s_i$

$\left(\begin{array}{cc}s.t.& {\sigma }_v^k+\frac{\partial {\sigma }_v}{\partial t}\left(t_i^k\right){\mathrm{\Delta t}}_i^k<{\sigma }_{\max },\forall v\in M,\end{array}\right)$

$t^{\min }<t_i^k+{\mathrm{\Delta t}}_i^k<t_i^{\max },i=\mathrm{1...}n,$

其中s_i表示顶点i邻接三角形面片面积和的控制曲面光滑程度， t^min由3D打印机打印最薄厚度确定，w_i表示其邻域顶点的权重，若与点i 邻接的顶点数为N_ring，则σ_max由打印所选取的材料决定。

等价于密闭模型的体积V_M。

通过迭代优化过程求解，在第k次迭代过程中，阶段1和阶段2的求解 分别为：

阶段1：

$\underset{{\mathrm{\Delta t}}_i^k}{\min }\underset{i=1}{\overset{n}{\Sigma }}({\alpha }_i^k+{\mathrm{\Delta \alpha }}_i^k)s_i$

$\left(\begin{array}{cc}s.t.& {\sigma }_v^k+\frac{\partial {\sigma }_v}{\partial t}\left(t_i^k\right){\mathrm{\Delta t}}_i^k<{\sigma }_{\max },\forall v\in M,\end{array}\right)$

${\alpha }_i^{\min }<{\alpha }_i^k+{\mathrm{\Delta \alpha }}_i^k<{\alpha }_i^{\max },i=\mathrm{1...}{n}_s,$

在阶段1求解过程中，固定β值不变，求解线性优化问题的解为 α^k+Δα^k，通过求解拉普拉斯方程，得到β^k。

阶段2：

$\underset{{\mathrm{\Delta t}}_i^k}{\min }\underset{i=1}{\overset{n}{\Sigma }}({\beta }_i^k+{\mathrm{\Delta \beta }}_i^k)s_i$

$\left(\begin{array}{cc}s.t.& {\sigma }_v^k+\frac{\partial {\sigma }_v}{\partial t}\left(t_i^k\right){\mathrm{\Delta t}}_i^k<{\sigma }_{\max },\forall v\in M,\end{array}\right)$

$t^{\min }<{\beta }_i^k+{\mathrm{\Delta \beta }}_i^k<{\beta }_i^{\max },i=\mathrm{1...}{n}_t,$

在阶段2的求解后，得到β^k+Δβ^k，此处同样通过求解拉普拉斯方程， 求得α^k+1，作为k+1次迭代的初值。

阶段1和阶段2依次进行，循环求解，直到约束条件满足(保证材料 不断裂，满足结构强度，保持所需要的打印材料最小并且保持生成曲面光 滑)，且优化目标V_M不在下降，则停止优化过程，得到最终结果t^F。

根据求解得到的网格M上每个顶点的厚度参数根据积分曲线上 的位置，向两侧各处取得生成曲面位置，即得到优化后的三维模型。

本发明的有益效果在于：

1、将整体模型划分成分片区域和过渡区域，减少了优化变量，使得 优化更为简便；

2、采用特殊的二次壳体单元和线性平面应力单元组合进行模拟仿真， 加快了计算速度；

3、使用敏感度分析技术将优化问题简化为线性规划，进一步采用迭 代优化技术，得到此优化问题的可行解；在不影响外观形状的前提下，将 材料消耗最小化，同时满足用户设计时的使用需求，包括受力情况等。

附图说明

图1为本发明的技术方案流程图。

图2为本发明初始模型示意图；

图3为划分区域示意图；

图4为优化后模型结果图。

具体实施方式

如图1所示，一种对于薄壳体3D打印厚度参数优化方法，包括以下 三个步骤，1、根据需要3D打印的三维网格模型几何结构，确定每个顶点 位置能够达到的最大厚度；2、根据模型的几何结构，对模型进行分片操 作，得到分片区域，对区域边界相交处进行扩展得到过渡区域,对于分片区 域中每一个独立的片采用一个厚度参数，对于过渡区域中每一个顶点采用 一个厚度参数；3、根据模型需要达到的使用需求，承受外力的情况和打 印所使用的材料，使用特定有限元模型模拟仿真和敏感度分析技术，通过 交替迭代技术，优化网格模型上每一个顶点需要达到的厚度并生成封闭的 三维模型。

下面结合实施例以及附图2-4对本发明进行详细说明。

现具体介绍本方法的三个步骤：

1)根据需要3D打印的三维网格模型的几何结构，确定每个顶点位置 能够达到的最大厚度的方法如下：

根据用户输入的三维网格，用M表示，求得其包围盒的在x,y,z轴的 尺寸为(a,b,c)。由于过大的厚度会造成有限元模型仿真的失真，根据整 体模型的尺寸估算薄壳状物体能够达到的最大厚度，则模型每个点最大厚 度设定为：

T_Max＝0.05*(a+b+c-min{a,b,c}-max{a,b,c})

最小厚度根据3D打印机能够打印的最薄厚度确定：

T_min＝0.5mm

根据模型网格几何特征对M上每一点进行厚度修正，采用averaged distancefunctions函数计算网格模型M的距离场函数d(X,M)，表示空间中 点X到网格模型M的距离，点x是位于模型M上的点，F(x,t)表示通过点x 的积分曲线，表示点x指向F(x,t)的向量，N(x)表示模型M上x点法 向方向，参见PengJianbo,Kristjansson,DanielZorin,Denis(2004). Interactivemodelingoftopologicallycomplexgeometricdetail.ACM TransactionsonGraphics(TOG)。

在曲面与法向相同方向一侧，最大厚度为：

同理在与法向相反一侧，最大厚度为：

对于M上点x根据几何形状计算的最大厚度T_G(x)表示为

$T_G\left(x\right)=2\min \{T_G^{-}\left(x\right),T_G^{+}\left(x\right),0.5T_{Max}\}$

此处采用averageddistancefunctions由于此函数的积分曲线能够根据 网格的几何特征改变积分方向，不同点间的积分曲线不会自交，有效避免 曲面凹凸变化造成的影响。而若采用法向方向构建等距曲面的方法则无法 避免在曲面凹处的自交情况。T_G(x)的取值是为了保证有限元仿真模型的正 确适用。

2)根据模型几何特征，采用fuzzycut方法对模型进行分片，参照Katz Sagi，TalAyellet，2003，Hierarchicalmeshdecompositionusingfuzzy clusteringandcuts，ACMTransactionsonGraphics(TOG)。分片结果如图 3所示。结合图描述一下

获得的分割区域的个数可以由用户指定或者程序自动做出判断。在划 分区域的基础上，用户可以自主调节，通过更改默认参数，扩大或者缩小 相应区域。

扩展分片之间的边界，输入模型的顶点数用n来表示，设定扩展边界的 迭代次数N_expand(至少进行一次)，在每次迭代过程中，将过渡区域中所有 顶点的单环邻域内顶点添加到过渡区域。经过N_expand迭代后，得到的过渡 区域表示为B_T，过渡区域中顶点个数为n_T，每一个顶点赋予一个厚度参数， 用向量β表示。β每个分量对应相应顶点的厚度参数，表示为在进行边界扩展后，剩下的部分被分割成彼此不相邻的部分，称为分块区 域表示为B_S，不相邻的分块数为n_S，每一个单独的分块区域采用同样的厚 度，用向量α表示。α每个分量对应相应分块区域的厚度，表示为

表示α厚度上限，表示α厚度参数 的下限。对于α_i的厚度参数变化范围由其分片区域包含的所有顶点集合 决定，即

$T_{\min }<{\alpha }_i<\min \left\{T_G\left(x\right)\right|x\in {\mathrm{sx}}_{{\alpha }_i}\}$

表示β厚度上限，表示β厚度参数 的下限。对于β_i的厚度参数变化范围由其对应的模型顶点决定，即

$T_{\min }<{\beta }_i<T_G\left(x_{{\beta }_i}\right)$

有此建立起α,β到模型M上每个顶点的厚度t的对应。即对于给定的 α,β可以求得t，反之亦然。

采用fuzzycut划分曲面方法能够有效的利用几何信息。分块区域的作 用在于保持模型整体特征的完整，降低优化部分的自由度，加快算法收敛， 过渡区域的作用在于光滑的链接不同的分块区域，不会在模型外形上观察 到由于厚度变化导致的明显差异。

3)图2为初始模型示意图，从图中可以看出，模型需承受的外力 为5N，本实施例中，打印所使用的材料的弹性模量用E表示，泊松比用ν 表示。

采用三角形二次壳体单元和线性平面应力单元组合模型，将整体刚度 矩阵可以表示为厚度参数α,β的多项式。采用敏感度分析技术，通过有限 元方程建立厚度参数α,β变化对模型每个顶点vonMises力σ变化的影响。

向量t表示模型M上每个顶点的厚度，则对于单元i,j,m的厚度t^e可以 表示为：

$t^e=\frac{t_i+t_j+t_m}{3}$

对于典型的三角形单元，节点i,j,m按逆时针顺序编号，使用面积坐标 表示表示，则其形状函数N^P可以表示为：

$N^P=\left(\begin{array}{cccccc}{N}_i^P& 0& N_j^P& 0& N_m^P& 0\\ 0& N_i^P& 0& N_j^P& 0& N_m^P\end{array}\right)$

其中，可以表示为：

$N_i^P=\frac{a_i+b_{i}x+c_{i}y}{2\Delta }$

其中，a_i,b_i,c_i计算公式如下：

$\left(\begin{array}{c}{a}_i=x_i{y}_m-x_m{y}_j\\ b_i=y_j-y_m\\ c_i=x_m-x_j\end{array}\right)$

其它系数可以通过下标i,j,m的顺序循环而求得，其中

$2\Delta =\det \left(\begin{array}{ccc}1& x_i& y_i\\ 1& x_j& y_j\\ 1& x_m& y_m\end{array}\right)$

由此可以计算矩阵B^P，

$B^P=\left(\begin{array}{ccc}{B}_i^P& B_j^P& B_m^P\end{array}\right)$

$B_i^P=\left(\begin{array}{cc}\frac{\partial N_i^P}{\partial x}& 0\\ 0& \frac{\partial N_i^P}{\partial y}\\ \frac{\partial N_i^P}{\partial y}& \frac{\partial N_i^P}{\partial x}\end{array}\right)=\frac{1}{2\Delta }\left(\begin{array}{cc}{b}_i& 0\\ 0& c_i\\ c_i& b_i\end{array}\right)$

根据材料和平面应力单元的性质，矩阵D^P表示为：

$D^P=\frac{E}{1-\nu }\left(\begin{array}{ccc}1& \nu & 0\\ \nu & 1& 0\\ 0& 0& \frac{1-\nu }{2}\end{array}\right)$

根据有限元理论，平面应力单元ijm的厚度用t^e表示，则刚度矩阵表 示为

$K_{ijm}^P=\int \int {\left(B^P\right)}^T{D}^P{B}^P{t}^{e}dxdy$

对于二次壳体单元，其计算公式与平面应力单元相仿，其形状函数N^S表示为：

$N^S=\left(\begin{array}{c}{P}_1-P_4+P_6+2(P_7-P_9)\\ -b_j(P_9-P_6)-b_m{P}_7\\ -c_j(P_9-P_6)-b_m{P}_7\\ P_2-P_5+P_4+2(P_8-P_7)\\ -b_m(P_7-P_4)-b_i{P}_8\\ -c_m(P_7-P_4)-c_i{P}_8\\ P_3-P_6+P_5+2(P_9-P_8)\\ -b_i(P_8-P_5)-b_j{P}_9\\ -c_i(P_8-P_5)-c_j{P}_9\end{array}\right)$

$P^T=\left(\begin{array}{c}{L}_i\\ L_j\\ L_m\\ L_i{L}_j\\ L_j{L}_m\\ L_m{L}_i\\ L_i^2{L}_j+\frac{1}{2}{L}_i{L}_j{L}_m\left(3\right(1-{\mu }_m)L_i-(1+3{\mu }_m)L_j+(1+{\mu }_m\left)L_m\right)\\ L_j^2{L}_m+\frac{1}{2}{L}_i{L}_j{L}_m\left(3\right(1-{\mu }_i)L_j-(1+3{\mu }_i)L_m+(1+3{\mu }_i\left)L_i\right)\\ L_m^2{L}_i+\frac{1}{2}{L}_i{L}_j{L}_m\left(3\right(1-{\mu }_j)L_m-(1+3{\mu }_j)L_j+(1+3{\mu }_j\left)L_j\right)\end{array}\right)$

其中L_i,L_j,L_m表示为，

$L_i=\frac{a_i+b_{i}x+c_{i}y}{2\Delta },L_j=\frac{a_j+b_{j}x+c_{j}y}{2\Delta },L_m=\frac{a_m+b_{m}x+c_{m}y}{2\Delta }$

其中l_i,l_j,l_m表示边i,j,m的长，则上式中μ_i,μ_j,μ_m表示为：

${\mu }_i=\frac{l_m^2-l_j^2}{l_i},{\mu }_j=\frac{l_i^2-l_m^2}{l_j},{\mu }_m=\frac{l_j^2-l_i^2}{l_m}$

由此可以计算矩阵B^S，

$B^S=\left(\begin{array}{cc}\frac{\partial }{\partial x}& 0\\ 0& \frac{\partial }{\partial y}\\ \frac{\partial }{\partial y}& \frac{\partial }{\partial x}\end{array}\right)\left(\begin{array}{c}\frac{\partial }{\partial x}\\ \frac{\partial }{\partial y}\end{array}\right)N^S$

根据材料和二次壳体单元的性质，矩阵D^S表示为：

$D^S=\frac{E{\left(t^e\right)}^3}{12(1-{\nu }^2)}\left(\begin{array}{ccc}1& \nu & 0\\ \nu & 1& 0\\ 0& 0& \frac{1-\nu }{2}\end{array}\right)$

根据有限元理论，平面应力单元ijm的厚度用t^e表示，则刚度矩阵表 示为

$K_{ijm}^S=\int \int \int {\left(B^S\right)}^T{D}^S{B}^{S}dxdydz$

故在此，每个单元由二次壳体单元和线性平面应力单元组合而成，其 单元ijm的刚度矩阵K_ele表示为：

$K_{ele}=\left(\begin{array}{ccc}{K}_{ijm}^P& 0& 0\\ 0& K_{ijm}^S& 0\\ 0& 0& K_z\end{array}\right)$

其中K_z＝1e-8，使得整个方程可求解。此处不同的K_z值会显著的影响 有限元仿真的效果。过大的K_z会使得整体元过硬，过小的K_z会使得整体方 程不可解。

上述刚度矩阵是假设ijm在xy平面上，在三维实体应用时，需要求解 旋转矩阵T_ijm，则在xyz坐标系下，

$K_{ele}^{global}=T_{ijm}{K}_{ele}{T}_{ijm}^T$

根据有限元理论，将所有元素的组装起来得到模型整体刚度矩阵 K_sys。

结合以上公式，K_sys可以显式的表示为t^e的多项式。由于t^e与α,β存在 映射关系，即等价于K_sys可以显式的表示为α,β的多项式。

根据有限元方程，在外力F作用下，根据整体刚度矩阵K_sys可以求解 模型每个顶点的位移U：

K_sysU＝F

根据敏感度分析原理，对于在模型M上顶点i处的厚度t_i变化对整体位 移产生的影响可以表示为：

$\frac{\partial F}{\partial t_i}=\frac{\partial K_{sys}}{\partial t_i}U+K_{sys}\frac{\partial U}{\partial t_i}$

由于外力保持不变，不随着厚度t_i变化，上式可简化为：

$\frac{\partial U}{\partial t_i}=-K_{sys}^{-1}\frac{\partial K_{sys}}{\partial t_i}U$

结合vonMises计算公式，可以根据位移对厚度参数的变化求得

若不采用敏感度分析方法，上述问题是非线性约束优化问题，求解较 慢且难以得到收敛解。若不采用交替迭代优化过程，则由于分块区域和过 渡区域中变量权重差别较大，且约束过度会造成优化求解部分不稳定。

采用敏感度分析方法将优化问题转化为线性规划问题，易于求解加速 了求解过程，且在有限元仿真过程中，减少的计算，加快了求解速度。

在实现过程中，K_sys可以显式的表示为α,β的多项式的性质显著提高 了计算速度。K_sys被存储为α,β的多项式，在优化过程需要对每次迭代求 解出的α,β重新计算K_sys，则只需将α,β乘以相应系数即可。此外，在计算 时也可以快速得到精确解。此性质成立的关键在于选取前述的二次壳 体单元和线性平面应力单元组合。若不采用此种元，则在计算过程中无法 得到此性质，会严重影响计算速度。

故建立优化目标函数，通过交替迭代优化过程，分别优化α,β。

优化目标是在模型承载指定外力的条件下，保证材料不断裂，保持所 需要的打印材料最小并且保持生成曲面光滑。即可用如下形式表达为：

$\underset{{\mathrm{\Delta t}}_i^k}{\min }\underset{i=1}{\overset{n}{\Sigma }}(t_i^k+{\mathrm{\Delta t}}_i^k)s_i$

$\left(\begin{array}{cc}s.t.& {\sigma }_v^k+\frac{\partial {\sigma }_v}{\partial t}\left(t_i^k\right){\mathrm{\Delta t}}_i^k<{\sigma }_{\max },\forall v\in M,\end{array}\right)$

$t^{\min }<y_i^k+{\mathrm{\Delta t}}_i^k<t_i^{\max },i=\mathrm{1...}n,$

其中s_i表示顶点i邻接三角形面片面积和的控制曲面光滑程度， t^min由3D打印机打印最薄厚度确定，w_i表示其邻域顶点的权重，若与点i 邻接的顶点数为N_ring，σ_max由打印材料确定。则等价 于密闭模型的体积V_M。

为了实现快速稳定的求得收敛解，通过迭代优化过程求解，在第k次迭 代过程中，阶段1和阶段2的求解分别为：

阶段1：

$\underset{{\mathrm{\Delta t}}_i^k}{\min }\underset{i=1}{\overset{n}{\Sigma }}({\alpha }_i^k+{\mathrm{\Delta \alpha }}_i^k)s_i$

$\begin{array}{cc}s.t.& {\sigma }_v^k\end{array}+\frac{\partial {\sigma }_v}{\partial t}\left(t_i^k\right){\mathrm{\Delta t}}_i^k<{\sigma }_{\max },\forall v\in M,$

${\alpha }_i^{\min }<{\alpha }_i^k+{\mathrm{\Delta \alpha }}_i^k<{\alpha }_i^{\max },i=\mathrm{1...}{n}_s,$

在阶段1求解过程中，固定β值不变，求解线性优化问题的解为 α^k+Δα^k，通过求解拉普拉斯方程，得到β^k。

阶段2：

$\underset{{\mathrm{\Delta t}}_i^k}{\min }\underset{i=1}{\overset{n}{\Sigma }}({\beta }_i^k+{\mathrm{\Delta \beta }}_i^k)s_i$

$\left(\begin{array}{cc}s.t.& {\sigma }_v^k+\frac{\partial {\sigma }_v}{\partial t}\left(t_i^k\right){\mathrm{\Delta t}}_i^k<{\sigma }_{\max },\forall v\in M,\end{array}\right)$

$t^{\min }<{\beta }_i^k+{\mathrm{\Delta \beta }}_i^k<{\beta }_i^{\max },i=\mathrm{1...}{n}_t,$

在阶段2的求解后，得到β^k+Δβ^k，此处同样通过求解拉普拉斯方程， 求得α^k+1，作为k+1次迭代的初值。

阶段1和阶段2依次进行，循环求解求解，直到约束条件满足，且优 化目标V_M不在下降，则停止优化过程，得到最终结果t^F。

若不采用迭代优化求解部分，则求解过程不收敛，无法得到可行解。

根据求解得到的网格M上每个顶点的厚度参数根据积分曲线上的 位置，像两侧各处取得生成曲面位置，即得到封闭的三维模型。