- Email: [email protected]
Simplified approach for steady thermal analysis of chips with variable power based on spatial autocorrelation
Applied Thermal Engineering 120 (2017) 347–357
Contents lists available at ScienceDirect
Applied Thermal Engineering journal homepage: www.elsevier.com/locate/apthermeng
Research Paper
Simplified approach for steady thermal analysis of chips with variable power based on spatial autocorrelation C. Pei ⇑, G.C. Fu, Y.H. Zhao, B. Wan, Y. Cheng School of Reliability and Systems Engineering, Beihang University, Beijing 100191, China
h i g h l i g h t s Approach for steady thermal analysis of chip with variable power is proposed. The determination of material properties and other factors can be tolerated. The highest temperature error in the case study is not more than 2%.
a r t i c l e
i n f o
Article history: Received 20 December 2016 Accepted 4 April 2017 Available online 5 April 2017 Keywords: Steady thermal analysis Spatial autocorrelation Undetermined coefficient method Power devices
a b s t r a c t Due to their growing power and shrinking size, thermal management and analysis have become one of the most critical concerns for power devices. Though finite element (FE) method is commonly employed in thermal analysis now, the uncertainty of material properties and the redundancy of repeat modeling still nag analyzers. In this paper, a simplified approach for steady thermal analysis of chips with variable power is proposed. Based on spatial autocorrelation and undetermined method, the steady temperature distribution could be directly determined with power distribution and example temperature measurement result, without identifying environment factors or material properties. The theory of the simplified approach is introduced at first in this paper. The characteristics of the simplified approach are discussed afterwards. Eventually, a case study is given to illustrate the procedure, of which the analytical result coincides with the actual measurement well. Ó 2017 Elsevier Ltd. All rights reserved.
1. Introduction As to the continual miniaturization and increasing power dissipation, thermal management has become a critical factor in the design process of microcircuits. Especially for high power microcircuits, i.e. insulated gate bipolar transistor (IGBT), high temperature, induced by heat flow, will cause parameter drifting or even failures. Numerous researches have been conducted on die-attach layer or its inner defects. Within the researches, the size [1,2], position [3], pattern [4], and even depth [5] of the die-attach defects have been discussed, and their influences on thermal [6], mechanical [7,8] and parameter-stable [9,10] performance have been also studied. In most of those researches, heat sources have been treated as an entity. However, the heat generation distribution affects thermal analysis results at the same time. Hot spots may be concealed by the assumption of uniform heating. For devices composed of plentiful repetition units, i.e. IGBTs or memories, the ⇑ Corresponding author. E-mail address: [email protected] (C. Pei). http://dx.doi.org/10.1016/j.applthermaleng.2017.04.010 1359-4311/Ó 2017 Elsevier Ltd. All rights reserved.
distribution of heat sources is also determined by the operating mode [11,12]. Meanwhile, current thermal analyses of devices are mainly based on equivalent thermal resistance network or finite element (FE) method [13–15]. No matter which method is employed, dissipation factors, including material property, dimension parameter, power distribution and boundary condition, should be determined in advance. The precision of the thermal analysis depends on the accuracy of the factors set in the analysis models to a great extent, which may be difficult to acquire sometimes. For example, the heat conductivity of the die-attach depends on both the material and the process, which will vary with each individual, especially in the case of porous die-attach pastes [16,17]. At the same time, whenever the operating mode switches, the analysis models have to be modified or rebuilt correspondingly, which will also increase the modeling work. To cope with the problems, an approach that could simplify the thermal analysis of chips with variable power will be proposed in this paper. With this approach, the steady temperature distribution could be acquired only with the information of power distribution
348
C. Pei et al. / Applied Thermal Engineering 120 (2017) 347–357
in different operating modes and the example temperature measurement result. The paper is arranged as below. Firstly, the theory of the simplified approach will be introduced. Subsequently, a series of typical power distributions will be given, and their thermal influences will be solved by the FE method and the proposed approach. Their temperature distributions will be compared to verify the accuracy of the approach. Meanwhile, the relationship between the undetermined coefficients and other dissipation factors will be discussed. Finally, a case study will be given. The steady temperature distribution measured by infrared microscopy and the analytical result based on this approach will be compared again to verify the approach. 2. Theory & methodology The typical heat dissipation route for encapsulated chips is illustrated in Fig. 1(a). As to the low heat conductivity of thermal molding compound, most of heat dissipates from the die-attach layer through the carrier to the mounting plane or the heat sink [18]. In a simplified form, the temperature of chip is determined by power distribution, power intensity, material thermal properties and boundary conditions, as Fig. 1(a) shows. In this paper, spatial analysis is utilized to depict the power distribution and identify the influence, which has been employed in other works on defect detection [19,20] or condition evaluation [21]. Among spatial analysis methods, spatial autocorrelation method, proposed by A. D. Cliff and J. K. Ord [22] in 1980s, is most promising and has been applied to wafer defect detection [23]. In these studies, the defective dies, fabricated on a single wafer, were summed as a wafer map to display the spatial dependence across the wafer. Based on it, spatial autocorrelation was extracted and spatial effect could be resolved. Similar to wafer map, chips with variable power will be also meshed into matrix-arranged lattices in this paper, as Fig. 2 illustrates. Let matrix Ma,b denote the chip. The power of unit (i, j) is expressed as X (i, j). The definition of the power of each unit is:
Xði; jÞ ¼
1; work 0; non work
ð1Þ
For simplification, the power of each unit is assumed to be equal and Xa,b is transformed into a 0–1 matrix. The distance was measured by counting the minimum steps between two units. In other words, the heat is only permitted to transfer between adjacent units. For any unit (h, i) and any other unit (j, k), their distance could be expressed as:
dh;i ðj; kÞ ¼ jh jj þ ji kj
ð2Þ
For unit (h, i), the units at d distance are separately summed as:
W 0h;i ðgÞ ¼ fXðj; kÞjdh;i ðj; kÞ ¼ g; Xðj; kÞ ¼ 0; g P 1g
ð3Þ
W 1h;i ðgÞ ¼ fXðj; kÞjdh;i ðj; kÞ ¼ g; Xðj; kÞ ¼ 1; g P 1g
ð4Þ
H1h;i ðgÞ ¼ jW 1h;i ðgÞj
ð5Þ
|| in Eq. (5) denotes the total number of set elements. Particularly, Hh,i(0) denotes the operating condition of unit (h, i) itself. According to Fourier’s law, the temperature rise of units is dependent on the thermal resistance, heat generation, and boundary heat flux. For each unit in the chip, excepting the heat generated interiorly and the heat flow exteriorly the chip bottom, the heat inflow (or outflow) from neighboring units should also be considered. Actually, the hot spots appear not only in the area with high heat intensity but also in the operation units ‘‘crowded” area. According to the first law of geography, ‘‘everything is related to everything else, but near things are more related than distant things” [24], the synergy influences between units should be considered and converge with distance. The simplified equivalent circuit of the dissipation route is illustrated in Fig. 1(b). According to superposition theorem, if temperature variation is limited and the thermal material properties remain stable, the steady temperature rise of the chip could be expressed as the linear superposition of influences of all heat sources (or operating units). Based on the hypothesis, the steady temperature matrix Ta,b could be expressed as below:
T ¼ AZ þ B
Fig. 1. Dissipation route and equivalent circuit: (a) Typical dissipation route for encapsulated chips (b) the simplified equivalent circuit.
Fig. 2. Chip meshed into matrix-arranged lattices.
ð6Þ
349
C. Pei et al. / Applied Thermal Engineering 120 (2017) 347–357
1 H11;1 ða þ b 1Þ H1;1 ð0Þ H11;1 ð1Þ C B B H2;1 ð0Þ H12;1 ð1Þ H11;2 ða þ b 1Þ C C B Z¼B C .. .. .. C B A @ . . . H1a;b ða þ b 1Þ ab;aþb1 Ha;b ð0Þ H1a;b ð1Þ 0
0 B B A¼B B @
Að0Þ Að1Þ .. .
3. Modeling
ð7Þ
1 C C C C A
ð8Þ
Aða þ b 1Þ Within the undetermined coefficients, A denotes the spatial correlation and the heat intensity. B denotes the environment factors, which is assumed to include boundary conditions related factor, i.e. boundary temperature and boundary heat exchange. As the power of each unit is normalized, not only matrix Z but also coefficient A affects the total power. The former represents the total power area, and the latter represents the heat intensity. Vector T denotes the temperature distribution. Obviously, with the power distribution of the chip and the example temperature measurement result, undetermined coefficient A and B in Eq. (6) could be solved. And if A and B known, the steady temperature distribution of the chip under any power distribution could be easily determined, during which the repetitive modeling and determination of material properties could be avoided. The flowchart of the proposed simplified approach is displayed in Fig. 3.
3.1. Object Before applying the proposed approach, the accuracy of this approach should to be confirmed, as well as its other characteristics. Several typical distributions are chosen here, based on which a series of numerical models are established to acquire reference value for the confirmation and discussions. To eliminate the effects of the package, a model only consisting of chip, carrier and heat sink is established, as Fig. 4 illustrates. For simplification, the dieattach layer supposed to be homogenous. The dimensions of the layers are listed in Table 1. The thermal properties are also listed in the table, which are acquired from [18]. The environment temperature of the thermal analysis is set to 25 °C, as well as the air. Only natural convection heat dissipation is considered in the model and the direction of gravity is set from the chip to the heatsink. Radiation dissipation is ignored in the model, as heat transfer by radiation is negligible under a low temperature like this case.
Table 1 Dimension and material of example chip. Layer
Material
Heat conductivity (W/m k)
Density (g/m3)
Chip Carrier Heat sink
GaAs Duralumin Duralumin
59.1 180 180
5310 2700 2700
Fig. 3. Flowchart of the Simplified approach.
Fig. 4. Structure diagram of numerical model.
350
C. Pei et al. / Applied Thermal Engineering 120 (2017) 347–357
3.2. Power distribution The chip is divided into 400 (20 ⁄ 20) lattices symbolically in its horizontal plane. As mentioned before, the power of each unit is supposed to be equal, which is set to 0.01 W in this part. The next step is to determine several typical distributions. In [23], four typical distributions are chosen. While the circle pattern is most common in wafer process, concentration is the most fatal case in chip dissipation. To take the pertinence in account, the circle pattern is substituted by center and edge concentration. Additionally, edge concentrated power distribution is also selected to investigate the edge effect, though it seldom appears in chips effect. Eventually, five distributions, including spatially homogeneous Bernoulli process, edge concentration, center concentration, cluster, and repetition, are established in the paper, which are displayed in Table 2. 3.3. Numerical settings Commercial thermal analysis software FloTHERM is employed in this paper. Finite volume method (FVM) is employed. Localized grids are toggled to reduce the amount of elements and assure the quality of the meshing. The elements of all the meshes are about 456,000. Multi grid method is utilized in solving to accelerate the
convergence. By default, if the residual error within the calculation reached 0.5% of total source value, the solution is deemed converged, as well as the velocity and pressure. 4. Results and discussions 4.1. Characterization of coefficients All the numerical models converge within 200 iterations. Their steady temperature distributions are shown in Fig. 5. As could be seen, the temperature distributions vary with the different power distributions. The highest temperature appears in the edge concentration distribution, which reveals that the edge effect is the most significant factor affecting the temperature rise. Besides, among the other four distributions, the highest and lowest temperatures both appear in the center concentration distribution. The uneven temperature distribution may be induced by power concentration. In other words, the concentration effect is another factor that should be considered. As mentioned in Section 2, the key procedure of the proposed approach is to determine underdetermined coefficients A and B in Eq. (6). Before utilizing the coefficients to predict the temperature distributions under different power distributions, their stabil-
Table 2 Typical power distributions. Percentage of operating unit
SHBP
Edge concentration
Center concentration
Cluster
Repetition
No. 30%
I
II
III
IV
V
No.
VI
VII
VIII
IX
X
20%
Fig. 5. Temperature distribution of by simulation: (a) distribution I, (b) distribution II, (c) distribution III, (d) distribution IV, (e) distribution V, (f) distribution VI, (g) distribution VII, (h) distribution VIII, (i) distribution IX, (j) distribution X.
351
C. Pei et al. / Applied Thermal Engineering 120 (2017) 347–357
ity should be validated beforehand. To achieve this goal, the spatial correlogram matrix of all the typical distributions are extracted. Then, A and B are solved with the thermal simulation analysis results by the least square method. Considering that the spatial correlogram matrix Z may not be full rank, underdetermined coefficient A is suggested to be truncated by a specific distance. Before taking this strategy, the influence of the truncation should also be analyzed in advance. As mentioned above, underdetermined coefficient A, correlation influences between units, is assumed to converge with distance. Moreover, it is assumed to keep stable under different truncations. Next, the assumptions will be proved. Taking distribution VII as an example, the fitting results under the different truncations are shown in Table 3. As could be seen, the coefficient A decreases with the distance. Then the relationship between coefficient A under 10order truncation and the distance is plotted in Fig. 6(a). An exponential function is employed, whose form is:
f ðgÞ ¼ ag b þ c
ð9Þ
The fitting curve is also included in Fig. 6(a). The goodness demonstrates that the data are well fitted to an exponential function. Besides, as symbol b denotes convergence rate, the convergence rate of A(g) is considered to be 2.25. Considering that Z(g) is proportional to distance g, their combined influence is considered to converge and the rate is 1.25. In other words, the influences of the coefficient A decrease with the distance, which means that the error induced by truncation could be ignored. Then the stability of the coefficients against the different truncation orders is investigated. Distribution VII is taken as an
example again, and the trend of A(0) values under the different truncation orders is plotted in Fig. 6(b). It could be observe that A(0) fluctuates at first. Then it stays stable briefly until the distance reaches the rank of matrix Z, which verifies the stability and uniqueness of the coefficients. Since the reasonability of the truncation has been proved, coefficient A and B in all the distributions by 4-order truncated are selected to compare, as shown in Table 4. The comparison reveals that coefficient B is proportional to the total power and A is immune to the power variation. Besides, the result demonstrates that nearly all the coefficients are independent from the power distribution. 4.2. Accuracy and improvement The average value of the coefficients, except distribution III and IV’s is chosen as the general value. Then the general value are substituted into Eq. (6) to predict the temperature. The residual error between the analysis and the simulation is displayed in Fig. 7. It could be found that the distribution of the residual error is independent from the distribution of the operating units and concentrates on the edge. The residual errors in distribution III and IV are mostly significant. This suggests that the dissipation of the edge units is much more complex than that of interior units. As the heat transfer coefficient of the chip edge is far less than the die-attach, the edge could be treated as an adiabatic boundary. Contrary to the conventional treatment in FEA, the chip is expanded to eliminate the edge influence, as Fig. 8(a) illustrates. The space of the chip is divided into 4 quadrants. Each quadrant
Table 3 Fitted coefficient A in distribution VII under different truncation orders. Truncation order
A(0)
A(1)
A(2)
A(3)
A(4)
A(5)
A(6)
A(7)
A(8)
A(9)
1 2 3 4 5 6 7 8 9 10
0.354 0.219 0.247 0.249 0.248 0.249 0.247 0.249 0.249 0.249
0.084 0.037 0.046 0.047 0.045 0.047 0.047 0.048 0.048
0.041 0.017 0.021 0.024 0.024 0.024 0.023 0.023
0.025 0.010 0.013 0.014 0.013 0.013 0.013
0.020 0.008 0.008 0.009 0.009 0.009
0.014 0.007 0.008 0.007 0.007
0.009 0.004 0.005 0.005
0.005 0.003 0.004
0.003 0.002
0.003
Fig. 6. Convergence of coefficient A: (a) coefficient A convergence with distance; (b) convergence of A(0) under different truncation orders.
352
C. Pei et al. / Applied Thermal Engineering 120 (2017) 347–357
Table 4 Comparison of coefficients with different distributions and total powers. Distribution No.
Total power (W)
B (°C)
Temperature rise (°C)
A(0)
A(1)
A(2)
A(3)
I II III IV V VI VII VIII IX X
0.8 1.2 0.8 1.2 0.8 1.2 0.8 1.2 0.8 1.2
40.634 46.784 40.971 46.608 40.509 46.567 40.568 47.109 40.548 46.661
15.634 21.784 15.971 21.608 15.509 21.567 15.568 22.109 15.548 21.661
0.257 0.258 0.765 0.447 0.250 0.257 0.248 0.245 0.246 0.247
0.048 0.047 0.338 0.102 0.049 0.053 0.047 0.051 0.048 0.047
0.020 0.018 0.216 0.042 0.025 0.027 0.024 0.025 0.023 0.025
0.009 0.008 0.073 0.026 0.017 0.015 0.013 0.015 0.014 0.015
Fig. 7. Residual error before expansion: (a) distribution I, (b) distribution II, (c) distribution III, (d) distribution IV, (e) distribution V, (f) distribution VI, (g) distribution VII, (h) distribution VIII, (i) distribution IX, (j) distribution X.
Fig. 8. Mirror expansion (a) schematic diagram (symbol ‘F’ represents the original distribution); (b) mirror expansion of distribution I.
is mirrored upside down, left and right. In this way, the original edge could be assured to be adiabatic. Distribution I is taken as an example. Its expanded distribution is shown in Fig. 8(b). The area outlined in yellow (top-left) is the original distribution and the area outlined in red (center) is the expanded distribution from the second quadrant. All the typical distributions are expanded as above. Then, the coefficients are solved again, and the convergence of coefficient A is shown in Fig. 9. New general coefficients are employed to
determine the residual errors by the mirror expansion. The results are illustrated in Fig. 10. By the comparison, it could be found that the residual errors decline from 0.2 to 0.1. Besides, the anomalous units located in the edge are also reduced. Nevertheless, though the anomaly caused by edge influence is relieved slightly, coefficient A in distribution III and IV are still larger than those of other distributions. Considering that power units are not likely to concentrate in the edge, the error induced by the edge influence may still be acceptable.
353
C. Pei et al. / Applied Thermal Engineering 120 (2017) 347–357
Firstly, distributions I, V, VII and IX are chosen. Higher environment temperature and fluid temperature, 35 °C, are set to analyze the influences. Mirror expansion method is employed in the calculation and the fitting results are displayed in Table 5. With the results, it could be observed that only coefficient B increases and the increment is close to the rise of the environment temperature. Oppositely, coefficient A is insensitive to the change of environment temperature.
Fig. 9. Convergence of new coefficient A in expanded distributions.
4.3. Relationship between coefficients and influencing factors 4.3.1. Environment influence In Section 4.1, the influences of heat sources (total power) on the coefficients have been investigated. Here, the relationship between the material properties or boundary conditions and the coefficients will be discussed further.
4.3.2. Material influence The influences of heat conductivity are also studied below. The heat conductivity of the chip and the carrier is doubled separately, which is depicted in Table 1. The new coefficient fitting results are displayed in Table 6. It could be observed that the increase of chip heat conductivity only affects coefficient A(0). By contrast, the variation of carrier heat conductivity will exert an influence on all the coefficients. Coefficient A nearly halves, except A(0). These differences could be explained by the different order of magnitudes of the resistances. The correlation resistance between units Rco in the Fig. 1(b), is about 600 K/W, while the exterior resistance Rout (Rtrans and Rbound) is about 160 K/W, of which nearly 90% is the inner resistance of carrier. Therefore, all of coefficient A are influenced by the heat conductivity of the carrier, as the influences of Rout are more global. By contrast, The variation of the ratio between A(0) and other coefficient A demonstrates that Rco only affect the heat flow between power units.
Fig. 10. Residual error after expansion: (a) distribution I, (b) distribution II, (c) distribution III, (d) distribution IV, (e) distribution V, (f) distribution VI, (g) distribution VII, (h) distribution VIII, (i) distribution IX, (j) distribution X.
Table 5 Comparison of coefficients under different environment temperature. Distribution No.
Coefficient
B
A(0)
A(1)
A(2)
A(3)
A(4)
A(5)
I
35 °C 25 °C Increment
50.615 40.499 10.120
0.250 0.247 0.003
0.048 0.048 0.000
0.024 0.024 0.000
0.015 0.015 0.000
0.010 0.010 0.000
0.006 0.006 0.000
V
35 °C 25 °C Increment
50.626 40.502 10.123
0.253 0.250 0.004
0.049 0.049 0.000
0.025 0.025 0.000
0.014 0.014 0.000
0.013 0.013 0.000
0.005 0.005 0.000
VII
35 °C 25 °C Increment
50.585 40.457 10.129
0.252 0.248 0.004
0.049 0.049 0.000
0.028 0.028 0.000
0.015 0.015 0.000
0.013 0.013 0.000
0.008 0.008 0.000
IX
35 °C 25 °C Increment
50.584 40.476 10.107
0.252 0.248 0.004
0.043 0.043 0.000
0.028 0.028 0.000
0.012 0.012 0.000
0.013 0.013 0.000
0.009 0.009 0.000
354
C. Pei et al. / Applied Thermal Engineering 120 (2017) 347–357
Table 6 Comparison of coefficients with different heat conductivity. Distribution No.
Coefficients
B
A(0)
A(1)
A(2)
A(3)
A(4)
A(5)
I
Original Chip doubled Carrier doubled
40.499 40.497 40.122
0.247 0.190 0.165
0.048 0.046 0.025
0.024 0.024 0.012
0.015 0.015 0.007
0.010 0.010 0.005
0.006 0.006 0.003
V
Original Chip doubled Carrier doubled
40.502 40.506 40.137
0.250 0.190 0.169
0.049 0.047 0.025
0.025 0.025 0.013
0.014 0.014 0.007
0.013 0.013 0.007
0.005 0.005 0.003
VII
Original Chip doubled Carrier doubled
40.457 40.462 40.097
0.248 0.188 0.169
0.049 0.047 0.025
0.028 0.028 0.014
0.015 0.015 0.008
0.013 0.013 0.007
0.008 0.008 0.004
IX
Original Chip doubled Carrier doubled
40.476 40.480 40.097
0.248 0.188 0.169
0.043 0.041 0.022
0.028 0.028 0.014
0.012 0.011 0.006
0.013 0.013 0.007
0.009 0.008 0.004
4.3.3. Resolution influence In some cases, in order to acquire a more accurate analysis result, the lattices are intended to be meshed smaller than the operating units. Nevertheless, the stability of the approach against different resolutions should be verified beforehand. The typical distributions and their simulation results are employed again. Bicubic interpolation is used to realize the zooming. Though the distribution is supposed to be binary in the former discussions, the approach could be to extend to non-integral spatial matrix actually. Back to the point, the binary matrix will be transformed into a real matrix during the zooming. In accordance with the change, Eqs. (3)–(5) and (7) are altered as below:
Hh;i ðgÞ ¼ 0
X
X h;i ðj; kÞjdh;i ðj; kÞ ¼ g
H1;1 ð0Þ H1;1 ð1Þ
ð9Þ
H1;1 ða þ b 1Þ
1
B H2;1 ð0Þ H2;1 ð1Þ H1;2 ða þ b 1Þ C C B C Z¼B .. . .. C B . A @ . . . Ha;b ða þ b 1Þ ab;aþb1 Ha;b ð0Þ Ha;b ð1Þ
ð10Þ
In Eq. (9), Xh,i(j, k) is a real number, in which the positive denotes heat source and the negative denotes cold source. Besides, Hh,i(0)
denotes the power of the lattice (h, i) itself. Taking distribution I as an example, the amount of the lattices shrinks from 400 to 324. With the zoomed temperature matrix and spatial matrix, the coefficients are fitted again. The results are displayed in Table 7. As could be observed, the ratio between the new coefficient A and original coefficient A is almost 1.1, which coincides with the expansion ratio. As coefficient A denotes the correlation, it is no wonder that it is proportional to the resolution. Comparatively, coefficient B stays stable, since it denotes the environment condition which will not vary with the resolution. It means that, the size of the lattices could be meshed smaller than the physical size of the units for a better accuracy, or could be enlarged to save calculation time. 4.3.4. Brief summary The relationship between coefficients and influencing factors is summarized in Table 8. In the table, it could be observed that coefficient B depends on the total power and the environment temperature, while coefficient A depends on material thermal properties and resolution (or the lattice size). Though the relationship is not necessary in the application of the proposed approach, it is deemed useful for the temperature prediction when these factors vary.
Table 7 Comparison of coefficients with different resolutions. Distribution No.
Coefficients
B
A(0)
A(1)
A(2)
A(3)
A(4)
A(5)
I
Original Zoomed Ratio
40.465 40.369 1.002
0.262 0.246 1.065
0.054 0.048 1.125
0.029 0.024 1.208
0.016 0.016 1.000
0.013 0.011 1.182
0.007 0.008 0.875
V
Original Zoomed Ratio
40.295 40.400 0.997
0.263 0.249 1.056
0.058 0.048 1.208
0.025 0.025 1.000
0.019 0.016 1.188
0.012 0.011 1.091
0.010 0.008 1.250
VII
Original Zoomed Ratio
40.277 40.328 0.999
0.259 0.243 1.066
0.054 0.048 1.125
0.028 0.025 1.120
0.015 0.015 1.000
0.013 0.011 1.182
0.008 0.008 1.000
IX
Original Zoomed Ratio
40.306 40.425 0.997
0.255 0.241 1.058
0.055 0.049 1.122
0.023 0.024 0.958
0.016 0.014 1.143
0.013 0.011 1.182
0.010 0.009 1.111
Table 8 Relationship between coefficients and influencing factors. Coefficient
Total power
Environment temperature
Heat conductivity of chip
Heat conductivity of die-attach or carrier
Resolution
B A(0) A except A(0)
d
d
s
s s
d d
‘’ denotes no influence, ‘d’ denotes proportional influence, ‘s’ denotes complex influence
C. Pei et al. / Applied Thermal Engineering 120 (2017) 347–357
5. Case study A microwave chip is chosen as an example, the top view of which is illustrated in Fig. 11(a). There are four main active regions on the chip, of which active region 1 is composed of transistor array. The power of active region 1 is the largest and the distribution depends on the chip’s operating modes. The chip is soldered to a low temperature co-fired ceramic carrier (LTCC) by a Mo-Cu eutectic solder. The actual heat conductivity of the solder is unknown, since it depends on the microstructure of die-attach greatly. To evaluate the thermal dissipation performance, infrared thermography is utilized. However, restricted by the testing condition, not all the operating modes could be realized during the measurement. One of the typical measure results is shown in Fig. 11(b). Particularly, the white box is an alignment tool in the measurement software and its color is irrelevant to temperature value, as well as the black box in Fig. 12(a). The power of the active regions is arranged as below: 8.25 W (N1), 0.715 W (N2 and N3), 0.8 W (N4). It could be seen that the power distribution of active region 1 is uniformly distributed in this mode. Besides, some low temperature regions (shows deep black in Fig. 12(a) and bright in (b)) are
355
ignored in the following calculation, since the real temperature is concealed by the passive components. To save calculation time, the chip is divided into 80 ⁄ 82 lattices. The coefficients in the case study are somewhat erratic, as shown in Fig. 12(a) and (b). By comparing Figs. 6(b) and 12(b), it could be found that the new oscillation period is longer than the former typical cases. Compared with the simulations, more lattices are set in the case study, so a bigger truncation order is recommended to be utilized. Afterwards, a new operating mode is shown in Fig. 13(a). The power of partial transistors increases from 0.093 W to 0.344 W, which generates a hot spot. The temperature distribution could be predicted by employing the distribution and the coefficients acquired in the former one (40-order truncated). For the sake of comparison, a real temperature measurement is also conducted and its image is converted into grey display. These two distributions and the residual errors are shown in Figs. 13 and 14. Except some low temperature regions caused by the passive components, the predicted result and the measured result coincide with each other well (the error of the highest temperature is 1.3 °C and the maximum error is less than 4.0 °C, considering a temperature rise of 50 °C).
Fig. 11. Top view of the chip in case study and measurement result: (a) top view, (b) temperature measurement (the white box is outlying).
Fig. 12. Coefficients convergence in case study: (a) coefficients convergence with distance; (b) convergence of A(0) under different truncation orders.
356
C. Pei et al. / Applied Thermal Engineering 120 (2017) 347–357
Fig. 13. Temperature distribution of measurement and prediction: (a) measurement result (the black box is outlying), (b) analysis result.
Acknowledgements The authors gratefully acknowledge the contribution of Prof. Suganuma, Osaka University for helping in the simulation and the case study. The work is supported by China Scholarship Council, (File No. 201606020080). Reference
Fig. 14. Residual error distribution in the case study.
6. Conclusion A simplified steady thermal analysis approach in this paper proposes, which is considered to be appropriate for devices with variable power. Though example power distribution and temperature measurement are demanded once to determine the unknown coefficients, the determination of material properties and boundary conditions could be relieved. Some characteristics of the approach are specified. The relationship between the coefficients in the approach and influencing factors is also determined, as well as the edge effect. In the case study, the approach is applied to an actual object. The case study show satisfactory results, of which the maximum error is not more than 5%. Compared to the conventional FE method, the approach proposed in this paper is more convenient and fast. Some shortcomings of the approach were also mentioned here. Firstly, the approach is based on linear superposition assumption, which means that the nonlinearity induced by temperature rise, i.e. the variations of material properties and boundary heat transfer coefficient, should be negligible. Meanwhile, in some cases, for instance, if the chip is covered by other passive components, the example temperature distribution might not be acquirable. Therefore, the proposed approach is considered not feasible for hybrid chips. In the future, the method to eliminate of edge effect will be studied. Meanwhile, the application of the approach to non-uniform material and transient thermal analysis is tended to be conducted.
[1] A.S. Fleischer, L. Chang, B.C. Johnson, The effect of die attach voiding on the thermal resistance of chip level packages, Microelectron. Reliab. 46 (2006) 794–804. [2] R. Liang, J. Zhang, S. Wang, Q. Chen, L. Xu, J. Dai, C. Chen, Investigation on thermal characterization of eutectic flip-chip UV-LEDs with different bonding voidage, IEEE Trans. Electron. Dev. 99 (2017) 1–6. [3] L. Ciamplini, M. Ciappa, P. Malberti, P. Regli, W. Fichtner, Modeling thermal effects of large contiguous voids in solder joints, Microelectr. J. 30 (1999) 1115–1123. [4] K.C. Otiaba, R.C. Bhatti, N.N. Ekere, S. Mallik, M.O. Alam, E.H. Amalu, M. Ekpu, Numerical study on thermal impacts of different void patterns on performance of chip-scale packaged power device, Microelectron. Reliab. 52 (2012) 1409– 1419. [5] K.C. Otiaba, M.I. Okereke, R.S. Bhatti, Numerical assessment of the effect of void morphology on thermo-mechanical performance of solder thermal interface material, Appl. Therm. Eng. 64 (2014) 51–63. [6] L. Yin, J. Zhang, P. Song, Y. Zhou, W. Yang, J. Zhang, The effects of void ratio in die attach layer on optical and thermal performances of high power light emitting diode, J. Nanoelectron. Optoe. 9 (2014) 1–6. [7] S.A. Paknejad, A. Mansourian, Y. Noh, K. Khtatba, S.H. Mannan, Thermally stable high temperature die attach solution, Mater. Design. 89 (2016) 1310–1314. [8] D.C. Katsis, J.D. vanWyk, A thermal, mechanical and electrical study of voiding in the solder die-attach of Power MOSFETs, IEEE T. Compon. Pack. Tecnol. 29 (2006) 127–136. [9] S.H. Tran, L. Dupont, Z. Khatir, Solder void position and size effects on electro thermal behavior of MOSFET transistor in forward bias conditions, Microelectron. Reliab. 54 (2014) 1921–1926. [10] L. Brunel, B. Lambert, P. Mezenge, Analysis of Schottky gate degradation evolution in AlGaN/GaN HEMTs during HTRB stress, Microelectron. Reliab. 53 (2013) 1450–1455. [11] Q. Li, X.Y. Zhu, Y.M. Xuan, Modeling heat generation in high power density nanometer scale GaAs/InGaAs/AlGaAs PHEMT, Int. J. Heat. Mass. Tran. 81 (2015) 130–136. [12] T. Niknam, M. Bornapour, A. Gheisari, Combined heat, power and hydrogen production optimal planning of fuel cell power plants in distribution networks, Energy Convers. Manage. 66 (2013) 11–25. [13] Y. Rahmani, D.D. Ganji, M.G. Bandpy, Analytical study of thermal spreading resistance in curved-edge heat spreader, Appl. Therm. Eng. 104 (2016) 527–533. [14] Y. Rahmani, H. Shokouhmand, Assessment of temperature-dependent conductivity effects on the thermal spreading/constriction resistance of semiconductors, J. Thermophys. Heat. Trans. 26 (2012) 638–643. [15] Y. Rahmani, H. Shokouhmand, A numerical study of thermal spreading/constriction resistance of silicon, ITherm (2012) 1–5. [16] T. Suzuki, T. Terasaki, Y. Kawana, D. Ishikawa, M. Nishimura, H. Nakao, K. Kurafuchi, Effect of manufacturing process on micro-deformation behavior of sintered-silver die-attach material, IEEE Trans. Device. Mat. Re. 16 (2016) 588–596.
C. Pei et al. / Applied Thermal Engineering 120 (2017) 347–357 [17] J. Carr, X. Mihlet, P. Gadaud, S.A.E. Boyer, G.E. Thompson, P. Lee, Quantitative characterization of porosity and determination of elastic modulus for sintered micro-silver joints, J. Mater. Process. Tech. 225 (2015) 19–23. [18] J. Sergent, A. Krum, Thermal Management Handbook: for Electronic Assemblies, McGraw-Hill, New York, 1998. [19] A.S.N. Huda, S. Taib, Application of infrared thermography for predictive/ preventive maintenance of thermal defect in electrical equipment, Appl. Therm. Eng. 61 (2013) 220–227. [20] Q.J. Tan, J.M. Dai, C.W. Bu, et al., Experimental study on debonding defects detection in thermal barrier coating structure using infrared lock-in thermographic technique, Appl. Therm. Eng. 107 (2016) 463–468.
357
[21] B. Gao, F. Yang, M.Y. Chen, et al., A temperature spectrum density distribution based condition evaluation method and application in IGBT, Appl. Therm. Eng. 106 (2016) 1440–1457. [22] A.D. Cliff, J.K. Ord, Spatial Processes: Models & Applications, Pion, New York, 1981. [23] Y.S. Jeong, S.J. Kim, M.K. Jeong, Automatic identification of defect patterns in semiconductor wafer maps using spatial correlogram and dynamic time warping, IEEE Trans. Semicond. Manuf. 21 (2008) 625–637. [24] W. Tobler, A computer movie simulating urban growth in the Detroit region, Econ. Geogr. 46 (1970) 234–240.
Contents lists available at ScienceDirect
Applied Thermal Engineering journal homepage: www.elsevier.com/locate/apthermeng
Research Paper
Simplified approach for steady thermal analysis of chips with variable power based on spatial autocorrelation C. Pei ⇑, G.C. Fu, Y.H. Zhao, B. Wan, Y. Cheng School of Reliability and Systems Engineering, Beihang University, Beijing 100191, China
h i g h l i g h t s Approach for steady thermal analysis of chip with variable power is proposed. The determination of material properties and other factors can be tolerated. The highest temperature error in the case study is not more than 2%.
a r t i c l e
i n f o
Article history: Received 20 December 2016 Accepted 4 April 2017 Available online 5 April 2017 Keywords: Steady thermal analysis Spatial autocorrelation Undetermined coefficient method Power devices
a b s t r a c t Due to their growing power and shrinking size, thermal management and analysis have become one of the most critical concerns for power devices. Though finite element (FE) method is commonly employed in thermal analysis now, the uncertainty of material properties and the redundancy of repeat modeling still nag analyzers. In this paper, a simplified approach for steady thermal analysis of chips with variable power is proposed. Based on spatial autocorrelation and undetermined method, the steady temperature distribution could be directly determined with power distribution and example temperature measurement result, without identifying environment factors or material properties. The theory of the simplified approach is introduced at first in this paper. The characteristics of the simplified approach are discussed afterwards. Eventually, a case study is given to illustrate the procedure, of which the analytical result coincides with the actual measurement well. Ó 2017 Elsevier Ltd. All rights reserved.
1. Introduction As to the continual miniaturization and increasing power dissipation, thermal management has become a critical factor in the design process of microcircuits. Especially for high power microcircuits, i.e. insulated gate bipolar transistor (IGBT), high temperature, induced by heat flow, will cause parameter drifting or even failures. Numerous researches have been conducted on die-attach layer or its inner defects. Within the researches, the size [1,2], position [3], pattern [4], and even depth [5] of the die-attach defects have been discussed, and their influences on thermal [6], mechanical [7,8] and parameter-stable [9,10] performance have been also studied. In most of those researches, heat sources have been treated as an entity. However, the heat generation distribution affects thermal analysis results at the same time. Hot spots may be concealed by the assumption of uniform heating. For devices composed of plentiful repetition units, i.e. IGBTs or memories, the ⇑ Corresponding author. E-mail address: [email protected] (C. Pei). http://dx.doi.org/10.1016/j.applthermaleng.2017.04.010 1359-4311/Ó 2017 Elsevier Ltd. All rights reserved.
distribution of heat sources is also determined by the operating mode [11,12]. Meanwhile, current thermal analyses of devices are mainly based on equivalent thermal resistance network or finite element (FE) method [13–15]. No matter which method is employed, dissipation factors, including material property, dimension parameter, power distribution and boundary condition, should be determined in advance. The precision of the thermal analysis depends on the accuracy of the factors set in the analysis models to a great extent, which may be difficult to acquire sometimes. For example, the heat conductivity of the die-attach depends on both the material and the process, which will vary with each individual, especially in the case of porous die-attach pastes [16,17]. At the same time, whenever the operating mode switches, the analysis models have to be modified or rebuilt correspondingly, which will also increase the modeling work. To cope with the problems, an approach that could simplify the thermal analysis of chips with variable power will be proposed in this paper. With this approach, the steady temperature distribution could be acquired only with the information of power distribution
348
C. Pei et al. / Applied Thermal Engineering 120 (2017) 347–357
in different operating modes and the example temperature measurement result. The paper is arranged as below. Firstly, the theory of the simplified approach will be introduced. Subsequently, a series of typical power distributions will be given, and their thermal influences will be solved by the FE method and the proposed approach. Their temperature distributions will be compared to verify the accuracy of the approach. Meanwhile, the relationship between the undetermined coefficients and other dissipation factors will be discussed. Finally, a case study will be given. The steady temperature distribution measured by infrared microscopy and the analytical result based on this approach will be compared again to verify the approach. 2. Theory & methodology The typical heat dissipation route for encapsulated chips is illustrated in Fig. 1(a). As to the low heat conductivity of thermal molding compound, most of heat dissipates from the die-attach layer through the carrier to the mounting plane or the heat sink [18]. In a simplified form, the temperature of chip is determined by power distribution, power intensity, material thermal properties and boundary conditions, as Fig. 1(a) shows. In this paper, spatial analysis is utilized to depict the power distribution and identify the influence, which has been employed in other works on defect detection [19,20] or condition evaluation [21]. Among spatial analysis methods, spatial autocorrelation method, proposed by A. D. Cliff and J. K. Ord [22] in 1980s, is most promising and has been applied to wafer defect detection [23]. In these studies, the defective dies, fabricated on a single wafer, were summed as a wafer map to display the spatial dependence across the wafer. Based on it, spatial autocorrelation was extracted and spatial effect could be resolved. Similar to wafer map, chips with variable power will be also meshed into matrix-arranged lattices in this paper, as Fig. 2 illustrates. Let matrix Ma,b denote the chip. The power of unit (i, j) is expressed as X (i, j). The definition of the power of each unit is:
Xði; jÞ ¼
1; work 0; non work
ð1Þ
For simplification, the power of each unit is assumed to be equal and Xa,b is transformed into a 0–1 matrix. The distance was measured by counting the minimum steps between two units. In other words, the heat is only permitted to transfer between adjacent units. For any unit (h, i) and any other unit (j, k), their distance could be expressed as:
dh;i ðj; kÞ ¼ jh jj þ ji kj
ð2Þ
For unit (h, i), the units at d distance are separately summed as:
W 0h;i ðgÞ ¼ fXðj; kÞjdh;i ðj; kÞ ¼ g; Xðj; kÞ ¼ 0; g P 1g
ð3Þ
W 1h;i ðgÞ ¼ fXðj; kÞjdh;i ðj; kÞ ¼ g; Xðj; kÞ ¼ 1; g P 1g
ð4Þ
H1h;i ðgÞ ¼ jW 1h;i ðgÞj
ð5Þ
|| in Eq. (5) denotes the total number of set elements. Particularly, Hh,i(0) denotes the operating condition of unit (h, i) itself. According to Fourier’s law, the temperature rise of units is dependent on the thermal resistance, heat generation, and boundary heat flux. For each unit in the chip, excepting the heat generated interiorly and the heat flow exteriorly the chip bottom, the heat inflow (or outflow) from neighboring units should also be considered. Actually, the hot spots appear not only in the area with high heat intensity but also in the operation units ‘‘crowded” area. According to the first law of geography, ‘‘everything is related to everything else, but near things are more related than distant things” [24], the synergy influences between units should be considered and converge with distance. The simplified equivalent circuit of the dissipation route is illustrated in Fig. 1(b). According to superposition theorem, if temperature variation is limited and the thermal material properties remain stable, the steady temperature rise of the chip could be expressed as the linear superposition of influences of all heat sources (or operating units). Based on the hypothesis, the steady temperature matrix Ta,b could be expressed as below:
T ¼ AZ þ B
Fig. 1. Dissipation route and equivalent circuit: (a) Typical dissipation route for encapsulated chips (b) the simplified equivalent circuit.
Fig. 2. Chip meshed into matrix-arranged lattices.
ð6Þ
349
C. Pei et al. / Applied Thermal Engineering 120 (2017) 347–357
1 H11;1 ða þ b 1Þ H1;1 ð0Þ H11;1 ð1Þ C B B H2;1 ð0Þ H12;1 ð1Þ H11;2 ða þ b 1Þ C C B Z¼B C .. .. .. C B A @ . . . H1a;b ða þ b 1Þ ab;aþb1 Ha;b ð0Þ H1a;b ð1Þ 0
0 B B A¼B B @
Að0Þ Að1Þ .. .
3. Modeling
ð7Þ
1 C C C C A
ð8Þ
Aða þ b 1Þ Within the undetermined coefficients, A denotes the spatial correlation and the heat intensity. B denotes the environment factors, which is assumed to include boundary conditions related factor, i.e. boundary temperature and boundary heat exchange. As the power of each unit is normalized, not only matrix Z but also coefficient A affects the total power. The former represents the total power area, and the latter represents the heat intensity. Vector T denotes the temperature distribution. Obviously, with the power distribution of the chip and the example temperature measurement result, undetermined coefficient A and B in Eq. (6) could be solved. And if A and B known, the steady temperature distribution of the chip under any power distribution could be easily determined, during which the repetitive modeling and determination of material properties could be avoided. The flowchart of the proposed simplified approach is displayed in Fig. 3.
3.1. Object Before applying the proposed approach, the accuracy of this approach should to be confirmed, as well as its other characteristics. Several typical distributions are chosen here, based on which a series of numerical models are established to acquire reference value for the confirmation and discussions. To eliminate the effects of the package, a model only consisting of chip, carrier and heat sink is established, as Fig. 4 illustrates. For simplification, the dieattach layer supposed to be homogenous. The dimensions of the layers are listed in Table 1. The thermal properties are also listed in the table, which are acquired from [18]. The environment temperature of the thermal analysis is set to 25 °C, as well as the air. Only natural convection heat dissipation is considered in the model and the direction of gravity is set from the chip to the heatsink. Radiation dissipation is ignored in the model, as heat transfer by radiation is negligible under a low temperature like this case.
Table 1 Dimension and material of example chip. Layer
Material
Heat conductivity (W/m k)
Density (g/m3)
Chip Carrier Heat sink
GaAs Duralumin Duralumin
59.1 180 180
5310 2700 2700
Fig. 3. Flowchart of the Simplified approach.
Fig. 4. Structure diagram of numerical model.
350
C. Pei et al. / Applied Thermal Engineering 120 (2017) 347–357
3.2. Power distribution The chip is divided into 400 (20 ⁄ 20) lattices symbolically in its horizontal plane. As mentioned before, the power of each unit is supposed to be equal, which is set to 0.01 W in this part. The next step is to determine several typical distributions. In [23], four typical distributions are chosen. While the circle pattern is most common in wafer process, concentration is the most fatal case in chip dissipation. To take the pertinence in account, the circle pattern is substituted by center and edge concentration. Additionally, edge concentrated power distribution is also selected to investigate the edge effect, though it seldom appears in chips effect. Eventually, five distributions, including spatially homogeneous Bernoulli process, edge concentration, center concentration, cluster, and repetition, are established in the paper, which are displayed in Table 2. 3.3. Numerical settings Commercial thermal analysis software FloTHERM is employed in this paper. Finite volume method (FVM) is employed. Localized grids are toggled to reduce the amount of elements and assure the quality of the meshing. The elements of all the meshes are about 456,000. Multi grid method is utilized in solving to accelerate the
convergence. By default, if the residual error within the calculation reached 0.5% of total source value, the solution is deemed converged, as well as the velocity and pressure. 4. Results and discussions 4.1. Characterization of coefficients All the numerical models converge within 200 iterations. Their steady temperature distributions are shown in Fig. 5. As could be seen, the temperature distributions vary with the different power distributions. The highest temperature appears in the edge concentration distribution, which reveals that the edge effect is the most significant factor affecting the temperature rise. Besides, among the other four distributions, the highest and lowest temperatures both appear in the center concentration distribution. The uneven temperature distribution may be induced by power concentration. In other words, the concentration effect is another factor that should be considered. As mentioned in Section 2, the key procedure of the proposed approach is to determine underdetermined coefficients A and B in Eq. (6). Before utilizing the coefficients to predict the temperature distributions under different power distributions, their stabil-
Table 2 Typical power distributions. Percentage of operating unit
SHBP
Edge concentration
Center concentration
Cluster
Repetition
No. 30%
I
II
III
IV
V
No.
VI
VII
VIII
IX
X
20%
Fig. 5. Temperature distribution of by simulation: (a) distribution I, (b) distribution II, (c) distribution III, (d) distribution IV, (e) distribution V, (f) distribution VI, (g) distribution VII, (h) distribution VIII, (i) distribution IX, (j) distribution X.
351
C. Pei et al. / Applied Thermal Engineering 120 (2017) 347–357
ity should be validated beforehand. To achieve this goal, the spatial correlogram matrix of all the typical distributions are extracted. Then, A and B are solved with the thermal simulation analysis results by the least square method. Considering that the spatial correlogram matrix Z may not be full rank, underdetermined coefficient A is suggested to be truncated by a specific distance. Before taking this strategy, the influence of the truncation should also be analyzed in advance. As mentioned above, underdetermined coefficient A, correlation influences between units, is assumed to converge with distance. Moreover, it is assumed to keep stable under different truncations. Next, the assumptions will be proved. Taking distribution VII as an example, the fitting results under the different truncations are shown in Table 3. As could be seen, the coefficient A decreases with the distance. Then the relationship between coefficient A under 10order truncation and the distance is plotted in Fig. 6(a). An exponential function is employed, whose form is:
f ðgÞ ¼ ag b þ c
ð9Þ
The fitting curve is also included in Fig. 6(a). The goodness demonstrates that the data are well fitted to an exponential function. Besides, as symbol b denotes convergence rate, the convergence rate of A(g) is considered to be 2.25. Considering that Z(g) is proportional to distance g, their combined influence is considered to converge and the rate is 1.25. In other words, the influences of the coefficient A decrease with the distance, which means that the error induced by truncation could be ignored. Then the stability of the coefficients against the different truncation orders is investigated. Distribution VII is taken as an
example again, and the trend of A(0) values under the different truncation orders is plotted in Fig. 6(b). It could be observe that A(0) fluctuates at first. Then it stays stable briefly until the distance reaches the rank of matrix Z, which verifies the stability and uniqueness of the coefficients. Since the reasonability of the truncation has been proved, coefficient A and B in all the distributions by 4-order truncated are selected to compare, as shown in Table 4. The comparison reveals that coefficient B is proportional to the total power and A is immune to the power variation. Besides, the result demonstrates that nearly all the coefficients are independent from the power distribution. 4.2. Accuracy and improvement The average value of the coefficients, except distribution III and IV’s is chosen as the general value. Then the general value are substituted into Eq. (6) to predict the temperature. The residual error between the analysis and the simulation is displayed in Fig. 7. It could be found that the distribution of the residual error is independent from the distribution of the operating units and concentrates on the edge. The residual errors in distribution III and IV are mostly significant. This suggests that the dissipation of the edge units is much more complex than that of interior units. As the heat transfer coefficient of the chip edge is far less than the die-attach, the edge could be treated as an adiabatic boundary. Contrary to the conventional treatment in FEA, the chip is expanded to eliminate the edge influence, as Fig. 8(a) illustrates. The space of the chip is divided into 4 quadrants. Each quadrant
Table 3 Fitted coefficient A in distribution VII under different truncation orders. Truncation order
A(0)
A(1)
A(2)
A(3)
A(4)
A(5)
A(6)
A(7)
A(8)
A(9)
1 2 3 4 5 6 7 8 9 10
0.354 0.219 0.247 0.249 0.248 0.249 0.247 0.249 0.249 0.249
0.084 0.037 0.046 0.047 0.045 0.047 0.047 0.048 0.048
0.041 0.017 0.021 0.024 0.024 0.024 0.023 0.023
0.025 0.010 0.013 0.014 0.013 0.013 0.013
0.020 0.008 0.008 0.009 0.009 0.009
0.014 0.007 0.008 0.007 0.007
0.009 0.004 0.005 0.005
0.005 0.003 0.004
0.003 0.002
0.003
Fig. 6. Convergence of coefficient A: (a) coefficient A convergence with distance; (b) convergence of A(0) under different truncation orders.
352
C. Pei et al. / Applied Thermal Engineering 120 (2017) 347–357
Table 4 Comparison of coefficients with different distributions and total powers. Distribution No.
Total power (W)
B (°C)
Temperature rise (°C)
A(0)
A(1)
A(2)
A(3)
I II III IV V VI VII VIII IX X
0.8 1.2 0.8 1.2 0.8 1.2 0.8 1.2 0.8 1.2
40.634 46.784 40.971 46.608 40.509 46.567 40.568 47.109 40.548 46.661
15.634 21.784 15.971 21.608 15.509 21.567 15.568 22.109 15.548 21.661
0.257 0.258 0.765 0.447 0.250 0.257 0.248 0.245 0.246 0.247
0.048 0.047 0.338 0.102 0.049 0.053 0.047 0.051 0.048 0.047
0.020 0.018 0.216 0.042 0.025 0.027 0.024 0.025 0.023 0.025
0.009 0.008 0.073 0.026 0.017 0.015 0.013 0.015 0.014 0.015
Fig. 7. Residual error before expansion: (a) distribution I, (b) distribution II, (c) distribution III, (d) distribution IV, (e) distribution V, (f) distribution VI, (g) distribution VII, (h) distribution VIII, (i) distribution IX, (j) distribution X.
Fig. 8. Mirror expansion (a) schematic diagram (symbol ‘F’ represents the original distribution); (b) mirror expansion of distribution I.
is mirrored upside down, left and right. In this way, the original edge could be assured to be adiabatic. Distribution I is taken as an example. Its expanded distribution is shown in Fig. 8(b). The area outlined in yellow (top-left) is the original distribution and the area outlined in red (center) is the expanded distribution from the second quadrant. All the typical distributions are expanded as above. Then, the coefficients are solved again, and the convergence of coefficient A is shown in Fig. 9. New general coefficients are employed to
determine the residual errors by the mirror expansion. The results are illustrated in Fig. 10. By the comparison, it could be found that the residual errors decline from 0.2 to 0.1. Besides, the anomalous units located in the edge are also reduced. Nevertheless, though the anomaly caused by edge influence is relieved slightly, coefficient A in distribution III and IV are still larger than those of other distributions. Considering that power units are not likely to concentrate in the edge, the error induced by the edge influence may still be acceptable.
353
C. Pei et al. / Applied Thermal Engineering 120 (2017) 347–357
Firstly, distributions I, V, VII and IX are chosen. Higher environment temperature and fluid temperature, 35 °C, are set to analyze the influences. Mirror expansion method is employed in the calculation and the fitting results are displayed in Table 5. With the results, it could be observed that only coefficient B increases and the increment is close to the rise of the environment temperature. Oppositely, coefficient A is insensitive to the change of environment temperature.
Fig. 9. Convergence of new coefficient A in expanded distributions.
4.3. Relationship between coefficients and influencing factors 4.3.1. Environment influence In Section 4.1, the influences of heat sources (total power) on the coefficients have been investigated. Here, the relationship between the material properties or boundary conditions and the coefficients will be discussed further.
4.3.2. Material influence The influences of heat conductivity are also studied below. The heat conductivity of the chip and the carrier is doubled separately, which is depicted in Table 1. The new coefficient fitting results are displayed in Table 6. It could be observed that the increase of chip heat conductivity only affects coefficient A(0). By contrast, the variation of carrier heat conductivity will exert an influence on all the coefficients. Coefficient A nearly halves, except A(0). These differences could be explained by the different order of magnitudes of the resistances. The correlation resistance between units Rco in the Fig. 1(b), is about 600 K/W, while the exterior resistance Rout (Rtrans and Rbound) is about 160 K/W, of which nearly 90% is the inner resistance of carrier. Therefore, all of coefficient A are influenced by the heat conductivity of the carrier, as the influences of Rout are more global. By contrast, The variation of the ratio between A(0) and other coefficient A demonstrates that Rco only affect the heat flow between power units.
Fig. 10. Residual error after expansion: (a) distribution I, (b) distribution II, (c) distribution III, (d) distribution IV, (e) distribution V, (f) distribution VI, (g) distribution VII, (h) distribution VIII, (i) distribution IX, (j) distribution X.
Table 5 Comparison of coefficients under different environment temperature. Distribution No.
Coefficient
B
A(0)
A(1)
A(2)
A(3)
A(4)
A(5)
I
35 °C 25 °C Increment
50.615 40.499 10.120
0.250 0.247 0.003
0.048 0.048 0.000
0.024 0.024 0.000
0.015 0.015 0.000
0.010 0.010 0.000
0.006 0.006 0.000
V
35 °C 25 °C Increment
50.626 40.502 10.123
0.253 0.250 0.004
0.049 0.049 0.000
0.025 0.025 0.000
0.014 0.014 0.000
0.013 0.013 0.000
0.005 0.005 0.000
VII
35 °C 25 °C Increment
50.585 40.457 10.129
0.252 0.248 0.004
0.049 0.049 0.000
0.028 0.028 0.000
0.015 0.015 0.000
0.013 0.013 0.000
0.008 0.008 0.000
IX
35 °C 25 °C Increment
50.584 40.476 10.107
0.252 0.248 0.004
0.043 0.043 0.000
0.028 0.028 0.000
0.012 0.012 0.000
0.013 0.013 0.000
0.009 0.009 0.000
354
C. Pei et al. / Applied Thermal Engineering 120 (2017) 347–357
Table 6 Comparison of coefficients with different heat conductivity. Distribution No.
Coefficients
B
A(0)
A(1)
A(2)
A(3)
A(4)
A(5)
I
Original Chip doubled Carrier doubled
40.499 40.497 40.122
0.247 0.190 0.165
0.048 0.046 0.025
0.024 0.024 0.012
0.015 0.015 0.007
0.010 0.010 0.005
0.006 0.006 0.003
V
Original Chip doubled Carrier doubled
40.502 40.506 40.137
0.250 0.190 0.169
0.049 0.047 0.025
0.025 0.025 0.013
0.014 0.014 0.007
0.013 0.013 0.007
0.005 0.005 0.003
VII
Original Chip doubled Carrier doubled
40.457 40.462 40.097
0.248 0.188 0.169
0.049 0.047 0.025
0.028 0.028 0.014
0.015 0.015 0.008
0.013 0.013 0.007
0.008 0.008 0.004
IX
Original Chip doubled Carrier doubled
40.476 40.480 40.097
0.248 0.188 0.169
0.043 0.041 0.022
0.028 0.028 0.014
0.012 0.011 0.006
0.013 0.013 0.007
0.009 0.008 0.004
4.3.3. Resolution influence In some cases, in order to acquire a more accurate analysis result, the lattices are intended to be meshed smaller than the operating units. Nevertheless, the stability of the approach against different resolutions should be verified beforehand. The typical distributions and their simulation results are employed again. Bicubic interpolation is used to realize the zooming. Though the distribution is supposed to be binary in the former discussions, the approach could be to extend to non-integral spatial matrix actually. Back to the point, the binary matrix will be transformed into a real matrix during the zooming. In accordance with the change, Eqs. (3)–(5) and (7) are altered as below:
Hh;i ðgÞ ¼ 0
X
X h;i ðj; kÞjdh;i ðj; kÞ ¼ g
H1;1 ð0Þ H1;1 ð1Þ
ð9Þ
H1;1 ða þ b 1Þ
1
B H2;1 ð0Þ H2;1 ð1Þ H1;2 ða þ b 1Þ C C B C Z¼B .. . .. C B . A @ . . . Ha;b ða þ b 1Þ ab;aþb1 Ha;b ð0Þ Ha;b ð1Þ
ð10Þ
In Eq. (9), Xh,i(j, k) is a real number, in which the positive denotes heat source and the negative denotes cold source. Besides, Hh,i(0)
denotes the power of the lattice (h, i) itself. Taking distribution I as an example, the amount of the lattices shrinks from 400 to 324. With the zoomed temperature matrix and spatial matrix, the coefficients are fitted again. The results are displayed in Table 7. As could be observed, the ratio between the new coefficient A and original coefficient A is almost 1.1, which coincides with the expansion ratio. As coefficient A denotes the correlation, it is no wonder that it is proportional to the resolution. Comparatively, coefficient B stays stable, since it denotes the environment condition which will not vary with the resolution. It means that, the size of the lattices could be meshed smaller than the physical size of the units for a better accuracy, or could be enlarged to save calculation time. 4.3.4. Brief summary The relationship between coefficients and influencing factors is summarized in Table 8. In the table, it could be observed that coefficient B depends on the total power and the environment temperature, while coefficient A depends on material thermal properties and resolution (or the lattice size). Though the relationship is not necessary in the application of the proposed approach, it is deemed useful for the temperature prediction when these factors vary.
Table 7 Comparison of coefficients with different resolutions. Distribution No.
Coefficients
B
A(0)
A(1)
A(2)
A(3)
A(4)
A(5)
I
Original Zoomed Ratio
40.465 40.369 1.002
0.262 0.246 1.065
0.054 0.048 1.125
0.029 0.024 1.208
0.016 0.016 1.000
0.013 0.011 1.182
0.007 0.008 0.875
V
Original Zoomed Ratio
40.295 40.400 0.997
0.263 0.249 1.056
0.058 0.048 1.208
0.025 0.025 1.000
0.019 0.016 1.188
0.012 0.011 1.091
0.010 0.008 1.250
VII
Original Zoomed Ratio
40.277 40.328 0.999
0.259 0.243 1.066
0.054 0.048 1.125
0.028 0.025 1.120
0.015 0.015 1.000
0.013 0.011 1.182
0.008 0.008 1.000
IX
Original Zoomed Ratio
40.306 40.425 0.997
0.255 0.241 1.058
0.055 0.049 1.122
0.023 0.024 0.958
0.016 0.014 1.143
0.013 0.011 1.182
0.010 0.009 1.111
Table 8 Relationship between coefficients and influencing factors. Coefficient
Total power
Environment temperature
Heat conductivity of chip
Heat conductivity of die-attach or carrier
Resolution
B A(0) A except A(0)
d
d
s
s s
d d
‘’ denotes no influence, ‘d’ denotes proportional influence, ‘s’ denotes complex influence
C. Pei et al. / Applied Thermal Engineering 120 (2017) 347–357
5. Case study A microwave chip is chosen as an example, the top view of which is illustrated in Fig. 11(a). There are four main active regions on the chip, of which active region 1 is composed of transistor array. The power of active region 1 is the largest and the distribution depends on the chip’s operating modes. The chip is soldered to a low temperature co-fired ceramic carrier (LTCC) by a Mo-Cu eutectic solder. The actual heat conductivity of the solder is unknown, since it depends on the microstructure of die-attach greatly. To evaluate the thermal dissipation performance, infrared thermography is utilized. However, restricted by the testing condition, not all the operating modes could be realized during the measurement. One of the typical measure results is shown in Fig. 11(b). Particularly, the white box is an alignment tool in the measurement software and its color is irrelevant to temperature value, as well as the black box in Fig. 12(a). The power of the active regions is arranged as below: 8.25 W (N1), 0.715 W (N2 and N3), 0.8 W (N4). It could be seen that the power distribution of active region 1 is uniformly distributed in this mode. Besides, some low temperature regions (shows deep black in Fig. 12(a) and bright in (b)) are
355
ignored in the following calculation, since the real temperature is concealed by the passive components. To save calculation time, the chip is divided into 80 ⁄ 82 lattices. The coefficients in the case study are somewhat erratic, as shown in Fig. 12(a) and (b). By comparing Figs. 6(b) and 12(b), it could be found that the new oscillation period is longer than the former typical cases. Compared with the simulations, more lattices are set in the case study, so a bigger truncation order is recommended to be utilized. Afterwards, a new operating mode is shown in Fig. 13(a). The power of partial transistors increases from 0.093 W to 0.344 W, which generates a hot spot. The temperature distribution could be predicted by employing the distribution and the coefficients acquired in the former one (40-order truncated). For the sake of comparison, a real temperature measurement is also conducted and its image is converted into grey display. These two distributions and the residual errors are shown in Figs. 13 and 14. Except some low temperature regions caused by the passive components, the predicted result and the measured result coincide with each other well (the error of the highest temperature is 1.3 °C and the maximum error is less than 4.0 °C, considering a temperature rise of 50 °C).
Fig. 11. Top view of the chip in case study and measurement result: (a) top view, (b) temperature measurement (the white box is outlying).
Fig. 12. Coefficients convergence in case study: (a) coefficients convergence with distance; (b) convergence of A(0) under different truncation orders.
356
C. Pei et al. / Applied Thermal Engineering 120 (2017) 347–357
Fig. 13. Temperature distribution of measurement and prediction: (a) measurement result (the black box is outlying), (b) analysis result.
Acknowledgements The authors gratefully acknowledge the contribution of Prof. Suganuma, Osaka University for helping in the simulation and the case study. The work is supported by China Scholarship Council, (File No. 201606020080). Reference
Fig. 14. Residual error distribution in the case study.
6. Conclusion A simplified steady thermal analysis approach in this paper proposes, which is considered to be appropriate for devices with variable power. Though example power distribution and temperature measurement are demanded once to determine the unknown coefficients, the determination of material properties and boundary conditions could be relieved. Some characteristics of the approach are specified. The relationship between the coefficients in the approach and influencing factors is also determined, as well as the edge effect. In the case study, the approach is applied to an actual object. The case study show satisfactory results, of which the maximum error is not more than 5%. Compared to the conventional FE method, the approach proposed in this paper is more convenient and fast. Some shortcomings of the approach were also mentioned here. Firstly, the approach is based on linear superposition assumption, which means that the nonlinearity induced by temperature rise, i.e. the variations of material properties and boundary heat transfer coefficient, should be negligible. Meanwhile, in some cases, for instance, if the chip is covered by other passive components, the example temperature distribution might not be acquirable. Therefore, the proposed approach is considered not feasible for hybrid chips. In the future, the method to eliminate of edge effect will be studied. Meanwhile, the application of the approach to non-uniform material and transient thermal analysis is tended to be conducted.
[1] A.S. Fleischer, L. Chang, B.C. Johnson, The effect of die attach voiding on the thermal resistance of chip level packages, Microelectron. Reliab. 46 (2006) 794–804. [2] R. Liang, J. Zhang, S. Wang, Q. Chen, L. Xu, J. Dai, C. Chen, Investigation on thermal characterization of eutectic flip-chip UV-LEDs with different bonding voidage, IEEE Trans. Electron. Dev. 99 (2017) 1–6. [3] L. Ciamplini, M. Ciappa, P. Malberti, P. Regli, W. Fichtner, Modeling thermal effects of large contiguous voids in solder joints, Microelectr. J. 30 (1999) 1115–1123. [4] K.C. Otiaba, R.C. Bhatti, N.N. Ekere, S. Mallik, M.O. Alam, E.H. Amalu, M. Ekpu, Numerical study on thermal impacts of different void patterns on performance of chip-scale packaged power device, Microelectron. Reliab. 52 (2012) 1409– 1419. [5] K.C. Otiaba, M.I. Okereke, R.S. Bhatti, Numerical assessment of the effect of void morphology on thermo-mechanical performance of solder thermal interface material, Appl. Therm. Eng. 64 (2014) 51–63. [6] L. Yin, J. Zhang, P. Song, Y. Zhou, W. Yang, J. Zhang, The effects of void ratio in die attach layer on optical and thermal performances of high power light emitting diode, J. Nanoelectron. Optoe. 9 (2014) 1–6. [7] S.A. Paknejad, A. Mansourian, Y. Noh, K. Khtatba, S.H. Mannan, Thermally stable high temperature die attach solution, Mater. Design. 89 (2016) 1310–1314. [8] D.C. Katsis, J.D. vanWyk, A thermal, mechanical and electrical study of voiding in the solder die-attach of Power MOSFETs, IEEE T. Compon. Pack. Tecnol. 29 (2006) 127–136. [9] S.H. Tran, L. Dupont, Z. Khatir, Solder void position and size effects on electro thermal behavior of MOSFET transistor in forward bias conditions, Microelectron. Reliab. 54 (2014) 1921–1926. [10] L. Brunel, B. Lambert, P. Mezenge, Analysis of Schottky gate degradation evolution in AlGaN/GaN HEMTs during HTRB stress, Microelectron. Reliab. 53 (2013) 1450–1455. [11] Q. Li, X.Y. Zhu, Y.M. Xuan, Modeling heat generation in high power density nanometer scale GaAs/InGaAs/AlGaAs PHEMT, Int. J. Heat. Mass. Tran. 81 (2015) 130–136. [12] T. Niknam, M. Bornapour, A. Gheisari, Combined heat, power and hydrogen production optimal planning of fuel cell power plants in distribution networks, Energy Convers. Manage. 66 (2013) 11–25. [13] Y. Rahmani, D.D. Ganji, M.G. Bandpy, Analytical study of thermal spreading resistance in curved-edge heat spreader, Appl. Therm. Eng. 104 (2016) 527–533. [14] Y. Rahmani, H. Shokouhmand, Assessment of temperature-dependent conductivity effects on the thermal spreading/constriction resistance of semiconductors, J. Thermophys. Heat. Trans. 26 (2012) 638–643. [15] Y. Rahmani, H. Shokouhmand, A numerical study of thermal spreading/constriction resistance of silicon, ITherm (2012) 1–5. [16] T. Suzuki, T. Terasaki, Y. Kawana, D. Ishikawa, M. Nishimura, H. Nakao, K. Kurafuchi, Effect of manufacturing process on micro-deformation behavior of sintered-silver die-attach material, IEEE Trans. Device. Mat. Re. 16 (2016) 588–596.
C. Pei et al. / Applied Thermal Engineering 120 (2017) 347–357 [17] J. Carr, X. Mihlet, P. Gadaud, S.A.E. Boyer, G.E. Thompson, P. Lee, Quantitative characterization of porosity and determination of elastic modulus for sintered micro-silver joints, J. Mater. Process. Tech. 225 (2015) 19–23. [18] J. Sergent, A. Krum, Thermal Management Handbook: for Electronic Assemblies, McGraw-Hill, New York, 1998. [19] A.S.N. Huda, S. Taib, Application of infrared thermography for predictive/ preventive maintenance of thermal defect in electrical equipment, Appl. Therm. Eng. 61 (2013) 220–227. [20] Q.J. Tan, J.M. Dai, C.W. Bu, et al., Experimental study on debonding defects detection in thermal barrier coating structure using infrared lock-in thermographic technique, Appl. Therm. Eng. 107 (2016) 463–468.
357
[21] B. Gao, F. Yang, M.Y. Chen, et al., A temperature spectrum density distribution based condition evaluation method and application in IGBT, Appl. Therm. Eng. 106 (2016) 1440–1457. [22] A.D. Cliff, J.K. Ord, Spatial Processes: Models & Applications, Pion, New York, 1981. [23] Y.S. Jeong, S.J. Kim, M.K. Jeong, Automatic identification of defect patterns in semiconductor wafer maps using spatial correlogram and dynamic time warping, IEEE Trans. Semicond. Manuf. 21 (2008) 625–637. [24] W. Tobler, A computer movie simulating urban growth in the Detroit region, Econ. Geogr. 46 (1970) 234–240.