Analytical relationships for imposing minimum length scale in the robust Topology Optimization formulation
PPREPRINT
Note on the minimum length scale and its defining parameters
Analytical relationships for Topology Optimization based on uniform manufacturinguncertainties.
Denis Trillet · Pierre Duysinx · Eduardo Fern´andez
Received: date / Accepted: date
Abstract
The robust topology optimization formula-tion that introduces the eroded and dilated versions ofthe design has gained increasing popularity in recentyears, mainly because of its ability to produce designssatisfying a minimum length scale. Despite its success invarious topology optimization fields, the robust formu-lation presents some drawbacks. This paper addressesone in particular, which concerns the imposition of theminimum length scale. In the density framework, theminimum size of the solid and void phases must be im-posed implicitly through the parameters that define thedensity filter and the smoothed Heaviside projection.Finding these parameters can be time consuming andcumbersome, hindering a general code implementationof the robust formulation. Motivated by this issue, inthis article we provide analytical expressions that ex-plicitly relate the minimum length scale and the pa-rameters that define it. The expressions are validatedon a density-based framework. To facilitate the repro-duction of results, MATLAB codes are provided.As a side finding, this paper shows that to obtainsimultaneous control over the minimum size of the solidand void phases, it is necessary to involve the 3 fields(eroded, intermediate and dilated) in the topology op-timization problem. Therefore, for the compliance min-imization problem subject to a volume restriction, theintermediate and dilated designs can be excluded fromthe objective function, but the volume restriction hasto be applied to the dilated design in order to involveall 3 designs in the formulation.
AuthorsDepartment of Aerospace and Mechanical Engineering,University of Li`ege, 4000 Li`ege, Belgium.Denis TrilletE-mail: [email protected]
Keywords
Length Scale · Robust Design · SIMP
Since the seminal work of Bendsøe and Kikuchi (1988),topology optimization has experienced huge advancesand nowadays is being massively adopted in the in-dustry (Pedersen and Allinger, 2006; Zhou et al, 2011;Zhu et al, 2016). Among the successful advancements,one can mention the famous density method under theSIMP interpolation scheme (Bendsøe, 1989). The well-known shortcomings of SIMP led to a succession of im-provements seeking to avoid the mesh dependency, thecheckerboard patterns and the presence of intermediatedensities. To date, one of the most effective approachesto dealing with the ill-effects of SIMP is the robust de-sign approach that brings the eroded and dilated ver-sions of the design (Sigmund, 2009). This method con-siders manufacturing errors that may incur in a uni-formly thinner or uniformly thicker component com-pared to the blueprint layout. Sigmund (2009) proposeda robust formulation that maximizes the performanceof the worst performing design among the eroded, di-lated and reference (intermediate) designs. This means,the formulation guarantees a design with good perfor-mance even if it is eventually eroded or dilated duringthe manufacturing process. Interestingly, the robust for-mulation yields a reference (intermediate) design thatfeatures minimum member size and minimum cavitysize (Wang et al, 2011).In addition to imposing minimum length scale intopology optimization, it has been seen that the ro-bust formulation provides a more stable convergencethan other projection or filtering strategies known sofar, which allows to reach almost discrete solutions in a r X i v : . [ m a t h . O C ] J a n Denis Trillet et al. the density method (Wang et al, 2011). The formula-tion has been recently applied in combination with ex-isting methods that allow to impose stress limits (daSilva et al, 2019; Silva et al, 2020), overhang angle con-straints (Pellens et al, 2018), maximum size restrictions(Fern´andez et al, 2020), and geometric non-linearities(Lazarov and Sigmund, 2011; Silva et al, 2020), not onlyin the density method but also in the level set method(Chen and Chen, 2011; Andreasen et al, 2020).Despite proving its effectiveness in several fields ofapplication (Wang et al, 2011b; Christiansen et al, 2015),the robust design approach presents some drawbackswith regard to its implementation. For example, themethod requires boundary treatments with respect tothe density filter, since the filtering region can be splitat the edges of the design domain affecting the imposedminimum size (Clausen and Andreassen, 2017; Kumarand Fern´andez, 2021). In addition, due to the erosionand dilation distances with respect to the intermedi-ate design, boundary conditions can be disconnectedfrom the designs involved in the formulation (Clausenand Andreassen, 2017). Another difficulty posed by themethod is that the minimum length scale must be im-plicitly imposed through the parameters defining thedensity filter and the Heaviside projection. This draw-back has been addressed in the literature using numer-ical (Wang et al, 2011) and analytical (Qian and Sig-mund, 2013) approaches. The numerical approach con-sists of applying the density filter and the Heaviside pro-jection to a 1D design, measuring the resulting lengthscale, and repeating the process several times with dif-ferent filter and projection parameters to subsequentlyconstruct a graph that relates the involved parameters.The analytical method consists of applying the densityfilter as a convolution integral in a continuous 1D de-sign. The solution of the integration yields explicit re-lationships of the projection/filtering parameters withthe obtained length scale (Qian and Sigmund, 2013).As the three field scheme is a scalar function, the re-lationships obtained from the 1D designs remain validfor 2D and 3D (Wang et al, 2011).The numerical and analytical relationships reportedby Wang et al (2011) and Qian and Sigmund (2013) areintended for the particular case where the minimum sizeof the void phase is set equal to that of the solid phase.This simplifies and facilitates the procedures to explic-itly relate the desired minimum length scale to the filterand projection parameters, but in this way, the mini-mum size of the solid phase cannot be different fromthe minimum size of the void phase. Nonetheless, dueto the rapid popularization of the robust formulation,it becomes necessary to find a method that allows toimpose different minimum sizes for each phase, in addi- tion to provide access to other geometric information.For example, the erosion and dilation distances with re-spect to the intermediate design have been required toimpose a maximum member size (Fern´andez et al, 2020)and to obtain designs tailored to the size of a deposi-tion nozzle (Fern´andez et al, 2021). To date, the erosionand dilation distances have been reported for specificminimum length scales, which hinders widespread ap-plicability of the above methods.The aim of this work is to provide a method to ob-tain the filter and projection parameters that imposeuser-defined minimum length scales, where the size ofthe solid phase can be defined differently from that ofthe void phase. In addition, for the desired minimumlength scales, the erosion and dilation distances are pro-vided. To this end, we extend the applicability of the an-alytical method proposed by Qian and Sigmund (2013)and the resulting relationships are validated using thenumerical method proposed by Wang et al (2011) and aset of 2D topology-optimized designs. To facilitate thereplication of results and the application of the methodsdiscussed in this paper, MATLAB codes are provided.The remainder of the document is developed as fol-lows. Section 2 introduces the robust topology opti-mization formulation and the case studies used to assesthe analytical expressions that relate the minimum len-gth scale to the parameters that impose it. Section 3presents the analytical method proposed by Qian andSigmund (2013) and describes the main contribution ofour work, which is the extension of the analytical ap-proach to allow choosing independent minimum sizesfor the solid phase and the void phase. Section 4 val-idates the analytical expressions using the numericalmethod proposed by Wang et al (2011). Section 5 dis-cusses the sources of error that are inherent to the an-alytical method. Section 6 assesses the analytical ex-pressions on 2D topology-optimized designs. Section 7gathers the final conclusions of this work while Section8 provides the codes that allow for the replication ofresults.
The analytical expressions that relate the minimumlength scale to the parameters that define it are devel-oped for the robust topology optimization formulationbased on the eroded, intermediate and dilated designs.This work considers the density approach based on theSIMP interpolation scheme (Bendsøe, 1989), even thoughthe proposed methodology can be applied to other topol-ogy optimization approaches with little efforts.Like most works in the literature, the eroded, inter-mediate and dilated designs that constitute the robust ote on the minimum length scale and its defining parameters 3 formulation are built using a three-field scheme (Sig-mund and Maute, 2013). The first field, denoted by ρ ,corresponds to the design variables. The second field,denoted by ˜ ρ , is obtained by a weighted average of thedesign variables within a circle of radius r fil . The thirdfield, denoted by ¯ ρ , is obtained by projecting the com-ponents of the filtered field towards 0 or 1. The filter andprojection functions are identical to those provided byWang et al (2011), however, these are reminded hereinfor the sake of clarity as they define the minimum lengthscale obtained in the optimized design.The filter of design variables, or the density filter(Bruns and Tortorelli, 2001; Bourdin, 2001), is definedas follows:˜ ρ i = N (cid:88) j =1 ρ j v j w ( x i − x j ) N (cid:88) j =1 v j w ( x i − x j ) , (1)where ρ j is the design variable associated to the element j and ˜ ρ i is the filtered variable associated to the element i . v j is the volume of the element j and w ( x i − x j ) is theweigh of ρ j in the definition of ˜ ρ i . As is common practicein the literature, the weighting function w ( x i − x j ) isdefined as a linear and decreasing function with respectto the distance between the elements i and j , as follows: w ( x i − x j ) = max (cid:18) , − (cid:107) x i − x j (cid:107) r fil (cid:19) , (2)where x i and x j represent the centroid of the elements i and j , respectively. It is recalled that r fil is the radiusof the density filter.To reduce the amount of intermediate densities pre-sent in the filtered field and to build the designs thatconstitute the robust formulation, the projected field isobtained with the following smoothed Heaviside func-tion (Wang et al, 2011):¯ ρ i = tanh ( βη ) tanh ( β (˜ ρ i − η ))tanh ( βη ) tanh ( β (1 − η )) (3)where β and η control the steepness and the thresholdof the projection, respectively. The eroded, intermedi-ate and dilated designs, denoted by ¯ ρ ero , ¯ ρ int and ¯ ρ dil ,are obtained from the smoothed Heaviside function ofEq. (3) by using the same β but with different thresh-olds η ero , η int and η dil , thus leading to ¯ ρ ero ( ˜ ρ , β, η ero ), ¯ ρ int ( ˜ ρ , β, η int ) and ¯ ρ dil ( ˜ ρ , β, η dil ).To assess the scope of the equations that provide thenecessary parameters to impose the minimum lengthscale, different topology optimization problems are sol-ved with variations in the desired minimum length scale. (a)(b) Fig. 1: (a) Heat sink and (b) force inverter design do-mains considered in this study.The minimum size of the optimized designs is measuredand compared against the intended values. In this work,the problems that are chosen to illustrate the develop-ments are the heat conduction and the non-linear com-pliant mechanism design. The former aims at minimiz-ing the thermal compliance subject to a volume restric-tion (Wang et al, 2011), whilst in the latter, the formu-lation of a force inverter is considered where an outputdisplacement is maximized for a given input force (Sig-mund, 1997). The design domains are shown in Fig. 1.According to the robust design approach, the topologyoptimization problems can be written as follows:min max ( c ( ¯ ρ ero ) , c ( ¯ ρ int ) , c ( ¯ ρ dil ))s.t.: v (cid:124) ¯ ρ dil ≤ V ∗ dil ( V ∗ int )0 ≤ ρ i ≤ , i = 1 , ..., N , (4)where c represents the thermal compliance of the heatsink or the output displacement of the force inverter, V ∗ dil is the upper bound of the volume restriction and Nis the total amount of design variables. As the designintended for manufacturing is the intermediate design,the upper bound of the volume constraint is scaled ac-cording to the user-defined limit V ∗ int . The optimizationproblem and the scaling process of the volume restric-tion are the same as described by Wang et al (2011)for the heat conduction problem. In the nonlinear forceinverter, the method to deal with the mesh distorsionappearing in low density elements is the same as ex-plained by Wang et al (2014). To avoid overextendingthe contents of the manuscript, the interested reader isrefered to the cited articles. Denis Trillet et al.
Comprehensive numerical tests have shown that the ro-bust formulation in Eq. (4) yields intermediate designsthat have the same topology as their eroded and di-lated projections (Wang et al, 2011). In other words,the erosion process does not destroy the solid membersof the optimized topology (intermediate), it just makesthem thinner. This implies that a minimum size of thesolid members in the intermediate design is imposedaccording to the erosion distance. Similarly, the dila-tion process does not close the cavities of the optimizedtopology (intermediate), it only makes them narrower.This implies that the minimum size of the void phaseis imposed according to the dilation distance.Considering a two-dimensional design domain, theminimum size of the solid phase is defined by the radius r min . Solid of the largest circle that can be circumscribedinto the smallest solid member of the topology, whilethe minimum size of the void phase is defined by theradius r min . Void of the largest circle that can be circum-scribed into the smallest cavity of the topology. It hasbeen shown that the minimum length scale, r min . Solid and r min . Void , are defined by the filter radius r fil andthe projection thresholds η ero , η int and η dil . Therefore,if specific values are to be imposed for the minimumlength scale, these must be implicitly imposed throughthe 4 parameters that define the density filter and theHeaviside projection.The dependence of the minimum length scale on theprojection and filtering parameters hurdles a generalimplementation of the robust formulation, specially incommercial codes, as there is no prompt approach tofind the relationship between the involved parametersand the desired length scale. Encouraged by this short-coming, in this section we present and extend the appli-cability of the analytical procedure proposed by Qianand Sigmund (2013), which is described below.To obtain an explicit relation between the minimumlength scale and the parameters that define it, Qianand Sigmund (2013) proposed to apply the three-fieldscheme over a uni-dimensional and continuous designdomain. For instance, to obtain the minimum size of thesolid phase, a design domain ρ ( x ) centered at the coor-dinate x m and containing solid elements on a stretch ofsize h is assumed, as shown in Fig. 2a. The continuous (a) 1D domain (b) Filtered Field(c) Projected Field (d) Filtered field Fig. 2: The three field scheme applied to a one-dimensional design domain in order to obtain the ana-lytical minimum length scale.form of the density filter reads as follows:˜ ρ ( x i ) = (cid:90) x i + r fil x i − r fil ρ ( x ) (cid:18) − | x i − x | r fil (cid:19) dx (cid:90) x i + r fil x i − r fil (cid:18) − | x i − x | r fil (cid:19) dx = (cid:90) x i + r fil x i − r fil ρ ( x ) r fil (cid:18) − | x i − x | r fil (cid:19) dx (5)where x i is any coordinate of x . For example, by choos-ing a filter radius greater than h/
2, the one-dimensionalfiltered field in Fig. 2b is obtained. The filtered field canbe projected using the smoothed Heaviside function ofEq. (3). To simplify the analysis, an infinite steepnessparameter ( β → ∞ ) is assumed. For example, for anHeaviside threshold η i = 0 .
2, the design of Fig. 2c isobtained. The size of the solid phase in the projected de-sign is defined by the length L i , which can be obtainedby finding x i from the equation ˜ ρ ( x i ) = η i . However,to this extent, the size of the solid phase is a function L i ( r fil , η i , h ) that depends on the assumed h value. Todetermine the value of h and to discover the explicitrelation between the size L of the projected field andthe filter and projection parameters, it is necessary torefer to the foundation of the robust formulation.An erosion projection removes solid material fromthe surface of the reference design. This operation en-larges the cavities and thins the structural memberspresent in the intermediate reference design (Sigmund, ote on the minimum length scale and its defining parameters 5 L ( η ero ) ≈ h by solving the following equation:˜ ρ ( x ero ) = (cid:90) h/ − h/ r fil (cid:18) − | x | r fil (cid:19) dx = η ero (6)where the reference system is placed at x ero for integra-tion. By solving the expression in Eq. (6), the size h isobtained as: h = 2 r fil (1 − (cid:112) − η ero ) (7)By obtaining the distance h that produces an erodedprojection of infinitesimal size, it is possible to relatethe minimum size to the filter radius for any projectionwhose threshold meets η i < η ero . A particular case isto choose η i = 0 .
5, which is the value chosen by Qianand Sigmund (2013) to define the intermediate design.However, to extend the analytical method we use anarbitrary η i such that η i < η ero .As mentioned above, the size of the solid phase L inthe projected field can be obtained by equating ˜ ρ ( x i )to η i . Since the condition of robustness is imposed on h (Eq. 7), the length corresponds to the minimum size,i.e. L = 2 r min . Solid . The expression that allows to relatethe minimum size to the filter radius is as follows:˜ ρ ( x i ) = (cid:90) x i + r fil x i − r fil ρ ( x ) r fil (cid:18) − | x i − x | r fil (cid:19) dx = η i (8)For the integration, it is convenient to place the ori-gin of the reference system at x i . To this end, 4 sit-uations must be considered to define the integrationlimits, which are summarized in Fig. 3. For example,if the length of the projected field L is greater than h and the filter radius is greater than ( L + h ) / ρ ( x i ) = (cid:90) L − L − h L − h r fil (cid:18) − | x | r fil (cid:19) dx = η i (9)Solving the integral of Eq. (9) leads to the followingexpression:˜ ρ ( x i ) = hr fil (1 − L r fil ) = η i (10)Finally, replacing Eq. (7) in (10):2 r min . Solid r fil = 2 − η i − √ − η ero (11) (a) L ≥ h and r fil ≥ L + h . (b) L < h and r fil ≥ L + h .(c) L ≥ h and r fil < L + h . (d) L < h and r fil < L + h . Fig. 3: Situations to be considered when defining thelimits of integration.Eq. (11) explicitly relates the minimum size of thesolid phase ( r min . Solid ), of a projected field defined by η i , with the filter radius ( r fil ) and the erosion threshold( η ero ). However, Eq. (11) is only valid for L > h and r fil ≥ ( L + h ) /
2. For implementation purposes, it ismore convenient to express the range of application interms of the projection thresholds. To this end, h canbe replaced from Eq. (7) and L from (11), which leadsto conditions depending only on η i and η ero , as follows: L > h = ⇒ η i < η ero − √ − η ero r fil ≥ ( L + h ) / ⇒ η i ≥ − √ − η ero − η ero (12)By repeating the procedure from Eq. (9) to (12) forthe 4 integration conditions shown in Fig. 3, a set ofequations is obtained which relate the filter and pro-jection parameters with the minimum size for any pro-jection threshold η i , provided that η i < η ero . The setof equations are summarized in the four first rows ofTable 1.To obtain the relationships that define the minimumsize of the void phase, the same procedure must be usedas for the solid phase, however, now starting from a one-dimensional design domain containing a cavity of size h .To avoid overextending the document with redundantinformation, this section is limited to presenting thefinal equations that define the minimum size of the voidphase. The expressions are summarized in the last 4rows of Table 1. Denis Trillet et al.
Table 1: The explicit relationship between the minimum length scale and the filter and projection parameters.
Phase Conditions Minimum size S o li d η i ≤ η ero − η i ≥ . r min . Solid r fil = 2 (cid:112) − η i − (cid:112) − η ero η i ≥ η ero − √ − η ero and η i > η ero − r min . Solid r fil = 2 √ η ero − η i η i < − √ − η ero − η ero and η i < . r min . Solid r fil = 4 − (cid:112) − η ero − (cid:112) η i η i ≥ − √ − η ero − η ero and η i < η ero − √ − η ero r min . Solid r fil = 2 − η i − √ − η ero V o i d η i ≥ η dil and η i ≤ . r min . Void r fil = 2 (cid:112) η i − √ η dil η i ≤ η d + 1 − √ η dil and η i < η dil r min . Void r fil = 2 √ η i − η dil η i ≥ √ η dil − η dil − η i > . r min . Void r fil = 4 − √ η dil − (cid:112) − η i η i < √ η dil − η dil − η i > η dil + 1 − √ η dil r min . Void r fil = 2 − − η i − √ η dil Having delivered the set of equations that expandthe scope of the method proposed by Qian and Sigmund(2013), the following section presents a methodology touse these equations.
In structural design, the minimum length scale con-trol is usually desired because of design requirementsor manufacturing limitations, hence in most cases, theminimum size of the solid and void phases are knownvalues established for the intermediate design. There-fore, for the set of equations presented in Table 1, theradii r intmin . Solid and r intmin . Void are assumed user-definedinput values. In this case, the projection threshold η i corresponds to the projection threshold η int ,hence the desired length scale for the intermediate de-sign is a function of the projection thresholds and ofthe size of the filter, namely, r intmin . Solid ( η int , η ero , r fil ) and r intmin . Void ( η int , η dil , r fil ). Given the number of unknowns( r fil , η ero , η int , η dil ), the system of equations in Table 1becomes indeterminate and the desired minimum lengthscale can be imposed through multiple combinations ofparameters. Nevertheless, such freedom of parametersselection can be reduced by considering the followingthree recommendations.Firstly, a number of advantages have been observedwhen defining the intermediate design with a threshold η int = 0 .
5. For instance, the projection features loweramounts of intermediate densities compared to those projections that use a threshold other than 0.5 (Xuet al, 2010; Wang et al, 2011; da Silva et al, 2019), anda threshold η int set to 0.5 provides the same size ranges(0.5) for the erosion and dilatation thresholds, whichis convenient for reducing rounding errors, since smalldifferences in projection thresholds could be insensitiveto the minimum size when using a coarse discretiza-tion of the design domain (Qian and Sigmund, 2013).Secondly, for a particular combination of thresholds,the filter size ( r fil ) can become considerably larger thanthe desired minimum length scale, which could signif-icantly increase computational requirements (Lazarovand Sigmund, 2011). This can be seen in Figs. 4 and6b. These figures show graphs that relate the filter size( r fil ) to the minimum size of solid or void phase andto the erosion or dilation threshold. These graphs showthat the closer η ero and η dil are to η int , the larger thefilter radius, which inevitably increases computationalrequirements. Under this observation, we recommendchoosing η ero ≥ .
75 and η dil ≤ .
25, thus it is ensuredthat r fil ≤ r intmin . Solid and r fil ≤ r intmin . Void . Thirdly, ero-sion and dilation thresholds too distant or too close tothe intermediate threshold increases oscillations of de-sign variables during the optimization process. In gen-eral, a good compromise is to choose 0 . ≤ η dil ≤ . . ≤ η ero ≤ . η int is set to 0.5, while thesecond and third observations limit the range of theerosion and dilation thresholds. Thus, for a user-defined ote on the minimum length scale and its defining parameters 7 Fig. 4: Graphical relationship between the minimumsize of the void phase, the filter radius and the dilationprojection.minimum length scale, it is possible to develop an algo-rithm that solves the system of equations consideringthe three observations. Here we propose an algorithmbased on graphic relationships, so that the reader caneasily find the desired parameters without the need toresort to a computational algorithm. Nonetheless, wealso provide as supplementary material a code writtenin MATLAB named
SizeSolution.m that performs theprocedure described below. It is important to note thatthe above observations are based on numerical testsconsidered for specific optimization problems formu-lated in the density approach. Therefore, it is possi-ble that under other topology optimization approachesor formulations the above observations are no longervalid. However, the proposed procedure can be appliedfor any other value of η int , or any other combination ofparameters that the user may consider convenient.Given that the ranges of application of the equationsin Table 1 are defined as a function of the projectionthresholds ( η i , η ero and η dil ), it is rather simple to con-struct graphs with respect to them. For instance, asshown in Fig. 5, to find the range of application of theequations defining the minimum size of the solid phase,it is simply necessary to know in which region of Fig. 5the projection η i falls.The first graph proposed in this work is shown inFig. 6a and gathers 4 parameters, the minimum lengthscale ( r intmin . Solid and r intmin . Void ) and the projection thresh-olds ( η ero and η dil ), provided that η int = 0 .
5. In thisgraph, the user can easily find the set of erosion and di-lation threshold that leads to the desired length scale.Then, the user can access the graph in Fig. 6b to obtainthe filter radius. For example, for the following mini-mum length scale, r intmin . Solid = 3 elements and r intmin . Solid = Fig. 5: Graphic representation of the applicability of theset of equations in Table 1 proposed for the solid phase.The denomination of the zone corresponds to the rownumber of Table 1.3 elements, the graph in Fig. 6a is accessed with a valueof 1.0 for the ordinate. According to the aforementionedobservations, the combination of thresholds [ η ero , η dil ]that can be chosen among others are:[ η ero , η dil ] = [0 . , . , (13a)[ η ero , η dil ] = [0 . , . , (13b)[ η ero , η dil ] = [0 . , . , (13c)[ η ero , η dil ] = [0 . , . . (13d)Arbitrarily, [0.75, 0.25] is selected, and from thegraph in Fig. 6b, it is obtained that r fil = 2 r intmin . Solid = 6elements.It should be noted that the graph in Fig. 6a is con-structed considering a resolution of 0.05 in the projec-tion thresholds. This is due to the fact that the dis-cretization of the filter radius in topology optimizationis generally coarse, and decimal numbers smaller than0.05 in the threshold value have usually a negligibleeffect on the minimum length scale of the optimizeddesign. However, if a better resolution is required, theuser can resort to the attached code
SizeSolution.m .As previously mentioned, the erosion and dilationdistances could be required, for instance, when imple-menting maximum size restrictions (Fern´andez et al,2020). These distances can be easily obtained from theequations of Table 1. As shown in Fig. 2d, the dilationdistance, denoted by t dil , is the offset distance betweenthe intermediate and dilated designs. Analogously, theerosion distance, denoted by t ero , is the offset distance Denis Trillet et al. (a) (b)(c) (d)
Fig. 6: Graphical relationships between the minimum length scale and η ero , η int , η dil , and r fil . The minimum lengthscale is defined for the intermediate design. The graphs are built using both, the analytical and the numericalmethod. Graphs (a) and (b) are designed to obtain the projection thresholds and the filter radius, respectively.Graphs (c) and (d) are designed to obtain the dilation and erosion distances, respectively. ote on the minimum length scale and its defining parameters 9 between the intermediate and eroded designs. There-fore: t dil = r dilmin . Solid ( η dil , η ero , r fil ) − r intmin . Solid ( η int , η ero , r fil ) t ero = r eromin . Void ( η ero , η dil , r fil ) − r intmin . Void ( η int , η dil , r fil )(14)where the minimum size of the solid phase in the dilateddesign ( r dilmin . Solid ) is obtained by choosing η i = η dil inthe equations of Table 1. Analogously, the minimumsize of the void phase in the eroded design ( r eromin . Void ) isobtained by choosing η i = η ero . In this way, analyticalexpressions can be obtained for the dilation and erosiondistances.It is noted from Eq. (14) that the erosion and di-lation depend on the 3 projection thresholds and onthe size of the filter, namely, t ero ( r fil , η ero , η int , η dil ) and t dil ( r fil , η ero , η int , η dil ). Therefore, if for some particularreason the user needs to impose the minimum lengthscale, and the erosion and dilation distances, then thefollowing system of equations must be solved: r intmin . Solid = r intmin . Solid ( r fil , η int , η ero ) r intmin . Void = r intmin . Solid ( r fil , η int , η dil ) t ero = r intmin . Solid ( r fil , η ero , η int , η dil ) t dil = r intmin . Solid ( r fil , η ero , η int , η dil ) (15)The system in Eq. (15) is determined and thereis only one combination of parameters [ r fil , η ero , η int , η dil ] that leads to the desired length scale [ r intmin . Solid , r intmin . Void , t ero , t dil ]. This can be illustrated with thefollowing example. For each combination of thresholds[ η ero , η dil ] given in Eq. (13), the erosion and dilationdistances [ t ero , t dil ] are provided, as follows:[ η ero , η dil ] = [0 . , . , [ t ero , t dil ] = [1 . , .
76] (16a)[ η ero , η dil ] = [0 . , . , [ t ero , t dil ] = [1 . , .
99] (16b)[ η ero , η dil ] = [0 . , . , [ t ero , t dil ] = [2 . , .
21] (16c)[ η ero , η dil ] = [0 . , . , [ t ero , t dil ] = [2 . , .
43] (16d)It is recalled that all combinations of thresholds[ η ero , η dil ] in Eq. (16) lead to the minimum length scale r intmin . Solid = r intmin . Void = 3 elements. However, they allresult in different erosion and dilation distances.To the authors’ knowledge, the need to impose spe-cific values for the erosion and dilation distances has notyet been claimed, so this manuscript is limited to pro-viding these values for a given combination of param-eters. The proposed graphs are shown in Fig. 6c and6d, which depend on the erosion and dilation thresh-olds and on the corresponding minimum size. For in-stance, considering [ η ero , η dil ] = [0.75, 0.25], the erosion and dilation distances are t ero = 0 . × t dil = 0 . × The analytical procedure developed in this article al-lows to quickly obtain the set of filtering and projectionparameters that imposes the desired minimum lengthscale. However, in practice, the length scale of the opti-mized design often differs from the desired values (Qianand Sigmund, 2013). This is mainly due to the fact thatthe analytical method assumes (i) a continuous domain,(ii) a perfect Heaviside projection ( β → ∞ ), and (iii)that the eroded and dilated fields project an infinitesi-mal minimum size in the solid and void phases, respec-tively. The assumptions (i), (ii), and (iii) are not metin topology optimization mainly because the design do-mains are discretized into finite elements and becausethe Heaviside projection is smoothed. To assess the er-ror introduced by these assumptions, in the followingwe compare the analytical method with the numericalmethod proposed by Wang et al (2011), which considersa discrete domain and a smoothed Heaviside function.The procedure for obtaining the minimum size us-ing the numerical method is analogous to the analyticalmethod but now using a one-dimensional domain dis-cretized into N elements. For example, to obtain theminimum size of the solid phase the three-field schemeis applied to a one-dimensional design domain ρ con-taining solid elements in a length h . Then, the size ofthe eroded, intermediate and dilated fields are mea-sured, and the length h is adjusted so that the resultingeroded field has an infinitesimal size (1 solid element).This process is repeated several times for different val-ues of η ero , and the resulting minimum size is normal-ized with respect to the size of the chosen filter. For im-plementation details regarding the numerical method,the reader is referred to the works of Wang et al (2011)and Fern´andez et al (2020), and to the attached MAT-LAB code named NumericalSolution.m .To validate the analytical and numerical methods,the latter is implemented using a Heaviside function at β = 500, a filter radius of 1000 elements and a designdomain discretized into 10 thousand elements, so thenumerical method can be considered continuous andunder similar assumptions than the analytical one. Therelationships obtained with the numerical method canbe seen in the graphs of Fig. 6. The agreement of re-sults from both methods allows us to validate the setof equations provided in Table 1.The three sources of error that affect the analyticalmethod are discussed below.0 Denis Trillet et al.
43] (16d)It is recalled that all combinations of thresholds[ η ero , η dil ] in Eq. (16) lead to the minimum length scale r intmin . Solid = r intmin . Void = 3 elements. However, they allresult in different erosion and dilation distances.To the authors’ knowledge, the need to impose spe-cific values for the erosion and dilation distances has notyet been claimed, so this manuscript is limited to pro-viding these values for a given combination of param-eters. The proposed graphs are shown in Fig. 6c and6d, which depend on the erosion and dilation thresh-olds and on the corresponding minimum size. For in-stance, considering [ η ero , η dil ] = [0.75, 0.25], the erosion and dilation distances are t ero = 0 . × t dil = 0 . × The analytical procedure developed in this article al-lows to quickly obtain the set of filtering and projectionparameters that imposes the desired minimum lengthscale. However, in practice, the length scale of the opti-mized design often differs from the desired values (Qianand Sigmund, 2013). This is mainly due to the fact thatthe analytical method assumes (i) a continuous domain,(ii) a perfect Heaviside projection ( β → ∞ ), and (iii)that the eroded and dilated fields project an infinitesi-mal minimum size in the solid and void phases, respec-tively. The assumptions (i), (ii), and (iii) are not metin topology optimization mainly because the design do-mains are discretized into finite elements and becausethe Heaviside projection is smoothed. To assess the er-ror introduced by these assumptions, in the followingwe compare the analytical method with the numericalmethod proposed by Wang et al (2011), which considersa discrete domain and a smoothed Heaviside function.The procedure for obtaining the minimum size us-ing the numerical method is analogous to the analyticalmethod but now using a one-dimensional domain dis-cretized into N elements. For example, to obtain theminimum size of the solid phase the three-field schemeis applied to a one-dimensional design domain ρ con-taining solid elements in a length h . Then, the size ofthe eroded, intermediate and dilated fields are mea-sured, and the length h is adjusted so that the resultingeroded field has an infinitesimal size (1 solid element).This process is repeated several times for different val-ues of η ero , and the resulting minimum size is normal-ized with respect to the size of the chosen filter. For im-plementation details regarding the numerical method,the reader is referred to the works of Wang et al (2011)and Fern´andez et al (2020), and to the attached MAT-LAB code named NumericalSolution.m .To validate the analytical and numerical methods,the latter is implemented using a Heaviside function at β = 500, a filter radius of 1000 elements and a designdomain discretized into 10 thousand elements, so thenumerical method can be considered continuous andunder similar assumptions than the analytical one. Therelationships obtained with the numerical method canbe seen in the graphs of Fig. 6. The agreement of re-sults from both methods allows us to validate the setof equations provided in Table 1.The three sources of error that affect the analyticalmethod are discussed below.0 Denis Trillet et al. r fil , r intmin . Solid , r dilmin . Solid , r intmin . Void and r eromin . Void ) are definedby a discrete number of elements. In this case, therounding error is ± δr , corresponds to 1 ele-ment in the radius, i.e. 2 elements in the diameter. Thatis, the vertical offset of the analytical curve due to therounding error is equal to 2 δr = 2 /
10 and to 2 δr = 2 / δr ) for a given analytical threshold( η ero or η dil ), or the margin of error of the analyticalthreshold ( δη ) for a given minimum size, as shown inFig. 7a.If the user of the robust approach needs to imposeprecisely the minimum size, he can choose the projec-tion thresholds ( η ero and η dil ) that lead to a larger filterradius, or he can use the numerical method (Wang et al,2011) using a discretization representative of the designdomain to be optimized.5.2 Perfect Heaviside projectionThe analytical method developed in this work assumesa perfect Heaviside function, which results in a pro-jected field ¯ ρ containing discrete densities (0 and 1). (a) N = 100, r fil = 10.(b) N = 200, r fil = 20 Fig. 7: Relationship between the minimum size of thesolid phase, the eroded threshold, and the filter radius.The dotted line represents the analytical relationshipincluding the rounding error.However, in topology optimization, the Heaviside func-tion is smoothed resulting in regions of intermediatedensities that lie on the surface of the optimized struc-ture. For this reason, in practice, once the optimizedsolution ¯ ρ int is obtained, a cut-off value ε is usuallyadded to the projected densities in order to get the op-timized structure. Therefore, the final design intendedfor manufacturing is obtained by projecting the design ¯ ρ int as follows:¯ ρ εi = (cid:40) , if ¯ ρ int( i ) ≥ ε , otherwise (17)Since the post-processed design ¯ ρ ε is the one in-tended for manufacturing, the analytical expressionsrelating the minimum length scale and the parametersthat define it must be elaborated for the field ¯ ρ ε . Aswill be explained hereafter, the equations provided inTable 1 can be easily adapted to consider the cut-offvalue ε . ote on the minimum length scale and its defining parameters 11 Fig. 8: Illustration of a transition region where pro-jected densities reach intermediate values.Consider Fig. 8 illustrating a transition zone be-tween the solid structure and the void region. In blueit is shown the filtered field ˜ ρ and in red the projectedfield ¯ ρ obtained with β = 32. In this illustration, thedensity cut-off is defined as ε = 0 .
95. A perfect Heav-iside projection ( β → ∞ ) would define the size of thesolid zone at the coordinate x i . However, the smoothedHeaviside projection defines the solid/void transition inthe coordinate x ε . The idea then is to find the projec-tion threshold η ε for which a perfect Heaviside projec-tion produces the same solid/void transition coordinatethan the one that is obtained with a smoothed Heavi-side projection with a threshold η i . Thus, the value ofthe filtered density present at the coordinate x ε mustbe found. To this end, from Fig. 8 it is observed that: ε = tanh ( βη i ) tanh ( β ( η ε − η i ))tanh ( βη i ) tanh ( β (1 − η i )) (18)Assuming β >
10 at the end of the optimizationprocess, which is a common practice when dealing withHeaviside projection, Eq. (18) can be simplified to:tanh ( β ( η ε − η i )) = 2 ε − η ε can be ob-tained: η ε = η i + 1 β atanh(2 ε −
1) (20)We recall that the filtered density η ε represents thevalue for which a perfect Heaviside function with thresh-old η ε and a smoothed Heaviside function with thresh-old η i result in the same size for the solid and voidphases. Therefore, to take into account the intermedi-ate densities resulting from a smoothed Heaviside func-tion, the shifting term atanh(2 ε − /β must be addedto the thresholds involved in the analytical equations(Table 1). Note that if β → ∞ or ε = 0 .
50, then theshifting term is zero and η ε = η i . (a)(b) Fig. 9: (a) Minimum length scale obtained with the an-alytical method assuming a perfect Heaviside projec-tion ( β → ∞ ), and with the numerical method usinga smoothed Heaviside projection ( β = 32) and threedifferent cut-off values ε . (b) The analytical curves ob-tained by adding the shifting term atanh(2 ε − /β tothe thresholds.To assess the error introduced by a smoothed Heav-iside projection, the minimum length scale obtainedwith the analytical and the numerical methods are com-pared. The analytical method is applied with β → ∞ ,so it does not consider the shifting term on the thresh-olds. For the numerical method we use β = 30, a filterradius equal to 200 and a one-dimensional domain dis-cretized into 2000 elements in order to avoid roundingerrors and isolate the effect of using a smoothed Heavi-side projection. In the numerical method, three cut-offvalues are used, ε = 0 . ε = 0 .
50 and ε = 0 .
99. Theresults are shown in Fig. 9a.The results show that the influence of using a smoothedHeaviside function is observed when ε (cid:54) = 0 . η ero > .
75. For this reason, we recommend to use adensity cut-off value equal to 0.5, as this avoids theneed to add the shifting term to the projection thresh-olds. Nevertheless, the user can easily adjust the ana-lytical curves for ε (cid:54) = 0 .
5, since it is only required to addatanh(2 ε − /β to the projection thresholds, as shownin the corrected curves in Fig. 9b.It is worth mentioning that probably the effect ofusing a smoothed Heaviside function is only seen for η ero ≥ .
75 because the projected field becomes morediscrete as the projection threshold is closer to 0.5 (Wanget al, 2011).5.3 Infinitesimal sizeIt is recalled that the analytical and numerical meth-ods are developed under the condition of robustnesswhich ensures that the manufactured design will fea-ture a good performance even if the blueprint designis uniformly thinned or thickened during the manufac-turing process. This condition dictates that the struc-tural members/cavities must be projected by at leastone solid/void element in the case of an eroded/dilateddesign.To obtain the minimum size, the analytical methodassumes that the size of the structural members/cavi-ties is infinitesimal ( ≈
0) in the eroded/dilated design,but in practice this does not occur. In topology op-timization, it is observed that the smallest structuralmembers/cavities in the eroded/dilated design are com-posed of one or a few elements, generally described withintermediate densities. Thus, in practice, the minimumsize of the solid phase in the eroded design ( r eromin . Solid )results in a discrete number of elements. This can beseen in Fig. 10, which is build using the numericalmethod. There, a filtered field and its eroded projec-tion under the robust condition are shown. It can beseen that the minimum size of the eroded design isnot infinitesimal, and therefore the minimum size ofthe solid phase in the intermediate and dilated designwould be larger than the value predicted by the analyt-ical method.The minimum size of the solid and void phases inthe eroded and dilated designs ( r eromin . Solid and r dilmin . Void )are defined by the size of the elements that discretizethe design space. In addition, the amount of elementswith intermediate densities defining the minimum sizedepends on the steepness of the smoothed Heavisidefunction ( β ), therefore, the error introduced by assum-ing an infinitesimal size in the robust condition is re-lated to the two sources of error mentioned previously.To isolate the effect of the infinitesimal size in the ro-bust condition and illustrate the error introduced by Fig. 10: Illustration in which an infinitesimal size is notreached for the eroded projection.the assumption of infinitesimal size, we make use ofthe numerical method. The numerical method is im-plemented in a discretized domain containing a largenumber of elements (10 elements) and using a largefilter size (10 elements) to simulate a continuous do-main. The steepness parameter of the smoothed Heav-iside projection is set as β = 512. Thus, the size anddensity of the elements are excluded from the analy-sis. The effect of not achieving an infinitesimal size inthe condition of robustness is intentionally imposed inthe numerical code. For this, the robustness conditionis considered satisfied if r eromin . Solid = α r fil . Consideringthat r eromin . Solid represents one finite element in topologyoptimization, and that representative values for r fil arebetween 2 and 10 elements, it is reasonable to considervalues for α ranging from 0.1 to 0.5. Taking into consid-eration the above, we build the graphical solutions thatrelate the minimum size in the solid phase ( r intmin . Solid ),the filter radius and the erosion threshold. The graphis shown in Fig. 11.The graph shows that the error of not achievingan infinitesimal size in the eroded design produces aminimum size of the solid phase bigger than the valuepredicted by the analytical method. The error relatedto the infinitesimal size is low in comparison to therounding error, so the latter would be the most rele-vant source of error from the analytical method thatassumes a continuous design domain.In summary, this section discussed the scope of theanalytical method (Qian and Sigmund, 2013) by com-paring it with the numerical one (Wang et al, 2011),both developed for a one-dimensional design domain.To this end, different sources of error were examined. Ingeneral, the errors can be controlled by mesh refinementor by correcting the projection thresholds according tothe cut-off ε value. In the following section, the analyt- ote on the minimum length scale and its defining parameters 13 Fig. 11: Effect of not reaching an infinitesimal size in thecondition of robustness r eromin . Solid = 0. This condition isassessed by imposing r eromin . Solid = α r fil .ical method is assessed using 2D-topology optimizationproblems. This section examines the reliability of the analyticalexpressions provided in Table 1. To this end, a set of 2Dtopology optimization problems are solved, from whichthe length scale is measured graphically and comparedwith the imposed values. Then, some designs obtainedwith maximum size constraints are provided to illus-trate the use of the erosion and dilation distances. Fi-nally, this section provides a remark regarding the sim-plified robust formulation, where the intermediate anddilated designs are removed from the objective function.6.1 Minimum length scaleThe minimum length scale is assessed using the heat ex-changer design problem described in Section 2. A set ofresults is generated from this design problem, which dif-fer in the desired minimum length scale and in the dis-cretization used. Specifically, three sets of solutions areobtained by discretizing the design domain into 100 × ×
200 and 400 ×
400 quadrilateral elements.In addition, three different length scales are prescribedfor each discretization, which are reported as the ratiobetween the minimum size of the solid phase and theminimum size of the void phase, i.e. r intmin . Solid /r intmin . Volid .The chosen ratios are 1 /
2, 1 / /
1. The minimumsize in the solid phase is the same in all scenarios andis defined as a physical dimension. In number of finiteelements, the radius that defines the minimum size of Fig. 12: Initial distribution of design variables consid-ered for the thermal compliance minimization problem.the solid phase ( r intmin ) is equal to 1, 2 and 4, for thediscretizations that use 100 , 200 and 400 elements,respectively. It is well known that the initial values ofdesign variables have a huge influence on the resultingtopology when it comes to thermal compliance mini-mization (Yan et al, 2018), hence, in order to facilitatethe comparison of results, we impose a base structureas a starting point, which is shown in Fig. 12. Beforepresenting the results, the procedure to obtain the min-imum length scale from the optimized designs is de-tailed.The minimum size of the solid phase is measuredgraphically by counting the number of finite elementsthat define the size of the thinnest structural branch.Similarly, the minimum size of the void phase is mea-sured by counting the elements in the radius of thelargest circumference that can be inscribed at the re-entrant corners of the design. For example, consider thedesign of Fig. 13c where the minimum length scales r intmin . Solid = min . Void int = 2 elements are imposed. Todetermine the real minimum size of the void phase, thelargest circle in Fig. 13a that falls into the re-entrantcorners of the design is identified. The corners analyzedare those that form a sharp angle between two struc-tural branches. Fig. 13b shows three representative re-entrant corners of the design depicted in Fig. 13c. Fromthere it is observed that the minimum size is given bya circle of radius 2 elements (zone C). Similarly, thelargest region that fits into the thinnest structural mem-bers is determined, as shown in Fig. 13d. There, theminimum size of the solid phase is given by a circle ofradius 1.5 elements (zone D).After describing the test case, the obtained resultsare presented. Table 2 contains the nine results gen-erated in this example (3 length scales × Fig. 13: Illustration of the minimum length scale mea-surement. In (a) the test regions and their radius sizesgiven in number of finite elements. In (b) the minimumsize of the void phase, (c) the optimized heat exchanger,and (d) the minimum size of the solid phase. The im-posed length scales are graphically shown in (c), at theupper left corner of the design.graphically next to each solution. The minimum size ofthe void phase is indicated in blue, while the minimumsize of the solid phase in magenta. The table also reportsthe 3 parameters required to impose the desired mini-mum length scales, i.e. η ero , η dil , and r fil . These param-eters have been obtained using the analytical methodimplemented in the MATLAB code provided with thispaper ( SizeSolution.m ). It is recalled that all the ex-amples assume η int = 0 . in the graphs of Fig. 14, while the valuesmeasured from the optimized designs are labeled as .The error bars in Fig. 14a illustrate the rounding errorpresent in the analytical method (see Section 5.1 fordefinition). Rounding error lines have not been plottedin Fig. 14b because the length scale ratios have differentrounding errors.As a general observation, we can point out fromFig. 14 that mesh refinement reduces the error betweenthe measured minimum size and the desired minimumsize (which is imposed through the set of parametersprovided by the analytical method). This is consistentwith the observations made in the previous section,where a continuous and uni-dimensional design domainwas used. In addition, the predicted error for the ana-lytical method (the error bars) are also consistent withthe measured sizes, which validates the scope and limi-tations of the analytical equations discussed in Section5. Regarding the graph in Fig. 14a, it can be men-tioned that in the nine 2D designs, the error betweenthe measured and desired minimum sizes is half a finiteelement. This error is relatively large in the coarse dis-cretization (50% error), and small in the fine discretiza-tion (12.5% error). However, despite the fact that theanalytical method is not exact, it seems to be accurateenough in the examples examined, since in all cases theerror is the same, half a finite element. It is interestingto note that the numerical method that assumes a dis-crete 1D domain, estimates a minimum size that differsby half a finite element with respect to the imposedvalue, which matches the measured error. However, thenumerical method provides a minimum size of the solidphase larger than the measured one. This is probablydue to the fact that the condition of robustness refersto one finite element and not to an infinitesimal size, asdiscussed in Section 5.3.On the other hand, the graph built for the voidphase (Fig. 14b) shows a different pattern. The mea-sured minimum size of the void phase is often equal tothe desired one and even bigger for some designs. Thiscould be explained by the simple fact that the chosentest case does not present a geometric singularity overthe minimum size of the void phase, therefore it is notpossible to guarantee that the smallest re-entrant cor-ner of the design will indeed correspond to the minimumsize imposed by the robust formulation. The represen-tative case might be that where r intmin . Void is chosen twicethe size r intmin . Solid (ratio 1/2). In such a case, the mea-sured error corresponds to half a finite element smallerthan the desired value, that is, the same result as forthe solid phase. ote on the minimum length scale and its defining parameters 15
Table 2: Optimized designs for the heat exchanger problem. The imposed volume constraint is V ∗ min = 20%. Mesh r intmin . Solid /r intmin . Void / / / using [ η ero , η dil ] = [0.65, 0.05] using [ η ero , η dil ] = [0.70, 0.30] using [ η ero , η dil ] = [0.80, 0.42]100 × r intmin . Solid = 1, r fil = 2 . r intmin . Solid = 1, r fil = 2 . r intmin . Solid = 1, r fil = 1 . × r intmin . Solid = 2, r fil = 5 . r intmin . Solid = 2, r fil = 4 . r intmin . Solid = 2, r fil = 3 . × r intmin . Solid = 4, r fil = 10 . r intmin . Solid = 4, r fil = 9 . r intmin . Solid = 4, r fil = 7 . ¯ ρ ero ) , c( ¯ ρ int ) , c( ¯ ρ dil ))s.t.: v (cid:124) ¯ ρ dil ≤ V ∗ dil ( V ∗ int )G ms ( ¯ ρ dil ) ≤ ≤ ρ i ≤ , i = 1 , ..., N , (21)where G ms is the maximum size constraint defined ex-actly as in (Fern´andez et al, 2020). Therefore, G ms rep-resents a p-mean aggregation function that gathers lo- Fig. 14: The graphs summarize the measured minimumsize of (a) the solid phase and (b) the void phase fromthe designs in Table 2. For each discretization, two setof results are provided, labeled as 2 D and 1 D . The 2 D represents the graphical measurement normalized withrespect to the imposed value. The 1 D represents theminimum size obtained from the numerical method nor-malized with respect to the imposed value.cal maximum size restrictions. As in the cited work, the p aggregation exponent is set at 100.Regarding the optimization parameters, these rep-resent the implementation of Wang et al (2014), i.e.the Heaviside parameter β is initialized at 1 and is in-creased by 1 every 20 iterations until a value β = 16 isreached. Then, 20 more iterations are carried out witha β = 32. The SIMP penalty parameter is set at 3. Wehave found that this setting of parameters works wellfor introducing maximum size restrictions into the non-linear force inverter formulated under the robust designapproach. To impose minimum and maximum length scales,the maximum size constraint G ms is applied on the di-lated design (Fern´andez et al, 2020). To do so, the re-gions where the local maximum size constraints are ap-plied must be scaled according to the dilation distance.For instance, if the desired maximum size in the inter-mediate design is defined by a circle of radius r intmax . Solid ,the maximum size constraint should be formulated forthe dilated design using a circle of radius: r dilmax . Solid = r intmax . Solid + t dil (22)Equation (22) explains the need of knowing the dila-tion distance when imposing maximum size restrictionsin the robust formulation. As mentioned previously, thisinformation can be easily obtained from the graphs gen-erated by the analytical equations. For example, in thefollowing, the half force inverter depicted in Fig. 1 issolved. The design domain is discretized into 200 × r intmin . Solid = 2 elements, while themaximum size is set as r intmax . Solid = 3 elements. Twodifferent values for the minimum size of the void phaseare chosen, r intmin . Void = 2 and r intmin . Void = 3 elements.The filter and projection parameters used to imposethe desired length scale are reported in Table 3. In allcases, the volume constraint is set to 25%.To obtain the dilation distance t dil , the projectionthresholds that define the minimum length scale mustbe defined. These values can be graphically obtainedfrom Fig. 6, or from the attached MATLAB code( SizeSolution.m ). For example, using this code, thesets of parameters in Table 3 are obtained. The dilationdistance is also reported there. For the minimum sizes r intmin . Void = 2 and r intmin . Void = 3 elements, the dilationdistances are t dil = 1 and t dil = 2 elements, respec-tively (the numbers have been rounded to the nearestinteger).The results are shown in Figs. 15 and 16. Each figurereports 2 optimized designs, which are obtained withand without maximum size restrictions. To facilitatethe visual comparison between results, the designs areplaced with respect to the symmetry axis and are re-ported in their deformed configuration. The imposedlength scale is also reported graphically next to eachdesign. As in the previous examples, the blue and ma-genta circles represent the minimum size of the voidTable 3: Parameters to impose r intmin . Solid = 2 on thedesign with β = 32. r intmin . Void r fil η ero η dil t dil t ero .
47 0 .
70 0 .
30 1 .
03 1 .
033 4 .
47 0 .
70 0 .
11 2 .
41 1 . Fig. 15: Nonlinear force inverter without (upper half)and with (lower half) maximum size constraints. Theminimum length scales are r intmin . Solid = r intmin . Void = 2elements. The maximum size is r intmax . Solid = 3 elements.Fig. 16: Nonlinear force inverter without (upper half)and with (lower half) maximum size constraints. Theminimum length scales are r intmin . Solid = 2 elements and r intmin . Void = 3 elements. The prescribed maximum sizeis r intmax . Solid = 3 elements.and solid phase, respectively. The black circle repre-sents the maximum size desired for the intermediatedesign (which is imposed through the dilated design) . The maximum size is actually imposed using an annularregion (Fern´andez et al, 2020), but for illustrative purposes,a circular region is drawn.
The results are consistent with the imposed lengthscales. In the 4 designs reported in Figs. 15 and 16,the minimum size of the solid and void phases are metwith a half-finite element of error (which agrees withthe analytical rounding error). Regarding the maximumsize, this also matches the imposed value. However, dueto the inherent drawbacks of the aggregation functionand the strong non-linearity of the force inverter de-sign problem, there are local regions where the imposedmaximum size restriction is not met, such as the hori-zontal bar in the design of Fig. 16.This force inverter test case shows the usefulness ofthe proposed equations and provided codes, since theyallow to quickly obtain the filter and projection param-eters, and the dilation distance that allow to impose thedesired length scales.6.3 Remark on the simplified robust formulationBefore concluding this manuscript we devote a discus-sion concerning the robust formulation. It is widelyknown that the robust formulation of the optimizationproblem can be simplified when the objective functionis monotonically dependent on the volume fraction. Themost widespread case in the literature is the complianceminimization problem, where the intermediate and di-lated fields have less compliance than the eroded de-sign and therefore can be removed from the objectivefunction without compromising the robustness of theformulation. In this case, the objective function is eval-uated only for the eroded design, with the intermedi-ate and dilated fields remaining solely for formulatingdesign constraints. When the volume restriction is theonly constraint included in the optimization problem,the constraint can be formulated using the intermedi-ate design, assuming that this design field is the oneintended for manufacturing. However, it has been men-tioned in the literature that evaluating the volume re-striction in the dilated design promotes convergence tobetter optimums since numerical instabilities are pre-vented. In this section we add another reason that hasnot been mentioned so far (to the best of the authors’knowledge). If the dilated design is not included in thevolume restriction, then it is not possible to ensure thecontrol over the minimum size of the void phase. Toexplain this statement, consider the following two ro-bust topology optimization formulations for the ther-mal compliance minimization problem: (cid:122) (cid:125)(cid:124) (cid:123) (cid:122) (cid:125)(cid:124) (cid:123) min ρ c ( ¯ ρ ero ) min ρ c ( ¯ ρ ero )s . t . : v (cid:124) ¯ ρ int ≤ V ∗ int , s . t . : v (cid:124) ¯ ρ dil ≤ V ∗ dil ( V ∗ int ) (23)0 ≤ ρ i ≤ ≤ ρ i ≤ V ∗ dil is scaled accord-ing to V ∗ int . From Eq. (23) it is clear that P.I( ¯ ρ ero , ¯ ρ int )and P.II( ¯ ρ ero , ¯ ρ int , ¯ ρ dil ), or alternatively, P.I( η ero , η int )and P.II( η ero , η int , η dil ). Therefore, P.I is not influencedby the dilation threshold nor the dilated design.We recall that the condition of robustness imposedfor the void phase involves the dilated design, which hasto project a cavity with at least one void element to bepresent in the 3 fields that constitute the robust formu-lation. The influence of the dilated design on the min-imum size of the void phase can be seen graphically inFig. 6a by considering a fixed erosion threshold and dif-ferent dilation thresholds. For example, for the set [ η dil , η ero ] consider the points [0 . , .
60] and [0 . , . .
3. These sets of parameters are resumed inTable 4. The results are summarized in Fig. 17.Clearly, the results obtained using P.II show a lengthscale consistent with the expected one, as the design inFig. 17a features bigger reentrant corners than the de-sign in Fig. 17b. However, for any value of η dil , theresult from P.I is always the same and corresponds tothat shown in Fig. 17d. The interpretation that can bemade of the P.I problem is that it uses a dilation thresh-old equal to the intermediate one, i.e. η dil = η int = 0 . η dil , η ero ] equal to [0 . , . The robust topology optimization formulation based onuniform manufacturing errors has gained increasing ac-ceptance in the topology optimization community. This Table 4: Parameters to compare the robust formulationand its simplified form. r intmin . Solid r intmin . Void r fil η e η d .
65 0 .
60 0 .
504 4 12 .
65 0 .
60 0 .
404 8 12 .
65 0 .
60 0 . η dil = 0 .
40. (b) P.II, η dil = 0 . η dil = 0 .
50. (d) P.I for any η dil . Fig. 17: Heat exchanger design problem using two dif-ferent variations of the robust formulation. P.I evaluatesthe volume constraint directly in the intermediate de-sign, while P.II does it through the dilated design. Here, η ero = 0 .
60 and r intmin . Solid = 4 elements.is mainly due to its ability to control the minimum sizeof both the solid and void phases, and its potential to becombined with other topology optimization approaches.Despite the increasing popularity of the formulation,no method was yet available to easily obtain the fil-ter and projection parameters that produce the desiredminimum length scales. This need encouraged us to fur-ther develop the analytical method proposed by Qianand Sigmund (2013). The scope and limitations of thismethod were assessed using the numerical method ofWang et al (2011) and a set of 2D design results fromtwo topology optimization problems, the thermal com-pliance minimization problem and the non-linear forceinverter.In addition to providing a fast and effective way toobtain the parameters that impose the desired mini- ote on the minimum length scale and its defining parameters 19 mum length scale, this work shows that to obtain si-multaneous control over the minimum sizes of the solidand void phases, it is necessary to involve the 3 fields(eroded, intermediate and dilated) in the robust topol-ogy optimization problem. For example, for the compli-ance minimization problem subject to a volume restric-tion, it is known that intermediate and dilated designscan be excluded from the objective function, but thevolume restriction has to be applied to the dilated de-sign in order to involve all 3 designs in the formulation.
This manuscript contains two MATLAB codes as sup-plementary material. The first is called
SizeSolution.m and provides a list of filter and projection parametersthat impose user defined minimum length scales. Thesecond is called
NumericalSolution.m and builds thegraphs in Fig. 6 using the numerical method proposedby Wang et al (2011).
Acknowledgements
The authors acknowledge the researchproject FAFIL (
Fabrication Additive par D´epˆot de Fil ), fundedby INTERREG and the European Regional Development Fund(ERDF).
Conflict of interest
On behalf of all authors, the corresponding author statesthat there is no conflict of interest.
References
Andreasen CS, Elingaard MO, Aage N (2020) Level settopology and shape optimization by density methodsusing cut elements with length scale control. Struc-tural and Multidisciplinary Optimization pp 1–23Bendsøe MP, Kikuchi N (1988) Generating optimaltopologies in structural design using a homogeniza-tion method. Computer methods in applied mechan-ics and engineering 71(2):197–224Bendsøe M (1989) Bendsoe, m.p.: Optimal shape designas a material distribution problem. structural opti-mization 1, 193-202. Structural Optimization 1:193–202Bourdin B (2001) Filters in topology optimization. In-ternational journal for numerical methods in engi-neering 50(9):2143–2158Bruns TE, Tortorelli DA (2001) Topology optimizationof non-linear elastic structures and compliant mecha-nisms. Computer methods in applied mechanics andengineering 190(26-27):3443–3459 Chen S, Chen W (2011) A new level-set based approachto shape and topology optimization under geometricuncertainty. Structural and Multidisciplinary Opti-mization 44(1):1–18Christiansen R, Lazarov B, Jensen J, Sigmund O (2015)Creating geometrically robust designs for highly sen-sitive problems using topology optimization: Acous-tic cavity design. Structural and MultidisciplinaryOptimization 52:737–754Clausen A, Andreassen E (2017) On filter boundaryconditions in topology optimization. Structural andMultidisciplinary Optimization 56(5):1147–1155da Silva GA, Beck AT, Sigmund O (2019) Topology op-timization of compliant mechanisms with stress con-straints and manufacturing error robustness. Com-puter Methods in Applied Mechanics and Engineer-ing 354:397 – 421Fern´andez E, Yang Kk, Koppen S, Alarc´on P, BauduinS, Duysinx P (2020) Imposing minimum and maxi-mum member size, minimum cavity size, and mini-mum separation distance between solid members intopology optimization. Computer Methods in Ap-plied Mechanics and Engineering 368:113,157Fern´andez E, Ayas C, Langelaar M, Duysinx P (2021)Topology optimization for large-scale additive man-ufacturing: Generating designs tailored to the depo-sition nozzle size (Under Review.)Kumar P, Fern´andez E (2021) A numerical schemefor filter boundary conditions in topology optimiza-tion on regular and irregular meshes (Under Review,arXiv:2101.01122v1)Lazarov BS, Sigmund O (2011) Filters in topology op-timization based on helmholtz-type differential equa-tions. International Journal for Numerical Methodsin Engineering 86(6):765–781Pedersen C, Allinger P (2006) Industrial Implementa-tion and Applications of Topology Optimization andFuture Needs, vol 137, Springer, pp 229–238Pellens J, Lombaert G, Lazarov B, Schevenels M (2018)Combined length scale and overhang angle control inminimum compliance topology optimization for addi-tive manufacturing. Structural and MultidisciplinaryOptimizationQian X, Sigmund O (2013) Topological design ofelectromechanical actuators with robustness towardover-and under-etching. Computer Methods in Ap-plied Mechanics and Engineering 253:237–251Sigmund O (1997) On the design of compliant mecha-nisms using topology optimization. Journal of Struc-tural Mechanics 25(4):493–524Sigmund O (2009) Manufacturing tolerant topology op-timization. Acta Mechanica Sinica 25(2):227–2390 Denis Trillet et al.