University of Southern California Dissertations and Theses

Deformable geometry design with controlled mechanical property based on 3D printing

(USC Thesis Other)

Deformable geometry design with controlled mechanical property based on 3D printing

PDF

Download a page range

Download transcript

Copy asset link

Request this asset

Transcript (if available)

Content University of Southern California Epstein Department of Industrial and System Engineering Deformable Geometry Design with Controlled Mechanical Property Based on 3D Printing by Yongqiang Li A Dissertation Presented to the FACULTY OF THE USC GRADUATE SCHOOL UNIVERSITY OF SOUTHERN CALIFORNIA In Partial Fulfillment of the Requirements for the Degree DOCTOR OF PHILOSOPHY (INDUSTRIAL AND SYSTEMS ENGINEERING) August 2014 Copyright 2014 Yongqiang Li i Acknowledgements First of all, I would like to acknowledge my advisor, Prof. Yong Chen. Without his insightful guide and persistent support I cannot complete my Ph.D. achievements. Second, I would like to say thanks to my wife, Lina Jiang. Without her persistent support, encouragement, tolerance and generosity of my anxiety from my pressure I cannot complete my Ph.D. work. Third, I would like to say thanks to my daughter Rebecca and Charlotte. Without their born I cannot enjoy the happiness of a father and have more courage to face the challenges in my Ph.D. work. In addition, I would like to say thanks to my parents, Zhonggao Li and CuiYing Peng; and my parents in law, Shuitao Jiang and Mei Kuang. Without their love and trust in a son, I cannot have the courage to complete my Ph.D. work. In the end, I would like to say thanks to my sister, Hongyong Li and my brother in law, Yunfei Chang. Without their care of my parents, I would have no enough time to complete my Ph.D. work. ii Table of Contents Acknowledgements ........................................................................................................................... i List of Tables ................................................................................................................................... iv List of Figures .................................................................................................................................. v Abstract ........................................................................................................................................ viii Chapter 1 Introduction .............................................................................................................. 1 1.1 Overview and Research Motivation ...................................................................................... 1 1.2 Statement of the Problem ....................................................................................................... 5 1.3 Background ............................................................................................................................ 7 1.4 Assumptions ......................................................................................................................... 15 Chapter 2 Design Deformable Geometries with Maximal Stiffness ....................................... 16 2.1 Problem Formulations and Our Approach ........................................................................... 16 2.2 Reduction and Growth Method on Principal Stress Line .................................................... 21 2.3 Principal Stress Line Computation ...................................................................................... 22 2.4 Domain Specification and Initial Structure Generation ....................................................... 25 2.5 Size Optimization of Beam Structures ................................................................................. 27 2.6 Topology Growth Based on Principal Stress Lines ............................................................. 30 2.7 Test Examples ...................................................................................................................... 42 2.8 Extensions ............................................................................................................................ 48 2.9 More Test Cases ................................................................................................................... 52 2.10 Summary ............................................................................................................................ 60 2.11 Appendix: Proof of Proposition 2.5.1.1 ............................................................................. 61 Chapter 3 Design Deformable Geometries with Target Stiffness ........................................... 65 3.1 Digital Material Cell and Digital Material ........................................................................... 65 3.2 Procedure of Digital Material Design .................................................................................. 73 3.3 Steps to Evaluate 2D Digital Material ................................................................................. 74 3.4 Steps to Interpolate 2D Digital Material .............................................................................. 77 3.5 Results of Digital Material Design ...................................................................................... 80 3.6 Case Study ......................................................................................................................... 103 3.7 Target Deformation Design Based on N Level stiffness Digital Material ......................... 104 3.7.1 Computational Framework for Desired Flexible Skin Design ................................... 104 3.7.2 Target Deformation Design Methodology ................................................................. 106 3.7.3 Case Study .................................................................................................................. 107 3.8 Summary ............................................................................................................................ 114 Chapter 4 Design Deformable Geometries with Target Compliance .................................... 116 4.1 Desired Compliance by Performance-tailored Meso-structures ........................................ 116 4.2 Formulation of Flexible Skin Design Problem .................................................................. 118 4.3 Computational Framework for Flexible Skin Design ........................................................ 120 iii 4.4 Optimization Methods for Flexible Skin Design ............................................................... 122 4.5 Flexibility Contribution Factors Methodology .................................................................. 125 4.6 Target displacements by Flexibility Contribution Factor Method ..................................... 130 4.7 Experimental Results and Analysis ................................................................................... 136 4.8 Summary ............................................................................................................................ 141 Chapter 5 Conclusion ............................................................................................................ 143 Bibliography ................................................................................................................................. 145 iv List of Tables Table 2.1: strain energy of growth cantilever at each iteration. ..................................................... 47 Table 2.2: strain energy of growth 3D cantilever at each iteration.. .............................................. 50 Table 2.3: strain energy of growth L-Shaped domain at each iteration... ...................................... 55 Table 2.4: strain energy of growth multi-loads bridge domain at each iteration... ........................ 57 Table 2.5: strain energy of growth convex L-Shape domain at each iteration............................... 57 Table 3.1: Evaluated Young’s modulus of all isotropic 2d digital cell patterns at level 1.............81 Table 3.2: Evaluated Young’s modulus of all anti-isotropic 2d digital cell patterns at level 1…..82 Table 3.3: Evaluated Young’s modulus of 2D digital cell with pattern 000000000 with varied count of cell rows…………………………………………………………………………………85 Table 3.4: Evaluated Young’s modulus of 2D digital cell with pattern 001011111 with varied count of cell rows…………………………………………………………………………………86 Table 3.5: 2D digital material basis patterns and their evaluated Young’s modulus…………….87 Table 3.6: Evaluated Young’s modulus of all interpolated patterns at the first span of level 1….88 Table 3.7: Evaluated Young’s modulus of all interpolated patterns from the second span of level 1………………………………………………………………………………………………......89 Table 3.8: Evaluated Young’s modulus of all interpolated patterns at the third span of level 1. ……………………………………………………………………………………………………90 Table 3.9: Evaluated Young’s modulus of all interpolated patterns at the fourth span of level 1. ……………………………………………………………………………………………………91 Table 3.10: Evaluated Young’s modulus of all interpolated patterns at the first span of level 3….. …………………………………………………..……………………………………………….93 Table 3.11: Isotropic 3D digital material basis patterns and their evaluated Young’s modulus…94 Table 3.12: transversely Isotropic 3D digital material basis patterns and their evaluated Young’s modulus………………………………………………………………………………………….. 96 Table 3.13: 3D digital material basis patterns and their evaluated Young’s modulus…………...96 Table 3.14: 3D isotropic material cell patterns – Level 2 from span 0 of BASIS 1 at Level 1….97 Table 3.15: 3D isotropic material cell patterns – Level 2 from span 1 of BASIS 1 at Level 1….98 Table 3.16: 3D isotropic material cell patterns – Level 2 from span 2 of BASIS 1 at Level 1. …………………………………………………………………………………………………..99 Table 3.17: 3D isotropic material cell patterns – Level 2 from span 3 of BASIS 1 at Level 1. …………………………………………………………………………………………………..100 Table 3.18: 3D isotropic material cell patterns – Level 2 from span 4 of BASIS 1 at Level 1. …………………………………………………………………………………………………..101 Table 3.19: Stiffness comparison of 2D digital material simulation result w/r physical experiment result……………………………………………………………………………………………..103 Table 4.1: Results of optimization methods (10-beam model)………………………………….133 Table 4.2: Results of optimization methods (bunny)……………………………………………141 v List of Figures Figure 2.22: L-Shaped domain…………………………………………………………………...54 Figure 2.23: The topology growth of an L-Shaped domain structure in 6 iterations…………….55 Figure 2.24: Multiple load bridge domain…………………………………………….…………57 Figure 2.25: The topology growth of multi-loads bridge domain in 5 iterations………………...57 Figure 2.26: convex L-Shaped domain…………………………………………………………..59 Figure 2.27: The topology growth of convex L-Shape domain in 7 iterations…………………..60 Figure 3.1: A 2D digital material cell……………………………………………………………65 Figure 3.2: Sequence ID of square elements in a 2D digital material cell……………………….66 Figure 3.3: Sequence of square elements in a mapped nine digit binary number………………..66 Figure 3.4: 2D digital cell 111101111……………………………………………………………66 Figure 3.5: 2D digital material bar FEM simulation setup with dimensions……………………..67 Figure 3.6: Target digital cell detail view………………………………………………………...68 Figure 3.7: A 3D digital material cell……………………………………………………………69 Figure 3.8: Sequence ID of cube elements in a 3D digital material cell…………………………70 Figure 3.9: Sequence ID of cube element groups in a mapped 4-digit binary number…………..70 Figure 3.10: Colored cube element groups as well as a mapped 4-digit binary number…………71 Figure 3.11: 2D digital material bar FEM simulation setup……………………………………...74 Figure 3.12: 2D digital cell at level zero…………………………………………………………75 Figure 3.14: 2D digital cell at level two………………………………………………………….76 Figure 3.15: 2D digital cell at level three………………………………………………………...76 Figure 3.16: 2D digital cell patterns at level zero………………………………………………...77 Figure 3.17: 2D digital cell patterns at level one…………………………………………………77 Figure 3.18: 2D digital cell patterns at level two…………………………………………………78 Figure 3.19: 2D digital cell 000000011 and its equivalents……………………………………...80 Figure 3.20: Sorted isotropic 2D digital cell patterns at level one with increment of stiffness…. 81 Figure 3.21: 2D digital cell with pattern 000000000……………………………………………..85 Figure 3.22: Graph of evaluated Young’s modulus with varied count of cell rows……………...86 Figure 3.24: Graph of evaluated Young’s modulus with varied count of cell rows……………...87 Figure 3.25: 2D digital material cell basis patterns with monotone increment of stiffness……...88 Figure 3.26: Soft and hard material patterns from the first span of level 1………………………88 Figure 3.27: Left: Graph of evaluated Young’s modulus of all interpolated patterns from the first span of level 1…………………………………………………………………………………….89 Figure 3.28: Soft and hard material patterns from the second span of level 1…………………...89 Figure 3.29: Left: Graph of evaluated Young’s modulus of all interpolated patterns from the second span of level 1…………………………………………………………………………….90 Figure 3.30: Soft and hard material patterns from the third span of level 1……………………...90 Figure 3.31: Left: Graph of evaluated Young’s modulus of all interpolated patterns from the third span of level 1…………………………………………………………………………………….91 Figure 3.32: Soft and hard material patterns from the fourth span of level 1…………………….91 Figure 3.33: Left: Graph of evaluated Young’s modulus of all interpolated patterns from the fourth span of level 1……………………………………………………………………………..92 Figure 3.34: Normalized graph of all four span curves at level 2………………………………..92 Figure 3.35: Soft and hard material patterns from the first span of level 2………………………93 Figure 3.36: Detail view of 2D cell patterns from the first span of level 2………………………93 vi Figure 3.37: Graph of evaluated Young’s modulus of interpolated patterns from the first span of level 3 ….………………………………………………………………………………………....94 Figure 3.38: Graph of evaluated Young’s modulus of all isotropic 3D cell patterns at level 1….95 Figure 3.39: 3D digital material cell basis patterns with monotone increment of stiffness……...97 Figure 3.40: Young’s modulus of all interpolated 3D cell patterns at the span 0 of level 2……..98 Figure 3.41: Young’s modulus of all interpolated 3D cell patterns at the span 1 of level 2……..99 Figure 3.42: Young’s modulus of all interpolated 3D cell patterns at the span 2 of level 2……100 Figure 3.43: Young’s modulus of all interpolated 3D cell patterns at the span 3 of level 2……101 Figure 3.44: Young’s modulus of all interpolated 3D cell patterns at the span 4 of level 2……102 Figure 3.45: Young’s modulus w/r 3D isotropic cell patterns at level 3………………………..102 Figure 3.46: stiffness comparison of 2D digital material simulation result w/r physical experiment result……………………………………………………………………………………………..103 Figure 3.47: Computational framework for desired flexible skin design……………………….105 Figure 3.48: Skin design domain with N-level printed digital material………………………...108 Figure 3.49: Target isosceles triangle like curve………………………………………………..109 Figure 3.50: Target isosceles triangle like curve is approximated by an isosceles triangle…….109 Figure 3.51: Approximated isosceles triangle is projected onto the skin design domain……….109 Figure 3.52: Stiff elements and design variables are assigned………………………………….110 Figure 3.53: Optimized material assignment result of isosceles triangle with given angle 150° . …………………………………………………………………………………………………..110 Figure 3.54: Optimized deformation result of isosceles triangle with given angle 150°………..111 Figure 3.55: Skin design domain with N-level printed digital material………………………...112 Figure 3.56: Target isosceles trapezoid like curve……………………………………………...112 Figure 3.57: Target isosceles trapezoid like curve is approximated by an isosceles trapezoid…112 Figure 3.58: Approximated isosceles trapezoid is projected onto the skin design domain……..113 Figure 3.59: Stiff elements and design variables are assigned………………………………….113 Figure 3.60: Optimized material assignment result of isosceles trapezoid with given angle 150° . …………………………………………………………………………………………………..114 Figure 3.61: Optimized deformation result of isosceles trapezoid with given angle 150°……...114 Figure 4.1: A feasibility study of achieving multiple design requirements……………………..117 Figure 4.2: A computational framework for flexible skin design………………………………120 Figure 4.3: An example of finite element analysis results generated by our automation tools…122 Figure 4.4: A test case of a 10-beam model……………………………………………………124 Figure 4.5: A beam example with negative contribution factor………………………………..124 Figure 4.6: Flexibility contribution factors of the 10-beam model……………………………..127 Figure 4.7: Flexibility contribution factors of the 10-beam model for different target points….127 Figure 4.8: Flexibility contribution factors of the 10-beam model for different load sizes……..129 Figure 4.9: Flexibility contribution factors of the 10-beam model for different link flexibility..129 Figure 4.10: Accumulated displacement with minimal and maximal spring ratios for the 10-beam model……………………………………………………………………………………………131 Figure 4.11: The optimized 10-beam models…………………………………………………...133 Figure 4.12: Optimized maximum displacement by GA and PSO……………………………...133 Figure 4.13: An example of computing a target displacement based on FCF and accumulated displacement…………………………………………………………………………………….135 Figure 4.14: Optimized distance from the target position for PSO……………………………..136 Figure 4.15: A complex test case: a bunny model………………………………………………137 Figure 4.16: Flexibility contribution factors for the bunny model……………………………...138 vii Figure 4.17: Flexibility contribution factors at a different target……………………………….138 Figure 4.18: Sorted flexibility contribution factors for the bunny model……………………….139 Figure 4.19: Accumulated displacement for the bunny model………………………………….139 Figure 4.20: The optimized bunny models by FCF method…………………………………….140 viii Abstract Deformable products exist broadly in different industries. These deformable products undergo deformation to achieve specific mechanical properties. Custom made insoles are applied to redistribute plantar pressure under feet to reduce stress related issues, such as diabetic ulceration. Also there are some other deformable products, i.e. cushioned chairs, custom made running shoes, light weighted space structure and so forth. The basic mechanical property under all these deformable products is stiffness with respect to deformation. Our research interest is to design these deformable geometries with controlled mechanical properties, and in order to do this, we investigated new design approaches and completed achievements to achieve mechanical property control in different research phases. First, we explored how to design complex geometries to achieve maximal stiffness in a single material. A PSL (principal stress line) based strategy is studied to design an optimal rigid Michell structure like topology. This Michell structure like topology obtained by PSL strategy is more refined than other approaches. We called this strategy reduction and growth approach. Second, we studied on how to design geometries to achieve target stiffness with multiple materials. Digital material design is investigated to design target stiffness in a given domain with a soft and a stiff material. Our digital material design method for desired stiffness is predictable and can be interpolated into more levels with higher resolution. Also we investigated how to apply these N-level digital materials to achieve controlled flexibility skin design. Simulation results for both digital material and flexibility skin design with N-level digital material have been shown and evaluated. A physical experiment was also investigated to verify the simulation results on digital material. ix Third, we investigated how to design geometries to achieve target compliance in a fixed topology of structure. A FCF (flexibility contribution factor) based method is provided to achieve the target compliance design of a geometry. This FCF method increases computation efficiency greatly compared to other optimal approaches. One challenge is how to fabricate these designed geometries with controlled stiffness,since they are either very complex geometries or contain multiple materials. Fortunately, 3D printing techniques have the unique capabilities to manufacture complex geometries and multiple materials. Nowadays, custom made products by 3D printing are going to main stream. We are expecting our research achievements to be applied with 3D printing and to impact the manufacturing industry. 1 Chapter 1 Introduction 1.1 Overview and Research Motivation Deformable products exist broadly in different industries. These deformable products undergo deformation to achieve specific mechanical properties. A cushioned chair is one of these deformable products to redistribute pressure, preventing a patient from sore or ulceration injuries (Chen and Wang 2008). A custom made insole is another example of these deformable products, which is employed to relieve the plantar peak pressure to reduce risks of stress related injuries, i.e. diabetic foot ulceration (Kä stenbauer, T. 1998). The basic mechanical property of all these deformable products is stiffness or compliance related to deformation. Our research interest is to design these deformable geometries with controlled mechanical property, i.e. maximum stiffness, target stiffness, target compliance and so forth. For the maximum stiffness design, people are interested in designing the most rigid topology of the structure given constant amount of material or using minimum amount of weight. Traditionally this maximum stiffness design problem is called topology optimization. Topology optimization is a classical subject in structural design. The study of fundamental properties of optimal grid-like continua was made by Michell (1904). In the last three decades, a considerable amount of work has been carried out on structural optimization. The two most well- known and developed structural optimization methods are , solid Isotropy Material with Penalization (SIMP), and Ground Truss Structure approach. More elaborated explanation of these methods are given in the following sections. Even though these two methods are well- established, limitations like checkboard issue (Bendsø e, M.P. and O. Sigmund 2002) and bars overlapping issue observed in our study exist. In order to address checkboard issue people presented filtering techniques. However, these filtering techniques prevent optimized structure 2 topology from being more refined to being an optimal Michell structure. Therefore, results from SIMP or ground truss structure strategy are not as elegant as the Michell type structure. Motivated by exploring optimization approach to design optimal Michell structure, we present an intuitive and heuristic approach. Our approach is based on principal stress line analysis of the given design domain. A comparison of our approach with the SIMP and ground truss structure methods is given in Figure 1.1. In evaluation, we believe our approach is more intuitive to understand and be able to obtain optimal Michell structure. It contains topology, geometry and size optimization. Strategy (a) SIMP (b) Ground Truss Structure (c) Principal Stress Line Step 1: Problem specification F F/4 F F/4 F F/4 Step 2: Design domain parameterization F F/4 F F/4 F F/4 Step 3: Optimization results F F/4 F F/4 Figure 1.1: A comparison of our PSL method with SIMP and Ground Truss Structure methods. One approach to design target stiffness is called N-level stiffness digital material design. Given two types of material, one is soft, another is stiff. We presented a method to design these N-level stiffness digital materials with key idea of cell pattern self propagation. These cell patterns consist of connected square elements are either filled with soft material or stiff material. With elements being filled with some cell pattern, the digital material shows target stiffness. Therefore, cluster of digital cells can be arrayed to form digital material. Figure 1.2 shows a digital cell and its corresponding digital material. The most significant advantage of this digital material is its predictability on design the pattern that achieves given stiffness. Another advantage 3 of our digital material design is its capability to propagate many levels of stiffness with high resolution. 33  Digital cell Digital material Hard material Soft material Figure 1.2: Digital cell and the digital material. To design target compliance, a wide variety of compliant mechanisms have been developed. The most well-known ones are springs and cantilever beams. They transfer certain amount of motion/energy from one point to another. However, most of these mechanisms are designed in the same size scale as the product component. Therefore, they can only be designed for limited types of displacements at limited number of positions. Figure 1.3:An example of using meso-structures to achieve heterogeneous material properties. In this thesis, we are considering a unique problem of using a big amount of compliant mechanisms in a smaller scale to achieve controllable flexibility in a product component. We call Ex1 Ey Ez X Y Z Ex2 Ex3 Ex4 E y = E z >> E x1 > E x2 > E x3 > E x4 R Y E Y R Z E Z L Spring L Inc r Spring R Spring E X 4 these compliant mechanisms in meso-scale level as meso-structures. They are the geometric arrangement of materials on a scale that is insignificant when compared to a product component in which they are used (Chen and Wang 2008). An example of such meso-structures (a set of springs and beams) is shown in Figure 1.3. Due to the large amount of meso-structures, a product component may achieve complex flexibility. That is, the component with meso-structures can behave differently indifferent directions (e.g. Ey and Ez >> Ex) and locations (e.g. Ex1> Ex2 > Ex3 > Ex4). Sufficiently large design space enable by the large amount of meso-structures may also allow us to accurately control such behaviors (i.e. achieving displacement at a given position in a given direction and amount under a loading condition). The characteristic length of our meso-structures is in the range of 1-10mm. It is mainly determined by a new type of manufacturing called Additive Manufacturing(AM) in 3D printing. Over the last twenty years, the AM processes have been developed based on the layer-based additive principle (Beaman, Barlow et al. 1996). Recent advances in material, process and machine development have enabled the AM process to evolve from prototype usage to final product manufacturing (rapid manufacturing) especially for aerospace and medical applications. The AM processes can directly fabricate parts from CAD models without part-specific tooling or fixtures. Consequently, truly complex 3-dimensional shapes in macro-and meso-scale levels (feature size > 0.1 mm) can be cost-effectively fabricated by them. This revolutionary capability allows complex meso-structures be designed with any configurations based on given design requirements. Therefore, in essence, we are using rigid materials and complex shapes that are enabled by the AM to achieve desired heterogeneous material properties, and consequently target flexibility, in a product component. 5 Motivated by observing deformable products applied broadly in industry and our daily life, my Ph.D. research is conducting deformable geometry design study with controlled mechanical properties. The deformable geometries designed in our study is not only for visualized in the computer or fabricated and put in a exhibit cabinet for show, but also is fabricated to achieve mechanical and functional property. In addition, the geometries are very complex and contain multiple materials potentially. Therefore, our study would generate more broad impact on the industry if these deformable geometries could be fabricated even they are complex and consist of multiple materials potentially. Fortunately, 3D printing techniques have these special capabilities to print complex geometries and multiple materials. In addition, 3D printing is going to main stream (Wohlers 2013). It is expected that my Ph.D. research is going to be applied in the manufacturing industry and generate broader impact. Our major achievements of this research are follows: (1) Investigated an intuitive method to achieve maximal stiffness design of geometries. We presented and developed reduction and growth method based on PSL(principal stress line) to achieve maximal stiffness design of geometries in a given domain; (2) Investigated design approaches to achieve target stiffness design of geometries. We presented and developed N-level stiffness digital material design approach with self propagation cell pattern to achieve controlled stiffness of a geometry with given multiple material; (3) Investigated a design approach to achieve target compliance design of geometries. We provided and developed FCF (flexibility contribution factor) approach to achieve target compliance design of meso-structure based geometries in single material. 1.2 Statement of the Problem In this thesis, we are going to investigate approaches to design the deformable geometries with controlled mechanical properties, with the achievements mentioned in Section 3. 6 Our work on deformable geometries design with controlled mechanical properties contains three phases in an ascending order by digging into the research problem more and more in-depth. The first one is deformable geometry design with maximal stiffness with given volume of a single material, i.e. minimum compliance problem or maximum stiffness problem. The second one is deformable geometry design with target stiffness with given multiple materials. And the third one is deformable geometry design with target flexibility with given tailored performance meso- structures in a single material. In this thesis, we are interested in how to distribute material in a given domain to design the Michel structure like topology to achieve the maximum stiffness or minimum compliance with given amount of material. Furthermore, we are interested in two types of deformable geometries design problems; one is to achieve deformable geometries design with target stiffness by investigating how to distribute multiple materials in a given domain; another is to investigate how to distribute meso-structures to achieve target compliance design in a given topology structure. Therefore, the deformable geometry design problem with respect to this thesis aims at achieving controlled mechanical property. In ascending order with respect to the depth of research, we are interested in addressing the design problem in three phases: Phase I: Deformable geometry design with maximal stiffness in a given domain. Given an amount of a single material, the problem in this phase is on how to find an optimal topology for material layout to maximize the stiffness in the domain with given constraints and loads. Phase II: Deformable geometry design with target stiffness by given multiple materials. This problem in this phase is on how to distribute multiple materials in a given domain to achieve deformable geometry design with target stiffness. 7 Phase III: Deformable geometry design with target compliance in a given topology of structures. The problem in this phase is on how to design deformable geometries by distributing meso-structures intelligently to achieve the target compliance in a given topology of structures. 1.3 Background Now-a-days 3D printing is going to mainstream(Wohlers 2013). 3D printing can fabricate custom products in small quantities economically (B. Berman 2012). More and more custom made products enter our daily life. Custom made sole and running shoes are being fabricated by 3D printing by New Balance, a sports manufacturer (Wohlers 2013). Team New Balance athlete, Jack Bolas, became the first ever track athlete to compete in custom made 3D printed plates, at the New Balance Games in January 2013 (Wohlers 2013). Customization of products by 3D printing is a trend in the future. Many injection molding products are designed and made with the identical model or limited types. These products are made from the same plastic or other materials. In practice when people sit down on a stool, chair, and cushion or run wearing running shoes; they have different needs to fit those products. Particularly these products come to deformation under stresses. The monotone or single material product can only provide people with boring deformation not fit every people. Straightforward, people would benefit from the custom made products to fit their needs better. A custom made insole is one of these customized products which can be fabricated by 3D printing.This custom made insole is an offloading device to relieve plantar pressure for avoiding and healing of plantar ulcers. There exists a causal relationship between diabetic foot ulceration and elevated plantar pressure (Kä stenbauer, T. 1998). And two approaches are provided to relieve plantar pressure for avoiding and healing of plantar ulcers. One is custom made insole, another is cushioned footwear. 8 Compared with planar insole, custom made insole can contact the foot totally. Researchers state that “total contact” insoles provide optimal offloading to release peak pressure (Bus SA 2004, Chen WP 2003, Lord M 1994). Chen and other authors employed three-dimensional finite element analysis to explore the effects of total contact insoles on the plantar pressure redistribution; and they found that two cases of total contact insoles in their study can both reduce high stresses at heel and metatarsal heads, redistributing the stress to the midfoot regions compared with the flat insoles (Chen WP 2003). Bus and other authors studied the effects of custom-made insoles on plantar pressures and load redistribution in neuropathic diabetic patients with foot deformity; and they found custom- made insoles decreased peak pressures and force-time integrals in the heel and first metatarsal head regions (Bus SA 2004). Owings and other authors studied to determine whether custom insoles tailored to contour of the barefoot pressure distribution and shape of a patient’s foot can reduce plantar pressures in the metatarsal head region to a bigger extent than conventional custom insoles. And they found the utilization of foot shape with barefoot plantar pressure measurement in designing custom insoles resulted in augmented off-loading of high-pressure regions under the forefoot(Owings 2008). A5513: “For diabetics only, multiple density insert, custom molded from model of patient’s foot, total contact with patient’s foot, including arch, base layer minimum of 3⁄16 inch material of shore A 35 durometer or higher, includes arch filler and other shaping material, custom fabricated”. This code is from the concern on the heterogeneity of insoles under Medicare coverage and their varied efficacies; and these specs stress the shape-based component of insole 9 design for a long-standing tenet of therapeutic shoemaking states that “total contact” insoles provided optimal offloading (Chen WP 2003, Bus SA 2004, Lord M 1994). However there is no standard definition of “total contact”. The current work only focus on the total contact before the insole is deformed. We believe there are research potentials on the total contact after the insole is deformed. Therefore, we believe the controlled flexibility skin design can contribute on custom insole design, or even can contribute the planar insole design. In recent years people are looking for more optimal way to customize the insole. Bus and other authors studied to assess metrics of in-shoe plantar pressure analysis to evaluate and optimize the pressure-reducing effects of diabetic therapeutic footwear. And their optimization is to modify the shoe or insole made from multidensity-layered materials (Bus SA 2011). Therefore it’s possible to employ flexible skin design to optimize an insole or even a shoe. Traditionally the custom made insoles were made from foam fabricated by CNC machining (Owings 2008). Bickel and other authors introduced a data-driven procedure for designing and fabricating multiple materials with desired deformation behavior; and they showed an insole example with complex heterogeneous materials fabricated by modern multi-material 3D printers (Bickel 2010). It is noticed that multiple materials play very important role in customized product by 3D printing, i.e. custom made insoles and custom printed optics (Willis 2012, Bickel 2010). Early in 2006, Hopkinson and other authors stated that there is potential to mix and grade materials in any combination that is desired, enabling materials with certain properties to be deposited where they are deserved with some of additive manufacturing processes(Hopkinson 2006). Composites are one category of multiple materials fabricated by 3D printing. They demonstrate superior mechanical properties to their base materials. Also 3D printing is a well- known technique to fabricate micrometer level of co-continuous structures or biological 10 composite topologies at micrometer resolution. It shows that one may employ computer geometric models to design superior mechanical properties of composites, such as fracture resistant composites, co-continuous composites materials with performance of stiffness, strength and energy dissipation (Leon 2013a, Leon 2013b, Wang 2011, Wang 2010). Topology optimization is a classical subject in structural design. Topology optimization of discretized and continuum structures are two broad categories in the structural optimization. We focus on discrete structure optimization in this paper. Introductions to truss topology problems can be found in (Topping, 1993; Kirsch, 1989; Kirsch, 1993; Bendsø e and Mota Soares, 1993b; Rozvany, et al., 1995a; Achtziger, 1997; Bendsø e and Sigmund, 2002). Hemp and Chan (Chan, 1960; Hemp, 1964; Hemp, 1973; Hemp and Chan, 1970), and in parallel efforts Dorn, Gomory and Greenberg (Dorn, et al., 1964) considered a ground structure to overcome the infeasibilities of Michell structures. With the development of numerical computational approaches, Dorn, Gomory and Greenberg were among the pioneers who developed a basis of a linear programming method (Dorn, et al., 1964; Haug and Arora, 1979). Given the same design domain, the analogous boundary conditions and external loads, they obtained the trusses coincide with the principal stresses directions of an optimal continuum structure (Achtziger, 1997). Nowadays the ground structure method is a well-known approach in the discrete topology optimization. Ground structure is composed of uniform spaced nodes connected with each other by boundary conditions and external loads or forces. The ground truss structure is thought to encompass the potential optimal structure. Introductions to ground structure approach can be found in (Topping, 1993; Achtziger, 1997; Bendsø e and Sigmund, 2002). The numerical computational theories on ground structure approach are mainly founded on minimization of compliance or maximization of stiffness. This objective function has been utilized in many literatures (Bendsø e and Sigmund, 2002; Achtziger, 1997; Bendsø e, 1995; Rozvany, et al., 11 1995a; Achtziger, et al., 1992; Bendsø e and Ben-Tal, 1993a; Svanberg, 1990; Svanberg, 1994; Taylor, 1969; Taylor and Rossow, 1977). In order to solve this objective function of minimization of compliance, linear or non linear programming techniques are developed (Achtziger, et al., 1992; Bendsø e and Ben-Tal, 1993a; Achtziger, 1997; Achtziger, et al., 2008; Achtziger and Stolpe, 2009). There are some other numerical computational approaches used to find the optimal truss structure from a ground truss structure (Hajela, et al., 1993; Xie and Steven, 1997). Node positions in a ground structure are to be optimized as well as topology and truss bar cross sectional size optimization. This further node positions optimization is called geometric approach(Topping, 1993). In geometric approach, the node coordinates are also considered to be variables as well as bar cross sectional size. Research on both topology and geometry of ground structures can be found in (Ben-Tal, et al., 1993; Achtziger, 2007). The complexity of ground structure approach is ) ( 2 n O , where n is the number of the nodes. When n is large, the ground structure is very dense; and the numerical computation of LP or nonlinear LP techniques is unstable; and some unreasonable structures are obtained (Burns, 2002). Rozvany, et al. (Rozvany, et al. 1995b) noticed that many literatures focused on truss structures rather than beam structures and pointed that the beam structures are more practical than truss structures. For truss structures pin-jointed by only considering the axial forces and buckling is unavoidable. Beam structures are rigid-jointed and considering the shear forces which are matched with the engineering purpose. Theories on beam structure optimization can be found in (Rozvany, 1994; Rozvany and Prager, 1976). However, these theories are targeted for the continuum or partial continuum beam structure optimization. They seem to be not simple to be applied in numerical computation. Some people combined the ground structure approach with 12 the continuum based material optimization method to design the beam structures with joints (Fredricson, et al., 2003; Fredricson, et al., 2003; Fredricson, 2005; Kim, et al., 2008). These beam structures with joints are optimized with their sizes and can be manufactured directly. They are called compliant assemblies. Basically their physical model is a 3D solid and either sensitivity analysis or SIMP (simple isotropic material with penalization (Bendsø e and Sigmund, 1999) are utilized to optimize the size of each beam. Their optimization methods are continuum based. Not only can composites achieve heterogeneous material properties, but also meso- structures can achieve heterogeneous material properties. Then another approach is to employ meso-structures to achieve heterogeneous material properties. The usage of meso-structures or different types of cellular structures in a product component is not new. The benefits of component design with cellular structures have been demonstrated in a wide variety of structural applications (Gibson and Ashby 1997), thermal applications (McDowell 1998) and biomedical applications (Ma and Elisseeff 2005). Some examples of such applications are shown in Figure 1.4. These structures are innovative periodic cellular materials that derive their outstanding mechanical and other physical performance from a structure of highly ordered unit cells rather than the properties of the parent materials. There are a lot of desired properties that can be achieved by cellular materials include high specific stiffness, high capacity for kinetic energy absorption, excellent vibration absorption and damping characteristics (Gibson and Ashby 1997). 13 Figure 1.4: Examples of cellular structure in various applications. For desired properties, the design of such cellular structures is highly related to the manufacturing processes that will be used to fabricate them. For example, most metal foams are made by stochastic manufacturing processes such as metal gas injection. Due to the lack of accurate control in those processes, the accurate shape modeling of cellular structures is not required or necessary. As another example, most corrugated panels or honeycomb structures are made by sheet metal forming processes (refer to Figure 1.4.a-bottom). Consequently, limited variations can be achieved in the corrugation configurations. Therefore, most studies on cellular structure design (Bendsoe and Kikuchi 1988; Bendsoe 1995; Gibson and Ashby 1997; Sigmund 1997; Sigmund 2000) are focused on the topology design of a single cellular structure and its properties. Based on such properties, the macro-level component property is studied by assuming a component has only a single type of cellular structures and all of them are identical. In comparison, we focus on the problem of how to use a combination of pre-defined meso-structures to achieve desired properties. In addition to multiple types of meso-structures that can be used in a component, each meso-structure can be tailored with a unique size. (a) Structural Applications (top: grc.nasa.gov; bottom: Kalpakjian & Schmid 2006) (b) Thermal Applications (McDowell 1998) (c) Biomedical Applications (www.microfab.com) 14 Another type of mesoscopic structures is the compliant structure, in which some elements are designed to be intentionally flexible. This flexibility allows the structure to function as a mechanism (i.e. compliant mechanism) (Howell 2001). The design of compliant mechanisms has been intensively studied. Two most popular design methods are the pseudo-rigid-body model (Jensen and Howell 2003; Jensen and Howell 2004) and topology optimization (Sigmund 1997; Joo, Kota et al. 2000; Joo and Kota 2004; Lu and Kota 2006). In the pseudo rigid-body model method, a flexible beam is approximated by a rigid, rotating beam and a torsional spring placed at the pin joint. The topology optimization method is more general. It studies how to distribute materials in a design domain such that an objective function is optimized. Some well-known approaches, such as SIMP and ground structure methods have been proposed for designing compliant mechanisms. The heterogeneous material properties of meso-structures can assist to design flexible skin. The opposite extremity of flexible skin design is to design a rigid structure with minimum compliance. It is observed that custom made insole or running shoes are deformable products. Motivated by designing these deformable geometries and fabricating them by 3D printing techniques, we are conducting this Ph.D. research study to design deformable geometries with controlled stiffness. For maximal stiffness geometry design, motivated by simplifying the complexity of ground structure and find the more Michell structure like topology structure, we proposed a principal stress line method for the optimal structure design. We also extended uniform strain energy density from truss structure to beam structure topology optimization. And for target stiffness geometry design, motivated by 3D printing’s capability to print multiple materials we study to design varied stiffness digital material with pattern predictable to meet the potential heterogeneous material distribution in the flexible skin design. 15 In addition, for target compliance geometry design, motivated by 3D printing’s capability to print complex geometries we study how to distribute meso-structures in a constant topology structure to achieve target compliance. 1.4 Assumptions In this research we have the following assumptions: (1) We assume all deformation in our study occurs in range of elasticity. If there were no special explanation, the deformation occurs in linear range of elasticity. (2) In our study on the component design with meso-structures, we assume all the rigid beams have the same size (we set r = 0.65mm in our study) since changing their sizes has a relatively small effect on the flexibility of a component’s skin. (3) In our FEA model, we did not consider extremity conditions such as yield failure, fracture failure or buckling. (4) In the beam structure study we assume a beam is an Euler-Bernuolli beam. An Euler- Bernuolli beam assumes the length of a beam if far bigger than dimension on the other two dimensions. 16 Chapter 5 Design Deformable Geometries with Maximal Stiffness Currently we provided and investigated reduction and growth method based on PSL(principal stress line) method to achieve maximal stiffness geometry design. This topology optimization approach is to mimic principal stress line to build and refine the optimized topology. The benefits of component design with cellular structures have been demonstrated in a wide variety of applications. The recent advances in additive manufacturing and high performance computing have enabled us to design a product component with adaptive cellular structures to achieve significantly better performance. However, designing a product component with such structures, especially its shape topology, poses significant challenges. Many approaches in topology optimization have been developed before for the purpose. In this proposal, we present a novel topology generation approach based on the principal stress line analysis for structural applications. We first discuss the computation of principal stress lines in a given design domain. Based on some well-known principles in biomimicry research, we then present a novel topology design method based on principal stress line analysis. Related algorithms are given for generating a shape topology that is adaptive to the given functional and shape requirements. Test cases are also presented to illustrate the presented method. 2.1 Problem Formulations and Our Approach 2.1.1 Problem Formulation The optimization problem we consider in the proposal is to minimize the compliance on the total structural volume of a beam structure under static loads and constraints with the cross- sectional areas, the nodal coordinates, and the beam connections as design variables. Hence all the three types of optimization problems are considered including (1) size optimization by 17 changing cross-sectional areas; (2) geometry optimization by changing the nodal coordinates; and (3) topology optimization by changing the beam connections. In form of finite element method the minimum of compliance problem can be written as in equation(2.1) (Eschenauer and Olhoff, 2001): ) ( min u U B χ  subject to , F Ku  (2.1) Where B is a specific volume of body;  denotes the design variable; U is the strain energy; K is the stiffness matrix; u is the displacements of all nodes; F is the external loads. 2.1.2 Michell Type continua and Principal Stress Lines Michell (1904) first studies the fundamental properties of optimal grid-like continua and contribute exact analytical solutions for several well-known truss structures. One of them is Michell-cantilever. Modern layout theory was founded by Prager and Rozvany (Prager and Rozvany 1977). Rozvany, et al. (Rozvany, et al. 1995a) states that Michell trusses have the maximum stiffness for a given volume. According to Strang and Kohn (Strang and kohn, 1983), the Michell problem is defined as the following equation: 0   div in  , f n  :  on  (2.2) The stress tensor is represented by the symmetric matrix as:        y xy xy x      , (2.3) and the principal stresses (eigenvalues) are ) 4 ) ( ( 2 1 2 2 2 , 1 xy y x y x            (2.4) Hence, one can imagine that the bars of a truss-like continuum are placed in the directions of principal stresses. That is, the Michell structure composes of infinite number of truss bars 18 perpendicular to each other at the intersection (Strang and Kohn, 1983; Prager, 1974). And the slope of the principal stresses with respect to the x-axis is given in (Hill, 1950) as: y x xy       2 2 tan (2.5) From (Pilkey, 2002), we have the stress strain relation as:                                       xy y x xy y x v v v v v v v E       2 2 1 0 0 0 1 0 1 ) 2 1 )( 1 ( (2.6) where E is the elastic or Young’s modulus; v is Poisson’s Ratio; x  , y  are strain in x and y direction respectively; xy  is shear stress in the xy plane. Substitute (2.5) using (2.6) we have y x xy       2 tan (2.7) Hence the principal stress lines(PSL) are defined as the two orthogonal families of curves whose directions at every point coincide with those of the minimum shear stress in a state of plane strain. Our approach is to mimic the PSL to obtain the topology of discrete structure of beams and geometry of nodes to solve the Michell type continua to achieve the optimal structure design. Let us discuss some properties of PSL first. 2.1.3 Principal Stress Line Analysis for Structural Optimization Since our study is focused on the linear deformation, from the conditions of elastic equilibrium, given the design domain, external forces and boundary conditions, in a finite element system of isotropic material we can show that the following properties hold obviously: Property 2.1.3.1 The displacement vector u is proportional to the external force vector F or stiffness matrix K . Property 2.1.3.2 The displacement vectoru is inversely proportional to the elastic modulus E. 19 From (Pilkey, 2002), we have y u x u y u x u x y xy y y x x                , , (2.8) Substitute (2.8) into (2.7) we have: y u x u x u y u y x y x             2 tan (2.9) Suppose 1 C is the scaling factor of the initial external force vector 0 F ; 2 C is the scaling factor of the initial elastic modulus 0 E . Then from 0 0 0 ,  K u F we have: , 0 1 0 2 F u K C C  (2.10) Compare it with 0 0 0 F u K  , we know: , ) ( 0 1 2 u u C C  (2.11) Then we have: y u x u x u y u y x y x            0 0 0 0 2 tan  (2.12) Hence, we have the following property. Property 2.1.3.3 The direction of the principal stress is not related to the scaling of the external forces, nor the material type for an isotropic material within the range of elastic deformation. As an example, we present the principal stress field of a simple cantilever beam as shown in Figure 2.1. For a single load with various sizes and different material properties, the principal stress lines are virtually the same. This suggests that their optimal topology and shape design should be identical, which is in accordance to the common sense in structural optimization. 20 (a) Steel, E=2.0e11(Pa), F=10(N) (b) Steel, E=2.0e11(Pa), F=100(N) (c) Aluminium, E=7.0e10(Pa), F=10(N) (d) Aluminium, E=7.0e10(Pa), F=100(N) Figure 2.1: The principal stress fields varied loads and material properties. Property 2.1.3.4 The principal stress field is mainly related to the topologic properties of the given structural design such as the position of external forces and the types of constraints. As an example, we present the principal stress field for a simple cantilever beam with different load and constraint positions. As shown in Figure 2.2(a), the external force F is moved from the middle to the bottom of the beam. As shown in Figure 2.2(b), two point constraints are used instead of a line constraints in Figure 2.2. Hence their optimal topology and shape design should be different, which is also in accordance to the common sense in structural optimization. 21 (a) Different load position (b) Different types of constraints Figure 2.2: Principal stress lines distribution in the design domain of the continuum domain. 2.2 Reduction and Growth Method on Principal Stress Line The principal stress line method is designed for solving the Michell type continua thru the minimum compliance or maximum stiffness of a beam structure under any given loads and constraints. The optimization of size, topology, and geometry of beam structure are all considered in our method. Not like SIMP and ground truss structure approach, our PSL structural optimization strategy does not distinct multiple load cases from single load case. For multiple load cases, our PSL strategy only loads all of them at one time. Therefore, our PSL strategy addresses those multiple load cases by single load strategy. We are not going to mention this multiple load case addressing strategy any more in the following discussion in our paper. If there is some case with more than one load appearing in our paper, this case could be either single load case or multiple load case. Based on a single loaded Michell type structure, the design process of the principal stress line method is shown in Figure 2.3. There are five major steps: (1) design domain specification, (2) principal stress lines computation, (3) initial structure generation, (4) size optimization, and (5) topology growth. 22 Section 7.2.4 Section 7.2.3 Section 7.2.6 Figure 2.3: An overview of the principal stress line method for beam structure design. For a given design domain specified by a user, the stress lines are computed and visualized based on the specified loads and constraints, the material information, and the maximum space that the resulting structure can occupy. An initial structure is then generated that connects all the constraints and loads. Based on the given beam structure, size optimization is then carried out, in which the cross section area of each beam is modified. Based on the principal stress lines and the designed beam structure, more principal stress lines will be identified in topology growth. The added principal stress lines will accordingly add more nodes to the beam structure and also modify the positions of some existing nodes(existing nodes positions are fixed from our methods, only the straight line shape are morphed with more principal stress line being added). Hence a new beam structure with different topology and shape will be generated to improve the structure performance. The process is repeated until no performance improvement is found or a limit on number of beams has been reached. Section 2.3~2.6 explain in detail each of the five steps in the above process. 2.3 Principal Stress Line Computation 23 Hegemier and Prager (Hegemier and Prager, 1969) utilizes the mathematical equations to compute the Michell truss bars based on the analytical computation of the directions of principal strains for a single loaded truss structure. The similar complementary slip lines are computed analytically in (Hill, 1950). In this paper, we develop a general numerical computation method to compute the directions of the principal stress lines. In addition to simple cases with a single load such as Michell truss bars, our numerical computation method can handle general cases with more complex boundary conditions and various external loads. 2.3.1 Computing Principal Stresses at Any Point We first analyze the given design domain with loads and constraints by using the Finite- Element Analysis (FEA) method. As shown in equation 2.3 and equation 2.4, the principal stresses and directions are well defined for any given point in the stress field of a continuum 2D solid. However, in the 2D plane static analysis results generated by a FEA system (we used COMSOL3.2: www.comsol.com), the finite elements are a set of triangles and all the post- processing data of normal stresses x  , y  and shear stress xy  are corresponding to triangle vertices. In order to query the post-processing data in any position in the design space, we compute all the stresses based on linear interpolation. That is, suppose the stresses at position P is needed and assume P is located in a triangle element 3 2 1 P P P  . The corresponding post- processing data (i.e. x  , y  or xy  ) for 1 P , 2 P and 3 P are 1 D , 2 D and 3 D separately. The post- processing data D at the position P can be interpolated as: 3 2 1 wD vD uD D    (2.13) where 3 1 2 1 3 2        P P P P PP PP u , 3 1 2 1 1 3        P P P P PP PP v , 3 1 2 1 2 1        P P P P PP PP w . Hence the principal stresses 1  , 2  and corresponding directions parameterized by  at P can be computed. 24 2.3.2 Generating 2D Principal Stress Lines in a Continuum Domain Based on the principal stresses computed for a given point in the design domain, we will have four principal directions (in  , 90   , 180   , and 270   ). Hence four principal stress lines can be generated beginning from a stress field point with the four principal directions. For simplicity and without loss of generality, we will discuss the principal stress line search strategy in one direction. As shown in Figure 2.4(a), from point 1 i P  , we can search one of its principal stress direction 1 i v  and compute point i P for a given approximation step  . Accordingly the resulted principal stress line is shown in Figure 2.4(b). Figure 2.4: principal stress line search based on stress field points. Hence, we have the following iterative equations:           1 0 , 1 1 0 0 i v P P i v P i i i    (2.14), Where 0 0 ,v P  are the starting position and the principal direction;  is the search step; i v  is one of the four principal directions at position i P closest to the search direction 1  i v  . The closest direction selection is finished by comparing the angle with the search direction 1  i v  . The inner product is used to do the comparing. Sufficient number of starting points is needed for generating principal stress lines that are dense enough to cover the entire design domain. We used two strategies in generating principal stress lines that are uniformly distributed in the design space. 1  i P i P 1  i P  1  i v  i v  1  i P i P 1  i P (a) (b) 25 (1) For simple design space, we add a guide curve and sample the curve uniformly. The sampling points are used as the starting points to generate the principal stress lines. (2) For complex design space, we used the design space boundary as the guide curve and sample curve uniformly. The sampling points are used as the starting points to generate the principal stress lines. Some examples of the accordingly generated stress lines are shown in Figure 2.1 and 2.2. 2.4 Domain Specification and Initial Structure Generation As shown in Figure 2.3, our method starts from the user specifying the loads and constraints, and also the maximum domain space that the resulting structure can occupy. Based on the computed principal stress lines, the user can identify the principal stress lines that connect the constraints and loads. We called such lines skeleton principal stress lines. An initial structure can then be generated from the skeleton principal stress lines. The generation of the initial structure can be performed automatically or by the user. For example, Figure 2.5 left shows the design domain for the Michell cantilever. The physics properties we used in the test are: H = 1.0 m; F = 1000 N; Young’s modulus E = 2.0e11 Pa; Poisson’s ratio  = 0.33; Density = 7850 3 / kg m ; Shape of a beam’s cross section = square; Given volume of the material = 0.015 3 m . Figure 2.5 right shows the computed principal stress lines. The skeleton principal stress lines are also shown to connect the load and constraints. Accordingly an initial structure is generated. 26 L=1.5H H F Figure 2.5: The design domain and initial structure of a Michell cantilever beam case. Notice after a candidate principal stress line is selected, we use a straight line connecting its two nodes in the beam structure instead of approximating the exact stress line by a set of polylines. This is due to the following theorem we called zero shear stress principle. Theorem 2.4.1: For two given nodes in a beam structure, a straight line connecting them will lead to a smaller compliance of the structure than any other curves to connect them. Proof: There are two main reasons for a straight line giving the minimum compliance of the structure: (1) A straight line segment between two adjacent nodes has the shortest distance. Therefore the consumed material of the beam is the least. For a constant volume of material given for the whole beam structure, consequently, the stiffness of the beam structure can be bigger for a straight line. (2) More importantly, the stiffness of a beam is dominated in the axial stress direction and the stiffness of a beam in the shear stress direction considerably smaller than the stiffness in the axial stress (refer to Section 2.5). Therefore the shear stress of a beam should be Skeleton Principal Stress Lines Initial Beam Structure 1 F 2 F 1 P 2 P 3 F Figure 2.6: Shear force on a curve. 27 reduced in order to achieve a smaller compliance. Suppose two adjacent nodes 1 P and 2 P are connected by a curve instead of a straight line segment as shown in Figure 2.6. This connection beam is either compressed or stretched. Without loss of generality, we assume the connection to be compressed. Hence two axial forces of 1 F and 2 F are along the tangential directions on both nodes of 1 P and 2 P . In order to achieve to the equilibrium of this connection beam, a force 3 F is generated by the internal strain deformation of the connection beam. Obviously 3 F is a shear force. 3 F will be zero if and only if the connection is a straight line. □ Hence in all the iterations during the topology growth (refer to Section 2.6), all the principal stress line will be approximated by straight lines between identified nodes. 2.5 Size Optimization of Beam Structures Once the topology of the optimized beam structure is fixed in each iteration in the growth process, the next step in our PSL strategy is to optimize the cross section size of a beam. In our beam structure design, we follow energy principles which present the basis size optimization of discrete and continuum structures. 2.5.1. Theoretical Basis – Axiom of Uniform Strain Energy Density Based on the assumption of the positive definiteness of the global stiffness matrix, Achtziger (Achtziger, 1997) proved that the strain energy density is uniform among the truss structure. However, no references were found on beam structures. Here we will extend some conclusions of topology optimization of truss structures in (Achtziger, 1997) to the topology optimization of beam structures. In order to employ structural optimization conclusions from truss structures we following two propositions For beam structures. 28 Proposition 2.5.1.1 Let i x is the volume of i-th beam in a planar beam structure, then the element stiffness matrix of the i-th beam can be written as 2 2 1 i i i i d i x x K K K   ; both 1 i K and 2 i K are positive semi-definite matrices. Under the classic Euler-Bernuolli assumption, the stiffness of a beam is dominated by 1 i K in the axial direction and we have i i d i x K K  (proof is given in Appendix B). Since the volume x has a linear relation to the stiffness matrix K , we know u x K u u u x ) ( ) , ( 2 1 T T f   is convex in x . Proposition 2.5.1.2 Let i i x K is the i-th beam’s element stiffness matrix in a truss or beam structure, then i K is a positive semi-definite matrix. Proof: Because i i T U  u K u 2 1 is the strain energy density of beam i for u  . 0  i U when the entire beam structure is within the elastic deformation range. According to the definition of a positive semi-definite matrix, we know stiffness matrix i K is a positive semi-definite matrix. □ From the optimum conditions for a truss structure in (Achtziger, 1997) we know:   ) ( ) ( max ) ( ) ( : 2 1 1 2 1 x u K x u x u K x u x x j T m j i T i new i    (2.15) For all m i , , 1   beams, ) ( ) ( 2 1 x u K x u i T is the strain energy density. From (2.15), it is easy to see the strain energy density is uniform when a beam structure achieves the equilibrium. 2.5.2. Computation Approach As shown in Section 2.5.1, the objective in size optimization is to achieve the maximal stiffness of the beam structure. Based on the method of uniform strain energy density, to 29 maximize stiffness also means to minimize the strain energy of the beam structures. In addition, the beam structure achieves the minimal potential energy when the strain energy density is uniform in the entire beam structure. We can define the strain energy density as strain energy per unit volume of a beam, denoted by  with unit of 3 / m J . Hence we can optimize the cross section size of each beam for a given beam structure similar to the Fully Stressed Design Stress-Ratio Method in truss structures (Topping, 1993). However, instead of stress ratio, we use the strain energy density ratio to modify the cross section area of beam members after each stress strain analysis of beam structures: 2 1 1 } { max             i m i i iold inew A A   (2.16) Notice we used the square root of the strain energy density ratio instead of the one order ratio in (2.16). This is because in our experiments, we found the one order ratio was not convergent well but square root of the ratio works very well. In stress-ratio method (Topping, 1993), it was found that the cross-section areas of some of truss bars will reduce to zero in a dense connected ground structure. In our study, the strain energy density ratio method will have the similar problem. In order to guarantee the numerical computational stability, we utilize the following strategy. Consider the ratio of the cross-section of a beam over the maximal cross-section of all beams: } { max 1 i m i i i A A     (2.17) 30 When the cross-section area ratio is small enough, i.e. threshold i    the cross-section area of the i-th beam is assigned a small number. In our study 6 10   e threshold  , and once a cross section area ratio is lower than this threshold value, a new value of area is assigned to this beam: } { max 1 i m i threshold i A A      (2.18) A beam with cross section area ratio lower than the threshold value is called a zero cross section beam. After cross section areas of all beams are resized by the strain energy density ratio method, all cross section areas of beams are scaled to fit the given volume V X (note zero cross section beams still have the threshold cross section area ratio). With the resized beam structure a new FEM stress strain calculation is carried out, in which the cross section area in the bigger strain energy density beams gets bigger and bigger( in a relative way), and the cross section area in the smaller strain energy density beams gets smaller and smaller. Steps 2 and 3 are repeated again and again, the strain energy distribution in the whole beam structure gets more and more uniform. This iteration process terminates when the strain energy of the entire beam structure reaches a plateau. Figure 2.7 shows an example of how size optimization works. (a) A beam structure with uniform cross section; (b) resulted beam structure after size optimization Figure 2.7: An test example of size optimization based on uniform strain energy density. 2.6 Topology Growth Based on Principal Stress Lines The theoretic Michell structure has infinite number of bars and each bar is infinitely small. Hence such exact analytical solution is not practical for a beam structure design. Instead we 31 would like to balance the number of beams and the performance of beam structure. For the purpose, we present a topology growth method based on the principal stress lines. From an initial structure, the main procedure of our method is: (1) Identify new principal stress lines that can reduce the approximation errors the most between the beam structure and the principal stress lines; (2) Use the identified principal stress lines to compute a set of intersection points as the nodes of beam structure (shape optimization); (3) Construct a candidate beam structure by connecting the computed nodes based on the principal stress lines (topology optimization). (4) Optimize the cross section size of each beam for the constructed beam structure (size optimization). The strain energy of the whole beam structure can be computed and utilized to measure the stiffness of the beam structure. The above procedure can be iterated until desired performance is reached or the maximum number of beams is exceeded. This is illustrated in Figure 2.8 based on a Michell cantilever beam case as shown in Figure 2.5. Iteration 0, strain energy= 4.1493mJ Uniform strain energy density distribution 32 Iteration 1, strain energy= 3.800 mJ Uniform strain energy density distribution Relation of Strain Energy and Iteration Times 0.00E+00 1.00E-03 2.00E-03 3.00E-03 4.00E-03 5.00E-03 6.00E-03 1 3 5 7 9 11 13 15 17 19 21 23 25 27 29 31 Iteration time Strain energy [J] Relation of Strain Energy and Iteration Times 0.00000000 0.00100000 0.00200000 0.00300000 0.00400000 0.00500000 0.00600000 1 4 7 10 13 16 19 22 25 28 31 34 37 40 43 46 49 52 55 58 Iteration time Strain energy [J] 33 Iteration 2, strain energy= 3.471 mJ Uniform strain energy density distribution Iteration 3, strain energy= 3.445 mJ Uniform strain energy density distribution Relation of Strain Energy and Iteration Times 0.00000000 0.00100000 0.00200000 0.00300000 0.00400000 0.00500000 1 5 9 13 17 21 25 29 33 37 41 45 49 53 57 61 65 69 73 77 81 85 89 Iteration time Strain energy [J] 34 Iteration 4, strain energy= 3.421 mJ Uniform strain energy density distribution Relation of Strain Energy and Iteration Times 0.00000000 0.00100000 0.00200000 0.00300000 0.00400000 0.00500000 1 7 13 19 25 31 37 43 49 55 61 67 73 79 85 91 97 103 Iteration time Strain energy [J] Relation of Strain Energy and Iteration Times 0.00000000 0.00100000 0.00200000 0.00300000 0.00400000 0.00500000 1 9 17 25 33 41 49 57 65 73 81 89 97 105 113 121 129 Iteration time Strain energy [J] 35 Iteration 5, strain energy= 3.407 mJ Uniform strain energy density distribution Iteration 6, strain energy= 3.393 mJ Uniform strain energy density distribution Relation of Strain Energy and Iteration Times 0.00000000 0.00100000 0.00200000 0.00300000 0.00400000 0.00500000 1 11 21 31 41 51 61 71 81 91 101 111 121 131 141 151 161 171 Iteration time Strain energy [J] 36 Iteration 7, strain energy= 3.390 mJ Uniform strain energy density distribution Relation of Strain Energy and Iteration Times 0.00000000 0.00100000 0.00200000 0.00300000 0.00400000 0.00500000 0.00600000 1 13 25 37 49 61 73 85 97 109 121 133 145 157 169 181 193 Iteration time Strain energy [J] Relation of Strain Energy and Iteration Times 0.00000000 0.00100000 0.00200000 0.00300000 0.00400000 0.00500000 0.00600000 1 18 35 52 69 86 103 120137 154 171188 205 222 239 256 273 290 Iteration time Strain energy [J] 37 Iteration 8, strain energy= 3.387 mJ Uniform strain energy density distribution Iteration 9, strain energy= 3.384 mJ Uniform strain energy density distribution Relation of Strain Energy and Iteration Times 0.00000000 0.00100000 0.00200000 0.00300000 0.00400000 0.00500000 0.00600000 1 18 35 52 69 86 103 120 137 154 171 188 205 222 239 256 273 290 Iteration time Strain energy [J] 38 Iteration 10, strain energy= 3.385 mJ Uniform strain energy density distribution Relation of Strain Energy and Iteration Times 0.000000 0.001000 0.002000 0.003000 0.004000 0.005000 0.006000 1 24 47 70 93 116 139 162 185 208 231 254 277 300 323 346 369 392 Iteration time Strain energy [J] Relation of Strain Energy and Iteration Times 0.000000 0.001000 0.002000 0.003000 0.004000 0.005000 0.006000 1 25 49 73 97 121 145 169 193 217 241 265289 313 337 361 385 409 Iteration time Strain energy [J] 39 Iteration 11, strain energy= 3.383 mJ Uniform strain energy density distribution Iteration 12, strain energy= 3.381 mJ Uniform strain energy density distribution Relation of Strain Energy and Iteration Times 0.00000000 0.00100000 0.00200000 0.00300000 0.00400000 0.00500000 0.00600000 1 25 49 73 97 121 145 169 193 217 241 265289 313337 361 385 409 Iteration time Strain energy [J] 40 Iteration 13, strain energy= 3.380 mJ Uniform strain energy density distribution Relation of Strain Energy and Iteration Times 0.00000000 0.00100000 0.00200000 0.00300000 0.00400000 0.00500000 0.00600000 0.00700000 1 31 61 91 121 151 181 211 241 271 301 331 361 391 421 451 481 511 Iteration time Strain energy [J] Relation of Strain Energy and Iteration Times 0.00000000 0.00100000 0.00200000 0.00300000 0.00400000 0.00500000 0.00600000 0.00700000 1 32 63 94 125 156 187 218249280 311 342373 404435 466497 528 Iteration time Strain energy [J] 41 Iteration 14, strain energy= 3.379 mJ Uniform strain energy density distribution Figure 2.8: The topology growth of a Michell cantilever beam structure in 15 iterations. The growth criteria to add one more PSL into the current beam structure is based on a heuristic way. The goal of PSL strategy to mimic the PSL is to approximate the PSL as much as possible. Therefore in each iteration we choose the right position to insert a new PLS to decrease Relation of Strain Energy and Iteration Times 0.00000000 0.00100000 0.00200000 0.00300000 0.00400000 0.00500000 0.00600000 0.00700000 1 37 73 109145 181 217 253289325 361397433 469505 541577 613 Iteration time Strain energy [J] Relation of Strain Energy and Growth Times 0.0000000000 0.0010000000 0.0020000000 0.0030000000 0.0040000000 0.0050000000 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 Growth time Strain energy 0.0001[J] 42 the approximate errors of those profiles of PSL. As shown in Figure 2. 8, based on the initial structure in iteration 0, the maximum errors from the principal stress lines and the related beams are 1a P and 1b P . Accordingly two new principal stress lines will be added which lead to a new node 1c P . Therefore, a new beam structure can be constructed in iteration 1. After size optimization, the strain energy of the new beam structure is 3.80 mJ, which is significantly reduced from that of the beam structure in the previous iteration (4.15 mJ). Hence the stiffness of the beam structure has been increased for the same amount of material. The process can be repeated. From the above topology growth process it can be noticed that as more principal stress lines are used, the beam structure is closer to the Michell structure. The strain energy of the entire structure is getting smaller. However, the performance improvements are getting smaller after several iterations. In our beam structure design, we use the PSL as a bridge between the continuum design domain and the discrete beam structure design. This is different from other methods such as SIMP, Ground truss Structure, and other expansion and growth methods for truss structure optimization (Kirsch, 1997; Martinez 2009). Compared to them, our method is dramatically more intuitive and easier to obtain more Michell type like structure. More importantly, based on our method, the user can easily terminate or add more beams in the designed beam structure. However, a limitation of our method is it can only be used for the minimum compliance or maximum stiffness structural optimization problem. 2.7 Test Examples To demonstrate the principal stress line method presented in this paper, two more examples, in addition to the Michell cantilever beam as shown in Section 2.6, were carried out: (1) a bridge structure, and (2) a bicycle frame. Similar to the Michell cantilever beam, the bridge structure is 43 another classical topology optimization test case. We also present a bicycle frame example, which is more complex with more than one load and constraints. 2.7.1 A Bridge Structure The design domain of a bridge structure is shown in Figure 2.9 on the left. The domain is fixed at the two end points of the bottom and a concentrated force is added in the middle of the bottom edge. The principal stress lines of the design domain based on the method discussed in section 2.3 are shown in Figure 9 on the right. L=H H F L=H Figure 2.9: a bridge like design domain and the principal stress lines. Based on the computed principal stress lines, the skeleton principal stress lines to connect the constraints and loads are shown in Figure 2.9 on the right. Notice the principal stress lines at two constraints have 45 degrees tangent with the horizontal axis. Accordingly, an initial structure can be generated which is also shown in the figure. The topology growth process of the structure is shown in Figure 2.10. After three iterations, the strain energy of the beam structure decreases from 1.33 mJ to 1.11 mJ. Hence the stiffness of the beam structure has been increased for the same amount of material used in the structure. Skeleton Principal Stress Lines Initial Beam Structure 44 Figure 2.10: The topology growth of a bridge structure in 3 iterations. 2.7.2 A Bicycle Frame Structure with Multiple Loads The design domain of a bicycle frame structure is shown in Figure 2.10 on the left. The domain is fixed at two points of the bottom. Multiple forces in different directions are applied on the top of the design domain. Accordingly the principal stress lines of the design domain are shown in Figure 2.11 on the right. 1 2 1 1.175 0.15 450N 350N length:[m] 0.125 350N Strain Energy = 1.327 mJ Iteration 0 Iteration 1 Iteration 2 Strain Energy = 1.135 mJ Strain Energy = 1.109 mJ Topology and shape of beam structures Beam structures after size optimization P 1a P 1b P 2a P 2b Skeleton Principal Stress Lines Initial Beam Structure 45 Figure 2.11. the design domain of a 2D bicycle frame Based on the computed principal stress lines, the skeleton principal stress lines to connect the constraints and loads are shown in Figure 2.11 on the right. Accordingly, an initial structure can be generated which is also shown in the figure. The topology growth process of the structure is shown in Figure 2.12. After five iterations, the strain energy of the beam structure decreases from 3.26 mJ to 3.12 mJ. Hence the stiffness of the beam structure has been increased for the same amount of material used in the structure. Figure 2.12: The topology growth of a bicycle frame structure in 5 iterations. Iteration 0 Iteration 1 Iteration 2 Strain Energy = 3.261 mJ Strain Energy = 3.194 mJ Strain Energy = 3.192 mJ Iteration 3 Iteration 4 Strain Energy = 3.147 mJ Strain Energy = 3.118 mJ Topology and shape of beam structures Beam structures after size optimization P 1a P 2a P 2b P 2c P 2d P 3a P 3b P 4a P 4b 46 Based on the strain energy of the whole beam structure, the performance improvements of the structures during the iterations are shown in Figure 2.13. 0.000 2.000 4.000 6.000 8.000 10.000 12.000 14.000 16.000 18.000 20.000 1 2 3 4 5 6 7 8 9 10 11 12 13 14 improvement(%) Iteration Stiffness Improvement vs. Iteration Cantilever Bridge Bicycle Figure 2.13: The structural performance increase of the three test examples. 2.7.3. Comparison We compared the cantilever optimized results with results from SIMP and ground truss structure approach. From SIMP From PSL strategy Figure 2.14: cantilever topology optimization results comparison between SIMP and PSL. From above figure we can see that the results from SIMP and PSL are different in the resolution of the beam structure meshes. The possible reason is like this. SIMP is a continuum based optimization strategy meshed by rectangle elements. The sensitivity analysis leads to the 47 material to be concentrated along the principal direction. However, the shape and connections of finite elements limits the material distribution along the real principal direction. Therefore it is hard for SIMP to get more refined optimized structure. While the topology of elements in PSL strategy is changed in next growth iteration. Therefore it is possible for PSL strategy to obtain more complex and refined Michel structure like result. In addition we interpreted the SIMP cantilever optimal result into a beam structure and optimized its size on the uniform strain energy density criteria shown in Figure 2.14. Topology sketch Strain Energy = 3.390mJ Figure 2.15: cantilever topology optimization from SIMP is interpreted into a beam structure. Table 2.1: strain energy of growth cantilever at each iteration. Growth # 1 2 3 4 5 6 7 8 9 Strain Energy (mJ) 4.149 3.799 3.471 3.445 3.421 3.407 3.390 3.387 3.384 Growth # 10 11 12 13 14 Strain Energy (mJ) 3.385 3.383 3.381 3.380 3.379 From Table 2.1, we can see the result from SIMP is at iteration 7 in PSL strategy. Then PSL strategy could achieve the more refined, more elegant and more optimal result in terms of strain energy performance. 48 Further we also computed the optimal truss structure from ground truss structure strategy. Here the strain energy problem formulation is utilized and solved in linear mathematical programming. Figure 2.16: cantilever topology optimization result from Ground Truss Structure. From Figure 2.16, it is noticed that there are some trusses are overlapped onto each other. Then it seems that the ground truss structure does not obtain the reasonable result for this cantilever case. 2.8 Extensions 2.8.1 Minimum Compliance of 3D Case Not only can the PSL strategy be applied in 2D domain, but also it can be applied in 3D domain. Here we have an example which shows how PSL is applied in 3D case successfully. The design domain with constraints and external forces of 3D cantilever is shown in Figure 2.17. 49 W W/2 0.866W 1.5W W/2 F Figure 2.17: 3D Cantilever Design domain And the PSL distribution is as shown in Figure 2.18. Figure 2.18: 3D Cantilever PSL distribution plot The PLS growth strategy was applied and we obtained the following results: 50 Figure 2.19: The growth process of 3D cantilever structure design by PSL strategy Table 2.2: strain energy of growth 3D cantilever at each iteration. Growth # 1 2 Strain Energy (mJ) 0.641 0.522 Iteration 0 0.641 mJ Iteration 1 0.522 mJ 51 Figure 2.20: The topology growth of a 3D Michell cantilever beam structure in 2 iterations. We can see strain energy in iteration 1 is smaller than that in iteration 0. 2.8.2 Minimum mass problem Still we use cantilever to show how PSL strategy can be applied in this minimum mass problem structural optimization. We continue to use the final optimized result from section 2.6. To solve the minimum mass problems, several types of constraints are configured. (I) Displacement Constraint From Property 2.1.3.1 and Proposition 2.5.1.1 we can see that volume of a beam structure is proportional to the displacement of nodes. Initial optimized topology with constraints Optimized Result Figure 2.21: The topology optimization results of 3D cantilever with given displacement constraint. Relation of Strain Energy and Growth Times 0.0000000000 0.0001000000 0.0002000000 0.0003000000 0.0004000000 0.0005000000 0.0006000000 0.0007000000 1 2 Growth time Strain energy[J] 52 The displacement at y direction at Point P is:Py = -6.758854e-6 (m). Suppose the displacement at P has upper bound like: |Py| < 5.0e-6(m).Then the volume of optimal structure can be optimized by scaling: ] [ 025 02 . 0 ] [ 015 . 0 0 . 5 6 7 . 6 3 3 0 0 m m V y y V upper     . The cross section area of all beams can be scaled by a coefficient 1.152 0  upper y y . (II) Stress Constraint From Property 2.1.3.1 and Proposition 2.5.1.1 we can see that volume of a beam structure is proportional to the admissible stress of a beam structure. Suppose max  is the maximal Von Mises stress. And the admissible stress range is: admission    max , Then the cross section area of all beams can be scaled by a coefficient admission   max From above minimum mass problem solving cases from different constraints, we can see how the PSL strategy is applied in this more practical engineering problem. First we formulate the minimum compliance problem and use PSL strategy to solve it and obtain an optimized beam structure solution. Second based on the constraints in minimum mass problem the optimized beam structure is scaled to satisfy the lower bound of those constraints. The PSL strategy can solve the minimum mass problem. 2.9 More Test Cases 2.9.1 L-Shaped Domain The design domain is L-Shaped as following figures shows. 53 H H 0.4H 0.4H 0.2H F (a) (b) Iteration 0 Strain Energy = 58.513 mJ Iteration 1 Strain Energy = 45.153 mJ 54 Iteration 2 Strain Energy = 38.681 mJ Iteration 3 Strain Energy = 36.142 mJ Iteration 4 Strain Energy = 35.633 mJ Iteration 5 Strain Energy = 35.419 mJ (c) Figure 2.22: L-Shaped domain (a), its PSL distribution plot (b) and The growth process of LShaped structure 55 design by PSL strategy (c) Table 2.3: strain energy of growth L-Shaped domain at each iteration. Growth # 1 2 3 4 5 6 Strain Energy (mJ) 58.513 45.153 38.681 36.142 35.633 35.419 Figure 2.23: The topology growth of an L-Shaped domain structure in 6 iterations. 2.9.2 Multiple loads onto a Bridge Structure The design domain is as following figure shows. F/3 F/3 F/3 H 3H/10 H/2 H/4 H/4 (a) 0 10 20 30 40 50 60 70 1 2 3 4 5 6 Strain energy (mJ) Iteration Relation of Strain Energy and Iteration Times 56 (b) PSL distribution plot Iteration 0, Strain Energy = 4.412 mJ Iteration 1, Strain Energy = 1.602 mJ Iteration 2, Strain Energy = 1.587 mJ Iteration 3, Strain Energy = 1.522 mJ 57 Iteration 4, Strain Energy = 1.491 mJ (c) Figure 2.24: domain (a), its PSL distribution plot (b) and The growth process of Multiple load bridge structure design by PSL strategy (c) Table 2.4: strain energy of growth multi-loads bridge domain at each iteration. Iteration # 1 2 3 4 5 Strain Energy (mJ) 4.412 1.602 1.587 1.522 1.491 Figure 2.25: The topology growth of multi-loads bridge domain in 5 iterations. 2.9.3 Single load onto a convex polygon domain Relation of Strain Energy and Growth Times 0.0000 0.0010 0.0020 0.0030 0.0040 0.0050 1 2 3 4 5 Growth time Strain energy[J] 58 H H 0.4H 0.4H 0.2H F (a) (b) Iteration 0 Strain Energy = 25.303 mJ Iteration 1 Strain Energy = 13.201 mJ 59 Iteration 2 Strain Energy = 13.078 mJ Iteration 3 Strain Energy = 11.493 mJ Iteration 4 Strain Energy = 10.485 mJ Iteration 5 Strain Energy = 10.191 mJ (c) Figure 2.26: convex L-Shaped domain (a), its PSL distribution plot (b) and The growth process of 60 convex L-Shaped domain design by PSL strategy (c) The growth process of convex LShaped domain design by PSL strategy. Table 2.5: strain energy of growth convex L-Shape domain at each iteration. Growth # 1 2 3 4 5 6 7 Strain Energy (mJ) 25.303 13.201 13.078 11.493 10.485 10.191 25.303 Figure 2.27: The topology growth of convex L-Shape domain in 7 iterations. 2.10 Summary Based on the energy principles and the uniform strain energy density proposition in truss discrete structures, we extended and proved the axiom of strain energy density to the beam structure. Then we utilize the axiom of strain energy density to conduct size optimization of a beam structure. In addition we presented a numerical method to compute the principal stress lines of the stress field in a continuum design domain. Further we proposed a reduction-expansion design method to approximate the principal stress lines in a continuum design domain to a discrete beam structure. Three examples are tested by using the reduction-expansion design method. The results showed the elegance of our design method. Compared to the ground truss structure approach, the reduction-expansion method based on principal stress lines are more predictable and can overcome the numerical computation instability and results of uncertainties when the ground like structure is highly dense. In addition our PSL strategy can work both on Relation of Strain Energy and Iteration Times 0.0000 0.0050 0.0100 0.0150 0.0200 0.0250 0.0300 1 2 3 4 5 6 7 Iteration times Strain energy [J] 61 single load case and multiple load cases by loading all loads at one time. And PSL strategy could be extended to 3D optimization from 2D optimization problem. Furthermore PSL strategy can be extended to minimum mass problem but not limited in this problem from minimum compliance problem. 2.11 Appendix: Proof of Proposition 2.5.1.1 Proposition 2.5.1.1 Let i x is the volume of i-th beam in a planar beam structure, then the element stiffness matrix of the i-th beam can be written as 2 2 1 i i i i d i x x K K K   ; both 1 i K and 2 i K are positive semi-definite matrices. Under the classic Euler-Bernuolli assumption, the stiffness of a beam is dominated by 1 i K in the axial direction and we have i i d i x K K  . Proof: According to Pilkey (Pilkey, 2002), the stiffness matrix for Axial Extension and Bending of a beam in a planar beam structure is found in (2.19):                                        l EI l EI l EI l EI l EI l EI l EI l EI l EA l EA l EI l EI l EI l EI l EI l EI l EI l EI l EA l EA e 4 6 0 2 6 0 6 12 0 6 12 0 0 0 0 0 2 6 0 4 6 0 6 12 0 6 12 0 0 0 0 0 2 2 2 3 2 3 2 2 2 3 2 3 K (2.19) In our study the square cross section of a beam is utilized(the cross section shape does not influence the conclusion here.). Then we have area moment of inertia 12 / 12 2 4 A a I I z y    ( the circular cross section has the similar relation), a is the side length. Then (2.19) can be written in form of: e i i e i i e i x x 2 2 1 K K K   , where i i i l A x  . (2.20) 62                         0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 1 2 1 l E e i K , (2.21)                             2 2 2 2 5 2 4 6 0 2 6 0 6 12 0 6 12 0 0 0 0 0 0 0 2 6 0 4 6 0 6 12 0 6 12 0 0 0 0 0 0 0 12 l l l l l l l l l l l l l E e i K . (2.22) Actually e i i x 1 K is the element stiffness matrix for a truss bar by only considering the axial loading forces; e i i x 2 2 K is the element stiffness matrix for a Bernoulli-Euler Beam when ignores the axial loading forces. From Example 2.10 in (Pilkey, 2002) the element stiffness matrix e i i x 2 2 K can be scaled and written as:                             4 6 0 2 6 0 6 12 0 6 12 0 0 0 0 0 0 0 2 6 0 4 6 0 6 12 0 6 12 0 0 0 0 0 0 0 3 2 2 l EI x e i i K (2.23) It is not difficult to verify that the eigenvalues of matrix (2.21) and (2.23) are non negative. Then both element stiffness matrix e i i x 1 K and e i i x 2 2 K are positive semi-definite. By applying the transform from local stiffness matrix to the global stiffness matrix: e e eT e G T K T K  , (2.24) where e T is the transform matrix from the local to the global coordinate system; and 63 e e                        1 0 0 0 0 0 0 cos sin 0 0 0 0 sin cos 0 0 0 0 0 0 1 0 0 0 0 0 0 cos sin 0 0 0 0 sin cos         T (2.25) Then apply the transform matrix e T onto the equation (2.29), we have: 2 2 1 2 2 1 i i i i e e i eT i e e i eT i d i x x x x K K T K T T K T K     . (2.26) Because e T is an orthogonal transform matrix; then 1 i i x K has the same eigenvalues with matrix e i i x 1 K ; and matrix 2 2 i i x K has the same eigenvalues with matrix e i i x 2 2 K . While all their eigenvalues are non negative, then both global element stiffness matrix 1 i i x K and 2 2 i i x K are positive semi-definite matrix. Also we can see                                                      4 6 0 2 6 0 6 12 0 6 12 0 0 0 0 0 0 0 2 6 0 4 6 0 6 12 0 6 12 0 0 0 0 0 0 0 ) ( 12 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 1 ) ( 3 l a Ea l a Ea e K (2.27) Under the classic Euler-Bernuolli assumption we have l a  ; let 0 ) ( 2   l a  . Now we have e i i e x l a Ea K K                                    3 2 0 6 2 0 2 0 2 0 0 0 1 0 0 1 6 2 0 3 2 0 2 0 2 0 0 0 1 0 0 1 ) (                 (2.28) 64 Where                                   3 2 0 6 2 0 2 0 2 0 0 0 1 0 0 1 6 2 0 3 2 0 2 0 2 0 0 0 1 0 0 1 2                 l E e i K (2.29) Since  is very small, and its variation from the cross section variable change is also very small and has very little influence on e i K in (2.29), we can assume e i K a constant matrix. Hence e K in (2.28) is a positive semi-definite matrix. Following above global transformation procedure, we have the global stiffness matrix of the i-th beam i i d i x K K  (2.30) where e e i eT i T K T K  . We also know i K is a positive semi-definite matrix. □ 65 Chapter 3 Design Deformable Geometries with Target Stiffness 3.1 Digital Material Cell and Digital Material We assume the deformation is linear elastic deformation to evaluate the Young’s modulus of digital materials. 3.1.1 2D Digital Material Cell Figure 3.1: A 33  2D digital material cell. A 2D digital material cell is a N x N square element. Each square element is similar to a pixel in an image. As a black and white image has either black or white pixels, a N x N square element 2D digital material cell has either hard or soft material. Our idea is to achieve different elastic material property of a cell by assigning either soft or hard material on a square element of a 2d digital material cell. To simplify the design and manufacturing N is set to be 3 in our study. Figure 3.1 shows a 2D digital material cell in our study. 3.1.2 The Pattern of a 2D Digital Material Cell A pattern of a 2D digital material cell is an ID or a finger print of a 2D digital material cell. The pattern is utilized to identify a 2D digital material cell uniquely. Since a square element is either soft material or hard material, a binary digit number 0 or 1 can be mapped to a square element. Totally a nine binary digit number is obtained to represent a pattern of a 2D digital material cell. Figure 3.2 shows an ID for each square element in a 2D digital cell. And Figure 3.3 66 shows the mapped 9 digit number w/r the corresponding square element ID in a 2D digital material cell. 7 8 9 4 5 6 1 2 3 Figure 3.2: Sequence ID of square elements in a 2D digital material cell. 1 2 3 4 5 6 7 8 9 Figure 3.3: Sequence of square elements in a mapped nine digit binary number. It can be seen thast an element is mapped to a position in a nine digit binary number. If the bit is 0, then the mapped element is soft material. Otherwise, the mapped element is hard material. For instance, a nine digit number of 111101111 maps the following pattern of 2d digital cell. Figure 3.4 shows how material is assigned into the 2D digital cell with an example pattern 111101111. Totally there are 2 9 different numbers. Then we obtain 512 possible patterns. Hard material Soft material Figure 3.4: 2D digital cell 111101111. 3.1.3 2D Digital Material Construction 2D digital material is constructed from 2D digital cells. If 2D digital material is looked as a 2d image, then the 2D digital cell is like a pixel. Therefore 2D digital material can be defined to be a rectangle array of 2D digital cells. 67 3.1.4 Elastic Modulus Evaluation of a 2D Digital Material For a bar, Young’s modulus E can be computed in the following equation FL E LA     (3.1) FEM simulation is employed to evaluate the elastic modulus of a 2D digital material. C2 C3 C4 C5 L6=6H H Fe = F / H (N/m) C6 C1 L5=5H L4=4H L3=3H L2=2H L1=H 1 l  2 l  3 l  4 l  5 l  6 l  H = N x h N 2 1 ... h Figure 3.5: 2D digital material bar FEM simulation setup with dimensions. Figure 3.5 is a 2D digital material bar domain with fixed left end and uniform loads on the right end. The bar only shows one row of cell here. In our study the rows number is monitored. In the results section we will see how rows number impact evaluation of Young’s modulus of digital material. 68 C2 C3 C4 C5 H Fe = F / H (N/m) C6 C1 L4=4H L3=3H 1 l 2 l 3 l 4 l 5 l 6 l C3 3 , 4 l 3 ,3 l 3 , 2 l 3 ,1 l 3:1 4 , 4 ,1 l 4 ,3 ,1 l 4 , 2 ,1 l 4 ,1,1 l 3:1 4 ,1, 2 l 4 ,1,3 l 4 ,1, 4 l 4 , 2 , 2 l 4 , 2 ,3 l 4 , 2 , 4 l 4 ,3 , 2 l 4 ,3 ,3 l 4 ,3 , 4 l 4 , 4 , 2 l 4 , 4 ,3 l 4 , 4 , 4 l 4 ,1,1 ( , ) hl P1 4 ,1, 2 2 2 ( , ) 3 h l P 4 ,1,3 3( , ) 3 h l P 4 ,1, 4 4 (0 , ) l P Figure 3.6: Target digital cell detail view. The displacements at nodes of middle cell elements are retrieved after the FEA simulation analysis running is done. Then the average Young’s modulus E can be retrieved by evaluating E at each node. All nodes are from the 4-th cell in Figure 3.5. And Figure 3.6 shows the detail view of these displacements of nodes at the 4-th cell. From equation (3.1) and average all these evaluated Young’s modulus of target nodes. It ends up with equation (3.2) 69   4 4 4 3 4, , 1 1 1 1 ( 1) 1 3 4 16 NN i ij i i j F i E E l l N h t             (3. 2) 3.1.5 3D Digital Material Cell Figure 3.7: A 333  3D digital material cell. A 3D digital material cell is a N x N x N cube element. Each cube element is similar to a voxel. A N x N x N voxel element of 3D digital material cell has either hard or soft material. Our idea is to achieve different elastic material property of a cell by assigning either soft or hard material on a voxel element of a 3d digital material cell. To simplify the design and manufacturing N is set to be 3 in our study. Figure 3.7 shows a 3D digital material cell used in our study. 3.1.6 The Pattern of a 3D Digital Material Cell A pattern of a 3D digital material cell is an ID or a finger print of a 2D digital material cell. This pattern is utilized to identify a 3D digital material cell uniquely. Since a cube element is either soft material or hard material, a binary digit number 0 or 1 can be mapped to a cube element. Totally a twenty-seven binary digit number is obtained to represent a pattern of a 3D digital material cell. 70 1 2 4 5 3 6 7 8 9 10 11 13 14 12 15 16 17 18 19 20 22 23 21 24 25 26 27 Figure 3.8: Sequence ID of cube elements in a 3D digital material cell. Figure 3.8 shows all sequence ID of cube elements in a 3D digital material cell. Since 27 bit number has 2 27 combinations, we do not plan to enumerate all possibilities. Therefore only Isotropic and transversely isotropic materials are investigated in our study. For an isotropic 3D digital material cell, it requires elements are distributed symmetrically in all x, y and z directions. Then 27 elements in a 3D digital cell can be classified into following four groups: (1,3,7,9,19,21,25,27), (2,8,20,26,4,6,22,24,10,12,18,16), (11,15,17,13,5,23), (14). Each group shares the common material, either soft or hard. Then a 4-digit number can be utilized to represent an isotropic 3D digital material cell as shown in Figure 3.9. 1 2 3 4 Figure 3.9: Sequence ID of cube element groups in a mapped 4-digit binary number. It can be seen that a group of elements are mapped to a position in a 4-digit binary number. If the bit is 0, then the mapped group of element is soft material. Otherwise, the mapped group of elements is hard material. Totally there are 2 4 different numbers. Then we obtain 16 patterns of 71 isotropic 3D digital material. Figure 3.10 shows the more intuitive colored cube elements groups as well as a mapped 4-digit binary number. Figure 3.10: Colored cube element groups as well as a mapped 4-digit binary number. For a transversely isotropic 3D digital material cell, it requires elements are distributed symmetrically in all xy, yz or zx planes. Then 27 elements in a 3D digital cell can be classified into following three groups as well as rotation: No rotation: (1,3,7,9,10,12,16,18,19,21,25,27), (2,4,6,8,11,13,15,17,20,22,24,26), (5,14,23) Rotation: (1,7,25,19,2,8,26,20,3,9,27,21), (4,16,22,10,5,17,23,11,6,18,24,12), (13,14,15) Each group shares the common material, either soft or hard. Then a 4-digit number can be utilized to represent a transversely-isotropic 3D digital material cell. The fourth digit represents whether this cell is rotated or not. Totally there are 2 4 different numbers. Then we obtain 16 patterns of transversely isotropic 3D digital material. 72 3.1.7 3D Digital Material Construction 3D digital material is constructed from 3D digital cells. If 3D digital material is looked as a 3d image, then the 3D digital cell is a 3d voxel like a 2d pixel. Therefore 2D digital material can be defined to be a 3-dimensional array of 3D digital cells. 3.1.8 Elastic Modulus Evaluation of 3D Digital Material The evaluation FEA simulation setup is very similar to that of 2D digital material. Let us write down the evaluation formula (3.3).   4 4 4 4 3 4, , , 1 1 1 1 1 ( 1) 1 3 4 16 N N N i i j k i i j k F i E E l l N h t               (3.3) 3.1.9 Customize a Digital Material Cell with Multi-levels Above digital material cell pattern only works for hard material and soft material. And there are some disadvantages or limitations: i). limited number of symmetric pattern is isotropic material. No more varied materials are available. ii). the resolution from soft material to hard material is very low. In order to address above issues or achieve high resolution of varied material between soft material and hard material, we presented a multi-level interpretation strategy. The idea is as follows: i). At level 0 given two materials, soft material with Young’s modulus Ea and hard material with Young’s modulus Eb; Ea < Eb. Obviously two basis patterns are: (0), (1) bb with ( (0)) ( (1)) E b E b  (3.4) ii). At level 1 eight symmetric or isotropic digital materials are available discussed in section 3.1.2. A series of 2D or 3D digital cell patterns are chosen to be basis patterns: (0,0), (0,1),..., (0, 1) b b b N  with ( (0,0)) ( (0,1)) ... ( (0, 1)) E b E b E b N     (3.5) 73 iii). At level 2, material pair of (soft material, hard material) is chosen from two next basis patterns at level 1,   (0, ), (0, 1) ,0 1 b i b i i N     . Then the pattern strategy is applied from section 3.1.2 to obtain another eight isotropic material patterns. Also another series of 2D or 3D digital cell patterns are chosen to be basis patterns: (0, ,0), (0, ,1),..., (0, , 1) b i b i b i N  with 0 1, ( (0, ,0)) ( (0, ,1)) ... ( (0, , 1)) i N E b i E b i E b i N        (3.6) iv). The interpretation process can be executed in a iterative way to deeper and deeper levels. The interpretation resolution to obtain a specific elastic modulus of material can be unlimited small close to zero. However, for the finite element simulation numerical computation error the levels are limited. 3.1.10 Monotone Increasing Rules on Basis Selection In order to make the basis interpretation idea work for multi-levels some rules have to obey: i). Patterns are sequenced with hard material ratio increasing in monotone way. ii). Patterns are sequenced with optimal structure increasing in monotone way. These two rules lead to another elegant property: the basis patterns are exactly the same sequenced ones from level one to the following deeper levels. 3.2 Procedure of Digital Material Design 3.2.1 Tools and Setup For 2D digital material, COMSOL(www.comsol.com) is employed to configure the finite element domain, loads and constraints conditions. And for 3D digital material, VEGA (http://run.usc.edu/vega/) is employed to conduct the similar setup. The overall simulation setup is discussed in section 3.1.4 either for 2D or 3D. 3.2.2 Verification 74 To verify whether the entire setup is correct or not, the entire domain is configured to be either soft or hard material. Then its Young’s modulus is evaluated to see if it is close to the Young’s modulus of the soft or the hard material. To verify if each element in a digital cell is configured to be soft or hard in a right way, we double check them by humane being in COMSOL. Since we have those patterns of 2D digital cells in our design sketches, so we randomly pick some patterns generated by our program and import them into COMSOL to visualize to see if they are configured correct or not. For 3d digital cell patterns, we utilized two ways to verify if the setup is configured correctly or not. One way is to output all elements property-soft or hard into a output text file for level one. Since there are 27 elements in a 3D digital cell at level one, we can employ this manually double check method for both isotropic and transversely isotropic 3d digital material. Another way is use test functions to configure each element to be soft or hard for level two and three of 3D digital material. The results from the test function are compared to results from the program. Once they are the same, we can conclude the verification is correct safely. 3.3 Steps to Evaluate 2D Digital Material Step 1: Tensile Experiment FEA Simulation Setup 6m 1m F = 1KN/m Digital cell Figure 3.11: 2D digital material bar FEM simulation setup. 75 Figure 3.11 shows the FEA simulation setup. The domain of this setup is composed of 5 by 30 digital cells. Step 2: Select Material or a Material pattern as a cell At level zero the digital cell is a 1 by 1 square element as shown in Figure 3.12. And we have two types of material, soft material S and hard material H. The Young’s modulus of soft material is 2.68 MPa. And the Young’s modulus of hard material is 2680MPa. Both of their Poisson ratio are 0.33.This cell at level zero can be filled with either material S or H. Hard material H Soft material S Figure 3.12: 2D digital cell at level zero. At level one the digital cell consists of 3 by 3 square elements. Each element can be filled with material S or H. Then we can get one pattern of cell at level one as shown in Figure 3.13. Hard material H Soft material S Figure 3.13: 2D digital cell at level one. At level two the digital cell consists of 9 by 9 square elements. Still we have two types of soft material S and H. But here S is one pattern of cell from level one. And H is another pattern of cell from level one. Recall that the digital cell at level two is composed of 3 by 3 square elements. Since each element is either filled with soft material or hard material from a digital cell at level one with 3 by 3 elements, then the digital cell at level two consists of 9 by 9 square elements as shown in Figure 3.14. 76 Soft material S Hard material H Figure 3.14: 2D digital cell at level two. At level three the digital cell consists of 27 by 27 square elements. Still we have two types of soft material S and H. But here S is one pattern of cell from level two. And H is another pattern of cell from level two. Recall that the digital cell at level three is composed of 3 by 3 square elements. Since each element is either filled with soft material or hard material from a digital cell at level two with 9 by 9 elements, then the digital cell at level three consists of 27 by 27 square elements as shown in Figure 3.15. Soft material S Hard material H Figure 3.15: 2D digital cell at level three. Iteratively, at level k (k>=1) the digital cell consists of 3 k by 3 k square elements. The cell at level k is 3 by 3 elements. Each element is either filled with soft material or hard material with cell pattern of 3 k-1 by 3 k-1 square elements at level k - 1. Step 3: Fill the Selected Cell into the Cell in the Simulation Domain 77 The simulation domain consists of 5 by 30 digital cells. And the selected cell is composed of 3 k by 3 k square elements at level k. Totally the domain contains 21 50 3 k   number of square elements. Step 4: Run FEA and Retrieve the Displacements of Elements The FEA simulation and analysis tool is running and the displacement of elements are retrieved. Step 5: Evaluate Young’s Modulus of 2D Digital Material Use equation in (3.3) then we have:   4 20 20 3 4, , 1 1 1 11 ( 1) 1 3 4 16 i ij i i j i E E l l t            (3.7) 3.4 Steps to Interpolate 2D Digital Material Step 1: Select Basis Material or Basis Material Pattern as a Cell At level zero the digital cell is a 1 by 1 square element. And we have two types of material, soft material b(0) and hard material b(1). The Young’s modulus of soft material is 2.68 MPa. And the Young’s modulus of hard material is 2680MPa. Totally we have two types of material b(0) and b(1), shown in Figure 3.16. Hard material b(1) Soft material b(0) Figure 3.16: 2D digital cell patterns at level zero. At level one the digital cell consists of 3 by 3 square elements. As discussed before, we have five basis isotropic patterns. Let us list all of them in Figure 3.17. b(0,0) b(0,1) b(0,2) b(0,3) b(0,4) Figure 3.17: 2D digital cell patterns at level one. 78 At level two the digital cell consists of 9 by 9 square elements shown in Figure 3.18. The cell pattern at level two is interpolated from two neighbor cell pattern at level one. Basically it can be seen that two 3 by 3 square 2D digital cell patterns are filled into the 3 by 3 digital cell. b(0,0,0) b(0,0,1) b(0,0,2) b(0,0,3) b(0,0,4) b(0,1,0) b(0,1,1) b(0,1,2) b(0,1,3) b(0,1,4) b(0,2,0) b(0,2,4) b(0,2,1) b(0,2,2) b(0,2,3) b(0,3,0) b(0,3,1) b(0,3,2) b(0,3,4) b(0,3,3) Figure 3.18: 2D digital cell patterns at level two. Here we visualize the digital basis cell patterns at level zero, one and two. We omit the visualization of digital basis cell pattern at level three and more since we have no enough space to 79 hold them. The visualization procedure can follow section 3.1.9 (customize digital material in multi-level) and Step 2 in section 3.3 (select material or a material pattern as a cell). Step 2: Evaluate the Young’s Modulus of Basis Patterns Once we basis patterns, we can follow the procedure in section 3.4(steps to evaluate 2D digital material) to evaluate Young’s modulus of digital material with basis pattern. The interpolated material property is evaluated level by level. Then more and more digital material with varied elastic property is retrieved in a monotone way. 3.4.1 Coincident Cell Pattern at Different Levels It’s interesting to notice those span end points at level k – 1 are also interpolated in level k. A question is whether this coincident cell pattern obtain the same young’s modulus evaluation from level k - 1 and level k. Intuitively this coincident cell pattern from k-1 is scaled 1/3 and duplicated by 9 copies and arrayed into the cell pattern at level k. At the first sight, the cell pattern at level k-1 is different from its 9 copies at level k for the dimension and topology in one single cell. However, recall that the young’s modulus evaluation is based on the entire beam bar simulation domain setup instead of one single cell. Therefore when inserting the single cell pattern at level k – 1 back into the setup domain. There would be cluster of cell pattern with the congruent pattern with those at level k. In addition, we assume the deformation is elastic deformation in our study. From the equilibrium equation: KU = F. K is dimensionless. Then the domain with elements in level k – 1 can be scaled into the same element size with that in level k without impacting its stiffness matrix K as well the digital material’s young’s modulus. We conclude the coincident cell pattern material’s young’s modulus at level k – 1 is the same with that at level k for their identical topology and dimensionless K in a clustered digital cells point of view. 80 3.5 Results of Digital Material Design 3.5.1 2D digital material patterns at level one Obviously, the 2D square cell shape has symmetry or mirror one about itself. Therefore for any pattern of nine digit number, it can rotate about itself by 90, 180, 270 degrees to obtain all derived symmetry ones. Also it can mirror itself about x axis or y axis to obtain all derived mirror ones. This original pattern, its derived symmetry ones and its derived mirror ones are grouped to one pattern supposed to have the same elastic material property. Certainly we have to evaluate Young’s modulus in both x and y directions to cover those non symmetric patterns. For instance, Figure 3.19 illustrates such a group of patterns. Figure 3.19: 2D digital cell 000000011 and its equivalents. 81 Once 512 patterns are simplified or grouped by finding the symmetry and mirror ones, 8 symmetric patterns and 94 non symmetric ones are defined. Table 3.1 and Figure 3.20 is the Young’s modulus evaluation from FEA for isotropic 2D digital material. Table 3.1: Evaluated Young’s modulus of all isotropic 2d digital cell patterns at level 1. E (Pa) Pattern Index Figure 3.20: Sorted isotropic 2D digital cell patterns at level one with increment of stifness. It is easy to see the increasing elastic property from soft material to hard material. Also we can see the increment in two next patterns is not uniform. Therefore it’s possible for us to choose the monotone increasing structure optimized and hard material ratio series of basis patterns. In addition Table 3.2 and shows 94 heterogeneous patterns as well as 8 isotropic patterns. 82 Table 3.2: Evaluated Young’s modulus of all anti-isotropic 2d digital cell patterns at level 1. 83 84 85 3.5.2 Impact of Number of Rows on Evaluation of Elastic Modulus In general the more finite elements there are, the more precise the simulation analysis is. We are interested in how rows of cells in the bar impact on evaluation of elastic modulus. Then an isotropic pattern 000000000 in Figure 3.21 is chosen to observe the impact shown in Table 3.3 and Figure 3.22. Figure 3.21: 2D digital cell with pattern 000000000. Table 3.3: Evaluated Young’s modulus of 2D digital cell with pattern 000000000 with varied count of cell rows. Count of cell rows 1 3 5 10 E (MPa) 2.697466 2.694189 2.693730 2.693425 86 Figure 3.22: Graph of evaluated Young’s modulus of 2D digital cell with pattern 000000000 with varied count of cell rows. Obviously it is shown that the curve reaches to a plateau with increment of count of cell rows. 5 is chosen in our study. In addition an anti-isotropic pattern 001011111 in Figure 3.23 is chosen to observe this impact as shown in Table 3.4 and Figure 3.24. Figure 3.23: 2D digital cell with pattern 001011111. Table 3.4: Evaluated Young’s modulus with varied count of cell rows. Count of cell rows 1 3 5 10 E (MPa) 136.882500 909.697300 989.326900 1023.461000 87 Figure 3.24: Graph of evaluated Young’s modulus with varied count of cell rows. The trend reaches plateau with increment of rows count of cells. The total curve trend is different from the isotropic one because of shear deformation caused by non-symmetric property of this anti-isotropic digital cell. 3.5.3 Pattern Basis on Monotone Increment Rules at Level One It is observed that there are eight isotropic patterns in level one. As shown in Table 3.5 and 3.25, a series of basis are to be picked out on the monotone increasing rules mentioned in previous section. Table 3.5: 2D digital material basis patterns and their evaluated Young’s modulus. Binary # 0000000 00 0000100 00 1010101 01 0101110 10 1111111 11 Cell pattern Ratio of Hard 0 1 5 5 9 88 Material(1/9) E( MPa ) 2.69373 0 3.37823 5 341.463 0 966.237 4 2697.46 6 Basis code b(0,0) b(0,1) b(0,2) b(0,3) b(0,4) b(0,0) b(0,1) b(0,2) b(0,3) b(0,4) Figure 3.25: 2D digital material cell basis patterns with monotone increment of stiffness. 3.5.4 Pattern Basis on Monotone Increment Rules at Level Two At level 2 basis cell patterns from level 1 are interpolated into cell patterns at level 2. Figure 3.26 shows the soft and hard material patterns at the first span of level 1 to be interpolated based on basis selected. Table 3.6 and Figure 3.27 show the interpolated results at the first span of level 1. b(0,0) b(0,1) Soft material Hard material Figure 3.26: Soft and hard material patterns from the first span of level 1. Table 3.6: Evaluated Young’s modulus of all interpolated patterns at the first span of level 1. Binary # 0000000 00 0000100 00 1010101 01 0101110 10 1111111 11 89 Cell pattern Filled pattern Ratio of Hard Material(1/81) 0 1 5 5 9 E( MPa ) 2.69364 6 2.75665 5 3.03136 4 3.05006 8 3.38240 8 Basis code b(0,0,0) b(0,0,1) b(0,0,2) b(0,0,3) b(0,0,4) b(0,0,0) b(0,0,1) b(0,0,2) b(0,0,3) b(0,0,4) b(0,0,0) b(0,0,1) b(0,0,2) b(0,0,3) b(0,0,4) Figure 3.27: Left: Graph of evaluated Young’s modulus of all interpolated patterns from the first span of level 1; right: normalized graph. Figure 3.28 shows the soft and hard material patterns at the second span of level 1 to be interpolated based on basis selected. Table 3.7 and Figure 3.29 show the interpolated results at the second span of level 1. Soft material Hard material b(0,1) b(0,2) Figure 3.28: Soft and hard material patterns from the second span of level 1. Table 3.7: Evaluated Young’s modulus of all interpolated patterns from the second span of level 1. Binary # 0000000 00 0000100 00 1010101 01 0101110 10 1111111 11 90 Cell pattern Filled pattern Ratio of Hard Material(1/81) 9 13 29 29 45 E( MPa ) 3.38240 8 3.91137 5 13.4260 1 126.509 3 343.255 0 Basis code b(0,1,0) b(0,1,1) b(0,1,2) b(0,1,3) b(0,1,4) b(0,1,0)b(0,1,1)b(0,1,2)b(0,1,3)b(0,1,4) b(0,1,0) b(0,1,1)b(0,1,2)b(0,1,3)b(0,1,4) Figure 3.29: Left: Graph of evaluated Young’s modulus of all interpolated patterns from the second span of level 1; right: normalized graph. Figure 3.28 shows the soft and hard material patterns at the third span of level 1 to be interpolated based on basis selected. Table 3.8 and Figure 3.29 show the interpolated results at the third span of level 1. Soft material Hard material b(0,2) b(0,3) Figure 3.30: Soft and hard material patterns from the third span of level 1. Table 3.8: Evaluated Young’s modulus of all interpolated patterns at the third span of level 1. Binary # 0000000 00 0000100 00 1010101 01 0101110 10 1111111 11 91 Cell pattern Filled pattern Ratio of Hard Material(1/81) 45 45 45 45 45 E( MPa ) 343.255 0 380.441 7 529.593 7 617.123 4 966.002 1 Basis code b(0,2,0) b(0,2,1) b(0,2,2) b(0,2,3) b(0,2,4) b(0,2,0)b(0,2,1)b(0,2,2)b(0,2,3)b(0,2,4) b(0,2,0) b(0,2,1)b(0,2,2) b(0,2,3) b(0,2,4) Figure 3.31: Left: Graph of evaluated Young’s modulus of all interpolated patterns from the third span of level 1; right: normalized graph. Figure 3.32 shows the soft and hard material patterns at the fourth span of level 1 to be interpolated based on basis selected. Table 3.9 and Figure 3.33 show the interpolated results at the fourth span of level 1. Soft material Hard material b(0,3) b(0,4) Figure 3.32: Soft and hard material patterns from the fourth span of level 1. Table 3.9: Evaluated Young’s modulus of all interpolated patterns at the fourth span of level 1. Binary # 0000000 0000100 1010101 0101110 1111111 92 00 00 01 10 11 Cell pattern Filled pattern Ratio of Hard Material(1/81) 45 49 65 65 100 E( MPa ) 966.0021 1009.004 1510.783 1687.922 2693.646 Basis code b(0,3,0) b(0,3,1) b(0,3,2) b(0,3,3) b(0,3,4) b(0,3,0) b(0,3,1) b(0,3,2) b(0,3,3) b(0,3,4) b(0,3,0) b(0,3,1) b(0,3,2) b(0,3,3)b(0,3,4) Figure 3.33: Left: Graph of evaluated Young’s modulus of all interpolated patterns from the fourth span of level 1; right: normalized graph. Apparently we can see the trend in each interpretation graph is strictly increasing. Elegantly the basis is exactly the same in each interpretation graph. The total graph is obtained when putting all four normalize integrated graphs into one single graph in Figure 3.34. Figure 3.34: Normalized graph of all four span curves at level 2. 0.00E+00 5.00E-01 1.00E+00 1.50E+00 0 2 4 6 93 3.5.5 Pattern Basis on Monotone Increment Rules at Level Three At level 3 basis cell patterns from level 2 are interpolated into cell patterns at level 3. Figure 3.35 shows the soft and hard material patterns at the first span of level 2 to be interpolated based on basis selected. Table 3.10 and Figure 3.37 show the interpolated results at the first span of level 2. And Figure 3.36 shows the detail view of interpolated 2D cell patterns from the first span of level 2 at level 3. b(0,0,0) b(0,0,1) Soft material Hard material Figure 3.35: Soft and hard material patterns from the first span of level 2. Table 3.10: Evaluated Young’s modulus of all interpolated patterns at the first span of level 3. Binary # 0000000 00 0000100 00 1010101 01 0101110 10 1111111 11 Cell pattern Ratio of Hard Material(1/729) 0 1 5 5 9 E( MPa ) 2.694007 2.701425 2.728280 2.72955 2.757052 Basis code b(0,0,0,0) b(0,0,0,1) b(0,0,0,2) b(0,0,0,3) b(0,0,0,4) b(0,0,0,0) b(0,0,0,1) b(0,0,0,2) b(0,0,0,3) b(0,0,0,4) Figure 3.36: Detail view of 2D cell patterns from the first span of level 2. 94 b(0,0,0,0) b(0,0,0,1) b(0,0,0,2) b(0,0,0,3) b(0,0,0,4) Figure 3.37: Graph of evaluated Young’s modulus of interpolated patterns from the first span of level 3. Also we can see the elastic property is increasing in a monotone way. 3.5.6 3D Digital Isotropic Material Patterns at Level One Recall that we employ a four digit number to represent the pattern of an isotropic 3D digital material. Like the 2D digital material FEA simulation setup, a 5 row in y direction by 5 row in z direction and 30 column in x direction array of 3D elements array is configured as a 3D digital material bar domain. Then the elastic modulus of each pattern is evaluated on the simulation constraints, loads and results. Table 3.11 shows all isotropic 3D digital material basis patterns and their evaluated Young’s modulus. Table 3.11: Isotropic 3D digital material basis patterns and their evaluated Young’s modulus. Binary number 0000 (0) 0001 (1) 0010 (2) 0011 (3) 0100 (4) 0101 (5) 0110 (6) 0111 (7) Cell pattern side view Center element 95 Ratio of Hard Material (1/27) 0 8 12 20 6 14 18 26 Ex=Ey= Ez ( MPa ) 2.708617 11.23852 1047.983 1790.737 346.5029 1246.248 1628.751 2590.925 Binary number 1000 (8) 1001 (9) 1010 (10) 1011 (11) 1100 (12) 1101 (13) 1110 (14) 1111 (15) Cell pattern side view Center element Ratio of Hard Material (1/27) 1 9 13 21 7 15 19 27 Ex=Ey= Ez ( MPa ) 3.671979 323.1294 1193.056 1999.265 370.5982 1312.160 1723.744 2708.610 Figure 3.38: Graph of evaluated Young’s modulus of all isotropic 3D cell patterns at level 1. 3.5.7 3D Digital Transversely Isotropic Material Patterns at Level One 96 For the transversely isotropic material patterns we have the following results as shown in Table 3.12. Table 3.12: transversely Isotropic 3D digital material basis patterns and their evaluated Young’s modulus. Binary number 0000 (0) 0001 (1) 0010 (2) 0011 (3) 0100 (4) 0101 (5) 0110 (6) 0111 (7) Cell pattern layer view Pole z z z z z z z z Ratio of Hard Material (1/27) 0 12 12 24 3 15 15 27 Ex( MPa ) 2.708617 8.535228 937.6360 2290.566 6.402255 1029.665 1123.096 2708.610 Binary number 0000 (0) 0001 (1) 0010 (2) 0011 (3) 0100 (4) 0101 (5) 0110 (6) 0111 (7) Cell pattern layer view Pole x x x x x x x x Ratio of Hard Material (1/27) 0 12 12 24 3 15 15 27 Ex( MPa ) 2.708617 1109.441 1096.883 2387.298 266.9560 1392.843 1476.475 2708.610 Transversely isotropic material is symmetric in a plane normal to the pole. 3.5.8 3D Digital Isotropic Material Pattern Basis at Level One The basis patterns at level one is selected out. And the results are shown in Table 3.13 and Figure 3.39. Table 3.13: 3D digital material basis patterns and their evaluated Young’s modulus. Binary number 1000 (0) 1000 (8) 1100 (12) 1101 (13) 1110 (14) 1111 (15) 97 Cell pattern side view Basis code b(0,0) b(0,1) b(0,2) b(0,3) b(0,4) b(0,5) Center element Ratio of Hard Material(1/27 ) 0 8 7 15 19 27 Ex=Ey=Ez ( MPa ) 2.708617 11.23852 370.5982 1312.160 1723.744 2708.610 Figure 3.39: 3D digital material cell basis patterns with monotone increment of stiffness. 3.5.9 3D Digital Isotropic Material Patterns at Level Two We show the results of 3D digital isotropic material of different patterns at level two as shown in Table 3.14~3.18 and Figure 3.40 ~3.44. Table 3.14: 3D isotropic material cell patterns – Level 2 from span 0 of BASIS 1 at Level 1 Binary number 1000 (0) 1000 (8) 1100 (12) 1101 (13) 1110 (14) 1111 (15) Cell pattern side view Basis code b(0,0,0) b(0,0,1) b(0,0,2) b(0,0,3) b(0,0,4) b(0,0,5) 98 Center element Hard Material Percentage(% ) 0 0.14 0.96 2.06 2.61 3.70 Ex=Ey=Ez ( MPa ) 2.71013 2.738148 2.930787 3.177479 3.361755 3.682196 Figure 3.40: Young’s modulus of all interpolated 3D cell patterns at the span 0 of level 2. Table 3.15: 3D isotropic material cell patterns – Level 2 from span 1 of BASIS 1 at Level 1 Binary number 1000 (0) 1000 (8) 1100 (12) 1101 (13) 1110 (14) 1111 (15) Cell pattern side view Basis code b(0,1,0) b(0,1,1) b(0,1,2) b(0,1,3) b(0,1,4) b(0,1,5) Center element Hard Material Percentage(% ) 3.70 6.53 9.47 16.05 19.34 26.93 Ex=Ey=Ez ( MPa ) 3.682196 3.900193 46.047111 48.977996 207.68405 370.32837 99 Figure 3.41: Young’s modulus of all interpolated 3D cell patterns at the span 1 of level 2. Table 3.16: 3D isotropic material cell patterns – Level 2 from span 2 of BASIS 1 at Level 1 Binary number 1000 (0) 1000 (8) 1100 (12) 1101 (13) 1110 (14) 1111 (15) Cell pattern side view Basis code b(0,2,0) b(0,2,1) b(0,2,2) b(0,2,3) b(0,2,4) b(0,2,5) Center element Hard Material Percentage(% ) 26.93 27.02 33.61 42.39 46.78 56.56 Ex=Ey=Ez ( MPa ) 370.32837 381.94620 548.56772 836.79870 996.16561 1316.8655 100 Figure 3.42: Young’s modulus of all interpolated 3D cell patterns at the span 2 of level 2. Table 3.17: 3D isotropic material cell patterns – Level 2 from span 3 of BASIS 1 at Level 1 Binary number 1000 (0) 1000 (8) 1100 (12) 1101 (13) 1110 (14) 1111 (15) Cell pattern side view Basis code b(0,3,0) b(0,3,1) b(0,3,2) b(0,3,3) b(0,3,4) b(0,3,5) Center element Hard Material Percentage(% ) 56.56 56.10 59.40 63.79 66.98 70.37 Ex=Ey=Ez ( MPa ) 1316.8655 1342.2890 1468.1874 1626.4207 1631.8966 1729.5892 101 Figure 3.43: Young’s modulus of all interpolated 3D cell patterns at the span 3 of level 2. Table 3.18: 3D isotropic material cell patterns – Level 2 from span 4 of BASIS 1 at Level 1 Binary number 1000 (0) 1000 (8) 1100 (12) 1101 (13) 1110 (14) 1111 (15) Cell pattern side view Basis code b(0,4,0) b(0,4,1) b(0,4,2) b(0,4,3) b(0,4,4) b(0,4,5) Center element Hard Material Percentage(% ) 70.37 71.47 78.05 86.83 91.22 100.00 Ex=Ey=Ez ( MPa ) 1729.5892 1756.9801 1960.15 2273.45 2408.02 2710.13 102 Figure 3.44: Young’s modulus of all interpolated 3D cell patterns at the span 4 of level 2. 3.5.10 3D Digital Isotropic Material Patterns at Level Three Figure 3.45 shows the monotone interpolation results at level three. Totally there are 126 interpretation points.(for the memory limitation, the simulation is 1 row by 1 row by 6 column). Figure 3.45: Young’s modulus w/r 3D isotropic cell patterns at level 3. There are three outliers at span b(0,3,3). The possible causes to these outliers are: 2.68 502.68 1002.68 1502.68 2002.68 2502.68 3002.68 0 20 40 60 80 100 120 140 Evaluated Young's Modulus(Pa) w/r Cell Patterns at Level 3 E(MPa) 103 1). (1626.422675, 1630.040429) span is very small. 2). 1row by 1row by 6column is very coarse to do FEA simulation. The numerical computation error is relatively big. 3.6 Case Study We conducted a physically tensile experiment to test the Young’s modulus for a 2 by 2 digital cell with four patterns. The material property refers to the experimental ones in simulation: Hard material: 213064121.6 (Pa) Soft material: 778786.4852 (Pa) Table 3.19: Stiffness comparison of 2D digital material simulation result w/r physical experiment result. pattern Hard material ratio (1/4) 4 3 2 1 0 Experiment E(Pa) 213064121.6 70661261 35089468.91 5954090 778786.4852 Simulation E(Pa) 214132762.3 115582661.8 29713667.22 1335410.639 782660.4193 Figure 3.46: stiffness comparison of 2D digital material simulation result w/r physical experiment result. 104 As shown in Table 3.19 and Figure 3.46 the trend of simulation results matches the trend of experimental results very well. 3.7 Target Deformation Design Based on N Level stiffness Digital Material We presented varied stiffness digital material approach to achieve the desired or target flexible/stiffness skin design. These digital materials are printed by 3D printer with combination of two different materials. One is stiff material, another is soft material. Since stiffness of digital materials required varies from the soft to the stiff. We presented an iterative strategy to interpolate the stiffness of digital material from two given materials. This digital material interpolation strategy is to be discussed in the subsequent sections. From the interpolated digital material library some of them are selected to serve the flexible skin design. These selected interpolated digital materials are called N levels. In addition this varied stiffness flexible skin design is based on N levels digital material. 3.7.1 Computational Framework for Desired Flexible Skin Design A large number of elements in a skin domain pose heavy computational challenges in approximating, analyzing, synthesizing, and optimizing a skin design to meet the given flexibility requirements. In this proposal, we only consider elastic deformations without special explanation. 105 Optimization Configuration of Skin Design Domain Description Constraints Loads Design Variables Skin Design Result Meet Requirements? NO YES CAD Generator 3D Printing FEA Figure 3.47: Computational framework for desired flexible skin design. The entire computational framework for desired flexible skin design is shown in Figure 3.47. The entire design flow is to configure skin domain, constraints, external loads and design variables first. Also the target design shape is described (not shown here). The design variables are to be determined by an optimization module. At the beginning design variables are assigned with initial values. Then the entire skin design domain with constraints, external loads and initial design variables are send to FEA analysis module to run, analyze and output the design result. Next the retrieved running skin design result is to compare with the given target design shape or pattern. If they match each other very well or the result is the most optimal solution, then the result is accepted and converted to a CAD generator to generate a STL geometry model. This resulted STL geometry model is sent to an AM machine to be fabricated in the end. Otherwise, the design variables would be modified by an optimization module using some criteria. Then another round of FEA simulation and analysis of the domain with updated design variables and other setting is executed. This design process is an iterative way until a most optimal solution is obtained or all possible cases in the design space are enumerated. 106 3.7.2 Target Deformation Design Methodology Observation: We observed that a continuous stiff edge is next to varied material with large difference of stiffness, this continuous stiff edge would be a straight line segment under uniform loads when the domain is configured properly. Based on this observation, the desired skin curve is simplified to be a polygon by approximation. Furthermore, this approximated polygon decreases the design solution space greatly. We will discuss this solution space simplification in the following paragraphs. Before we utilize above framework to design the desired skin we presented the following steps or methodology to implement our idea. The N level printed digital material vector is an array holding all digital material from soft to stiff. GIVEN: target skin curve T C ; a 2D domain being meshed into D elements; distributed loads; both ends fixed; N level printed digital material vector   12 , , , dm N V E E E  . FIND: Digital Material i E at j th  element SATISFY: 01 ( , , , ) k CP     , 1,2, , iN  , 1,2, , jD  , 1,2, , kM  MINIMIZE: T CC  Where k  is the internal angle at k th  joint; P is a polygon shape function. Step 1: Configure the 2D domain constrained two ends, distributed loads, N level printed digital material vector dm V and target skin shape T C . Step 2: The given target skin curve is approximated by a polygon with N edges. Curve 01 ( , , , ) k CP     is defined. Step 3: This N-edge polygon is projected onto a 2D bar domain with two ends fixed. Then this bar 2D domain is divided into N chained sub-domains. These N chained sub-domains contain 107 some amount of elements separately. These elements are decomposed into two sets. One set consists of elements with respect to edges on the polygon; another set consists of elements are responding to vertices or joints on the polygon. Step 4: Those elements of links are assigned with most stiff digital material. And those elements of joints are set to be design variables. Each design variable varies among N level digital materials. Therefore, assume there are M design variables in this 2D bar domain; then the design solution space is N M . Step 5: Convert design variables and target polygon shape into an optimization problem. Then framework to design desired skin is employed to run and find the solution for this optimization problem. It is noticed that the design solution space is N M in our methodology. N is a constant, M is the polygon joint number. When M is as small as possible, N M << D M . Therefore, our desired skin design methodology decreases the design solution space greatly. There are some limitations in our methodology. Since the deformation is limited in some range in practice the desired skin curve is not arbitrary large in y extent. It is noticed that the distributed loads are not design variables in our methodology, we assume they are uniform. Therefore the desired skin curve is supposed to be single y-value one. This single y-value curve means that there is only one single intersection point between the desired skin curve and a y- parallel ray shooting from x axis. In addition, the desired skin curve is supposed to be tangent to x axis at both ends of the 2D domain. And this desired skin curve is not expected to be concave in the middle of the curve. 3.7.3 Case Study 108 The following case is a simplified prototype from the stool, chair, cushion or running shoe design. We simplify these applications into a 2D strip domain issue. And the external loads are distributed onto this domain. We assume these loads are uniform. 3.7.3.1 Isosceles triangle like curve skin design Step 1: Figure 3.48 shows the 2D domain configuration with both ends constrained. Distributed loads and N-level printed digital material vector is selected properly. The N-level printed digital material vector is obtained and selected carefully from the digital material design discussed in the following sections. 1 1000 soft E E 496.5[Pa] 0.005[atm]  1.4 1.6 3 7.7 23 235 48mm Thickness t = 6mm 0.8mm ] MPa [ 2.68  soft E 33 . 0   TBD TBD TBD TBD TBD TBD TBD TBD TBD TBD TBD TBD TBD TBD TBD TBD TBD TBD TBD TBD TBD TBD TBD TBD Figure 3.48: Skin design domain with N-level printed digital material. Here the dimension of the 2D domain is 48 0.8 6( ) mm mm mm     with two ends fixed. The distributed load is 496.5Pa. And the printed varied stiffness digital material vector is   1 2 8 , , , dm V E E E  ; where 1 2 3 4 , 1.4 , 1.6 , 3 soft soft soft soft E E E E E E E E     5 6 7 8 7.7 , 23 , 235 , 1000 soft soft soft soft E E E E E E E E     . And the target skin curve is given in the Figure 3.49. 109 Figure 3.49: Target isosceles triangle like curve. Step 2: The given target skin curve is approximated by an open polygon with two edges. In fact this approximation is an isosceles triangle with biggest internal angel 150° . The approximation is shown in the Figure 3.50. Figure 3.50: Target isosceles triangle like curve is approximated by an isosceles triangle. 0 ( 150 ) CP   Step 3: The approximated isosceles triangle is projected onto the skin design domain as shown in Figure 3.51. Totally there are 24 D  elements in this skin domain. TBD TBD TBD TBD TBD TBD TBD TBD TBD TBD TBD TBD TBD TBD TBD TBD TBD TBD TBD TBD TBD TBD TBD TBD Figure 3.51: Approximated isosceles triangle is projected onto the skin design domain. 110 Step 4: Those elements in the skin domain with respect to the links or edges are filled with most stiff digital material. And those elements of joints are set to be design variables. The illustration is shown in Figure 3.52. TBD TBD TBD TBD TBD 1 1000 soft E E 1.4 1.6 3 7.7 23 235 Figure 3.52: Stiff elements and design variables are assigned. Here levels of varied stiffness digital materials N=8; and the count of design variables M=2(considering symmetry). The entire design solution space is 8 2 = 64. Compared to the entire solution space 8 24 , 8 2 <<8 24 . Step 5: It ends up with the following optimization problem. Target Design: ) ) 3 5 ( max( ) 6 5 min( 2 2            s.t. 3 int 1 int jo jo E E  } , , { , , 8 2 1 3 int 2 int 1 int E E E E E E jo jo jo   (3.8) Where 3 int 2 int 1 int , , jo jo jo E E E are with respect to the Young’s modulus of three joint elements. The optimization module in above target skin shape design has many choices, i.e. brute force, GA, PSO and so on. Here we chose GA algorithm to find the solution for this optimization problem. In order to translate this problem into expression by GA we used a 6-digit binary number to represent the Young’s modulus ID of joint-1 and joint-2. Each joint has three binary 111 digits to represent up to 8 Young’s modulus from the digital material library. The population in each generation is set to be 32; generation is set to be 32; mutation rate is set to be 0.001. The results of gnome is 000001, Since 3 int 1 int jo jo E E  , we obtained the Young’s modulus of these three TBD elements shown in the Figure 3.53: Joint 1 Joint 2 Joint 3   8 . 1 5 1  1 1000 soft E E 1.4 1.6 3 7.7 23 235 Figure 3.53: Optimized material assignment result of isosceles triangle with given angle 150° . The displacement at y direction is plotted as shown in Figure 3.54: Figure 3.54: Optimized deformation result of isosceles triangle with given angle 150° . It is easy to verify that the error to the target isosceles triangle is (151.8-150)*100/150 = 1.2 (%). 3.7.3.2 Isosceles trapezoid like curve skin design This example shows how to design an isosceles trapezoid like curve skin. For simplicity only key words in each steps and figures are listed here. Step 1: 2D domain configuration, N-level printed varied stiffness digital material vector and target skin shape as shown in Figure 3.55. 112 1 1000 soft E E 496.5[Pa] 0.005[atm]  1.4 1.6 3 7.7 23 235 48mm Thickness t = 6mm TBD 0.8mm ] MPa [ 2.68  soft E 33 . 0   TBD TBD TBD TBD TBD TBD TBD TBD TBD TBD TBD TBD TBD TBD TBD TBD TBD TBD TBD TBD TBD TBD TBD TBD Figure 3.55: Skin design domain with N-level printed digital material. And the target skin curve is given in Figure 3.56. Figure 3.56: Target isosceles trapezoid like curve. Step 2: Approximate target skin curve into a polygon as shown in Figure 3.57. Figure 3.57: Target isosceles trapezoid like curve is approximated by an isosceles trapezoid. Step 3: Project the approximated polygon onto the 2D domain shown in Figure 3.58. 113 TBD TBD TBD TBD TBD TBD TBD TBD TBD TBD TBD TBD TBD TBD TBD TBD TBD TBD TBD TBD TBD TBD TBD TBD TBD Figure 3.58: Approximated isosceles trapezoid is projected onto the skin design domain. Step 4: Assign stiff materials to edge elements and design variable to joint elements shown in Figure 3.59. 1 1000 soft E E 1.4 1.6 3 7.7 23 235 TBD TBD TBD TBD TBD TBD TBD Figure 3.59: Stiff elements and design variables are assigned. Step 5: It ends up solving the following optimization problem. Target Design: ) ) 3 5 ( max( ) 6 5 min( 2 2            s.t. 4 int 1 int jo jo E E  3 int 2 int jo jo E E  } , , { , , , 8 2 1 4 int 3 int 2 int 1 int E E E E E E E jo jo jo jo   (3.9) 114 The results of gnome is 010110, we obtained the Young’s modulus of these three TBD elements shown in Figure 3.60: 1 1000 soft E E 1.4 1.6 3 7.7 23 235   149.2  Joint 1 Joint 2 Joint 3 Joint 4 Figure 3.60: Optimized material assignment result of isosceles trapezoid with given angle 150° . The displacement at y direction is plotted as shown in Figure 3.61: Figure 3.61: Optimized deformation result of isosceles trapezoid with given angle 150° . 3.8 Summary We present the digital cell representation consists of N by N elements to design varied isotropic or anti-isotropic material in 2D. And we present the digital cell representation consists of N by N by N elements to design varied isotropic or transversely-isotropic material in 3D. It shows this digital cell approach is very fundamental and simple to be implemented. And we build the mathematical model and simulation setup to run, analyze and evaluate the Young’s modulus of a digital material either 2D or 3D. This computational framework for N-level digital 115 heterogeneous material is a powerful tool helping us study the material stiffness property. We creatively present multi-levels to interpolate a digital cell with digital cell pattern from previous level. This iterative interpolation strategy or scheme makes a digital cell to multiply by it. And this scheme enables to interpolate high resolution of digital material from a soft material and a hard material. We present monotone material property increasing rules to choose basis patterns at each level to enable the valid interpolation of two neighbor digital material patterns. These two rules are shown to be self-closed to guide to design monotone increasing N-level digital material. And the trend of physical experimental design matches simulation results very well. We presented a design method to design desired flexible surface based on N-level digital materials. This design approach decreases the entire solution space in a gigantic way based on limited number of joints design variables. 116 Chapter 4 Design Deformable Geometries with Target Compliance As mentioned in the previous section. One approach is to distribute meso-structures to achieve target compliance geometry design. Here meso-structures are made from a single stiff material. Our research interest is to employ single stiff material to design desired compliance. Recent advances in sold freeform fabrication (AM) present tremendous design freedom for a product design with complex geometries. In this paper, we consider the problem of using rigid materials to design flexible surface of a product component for AM. A design strategy based on the combination of well-defined meso-structures is presented to achieve desired heterogeneous material properties, and consequently desired flexibility in target directions and positions. We present our computational framework to automate the design optimization process. Due to the dramatically increased design space, a brute force or traditional design optimization method such as the genetic algorithm (GA) and particle swam optimization (PSO) is not efficient. We present a design method based on the idea of analyzing the flexibility of each link for given meso- structures. Two experimental examples are presented to demonstrate its usage in generating the maximum/minimum and target displacements. We also present its comparison with the GA and PSO methods. 4.1 Desired Compliance by Performance-tailored Meso-structures In our method, we are trying to tailor each meso-structure by adjusting its type and parameters to achieve desired properties of a product component. Such performance-tailored meso-structures provide an alternative way to design product components with desired flexibility. A component with such performance-tailored meso-structures will definitely be more complex. However, such a complex design may be beneficial since as stated in the Law of Requisite Variety(Ashby 1956), the more states (variety) a system is capable of possessing, the more 117 requirements it is capable of achieving and the more disturbances it is capable of handling. By observing biological systems, Ashby concluded that in order for a system to be capable of surviving in a complex environment, the level of complexity of the system must increase to at least the same level as its environment. Thus an increase in the complexity of the system is critical for its robustness and adaptability. Similarly a component design must be more complex if more design specifications are required to define it. The introduction of performance-tailored meso-structures in a product component can dramatically increase the design freedoms that are available for a designer. Consequently, we believe performance-tailored meso-structures provide a potential solution to achieve heterogeneous material properties and consequently multi-functionality in a component. Figure 4.1: A feasibility study of achieving multiple design requirements (i.e. pushing up/down and rotating) by meso-structures. A feasibility study based on the CAD model of a bunny has been performed (refer to Figure 4.1). We used two types of meso-structures (rigid beams and flexible springs) to design a flexible 118 skin. We designed the CAD model and built it by using a AM process, selective laser sintering (SLS). Based on the built physical model, we believe the combination of such a large number of meso-structures provides us sufficient design freedom to achieve multiple types of movements. For example, we can push up and down the bunny’s head as shown in Figure 4.1. We can also rotate its head in various directions. By changing some meso-structures, we can achieve different flexibility requirements in different locations. We may also add certain degree of robustness in the design such that the desired flexibility will not be affected even if one spring or beam is broken. Therefore, by using an extra layer of meso-structures in a component design, we will lead to a higher level of complexity, which correlates to a higher level of variety. The intelligently design of such combination of meso-structures may allow a component design to satisfy complex design requirements. 4.2 Formulation of Flexible Skin Design Problem In this proposal, we investigate the flexible skin design problem based on performance- tailored meso-structures. Assume constraints and n loads acting on a component skin are given. An external load can be either a force or a moment. A product’s skin is usually defined by a boundary representation (B-rep). Therefore, such skin can easily be approximated by a set of nodes and links. In this paper, we assume a set of nodes N i and links L j between N i have been given. We will only consider L j as the candidates for performance-tailored meso-structures. Further optimizing the positions of N i and topology of L j (i.e. node connections) can give us an even bigger design space for a flexible skin design. There are many different types of meso-structures that can be considered including the research results generated by the topology optimization (Wang 2005; Chen 2007). In this 119 proposal, we just use two types of meso-structures, a rigid beam B j and a spring S j , since they are simple and easy to understand. Yet a large combination of such simple elements can have sufficient capability to achieve desired flexibility. This is similar to the black/white pixels used in a digital image. Even though we have only two bits (0/1), their combinations of a large set of bits can draw a wide variety of images. In our study, we assume all the rigid beams will have the same size (we set r = 0.65mm in our study) since changing their sizes has a relatively small effect on the flexibility of a component’s skin. A spring can be defined by a set of parameters: radius of spring wire r, outer radius of spring R, lead of coil L coil , and total spring length L spring . To simplify the problem, we fix r, R, and L coil based on our experience of the target manufacturing process (for SLS process, we set r = 0.65mm, R = 1.35mm, L coil = 2.1mm) and only change the spring length L spring . Suppose the length of a link between two nodes is L Beam . We define the spring length ratio  = L spring / L Beam . Therefore, we can tailor the flexibility of a spring by changing its . In our study, we set the range of the ratio  from 1 / 3 to ½. The inclusion of more spring variables will only give us more design freedoms without affecting our design method. Therefore, the flexible skin design problem considered in this paper aims at using a combination of beams B j (fixed) and springs S j (with some variations) to achieve desired flexibility at target positions. Due to the tremendous design freedom, a design method based on the meso-structures may provide novel designs for various design problems such as: (1) Type I: the minimum or maximum displacements; (2) Type II: target displacements. These two types of flexible skin design problems have been mentioned in the previous chapter. 120 In this proposal, we present an optimization method based on the flexibility analysis of the links in a given skin model. The computational framework and related optimization techniques are discussed for the aforementioned design problems. 4.3 Computational Framework for Flexible Skin Design A large number of meso-structures pose significant challenges in analyzing, synthesizing, and optimizing a component’s skin design to meet the given flexibility requirements. In this research, we will only consider elastic deformations under static loads. Although a wealth of literatures exist on analyzing structures with beams and springs under such conditions, we decide to use a commercial finite element analysis (FEA) system for computing the deformation of meso-structures under loads. A well-accepted FEA system can provide us general and reliable simulation results for a wide variety of displacements and loads. Figure 4.2: A computational framework for flexible skin design. We use COMSOL Multiphysics 3.2 (www.comsol.com) as the FEA tool in our computational framework. Based on a provided programming language, COMSOL script, we develop a set of automation tools for the analysis process. The computational framework we set up for the flexible skin design is shown in Figure 4.2. For a given component skin, we first define a script design pattern which incorporates design requirements such as base material Refined Design Configuration FEA Result Analyzer Analysis Script Generator Script (.M) SFF Machines CAD Generator COMSOL Script COMSOL Multiphysics Post-processing Data Design Optimization Initial Design Configuration Script Design Pattern Flexibility Requirements Product Skin Design 121 properties, supports, and loads. According to the design pattern, a software module, analysis script generator, will automatically generate a script file (.M) for a meso-structure configuration. The script file will be analyzed by COMSOL Multiphysics. We run the simulation in the batch mode such that the FEA analysis can be performed in the background without any user interactions. After the simulation, COMSOL Multiphysics will save post-processing data in a specified format at a given file path. Accordingly, we developed another software module, FEA result analyzer, to extract target displacements related to the flexibility requirements. The analyzed results will be sent to a design optimization module, which will be discussed in more details in the next section. Based on the FEA results, the design optimization module will determine how to refine the configurations of meso-structures for further analysis. Finally, an optimized design configuration will be generated for achieving the desired flexibility requirements. A software module, CAD model generator, will then construct a valid polygonal model (.STL) for the skin with meso-structures based on the design configuration (Chen and Wang 2008). The model can then be built in a AM system for physical models. In the analysis of meso-structures consisting of beams and springs, we use 3-dimensional Euler beams in the FEA model. This is based on the consideration that both bending and torsion in addition to axial loading exist in such meso-structures. For a spring that is mathematically well-defined, we use a large number of small beams to approximate it. Therefore, the FEA model we send to COMSOL Multiphysics actually contains a set of 3-dimensional Euler beams. In order to ensure the flexibility analysis results are accurate, we select various test cases including examples provided in textbooks and by using another FEA software system (ANSYS) for verification. In all the test cases, the FEA simulation results generated by our automation tools are satisfactory with small differences (generally less than 5%). In our FEA model, we did not 122 consider conditions such as yield failure, fracture failure or buckling. They can be added by checking the stresses of each Euler beam based on the FEA results. (a) (b) Figure 4.3: An example of finite element analysis results generated by our automation tools. We also test the simulation time of our automation tools. Not surprisingly, the running time for generating the analysis results is proportional to the number of Euler beams in a FEA model. The running time is generally satisfactory. For example, two FEA results generated by our automation tools are shown in Figure 4.2. A color scheme is also given to show the deformations. In Figure 4..(a), the bunny model has all rigid beams (beam number is 220); in Figure 4.3.(b), some of the beams have been replaced by flexible springs which are approximated by smaller Euler beams (beam number is 1,750). The supports and loads are shown in Figure 4.15. It takes around 22 and 31 seconds respectively for generating the analysis results by using a PC with Pentium D 3.2GHz and 2GB DRAM. The same base material (Nylon) and amount of forces are used in both simulations. It can be seen from the results that the maximum displacements will significantly increase (by ~10 times) after replacing some beams with springs. 4.4 Optimization Methods for Flexible Skin Design 123 The flexible skin design in our work is a complex combinatorial problem, in which the computational expense exponentially increases with the problem scale. Worst of all, we have no explicit objective function related to design variables. Consequently, this type of problem is usually an NP hard problem and will need some type of heuristics in its optimization process. To perform the optimization of a flexible skin design, we investigate various optimization approaches including the state-of-the-art meta-heuristics methods such as Genetic Algorithm(Holland 1992) and Particle Swarm Optimization(Kennedy and Eberhart 1995; Kennedy, Eberhart et al. 2001). (1) The genetic algorithm (GA) method is based on the natural selection by mimicking the evolution process of biological systems. GA first randomly generates a generation of individuals, then removes bad individuals and produces better and better ones. GA is a probabilistic method and only relies on the objective value; it has no requirement on the properties about the objective function like continuity, differentiability, convexity, etc. (2) The particle swarm optimization (PSO) is another kind of evolutionary computing technology to mimic the social behavior of animals such as fish schooling and bird flocking. The biologic colonies demonstrate that intelligence generated from complex activities such as cooperation and competition among individuals can provide efficient solutions for various optimization problems. Based on the investigation, we found the aforementioned meta-heuristics methods are both population based methods. They can solve small scale problem efficiently. However, the quality and efficiency of the solution is much worse for more complex models and design requirements (refer to a comparison of the test results as shown in Table 4.1 and 4.2). For several simple design cases, we use the computational framework as shown in Figure 4.2 to analyze all the combinations of rigid beams and flexible springs. Such a brute force method 124 may take days for even a small scale problem (e.g. for a structure with only 20 links, the total combination number of beams and springs would be 2 20 = 1,048,576). However, the exhausted search can provide us a benchmark for a comparison with the optimization results of other methods. An analysis of the exhausted search results also provides us some interesting findings. (1) The flexibility of each link will have a different impact on the desired flexibility at a target position. That is, some links will have bigger impacts to the target displacements while some other links will have neglectable impacts. (2) The flexibility of a link may inversely affect the flexibility at a target position. That is, when replacing some beams by springs for those links, the target displacement will decrease instead of increase. (3) It is important to determine the flexibility of each link such that a desired heterogeneous material property can be achieved for the target flexibility in a skin design. Figure 4.4: A test case of a 10-beam model. Figure 4.5: A beam example with negative contribution factor. An example of a simple 10-beam structure (B 0 ~ B 9 ) is shown in Figure 4.4 to demonstrate such a phenomenon. Suppose four beams B 4 ~ B 7 are constrained by a bottom plate. An external force is added at node N 1 along the direction of beam B 8 , and we are interested in the displacement of node N 2 . The analysis results indicate that by replacing B i by S i for each link, the displacement at N 2 will generally increase except for B 8 and B 9 . That is, the replacement of B 8 by Target Point External force B 7 B 4 B 5 B 6 B 8 N 1 N 2 B 9 B 1 B 2 B 0 B 3 N 3 N 4 Disp = 0.782mm Disp = 0.768mm 125 a spring will reduce the displacement at N 2 instead of increasing it (refer to the results as shown in Figure 4.5). A possible explanation is that the spring at B 8 will store some deformation energy, and, consequently, reduce the displacement at N 2 . Beam B 9 is not connected to any constraint or external load. Therefore, B 9 will have no effect on the displacement of N 2 . The flexibility analysis of each link can provide us helpful heuristic in the flexible skin design. In the next section, we will present a technique called flexibility contribution factors (FCF) for such analysis. Correspondingly, a new heuristic method based on FCF is developed and presented in Section 4.6 for the design problems of achieving target displacements under loads. 4.5 Flexibility Contribution Factors Methodology In the Type I design problem as shown in 4.2, we are interested in the size of a displacement without considering its direction. This is common in applications such as vibration controls, in which, we are interested in the displacement of a target point due to a vibration source. It is well known that using some flexible materials can reduce the broadcasting of vibration. In this study, we further demonstrate that a combination of rigid beams and flexible springs can reduce or magnify the displacements of one or multiple target points. Although we only consider static loads in our study, we believe our method can be extended to dynamic loads that are typically required by vibration controls. 4.5.1 Definition of Flexibility Contribution Factor In our study, a link L i can be either a beam (B i ) or a spring (S i ). Its stiffness can be measured by the well-known spring constant K Li = F/x as shown in Hooke’s law. We can also define C Li = 1/K Li = x/F as a measurement of the flexibility of L i . A quantitative study indicates that a beam in our study has K Li  70,000 N/m; while a spring with the spring length ratio  = 1 / 3 126 and ½ has K Li  2,400 N/m and 1,600 N/m respectively. Therefore, stiffness K Bi >>K Si while flexibility C Bi <<C Si . By changing a link from a beam to a spring, its flexibility C Li increases dramatically (30~40 times). This is because different physical principles are dominant in a beam and a spring respectively. However, for a spring with different , its flexibility change will be much smaller (~30%). Therefore, in the flexibility analysis of a link, we will focus on the load- deformation behavior difference between a beam (B i ) and a spring (S i ). Suppose all the links of a 3D skin model are represented as rigid beams B i which have flexibility C Bi . For given loads and constraints, the displacement of a target point T k has a size of 1 () k i n T B B B C C C  . By replacing B i by a spring S i , the flexibility of the link is increased to C Si . Due to such a single replacement in the skin model, the displacement of T k will have a size of 1 () k i n T B S B C C C  under the same loads and constraints. Suppose the displacement size of  Tk is the output that we are interested in the flexible skin design. We can define the output-to- flexibility amplification (deamplification) relationship based on () k i i T B S CC   = 1 () k i n T B S B C C C  - 1 () k i n T B B B C C C  , which is denoted as the flexibility contribution factor of a link L i (denoted as FCF Li ). Therefore FCF Li is a qualifier that attaches a number to a specific link for characterizing the link by its capacity of creating desired target displacements. Based on the flexibility contribution factor, we can classify all the links into three categories: (1) Positive contribution factors: FCF Li > 0; (2) Negative contribution factors: FCF Li < 0; (3) Neutral contribution factors: FCF Li  0. We can further sort FCF Li based on their sizes for the same amount of flexibility changes of links (i.e. C Si - C Bi ). For example, for the 10-beam model as shown in Figure 4.4, we use a spring 127 with length ratios  = 1 / 3 to compute FCF Li . The computed flexibility contribution factors are shown in Figure 4.6. In the figure, we use various gray scale levels to represent the size of FCF Li . A link L i is drawn in a darker line if its |FCF Li | is larger; and it is drawn in a dotted line (e.g. B 8 ) if its FCF Li < 0. 4.5.2 Properties of Flexibility Contribution Factors Figure 4.6: Flexibility contribution factors of the 10-beam model. We study the flexibility contribution factors and identify some of its properties. (1) Different target points. Figure 4.7: Flexibility contribution factors of the 10-beam model for different target points. Target Point Negative contribution factor Positive contribution factor Target Point Target Point Negative contribution factor Positive contribution factor 128 The flexibility contribution factor of a link is closely related to the location of a target point. Their values including their signs and relative sizes will be dramatically different for different target points. For example, for the 10-beam model as shown in Figure 4.6, if we change the target point from N 2 to N 3 or N 4 , we will get different sets of flexibility contribution factors, which are shown in Figure 4.7. Notice there are no negative contribution factors among them. (2) Different skin topology, load positions, load directions and constraints. Similar to different target points, the flexibility contribution factors are closely related to the topology of links, the positions and directions of loads and constraints. (3) Different load sizes and base materials. Assuming that n loads act on an elastic body, the deformation at a specific location T k can be found as 1 k n T j j j CF     , where C j are flexibility coefficients. Based on a principle for the linearly elastic bodies called principle of linear superposition (Lobontiu 2003), the deformation at any point in an elastic body under the action of a load system is equal to the sum of deformations produced by each load when acting separately. Therefore, we only need to consider one load in studying its properties. The test results on the 10-beam model indicate that, for different load sizes or base materials (e.g. Nylon or steel), (1) the sign of FCF Li will not change (i.e. the links with positive/neutral/negative contribution factors will remain in the same sets); (2) the ranked orders of FCF Li based on their sizes will not change. For the 10-beam model, an example of the flexibility contribution factors for different load sizes is shown in Figure 4.8. In the figure, the number of beams in the X axis refer to B 0 ~ B 9 as shown in Figure 4.4. 129 Figure 4.8: Flexibility contribution factors of the 10-beam model for different load sizes. (4) Different link flexibility (C Si ). In computing FCF Li , we can use a spring with different flexibility C Si to replace each link. However, the signs of FCF Li and the ranked orders based on their sizes will not change for different C Si . For example, in the 10-beam model, we use springs with three different length ratios (  = 1 / 3 , 2 / 5 and ½) as C Si to replace C Bi for computing FCF Li . The computed flexibility contribution factors are shown in Figure 4.96. In the figure, the number of beams in the X axis refer to B 0 ~ B 9 as shown in Figure 4.4. Figure 4.9: Flexibility contribution factors of the 10-beam model for different link flexibility. 130 4.6 Target displacements by Flexibility Contribution Factor Method Based on FCF Li , we present a new constructive heuristics method, Flexibility Contribution Factor (FCF) Method. 4.6.1 Flexibility Contribution Factor Method Our FCF method contains two major steps for analyzing the relation between the flexibility of a link and a target displacement. (1) Computing FCF Li . The equations and related descriptions for computing the flexibility contribution factor FCF Li have been presented in Section 4.5.1. As the result, we will have three sets of links, positive, negative and neutral contribution factors. In addition, all the links in each set can be easily sorted based on their sizes. The computation complexity of computing FCF Li is linear ( ) (n O ). Therefore, the computation time is satisfactory. For example, for a bunny model as shown in Figure 4.64 which has 220 links, it takes ~1.2 hour to compute FCF Li for all the links. (2) Computing accumulated displacement of FCF Li . For a desired flexibility in a target position, we can determine the order for sorting FCF Li . For example, for achieving the maximum displacement, FCF Li can be sorted from the biggest to the smallest considering their signs; while for achieving the minimum displacement they can be sorted in the reverse order. Therefore, based on the sorted FCF Li , we can compute the accumulated displacement of FCF Li by using the following steps: sort all the links according to their FCF Li set all the links as beams for each link L i in the sorted order replace beam B i by a spring S i calculate  Tk as the accumulated displacement of L i end for 131 Figure 4.10: Accumulated displacement with minimal and maximal spring ratios for the 10-beam model. For FCF Li as shown in Figure 4.10, the computed accumulated displacements are shown in Figure 4.12. The test results indicate that the displacement  Tk keeps increasing with each replacement of a beam with positive contribution factor; while  Tk keeps decreasing with each replacement of a beam with negative contribution factor. In addition, among the positive contribution factors, the displacement increases are bigger for the first few replacements. A test of a more complex model (refer to Section 4.7) indicates a similar trend. Therefore, the flexibility contribution factor method is effective. The accumulated displacement of FCF Li is a useful tool for the flexibility analysis of links. Its computation is also efficient. The complexity of computing the accumulated displacement is linear ( ) (n O ). Even thought the sorting of FCF Li is ) lg ( n n O , it generally takes a short time (<<1 second) since all related computations are rather simple. 4.6.2 Flexible Skin with Multiple Maximum/Minimum Displacements Assuming that the deformation at a location T k under a given loading condition is  Tk . For a multi-objective optimization problem with multiple target displacements, the math formulation of the maximum/minimum displacements is given as follows: 132 GIVEN: Nodes N i , links L j , targets T k , loads, and constraints FIND: Beams B j , springs S j with length ratio  Sj . SATISFY: L j = B j S j , B j S j = 0  Tk = disp Tk (B j , S j )  min   Sj   max MAXIMIZE: 0 k m kT k    . Suppose appropriate weights  k have been assigned for the multiple objectives with 0 | | 1 m k k     . In addition,  k > 0 for the maximum displacements; and  k < 0 for the minimum displacements since the objective in the math formulation is set as maximize. The steps of our design method based on FCF Li are given as follows. Step 1: Compute FCF Lj (  Tk ) for each location T k ; Step 2: Compute FCF Lj = 0 () k m k Lj T k FCF    for each link L j ; Step 3: Sort FCF Lj based on its size; Step 4: Set L j as S j with  Sj =  max if FCF Lj > 0; otherwise set L j as B j . We apply the FCF method to optimize the maximum/minimum displacement of node N 2 in the 10-beam model. After computing FCF Lj (  N2 ), (1) for the maximum displacement of N 2 , its weight = 1. We simply replace all the links whose FCF Lj (  N2 ) > 0 with springs  Sj =  max . Consequently an optimized result can be generated as shown in Figure 4.11(a); (2) for the minimum displacement of N 2 , its weight = -1. We can replace all the links whose FCF Lj (  N2 ) < 0 with springs  Sj =  max . The related result is shown in Figure 4.11(b). 133 Figure 4.11: The optimized 10-beam models (maximum vs. minimum displacements). We also use a hybrid optimization method by integrating GA and PSO to calculate the maximum displacement of N 2 . We first use GA by fixing  Sj = 1 / 3 to try different combinations of beams and springs. Based on the resulted combination, we then use PSO to optimize  Sj of all the springs. The search space of  Sj is set as ( 1 / 3 , 1 / 2 ). The converging curve of the maximal displacement according to the evolutionary process is shown in Figure 4.12. It turns out that the FCF method can generate a better result with much shorter time. The running results are given in Table 4.1. Figure 4.12: Optimized maximum displacement by GA and PSO. Table 4.1: Results of optimization methods (10-beam model). Optimization method FCF GA GA+PSO Disp. = 1.205mm Disp. = 0.766mm (a) (b) 134 Max displacement (mm) 1.20 0.92 1.18 CPU time 2min 50min 90min The displacement of N 2 in the 10-beam model without any springs is 0.782mm. Therefore, the replacement of B 8 by a spring S 8 will reduce its displacement; hence the flexibility of S 8 will contribute to the damping of the vibration from N 1 to N 2 . When we increase  S8 from 1 / 3 to 1 / 2 ,  N2 decreases from 0.768mm to 0.766mm. Ideally if we keep on increasing the flexibility of S 8 , the displacement  N2 will keep decreasing. When the flexibility of S 8 is set as infinity, the link between N 1 and N 2 is essentially deleted. This leads to a new model with only 9 beams. We analyze the 9-beam model (with B 8 deleted). The displacement  N2 = 0.618mm, which verifies our analysis. 4.6.3 Flexible Skin with a Target Displacement A more challenging optimization problem is to achieve a target displacement d Tk which is not necessarily the maximum or minimum displacement at T k . The math formulation of the target displacement problem is given as follows: GIVEN: Nodes N i , links L j , target T k and d Tk , loads, and constraints FIND: Beams B j and springs S j with length ratio  Sj . SATISFY: L j = B j S j , B j S j = 0  Tk = disp Tk (B j , S j )  min   Sj   max MINIMIZE: || kk TT d   . The steps of our design method based on FCF Li for achieving the target displacements d Tk are given as follows. 135 Step 1: Compute FCF Lj (  Tk ) for T k ; Step 2: Compute the accumulated displacements (Acc-Disp) of  Tk based on springs with  Sj =  min and  max respectively (refer to Figure 4.71); Step 3: Identify a minimum set of links (L m ) such that d Tk is between the accumulated displacements Acc-Disp min and Acc-Disp max (i.e. Acc-Disp min d Tk Acc-Disp max ); Step 4: Use particle swarm optimization (PSO) method to optimize the spring settings of S m for a minimum difference between  Tk and d Tk . We apply the FCF method to optimize the target displacement of node N 2 in the 10-beam model. Suppose the target displacement d N2 is 1.0mm. Based on FCF Lj (  N2 ), we can compute the accumulated displacements Acc-Disp min and Acc-Disp max , which are shown in Figure 4.13(a). For d N2 = 1.0mm, we can identify a set of candidates with the minimal number of links (B 4 , B 6 , B 7 ) such that Acc-Disp min d Tk Acc-Disp max . Therefore for a corresponding structure as shown in Figure 4.13(b), it has sufficient design space for achieving the target displacement. We further use PSO to optimize the spring length ratios of the identified springs. The optimization result based on PSO is shown in Figure 4.14. After around 10 iterations, the difference between  Tk and d Tk converges to around 0. Figure 4.13: An example of computing a target displacement based on FCF and accumulated displacement. B 4 B 6 B 7 L m (a) (b) 136 Figure 4.14: Optimized distance from the target position for PSO. Therefore, in essence, we use FCF method to select a set of candidate links such that the given target displacement can be achieved based on the selected links. However, compared to the full optimization of all the links, our method can significantly reduce the number of variables required in an optimization method. For example, in the case as shown in Figure 4.13, only 3 springs, instead of 10, will be considered in PSO; in the case as shown in Figure 4.20(b), only 4 springs, instead of 220, will be optimized. Consequently the optimization process can be run much faster, especially for complex models. 4.7 Experimental Results and Analysis A Flexibility Design Exploration Testbed has been implemented using C++ programming language with Microsoft Visual C++ 2005 complier. In addition to the FEA (COMSOL Multiphysics) and the automation tools as shown in Figure 4.2, we also integrate various optimization modules based on the methods including flexibility contribution factor, greedy, GA and PSO in our Testbed. The optimized result can be dynamically shown in the Testbed to assist the data analysis. 137 Figure 4.15: A complex test case: a bunny model. We present the test results of a more complex model in this section. As shown in Figure 4.15, a bunny model with 220 beams is considered. Suppose the model is constrained at its bottom. Two forces (F 1 and F 2 ) are provided to twist its head. The target displacement we are interested in is at the end of one ear (T 1 ). Based on the procedure described in Section 4.5.1, we can compute the flexibility contribution factors FCF Li for all the links (we set springs at  = 1 / 3 ). The computed FCF Li are shown in Figure 4.16. In the figure, a link is drawn in a darker line if its |FCF Li | is bigger; in addition, a link with negative contribution factors is drawn in a dotted line. T 1 138 Figure 4.16: Flexibility contribution factors for the bunny model. We also compute the flexibility contribution factors FCF Li for another target point T 2 . The results are shown in Figure 4.17. They are quite different from those as shown in Figure 4.16. Figure 4.17: Flexibility contribution factors at a different target. The sorted FCF Li and related accumulated displacements for target point T 1 are shown in Figure 4.18 and Figure 4.19 respectively. Target Point Target Point Negative contribution factor Positive contribution factor T 1 T 2 139 Figure 4.18: Sorted flexibility contribution factors for the bunny model. Figure 4.19: Accumulated displacement for the bunny model. For the maximal displacement at T 1 based on the same spring setting (  = 1 / 3 ), we use the GA, greedy and FCF methods to compute the optimized flexible skin respectively. GA method: We use a binary string to represent the parameters: bit “1” denotes the beam is changed with a spring and bit “0” denotes the beam is untouched. For the genetic operators, we choose uniform crossover and one point mutation operator, roulette wheel selector. The probability of crossover is set at p c =0.9; the probability of mutation is set at p m =0.01; the population size N ga =30; and the maximal cut off generation is MaxGen ga =100. 140 Greedy method: In this method, we calculate the flexibility contribution factors for all the links that are still set as beams. Base on the computed results, we identify the best link L i among the candidates whose replacement by a spring S i can give us the maximum displacement at T 1 . After L i is set as S i , we update the FEA model accordingly and remove L i from the candidate links. For the remaining candidates, we repeat the aforementioned process by recalculating their flexibility contribution factors under the new settings. The optimization process will continue until no replacement of S i can increase the displacement. Therefore, unlike the FCF method, in the greedy method we will update the flexibility contribution factors in each iteration. The updating scheme will dramatically increase the computation time. However, the interactions between FCF Li will be considered in the greedy method. FCF method: Based on the steps described in Section 4.6.2, we compute FCF Lj (  T1 ) and replace all the links whose FCF Lj (  T1 ) > 0 with springs. The model generated by the FCF method respect to the maximum displacement is shown in Figure 4.20(a), which leads to  T1 = 3.24mm. Figure 4.20: The optimized bunny models by FCF method (maximum vs. minimum displacements). The running results by the three aforementioned optimization methods are shown in Table 4.. The results illustrate that the flexibility contribution factor (FCF) method is both effective and efficient. Disp. = 3.24mm Disp. = 0.08mm Target Point Target Point (a) (b) 141 Table 4.2. Results of optimization methods (bunny). Optimization method FCF GA Greedy Max displacement (mm) 3.24 2.74 3.30 CPU time 1h 25h 100h We also use the FCF method as described in Section 4.6.3 for a target displacement of 0mm at T 1 . The model generated by considering the first four links which have the biggest negative contribution factors is shown in Figure 4.20(b). The displacement is 0.08mm, a significant decrease compared to 0.219mm for the all-beam model. The minimum displacement can be further reduced by using PSO method to refine the spring settings. 4.8 Summary Ever-more challenging design requirements require the usage of heterogeneous material properties within a product component. We illustrate an important design strategy by using performance-tailored meso-structures to achieve desired heterogeneous material properties. A product component with complex meso-structures can be fabricated by solid freeform fabrication processes. In a flexible skin design with only two types of meso-structures (beams and springs), both meso-structures are quite simple. However, a large number of their combinations require better design exploration methods. To effectively get a good solution for the desired flexibility, we proposed a constructive heuristics method based on the flexibility analysis of all the links. Related techniques for such purpose are presented including flexibility contribution factors and accumulated displacement. Compared with other meta-heuristics methods such as genetic algorithm and particle swam optimization, the flexibility contribution factor method can 142 dramatically improve the efficiency without loss of the quality. Moreover, it can provide us a group of candidate links for meso-structure design. 143 Chapter 5 Conclusion This research, which is focused on controlled flexibility skin design, is expected to make the following contributions: We present our computational framework to automate the design optimization process of flexible skin design based on beam/spring switch meso-structures. Based on the energy principles and the uniform strain energy density proposition in truss discrete structures, we extended and proved the axiom of strain energy density to the beam structures. We present a numerical method to compute the principal stress lines of the stress field in a continuum design domain. We propose a reduction-expansion design method to approximate the principal stress lines(PSL) in a continuum design domain to a discrete Michell-like topology optimized beam structure, which has minimum compliance with given volume of material. We also extend our PSL topology optimization method to maximal permission stress, given displacement and 3D case. We present a design method to design desired flexible skin based on N-level stiffness digital materials. This design method decreases the entire solution space in a gigantic way based on limited number of joints design variables. We present the digital cell representation scheme consists of N by N elements to design varied isotropic or anti-isotropic material in 2D. We present the digital cell representation consists of N by N by N elements to design varied isotropic or transversely-isotropic material in 3D. We build the mathematical model and simulation setup to run, analyze and evaluate Young’s modulus of a digital material either 2D or 3D. 144 We creatively present multi-levels to interpolate a digital cell with digital cell pattern from the previous level. This iterative interpolation strategy or scheme makes a digital cell to multiply by it. And this scheme enables to interpolate high resolution of digital material from a soft material and a hard material. We present monotone material property increasing rules to choose basis patterns at each level to enable the valid interpolation of two neighbor digital material patterns. A Flexible Contribution Factor design strategy based on the combination of well-defined meso-structures is presented to achieve desired heterogeneous material properties, and consequently desired flexibility in target directions and positions. We prove the correctness of the optimum condition in a beam structure that the entire beam structure achieves equilibrium and optimum of material distribution when strain energy density is uniform. 145 Bibliography 1. Achtziger, W. and M. Stolpe, "Global optimization of truss topology with discrete bar areas - Part II: Implementation and numerical results," Computational Optimization and Applications, Vol. 44, pp. 315-341, Nov. 2009. 2. Achtziger, W., "Topology optimization of discrete structures: An introduction in view of computational and nonsmooth aspects," Topology Optimization in Structural Mechanics, pp. 57-100. Springer, Vienna, 1997. 3. Achtziger, W., "On simultaneous optimization of truss geometry and topology," Structural and Multidisciplinary Optimization, Vol. 33, pp. 285-304, Jan. 2007. 4. Achtziger, W., M.P. Bendsø e, A. Ben-Tal, and J. Zowe, "Equivalent Displacement Based Formulations for Maximum Strength Truss Topology Design," Impact of Computing in Science and Engineering, Vol. 4, pp. 315-345, Dec. 1992. 5. Achtziger, W. and M. Stolpe, "Global optimization of truss topology with discrete bar areas - Part I: Theory of relaxed problems, " Computational Optimization and Applications, Vol. 40, pp. 247–280, Nov. 2008. 6. Ashby, W. R. (1956). An Introduction to Cybernetics. London, UK, Chapman & Hall Ltd. 7. B. Berman, "3-D printing: The new industrial revolution", Business Horizons, Elsevier, Vol. 55, 2012, pp. 155-162. 8. Beaman, J. J., J. W. Barlow, et al. (1996). Solid Freeform Fabrication: A New Direction in Manufacturing, Springer. 9. Bendsø e, M. P. (1995). Optimization of Structural Topology, Shape, and Material. Berlin Heidelberg, Springer-Verlag. 146 10. Bendsø e, M.P. and A. Ben-Tal, "Truss topology optimization by a displacement based optimality criterion approach," Optimization of Large Structural Systems Vol I, pp. 139-155. Kluwer, Dordrecht, 1993a. 11. Bendsø e, M. P. and N. Kikuchi (1988). "Generating Optimal Topologies in Structural Design Using a Homogenization Method." Computer Methods in Applied Mechanics and Engineering(71): 197-224. 12. Bendsø e, M.P. and C.A. Mota Soares, eds., Topology Design of Structures, Kluwer Academic Publishers, Dordrecht 1993b. 13. Bendsø e, M.P. and O. Sigmund, "Material Interpolation Schemes in Topology Optimization, " Archive of Applied Mechanics, Vol. 69, pp. 635-654, Nov. 1999. 14. Bendsø e, M.P. and O. Sigmund, Topology optimization: theory, methods and applications. Springer, Berlin, 2002. 15. Bendsø e, M.P., Methods for Optimization of Structural Topology, Shape and Material, Springer Verlag, New York 1995. 16. Bendsø e, M.P. and O. Sigmund, Topology optimization: theory, methods and applications. Springer, Berlin, 2002. 17. Ben-Tal, A., M. Kočvara and J. Zowe, "Two nonsmooth approaches to simultaneous geometry and topology design of trusses," Topology optimization of structures, pp. 31-42, Kluwer Academic Publishers, Dordrecht, The Netherlands, 1993. 18. Bickel, Bernd, et al. "Design and fabrication of materials with desired deformation behavior." ACM Transactions on Graphics (TOG) 29.4 (2010): 63. 19. Burns, S.A., ed., Recent advances in optimal structural design, Institute of the American Society of Civil Engineers, Reston, VA, 2002. 147 20. Bus SA, Ulbrecht JS, Cavanagh PR: Pressure relief and load redistribution by custom-made insoles in diabetic patients with neuropathy and foot deformity. ClinBiomech 19:629–638, 2004. 21. Chan, H.S.Y., "The Design of Michell Optimum Structures", College of Aeronautics, Cranfield, United Kingdom, Report No. 142, Dec., 1960. 22. Chen, Y. (2007). "3D Texture Mapping for Rapid Manufacturing." Computer-aided Design and Application 4(6): 761-771. 23. Chen, Y. and S. Wang (2008). "Computer-aided Product Design with Performance-Tailored Mesostructures." Computer-aided Design and Application 5(1-4): 565-575. 24. Chen WP, Ju CW, Tang FT: Effects of total contact insoles on the plantar stress redistribution: a finite element analysis. Clin Biomech 18:S17–S24, 2003. 25. COMSOL,www.comsol.com 26. Dorn, W.S., R.E. Gomory, and H.J. Greenberg, "Automatic design of optimal structures," Journal de Mecanique, Vol. 3, pp. 25-52, 1964. 27. Eschenauer H.A. and N. Olhoff, "Topology optimization of continuum structures: a review, " Applied Mechanics Reviews, Vol. 54(4), pp. 331–389, Jul. 2001. 28. Fredricson, H., "Topology Optimization of Frame Structures—Joint Penalty and Material Selection," Structural and Multidisciplinary Optimization, Vol. 30, pp. 193-200, Apr. 2005. 29. Fredricson, H., T. Johansen, A. Klarbring, and J. Petersson, "Topology Optimization of Frame Structures With Flexible Joints," Structural and Multidisciplinary Optimization, Vol. 25, pp. 199-214, Feb. 2003. 30. Gibson, L. J. and M. F. Ashby (1997). Cellular Solids: Structure and Properties, Cambridge University Press. 148 31. Hajela, P., E. Lee, and C.Y. Lin, "Genetic algorithms in structural topology optimization," Topology Design of Structures, pp. 117-134, Kluwer Academic Publishers, Dordrecht, 1993. 32. Haug, E.J. and J.S. Arora, Applied Optimal Design, Wiley&Sons, New York 1979. 33. Hegemier, G.A. and W. Prager, "On Michell trusses," International Journal of Mechanical Sciences, Vol. 11(2), pp. 209-215, Feb. 1969. 34. Hemp, W.S. and H.S.Y. Chan, "Optimum Design of Pin-Jointed Frameworks", Aeronautical Research Council Reports and Mem. No. 3632, Her Majesty’s Stationery Office, London, United Kingdom, 1970. 35. Hemp, W.S., "Studies in the Theory of Michell Structures," Proceedings International Congress of Applied Mechanics, Munich, West Germany, 1964. 36. Hemp, W.S., Optimum Structures, Clarendon Press, Oxford, United Kingdom, 1973. 37. Hill, R., The Mathematical Theory of Plasticity, Oxford University Press, Oxford, 1950. 38. Holland, J. H. (1992). Adaption in Nature Artificial Systems, MIT Press. 39. Hopkinson, Neil, Richard Hague, and Philip Dickens, eds. Rapid manufacturing: an industrial revolution for the digital age. John Wiley & Sons, 2006. 40. Howell, L. (2001). Compliant Mechanisms, Wiley. 41. Jensen, B. and L. Howell (2003). "Identification of Compliant Pseudo-Rigid-Body Four-Link Mechanism Configurations Resulting in Bistable Behavior." ASME Journal of Mechanical Design 125: 701-708. 42. Jensen, B. and L. Howell (2004). "Bistable Configurations of Compliant Mechanism Modeled Using Four-Links and Translational Joints." ASME Journal of Mechanical Design 126(4): 657-666. 43. Jernej Barbic, Fun Shing Sin, Daniel Schroeder: Vega FEM Library. 2012. http://www.jernejbarbic.com/vega 149 44. Joo, J. and S. Kota (2004). "Topological Synthesis of Compliant Mechanisms Using Nonlinear Beam Elements." Mechanical Based Design of Structures and Machines 32(1): 17- 38. 45. Joo, J., S. Kota, et al. (2000). "Topological Synthesis of Compliant Mechanisms Using Linear Beam Elements." Mechanical Structures and Machines 28(4): 245-280. 46. Kä stenbauer, T., et al. "Running shoes for relief of plantar pressure in diabetic patients." Diabetic medicine 15.6 (1998): 518-522. 47. Kennedy, J. and R. C. Eberhart (1995). Particle Swarm Optimization. Proceedings of IEEE International Conference on Neutral Networks, Perth, Australia. 48. Kennedy, J., R. C. Eberhart, et al. (2001). Swarm Intelligence. San Francisco, Morgan Kaufman. 49. Kim, M.J., G.W. Jang, and Y.Y. Kim, "Application of a Ground Beam-Joint Topology Optimization Method for Multi-Piece Frame Structure Design, " Journal of Mechanical Design, Vol. 130, DOI: 10.1115/1.2936930, Aug. 2008 50. Kirsch, U., "Optimal Topologies of Truss Structures," Computer Methods in Applied Mechanics and Engineering, vol. 72, pp. 15-28, 1989. 51. Kirsch, U., "Fundamental Properties of Optimal Topologies," Topology Design of Structures, pp. 3-18, Kluwer Academic Publishers, Dordrecht, 1993. 52. Kirsch, U., "Reduction and expansion processes in topology optimization," Topology optimization in structural mechanics, pp. 197-206, Springer, New York, 1997. 53. Leon S. Dimas and Markus J. Buehler, tough and stiff composites with simple building blocks, J. Mater Res., 2013a, 28, 1295. 54. Leon S. Dimas, Graham H. Bratzel, Ido Eylon, and Markus J. Buehler, tough composites inspired by mineralized natural materials, Adv. Funct. Mater. 2013b, 23, 4629-4638. 150 55. Lobontiu, N. (2003). Compliant Mechanisms: Design of Flexure Hinges, CRC Press. 56. Lord M, Hosein R: Pressure redistribution by molded inserts in diabetic footwear: a pilot study. J Rehab Res Develop 31:214–221, 1994. 57. Lu, K. and S. Kota (2006). "Topology and Dimensional Synthesis of Compliant Mechanism Using Discrete Optimization." ASME Journal of Mechanical Design 128(1080-1091). 58. Ma, P. X. and J. Elisseeff (2005). Scaffolding in Tissue Engineering, CRC Press. 59. Martí nez, M., Garcí a-Nieto, S., Sanchis, J., & Blasco, X. (2009). Genetic algorithms optimization for normalized normal constraint method under Pareto construction. Advances in engineering software, 40(4), 260-267. 60. McDowell, D. L. (1998). New Directions in Materials Design Science and Engineering. Atlanta, CA, Georgia Institute of Technology and Morehouse College. 61. Owings, Tammy M., et al. "Custom therapeutic insoles based on both foot shape and plantar pressure measurement provide enhanced pressure relief."Diabetes Care 31.5 (2008): 839-844. 62. Pilkey, W.D., Analysis and design of elastic beams, John Wiley &Sons, Inc., 2002. 63. Rozvany, G.I.N. and W. Prager, "Optimal design of partially discretized grillages," Journal of the Mechanics and Physics of Solids, Vol. 24, pp. 125-136, Jun. 1976. 64. Prager, W., "A note on discretized Michell structures," Computer Methods in Applied Mechanics and Engineering, Vol. 3(3), pp. 349-355, May, 1974. 65. Rozvany, G.I.N. and W. Prager, "Optimal design of partially discretized grillages," Journal of the Mechanics and Physics of Solids, Vol. 24, pp. 125-136, Jun. 1976. 66. Rozvany, G.I.N., "Optimum layout of grillages: allowance for the cost of supports and optimization of support locations," Mechanics of Structures and Machines, Vol. 22, p.49-72, 1994. 151 67. Rozvany, G.I.N., M.P. Bendsø e, and U. Kirsch, "Layout Optimization of Structures," Applied Mechanics Reviews, Vol. 48(2), pp. 41-119, Feb. 1995a. 68. Rozvany, G.I.N., "What is meaningful in topology design? An engineer’s viewpoint," Advances in structural optimization, pp. 149-185, Kluwer Academic Publishers, Dordrecht, 1995b. 69. Rozvany, G.I.N., M.P. Bendsø e, and U. Kirsch, "Layout Optimization of Structures," Applied Mechanics Reviews, Vol. 48(2), pp. 41-119, Feb. 1995a. 70. Sigmund, O. (1997). "On the Design of Compliant Mechanisms Using Topology Optimization." Mechanical Structures and Machines 25: 493-524. 71. Sigmund, O. (2000). "A New Class of Extremal Composites." Journal of Mechanics and Physics of Solids(48): 397-428. 72. Strang, G. and R.V. Kohn, "Hencky-Prandtl nets and constrained Michell trusses," Computer methods in Applied Mechanics and Engineering, Vol. 36, pp. 207-222, 1983. 73. Svanberg, K., "Optimal truss sizing based on explicit Taylor series expansions," Structural and Multidisciplinary Optimization, Vol. 2, pp. 153-162, 1990. 74. Svanberg, K., "Global convergence of the stress ratio method for truss sizing," Structural and Multidisciplinary Optimization, Vol. 8, pp. 60-68, 1994. 75. Taylor, J.E. and M.P. Rossow, "Optimal Truss Design Based on an Algorithm Using Optimality Criteria," International Journal of Solids and Structures, Vol. 13, pp. 913-923, 1977. 76. Taylor, J. E., "Maximum Strength Elastic Structural Design, Proceedings of the ASCE," Journal of the Engineering Mechanics Division, Vol. 95, pp. 653-663, 1969. 77. Topping, B.H.V. "Topology Design of Discrete Structures," Topology Design of Structures, pp. 517-534, Kluwer Academic Publishers, Dordrecht, 1993. 152 78. Wang, H. V. (2005). A Unit Cell Approach for Lightweight Structure and Compliant Mechanism. School of Mechanical Engineering. Atlanta, GA, Georgia Institue of Technology. 79. Dr. Lifeng Wang, Jacky Lau, Nicholas V Soane, Matthew J. Rosario, and Prof. Mary C. Boyce, Mechanical Behavior of Co-Continuous Polymer Composites, Proceedings of the SEM Annual Conference, June 7-10, 2010 Indianapolis, Indiana USA. 80. Lifeng Wang , Jacky Lau , Edwin L. Thomas , and Mary C. Boyce, Co-Continuous Composite Materials for Stiffness, Strength, and Energy Dissipation, Adv. Mater. 2011, 23, 1524–1529. 81. Willis, Karl, et al. "Printed optics: 3d printing of embedded optical elements for interactive devices." Proceedings of the 25th annual ACM symposium on User interface software and technology. ACM, 2012. 82. Wohlers, Terry, and President Tim Caffrey. "Additive Manufacturing: Going Mainstream." MANUFACTURING ENGINEERING 150.6 (2013): 67-73. 83. Xie, Y.M. and G.P. Steven, Evolutionary structural optimization, Springer, New York, 1997.

Abstract (if available)

Abstract Deformable products exist broadly in different industries. These deformable products undergo deformation to achieve specific mechanical properties. Custom made insoles are applied to redistribute plantar pressure under feet to reduce stress related issues, such as diabetic ulceration. Also there are some other deformable products, i.e. cushioned chairs, custom made running shoes, light weighted space structure and so forth. ❧ The basic mechanical property under all these deformable products is stiffness with respect to deformation. Our research interest is to design these deformable geometries with controlled mechanical properties, and in order to do this, we investigated new design approaches and completed achievements to achieve mechanical property control in different research phases. ❧ First, we explored how to design complex geometries to achieve maximal stiffness in a single material. A PSL (principal stress line) based strategy is studied to design an optimal rigid Michell structure like topology. This Michell structure like topology obtained by PSL strategy is more refined than other approaches. We called this strategy reduction and growth approach. ❧ Second, we studied on how to design geometries to achieve target stiffness with multiple materials. Digital material design is investigated to design target stiffness in a given domain with a soft and a stiff material. Our digital material design method for desired stiffness is predictable and can be interpolated into more levels with higher resolution. ❧ Also we investigated how to apply these N‐level digital materials to achieve controlled flexibility skin design. Simulation results for both digital material and flexibility skin design with N‐level digital material have been shown and evaluated. A physical experiment was also investigated to verify the simulation results on digital material. ❧ Third, we investigated how to design geometries to achieve target compliance in a fixed topology of structure. A FCF (flexibility contribution factor) based method is provided to achieve the target compliance design of a geometry. This FCF method increases computation efficiency greatly compared to other optimal approaches. ❧ One challenge is how to fabricate these designed geometries with controlled stiffness,since they are either very complex geometries or contain multiple materials. ❧ Fortunately, 3D printing techniques have the unique capabilities to manufacture complex geometries and multiple materials. Nowadays, custom made products by 3D printing are going to main stream. We are expecting our research achievements to be applied with 3D printing and to impact the manufacturing industry.

Linked assets

University of Southern California Dissertations and Theses

Conceptually similar

PDF

Deformation control for mask image projection based stereolithography process

PDF

Hybrid vat photopolymerization: methods and systems

PDF

Fabrication-aware machine learning for accuracy control in additive manufacturing

PDF

Multi-scale biomimetic structure fabrication based on immersed surface accumulation

PDF

Hybrid vat photopolymerization processes for viscous photocurable and non-photocurable materials

$Metasurfaces in 3D applications: multiscale stereolithography and inverse design of diffractive optical elements for structured light$

PDF

Metasurfaces in 3D applications: multiscale stereolithography and inverse design of diffractive optical elements for structured light

PDF

3D printing of polymeric parts using Selective Separation Shaping (SSS)

PDF

Machine learning-driven deformation prediction and compensation for additive manufacturing

PDF

Energy control and material deposition methods for fast fabrication with high surface quality in additive manufacturing using photo-polymerization

PDF

Correspondence-based analysis of 3D deformable shapes: methods and applications

PDF

Statistical modeling and machine learning for shape accuracy control in additive manufacturing

PDF

Nanostructure interaction modeling and estimation for scalable nanomanufacturing

PDF

Data-driven 3D hair digitization

PDF

Modeling and analysis of nanostructure growth process kinetics and variations for scalable nanomanufacturing

PDF

Feature-preserving simplification and sketch-based creation of 3D models

PDF

3D printing and compression testing of biomimetic structures

PDF

The extension of selective inhibition sintering (SIS) to high temperature alloys

PDF

Accurate 3D model acquisition from imagery data

PDF

Fabrication of ultrasound transducer and 3D-prinitng ultrasonic device

PDF

3D inference and registration with application to retinal and facial image analysis

Asset Metadata

Creator Li, Yongqiang (author)

Core Title Deformable geometry design with controlled mechanical property based on 3D printing

School Viterbi School of Engineering

Degree Doctor of Philosophy

Degree Program Industrial and Systems Engineering

Publication Date 07/15/2014

Defense Date 05/14/2014

Publisher University of Southern California (original), University of Southern California. Libraries (digital)

Tag 3D printing,deformation,digital material,geometry,OAI-PMH Harvest,stiffness,tailored performance meso‐structure,topology optimization

Format application/pdf (imt)

Language English

Contributor Electronically uploaded by the author (provenance)

Advisor Chen, Yong (committee chair), Huang, Qiang (committee member), Jin, Yan (committee member)

Creator Email actionafterthought@gmail.com,yongqial@usc.edu

Permanent Link (DOI) https://doi.org/10.25549/usctheses-c3-440296

Unique identifier UC11287208

Identifier etd-LiYongqian-2685.pdf (filename),usctheses-c3-440296 (legacy record id)

Legacy Identifier etd-LiYongqian-2685.pdf

Dmrecord 440296

Document Type Dissertation

Format application/pdf (imt)

Rights Li, Yongqiang

Type texts

Source University of Southern California (contributing entity), University of Southern California Dissertations and Theses (collection)

Access Conditions The author retains rights to his/her dissertation, thesis or other graduate work according to U.S. copyright law. Electronic access is being provided by the USC Libraries in agreement with the a...

Repository Name University of Southern California Digital Library

Repository Location USC Digital Library, University of Southern California, University Park Campus MC 2810, 3434 South Grand Avenue, 2nd Floor, Los Angeles, California 90089-2810, USA