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

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

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

3、模型。因此，本报告目的在于通过ABAQUS的用户子程序接口UMAT,选择合适的算法，将Anand粘塑性本构模型引入ABAQUS中，以便后续的研究。2Anand本构方程统一型粘塑性Anand本构模型有两个基本特征：（1）在应力空间没有明确的屈服面，故在变形过程中不需要加载/卸载准则,塑性变形在所有非零应力条件下产生。（2）采用单一内部变量描述材料内部状态对塑性流动的宏观阻抗。内部变量（或称变形阻抗）用S标记，具有应力量纲。粘塑性Anand模型的流动方程采用双曲蠕变规律对材料的率相关性与温度相关性进行预测，如下式：Q=C:-p-Ia（TT0）i3s:s28p=3Ae（-詈）p2sign（1-2S*

4、QeRT2.l8p:8pS*=53PPA式中，q为Cauchy应力，s为偏应力，C为弹性张量，乞为总应变，为热膨胀系数,gp为非弹性应变速率，A为常数，Q为激活能，m为应变敏感指数，g为应力乘子，R为气体常数，T为绝对温度，T0为参考温度，h0为形变硬化-软化常数，a为与硬化-软化相关的应变敏感系数，S*为变量饱和值，S为系数，n为指数。粘塑性Anand本构方程中，共有9个材料参数:A,Q,g,m,n,h0，S,a以及初始形变阻抗S0。为真实模拟钎焊材料内部损伤变化，引入损伤，演化率如下式：DD=s-sp-Q(T-70):C:E_2p_(T_T0)加入损伤的Anand模型方程如下：Q=(1-D

5、)-C:S_Sp-血(T-T0)D=s_sp_Ia(T_T0):C:s&p1-少(T-T0)*S*,1She10*0SA=SI3s:ssinh(g2)Sa-(iS)sign(1_*)S_Q_eRT3s:s2(1)(2)(3)(4)(5)3.算法与计算流程计算(2)，(3)，(4)式，主要有三种数值算法，向前显式Euler方法，向后隐式Euler方法，中点法。向前显式Euler方法是条件稳定，具有一阶精度；向后隐式Euler方法和中点法是绝对稳定，分别具有一阶和二阶精度，它们均为隐式方法，需要利用迭代解隐式方程。迭代主要有四种方法：普通迭代，牛顿法，弦位法，抛物线法。本报告中主要采用数值绝对稳定

6、的向后Euler方法和中点法两种数值算法，迭代采用普通迭代和弦位法，进行试算比较方法的优劣。3.1向后隐式Euler方法+普通迭代向后隐式Euler计算公式为：y(x0)=y0yn+1=yn+hf(xn+Tyn+1)向后Euler方法是隐式方法，计算yn+1时要解隐式方程，通常用到迭代法。例如，先用向前显式Euler方法的计算结果作为初值，再作迭代，计算格式为：(0)n+1=yn+hf(xn，yn)(k+1)n+1yn+hf(xn+1,yn(k+)1普通迭代的格式为：xk+1r(xk)，判别迭代过程收敛的条件为：xnj)-xn%“采用上述算法，(1)-(5)式数值计算的离散格式可以表述如下：Q

7、(n+1)=(1-D(n+1)C:&(n+1)-(p+1)D(n+1)=D(n)+A-&(n+1)-8P+1)(-)8(n+1)丸屮+3AtAeRT(n+1)P2-Ia(T(n+1)-70)-Ia(T(n+1)-70):C:(n+1)-P+1)-Ia(T丄m,3s(n+1):s(n+1)Vo*sinhG2()S(n+1)s(n+1)(n+1).2s(n+1)：s(n+1)S(n+1)=S(n)+h0S(n+1)1-S*(n+1)asign(1S(n+1)严)2(8(n+1)-即)：(8(n+1)-8聘)*(n+1)人S=SA+Dy即皿+1)-8即)eRTAtA6)(7)(8)(9)10)根据上

8、述计算格式，UMAT子程序的计算流程为：读取由ABAQUS传递给UMAT子程序的q(n),E(n),A_,D(n),&(n),S(n),T(n+1)p和AT，作为计算初值；采用迭代法，联立方程(6)-(10)式，求解Q(n+1),D(n+1)，首+1)及S(n+1)；更新应力及全部状态变量，更新Jacobian矩阵。其中，迭代法的计算流程具体如下：(1)迭代循环开始，针对q(n+1),D(n+1),&器+1)及S(n+1)赋计算初值;(2)将方程(6)-(9)式写成形如y二f(x,y)的迭代格式，由第k步的k+1kkQ(n+1),D(n+1),8卩+1)及S(n+1)计算第k+1步的Q(n+1

9、),D(n+1)，&(n+1)及S(n+1)；(3)分析比较第k步与第k+1步的Q(n+1),D(n+1),8(+1)及s(n+1)，若它们之间的差满足精度要求，结束循环；否则，继续循环。若循环次数大于预定最大循环次数时，迭代失败。向后Euler方法具有绝对数值稳定性，误差具有一阶精度。虽然是绝对稳定的，但是迭代步长仍要受到一定限制。3.2中点法+普通迭代中点法计算公式为：y(x0)=y0hyn+1二yn+2(f(Xn，yn)+f(Xn+1,yn+1)中点法也是隐式方法计算yn+1时要用到迭代法解隐式方程。先用向前显式Euler方法的计算结果作为初值，再作迭代，同样采用普通迭代方法，计算格式为

10、y(0)1=y+hf(x,y)n+1nnn(k+1),h(),r(k)yn+1=yn+2(f(Xn，yn)+f(xn+1，yn+1)采用上述算法，(1)-(5)式数值计算的离散格式可以表述如下：d(n+1)=(1-d(n+1)C:a(n+1)-ap+1)-Ia(T(n+1)-T0)D(n+1)=D(n)+Al(f(n)+f(n+1)f(n)=a(n)-a护)-Ia(T(n)a(n+1)=a(n)Pp弋(g(n)+g(n+1)3AeRT(2sinh(g(11)(12)-T0):C:a(n)-a(，)-Ia(T(n)-T0)(n),2s(n):s(n)(13)(14)(15)(16)S(n+1)=

11、S(n)弋(h(n)+h(n+1)h(n)=Jh0S(n)1一S*(n)asign(1一S(n)R)*3隹(n+1)-a(n):(a(n+1)-a(n)(17)2(a(n+1)-a(n)：(a(n+1)-a(n)Qs*(n)打(18)!3IPP八、PP丿RT(n)AtAe根据上述计算格式，UMAT子程序的计算流程为：读取由ABAQUS传递给UMAT子程序的(n),卫(n),Aa,D(n)，a(n),S(n),T(n+1)和人卩，作为计算初值，并计算f(n)g(n)h(n)。采用迭代法，联立方程(11)-(18)式，求解(n+1),D(n+1)，aP+D及S(n+D；更新应力及全部状态变量，更新

12、Jacobian矩阵。其中，迭代法的计算流程具体如下：迭代循环开始，针对(n+1),D(n+1)，a(n+1)及S(n+1)赋计算初值；将方程(11)(12)(14)(16)式写成形如y=f(x,y)的迭代格式，由第k步的k+1kk- - -丁(n+1),D(n+1),&(n+1)及s(n+D计算第k+1步的工(n+1)D(n+1),&(n+1)及S(n+1)；(3)分析比较第k步与第k+1步的丁(n+1)D(n+1),(n+1)及S(n+1)，若它们之间的差满足精度要求，结束循环；否则，继续循环。若循环次数大于预定最大循环次数时，迭代失败。中点法也具有绝对数值稳定性。它的迭代收敛条件比用向后

13、Euler方法的迭代步长可以大一倍，误差具有二阶精度，比向后Euler高一阶。但中点法每积分一步需要计算两次函数值,精度的提高时以增加计算量为代价的。3.3中点法+弦位法迭代数值计算同样采用中点法，同3.2中的离散格式(11)-(18)。弦位法迭代格式为：/(xk)xk+1=xkf(xk)-f(xk_1)j)根据弦位法迭代格式，UMAT子程序的迭代格式表示为:k(o丿_n(_)珥+1_珥_k()k(jk_珥-1)kk_1k()Dk+1_Dk_丁败_Dk_1丿Sk+1_Skk(%)总气)_k(k_1丿kSJ(19)(20)(21)(22)UMAT子程序中迭代的计算流程为:(1)迭代循环开始，针对

14、丄(n+1),D(n+1)，&P+1)及S(n+1)赋计算初值；将方程(11)(12)(14)(16)式写成形如式(19)-(22)的迭代格式，由第k步的工(n+1)D(n+1)，&(n+1)及S(+1)计算卩(珥)，再利用式(20)-(22)计算第k+1步的丁(n+1)D(n+1)，&(n+1)及S(+1)；分析第k+1步的p(1)，若它的值满足精度要求，结束循环；否则，继续循环。若循环次数大于预定最大循环次数时，迭代失败。弦位法比普通迭代收敛得快，但是计算工作量要增多，需要计算前两步的值。具体的优劣要针对不同的模型来定。3.4应力、应变增量的Jacobian矩阵- - #-采用隐式数值算法

15、，在UMAT子程序中，需要更新应力、应变增量的Jacobian矩阵。推导后得到一致的应力、应变增量的Jacobian矩阵，如下式：-(1-D)C:念1000001000J=(1-D00100000100000100000(1-D)C-3At2b(n+1)eq00000+3At21eq如3九九000九业3九000九九业3000000卩00000000)00卩00卩000000卩00卩0000000卩九+2卩九九九九+2卩九九九九+2卩0000000004.单元测试：根据上述算法和流程，编写了Anand粘塑性本构模型的UMAT子程序，分别为：KANAND.FOR（向后隐式Euler方法+普通迭代）

16、、KMIDDLE.FOR（中点法+普通迭代）、KMIDDLE-XUE.FOR（中点法+弦位法迭代）。三个程序具有相同的状态变量，共8个，即非弹性应变张量的6个分量、非弹性形变阻抗S和损伤值D。程序中涉及15个相同的材料参数，它们依次为：Youngs模量，Possion比，热膨胀系数，参考温度，初始损伤，损伤阈值以及式（2）-（5）中的参数A,Q,g,m,n,h0,j,a以及初始形变阻抗S0。为了测试程序的正确性，对实际模型进行有限元分析。分析采用三维实体单元，单元类型C3D8：An8-nodelinearbrick。单位取毫米，吨，秒制。材料参数如下：表1：Pb90Sn10焊料的粘塑性Anan

17、d方程的材料参数E/MPaX104Va/K-1X10-5T/K0A/s-1,X107QJ/molmh0MPaASMPanX10-3aS0MPaD0D阈值1.620.382.782933.251558370.143178772.734.383.7315.09暂取0.08暂取1三个程序中，程序KMIDDLE-XUE.FOR（中点法+弦位法迭代）暂有收敛问题未解决,故本报告只进行前两程序的测试。分别进行正方体模型和长方体模型测试。4.1正方体测试正方体模型边长取0.02，位移加载2E-007，时间5s。经测试，所编程序在拉伸、剪切和热加载测试中具有统一的正确性，即程序如果可以进行一种载荷，就可以进行

18、其他载荷加- - -载，反之不能加载。故本报告主要利用拉伸载荷来测试程序在不同网格数下的正确性。图1是网格边13个种子的正方体模型单向拉伸时得到的变形情况。图1网格边13个种子的正方体模型(a)向后Euler方法(b)中点法图2应力-应变曲线a曲5.C*4.W3.的TiirneDmax.BACK=8.000126912659419E-002Dmax.MIDDLE=8.000126912659419E-002图3损伤-时间曲线及损伤云图图2(a)是采用向后Euler方法计算单向拉伸时的应力应变曲线，图2(b)是采用中点法得出的曲线，两者都正确反映了粘塑性材料的单拉变形特征。图3是采用两种方法得出

19、的损伤-时间曲线和损伤云图，因完全相同只列出一个，从输出数据读出损伤最大值如图中所示，数值相同。从以上分析可以看出采用两种方法在图1所示网格数下可以计算，并得出完全相同结果。图4是网格边15个种子的正方体模型单向拉伸时得到的变形情况。图4网格边15个种子的正方体模型(b)中点法(a)向后Euler方法02(KW1WIM.WMOTtow图5应力-应变曲线Dmax.BACK=8.000126811290798E-002图6向后Euler方法计算的损伤-时间曲线及损伤云图图5(a)是采用向后Euler方法计算单向拉伸时的应力应变曲线，图2(b)是采用中点法得出的曲线，后者在计算中出现错误停止，图中为

20、计算完成了的曲线。对两者已加载的应变部分比较，得出的曲线相近。图6是采用Euler方法得出的损伤-时间曲线和损伤云图，从输出数据读出损伤最大值如图中所示，与图3中得到的数值相近。为找出中点法在计算中出现错误的原因和出现错误时模型的状态，绘出以下图示。Di砧IXW1- - - -0.000.501.001.502.W2.50100350Timso.no.60tootso2.n2.503.003.5a4.00fl.QQxJMJuc42-a-ts匸上d.xmUJTima(a)E-TIME(b)E:Max.Principal-TIMEs-tfl宀Dfyse5_sLLj0.000.501.W1502.0

21、02.50工轴3.S机侗TimeD.DflmD.ga.mfawclki.ro1jsue砂常mTime(c)SandS.Mises-TIME(d)S:Max.Principal-TIMEUHrdHLbUMHlHLb=DfiMfiGEUfiRIfiLE=DfiMfiGEUfiRIfiLE=DfiMfiGEUfiRIfiLE=DfiMfiGEUfiRIfiLE=DfiMfiGEUfiRIfiLE=整UUUUV4Z33Z1UUZg0D00m4272S710E-002g0D0011290g333SE-92g0D901120i85ifc9E-002g0D0J47SS&1E-92g0D90112713290

22、2iE-02ITERATIONTOSOLUESTRESSFAILSDAMAGEUARIALE=NaNITERATIONTOSQLUESTRESSFAIL.|DAMAGEUARIALE=NaNITERATIONTOSOLUESTRESSFAILSDAMAGEUARIfiLE=NdN(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