Accepted Manuscript A new multi-material level set topology optimization method with the length scale control capability Jikai Liu, Yongsheng Ma PII: DOI: Reference: S0045-7825(17)30679-5 https://doi.org/10.1016/j.cma.2017.10.011 CMA 11638 To appear in: Comput. Methods Appl. Mech. Engrg. Received date : 17 April 2017 Revised date : 2 October 2017 Accepted date : 5 October 2017 Please cite this article as: J. Liu, Y. Ma, A new multi-material level set topology optimization method with the length scale control capability, Comput. Methods Appl. Mech. Engrg. (2017), https://doi.org/10.1016/j.cma.2017.10.011 This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain. A new multi-material level set topology optimization method with the length scale control capability Jikai Liu1,2,*, Yongsheng Ma2 1 Department of Mechanical Engineering and Materials Science, University of Pittsburgh, Pittsburgh, PA 15261, United States 2 Department of Mechanical Engineering, University of Alberta, Edmonton, Canada *Corresponding author. Email: [email protected] Abstract Multi-material level set topology optimization methods conventionally rely on the overlap of multiple level set functions to realize the multi-material structural representation. However, this representation may produce redundant material phases and the signed distance information is no longer available within the individual material regions. To fix these issues, a new multi-material level set topology optimization method is proposed in this paper, where m level set functions represent m material phases plus the void, and the signed distance information is straightforwardly available in each material phase. To be specific, each level set function corresponds to a material phase and the overlapping areas are filled with an artificial material type, the property of which is weaker than any of the involved material types. This artificial material plays the role of penalizing the overlap areas to vanish, which is different from the close-to-zero property material type for the void. Hence, the optimization process starts with any multilevel set interpolated input and the overlapping areas will gradually vanish. Based on this new method, we have successfully realized the component length scale control on multi-material structures and innovatively, an approach has been proposed to realize the uniform component thickness control without need of pre-specifying the thickness target. Keywords: Multi-material; Level set; Topology optimization; Length scale control 1. Introduction Multi-material structure is a type of composite, and because of the combined multiple material properties, it has the potential to produce creative structural design with multi-functionality. On the other hand, because of the enlarged design space, it is impractical to perform the multi-material structural design through the conventional trial-and-error approach. Therefore, multi-material topology optimization has attracted a great deal of attention. Also, owing to the rapid development of additive manufacturing, physically producing the multi-material component at a reasonable cost is no longer a major issue. So far, multi-material topology optimization has been implemented under the main topology optimization frameworks, including SIMP (Solid Isotropic Material with Penalization) [1], ESO (Evolutionary Structural Optimization) [2], and level set [3,4], and a brief literature survey is presented in the rest of this section. 1 The general multi-material interpolation under the SIMP framework was summarized by Bendsoe and Sigmund in [5]. In [6,7], the authors applied the multi-material SIMP method to design two-material meta-materials for extreme physical properties; and later, they [8] extended the method to design multiphysics actuators with two types of materials. The general way of two-material interpolation is to assign ) for two design variables to material element: one ( ) for the solid material fraction and the other ( the material composition fraction of material 1. For instance, the elastic tensor is interpolated through Eq. (1). 1 (1) and are the elasticity tensors of material type 1 and 2, respectively; is the penalization where term. Hvejsel and Lund [9] generalized the SIMP method to make it applicable to interpolate any number of predefined materials with isotropic or anisotropic properties. Taheri and Suresh [10] recently combined the SIMP method and iso-geometric analysis (IGA) to design multi-material and functionally graded structures. And Luo et al. [11] applied the SIMP method to design the wrinkle-free cable-suspended membrane structures. Lund and Stegmann [12,13] proposed the DMO (discrete material optimization) method to optimize the discrete fiber orientation selection for laminate composites, which can be treated as a new material interpolation scheme; see Eq. (2). 1 1 (2) This method has been widely adopted for laminate composite design. Lund [14] employed the DMO for multi-material shell laminate design to maximize the buckling load factor. Gao and Zhang [15] compared the SIMP and DMO interpolation schemes through studying the mass-constrained multi-material topology optimization problems. Later, several modified multi-orientation interpolations were proposed such as the Shape Function with Penalization (SFP) [16,17] and the Bi-value Coding Parameterization (BCP) [18,19], which demonstrated the advantages of reduced design variables. Level set method is another topology optimization framework, where multi-material problems have been extensively studied. Wang and Wang [20] developed the ‘color’ level set method for multi-material interpolation. The advantage of this method is that, it only requires level set functions to represent 2 material phases. Mei and Wang [21] presented a similar idea to perform multi-material level set topology optimization. Later, the ‘color’ level set method was adopted for different topology optimization problems, including design of compliant mechanisms [22], heterogeneous objects [23], heat conduction structures [24], and stress-constrained structures [25]. In [25,26], the authors demonstrated rigorous shape derivatives for multi-material topology optimization problems. Recently, the continuous material property transition at the material/material interface was considered by Vermaak et al. [27] and the cohesive model was addressed by Liu et al. [28]. Other than the ‘color’ level set method, Wang et al. [29] proposed the Multi-Material Level Set (MMLS)。 This method represents 1 material phases by level set functions through the whether or not logic. And it guarantees no redundant material phases. The MMLS was inherited by Cui et al. [30], who conducted the multi-material level set topology optimization by solving the reaction diffusion equation; 2 and also by Liu and Ma [31] for a novel perforated pipeline design, where the solid and perforated shell elements were treated as two different material types. Noted that, the ‘degeneration’ process demonstrated in [31] is only a special case, and the ‘color’ level set and MMLS are essentially different, where the former employs the table-like indicator to combine the level set functions and the latter uses the ‘whether or not’ logic to remove the redundancy material phases. Wang et al. [32] applied the MMLS to design multi-phase mechanical meta-materials. Other level set related methods including the piecewise constant level set method [33] and the reconciled level set method (RLSM) [34,35] are also effective in deriving multi-material topological design. In addition, level set topology optimization with multiple embedded components belongs to another type of multi-material problem [36]. In addition, multi-material topology optimization has also been performed under other frameworks, e.g., BESO (Bidirectional Evolutionary Structural Optimization) [37,38] and Phase field [39,40]. In summary of the level set-based methods, length scale control has rarely been studied for multi-material problems. Given the reason, overlap of the multiple level set functions causes the signed distance information lost [35], where however, length scale control on level set method relies heavily on the signed distance characteristic [41–43]. Therefore, we propose a new multi-material level set topology optimization method, where m level set functions represent m material phases plus the void, and the signed distance information is straightforwardly available in each material phase. To be specific, each level set function corresponds to a material phase and overlapping of the level set functions is allowed but penalized by filling the overlapping areas with an artificial material type. The property of this artificial material is weaker than any of the formally involved material types. In this way, existence of the overlapping areas means low effectiveness in material use because of the material property weakening. Apart from that, the material consumption is double counted within the overlapping area, which aggravates its ineffectiveness. Hence, the overlapping areas will gradually vanish during the optimization process. Based on this new method, we have successfully realized the component length scale control on multi-material structures and innovatively, an approach has been proposed to realize the uniform component thickness control without need of pre-specifying the thickness target. 2. Multi-material level set interpolation 2.1 Brief introduction to level set method Level set method, proposed by Osher and Sethian (1988), is an effective implicit interface modeling and tracing method, which has demonstrated the success in several engineering fields, including image processing, computational fluid dynamics, and structural optimization, etc. Specifically in structural optimization, the level set function Φ : ⟼ models the design domain in an implicit manner, through Eq. (3). 0, 0, 0, ∈ / ∈ ∈ / 2 3 represents the design domain, where ∈ and ∂Ω represents the material/void interface. (3) ∈ 2 3 is the material domain, 3 The Heaviside and Dirac delta functions are generally adopted to formulate the numerical region and boundary integration, respectively; see Eq. (4) and (5). Φ Φ Φ 1, 0, Φ , Φ Φ Φ 0 0 Φ (4) Φ (5) 1 Then, interior and boundary of the material domain can be represented by: | | 1 (6) 0 (7) 2.2 Multi-material level set interpolation Normally, multiple level set functions are combined to interpolate the material phases, because a single level set function can only distinguish two phases through the positive and negative signs. Noted that Saye and Sethian [44] recently proposed the non-signed distance level set method; however, it has yet been extended to address structural optimization problems. There are two representative ways to realize the multi-level set interpolation: the ‘Color’ Level Set (CLS) [20] and the Multi-Material Level Set (MMLS) [29]. Figure 1. Multi-material level set schemes The CLS has the characteristic that only level set functions are required to represent 2 material |Φ phases. As presented in Fig. 1a, there are two level set domains: Ω 0 1, 2 , four , Φ 1, 2, 3, 4 . Taking | Φ material phases: elastic properties for example, the interpolation is: 1 Φ Φ Φ Φ 1 Φ Φ Φ 1 Φ (8) 4 in which , , , and are the elastic tensors of material type 1, 2, 3, and 4, respectively. Differently, MMLS represents 1 material types by level set functions. As presented in Fig. 1b, the |Φ 1, 2, 3 . And two level set domains: Ω 0 1, 2 form three material phases: the specific interpolation is shown in Eq. (9). This interpolation method effectively eliminates the redundant material phases. 1 Φ Φ Φ Φ 1 ⋃ With either the CLS or MMLS, both the completeness are satisfied. Φ and uniqueness (9) ∩ ∅ 2.3 New multi-material level set interpolation Given the single-material level set topology optimization, the material domain is modeled by the signed distance field, where abundant geometric information is accessible through mathematical calculations, e.g. the length scale evaluation [41,42]. However, overlap of the level set fields causes loss of this information. Performing geometric manipulations such as length scale control is no longer a trivial task. To address this issue, this paper presented a novel multi-material level set interpolation, where the overlapping areas are filled with an artificial weak material type and would be suppressed during the evolution process. Consequently in the optimization result, each material phase is represented by a single level set function and therefore, the signed distance-based geometric information are well maintained in each material phase, which enables the independent length scale control. To be specific, level set functions are employed to represent the new interpolation is demonstrated in Eq. (10). Φ Φ 1 Φ Φ Φ material types plus a void phase. Then, 1 Φ (10) represents the elastic tensor of the artificial weak material type. Note that, this artificial In Eq. (10), weak material type is different from the close-to-zero properties applied to the void phase. The purpose of Φ . adopting this artificial weak material phase is to suppress the overlapping area Φ Given its working principle, the material volume fraction constraints as demonstrated in Eq. (11) are employed as part of the topology optimization problem formulation. Φ Ω (11) Φ Ω and are the upper bounds of the material volume fractions of phase 1 and 2, respectively. Therefore, it is obvious that eliminating the overlapping area always leads to better compliance-type structural 5 performance (e.g. structural compliance and thermal compliance), as both material type 1 ( ) and type 2 , and the overlapping area consumes double material ( ) are stiffer than the artificial material type ( volumes. 3. Optimization problem formulation and solution Based on the multi-material level set interpolation in Eq. (10), the compliance-minimization topology optimization problem is formulated below: ,Φ ,Φ . Φ 1 Φ Ω Φ 1 Φ Ω Φ Φ Ω , ,Φ ,Φ . . Φ ,∀ ∈ Ω (12) Φ , ,Φ ,Φ 1 Φ Φ 1 Φ ∙ where Ω , , Φ is the energy bilinear form and Φ is the displacement and Ω Ω Φ Ω Ω ∙ is the load linear form; is the space of kinematically admissible displacement field. traction force. Φ is the test vector, and is the body force and is the boundary represents the strain. 6 in To solve this problem, sensitivity analysis is required to calculated the boundary propagating speed direction of n = - . And design update is performed by solving the Hamilton-Jacobi equation as | | demonstrated in Eq. (13). Φ | Φ | (13) Specifically, the Lagrange multiplier method is employed and the Lagrangian function is constructed as: , ,Φ ,Φ ,Φ ,Φ where is the Lagrange multiplier and Ω Φ ∈ (14) is the adjoint displacement field. It is noted, the adjoint sensitivity analysis subject to multi-level set interpolation is widely studied [25,29,31], and therefore would not be further specified here. Only the results are demonstrated. 1 | Ω | Φ Φ Φ Φ | Φ Φ Φ 1 | Ω (15) Φ Φ Φ Then, by following Eq. (16), (16) could be guaranteed to change in the descent direction, as shown in Eq. (17), ′ Φ | Φ | Ω During the derivations, it is assumed that domain D is fixed. Φ , , , and | Φ | Ω 0 (17) are spatially invariant and the design In addition, the augmented Lagrange multiplier is adopted to address the material volume fraction constraints as: Φ Ω (18) 0 It is iteratively updated through Eq. (18) in which 1 indicates the iteration number. 7 4. Length scale control Length scale control has been a long-lasting issue tightly bonded to the part manufacturability, e.g., very small components are difficult to manufacture and small component length scale gradient is preferred by cooling-included manufacturing processes. Therefore, length scale control has been extensively studied during the last two decades, and approaches have been proposed under both the density-based [45–53] and level set[41–43,54–60] frameworks. Literature surveys on length scale control can be found in [61,62]. On the other hand, a gap can be identified that length scale control has rarely been performed based on multi-material topology optimization. Given the reasons for level set method, length scale control conventionally relies heavily on the signed distance information [41–43] but overlap of the level set function of the existing multi-material methods causes troubles in accessing the signed distance information in each material phase. This limitation motives us to represent each material phase by a single level set function and apply penalization to vanish the overlapping areas, as specified in the Section 2. The component length scale was previously defined in [41,42,57] and this concept is similarly defined in this work. The component length scale at a boundary point is defined by its shortest distance to the structural skeleton and the structural skeleton, also named the medial axis, is constructed by the set of centers of the maximally inscribed circles in 2D (and balls in 3D). Specifically, the length scale control functional proposed in [56] is adopted, as shown below: Φ Φ Φ 2 ,0 ; The notations: where Φ 2 Ω (19) ,0 is the designated length scale target of the ith material phase. The physical meaning of this functional is trivial to explain. The material areas of either larger or smaller than the designated length scale will be penalized. Put this control functional into the objective function in Eq. (12) and there we get: . ,Φ ,Φ Φ (20) where is the weighting factor balancing the length scale control effect and the structural compliance. It is noted that is applied as a non-zero value only after N iterations of the regular optimization loops, in order to let the structural topology freely evolve to derive a close-to-optimal result. Sensitivity analysis on the length scale control functional was previous studied in [56], and derivative of the thickness-control functional is presented in Eq. (21). Φ Φ 2 Φ 2 Φ Φ′ Ω (21) 8 2 Φ 2 Φ 2 Φ Φ′ Ω 2 formed into In Eq. (221), the secoond term on the right sidde is in fieldd integrationn and needs to be transfo boundaryy integrationn for the lateer Hamilton--Jacobi equaation based ddesign updatte. As the baasis of this transform mation, a corrollary is cited from [26]: Corollarry 1. For a 2D D integrable function f Φ Φ as demonsstrated in Fig. 2, Φ Ω 1 Γ (22) ∩ Two deffinitions are ppresented bellow to interprret this corolllary. ∈ Definitioon 1. For anyy of on Ω. When Π | | iss the set of prrojections ,Π ∈ Ω, | ≔ ∈ | reducees to a singlee point, it is ccalled the proojection of ontoo Ω. ≔ Definitioon 2. For anyy ∈ Ω, the ray eemerging from m . ∈ , is Figure 2 Scchematic plott of corollaryy 1 In Corolllary 1, representss the distancce which eqquals to Φ in case thhat the signeed distance regulatioon is satisfiedd. By applyinng Corollary 1, the sensitiivity result iss adapted intoo: ′ Φ Φ′ Ω Φ 2 Φ ∩ Φ 2 2 2 Φ (21) 2 2 1 Φ where Ω is the ith maaterial phase and Ω is itts boundary. 9 The analytical solution of Eq. (21) can be derived as shown below: 2 3 1 2 3 1 4 2 2 , 2 12 4 , (22) 2 In Eq. (22), represents the transient length scale value measured at the boundary point , and represents the related curvature value. Calculation of is straightforward [43] but lacks accuracy [57]. It is noted that, by ignoring the term, similar length scale control effects can be achieved where the updated sensitivity result is shown below: , 4 2 , 4 (23) 2 The function of Eq. (19) is to constrain the component length scale approaching a constant thickness. In fact, this function can be adapted to provide versatile length scale control effects, e.g. the maximum length scale only, the minimum length scale only, and concurrently both. Specifically, Eq. (19) is changed into: Φ where and Φ Φ 2 2 Φ Ω (24) are the upper and lower bounds of the length scale control. Then, the sensitivity result is changed into: ′ Φ Φ′ Ω , 4 0, 2 2 (25) 2 4 , 2 Another point worth noticing is that, there involves several max/min operations in the length scale control functional as shown in Eq. (19) which seems to be a non-desirable numerical feature. However, the sensitivity result is shown in Eq. (21) and later being simplified into Eq. (23). Eq. (23) is of C1 continuity even though it is a piecewise function. Hence, the involvement of the max/min operations in Eq. (19) does not harm the numerical solution. So far, the only unknown is then the transient length scale . In this paper, we built an auxiliary level set field to calculate the . Specifically, is calculated by measuring the shortest distance from the boundary point to the structural skeleton [41,42,57]. Hence, the structural skeleton is captured through the following procedures [42]. 10 (1) Identify the full skeleton Ω based on Eq. (26). Ω | Φ min ,Φ Φ 0 (26) where is a small positive number to guarantee the numerical computation robustness which is assigned the value 0.1 in this work. Laplacian of the level set function is calculated according to Eq. (27). Φ Φ Φ , Φ, , (2) Delete the skeleton points belonging to the corner area ̅ Ω Φ, 4Φ , (27) Ω , i.e., Ω / Ω (28) where, Ω in which, and |0.5 0.5 | Φ Φ ∙ Φ |∙| Φ | ,Φ (29) represent backward and forward, respectively. ∈ 0,1 works based on the fact that the backward and forward gradients employ opposite signs. is a positive number which determines the size of the corner area. The purpose of deleting these skeleton points is to ensure the proper length scale measure around the corner areas. Once the structural skeleton has been identified, the points located in the skeleton area ̅ are assigned the level set value of zero. Then, by solving Eq. (30), the signed distance field is re-initialized and the auxiliary level set field is created. The auxiliary level set value of any boundary point indicates its shortest distance to the structural skeleton; in other words, the value is obtained. | Φ | 1 (30) Going through the procedures to identify the structural skeleton and create the auxiliary level set field takes extra computational time. However, the impact is limited since there only involves calculation of the local finite difference information and Eq. (30) can be efficiently solved through the highly efficient fast marching method [63]. The extra procedures add around 14 percent of the overall computational time for 2D examples and only 1.5 percent for 3D examples. 5. Numerical examples In this section, several numerical examples are studied to prove the effectiveness of the proposed multimaterial level set topology optimization method and the associated length scale control technique. In all 0.5E , where E the following examples, the artificial weak material has the elastic tensor E indicates the elastic tensor of the weak material; and the material properties in the void area is E E ). 10 E 11 5.1 Three point-loaded beam Figure 3. Initial design domain and boundary conditions (size:150*75) In this example, the beam structure is optimized. The initial design domain and boundary conditions are demonstrated in Fig. 3, where three point forces of total magnitude 1 are loaded at the bottom, the left bottom corner is fully clamped and the right bottom corner is only fixed in the vertical direction. Two compositing materials will be considered: the strong material has a Young’s modulus of 1.30 and a Poisson’s ratio of 0.4, while the weak material has the Young’s modulus of 0.65 and the same Poisson’s ratio. Besides, the maximum volume fraction of the strong material is 0.25 and that of the weak material is 0.15. Correspondingly, the optimization process without length scale control is demonstrated in Fig. 4 and the convergence histories are shown in Fig. 5. The objective value converges at 5.67. The two materials form different components of the structure since they have the comparable properties. In addition, the strong material mainly distribute around the loading areas and the clamped areas where the deformation tends to be larger. Figure 4. Optimization process without length scale control (Blue: strong material; Green: weak material; Yellow: the artificial material type) 12 (a) (b) Figure 5. Convergence histories (a) Evolution history of the objective value (Final value = 5.67) (b) Evolution history of the volume fractions Other than that, we further tested this example by imposing different length scale limits: T = 6 and 8 respectively, and it is noted that, the weight factor in this case is defined in Eq. (31). ∙ 5 ∙ 10 (31) where iterNum is short for the optimization process iteration number, is the iteration number when the length scale control starts to take effect and H means the explicit Heaviside function. For the skeleton identification, 2 0.5, 3 2.5. Correspondingly, the optimization results with length scale control are demonstrated in Fig. 6. It can be concluded from the results that, i) the length scale targets have been well satisfied even though multiple material phases are involved; ii) imposing the length scale control compromises the structural compliance because of the constrained design space; and iii) the volume fractions of material phase 1 and 2 do not meet the upper bounds. More specifically, for the T = 6 case, a set of different values are explored to see the effect of choosing different formulations for . We can see from Fig. 6b and 6c that, the optimization results with 80 and 100 have the identical topology while only the shape is slightly different; comparatively, the optimization result with 120 has a different structural topology. For the reason, the length scale 13 control effort actually has the effect of blocking topologically changes since some structural component needs to approach zero thickness to make the topological change happen. Hence, the value should not be too small, since it may block most of the necessary topological changes. In addition, the objective values of Fig. 6(b-d) are at the same level and there is no major difference. Hence, it is recommended that can be selected within a range (such as 80-120 for this case). A too large is not suggested since it will reduce the optimization efficiency. In summary, there exist some rules of the value selection; however, the most appropriate value is case specific and the determination relies on some trials and errors, which can be summarized as a limitation of the current method. Other than that, the same limitation exists in choosing the value (e.g., 5 ∙ 10 ) which also relies on some trials. In addition, imposing length scale control will inevitably result in numerous local optima, which may explain why the results do not approach the volume fraction upper limits. The uniform length-scale constraint is very strict and results is a much smaller feasible design space, compared with only constraining the maximum or minimum length scale. Within this smaller feasible design space, the optimal design only has the strict uniform length scale constraint active. Another possible reason is the initial design dependency feature of the level set method [41]. (a) (b) (c) (d) Figure 6. Optimization results with length scale control: (a) T = 8, 100, objective value = 6.69, volume fraction of the strong material is 0.232 and that of the weak material is 0.150; (b) T = 6, 80, objective value = 11.66, volume fraction of the strong material is 0.208 and that of the weak material is 0.102; (c) T = 6, 100, objective value = 11.59, volume fraction of the strong material is 0.207 and that of the weak material is 0.100; (d) T = 6, 120, objective value = 11.82, volume fraction of the strong material is 0.174 and that of the weak material is 0.102 14 m 5.2 Two-bar cantilevver problem Figgure 7. Initiall design dom main and bounndary conditiions (size:150*75) mized. The initial design domain d and In this exxample, the ttwo-bar cantiilever structuure is topologgically optim boundaryy conditions are demonsstrated in Figg. 7, where a point forcce of magnittude 1 is loaaded at the bottom ccenter in the horizontal ddirection, andd the top edgge is fully clamped. Twoo compositingg materials will be cconsidered: thhe strong maaterial has a Y Young’s moddulus of 1.300 and a Poissson’s ratio off 0.4, while the weakk material haas the Youngg’s modulus oof 0.65 and tthe same Poissson’s ratio. Besides, thee maximum volume ffraction of thhe strong matterial is 0.15 and that of thhe weak matterial is 0.1. The optimization prooblem is studdied from twoo starts, as deemonstrated in Fig. 8(a, cc) and corresspondingly, the resullts without leength scale coontrol are preesented in Fiig. 8(b, d). T The convergennce histories are shown in Fig. 9. It can be cooncluded thatt the optimization result ccan be very dependent d onn the initial design since there arre clear disttinctions bettween the ttwo optimization resultss. Besides, the trend of o material distributiion is similarr to the last eexample wheere the strongg material m mainly distribuutes around the t loading and clam mped areas. The optimizzation result is different from [29] which w means the optimizzed designs could bee just local opptima. (a) (b) 15 (c) (d) Figure 8. Optimization results without length scale control (Blue: strong material; Green: weak material; Yellow: the artificial material type) (a) Start 1 of optimization process (b) Optimization result with start 1 (c) Start 2 of optimization process (d) Optimization result of start 2 (a) (b) 16 (c) (d) Figure 9. Convergence histories (a) Evolution history of the objective value with start 1 (Final value = 9.29) (b) Evolution history of the volume fractions with start 1 (c) Evolution history of the objective value with start 2 (Final value = 9.40) (d) Evolution history of the volume fractions with start 2 Then, we further tested this example by applying the length scale control and two schemes are studied: i) T = 8 for the strong material phase and T = 10 for the weak material phase; and ii) T = 10 for the strong material phase and T = 8 for the weak material phase. Start 1 as shown in Fig. 8 is used as the input. Noted that, the weight factor in this case is defined in Eq. (32). For the skeleton identification, 2 0.5, 3 4.5. 150 ∙ 2.5 ∙ 10 (32) The optimization results are demonstrated in Fig. 10. It can be concluded from the results that, i) the length scale targets have been satisfied even though multiple material phases are involved; ii) imposing the length scale control compromises the structural compliance because of the constrained design space; and iii) the volume fractions of the strong and weak material phases are only 0.070 and 0.100 for scheme 1 and 0.146 and 0.032 for scheme 2, which means that they do not meet the upper bounds. The convergence histories are plotted in Fig. 11. 17 (a) (b) Figu ure 10. Optim mization resuults with lenggth scale conttrol (a) Optim mization resuult of schemee 1 and objective vaalue = 14.23 (b) Optimizaation result of scheme 2 aand objective value = 12.000 (a) (b) (c) (d) Figure 11. Converggence historiees (a) Evolutiion history of the objectivve value of sccheme 1 (b) Evolution me fractions oof scheme 1 (c) Evolutionn history of thhe objective value of scheme 2 (d) history of the volum Evolutiion history off the volume fractions of scheme 2 18 m 5.3 3D ccube problem In this exxample, the 3D cube prooblem is studdied. The inittial design doomain and bboundary connditions are demonsttrated in Fig. 12, where a point force of magnitudde 5 is loadedd at the top ccenter and thhe four foot corners are fully claamped. Twoo compositinng materials will be connsidered: the strong mateerial has a Young’ss modulus off 1 and a Poissson’s ratio oof 0.3, whilee the weak material m has thhe Young’s m modulus of 0.5. Besiides, the maxximum volum me fraction oof the strong material is 00.15 and that of the weak material is 0.1. Figu ure 12. Input material disttribution andd boundary coonditions (Bllue: strong m material; Yelloow: the artifficial materiaal type) The inpuut material diistribution is also demonsstrated in Figg. 12, and thee optimizatioon result withhout length scale conntrol is show wn in Fig. 133. The multi--material disttribution is reasonable r since the stronng material forms the main load ppath while thhe weak mateerial forms thhe stiffening structure. 19 Figure 13. The non-constrained optimization result with objective value of 16.25 (Blue: strong material; Green: weak material; Yellow: the artificial material type) To perform the length scale control, we proposed a new functional in this example, where the length scale target need not be pre-specified; see below: ∗ Φ Φ ∗ Φ 2 Φ 2 Ω (33) ∗ 2 Φ ∗ Φ Ω/ Φ Ω Φ ∗ is the auxiliary level set field for material phase i. Physically, ∗ represents the length scale average at the specific optimization iteration. This upgraded length scale control functional penalizes the local length scales to approach an averaged value for each material phase, and therefore, uniform thickness distribution in each material phase could be expected, even though no length scale target is pre-specified. With this upgraded functional, the optimization problem is re-studied and we tested two sets of 0.6, 3.5) is shown in Fig. 14 and the result with parameters of Eq. (29). The result of using ( 0.8, 3.5) is presented in Fig. 15. The weight factor in this case is defined in Eq. (34). ( 30 ∙ 5 ∙ 10 (34) From the results, we can conclude that the uniform length scale has largely been realized within each material phase, even though local violations exist because of the multi-objective nature of the problem. On the other hand, some minor branches of the structural skeleton are captured with the parameter set ( 2 0.6, 3 3.5), and consequently, we can see the embossment in the strong material phase as shown in Fig. 14. These local structural details contribute better interactions between the two material phases, but it is non-intuitive to judge whether or not and how much they have violated the uniform length scale requirement. Comparatively, by raising the value of 2 0.8 to enforce the structural skeleton screening, these minor structural skeleton branches are removed from the identification result; see Fig. 15 where the embossment feature has been eliminated from the strong material phase. In addition, by applying the length scale control, the structural performance is slightly compromised compared to the original result. Noted that in Fig. 14 and 15, we used the open-source code from [64] to extract the iso-surface model and export the STL File. 20 (a) (b) (c) Figu ure 14. The ooptimization result with leength scale ccontrol with oobjective valuue of 17.10 ((a) The optim mization resultt with lengthh scale controol (b) Materiaal distributionn of the stronng and weak material phases, reespectively (cc) The extraccted iso-levell set surface 21 (a) (b) (c) Figu ure 15. The ooptimization result with leength scale ccontrol with oobjective valuue of 17.27 ((a) The optim mization resultt with lengthh scale controol (b) Materiaal distributionn of the stronng and weak material phases, reespectively (cc) The extraccted iso-levell set surface 22 6. Conclusion This paper presents a new multi-material level set topology optimization method, where each solid material phase is represented by a single level set function in the converged result. Therefore, the technical merits would be the no redundancy material phases and the well maintained signed distance field within each material phase. Because of the latter, length scale control has been realized on multimaterial level set problems, which to the best of the authors’ knowledge has never been achieved before, especially under the level set framework. Another contribution is that, we proposed a new length scale control functional through which the length scale target need not be pre-specified to realize the uniform component thickness. So far, the proposed method has been proved effective for compliance optimization; however, the extension to other types of optimization problems is still questionable which could be a limitation. It can be foreseen that, the proposed method is suitable to problems where the objective is inversely proportional to the material volume fraction and the material properties, e.g., structural compliance optimization and thermal compliance optimization. Given a more general scenario, the applicability is limited. For instance, for compliant mechanism design, the artificial weak material may form the hinges; and for stress-constrained problems, the artificial weak material may accumulate at the re-entrant corners to relieve the stress concentration. On the other hand, even though there is the limitation, the proposed method has its potential to be extended to even diversified optimization problems. Such as the compliant mechanism design, the objective function is often configured with multiple terms including the compliance term to eliminate potential hinges [65]. This strategy can be adopted by our method to design compliant mechanisms while avoiding the artificial weak material formed hinges. Hence, we will spend the effort to improve the proposed method for more general applicability. For future work, improvement of the multi-phase length scale control accuracy will be the main focus, as the structural skeleton identification is sensitive to the parameter changes and also the boundary perturbation. Acknowledgements The authors would like to acknowledge the support from China Scholarship Council (CSC) and Natural Sciences and Engineering Research Council of Canada (NSERC). Reference [1] M.P. Bendsøe, O. Sigmund, Topology Optimization, Springer Berlin Heidelberg, Berlin, Heidelberg, 2004. http://link.springer.com/10.1007/978-3-662-05086-6 (accessed May 26, 2016). [2] Y.M. Xie, G.P. Steven, A simple evolutionary procedure for structural optimization, Comput. Struct. 49 (1993) 885–896. doi:10.1016/0045-7949(93)90035-C. [3] G. Allaire, F. Jouve, A.-M. Toader, Structural optimization using sensitivity analysis and a level-set method, J. Comput. Phys. 194 (2004) 363–393. doi:10.1016/j.jcp.2003.09.032. [4] M.Y. Wang, X. Wang, D. Guo, A level set method for structural topology optimization, Comput. Methods Appl. Mech. Eng. 192 (2003) 227–246. doi:10.1016/S0045-7825(02)00559-5. 23 [5] [6] [7] [8] [9] [10] [11] [12] [13] [14] [15] [16] [17] [18] [19] [20] [21] [22] [23] [24] M.P. Bendsøe, O. Sigmund, Material interpolation schemes in topology optimization, Arch. Appl. Mech. 69 (1999) 635–654. doi:10.1007/s004190050248. O. Sigmund, S. Torquato, Design of materials with extreme thermal expansion using a three-phase topology optimization method, J. Mech. Phys. Solids. 45 (1997) 1037–1067. doi:10.1016/S00225096(96)00114-7. L.V. Gibiansky, O. Sigmund, Multiphase composites with extremal bulk modulus, J. Mech. Phys. Solids. 48 (2000) 461–498. doi:10.1016/S0022-5096(99)00043-5. O. Sigmund, Design of multiphysics actuators using topology optimization – Part II: Two-material structures, Comput. Methods Appl. Mech. Eng. 190 (2001) 6605–6627. doi:10.1016/S00457825(01)00252-3. C.F. Hvejsel, E. Lund, Material interpolation schemes for unified topology and multi-material optimization, Struct. Multidiscip. Optim. 43 (2011) 811–825. doi:10.1007/s00158-011-0625-z. A.H. Taheri, K. Suresh, An isogeometric approach to topology optimization of multi-material and functionally graded structures, Int. J. Numer. Methods Eng. (2016) n/a-n/a. doi:10.1002/nme.5303. Y. Luo, Y. Niu, M. Li, Z. Kang, A multi-material topology optimization approach for wrinkle-free design of cable-suspended membrane structures, Comput. Mech. (2017) 1–14. doi:10.1007/s00466017-1387-2. J. Stegmann, E. Lund, Discrete material optimization of general composite shell structures, Int. J. Numer. Methods Eng. 62 (2005) 2009–2027. doi:10.1002/nme.1259. E. Lund, J. Stegmann, On structural optimization of composite shell structures using a discrete constitutive parametrization, Wind Energy. 8 (2005) 109–124. doi:10.1002/we.132. E. Lund, Buckling topology optimization of laminated multi-material composite shell structures, Compos. Struct. 91 (2009) 158–167. doi:10.1016/j.compstruct.2009.04.046. T. Gao, W. Zhang, A mass constraint formulation for structural topology optimization with multiphase materials, Int. J. Numer. Methods Eng. 88 (2011) 774–796. doi:10.1002/nme.3197. M. Bruyneel, SFP—a new parameterization based on shape functions for optimal material selection: application to conventional composite plies, Struct. Multidiscip. Optim. 43 (2011) 17–27. doi:10.1007/s00158-010-0548-0. M. Bruyneel, P. Duysinx, C. Fleury, T. Gao, Extensions of the Shape Functions with Penalization Parameterization for Composite-Ply Optimization, AIAA J. 49 (2011) 2325–2329. doi:10.2514/1.J051225. T. Gao, W. Zhang, P. Duysinx, A bi-value coding parameterization scheme for the discrete optimal orientation design of the composite laminate, Int. J. Numer. Methods Eng. 91 (2012) 98–114. doi:10.1002/nme.4270. T. Gao, W.H. Zhang, P. Duysinx, Simultaneous design of structural layout and discrete fiber orientation using bi-value coding parameterization and volume constraint, Struct. Multidiscip. Optim. 48 (2013) 1075–1088. doi:10.1007/s00158-013-0948-z. M.Y. Wang, X. Wang, “Color” level sets: a multi-phase method for structural topology optimization with multiple materials, Comput. Methods Appl. Mech. Eng. 193 (2004) 469–496. doi:10.1016/j.cma.2003.10.008. Y. Mei, X. Wang, A level set method for structural topology optimization with multi-constraints and multi-materials, Acta Mech. Sin. 20 (2004) 507–518. doi:10.1007/BF02484273. M.Y. Wang, S. Chen, X. Wang, Y. Mei, Design of Multimaterial Compliant Mechanisms Using Level-Set Methods, J. Mech. Des. 127 (2005) 941–956. doi:10.1115/1.1909206. M.Y. Wang, X. Wang, A level-set based variational method for design and optimization of heterogeneous objects, Comput.-Aided Des. 37 (2005) 321–337. doi:10.1016/j.cad.2004.03.007. C. Zhuang, Z. Xiong, H. Ding, Topology optimization of multi-material for the heat conduction problem based on the level set method, Eng. Optim. 42 (2010) 811–831. doi:10.1080/03052150903443780. 24 [25] X. Guo, W. Zhang, W. Zhong, Stress-related topology optimization of continuum structures involving multi-phase materials, Comput. Methods Appl. Mech. Eng. 268 (2014) 632–655. doi:10.1016/j.cma.2013.10.003. [26] G. Allaire, C. Dapogny, G. Delgado, G. Michailidis, Multi-phase structural optimization via a level set method, ESAIM Control Optim. Calc. Var. 20 (2014) 576–611. [27] N. Vermaak, G. Michailidis, G. Parry, R. Estevez, G. Allaire, Y. Bréchet, Material interface effects on the topology optimizationof multi-phase structures using a level set method, Struct. Multidiscip. Optim. 50 (2014) 623–644. doi:10.1007/s00158-014-1074-2. [28] P. Liu, Y. Luo, Z. Kang, Multi-material topology optimization considering interface behavior via XFEM and level set method, Comput. Methods Appl. Mech. Eng. 308 (2016) 113–133. doi:10.1016/j.cma.2016.05.016. [29] Y. Wang, Z. Luo, Z. Kang, N. Zhang, A multi-material level set-based topology and shape optimization method, Comput. Methods Appl. Mech. Eng. 283 (2015) 1570–1586. doi:10.1016/j.cma.2014.11.002. [30] M. Cui, H. Chen, J. Zhou, A level-set based multi-material topology optimization method using a reaction diffusion equation, Comput.-Aided Des. 73 (2016) 41–52. doi:10.1016/j.cad.2015.12.002. [31] J. Liu, Y. Ma, Design of pipeline opening layout through level set topology optimization, Struct. Multidiscip. Optim. 55 (2017) 1613–1628. doi:10.1007/s00158-016-1602-3. [32] Y. Wang, Z. Luo, N. Zhang, T. Wu, Topological design for mechanical metamaterials using a multiphase level set method, Struct. Multidiscip. Optim. (2016) 1–16. doi:10.1007/s00158-0161458-6. [33] P. Wei, M.Y. Wang, Piecewise constant level set method for structural topology optimization, Int. J. Numer. Methods Eng. 78 (2009) 379–402. doi:10.1002/nme.2478. [34] S. Chen, S. Gonella, W. Chen, W.K. Liu, A level set approach for optimal design of smart energy harvesters, Comput. Methods Appl. Mech. Eng. 199 (2010) 2532–2543. doi:10.1016/j.cma.2010.04.008. [35] P. Vogiatzis, S. Chen, X. Wang, T. Li, L. Wang, Topology optimization of multi-material negative Poisson’s ratio metamaterials using a reconciled level set method, Comput.-Aided Des. 83 (2017) 15–32. doi:10.1016/j.cad.2016.09.009. [36] Z. Kang, Y. Wang, Y. Wang, Structural topology optimization with minimum distance control of multiphase embedded components by level set method, Comput. Methods Appl. Mech. Eng. 306 (2016) 299–318. doi:10.1016/j.cma.2016.04.001. [37] K. Ghabraie, An improved soft-kill BESO algorithm for optimal distribution of single or multiple material phases, Struct. Multidiscip. Optim. 52 (2015) 773–790. doi:10.1007/s00158-015-1268-2. [38] X. Huang, Y.M. Xie, Bi-directional evolutionary topology optimization of continuum structures with one or multiple materials, Comput. Mech. 43 (2009) 393. doi:10.1007/s00466-008-0312-0. [39] S. Zhou, M.Y. Wang, Multimaterial structural topology optimization with a generalized Cahn– Hilliard model of multiphase transition, Struct. Multidiscip. Optim. 33 (2007) 89. doi:10.1007/s00158-006-0035-9. [40] R. Tavakoli, Multimaterial topology optimization by volume constrained Allen–Cahn system and regularized projected steepest descent method, Comput. Methods Appl. Mech. Eng. 276 (2014) 534–565. doi:10.1016/j.cma.2014.04.005. [41] X. Guo, W. Zhang, W. Zhong, Explicit feature control in structural topology optimization via level set method, Comput. Methods Appl. Mech. Eng. 272 (2014) 354–378. doi:10.1016/j.cma.2014.01.010. [42] Q. Xia, T. Shi, Constraints of distance from boundary to skeleton: For the control of length scale in level set based structural topology optimization, Comput. Methods Appl. Mech. Eng. 295 (2015) 525–542. doi:10.1016/j.cma.2015.07.015. [43] J. Liu, H. Yu, Y. Ma, Minimum void length scale control in level set topology optimization subject to machining radii, Comput.-Aided Des. 81 (2016) 70–80. doi:10.1016/j.cad.2016.09.007. 25 [44] R.I. Saye, J.A. Sethian, The Voronoi Implicit Interface Method for computing multiphase physics, Proc. Natl. Acad. Sci. 108 (2011) 19498–19503. doi:10.1073/pnas.1111557108. [45] T.A. Poulsen, A new scheme for imposing a minimum length scale in topology optimization, Int. J. Numer. Methods Eng. 57 (2003) 741–760. doi:10.1002/nme.694. [46] J.K. Guest, J.H. Prévost, T. Belytschko, Achieving minimum length scale in topology optimization using nodal design variables and projection functions, Int. J. Numer. Methods Eng. 61 (2004) 238– 254. doi:10.1002/nme.1064. [47] J.K. Guest, Topology optimization with multiple phase projection, Comput. Methods Appl. Mech. Eng. 199 (2009) 123–135. doi:10.1016/j.cma.2009.09.023. [48] M. Schevenels, B.S. Lazarov, O. Sigmund, Robust topology optimization accounting for spatially varying manufacturing errors, Comput. Methods Appl. Mech. Eng. 200 (2011) 3613–3627. doi:10.1016/j.cma.2011.08.006. [49] O. Sigmund, Manufacturing tolerant topology optimization, Acta Mech. Sin. 25 (2009) 227–239. doi:10.1007/s10409-009-0240-z. [50] O. Sigmund, Morphology-based black and white filters for topology optimization, Struct. Multidiscip. Optim. 33 (2007) 401–424. doi:10.1007/s00158-006-0087-x. [51] F. Wang, B.S. Lazarov, O. Sigmund, On projection methods, convergence and robust formulations in topology optimization, Struct. Multidiscip. Optim. 43 (2011) 767–784. doi:10.1007/s00158-0100602-y. [52] M. Zhou, B.S. Lazarov, F. Wang, O. Sigmund, Minimum length scale in topology optimization by geometric constraints, Comput. Methods Appl. Mech. Eng. 293 (2015) 266–282. doi:10.1016/j.cma.2015.05.003. [53] W. Zhang, W. Zhong, X. Guo, An explicit length scale control approach in SIMP-based topology optimization, Comput. Methods Appl. Mech. Eng. 282 (2014) 71–86. doi:10.1016/j.cma.2014.08.027. [54] S. Chen, M.Y. Wang, A.Q. Liu, Shape feature control in structural topology optimization, Comput.Aided Des. 40 (2008) 951–962. doi:10.1016/j.cad.2008.07.004. [55] J. Luo, Z. Luo, S. Chen, L. Tong, M.Y. Wang, A new level set method for systematic design of hinge-free compliant mechanisms, Comput. Methods Appl. Mech. Eng. 198 (2008) 318–331. doi:10.1016/j.cma.2008.08.003. [56] J. Liu, Y. Ma, J. Fu, K. Duke, A novel CACD/CAD/CAE integrated design framework for fiberreinforced plastic parts, Adv. Eng. Softw. 87 (2015) 13–29. doi:10.1016/j.advengsoft.2015.04.013. [57] G. Allaire, F. Jouve, G. Michailidis, Thickness control in structural optimization via a level set method, Struct. Multidiscip. Optim. 53 (2016) 1349–1382. doi:10.1007/s00158-016-1453-y. [58] W. Zhang, D. Li, J. Zhang, X. Guo, Minimum length scale control in structural topology optimization based on the Moving Morphable Components (MMC) approach, Comput. Methods Appl. Mech. Eng. 311 (2016) 327–355. doi:10.1016/j.cma.2016.08.022. [59] Y. Wang, L. Zhang, M.Y. Wang, Length scale control for structural optimization by level sets, Comput. Methods Appl. Mech. Eng. 305 (2016) 891–909. doi:10.1016/j.cma.2016.03.037. [60] J. Liu, L. Li, Y. Ma, Uniform thickness control without pre-specifying the length scale target under the level set topology optimization framework, Adv. Eng. Softw. (2018). doi:10.1016/j.advengsoft.2017.09.013. [61] J. Liu, Y. Ma, A survey of manufacturing oriented topology optimization methods, Adv. Eng. Softw. 100 (2016) 161–175. doi:10.1016/j.advengsoft.2016.07.017. [62] B.S. Lazarov, F. Wang, O. Sigmund, Length scale and manufacturability in density-based topology optimization, Arch. Appl. Mech. 86 (2016) 189–218. doi:10.1007/s00419-015-1106-4. [63] J. Sethian, Fast Marching Methods, SIAM Rev. 41 (1999) 199–235. doi:10.1137/S0036144598347059. [64] P. Vogiatzis, S. Chen, C. Zhou, An Open Source Framework For Integrated Additive Manufacturing And Level-Set Based Topology Optimization, ASME J. Comput. Inf. Sci. Eng. (2017). 26 [65] Q. Xia, T. Shi, Topology optimization of compliant mechanism and its support through a level set method, Comput. Methods Appl. Mech. Eng. 305 (2016) 359–375. doi:10.1016/j.cma.2016.03.017. Figure caption list: Figure 1. Multi-material level set schemes Figure 2 Schematic plot of corollary 1 Figure 3. Initial design domain and boundary conditions (size:150*75) Figure 4. Optimization process without length scale control (Blue: strong material; Green: weak material; Yellow: the artificial material type) Figure 5. Convergence histories (a) Evolution history of the objective value (Final value = 5.67) (b) Evolution history of the volume fractions Figure 6. Optimization results with length scale control: (a) T = 8, 100, objective value = 6.69, volume fraction of the strong material is 0.232 and that of the weak material is 0.150; (b) T = 6, 80, objective value = 11.66, volume fraction of the strong material is 0.208 and that of the weak material is 0.102; (c) T = 6, 100, objective value = 11.59, volume fraction of the strong material is 0.207 and that of the weak material is 0.100; (d) T = 6, 120, objective value = 11.82, volume fraction of the strong material is 0.174 and that of the weak material is 0.102 Figure 7. Initial design domain and boundary conditions (size:150*75) Figure 8. Optimization results without length scale control (Blue: strong material; Green: weak material; Yellow: the artificial material type) (a) Start 1 of optimization process (b) Optimization result with start 1 (c) Start 2 of optimization process (d) Optimization result of start 2 Figure 9. Convergence histories (a) Evolution history of the objective value with start 1 (Final value = 9.29)(b) Evolution history of the volume fractions with start 1(c) Evolution history of the objective value with start 2 (Final value = 9.40) (d) Evolution history of the volume fractions with start 2 Figure 10. Optimization results with length scale control (a) Optimization result of scheme 1 and objective value = 14.23 (b) Optimization result of scheme 2 and objective value = 12.00 Figure 11. Convergence histories (a) Evolution history of the objective value of scheme 1 (b) Evolution history of the volume fractions of scheme 1 (c) Evolution history of the objective value of scheme 2 (d) Evolution history of the volume fractions of scheme 2 Figure 12. Input material distribution and boundary conditions (Blue: strong material; Yellow: the artificial material type) Figure 13. The non-constrained optimization result with objective value of 16.25 (Blue: strong material; Green: weak material; Yellow: the artificial material type) 27 Figure 14. The optimization result with length scale control with objective value of 17.10 (a) The optimization result with length scale control (b) Material distribution of the strong and weak material phases, respectively (c) The extracted iso-level set surface Figure 15. The optimization result with length scale control with objective value of 17.27 (a) The optimization result with length scale control (b) Material distribution of the strong and weak material phases, respectively (c) The extracted iso-level set surface 28 Highlights Develop a new multi-material level set topology optimization method which eliminates the overlapping of level set functions for multi-material interpolation. Maintain the signed distance characteristic within each material phase. Realize independent component length scale control for each material phase. Propose a new functional which realizes the uniform component length scale control without need of pre-specifying the length scale target. 29

1/--страниц