[YuL4]
图6 CBGA焊点的有限元模型
表2 有限元模型的材料参数
2.2 稳健性分析
仿真和实验数据的变化趋势一致性直接影响寿命模型的准确度。为使仿真结果和实验数据获得稳定的拟合匹配关系，研究提取了仿真结果中每颗焊点的应变能密度、等效塑性应变的最大值增量和平均值增量4种指标，对它们的变化趋势和稳健性进行了分析。其中，指标的平均值是单元体积加权平均值。以平均塑性应变能密度增量为例，焊点的仿真结果分布趋势及加权平均值求解方法如图7所示。
[YBC5]
图7 塑性应变能密度仿真结果
2.2.1 单元类型的选择
单元与结点是有限元方法中的核心要素，应力、应变、能量等一切数值计算都离不开单元类型的定义。在CBGA板级封装的温循过程中，焊球受到上侧陶瓷基板和下侧PCB的双重作用，再考虑到其自身的热胀冷缩效应，焊球处于拉伸-压缩-剪切的耦合状态，其受力特征比较复杂，对有限元模型的求解准确性要求较高。
在单元属性方面，可能影响仿真结果的主要有单元阶次、积分点数量等因素。单元阶次主要有一阶和二阶2种，二阶单元相比一阶单元的节点数量更多，求解精度更高，但计算效率较低。积分点形式主要有全积分和减缩积分2种，全积分形式的求解精度更高，但在配合一阶单元使用时存在“剪切自锁”现象，导致对弯曲的描述失真；减缩积分可避免这种现象，但积分点的减少又导致了“沙漏”问题。为确保有限元仿真的稳健性，研究使用5类六面体单元对焊球的状态进行了仿真对比分析。5类单元属性详见表3。其中，单元类型1和类型2为二阶单元，类型3和类型4为一阶单元，后缀字母R的为减缩积分单元，类型5是在一阶减缩积分单元C3D8R的基础上进行了沙漏控制强化处理。
表3 5类六面体单元类型属性概况
（a） 平均塑性应变能密度增量[YuL6] 随中心距变化曲线
（b） 平均等效塑性应变增量随中心距变化曲线
（c） 最大塑性应变能密度增量随中心距变化曲线
（d） 最大等效塑性应变增量随中心距变化曲线
图8 不同单元类型的仿真结果稳健性分析
图8为不同单元类型的仿真结果稳健性分析。通过对比可以发现，平均塑性应变能密度增量受单元类型的影响最小，使用5种单元类型能获得几乎一致的仿真结果，并且随着焊点中心距的增大，其位置越来越靠近基板边缘，平均塑性应变能密度增量呈单调上升的趋势，符合焊点可靠性的空间分布特征。相比之下，其余3种表征因子随中心距的变化趋势一致性相对较差，受单元类型的影响都很大。
2.2.2 网格收敛性分析
网格划分是有限元方法的重要内容，网格密度是否合理直接关系到仿真精度和可复现性。下面针对焊球单元基本尺寸分别为0.015 mm、0.030 mm和0.060 mm时的仿真结果进行了对比分析，见图9。
[YuL7] [YuL8]
（a） 不同网格密度的焊球模型
（b） 平均塑性应变能密度增量随中心距变化的曲线
（c） 平均等效塑性应变增量随中心距变化的曲线
（d） 最大塑性应变能密度增量随中心距变化的曲线
（e） 最大等效塑性应变增量随中心距变化的曲线
图9 网格密度对比
由图9可以看到，塑性应变能密度的平均值增量受网格密度的影响最小，在3种不同的网格密度下的仿真结果几乎完全一致，说明该表征因子的网格收敛性很好。此外，塑性应变的最大值增量受到的影响也相对较小，而其余2种表征因子受网格密度的影响较大。因此，从网格收敛性的角度看，应变能密度平均值增量仿真结果不易受单元尺寸影响，能获得更稳健的表征效果。
2.2.3 分析区域的选择
下面对于可靠性表征因子的提取位置进行了对比研究。分别选取了焊球在靠近陶瓷侧的薄层、靠近PCB侧的薄层、整个焊球3种区域，对塑性应变能密度、等效塑性应变的平均值增量和最大值增量几个表征因子进行分析。中心距越大（越靠近边缘）的焊点越容易发生失效，可靠性越低，焊点的可靠性和中心距呈负相关。由图10可以看出，焊球全域范围的应变能密度平均值增量和等效塑性应变平均值增量2个指标状态和中心距的相关关系最明确，能更准确地表征焊点可靠性特征，而陶瓷侧或PCB侧的薄层区域内的结果指标变化趋势比较离散，与焊点的可靠性规律不相符，因此采用全域内的指标平均值增量更为合理。
[YuL9]
（a） 区域示意图
（b） 平均塑性应变能密度增量随中心距变化的曲线
（c） 平均等效塑性应变增量随中心距变化的曲线
（d） 最大塑性应变能密度增量随中心距变化的曲线
（e） 最大等效塑性应变增量随中心距变化的曲线
图10 分析区域对比
焊球整体的应变能密度平均值增量能适应不同单元类型和网格密度，并且与焊点的可靠性空间分布特征一致，因此本研究选取该指标作为CBGA焊点的可靠性表征因子。
3 寿命模型
3.1 模型建立
结合仿真表征和试验标定，研究团队建立了焊点的疲劳寿命预测模型，从而可以根据仿真结果预测焊点在温度循环条件下的寿命。目前，电子封装领域对BGA焊点寿命的预测方法主要有4种模型[22]，包括Coffin-Manson模型、Darveaux模型、断裂力学模型和损伤线性叠加模型，详见表4。CM模型用塑性应变来表征焊点失效机理；Darveaux模型则提出形变能量是焊点疲劳寿命的关键因素；断裂力学模型理论以裂纹扩展速率相关的参数来表征焊点失效；损伤线性叠加模型则认为稳态蠕变是失效的主要原因。根据前文中仿真模型稳健性部分的研究结论，平均应变能密度增量不易受有限元模型中各种因素的影响，有限元仿真结果更为稳定，因此本研究选取Darveaux模型，并以平均塑性应变能密度增量ΔWave作为表征因子。
表4 4类寿命预测模型
Darveaux模型是基于能量法的寿命模型，主要通过测量焊点的裂纹扩展率，建立与裂纹扩展相关的裂纹萌生和扩展的寿命方程，2个方程都用单次温度循环中的平均塑性应变能密度增量来表达。其中，裂纹萌生的温度循环次数N0满足：
[YuL10] （1）
裂纹扩展过程需要的温度循环次数N1满足：
N1=a(da/dN) （2）
分子a为裂纹扩展总长，分母为裂纹扩展速率，也和能量相关：
[YuL11] （3）
N0与N1之和为焊点总寿命：
Nf = N0 + N1 （4）
由于裂纹扩展的速度很快，温度循环次数一旦达到N0，则萌生裂纹，往往在几次循环后焊点就完全断裂，因此总寿命Nf≈N0。
将同一CBGA板级封装结构的实验寿命数据与仿真所得数据ΔWave进行对比，如图11所示。所有焊点根据其到结构中心的距离从小到大排序，结果显示焊点的MTTF呈现出逐渐下降的趋势，而仿真得到的ΔWave值则呈现逐渐增长的趋势[YBC12] 。虽然2组数据变化方向相反，但变化幅度相似，显示出较强的负相关性。
图11 不同位置的焊点平均寿命和平均塑性应变能密度增量随中心距变化的曲线
通过拟合每种中心距焊点的试验数据与仿真数据，焊点寿命的Darveaux能量模型拟合曲线如图12所示，可以看出拟合优度较高，数据点均分布在拟合曲线周围。通过这样的方式，利用Darveaux模型将试验数据与仿真表征有机结合，得到模型参数K1=226.12，K2=-0.89。通过这些参数建立了仿真表征和实验数据之间的桥梁，在CBGA产品研发阶段，可以利用该模型评估相同制备和封装工艺下的焊点可靠性，提高产品研发效率。
[YuL13] [YuL14] [YuL15]
图12 焊点寿命的Darveaux能量模型拟合曲线
3.2 寿命预测案例
研究团队尝试使用寿命模型预测某型新研发产品的焊点寿命。建立了该封装产品的板级有限元模型，设定焊球直径为0.5 mm，并进行了高低温循环的有限元仿真，温度范围为-55~100 ℃。仿真所得的代表产品板级焊点的应变能密度分布云图见图13，可以看出，越靠近结构外侧的焊球，其塑性应变能越高。
[YuL16]
图13 代表产品板级焊点的应变能密度分布云图
随后，提取了关键电路焊球的可靠性表征因子（即ΔWave），其值为0.25×10-3 J/mm3。将所得的可靠性表征因子带入前文所述的寿命模型，计算得到该CBGA封装焊点的MTTF为779次温度循环。
采用第2节提出的方法对样品进行了温度循环试验，并使用Weibull分布对实验数据进行拟合分析，预测案例的焊点失效概率图见图14。分析得到该产品的关键焊点MTTF为875，其95%置信区间为(725, 1 057)。预测值落在了95%置信区间内，略低于估计值。验证结果表明，所建立的寿命评估方法具有较高的可行性。
图14 预测案例的焊点寿命验证结果
4 结论
研究团队采用有限元法对CBGA焊点的寿命进行了深入预测分析，为电子封装领域的可靠性设计提供了坚实的理论基础。通过综合仿真分析和实验数据，成功识别出焊点平均塑性应变能密度增量ΔWave作为关键的可靠性表征因子，并据此优化了模型的稳健性。通过将实验实测数据与仿真结果有机融合，并集成至Darveaux寿命模型中，建立了一套有效的焊点寿命预测方法。研究结果表明，该方法能够准确预测CBGA焊点的寿命，为同类封装产品的可靠性设计和评估提供了有力的技术支持。
在实际应用中，多种因素如植球工艺、基板表面焊盘强度、焊球材料本构准确性等均可能影响焊点寿命的预测准确性。因此，在试验标定过程中应尽量采用与预测对象相同的工艺条件和试验环境，以获得更高的预测准确度。同时，这些影响因素的作用机理也有待未来进一步研究，以推动CBGA焊点寿命预测技术的发展。
参考文献：
[1] 朱淼, 杨雪霞, 孙艳玺, 等. 热循环条件下BGA叠层焊点可靠性研究及参数优化[J]. 电子元件与材料, 2024, 43(7): 892-898.
[2] 张惟斌, 申坤, 姚颂禹, 等. 基于Anand模型的BGA封装热冲击循环分析及焊点疲劳寿命预测[J]. 电子与封装, 2024, 24(3): 030207.
[3] LIM N R E G, UBANDO A T, GONZAGA J A, et al. Finite element analysis on the factors affecting die crack propagation in BGA under thermo-mechanical loading[J]. Engineering Failure Analysis, 2020, 116: 104717.
[4] XIE D J, ZHANG A, SHIRANGI H, et al. Duplicable and effective: a new drop test for BGA assemblies[C]// 2015 IEEE 65th Electronic Components and Technology Conference (ECTC), San Diego, CA, 2015: 417-424.
[5] MAO M H, WANG W W, LU C H, et al. Machine learning for board-level drop response of BGA packaging structure[J]. Microelectronics Reliability, 2022, 134: 114553.
[6] LIU F, MENG G, ZHAO M, et al. Experimental and numerical analysis of BGA lead-free solder joint reliability under board-level drop impact[J]. Microelectronics Reliability, 2009, 49(1): 79-85.
[7] SYED M M S M Z, ANI F C, RAMLI M R, et al. Effect of moisture content on crack formation during reflow soldering of ball grid array (BGA) component[C]// Advances in Robotics, Automation and Data Analytics. Cham: Springer International Publishing, 2021: 309-314.
[8] WANG J, NIU Y L, SHAO S, et al. A comprehensive solution for modeling moisture induced delamination in electronic packaging during solder reflow[J]. Microelectronics Reliability, 2020, 112: 113791.
[9] HENG E, ABDULLAH M Z. Influence of moisture concentration and hydrophobic material on induced stress in FCBGA package under reflow[J]. Soldering & Surface Mount Technology, 2022, 34(2): 103-115.
[10] 张跃平. BGA封装结构温循载荷下力学响应研究[D]. 北京: 中国运载火箭技术研究院, 2021.
[11] 杨振涛, 余希猛, 张俊, 等. 一种基于倒装芯片的超宽带BGA封装差分传输结构[J]. 半导体技术, 2024, 49(1): 91-96.
[12] 李志丹, 黄勇, 陈世金, 等. 印制电路板BGA焊盘掉点失效案例分析[J]. 印制电路信息, 2015, 23(6): 61-62.
[13] NOURANI A, AKBARI S, SPELT J K. Fracture load prediction of BGA solder joints: Cohesive zone modeling and experimental verification[J]. International Journal of Solids and Structures, 2016, 90: 30-44.
[14] KIM J W, JUNG S B. Experimental and finite element analysis of the shear speed effects on the Sn–Ag and Sn–Ag–Cu BGA solder joints[J]. Materials Science and Engineering: A, 2004, 371(1/2): 267-276.
[15] HU W, LI Y, SUN Y, et al. A model of BGA thermal fatigue life prediction considering load sequence effects[J]. Materials (Basel), 2016, 9(10): E860.
[16] LUO X B, MAO Z M, LIU J, et al. An analytical thermal resistance model for calculating mean die temperature of a typical BGA packaging[J]. Thermochimica Acta, 2011, 512(1/2): 208-216.
[17] 黄春跃, 黄根信, 梁颖, 等. 面向完整传输路径的BGA焊点信号完整性分析及优化[J]. 焊接学报, 2019, 40(3): 25-31.
[18] 夏卓杰, 张亮, 熊明月, 等. 有限元数值模拟在BGA/QFP/CCGA器件焊点可靠性研究中的应用[J]. 电子与封装, 2020, 20(2): 020101.
[19] LIU F, LU Y, WANG Z, et al. Numerical simulation and fatigue life estimation of BGA packages under random vibration loading[J]. Microelectronics Reliability, 2015, 55(12): 2777-2785.
[20] JIANG T B, DU C, XU L H. Finite element analysis and fatigue life prediction of BGA mixed solder joints[C]//High Density Design Packaging and Microsystem Integration, 2007 International Symposium on, 2007: 1-6.
[21] 田艳红, 贺晓斌, 杭春进. 残余应力对混合组装BGA热循环可靠性影响[J]. 机械工程学报, 2014, 50(2): 86-91.
[22] 朱玉丹. BGA焊点的热可靠性研究[D]. 上海: 上海交通大学, 2017.
|
