Generation of Efficient Iso-Planar Printing Path for Multi-Axis FDM Printing

The emerging multi-axis fused deposition modeling (FDM) printing process is a powerful technology for fabricating complicated 3D models that otherwise would require extensive support structures or suffer the severe stair-case effect if printed on a conventional three-axis FDM printer. However, because of the addition of two rotary axes which enables the printing nozzle to change its orientation continuously, and the fact that the printing layer is now curved, determining how a nozzle printing path to cover the layer becomes a non-trivial issue, since the rotary axes of the printer in general have a much worse kinematic capacity than the linear axes. In this paper, specifically targeting robotic printing, we first propose an efficiency indicator called the material deposition rate which considers both the local geometry of the layer surface and the kinematic capacities of the printer. By maximizing this indicator globally, a best drive plane direction is found, and then the classic iso-planar method is adopted to generate the printing path for the layer, which not only upholds the specified printing quality but also strives to maximize the kinematic capacities of the printer to minimize the total printing time. Preliminary experiments in both computer simulation and physical printing are carried out and the results give a positive confirmation on the proposed method.


Introduction
Fused deposition modeling (FDM) is one of the most popular types of additive manufacturing (AM), wherein a nozzle mechanically extrudes molten thermoplastic or metal materials layer by layer to generate three-dimensional surfaces or objects. Owing to its simplicity and flexibility, FDM printing has been widely used in various fields such as mechanical engineering, aerospace, and even biological engineering [1]. Traditionally, three-axis AM systems with three degree-of-freedom (DOF) are used for printing simplefeatured objects, in which the part is accumulated in planar layers along a fixed nozzle orientation. This feature leads to an easy-handle fabrication process, while it also brings problems such as the stair-case effect, the need for support for overhang features, weak mechanical strength, and restriction of printing motions [2]. The advent multi-axis AM provides a considerably significant solution to resolve these problems. In general, for the multi-axis AM with additional DOFs, the nozzle axis does not have to align with the z-axis, and neither do adjacent printing paths have to lie in a same plane [3]. By carefully adjusting the nozzle orientation with respect to the surface normal, the defects such as the stair-case effect, local gouging problem, and the need for support that existed in three-axis printing can be eliminated, at least in theory.
Research in FDM AM is primarily focused on the following three principal issues: (1) the deposited materials, (2) the mechanism of hardware, and (3) the tool path planning and process control. Among the three, tool path planning/process control can be regarded details of our iso-planar printing path generation method based on the propo tor. The results of the experiments and discussion are reported in Section 4, f the conclusion in Section 5.

Material Deposition Rate
Similar to five-axis milling, when planning a multi-axis printing path, both location (NL) curves and the associated nozzle orientation need to be prop mined. In this section, we described in detail how to determine the effecti width and the kinematic-constrained feed rate, which subsequently leads to th indicator to evaluate the printing efficiency, i.e., the material deposition rate a of the feed direction .

Effective Printing Width
As shown in Figure 2, numerically, a single NL curve (e.g., ( )) identif zle movement during a printing process, which typically is represented as a se ple points { 1 , 2 , … , } with respect to the workpiece coordinate system (WC the workpiece (i.e., the layer surface S) is defined. At the same time, a standard --is naturally defined at each point on ( ) where is an angular that determines the drive plane of this NL curve (see Section 2.3), while and feed direction and normal direction of the surface at , respectively, and product of the former two.

Material Deposition Rate
Similar to five-axis milling, when planning a multi-axis printing path, both the nozzle location (NL) curves and the associated nozzle orientation need to be properly determined. In this section, we described in detail how to determine the effective printing width and the kinematic-constrained feed rate, which subsequently leads to the proposed indicator to evaluate the printing efficiency, i.e., the material deposition rate as a function of the feed direction θ.

Effective Printing Width
As shown in Figure 2, numerically, a single NL curve (e.g., C P (θ)) identifies the nozzle movement during a printing process, which typically is represented as a series of sample points {p 1 , p 2 , . . . , p m } with respect to the workpiece coordinate system (WCS) wherein the workpiece (i.e., the layer surface S) is defined. At the same time, a standard local frame f p -n p -k p is naturally defined at each point p on C P (θ) where θ is an angular parameter that determines the drive plane of this NL curve (see Section 2.3), while f p and n p are the feed direction and normal direction of the surface at p, respectively, and k p is the cross product of the former two. In the realm of FDM printing, some previous research [16][17][18][19] studied various acteristics of the FDM process, and many mathematical models were proposed f modeling of the geometry of the deposed material. Of all the proposed methods for eling the cross-sectional shape of the deposited filaments, the ellipse model, whic originally proposed by Ahn et al. [17] and subsequently extended by Angelo1 et a with experimental validation, is commonly adopted. Therefore, the ellipse model adopted by us in this paper, which represents the cross-section of adjacent NL cur two intersected ellipses. It should be mentioned that, for a genuine curved layer concave and convex cases should be assumed, thus the two different scenarios of ad filaments as shown in Figure 3. Considering that the modeling ellipse is very small the printing layer is sufficiently smooth, the curve intersected by theplane of and the layer surface can be approximated as an arc passing through . The sam proximation can also be applied to the intersection with the previous layer. Thu cross-sectional ellipse can be parametrically formulated in theplane as: where is the angular parameter measured from the -axis, with and ℎ being th jor-axis length and minor-axis length, respectively. In this paper, we assume that s is to the nozzle diameter and h is the layer thickness at p. In the realm of FDM printing, some previous research [16][17][18][19] studied various characteristics of the FDM process, and many mathematical models were proposed for the modeling of the geometry of the deposed material. Of all the proposed methods for modeling the cross-sectional shape of the deposited filaments, the ellipse model, which was originally proposed by Ahn et al. [17] and subsequently extended by Angelo1 et al. [19] with experimental validation, is commonly adopted. Therefore, the ellipse model is also adopted by us in this paper, which represents the cross-section of adjacent NL curves by two intersected ellipses. It should be mentioned that, for a genuine curved layer, both concave and convex cases should be assumed, thus the two different scenarios of adjacent filaments as shown in Figure 3. Considering that the modeling ellipse is very small while the printing layer is sufficiently smooth, the curve intersected by the n p -k p plane of point p and the layer surface can be approximated as an arc passing through p. The same approximation can also be applied to the intersection with the previous layer. Thus, the cross-sectional ellipse can be parametrically formulated in the n p -k p plane as: where ρ is the angular parameter measured from the k p -axis, with s and h being the majoraxis length and minor-axis length, respectively. In this paper, we assume that s is equal to the nozzle diameter and h is the layer thickness at p.  Figure 2. A printing nozzle location (NL) curve defined in WCS and its local frame.
In the realm of FDM printing, some previous research [16][17][18][19] studied various characteristics of the FDM process, and many mathematical models were proposed for the modeling of the geometry of the deposed material. Of all the proposed methods for modeling the cross-sectional shape of the deposited filaments, the ellipse model, which was originally proposed by Ahn et al. [17] and subsequently extended by Angelo1 et al. [19] with experimental validation, is commonly adopted. Therefore, the ellipse model is also adopted by us in this paper, which represents the cross-section of adjacent NL curves by two intersected ellipses. It should be mentioned that, for a genuine curved layer, both concave and convex cases should be assumed, thus the two different scenarios of adjacent filaments as shown in Figure 3. Considering that the modeling ellipse is very small while the printing layer is sufficiently smooth, the curve intersected by theplane of point and the layer surface can be approximated as an arc passing through . The same approximation can also be applied to the intersection with the previous layer. Thus, the cross-sectional ellipse can be parametrically formulated in theplane as: where is the angular parameter measured from the -axis, with and ℎ being the major-axis length and minor-axis length, respectively. In this paper, we assume that s is equal to the nozzle diameter and h is the layer thickness at p. Referring to [11,15], to determine how far away the two adjacent printing paths should be placed, we recall that the cusp height is commonly used for gauging the surface accuracy in CNC machining. However, considering the unique characteristics of FDM, we instead opt to adopt a cusp-area criterion for a better measurement of printing quality, based on the observation that, when the melted material is extruded, it will be squeezed as its shape is bounded by its neighbor's re-solidified filaments. Therefore, there is a cusp-area between the two adjacent ellipses, which is denoted as the overlapping area A o (p) (i.e., the region shaded in black in Figure 3). Besides, there are also two enclosed gap areas A g (p) (the shaded in blue in Figure 3) that are sandwiched between the adjacent ellipses and the previous layer surface as well as between the ellipses and the current layer surface, respectively. As such, the cusp-area equilibrium criterion-A g (p) = A o (p)-is thus proposed to help determine the effective printing width w e f f which is the distance between the two intersection points between an ellipse and its two neighboring ellipses. The intersection position satisfying the cusp-area equilibrium criterion can be numerically obtained by applying the Gauss-Green formula to the segmented area of ellipses [20], which will be designated by the angle ρ * (see Figure 3). By combining with Equation (1), the effective printing width w e f f at point p can then be expressed as: It can be seen that w e f f is highly related to the local curvature at point p in the n p -k p plane; along a different feed direction θ, the value of w e f f will change accordingly to satisfy the cusp-area equilibrium criterion. This implies that w e f f is a function of θ, which can be written as w e f f (θ).

Kinematic-Constrained Feed Rate
To better determine the nozzle feed rate for upholding the printing quality, not only should the printing path be carefully designed, but so should the motion of the robotic manipulator. During the transformation from a printing task defined in the Cartesian space to the configuration space of the robotic manipulator, the kinematic redundancy can be solved by taking the printing task as the priority over the degree redundancy of the robotic manipulator. In this sub-section, we will first introduce the classic kinematic redundancy problem and the singularity avoidance problem, and then describe in detail how to optimize the nozzle feed rate of printing.

Kinematic Redundancy and Singularity Avoidance
As the self-rotation of the nozzle-axis is irrelevant to printing, a general printing task normally requires five DOFs with 3-DOF for positioning and 2-DOF for identifying the orientation of the nozzle. As a result, when an industrial 6-DOF robotic manipulator is applied for printing, such as the frequently used UR5 manipulator shown in Figure 4, the extra DOF will result in a functional redundancy in robotic motion planning (for more details about the kinematic redundancy, please refer to [21][22][23]).
By default, the printing path is firstly defined in the WCS. According to the geometry information of the workpiece, each NL point on the printing path can be expressed as a vector (x w , y w , z w , i, j, k) defining its point location and the associated orientation. In addition, another two frames are introduced to specify the coordinate transformation: the nozzle coordinate system (NCS) and the base coordinate system (BCS) of the manipulator. As shown in Figure 4, the NCS is defined at the tip of the nozzle, while the WCS is fixed on the end-effector of the manipulator. Since the central point of the nozzle is taken as the origin of NCS and its z n -axis is set exactly coincident with the NL orientation (see the auxiliary graph in Figure 4), the NCS can thus be expressed as (x n , y n , z n , α, β, γ) where the Euler angles α, β, and γ indicate the axle rotation with respect to the WCS. The transformation from WCS to NCS can be identified by a homogeneous transformation matrix T W N as: where trans() and rot() are the standard translation and rotation matrices, respectively. Thanks to the consistency of nozzle orientation defined in the NCS and WCS, a simplified relation between (i, j, k) and (α, β, γ) can be extracted from Equation (3) as:  By default, the printing path is firstly defined in the WCS. According to the geometry information of the workpiece, each NL point on the printing path can be expressed as a vector ( , , , i, j, k) defining its point location and the associated orientation. In addition, another two frames are introduced to specify the coordinate transformation: the nozzle coordinate system (NCS) and the base coordinate system (BCS) of the manipulator. As shown in Figure 4, the NCS is defined at the tip of the nozzle, while the WCS is fixed on the end-effector of the manipulator. Since the central point of the nozzle is taken as the origin of NCS and its -axis is set exactly coincident with the NL orientation (see the auxiliary graph in Figure 4), the NCS can thus be expressed as ( , , , α, β, γ) where the Euler angles α, β, and γ indicate the axle rotation with respect to the WCS. The transformation from WCS to NCS can be identified by a homogeneous transformation matrix as: where () and () are the standard translation and rotation matrices, respectively. Thanks to the consistency of nozzle orientation defined in the NCS and WCS, a simplified relation between (i, j, k) and (α, β, γ) can be extracted from Equation (3) as: This relation reveals the fact that the angle parameter ∈ [− , ] in NCS can be arbitrary, which is the cause of kinematic redundancy.
To facilitate the kinematic analysis of robot manipulators, the configuration Θ is used to denote the joint variables of the manipulator with respect to BCS. Now that the redundant implies infinite joint solutions for the end-effector posture, it can be capitalized for a better planning of the joint configuration under certain specific constraints, e.g., to avoid the singularities. As it is well known in robotic path planning, a good path should not only bypass any singularities but also preferably stay away from them, as near a singular configuration, even a small change in the operational space on the printing path would cause a drastic change in the joint motion in the configuration space in order to maintain the desired trajectory [16]. Therefore, singularity avoidance is chosen by us as one objec- This relation reveals the fact that the angle parameter γ ∈ [−π, π] in NCS can be arbitrary, which is the cause of kinematic redundancy.
To facilitate the kinematic analysis of robot manipulators, the configuration Θ is used to denote the joint variables of the manipulator with respect to BCS. Now that the redundant γ implies infinite joint solutions for the end-effector posture, it can be capitalized for a better planning of the joint configuration under certain specific constraints, e.g., to avoid the singularities. As it is well known in robotic path planning, a good path should not only bypass any singularities but also preferably stay away from them, as near a singular configuration, even a small change in the operational space on the printing path would cause a drastic change in the joint motion in the configuration space in order to maintain the desired trajectory [16]. Therefore, singularity avoidance is chosen by us as one objective to optimize the kinematic performance. To evaluate a singular state of manipulator configuration, Huo and Baron [24] proposed a manipulation performance index w PS , called the parameter of singularity (PS), which is a function of angle γ and defined as: where scalars σ 1 , σ 2 , . . . , σ m are the eigenvalues of the geometric Jacobian matrix sorted in the descending order (more details can be found in [24]). As w PS varies with angle γ, different configurations will lead to different w PS 's, and the smaller it is, the better conditioning is the manipulator. Therefore, this measure is chosen by us as the criterion for evaluating and, thus, to effectively avoid the singularities. For all possible configurations, the optimization objective aims at minimizing w PS to physically place the manipulator as far away as possible from any singular state. Specifically, for a single nozzle posture (x w , y w , z w , i, j, k) defined by its position and orientation, a minimized singularity pa-rameter w PS, min can be resolved by solving the minimization problem of the following form: With γ * identifying the posture of the end-effector (x n , y n , z n , α, β, γ * ), we can then finally determine the joint configuration Θ * expressed as a set of six joint angles (q 1 , q 2 , q 3 , q 4 , q 5 , q 6 ) by the inverse kinematic transformation (IKT) which is a standard computation procedure in robotics [25].

Kinematic-Constrained Feed Rate Optimization
For the robotic arm, it is physically constrained by its configuration limits, such as the velocity, acceleration and jerk, which are also called the kinematic capacities. To make the robotic arm operate within its kinematic capacities, motion planning considering the kinematic constraints is important, as a stable end-effector with less vibration is crucial for maintaining a smooth and steady nozzle path. For this purpose, the feed rate, i.e., the sweeping speed of the nozzle tip, should be adjusted/controlled carefully to assure that the motion of robotic manipulator is always within its kinematic capacities.
Specifically, to calculate the joint velocity, acceleration and jerk at a single NL point that lies on an NL curve, the continuous feed rate scheduling problem is first converted to a discrete form by sampling the NL curve into equal-spaced sample points. The arc length between each pair of adjacent sample points is approximated by the chord-length l c , which is assumed to be sufficiently small. As aforementioned, for any sample NL point p, its corresponding configuration can be determined by choosing from appropriate angle parameters according to the stipulated secondary criterion: the singularity avoidance. Let p 0 and { p 1 , p 2 } be the preceding and succeeding sample points of p, respectively, and Θ 0 , Θ, Θ 1 , and Θ 2 their corresponding joint configurations, with each configuration represented by six joint angles (q 1 , q 2 , q 3 , q 4 , q 5 , q 6 ). Referring to [26], for each point p, its velocity, acceleration and jerk of those corresponding six joints can be numerically derived by: where f (θ) is the feed rate at point p, which is a function of feed direction θ with respect to the local geometry, and l c denotes the chord-length. To satisfy the physical constraints, each joint should be constrained by its maximum values: where v j,max , a j,max , and j j,max denote respectively the maximum (angular) velocity, acceleration, and jerk of joint j, all of which are constants. Combining Equation (7) with Equation (8), the maximized feed rate under the velocity, acceleration and jerk constraints along feed direction θ i can then be identified, which will be denoted as f max, v (θ i ), f max, a (θ i ), and f max, j (θ i ) respectively. Naturally, the most restrictive of the above three should be regarded as the maximum allowed feed rate at p in the θ i direction, i.e., f * (θ i ): By now, the most conservative feed rate under the physical constraints of the robotic manipulator considering the singularity-avoidance has been established. Furthermore, combining with the effective printing width constraint (Section 2.1), we next proposed an indicator of printing efficiency that will synthesize on both these two feed-direction sensitive functions.

Indicator-Determined Preferable Feed Direction
From the perspective of printing, as the surface layer to be printed is given with a specific geometry and thickness, it naturally leads to a fixed volume to be filled by the deposited material. Then, if the volume deposited per unit time is larger, a shorter time will be needed to infill this volume. Since the cross-section of a filament is regarded as an ellipse, the printed volume along a single printing path is thus that of an elliptic cylinder. Considering this fact, an indicator named the material deposition rate v e (θ) is proposed here to evaluate the printing efficiency, which represents the printed volume per unit time. Explicitly, v e (θ i ) at any NL point is defined as: where h is the thickness of the current printing layer, w e f f (θ i ) is the effective width, and f * (θ i ) is the most conservative feed rate determined by the physical constraints and singularity avoidance along the feed direction θ i .
To construct the indicators of the sample points on the layer surface S and find out the representative optimal feed directions for the whole surface, we propose a samplingbased strategy. For any sample point on S, the potential feed direction in the range of [0, π] is first discretized into n equal interval radians, i.e., O = {θ 1 = 0, θ 2 , . . . , θ n = π}. Assuming that S is given as a triangulated mesh, its vertices are chosen as the sampling points for evaluation by the proposed indicator and collectively called the sampling point set V = {s 1 , s 2 , . . . , s m }. At each sample point s i , the effective printing width w e f f (θ i ) and the optimized feed rate f * (θ i ) for any θ i ∈ [0, π] can be found as described in the previous sections. Then, by Equation (10), the material deposition rate on S is uniquely defined as a continuous scalar function of the feed direction θ i at point s i . However, although there is an optimal feed direction θ i * and the corresponding maximum v e (θ i * ) for each s i , our objective is to find a single feed direction (as associated with a single drive plane) that will maximize the total material deposition rate for the entire S. For this purpose, the following discrete optimization formula is proposed: where 1 m ∑ m j=1 v e (θ i ) denotes the average material deposition rate at θ i among m sampling points. The single solution to Equation (11), i.e., θ * , will then be taken as the preferred printing direction for the entire layer surface S. Below, Algorithm 1 gives the pseudo-codes of this procedure for finding angle θ * .
With the drive plane direction found (as identified by the single angle θ * from Equation (11)), we next proceed to describe how to generate an iso-planar nozzle path to print the entire S.  (6)  4 Θ * ← IKT(x n , y n , z n , α, β, γ * ) referring to [25] (10) 8 θ * ← max 1/m(∑ m j=1 v e,j (θ i )) by Equation (11) 9 return θ * Figure 5 gives an illustration of NL curve generation based on the iso-planar method. Basically, the NL curves C P,i (θ * ) are the intersections between the layer surface S and a series of parallel driving planes. For each drive plane Q, it is perpendicular to the plane x w o w y w and with a normal (−sin(θ * ), cos(θ * ), 0) in x w y w z w .

Optimized Iso-Planar Path Generation
printing direction for the entire layer surface S. Below, Algorithm 1 gives the pseudocodes of this procedure for finding angle * .

Algorithm 1 Optimized feed direction determination
Input: A triangulated mesh surface with vertices set Output: Optimized feed angle * //* Refer to Sections 2.1 and 2.2 1 for ∀ ∈ and ∀ ∈ [0, ] do 2 ( ) ← ℎ( , , , ℎ * ) by Equation (2) 3 * ← ( * ) ← ( , , ( )) by Equation (6)  } by Equation (11)  9 return * With the drive plane direction found (as identified by the single angle * from Equation (11)), we next proceed to describe how to generate an iso-planar nozzle path to print the entire S.  Let Q 1 , Q 2 . . . Q k denote these yet-to-be-determined successive drive planes, each inducing an NL path, e.g., C P,i (θ * ) in Figure 5, which is sampled into a discrete form. The very first drive plane Q 1 is readily defined-it is simply an extreme plane of S with a normal vector of (−sin(θ * ), cos(θ * ), 0). The approximation as introduced in [6] provides a way to determine the sample density of an NL curve, i.e., the forward step l f between any two adjacent points p i and p i+1 is approximated by the chord length satisfying the chord-error constraint δ * c , which can be explicitly expressed as:

Optimized Iso-Planar Path Generation
where r f is the curvature radius of S along the direction of θ * at point p i and δ * c is the specified maximal chord-error.
For the feed-material flow, it is achieved by the mechanics of the motor controllers of the 3D printer. In our experiment setting, the motor control of the extrusion system works by using piecewise linear segments to approximate a smooth curve. As shown in Figure 6, to synchronize the nozzle motion with the material extrusion, the volume-equilibrium between the deposited material and the feed material is used to determine the feed step l m . Specifically, this volume equilibrium, of which the left side is the volume of the deposited material from p i to p i+1 and the right one represents the cylinder volume of feed material from e i to e i+1 , is given as: where the nozzle diameter s, the layer thickness h, and the radius of feed material r m are user-determined and l f is the forward step computed by Equation (12). Accordingly, the feed step of the supplied filament l m is fully decided.

+1
where is the curvature radius of S along the direction of * at point and * is the specified maximal chord-error. For the feed-material flow, it is achieved by the mechanics of the motor controllers of the 3D printer. In our experiment setting, the motor control of the extrusion system works by using piecewise linear segments to approximate a smooth curve. As shown in Figure  6, to synchronize the nozzle motion with the material extrusion, the volume-equilibrium between the deposited material and the feed material is used to determine the feed step . Specifically, this volume equilibrium, of which the left side is the volume of the deposited material from to +1 and the right one represents the cylinder volume of feed material from to +1 , is given as: where the nozzle diameter , the layer thickness ℎ, and the radius of feed material are user-determined and is the forward step computed by Equation (12). Accordingly, the feed step of the supplied filament is fully decided.

Feed material
Deposited material +1 +1 h Figure 6. Illustration of material extrusion synchronization.
As for the nozzle orientation assignment, we enforce that the nozzle axis should align with the surface normal at the NL point, lest the local gouging between the printed layer and the nozzle may occur. In such a way, a single NL curve containing the nozzle location and the orientation can be well determined.
After the initial NL curve is decided, its adjacent path can be expanded accordingly. Akin to carefully controlling the sampling density of an individual NL curve, to uphold As for the nozzle orientation assignment, we enforce that the nozzle axis should align with the surface normal at the NL point, lest the local gouging between the printed layer and the nozzle may occur. In such a way, a single NL curve containing the nozzle location and the orientation can be well determined.
After the initial NL curve is decided, its adjacent path can be expanded accordingly. Akin to carefully controlling the sampling density of an individual NL curve, to uphold the printing quality, the density of NL curves should also be strictly controlled. Let w s denote the side-step distance between two arbitrary adjacent drive planes. As discussed in Section 2.1, each point on an NL curve C P (θ * ) has a corresponding w e f f value, even along the same printing direction. To determine the side-step w s for the current NL curve, the minimal one w * s among all the m sample points on the NL curve will be taken as the most conservative side-step, i.e., Starting from the first drive plane Q 1 and its corresponding NL curve, by applying Equation (14), we obtained the second drive plane Q 2 and its NL curve, then the third Q 3 and its NL curve, . . . ., until the last one Q k and its NL curve when the opposite extreme edge of S was reached. As the side-step w * s (θ * ) is adaptively decided, the material deposition becomes more even and regular, which is opposite to the conventional equaldistance side-step approach. Algorithm 2, given below, summarizes this procedure.

Results and Discussion
To validate the proposed multi-axis nozzle path planning methodology, the prescribed algorithms have been implemented by us using MATLAB. In addition to computer simulation, as shown in Figure 7, for carrying out physical printing experiments, we built a prototype of a robotic manipulator-based multi-axis FDM printing system which integrates a popular 6-DOF UR5 robotic manipulator with a nozzle extruder. The processing parameters used in the experiments are listed in the following Table 1: Nozzle diameter = 1 mm Joint velocity limit v i,max = π 6 rad/s Joint acceleration limit a i,max = π 2 rad/s 2 Joint jerk limit j i,max = π rad/s 3 Maximal chord-error for forward step δ * c = 0.5 mm Major-axis length of filament ellipse s = 1 mm Minor-axis length of filament ellipse h = 0.75 mm In total, there are four printing surfaces selected for testing the proposed algorithms, i.e., a saddle surface (Figure 8a), a bicycle seat surface (Figure 8b), and two freeform surfaces with concave and convex regions (Figure 8c,d). For these test surfaces, which are originally given in CAD NURBs format, they are first discretized into triangulated meshes with the help of the software HyperMesh. During the process of triangulation, the specified maximal chord-error δ * c is used as a control variable of point sampling, so to ensure the accuracy of the mesh models. Then, the triangulated meshes as shown in Figure 8 are used as the input to the proposed algorithms, and the vertices of the meshes are used to calculate the material deposition rate.  Nozzle diameter = 1 mm Joint velocity limit , = 6 ⁄ Joint acceleration limit , = 2 2 ⁄ Joint jerk limit , = π 3 ⁄ Maximal chord-error for forward step * = 0.5 mm Major-axis length of filament ellipse s = 1 mm Minor-axis length of filament ellipse h = 0.75 mm In total, there are four printing surfaces selected for testing the proposed algorithms, i.e., a saddle surface (Figure 8a), a bicycle seat surface (Figure 8b), and two freeform surfaces with concave and convex regions (Figure 8c,d). For these test surfaces, which are originally given in CAD NURBs format, they are first discretized into triangulated meshes with the help of the software HyperMesh. During the process of triangulation, the specified maximal chord-error * is used as a control variable of point sampling, so to ensure the accuracy of the mesh models. Then, the triangulated meshes as shown in Figure 8 are used as the input to the proposed algorithms, and the vertices of the meshes are used to calculate the material deposition rate.  After running Algorithm 1 using the parameters of Table 1, the ( ) among different feed angle for the four test models are computed and displayed in Figure 9, with both the optimized angle and its corresponding deposition rate annotated. For S1, as shown in Figure 9a, the indicator indicates several distinct local maxima, and the best of which oc- After running Algorithm 1 using the parameters of Table 1, the v e (θ) among different feed angle θ for the four test models are computed and displayed in Figure 9, with both the optimized angle and its corresponding deposition rate annotated. For S 1 , as shown in Figure 9a, the indicator indicates several distinct local maxima, and the best of which occurs at θ * = 158 • (0.878π rad) with a maximum v e (θ * ) = 0.599 mm 3 /s. For the second model, S 2 shown in Figure 9b, since the surface is symmetric about the x-axis, the indicator v e (θ) shows a symmetry about θ = 90 • (barring the inaccuracy due to the meshing and numerical errors), and a global maximum v e (39 • (0.217π rad)) is identified. For the two freeform surfaces S 3 and S 4 (Figure 9c,d), the optimized angles are 22 • and 158 • , respectively, while there is little difference between the deposition rates with various angles, especially for S 4 . This is because, for these two surfaces, for any direction, the intersection curve between the corresponding drive plane and the surface always has portions of large curvature, which makes it difficult for Algorithm 1 to find a distinguishingly better direction angle. With the optimal angle * found by Algorithm 1, Algorithm 2 (see Section 3) is then executed, which iteratively generates the NL paths for the given curved surfaces. Finally, these NL curves are connected in a zig-zag pattern to form the continuous printing paths, as shown in Figure 10. With the optimal angle θ * found by Algorithm 1, Algorithm 2 (see Section 3) is then executed, which iteratively generates the NL paths for the given curved surfaces. Finally, these NL curves are connected in a zig-zag pattern to form the continuous printing paths, as shown in Figure 10. Among the four surface models, two are selected further for physical printing tests, i.e., the saddle surface 1 and the bicycle seat surface 2 . For these two surfaces, the printing paths generated from Algorithm 2 are first converted into the G-codes of the homebuilt multi-axis FDM printing system and then executed. For comparison purposes, a benchmark printing path is also executed for each of the two surfaces. To facilitate the printing, a support-base is pre-fabricated beforehand for each of the two surfaces. Figure  11 displays some snapshots of the two printing processes, wherein the red surfaces are the partially printed layers while the support-bases are in black. Among the four surface models, two are selected further for physical printing tests, i.e., the saddle surface S 1 and the bicycle seat surface S 2 . For these two surfaces, the printing paths generated from Algorithm 2 are first converted into the G-codes of the homebuilt multi-axis FDM printing system and then executed. For comparison purposes, a benchmark printing path is also executed for each of the two surfaces. To facilitate the printing, a support-base is pre-fabricated beforehand for each of the two surfaces. Figure 11 displays some snapshots of the two printing processes, wherein the red surfaces are the partially printed layers while the support-bases are in black. Figure 12 displays the final printed surfaces S 1 and S 2 , respectively, using both ours and the benchmarking method. In Figure 12a,d which show the printed parts using the proposed method, the black lines mark the optimized feed directions of the printing paths, with the optimized angles tagged by θ * , and w s denoting the side step between the two involved adjacent NL curves. For a more detailed examination, the regions circled in green are magnified, as shown Figure 12b,e, respectively. As a whole, the re-solidified material is well deposited and evenly arranged onto the support-base, which owes much to the adaptive determination of both the side step and forward step. However, the zoomed pictures also indicate that over-extrusion or under-extrusion of the deposited material may not be fully eliminated, for which we offer the following explanations. First, as aforementioned in Section 3, the motor control of the extrusion system in our setup is based on the piecewise linear approximation of a smooth NL curve, which may induce certain inaccuracy. Secondly, the current adaptive side-step decision is based on a simplified ellipse model, which may incur inaccuracy in the actual amount of the deposited material.
Thirdly, the controller of the UR5 arm may automatically alter the scheduled feed rate during the printing process whenever it deems necessary, which also affects the deposition rate. Regarding the comparison, Figure 12c,f shows photos of the printed S 1 and S 2 using the benchmarking method which is still the iso-planar method except that the drive plane direction is arbitrary. Note that, except for the different directions of the drive planes, both the proposed method and the benchmarking method use the identical method to determine the side-step w s and the feed rate. Therefore, as the photos show, both the proposed method and the benchmarking method are similar in terms of printing quality.  Figure 12 displays the final printed surfaces 1 and 2 , respectively, using both ours and the benchmarking method. In Figure 12a,d which show the printed parts using the proposed method, the black lines mark the optimized feed directions of the printing paths, with the optimized angles tagged by * , and denoting the side step between the two involved adjacent NL curves. For a more detailed examination, the regions circled in green are magnified, as shown Figure 12b,e, respectively. As a whole, the re-solidified material is well deposited and evenly arranged onto the support-base, which owes much to the adaptive determination of both the side step and forward step. However, the zoomed pictures also indicate that over-extrusion or under-extrusion of the deposited material may not be fully eliminated, for which we offer the following explanations. First, as aforementioned in Section 3, the motor control of the extrusion system in our setup is based on the piecewise linear approximation of a smooth NL curve, which may induce certain inaccuracy. Secondly, the current adaptive side-step decision is based on a simplified ellipse model, which may incur inaccuracy in the actual amount of the deposited material. Thirdly, the controller of the UR5 arm may automatically alter the scheduled feed rate during the printing process whenever it deems necessary, which also affects the deposition rate. Regarding the comparison, Figure 12c,f shows photos of the printed S1 and S2 using the benchmarking method which is still the iso-planar method except that the drive plane direction is arbitrary. Note that, except for the different directions of the drive planes, both the proposed method and the benchmarking method use the identical method to determine the side-step and the feed rate. Therefore, as the photos show, both the proposed method and the benchmarking method are similar in terms of printing quality. The real printing time and the total printing length have been recorded by the controller of our homebuilt multi-axis printer and they are listed in Table 2. The results indicate that the proposed method has reduced the total printing time against the benchmarking method by 10.8% and 35.5%, respectively, for S 1 and S 2 . For these improvements, we offer several plausible explanations. The first and also the primary one is that, by considering the effective printing volume as a global optimization objective, the proposed method has effectively determined an optimized drive plane angle to minimize the total printing time. Secondly, as revealed by the total path length data shown in Table 2, the printing paths generated by the proposed method are shorter than that by the benchmark method, which is chiefly due to the fact that, as the proposed indicator involves the effective printing width which is related to the local curvature of the printing layer, a higher deposition rate is always sought locally, which in turn will reduce the total printing length (since the total printing volume is fixed). Finally, it is noted that the proposed indicator considers not only the local curvature of the printing layer but also the kinematic capacities of the robot arm. Altogether, these factors contribute to a significant improvement in printing efficiency by the proposed method. The real printing time and the total printing length have been recorded by the controller of our homebuilt multi-axis printer and they are listed in Table 2. The results indicate that the proposed method has reduced the total printing time against the benchmarking method by 10.8% and 35.5%, respectively, for 1 and 2 . For these improvements, we offer several plausible explanations. The first and also the primary one is that, by considering the effective printing volume as a global optimization objective, the proposed method has effectively determined an optimized drive plane angle to minimize the total printing time. Secondly, as revealed by the total path length data shown in Table 2, the printing paths generated by the proposed method are shorter than that by the benchmark method, which is chiefly due to the fact that, as the proposed indicator involves the effective printing width which is related to the local curvature of the printing layer, a higher deposition rate is always sought locally, which in turn will reduce the total printing length (since the total printing volume is fixed). Finally, it is noted that the proposed indicator considers not only the local curvature of the printing layer but also the kinematic capacities of the robot arm. Altogether, these factors contribute to a significant improvement in printing efficiency by the proposed method.

Conclusions
The primary objective of this paper is to develop a multi-axis printing path generation algorithm for an arbitrary curved freeform layer, such that the kinematic capacities of the 6-DOF robot arm (which enables the rotation of the nozzle) can be fully utilized while the required printing quality in specific tolerance is upheld. Our solution is to first establish a printing efficiency indicator that considers both the local geometry of the layer surface and the kinematic capacities of the robot arm. The classic iso-planar method is then adopted to generate a printing path; wherein, the normal of the parallel drive planes is exactly the one that globally maximizes the efficiency indicator and the side-step between adjacent drive planes is adaptively adjusted to maintain the required printing quality. We have carried out physical printing experiments of the proposed method and compared with some benchmark printing paths, and the preliminary results give a positive confirmation of the proposed method.
Future work will focus on three topics. First, while the iso-planar type of printing path is prevailing in practice, other types do exist, such as the parallel-contour type and its various variations. Thus, it is natural to inquire if similar ideas of this paper can also be applied or extended to some of those non-iso-planar printing path patterns. Secondly, at present, the nozzle orientation is mechanically assigned-it is simply aligned with the surface normal at the NL point. Since the change of the nozzle orientation has a huge impact on the joint motions of the robot arm, it is plausible to include the nozzle orientation as a variable to optimize, tending to all the constraints such as collision-avoidance. Finally, our proposed printing path generation method is restricted to a given printing layer, while the printing layers are generated independently of how the printing paths are generated on them. Aiming at further improving the printing efficiency, it will be interesting to see if these two processes can be coupled together.

Conflicts of Interest:
The authors declare no conflict of interest.

Nomenclature θ
Angle used to indicate the feed direction C P (θ) A single nozzle location curve, a function of θ f p Feed direction at point p n p Normal direction of the surface at point p ρ Angular parameter defined in the n p -k p plane h Minor-axis length of filament ellipse (equal to the layer thickness) s Major-axis length of filament ellipse (equal to the nozzle diameter) w e f f Effective printing width A g (p) Enclosed gap area A o (p) Overlapping area WCS Workpiece coordinate system NCS Nozzle coordinate system BCS Base coordinate system of the manipulator T W N Coordinate transformation from WCS to NCS (x w , y w , z w , i, j, k) A single nozzle location point defined by position (x w , y w , z w ) and its associated orientation (i, j, k) with respect to WCS (x n , y n , z n , α, β, γ) NCS frame described with respect to WCS in terms of translation vector (x n , y n , z n ) and the Euler angles represented rotation vector (α, β, γ) w PS (γ) Parameter of singularity, a function of γ Θ = (q 1 , q 2 , q 3 , q 4 , q 5 , q 6 ) Joint configurations of the robotic manipulator V = (v 1 , v 2 , . . . , v 6 ) Joint velocities of the robotic manipulator A = (a 1 , a 2 , . . . , a 6 Joint accelerations of the robotic manipulator J = (j 1 , j 2 , . . . , j 6 ) Joint jerks of the robotic manipulator v e (θ) Material deposition rate, a function of θ f (θ) Feed rate, a function of θ l f Forward step between adjacent NL points w s Side step between adjacent NL paths Q Drive plane S Mesh surface V Vertices set of the meshed surface