CN114117864A - 自适应时间步长有限元法在电子器件热仿真中的应用方法 - Google Patents

Abstract

自适应时间步长有限元法在电子器件热仿真中的应用方法，属于热场仿真领域。包括步骤：1)构建仿真的几何模型，选定模型中的材料参数；2)对仿真的几何模型进行子区域分割、非相容网格剖分；3)读取子区域网格，对三维模型进行预处理，设定边界条件及材料参数，选取基函数阶数，生成子区域系统矩阵，得到空间离散后的含时矩阵方程组；4)对时间进行离散，自适应时间迭代，求解温度场量，绘制温度分布，并计算求解误差；5)若温度场量的结果不收敛，或误差较大，则通过调整子区域交界面的稳定项系数来修正结果；同时检查设定的迭代误差，重复步骤3)～4)，至获得准确的温度场分布。可减少迭代次数及模型运行时间，提高仿真效率。

Description

自适应时间步长有限元法在电子器件热仿真中的应用方法

技术领域

本发明属于热场仿真领域，涉及到区域分解方法，有限元方法和自适应时间步长相结合的一种自适应时间步长有限元法在电子器件热仿真中的应用方法。

背景技术

现代集成电路朝着超大规模方向发展，芯片的特征尺寸不断逼近理论极限，堆叠层数也不断上升。由于芯片体积有限导致芯片的单位面积产热率不断增加，从而器件工作温度升高。芯片温度过高导致系统可靠性降低，同时减少芯片的使用寿命。芯片的散热问题成为当前的难点之一，并随着硬件设备的性能提升而亟需解决。因此能够快速，高效的对芯片器件建模仿真来指导其散热设计意义重大。

对芯片器件的电磁场与热场仿真的商业软件大多采用有限元方法，但是三维集成电路模型往往是多尺度的，如硅通孔的绝缘层、场效应管沟道之间的间隔相对硅基衬底的尺寸而言相差几个数量级。多尺度的芯片模型导致网格的自由度过大，从而影响系统矩阵求解效率。因此为了减少计算量，仅仅对芯片的局部进行仿真或忽略其中的细小结构进而简化模型。尽管如此，在局部结构仿真和简化时，仍然需要不断提高算法的计算效率，减少计算时间。区域分解法通过采用灵活的子区域非相容网格，而更适用于多尺度模型的多物理场耦合计算。本发明主要围绕如何快速对芯片热分析进行建模和计算，主要发明技术在于：(1)运用基于内罚方法的传输条件，保证温度与热通量在不连续交界面处的连续性；最终生成子区域的系统矩阵并求解，而无需组装整个系统的稀疏矩阵。(2)依据校正项和预测项的温度差以自适应控制时间步长。在初始时刻温度梯度较大时，采用小迭代步，在温度趋近稳定时，采用大迭代步。进而自适应控制步长，减少计算迭代次数，节省计算时间。最后，本发明通过研究一个典型系统级微电子器件热分析，验证恒温、对流等常用边界条件下，新技术的数值精度、计算效率以及稳定性。

发明内容

本发明的目的是为了提高微电子器件在固体传热仿真求解时的运行速度，提供采用基于区域分解的有限元方法减少网格自由度，结合自适应时间步长，减少时间迭代次数，从自由度和时间迭代步两个方面提高仿真效率的一种自适应时间步长有限元法在电子器件热仿真中的应用方法。

本发明包括以下步骤：

1)构建仿真的几何模型，选定模型中的材料参数；

2)对仿真的几何模型进行子区域分割、非相容网格剖分；

3)读取子区域网格，对三维模型进行预处理，设定边界条件及材料参数，选取基函数阶数，生成子区域系统矩阵，得到空间离散后的含时矩阵方程组；

4)对时间进行离散，自适应时间迭代，求解温度场量，绘制温度分布，并计算求解误差；

5)若温度场量的结果不收敛，或误差较大，则通过调整子区域交界面的稳定项系数来修正结果；同时检查设定的迭代误差，重复步骤3)～4)的过程，直到获得准确的温度场分布。

在步骤1)中，所述构建仿真的几何模型是通过确定计算区域、边界条件、热导系数、比热容、材料密度等构建三维模型。

在步骤2)中，所述对仿真的几何模型进行子区域分割、非相容网格剖分，子区域网格尺寸依照子区域大小和材料参数选取，生成多尺度非相容网格。

在步骤3)中，所述空间离散后的含时矩阵方程组包括：

(1)控制方程为三维傅里叶热传导方程：

常用的边界条件：

交界面处的连续性条件：

这里的i与j表示不同的子区域；本发明中采用基于内罚的方式以满足交界面处的连续条件；不同子区域间满足热传导方程表达式为：

其残差表示为：

其中，δ是温度项的稳定系数，不同区域的残差之和为0；权函数取为区域的基函数，采用伽辽金方法，则残差加权变分后得到：

其中，W为权函数。c₁,c₂是常数项系数。式(7)中(·，·)表示子区域上积分计算。<·，·>表示是交界面上的积分；考虑格林恒等式变换，并设置常数项c₁＝c₂＝1抵消微分项

进而简化计算；同理，N个子区域的弱形式可以理解为在N个子区域内求解T∈H¹(Ω₁,Ω₂,..Ω_N)满足：

将第i个区域的表达式总写成矩阵形式为：

其中：

公式(9)即为子区域i空间离散后的系统矩阵。

(2)系统矩阵的计算需要将物理单元变换到参考单元进行积分计算：

物理单元坐标(x,y,z)的积分运算转换成参考单元坐标(ξ,η,ζ)的积分，基函数运算时，采用如下变换：

其中，J为雅可比矩阵。

在步骤3)第(1)部分中，交界面处的稳定项系数δ，其取值正比于δ～O(k/h)，其中k为热导率，h为离散网格的特征尺寸,δ最终的取值范围为10⁴～10⁶。

在步骤3)第(2)部分中，子区域的系统矩阵采用高斯数值积分法：

式中

是四面体单元的面积坐标，W_i是积分点对应的权重，Δ^e为四面体单元体积。

在步骤4)中，所述对时间进行离散，采用自适应时间步长对子区域的时间导数进行离散，离散的策略为：

流程中时间更新公式：

与现有技术相比，本发明的优点如下：

一、本发明采用解耦合的子区域系统矩阵求解，相较于DDM中整体系统求解能够更快速的达到收敛；

二、本发明通过计算预测项与校正项误差，采用初始时刻密集采样，稳态时稀疏采样的自适应时间步长方法进行含时矩阵求解，极大减少迭代次数及模型运行时间。

附图说明

图1是本发明三维堆叠芯片互联模型结构示意图。

表1是本发明实施例的各层材料参数；

图2是本发明实施例的CPU与芯片存储层的发热区域示意图；

图3是本发明中交叠网格的划分：两个不同区域交叠，形成的多边形区域划分为小的三角形；

图4是本发明中求解交界面处的积分流程示意图；

图5是本发明示例在点(-0.0056,-0.0060,-1.0e-4)m计算的温度随时间变化图；

图6是根据本发明实施例计算在x＝1e-6m处的温度截面图；

图7是根据本发明实施例计算在时间为12s的三维模型温度分布图；

具体实施方式

以下实施例将结合附图对本发明的技术方案作进一步说明。

参见图1中三维堆叠芯片互联模型结构模型，本发明实施例包括以下步骤：

S1：构建具有六层结构的三维堆叠芯片互联几何仿真模型；

S2：对仿真的几何模型沿着z切面进行子区域分割，并剖分为密集度不同的网格。记录不同子区域中结点坐标、单元、边界以及面的全局编号信息等；

S3：读取不同子区域网格，对三维模型算法进行预处理，预处理的程序包含对每个区域内的网格基函数进行选定，同时设定热绝缘边界及辐射边界条件，建立子区域系统矩阵方程；

S4：循环计算每个子区域单元的质量矩阵，刚度矩阵以及子区域交界面矩阵，将单元矩阵按照全局结点编号进行组装，最后形成子区域含时系统矩阵；

S5：采用自适应步长离散含时系统矩阵，并运用稀疏矩阵求解器Pardiso得到温度分布；

S6：利用温度系数插值得到整个模型的温度场分布，上述步骤S5中得到的温度场量如果不收敛，或者偏离实际较大，则通过调整控制非相容网格交界面的稳定项系数来修正结果。重复步骤S3～S5的过程，直到获得比较准确的温度场分布。

步骤S3中，关于子区域系统矩阵方程的建立，具体步骤为：

①控制方程为三维热传导方程：

常用的边界条件为：

位于子区域交界面处的连续性条件为：

其中，i与j表示不同的网格区域。本发明中采用基于内罚的方式以满足交界面处的连续条件。不同子区域间满足热传导方程表达式为：

其残差表示为：

其中，δ是温度项的稳定系数，不同区域的残差之和为0。采用伽辽金方法，权函数取值为区域的基函数，则残差加权变分后得到：

其中，W为权函数，c₁,c₂是常数项系数。式(1-42)中(·，·)表示子区域上积分计算。<·，·>表示是交界面上的积分。考虑格林恒等式变换，并设置常数项c₁＝c₂＝1抵消微分项

进而简化计算。同理，N个子区域的弱形式可以理解为在N个子区域内求解T∈H¹(Ω₁,Ω₂,..Ω_N)满足：

将第i个区域的表达式总写成矩阵形式为：

其中：

公式(9)即为子区域i空间离散后的系统矩阵。

②系统矩阵的计算需要将物理单元变换到参考单元进行积分计算：

物理单元坐标(x,y,z)的积分运算转换成参考单元坐标(ξ,η,ζ)的积分，基函数运算时，采用如下变换：

其中J为雅可比矩阵。

在步骤S4中，计算子区域交界面矩阵时，需要先找出所有位于交界面处面的单元信息。由于两个不同的子区域网格非相容，形成的多边形无法直接矩阵积分计算。因此本发明中将交界面上的多边形划分成多个小三角形，根据图3中所示，子区域i与子区域j在交界面处含有多边形交集，将多边形区域分为若干个小三角形。如图3多边形交集划分为小三角形1，小三角形2，小三角形3。记录这些公共小三角形面在子区域i、子区域j的单元全局编号及其对应面的编号信息。在进行第i个子区域的积分计算时，依照交界面矩阵计算流程图4，将该子区域交界面上的小三角形进行遍历，对所有的小三角形进行面积分，获得交界面单元矩阵，最后依据其全局编号信息组装成交界面矩阵，完成整个系统矩阵的计算。

在步骤S5中，对时间进行离散时，采用预测项与校正项的差异以自适应调整单步时间步长进而获得动态采样。本发明采用后向欧拉公式作为校正项，前向欧拉公式作为预测项,在实际的数值仿真过程中，误差值的大小影响步长更新的快慢。误差较大时，得到大时间步长进行更新；小于预测误差时，采用小时间步长进行时间迭代。进而自适应控制时间步长，减少整体计算时间，提高计算效率。

以上为本发明的具体实施流程。

最后本发明给出一个典型系统级微电子器件热分析结果来验证算法的可靠性。三维堆叠芯片互连模型拥有多层互连结构，采用本发明提出的区域分解与自适应步长耦合方法进行热分析求解。图1展示了这个具有六层结构的堆叠芯片互连模型，其中芯片最底部的C4表示系统级芯片中的基板材料。设置该层长宽为为12mm，层高为0.2mm。两个C4聚合层中间是CPU芯片处理层，长宽与底部基底设置相同，长宽同为12mm，高度为0.5mm。CPU层中间的圆柱形表示芯片互连的硅通孔。TSV能够在层之间进行互连，取代引线键合为电路新一代芯片馈电技术。TSV近似为直径0.2mm，高度与层高相同的铜柱。与CPU相同的存储芯片层长宽与CPU层相同，作为器件内存管理芯片。存储芯片层同样也需要TSV来传递电路信息。在存储芯片层之上是热扩散层，用来将底部芯片层的热源进行传导扩散。顶部是热栅指结构，用来对流散热和辐射散热。

表1

材料	W/(m<sup>2</sup>·K)	C(J/[kg·K])	ρ(kg/m<sup>3</sup>)
				热扩散层	400	385	8933
散热片	220	896	2707
				芯片	135	704	2330
通孔	400	385	8933
				基板	50	180	9290

表1显示了具有多层结构的芯片互连系统的每层材料参数取值。TSV的层间布局为：CPU层具有9*5个TSV组成的TSV阵列，芯片存储层具有6*4的TSV阵列结构，而热扩散层有3*4的TSV用来导热。模型将CPU芯片发热层与芯片存储层划分为4个大小相同的部分。图2中给出了热源部分的CPU芯片发热层和芯片存储层的划分。两个芯片层结构等效为体热源，四等分的热源取值不同。分别遵循：

同时Q_c＝|cos(2πy/3l₁)|，

且Q_m＝|cos(3πy/3l₁)|。热扩散层也等效成热耗率为P＝10W的体热源。顶部热栅指的表面作为对流散热和辐射散热的边界。等效的边界传热系数为h＝500W/(m²·K)，表面辐射率为0.93W/m²。模型初始温度为300K，环境温度为300K。该模型被分成6个不同子区域；子区域间是非相容网格。

图5给出商业软件COMSOL和本发明方法计算得到的观察点(-0.0056，-0.0060，-1.0e-4)m随时间温度变化分布情况，从图中可以看到：本发明的区域分解自适应步长计算技术与COMSOL计算结果一致，吻合很好，整体误差为1.0％。图6给出芯片互连模型在x＝1e-6m处的切面温度分布。YZ切面的温度分布中底部的分布比较相近，顶部温度较低，且结构的层与层之间存在温度差异，同层之间温度分布相同。越靠近顶部热栅指温度变化越大。此外，因给定的热源是根据空间位置变化的，芯片发热层的中间部分温度聚集，边缘部分温度略低。图7给出三维芯片互连模型在时间为12秒的温度分布图，温度单位为K。

上述实施例仅为本发明的较佳实施例，不能被认为用于限定本发明的实施范围。凡依本发明申请范围所作的均等变化与改进等，均应仍归属于本发明的专利涵盖范围之内。

Claims (7)

1.自适应时间步长有限元法在电子器件热仿真中的应用方法，其特征在于包括以下步骤：

1)构建仿真的几何模型，选定模型中的材料参数；

2)对仿真的几何模型进行子区域分割、非相容网格剖分；

3)读取子区域网格，对三维模型进行预处理，设定边界条件及材料参数，选取基函数阶数，生成子区域系统矩阵，得到空间离散后的含时矩阵方程组；

4)对时间进行离散，自适应时间迭代，求解温度场量，绘制温度分布，并计算求解误差；

5)若温度场量的结果不收敛，或误差较大，则通过调整子区域交界面的稳定项系数来修正结果；同时检查设定的迭代误差，重复步骤3)～4)的过程，直到获得准确的温度场分布。

2.如权利要求1所述自适应时间步长有限元法在电子器件热仿真中的应用方法，其特征在于在步骤1)中，所述构建仿真的几何模型是通过确定计算区域、边界条件、热导系数、比热容、材料密度构建三维模型。

3.如权利要求1所述自适应时间步长有限元法在电子器件热仿真中的应用方法，其特征在于在步骤2)中，所述对仿真的几何模型进行子区域分割、非相容网格剖分，子区域网格尺寸依照子区域大小和材料参数选取，生成多尺度非相容网格。

4.如权利要求1所述自适应时间步长有限元法在电子器件热仿真中的应用方法，其特征在于在步骤3)中，所述得到空间离散后的含时矩阵方程组包括以下：

(1)控制方程为三维傅里叶热传导方程：

常用的边界条件：

T|_Γ1＝T₀(r,t) 恒温边界

绝缘边界

对流边界

交界面处的连续性条件：

T⁽ⁱ⁾＝T^(j)

这里的i与j表示不同的子区域；采用基于内罚的方式以满足交界面处的连续条件；不同子区域间满足热传导方程表达式为：

其残差表示为：

其中，δ是温度项的稳定系数，不同区域的残差之和为0；权函数取为区域的基函数，采用伽辽金方法，则残差加权变分后得到：

其中，W为权函数；c₁,c₂是常数项系数；式(7)中(·，·)表示子区域上积分计算；<·，·>表示是交界面上的积分；考虑格林恒等式变换，并设置常数项c₁＝c₂＝1抵消微分项

进而简化计算；同理，N个子区域的弱形式可以理解为在N个子区域内求解T∈H¹(Ω₁,Ω₂,..Ω_N)满足：

将第i个区域的表达式总写成矩阵形式为：

其中：

公式(9)即为子区域i空间离散后的系统矩阵；

(2)系统矩阵的计算需要将物理单元变换到参考单元进行积分计算：

物理单元坐标(x,y,z)的积分运算转换成参考单元坐标(ξ,η,ζ)的积分，基函数运算时，采用如下变换：

其中，J为雅可比矩阵。

5.如权利要求4所述自适应时间步长有限元法在电子器件热仿真中的应用方法，其特征在于在步骤3)第(1)部分中，交界面处的稳定项系数δ，其取值正比于δ～O(k/h)，其中k为热导率，h为离散网格的特征尺寸,δ最终的取值范围为10⁴～10⁶。

6.如权利要求4所述自适应时间步长有限元法在电子器件热仿真中的应用方法，其特征在于在步骤3)第(2)部分中，系统矩阵采用高斯数值积分法：

式中

是四面体单元的面积坐标，W_i是积分点对应的权重，Δ^e为四面体单元体积。

7.如权利要求1所述自适应时间步长有限元法在电子器件热仿真中的应用方法，其特征在于在步骤4)中，所述对时间进行离散，采用自适应时间步长对子区域的时间导数进行离散，离散的策略为：

流程中时间更新公式：

Images