用于分析电子封装结构的有限元-边界元耦合方法研究

马冲 公颜鹏 侯传涛 秦飞

（1电子封装技术与可靠性研究所，北京工业大学，北京 100124；2 可靠性与环境工程技术重点实验室，北京强度环境研究所，北京 100076）

0 引言

芯片封装是指利用膜技术及微细加工技术，将芯片及其他要素在框架或基板上布置、粘贴固定及连接，引出接线端子并通过可塑性绝缘介质灌封固定，构成整体立体结构的工艺[1]。封装结构涉及多种不同性能的材料，而且在界面处存在几何不连续性。在服役过程中，温度场、湿度场等恶劣环境会使封装结构产生多种失效问题，例如，空洞、开裂、翘曲和分层等。随着封装密度增加和功能的多样化，具有多物理场、多尺度特征的电子封装结构的可靠性问题正成为设计的关键。

由于易实现、精度高、无环境限制等特点，数值模拟技术正成为电子封装结构设计领域的主流方案。有限元法（FEM）是一种常用的获得各种工程问题近似解的数值分析技术。目前，商业有限元软件ABAQUS、ANSYS等被广泛用于电子封装结构的可靠性分析[2-7]。

实际上，由于电子封装结构的多尺度特征，有限元模型往往需要被离散成大量的单元，以保证数值结果的合理性。特别是对于一些具有多尺度特征的复杂模型，会耗费大量计算时间和成本。边界元法（BEM）是用以求解工程和科学问题的常用数值分析方法之一，其特点是降低求解问题的自由度数和只需要边界离散化。有限元法和边界元法各有特点，因此发挥各自特长的有限元-边界元耦合方法一直受到研究者的重视。Liu等人[8]提出了一种用于解决动态弹塑性问题的有限元-边界元耦合方法。Fernandes 等人[9]提出了一种分析非均匀材料板弯曲问题的有限元-边界元耦合方法。Bonari 等人[10]提出了一种用于求解接触问题的有限元-边界元耦合方法，其中有限元用于处理宏观模型、边界元用于处理微观粗糙模型。Dong 等人[11]利用边界元法对力载荷作用下的集成电路封装进行应力计算。基于模型几何和加载的对称性质，得到相应的基本解。Neto 等人[12]提出了一种将等几何边界元法（IGABEM）与有限元软件相结合的方法。此研究工作为弹性力学中的三维IGABEM公式提出了网格自适应方法，提供了精确的几何表示和力学场描述，有助于BEM和CAD方案的耦合。Qin等人[13]提出一种新颖的有限元-边界元耦合方法来研究多尺度结构二维稳态传热问题，其中边界元用于处理不重要或线性区域。Gong等人[14]提出了一种能够研究电子封装中多尺度结构传热问题的等几何边界元方法。数值结果与解析解和有限元解的比较验证了此方法的有效性。虽然边界元法可以减少单元的数量，但对于非线性、非均匀问题的分析，该方法还存在很多问题。

结合边界元法和有限元法各自的优势，本文提出了一种用于分析电子封装结构的有限元-边界元耦合方法，将自编写的BEM代码集成到商业软件ABAQUS中。在耦合方法中，将整个模型域划分为有限元（FE）域和边界元（BE）域。有限元法（借助ABAQUS）用于分析具有非线性或非均匀特性的子域，而具有线性特征或不重要的域则由BEM代码求解。边界元域被视为有限元法的一个超单元，它的有效刚度和有效节点力可通过BEM代码求解得到。通过子程序UEL将超单元的等效量装配到整体模型的有限元公式。通过这种方法，实现用ABAQUS进行电子封装结构的耦合方法分析。

1 有限元-边界元耦合理论

有限元法与边界元法的耦合方法需要把整个求解域分成两部分，一部分以有限元法进行离散，另一部分以边界元法进行离散。求解方法基本上分为两类：一类是先求解边界元域，把边界元区域看作是有限元法的一个超单元，通过边界条件和有限元方程耦合起来求解；另一类是先求解有限元域，然后将有限元解作为边界元域的边界条件，再用边界元法求解[15]。本文采用第一类求解方法。

对于先求解边界元域的方法，边界元域的边界面力应转化为有限元中的等效节点力。如图1所示的模型，由一个BE域和一个FE域组成。对于BE域，其界面的位移{uc}和面力{tc}之间的关系可以表示为[15]

图1 边界元域与有限元域的界面力Fig.1 Interface forces between FE domain and BE domain

式中，{tc0}是所有界面位移为零时的面力向量，[KBE]是边界元域的伪刚度矩阵。

对于FE域，界面位移{uc}和界面节点力{Fc}之间的关系为

式中，{Fc0}是界面位移为零时的节点力向量，[KFE]是有限元域界面节点的刚度矩阵。

在界面的某一节点上施加x方向的任意虚位移，则有限元域界面节点力做的功为

式中，n为单元e的节点数。根据守恒条件，节点力做的功应与面力做的功相等

式中，Γ为积分域边界（如图1所示）。将面力tx和虚位移xuδ表示为插值形式

其中，Ni和Nj为界面处的形函数。

由公式(3)-(6)，可以得到

同理，施加y向虚位移，可得到

因此，对公式(1)左乘矩阵[M]，得到等效节点力的表达式

其中，由[Me]组合而成的矩阵[M]表示整个界面的变换矩阵，矩阵[Me]可以通过下式获得

由公式(9)可知，可以通过转换矩阵[M]建立有限元法与边界元法的关系。将公式（9）带入公式（2），即可形成完整的耦合方程。

2 耦合算法的实现

本文所开发的有限元法-边界元法耦合算法是通过子程序UEL把边界元程序集成到有限元软件ABAQUS中。图2展示了将自编写的边界元程序与ABAQUS耦合起来的过程：

图2 耦合算法实现流程图Fig.2 The flowchart of coupling scheme

1）根据模型的几何信息、材料特性等参数将模型划分为有限元子域和边界元子域。其中有限元子域需通过ABAQUS /CAE处理，而边界元子域通过边界元程序进行计算。

2）为了使边界元程序能与ABAQUS软件耦合在一起，需在ABAQUS/CAE中以标准操作建立两个模型：有限元域模型和等效边界元模型。等效边界元部分是将边界元域等效为与界面相重合的几何模型， 即接触界面。需要注意的是应该保证等效边界元部分所划分的节点与有限元域界面部分的节点保持一致，然后两个模型通过绑定约束“ Tie ”组装起来。

3）ABAQUS/CAE完成两个部件的前处理后，输出存储模型信息的初始input文件。

4）对初始input文件进行修改，将等效边界元部分修改为用户自定义单元，重新定义其几何信息、物理属性和边界条件等参数。

5）将修改后的input文件提交计算，同时调用包含自编写边界元程序的子程序UEL来实现边界元程序与ABAQUS的耦合。

6）通过UEL得到用户自定义单元（即边界元域界面处）的等效刚度矩阵和等效节点载荷，并传递给ABAQUS，此时有限元域便可通过ABAQUS求解器进行求解得到有限元域的所有物理量。

3 数值算例

本节通过两个算例来实现该算法对常见结构的分析并探索其在实际问题中的应用，进一步探究该算法的性能。

3.1 具有复杂界面的矩形板结构

在这个算例中，我们选择带有复杂界面形状的矩形板结构，该模型由两个部分组成，其模型示意图如图3所示。该模型的几何参数为a= 1 mm，l1= 0.2 mm，l2= 1 mm。边界条件为悬臂梁左侧固定，右侧有沿x方向分布的均布载荷q= 1000 N/mm。该模型两个域的弹性模量为E= 1000 MPa，泊松比为v= 0.2。耦合方法对模型求解时，左侧区域用有限元法离散，右侧区域用边界元法离散。为了验证该耦合算法的精度，在此选用商业软件ABAQUS的细化网格模型的计算结果作为参考解。耦合算法处理的模型被离散成1个自定义单元和800个四节点单元，而细化的有限元模型被离散成5200个四节点单元。针对该矩形板模型，从提交作业到分析完成的过程，耦合方法用时24 s，有限元法用时50 s。图4 (a) 和 (b) 所示为该耦合方法和有限元法得到的左侧部分的von Mises应力云图。从图中可以看出该耦合算法对于弹性问题的求解有着较好的精度。

图3 带有复杂界面的多尺度结构示意图Fig.3 Multiscale structure with complex interfaces

图4 两种方法得到的左侧结构的von Mises应力分布云图Fig.4 von Mises stress contours computed by the present method and FEM

3.2 电子封装TSV-Cu结构

此模型为三维电子封装中较为先进的硅通孔（Through silicon via, TSV）结构。TSV-Cu 技术是在硅片上通过干法刻蚀技术在硅基体中刻蚀孔洞，然后在其中电镀填充铜来实现晶体管的机械支撑和电气互连等目的。后续需要通过CMP技术去除沉积在硅表面的铜层，并实现晶圆背面TSV结构的暴露及结构正、背面的平坦化[16]。

在此算例中，对简化的TSV结构模型（如图5所示）进行力学分析。其中TSV-Cu 部分的直径为10 μm，深度为100 μm，电镀Cu淀积层的厚度为δ= 6μm；硅部分的几何尺寸分别为a= 630 μm，b= 250 μm。硅部分的底部被固定，在淀积层的上表面施加剪切载荷t= 0.25 N/m。TSV-Cu部分的弹性模量ECu= 155 GPa，泊松比v= 0.3，屈服强度为47.91MPa；Si部分的弹性模量ESi= 131 GPa，泊松比v= 0.25。

图5 TSV模型的示意图Fig.5 Cross-section of TSV

在耦合算法计算过程中，TSV-Cu部分以有限元法划分网格，而硅部分以边界元法划分边界 网格，如图6所示的网格示意图，在TSV - Cu区域使用四节点四边形单元，而在硅部分使用二次边界单元进行离散。耦合算法处理的模型被离散成1个自定义单元和7434个四节点单元，而细化的有限元模型被离散成161934个四节点单元。由于硅部分只需在边界部分进行离散，所以该耦合方法对于这种跨尺度结构可以大幅度减少单元数量。对于硅通孔模型，从提交作业到分析完成的过程，耦合算法用时82 s，有限元法用时102 s。

图6 TSV结构的网格示意图Fig.6 Meshes used in the TSV model

图7展示了该耦合算法和有限元法计算得到的位于路径L上的节点处的von Mises应力结果，有限元法的计算结果作为参考结果。从图中可以清晰的看出该耦合算法对于非线性问题的求解也有着较好的精度。

图7 路径L上各点的von Mises应力分布值Fig.7 The von Mises stress at points along the contour L

4 结论

本文提出了一种将商业软件ABAQUS和自编BEM代码结合的耦合方法。在当前的耦合方法中，根据结构的几何特征、材料属性等，将整个模型分为BE域和FE域。自编 BEM 代码的功能通过 ABAQUS的子程序 UEL 来实现，以此获得用户自定义单元的有效刚度和有效节点力。

数值算例说明了该算法在求解弹性问题和非线性问题时均有着较高的精度，能大幅降低单元数量并减少计算时间。