Convexity and Surface Quality Enhanced Curved Slicing for Support-Free Multi-Axis Fabrication

Abstract

In multi-axis fused deposition modeling (FDM) printing systems, support-free curved layer fabrication is realized by continuous transition of the printer nozzle orientation. However, the ability to print 3D models with complex geometric (e.g., high overhang) and topological (e.g., high genus) features is often restricted by various manufacturability constraints inherent to a curved layer design process. The crux in a multi-axis printing process planning pipeline is to design feasible curved layers and their tool paths that will satisfy both the support-free condition and other manufacturability constraints (e.g., collision-free). In this paper, we propose a volumetric curved layer decomposition method that strives to not only minimize (if not prevent) collision-inducing local shape features of layers, but also enable adaptive layer thickness to comply with a new volumetric error-based surface quality criterion. Our method computes an optimal Radial Basis Functions (RBF) field to modify the fabrication sequence field, from which, the iso-surface layers are extracted to design the corresponding multi-axis printing tool paths. A method to fine-tune variable nozzle orientations on each resulting curved layer is then proposed to overcome possible collisions in high-genus geometries. To validate the concept and exhibit its potential, several support-free fabrication experiments and comparisons with the conventional geodesic field-based slicing are presented, and the results give a preliminary confirmation of the feasibility and advantages of the proposed method.

1. Introduction

In contrast to subtractive manufacturing (i.e., machining), Additive Manufacturing (AM) offers a better means for fabricating complex 3D freeform models [1]. However, the conventional planar-layered material deposition (with a fixed build direction) employed in AM constrains its design space. In Fused Deposition Modelling (FDM), this restricted fabrication nature induces manufacturability constraints such as support structure requirement in overhang regions and degraded surface quality due to the staircase effect [1,2]. Despite numerous efforts [3,4,5] on optimising certain process parameters (e.g., build orientation, layer thickness) in the conventional (3-axis) FDM printer setting, the overall impact of the constraints remains inevitable in 3D models that contain complex geometric (e.g., large overhang, high curvature) and topological (e.g., high-genus) features. The emerging multi-axis AM technology [6,7] has introduced the concept of employing high degrees-of-freedom (DOF) printer systems to enable the printing nozzle to continuously change its orientation during fabrication. Wu et al. [8] exploited this rotational flexibility to realise support-free fabrication of 3D models by printing (with planar layers) segmented parts of a 3D model along different build orientations. Inspired by the potential of support-free printing, several notable efforts [9,10,11] on model-segmentation-based (i.e., the 3+2-axis) multi-axis printing were proposed in recent years. However, due to the limitations inherent to model segmentation techniques (e.g., relying on planar slices, fixed build direction for each segmented part, etc.), a complete elimination of support structures remains impossible (especially for 3D models with both overhanging and high-genus structures) for these so-called 3+2-axis multi-axis printing methods [10,11].

To overcome the limitations in both 3-axis and 3+2-axis multi-axis AM methods, researchers have now endeavoured to develop 3D printing with non-planar printing layers. It is evident that non-planar (i.e., curved) printing layers are effective in minimising the staircase effect in curved 3D models [12,13]. Dai et al. [14] presented a curved layer decomposition method, which employs continuous nozzle orientation change (i.e., continuous multi-axis printing) to realise non-planar print paths to achieve support-free printing. As their method relies on the convexity of decomposed surface layers (to avoid collision during fabrication), it nevertheless suffers from severe surface artifacts and high computational load. Recently, some methods of curved layer decomposition for support-free continuous multi-axis fabrication based on the use of geodesic distance field were introduced [15,16,17]. However, those by Xu et al. [15] and Li et al. [16] showed that highly concave curved layers are inevitable when decomposed based on the geodesic distance field. Thus, local gouging (i.e., the collision between the printer nozzle and the printing layers) remains a serious concern (see Figure 1), restricting their algorithms to considerably simpler 3D models. In contrast to the conventional multi-axis printer setting (i.e., the nozzle orientation remains fixed while the workpiece changes its orientation), Mitropoulou et al. [17] deposited material (with a mobile tool) along non-planar paths on a fixed workpiece to realise support-free fabrication of large-scale single-shell surfaces. However, their method offers no direct solution to prevent local gouging and collisions that could occur between the nozzle and the in-process workpiece. Mainly concerning strengthening a printed part, Fang et al. [18] presented a strength-aware curved layer decomposition method, which, though, did not consider the support-free fabrication requirement.

A recent work by Xie et al. [19] showed how convex ellipsoidal layers could be utilized to overcome the local-gouging problem while upholding the support-free requirement. In fact, some quite complex 3D geometries (e.g., a genus-3 model) were successfully fabricated by their method. However, as the skeleton is used in their method, full covering of a complex 3D geometry with ellipsoidal slices becomes a serious problem (as typically, it is extremely difficult to stably compute the skeleton of a complex 3D model, particularly near the model boundary). In addition, their method suffers from risks of collision during fabrication of highly convex surface layers in proximity (see Figure 2). Recently, a vector-field-based curved slicing method was proposed by Li et al. [20] to realise support-free multi-axis fabrication. Unfortunately, their strategy to eliminate concave (i.e., local gouging-inducing) shape features in curved layers solely depends on minimising the mean curvature, which is an extrinsic property that contains little information of the local shape properties of a given curved surface layer. Therefore, the resultant shapes of most decomposed layers tend to be undesirable (highly planar) rather than being convexity-enhanced. This caveat complicates (if not prevents) maintaining the support-free requirement and, as a result, low printing surface quality becomes a serious concern.

From the above review, it is evident that manufacturability constraints native to curved layer decomposition techniques have caused loss of generality in handling complex geometrical features. For example, determining the severity of concavity in decomposed layers has never been well studied before to evaluate the manufacturability of a given curved surface layer. Even if the concavity were eliminated (i.e., when convexity is guaranteed), constraints such as collisions due to highly convex layers in proximity of one another remains inevitable in the existing methods [14,19]. On the other hand, improvement of the surface quality of printed models has been rarely (if ever) studied in the domain of multi-axis printing.

In this paper, we present a new curved layer decomposition method that employs a multi-objective optimisation model to provide a generalised solution to the above-explained challenges (i.e., evaluating the severity of concave features, collision in high-convex layers, and unawareness of surface quality control). Our method employs the powerful algebraic tool Radial Basis Functions (RBFs) to compute an optimised distance field, which focuses on minimising (if not completely eradicating) unsafe regions within the decomposed curved layers and improving the surface quality based on a volumetric error model, while ensuring the support-free fabrication requirement for a range of complex 3D models. As a result of our effective unsafe shape feature-minimising strategy, we have successfully fabricated some complex 3D models that would fail under conventional geodesic field-based slicing methods [15,16]. We will show that allowing adaptive layer thickness in curve-layer design has potential to improve the surface quality of fabrication. Finally, we propose an experimental approach to fine-tune the nozzle orientation vector on a multi-axis tool path to realise collision-free curved-layer fabrication of complex 3D models that contain both high-genus and large overhanging structures.

2. RBF Field-Based Curved Layer Decomposition

2.1. Overview

Given an input 3D model represented as a tetrahedral mesh $M$ containing $V$ vertices, our goal is to decompose M into a set of curved surface layers that can ensure an accurate approximation of the input model geometry while enabling fabrication of complex shape features (e.g., high-genus structure, large overhang, etc.) in a support-free manner. Therefore, it is evident that our goal demands for multiple objectives to be satisfied during the non-planar layer decomposition stage. Optimisation models of multiple objectives are generally non-linear [14,18]; thus, the computational cost and the complexity of the optimisation algorithm could rise to an exceptional level [18].

Our process-planning pipeline employs the RBF interpolation (Section 2.2) technique to compute a weighted fabrication sequence field (Section 2.3) for the input 3D model. To optimise (i.e., compute the optimal RBF design parameters) the fabrication sequence field, we first introduce several mathematical models (Section 2.5) to evaluate the impact of fabrication constraints and to determine the manufacturability on each extracted iso-surface layer (Section 2.4). Ideally, the optimisation phase iteratively modifies the shape of each problematic (i.e., containing unsafe regions) curved surface layer while respecting the given thresholds (e.g., the overhang angle criterion, the layer thickness, etc.). Finally, a (printing) tool path on each decomposed curved layer is planned to realise continuous multi-axis printing.

2.2. RBF Interpolation

In this section, we give a brief introduction of the RBF interpolation technique to make this paper self-contained. The RBF theory emerged in topographical data representation for the first time [21]. Given a certain metric space (e.g., the Euclidean distance), RBFs have a potential to generate a node-wise response, which mainly relies on the distance between a given node $P \left(P\in V\right)$ and the centre point (RBF centre) of each overlapping RBF field. Since the overall response of overlapping RBF fields is simply formulated as a linear combination, the interpolated scalar value at $P$, $O\left(P\right)$, can be computed by the following equation [21]:

$$O\left(P\right)={\displaystyle \sum }_{k=1}^n{W}_k\mathsf{\Phi }(\parallel P-C_k\parallel )$$

(1)

in which $\mathsf{\Phi } , C_k$ $,W_k$ represent the RBF kernel, the RBF centre location, and the RBF centre weight, respectively, in a domain comprising n RBF centres, while $\parallel .\parallel$ indicates the L1 norm. From Equation (1), it is evident that the RBF interpolated values also depend on the following design parameters,

• The number of RBF centres
• The RBF centre weights
• The RBF centre locations
• The RBF kernel type

2.3. RBF Weighted Fabrication Sequence Field

One of the challenges in approximating a 3D model with non-planar layers is to automatically induce a natural design for the curved layers that are stacked in a sequential order [15,16]. As the first step of our curved layer decomposition method (Section 2.4), we employ the heat geodesics method proposed by Crane et al. [22] to compute an initial heat gradient vector field $H$ (see Figure 3a), which can help induce a geodesic distance field in different kinds of discretized domains (in our case, a volumetric mesh). It is important to note that the base vertices of the input 3D model are set as multiple heat sources for this computation, thus propagating a heat flow from the base to the top of the model. As the result of this step, each tetrahedron of $M$ will be assigned with a unique discrete gradient vector by the end of this step [22].

Next, we compute the RBF weight field $R_X$, which will be used to modify the original heat propagation patterns inside $M$. Ensuring the generality of the function domain of RBF interpolation is crucial. Therefore, for each tetrahedron of $M$, its centroid is considered as a node of the RBF function domain. The coordinates of the centroid can be computed by averaging the sum of the vertex coordinates of each tetrahedron [23]. To guarantee the smoothness of weights spreading over the function domain ($M$) and to preserve their magnitude within an appropriate range, we choose the infinitely smooth Multiquadric (MQ) RBF kernel [21] for the interpolation. The MQ kernel can be represented as

$$\mathsf{\Phi }=\sqrt{1+{\left(\epsilon r\right)}^2}$$

(2)

in which $\epsilon$ and $r$ stand for the shape parameter of the MQ kernel and Euclidean distance between a given RBF centre and each node of the function domain (i.e., $\parallel P-C_k\parallel$), respectively. Even though $\epsilon$ has a slight influence over the interpolated RBF field values (i.e., changing $\epsilon$ will result in a negligible difference in the RBF field values) [21], we observe that it is considerably less impactful than the design parameters listed in Section 3.2. Therefore, keeping $\epsilon$ fixed (in our case $\epsilon =1$), we compute the RBF weights by interpolating the functions defined at RBF centres. The interpolated RBF weight field $R_X$ can be viewed as a point cloud (Figure 3b) that comprises a unique scalar value at every node of the function domain. Finally, to modify the original heat propagation pattern (i.e., modify the original geodesic distance field with the RBF weights), the weighted gradient vector field $H_W$ is derived as the product of the initial gradient vector field and the RBF weight scalar field, as follows:

$$H_W=H\times R_X$$

(3)

To extract the RBF weighted distance field (i.e., the fabrication sequence field) (see Figure 3d), we solve the Poisson equation on $H_W$ to convert it into a scalar field, which assigns a scalar distance value $\phi$ at each $V$ of the mesh by the following equation:

$$\Delta \phi =\nabla H_W$$

(4)

2.4. Iso-Surface Extraction

In this step, we decompose the fabrication sequence field into a sequence of curved surface layers, on each of which, the manufacturability will be evaluated, and a multi-axis tool path will be designed. Given an iso-value F of the fabrication sequence field, we march through every tetrahedron and check for any intersection by the corresponding iso-surface S (i.e., check if iso-value F lies between the extremal fabrication field values at the four vertices of the current tetrahedron). If such intersections exist, the iso-points (V_iso) of S are computed by linearly interpolating the vertex coordinates of the tetrahedron by the following equation:

$$V_{iso}=\frac{F-D_i}{D_j-D_i}\left(V_j-V_i\right)+V_i$$

(5)

in which $V_j \mathrm{and} V_i$ represent the end vertex coordinates of a candidate edge while $D_j \mathrm{and} D_i$ stand for the field-values of $V_j$ and $V_i$ respectively. While there are many possibilities for occurrence of the three intersecting points (see Figure 4), the extremely rare case of having four intersecting points of an iso-surface in one tetrahedron is handled by splitting the quadrilateral into two triangles across the diagonal (see Figure 5). For example, in the first case (Figure 4a) of three intersecting points in Figure 4 (assume iso-value F is given), the coordinates of V_iso can be computed by replacing V_i, V_j, D_i, and D_j of Equation (5) with V_r, V_s, D_r, and D_s, respectively. The exact similar interpolation is performed along V_r, V_q and V_r, V_p edges. In the case of four intersecting points (Figure 5), the above interpolation is performed twice along the corresponding edges of both triangles. Once all these intersection points inside model $M$ are connected, the iso-surface layer that corresponds to iso-value F can be automatically extracted in the form of a triangular mesh. However, the mesh connectivity of the extracted surface layers would mostly remain irregular, leading to numerical degeneracies in the later tool path design phase. Therefore, we apply a remeshing technique proposed by Botsch et al. [24] to ensure the required regular mesh connectivity within all the decomposed curved layers.

2.5. Manufacturability Constraint Evaluation Models

In this section, we introduce mathematical models that evaluate the manufacturability on each decomposed curved surface layer (i.e., to identify whether a given curved layer qualifies for fabrication based on the pre-defined criteria).

2.5.1. High-Risk Local Shape Feature Detection

As is well known, concave features (i.e., local areas where inward pointing vertex normal vectors exist) in non-planar surface layers are the most prominent cause behind local gouging [15,16]. In addition to constraining the complexity of manufacturable 3D geometries, local gouging also prevents safe fabrication on multi-axis printer systems. Our approach is to first identify the existing local concavities and to compute the severity of total local concave features in each curved layer to determine whether the curved layer could be fabricated in a gouging-free manner. For this purpose, we implement the surface shape analysis model as introduced by Koenderinek et al. [25]. Both the Shape Index ($S_V$) and Curvedness Index ($C_V$) measures presented in this model are computed by employing the discrete principal curvatures ($K_1,K_2$) of the surface layer mesh, as follows:

$$S_V=\frac{2}{\pi }ta{n}^{-1}\frac{K_2+K_1}{K_2-K_1}$$

(6)

$$C_V=\sqrt{\frac{K_1{}^2+K_2{}^2}{2}}$$

(7)

$S_V$, which ranges from −1 to 1, distinguishes between vertices located in convex/concave local shape features. Given a set of layers with concave features, we observe that the degree of curvature measured by $C_V$ is an ideal indicator of the severity of local concavities. Since the first objective of our optimisation is to improve the convexity of decomposed layers, it is crucial to integrate both the above measures to define certain criteria that could filter out the unsafe features that are prone to local gouging. Once both $S_V \mathrm{and} C_V$ are evaluated on a given curved surface layer (see Figure 6), it is conceivable that, for local concave features, their $S_V$ is generally negative. Through our tests, it is found that in highly curved regions of such a concave surface layer the value of $C_V$ is typically greater than 5. Therefore, by following an experimental approach, we consider a vertex to be in a high-risk (i.e., local-gouging inducing) region if it fulfils the below criteria:

$$S_V<0.2 and C_V>4$$

(8)

The total number of at-high-risk vertices, $H_K$(which fulfill the above criteria), will serve as an indicator for the optimisation model. For instance, on a totally convex curved surface layer, $H_K$ would be zero. Thus, minimising $H_K$ is set as one of our objectives for optimising the fabrication sequence field to improve the convexity within the decomposed curved surface layers and to avoid local gouging.

2.5.2. Volumetric Error Approximation

In this section, we introduce a new mathematical model that approximates the volumetric error between two adjacent curved surface layers (i.e., volumetric deviation from the original 3D representation). The discretisation effect of both planar and non-planar slicing procedures causes compromised surface quality of the final printed part. While surface quality is directly influenced by the layer thickness distribution, curvedness of surface profile (i.e., surface angle), and the part orientation, numerous adaptive slicing methods (in 3-axis AM) based on the cusp-height measure [26,27] and volumetric error [28,29] have been proposed to uphold the printed surface quality. Volumetric error approximates the deviation of the sliced model and is considered superior to cusp-height for adaptive slicing techniques [29]. However, in the case of a curved layer, surface quality improvement has rarely been studied, mainly due to two reasons: the unpredictable nature of curved layer decomposition and the difficulty in handling numerous objectives to satisfy both the support-free and collision-free requirements.

It is important to note that the geometric deviation occurs (or is measured) along the boundary of each decomposed layer [4,29]. Our volumetric error approximation model starts with finding the boundaries (i.e., the boundary vertices) of each curved surface layer. For this task, all the vertices that form non-repeated edges (i.e., boundary edges) are extracted as boundary vertices. We compute the volumetric error at each boundary vertex since the local layer thickness is adaptively changed in the optimisation phase. If the total number of boundary vertices in adjacent layers $S_K \mathrm{and} S_{K+1}$ are given as $m \mathrm{and} n$, we consider $m (m<n)$ as the total number of corresponding vertices. In our case, it is crucial to form the correspondence between the boundary vertices of two adjacent layers $S_K \mathrm{and} S_{K+1}$ to approximate the triangular cross sections formed by the deviation between the two [29]. To find the corresponding vertex pairs between the boundary vertices $\left(B_{k,j=1,\dots ,m},B_{k+1,j=1,\dots ,m}\right)$ of $S_K \mathrm{and} S_{K+1}$, we first calculate the Euclidean distance from each vertex $B_{k,j}$ of the current layer to all the boundary vertices $B_{k+1,j}$ of the adjacent layer. We label the closest vertex pairs as the corresponding vertices for any two adjacent curved surface layers. The next step is to form triangular cross sections for each corresponding vertex pair.

If the boundary lengths of $S_K \mathrm{and} S_{K+1}$ are given as $R_K$ and $R_{K+1}$ ($R_K$ > $R_{K+1}$), we find the closest point $V_{K,j}$ located on $S_K$ (the longer boundary length layer) to form the triangular cross section for each corresponding vertex pair. Thus, $d=\parallel V_{K,j}-B_{k+1,j}\parallel$ corresponds to the local layer thickness $h$ (i.e., the height of the triangular cross section) that ultimately influences the approximated volumetric error. Once all the cross sections are formed, next, we compute the edges (by considering the Euclidean distance) to connect each cross section. This process forms oblique prisms that approximate the volumetric errors of individual segments between the adjacent curved layers (Figure 7b). The volume of an oblique prism $Vol\left(P_k\right)$ can be easily calculated as:

$$Vol\left(P_k\right)=A\times H$$

(9)

in which $A$ and $H$ $\left(=\frac{l_1+l_2+l_3}{3}\right)$ represent the area of a triangular cross section and the average length of the prism, respectively. Thus, we approximate the total volumetric error $H_V$ between two adjacent curved layers as the total volume of all the surrounding oblique prisms:

$$H_V={\sum }_{k\in m}Vol\left(P_k\right)$$

(10)

It is important to note that the shape and local thickness of curved layers vary (often significantly) in our layer decomposition algorithm, which is one of the significances of the non-uniform weighted heat gradient vector field. By exploiting this ability, our method focuses on optimising the fabrication sequence field in a way that minimises $H_V$.

2.5.3. Boundary Overhang Angle

To achieve support-free fabrication for 3D models with large overhang areas, the angle ${\theta }_B$ between the nozzle orientation vector at each boundary vertex and the surface normal vector of the curved surface layer should be maintained [15] under a pre-determined threshold (in our case, ${\theta }_B<135°$). Since it is clear that the shape of a curved surface layer has a direct impact on the computed overhang angle, any modification to the fabrication sequence field must respect the support-free overhang angle criterion. For this purpose, we evaluate the overhang angles for each decomposed surface layer to ensure that the support-free condition is met throughout. While the nozzle orientation vector computation (Section 3.2) could be modelled in different ways, the surface normal vector at any vertex of any curved layer is considered to be the weighted mean of the face normal vectors of the incident triangles of the layer [30]. Once both the nozzle orientation vector $P_B$ and surface normal vector $R_{B }$ are known, the angle ${\theta }_B$ can be computed by the following equation:

$${\theta }_B=arccos\left[\frac{P_B.R_{B }}{\left|P_B\right|\left|R_{B }\right|}\right]$$

(11)

3. Multi-Objective Optimisation and Multi-Axis Fabrication Path Generation

To find an optimal solution for the multiple objectives in our curved layer decomposition stage, we employ a multi-objective genetic algorithm (GA) optimisation model to search for an effective result. The optimisation variable setting (Section 3.1.1), objective functions (Section 3.1.2) and definition of constraints (Section 3.1.3) in our model are described in this section. Once the optimal fabrication sequence field is computed within the time and other pre-defined constraints, we then design a multi-axis tool path with a variable nozzle orientation vector (Section 3.2) that strives to avoid collisions during the fabrication.

3.1. Multi-Objective GA Model for Field Optimisation

Conceivably, a global optimisation solution (i.e., minimising the total volumetric error or high-risk region vertices of all the decomposed curved layers) would not be a local optimal solution to individual curved layers (e.g., maximizing the convexity). In fact, such a global approach could cause overcompensation in a non-linear optimisation problem. Therefore, to ensure improvement of shape features within each decomposed curved layer, we opt for a local (i.e., layer-based control) optimisation approach (Section 3.1.3).

3.1.1. Optimisation Variable Setting

Since our optimisation solely depends on finding the optimal RBF scalar weights interpolated over the input 3D geometry, the following impactful RBF field design parameters will serve as the optimisation variables:

• Number of RBF centres (N_B)
• RBF centre locations
• RBF centre weights

3.1.2. Objective Functions

The purpose of optimising the fabrication sequence field is to iteratively improve the convexity (i.e., minimising the concavity), while at the same time, minimising the volumetric error within each decomposed curved surface layer through adaptively changing the layer thickness and the shape of surface layers under the support-free overhang angle condition. Thus, the objective functions aim to minimise $H_k$ and $H_V$ at each individual curved layer (Section 2.5). Based on the evaluation score of these two objectives, the fabrication sequence field values are modified until the stopping criteria are met. For a 3D geometry which consists of both concave surface layer features and a high-genus topology (which usually induces high $H_V$), the above two objectives are combined with user-provided weight parameters as $W_1\times H_V+W_2\times H_K$.

3.1.3. Constraints for Optimisation

We set several constraints in optimising the fabrication sequence field to search for a high-quality outcome while considering the computational load and the printer hardware constraints. As explained in Section 2.5.3, the first constraint to ensure support-free fabrication is the overhang angle criterion (${\theta }_B<135°$). Next, the layer thickness range $\left[h_{min} , h_{max}\right]$ is set between the values of 0.1 mm and 0.9 mm to ensure the material extrusion within the printer nozzle size specifications (this depends on the specific hardware). In terms of the RBF interpolation, the maximum number of RBF centres is restricted to 10, which is chiefly due to the concern on computational load. Nevertheless, these parameters can be defined by the user, depending on their requirements and the available processing power.

Finally, we introduce two layer-based constraints, $H_V\le K_1$ and $H_K\le K_2$, to gauge susceptible curved layers with highly concave shape features or causing a high volumetric error. In our case, $K_1=6$ (the maximum number of high-risk region vertices) and $K_2=200$ mm³ (the maximum allowable volumetric error), which are experimentally determined. Finally, to influence the fabrication field values to comply with the above layer-based constraints, we have also included another objective of minimising the total number of erroneous layers $S_E$ and added a penalty score depending on the $S_E$ for each iteration in the optimisation.

3.2. Multi-Axis Tool Path Design

A (printing) tool path for curved layer fabrication on a high-DOF system is usually defined by a position vector, a corresponding nozzle orientation vector, and the accumulated feed value [14,15]. For the position vector that defines the tool path pattern, we adopt the popular parallel-contour pattern. Given a curved surface layer as a triangular mesh, we first extract the boundary vertices (as explained in Section 2.5.2) again and by setting them as the source vertices. The exact surface geodesic field [31] is then computed on the triangular mesh. Then, the intersecting points that form a g number (g = 5 in our case) of iso-contours are computed by linearly interpolating the geodesic distance values assigned at the mesh vertices. Finally, the last point of each contour is connected to the first point of the next contour to form a continuous printing path on each curved layer.

The conventional planning of nozzle orientation to guarantee non-warping fabrication has been to align the nozzle orientation vector along the layer surface normal vector at each vertex of the curved layer. However, sometimes this is susceptible to collisions, regardless of the improved convexity, especially in high-genus 3D geometries with multiple branches where other curved layers exist in the proximity (Figure 8a). To overcome this bottleneck, we propose a new method to design the nozzle orientation vector. Specifically, both the vertex surface normal vector $P_{N }$ and the average surface normal vector (i.e., a fixed vector) $P_{avg}$ at each vertex of a curved surface layer are used to derive the fine-tuned nozzle orientation $P_V$ through the following equation:

$$P_V=\alpha \times P_{avg}+\beta \times P_{N }$$

(12)

in which $\alpha$ and $\beta$ represent the experimentally determined fine-tuning parameters for the input 3D model depending on the geometric complexity (Section 4.4). The motivation of this is that, with a careful tuning, $P_V$ eliminates most of the inward pointing nozzle orientation vectors (Figure 8b) in the proximity, thus preventing (if not eliminating) possible (global) collisions in high-genus 3D geometries (Section 4.4). Note that, while set by the user, the selected $\alpha$ and $\beta$ must ensure that $P_V$ complies with the support-free overhang angle criterion, as defined in Equation (11).

4. Results and Discussion

In this section, the fabrication results generated by our curved layer decomposition method are presented. Our process-planning framework is implemented in MATLAB and the multi-objective optimisation model is also based on the in-built multi-objective genetic algorithm functions of MATLAB. First, we conduct several layer-only fabrication experiments to provide a comparison between our convexity-enhanced layers and the conventional curved layers generated by the geodesic field-based methods [15,16,17]. We show that the layers extracted from the geodesic field often fail due to local gouging, while that generated by our RBF-based field are rid of local gouging. Furthermore, it is important to note that all the selected models in our experiments are too complex to realise a full-model printing (in a safe no-local-gauging and collision-free manner) under the conventional geodesic-field-based curved slicing method and the conventional nozzle orientation settings. Finally, to demonstrate the potential of our method in surface quality enhancement through the proposed volumetric error-based adaptive slicing strategy, simulation results are presented and comparisons are made with both the planar slicing and geodesic field-based slicing methods.

4.1. The Fabrication System

The (homebuilt) multi-axis printer system (Figure 9) employed for our experiments comprises a fixed nozzle with an $n_D$ = 1 mm diameter and a mobile workpiece attached to the end effector of a 6-axis UR5 robot arm. Thus, the workpiece orientation is continuously altered according to the nozzle orientation vectors of a tool path. Meanwhile, the nozzle deposits fused filament (1.75 mm standard PLA) at the position vector coordinates of the tool path with a variable accumulated feed value $(f_p)$ to realise adaptive layer thickness control, so to fabricate decomposed curved layers from our method. The varying rate $f_p$ is computed by the following equation:

$$f_p=\frac{\mu \times \left(l_e\times d_n\times h\right)}{\pi r_m{}^2}$$

(13)

where $l_e,d_n, \mathrm{and} h$ represent the Euclidean distance between any two adjacent points of the tool path, the nozzle diameter, and the local layer thickness, respectively; $r_m$ stands for the filament radius (fixed); while $\mu$ (in our case $\mu$ = 0.7) is a pre-assigned shape parameter that improves the accuracy of $f_p$ by complying to the elliptic shape of the deposited material.

It is important to emphasize that the main purpose of our inexpensive homebuilt multi-axis printer is to assess the feasibility of the developed multi-step printing process, which suffers from a few hardware limitations that, though, could be improved. For example, the positioning error of the UR5 robot arm [14], which causes a considerable impact on the printing surface quality, can be significantly reduced by employing a high-precision 5-axis CNC table [18]. In addition, a higher quality nozzle (e.g., 0.2 mm radius nozzle instead of the current 1mm radius nozzle) would be an ideal candidate to showcase the true potential of our surface-quality-improvement method.

4.2. Simulation Results

In this section, we provide the simulation results for the selected 3D geometries (Figure 10). The input 3D models are decomposed into the traditional planar layers, the conventional geodesic field-based curved layers, and the RBF-based (our method) curved layers. From the test data given in

Table 1. Comparison of simulation statistics for different layer decomposition methods (V_O—number of overhang vertices, V_U—number of high-risk region vertices, VE—total volumetric error in mm³).

3D Model	Number of Tetrahedrons	Layer Decomposition Method
		Planar			Geodesic			RBF-Based (Our Method)
		VO	VU	VE	VO	VU	VE	VO	VU	VE
Kitten	100,045	435	-	5.89	115	1539	3.42	86	900	2.21
2-genus	105,911	580	-	6.2	302	1150	4.3	166	592	3.8
Nintendo	135,913	996	-	7.85	502	1288	4.71	382	992	4.02
3-genus	138,746	289	-	4.12	153	1473	1.8	81	837	1.05

, it is evident that the curved layers generated from our RBF-based method exhibit much enhanced convexity (i.e., lower number of high-risk region vertices—V_U) and contain the lowest approximated volumetric error (VE), while at the same time, achieving a higher compliance to the support-free requirement in contrast to geodesic field-based method.

To visualise the significance of convexity-improvement of our method, we compare the extracted curved layers at some critical regions (i.e., concavity inducing) from both our RBF-based fabrication sequence field and the geodesic distance field for the Kitten and two-genus model. For the Kitten model, the $S_V$ and $C_V$ plots of the layer shown in Figure 11 indicate high concave features (max $C_V$ of 8), posing a risk of local gouging. In contrast, the curved layer from our method, as shown in Figure 12, achieves better convexity (the max $C_V$ is reduced to 5). Figure 13 and Figure 14 depict the similar comparison plots for the two-genus model. Again, by examining the $S_V$ and $C_V$ plots in Figure 13 and Figure 14 respectively, our RBF-based method fares much better in terms of these two important metrics. On a different note, as revealed in Figure 12 and Figure 14, our method naturally generates adaptive layer thickness depending on the local features and volumetric error measure, in contrast to the uniform layer thickness in the geodesic distance-based field method [15].

Naturally, as the number of RBF centres N_B directly controls the span of design space in our optimization model, it greatly impacts the quality of the optimised output.

Table 2. Test result of the relationship between the optimality and the number of RBF centres N_B.

3D Model	Number of RBF Centres (NB)
	2–3			4–6			8–10
	VO	VU	VE	VO	VU	VE	VO	VU	VE
Kitten	102	979	2.32	98	941	2.29	86	900	2.11
2-genus	252	744	5.9	210	682	5.06	166	592	3.8
Nintendo	411	1019	6.85	382	992	4.12	371	921	4.08
3-genus	123	993	1.58	112	941	1.45	81	837	1.05

provides the test data of a comparison experiment of our method with three different ranges of N_B, which shows a clear monotone positive trend between the optimality and N_B. However, the mitigating factors must be considered when selecting a larger N_B, e.g., the increasing computational load, the numerical stability, and the point of diminishing returns.

4.3. Layer-Only Comparison with the Geodesic Field-Based Slicing

The above simulation results indicate that the curved layers decomposed from the geodesic distance field are susceptible to a greater risk for local gouging in all our tested 3D models. To better exhibit the local gouging problem in the geodesic field-based curved layers and the advantages of our convexity-enhanced curved layers, we have fabricated the layers shown in Figure 11, Figure 12, Figure 13 and Figure 14 from both the geodesic distance field and our optimised RBF-based field. Figure 15 shows the fabrication comparison for the Kitten model. The geodesic field-based layer of the Kitten model from Figure 11 is confirmed to incur local-gouging (Figure 15b,c) and therefore, the fabrication could be considered as a failure (as the nozzle must bypass the location to avoid damage to the nozzle). Figure 15a displays the fabrication result of the convexity-enhanced layer (Figure 12) from our method, which incurs no local-gouging. Similarly, for the two-genus model, Figure 16b,c reveal the failed fabrication of the curved layer (Figure 13) by the geodesic distance field, while Figure 16a shows the success by our method on the curved layer (Figure 14).

From Figure 15 and Figure 16, it is evident that in the presence of high concave local features, the quality of the printed part would deteriorate if printing is allowed to continue, regardless of local gouging, as cracks or defects would result on the layer. In addition, it is also highly unsafe for the nozzle. Therefore, local gouging must be prevented.

4.4. Physically Fabricated 3D Models

Finally, we show the physically fabricated models by our RBF-based multi-axis printing process planning system. The nozzle orientation tuning parameters ($\alpha$ and $\beta$) and the total fabrication time (in hours) for each model are tabulated in

Table 3. Comparison of the impact on the number of RBF centres employed for optimal gradient vector field computation.

3D Models	Nozzle Orientation Vector Tuning Parameters		Fabrication Time (Hours)
	$\mathit{\alpha }$	$\mathit{\beta }$
Kitten	0.35	0.65	7.1
2-genus	0.61	0.39	3.9
Nintendo	0.4	0.6	7.8
3-genus	0.42	0.58	4.7

, while Figure 17, Figure 18, Figure 19 and Figure 20, respectively, show the photos of the four fabricated models, i.e., the Kitten, two-genus, Nintendo, and three-genus. Figure 21 depicts the support-structures required to print the overhang areas of the selected models in a 3-axis setting (i.e., with planar slicing).

For the Kitten model, Figure 17a,b illustrate how local gouging is effectively avoided on convexity-improved curved layers during the fabrication process, while Figure 17c displays the effect of varying nozzle orientation parameters towards fabricating highly convex (potential risk of local gouging with a fixed nozzle orientation) layers, and Figure 17d displays the successfully fabricated Kitten model.

In terms of dealing with high-genus topology structures, Figure 18a and Figure 20c indicate the significance of $\alpha$ and $\beta$ (the nozzle orientation tuning parameters) in realising safe fabrication of curved layers of 2-genus and 3-genus models. Fabricating such layers would be impossible with a vertex normal vector-based nozzle orientation. Figure 18b depicts how the convexity-improved layers of our method eliminate local gouging to realise successful fabrication of the slanted two-genus model (see Figure 18d). With the similar convexity-enhancement (e.g., Figure 20a) and adaptive thickness control through our fully automatic process planning, the three-genus model (see Figure 20d) is successfully fabricated in a support-free manner while maintaining a desirable finish surface quality (i.e., reducing the stair-case effect by adjusting the nozzle orientation when printing the part boundary). It is worthwhile to note that Xie’s method [19] relied on manual processing to overcome the complexity of approximating a three-genus model with curved layers.

On the other hand, the Nintendo model contains a significant amount of overhang areas (see Figure 21d); thus, it poses additional challenges for the optimisation phase and the fine-tuning process of the nozzle orientation. However, as depicted in Figure 19b, the tuned nozzle orientation vectors by our algorithm have clearly handled the highly convex layers (in proximity) without collision and achieved a continuous fabrication of the entire model (Figure 19d). In summary, our proposed RBF-based curved layer decomposition algorithm has effectively upheld the numerous manufacturability constraints in multi-axis printing and successfully fabricated the four selected models containing various complex geometric features, while achieving high surface quality (i.e., reducing the stair-case effect by adjusting the nozzle orientation when printing the part boundary surface).

4.5. Limitations

Despite the promising fabrication results (see Figure 22), there are a few limitations in our RBF-based volumetric curved layer decomposition system. First, though not the focus of this work, the hardware (i.e., the homebuilt printer) used in the physical printing experiments is quite crude and inaccurate. E.g., it is even not equipped with a cooling fan to prevent material dragging during the fabrication (as evidently seen in Figure 19a). Algorithm-wise, the optimisation model is currently unaware of generating curved layers that conform to the shape of the input model’s starting base. As a result, when curved layers are decomposed from the bottom (i.e., the base) of the model, there is a slight distortion (e.g., the two-genus model) in contrast to the original 3D model representation (see Figure 14). Additionally, due to the non-linear mapping in the conversion of the nozzle orientation vectors from the workpiece coordinate system into the corresponding joint angles of the robot arm (the machine coordinate system), the issue of material under-extrusion becomes inevitable in some regions of the fabricated models (e.g., the two-genus model). While this issue also affects the surface quality, it could be potentially eliminated by the singularity-aware motion planning method proposed by Zhang et al. [32]. In addition, currently our optimisation model is not focused on improving the mechanical strength properties of decomposed layers. Although formulating the mechanical strength properties as objectives (e.g., see Fang et al. [18]) could be an immediate solution, its impact on maintaining the support-free fabrication requirement for all the decomposed curved layers is yet unknown and the increased computational complexity of the optimisation algorithm is also a serious concern.

5. Conclusions and Future Work

The challenging nature of decomposing a 3D model into printable curved layers (i.e., support-free, no local gouging, and global collision-free) in multi-axis printing restricts the achievable geometric complexity. In this paper, we introduce a computational framework that employs the RBF interpolation technique to compute an optimised fabrication sequence field that is focused on improving the convexity of decomposed curved layers and the resulting surface quality under the support-free fabrication requirement. By evaluating the impact of manufacturability constraints on each decomposed layer using several fitting mathematical models, we identify the layers that require further shape modification to comply with the pre-defined standards. Finally, based on the evaluation result, the GA optimisation model iteratively optimises the shapes of the curved layers by changing the positions and the weights of the RBF centres inside the model, with an aim of realising support-free and gouging-free curved layer fabrication. While our process planning is completely autonomous, the optimisation parameters are also flexible to be set as user defined inputs; thus, users can be given the freedom to select the desired accuracy and the corresponding computational load. We also present a method to fine-tune the nozzle orientation vectors to avoid possible collisions in highly curved (both convex and concave) areas of the decomposed layers.

Our process-planning pipeline presents a new path towards achieving safer support-free curved layer multi-axis fabrication. The physical fabrication results of our curved layer decomposition method are promising as large overhang areas, concave shape features and highly convex layers in proximity of various complex 3D models have been printed without global collision or local gouging. For the future work, our aim is to integrate a global 3D shape analysis (in addition to the current method of analysing local geometric properties of curved layers) into the process-planning pipeline to automatically optimize the nozzle orientations that are aware of the topology of the input model. This would not only eliminate the current experimental approach in computing nozzle orientations, but also increase both the computing efficiency and quality of process planning. Finally, to further improve the printing quality, we hope to include more objectives (e.g., the printing time, the mechanical strength properties) into the current multi-objective optimisation model under the support-free multi-axis printing framework.