A Thermal Impedance Model for IGBT Modules Considering the Nonlinear Thermal Characteristics of Chips and Ceramic Materials

Abstract

The traditional method of calculating junction temperature does not consider the dependence of a material’s thermal conductivity on temperature, in which the thermal conductivity changes with temperature. However, with an increase in junction temperature, the temperature sensitivity (TS) will have a more significant impact on the actual temperature of chips. This study established an improved IGBT equivalent thermal impedance model that considers the nonlinear characteristics of the TS of chips and ceramic materials. The Fourier series analysis method was used to obtain the heat flux density curve, and then the heat diffusion angles of each layer were solved. Moreover, iterations were performed until the thermal conductivity and temperature of the chip and ceramic layers matched the nonlinear characteristics of the TS. When the power loss was less than 200 W, the maximum error of the junction temperature calculated by the proposed method considering TS was 3%, while the maximum error of the method without considering TS was 9.5%. Compared with the finite element simulation, the proposed method has a faster solving speed and high accuracy. The proposed method only requires the input material parameters, size parameters, and boundary conditions to solve the junction temperature, which has strong practicality and high accuracy.

1. Introduction

Insulated-gate bipolar transistor (IGBT) modules play a significant role in modern power electronics technology. They are widely used in various industrial fields, including electric vehicles [1,2], aerospace [3], and wind power generation [4,5]. IGBT devices are widely used in rail transport. However, owing to the high probability of failure, research related to fault diagnosis is more extensive, which is an important reason for conducting device reliability studies [6,7,8]. However, as the power level and power density of IGBT modules increase, the reliability of IGBT modules is increasingly affected by the device junction temperature [9,10]. It is important to accurately calculate the junction temperature of IGBT modules [11].

The equivalent RC circuit of an IGBT module includes a partial fraction circuit and a continued fraction circuit. The partial fraction circuit, also known as the Foster model, obtains the dynamic thermal impedance curve through simulations or experiments, and the RC thermal characteristic parameters of each order are obtained by exponential fitting. However, the Foster model cannot reflect the temperature distribution inside the device [12], nor can it consider the actual situation where thermal resistance varies owing to internal factors [13]. The continued fraction circuit, also known as the Cauer model, reflects the actual values of thermal resistance and thermal capacitance of the physical layers and can predict the temperature of each layer. Reference [14] improved the accuracy of the Cauer model by subdividing it into layers. References [15,16] established an improved Cauer model that considered chip solder voids. Reference [17] proposed a method for estimating the thermal network parameters using IGBT module cooling curves to obtain the RC parameters of the Cauer model.

The most important step in solving the Cauer RC network is to obtain the heat diffusion angle. The thermal impedance of each layer can be precisely calculated only when the obtained heat diffusion angle is accurate. References [14,16,18,19] consider the heat diffusion angle as a fixed value to calculate the effective heat conduction area. References [20,21,22] consider the heat diffusion angle as the arctangent value of the ratio of the thermal conductivity of this layer to that of the next layer. However, the accuracy of the heat diffusion angle obtained using these two methods is insufficient. Reference [23] obtained the heat diffusion angle by obtaining the heat flux density curve using the finite element method (FEM), but this method is time-consuming.

IGBT modules typically have multiple chips, and when power is applied to two or more chips, there may be overlapping areas of heat flux diffusion, resulting in thermal coupling effects. However, the thermal coupling effect is usually not considered in traditional Cauer thermal network models, which can lead to lower junction temperatures. Some scholars have studied the thermal coupling effect. Reference [24] mainly analyzed the thermal coupling effect of the cooling system and believed that the case-to-ambient coupling thermal resistance was the main coupling thermal resistance. However, the coupling thermal resistance was derived from the junction temperature, and an essential formula for calculating the coupling thermal resistance was not provided. Reference [25] used electrical phenomena to obtain a thermal coupling model but had high requirements for the equipment. References [26,27] used the FEM to observe the thermal coupling effect, and reference [28] used the FEM to obtain the transient thermal impedance curve before establishing the thermal coupling effect. They took a long time to obtain the thermal coupling model using the FEM. Reference [29] obtained a thermal coupling model by establishing a discrete thermal network model; however, the extraction process of the Cauer model is complex. Reference [30] established a thermal coupling model that decouples internal thermal coupling from cooling system thermal coupling but does not determine the thermal coupling region.

However, the nonlinear characteristics of the TS of materials are usually ignored in existing Cauer thermal network models, which can result in the Cauer thermal network model not being in line with the actual situation of the IGBT modules, and the accuracy of the thermal network model deviates from the actual situation. As the junction temperature increases, the TS of the materials has a more significant impact on the actual temperature of the chip. References [23,31] considered the temperature-dependent characteristics of the chip and ceramic layers and established a thermal network model. Reference [32] considered the temperature dependence of packaging materials. However, these studies have a slow solving process owing to the use of the finite element method. Reference [33] obtained a transient thermal model using a finite element thermal simulation method with physical parameters varying with temperature; however, this study did not consider the temperature-dependent characteristics of the ceramic layer.

Table 1. Summary of characteristics of existing research.

	Method	Accuracy	Speed
Heat diffusion angle	Fixed value [14,16,18,19]	low	fast
	The inverse tangent value of the ratio of thermal conductivity [20,21,22]	low	fast
	Heat flux density curve by the FEM [23]	high	slow
	Proposed method	high	fast
	Reference	Shortcoming
Thermal coupling effects	[24,30]	thermal coupling region resistance is not provided
	[25,29]	difficult to implement
	[26,27,28]	slow calculation
	Reference	Accuracy	Speed
Temperature sensitivity	[23,31,32]	high	slow
	[33]	low	slow
	proposed method	high	fast

summarizes the characteristics of existing research.

Considering the existing research strategies and existing problems mentioned above, this article proposes a thermal impedance model for IGBT modules considering the nonlinear characteristics of the TS of materials using the Fourier series analysis method. First, Fourier Transform was used to analyze the heat flux density curve, and the heat diffusion angles were obtained from the heat flux density curve, which can be used to calculate the effective heat conduction area and obtain the thermal impedance of each layer. Then, considering the TS of the ceramic and chip layer materials, the thermal conductivities of the chip and ceramic layers were corrected to obtain the corrected thermal resistance and thermal capacitance of each layer. In addition, the thermal coupling effect was considered. Finally, the accuracy of the thermal impedance model established in this study is verified through a comparative analysis with the results of the FEM.

2. Single-Chip Thermal Impedance Model

2.1. Equivalent Circuit Model of the IGBT Module

The IGBT module consisted of a seven-layer structure [34,35]. The cross-section of a typical IGBT module is shown in Figure 1.

Assuming that the heat flux generated by the IGBT chip diffuses at a fixed angle in different layers, the area of the heat flux diffusion and the corresponding circuit structure are shown in Figure 2. In this system, $k_{chip}$,$k_i$ (i = 1,2,3,4,5,6) represents the thermal conductivity of each layer, $c_{chip}$,$c_i$ (i = 1,2,3,4,5,6) represents the specific heat capacity of each layer, $t_{chip}$,$t_i$ (i = 1,2,3,4,5,6) represents the thickness of each layer, ${\phi }_i$ represents the heat diffusion angle of each layer, $T_{chip},T_i$ represents the temperature of the upper surface of each layer vertically from the center of the chip downwards, $R_{chip},R_i$ represents the thermal resistance of each layer, $C_{chip},C_i$ represents the thermal capacitance of each layer, $T_a$ represents the temperature of the environment, and Q represents the heat flux, that is, the power loss generated by the IGBT chip. Consider a thin layer with a thickness $∆z$ in the heat flux diffusion area. The length of the j-1th thin layer is represented by $l_{j-1}$, the width is represented by $w_{j-1}$, and the effective heat conduction area is represented by $A_{j-1}$. The length of the jth thin layer is represented by $l_j$, the width is represented by $w_j$, and the effective heat conduction area is represented by $A_j$.

The effective heat conduction area $A_j$ of the jth thin layer is given by

$$\left\{\begin{array}{l}{l}_j=l_{j-1}+2\mathsf{\Delta }z\tan {\phi }_i,w_j=w_{j-1}+2\mathsf{\Delta }z\tan {\phi }_i\\ A_j=l_j{w}_j\\ l_0=c,w_0=d\end{array}\right.,$$

(1)

where c is the length of the IGBT chip, d is the width of the IGBT chip, j is the order of the thin layer with a thickness $∆z$ from top to bottom, and i (i = 1,2,3,4,5,6) is the order of different material layers from the chip solder layer to the baseplate.

The thermal resistance $∆R_j$ of the jth thin layer is calculated by [36]

$$\mathsf{\Delta }{R}_j={\displaystyle \frac{\mathsf{\Delta }z}{k_i{A}_j}}.$$

(2)

The thermal capacitance $∆C_j$ of the jth thin layer is calculated by [36]

$$\mathsf{\Delta }{C}_j=c_i{\rho }_i{A}_j\mathsf{\Delta }z,$$

(3)

where ${\rho }_i$ is the material density of each layer.

The convective heat transfer resistance is

$$R_{conv}={\displaystyle \frac{1}{h{A}_s}},$$

(4)

where h is the convective heat transfer coefficient and $A_s$ is the effective convective heat transfer area; h can first be calculated theoretically using empirical formulas, then verified through simulation, and finally calibrated through steady-state experiments.

The steady-state junction temperature $T_{chip}$ of the chip is calculated by [37]

$$T_{chip}=(R_{chip}+R_1+R_2+R_3+R_4+R_5+R_6+R_{conv})Q+T_a.$$

(5)

The transient junction temperature $T_{chip}\left(t\right)$ of the chip is calculated by

$$T_{chip}(t)=L^{-1}({\displaystyle \frac{1}{s}}{Z}_{th}(s)Q)+T_a,$$

(6)

where $L^{-1}$ is the Inverse Laplace Transform of L, and $Z_{th}\left(s\right)$ is the total thermal impedance of the Cauer model.

2.2. Solving the Heat Diffusion Angle of Each Layer in the IGBT Module

The heat generated by the IGBT chips does not transfer vertically downwards but undergoes lateral diffusion [38]. The lateral diffusion of heat leads to a smaller vertical heat flux density component $q_z$. Because the heat remains constant, the increased effective heat conduction area generates a heat diffusion angle. Therefore, the change in the heat flux density generates the heat diffusion angle, and there is a relationship between the heat flux density and heat diffusion angle. Heat diffusion angles can be obtained from the heat flux density curve.

The effective heat conduction area $A_j$ of the jth thin layer is calculated by [31]

$$A_j={\displaystyle \frac{Q}{q_z(z)}},$$

(7)

where $q_z\left(z\right)$ is the magnitude of the component of the heat flux vector in the z-axis. The effective heat conduction area can be approximately circular [23], and a schematic diagram of a thin layer with a thickness $∆z$ is shown in Figure 3; $r_{j,1}$ represents the radius of the upper surface of the thin layer, and $r_{j,2}$ represents the radius of the lower surface of the thin layer. The tangent value of the heat diffusion angle of the jth thin layer can be calculated by

$$\tan {\phi }_j={\displaystyle \frac{r_{j,2}-r_{j,1}}{\mathsf{\Delta }z}}.$$

(8)

Let r(z) be the curve of the effective heat conduction radius. According to Equation (8), the tangent of the heat diffusion angle is the slope of r(z), which is given by

$$r(z)=\sqrt{A(z)/\pi }.$$

(9)

The heat flux density component $q_z\left(z\right)$ and effective heat conduction radius r(z) from the center of the upper surface of the chip solder layer vertically downwards to the bottom surface of the baseplate are shown in Figure 4. By linearly fitting the r(z) curves corresponding to each layer separately, that is, $r=a_{1}z+a_0$, the tangent value of the heat diffusion angle of each layer is $a_1$.

2.3. Temperature-Sensitive Properties of Materials

The TS of each layer of material is shown in Figure 5. The thermal conductivity of Cu changes very little with temperature, indicating the low sensitivity of the thermal conductivity of copper to temperature. The thermal conductivities of Si and Al₂O₃ vary greatly with temperature, and the rate of change in thermal conductivity at 300 °C relative to that at 25 °C can even exceed 50%. The thermal conductivities of Si and Al₂O₃ are highly sensitive to temperature. As the junction temperature increased, the difference between the temperature of each layer material and the set temperature increased. The difference between the thermal conductivity and initial set thermal conductivity increased. The TS of the material has a more significant impact on the actual junction temperature. Therefore, it is necessary to consider the TS of materials with a high TS to make the calculated junction temperature more accurate. This study considered the TS of the ceramic and chip layers.

Although the datasheet of the devices does not provide material-related parameters, it does provide information about what the material of each layer is. After knowing what the material of each layer is, we can obtain the material-related parameters by searching the books and literature. The fitting formula for the relationship between the thermal conductivity and temperature of Si and Al₂O₃ is described as [23]

$$\left\{\begin{array}{l}k=438056{(T+273.15)}^{-1.4},\mathrm{Si}\\ k=40250{(T+273.15)}^{-1.264}{,\mathrm{Al}}_2{\mathrm{O}}_3\end{array}\right..$$

(10)

The fitting formula for the relationship between the specific heat capacity and the temperature of Si and Al₂O₃ is described as [23]

$$\left\{\begin{array}{l}c=6\ast {10}^{-6}\ast T^3-0.0044\ast {\mathrm{T}}^2+1.3946\ast \mathrm{T}+673.43,\mathrm{Si}\\ c=4\ast {10}^{-6}\ast T^3-0.0049\ast {\mathrm{T}}^2+2.4163\ast \mathrm{T}+693.29,{\mathrm{Al}}_2{\mathrm{O}}_3\end{array}\right..$$

(11)

When considering the TS of chips and ceramic materials, it is necessary to adjust the thermal conductivity and the thermal resistance of the ceramic and chip layers. The center-point temperature of the chip layer can be selected to represent the temperature of the chip layer, and the center-point temperature of the ceramic layer can be selected to represent the temperature of the ceramic material. When $k_{chip}$ and $k_3$ satisfy Equation (10) with the temperature at the center points of the chip and ceramic layers, respectively, $k_{chip}$ and $k_3$ are the corrected thermal conductivities that consider the TS of the chip and ceramic.

2.4. Solving Heat Flux Density Using Fourier Series Analysis Method

Figure 6 shows the structure of the multilayer IGBT modules. The heat flux generated by the chip layer, which is not shown in the figure, is used as the boundary condition. In the figure, c is the length of the chip on the x-axis, d is the length of the chip on the y-axis, and a and b are the lengths of the chip solder layer to the baseplate along the x- and y-axes, respectively. $(X_c,Y_c,0)$ is the coordinate of the center of the lower surface of the chip, and $T_f$ is the temperature of the cooling fluid.

The three-dimensional steady-state heat conduction differential equation is the Laplace Equation [39]:

$${\displaystyle \frac{{\partial }^2{\theta }_i}{\partial x^2}}+{\displaystyle \frac{{\partial }^2{\theta }_i}{\partial y^2}}+{\displaystyle \frac{{\partial }^2{\theta }_i}{\partial z^2}}=0,$$

(12)

where ${\theta }_i$ is the temperature rise of the i-th layer relative to the cooling fluid and ${\theta }_i=T_i-T_f$. The method of separating variables can be used to solve Equation (12) [40,41,42], and the general solution is

$$\begin{array}{l}{\theta }_i(x,y,z)=A_{i0}+B_{i0}z\\ +{\displaystyle \sum _{m=1}^{\infty }\cos ({\lambda }_{m}x)[A_{im}\cosh ({\lambda }_{m}z)+B_{im}\sinh ({\lambda }_{m}z)]}\\ +{\displaystyle \sum _{n=1}^{\infty }\cos ({\delta }_{n}y)[A_{in}\cosh ({\delta }_{n}z)+B_{in}\sinh ({\delta }_{n}z)]}\\ +{\displaystyle \sum _{m=1}^{\infty }{\displaystyle \sum _{n=1}^{\infty }\cos ({\lambda }_{m}x)\cos ({\delta }_{n}y)[A_{imn}\cosh ({\beta }_{mn}z)+B_{imn}\sinh ({\beta }_{mn}z)]}},\end{array}$$

(13)

where ${\lambda }_m$ and ${\delta }_n$ are the eigenvalues of the Fourier expansion in the x and y directions, respectively, ${\beta }_{mn}$ is the eigenvalue of the cross-coupled Fourier expansion; ${\lambda }_m=m\pi /a$, ${\delta }_n=n\pi /b$; and ${\beta }_{mn}=\sqrt{{\lambda }_m^2+{\delta }_n^2}$. $A_{i0}$ and $B_{i0}$ are the zero-order Fourier coefficients; $A_{im}$ and $B_{im}$ are the mth Fourier coefficients that represent the temperature change in the x-direction under constant z; $A_{in}$ and $B_{in}$ are the nth Fourier coefficients that represent the temperature change in the y-direction under constant z; and $A_{imn},B_{imn}$ is the Fourier cross-coupling coefficient. The heat flux density of each layer can be obtained by

$$\begin{array}{l}{q}_{zi}(x,y,z)=-k_i{\displaystyle \frac{d{\theta }_i(x,y,z)}{dz}}=-k_i(B_{i0}+{\displaystyle \sum _{m=1}^{\infty }{\lambda }_m\cos ({\lambda }_{m}x)[A_{im}\sinh ({\lambda }_{m}z)+B_{im}\cosh ({\lambda }_{m}z)]}\\ +{\displaystyle \sum _{n=1}^{\infty }{\delta }_n\cos ({\delta }_{n}y)[A_{in}\sinh ({\delta }_{n}z)+B_{in}\cosh ({\delta }_{n}z)]}\\ +{\displaystyle \sum _{m=1}^{\infty }{\displaystyle \sum _{n=1}^{\infty }{\beta }_{mn}\cos ({\lambda }_{m}x)\cos ({\delta }_{n}y)[A_{imn}\sinh ({\beta }_{mn}z)+B_{imn}\cosh ({\beta }_{mn}z)]}}).\end{array}$$

(14)

Assuming that the heat generated by the chip is uniformly applied to the lower surface of the chip and that the heat flux density is applied as a boundary condition to the upper surface of the chip solder, the four sides of the entire module and the upper surface of the chip solder layer that is not in contact with the chip are insulated, and the bottom surface of the baseplate is subjected to convective boundary conditions. When the contact between the layers is perfect, the heat flux density and temperature between the layers are continuous. Therefore, the boundary conditions for the IGBT modules can be written as [40]

$$\left\{\begin{array}{l}{\left.{\displaystyle \frac{\partial {\theta }_1}{\partial z}}\right|}_{z=0,contact\mathrm{area}}=-{\displaystyle \frac{Q}{k_{1}cd}},{\left.{\displaystyle \frac{\partial {\theta }_1}{\partial z}}\right|}_{z=0,\mathrm{not}\mathrm{contact}\mathrm{area}}=0\\ {\left.k_i{\displaystyle \frac{\partial {\theta }_i}{\partial z}}\right|}_{z=z_i}={\left.k_{i+1}{\displaystyle \frac{\partial {\theta }_{i+1}}{\partial z}}\right|}_{z=z_i},{\theta }_i(x,y,z_i)={\theta }_{i+1}(x,y,z_i)\\ {\left.k_6{\displaystyle \frac{\partial {\theta }_6}{\partial z}}\right|}_{z=z_6}={\left.-h{\theta }_6\right|}_{z=z_6}\end{array}\right.$$

(15)

where $z_i=\sum _{g=1}^i{t}_g$.

Based on the boundary conditions on the upper surface of the first layer, it can be concluded that [40]

$$\left\{\begin{array}{l}{B}_{1m}={\displaystyle \frac{-4Q\cos ({\lambda }_m{X}_c)\sin ({\lambda }_{m}c/2)}{abc{k}_1{{\lambda }_m}^2}}\\ B_{1n}={\displaystyle \frac{-4Q\cos ({\delta }_n{Y}_c)\sin ({\delta }_{n}d/2)}{abd{k}_1{{\delta }_n}^2}}\\ B_{1mn}={\displaystyle \frac{-16Q\cos ({\lambda }_m{X}_c)\sin ({\lambda }_{m}c/2)\cos ({\delta }_n{Y}_c)\sin ({\delta }_{n}d/2)}{abcd{k}_1{\lambda }_m{\delta }_n{\beta }_{mn}}}\\ B_{10}={\displaystyle \frac{-Q}{k_{1}ab}}\end{array}\right.$$

(16)

Based on the bottom boundary conditions of the module, it can be concluded that [40]

$$B_{6p}(\xi )=-{\sigma }_6(\xi )A_{6p}$$

(17)

where $\xi$ can be replaced by ${\lambda }_m{,\delta }_n{,\beta }_{mn}$, respectively, and the subscripts p of the corresponding variables are m, n, and mn, respectively. Then, ${\sigma }_6$ is given by [40]

$${\sigma }_6={\displaystyle \frac{\xi \sinh (\xi z_6)+h/k_6\cosh (\xi z_6)}{\xi \cosh (\xi z_6)+h/k_6\sinh (\xi z_6)}}.$$

(18)

Based on the boundary conditions between each layer, it can be concluded that [40]

$${\sigma }_i={\displaystyle \frac{-B_{ip}}{A_{ip}}}={\displaystyle \frac{{\sigma }_{i+1}{\varphi }_i-{\psi }_i}{{\sigma }_{i+1}{\psi }_i-{\rho }_i}},$$

(19)

where

$$\left\{\begin{array}{l}{\varphi }_i={\cosh }^2(\xi {\displaystyle \sum _{j=1}^i{t}_j})-k_i/k_{i+1}{\sinh }^2(\xi {\displaystyle \sum _{j=1}^i{t}_j})\\ {\psi }_i=\sinh (\xi {\displaystyle \sum _{j=1}^i{t}_j})\cosh (\xi {\displaystyle \sum _{j=1}^i{t}_j})(1-k_i/k_{i+1})\\ {\rho }_i={\sinh }^2(\xi {\displaystyle \sum _{j=1}^i{t}_j})-k_i/k_{i+1}{\cosh }^2(\xi {\displaystyle \sum _{j=1}^i{t}_j})\end{array}\right..$$

Once ${\sigma }_6\left(\xi \right)$ is calculated, further calculations can be performed to obtain ${\sigma }_5\left(\xi \right)$, ${\sigma }_4\left(\xi \right)$, ${\sigma }_3\left(\xi \right)$, ${\sigma }_2\left(\xi \right)$, and ${\sigma }_1\left(\xi \right)$. By combining Equations (16) and (19), $A_{1m},A_{1n}$, and $A_{1mn}$ can be obtained. And $A_{ip},B_{ip}(i=\mathrm{2,3},\mathrm{4,5},6)$ can be sequentially calculated from [40]

$$\left\{\begin{array}{l}{A}_{(i+1)p}=A_{ip}{\varphi }_i+B_{ip}{\psi }_i\\ B_{(i+1)p}=-(A_{ip}{\psi }_i+B_{ip}{\rho }_i)\end{array}\right..$$

(20)

The zero-order Fourier coefficients can be calculated by [40]

$$\left\{\begin{array}{l}{B}_{i0}=-{\displaystyle \frac{Q}{k_{i}ab}}\\ A_{60}={\displaystyle \frac{Q}{ab}}({\displaystyle \frac{1}{h}}+{\displaystyle \frac{z_6}{k_6}})\\ A_{i0}={\displaystyle \frac{Q}{ab}}({\displaystyle \frac{1}{h}}+{\displaystyle \sum _{j=i+1}^{6-i}\frac{t_j}{k_j}}+{\displaystyle \frac{z_i}{k_i}}),i=1,2,3,4,5\end{array}\right..$$

(21)

3. Thermal Resistance Model Considering Multi-Chip Thermal Coupling

3.1. Improved Calculation Method of Coupling Thermal Resistance

Thermal coupling effects must be considered in IGBT modules with multiple chips. When heat is generated by multiple chips, there may be overlapping areas of heat diffusion among the chips, resulting in thermal coupling effects (TCEs). The effective heat conduction area is equivalent to a circle to facilitate the calculation. Figure 7 shows the thermal coupling region. The overlapping region is the thermal coupling region, and the corresponding thermal resistance is the thermal coupling resistance. The thermal coupling resistance increases the total thermal resistance, thereby increasing the junction temperature. $D_{pq}$ is the distance between the center of chip p and the center of chip q.

Figure 8 shows the thermal coupling x-y plane cross-section of the jth thin layer: $r_p$ represents the effective heat conduction radius of chip p in the jth thin layer, and $r_q$ represents the effective heat conduction radius of chip q in the jth thin layer.

The cosine theorem can be used to obtain the angles, and ${\alpha }_p$ and ${\alpha }_q$ are given by

$$\left\{\begin{array}{l}{\alpha }_p={\cos }^{-1}({\displaystyle \frac{r_p^2(z)+D_{pq}^2-r_q^2(z)}{2r_p(z)D_{pq}}})\\ {\alpha }_q={\cos }^{-1}({\displaystyle \frac{r_q^2(z)+D_{pq}^2-r_p^2(z)}{2r_q(z)D_{pq}}})\end{array}\right..$$

(22)

The thermal coupling area can be calculated by

$$\begin{array}{c}{A}_{pq,j}=({\stackrel{⌢}{A}}_{O_{p}AB}-A_{O_{p}AB})+({\stackrel{⌢}{A}}_{O_{q}AB}-A_{O_{q}AB})\\ =[{\displaystyle \frac{2{\alpha }_p}{2\pi }}\pi {r_p}^2(z)-{\displaystyle \frac{1}{2}}{r_p}^2(z)\sin (2{\alpha }_p)]\\ +[{\displaystyle \frac{2{\alpha }_q}{2\pi }}\pi {r_q}^2(z)-{\displaystyle \frac{1}{2}}{r_q}^2(z)\sin (2{\alpha }_q)].\end{array}$$

(23)

The effective heat conduction area of the jth thin layer can be obtained using Equation (7), and the effective heat conduction radius curve can be obtained using Equation (9). When the distance between two chips is too large and satisfies ${(r}_p\left(z_6\right)+r_q\left(z_6\right))\le D_{pq}$, there is no heat diffusion overlap, and no TCE occurs. When the distance between the two chips is relatively close and satisfies ${(r}_p\left(z_6\right)+r_q\left(z_6\right))>D_{pq}$, a heat flux diffusion overlap area will occur, resulting in TCEs.

As shown in Figure 8, the starting layer of the thermal coupling region is the layer where the sum of the effective heat conduction radius of chip p and the effective heat conduction radius of chip q is equal to the distance between the center points of the two chips, that is, when $r_{p,j-1}+r_{q,j-1}\le D_{pq},r_{p,j}+r_{q,j}>D_{pq}$ is satisfied, and the jth thin layer is the starting layer of the thermal coupling region.

Compared with the methods in [43], the proposed method can directly obtain the effective heat conduction radius of the two chips from the heat flux density curve. The starting position of the thermal coupling region can be determined based on the relationship between the effective heat conduction radius of the two chips and the distance between them. Compared to reference [43], the proposed method has a simpler and clearer calculation process.

The coupling thermal resistance of chip q to chip p includes the coupling heat conduction resistance and coupling convective heat transfer resistance. The coupling heat conduction resistance of chip q to chip p in the jth thin layer can be calculated by [43]

$$\mathsf{\Delta }{R}_{pq,j}={\displaystyle \frac{A_{pq,j}}{A_{q,j}}}\mathsf{\Delta }{R}_q,$$

(24)

where $A_{pq,j}$ represents the thermal coupling area of the jth thin layer, $A_{q,j}$ represents the effective heat conduction area of the jth thin layer of chip q, and $∆R_q$ represents the thermal resistance of the jth thin layer of chip q.

The coupling convective heat transfer resistance of chip q to chip p is given by

$$R_{pq\_conv}={\displaystyle \frac{A_{pq}(z_6)}{A_q(z_6)}}{R}_{q\_conv},$$

(25)

where $A_{pq}\left(z_6\right)$ represents the coupled convective heat transfer area of chip q to chip p, $A_q\left(z_6\right)$ represents the effective convective heat transfer area of chip q, and $∆R_{q\_conv}$ represents the convective heat transfer resistance of chip q.

3.2. Program Implementation Steps for Calculating Multi-Chip Junction Temperature

The improved calculation process for the junction temperature considering the thermal characteristics of the chips and ceramic materials is described as follows:

Step 1: The boundary conditions, material parameters, and size parameters of the IGBT module are entered.

Step 2: The Fourier coefficients are solved using Equations (16)–(21), and the path perpendicular to the center of the chip downwards is taken. The heat flux density curve $q_z\left(z\right)$ from the upper surface of the chip solder layer to the bottom of the baseplate is obtained using Equation (14).

Step 3: The effective heat conduction area curve A(z) from the upper surface of the chip solder layer to the bottom of the baseplate is calculated using Equation (7), and the effective heat conduction radius curve r(z) is calculated using Equation (9). Linear fitting is used to obtain the tangent values of the heat diffusion angle from the chip solder layer to the baseplate.

Step 4: The thermal resistance $R_i(i=\mathrm{1,2},\mathrm{3,4},\mathrm{5,6})$ and thermal capacitance $C_i(i=\mathrm{1,2},\mathrm{3,4},\mathrm{5,6})$ of each layer are derived using Equations (1)-(3). Assuming that the effective heat conduction area of the chip layer is the size area of the chip, the thermal resistance $R_{chip}$ and thermal capacitance $C_{chip}$ of the chip layer are obtained from Equations (2) and (3), respectively, and the convective heat transfer resistance $R_{conv}$ is obtained from Equation (4).

Step 5: Given the initial thermal conductivities of the chip and ceramic, the temperatures $T_{chip\_con}$ and $T_{3\_con}$ corresponding to $k_{chip}$ and $k_3$ are calculated with Equation (10). $T_{chip\_ct}$ and $T_{3\_ct}$ are the temperatures at the center of the chip and ceramic layers, respectively. $\left|T_{chip\_con}-T_{chip\_ct}\right|$, $\left|T_{3\_con}-T_{3\_ct}\right|$ are calculated, and when the error does not meet the requirements, new values are assigned to $k_{chip}$ and $k_3$.

Step 6: Steps 2, 3, 4, and 5 are repeated until $\left|T_{chip\_con}-T_{chip\_ct}\right|<{\epsilon }_1\&\left|T_{3\_con}-T_{3\_ct}\right|<{\epsilon }_2$, and then the iteration stops. At this point, the obtained $k_{chip}$ and $k_3$ are matched with the corresponding material temperature, and a thermal impedance model considering the TS of the material is obtained.

Step 7: The starting layer of the thermal coupling region is determined, the thermal coupling area is obtained from Equation (23), the thermal coupling resistance is obtained from Equations (24) and (25), and the junction temperature is obtained considering TCEs.

A program flowchart for calculating the multichip junction temperature is shown in Figure 9.

4. Simulation Verification

To validate the proposed method, the SKM50GB12T4 IGBT module was selected, and its ANSYS16.0 simulation model was built for steady-state and transient thermal simulations. Figure 10 shows the simulation model of the constructed single chip. Unlike the actual IGBT module model, the upper copper layer in Figure 10 does not have a pattern, and the area of the solder layer is larger than that of the chip. Although the simulation model of the IGBT module has been simplified, it is known from reference [40] that the simplified structure matches the junction temperature obtained from the actual structure very well. The parameters of each layer of the IGBT module, including the size of the IGBT module, thickness, specific heat capacity, thermal conductivity, and position of the chip, are listed in

Table 2. IGBT module parameters.

Parameter	Value	Parameter	Value
a	30.3 mm	$c_2$$,c_4$$,c_6$	$384\mathrm{J}/(\mathrm{k}\mathrm{g}·°\mathrm{C})$
b	28 mm	$c_3$	$\left(11\right)\mathrm{J}/(\mathrm{k}\mathrm{g}·°\mathrm{C})$
c	7.2 mm	$k_{chip}$	$\left(10\right)\mathrm{W}/(\mathrm{m}·°\mathrm{C})$
d	6.75 mm	$k_1$$,k_5$	$54\mathrm{W}/(\mathrm{m}·°\mathrm{C})$
$t_{chip}$	0.15 mm	$k_2$$,k_4$$,k_6$	$390\mathrm{W}/(\mathrm{m}·°\mathrm{C})$
$t_1$	0.1 mm	$k_3$	$\left(10\right)\mathrm{W}/(\mathrm{m}·°\mathrm{C})$
$t_2$$,t_4$	0.3 mm	$X_c$	11.25 mm
$t_3$	0.38 mm	$Y_c$	14.485 mm
$t_5$	0.12 mm	$X_{c2}$	25 mm
$t_6$	2.8 mm	$Y_{c2}$	22 mm
$c_{chip}$	$\left(11\right)\mathrm{J}/(\mathrm{kg}·°\mathrm{C})$	$T_a$	25 °C
$c_1$$,c_5$	$230\mathrm{J}/(\mathrm{kg}·°\mathrm{C})$

.

The results obtained from simulation can be used to calculate the heat conduction resistance $R_{sim}$ and convective heat transfer resistance $R_{conv\_sim}$:

$$\left\{\begin{array}{l}{R}_{sim}={\displaystyle \frac{T_{chip\_sim}-T_a}{Q}}\\ R_{conv\_sim}={\displaystyle \frac{T_c-T_a}{Q}}\end{array}\right.$$

(26)

where $T_{chip\_sim}$ is the junction temperature obtained through the simulation and $T_c$ is the case temperature.

Figure 11 shows the thermal resistance of the chip and ceramic layers, and the tangent of the heat diffusion angle of the ceramic layer is obtained by the FEM considering the TS, the proposed method considering the TS, and the method without considering the TS at different power losses. Because device suppliers typically do not provide information on active and passive areas, it can be assumed that the heat is evenly distributed within the chip and is transferred vertically downwards with a heat diffusion angle of 0°. The influence of temperature on the tangent value of the heat diffusion angle of the chip layer is not shown.

As shown in Figure 11, when the value of power loss is low, the results obtained by the three methods are relatively close. However, as the power loss increases, the difference between the results obtained by the method without considering TS and the FEM becomes increasingly significant. In contrast, the result obtained by the proposed method is still very close to that of the FEM, and their error is very small. This indicates that as the power loss increases, the TS of the material has an increasingly significant impact on the heat diffusion angle of the ceramic layer and the thermal resistance of the chip and ceramic layers. This is because as the material temperature increases, the values of thermal conductivity of the ceramic and chip layers deviate more from those at the predetermined temperature. The decrease in the thermal conductivity results in an increase in the thermal resistance and a decrease in the tangent value of the heat diffusion angle.

Under different boundary conditions of power loss, the junction temperatures obtained by the FEM considering the TS, the proposed method considering the TS, and the method without considering the TS are shown in Figure 12. This shows that as the power loss increases, the error of the method without considering TS compared to the FEM considering TS increases, even exceeding 9% when the power loss is 200 W. When the power loss is less than 200 W, the maximum error of the method proposed in this study does not exceed 4%, which is more than 5% lower than that of the method that does not consider the TS.

This indicates that considering the TS can significantly improve the accuracy of the junction temperature. The results obtained using the proposed method are extremely close to those obtained using the FEM, which verifies that the proposed method is effective and accurate.

The proposed method calculates a lower junction temperature than the FEM because using the midpoint temperature of the ceramic layer and chip layer to represent the temperature of the entire layer is not sufficiently accurate, resulting in a higher thermal conductivity, a larger heat diffusion angle, and a lower thermal resistance. This leads to a lower junction temperature.

To increase the accuracy of the calculated junction temperature, the chip and ceramic layers can be refined into multiple thin layers. The temperature of each thin layer and the thermal conductivity of each thin layer can be determined iteratively. At this point, the temperature and thermal conductivity of each point in the chip and ceramic layers match the fitting equation very well. Although this method can improve accuracy, it introduces an iterative calculation process to determine the thermal conductivity of each thin layer, which significantly reduces the calculation speed. Owing to the poor linearity of the effective heat conduction radius curve of the baseplate, the baseplate can be divided into multiple layers, and linear fitting can be performed on each layer separately. This will result in a more accurate heat diffusion angle of the baseplate, thereby improving the accuracy of the thermal impedance model.

When the power loss is 120 W, under different boundary conditions of convective heat transfer, the junction temperatures derived from the FEM considering TS, the proposed method considering TS, and the method without considering TS are shown in Figure 13. Under different convective heat transfer conditions, the junction temperature obtained by the proposed method considering TS is very close to that obtained by the FEM, whereas the junction temperature obtained by the method without considering TS is significantly different from that obtained by the FEM. This verifies that the proposed method is effective under different boundary conditions of convective heat transfer.

Figure 14 shows the transient junction temperatures obtained by the FEM considering the TS, the proposed method considering the TS, and the method without considering the TS when the power loss is 120 W. Compared with the FEM, the proposed method obtains a faster increase in the junction temperature. Although the junction temperature rise process is somewhat different, the final steady-state junction temperature obtained is almost the same. Because there are no spikes during the process of junction temperature rise, the difference in the temperature rise of each junction has little effect on the reliability analysis of the chip. However, the steady-state junction temperature obtained by the method without considering TS differs significantly from that derived using the FEM.

By increasing the values of the thermal capacitance of each layer obtained using the proposed method, the transient junction temperatures compared with the FEM were obtained when the values of the thermal capacitance of each layer obtained using the proposed method were increased by 10%, 30%, and 50%, as shown in Figure 15. The larger the proportion of increase in the value of thermal capacitance, the slower the rate of junction temperature rise, and the closer the junction temperature curve obtained by the proposed method is to that obtained by the FEM. This indicates that the thermal capacitance of each layer is a key factor that affects the rate of increase of the junction temperature. The larger the value of the thermal capacitance in each layer, the slower the rise in the junction temperature and the longer it takes to reach a steady state.

Using the same computer, the CPU model is Intel (R) Xeon (R) Gold6426Y, which is a 32-core processor with 64 GB of running memory. The number of nodes in the FEM is 140,000, and the solution time is as long as 7.5 min. However, the solution time of the proposed method is 51 s, which reduces the solution time by 88.7% and significantly accelerates the solution speed.

When the power loss is applied to chip p, Equation (26) can be used to obtain the heat conduction resistance $R_{p\_sim}$ and convective heat transfer resistance $R_{p\_conv\_sim}$ of chip p. When power loss is applied to chips p and q, Equation (26) can be used to obtain the heat conduction resistance $R_{p\_sim2}$ and convective heat transfer resistance $R_{p\_conv\_sim2}$ of chip p. The coupling heat conduction resistance $R_{pq\_sim}$ and coupling convective heat transfer resistance $R_{pq\_conv\_sim}$ of chip q to chip p through simulation are given by

$$\left\{\begin{array}{l}{R}_{pq\_sim}=R_{p\_sim2}-R_{p\_sim}\\ R_{pq\_conv\_sim}=R_{p\_conv\_sim2}-R_{p\_conv\_sim}\end{array}\right.$$

(27)

Based on considering the TS of a single chip and considering the TCE between chips, Figure 16 shows the multi-chip thermal simulation model. The power loss of both chips was 120 W, and Table 2 lists the coordinates of the two chips. It can be seen from Figure 16 that the junction temperature of one of the chips was 139.51 °C.

Table 3. Junction temperatures obtained by different methods.

Method	Junction Temperature (°C)	Error (%)
FEM	139.51	0
Proposed method considering TS	111.15	20.33
Proposed method considering TS and TCE	128.72	7.73

shows the junction temperatures obtained using the different methods. It shows that based on considering the TS mentioned above, considering the TS and TCEs between chips can reduce the error by 12.6%. Considering the TCEs can greatly enhance accuracy. The error of the proposed method considering the TS and TCEs is only 7.73%, which verifies that the proposed method is highly accurate. The starting layer of the thermal coupling region was located in the baseplate in this simulation. The coupling heat conduction resistance calculated using the proposed method was 0.00083 (K/W), and the coupling convective heat transfer resistance was 0.146 (K/W). This indicates that the coupling convective heat transfer resistance contributes significantly to the junction temperature, whereas the contribution of the coupling heat conduction resistance to the junction temperature can be ignored. The reason for the lower junction temperature calculated using the proposed method may be the simplification of the shape of the effective heat conduction area. Actually, the shape of the effective heat conduction area is similar to a rectangle, with curvature at all four corners. This study simplified the shape of the effective heat conduction area into a circle to facilitate the calculation. This leads to a smaller calculated coupled heat conduction area, resulting in a lower calculated coupled thermal resistance.

5. Conclusions

This study proposes a thermal impedance model that considers the nonlinear thermal characteristics of the TS of chips and ceramic materials. Using the Fourier series method to analyze the heat flux density, the effective heat conduction radius curve was obtained, and then the heat diffusion angle was obtained through linear fitting. Considering the nonlinear thermal characteristics of the material, the modified values of the thermal resistance and thermal capacitance were obtained, thus obtaining a single-chip thermal impedance model that considers the nonlinear thermal characteristics of the material. On this basis, a thermal resistance model for multiple chips was obtained considering the TCEs of multiple chips.

To verify the reliability of the method proposed in this study, a finite element model was constructed and simulated. Through simulation verification, the following can be taken to be the conclusions of this study:

• When the power loss was less than 200 W, compared with the FEM, the maximum errors in the thermal resistance of the ceramic layer, the thermal resistance of the chip layer, and the tangent of the heat diffusion angle of the ceramic layer obtained by the method proposed in this study were 6.73%, 1%, and 5.5%, respectively. However, compared with the FEM, the maximum errors in the thermal resistance of the ceramic layer, thermal resistance of the chip layer, and tangent of the heat diffusion angle of the ceramic layer obtained by the method without considering the TS of the material were 30%, 25.1%, and 17.34%, respectively. After considering the TS of the ceramic and chip materials, the accuracy of the thermal resistance of the ceramic and chip layers and the accuracy of the tangent value of the ceramic layer heat diffusion angle were greatly improved.
• When the power loss is less than 200 W, the error of the method without considering TS compared to the FEM considering TS increases, even exceeding 9% when the power loss is 200 W. When the power loss is less than 200 W, the maximum error of the method proposed in this study does not exceed 4%, which is more than 5% lower than that of the method that does not consider the TS. This indicates that the nonlinear thermal characteristics of the chip and ceramic materials can affect the heat diffusion angle and thermal resistance of the chip and ceramic layers, thereby affecting the junction temperature. Considering the TS of the materials can significantly improve the accuracy of the junction temperature.
• Compared to the FEM, the junction temperature obtained by the proposed method increased faster. Although there are some differences in the process of junction temperature rise, the final steady-state junction temperature obtained is almost the same. The rise in junction temperature becomes slower by increasing the thermal capacitance of each layer, indicating that thermal capacitance affects the speed of the junction temperature rise.
• Based on considering the TS mentioned above, considering the TS and TCEs between chips can reduce the error by 12.6%, and the error of the proposed method considering TS and TCEs is only 7.73%. This indicates that the TCEs between chips will greatly increase the junction temperature of the chips, and considering the TCEs can greatly improve the accuracy.
• Compared to the FEM, the solution time of the proposed method was reduced by 88.7%. In addition, the FEM has a complex modeling process, whereas the proposed method requires only the input of the size, material, and boundary parameters of the IGBT modules.

Compared to [23,31], the proposed method has a similar accuracy but a much faster solving speed. Compared to reference [33], the proposed method not only has a higher accuracy but also has a much faster solving speed. Compared to references [25,29], the proposed method not only has high accuracy but is also easier to implement.

Overall, the method proposed in this article not only has high computational accuracy but also has a fast solving speed and simple implementation process, making it highly practical.

A future research agenda is to conduct experimental verification, which can better validate the method proposed in this study.