Time-domain viscoelastic constitutive model based on concurrent fitting of frequency-domain characteristics

MR-11731; No of Pages 9 Microelectronics Reliability xxx (2015) xxx–xxx

Contents lists available at ScienceDirect

Microelectronics Reliability journal homepage: www.elsevier.com/locate/mr

Time-domain viscoelastic constitutive model based on concurrent fitting of frequency-domain characteristics Tz-Cheng Chiu ⁎, Bo-Sheng Lee, Dong-Yi Huang, Yu-Ting Yang 1, Yi-Hsiu Tseng 1 Department of Mechanical Engineering, National Cheng Kung University, Tainan 701, Taiwan, ROC

a r t i c l e

i n f o

Article history: Received 10 May 2015 Received in revised form 30 June 2015 Accepted 9 July 2015 Available online xxxx Keywords: Viscoelasticity Storage modulus Relaxation modulus Interconversion Poisson's ratio Warpage

a b s t r a c t A numerical procedure for constructing the multiaxial viscoelastic model for polymeric packaging material over a wide range of temperature is presented. By using the proposed best-fitting procedure, experimentally measured frequency-domain Young's and shear storage moduli are used to calculate the time-domain bulk and shear relaxation moduli which describe the three-dimensional constitutive behavior of a viscoelastic solid. The numerical procedure incorporates restrictions that ensure that the derived time-domain material function is physics compatible. The proposed procedure was applied to construct the viscoelastic constitutive models of epoxy molding compounds (EMCs), and compared to results obtained by using approximate-formula based direct conversion procedure. It was shown that, without using the proposed procedure, the directly calculated time-dependent Poisson's ratio oscillates significantly in the rubbery regime and is physically inadmissible. To validate the constitutive model constructed by using the proposed procedure, a numerical finite element model that incorporates the viscoelastic constitutive model of the EMC was applied to simulate warpage of an overmolded package under the solder reflow process and compared to experimental shadow Moiré measurements. © 2015 Elsevier Ltd. All rights reserved.

1. Introduction In organic microelectronic assemblies the viscoelastic nature of polymeric packaging materials plays an important role in the overall deformations and stresses because the process temperatures are typically around or higher than the glass transition temperatures (Tg) of the polymeric constituents [1]. Evaluations of the influences of these timedependent behaviors on the thermomechanical characteristics of the assemblies can be performed by using numerical finite element procedures. In these numerical models, viscoelastic behaviors of the polymeric materials are described by using time-dependent material functions. Within the context of linear viscoelasticity, any two of the relaxation (or creep) functions corresponding to the uniaxial, shear, bulk and Poisson's responses are sufficient to define the three-dimensional constitutive behavior of the viscoelastic solid, and the remaining material functions can be obtained from the two selected functions. Experimental measurements of the viscoelastic behaviors are typically performed either under quasi-static relaxation (or creep) conditions or under time harmonic oscillating conditions. Aside from the standard quasi-static uniaxial or shear tests, tools such as dynamic mechanical analysis (DMA) instrument under cyclic tensile or bending conditions [2], rheometer under cyclic torsion condition [3], and high pressure dilatometer under compressive creep condition [4] are often used for the viscoelasticity characterization. Among the various tools ⁎ Corresponding author. Central Development Engineering, ASE Group Chung-Li, Taoyuan 320, Taiwan, ROC.

1

used for measuring viscoelastic properties, the DMA instrument and rheometer are the most widely adopted because the test durations required in these time harmonic oscillating setups are much shorter than those required in the quasi-static creep or relaxation tests. The viscoelastic material functions obtained from these oscillation tests are typically in the forms of Young's/shear storage and loss moduli as functions of oscillating frequency and temperature. However, in most numerical finite element procedures used for mechanical simulations, the viscoelastic formulations are based on the shear and bulk relaxation functions. As a consequence, additional post-processes are required to convert the Young's and shear storage moduli in the frequency domain to the storage and shear relaxation moduli in the time domain. The post-process would involve converting the frequency-domain viscoelastic functions to the corresponding time-domain functions, and calculating the viscoelastic bulk relaxation modulus from the uniaxial and shear responses. The viscoelastic Poisson's ratio, while is not required for the finite element simulation, should also be calculated in the post-process to assess the accuracy of the multiaxial viscoelastic constitutive behavior described by the corresponding bulk and shear relaxation moduli. Interconversions of the linear viscoelastic material functions between their values in the frequency and time domains can be achieved by using approximate analytical formulas [5,6] or numerical procedures [7,8]. In general these procedures were shown to give relatively good estimations. On the other hand, the calculation of the bulk and the Poisson's ratio functions are more complicated. It was highlighted by Tschoegl et al. [9] that the smallest specimen-to-specimen property

http://dx.doi.org/10.1016/j.microrel.2015.07.031 0026-2714/© 2015 Elsevier Ltd. All rights reserved.

Please cite this article as: T.-C. Chiu, et al., Time-domain viscoelastic constitutive model based on concurrent fitting of frequency-domain characteristics, Microelectronics Reliability (2015), http://dx.doi.org/10.1016/j.microrel.2015.07.031

T.-C. Chiu et al. / Microelectronics Reliability xxx (2015) xxx–xxx

variation and differences in the environment variables would result in significant variation or even physically contradicting Poisson's ratio values, and it is required to measure the two independent viscoelastic functions simultaneously on the same specimen in a single experiment to minimize this issue. While the protocol of simultaneously measuring two material functions are essential for the accuracy of the viscoelastic Poisson's ratio, it is very difficult to accomplish in practice, and the reported data are scarce [9]. As an alternative to this rigorous experimental protocol, a numerical best-fitting procedure with restriction on Poisson's ratio's acceptability is proposed to post-process the separately measured frequency-domain Young's and shear storage moduli for determining a set of time-domain bulk and shear relaxation moduli that would meet the requirements for physically admissible multiaxial viscoelastic constitutive model. In this paper multiaxial viscoelastic constitutive models obtained by using the proposed moduli interconversion procedure are presented for EMCs of two different filler percentages. Experimental characterizations for the frequency dependent Young's and shear storage moduli under various temperatures were first performed by using a DMA instrument and a rheometer, respectively. Mastercurves of the experimental Young's and shear storage moduli were constructed based on the time-temperature-superposition principle. A numerical program based on the proposed viscoelastic function interconversion procedure was then applied to convert the Young's and shear storage moduli to the corresponding timedomain viscoelastic material functions. The mathematical models of the bulk and shear relaxation moduli in the forms of Prony series were implemented in finite element models for simulating the warpage of electronic package under surface mount solder reflow condition, and compared to shadow Moiré experimental measurements to evaluate the accuracy of the viscoelastic model. 2. Multiaxial viscoelastic model

t −

0

κ ðt−τ Þ

dεkk ðτ Þ dτ þ dτ

t −

0

2μ ðt−τÞ

dei j ðτ Þ dτ; i; j ¼ 1; 2; 3; dτ ð1Þ

where σ ij is the stress, εkk is the sum of normal strains (εkk = ε 11 + ε22 + ε33 ), eij is the deviatoric strain, δ ij is the Kronecker delta, and κ and μ are the bulk and shear relaxation moduli, respectively. Under uniaxial loading condition, the multiaxial constitutive Eq. (1) reduces to Z σ 11 ðt Þ ¼

t −

0

Eðt−τÞ

dε11 ðτÞ dτ; dτ

  σ 11 ðt Þ ¼ E0 ðωÞ þ iE″ðωÞ ε0 eiωt ;

ð3Þ

where ε0eiωt denotes the strain excitation, E′(ω) and E″(ω) are the Young's storage and loss moduli, respectively, and are given by E0 ðωÞ ¼ E∞ þ ω Z E″ðωÞ ¼ ω

0

∞

0

140 C

15

150 C

10

160 C

5

170 C

220 C

0 0

20

40

60

80

100

120

(1/s) Fig. 1. Young's storage modulus of Compound A.

where E∞ is the value of E as t → ∞. The inverse relationship of (4) can be expressed as [10] Eðt Þ ¼

2 π

Z

∞

0

E0 ðωÞ

sin ωt dω: ω

ð6Þ

For the case of shear loading such as the time-harmonic shear strain excitation used in torsional DMA, the shear stress-shear strain and the related time-frequency domain moduli relationships can be obtained by replacing the normal stress, normal strain, and Young's modulus terms in Eqs. (2)–(6) with the shear stress, shear strain, and shear modulus terms, respectively. The relationship between shear relaxation and storage moduli can be therefore expressed as 2 π

Z

∞

0

μ 0 ðωÞ

sin ωt dω: ω

ð7Þ

2.1. Direct interconversion Because the viscoelastic stress–strain relationships such as Eqs. (1) and (2) are in the forms of convolution integrals, relationships between the material functions cannot be expressed in simple equations. Interconversions of these materials functions including Young's relaxation modulus, viscoelastic Poisson's ratio, bulk relaxation modulus and shear relaxation modulus are typically conducted in the Laplace transform domain, in which the relationships between the material functions are in the forms of algebraic equations. Consequently, a straight-forward approach for determining the time-domain bulk relaxation modulus and viscoelastic Poisson's ratio from the frequency-

ð2Þ

where E is the Young's relaxation modulus. When the viscoelastic material is subjected to uniaxial time harmonic strain excitation such as in the DMA, the steady-state stress response can be written as

Z

20

20

40 30 10 0 10 20 40 50

15

' (GPa)

σ i j ðt Þ ¼ δi j

Z

40 C 60 C 80 C 90 C 100 C

25

μ ðt Þ ¼

The multiaxial constitutive behavior of an isotropic linearviscoelastic material can be expressed in a convolution integral form given by Z

30

' (GPa)

2

C C C C C C C C

10 130 C

140 C

5

150 C ∞

½Eðt Þ−E∞  sinðωt Þdt;

ð4Þ

0

270 C

0

50

100

150

200

250

300

(1/s) ½Eðt Þ−E∞  cosðωt Þdt;

ð5Þ

Fig. 2. Young's storage modulus of Compound B.

Please cite this article as: T.-C. Chiu, et al., Time-domain viscoelastic constitutive model based on concurrent fitting of frequency-domain characteristics, Microelectronics Reliability (2015), http://dx.doi.org/10.1016/j.microrel.2015.07.031

T.-C. Chiu et al. / Microelectronics Reliability xxx (2015) xxx–xxx

25

10

30 C 40 C

8

140 C 150 C

6

storage modulus (GPa)

' (GPa)

12

160 C

4 170 C

2 0

3

260 C 270 C

0

20

40

60

80

100

20

E' 15 10

'

5 0 10-20

120

10-15

10-10

domain Young's and shear storage moduli is by: (i) converting the Young's and shear storage moduli to their corresponding relaxation moduli in the Laplace transform domain, (ii) calculating the bulk relaxation modulus and viscoelastic Poisson's ratio in the Laplace transform domain, and (iii) inverting these Laplace-domain functions to obtain their time-domain values. This procedure can be easily implemented by using the approximate analytical formulas proposed by Schapery and Park [5]. The conversion of the frequency-domain Young's storage modulus to the Laplace-domain Young's relaxation modulus can be performed by using the approximate formula [5]:  E0 ðωÞ  ; cosðnπ=2Þω¼s

ð8Þ

where Ẽ denotes the Young's relaxation modulus in the Laplace–Carson transform domain, and is defined as ∞

0

Eðt Þe−st dt;

ð9Þ

in which Ē(s) is the standard Laplace transform of the Young's relaxation modulus. In Eq. (8), n is the slope of the logarithmic E′ — logarithmic ω mastercurve, i.e.,

1010

1015

 d log10 E0 ðωÞ : n¼ dð log10 ωÞ

Calculation of the shear relaxation modulus in the Laplace–Carson transform domain from the experimentally measured shear storage modulus can also be done by using the same Eqs. (8) and (10), with E replaced by μ. By following to the elastic–viscoelastic correspondence principle, the bulk relaxation modulus and the viscoelastic Poisson's ratio in the Laplace–Carson transform domain can be expressed as κ~ ðsÞ ¼

~ ðsÞ ¼ ν

ð10Þ

~ðsÞ μ~ ðsÞ E

~ðsÞ 9μ~ ðsÞ−3E

;

ð11Þ

~ðsÞ E −1; 2 μ~ ðsÞ

ð12Þ

respectively. The time-domain bulk relaxation modulus can be obtained by using an approximate formula for the inverse Laplace–Carson transform given by [5] κ ðt Þ ≃

 κ~ ðsÞ  ; Γ ð1−mÞt¼1=s

ð13Þ

where Γ denotes the Gamma function, and m is given by



m¼

dð log10 κ~ ðsÞÞ : dð log10 sÞ

ð14Þ

25

5

shift factor WLF model (>120 ° C) WLF model (<120 ° C)

20 40 C

4

15 120 C 130 C

3 2

log(aT)

' (GPa)

105

Fig. 5. Young's and shear storage modulus mastercurves for Compound A, reference temperature: 140 °C.

Fig. 3. Shear storage modulus of Compound A.

Z ~ðsÞ ¼ sEðsÞ ¼ s E

100

(1/s)

(1/s)

~ðsÞ ≃ E

10-5

10 5 0

140 C

-5

1 0

150 C 270 C

0

20

40

60

80

100

(1/s) Fig. 4. Shear storage modulus of Compound B.

120

-10 -15

0

50

100

150 (°C)

200

250

300

Fig. 6. Temperature shift factors for storage modulus mastercurves of Compound A, reference temperature: 140 °C.

Please cite this article as: T.-C. Chiu, et al., Time-domain viscoelastic constitutive model based on concurrent fitting of frequency-domain characteristics, Microelectronics Reliability (2015), http://dx.doi.org/10.1016/j.microrel.2015.07.031

4

T.-C. Chiu et al. / Microelectronics Reliability xxx (2015) xxx–xxx

16

Table 1 WLF shift function parameters for the EMCs.

storage modulus (GPa)

14 12

E'

Compound A

10

Above 120 °C Below 120 °C Above 120 °C Below 120 °C

Compound B

8 6

Tref (°C)

C1

C2

140 167 120 120

71.9 1460 41.6 210

324.3 15,427 148.6 1625

'

4 where the subscript g denotes the glassy modulus value, w and τ are the weighting constant and relaxation time, respectively, and

2 0 10-20

10-10

100

1010

1020

1030

w∞ ¼ 1−

(1/s)

X

wi :

ð18Þ

Fig. 7. Young's and shear storage modulus mastercurves for Compound B, reference temperature: 120 °C.

By substituting Eq. (17) into Eq. (4) and observing E∞ = EgwE∞, the Prony series representation of the Young's storage modulus is given by

The approximate formula (13) can also be applied to estimate other viscoelastic material functions from their values in the Laplace–Carson transform domain.

E0 ðωÞ ¼ Eg wE∞ þ

2.2. Numerical best-fitting An issue that may arise in the direct interconversion for viscoelastic material functions is that small variations in the material functions from separate characterization experiments may lead to significant errors in the derived material functions [9]. A general approach that may be able to avoid this issue is to express the to-be-derived material functions in a prescribed mathematical form, and to assign the physics requirements as the restriction conditions in the experimental data best fitting process such that a coherent set of material constants for the prescribed mathematical model can be obtained. The mathematical functions typically used for representing the viscoelastic relaxation moduli are the Prony series: " κ ðt Þ ¼ κ g wκ∞ þ

M X i¼1

" μ ðt Þ ¼ μ g wμ∞ þ

N X

t μ τi

wEi exp −

t τEi

μ

i¼1

" Eðt Þ ¼ Eg wE∞ þ

P X i¼1

wEi 

i¼1

ω2  E 2

1=τi

# þ ω2

:

ð19Þ

Similarly, the Prony series representation of the shear storage modulus can also be expressed in the same mathematical form as Eq. (19), with the model constants replaced with those given in Eq. (16). For particulate composite such as the EMCs, typically the viscoelastic Poisson's ratio falls in the range of 0 and 0.5, and is a monotonically increasing function of time [9]. A Prony series representation of the viscoelastic Poisson's ratio that would satisfy these requirements can be expressed as   Q X ξ νðξÞ ¼ ν∞ − ν i exp − ν ; τi i¼1

ð15Þ

!#

Q X νi ¼ νi N 0; ν ∞ − i¼1

;

ð16Þ

ð20Þ

;

ð17Þ

!#

! Eg −1 N 0; ν ∞ ¼ 2μ g

! Eg wE∞ μ −1 ≤ 0:5: ð21Þ 2μ g w∞

Jerabek et al. [11] showed that, without loss of generality, the same number of Prony series terms and relaxation time constants can be used in the relaxation moduli (16), (17) and the viscoelastic Poisson's ratio (Eq. (18)), i.e., N = P = Q and τμi = τEi = τνi = τi. By taking these assumptions, converting Eqs. (16), (17) and (20) into

30

25

shift factor WLF model (>120 ° C) WLF model (<120 ° C)

20

P X

where ν∞ is the large-time value of the Poisson's ratio, νi and τνi are material constants, and

 # t κ wi exp − κ ; τi

wi exp −

"

20

Compound A Compound B

(GPa)

log(aT)

10 0

15

-10

10

-20

5

-30

0

100

200

300

(°C)

0 10-30

10-20

10-10

100

1010

1020

(s) Fig. 8. Temperature shift factors for storage modulus mastercurves of Compound B, reference temperature: 120 °C.

Fig. 9. Bulk relaxation modulus mastercurve from direct interconversion.

Please cite this article as: T.-C. Chiu, et al., Time-domain viscoelastic constitutive model based on concurrent fitting of frequency-domain characteristics, Microelectronics Reliability (2015), http://dx.doi.org/10.1016/j.microrel.2015.07.031

T.-C. Chiu et al. / Microelectronics Reliability xxx (2015) xxx–xxx

5

25

relaxation modulus (GPa)

0.50

0.25

0.00

-0.25

Compound A Compound B

-0.50 -20 10

10-10

100

1010

1020

20

E 15 10 5 0 10-15

1030

s (1/s)

10-10

10-5

100

105

1010

1015

1020

(s)

Fig. 10. Laplace–Carson-transform-domain Poisson's ratio mastercurve from direct interconversion.

the Laplace–Carson transform domain, and then substituting into Eq. (12), it gives

Fig. 12. The fitted relaxation modulus mastercurve of Compound A, reference temperature 140 °C.

3. Constitutive model of EMC 3.1. Experimental characterization

N X ν∞ − i¼1

νi s − s þ 1=τi

Eg E w −wμ∞ 2μ g ∞

! þ μ wi

N X i¼1 N X

þ

i¼1

Eg E μ w −wi 2μ g i

!

s s þ 1=τi

μ

wi s s þ 1=τi

¼ 0:

ð22Þ

By substituting the collocation points s = 1/τ1, 1/τ2, …, 1/τN into Eq. (22), an equivalent set of N simultaneous equations can be used to determine ν1, ν2, …, νN in terms of wE1, wE2, …, wEN and wμ1, wμ2, …, w-N μ , and then Eq. (21) can be rewritten as N + 2 restriction equations for wE1, wE2, …, wEN and wμ1, wμ2, …, wμN. Consequently, the task of developing a coherent set of viscoelastic material functions can be achieved by obtaining the Prony series constants of Eqs. (16) and (17) such that the experimentally measured Young's and shear storage moduli are least square fitted by the Prony series representations in the frequency domain (e.g., (Eq. (19)) for E′(ω)) with the N + 2 restriction equations generated from Eqs. (21) to (22). A numerical implementation of this procedure was developed by using sequential-quadratic-programming based optimization method in MATLAB [13].

0.5

storage modulus (GPa)

25 20

Experiment Fitted

0.4

E'

15

0.3

'

10

0.2 0.1

5 0 10-15

In this study two commercial EMCs are investigated: Compound A is 88 wt.% silica filled and is used in typical overmolded laminate package, Compound B is a 70 wt.% filled compound and is used for encapsulating ultra-thin packages. Measurements of the Young's storage moduli for both compounds were conducted by using a TA Instruments Q800 DMA tester. As cured EMC specimens with cross-sectional dimensions of 2 mm × 0.4 mm and a gauge length of 20 mm were characterized under time-harmonic tension loads. The tests were conducted under frequency sweep between 0.1 and 40 Hz at constant temperatures between −60 and 270 °C. Experimental results of the Compounds A and B Young's storage moduli are shown in Figs. 1 and 2, respectively. Shear characterizations of the EMCs were conducted by using an Anton Paar MCR-302 rheometer. As cured specimens of gauge dimension 45 mm × 6 mm × 2 mm were tested under time harmonic torsion configuration over the frequency and temperature ranges from 0.1 to 100 rad/s and from − 70 to 270 °C, respectively. The experimental shear storage moduli at various temperatures for Compounds A and B are shown in Figs. 3 and 4, respectively. Comparing the experimental results for Compound A in Figs. 1 and 3, it can be seen that the temperature ranges of glassy to rubbery transition observed from the tension and torsion DMA tests are in good agreement. Similar observation can also be made for the Compound B results (Figs. 2 and 4).

10-10

10-5

100

105

1010

1015

(1/s) Fig. 11. The fitted storage modulus mastercurve of Compound A, reference temperature 140 °C.

0.0 -15 10

10-10

10-5

100

105

1010

1015

1020

(s) Fig. 13. The fitted viscoelastic Poisson's ratio mastercurve of Compound A, reference temperature 140 °C.

Please cite this article as: T.-C. Chiu, et al., Time-domain viscoelastic constitutive model based on concurrent fitting of frequency-domain characteristics, Microelectronics Reliability (2015), http://dx.doi.org/10.1016/j.microrel.2015.07.031

6

T.-C. Chiu et al. / Microelectronics Reliability xxx (2015) xxx–xxx

Assuming that the EMCs are thermorheologically simple, the mastercurves of the Young's and shear storage moduli for the EMCs can be constructed based on the time-temperature superposition principle. Shown in Figs. 5 and 6 are the storage modulus mastercurves and the corresponding temperature shift factors, respectively, for Compound A. The reference temperature selected for the time-temperature superposition is 140 °C, which is around the Tg of the compound. The temperature shift factors as shown in Fig. 6 are defined such that Z ξ¼

0

t

aT ðT Þdτ;

ð23Þ

where ξ is the reduced time and aT is the temperature shift factor. Under constant temperature conditions during the DMA experiments, the integral form of Eq. (15) can be reduced to ξ = aT(T)t, and similarly, the reduced frequency for the storage modulus mastercurves (Fig. 5) can be expressed as ϖ¼

ω : aT ðT Þ

ð24Þ

The storage modulus mastercurves and the corresponding temperature shift factors for Compound B are shown in Figs. 7 and 8, respectively. The reference temperature used for the superposition is 120 °C, which is around the Tg of Compound B. Mathematical representation of the shift factors shown in Figs. 6 and 8 is expressed by using the Williams, Landel, and Ferry (WLF) equation, given by log10 aT ¼

C 1 ðT−T re f Þ C 2 þ ðT−T re f Þ

ð25Þ

where Tref is the reference temperature, C1 and C2 are model constants. The fitted WLF models are also shown in Figs. 6 and 8, and the models constants for Compounds A and B are given in Table 1. 3.2. Multiaxial viscoelastic model development Conversions of the frequency-domain experimental storage modulus mastercurves to the time-domain material functions to be used in the multiaxial constitutive Eq. (1) were first performed by the procedures described in Section 2.1. Note that the interconversions were performed on the mastercurves and (t, ω) are replaced by (ξ, ϖ). Using Eqs. (8) and (10), the experimentally measured storage moduli were first converted to the relaxation moduli in the Laplace–Carson transform domain. Substituting the Laplace–Carson-domain relaxation moduli into Eqs. (11) and (13), the bulk relaxation moduli of the EMCs were obtained and are shown in Fig. 9. It can be seen from Fig. 9 that the calculated bulk moduli exhibit the proper relaxation over time behavior, but display some fluctuations in both small- and large-time regimes. The fluctuating trend is more significant in the Laplace–Carson-domain viscoelastic Poisson's ratios, which were calculated from Eq. (12) and are shown in Fig. 10. It also can be seen from Fig. 10 that, for small s values (corresponding to large time), the Poisson's ratios oscillate and become negative, which is physically inadmissible for the particulate composites considered in this study. The peculiar trend observed in Fig. 10 may be attributed to that the viscoelastic responses (Young's and shear storage moduli) were not measured simultaneously on the same specimen, as cautioned by Lu et al. [12]. To further investigate into this issue, multiple DMA runs using different specimens having the same geometry were performed to check repeatability of the experimental results. The results show some variations in the storage modulus values, mainly at temperatures slightly higher than Tg. These differences were not significant considering the overall response mastercurve, but could lead to significant changes in the viscoelastic Poisson's ratio, as can be seen from Eq. (12). To overcome this difficulty, the numerical procedure described in Section 2.2 was applied for constructing the multiaxial constitutive

Table 2 Prony series constants for the bulk modulus of Compound A, M = 17. i

wκi

τκi (s)

i

wκi

1 2 3 4 5 6 7 8 9

0.0255 0.0380 0.0455 0.0500 0.0684 0.0260 0.0812 0.0734 0.1372

1.0E−11 1.0E−08 1.0E−05 1.0E−03 1.0E−01 1.0E+00 1.0E+01 1.0E+02 1.0E+03

10 0.1432 11 0.1207 12 0.0756 13 0.0235 14 0.0103 15 0.0103 16 0.0103 17 0.0011 κg = 11.1 GPa

τκi (s) 1.0E+04 1.0E+05 1.0E+06 1.0E+07 1.0E+08 1.0E+09 1.0E+11 1.0E+13

model. Shown in Figs. 11–13 are the fitted viscoelastic material functions of the Compound A, with the corresponding Prony series constants given in Tables 2 and 3. It can be seen from Fig. 11 that the fitted storage modulus mastercurves agree very well to the experimental data in the glassy regime, but are slightly different in the glass transition and rubbery regimes. However, the slightly compromised fitting to the storage moduli gives much better time-domain relaxation moduli (Fig. 12) and viscoelastic Poisson's ratio (Fig. 13) than the results of the direct conversion. The fitted viscoelastic functions for Compound B are shown in Figs. 14–16, with the corresponding Prony series constants given in Tables 4 and 5. Similar to the results for Compound A, the fitted models for Compound B agree well to the experimental data shown Fig. 14, and the fitted time-domain materials functions shown in Figs. 15 and 16 are more realistic. From the comparison of the viscoelastic Poisson's ratios shown in Figs. 13 and 16, it can be seen that the filler percentage has a strong influence on the Poisson's ratio. The characteristic of the Compound A (88% filler) is dominated by the behavior of the silica, which has a low Poisson's ratio, and the compound's Poisson's ratio stays at relatively low values over time. On the other hand, the epoxy resin in Compound B (70% filler) appears to have a stronger influence on the overall compound behavior, and the Poisson's ratio value is much higher compared to the one of Compound A. 3.3. Package warpage simulation To evaluate the multiaxial viscoelastic model, the bulk and shear relaxation moduli of Compound A were implemented in the thermomechanical model of an overmolded ball grid array (BGA) package. The numerical model was developed by using finite element software ANSYS [14] and applied to simulate warpage response of the package under solder joint reflow temperature range (between 30 and 260 °C), and compared to experimental shadow Moiré measurements. The BGA package considered has an overall dimension of 20.2 mm × 20.2 mm × 1.06 mm, and contains a Si die of 5.3 mm × 4.7 mm × 0.2 mm. Thicknesses of the laminate substrate

Table 3 Prony series constants for the shear modulus of Compound A, N = 17. i

wμi

τμi (s)

i

wμi

1 2 3 4 5 6 7 8 9

0.0291 0.0429 0.0451 0.0498 0.0675 0.0264 0.0794 0.0739 0.1348

1.0E−11 1.0E−08 1.0E−05 1.0E−03 1.0E−01 1.0E+00 1.0E+01 1.0E+02 1.0E+03

10 0.1431 11 0.1185 12 0.0759 13 0.0224 14 0.0118 15 0.0204 16 0.0127 17 0.0027 μg = 9.60 GPa

τμi (s) 1.0E+04 1.0E+05 1.0E+06 1.0E+07 1.0E+08 1.0E+09 1.0E+11 1.0E+13

Please cite this article as: T.-C. Chiu, et al., Time-domain viscoelastic constitutive model based on concurrent fitting of frequency-domain characteristics, Microelectronics Reliability (2015), http://dx.doi.org/10.1016/j.microrel.2015.07.031

T.-C. Chiu et al. / Microelectronics Reliability xxx (2015) xxx–xxx

16

storage modulus (GPa)

14 12

Table 4 Prony series constants for the bulk modulus of Compound B, M = 17.

Experiment Fitted

10

7

E'

8 6

'

4

i

wκi

τκi (s)

i

1 2 3 4 5 6 7 8 9

0.0164 0.0299 0.0209 0.0164 0.025 0.0213 0.0208 0.0182 0.0298

1.0E−21 1.0E−19 1.0E−17 1.0E−15 1.0E−13 1.0E−11 1.0E−09 1.0E−07 1.0E−05

10 0.0188 11 0.1147 12 0.1699 13 0.2836 14 0.124 15 0.0295 16 0.0141 17 0.0087 κg = 14.4 GPa

2 0 10-20

10-10

100

1010

1020

1030

(1/s) Fig. 14. The fitted storage modulus mastercurve of Compound B, reference temperature 120 °C.

relaxation modulus (GPa)

16 14 12

E

10

wκi

τκi (s) 1.0E−03 1.0E−01 1.0E+01 1.0E+03 1.0E+05 1.0E+07 1.0E+09 1.0E+11

Table 5 Prony series constants for the shear modulus of Compound B, N = 20. i

wμi

1 0.0107 2 0.0213 3 0.0293 4 0.0206 5 0.0162 6 0.025 7 0.0209 8 0.0207 9 0.0178 10 0.0301 μg = 4.75 GPa

τμi (s)

i

wμi

τμi (s)

1.0E−23 1.0E−21 1.0E−19 1.0E−17 1.0E−15 1.0E−13 1.0E−11 1.0E−09 1.0E−07 1.0E−05

11 12 13 14 15 16 17 18 19 20

0.03779 0.02332 0.02072 0.01981 0.01986 0.02892 0.02698 0.01614 0.02442 0.01664

0.0185 0.1132 0.168 0.2813 0.1226 0.0295 0.013 0.0108 0.0019 0.0104

8 6 4 2 0 10-30

10-20

10-10

100

1010

1020

(s) Fig. 15. The fitted relaxation modulus mastercurve of Compound B, reference temperature 120 °C.

and mold compound above die are 0.36 mm and 0.475 mm, respectively. Shown in Fig. 17 is the quarter finite element model used for the warpage simulation. The model was constructed by using quadratic solid elements with model boundary conditions such that zero out-of-plane displacements were assigned on the package

0.5 0.4

symmetric planes, and the bottom center point was fixed to prevent free-body motion. Aside from the EMC, the constitutive behavior of the package constituents including the Si die, die attach adhesive and the laminate substrate were modeled as linear elastic with temperature-dependent elastic moduli and coefficients of thermal expansion (CTEs). Material constants of the linear thermomechanical models are given in Table 6. Shown in Fig. 18 are the warpage simulation results compared to experimental shadow Moiré measurements. The experimental warpage values were collected on 15 package samples by using the AkroMetrix TherMoiré PS400 instrument. Note that the sign of the warpage value is such that, viewing from the substrate side, the concave and convex deformations correspond to positive and negative warpages, respectively. Also included in Fig. 18 are the simulation results based on a simple temperature dependent linear elastic constitutive model for the Compound A. It can be seen from Fig. 18 that the numerical results based on the viscoelastic model agree well to the shadow Moiré measurements. It validates the multiaxial EMC constitutive model developed in this study. When the viscoelastic EMC model in the finite element model is replaced with an elastic one, the simulation gives a different warpage evolution trend. In the elastic model, the error in the relative warpage trend between 130 °C and 180 °C can be attributed to the significant change in the thermomechanical properties of the EMC near its Tg. As temperature cools down from above to below the

0.3 0.2 0.1 0.0 -30 10

10-20

10-10

100

1010

1020

(s) Fig. 16. The fitted viscoelastic Poisson's ratio mastercurve of Compound B, reference temperature 120 °C.

Fig. 17. Quarter finite element model for the overmolded BGA package, 8064 elements with 34,513 nodes.

Please cite this article as: T.-C. Chiu, et al., Time-domain viscoelastic constitutive model based on concurrent fitting of frequency-domain characteristics, Microelectronics Reliability (2015), http://dx.doi.org/10.1016/j.microrel.2015.07.031

8

T.-C. Chiu et al. / Microelectronics Reliability xxx (2015) xxx–xxx

200

Table 6 Linear thermomechanical constants of the BGA package.

shadow Moire (t), constant E(t), constant E(t) and (t) E(t), constant (t), constant

Si die Die attach adhesive

Laminate substrate

EMC

Young's modulus (GPa)

Poisson's ratio

CTE (ppm/°C)

150 0.821 (11 °C) 0.499 (21 °C) 0.313 (31 °C) 0.155 (50 °C) 0.06 (100 °C) 0.048 (146 °C) 23.7 (−60 °C) 19.3 (90 °C) 14.9 (170 °C) 12.4 (250 °C) 10.1 (290 °C) 22.5 (40 °C) 21 (110 °C) 19.4 (130 °C) 17.3 (140 °C) 3.25 (170 °C) 1.43 (190 °C) 1.12 (220 °C)

0.17 0.4

2.9 95 (−40 °C) 95 (17 °C) 220 (27 °C) 220 (260 °C)

0.19

warpage ( m)

150

17.5 (−40 °C) 17.5 (50 °C) 17 (200 °C) 17 (260 °C)

0.17

7.8 (−80 °C) 7.8 (45 °C) 9.77 (105 °C) 13.7 (130 °C) 29.9 (155 °C) 37.2 (260 °C)

EMC glass transition temperature, the elastic EMC model cannot account for the time-dependent material relaxation, and, as a result, would overestimate the thermal stress and warpage change. An often adopted engineering practice when only one independent viscoelastic function is available through experiments is to assume that either the bulk modulus or the Poisson's ratio is time- and temperature-independent. For example, the multiaxial viscoelastic model based on tensile DMA characterization and an assumed constant bulk modulus can be described by (μ(t), κ = κg), in which the shear relaxation modulus can be obtained by first converting the Prony series of ~ κ~ −EÞ ~ (in E′(ω) to the Laplace–Carson domain, calculating μ~ ¼ 3κ~ E=ð9 ~ which κ ¼ κ g ), and then fitting for the Prony series constants of μ(t) in the Laplace–Carson domain. To investigate the effectiveness of these assumptions, multiaxial viscoelastic models of Compound A were constructed from the viscoelastic Young's or shear modulus mastercurve (Fig. 12) with the assumption of either the bulk modulus or Poisson's ratio remains constant at its glassy value over time, and then implemented in the finite element model to simulate package warpage. Comparisons of the numerical results based on these assumptions and the experimental data are shown in Fig. 19, from which it can be seen that the warpage predictions based on the constant Poisson's ratio are very close to the one based on the comprehensive viscoelastic model; while the prediction based on constant bulk modulus deviates slightly more. The less accuracy for the predictions based on constant bulk modulus assumption can simply attributed to that the viscoelastic behavior of Compound A (Figs. 12 and 13) is closer to the constant-ν assumption

0 -50

-150

0

50

100

shadow Moire viscoelastic model elastic model

-50 -100

250

300

shadow Moire E(t) and (t) E(t), constant (t), constant E(t), constant (t), constant

150

0

200

than to the constant-κ assumption. Nevertheless, the warpage trends predicted with the constant bulk modulus assumption are still reasonable good if only one of the Young's and shear storage modulus mastercurves is available for constructing the multiaxial constitutive model. Another typically made assumption in modeling the viscoelastic behavior of EMCs is to assume a glassy Poisson's ratio in the range of 0.3 and 0.4, which are true for many engineering materials. By assuming a glassy Poisson's ratio of 0.35, multiaxial constitutive models of Compound A based on the viscoelastic Young's or shear modulus with the assumption of either ν or κ remains constant at its glassy value were also developed and implemented in finite element warpage simulation. Results of the warpage simulations based on these models are shown in Fig. 20. By comparing Figs. 19 and 20, it can be seen that the warpage trends predicted with the assumed Poisson's ratio deviate a lot more from the experiment data. In addition, from Fig. 20 it can be seen that the simulation predictions based on multiaxial models constructed by using only viscoelastic Young's modulus mastercurve are more accurate compared to models based on viscoelastic shear modulus mastercurve. The better warpage prediction accuracy for the Young's-modulus-mastercurve based multiaxial model may be attributed to that the package warpage during thermal excursion is mainly associated to the normal strains and stresses in the package in-plane directions, and the multiaxial model based on viscoelastic Young's modulus inherently gives a better prediction of this deformation response than the viscoelastic-shear-modulus based model would provide.

warpage ( m)

50

150 (°C)

Fig. 19. Effects of viscoelastic model assumptions on warpage simulation results.

200

100

warpage ( m)

50

-100

150

100 50 0 -50 -100

-150 -200

100

assuming (0) = 0.35 -150 0

50

100

150 (°C)

200

250

300

Fig. 18. Comparisons of warpage simulation predictions and experimental data.

0

50

100

150 (°C)

200

250

300

Fig. 20. Effect of viscoelastic model with assumed glassy Poisson's ratio on warpage simulation results.

Please cite this article as: T.-C. Chiu, et al., Time-domain viscoelastic constitutive model based on concurrent fitting of frequency-domain characteristics, Microelectronics Reliability (2015), http://dx.doi.org/10.1016/j.microrel.2015.07.031

T.-C. Chiu et al. / Microelectronics Reliability xxx (2015) xxx–xxx

4. Conclusions A numerical procedure was developed for constructing the timedomain multiaxial linear viscoelastic constitutive model from separately measured frequency-domain Young's and shear storage moduli. Possible physics incompatibility in the directly converted time-domain material functions is prevented with the incorporation of optimization procedure to best fit the separately measured material functions under the restriction of Poisson's ratio being a monotonic increasing function of time. Multiaxial viscoelastic model constructed based on the proposed procedure was implemented in warpage simulation of an overmolded package and validated with shadow Moiré measurements. From additional warpage simulations based on only either the Young's or shear storage modulus and the assumption of either constant bulk modulus or constant Poisson's ratio, it was found that the viscoelastic models based only on the tension DMA characterization predict package warpage deformation better than the one based only on the torsion DMA characterization. Acknowledgments This study was sponsored by the National Science Council of ROC under the grant NSC 101-2221-E-006-023-MY2 and by ASE Group. The authors would also like to thank the Group R&D of ASE for equipment and material support.

9

[2] T.C. Chiu, J.L. Gung, H.W. Huang, Y.S. Lai, Effects of curing and chemical aging on warpage — characterization and simulation, IEEE Trans. Device Mater. Reliab. 11 (2011) 339–348. [3] L.J. Ernst, C. van't Hof, D.G. Yang, M.S. Kiasat, G.Q. Zhang, H.J.L. Bressers, J.F.J. Caers, A.W.J. den Boer, J. Janssen, Determination of visco-elastic properties during the curing process of underfill materials, Proceedings of the 50th Electronic Components and Technology Conference 2000, pp. 1070–1077. [4] M. Sadeghinia, K.M.B. Jansen, L.J. Ernst, Characterization of the viscoelastic properties of an epoxy molding compound during cure, Microelectron. Reliab. 52 (2012) 1711–1718. [5] R.A. Schapery, S.W. Park, Methods of interconversion between linear viscoelastic material functions. Part II—an approximate analytical method, Int. J. Solids Struct. 36 (1999) 1677–1699. [6] I. Emri, B.S. von Bemstorff, R. Cvelbar, A. Nikonov, Re-examination of the approximate methods for interconversion between frequency- and time-dependent material functions, J. Non-Newtonian Fluid Mech. 129 (2005) 75–84. [7] S.W. Park, R.A. Schapery, Methods of interconversion between linear viscoelastic material functions. Part I—a numerical method based on Prony series, Int. J. Solids Struct. 36 (1999) 1653–1675. [8] J. Luk-Cyr, T. Crochon, C. Li, M. Lévesque, Interconversion of linearly viscoelastic material functions expressed as Prony series: a closure, Mech. Time Depend. Mater. 17 (2013) 53–82. [9] N.W. Tschoegl, W.G. Knauss, I. Emri, Poisson's ratio in linear viscoelasticity—a critical review, Mech. Time Depend. Mater. 6 (2002) 3–51. [10] N.W. Tschoegl, The Phenomenological Theory of Linear Viscoelastic Behavior, Springer-Verlag, Berlin, 1989. [11] M. Jerabek, D. Tscharnuter, Z. Major, K. Ravi-Chandar, R.W. Lang, Relaxation behavior of neat and particulate filled polypropylene in uniaxial and multiaxial compression, Mech. Time Depend. Mater. 14 (2010) 47–68. [12] H. Lu, X. Zhang, W.G. Knauss, Uniaxial, shear and Poisson relaxation and their conversion to bulk relaxation, Polym. Eng. Sci. 37 (1997) 1053–1064. [13] MATLAB Release 2012b, The MathWorks, 2012. [14] ANSYS Version 15, Ansys, 2013.

References [1] W. Lin, M.W. Lee, PoP/CSP warpage evaluation and viscoelastic modeling, Proceedings of the 58th Electronic Components and Technology Conference 2008, pp. 1576–1581.

Please cite this article as: T.-C. Chiu, et al., Time-domain viscoelastic constitutive model based on concurrent fitting of frequency-domain characteristics, Microelectronics Reliability (2015), http://dx.doi.org/10.1016/j.microrel.2015.07.031