CHARACTERIZATION AND PREDICTION OF MATERIAL RESPONSE OF MICRO AND NANO-UNDERFILLS FOR FLIP CHIP DEVICES Except where reference is made to the work of others, the work described in this dissertation is my own or was done in collaboration with my advisory committee. This dissertation does not include proprietary or classified information. ________________________ Muhammad Saiful Islam Certificate of Approval: _______________________________ Pradeep Lall, Co-Chair Associate Professor Mechanical Engineering ______________________________ Jeffrey C. Suhling, Co-Chair Quina Distinguished Professor Mechanical Engineering _____________________________ Richard C. Jaeger Distinguished University Professor Electrical and Computer Engineering ____________________________ Roy W. Knight Assistant Professor Mechanical Engineering ___________________________ Stephen L. McFarland Dean Graduate School CHARACTERIZATION AND PREDICTION OF MATERIAL RESPONSE OF MICRO AND NANO-UNDERFILLS FOR FLIP CHIP DEVICES Muhammad Saiful Islam A Dissertation Submitted to the Graduate Faculty of Auburn University in Partial Fulfillment of the Requirements of the Degree of Doctor of Philosophy Auburn, Alabama May 11, 2006 iii CHARACTERIZATION AND PREDICTION OF MATERIAL RESPONSE OF MICRO AND NANO-UNDERFILLS FOR FLIP CHIP DEVICES Muhammad Saiful Islam Permission is granted to Auburn University to make copies of this dissertation at its discretion, upon request of individuals or institutions and at their expense. The author reserves all publication rights. ______________________________ Signature of Author ______________________________ Date of Graduation iv VITA Muhammad Saiful Islam, son of Muhammad Abdul Malek was born on July 01, 1973 in Pirojpur, Bangladesh. He graduated from Adamjee Cantonment College, in Dhaka, Bangladesh, in 1990. He was admitted to Bangladesh University of Engineering and Technology (BUET) in Dhaka, Bangladesh, and graduated with the degree of Bachelor of Science in Mechanical Engineering in July 1997. After graduation, he joined the Department of Mechanical Engineering at BUET in September 1997 as a lecturer. He served two years in the Department of Mechanical Engineering at BUET as a faculty member and then came to the US in 1999 for higher studies. He received his Master of Science in Mechanical Engineering degree in May 2005 from Auburn University, Auburn, Alabama. v DISSERTATION ABSTRACT CHARACTERIZATION AND PREDICTION OF MATERIAL RESPONSE OF MICRO AND NANO-UNDERFILLS FOR FLIP CHIP DEVICES Muhammad Saiful Islam Doctor of Philosophy, May 11, 2006 (M. S., Auburn University, 2005) (B. S., Bangladesh University of Engineering and Technology, 1997) 276 Typed Pages Directed by Pradeep Lall and Jeffrey C. Suhling Silica particles are used as a filler material in electronic underfills to reduce coefficient of thermal expansion of the underfill-epoxy matrix. In traditional underfills, the size of silica particles is in the micrometer range. Reduction in particle sizes into the nanometer range has the potential of attaining higher volume fraction particle loading in the underfills and greater control over underfill properties for higher reliability applications. Presently, no-flow underfills have very low or no filler content because micron-size filler particles hinder solder joint formation. Nano-silica underfills have the potential of attaining higher filler loading in no-flow underfills without hindering solder interconnect formation [1]. vi In this research, property prediction models based on representative volume element (RVE) and modified random spatial adsortion have been developed. A size distribution also applied to the spherical particles to generate random diameter of the filler particles. The models can be used for development of nano-silica underfills with desirable thermo-mechanical properties. Elastic Modulus, coefficient of thermal expansion, bulk modulus, Poisson?s ratio and viscoelastic properties have been predicted. Temperature dependent thermo-mechanical properties of nano-silica underfills and micro-silica underfill have been evaluated and correlated with the model prediction in a temperature range of ?175?C to +150 ?C. Properties investigated include temperature dependent stress-strain, creep and stress relaxation behavior. The trade-offs between using nano-fillers instead of micron-fillers on thermo-mechanical properties and reliability has been benchmarked. vii ACKNOWLEDGEMENT The author would like to thank his advisors Dr. Pradeep Lall, Dr. Jeffrey C. Suhling and the other internal and external committee members for their continuous guidance and encouragement. The author also thanks the NSF Center for Advanced Vehicle Electronics (CAVE); the NASA Center for Space Power, and Advanced Electronics (CSPAE); and Auburn University; for supporting and funding this project. Very special thanks go to his daughter Sharia Islam Maisha and his wife Mehroze Sultana for their support to this effort. Many thanks go to his parents Muhammad Abdul Malek and Monoara Begum, his father-in-law Md. Mahboob Ali and his mother-in-law Shahanara Begum for their inspirations. viii Style manual and journal used: Guide to Preparation and Submission of Theses and Dissertations Computer software used: Microsoft Word 2000, Microsoft Excel, Sigma Plot, ANSYS etc. ix TABLE OF CONTENTS LIST OF FIGURES xii LIST OF TABLES xx CHAPTERS 1 INTRODUCTION 1 1.1 Introduction 1 1.2 Research Objective 13 2 LITERATURE REVIEW 15 2.1 Literature Review on Underfill Materials 15 3 CHARACTERIZATION OF UNDERFILL MATERIAL PROPERTIES 19 3.1 Characterization of Nano Silica Underfill 19 3.1.1 Micro Tension/Torsion Testing Machine 19 3.1.2 Uniaxial Testing of Nano Underfill 1 (NUF1) 21 3.1.3 Uniaxial Testing of Nano Underfill 2 (NUF2) 35 3.1.4 CTE Measurement 38 3.2 Characterization of Macro Silica Underfill 47 3.2.1 Effect of Extreme Low Temperatures 47 3.2.2 Temperature Dependent Creep Curves 51 4 Finite Element Simulation of Uniaxial Test Specimen 53 4.1 Effect of Length to Width Aspect Ratio 53 4.2 Nonlinear Simulation of Temperature Dependent Stress Strain Curves 76 4.3 Error Estimation by Multi-mesh Extrapolation 82 5 PREDICTION OF MATERIAL PROPERTIES MACRO-SILICA UNDERFILLS 88 5.1 Underfill Constituent Properties of Macro-Silica Underfills 88 x 5.2 Prediction of Material Properties of Macro Silica Underfills 88 5.2.1 Unit Cell Generation 88 5.2.2 Prediction of Coefficient of Thermal Expansion by 1 Filler Model 90 5.2.3 Prediction of Coefficient of Thermal Expansion by 8 Filler Model 97 5.2.4 Prediction of Elastic Modulus by 1 Filler Model 102 5.2.5 Prediction of Elastic Modulus by 8 Filler Model 102 6 PREDICTION OF MATERIAL PROPERTIES OF NANO-SILICA UNDERFILLS BY CONSTANT SIZE FILLER MODELS 110 6.1 Underfill Constituent Properties of Nano-Silica Underfills 110 6.2 Prediction of Material Properties of Nano Silica Underfills 112 6.2.1 Unit Cell Generation 112 6.2.2 Isotropy of the Unit Cell 115 6.2.3 Randomness of Filler Distribution 116 6.2.4 Finite Element Model of Unit Cell 120 6.2.5 Prediction of Coefficient of Thermal Expansion (CTE) 120 6.2.6 Repeatability of Unit Cell 128 6.2.7 Prediction of Elastic Modulus 128 7 PREDICTION OF MATERIAL PROPERTIES OF NANO-SILICA UNDERFILLS BY RANDOM SIZE FILLER MODELS 132 7.1 Unit Cell Generation 132 7.2 Finite Element Model of Unit Cell 137 7.3 Prediction of Coefficient of Thermal Expansion (CTE) 138 7.4 Prediction of Elastic Modulus 138 7.5 Prediction of Bulk Modulus 145 7.6 Prediction of Poisson?s Ratio 147 7.7 Theoretical Background of Viscoelastic Model 147 7.8 ANSYS Input Constants for Viscoelastic Material 155 7.9 Prediction of Relaxation Behavior 171 8 CONCLUSIONS 177 REFERENCES 178 APPENDICES 183 A FILLER DISTRIBUTION IN DIFFERENT UNIT CELLS 184 B MATLAB PROGRAM FOR RANDOM COORDINATE GENERATION-I 217 xi C MATLAB PROGRAM FOR RANDOM COORDINATE GENERATION-II 224 D MATLAB PROGRAM FOR RANDOM COORDINATE GENERATION-III 236 E MATLAB PROGRAM FOR RANDOM COORDINATE GENERATION-IV 243 LIST OF FIGURES 1.1 Cross-Sectional Schematic of a Flip Chip 3 1.2 Shear Deformation in Solder Joint Due to Heating of the Flip Chip Assembly without Underfill 5 1.3 Shear Deformation in Solder Joint Due to Cooling of the Flip Chip Assembly without Underfill 6 1.4 Figure 1.4: Warping During Heating of the Flip Chip Assembly with Underfill 8 1.5 Warping During Cooling of the Flip Chip Assembly with Underfill 9 3.1 MT-200 Testing System with Environmental Chamber 20 3.2 Typical Stress-Strain Curves of Underfill NUF1. (T = 50 ?C, = 0.001 sec ?& -1 ) 23 3.3 Hyperbolic Tangent Model Fit to Typical Underfill Stress-Strain Data (T = 50 ?C, = 0.001 sec?& -1 ) 24 3.4 Temperature Dependent Average Stress-Strain Curves of Underfill NUF1 ( = 0.001 sec? & -1 ) 26 3.5 Temperature Dependent Elastic Modulus of NUF1 (?& = 0.001 sec -1 ) 27 3.6 Elastic Modulus of a Nano-Filler Underfill and a Micro Filler Underfill (?& = 0.001 sec -1 ) 28 3.7 Cryogenic Test Chamber for MT 200 Testing System 30 3.8 Temperature Dependent Stress Strain Curves at Extreme Low Temperatures (? = 0.001 sec & -1 ) 31 3.9 Figure 3.9: Temperature Dependent Elastic Modulus in Cryogenic Temperatures (? = 0.001 sec& -1 ) 32 xii 3.10 Temperature Dependent Creep Data of NUF1 33 3.11 Temperature Dependent Stress Relaxation Test Data of Underfill NUF1 34 3.12 Temperature Dependent Stress Strain Curves (?& = 0.001 sec -1 ) 36 3.13 Temperature Dependent Elastic Modulus (?& = 0.001 sec -1 ) 37 3.14 Elastic Modulus of Nano-Filler Underfills and a Micro Filler Underfill (?& = 0.001 sec -1 ) 39 3.15 Temperature Dependent Creep Data of NUF2 40 3.16 Underfill Specimen with Strain Gage 42 3.17 Strain Gage on Reference Material. (Titanium-Silicate Bar,TSB) 43 3.14 Strain Gage Reading from the Gage on the Reference Material 44 3.19 Strain Gage Reading from the Gage on the Test Specimen 45 3.20 Thermal Strain on Underfill Specimen at Different Temperature 46 3.21 Underfill Stress-Strain Curves for Material UF2 (T = -175 to 150 o C, = 0.001 sec ?& -1 ) 49 3.22 Elastic Modulus vs. Temperature for Material UF2 (T = -175 to 150 o C, ?& = 0.001 sec -1 ) 50 3.23 Temperature Dependent Underfill Creep Data (Stress Level 20 MPa) 52 4.1 FEA Model of Symmetric Half of the Tensile Specimen. Width to Length Aspect Ratio = 40 54 4.2 FEA Mesh of the Model 55 4.3 FEA Mesh of the Model. Specimen Length, L = 6 mm, Width, W = 3 mm and Length to Width Aspect Ratio = 2 57 4.4 Axial Stress Distribution (? x ). Length to Width Aspect Ratio = 2 58 4.5 Axial Stress (? x ) Value along the Right Most Vertical Path. Length to Width Aspect Ratio = 2 59 xiii xiv 4.6 Comparison of FEA Output with the Input Experimental Data, Length to Width Aspect Ratio = 2 60 4.7 FEA Mesh of the Model. Specimen Length, L = 18 mm, Width, W = 3 mm and Length to Width Aspect Ratio = 6 61 4.8 Axial Stress Distribution (? x ). Length to Width Aspect Ratio = 6 62 4.9 Axial Stress (? x ) Value along the Right Most Vertical Path. Length to Width Aspect Ratio = 6 63 4.10 Comparison of FEA Output with the Input Experimental Data. Length to Width Aspect Ratio = 6 64 4.11 FEA Mesh of the Model. Specimen Length, L = 60 mm, Width, W = 3 mm and Length to Width Aspect Ratio = 20 66 4.12 Axial Stress Distribution (? x ). Length to Width Aspect Ratio = 20 67 4.13 Axial Stress (? x ) Value Along the Right Most Vertical Path. Length to Width Aspect Ratio = 20 68 4.14 Comparison of FEA Output with the Input Experimental Data. Length to Width Aspect Ratio = 20 69 4.15 FEA Mesh of the Model. Specimen Length, L = 120 mm, Width, W = 3 mm and Length to Width Aspect Ratio = 40 70 4.16 Axial Stress Distribution (? x ). Length to Width Aspect Ratio = 40 71 4.17 Axial Stress (? x ) Value along the Right Most Vertical Path. Length to Width Aspect Ratio = 40 72 4.18 Comparison of FEA Output with the Input Experimental Data. Length to Width Aspect Ratio = 40 73 4.19 Effect of Length to Width Aspect Ratio on Elastic Modulus 74 4.20 Percent Difference in Elastic Modulus Values from the Elastic Modulus Obtained for Length to Width Aspect Ratio 40 75 4.21 Axial Stress along Top Horizontal Path. Length to Width Aspect Ratio = 40 77 4.22 Axial Stress along Mid Horizontal Path. Length to Width Aspect Ratio = 78 xv 4.23 Axial Stress along Bottom Horizontal Path. Length to Width Aspect Ratio=40 79 4.24 Comparison of Stress-Strain Curves Generated by FEA with Input Experimental Data. Length to Width Aspect Ratio = 40 80 4.25 FEA Prediction of Stress-Strain Curves at 110 ?C. Length to Width Aspect Ratio = 40 81 4.26 Original Mesh Before Refinement. Total Number of Element = 720. Length to Width Aspect Ratio = 40 83 4.27 Mesh Refinement. Total Number of Element = 2880. Length to Width Aspect Ratio = 40 84 4.28 More Mesh Refinement. Total Number of Element = 11520. Length to Width Aspect Ratio = 40 85 4.29 Axial Stress vs. h Plot. Calculations are done for the Mid Node of the Rightmost Vertical Edge 86 5.1 One Filler Unit Cell 91 5.2 Eight Filler Unit Cell 92 5.3 Finite Element Model of One Filler Unit Cell 93 5.4 Applied Boundary Condition in One Filler Unit Cell 94 5.5 Thermal Expansion in x Direction 95 5.6 U x Displacement of the Deformed Surface at X = 2.3 ?m from Center [Right Most Surface] 96 5.7 Finite Element Model of Eight Filler Unit Cell 98 5.8 Applied Boundary Condition in Eight Filler Unit Cell 99 5.9 Thermal Expansion in x Direction 100 5.10 U x Displacement of the Deformed Surface at X = 4.6 ?m from Center [Right Most Surface] 101 5.11 Applied Boundary Condition 103 xvi 5.12 Displacement in y Direction 104 5.13 Stress Strain Curve at the Top Surface of the Unit Cell 105 5.14 Applied Boundary Condition in Eight Filler Model 106 5.15 Displacement in y Direction 108 5.16 Stress Strain Curve Calculated at the top Surface of the Unit Cell 109 6.1 Centroids for Acceptable Nano-Particle Distribution 117 6.2 Plot of Radial Distribution Function 119 6.3 Isotropic View of Filler Distribution, ? = 0.20 121 6.4 Front View of Filler Distribution, ? = 0.20 122 6.5 RHS View of Filler Distribution, ? = 0.20 123 6.6 FEA Mesh, ? = 0.20 124 6.7 Isotropic View of Filler Distribution, ? = 0.38 125 6.8 Prediction of x-CTE by finite-element analysis of Unit Cell 126 6.9 Prediction of x, y, z-CTE by finite-element analysis of Unit Cell 127 6.10 Prediction of Elastic Modulus by FEA Analysis of Unit Cell 130 7.1 Isometric View of the Size and Position Distribution of the Filler Particles, ? = 0.20 139 7.2 Front View of the Size and Position Distribution of the Filler Particles, ? = 0.20 140 7.3 Right Side View of the Size and Position Distribution of the Filler Particles, ? = 0.20 141 7.4 Right Side View of the Size and Position Distribution of the Filler Particles, ? = 0.39 142 7.5 FEA Mesh of Unit Cell (Isometric View) 143 xvii 7.6 Prediction of CTE by finite-element analysis of Unit Cell 144 7.7 Prediction of Elastic Modulus by FEA Analysis of Unit Cell 146 7.8 Prediction of Bulk Modulus by FEA Analysis of Unit Cell 148 7.9 Prediction of Poisson?s Ratio by FEA Analysis of Unit Cell 149 7.10 Maxwell Model 151 7.11 Stress Relaxation of Maxwell Element under Constant Strain 152 7.12 Generalized Maxwell Model by Connecting Elements in Series 153 7.13 Generalized Maxwell Model by Connecting Elements in Parallel 154 7.14 Kelvin Model 156 7.15 Generalized Kelvin Model by Connecting the Elements in Parallel 157 7.16 Generalized Kelvin Model by Connecting the Elements in Series 158 7.17 Stress Relaxation Data Used for Calculating ANSYS Viscoelastic Constant 160 7.18 Log-Log Plot of the Relaxation Modulus versus Time Data 161 7.19 Master Relaxation Modulus at 25 ?C 162 7.20 Temperature Dependent Shift Factor 163 7.21 Generalized Parallel Maxwell Model with a Free Spring in Parallel 167 7.22 Prediction of Elastic Modulus Relaxation 173 7.23 Prediction of Shear Modulus Relaxation 175 7.24 Prediction of Bulk Modulus Relaxation 176 A.1 Random Distribution of Constant Size Fillers, ? = 0.10 [Isotropic View] 185 A2 Random Distribution of Constant Size Fillers, ? = 0.10 [Front View] 186 A3 Random Distribution of Constant Size Fillers, ? = 0.10 [RHS View] 187 xviii A.4 Random Distribution of Constant Size Fillers, ? = 0.10 [Top View] 188 A.5 Random Distribution of Constant Size Fillers, ? = 0.20 [Isotropic View] 189 A.6 Random Distribution of Constant Size Fillers, ? = 0.20 [Front View] 190 A.7 Random Distribution of Constant Size Fillers, ? = 0.20 [RHS View] 191 A.8 Random Distribution of Constant Size Fillers, ? = 0.20 [Top View] 192 A.9 Random Distribution of Constant Size Fillers, ? = 0.25 [Isometric View] 193 A.10 Random Distribution of Constant Size Fillers, ? = 0.25 [Front View] 194 A.11 Random Distribution of Constant Size Fillers, ? = 0.25 [RHS View] 195 A.12 Random Distribution of Constant Size Fillers, ? = 0.25 [Top View] 196 A.13 Random Distribution of Constant Size Fillers, ? = 0.38 [Isometric View] 197 A.14 Random Distribution of Constant Size Fillers, ? = 0.38 [Front View] 198 A.15 Random Distribution of Constant Size Fillers, ? = 0.38 [RHS View] 199 A.16 Random Distribution of Constant Size Fillers, ? = 0.38 [Top View] 200 A.17 Random Distribution of Random Size Fillers, ? = 0.10 [Isometric View] 201 A.18 Random Distribution of Random Size Fillers, ? = 0.10 [Front View] 202 A.19 Random Distribution of Random Size Fillers, ? = 0.10 [RHS View] 203 A.20 Random Distribution of Random Size Fillers, ? = 0.10 [Top View] 204 A.21 Random Distribution of Random Size Fillers, ? = 0.20 [Isometric View] 205 A.22 Random Distribution of Random Size Fillers, ? = 0.20 [Front View] 206 A.23 Random Distribution of Random Size Fillers, ? = 0.20 [RHS View] 207 A.24 Random Distribution of Random Size Fillers, ? = 0.20 [Top View] 208 A.25 Random Distribution of Random Size Fillers, ? = 0.25 [Isometric View] 209 xix A.26 Random Distribution of Random Size Fillers, ? = 0.25 [Front View] 210 A.27 Random Distribution of Random Size Fillers, ? = 0.25 [RHS View] 211 A.28 Random Distribution of Random Size Fillers, ? = 0.25 [Top View] 212 A.29 Random Distribution of Random Size Fillers, ? = 0.39 [Isometric View] 213 A.30 Random Distribution of Random Size Fillers, ? = 0.39 [Front View] 214 A.31 Random Distribution of Random Size Fillers, ? = 0.39 [RHS View] 215 A.32 Random Distribution of Random Size Fillers, ? = 0.39 [Top View] 216 xx LIST OF TABLES 1.1 CTE Values of Flip Chip Packaging Materials 4 3.1 Basic Underfill Properties as Provided by Vendor 48 4.1 Value of h and Axial Stress for Different Number of Element 87 5.1 Material Properties for Epoxy Matrix and Micro-Silica Particle 89 6.1 Material Properties for Epoxy Matrix and Nano-Silica Particle 111 6.2 Moment of Inertia of Acceptable Nano-Particle Distributions 118 6.3 Calculation of CTE Value to Check Repeatability of the Unit Cell 129 7.1 Prony Viscoelastic Shear Response 169 7.2 Prony Viscoelastic Volumetric Response 170 1 CHAPTER 1 INTRODUCTION 1.1 Introduction Packaging of electronic chips are performed to provide input/output signal connections, power connections, means of heat dissipation, structure to support and protect the chip etc. After manufacturing, the wafer is sliced into small chips. These chips are then packaged into chip carriers, known as first level packaging. In the second level packaging chip carriers are attached to a printed circuit board (PCB) or to a card. At last third level packaging is performed by placing the printed circuit board or the card on to a motherboard. Recent trends in electronic industries are to make smaller, cost effective, reliable and better performing devices. Flip chip, ball grid array (BGA) and chip scale package (CSP) are some examples of recent trend in IC packaging. Flip chip is one of the most attracted first level interconnection technology. In flip chip assembly the semiconductor chip is interconnected to a substrate by inverting the die (i.e. flipping the die so that the active side or the I/O side is in face down condition) and making electrical contacts and mechanical connections between the pads on the die and the pads on the substrate, circuit boards, or carriers, by means of conductive bumps. IBM first introduced flip chip in early 1960?s in the form of Controlled Collapse Chip Connection (C4). A cross-sectional 2 schematic of a flip chip is shown in Figure 1.1. The driving forces towards using flip chip are shortest electrical connection, area array connections, smallest size and weight, high speed production, lowest cost, etc. In flip chip the die is inverted and the I/O connection pads are directly connected with the conductive bump. This reduces the electrical path and hence increases the electrical performance. Since flip chip allows area array connections instead of only perimeter connections, it provides the possibility for more I/O connections per chip. Flip chip is the smallest and lightest packaging possible. This makes flip chip more suitable for portable electronics. In wire bonding, only one connection can be done at a time. In flip chip, all of the bump connections can be done simultaneously by a reflow process. Thus, all of the chips can be connected simultaneously on a printed circuit board. Flip chips are typically used in most cost-dominated applications. A major limitation of flip chip on laminate is it requires an underfill encapsulant to achieve acceptable thermal cycling reliability and impact resistance. In a flip chip without underfill, the mismatch in coefficient of thermal expansion (CTE) between the silicon chip and the substrate affects the reliability greatly. The CTE values of the materials used in a flip chip on laminate are shown in Table 1.1. The printed circuit board (PCB) has much higher CTE than the silicon chip. When a flip chip without underfill is heated during thermal cycling or thermal shock, the PCB expands more than the silicon chip. Thus the solder joints between them suffer sever shear deformation (Figure 1.2), which leads to early failure of the connections. When the flip chip is cooled down, the PCB contracts (Figure 1.3) more than silicon chip Silicon Chip Printed Circuit Board Copper Pad Passivation Solder mask Solder Joint Underfill Figure 1.1: Cross-Sectional Schematic of a Flip Chip 3 Table 1.1: CTE Values of Flip Chip Packaging Materials Material CTE (10 -6 1/?C) Silicon 2.3 Ceramic Substrate 4-6 FR4/Glass PCB 15-20 Underfill 25-30 63/37 Solder 21 Copper 17 4 At T = T 1 , Flip Chip Assembly At T= T 2 > T 1 , Shear Deformation in Solder Joint Due to Higher Expansion of PCB Figure 1.2: Shear Deformation in Solder Joint Due to Heating of the Flip Chip Assembly without Underfill 5 At T = T 1 , Flip Chip Assembly At T= T 2 < T 1 , Shear Deformation in Solder Joint Due to Higher Contraction of PCB Figure 1.3: Shear Deformation in Solder Joint Due to Cooling of the Flip Chip Assembly without Underfill 6 7 and causes the opposite sign of shear deformation on the solder joints. Underfills have been used as a supplemental restraint mechanism to enhance the reliability for flip-chip devices and chip-scale packages [1] in a wide variety of applications including portable consumer electronics like cellular phones [2], laptops [3], under-the-hood electronics [4], microwave applications [5], system in package (SIP) [6], high-end workstations [7], and several other high performance applications. Underfills compensate for the mismatch in coefficient of thermal expansion between silicon and the printed circuit board and more evenly distribute and minimize the solder joint strains, thus improving thermal cycling fatigue life. Underfill encapsulants are mainly non- conductive epoxies with silica or alumina particles as fillers. In the final stage of the flip chip assembly, liquid underfills are dispensed along the side of the die. After dispensing, it flows in the gap between the die and the substrate around the solder joints by capillary action. The underfill is then cured according to the recommendations of the manufacturer. After curing, the underfill creates a mechanical bond between the silicon die and PCB. When this flip chip assembly with underfill is heated or cooled, this bond reduces the differential deformation between the die and the substrate, and thus reduces the shear strain imposed on the solder joints. In addition, the whole assembly warps during temperature change (Figure 1.4 and 1.5). Epoxy is a common ingredient in the underfill material used because of its desirable characteristics like corrosion resistance, good adhesion, in addition to physical and electrical properties. However, epoxies by themselves possess a high coefficient of thermal expansion (above 80 ppm/?C), making them unable to meet the very first At T = T 1 , Flip Chip Assembly with Underfill At T= T 2 > T 1 , Warping Due to Higher Expansion of PCB Figure 1.4: Warping During Heating of the Flip Chip Assembly with Underfill 8 At T = T 1 , Flip Chip Assembly with Underfill At T = T 2 < T 1 , Warping Due to Higher Contraction of PCB Figure 1.5: Warping During Cooling of the Flip Chip Assembly with Underfill 9 10 requirement of a good underfill material. For this reason, the epoxies are filled with filler particles that decrease the CTE of the adhesive. This ensures lower thermal mismatch between the materials in the flip chip assembly. The primary reason of adding fillers are to control the thermomechanical, elastic and time dependent properties of underfill. Common filler particles include, silica, and ceramics. There are several other key requirements for underfill materials [8]. Relatively high elastic modulus is also desired for underfill materials. If the modulus is high, the tendency for flip chip assembly deformation is less. The elastic modulus of most underfill materials depends on temperature and rate of deformation. Good underfills should have relatively low cure temperatures. For example, the cure temperature of an underfill should not be more than the solder reflow temperature. Again a very high cure temperature may cause damage to the other components. Fast curing of underfills is also desired to minimize assembly times, while lower viscosity materials will ensure faster flow during dispense. Resistance to moisture absorption is another important characteristic of underfill materials. In addition, a relatively high glass transition temperature (T g ) is generally desirable. At temperatures near or above its glass transition temperature, a material looses its mechanical stiffness and strength, and becomes very soft. In addition, the CTE of a material greatly increases. In order to ensure mechanical reliability of flip chip assemblies at elevated temperatures, use of an underfill with a high T g is typically required. In traditional underfills, the size of silica particles is in the micrometer range. Reduction in particle sizes into the nanometer range has the potential of attaining higher volume fraction particle loading in the underfills and greater control over underfill 11 properties for higher reliability applications. Presently, no-flow underfills have very low or no filler content because micron-size filler particles hinder solder joint formation. Nano-silica underfills have the potential of attaining higher filler loading in no-flow underfills without hindering solder interconnect formation [9, 10] Nano-silica imparts the same modulus enhancement and CTE reduction to adhesives as micro-sized silica particles. However, nano-silica does not settle in an underfill formulation and does not interfere with the solder interconnection process, unlike micron-sized particles. Low filler-loading in no-flow underfills causes the coefficient of thermal expansion to be higher than capillary flow underfills. The no-flow underfills are largely unfilled or have very low filler loading because micron-sized filler particles interfere with the solder interconnection process [9, 10]. No-flow underfills are an attractive alternative to the use of capillary flow underfills, which typically require a post-reflow batch cure operation adding to production cycle time. In addition, nano-silica particles can achieve much higher volume fraction loading than micron-sized particles, thus providing greater control over underfill properties. Underfill formulation including volume fraction, size, distribution and material properties of the filler particles, influence the elastic modulus, coefficient of thermal expansion, and mechanical deformation behavior and determine the thermo-mechanical reliability of flip-chip devices. There is need for property-prediction techniques for formulation of the underfills with desirable thermo-mechanical characteristics. In this study methodologies for property prediction have been investigated. Models have been developed for prediction of underfill properties. The models are based on constituent component properties and enable the prediction of effective equivalent properties of statistically isotropic 12 composites formed by random distribution of spherical filler particles. Drugan [11, 12] showed that the representative volume element (RVE) is an effective technique for prediction of elastic composites. Segurado [13] also demonstrated the RVE with modified random sequential adsorption algorithm (RSA) as a reliable approach to estimate the equivalent properties. In the present study an algorithm similar to random sequential adsortion algorithm was developed to generate statistically isotropic cubic unit cells of underfill containing up to 38% nano fillers. Developed unit cell has been analyzed and the elastic modulus and coefficient of thermal expansion of the underfill were computed by using RVE implementation in implicit-finite elements. Reliable, consistent, and comprehensive material property data are needed for underfills for the purpose of mechanical design, reliability assessment, and process optimization of electronic packages. For example, it is common to use finite element models to simulate the reliability of flip chip on laminate assemblies subjected to thermal cycling (e.g. -40 to 125 ?C). Underfill mechanical property data are scarce on vendor datasheets, and typically only the room temperature elastic modulus is available. Given the known temperature dependent viscoplastic nature of filled epoxy encapsulants such as flip chip underfills, one must naturally doubt the accuracy of a finite element simulation based on only room temperature elastic properties. In this present work, temperature dependent thermo-mechanical properties of nano-underfills were evaluated in a temperature range of ? 175?C to +150 ?C. Underfills in electronic application are used at small scales typical of the gaps between the flip-chip device and the printed circuit assembly. Material properties in the present effort have been validated versus properties measured on micro-mechanical test structures. The test 13 structures are uniaxial-tension specimen with 75-125 ?m (3-5 mil) thickness, typical of the gaps between the silicon chip and the printed circuit board. The details of the specimen preparation procedure were presented in prior reference by the author [19]. A three parameter hyperbolic tangent empirical relation was used accurately to model the temperature dependent nonlinear stress-strain data. Strain gages were used to measure the CTE of cured underfill. Temperature dependent material behaviors of micro-silica underfills were also evaluated. 1.2 Research Objective ? Use micro tension/torsion testing machine to evaluate temperature dependent stress-strain, creep and relaxation behavior of nano-silica and micro-silica underfills ? Provide the basic underfill mechanical property data for future finite element analysis of flip chip ? Apply strain gage method to measure the coefficient of thermal expansion of nano-silica underfills ? Develop testing chamber for cryogenic temperature (-180 to 0 ?C) thermomechanical testing ? Develop finite element models to evaluate effect of aspect ratio on the uniaxial tensile specimen ? Develop finite element models to predict temperature dependent nonlinear stress- strain curves of underfills 14 ? Develop algorithm to generate random coordinates of nano filler particles to form statistically isotropic unit cell of nano-underfill ? Perform finite element analysis to predict effective elastic modulus, poisson?s ratio, shear modulus, bulk modulus, etc. of nano underfill ? Perform finite element analysis to predict effective creep and relaxation behavior of nano underfill 15 CHAPTER 2 LITERATURE REVIEW 2.1 Literature Review on Underfill Materials Recent trends of electronic industries towards making smaller, cost effective, reliable and better performing devices increases the necessity of using flip chip, ball grid array, chip scale packages etc. in various consumer electronics. Pascariu [3] demonstrated that in recent years the use of flip chip technology is increasing very rapidly. The author explained that the main reason behind the increased use of flip chip technology is the noticeable expansion of the demand for the internet, cellular phones, PDA?s, desktop computers, laptops, digital camcorder, digital cameras etc. A range of products including thousands ball grid array products uses flip chip technology. Use of wide range of organic substrates enabled development of various different applications of flip chip technology. Bedinger [5] described that the use of flip chip and ball grid array technology are very important in microwave power amplifier. The author also demonstrated that due to small size and mass for a flip chip thermal cycling is the primary mood of failure and vibration, mechanical shock etc. are not typically the mood of failure. Mismatch in the coefficient of thermal expansion is the main reason behind the thermal cycle failure. In the present work it was also explained that underfill can improve the solder joint fatigue 16 life significantly. This increase in life some can be practically achievable in the order of 8 to 9. Also adding of underfill does not change or degrade the performance of the device considerably for several flip chip technology. Liji [1] demonstrated that underfill provides a supplemental restraint mechanism to enhance the reliability of flip chip, CSP etc. The author found that use of underfill increases the life of solder joint interconnects in some PBGA samples in the order of 5 to 6. Author also discussed that although use of underfill increases the reliability, new failure moods are also introduced. For example, delamination at the underfill to die interface, PCB cracking etc. The author performed finite element study to see the stress situation in the underfilled PBGA packages. Obtained result showed that use of underfill reduces the maximum value of the solder joint stresses. Sillanpaa [2] analyzed the reliability of the flip chip assemblies. The author tested the assemblies for potential damage to the flip chips and solder interconnects due to the application of thermal cycling and mechanical shock by dropping. The failure analysis was also performed. The author performed finite element analysis to simulate the drop test. During the drop tests it was also found that the system is sensitive to the quality of underfill in drop tests. Ray [7] and coauthors developed a multilayer ceramic substrate to package a 335 mm 2 silicon chip which has about 15 millions transistors. The authors discussed the chip and package attributes in their paper. In their study a highly conductive thermal paste was developed by the author to provide efficient thermal path from the flip chip to the aluminum cap. The author achieved excellent reliability by column grid array connection to the second level organic card. 17 Liu and coauthors [9] introduced a series of single pass reflow encapsulants that requires no post curing after solder bump reflow. These materials first serve as a flux to activate soldering and then being cured during the reflow process and establish a mechanical connection between the die and the substrate. Their study employed DSC, DMA, TMA and shear test technique to investigate the properties of the material. The author concluded that the developed materials can reduce the assembly line throughput and cost as it eliminates the underfilling and post curing process. Shi and coworkers [10] investigated the filler interference mechanism of no flow underfill. The author demonstrated that the no-flow underfills are largely unfilled or have very low filler loading because micron-sized filler particles interfere with the solder interconnection process. Nano-silica underfills have the potential of attaining higher filler loading in no-flow underfills without hindering solder interconnect formation. Nano-silica underfill is a very new concept in underfill family. In these underfills nano-silica particles are dispersed in an epoxy matrix. Adding nano-silica particle instead of micro-silica particle will help to achieve better control of the underfill material property. These underfills have the potential of achieving higher volume fraction also. Very few studies have been performed recently on nano-silica underfills. In recent days development and application of nano-silica underfills are very important to achieve the success in electronic packaging and manufacturing industries. Knowledge of constitutive behavior of nano-silica underfill is very critical too. No study has been performed yet to characterize the thermomechanical properties of nano-silica underfill. The tradeoffs for moving towards nano-silica underfill instead of using the micro-silica underfill is not well understood too. 18 There is need for property-prediction techniques for formulation of the underfills with desirable thermo-mechanical characteristics. Finite element and analytical techniques need to be developed to predict properties of nano-silica underfills. Thermomechanical, elastic, viscoelastic, viscoplastic properties of nano-silica underfills need to be studied. Property prediction technique will provide the vendors an early idea about the thermomechanical behavior of new underfills before actually making it and testing it in the laboratory. This will also provide vendors an early idea about the reliability status of the electronic packages with any new underfill before actually making the underfill. This technique will also suggest necessary changes in the constituent component to obtain the desired underfill to optimize the reliability of electronic packages. Drugan [11, 12] showed that the representative volume element (RVE) is an effective technique for prediction of elastic composites. Segurado [13] also demonstrated the RVE with modified random sequential adsorption algorithm (RSA) as a reliable approach to estimate the equivalent properties. Segurado [13] also demonstrated that using a RVE with few dozen of fillers is a very good approach to estimate the equivalent properties of a particle reinforced composite. In the above study the author used modified random sequential adsorption algorithm to calculate the random position of spheres. CHAPTER 3 CHARACTERIZATION OF UNDERFILL MATERIAL PROPERTIES 3.1 Characterization of Nano Silica Underfill 3.1.1 Micro Tension/Torsion Testing Machine A MT-200 tension/torsion thermo-mechanical test system from Wisdom Technology, Inc., as shown in Figure 3.1, has been used to test the samples in this study. The system provides an axial displacement resolution of 0.1 micron and a rotation resolution of 0.001 o . Testing can be performed in tension, shear, torsion, bending, and in combinations of these loadings, on small specimens such as thin films, solder joints, gold wire, fibers, etc. Cyclic (fatigue) testing can also be performed at frequencies up to 5 Hz. In addition, a universal 6-axis load cell is utilized to simultaneously monitor three forces and three moments/torques during sample mounting and testing. An environmental chamber provides a temperature range capability of approximately ?180 to +300 o C. For uniaxial testing with the MT-200, forces and displacements are measured. The axial stress and axial strain are calculated from the applied force and measured cross-head displacement using LL L A F ? = ? =?=? (3.1) 19 Figure 3.1: MT-200 Testing System with Environmental Chamber 20 where is the uniaxial stress, ? ? is the uniaxial strain, F is the measured uniaxial force, A is the original cross-sectional area, ? is the measured cross-head displacement, and L is the specimen gage length (initial length between the grips). 3.1.2: Uniaxial Testing of Nano Underfill 1 (NUF1) Uniaxial testing specimens of the nano-underfill were developed. The specimen preparation procedure is presented elsewhere [19]. With the developed test specimen; tensile, creep and relaxation tests have been performed in a wide temperature range. The nano-underfill material for detailed material testing is NUF1, and has a cure time of 30 min at 150?C. The underfill has a volume fraction in the neighborhood of ? = 0.22. The glass transition temperature of this material is 155.86?C (DMA Method). The thickness of the cured uniaxial specimen was 75-125 ?m (3-5 mil). The typical length and width of a specimen are 90 mm and 3 mm respectively. In all uniaxial tests, the effective test length of the uniaxial specimen was 60 mm. Stress-Strain Data Figure 3.2 shows typical stress-strains curves for underfill NUF1 at 50 ?C, and a strain rate of = 0.001 sec?& -1 . The observed variation in the data between different tests is typical for cured polymeric materials. The elastic modulus E is the slope of the initial linear portion of the stress-strain curves. At 50?C and ?& = 0.001 sec -1 , the value for this underfill was measured to be E = 3.74 GPa. This value was found by averaging the results from 5 tests. An empirical three-parameter hyperbolic tangent model has been used to model the observed nonlinear underfill stress-strain data. Such a model has been 21 used historically to model the stress-strain curves of Cellulosic materials [17,33]. The general representation of the hyperbolic tangent empirical relation is ?+?=?? 321 C)Ctanh(C)( (3.2) where C 1 , C 2 , and C 3 are material constants. Differentiation of equation (3.2) gives an expression for the initial (zero strain) elastic modulus, 321 CCCE += (3.3) Likewise, constant C 3 represents the limiting slope of the stress-strain curve at high strains. For a given set of experimental data, constants C 1 , C 2 and C 3 are determined by performing a nonlinear regression analysis of equation (3.2) through experimental data points. Based on the results from reference [33], the data from all of the stress-strain curves in a set should be fit simultaneously in order to obtain the best set of hyperbolic tangent model material constants. The stress-strain data shown in Figure 3.2 were fit with the hyperbolic tangent model using a nonlinear regression analysis. Results from this calculation were C 1 = 53.60 MPa, C 2 = 68.67, and C 3 = 10.68 MPa. The data from Figure 3.2 and the best fit of the empirical model are shown together in and Figure 3.3. Excellent correlation is observed. Such results were typical for all of temperatures at which testing was performed. 22 Strain 0.000 0.005 0.010 0.015 0.020 0.025 S t re ss (M P a ) 0 10 20 30 40 50 60 Figure 3.2: Typical Stress-Strain Curves of Underfill NUF1. (T = 50 ?C, = 0.001 sec? & -1 ) 23 St rain 0.000 0.005 0.010 0.015 0.020 0.025 0 10 20 30 40 50 60 Model Stress Figure 3.3: Hyperbolic Tangent Model Fit to Typical Underfill Stress-Strain Data (T = 50 ?C, ?& = 0.001 sec -1 ) 24 Typical variation of the stress-strain curves of the tested underfill with temperature is shown in Figure 3.4. The strain rate for these tests was = 0.001 sec?& -1 . The total test time (to failure) of a typical tensile test is less than 5 seconds for room temperature experiments. Tests were performed at T = 25, 50, 75, 100, 125 and 150?C. The curves shown in Figure 3.4 are the hyperbolic tangent empirical fits to the multiple curves measured at each temperature, and therefore represent an average measured stress- strain curve at each temperature. The stress-strain curves in Figure 3.4 illustrate considerable softening and viscoplastic type behavior as the temperature is increased. Elastic modulus decreases approximately linearly with the increase of temperature from 25 to 125 ?C. Figure 3.5 shows temperature dependent elastic modulus of underfill NUF1. The value of elastic modulus at 25?C is in the neighborhood of 4.65 GPa. This correlates well with the predicted value of 4.4 GPa from the unit cell model in chapter 6. For temperature, T = 125 ?C elastic modulus decreases dramatically and at T = 150 ?C elastic modulus is almost near to zero. This is typical as underfill approaches near to its glass transition temperature. Figure 3.6 shows a comparison of elastic modulus of the nano-underfill NUF1 versus micron-filler underfill UF1. UF1 is has micro size particles as fillers. The volume fraction of filler particle in UF1 is almost double than NUF1. Figure 3.6 shows that due to higher filler concentration UF1 has higher elastic modulus than NUF1. Glass transition temperature of UF1 is 150 ?C and its modulus drops greatly after 100 ?C but NUF1 has glass transition temperature as 155.86 ?C and its modulus drops greatly after 125 ?C. Measurements of accurate mechanical properties at extreme low temperatures are very important for various applications including the solar system exploration missions 25 Strain 0.000 0.005 0.010 0.015 0.020 0.025 0.030 0.035 St res s ( M P a ) 0 10 20 30 40 50 25 ? C 50 ? C 75 ? C 125 ? C 150 ? C 100 ? C Figure 3.4: Temperature Dependent Average Stress-Strain Curves of Underfill NUF1 (? & = 0.001 sec -1 ) 26 0 1 2 3 4 5 6 0 25 50 75 100 125 150 175 200 Temperature, T (?C) E l a s tic M o d u lu s , E (G P a ) Figure 3.5: Temperature Dependent Elastic Modulus of NUF1 ( = 0.001 sec? & -1 ) 27 0 2 4 6 8 10 12 0 25 50 75 100 125 150 175 200 Temperature, T (?C) E l as t i c M odul us , E ( G P a NUF1 UF1 Figure 3.6: Elastic Modulus of a Nano-Filler Underfill and a Micro Filler Underfill (? & = 0.001 sec -1 ) 28 29 by NASA. In this present study tests has been performed at very extreme low temperatures down to -175 ?C. A newly developed environmental chamber is used with the MT-200 testing system in this purpose. The environmental chamber shown in Figure 3.7 is used for the cryogenic tests and it generates the cooling effect by re-circulating cold nitrogen gas from a liquid nitrogen source. Figure 3.8 shows stress strain curves at -50, - 100 and -175 ?C. The measured data have been plotted together with the room temperature and higher stress-strain curves shown earlier in Figure 3.4. The test data shows that underfill NUF1 became more linear elastic as the temperature decreases to cryogenic temperature. Figure 3.9 shows the elastic modulus of underfill NUF1 at extreme low temperatures. The elastic modulus also increases linearly with decrease in temperature from the room temperature to negative temperatures. At -175 ?C elastic modulus of NUF1 became almost double than the room temperature value. Preliminary Creep Data Preliminary creep curves for the underfill NUF1 are shown in Figure 3.10. All the tests were performed at a constant stress level of 10 MPa. Figure 3.10 illustrates the strong influence of temperature upon the deformation under constant load. At high temperatures (near the T g ) the creep compliance is greatly increased. Preliminary Relaxation Data Preliminary stress relaxation test curves for the underfill NUF1 are shown in Figure 3.11. All the tests are performed for a constant strain level which is 1% in this case. From test data it is observed that the stress relaxation rate increases with the Figure 3.7: Cryogenic Test Chamber for MT 200 Testing System 30 Strain 0.000 0.005 0.010 0.015 0.020 0.025 0.030 0.035 St ress (M Pa) 0 20 40 60 80 25 ? C 31 50 ? C 75 ? C 125 ? C 150 ? C -50 ? C -100 ? C -175 ? C 100 ? C Figure 3.8: Temperature Dependent Stress Strain Curves at Extreme Low Temperatures (? & = 0.001 sec -1 ) 0 2 4 6 8 10 -200 -150 -100 -50 0 50 100 150 200 Temperature, T (?C) El a s t i c M o d u l u s , E ( G P a ) .. . Figure 3.9: Temperature Dependent Elastic Modulus in Cryogenic Temperatures (?& = 0.001 sec -1 ) 32 Time (min) 0 20 40 60 80 100 120 C r ee p S t r a in 0.000 0.005 0.010 0.015 0.020 0.025 0.030 125 ? C 33 100 ? C 75 ? C 25 ? C 150 ? C Figure 3.10: Temperature Dependent Creep Data of NUF1 Time (min) 0 20 40 60 80 100 120 St ress (M Pa) 0 5 10 15 20 25 30 35 25 ? C 100 ? C 125 ? C 75 ? C Figure 3.11: Temperature Dependent Stress Relaxation Test Data of Underfill NUF1 34 increase of temperature. All the tests were stopped after 100 minutes as it almost reaches to a very smaller rate of change of stress after this time. 3.1.3: Uniaxial Testing of Nano Underfill 2 (NUF2) Uniaxial testing specimens of the nano-underfill were developed. The specimen preparation procedure is presented elsewhere [19]. With the developed test specimen; tensile and creep tests have been performed in a wide temperature range. The nano- underfill material for detailed material testing is NUF2, and has a cure time of 30 min at 150?C. The underfill has a volume fraction in the neighborhood of ? = 0.11. The glass transition temperature of this material is 98.24 ?C (DMA Method). The thickness of the cured uniaxial specimen was 75-125 ?m (3-5 mil). The typical length and width of a specimen are 90 mm and 3 mm respectively. In all uniaxial tests, the effective test length of the uniaxial specimen was 60 mm. Stress-Strain Data Typical variation of the stress-strain curves of the tested underfill with temperature is shown in Figure 3.12. The strain rate for these tests was = 0.001 sec?& -1 . The total test time (to failure) of a typical tensile test is less than 5 seconds for room temperature experiments. Tests were performed at T = 25, 50, 75, 100, 125 and 150?C. The stress-strain curves in Figure 3.12 illustrate considerable softening and viscoplastic type behavior as the temperature is increased. Elastic modulus decreases approximately linearly with the increase of temperature from 25 to 90 ?C. Figure 3.13 shows temperature dependent elastic modulus of underfill 35 Strain, ? 0.00 0.02 0.04 0.06 0.08 0.10 S t re ss, ? (M P a ) 0 5 10 15 20 25 30 35 25 ? C 50 ? C 75 ? C 90 ? C 100 ? C 125 ? C 150 ? C Figure 3.12: Temperature Dependent Stress Strain Curves (? = 0.001 sec & -1 ) 36 Temperature, t ( ? C) 0 25 50 75 100 125 150 175 El as tic M o d u lu s , E ( G Pa ) 0.0 0.5 1.0 1.5 2.0 2.5 3.0 3.5 Figure 3.13: Temperature Dependent Elastic Modulus (? & = 0.001 sec -1 ) 37 38 NUF2. The value of elastic modulus at 25?C is in the neighborhood of 2.9 GPa. This correlates well with the predicted value of 3.25 GPa from the unit cell model [21]. At temperature, T = 125 ?C elastic modulus is almost near to zero. This is typical as underfill temperature exceeds its glass transition temperature. Figure 3.14 shows a comparison of elastic modulus of the nano-underfill NUF1, NUF2 versus micron-filler underfill UF1. UF1 is has micro size particles as fillers. The volume fraction of filler particle in UF1 is almost double than NUF1 and almost four times than NUF2. Figure 3.14 shows that due to higher filler concentration UF1 has higher elastic modulus than NUF1 and NUF2. Glass transition temperature of UF1 is 150 ?C and its modulus drops greatly after 100 ?C but NUF1 has glass transition temperature as 155.86 ?C and its modulus drops greatly after 125 ?C. As glass transition temperature of NUF2 is 98.24 ?C its property degrades a lot after 90 ?C. Preliminary Creep Data Preliminary creep curves for the underfill NUF2 are shown in Figure 3.15. All the tests were performed at a constant stress level of 10 MPa. Figure 3.15 illustrates the strong influence of temperature upon the deformation under constant load. At high temperatures (near the T g ) the creep compliance is greatly increased 3.1.4: CTE Measurement Coefficient of thermal expansion of underfill is the most critical thermo- mechanical property of underfill. Accurate measurement of CTE of the cured underfill is challenging. In this present work strain gage method was applied to measure the CTE of 0 2 4 6 8 10 12 0 25 50 75 100 125 150 175 200 Temperature, T (?C) E l ast i c M odu l u s, E (G P a ) AAA-01 39 NUF2 NUF1 E l ast i c M odu l u s, E (G P a ) Figure 3.14: Elastic Modulus of Nano-Filler Underfills and a Micro Filler Underfill (? & = 0.001 sec -1 ) Time, t (min) 0 20 40 60 80 100 120 Cre ep S t rain 0.000 0.002 0.004 0.006 0.008 0.010 0.012 0.014 0.016 0.018 90 ? C 75 ? C 50 ? C 25 ? C Figure 3.15: Temperature Dependent Creep Data of NUF2 40 cured underfill. Strain gage technique is a very accurate method for CTE measurement [29]. Hence the specimen for this method should be very stiff, in this case much thicker underfill specimens than uniaxial test specimens were developed. A 25 mm x 15 mm x 1 mm specimen was cast and cured at 150?C for 30 minutes. A strain gage was placed on the cured specimen. A special type of CTE measurement strain gage was chosen for this purpose. A similar strain gage was placed on a reference material. The reference material was a Titanium-Silicate Bar (TSB) bar specially made for using in CTE measurement. Figure 3.16 shows the photograph of the underfill specimen with strain gage. Figure 3.17 shows strain gage on reference material. The underfill specimen and the TSB bar were heated in an oven from room temperature to the 120 ?C of the underfill, NUF1. Figure 3.18 and Figure 3.19 shows the strain gage readings at different temperatures for the gages placed on the reference material and underfill specimen respectively. Figure 3.20 shows the thermal strain developed in the underfill specimen at different temperatures. The CTE was calculated by using the following equation [29], T EE RS RS ? ? =??? (3.4) where and are CTE of the test specimen and reference material respectively. E S ? R ? S and E R are strain gage output of the gages on the test specimen and the reference material respectively and ?T is temperature difference. The calculated value of CTE of the underfill NUF1 is 39 ppm/?C. The value correlates reasonably with the predicted CTE value of 45 ppm/?C from the unit cell model presented later in chapter 6. 41 Figure 3.16: Underfill Specimen with Strain Gage 42 Figure 3.17: Strain Gage on Reference Material. (Titanium-Silicate Bar,TSB) 43 -2500 -2000 -1500 -1000 -500 0 0 2040608010120140 Temperature (?C) S t r a i n G a ge R eadi ng .... Figure 3.18: Strain Gage Reading from the Gage on the Reference Material. 44 0 200 400 600 800 1000 1200 1400 1600 0 20 40 60 80 100 120 140 Temperature (?C) S t r a i n G a ge R eadi ng .... Figure 3.19: Strain Gage Reading from the Gage on the Test Specimen 45 0 0.0008 0.0016 0.0024 0.0032 0.004 0 25 50 75 100 125 150 Temperature (?C) T h e r m a l S t r a in .. Figure 3.20: Thermal Strain on Underfill Specimen at Different Temperature 46 3.2 Characterization of Macro Silica Underfill Reliable, consistent, and comprehensive material property data are needed for underfills for the purpose of mechanical design, reliability assessment, and process optimization of electronic packages. For example, it is common to use finite element models to simulate the reliability of flip chip on laminate assemblies subjected to thermal cycling (e.g. -40 to 125 ?C). In this section material characterization of several macro underfills are performed via thermomechanical testing. The six-axis micro tension/torsion machine has been used to perform tensile, creep tests. 3.2.1: Effect of Extreme Low Temperatures Accurate evaluation of the stress-strain behavior of underfill materials at extreme of low temperatures is very important for applications in space environments such as planned for several upcoming NASA solar system exploration missions. There have been very few studies on the low temperature behavior of microelectronic encapsulants. In the present section, the stress-strain testing has been extended down to -185 ?C by using the newly developed extreme cold temperature test chamber shown in Figure 3.7. With the developed system shown in Figure 3.7, stress-strain and creep tests have been performed on several underfill materials down to cryogenic temperatures. For example, stress-strain curves of underfill encapsulant UF1 (curing condition and basic properties listed in Table 3.1) are shown in Figure 3.21 for the temperature range of -175 to +150 o C. For these experiments, the strain rate was ?& = 0.001 sec -1 . The measured elastic moduli of the same underfill encapsulant are plotted over the same temperature range in Figure 3.22. As can be seen from the data, the stress-strain curves become linear elastic 47 Table 3.1: Basic Underfill Properties as Provided by Vendor Underfill Material g T ? (1/ o C) (Below Tg) Filler Particle Diameter Filler Content Recommended Cure Conditions UF1 150 o C 25 x 10 -6 1-10 ?m 50% 150 o C / 30 min 48 Strain, ? 0.000 0.004 0.008 0.012 0.016 0.020 Stress, ? (MPa) 0 20 40 60 80 100 120 -175 ? C -125 ? C -100 ? C -75 ? -25 ? C 0 ? C +25 ? C +50 ? C +75 ? C +100 ? C +125 ? C +150 ? C C -50 ? C Figure 3.21 - Underfill Stress-Strain Curves for Material UF2 (T = -175 to 150 o C, ?& = 0.001 sec -1 ) 49 0 5 10 15 20 25 -200 -175 -150 -125 -100 -75 -50 -25 0 25 50 75 100 125 150 175 Temperature, T (?C) E l as t i c M odu l u s , E ( G P a ) Figure 3.22 - Elastic Modulus vs. Temperature for Material UF2 (T = -175 to 150 o C, ?& = 0.001 sec -1 ) 50 51 with much higher stiffness as the temperature decreases. With constant load tests, we have verified that creep is negligible for temperatures at or below T = -50 o C. The material is highly viscoplastic at room temperature and above. At the cryogenic temperature T = -175 o C, the value of the elastic modulus is almost double that of the room temperature value. For temperatures from -175 to +100 o C, the elastic modulus decreases linearly with a constant slope. Above T = +100 o C, the elastic modulus decreases at a faster rate as the material enters the glass transition region. The vendor reported the T g of this material is +150 o C (DMA method). 3.2.2: Temperature Dependent Creep Curves Complete characterization of the uniaxial constitutive behavior of underfills requires the measurements of time dependent creep or stress relaxation curves. Temperature dependent preliminary creep curves of UF2 for the underfill under consideration are shown in Figures 3.23. All the tests were done at 20 MPa. Figure 25 illustrates the strong influence of temperature. At high temperatures (near the Tg) the creep compliance is greatly increased. Time (Sec) 0 2000 4000 6000 8000 10000 Cr eep St r a in 0.000 0.005 0.010 0.015 0.020 0.025 0.030 25 ? C 110 ? C 90 ? C 100 ? C 75 ? C Figure 3.23: Temperature Dependent Underfill Creep Data (Stress Level 20 MPa) 52 53 CHAPTER 4 FINITE ELEMENT SIMULATION OF UNIAXIAL TEST SPECIMEN In this chapter detail Finite Element Analysis (FEA) simulation of uniaxial tensile tests has been performed by using a general purpose software ANSYS. Figure 4.1 shows a general description of the model geometry, boundary conditions and loading. Due to the symmetry, half of the uniaxial specimen has been modeled. The specimen chosen for this study has length to width aspect ratio 40, length 120 mm and width 3 mm. For the FEA study half of the specimen has been modeled with length 60 mm and width 3 mm. The actual tests are displacement controlled and one end of the specimen is attached with the stationary stage of the tester and the other end of the specimen is attached with the moving stage. In Figure 4.1 the left end of the specimen represents the stationary stage of the tester and fixed boundary condition is applied at the left end of the geometry. The right end represents the middle of the specimen and prescribed displacement is applied here. Figure 4.2 shows the details of the FEA mesh used in the model. In this FEA study the measured temperature dependent stress-strain curves from the real tensile tests have been used as input material property. 4.1: Effect of Length to Width Aspect Ratio For uniaxial tests, stress uniformity throughout the specimen is a critical factor. Earlier it was stated that desired value for length to width aspect ratio is 10 - 20 to yield a Prescribed Displacement Fixed Boundary Condition Figure 4.1: FEA Model of Symmetric Half of the Tensile Specimen. W Aspect Ratio = 40. 54 idth to Length Figure 4.2: FEA Mesh of the Model 55 56 reasonably pure uniaxial stress state in the specimen. Finite element analysis has been performed with the aspect ratios 2, 6, 20 and 40. The analysis temperature was 25 ?C for this part of the study. Figure 4.3, 4.4, 4.5, 4.6 represents the results for the aspect ratio 2. The FEA mesh of the model for aspect ratio 2 is shown in Figure 4.3. The element size was chosen as 0.5 mm X 0.5 mm square and kept same for all the models with different aspect ratios. Figure 4.4 shows the axial stress distribution along the length of the specimen. Stress in this case is not uniform at all throughout the whole specimen. This large nonuniformity in the axial stress value demonstrates that the length to width aspect ratio 2 is a worst choice for uniaxial test specimen. At this case value of stress at the rightmost vertical path of the model is not also constant, Figure 4.5. As symmetric half of the specimen was modeled this rightmost vertical path represents the middle of the specimen where stress must have a constant value in a good acceptable uniaxial tensile test. Figure 4.6 shows a comparison between the input experimental stress-strain curve and the stress-strain curve generated by the FEA output. The generated stress-strain curve from FEA output is much different than the input experimental stress-strain curve. This also supports the idea that a very small length to width aspect ratio like 2 is a very poor choice for the uniaxial tensile specimen. Figure 4.7 to 4.10 shows the simulation results for the aspect ratio 6. Axial stress distribution in Figure 4.8 shows that stress is not uniform through a large portion of the specimen. Non uniformity of the stress spreads towards approximately 40% of the total length of the specimen. The non uniformity of the stress starts at the edge of the specimen where fixed boundary condition is applied and spreads towards the center of the specimen. At the edge, the application of the fixed boundary condition prevents the Figure 4.3: FEA Mesh of the Model. Specimen Length, L = 6 mm, Width, W = 3 mm and Length to Width Aspect Ratio = 2 57 Figure 4.4: Axial Stress Distribution (? x ). Length to Width Aspect Ratio = 2 58 3.00E+07 3.10E+07 3.20E+07 3.30E+07 3.40E+07 3.50E+07 0 0.0005 0.001 0.0015 0.002 0.0025 0.003 0.0035 Vertical Position (m) [zero at top] St re s s (Pa ) Stress Along Rightmost Vertical Path Figure 4.5: Axial Stress (? x ) Value along the Right Most Vertical Path. Length to Width Aspect Ratio = 2 59 0 10 20 30 40 0 0.0005 0.001 0.0015 0.002 0.0025 0.003 0.0035 0.004 Strain Stress (MPa) Input Experimental Data FEA Output Figure 4.6: Comparison of FEA Output with the Input Experimental Data, Length to Width Aspect Ratio = 2 60 Figure 4.7: FEA Mesh of the Model. Specimen Length, L = 18 mm, Width, W = 3 mm and Length to Width Aspect Ratio = 6 61 Figure 4.8: Axial Stress Distribution (? x ). Length to Width Aspect Ratio = 6 62 3.00E+07 3.10E+07 3.20E+07 3.30E+07 3.40E+07 3.50E+07 0 0.0005 0.001 0.0015 0.002 0.0025 0.003 0.0035 Vertical Position (m) [zero at top] Stre s s (Pa ) Stress Along Rightmost Vertical Path Figure 4.9: Axial Stress (? x ) Value along the Right Most Vertical Path. Length to Width Aspect Ratio = 6 63 0 10 20 30 40 0 0.0005 0.001 0.0015 0.002 0.0025 0.003 0.0035 0.004 Strain S t re s s (MP a ) Input Experimental Data FEA Output Figure 4.10: Comparison of FEA Output with the Input Experimental Data. Length to Width Aspect Ratio = 6 64 65 expected contraction due to the poisons ratio effect in the lateral direction and aids in the generation of the stress non uniformity. During the real experiment the edge of the specimen also can not deform laterally due to the grip, hence Figure 4.8 represents the stress distribution of a real tensile test. Uniaxial specimen with length to width aspect ratio 6 will be also not acceptable as a good test specimen. The value of stress at the rightmost vertical path is almost constant, Figure 4.9. Figure 4.10 shows that generated stress-strain curve by FEA output is almost same as the experimental stress-strain curve. The results for the aspect ratio 20 and 40 are summarized in Figure 4.11 to 4.14 and Figure 4.15 to 4.18 respectively. Both of these cases have very little non uniformity in axial stress distribution. Non uniformity of the stress spreads towards approximately 5% of the total length of the specimen in both cases. In both cases stress at the rightmost vertical edge have constant value and stress-strain curves from experimental input data and FEA output are exactly same in both cases. Figure 4.19 shows a plot demonstrating the change in elastic modulus value with the change in aspect ratio. The elastic modulus values in this plot are the initial slopes of the stress-strain curves generated by the FEA output data. The stress-strain curves were measured at the rightmost end of the model geometry which represents the center of the uniaxial specimen. The input material properties were same in all cases. This plot shows that at first the elastic modulus decreases with the increase of length to width aspect ratio of a thin uniaxial tensile specimen and then remains almost constant for any length to width aspect ratio which is 20 or more. Figure 4.20 shows the change in percent difference in modulus value with the increase in aspect ratio. Here the modulus value for the specimen with length to width aspect ratio 40 was taken as the base value. These results demonstrate that uniaxial Figure 4.11: FEA Mesh of the Model. Specimen Length, L = 60 mm, Width, W = 3 mm and Length to Width Aspect Ratio = 20 66 Figure 4.12: Axial Stress Distribution (? x ). Length to Width Aspect Ratio = 20 67 3.00E+07 3.10E+07 3.20E+07 3.30E+07 3.40E+07 3.50E+07 0.0000 0.0005 0.0010 0.0015 0.0020 0.0025 0.0030 0.0035 Vertical Position (m) [zero at top] Stre s s (Pa ) Stress Along Rightmost Vertical Path Figure 4.13: Axial Stress (? x ) Value Along the Right Most Vertical Path. Length to Width Aspect Ratio = 20 68 0 10 20 30 40 0 0.0005 0.001 0.0015 0.002 0.0025 0.003 0.0035 0.004 Strain S t re s s (MP a ) Input Experimental Data FEA Output Figure 4.14: Comparison of FEA Output with the Input Experimental Data. Length to Width Aspect Ratio = 20 69 Figure 4.15: FEA Mesh of the Model. Specimen Length, L = 120 mm, Width, W = 3 mm and Length to Width Aspect Ratio = 40 70 Figure 4.16: Axial Stress Distribution (? x ). Length to Width Aspect Ratio = 40 71 3.00E+07 3.10E+07 3.20E+07 3.30E+07 3.40E+07 3.50E+07 0 0.0005 0.001 0.0015 0.002 0.0025 0.003 0.0035 Vertical Position (m) [zero at top] Stre s s (Pa ) Stress Along Rightmost Vertical Path Figure 4.17: Axial Stress (? x ) Value along the Right Most Vertical Path. Length to Width Aspect Ratio = 40 72 0 10 20 30 40 0 0.0005 0.001 0.0015 0.002 0.0025 0.003 0.0035 0.004 Strain Stress (MPa) Input Experimental Data FEA Output 73 Figure 4.18: Comparison of FEA Output with the Input Experimental Data. Length to Width Aspect Ratio = 40 10.2 10.3 10.4 10.5 10.6 10.7 10.8 0 5 10 15 20 25 30 35 40 45 Length to Width Aspect Ratio E l as t i c M odul us by F E A D a t a ( G P a ) . Figure 4.19: Effect of Length to Width Aspect Ratio on Elastic Modulus. 74 0 0.5 1 1.5 2 2.5 0 5 10 15 20 25 30 35 40 45 Length to Width Aspect Ratio % D i f f e r e nce i n V a l u e o f E l ast i c M o d u l u s . Figure 4.20: Percent Difference in Elastic Modulus Values from the Elastic Modulus Obtained for Length to Width Aspect Ratio 40 75 76 specimen with a large length to width aspect ratio like 20 or more is a very good choice for the uniaxial tensile tests. Again for thin film testing too long specimen may produce erroneous results due to the bending effect by self weight while a horizontal test setup like the thermo-mechanical test system of this present study is been used. Considering all of this, length to width aspect ratio 20 will be the best choice for a uniaxial tensile specimen. Checking for Stress Uniformity The axial stress value along some longitudinal paths of the model has been plotted to check the stress uniformity through out the specimen. This study has been done with the model of length to width aspect ratio 40. Figure 4.21 to 4.23 shows that stress is constant along the longitudinal paths (horizontal paths in the figures) of the specimen. 4.2: Nonlinear Simulation of Temperature Dependent Stress Strain Curves In this part detailed comparison of temperature dependent stress-strain curves generated by FEA and experimental data has been performed. FEA model with length to width aspect ratio 40 has been chosen for this purpose. The simulations are performed at temperature 125 ?C. Figure 4.24 shows that experimental stress-strain curves at all temperatures are exactly same as the stress-strain curves generated by FEA output. Due to this nice similarity between the experimental and FEA stress-strain curves it will be possible to make FEA prediction of the stress-strain curves at some temperature where no test has been done. Figure 4.25 shows such a prediction at 110 ?C. 0.00E+00 2.00E+06 4.00E+06 6.00E+06 8.00E+06 1.00E+07 1.20E+07 1.40E+07 1.60E+07 1.80E+07 2.00E+07 0 0.01 0.02 0.03 0.04 0.05 0.06 0.07 Position Along Top Horizontal Line (m) [Zero at Left] X - S t re s s (P a ) Stress Along Top Horizontal Path Figure 4.21: Axial Stress along Top Horizontal Path. Length to Width Aspect Ratio = 40 77 0.00E+00 2.00E+06 4.00E+06 6.00E+06 8.00E+06 1.00E+07 1.20E+07 1.40E+07 1.60E+07 1.80E+07 2.00E+07 0 0.01 0.02 0.03 0.04 0.05 0.06 0.07 Position Along Mid Horizontal Line (m) [Zero at Left] X - St r e ss (Pa ) Stress Along Mid Horizontal Path Figure 4.22: Axial Stress along Mid Horizontal Path. Length to Width Aspect Ratio = 40 78 0.00E+00 2.00E+06 4.00E+06 6.00E+06 8.00E+06 1.00E+07 1.20E+07 1.40E+07 1.60E+07 1.80E+07 2.00E+07 0 0.01 0.02 0.03 0.04 0.05 0.06 0.07 Position Along Bottom Horizontal Line (m) [Zero at Left] X - S t re s s (P a ) Stress Along Bottom Horizontal Path Figure 4.23: Axial Stress along Bottom Horizontal Path. Length to Width Aspect Ratio=40 79 0 5 10 15 20 25 30 35 40 0 0.002 0.004 0.006 0.008 0.01 0.012 0.014 Strain St re s s (M Pa ) FEA Output Input Experimental Data Si 3 25 ?C 50 ?C 75 ?C 100 ?C 125 ?C 150 ?C Figure 4.24: Comparison of Stress-Strain Curves Generated by FEA with Input Experimental Data. Length to Width Aspect Ratio = 40 80 0 5 10 15 20 25 30 35 40 0 0.002 0.004 0.006 0.008 0.01 0.012 0.014 Strain St r e s s ( M Pa ) Experimenta Input Datal FEA Output 125 C (E i t l) Figure 4.25: FEA Prediction of Stress-Strain Curves at 110 ?C. Length to Width Aspect Ratio = 40 81 4.3: Error Estimation by Multi-mesh Extrapolation Error in axial stress value due to mesh refinement is evaluated by Multimesh extrapolation [14]. Figure 4.26 to 4.28 shows mesh refinement. In each refinement, nodes and interelement boundaries of the coarser mesh are preserved, while new nodes, elements, and interelement boundaries are added. A plot of axial stress vs. h q is shown in Figure 4.29. Here q = 1 and h is defined by eq. (4.1). 82 )(ElementofNumber 1 h = (4.1) Calculations are done for the mid node at the rightmost vertical path. Value of h and Axial Stress for Different Number of Element are given in Table 4.1. Calculated error in axial stress value for the model with 720 elements is found as 0.169%. This small error ensures that mesh picked for the models discussed before are refine enough. These Nodal Positions will be Preserved During Refinement Figure 4.26: Original Mesh Before Refinement. Total Number of Element = 720. Length to Width Aspect Ratio = 40 83 Mesh Refinement Preserved Nodal Positions from Coarser Mesh Figure 4.27: Mesh Refinement. Total Number of Element = 2880. Length to Width Aspect Ratio = 40 84 Preserved Nodal Positions from Coarser Mesh More Mesh Refinement Figure 4.28: More Mesh Refinement. Total Number of Element = 11520. Length to Width Aspect Ratio = 40 85 Stress = -818784 h + 16583000 0 2500000 5000000 7500000 10000000 12500000 15000000 17500000 20000000 0 0.005 0.01 0.015 0.02 0.025 0.03 0.035 0.04 h Axial Stress (Pa) Figure 4.29: Axial Stress vs. h Plot. Calculations are done for the Mid Node of the Rightmost Vertical Edge 86 Table 4.1: Value of h and Axial Stress for Different Number of Element 165760000.0093169511520 165668000.01863392880 165528000.0372678720 Axial Stress (Pa)h Number of Elements 87 88 CHAPTER 5 PREDICTION OF MATERIAL PROPERTIES MACRO-SILICA UNDERFILLS 5.1: Underfill Constituent Properties of Macro-Silica Underfills Macro silica underfills are formed by dispersing micron sized fillers in epoxy matrix. In this chapter Finite Element Models have been generated to predict coefficient of thermal expansion and elastic modulus of macro silica underfills. The property prediction has been performed based on the constituent properties of underfill. Underfill constituent properties are given in Table 5.1. The underfill material chosen for the prediction study has 49% filler particles by volume. 5.2: Prediction of Material Properties of Macro Silica Underfills 5.2.1: Unit Cell Generation Unit cells are generated to predict the underfill material properties by symmetric distribution of constant sized spherical fillers in epoxy matrix. In real underfill randomly sized fillers particles are distributed randomly in an epoxy matrix. Unit cells with symmetrically distributed and constant sized spherical filler particles are considered to avoid complexity of the finite element models. Unit cells with random sized and randomly distributed filler particles will be discussed in chapter six. Two types of unit cells are developed. First type of unit cell has one filler at the center of an epoxy cube. 89 Table 5.1 Material Properties for Epoxy Matrix and Micro-Silica Particle Elastic Modulus (GPa) Poisson Ratio CTE (ppm/?C) Filler 77.8 0.19 0.5 Epoxy 2.5 0.40 62.46 The length of the epoxy cube is 4.6 ?m and the diameter of the filler particle is 4.5 ?m. The volume of the cube is 97.34 ?m 3 and the volume of the filler particle is 47.71 ?m 3 .The volume fraction of the filler is 49% and was calculated by equation (5.1). The CubetheofVolume FillertheofVolume FractionVolume =?, (5.1) other unit cell has total eight fillers distributed symmetrically. The length of the cube is 9.2 ?m and the diameter of the particle is 4.5 ?m. The volume fraction of the filler was 49% and also calculated by equation (5.1). Figure 5.1 and 5.2 shows the filler distribution in the one filler and eight filler unit cell respectively. 5.2.2: Prediction of Coefficient of Thermal Expansion by 1 Filler Model Finite element model was developed by the generated unit cells to predict the coefficient of thermal expansion. Figure 5.3 shows the finite element model of one filler cell. The origin is at the center of the unit cell. Symmetric boundary condition was applied as shown in Figure 5.4. The reference temperature was 25 ?C and the maximum temperature was applied as 125 ?C. Figure 5.5 shows the thermal expansion in x direction due to rise in the temperature from 25 ?C to 125 ?C. Figure 5.6 shows the surface plot of the displacement values of the nodes on the rightmost surface. Since the filler particle is more rigid than the epoxy, the nodes at the center have less displacement and the nodes away from the center have more displacement. This results a cone shape 90 Epoxy Cube Side Length = 4.6 ?m CTE = 62.46 X 10 -6 /?C Volume = 97.336 m 3 Filler Particle Diameter = 4.5 ?m CTE = 0.5 X 10 -6 /?C Volume = 47.713 m 3 Figure 5.1: One Filler Unit Cell 91 Epoxy Cube Side Length = 9.2 ?m CTE = 0.5 X 10 -6 /?C Volume = 778.688 ?m3 Filler Particle Diameter = 4.5 ?m CTE = 62.46 X 10 -6 /?C Volume = 8 X 47.713 ?m 3 Figure 5.2: Eight Filler Unit Cell 92 Figure 5.3: Finite Element Model of One Filler Unit Cell 93 Figure 5.4: Applied Boundary Condition in One Filler Unit Cell 94 Figure 5.5: Thermal Expansion in x Direction 95 Cg = 0.0068 ?m from Un-deformed Surface 0.000 0.002 0.004 0.006 0.008 0.010 -3 -2 -1 0 1 2 3 -3 -2 -1 0 1 2 D is pl ac e m e nt ( ? m ) Y C o - o r d i n a t e ( ? m ) X C o -o rd in a te (? m ) Maximum Displacement = 0.009216 ?m z y x Figure 5.6: U x Displacement of the Deformed Surface at X = 2.3 ?m from Center [Right Most Surface] 96 97 surface after deflection. The displacements are symmetric about the center axes as the cube is symmetric about that axes. The distance of the C g of the deflected surface from the neutral surface is 0.0068 ?m. Based on the C g the deflected surface the coefficient of thermal expansion of the underfill is UF1 is 29.53 ppm. The published value of the C g by the vendor is 25 ppm, which well correlates with the prediction. 5.2.3: Prediction of Coefficient of Thermal Expansion by 8 Filler Model Figure 5.7 shows the finite element model of one filler cell. The origin is at the center of the unit cell. Symmetric boundary condition was applied as shown in Figure 5.8. The reference temperature was 25 ?C and the maximum temperature was applied as 125 ?C. Figure 5.9 shows the thermal expansion in x direction due to rise in the temperature from 25 ?C to 125 ?C. Figure 5.10 shows the surface plot of the displacement values of the nodes on the rightmost surface. Since the filler particle is more rigid than the epoxy, the nodes near the filler particles have less displacement and the nodes away from the filler particles have more displacement. This results the shape of the deflected surface shown in Figure 5.10. The displacements are symmetric about the center axes as the cube is symmetric about that axes. The maximum deflection occurs at the center of the right most surface as there are no filler at the center of the eight filler unit cell. The distance of the C g of the deflected surface from the neutral surface is 0.013 ?m. Based on the C g the deflected surface the coefficient of thermal expansion of the underfill is UF1 is 28.81 ppm. The published value of the C g by the vendor is 25 ppm, which well correlates with the prediction. Figure 5.7: Finite Element Model of Eight Filler Unit Cell 98 Figure 5.8: Applied Boundary Condition in Eight Filler Unit Cell 99 Figure 5.9: Thermal Expansion in x Direction 100 0.008 0.010 0.012 0.014 0.016 0.018 0.020 0.022 -6 -4 -2 0 2 4 6 -6 -4 -2 0 2 4 D i s p la c e m en t ( ? m ) Y C o - o r d i n a t e ( ? m ) Z C o -o rd in a te (? m ) Cg = 0.013254 ?m from Un-deformed Surface Maximum Displacement = 0.020915 ?m z y x Figure 5.10: U x Displacement of the Deformed Surface at x = 4.6 ?m from Center [Right Most Surface] 101 102 5.2.4: Prediction of Elastic Modulus by 1 Filler Model Finite element models have also developed to predict the elastic modulus of the micro-silica underfill. Figure 5.11 shows the applied boundary conditions in the model. The origin is at the center of the cube. The x displacement is zero at x equal to zero, the z displacement is zero at z equal to zero and the y displacement is zero at y equal to negative 2.3 ?m. The fixed boundary condition is applied at the center node of the bottom surface. The displacement is coupled in the y direction at the top surface at y equal to 2.3 ?m. A compressive pressure load is applied on the nodes at the top surface at y equal to 2.3 ?m. The load was applied in a ramped manner. Figure 5.12 shows the y deflection contour in the unit cell. The stress strain curve was calculated at the top surface of the unit cell and shown in Figure 5.13. Only linear elastic properties for both epoxy and filler particles were entered as input material models. This causes the calculated stress strain curve to be linear elastic all the way. The predicted elastic modulus was 13.64 GPa, which correlates well with the measured value of the elastic modulus 10.43 GPa. 5.2.5: Prediction of Elastic Modulus by 8 Filler Model Elastic modulus was also predicted by finite element analysis of the eight filler unit cell. Figure 5.14 shows the applied boundary condition. The origin is at the center of the cube. The x displacement is zero at x equal to zero, the z displacement is zero at z equal to zero and the y displacement is zero at y equal to negative 4.6 ?m. Fixed boundary condition was applied at the center node of the bottom surface. The displacement is coupled in the y direction at the top surface at y equal to 4.6 ?m. x z y U x = 0 @ x = 0 U z = 0 @ z = 0 U y = 0 @ y = -2.3 ?m Figure 5.11: Applied Boundary Condition 103 Figure 5.12: Displacement in y Direction 104 ? = 13639 ? + 1E-06 0 5 10 15 20 25 30 35 40 0 0.0005 0.001 0.0015 0.002 0.0025 0.003 Strain A v er ag e E l em en t a l S t r e s s ( M P a ) Figure 5.13: Stress Strain Curve at the Top Surface of the Unit Cell 105 x z y U x = 0 @ x= 0 U y = 0 @ y= - 4.6 ?m U z = 0 @ z = 0 Figure 5.14: Applied Boundary Condition in Eight Filler Model 106 107 A compressive pressure load is applied on the nodes at the top surface at y equal to 2.3 ?m. The load was ramped from zero to its maximum value. Figure 5.15 shows the y deflection contour in the unit cell. The stress strain curve was calculated at the top surface of the unit cell and shown in Figure 5.16. Only linear elastic properties for both epoxy and filler particles were entered as input material models. This causes the calculated stress strain curve to be linear elastic all the way. The predicted elastic modulus was 11.33 GPa, which correlates well with the measured value of the elastic modulus 10.43 GPa. The predicted elastic modulus value from the eight cell model is more close to the measured value of the elastic modulus. Figure 5.15: Displacement in y Direction 108 ? = 11325 ? + 9E-06 0 5 10 15 20 25 30 35 40 0 0.0005 0.001 0.0015 0.002 0.0025 0.003 Strain A v er ag e E l em en t a l S t r e s s ( M P a ) Figure 5.16: Stress Strain Curve Calculated at the top Surface of the Unit Cell 109 110 CHAPTER 6 PREDICTION OF MATERIAL PROPERTIES OF NANO-SILICA UNDERFILLS BY CONSTANT SIZE FILLER MODELS 6.1 Underfill Constituent Properties of Nano-Silica Underfills Nano-silica underfills are particle reinforced composite materials in which nano silica fillers are distributed randomly in an epoxy matrix. In this present chapter methodologies for property prediction have been investigated. Models in this chapter is developed based on constituent component properties and able to predict the effective equivalent properties of statistically isotropic unit cell of nano-underfill formed by random distribution of spherical nano filler particles. The properties of constituent components are given in table 6.1. Since elastic modulus, shear modulus, bulks modulus, Poisson?s ratio, etc. are determined from the slope of initial linear portion of the corresponding curves, both epoxy and silica fillers are assumed to be linear elastic during predicting the above properties of the nano underfill. The epoxy was assumed viscoelastic and the silica fillers are assumed as linear elastic to develop the models to predict creep and relaxation behavior. In the present study an algorithm similar to random sequential adsortion algorithm was developed to generate statistically isotropic cubic unit cells of underfill containing up to 38% nano fillers. The models can be used for development of nano-silica underfills with desirable thermo-mechanical properties. 111 Table 6.1 Material Properties for Epoxy Matrix and Nano-Silica Particle Elastic Modulus (GPa) Poisson Ratio CTE (ppm/?C) Filler 77.8 0.19 0.5 Epoxy 2.5 0.40 62.46 6.2 Prediction of Material Properties of Nano Silica Underfills 6.2.1 Unit Cell Generation The modulus of elasticity and the coefficient of thermal expansion of nano underfills were predicted by the implicit finite-element models of three-dimensional cubic unit cells. The unit cell was generated by randomly distributing spherical fillers in an epoxy matrix. The volume of the cube is L 3 , N is the total number of particles and r is the radius of the spherical particle. Volume fraction (?) of the filler was determined as the ratio of total volume of the sphere to the volume of the cube, 3 3 3 4 L rN ? ? ? ? ? ? = ? ? (6.1) Volume fraction of the cube was controlled by varying the total sphere number N as required by the value of L. The radius of the sphere r is kept same for the analysis of all volume fractions. The fillers should be distributed in such way that the unit cell should be isotropic i.e. equivalent in all direction and it should be quite suitable for generating good finite element mesh. An algorithm based on modified random sequential adsortion (RSA) has been developed to generate the random center coordinates of the nano-silica particles in the underfill [13]. According to this algorithm all the accepted random coordinates of the particles (for ? = 0 to 0.25) pass the following conditions. (a) If the particle surface touches the surface of the cube, or if they are very close, it may not be possible to mesh or the generated finite element mesh will be 112 distorted or some time meshing may not be possible at all. To avoid these, the particles were kept inside of the cube at some minimum distance d 2 from the surface of the cube. r1.0rd 2 += (6.2) To fulfill the above condition, center coordinates of the i th particle must pass the following check 3,2,1jdLx 3,2,1jdx 2 i j 2 i j =?? =? (6.3) (b) If two adjacent particles overlap each other they will violate the rigid sphere condition. In addition, two adjacent particles can not touch each other. To fulfill this condition center coordinates of the i th particle must pass the following check r07.2d 1 ?= (6.4) ( ) )1i(.....,,1kdxx 1 ki ?=?? vv (6.5) The modified RSA method has been used to achieve volume fraction higher than ? = 0.25. In this study, the unit cells with volume fraction, ? , such that 0 < ? < 0.25, all the particle centers were kept inside the cube, no particle overlaps the outer surfaces of the cube. For unit cell volume fractions more than 0.25, the particles were allowed to 113 overlap the surface boundary of the cube. If any particle overlapped the surface, the portion of the particle outside the cube was carefully cut into several sections and copied at a suitable position on the opposite surface of the cube. First, a cubic cell with volume fraction less than 0.25 was generated by fulfilling the above two conditions and then the cell was compressed in to a smaller size while keeping the size of the spherical particle same. The size of the sphere and total number of the sphere does not change during the operation, but the cell is compressed, giving a higher volume fraction of filler content. Length of the cubic cell is compressed first, by multiplying it by a user defined shrinking factor c f (< 1). The new length of the cube is given by LcL fn ?= (6.6) The old position of the center of i th particle is given by i x r then new position of the particle will be given by i n x r f ii n cxx ?= rr (6.7) Once the particles have moved, they will be allowed to overlap the surface of the cube but no particle center will be outside of the cube. If the new coordinates of the i th particle falls outside of the cube then it moved back in the cube at a random position between the surface and at a inward distance (r-?). Where ? is a user defined constant and ? > 0.1r. i n x r 114 If the particle center sits at a distance r from the surface of the cube then the particle surface will touch the cube surface and meshing will not be possible. If the new coordinates of the ith particle i n x r falls inside of the cube and at a distance r from the surface then it will be moved at a random position between the surface and at a inward distance (r-?). Where ? is a user defined constant and ? > 0.1r If the particle overlaps any other previously accepted position then it will be moved to a new random position in a random direction. The new position will be given by ??= i n i n xx rr (6.8) where ? is a small random number. Smaller ? ensures faster convergence. The new position of the particle will be accepted only if all above conditions are fulfilled. Since the algorithm involves random movement and random positioning each of the iterations may not produce an acceptable distribution. The algorithm also counts the number of iterations and stops the program if number of iteration exceeds a certain number. 6.2.2 Isotropy of the Unit Cell Once a valid distribution of the filler particles has been created, the algorithm also calculates the centroid and moment of inertia of the distributed particles. These quantities are calculated to check the isotropy of the distribution. The distribution with centroids at a position in the neighborhood of 0.5L are accepted, where L is the length of 115 the side of the cube. An isotropic distribution will have identical moment-of-inertia for the nano-particles about three orthogonal axes. Figure 6.1 shows a plot of the coordinates of the centroids of the accepted distribution for different volume fractions. Table 6.2 shows the values of moment of inertia for the accepted distributions for several volume fractions. In almost all cases, moments of inertia about all three axes have very close values. This indicates that the generated unit cell is isotropic in all directions. The unit cell for volume fraction 0.38 has a smaller moment of inertia due to the compression of the cell Size. 6.2.3 Randomness of Filler Distribution Randomness of filler distribution also has been checked to ensure that no periodicity exists in the filler distribution. A radial distribution function of filler centroids, g(R), is used for this purpose [48, 49]. The function is given as dR )R(dK NR4 V )R(g 2 ? = (6.9) where V is the volume considered, N is the total number of filler particles in the considered volume, K(R) is the average number of filler centroids within a radius R from an arbitrary filler. Figure 6.2 shows a plot of g(R) for a distribution of filler volume fraction 0.20. The distributions are not periodic as very little presence of spikes of maxima and minima are seen. 116 Volume Fraction of Filler Particles 0.0 0.1 0.2 0.3 0.4 0.5 Pos i tion of Centr o id Along L 0.00 0.25 0.50 0.75 1.00 Position of X Centroid Position of Y Centroid Position of Z Centroid Figure 6.1: Centroids for Acceptable Nano-Particle Distribution 117 118 Table 6.2: Moment of Inertia of Acceptable Nano-Particle Distributions. Volume Fraction Moment of Inertia About Axis 1 Moment of Inertia About Axis 2 Moment of Inertia About Axis 3 0.10 100105 85312 86065 0.20 208298 195363 210952 0.25 278356 259267 254355 0.38 168768 179811 186471 r/d 12345678910 Radial Dis t ribution Func tion, g(R) 0 10 20 30 40 50 Figure 6.2: Plot of Radial Distribution Function 119 120 6.2.4 Finite Element Model of Unit Cell Finite element models of the cubic unit cells have been created to predict the CTE and elastic modulus of the nano underfill. In these models, spherical silica fillers were generated in an epoxy cube. The length of the epoxy cube is equal to the length of the unit cell considered in the unit cell algorithm. The filler particles were generated by using the center coordinates and radius. For the models having volume fraction of filler more than 0.25, fillers overlapping the cube surface were carefully cut into pieces and copied at a suitable place onto the other side of the cube. The model geometry was meshed by 8 node tetrahedral brick element. Figure 6.3 to Figure 6.5 shows representative models of the unit cell and filler distribution. More filler distribution Figure 6.6 shows example of the mesh generated in the FEA models. Figure 6.7 shows the isotropic distribution of particles with volume fraction, ? = 0.38. 6.2.5 Prediction of Coefficient of Thermal Expansion (CTE) The developed finite element model was used to predict the coefficient of thermal expansion of the nano-underfill for different volume fraction of filler particles. Symmetric boundary conditions were used at the cube faces at x, y, and z equal to zero. The degrees of freedom were coupled at the faces at x, y, z equal to L. The temperature of the model was raised to a user defined uniform temperature. This ensures that the cube faces will not be distorted after deformation. Calculated CTE values are presented in Figure 6.8 and Figure 6.9. In Figure 6.8, CTE in the direction of x-axis has been plotted for different volume fraction of the filler particles. Figure 6.8 shows that CTE of the underfill decreases linearly with the increase of volume fraction of the filler particles. Figure 6.3: Isotropic View of Filler Distribution, ? = 0.20 121 Figure 6.4: Front View of Filler Distribution, ? = 0.20 122 Figure 6.5: RHS View of Filler Distribution, ? = 0.20 123 Figure 6.6: FEA Mesh, ? = 0.20 124 Figure 6.7: Isotropic View of Filler Distribution, ? = 0.38 125 0 10 20 30 40 50 60 70 0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5 Volume Fraction of Filler Particle C T E ( ppm / ? C ) Figure 6.8: Prediction of x-CTE by finite-element analysis of Unit Cell 126 0 10 20 30 40 50 60 70 0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5 Volume Fraction of Filler Particle C T E ( ppm / ? C ) CTE in x Direction CTE in y Direction CTE in z Direction Figure 6.9: Prediction of x, y, z-CTE by finite-element analysis of Unit Cell 127 128 In Figure 6.9 calculated CTE values in x, y, z directions were plotted together on the same graph. This graph shows that CTE in all directions is almost identical which once again proves that the generated unit cells are isotropic. 6.2.6 Repeatability of Unit Cell All the results that presented in this chapter calculated from the analysis of one unit cell in each case. Due to the excellent repeatability of the results between different unit cells calculating results from one unit cell gives enough accuracy. In this section repeatability of calculated quantities from different unit cells have been checked. Five different unit cells with same volume fraction of fillers have been chosen. CTE has been calculated by the Finite Element Analysis of these five unit cells while same thermal load and boundary conditions have been applied. In each case obtained coefficient of thermal expansion value is almost exactly same. Results have been summarized in Table 6.3. 6.2.7 Prediction of Elastic Modulus The elastic modulus of nano underfill was also predicted for different volume fractions. For this analysis, different boundary conditions than the Coefficient of Thermal Expansion (CTE) models have been used. Displacement in the x-direction is zero at the cube surface at x equal to zero. Fixed boundary condition is applied at one of the point at the mid position of the surface at x equal to zero. Tensile load was applied at the surface at x equal to L in the form of uniform pressure. A suitable small value of the pressure load will give better results because the modulus is defined as the initial slope of the linear part of stress-strain curves. Figure 6.10 shows the elastic modulus vs. volume 129 Table 6.3: CTE Value from Several Different Unit Cell Model CTE (ppm/?C) Average CTE (ppm/?C) % Difference from Average CTE Unit Cell #1 47.1 0.1 Unit Cell #2 47.0 0.2 Unit Cell #3 46.8 0.6 Unit Cell #4 46.6 1.0 Unit Cell #5 48 47.1 2.0 0 2 4 6 8 10 12 0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5 Volume Fraction of Filler Particle El a s t i c M o d u lu s , E ( G Pa ) Figure 6.10: Prediction of Elastic Modulus by FEA Analysis of Unit Cell 130 131 fraction of filler particle. Elastic modulus increases exponentially with the increase of volume fraction of the fillers particles. 132 CHAPTER 7 PREDICTION OF MATERIAL PROPERTIES OF NANO-SILICA UNDERFILLS BY RANDOM SIZE FILLER MODELS 7.1 Unit Cell Generation In this chapter, random sequential adsorption has been used to develop finite- element models for nano-underfill unit cells. Widom [43] showed that a random configuration could be produced by random sequential addition of spheres volume subject to condition of no overlap. The addition may be achieved by random addition of particles or through perturbing a regular array of particles through random movements. [13, 44, 45, 46, 47]. Application of technique for analysis of time-dependent properties of nano- underfills investigated in this paper is new. Previously authors have investigated the application of RSA to prediction of elastic properties of nano-underfills [21]. In this chapter, volume fractions between 10 to 40 percent have been investigated, well short of the jamming limit of the particles. The unit cell of 10 percent filler volume fraction had 24 particles, 20 percent filler volume fraction had 48 particles, 25 percent filler volume fraction had 60 particles and 38 percent filler volume fraction had 48 particles. In real underfill materials all the filler particles do not have same diameter. In this study, particles sizes in unit cell have been randomly distributed in both size and location. Ability of using Random sequential adsorption for time-dependent property prediction in distribution of particle diameters has been investigated. The size distribution of the particles has been specified by a mean and standard deviation. Previous studies have focused on elastic property prediction with identical size particles [13, 21]. Bulk modulus, shear modulus and viscoelastic properties of nano underfills have been predicted by the implicit finite-element models of three-dimensional cubic unit cells. The unit cell was generated by distributing spherical fillers in an epoxy matrix. In real underfill all the spherical fillers do not have same diameter instead they follow a size distribution function. To simulate the real behavior of a nano-underfill, size distribution of filler particles is introduced. A new algorithm is developed which generates the random location of the filler center as well as a random diameter for the filler particles. The size randomization was performed by assuming a Gaussian distribution with ?3? . The median diameter was 50 nm and the minimum and maximum diameters were 34 and 66 respectively. The length of the new unit cell is 10 times than the median diameter of the filler particle. The volume of the cube is L 3 , N is the total number of particles and is the random radius of the spherical particle where l r ...3,2,1=l etc. The random radius was generated by using a built-in MATLAB function in such way that it follows the following conditions. l r max min rr andrr l l ? ? (7.1) After creating each successful sphere with random size and location the volume 133 134 )fraction ( N ? of the filler concentration was determined by the ratio of total volume of the sphere to the volume of the cube, 3 1 3 4 L r N l l N ? = ? ? ? ? ? ? = ? ? (7.2) The program stops executing when the volume fraction ( ) N ? becomes equal or slightly greater than the user defined volume fraction value ? i.e. the program stops when the following condition is passed. ?? ? N (7.3) The fillers should be distributed in such way that the unit cell should be isotropic i.e. equivalent in all direction and it should be quite suitable for generating good finite element mesh. An algorithm based on modified random sequential adsortion (RSA) [13] has been used to generate the random center coordinates of the nano-silica particles which have random diameters too. According to this algorithm all the accepted random coordinates of the particles (for ? = 0 to 0.25) pass the following conditions. (a) If the particle surface touches the surface of the cube, or if they are very close, it may not be possible to mesh or the generated finite element mesh will be distorted or some time meshing may not be possible at all. To avoid these, the center of the i th particle with radius was kept inside of the cube at some minimum distance l r l d 2 from the surface of the cube. (7.4) ll l rrd 1.0 2 += To fulfill the above condition, center coordinates of the i th particle with radius must pass the following check l r 3,2,1 3,2,1 2 2 =?? =? jdLx jdx li j li j (7.5) (b) If two adjacent particles overlap each other they will violate the rigid sphere condition. In addition, two adjacent particles can not touch each other. To fulfill this condition center coordinates of the i th particle which has radius must pass the following check l r ( ) )1(.....,,107.1 1 ?=+= ikrrd kl l (7.6) ( ) )1(.....,,1 1 ?=?? ikdxx lki vv (7.7) The modified RSA method has been used to achieve volume fraction higher than ? = 0.25. In this study, the unit cells with volume fraction, ? , such that 0 < ? < 0.25, all the particle centers were kept inside the cube, no particle overlaps the outer surfaces of the cube. For unit cell volume fractions more than 0.25, the particles were allowed to 135 overlap the surface boundary of the cube. If any particle overlapped the surface, the portion of the particle outside the cube was carefully cut into several sections and copied at a suitable position on the opposite surface of the cube. First, a cubic cell with volume fraction less than 0.25 was generated by fulfilling the above two conditions and then the cell was compressed in to a smaller size while keeping the size of the spherical particle same. The size of the sphere and total number of the sphere does not change during the operation, but the cell is compressed, giving a higher volume fraction of filler content. Length of the cubic cell is compressed first, by multiplying it by a user defined shrinking factor c f (< 1). The new length of the cube is given by LcL fn ?= (7.8) The old position of the center of i th particle which has radius is given by l r i x r then new position of the particle will be given by i n x r f ii n cxx ?= rr (7.9) Once the particles have moved, they will be allowed to overlap the surface of the cube but no particle center will be outside of the cube. If the new coordinates of the i th particle falls outside of the cube then it moved back in the cube at a random position between the surface and at an inward distance i n x r ( )?? l r , where ? is a user defined constant. 136 If the particle center sits at a distance from the surface of the cube then the particle surface will touch the cube surface and meshing will not be possible. If the new coordinates of the i l r th particle i n x r falls inside of the cube and at a distance from the surface then it will be moved at a random position between the surface and at an inward distance( l r )?? l r , where ? is a user defined constant. If the particle with radius overlaps any other previously accepted particle then it will be moved to a new random position in a random direction. The new position will be given by l r ??= i n i n xx rr (7.10) where ? is a small random number. Smaller ? ensures faster convergence. The new position of the particle will be accepted only if all above conditions are fulfilled. Above algorithm involves generation of random size of sphere, random movement and random positioning, hence each of the iterations may not produce an acceptable distribution. The algorithm also counts the number of iterations and stops the execution if the following condition is not full filled, here i is number of iteration and ? is a user defined constant. Generally ? is a very big number. ??i (7.1) 7.2 Finite Element Model of Unit Cell Finite element models of the cubic unit cells have been created to predict the bulk modulus, Poisson?s ratio and and viscoelastic properties of the nano underfill. In these 137 138 models, spherical random sized silica fillers were generated in an epoxy cube. The length of the epoxy cube is equal to the length of the unit cell considered in the unit cell algorithm. The filler particles were generated by using the center coordinates and radius. The model geometry was meshed by 8 node tetrahedral brick element. Figure 7.1 to Figure 7.3 shows representative models of the unit cell and filler distribution for ? = 0.20. Figure 7.4 shows the isotropic distribution of particles with volume fraction, ? = 0.39. Figure 7.5 shows example of the mesh generated in the FEA models. 7.3 Prediction of Coefficient of Thermal Expansion (CTE) The developed finite element model with random diameter filler particles was used to predict the coefficient of thermal expansion of the nano-underfill for different volume fraction of filler particles. Symmetric boundary conditions were used at the cube faces at x, y, and z equal to zero. The degrees of freedom were coupled at the faces at x, y, z equal to L. This ensures that the cube faces will not be distorted after deformation. The temperature of the model was raised to a user defined uniform temperature. Calculated CTE values are presented in Figure 7.6. In Figure 7.6, CTE in the direction of x-axis has been plotted for different volume fraction of the filler particles. Figure 7.6 shows that CTE of the underfill decreases linearly with the increase of volume fraction of the filler particles, and that the same results as Figure 6.8 were obtained. 7.4 Prediction of Elastic Modulus The elastic modulus of nano underfill was also predicted for different volume fractions by performing finite element analysis of the models with random fillers. For Figure 7.1: Isometric View of the Size and Position Distribution of the Filler Particles, ? = 0.20 139 Figure 7.2: Front View of the Size and Position Distribution of the Filler Particles, ? = 0.20 140 Figure 7.3: Right Side View of the Size and Position Distribution of the Filler Particles, ? = 0.20 141 Figure 7.4: Right Side View of the Size and Position Distribution of the Filler Particles, ? = 0.39 142 Figure 7.5: FEA Mesh of Unit Cell (Isometric View) 143 0 10 20 30 40 50 60 70 0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5 Volume Fraction of Filler Particle CT E ( p p m / ? C) Figure 7.6: Prediction of CTE by finite-element analysis of Unit Cell 144 this analysis, different boundary conditions than the CTE models have been used. Displacement in the x-direction is zero at the cube surface at x equal to zero. Fixed boundary condition is applied at couple of nodes at the mid position of the surface at x equal to zero. Tensile load was applied at the surface at x equal to L in the form of uniform pressure. A suitable small value of the pressure load will give better results because the modulus is defined as the initial slope of the linear part of stress-strain curves. Figure 7.7 shows the elastic modulus vs. volume fraction of filler particle. Elastic modulus increases exponentially with the increase of volume fraction of the fillers particles, and very similar results to Figure 6.10 were obtained. 7.5 Prediction of Bulk Modulus The bulk modulus of the nano-underfill was calculated by finite element analysis of the above unit cell. A fixed boundary condition is applied on the node at the origin. The origin is located at one corner of the unit cell. The x, y and z deflections are kept zero on the surfaces located at x = 0, y = 0 and z = 0 respectively. A uniform tensile loading in the form of hydrostatic pressure is applied on the surfaces located at x = L, y = L and z = L respectively. The degrees of freedom were coupled in the corresponding directions at x = L, y = L and z = L respectively. The loading causes an increase in the volume of the cube without changing the shape at all. As the unit cells were isotropic it generated a uniform change in length in all directions. Bulk modulus values were calculated by using the following formula V V k ? = ? (7.12) 145 0 2 4 6 8 10 12 0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5 Volume Fraction of Filler Particle E l ast i c M odul us , E ( G P a ) Figure 7.7: Prediction of Elastic Modulus by FEA Analysis of Unit Cell 146 147 Where, k = bulk modulus, ? = value of applied hydrostatic stress, V = original volume, ?V = change in volume. The analysis was performed for unit cells with ? = 0, 0.1, 0.2, and 0.25. Figure 7.8 shows the predicted bulk modulus as a function of volume fraction of the filler particle. Here the bulk modulus of the epoxy is 4.16 GPa and it increases exponentially with increases of volume fraction of the filler particles. 7.6 Prediction of Poisson?s Ratio Finite element analysis of the unit cell was also performed to predict the Poisson?s ratio of the nano-underfill. As like the bulk modulus, new unit cells in which size distribution of the filler particles were applied have chosen for this purpose. Displacement in the x-direction is zero at the cube surface at x equal to zero. Fixed boundary condition is applied at one of the point at the mid position of the surface at x equal to zero. Tensile load was applied on the surface at x equal to L in the form of uniform pressure. Nodal degrees of freedom on surfaces at y = 0 and L, z = 0 and L are coupled at corresponding directions. Loading causes extension in x direction while contraction in y and z direction. Poison?s ratio was calculated by using the extension and contraction values. Figure 7.9 shows that Poisson?s ratio decreases nonlinearly with increase of filler volume fraction. 7.7 Theoretical Background of Viscoelastic Model Prediction of time dependent properties for the underfill materials is very important. Commercial underfills show creep behavior, relaxation behavior and strain rate dependence. Silica particles used in underfills as fillers show linear elastic behavior. 4.00 4.80 5.60 6.40 7.20 8.00 0 0.1 0.2 0.3 0.4 0.5 Volume Fraction of Filler Particles B u l k M o dul us , k ( G P a ) Figure 7.8: Prediction of Bulk Modulus by FEA Analysis of Unit Cell 148 0.25 0.3 0.35 0.4 0.45 0 0.1 0.2 0.3 0.4 0.5 Volume Fraction of Filler Particles Po i s s o n ' s R a t i o Figure 7.9: Prediction of Poisson?s Ratio by FEA Analysis of Unit Cell 149 150 Viscoelastic property of the epoxy matrix is responsible for the time dependent behavior of underfill. A material will be viscoelastic if it has both elastic and viscous behavior. These materials exhibit instantaneous elasticity during loading which follows a slow and continuous increase of strain at a decreasing rate. During unloading a continuously decreasing strain is observed which follows an instantaneous elastic recovery. Viscoelastic materials also have the stress and strain rate dependence i.e. the longer the time to reach the final value of stress at a constant rate stressing, the larger is the corresponding strain. On the other hand the longer the time to reach the final value of strain at a constant rate straining, the smaller is the corresponding stress. For the time dependent analysis, a generalized Maxwell material model is assumed to represent the material behavior of the viscoelastic epoxy matrix in underfill. The Maxwell model consists by connecting a linear spring element to a linear viscous dashpot element in series [15]. Figure 7.10 and 7.11 shows the configuration and stress relaxation behavior under constant strain situation of a Maxwell element. The behavior of generalized Maxwell model formed by connecting several Maxwell elements in series becomes identical to the behavior of a single Maxwell element as described in Figure 7.12. On the other hand the generalized Maxwell model formed by connecting several Maxwell elements in parallel represents instantaneous elasticity, delayed elasticity with various retardation times, stress relaxation with various relaxation times and also viscous flow, Figure 7.13. The generalized Maxwell model formed by parallel elements is more convenient than generalized Kelvin model formed by parallel elements, which will be discussed later, to predict the stress associated with a prescribed strain variation since the same prescribed ? ? ? 2 ? 1 ? R ? Figure 7.10: Maxwell Model 151 ? t 0 ? (a) Applied Constant Strain t 0 ?? R= ? (b) Stress Response Figure 7.11: Stress Relaxation of Maxwell Element under Constant Strain 152 k 1 ? 1 1 N 1i i k 1 k ? = ? ? ? ? ? ? ? ? = ? 1 N 1i i 1 ? = ? ? ? ? ? ? ? ? = ? ? ? k 2 ? 2 k i ? i Figure 7.12: Generalized Maxwell Model by Connecting Elements in Series 153 ? 2 ? i? 1 k i k 2 k 1 Figure 7.13: Generalized Maxwell Model by Connecting Elements in Parallel 154 155 strain is applied to each individual element, and also the resulting stress is the sum of the individual contributions. Kelvin model is formed by connecting a spring element parallel to a dashpot element as shown in Figure 7.14. Generalized Kelvin model formed by connecting several Kelvin elements in parallel produces behavior as an equivalent Kelvin model, as shown in Figure 7.15. The generalized Kelvin model in series, Figure 7.16, is more convenient for viscoelastic analysis in cases where the stress history prescribed since same prescribed stress is applied to each individual element, and the resulting strain is the sum of the individual strain in each element. Again generalized Maxwell model in parallel is more convenient in cases where the strain history is prescribed since same prescribed strain is applied to each individual element, and also resulting stress is the sum of the individual contribution. 7.8 ANSYS Input Constants for Viscoelastic Material The finite element analysis software ANSYS was used to compute the predicted time dependent behavior of the developed unit cell model. In ANSYS a viscoelastic material can be defined by two methods [16]. The first method requires defining total of 95 constants. The list of the constants and their explanations are available in ANSYS documentation. Although the available automatic curve fitting procedure in ANSYS can produce some of the above 95 constants but defining total 95 constant is a very complicated method. The other available method requires defining fewer numbers of constants. In this method the required constants are Williams-Landel-Ferry (WLF) shift function constants and prony series constants of volumetric and shear response. ? ? ? 2 ? 1 ? k Figure 7.14: Kelvin Model 156 ? i ? 2 ? 1 k i k 2 k 1 ? = = N 1i i k k ? = = N 1i i ?? Figure 7.15: Generalized Kelvin Model by Connecting the Elements in Parallel 157 k 1 ? 1 k 2 ?2 k i ? i Figure 7.16: Generalized Kelvin Model by Connecting the Elements in Series 158 WLF shift functions are calculated from the elastic modulus relaxation data and time- temperature superposition method is used in this purpose [15]. Figure 7.17 shows the stress relaxation data that will be used in calculation of the viscoelastic constants. Figure 7.18 shows the log-log plot of the temperature dependent relaxation modulus versus time data obtained from the data shown in Figure 7.17. The curve at 25?C was taken as the reference curve and other curves at 75?C, 100?C and 125?C were shifted sideways parallel to the time axis to an appropriate distances and the single master stress relaxation curve was formed, Figure 7.19. After sifting by appropriate amount log stress relaxation modulus curves shown in Figure 7.18, needed to be extrapolated to get the single continuous master curve without any break, Figure 7.19. The magnitude of the total shift is plotted against temperature and given in Figure 7.20. Williams-Landel-Ferry (WLF) constants were calculated by substituting the temperature dependent shift function a T and temperature value into equation 7.13 and then performing nonlinear regression of that equation. Equation 7.13 has been originally proposed by Williams, Landel and Ferry and it gives the relation between the shift factor and the temperature. Calculated Williams- Landel-Ferry constants are C 1 = T 0 = 25 ?C, C 2 = -42.6, C 4 = 517 ?C () ( ) o o T TTC TTC aLog ?+ ?? = 4 2 (7.13) ANSYS constants to represent volumetric response and shear response were calculated by using Prony series. Prony series is derived from the solution of generalized Maxwell model in parallel. Stress-strain relation of spring and dashpot in Figure 7.10 can be given as 159 Time (min) 0 204060801012 S t re s s (M P a ) 0 5 10 15 20 25 ? C 100 ? C 125 ? C 75 ? C Strain Level = 1% 0 Figure 7.17: Stress Relaxation Data Used for Calculating ANSYS Viscoelastic Constant 160 Log Time, sec 0.00.51.01.52.02. Lo g St re s s R e la xa tion Mo du lu s E( t) , G P a -0.8 -0.6 -0.4 -0.2 0.0 0.2 0.4 25 ? C 75 ? C 125 ? C 100 ? C 5 Figure 7.18: Log-Log Plot of the Relaxation Modulus versus Time Data 161 Log Reduced Time Scale for T = 25 ? C, sec 02468 L og St r e ss Rel a xat i on M odul us E ( t ) , G P a -0.8 -0.6 -0.4 -0.2 0.0 0.2 0.4 Master Stress Relaxation Modulus at T = 25 ? C 10 Figure 7.19: Master Relaxation Modulus at 25 ?C 162 Temperature, T ( ? C) 0 25 50 75 100 125 150 S h i f t Fa ct or , Lo g a T -10 -8 -6 -4 -2 0 2 Figure 7.20: Temperature Dependent Shift Factor 163 2 ?? k= (7.14) 1 ? = ??? (7.15) Total strain, 21 ??? += (7.16) and strain rate ??? += 2 1 ??? (7.17) Inserting Equation 7.15 and time derivative of Equation 7.14 in Equation 7.17 will give ? ?? ? += ? ? k (7.18) or, ? ??? += dt d kdt d 1 (7.19) Taking Laplace transform of Equation 7.19 () () () () () ? ? ???? s ss k ss ? ?? + ? ? ? ? ? ? ?=? 0 1 0 Since () () 000 ==?? , () () () ? ? ?? s ss k ss ? ?? + ? ? ? ? ? ? = 1 Rearranging gives, () ()s k s sk s ?? + = ? ? ? (7.20) 164 Assuming applied strain is a step function and ( ) ( )tHt 0 ?? = Then, () s s 1 0 ?? = ? Substituting in Equation 7.20 gives, () s k s sk s 0 ? ? ? + = ? or, () ? ? ? k s k s + = ? 0 (7.21) Taking inverse Laplace transform of Equation 7.21 gives, () ? ? ? ? ? ? ? ? ?= t k kt ? ?? exp 0 Assuming k ? ? = gives, () ? ? ? ? ? ? ?= ? ?? t kt exp 0 (7.2) For generalized Maxwell model several Maxwell elements are connected in parallel as shown in Figure 7.13. Since a set of parallel Maxwell elements are in parallel connection, the stress relaxation response will be obtained by summing the stress relaxation response of the individual elements. For the combination 165 () ? = ? ? ? ? ? ? ? ? ?= N i i i t kt 1 0 exp ? ?? (7.23) Prony series representation of the elastic modulus can be obtained by dividing the stress response equation by applied constant strain. Derived prony series will have the following form () ? = ? ? ? ? ? ? ? ? ?= N i i i t ktE 1 exp ? (7.24) Now, a generalized parallel Maxwell model with a free spring in parallel of modulus , Figure 7.21, will give the modulus relaxation response to a constant strain ? k 0 ? as shown in Equation 7.25, () ? = ? ? ? ? ? ? ? ? ? ?+= N i i i t kktE 1 exp ? (7.25) Denoting k?s by E?s will give () ? = ? ? ? ? ? ? ? ? ? ?+= N i i i t EEtE 1 exp ? (7.26) ANSYS allows maximum ten Maxwell elements to approximate the relaxation function. Equation 7.27 and 7.28 are used to approximate the relaxation of shear modulus and bulk modulus. 166 k ? ? 2 ? i? 1 k i k 2 k 1 Figure 7.21: Generalized Parallel Maxwell Model with a Free Spring in Parallel 167 () ? = ? ? ? ? ? ? ? ? ? ?+= G n i G i i t GGtG 1 exp ? (7.27) () ? = ? ? ? ? ? ? ? ? ? ?+= k n i k i i t kktk 1 exp ? (7.28) Where, ComponentProny Each for Time Relaxation and Relaxationk eApproximat toElements Maxwell ofNumber n RelaxationG eApproximat toElements Maxwell ofNumber n (GPa) ModulusBulk Final k (GPa) ModulusShear Final G k i G i k G = = = = = ? ? ?? For the present study, total 5 Maxwell elements are considered. Shear modulus and bulk modulus relaxation data were calculated from Elastic modulus relaxation data by using Equations 7.29 and 7.30. Temperature dependent prony constants were calculated by the nonlinear regression of equation 7.27 and 7.28. Calculated prony series constants are given in Table 7.1 and 7.2. () ( ) ()?+ = 12 tE tG (7.29) () ( ) ()?213 ? = tE tk (7.30) 168 Table 7.1: Prony Viscoelastic Shear Response 25 ?C 75 ?C 100 ?C 125 ?C G 1 0.0059 0.0208 0.1269 0.0215 G 1 ? 9.39 17.55 0.9385 2.92 G 2 0.0196 0.0278 0.0256 0.1146 G 2 ? 10.91 35.43 13.91 0.3227 G 3 0.0646 0.0887 0.0374 0.0396 G 3 ? 1.09 1.60 12.85 23.02 G 4 0.0195 0.0294 0.0376 0.0212 G 4 ? 12.12 32.07 12.85 2.94 G 5 0.018 0.0298 0.0093 0.0188 G 5 ? 32.23 34.81 66.02 2.73 169 Table 7.2: Prony Viscoelastic Volumetric Response 25 ?C 75 ?C 100 ?C 125 ?C k 1 0.0972 0.0922 0.1148 0.07078 k 1 ? 16.09 31.15 13.24 2.84 k 2 0.0543 0.0804 0.4311 0.1345 k 2 ? 16.12 31.15 0.9412 23.06 k 3 0.2078 0.0922 0.1146 0.0709 k 3 ? 1.24 31.15 13.24 2.84 k 4 0.021 0.0919 0.1132 0.0679 k 4 ? 1.27 31.15 13.21 2.93 k 5 0.0429 0.3049 0.0288 0.3894 k 5 ? 16.03 1.71 58.52 0.3223 170 171 ANSYS provides some automatic curve fitting method to calculate the WLF, shear and bulk relaxation constant. In this method only the raw relaxation data need to be entered in the program and all the constants are calculated automatically by curve fitting procedure. This is a very easy and convenient way of calculating some of the input viscoelastic constants, however during calculating the WLF constants this method sets the reference temperature automatically by its own. Since the selection of reference temperature is not in user?s control, the defined reference temperature may not be useful for the user. For this reason this automatic curve fitting procedure was not used to calculate the required viscoelastic constants for ANSYS. 7.9 Prediction of Relaxation Behavior FEA modeling of the unit cell has developed to predict the relaxation behavior as a function of time. Boundary conditions were applied at the cube surface at y equal to zero. Displacement in the y-direction is zero at the cube surface at y equal to zero. Fixed boundary condition is applied at couple of nodes at the mid position of the surface at y equal to zero. Tensile loading was applied at the surface at y equal to L in the form of uniform nodal displacement. Step displacement loading was chosen instead of ramp loading to make sure instantaneous application of strain as shown in Figure 7.11(a). The value of applied instantaneous strain was 1% and it was kept constant throughout the solution. Filler particles were assumed to be linear elastic and only linear elastic material properties were applied for the filler material. Viscoelastic material behavior of the parallel Maxwell model shown in Figure 7.21 was applied for the epoxy matrix. In order to define the viscoelastic properties of epoxy WLF constant data, Prony viscoelastic shear 172 response data, Prony viscoelastic volumetric response data of epoxy were applied. In ANSYS for any viscoelastic material, linear elastic material properties also need to be applied along with the viscoelastic properties. The elastic modulus value of the epoxy is 2.5 GPa, which is the initial slope of the stress-strain curve of a tensile test. The above elastic modulus value will not be appropriate for viscoelastic analysis, since the applied strain loading is instantaneous. During actual relax testing the stress-strain value starts from zero and ramps up to the applied strain level and then relaxing starts. Since the epoxy matrix is not a linear elastic material, after starting the tests the instantaneous modulus or the ratio of instantaneous stress to instantaneous strain continuously decreases and at the point where relaxing starts this value becomes much lower than the initial elastic modulus. In this present FEA analysis this initial ramping is ignored and instantaneous strain and stress is assumed to be applied at time equal to zero. During calculating the Prony constants the initial ramp in the stress-strain data was also ignored and only the usefull relax portion of the data was utilized. So using of this initial elastic modulus value, 2.5 GPa, will considerably over predict the relaxation behavior of the underfill. In this case the actual instantaneous modulus at the time of starting of the relaxation should be used. For example at the point where relaxation starts the ratio of the instantaneous stress to the instantaneous strain for the 25 ?C curve in Figure 7.17 is 1.84 GPa. Hence the present FEA modeling is performed to predict the relaxation behavior at 25 ?C this 1.84 GPa was applied as the elastic modulus value of the viscoelastic epoxy matrix. Nonlinear solutions were performed in this case with several time steps. The nodal stress response of the surface at y = L were computed and plotted in Figure 7.22. In Figure 7.22 the prediction for the relaxation behavior of the underfill 0 0.5 1 1.5 2 2.5 3 3.5 2040608010 Time, t (min) St r e s s R e l a x a t i o n M o d u l u s , E ( G Pa ) Experimental-20% Filler Volume Fraction FEA Prediction-20% Filler Volume Fraction Experimental-10% Filler Volume Fraction FEA Prediction-10% Filler Volume Fraction Epoxy Only 120 Figure 7.22: Prediction of Elastic Modulus Relaxation 173 174 which has 20% nano-silica fillers is presented. In Figure 7.22, the relaxation of instantaneous modulus is plotted. Figure 7.22 shows that the predicted value for 20% volume fraction correlates very well with the experimentally measured relaxation behavior. For the 10% filler volume fraction a small deviation is observed between the predicted and experimental values. The reason behind this is the experimental curve at 10% filler volume fraction was generated with some unqualified specimen since the supply of the nano underfill was very limited. Further testing needed to be done at 10% volume fraction level in order to validate the predicted data properly at this volume fraction level. In the same figure the relaxation behavior of the epoxy only data is also plotted. Figure 7.23 and 7.24 shows the prediction of the shear modulus and bulk modulus relaxation. Shear modulus and bulk modulus were calculated from the elastic modulus by using equations 7.29 and 7.30. Like the instantaneous modulus relaxation same types of behavior were observed for the shear and bulk modulus relaxation. 0 0.2 0.4 0.6 0.8 1 1.2 1.4 0 20 40 60 80 100 120 Time, t (min) S h ear R e l a xa t i on M odul us, G ( G P a Experimental-20% Filler Volume Fraction FEA Prediction-20% Filler Volume Fraction Experimental-10% Filler Volume Fraction FEA-10% Filler Volume Fraction Epoxy Only Figure 7.23: Prediction of Shear Modulus Relaxation 175 0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 2040608010 Time, t (min) B u l k R e l a xa t i o n M o d u l u s, k ( G P a Experimental-20% Filler Volume Fraction FEA Prediction-20% Filler Volume Fraction Epoxy Only Experimental-10% Filler Volume Fraction FEA-10% Filler Volume Fraction 120 Figure 7.24: Prediction of Bulk Modulus Relaxation 176 177 CHAPTER 8 CONCLUSIONS In this present study, a methodology has been developed to create statistically isotropic cubic unit cells of underfill materials containing randomly distributed, random sized micro and nano spherical fillers. Finite element models based on random sequential addition algorithm have been developed and applied successfully for prediction of linear and non-linear material properties of underfill materials. Model predictions of properties have been validated with experimental data. The multi-filler models with randomly oriented, random sized spherical filler particles predict the properties more accurately than one or eight filler models with periodic filler distribution. Multi-mesh extrapolation has been used to model macro behavior in uniaxial test. Finite element results indicate that the length-to-width aspect ratio of thin underfill specimen of 20 or more is required for accurate uniaxial testing of nano and micro- underfills. Mesh refinement has been shown to reduce uncertainty in predicted stress- strain behavior in the neighborhood of 0.2%. A cold chamber has been developed newly which can be used with mechanical testing machine to perform uniaxial tests at very low temperature up to -180?C. A newly developed specimen preparation technique has been implemented to generate very thin (3 to 5 mil) specimen of nano and micro underfill materials. Temperature dependent stress- strain curves showed that both nano and micro underfill become viscoplastic at very high temperatures and linear elastic at very low temperatures. REFERENCES 1 Liji, Z., Li, W., Xiaoming, X., Kempe, W., ?An Investigation on Thermal Reliability of Underfilled PBGA Solder Joint,? IEEE Transactions on Electronics Packaging Manufacturing, vol. 25, no. 4, October 2002 2 Sillanpaa, M., Okura, J.H., ?Flip Chip on Board: Assessment of Reliability in Cellular Phone Application?, IEEE-CPMT Vol.27, Issue:3, pp. 461 ? 467, Sept.2004 3 Pascariu G., Cronin P, Crowley D, ?Next-generation Electronics Packaging Using Flip Chip Technology?, Advanced Packaging, Nov.2003 4 Jung, E., Heinricht, K. Kloeser, J. Aschenbrenner, R. Reichl, H., Alternative Solders for Flip Chip Applications in the Automotive Environment, IEMT- Europe, Berlin, Germany, pp.82-91, 1998 5 Bedinger, John M., ?Microwave Flip Chip and BGA Technology?, IEEE MTT-S International Microwave Symposium Digest, v 2, pp 713-716, 2000 6 Van den Crommenacker, J., ?The System-in-Package Approach?, IEEE Communications Engineer, Vol 1, No. 3, pp. 24-25, June/July, 2003 7 Ray, S.K., Quinones, H., Iruvanti, S., Atwood, E., Walls, L., ?Ceramic Column Grid Array (CCGA) Module for a High Performance Workstation Application?, Proceedings - Electronic Components and Technology Conference, pp 319-324, 1997 8 Lau, J. H., and Chang, C., ?Characterization of Underfill Materials for Functional Solder Bumped Flip Chips on Board Applications,? IEEE Transactions on Components and Packaging Technologies, Vol. 22(1), pp. 111- 119, 1999 9 Liu, J.; Kraszewshi, R.; Lin, X.; Wong, L.; Goh, S. H.; Allen, J., ?New Developments in Single Pass Reflow Encapsulant for Flip Chip Application,? Proc International Symposium on Advanced Packaging Materials, Atlanta, GA, pp. 74-79, March 2001 178 10 Shi, S. H.; Wong, C. P., ?Recent Advances in the Development of No-Flow Underfill Encapsulants ? A Practical Approach towards the Actual Manufacturing Application?, IEEE Transactions on Electronics Packaging Manufacturing, Vol. 22, pp. 331-339, 1999 11 Drugan, W. J., ?Micromechanics-Based Variational Estimates for a Higher- Order Nonlocal Constitutive Equation and Optimal Choice of Effective Moduli of Elastic Composites,? J. Mech. Phys. Solids, 48, pp 1359-1387, 2000 12 Drugan, W. J., Wills, J. R., ?A Micromechanics-Based Nonlocal Constitutive Equation and Estimates of the Representative Volume Element Size for Elastic Composites,? J. Mech. Phys. Solids, 44, pp 497-524, 1996 13 Segurado, J., Llorca, J., ?A Numerical Approximation to Elastic Properties of Sphere-Reinforced Composites,? Journal of Mechanics and Physics of Solids, Vol. 50, pp. 2107-2121, 2002 14 Cook, R. D., ?Concepts and Applications of Finite Element Analysis?, 4 th Edition, 2001 15 Findley, W. N., Lai, J. S., Onaran, K., ?Creep and Relaxation of Nonlinear Viscoelastic Material,? Dovar Publications, Inc., New York, 1976 16 ANSYS Documentation, ?Viscoelastic Material Constants,? Section 2.5.3, ANSYS Element Reference Manual, ANSYS Release 9.0 Documentation, 2005 17 Andersson, O. and Berkyto, E., ?Some Factors Affecting the Stress Strain Characteristics of Paper,? Svensk Papperstidning, Vol. 54(13), pp. 437-444, 1951 18 Bo, C., Li, W., Qun, Z., Xia, G., Xiaoming, X., Kempe, W., ?Reliability and New Failure Modes of Encapsulated Flip Chip on Board Under Thermal Shock Testing,? Electronic Packaging Technology Proceedings, ICEPT, pp 416 ? 421, 2003 19 Islam, M. S., Suhling J. C., Lall, P., ?Measurement of the temperature dependent constitutive behavoir of underfill encapsulants,? Intersociety Conference on Thermal and Thermomechanical Phenomena in Electronic Systems, ITHERM , Vol. 2 , pp 145 - 152, June 2004 20 Islam, S., Suhling, J. C., Lall, P., ?Measurement of the Temperature Dependent Constitutive Behavior of Underfill Encapsulants,? IEEE Transactions on Components and Packaging Technologies, Vol. 28, No. 3, pp. 467-476, September 2005 179 21 Lall, P., Islam, S., Suhling, J. C., Tian, G., ?Nano-Underfills for High-Reliability Applications in Extreme Environments,? The 55 th Electronic Components and Technology Conference (ECTC), 2005 22 Rahim, M. K., Suhling, J. C., Copeland, S., Islam, M. S., Jaeger, R. C., Lall, P., Johnson, R. W., ?Die Stress Chracterization in Flip Chip Laminate Assemblies,? IEEE Transactions on Components and Packaging Technologies, Vol. 28, No. 3, pp. 415-429, September 2005 23 Islam, M. S., Suhling, J. C., Lall, P., ?Measurement of the Constitutive Behavior of Underfill Encapsulants,? Proceedings of IPACK03, International Electronic Packaging Technical Conference and Exhibition, 2003 24 Abdel-Hady, H., Ma, H., Suhling, J. C., Islam, M. S., Lall, P., ?Measurement of the Constitutive Behavior of Lead Free Solders? SEM, 2004 25 Islam, M. S., Suhling J. C., Lall, P., Xu, B., Johnson, R. W, ?Measurements and Modeling of the Temperature Dependent Material Behavior of Underfill Encapsulants,? The 53 rd Electronic Components and Technology Conference (ECTC), 2003 26 Suhling J. C., Islam, M. S., Lall, P., Xu, B., Johnson, R. W, ?Accurate Measurements of the Material Behavior of Underfill Encapsulants,? IMAPS, 2004 27 Islam, S., Suhling J. C., Lall, P., ?Measurement of Underfill to Soldermask Adhesion Strength?, SEM, 2003 28 Rahim, K, Islam, S., Lin, C., Ma, H., Suhling J. C., Jaeger, R. C., Lall, P., ?Packaging Stress Characterization at Low Temperatures, IMAPS, 2005 29 TN 513-1, Measurement Group, Technical Note, 2004 30 Pyrz, R., ?Correlation of Microstructure Variability and Local Stress Field in Two-Phase Materials,? Mater. Sci. Eng., Mater. Sci. Eng., A 177, pp 253-259, 1994a 31 Pyrz, R., Quantitative Description of the Microstructure of Composites. Part I: Morphology of Unidirectional Composite Systems, Composites Science and Technology, Vol. 50, pp. 197 ? 208, 1994b 32 Van den Crommenacker, J., ?The System-in-Package Approach?, IEEE Communications Engineer, Vol 1, No. 3, pp. 24-25, June/July, 2003 180 33 Yeh, K. C., Considine, J. M., Suhling, J. C., ?The Influence of Moisture Content on the Nonlinear Constitutive Behavior of Cellulosic Materials,? Proceedings of the 1991 International Paper Physics Conference, TAPPI, pp. 695-711, Kona, Hawaii, September 22-26, 1991 34 Qian Z., Lu M., and Liu S., ? Constitutive Modeling of Polymer Films Based on A Molecular Chain Network Theory? Advances in Electronic Packaging, ASME EEP-vol. 19-2, 1997a, pp. 1291-1298 35 Merton M. M., Mahajan R. L., and Nikmanesh N., ?Alternative Curing Methods for FCOB Underfill? Advances in Electronic Packaging, ASME EEP, vol. 19-1, 1997, pp. 291-300 36 Qian, Z., Wang, J., Yang, J., Liu, S., ?Visco-Elastic-Plastic Properties and Constitutive Modeling of Underfills?, IEEE Transactions on Components and Packaging Technology?, vol. 22, No. 2, pp. 152-157, June 1999 37 Shi, S. H., Yao, Q. Qu, J., Wong, C. P., ?Study on The Correlation of Flip-Chip Reliability with Mechanical Properties of No-Flow Underfill Materials?, Proceedings of International Symposium on Advanced Packing Materials, pp. 271-277, March 2000 38 Schubert A., 1997, ?Thermo-Mechanical Reliability Analysis of Flip-Chip Assemblies by Combined Microdac and the Finite Element Method?, Advances in Electronic Packaging, ASME EEP-Vol. 19-2, 1647-1654 39 M. Lu, W. Ren and S. Liu, ?A Multi-Axial Thermo-Mechanical Fatigue Tester for Electronic Packaging Material?, EEP-Vol. 17, ?ASME Symposium of Censoring, monitoring, modeling and verification of Electronic Packaging, Nov. 17-21, 1996, Atlanta, pp 87-92 40 M. Lu, Z. Qian and S Liu, ?Unified Multi-Axial Submicron Fatigue Tester for Miniaturized Specimen?, EEP-Vol. 19-1, Advances in Electronic Packaging- 1997, ASME 1997, Volume 1, pp 1145-1150 41 Ren, W., ?Thermo-Mechanical Properties of Packaging Materials and Their Applications to Reliability Evaluation for Electronic Packages?, Ph. D. Thesis, Department of mechanical engineering, Wayne State University, Detroit, Michigan, 2000 42 Merton M. M., Mahajan R. L., and Nikmanesh N., 1997, ?Alternative Curing Methods for FCOB Underfill?, Advances in Electronic Packaging, ASME EEP- Vol. 19-1, 291-300 181 43 Widom, B., ?Random Sequential Addition of Hard Spheres to a Volume,? The Journal of Chemical Physics, Vol. 44, No. 10, pp. 3888-3894, 1966 44 Rintoul, M. D., Torquato, S., ?Reconstruction of the Structure of Dispersions,? Journal of Colloid and Interface Science, 186, pp. 467-476, 1997 45 Bohm, H.J., Han, W., ?Comparisons between three-dimensional and two- dimensional multi-particle unit cell models for particle reinforced MMCs.,? Model. Simul. Mater. Sci. Eng., 9, pp. 47?65, 2001 46 Gusev, A.A., ?Representative volume element size for elastic composites: a numerical study,? J. Mech., Phys. Solids, 45, 1449?1459, 1997 47 Gusev, A.A., Hine, P.J., Ward, I.M., ?Fiber packing and elastic properties of transversely random unidirectional glass/epoxy composite,? Comp. Sci. Technol. 60, 535?541, 2000 48 Pyrz, R., ?Correlation of Microstructure Variability and Local Stress Field in Two-Phase Materials,? Material Science and Engineering, A177, pp. 253-259, 1994 49 Pyrz, R., ?Quantitative Description of the Microstructure of Composites. Part I: Morphology of Unidirectional Composite Systems,? Composite Science and Technology,50, pp. 197-208, 1994 182 183 APPENDICES 184 APPENDIX A FILLER DISTRIBUTION IN DIFFERENT UNIT CELLS In the property prediction effort several unit cells have been developed to predict the properties of nano-silica underfill as a function of filler volume fractions. Detail view of the filler distribution in all unit cells are given in this section. Figure A.1: Random Distribution of Constant Size Fillers, ? = 0.10 [Isotropic View] 185 Figure A.2: Random Distribution of Constant Size Fillers, ? = 0.10 [Front View] 186 Figure A.3: Random Distribution of Constant Size Fillers, ? = 0.10 [RHS View] 187 Figure A.4: Random Distribution of Constant Size Fillers, ? = 0.10 [Top View] 188 Figure A.5: Random Distribution of Constant Size Fillers, ? = 0.20 [Isotropic View] 189 Figure A.6: Random Distribution of Constant Size Fillers, ? = 0.20 [Front View] 190 Figure A.7: Random Distribution of Constant Size Fillers, ? = 0.20 [RHS View] 191 Figure A.8: Random Distribution of Constant Size Fillers, ? = 0.20 [Top View] 192 Figure A.9: Random Distribution of Constant Size Fillers, ? = 0.25 [Isometric View] 193 Figure A.10: Random Distribution of Constant Size Fillers, ? = 0.25 [Front View] 194 Figure A.11: Random Distribution of Constant Size Fillers, ? = 0.25 [RHS View] 195 Figure A.12: Random Distribution of Constant Size Fillers, ? = 0.25 [Top View] 196 Figure A.13: Random Distribution of Constant Size Fillers, ? = 0.38 [Isometric View] 197 Figure A.14: Random Distribution of Constant Size Fillers, ? = 0.38 [Front View] 198 Figure A.15: Random Distribution of Constant Size Fillers, ? = 0.38 [RHS View] 199 Figure A.16: Random Distribution of Constant Size Fillers, ? = 0.38 [Top View] 200 Figure A.17: Random Distribution of Random Size Fillers, ? = 0.10 [Isometric View] 201 Figure A.18: Random Distribution of Random Size Fillers, ? = 0.10 [Front View] 202 Figure A.19: Random Distribution of Random Size Fillers, ? = 0.10 [RHS View] 203 Figure A.20: Random Distribution of Random Size Fillers, ? = 0.10 [Top View] 204 Figure A.21: Random Distribution of Random Size Fillers, ? = 0.20 [Isometric View] 205 Figure A.22: Random Distribution of Random Size Fillers, ? = 0.20 [Front View] 206 Figure A.23: Random Distribution of Random Size Fillers, ? = 0.20 [RHS View] 207 Figure A.24: Random Distribution of Random Size Fillers, ? = 0.20 [Top View] 208 Figure A.25: Random Distribution of Random Size Fillers, ? = 0.25 [Isometric View] 209 Figure A.26: Random Distribution of Random Size Fillers, ? = 0.25 [Front View] 210 Figure A.27: Random Distribution of Random Size Fillers, ? = 0.25 [RHS View] 211 Figure A.28: Random Distribution of Random Size Fillers, ? = 0.25 [Top View] 212 Figure A.29: Random Distribution of Random Size Fillers, ? = 0.39 [Isometric View] 213 Figure A.30: Random Distribution of Random Size Fillers, ? = 0.39 [Front View] 214 Figure A.31: Random Distribution of Random Size Fillers, ? = 0.39 [RHS View] 215 Figure A.32: Random Distribution of Random Size Fillers, ? = 0.39 [Top View] 216 217 APPENDIX B MATLAB PROGRAM FOR RANDOM COORDINATE GENERATION-I Following MATLAB program was used to generate random center coordinates of the constant size filler particles and for volume fraction up to 25 percent. %********************************************************************** %* Random Coordinate Generation * %* x1, x2, x3 are three coordinates of the center of the spherical particles * %********************************************************************** %*********************** Defining Number of Particles ********************* % Total Number of Particle = n * %********************************************************************** diary('random_coordinates_20%_volume_fraction.txt') vv=20; v=0; n=1; i=2; tvs=0; mmix=0; mmiy=0; mmiz=0; 218 x1vs=0; x2vs=0; x3vs=0; h=250; x1(1)=0.1; x2(1)=0.1; x3(1)=0.1; r=25; s1=2.07*r; s2=r+(0.1*r); fprintf('Random Coordinates of Filler Particles\n'); while v<=vv % Start of while #1 [n= no of sphere] a=rand(1,1); % Generating Real Number Between 0-1 b=rand(1,1); c=rand(1,1); d=randint(1,1,[0,h]); % Generating Integer Number Between 0-h e=randint(1,1,[0,h]); f=randint(1,1,[0,h]); p1=abs(a+d); % Generating real number p1, p2 and p3 p2=abs(b+e); p3=abs(c+f); 219 %******************* Checking for Acceptable Position ****************** if p1>=s2 % Checking that center is enough away from surface m1=1; % Start of if #1 else m1=0; continue % Go to the top of while loop with i=i end % End of if #1 if p2>=s2 % Start of if #2 m1=1; else m1=0; continue % Go to the top of while loop with i=i end % End of if#2 if p3>=s2 % Start of if #3 m1=1; else m1=0; continue % Go to the top of while loop with i=i end % End of if#3 if abs(p1-h)>=s2 % Start of if #4 m1=1; 220 else m1=0; continue % Go to the top of while loop with i=i end % End of if #4 if abs(p2-h)>=s2 % Start of if #5 m1=1; else m1=0; continue % Go to the top of while loop with i=i end % End of if #5 if abs(p3-h)>=s2 % Start of if #6 m1=1; else m1=0; continue % Go to the top of while loop with i=i end % End of if#6 j=1; while j=s1 % Start of if #7 m=1; j=j+1; else m=0; j=i; end % End of if #7 end % End of while #2 if m==1 % Start of if #8 x1(i)=p1; x2(i)=p2; x3(i)=p3; fprintf('wpave,%f,%f,%f\n',0,0,x3(i)); fprintf('sph4,%f,%f,%f\n',x1(i),x2(i),r); % Calculation of Volume vs=(4*pi*(r^3)/3); tvs=tvs+vs; % Calculation of Centroid x1vs=x1vs+(4*pi*(r^3)/3)*x1(i); x2vs=x2vs+(4*pi*(r^3)/3)*x2(i); x3vs=x3vs+(4*pi*(r^3)/3)*x3(i); 222 % Calculation of Inertia mmix=mmix+(1*(((2/5)*(r*r))+((x1(i)-h/2)*(x1(i)-h/2)))); mmiy=mmiy+(1*(((2/5)*(r*r))+((x2(i)-h/2)*(x2(i)-h/2)))); mmiz=mmiz+(1*(((2/5)*(r*r))+((x3(i)-h/2)*(x3(i)-h/2)))); i=i+1; else i=i; end % End of if #8 vc=h^3; v=100*(tvs/vc); end % end of while#1 x1bar=x1vs/tvs; x2bar=x2vs/tvs; x3bar=x3vs/tvs; fprintf('Length of One Side of the Cube, h = %d\n',h); fprintf('Volume of One Sphere = %f\n',vs); fprintf('Total Volume of Sphere = %f\n',tvs); fprintf('Volume of cube, vc = %f\n',vc); fprintf('Volume Fraction of Filler Particle = %f\n',v); fprintf('X Centroid = %f\n',x1bar); fprintf('Y Centroid = %f\n',x2bar); fprintf('Z Centroid = %f\n',x3bar); fprintf('Mass Moment of Inertia wrt x Centroidal Axis= %f\n',mmix); 223 fprintf('Mass Moment of Inertia wrt y Centroidal Axis= %f\n',mmiy); fprintf('Mass Moment of Inertia wrt z Centroidal Axis= %f\n',mmiz); fprintf('\nTotal %d particle is created \n',(i-2)); ii=2; while ii<=(i-1) fprintf('Sphere No,%d,x1,%f,x2,%f,x3,%f\n',ii,x1(ii),x2(ii),x3(ii)); ii=ii+1; end fprintf('\nThanks! for using this program\n') 224 APPENDIX C MATLAB PROGRAM FOR RANDOM COORDINATE GENERATION-II Following MATLAB program was used to generate random center coordinates of the constant size filler particles and for volume greater than to 25 percent. %********************************************************************** %* %Random Coordinate Generation * %***** x1, x2, x3 are three coordinates of the center of the spherical particles * %********************************************************************** %*********************** Defining Number of Particle ********************* % Total Number of Particle = n * %********************************************************************** diary('random_coordinate_volume_fraction_more_than_25_percent.txt') fprintf('SSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSS\n'); fprintf('SSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSS\n'); fprintf('Start of a New Run\n'); vv=20; v=0; n=1; i=2; 225 vs=0; tvs=0; tvsn=0; mmix=0; mmiy=0; mmiz=0; x1vs=0; x2vs=0; x3vs=0; mmixn=0; mmiyn=0; mmizn=0; x1vsn=0; x2vsn=0; x3vsn=0; h=250; r=25; s1=2.07*r; s2=r+(0.1*r); x1(1)=0.01; x2(1)=0.01; x3(1)=0.01; fprintf('Random Coordinates of Filler Particles\n'); 226 while v<=vv % Start of while #1 [n= no of sphere] a=rand(1,1); % Generating Real Number Between 0-1 b=rand(1,1); c=rand(1,1); d=randint(1,1,[0,h]); % Generating Integer Number Between 0-h e=randint(1,1,[0,h]); f=randint(1,1,[0,h]); p1=abs(a+d); % Generating real number p1, p2 and p3 p2=abs(b+e); p3=abs(c+f); %********************** Checking for Acceptable Position ***************** if p1>=s2 % Checking that center is enough away from surface m1=1; % Start of if #1 else m1=0; continue % Go to the top of while loop with i=i end % End of if #1 if p2>=s2 % Start of if #2 m1=1; else m1=0; continue % Go to the top of while loop with i=i end % End of if #2 227 if p3>=s2 % Start of if #3 m1=1; else m1=0; continue % Go to the top of while loop with i=i end % End of if #3 if abs(p1-h)>=s2 % Start of if #4 m1=1; else m1=0; continue % Go to the top of while loop with i=i end % End of if #4 if abs(p2-h)>=s2 % Start of if #5 m1=1; else m1=0; continue % Go to the top of while loop with i=i end % End of if #5 if abs(p3-h)>=s2 % Start of if #6 m1=1; else m1=0; continue % Go to the top of while loop with i=i 228 end % End of if #6 %*** Checking Distance From Previously Accepted Particles *** j=1; while j=s1 % Start of if #7 m=1; j=j+1; else m=0; j=i; end % End of if #7 end % End of while #2 %******************** Accepting Position *************************** if m==1 % Start of if #8 x1(i)=p1; x2(i)=p2; x3(i)=p3; fprintf('wpave,%f,%f,%f\n',0,0,x3(i)); fprintf('sph4,%f,%f,%f\n',x1(i),x2(i),r); vs=(4*pi*(r^3)/3); tvs=tvs+vs; 229 vc=h^3; v=100*(tvs/vc); x1vs=x1vs+(4*pi*(r^3)/3)*x1(i); x2vs=x2vs+(4*pi*(r^3)/3)*x2(i); x3vs=x3vs+(4*pi*(r^3)/3)*x3(i); mmix=mmix+(1*(((2/5)*(r*r))+((x1(i)-h/2)*(x1(i)-h/2)))); mmiy=mmiy+(1*(((2/5)*(r*r))+((x2(i)-h/2)*(x2(i)-h/2)))); mmiz=mmiz+(1*(((2/5)*(r*r))+((x3(i)-h/2)*(x3(i)-h/2)))); i=i+1; else i=i; end % End of if #8 end % End of while #1 x1bar=x1vs/tvs; x2bar=x2vs/tvs; x3bar=x3vs/tvs; fprintf('\n') fprintf('Length of One Side of the Cube, h = %f\n',h); fprintf('Volume of One Sphere, vs = %f\n',vs); fprintf('Total Volume of Sphere = %f\n',tvs); fprintf('Volume of cube, vc = %f\n',vc); fprintf('Volume Fraction of Filler Particle = %f\n',v); fprintf('X Centroid = %f\n',x1bar); 230 fprintf('Y Centroid = %f\n',x2bar); fprintf('Z Centroid = %f\n',x3bar); fprintf('Mass Moment of Inertia wrt x Centroidal Axis= %f\n',mmix); fprintf('Mass Moment of Inertia wrt y Centroidal Axis= %f\n',mmiy); fprintf('Mass Moment of Inertia wrt z Centroidal Axis= %f\n',mmiz); fprintf('\nTotal %d particle is created \n',(i-2)); fprintf('\n') zz=2; while zz<=(i-1) fprintf('Sphere No,%d,x1,%f,x2,%f,x3,%f\n',zz,x1(zz),x2(zz),x3(zz)); zz=zz+1; end %******************** Shirinking the Cube ****************** fprintf('*********************************************\n'); fprintf('*********************************************\n'); fprintf('Random Coordinates of Fillers After Shrinking\n'); sc=0.7937; hn=h*sc; ii=2; % Useful coordinates are from i=2 to i=i-1 x1new(1)=0.01; x2new(1)=0.01; x3new(1)=0.01; while ii<=(i-1) % Start of while # 3 231 stopcount=1; x1n=x1(ii)*sc; x2n=x2(ii)*sc; x3n=x3(ii)*sc; %******************* Checking for Acceptable Position ****************** %%*** and Checking Distance from Previously Accepted Particles ************* jj=1; while jjhn % Start of if #41 x1n=hn-randint(1,1,[0,(r-5)]); else 232 end % End of if#41 if x2n>hn % Start of if #51 x2n=hn-randint(1,1,[0,(r-5)]); else end % End of if #51 if x3n>hn % Start of if #61 x3n=hn-randint(1,1,[0,(r-5)]); else end % End of if #61 dist_iijj=sqrt(((x1n-x1new(jj))*(x1n-x1new(jj)))+((x2n-x2new(jj))*(x2n- x2new(jj)))+((x3n-x3new(jj))*(x3n-x3new(jj)))); if abs(dist_iijj)>=s1 % Start of if #71 mm=1; jj=jj+1; else mm=0; jj=1; u11=randint(1,1,[1,10]); u22=randint(1,1,[1,10]); u33=randint(1,1,[1,10]); uval=sqrt((u11*u11)+(u22*u22)+(u33*u33)); u1=(u11/uval); u2=(u22/uval); 233 u3=(u33/uval); x1n=x1n+(randint(1,1,[-10,10]))*u1; % Random moving of sphere x2n=x2n+(randint(1,1,[-10,10]))*u2; x3n=x3n+(randint(1,1,[-10,10]))*u3; end % End of if #71 stopcount=stopcount+1; if stopcount >= 500000 % Stopcount if ii=i; jj=ii; fprintf('Stopcount = %d\n',stopcount); else end % End of stopcount if end % End of while #4 %********************************************************************** if mm==1 % Start of if #81 x1new(ii)=x1n; x2new(ii)=x2n; x3new(ii)=x3n; fprintf('wpave,%f,%f,%f\n',0,0,x3new(ii)); fprintf('sph4,%f,%f,%f\n',x1new(ii),x2new(ii),r); vs=(4*pi*(r^3)/3); tvsn=tvsn+vs; x1vsn=x1vsn+(4*pi*(r^3)/3)*x1new(ii); 234 x2vsn=x2vsn+(4*pi*(r^3)/3)*x2new(ii); x3vsn=x3vsn+(4*pi*(r^3)/3)*x3new(ii); mmixn=mmixn+(1*(((2/5)*(r*r))+((x1new(ii)-hn/2)*(x1new(ii)-hn/2)))); mmiyn=mmiyn+(1*(((2/5)*(r*r))+((x2new(ii)-hn/2)*(x2new(ii)-hn/2)))); mmizn=mmizn+(1*(((2/5)*(r*r))+((x3new(ii)-hn/2)*(x3new(ii)-hn/2)))); ii=ii+1; else ii=ii; end % End of if #81 end % End of while# 3 vcn=hn^3; vn=100*(tvsn/vcn); fprintf('\n') fprintf('\n') fprintf('New Length of One Side of the Cube, hn = %f\n',hn); fprintf('Volume of Sphere, vs = %f\n',vs); fprintf('Total New Volume of Sphere, tvsn = %f\n',tvsn); fprintf('New Volume of cube, vcn = %f\n',vcn); fprintf('New Volume Fraction of Filler Particle = %f\n',vn); fprintf('\n') fprintf('\nTotal %d particle is created \n',(ii-2)); x1barn=x1vsn/tvsn; x2barn=x2vsn/tvsn; 235 x3barn=x3vsn/tvsn; fprintf('New X Centroid = %f\n',x1barn); fprintf('New Y Centroid = %f\n',x2barn); fprintf('New Z Centroid = %f\n',x3barn); fprintf('New Mass Moment of Inertia wrt x Centroidal Axis= %f\n',mmixn); fprintf('New Mass Moment of Inertia wrt y Centroidal Axis= %f\n',mmiyn); fprintf('New Mass Moment of Inertia wrt z Centroidal Axis= %f\n',mmizn); fprintf('\n') zzz=2; while zzz<=(i-1) fprintf('Sphere No,%d,x1new,%f,x2new,%f,x3new,%f\n', zzz,x1new(zzz), x2new(zzz), x3new(zzz)); zzz=zzz+1; end fprintf('\n') fprintf('\nThanks! for using this program\n') 236 APPENDIX D MATLAB PROGRAM FOR RANDOM COORDINATE GENERATION-III Following MATLAB program was used to generate random center coordinates of the random size filler particles and for volume fraction up to 25 percent. %********************************************************************** %* %Random Coordinate Generation * %* x1, x2, x3 are three coordinates of the center of the spherical particles * % Total Number of Particle = n * %********************************************************************** diary('random_coordinate_random_radius_20_percent_volume_fraction.txt') vv=20; v=0; n=1; i=2; tvs=0; mmix=0; mmiy=0; mmiz=0; x1vs=0; x2vs=0; x3vs=0; 237 h=250 ; x1(1)=0.1; x2(1)=0.1; x3(1)=0.1; radius(1)=0.1; fprintf('Random Coordinates of Filler Particles\n'); fprintf('\n') while v<=vv % Start of while #1 [n= no of sphere] %********************* Defining Radius ****************************** r1=(rand(1,1)); r2=(randint(1,1,[17,33])); r=(r1+r2); % Random radius s2=r+(0.1*r); %********************** Center Coordinate ************************* a=rand(1,1); % Generating Random Real Number Between 0-1 b=rand(1,1); c=rand(1,1); d=randint(1,1,[0,h]); % Generating Random Integer Number Between 0-h e=randint(1,1,[0,h]); f=randint(1,1,[0,h]); p1=abs(a+d); % Generating Real Number p2=abs(b+e); p3=abs(c+f); 238 %******************* Checking for Acceptable Position ****************** if p1>=s2 % Checking that center is enough away from surface m1=1; % Start of if #1 else m1=0; continue % Go to the top of while loop with i=i end % End of if#1 if p2>=s2 % Start of if #2 m1=1; else m1=0; continue % Go to the top of while loop with i=i end % End of if#2 if p3>=s2 % Start of if #3 m1=1; else m1=0; continue % Go to the top of while loop with i=i end % End of if #3 if abs(p1-h)>=s2 % Start of if #4 m1=1; else m1=0; 239 continue % Go to the top of while loop with i=i end % End of if #4 if abs(p2-h)>=s2 % Start of if #5 m1=1; else m1=0; continue % Go to the top of while loop with i=i end % End of if#5 if abs(p3-h)>=s2 % Start of if #6 m1=1; else m1=0; continue % Go to the top of while loop with i=i end % End of if #6 j=1; while j=s1 % Start of if #7 m=1; j=j+1; 240 else m=0; j=i; end % End of if #7 end % End of while #2 if m==1 % Start of if #8 x1(i)=p1; x2(i)=p2; x3(i)=p3; radius(i)=r; fprintf('wpave,%f,%f,%f\n',0,0,x3(i)); fprintf('sph4,%f,%f,%f\n',x1(i),x2(i),r); % Calculation of Volume vs=(4*pi*(r^3)/3); tvs=tvs+vs; % Calculation of Centroids x1vs=x1vs+(4*pi*(r^3)/3)*x1(i); x2vs=x2vs+(4*pi*(r^3)/3)*x2(i); x3vs=x3vs+(4*pi*(r^3)/3)*x3(i); % Calculation of Inertia mmix=mmix+(1*(((2/5)*(r*r))+((x1(i)-h/2)*(x1(i)-h/2)))); mmiy=mmiy+(1*(((2/5)*(r*r))+((x2(i)-h/2)*(x2(i)-h/2)))); mmiz=mmiz+(1*(((2/5)*(r*r))+((x3(i)-h/2)*(x3(i)-h/2)))); 241 i=i+1; else i=i; end % End of if #8 vc=h^3; v=100*(tvs/vc); end % End of while #1 x1bar=x1vs/tvs; x2bar=x2vs/tvs; x3bar=x3vs/tvs; fprintf('\n') fprintf('Length of One Side of the Cube, h = %d\n',h); fprintf('Volume of One Sphere = %f\n',vs); fprintf('Total Volume of Sphere = %f\n',tvs); fprintf('Volume of cube, vc = %f\n',vc); fprintf('Volume Fraction of Filler Particle = %f\n',v); fprintf('X Centroid = %f\n',x1bar); fprintf('Y Centroid = %f\n',x2bar); fprintf('Z Centroid = %f\n',x3bar); fprintf('Mass Moment of Inertia wrt x Centroidal Axis= %f\n',mmix); fprintf('Mass Moment of Inertia wrt y Centroidal Axis= %f\n',mmiy); fprintf('Mass Moment of Inertia wrt z Centroidal Axis= %f\n',mmiz); 242 fprintf('\nTotal %d particle is created \n',(i-2)); ii=2; while ii<=(i-1) fprintf('Sphere No,%d,x1=%f,x2=%f,x3=%f,r=%f\n',ii,x1(ii),x2(ii),x3(ii),radius(ii)); ii=ii+1; end fprintf('\nThanks! for using this program\n') 243 APPENDIX E MATLAB PROGRAM FOR RANDOM COORDINATE GENERATION-IV Following MATLAB program was used to generate random center coordinates of the random size filler particles and for volume fraction more than 25 percent. %********************************************************************** %* %Random Coordinate Generation * %* x1, x2, x3 are three coordinates of the center of the spherical particles * %********************************************************************** %*********************** Defining Number of Particle ********************* % Total Number of Particle = n * %********************************************************************** diary('random_coordinate_random_radius_40_percent_volume_fraction.txt') vv=20; v=0; n=1; i=2; vs=0; tvs=0; tvsn=0; 244 mmix=0; mmiy=0; mmiz=0; x1vs=0; x2vs=0; x3vs=0; mmixn=0; mmiyn=0; mmizn=0; x1vsn=0; x2vsn=0; x3vsn=0; h=250; x1(1)=0.1; x2(1)=0.1; x3(1)=0.1; radius(1)=0.1; fprintf('Random Coordinates of Filler Particles\n'); fprintf('\n') while v<=vv % Start of while #1 [n= no of sphere] %********************* Defining Radius *********************************** r1=(rand(1,1)); r2=(randint(1,1,[17,33])); 245 r=(r1+r2); % Random radius between 17 and 33 s2=r+(0.1*r); %********************** Center Coordinate ************************* a=rand(1,1); % Generating Real Number Between 0-1 b=rand(1,1); c=rand(1,1); d=randint(1,1,[0,h]); % Generating Integer Number Between 0-h e=randint(1,1,[0,h]); f=randint(1,1,[0,h]); p1=abs(a+d); % Generating Real Numbers p2=abs(b+e); p3=abs(c+f); % %******************* Checking for Acceptable Position ****************** % if p1>=s2 % Checking that center is enough away from surface m1=1; % Start of if #1 else m1=0; continue % Go to the top of while loop with i=i end % End of if #1 if p2>=s2 % Start of if #2 m1=1; 246 else m1=0; continue % Go to the top of while loop with i=i end % End of if #2 if p3>=s2 % Start of if #3 m1=1; else m1=0; continue % Go to the top of while loop with i=i end % End of if #3 if abs(p1-h)>=s2 % Start of if #4 m1=1; else m1=0; continue % Go to the top of while loop with i=i end % End of if #4 if abs(p2-h)>=s2 % Start of if #5 m1=1; else m1=0; continue % Go to the top of while loop with i=i end % End of if #5 if abs(p3-h)>=s2 % Start of if #6 247 m1=1; else m1=0; continue % Go to the top of while loop with i=i end % End of if #6 j=1; while j=s1 % Start of if #7 m=1; j=j+1; else m=0; j=i; end % End of if #7 end % End of while #2 if m==1 %start of if #8 x1(i)=p1; x2(i)=p2; x3(i)=p3; 248 radius(i)=r; fprintf('wpave,%f,%f,%f\n',0,0,x3(i)); fprintf('sph4,%f,%f,%f\n',x1(i),x2(i),r); % Calculation of Volume vs=(4*pi*(r^3)/3); tvs=tvs+vs; % Calculation of Centroid x1vs=x1vs+(4*pi*(r^3)/3)*x1(i); x2vs=x2vs+(4*pi*(r^3)/3)*x2(i); x3vs=x3vs+(4*pi*(r^3)/3)*x3(i); % Calculation of Inertia mmix=mmix+(1*(((2/5)*(r*r))+((x1(i)-h/2)*(x1(i)-h/2)))); mmiy=mmiy+(1*(((2/5)*(r*r))+((x2(i)-h/2)*(x2(i)-h/2)))); mmiz=mmiz+(1*(((2/5)*(r*r))+((x3(i)-h/2)*(x3(i)-h/2)))); i=i+1; else i=i; end % End of if #8 vc=h^3; v=100*(tvs/vc); end % End of while#1x1bar=x1vs/tvs; x2bar=x2vs/tvs; x3bar=x3vs/tvs; 249 fprintf('Length of One Side of the Cube, h = %d\n',h); fprintf('Volume of One Sphere = %f\n',vs); fprintf('Total Volume of Sphere = %f\n',tvs); fprintf('Volume of cube, vc = %f\n',vc); fprintf('Volume Fraction of Filler Particle = %f\n',v); fprintf('\n') fprintf('X Centroid = %f\n',x1bar); fprintf('Y Centroid = %f\n',x2bar); fprintf('Z Centroid = %f\n',x3bar); fprintf('\n') fprintf('Mass Moment of Inertia wrt x Centroidal Axis= %f\n',mmix); fprintf('Mass Moment of Inertia wrt y Centroidal Axis= %f\n',mmiy); fprintf('Mass Moment of Inertia wrt z Centroidal Axis= %f\n',mmiz); fprintf('\n') fprintf('\nTotal %d particle is created \n',(i-2)); fprintf('\n') fprintf('\n') ii=2; while ii<=(i-1) fprintf('Sphere No,%d,x1=%f,x2=%f,x3=%f,r=%f\n',ii,x1(ii),x2(ii),x3(ii),radius(ii)); ii=ii+1; end 250 %******************** Shirinking the Cube ****************** fprintf('\n') fprintf('*********************************************\n'); fprintf('*********************************************\n'); fprintf('\n') fprintf('Random Coordinates of Fillers After Shrinking\n'); fprintf('\n') sc=0.7936; hn=h*sc; ii=2; x1new(1)=0.01; x2new(1)=0.01; x3new(1)=0.01; while ii<=(i-1) % Start of while # 3 % Useful coordinates are from i=2 to i=i-1 stopcount=1; x1n=x1(ii)*sc; x2n=x2(ii)*sc; x3n=x3(ii)*sc; rad=radius(ii); 251 %******************* Checking for Acceptable Position ****************** %%*** and Checking Distance From Previously Accepted Particles ************* jj=1; while jjhn % Start of if #41 x1n=hn-randint(1,1,[0,round(rad-5)]); else end % End of if#41 if x2n>hn % Start of if #51 x2n=hn-randint(1,1,[0,round(rad-5)]); 252 else end % End of if #51 if x3n>hn % Start of if #61 x3n=hn-randint(1,1,[0,round(rad-5)]); else end % End of if#61 %********************** Condition 2 in Paper ********************** if x1n==round(rad) x1n=randint(1,1,[0,round(rad-5)]); else end if x2n==round(rad) x2n=randint(1,1,[0,round(rad-5)]); else end if x3n==round(rad) x3n=randint(1,1,[0,round(rad-5)]); else end if x1n==round(hn-rad) x1n=hn-randint(1,1,[0,round(rad-5)]); else end 253 if x2n==round(hn-rad) x2n=hn-randint(1,1,[0,round(rad-5)]); else end if x3n==round(hn-rad) x3n=hn-randint(1,1,[0,round(rad-5)]); else end %******************* Condition 3 in Paper ******************************* dist_iijj=sqrt(((x1n-x1new(jj))*(x1n-x1new(jj)))+((x2n-x2new(jj))*(x2n- x2new(jj)))+((x3n-x3new(jj))*(x3n-x3new(jj)))); s1n=1.07*(rad+radius(jj)); if abs(dist_iijj)>=s1n % Start of if #71 mm=1; jj=jj+1; else mm=0; jj=1; u11=randint(1,1,[1,10]); u22=randint(1,1,[1,10]); u33=randint(1,1,[1,10]); 254 uval=sqrt((u11*u11)+(u22*u22)+(u33*u33)); u1=(u11/uval); u2=(u22/uval); u3=(u33/uval); x1n=x1n+(randint(1,1,[-10,10]))*u1; % Random moving of sphere x2n=x2n+(randint(1,1,[-10,10]))*u2; x3n=x3n+(randint(1,1,[-10,10]))*u3; end % End of if #71 stopcount=stopcount+1; if stopcount >= 500000 % Stopcount if ii=i; jj=ii; fprintf('Stopcount = %d\n',stopcount); else end % End of stopcount if end % End of while #4 if mm==1 % Start of if #81 x1new(ii)=x1n; x2new(ii)=x2n; x3new(ii)=x3n; fprintf('wpave,%f,%f,%f\n',0,0,x3new(ii)); fprintf('sph4,%f,%f,%f\n',x1new(ii),x2new(ii),rad); 255 vs=(4*pi*(rad^3)/3); tvsn=tvsn+vs; x1vsn=x1vsn+(4*pi*(rad^3)/3)*x1new(ii); x2vsn=x2vsn+(4*pi*(rad^3)/3)*x2new(ii); x3vsn=x3vsn+(4*pi*(rad^3)/3)*x3new(ii); mmixn=mmixn+(1*(((2/5)*(rad*rad))+((x1new(ii)-hn/2)*(x1new(ii)-hn/2)))); mmiyn=mmiyn+(1*(((2/5)*(rad*rad))+((x2new(ii)-hn/2)*(x2new(ii)-hn/2)))); mmizn=mmizn+(1*(((2/5)*(rad*rad))+((x3new(ii)-hn/2)*(x3new(ii)-hn/2)))); ii=ii+1; else ii=ii; end % End of if #81 end % End of while # 3 vcn=hn^3; vn=100*(tvsn/vcn); fprintf('\n') fprintf('New Length of One Side of the Cube, hn = %f\n',hn); fprintf('Volume of Sphere, vs = %f\n',vs); fprintf('Total New Volume of Sphere, tvsn = %f\n',tvsn); fprintf('New Volume of cube, vcn = %f\n',vcn); 256 fprintf('New Volume Fraction of Filler Particle = %f\n',vn); fprintf('\n') fprintf('\nTotal %d particle is created \n',(ii-2)); x1barn=x1vsn/tvsn; x2barn=x2vsn/tvsn; x3barn=x3vsn/tvsn; fprintf('New X Centroid = %f\n',x1barn); fprintf('New Y Centroid = %f\n',x2barn); fprintf('New Z Centroid = %f\n',x3barn); fprintf('\n') fprintf('New Mass Moment of Inertia wrt x Centroidal Axis= %f\n',mmixn); fprintf('New Mass Moment of Inertia wrt y Centroidal Axis= %f\n',mmiyn); fprintf('New Mass Moment of Inertia wrt z Centroidal Axis= %f\n',mmizn); fprintf('\n') zzz=2; while zzz<=(i-1) fprintf('Sphere No,%d,x1new,%f,x2new,%f,x3new,%f\n',zzz,x1new(zzz),x2new(zzz),x3new(zzz)); zzz=zzz+1; end fprintf('\nThanks! for using this program\n')