A stress-based criterion to identify and control intersections in 2D compliance minimization topology optimization

Topology optimization typically generates designs that exhibit significant geometrical complexity, which can pose difficulties for manufacturing and assembly. The number of occurrences of an important design feature, in particular intersections, increases with geometrical complexity. Intersections are essential for load transfer in many engineering structures. For certain upcoming manufacturing processes, such as direct metal deposition, the size of an intersection plays a role. During metal deposition, slim intersections are more prone to manufacturing defects than bulkier ones. In this study, a computationally tractable methodology is proposed to both control occurrence and size of intersections in topology optimization. To identify intersections, a stress-based quantity is proposed, denoted as Intersection Indicator. This quantity is based on the local degree of multi-axiality of the stress state, and identifies material points at intersections. The proposed intersection indicator can identify intersections in both single as well as multi-load case problems. To detect the relative size of intersections, the average density in the vicinity of an intersection is used to penalize or promote intersection sizes of interest. The corresponding sensitivity analysis involves solving a set of adjoint equations for each load case. Numerical 2D experiments demonstrate a controllable reduction of penalized slim intersections compared to the designs obtained from conventional compliance minimization. The overall geometrical complexity of the design is reduced due to the promotion of bulkier intersections which leads to an increase in compliance. The designs obtained are more suitable for manufacturing processes such as direct metal deposition.


Introduction
Topology Optimization (TO) is an effective method to generate early phase designs in a variety of engineering applications such as aerospace, biomedical, optics and maritime (Bendsoe and Sigmund 2013). Designs obtained by TO, although having a superior mechanical performance, can be difficult or impossible to manufacture due to their geometrical complexity (Lazarov et al. 2016). Additional manufacturing constraints are hence usually integrated into TO to generate designs which are suitable for processes such as casting (Xia et al. 2010;Zhou and Zhang 2019), machining (Sigmund 2009;Langelaar 2019) and Additive Manufacturing (AM) (Langelaar 2016;Gaynor and Guest 2016; Van de Ven et al. 2020). Many other studies have also provided methods to specifically reduce the geometrical complexity of TO designs. For instance, controlling minimum length scale in material and void regions. Length scale control can be achieved by filters or projection methods (Bourdin 2001;Guest et al. 2004;Zhou et al. 2015) and skeleton-based methods (Zhang et al. 2014;Xia and Shi 2015;Zhang et al. 2016Zhang et al. , 2017a. Moreover, geometrical complexity can also be controlled by restricting the number of structural components, in, e.g., the moving morphable component framework (Zhang et al. 2017b;Wein et al. 2020). In all mentioned methods a certain amount of Responsible Editor: Julián Andrés Norato Topical Collection: 14th World Congress of Struct Multidisc Optim. Guest Editors: Y. Noh, E. Acar, J. Carstensen, J. Guest, J. Norato, and P. Ramu. mechanical performance is typically sacrificed in order to achieve a manufacturable design.
Here we focus on generating designs by TO while controlling an important geometrical aspect of a structural design, namely intersections. The number of occurrences of intersections increase with the increase in geometrical complexity. Intersections are critical from both a performance and manufacturing point of view. However, how to explicitly control intersections is an open question in the literature. An intersection is a connection between two or more structural members in a design where a major change in the stress state occurs (Ambrozkiewicz and Kriegesmann 2018). Consequently, a multi-axial stress state is encountered at intersections (Bendsøe and Haber 1993). Intersections can be susceptible to failure due to stress multi-axiality (Clausmeyer et al. 1991). Moreover, from a manufacturing perspective, intersections can be associated with high cost. For structures produced by assembly operations, intersections may imply joining through riveting, bolting or welding operations (Megson 2019). Consequently, reducing the number of intersections in designs will lead to less joining operations, and lower cost. Intersections also pose a problem in certain AM methods. A prominent example is Direct Metal Deposition (DMD). DMD is a metal AM variant which is rapidly maturing and used currently to produce large functional structures such as ship propellers and aircraft components (Lockett et al. 2017). In DMD, simultaneous material melting and deposition occurs along the deposition lines. It is difficult to avoid overlapping of the deposition paths at intersections. This leads to extra material deposition at intersections which in turn cause geometrical deviations (Mehnen et al. 2014). However, when the size of an intersection is sufficiently large, overlapping of the deposition paths can be avoided. Hence, it is of interest to control not only the number but also the size of intersections.
Intersection detection in density based TO problem has been investigated through an image-based method (Gamache et al. 2018) and a stress-based method (Ambrozkiewicz and Kriegesmann 2018). In the image-based method, the skeleton of the structure is used to identify the structural members intersecting at nodes. These nodes are considered as intersections. However, the algorithm used for detection is not differentiable and therefore design sensitivities can not be calculated. This poses a problem in using this approach in gradient-based optimization (Gamache et al. 2018). Detection of intersections through a stress-based parameter provides a possibility for design sensitivity calculation. Stresses and sensitivities of stress-based quantities are readily computable for mechanical TO problems. Ambrozkiewicz and Kriegesmann (2018) used principal stress ratios to detect intersections and straight members connected to it, in order to evaluate the load path in the design. However, specific control on intersections during optimization was not performed and the use of these ratios was limited to detection in single load case problem.
In this paper, the control of intersections is studied using stress-state information in the context of the classical density based TO problem of compliance minimization (Bendsøe 1989). Controlling the number and/or size of intersections that are typically generated in TO may lead to major topological changes. To illustrate this, consider the schematic of a structural baseline design, shown in Fig. 1a. The intersections are marked by red dashed circles. The designs shown in Fig. 1b and c have the same material volume as the baseline design, but the intersections are larger and fewer. A significant reduction in the geometrical complexity is observed from Fig. 1a to b and from Fig. 1b to c. The aim of this paper is to present a computationally tractable approach for identifying and controlling intersections as a function of their size, in 2D density-based compliance minimization for single and multi load case problems. 2D problems are relevant for industrial DMD applications, because most of the structural components produced by DMD are manufactured by stacking 2D layers on top of each other. A discussion on challenges to extend the method to 3D is given at the end of the paper. Note that controlling intersections is shown for compliance minimization TO problems, however, the stress multi-axiality at intersections can also be used in many other applications.
In order to achieve the above mentioned aim, a mathematical formulation is developed to identify the intersections in the density based minimum compliance TO. There are two aspects to consider: detection and discriminating intersections depending on their relative size. In order to detect intersections, the stress state generated in the design due to the applied loads and boundary conditions is evaluated. The key idea is that at an intersection, the stress state Fig. 1 Schematic representation : a Complex baseline design with marked intersections (red dashed circles). b Design modification with reduced geometrical complexity and relatively thicker intersections than the baseline design. c Design modification with further reduced geometrical complexity and an even thicker intersection than in (b) exhibits an higher level of multi-axiality compared to uniaxially loaded design features. Therefore, a stress-based parameter is defined to measure the local degree of stress multi-axiality, making use of the principal stress ratios. An absolute value of this ratio close to 1 indicates a multi-axial state that occurs typically in the vicinity of an intersection. The stress-based parameter is selected because the sensitivity with respect to the design variables can be calculated easily. The mathematical description of the proposed Intersection Indicator is given for both single load and multi load problems. The latter is more challenging, because stress fields produced in response to each load system are different.
Once the intersections can be identified, a further step is to control the occurrence as a function of the relative size of the intersections. In order to identify relative size, the average density in a local circular domain at the location of an intersection is used as a measure of the size of the intersection. Consequently, intersections slimmer (or bulkier) in size can be penalized more, whereas bulkier (or slimmer) can be allowed in the design domain using a weight function. To illustrate how these measures can be combined, we propose an objective function aiming to minimize relatively slim intersections, as this is of particular interest for the DMD process.
During optimization, if the objective is dominated equally by compliance and intersection reduction, then black and white design will be realized. However, if the objective is dominated by intersection minimization, designs with considerable gray regions may result. This problem arises because, the intersection identification is indifferent to the magnitude of the local stress but solely depends on the local principal stress ratio. For this, a threshold stress level is introduced to distinguish between elements with similar principal stress ratios but different stress levels. Moreover, since the intersection indicator is a function of the principal stress ratio, if one or both principal stresses vanish during optimization the intersection indicator will be ill-defined. By introducing the threshold stresses into the formulation, this problem is also resolved. In general for stress-based TO, problems related to the local nature of the stresses and possible stress singularity are of concern (Kirsch 1990;Duysinx and Bendsøe 1998;Bruggi 2008). However, these problems are encountered when the stresses are constrained in the design domain. In our case, no constraints are imposed on the stresses, hence, reported problems are not encountered in our work. The effectiveness of our formulation is demonstrated, and the influence of various numerical parameters is studied based on several 2D compliance minimization problems.
The remainder of the paper is organized as follows: Sect. 2 presents the formulation of the Intersection Indicator, influence of threshold stress levels and intersection size estimation. The topology optimization problem formulation to minimize slim intersections and influence of threshold stress during optimization are presented in Sect. 3. Description of the numerical examples and resulting optimized designs and discussion are given in Sect. 4. A brief discussion on challenges to extend the method to 3D problems are outlined in Sect. 5. Finally, Sect. 6 provides the final conclusions.

Optimization problem
For completeness, we briefly summarize the optimization problem considered: In Eq. (1), c is the structural compliance, e represents the element nodal degrees of freedom, is the element stiffness matrix, and E e is the Young's modulus which is scaled using the filtered density ̃e as E 0 and E min are the Young's modulus of material and void, respectively, p is the SIMP penalization exponent and N is the total number of elements in design domain N . Eq. (2) defines the density filter applied to the design variable e at position with element volume v e (Bruns and Tortorelli 2001). Eq. (3) represents the weight ( t i ) calculation for the density filter. min e is the local circular region in which the filter is effective with a radius r min . Eq. (4) is the linear elastic state equation, where is the global stiffness matrix and, and are the global nodal degrees of freedom and nodal loads, respectively. Eq. (5) represents the volume constraint, and Eq. (6) bounds the density variable for all elements in the design domain N . The 2D problem is solved using (4) s.t.

Single load case
In order to identify intersections in compliance minimization TO, multi-axiality of the stress state is exploited. In compliance optimization the TO generates structural features typically under uniaxial stress to bear the load, such that the principal stress direction align in accordance with the orientation of a member (Bendsøe and Haber 1993;Ambrozkiewicz and Kriegesmann 2018). Since an intersection is where members with different orientation meet, the stress tensor should exhibit a multi-axial characteristic for single load case. This notion is illustrated using an optimized cantilever beam layout shown in Fig. 2. The optimization is carried out using the default parameters shown in Table 1, for the design domain depicted in Fig. 2a along with loading and boundary conditions. For the design obtained by standard compliance minimization, as depicted in Fig. 2b, two principal stresses I and II are calculated and the directions of the dominant principal stresses, i.e. either I or II , for the material region are realized by line segments in Fig. 2c. The length of the line segments in Fig. 2c are scaled with the magnitude of the principal stresses. As expected, the stress state is primarily uniaxial in the members except at the intersections. It remains to quantify stress multi-axiality in a design. The stresses in each element can be computed as follows: Note that the displacement field is already calculated in each optimization iteration. Stresses can therefore be obtained for a low computational cost. Here, e is the stress state of element e calculated at its centroid, e is the strain-displacement matrix and e is the constant constitutive matrix for isotropic linear elasticity. A plane stress condition is considered, i.e. e is given by [ xx , yy , xy ] T using Voigt notation. Next, the principal stresses for element e are calculated using: To measure stress biaxiality, we introduce a quantity R e composed of the ratios of these principal stresses The above function is chosen because it involves both ratios I ∕ II and II ∕ I . The stresses are squared because principal stresses could be both negative as well as positive and here we are only interested in the magnitude of stresses. For a biaxial stress state with principal stresses of roughly equal magnitude we obtain R e ≈ 2 . In cases where the magnitude of the principal stresses differ significantly, including uniaxial stress states, it follows that R e >> 2 . In the void regions, where the magnitude of the principal stresses is relatively small, R e may give false positives. Therefore, the filtered design variables ̃e can be used to avoid incorrect intersection detection in void regions. Thus, an Intersection Indicator ( e ) is proposed as: The operation 1∕ log 2 (.) is applied on R e to map the value of e between 0 and 1. As a result, three distinct settings can be identified: biaxial stress states in the solid ( e ≈ 1 ), uniaxial stress states in the solid ( 0 < e ≪ 1 ) and void regions ( e ≈ 0 ). Contours of this Intersection Indicator of the cantilever beam presented in Fig (9) regions. Thus, through exploiting biaxiality of the stress state intersections can be accurately identified.

Influence of threshold stress on Intersection Indicator
As mentioned previously, e is the intersection indicator which is a function of filtered density variable ̃e and stress multiaxiality measure R e . The intersection indicator is used for optimization but when the design focus is primarily on intersection reduction, gray regions in mildly stressed regions remain. This is because the multi-axiality criterion R e is indifferent to stress levels. Therefore, the optimizer tends to reduce the multi-axiality of the local stress state instead of reducing the density of the element. Two elements which are at different stress levels but with comparable principal stress ratios will give a similar value of R e . Therefore, it is desirable to distinguish between elements at different stress levels during optimization. Although void regions that are not load bearing are filtered out by the ̃e dependence of Eq. (11), the indifference of stress multi-axiality to differing stress levels can result in gray densities remaining in the design. We introduce a threshold stress to distinguish between points at different stress levels to mitigate the complications discussed above. Moreover, to define intersection indicator for the cases where principal stresses may vanish, i.e., I ≈ 0 or II ≈ 0 or both, the threshold stress is used. The multi-axiality criterion, R e , is thus modified to account for multi-axiality of the stress state only above a threshold stress level 0 : For 0 = 0 , H e is equal to R e . The effect of threshold stresses on the multi-axiality criterion is shown in Fig. 3. For a nonzero 0 , the multi-axiality criterion will behave almost identical for points in the domain where stress levels are significantly above 0 (compare cases | 0 ∕ I | = 0.05 and | 0 ∕ I | = 0 ). However, regions where both principal stresses are significantly smaller than the threshold stress level (see for example the case | 0 ∕ I | = 10 ) will be considered artificially as biaxially stressed regions even when | II | is very different from | I | . Consequently, for mildly stressed areas the multi-axiality criterion becomes insensitive to the modification of the stress state. Replacing R e with H e in Eq. (11) provides the modified Intersection Indicator ( I e ): The effect of 0 is demonstrated in Fig. 4. In this figure, 1∕ log 2 (H e ) and I e are plotted for values of threshold stresses equal to 0, 10 MPa and 100 MPa for the design shown in For 0 = 0 , in Fig. 4a it can be seen that the low-stressed regions, which are the void regions, exhibit both uniaxiality as well as biaxiality of the stress state. Since the stress levels are insignificant in these regions, the notion of multiaxiality of the stress state does not signify presence of any relevant design features at these locations. However, in the solid region the multiaxiality parameter clearly identifies the straight members and intersections via an uniaxial or biaxial stress state, respectively. Note that the magnitude of principal stresses in the void areas are orders of magnitude smaller than those in the solid area, but still the multiaxiality criterion does not differentiate between them. Scaling the multiaxiality criterion with the local density yields the intersection indicator shown in Fig. 4d. For the case 0 = 10 MPa , the multiaxiality criterion is shown in Fig. 4b. It can be observed that the void regions in which the stress levels are negligible are now artificially identified as being in the biaxial stress state. Scaling with the filtered design variable removes these artificial biaxial regions. The resulting intersection indicator plot shown in Fig. 4e, shows no sign of artificial biaxial regions and no effect in the material regions. Hence, through the introduction of the threshold stresses, for mildly stressed regions multi-axiality becomes constant and thus cannot be influenced by the optimizer. The effect of threshold stress during optimization for intersection minimization of the same problem is discussed in detail in Section 3.2. More importantly, most solid regions experience stress levels far above the threshold stress levels, thus the effect of the threshold stress is not visible in the solid region. This observation will be used in Section 3.2 to select the threshold stress for different considered problems. A higher threshold stress value, like 0 = 100 MPa , affects both the multiaxiality criterion and the intersection indicator, such that the material regions are also considered in the artificial biaxial stress state as shown in Fig. 4c and f, which should be avoided. Therefore, threshold stresses should be lower than the principal stresses in the solid region.

Multi-load case
The Intersection Indicator defined above can be easily extended to multi-load case problems. First, a slight modification in notation is required for a multi-load case problem. Analogous to Eq. (10), for a total of M load cases with Fig. 3 The effect of the threshold stress on the multi-axiality criterion. The x-axis represent the stress ratios assuming | II | ≤ | I | , and y-axis represent the corresponding value of the multi-axiality criterion. The multi-axiality criterion is plotted for various ratios of threshold stress and first principal stress e ) corresponding to the m th load case. Fig. 5a shows the design domain of a cantilever beam with prescribed boundary and loading conditions. Note that, P 1 and P 2 are applied as separate load cases, corresponding to m = 1 and m = 2 , respectively. Through standard topology optimization for multi-load case problem, the design shown in Fig. 5b is obtained. The default parameters used for optimization are given in Table 1. The average | I | and | II | in the domain corresponding to only load P 1 are 26 MPa and 24 MPa , respectively, and corresponding to only load P 2 are 9 MPa and 7 MPa , respectively. Therefore, the threshold stresses (m) 0 corresponding to only load P 1 and P 2 are selected as 10 MPa and 3 MPa , respectively, which are lower than both average principal stress values. The biaxiality criterion 1∕ log 2 (H (m) e ) for load case m = 1 and 2 are shown in Fig. 5c and d, respectively. The corresponding intersection indicator I (m) e for load case m = 1 and 2 are shown in Fig. 5e and f, respectively. It is evident that through I (1) e all intersections of the design can be identified. However, I (2) e shows several false positives in the structural members (for example, Point d 2 in Fig. 5f). This is due to the fact that part of the structure on the right side of the load P 2 is not load-bearing for the load case with P 2 only. Therefore, the stresses generated on the right side are negligible. Thus, false multi-axiality is detected in the unloaded part of the structure. To account for all the load cases and mitigate false positives, an aggregation scheme is applied considering the following statements: -If an Element e has an uniaxial stress state and is load bearing in all load cases, then for all load cases, | (m) I | ≫ | (m) II | or vice versa, leading to high H (m) e values in each case. In the example problem depicted in Fig. 5, Point d 1 is load bearing for both load P 1 and P 2 . The values of H (m) e for Point d 1 in both load cases are considerably higher than 2 which can be inferred by Fig. 5c- stress state for several loadcases which is not load bearing in other loadcases. Then, for loadcases in which it is load bearing | (m) I | ≫ | (m) II | or vice versa, again leading to high H (m) e values in each case. However, for loadcases in which it is not loaded | (m) e will be approximately equal to 2. In the example problem, Point d 2 is uniaxially loaded for loadcase P 1 but is not loaded for P 2 . Therefore, the value of H (1) e is high and H (2) e ≈ 2 as shown in Fig. 5c e over all the loadcases will hence result into a value approximately equal to 2M.
Motivated by above observations to formulate the intersection indicator for multi load case problems, the biaxiality in the structure is aggregated as follows: Through this definition the uniaxially loaded elements in at least one loadcase will exhibit value of H e ≫ 2 . Moreover, the biaxially loaded elements at least in one loadcase will exhibit value of H e ≈ 2 . Now, substituting Eq. (16) in Eq. (13) provides intersection indicator ( I e ) for the multiload case problem. For the example problem, the plot of 1∕ log 2 (H e ) and corresponding intersection indicator is shown in Fig. 5g and h, respectively. It can be observed in Fig. 5g that through aggregation the uniaxially and biaxially loaded elements in the solid are clearly identified. The void regions are considered in artificially biaxial state which are then filtered out through filtered density as shown in corresponding intersection indicator plot Fig. 5h. This formulation is also applied to problems with more than 2 load cases to check the scope of the approach. It is observed that through proposed aggregation, intersections can be detected in multi-load case problems with more than two load cases. We note in passing, for M = 1 , Eq. (16) reduces to Eq. (12).

Intersection Size
As discussed in the Introduction, in certain applications one may require to control the occurrence of intersections depending on their sizes. Consider the two intersections, shown in Fig. 6a and b. To evaluate the intersection size, assume a circular domain, int e , centered at Element e (indicated by a red dot in Fig. 6). Averaging the filtered densities in this local domain provides information on the size of the intersection at Element e. For a relatively bulky intersection, shown in Fig. 6b, the average filtered density is higher than in case of a relatively slim intersection, shown in Fig. 6a. However, the intersection shown in Fig. 6c, which happens to be an intersection as slim as the intersection depicted in Fig. 6a, has a higher average density in the local region int e compared to that of the intersection shown in Fig. 6a. Note that, the farther the material is located from the center of Multi-load case Cantilever Beam : a design domain with specified size, loading and boundary conditions, b design ( ̃e ) obtained by solving standard topology optimization compliance minimization problem, c and d represents the biaxiality criterion ( 1∕ log 2 (H (m) e ) ) for each load P 1 and P 2 , respectively. e and f shows the Intersection Indi-cator ( I (m) e ) field of the corresponding design considering load case P 1 and P 2 , respectively. (g) represents the aggregated biaxiality criterion ( 1∕ log 2 (H e ) ). H e is calculated using Eq. (16). (h) is the corresponding intersection indicator field obtained by aggregating the biaxiality of the stress state experienced by the structure by load P 1 and P 2 the circular domain at Element e, the lower the probability this material is part of the intersection. To account for this, a weighted averaging scheme is proposed, such that the contribution of the filtered density of a particular element to the local averaged density at a point of interest depends on the distance between the element and the point of interest. The weighted density filter (Bruns and Tortorelli 2001) is used to calculate the local weighted average of the filtered density(̂e): Here, r int is the radius of the local circular domain int .
A straight member, as shown in Fig. 6d, has values of ̂e equivalent to that of a bulky intersection. To distinguish between a straight member and an intersection, this size measure ̂e must be combined with the Intersection Indicator (I e ).
Inspired by the requirements in the DMD process, we consider the case of allowing bulky intersections, while suppressing slim ones. Thus, to penalize slim intersections, weight factors are calculated as function of intersection size (18) s e = 1 −̂e, Here, ŝ e is the weight assigned to the intersections at Element e and S e is the weighted intersection indicator that is higher for slim intersections in a design. These weights are chosen as an example for the DMD application, however, the weights can be chosen in other ways to target different intersection sizes for other applications.
Recall that the minimum member size in TO is set to be 2r min therefore, the minimum intersection size is 2 √ 2r min . Thus, to detect the size of an intersection, we require r int > √ 2r min . Consequently, intersections larger than the chosen value r int are not detected. To illustrate this with an example, we revisit the cantilever problem shown in Fig. 2. The filter given in Eq. (17) is applied to the filtered density field shown in Fig. 2b for three different values, where = r int ∕r min , as shown in Fig. 7. The chosen values of are 2, 4 and 6. As observed from Fig. 7a-c, increasing the radius increases ŝ e in entire domain, and thus, a greater number of thin intersections are identified below the domain size int e . This is evident from Fig. 7d-f. Note that the weighted averaging in the local domain provides a relative size measure which can distinguish between bulky and slim intersections.

Multi-Objective Formulation : Compliance and Intersection Objective
So far it has been established that the presence of intersections can be identified and their relative size can be determined. It remains to control the occurence of intersections in TO. The intersection indicator defined in the previous section only gives an indication of intersection presence. It does not provide insight on the actual number of intersections. Moreover, as stated in the Introduction, intersections facilitate load transfer, therefore their presence is paramount for mechanical performance. However, from a manufacturing point of view, such as for DMD, slim intersections are problematic. Therefore, the number of slim intersections should be restricted to improve manufacturability of through DMD. Therefore, there is a trade off between performance and manufacturability. Consequently, for optimization we adopt a multi-objective approach where a total objective function is defined as a combination of compliance of the structure and a function suppressing slim intersections. Firstly, we define the function that will be minimized to achieve fewer slim intersections for DMD application.  The proposed function which is also termed as intersection objective I is defined as: The numerator of the intersection objective I is the summation of the weighted intersection indicator S e over the entire domain. Given Eq. (18), the contribution from relatively slim intersections will be higher, whereas bulky intersections and straight members will only contribute marginally. Thus, minimizing the numerator will minimize thin intersections and promote bulky intersections or straight members in a design. The denominator of the intersection objective is the design volume normalized by the allowed volume in the design domain V 0 . Upon omitting the term, the optimizer removes material from the intersection location and the optimization tends towards a trivial solution of no material. Therefore, the denominator term promotes material in the design domain.
In Eq. (21), O is the total objective function, is the weight assigned to the compliance objective, consequently, (1 − ) is the weight assigned to the intersection objective. The parameter is a tuning parameter which can be reduced to emphasis on manufacturability of the design through DMD. Both the objectives are normalized by their respective values of the compliance ( c * ) and intersection ( I * ) calculated for the design obtained by the standard compliance minimization problem. Superscripts ( * ) are used to indicate values that are calculated on the converged standard compliance minimization design. The constant factor is introduced to scale the value of the objective such that it ranges between 1 and 100 for the stability of the MMA optimizer, see (Svanberg 1987). The sensitivity analysis of the objective and constraint functions are determined using the adjoint method and detailed in Appendix A.

Effect of threshold stress levels on Intersection Minimization and Selection Strategy
As mentioned in Section 2.2.2, the effect of the chosen threshold stress during optimization is crucial. To demonstrate its effect, initially the threshold stress is assumed to be zero ( 0 = 0 ). The multi-objective formulation is applied to Problem 1 with = 0.2 , i.e, the contribution of intersection objective I is dominant. The pure compliance minimization design previously shown in Fig. 2b is chosen as the reference case for the calculation of c * and I * . The value of = 6 is chosen because it targets every intersection in a compliance minimization design, shown in Fig. 7f. The results obtained after 100 iterations are presented in Fig. 8. In the design shown in Fig. 8a, it can be observed gray regions d, e and f represent the corresponding weighted intersection indicator S e =ŝ e I e representing the intersections smaller than the threshold defined by r int . The red and black circles represent the size of the r min and r int , respectively. As radius r int increases more number of intersecting features are identified are present at locations which experience low stress levels, such as Region 1 and 2 (indicated by red circles) and an arch shaped design feature is observed with high material density in Region 3. The corresponding objective values are c∕c * = 2.01 , I∕I * = 0.44 and O∕O * = 0.75 , see Eq. (21). It shows that the intersection objective and total objective have reduced, as desired. However, the obtained result does not represent a desired, manufacturable geometry. The features formed in the low stresses regions demonstrate that the effect of reducing the local multi-axial stress state dominates removal of local densities in these regions, as shown in Fig. 8b, leading to the depicted undesired result. In order to facilitate only removal of local densities from low stresses regions, the effect of reducing local multi-axial stress state should be nullified.
To avoid generation of the design features in regions where the stress is negligible during optimization, the threshold stress, mentioned in Section 2.2.2, is introduced. The stress response from the pure compliance minimization reference design is used to aid in the selecting the threshold stress level. Recalling that for the compliance minimization design in Fig. 2b average and maximum value of both | I | and | II | in the domain are 25 MPa and 618 MPa , respectively. Also, it is shown previously that for 0 = 10 MPa , which is below the average stress value, does not affect the intersection indicator of the converged compliance design (see Fig. 4d and e). Therefore, for the intersection minimization problem, 0 = 10 MPa is selected because in the reference compliance minimization case, the stresses in the solid region are significantly above 10 MPa . The result of intersection minimization after 100 optimization iterations considering 0 = 10 MPa is shown in Fig. 8c. The design does not contain gray design features in the lowly stressed region, as they are artificially considered in the biaxial stress state which can be observed in Fig. 8d. Since, there is no possibility to reduce the local multi-axial stress state in the low stressed region, the intersection objective is reduced by removing the local densities in these regions, which promotes convergence to a manufacturable designs. The corresponding objective values are c∕c * = 1.34 , I∕I * = 0.35 and O∕O * = 0.55.
It remains to determine the value of the threshold stress for general cases. The key aspect is that low-stressed regions should be discouraged to affect the detection of stress multiaxiality. Since the intersection minimization problem is solved simultaneously with the compliance minimization problem, the stress levels found in the solid regions of the compliance minimization design are considered as representative. Therefore, the threshold stress value should be selected such that it is smaller than the stress levels present in the solid regions of the compliance minimization design. Now, the stress levels in the solid regions are either dictated by I , II or both. Thus, the measure of stress levels should be combination of the individual principal stress components. This motivates the selection of Von Mises stresses as a measure of stress levels in solid regions and the threshold stress should be smaller than the Von Mises stress field in the solid regions of compliance minimization design. The compliance minimization design consists of void (white) and material (black) regions separated by interface (gray) regions. The stress levels experienced by interface (gray) regions fall between those in void and material regions. In this paper, the following procedure is adopted to select the threshold stresses for single and multi loadcase problems: 1. Apply the load and boundary conditions to the converged design obtained from the standard compliance minimization problem and determine the principal stresses at each element in the domain using Eq. (8) and Eq. (9). 2. Calculate the Von Mises stress ( (vm) ) at each element in the domain using √ 2 I − I II + 2 II . 3. To identify the Von Mises stress in the gray regions, stresses which corresponds to filtered design variables 0.1 <̃e < 0.5 are selected. Due to SIMP penalization, the stress levels corresponding to the value of ̃e = 0.1 and ̃e = 0.5 will be in the order of 0.1% and 12.5% compared to the stress levels in the solid region ( ̃e = 1 ), respectively. Therefore, the stress levels in the selected range will be significantly lower than the stresses in the solid region. 4. To select the threshold stresses as low as possible, sort the Von Mises stresses of the selected filtered density range in ascending order. Pick the stress value found at an index closest to 90% of the length of the array. The 90% value is selected as threshold so that most of the lowly stressed regions in a design are artificially considered in biaxial stress state.
The proposed algorithm used to systematically select threshold stresses presented a good performance distinguishing between the lowly and significantly stressed regions during optimization for all problems discussed in this paper. Nevertheless, note that the proposed algorithm is only one option and choosing the best algorithm is beyond the scope of this work.

Results
In this section, results obtained by applying the proposed multi-objective formulation on various single and multi load problems are presented. The design domains and boundary conditions of the 5 considered test problems are given in Fig. 9. Each problem is labeled with a problem number for reference. All the design domains have equal dimensions.
Recall that the default optimization settings are given in Table 1. Optimization results of all the problems and effect of different parameters on the designs are shown in Section 4.1. An overall discussion is given in Section 4.2.

Effect of
The designs obtained for all of the test problems defined in Fig. 9 are shown in Fig. 10. Results obtained for = 6 are shown. Here, r int is chosen to be sufficiently large to target every intersection in the respective compliance minimization design, further effects of r int on the design are discussed in Section 4.1.2. Converged designs obtained for different values of and their corresponding intersection and compliance objective values are shown in Fig. 10. Designs for = 1.0 correspond to standard compliance minimization, the reference design for each of the problems. It is evident from the designs that, as the weight on the intersection objective increases, the number of intersections reduces in the design.
Consistently the compliance objective is increasing while the intersection objective is decreasing. In all the problems, it is observed that the minimization of intersections comes at the cost of inferior stiffness performance.

Effect of r min and r int
Recall that the proposed method utilizes two local circular domains, determined by the radii, r min and r int , that are anticipated to affect the optimized design. Increasing r min leads to an increase in the minimum member size and, thus, leading to fewer intersections. In order to isolate the effect of r int , r min is kept constant and r int is varied. In Fig. 11, designs for Problem 1 are illustrated for three selected values of = r int ∕r min equal to 2 ,4 and 6, while keeping r min = 2.5 mm and = 0.6 . The values of compliance and intersection objective given in Fig. 11 are normalized by the reference design shown in Fig. 2b. For all the designs, the intersections are clearly identified by the intersection indicator as shown in the corresponding plot of I e , see Fig. 11. During optimization, if an intersection size reaches the size of the local domain with radius r int , then the contribution of that intersection to the intersection objective will be negligible and therefore, the intersection will not be reduced. This effect is clearly visible in the design for = 2 , as shown in Fig. 11. The weights ŝ e , which is inversely proportional to intersection size (see Eq. (18)), corresponding to each intersection are close to 0. It suggests that the intersection size is equivalent to the local domain with radius r int . Also the contribution of the intersections in the intersection objective is close to 0 as evident from the S e plots.
As r int increases, bulkier intersections are preferred which can be seen in the design for = 4 , see Fig. 11. The weight ŝ e corresponding to the bulky intersection is close to 0 as can Page 13 of 19 307 be seen in the ŝ e plots. Furthermore, the contribution from bulky intersections to the intersection objective is also close to 0 as shown in the S e plots.
Further increasing r int , i.e. for = 6 , relatively thicker intersections are observed than for = 2 . However, a slim intersection is also observed which is not removed during optimization. Plot of weight ŝ e shows that the bulky intersection has relatively lower value of ŝ e than the slim intersection. The contribution from slim intersection to the intersection objective is also higher as shown in S e plot. It means that the slim intersections are detected, however, it is not removed during optimization. This happens because, firstly, the slim intersections are minimized which does not mean that slim intersections can not exist in a design. Secondly, using = 0.6 means that the compliance objective function has a greater influence, therefore, the slim feature is not removed completely.

Discussions on Results
The results show that by minimizing the slim intersections, the topology of the compliance minimization design has been simplified. Moreover, calculation of the intersection indicator is computationally efficient because at every iteration step displacement information is available, and calculation of stresses based on this requires only element level linear operations (see Eq. (8)). However, for minimizing the intersection objective an additional set of adjoint equations has to be solved which increases its computational cost. For the multi-load case problem, the number of adjoint equations scales with the number of load cases in a single problem because stress fields from different load cases are aggregated together. The convergence behavior of normalized compliance and intersection objective with respect to iteration step of Problem 1 with = 6 are shown in Figs. 12 and 13, respectively. The convergence plots show smooth convergence. However, the number of iterations required was quite high compared to standard compliance minimization problem. This is due to the fact that the sensitivity of the compliance minimization problem is always negative with respect to every design variable, but for the intersection objective it is not the case, thus the convergence gets slower compared to compliance minimization.

Extension to 3D Problems
For 3D problems, again there are two aspects to be considered, identification and control of 3D intersections. In contrast to the 2D problems, it was observed, by examining 3D TO results, one may identify following types of intersections: -Straight beams intersecting with each other.
-Straight beams intersecting with a plate type structure.
-Two or more plate like structures intersecting with each other.
To identify above mentioned 3D intersections certain 'signatures' in terms of principal stress ratios have been recognized. However, defining an indicator, similar to the 2D case, to reliably identify these features proves more challenging due to the larger variety of intersection types. Further investigation is needed to develop an intersection indicator for 3D problems. Moreover, controlling 3D intersections to generate designs for certain manufacturing process such as DMD have two important aspects. Firstly, 3D intersections need not be problematic for the DMD process because mostly the production happens in a layer-by-layer manner. Thus, it is possible that within the deposition layer intersections are not encountered and the part could be produced easily. A rather challenging problem left for future research is to identify the intersections within the deposition layers and then minimize them for DMD related applications. This could be achieved by projecting the stress components to the plane of deposition. Secondly, for DMD processes, designs with high geometric complexity could also pose a problem during manufacturing. Hence, once a reliable 3D intersection indicator is developed, then, based on their relative size, intersections can be controlled to achieve minimal geometric complexity, similar to the 2D case.

Conclusions
In this study a novel approach to reduce geometric complexity of designs in 2D compliance minimization topology optimization is presented. The reduction in the geometric complexity is attained through identifying and controlling intersections by means of the local stress state for single and multi load case problems. The reduction of geometric complexity of a design improves manufacturability through both conventional as well as additive manufacturing processes, such as DMD.
The proposed Intersection Indicator is a function of the ratio of local principal stresses. Therefore, the design sensitivities of the indicator can be easily evaluated by solving a set of adjoint equations per loadcase. This renders the indicator suitable for gradient-based optimization. To prevent convergence problems due to spurious stress ratios in low-stressed regions, a threshold stress is introduced as well as a procedure to select an appropriate value. To control intersections, a multi-objective formulation is proposed comprising of a compliance and intersection objective. A provision to emphasize intersections according to their size has also been presented, and demonstrated in the context of the DMD process, where slim intersections are to be avoided. Controlled by the relative weight of the intersection objective, different levels of design simplification are achieved, both for single as well as multi-load problems. As expected, increasing restrictions on design complexity leads to increased compliance.
Smooth convergence is observed without a continuation scheme. However, the number of iterations required for convergence increases when the relative weight of the intersection objective increases. Implementing the proposed formulation in an existing topology optimization framework, such as the 88-line MATLAB TO code presented by Andreassen et al. (2011), is straightforward. It involves calculation of principal stresses, threshold stresses and sensitivities of intersection objective.
Although the proposed approach proved effective in various test problems, it has limitations. First, it does not provide control on the exact number of intersections. Second, the stress-based indicator relies on characteristics of optimal designs obtained in compliance minimization. Its extension to other problems such as compliant mechanism or thermal optimization will require developing a relation between geometrical aspects of intersections and the 'physics' observed at intersections. Currently, the presented intersection indicator is therefore restricted to 2D compliance minimization problems. Extension to a 3D setting is a topic for future work, where stress-based recognition of the variety of intersection types that can occur in 3D is the main challenge. Fig. 10 Designs obtained after applying compliance and intersection minimization for considered problems. Arrows indicate designs for the corresponding data point. The parameter associated to a design is given below the design. In all the test problems, as the weight on the intersection objective increases slim intersections are reduced and bulky intersections are preferred, consequently, geometrical complexity is reduced ◂ Fig. 11 Effect of r int on Problem 1 : Cantilever Beam. The red and black circles indicate the size of the local circular domain with radius r min and r int , respectively. As the radius r int increases the size of the intersections promoted in the design domain increases. Therefore, the size of the intersections also increases in the converged designs with increase of r int The sensitivity of the intersection objective is calculated by applying the quotient, product and chain rules of derivatives. The sensitivity calculation is shown for a single loadcase problem for clarity and it can be extended to multi loadcase problem with minor modification. Firstly, design differentiation of Eq. (20),yields: The term ̃ẽ j = ej , where ej is the Kronecker delta function, thus the above equation reduces to: Now the product rule is applied to the later part of the above equation using Eq. (18), as shown below: In the above equation, the term ̂ẽ j , can be derived from Eq.
(17) as shown below: Now, since filtered densities are independent variables thus, the above equation can be written as: To calculate the term with I ẽ j in Eq. (29), the product rule is applied to Eq. (13). With some simplifications and using Eq. (8), Eq. (9), Eq. (12) and Eq. (13) the term can be reduced to the following form: where, The sensitivity calculation of individual terms is given subsequently as where, where,