基于RVE-子模型法的多尺度封装结构分析方法

孙国立 公颜鹏 侯传涛 李尧 秦飞

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

0 引言

随着集成电路发展进入后摩尔定律时代，集成电路产业已经从传统的不断减小结构尺寸的方法转向通过先进封装技术提升系统性能[1,2]。近年来，系统级封装、异质集成、3D封装等先进封装技术得到越来越多的关注[3,4]。封装后的芯片以封装模块的形式服役，目前，载荷环境愈发复杂[5]。当封装模块处在恶劣的工作环境时会导致模块出现可靠性问题[6]，尤其对于弹载的电子产品，其可靠性对飞行安全至关重要[7]。因此对完整的封装结构进行仿真建模分析尤为重要。但是，完整的封装结构中会存在微米尺寸的RDL和TSV-Cu，同时也存在毫米尺寸的基板，甚至会在界面处形成了相互镶嵌的微观组织[8]。因此，封装模型具有典型的结构多尺度特征[9]。虽然计算机硬件和软件水平已有长足发展，但结构多尺度会导致计算规模迅速增加，耗费大量计算资源。同时计算过程中各类误差不断累积，不仅影响计算精度，计算的收敛性也面临挑战。因此，亟待发展一种多尺度分析方法，实现几何多尺度结构的快速精确分析。

目前，针对多尺度封装结构的分析，最简单常用的方法是体积等效法，该方法可以很好地估计两相复合材料的等效热膨胀系数[10-12]和等效弹性模量[10]。根据该均匀化方法，Lee研究发现当铜通孔占总体积14%以下时，等效杨氏模量接近于硅，而当铜的含量为88%时弹性模量接近于铜[13]。Cheng 和 Shen使用有限元计算的热应变来估计等效的热膨胀系数[11]。这种方法只能粗略的计算平面内的等效参数。Che等用同样的方法根据TSV-Cu在晶圆上的分布，得到铜/硅胞体的线弹性力学参数[14]。但以上方法并没有考虑TSV-Cu塑性和胞体的边界条件的影响，建立的等效模型与真实模型仍存在一定误差。一些学者通过理论推导的方式得到解析的等效参数。Chen等人对 TSV转接板的平面内和平面外等效机械性能的理论进行了详细阐述并验证其有效性[15]。在Chen基础上，Ye等人通过在计算铜和烧结银的等效力学参数时引入几何参数，对其计算理论进行了修正，数值结果表明全网格模型与等效模型应力分布一致[16]。但是，由于封装结构的复杂度越来越高，这种方法的效率较低。因此，亟需发展一种封装结构多尺度分析方法。由于三维编织材料的细观结构呈周期性排布，因此利用代表性体积单元（Representative Volume Element，RVE）研究三维编织材料是一种有效且常用的方法[17]。与编织材料类似，几何周期性也同样普遍存在于封装结构中。

本文基于代表性体积单元法和子模型技术，将两种方法耦合，发展了RVE-子模型法。建立多尺度封装结构的等效模型，在保证计算精度的情况下，达到减小计算规模和降低计算资源的目的。基于该方法本文对三点弯仿真建立等效模型，通过对关键位置各应力、应变分量的对比，验证了本文方法的有效性。

1 RVE-子模型法理论

多尺度结构建模的中心思想是尺度分离。对于周期性结构，利用直接平均理论计算单个RVE的本构，并将其赋给所建立的均质等效模型。再通过子模型技术，将重点关注的位置建立精细模型（子模型），得到关键区域的应力应变场。该方法可以对具有结构多尺度的非线性周期性结构进行分析。下面我们以TSV-Cu结构为例介绍本文算法。

1.1 RVE分析法

图1为建立等效模型的示意图，在实际的封装结构中，TSV-Cu直径在30-100 μm之间，在芯片或转接板上呈周期性分布，宏观结构可由单个胞体阵列得到。根据其结构特点，在某一区域内定义一个单胞作为RVE，并将RVE的力学参数作为整个周期性结构的力学参数。

图1 建立等效模型示意图Fig.1 Schematic diagram of the established equivalent model

RVE的边界条件对获得准确本构具有重要影响，研究表明施加周期性边界条件使得单个RVE变形协调，相邻RVE界面之间不会分离，可以保证相邻RVE界面处的应力连续[18]。式（1）为任意RVE周期性边界条件。

在ABAQUS中实现对RVE划分周期性网格，并采用多点约束的方式施加周期性边界条件。式（2）为施加的约束方程，使得相对面上节点在某一方向的位移差为一常数。计算中，需要对RVE每个面上相对节点之间自由度建立约束，由于节点数量过于庞大，本文通过 Python脚本实现线性多点约束。

式中，Cn为根据节点和参考点的位置确定的系数，n为线性多点约束方程的项数（本文取3），是位移变量，i表示方向（=1,2,3i），P表示节点的位置。

对RVE进行单向拉伸，并利用公式（3）和（4）所示的直接平均方法得到如图4所示的RVE的应力应变曲线，即等效模型的本构。

式中N为RVE模型的单元总数，Vm是RVE中第m个单元的体积，σm是RVE中第m个单元的应力，mε是RVE中第m个单元的应变。

1.2 子模型技术

子模型方法又称特定边界位移法，特定边界就是子模型从全模型中切割出的边界。在同一坐标系中对部分区域进行精细建模，并将其作为子模型，其边界条件即为全模型边界的位移值。子模型技术基于圣维南原理。下简要说明子模型的有限元理论。式（5）为基于弹性力学最小位能原理的有限元表达式[19]。

式中，[K]为整体刚度矩阵，[D]为待求位移向量，[F]为施加的载荷。

由于子模型的边界条件是整体模型切割位移值，因此部分位移解已知。假设位移向量[D]已知的部分为[D1]，待求解的为[D2]。则式（5）可以写为

将式（6）展开得到公式（7）

由（7）式可知，已知的位移向量[]1D转变为载荷项，可求得未知的[D2]。

本文将等效模型得到的位移场作为子模型边界条件，利用节点位移和单元应力，得到关注位置的应力应变场。

2 有限元求解

本文以ABAQUS软件为实验平台，对RVE进行单轴拉伸，采用直接均匀理论得到每一个子步下平均应力和平均应变值[20]。图2利用RVE-子模型法建立计算宏观位移场和关注位置的应力应变相应的流程图。求解过程通过商业软件ABAQUS和Python脚本实现。假设RVE为各向同性材料，且由RVE在x向和y向阵列得到的非线性复合板可以等效为均质板。首先，对RVE施加周期性边界条件并进行单轴拉伸得到其应力应变曲线。并将获得的RVE本构赋给建立的等效模型，求解得到等效模型的位移场，然后，与全网格模型的结果进行对比，若误差较小则将等效模型的位移作为子模型的边界条件，可得到关注位置的应力应变场。若误差较大，则重新选取RVE后再次进行求解。

图2 RVE-子模型法计算流程Fig.2 RVE-Submodel method calculation process

3 算例验证

3.1 三点弯模型

如图3所示，考虑一种非线性周期性结构的三点弯模型，模型由15个单胞组成，每个单胞的尺寸均为111×× mm，中间支撑圆柱的半径为0.08 mm，整个模型的最大和最小尺寸相差两个数量级，具备一定的多尺度特征。为体现RVE-子模型在处理非线性问题上的优势，设置复合板的基体为TSV-Cu，其本构如图4所示[21]，增强体为硅，弹性模量为131 GPa，泊松比为0.3[22]。由于基体材料为非线性本构，因此本文假设等效模型的也为非线性结构，通过第2节所述方法得到等效模型的非线性本构。有限元模型的两个支撑圆柱和中间压头圆柱均被约束为刚体，两端的支撑为固定约束，中间压头施加向下的0.3 mm位移。三个圆柱压头与中间的复合板之间采用绑定约束。全网格模型和基于RVE-子模型法建立的等效模型参数对比表1所示。为排除计算机硬件的影响，采用同一台计算机和相同的单元类型（C3D8R）对两种模型进行划分网格，得到的全网格模型的单元数为835165，而等效模型的单元数为371925，等效模型单元约为全网格模型的一半。本算例假设关注区域为图3所示的RVE区域，因此，需要对该位置建立精细模型并将等效模型的计算结果作为子模型的边界条件。

图3 三点弯模型Fig.3 Three-point bend model

图4 模拟单轴拉伸下的应力应变曲线Fig.4 Stress-strain curve under uniaxial tensile

表1 有限元模型参数Table 1 Parameters used in the finite element model

3.2 等效模型的有效性验证

为验证本文RVE-子模型法所建立的等效模型的有效性。首先对比了相同载荷下两种模型的宏观位移场，其次对比了关注位置的全网格模型和子模型的计算结果，最后比较了两种模型在相同计算机配置下的计算时间。图5为模拟三点弯曲实验下实际模型和等效均质模型位移云图。由图可知两种模型得到的最大位移相差仅7.6 μm，并且两种方法计算结果的趋势一致。

图5 全网格模型与等效模型的位移云图对比Fig.5 Comparisons of displacement obtained by the full mesh model and present method

图6 给出了本文RVE-子模型法得到的关注区域位移值与全网格模型得到的值对比情况，ux、uy、分别为关注区域x、y、z三个方向的最大位移值和最大合位移。从图中可以看出两种方法得到位移值基本一致，最大相差0.3 mm，结果吻合较好。

图6 关注区域各位移分量对比(ux、uy、uz、u)Fig.6 Comparisons of displacement (ux、uy、uz、u)

图7 为关注位置的应力分量对比，由图可知本文算法得到的最大值与参考值趋势相同。其中子模型得到的关注位置最大von Mises应力为1464 MPa，参考值为1114 MPa，两值相差较小。因此，可为材料失效分析提供依据。x方向正应力模拟值为-786.3 MPa，与参考值的误差仅为5.02%。许多经典疲劳理论是基于应变理论的（如Manson-Coffin等），因此应变对于结构可靠性有重要影响。关注位置的应变分量对比在图8给出，本文数值解和参考值在同一水平，与图7对比可知，应力大的分量应变也相应大，这与力学理论一致。两种模型的计算时间如表2所示。与全模型相比，等效模型在计算结果精度一致的情况下，计算时间为40.85 min，比全模型减少了78.57%。等效模型和子模型的计算时间之和也仅有全模型的五分之一，这进一步体现了RVE-子模型法处理结构多尺度的优势。

表2 不同模型的计算时间Table 2 Computation time for different models

图7 关注区域各应力分量对比Fig.7 Comparison of stress components in the area of interest

图8 关注区域各应变分量对比Fig.8 Comparison of strain components in the area of interest

结果证明本文提出的方法用于研究多尺度封装结构的可行性。该方法可以降低计算规模，提高计算效率。

4 结论

针对封装结构中存在大量多尺度结构的问题，本文将代表性体积单元法和子模型技术相结合发展了RVE-子模型法。利用ABAQUS可以解决封装结构中大量周期结构多尺度建模的困难，实现宏观和微观的双尺度交互。文章给出了该方法的具体计算过程，并通过三点弯仿真验证了方法的有效性。结果表明， RVE-子模型法节省了大量计算时间，加快了计算进程，为解决先进封装结构多尺度建模问题提供一种新方法。