Minimizing form errors in additive manufacturing with part build orientation: An optimization method for continuous solution spaces

2.1 Related work

Part build orientation is a key factor in AM as it is easily manipulated and has a clear influence on final properties. The effect of part build orientation has been the subject of many research efforts but remains an open issue [7]. Methods for determining the optimal part build orientation can be largely divided into two groups: (i) those evaluating a set of candidate solutions with respect to an objective function and (ii) those mathematically describing the solution space to explore this continuous space for the optimal solution. A selection of the most relevant related methods is showcased in Table 1. The interested reader may refer to refs [7,8] for recent reviews on the orientation problem.

The first group is the most common and starts by identifying candidate orientations. This process is either based on a set of rules (e.g., flat surfaces of the convex hull) or by discretization of the solution space, i.e., certain intervals of rotation about one or more predefined axes. In the next step, either an exhaustive search is performed where all candidates are considered or a guided search is conducted, e.g., by a genetic algorithm (GA). These discrete methods have certain advantages, but they inherently fail to consider the entire (continuous) solution space and therefore risk missing good orientations. In addition, any attempt to refine the search space by including more candidate orientations inevitably increases the computation time.

The second group, on the other hand, grants access to the entire solution space, where the complexity of this solution space generally correlates with the complexity of the geometry. This category is not as well explored, perhaps due to the simplicity of discrete solution spaces or the ability of discrete approaches to handle discontinuous functions. Nevertheless, a continuous solution space may facilitate more nuanced objective functions and complex solution spaces. Continuous solution spaces can be explored by evolutionary algorithms (EAs), or gradient-based methods where knowledge of the topology of the solution space is exploited in an iterative search for the global optimum.

The approach described herein constitutes a hybrid of the two methods outlined earlier: The solution space is described with a differentiable function, which is used to identify critical points. The critical points will then constitute the candidate orientations in the final evaluation, which identifies the global optimum. The novelty of this approach lies in the combination of simplicity from generalizing surface types, and the flexibility brought forward by the general framework that can be populated with any objective function.

Cheng et al. [9] considered all flat surfaces as candidates for determining the optimal part build orientation and evaluated every candidate orientation with respect to the accuracy, build time, and stability. A similar approach was proposed for the minimization of cost by Alexander et al. [5] who included the cost of postprocessing and build time in the cost calculations. Zhang et al. [16] generated a set of candidate orientations from every shape feature and evaluated them for several attributes including surface roughness and support volume. Similar methods have also been proposed with facet clustering, where groups of facets are considered collectively based on how they are affected by part build orientation [20,22,21]. All these methods benefit from the generalization of how different surfaces are affected by the build direction. While various methods are employed in these studies in the search for an optimal solution, none of them consider a continuous solution space generated from the geometry.

Zhang and Li [13] proposed a discretization of the solution space (i.e., a unit sphere), and let each facet of the STL file promote the two orientations parallel to the facet normal, as well as the great circle corresponding to the perpendicular of the normal vector. The authors describe methods to make the selection using an exhaustive search for minimizing VE, as well as a GA for combined optimization of VE and part height [13,14]. The authors argue for the use of GA over exhaustive search when the number of discretized points becomes large due to time concerns. Discretization allows for controlling the resolution of the solution space and therefore also the computational cost. However, the approach remains oblivious to any effects other than those of the predefined points in the solution space.

Das et al. [15,17] proposed a method for minimizing form errors and support structures by formulating a minimization problem to be solved in MATLAB with gradient-based optimization. The method is based on the one-dimensional tolerance maps of Paul and Anand [23,24] for cylindricity and Arni and Gupta [25] for flatness – both theoretically derived from the staircase effect. A similar approach was proposed by Chowdhury et al. [19] who continued to compensate the geometry for any expected deviations still present after optimization. However, both of these proposed methods involve a risk of convergence to local optima due to the gradient-based approach. The convergence to local optima is avoided in the study by Budinoff and McMains [18] who performs an exhaustive search of one-degree increments to identify feasible regions from which the final orientation may be selected.

The existing body of literature describes various approaches to the orientation problem with various benefits and drawbacks. Exploiting higher level information about local topology facilitates the generation of candidate orientations. This information can also be used to construct continuous solution spaces for achievable tolerances and other objectives. The method presented herein benefits from the generalization of part geometry to reduce the number of parameters in optimization, while also accessing the entire solution space for mathematical analysis.

2.2 Mathematical foundations

There is no shortage of mathematical formulations of the orientation problem in the academic literature. The plethora of formulations arises from the subtle differences in AM technologies, which have varying parameters with different effects. The number of formulations is further amplified by the deviating scope and objectives of previous works. For instance, the VE is a common measure of accuracy in AM. However, the calculation of VE is based on fundamental assumptions regarding the surface profile, typically assuming right-angled steps. However, the real surface will be filleted as various effects will round off the corners and hence throw off the theoretical models [26,27].

Nevertheless, the effect of part build orientation on surface accuracy is indisputable in the current AM systems, although of less concern when the layers are thinner [28]. Therefore, the modeling of quality as a function of part build orientation is warranted, and mathematical descriptions of orientations are necessary.

According to Euler’s rotation theorem, any orientation of a rigid body in R 3 can be described as a sequence of three basic rotations ( α , β , γ ) , where no two subsequent rotations are performed about the same axis. The basic rotations are performed about the x -, y -, and z -axis individually, and a rotation of θ degrees about the respective axes can be calculated using the 3 × 3 rotation matrices R x , R y , and R z where

The rotation matrices of equations (2)–(4) will rotate any column vector about the respective axis by the angle θ . The direction of the rotation is determined by the right-hand rule, i.e., counterclockwise as seen from the positive end toward the origin. According to Euler, any orientation can be achieved by three successive rotations. However, because matrix multiplication is non-commutative, the sequence of rotations influences the final orientation of the body. Consequently, when three rotations are performed in succession following Euler’s rules, there are still 12 possible combinations divided into two distinct groups:

1. Proper Euler angles where the first and last rotations are performed about the same axis; and

2. Tait-Bryan angles where all three rotations are performed about unique axes.

In the current work, the notation for rotation sequence is simply the axes in the order of which the rotations are performed with a subscript further emphasizing the sequence. For instance, the most common rotation following a proper Euler angle sequence is Z 1 X 2 Z 3 , which means first rotation about the z -axis, second rotation about the (new) x -axis, and finally another rotation about the (now current) z -axis. Table 2 tabulates all 12 rotation sequences following this notation.

The range of rotations is limited to the unit circle as rotations of more than 2 π radians make little sense. Hence, the range of α and γ can be restricted to [ 0 , 2 π ) . Conversely, the range of β need not exceed π radians because larger rotations will only repeat previous spaces; hence, the range of β can be reduced to [ 0 , π ) . In consequence, the orientation problem in AM is bound to the finite space, where ( α , γ ) ∈ [ 0 , 2 π ) and β ∈ [ 0 , π ) .

3.1 Mathematical description of part build orientation

Consider a part P = { F 1 , F 2 , … , F n } . These features may be of different sizes and feature types, i.e., planes, cylinders, spheres, cones, and tori. The quality of each feature type can be modeled as a function of some parameters Q F = f ( x 1 , x 2 , … , x n ) . If a consistent measure of quality is applied (e.g., surface roughness Ra), the mean quality of the entire part Q part may be calculated as follows:

where n is the number of features, Q i and A i are the quality and area of the ith feature, respectively, and A part is the total surface area of the part.

While the area may be a suitable weight factor in many cases, it is also possible to introduce a separate weight factor that enables the prioritization of features. This weight factor can be applied on any level, e.g., certain feature types may be ignored, or individual features may attain a higher priority. Furthermore, the fundamental assumption is made that the quality of an additively manufactured surface can be described as a function of its orientation with respect to the build direction.

The orientation of a feature F with respect to the build direction z ˆ may be described as two Euler angles ( α , β ) necessary for the rotational transformation of z ˆ to v → where v → is the feature vector of F . In the current work, ( α , β ) represents rotations about the x - and y -axis respectively, both with reference to the original reference frame (i.e., Tait-Brian angles X 1 Y 2 Z 3 ).

The angle θ between two arbitrary vectors u → and v → is found by:

By utilizing unit vectors, ∣ u ∣ → ⋅ ∣ v ∣ → evaluates to 1 and reduces the expression to

By inserting z ˆ = [ 0 0 1 ] for u → , the aforementioned expression can be further reduced because

In this context, z will be the z -component of the feature vector v → , denoted v → z . However, we want to express the angle θ as a function of the part’s orientation in 3D space to enable optimization of orientation. For this purpose, the orientation vector can be expressed in terms of Tait-Bryan angles ( Z 1 Y 2 X 3 ). These rotations are commonly referred to as yaw, pitch, and roll in engineering applications and describes the orientation of a rigid body with respect to the world coordinate system (WCS). The feature vector can then be derived from the rotation matrices in equations (2)–(4), and the sequential rotations about the x -, y -, and z -axis can then be performed in a single operation by multiplying the matrices as follows:

This means that the orientation vector v → of a feature may be expressed as follows:

which means

When v → z from equation (11) is inserted in equation (7), the final expression simply becomes

To enable the mathematical description of the entire geometry in a single expression, the orientation of each feature is described relative to a common coordinate frame, i.e., the WCS. The geometry may be regarded as a rigid body, which means that the relation between the surfaces remains constant and any transformation acts on all surfaces equally. In consequence, the orientation of all surfaces may be collectively calculated from the same values of α , β , and γ using equation (10), where x , y , and z are the components of the feature’s initial orientation vector. Accordingly, equation (12) is valid for all surfaces with α and β as the only variables.

3.2 Finding the optimal part build orientation

Based on the aforementioned theory, it is possible to determine the optimal part build orientation mathematically by evaluating the critical points of the objective function. The critical points of a function f ( x , y ) are found, where

Typically, there will be several solutions to equation (13), and each solution needs to be evaluated separately. These solutions will represent points, edges, and perhaps even entire areas in the 2D solution space. Provided the formalization in the previous section, the part build orientation can be described as a set of two rotations α and β . If a function is based on these two rotations, α and β will replace x and y in equation (13) in the search for critical points. The solution space will be the surface of a unit sphere, where α ∈ [ 0 , 2 π ) and β ∈ [ 0 , π ) . A point on this surface will correspond to a single unique orientation, while a line will correspond to a range of orientations. Note that the entire unit sphere is accessible already with only half a rotation of β , and still, a full rotation of β is used in the visualizations of this article to make their analysis more intuitive.

The fundamental assumption remains that the quality of a surface can be modeled as a function of its orientation with respect to the build orientation and that this function is differentiable. A conditional function (such as the one in equation (1)) introduces certain challenges to this method. However, such discontinuities would represent edges and areas in the 2D solution space that could be added to the list of candidate orientations.

Each feature of the geometry will add a term to the objective function, and each term will typically add one or more candidate orientations to the list. The exact number of additional candidates depends on the mathematical model as higher order functions will yield more candidates. It is therefore beneficial to limit or minimize the number of terms to avoid excessive computations. A simple way to minimize the number of terms is to join similar terms, e.g., two features of the same type and orientation can be combined into a single term of the objective function. Other measures include a manual selection of significant features and automatic filtering of features based on type, size, etc.

As surfaces are affected differently by build direction, separate models for each feature type are necessary to enable proper evaluation. This is easily implemented by inserting the appropriate expression for F Quality into equation (5) for each feature.

4.1 Constructing the mathematical models

The following models are based on the orientation rules described by Frank and Fadel [32], stating that cylinders should be oriented with the axis parallel to the build direction, while planes can be oriented either parallel or perpendicular to the build direction. Up-facing and down-facing surfaces are treated equally in the examples to keep the objective functions simple. This can be naturally incorporated in the objective function to account for any additional effects, e.g., overcure, support structures, and so on.

The central assumption in the current work is that the accuracy of a feature can be modeled as a function of its angle to the build direction ( θ ), which in turn is a function of the part’s orientation as described in equation (12). Also, we are modeling deviations from nominal geometry, which can be considered a cost; hence, a low cost indicates a high fitness of a given orientation. However, the angle will not be sufficient to evaluate the fitness of a certain orientation. Consider, for instance, a cylinder oriented at a 4 5 ∘ angle from the build direction ( θ = 4 5 ∘ ). The same cylinder oriented at θ = 13 5 ∘ would yield the same result, but the angle θ is quite different. Clearly, the objective function should be more sophisticated to incorporate this behavior.

Because the angle θ always will be in the interval [ 0 ∘ , 18 0 ∘ ], the sine of the angle will provide three desirable properties of an objective function: (i) the result is always a number between 0 and 1, (ii) the function is minimized at vertical orientations and maximized at horizontal orientations, and (iii) the function is periodic and symmetric. This study employs this function as an expression for the quality of cylinders:

where θ is the angle between the feature normal vector and the build direction. This can be formulated as a function of ( α , β ) by inserting equation (12) for θ . With this substitution for θ , equation (14) may be written as follows:

Planes, on the other hand, require some additional configurations as we must consider both vertical and horizontal orientations as positive. To incorporate this new behavior, equation (14) is multiplied by cos 2 θ . This term ensures that planes are positively evaluated at both vertical and horizontal orientations. In addition, the function is more sensitive to minor changes when a plane is horizontal than when the plane is vertical. This fit well together with the VE being large when the plane is close to horizontal, while not being as prominent in close-to-vertical orientations.

Finally, the function is normalized to facilitate comparison with other feature types, and so on. The normalization is achieved by introducing a divisor equal to the maximum of the function. The maximum is easily obtained by derivation and reveals a normalization factor of 2.598 when rounded to three decimal points. This yields the following expression for the quality of planes:

where θ is the angle between the feature normal vector and the build direction. As with equation (14), the expression for planes may also be rewritten as a function of ( α , β ) by inserting equation (12) for θ . The result may be formulated as follows:

The solution spaces of equations (15) and (17) are illustrated in Figure 4, where a 3D graph and a contour plot are presented for each of the equations.

Figure 4

Solution spaces for planes and cylinders α , β ∈ [ 0 , 2 π ] from equations (15) and (17), respectively. (a) 3D-graph for cylinders. (b) Contour for cylinders. (c) 3D graph for planes. (d) Contour for planes.

4.2 Case 1: A simple geometry

The first case study is the simple geometry reconstructed from the study by Cheng et al. [9] and presented in Figure 3. This geometry provides a gentle introduction to the method by enabling a step-wise analysis of the geometry and the accuracy model. The first step of the method is to obtain a numeric description of the geometry in the appropriate format. Table 3 provides the positions and orientations of the features defined relative to the PCS as illustrated in Figure 3.

The data from Table 3 are inserted into equations (15) and (17) one row at a time as follows:

1. The feature type determines which equation to use (equation (15) for cylinders or equation (17) for planes)

2. The variables x , y , and z are substituted with E x , E y , and E z from Table 3

3. The expression is multiplied with the relative surface area of the feature (final column of Table 3)

Following the progression mentioned earlier, we start with the first row from Table 3 and perform the following steps: (1) The feature is identified as a planar type, and we, therefore, use equation (17). (2) The values of E x , E y , and E z are inserted for x , y and z , respectively. For the first feature, this means that 0 is inserted for x and z , while 1 is inserted for y . (3) Finally, the weight factor is introduced by multiplying by the relative surface area as found in the final column of Table 3, namely, 0.093 (9.3%). These steps are displayed in the calculations of equation (18).

When all the data from Table 3 are inserted into equations (15) and (17), all the equations may be collected in a single expression for the entire geometry as follows:

This yields the solution space illustrated in Figure 5. The solution space reflects the symmetry and regularity of the geometry as the orientations where feature vectors align with the build direction are clear.

Figure 5

Solution space for case 1 based on equation (19). (a) 3D-graph for case 1. (b) Contour for case 1.

In the next step, the partial derivatives are calculated and evaluated according to equation (13). Figure 6 shows the graphs of the partial derivatives where the dashed lines of Figure 6(a) and (b) correspond to the contour lines where the derivative evaluates to zero or are undefined.

Figure 6

Graphs displaying the partial derivatives of equation (19). (a) Contour for ∂ ∂ α . (b) Contour for ∂ ∂ β .

The critical points are found where both derivatives either evaluate to zero or are undefined. This can be displayed graphically by plotting the dashed lines of Figure 6(a) and (b) in a single figure as shown in Figure 7. Finally, individual evaluation of these points and edges must be conducted to identify the global optimum as tabulated in Table 4.

Figure 7

Contour lines for the partial derivatives of equation (19).

The evaluation of α , β ∈ [ 0 , 2 π ) reveals four solutions with equal cost. However, these four solutions are pairwise identical, i.e., ( 27 0 ∘ , 18 0 ∘ ) is the same as ( 9 0 ∘ , 0 ∘ ), and ( 9 0 ∘ , 18 0 ∘ ) is the same as ( 27 0 ∘ , 0 ∘ ). Moreover, these two unique orientations are polar opposites, corresponding to the object lying on its left or right side as exemplified in Figure 8(a). This evaluation of the optimal orientation is also consistent with previous assessments of the same geometry [9,11,16].

Figure 8

Optimal orientations identified for case 1. (a) Identified optimum ( 9 0 ∘ , 0 ∘ ). (b) Second-best orientation ( 0 ∘ , 9 0 ∘ ).

Another orientation achieving a low cost is the upright position which is achieved for any value of α when β is 9 0 ∘ or 27 0 ∘ . This constitutes two edges in the solution space with identical solutions. These orientations correspond to the front or the back facing upwards, which aligns the largest cylinder with the build direction as displayed in Figure 8(b). Note that the up-facing and down-facing surfaces are not differentiated by the objective function. Incorporating this behavior would yield different results where repetition of the solution space is avoided.

4.3 Case 2: Example with more complex geometry

The second case is an original design created for the sole purpose of demonstrating the applicability of the proposed method on a slightly more complex geometry. The part displayed in Figure 9 is a joint with connectors in various directions. Forty-four features may be identified where 15 are cylindrical, and the remaining 29 are planes. A numeric description of the geometry is presented in Table 5.

Figure 9

The joint designed for case 2.

Using the functions for planes and cylinders from equations (15) and (17) populated with the data of Table 5, the solution space of Figure 10 is obtained. Because many features share orientation vectors, the overall objective function can be simplified to contain a minimal number of terms.

Figure 10

Solution space for case 2 based on equation (20). (a) 3D graph for case 2. (b) Contour for case 2.

The solution space for the second case is also quite regular as demonstrated in Figure 11 where the partial derivatives of equation (20) are displayed. The regularity of these plots reflects the redundancy in the domain α , β ∈ [ 0 , 2 π ] .

Figure 11

Graphs displaying the partial derivatives of equation (20). (a) Contour for ∂ ∂ α . (b) Contour for ∂ ∂ β .

As demonstrated in case 1, the critical points are found where both partial derivatives are either zero or undefined, which corresponds to the intersections of red and black lines in Figure 12. All critical points are tabulated in Table 6 where the minima are highlighted in bold text.

Figure 12

Contour lines for the partial derivatives of equation (20).

For the second case study, the optimal orientation is achieved by 9 0 ∘ rotation about the x -axis as displayed in Figure 13(a). The same cost is also observed for three other combinations of rotations as displayed in Table 6. Due to redundancy in the solution space, these four combinations of rotations only correspond to two unique orientations in exactly opposite directions. With the hexagonal protrusions oriented parallel to the build direction, 25 out of the 29 planes are in a favorable orientation, while only two of the 15 cylinders are vertical.

Figure 13

Optimal orientations identified for case 2. (a) Identified optimum ( 9 0 ∘ , 0 ∘ ). (b) Alternative orientation ( 0 ∘ , 6 3 ∘ ).

Clearly, the large number of planes favors the orientation at ( 9 0 ∘ , 0 ∘ ). However, if cylinders are given a higher priority than planes, this orientation may no longer be as favorable. Figure 14 compares the effect of assigning a weight factor to cylindrical features for three orientations. A point of intersection is identified at a weight factor of 3.4, where the orientation ( 0 ∘ , 6 3 ∘ ) becomes more favorable than the previous best at ( 9 0 ∘ , 0 ∘ ) (see Figure 13(b)). This analysis indicates that for this part, cylinders should be at least 3.4 times as important relative to planes before the orientation at ( 0 ∘ , 6 3 ∘ ) is selected.

Figure 14

Comparison of three orientations for case part 2 when an additional weight factor is introduced for cylinders.

Both case studies exhibit repetitive symmetric patterns in the evaluation table. This symmetry arises from two sources: (i) as stated in Section 2.2, the domain of β can be constrained to [ 0 , π ) as rotations exceeding this range will only repeat previous solutions, and (ii) the objective function of these case studies takes no regard of the difference between up-facing and down-facing surfaces. One may also question the alignment of many critical points with right-angle orientations. This is, however, a result of the initial orientation of the part, which originates from the design phase. If the initial orientation was different, the patterns observed in the plots would change due to the projection, but the solution space would remain unaltered.

A more complex geometry with more features would yield a solution space with more local extremes where the complexity of the solution space reflects the complexity of the geometry. Naturally, the objective functions for each surface type also contribute toward the topology of the solution space. Consequently, different objective functions, e.g., those obtained through experiments, would give different results to those reported earlier.