CHARACTERIZATION AND PREDICTION OF MATERIAL RESPONSE
OF MICRO AND NANOUNDERFILLS 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, CoChair
Associate Professor
Mechanical Engineering
______________________________
Jeffrey C. Suhling, CoChair
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 NANOUNDERFILLS 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 NANOUNDERFILLS 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 NANOUNDERFILLS 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 underfillepoxy 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, noflow underfills have very low or no filler content because
micronsize filler particles hinder solder joint formation. Nanosilica underfills have the
potential of attaining higher filler loading in noflow 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 nanosilica underfills with
desirable thermomechanical properties. Elastic Modulus, coefficient of thermal
expansion, bulk modulus, Poisson?s ratio and viscoelastic properties have been predicted.
Temperature dependent thermomechanical properties of nanosilica underfills and
microsilica 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 stressstrain, creep and stress relaxation behavior. The tradeoffs between
using nanofillers instead of micronfillers on thermomechanical 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 fatherinlaw Md. Mahboob Ali and his motherinlaw
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 Multimesh Extrapolation 82
5 PREDICTION OF MATERIAL PROPERTIES MACROSILICA
UNDERFILLS
88
5.1 Underfill Constituent Properties of MacroSilica 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 NANOSILICA
UNDERFILLS BY CONSTANT SIZE FILLER MODELS
110
6.1 Underfill Constituent Properties of NanoSilica 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 NANOSILICA
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 GENERATIONI
217
xi
C MATLAB PROGRAM FOR RANDOM COORDINATE GENERATIONII
224
D MATLAB PROGRAM FOR RANDOM COORDINATE GENERATIONIII
236
E MATLAB PROGRAM FOR RANDOM COORDINATE GENERATIONIV 243
LIST OF FIGURES
1.1 CrossSectional 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 MT200 Testing System with Environmental Chamber
20
3.2
Typical StressStrain Curves of Underfill NUF1. (T = 50 ?C, = 0.001
sec
?&
1
)
23
3.3 Hyperbolic Tangent Model Fit to Typical Underfill StressStrain Data (T
= 50 ?C, = 0.001 sec?&
1
)
24
3.4 Temperature Dependent Average StressStrain 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 NanoFiller 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 NanoFiller 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. (TitaniumSilicate 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 StressStrain 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 StressStrain Curves Generated by FEA with Input
Experimental Data. Length to Width Aspect Ratio = 40
80
4.25 FEA Prediction of StressStrain 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 NanoParticle 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 xCTE by finiteelement analysis of Unit Cell
126
6.9 Prediction of x, y, zCTE by finiteelement 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 finiteelement 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 LogLog 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 MicroSilica Particle
89
6.1 Material Properties for Epoxy Matrix and NanoSilica Particle
111
6.2 Moment of Inertia of Acceptable NanoParticle 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 crosssectional
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 costdominated 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: CrossSectional 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
46
FR4/Glass PCB
1520
Underfill
2530
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 flipchip devices and chipscale packages [1] in a wide variety of
applications including portable consumer electronics like cellular phones [2], laptops [3],
underthehood electronics [4], microwave applications [5], system in package (SIP) [6],
highend 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, noflow underfills have very low
or no filler content because micronsize filler particles hinder solder joint formation.
Nanosilica underfills have the potential of attaining higher filler loading in noflow
underfills without hindering solder interconnect formation [9, 10]
Nanosilica imparts the same modulus enhancement and CTE reduction to
adhesives as microsized silica particles. However, nanosilica does not settle in an
underfill formulation and does not interfere with the solder interconnection process,
unlike micronsized particles. Low fillerloading in noflow underfills causes the
coefficient of thermal expansion to be higher than capillary flow underfills. The noflow
underfills are largely unfilled or have very low filler loading because micronsized filler
particles interfere with the solder interconnection process [9, 10]. Noflow underfills are
an attractive alternative to the use of capillary flow underfills, which typically require a
postreflow batch cure operation adding to production cycle time. In addition, nanosilica
particles can achieve much higher volume fraction loading than micronsized 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 thermomechanical reliability of flipchip devices.
There is need for propertyprediction techniques for formulation of the underfills
with desirable thermomechanical 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 implicitfinite 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 thermomechanical properties of
nanounderfills 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 flipchip
device and the printed circuit assembly. Material properties in the present effort have
been validated versus properties measured on micromechanical test structures. The test
13
structures are uniaxialtension specimen with 75125 ?m (35 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 stressstrain data. Strain gages were used to measure
the CTE of cured underfill. Temperature dependent material behaviors of microsilica
underfills were also evaluated.
1.2 Research Objective
? Use micro tension/torsion testing machine to evaluate temperature dependent
stressstrain, creep and relaxation behavior of nanosilica and microsilica
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
nanosilica 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 nanounderfill
? 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 noflow underfills are largely unfilled or have
very low filler loading because micronsized filler particles interfere with the solder
interconnection process. Nanosilica underfills have the potential of attaining higher
filler loading in noflow underfills without hindering solder interconnect formation.
Nanosilica underfill is a very new concept in underfill family. In these underfills
nanosilica particles are dispersed in an epoxy matrix. Adding nanosilica particle
instead of microsilica 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 nanosilica underfills. In recent days
development and application of nanosilica underfills are very important to achieve the
success in electronic packaging and manufacturing industries. Knowledge of constitutive
behavior of nanosilica underfill is very critical too. No study has been performed yet to
characterize the thermomechanical properties of nanosilica underfill. The tradeoffs for
moving towards nanosilica underfill instead of using the microsilica underfill is not well
understood too.
18
There is need for propertyprediction techniques for formulation of the underfills
with desirable thermomechanical characteristics. Finite element and analytical
techniques need to be developed to predict properties of nanosilica underfills.
Thermomechanical, elastic, viscoelastic, viscoplastic properties of nanosilica 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 MT200 tension/torsion thermomechanical 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 6axis 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 MT200, forces and displacements are measured. The axial
stress and axial strain are calculated from the applied force and measured crosshead
displacement using
LL
L
A
F ?
=
?
=?=?
(3.1)
19
Figure 3.1: MT200 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 crosssectional area, ? is the measured crosshead 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 nanounderfill 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
nanounderfill 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 75125 ?m (35 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.
StressStrain Data
Figure 3.2 shows typical stressstrains 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 stressstrain 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 threeparameter hyperbolic tangent model has been
used to model the observed nonlinear underfill stressstrain data. Such a model has been
21
used historically to model the stressstrain 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 stressstrain 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 stressstrain
curves in a set should be fit simultaneously in order to obtain the best set of hyperbolic
tangent model material constants.
The stressstrain 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 StressStrain 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 StressStrain Data (T = 50
?C, ?& = 0.001 sec
1
)
24
Typical variation of the stressstrain 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 stressstrain 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
nanounderfill NUF1 versus micronfiller 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 StressStrain 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 NanoFiller 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 MT200 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 recirculating 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 stressstrain 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 nanounderfill 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 75125 ?m (35 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.
StressStrain Data
Typical variation of the stressstrain 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 stressstrain 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 nanounderfill NUF1,
NUF2 versus micronfiller 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
)
AAA01
39
NUF2
NUF1
E
l
ast
i
c M
odu
l
u
s,
E
(G
P
a
)
Figure 3.14: Elastic Modulus of NanoFiller 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 TitaniumSilicate 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. (TitaniumSilicate 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 sixaxis micro tension/torsion machine has
been used to perform tensile, creep tests.
3.2.1: Effect of Extreme Low Temperatures
Accurate evaluation of the stressstrain 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 stressstrain 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, stressstrain and creep tests have been
performed on several underfill materials down to cryogenic temperatures. For example,
stressstrain 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 stressstrain 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
110 ?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 StressStrain 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 stressstrain 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 stressstrain curve
and the stressstrain curve generated by the FEA output. The generated stressstrain
curve from FEA output is much different than the input experimental stressstrain 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
stressstrain curve by FEA output is almost same as the experimental stressstrain 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 stressstrain 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 stressstrain curves generated by the FEA
output data. The stressstrain 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 thermomechanical 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 stressstrain 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 stressstrain curves at all
temperatures are exactly same as the stressstrain curves generated by FEA output. Due
to this nice similarity between the experimental and FEA stressstrain curves it will be
possible to make FEA prediction of the stressstrain 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 StressStrain 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 StressStrain Curves at 110 ?C. Length to Width Aspect
Ratio = 40
81
4.3: Error Estimation by Multimesh 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 MACROSILICA UNDERFILLS
5.1: Underfill Constituent Properties of MacroSilica 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 MicroSilica 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 Undeformed 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 Undeformed 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
microsilica 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 ? + 1E06
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 ? + 9E06
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 NANOSILICA
UNDERFILLS BY CONSTANT SIZE FILLER MODELS
6.1 Underfill Constituent Properties of NanoSilica Underfills
Nanosilica 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 nanounderfill 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 nanosilica underfills with desirable thermomechanical properties.
111
Table 6.1 Material Properties for Epoxy Matrix and NanoSilica 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 finiteelement models of threedimensional
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
nanosilica 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 momentofinertia for
the nanoparticles 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 NanoParticle Distribution
117
118
Table 6.2: Moment of Inertia of Acceptable NanoParticle 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 nanounderfill 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 xaxis 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 xCTE by finiteelement 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, zCTE by finiteelement 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 xdirection 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 stressstrain 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 NANOSILICA
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 nanounderfill 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 timedependent properties of nano
underfills investigated in this paper is new. Previously authors have investigated the
application of RSA to prediction of elastic properties of nanounderfills [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 timedependent 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 finiteelement models of threedimensional 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 nanounderfill, 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 builtin 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 nanosilica 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 nanounderfill 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
xaxis 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 finiteelement analysis of Unit Cell
144
this analysis, different boundary conditions than the CTE models have been used.
Displacement in the xdirection 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 stressstrain
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 nanounderfill 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 nanounderfill. 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 xdirection 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 WilliamsLandelFerry (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 loglog 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. WilliamsLandelFerry (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
LandelFerry 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. Stressstrain 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: LogLog 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 ydirection 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 stressstrain 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 stressstrain 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 stressstrain 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
)
Experimental20% Filler Volume Fraction
FEA Prediction20% Filler Volume Fraction
Experimental10% Filler Volume Fraction
FEA Prediction10% Filler Volume Fraction
Epoxy Only
120
Figure 7.22: Prediction of Elastic Modulus Relaxation
173
174
which has 20% nanosilica 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
Experimental20% Filler Volume Fraction
FEA Prediction20% Filler Volume Fraction
Experimental10% Filler Volume Fraction
FEA10% 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
Experimental20% Filler Volume Fraction
FEA Prediction20% Filler Volume Fraction
Epoxy Only
Experimental10% Filler Volume Fraction
FEA10% 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 nonlinear material properties of underfill materials. Model predictions of properties
have been validated with experimental data. The multifiller models with randomly
oriented, random sized spherical filler particles predict the properties more accurately
than one or eight filler models with periodic filler distribution.
Multimesh extrapolation has been used to model macro behavior in uniaxial test.
Finite element results indicate that the lengthtowidth 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?, IEEECPMT Vol.27, Issue:3, pp. 461 ? 467,
Sept.2004
3 Pascariu G., Cronin P, Crowley D, ?Nextgeneration 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.8291, 1998
5 Bedinger, John M., ?Microwave Flip Chip and BGA Technology?, IEEE MTTS
International Microwave Symposium Digest, v 2, pp 713716, 2000
6 Van den Crommenacker, J., ?The SysteminPackage Approach?, IEEE
Communications Engineer, Vol 1, No. 3, pp. 2425, 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 319324,
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. 7479, March 2001
178
10 Shi, S. H.; Wong, C. P., ?Recent Advances in the Development of NoFlow
Underfill Encapsulants ? A Practical Approach towards the Actual
Manufacturing Application?, IEEE Transactions on Electronics Packaging
Manufacturing, Vol. 22, pp. 331339, 1999
11 Drugan, W. J., ?MicromechanicsBased Variational Estimates for a Higher
Order Nonlocal Constitutive Equation and Optimal Choice of Effective Moduli
of Elastic Composites,? J. Mech. Phys. Solids, 48, pp 13591387, 2000
12 Drugan, W. J., Wills, J. R., ?A MicromechanicsBased Nonlocal Constitutive
Equation and Estimates of the Representative Volume Element Size for Elastic
Composites,? J. Mech. Phys. Solids, 44, pp 497524, 1996
13 Segurado, J., Llorca, J., ?A Numerical Approximation to Elastic Properties of
SphereReinforced Composites,? Journal of Mechanics and Physics of Solids,
Vol. 50, pp. 21072121, 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. 437444,
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. 467476,
September 2005
179
21 Lall, P., Islam, S., Suhling, J. C., Tian, G., ?NanoUnderfills for HighReliability
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. 415429, 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 AbdelHady, 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 5131, Measurement Group, Technical Note, 2004
30 Pyrz, R., ?Correlation of Microstructure Variability and Local Stress Field in
TwoPhase Materials,? Mater. Sci. Eng., Mater. Sci. Eng., A 177, pp 253259,
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 SysteminPackage Approach?, IEEE
Communications Engineer, Vol 1, No. 3, pp. 2425, 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. 695711, Kona,
Hawaii, September 2226, 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
EEPvol. 192, 1997a, pp. 12911298
35 Merton M. M., Mahajan R. L., and Nikmanesh N., ?Alternative Curing Methods
for FCOB Underfill? Advances in Electronic Packaging, ASME EEP, vol. 191,
1997, pp. 291300
36 Qian, Z., Wang, J., Yang, J., Liu, S., ?ViscoElasticPlastic Properties and
Constitutive Modeling of Underfills?, IEEE Transactions on Components and
Packaging Technology?, vol. 22, No. 2, pp. 152157, June 1999
37 Shi, S. H., Yao, Q. Qu, J., Wong, C. P., ?Study on The Correlation of FlipChip
Reliability with Mechanical Properties of NoFlow Underfill Materials?,
Proceedings of International Symposium on Advanced Packing Materials, pp.
271277, March 2000
38 Schubert A., 1997, ?ThermoMechanical Reliability Analysis of FlipChip
Assemblies by Combined Microdac and the Finite Element Method?, Advances
in Electronic Packaging, ASME EEPVol. 192, 16471654
39 M. Lu, W. Ren and S. Liu, ?A MultiAxial ThermoMechanical Fatigue Tester
for Electronic Packaging Material?, EEPVol. 17, ?ASME Symposium of
Censoring, monitoring, modeling and verification of Electronic Packaging, Nov.
1721, 1996, Atlanta, pp 8792
40 M. Lu, Z. Qian and S Liu, ?Unified MultiAxial Submicron Fatigue Tester for
Miniaturized Specimen?, EEPVol. 191, Advances in Electronic Packaging
1997, ASME 1997, Volume 1, pp 11451150
41 Ren, W., ?ThermoMechanical 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. 191, 291300
181
43
Widom, B., ?Random Sequential Addition of Hard Spheres to a Volume,? The
Journal of Chemical Physics, Vol. 44, No. 10, pp. 38883894, 1966
44
Rintoul, M. D., Torquato, S., ?Reconstruction of the Structure of Dispersions,?
Journal of Colloid and Interface Science, 186, pp. 467476, 1997
45
Bohm, H.J., Han, W., ?Comparisons between threedimensional and two
dimensional multiparticle 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
TwoPhase Materials,? Material Science and Engineering, A177, pp. 253259,
1994
49 Pyrz, R., ?Quantitative Description of the Microstructure of Composites. Part I:
Morphology of Unidirectional Composite Systems,? Composite Science and
Technology,50, pp. 197208, 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 nanosilica 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 GENERATIONI
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 01
b=rand(1,1);
c=rand(1,1);
d=randint(1,1,[0,h]); % Generating Integer Number Between 0h
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(p1h)>=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(p2h)>=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(p3h)>=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',(i2));
ii=2;
while ii<=(i1)
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 GENERATIONII
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 01
b=rand(1,1);
c=rand(1,1);
d=randint(1,1,[0,h]); % Generating Integer Number Between 0h
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(p1h)>=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(p2h)>=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(p3h)>=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',(i2));
fprintf('\n')
zz=2;
while zz<=(i1)
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=i1
x1new(1)=0.01;
x2new(1)=0.01;
x3new(1)=0.01;
while ii<=(i1) % 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=hnrandint(1,1,[0,(r5)]);
else
232
end % End of if#41
if x2n>hn % Start of if #51
x2n=hnrandint(1,1,[0,(r5)]);
else
end % End of if #51
if x3n>hn % Start of if #61
x3n=hnrandint(1,1,[0,(r5)]);
else
end % End of if #61
dist_iijj=sqrt(((x1nx1new(jj))*(x1nx1new(jj)))+((x2nx2new(jj))*(x2n
x2new(jj)))+((x3nx3new(jj))*(x3nx3new(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',(ii2));
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<=(i1)
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 GENERATIONIII
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 01
b=rand(1,1);
c=rand(1,1);
d=randint(1,1,[0,h]); % Generating Random Integer Number Between 0h
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(p1h)>=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(p2h)>=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(p3h)>=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',(i2));
ii=2;
while ii<=(i1)
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 GENERATIONIV
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 01
b=rand(1,1);
c=rand(1,1);
d=randint(1,1,[0,h]); % Generating Integer Number Between 0h
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(p1h)>=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(p2h)>=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(p3h)>=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',(i2));
fprintf('\n')
fprintf('\n')
ii=2;
while ii<=(i1)
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<=(i1) % Start of while # 3
% Useful coordinates are from i=2 to i=i1
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=hnrandint(1,1,[0,round(rad5)]);
else
end % End of if#41
if x2n>hn % Start of if #51
x2n=hnrandint(1,1,[0,round(rad5)]);
252
else
end % End of if #51
if x3n>hn % Start of if #61
x3n=hnrandint(1,1,[0,round(rad5)]);
else
end % End of if#61
%********************** Condition 2 in Paper **********************
if x1n==round(rad)
x1n=randint(1,1,[0,round(rad5)]);
else
end
if x2n==round(rad)
x2n=randint(1,1,[0,round(rad5)]);
else
end
if x3n==round(rad)
x3n=randint(1,1,[0,round(rad5)]);
else
end
if x1n==round(hnrad)
x1n=hnrandint(1,1,[0,round(rad5)]);
else
end
253
if x2n==round(hnrad)
x2n=hnrandint(1,1,[0,round(rad5)]);
else
end
if x3n==round(hnrad)
x3n=hnrandint(1,1,[0,round(rad5)]);
else
end
%******************* Condition 3 in Paper *******************************
dist_iijj=sqrt(((x1nx1new(jj))*(x1nx1new(jj)))+((x2nx2new(jj))*(x2n
x2new(jj)))+((x3nx3new(jj))*(x3nx3new(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',(ii2));
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<=(i1)
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')