从上个世纪80年代第一家3D打印设备公司诞生到现在, 随着材料与工艺技术的不断进步, 3D打印技术愈发成熟并逐渐流行起来. 方便个人定制、适用于任意形状、快速成型等突出优点使3D打印有力的推动了个性化产品的设计制造^[1]、机械结构设计^[2]、三维重建^[3]等领域的发展.

尽管3D打印机器越来越普及, 打印材料的价格依然远高于传统制造方式所用的原料, 如常用的光敏树脂材料价格高达2000元/千克. 这一定程度上制约了3D打印技术的发展和推广. 因此, 如何有效降低模型的实际打印体积, 节约打印成本是3D打印技术面临的一个重要挑战^[4,5]. 针对节材打印这一问题, 国内外研究者提出了许多解决方案. Wang等^[6]受到建筑中的桁架结构启发提出一种“蒙皮-钢架”的轻质结构, 并基于该结构进行模型的优化. 这种方法能在有效降低打印成本同时, 保证打印结果满足各种约束条件, 如物理强度、平衡性等. Lu等^[7]受到自然界中的一种介观结构——蜂窝结构的启发, 提出一种内部掏空优化方法. 该方案基于应力分布来对模型空间进行自适应Voronoi分割得到大小不同的Voronoi单元, 通过掏空Voronoi单元构造出类似于蜂窝的内部结构. 其优化结果也能够一定程度上降低打印体积, 同时有着很好的物理强度. Zhang等^[8]受到物体的中轴和骨骼结构的启发, 提出一种中轴树结构作为物体内部结构型. 这种结构型主要由中轴结构、边界框架以及一组连接杆件这三部分组成. 徐文鹏等^[9]借鉴传统渐进结构优化方法, 给出一种拓扑优化算法. Zhang与徐文鹏的方法都很好地反应了模型受力的传递路线并且大幅节省了材料的使用.

另外, 考虑到打印出的模型在现实场景下可能受到的各种外力, 在节省材料的同时, 也要保证模型结构的力学稳定性, 使得打印结果满足一定的强度要求. 这一方面可以通过采用具有更好力学性能的材料进行打印, 如具有点阵结构的新型复合材料^[10,11]; 另一方面则可以对模型做进行几何结构优化. 本文主要讨论后者. 在打印模型之前对模型进行准确的力学稳定性判断, 并对存在强度问题的部位进行合理的结构优化可以归结为模型的强度优化问题^[4,5]. 对此类问题 Stava等^[12]给出一种解决方案, 该方案通过自动检测薄弱区域并对模型进行内部挖洞、局部加厚以及加支撑这三种方式来提高模型的强度, 调整后的模型可以承受用户指定的载荷. Zhou等^[13]依据模型的几何特征和材料属性直接计算出模型的最不利载荷, 同时识别出模型上最容易被破坏的区域. Zhao等^[14]提出了一种使用厚度参数来构造满足给定力学稳定性约束的壳模型的方法. 该方法依据模型的边界面和用户指定的外力, 根据应力约束来优化厚度参数以挤出内外表面, 最终构造出满足约束的壳模型.

上述研究成果在节材打印、强度优化方面取得了很好的效果, 为本文提供了宝贵的经验和参考. 同时, 我们注意到, 节材打印的相关研究^[6–8]都采用了一些特定的结构型式来生成内部稀疏的空间结构, 以此来减小打印体积. 这样带来的一个问题是模型中会存在着大量的悬空结构, 如Lu等^[7]的蜂窝结构、Zhang等^[8]的中轴树结构. 这种模型需要使用无需支撑的打印技术如3DP、SLS或者使用可溶解的支撑材料才能打印完成, 一定程度上提高了打印的硬件门槛和操作难度. 而关于强度优化的相关研究^[12–14]均是通过对模型做变形来解决的, 虽然能有效的提高模型强度, 但会导致模型外观的变化, 影响使用体验.

针对当前研究中存在的上述两个问题, 本文提出一种基于应力分布的壳模型构造优化方案. 本方案首先构造出均匀厚度的初始壳模型, 然后进行基于应力分布的壳模型优化迭代, 得到具有极小体积的壳模型. 经过多组实验验证, 优化后模型的体积为原模型的17.2%~24.4%, 节材效果显著. 此外, 本文方案生成的壳模型还具有如下优点: 不改变模型的外观; 使用普通的桌面级打印设备即可打印, 如最常见的FDM打印机; 添加的支撑材料方便去除.

1 优化方案概述

给定一个待打印的实心三维模型M, 使用打印材料m, 指定载荷F, 屈服强度σ_y, 最小打印厚度t_min. 本文所要解决的问题是如何通过优化 $M$ 的几何形状来减小模型的打印体积. 同时, 需要满足如下约束条件, 模型的外观不变、在载荷的作用下不发生屈服, 以及模型的最小厚度t_m不低于最小打印厚度. 所以优化问题可以描述为:

$\left\{\begin{array}{l} \min Vol(M) \\ {\rm{s.t.}} \;\;\;{\sigma _{\max }} < {\sigma _y},\;\;{t_m} \ge {t_{\min }} \\ \end{array} \right.$

(1)

其中, 第一个约束条件为屈服强度约束条件, 用来保证模型可以承载指定的载荷而不发生破坏. 依据第四强度理论设置, 即当构件各部位的von Mises等效应力(见第2节)的最大值σ_max达到材料的屈服强度σ_y时, 材料就会开始屈服, 发生塑性形变, 外力撤销后无法恢复原来形状.

针对上述优化问题, 本文给出一种基于应力分布的壳模型构造优化方案, 概述如下: 首先构造与M对应的均匀厚度空壳模型M₀, 对M₀进行应力分析, 如果不满足屈服强度约束, 则均匀增加M₀的厚度, 直到满足约束为止; 然后基于M₀的应力分布情况, 自适应地调整M_i各位置的厚度得到M_i+1, 直到迭代结果M_k达到某个约束条件时停止, M_k–1即是具有极小打印体积的优化后的壳模型. 由于模型外观不变, 即模型的外表面R不变 (注: 外表面R围成实心模型M), 因此优化壳的厚度也就等价于优化内表面S的形状.

2 应力分析

应力描述了连续介质内部各部分之间相互作用力的强度, 在连续介质力学中将其定义为单位面积所承受的作用力. 本文的壳模型构造优化方案一方面借助应力分析结果指导迭代过程, 调整模型厚度, 另一方面也使用应力分布信息验证优化结果, 判断是否还可以继续优化.

由于绝大多数桌面级3D打印机仅支持塑料(PLA、ABS、尼龙等)作为打印材料, 因此本文只考虑塑性材料作为打印材料的情况. 对塑性材料的屈服情况的判断, 一般使用材料力学的第四强度理论^[15]. 该理论认为引起材料屈服破坏的主要因素是畸变能密度, 其强度条件见公式(2), 其中σ₁、σ₂、σ₃为主应力, σ_y为材料的屈服强度.

$\sqrt {\frac{1}{2}\left[{{({\sigma _1} - {\sigma _2})}^2} + {{({\sigma _2} - {\sigma _3})}^2} + {{({\sigma _3} - {\sigma _1})}^2}\right]} < {\sigma _y}$

(2)

von Mises是一种屈服准则, 表示在一定的变形条件下, 当受力物体内一点的等效应力达到某一值时, 该点就开始进入塑性状态, 产生屈服. 屈服准则的值即为von Mises等效应力, 计算公式如下:

${\sigma _v} = \sqrt {\frac{1}{2}\left[ {{{\left( {{\sigma _1} - {\sigma _2}} \right)}^2} + {{\left( {{\sigma _2} - {\sigma _3}} \right)}^2} + {{\left( {{\sigma _1} - {\sigma _3}} \right)}^2}} \right]} $

(3)

由式(2)、式(3)可知, von Mises等效应力遵循第四强度理论, 所以本文使用von Mises等效应力来度量模型各个部位的应力水平, 并通过与材料的屈服强度进行比较来判断模型是否发生屈服.

为了给构件必要的强度储备, 材料力学中引入了许用应力^[15]的概念, 即将材料的屈服强度除以一个大于1的安全因数n, 作为构件应力不可超过的值, 通常对于塑性材料n取1.1~2, 此时屈服强度条件约束为:

${\sigma _{\max }} \le \frac{{{\sigma _y}}}{n}$

(4)

同理, 本文取不大于1的安全系数s₁作为屈服强度的系数, 则引入安全系数的屈服强度约束条件为:

${\sigma _{\max }} \le {s_1}{\sigma _y}$

(5)

本文中使用PLA材质时将s₁设置为0.7, 使用PC-ABS材质时将s₁设置为1(为了与文献[7]中Lu方案进行对比). 若想给模型更高的强度储备, 可以将s₁设置为更小的值(一般不小于0.5). 此外, 为了消除打印机的打印误差带来的影响, 本文还引入了最小打印厚度对应的安全系数s₂, 实验中将s₂设为1.2. 至此, 引入了安全系数的优化问题表示为:

$\left\{\begin{array}{l} \min Vol(M) \\ {\rm{s.t.}}\;\;\;{\sigma _{\max }} \le {s_1}{\sigma _y},\;\;{t_m} \ge {s_2}{t_{\min }} \end{array}\right. $

(6)

为了实现对模型应力分布的高效准确计算, 本文借助机械设计软件SolidWorks的Simulation工具包进行基于有限元的静力分析. 依次进行模型的四面体化, 物理环境的配置(包括所受外力F, 材料属性m, 固定位置等), 以及基于有限元的静力分析. 最后将四面体化的结果和应力分布结果导出, 作为模型优化程序的输入. 图1展示了此应力分析步骤的工作流程. 需要说明的是, 虽然此工具能够得到精确程度较高的应力分布结果, 但计算效率不高, 本文中所有模型的应力分析时间在3~10 min/次.

3 壳模型的构造与优化 3.1 均匀厚度壳模型的构造

利用三维建模软件创建的待打印模型往往只有一个封闭的外表面R, 该外表面一般由一组三角面片连接组成. 因此如果不对模型进一步处理, 直接打印会得到一个实心的模型, 材料耗费会非常大. 所以一般在打印前需要对模型进行空壳处理, 即计算出一层内表面用来与外表面共同构成封闭的空壳模型. 对于凸多面体, 可以直接将外表面顶点沿着法线方向向内偏置一段距离来得到内表面. 但对于一般非凸多面体, 简单的向内偏置外表面会导致三角面片的自交. Jung等^[16]提出的自交检查和去除方法可以消除三角面片自交的情况, 但壳模型却不再厚度可控. 为此, 本文采用了基于模型体素化表示的等值面提取算法来构造内表面(见算法1).

算法1. 基于体素化的均匀厚度壳模型构造算法

输入: 模型外表面R, 指定厚度为t.

输出: 模型内表面S (与R围成均匀厚度壳模型(S, R)).

1. 以a为边长, 体素化R所在包围盒, 得到一组体素, 记为Voxel.

2. FOR每个在Voxel中的体素vox.

3. 判断vox相对R的位置(分为内部, 外部, 相交).

4. 以R为边界, Voxel中体素为离散元素构建有向距离场D_c.

5. 以0.25a为边长, 构建细粒度有向距离场D_f.

6. 取D_c中具有距离值t/a的一层体素Voxel_c.

7. FOR每个在Voxel_c中的体素vox_c.

8. 获取8个顶点在D_f中的距离.

9. 对Voxel_c应用Marching Cubes算法^[17]提取内表面S.

10. RETURN S.

图2展示了算法1的相关流程. 此算法不但可以避免三角面片自交, 还能够提取出过渡非常平滑的内表面, 见图2(c). 需要说明的是, 在算法1第2、3行的实现上, 本文采用了由粗到细的体素化方式进行加速, 大幅提升了算法效率. 为了得到足够精细的内表面, 需要保证体素尺寸足够小, 因此就需要处理数目庞大的体素, 以尺寸为10×20×15 cm³的Kitten模型为例, 若体素边长为1 mm, 则需要计算 3e⁹个体素的位置, 这涉及到大量的相交检测判断. 在 Intel Core i5-4590 @ 3.30 GHz 平台下, 直接的对细粒度的体素进行位置判断的平均耗时为1.9分钟. 为了减少时间消耗, 我们采用由粗到细的体素化方式进行加速. 具体做法是, 先进行粗粒度的体素化, 体素边长设为 4 mm. 对这些大体素进行位置判断, 删除在模型外部的大体素, 由模型内部的大体素细分出的小体素(边长1 mm)一定在模型内部, 所以只需对那些由模型边界上的大体素细分出的小体素进行位置判断, 从而大大减少了相交检测判断次数. 且所有在模型内部或边界上的边长1 mm的小体素的相对位置均被正确判断, 不会影响到有效信息的完整性. 使用这个新的方式对该模型进行处理, 平均耗时缩短为15秒, 为优化前的13.2%, 效率提升明显.

3.2 基于应力分布的壳模型构造与优化

为了保证模型打印完成之后具有一定的力学稳定性(即可以承受给定的载荷), 需要模型上各单元的von Mises等效应力不超过材料的屈服强度. 所以当模型简单处理为均匀厚度的壳模型时, 为保证所有单元的应力值在安全范围内, 对于那些本身应力值较低的区域需要和应力值高的区域具有相同厚度. 这会导致材料的浪费. 这一问题在应力分布不均的模型身上会突显出来. 故而对于算法1构造出的均匀厚度模型, 需要进一步去调整内表面形状, 最大限度的减小模型体积.

本文提出一种基于应力分布自适应调整内表面形状的方法. 具体做法是让顶点沿着法线方向外扩, 应力大的顶点具有较小外扩速度, 应力小的顶点具有较大外扩速度. 图3给出了基于此方法进行一次内表面顶点外扩的示意. 其中顶点a, b具有较小应力, 顶点c, d具有较大应力, 外扩后的形状如图3(c)所示.

基于以上思路, 我们给出每次迭代时顶点v的外扩距离计算公式, 见式(7). 其中, d_o为顶点v到模型外表面的距离, σ_v为v的应力值, s₁、s₂分别为屈服强度、最小打印厚度的安全系数(见第2节). 公式第一项用来保证最小打印厚度约束, 第二项用来计算顶点在当前应力情况下的外扩距离, 其中p用来控制每次外扩的最大距离.

$d = \min \left\{ {{d_o} - {s_2}{t_{\min }},{\rm{ }}p\max \left\{ {0,{\rm{ }}1 - \frac{{{\sigma _v}}}{{{s_1}{\sigma _y}}}} \right\}} \right\}$

(7)

基于式(7), 我们提出基于应力分布的壳模型构造与优化算法, 见算法2.

算法2. 基于应力分布的壳模型构造与优化算法

输入: 外表面R, 屈服强度σ_y, 最小打印厚度t_min.

输出: 壳模型(S, R), S为构造的内表面.

1. 对R进行体素化并构建距离场D_c、D_f.

2. 利用Marching Cubes算法按厚度t提取S₀.

3. WHILE (S₀, R)最大应力σ_max>s₁σ_y.

4. t←t+t₀.

5. 利用Marching Cubes算法按厚度t提取S₀.

6. WHILE (S_i, R)最大应力σ_max≤s₁σ_y.

7. FOR 每个S_i上的顶点v.

8. 利用式(7) 计算外扩距离d_v.

9. IF $\scriptstyle \forall $ d_v==0.

10. BREAK.

11. FOR 每个S_i上的顶点v

12. 沿顶点法线方向向外移动d_v距离.

13. 对S_i应用算法1按厚度0提取S_i+1.

14. i←i+1.

15. IF $\scriptstyle \forall $ d_v==0.

16. S←S_i.

17. ELSE

18. S←S_i–1.

19. RETURN (S, R).

算法2分为两个阶段: 前5行进行均匀厚度壳模型的迭代构造; 后面则是进行基于应力的壳模型的迭代优化. 此算法可以构造出具有极小体积的非均匀厚度壳模型, 从而大幅减少打印材料的成本.

注意这里得到的不是具有全局最小体积值的模型. 因为在每次调整内表面时, 各顶点移动方向均是向模型外部, 模型体积单调减小, 而载荷保持不变, 各单元的应力会逐渐增大, 直到某一单元的应力达到材料的屈服强度时为止. 若此时允许那些达到屈服强度约束的顶点向内移动, 使得其应力回到屈服强度约束之下, 则可以固定这些顶点, 继续向外扩展其他顶点, 这会进一步减少模型体积. 然而这种增加顶点移动自由度的做法会极大地增加算法的复杂程度, 考虑到每次迭代过程中应力分析耗时较大的特点, 如此进行下去会非常耗时, 可行性不强. 因此, 本文优化算法不允许顶点向内移动, 而是让应力较小的顶点具有较大的外扩速度, 好处是可以很快收敛生成具有极小体积的壳模型, 保证方案的可用性.

3.3 算法实例分析

这里以Kitten模型为例, 展示算法2的迭代过程与优化效果. 初始Kitten模型包围盒尺寸为6.5×10×6 cm³, 最小打印厚度为1 mm, 施加载荷为 1000 N, 这里使用的打印材料是PLA, 其相关物理属性见表1. 另外, 式(4) 中相关系数为: s₁=0.7, s₂=1.2, p=0.002.

首先构造均匀厚度壳模型, 初始厚度为5 mm. 然后对其进行应力分析(见图4). 上方箭头表示载荷, 底部箭头表示模型的固定位置. 此时最大应力为22.7 MPa, 小于屈服强度约束s₁σ_y (=49 MPa).

接着对初始壳模型进行优化迭代, 得到最终结果. 图5展示模型优化过程中剖面的变化. 结合图4中的应力分布信息, 可以发现优化后模型在初始应力较高的区域(方框标识), 厚度较大, 在应力较低的区域(圆圈标识), 厚度较小. 图6展示了优化过程中体积、应力的变化, 可见随着模型体积的降低, 模型的平均应力在逐渐增加, 直至模型的最大应力超过屈服强度约束. 经计算, 优化后Kitten模型的体积为输入模型的22.4%, 节材效果显著.

4 实验与讨论

本文算法由一个待打印的实心3D网格模型开始, 对其进行内壳模型的构造, 并基于应力分布进行内壳优化, 最终得到的内表面与输入的外表面围成体积极小的壳模型. 本文的实验环境为Windows操作系统, 处理器为Intel Core i5-4590 @ 3.30 GHz, 内存为8 GB. 本文算法2的壳模型优化部分需要迭代5~9次得到最终结果, 而每次迭代时应力分析过程是影响算法速度的瓶颈, 一般在3~10分钟/次. 本文进行的所有实验中优化时间最多为1.5小时, 平均为35分钟. 为了充分验证算法的有效性, 我们在对多个经典模型进行优化. 图7和表2展示了优化结果的信息. 从表2中模型的优化数据来看, 使用本文算法优化后壳模型的体积仅为输入模型的17.2%~24.4%, 大幅降低了材料的用量.

我们与同样基于应力分布的Lu等在文献[7]中提出的节材优化方案做了3组对比. 两个方案均使用PC-ABS作为打印材料, 屈服强度约束取41 MPa. 对比情况见表3和图8. 可以发现, Lu等方案的优化结果虽然具有较小的最大应力, 但模型体积为原来的26.0%~43.6%, 远高于本文的19.8%~24.4%. 这主要因为Lu等方案的优化结果受到蜂窝结构这种特定内部结构型式的限制, 无法进一步减小模型体积.

我们使用基于FDM打印技术的3DTALK MINI打印机打印制造了3组PLA材质的优化模型及剖面(见图9). 模型高度均为9~10 cm. 打印优化模型需要3~3.6小时, 打印剖面模型需要1.4~2小时. 在模型底部预留小孔, 方便打印完成之后使用工具取出打印支撑. 为了验证优化后模型的力学稳定性, 本文使用英斯特朗电子动静态疲劳试验机对打印出的模型进行压力测试, 见图10. 以Bimba模型为例, 观察机器数据采样系统给出的折线图10(b)可以发现, 随着压力的增加, 模型形变大小呈线性增长, 并在载荷达到1000 N时, 形变为1.2 mm. 这说明模型处在弹性形变状态, 没有发生屈服, 表明优化结果模型可以承受指定载荷而不被破坏.

5 结论与展望

本文针对3D打印领域中的材料节省问题提出了一种基于应力分布的壳模型构造优化算法. 该算法能够生成体积极小的非均匀厚度壳模型. 在初始壳模型的构造过程中, 先将模型体素化表示, 采用由粗到细的体素化方式加速体素相对位置的判断, 接着构建距离场, 利用Marching Cubes算法提取初始内壳; 在内壳的优化过程中, 基于应力分布的内壳外扩过程可以自适应降低模型各位置厚度, 最终得到具有极小体积的壳模型.

本文实验模型的优化结果体积为原来体积的17.2%~24.4%, 显著降低了材料消耗. 此外, 本文方法还具有不改变模型外观、内部支撑容易去除、适用于多种类型打印机的突出优点. 最后, 我们通过压力测试验证了模型的力学稳定性, 进一步验证了本文算法的有效性.

本文工作还存在很多改进之处. 首先, 基于有限元方法的应力分析方法虽然计算准确, 但效率偏低, 影响了算法的执行速度. 因此, 提高有限元分析的计算效率、减少应力分析所需的时间是让本文方法能够更加高效应用的关键. 其次, 本文方法仅支持能够在打印机中一次打印完成的模型, 对于那些超出打印机尺寸的模型,需要先进行模型的分割^[18,19], 然后在进行体积优化. 这涉及到多个部件模型组装后的应力分析以及协同优化, 很具有挑战性, 也是我们下一步要研究的主要内容. 再次, 本文的约束条件包括屈服应力、最小打印厚度, 以及外表面保持不变, 而在实际应用中, 可能会有其他的要求, 比如模型的平衡性, 即如何让模型以某个姿态稳定放置或稳定运动^[20,21]. 因此, 下一步考虑从这个角度对本文优化方案进行扩展.