Investigations of Surface Instabilities in Soft Materials and Phase Transformations in Additive Manufacturing by Qi Zhang B.Eng., Peking University, 2013 a dissertation submitted in partial fulfillment of the requirements for the Degree of Doctor of Philosophy in the School of Engineering at Brown University Providence, Rhode Island may 2019 c Copyright by Qi Zhang, 2019. All rights reserved. This dissertation by Qi Zhang is accepted in its present form by the School of Engineering as satisfying the dissertation requirement for the degree of Doctor of Philosophy Date Professor Janet Blume, Advisor Recommended to the Graduate Council Date Professor Kyung-Suk Kim, Reader Date Doctor Victor Oancea, Reader Approved by the Graduate Council Date Professor Andrew G. Campbell Dean of the Graduate School iii Vita Qi Zhang was born on December 19, 1992 in Weihui, Henan Province, China. He attended Peking University in 2009, and received Bachelor degree in Engineering Structure Analysis in 2013. Following, he entered Solid Mechanics Research Group at Brown University, and received the James R. Rice Graduate Fellowship in Solid Mechanics for the 2013-2014 academic year. His research focused on the theoretical derivations and computational modeling of surface instabilities in compressible soft materials, and simulations of metallurgical phase transformations in additive manu- facturing process. iv Acknowledgements First of all I would like to express my sincere gratitude to Professor Janet Blume for accepting me as her student and offering valuable guidance and tremendous help during my years at Brown. She is a brilliant scholar and also a wonderful mentor. She taught me how to be a capable researcher and showed me a broader picture whenever I got stuck in my own head. Second, I would like to thank Professor Kyung-Suk Kim and Dr. Victor Oancea, who spent a lot of time reading and offering me suggestions for my thesis writing. Discussions with Professor Kim inspired my research work on surface wrinkling and I also benefited a lot from his courses. I would also like to thank and acknowledge all faculty at Brown Solid Mechanics Group. Achieving this degree has been a long and maturing road, and I truly feel privileged that I was able to complete my academic journey at such an amazing group with Professor Bower, Professor Gao, Professor Henann and Professor Kesari. While studying at Brown, I have made numerous friends. They not only provided support and advices for my research, but also offered help and fun in life. To Xin, Teng, Sheng, Kai, Odysseas, Daniel, Bo, Xing, Yue, Zhi and many others: I’ll never forget our scientific discussions and happy hours. Special thanks to Guijin, my five- year roommate, and thanks to the time when we struggled and made achievements. Finally, and with my deepest emotion, I would like to thank my parents, Wenbin Zhang and Guangqin Yan, for their endless support and understanding. One last thank you I would like to give is to Xiaoguai Li, the one and the only in my heart. Without their love and care, I would have never completed this thesis. v To my parents. vi Contents Signature Page . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . iii Vita . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . iv Acknowledgements . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . v List of Tables . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . x List of Figures . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xi 1 Introduction 1 1.1 Surface instabilities in soft materials . . . . . . . . . . . . . . . . . . 1 1.2 Metallurgical phase transformations in additive manufacturing process 4 1.3 Methodology and approaches . . . . . . . . . . . . . . . . . . . . . . 6 1.4 Outline of the thesis . . . . . . . . . . . . . . . . . . . . . . . . . . . 7 2 Wrinkling in Generalized Blatz-Ko Materials I: the Bilayer System and Finite Slabs 10 2.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10 2.2 Generalized Blatz-Ko material . . . . . . . . . . . . . . . . . . . . . . 12 2.2.1 Uniaxial loading test . . . . . . . . . . . . . . . . . . . . . . . 13 2.3 Results and discussion: surface wrinkles in compressible Blatz-Ko ma- terials . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15 2.3.1 f = 1 (compressible Neo-Hookean material) . . . . . . . . . . 15 2.3.2 f =0. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22 vii 2.3.3 0≤f ≤1 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31 2.4 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 32 3 Wrinkling in Generalized Blatz-Ko Materials II: Instability and Imperfection-sensitivity 34 3.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 34 3.2 Stability analysis of wrinkles in finite slabs . . . . . . . . . . . . . . . 36 3.2.1 Linear perturbation analysis . . . . . . . . . . . . . . . . . . . 37 3.2.2 Second-order perturbation analysis for nonlinear stability of a wrinkle . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 40 3.3 Finite element analysis of post-bifurcation wrinkling and creasing . . 44 3.3.1 Critical wrinkling state . . . . . . . . . . . . . . . . . . . . . . 45 3.3.2 Fourier transformation analysis of post-wrinkling deformations 47 3.3.3 Sensitivity to the perturbation: wrinkling and creasing . . . . 49 3.4 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 51 4 A Metallurgical Phase Transformation Framework Applied to Ad- ditive Manufacturing Processes 53 4.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 53 4.2 Finite element simulation of powder-bed fusion additive manufacturing 55 4.2.1 Overall modeling considerations . . . . . . . . . . . . . . . . . 57 4.2.2 AM predictive simulations: spatial and time multiscale consid- erations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 59 4.3 An FE-based generic metallurgical phase transformation Framework . 60 4.3.1 Kinetics of solid-phase transformations . . . . . . . . . . . . . 62 4.3.2 Modeling framework . . . . . . . . . . . . . . . . . . . . . . . 66 4.4 Demonstration of the framework for two different alloys . . . . . . . 68 4.4.1 Ti-6Al-4V . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 69 viii 4.4.2 Steel 5140 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 76 4.5 Mechanical properties prediction for Ti-6Al-4V . . . . . . . . . . . . . 78 4.5.1 An empirically-driven plasticity model for Ti-6Al-4V alloys . . 79 4.5.2 Ramberg-Osgood plasticity model . . . . . . . . . . . . . . . . 80 4.5.3 Microstructure and plastic yield behavior . . . . . . . . . . . . 81 4.5.4 As-printed microstructure . . . . . . . . . . . . . . . . . . . . 82 4.5.5 Microstructure size evolutions during heat treatment . . . . . 84 4.6 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 84 5 Validation and Calibration of The Metallurgical Phase Transforma- tion Framework in Ti-6Al-4V Alloy 86 5.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 86 5.2 Experiments and measurements . . . . . . . . . . . . . . . . . . . . . 88 5.2.1 Column printing with in-situ thermocouple measurements . . 88 5.2.2 Cruciform printing and tensile testing . . . . . . . . . . . . . . 90 5.3 Calibration and validation for Ti-6Al-4V . . . . . . . . . . . . . . . . 97 5.3.1 Column printing and comparison with thermocouple data . . . 97 5.3.2 Cruciform printing model . . . . . . . . . . . . . . . . . . . . 99 5.4 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 106 6 Conclusion 108 Bibliography 111 ix List of Tables 4.1 Latent heat of Ti-6Al-4V. . . . . . . . . . . . . . . . . . . . . . . . . 58 4.2 Other thermal properties of Ti-6Al-4V. . . . . . . . . . . . . . . . . 59 4.3 JMA coefficients for the transformation (α0 → β + αW ). . . . . . . . 71 4.4 Latent heat of Steel 5140. . . . . . . . . . . . . . . . . . . . . . . . . 76 4.5 Thermal properties of Steel 5140. . . . . . . . . . . . . . . . . . . . . 76 4.6 Phase transformations considered in the metallurgical model of Steel 5140. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 76 5.1 Process parameters for the thermal validation exercise involving the thermocouple substrate. . . . . . . . . . . . . . . . . . . . . . . . . . 89 5.2 Volume fraction predictions . . . . . . . . . . . . . . . . . . . . . . . 101 x List of Figures 2.1 Transverse stretch ratio λ2 and volume change as a function of λ1 . . 14 2.2 The first Piola-Kirchoff stress as a function of stretch ratio predicted by the generalized Blatz-Ko material model. . . . . . . . . . . . . . . 14 2.3 Schematic diagram of bilayer system wrinkling. (a) Undeformed state of a soft substrate with width l and thin film with width λps l ; (b) The substrate is pre-stretched with stretch ratio λps ; (c) The thin film layer is firmly bonded to the substrate; (d) The bilayer system is compressed to λ1 l and wrinkle pattern forms on the surface. . . . . . . . . . . . 19 2.4 Relationship between the critical normalized wave number k¯ and λ1 of substrate for different Poisson’s ratio. (c) matches Figure 2a in [92]. 21 2.5 (a, b, c) The first-order perturbation prediction of wavelengths for different Poisson’s ratios. Solid lines show the normalized wavelength in deformed configuration L∗ /h = λ1 L/h = 2πλ1 /kh with respect to compressive strain C L for first-order-predicted post-buckling wrinkles; Dashed lines indicate the critical buckling points. (d, e, f) Effect of pre-stretch on wrinkling behaviors. (c) matches Figure 2b in [92]. . . 23 2.6 Undeformed configuration of finite substrate with width l and thick- ness d. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 26 xi 2.7 Wrinkling eigenmodes in a finit slab with initial coordinates (0 ≤ X1 ≤ 5, 0 ≤ X2 ≤ 1). Solid lines show wrinkle patterns, while dash lines are unperturbed deformations with the same stretch ratio. . . . . . . . . . . . . . . . 29 2.8 Characterization of wrinkle patterns for f = 0. (a) in region I, value of α2 is positive, which corresponds to wrinkling in semi-infinite substrate with a thin film layer; (b) in region II, value of α2 is negative and therefore α is imaginary. This corresponds to wrinkling in a finite slab; (c) in region III, α2 is a complex number and might leads to more complex perturbation forms. . . . . . . . . . . . . . . . . . . . . . . 30 2.9 Characterization of wrinkle patterns for different f . . . . . . . . . . 32 (2) 3.1 λ1 as a function of the stretch ratio λ1 . . . . . . . . . . . . . . . . . 44 3.2 FE model of the finite slab with d=1 and L=5. . . . . . . . . . . . . 45 3.3 Wrinkling deformation of the material with (f = 0&ν = 0.2) around critical point ( = 0.6437). Contour plot demonstrates the deformed maximum principal nominal strain, max p , distribution. (a)(b) prior to the onset of wrinkles ( = 0.4343, 0.5943), the perturbation is first localized near the surface and decays as away from the surface; (c) beyond the critical point ( = 0.7143), the non-decaying wrinkling de- formation, characterized by evenly arranged rectangular cells, develops. 46 3.4 Critical wrinkling strain in materials with f = 0 and 0.5. Solid line is the theoretical predicted critical stretch ratio; markers show strain for uniform deformation to wrinkling state transition determined by FEM simulations, the error bars are due to the displacement-controlled loading steps ∆. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 47 xii 3.5 Post-buckling deformation of Blatz-Ko material at a stretch ratio λ1 = 0.2941. (a) ABAQUS simulation result; (b) Fourier transformation analysis. X and Y axis show frequencies of wrinkling deformations in the X1 &X2 directions, respectively. Solid lines correspond to theoret- ical predictions and dots show FFT calculation of the finite element simulation output;(c) recover wrinkling deformation with the first sev- eral eigenmodes; (d) plot of fy for fx equals the frequency of initial perterbation. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 48 3.6 Imperfection sensitivity test for the material with parameter (f = 0&ν = 0.2). The perturbation is set to be up2 = ζ cos 2πkX  L 1 , where k = 4, 6, 8 and ζ = 0.0002. For all the imperfection cases, finite element simulation results at compressive strain  = 0.7143 are considered. Contour plots of maximum principal nominal strain (max p ) distribution and corresponding deformation frequency analysis results are plotted. . . 50 3.7 Contour plots of wrinkles and creases in a 2D square model. Unde- formed dimension is 1 × 1. (a) ζ = 10−4 ,  = 0.6812; (b) ζ = 10−3 ,  = 0.3959. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 51 3.8 Critical stretch ratio for creasing in the generalized Blatz-Ko material. 51 4.1 Ti-6Al-4V thermal properties: conductivity and specific heat capacity. 58 4.2 Diagram of physical state changes between raw/liquid/solid (RLS) states. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 61 4.3 An example time-temperature-transformation (TTT) diagram. . . . 63 4.4 The algorithm to model a diffusional transformation. . . . . . . . . 64 4.5 The algorithm to model a martensitic transformation. . . . . . . . . 66 4.6 General algorithm to model the generic metallurgical phase transfor- mation framework. . . . . . . . . . . . . . . . . . . . . . . . . . . . 67 4.7 Equilibrium β phase fraction. . . . . . . . . . . . . . . . . . . . . . . 71 xiii 4.8 TTT-diagram for diffusional phase transformations β → αgb and β → αW in Ti-6Al-4V. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 72 4.9 Ti-6Al-4V metallurgical phase transformations and α-lathe width pre- diction one element analysis 1. . . . . . . . . . . . . . . . . . . . . . 73 4.10 Ti-6Al-4V metallurgical phase transformation one element analysis 2. 74 4.11 The solidification map for β grain morphology prediction. . . . . . . 75 4.12 Steel 5140 metallurgical phase transformations TTT-diagram. . . . 77 4.13 Steel metallurgical phase transformations. . . . . . . . . . . . . . . . 78 4.14 Diagram illustrating the approach to the empirically-driven plasticity model for Ti-6Al-4V. . . . . . . . . . . . . . . . . . . . . . . . . . . 80 5.1 Image of the thermocouple substrate (left) and an illustration (right) showing the location of the blind holes for thermocouple placement. 89 5.2 Diagram (left) showing the naming convention for the cruciform spec- imens and a photograph (right) of the as-built cruciform specimens. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 90 5.3 Heat treatment of CRUC2. . . . . . . . . . . . . . . . . . . . . . . 91 5.4 Photographs of tensile specimens. . . . . . . . . . . . . . . . . . . . 92 5.5 Stress-strain curves from the cruciform tensile test specimens. . . . 93 5.6 Grain orientation maps from heat treated and as-deposited specimens. 94 5.7 Back-scatter SEM images of heat treated and as-deposited specimens. 94 5.8 Backscatter SEM images and elemental maps for the heat treated spec- imen: Backscatter image (left); V chemical map (right). . . . . . . . 95 5.9 XRD patterns obtained in Bragg-Brentano configuration for as de- posited and heat treated specimens. Black: as-deposited, Red: heat treated. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 96 5.10 FE model of column printing experiments. . . . . . . . . . . . . . . 98 xiv 5.11 Column printing, thermo-couples temperature experiment measure- ments vs. FEA predictions at cube 5. top) first 20 layers printing, the simulation data was taken from the TC 0.0625mm below the cube; bottom) printing of the second layer. Blue solid line shows the ther- mocouple recorded temperatures, and all other lines represent finite element predictions at different depths beneath the cube, as shown in the table. The closer the position is to the cube, the higher the temperature is predicted. . . . . . . . . . . . . . . . . . . . . . . . . 99 5.12 Geometry (left) and mesh (right) of the SLM process on one dog bone specimen in the cruciform model . . . . . . . . . . . . . . . . . . . . 100 5.13 Temperature evolutions: snapshot during a lasering sequence. Gray contourplots represent temperature as above melting. . . . . . . . . 101 5.14 Histogram of measured αlath thickness from image processing of EBSD data in several samples. . . . . . . . . . . . . . . . . . . . . . . . . . 102 5.15 FE model of the rectangular sample for virtual tensile tests. . . . . . 104 5.16 Stress-strain curves of the rectangular samples under uniaxial tension. 105 5.17 Heat transfer analysis of dogbone specimen. a) FEA model; b) lath thickness contour plots of both as-built sample and heat treated sam- ple. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 105 5.18 Tensile tests of dogbone specimens. a) necking behavior captured in FE simulations for both as-built (upper plot) and heat-treated (lower plot) samples; b) comparison of stress-strain curves between finite ele- ment predictions with experimental results. Properties of both as-built sample and heat treated sample are included. . . . . . . . . . . . . . 106 xv Chapter 1 Introduction 1.1 Surface instabilities in soft materials Soft materials, such as rubbers, elastomers, foams, and living tissues are often com- pressed by mechanical loading, swelling, or growth. When the compressive strain reaches a critical point, instability forms, e.g., wrinkles or creases, are found on the surface of soft materials [35, 22, 68, 110]. Surface instabilities are widely observed in nature, such as undulations on the skin, folding of brains, and morphogenesis of plants [64, 85, 7, 15, 69]. Among those patterns, wrinkles and creases are two widely observed types of instability. Wrinkles are characterized by a field of infinitesimal strain superimposed upon a large strain in a finite region, whereas a crease deviates from the flat state by a field of large strain in an infinitesimal region [48]. Surface in- stabilities were originally seen as phenomena to be avoided, but more recently, various surface patterns have been seen engineering applications, such as stretchable electron- ics, measuring the mechanical properties of materials, tune adhesion, and fabricating surfaces [12, 52, 24, 35, 70, 21]. It is important to understand how those patterns are formed and what is the critical condition for them to develop. Investigations of the post-wrinkling process can lead to ways to control the surface patterns. As 1 soft materials are compressed beyond the critical point, those instabilities develop and transit to advanced Ruga phases, such as, multi-mode wrinkles, folds, and ridges [28, 112, 14, 42, 54, 81, 92]. Much theoretical and experimental effort has been put forward to characterize the critical conditions for wrinkling, wrinkle patterns as well as post-buckling process [12, 52, 24, 35]. In one early study, Biot derived a critical compressive strain 0.456 for wrinkling in a semi-infinite incompressible solid under plain-strain compression, using a bifurcation analysis [10]. A more complete analysis of elasticity for buckling was provided by Shield et al. [89], which presented the exact formulation of the first-order perturbation analysis for the wrinkle formation with incrementally small- strain deformation of an isotropic substrate. However, it has been shown by Cao et al. [18] that wrinkling in semi-infinite homogeneous incompressible neo-Hookean solids is highly unstable and imperfection-sensitive, while creases are stable and could occur at a strain lower than the Biot critical strain. The critical strain for creasing in incompressible neo-Hookean material was measured to be 0.35 ± 0.07 by bending elastomer surface experimentally [34], and it was further verified numerically to be 0.35 [48]. Wrinkles are stabilized when a stiff thin-film is firmly attached to the substrate, and the wrinkling in the bilayer system has been studied extensively [19, 112, 28]. It has been shown that the critical compressive strain is influenced by the pre-stretch on the substrate, the film/substrate modulus ratio, and the variation in substrate properties [92, 19, 28]. A stiffer and thinner film tends to form smooth wrinkles, while a softer and thicker film tends to form creases [5]. For bilayer systems with modest stiffness ratio between the film and the substrate, the bifurcation strain is relatively large. For stiffer films, the bifurcation strain is relatively small and it is also applicable to any linear elastic film material [19]. A computational study was performed by Tallinen et al. [94] to investigate how the stiffness difference between 2 two layers influences the transition from wrinkles to creases or folds. Kim et al. [28] showed that wrinkles can transition to creases for graded samples with an elastic stiffness that decreases with distance from the free surface. Later, they studied the evolution of the surface instabilities in bilayers and constructed a phase diagram of three qualitatively different types of post-wrinkling bifurcations [112]: (i) wrinkles that transit directly to creases, (ii) wrinkling-period doubling-creasing transition, and (iii) wrinkles-period doubling-period quadrupling-folding transition. It was shown that as the pre-stretch of the substrate is increased, wrinkles are destabilized and period doubling occurs earlier at a lower strain [4], and if the substrate is subject to a sufficiently large pre-stretch, ridge mode occurs through a different ruga phase path [17]. To date, the theoretical and computational efforts in analyzing the surface instabili- ties are mostly restricted to the film-substrate systems of incompressible hyperelastic materials, e.g. neo-Hookean model. This thesis seeks to explore the effects of sub- strate compressibility on the formation of surface wrinkles and shows a new wrinkling type that is only admissible in a finite slab (rather than in half-spaces). The mate- rial is modeled by a generalized Blatz-Ko constitutive relation. This material model (introduced in Chapter 2) is widely used for compressible isotropic nonlinearly elas- tic solids. It has 3 adjustable parameters, which may be specified to reduce the form to represent commonly used incompressible and compressible materials such as neo-Hookean, Mooney-Rivlin and Blatz-Ko models [44]. For instance, the Blatz-Ko material, which has been arrived at by experiment [11] and whose deformations have been studied previously [8, 44], is obtained from the considered class of materials by specifying two parameters involved in the definition of the class. This generalized model provided more flexibility and wider choices of material models. Therefore, it is adopted in this thesis to investigate the surface instabilities in compressible soft materials. 3 1.2 Metallurgical phase transformations in addi- tive manufacturing process Additive manufacturing (AM) is the process of joining materials to make objects from 3D model data, usually layer upon layer, as opposed to subtractive manufacturing methodologies, such as traditional machining [78]. The development of additive man- ufacturing started in the 1980s [62], and significant progress has been made since then [46, 9]. It opens up the possibility of producing parts that are characterized by complex geometries and are thus difficult to build using conventional manufacturing techniques [102]. Although the AM process is still hampered by low productivity, poor quality and uncertainty of final part mechanical property [9], it is reasonable to expect that it will play an increasingly important role in manufacturing. First of all, it enables engineers to make efficient use of materials and obtain products with minimal waste [66, 78, 61]. Since there is much less manufacturing constraint, e.g. less requirement for post-processing, for AM than conventional processes, AM can be applied to increase design freedom and achieve topology optimizations [13, 109]. In addition, AM has the ability to produce a consolidated part that was previously generated by many separated parts, and thus reduce energy and natural resource costs [23]. Therefore, it is also environmentally friendly. Due to the previous advan- tages, AM has seen many applications in aerospace, automotive and biomedical area [39, 97, 91, 36]. Metal powder bed fusion (PBF) based AM processes are capable of fabricating metal- lic parts directly from many engineering alloy powders, such as stainless steel, tita- nium alloys, etc. [38]. Selective laser melting (SLM) and electron beam melting (EBM) are the most commonly used PBF based AM processes [63, 96, 76]. During the printing process, metallic AM parts/alloys undergoes a complex thermal history and complicated metallurgical phase transformations which introduce complexities to 4 the analysis of microstructural evolution and properties not typically found in con- ventional processes [32]. While significant progress has been made in the last few years, the reliability of AM parts is often less than desirable as they suffer from man- ufacturing defects and hence subpar strength and fatigue life [106]. To address this challenge, there is an increasing interest towards simulation of manufacturing pro- cesses to provide insight into the process and help accelerate progress in raising the quality of AM parts, including predicting thermal evolutions, part distortions, and residual stresses [82, 113, 77, 84]. Material modelling is a crucial part in these mod- els and is particularly challenging for the parts which experience repeated melting, solidification, and metallurgical phase transformations [75]. The commercial finite element software ABAQUS provides a general simulation framework that is capable of performing thermal and mechanical analysis during printing processes which ac- counts for material deposition and moving energy sources [1]. In this thesis, a generic metallurgical phase transformation framework applicable to arbitrary metal alloys is introduced. The user systematically defines all possible transformations that can take place via a parent-children paradigm (e.g., austenite to martensite); the temperature conditions when transformations can occur with associated transformation diagrams; and whether the transformations are reversible, diffusional, or non-diffusional. The framework calibrates the parameters in evolution equations, such as the Johnson- Mehl-Avrami-Kolmogorov equation [105, 6, 60] or the Koistinen-Marburger equation [58], from input and predicts temperature evolutions, phase transformations during printing and heat treatment, resulting microstructure, and mechanical properties. Titanium alloys have been key metals for the aircraft industry because of the com- bination of excellent mechanical properties with light weight [67]. Additively man- ufactured Ti-6Al-4V parts are characterized by excellent fracture resistance, fatigue behavior and corrosion resistance [38, 90]. These properties are very sensitive to the microstructure of the material, and therefore understanding and control of the mi- 5 crostructure evolution is of great importance, as it is in all titanium manufacturing processes [104, 98, 25]. In this study, a comprehensive Ti-6Al-4V metallurgical model [74] is adopted and implemented in the generic framework mentioned above. The model accounts for complicated phase transformations involved in the printing pro- cess and heat treatment process, and predicts microstructure of the processed parts. In addition, this work adopts existing empirical plasticity models [27, 40, 72, 49] for Ti-6Al-4V alloys that connect elastic behavior, yield strength/strain, and ultimate strength/ductility limits to phase volume fraction content and grain size information. The models utilize the predictions described above to estimate, in the context of a Ramberg-Osgood plasticity model, the hardening behavior of the as-printed and heat treated alloys. 1.3 Methodology and approaches The studies in this thesis involve both theoretical derivations and numerical analy- sis. The simulations are mainly carried out by finite element analysis (FEA) using commercial software package ABAQUS/Standard [41]. For simulations of surface instability in soft materials, ABAQUS/standard provides an automatic mechanism to stabilize the unstable quasi-static problem through the automatic addition of volume-proportional damping to the model [41]. In the FEA process, a viscous force of the form Fv = cM ∗ v is added to the global equilibrium equations P − I − Fv = 0, 6 where M ∗ is an artificial mass matrix calculated with unity density, c is a damping factor, v = ∆u/∆t is the vector of nodal velocities, and ∆t is the increment of time. New physics-based FEM formulations and models implemented in ABAQUS were used to simulate the 3D printing process [1]. Accurate physics-based computer sim- ulations can provide significant insight into the process and can help guide the opti- mization of process parameters and increase the quality of additively manufactured parts. The ABAQUS FE solver offers functionalities, including machine information (e.g., powder recoating sequence, laser scan path, material deposition of the printing head, etc.) and progressive element activation, to provide accurate, scalable, and predictive simulation of part-level 3D printing processes. 1.4 Outline of the thesis In this thesis, I will first investigate a class of surface instabilities in the generalized Blatz-Ko materials with theoretical derivations and finite element simulations. Then, I will develop a generic user subroutine for AM metallurgical phase transformations in ABAQUS and connect the microstructure and mechanical properties of additively manufactured Ti-6Al-4V parts. This thesis consists of six chapters and is organized as follows. Chapter 1 (this chapter) is an overview introduction to the surface instabilities of soft materials that have been observed in nature and recent studies, and of the finite element modelling of material evolution in AM process. The challenges and unsolved problems in those fields are discussed. The methodology of finite element simulation is introduced. The research motivation and goal are also presented. In Chapter 2, the effect of substrate compressibility on surface wrinkling in soft materials is analyzed theoretically. The generalized BlatzKo model, which may re- duce to commonly used incompressible and compressible materials is employed. The 7 wrinkling behavior of a thin film bonded to a soft compressible substrate is also an- alyzed. Subsequently, by performing a linear perturbation analysis, we also find a bi-sinusoidal wrinkling pattern in slabs of finite thickness, which is admissible in a subclass of the materials considered. In this chapter, I focus on the characteriza- tion of critical wrinkling pattern determined by material compressibility and other properties. In Chapter 3, the stability of the bi-sinusoidal wrinkles in a finite slab are studied. I first investigated the stability of the wrinkles using Koiters method. Then, extensive finite element simulations were performed to validate the theoretically determined critical condition of this wrinkling deformation. The post-wrinkling deformation was also analyzed and compared with theoretical predictions though Fourier transforma- tion analysis. In Chapter 4, a metallurgical phase transformation modeling framework for additive manufacturing process is presented. The framework accounts for the phase transfor- mations (physical state changes) involving raw/liquid/solid material states followed by metallurgical solid state phase transformations induced by either rapid heating or cooling events associated with typical 3D printing sequences but also with slower- rate temperature evolutions as associated with heat treatment applications. Two detailed exemplifications are also given in this chapter. Subsequently, an empirical plasticity model for Ti-6Al-4V alloy is utilized to connect microstructure predictions with mechanical properties. In this chapter, I focus on the development and general application of the phase transformation framework. In Chapter 5, the metallurgical phase transformation framework is employed to model the AM processes of Ti-6Al-4V alloy. By modeling the 3D printing process, heat- treated process and tensile tests, the predictive ability of the numerical models is assessed. Experimental models are introduced and numerical simulation predictions are compared with the temperature measurements, EBSD/XRD microstructural ex- 8 aminations, and tensile tests on dogbone specimens. In this chapter, I focus on the validation of the generic framework and calibration of the empirical mechanical mod- els. Chapter 6 concludes the entire thesis with major scientific findings and discussions. 9 Chapter 2 Wrinkling in Generalized Blatz-Ko Materials I: the Bilayer System and Finite Slabs 2.1 Introduction Wrinkling in rubber-like materials has been studied extensively. In one early study, Biot derived a critical compressive strain 0.456 for wrinkling in a semi-infinite incom- pressible solid under plain-strain compression, using a bifurcation analysis [10]. Sur- face instabilities were originally seen as phenomena to be avoided, but more recently, wrinkling in thin-film-substrate systems has been seen to have various applications, such as stretchable electronics, measuring the mechanical properties of materials and fabricating surfaces [12, 52, 24, 35]. Much theoretical and experimental effort has been put forward to characterize the critical conditions for wrinkling, wrinkle pat- terns as well as post-buckling process [45, 18, 108]. It has been shown that the critical compressive strain is influenced by the pre-stretch on the substrate, the film/substrate modulus ratio, and the variation in substrate properties [92, 19, 28]. Wrinkling in 10 an isotropic half-space under biaxial plane stress has been studied for compressible elastic materials [101]. This study seeks to explore the effects of substrate compressibility on the formation of surface wrinkles and introduces a new wrinkling type that is only admissible in a finite slab. The substrate material is modeled by a generalized Blatz-Ko constitutive rela- tion. This material model (Equation (2.3)) is widely used for compressible isotropic nonlinearly elastic solids. It has 3 adjustable parameters, which may be specified to reduce the form to represent commonly used incompressible and compressible mate- rials such as neo-Hookean, Mooney-Rivlin and Blatz-Ko models [44](see Section 2.2). Section 2.3 of this paper provide a detailed analysis of wrinkling of the generalized Blatz-Ko material. Wrinkling of compressible neo-Hookean materials (found by setting the parameter f in Equation (2.3) equal to 1) is analyzed in Section 2.3.1, which matches results in [92] for the incompressible case. Effects of the substrate pre-stretch and Poisson’s ratio (0 < ν < 0.5) are discussed. In Section 2.3.2 the wrinkling state for the material with the parameter f = 0 with varying compressibility is discussed. We derive conditions for wrinkling in a finite slab in this case. The theoretical analyses reveal that the compressibility as well as the parameter f (related to the void ratio of the sponge material) play important roles in the wrinkles and will be summarized in Section 2.4. It has been shown that these wrinkles are unstable in a half-space homogeneous incompressible neo-Hookean solid [18] and creases are more prone to develop [43].This chapter does not address the stability of the wrinkles in finite slabs. However, it will be discussed in Chapter 3. 11 2.2 Generalized Blatz-Ko material The strain-energy density per unit undeformed volume for a homogeneous, isotropic elastic compressible material has the form W = W (I1 , I2 , I3 ) , (2.1) where I1 , I2 , I3 are the principal invariants of the deformation tensor. B = FFT where F is the deformation gradient tensor: I1 = trB = λ21 + λ22 + λ23 , 1 2 I1 − tr B2 = λ21 λ22 + λ22 λ23 + λ23 λ21 ,  I2 = (2.2) 2 I3 = detB = λ21 λ22 λ23 , where λi (i = 1, 2, 3) are the principal stretches. Now we consider a specific form of Equation (2.1), known as the generalized Blatz-Ko material, proposed by Blatz and Ko [11] to describe a foam-rubber material:   µ 1 (1 − 2ν) −ν/(1−2ν) W = f I1 − 1 − + I3 2 ν ν   (2.3) µ I2 1 (1 − 2ν) ν/(1−2ν) + (1 − f ) −1− + I3 , 2 I3 ν ν where the constants µ, ν are the shear modulus and Poisson’s ratio at infinitesimal deformations and 0 ≤ f ≤ 1 is a constant related to the volume fraction of voids in the foam-rubber material. 12 Two specific cases of (2.3) have been extensively studied. When f = 0, (2.3) reduces to   µ I2 1 (1 − 2ν) ν/(1−2ν) W = −1− + I3 (2.4) 2 I3 ν ν which characterizes foamed polyurethane elastomers. When f = 1, it becomes   µ 1 (1 − 2ν) −ν/(1−2ν) W = I1 − 1 − + I3 (2.5) 2 ν ν which models solid polyurethane rubbers.[8] As ν goes to 1/2, Equation (2.3) yields µ µ W = f (I1 − 3) + (1 − f ) (I2 − 3) (2.6) 2 2 Which is known as Mooney-Rivlin model for incompressible materials. It’s also shown that the particular foam rubber material in [11] can be characterized by a specific form of (2.4) when Poisson’s ratio ν = 1/4,   µ I2 1/2 W = + 2I3 − 5 (2.7) 2 I3 which is often called the Blatz-Ko material. Similarly, (2.5) may be viewed as a com- pressible generalization of the neo-Hookean incompressible material W = µ (I1 − 3) /2 [44]. These two models (2.5) & (2.7) will be further discussed in Sections 2.3.1 and 2.3.2 to investigate how substrate compressibility influences wrinkling behavior of soft materials. 2.2.1 Uniaxial loading test To illustrate the properties of these materials, we show their behavior in uniaxial tension. With the longitudinal stretch λ1 , material (2.4) and (2.5) have similar defor- mation response but different stresses. Figure 2.1 shows the transverse stretch λ2 and 13 V as a function of λ1 for the generalized Blatz-Ko material, independent of parameter f . However, Poisson’s ratio & f has an effect on the first Piola-Kirchoff stress curve. (a) λ1 & λ2 (b) λ1 & Volume ratio Figure 2.1: Transverse stretch ratio λ2 and volume change as a function of λ1 . Figure 2.2: The first Piola-Kirchoff stress as a function of stretch ratio predicted by the generalized Blatz-Ko material model. 14 2.3 Results and discussion: surface wrinkles in compressible Blatz-Ko materials 2.3.1 f = 1 (compressible Neo-Hookean material) Wrinkling of the semi-infinite bilayer system as shown in Figure 2.3 is analyzed in this section. The substrate is represented as a compressible neo-Hookean material and the thin film is an incompressible neo-Hookean solid. We show the eigenmodes of the deformation in the soft substrate and how compressibility and pre-stretch of the substrate affect the onset of the wrinkles as well as the evolution of the eigenmodes with further compression. 2.3.1.1 Analysis of the substrate deformation The elastic potential is given by Equation (2.5), and thus the first Piola-Kirchhoff stress can be derived as, ∂W  − ν  T= = µ F − I3 1−2ν F−T . (2.8) ∂F We consider a half space compressible substrate in (−∞ < X1 < ∞, −∞ < X2 < 0), ∂x1 which undergoes a homogeneous plane-strain deformation with stretches λ1 = ∂X1 , ∂x2 λ2 = ∂X2 , λ3 = 1. xi and Xi represent the deformed and undeformed coordinates respectively. Then the deformation gradient is given by    λ1    F0 =   λ2 .  (2.9)   1 15 Combining Equation (2.8) & (2.9) gives the first Piola-Kirchhoff stress components n o 0 0 2 2 2 ν/(2ν−1)−1 as T12 = 0 and T22 = µ λ2 − λ1 λ2 (λ1 λ2 ) . Then the unperturbed solution −ν/(1−ν) 0 0 is solved as λ2 = λ1 from the boundary condition T12 = T22 = 0. Now we consider applying an infinitesimal perturbation of the form:  x1 = λ1 X1 + δg1 (X1 , X2 )       x2 = λ2 X2 + δg2 (X1 , X2 ) (2.10)      x3 = X 3  where δ  1 is a perturbation parameter and g1 and g2 are arbitrary functions to be ∂gi determined. With gi,j = ∂Xj , the deformation gradient is    λ1 + δg1,1 δg1,2 0    F=  δg 2,1 λ2 + δg 2,2 0 ,  (2.11)   0 0 1 and Equation (2.8) gives   2/(ν−1) 1/(ν−1) T11   λ1 + λ21 (1 − 2ν) g1,1 + 2νλ1 g2,2 (1+ν)/(−1+ν) + O δ2 ,  = λ1 − λ1 +δ µ 1 − 2ν T12  1/(ν−1)  g2,1 + O δ 2 ,  = δ g1,2 + λ1 µ T21  1/(ν−1)  g1,2 + g2,1 + O δ 2 ,  = δ λ1 µ   1/(ν−1) T22 2 νλ 1 g 1,1 + (1 − ν) g2,2 + O δ2 .  =δ µ 1 − 2ν (2.12) 16 With the above results, the equilibrium equations are given, up to the first order, as   2/(ν−1) 1/(ν−1) T11,1 + T12,2 1 − 2ν + λ1 g1,11 + (1 − 2ν) g1,22 + λ1 g2,12 + O δ 2 = 0,  =δ µ 1 − 2ν 1/(ν−1) T21,1 + T22,2 λ1 g1,12 + (1 − 2ν) g2,11 + 2 (1 − ν) g2,22 + O δ2 = 0  =δ µ 1 − 2ν (2.13) Following Sun et al. [92], we assume that g1 and g2 have the form of   g1 = C1 sin (kX1 ) exp (αkX2 )  (2.14)  g2 = C2 cos (kX1 ) exp (αkX2 )  where α > 0 and (C1 , C2 ) are coefficients to be determined. k = 2π/L is the wrinkle wave number with L the undeformed configuration wavelength of the perturbation. We then plug it back into equilibrium equations (2.13) and solve for α. We find two r  2/(−1+ν) eigenvalues, α = 1 and α = 1 − 2ν + λ1 / (2 − 2ν), which reduce to α = 1 and α = 1/λ21 (results from [92]) if ν = 1/2. The corresponding two eigenmodes are r  1/(−1+ν) 1/(1−ν) 2/(−1+ν) C2 = −λ1 C1 and C2 = −λ1 1 − 2ν + λ1 / (2 − 2ν)C1 . Combining these two eigensolutions gives the first Piola-Kirchhoff stress components as 17 T11    n   (1+ν)/(−1+ν) 1/(−1+ν) 1/(1−ν)   = λ1 − λ1 − kδ A λ1 + λ1 exp (kX2 )     µ  s  s  2/(−1+ν) 2/(−1+ν)  1 − 2ν + λ1 1 − 2ν + λ1     cos kX1 + O δ 2  +2B exp kX2        2 − 2ν 2 − 2ν    T12   n   1/(1−ν) 1/(−1+ν) = −kδ A λ1 + λ1 exp (kX2 )        µ   s  2/(−1+ν)  1 − 2ν + λ1    1/(−1+ν) sin kX1 + O δ 2  +2Bλ1 exp kX2          2 − 2ν     T21   = −kδ {2A exp (kX2 )  µ   s  2/(−1+ν)  1 − 2ν + λ1     2/(−1+ν)  sin kX1 + O δ 2  +B 1 + λ1 exp kX2         2 − 2ν    n   o  2/(−1+ν) kδ 2A 1 − 2ν + λ1     T22 exp (kX2 ) cos kX1 =    µ 2/(−1+ν) 1 − 2ν + λ1       q  q    2/(−1+ν) 2/(−1+ν)  2/(−1+ν) 1−2ν+λ1 1−2ν+λ1     kδ 2B (−1 + ν) 1 + λ1 2−2ν exp kX2 2−2ν cos kX1 −   2/(−1+ν)  1 − 2ν + λ1       + O δ2    (2.15) where A and B are coefficients to be determined by boundary conditions. The dis- placement components at the interface between substrate and stiff layer are obtained as  s ! 1/(1−ν) 1/(−1+ν) 2 − 2ν kx1 + O δ2    u1 (x1 , 0) = −δ Aλ1 + Bλ1 sin   2/(−1+ν) λ1 1 − 2ν + λ1  u2 (x1 , 0) = δ (A + B) cos kx1 + O δ 2 .     λ1 (2.16) 18 The Cauchy stress components at the interface can also be derived as  σ12 ν/(1−ν) n  2/(−1+ν) o kx1 2  = −kδλ 2A + B 1 + λ sin + O δ   1 1 µ λ1        r   √   2/(−1+ν) 2/(−1+ν)  kδ 2A (1 − ν) 1 − 2ν + λ1 + 2B (1 − ν) 1 + λ1 cos kx 1    σ22 λ1 = r  µ  2/(−1+ν)  λ (1 − ν) 1 − 2ν + λ     1 1      + O δ2 .    (2.17) 2.3.1.2 Plate-theory analysis of bilayer system Following [92], an undeformed film is bonded to the substrate with pre-stretch λps and the bilayer system subsequently undergoes plane strain compression (Figure 2.3). The final stretch for the substrate is denoted by λ1 and thus the compressive strain of thin film is C L = (λps − λ1 ) /λps . Figure 2.3: Schematic diagram of bilayer system wrinkling. (a) Undeformed state of a soft substrate with width l and thin film with width λps l ; (b) The substrate is pre-stretched with stretch ratio λps ; (c) The thin film layer is firmly bonded to the substrate; (d) The bilayer system is compressed to λ1 l and wrinkle pattern forms on the surface. 19 For this film-substrate system, the stress components σ12 & σ22 and displacements u1 & u2 have to be continuous at the interface. Therefore, expressions in (2.16) and (2.17) also describe the stress and displacement in the film at the interface. The relationship between the stresses and displacements from plate-theory analysis [89] is given as,  2 4 h ∂ u2 (x1 , 0) h ∂ 3 u1 (x1 , 0) 2  EL h C ∂ u2 (x1 , 0)   − + L + σ22 (x1 , 0) = 0 (1 − νL2 ) 3 ∂x41 ∂x31 ∂x21   2  2 ∂ u1 (x1 , 0) h ∂ 3 u2 (x1 , 0)   EL h   − − σ21 (x1 , 0) = 0 (1 − νL2 ) ∂x21 2 ∂x31 (2.18) where C L is the thickness average of 11 in the thin film layer, and (EL , νL ) are Young’s Modulus and Poisson’s ratio of the layer. Substituting Equations (2.16) and (2.17) into the above equations gives     M11 M12   A     = 0. (2.19) M21 M22 B where  s    1 − 2ν + λ2/(−1+ν) (2−ν)/(1−ν) ¯ 2 M11 = 2k¯3 − 3λ1 k − 6λ21 C ¯ 3 1  L k + 12λ1 µ ¯   2 − 2ν       s s  2/(−1+ν) 2/(−1+ν)   ¯ 3 1 − 2ν + λ1 ν/(−1+ν) ¯ 2 2 C¯ 1 − 2ν + λ1 M12 = 2k − 3λ1 k − 6λ1 L k   2 − 2ν 2 − 2ν         2/(−1+ν) + 6λ31 1 + λ1      µ ¯ s   1 − 2ν + λ2/(−1+ν) (1+ν)/(−1+ν) ¯ 2 (−1+2ν)/(−1+ν) (−2+3ν)/(−1+ν)  1 M21 = λ1 k − 2λ1 − 4λ1 µ ¯   2 − 2ν       s 2/(−1+ν)  1 − 2ν + λ1   (1+ν)/(−1+ν) ¯ 2 (1+2ν)/(−1+ν) ¯ M22 = λ1 k − 2λ1 k        2 − 2ν   s 2/(−1+ν) 1 − 2ν + λ1      (−2+3ν)/(−1+ν) 3ν/(−1+ν)   − 2 λ1 + λ1 µ ¯ 2 − 2ν (2.20) 20 with k¯ = kh and µ ¯ = µ (1 − νL2 ) /EL . Then det M = 0 gives,  q  2/(−1+ν) q ( 2 λ2/(−1+ν) − 1−2ν+λ1 k¯4  2/(−1+ν) 1−2ν+λ1 2/(−1+ν) 1 2−2ν 8 −1 + λ1 2−2ν ¯k¯3 µ C L = + λ21 λ1   s  s  2/(−1+ν) 2/(−1+ν) 1/(1−ν) λ4/(−1+ν) 2/(−1+ν)  1 − 2ν + λ1 1 − 2ν + λ1 + 12λ1 1 + 3λ1 1− −  ¯k¯2 µ  2 − 2ν 2 − 2ν  s  ),( 2/(−1+ν)  2/(−1+ν)    2/(−1+ν)  1 − 2ν + λ1 + 24λ1 −1 + λ1 ¯k¯ + 48λ21 1 + λ1 µ −1 + µ ¯2 24k¯ 2 − 2ν  s  s ) 2/(−1+ν) 2/(−1+ν) λ2/(−1+ν) 1 − 2ν + λ1   k¯ + −λ1 + λ(1+ν)/(−1+ν)  1 − 2ν + λ1 1 − 1 µ ¯ . 2 − 2ν 2 − 2ν (2.21) This form of C L indicates that it is always positive and has a minimum in the region of k¯ > 0, λ1 > 0, which represents the critical buckling point for a particular bilayer system with normalized shear modulus µ ¯ and substrate Poisson’s ratio ν. Figure 2.4 shows plots of the substrate stretch ratio λ1 and buckling wave number k¯ at which the layer critically buckles for different µ ¯ and Poisson’s ratio of the substrate. (a) ν = 0 (b) ν = 0.25 (c) ν = 0.5 Figure 2.4: Relationship between the critical normalized wave number k¯ and λ1 of substrate for different Poisson’s ratio. (c) matches Figure 2a in [92]. With C L = (λps − λ1 ) /λps , Equation (2.21) gives the first-order perturbation pre- dicted wave number k¯ as a function of nominal compressive strain for given Poisson’s ratio ν, pre-stretch λps and shear modulus ratio µ ¯. Solid line in Figure 2.5c shows 21 ¯ = 1 × 10−5 and λps = 1.5 where the substrate the roots of Equation (2.21) with µ Poisson’s ratio is taken to be 0.5. We plot the wavelength in the deformed configu- ration (L∗ = λ1 L = 2πλ1 /k) as a function of the nominal compressive strain C L near the buckling point. As shown in the figure, there is a bifurcation at the critical buckling point for the bilayer system. As the nominal compressive strain increases, the wavelength of the primary buckling mode decreases and the wavelength of the secondary buckling mode increases rapidly. As pointed out in [92], the bifurcation phenomenon was verified with computational simulation in [30] & [47] and observed in experiments in [89]. As shown in the first three subfigures of Figure 2.5, increased substrate compressibility leads to the decrease of critical buckling strain. However, the wavelength at the critical buckling point almost remains the same. Note that the linear perturbation analysis has limited accuracy beyond the buckling point. A higher order perturbation analysis is needed to analyze the post-buckling behavior. The role of the substrate pre-stretch is also investigated. Pre-elongation of the sub- strate influences the wrinkles in a similar way as does increased incompressibility of the substrate material, as shown in the last three subfigures of Figure 2.5. It in- creases the critical buckling strain and decreases the wavelength of the secondary buckling mode, which is similar to the effect of replacing the substrate with a more incompressible material. 2.3.2 f =0 We then study the wrinkling in another subclass of the generalized Blatz-Ko material where the parameter f = 0. The strain-energy potential is given by Equation (2.4), and the first Piola-Kirchhoff stress can be derived as, I1 F − BF − I2 F−T   ∂W ν/(1−2ν) −T T= =µ + I3 F . (2.22) ∂F I3 22 (a) ν = 0 (b) ν = 0.25 (c) ν = 0.5 (d) ν = 0 (e) ν = 0.25 (f) ν = 0.5 Figure 2.5: (a, b, c) The first-order perturbation prediction of wavelengths for differ- ent Poisson’s ratios. Solid lines show the normalized wavelength in deformed configu- ration L∗ /h = λ1 L/h = 2πλ1 /kh with respect to compressive strain CL for first-order- predicted post-buckling wrinkles; Dashed lines indicate the critical buckling points. (d, e, f) Effect of pre-stretch on wrinkling behaviors. (c) matches Figure 2b in [92]. 0 0 Following the same procedure as in Section 2.3.1, boundary condition T12 = T22 =0 −ν/(1−ν) gives λ2 = λ1 . we assume the deformed configuration of the half substrate to be, 23  x1 = λ1 X1 + δh1 (X1 , X2 )       x2 = λ2 X2 + δh2 (X1 , X2 ) (2.23)      x 3 = X3  where δ  1 is a perturbation parameter and h1 and h2 are arbitrary functions to be determined. Then the equilibrium equations are given as, T11,1 + T12,2 δ  3  = λ2 3 − 6ν + λ21 (λ1 λ2 )2ν/(1−2ν) (−1 + 4ν) h1,11 + λ21 λ2 (1 − 2ν) h1,22 µ 1 − 2ν    −λ1 (λ1 λ2 )(2−2ν)/(1−2ν) (1 − 4ν) + λ21 + λ22 (−1 + 2ν) h2,12 + O δ 2 = 0,   T21,1 + T22,2 δ  3  = λ1 3 − 6ν + λ22 (λ1 λ2 )2ν/(1−2ν) (−1 + 4ν) h2,22 + λ1 λ22 (1 − 2ν) h2,11 µ 1 − 2ν    (2−2ν)/(1−2ν) 2 2 (1 − 4ν) + λ1 + λ2 (−1 + 2ν) h1,12 + O δ 2 = 0   −λ2 (λ1 λ2 ) (2.24) 2.3.2.1 ν =1/2 When ν = 1/2, the first Piola-Kirchoff stress is T = µ (I1 F − BF) − pF−T . (2.25) where p is the hydrostatic pressure. By assuming   h1 = C1 sin (kX1 ) exp (αkX2 )  (2.26)  h2 = C2 cos (kX1 ) exp (αkX2 )  and p = p0 + δp1 cos (kX1 ) exp (αkX2 ), two eigenvalues were determined to be α1 = 1 and α2 = 1/λ21 using equilibrium equations (2.24). The two eigenmodes gives C2 = −C1 and C2 = −λ−2 1 C1 . 24 Combining these two eigensolutions gives the displacements components on X2 = 0 as  2 kx1 + O δ2    u1 (x1 , 0) = −δ Aλ1 + B sin   λ1 (2.27) kx  u2 (x1 , 0) = δ (A + B) cos 1 + O δ 2    λ1 where A and B are coefficients to be determined by boundary conditions. The Cauchy stress components at the interface between substrate and layer can also be derived as     σ12 1 kx1  2   µ = −2kδ 2λ1 A + λ1 + λ3 B sin λ + O δ   1 1     (2.28) σ 2A 1 kx1  22 = 2kδλ1 + O δ2 .    + λ1 + 3 B cos µ λ1 λ1 λ1 These are identical to the results for incompressible neo-Hookean material (f = 1, ν = 1/2) as provided in [92]. From the stress and displacement components above, the same buckling strain is derived, 2 ¯3 2 ¯2 λ1 k¯4 + 4 (1 + λ21 ) µ ¯k¯ + 24λ31 (1 + λ41 ) µ ¯k + 12λ41 (1 + λ21 ) µ ¯k + 6λ1 (1 − λ21 ) µ ¯2 C L = ¯ 2 µ ¯  12kλ 1 ¯ + λ1 k + λ1 µ ¯ (2.29) Thus the onset of wrinkling and wavelength bifurcation in these two materials are the same so they would have similar wrinkling behavior. In fact, since the invariants I1 and I2 are equal in plane strain problems, the strain energy of Mooney-Rivlin model (Equation (2.6)) degenerates to W = µ (I1 − 3) /2, which describes the incompressible neo-Hookean material. Therefore, these results above hold not only for f = 0 and f = 1, but also for any incompressible material with (0 ≤ f ≤ 1), i.e. the Mooney- Rivlin material. 25 2.3.2.2 ν =1/4 (Blatz-Ko material) When ν = 1/4, the material becomes the well-known Blatz-Ko material. In this d case, solving the equilibrium equation (2.24) with the perturbation assumption (2.26) l leads to an imaginary α, and therefore imaginary displacement. We therefore set the perturbation on displacements to be   h1 = C1 sin (kX1 ) sin (αkX2 )  (2.30)  h2 = C2 cos (kX1 ) cos (αkX2 ) .  d l Figure 2.6: Undeformed configuration of finite substrate with width l and thickness d. Since this deformation has no decay in the X2 direction, it could not represent a deformation in a semi-infinite substrate. However, the deformation could take place in a slab of finite thickness d with traction free surfaces on X2 = 0 & X2 = d and the stretches λ1 & λ2 , as shown in Figure 2.6. From the equilibrium equation, eigenvalues of α are found as,  r ! 1 −16/3 −8/3 −16/3 16/3  8/3 16/3 2 α12 = − 8λ1 − λ1 −36λ1 + 1 − 8λ1 + λ1    1 + λ1   6 r !   2 1 −16/3 −8/3 −16/3 16/3  8/3 16/3 2  α2 =   1 + λ1 − 8λ1 + λ1 −36λ1 + 1 − 8λ1 + λ1 , 6 (2.31) d l 26 λ1 l (a) (b) and deformations can be derived as,  ! 2 8/3 2 8/3  1 + 3α λ L 1 + 3α λ L  x1 = λ1 X1 + δ A −4/3 1 14/3 sin (α1 kX2 ) + B −4/3 2 14/3 sin (α2 kX2 ) sin (kX1 )   λ1 + λ1 α1 λ1 + λ1 α 2   x2 = λ−1/3 X2 + δL (A cos (α1 kX2 ) + B cos (α2 kX2 )) cos (kX1 ) .   1 (2.32) In this case, eigenvalues of α are real numbers and the displacements are meaningful only if λ1 ≤ 0.3724. Therefore, deformations of this wrinkle form may be admissible when the material is compressed beyond that point. The two eigenmodes also provide the first Piola-Kirchhoff stress as,    T11  −1/3  8/3  8/3  sin (α kX ) 1 2 = λ1 − λ−3 + 2πδ A 3 − α12 λ1 −8 + λ1   1 µ α1          sin (α kX )  cos (kX )  2 8/3 8/3 2 2 1  + O δ2  +B 3 − α2 λ1 −8 + λ1    α2   8/3 8/3    λ1 1 + λ1   T12 o sin (kX )   n      2 16/3 2 16/3 1    µ = 2πδ A 3α1 λ1 − 1 cos (α1 kX2 ) + B 3α2 λ1 − 1 cos (α2 kX2 ) 8/3 16/3 λ1 + λ1  + O δ2        λ14/3 sin (kX1 )   T21 3α12 2 + O δ2      = 2πδ A − 1 cos (α1 kX2 ) + B 3α2 − 1 cos (α2 kX2 )    µ 8/3 1 + λ1        4/3 T22  sin (α1 kX2 )  sin (α2 kX2 ) λ1 cos (kX1 )   = 2πδ A 1 − 3α12 + B 1 − 3α22 2      µ 8/3 + O δ . α1 α2 1 + λ1 (2.33) Plugging the first Piola-Kirchhoff stress into boundary conditions (T12 = T22 = 0 on X2 = 0, d) yields,       16/3 16/3 3α12 λ1 −1 3α22 λ1 −1   0   O (δ)   A               3α2 λ16/3 − 1 cos (α kd) 3α2 λ16/3 − 1 cos (α kd) =  0 + O (δ)    .  1 1 1 2 1 2     B     (1 − 3α12 ) sin(α α1 1 kd) (1 − 3α22 ) sin(α α2 2 kd) 0 O (δ) (2.34) 27 Note that T22 vanishes when X2 = 0, which induces no constraint on parameters A and B. Equation (2.34) requires α1 kd and α2 kd to be integer multiples of π, and therefore the wave number k = n1 π/ (α1 d) = (n1 + 2n2 ) π/ (α2 d), where n1 and n2 are positive integers. Different values of n1 and n2 , along with the thickness d, determine different values of wave number k and wave length L, which correspond to different wrinkle modes. From the above equation and Equation (2.31), multi-eigenmodes with identical (n1 + 2n2 ) /n1 correspond to the same stretch ratio λ1 . The maximum value λcr 1 corresponds to the critical wrinkling stretch ratio. This maximum is found to be 0.3724 as (n1 + 2n2 ) /n1 approaches 1. Assuming d = 1, nine eigenmodes asso- ciated with three stretch ratios are plotted in Figure 2.7, with initial coordinates (0 ≤ X1 ≤ 5, 0 ≤ X2 ≤ 1). 28 (a) α1 kd = α2 kd/3 = π (b) α1 kd = α2 kd/3 = 2π (c) α1 kd = α2 kd/3 = 3π (d) α1 kd = α2 kd/5 = π (e) α1 kd = α2 kd/5 = 2π (f) α1 kd = α2 kd/5 = 3π (g) α1 kd = α2 kd/7 = π (h) α1 kd = α2 kd/7 = 2π (i) α1 kd = α2 kd/7 = 3π Figure 2.7: Wrinkling eigenmodes in a finit slab with initial coordinates (0 ≤ X1 ≤ 5, 0 ≤ X2 ≤ 1). Solid lines show wrinkle patterns, while dash lines are unperturbed deformations with the same stretch ratio. 2.3.2.3 0 ≤ ν ≤ 1/2 The compressible Blatz-Ko and the incompressible Mooney-Rivlin materials exhibit different wrinkle patterns. The perturbation form with a sinusoidal variation in the X1 direction of the incompressible material has an exponential decay in X2 direction, while for the case of ν = 1/4, the perturbation is characterized by sinusoidal curves in the X2 direction. Wrinkling of the material with f = 0, unlike that of f = 1, is significantly dependent on compressibility, which was verified by further investigation on the eigenvalue α (Figure 2.8). From Equation (2.26), when α is real-valued, the perturbation form includes exponential decay perpendicular to the surface while 29 imaginary values of α correspond to sinusoidal variations perpendicular to the surface. Complex values of α do not correspond to a perturbation field of the form considered. It might lead to complex perturbation forms but this was not investigated in this paper. Figure 2.8: Characterization of wrinkle patterns for f = 0. (a) in region I, value of α2 is positive, which corresponds to wrinkling in semi-infinite substrate with a thin film layer; (b) in region II, value of α2 is negative and therefore α is imaginary. This corresponds to wrinkling in a finite slab; (c) in region III, α2 is a complex number and might leads to more complex perturbation forms. I) Wrinkling in a semi-infinite solid bonded to a thin film. For the pertur- bation form considered, wrinkling in a bilayer system is possible for certain values of ν & λ1 , that fall in region I of the ν − λ1 graph in Figure 2.8. In this region, wrinkling behavior is similar to that of generalized neo-Hookean material (parameter f = 1). The critical buckling strain is determined by the mechanical properties of substrate/thin film as well as the pre-stretch ratio. With lower pre-stretch ratio and lower Poisson’s ratio, the critical buckling strain decreases and thus wrinkling is more easily induced. Moreover, as the material becomes more compressible, the pre-stretch turns out to play a less important role in wrinkling. The most significant effect of compressibility is the transition of wrinkling patterns. As shown in Figure 2.8, the borderline between region I and III marks the stretch ratio where α becomes a com- plex number. As a result, a wrinkle pattern of the form considered is not possible in region II. This line climbs up with Poisson’s ratio decreasing and intersects λ1 = 1 30 at ν = 0.4375. This indicates wrinkling of bilayer system only exists in the material with Poisson’s ratio 0.4375 < ν ≤ 0.5 given no pre-stretch. With higher pre-stretch of the substrate, this wrinkle form is not admissible at lower Poisson’s ratio. II) Wrinkling in a finite slab. Wrinkling deformations of the form considered are characterized by sinusoidal waves in both X1 & X2 directions. As shown in Figure 2.8 this perturbation form is possible only for values of ν & λ1 that fall in region II. The borderline between II and III represents the stretch ratio at which α changes into an imaginary number. Because the perturbation does not decay in the X2 direction, the deformation is not admissible in a semi-infinite solid. However, in a finite slab with traction free surfaces under uniaxial compression, this wrinkle pattern is admissible for all values of ν < 0.5. The borderline between regions II & III serves as the stretch ratio beyond which wrinkling of the form considered is admissible. To summarize, compressibility influences wrinkling of the material with f = 0 in many ways, particularly in wrinkling onset. Wrinkling in a finite slab (without a thin film on top of the surface) occurs only in compressible materials (region II). Wrinkling in semi-infinite bilayer systems (region I) is affected by the substrate compressibility as well. This phenomenon provide an option for wrinkling control in various applications such as stretchable electronics and microfluidic systems. 2.3.3 0≤f ≤1 Theoretically the two wrinkle patterns (in the finite and semi-infinite geometries) could be induced for materials with parameter f between 0 and 1. Figure 2.9 shows the characterization of wrinkle forms for materials with different values of f . Similar to that of f = 0, wrinkling in the bilayer system is possible for values of ν & λ1 in region I and wrinkling in a finite-slab occurs only for values of ν & λ1 in region II. No wrinkling of the type considered occurs in region III. As f gets closer to 1, region I expands to lower Poisson’s ratio and higher compressive strain whereas regions II 31 and III shrink and finally vanish (at f = 1). Therefore, the parameter f influences wrinkling in the generalized Blatz-Ko material in two ways. The critical strain for wrinkling in a finite slab is increased as f goes to 1. This form of wrinkling is not possible at f = 1. Wrinkling in a semi-infinite solid bonded with a thin film is not admissible for the material with high compressibility when f is close to 0. (a) f = 0 (b) f = 0.5 (c) f = 0.9 Figure 2.9: Characterization of wrinkle patterns for different f . 2.4 Conclusion Through theoretical analysis, we investigate the effects of material compressibility on certain types of wrinkle formation in a class of non-linearly elastic materials. Two wrinkle patterns are discussed and the effects of the compressibility, parameter f as well as pre-stretch are studied. A bilayer system consists of an incompressible thin film and a compressible semi- infinite substrate. When it is uniaxially compressed, wrinkling of the form considered is induced at the critical buckling strain B . Firstly, the substrate is taken to be a compressible neo-Hookean material (f = 1) to determine the effect of Poisson’s ratio and pre-stretch. We find that increased incompressibility in the substrate, as well as pre-stretch, slightly increases B . Secondly, for materials with (0 ≤ f < 1), wrinkling of the bilayer system is suppressed for more compressible materials. This effect is maximized when the parameter f = 0, where wrinkle patterns are not possible for values of ν ≤ 0.4375. Therefore, by using generalized Blatz-Ko material with various 32 compressibility, surface instabilities of the type studied can be suppressed. This could have applications in settings where adhesion control is important. For certain ranges of the parameters (f and ν), wrinkling is admissible in a slab of finite thickness but infinite horizontal extent (without a thin film). The perturbation form is sinusoidal in both X1 and X2 directions, and traction free surface conditions may be enforced on the upper and lower surfaces. The critical buckling strain in these materials is found to be higher than that of the wrinkling in a semi-infinite system and it is found to be more significantly influenced by the compressibility of material. The wrinkling formation is admissible in the compressible generalized Blatz-Ko material with (0 ≤ f < 1). As f gets closer to 0, wrinkles could be induced in materials with higher incompressibility. However, as [18] pointed out, wrinkling in a half-space homogeneous incompressible neo-Hookean solid is highly unstable. It remains for future work to explore the stability of wrinkling patterns in finite slabs, with further experimental and computational study. 33 Chapter 3 Wrinkling in Generalized Blatz-Ko Materials II: Instability and Imperfection-sensitivity 3.1 Introduction In the previous chapter, the effects of material compressibility on the formation of surface wrinkles were investigated. A linear bifurcation analysis of the generalized Blatz-Ko material was performed. Apart from the widely studied wrinkling form in the incompressible neo-Hookean material, a new wrinkling pattern, characterized by sinusoidal displacements in both horizontal and vertical directions, was predicted theoretically. For certain ranges of the parameters in the generalized Blatz-Ko model, this pattern is admissible in a slab of finite thickness. It has been shown that the wrinkling in the semi-infinite homogeneous incompressible neo-Hookean solid is highly unstable but it could be stabilized and occur in experiments when a stiff thin-film is firmly attached to the solid [18]. It was found that wrinkling in the bilayer system is influenced by the pre-stretch on the substrate, the film/substrate stiffness ratio, 34 and the variation in substrate properties [92, 19, 28]. In this section, the stability of the bi-sinusoidal wrinkling patterns in compressible materials is investigated. If it is observed, does the wrinkling deformation match the theoretical prediction? If it is not stable, what method can be employed to stabilized the instability? To answer the previous questions, we perform a second-order displacement pertur- bation analysis with the Koiter method [59]. The elastic system is completely de- scribed by its potential energy (u, λ1 ), where u is the admissible displacement and λ1 is the stretch ratio parameter. For a given stretch, the corresponding equilib- rium equation is found, by taking the Gateau derivative of  with respect to u, to be ,u (u, λ1 )δu = 0 . At the critical wrinkling point, there exists such admissible δu that ,uu (u, λ1 )δuδu = 0. Employing the Lyapunov-Schmidt-Koiter decomposi- tion method, one can write the Taylor series expansion of (u, λ1 ) around the critical point. By considering up to the second order deformation, stability of the wrinkling in finite slab is analyzed in this chapter. The finite element method (FEM) is employed to analyze the wrinkles and their sta- bility. A two dimensional slab with finite thickness was modeled and the slab is com- pressed in X1 direction with periodic boundary condition applied. The generalized Blatz-Ko material with various parameters (f &ν) are considered. In this Chapter, we demonstrate the critical wrinkling strain, employing extensive finite element analyses, and compare the numerical and theoretical predictions. In addition, post-wrinkling process is also investigated. While the wrinkling in bilayer systems may grow and transition to other instability patterns [112], the wrinkle amplitude in the finite slab does not become dramatically larger. By performing Fourier transform analyses, we calculated the dominant frequencies of the wrinkling deformation, and compare those frequencies with the theoretical predictions from the preceding chapter. In addition, imperfection sensitivity is studied through finite element simulations. Creasing is observed when an initial imperfection of large amplitude is applied. In 35 contrast to wrinkles, creasing is characterized by localized deformation. The surface folds into a self-contact region with a sharp tip [16]. The critical strain for creasing in incompressible neo-Hookean material was measured to be 0.35 ± 0.07 by bending elastomer surface experimentally [34], and it was further verified numerically to be 0.35 [48]. In this chapter, the critical creasing strain for the generalized Blatz-Ko (f = 0, 0.5&1) slab is investigated via finite element simulations. Elastic instability of various kinds has applications in engineering field from stretch- able electronics [52] to control surface chemistry [53]. By understanding the critical conditions and imperfection sensitivity of wrinkles and creases, one can control the instability types of the soft materials and forecast formation of such patterns in na- ture. This chapter is organized as follows: We study the stability of the wrinkling in finite slabs through a second-order perturbation analysis via Koiter’s method in Section 3.2. We validate theoretical wrinkle formation through finite element simulations in Section 3.3. In particular, we employ ABAQUS/standard to simulate the initiation and growth of the wrinkles and creases. The post-wrinkling deformation was an- alyzed through a Fourier transformation analysis in Matlab and compared against theoretical prediction. In addition, the critical condition for creasing was also studied numerically. 3.2 Stability analysis of wrinkles in finite slabs Using a linear perturbation analysis, we examined bifurcation of a Blatz-Ko finite slab to a bi-directional sinusoidal wrinkling pattern under plane-strain uniaxial com- pression. The critical compressive strain for the wrinkle to develop was determined, however, whether the wrinkle form is stable or not remained hidden. Here, we in- vestigate the stability of these sinusoidal wrinkles by carrying out Koiter’s nonlinear 36 bifurcation analysis with second-order perturbation solutions. The first order defor- mation is also provided here in a systematic fashion. 3.2.1 Linear perturbation analysis The strain energy density of the generalized Blatz-Ko material is given by equation (2.3). For the simplicity of the problem and without loss of the generality, here we present the solution for the special Blatz-Ko material (f = 0, ν = 1/4). From the elastic potential, the first Piola-Kirchhoff stress can be derived as equation (2.22). Given an initial (zeroth-order) uniform stretch deformation with gradient of (F11 = λ1 , F22 = λ2 and F12 = F21 = 0), the traction-free boundary condition 0 −1/3 T22 = 0 on the slab surfaces indicate λ2 = λ1 . To investigate stability of the wrinkling pattern, perturbation analysis is performed. The displacement is assumed in a similar way as described in Chapter 1, except that it is now expanded up to the second order as, u(X1 , X2 , λ1 ) = uc + ζu(1) + ζ 2 u(2) + O(ζ 3 ), (3.1) where uc is the displacement field at the critical wrinkling point, u(1) and u(2) are the first and the second order eigenfunction displacement fields, and ζ  1 is the perturbation factor. 37 The equilibrium equation is given by, T + T12,2 δ   11,1  = 4 (1) 4/3 (1) 8/3 (1) [3u1,11 + (λ41 + λ1 )u2,12 + λ1 u1,22 ]     µ λ1 δ 2 h (1)     (1) (1) 4/3 (1) 8/3 (1) + 5 u1,1 −12u1,11 − λ41 u2,12 − 3λ1 u2,12 − 2λ1 u1,22        λ1    (1) 4 (1) 4/3 (1) 8/3 (1) 4 (1)  + u1,2 −λ1 u2,11 − 3λ1 u2,11 − 4λ1 u1,12 − 2λ1 u2,22          (1) 8/3 (1) (1) 4/3 (1) 8/3 (1) 16/3 (1)  + u2,1 −2λ1 u2,11 − 2λ41 u1,12 − 6λ1 u1,12 − λ1 u2,22 − 3λ1 u2,22          (1) 16/3 (1) 8/3 (1) (1)  + u2,2 −3λ1 u2,12 − λ1 u2,12 − 2λ41 u1,22        i (2) 7/3 (2) 11/3 (2)  + 3λ1 u1,11 + (λ51 + λ1 )u2,12 + λ1 u1,22         =0    T21,1 + T22,2 δ 4/3 (1) 8/3 (1) (1)   = 8/3 [λ1 u2,11 + (1 + λ1 )u1,12 + 3λ41 u2,22 ]     µ λ1  2 δ   h   (1) 4/3 (1) 8/3 (1) (1) u1,1 −2λ1 u2,11 − λ1 u1,12 − 3u1,12     + 11/3    λ1      (1) (1) 8/3 (1) 4/3 (1) (1) 8/3 (1) + u1,2 −3u1,11 − λ1 u1,11 − 2λ1 u2,12 − 6λ41 u2,12 − 2λ1 u1,22           (1) 4/3 (1) 8/3 (1) 4/3 (1) 4 (1) + u2,1 −2λ1 u1,11 − 4λ1 u2,12 − λ1 u1,22 − 3λ1 u1,22           (1) 16/3 (1) 4 (1) 8/3 (1) 4/3 (1) + u2,2 −12λ1 u2,22 − 3λ1 u1,12 − 2λ1 u2,11 − λ1 u1,12        i  7/3 (2) 11/3 (2) (2) 5 (2)     + λ1 u2,11 + λ1 u1,12 + λ1 u1,12 + 3λ1 u2,22      =0 (3.2) The first-order equations are given as,   a1 u(1)  (1) (1) 1,22 + b1 u2,12 + c1 u1,11 = 0 (3.3)  a2 u(1) + b2 u(1) + c2 u(1) = 0  2,22 1,12 2,11 38     8/3 4/3 8/3 4/3 where a1 = λ1 , b1 = λ1 + λ41 , c1 = 3, a2 = 3λ41 , b2 = 1 + λ1 , c2 = λ1 . The displacements can be solved through a function Φ (X1 , X2 ) as, 1  (1)  u1 = − Φ,12   a1 c 1   (3.4) (1) 1 1 1  u2 = Φ,11 + Φ,22 .   b 1 a1 c1 Then the first equilibrium equation vanishes and the second equation becomes a fourth order partial differential equation,   c2 a2 b2 c2 a2 Φ,1111 + − + Φ,1122 + Φ,2222 = 0 (3.5) a1 b 1 a1 b 1 a1 c 1 b 1 c 1 b1 c 1 The general solution has the form, 4 a1 c 1 X Φ=− 2 Ci F1 cos (kX1 ) exp (αi kX2 ) (3.6) αk i=1 8/3 1+3α2 λ1 L where F1 = −4/3 4/3 α , k, L are the wavenumber and wavelength in the undeformed λ1 +λ1 configuration, and the coefficients Ci are determined from the boundary conditions. Since the partial differential equation given above leads to imaginary α, and therefore imaginary displacements. Therefore we set the X2 part of Φ to be sinusoidal functions as follows, 2 4 a1 c 1 X a1 c 1 X Φ=− 2 Ci F1 cos (kX1 ) cos (αi kX2 ) − Ci F1 cos (kX1 ) sin (αi kX2 ) αk i=1 αk 2 i=3 (3.7) 39 Enforce T12 = 0 on the surface and C3 , C4 are determined to be 0. Then the first order displacement field is obtained as,  2 X  (1)  u1 =   Ci F1 sin (kX1 ) sin (αi kX2 )  i=1 2 (3.8)  X  (1)  u2 = Ci F2 cos (kX1 ) cos (αi kX2 ) .   i=1 and F2 = L. C1 , C2 can be further determined in a similar way as the preceding chapter. However, since we are interested whether the present wrinkling form is stable, equation (3.8) is used for the following calculations. 3.2.2 Second-order perturbation analysis for nonlinear sta- bility of a wrinkle In this section, nonlinear stability of the wrinkles is examined by second-order per- turbation analysis. Substituting the first order displacement field into the second order equilibrium equation, we can get a system of inhomogeneous partial differential equations for the second-order displacement fields.   a1 u(2)  (2) (2) 3 1,22 + b1 u2,12 + c1 u1,11 = k (A1 + B1 cos (2αkX2 )) sin (2kX1 ) (3.9)  a2 u(2) + b2 u(2) + c2 u(2) = k 3 (A2 + B2 cos (2kX1 )) sin (2αkX2 )  2,22 1,12 2,11 1 5 where A1 = −3F12 /λ1 − 3F1 F2 αλ13 /2 + (F2 2 + F1 2 α2 )λ13 /2 − F1 F2 αλ31 /2, B1 = 1 8 5 8 −3F1 F2 αλ13 − F1 F2 (α + α3 )λ31 + 3F1 2 (2 + α2 λ13 )/2λ1 + F2 2 λ13 (1 + α2 + 3α2 λ13 )/2, 1 4 8 A2 = αλ13 (F1 F2 α − (F2 2 + F1 2 α2 )λ13 + 3F1 F2 αλ13 + 6F2 2 α2 λ41 )/2, and B2 = 8 8 4 8 8 (3F2 2 αλ13 (1 + 2α2 λ13 ) − 2F1 F2 λ13 (1 + α2 + 3α2 λ13 ) + F1 2 α(3 + (1 + α2 )λ13 )/2λ1 . 40 Denote that  (2)  u1 = v1   k2  B1  (3.10) (2) u = v − A1 X2 + sin(2αkX2 ) cos(2kX1 ).  2   2 2b1 2αk The equations become,     a1 v1,22 + b1 v2,12 + c1 v1,11 = 0    a2 v2,22 + b2 v1,12 + c2 v2,11 = k 3 (A2 + B2 cos (2kX1 )) sin (2αkX2 )  k2  2   c    2 − 4k A1 c2 X2 + 2B1 k + 2αa2 sin(2αkX2 ) cos(2kX1 )   2b1 α (3.11) Follow similar procedure as in the first-order solution, by assuming 1   v1 = − Ψ,12   a1 c 1   (3.12) 1 1 1  v2 = Ψ,11 + Ψ,22 .   b 1 a1 c1 we can derive a fourth-order inhomogeneous partial differential equation as,   c2 a2 b2 c2 a2 Ψ,1111 + − + Ψ,1122 + Ψ,2222 a1 b 1 a1 b 1 a1 c 1 b 1 c 1 b1 c 1 = k 3 (A2 + B2 cos (2kX1 )) sin (2αkX2 ) (3.13) k2  2 c 2   − 4k A1 c2 X2 + 2B1 k + 2αa2 sin(2αkX2 ) cos(2kX1 ). 2b1 α The solution to the above equation is determined by combining a particular solution and the general solution to homogeneous equations as follows, Ψ(X1 , X2 ) =F3 sin(2αkX2 ) + F4 X2 cos(2kX1 ) + F5 X2 cos(2kX1 ) cos(2αkX2 ) 6 8 a1 c1 X ¯  ¯  a1 c 1 X ¯ 1 sin αi kX  ¯ 2 .  − C i F 1 cos kX1 cos αi kX 2 − C i F1 cos kX αk¯2 i=5 αk¯2 i=7 (3.14) 41 b1 c1 A2 a1 c1 (αb1 B2 −B1 (α2 a1 +c2 )) where the parameters are given as F3 = 16α4 a2 k , F4 = − A18a1 , F5 = 16α2 (2α2 a1 a2 +a1 c2 +a2 c1 −b1 b2 ) , and Ci and k¯ are coefficients to be determined by boundary conditions. Therefore, the second-order displacements can be expressed as,  (2) 2k u1 = (F4 + F5 (cos (2αkX2 ) − 2αkX2 sin (2αkX2 ))) sin (2kX1 )   a1 c 1     k  (2) u2 = − 16a1 α3 kF3 sin (2αkX2 ) + (2c1 αk(Aa1 + 8F4 )X2 + a1 (Bc1    4a b c 1 1 1 α   +16α2 F5 sin (2αkX2 ) + 16αk(c1 + a1 α2 )X2 F5 cos (2αkX2 ) cos (2kX1 )      (3.15) With above results, the displacement field and stretch ratio in the post-buckling process can be approximated by the following Taylor expansion up to the second order,   u(X1 , X2 , λ1 ) = uc + ζu(1) + ζ 2 u(2) + O(ζ 3 )  (3.16)  λ1 = λc + ζλ(1) + ζ 2 λ(2) + O(ζ 3 )  1 1 1 (1) (2) where the parameters λ1 and λ1 in the λ1 expansion determine the stability of the wrinkling form considered at the initial post-bifurcation process. By Koiter’s theory, (1) (2) the wrinkle is supercritical or stable if λ1 = 0 and λ1 < 0 , and it is subcritical or (1) (2) (1) unstable if λ1 = 0 and λ1 > 0 . If λ1 6= 0 , the wrinkle is transcritical. Here, we present a nonlinear stability analysis by calculation of the energy variation in the initial post-bifurcation state around critical compressive strain. Let W (u) be the strain energy per unit volume of the deformed material, then Z = (W (u) − W (uc ))dV (3.17) V is the energy change relative to the uniform state associated with imposed λ1 in the Blatz-Ko slab. The variational form of the equilibrium equation ,u (u, λ1 )δu = 0 , 42 along the post-bifurcation equilibrium path can be expanded as follows, 1 c ,u (u, λ1 )δu =c,u δu + c,uu δuδu + c,uλ1 δu∆λ1 + ,uuu δuδuδu + 2c,uuλ1 δuδu∆λ1 2  1 c +c,uλ1 λ1 δu(∆λ1 )2 +  δuδuδuδu + 3c,uuuλ1 δuδuδu∆λ1 6 ,uuuu +3c,uuλ1 λ1 δuδu(∆λ1 )2 + c,uλ1 λ1 λ1 δu(∆λ1 )3 + O(ζ 4 )   (1)  1  (2) =c,u δu + ζ c,uu u(1) + c,uλ1 λ1 δu + ζ 2 2c,uu u(2) + 2c,uλ1 λ1 + c,uuu (u(1) )2 2 (1) (1)  1  (1) +2c,uuλ1 λ1 u(1) + c,uλ1 λ1 (λ1 )2 δu + ζ 3 6c,uuu u(1) u(2) + 6c,uuλ1 λ1 u(2) 6 (2) (1) (2) (1) +6c,uuλ1 λ1 u(1) + 6c,uλ1 λ1 λ1 λ1 + c,uuuu (u(1) )3 + 3c,uuuλ1 λ1 (u(1) )2  (1) (1) +3c,uuλ1 λ1 (λ1 )2 u(1) + c,uλ1 λ1 λ1 (λ1 )3 δu + O(ζ 4 ) = 0 (3.18) where the superscript ‘c’ means that the functions are evaluated at ζ = 0. In the above equation, ,c,u u(1) implies the variation of the functional caused by variation of the displacement field δu = u(1) . From the equilibrium equation, c,uλ1 δu = c,uλ1 λ1 δu = c,uλ1 λ1 λ1 δu = 0 along the post-buckling path. In addition, c,uu u(1) u(2) vanishes because of the orthogonality between the first- and the second-order displacement fields. By choosing δu = ηu(1) with a small variational parameter η in the above equation, the (1) (2) coefficients λ1 and λ1 are determined as, 1 c (1) ,uuu (u(1) )3 λ1 = − 2 c ,uuλ1 (u(1) )2   1 c (1) 1 c (1)  6 ,uuuu (u(1) )4 + λ1  2 ,uuuλ1 (u(1) )3 + c,uuλ1 u(1) u(2) + 12 c,uuλ1 λ1 (u(1) )2 λ1 + c,uuu (u(1) )2 u(2) (2) λ1 =− c,uuλ1 (u(1) )2 (3.19) By substituting the wrinkling displacement fields from equation (3.8) and (3.15) and (1) (2) integrating the variations per wavelength, λ1 is solved to be 0 and the value of λ1 (2) is plotted in Figure 3.1. As shown in the Figure 3.1, λ1 is negative, which indicates 43 a supercritical bifurcation. Therefore, the wrinkling of the form considered is stable. (a) α1 (b) α2 (2) Figure 3.1: λ1 as a function of the stretch ratio λ1 . 3.3 Finite element analysis of post-bifurcation wrinkling and creasing A finite element model has been employed to simulate uniaxial compression of the generalized Blatz-Ko slabs. The critical wrinkling pattern and the post-bifurcation evolution of the wrinkles are shown for the compressible materials with different values of the parameters f &ν. Since the critical and post-bifurcation displacement fields are independent of the shear modulus, µ is taken to be 1 Mpa for all simulation cases. Plane strain finite element simulations have been performed with ABAQUS/standard. Figure 3.2 shows a schematic of the finite element model in its initial configuration. The simulation is of a rectangular slab with thickness d = 1 and length L = 5d, where the out of plane depth is 1. Traction free boundary conditions are enforced on the top and bottom surfaces, and periodic boundary conditions are applied at the lateral sides. The Blatz-Ko material model has been implemented for ABAQUS simulation through a user subroutine, UHYPER, defining the strain-energy potential and its derivatives with respect to invariants [41]. The mesh size is taken to be small 44 enough to ensure it doesnt influence the simulation results. Bilinear plane strain element CPE4 is adopted to simulate the compressible materials and CPE4H is used for incompressible solids. Figure 3.2: FE model of the finite slab with d=1 and L=5. Prior to compression, a sinusoidal vertical displacement perturbation is introduced to the top surface to trigger an instability. The perturbation is expressed as,   2πkX1 up2 = ζ cos (3.20) L where k = 2, 4, 6, 8 . . . and ζ is the imperfection amplitude. Then the model is compressed by small incremental displacement-controlled loading on the right face. Finite element simulations show that wrinkling in finite slabs of Blatz-Ko material, distinct from that in half infinite neo-Hookean solid, is characterized by sinusoidal displacements in both horizontal and vertical directions. 3.3.1 Critical wrinkling state Under isoperiodic compression, the slab is first uniformly deformed with perturbation localized on the top surface. At the bifurcation point, sinusoidal wrinkling develops throughout the model. Finite element simulations verify that the wrinkling form considered occurs in a subclass of the generalized Blatz-Ko material, as theoretically predicted previously. To illustrate the wrinkling process, here we present analysis of the material with (f = 0&ν = 0.2) as an example. The initial perturbation is chosen 45 to be up2 = 0.001 cos 4πX1  L . Figure 3.3 shows the finite element simulation results near bifurcation point. Figure 3.3: Wrinkling deformation of the material with (f = 0&ν = 0.2) around critical point ( = 0.6437). Contour plot demonstrates the deformed maximum principal nominal strain, max p , distribution. (a)(b) prior to the onset of wrinkles ( = 0.4343, 0.5943), the perturbation is first localized near the surface and decays as away from the surface; (c) beyond the critical point ( = 0.7143), the non-decaying wrinkling deformation, characterized by evenly arranged rectangular cells, develops. The critical wrinkling strain that corresponds to the onset of wrinkles was determined from the finite element simulation results to be 0.6373 ± 0.01. The error range is due to the displacement-controlled loading increment ∆ = 0.02. This critical strain is in good agreement with the previous theoretical predictions by solving for the eigenvalue α, according to which crit = 0.6437. It is also found to be independent of the imperfection wavelength and amplitude, as discussed in Section 3.3.3. Similar wrinkling behaviors were observed for the generalized Blatz-Ko material with different 46 parameters. Figure 3.4 shows the critical wrinkling stretch ratio for the material with parameter values of f = 0 and f = 1/2 and with various levels of material compressibility. The figure confirms that the wrinkling forms are predicted by finite element analysis and align with the theoretical derivation. (a) f = 0 (b) f = 0.5 Figure 3.4: Critical wrinkling strain in materials with f = 0 and 0.5. Solid line is the theoretical predicted critical stretch ratio; markers show strain for uniform deformation to wrinkling state transition determined by FEM simulations, the error bars are due to the displacement-controlled loading steps ∆. 3.3.2 Fourier transformation analysis of post-wrinkling de- formations Beyond the critical strain, the slab is compressed further to study the post-bifurcation behavior and to investigate whether wrinkling of the form considered is stable. In the FE simulation, the slab is deformed to a compressive ratio that is much larger than the critical strain and the wrinkling pattern remains. To analyze wrinkling deformations in the post-buckling processes, displacements of the finite slab are extracted from the ABAQUS output database with python scripts. Fourier transform analysis in MATLAB is then employed to capture the frequency of the deformation in both horizontal and vertical directions. Theoretical predictions show that the wrinkling deformation is characterized by sinusoidal functions in both horizontal and vertical directions. For each sinusoidal eigenmode, the wave number 47 is k/L and αk/L in X1 and X2 directions, respectively. Therefore, the ratio between the frequencies in two directions fy /fx is determined to be α. (a) (b) (c) (d) Figure 3.5: Post-buckling deformation of Blatz-Ko material at a stretch ratio λ1 = 0.2941. (a) ABAQUS simulation result; (b) Fourier transformation analysis. X and Y axis show frequencies of wrinkling deformations in the X1 &X2 directions, respectively. Solid lines correspond to theoretical predictions and dots show FFT calculation of the finite element simulation output;(c) recover wrinkling deformation with the first several eigenmodes; (d) plot of fy for fx equals the frequency of initial perterbation. The simulation results for materials with different values of the parameters f &ν are investigated. With the two dimensional fast Fourier transformation, we were able to get the dominant frequencies. In the frequency domain, only points with a value of higher than 10% of the highest one were taken into account. As shown in Fig- ure 3.5, those first several eigenmodes are sufficient to represent the overall wrinkling deformation. Its shown that the frequency fx of the wrinkling deformation in the horizontal direction is associated with the frequency of perturbation, especially in the initial-post-bifurcation state. However, the ratio of frequencies in the horizon- 48 tal and vertical directions are in good agreement with theoretical predictions. As before, we present analysis of classical Blatz-Ko material (f = 0&ν = 1/4) as an example. The initial perturbation is set as up2 = 0.0001 cos 4πX  L 1 . The critical wrin- kling strain is determined to be 0.6303 ± 0.01, where the theoretical value is 0.6276. Post-buckling analysis reveals that the wrinkling mode remains stable as the overall compressive strain increases. Fourier Transformation of vertical displacements are performed and the dominate frequencies (at  = 0.7059) with highest amplitude are plotted in Figure 3.5. Linear perturbation analysis demonstrates that the first order vertical displacement is δL (A cos(α1 kX2 ) + B cos(α2 kX2 )) cos(kX1 ) . Hence the fre- quency ratio is theoretically predicted as α1 = 2.1067 and α2 = 12.4093. As shown in the figure, the simulation results locate around the theoretical-predicted line. More- over, this conclusion is held for all post-buckling steps for all admissible materials . 3.3.3 Sensitivity to the perturbation: wrinkling and creasing Finite element analysis was performed to investigate the sensitivity of the wrin- kling form considered to the initial perturbation. The perturbation form is given as up2 = ζ cos 2πkX  L 1 . As mentioned in the previous section, the critical wrinkling strain remains the same for different wavenumbers (k = 2, 4, 6, 8) in the initial pertur- bation setting. It is found that changing the initial wavenumber will induce different wrinkling deformations, as shown in Figure 3.6. The dominant frequency in X1 di- rection is determined from the imperfection form applied. However, the frequency ratios of fy /fx are always in agreement with theoretically predicted α. The wrinkling in the finite slab is found to be sensitive to the imperfection amplitude. Although the wrinkling pattern remains the same when the amplitude ζ is within a certain range, creases will develop if a larger ζ (typically an order of 10−3 d) is used. As shown in Figure 3.7, a 2D square model was used to simulate the wrinkles and creases 49 Figure 3.6: Imperfection sensitivity test for the material with parameter (f = 0&ν = 0.2). The perturbation is set to be up2 = ζ cos 2πkX L 1 , where k = 4, 6, 8 and ζ = 0.0002. For all the imperfection cases, finite element simulation results at compressive strain  = 0.7143 are considered. Contour plots of maximum principal nominal strain (max p ) distribution and corresponding deformation frequency analysis results are plotted. on the surface. The slab is compressed in X1 direction with vertical displacement on the bottom surface fixed. When an imperfection of 10−4 d is applied, wrinkles are observed. While an imperfection of 10−3 d is applied, creases are observed. The critical strain for creasing in the generalized Blatz-Ko material is determined from the finite element simulation and shown in Figure 3.8. Distinct from the wrinkling strain, the critical strain for creasing is not significantly influenced by the parameter f and the compressibility of the material. As the material becomes more compressible 50 Figure 3.7: Contour plots of wrinkles and creases in a 2D square model. Undeformed dimension is 1 × 1. (a) ζ = 10−4 ,  = 0.6812; (b) ζ = 10−3 ,  = 0.3959. Figure 3.8: Critical stretch ratio for creasing in the generalized Blatz-Ko material. and f gets closer to 1, the critical stretch ratio becomes smaller and thus it requires a larger compressive strain to induce the creasing. 3.4 Conclusion In this chapter, I explore the stability of the wrinkling pattern in finite slabs theoreti- cally and numerically. A nonlinear bifurcation analysis via Koiter method is employed to analyze the wrinkle deformations. The bifurcated displacements are solved in a systematically to the second order. By calculation of the energy variation in the 51 initial post-bifurcation state near a critical compressive strain, it is shown that the critical bifurcation point is supercritical and the wrinkles are stable. Extensive finite element simulations are also performed to validate and investigate the post-buckling behaviors. FE analyses demonstrate that the critical wrinkling strain is in agreement with the preceding derivations. A Fourier transformation analysis shows that the wrinkling deformations that emerge in the FE analyses follow the theoretically pre- dicted sinusoidal functions and frequencies. The sinusoidal periodicity is dependent on the initial perfection settings. It should be noted that the wrinkle amplitude does not grow dramatically as in the film-substrate case. In addition, the instability type is not sensitive to the imperfection amplitude on the order of 10−4 . However, a larger imperfection amplitude will induce creasing at a lower strain. The results may inspire designing new systems with various surface modes and lead to a new way for control of the instabilities. 52 Chapter 4 A Metallurgical Phase Transformation Framework Applied to Additive Manufacturing Processes 4.1 Introduction A customizable general simulation framework was presented and validated for a wide spectrum of additive manufacturing (AM) processes (laser and electron beam powder bed fabrication, direct energy deposition, arc welding, polymer extrusion, ink jetting) as implemented on the Dassault Systemes 3DX Platform based on new Finite Element (FE) technology implemented in Abaqus [106]. The framework allows for: 1) arbitrary meshes of CAD representations; 2) exact specification in time and space of processing conditions (e.g., powder addition, laser trajectories, dwell times, etc.); 3) precise tracking of the progressive addition of raw material to each element in the mesh via complex geometric computations; 4) precise integration of the moving energy 53 sources (e.g., laser, electron beams, arc welds, high temperature polymer extrusion) and; 5) automatic computation of the continuously evolving convection and radiation surfaces. In this chapter, in conjunction with the general simulation framework above, a generic metallurgical phase transformation framework applicable to arbitrary metal alloys is introduced. The framework accounts for phase transformations (physical state changes) from stock feed or raw material (e,g., powder) via melt- ing/vaporization/solidification followed by metallurgical solid state phase transfor- mations induced by either rapid heating or cooling events associated with typical 3D printing sequences but also with slower-rate temperature evolutions as associ- ated with heat treatment applications. The user systematically defines all possible transformations that can take place via a parent-children paradigm (e.g., austenite to martensite); the temperature conditions when transformations can occur with associated TTT diagrams; and whether the transformations are reversible, diffu- sional, or non-diffusional. The framework predicts temperature evolutions, calibrates JMA/KM models, and predicts phase transformations during printing and heat treatment. Grain growth and grain aspect ratio assessment models are also included. The quality of additively manufactured parts is often hindered by the inherent nature of the manufacturing process. It is not uncommon that the mechanical properties of AM parts are inferior to those manufactured via classical manufacturing techniques and therefore it is desirable to leverage numerical models for the reliable prediction of mechanical properties of printed parts. In this work we adopt existing empirical plasticity models [27, 40, 72, 49] for a Ti-6Al-4V alloy that connect elastic behavior, yield strength/strain, and ultimate strength/ductility limits to phase volume fraction content and grain size information. The models utilize the predictions described above to estimate, in the context of a Ramberg-Osgood plasticity model, the hardening behavior of the as-printed and heat-treated alloys. 54 The chapter is organized as follows. In Section 4.2, a brief review of the general FE-based framework utilized for thermo-mechanical analysis of the printing process is performed. In Section 4.3, the generic metallurgical phase transformation model for any alloy is presented, and two detailed exemplifications are given in Section 4.4. Section 4.5 describes the empirical plasticity models that have been used. Finally, in Section 4.6, conclusions are drawn. 4.2 Finite element simulation of powder-bed fu- sion additive manufacturing New physics-based FEM formulations and models implemented in software Apps from DASSAULT SYSTEMES R2018x [1] were used in this study. It is not the focus of this paper thesis to describe these features in detail, but key aspects of the framework include: • Meshing: the geometric shape of the part whose 3D printing process is to be modeled is first discretized with finite elements. Arbitrary mesh densities can be used as generated by existing pre-processors. This overcomes the challenge of producing uniform, voxel-based or layer-conforming finite element meshes which are difficult, if not impractical, to generate for complex, topology optimised geometries that are often of interest to manufacture using powder-bed fusion processes. • Machine information: the 3D printing machine related information (e.g., powder recoating sequence, laser scan path, material deposition of the printing head, arc weld power, etc.) is pre-processed with no loss of accuracy from actual data as used by the physical machines. 55 • A path intersection module: the path intersection algorithm is used to sweep through the finite element mesh with the laser scanning sequence (or tool path) and associated heat source configurations. The intersection can be based on either the original shape of the part or the current shape of the predicted de- formed/distorted shape of the part during the analysis. • Progressive element activation: at any given point during the simulation, any finite element could be completely filled with matter, partially filled with matter, or completely empty. The calculation method precisely keeps track of this evolution, keeps track of the mass inventory and distribution to account for the addition of the material during the printing. This is a significant step-change in AM process simulation technology from conventional element birth techniques. • Progressive heating computations: at any point in time, heat bursts are com- puted by taking into account the actual path and power distribution of the heat source (the laser beam in this case). An arbitrary number of heating events (characterized as a sequence of heat fluxes at given locations) are computed for an accurate representation of the heating source in both time and space. • Progressive cooling via convection and radiation: partial facet areas are com- puted to allow for a precise assessment of cooling-related heat fluxes, regardless of the finite element discretization. Radiation and convective heat transfer can be modeled on a continuously evolving surface that reflects the current shape of the part at any given point during the print. As indicated previously, this represents a step-change over confirming element activation methods. • Openness: a comprehensive set of new APIs are available to access information computed inside the code. While solutions for all major AM processes are provided, users can tailor physics modeling needs corresponding the particular process, machine environment, etc.. 56 4.2.1 Overall modeling considerations The Apps are based on the FE solver Abaqus and the usual transient heat transfer and quasi-static stress divergence balance equations are leveraged. For brevity they will not be reviewed here but the details are outlined in the softwares User Manuals [1]. In this study, the laser heat source is considered as a distributed volumetric heat flux and modelled by Goldaks semi-ellipsoid model [37]. In the local x-y-z frame of the laser, where x-axis aligns with the direction of travel, the laser heat source is described as " !# 2P η (x + vx t)2 (y)2 (z)2 Q (x, y, z, t) = ¯ √ exp − + ¯2 + 2 (4.1) a ¯b¯ cπ π ¯2 a b c¯ where P is the nominal power, η is absorption coefficient and vx is the travel speed ¯, ¯b, and c¯ are the dimensions of the ellipsoid along of the laser; t is the time; and a x-, y- and z-axis, respectively. The dimensions can be defined by laser parameters, including radius r, eccentricity e and penetration depth d, as r ¯ = er, ¯b = , c¯ = d a (4.2) e Thermal radiation qrad = εσ (T − Tz )4 − (T∞ − Tz )4  (4.3) and surface convection qconv = h (T − T∞ ) (4.4) of free surfaces during manufacturing processes are considered, where  is the material emissivity, h is convective coefficient, σ = 5.67 × 10−8 W/(m2 · ◦ C4 ) is the Stefan- Boltzmann constant, Tz = −273.15 ◦ C is the absolute zero temperature, and T∞ is the ambient temperature. 57 For the stress analysis, the temperatures are read from the heat transfer analysis. Similar techniques as described above are used to add material during the stress analysis. As for the material properties, temperature-dependent material models for conduc- tivity, specific heat, latent heat for melting and vaporization [31], thermal expansion, and elasticity from were chosen as shown in Figure 4.1. Figure 4.1: Ti-6Al-4V thermal properties: conductivity and specific heat capacity. Latent heat of fusion 2.86 × 105 J/kg Solidus temperature 1604 ◦ C Liquidus temperature 1650 ◦ C Latent heat of vaporization 9.83 × 106 J/kg Liquidus temperature 3290 ◦ C Vaporized temperature 3390 ◦ C Table 4.1: Latent heat of Ti-6Al-4V. For the purpose of this study a Ramberg-Osgood plasticity model with isotropic strain hardening was also used for the mechanical response. This model and its calibration is discussed extensively in Section 4.5. 58 Density ρ 4420kg/m3 Emissivity  0.25 2 ◦ Convection Coefficient h 18W/(m · C) Table 4.2: Other thermal properties of Ti-6Al-4V. 4.2.2 AM predictive simulations: spatial and time multiscale considerations The gigantic multiscale problem in both space and time associated with additive manufacturing process simulation is well discussed in the literature [79, 55]. The multi-physics relevant for sub-millimeter melt pools combined with local cooling rates of thousands of degrees Celsius per second contrast with the typical dimensions of parts being printed accompanied with printing times of hours/days. When it comes to predicting far field temperature evolutions (away from the areas where active melting/fusion occurs), during print and post removal from build plate distortions, as well as for residual stress predictions, averaged/lumped techniques of various degrees can be employed with reasonable reliability [106]. Macroscale phenomenological temperature-dependent elasto-plastic constitutive models appear to be reliable and computationally efficient for effective simulations when it comes to predicting the distortions and stresses in printed metal parts. Examples are given in the sections below. However, in order to capture the rapidly evolving temperature evolutions and high temperature gradients associated with the melt pools and their vicinity, special care must be exercised for both the meshes being used as well as the time marching strategies. A coarse mesh is employed in regions away from the zone of interest while a fine mesh with a 60µm characteristic dimension was used in the areas where metallurgical phase transformations computations are performed. After a mesh convergence study, it was found that this is the coarsest mesh (and hence most economical) that allows 59 for accurate computation of thermal gradients. In addition, a coarse time stepping scheme is adopted in regions away from the zone of interest or for most of the time associated with the recoating process when temperature evolution rates are low. By contrast, the time stepping in the regions of interest during laser scanning is adopted to be close to 1ms. After a time stepping convergence study, it was found that this is the coarsest time step that allows for accurate computation of heating and cooling rates associated with the vicinity of the melt pools. 4.3 An FE-based generic metallurgical phase transformation Framework Xie et al. [106] implemented a phenomenological metallurgical phase transformation model for Ti-6Al-4V alloy [26] within the general simulation framework of additive manufacturing processes, and validated the predicted porosity (unfused areas) and phase content against experiment measurements. In this chapter, the model for Ti- 6Al-4V is extended to a generic framework applicable to arbitrary metal alloys. First, the framework considers physical state changes from a stock feed (i.e. a raw material, such as a powder) via melting, vaporization, and solidification. During an additive manufacturing process, raw (unfused) materials are melted to a liquid state by an intense moving heat source, and then undergoes cooling and solidification when the heat source is removed. The solidified (fused) material can be re-melting when the heat source passes by a nearby location. As shown in Figure 4.2, transitions between raw, liquid, and solid states of an alloy material are considered as govern simply by temperature and occurring in a temperature range bounded by the solidus and liquidus temperature of the alloy. Thus, melt pool size and unfused zones can be predicted. Thermal and mechanical properties as a function of the physical state of the materials can also be considered in FE simulations. It should be noted that 60 this method does not consider porosity that might arise due to melt pool turbulence and the entrapment (or release) of gases during processing. To approximate such phenomena, semi-empirical models have been used for casting simulations [93, 56] or computational fluid dynamics studies for laser welding (with applicability to laser- based AM) [111] have previously been used. Secondly, the framework accounts for metallurgical solid-state phase transformations induced by heating and cooling events associated with rapid manufacturing processes and/or slower heat treatment applications. Modeling of metallurgical solid phase transformations is discussed below. The onset and the speed of transformations are associated with the alloys time-temperature-transformation (TTT) diagram that can be measured experimentally or generated by commercial software such as JMatPro and ThermoCalc. As a result, distribution and evolutions of the solid-phase com- position that further influences thermal and mechanical behavior can be assessed. The framework is based on phenomenological (internal state variable) metallurgical phase transformation models and, therefore, it is not a replacement for JMatPro or ThermoCalc or discrete-level methods. Figure 4.2: Diagram of physical state changes between raw/liquid/solid (RLS) states. 61 4.3.1 Kinetics of solid-phase transformations A solid-phase transformation can be written in a form of p1 + p2 + . . . + pnp → c1 + c2 + . . . + cnc (4.5) where np and nc are numbers of the parent and children phases, respectively; pi and ci denote the i th parent phase, and the i th child phases, respectively. The change of volume fraction of a parent phase during a transformation is assumed to be proportional to its current volume fraction. A similar assumption is made for children phases. With these assumptions, a transformation can be simplified as p→c (4.6) It consists of one master parent phase and one master child phase. The volume fractions of the master phases are Pnp Pnc fp = i=1 fp i fc = i=1 f ci (4.7) Two types of transformations are considered, diffusional transformation and marten- sitic (non-diffusional) transformation. The kinetics of a diffusional transformation under isothermal conditions can be described by the Johnson-Mehl-Avrami (JMA) model [26]. fc = ftot (1 − exp [−ktn ]) where ftot = fp + fc is the total volume frac- tion of all parent and children phases participated in the transformation; k and n are temperature-dependent coefficients that can be calibrated from TTT diagrams. Figure 4.3 shows an example TTT diagram that plots a diffusional transformation (D) and a martensitic transformation (M). A transformation is usually represented by at least two curves in a TTT diagram: a start curve (e.g., Ds , Ms in Figure 4.3) that marks the beginning of the transformation and a finish curve (e.g., Df , Mf in 62 Figure 4.3) where the transformation is completed. A curve in a TTT diagram rep- resents the time (as a function of temperature) that a transformation has completed a certain percentage. For instance, the start curve of the diffusional transformation Ds can be represented by the time, tD s (T ), when the transformation begins and the corresponding volume fraction of the children phases, fsD . Thus, the JMA coefficients of a diffusional transformation can be calibrated from TTT curves    log log 1 − fsD / log 1 − ffD n (T ) =  D  (4.8) log tD s (T ) /tf (T )  log 1 − fsD k (T ) = − n(T ) (4.9) (tD s (T )) Figure 4.3: An example time-temperature-transformation (TTT) diagram. Additive manufacturing processes involve rapid heating and cooling events that are considered as non-isothermal conditions. Therefore, an extended JMA model [26] that considers continuous temperature variation as a series of small consecutive isothermal steps is used fp (T + ∆T ) = ftot 1 − (1 − exp [−k (ξ + ∆t)n ]) 1 − fp,0 eq eq  , if fp > ftot fp,0 (4.10) 63 eq where ∆t is the time step size; ∆T is the temperature change since the last step; fp,0 are temperature-dependent equilibrium volume fractions of all parent phases in the mixture that only contains parent and children phases of the transformation; ξ is the fictitious time which would provide the same amount of transformation, resulted in volume fraction, fp (T ), at a constant temperature T + ∆T , defined as " eq !# n1 1 fp (T ) − ftot fp,0 ξ = − ln eq  (4.11) k ftot 1 − fp,0 If the diffusional transformation is reversible (p ↔ c), the reverse transformation c → p can be described as fp (T + ∆T ) = ftot (1 − exp [−k (ξ + ∆t)n ]) fp,0 eq , eq if fp < ftot fp,0 (4.12) Above all, the logic to model a diffusional transformation is summarized in Figure 4.4. Figure 4.4: The algorithm to model a diffusional transformation. 64 A martensitic transformation solely depends on temperature, as described by the Koistinen-Marburger (KM) model [26], fc (T ) = fc (To ) + (fp (To ) − fpr ) [1 − exp [−γ (Ms − T )] (4.13) where To is the temperature at the beginning of the current cooling cycle; the marten- sitic start temperature Ms and the transformation coefficient γcan be calibrated from TTT diagrams. Since a martensitic transformation is independent of time, this type of transformation is usually represented by at least two horizontal lines in the TTT diagram, as shown in Figure 4.5. For example, the start curve of the martensitic transformation can be represented by the temperature Ms where the transformation begins and the corresponding volume fraction of the children phases fsM . Thus, the transformation coefficient γ can be computed as    log 1 − fsM / 1 − ffM γ=− (4.14) Mf − Ms Under arbitrary thermal histories, a backward Euler integration of the differential form of the original KM model [106] is used fc (T ) − γ∆T (fp (T0 ) − fpr + fp (T0 )) fc (T + ∆T ) = (4.15) 1 − γdT The transformation occurs if the initial volume fraction of parent phases is higher than the retained amount of parent phases, namely, fp (T0 ) > fpr . The algorithm to model a martensitic transformation is shown in Figure 4.5. After the volume fraction of the master phases of a transformation are computed, the volume fraction of parent and children phases can be updated as fpi (T ) fci (T ) fpi (T + ∆T ) = fp (T + ∆T ) fp (T ) , fci (T + ∆T ) = fc (T + ∆T ) fc (T ) (4.16) 65 Figure 4.5: The algorithm to model a martensitic transformation. 4.3.2 Modeling framework The generic metallurgical phase transformation framework is implemented via user subroutine USDFLD that can be called by FE solver Abaqus/Standard at each inte- gration point in each increment during simulations. Results, including physical state of the material and volume fractions of solid phases, are stored as field variables for visualization and can be further used to evaluate thermal and mechanical properties as a function of those quantities. In each increment, the physical state of the material is first assessed. If the material is in the solid state, metallurgical phase transforma- tions are evaluated by following the general algorithm shown in Figure 4.6. The computation loops over all possible solid-phase transformations. The transforma- tions satisfying the condition of temperature and/or temperature rate for occurring will be computed and the volume fractions of its parent and children phases will be updated. Two types of transformations, diffusional and martensitic transformation, 66 are modeled in separate subroutines. In this framework, one can model any number of phases and transformations as desired. Priority of transformations is considered as the same as the input sequence of transformations if any two transformations can be concurrent. Future work will include considerations of the coupling of concurrent transformations. Figure 4.6: General algorithm to model the generic metallurgical phase transforma- tion framework. The framework allows the user to systematically define and to model phase transfor- mations in an alloy during arbitrary additive manufacturing process. Modeling pa- rameters are associated with the alloys phase diagram and TTT diagram, including solidus and liquidus temperature, number, and name of solid phases and solid-phase transformations, phase fraction upon solidification; and for each transformation, num- bers and names of parent and children phases, the transformation type, reversibility, temperature conditions when transformations can occur, curves in TTT diagrams for calibrating the JMA and KM models or calibrated model coefficients, etc. To- gether with the general simulation framework of additive manufacturing processes, the framework can calibrate the JMA and KM models and predicts phase transforma- 67 tions during print and heat treatment. Influences of the microstructure composition on thermal and mechanical material behaviors can also be considered in simulations. This is particularly useful, for example, when considering how microstructural fea- tures evolve as a function of the number of times that the powder has been recycled. In the case of Ti-6Al-4V, it is known that oxygen uptake during multiple powder uses leads to a loss of ductility after multiple re-uses [95]. In this case, the results of Malinov et al. [71] (an artificial neural network that can generate Ti-6Al-4V TTT diagrams for different oxygen content) could be used in conjunction with the general metallurgical modelling framework to simulate the effect of oxygen concentration on microstructural evolutions and therefore on mechanical properties. 4.4 Demonstration of the framework for two dif- ferent alloys The generic framework described above has been implemented in the commercial soft- ware ABAUS/Standard as a user subroutine. The framework enables users to define metallurgical phase transformations in any metal alloys through parameter/property tables, and to run finite element simulations of additive manufacturing processes. Two examples are given below with basic validation in one element tests. It should be noted that although the majority of the metallurgical phase transformations can be described by JMA/KM equations, some other models are also employed to pre- dict phase transformations, such as a parabolic model. The generic framework allows users to define those transformation models through minimal changes in the user sub- routine and input files. Moreover, one may also add features to predict non-generic microstructure morphologies as discussed in the Ti-6Al-4V simulation. 68 4.4.1 Ti-6Al-4V Titanium alloys and Ti-6Al-4V in particular have been key metals for the aircraft industry because of the combination of excellent mechanical properties with light weight [67]. These properties are very sensitive to the microstructure of the material, and therefore understanding and control of the microstructure evolution is of great importance in all titanium manufacturing processes [74]. During solidification from the liquid phase, Ti-6Al-4V first solidifies with a body centered cubic (β) crystal structure giving rise to a columnar beta grain structure. During slow cooling from the β phase region these grains partially transform to the hexagonally close packed (α) phase leaving an equilibrium microstructure containing a mixture of both α and β material. The morphology and distribution of the α phase which forms during this process is strongly dependent upon the cooling rate and processing conditions; During slow cooling the β transforms to a Widmanst¨atten microstructure characterized by laths of α, even slower cooling or annealing close to the beta transus can lead to the formation of a more equiaxed alpha phase which can form on β grain boundaries. Very high cooling rates, such as those associated with the SLM process, lead to the formation of non-equilibrium metastable martensitic phases the most common of which, α0 , has a hexagonal crystal structure similar to α phase but forms by a diffusionless (Martensitic) transformation. α0 has two morphologies an acicular lath morphology similar to Widmanst¨atten α and massive α0 which differs primarily with respect to the crystallographic orientations of groups of α0 laths. Both α0 morphologies are very often difficult to distinguish from the Widmanst¨atten α phase. In SLM processes, the high cooling rates leads to the almost complete transformation of the columnar β grains to α0 [107] and subsequent heat treatment, below the β transus temperature, is required to the decomposition the α0 phase into a mixture of the equilibrium α and β phases. 69 A phenomenological metallurgical phase transformation model for Ti-6Al-4V alloy was implemented to perform thermal and mechanical analysis of additive manufac- turing process by Xie et al. [106]. In this thesis, the same thermal properties, such as conductivity, latent heat and convection coefficient, are adopted. However, a more comprehensive solid phase transformation model [74] is implemented to improve the quality of mechanical property predictions. The α0 -phase in the previous model is replaced by αm -phase which consists of martensite α and massive α. The previ- ous α-phase is further divided into the equiaxed grain boundary α-phase, αgb , and Widmanst¨atten α-phase, αW [75]. Five metallurgical phase transformations between the hexagonal close-packed (hcp) αgb , αW , αm and the body-centered cubic (bcc) β-phase have been incorporated in the generic parent-children framework. During solidification, the material is consid- ered to be fully β-phase. Due to the subsequent thermal history, one martensitic transformation (β → αm ) and four diffusional transformation (β → αgb , β → αw , αm → β + αw , αgb + αw → β) are considered for calculation. The martensitic phase transformation occurs when the cooling rate is larger than 410 ◦ C/s and the massive phase transformation occurs when cooling rate is between 20 ◦ C/s and 410 ◦ C/s. Since the two diffusionless formations have similar chemical compositions and crystal structure, they are assigned one state variable (fαm ) in the model [2]. The overall transformation (β → αm ) is described by the KM model, with the transformation start temperature Ms = 851 ◦ C and γ = 0.005 ◦ C−1 . The αm - phase has a non-stable equilibrium crystal structure and can further transform to αw and β when the sample is annealed at medium high temperature. The transformation (αm → β + αw ) is described by the JMA model where the coefficients are given in Table 4.3 [74]. Fraction change of the children phase is proportional to its current equilibrium phase fraction. 70 Temperature [ ◦ C] k1 [s−n2 ] n1 [−] 400 0.0192 0.667 500 0.0147 1.106 700 0.0246 1.252 800 0.0307 1.326 Table 4.3: JMA coefficients for the transformation (α0 → β + αW ). Figure 4.7: Equilibrium β phase fraction. At a slower cooling rate from β phase, the alloy forms grain boundary α-phase and Widmanst¨atten α-phase via diffusion controlled transformations. These transforma- tions (β → αgb , β → αw ) are described by JMA model with two different sets of coefficients n and k. The coefficients are calibrated from the TTT-diagram as shown in Figure 4.8 [50]. Although the same JMA model is also applicable for the α → β transformation, a parabolic growth model [75] is chosen to describe the dissolution of α to β phase (αgb + αw → β) for the purpose of additive manufacturing process 71 simulation.  √ 1 − f eq (t + t)f 0 < t + t∗ < tcrit   β diss (T ) t + t∗ fαw +αgb (t + ∆t) = (4.17) 1 − fβeq (t + t) = fαeqw +α (t + t) ∗   gb t + t > tcrit where fdiss (T ) = 2.2 × 10−31 T 9.89 with T in kelvin [50, 51], tcrit = (fdiss (T ))−2 ,  2 fβ (t) t∗ = f eq (t+∆t)f diss (T ) and the equilibrium phase fraction fβeq is shown in Figure 4.7. β It is assumed that αw transforms to β first and αgb is transformed only when the αw phase has been extinguished [75]. Figure 4.8: TTT-diagram for diffusional phase transformations β → αgb and β → αW in Ti-6Al-4V. Two simple one element tests were performed to validate the metallurgical model. For the first analysis, two heating-cooling cycles between 26 ◦ C and 1826 ◦ C are assigned. As shown in Figure 4.9, the alloy is fully β-phase upon solidification. During rapid cooling, the β-phase transforms to αm and the diffusional transformations occur slowly upon heating. It is noted that, different from the result in previous model [106], the martensitic and massive phase transformation also occurs at a medium rapid cooling 72 rate ( 20 ◦ C/s to 410 ◦ C/s) and the transformation speed is a little slower. For the second analysis, the applied temperature history is taken from a node in the transient heat transfer analysis of the experiment-based column printing model (as discussed in Section 5.3). Printing process and heat treatment are both included, as shown in Figure 4.10. In the printing process, the αm -phase is formed during rapid cooling process from the highest temperature. In the subsequent slow cooling process, the diffusional phase transformation occurred is negligible because of extremely low transformation speed at low temperature. As a result, 78% of αm is predicted at the end of printing. In the heat treating process, αm transforms to (αw + β) through diffusional transformation and the remaining β phase from solidification also partially transforms to αw . The phase fraction of αw is predicted to be 90% at the end of heat treatment. Figure 4.9: Ti-6Al-4V metallurgical phase transformations and α-lathe width pre- diction one element analysis 1. 73 (a) 3D printing process (b) Heat treatment process Figure 4.10: Ti-6Al-4V metallurgical phase transformation one element analysis 2. 4.4.1.1 Microstructure morphology for Ti-6Al-4V The mechanical properties of Ti-6Al-4V are associated with microstructures and com- position of solid phases. The β grain structure is known to play a significant effect on fatigue performance [29, 106], and the Widmanst¨atten α lath width has been identified to have the largest effect on the strength properties. Therefore, it is impor- tant to calculate the α lathe width, β grain size and morphology in the FE analysis. These features are added to the metallurgical-microstructural-mechanical model with minimal changes in the generic user subroutine. β Grain Size and Grain Morphology Following Xie et al. [106], the (solidifica- tion) β Grain Size is updated by a parabolic growth law [57] when the temerature is between the β-transus temperature and the solidus temperature.   Q dnt+∆t = dnt + ∆t k exp − (4.18) R (T − TZ ) where d is the grain size at time t; d0 = 0 is the initial grain size, n = 2, k = 2.2m2 /s are coefficients, and Q = 251kJ/mol is an activation energy of β-phase of Ti-6Al-4V; R = 8.314J/(mol · ◦ C) is the universal gas constant. The β grain morphology is predicted by an experimentally measured G-R lookup map as shown in Figure 4.11 [57]. Given the magnitude of thermal gradient G and 74 T˙ solidification rate R = G , the β morphology can be determined to be equiaxed, columnar or their mixture. Figure 4.11: The solidification map for β grain morphology prediction. α Lath Width The widmanst¨atten α lathe width is assumed to be determined by temperature by an empirical Arrhenius equation as [75], teq lath = klath e −Rlath /T . (4.19) And the incremental calculation of the lath size is given as [75], fαw (t) tlath (t + ∆t) = (tlath − teq lath ) × + teq lath . (4.20) fαw (t + ∆t) The initial condition t0lath is set to be 1µm. According to Murgau et al. [74], the parameters are first chosen to be Arrhenius prefactor klath = 1.42µm and activation temperature Rlath = 294K. However, the suggest values seem to be unable to pre- dict the expected lath width. The calibration of the model parameters is discussed extensively in Chapter 5. 75 4.4.2 Steel 5140 The generic framework was also applied to a simplified metallurgical model of Steel 5140 to validate. Steel 5140 is a low carbon, low alloy chromium, molybdenum, nickel case hardening steel and has good weldability. Thermal properties are given in Table 4.5 [80]. Apart from the raw/liquid/solid material states, five metallurgical phases, Austenite(A), Bainite(B), Ferrite(F), Pearlite(P) and Martensite(M), were included for calculation. The alloy was considered as fully austenite upon solidification. Eight solid phase transformations between Austenite and the other three phases were con- sidered as shown in Table 4.6. Latent heat of fusion 2.70 × 105 J/kg Solidus temperature 1750 ◦ C Liquidus temperature 1800 ◦ C Latent heat of vaporization 7.60 × 106 J/kg Liquidus temperature 2910 ◦ C Vaporized temperature 3010 ◦ C Table 4.4: Latent heat of Steel 5140. Solid Liquid 3 Density [kg/m ] 4420 7800 Conductivity [W/(m · ◦ C)] 45 35 Specific heat capacity [mJ/(ton · ◦ C)] 650 840 Table 4.5: Thermal properties of Steel 5140. Transformation Type Occurring Conditions 1 A→F Bs < T < A1 (Cooling) 2 A→P Diffusional Bs < T < A1 (Cooling) 3 A→B Ms < T < Bs (Cooling) 4 A→M Martensitic Mf < T < Ms (Cooling) 5 F→A 6 P→A Diffusional A1 < T < A3 (Heating) 7 B→A 8 M→A Table 4.6: Phase transformations considered in the metallurgical model of Steel 5140. 76 During cooling, the martensitic transformation (A→M) was described by KM model and the other transformations (A→B, A→F and A→P) were described by JMA model. The transformation temperature bounds and the JMA/KM coefficients were calibrated from the TTT diagram as shown in Figure 4.12 [88]. During heating, all the phase transformations were modeled by the JMA equation. Since the phase trans- formations occur very fast [80], coefficients were chosen such that the transformation speeds are ten times of the cooling cases. Figure 4.12: Steel 5140 metallurgical phase transformations TTT-diagram. One element validation test was also conducted. As shown in Figure 4.13, the ma- terial undergoes heating-cooling cycles of different cooling rates. At slow cooling rate (1 ◦ C/s), the alloy transforms to mainly Ferrite and Pearlite; at medium cool- ing rate (10 ◦ C/s), a small amount of Banite phase was predicted; at high cooling rate (100 ◦ C/s), the martensitic phase transformation occurs and only negligible dif- 77 fusional formations were predicted. For the purpose of validation test, no quantita- tive conclusions are drawn here. However, with a more precise metallurgical model (TTT-diagram), the generic subroutine is expected to be applied to the simulation of arbitrary temperature histories in manufacturing process. (a) Cooling Rate = 1 ◦ C/s (b) Cooling Rate = 10 ◦ C/s (c) Cooling Rate = 100 ◦ C/s Figure 4.13: Steel metallurgical phase transformations. 4.5 Mechanical properties prediction for Ti-6Al- 4V The microstructure of Ti-6Al-4V is sensitive to the thermal history experienced, and different microstructures can be obtained under controlled thermal profiles [104, 65]. These microstructures are often quite unlike those of traditional wrought and cast materials and because the mechanical properties of additively manufactured parts depend on the underlying microstructure, the properties of SLM Ti-6Al-4V parts also have properties that are significantly different to conventionally cast or wrought products. For example, the yield and tensile strengths in laser-based AM components are gen- erally higher than those in annealed material and in the same range as age-hardened Ti-6Al-4V. This is likely due to the presence of a fine-grained microstructure in laser- based AM components [103]. On the other hand, the ductility of SLM Ti-6Al-4V components tends to range from 6% to 11% (depending on orientation with respect 78 to the build direction) which is lower than the 12 − 17% elongations observed in wrought conditions [20]. 4.5.1 An empirically-driven plasticity model for Ti-6Al-4V alloys Whilst there has been recent progress in applying CALPHAD-informed thermal mod- elling for cyclic thermal heating; modified phase transformation modelling considering superheating and supercooling effects on microstructure; and grain texture modelling based on melt-pool solidification [83], these approaches require significant compu- tational expense, expansive experimental validation exercises, and are still at low Technology Readiness Levels (e.g., not suitable for industrial applications). The ben- efit of these approaches is that unlike conventional JMAK or internal state variable approaches, phase field methods provide quantitative, rather than qualitative pre- dictions of the microstructure. In an attempt to try to leverage the computational benefits of internal state variable approaches whilst achieving the accuracy of phase field methods, the approach taken in this work was to generate a map from microstruc- tural features to mechanical properties that makes use of targeted mechanical testing, metallographic examinations, and the unique characteristics of Ti-6Al-4V SLM parts. The result is an empirically-driven plasticity model. Figure 4.14 illustrates the approach taken in this work to define an empirically-driven plasticity model for SLM Ti-6Al-4V parts 79 Figure 4.14: Diagram illustrating the approach to the empirically-driven plasticity model for Ti-6Al-4V. 4.5.2 Ramberg-Osgood plasticity model The stress-strain behavior at each material point in the finite element model was represented by a Ramberg-Osgood equation of the form:  n σ Ar σY σ ε= + (4.21) E E σY where  is the true total strain; σ is the true stress; E is the elastic modulus, Ar is a coefficient, σY is the yield strength at a strain level related to Ar , and n is the strain hardening exponent. Ti-6Al-4V exhibits continuous yielding behaviour with limited strain hardening and therefore a Ramberg-Osgood generally provides a reasonable fit to experimental stress-strain curves of Ti-6Al-4V samples. In this work, the yield stress was assumed to be defined as the 0.2% offset yield stress occurring at a strain of Y = 0.002 + σY /E, and EY Ar = −1 (4.22) σY log [(Eεu /σu ) − 1] − log [(EεY /σY ) − 1] n= +1 (4.23) log (σu /σY ) 80 where σu is the ultimate tensile strength. The elastic modulus was defined using a linear mixture rule based on the volume fractions of α, α0 , and β at each integration point. α and α0 have similar elastic moduli, higher than the elastic modulus of the β phase, with E = 117, 114 and 82 GPa, respectively [27]. 4.5.3 Microstructure and plastic yield behavior The microstructures of SLM builds are dominated by martensitic α0 and αw and so the yield strength of the material is assumed to be largely determined by the strengths of these phases. The yield strength of the αw phase can be defined using a Hall-Petch relationship of the form: K σY = σ0 + √ (4.24) d where σ0 is the intrinsic strength term, K is the strengthening coefficient, and, for the purposes of this work, d is the thickness of αw lamella. The Hall-Petch relationship is a simplistic and empirical relationship which traditionally approximates the yield strength materials as a function of grain size (d). An augmented Hall-Petch relation- ship, with grain size substituted for lamellar width, αw , has been successfully applied to Ti-6Al-4V SLM, electron beam melting and direct energy deposition materials, providing a good fit to Vickers hardness and yield strength measurements in heat treated materials, in which the as built α0 microstructure has transformed primarily αw lamellar [33, 3, 73, 40]. Measurement of αw lath width is relatively straightfor- ward using a linear intercept method and hence the availability of literature data to validate the applicability of the relationship between lath width and yield strength however, this relationship is more difficult to apply directly to as-built SLM materi- als which are primarily composed of martensitic α0 which forms over a range of very fine length scales, of 1µm and below, meaning defining an average α0 lath width is 81 more challenging in as-built material. This makes the augmented Hall-Petch rela- tionship more difficult to validate for as-built components. Nevertheless, a review of published data related to Hall-Petch coeffficients was undertaken and compared with nano-indentation hardness measurements (using a Berkovitch indenter) performed on Ti-6Al-4V SLM samples built for this work and suitable Hall-Petch coefficients were determined. Finally, with regards to the ultimate tensile strength and ductility (assumed for the purposes of the Ramberg-Osgood relationship to be the true total strain at the ulti- mate tensile strength), σu was defined as a linear function of the yield strength [72], and u was defined as a hyperbolic tangent function fit to published data [3, 49, 72]. The plasticity model was implemented as an ABAQUS user subroutine and will be presented detailly in the next chapter. 4.5.4 As-printed microstructure The definitions of the parameters above principally depend on the phase constitu- tion (volume fractions of α, α0 and β) and the α-lathe width. The generic phase transformation framework described in Section 4.3 enables the prediction of phase constitution, independent of the thermal processing conditions analyzed; that is, the methods presented in this chapter are applicable for other additive processes such as direct energy deposition and electron beam melting. Therefore, all that is required to define the Ramberg-Osgood parameters is a method for predicting the α-lathe width in the as-built and heat-treated condition. However, this is challenging because the as-built microstructure of SLM parts is pre- dominantly comprised of martensitic α0 . Presently, phase field methods are one of the few approaches that can be used to simulate martensitic transformations [100]; however, these methods are presently not well-suited to address computational do- 82 mains on the millimeter scale. Instead, an empirical solution was implemented for the current work. Based on the observation of an entirely martensitic microstructure, it was determined that a pragmatic approach would be to initialize the α-lathe widths in the as-built condition based on physical measurements of lathe widths from the beam and cube specimens. For a subset of samples, manual and semi-automatic measurements of alpha lathes were undertaken (using image segmentation software), and a histogram of alpha lathe widths was created. This histogram was then used to initialize the lathe widths in the finite element model after SLM processing by using it as a probability distribution function. Predicting grain nucleation and grain growth during 3D printing in PBF processes is a field of rather intense research activity. Methods such as phase field and cellular automata are used rather extensively in making numerical assessments of preferred grain morphology such as growth orientations, grain aspect ratios, etc.. All meth- ods require accurate resolution of the temperature evolutions (including temperature gradients) which often relies on a much finer time incrementation than the 1msec time stepping adopted in this work for resolving the temperature fields in regions of interest. Often, because of their expense, these methods are used for predicting grain morphology/growth in spatial regions that are very small (smaller than 1mm3 ). In this work we have not aimed at predicting grain growth in during the printing process primarily because of the expense associated with the computations in a rea- sonably large sample. Instead, we have adopted a more practical approach and we are initializing the thickness of as printed samples in a statistical fashion from im- age processing information from EBSD measurements. While the simulation predicts the temperature fields in the sample, at the end of the print simulation the software initializes the thickness of αlath according to this distribution in a spatially random fashion, for arbitrarily large samples. 83 4.5.5 Microstructure size evolutions during heat treatment It should be noted that the Hall-Petch relationship is being used in an empirical sense in this work. The relationship has proven useful for the purposes of this study, but it is only an approximation. The Hall-Petch relationship represents grain-boundary strengthening, but the coarsening of α-lathes during heat treatment (and therefore the softening of the yield stress through the application of a Hall-Petch relationship in this work) does not necessarily represent the actual mechanism that is resulting in the decrease in strength. Sallica-Leva et al. [87] argue that before α0 can decompose into α, the β phase (nano β particles) must first form on the α0 grain boundaries. This allows vanadium (and other β stabilisers) to diffuse into the β which in turn allows the α0 to transform into α. The authors note that the strength drop during this process is associated not only the coarsening of α lamella, but also with the diffusion of vanadium from α0 into β grains. Nevertheless, because of the difficulty in distinguishing α0 and α phases using analytical methods, it is likely that most published data on Hall-Petch relationships for additively manufactured Ti-6Al-4V parts include (or failure to distinguish) both α and α0 phases in the fitted parameters. 4.6 Conclusion While recent progress in improving the understanding of additive manufacturing pro- cesses, the process is often hindered by variable strength parts with inferior fatigue life when compared to traditional manufacturing processes such as forging and casting. The complex multiscale-multiphysics aspects of metallurgical phase transformations are at the very center of this problem. In this chapter we have introduced a general framework for predicting volume fractions of metallurgical phase transformations as an add-on to a previously introduced general FE-based framework for predicting tem- 84 perature evolutions, distortions, and residual stresses in printed parts. The modeling framework was demonstrated for two different alloys, Ti-6Al-4V and Steel 5140: • Computes detailed temperature evolutions and gradients associated with rapid heating and cooling. • Requires as input the collection of all relevant metallurgical phase transforma- tions that can take place, diffusional or non-diffusional. • Can be calibrated upfront or on-the-fly from temperature transformation dia- grams. • Predicts volume fractions of phases post printing or post heat treatment. Numerical models have little value unless they provide reasonable predictions when compared to physical test measurements. The focus of next chapter is to comprehen- sively validate the proposed framework. Nevertheless plenty of challenges remain as far as predictive models are concerned from a quantitative (not just qualitative) perspective. Most notably, the current phenomenological models for lamellar thickness predictions have been shown to be somewhat limited as far as predicting mechanical properties are concerned. Work is underway to assess other techniques (including cellular automata grain growth). 85 Chapter 5 Validation and Calibration of The Metallurgical Phase Transformation Framework in Ti-6Al-4V Alloy 5.1 Introduction In previous chapter, a customizable general simulation framework was documented for wide spectrum of additive manufacturing (AM) processes (laser and electron beam powder bed fabrication, direct energy deposition, arc welding, polymer extrusion, ink jetting) as implemented on the Dassault Systemes 3DX Platform based on new Finite Element (FE) technology implemented in Abaqus. The general simulation framework was complemented by a generic metallurgical phase transformation framework appli- cable to arbitrary metal alloys. Grain growth and grain aspect ratio assessment models are also included. 86 The quality of additively manufactured parts is often hindered by the inherent nature of the manufacturing process. It is not uncommon that the mechanical properties of AM parts are inferior to those manufactured via classical manufacturing techniques and therefore it would be desirable if numerical models for the prediction of mechan- ical properties of printed parts could be reliably. In this chapter we adopt existing empirical plasticity models [27, 40, 72, 49] for a Ti-6Al-4V alloy that connect yield strength/strain and ultimate strength/ductility limits to phase volume fraction con- tent and grain size information. The models utilize the predictions described above to estimate in the context of a Ramberg-Osgood plasticity model the hardening behavior of the as-printed and heat treated alloys. Experimental work is always necessary in order to gain insight into the physics gov- erning the processes being studied. Since numerical models are only valuable as long as they are reasonably predictive, a number of targeted experiments were conducted with the dual scope of understanding solidification behavior and to verify and validate the numerical predictions. We performed a number of validation exercises relying on a minimal amount of calibration activity. First, the overall thermal balance of the printing process was studied by recording the temperature evolutions at thermocou- ples placed in the build plate while the printing process is ongoing. After gaining con- fidence that the thermal evolutions (heating and cooling rates and thermal gradients) were predicted with reasonable accuracy, experiments were conducted to measure via SEM, EBSD and XRD techniques the post-printing and post-heat treatment the volume fraction of phases present as well as thickness of the lamellar structures de- veloped. Lastly, printed dog bone specimens were utilized to analyze plastic behavior and to assess the predictiveness of the models for mechanical properties. The chapter is organized as follows. In section 5.2 the experimental work that was conducted is presented while in section 5.3 we assess the predictiveness of the numer- 87 ical models. In section 5.4 conclusions are drawn to outline the success but also the remaining challenges. 5.2 Experiments and measurements (All experiments were conducted by TWI. Ltd, UK.) To support the calibration and validation of the generic metallurgical phase trans- formation framework, several Ti-6Al-4V SLM samples were manufactured and char- acterized. For all experimental trials, process parameters, build and recoat times, chamber conditions and powder characteristics were recorded for direct translation into the modelling framework. 5.2.1 Column printing with in-situ thermocouple measure- ments Tripathy et al. [99] demonstrated the validation of residual stress predictions using the additive manufacturing simulation framework in Abaqus. However, part-level residual stress and distortion predictions can be accurate without capturing the local temperature transients that are required for accurate metallurgical phase transforma- tions modelling. Therefore, a set of tests were defined for the present work whereby local temperature transients could be measured experimentally. To complement these tests, a bespoke substrate was manufactured by drilling a series of blind holes parallel to the build surface as shown in figure 5.1. The holes were position 0.50mm below the build surface and Type-K thermocouples were placed at the end of each hole, with the tip of the thermocouple pointing and in contact with the build surface. Five nominally identical cubes with a square cross-section having side length 10mm and a total build height of 13mm were built on top of 5 different thermocouples. The cubes were placed sufficiently far from one another to ensure 88 Figure 5.1: Image of the thermocouple substrate (left) and an illustration (right) showing the location of the blind holes for thermocouple placement. thermal isolation in the thermocouple measurements. A Renishaw AM 250 SLM machine was used with the process parameters shown in table 5.1. Parameter Hatch Border Power (W) 200 160 Exposure time (µs) 70 30 Point distance (µs) 60 20 Hatch spacing (µm) 95 N/A Beam diameter (µm) 75 Layer thickness (µm) 60 Border distance offset of 40um with one boundary contour; powder enters the chamber at room temperature, but the substrate pre-heat Notes: temperature was 100oC. Scan strategy for the in-fill was parallel, alternating scan lines, with an alternation of 67o every layer. Table 5.1: Process parameters for the thermal validation exercise involving the ther- mocouple substrate. During SLM processing of the first 20 layers, the thermocouples recorded tempera- tures at 100Hz. The entire cross section (slice) of one square was built before process- ing of a subsequent cube began. The intention was to have five nominally identical temperature records that could be aligned using a time shift factor. However, due to the inability to accurately control the placement of the thermocouple tip, coupled 89 with the highly localized temperature transients during SLM processing of the first layers of the build, some differences between measurements could be observed. 5.2.2 Cruciform printing and tensile testing Having validated the heat transfer model using the thermocouple substrate, a new set of experiments was devised to both calibrate and validate the phase transformation model and subsequent structure-property map. Two cruciform specimens were built using the same powder as all previous exercises and the same machine. The cruciform had total build height 40mm, wall thickness 4mm, and length 140mm. The process parameters for the cruciform parts were similar to those of the cubes for the ther- mal validation except that four boundary contours were used and the layer-to-layer rotation of scan vectors was 90 ◦ . 5.2.2.1 Tensile specimen preparation Following SLM processing of the cruciform specimens, the parts were removed from the substrate. A diagram and photograph of the cruciform parts are shown in fig- ure 5.2. Of the two cruciform, CRUC1 was left in the as-built condition, and CRUC2 Figure 5.2: Diagram (left) showing the naming convention for the cruciform speci- mens and a photograph (right) of the as-built cruciform specimens. was heat treated. The heat treatment involved using molybdenum heating elements 90 under vacuum (10−5 mbar) condition with 10 ◦ C/min heating, holding at 750 ◦ C for two hours and cooling by convection for 24 hours (see figure 5.3). This heat treatment is similar to a typical mill annealing heat treatment applied to wrought products. Figure 5.3: Heat treatment of CRUC2. 5.2.2.2 Tensile testing of cruciform specimens From each cruciform, 8 flat, dogbone tensile specimens were machined: 4 in the vertical build direction, and 4 in the horizontal orientation as shown in figure 5.4. The specimen dimensions were in accordance with BS EN 2002-1 (Aerospace series metallic materials tensile testing) which has a 3mm wide gauge, 12mm gauge length, 15mm central parallel length, and grip ends that are 9mm wide. The photographs of the tensile specimens are shown in figure 5.4. In the as-built condition, the specimens are characterized by 45 ◦ failures with little necking, while in the heat treated condition more of ductile cup-cap failures with clear necking were observed. The stress-strain 91 Figure 5.4: Photographs of tensile specimens. curves from the cruciform tensile test specimens are shown in figure 5.5. The following observations can be made regarding the stress-strain curves from the cruciform parts: • Whilst SLM processing is known to introduce anisotropy and texture, the stress- strain curves show only small differences in the yield and hardening behavior for horizontal and vertical orientation specimens in the same condition (as-built or heat treated). 92 Figure 5.5: Stress-strain curves from the cruciform tensile test specimens. • The heat treated specimens exhibit a lower yield point and ultimate tensile strength compared with the as-built specimens, as expected. • The horizontal specimens, in general, exhibit larger fracture strains (and uni- form elongation limits) than the vertical specimens. 5.2.2.3 Metallographic examination of cruciform specimens Following tensile testing, the remaining plates from the cruciform parts (in addi- tion to the gripped and unstrained portions of the tensile specimens) were subjected to metallographic examination. As-deposited and heat treated specimens were, sec- tioned, metallographically prepared and examined by back-scatter scanning electron microscopy (SEM) in a Zeiss Sigma SEM. Qualitative chemical area maps were ob- tained using Oxford Instruments energy-dispersive X-ray (EDX) spectroscopy equip- ment in the SEM. Additionally, grain structure was quantified using EBSD analysis performed on 115 × 115µm areas using an accelerating voltage of 30kV and a step size of 0.05µm. Analysis was completed using Oxford Instruments Aztec software and grain reconstruction was achieved using a critical misorientation angle of 10 degrees 93 and grain completion to 2 degrees. The average width of laths was determined by dividing the feret length of a grain, by the grain area. Grain orientation EBSD maps from the as-deposited and heat treated specimens are shown in figure 5.6. Backscatter SEM images of the two specimens are shown in figure 5.7, which show a, qualitatively, larger lath width for the heat treated speci- men. Features showing as white in the backscatter image were evident for the heat treated specimen, but were absent in the as-deposited specimen. Brighter features, in backscatter images, are indicative of higher atomic number elements being present. EDX elemental mapping was conducted on the heat treated specimen to determine the chemical makeup of the white features. As can be seen from figure 5.8, the white features were rich in V and deficient in Al. (a) As-deposited (b) Heat-treated Figure 5.6: Grain orientation maps from heat treated and as-deposited specimens. (a) As-deposited (b) Heat-treated Figure 5.7: Back-scatter SEM images of heat treated and as-deposited specimens. 94 Figure 5.8: Backscatter SEM images and elemental maps for the heat treated speci- men: Backscatter image (left); V chemical map (right). 5.2.2.4 X-ray diffraction analysis X-ray diffraction spectra in Bragg-Brentano, i.e. symmetric, configuration were ob- tained from as deposited and heat treated specimens, using a Bruker D8 Advance, equipped with a Cu tube X-ray source. Diffraction patterns were obtained between 30 and 130 degrees in θ. The purpose of the analysis was to determine differences in the phases present in the as-deposited and heat treated specimens. Polished metal- lographic sections, used previously for metallographic analysis were illuminated with X-rays. Strain free lattice parameter, d0 , measurements, for the α/ α0 phases, of the as- deposited and heat-treated specimens, were conducted in order to determine if peak shifts were related to changes in the residual stress state or related to the chemistry of the phases. These features are convoluted when using Bragg-Brentano geometry. The strain-free lattice parameter can be determined from data generated in a stan- dard sin2 ψ surface residual stress measurement. Such a measurement involves the determination of the d-spacing for a particular diffraction peak at various tilt angles, ψ. A plot of d-spacing v.s. sin2 ψ should be linear with a gradient proportional to the residual stress, σ, i.e: ∂dψ 1 2 = S2hkl σd0 (5.1) ∂ sin ψ 2 95 d0 , can then be determined by interpolation of the d-spacing v.s. sin2 ψ plot to the strain free direction of the biaxial stress state: s ∗ −2S1hkl ψhkl = sin−1 1 hkl (5.2) S 2 2 The diffraction elastic constants S1hkl and 12 S2hkl were inferred from single crystal elastic constants for titanium using the Kroner model in the IsoDEC software package. d0 measurements were also obtained by illuminating a polished metallographic specimen with X-rays. It is important to note that although such sectioning of the component would alter the residual stress state, the strain-free lattice parameter would of course be unaffected. The remaining residual stress as well as d0 were determined from examination of the 11¯24. Figure 5.9: XRD patterns obtained in Bragg-Brentano configuration for as deposited and heat treated specimens. Black: as-deposited, Red: heat treated. XRD patterns in Bragg-Brentano configuration are shown in figure 5.9. Small β titanium peaks were evident in the heat treated specimen, but were absent in the as- deposited. Additionally there was a peak shift to lower 2θ values after heat treatment, as well as a narrowing of the α titanium peaks. 96 5.2.2.5 Discussion on EBSD and XRD analysis EBSD results indicated changes in lath width with heat treatment. This suggestion is reinforced by the backscatter images which showed an increased grain size for the heat treated specimen; however, it is challenging to quantify this for the present work as the indexing rate was quite low (70%) in the as-deposited specimen and therefore the statistics were skewed when comparing as-built to heat-treated measurements. The β titanium phase was not detected, in either of the specimens, through EBSD analysis. However, features with high levels of β stabilising elements, i.e. V, and lower than back-ground levels of α stabilising elements, i.e. Al, were detected through backscatter imaging and qualitative chemical analysis. The presence of β titanium was also confirmed by XRD analysis. The sin2 ψ measurements suggested that α0 had decomposed to α+β during heat treatment. After heat treatment, the α lattice parameters were much closer to what would be expected of Ti-6Al-4V [86]. Peak narrowing was also evident for the heat treated specimen. Such peak narrowing is indicative of a reduction in deformation, i.e. type III residual stress, and/or an increase in grain size. It is likely that both of these factors would have played a role in the peak narrowing evident in the heat treated specimen. 5.3 Calibration and validation for Ti-6Al-4V 5.3.1 Column printing and comparison with thermocouple data The Ti-6Al-4V column printing experiment was modeled with the generic frame- work using ABAQUS. Figure 5.10 shows the finite element model used in the simula- tions. Twenty layers (60µm per layer) of five 10mm × 10mm cubes were printed on a 97 226mm × 248mm × 15mm substrate. A uniform fine mesh of 100µm × 100µm × 60µm was used for one cube (referred as cube 5) as shown in figure 5.10, and a coarse mesh of 1mm × 1mm × 1mm was used for the other cubes (referred as cube 1-4). The substrate was also partitioned into a coarse mesh (10mm × 10mm × 10mm) part and a fine mesh part underneath cube 5 to perform the analysis. The laser beam profile Figure 5.10: FE model of column printing experiments. was described by the Goldak model with laser spot radius r = 50µm and penetra- tion depth h = 60µm. The laser power absorption coefficient was 55%. Ambient temperature was taken to be 26 ◦ C and the bottom surface of the substrate was kept at 100 ◦ C during the heat transfer analysis. The time increment was selected to be 0.4s. As shown in figure 5.11, the temperature histories at positions with several different depths underneath the cube 5 were predicted during first 20 layers print- ing. Although there are some differences between the finite element predictions of temperature transients and the thermocouple measurements, the overall agreement with regards to the decay of peak temperatures and heating rates are satisfactory. The differences in cooling rates are challenging to address due to the experimental errors incurred with respect to uncertainty in thermocouple placement and surround- ing powder bed temperature. Nevertheless, this exercise demonstrates the validation of the heat transfer model. 98 Figure 5.11: Column printing, thermo-couples temperature experiment measure- ments vs. FEA predictions at cube 5. top) first 20 layers printing, the simulation data was taken from the TC 0.0625mm below the cube; bottom) printing of the second layer. Blue solid line shows the thermocouple recorded temperatures, and all other lines represent finite element predictions at different depths beneath the cube, as shown in the table. The closer the position is to the cube, the higher the temperature is predicted. 5.3.2 Cruciform printing model 5.3.2.1 Model setup Following the experiments presented in section 5.2.2, the SLM processing of the cru- ciform specimens and the tensile tests on one dog bone sample are simulated within the additive manufacturing framework in Abaqus. Figure 5.12 shows the geometry and mesh of the finite element model for the heat transfer analysis. To simplify the manufacturing process, only the reduced section of the CRUC1-EAST-H1 dog bone specimen (figure 5.2) is modeled with fine mesh (60 × 58.8 × 60µm elements). The layers below this section and above the substrate are modeled with coarse mesh of element size 1.2 × 1 × 1.15mm. The overall size of the finely computed section is 99 6 × 4 × 0.24mm. This corresponds to four layers of laser melting. The same laser scan paths as in the experiments are generated for the simulation, where each layer con- tains 4 boundary contours, and hatch-type infill with 90 ◦ alternation angle between layers. Figure 5.12: Geometry (left) and mesh (right) of the SLM process on one dog bone specimen in the cruciform model After the SLM process is simulated in the heat transfer analysis, part of the finite element model, together with their local mechanical properties, are imported into the stress analysis for the prediction of tensile tests. To avoid the boundary effects resulted from the simplified heat transfer model, only a 3.96 × 2 × 0.12mm material block within the previous finely meshed part is imported into the new analysis. This material block represents the gage length section of the dog bone specimen in real experiments. 5.3.2.2 Metallurgical transformation predictions A snapshot of a temperature evolutions during the printing of the gage length section of the dog bone specimen is shown in figure 5.13. The detailed temperature histories are driving the phase transformations as described/calibrated above. Phase fraction volumes are predicted at each integration point as driven by the rapidly evolving temperatures, melting/solidification/remelting patterns. Volume fractions are then averaged for the entire specimen to compute an average volume fraction at the of the print and at the end of the heat treatment history, as shown in table 5.2. 100 Figure 5.13: Temperature evolutions: snapshot during a lasering sequence. Gray contourplots represent temperature as above melting. β α0 α As-built (Simulation) 0.2% 95.58% 4.212% As-built (Experiment) None observed 93.87% 6.13% (by morphology) Heat-treated (Simulation) 11.8% 0.1% 88.04% Heat-treated (Experiment) 4.3% Cannot distinguish between two (95.7%) Table 5.2: Volume fraction predictions 5.3.2.3 As-printed grain morphology Predicting grain nucleation and grain growth during 3D printing in PBF processes is a field of rather intense research activity. Methods such as phase field and cellular automata are used rather extensively in making numerical assessments of preferred grain morphology such as growth orientations, grain aspect ratios, etc.. All meth- ods require accurate resolution of the temperature evolutions (including temperature gradients) which often relies on a much finer time incrementation than the 1ms time stepping adopted in this work for resolving the temperature fields in regions of in- terest. Often, because of their expense, these methods are used for predicting grain morphology/growth in spatial regions that are very small (smaller than 1mm3 ). In this work we have not aimed at predicting grain growth in during the printing process primarily because of the expense associated with the computations in a rea- 101 sonably large sample. Instead, we have adopted a more practical approach and we are initializing the thickness of as printed samples in a statistical fashion from image processing information from EBSD measurements. Figure 5.14 shows the histogram used. While the simulation predicts the temperature fields in the sample, at the end of the print simulation the software initializes the thickness of αlath according to this distribution in a spatially random fashion, for arbitrarily large samples. Figure 5.14: Histogram of measured αlath thickness from image processing of EBSD data in several samples. 5.3.2.4 Grain thickness evolution during heat treatment It should be noted that the Hall-Petch relationship is being used in an empirical sense in this work. The relationship has proven useful for the purposes of this study, but it is only an approximation. The Hall-Petch relationship represents grain-boundary strengthening, but the coarsening of α-lathes during heat treatment (and therefore the softening of the yield stress through the application of a Hall-Petch relationship in this work) does not necessarily represent the actual mechanism that is resulting in the decrease in strength. Sallica-Leva et al. [87] argue that before α0 can decompose 102 into α, the β phase (nano β particles) must first form on the α0 grain boundaries. This allows vanadium (and other β stabilisers) to diffuse into the β which in turn allows the α0 to transform into α. The authors note that the strength drop during this process is associated not only the coarsening of α lamella, but also with the diffusion of vanadium from α0 into β grains. Nevertheless, because of the difficulty in distinguishing α0 and α phases using analytical methods, it is likely that most published data on Hall-Petch relationships for additively manufactured Ti-6Al-4V parts include (or failure to distinguish) both α and α0 phases in the fitted parameters. 5.3.2.5 Mechanical properties The plasticity model was implemented as an ABAQUS user subroutine. Given the phase fractions (fα , fα0 , fβ ) and α lath thickness tlath , mechanical properties were calculated as follows, • Youngs modulus: E = Eα × fα + Eα0 × fα0 + Eβ × fβ , σy • Yield stress and yield strain: σy = σ0 + √ K , y = 0.002 + , tlath E • Ultimate tensile stress: σu = 1.1046σy − 13.858, • Uniform elongation limit: u = 10.00 + 1.9 tanh (11.00tlath − 7.0), where Eα = 117GPa, Eα0 = 114GPa, Eβ = 82GPa [27, 40, 72, 49]. The yielding and ultimate points were converted into true stress and true strain, and a Ramberg- Osgood power-law equation was fitted from those true values.  n σ Ar σY σ ε= + (5.3) E E σY EY log [(Eεu /σu )−1]−log [(EεY /σY )−1] where Ar = σY − 1 and n = log (σu /σY ) + 1. 103 Virtual tensile tests of Ti-6Al-4V SLM parts were performed with the plasticity model. A minimal calibration in the empirical Hall-Petch coeffficients σ0 and K was made to match the experimental stress-strain curves. Rectangular sample The rectangular sample taken from the metallurgical- microstructrue prediction was imported to perform the mechanical analysis. The tensile tests of both as-built and heat-treated samples were performed. As shown in figure 5.15, a uniform mesh of 60µm × 60µm × 60µm was used in the model. The rectangular sample was elongated with the generalized plane stress condition. The left surface was fixed in X direction, and a displacement loading of 0.4 mm in X direction was applied to the right surface in 100 time increments. The top and bottom surfaces were kept flat during the analysis. Figure 5.15: FE model of the rectangular sample for virtual tensile tests. F F l−l0 The engineering stress σ = A0 = 2mm×0.12mm and the engineering strain  = l0 = l−3.96mm 3.96mm were calculated to calibrate with the experimental results. A final choice of σ0 = 691MPa and K = 420MPa · mm1/2 were taken and the corresponding stress- strain curves are shown in figure 5.16. 104 Figure 5.16: Stress-strain curves of the rectangular samples under uniaxial tension. Dogbone experiments The FEA stress-strain curves of rectangular samples agree well with the experimental results up to ultimate tensile limiting point, however, the post-ductility limit necking behavior was not captured. Therefore, some additional simulations were performed with the dog bone specimen. As shown in figure 5.17, a fine mesh of 60µm × 60µm × 60µm was used for the center part, and a coarse mesh of 0.5mm × 0.5mm × 0.5mm was used for the rest. (a) (b) Figure 5.17: Heat transfer analysis of dogbone specimen. a) FEA model; b) lath thickness contour plots of both as-built sample and heat treated sample. For the as printed specimen, metallurgical phase fractions were taken from simulation results for the rectangular sample in section 5.3.2.2. Similar as before, the α lathe width was initialized from the EBSD measurement in a statistical fashion. The cali- brated plasticity model was used for center part only and the side parts were modeled 105 with elastic relationship. Plane stress boundary condition was considered in the sim- ulation and a displacement loading in X direction was applied to the right surface. During the analysis, necking will occur in middle of the specimen at an elongation of 6%. The stress-strain curve of center area was plotted against the experimental results as shown in figure 5.18. To model the heat treated dog bone specimen, a heat transfer analysis was per- formed followed by the mechanical analysis. The heat transfer simulation predicted the metallurgical phase fractions and microstructure morphology from the as printed specimen. Then the same virtual tensile test was performed and the stress-strain curve of center area is also shown in figure 5.18. (a) (b) Figure 5.18: Tensile tests of dogbone specimens. a) necking behavior captured in FE simulations for both as-built (upper plot) and heat-treated (lower plot) samples; b) comparison of stress-strain curves between finite element predictions with experimen- tal results. Properties of both as-built sample and heat treated sample are included. 5.4 Conclusion The AM process simulation framework described in the previous chapter for predicting volume fractions of metallurgical phase transformations was validated for tempera- ture; melt pool sizes; unfused powders; grain morphology; and mechanical properties predictions. 106 Numerical models have little value unless they provide reasonable predictions when compared to physical test measurements. Temperature measurements and volume fractions metallographic assessments via EBSD/XRD techniques of printed and heat- treated samples were performed. Numerical predictions were compared against test measurement with reasonable success. The models were calibrated/validated against macro-scale dogbone tests. The dog- bones were cutout from 3D printed walls and comparisons for as-printed and post heat treatments were executed again with reasonable success. The models are also able to capture post-ductility limit necking behavior. The combined numerical and experimental methodology introduced in this chapter and the corresponding modelling framework paper provide a pragmatic approach to address process-structure-property-performance relationships in additive manufactur- ing. Whilst the proposed approach is not derived entirely from first principles (like emerging integrated computational materials engineering approaches), the approach is valuable in that it can be easily calibrated with a small number of experimental trials. Once the model is calibrated, the end user has access to powerful metallurgical modelling tools that can help inform process parameter optimization in reasonable computational run times. Therefore, whilst still a continuum-level internal state vari- able approach, the method can directly link processing conditions to microstructural features and then link these features to mechanical performance metrics which can be explored and optimized iteratively. Nevertheless plenty of challenges remain as far as predictive models are concerned from a quantitative (not just qualitative perspective). Most notably, the current phenomenological models for lamellar thickness predictions have been shown to be somewhat limited as far as predicting mechanical properties are concerned. Work is underway to assess other techniques (including cellular automata grain growth). 107 Chapter 6 Conclusion Employing theoretical derivations and finite element simulations, I have investigated the surface instabilities in generalized Blatz-Ko material and developed a metallurgi- cal phase transformation framework for modeling of additive manufacturing process. The scientific findings for each chapter are summarized as follows. In Chaper 2, a linear perturbation analysis was performed to investigate the effects of material compressibility on the surface wrinkling of generalized Blatz-Ko materials. For the wrinkling in bilayer systems, it is found that increased incompressibility in the substrate, as well as pre-stretch, slightly increases the critical wrinkling strain. In addition, this wrinkling type is suppressed for a subclass of generalized Blatz-Ko material, and this effect is maximized for the material with parameter f = 0. Our theoretical study shows that for certain ranges of parameters (f and ν), wrinkling is admissible in a slab of finite thickness but infinite horizontal extent (without a thin film). The wrinkling deformation is characterized by sinusoidal displacements in both horizontal and vertical directions. Moreover, the critical wrinkling strain is found to be higher than that of the wrinkling in a bilayer system and it is found to be more significantly influenced by the compressibility of material. Our study opens 108 up a promising prospect of designing stretchable electronics and controlling surface instabilities. In Chapter 3, I have explored the stability of the wrinkling in finite slabs though theoretical analysis and finite element simulations. I performed a second-order per- turbation analysis with the Koiter method and it was shown that the wrinkling of the form considered is stable. Extensive finite element simulations are also performed to validate and investigate the post-buckling behaviors. FE analyses demonstrate that the critical wrinkling strain is in agreement with the preceding derivations. A Fourier transformation analysis shows that the wrinkling deformations that emerge in the FE analyses follow the theoretically predicted sinusoidal functions and frequencies. The sinusoidal periodicity is found to be dependent on the initial perfection settings. In Chapter 4, we have introduced a general framework for predicting volume frac- tions of metallurgical phase transformations as an add-on to a previously introduced general FE-based framework for predicting temperature evolutions, distortions, and residual stresses in printed parts. The framework accounts for phase transforma- tions (physical state changes) from stock feed or raw material (e,g., powder) via melting/vaporization/solidification followed by metallurgical solid state phase trans- formations induced by either rapid heating or cooling events associated with typical 3D printing sequences but also with slower-rate temperature evolutions as associated with heat treatment applications. The modeling framework was demonstrated for two different alloys, Ti-6Al-4V and Steel 5140. In summary, the framework is applicable for modelling additive manufacturing process of arbitrary alloys and connect mechan- ical properties to phase volume fraction content and microstructure morphologies. In Chapter 5, several Ti-6Al-4V SLM samples were manufactured and characterized to validate and calibrate the generic metallurgical phase transformation framework. Temperature measurements and volume fractions metallographic assessments via SEM/EBSD/XRD techniques of printed and heat-treated samples were performed. 109 Numerical predictions were compared against test measurement with reasonable suc- cess. An empirical plasticity model was adopted to predict mechanical properties of the samples. The virtual tensile tests of dogbone specimens were performed and the models are able to predict the stress-strain curve up to ultimate strength and capture the post-ductility limit necking behavior. The framework has shown capability of producing reliable predictions of volume fractions, grain morphology, and mechanical properties in printing process and heat treatment process, and has seen applications in optimization of additively manufactured parts. 110 Bibliography [1] Abaqus 2018 user manual. https://www.3ds.com/products-services/ simulia/support/documentation/, 2018. [2] T Ahmed and HJ Rack. Phase transformations during cooling in α+ β titanium alloys. Materials Science and Engineering: A, 243(1-2):206–211, 1998. [3] Sinan Saadi Al-Bermani. An investigation into microstructure and microstruc- tural control of additive layer manufactured Ti-6Al-4V by electron beam melting. PhD thesis, University of Sheffield, 2011. [4] Anesia Auguste, Lihua Jin, Zhigang Suo, and Ryan C Hayward. The role of substrate pre-stretch in post-wrinkling bifurcations. Soft Matter, 10(34):6520– 6529, 2014. [5] Anesia Auguste, Lihua Jin, Zhigang Suo, and Ryan C Hayward. Post-wrinkle bifurcations in elastic bilayers with modest contrast in modulus. Extreme Me- chanics Letters, 11:30–36, 2017. [6] Melvin Avrami. Granulation, phase change, and microstructure kinetics of phase change. iii. The Journal of chemical physics, 9(2):177–184, 1941. [7] PV Bayly, LA Taber, and CD Kroenke. Mechanical forces in cerebral cortical folding: a review of measurements and models. Journal of the mechanical behavior of biomedical materials, 29:568–581, 2014. [8] Millard F Beatty. Topics in finite elasticity: hyperelasticity of rubber, elas- tomers, and biological tissueswith examples. Applied Mechanics Reviews, 40(12):1699–1734, 1987. [9] H Bikas, Panagiotis Stavropoulos, and George Chryssolouris. Additive manufac- turing methods and modelling approaches: a critical review. The International Journal of Advanced Manufacturing Technology, 83(1-4):389–405, 2016. [10] MA Biot. Surface instability of rubber in compression. Applied Scientific Re- search, Section A, 12(2):168–182, 1963. [11] Paul J Blatz and William L Ko. Application of finite elastic theory to the deformation of rubbery materials. Transactions of the Society of Rheology, 6(1):223–251, 1962. 111 [12] Ned Bowden, Scott Brittain, Anthony G Evans, John W Hutchinson, and George M Whitesides. Spontaneous formation of ordered structures in thin films of metals supported on an elastomeric polymer. Nature, 393(6681):146– 149, 1998. [13] D Brackett, I Ashcroft, and R Hague. Topology optimization for additive man- ufacturing. In Proceedings of the solid freeform fabrication symposium, Austin, TX, volume 1, pages 348–362. S, 2011. [14] Fabian Brau, Pascal Damman, Haim Diamant, and Thomas A Witten. Wrinkle to fold transition: influence of the substrate response. Soft Matter, 9(34):8177– 8186, 2013. [15] Silvia Budday, Charles Raybaud, and Ellen Kuhl. A mechanical model predicts morphological abnormalities in the developing human brain. Scientific reports, 4:5644, 2014. [16] Shengqiang Cai, Dayong Chen, Zhigang Suo, and Ryan C Hayward. Creasing instability of elastomer films. Soft Matter, 8(5):1301–1304, 2012. [17] Changyong Cao, Hon Fai Chan, Jianfeng Zang, Kam W Leong, and Xu- anhe Zhao. Harnessing localized ridges for high-aspect-ratio hierarchical pat- terns with dynamic tunability and multifunctionality. Advanced Materials, 26(11):1763–1770, 2014. [18] Yanping Cao and John W Hutchinson. From wrinkles to creases in elastomers: the instability and imperfection-sensitivity of wrinkling. In Proc. R. Soc. A, volume 468, pages 94–115. The Royal Society, 2012. [19] Yanping Cao and John W Hutchinson. Wrinkling phenomena in neo-hookean film/substrate bilayers. Journal of applied mechanics, 79(3):031019, 2012. [20] Beth E Carroll, Todd A Palmer, and Allison M Beese. Anisotropic tensile behavior of ti–6al–4v components fabricated with directed energy deposition additive manufacturing. Acta Materialia, 87:309–320, 2015. [21] Edwin P Chan, Erica J Smith, Ryan C Hayward, and Alfred J Crosby. Surface wrinkles for smart adhesion. Advanced Materials, 20(4):711–716, 2008. [22] Xi Chen and Jie Yin. Buckling patterns of thin films on curved compliant substrates with applications to morphogenesis and three-dimensional micro- fabrication. Soft Matter, 6(22):5667–5680, 2010. [23] Chen Chu, Greg Graf, and David W Rosen. Design for additive manufacturing of cellular structures. Computer-Aided Design and Applications, 5(5):686–696, 2008. 112 [24] Jun Young Chung, Adam J Nolte, and Christopher M Stafford. Surface wrin- kling: A versatile platform for measuring thin-film properties. Advanced Mate- rials, 23(3):349–368, 2011. [25] Peter C Collins, B Welk, T Searles, J Tiley, JC Russ, and HL Fraser. De- velopment of methods for the quantification of microstructural features in α+ β-processed α/β titanium alloys. Materials Science and Engineering: A, 508(1- 2):174–182, 2009. [26] Ant´onio Crespo. Modelling of heat transfer and phase transformations in the rapid manufacturing of titanium components. In Convection and Conduction Heat Transfer. InTech, 2011. [27] Antonio Crespo, Augusto Deus, and Rui Vilar. Modeling of phase transforma- tions and internal stresses in laser powder deposition. In XVII International Symposium on Gas Flow, Chemical Lasers, and High-Power Lasers, volume 7131, page 713120. International Society for Optics and Photonics, 2009. [28] Mazen Diab, Teng Zhang, Ruike Zhao, Huajian Gao, and Kyung-Suk Kim. Ruga mechanics of creasing: from instantaneous to setback creases. In Proceed- ings of the Royal Society of London A: Mathematical, Physical and Engineering Sciences, volume 469, page 20120753. The Royal Society, 2013. [29] Matthew J Donachie. Titanium: a technical guide. ASM international, 2000. [30] Kirill Efimenko, Mindaugas Rackaitis, Evangelos Manias, Ashkan Vaziri, L Ma- hadevan, and Jan Genzer. Nested self-similar wrinkling patterns in skins. Nature materials, 4(4):293–297, 2005. [31] Zhiqiang Fan and Frank Liou. Numerical modeling of the additive manufac- turing (am) processes of titanium alloy. In Titanium Alloys-Towards Achieving Enhanced Properties for Diversified Applications. InTech, 2012. [32] William E Frazier. Metal additive manufacturing: a review. Journal of Mate- rials Engineering and Performance, 23(6):1917–1928, 2014. [33] Haize Galarraga, Robert J Warren, Diana A Lados, Ryan R Dehoff, Michael M Kirka, and Peeyush Nandwana. Effects of heat treatments on microstructure and properties of ti-6al-4v eli alloy fabricated by electron beam melting (ebm). Materials Science and Engineering: A, 685:417–428, 2017. [34] AN Gent and IS Cho. Surface instabilities in compressed or bent rubber blocks. Rubber Chemistry and Technology, 72(2):253–262, 1999. [35] Jan Genzer and Jan Groenewold. Soft matter with hard skin: From skin wrin- kles to templating and material characterization. Soft Matter, 2(4):310–323, 2006. 113 [36] J Giannatsis and V Dedoussis. Additive fabrication technologies applied to medicine and health care: a review. The International Journal of Advanced Manufacturing Technology, 40(1-2):116–127, 2009. [37] John Goldak, Aditya Chakravarti, and Malcolm Bibby. A new finite element model for welding heat sources. Metallurgical transactions B, 15(2):299–305, 1984. [38] Haijun Gong, Khalid Rafi, Hengfeng Gu, Thomas Starr, and Brent Stucker. Analysis of defect generation in ti–6al–4v parts made using powder bed fusion additive manufacturing processes. Additive Manufacturing, 1:87–98, 2014. [39] Nannan Guo and Ming C Leu. Additive manufacturing: technology, applica- tions and research needs. Frontiers of Mechanical Engineering, 8(3):215–243, 2013. [40] Brian J Hayes, Brian W Martin, Brian Welk, Samuel J Kuhr, Thomas K Ales, David A Brice, Iman Ghamarian, Andrew H Baker, Christina V Haden, D Gary Harlow, et al. Predicting tensile properties of ti-6al-4v produced via directed energy deposition. Acta Materialia, 133:120–133, 2017. [41] Hibbett, Karlsson, and Sorensen. ABAQUS/standard: User’s Manual, vol- ume 1. Hibbitt, Karlsson & Sorensen, 1998. [42] Douglas P Holmes and Alfred J Crosby. Draping films: A wrinkle to fold transition. Physical review letters, 105(3):038303, 2010. [43] Wei Hong, Xuanhe Zhao, and Zhigang Suo. Formation of creases on the surfaces of elastomers and gels. Applied Physics Letters, 95(11):111901, 2009. [44] CO Horgan. Remarks on ellipticity for the generalized blatz-ko constitu- tive model for a compressible nonlinearly elastic solid. Journal of elasticity, 42(2):165–176, 1996. [45] R Huang and Z Suo. Instability of a compressed elastic film on a viscous layer. International Journal of Solids and Structures, 39(7):1791–1802, 2002. [46] Samuel H Huang, Peng Liu, Abhiram Mokasdar, and Liang Hou. Additive manufacturing and its societal impact: a literature review. The International Journal of Advanced Manufacturing Technology, 67(5-8):1191–1203, 2013. [47] Wilhelm TS Huck. Artificial skins: hierarchical wrinkling. Nature materials, 4(4):271–272, 2005. [48] Lihua Jin, Anesia Auguste, Ryan C Hayward, and Zhigang Suo. Bifurcation diagrams for the formation of wrinkles or creases in soft bilayers. Journal of Applied Mechanics, 82(6):061008, 2015. 114 [49] Sujoy Kumar Kar. Modeling of mechanical properties in alpha/beta-titanium alloys. PhD thesis, The Ohio State University, 2005. [50] Shawn Michael Kelly. Thermal and microstructure modeling of metal deposition processes with application to Ti-6Al-4V. PhD thesis, Virginia Tech, 2004. [51] SM Kelly, SS Babu, SA David, T Zacharia, and SL Kampe. A microstruc- ture model for laser processing of ti-6al-4v. In Trends In Welding Research: Proceedings of the 7Th International Conference, page 65. ASM International, 2006. [52] Dahl-Young Khang, Hanqing Jiang, Young Huang, and John A Rogers. A stretchable form of single-crystal silicon for high-performance electronics on rubber substrates. Science, 311(5758):208–212, 2006. [53] Jungwook Kim, Jinhwan Yoon, and Ryan C Hayward. Dynamic display of biomolecular patterns through an elastic creasing instability of stimuli- responsive hydrogels. Nature materials, 9(2):159, 2010. [54] Pilnam Kim, Manouk Abkarian, and Howard A Stone. Hierarchical fold- ing of elastic membranes under biaxial compressive stress. Nature materials, 10(12):952, 2011. [55] W King, AT Anderson, RM Ferencz, NE Hodge, C Kamath, and SA Khairal- lah. Overview of modelling and simulation of metal powder bed fusion process at lawrence livermore national laboratory. Materials Science and Technology, 31(8):957–968, 2015. [56] YK Ko, V Sahai, RA Overfelt, and JT Berry. Prediction of porosity in cast equiaxed alloy 718. Technical report, Minerals, Metals and Materials Society, Warrendale, PA (United States), 1995. [57] Pamela A Kobryn and SL Semiatin. Microstructure and texture evolution dur- ing solidification processing of ti–6al–4v. Journal of Materials Processing Tech- nology, 135(2-3):330–339, 2003. [58] DP Koistinen. A general equation prescribing the extent of the austenite- martensite transformation in pure iron-carbon alloys and plain carbon steels. acta metallurgica, 7:59–60, 1959. [59] Warner Tjardus Koiter. The stability of elastic equilibrium. Technical report, STANFORD UNIV CA DEPT OF AERONAUTICS AND ASTRONAUTICS, 1970. [60] Andrei Nikolaevich Kolmogorov. On the statistical theory of the crystallization of metals. Bull. Acad. Sci. USSR, Math. Ser, 1:355–359, 1937. [61] J-P Kruth, Ming-Chuan Leu, and Terunaga Nakagawa. Progress in additive manufacturing and rapid prototyping. Cirp Annals, 47(2):525–540, 1998. 115 [62] Jean-Pierre Kruth. Material incress manufacturing by rapid prototyping tech- niques. CIRP Annals-Manufacturing Technology, 40(2):603–614, 1991. [63] Jean-Pierre Kruth, Ludo Froyen, Jonas Van Vaerenbergh, Peter Mercelis, Mar- leen Rombouts, and Bert Lauwers. Selective laser melting of iron-based powder. Journal of materials processing technology, 149(1-3):616–622, 2004. [64] M K¨ucken and AC Newell. A model for fingerprint formation. EPL (Europhysics Letters), 68(1):141, 2004. [65] S Lampman. Wrought titanium and titanium alloys. Metal handbook, 2:592– 633, 1990. [66] Gideon N Levy, Ralf Schindel, and Jean-Pierre Kruth. Rapid manufacturing and rapid tooling with layer manufacturing (lm) technologies, state of the art and future perspectives. CIRP Annals-Manufacturing Technology, 52(2):589– 609, 2003. [67] Christoph Leyens and Manfred Peters. Titanium and titanium alloys: funda- mentals and applications. John Wiley & Sons, 2003. [68] Bo Li, Yan-Ping Cao, Xi-Qiao Feng, and Huajian Gao. Mechanics of mor- phological instabilities and surface wrinkling in soft materials: a review. Soft Matter, 8(21):5728–5745, 2012. [69] Haiyi Liang and L Mahadevan. Growth, geometry, and mechanics of a blooming lily. Proceedings of the National Academy of Sciences, 108(14):5516–5521, 2011. [70] Pei-Chun Lin, Shilpi Vajpayee, Anand Jagota, Chung-Yuen Hui, and Shu Yang. Mechanically tunable dry adhesive from wrinkled elastomers. Soft Matter, 4(9):1830–1835, 2008. [71] S Malinov, W Sha, and Z Guo. Application of artificial neural network for pre- diction of time–temperature–transformation diagrams in titanium alloys. Ma- terials Science and Engineering: A, 283(1-2):1–10, 2000. [72] Justin Mezzetta. Process-property relationships of Ti6Al4V fabricated through selective laser melting. PhD thesis, Ph. D. thesis, McGill University, Montreal, QC, Canada, 2016. [73] Scott Mitzner, Stephen Liu, Marcia S Domack, and Robert A Hafley. Grain refinement of freeform fabricated ti-6al-4v alloy using beam/arc modulation. 2012. [74] C Charles Murgau. Microstructure Model for Ti-6Al-4V used in Simulation of Additive Manufacturing. PhD thesis, Doctoral thesis, Lule˚ a University of Technology, 2016. 116 [75] C Charles Murgau, Robert Pederson, and Lars-Erik Lindgren. A model for ti– 6al–4v microstructure evolution for arbitrary temperature changes. Modelling and Simulation in Materials Science and Engineering, 20(5):055006, 2012. [76] Lawrence E Murr, Sara M Gaytan, Diana A Ramirez, Edwin Martinez, Jen- nifer Hernandez, Krista N Amato, Patrick W Shindo, Francisco R Medina, and Ryan B Wicker. Metal fabrication by additive manufacturing using laser and electron beam melting technologies. Journal of Materials Science & Technology, 28(1):1–14, 2012. [77] Pulin Nie, OA Ojo, and Zhuguo Li. Numerical modeling of microstructure evolution during laser additive manufacturing of a nickel-based superalloy. Acta Materialia, 77:85–95, 2014. [78] ASTM Committee F42 on Additive Manufacturing Technologies and ASTM Committee F42 on Additive Manufacturing Technologies. Subcommittee F42. 91 on Terminology. Standard terminology for additive manufacturing technolo- gies. ASTM International, 2012. [79] Deepankar Pal, Nachiket Patil, Kai Zeng, and Brent Stucker. An integrated approach to additive manufacturing simulations using physics based, coupled multiscale process modeling. Journal of Manufacturing Science and Engineer- ing, 136(6):061022, 2014. [80] W Piekarska, M Kubiak, and A Bokota. Numerical simulation of thermal phe- nomena and phase transformations in laser-arc hybrid welded joints. Archives of Metallurgy and Materials, 56(2):409–421, 2011. [81] Luka Pocivavsek, Robert Dellsy, Andrew Kern, Sebasti´an Johnson, Binhua Lin, Ka Yee C Lee, and Enrique Cerda. Stress and fold localization in thin elastic membranes. Science, 320(5878):912–916, 2008. [82] Remi Ponche, Olivier Kerbrat, Pascal Mognol, and Jean-Yves Hascoet. A novel methodology of design for additive manufacturing applied to additive laser manufacturing process. Robotics and Computer-Integrated Manufactur- ing, 30(4):389–398, 2014. [83] Patcharapit Promoppatum, Shi-Chune Yao, P Chris Pistorius, Anthony D Rol- lett, Peter J Coutts, Frederick Lia, and Richard Martukanitz. Numerical mod- eling and experimental validation of thermal history and microstructure for additive manufacturing of an inconel 718 product. Progress in Additive Manu- facturing, pages 1–18, 2018. [84] Narendran Raghavan, Ryan Dehoff, Sreekanth Pannala, Srdjan Simunovic, Michael Kirka, John Turner, Neil Carlson, and Sudarsanam S Babu. Numerical modeling of heat-transfer and the influence of process parameters on tailoring the grain morphology of in718 in electron beam additive manufacturing. Acta Materialia, 112:303–314, 2016. 117 [85] David P Richman, R Malcolm Stewart, John W Hutchinson, and Verne S Caviness. Mechanical model of brain convolutional development. Science, 189(4196):18–21, 1975. [86] J Romero, MM Attallah, M Preuss, M Karadge, and SE Bray. Effect of the forging pressure on the microstructure and residual stress development in ti– 6al–4v linear friction welds. Acta Materialia, 57(18):5582–5592, 2009. [87] E Sallica-Leva, R Caram, AL Jardini, and JB Fogagnolo. Ductility improve- ment due to martensite α decomposition in porous ti–6al–4v parts produced by selective laser melting for orthopedic implants. Journal of the mechanical behavior of biomedical materials, 54:149–158, 2016. [88] N Saunders, Z Guo, X Li, AP Miodownik, and J Ph Schill´e. The calculation of ttt and cct diagrams for general steels. JMatPro Software Literature, 2004. [89] TW Shield, KS Kim, and RT Shield. The buckling of an elastic layer bonded to an elastic substrate in plane strain. Journal of Applied Mechanics, 61(2):231– 235, 1994. [90] Winston O Soboyejo and TS Srivatsan. Advanced structural materials: proper- ties, design optimization, and applications. CRC press, 2006. [91] Yuhua Song, Yongnian Yan, Renji Zhang, Da Xu, and Feng Wang. Manufacture of the die of an automobile deck part based on rapid prototyping and rapid tooling technology. Journal of materials processing technology, 120(1-3):237– 242, 2002. [92] Jeong-Yun Sun, Shuman Xia, Myoung-Woon Moon, Kyu Hwan Oh, and Kyung- Suk Kim. Folding wrinkles of a thin stiff layer on a soft substrate. In Proceedings of the Royal Society of London A: Mathematical, Physical and Engineering Sciences, volume 468, pages 932–953. The Royal Society, 2012. [93] PK Sung, DR Poirier, SD Felicelli, EJ Poirier, and A Ahmed. Simulations of microporosity in in718 equiaxed investment castings. Journal of crystal growth, 226(2-3):363–377, 2001. [94] Tuomas Tallinen and John S Biggins. Mechanics of invagination and folding: Hybridized instabilities when one soft tissue grows on another. Physical Review E, 92(2):022720, 2015. [95] HP Tang, M Qian, N Liu, XZ Zhang, GY Yang, and J Wang. Effect of powder reuse times on additive manufacturing of ti-6al-4v by selective electron beam melting. Jom, 67(3):555–563, 2015. [96] Lore Thijs, Frederik Verhaeghe, Tom Craeghs, Jan Van Humbeeck, and Jean- Pierre Kruth. A study of the microstructural evolution during selective laser melting of ti–6al–4v. Acta materialia, 58(9):3303–3312, 2010. 118 [97] Charles L Thomas, Thomas M Gaffney, Srinivas Kaza, and Cheol H Lee. Rapid prototyping of large scale aerospace structures. In Aerospace Applications Con- ference, 1996. Proceedings., 1996 IEEE, volume 4, pages 219–230. IEEE, 1996. [98] J Tiley, T Searles, E Lee, S Kar, Russ Banerjee, JC Russ, and HL Fraser. Quantification of microstructural features in α/β titanium alloys. Materials Science and Engineering: A, 372(1-2):191–198, 2004. [99] Sakya Tripathy, Charlie Chin, Tyler London, Utkarsha Ankalkhope, and Vic- tor Oancea. Process modeling and validation of powder bed metal additive manufacturing. [100] Yuhki Tsukada, Yasuhiro Kojima, Toshiyuki Koyama, and Yoshinori Murata. Phase-field simulation of habit plane formation during martensitic transforma- tion in low-carbon steels. ISIJ International, 55(11):2455–2462, 2015. [101] SA Usmani and MF Beatty. On the surface instability of a highly elastic half- space. Journal of Elasticity, 4(4):249–263, 1974. [102] Mohammad Vaezi, Srisit Chianrabutra, Brian Mellor, and Shoufeng Yang. Mul- tiple material additive manufacturing–part 1: a review: this review paper covers a decade of research on multiple material additive manufacturing technologies which can produce complex geometry parts with different materials. Virtual and Physical Prototyping, 8(1):19–50, 2013. [103] Bey Vrancken, Lore Thijs, Jean-Pierre Kruth, and Jan Van Humbeeck. Heat treatment of ti6al4v produced by selective laser melting: Microstructure and mechanical properties. Journal of Alloys and Compounds, 541:177–185, 2012. [104] Gerhard Welsch, Rodney Boyer, and EW Collings. Materials properties hand- book: titanium alloys. ASM international, 1993. [105] Johnson William and Robert Mehl. Reaction kinetics in processes of nucleation and growth. Transactions of the Metallurgical Society of AIME, 135:416–442, 1939. [106] Jiawen Xie, V. Oancea, and Juan A Hurtado. Phase transformations in metals during additive manufacturing processes. NAFEMS World Congress 2017, 2017. [107] Jingjing Yang, Hanchen Yu, Jie Yin, Ming Gao, Zemin Wang, and Xiaoyan Zeng. Formation and control of martensite in ti-6al-4v alloy produced by selec- tive laser melting. Materials & Design, 108:308–318, 2016. [108] Jianfeng Zang, Xuanhe Zhao, Yanping Cao, and John W Hutchinson. Localized ridge wrinkling of stiff films on compliant substrates. Journal of the Mechanics and Physics of Solids, 60(7):1265–1279, 2012. 119 [109] Tom´as Zegard and Glaucio H Paulino. Bridging topology optimization and ad- ditive manufacturing. Structural and Multidisciplinary Optimization, 53(1):175– 192, 2016. [110] Qi Zhang and Janet A Blume. Surface wrinkling in generalized blatz–ko mate- rials. Extreme Mechanics Letters, 11:68–76, 2017. [111] Haiyan Zhao, Lugui He, Wen Chong Niu, and Bin Zhang. Investigation on porosity suppression in deep-penetration laser welding by using computational fluid dynamics. Journal of Laser Applications, 28(3):032011, 2016. [112] R Zhao, T Zhang, M Diab, H Gao, and K-S Kim. The primary bilayer ruga- phase diagram i: Localizations in ruga evolution. Extreme Mechanics Letters, 4:76–82, 2015. [113] Jianhua Zhou, Yuwen Zhang, and JK Chen. Numerical simulation of random packing of spherical particles for powder-based additive manufacturing. Journal of manufacturing science and engineering, 131(3):031004, 2009. 120