Robotic FDM for free-form fabrication: evaluating adaptive non-planar slicing with different contour methods

Additive manufacturing (AM) is a family of processes that builds 3D objects with materials like polymers, ceramics and alloys (Anerao et al., 2024; Bhuvanesh Kumar and Sathiya, 2021; Chen et al., 2022; De Assis and Rampazo, 2024; Pascual-González et al., 2022; Shahrubudin et al., 2019; Solomon et al., 2021) and for different application fields from automotive, aerospace and civil applications to jewelry, musical and pharmaceutical (Bacciaglia et al., 2019; Breseghello and Naboni, 2022; Fabrizio et al., 2022; Fatma et al., 2021; Joshi and Sheikh, 2015; Melocchi et al., 2016; Mohammadi Zerankeshi et al., 2022). It adds material layer by layer, allowing the object fabrication with a single process. Among AM processes, fused deposition modeling (FDM) is the most common and widely adopted for fast prototyping (Salifu et al., 2022). The slicing process outputs the deposition toolpath, taking as input the mesh representation of the object, commonly as triangular mesh and some process parameters. Conventional slicing works by cutting the mesh object with parallel planes to find the contour toolpath, which is also the starting point for filled components. Conventional approaches are very suitable to realize objects without curved features. Nevertheless, this approach, also known as 2D-slicing, has significant limitations for free-form shapes that present curvatures and undercuts. In real components, it causes the staircase effect and the need for support structures, affecting the surface quality and process material efficiency. In fact, on overhanging sides with strong curvatures, the layer distance on the contour, computed as Hausdorff distance, increases, leading to the need for support structures (Xu et al., 2019). This layer thickness variation (LTV) is present on the boundaries of overhanging sides with 2D-slicing methods and it is null when no curvatures or overhangs are present. Previous works (Insero et al., 2022, 2023) have discussed this issue of 2D-slicing and planar multi-directional slicing approaches in free-forms. The defects of 2D-slicing methods can be avoided in free-forms using new slicing approaches, such as non-planar (NP) slicing techniques. NP slicing allows supportless fabrication depositing material outside planes and in the entire 3D space (Nayyeri et al., 2022). Such freedom in deposition is enabled only by increasing the number of degrees of freedom (DoFs) between the tool head and the component to fabricate. The introduction of anthropomorphic robots in AM can enable NP slicing providing six or more DoFs to the motion (Pires et al., 2022). Besides, anthropomorphic robots are a compact solution providing large-scale fabrication with a relatively small machine (Bhatt et al., 2020; Rebaioli et al., 2019). NP slicing techniques can rely on optimization approaches acting on geometrical parameters, e.g. cusp height and layer thickness (LT). Several works have been done in the field. Some of the most relevant for the authors are (Chalvin et al., 2019; Feng et al., 2021; Huang and Singamneni, 2015; Jin et al., 2017; Mitropoulou et al., 2020, 2022; Singamneni et al., 2012; Xu et al., 2019; Zhao et al., 2018). Nevertheless, the majority of works are restricted to thin-walled and/or parametric geometries that are not fully representative of requirements in AM parts. In (Chalvin et al., 2019), the authors formalize the problem of LTV for the thin-walled toroidal geometry applying geometrical and optimization methods to minimize the LTV using a parametric representation of the geometry. Related works of (Mitropoulou et al., 2020, 2022; Xu et al., 2019) have proposed a scalar-field approach to generate NP contours on objects with a mesh representation. Contours are generated as isolines of the geodesic scalar field computed on the object surface with the algorithm introduced in (Mitchell et al., 1987). On the other hand, Chalvin et al. (2019) and Xu et al. (2019) show that contours generated by both optimization and iso-geodesic methods provide 3D curves with two issues. From one side, the contour presents sharp edges that can create issues with the robotic motions due to fast tool orientation changes, a reduced tool central point (TCP) speed and possible tool collisions. On the other side, over-extrusion problems are still present due to the LTV in curved layers. To solve this problem an adaptative approach applied to NP slicing has been suggested by (Xu et al., 2019). Table 1 reports a collection of recent and relevant algorithms regarding NP slicing. The works are categorized according to the type of slicing approach, i.e. iso-parametric, optimization according to a specific functional and geodesic approach. Within that categorization, methods related to thin-walled fabrication (Contour) are distinguished from filled fabrication. Despite several investigations being done on NP slicing algorithms, to the authors’ knowledge, the LTV behavior on the inner side has not received enough attention both on the real fabrication and on the theoretical toolpath. In a previous work (Insero et al., 2023), it has been shown the LTV in the inner side for a quarter of the torus, using an iso-v parametrization for the contour generation and from which the inner side has been reconstructed using a projection method. In the present work, the authors contribute to:

• Evaluate, for the first time, different NP slicing techniques on the same geometry, in terms of LTV. The comparison is done on the quarter of torus considering the limits of a parametric representation and moving on a mesh representation that is more representative of a real object.

• Present a method for the reconstruction of the inner side of a layer with the complexity of curvature using B-spline patches.

• Present a new adaptive slicing method based on sublayers, to reduce LTV issues, for the different NP slicing methods.

• Compare, for the first time, different NP slicing techniques from theoretical and experimental sides.

In detail, in Section 2, the methodological approach and the geometry modeling will be presented. The contouring methods applied for the parametric and mesh-based representation are presented in Section 3. Then, Section 4 presents inner surface reconstruction methods and the infill strategy. Section 5 describes the adaptive approach based on sublayer. In Section 6, the setup involved is detailed. Section 7 discusses the results. To evaluate the possibility of the presented method being generalized to other complex shapes, a bifurcated and a bifurcated chain case studies are shown in Section 8. Finally, the conclusions of the work are drawn in Section 9.

The flow of the research is depicted in Figure 1. The primary goal of the NP slicing algorithms is to optimize the LTV, i.e. reach a constant LT across the component. As discussed in the introduction, consistent LTV reduces anisotropy, minimizes staircase effects, improves mechanical performance and promotes supportless fabrication. This performance evaluation, in terms of LTV, uses a representative object: a quarter torus. The quarter of torus offers both overhang and constant curvature making it ideal as a representative object of free-form, including the drawbacks of curvature and overhangs. Besides, it has been considered because it can be described with a parametric (Chalvin et al., 2019; Insero et al., 2022, 2023) and a mesh representation. Quarter of the torus is represented as a revolution geometry that uses two angles (u, v), the constant curvature radius R and the torus radius r. Three contouring algorithms are applied to the selected geometry. The NP slicing is composed of contouring stage and infill stage. For the contour stage, the choice between the possible approaches depends on the input type. For the parametric representation, two approaches have been used: an approximated solution, referred to us as iso-v, and a numerical optimization method. All of them perform contouring on the external surface of the input model. The third approach can be applied to discrete 3D objects represented as triangular meshes, enabling a generalization of the NP slicing to general CAD models. In a real-case scenario, CAD objects are used in an approximated form with surface mesh made of triangular elements (i.e. STL and OBJ formats) (King et al., 2021; Stroud and Xirouchakis, 2000). Iso-geodesic method applies to the mesh discretized representation. Following the contour generation stage, the infill stage is carried out. First, the contour’s internal area is reconstructed with B-spline patches. B-spline patches are chosen as interpolation methods due to their flexibility and capability of providing smooth surfaces. Then, a bidirectional linear pattern has been chosen for the toolpath. Bidirectional linear patterns are commonly adopted to produce filled components in different AM processes, moving from FDM conventional 2D-slicing to other AM processes, such as laser metal deposition (Chen et al., 2022; Dinovitzer et al., 2019; Donadello et al., 2022; Pant et al., 2023; Ribeiro et al., 2020). After, the LT is evaluated by computing the minimum distance between each point of the next layer and the previous one in both the contour and inner side. The thickness is theoretically estimated because all three methods operate on the contour distance but not in the infill. If the theoretical LTV overcomes a threshold imposed by component geometry, then an adaptive slicing strategy is activated. The adaptive slicing adjusts the nominal LT and iterates the slicing process with the adjusted thickness. Finally, to evaluate the real LTV generated in the process, the components are cross-sectioned and the LT are measured with an optical microscope, and then compared with the LTV of the theoretical toolpath.

This section presents the three contouring algorithms applied in this work. Section 3.1 presents the two algorithms for parametric representation. First, the approximated method, iso-v, and then, the numerical optimization of the LTV is shown. Both methods provide a geometry parametrization according to the fixed constraint. This mathematical approach requires geometry parametric representation that is not always available. Section 3.2 presents the mesh algorithm. It is based on the iso-geodesic curves and it can be applied to triangle-meshed objects (i.e. STL, OBJ formats).

The torus can be represented by equation (1) (Chalvin, 2021; Chalvin et al., 2019; Insero et al., 2023). Where R and r represent the curvature and torus radii respectively. Whereas v and u are the angular parameter of the circular section and the revolution angle, respectively. The v parameter has been discretized equally spaced in the interval [0; 2π]. Considering a number of contour points equal to 100 and the radius r = 13 mm, the arc length for the first layer is $0.82 mm$ :

The first contouring method exploits the parametric representation providing a rough approximation of LTV optimization. Considering a nominal LT, Δh, the layer distance between two contiguous layers is approximatively considered as the Euclidean distance between two points at the same angular position v. This distance, with iso-v parameter, can provide a geometry parametrization where the u parameter is a function of the nominal LT Δh, the layer index n and v. Considering that u₀ = 0 (i.e. starting with angle zero for the first layer) and applying the Euclidean distance to equation (1) between two points with same v. parameter and consecutive u_n, we can achieve equation (2), which expresses this parametrization. Considering the curvature of the geometry, the real minimum LT cannot be the distance between two points with the same v, as reported in Figure 2. Figure 2 shows iso-v points in two contiguous layers. In that representation, it is evident that the distance between the point P_n, vi (i.e. the red point in Figure 2) and the iso-v point of the previous layer P_n−1, vi (i.e. the blue point in Figure 2) is not the minimum. In fact, the real closest point is the green one in Figure 2 (Chalvin et al., 2019). This means that the iso-v contouring method introduces an appreciable error in the theoretical calculation of LTV:

To achieve a global constancy of the desired nominal LT (Δh), an optimization approach has been introduced starting from the second layer. The objective function is the error between the actual layer inter-contour distance, δh and the desired one Δh. δh is defined as the minimum distance between the considered point (⁠$δu,vi$⁠) and the lower layer (Chalvin, 2021; Chalvin et al., 2019). The problem statement is expressed in equation (3) where the objective function to minimize is ϕ, which is the squared error between the desired Δh and the actual layer inter-contour distance. The function ϕ is subjected to an inequality constraint ψ_i that controls the contouring algorithm to construct the geometry from the base (u₀ = 0) to an increasing u (i.e. the layering process will proceed to construct the entire geometry). The actual LT $δhun,vi$ of a point i of the layer n is computed as the minima height of the triangle, which has as vertexes the considered points and the points of each segment of the previous layer as shown in Figure 2. The problem results in a non-linear programming problem with a quadratic objective function, which can be solved using the sequential least square quadratic programming algorithm (Nocedal and Wright, 2006) and the minimize function in the Scipy library for Python 3 (Virtanen et al., 2020). The parameters chosen for the optimization are the maximum number of iterations and the tolerance of the objective function. The maximum iteration number is set to 1000, but it is never reached. The tolerance on the objective function is set to 10^{−8. $pn-1$} is the number of points of the layer $n-1$⁠:

Most applications, however, do not use parametric object representation as slicing input. In real-case scenarios, input models consist of a triangulated surface representation (Ding et al., 2016; King et al., 2021). To broaden the problem to generic STL objects, a NP slicing contouring is presented based on computing a parameter called geodesic distance, defined as the shortest path length embedded on a surface between two points. Computing the minimum geodesic distance between all the triangle vertexes of the object with respect to a fraction of vertexes, identified as base vertexes, allows to definition of a scalar field on the geometry. The computation of the minimum geodesic distance in any mesh vertex can be achieved using an implementation based on the Mitchel–Mount–Papadimitriou (MMP) algorithm (Mitchell et al., 1987). The MMP algorithm (also known as exact geodesic algorithm) has been implemented by Surazhsky et al. (2005), providing an implementation that computes faster geodesic field on triangle meshes. To generate the NP-contours, isolines are extracted from the scalar field, that is interpolated on the edges between vertices. The NP-contours are taken as isolines equally spaced by the nominal LT Δh value (Mitropoulou et al., 2020; Xu et al., 2019). Isoline’s points (i.e. P_n) are all points $P∈∂Ω$ that have the geodesic distance multiple of the nominal LT value Δh as shown in equation (4). Figure 3 shows a schematization of geodesic contouring. The left figure schematizes the difference between the geodesic distance $dg*$ and the Euclidean distance d_e computed between two points P_b,i and P_d,i on the generic surface $∂Ω$⁠. The geodesic scalar field uses the vertexes on the boundary curve I_b as starting points. Then, the geodesic distance value assigned to each vertex is the minima distance $dg*(Ib, Pd,i)$⁠. A depiction of the generated scalar field for the quarter of torus is shown in Figure 3 on the right. The scalar field is computed between all the vertices starting from the bottom. The colors represent the growing distance from the base. The black lines are isolines equally spaced of the scalar field. As highlighted by Mitropoulou et al. (2020), geodesic field suffers from the triangle meshes’ irregularities. For the purposes of this work, discretization issues are considered out of scope. To avoid discretization issues, higher discretization quality settings have been used to generate the STL geometry:

This section will introduce the inner surface reconstruction method using the contour points provided by the presented NP slicing methods. First, the B-spline surface interpolation is used to recover an interpolation per layer. Then, the bidirectional linear toolpath is generated projecting the layer on a plane and using a grid for the line discretization. The layer’s interpolated points are the ones inside the contour. The contour must be a closed curve with a single area to be interpolated. The required inputs for the infill generation are the contour points and the degree parameters along the two main directions.

Compared to the conventional planar case, where the inner side lays on a plane. In the curved layers, the inner side can be any surface bounded by the contour. Hence, a suitable surface on which to embed the toolpath must be uniquely determined. A novel approach is presented to reconstruct internal layer information based on a B-spline surface interpolation starting from contour points. B-spline surfaces are composed of two basis functions (one per direction), each of them is characterized by a degree (m and n respectively), which needs to be tuned according to layers curvature to achieve the best fitting result. A general formulation of a B-spline patch with degrees (m, n) is defined in equation (6). The points used for the interpolation are called control points (P_i,j). The two basis functions (B_i,m and B_j,n) can be recursively computed using the Cox-deBoor algorithm (Cox, 1972; De Boor, 1972) reported in equation (5). The basis function B_i,p (t) of degree p is recursively computed with basis function with order p − 1 in each node segment [t_i, t_i+1]. B-splines provide smooth surfaces and choose the polynomial degree along the two directions independently computing different basis functions. The B-spline surface s(u, v) in two generic directions (u, v) is computed as the combination of the basis function along the two directions. The contour points have been interpolated using the B-spline implementation provided by SciPy (bisplrep and bisplev) based on (Dierckx, 1981, 1996). Then, a discretized grid is defined around the contour points. To extrapolate the points inside the contour region, first the contour is projected on the horizontal plane. Then, the grid points exceeding the projected curve are discarded (i.e. assigning NaN to the grid point). To find the grid points inside the projection, a polygon is constructed from the projection, and then, each grid point is checked whether it is inside the polygon using the Shapely library (Gillies et al., 2022). Step 1 depicted in Figure 4 represents the output of the B-spline interpolation. It shows a planar grid embedding the layer projected. From this, the linear bidirectional infill is generated:

An alternating bidirectional infill strategy is adopted for filling each layer after surface interpolation. Infill points are extracted from the surface to form parallel lines placed at a fixed distance, named hatching distance (H_d), selected on the base of the AM technique and the related process parameters. Each line is then linearly interpolated and resampled using a discretization factor η according to the line length. The η parameter adjusts the tradeoff between path accuracy and the blending effect in the robot speed which is determined experimentally. A picture of the algorithm to gather the points to construct the infill lines is depicted in Figure 4. To achieve an infill direction along axis k, in Step 1, the algorithm scrolls the j-axis until it finds the first not-NaN value (i.e. a value that is inside the projection of the contour on the plane) in the k-axis (i.e. P_k5,j0 in the picture), this point is considered as the starting point of the line, Figure 4a. In Step 2, the algorithm is scrolling along the j-axis choosing the j-index that has the point with the distance closer to the hatching distance concerning the previous point (within a tolerance ε), Figure 4b. The j-axis is scrolled until it reaches the maximum grid index, as shown in Step 3 , Figure 4c. Knowing all the j-indexes spaced of H_d, the boundary line points are found at the margin of the k-axis as shown in Steps 4 and 5, Figure 4f-e. The complete line considers the k-indexes between the two-line margins reduced by the factor η as shown in Step 6, Figure 4d. Collected points are then organized to form a linear bidirectional printing path, alternating the printing direction of 90° each layer, which will be followed by the robot during the process. To always guarantee a good deposition during the printing, the normal to the surface in the toolpath has been used. Its orientations are locally calculated from local building directions coincident with the surface normal of each infill point computing using the Frenet reference frame.

The surface normal for each infill point is computed as the binormal versor of the Frenet reference frame. The tangent vector T_j,k is computed as the normalized difference between two consecutive points on the infill line direction (i.e. if the line direction is along the j-axis, the tangent vector T_j,k is computed by normalizing the difference between P_j+1,k and P_j,k). While the normal direction is computed similarly with the previous parallel line. The versors computation is shown in equation (7):

The Euler angles are determined considering the orientation of the binormal with respect to the vertical axis. The rotation along the z-axis is not considered. The two angles are then computed as the following projecting angle around the z-axis:

The surface reconstruction and infill generation return a filled representation of the layers. Compared to the contour, the surface generation and infill path do not consider the layer distance between consecutive layers. As so far considered for the contour distance, the inter-layer distance (τ_n,i) must be computed for any inner side point i of layer n. Inter-layer distance is defined as the minimum Euclidean distance between the selected point i on layer n and all the points belonging to the layer n − 1 shifted in the grid along the j, k-axes, as shown in equation (9), where n_G is the grid discretization factor. τ_n,i should coincide with the nominal LT, Δh. The theoretical LT on the contour and in the inter-layer has been evaluated with equation (9) for all three contouring methods. In the results section, the theoretical results are reported for all the cases:

A previous work (Insero et al., 2023) shows that the LT on the inner side reduces significantly if it is compared to the contour. In case the LT reduces significantly, over-extrusion defects appear. In particular, it can be noticed that as the layer deforms, the layer’s projection becomes closer to the Z-coordinate due to the geodesic distance definition. Figure 5 shows a schematization of this phenomenon. A tentative approach to solve this problem is to increase the LT (δ) used for slicing if the inter-layer distance error falls under a specified limit (τ_min) after the surface reconstruction phase. This new LT δ will be greater than the nominal Δh used for slicing the contour. We refer to this method as adaptive slicing. The adaptive condition is enabled whenever the internal LT drops below the τ_min limit. When enabled, it increases the distance between contours, iterating the contour slicing phase, until the mean value of τ in the entire layer becomes equal to the desired LT value (δ_adp). Besides, it is guaranteed that τ is punctually greater than τ_min. Although increasing the layer height δ_adp ensures having τ over the minimum limit, it implies also increasing its maximum value. If a small inter-layer distance generates layers overlapping, a higher value means a lack of adherence between two layers. To limit the issue two logical solutions are possible. The first one consists of varying the material feed rate. Theoretically, the first solution allows us to achieve the best result, however, it is not reliable. Varying the material feed rate means varying either the robot speed or the extrusion speed. Those cannot be modified instantaneously during the process, due to the entire system dynamics. Hence, the deposition process does not fit well for fast interruptions or fast material variations. The second solution, instead, consists of filling the regions in which the τ is too large using sublayers between two layers. These are smaller in size than the original ones. They are placed just in portions where the thickness exceeds the τ_max value, resulting in an additional material where it is required. The main drawback of using sublayers concerns the thickness variability. This method does not guarantee a uniform thickness within a layer. It allows to constrain its value between a minimum and a maximum. The minimum is guaranteed by the adaptive condition while the maximum is by the sublayer’s presence. Figure 6 shows a schematization between the nonadaptive, adaptive on the contour, and the last one, the addition of the sublayer. In the nonadaptive case, the minimum distance between the layers n + 1 and n is depicted as min τ_{i,n +1}. Then, the adaptive condition allows to define a newer δ_adapt for the slicer to get a minimum distance on the layer, which is greater than τ_min. Finally, for the adaptive with sublayers case, the points that have a τ_{n + 1,i} greater than τ_max on the layer $ln+1adapt$ are detected as critical points to define sublayers. Critical points on the surface layer are projected along their normal direction for a distance equal to half of their distance from the previous layer (orange arrows in Figure 6). This enables the placement of a new point between the two layers, which is again interpolated through a B-spline surface, generating the sublayers (as shown in the purple line in Figure 6). This time, the introduced points are not representative of the layer contour but contain internal points providing additional constraints for interpolation. The adaptive slicing with sublayer has been applied only to the most critical slicing methods namely, iso-geodesic and error optimization approaches, whereas the iso-v method has been used as a reference.

The robot-assisted additive manufacturing (RAAM) is a robotic cell consisting of a six-axis industrial anthropomorphic robot (IRB 2400/16, ABB) with an IRC5 controller. RAAM system is laid out in the fixed extruder-movable bed configuration. The layout choice is motivated by generated tool paths. They involve high orientation angles that reduce the pellet flow rate in the extruder. The robot is equipped with a squared 3D-printed plate with sides of 60 mm and 10 mm thickness. The bed plate is covered with blue tape (Scotch Blue Tape 2090 3M) and the spray glue for AM (3DLAC) is added to enhance print adhesion. An FDM pellet extruder (Mahor xyz v3) is fixed on an aluminum frame at a height of about 1,300 mm with respect to the robot base reference frame. The extruder is controlled by a microcontroller ATmega2560 with Ethernet onboard (EtherMega, Freetronic) and a PCB shield (RAMPS v1.6, BIGTREETECH) that allows a voltage supply of 12V for the heat cartridge and stepper motor drive. Figure 7 shows a schematization of the RAAM setup. The robot controller uses 3 RAPID tasks:

1. T_ROB manages the robot motion and receives the print information from a PC;

2. T_EXTR controls the FDM extruder; and

3. T_PAR collects robot information during the prints.

It saves robot poses from the robot, TCP speeds (V_TCP) and time with a sampling rate of 10 Hz. The T_ROB task receives messages over transmission control protocol/internet protocol (TCP/IP). Each message contains the pose information, i.e. x-y-z coordinates of the movable bed, the orientation with the Euler XYZ convention, TCP speed and the extruder screw speed. Blended motion is achieved using a pose buffer in the T_ROB task continuously received from the PC runtime the robot motion. The primitive motion used is linear, i.e. MoveL, followed by the continuous path parameter (zonedata) z0. The computed trajectory is reduced to a pose matrix including V_TCP and extruder command. Each pose is sent one by one to the robot controller filling a pose buffer in the T_ROB task. The FDM material used in this work is PLA with a globular shape with a diameter ranging from 3.5 mm to 5.0 mm (sold by FiloAlfa). Process parameters involved in the experimental campaign are collected and summarized in Table 2.

The printing strategy employed in this work consists of printing only the infill path of the geometry using the bidirectional line strategy to highlight the LT. Three components are fabricated, one per NP slicing method discussed, i.e. iso-v, optimization method and iso-geodesic. The approximated iso-v method is filled using the B-spline surface reconstruction strategy and it is used as a reference case. The optimization and iso-geodesic methods require the use adaptive approach in the filling. All the parameters used to generate the torus, the slicing parameters and the infill parameters are described in Table 3. The printed components are then manually cut and grinded. The cross-sections are analyzed under the optical microscope to find a match with theoretical results in terms of inter-layer distance by measuring distances between centers of two consecutive layers. Cross-sections are divided into different regions and analyzed in couples (on the inner curvature and the outer curvature). Eight measures per region are taken to compute mean µ and standard deviation σ.

The results section is divided into the following. Section 7.1 presents the theoretical results of contour case comparison. Section 7.2 shows the theoretical results of the intra-layer distances in filled components. Section 7.3 presents the results after the implementation of the adaptive approach. Section 7.3.1 present the theoretical results, whereas Section 7.3.2 presents the experimental results and compares with the theoretical ones, validating the modeling section.

The three contouring methods provide three different toolpaths to compare. First, the analytical-based method will be compared, i.e. the iso-v and optimization methods. Then, the iso-geodesic and optimization approaches. Figure 8 shows a side view of the toolpath between iso-v and optimization methods. It is possible to notice that the optimization and iso-geodesic approaches provide an increased layer curvature variation compared to iso-v. In both analytical cases, the upper layers present an open curve due to the fact those points are not present in the analytical representation. To apply the infill method, a straight line will be used to close the contour before applying the surface reconstruction method (Insero et al., 2023). Instead, for the iso-geodesic method, the contours are automatically closed. This is due to the meshed geometry, which represents the object as a closed surface. Leaving aside the closure of the top layers, the iso-geodesic method and optimized methods provide almost the same toolpath. This is demonstrated in Figure 9, where the Hausdorff distance is used as a comparison method between the layers (Alt et al., 2004). A greater Hausdorff distance means a greater difference between the two objects. A greater Hausdorff distance is present in the closure part. This difference is highlighted in Figure 10, where the two toolpaths are shown. The contour layer distance has been measured considering the minimum Euclidian distance between all the points of the previous layer. The distribution per layer of the distances of all the methods is shown in Figure 11. It shows the τ error percentage between the nominal layer height and the effective layer distance. The τ in the iso-v case shows a reduction of layer distance of 46%, whereas a 0% variation in the optimization case, as demonstrated in Chalvin et al. (2019). For the iso-geodesic case, the τ is almost 0% excluding discontinuity points in the closures and in the sharper edges on the back of the layer. This variation is mainly due to the quality of the discretization of the mesh object. The different contour closures and the discretization effect in the iso-geodesic will affect the LTs due to a different interpolated surface.

The B-spline interpolation method has been applied to the three contours. From the B-spline surface reconstruction, the toolpath generated is extracted. Figure 12 shows the result of the B-spline surface reconstruction algorithm and the toolpath generated for two consecutive layers. In light blue, it is possible to notice the interpolating surface, in dark blue the discrete points are used to create the path, and then they are ordered to represent the toolpath shown as a red line. The LT has been computed on the sampled point of the toolpath for all three cases. Figure 13 shows the trends of the LTV τ with respect to the nominal Δh = 0.9 mm. The iso-v case shows a monotonic trend where the error value grows almost linearly up to −46% following the contour case. For the optimization case, the LTV reduces up to −94%. In particular, the LT is reduced as the number of layers increases located mainly on the external overhang (external side of the curvature). Instead, iso-geodesic methods show a more “robust” trend for the minimum values lying almost at 0%, but with a less uniform variation compared to the optimization case. This is due to the discontinuities on the original contours which affects the interpolation. For both cases, the maximum values are double concerning the iso-v method. The optimization and iso-geodesic reach maximum values of −94% and −98%, respectively. These results mean the two methods provide a toolpath where the nozzle tip can collide with the previous layer. This phenomenon is an indirect effect of setting constant distance between points on layer contours. That induces a process through which the growing layer direction passes from completely vertical to almost horizontal along the curvature. This effect is particularly emphasized in correspondence of the external curvature region, as shown in Figure 5. These intersections aim to figure out an adaptive approach that changes the contours according to the actual τ computed. As shown, τ becomes critical in optimization and iso-geodesic methods. To limit τ variations and particularly the maximum, the adaptive approach has been applied to these two cases.

The torus geometry has been fabricated using all the three presented contouring methods. All the methods use the B-spline approach for reconstructing the inner side. Optimization and iso-geodesic use the adaptive method with sublayers, due to their reduced LT produced by the B-spline approach. Figure 14 shows the layer representation of the three presented approaches. For the iso-geodesic case, it is also possible to notice how the layer’s closure affects the shape of the B-spline reconstruction, which creates a nonuniform distribution on the lateral surface. The difference in the top part between the iso-geodesic and optimization is also due to the coarsening of the mesh discretization in the iso-geodesic case compared to the numerical generated contour for the optimization case.

Figure 15 shows the theoretical layer distance per each toolpath point. The iso-v method without any adaptive approach shows an LTV which at most arrives at −46% and with the upper limit of 0%, meaning that only over-extrusion phenomena can appear. The iso-v case also shows a monotonic behavior of LTV, as clear from the linear behavior of the mean LTV shown in Figure 13. For optimization and iso-geodesic methods, the LTV behavior is almost the same. Both show the same LTV bounded between 33% and −33%. This is expected because both methods return the same toolpath. Compared to iso-v, the lower limit is higher, reducing over-extrusion probability. The main difference between optimization and iso-geodesic is that the layer closure in iso-geodesic is sharper and noisier as shown in Figure 10. This causes less uniformity in the deposition from the LTV and discontinuity in tool orientation during deposition compared to the optimization case.

For each method, one part is fabricated with the parameters shown in Table 2. Figures 16–18 show the fabricated components for the three cases: iso-v and the two contour methods with adaptive slicing and sublayers. The components are then cut and grinded as explained in Section 4.2. For each method, quarters of torus cross-sections are shown in Figures 19–21. In the Figures are highlighted the real interlayer distance behavior and microscopic-analyzed regions are compared with the theoretical ones. For each sampled region layer, distance measurements are collected with an optical microscope and compared with the expected layer distance from analytical computation. The measurement ranges, mean and standard deviation are summarized for each case in Tables 4–6. The iso-v case LT measures are following the same analytical results for the region analyzed. In regions H-I, over-extrusion happens, making it impossible to collect measures. Figure 19 shows Region I where the cross section of the component is opaque and uniform not allowing to see the differences between the toolpath lines. Those regions expect a LT of 0.4 mm at most. For the other two cases, the over-extrusion effect does not appear for bigger regions. Compared to the iso-v case, looking from the bottom, the LT error is lower up to E-F regions. Besides, the theoretical results mostly match the experimental ones. From a punctual LTV point of view, in the iso-v case, the LT ranges between −46% and 0%. Instead, adaptive methods provide both positive and negative variation. Adaptive error optimization has a maximum layer distance of 1.2 mm (corresponding to an error of +33%) and a minimum of 0.6 mm (corresponding to −33%). Similarly, the iso-geodesic method leads to the same LTV variation, −33% for the minimum and +33% for the maximum. Optimization and iso-geodesic methods reduce the extrema values of LTV. On the other hand, both methods have positive and negative variations providing an overall variation (ΔLTV_max) of 66%. Besides, in adaptive optimization and iso-geodesic cases, the theoretical LT changes adaptively to be on average closer to the desired LT (0.9 mm). Table 7 summarizes the main results obtained from the comparison between slicing algorithms. In particular, summarized, the total number of layers (NL), LTV and ΔLTV_max indicate the maximum internal interlayer distance variation.

This work compares three NP-contour algorithms and proposes an NP-adaptive slicing approach with an infill strategy using the B-spline method. The aim is to understand the LTV behavior on the contour and the infill. To show the capability of the geodesic algorithm with the infill algorithm to process more complex geometry, we present two case studies: a tubular bifurcation (Y-case) and a single chain (O-case). These geometries have been sliced using the iso-geodesic method for the contour. The iso-geodesic allows to be applied on STL object allowing to fabricate CAD and real objects. The adaptive approach with sub-layers was applied to fill the contour. The presented Y-case has been then fabricated showing a general and simple strategy to manage bifurcation similar to (Mitropoulou et al., 2020). The O-case is presented theoretically to discuss the geometry complexity and the possible strategies that can be used.

The Y-case study is bifurcated geometry based on a Sweep feature of circular contour along a path. The proposed algorithm for filling is not able to handle directly bifurcated geometry, due to the impossibility of differentiating the two surfaces of the bifurcation (i.e. the two branches’ sections). Those surfaces should be handled independently. Hence, a volume decomposition strategy is required to handle this case. The component has been cut to have two branches (i.e. a main branch and a secondary branch) using CAD. Then, the main branch (i.e. the one with the base on the base plane) has been processed as previously discussed for the regular toroidal shape. For the secondary branch, it is necessary to identify the vertices on the support face that can be used as the sources for the geodesic field computation. This can be achieved considering that the orientation of triangles’ normal is perpendicular to the cut. From those triangles, the vertices are used as source points for the geodesic scalar field computation. After the field has been computed independently on both branches, the infill can be applied to both. The geodesic field of the object is shown in Figure 22. The two components have two different color bars according to their scalar field. The lines represent the isolines of the geodesic field (i.e. the layers before being processed by infill strategy). The black lines for the main branch, the white ones for the secondary branch. The 3D-printed Y-case is shown in Figure 23. As already shown in Figure 22, the second branch results almost multi-directional. It is printed with an almost constant tool orientation of about 27 degrees concerning the vertical direction. It can be noticed also from the isolines that the volume decomposition and the different slicing can cause anisotropy between the two parts and tool repositioning issues. Anisotropy is more evident due to the fact that a filled volume is fabricated. Those effects can be mitigated with a more suitable volume decomposition strategy, a different deposition strategy that can deposit materials switching between the different branches and, joined with a collision avoidance strategy able to consider global and local collisions (Jayakody et al., 2022).

A more complex geometry is shown using a combination of volume decomposition and the proposed method. In this case, a bifurcation with a joint of the bifurcation is proposed. As in the previous case study, the part needs to be appropriately cut to achieve a single area between the contour to be interpolated. Besides, in this case, it is more detrimental to consider the foreseen tool angles (mainly to avoid collision and considering gravity) and a feasible strategy to deposit material avoiding collisions. Feasible strategies include orienting properly the tool considering the tool geometry (i.e. tool size and nozzle angles) and depositing layers switching between parallel branches. The volume decomposition strategy should have as input the geometric constraints of the tool head. Figure 24 shows the iso-geodesic slicing for the geometry. The white lines are the iso-geodesic contours used for slicing. As for the previous case, each single portion uses an independent geodesic field. The cuts are chosen in the saddle points of the bifurcation and the cuts on the top side are chosen with a small angle to reduce the probability of collision between the tool and the two branches. However, at present, the slicing procedure of such complex shapes requires cutting manually the geometry and analyzing the toolpath in simulation to check the feasibility. Furthermore, it is possible to notice that more branches and joints are present, and the number of cuts increases, increasing the method complexity. To address such complexity, it is required to integrate more complex strategies in the slicing procedure that consider the collision avoidance between the tool head, previous layers and bed plate. Besides, a strategy to fabricate sub-portions of the single branches can be adopted to mitigate collisions and being more flexible in the volume decomposition.

This work provides a comprehensive overview of the main benefits and challenges of NP slicing for filled free-form fabrication. It compares three different contour methods for slicing the quarter of the torus. The contour methods used both an analytical and mesh-based representation as input geometry. The methods have the goal of reducing the LTV. To realize filled fabrication, a B-spline patch method is used for reconstructing the inner side. Then, an adaptive slicing approach is performed for the cases where LTV decreases significantly. Furthermore, a comparison between numerical and experimental cases is provided using the LTV. Finally, two bifurcation case studies have been used to test on a more complicated case the proposed algorithm. The main outcomes from the comparison between results coming from presented NP slicing approaches can be summarized in:

• Best performance in terms of distance between layers contours achieved by iso-geodesic and error optimization slicing algorithms, achieving zero error against 46% in reduction obtained from the iso-v.

• Successful supportless printing is achieved for all strategies.

• Minimum absolute internal interlayer distance error is achieved by the iso-v slicing with a maximum reduction of 46% against 66% obtained by adaptive error optimization and adaptive iso-geodesic.

• Matching between theoretical and real results in terms of internal distance between layers.

• Error optimization method provides similar results to iso-geodesic enabling constant LT on contour for general mesh discretized objects.

• More complex geometries such as bifurcated shapes should be manually handled cutting manually the geometry.

• Process planning is more relevant in the case of more complex geometry.

The authors acknowledge the MSc. student Vincenzo Orlando that during his thesis helped in the design of the path planning, algorithms comparison and experimental campaign.

Declarations: The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Conflict of interests: The authors declare that they have no competing interests.