CN113139309A - 封装单元bga板级互连交变温度载荷下数值仿真方法 - Google Patents

Abstract

本发明公开了一种封装单元BGA板级互连交变温度载荷下数值仿真方法，包括以下步骤：S1、简化并绘制带宽射频单元BGA板级互联三维几何模型；S2、对几何模型中力学行为高保真的非关键区域BAG焊点等效处理；S3、对几何模型设置求解物理场、定义边界条件和温度载荷，并赋值材料属性；S4、进行适用于应力突变条件的结构化网络控制，保证非线性求解的收敛性；S5、通过宏观‑微观多尺度单元迭代实现关键BGA焊点的精确求解；S6、对关键单元误差控制下的特征循环寿命计算。本发明克服了传统温度循环(温度冲击)试验研究评估BGA板级互连可靠性成本高、周期长的缺点，减少了可靠性设计对实物验证的依赖性。

Description

封装单元BGA板级互连交变温度载荷下数值仿真方法

技术领域

本发明属于电子封装板级互联技术领域，具体涉及一种封装单元BGA板级互连交变温度载荷下数值仿真方法。

背景技术

封装单元采用BGA阵列集成实现射频和数字信号的传输，由于板级互连材料体系之间存在的固有热膨胀系数差异，产品在经历交变温度载荷(温度循环或温度冲击)后的由于热应力导致的互连焊点过应力开裂、塑性应变失效或热疲劳失效时有发生，最终导致互连可靠性问题。因此，需要通过正向的理论计算及模拟仿真支撑板级BGA互连的可靠性设计及特征循环寿命预估。

传统的基于试验研究评估BGA互连可靠性的方法周期长、成本高，且重点集中在互连焊点的热疲劳寿命预测上。中国专利CN107203666A公布了一种BGA焊点热疲劳寿命预测方法及系统，通过理论计算及试验测试数据相结合的方法预测疲劳寿命。在前述专利的基础上，中国专利CN104820781B从失效物理技术领域出发，公布了一种不同温度循环载荷顺序加载影响的BGA焊点热疲劳寿命预测方法，考虑了温度载荷的施加顺序带来的影响。上述专利均是采用理论计算公式和试验数据拟合得出最终的热疲劳寿命预估模型，并不涉及阵列BGA板级互连的高精度数值仿真方法。

部分研究文献中，涉及到BGA板级互连温度循环的热力耦合仿真，但由于行业领域的限制，此类仿真并未提出力学行为高保真的BGA焊点等效处理、适用于应力突变下的非线性求解有限元网格控制、宏观-微观多尺度单元迭代精细化求解以及关键单元误差控制下的特征循环寿命计算等方法，导致仿真结果的一致性、稳定性和计算精度无法保证。

发明内容

针对现有技术中的上述不足，本发明提供的一种封装单元BGA板级互连交变温度载荷下数值仿真方法解决了BGA板级互连温度循环的热力耦合仿真时仿真结果的一致性、稳定性和计算精度无法保证的问题。

为了达到上述发明目的，本发明采用的技术方案为：一种封装单元BGA板级互连交变温度载荷下数值仿真方法，包括以下步骤：

S1、简化并绘制带宽射频单元BGA板级互联三维几何模型；

S2、对几何模型中力学行为高保真的非关键区域BAG焊点等效处理；

S3、对几何模型设置求解物理场、定义边界条件和温度载荷，并赋值材料属性；

S4、进行适用于应力突变条件的结构化网络控制，保证非线性求解的收敛性；

S5、通过宏观-微观多尺度单元迭代实现关键BGA焊点的精确求解；

S6、对关键单元误差控制下的特征循环寿命计算。

进一步地：所述步骤S2中BAG焊点等效为圆柱形，其直径和原始BAG焊点的直径相同；通过模拟单个BGA焊点的剪切和拉伸实验，确定合适的弹性模量E^equiv和屈服强度

使得等效焊点和原始焊点在受到剪切力和拉伸力的变形相同，保证等效焊点和初始的复杂结构圆鼓形焊点具有相同的非线性机械属性。

进一步地：所述弹性模量E^equiv及屈服强度

通过多轮模拟迭代调整并确定，使得在相同轴向位移、剪切位移下，等效模型与原始模型之间的轴向力σ_normal、剪切力σ_shear平均误差都小于10％，只在负载最大的四个最外侧角落焊点处采用真实的圆鼓形BGA焊球。

进一步地：述最外侧角落的BGA焊球的网格，邻近焊盘区域的一层单元厚度≤30um。

进一步地：所述步骤S4中结构化网络控制具体为：封装单元BGA板级互连有限元模型中各个部件界面处理为共用面，采用共节点互连方式，避免接触带来的非线性问题；所述金属围框及盖板采用六面体8节点低阶单元，封装单元封装基板、系统母板采用六面体8节点低阶单元，阵列BGA焊球、焊盘采用六面体20节点高阶单元。

进一步地：所述单元的特征值为：单元长宽比平均值≤5，单元最大顶角平均值≤120°，单元雅各比平均值≤2，单元翘曲因子平均值＜0.1。

进一步地：所述步骤S5中宏观-微观多尺度单元迭代为对封装单元BGA板级互连建立宏观粗糙模型，通过力学行为高保真的BGA焊点模型等效方法以及应力突变条件下的结构化网格控制方法完成仿真计算，通过模型切割和边界再生成，将计算出的热变形作为边界条件入到更加详细的单个BGA焊点微观精细模型，实现关键BGA焊点的精细化数值仿真计算。

进一步地：所述单个BGA焊点微观精细有限元模型，经过单元尺寸缩放求解，使得BGA焊点的最终求解结果不依赖于单元尺寸，求解精度和稳定性得到保证，精度不低于90％。

进一步地：所述步骤S6中特征循环寿命计算为对单个BGA焊点微观精细模型的数值仿真结果，根据焊点的应力应变分布梯度，提取靠近焊盘且应力应变最大的一层单元的平均塑性应变能密度ΔW_ave，采用基于能量的Darveaux经验方程计算特征循环寿命；所述一层单元的厚度≤30um。

进一步地：特征循环寿命计算的具体步骤为：

A.对所关注单元建立单元组，对单元组内的各分析单元进行编号，提取选定单元在某个温度循环后累积的塑性应变能密度ΔW_i，并乘以对应单元编号下的单元体积vi，通过线性叠加得到整个单元组的应变能ΔW_layer；

B.通过整个单元组的应变能ΔW_layer除以单元组的总体积V_layer，得到该循环内的单元的累积平均塑性应变能密度；

C.重复步骤A和步骤B，得到第n个稳定循环后的累积平均塑性应变能密度，减去第n-1个稳定循环后的累积平均塑性应变能密度，得到该单元组单个稳定循环内的平均塑性应变能密度ΔW_ave，并采用基于能量的Darveaux经验方程计算特征循环寿命。

本发明的有益效果为：本发明所提出的交变温度载荷下封装单元BGA板级互连数值仿真方法，实现了高效、精准、稳定的仿真过程，克服了传统温度循环(温度冲击)试验研究评估BGA板级互连可靠性成本高、周期长的缺点，减少了可靠性设计对实物验证的依赖性。与现有技术相比，本发明所提出的数值仿真方法有以下优点：

(1)通过力学行为高保真的BGA焊点模型等效方法，模拟单个BGA焊点的剪切和拉伸实验，确定合适的弹性模量E^equiv和屈服强度

保证等效圆柱焊点和初始的复杂结构圆鼓形焊点具有相同的非线性力学行为，将非关键区域BGA采用等效圆柱焊点替代，在保证精度的同时大幅提高了仿真计算效率。

(2)通过宏观-微观多尺度单元迭代法，实现了对危险BGA焊点精准、稳定的数值仿真计算，避免了宏观全模型加密网格所带来的仿真计算量猛增的问题，降低了对仿真计算硬件资源的苛刻要求。

(3)提出了应力突变条件下的结构化网格控制方法，保证非线性迭代求解的收敛性，通过对关键单元计算结果加权平均处理，完成交变温度载荷下互连焊点的特征循环寿命精确测算，避免了实物样件进行温度循环试验耗时长、成本高的问题。

附图说明

图1是本发明的实施流程图；

图2是某封装单元BGA板级互连示意图；

图3是简化后的封装单元BGA板级互连1/4对称三维几何模型图；

图4是仿真时各部件所选材料模型图；

图5是封装单元BGA板级互连1/4对称模型约束施加图；

图6是非关键区域BGA焊点等效过程示意图；

图7是原始焊点和等效焊点轴向响应对比曲线图；

图8是原始焊点和等效焊点切向响应对比曲线图；

图9是某时刻封装单元BGA板级互连BGA焊点的应力分布图；

图10是某时刻封装单元BGA板级互连BGA焊点的应变分布图；

图11是宏观-微观多尺度单元迭代处理过程示意图；

图12是危险BGA焊点靠近封装基板侧的一层单元示意图；

图13是危险BGA焊点等效应力随时间变化曲线图；

图14是危险BGA焊点累积塑性应变随时间变化曲线图；

图15是危险BGA焊点累积塑性应变能密度随时间变化曲线图。

其中：1、金属围框和盖板；2、封装单元基板；3、系统母板；4、系统母板上的焊盘；5、BGA焊球；6、封装单元基板上的焊盘。

具体实施方式

下面对本发明的具体实施方式进行描述，以便于本技术领域的技术人员理解本发明，但应该清楚，本发明不限于具体实施方式的范围，对本技术领域的普通技术人员来讲，只要各种变化在所附的权利要求限定和确定的本发明的精神和范围内，这些变化是显而易见的，一切利用本发明构思的发明创造均在保护之列。

如图1所示，一种封装单元BGA板级互连交变温度载荷下数值仿真方法，包括以下步骤：

S1、简化并绘制带宽射频单元BGA板级互联三维几何模型；

S2、对几何模型中力学行为高保真的非关键区域BAG焊点等效处理；

S3、对几何模型设置求解物理场、定义边界条件和温度载荷，并赋值材料属性；

S4、进行适用于应力突变条件的结构化网络控制，保证非线性求解的收敛性；

S5、通过宏观-微观多尺度单元迭代实现关键BGA焊点的精确求解；

S6、对关键单元误差控制下的特征循环寿命计算。

本发明以一种封装单元BGA板级集成产品为例，用数值仿真方法分析其温度循环试验过程。图1为本发明的实施流程图，本示例的实施采用ANSYS Workbench软件实现，其具体实现过程包含以下步骤：

(1)建立有限元模型

图2是某封装单元BGA板级互连产品的截面示意图。

1)建立封装单元BGA板级互连温度循环仿真的几何模型

A.几何模型简化：图3为简化去除射频封装单元中的芯片、元器件、粘接胶、键合金丝等结构要素后的1/4几何模型，图中包括：1—金属围框和盖板，2—封装单元基板，3—系统母板，4—系统母板上的焊盘，5—阵列BGA焊球，6—封装单元基板上的焊盘。几何模型不考虑制造过程中所造成的残余应力及应变，忽略不同材料间的分层与预裂纹影响。

B.BGA互连区域模型等效处理：本示例中互连BGA焊球直径d＝0.6mm，高度H＝0.3mm，将非关键区域的BGA焊点采用同样直径、同样高度的圆柱代替。图6是非关键区域BGA焊点等效过程示意图。通过模拟单个焊点施加剪切力、拉伸力的过程，使得等效焊点和原始焊点在受到剪切力和拉伸力的变形相同。图7是原始焊点和等效焊点轴向响应对比曲线，图8是原始焊点和等效焊点切向响应对比曲线。通过本实施例，轴向力响应误差控制在7％以内，切向力响应误差控制在10％以内，最终确定等效焊点的弹性模量E^equiv＝86938MPa，屈服强度

2)材料模型构建及求解物理场的设置

A.材料模型构建及设置：在ANSYS Workbench工作台中，调用Engineering data模块，对封装单元封装基板材料和系统母板材料，建立多温度点下的线弹性材料模型；对金属围框和盖板、封装单元基板和系统母板上的焊盘材料，建立多线性随动强化弹塑性材料模型；对阵列BGA焊球，建立多温度点下的粘塑性本构材料模型。材料模型最高温度85℃，最低温度-55℃，温度间隔为15℃。图4为本实施例仿真时各部件所选材料模型。

B.求解物理场设置：设置物理场为“热力耦合”多物理场。

3)定义载荷和边界条件

A.温度载荷：以封装单元BGA板级互连的交变温度载荷试验典型参数为输入，设置参考温度T_ref＝25℃，初始温度T_start＝220℃，回流冷却速率ν＝4℃/s，最高温度T_max＝85℃，最低温度T_min＝-55℃，保温时间t_dwell＝30min，温度转换时间t_ramp＝2min，仿真中共进行5个循环，采用实测温度曲线对交变温度曲线高低温转换拐点邻近区域的具体温度值进行平顺化修正处理。图10是仿真时的温度载荷曲线示意图。

B.约束施加：本实施例中采用1/4等效对称模型，两个对称面施加对称约束，使X向和Y向位移为0；底面中心节点施加Z向约束，使Z向位移为0；保证三个方向旋转自由度受约束。图5是封装单元BGA板级互连1/4对称模型约束施加图。

4)网格划分：图9为封装单元BGA板级互连1/4对称模型有限元网格划分图，其中金属围框及盖板采用六面体8节点低阶单元，网格尺寸1.0mm，封装单元基板、系统母板采用六面体8节点低阶单元，网格尺寸1.0mm，关键区域BGA焊球、焊盘采用六面体20节点高阶单元，网格尺寸0.1mm，非关键区域采用等效圆柱焊点，单元尺寸为0.5mm。总节点数858318个，单元数257161个。为了适用应力突变条件下的非线性求解，本实施例控制单元长宽比平均值4.73，单元最大顶角平均值104°，单元雅各比平均值1.3，单元翘曲因子平均值为0.02。

(2)封装单元BGA板级互连仿真计算及后处理

1)有限元求解：针对热力耦合多物理场分析，采用瞬态牛顿-辛普森(Newton-Raphson Method)非线性迭代计算方法求解，采用ANSYS内置的瞬态直接耦合PCG求解器，通过预置收敛准则控制收敛精度，根据本实施例中的温度载荷曲线，划分载荷区间为24个载荷步，采用固定间距的载荷步长，温度转换过程的载荷步，时间步长为5s，保温过程的载荷步，时间步长为30s。

2)结果后处理：通过后处理提取并查看以下结果

A.BGA焊点的等效应力：添加三维绘图，选择阵列BGA焊点，输入表达式etable,eseqv,s,eqv，并选择时间为17363s，绘制出在温度循环17363s时的BGA焊点等效应力分布图，见图11。可见，最外侧角落处的BGA焊点所受应力最大，且靠近封装单元基板一侧。

B.BGA焊点的等效塑性应变：添加三维绘图，选择阵列BGA焊点，输入表达式etable,epleqv,eppl,eqv，并选择时间为19163s，绘制出在温度循环19163s时的BGA焊点等效塑性应变分布图，见图12。可见，最外侧角落处的BGA焊点所受应变最大，且靠近封装基板一侧。

3)宏观-微观多尺度单元迭代下的BGA焊点精细求解：针对前述的计算结果，通过模型切割和边界再生成，将计算出的热变形(关键焊点邻近区域轴向变形和剪切变形)作为边界条件入到更加详细的单个BGA焊点微观精细模型，实现关键BGA焊点的精细化数值仿真计算。图13是宏观-微观多尺度单元迭代处理过程示意图。进行的单元尺寸无关性验证如下：

A.宏观粗糙模型中相同位置BGA焊点单元尺寸为0.1mm，上下基板的单元尺寸为1.0mm，等效应力26.974MPa，等效塑性应变0.032324；

B.微观BGA焊点精细模型单元尺寸0.1mm，上下基板的单元尺寸为1.0mm，等效应力26.69MPa，等效塑性应变0.032289；

C.微观BGA焊点精细模型单元尺寸0.05mm，上下基板的单元尺寸为0.25mm，等效应力30.166MPa，等效塑性应变0.052132；

D.微观BGA焊点精细模型单元尺寸0.05mm，上下基板的单元尺寸为0.1mm，等效应力30.166MPa，等效塑性应变0.052122；

可见，单元足够精细时，危险区域BGA焊点的求解结果更加一致和稳定，不再依赖单元尺寸。

4)关键单元误差控制下的BGA焊点特征循环寿命计算：

针对3)中切割边界位移处理后的计算结果，采用底层命令提取靠近封装基板焊盘侧的一层单元的加权平均处理结果，图14是危险BGA焊点靠近封装基板侧的一层单元示意。根据前述网格划分的要求，该层单元的厚度为20um。绘制累积塑性应变能密度△W随时间变化的曲线，如图15所示，计算得最后一个温度循环稳定后的平均塑性应变能密度△W_ave＝5.788MPa。

采用基于能量的Darveaux经验方程计算特征循环寿命，方程如下：

特征循环寿命：

其中，a为焊点在界面处的直径，实际产品为0.52mm。K₁～K₄为BGA封装相关系数，K₁＝17.95，K₂＝-1.62，K₃＝1.49*10-3，K₄＝1.04。

带入方程，求得特征循环寿命N_w＝57次。

Claims (10)

1.一种封装单元BGA板级互连交变温度载荷下数值仿真方法，其特征在于，包括以下步骤：

S1、简化并绘制带宽射频单元BGA板级互联三维几何模型；

S2、对几何模型中力学行为高保真的非关键区域BAG焊点等效处理；

S3、对几何模型设置求解物理场、定义边界条件和温度载荷，并赋值材料属性；

S4、进行适用于应力突变条件的结构化网络控制，保证非线性求解的收敛性；

S5、通过宏观-微观多尺度单元迭代实现关键BGA焊点的精确求解；

S6、对关键单元误差控制下的特征循环寿命计算。

2.根据权利要求1所述的封装单元BGA板级互连交变温度载荷下数值仿真方法，其特征在于，所述步骤S2中BAG焊点等效为圆柱形，其直径和原始BAG焊点的直径相同；通过模拟单个BGA焊点的剪切和拉伸实验，确定合适的弹性模量E^equiv和屈服强度

使得等效焊点和原始焊点在受到剪切力和拉伸力的变形相同，保证等效焊点和初始的复杂结构圆鼓形焊点具有相同的非线性机械属性。

3.根据权利要求2所述的封装单元BGA板级互连交变温度载荷下数值仿真方法，其特征在于，所述弹性模量E^equiv及屈服强度

通过多轮模拟迭代调整并确定，使得在相同轴向位移、剪切位移下，等效模型与原始模型之间的轴向力σ_normal、剪切力σ_shear平均误差都小于10％，只在负载最大的四个最外侧角落焊点处采用真实的圆鼓形BGA焊球。

4.根据权利要求3所述的封装单元BGA板级互连交变温度载荷下数值仿真方法，其特征在于，所述最外侧角落的BGA焊球的网格，邻近焊盘区域的一层单元厚度≤30um。

5.根据权利要求1所述的封装单元BGA板级互连交变温度载荷下数值仿真方法，其特征在于，所述步骤S4中结构化网络控制具体为：封装单元BGA板级互连有限元模型中各个部件界面处理为共用面，采用共节点互连方式，避免接触带来的非线性问题；所述金属围框及盖板采用六面体8节点低阶单元，封装单元封装基板、系统母板采用六面体8节点低阶单元，阵列BGA焊球、焊盘采用六面体20节点高阶单元。

6.根据权利要求5所述的封装单元BGA板级互连交变温度载荷下数值仿真方法，其特征在于，所述单元的特征值为：单元长宽比平均值≤5，单元最大顶角平均值≤120°，单元雅各比平均值≤2，单元翘曲因子平均值＜0.1。

7.根据权利要求1所述的封装单元BGA板级互连交变温度载荷下数值仿真方法，其特征在于，所述步骤S5中宏观-微观多尺度单元迭代为对封装单元BGA板级互连建立宏观粗糙模型，通过力学行为高保真的BGA焊点模型等效方法以及应力突变条件下的结构化网格控制方法完成仿真计算，通过模型切割和边界再生成，将计算出的热变形作为边界条件入到更加详细的单个BGA焊点微观精细模型，实现关键BGA焊点的精细化数值仿真计算。

8.根据权利要求7所述的封装单元BGA板级互连交变温度载荷下数值仿真方法，其特征在于，所述单个BGA焊点微观精细有限元模型，经过单元尺寸缩放求解，使得BGA焊点的最终求解结果不依赖于单元尺寸，求解精度和稳定性得到保证，精度不低于90％。

9.根据权利要求7所述的封装单元BGA板级互连交变温度载荷下数值仿真方法，其特征在于，所述步骤S6中特征循环寿命计算为对单个BGA焊点微观精细模型的数值仿真结果，根据焊点的应力应变分布梯度，提取靠近焊盘且应力应变最大的一层单元的平均塑性应变能密度ΔW_ave，采用基于能量的Darveaux经验方程计算特征循环寿命；所述一层单元的厚度≤30um。

10.根据权利要求9所述的封装单元BGA板级互连交变温度载荷下数值仿真方法，其特征在于，所述特征循环寿命计算的具体步骤为：

A.对所关注单元建立单元组，对单元组内的各分析单元进行编号，提取选定单元在某个温度循环后累积的塑性应变能密度ΔW_i，并乘以对应单元编号下的单元体积vi，通过线性叠加得到整个单元组的应变能ΔW_layer；

B.通过整个单元组的应变能ΔW_layer除以单元组的总体积V_layer，得到该循环内的单元的累积平均塑性应变能密度；

C.重复步骤A和步骤B，得到第n个稳定循环后的累积平均塑性应变能密度，减去第n-1个稳定循环后的累积平均塑性应变能密度，得到该单元组单个稳定循环内的平均塑性应变能密度ΔW_ave，并采用基于能量的Darveaux经验方程计算特征循环寿命。

Images