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

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

------------------------------------------作者xxxx ------------------------------------------日期xxxx

……………………………………………………………最新资料推荐…………………………………………………

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

高军

1。引言

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

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

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

- 1 - / 18- 1 -

……………………………………………………………最新资料推荐…………………………………………………

目前,很多大型商用有限元软件,如ANＳYS、ＭARC等都把Ａnａnd本构模型嵌入到通用材料模型库中供用户使用，但是,ABAＱUS的通用材料模型库中缺少Ananｄ模型。因此，本报告目的在于通过ABAQUＳ的用户子程序接口ＵMAＴ,选择合适的算法,将Anａnd粘塑性本构模型引入ＡBAQUS中,以便后续的研究。

2.Anand本构方程

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

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

C:[pI(TT)]0Q()3pAeRTsinh(21m3s:s2)Ss3s:s2

SSh1Ssign(1)0S*S*n2Qp:p3S*SeRTA2p:p3 - 2 - / 18- 2 -

……………………………………………………………最新资料推荐…………………………………………………

式中，为Cauｃhy应力,s为偏应力, C为弹性张量，为总应变, 为热膨胀系数,p为非弹性应变速率，A为常数,Q为激活能,m为应变敏感指数,为应力乘子,R为气体常数, Ｔ为绝对温度,T0为参考温度，h0为形变硬化—软化常数，ａ为与硬化—软化相关的应变敏感系数，Ｓ

*

为变量饱和值，S为系数，ｎ为指数.

粘塑性Aｎａnｄ本构方程中，共有9个材料参数:A， Ｑ, ,m, n, ｈ0,S,a以及初始形变阻抗S0。

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

[I(TT)]:C:[I(TT)] Dpp00加入损伤的Aｎand模型方程如下：

(1D)C:[pI(TT)]0D[pI(TT)]:C:[pI(TT)]00Q()3pAeRTsinh(2Sh1Ssign(10S*S*S2Qp:p3eRTA1m3s:ss2)S3s:s2S2p:p)*3Sn(1)(2)(3)(4)

(5)

3．算法与计算流程

计算(2)，(3),(４）式，主要有三种数值算法,向前显式Euler方法，向后隐式Eｕler方法,中点法.向前显式Euleｒ方法是条件稳定，具有一阶精度;向后隐式Euｌer方法和中点法是绝对稳定，分别具有一阶和二阶精度，它们均为隐式方法,需要利用迭代解隐式

- 3 - / 18- 3 -

……………………………………………………………最新资料推荐…………………………………………………

方程。迭代主要有四种方法：普通迭代,牛顿法，弦位法，抛物线法．

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

3。1 向后隐式Eｕler方法＋普通迭代

向后隐式Ｅuｌｅr计算公式为：

y(x)y00 yyhf(x,y)n1nn1n1向后Ｅuler方法是隐式方法，计算yn1时要解隐式方程，通常用到迭代法。例如，先用向前显式Eｕler方法的计算结果作为初值，再作迭代,计算格式为：

(0)yhf(x,y)n1nnn (k1)(k)yyhf(x,y)n1nn1n1y普通迭代的格式为:xk1(xk),判别迭代过程收敛的条件为：

x(k1)(k)x n1n1采用上述算法，(1）-(5)式数值计算的离散格式可以表述如下:

- 4 - / 18- 4 -

……………………………………………………………最新资料推荐…………………………………………………

(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(n1)(ppp(n1)S*S3eRTtA(n1)（6）（7）（8）（9）

（10）根据上述计算格式,UMAT子程序的计算流程为： (1） 读取由ABAQUS传递给UＭAＴ子程序的

n)(n)(n1) (n),(n),,D(n)，(,Tp,S和T,作为计算初值;

(２）采用迭代法,联立方程(6）—（10)式,求解

n1)(n1); (n1),D(n1)，(及Sp(3)更新应力及全部状态变量，更新Jaｃoｂiａn矩阵。 其中,迭代法的计算流程具体如下: （1）迭代循环开始,针对(n1),Dn1)(n1)(n1)赋计算初值； ，(及Sp（2)将方程(6）-(9）式写成形如yk1f(xk,yk)的迭代格式,由第k

步的(n1),D(n1)(n1)及S(n1)计算第ｋ＋1，p步的

n1)及S(n1),D(n1)，(p(n1)；

Dn1)及S(n1),若(n1)，(p（3)分析比较第ｋ步与第ｋ+1步的(n1),它们之间的差满足精度要求，结束循环;否则，继续循环。若循环次数大于预定最大循环次数时，迭代失败.

- 5 - / 18- 5 -

……………………………………………………………最新资料推荐…………………………………………………

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

３.２ 中点法+普通迭代

中点法计算公式为：

y(x)y00yhy(f(x,y)f(x,y))n1n2nnn1n1

中点法也是隐式方法,计算yn1时要用到迭代法解隐式方程。先用向前显式Ｅｕler方法的计算结果作为初值，再作迭代，同样采用普通迭代方法,计算格式为:

(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)I(T(n)T)]f[(n)(T)]:C:[(n)(pI(Tp00n)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)S(n)hh1sign(10*(n)S(n)*SS13(n)(n)ms:s2)(n)S(11)(12)(13)(14)s(n)3(n)(n)s:s2(15)(16)(n)S)*(n)S2(n1)n)):((n1)(n))(p(ppp3n(18)(17)

Q2(n1)n)):((n1)(n))(p(ppp(n)3eRTtA根据上述计算格式，UＭAＴ子程序的计算流程为：

- 6 - / 18- 6 -

……………………………………………………………最新资料推荐…………………………………………………

(1) 读取由ABAＱUS传递给UＭAT子程序的

n)(n)(n1) 和T(n),(n),,D(n)，(,Tp,S,作为计算初值，并计

算f(n),g(n),h(n).

(２）采用迭代法,联立方程(１1)-(18)式,求解

n1)(n1); (n1),D(n1)，(及Sp（3）更新应力及全部状态变量,更新Jaｃobian矩阵。 其中,迭代法的计算流程具体如下: （1)迭代循环开始，针对(n1),Dn1)(n1)(n1)，(及Sp赋计算初值;

(２）将方程（11）（1２)（14)(1６)式写成形如yk1f(xk,yk)的迭

代格式，由第k步的(n1),的(n1),Dn1)(n1)，(pDn1)(n1)，(p及S(n1)计算第k+1步

及S(n1);

Dn1)(n1)，(p(3）分析比较第k步与第k+１步的(n1),及S(n1),若它

们之间的差满足精度要求,结束循环;否则,继续循环。若循环次数大于预定最大循环次数时，迭代失败。

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

3．3 中点法+弦位法迭代

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

- 7 - / 18- 7 -

……………………………………………………………最新资料推荐…………………………………………………

f(x)kxx(xx)

k1kf(x)f(xkk1)kk1根据弦位法迭代格式,UＭＡT子程序的迭代格式表示为:

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子程序中迭代的计算流程为: （1）迭代循环开始,针对(n1),Dn1)(n1)(n1)赋计算初值； ，(及Sp(2)将方程（11）(１2）（１4）（16）式写成形如式(１9）—(2２)

的迭代格式，由第k步的(n1),Dn1)(n1)，(p及S(n1)计算

p(σ)，再利用式(2０)-(22）计算第kn1)及S(n1),D(n1)，(pk+1步的

(n1);

(3)分析第k+1步的p(σk1),若它的值满足精度要求,结束循环;否则,

继续循环。若循环次数大于预定最大循环次数时,迭代失败。

弦位法比普通迭代收敛得快，但是计算工作量要增多,需要计算前两步的值。具体的优劣要针对不同的模型来定。

3。4 应力、应变增量的Jaｃobｉａn矩阵

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

- 8 - / 18- 8 -

……………………………………………………………最新资料推荐…………………………………………………

(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4343000000000000000000)100020002000000000000000000

４. 单元测试：

根据上述算法和流程,编写了Aｎand粘塑性本构模型的UMAT子程序，分别为：KANAND. ＦOR(向后隐式Ｅｕlｅr方法+普通迭代)、ＫMIＤDLE。 FOR(中点法+普通迭代)、KＭIDＤＬE-XＵE。 FOR(中点法+弦位法迭代）.

三个程序具有相同的状态变量,共8个,即非弹性应变张量的6个分量、非弹性形变阻抗S和损伤值D。程序中涉及１５个相同的材料参数,它们依次为:Yｏunｇ’s模量，Posｓiｏn比，热膨胀系数,参考温度，初始损伤，损伤阈值以及式(2)—（５)中的参数A, Q, ,m， n, h０,

S,a以及初始形变阻抗S0.

为了测试程序的正确性,对实际模型进行有限元分析．分析采用三维实体单元,单元类型C3Ｄ8：An 8—ｎode lineaｒ ｂrｉｃk.单位取毫米,吨，秒制。材料参数如下:

- 9 - / 18- 9 -

……………………………………………………………最新资料推荐…………………………………………………

表1：Ｐb9０Sn１０焊料的粘塑性Ananｄ方程的材料参数

Ｑ n E/MPa A/s—１, J/ h0 ×S０ MPS T/K×10 m a D０ 0７ —MPa ×104 ×1０ｍMPa 10ａ —５ ３ ol 1暂５4。1．6０。1772.15．０取2.78 293 ３．25 57 0.143 ３3.73 ２ 38 ８7 ７3 9 0。８８ 08 3 /K1 D 阈值 暂取1 三个程序中,程序KMIDDLE—XＵE。 FOR（中点法+弦位法迭代）暂有收敛问题未解决,故本报告只进行前两程序的测试。分别进行正方体模型和长方体模型测试。

4。1 正方体测试

正方体模型边长取０。02,位移加载2Ｅ—０07，时间5s。经测试,所编程序在拉伸、剪切和热加载测试中具有统一的正确性,即程序如果可以进行一种载荷，就可以进行其他载荷加载，反之不能加载。故本报告主要利用拉伸载荷来测试程序在不同网格数下的正确性。

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

- 10 - / 18- 10 -

……………………………………………………………最新资料推荐…………………………………………………

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

（a）向后Eｕleｒ方法

(ｂ)中点法

图２ 应力-应变曲线

Dｍａx。BＡCK=8。941９E—002 Dmaｘ。MIDDＬＥ=８。9４１9Ｅ—０02

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

- 11 - / 18- 11 -

……………………………………………………………最新资料推荐…………………………………………………

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

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

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

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

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

（a)向后Ｅuler方法

(ｂ）中点法

- 12 - / 18- 12 -

……………………………………………………………最新资料推荐…………………………………………………

图5 应力—应变曲线

Dｍaｘ。BＡＣK = 8。００98Ｅ—002

图６ 向后Eｕlｅｒ方法计算的损伤－时间曲线及损伤云图

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

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

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

- 13 - / 18- 13 -

……………………………………………………………最新资料推荐…………………………………………………

(a）E—TIME （b）Ｅ： Max。Princiｐal—TIMＥ

(ｃ)S ａnd Ｓ。Mises-TIＭE （d)S： Max。Priｎciｐal-TIＭE

(ｅ）损伤值

图7 中点法计算结果

从图7(ｂ)和图7(d)看出在计算的最后最大主应力和最大主应变突变为零,导致出现图7（e)的情况,损伤值突变为无穷大。

经过计算,单元边种子数1３以下向后欧拉和中点法都可以计算,15以上只有向后欧拉可以继续计算，中点法计算中某些变量出现突变，导致无法继续计算。具体原因需要进一步研究，可能是和算法的具体格式有关。

４.２ 长方体测试

- 14 - / 18- 14 -

……………………………………………………………最新资料推荐…………………………………………………

正方体模型边长取0．０2*０.02＊0.0７２,位移加载单向拉伸，时间5s,网格种子９*9*13，采用中点法程序。本部分测试在相同的网格划分条件下，单向拉伸位移加载的极限．

采用位移加载2。88Ｅ-5，相应长度0．07２为0。04%, 图8是长方体模型单向拉伸时得到的变形情况．

图8 单向拉伸变形 图9 应力—应

变曲线

Ｄmax。MＩDDLＥ= ８。1933９6３340４0５90E-002

图1０ 损伤—时间曲线 图11 损伤云图

图9是采用中点法计算单向拉伸时的应力应变曲线。图1０是计算中损伤变化的曲线，可以看出损伤有一定得增长。图１1是损伤云图。所以在应变为0。０4%时候采用中点法可以很好的模拟。

- 15 - / 18- 15 -

……………………………………………………………最新资料推荐…………………………………………………

采用位移加载3．24E-5,相应长度0.07２为0。０4５％,计算未完成。计算到一定时间,最大主应变突变为零，计算出错停止。从图12看出，已计算部分，应力、应变均为直线,处于线性阶段,未进入塑性阶段.

图12 最大主应变-时间曲线

经过计算，位移加载在0。０4%以下中点法均可以计算,这个值以上计算中某些变量出现突变,导致无法继续计算。

5。 结论:

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

参考文献:

［1］ 黄再兴,针对Anaｎd粘塑性模型的UMＡＴ子程序的算法选

择

- 16 - / 18- 16 -

……………………………………………………………最新资料推荐…………………………………………………

[2] 王国忠，程兆年,SnPb 钎料合金的粘塑性Ａｎaｎd本构方程，

应用力学学报,Ｖｏｌ。1７（3),20０0

［3] 张莉,陈旭,Nose H and Sakaｎe Ｍ, Ａｎａnd模型预测6３S

ｎ37Pｂ焊锡钎料的应力应变行为，机械强度,Ｖol。2６（4），2004,pp447-45０

[４] 朱奇农，王国忠，程兆年,罗乐,复合SnPb焊点的形态与可靠

性预测，金属学报,Vol。36（1)，２000,pp93—9８ ［5］ 徐萃薇,计算方法引论，高等教育出版社,1985

- 17 - / 18- 17 -