Anand粘塑性模型的UMAT子程序及验证.docx

1、Anand粘塑性模型的UMAT子程序及验证Anand粘塑性模型的UMAT子程序及验证Anand粘塑性模型的UMAT子程序及验证高军1.引言电子封装及其组件在工艺或者服役过程中, 由于功率耗散和环境温度的周期变化, 会因为电子印制电路板、芯片和焊点的热膨胀失配，在合金钎焊焊点处产生交变的应力应变, 导致焊点的电、热或者机械失效。焊点的热循环失效(可靠性)是电子封装及组装技术中的关键问题之一, 受到了人们的普遍关注。焊点体积细小, 应力应变很复杂。为了准确模拟焊点在服役条件下的应力应变响应, 对可靠性进行评估, 必须建立合理有效的描述钎焊合金材料力学响应的本构方程。SnPb基焊锡钎料广泛应用于电子

2、封装领域，作为电的连接和机械的连接。对于钎料的力学性能的试验和本构模型，许多学者都进行了研究。通常SnPb基焊锡钎料具有很强的温度和加载速率的相关性，应该采用统一型粘塑性本构模型描述SnPb钎料的变形行为。在统一型粘塑性本构模型中，应用最广泛的是Anand模型。具有形式简单，模型参数少等特点，在电子焊点的寿命预测中广泛应用。它采用与位错密度、固溶体强化以及晶粒尺寸效应等相关的单一内部变量S描述材料内部状态对塑性流动的宏观阻抗，可以反映粘塑性材料与应变速度、温度相关的变形行为，以及应变率的历史效应、应变硬化和动态回复等特征。 目前，很多大型商用有限元软件，如ANSYS、MARC等都把Anand本

3、构模型嵌入到通用材料模型库中供用户使用，但是，ABAQUS的通用材料模型库中缺少Anand模型。因此，本报告目的在于通过ABAQUS的用户子程序接口UMAT，选择合适的算法，将Anand粘塑性本构模型引入ABAQUS中，以便后续的研究。2.Anand本构方程统一型粘塑性Anand本构模型有两个基本特征：(1) 在应力空间没有明确的屈服面, 故在变形过程中不需要加载/卸载准则, 塑性变形在所有非零应力条件下产生。(2) 采用单一内部变量描述材料内部状态对塑性流动的宏观阻抗。内部变量(或称变形阻抗) 用S标记, 具有应力量纲。粘塑性Anand模型的流动方程采用双曲蠕变规律对材料的率相关性与温度相关

4、性进行预测，如下式：式中，为Cauchy应力，为偏应力, 为弹性张量，为总应变, 为热膨胀系数，为非弹性应变速率，A为常数，Q为激活能，m为应变敏感指数，为应力乘子，为气体常数， T为绝对温度，T0为参考温度，h0为形变硬化-软化常数，a为与硬化-软化相关的应变敏感系数，S*为变量饱和值，为系数，n为指数。粘塑性Anand本构方程中，共有9个材料参数:A, Q, ,m, n, h0，,a以及初始形变阻抗S0。为真实模拟钎焊材料内部损伤变化，引入损通迭代和弦位法，进行试算比较方法的优劣。3.1 向后隐式Euler方法+普通迭代向后隐式Euler计算公式为：向后Euler方法是隐式方法，计算时要解

5、隐式方程，通常用到迭代法。例如，先用向前显式Euler方法的计算结果作为初值，再作迭代，计算格式为：普通迭代的格式为：，判别迭代过程收敛的条件为：采用上述算法，（1）-（5）式数值计算的离散格式可以表述如下：根据上述计算格式，UMAT子程序的计算流程为：（1） 读取由ABAQUS传递给UMAT子程序的 和，作为计算初值； （2）采用迭代法，联立方程（6）-（10）式，求解；（3）更新应力及全部状态变量，更新Jacobian矩阵。其中，迭代法的计算流程具体如下：（1）迭代循环开始，针对赋计算初值；（2）将方程（6）-（9）式写成形如的迭代格式，由第k步的及计算第k+1步的及；（3）分析比较第k步

6、与第k+1步的及，若它们之间的差满足精度要求，结束循环；否则，继续循环。若循环次数大于预定最大循环次数时，迭代失败。向后Euler方法具有绝对数值稳定性，误差具有一阶精度。虽然是绝对稳定的，但是迭代步长仍要受到一定限制。3.2 中点法+普通迭代中点法计算公式为：中点法也是隐式方法，计算时要用到迭代法解隐式方程。先用向前显式Euler方法的计算结果作为初值，再作迭代，同样采用普通迭代方法，计算格式为：采用上述算法，（1）-（5）式数值计算的离散格式可以表述如下：根据上述计算格式，UMAT子程序的计算流程为：（1） 读取由ABAQUS传递给UMAT子程序的 和，作为计算初值，并计算。（2）采用迭代

7、法，联立方程（11）-（18）式，求解；（3）更新应力及全部状态变量，更新Jacobian矩阵。其中，迭代法的计算流程具体如下：（1）迭代循环开始，针对赋计算初值；（2）将方程（11）（12）（14）（16）式写成形如的迭代格式，由第k步的及计算第k+1步的及；（3）分析比较第k步与第k+1步的及，若它们之间的差满足精度要求，结束循环；否则，继续循环。若循环次数大于预定最大循环次数时，迭代失败。中点法也具有绝对数值稳定性。它的迭代收敛条件比用向后Euler方法的迭代步长可以大一倍，误差具有二阶精度，比向后Euler高一阶。但中点法每积分一步需要计算两次函数值，精度的提高时以增加计算量为代价的。

8、3.3 中点法+弦位法迭代数值计算同样采用中点法，同3.2中的离散格式（11）-（18）。弦位法迭代格式为： 根据弦位法迭代格式，UMAT子程序的迭代格式表示为：UMAT子程序中迭代的计算流程为：（1）迭代循环开始，针对赋计算初值；（2）将方程（11）（12）（14）（16）式写成形如式（19）-（22）的迭代格式，由第k步的及计算，再利用式（20）-（22）计算第k+1步的及；（3）分析第k+1步的，若它的值满足精度要求，结束循环；否则，继续循环。若循环次数大于预定最大循环次数时，迭代失败。弦位法比普通迭代收敛得快，但是计算工作量要增多，需要计算前两步的值。具体的优劣要针对不同的模型来定。3

9、.4 应力、应变增量的Jacobian矩阵 采用隐式数值算法，在UMAT子程序中，需要更新应力、应变增量的Jacobian矩阵。推导后得到一致的应力、应变增量的Jacobian矩阵，如下式：4. 单元测试：根据上述算法和流程，编写了Anand粘塑性本构模型的UMAT子程序，分别为：KANAND. FOR（向后隐式Euler方法+普通迭代）、KMIDDLE. FOR（中点法+普通迭代）、KMIDDLE-XUE. FOR（中点法+弦位法迭代）。三个程序具有相同的状态变量，共8个，即非弹性应变张量的6个分量、非弹性形变阻抗S和损伤值D。程序中涉及15个相同的材料参数，它们依次为：Youngs模量，P

10、ossion比，热膨胀系数，参考温度，初始损伤，损伤阈值以及式（2）-（5）中的参数A, Q, ,m, n, h0, ,a以及初始形变阻抗S0。为了测试程序的正确性，对实际模型进行有限元分析。分析采用三维实体单元，单元类型C3D8：An 8-node linear brick。单位取毫米，吨，秒制。材料参数如下：表1：Pb90Sn10焊料的粘塑性Anand方程的材料参数E/MPa10410-5A/s-1, 107QJ/molmh0 MPaMPan10-3aS0 MPaD0D阈值1.620.382.782933.251558370.143178772.734.383.7315.09暂取0.08暂

11、取1三个程序中，程序KMIDDLE-XUE. FOR（中点法+弦位法迭代）暂有收敛问题未解决，故本报告只进行前两程序的测试。分别进行正方体模型和长方体模型测试。4.1 正方体测试正方体模型边长取0.02，位移加载2E-007，时间5s。经测试，所编程序在拉伸、剪切和热加载测试中具有统一的正确性，即程序如果可以进行一种载荷，就可以进行其他载荷加载，反之不能加载。故本报告主要利用拉伸载荷来测试程序在不同网格数下的正确性。图1是网格边13个种子的正方体模型单向拉伸时得到的变形情况。图1 网格边13个种子的正方体模型 （a）向后Euler方法 （b）中点法图2 应力-应变曲线 Dmax.BACK=8.

12、000126912659419E-002Dmax.MIDDLE=8.000126912659419E-002图3 损伤-时间曲线及损伤云图图2（a）是采用向后Euler方法计算单向拉伸时的应力应变曲线，图2（b）是采用中点法得出的曲线，两者都正确反映了粘塑性材料的单拉变形特征。图3是采用两种方法得出的损伤-时间曲线和损伤云图，因完全相同只列出一个，从输出数据读出损伤最大值如图中所示，数值相同。从以上分析可以看出采用两种方法在图1所示网格数下可以计算，并得出完全相同结果。图4是网格边15个种子的正方体模型单向拉伸时得到的变形情况。图4 网格边15个种子的正方体模型 （a）向后Euler方法 （b

13、）中点法图5 应力-应变曲线 Dmax.BACK = 8.000126811290798E-002图6 向后Euler方法计算的损伤-时间曲线及损伤云图图5（a）是采用向后Euler方法计算单向拉伸时的应力应变曲线，图2（b）是采用中点法得出的曲线，后者在计算中出现错误停止，图中为计算完成了的曲线。对两者已加载的应变部分比较，得出的曲线相近。图6是采用Euler方法得出的损伤-时间曲线和损伤云图，从输出数据读出损伤最大值如图中所示，与图3中得到的数值相近。为找出中点法在计算中出现错误的原因和出现错误时模型的状态，绘出以下图示。 （a）E-TIME （b）E: Max.Principal-TIM

14、E （c）S and S.Mises-TIME （d）S: Max.Principal-TIME（e）损伤值图7 中点法计算结果 从图7（b）和图7（d）看出在计算的最后最大主应力和最大主应变突变为零，导致出现图7（e）的情况，损伤值突变为无穷大。经过计算，单元边种子数13以下向后欧拉和中点法都可以计算，15以上只有向后欧拉可以继续计算，中点法计算中某些变量出现突变，导致无法继续计算。具体原因需要进一步研究，可能是和算法的具体格式有关。4.2 长方体测试正方体模型边长取0.02*0.02*0.072，位移加载单向拉伸，时间5s，网格种子9*9*13，采用中点法程序。本部分测试在相同的网格划分条

15、件下，单向拉伸位移加载的极限。采用位移加载2.88E-5，相应长度0.072为0.04%， 图8是长方体模型单向拉伸时得到的变形情况。 图8 单向拉伸变形 图9 应力-应变曲线 Dmax.MIDDLE= 8.193396334040590E-002图10 损伤-时间曲线 图11 损伤云图图9是采用中点法计算单向拉伸时的应力应变曲线。图10是计算中损伤变化的曲线，可以看出损伤有一定得增长。图11是损伤云图。所以在应变为0.04%时候采用中点法可以很好的模拟。 采用位移加载3.24E-5，相应长度0.072为0.045%，计算未完成。计算到一定时间，最大主应变突变为零，计算出错停止。从图12看出，已计算部分，应力、应变均为直线，处于线性阶段，未进入塑性阶段。 图12 最大主应变-时间曲线经过计算，位移加载在0.04%以下中点法均可以计算，这个值以上计算中某些变量出现突变，导致无法继续计算。5. 结论：计算过程中，到达某个点之后应力、应变中的某个量突变为零，这个点不是应力