冰雪颗粒3D打印的材料填充算法

陈贵亮， 贾镇涛， 杨冬， 李子浩， 杜文博

(河北工业大学机械工程学院， 天津 300401)

3D打印材料输出和粘接的方法依所处理的材料的物理属性而不同。目前常见的方法包括：光固化技术[1]、熔融沉积成型技术[2]、片材堆积制造技术、粉末粘接成型技术[3]、选择性激光烧结技术[4]、激光直接烧结金属技术等。所用材料涉及金属、陶瓷、塑料、凝胶、纸张、水泥等，应用范围包括大到建筑构件、航空航天大型腔体、支撑零件的制造，小到亚微米级的微机构零件的成型[5]，包括对生物材料的打印[6]。将上述成型方法中所使用的材料替换，也可实现其他方向的应用。文献[7]中将激光烧结技术中的金属粉末替换为石英或石膏粉末即可打印出类岩石样品。而目前成熟的3D打印技术并不适用于冰雪材料的打印。

冰雪运动是中国北方常见的运动与娱乐活动，其中的冰雪景观的建筑和场景的修补一直由人工实施，存在着劳动强度大、效率低和对人员技术要求高的问题。北京冬奥会的即将召开，必将激发中国冬季运动的发展，因此对于冰雪景观建筑和修补的需求也会扩大。为满足这种需求，秉承科技冬奥的精神，提出了针对冰雪材料的机器人3D打印技术的研究课题。目前，机器人技术与建筑材料成型技术相结合的机器人化建筑3D打印在建筑领域已有成功的应用实例，如用混凝土结合其他材料，打印建筑景观、复原文物展览[8]。文献[9]研究了利用机器人移动打印混泥土建筑构件的方法。

冰雪在松散情况下质地松软、强度低、易结块，受压缩时水分容易析出，可形成具有较高强度的块[10]，液态水流动性强，难于迅速凝结成型为足够强度的实体，但实体的冰雪块遇到水分容易相互粘接。文献[11]中研究了微小冰结构的3D打印成型方法，是在液氮氛围中使微米级液滴沾附凝结从而成型的方法；文献[12]也针对液态材料的微结构凝结成型方法进行了研究，但均是在微尺度上有效，在宏观尺度上这种方法不能实现。文献[13]中将雪经过喷嘴喷出，然后堆积成型，但是未压实的冰雪材料松散而机械强度差，不能堆积过多，否则难以形成理想的形状。文献[14]使用液态水在冷平面上冻结成冰，最终层层叠加成为三维实体，液态水的流动性较强，因而难以成型。从设备角度考虑，多数3D打印均采用结构封闭的框架结构。文献[15]设计了一种悬臂式3D打印机构，相对于传统打印机构，具有高灵活性，可以很好地保证打印精度。但是如用于户外，打印形体较大的作品，则需要较大的设备，显然不利于在户外多种环境下灵活的作业。

分析可知，实现户外环境的冰雪3D打印，在设备上应该采用可灵活移动、结构开放的机器人系统，在材料上使用制备的坚固冰雪颗粒材料通过排列，以水膜粘结，可具有较好的材料机械强度，能满足成型的效率需求。材料间用于粘接的水膜薄，在户外寒冷温度下容易凝结。以此形成冰雪雕塑的基本形体再辅以整体加强和外观修饰手段，可以成为一种较实用的户外场所大型冰雪雕塑3D打印技术。

轨迹规划是3D打印的重要问题，需根据材料特性研究适合的方法[16-17]，进行轨迹规划的前提则是对材料填充的规划计算。目前对于颗粒化冰雪填充的成熟方法尚鲜见报道，常见的线成面、面成体的填充方法不能直接应用于颗粒填充。针对这个问题，开展颗粒化冰雪的3D打印材料填充的快速规划计算研究。

1 冰雪颗粒3D打印系统

图1 系统方案Fig.1 System scheme

根据我们设计的冰雪颗粒3D打印方案，如图1所示，系统包含如下几个模块：移动平台、机器人、冰雪3D打印输出机构、水雾喷淋系统、冰雪颗粒材料的储存/输送系统、机器人控制系统，另外还有冰雪颗粒的制备成型和固化系统。其大致工艺过程是利用机器人进行颗粒坐标定位、打印机构接受输送系统的颗粒进行定周期的间歇打印，同周期内配合以水雾喷淋系统对颗粒粘结部位的间歇性水雾喷淋，作为颗粒的黏接剂。

整个工艺过程的实现，要求计算机能够预先对拟打印的数字模型根据颗粒的尺寸进行切片、逐层的颗粒排列填充。

2 填充模型和模型切片处理

由于是对等直径球粒的打印，不仅要求同层内球粒排列紧密有序，还要求相邻层间凸凹互补，才能保证整个形体的稳固致密，拥有最小的孔隙率。

2.1 填充模型

把三维模型中，填充的球体总体积与模型体积的比例称为密度。等直径球体堆积与一个体区域中能够达到的最大密度为74%。能实现密度为74%的堆积方法有两种：一是六方紧密堆积，其中相隔层的球位置正对，记为ABABAB；二是面心立方紧密堆积，其中每相隔两层的球位置正对，记为ABCABC。其中六方紧密堆积方式较为简单，因此，冰雪颗粒3D打印采用此堆积方式。依据此方法，设图2(a)为某一层的球体的紧密排列情况。任意球A共与6个球相邻接触，每3个球可围成一个空隙，其中一半是尖角向下的B空隙，另一半是尖角向上的C空隙。相邻的向上一层的球按照对应放置于尖角向下的B空隙之上，每个下层球可以同时与3个上层球相互接触，围成如图2(b)的四面体空隙。按此方法重复进行，层层累积，可以堆出图3的效果。

图2 颗粒两层分布情况Fig.2 Particle distribution in two layers

图3 颗粒堆积效果Fig.3 Particle accumulation effect

第一层网格直线表达式为

(1)

第二层网格直线表达式为

(2)

式中：i为网格水平线的取值范围；j为网格斜线的取值范围；yi1为第一层网格水平线纵坐标；与第一层网格斜线相交求得的球心横纵坐标为xj1和yj1；yi2为第二层网格水平线纵坐标；与第二层网格斜线相交求得的球心横纵坐标为xj2和yj2。

切片层厚度的计算公式为

(3)

2.2 模型数据切片处理

模型被存为STL格式后可由MATLAB软件的fread函数直接读取，读出的数据为表达包围实体模型的外表曲面的三角面片。按照切片的厚度，在模型的坐标系中设置一系列等距且于地面(x-o-y平面)平行的平面作为切片平面。逐个切平面检索与之相交的三角面片，并且求这些面片与该切平面的交点，按顺序连接这些交点成为若干个封闭的轮廓，此轮廓即形成了对应层切平面的填充轮廓。

2.2.1 模型数据的读取

为方便后面处理，将读出的三角面片数据存入一个数组矩阵tri_facet，矩阵的一行表示一个三角面片，由3个顶点坐标和面片平面的法相矢量共12个元素组成，如第i行表示第i个三角面片的参数，可表示为

tri_facet(i,:)=[x1i,y1i,z1i,x2i,y2i,z2i,x3i,y3i,z3i,nix,niy,niz]

(4)

矩阵的行数等于三角面片的总数目。

2.2.2 切片处理并建立轮廓

对三角面片数据进行切片并生成切片面填充轮廓需要：①对三角面片按坐标排序；②检索匹配切片平面和三角面片；③求各层切片平面与三角面片的交点；④连接交点为填充轮廓曲线。具体步骤如下。

步骤1对三角面片按坐标高低排序。提取tri_facet矩阵中每行数据的3个z向顶点坐标的最大值zimax和最小值zimin，将三角面片按zimin由小到大进行升序排序，若有三角面片的zimin相同，则对zimax进行比较，zimax值小的序号靠前。经过排序的面片存放于新的矩阵tri_facet_ord中，便于下一步的切平面匹配。

步骤2检索匹配切片平面和三角面片。这一步是为每一切平面匹配与之相交的三角面片，是一个检索过程。选择检索区域分割系数，经比较，最终确定采用黄金分割即0.618法速度较快。

根据切片层厚度和模型的高度确定切片层的数量和各层球心的高度(z值)，将其由低到高存数组lay_z中。令切片层最小序号min=1，最高层序号max可由size函数求得数组的列数来表示，即

max=size(lay_z,2)

(5)

分割检索区域的切平面层序号gold按式(6)计算，用floor函数将计算结果向负无穷大方向取整。分割检索区域的切平面高度值表达式见式(7)。

gold=min+floor[0.618(max-min)]

(6)

zseg=lay_z(gold)

(7)

以zseg与步骤1中排序后的三角面片矩阵tri_facet_ord中的数值做比较，可按式(8)判断三角面片与分割面的匹配关系。

(8)

式(8)中：zimin为第i个三角面片3个顶点坐标中z坐标的最小值；zimax为第i个三角面片3个顶点坐标中z坐标的最大值。

对于满足第i个面片与分割面相交情况，建立新的矩阵，将面片i记录为与gold层相交的面片；对于不相交、位于分割面上方情况，下一步对区间[gold+1,max]的层进行黄金分割检索；对于不相交、位于分割面下方情况，下一步对[min,gold-1]区间的切片层进行黄金分割检索。重复执行这个检索过程，直到切片平面检索序号矩阵中元素为空，即所有切片平面均匹配到与之相交的三角面片。

步骤3求各层切片平面与三角面片的交点。按顺序为每个切片平面层求解与之相交的三角面片的交点作为描绘轮廓的节点。

第i层切片平面的高度zi表达式为

(9)

步骤4连接交点为填充轮廓曲线。为避免节点过近或重复，需要设置一个阈值，两节点间距离小于该值即可认为二者重复，需删掉一个点。将每层删除冗余信息后的交点以此层原点为基准，根据三角面片的毗邻关系，按逆时针排序，并将相邻节点依次用线段相连，形成连续的平面轮廓线，在同一层内的不同轮廓曲线之间用NaN(非数值点)进行分隔，避免不同轮廓之间产生无效链接导致轮廓出现错误。

文献[18]中，对于相同问题采用二分法作为检索匹配方法，为比较其与上述方法的效率，针对各模型分别采用两种方法进行了对比，结果如表1所示。可以看出，在所有计算示例中黄金分割法所用时间(计算10次取平均值)均略低于二分法，这个趋势在三角面片数较多或者分层数较多时更明显一点。

表1 黄金分割和二分法进行匹配计算的时间比较Table 1 Time comparison of matching calculation between golden section and dichotomy

由于一层中有可能存在相互独立的多个轮廓岛或者相互为内外轮廓的复合轮廓，需要判别已确定层中材料填充的内部区域，具体可参考文献[19]，采用与之类似算法。经判断整理的区域数据按节点的x、y坐标值写为两个行矩阵分别为xv和yv。

3 切片的颗粒材料填充算法

多数3D打印采用线材输出，其填充算法是实现连续的直线或曲线在区域内的填充。采用颗粒材料实现冰雪3D打印，需要按照的等径球粒排列方法实现在区域内对点的阵列填充，而且相邻切片层的球粒需要互补错开排列。如果逐层按区域计算球粒中心的排放位置，不仅计算量较大，而且需要频繁进行边界判断，影响该工作效率。为简化计算，获得较好的填充规划效率和可靠性，提出利用平面排列网格模板进行区域筛选的规划方法。

设置两个可包容所有层轮廓的水平平面区域，分别在其内设置直线网格，求出网格交点作为区域填充的种子点，然后逐层与层内填充轮廓做比较，保留轮廓内的种子点作为填充点即可。

假设颗粒半径为10 mm，第一层网格的直线方程可写为式(10)，第二层为式(11)，其中i、j分别为与x轴平行直线和斜线的序号。两层的交点x、y坐标分别为x1、y1和x2、y2，分别表示奇数层和偶数层的种子点集合。

(10)

(11)

判断切片平面层中区域内的点。可以直接调用MATLAB中的inpolygon函数来判断网格交点是否在复合区域内。

当层数为奇数时，判断程序为

in1=inpolygon(x1,y1,xv,yv)；

当层数为偶数时，判断程序为

in2=inpolygon(x2,y2,xv,yv)；

返回值为1的点即为复合区域内的填充点。一层中所有点存于矩阵PointIn的对应行中，当对所有层均筛选完成后，即获得所有层的填充区域中的填充点，即完成了填充规划，以此为基础即可继续进行轨迹规划。

4 实例检验

为了验证本文算法的准确性，以图4所示零件对该切片筛点算法进行测试。模型的尺寸为600 mm(长)×400 mm(宽)×400 mm(高)，三角面片个数为952，分层厚度经计算为16.329 9 mm，对之进行处理。图5为第1层和第15层的球粒中心填充规划情况。可以看出，该算法能够精确筛选出多轮廓组成的复合区域内部的填充点。

图6为对各层填充点进行组合，形成的整个模型的球粒规划，该算法能够准确地以六方紧密堆积原理，将三维实体以点为单位进行离散填充规划。

图4 颗粒填充规划模型Fig.4 Particle filling programming model

图5 切平面颗粒填充效果Fig.5 Particle filling effect of the layer tangent plane

为进一步验证本文算法的可靠性，对图7展示的兔子与熊的模型进行填充计算。可以看出，本文算法对于具有复杂曲面和多重轮廓的实体进行颗粒填充规划也有效。其中，图7(a)模型尺寸为576.6 mm(长)×576.6 mm(宽)×1 185 mm(高)，填充点数量为 20 612，计算所需时间为 2.04 s；图7(b) 模型尺寸为1 873.2 mm(长)×1 585.5 mm(宽)×1 608.4 mm(高)，填充点数量为226 451，计算所需时间为15.18 s。该算法可以每秒筛选出超过104个填充点。

图6 三维填充效果Fig.6 Three-dimensional filling effect

图7 模型颗粒填充效果Fig.7 Model and particle filling effect

5 结论

根据冰雪材料3D打印机器人课题，针对冰雪材料的物理特性，提出了颗粒材料3D打印的新方法，对打印模型的数据处理及材料的颗粒填充算法进行了研究。得出如下结论。

(1)针对三角面片的排序和切片层平面与三角面片的匹配问题，提出采用黄金分割法作为区域分割搜索的方法，可以提高匹配速度。对于切片平面内颗粒材料区域填充的问题，确定采用六方紧密堆积方式，提出采用预设切片平面层的球粒填充模板，使之与切片平面轮廓进行比较的方法，保留其中轮廓内的点作为填充点，而逐层实现颗粒材料对三维实体的填充规划。这种方法简单直观，容易编程。通过实例计算可证明这种方法对于具有多重复合区域和具有复杂表面的模型的填充规划均有效。

(2)所研究的填充方法是根据颗粒材料的特点进行的填充规划，该工作为后续的路径规划奠定了基础。由于颗粒材料可以通过简单预制成为具有理想强度的材料，水分为粘接材料通常凝固速度较快，可以避免有些材料因为需要较长干透时间而影响打印效率的情况，该方法也可推广到其他材料的打印中。因此本文研究方法应该具有更广的应用前景。