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

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

高军

1.引言

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

SnPb基焊锡钎料广泛应用于电子封装领域，作为电的连接和机械的连接。对于钎料的力学性能的试验和本构模型，许多学者都进行了研究。通常SnPb基焊锡钎料具有很强的温度和加载速率的相关性，应该采用统一型粘塑性本构模型描述SnPb钎料的变形行为。

在统一型粘塑性本构模型中，应用最广泛的是Anand模型。具有形式简单，模型参数少等特点，在电子焊点的寿命预测中广泛应用。它采用与位错密度、固溶体强化以及晶粒尺寸效应等相关的单一内部变量S描述材料内部状态对塑性流动的宏观阻抗，可以反映粘塑性材料与应变速度、温度相关的变形行为，以及应变率的历史效应、应变硬化和动态回复等特征。

目前，很多大型商用有限元软件，如ANSYS、MARC等都把Anand本构模型嵌入到通用材料模型库中供用户使用，但是，ABAQUS的通用材料模型库中缺少Anand模型。因此，本报告目的在于通过ABAQUS的用户子程序接口UMAT，选择合适的算法，将Anand粘塑性本构模型引入ABAQUS中，以便后续的研究。

2.Anand本构方程

统一型粘塑性Anand本构模型有两个基本特征：(1) 在应力空间没有明确的屈服面, 故在变形过程中不需要加载/卸载准则, 塑性变形在所有非零应力条件下产生。(2) 采用单一内部变量描述材料内部状态对塑性流动的宏观阻抗。内部变量(或称变形阻抗) 用S标记, 具有应力量纲。

粘塑性Anand模型的流动方程采用双曲蠕变规律对材料的率相关性与温度相关性进行预测，如下式：

??C:[???p?I?(T?T)]01?m3Q?s:s?(?)??3s?p?AeRT?sinh(?2?)??2S??3s:s ??2?????SS2?????h1??p:??pSsign(1?)???0**3SS????n?2?Q?p:??p?????3*S?S?eRT?A??????

- 1 -

式中，?为Cauchy应力，s为偏应力, C为弹性张量，?为总应变, ?为热膨胀系数，

?p为非弹性应变速率，A为常数，Q为激活能，m为应变敏感指数，?为应力乘子，R为气??体常数， T为绝对温度，T0为参考温度，h0为形变硬化-软化常数，a为与硬化-软化相关的应

?变敏感系数，S为变量饱和值，S为系数，n为指数。

*

?粘塑性Anand本构方程中，共有9个材料参数:A, Q, ?,m, n, h0，S,a以及初始形变阻抗

S0。

为真实模拟钎焊材料内部损伤变化，引入损伤，演化率如下式：

??[????I?(T?T)]:C:[????I?(T?T)] Dpp00加入损伤的Anand模型方程如下：

??(1?D)?C:[???p?I?(T?T)]0??[????I?(T?T)]:C:[????I?(T?T)]Dpp00Q?(??)?3?p?AeRT?sinh(??2???1?m3s:s?2)??S???(1)(2)s3s:s2(3)???SS?????h1?Ssign(1?)??0?S*S*???n?2?Q?p:??p????3?*S?S?eRT?A??????2?:??p?3p(4)

(5)

3.算法与计算流程

计算（2），（3），（4）式，主要有三种数值算法，向前显式Euler方法，向后隐式Euler方法，中点法。向前显式Euler方法是条件稳定，具有一阶精度；向后隐式Euler方法和中点法是绝对稳定，分别具有一阶和二阶精度，它们均为隐式方法，需要利用迭代解隐式方程。迭代主要有四种方法：普通迭代，牛顿法，弦位法，抛物线法。

本报告中主要采用数值绝对稳定的向后Euler方法和中点法两种数值算法，迭代采用普通迭代和弦位法，进行试算比较方法的优劣。

3.1 向后隐式Euler方法+普通迭代

向后隐式Euler计算公式为：

y(x)?y00 y?y?hf(x,y)n?1nn?1n?1向后Euler方法是隐式方法，计算y

n?1时要解隐式方程，通常用到迭代法。例如，先用向

- 2 -

前显式Euler方法的计算结果作为初值，再作迭代，计算格式为：

(0)?y?hf(x,y)n?1nnn (k?1)(k)y?y?hf(x,y)n?1nn?1n?1y普通迭代的格式为：x(k?1)(k)?x?? ??(x)，判别迭代过程收敛的条件为：xk?1kn?1n?1采用上述算法，（1）-（5）式数值计算的离散格式可以表述如下：

?(n?1)?(1?Dn?1)?I?(T(n?1)?T)])C:[?(n?1)??(p0(n?1)(n)n?1)?I?(T(n?1)?T)]:C:[?(n?1)??(n?1)?I?(T(n?1)?T)]D?D??t?[?(n?1)??(pp001Q?3(n?1)(n?1)?m(??)?s:s?(n?1)3s(n?1)(n?1)2(n)RT?????p??t?Aesinh(?)?(n?1)2??3(n?1)(n?1)pSs:s??2???(n?1)?(n?1)??SS?2(n?1)(n?1)(n)n)):(?(n?1)??(n))S?S??h1?sign(1?)??(?p??(ppp0*(n?1)*(n?1)?3?SS??nQ?2(n?1)?(n)(n?1)(n)??p):(?p??p)?(n?1)?(n?1)??3(?p*?S?S?eRT?tA??????(n?1)（6）（7）（8）（9）

（10）根据上述计算格式，UMAT子程序的计算流程为：

n)(n)(n?1) （1） 读取由ABAQUS传递给UMAT子程序的?(n),?(n),??,D(n)，?(,S,Tp和?T，作为计算初值；

n?1)(n?1)；（2）采用迭代法，联立方程（6）-（10）式，求解?(n?1),D(n?1)， ?(及Sp（3）更新应力及全部状态变量，更新Jacobian矩阵。

其中，迭代法的计算流程具体如下：

n?1)(n?1)赋计算初值； （1）迭代循环开始，针对?(n?1),D(n?1)，?(及Sp（2）将方程（6）-（9）式写成形如yk?1?f(xk,yk)的迭代格式，由第k步的

?(n?1),D(n?1)及S(n?1)计算第k+1步的(n?1)(n?1)n?1)及(n?1)，?p?,D，?(pS(n?1)；

n?1)及S(n?1)，若它们之间的差（3）分析比较第k步与第k+1步的?(n?1),D(n?1)，?(p满足精度要求，结束循环；否则，继续循环。若循环次数大于预定最大循环次数时，

迭代失败。

向后Euler方法具有绝对数值稳定性，误差具有一阶精度。虽然是绝对稳定的，但是迭代步长仍要受到一定限制。

- 3 -

3.2 中点法+普通迭代

中点法计算公式为：

y(x)?y00yh?y?(f(x,y)?f(x,y))n?1n2nnn?1n?1

中点法也是隐式方法，计算yn?1时要用到迭代法解隐式方程。先用向前显式Euler方法的

计算结果作为初值，再作迭代，同样采用普通迭代方法，计算格式为：

(0)?y?hf(x,y)n?1nnn

h(k?1)(k)y?y?(f(x,y)?f(x,y)n?1n2nnn?1n?1y采用上述算法，（1）-（5）式数值计算的离散格式可以表述如下：

n?1)?I?(T(n?1)?T)]?(n?1)?(1?D(n?1))C:[?(n?1)??(p0?t(n)(n?1)(f?f)2(n)(n)(n)n)n)f?[?(n)??(?T)]:C:[?(n)??(?T)]p?I?(Tp?I?(T00n)?t(n)?g(n?1))?(n?1)??(p?2(gpD(n?1)?D(n)?Q?(??)?(n)3?sinh(?g(n)?AeRT2???(n?1)(n)?t(n)(n?1)S?S?(h?h)2??(n)(n)?SS?(n)?h??h1?sign(1?)??0*(n)*(n)??SS?????(n)S*?S????13(n)(n)?ms:s?2)??(n)?S??(11)(12)(13)(14)s(n)3(n)(n)s:s2(15)(16)2(n?1)n)):(?(n?1)??(n))(?p??(ppp3??????n(18)(17)

Q2(n?1)n)):(?(n?1)??(n))(?p??(?(n)ppp3eRT?tA根据上述计算格式，UMAT子程序的计算流程为：

n)(n)(n?1) （1） 读取由ABAQUS传递给UMAT子程序的?(n),?(n),??,D(n)，?(,Tp,S和?T，作为计算初值，并计算

f(n),g(n),h(n)。

n?1)(n?1)；（2）采用迭代法，联立方程（11）-（18）式，求解?(n?1),D(n?1)， ?(及Sp（3）更新应力及全部状态变量，更新Jacobian矩阵。

其中，迭代法的计算流程具体如下： （1）迭代循环开始，针对?(n?1),Dn?1)(n?1)(n?1)赋计算初值； ，?(及Sp（2）将方程（11）（12）（14）（16）式写成形如yk?1?f(xk,yk)的迭代格式，由第k步的

- 4 -

n?1)及S(n?1)计算第k+1步的(n?1)n?1)及(n?1)?(n?1),D(n?1)，?(?,D，?(ppS(n?1)；

Dn?1)及S(n?1)，若它们之间的差满(n?1)，?(p（3）分析比较第k步与第k+1步的?(n?1),足精度要求，结束循环；否则，继续循环。若循环次数大于预定最大循环次数时，迭

代失败。

中点法也具有绝对数值稳定性。它的迭代收敛条件比用向后Euler方法的迭代步长可以大一倍，误差具有二阶精度，比向后Euler高一阶。但中点法每积分一步需要计算两次函数值，精度的提高时以增加计算量为代价的。

3.3 中点法+弦位法迭代

数值计算同样采用中点法，同3.2中的离散格式（11）-（18）。弦位法迭代格式为：

f(x)kx?x?(x?x)

k?1kf(x)?f(xkk?1)kk?1根据弦位法迭代格式，UMAT子程序的迭代格式表示为：

tryp(σ)?σ?σkkkk(σ)nσ?σ?(σ?σ)k?1kk(σ)?k(σk?1)kkk?1k(σ)nD?D?(D?D)k?1kk(σ)?k(σk?1)kkk?1k(σ)nS?S?(S?S)k?1kk(σ)?k(σk?1)kkk?1(19)(20)

(21)(22)UMAT子程序中迭代的计算流程为：

n?1)(n?1)赋计算初值； （1）迭代循环开始，针对?(n?1),D(n?1)，?(及Sp（2）将方程（11）（12）（14）（16）式写成形如式（19）-（22）的迭代格式，由第k步的

n?1)及S?(n?1),D(n?1)，?(p(n?1)计算

p(σ)，再利用式（20）-（22）计算第k+1

k步的?(n?1),Dn?1)及S(n?1)； (n?1)，?(p（3）分析第k+1步的p(σk?1)，若它的值满足精度要求，结束循环；否则，继续循环。若

循环次数大于预定最大循环次数时，迭代失败。

弦位法比普通迭代收敛得快，但是计算工作量要增多，需要计算前两步的值。具体的优

劣要针对不同的模型来定。

3.4 应力、应变增量的Jacobian矩阵

- 5 -

采用隐式数值算法，在UMAT子程序中，需要更新应力、应变增量的Jacobian矩阵。推导后得到一致的应力、应变增量的Jacobian矩阵，如下式：

?(n?1)????3??sp?(1?D)C??t?(1?D)C:???2?(n?1)???eq?1?0??0?J???1?D?(??0?0??0????2?????????0?0???0001001000000?4??3000??000?????(n?1)??000?3p????t100?2?(n?1)??eq?0010????0001???0??000??000????2?000??0?00?00?0??000????4?3??4?30000000000?000?00?0?0??0??)?10??0??0???????2??000

4. 单元测试：

根据上述算法和流程，编写了Anand粘塑性本构模型的UMAT子程序，分别为：KANAND. FOR（向后隐式Euler方法+普通迭代）、KMIDDLE. FOR（中点法+普通迭代）、KMIDDLE-XUE. FOR（中点法+弦位法迭代）。

三个程序具有相同的状态变量，共8个，即非弹性应变张量的6个分量、非弹性形变阻抗S和损伤值D。程序中涉及15个相同的材料参数，它们依次为：Young’s模量，Possion比，热膨胀系数，参考温度，初始损伤，损伤阈值以及式（2）-（5）中的参数A, Q, ?,m, n, ,a以及初始形变阻抗S0。

为了测试程序的正确性，对实际模型进行有限元分析。分析采用三维实体单元，单元类型C3D8：An 8-node linear brick。单位取毫米，吨，秒制。材料参数如下：

表1：Pb90Sn10焊料的粘塑性Anand方程的材料参数 ??/K?1 A/s-1, Q h0 n S0 S ? ?T/K m a D0 0MPa MPa ×10-3 MPa ×10-5 ×107 J/mol 0.38 2.78 293 3.25 15583 7 0.143 1787 72.73 4.38 3.73 15.09 暂取0.08 h0,

?SE/MPa ×104 1.62 D 阈值 暂取1

三个程序中，程序KMIDDLE-XUE. FOR（中点法+弦位法迭代）暂有收敛问题未解决，故本报告只进行前两程序的测试。分别进行正方体模型和长方体模型测试。

4.1 正方体测试

正方体模型边长取0.02，位移加载2E-007，时间5s。经测试，所编程序在拉伸、剪切和热加载测试中具有统一的正确性，即程序如果可以进行一种载荷，就可以进行其他载荷加

- 6 -

载，反之不能加载。故本报告主要利用拉伸载荷来测试程序在不同网格数下的正确性。

图1是网格边13个种子的正方体模型单向拉伸时得到的变形情况。

图1 网格边13个种子的正方体模型

（a）向后Euler方法 （b）中点法

图2 应力-应变曲线

Dmax.BACK=8.000126912659419E-002 Dmax.MIDDLE=8.000126912659419E-002

图3 损伤-时间曲线及损伤云图

图2（a）是采用向后Euler方法计算单向拉伸时的应力应变曲线，图2（b）是采用中点法得出的曲线，两者都正确反映了粘塑性材料的单拉变形特征。

图3是采用两种方法得出的损伤-时间曲线和损伤云图，因完全相同只列出一个，从输出数据读出损伤最大值如图中所示，数值相同。

从以上分析可以看出采用两种方法在图1所示网格数下可以计算，并得出完全相同结果。

图4是网格边15个种子的正方体模型单向拉伸时得到的变形情况。

- 7 -

图4 网格边15个种子的正方体模型

（a）向后Euler方法 （b）中点法

图5 应力-应变曲线

Dmax.BACK = 8.000126811290798E-002

图6 向后Euler方法计算的损伤-时间曲线及损伤云图

图5（a）是采用向后Euler方法计算单向拉伸时的应力应变曲线，图2（b）是采用中点法得出的曲线，后者在计算中出现错误停止，图中为计算完成了的曲线。对两者已加载的应变部分比较，得出的曲线相近。

图6是采用Euler方法得出的损伤-时间曲线和损伤云图，从输出数据读出损伤最大值如图中所示，与图3中得到的数值相近。

为找出中点法在计算中出现错误的原因和出现错误时模型的状态，绘出以下图示。

- 8 -

（a）E-TIME （b）E: Max.Principal-TIME

（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，采用中点法程序。本部分测试在相同的网格划分条件下，单向拉伸位移加载的极限。

采用位移加载2.88E-5，相应长度0.072为0.04%， 图8是长方体模型单向拉伸时得到的变形情况。

- 9 -

图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. 结论：

- 10 -

计算过程中，到达某个点之后应力、应变中的某个量突变为零，这个点不是应力、应变的一个固定值，具体原因需要进一步研究，可能是和算法的具体格式有关。

参考文献：

[1] 黄再兴，针对Anand粘塑性模型的UMAT子程序的算法选择 [2] 王国忠，程兆年，SnPb 钎料合金的粘塑性Anand本构方程，应用力学学报，Vol.17（3），

2000

[3] 张莉，陈旭，Nose H and Sakane M, Anand模型预测63Sn37Pb焊锡钎料的应力应变行

为，机械强度，Vol.26（4），2004，pp447-450

[4] 朱奇农，王国忠，程兆年，罗乐，复合SnPb焊点的形态与可靠性预测，金属学报，

Vol.36(1),2000,pp93-98

[5] 徐萃薇，计算方法引论，高等教育出版社，1985

- 11 -

计算过程中，到达某个点之后应力、应变中的某个量突变为零，这个点不是应力、应变的一个固定值，具体原因需要进一步研究，可能是和算法的具体格式有关。

参考文献：

[1] 黄再兴，针对Anand粘塑性模型的UMAT子程序的算法选择 [2] 王国忠，程兆年，SnPb 钎料合金的粘塑性Anand本构方程，应用力学学报，Vol.17（3），

2000

[3] 张莉，陈旭，Nose H and Sakane M, Anand模型预测63Sn37Pb焊锡钎料的应力应变行

为，机械强度，Vol.26（4），2004，pp447-450

[4] 朱奇农，王国忠，程兆年，罗乐，复合SnPb焊点的形态与可靠性预测，金属学报，

Vol.36(1),2000,pp93-98

[5] 徐萃薇，计算方法引论，高等教育出版社，1985

- 11 -

本文来源：https://www.bwwdw.com/article/n7uf.html