CN105313336B - 一种薄壳体3d打印优化方法 - Google Patents

Abstract

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

Description

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

技术领域

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

背景技术

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

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

近几年有一些研究者专注研究如何增强打印物体的结构强度。STAVA和VANEK提出通过改善物体形状和增加支撑来增加物体的结构强度，同时减少打印消耗的材料，参见STAVA,O.,VANEK,J.,BENES,B.,CARR,N.,AND MˇECH,R.：Stress relief:Improvingstructural strength of 3d printable，566，2012，ACM Trans.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.,AND LIU,X.2013.Cost-effectiveprinting of 3d objects with skin-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)表示为

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

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

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

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

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

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

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

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

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

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

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

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

K_sysU＝F

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

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

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

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

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

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

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

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

阶段1：

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

阶段2：

在阶段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点法向方向，参见Peng Jianbo,Kristjansson,Daniel Zorin,Denis(2004).Interactive modeling of topologically complex geometric detail.ACMTransactions on Graphics(TOG)。

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

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

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

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

2)根据模型几何特征，采用fuzzy cut方法对模型进行分片，参照Katz Sagi，TalAyellet，2003，Hierarchical mesh decomposition using fuzzy clustering and cuts，ACM Transactions on Graphics(TOG)。分片结果如图3所示。结合图描述一下

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

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

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

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

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

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

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

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

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

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

其中，可以表示为：

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

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

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

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

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

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

其中L_i,L_j,L_m表示为，

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

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

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

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

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

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

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

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

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

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

K_sysU＝F

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

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

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

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

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

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

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

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

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

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

阶段1：

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

阶段2：

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

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

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

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

Claims (1)

1.一种薄壳体3D打印优化方法，其特征在于：包括以下步骤：

(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上每一点进行厚度修正，采用averageddistancefunctions函数计算网格模型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)=2min\{T_G^{-}\left(x\right),T_G^{+}\left(x\right),0.5T_{Max}\};$

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

其中，采用fuzzy cut方法对模型进行分片，具体为：

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

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

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

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

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

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

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

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

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

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

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

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

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

K_sysU＝F

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

$\frac{\&part;F}{\&part;t_i}=\frac{\&part;K_{sys}}{\&part;t_i}U+K_{sys}\frac{\&part;U}{\&part;t_i}$

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

$\frac{\&part;U}{\&part;t_i}=-K_{sys}^{-1}\frac{\&part;K_{sys}}{\&part;t_i}U$

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

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

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

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

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

$t^{min}<t_i^k+{\mathrm{\&Delta;t}}_i^k<t_i^{\max },i=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;}}({\mathrm{\&alpha;}}_i^k+{\mathrm{\&Delta;\&alpha;}}_i^k)s_i$

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

${\mathrm{\&alpha;}}_i^{min}<{\mathrm{\&alpha;}}_i^k+{\mathrm{\&Delta;\&alpha;}}_i^k<{\mathrm{\&alpha;}}_i^{max},i=1...n_s,$

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

阶段2：

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

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

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

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

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

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