Multi-Axis Support-Free Printing of Freeform Parts with Lattice Infill Structures
11 Multi-Axis Support-Free Printing of Freeform Parts with Lattice Infill Structures
Yamin Li, Kai Tang*, Dong He, Xiangyu Wang
Department of Mechanical and Aerospace Engineering, Hong Kong University of Science and Technology, Clear Water Bay, Kowloon, Hong Kong *Corresponding author. E-mail: [email protected] (Yamin Li), [email protected] (Kai Tang)
Abstract : In additive manufacturing, infill structures are commonly used to reduce the weight and cost of a solid part. Currently, most infill structure generation methods are based on the conventional 2.5-axis printing configuration, which, although able to satisfy the self-supporting condition on the infills, suffer from the well-known stair-case effect on the finished surface and the need of extensive support for overhang features. In this paper, based on the emerging continuous multi-axis printing configuration, we present a new lattice infill structure generation algorithm, which is able to achieve both the self-supporting condition for the infills and the support-free requirement at the boundary surface of the part. The algorithm critically relies on the use of three mutually orthogonal geodesic distance fields that are embedded in the tetrahedral mesh of the solid model. The intersection between the iso-geodesic distance surfaces of these three geodesic distance fields naturally forms the desired lattice of infill structure, while the density of the infills can be conveniently controlled by adjusting the iso-values. The lattice infill pattern in each curved slicing layer is trimmed to conform to an Eulerian graph so to generate a continuous printing path, which can effectively reduce the nozzle retractions during the printing process. In addition, to cater to the collision-free requirement and to improve the printing efficiency, we also propose a printing sequence optimization algorithm for determining a collision-free order of printing of the connected lattice infills, which seeks to reduce the air-move length of the nozzle. Ample experiments in both computer simulation and physical printing are performed, and the results give a preliminary confirmation of the advantages of our methodology.
Key words : additive manufacturing; support-free; self-supporting, lattice infill; multi-axis printing, FDM
1. Introduction
Additive manufacturing (AM) technologies have brought a significant change to manufacturing, which makes it possible to fabricate parts with extremely complex features that are otherwise impossible to be produced by traditional means such as machining. Fused deposition modeling (FDM) is one of the most commonly-used types of AM, owing to its low cost and simplicity, which builds a part layer by layer by extruding a molten filament [1] on the layers. Most FDM systems are of the 2.5-axis printing configuration, namely, the geometric model of the part is first sliced into a series of parallel planar layers and then the material is deposited layer by layer along a fixed direction (Z+). The biggest limitation of this configuration is that support structure is usually needed in order to prevent the collapse of material when printing overhang features, which not only causes the extra cost of the printing material and printing time but also subjects the part surface to potential damage when the support is being removed eventually. Moreover, due to the nature of parallel slicing, parts built under 2.5-axis configuration are inherently susceptible to the staircase effect and the consequent poor surface quality. Infill structures are commonly used in AM to reduce the weight and cost of a solid part. One critical requirement for the design of infill structures is that they must be self-supporting, as it would be extremely difficult or even impossible to remove any interior support after the part is printed. Dong et al. [2] and Tamburrino et al. [3] reviewed the properties and the modeling processes of lattice infill structures in additive manufacturing. Wu et al. [4] used the idea of adaptive rhombic grids to generate infill structures to satisfy the manufacturing requirements on both the overhang-angle and wall-thickness, and the generated infill structures exhibit improved properties of both high stiffness and static stability. Similarly, Lee et al. [5] proposed a method for generating support-free elliptic voids by constructing a Voronoi diagram of ellipses, which aims at not only avoiding the need of interior supports but also achieving better mechanical properties than Wu’s rhombic infill structures [4]. Kuipers et al. [6] proposed a new self-supporting infill structure called
CrossFill , which has the advantage that the extrusion printing paths are continuous and free of self-overlap. Wang et al. [7] presented a support-free hollowing algorithm based on the offsetting operation which can reduce more volume of material than Wu’s method [4]. Yang et al. [8] also proposed a hollow-to-fill algorithm based on the voxel model to guarantee the support-free property of inner surfaces for shape optimization. Similar voxel model-based hollowing methods can also be found in [9–12]. Gupta and Krishnamoorthy [13] developed a framework that can guarantee a sparse infill pattern according to a given arbitrary polygonal mesh; in addition, it guarantees the existence of a single, continuous and crossover-free tool path in each layer by using a novel Euler transformation. Additionally, there are also some infill structure generation methods based on topology optimization [14–16] with self-supporting as a constraint. Notwithstanding their richness, all the above infill structure generation methods are based on the conventional 2.5-axis printing configuration, which means that, although the generated infill structures in the interior of the part are self-supporting, exterior supports are still required when printing the surface of overhang features. To reduce the exterior support, many approaches have been proposed, although they are still limited to the 2.5-axis printing configuration. Hu et al. [17] proposed an orientation-driven shape optimizer which could considerably slim down the support by finding an optimal building direction. Zhou et al. [18] developed a tree-like support structure generation method based on the Lindenmayer system and an Octree. Vanek et al. [19] presented an optimization framework for the reduction of support, which tries to minimize the support material while providing sufficient support. Tricard et al. [20] proposed an interior rib-like support structure to build hollowed parts. Nevertheless, although in certain cases the support volume could be significantly reduced by these methods, they are not able to fully eliminate the need of exterior support, simply due to the nature of 2.5-axis printing configuration. As a viable solution to the shortcomings of 2.5-axis printing configuration, the 3+2-axis configuration allows a part to be printed with a finite number of different building directions, although for each building direction the printing configuration is still of the 2.5-axis type. Gao et al. [21] proposed a method to properly decompose a model into several sub-parts, which can then be printed consistently with different building directions with a much reduced amount of supporting material. Wei et al. [22] developed a skeleton-based algorithm for partitioning a model into the least number of sub-parts, which can (at least in theory) totally eliminate the need of support. Wu et al. [23] gave a support-effective volume decomposition algorithm that can minimize the surface area of regions with large overhangs. Based on the 3+2-aixs configuration, Bhatt et al. [24] developed an algorithm to print accurate thin-shell parts with no support. Similar research can also be found in [25–31]. Nonetheless, all these improvements are based on the 3+2-axis configuration, which become less effective for freeform parts with complex features and thus lack generality. The emerging continuous multi-axis printing configuration is perhaps the ultimate solution to the support-free requirement. Basically, on a multi-axis printing platform, we can not only slice the model into non-planar (i.e., curved) layers, but also continuously adjust the nozzle orientation to align with the layer normal, so to restrict the overhang angle below a threshold, thus achieving a total support-free printing (at least in theory). Dai et al. [32] proposed a curved layer decomposition method for multi-axis printing based on the voxel model. This method is considered to be general for printing freeform parts; however, it demands a huge computational cost due to the nature of voxelization. Xu et al. [33] recently presented a curved layer-decomposition algorithm, which first establishes a geodesic distance field on the part surface, then generates a set of closed iso-geodesic contours on the part surface, and finally fills these 3D contours into curved layers. However, one critical problem with their algorithm is that, as the curved layers are obtained by hole-filling the boundary loops, they are easy to intersect with each other when the slicing distance is small. Xu’s method uses the geodesic distance field as an intrinsic indicator to slice the part, which shows good generality and robustness. In terms of the calculation of geodesic field, Crane et al. [34] proposed an efficient and robust method for computing geodesics in a Riemannian manifold. Facilitated by a geodesic distance field embedded on the solid model, the part can be decomposed into a set of curved layers [35,36] that might be suitable for multi-axis printing. Additionally, Etienne et al. [37] proposed a curved layer generation method for printing a part on a 3-axis printers, which can effectively eliminate the staircase effect. Their method first uses planar planes to slice a deformed model of the part and then maps the planar layers back to the original model to obtain the corresponding curved layers. Although the above curved-layer decomposition methods based on continuous multi-axis printing configuration can significantly reduce or even, in most cases, completely eliminate the need of exterior support, they are all for printing the entire solid volume. Actually, to the best of our knowledge, due to the newness of continuous multi-axis printing configuration, so far there has no published reports on how to automatically design an infill structure and generate a printing path for an arbitrary freeform solid such that no support will be required for either the boundary or the infill structure. In this paper, under the continuous multi-axis printing configuration, we present a new methodology for automatically generating an infill structure as well as the accompanying multi-axis printing path for an arbitrary freeform part, which will be support-free for printing both the infill structure and the boundary surface of the part. The outline of the methodology is given next while its details will be presented in the ensuing sections. First, a new geodesic distance field (GDF) based curved-layer slicing algorithm is proposed. Rather than using only the iso-geodesic contours on the part’s surface as the boundary loops of the curved layers as did in [33], we compute the geodesic distance field directly inside the 3D volume of the part, which provides a more natural volume decomposition. Referring to Figure 1, for a three-dimensional manifold (i.e., a watertight a solid), from the base of the part, we can define locally parallel geodesics that will fill the entire manifold, and the iso-geodesic distance surfaces (IGDSs) of this 3D field naturally decompose the whole part. Let us call this geodesic distance field the 𝛾 -GDF. Because IGDSs are always perpendicular to the geodesics, the overhang angle at the boundary of any IGDS of 𝛾 -GDF will be significantly reduced, making it possible to print the part without any support. As IGDSs never intersect each other, the tangling issue of potential intersections between close-by curved layers as filled contours is now conveniently averted. Then, based on 𝛾 -GDF, a new lattice infill structure generation method is proposed, which will enable support-free printing for both the part boundary surface and the generated infills. Assuming that the printing always begins from the base (i.e., the bottom), a series of IGDSs of 𝛾 -GDF, called 𝛾 -IGDSs (as shown in Figure 1 (b)), can be constructed which naturally decompose the entire part. These IGDSs are exactly the sought curved layers on which lattice infills will be planned. To generate the interior lattice infills, two other types of GDF, called 𝛼 -GDF and 𝛽 -GDF respectively, are established, which satisfy the desirable orthogonal property – the geodesics of the three fields are mutually orthogonal to each other. Next, the other two clusters of IGDSs – i.e., the 𝛼 -IGDSs and the 𝛽 –IGDSs – are generated, and the three clusters of IGDSs (i.e., the 𝛾 -IGDSs, the 𝛼 -IGDSs, and the 𝛽 -IGDSs) are orthogonal to each other, as illustrated in Figure 1 (b). For each 𝛾 -IGDS, a lattice infill pattern is formed by the intersection lines (isolines) between this 𝛾 -IGDS and the 𝛼 -IGDSs and the 𝛽 -IGDSs. Because the three clusters of IGDSs are orthogonal to each other, the generated lattice pattern on each 𝛾 -IGDS is assured of self-supporting. In the rest of the paper, the symbols “ 𝛾 -”, “ 𝛼 -” and “ 𝛽 -”, or the superscripts 𝛾 , 𝛼 and 𝛽 will be used to differentiate similar elements regarding the three GDFs. Finally, to cater to the collision-free requirement and also to improve the printing efficiency, we propose a printing sequence optimization method that aims at effectively reducing the air-move length of the nozzle while upholding the collision-free constraint. The rest of the paper is organized as follows. In Section 2, the detailed algorithm of curved layer slicing and lattice infill structure generation is presented. Then Section 3 gives the details of the printing sequence optimization method and the printing path planning method for the lattice infill patterns. In Section 4, to validate the proposed methodology, we report the results of both computer simulation and physical printing experiments on several representative freeform parts. Finally, in Section 5, we conclude the paper and offer some discussions. (a) (b) Figure 1
Illustration of the infill structure generation method: (a) the three geodesic distance fields; (b) the lattice structure formed by the intersection of the three clusters of IGDSs.
2. Geodesic distance field-based slicing and infill generation
In this section, we present our GDF-based algorithms for curved layer slicing and generating a self-supporting lattice infill structure. Specifically, in Section 2.1, the details of computing the three mutually orthogonal 3D GDFs on a tetrahedral mesh are given. Then, in Section 2.2, the algorithms for curved layer slicing and the generation of a lattice infill pattern in each layer are presented. Finally, in Section 2.3, the algorithm for the generation of a skeleton tree of the connected lattice infill patterns is proposed, which will facilitate the printing sequence optimization to be presented in Section 3.
We assume that the given freeform model is represented by a tetrahedral mesh M ( V , E , F , T ), where V , E , F and T are the collections of vertices, edges, faces, and tetrahedrons respectively. First, a 3D GDF embedded on the tetrahedral mesh M is established, where the field value at any vertex is its geodesic distance to the specified bottom of the model, and this GDF is named as 𝛾 -GDF (see Figure 1 (a)). We apply the Crane’s heat method on the tetrahedral mesh to calculate 𝛾 -GDF by setting the bottom as the heat source. In their method, the heat flow equation 𝑢̇ = ∆𝑢 is firstly solved discretely for a fixed time t ; then, the gradient vector field X can be calculated by 𝑋 = −∇𝑢/|∇𝑢| ; and finally the GDF can be determined by solving the Poisson equation ∆𝜙 = ∇ ∙ 𝑋 . First, to solve the heat diffusion equation 𝑢̇ = ∆𝑢 on the tetrahedral mesh, it is rewritten in a discrete form as (𝐈 − 𝑡𝐕 −1 𝐋 𝑐 )𝒖 𝑡 = 𝒖 (1) where 𝐈 is the identity matrix, 𝒖 is the initial temperature field vector, 𝒖 𝑡 is the temperature field vector at moment t , 𝐕 ∈ ℝ 𝑛×𝑛 is a diagonal matrix containing the vertex volumes, and 𝐋 𝑐 ∈ ℝ 𝑛×𝑛 is the Laplacian matrix. The detailed values of matrix 𝐋 𝑐 can be found in [34]. For the initial temperature field vector 𝒖 , the temperatures at the bottom vertices are set as 1, while the temperatures at other vertices are set as 0. The appropriate time step t can be set as 𝑡 = ℎ , where h is the average length of edges [34]. Then, the temperature field 𝒖 𝑡 at moment t can be calculated by solving Eq. (1). For the k th tetrahedron 𝑇 𝑘 , the temperature scalar field inside its volume is defined as a piecewise linear function 𝑢 𝑘 (𝑥) = ∑ 𝜑 𝑖𝑘 (𝑥)𝑢 𝑖𝑘4𝑖=1 , with 𝜑 𝑖𝑘 being the piecewise linear basis function that is valued 1 at vertex 𝑣 𝑖 and 0 at all other vertices, and 𝑢 𝑖𝑘 being the temperature value at vertex 𝑣 𝑖 . The discrete temperature gradient inside the tetrahedron can then be expressed as ∇𝑢 𝑘 = ∑ 𝑢 𝑖𝑘 ∇𝜑 𝑖𝑘4𝑖=1 (2) It should be noted that ∇𝜑 𝑖𝑘 is simply the vector orthogonal to face 𝑓 𝑖𝑘 and opposite to vertex 𝑣 𝑖 in tetrahedron 𝑇 𝑘 , pointing towards vertex 𝑣 𝑖 and with a magnitude of |∇𝜑 𝑖𝑘 | = 𝑎𝑟𝑒𝑎(𝑓 𝑖𝑘 )3|𝑇 𝑘 | [38], where |𝑇 𝑘 | denotes the volume of tetrahedron 𝑇 𝑘 . The gradient vector field 𝒈 𝑘𝛾 of 𝛾 -GDF can be obtained by normalizing ∇𝑢 𝑘 , i.e., 𝒈 𝑘𝛾 =∇𝑢 𝑘 /|∇𝑢 𝑘 | . In order to generate the other two vector fields 𝒈 𝑘𝛼 and 𝒈 𝑘𝛽 which are orthogonal to 𝒈 𝑘𝛾 , a reference vector 𝒓 is introduced, which can usually be set as 𝒓 = ( ) . Then, 𝒈 𝑘𝛼 and 𝒈 𝑘𝛽 can be calculated by {𝒈 𝑘𝛼 = 𝒓×𝒈 𝑘𝛾 |𝒓×𝒈 𝑘𝛾 | 𝒈 𝑘𝛽 = 𝒈 𝑘𝛼 × 𝒈 𝑘𝛾 (3) wherein the vectors 𝒈 𝑘𝛾 , 𝒈 𝑘𝛼 and 𝒈 𝑘𝛽 are all mutually orthogonal with each other. The integrated divergence of the gradient field associated with vertex 𝑣 𝑖 can then be written as (𝐷𝑖𝑣𝒈)(𝑣 𝑖 ) = ∑ ∇𝜑 𝑖𝑘 ⋅ 𝒈|𝑇 𝑘 | 𝑇 𝑘 ∈𝑁(𝑖) (4) where N ( i ) is the set of vertices immediately adjacent to vertex 𝑣 𝑖 . Finally, the three GDFs, i.e., the 𝛾 -GDF 𝝓 𝛾 , the 𝛼 -GDF 𝝓 𝛼 , and the 𝛽 -GDF 𝝓 𝛽 for all the vertices can be obtained by solving the following discrete Poisson equation 𝐋 𝑐 𝝓 = 𝒃 (5) where 𝒃 is the divergence field vector of the gradient field. Once the three orthogonal GDFs embedded on the tetrahedral mesh M ( V , E , F , T ) are obtained, the curved layers and the lattice infill pattern in each layer can be constructed. The 𝛾 -GDF is used to slice the whole part, namely, the curved layers are generated by interpolating a number of 𝛾 -IGDSs, and each layer is sandwiched between two adjacent 𝛾 -IGDSs. Let 𝛹 𝛾 = {𝜙 , 𝜙 , . . . , 𝜙 𝑖𝛾 , . . . } be the set of sampling geodesic distances for 𝛾 -GDF. For each 𝜙 𝑖𝛾 , a 𝛾 -IGDS 𝑆 𝑖𝛾 can be defined. Similarly, let 𝛹 𝛼 = {𝜙 , 𝜙 , . . . , 𝜙 𝑗𝛼 , . . . } and 𝛹 𝛽 = {𝜙 , 𝜙 , . . . , 𝜙 𝑗𝛽 , . . . } be the sets of sampling geodesic distances for the 𝛼 -GDF and 𝛽 -GDF respectively, and the corresponding 𝛼 -IGDSs and 𝛽 -IGDSs can also be interpolated . The lattice infill pattern in 𝑆 𝑖𝛾 can be formed by the intersection lines (i.e., the 𝛼 -isolines) between 𝑆 𝑖𝛾 and the 𝛼 -IGDSs, as well as the intersection lines (i.e., the 𝛽 -isolines) between 𝑆 𝑖𝛾 and the 𝛽 -IGDSs. Because the 𝛾 -, 𝛼 -, and 𝛽 -IGDSs are all mutually orthogonal to each other, the generated lattice infill pattern in 𝑆 𝑖𝛾 is assured of self-supporting if the nozzle orientation is aligned with the normal direction of 𝑆 𝑖𝛾 . Figure 2 illustrates an lattice infill pattern 𝐺 𝑖 (𝑉 𝑖 , 𝐸 𝑖 ) at 𝑆 𝑖𝛾 , which is an undirected graph whose vertex set 𝑉 𝑖 contains the intersection vertices between the isolines and faces of the tetrahedral mesh, the intersection vertices between the 𝛼 -isolines and the 𝛽 -isolines, and the interpolation vertices in terms of 𝜙 𝑖𝛾 at the mesh boundary, and the edge set 𝐸 𝑖 is a collection of edges between the vertices. Algorithm 1 shows the pseudocodes for the generation of 𝐺 𝑖 (𝑉 𝑖 , 𝐸 𝑖 ) . Specifically, Steps 1-12 find all the intersection vertices between the mesh faces and the 𝛼 -isolines, wherein a function called GenGxVertex ( 𝐹 𝑘 , 𝜙 𝑖𝛾 , 𝜙 𝑗𝛼 , & 𝑉 𝛼 ) is used to judge whether there exists an intersection vertex 𝑉 𝛼 between the triangle face 𝐹 𝑘 and an 𝛼 -isoline in terms of 𝜙 𝑗𝛼 – it returns true if an intersection vertex is found and false otherwise. Figure 3 illustrates the intersection vertex between an 𝛼 -isoline in terms of 𝜙 𝑗𝛼 and a face of the tetrahedral mesh. Because the 𝛼 -isolines are on 𝑆 𝑖𝛾 which has a 𝛾 -geodesic distance of 𝜙 𝑖𝛾 , there must exist two interpolation vertices in terms of 𝜙 𝑖𝛾 at the two edges of the triangle face respectively; assuming these two vertices are P and Q at edge AB and AC respectively (see Figure 3), the following condition must be true: {(𝜙 𝑎𝛾 − 𝜙 𝑖𝛾 )(𝜙 𝑏𝛾 − 𝜙 𝑖𝛾 ) < 0(𝜙 𝑎𝛾 − 𝜙 𝑖𝛾 )(𝜙 𝑐𝛾 − 𝜙 𝑖𝛾 ) < 0 (6) where 𝜙 𝑎𝛾 , 𝜙 𝑏𝛾 and 𝜙 𝑐𝛾 are 𝛾 -geodesic distances at vertex A , B and C respectively. Then, the coordinates of vertex P and Q , i.e., 𝑣 𝑃 and 𝑣 𝑄 , as well as the corresponding 𝛼 -geodesic distances can be calculated by { 𝑣 𝑃 = (|𝜙 𝑎𝛾 − 𝜙 𝑖𝛾 |𝑣 𝐵 + |𝜙 𝑏𝛾 − 𝜙 𝑖𝛾 |𝑣 𝐴 )/|𝜙 𝑎𝛾 − 𝜙 𝑏𝛾 |𝑣 𝑄 = (|𝜙 𝑎𝛾 − 𝜙 𝑖𝛾 |𝑣 𝐶 + |𝜙 𝑐𝛾 − 𝜙 𝑖𝛾 |𝑣 𝐴 )/|𝜙 𝑎𝛾 − 𝜙 𝑐𝛾 |𝜙 𝑝𝛼 = (|𝜙 𝑎𝛾 − 𝜙 𝑖𝛾 |𝜙 𝑏𝛼 + |𝜙 𝑏𝛾 − 𝜙 𝑖𝛾 |𝜙 𝑎𝛼 )/|𝜙 𝑎𝛾 − 𝜙 𝑏𝛾 |𝜙 𝑞𝛼 = (|𝜙 𝑎𝛾 − 𝜙 𝑖𝛾 |𝜙 𝑐𝛼 + |𝜙 𝑐𝛾 − 𝜙 𝑖𝛾 |𝜙 𝑎𝛼 )/|𝜙 𝑎𝛾 − 𝜙 𝑐𝛾 | (7) where 𝑣 𝐴 , 𝑣 𝐵 and 𝑣 𝐶 are the coordinates of vertex A , B and C respectively, and 𝜙 𝑎𝛼 , 𝜙 𝑏𝛼 and 𝜙 𝑐𝛼 are the 𝛼 -geodesic distances at vertex A , B and C respectively. The intersection vertex between the triangle face and the 𝛼 -isoline in terms of 𝜙 𝑗𝛼 must be on edge PQ , and the following condition holds: (𝜙 𝑝𝛼 − 𝜙 𝑗𝛼 )(𝜙 𝑞𝛼 − 𝜙 𝑗𝛼 ) < 0 (8) Finally, the coordinate 𝑣 𝑁 of the intersection vertex N can be calculated by 𝑣 𝑁 = (|𝜙 𝑝𝛼 − 𝜙 𝑗𝛼 |𝑣 𝑄 + |𝜙 𝑞𝛼 − 𝜙 𝑗𝛼 |𝑣 𝑃 )/|𝜙 𝑝𝛼 − 𝜙 𝑞𝑎 | (9) Function GenGxVertex ( 𝐹 𝑘 , 𝜙 𝑖𝛾 , 𝜙 𝑗𝛼 , & 𝑉 𝛼 ) first uses Eq. (6) and Eq. (8) to judge whether there exists an intersection vertex at the triangle face 𝐹 𝑘 , and then, if it exists, uses Eq. (7) and Eq. (9) to calculate the coordinate of the intersection vertex 𝑉 𝛼 . After the calculation of all the intersection vertices on the 𝛼 -isoline in terms of 𝜙 𝑗𝛼 , the corresponding edges are constructed by traversing all the tetrahedrons, i.e., if two faces of a tetrahedron contain an 𝛼 -vertex respectively, an 𝛼 -edge is defined to connect these two 𝛼 -vertices (as stipulated in Steps 7-11 in Algorithm 1). Similarly, Steps 13-24 are used to find all the intersection vertices between the mesh faces and the 𝛽 -isolines, wherein the function GenGyVertex ( 𝐹 𝑘 , 𝜙 𝑖𝛾 , 𝜙 𝑗𝛽 , & 𝑉 𝛽 ) is used to judge whether there exists an intersection vertex 𝑉 𝛽 between the triangle face 𝐹 𝑘 and the 𝛽 -isoline in terms of 𝜙 𝑗𝛽 – it returns true if an intersection vertex is found and false otherwise, and the corresponding 𝛽 -edges are constructed by traversing all the tetrahedrons. 0 Figure 2
Illustration of the lattice infill pattern 𝐺 𝑖 (𝑁 𝑖 , 𝐸 𝑖 ) at 𝑆 𝑖𝛾 Figure 3
Illustration of the intersection vertex between an 𝛼 -isoline of 𝜙 𝑗𝛼 and a face of the tetrahedral mesh After the generation of all the 𝛼 -edges and the 𝛽 -edges, the intersection vertices between them can be obtained by traversing all the tetrahedrons. As shown in Figure 4, a tetrahedron may contain several 𝛼 -edges and 𝛽 -edges, and some of them may intersect with each other. Take as an example the intersection vertex P between the 𝛼 -edge AB and the 𝛽 -edge CD shown in Figure 4, the following condition must be held: {(𝜙 𝑎𝛼 − 𝜙 𝑐𝑑𝛼 )(𝜙 𝑏𝛼 − 𝜙 𝑐𝑑𝛼 ) < 0(𝜙 𝑐𝛽 − 𝜙 𝑎𝑏𝛽 )(𝜙 𝑑𝛽 − 𝜙 𝑎𝑏𝛽 ) < 0 (10) where 𝜙 𝑎𝛼 , 𝜙 𝑏𝛼 and 𝜙 𝑐𝑑𝛼 are the 𝛼 -geodesic distances of node A , node B , and the 𝛼 -edge CD respectively, while 𝜙 𝑐𝛽 , 𝜙 𝑑𝛽 and 𝜙 𝑎𝑏𝛽 are the 𝛽 -geodesic distances of node C , node D , and the 𝛽 -edge AB respectively. Then the coordinate 𝑣 𝑃 of the intersection vertex P can be calculated by: 𝑣 𝑃 = (|𝜙 𝑐𝑑𝛼 − 𝜙 𝑎𝛼 |𝑣 𝐵 + |𝜙 𝑐𝑑𝛼 − 𝜙 𝑏𝛼 |𝑣 𝐴 )/|𝜙 𝑎𝛼 − 𝜙 𝑏𝑎 | 𝑜𝑟. 𝑣 𝑃 = (|𝜙 𝑎𝑏𝛽 − 𝜙 𝑑𝛽 |𝑣 𝐶 + |𝜙 𝑎𝑏𝛽 − 𝜙 𝑐𝛽 |𝑣 𝐷 ) /|𝜙 𝑐𝛽 − 𝜙 𝑑𝛽 | (11) where 𝑣 𝐴 , 𝑣 𝐵 , 𝑣 𝐶 , and 𝑣 𝐷 are the coordinates of vertex A , B , C , and D , respectively. Inside the tetrahedron 𝑇 𝑘 , all the possible intersection points between the 𝛼 -edges and the 𝛽 -edges should be first calculated according to Eq. (10) and Eq. (11). Then, the corresponding edges will be segmented by these intersection vertices, and new edges should be defined according to these intersection vertices. The steps in Algorithm 1 from the 25 th row to the 29 th are used to find the intersection vertices between the 𝛼 -isolines and the 𝛽 -isolines by traversing all the tetrahedrons, where function GenGxyVertices ( 𝑇 𝑘 , & VL ) is used to calculate the intersection vertices VL in the tetrahedron 𝑇 𝑘 (it returns true if an intersection vertex is found, and false otherwise). Figure 4
Illustration of the intersection vertices between the 𝛼 -edges and 𝛽 -edges inside a tetrahedron In Algorithm 1, the steps from the 30 th to the 40 th row are used to calculate the interpolation vertices in terms of 𝜙 𝑖𝛾 at the boundary of the tetrahedral mesh, as well as to define the corresponding boundary edges for 𝐺 𝑖 (𝑁 𝑖 , 𝐸 𝑖 ) . The interpolation vertices can be obtained by traversing all the boundary edges. Specifially, function GenBoundVertex ( 𝐸 𝑘 , & 𝑉 ) is used to judge whether there is an interpolation vertex V in terms of 𝜙 𝑖𝛾 at the boundary edge 𝐸 𝑘 , which returns true if an interpolation vertex is found, and false otherwise. Take the interpolation vertex A at the boundary edge MN shown in Figure 5 as an example, the following condition must be satisfied: (𝜙 𝑚𝛾 − 𝜙 𝑖𝛾 )(𝜙 𝑛𝛾 − 𝜙 𝑖𝛾 ) < 0 (12) where 𝜙 𝑚𝛾 and 𝜙 𝑛𝛾 are the 𝛾 -geodesic distances at vertex M and N respectively, and the coordinate 𝑣 𝐴 of the interpolation vertex A can be calculated by 𝑣 𝐴 = (|𝜙 𝑚𝛾 − 𝜙 𝑖𝛾 |𝑣 𝑁 + |𝜙 𝑛𝛾 − 𝜙 𝑖𝛾 |𝑣 𝑀 )/|𝜙 𝑚𝛼 − 𝜙 𝑛𝑎 | (13) 2 where 𝑣 𝑀 and 𝑣 𝑁 are the coordinates of vertex M and N respectively. After the calculation of all the boundary vertices, the corresponding boundary edges can also be established by traversing all the boundary faces. Take the boundary face MNL shown in Figure 5 as an example. Vertex A and B are the interpolation vertices in terms of 𝜙 𝑖𝛾 at edge MN and ML respectively, with which the edge AB can be defined. However, if there are some intersection vertices between face MNL and the 𝛼 -edges or the 𝛽 -edges, such as point C and E shown in Figure 5, edge AB will be broken by these intersection vertices and new edges (i.e., edge AE , EC and CB ) should be defined between vertex A and B . The above procedure to define the boundary edges in 𝐺 𝑖 (𝑁 𝑖 , 𝐸 𝑖 ) is realized by the steps in the 35 th to the 40 th row in Algorithm 1, where function
GenBoundEdges ( 𝐹 𝑘 , & EL ) is used to find all the boundary edges at face 𝐹 𝑘 . Finally, all the edges inside the tetrahedrons are saved in the edge list 𝐸 𝑖 of 𝐺 𝑖 (𝑉 𝑖 , 𝐸 𝑖 ) (i.e., by the steps in the 41 st to the 43 rd row in Algorithm 1) . Figure 5
Illustration of the interpolation points of 𝜙 𝑖𝛾 at the boundary edges Algorithm 1 Generation of the lattice infill pattern 𝐺 𝑖 at 𝑆 𝑖𝛾 Input: the tetrahedral mesh M ( V , E , F , T ) of the part with three embedded GDFs ( 𝜙 𝛾 , 𝜙 𝛼 and 𝜙 𝛽 ) , 𝜙 𝑖𝛾 , 𝛹 𝛼 = {𝜙 , 𝜙 , . . . , 𝜙 𝑖𝛼 , . . . } and 𝛹 𝛽 = {𝜙 , 𝜙 , . . . , 𝜙 𝑖𝛽 , . . . } Output: the lattice infill pattern 𝐺 𝑖 (𝑉 𝑖 , 𝐸 𝑖 ) at 𝛾 -IGDS 𝑆 𝑖𝛾 for each 𝜙 𝑗𝛼 in 𝛹 𝛼 for each face 𝐹 𝑘 in F if GenGxVertex ( 𝐹 𝑘 , 𝜙 𝑖𝛾 , 𝜙 𝑗𝛼 , & 𝑉 𝛼 ) 4 Put 𝑉 𝛼 into 𝑉 𝑖 , 𝐹 𝑘 𝑉 𝛼 end end for each tetrahedron 𝑇 𝑘 in T if two faces of 𝑇 𝑘 contain a node 𝑉 and 𝑉 respectively
9 Build an 𝛼 -edge 𝐸 𝛼 of 𝑉 and 𝑉 , 𝑇 𝑘 ← 𝐸 𝛼 end end end for each 𝜙 𝑖𝛽 in 𝛹 𝛽 for each face 𝐹 𝑘 in F if GenGyVertex ( 𝐹 𝑘 , 𝜙 𝑖𝛾 , 𝜙 𝑖𝛽 , & 𝑉 𝛽 ) Put 𝑉 𝛽 into 𝑉 𝑖 , 𝐹 𝑘 𝑉 𝛽 end end for each tetrahedron 𝑇 𝑘 in T if two faces of 𝑇 𝑘 contain a nodes 𝑉 and 𝑉 respectively 21 Build a 𝛽 -edge 𝐸 𝛽 of 𝑉 and 𝑉 , 𝑇 𝑘 ← 𝐸 𝛽 end end end for each tetrahedron 𝑇 𝑘 in T if GenGxyVertices ( 𝑇 𝑘 , & VL ) 27 Put all the vertices in VL into 𝑉 𝑖 end end for each edge 𝐸 𝑘 in E if GenBoundVertex ( 𝐸 𝑘 , & 𝑉 ) .and. 𝐸 𝑘 is a boundary edge Put V into 𝑉 𝑖 , 𝐸 𝑘 ← V end end for each face 𝐹 𝑘 in F if 𝐹 𝑘 is a boundary face GenBoundEdges ( 𝐹 𝑘 , & EL ) Put all the edges in EL into 𝐸 𝑖 end end for each tetrahedron 𝑇 𝑘 in T Put all the edges inside 𝑇 𝑘 into 𝐸 𝑖 end Take the Y model shown in Figure 6 (a) as an example. The tetrahedral model contains 4014 vertices and 17260 tetrahedrons. The vertices at the bottom are set as the heat source to calculate the 𝛾 -GDF, as well as the 𝛼 -GDF and the 𝛽 -GDF. The maximum 𝛾 -geodesic distance is 44.29 mm, and the whole part is sliced into 44 layers with the 𝛾 -geodesic interval set to be 1 mm ( 𝛹 𝛾 = {1,2,3. . .44} ). At each layer, the lattice infill pattern is generated by Algorithm 1 (the maximum 𝛼 -geodesic distance and 𝛽 -geodesic distance are 24.17 mm and 10.11 mm respectively, and the geodesic intervals for these two fields are all set as 2 mm, i.e., 𝛹 𝛼 = {2,4,6. . .24} , 𝛹 𝛽 = {2,4,6. . .10} ). The generated lattice infill structures are shown in Figure 6 (b-c), and Figure 6 (d) shows the lattice infill pattern at the 24 th layer. It can be seen that the lattice infill pattern in each layer is self-supporting because both the 𝛼 -GDF and the 𝛽 -GDF are orthogonal to the 𝛾 -GDF. 4 (a) (b) (c) (d) Figure 6
Generation of lattice infill structures for the Y model: (a) tetrahedral mesh of the model; (b) the generated curved layers; (c) lattice infill structures from the 1 st to the 24 th layer; (d) lattice infill pattern in the 24 th layer. Referring to Figure 7, the overhang angle 𝜃 at the part boundary is defined as the angle between the nozzle orientation and the normal direction n of the boundary surface, which measures the degree of danger of material falling. When angle 𝜃 is larger than a threshold (e.g., 135°), external support structure will be required in order to avoid the falling of material. Figure 8 (a) and (b) show the overhang angles of the Y model under the traditional 2.5-axis planar slicing configuration and our multi-axis GDF based curved-layer slicing method, respectively. As clearly seen, under our method, by aligning the nozzle orientation with the gradient direction of 𝛾 -GDF, angle 𝜃 is kept at or near 90 o , thus avoiding the need of extra support. 5 (a) (b) Figure 7
The overhang angle: (a) traditional 2.5-axis planar slicing; (b) multi-axis curved layer slicing. (a) (b)
Figure 8
Overhang angle of the Y model: (a) traditional 2.5-axis slicing method; (b) our multi-axis curved layer slicing method. The generated lattice infill pattern in each layer resembles an undirected graph that may contain several connected sub-graphs, which can be identified by using the DFS (depth first search) algorithm. For all the connected sub-graphs, we define a tree data structure simply called a skeleton tree that identifies the topological relationship between them. As illustrated in Figure 9, on the skeleton tree, each node represents a connected sub-graph, and every pair of adjacent sub-graphs are corresponded by an edge between the two representative nodes on the tree. For any node 𝑎 on the skeleton tree, if node 𝑏 is connected to 𝑎 by an edge on the tree, we call 𝑏 an upper-node of 𝑎 if 𝑏 ’s 𝛾 -geodesic distance is larger than 𝑎 ’s, and a lower-node otherwise. For example, on the tree in Figure 9, node 𝐺 has two upper-nodes, i.e., 𝐺 and 𝐺 , and only one lower-node, 𝐺 . The skeleton tree is constructed from the bottom towards the top, and the corresponding algorithm is given in Algorithm 1, wherein function AreTwoGraphsAdjacent ( 𝐺 𝑖,𝑗 , 𝐺 𝑖+1,𝑘 ) is used to judge whether any two nodes 𝐺 𝑖,𝑗 and 𝐺 𝑖+1,𝑘 are adjacent to each other (i.e., to be connected by an edge). Referring to Figure 10, to judge whether 𝐺 𝑖,𝑗 and 𝐺 𝑖+1,𝑘 are adjacent, we can first randomly select an edge on the triangular mesh (the boundary of the part) that intersects the boundary of 𝐺 𝑖,𝑗 ; then, from this edge, we trace out a geodesically steepest ascending path on the triangular mesh. Similarly, we can also trace out a geodesically steepest descending path from 𝐺 𝑖+1,𝑘 . The two nodes are adjacent to each other if at least one of the two paths go through both. Figure 9
Generating the skeleton tree of the connected lattice infill patterns Algorithm 2 Construction of a skeleton tree
Input:
All the connected sub-graphs of the decomposed layers: 𝐺 𝑖 {𝐺 𝑖,1 , 𝐺 𝑖,2 , . . . , 𝐺 𝑖,𝑗 , . . . } , i = 1, 2, …, n . Output:
Skeleton tree of the connected lattice infill patterns 1 /* n = the number of layers */ for i = 1 : n -1 3 for each graph 𝐺 𝑖,𝑗 of i th layer do for each IGDS 𝐺 𝑖+1,𝑘 of ( i +1)th layer do if AreTwoGraphsAdjacent ( 𝐺 𝑖,𝑗 , 𝐺 𝑖+1,𝑘 ) then
6 Construct an edge between the two graph nodes 7 end end end end Figure 10
Judge whether the two sub-graphs are adjacent to each other
3. Printing process planning
We have now decomposed the whole part into a number of curved layers according to the 𝛾 -IGDS and constructed the lattice infill pattern in each layer, as well as a skeleton tree that identifies the topological relationship of the connected sub-graphs of the lattice infill patterns. In this section, we will first present our printing sequence generation method (Section 3.1) and then describe how the printing paths are planned (Section 3.2). The strict increasing geodesic order of the skeleton tree has already defined a partial order of printing – any node must be printed before its upper-node(s). However, since printing is a time-continuous process, we must convert this partial ordering into a total ordering of traversal of the nodes, i.e., a single sequence of nodes to print, called a printing sequence. Hereafter we will interchangeably use the terms a “node” and a connected sub-graph (of a layer). The following criterion must be satisfied for any valid printing sequence:
Criterion 1 : a sub-graph can only be printed if its lower sub-graph(s) have already been printed. 8
Generally, there are two traversal strategies of a skeleton tree, i.e. the layer priority traversal (
LPT ) and the depth priority traversal (
DPT ). Refer to Figure 11, the layer priority traversal strategy traverses the skeleton tree layer by layer from bottom-up, which tends to avoid the collision to the utmost, for that the sub-graphs of each layer share the same 𝛾 -geodesic distance. On the other hand, the depth priority traversal strategy traverses the skeleton tree along branches in priority, which favours minimizing the air-movement of the nozzle. However, as shown in Figure 12, under DPT, if a branch grows too deep, it may cause collisions when printing other branches. (a) (b) Figure 11
Traversal strategies of a skeleton tree: (a) layer priority traversal; (b) depth priority traversal.
Figure 12
Illustration of possible collisions during a depth priority traversal
In order to reduce the air-move path length while ensuring no collisions, we propose an optimization strategy which seeks a compromise between LPT and DPT. First of all, the collision check between the nozzle and the sub-graphs must be modelled. In this paper, the shape of the nozzle is simplified by its bounding cone, as shown in Figure 13 (a). Admittedly, this simplification is too conservative; however, because collision check is not the main topic of this paper, we opt for this simplification to implement our algorithm. When the nozzle cone sweeps along the boundary curve of 9 a sub-graph with its orientation coincident with the surface normal direction n , the envelope of motion will be a ring-like ruled surface 𝑆(𝑢, 𝑣) = 𝑷(𝑢) + 𝒗𝒌(𝑢) , where
𝑷(𝑢) is an arbitrary point on the boundary curve of the sub-graph, and the unit vector 𝒌(𝑢) of the generator can be obtained by rotating the normal vector n at 𝑷(𝑢) around the tangent vector 𝝉(𝑢) of the boundary curve with the nozzle angle α , as shown in Figure 13 (a). Specifically, to construct the triangular mesh of the ruled surface, we first place a few sample points on the boundary curve of the sub-graph, then calculate their generators, and finally connect the generators as triangles. The upper and lower holes of the ring-like ruled surface should be filled to approximate the envelope volume of the cone over the entire sub-graph. In this paper, we adopt the advancing front mesh (AFM) technique [39] to fill the holes, which is robust and simple. To determine whether there is a potential collision when printing a sub-graph, we only need to check whether there are intersections between other sub-graphs and this envelope volume. For each sub-graph, we can calculate all the potential collision sub-graphs (PCG) (i.e., other sub-graphs that intersect the envelope volume of this sub-graph). Take the part shown in Figure 13 (b) as an example, the potential collision sub-graphs for sub-graph 𝐺 will be 𝐺 , 𝐺 , 𝐺 , 𝐺 and 𝐺 . The detailed procedures for calculating the PCGs of each sub-graph are given in Algorithm 3, where function
CollisionCheck ( G i , G j ) is used to judge whether there are intersections between sub-graph G j and the envelope volume of sub-graph G i – it returns true if an intersection is identified and false otherwise. (a) (b) Figure 13
Illustration of collision check between the nozzle and a sub-graph: (a) ring-like surface generated by sweeping the nozzle along the boundary of the sib-graph; (b) envelope volume of the nozzle motion. Algorithm 3 Calculation of the potential collision sub-graphs of each sub-graph
Input: the graph list { G , G , …, G i , …} and the nozzle cone angle α Output:
PCGs list for each graph 1 integer k = total number of sub-graphs 2 for i = 1: k do for j = 1: k do if i != j then if CollisionCheck ( G i , G j ) then Put G j into the PCGs list of G i . 7 end end end end Facilitated by the PCGs of each sub-graph, we propose a greedy traversal ( GT ) algorithm that can generate a collision-free printing sequence with a shorter air-move path length than that of the layer priority traversal. Besides Criterion 1, another criterion must be satisfied during the traversal: Criterion 2 : a sub-graph can only be printed if the PCGs of all the unprinted sub-graphs exclude this sub-graph. Specifically, in Algorithm 4 below, function
UpdateNodeCadidates ( ST ) is used to find all the printable candidate nodes which satisfy both Criterion 1 and Criterion 2 according to the current skeleton tree ST . And function SelectANode ( NC , N c ) is used to update the current node N c from the candidate node list NC . It will select the upper-node(s) of N c in priority. However, if the upper-node(s) are not included in NC , the node which is nearest to N c will be selected. To summarize, we traverse the skeleton tree along the branches in priority unless a potential collision is encountered. Algorithm 4 Printing sequence optimization algorithm
Input:
The skeleton tree ST of all the sub-graphs Output:
Printing sequence list PQ of the sub-graphs 1 Node candidates list NC = UpdateNodeCadidates ( ST ) 2 Current node N c = either node in NC while NC != ∅ do Label N c as printed Put N c into PQ NC = UpdateNodeCadidates ( ST ) N c = SelectANode ( NC , N c ) end For any node (sub-graph) in the skeleton tree, a traversal printing path needs to be determined. As illustrated in Figure 14 (a), for a connected sub-graph
𝐺(𝑉, 𝐸) , the number of intersection vertices ( 𝑣 , 𝑣 , . . . , 𝑣 𝑘 , . .. ) between the boundary curve and the isolines is always even and the degrees of these vertices are all 3, while the degrees of other vertices are either 2 or 4. To avoid excessive tool retractions, the sub-graph can be transformed into an Eulerian graph by properly trimming the boundary curve. As shown in Figure 14 (b), the boundary edges between the intersection vertex 𝑣 𝑘 and 𝑣 𝑘+1 ( k = 2, 4, 6, …) are deleted, so that the degrees of all the vertices become even and the new graph must contain an Eulerian tour. In this paper, the well-known Fleury’s algorithm is used to find an Eulerian tour in a connected graph. However, to avoid crossovers on the printing path, the printing path will turn at any vertex whenever its degree is 4, i.e., at the intersection vertices between the 𝛼 -isolines and the 𝛽 -isolines, as shown in Figure 15. To provide supports at places where the boundary edges are deleted, a support perimeter around the trimmed graph is added, which though is jagged at these places, as shown in Figure 14 (c) and (d). The offset distance between the boundary curve and the support perimeter is set to be the path width w , and the tooth length l can be adjusted to a proper value, e.g., l = 2 w . The nozzle orientation at each vertex is set to be coincident with the layer surface normal vector of the corresponding 𝛾 -IGDS. Because IGDSs are always perpendicular to the geodesics, the gradient vector of the geodesic field at the vertex can be directly used as the nozzle orientation. Figure 14
Printing path for a lattice infill pattern: (a) lattice infill pattern; (b) trimmed pattern that contains an Eulerian tour; (c) printing path for the lattice infill pattern; (d) generation of the support perimeter. Figure 15
Avoidance of crossovers
Due to the nature of curved layer slicing, the layer thickness h will no longer be a constant along the printing path. Refer to Figure 16, the intersection of the extruded filament is simplified to be a rectangle of h × w , where w is the path width; then, the following mass conservation equation should be satisfied during the printing process: 𝜋𝑟 𝑚2 𝑓 𝑚 = 𝜇𝑤ℎ𝑓 𝑝 (14) where r m is the radius of the original filament, f m is the feed rate of the filament, f p is the feed rate of the nozzle, and 𝜇 is a correction coefficient which is smaller than 1 and can be determined by experiments. The layer thickness ℎ 𝑘 at vertex 𝑣 𝑘 of the i th sub-graph 𝐺 𝑖 (𝑉 𝑖 , 𝐸 𝑖 ) can be calculated by finding the shortest distance of this vertex to the previous ( i- 𝐺 𝑖−1 (𝑉 𝑖−1 , 𝐸 𝑖−1 ) : {ℎ 𝑘 = 𝐹𝑖𝑛𝑑𝑆ℎ𝑜𝑟𝑡𝑒𝑠𝑡𝐷𝑖𝑠𝑡𝑎𝑛𝑐𝑒(𝑣 𝑘 , 𝐺 𝑖−1 (𝑉 𝑖−1 , 𝐸 𝑖−1 ))𝑖𝑓 ℎ 𝑘 ≥ 𝜆𝜙 𝑖 , ℎ 𝑘 = 𝜆𝜙 𝑖 (15) wherein function 𝐹𝑖𝑛𝑑𝑆ℎ𝑜𝑟𝑡𝑒𝑠𝑡𝐷𝑖𝑠𝑡𝑎𝑛𝑐𝑒(𝑣 𝑘 , 𝐺 𝑖−1 (𝑉 𝑖−1 , 𝐸 𝑖−1 )) traverses all the edges of 𝐺 𝑖−1 (𝑉 𝑖−1 , 𝐸 𝑖−1 ) to find the shortest distance. To avoid a too large layer thickness, ℎ 𝑘 should not exceed a threshold 𝜆𝜙 𝑖 , where 𝜙 𝑖 is the 𝛾 -geodesic interval between 𝐺 𝑖 (𝑉 𝑖 , 𝐸 𝑖 ) and 𝐺 𝑖−1 (𝑉 𝑖−1 , 𝐸 𝑖−1 ) , and 𝜆 is a coefficient which is larger than 1 (1.5 in our tests). Additionally, at the first layer, ℎ 𝑘 can be set as the z-coordinate 𝑧 𝑘 of vertex 𝑣 𝑘 . For the Y model, Figure 17 (a) shows the distribution of layer thickness deviation e at different layers when the geodesic distance interval 𝜙 𝑖 is 1mm ( e = ( ℎ 𝑘 - 𝜙 𝑖 ) / 𝜙 𝑖 ), with a maximum percentage deviation of 35%, and Figure 17 (b) shows the statistical results. Although the overhang angle at the boundary of the part is considerably reduced, the layer thickness is nonuniform due to the intrinsic nonuniform distribution of geodesics. 3 Figure 16
Illustration of the printing parameters
Figure 17
Illustration of layer thickness of the lattice infill patterns: (a) distribution of layer thickness deviation at different layers; (b) statistical chart of layer thickness deviation at different layers.
To avoid possible collisions when the nozzle moves from one connected sub-graph to another, in this paper, we make use of the safe box method (as reported in our recent works [40,41]). As schematically shown in Figure 18, under the safe box paradigm, the nozzle first air-moves from the current sub-graph to one of the safe planes (which are outside the current in-process workpiece), then moves on this safe plane, crosses the obstacle (i.e., the in-process workpiece), and finally approaches another sub-graph. 4
Figure 18
Safe box of the in-process workpiece and planning of the collision-free air-move of the nozzle
4. Experiments and discussion
We have implemented the proposed methodology of automatically generating the lattice infill structures and printing path for multi-axis support-free printing of an arbitrary freeform part in C++ and run the computer program on a laptop with an Intel i7 CPU. In addition, for the purpose of physical validation, as shown in Figure 19, we have built a simple multi-axis FDM printer, which is composed of a 6DOF robot arm (UR5) and a three-axis filament feed system. The robot enables the in-process workpiece to realize any desirable posture with respect to the nozzle, while at the same time the robot and the filament feed rates are synchronously controlled to ensure that Eq. (14) is always satisfied during the printing process. Five freeform parts with complicated structures such as with overhangs and of non-zero genus numbers were chosen for the test and were also physically printed. In this section we report the experimental results of both computer simulation and physical printing, and together with our discussion.
Figure 19
Homebuilt multi-axis robot printing system Figure 20 depicts the tetrahedral models of the five freeform parts that are printed by using the proposed lattice infilling method, i.e., the Y model, the spiral model, the bunny, the kitten, and the propeller. Table 1 lists the printing parameters of the five parts. The 𝛾 -geodesic distance interval for the curved layer slicing is set to be 0.6 mm. The lattice width for the first four parts are set to be 4 mm, 6 mm, 6mm, and 6 mm respectively, while the propeller is printed in two steps – in the first step, the cylinder is printed with the lattice width set to be 8 mm, and in the second step the three blades are printed with the lattice width set to be 3 mm. Figure 21 and Figure 22 show the actual and simulated printing processes of the five parts, respectively, and Figure 23 shows some cross-sections of three printed parts. All the five freeform parts are successfully printed without any supports, for both the part boundary surface and the interior infills. Due to the intrinsic nature of geodesic distance field, the overhang angle at the part boundary surface is considerably reduced, making it possible to print a part without any exterior supports. As already described, the lattice infill structures are formed by the intersections between the three clusters of orthogonal IGDSs (i.e., the 𝛾 - IGDSs, the 𝛼 -IGDSs and the 𝛽 -IGDSs.), which guarantees that the generated infill structures are self-supporting. In terms of the computational cost, as our volume decomposition method is a combination of several computational processes, Table 2 lists the time complexities of these processes and the actual amounts of running time of the first four tests. The calculation of the geodesic distance field involves solving a linear system, so the time complexity is only O ( n ), where n is the number of mesh vertices. The time complexity of the lattice infill structures’ generation is O ( n*m ), where n and m are respectively the numbers of sliced layers and mesh vertices, as the time complexity of Algorithm 1 is O ( m ), and Algorithm 1 will be executed n times to generate the lattice infill for each layer. After the trimming of the lattice infill in each layer, the well-known Fleury’s algorithm is used to find an Eulerian tour, whose time complexity is O ( n ), where n is the number of vertices in the sub-graph. The time complexity of layer thickness calculation (Eq. (15)) is O ( n*m ), where n is the number of vertices of the current sub-graph, and m is the number of edges of the previous sub-graph. By using a kd -tree, the time complexity can be reduced to O ( n* log m ). Table 1 Printing parameters of the five parts Part Number of tetrahedrons Time for printing path generation (s) Path length (mm) Actual printing time (min) Number of layers Width of the lattice (mm) Y 17148 33 49765 130 145 4 Spiral 31152 74 71139 171 401 6 Bunny 66467 203 194260 400 206 6 Kitten 58146 120 88338 232 161 6 Propeller 83314 145 125425 327 Cylinder: 183 Blade: 83 Cylinder: 8 Blade: 3 Table 2 Time complexities of the algorithms and the running time of the three tests Process Algorithm Time complexity Running time (s) Geodesic distance field generation Sec. 2.1 O ( n ) Y: 6; Spiral: 23; Bunny: 48; Kitten: 39 Lattice infill generation Algorithm 1 O ( n*m ) Y: 5; Spiral: 17; Bunny: 26; Kitten: 16 Finding a Eulerian tour Fleury’s algorithm O ( n ) Y: 14; Spiral: 23; Bunny: 77; Kitten: 39 Layer thickness calculation Eq. (15) O ( n*m ) Y: 4; Spiral: 6; Bunny: 36; Kitten: 15
Figure 20
Five freeform parts to print: (a) Y model; (b) spiral model; (c) bunny; (d) kitten; (e) propeller. (a) (b) (c) (d) (e) Figure 21
Multi-axis support-free printing with lattice infill structures: (a) Y model; (b) spiral model; (c) bunny; (d) kitten; (e) propeller. (a) (b) (c) (d) (e) Figure 22
Simulation of the printing processes: (a) Y model; (b) spiral model; (c) bunny; (d) kitten; (e) propeller. (a) (b) (c) Figure 23
Cross-sections of the printed parts: (a) Y model; (b) bunny; (c) kitten.
The lattice width (i.e., geodesic intervals for 𝛼 -GDF and 𝛽 -GDF) can be set to be a constant, which will result in even lattices (see in Figure 21 and Figure 24 (a)). However, to print a part with graded material properties (e.g., variable Young’s modulus), the lattice width can also be adjusted adpatively, namely, by increasing the density of infill lattices at certain places, as shown in Figure 24 (b). For example, the “roof” region of a printed part is typically the weakest, as the overhang angle 𝜃 there is very small. Figure 25 shows the distribution of the overhang angle 𝜃 of the bunny model, where the blue regions identify the “roofs” that are susceptible to material collapse and the density of the lattices in these places should be increased. (a) (b) Figure 24
Even and non-even lattices Figure 25
Distribution of the angle between the surface normal vector and the printing orientation
Next, we report the experimental results of three different printing sequence traversal algorithms, i.e., the benchmarking
LPT and
DPT , and our Algorithm 4 (A4), on a tree-structured part with three branches. As shown in Figure 26, the part is automatically decomposed into 151 infilling layers by our Algorithm 1 with the geodesic distance interval set at 0.6 mm, and the total number of the connected sub-graphs of the infilling layers is 311. Table 3 and Figure 27 show the simulation results of different cases. When the nozzle angle (which is denoted by NA that measures the size of the nozzle) is 75°, the printing sequence generated by the
LPT is collision-free, which requires 162 retractions and the total air-move path length is 2382 mm. However, the
DPT fails to generate a collision-free printing sequence, although the number of retractions (only 2) and the total air-move path length (only 101 mm) would be ideal. The proposed Algorithm 4 successfully plans a collision-free printing sequence with fewer retractions (24) and a shorter path length (654 mm) as compared with those of
LPT . As expected, the number of nozzle retractions and the air-move path length are inversely related with the size of the nozzle (see Figure 27). Because the calculation of PCGs for each sub-graph involves collision check, the time complexity of A4 is the largest, i.e., O ( k * m ), where k and m are the number of sub-graphs and the average number of vertices of each sub-graph, respectively. Figure 28 shows some snapshots of the actual printing processes of the A4 and LPT printing paths when the nozzle angle is 45°. (Note that we did not compare with the
DPT path since it failed to print due to the unresolvable collisions.) The comparison results clearly confirm that, when compared to the
LPT path, the filament drag problem is mitigated considerably by our A4 path owing to its significantly reduced nozzle retractions, thus leading to a much higher printing quality. 1 (a) (b)
Figure 26
A three-branch model and its IGDSs: (a) model; (b) IGDSs. Table 3 Simulation results of the three traversal algorithms Algorithm Number of retractions Air-move path length (mm) Is collision-free or not Running time (s) LPT, NA: 75° 162 2382 Yes 0 DPT, NA: 75° 2 101 No 0 A4, NA: 75° 24 654 Yes 271 A4, NA: 60° 17 371 Yes 177 A4, NA: 45° 7 227 Yes 277 A4, NA: 30° 5 146 Yes 154 A4, NA: 15° 4 125 Yes 155 A4, NA: 1° 3 108 Yes 242 (a) (b) (c) (d) (e) (f) Figure 27
Printing sequences generated by A4 under different nozzle angles: (a) 75°; (b) 60°; (c) 45°; (d) 30°; (e) 15°; (f) 1°. (a) (b) Figure 28
Actual printing processes when the nozzle angle is 45°: (a) A4, printing path length is 39946 mm, printing time is150 min; (b) LPT, printing path length is 42101 mm, printing time is 155 min.
5. Conclusion
This paper is motivated by the need of an automatic infill structure generation method for an arbitrary freeform solid part, that will ensure that both the generated infills and the given part boundary can be printed without any support under the continuous multi-axis printing configuration. Towards this objective, three mutually orthogonal geodesic distance fields embedded in the volume of the part are established. The iso-geodesic distance surfaces (IGDSs) of these three fields naturally form the curved layers of the part and the lattice infill structure, which, by aligning the nozzle orientation with the surface normal of the curve layers, guarantee that the overhang angles at both the part surface boundary and the infills are within the self-support range and thus eliminate the need of extra support. To avoid excessive nozzle retractions when printing the infills, the lattice infill pattern in each layer is first trimmed to an Eulerian graph and then a continuous printing path is constructed by using Fleury’s algorithm. In addition, we present a printing sequence optimization algorithm for establishing a total ordering of the connected lattice infills which, while respecting the collision-free requirement, tries to minimize the air-move path length of the nozzle. The results of both computer simulation and physical printing experiments have given a positive confirmation of the proposed methods. Regarding the limitations and future research, the sub-graphs of the generated infills by our method cannot always maintain the convexity if the part has a complicated topology, which may cause local collisions when the nozzle angle is large. It is conceivable that, even under the most conservative LPT, there can be cases when the collision simply cannot be avoided on a skeleton tree. One solution to this problem is using a slender nozzle to reduce the potential of local collision. On the other hand, as there are many ways to decompose a solid and generate curved slicing layers (e.g., [17] and [18]), a solid that fails our 3D geodesics-based curved-layer slicing and printing sequencing algorithm might well be printable without collisions under other strategies of slicing and printing sequencing. One plausible solution is that, rather than given a base, we freely select a base (including both its location and the size) on the boundary of the solid to define the γ-geodesic distance field and the corresponding lattice infill patterns so that their corresponding skeleton tree will be printable at least under LPT. Alternatively, for a given solid we could combine the proposed 3D geodesics-based volume decomposition with other types of decomposition and come up with a different set of curved slicing layers and their printing sequence that are better in dealing with the collision. All these will be our future research topics. 4
6. Acknowledgement
This work is supported by Hong Kong RGC-GRF/16200819.
References [1] W. Gao, Y. Zhang, D. Ramanujan, K. Ramani, Y. Chen, C.B. Williams, C.C.L. Wang, Y.C. Shin, S. Zhang, P.D. Zavattieri, The status, challenges, and future of additive manufacturing in engineering, CAD Comput. Aided Des. 69 (2015) 65–89. https://doi.org/10.1016/j.cad.2015.04.001. [2] F. Tamburrino, S. Graziosi, M. Bordegoni, The design process of additively manufactured mesoscale lattice structures: A review, J. Comput. Inf. Sci. Eng. 18 (2018) 1–16. https://doi.org/10.1115/1.4040131. [3] G. Dong, Y. Tang, Y.F. Zhao, A survey of modeling of lattice structures fabricated by additive manufacturing, J. Mech. Des. Trans. ASME. 139 (2017) 1–13. https://doi.org/10.1115/1.4037305. [4] J. Wu, C.C.L. Wang, X. Zhang, R. Westermann, Self-supporting rhombic infill structures for additive manufacturing, CAD Comput. Aided Des. 80 (2016) 32–42. https://doi.org/10.1016/j.cad.2016.07.006. [5] M. Lee, Q. Fang, Y. Cho, J. Ryu, L. Liu, D.S. Kim, Support-free hollowing for 3D printing via Voronoi diagram of ellipses, CAD Comput. Aided Des. 101 (2018) 23–36. https://doi.org/10.1016/j.cad.2018.03.007. [6] T. Kuipers, J. Wu, C.C.L. Wang, CrossFill: Foam Structures with Graded Density for Continuous Material Extrusion, CAD Comput. Aided Des. 114 (2019) 37–50. https://doi.org/10.1016/j.cad.2019.05.003. [7] W. Wang, Y. Liu, S. Member, J. Wu, S. Tian, C.C.L. Wang, S. Member, L. Liu, X. Liu, Support-Free Hollowing, 24 (2018) 2787–2798. [8] Y. Yang, S. Chai, X.M. Fu, Computing interior support-free structure via hollow-to-fill construction, Comput. Graph. 70 (2018) 148–156. https://doi.org/10.1016/j.cag.2017.07.005. [9] Y. Xie, X. Chen, Support-free interior carving for 3D printing, Vis. Informatics. 1 (2017) 9–15. https://doi.org/10.1016/j.visinf.2017.01.002. 5 [10] J. Wang, J. Dai, K.S. Li, J. Wang, M. Wei, M. Pang, Cost-effective printing of 3D objects with self-supporting property, Vis. Comput. 35 (2019) 639–651. https://doi.org/10.1007/s00371-018-1493-y. [11] J. Lee, K. Lee, Block-based inner support structure generation algorithm for 3D printing using fused deposition modeling, Int. J. Adv. Manuf. Technol. (2017). https://doi.org/10.1007/s00170-016-9239-3. [12] T. Lee, J. Lee, K. Lee, Extended block based infill generation, Int. J. Adv. Manuf. Technol. 93 (2017) 1415–1430. https://doi.org/10.1007/s00170-017-0572-y. [13] P. Gupta, B. Krishnamoorthy, G. Dreifus, Continuous toolpath planning in a graphical framework for sparse infill additive manufacturing, Comput. Des. (2020). https://doi.org/10.1016/j.cad.2020.102880. [14] J. Wu, C. Dick, R. Westermann, A System for High-Resolution Topology Optimization, IEEE Trans. Vis. Comput. Graph. 22 (2016) 1195–1208. https://doi.org/10.1109/TVCG.2015.2502588. [15] J. Wu, N. Aage, R. Westermann, O. Sigmund, Infill Optimization for Additive Manufacturing-Approaching Bone-Like Porous Structures, IEEE Trans. Vis. Comput. Graph. 24 (2018) 1127–1140. https://doi.org/10.1109/TVCG.2017.2655523. [16] M. Langelaar, Topology optimization of 3D self-supporting structures for additive manufacturing, Addit. Manuf. (2016). https://doi.org/10.1016/j.addma.2016.06.010. [17] K. Hu, S. Jin, C.C.L. Wang, Support slimming for single material based additive manufacturing, CAD Comput. Aided Des. (2015). https://doi.org/10.1016/j.cad.2015.03.001. [18] Y. Zhou, H. Lu, Q. Ren, Y. Li, Generation of a tree-like support structure for fused deposition modelling based on the L-system and an octree, Graph. Models. 101 (2019) 8–16. https://doi.org/10.1016/j.gmod.2018.12.003. [19] J. Vanek, J.A.G. Galicia, B. Benes, Clever support: Efficient support structure generation for digital fabrication, Eurographics Symp. Geom. Process. (2014). https://doi.org/10.1111/cgf.12437. [20] T. Tricard, F. Claux, S. Lefebvre, Ribbed Support Vaults for 3D Printing of Hollowed Objects, Comput. Graph. Forum. 39 (2019) 147–159. https://doi.org/10.1111/cgf.13750. 6 [21] Y. Gao, L. Wu, D.M. Yan, L. Nan, Near support-free multi-directional 3D printing via global-optimal decomposition, Graph. Models. 104 (2019) 101034. https://doi.org/10.1016/j.gmod.2019.101034. [22] X. Wei, S. Qiu, L. Zhu, R. Feng, Y. Tian, J. Xi, Y. Zheng, Toward Support-Free 3D Printing: A Skeletal Approach for Partitioning Models, IEEE Trans. Vis. Comput. Graph. 24 (2018) 2799–2812. https://doi.org/10.1109/TVCG.2017.2767047. [23] C. Wu, C. Dai, G. Fang, Y.J. Liu, C.C.L. Wang, General Support-Effective Decomposition for Multi-Directional 3-D Printing, IEEE Trans. Autom. Sci. Eng. (2020). https://doi.org/10.1109/TASE.2019.2938219. [24] P.M. Bhatt, R.K. Malhan, P. Rajendran, S.K. Gupta, Building free-form thin shell parts using supportless extrusion-based additive manufacturing, Addit. Manuf. 32 (2020) 101003. https://doi.org/10.1016/j.addma.2019.101003. [25] E. Karasik, R. Fattal, M. Werman, Object Partitioning for Support-Free 3D-Printing, Comput. Graph. Forum. 38 (2019) 305–316. https://doi.org/10.1111/cgf.13639. [26] H. ming Zhao, Y. He, J. zhong Fu, J. jiang Qiu, Inclined layer printing for fused deposition modeling without assisted supporting structure, Robot. Comput. Integr. Manuf. 51 (2018) 1–13. https://doi.org/10.1016/j.rcim.2017.11.011. [27] H. Liu, L. Liu, D. Li, R. Huang, N. Dai, An approach to partition workpiece CAD model towards 5-axis support-free 3D printing, Int. J. Adv. Manuf. Technol. 106 (2020) 683–699. https://doi.org/10.1007/s00170-019-04495-3. [28] H. Zhang, Z. Sun, Q. Hu, H. Lammer, Research and implementation of integrated methods of unsupported printing and five-axis CNC machining technology, Teh. Vjesn. 26 (2019) 1267–1274. https://doi.org/10.17559/TV-20180313125041. [29] M. Wang, H. Zhang, Q. Hu, D. Liu, H. Lammer, Research and implementation of a non-supporting 3D printing method based on 5-axis dynamic slice algorithm, Robot. Comput. Integr. Manuf. (2019). https://doi.org/10.1016/j.rcim.2019.01.007. [30] C. Wu, C. Dai, G. Fang, Y.J. Liu, C.C.L. Wang, RoboFDM: A robotic system for support-free fabrication using FDM, in: Proc. - IEEE Int. Conf. Robot. Autom., 2017. https://doi.org/10.1109/ICRA.2017.7989140. 7 [31] K. Xu, L. Chen, K. Tang, Support-Free Layered Process Planning Toward 3 + 2-Axis Additive Manufacturing, IEEE Trans. Autom. Sci. Eng. (2019). https://doi.org/10.1109/TASE.2018.2867230. [32] C. Dai, C.C.L. Wang, C. Wu, S. Lefebvre, G. Fang, Y.J. Liu, Support-free volume printing by multi-axis motion, ACM Trans. Graph. 37 (2018). https://doi.org/10.1145/3197517.3201342. [33] K. Xu, Y. Li, L. Chen, K. Tang, Curved layer based process planning for multi-axis volume printing of freeform parts, CAD Comput. Aided Des. (2019). https://doi.org/10.1016/j.cad.2019.05.007. [34] K. Crane, C. Weischedel, M. Wardetzky, Geodesics in heat: A new approach to computing distance based on heat flow, ACM Trans. Graph. (2013). https://doi.org/10.1145/2516971.2516977. [35] T. Nguyen, K. Karčiauskas, J. Peters, C1 finite elements on non-tensor-product 2d and 3d manifolds, Appl. Math. Comput. (2016). https://doi.org/10.1016/j.amc.2015.06.103. [36] A.G. Belyaev, P.A. Fayolle, On Variational and PDE-Based Distance Function Approximations, Comput. Graph. Forum. (2015). https://doi.org/10.1111/cgf.12611. [37] J. Etienne, N. Ray, D. Panozzo, S. Hornus, C.C.L. Wang, J. Martínez, S. McMains, M. Alexa, B. Wyvill, S. Lefebvre, Curvislicer: Slightly curved slicing for 3-axis printers, ACM Trans. Graph. (2019). https://doi.org/10.1145/3306346.3323022. [38] Y. Tong, S. Lombeyda, A.N. Hirani, M. Desbrun, Discrete multiscale vector field decomposition, in: ACM SIGGRAPH 2003 Pap. SIGGRAPH ’03, 2003. https://doi.org/10.1145/1201775.882290. [39] P.L. George, E. Seveno, The advancing‐front mesh generation method revisited, Int. J. Numer. Methods Eng. (1994). https://doi.org/10.1002/nme.1620372103. [40] Y. Li, L. Zeng, K. Tang, C. Xie, Orientation-point relation based inspection path planning method for 5-axis OMI system, Robot. Comput. Integr. Manuf. (2020). https://doi.org/10.1016/j.rcim.2019.101827. [41] Y. Li, K. Tang, L. Zeng, A Voxel Model-Based Process-Planning Method for Five-Axis Machining of Complicated Parts, J. Comput. Inf. Sci. Eng. 20 (2020) 1–14. https://doi.org/10.1115/1.4046589. 8[31] K. Xu, L. Chen, K. Tang, Support-Free Layered Process Planning Toward 3 + 2-Axis Additive Manufacturing, IEEE Trans. Autom. Sci. Eng. (2019). https://doi.org/10.1109/TASE.2018.2867230. [32] C. Dai, C.C.L. Wang, C. Wu, S. Lefebvre, G. Fang, Y.J. Liu, Support-free volume printing by multi-axis motion, ACM Trans. Graph. 37 (2018). https://doi.org/10.1145/3197517.3201342. [33] K. Xu, Y. Li, L. Chen, K. Tang, Curved layer based process planning for multi-axis volume printing of freeform parts, CAD Comput. Aided Des. (2019). https://doi.org/10.1016/j.cad.2019.05.007. [34] K. Crane, C. Weischedel, M. Wardetzky, Geodesics in heat: A new approach to computing distance based on heat flow, ACM Trans. Graph. (2013). https://doi.org/10.1145/2516971.2516977. [35] T. Nguyen, K. Karčiauskas, J. Peters, C1 finite elements on non-tensor-product 2d and 3d manifolds, Appl. Math. Comput. (2016). https://doi.org/10.1016/j.amc.2015.06.103. [36] A.G. Belyaev, P.A. Fayolle, On Variational and PDE-Based Distance Function Approximations, Comput. Graph. Forum. (2015). https://doi.org/10.1111/cgf.12611. [37] J. Etienne, N. Ray, D. Panozzo, S. Hornus, C.C.L. Wang, J. Martínez, S. McMains, M. Alexa, B. Wyvill, S. Lefebvre, Curvislicer: Slightly curved slicing for 3-axis printers, ACM Trans. Graph. (2019). https://doi.org/10.1145/3306346.3323022. [38] Y. Tong, S. Lombeyda, A.N. Hirani, M. Desbrun, Discrete multiscale vector field decomposition, in: ACM SIGGRAPH 2003 Pap. SIGGRAPH ’03, 2003. https://doi.org/10.1145/1201775.882290. [39] P.L. George, E. Seveno, The advancing‐front mesh generation method revisited, Int. J. Numer. Methods Eng. (1994). https://doi.org/10.1002/nme.1620372103. [40] Y. Li, L. Zeng, K. Tang, C. Xie, Orientation-point relation based inspection path planning method for 5-axis OMI system, Robot. Comput. Integr. Manuf. (2020). https://doi.org/10.1016/j.rcim.2019.101827. [41] Y. Li, K. Tang, L. Zeng, A Voxel Model-Based Process-Planning Method for Five-Axis Machining of Complicated Parts, J. Comput. Inf. Sci. Eng. 20 (2020) 1–14. https://doi.org/10.1115/1.4046589. 8