DESIGN OPTIMIZATION OF ADDITIVELY MANUFACTURED LATTICE STRUCTURES ACCOUNTING FOR MANUFACTURING-INDUCED MATERIAL ORTHOTROPY AND MULTIPHYSICS LOADING By Ming Liu A DISSERTATION Submitted to Michigan State University in partial fulfillment of the requirements for the degree of Mechanical Engineering — Doctor of Philosophy 2022 ABSTRACT Additively manufactured (AM) lattice structures are designed and optimized while accounting for the effects of manufacturing-induced material orthotropy, the resolution limitations of AM, and multiphysics (structural and thermal) loading conditions. A model is formulated for general 3D orthotropy in each member of an AM lattice structure, with principal material coordinates associated with the AM build direction. This model is used within a bi-level design optimization strategy to simultaneously optimize the radius of each lattice member and the manufacturing build direction. Design study results illustrate the importance of accounting for material anisotropy in the design of AM lattice structures and the benefits of simultaneously optimizing the lattice and the manufacturing build direction. To account for AM resolution limitations, a model-based approach is implemented that accounts for geometrical tolerances while minimizing the geometric differences between the as- designed structure and the as-manufactured one. The influences of both material orthotropy and manufacturing resolution are studied and compared, demonstrating that both printing direction and process resolution play crucial roles in guaranteeing the performance and manufacturability of the as-manufactured lattice structure. For lattice design under multiphysics fields, a method is proposed to generate lattice structure designs under multiphysics loading, including stress, heat conduction, and heat convection. For a given set of loading and boundary conditions, a lattice trajectory is generated with an orientation vector selection algorithm based on the principal vectors of each field. Considering the design bias for a given physics field, weighting factors are assigned for each field to generate the lattice layout. A set of design studies are conducted to demonstrate the impact of different combinations of weighting factors on the lattice layout. It is shown that the proposed approach is capable of generating efficient lattice designs that represent effective trade-offs between mechanical and thermal load-carrying capabilities. ACKNOWLEDGEMENTS I would like to express my gratitude to my advisor, Dr. Ronald Averill, for his guidance, trust, and support during my Ph.D. study. Thank you for not only teaching me what research is and how to solve problems, but also for encouraging me and pushing me forward to try something new. Also, I could not have taken this journey without my defense committee. Thank you, Dr. Kalyanmoy Deb, for teaching me the foundation of multi-objective optimization and generously offering invaluable advice and knowledge in my research. Thank you, Dr. Erik Goodman, for sharing knowledge and demonstrating passion and diligence for the work. Thank you, Dr. Alejandro Diaz, for providing expertise and suggestions to improve my research and educating me to be meticulous with research. I am extremely grateful to have this opportunity to work with each of you. Not only did your intelligence and expertise provide invaluable advice and make a significant impact on me, but also your diligence, work ethic, and personal characteristics inspired me to become someone like you. I would also like to thank Matthew Lee Ryerkerk, Abhiroop Ghosh, Qiren Gao, Zhenxiang Jiang, and Erva Ulu for their support and help along the way. Thanks to all of our collaborators for their contributions to the TRADES project. Thanks also go to all professors and instructors in the Mechanical Engineering department at Michigan State University. I am grateful for the opportunity to learn new topics, to experience how to lead labs, and for all the wonderful classmates and colleagues I've met along the way. My appreciation also goes out to my family, especially my parents and my grandmother. Without their tremendous understanding and encouragement over the past few years, it would be impossible for me to complete my study. I would also like to express my deepest memory of my iv grandmother. I am particularly grateful to my grandmother, who imbued in me the spirit that has empowered me to reach this point. I am thankful for a lot of friends. Thanks to them for their company, loving-kindness, hospitality, and encouragement, which supported me and made me feel like I was with family all the time. This research was developed with funding from the Defense Advanced Research Projects Agency (DARPA). The views, opinions, and/or findings expressed are those of the author and should not be interpreted as representing the official views or policies of the Department of Defense or the U.S. Government. v TABLE OF CONTENTS LIST OF ABBREVIATIONS .................................................................................................... vii Chapter 1 Introduction................................................................................................................ 1 1.1 Motivation ..................................................................................................................... 1 1.2 Dissertation Outline ...................................................................................................... 2 Chapter 2 Background ................................................................................................................ 4 2.1 Anisotropic Material Properties in AM Specimens ........................................................ 4 2.2 Lattice Structure Design .................................................................................................. 5 2.3 Lattice Structure Optimization ........................................................................................ 5 2.4 Manufacturability of AM Process ................................................................................... 6 2.5 Multiphysics Lattice Design............................................................................................ 8 Chapter 3 Effect of Manufacturing-Induced Anisotropy on the Design of Additively Manufactured Lattice Structures ........................................................................... 10 3.1 Methodology ................................................................................................................. 10 3.2 Model Validation........................................................................................................... 22 3.3 Studies and Results........................................................................................................ 26 Chapter 4 Optimizing the Manufacturability and Performance of Additively Manufactured Lattice Structures ........................................................................... 42 4.1 Methodology ................................................................................................................. 42 4.2 Studies and Results........................................................................................................ 45 Chapter 5 Multiphysics Design Optimization of Lattice Structures ..................................... 59 5.1 Intermediate Methods and Results ................................................................................ 59 5.2 Methodology ................................................................................................................. 70 5.3 Studies and Results........................................................................................................ 78 Chapter 6 Summary, Conclusions, and Future Work............................................................ 89 6.1 Summary and Conclusions ............................................................................................ 89 6.2 Future Work .................................................................................................................. 90 BIBLIOGRAPHY ....................................................................................................................... 92 vi LIST OF ABBREVIATIONS FEA Finite Element Analysis FEM Finite Element Method AM Additive manufacturing GCS Global coordinate system LCS Local Coordinate System MCS Material Coordinate System TO Topology Optimization FDM Fused Deposition Modeling SLM Selective Laser Melting vii Chapter 1 Introduction 1.1 Motivation Additive manufacturing (AM) allows much greater design freedom compared to subtractive manufacturing strategies [1], and AM is well suited to producing lattice structures. These structures are sometimes constructed with repeating unit cells, within which differently oriented struts are interconnected [2], and advanced design strategies allow the lattice member orientations to be locally optimized to align with dominant stress or thermal fields [3]–[6]. Compared to traditionally manufactured solid structures, lattice structures have many potential advantages, such as high strength and stiffness-to-weight ratios, efficient heat dissipation due to the high porosity and large surface area, tunable energy absorption [7], and more. The mechanical properties of a lattice structure can be tailored to meet special system requirements and/or local performance requirements that may vary from region to region within a structure. Considering these performance advantages and their material and energy-saving properties, lattice structures hold great potential for enabling novel designs of advanced structures, and for re-designing many existing systems [8]. Due to the mechanism of AM process, some forms of AM induce local material anisotropy. For example, the laser powder bed fusion process produces parts with different mechanical properties in the build direction compared to those properties in the printing plane [9]. The printing direction controls the principal directions of the resulting orthotropy and thus can have a 1 significant effect on the performance of the AM structure. Several experimental studies have demonstrated the material strength dependence on the build direction [10], [11]. The design, manufacturing, and optimization of AM lattice structures have been investigated by many researchers. However, material anisotropy has seldom been considered. The performance of AM structures is also known to be affected by the printing direction (or the orientation of the part). Previous studies [12], [13] show that the printing direction can have a significant effect on performance, with the strength of parts varying with printing direction [14], [15]. The principal material directions of the manufacturing-induced material orthotropy are associated with the printing direction [16]. It is necessary to take into consideration this material orthotropy in the design optimization of AM lattice structures [17]. Advanced design strategies allow the lattice member orientations to be locally optimized to align with dominant stress or thermal fields [5], [6]. Compared to traditionally manufactured solid structures, lattice structures have many potential advantages such as high strength- or stiffness-to-weight ratios, efficient heat dissipation due to the high porosity and large surface area, tunable energy absorption [7], and more. The mechanical properties of a lattice structure can be tailored to meet special system requirements and/or local performance requirements that may vary from region to region within a structure. Considering these performance advantages and their material and energy-saving properties, lattice structures hold great potential for enabling novel designs of advanced structures, and for re-designing many existing systems [8]. 1.2 Dissertation Outline The remainder of the dissertation is structured as follows. A background and literature review of existing studies on lattice structure design and optimization is provided in Chapter 2. The study on the effect of manufacturing-induced anisotropy on the design of additively 2 manufactured lattice structure optimization is presented in Chapter 3. Optimization studies of lattice structure incorporated both geometry correction and material orthotropy are presented in Chapter 4. In Chapter 5, the study of design optimization of lattice structure layout under thermo-mechanical loadings is demonstrated. Finally, conclusions are drawn, and future work directions are provided in Chapter 6. 3 Chapter 2 Background 2.1 Anisotropic Material Properties in AM Specimens The anisotropic mechanical properties of various AM materials and processes have been investigated broadly. Ahn et al. [18] investigated the anisotropic material properties for fused deposition modeling ABS structures with tensile and compressive tests for different raster direction specimens, and the results showed that both tensile and compressive strengths for transverse specimens are lower than those for axial specimens. In addition to tension and compression testing, Ziemian et al. [19] conducted flexural testing, impact testing, and fatigue testing for ABS specimens, finding that the tensile and bending yield strengths are largest for a 0-degree raster orientation, and a 45-degree raster orientation is weakest for compression. With the constitutive material models developed for fused deposition modeling printed parts using a numerical homogenization procedure, Somireddy et al. [20] investigated the influence of AM process parameters and printing directions on the material properties. For 3D printed ABS and polycarbonate (PC), Cantrell et al. [21] characterized both tensile and shear modulus and yield strength, finding that anisotropic properties varied as the raster orientation and build orientation changed. For the directed energy deposition additive manufacturing process, Carroll et al. [22] studied the anisotropic tensile behavior of Ti-6Al-4V material components. The anisotropy was also characterized for laser powder bed fusion processed materials [9]. Within this study, the material anisotropy was calculated based on the classic laminate theory for stress analysis. 4 2.2 Lattice Structure Design Several investigations have focused on new design strategies for AM lattice structures. Assuming a periodic distribution of unit cells in regions of the lattice structure, Tao et al. [23] reviewed design methods for the lattice structure design process, which could be classified into unit cell design and pattern design. For the unit cell design, three methods were summarized: primitive-based method, implicit surface method, and topology optimization method. For the pattern design, direct patterning, conformal patterning, and topology optimization were utilized for distributing the unit cells to the structure. The phases of the design process of additively manufactured mesoscale lattice structures were also reviewed by Tamburrino et al. [24]. A novel approach was investigated by Daynes et al. [4] to generate lattice structures aligning with the strain trajectories, yielding lattice structures with better stiffness and strength. The design of the lattice structure focused not only on the geometry but also on the material properties, manufacturing testing, and analysis. The material anisotropy was accounted for in the lattice design, and new methods were proposed for designing the lattice structure with the prescribed anisotropic properties [25]. 2.3 Lattice Structure Optimization A large number of lattice optimization studies have been conducted based on topology optimization. The hexagonal unit cells have been distributed with the energy-based homogenization filter [26]. The topology optimization has also been applied with other methods for AM lattice design. The cooling channel design is based on lattice structure topology optimization and the movable design-dependent feature [27]. For AM lattice structure fabrication, the Bidirectional Evolutionary Structural Optimization (BESO) based design method was also got employed [28]. 5 From the material perspective, Stanković et al. [29] conducted optimization on the anisotropic material properties of AM lattice structures. Using the formulated anisotropy material properties related to the build orientation, the effect of anisotropy on the performance of the lattice structure was studied. Additional investigations studied the effect of anisotropy with the inkjet 3D printing model, in which they implemented an algorithm for modeling anisotropic compression and tension states into the Generalized Optimality Criteria (GOC). The study demonstrates the importance of considering the anisotropy effect and printing direction influence within the lattice structure design [30]. 2.4 Manufacturability of AM Process Due to the effects of AM process parameters [31], [32] such as manufacturing resolution and minimum printable feature size, the shape of as-manufactured parts can be different than as- designed parts [33], [34]. These AM process parameters can also introduce accuracy errors to the manufactured lattice structure[35]. For example, at a sharp corner of the structure, as the minimum printable feature size of the process increases the variation between the finished product and the designed model also increases. Lattice structures that have small beam members or small angles at beam intersections also present a challenge to AM, which may not be able to accurately produce these geometries due to lack of resolution. Each AM process has a typical resolution limit [36]. For fused deposition melting (FDM), the resolution is normally around 0.3mm [37], which is used in this study. Stereolithography (SLA) has a higher resolution than FDM at 0.05mm. For metal additive manufacturing processes, the minimum printing feature size for Selective Laser Sintering (SLS) is 0.762mm, and for Electron Beam Melting (EBM) is 0.1mm. 6 To account for these geometric limitations associated with an AM process, a manufacturability analysis based on model correction [38] is used here within the design process. The model correction analysis first slices the designed AM geometry into voxels. Then, the critical features affected by the given AM process are identified by considering the resolution or the minimum printable feature of the AM process. Modifications are based on the topology of the original design geometry and the toolpath of the print head. The performance of AM structures is also known to be affected by the printing direction (or the orientation of the part). Previous studies [12], [13] show that the printing direction can have a significant effect on performance, with the strength of parts varying with printing direction [14], [15]. The principal material directions of the manufacturing-induced material orthotropy are associated with the printing direction [16]. It is necessary to take into consideration this material orthotropy in the design optimization of AM lattice structures [17]. Additive manufacturing (AM), as a newly developing technology, is raising more possibilities and ideas for the manufacturing area. The complexity of the fabricated part also increases for more possibilities and adaptivity. For the manufacturability of the structure fabricated by the additive manufacturing process, the printing direction (or the orientation of the model) for the process is a crucial setting. From the studies [11], [39], the printing direction has an undeniable effect on the performance of the fabricated structure. With the mechanism of additive manufacturing, the strength of the fabricated parts would vary with different printing directions and the performance of the structure could change accordingly[18]. The material orthotropy gets introduced and is correlated to the printing direction. For the manufacturability of additively manufactured structures, it is necessary to take into consideration the material orthotropy of various printing directions [40]. The former study 7 includes the additive manufacturing (AM) process that introduced material orthotropy in the optimization of the additively manufactured lattice structure. From another perspective of lattice structure manufacturability, the restriction of the AM processes could also introduce errors to the manufactured parts[41]. Like the sharp corner of the structure, as the resolution of the process decreased, the variation between the finished product and the designed model would be larger. Especially for the lattice structure, the size of beams could be small and the intersections between beams would have various angles, some of the lattice beams may lose the material strength since the size is smaller than the resolution, or the fabricated corner would lose its geometry accuracy. 2.5 Multiphysics Lattice Design The development of Additive Manufacturing (AM) technology creates new possibilities for complex yet manufacturable geometries, wherein the scale of the manufactured part varies from the microscale unit cell to macroscale structures, and the manufacturing materials change from a single material to multiple materials [42]. Lattice structures have the advantages of a high strength-to-weight ratio, versatility, and broad applications, which have attracted much attention in the AM field of research. Multifunctional lattice structures [43] were designed for support structures [44], heat dissipation [45], energy absorption [46], negative Poisson's ratio [47], bioimplants [48], and aerospace components [49]. The design and optimization of lattice structures could be categorized into periodic lattice structures and stochastic lattice structures[50]. The research on periodic lattice structure is focused on the properties of the different lattice unit cells to improve the performance of the lattice structure[51]. Due to the repetitive property of the lattice structure, a large number of studies applied topology optimization (TO) to design and optimize the design with lattice unit 8 cells [52], [53]. TO method could minimize the total material and maximize the stiffness of the design part by changing the material distribution of the meshed design domain according to the loading conditions. For lattice structure, the lattice unit cell could be regarded as the meshed element which could be assigned to the remaining region of the TO domain. For ununiform lattice structures, another lattice design approach is based on the load path or trajectory of the prescribed field[54]–[57]. With a preprocessed Finite Element Analysis (FEA) of the given design domain, a trajectory based on the principal field could be obtained. Gao et al.,[5] developed an algorithm to generate a stress trajectory with principal stress values and vectors for a model under certain mechanical loading conditions and aligned the lattice beams along the stress trajectory to optimize the loading carrying efficiency of the lattice layout. Lattice structure was also investigated under certain thermal conditions. The thermal conductivity of lattice structures under different metallic AM processes was studied[58]. The lattice structure was also designed for AM heat conduction[27], enhancing thermal dissipation[59]. Jiang et al., [6] applied the field-aligned lattice algorithm [5] within the thermal field to generate the lattice layout aligning the heat flux trajectory, and optimized the lattice beam radii with the Method of Moving Asymptotes (MMA) [60] to minimum compliance. As a promising multifunctional structure, lattice structure could be beneficial for applications in multi-physics fields, which is one of the research areas that need to draw more attention. The uniform cellular material composed of thermoelastic structures was optimized for under mechanical and thermal loads[61], [62]. With the multi-physics topology optimization method, the graded porous structures were optimized for structural and thermal performance [63]. 9 Chapter 3 Effect of Manufacturing-Induced Anisotropy on the Design of Additively Manufactured Lattice Structures In this chapter, the material anisotropy induced by the AM process is modeled within each lattice member, and a bi-level optimization strategy is employed to simultaneously optimize the lattice and the AM build direction. Application examples demonstrate the influence of material anisotropy on the resulting lattice structure design as well as the importance of simultaneously optimizing the build direction and the lattice properties. 3.1 Methodology In the current model, three different coordinate systems are used to fully describe the orthotropic material behavior in a lattice member: 1. Global coordinate system (GCS): The global coordinate system is defined at the lattice structure level. The global coordinate axes are denoted by capital letters X, Y, Z. 2. Local coordinate system (LCS): A local coordinate system is defined for each beam element within the lattice structure. It may vary from beam to beam according to its orientation. The local coordinate axes are denoted by lowercase letters x, y, z, where x is the centroidal axis along the length of the beam, and the y and z axes are in the plane of the beam cross-section that is perpendicular to the x-axis. 3. Material coordinate system (MCS): The material coordinate system is defined within each layer of the material buildout. As the depositing direction of the nozzle changes, the orientation of the principal material coordinates for each beam will change as well. The 10 material coordinate axes are denoted by numbers 1, 2, 3, where axis 3 is defined as the nozzle’s depositing direction, which is normal to the 3D printing plane and often referred to as the “printing direction.” Note that material coordinate 3 is often, but not necessarily, defined globally for the entire part being manufactured. 3 Y 1 2 (MCS) x y (LCS) X (GCS) z Z Figure 3.1: Different coordinate systems shown for an AM beam. 3.1.1 Orthotropic Material Properties of 3D lattices 3.1.1.1 Stress-strain relations for an orthotropic material in the MCS From generalized Hooke’s law, the strain-stress relations can be expressed as: {𝜀𝜀} = [𝑆𝑆]{𝜎𝜎} (1) where [𝑆𝑆] is the compliance matrix, {𝜀𝜀} is the strain vector and {𝜎𝜎} is the stress vector. For a transversely isotropic material [64], (1) can be expressed in principal material coordinates as: 11 𝜀𝜀1 𝑆𝑆11 𝑆𝑆12 𝑆𝑆13 0 0 0 𝜎𝜎 ⎧ 𝜀𝜀2 ⎫ ⎡𝑆𝑆12 𝑆𝑆11 𝑆𝑆13 ⎤ ⎧ 𝜎𝜎1 ⎫ 0 0 0 ⎪ 𝜀𝜀 ⎪ ⎢𝑆𝑆 ⎥⎪ 2 ⎪ 3 ⎢ 13 𝑆𝑆13 𝑆𝑆33 0 0 0 ⎥ 𝜎𝜎3 𝛾𝛾23 = ⎢ 0 0 0 𝑆𝑆44 0 0 ⎥ ⎨𝜏𝜏23 ⎬ (2) ⎨ ⎬ ⎪𝛾𝛾31 ⎪ ⎢ 0 0 0 0 𝑆𝑆44 0 ⎥ ⎪𝜏𝜏31 ⎪ ⎩𝛾𝛾12 ⎭ ⎣ 0 0 0 0 0 𝑆𝑆66 ⎦ ⎩𝜏𝜏12 ⎭ The components of the compliance matrix can be expressed in terms of the engineering material constants as: 1 1 𝜗𝜗12 𝜗𝜗31 1 𝑆𝑆11 = , 𝑆𝑆33 = , 𝑆𝑆12 = − , 𝑆𝑆13 = − , 𝑆𝑆44 = , 𝑆𝑆66 = 2(𝑆𝑆11 − 𝑆𝑆12 ) (3) 𝐸𝐸1 𝐸𝐸3 𝐸𝐸1 𝐸𝐸3 𝐺𝐺13 where 𝐸𝐸1 , 𝐸𝐸3 , 𝐺𝐺13 are transversely isotropic material moduli. 3.1.1.2 Stress-strain relations for an orthotropic material in terms of LCS Transforming the constitutive relations (Equation 1) from the printing material coordinate system (MCS) 1, 2, 3 into the beam local coordinate system (LCS) x, y, z relies on the direction cosines in Table 3.1. Table 3.1: Direction cosines between MCS and LCS. 1 2 3 x 𝑙𝑙1 𝑚𝑚1 𝑛𝑛1 y 𝑙𝑙2 𝑚𝑚2 𝑛𝑛2 z 𝑙𝑙3 𝑚𝑚3 𝑛𝑛3 In Table 3.1, the term 𝑙𝑙1 , 𝑙𝑙2 , 𝑙𝑙3 , 𝑚𝑚1 , 𝑚𝑚2 , 𝑚𝑚3 , 𝑛𝑛1 , 𝑛𝑛2 , 𝑛𝑛3 denotes the cosines of the angles between axes in the two coordinate systems. For example: 𝑚𝑚3 =cos (2, z). The relationship between the stress tensor σ′ in the beam LCS and stress tensor σ in the MCS is: {σ′ } = [A] ∗ {σ} ∗ [A]T (4) 12 where σ′ and σ are expressed in matrix form and A is the rotation matrix generated from the transformation Table 3.1: 𝑙𝑙1 𝑚𝑚1 𝑛𝑛1 [A] = �𝑙𝑙2 𝑚𝑚2 𝑛𝑛2 � (5) 𝑙𝑙3 𝑚𝑚3 𝑛𝑛3 Thus, (4 ) can be written as: 𝜎𝜎x 𝜏𝜏xy 𝜏𝜏xz 𝑙𝑙1 𝑚𝑚1 𝑛𝑛1 𝜎𝜎1 𝜏𝜏12 𝜏𝜏13 𝑙𝑙1 𝑙𝑙2 𝑙𝑙3 �𝜏𝜏yx 𝜎𝜎y 𝜏𝜏yz � = �𝑙𝑙2 𝑚𝑚2 𝑛𝑛2 � �𝜏𝜏21 𝜎𝜎2 𝜏𝜏23 � �𝑚𝑚1 𝑚𝑚2 𝑚𝑚3 � (6) 𝜏𝜏zx 𝜏𝜏zy 𝜎𝜎z 𝑙𝑙3 𝑚𝑚3 𝑛𝑛3 𝜏𝜏31 𝜏𝜏32 𝜎𝜎3 𝑛𝑛1 𝑛𝑛2 𝑛𝑛3 The relationship between the stress in the beam LCS and the corresponding stress state in MCS can be simplified as: 𝜎𝜎x 𝜎𝜎1 ⎧ 𝜎𝜎y ⎫ ⎧ 𝜎𝜎2 ⎫ ⎪ 𝜎𝜎 ⎪ ⎪ 𝜎𝜎 ⎪ z −1 3 𝜏𝜏yz = [𝑇𝑇] 𝜏𝜏23 (7) ⎨ ⎬ ⎨ ⎬ ⎪𝜏𝜏zx ⎪ ⎪𝜏𝜏31 ⎪ ⎩𝜏𝜏xy ⎭ ⎩𝜏𝜏12 ⎭ where T is the transformation matrix: 2 ⎡ 𝑙𝑙1 2 𝑚𝑚1 2 𝑛𝑛1 2 2𝑚𝑚1 𝑛𝑛1 2𝑙𝑙1 𝑛𝑛1 2𝑙𝑙1 𝑚𝑚1 ⎤ ⎢ 𝑙𝑙2 𝑚𝑚2 2 𝑛𝑛2 2 2𝑚𝑚2 𝑛𝑛2 2𝑙𝑙2 𝑛𝑛2 2𝑙𝑙2 𝑚𝑚2 ⎥ ⎢ 2 𝑚𝑚3 2 n3 2 2𝑙𝑙3 𝑚𝑚3 ⎥ [𝑇𝑇]−1 = ⎢ 𝑙𝑙3 2𝑚𝑚3 𝑛𝑛3 2𝑙𝑙3 𝑛𝑛3 ⎥ (8) ⎢𝑙𝑙2 𝑙𝑙3 𝑚𝑚2 𝑚𝑚3 𝑛𝑛2 𝑛𝑛3 𝑚𝑚3 𝑛𝑛2 + 𝑚𝑚2 𝑛𝑛3 𝑙𝑙3 𝑛𝑛2 + 𝑙𝑙2 𝑛𝑛3 𝑙𝑙3 𝑚𝑚2 + 𝑙𝑙2 𝑚𝑚3 ⎥ ⎢𝑙𝑙1 𝑙𝑙3 𝑚𝑚1 𝑚𝑚3 𝑛𝑛1 𝑛𝑛3 𝑚𝑚3 𝑛𝑛1 + 𝑚𝑚1 𝑛𝑛3 𝑙𝑙3 𝑛𝑛1 + 𝑙𝑙1 𝑛𝑛3 𝑙𝑙3 𝑚𝑚1 + 𝑙𝑙1 𝑚𝑚3 ⎥ ⎣𝑙𝑙1 𝑙𝑙2 𝑚𝑚1 𝑚𝑚2 𝑛𝑛1 𝑛𝑛2 𝑚𝑚2 𝑛𝑛1 + 𝑚𝑚1 𝑛𝑛2 𝑙𝑙2 𝑛𝑛1 + 𝑙𝑙1 𝑛𝑛2 𝑙𝑙2 𝑚𝑚1 + 𝑙𝑙1 𝑚𝑚2 ⎦ For small strains undergoing elastic deformation, the transformation equation is similar in form to the stress transformation equation: 𝜀𝜀x 𝜀𝜀1 ⎧ 𝜀𝜀y ⎫ ⎧ 𝜀𝜀2 ⎫ ⎪ 𝜀𝜀z ⎪ ⎪ 𝜀𝜀3 ⎪ ⎪ ⎪ ⎪ ⎪ 𝛾𝛾yz 𝛾𝛾23 2 = [𝑇𝑇]−1 2 (9) ⎨𝛾𝛾xz ⎬ ⎨𝛾𝛾31 ⎬ ⎪2⎪ ⎪2⎪ ⎪𝛾𝛾xy ⎪ ⎪𝛾𝛾12 ⎪ ⎩2⎭ ⎩2⎭ 13 The transformation between the engineering strain vector and the tensor strain vector could be multiplied by the Reuter matrix R 𝜀𝜀x 𝜀𝜀1 𝜀𝜀x ⎧ 𝜀𝜀y ⎫ 𝜀𝜀1 ⎧ 𝜀𝜀2 ⎫ ⎧ 𝜀𝜀y ⎫ ⎪ 𝜀𝜀z ⎪ ⎧ 𝜀𝜀2 ⎫ ⎪ 𝜀𝜀3 ⎪ ⎪ 𝜀𝜀 ⎪ ⎪ ⎪ ⎪ 𝜀𝜀 ⎪ ⎪ ⎪ z 𝛾𝛾yz 𝛾𝛾23 3 𝛾𝛾yz = [R] 2 , 𝛾𝛾 = [R] 2 (10) ⎨ ⎬ ⎨𝛾𝛾zx ⎬ ⎨ 23 ⎬ ⎨𝛾𝛾31 ⎬ ⎪ ⎪𝛾𝛾zx ⎪𝛾𝛾31 ⎪ ⎪2⎪ ⎪2⎪ ⎩𝛾𝛾xy ⎭ ⎪𝛾𝛾xy ⎪ ⎩𝛾𝛾12 ⎭ ⎪𝛾𝛾12 ⎪ ⎩2⎭ ⎩2⎭ where 1 0 0 0 0 0 ⎡0 1 0 0 0 0⎤ ⎢ ⎥ 0 0 1 0 0 0⎥ [R] = ⎢⁠ (11) ⎢0 0 0 2 0 0⎥ ⎢0 0 0 0 2 0⎥ ⎣0 0 0 0 0 2⎦ According to (1), (2), (7), (9), and (10), the relationship between the stress and engineering strain in the beam local coordinate system (LCS) becomes 𝜀𝜀x 𝜎𝜎x 𝜎𝜎x ⎧ 𝜀𝜀y ⎫ ⎧ 𝜎𝜎y ⎫ ⎧ 𝜎𝜎y ⎫ ⎪ 𝜀𝜀 ⎪ ⎪ 𝜎𝜎 ⎪ ⎪ 𝜎𝜎 ⎪ z −1 −1 z ̅ z 𝛾𝛾yz = [R][T] [R] [S][T] 𝜏𝜏yz = [𝑆𝑆] 𝜏𝜏yz (12) ⎨ ⎬ ⎨ ⎬ ⎨ ⎬ ⎪𝛾𝛾zx ⎪ ⎪𝜏𝜏zx ⎪ ⎪𝜏𝜏zx ⎪ ⎩𝛾𝛾xy ⎭ ⎩𝜏𝜏xy ⎭ ⎩𝜏𝜏xy ⎭ Expanding equation (12): ���� 𝑆𝑆 ���� 𝑆𝑆12 ���� 𝑆𝑆13 ���� 𝑆𝑆14 ���� 𝑆𝑆15 ���� 𝑆𝑆16 𝜎𝜎 𝜀𝜀x ⎡ 11 ⎤ x ⎧ 𝜀𝜀y ⎫ ⎢���� 𝑆𝑆21 ���� 𝑆𝑆 22 ���� 𝑆𝑆23 ���� 𝑆𝑆24 ���� 𝑆𝑆25 ���� 𝑆𝑆26 ⎥ ⎧ 𝜎𝜎y ⎫ ⎪ 𝜀𝜀 ⎪ ⎢���� ���� ���� ���� ���� ���� ⎪ ⎪ z 𝑆𝑆31 𝑆𝑆32 𝑆𝑆33 𝑆𝑆34 𝑆𝑆35 𝑆𝑆36 ⎥ 𝜎𝜎z 𝛾𝛾yz = ⎢���� ���� ���� ���� ���� ���� ⎥ 𝜏𝜏yz (13) ⎨ ⎬ ⎢𝑆𝑆41 𝑆𝑆42 𝑆𝑆43 𝑆𝑆44 𝑆𝑆45 𝑆𝑆46 ⎨ ⎥ ⎪𝜏𝜏zx ⎬ ⎪𝛾𝛾zx ⎪ ���� 𝑆𝑆51 ���� 𝑆𝑆52 ���� 𝑆𝑆53 ���� 𝑆𝑆54 ���� 𝑆𝑆55 ���� 𝑆𝑆56 ⎥ ⎪ ⎩𝛾𝛾xy ⎭ ⎢���� ���� ���� ���� ���� ���� ⎩ 𝜏𝜏xy ⎭ ⎣𝑆𝑆61 𝑆𝑆62 𝑆𝑆63 𝑆𝑆64 𝑆𝑆65 𝑆𝑆66 ⎦ where [𝑆𝑆̅] is used to compute the effective material properties in the local coordinate system (LCS) for the beam finite element model. 14 3.1.2 Finite Element Analysis Each lattice member is represented as a single 1D beam finite element based on classical (Euler-Bernoulli) beam theory [65]. A finite element model is developed in Abaqus [66] using the B33 element type, which uses a cubic interpolation of the transverse displacement components. We assume the stress of the beam caused by tension, bending, and torsion are decoupled. 3.1.3 Beam Strain and Stress To predict failure in the orthotropic material within each lattice member, the stress state in the material coordinate system is required. Starting with the displacement and rotation fields in the global coordinates, the following procedure is used. This process is programmed in MATLAB [67]. Step 1: Extract the nodal displacement and rotation fields from the Abaqus simulation results and calculate the strain components at various section points on the circumference of each element node. Note that the beam strains are exact throughout each element when external loads are applied only at lattice joints, neglecting the effects of the joint geometry. Strains are calculated at nodes of beam elements, where they are the largest. Since the joint geometry is not considered, these strains are considered to be nominal but sufficiently accurate for the current purpose of the overall design. The nodal displacements and rotations of the beam are originally in the global coordinate system (GCS), and these must be converted to the beam local coordinate system (LCS). Each beam element has 2 nodes, with 6 degrees of freedom at each node (3 translations and 3 15 rotations). The vector of element displacements and rotations in the global coordinate system (GCS) is: 𝑢𝑢𝑋𝑋𝑋𝑋 ⎧𝑢𝑢𝑌𝑌𝑌𝑌 ⎫ ⎪𝑢𝑢𝑍𝑍𝑍𝑍 ⎪ ⎪ 𝜃𝜃 ⎪ ⎪ 𝑋𝑋𝑋𝑋 ⎪ 𝜃𝜃 ⎪ 𝑌𝑌𝑌𝑌 ⎪ 𝜃𝜃 {𝑈𝑈} = 𝑢𝑢𝑍𝑍𝑍𝑍 (14) ⎨ 𝑋𝑋𝑋𝑋 ⎬ ⎪𝑢𝑢𝑌𝑌𝑌𝑌 ⎪ ⎪𝑢𝑢𝑍𝑍𝐵𝐵 ⎪ ⎪𝜃𝜃𝑋𝑋𝑋𝑋 ⎪ ⎪ 𝜃𝜃𝑌𝑌𝑌𝑌 ⎪ ⎩ 𝜃𝜃𝑍𝑍𝑍𝑍 ⎭ where the displacement components in the GCS of node A are (𝑢𝑢𝑋𝑋𝑋𝑋 , 𝑢𝑢𝑌𝑌𝑌𝑌 , 𝑢𝑢𝑍𝑍𝑍𝑍 ), the displacement components in the GCS of node B are (𝑢𝑢𝑋𝑋𝑋𝑋 , 𝑢𝑢𝑌𝑌𝑌𝑌 , 𝑢𝑢𝑍𝑍𝑍𝑍 ), the rotation components in the GCS for node A are ( 𝜃𝜃𝑋𝑋𝑋𝑋 , 𝜃𝜃𝑌𝑌𝑌𝑌 , 𝜃𝜃𝑍𝑍𝑍𝑍 ) and the rotation components in the GCS for node B are (𝜃𝜃𝑋𝑋𝑋𝑋 , 𝜃𝜃𝑌𝑌𝑌𝑌 , 𝜃𝜃𝑍𝑍𝑍𝑍 ). The length of the beam element is: L = √𝑑𝑑𝑑𝑑 2 + 𝑑𝑑𝑑𝑑 2 + 𝑑𝑑𝑑𝑑 2 (15) 𝑑𝑑𝑑𝑑 𝑑𝑑𝑑𝑑 𝑑𝑑𝑑𝑑 where 𝑑𝑑𝑑𝑑 = 𝑋𝑋𝐵𝐵 − 𝑋𝑋𝐴𝐴 , 𝑑𝑑𝑑𝑑 = 𝑌𝑌𝐵𝐵 − 𝑌𝑌, 𝑑𝑑𝑑𝑑 = 𝑍𝑍𝐵𝐵 − 𝑍𝑍𝐴𝐴 . Let l = ,m= ,n= . 𝐿𝐿 𝐿𝐿 𝐿𝐿 When D = √𝑙𝑙2 + 𝑚𝑚2 > 0, we define the pointwise rotation matrix r as: 𝑙𝑙 𝑚𝑚 𝑛𝑛 −𝑚𝑚 𝑙𝑙 [𝑟𝑟] = � 0� (16.1) 𝐷𝐷 𝐷𝐷 −𝑙𝑙∗𝑛𝑛 −𝑚𝑚∗𝑛𝑛 𝐷𝐷 𝐷𝐷 𝐷𝐷 When D = √𝑙𝑙2 + 𝑚𝑚2 = 0, the rotation matrix r is: 0 0 𝑠𝑠𝑠𝑠𝑠𝑠𝑠𝑠(𝑛𝑛) [𝑟𝑟] = �−1 0 0 � (16.2) 0 −𝑠𝑠𝑠𝑠𝑠𝑠𝑠𝑠(𝑛𝑛) 0 Assembling the 12x12 rotation matrix for each beam element as: 16 𝑟𝑟 0 0 0 [𝑅𝑅] = � 0 𝑟𝑟 0 0� (16.3) 0 0 𝑟𝑟 0 0 0 0 𝑟𝑟 the displacement vector in the local coordinate system {𝑢𝑢} can be defined in terms of the global displacement vector {𝑈𝑈}: {𝑢𝑢} = [𝑅𝑅]{𝑈𝑈} (17) For the beam model, the “centroidal axis” coincides with the axis where the strain due to bending is zero, and the strains caused by tension F, torsion T, and orthogonal bending moments 𝑀𝑀𝑀𝑀 and 𝑀𝑀𝑀𝑀 are decoupled (shear due to bending is neglected). These loading conditions are illustrated in Figure 3.2 (loads are only applied at nodes). Figure 3.2: The loadings on the beam, where M is bending moment, F is axial load, and T is torsion. The strain components are calculated for bending in both XY and XZ planes. The calculation steps are similar for each plane; thus, we only demonstrate the steps of calculating strains for bending in the XY plane. For the Euler Bernoulli beam element, the shape functions are cubic polynomials: 𝐿𝐿3 −3𝐿𝐿𝑥𝑥 2 +2𝑥𝑥 3 𝐿𝐿3 𝑥𝑥−2𝐿𝐿2 𝑥𝑥 2 +𝐿𝐿𝑥𝑥 3 3𝐿𝐿𝑥𝑥 2 −2𝑥𝑥 3 −𝐿𝐿2 𝑥𝑥 2 +𝐿𝐿𝑥𝑥 3 𝑁𝑁1 = 𝑁𝑁2 = 𝑁𝑁3 = 𝑁𝑁4 = (18) 𝐿𝐿3 𝐿𝐿3 𝐿𝐿3 𝐿𝐿3 From the shape functions, the strain-displacement matrix for bending in the XY plane is: 17 1 {𝐵𝐵1 }𝑇𝑇 = [−6𝐿𝐿 + 12𝑥𝑥, −4𝐿𝐿2 + 6𝐿𝐿𝐿𝐿, 6𝐿𝐿 − 12𝑥𝑥, −2𝐿𝐿2 + 6𝐿𝐿𝐿𝐿] 𝑥𝑥 = [0, 𝐿𝐿] (19) 𝐿𝐿3 The strain displacement matrix for the axial load is: 1 {𝐵𝐵2 }𝑇𝑇 = [−1,1] (20) 𝐿𝐿 The strain displacement matrix for torsion is: 1 {𝐵𝐵3 }𝑇𝑇 = [−1,1] (21) 𝐿𝐿 Figure 3.3: The displacement and rotation components in the XY plane. From Figure 3.2, the axial strain at the beam element node is the combination of strains due to bending in the XY plane and axial load. Using equations (20) and (21), the axial strain is at the top and bottom positions on the y-axis are: 𝜀𝜀𝑥𝑥𝑥𝑥 = ±𝑦𝑦 {𝐵𝐵1 }𝑇𝑇 {𝑢𝑢1 } + {𝐵𝐵2 }𝑇𝑇 {𝑢𝑢2 } {𝑢𝑢1 } = {𝑢𝑢𝑌𝑌𝑌𝑌 𝜃𝜃𝑍𝑍𝑍𝑍 𝑢𝑢𝑌𝑌𝑌𝑌 𝜃𝜃𝑍𝑍𝑍𝑍 }𝑇𝑇 (22) {𝑢𝑢2 } = {𝑢𝑢𝑋𝑋𝑋𝑋 𝑢𝑢𝑋𝑋𝑋𝑋 }𝑇𝑇 where y=r (-r) at the top (bottom) of the section. From Figure 3.3, the shear strain at the beam element node results from the torsion. With Equation (22), the shear strain component is: 𝜀𝜀𝑥𝑥𝑥𝑥 = ±𝑟𝑟 {𝐵𝐵3 }𝑇𝑇 {𝑢𝑢3 } (23) {𝑢𝑢3 } = {𝜃𝜃𝑋𝑋𝑋𝑋 𝜃𝜃𝑋𝑋𝑋𝑋 }𝑇𝑇 18 where +r is the counterclockwise rotation around the x-axis, -r is the clockwise rotation around the x-axis. From the equations above, we can calculate the 𝜀𝜀𝑥𝑥𝑥𝑥 strain components at the beam ends for bending in the XY plane. Similar computations are performed for bending in the XZ plane. After we obtain all strain components for bending in both planes, along with torsion and pure axial deformation, we inspect the strain (and stress) state at 36 section points distributed on the beam cross-section (Figure 3.2), equally spaced by 10° . Step 2: Stress in the Material Coordinate System Failure criteria are expressed and evaluated in the MCS, so it is necessary to compute the stress state at each point of interest in the MCS. This can be achieved in the following method. Compute the strain in the beam LCS, compute the stress state in the beam LCS, then transform the stress state into the MCS. The flowchart in Figure 3.4 illustrates the stress calculation process. Figure 3.4: The process to calculate stress in the MCS. The strain state in the beam LCS is calculated using (13). To calculate the strain in the MCS, we start from the plane stress assumption used in classical beam theory: 𝜎𝜎y = 𝜎𝜎z = 𝜏𝜏yz = 0 (24) Thus, at the section points where the strain components were computed formerly, the stress in the local coordinate system (LCS) becomes 𝜎𝜎x 𝜎𝜎x ⎧ 𝜎𝜎y ⎫ ⎧ 0 ⎫ ⎪ 𝜎𝜎 ⎪ ⎪ 0 ⎪ z {𝜎𝜎𝐿𝐿𝐿𝐿𝐿𝐿 } = 𝜏𝜏yz = 0 (25) ⎨ ⎬ ⎨ ⎬ ⎪𝜏𝜏zx ⎪ ⎪𝜏𝜏zx ⎪ ⎩𝜏𝜏xy ⎭ ⎩𝜏𝜏xy ⎭ 19 Then (13) turns to ���� 𝑆𝑆 ���� 𝑆𝑆12 ���� 𝑆𝑆13 ���� 𝑆𝑆14 ���� 𝑆𝑆15 ���� 𝑆𝑆16 𝜎𝜎x 𝜀𝜀x ⎡ 11 ⎤ ⎧ 𝜀𝜀y ⎫ ⎢���� 𝑆𝑆21 ���� 𝑆𝑆22 ���� 𝑆𝑆23 ���� 𝑆𝑆24 ���� 𝑆𝑆25 ���� 𝑆𝑆26 ⎥ ⎧ 0 ⎫ ⎪ 𝜀𝜀 ⎪ ⎢���� ���� ���� ���� ���� ���� ⎪0⎪ z 𝑆𝑆31 𝑆𝑆32 𝑆𝑆33 𝑆𝑆34 𝑆𝑆35 𝑆𝑆36 ⎥ 𝛾𝛾yz = ⎢���� ���� ���� ���� ���� ���� ⎥ 0 (26) ⎨ ⎬ ⎢𝑆𝑆41 𝑆𝑆42 𝑆𝑆43 𝑆𝑆44 𝑆𝑆45 𝑆𝑆46 ⎨ ⎥ ⎪𝜏𝜏zx ⎬ ⎪𝛾𝛾zx ⎪ ���� 𝑆𝑆51 ���� 𝑆𝑆52 ���� 𝑆𝑆53 ���� 𝑆𝑆54 ���� 𝑆𝑆55 ���� 𝑆𝑆56 ⎥ ⎪ ⎩𝛾𝛾xy ⎭ ⎢���� ���� ���� ���� ���� ���� ⎩ 𝜏𝜏xy ⎭ ⎣𝑆𝑆61 𝑆𝑆62 𝑆𝑆63 𝑆𝑆64 𝑆𝑆65 𝑆𝑆66 ⎦ Since the effective material properties were used within the 1D beam element, defining effective material properties 1 1 1 1 1 1 𝐸𝐸𝑥𝑥 = ����� 𝐸𝐸𝑦𝑦 = ����� 𝐸𝐸𝑧𝑧 = ����� 𝐺𝐺𝑦𝑦𝑦𝑦 = ����� 𝐺𝐺𝑥𝑥𝑥𝑥 = ����� 𝐺𝐺𝑥𝑥𝑥𝑥 = ����� (27) 𝑆𝑆11 𝑆𝑆22 𝑆𝑆33 𝑆𝑆44 𝑆𝑆55 𝑆𝑆66 1 1 1 𝜗𝜗𝑦𝑦𝑦𝑦 = −𝐸𝐸𝑦𝑦 ����� 𝜗𝜗𝑥𝑥𝑥𝑥 = −𝐸𝐸𝑥𝑥 ����� 𝜗𝜗𝑥𝑥𝑥𝑥 = −𝐸𝐸𝑥𝑥 ����� 𝑆𝑆32 𝑆𝑆32 𝑆𝑆12 the effective Young’s modulus and poison ratio are used for calculating the stress in the local coordinate system (LCS). 𝜎𝜎x = 𝐸𝐸𝑥𝑥 ∗ 𝜀𝜀x (28.1) 𝐸𝐸𝑥𝑥 𝜏𝜏zx = ∗ 𝛾𝛾zx (28.2) 2(1+𝜗𝜗𝑥𝑥𝑥𝑥 ) 𝐸𝐸𝑥𝑥 𝜏𝜏xy = ∗ 𝛾𝛾xy (28.3) 2(1+𝜗𝜗𝑥𝑥𝑥𝑥 ) Thus, all strain and stress components in the beam local coordinate system (LCS) are known. For the stress in the material coordinate system (MCS), with (2), we have: {𝜎𝜎𝑀𝑀𝑀𝑀𝑀𝑀 } = [𝑆𝑆]−1 {𝜀𝜀𝑀𝑀𝑀𝑀𝑀𝑀 } (29) Then the stress state can be transformed into the MCS by (7) to obtain: {𝜎𝜎𝑀𝑀𝑀𝑀𝑀𝑀 } = [𝑇𝑇]{𝜎𝜎𝐿𝐿𝐿𝐿𝐿𝐿 } (31) 3.1.4 Failure Criteria Within the optimization study, three failure criteria are used as constraints: 20 𝑀𝑀𝑀𝑀𝑀𝑀 𝑠𝑠𝑠𝑠𝑠𝑠𝑠𝑠𝑠𝑠𝑠𝑠 𝑖𝑖𝑖𝑖 𝑝𝑝𝑝𝑝𝑝𝑝𝑝𝑝𝑝𝑝𝑝𝑝𝑝𝑝𝑝𝑝 𝑝𝑝𝑝𝑝𝑝𝑝𝑝𝑝𝑝𝑝 𝜎𝜎1 𝐹𝐹𝐹𝐹𝐹𝐹𝐹𝐹𝐹𝐹𝐹𝐹𝐹𝐹 𝑖𝑖𝑖𝑖𝑖𝑖𝑖𝑖𝑖𝑖 𝛼𝛼1 = = (32.1) 𝑆𝑆𝑆𝑆𝑆𝑆𝑆𝑆𝑆𝑆𝑆𝑆𝑆𝑆ℎ 𝑖𝑖𝑖𝑖 𝑝𝑝𝑝𝑝𝑝𝑝𝑝𝑝𝑝𝑝𝑝𝑝𝑝𝑝𝑝𝑝 𝑝𝑝𝑝𝑝𝑝𝑝𝑝𝑝𝑝𝑝 𝑆𝑆1 𝑀𝑀𝑀𝑀𝑀𝑀 𝑠𝑠𝑠𝑠𝑠𝑠𝑠𝑠𝑠𝑠𝑠𝑠 𝑣𝑣𝑣𝑣𝑣𝑣𝑣𝑣𝑣𝑣𝑣𝑣𝑣𝑣𝑣𝑣 𝑡𝑡𝑡𝑡 𝑝𝑝𝑝𝑝𝑝𝑝𝑝𝑝𝑝𝑝𝑝𝑝𝑝𝑝𝑝𝑝 𝑝𝑝𝑝𝑝𝑝𝑝𝑝𝑝𝑝𝑝 𝜎𝜎3 𝐹𝐹𝐹𝐹𝐹𝐹𝐹𝐹𝐹𝐹𝐹𝐹𝐹𝐹 𝑖𝑖𝑖𝑖𝑖𝑖𝑖𝑖𝑖𝑖 𝛼𝛼2 = = (32.2) 𝑆𝑆𝑆𝑆𝑆𝑆𝑆𝑆𝑆𝑆𝑆𝑆𝑆𝑆ℎ 𝑣𝑣𝑣𝑣𝑣𝑣𝑣𝑣𝑣𝑣𝑣𝑣𝑣𝑣𝑣𝑣 𝑡𝑡𝑡𝑡 𝑝𝑝𝑝𝑝𝑝𝑝𝑝𝑝𝑝𝑝𝑝𝑝𝑝𝑝𝑝𝑝 𝑝𝑝𝑝𝑝𝑝𝑝𝑝𝑝𝑝𝑝 𝑆𝑆3 𝑀𝑀𝑀𝑀𝑀𝑀 𝑠𝑠ℎ𝑒𝑒𝑒𝑒𝑒𝑒 𝑠𝑠𝑠𝑠𝑠𝑠𝑠𝑠𝑠𝑠𝑠𝑠 𝑖𝑖𝑖𝑖 𝑝𝑝𝑝𝑝𝑝𝑝𝑝𝑝𝑝𝑝𝑝𝑝𝑝𝑝𝑝𝑝 𝑝𝑝𝑝𝑝𝑝𝑝𝑝𝑝𝑝𝑝 𝜎𝜎13 𝐹𝐹𝐹𝐹𝐹𝐹𝐹𝐹𝐹𝐹𝐹𝐹𝐹𝐹 𝑖𝑖𝑖𝑖𝑖𝑖𝑖𝑖𝑖𝑖 𝛼𝛼3 = = (32.3) 𝑆𝑆ℎ𝑒𝑒𝑒𝑒𝑒𝑒 𝑠𝑠𝑠𝑠𝑠𝑠𝑠𝑠𝑠𝑠𝑠𝑠𝑠𝑠ℎ 𝑖𝑖𝑖𝑖 𝑝𝑝𝑝𝑝𝑝𝑝𝑛𝑛𝑡𝑡𝑡𝑡𝑡𝑡𝑡𝑡 𝑝𝑝𝑝𝑝𝑝𝑝𝑝𝑝𝑝𝑝 𝑆𝑆13 For failure index 𝛼𝛼1 , the maximum stress is calculated as the principal stress within the printing plane (1st -2nd material plane). For the failure index 𝛼𝛼2 , the stress is calculated in the plane vertical to the printing plane (the 3rd material plane). The failure index 𝛼𝛼3 considers the shear stress between each layer of material deposited. 3.1.5 Material Properties For the study, we used the ABS material of the FDM process as an example. The orthotropic material properties of ABE material listed in Table 3.2 [21], [68] were applied in the subsequent studies. To test the reliability of the optimized lattice structure based on isotropic material properties, the isotropic material properties of ABS material were also listed in Table 3.2 for comparison studies. 21 Table 3.2: Material properties for ABS Orthotropic material properties Strength Young’s modulus (MPa) Compression (MPa) Tension (MPa) 𝐸𝐸1 2,078 28.83 25.51 𝐸𝐸3 16,48 29.48 14.35 𝐺𝐺12 E1/2/(1+v12) 𝐺𝐺13 686 𝐺𝐺23 686 𝛾𝛾12 0.3 𝛾𝛾13 0.315 𝛾𝛾23 0.315 Isotropic material properties Young’s modulus (MPa) Strength (MPa) E 2200 31 3.2 Model Validation The beam model formulation and implementation were validated by comparing the model results with those of a 3D solid beam with orthotropic material properties. For this single-beam validation model, the beam radius was set to 5mm, and the length was 100mm. For both models, one end of the beam was clamped, and the forces in x and y directions 𝐹𝐹𝑥𝑥 and 𝐹𝐹𝑦𝑦 , and torsion in the x-axis 𝑇𝑇𝑥𝑥 were applied at the opposite end of the beam. Figures 3.5 and 3.6 illustrate the models for the 1D beam element and the 3D solid beam. 22 Figure 3.5: 1D beam element model used for validation. Figure 3.6: 3D solid beam model used for validation. For deflection and stress components of the end of the lattice beam, the results were tested for different manufacturing orientations within the first quadrant, as shown in Figure 3.7. Recall that the principal MCS is aligned with the printing direction. The printing direction is defined by 2 angles in the spherical coordinate system (𝜃𝜃, 𝜑𝜑), where 𝜑𝜑 is the angle between the printing direction and the XY plane of the global coordinate system, and the 𝜃𝜃 is the angle between the reflection of the printing direction on the XY plane and the X coordinate of the global coordinate 𝜋𝜋 𝜋𝜋 system. The range of printing angles considered was: 𝜑𝜑 = [− , ], 𝜃𝜃 = [0,2𝜋𝜋]. 2 2 The ranges of the print direction angles were discretized in 15-degree increments for both φ and θ. For each build direction (set of angles) the ratio of results from the two models (1D and 3D) was calculated. A ratio near 1 indicated good agreement between the results. 23 Figure 3.7: Definition of the angles that define the print direction. 𝐸𝐸1 The validation study was performed for three different E1 to E3 ratios: ≈ 1, 10, and 100. 𝐸𝐸3 In all three cases, the ratios of deflection and stress components in the principal MCS were 𝐸𝐸1 compared favorably. For brevity, only the results for ≈ 10 (E1=2,078 MPa, E3=210 MPa) are 𝐸𝐸3 shown here in Figures 3.8 and 3.9. It is well-known that stress components having a small contribution to the overall strain energy may not be accurately predicted in a finite element model, where the objective is the minimization of total potential energy. Thus, when the value of a stress component was less than 1% of the maximum stress value among all other stress components for that print direction, the stress ratio results for these small stress components were not included in the plots. At some locations, the values of some stress components are very small compared to the other components. While the difference between the two computed values is also very small, the ratio of these two small values can be large, falsely implying a significant error. We have chosen to neglect these components where this numerical issue occurs. 24 Figure 3.8: Ratio of beam deflections predicted by 1D and 3D orthotropic beam models. Figure 3.9: Ratio of beam stress components (S11, S22, S33, S12) predicted by 1D and 3D orthotropic beam models. Components S13 and S23 were less than 1% of the maximum stress value at all points, so these plots were omitted. 25 3.3 Studies and Results In an additively manufactured lattice structure, the principal directions of material orthotropy are related to the build direction, and the orientation of lattice members relative to those principal material directions can have a significant effect on the structural performance of each lattice member as well as the overall structure. In other words, the design of each lattice member (e.g., its minimum radius) depends on the build direction. For this reason, it is necessary to simultaneously optimize the manufacturing build direction of the overall structure and the radius of each lattice member. In the current study, we consider the build (or printing) direction as well as the radius of each lattice member as design variables to be optimized. The goal of the study is to generate the lightest lattice structure (minimize volume) that also satisfies constraints on strength under a prescribed loading condition. Due to the strong coupling of beam radii to a particular build direction, a bi-level optimization strategy is used here to increase the efficiency of the design exploration. In the 1st optimization level (the upper level), the objective is to minimize the material volume of the lattice structure by modifying the printing direction. In the 2nd optimization level (the lower or inner level), the objective is to minimize the material volume of the lattice structure for a given printing direction by modifying the beam radii while accounting for material orthotropy in stiffness and strength. The optimization process is shown in the flowchart in Figure 3.10. 26 Figure 3.10: Flowchart of the bi-level optimization process. Procedures within the solid line belong to the 2nd (lower) level optimization, while procedures within the dotted line belong to the 1st (upper) level optimization. The software package(s) used for each procedure are listed beside the blocks. The 1st (upper) level optimization defines a proposed printing direction by assigning values to θ, φ. Then, for that printing direction, optimization is performed on the radii variables to minimize the material volume within the 2nd (lower) level, while satisfying constraints on orthotropic material strength. The failure criteria are defined in terms of failure indices 𝛼𝛼𝑖𝑖 , as defined herein. The bi-level optimization returns the optimized radii and the printing angles for the optimal solutions with minimum volume. The software packages HEEDS MDO, Abaqus, Python, and MATLAB are utilized in various steps within this study. HEEDS MDO [69] is a multi-disciplinary design exploration software package. It is used here to perform design space exploration at both levels and to control the overall process automation. Abaqus[66], [70] is used to perform the structural finite element 27 analysis of each lattice design. Internal code written in the Python [70], [71] scripting language is used to extract simulation results from the Abaqus output (*.odb) file and to send it to MATLAB for post-processing. Internal code written in MATLAB [67] is used to compute orthotropic material properties and effective beam properties for a set of given printing angles and beam orientation, and to post-process the stress state and evaluate failure criteria in the material coordinate system. 3.3.1 Bi-level 1D Beam Optimization To illustrate the effect of material orthotropy on the design of an additively manufactured lattice structure, several optimization studies have been performed. The first one is a simple 1D single beam model, as shown in Figure 3.11. Figure 3.11: 1D beam model clamped at one end and loaded at the opposite end. The beam is fixed (clamped) at one end and loaded at the other end by a torque along the X- axis and forces in the X and Y directions, producing axial, bending, shear, and torsional loading on the beam. At the 1st optimization level (the upper level), the study sweeps along the printing angles 𝜃𝜃 ∈ [0° , 90° ] and 𝜑𝜑 ∈ [0° , 90° ] with an increment of 10° . At the 2nd optimization level (the lower level), the objective is to minimize the material volume of the beam for a given printing direction by modifying the beam radius while accounting for material orthotropy in stiffness and strength. The optimization problem statement for the 1D beam element is: 28 Optimization statement: 1st (upper) level sweep (not optimization): Objective: Minimize the volume of the beam Constraints: Failure indices 𝛼𝛼𝑖𝑖 ≤ 1 (i =1,2,3) Variables: Printing angles (θ, φ) 𝜃𝜃 ∈ [0° , 90° ] and 𝜑𝜑 ∈ [0° , 90° ], the increment is 10° 2nd (lower) level optimization (for a fixed set of printing angles θ, φ): Objective: Minimize the volume of the beam Constraints: Failure indices 𝛼𝛼𝑖𝑖 ≤ 1 (i =1,2,3) Variables: The radius of the beam 𝑟𝑟 ∈ [4, 6], the increment is 0.05mm Results The sweep method is employed to demonstrate the influence of the printing angles on the beam volume. For the 1D beam optimization study, we only choose the first quadrant, in which printing angle 𝜃𝜃 ∈ [0° , 90° ] and 𝜑𝜑 ∈ [0° , 90° ], to display the results. In Figure 3.12, the color of each square as well as its height on the volume axis represent the volume for the given printing angles 𝜃𝜃 and 𝜑𝜑. 29 Figure 3.12: Bi-level 1D beam optimization result. For this 1D beam bi-level optimization study, the maximum volume design occurs at a build direction defined by (θ, φ) = (90, 0) and (θ, φ) = (90, 90), where the volume is 7390 𝑚𝑚𝑚𝑚3 and shown in yellow on the left and right corners in Figure 3.12. The minimum volume design occurs at the build direction (30, 30), where the volume is 5670 𝑚𝑚𝑚𝑚3 and shown at the concave part in Figure 3.12. The difference between designs with maximum volume and minimum volume is about 23.3%. This significant difference and the strong variations in allowable volume shown in Figure 3.12 illustrate the importance of including build direction and material orthotropy when designing an AM lattice. Comparison with the isotropic material lattice To further demonstrate the influence of orthotropic material strengths on the performance of the lattice structure, an optimization study with the isotropic material property was performed. In this isotropic case, the isotropic material properties were used (Table 3.2) and the isotropic material strength was considered as the constraint. Within the same range of the radii, an optimization study was performed to obtain the optimal design for the isotropic material beam. 30 Using the radii from the optimal design using isotropic materials, a study using the sweep method was performed by varying the printing angle (θ, φ). Using the corresponding orthotropic material properties, it was found that for all printing angles, the optimized isotropic design violated the orthotropic material strength constraints. Thus, the optimized beam based on isotropic material properties would likely fail, again illustrating the risk of ignoring the manufacturing orientation-induced orthotropy. 3.3.2 Bi-level BCC Lattice Structure Optimization To illustrate the influence of printing direction and the corresponding material orthotropy on the design of a lattice structure, an optimization study of a BCC lattice structure was performed. A 2 × 1 × 1 BCC lattice structure was considered. It is fixed at one end, and loads were applied on the other end in the negative Y-axis, as seen in Figure 3.13. Figure 3.13: 3D bcc lattice structure. A bi-level optimization strategy is used. In the 1st optimization level (the upper level), the study sweeps along the printing angles 𝜃𝜃 ∈ [0° , 90° ] and 𝜑𝜑 ∈ [0° , 90° ] with an increment of 15° to show the printing direction’s effect on the volume of the structure. At the 2nd optimization level (the lower level), the objective is to minimize the material volume for the given printing 31 direction by modifying the members' radii, while accounting for material orthotropy in stiffness and strength. The optimization statement is shown below. Optimization statement: 1st (upper) level sweep (not optimization): Variables: Printing angles (θ, φ) 𝜃𝜃 ∈ [0° , 90° ] and 𝜑𝜑 ∈ [0° , 90° ], the increment is 15° 2nd (lower) level optimization (for a fixed set of printing angles θ, φ): Objective: Minimize the material volume of the BCC lattice structure Constraints: Failure indices 𝛼𝛼𝑖𝑖 ≤ 1 (i =1,2,3) Variables: Radii of lattice members 𝑟𝑟𝑖𝑖 (i=1 to number of lattice members) 𝑟𝑟𝑖𝑖 ∈ [0.05,0.1], the increment is 0.01mm Results The sweep method was used for the 1st (upper) level of the BCC lattice structure optimization study with an increment of 15° for each printing angle, and the first quadrant with printing angle 𝜃𝜃 ∈ [0° , 90° ] and 𝜑𝜑 ∈ [0° , 90° ] was chosen to demonstrate the volume for each printing angle. 32 Figure 3.14: Bilevel BCC lattice structure optimization result. From the sweeping study of the printing angle for the BCC lattice structure (Figure 3.14.), the maximum volume of the lattice structure occurs at the printing angle (θ, φ) = (90° , 15° ), which is 0.34398 𝑚𝑚𝑚𝑚3 . The minimum volume of 0.2788 𝑚𝑚𝑚𝑚3 occurs at the printing angle (θ, φ) = (30° , 60° ), (45° , 75° ). Comparing the maximum and the minimum volume of the BCC lattice structure, there is a difference of 18.9%. From the results, the printing angles do have a significant impact on the volume of the AM BCC lattice structure. Comparison with the isotropic material lattice With the aim of comparison and validation, an optimization study using the isotropic material properties was performed. Using the radii from this optimization study, a sweeping study was then performed on the printing angles, using orthotropic material properties. This study demonstrated that the optimal isotropic design was no longer feasible under the orthotropic material strength constraints for all printing angles. This further confirms the necessity for considering manufacturing orientation-induced orthotropy. 33 3.3.3 3D Bracket Lattice Structure Optimization To demonstrate the influence of orthotropic material strength effects on a medium-scale lattice structure, a 3D bracket with 1,369 lattice members was considered. All members were assumed to have the same radius, so there was only one design variable. Figure 3.15: The 3D bracket lattice model. As shown in Figure 3.15, this model was a hanger bracket. Its outer dimensions for were: length in x direction = 300 mm, width in z direction = 200 mm, depth in y = 20 mm. The two mounting holes on the left of the bracket were fixed, and a vertical load (in the z-direction) was applied to the mounting hole on the right side. Problem statement: 1st (upper) level optimization: Objective: Minimize the material volume of the 3D bracket lattice structure Constraints: Failure indices 𝛼𝛼𝑖𝑖 ≤ 1(i =1,2,3) (Values of 𝛼𝛼𝑖𝑖 are obtained from Level 2) Variables: Printing angles (θ, φ) 34 𝜃𝜃 ∈ [0° , 90° ] and 𝜑𝜑 ∈ [0° , 90° ], the increment is 5° 2nd (lower) level optimization (for a fixed set of printing angles θ, φ): Objective: Minimize the material volume of the 3D bracket lattice structure Constraints: Failure indices 𝛼𝛼𝑖𝑖 ≤ 1 (i =1,2,3) Variables: Radii of lattice members 1 ≤ 𝑟𝑟𝑖𝑖 ≤ 3, the increment is 0.05mm Results Multiple solutions were found through the bilevel optimization of this 3D lattice bracket. As seen in the parallel plot for the bi-level optimization study (Figure 3.16), several different printing angles are possible for designs with the same minimum volume. This is illustrated more clearly in Figure 3.17, where only the minimum volume solutions are shown. For the optimal design of the 3D lattice bracket, the minimum volume is 474,400 𝑚𝑚𝑚𝑚3 ; however, the printing angles for the optimal solutions vary, with 𝜑𝜑 ranging from [−35° , 25° ], and 𝜃𝜃 ranging between [ 95° , 120° ]. This allows the designer and manufacturer to select the most convenient or economical printing angle from among the optimal solutions. One of these optimal solutions from the bi-level optimization is shown in Figure 3.18. 35 Figure 3.16: Parallel coordinates plot for bi-level optimization for 3D bracket lattice model. Figure 3.17: Parallel coordinate plot for optimal solutions for 3D bracket lattice model. Figure 3.18: Optimized 3D lattice bracket for 3D bracket lattice model. 36 Comparison between the optimized isotropic material lattice hanger and the orthotropic material strength constraints An optimization study was also conducted using the isotropic material properties in Table 3.2. The young’s modulus for isotropic material is higher than 𝐸𝐸1 , the minimum volume for the isotropic material lattice was 394,380 𝑚𝑚𝑚𝑚3 . When this design was evaluated using orthotropic material strength constraints, it violated the constraints for all printing angles. 3.3.4 Bottle Opener Lattice Structure Optimization The final example application is a lattice bottle opener, shown in Figure 3.19. The overall dimensions of the model are 80 mm × 20 mm × 15 mm. The model is symmetric about the mid- plane, but the full structure is represented in the model. For the FEA simulation of this model, boundary conditions representing the levering mechanism of a real bottle opener were used. Region A in the model was acting as the fulcrum, in which the movement was fixed on the cap (𝑈𝑈𝑈𝑈 = 𝑈𝑈𝑈𝑈 = 0, with 𝑈𝑈𝑈𝑈 = 0 at the center). Before the cap is removed, the tooth or lip B remained still (𝑈𝑈𝑈𝑈 = 𝑈𝑈𝑈𝑈 = 0, with 𝑈𝑈𝑈𝑈 = 0 at the center). A distributed force was exerted on the edge of handle C, with a total magnitude of 36N in the y- direction. B 37 Figure 3.19: Boundary conditions for bottle opener lattice structure. Optimization statement: 1st (upper) level optimization: Objective: Minimize the material volume of the bottle opener lattice structure Constraints: Failure indices 𝛼𝛼𝑖𝑖 ≤ 1, i =1,2,3 (Values of 𝛼𝛼𝑖𝑖 are obtained from 2nd level) Variables: Printing angles (θ, φ), 𝜑𝜑 ∈ [−90° , 90° ] the increment is 5° ; 𝜃𝜃 = 90° 2nd (lower) level optimization (for a fixed set of printing angles θ, φ): Objective: Minimize the material volume of the bottle opener lattice structure Constraints: Failure indices 𝛼𝛼𝑖𝑖 ≤ 1 (i =1,2,3) Variables: Radii of lattice members 1 ≤ 𝑟𝑟𝑖𝑖 ≤ 3.6 the increment is 0.1mm Settings for the optimization study: For the 1st (upper) level optimization, due to the symmetry, the printing angle 𝜃𝜃 is fixed at 90° , so only the printing angle 𝜑𝜑 is varied. The number of iterations (the optimization budget) at this level was set to 16. At the 2nd (lower) level optimization, there are 102 lattice beams in total to be optimized. Due to the symmetry of the model, the number of optimization variables was reduced to 51. Based on initial experiments, it was found that 12,000 evaluations were needed to achieve nearly 38 converged results at this level. While it cannot be easily proven, these experiments suggested that the design space is non-convex, likely due to the complexity of the orthotropic material strength constraints. Results Because the optimization approach is stochastic, three sets of optimization studies were conducted, each starting with a different baseline design. These three baseline designs had uniform initial lattice beam radii (volume) of 1.5 mm (9593.3 𝑚𝑚𝑚𝑚3 ), 2.5 mm (26648 𝑚𝑚𝑚𝑚3 ), and 3.5 mm (52230 𝑚𝑚𝑚𝑚3 ). For the three sets of optimization results, the trends for different printing angles 𝜑𝜑 were similar (Figure 3.20). When 𝜑𝜑 = 0° or −30° , the volume of the lattice structure was the highest. The optimal solution was at 𝜑𝜑 = 50° with a minimum volume of 6949.8 𝑚𝑚𝑚𝑚3 . These results illustrate that, when material orthotropy is considered, the printing angle has a significant impact on the volume of the lattice structure. Multi Baseline Bi-level Optimization results 12,000 10,000 8,000 volume 6,000 6949.8 4,000 1.5mm baseline 2,000 2.5mm baseline 3.5mm baeline 0 -90 -75 -60 -45 -30 -15 0 15 30 45 60 75 90 Printing Angle: 𝜑𝜑 Figure 3.20: Bi-level optimization results from 3 baseline designs: 1.5 mm, 2.5 mm, and 3.5 mm. 39 The optimized bottle opener lattice structure is shown in Figure 3.21. The thickest members are located at the lever mechanism, where the stress concentration occurs. The beam members are very thick in that region, so beam theory may not represent the behavior adequately. Nevertheless, the results are useful from an intermediate design stage perspective. Figure 3.22 shows the optimal build orientation for additive manufacturing. Figure 3.21: Optimal lattice structure from 2.5mm baseline design. As shown in Figure 3.20, near the optimal build angle of 𝜑𝜑 = 50° , several printing angles generate good designs, namely at 𝜑𝜑 = 45° , 50° , 55° . To further investigate the results at these printing angles, a convergence study of these points was performed. Figure 3.22: The optimal printing orientation (Z) of the bottle opener lattice. 40 Using each optimized solution as baseline designs for 𝜑𝜑 = 45° , 50° , 55° , more search evaluations were performed within a reduced range of radii containing the optimized designs. Assuming that the design space is convex within this limited range, a gradient-based optimization algorithm, quadratic programming (QP), was applied to search for the local minimum. Results with Limited Radii Range 8200 8000 7800 Volume 7600 1.5mm 7400 2.5mm 7200 3.5mm 7000 6949.8 6800 6,932.30 40 45 50 55 60 Printing Angle: Phi Figure 3.23: The optimal results from the convergence study. As shown in Figure 3.23, after 3,000 additional search evaluations in the region of each optimized design, the minimum volume of all three designs for 𝜑𝜑 = 50° converge to approximately 6932 𝑚𝑚𝑚𝑚3 . For other print angles, the solutions did not change. These results imply that 𝜑𝜑 = 50° is the best choice for printing angle. 41 Chapter 4 Optimizing the Manufacturability and Performance of Additively Manufactured Lattice Structures To account for these geometric limitations associated with an AM process, a manufacturability analysis based on model correction [38] is used here within the design process. The model correction analysis first slices the designed AM geometry into voxels. Then, the critical features affected by the given AM process are identified by considering the resolution or the minimum printable feature of the AM process. Modifications are based on the topology of the original design geometry and the toolpath of the print head. In this chapter, we evaluate the manufacturability of the AM lattice structure with both the material orthotropy and manufacturing resolution limitations. First, the theory based on material orthotropy analysis and manufacturability-oriented model correction is explained. Then, the effect of geometry influence is studied based on the optimization of the example cases. Last, a comprehensive manufacturability analysis with bi-level optimization studies based on geometry effects and material orthotropy is conducted. 4.1 Methodology 4.1.1 Design of As-Manufactured Lattice Structures In this study, the originally designed geometry and material of the lattice structure compose what is called the as-designed lattice. Because AM processes have material and resolution limitations, there may be significant differences between the original as-designed model and the final manufactured part. To account for this, a model correction was introduced into the design process to generate an optimized as-manufactured design. 42 In a previous study [17], additively manufactured lattice structures were optimized while accounting for manufacturing-induced material orthotropy. A generalized three-dimensional orthotropic material model was introduced into each lattice beam, with the principal material coordinates dependent on the manufacturing build direction or other manufacturing features. It was shown that not accounting for this orthotropy in strength and stiffness led to very unconservative designs. To account for as-manufactured geometry, we aim to analyze as-designed lattice structures and make necessary geometric modifications to solve any manufacturability issues for a selected AM approach and process parameters associated with it. We target to alleviate problems such as thin features, sharp corners, and topological changes caused by these issues. We adopt a similar approach to [38] in generating a corrected model. The corrected model accurately approximates the as-manufactured part, and it can be printed using the indented AM process without any failure. This approach allows us to evaluate the manufacturability of each candidate lattice design and measure the minimal amount of change required to make the design manufacturable quantitatively. Given an input 3D lattice model and minimum printable feature size dictated by the selected AM process, we slice the shape in intervals equivalent to the print layer height and rasterize them to generate images representing 2D slices, 𝑆𝑆𝑖𝑖 . Each slice is then processed to construct a meso- skeleton, 𝑀𝑀𝑖𝑖 — the maximal area where a print head can be positioned during the printing process— that is topologically equivalent to the corresponding slice. Dilation of a meso-skeleton with a circular structuring element, 𝐹𝐹 corresponding to the minimum printable feature results in a corrected slice 𝐶𝐶𝑖𝑖 , i.e., 𝐶𝐶𝑖𝑖 = 𝑀𝑀𝑖𝑖 ⨁ 𝐹𝐹, denoting a morphological dilation operation. The union of the corrected slices constitutes the corrected model. 43 In order to compute the meso-skeleton of a slice, we perform a topology preserving thinning operation described in [38]. Starting from a binary image representing a slice 𝑆𝑆𝑖𝑖 , contour pixels that do not contribute to the topology are removed iteratively. Deletion or retention of a contour pixel is determined based on the configuration of the pixel in its local 8-connected neighborhood. The thinning process continues until there is no removable contour pixel left such that (1) removal of any contour pixel changes the topology or (2) shrinks the meso-skeleton beyond the erosion of the slice 𝑆𝑆𝑖𝑖 with 𝐹𝐹 , i.e., 𝑆𝑆𝑖𝑖 ⊖ 𝐹𝐹 . While the former condition guarantees that the topology of the meso-skeleton is identical to the input slice, the latter one prevents unnecessary iterations as the erosion serves as an upper bound for the maximal allowable region that the print head can traverse during the printing process. In model correction, the rasterization resolution and structuring element size play an important role in the accurate estimation of the as-manufactured shape. Suppose the manufacturing resolution is defined by the minimum printable feature size that is represented by a circle of diameter d. We discretize d by 𝑛𝑛𝑑𝑑 pixels to approximate the circular structuring element 𝐹𝐹. Then, the rasterization resolution of slices is selected as 𝑛𝑛𝑆𝑆 = 𝑙𝑙 × 𝑛𝑛𝑑𝑑 / 𝑑𝑑 where 𝑙𝑙 is the largest dimension of the input object bounding box. Expectedly, the selection of a large 𝑛𝑛𝑑𝑑 value results in a more accurate approximation of the as-manufactured shape. Yet, larger 𝑛𝑛𝑑𝑑 increases the computational cost, as well as the memory usage as the number of iterations in the thinning process and the size of each image representing the slices increase simultaneously. In our experiments, we observed that 𝑛𝑛𝑑𝑑 ≥ 5 is sufficient for our purposes. 44 4.2 Studies and Results We first investigated the influence of voxel size on the final volume of the as-manufactured lattice structure. We conducted convergence tests with a different number of voxels for a given AM process resolution. Based on the convergence tests, we optimized the printing direction of AM lattice structure models with different voxel sizes to get the minimum material difference between the as- designed lattice and the as-manufactured lattice. Last, we combined manufacturability-oriented model correction and the AM-induced material orthotropy within the evaluation of potential AM lattice structures to optimize the lattice structure for a specific AM process minimum printable feature size and material orthotropic strengths. In the study, we are ignoring the possible need for support structures in the manufacturing process. 4.2.1 Convergence Test for Different Voxel Sizes To analyze the result of a study based on manufacturability-oriented model correction, there are several parameters needed. Voxel size (measured in mm) is an input parameter that can be adjusted according to the minimum printable feature size of the given AM process. The absolute difference is the absolute number of voxel changes throughout the model correction process. Thus, absolute difference = number of voxels added + number of voxels removed (1) The absolute material change is: absolute material change = absolute difference * 𝑣𝑣𝑣𝑣𝑣𝑣𝑣𝑣𝑣𝑣 𝑠𝑠𝑠𝑠𝑠𝑠𝑠𝑠 3 (2) 45 For the final modified geometry, the total volume of the as-manufactured model is calculated by: total volume = (original voxels + added voxels – removed voxels) * 𝑣𝑣𝑣𝑣𝑣𝑣𝑣𝑣𝑙𝑙 𝑠𝑠𝑠𝑠𝑠𝑠𝑠𝑠 3 (3) To illustrate the influence of voxel size on the final volume of the as-manufactured lattice structure, we optimized a single beam of circular cross-section with different voxel sizes. The length of the lattice beam was 10mm and the radius was 0.25mm. The Selective Laser Melting (SLM) with 0.1mm resolution was selected as the AM process for the model. Figure 4.1: Single lattice beam We conducted optimization studies based on different voxel sizes, try to compare the final volume of the as-manufactured lattice beam between different voxel sizes and printing angles to test what voxel size would be reasonable, and got the results converged. Optimization statement Objective: Minimize the absolute volume difference (voxelized volume) Variables: Rotation angles of the lattice beams (𝜑𝜑, 𝜃𝜃) 𝜑𝜑 ∈ [−90°, 90°], 𝜃𝜃 ∈ [−90°, 90°] 𝜑𝜑, 𝜃𝜃 define the print angle for the SLM process and represent the rotation of the geometry about the z and x-axes, respectively. 46 (a) (b) Figure 4.2: Definition of rotation angles 𝜑𝜑, 𝜃𝜃. Allowing a 5-degree increment of the rotation angle, the optimization results are shown in Table 4.1. With different voxel sizes, the absolute difference and absolute material change vary, but the total volume of the model changes very little, with the difference between the maximum volume and minimum volume being 0.2% (Figure 4.3). Thus, even the largest voxel size considered (4mm) is sufficiently small to represent the geometry for the purposes of this study. Table 4.1: Results with optimization of 5-degree incremental rotation angle Absolute Absolute Voxel size Total volume Material difference (mm) Change 4 0.05015 0.0005 2.2258 336 0.02503 0.00527 2.23057 Same AM 432 0.01667 0.002 2.2273 1 D beam resolution with 656 0.0125 0.00128 2.22658 different voxel sizes 955 0.01 0.00096 2.22626 715 0.00834 0.00041 2.22571 2992 0.00714 0.00109 2.22639 8662 0.00625 0.00211 2.22741 47 Total Volume 2.23057 2.231 2.23 Total Volume 2.229 2.2273 2.22741351 2.228 2.22658 2.22626 2.22639 2.227 2.2258 2.22571 2.226 2.225 2.224 2.223 Voxel Size Figure 4.3:The total volume of optimized beams allowing a 5-degree increment of the rotation angles. The above study was repeated using a 1-degree increment of the rotation angles. In Table 4.2, with different voxel sizes, the total volume is almost the same, the difference between the maximum and minimum volume is less than 0.18% for this study (Figure 4.4). Table 4.2: Results with optimization of 1-degree incremental rotation angle Absolute Absolute Total volume Voxel size Material difference Change 4 0.05014 0.0005 2.2258 294 0.02503 0.00461 2.22991 1D beam (1 degree) Same AM 394 0.01667 0.00183 2.22713 resolution with 546 0.0125 0.00107 2.22637 different voxel sizes 711 0.01 0.00071 2.22601 398 0.00834 0.00023 2.22553 2829 0.00715 0.00103 2.22633 7635 0.00625 0.00186 2.22716 48 Total volume 2.24 2.235 2.22991 Total Volume 2.23 2.2258 2.22713 2.22637 2.22601 2.22633 2.22716 2.22553 2.225 2.22 2.215 2.21 voxel size Figure 4.4: Total volume of optimized beams allowing a 1-degree increment of the rotation angles. These results suggest that as long as the voxel size is less than half of the minimum printable feature size, the final volume of the modified lattice model is predicted with sufficient accuracy. To decrease the computational time for subsequent optimization studies, we will use this rule to set the voxel size. 4.2.2 Optimization Studies Based on the Manufacturability-Oriented Model Correction The influential parameters for the manufacturability of AM processes are material properties and geometric features. We have already performed some optimization studies on the orthotropic material strength influence on the AM lattice structure performance [reference]. Here, we are going to employ a model correction approach to analyze the manufacturability of the lattice geometries and the restriction of the AM process parameters [38]. The study aims to optimize the build orientation of the designed lattice structure so that the as-manufactured lattice structure will have the least amount of volume difference compared to the as-designed structure. The model correction approach is a voxel-based analysis method considering the manufacturability of AM processes, which slices the 3D geometry into 2D images and voxelizes 49 the 2D plane for manufacturability checking. The most crucial input parameters are the maximum number of voxels in the largest dimension and the number of voxels that constitute the smallest manufacturing size of the specific AM process. Due to the different dimensions of designs and different AM processes, these parameters should be carefully tuned for the optimization study. 4.2.2.1 Optimization Settings For the optimization study, the optimization algorithm SHERPA within the HEEDS multi- disciplinary optimization software [reference] has been used. . The FDM process was selected as the AM process for lattice manufacturing. The minimum printable feature of FDM was selected to be 0.3 mm. The number of optimization evaluation iterations was chosen to be 300. The study was conducted for various numbers of voxel sizes for further comparison and validation. (a) (b) Figure 4.5: (a) 1D beam, (b) 3D bcc lattice. Two lattice models were selected for demonstrating the influence of the orientation of the lattice beams. Model 1 is a single 1D beam with a circular cross-section of a radius of 0.5mm and a length of 10mm, Figure 4.5(a). Model 2 is a 3D lattice structure composed of two BCC unit cells, Figure 4.5(b). The dimensions of model 2 are 20 mm x 10 mm x 10 mm, and the radii are 0.5 mm. 50 For each of these structures, the optimization problem statement is the same as follows: Objective: Minimize the absolute volume difference (voxelized volume) Variables: Rotation angles of the lattice beams (𝜑𝜑, 𝜃𝜃) 𝜑𝜑 ∈ [−90°, 90°], 𝜃𝜃 ∈ [−90°, 90°] 𝜑𝜑, 𝜃𝜃 represents the rotation of the geometry about the z and x-axes, respectively. 4.2.2.2 Results As seen in Table 4.3, the optimal solutions may arrive at the different printing orientations with different sizes of voxelization. Table 4.3.a: Optimal solutions for 1D beam orientation optimization Voxel size (𝑚𝑚𝑚𝑚3 ) Optimal Orientation (𝜑𝜑, 𝜃𝜃) 0.05 (−75°, 30°) 0.025 (80°, 25°) Table 4.3.b: Optimal solutions for 3D lattice orientation optimization Voxel size (𝑚𝑚𝑚𝑚3 ) Optimal Orientation (𝜑𝜑, 𝜃𝜃) 0.05 (50°, −5°) 0.025 (90°, −70°) For model 1, when the voxel size is 0.05𝑚𝑚𝑚𝑚3 , the optimal orientation with minimum volume difference is [-75, 30], when the voxel size changes to 0.025mm^3, the optimal orientation with minimum volume difference is [80, 25]. For model 2, when the voxel size is 0.05𝑚𝑚𝑚𝑚3 , the optimal orientation with minimum volume difference is [50, -5] when the voxel size changes to 0.025𝑚𝑚𝑚𝑚3 , the optimal orientation with minimum volume difference is [90, -70]. 51 Based on the results, there are some possible hypotheses made: 1. With only the absolute difference of the model volume as the optimization objective, the optimal solution may not be reliable. 2. the studies are not converged yet, and the modified volume of the as-manufactured lattice structure may not be minimized through optimization. 3. since the 1D beam is simple, the volume difference due to different orientations is quite small, and different angles provide a suitable solution within the allowed variations 4. Because the voxel size compared to the beam size is small enough, which results in similar results between optimal orientations. Additional studies are needed to further investigate on this problem. 4.2.3 Optimization Studies Based on Both Geometry Correction and Material Orthotropy In addition to considering AM effects on material orthotropy and performance, we also take geometric manufacturability into account. Following the optimization of lattice members’ radii for a given build direction, we employ a separate model [38] to analyze the manufacturability of that design for an assumed set of manufacturing process parameters. This model determines where and how much material must be added or removed for a specific AM process in order to produce the design. The result of this is called the as-manufactured design. The mass (or volume) of the as-manufactured design will be slightly different from that of the as-designed lattice. The build direction will impact the design in terms of both material and geometrical effects. There will be a strong coupling between build direction-induced material orthotropy and the radii of lattice members needed to meet the strength criteria, as shown in a previous section. Similarly, there is often a strong coupling between build direction and the material volume of the as- manufactured part, as shown in [38]. 52 Due to the complex design space resulting from this type of coupling, a bi-level optimization strategy is used here to increase the efficiency of the design exploration. In the 1st optimization level (the upper level), the objective is to minimize the material volume of the as-manufactured lattice by modifying the printing direction. The as-manufactured design is obtained by performing a manufacturability analysis of the as-designed lattice, which is obtained in the 2nd optimization level. In the 2nd optimization level (the lower level), the objective is to minimize the material volume of the as-designed lattice for a given printing direction by modifying the member radii to account for material orthotropy in stiffness and strength. Figure 4.6: Flowchart for bi-level optimization with model correction and material orthotropy The complete problem statement for the bi-level lattice optimization study is: 1st (upper) level optimization: Objective: Minimize the material volume of the as-manufactured lattice structure As-manufactured total volume=(original voxels + added voxels - removed voxels) * 𝑣𝑣𝑣𝑣𝑣𝑣𝑣𝑣𝑣𝑣 𝑠𝑠𝑠𝑠𝑠𝑠𝑠𝑠 3 Constraints: 53 Failure indices for orthotropic material strengths 𝛼𝛼𝑖𝑖 ≤ 1 (i =1,2…5) (Values of 𝛼𝛼𝑖𝑖 are obtained from Level 2, and the constraints are defined in [17]) Variables: Printing angles θ, φ −90° ≤ 𝜑𝜑 ≤ 90° 0° ≤ 𝜃𝜃 ≤ 180° 2nd (lower) level optimization (for a fixed set of printing angles θ, φ): Objective: Minimize the material volume of the as-designed lattice structure Constraints: Failure indices for material strengths constraints 𝛼𝛼𝑖𝑖 ≤ 1 (i =1,2…5) Variables: Radii of lattice members 𝑟𝑟𝑖𝑖 (i=1 to number of lattice members) The detailed workflow of the bi-level optimization is illustrated in Figure 4.7. To achieve coherence with the former material orthotropy study, the FDM process is selected. For a good voxelization of the lattice structure with a given minimum printable feature size, the radius of the lattice beam would be equal to or greater than 3 times the voxel size. 54 Figure 4.7: Flowchart of the bi-level optimization process. Procedures within the solid line belong to the 2nd (lower) level optimization, while procedures within the dotted line belong to the 1st (upper) level optimization. The software package(s) used for each procedure are listed beside the blocks. The 1st (upper) level optimization determines the proposed printing direction first by assigning values to θ, φ. Then, for that printing direction, optimization is performed on the radii variables to minimize the material volume of the originally designed lattice within the 2nd (lower) level, while satisfying constraints on material strength. The failure criteria are defined in terms of failure indices 𝛼𝛼𝑖𝑖 , as defined in [reference of previous paper]. Thereafter, the parameters of the optimal as-designed lattice are returned to the 1st level, where a manufacturability analysis is performed to produce the as-manufactured lattice. This process is repeated until an optimal printing direction and associated as-manufactured design are obtained. 55 HEEDS, Mithril, Abaqus, Python, Manufacturability Oriented Model Correction, and MATLAB are software tools utilized in various procedures of this study. HEEDS [72] is used to perform the design space search and to automate the entire optimization process between the software packages mentioned above. Mithril is a generative design package for lattice structures [need a reference]. It allows a very large design to be represented as a few lines of code (a script). In the current optimization process, each as-designed lattice model emerging from the lower level is converted into the format of a Mithril output file (*.mtl), which is used as the input for the manufacturability- oriented model correction analysis. Abaqus [66] is a finite element analysis software package that is used to perform the structural analysis of each lattice design. Python [71] is used as the scripting language to extract simulation results from the Abaqus output (*.odb) file and to send them to MATLAB for post- processing. MATLAB [73] is employed to compute orthotropic material properties for given printing angles and post-processing the stress state at the material coordinate system, as described in [17]. 4.2.3.1 3D BCC Lattice Structure The 3D bcc lattice structure is composed of two BCC unit cells. The dimensions of the model are 20 mm * 10 mm * 10 mm, and the radii are 0.5 mm. For the loading condition, a force is applied in the negative z-direction on the lattice nodes of the right surface. On the left surface, the nodes are fixed, as shown in Figure 4.8. 56 Figure 4.8: 3D bcc lattice structure optimized results The optimal orientation of the lattice structure was found to be (𝜑𝜑, 𝜃𝜃) = (15°,0°). This orientation is shown in Figure 4.8. From the figure above, the optimal orientations of the as- designed lattice structure to the as-manufactured lattice are quite different which shows that both material orthotropy and AM process resolution are crucial for the manufacturability analysis of the lattice structure. 4.2.3.2 3D Bracket Lattice The 3D bracket lattice structure shown in Figure 4.9 was optimized using the same process. The overall dimensions of the bracket were 300 mm * 200 mm * 20 mm. The two mounting holes on the left of the bracket were fixed, and a vertical load (in the negative z-direction) was applied to the mounting hole on the right side. Figure 4.9: Bracket lattice structure optimized results 57 The optimal orientation of the bracket lattice structure was found to be (𝜑𝜑, 𝜃𝜃) = (-30°,95°) (see Figure 4.9). The original as-designed volume is 313815 𝑚𝑚𝑚𝑚3 and the as-manufactured volume of the lattice is 469,239 𝑚𝑚𝑚𝑚3 . The difference is 33%. Z Y X (a) (b) Figure 4.10: Bracket lattice structure built-in XY plane (a), XZ plane (b) Recalling that the Z direction is the printing head deposit direction, Figure 4.10 compares the optimal orientation of the as-manufactured lattice (Figure 4.10(a)) to the build orientation along Z axis (Figure 4.10(b)). Accounting for material orthotropy and AM resolution limitations, the optimal printing direction yields a lattice volume that is 10.76% lower than the direction shown in Figure 4.10(b). 58 Chapter 5 Multiphysics Design Optimization of Lattice Structures Within this study, we are going to generate lattice layouts under both mechanical and thermal loads with the field-aligned lattice algorithm. Based on the preference of different physic fields, a weight factor could be assigned within the load trajectory selection algorithm. The lattice beams would also be optimized on radii to achieve minimum mechanical and thermal compliance. This study is composed of three sections. First, we present an intermediate method and studies on this topic. Second, we introduce the method of lattice trajectory generation, which includes the algorithm of computing the principal vectors between multi-physics fields. Last, based on the multi-physics trajectory selection algorithm, we conduct a comparative study based on different weighting factors for the stress and the heat conduction field, to illustrate the trade- off of lattice layout between the different importance of multi-physics fields. 5.1 Intermediate Methods and Results 5.1.1 The Trajectory Method In order to estimate the output fields for each type of loading, we initially assume that the domain is solid and uniform. For each physics field (e.g., structural or thermal), we use finite element analysis (FEA) to predict the distribution of the responses (e.g., stress field or heat flux). The principal vectors are computed at each nodal point, and these vectors are used in the creation of the lattice trajectories, as described below. 59 5.1.2 The Multifield Trajectory The lattice structure is designed to be applied to different scenarios, stress, heat conduction, and heat convection. Within earlier studies, the lattice layout aligned with the stress field and thermal field were studied separately [5], [6]. In this study, we propose a method to generate the lattice layout considering both the stress field and heat flux field applied simultaneously. The lattice layout is based on the trajectory of the stress and heat flux. To generate the trajectory for each field, the stress and heat flux values within a finely meshed domain are employed. The stress and heat flux values are obtained from Abaqus FEA. Then based on the start point and the weight factor for each field, the trajectory of the lattice is drawn based on the vectors for each field. 5.1.3 Multiphysics Fields Vector Selection Algorithm Input: i = number of physic fields 𝛾𝛾𝑖𝑖 : user-defined global weights for each physics field 𝑢𝑢𝑖𝑖 : local computed weight based on local physics field values at the point 𝑆𝑆1 , 𝑆𝑆2 , 𝐻𝐻1 , 𝐻𝐻2 , 𝑅𝑅1 ,(𝑅𝑅2 ): principal field vectors at the point While i: 𝑢𝑢1 *𝛾𝛾1 *𝑆𝑆1 +𝑢𝑢2 *𝛾𝛾2 *𝐻𝐻1 𝑣𝑣1 = (based on the angle between two vectors) 𝑢𝑢1 *𝛾𝛾1 + 𝑢𝑢2 *𝛾𝛾2 𝑢𝑢1 *𝛾𝛾1 *𝑆𝑆2 +𝑢𝑢2 *𝛾𝛾2 *𝐻𝐻2 𝑣𝑣2 = 𝑢𝑢1 *𝛾𝛾1 + 𝑢𝑢2 *𝛾𝛾2 New weight for 𝑣𝑣1 and 𝑣𝑣2 is n= (𝑢𝑢1 *𝛾𝛾1 + 𝑢𝑢2 *𝛾𝛾2 ) 𝑛𝑛*𝑣𝑣1 +𝑢𝑢3 *𝛾𝛾3 *𝑅𝑅1 (𝑢𝑢1 *𝛾𝛾1 + 𝑢𝑢2 *𝛾𝛾2 ) *𝑣𝑣1 +𝑢𝑢3 *𝛾𝛾3 *𝑅𝑅1 𝑢𝑢1 *𝛾𝛾1 *𝑆𝑆1 +𝑢𝑢2 *𝛾𝛾2 *𝐻𝐻1 +𝑢𝑢3 *𝛾𝛾3 *𝑅𝑅1 𝑚𝑚1 = = = 𝑛𝑛+𝑢𝑢3 *𝛾𝛾3 (𝑢𝑢1 *𝛾𝛾1 + 𝑢𝑢2 *𝛾𝛾2 )+𝑢𝑢3 *𝛾𝛾3 𝑢𝑢1 *𝛾𝛾1 + 𝑢𝑢2 *𝛾𝛾2 +𝑢𝑢3 *𝛾𝛾3 𝑛𝑛*𝑣𝑣2 +𝑢𝑢3 *𝛾𝛾3 *𝑅𝑅2 (𝑢𝑢1 *𝛾𝛾1 + 𝑢𝑢2 *𝛾𝛾2 ) *𝑣𝑣2 +𝑢𝑢3 *𝛾𝛾3 *𝑅𝑅2 𝑢𝑢1 *𝛾𝛾1 *𝑆𝑆1 +𝑢𝑢2 *𝛾𝛾2 *𝐻𝐻1 +𝑢𝑢3 *𝛾𝛾3 *𝑅𝑅2 𝑚𝑚2 = = = 𝑛𝑛+𝑢𝑢3 *𝛾𝛾3 (𝑢𝑢1 *𝛾𝛾1 + 𝑢𝑢2 *𝛾𝛾2 )+𝑢𝑢3 *𝛾𝛾3 𝑢𝑢1 *𝛾𝛾1 + 𝑢𝑢2 *𝛾𝛾2 +𝑢𝑢3 *𝛾𝛾3 New weight for 𝑚𝑚1 and 𝑚𝑚2 is n= (𝑢𝑢1 *𝛾𝛾1 + 𝑢𝑢2 *𝛾𝛾2 + 𝑢𝑢3 *𝛾𝛾3 ) 60 End The weighted function terms must satisfy the condition: 𝛾𝛾1 + 𝛾𝛾2 + ⋯ + 𝛾𝛾𝑖𝑖 = 1 The local weight functions are calculated as follows: 𝑢𝑢1 = Max(Abs(Sig1), Abs(Sig2)) / Max(Abs[Globalsig1],Abs[Globalsig2]) 𝑢𝑢2 = Abs(heat flux value)/ Max (abs (global heat flux values)) 𝑢𝑢𝑖𝑖 = Abs(R field value)/ Max (abs (global R field values)) Where, 𝑢𝑢1 computed weight based on local principal stress values at the point 𝑢𝑢2 computed weight based on local heat flux value at the point 𝑢𝑢𝑖𝑖 computed weight based on local R field values at the point a. Angle Definition S2: 2nd principal stress direction 𝛽𝛽 H1: heat flux direction 𝛼𝛼 S1: 1st principal stress direction Figure 5.1: the angle definition of the principal vectors S1= 1st principal stress direction S2=2nd principal stress direction H1= heat flux direction H2=direction orthogonal to heat flux direction H1 𝛼𝛼 angle: from 1st principal stress direction S1 to heat flux direction H1 61 𝛽𝛽 angle: between 2nd principal stress direction S2 and heat flux direction H1 b. Vector Selection 1. 𝛼𝛼=[0, 𝜋𝜋/2 , 𝜋𝜋 , -𝜋𝜋/2] H1 S1 S1 H1 (a) (b) H S 𝛼𝛼 𝛼𝛼 S H (c) (d) Figure 5.2: The angle 𝜶𝜶 between S1 and H1 equals 𝝅𝝅 (a), 0 (b), 𝝅𝝅/2 (c), −𝝅𝝅/2 (d) In this case, select 2 principal stress directions S1, S2, or heat flux direction H1 and its orthogonal direction H2 as the lattice trajectory vectors, depending on the weighting factor of each field: V1= S1 or H1 V2= S2 or H2 2. -𝜋𝜋/4 ≤ 𝛼𝛼 ≤ 𝜋𝜋/4 [𝛼𝛼 ≠ 0] 62 (a) (b) Figure 5.3: The angle between S1 and H1, 𝟎𝟎 < 𝜶𝜶 ≤ 𝝅𝝅/4 (a), -𝝅𝝅/4 ≤ 𝜶𝜶 < 0 (b). When −𝜋𝜋/4 ≤ 𝛼𝛼 ≤ 𝜋𝜋/4 [𝛼𝛼 ≠ 0], select the 2 principal vectors based on the following equations: w1* γ1*S1+ w2*γ2*H1 w1* γ1*S2+ w2*γ2*H2 V1= V2= w1* γ1+ w2*γ2 w1* γ1+ w2*γ2 3. 𝜋𝜋/2 < |𝛼𝛼|< 𝜋𝜋/4 (a) (b) Figure 5.4: The angle between S1 and H1. (a) 𝝅𝝅/2 < |𝜶𝜶|< 𝝅𝝅/4, 𝝅𝝅/𝟒𝟒 < 𝜶𝜶 ≤ 𝝅𝝅/𝟐𝟐 & 𝜷𝜷 < 𝝅𝝅/2, (b) -𝝅𝝅/2 < 𝜶𝜶 < -𝝅𝝅/4 & 𝜷𝜷 > 𝝅𝝅/2. When 𝜋𝜋/4 < 𝛼𝛼 ≤ 𝜋𝜋/2 & 𝛽𝛽 < 𝜋𝜋/2, Cos(𝛽𝛽)>0, then w1* γ1*S2+ w2*γ2*H1 w1* γ1*(-S1)+ w2*γ2*H2 V1= V2= w1* γ1+ w2*γ2 w1* γ1+ w2*γ2 When -𝜋𝜋/2 < 𝛼𝛼 < -𝜋𝜋/4 & 𝛽𝛽 < 𝜋𝜋/2, Cos(𝛽𝛽)<0, then w1* γ1*(-S2)+ w2*γ2*H1 w1* γ1*S1+ w2*γ2*H2 V1= V2= w1* γ1+ w2*γ2 w1* γ1+ w2*γ2 4. 𝜋𝜋/2 < |𝛼𝛼 |<3𝜋𝜋/4 63 (a) (b) Figure 5.5: The angle between S1 and H1 𝝅𝝅/2 < |𝜶𝜶 |<𝟑𝟑𝝅𝝅/4, 𝜋𝜋/2 < 𝛼𝛼 ≤ 3𝜋𝜋/4 & 𝛽𝛽 < 𝜋𝜋/2 (a), -3𝜋𝜋/4 < 𝛼𝛼 < 𝜋𝜋/2 & 𝛽𝛽 > 𝜋𝜋/2 (b). When 𝜋𝜋/2 < 𝛼𝛼 ≤ 3𝜋𝜋/4 & 𝛽𝛽 < 𝜋𝜋/2, shown in Figure 5.5(a), Cos(𝛽𝛽)>0, the vector is w1* γ1*|S2|+ w2*γ2*H1 w1* γ1*(-S1)+ w2*γ2*H2 V1= V2= w1* γ1+ w2*γ2 w1* γ1+ w2*γ2 When -3𝜋𝜋/4 < 𝛼𝛼 < 𝜋𝜋/2 & 𝛽𝛽 > 𝜋𝜋/2, shown in Figure 5.5(b), Cos(𝛽𝛽)<0, the vector is w1* γ1*(-S2)+ w2*γ2*H1 w1* γ1*S1+ w2*γ2*H1 V1= V2= w1* γ1+ w2*γ2 w1* γ1+ w2*γ2 5. 3𝜋𝜋/4 <|𝛼𝛼 |< 𝜋𝜋 (a) (b) Figure 5.6: The angle between S1 and H1 is 3𝝅𝝅/4 < 𝜶𝜶 ≤ 𝝅𝝅 (a), -3𝝅𝝅/4 < 𝜶𝜶 < 𝝅𝝅 (b). In this case, select 2 lattice trajectory vectors: w1* γ1*(-S1)+ w2*γ2*H1 w1* γ1*(-S2)+ w2*γ2*H2 V1= V2= w1* γ1+ w2*γ2 w1* γ1+ w2*γ2 64 5.1.4 Studies and Intermediate Results 5.1.4.1 Model Settings The model used for the lattice trajectory generation is the plate with a hole. The dimensions of the plate are 40 mm* 20 mm*2 mm, in length, height, and thickness respectively. The model is shown in Figure 5.7. The origin is at the center of the circle and the radius of the hole is 5 mm. Figure 5.7: Plate with a hole model For the stress field, we use the same loading and boundary conditions as the study [5], shown in Figure 5.8. The left edge of the plate is fixed, and the force is applied to the right edge of the plate. For the material properties of the given model, the elasticity is 1 and the poison ratio is 0.33. Figure 5.8: Load and boundary conditions for stress field For the thermal field, we adopt a similar surface-to-point problem as [6]. Within the design domain, heat flux q = 1/12 is applied at the left edge, top half, and bottom half of the plate. At 65 the right edge, a constant temperature To = 0 is set at the 2 mm region in the middle. The load and boundary conditions for the heat conduction field are shown in Figure 5.9. The heat conductivity for the model is set to 1. 2mm Heat Flux 𝑞𝑞𝑒𝑒 Temperature 𝑇𝑇0 =0 Figure 5.9: Load and boundary conditions for heat conduction field 5.1.4.2 Optimization Problem Statement The optimization for the lattice structure would be set into bi-level optimization, at the 1st level, the lattice trajectory would be optimized on the input parameters of the field-aligned lattice algorithm, which includes the start point, the gap between seed points, the radius of merge circle and spacing constant. At the 2nd level, for a specific lattice layout, the radii would be optimized by MMA to minimize the weighted sum compliance of both mechanical and thermal compliance. 1st level optimization: 𝐶𝐶𝑚𝑚𝑚𝑚𝑚𝑚ℎ𝑎𝑎𝑎𝑎𝑎𝑎𝑎𝑎𝑎𝑎𝑎𝑎 𝐶𝐶𝑡𝑡ℎ𝑒𝑒𝑒𝑒𝑒𝑒𝑒𝑒𝑒𝑒 Objective: minimize compliance C = 𝑤𝑤1 + 𝑤𝑤2 𝐶𝐶𝑚𝑚_𝑜𝑜𝑜𝑜𝑜𝑜𝑜𝑜𝑜𝑜𝑜𝑜𝑜𝑜𝑜𝑜 𝐶𝐶𝑡𝑡_𝑜𝑜𝑜𝑜𝑜𝑜𝑜𝑜𝑜𝑜𝑜𝑜𝑜𝑜𝑜𝑜 Constraints: combined compliance C of the lattice structure is larger than 0(>1e-5) Variable: 𝑥𝑥 and 𝑦𝑦 coordinate of the start point (𝑥𝑥 2 + 𝑦𝑦 2 > 52 ) 𝑀𝑀𝑀𝑀𝑀𝑀𝑀𝑀𝑀𝑀𝑀𝑀, 𝑀𝑀𝑀𝑀𝑀𝑀𝑀𝑀𝑀𝑀𝑀𝑀 (maximum and minimum allowable distance between the seed points) 66 the radius of the merged circle: 𝑅𝑅 spacing constant: 𝛼𝛼 2nd level optimization: (with the optimized lattice layout) Objective: 𝐶𝐶𝑚𝑚𝑚𝑚𝑚𝑚ℎ𝑎𝑎𝑎𝑎𝑎𝑎𝑎𝑎𝑎𝑎𝑎𝑎 𝐶𝐶𝑡𝑡ℎ𝑒𝑒𝑒𝑒𝑒𝑒𝑒𝑒𝑒𝑒 Minimize compliance C= 𝑤𝑤1 + 𝑤𝑤2 𝐶𝐶𝑚𝑚_𝑜𝑜𝑜𝑜𝑜𝑜𝑜𝑜𝑜𝑜𝑜𝑜𝑜𝑜𝑜𝑜 𝐶𝐶𝑡𝑡_𝑜𝑜𝑜𝑜𝑜𝑜𝑜𝑜𝑜𝑜𝑜𝑜𝑜𝑜𝑜𝑜 Constraints: The total volume of the lattice structure is less than 20% 𝑉𝑉𝑡𝑡𝑡𝑡𝑡𝑡𝑡𝑡𝑡𝑡 Variable: Radii of lattice members: 𝑟𝑟𝑖𝑖 : 𝑟𝑟𝑖𝑖𝑚𝑚𝑚𝑚𝑚𝑚 ≤ 𝑟𝑟𝑖𝑖 ≤ 𝑟𝑟𝑖𝑖𝑚𝑚𝑚𝑚𝑚𝑚 (𝑖𝑖 = 1, 2, … 𝑡𝑡𝑜𝑜𝑜𝑜𝑜𝑜𝑜𝑜 𝑛𝑛𝑛𝑛𝑛𝑛𝑛𝑛𝑛𝑛𝑛𝑛 𝑜𝑜𝑜𝑜 𝑙𝑙𝑙𝑙𝑙𝑙𝑙𝑙𝑙𝑙𝑙𝑙𝑙𝑙 𝑚𝑚𝑚𝑚𝑚𝑚𝑚𝑚𝑚𝑚𝑚𝑚𝑚𝑚) where 𝑟𝑟𝑖𝑖𝑚𝑚𝑚𝑚𝑚𝑚 and 𝑟𝑟𝑖𝑖𝑚𝑚𝑚𝑚𝑚𝑚 represent the minimum and maximum allowable radius for the 𝑖𝑖𝑡𝑡ℎ lattice member 5.1.4.3 Weighting Factor Comparison Studies a. Mechanical field weighting factor =1, thermal field weighting factor =0. To begin the comparison study, I start with a pure stress field, in which the weighting factor for the stress field is set to 1. The lattice layout of the optimization is shown in Figure 5.10. Figure 5.10: Lattice layout for mechanical weighting factor =1 67 The temperature distribution of the optimal solution is shown in Figure 5.11, the unit of the temperature is K. As a result, the thermal compliance is 224.2889. Figure 5.12 is the stress plot of the lattice layout, and the stress compliance is 0.1804. Figure 5.11: Temperature distribution of the lattice layout Figure 5.12: The stress distribution of the lattice layout which is calculated at the midpoint of each lattice beam, color shows the stress selected based on the maximum absolute value of the compression and tension stress values. b. Mechanical field weighting factor =0, thermal field weighting factor =1. To compare with the former study results, another study with a thermal field weighting factor equal to 1 is also conducted. The optimal lattice layout is shown below. 68 Figure 5.13: Lattice layout for thermal weighting factor =1 Figure 5.14: Temperature distribution of the lattice layout Figure 5.15: The stress distribution of the lattice layout which is calculated at the midpoint of each lattice beam, color shows the stress selected based on the maximum absolute value of the compression and tension stress values. 69 The temperature distribution of the optimal solution is shown in Figure 5.14, the unit of the temperature is K. As a result, the thermal compliance is 674. Figure 5.15 is the stress plot of the lattice layout, and the stress compliance is 2.0526. 5.2 Methodology To generate a lattice layout for the solid domain, typically, the topology optimization method is used to optimize the material layout for the solid domain, then distribute lattice unit cells to the optimized layout or mesh element. Within this study, a novel algorithm is developed to generate the lattice trajectory directly considering multi-physics fields. The algorithm was invented to address the following problems: 1. Combine multiple physics fields, considering the characteristics/properties of each field. For example, within the mechanical field, the stress components are represented by 2D tensors; but in the thermal field, the heat flux is 1D vectors. How to combine 2D tensor and 1D vector is the problem that should be addressed by the algorithm. 2. Account for the importance of each field at local points. Since multi-physics fields could introduce different components within the domain, they will have different values and even differences in magnitude. 3. Implement the designer’s preference and input within the layout generation. For a lattice design, the designer may have some requirements or preferences for each field. This algorithm can take consider of designer’s bias. 4. Be generic and could be applied to multiple physics fields. With this criterion in mind, we want to develop the algorithm as a generic method that could be applied to different physics fields rather than be tuned for each problem. 70 5.2.1 The Trajectory Method From studies [54]–[57], the concept and theory of the lattice design based on the trajectory of the principal field were introduced. Gao et al.,[5] proposed the method to generate lattice with a principal stress trajectory. For the mechanical field, the way to generate a stress trajectory is: first, conduct a Finite Element Analysis (FEA) for a uniform solid domain, with the given boundary conditions and mechanical loads. Since the interpolation of stress values will be used in the following steps, a really fine mesh is strongly recommended for the analysis. Then, with the stress components generated from the FEA study at each point, the principal stress vectors could be computed at each point. The 1st principal stress trajectory would go along the 1st principal stress vectors and the trajectory orthogonal to the 1st principal stress trajectory would go along the 2nd principal stress vectors. For the thermal field, the same trajectory idea could be applied to heat flux vectors. With a thermal analysis of the solid domain, the heat flux vectors could be computed. With the same algorithm from [5], the heat flux trajectory could be generated. After the trajectories are generated for the whole domain, we would distribute the lattice elements along the trajectories to convert the trajectories to lattice layout for further studies. 71 Figure 5.16: Flowchart for lattice trajectory method. 5.2.2 The Multiphysics Field Trajectory The lattice structure is designed to be applied to different physics fields at the same time, including mechanical, heat conduction, heat convection, vibration, etc. Within earlier studies, the lattice layout aligned with the principal stress trajectories and heat flux trajectories is studied separately[5], [6]. In this study, we propose a method to generate the lattice layout considering both the principal stress field and heat flux field. The lattice layout is based on the trajectory of the principal stress and heat flux. To generate the trajectory for each field, the stress and heat flux values within a finely meshed domain are employed. The stress and heat flux values are obtained from Abaqus FEA. Then based on the start point and the weight factor for each field, the trajectory of the lattice is drawn based on the vectors for each field. 72 5.2.3 Multiphysics Fields Pseudo Principal Vector Computation Algorithm Since the trajectory generation is based on the principal vectors at each point, the first problem we want to solve is how to generate the pseudo principal vectors with two physics fields. As we know, mechanical fields and thermal fields have different values and properties. Within the same 2D solid domain, the stress components from the mechanical field are represented by 2D tensors. However, for the thermal field, the heat flux is a 1D vector. The algorithm to compute the pseudo principal vector for 2 physics fields is discussed below. 5.2.3.1 Vector Transformation In the global XY coordinate systems, the mechanical/stress field has 3 stress components, Sx, Sy, and Sxy, the heat flux field has heat flux H, shown in Figure 5.17. Figure 5.17: Stress tensor and heat flux vector in the global coordinate system Stress components in XY coordinate could be represented in the 2D tensor: 𝑆𝑆𝑥𝑥 𝑆𝑆𝑥𝑥𝑥𝑥 𝑆𝑆 = � � (1) 𝑆𝑆𝑥𝑥𝑥𝑥 𝑆𝑆𝑦𝑦 73 In order to combine the 1D vector with a 2D tensor, we want to transform the heat flux vector into a 2D tensor within XY coordinates. 𝐻𝐻1 is the 1D heat flux vector, which could be regarded as a 2D tensor with only one non-zero eigenvalue in principal coordinates: 𝐻𝐻1 0 𝐻𝐻1 = � � (2) 0 0 H in the XY coordinate would then be: 𝐻𝐻𝑥𝑥 𝐻𝐻𝑥𝑥𝑥𝑥 𝐻𝐻 = � � (3) 𝐻𝐻𝑥𝑥𝑥𝑥 𝐻𝐻𝑦𝑦 From the 2nd order tensor transformation equation, the principal tensor 𝐻𝐻1 𝐻𝐻1 = 𝑄𝑄 ∗ 𝐻𝐻 ∗ 𝑄𝑄𝑇𝑇 (4) 𝐻𝐻1 0 cos 𝜃𝜃 −sin 𝜃𝜃 𝐻𝐻𝑥𝑥 𝐻𝐻𝑥𝑥𝑥𝑥 cos 𝜃𝜃 sin 𝜃𝜃 � �=� �� �� � (5) 0 0 sin 𝜃𝜃 cos 𝜃𝜃 𝐻𝐻𝑥𝑥𝑥𝑥 𝐻𝐻𝑦𝑦 −sin 𝜃𝜃 cos 𝜃𝜃 The 𝜃𝜃 is the angle between the heat flux vector 𝐻𝐻1 and X axis. Then, the heat tensor H in the XY coordinate is 𝐻𝐻 = 𝑄𝑄𝑇𝑇 ∗ 𝐻𝐻1 ∗ 𝑄𝑄 (6) 𝐻𝐻𝑥𝑥 𝐻𝐻𝑥𝑥𝑥𝑥 cos 𝜃𝜃 sin 𝜃𝜃 𝐻𝐻1 0 cos 𝜃𝜃 − sin 𝜃𝜃 � �=� �� �� � 𝐻𝐻𝑥𝑥𝑥𝑥 𝐻𝐻𝑦𝑦 −sin 𝜃𝜃 cos 𝜃𝜃 0 0 sin 𝜃𝜃 cos 𝜃𝜃 𝐻𝐻1 cos 𝜃𝜃 cos 𝜃𝜃 −𝐻𝐻1 cos 𝜃𝜃 sin 𝜃𝜃 =� � (7) −𝐻𝐻1 cos 𝜃𝜃 sin 𝜃𝜃 𝐻𝐻1 sin 𝜃𝜃 sin 𝜃𝜃 Now, with both stress components and heat flux components in the 2D tensor format, we could superimpose them together. 5.2.3.2 Normalization Considering the values of components from each physics field, normalization is needed for the combination of components. 74 𝑆𝑆𝑣𝑣𝑣𝑣𝑣𝑣 = 𝐴𝐴𝐴𝐴𝐴𝐴(�𝜎𝜎𝑥𝑥 2 + 𝜎𝜎𝑦𝑦 2 − 𝜎𝜎𝑥𝑥 𝜎𝜎𝑦𝑦 + 3𝜏𝜏𝑥𝑥𝑥𝑥 2 ) (8) 𝐻𝐻𝑎𝑎𝑎𝑎𝑎𝑎 = 𝐴𝐴𝐴𝐴𝐴𝐴(𝐴𝐴𝐴𝐴𝐴𝐴(𝐻𝐻1 )) (9) In order to consider the sign of each component, the normalization factors were chosen to be positive values. The normalization factor for the stress components is the average of all von- mises stress within the domain. For the heat flux, since it is a 1D vector, the normalization factor is the average of absolute values of heat flux. For each component, the normalized components are: 𝑆𝑆𝑥𝑥 𝑆𝑆𝑦𝑦 𝑆𝑆𝑥𝑥𝑥𝑥 𝑢𝑢𝑥𝑥 = 𝑢𝑢𝑦𝑦 = 𝑢𝑢𝑥𝑥𝑥𝑥 = (10) 𝑆𝑆𝑣𝑣𝑣𝑣𝑣𝑣 𝑆𝑆𝑣𝑣𝑣𝑣𝑣𝑣 𝑆𝑆𝑣𝑣𝑜𝑜𝑜𝑜 𝐻𝐻𝑥𝑥 𝐻𝐻𝑦𝑦 𝐻𝐻𝑥𝑥𝑥𝑥 𝑣𝑣𝑥𝑥 = 𝑣𝑣𝑦𝑦 = 𝑣𝑣𝑥𝑥𝑥𝑥 = 𝐻𝐻𝑎𝑎𝑎𝑎𝑎𝑎 𝐻𝐻𝑎𝑎𝑎𝑎𝑎𝑎 𝐻𝐻𝑎𝑎𝑎𝑎𝑎𝑎 5.2.3.3 Designer Preference With a multi-physics problem, the designer may have a design preference for each field. The weighting factor defined with the designer’s preference is as follows: 𝛾𝛾1 + 𝛾𝛾2 + ⋯ + 𝛾𝛾𝑖𝑖 = 1 (11) 𝛾𝛾1 : Given weight for mechanical field 𝛾𝛾2 : Given weight for thermal field 𝛾𝛾𝑖𝑖 : Given weight for any physics field i 5.2.3.4 Weighted Sum Function We use the weighted sum function in the global coordinate: 𝑤𝑤𝑥𝑥 = 𝑢𝑢𝑥𝑥 *𝛾𝛾1 +𝑣𝑣𝑥𝑥 *𝛾𝛾2 (12) 𝑤𝑤𝑦𝑦 = 𝑢𝑢𝑦𝑦 *𝛾𝛾1 +𝑣𝑣𝑦𝑦 *𝛾𝛾2 (13) 𝑤𝑤𝑥𝑥𝑥𝑥 = 𝑢𝑢𝑥𝑥𝑥𝑥 *𝛾𝛾1 + 𝑣𝑣𝑥𝑥𝑥𝑥 *𝛾𝛾2 (14) 75 5.2.3.5 Tensor Combination The combined tensor is: 𝑤𝑤𝑥𝑥 𝑤𝑤𝑥𝑥𝑥𝑥 𝑊𝑊 = [𝑤𝑤 𝑤𝑤𝑦𝑦 ] (15) 𝑥𝑥𝑥𝑥 5.2.3.6 Pseudo Principal Vector in 2D Calculate the pseudo principal vector and principal values by finding the eigenvector and eigenvalues from the combined tensor: 𝑤𝑤𝑥𝑥 +𝑤𝑤𝑦𝑦 𝑤𝑤𝑥𝑥 −𝑤𝑤𝑦𝑦 2 𝑤𝑤1,2 = ± �� � + 𝑤𝑤𝑥𝑥𝑥𝑥 2 (16) 2 2 With the pseudo principal vectors computed, the new trajectories and lattice layout could be generated. 5.2.4 Multiphysics Fields Lattice Optimization Within this study, we are going to use the optimization method to test and demonstrate the ability of the algorithm and the influence of the designer’s preference on the design of lattice layout. We use multi-objective optimization studies to investigate the design parameters and weighting factors’ impact on the lattice layout. One objective is to minimize thermal compliance, and another is to minimize the mechanical compliance of the structure. Since for each lattice layout, the radii of the lattice beams would also impact the objectives, the 2nd level of optimization was adapted to minimize the weighted sum of both objectives by MMA [60]. The weighted sum approach works for convex pareto fronts, but may not work well for non-convex fronts. 76 1st level optimization: Objective: minimize mechanical compliance 𝐶𝐶𝑚𝑚𝑚𝑚𝑚𝑚ℎ𝑎𝑎𝑎𝑎𝑎𝑎𝑎𝑎𝑎𝑎𝑎𝑎 minimize mechanical compliance 𝐶𝐶𝑡𝑡ℎ𝑒𝑒𝑒𝑒𝑒𝑒𝑒𝑒𝑒𝑒 Constraints: Each compliance C of the lattice structure is larger than 0(>1e-5) Variable: 𝑥𝑥 and 𝑦𝑦 coordinate of the start point 𝑀𝑀𝑀𝑀𝑀𝑀𝑀𝑀𝑀𝑀𝑀𝑀, 𝑀𝑀𝑀𝑀𝑀𝑀𝑀𝑀𝑀𝑀𝑀𝑀 (maximum, minimum allowable distance between the seed points) the radius of the merged circle: 𝑅𝑅 spacing constant: 𝛼𝛼 the weighting factor for each field: 𝛾𝛾1 , 𝛾𝛾2 2nd level optimization: (with the optimized lattice layout) Objective: 𝐶𝐶𝑚𝑚𝑚𝑚𝑚𝑚ℎ𝑎𝑎𝑎𝑎𝑎𝑎𝑎𝑎𝑎𝑎𝑎𝑎 𝐶𝐶𝑡𝑡ℎ𝑒𝑒𝑒𝑒𝑒𝑒𝑒𝑒𝑒𝑒 Minimize compliance C= 𝛾𝛾1 + 𝛾𝛾2 (𝛾𝛾2 = 1 − 𝛾𝛾1 ) 𝐶𝐶𝑚𝑚_𝑜𝑜𝑜𝑜𝑜𝑜𝑜𝑜𝑜𝑜𝑜𝑜𝑜𝑜 𝐶𝐶𝑡𝑡_𝑜𝑜𝑜𝑜𝑜𝑜𝑜𝑜𝑜𝑜𝑜𝑜𝑜𝑜 𝐶𝐶𝑚𝑚_𝑜𝑜𝑜𝑜𝑜𝑜𝑜𝑜𝑜𝑜𝑜𝑜𝑜𝑜 , 𝐶𝐶𝑡𝑡_𝑜𝑜𝑜𝑜𝑜𝑜𝑜𝑜𝑜𝑜𝑜𝑜𝑜𝑜 are the optimal solution from the single optimization study Constraints: The total volume of the lattice structure is less than 10% 𝑉𝑉𝑡𝑡𝑡𝑡𝑡𝑡𝑡𝑡𝑡𝑡 Variable: Radii of lattice members: 𝑟𝑟𝑖𝑖𝑚𝑚𝑚𝑚𝑚𝑚 ≤ 𝑟𝑟𝑖𝑖 ≤ 𝑟𝑟𝑖𝑖𝑚𝑚𝑚𝑚𝑚𝑚 𝑖𝑖 = 1, 2, … 𝑡𝑡𝑡𝑡𝑡𝑡𝑡𝑡𝑡𝑡 𝑛𝑛𝑛𝑛𝑛𝑛𝑛𝑛𝑛𝑛𝑛𝑛 𝑜𝑜𝑜𝑜 𝑙𝑙𝑙𝑙𝑙𝑙𝑙𝑙𝑙𝑙𝑙𝑙𝑙𝑙 𝑚𝑚𝑚𝑚𝑚𝑚𝑚𝑚𝑚𝑚𝑚𝑚𝑚𝑚 where 𝑟𝑟𝑖𝑖𝑚𝑚𝑚𝑚𝑚𝑚 and 𝑟𝑟𝑖𝑖𝑚𝑚𝑚𝑚𝑚𝑚 represent the minimum and maximum allowable radius for the 𝑖𝑖𝑡𝑡ℎ lattice member. 77 5.3 Studies and Results Based on the method and the algorithm discussed in the sections above, we deployed two examples to demonstrate the approach. One is a cantilever beam model, and another is the plate with a hole model. 5.3.1 Cantilever Beam Example First, an example of a cantilever beam is deployed for demonstrating the approach. 5.3.1.1 Model Settings The dimensions of the plate are 40 mm* 20 mm*1 mm, in length, height, and thickness respectively. The representation of the cantilever beam is shown in Figure 5.18. Figure 5.18: Load and boundary conditions for mechanical field and heat conduction field To show the contrast of lattice layout due to two physics fields, we set the loading and boundary conditions in Figure 5.18 so as to maximize the trade-off between designs that are ideal for mechanical versus thermal loading conditions. 78 For the mechanical field, we use the loading and boundary conditions as the cantilever beam, as shown in Figure 5.18. The left edge of the plate is fixed, and the uniform distributed force F=5 N/mm is applied on the top edge of the beam. For the material properties of the given model, the modulus of elasticity is 2078 MPa and the poison ratio is 0.33. For the thermal field, the loading and boundary conditions are assigned to the design domain as shown in Figure 5.18. The heat flux 𝑞𝑞𝑒𝑒 = 1 W/𝑚𝑚𝑚𝑚2 is applied at the center of the left edge. At the right edge, a constant temperature 𝑇𝑇0 = 0 ℃ = 273 K is set at two 4 mm regions at two corners of the right edge. The heat conductivity for the model is set to 0.1 𝑊𝑊/𝑚𝑚 ∙ 𝐾𝐾. 5.3.1.2 Optimization Problem Statement A bi-level optimization procedure was used. In the 1st level, the lattice trajectory was optimized on the orientation input parameters of the field-aligned lattice algorithm, which includes the start point, the gap between seed points, the radius of the merging circle, and the spacing constant. In the 2nd level, for a specific lattice layout, the radii were optimized by MMA to minimize the weighted sum compliance of both mechanical and thermal compliance. 1st level optimization: Objective: minimize mechanical compliance 𝐶𝐶𝑚𝑚𝑚𝑚𝑚𝑚ℎ𝑎𝑎𝑎𝑎𝑎𝑎𝑎𝑎𝑎𝑎𝑎𝑎 minimize mechanical compliance 𝐶𝐶𝑡𝑡ℎ𝑒𝑒𝑒𝑒𝑒𝑒𝑒𝑒𝑒𝑒 Constraints: Each compliance C of the lattice structure is larger than 0(>1e-5) Variable: 𝑥𝑥 and 𝑦𝑦 coordinate of the start point 𝑀𝑀𝑀𝑀𝑀𝑀𝑀𝑀𝑀𝑀𝑀𝑀, 𝑀𝑀𝑀𝑀𝑀𝑀𝑀𝑀𝑀𝑀𝑀𝑀 (maximum, minimum allowable distance between the seed points) 79 the radius of the merged circle: 𝑅𝑅 spacing constant: 𝛼𝛼 the weighting factor for the mechanical field: 𝛾𝛾1 2nd level optimization: (with the optimized lattice layout) Objective: 𝐶𝐶𝑚𝑚𝑚𝑚𝑚𝑚ℎ𝑎𝑎𝑎𝑎𝑎𝑎𝑎𝑎𝑎𝑎𝑎𝑎 𝐶𝐶𝑡𝑡ℎ𝑒𝑒𝑒𝑒𝑒𝑒𝑒𝑒𝑙𝑙 Minimize compliance C= 𝛾𝛾1 + 𝛾𝛾2 (𝛾𝛾2 = 1 − 𝛾𝛾1 ) 𝐶𝐶𝑚𝑚_𝑜𝑜𝑜𝑜𝑜𝑜𝑜𝑜𝑜𝑜𝑜𝑜𝑜𝑜 𝐶𝐶𝑡𝑡_𝑜𝑜𝑜𝑜𝑜𝑜𝑜𝑜𝑜𝑜𝑜𝑜𝑜𝑜 𝐶𝐶𝑚𝑚_𝑜𝑜𝑜𝑜𝑜𝑜𝑜𝑜𝑜𝑜𝑜𝑜𝑜𝑜 , 𝐶𝐶𝑡𝑡_𝑜𝑜𝑜𝑜𝑜𝑜𝑜𝑜𝑜𝑜𝑜𝑜𝑜𝑜 are the optimal solution from the single optimization study Constraints: The total volume of the lattice structure is less than 10% 𝑉𝑉𝑡𝑡𝑡𝑡𝑡𝑡𝑡𝑡𝑡𝑡 Variable: Radii of lattice members: 𝑟𝑟𝑖𝑖𝑚𝑚𝑚𝑚𝑚𝑚 ≤ 𝑟𝑟𝑖𝑖 ≤ 𝑟𝑟𝑖𝑖𝑚𝑚𝑚𝑚𝑚𝑚 𝑖𝑖 = 1, 2, … 𝑡𝑡𝑡𝑡𝑡𝑡𝑡𝑡𝑡𝑡 𝑛𝑛𝑛𝑛𝑛𝑛𝑛𝑛𝑛𝑛𝑛𝑛 𝑜𝑜𝑜𝑜 𝑙𝑙𝑙𝑙𝑙𝑙𝑙𝑙𝑙𝑙𝑙𝑙𝑙𝑙 𝑚𝑚𝑚𝑚𝑚𝑚𝑚𝑚𝑚𝑚𝑚𝑚𝑚𝑚 where 𝑟𝑟𝑖𝑖𝑚𝑚𝑚𝑚𝑚𝑚 and 𝑟𝑟𝑖𝑖𝑚𝑚𝑚𝑚𝑚𝑚 represent the minimum and maximum allowable radius for the 𝑖𝑖𝑡𝑡ℎ lattice member 5.3.1.3 Weighting Factor Comparison Studies After testing the impact of different weighting factors on the layout of the lattice structure, the results are summarized in the following table of different combinations of weighting factors. From the lattice layouts shown in Table 5.1, it is clear that as the weighting factor for mechanical loading changes from 1 to 0, the lattice layout varies from preferring mechanical loading to biasing the thermal loading. The mechanical compliance increases and the thermal compliance decreases. There is a strong tradeoff of the lattice designs between mechanical compliance and thermal compliance. 80 Table 5.1: Lattice layout under different weighting factors Mechanical The weighting compliance (mm/N) / factor for Optimized Lattice layout Thermal compliance mechanical field (103 𝐾𝐾/𝑊𝑊) 1.1227/ 1 37414 1.9194/ 0.7 9776.8 2.3036/ 0.5 8072.3 8.58/ 0.3 8707.3 114500/ 0 4833.8 81 In order to investigate the tradeoff of the lattice design, we used multi-objective optimization to systematically investigate the impact of the weighting factor and design parameters on the design of the lattice layout. We used the bi-level optimization strategy discussed in subsection 5.3.1.2. Within the 2nd level optimization, we used MMA to minimize the weighted sum of two objectives. To help with comparing the values for each objective within the weighted sum function, normalization factors are necessary for normalizing each objective. We deployed two single objective optimization studies for simply mechanical loading conditions and simply thermal loading conditions to find the normalization factors for each objective. The optimal values for each study are shown in the table below. Table 5.2: The optimal solution for single objective optimization Single optimization study Mechanical compliance Thermal compliance Mechanical 1.1227 37414 Thermal 114500 4833.8 With the normalization factors for mechanical compliance and thermal compliance, we conducted the bi-level optimization with three sets of baseline designs to test the consistency of the results. The multi-objective optimization study results are presented in the figure below. 82 Figure 5.19: Pareto fronts from 3 sets of bilevel optimization studies With three sets of bilevel optimization studies, the Pareto fronts show the same tendency in the distribution of designs. And the weight factor of mechanical is decreasing as the designs move from the upper left corner to the lower right corner. In Figure 5.19, two optimal solutions from two single objective optimization studies are also shown. The plot demonstrates that two extreme cases are distributed in the two corners. Most designs in the knee region could be more interesting to explore the tradeoff between each objective. Besides three sets of optimization studies, we also did a convergence test ahead to choose a reasonable number for all optimization studies. The hypervolume for different evaluation numbers is calculated in the table and figure below. It’s clear that the study gets converged at 800 evaluations. 83 Table 5.3: Hypervolume for optimization evaluation numbers Evaluation numbers Hypervolume 50 0.2953 100 0.4871 200 0.5401 400 0.4820 600 0.4872 800 0.5564 1000 0.4872 Hypervolume 0.6 0.5564 0.5401 0.4871 0.482 0.4872 0.4872 0.5 0.4 Hypervolume 0.2953 0.3 0.2 0.1 0 50 100 200 400 600 800 1000 Evaluations Number Figure 5.20: Hypervolume for different evaluation numbers. 84 5.3.2 Plate with a Hole Example Second, an example of a plate with a hole using the approach is presented below. (The study hasn’t been completed yet, some preliminary results are shown in this section) 5.3.2.1 Model Settings The model used for the lattice trajectory generation is the plate with a hole. The dimensions of the plate are 40 mm* 20 mm*1 mm, in length, height, and thickness respectively. The model is shown in Figure 5.21. The origin is at the center of the circle and the radius of the hole is 5 mm. Figure 5.21: Plate with a hole model. For the mechanical field, we use the loading and boundary conditions shown in Figure 5.22. The left edge of the plate is fixed, and the uniform distributed force F=10 N/mm is applied to the right edge of the plate. For the material properties of the given model, the modulus of elasticity is 2078 MPa and the poison ratio is 0.33. For the thermal field, the loading and boundary conditions are assigned to the design domain as shown in Figure 5.21. The heat flux 𝑞𝑞𝑒𝑒 = 1 W/𝑚𝑚𝑚𝑚2 is applied to the center of the bottom edge. A constant temperature 𝑇𝑇0 = 0 ℃ = 273 K is set to left, right and top edge. The circle and the bottom edge are adiabatic. The heat conductivity for the model is set to 0.1 𝑊𝑊/𝑚𝑚 ∙ 𝐾𝐾. 85 Figure 5.21: Load and boundary conditions for mechanical field and heat conduction field 5.3.2.2 Optimization Problem Statement A bi-level optimization procedure was used. In the 1st level, the lattice trajectory was optimized on the orientation input parameters of the field-aligned lattice algorithm, which includes the start point, the gap between seed points, the radius of the merging circle, and the spacing constant. In the 2nd level, for a specific lattice layout, the radii were optimized by MMA to minimize the weighted sum compliance of both mechanical and thermal compliance. 1st level optimization: Objective: minimize mechanical compliance 𝐶𝐶𝑚𝑚𝑚𝑚𝑚𝑚ℎ𝑎𝑎𝑎𝑎𝑎𝑎𝑎𝑎𝑎𝑎𝑎𝑎 minimize mechanical compliance 𝐶𝐶𝑡𝑡ℎ𝑒𝑒𝑒𝑒𝑒𝑒𝑒𝑒𝑒𝑒 Constraints: Each compliance C of the lattice structure is larger than 0(>1e-5) Variable: 𝑥𝑥 and 𝑦𝑦 coordinate of the start point 𝑀𝑀𝑀𝑀𝑀𝑀𝑀𝑀𝑀𝑀𝑀𝑀, 𝑀𝑀𝑀𝑀𝑀𝑀𝑀𝑀𝑀𝑀𝑀𝑀 (maximum, minimum allowable distance between the seed points) 86 the radius of the merged circle: 𝑅𝑅 spacing constant: 𝛼𝛼 the weighting factor for the mechanical field: 𝛾𝛾1 2nd level optimization: (with the optimized lattice layout) Objective: 𝐶𝐶𝑚𝑚𝑚𝑚𝑚𝑚ℎ𝑎𝑎𝑎𝑎𝑎𝑎𝑎𝑎𝑎𝑎𝑎𝑎 𝐶𝐶𝑡𝑡ℎ𝑒𝑒𝑒𝑒𝑒𝑒𝑒𝑒𝑒𝑒 Minimize compliance C= 𝛾𝛾1 + 𝛾𝛾2 (𝛾𝛾2 = 1 − 𝛾𝛾1 ) 𝐶𝐶𝑚𝑚_𝑜𝑜𝑜𝑜𝑜𝑜𝑜𝑜𝑜𝑜𝑜𝑜𝑙𝑙 𝐶𝐶𝑡𝑡_𝑜𝑜𝑜𝑜𝑜𝑜𝑜𝑜𝑜𝑜𝑜𝑜𝑜𝑜 𝐶𝐶𝑚𝑚_𝑜𝑜𝑜𝑜𝑜𝑜𝑜𝑜𝑜𝑜𝑜𝑜𝑜𝑜 , 𝐶𝐶𝑡𝑡_𝑜𝑜𝑜𝑜𝑜𝑜𝑜𝑜𝑜𝑜𝑜𝑜𝑜𝑜 are the optimal solution from the single optimization study Constraints: The total volume of the lattice structure is less than 10% 𝑉𝑉𝑡𝑡𝑡𝑡𝑡𝑡𝑡𝑡𝑡𝑡 Variable: Radii of lattice members: 𝑟𝑟𝑖𝑖𝑚𝑚𝑚𝑚𝑚𝑚 ≤ 𝑟𝑟𝑖𝑖 ≤ 𝑟𝑟𝑖𝑖𝑚𝑚𝑚𝑚𝑚𝑚 𝑖𝑖 = 1, 2, … 𝑡𝑡𝑡𝑡𝑡𝑡𝑡𝑡𝑡𝑡 𝑛𝑛𝑛𝑛𝑛𝑛𝑛𝑛𝑛𝑛𝑛𝑛 𝑜𝑜𝑜𝑜 𝑙𝑙𝑙𝑙𝑙𝑙𝑙𝑙𝑙𝑙𝑙𝑙𝑙𝑙 𝑚𝑚𝑚𝑚𝑚𝑚𝑚𝑚𝑚𝑚𝑚𝑚𝑚𝑚 where 𝑟𝑟𝑖𝑖𝑚𝑚𝑚𝑚𝑚𝑚 and 𝑟𝑟𝑖𝑖𝑚𝑚𝑚𝑚𝑚𝑚 represent the minimum and maximum allowable radius for the 𝑖𝑖𝑡𝑡ℎ lattice member 5.3.2.3 Weighting Factor Comparison Studies After testing the impact of different weighting factors on the layout of the lattice structure, the results are summarized in the following table of different combinations of weighting factors. (The study hasn’t been completed yet, some preliminary results are shown in this section) 87 Table 5.4: Lattice layout under different weighting factors The weighting factor for Optimized Lattice layout mechanical field 1 0 88 Chapter 6 Summary, Conclusions, and Future Work 6.1 Summary and Conclusions Chapter 3 A method to calculate orthotropic material properties in an AM process for 3D lattice structures has been developed and achieved using a bi-level design optimization approach. The optimization results have demonstrated the importance of accounting for material orthotropy within the lattice structure design and optimization process. Investigations on the effect of the printing direction have demonstrated that a feasible lattice made with the optimal printing direction could result in as much as a 20% reduction in volume compared to one made with the worst printing direction. Designs that ignored manufacturing- induced material orthotropy have been found to be unreliable as well. Further, since the principal material coordinates are often associated with the printing direction, a simultaneous optimization of the printing direction and lattice structure has been found to be important. Chapter 4 In this study, a lattice design methodology was developed that considers both material orthotropy and AM process resolution. Material orthotropy impacts the structural performance of the AM lattice structure and thus must be included in the design, while the minimum printing feature size of AM process influences the geometry of the final additively manufactured model. Both aspects are considered to produce and optimized as-manufactured lattice structure. 89 Chapter 5 Within the study in Chapter 5, we generate optimized lattice layouts under multiphysics loadings, with both mechanical and thermal loadings. The lattice layout algorithm combines multi-physics fields in a generalized way. This new field can then be used to generate the lattice layout. With the new algorithm, the new vector could be computed regardless of tuning for different boundary conditions. Based on different weighting factors for the stress and the heat conduction field, the trade-off of lattice layout between the different importance of multi-physics fields is presented. 6.2 Future Work With some preliminary results from multi-physics lattice optimization studies, future work would focus on the following studies: 1. Design and optimize lattice structures for different models with various loads, and boundary conditions, including heat convection. The current study for multi-physics fields lattice is only based on the loading conditions similar to former studies [5], [6] for validating the results. Further studies could apply the algorithm to more physics fields, including heat convection, thermal stress, etc. Different design domains could also be used. 2. Extend the lattice layout generation algorithm to 3D. Within the study, we only consider the 2D domain with all examples. To generate lattice structures for the 3D space, the algorithm should be extended to 3D space and consider more dimension combinations of multi-physics fields. 90 3. Collaborate with the geometric and material orthotropic manufacturability analysis within the multi-physics lattice design and optimization. Combining the manufacturability analysis within the design of multi-physics lattice structure could be investigated in the future. For multi-physics lattice structures, it is crucial to consider the manufacturing of lattice structures to be practical. Including both material orthotropy and AM process restrictions on multi-physics lattice design, the solution could be more realistic. 91 BIBLIOGRAPHY [1] M. Attaran, “The rise of 3-D printing: The advantages of additive manufacturing over traditional manufacturing,” Business Horizons, vol. 60, no. 5, pp. 677–688, Sep. 2017, doi: 10.1016/j.bushor.2017.05.011. [2] M. Scheffler and P. Colombo, Cellular Ceramics: Structure, Manufacturing, Properties and Applications. John Wiley & Sons, 2006. [3] S. Daynes, S. Feih, W. F. Lu, and J. Wei, “Optimisation of functionally graded lattice structures using isostatic lines,” Materials & Design, vol. 127, pp. 215–223, Aug. 2017, doi: 10.1016/j.matdes.2017.04.082. [4] S. Daynes, S. Feih, W. F. Lu, and J. Wei, “Design concepts for generating optimised lattice structures aligned with strain trajectories,” Computer Methods in Applied Mechanics and Engineering, vol. 354, pp. 689–705, Sep. 2019, doi: 10.1016/j.cma.2019.05.053. [5] Q. Gao, R. C. Averill, A. Diaz, K. Deb, and E. D. Goodman, “Towards Stress-Aligned Lattice Design,” submitted for publication 2022. [6] Z. Jiang, Q. Gao, R. C. Averill, K. Deb, and D. G. Erik, “Field-Aligned Lattice Design for Thermal Problems,” submitted for publication 2022. [7] W. P. Syam, W. Jianwei, B. Zhao, I. Maskery, W. Elmadih, and R. Leach, “Design and analysis of strut-based lattice structures for vibration isolation,” Precision Engineering, vol. 52, pp. 494–506, Apr. 2018, doi: 10.1016/j.precisioneng.2017.09.010. [8] J. Plocher and A. Panesar, “Review on design and structural optimisation in additive manufacturing: Towards next-generation lightweight structures,” Materials & Design, vol. 183, p. 108164, Dec. 2019, doi: 10.1016/j.matdes.2019.108164. [9] E. W. Hovig, A. S. Azar, F. Grytten, K. Sørby, and E. Andreassen, “Determination of Anisotropic Mechanical Properties for Materials Processed by Laser Powder Bed Fusion,” Advances in Materials Science and Engineering, vol. 2018, p. e7650303, Nov. 2018, doi: 10.1155/2018/7650303. [10] P. Maroti et al., “Printing orientation defines anisotropic mechanical properties in additive manufacturing of upper limb prosthetics,” Mater. Res. Express, vol. 6, no. 3, p. 035403, Dec. 2018, doi: 10.1088/2053-1591/aaf5a9. [11] J. Mueller and K. Shea, “The effect of build orientation on the mechanical properties in inkjet 3D-printing,” Aug. 2015. [12] W. Sun, Y. Ma, W. Huang, W. Zhang, and X. Qian, “Effects of build direction on tensile and fatigue performance of selective laser melting Ti6Al4V titanium alloy,” International Journal of Fatigue, vol. 130, p. 105260, Jan. 2020, doi: 10.1016/j.ijfatigue.2019.105260. 92 [13] R. Wauthle et al., “Effects of build orientation and heat treatment on the microstructure and mechanical properties of selective laser melted Ti6Al4V lattice structures,” Additive Manufacturing, vol. 5, pp. 77–84, Jan. 2015, doi: 10.1016/j.addma.2014.12.008. [14] S. Ahn, M. Montero, D. Odell, S. Roundy, and P. K. Wright, “Anisotropic material properties of fused deposition modeling ABS,” Rapid Prototyping Journal, vol. 8, no. 4, pp. 248–257, Jan. 2002, doi: 10.1108/13552540210441166. [15] P. Maroti et al., “Printing orientation defines anisotropic mechanical properties in additive manufacturing of upper limb prosthetics,” Mater. Res. Express, vol. 6, no. 3, p. 035403, Dec. 2018, doi: 10.1088/2053-1591/aaf5a9. [16] M. Rybachuk, C. Alice Mauger, T. Fiedler, and A. Öchsner, “Anisotropic mechanical properties of fused deposition modeled parts fabricated by using acrylonitrile butadiene styrene polymer,” Journal of Polymer Engineering, vol. 37, no. 7, pp. 699–706, Aug. 2017, doi: 10.1515/polyeng-2016-0263. [17] M. Liu, R. C. Averill, K. Deb, A. Diaz, and E. D. Goodman, “Effect of Manufacturing- induced Anisotropy on the Design of Additively Manufactured Lattice Structures”. [18] S. Ahn, M. Montero, D. Odell, S. Roundy, and P. K. Wright, “Anisotropic material properties of fused deposition modeling ABS,” Rapid Prototyping Journal, vol. 8, no. 4, pp. 248–257, Jan. 2002, doi: 10.1108/13552540210441166. [19] C. Ziemian, M. Sharma, and S. Ziemi, “Anisotropic Mechanical Properties of ABS Parts Fabricated by Fused Deposition Modelling,” in Mechanical Engineering, M. Gokcek, Ed. InTech, 2012. doi: 10.5772/34233. [20] M. Somireddy, A. Czekanski, and C. V. Singh, “Development of constitutive material model of 3D printed structure via FDM,” Materials Today Communications, vol. 15, pp. 143–152, Jun. 2018, doi: 10.1016/j.mtcomm.2018.03.004. [21] J. Cantrell et al., “Experimental Characterization of the Mechanical Properties of 3D-Printed ABS and Polycarbonate Parts,” p. 19. [22] B. E. Carroll, T. A. Palmer, and A. M. Beese, “Anisotropic tensile behavior of Ti–6Al–4V components fabricated with directed energy deposition additive manufacturing,” Acta Materialia, vol. 87, pp. 309–320, Apr. 2015, doi: 10.1016/j.actamat.2014.12.054. [23] W. Tao and M. C. Leu, “Design of lattice structure for additive manufacturing,” in 2016 International Symposium on Flexible Automation (ISFA), Aug. 2016, pp. 325–332. doi: 10.1109/ISFA.2016.7790182. [24] F. Tamburrino, S. Graziosi, and M. Bordegoni, “The Design Process of Additively Manufactured Mesoscale Lattice Structures: A Review,” J. Comput. Inf. Sci. Eng, vol. 18, no. 4, Dec. 2018, doi: 10.1115/1.4040131. 93 [25] S. Xu, J. Shen, S. Zhou, X. Huang, and Y. M. Xie, “Design of lattice structures with controlled anisotropy,” Materials & Design, vol. 93, pp. 443–447, Mar. 2016, doi: 10.1016/j.matdes.2016.01.007. [26] Y. Du, H. Li, Z. Luo, and Q. Tian, “Topological design optimization of lattice structures to maximize shear stiffness,” Advances in Engineering Software, vol. 112, pp. 211–221, Oct. 2017, doi: 10.1016/j.advengsoft.2017.04.011. [27] L. Cheng, J. Liu, X. Liang, and A. C. To, “Coupling lattice structure topology optimization with design-dependent feature evolution for additive manufactured heat conduction design,” Computer Methods in Applied Mechanics and Engineering, vol. 332, pp. 408–439, Apr. 2018, doi: 10.1016/j.cma.2017.12.024. [28] Y. Tang, A. Kurtz, and Y. F. Zhao, “Bidirectional Evolutionary Structural Optimization (BESO) based design method for lattice structure to be fabricated by additive manufacturing,” Computer-Aided Design, vol. 69, pp. 91–101, Dec. 2015, doi: 10.1016/j.cad.2015.06.001. [29] T. Stanković, J. Mueller, and K. Shea, “Optimization for Anisotropy in Additively Manufactured Lattice Structures,” in Volume 2A: 42nd Design Automation Conference, Charlotte, North Carolina, USA, Aug. 2016, p. V02AT03A027. doi: 10.1115/DETC2016- 59494. [30] T. Stanković, J. Mueller, and K. Shea, “The effect of anisotropy on the optimization of additively manufactured lattice structures,” Additive Manufacturing, vol. 17, pp. 67–76, Oct. 2017, doi: 10.1016/j.addma.2017.07.004. [31] H. Gonabadi, A. Yadav, and S. J. Bull, “The effect of processing parameters on the mechanical characteristics of PLA produced by a 3D FFF printer,” Int J Adv Manuf Technol, vol. 111, no. 3, pp. 695–709, Nov. 2020, doi: 10.1007/s00170-020-06138-4. [32] S. Tsopanos et al., “The Influence of Processing Parameters on the Mechanical Properties of Selectively Laser Melted Stainless Steel Microlattice Structures,” J. Manuf. Sci. Eng, vol. 132, no. 4, Aug. 2010, doi: 10.1115/1.4001743. [33] S. Chowdhury, K. Mhapsekar, and S. Anand, “Part Build Orientation Optimization and Neural Network-Based Geometry Compensation for Additive Manufacturing Process,” Journal of Manufacturing Science and Engineering, vol. 140, no. 3, Dec. 2017, doi: 10.1115/1.4038293. [34] A. Alafaghani and A. Qattawi, “Investigating the effect of fused deposition modeling processing parameters using Taguchi design of experiment method,” Journal of Manufacturing Processes, vol. 36, pp. 164–174, Dec. 2018, doi: 10.1016/j.jmapro.2018.09.025. [35] G. A. O. Adam and D. Zimmer, “On design for additive manufacturing: evaluating geometrical limitations,” Rapid Prototyping Journal, vol. 21, no. 6, pp. 662–670, Jan. 2015, doi: 10.1108/RPJ-06-2013-0060. 94 [36] X. Wang, M. Jiang, Z. Zhou, J. Gou, and D. Hui, “3D printing of polymer matrix composites: A review and prospective,” Composites Part B: Engineering, vol. 110, pp. 442– 458, Feb. 2017, doi: 10.1016/j.compositesb.2016.11.034. [37] A. Dey and N. Yodo, “A Systematic Survey of FDM Process Parameter Optimization and Their Influence on Part Characteristics,” Journal of Manufacturing and Materials Processing, vol. 3, no. 3, Art. no. 3, Sep. 2019, doi: 10.3390/jmmp3030064. [38] E. Ulu, N. G. Ulu, W. Hsiao, and S. Nelaturi, “Manufacturability Oriented Model Correction and Build Direction Optimization for Additive Manufacturing,” arXiv:1909.11620 [cs], Sep. 2019, Accessed: Nov. 02, 2019. [Online]. Available: http://arxiv.org/abs/1909.11620 [39] J. Mueller, K. Shea, and C. Daraio, “Mechanical properties of parts fabricated with inkjet 3D printing through efficient experimental design,” Materials & Design, vol. 86, pp. 902– 912, Dec. 2015, doi: 10.1016/j.matdes.2015.07.129. [40] W. E. Frazier, “Metal Additive Manufacturing: A Review,” J. of Materi Eng and Perform, vol. 23, no. 6, pp. 1917–1928, Jun. 2014, doi: 10.1007/s11665-014-0958-z. [41] M. K. Thompson et al., “Design for Additive Manufacturing: Trends, opportunities, considerations, and constraints,” CIRP Annals, vol. 65, no. 2, pp. 737–760, 2016, doi: 10.1016/j.cirp.2016.05.004. [42] L. J. Gibson, “Cellular Solids,” MRS Bulletin, vol. 28, no. 4, pp. 270–274, Apr. 2003, doi: 10.1557/mrs2003.79. [43] H. N. G. Wadley, “Multifunctional periodic cellular metals,” Phil. Trans. R. Soc. A., vol. 364, no. 1838, pp. 31–68, Jan. 2006, doi: 10.1098/rsta.2005.1697. [44] A. Hussein, L. Hao, C. Yan, R. Everson, and P. Young, “Advanced lattice support structures for metal additive manufacturing,” Journal of Materials Processing Technology, vol. 213, no. 7, pp. 1019–1026, Jul. 2013, doi: 10.1016/j.jmatprotec.2013.01.020. [45] S. Das and A. Sutradhar, “Multi-physics topology optimization of functionally graded controllable porous structures: Application to heat dissipating problems,” Materials & Design, vol. 193, p. 108775, Aug. 2020, doi: 10.1016/j.matdes.2020.108775. [46] I. Maskery, N. T. Aboulkhair, A. O. Aremu, C. J. Tuck, and I. A. Ashcroft, “Compressive failure modes and energy absorption in additively manufactured double gyroid lattices,” Additive Manufacturing, vol. 16, pp. 24–29, Aug. 2017, doi: 10.1016/j.addma.2017.04.003. [47] Z. Chen, Z. Wang, S. Zhou, J. Shao, and X. Wu, “Novel Negative Poisson’s Ratio Lattice Structures with Enhanced Stiffness and Energy Absorption Capacity,” Materials (Basel), vol. 11, no. 7, p. E1095, Jun. 2018, doi: 10.3390/ma11071095. [48] P. Heinl, L. Müller, C. Körner, R. F. Singer, and F. A. Müller, “Cellular Ti–6Al–4V structures with interconnected macro porosity for bone implants fabricated by selective 95 electron beam melting,” Acta Biomaterialia, vol. 4, no. 5, pp. 1536–1544, Sep. 2008, doi: 10.1016/j.actbio.2008.03.013. [49] C. Wang et al., “Multi-scale design and optimization for solid-lattice hybrid structures and their application to aerospace vehicle components,” Chinese Journal of Aeronautics, vol. 34, no. 5, pp. 386–398, May 2021, doi: 10.1016/j.cja.2020.08.015. [50] C. Pan, Y. Han, and J. Lu, “Design and Optimization of Lattice Structures: A Review,” Applied Sciences, vol. 10, no. 18, Art. no. 18, Jan. 2020, doi: 10.3390/app10186374. [51] D. Raymont, C. Yan, A. Hussein, and P. Young, “Design and additive manufacturing of cellular lattice structures,” presented at the Innovative developments in virtual and physical prototyping, Oct. 2011. doi: 10.1201/b11341-40. [52]J. Liu et al., “Current and future trends in topology optimization for additive manufacturing,” Struct Multidisc Optim, vol. 57, no. 6, pp. 2457–2483, Jun. 2018, doi: 10.1007/s00158-018- 1994-3. [53] J. Zhu, H. Zhou, C. Wang, L. Zhou, S. Yuan, and W. Zhang, “A review of topology optimization for additive manufacturing: Status and challenges,” Chinese Journal of Aeronautics, vol. 34, no. 1, pp. 91–110, Jan. 2021, doi: 10.1016/j.cja.2020.09.020. [54] S. Daynes, S. Feih, W. Lu, and J. Wei, Design concepts for generating optimised lattice structures aligned with strain trajectories, vol. 354. 2019. doi: 10.1016/j.cma.2019.05.053. [55] D. W. Kelly and M. W. Tosh, “Interpreting load paths and stress trajectories in elasticity,” Engineering Computations, vol. 17, no. 2, pp. 117–135, Jan. 2000, doi: 10.1108/02644400010313084. [56] A. Y. Tamijani, K. Gharibi, M. H. Kobayashi, and R. M. Kolonay, “Load paths visualization in plane elasticity using load function method,” International Journal of Solids and Structures, vol. 135, pp. 99–109, Mar. 2018, doi: 10.1016/j.ijsolstr.2017.11.013. [57] R. Weijermars, “Mapping stress trajectories and width of the stress-perturbation zone near a cylindrical wellbore,” International Journal of Rock Mechanics and Mining Sciences, vol. 64, pp. 148–159, Dec. 2013, doi: 10.1016/j.ijrmms.2013.08.019. [58] S. Catchpole-Smith, R. R. J. Sélo, A. W. Davis, I. A. Ashcroft, C. J. Tuck, and A. Clare, “Thermal conductivity of TPMS lattice structures manufactured via laser powder bed fusion,” Additive Manufacturing, vol. 30, p. 100846, Dec. 2019, doi: 10.1016/j.addma.2019.100846. [59] B. Vaissier, J.-P. Pernot, L. Chougrani, and P. Véron, “Parametric design of graded truss lattice structures for enhanced thermal dissipation,” Computer-Aided Design, vol. 115, pp. 1–12, Oct. 2019, doi: 10.1016/j.cad.2019.05.022. [60] K. Svanberg, “The method of moving asymptotes—a new method for structural optimization,” International Journal for Numerical Methods in Engineering, vol. 24, no. 2, pp. 359–373, 1987, doi: 10.1002/nme.1620240207. 96 [61] J. Deng, J. Yan, and G. Cheng, “Multi-objective concurrent topology optimization of thermoelastic structures composed of homogeneous porous material,” Struct Multidisc Optim, vol. 47, no. 4, pp. 583–597, Apr. 2013, doi: 10.1007/s00158-012-0849-6. [62] J. Yan, G. Cheng, and L. Liu, “A Uniform Optimum Material Based Model for Concurrent Optimization of Thermoelastic Structures and Materials,” Int. J. Simul. Multidisci. Des. Optim., vol. 2, no. 4, Art. no. 4, Dec. 2008, doi: 10.1051/ijsmdo/2008035. [63] E. M. Dede, “Multiphysics optimization, synthesis, and application of jet impingement target surfaces,” in 2010 12th IEEE Intersociety Conference on Thermal and Thermomechanical Phenomena in Electronic Systems, Jun. 2010, pp. 1–7. doi: 10.1109/ITHERM.2010.5501408. [64] R. M. Jones, Mechanics Of Composite Materials. CRC Press, 2018. [65] O. A. Bauchau and J. I. Craig, “Euler-Bernoulli beam theory,” in Structural Analysis, O. A. Bauchau and J. I. Craig, Eds. Dordrecht: Springer Netherlands, 2009, pp. 173–221. doi: 10.1007/978-90-481-2516-6_5. [66] “Abaqus 2016 Documentation.” http://130.149.89.49:2080/v2016/index.html (accessed Jul. 20, 2020). [67] “MATLAB - MathWorks.” https://www.mathworks.com/products/matlab.html (accessed Oct. 26, 2021). [68] J. Rodriguez, J. Thomas, and J. Renaud, “Mechanical behavior of acrylonitrile butadiene styrene fused deposition materials modeling,” Rapid Prototyping Journal, vol. 9, pp. 219– 230, Oct. 2003, doi: 10.1108/13552540310489604. [69] “HEEDS Design Exploration Software.” https://www.redcedartech.com/ (accessed Jul. 20, 2020). [70] “Abaqus Scripting Reference Guide (2016).” http://130.149.89.49:2080/v2016/books/ker/default.htm (accessed Jul. 20, 2020). [71] “Welcome to Python.org,” Python.org. https://www.python.org/doc/ (accessed Mar. 16, 2022). [72] “Simcenter HEEDS,” Siemens Digital Industries Software. https://www.plm.automation.siemens.com/global/en/products/simcenter/simcenter- heeds.html (accessed Jul. 20, 2020). [73] “MATLAB Documentation.” https://www.mathworks.com/help/matlab/ (accessed Mar. 16, 2022). 97