A new overhang constraint for topology optimization of self-supporting structures in additive manufacturing

This work falls within the scope of computer-aided optimal design, and aims to integrate the topology optimization procedures and recent additive manufacturing technologies (AM). The elimination of scaffold supports at the topology optimization stage has been recognized and pursued by many authors recently. The present paper focuses on implementing a novel and specific overhang constraint that is introduced inside the topology optimization problem formulation along with the regular volume constraint. The proposed procedure joins the design and manufacturing processes into a integrated workflow where any component can directly be manufactured with no requirement of any sacrificial support material right after the topology optimization process. The overhang constraint presented in this work is defined by the maximum allowable inclination angle, where the inclination of any member is computed by the Smallest Univalue Segment Assimilating Nucleus (SUSAN), an edge detection algorithm developed in the field of image analysis and processing. Numerical results on some benchmark examples, along with the numerical performances of the proposed method, are introduced to demonstrate the capacities of the presented approach.


Introduction
A topology optimization problem consists in the formulation of a lay-out problem, where the goal is to find the optimal distribution of material in a specific region, according to the applied loads, the possible support conditions, the volume of the structure to be constructed and possibly some additional design restrictions. The first topology optimization method (the so called homogenization method) was proposed by Bendsøe and Kikuchi (1988). In this method porous microstructures are assumed to relax the optimization problem since the binary discrete setting for structural compliance designs is known to be ill-posed. However, the determination and evaluation of optimal microstructures and their orientations was cumbersome if not irresolvable, and the resulting structures could not be built as there was not associated length scale to these microscopic materials. Alternatively, density based methods were proposed, where the design variables are the elemental densities representing solid material or void. In order to approximate the solution to a discrete solution, material properties are scaled according to the density value of each element. One can find several approaches as the Solid Isotropic Material with Penalization approach (Bendsoe 1989), the Rational Approximation of material property approach (Stolpe and Svanberg 2001) and the SINH method (Bruns 2005) in the literature.
At present the SIMP method can be stated as the most popular topology optimization method, but apart from this approach there is an extensive amount of other alternative methods. Some are based on more heuristic concepts, such as the evolutionary methods, which have undergone great development over the past decades. Among all its variants, the Bi-directional Evolutionary Structural Optimization (BESO) method and the Sequential Element Rejection and Admission (SERA) approaches have been extensively investigated (Huang and Xie 2010;Querin et al. 2017). Some other methods explore shape and boundary variation and can be considered relatively new, like Topological Derivative methods (Sokolowski and Zochowski 1999) and Level Set methods (Allaire et al. 2004;Wang et al. 2003). Other recent approaches that recall shape and size optimization procedures for topology optimum design are the Moving Morphable Conponent (MMC) and the Moving Morphable Void (MMV) methods, introduced in Guo et al. (2014) and Norato et al. (2015), where the topology optimization problem can be solved based on explicit geometry description. These methods were further improved in  and (Gou et al. 2016). A comprehensive review of these methods can be found in the monograph by Bendsoe and Sigmund (2003), the general survey by Deaton and Grandhi (2014) and the comparative review by Sigmund and Maute (2013).
Even if topology optimization has shown to be a powerful technique that allows for increasingly efficient designs, because of the complexity and intricacy of the solutions obtained, it was often constrained to research and theoretical studies (Zegard and Paulino 2016). Recently, Additive Manufacturing is filling the gap between topology optimization theory and its application. AM is a rapidly evolving field that first emerged in the 1980s, and by now, because of the important technological developments, one can directly manufacture end-used parts, opening the design process to new challenges (materials, shape and internal structures) and enhancing the design freedom (Ponche et al. 2014). Hitherto, along with technological development, a large variety of AM processes have been proposed, for both metallic and plastic materials. However, and although AM overcomes the limits currently imposed by the conventional manufacturing technologies, the processes based on powder beds and polymer materials still present construction limitations and there are some technical constraints that must be satisfied in order to generate consistent geometries (Doubrovski et al. 2011). These constraints include the minimum member size, the overhanging distance and the member inclination (Chang and Tang 2001;Leary et al. 2013). Extremely thin members cannot be properly printed if that thickness is below the minimum allowable thickness of the tool. The fabrication overhang angle is another example of a rule which is of great importance to ensure that the part will not collapse during the layer adding process when additive manufacturing is used. Generally the inclination of a member is described as the angle between the base plate and the down facing contour. This is similar to considering the angle between the material growing direction and the vector normal to the contour. In this work, the member inclination referred to the material growing direction will define the overhang and whether it is properly supported or not. The theoretical threshold for the overhanging angle from which the member is considered non self supported is commonly considered to be 45 degrees (Daniel 2009), though there are some authors that propose that not only the overhanging angle but the combination of overhanging angle and overhanging length should be considered (Vayre et al. 2012). Other authors propose different values for that angle, ranging from 20 to 45 degrees (Leary et al. 2014). A structure satisfying such an overhang angle constraint is called self-supporting. For a non self-supporting structure to be manufactured, its geometry has to be modified or additional support structures need to be generated. Such support structures are formed during the fabricating process so that the primary object can be manufactured layer by layer without collapsing. This is very time-consuming and a waste of material, rising the issue of further post-processing activities to remove the unwanted supports.
Modifying the geometry manually in order to fulfill the overhang constraint may ultimately reduce the structure's physical performance and lead to structures with a compliance that is far from the optimum support-needed structure. Since the idea of coupling topology optimization and AM was formulated, many different considerations have been proposed but generally we can group them in three main groups, depending on the level of engagement between the two technologies. In the first group we find the standard topology optimization process, where once the structure is optimized, the supporting material is added where it is required during the manufacturing process. A second level engagement (Leary et al. 2014;Hu et al. 2015) gathers strategies that consist in adapting the optimized design to AM processes. Once the topology optimization is over, the generated geometry may or not be suitable for AM. In case the use of scaffold structures is not an option, the designer proceeds to introduce some manual changes on the geometry adapting it to the fabrication process. These engagement strategies include a critical step where designer's experience plays an important role. This step is where the designer introduces manually the changes on the post topology optimization geometry and creates some variations of that geometry in order to avoid bad supported overhanging members. Hence, there is a level of interference and uncertainty with the new level of optimality and deterioration of the structures' physical performance. Finally, the third level engagement considers the strategies involving a total integration of topology optimization and AM processes by introducing the overhang constraint into the design process, which enables a final result that can be directly manufactured and shows and acceptable compliance ratio with respect to the non-supported reference optimum structure. The main advantage of this approach is a total integration of topology optimization and design process with the AM procedure.
This paper progresses in this line of work and presents a new topology optimization approach that helps to design optimal self-supporting structures, which are ready to be fabricated via additive manufacturing without the usage of manual modifications or additional support structures, so that their associated optimal compliance is close to that of the reference structure. A strategy to solve the issue of designing a completely self-supporting structure via topology optimization was presented by Brackett et al. (2011), where the overhang angle was included in the optimization problem but did not produce a complete self-supporting structure. Gaynor and Guest (2016) introduced a wedge-shaped filter during the optimization process to obtain self-supporting structures. Recently, Langelaar presented a novel self-supporting filter into the topology optimization process, which is achieved via building smooth approximations to the minimum and maximum functions (Langelaar 2016a, b). Also very recently, Qian (2017) has proposed a density gradient based approach for undercut and overhang angle control in topology optimization. Recently, Guo et al. (2017) introduced a novel formulation for AM-oriented topology optimization where the corresponding problem can be solved in a geometrically explicit way and simultaneous optimization of the selfsupporting structural topology and the angle of working plane is achieved for the first time. Theoretical issues associated with AM-oriented topology optimization are also discussed in this work, which are important to examine the effectiveness of different approaches.
The approach proposed in this paper follows a different way and develops a novel global inequality constraint for the overhanging angle which is included explicitly in the topology optimization formulation problem. This constraint can be easily added to the problem formulation as an additional inequality constraint along with the regular volume constraint for compliance minimization problems, and coupled with traditional mathematical programming optimization algorithms, like the Method of Moving Asymptotes (Svanberg 1987). The constraint that is proposed in this paper refers to the ratio between acceptable and not acceptable material members, and is based on the evaluation of self-supported "contours", instead of the traditional formulation of supported and unsupported "elements". The methodology used to identify and control the inclination of members is based on an edge detection algorithm developed by Smith and Brady (1997) in the field of image processing, known as the Smallest Univalue Segment Assimilating Nucleus (SUSAN). The implemented procedure allows the designer to specify not only any critical overhang angle and printing direction, but also the tightness of the constraint through a harshness parameter that controls the ratio of selfsupported contours. The capacities of this approach are demonstrated with extensive numerical examples and obtained results are compared with previous studies.
Finally, making a comparison of the proposed method with other methods existing in the literature, the following aspects could be highlighted. Most of previous works need the incorporation of filters that must be included in the traditional density-based topology optimization procedures (Serphos 2014), or they introduce some sort of additional constraints to control several problems like lateral overhangs and the grayness level of the density field (Qian 2017). On the contrary, the method proposed here needs only a single overhang constraint, defined as a relative value instead of an absolute value, which allows controlling the impact of the restriction easily and makes it straightforward to choose an adequate value for it. Concerning the SUSAN algorithm adopted for edge detection, it must be said that the choice of this algorithm was not casual, since it has been demonstrated that it is extremely efficient for contour detection. The inclination of the members is detected precisely and locally along the full contour of the structure, while other methods may require introducing a series of auxiliary constraints to efficiently detect not self-supported edges and singular parts (Guo et al. 2014). Its extension to 3D problems is also straightforward, as it can be checked in the work by Walter et al. (2009aWalter et al. ( , 2009b.

Edge detection and overhang constraint
In this section the contour detection algorithm and the selfsupporting constraint are formulated, which will be integrated within the modified SIMP framework in terms of the element density, also known as RAMP parameterization. First we shall develop a tool based on edge detection capable of finding forms and shapes in a grey-scale topology optimization model that will be used to evaluate the amount of feasible and nonfeasible contours during the topology optimization iterations. Following on from this, we will present the proposed explicit overhang constraint and use it for analysis and sensitivity analysis of the manufacturability of the structure. A brief description of the well known Heaviside projection technique for density filtering is also included at the end of the section, since it will be used later for the sensitivity analysis.

Contour evaluation algorithm
In order to obtain a magnitude that describes the overhang situation of the structure, we will start from an edge detection algorithm that is typically used in image processing, but here the pixel intensity is substituted by the element material density. The adopted algorithm is called Smallest Univalue Segment Assimilating Nucleus also known as SUSAN. The performance of this algorithm is based on sweeping images (domains) with circular masks along them and analyzing the image intensity gradient created by the pixels covered by the mask. It counts the value that characterizes the similarities between each pixel of the image and its neighborhood (Mokrzycki and Samko 2009). This edge detection algorithm is adopted here since it has already shown its reliability and effectiveness in the field of image processing. The information about edge location that this method reports is as close as possible to the correct position, the number of false negatives and false positives is minimum and it is fast enough to be usable in an optimization system that uses iterative evaluation of functions (Fynbo et al. 2001). The masks are located so that the element in the center of the mask matches the element being analyzed. In order to adapt the algorithm to the design space the mask used will be squared with a size of 3 × 3 (see Fig. 1). Higher size masks could be chosen but it is not recommended since they generate a blurring of the spatial information which may produce inaccurate results. The contour detection algorithm works calculating the density gradient to detect geometric boundaries. The inclination of this gradient is of great importance as the overhanging angle is stated as the deviation between this gradient and structure's growing direction. The density gradient offers some unique properties that are helpful and that are used to facilitate the process: the density gradient joints the geometrical center and the gravity center of the mask pointing always to the gravity center, the gradient vector is always normal to the boundaries, and the module of the vector defines the existence of a possible boundary. Quick variations of the density are presented by high values of the density gradient module and represent the closeness to a border.
To obtain the density gradient the coordinates of vector v cg (x cg , y cg ) which represents the gravity center of the actual 3 × 3 mask are calculated first in the xy reference system. This vector is calculated using (1) and (2) from mechanics in the following way: where x i and y i are the local coordinates of the center of the elements inside the mask and ρ i represents the density of each element inside the mask. Coordinates x i and y i are calculated with respect to a local xy coordinate system whose origin matches the geometric center of the mask. Here we have considered a regular quadratic grey map so the surface area A i of the finite elements in the mesh can be left out of the equation. The information supplied by this vector is a direction and a weight, i.e., the length of the vector for each element, where the weight is a measure of the element's quality as an edge. Masks in areas of uniform density will have centers of gravity close to the geometrical centre, and the weight will therefore be small. As it will be shown in the next chapter, comparing the direction given by the gradient with the growing direction of additive printing, it is possible to compute both the deviation of the contour inclination from the limit value of the overhanging angle and to quantify a mean value of the allowable and non-allowable edges forming the topology of the structure.
Finally and concerning structure edge detection, it must be noted that an expanded dummy finite element mesh is introduced in the background, so that we can take into account the outer boundaries of the rectangular design domain (Garaigordobil et al. 2017). This dummy mesh consists of the regular mesh that defines the density distribution and two additional rows and columns of elements bordering the boundary of the design domain. The density of these dummy elements is taken to be zero as they always correspond to the outside of the design region, except at the lower boundary where the plate of the printing machine is located and density is equal to one. Obviously, this auxiliary mesh is for contour detection purposes only and it is not considered during the finite element analysis of the structure. Figure 2 shows how effective the contour extraction based on this technique is when applied to the geometry obtained from the solution of a topology optimization problem. Light areas correspond to elements that are close to the border of the structure, while dark elements represent material or void areas which are located far from the edges of the structure.

Overhang constraint
The main objective of this work is to develop an efficient and versatile constraint for controlling the overhanging angle or member inclination. Generally the overhanging angle is presented with respect to the platform or base-plate, however, in the following pages, seeking for a better understanding of the calculation, this angle will be calculated respect to the growing direction. Here, the overhang angle α is the relative angle between the growing direction and the edge gradient vector v cg . Since the gradient vector is always normal to the edge of the member, the angle α will represent also the physical inclination of any boundary referred to the base-plate. Hence, any member is considered self supported if the overhang angle α overcomes a predefined threshold angle ψ (see Fig. 3).
In order to build the proposed overhang constraint, the density gradient v cg of each mask will be compared with the growing direction, a direction that is defined by the designer with the unitary vector r that is normal to the plate. The deviation between these two vectors is the angle previously named as α. This angle is compared with the limit value of the overhanging angle defined by the designer, which stands as ψ. The designer is free to choose the value of ψ, that is, the critical value for the member inclination α, among the different values proposed in the bibliography which depends on the material and the manufacturing process used. Hence, the application of the procedure described here for any value of the overhang angle is straightforward, as it will be shown later. These two angles may be compared to calculate the deviation between the gradient vector v cg and the allowable vector v ψ for each mask using the function named φ(ρ), which depends on the density distribution vector ρof the mask and should be non positive: As it can be seen in Fig. 3, when the value of function φ is larger than zero, the contour defined by the material densities associated to the central element in that mask is said to be not supported. However, in order to avoid potential indeterminations when |v cg | = 0, (3) will be reformulated here reconsidering a different equation. Rather than comparing angles, we shall compare the vertical projection of the density gradient (y cg ) and the vertical projection bound (y ψ ) of vector v ψ , building a new function to evaluate the deviation between admissible and inadmissible gradient vectors for the contours associated to each mask: Effectively, we can use any of these equations for our purposes, as it will be shown immediately. After substituting the components of vectors r and v cg and operating we can get: Since (5) is compared to 0 so that bad supported elements are detected, we can simplify the expression considering only the numerator in the following way: For a growing direction perpendicular to a horizontal base plate, that is, r = (0,1) T , we get: Rearranging terms and calculating the square of right and left hand side terms, after some trigonometric substitutions the following equation is obtained, where the negative root is not considered because it would result in an imaginary solution.
A negative value of the y cg coordinate always corresponds to an admissible edge, therefore we will always consider the positive value of the right hand side of (7): It can be seen in Fig. 3 that the right hand side of (9) represents the maximum allowable vertical projection of the gradient vector, previously denoted by y ψ . Therefore (3) and (4) are equivalent and both of them can be used to decide if a contour is admissible or not. This expression will be of great importance later for the definition of the global constraint function proposed in this work, so it will be normalized in the following way: It can be seen in (10) that now there is not any issue with this alternative expression because the term in the denominator that could became null has been eliminated. Once again, when φ is positive the central element in that mask is said to be not supported because the vertical projection of the contour's gradient is larger than the allowed projected vertical length.
Moreover, it can be demonstrated easily that (10) corresponds to the cross product of vectors v cg and the unitary vector v ψ , considering both vectors located in the first quadrant, since we want the same behavior for vector that are symmetric with respect to the vertical growing direction. The (10) will be used to classify the contours associated to each of the m masks in the design domain. If the angle α overcomes the angle ψ, then that mask is said to be properly supported and it will be verified that φ m ≤ 0. On the other hand, if α is smaller than ψ a non appropriate contour is detected and φ m > 0, so it is possible to divide the masks into two groups, φ − and φ + , where the absolute value of the sum in each of the groups represents an amount of appropriate and non appropriate contours in the structure, respectively, that is: where M corresponds to the total number of masks in the design domain and η m represents a weighting factor defined as the sum of the elements' densities in the mask (13). This factor is a constant for each mask and for each iteration and when multiplied with φ m reflects the quality as an edge of the grey scaled contour associated to the corresponding mask, since elements with high densities produce better defined shapes than elements with low ratios of material.
Once the previous classification is done, we can introduce the new global constraint that is proposed in this work. This constraint is defined as the ratio between the value of self supported contours (φ − ) and the total amount of admissible and inadmissible contours (φ − + φ + ), which should be generally as close as possible to one in order to obtain a solution free of non supported contours. The value of this ratio will be of great help to control the support-free manufacturability of the final optimum solution, so it will be considered as a tunable harshening parameter named Φ 0 depicting the "overhang rate" of the structure. Changes of this parameter have a direct effect on the final optimum solution, as it is described in the section that covers the numerical application examples at the end of the paper. (14) represents this additional global constraint that will be added to the general topology optimization formulation problem.

Density filtering and heaviside projection
A common problem in topology optimization of continuum structures is the occurrence of checkerboard patterns, or regions of alternating solid and void elements, in the final solution. This phenomenon of alternating presence of solid and void elements ordered in a checkerboard like fashion makes the interpretation of optimal material distribution and geometric extraction for manufacturing difficult. Diaz and Sigmund (1995) and Jog and Haber (1996) showed that checkerboards are not optimal but rather result from numerical instabilities. A large amount of literature exists on preventing checkerboards patterns and mesh dependence. Popular approaches are to restrict the design space placing a constraint on the total perimeter of the structure so that solution exists for the original continuum problem (Haber et al. 1996) or to filter the values of sensitivities (Sigmund and Peterson 1998), which have been used extensively in a significant amount of works. In this work we will adopt a Heaviside projection method (Bourdin 2001), since density filters have demonstrated to be more robust than sensitivity filtering schemes. First, an intermediate field of filtered design variablesρ e is obtained using a weighted average function in the following way: where ω i is the weight factor, S e is the set of i elements in the domain of influence of element e, that is, the set where the center-to-center distance to element e is smaller than the filter radius r min : The smoothing operation above limits the space of possible designs and solves the mesh dependency problem. However, the smoothness of the designs also implies gray transition zones at the interface between material and void phases (Guest et al. 2004). In order to avoid these gray zones, a projection step is added to the density filter where intermediate densities are projected by means of a regularized Heaviside function. In this work a continuous approximation of the Heaviside function based on the hyperbolic tangent function is used (Wang et al. 2011), where the values ofρ e below T are projected to 0 and the values above are projected to 1: where β is a scaling parameter which controls the steepness of the continuous approximation of the Heaviside function and T is the threshold parameter of the Heaviside function. The projected densities ρ e are referred to as the physical densities and will substitute the regular density design variables; therefore we will always present the projected density field ρ e rather than the original density field ρ e as the solution to the optimization problem. Similarly, the objective and constraint functions will be computed usign physical densities, provided that the variable ρ e is replaced with ρ e . It must be noted that the sensitivities in the topology optimization problem become increasingly ill-conditioned for high values of the parameter β and a continuation scheme on the parameter β is usually necessary during the optimization process (Andreassen et al. 2011). When working with large values of β using the default MMA approach as we do in this investigation, the sharpness of that approximation to the step function may generate aggressive oscillations. Therefore, there are some adjustments that should be considered when a continuation scheme is applied, since there is an iteration from which its value is great enough to generate these oscillations. The required parameter adjustments that should be included in the code when the MMA algorithm is adopted can be found in the paper by Guest et al. (2011), which basically consist in modifying the distance of the asymptotes from the current point in terms of the parameter β. Similarly, we will also make use of the continuation scheme proposed for the parameter p in the same paper.

Problem formulation and sensitivity analysis
This paper considers the maximum stiffness topology design where a new global overhang constraint is included in order to control and minimize the need of scaffold structures in additive manufacturing of optimized designs. This problem will be formulated in the following way: s:t: where u is the nodal vector of displacements, K is the global stiffness matrix that depends on element's E e Young modulus, F represents the nodal force vector, and k 0 . is the element stiffness matrix for a solid element. V(ρ) and V 0 correspond to the volume fraction of the actual iteration and the maximum allowable volume fraction, respectively. Finally,Φ ρ ð Þ represents the overhang constraint and corresponds to the ratio between the values of allowable and total number of contours, which should be larger than a specified value Φ 0 .
The RAMP method that converts intermediate densities less efficient in the optimization will be used. Each element is assigned a density and its Young's modulus is still given by the well known density-stiffness interpolation scheme, but replacing the regular density variable ρ e with the physical density ρ e .
where p is the penalization parameter, E 0 is the Young's modulus of the solid isotropic phase and E min represents the modulus of the void material, which is taken to be 10 −9 in this work. In the following, the approach for the sensitivity analysis of the functions depicted in the topology optimization problem formulation is described. In order to compute the derivatives of any function f(ρ) involved in the problem with respect to the density of each element, we will apply the chain rule, so that the filtering and projection operations are considered. The sensitivities of the compliance and the volume and overhang constraints will therefore be calculated using the following rule: The second and third terms in the right side of the equation are obtained for any f(ρ) function as The derivative of the first term on the right hand side of (24) will vary depending on the nature of function f(ρ) considered in each case. It is well known that the derivative of the compliance with respect to variable ρ e can be easily obtained in the following way: The derivatives of the volume constraint with respect to the projected densities for substitution in the chain rule of (24) can be obtained in the following way: where v e represents the volume of finite element e. Finally, let us consider the sensitivity computation of the overhang constraint function with respect to ρ e . First we should normalize the inequality constraint shown in (14) in the following way: where a different name has been adopted to avoid confusion. Taking derivatives in (29) and omitting density dependence for simplicity we have: The sensitivities of the allowable and non-allowable contours, φ − and φ + , can be calculated taking derivatives of expressions (11) and (12): Finally, it is necessary to compute the derivatives of the contour values of each mask φ m for substitution in (31) and (32). Recalling the function built in (10) for contour evaluation, the derivatives can be obtained with the following expression: where the derivatives of the coordinates of the gravity center can be easily obtained with the relations given in (1) and (2) and remembering that we are using the physical density variables for contour evaluation once the smoothing and projection filters were applied: Once the sensitivities of objective and constraint functions are obtained, the topology optimization problem formulated in (18)-(22) is solved using the Method of Moving Asymptotes (MMA) developed by Svanberg (1987). MMA minimizes a set of sequential convex approximations of the original functions enabling the solution of problems with more than one constraint and is well known to be very efficient in the field of structural topology optimization. Figure 4 shows the flowchart of the proposed procedure that summarizes the different loops considered and the sequence of their application. Once a suitable reference domain is chosen, the basic data for structural analysis should be defined, including loads, boundary conditions, fixed void and solid regions, finite element mesh, etc. Then the topology optimization problem is defined for computation of the optimal material distribution over the reference domain, starting with an initial homogeneous distribution of material. Volume and overhang constraint bounds should be specified, as well as the normalized growing direction and the critical allowable inclination of members for additive manufacturing. After variables are smoothed and projected in order to obtain the physical density field, the contour evaluation loop starts, where the gradients for all the masks in the design domain are calculated and the overhang constraint value is obtained. For this distribution, the volume and compliance are computed by the finite element method. Finally, sensitivity analysis is performed and densities are updated with the help of the method of moving asymptotes. This loop is repeated until only a marginal improvement of the compliance is achieved over the last design. Once the optimization problem is over and in order to send the component to the 3D printing machines, the optimal material distribution should be interpreted so that the designer can build a CAD representation of the shape. Depending on the harshness parameter selected for the overhang constraint, there may be a little amount of elements that are badly supported and violate the overhang constraint, which could be fixed in this stage. However, this is not always a required step.

Numerical examples
To illustrate the functionality of the developed optimization procedure and the characteristics of the overhang constraint parameters, several examples are discussed in this chapter. In all cases a Young's modulus E 0 = 1 and a lower bound E min = 10 −9 are used in the RAMP law. As it was mentioned in the paper, the common practice of a continuation method will be applied on the RAMP exponent p and the Heaviside parameter β. The penalization parameter is initially set to 10 and increased by 2 every continuation step until reaching a maximum penalization of 18 while the Heaviside parameter is initially set to 5 and increased by 5 until a maximum magnitude of 25 (Gaynor and Guest 2016;Guest et al. 2011). Finally, the filter radius in the density filter is chosen equal to 4 times the size of the element and the threshold parameter T of the Heaviside function is set to 0.5 unless expressly indicated. The values adopted for these parameters have shown to be effective and lead to a very accurate geometry definition. It should be noted that the value of the filter radius does not produce particular effects on the proposed strategy and its influence is similar to traditional topology optimization. Although it is actually very important in other strategies in the literature where large filter values are required to regularize oscillations, no remarkable effect is present in this case. The printing direction is always considered vertical and upwards. It must be noted that although scaffold structures can be avoided in any 2D example if the structure layout is placed in the horizontal plane of the printing machine, they are still valid to demonstrate the proposed approach's ability in handling overhang constraints if required.

Two bar structure example
First this benchmark example is introduced as an academic application to understand the effectiveness of the previously described procedure for both boundary detection and optimization of the overhanging angle. The design domain and boundary conditions for the topology optimization problem are shown in Fig. 5. The design domain is a rectangular area with a width W and height H=W/3. A concentrated unitary load F is applied at the middle point of the top edge while the bottom edge of the design domain is clamped. The design domain is discretized with 120 × 40 equally-sized square finite elements. The maximum volume fraction for the minimum compliance problem is chosen equal to 0.2, and the ratio of admissible contours in the overhang constraint is set to 0.97. In order to check the ability of the proposed algorithm to find optimum topologies without members that exceed the   allowable inclinations for subsequent additive manufacturing, four different angles are considered in the overhang constraint: 45°, 60°, 80°and 90°. The analytic solution for this topology optimization problem is already known and corresponds to a truss-like structure with two bars forming 45°with the base of the reference domain. This means that if the allowable angle for member inclination is equal or less than 45°, we always obtain the result shown in Fig. 6a, because the optimum shape always fulfils the overhang constraint. When higher angles are introduced for the overhang constraint, this shape becomes inadmissible for additive manufacturing and the algorithm redirects the optimization process to find again the stiffest two bar truss-like structure with admissible inclination of members. As it can be seen in Fig. 6b, c and d, the slope of the inner walls grows gradually in order to fulfill the admissible angle specified in the overhang constraint, while the slope of the outer contour is free to adopt any shape since it does not have any influence in the formation of scaffold structures. In this sense, the bound of 90°is particularly interesting because in this case there is not allowed any hole inside the structure, since it would necessarily lead to a non-allowable contour. Actually, for a limit angle of 90°every density gradient should be horizontal or be pointing downwards, hence, y cg must be equal o lower than zero.

Cantilever beam
The next example solves the problem of optimizing the material distribution with an overhang restriction in the cantilever beam shown in Fig. 7. The design domain has a height H and length L = 2H. The left edge of the domain is clamped and a unit point load is applied at the center of the right edge. The maximum volume fraction allowed is 50% of the volume of the full design domain and it is discretized using 160 × 80 unit sized finite elements. The ratio of admissible contours in the overhang constraint is 0.97. In Fig. 8 the overhang boundaries that would eventually need some sacrificial support material are highlighted with dotted lines, assuming that the minimum allowable inclination of members so that they can be manufactured by AM processes is 45°. This exercise is very appropriate to demonstrate the capacities of the algorithm to deal with overhang constraints due to its complexity. It contains several tricky boundaries with remarkable lengths some of which are virtually horizontal. Figure 9 show the optimum topologies for the optimization problem with 45°, 60°and 90°o verhang constraints, respectively. It is also used to evaluate the impact of the objective volume fraction and the grayness level of the results.
It is worth to mention that the more severe the constraint is, that is, the larger the allowable angle for member inclination is, the greater the number of members and holes the optimum structure shows. Figure 9 shows that these members branch from the base to the dome of the structure, satisfying the prescribed overhang constraint while they carry the load and give support to the elements above. When the definitely restrictive limit angle of 90°is reached, where the presence of any hole involves a violation of the constraint, the implemented algorithm works coherently and yields a continuum plate filled with material. It is also noticeable the way the holes are shaped. Looking at the up-facing contours of the holes, we can appreciate that they form curves and are rarely straight segments, however, the down-facing contours, which are those subjected to the overhang constraint, posses a more straight segment nature. This behavior complies with the overhang limitation, as a curved down-facing contour could have somewhere along the member a tangent with an inclination that is inadequate and violates the constraint. For that situation to be avoided, these contours are virtually   Concerning the grayness of the optimum design, it can be seen that the obtained solutions are virtually black-andwhite except for a very low fraction of elements. The grayness level (GL) can be computed with (36) as suggested by Sigmund (2007), where n is the total number of elements, and the values obtained for the optimized topologies shown in Fig. 9 are 0.036, 0.058 and 0.012, respectively. These values are very close to zero and confirm that the topologies are almost 0-1 solutions.
In order to analyze the effect of the volume fraction on the final solution, the cantilever example was solved again with two different volume ratios, representing relatively low and a high values. In both cases a critical angle of 45°was considered. When a low target volume is specified, frequently great voids are created in the non-constrained problem. As big holes will generally result in not self supported areas and the available material to introduce supporting members is low, solutions may differ markedly from topologies obtained for larger volume fractions. Figure 10a show the optimum topology when the volume fraction is taken to be 0.35, where it can be seen that material is more localized at the bottom of the design domain, while the shape of the holes recall the water-drop geometry used frequently when designing for additive manufacturing ). On the contrary, for larger volumes (0,65) the algorithm is free to place the material on the top and generate a straight member that provides the structure with greater stiffness (see Fig. 10b).
A key concept in AM processes is the positioning of the structure to be manufactured in the printing machine, since depending on that positioning, the amount of non-supported members in the structure changes, as well as the mechanical properties that the final structure will present (Hofland et al. 2017). In case of the cantilever beam shown in this example, a vertical design domain may be more appropriate than the horizontal layout as the non-penalized optimum topology naturally shows members with inclinations close to 45°, suitable to be manufactured directly without the need of support material (see Fig. 11a). Moreover, the optimal topology that is obtained with the vertical setup preserves the symmetry axis of the nonconstrained regular solution. In the following it will be analyzed the geometry changes that the optimum solution suffers when different constraint angles are considered with the vertical positioning of the design domain.
Comparing the geometries generated for 45°and 60°overhang constraints, there can be noticed some small but meaningful variations in the geometry, where members sprout and die in different areas depending on the constraint angle. Effectively, for the 60°constraint some key points are moved downwards and more material is placed around the central members so that larger slopes can be introduced in the down-facing contours (see Fig. 11b). When the minimum allowable inclination is set to 80°, the topology of the optimum solution changes dramatically and almost vertical walls are built in order to fulfill the overhang constraint, as Fig. 11c shows. As in the previous horizontal layout case, a restriction of 90°implies again that no holes can be introduced in the domain because an inadmissible gradient would be created. Therefore a solid domain is created, as it is shown Fig. 11d. It is also noticeable that this solution shows a curved outer shape that matches approximately the shape of a solid cantilever beam with maximum stiffness.

MBB beam
Let us consider now the classical MBB topology optimization problem in Fig. 12. The goal of this example is to show the effect of the parameter Φ 0 , which represents the minimum ratio of self supported contours in the structure, as it was mentioned in the theoretical part of the paper. The rectangular domain is supported in the lower right and left corners and the load is applied vertically at the middle point of the upper edge. Taking advantage of the symmetry, only half of the domain is analyzed, applying a symmetric boundary condition along the left edge of the structure. The finite element mesh is uniform with 156 × 52 elements of unit dimension and the volume Fig. 9 Optimum topologies for different critical angles a 45°b 60°a nd c 90°A new overhang constraint for topology optimization of self-supporting structures in additive manufacturing restriction is taken to be 50% of the domain. In this example it is specified an overhang constraint of 45°and three different values of the parameter Φ 0 will be considered: 0.87, 0.97 and 0.99. Figure 13a corresponds to the optimized MMB topology when no overhang constraint is considered. The value of the parameter Φ 0 for this solution is 0.865, which means that for the problem without overhang constrains, the ratio of well supported masks of the MBB example is 86,5%. Increasing just a bit the value of Φ 0 creates some meaningful changes in the inclinations of some of the members while the rest of the geometry remains very alike. This evolution can be noticed in Fig. 13a, b, c and d where overhangs are gradually corrected as the constraint is tightened with higher values of Φ 0 .
Comparing Fig. 13a and b it can be noticed how elements and connection points are moved when Φ 0 is increased. The most meaningful modifications are the relocation of the knots and the shortening of the horizontal boundaries. For the middle knot, the algorithm divides and moves it up and downwards giving a more appropriate inclination to the boundaries of the hole beneath and introducing a new member while not sacrificing other boundaries. Also it can be seen that the non acceptable overhangs are reduced by making the horizontal down faced contours shorter. Another noticeable change is the slight modification of the rest of the knots, so that the amount of badly inclined members is reduced by modifying their inclination. For tighter values of the parameter Φ 0 , it can be seen that the number of holes in the structure decreases and their geometry changes revealing triangular shapes in most cases (see Fig. 13c and d). This behavior is not surprising since triangle-shaped holes can manage and reduce effectively the over-inclination of contours, so they turn out to be a natural choice for minimizing the number of badly supported elements. In fact, the top corner of the triangular hole is the only area where problematic contours can exist, since the top element could lead to a false-positive of the overhang constraint due to a single up pointing density gradient. The role of the parameter Φ 0 becomes very significant in this case because setting a value lower than unity helps the algorithm to skip false-positive contours and converge correctly. Concerning the value of the objective function, despite the deviation of the solutions with respect to the non-constrained optimum topology, the final geometries show compliance values that are close to the original solution of the MBB problem; obviously, the non-constrained optimization results always in a better structural efficiency.
Theoretically speaking, the ratio of the value objective function associated with the optimized design considering selfsupport property and that associated with the optimized design without considering self-support property should approach one . Actually, this measure can be used to evaluate the performances of optimized designs obtained by numerical approaches, so we will use the well known MBB problem in this section to make some comparisons with use of this measure. Assuming the load is 1 N and the sides of elements are equal to 1 mm, the values of compliance for the solutions shown in Fig.  13 are: a) 203.81 N.mm b) 208.33 N.mm c) 227.05 N.mm and d) 230.95 N.mm. The value of the compliance increases as the overhang restriction is tightened, as it was expected. The ratios with respect to the compliance that corresponds to the optimized design without considering self-support design are: a) 1.00 b) 1.02 c) 1.11 and d) 1.13. These values are close to one so we can consider that the performance of the optimized designs is acceptable and the proposed approach is an effective alternative for self-supported design.

Bridge structure
Finally we will consider the example shown in Fig. 14. This structure is another major challenge for the proposed algorithm as the non-constrained optimum topology shows several  boundaries where the inclination exceeds 45°. The rectangular design domain is discretized using 160 × 80 finite elements and it is supported in the lower right and left corners. A load is applied vertically at the middle point of the lower edge. The volume restriction is 30% of the domain and overhang constraints of 45°and 60°are introduced.
This example is used also to evaluate the effect of the parameter T, the threshold of the Heaviside Projection that defines the value of the design variable where the jump of the approximated step function takes place. Two different values of T are used, representing low and high values of the parameter, T = 0.01 and T = 0.8, respectively. Theoretically, low values would guide the design through initially denser domain, opposite to high values. Figure 15 is a source of information for discussing these effects on final designs. It can be appreciated in Fig. 15c, where a very low value of 0.01 is taken for T, that a vertical solid member is introduced in substitution of the central hole giving support to the horizontal boundary at the top of the structure. The rest of the members sprout from this central member with an inclination according to the constraint introduced. It is also worth to mention that the height of the structure is reduced so that more material can be placed in the inner zone and also to generate the vertical center member. This generates shorter and thicker members that make the topology ready to be fabricated via additive manufacturing but increases the value of the compliance losing part of the original structural efficiency.
On the other hand, Fig. 15b shows a bifurcation of the previously mentioned vertical member forming one vertical and two "v" shaped thinner members. This solution is closer to the unconstrained classical optimum design and provides the structure with greater stiffness than the previous solution. The choice for a value of T seems to be unimportant from a manufacturing point of view and for every value designers can reach perfectly printable topologies. However, considering the mechanical behavior, a different value can provide the structure with greater stiffness.

Conclusions
The method presented in this paper offers direct control over the inclination of structural members in topology optimization problems, so that self supporting, print-ready designs can be generated for additive manufacturing. The proposed approach presents some new features compared with previous works in this area. The inclination of members is calculated by an efficient edge detection algorithm developed in the field of image processing and applied here to evaluate the overhang angle of contours. This paper proposes also a new constraint function to control the amount of self-supporting contours, defined as the ratio between the admissible and the total amount of contours. That constraint can be easily added to the conventional volume requirement of density-based topology optimization procedures and combined with any general purpose optimizer like the popular MMA method. This technique has been combined with the commonly applied Heaviside projection to enhance black and white designs. Results obtained from the numerical examples suggest that the methodology developed allows the designer to efficiently control the amount of self-supporting members so that unprintable designs with infeasible overhanging sections can be banned from the design space, producing self-supporting topologies that are ready for additive printing. Moreover, the examples included in the paper demonstrate that the proposed strategy is effective at generating designs that can be printed without additional supports for any critical overhang angle defined by the user. Even if important aspects such as deformation of members due to overheating and residual stresses have not been included, authors expect that the approach suggested in this work could be of considerable practical value as a first approximation in early design stages. Future works include an implementation of the method for 3D structures optimization, where it would be especially valuable.