李继生+王勐+黄战武

摘 要： 随着电子产品的尺寸和重量日趋变小，经典Fourier理论已经不能很好地解决实际问题。为了使理论结果尽可能准确地反映实际情况，需要用到Non?Fourier导热模型。首先，建立了三维多芯片组件的Fourier和Non?Fourier导热模型；其次，采用有限差分方法求解相应模型的传热方程，得到了温度分布和温度响应；最后，使用有限元热分析软件ICEPAK建立了相应的热分析模型并进行计算。实验结果表明，与经典的Fourier模型相比，Non?Fourier模型更加接近实际温度，且温度场进入稳态的时间较长，热耦合的现象也更强。

关键词： 多芯片组件； Fourier； Non?Fourier； 有限差分法

中图分类号： TN710?34； TP311 文献标识码： A 文章编号： 1004?373X（2014）03?0117?05

Thermal analysis and research of multi?chip component substrate

based on Non?Fourier thermal conductivity model

LI Ji?sheng 1， WANG Meng2， HUANG Zhan?wu2

（1. Shaanxi Changling Photovoltaic Electric Co.， Ltd.， Baoji 721006， China； 2. Institute of Electronic CAD， Xidian University， Xian 710071， China）

Abstract： With the reduction in size and weight of electronic products， classical Fourier theory can′t solve practical problems well. To make the theoretical results close to the fact， Non?Fourier model should be used. Firstly， the Fourier model and the Non?Fourier model for three?dimension multi?chip module （MCM） are built separately. Secondly， finite difference method （FDM） is used to solve the corresponding equation of heat transfer and get distribution and response of temperature. Finally， the relevant thermoanalysis model is built up by finite element analysis software of ICEPAK and calculated. The result shows that the temperature of Non?Fourier model is closer to the practical situation than classical Fourier model. The time is longer for temperature field to enter the stable state， and the thermo?coupled phenomenon is stronger.

Keywords： MCM； Fourier； Non?Fourier； FDM

0 引 言

随着便携式电子产品小型化、多功能化和环境保护的要求，片上系统（System on Chip，SoC）等高密度三维立体（3D）堆叠封装技术得到了迅速发展，这种技术在减小芯片尺寸和实现紧密集成度等方面可以满足目前电子行业的要求[1?2]。

现代电子技术迅猛发展，要求电子整机朝着短、小、轻、薄及高可靠、高速、高性能和低成本的方向发展[3]。热失效是电子设备的主要失效形式，其失效率随温度增加将呈指数增长趋势[4]。在设计阶段，热分析非常关键，通过其可以发现电子设备的热缺陷[5?6]。

当前电子设备的集成度不断提高，器件的特征尺寸也愈来愈小。当尺度（空间尺度和时间尺度）微细化后，将会出现传热的尺度效应，此时的传热情况与常规尺度下不同。对于多芯片组件，更多的使用经典Fourier理论去分析解决问题，但尺度微细化后的热量传播包含了不同于扩散的波动传播机制，这是以经典Fourier定律为基础而建立起来的瞬态传导理论不能合理解释的。但目前关于电子封装的Non?Fourier效应研究较少，相关的数据资料也比较有限。本文尝试以三维多芯片组件为模型，分别用Fourier和Non?Fourier分析方法研究其热效应问题，并对实验结果进行了分析对比。

1 多芯片组件导热模型

1.1 多芯片组件模型

本文对三维多芯片组件进行分析，在直角坐标系中，建立如图1所示的三维多芯片组件模型。基板模型呈长方体状，其几何尺寸为：[xa×yb×zc=70×70×5 mm3][Z=0]表面标有编号1～6的位置布置有6个芯片，在此忽略芯片的实际尺寸，只在基板表面有芯片的位置考虑有若干热源的情况。在热源位置内，热源的热耗散均匀分布，与坐标位置和时间均无关。

图1 三维多芯片组件模型

1.2 热传导问题的Fourier分析

在《热的分析理论》中[7]，法国科学家约瑟夫·傅里叶（Joseph Fourier）结合实验观察，阐述了把温度梯度和热流联系起来的基本定律：对于均匀的各向同性固体，傅里叶定律可表示为：

[q（r，t）=-k?T（r，t）] （1）

式中：热流密度[q（r，t）]表示单位等温面上，在温度降低方向上单位时间内的热流量（单位：W·m2），[k]为材料的导热系数（单位：W/（m·℃）），[T（r，t）]为温度分布（单位：℃）。由于[q（r，t）]指向温度降低方向，式（1）中的负号就使热流量成为一个正的量。在直角坐标系中，式（1）被改写为：

[q（x，y，z，t）=-ik?T?x-jk?T?y-kk?T?z] （2）

式中[i，j，k]分别为沿[x，y，z]方向的单位向量。

对一个很小的控制体[V]建立热量平衡方程，可描述如下[8]：

[单位时间内通过V的边界进入的热量+单位时间内V内产生的热量=单位时间内V内能量的累积]

对于上式中的每一项，以其各自表达式表示，可得式（3）：

[-Aq?ndA+-Vg（r，t）dV=-VρCp?T（r，t）?tdV] （3）

对式（3）进行转化，可得静止均匀物体内含有热源的各向同性物质的热传导方程，即：

[?[k?T（r，t）]+g（r，t）=ρCp?T（r，t）?t] （4）

1.3 热传导问题的Non?Fourier分析

随着工程技术的快速发展，热传导领域的瞬态过程日益受到人们的重视[9?10]。对于稳定温度分布的无限大介质，突然出现的热扰动会对介质的稳定状态造成破坏，历经松弛过程之后，介质将处于新的稳定状态。在温度场的重建过程中，松弛时间以[τ0]表示，热量传播速度以[Ch]表示，从热扰动的所在点算起，热量传播深度为[δh=Chτ0。]由于热扩散率的量纲m2/s=（m/s）·m，因此可得[δh=αCh，]则热量传播速度的表达式（以松弛时间和热扩散率表示）如下：

[Ch=ατ0=kρCτ0] （5）

修正经典的Fourier定律，可得：

[q=-kgradt-τ0?q?t] （6）

其中，[?q?t]为在温度梯度截面上热流密度对时间的变化率，[τ0?q?t]为松弛时间间隔内截面上热流密度的改变量。

当热量传播速度为有限值时，整理后的热传导微分方程如下：

[τ0?2T?t2+?T?t-τ0ρC?Q?t=α?2T?x2+?2T?y2+?2T?z2+QρC] （7）

若内热源[Q]与时间无关，式（7）可简化为：

[τ0?2T?t2+?T?t=α?2T?x2+?2T?y2+?2T?z2+QρC] （8）

2 Non?Fourier差分方程的求解

相比于Fourier导热微分方程，Non?Fourier导热微分方程的求解更为困难。对于Non?Fourier导热微分方程，仍然使用交替方向隐式法，得到如下差分格式：

3 温度场计算及结果

首先建立Fourier和Non?Fourier导热微分方程相应的差分格式，然后通过求解得到的关于网格节点温度的方程组，即可求得所有节点温度值。

图2分别是Fourier和Non?Fourier分析中，顶面最高温度随时间变化的温度曲线图。在开始的100 s内，Fourier和Non?Fourier分析对应的顶面最高温度上升趋势相当。100 s以后，Fourier分析对应的温度上升速度逐步放缓，并迅速进入到温度场稳定状态，在600 s后稳定在大约70 ℃。Non?Fourier分析对应的温度曲线则比较陡峭，在100 s之后温度上升的速度仍然比较快，在1 200 s左右才开始逐渐进入温度场的稳定状态，稳定之后温度大约为105 ℃。在整个过程中，Non?Fourier分析对应的顶面最高温度要高于Fourier分析，两种分析相应的温度差也不断增大。Fourier分析可以看成Non?Fourier分析的特殊情况，即松弛时间为零的情况，在相同条件下，其顶面最高温度的相应曲线落在图2两条曲线之间。

图2 顶面最高温度曲线图

图3为顶面温度场的相应最大温度差随时间变化的曲线图。短暂的时间过后，两种分析的最大温度差都迅速稳定下来，远快于整个温度场达到稳定状态所用的时间。相比于Fourier分析， Non?Fourier分析的最大温度差达到稳定所需的时间要略长一些。

图2和图3出现的这些差别，是由于Fourier分析以热量传播速度无限大为前提，同时热量被整个物体的全部质点不同程度地吸收；而Non?Fourier分析考虑了热量的传播速度，进入某一截面的热流一部分以热传导的方式向物体内部传导，而另一部分用在了该截面分子松弛过程所导致的热量改变。

图3 顶面最大温差曲线图

应用有限元热分析软件Icepak建立三维多芯片组件模型，如图4所示。

图4 Icepak模型图

由图4可知，软件仿真和通过计算所得的温度场分布大致相同，但Non?Fourier模型更加接近。

分析在某些关键点和路径上的温度分布状况，在基板上选取下述关键点和路径。

路径1：[y=21，z=0，x∈（0，70）]

路径2：[y=35，z=0，x∈（0，70）]

路径3：[y=56，z=0，x∈（0，70）]

路径4：[x=42，y=42，z∈（0，5）]

路径5：[x=21，y=21，z∈（0，5）]

路径6：[x=14，y=56，z∈（0，5）]

路径7：[x=y，z=0，x∈（0，70），y∈（0，70）]

路径8：[x=70-y，z=0，x∈（0，70），y∈（0，70）]

图5分别为在Fourier和Non?Fourier分析两种情况下，沿[x]方向路径1、2、3的温度曲线图。由图5可知，除温度值不同，及热源6位置与其附近区域的温度值在曲线上略有不同外，两种分析情况中的路径曲线走势基本是一致的。

图5 路径1、2、3的温度曲线图（1 300 s）

图6分别为在Fourier分析和Non?Fourier分析两种情况下，沿[z]方向路径4、5、6的温度曲线图。由图6可知，这三条路径上的温度从顶面到底面逐渐降低，且温度变化不大，变化规律基本一致。这是因为沿[z]方向的尺寸很小，并且仅有底面对外传热。

图7的路径7、8是顶面上的两条对角线。路径8的温度曲线左右对称，越靠近路径8的中点，温度越高。整个温度场的分布关于路径7对称。这是因为基板表面布置的热源在也是以路径7为对称轴对称分布的，而且发热量在对称轴附近是最大的。因而热源的分布状况对温度场的影响很大，所以基板上芯片与其他热源的位置分布就显得尤为重要，要尽量分散热源，避免过于集中而导致局部温度过高。

图6 路径4、5、6的温度曲线图（1 300 s）

图7 路径7、8的温度曲线图（1 300 s）

路径7贯穿热源1和热源6所在的位置，Fourier分析中路径8的温度曲线中部比Non?Fourier分析的曲线变化更为剧烈，Fourier分析中曲线的向上凸起程度不及Non?Fourier分析。同时，与热源位置处的温度相比，热源1和6之间相同区域的温度都比较低，而Non?Fourier分析中没有出现这一情况。

由此可见，Non?Fourier分析中的热耦合现象明显强于Fourier分析。

对温度场进行分析可得，使用Fourier分析和Non?Fourier分析得到的结果存在着很多差异。Non?Fourier导热模型的温度上升速度快于Fourier分析，温度曲线比Fourier导热模型的更为陡峭，进入稳定状态的时间远迟于Fourier分析且进入稳定状态后温度值较高，Non?Fourier模型的温度场最大温差要大于Fourier模型，热耦合现象也较强。

4 结 语

本文分别建立了Fourier分析和Non?Fourier分析导热微分方程的有限差分格式，通过求解相应代数方程组，得到两种分析情况下的温度场分布；分析比较了使用Fourier和Non?Fourier分析对基板温度场分布所产生的不同影响，并在热分析软件Icepak中对模型进行了验证。由于对Non?Fourier分析的数学处理相对比较困难，本文研究的是相对简单的基板模型，没有考虑芯片与基板的连接、芯片的体积、以及基板上的布线结构等问题，与多芯片组件结构的实际情况有一定的区别，还有待进一步改进。

参考文献

[1] LIU Shu?hua， CAO Li?qiang， LI Jun， et al. Study of mixed?signal crosstalk in 3?D package [C]// Proceedings of 2010 11th International Conference on Electronic Packaging Technology& High Density. Xian， China： ICEPT?HDP， 2010： 559?562.

[2] RICKERT P， KRENIK W. Cell phone integration SiP， SoC， and PoP [J]. IEEE Design & Test of Computers， 2006， 23（3）： 188?195.

[3] KUO W S， TZENG Y L， CHEN E， et al. POP package （Cavity BGA） warpage improvement and stress characteristic analyses [C]// Proceedings of 2007 International Microsystems， Packaging， Assembly and Circuit Technology. Orlando， FL： IMPACT， 2007： 342?345.

[4] 于慈远，于湘珍，杨为民.电子设备热分析/热设计/热测试技术初步研究[J].微电子学，2000，30（5）：334?337.

[5] WEISS J， LANGHORST F. Thermal analysis of electronic packaging [J]. Australian Electronics Engineering， 1988， 21（8）： 62?64.

[6] SWAGER A W. Thermal images provide reliability clues [J]. END， 1990， （1）： 47?58.

[7] 王进春.热传导中的非傅里叶效应研究[J].中国电力教育，2009（3）：144?145.

[8] BERGLES A E. Heat transfer in electronic and microelectronic equipment [M]. Troy， NY： Hemisphere Publishing Corporation， 1990.

[9] 蒋方明，刘登瀛.非傅里叶导热的最新研究进展[J].力学进展，2002，32（1）：128?140.

[10] 张浙，刘登瀛.非傅里叶热传导研究进展[J].力学进展，2000，30（3），446?456.

路径2：[y=35，z=0，x∈（0，70）]

路径3：[y=56，z=0，x∈（0，70）]

路径4：[x=42，y=42，z∈（0，5）]

路径5：[x=21，y=21，z∈（0，5）]

路径6：[x=14，y=56，z∈（0，5）]

路径7：[x=y，z=0，x∈（0，70），y∈（0，70）]

路径8：[x=70-y，z=0，x∈（0，70），y∈（0，70）]

图5分别为在Fourier和Non?Fourier分析两种情况下，沿[x]方向路径1、2、3的温度曲线图。由图5可知，除温度值不同，及热源6位置与其附近区域的温度值在曲线上略有不同外，两种分析情况中的路径曲线走势基本是一致的。

图5 路径1、2、3的温度曲线图（1 300 s）

图6分别为在Fourier分析和Non?Fourier分析两种情况下，沿[z]方向路径4、5、6的温度曲线图。由图6可知，这三条路径上的温度从顶面到底面逐渐降低，且温度变化不大，变化规律基本一致。这是因为沿[z]方向的尺寸很小，并且仅有底面对外传热。

图7的路径7、8是顶面上的两条对角线。路径8的温度曲线左右对称，越靠近路径8的中点，温度越高。整个温度场的分布关于路径7对称。这是因为基板表面布置的热源在也是以路径7为对称轴对称分布的，而且发热量在对称轴附近是最大的。因而热源的分布状况对温度场的影响很大，所以基板上芯片与其他热源的位置分布就显得尤为重要，要尽量分散热源，避免过于集中而导致局部温度过高。

图6 路径4、5、6的温度曲线图（1 300 s）

图7 路径7、8的温度曲线图（1 300 s）

路径7贯穿热源1和热源6所在的位置，Fourier分析中路径8的温度曲线中部比Non?Fourier分析的曲线变化更为剧烈，Fourier分析中曲线的向上凸起程度不及Non?Fourier分析。同时，与热源位置处的温度相比，热源1和6之间相同区域的温度都比较低，而Non?Fourier分析中没有出现这一情况。

由此可见，Non?Fourier分析中的热耦合现象明显强于Fourier分析。

对温度场进行分析可得，使用Fourier分析和Non?Fourier分析得到的结果存在着很多差异。Non?Fourier导热模型的温度上升速度快于Fourier分析，温度曲线比Fourier导热模型的更为陡峭，进入稳定状态的时间远迟于Fourier分析且进入稳定状态后温度值较高，Non?Fourier模型的温度场最大温差要大于Fourier模型，热耦合现象也较强。

4 结 语

本文分别建立了Fourier分析和Non?Fourier分析导热微分方程的有限差分格式，通过求解相应代数方程组，得到两种分析情况下的温度场分布；分析比较了使用Fourier和Non?Fourier分析对基板温度场分布所产生的不同影响，并在热分析软件Icepak中对模型进行了验证。由于对Non?Fourier分析的数学处理相对比较困难，本文研究的是相对简单的基板模型，没有考虑芯片与基板的连接、芯片的体积、以及基板上的布线结构等问题，与多芯片组件结构的实际情况有一定的区别，还有待进一步改进。

参考文献

[1] LIU Shu?hua， CAO Li?qiang， LI Jun， et al. Study of mixed?signal crosstalk in 3?D package [C]// Proceedings of 2010 11th International Conference on Electronic Packaging Technology& High Density. Xian， China： ICEPT?HDP， 2010： 559?562.

[2] RICKERT P， KRENIK W. Cell phone integration SiP， SoC， and PoP [J]. IEEE Design & Test of Computers， 2006， 23（3）： 188?195.

[3] KUO W S， TZENG Y L， CHEN E， et al. POP package （Cavity BGA） warpage improvement and stress characteristic analyses [C]// Proceedings of 2007 International Microsystems， Packaging， Assembly and Circuit Technology. Orlando， FL： IMPACT， 2007： 342?345.

[4] 于慈远，于湘珍，杨为民.电子设备热分析/热设计/热测试技术初步研究[J].微电子学，2000，30（5）：334?337.

[5] WEISS J， LANGHORST F. Thermal analysis of electronic packaging [J]. Australian Electronics Engineering， 1988， 21（8）： 62?64.

[6] SWAGER A W. Thermal images provide reliability clues [J]. END， 1990， （1）： 47?58.

[7] 王进春.热传导中的非傅里叶效应研究[J].中国电力教育，2009（3）：144?145.

[8] BERGLES A E. Heat transfer in electronic and microelectronic equipment [M]. Troy， NY： Hemisphere Publishing Corporation， 1990.

[9] 蒋方明，刘登瀛.非傅里叶导热的最新研究进展[J].力学进展，2002，32（1）：128?140.

[10] 张浙，刘登瀛.非傅里叶热传导研究进展[J].力学进展，2000，30（3），446?456.

路径2：[y=35，z=0，x∈（0，70）]

路径3：[y=56，z=0，x∈（0，70）]

路径4：[x=42，y=42，z∈（0，5）]

路径5：[x=21，y=21，z∈（0，5）]

路径6：[x=14，y=56，z∈（0，5）]

路径7：[x=y，z=0，x∈（0，70），y∈（0，70）]

路径8：[x=70-y，z=0，x∈（0，70），y∈（0，70）]

图5分别为在Fourier和Non?Fourier分析两种情况下，沿[x]方向路径1、2、3的温度曲线图。由图5可知，除温度值不同，及热源6位置与其附近区域的温度值在曲线上略有不同外，两种分析情况中的路径曲线走势基本是一致的。

图5 路径1、2、3的温度曲线图（1 300 s）

图6分别为在Fourier分析和Non?Fourier分析两种情况下，沿[z]方向路径4、5、6的温度曲线图。由图6可知，这三条路径上的温度从顶面到底面逐渐降低，且温度变化不大，变化规律基本一致。这是因为沿[z]方向的尺寸很小，并且仅有底面对外传热。

图7的路径7、8是顶面上的两条对角线。路径8的温度曲线左右对称，越靠近路径8的中点，温度越高。整个温度场的分布关于路径7对称。这是因为基板表面布置的热源在也是以路径7为对称轴对称分布的，而且发热量在对称轴附近是最大的。因而热源的分布状况对温度场的影响很大，所以基板上芯片与其他热源的位置分布就显得尤为重要，要尽量分散热源，避免过于集中而导致局部温度过高。

图6 路径4、5、6的温度曲线图（1 300 s）

图7 路径7、8的温度曲线图（1 300 s）

路径7贯穿热源1和热源6所在的位置，Fourier分析中路径8的温度曲线中部比Non?Fourier分析的曲线变化更为剧烈，Fourier分析中曲线的向上凸起程度不及Non?Fourier分析。同时，与热源位置处的温度相比，热源1和6之间相同区域的温度都比较低，而Non?Fourier分析中没有出现这一情况。

由此可见，Non?Fourier分析中的热耦合现象明显强于Fourier分析。

对温度场进行分析可得，使用Fourier分析和Non?Fourier分析得到的结果存在着很多差异。Non?Fourier导热模型的温度上升速度快于Fourier分析，温度曲线比Fourier导热模型的更为陡峭，进入稳定状态的时间远迟于Fourier分析且进入稳定状态后温度值较高，Non?Fourier模型的温度场最大温差要大于Fourier模型，热耦合现象也较强。

4 结 语

本文分别建立了Fourier分析和Non?Fourier分析导热微分方程的有限差分格式，通过求解相应代数方程组，得到两种分析情况下的温度场分布；分析比较了使用Fourier和Non?Fourier分析对基板温度场分布所产生的不同影响，并在热分析软件Icepak中对模型进行了验证。由于对Non?Fourier分析的数学处理相对比较困难，本文研究的是相对简单的基板模型，没有考虑芯片与基板的连接、芯片的体积、以及基板上的布线结构等问题，与多芯片组件结构的实际情况有一定的区别，还有待进一步改进。

参考文献

[1] LIU Shu?hua， CAO Li?qiang， LI Jun， et al. Study of mixed?signal crosstalk in 3?D package [C]// Proceedings of 2010 11th International Conference on Electronic Packaging Technology& High Density. Xian， China： ICEPT?HDP， 2010： 559?562.

[2] RICKERT P， KRENIK W. Cell phone integration SiP， SoC， and PoP [J]. IEEE Design & Test of Computers， 2006， 23（3）： 188?195.

[3] KUO W S， TZENG Y L， CHEN E， et al. POP package （Cavity BGA） warpage improvement and stress characteristic analyses [C]// Proceedings of 2007 International Microsystems， Packaging， Assembly and Circuit Technology. Orlando， FL： IMPACT， 2007： 342?345.

[4] 于慈远，于湘珍，杨为民.电子设备热分析/热设计/热测试技术初步研究[J].微电子学，2000，30（5）：334?337.

[5] WEISS J， LANGHORST F. Thermal analysis of electronic packaging [J]. Australian Electronics Engineering， 1988， 21（8）： 62?64.

[6] SWAGER A W. Thermal images provide reliability clues [J]. END， 1990， （1）： 47?58.

[7] 王进春.热传导中的非傅里叶效应研究[J].中国电力教育，2009（3）：144?145.

[8] BERGLES A E. Heat transfer in electronic and microelectronic equipment [M]. Troy， NY： Hemisphere Publishing Corporation， 1990.

[9] 蒋方明，刘登瀛.非傅里叶导热的最新研究进展[J].力学进展，2002，32（1）：128?140.

[10] 张浙，刘登瀛.非傅里叶热传导研究进展[J].力学进展，2000，30（3），446?456.