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

-1-Anand粘塑性模型的UMAT子程序及验证高军1.引言电子封装及其组件在工艺或者服役过程中,由于功率耗散和环境温度的周期变化,会因为电子印制电路板、芯片和焊点的热膨胀失配，在合金钎焊焊点处产生交变的应力应变,导致焊点的电、热或者机械失效。焊点的热循环失效(可靠性)是电子封装及组装技术中的关键问题之一,受到了人们的普遍关注。焊点体积细小,应力应变很复杂。为了准确模拟焊点在服役条件下的应力应变响应,对可靠性进行评估,必须建立合理有效的描述钎焊合金材料力学响应的本构方程。SnPb基焊锡钎料广泛应用于电子封装领域，作为电的连接和机械的连接。对于钎料的力学性能的试验和本构模型，许多学者都进行了研究。通常SnPb基焊锡钎料具有很强的温度和加载速率的相关性，应该采用统一型粘塑性本构模型描述SnPb钎料的变形行为。在统一型粘塑性本构模型中，应用最广泛的是Anand模型。具有形式简单，模型参数少等特点，在电子焊点的寿命预测中广泛应用。它采用与位错密度、固溶体强化以及晶粒尺寸效应等相关的单一内部变量S描述材料内部状态对塑性流动的宏观阻抗，可以反映粘塑性材料与应变速度、温度相关的变形行为，以及应变率的历史效应、应变硬化和动态回复等特征。目前，很多大型商用有限元软件，如ANSYS、MARC等都把Anand本构模型嵌入到通用材料模型库中供用户使用，但是，ABAQUS的通用材料模型库中缺少Anand模型。因此，本报告目的在于通过ABAQUS的用户子程序接口UMAT，选择合适的算法，将Anand粘塑性本构模型引入ABAQUS中，以便后续的研究。2.Anand本构方程统一型粘塑性Anand本构模型有两个基本特征：(1)在应力空间没有明确的屈服面,故在变形过程中不需要加载/卸载准则,塑性变形在所有非零应力条件下产生。(2)采用单一内部变量描述材料内部状态对塑性流动的宏观阻抗。内部变量(或称变形阻抗)用S标记,具有应力量纲。粘塑性Anand模型的流动方程采用双曲蠕变规律对材料的率相关性与温度相关性进行预测，如下式：nTRQeAppSSppSSsignSShSsssmSssTRQAepTTIpC:32*:32)*1(*10:231):23sinh()(23)]0([:-2-式中，为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。为真实模拟钎焊材料内部损伤变化，引入损伤，演化率如下式：)]0([::)]0([TTIpCTTIpD加入损伤的Anand模型方程如下：)5(:32*)4(:32)*1(*10)3(:231):23sinh()(23)2()]0([::)]0([)1()]0([:)1(nTRQeAppSSppSSsignSShSsssmSssTRQAepTTIpCTTIpDTTIpCD3.算法与计算流程计算（2），（3），（4）式，主要有三种数值算法，向前显式Euler方法，向后隐式Euler方法，中点法。向前显式Euler方法是条件稳定，具有一阶精度；向后隐式Euler方法和中点法是绝对稳定，分别具有一阶和二阶精度，它们均为隐式方法，需要利用迭代解隐式方程。迭代主要有四种方法：普通迭代，牛顿法，弦位法，抛物线法。本报告中主要采用数值绝对稳定的向后Euler方法和中点法两种数值算法，迭代采用普通迭代和弦位法，进行试算比较方法的优劣。3.1向后隐式Euler方法+普通迭代向后隐式Euler计算公式为：)1,1(10)0(nynxhfnynyyxy向后Euler方法是隐式方法，计算1ny时要解隐式方程，通常用到迭代法。例如，先用向-3-前显式Euler方法的计算结果作为初值，再作迭代，计算格式为：))(1,1()1(1),()0(1knynxhfnyknynynxhfnyny普通迭代的格式为：)(1kxkx，判别迭代过程收敛的条件为：)(1)1(1knxknx采用上述算法，（1）-（5）式数值计算的离散格式可以表述如下：）（）（）（）（）（10)1())()1((:))()1((32)1(*9))()1((:))()1((32))1(*)1(1()1(*)1(10)()1(8)1(:)1(23)1(1))1()1(:)1(23sinh())1((23)()1(7)]0)1(()1()1([::)]0)1(()1()1([)()1(6)]0)1(()1()1([:))1(1()1(nnTRQetAnpnpnpnpSnSnpnpnpnpnSnSsignnSnShnSnSnsnsnsmnSnsnsnTRQAetnpnpTnTInpnCTnTInpntnDnDTnTInpnCnDn根据上述计算格式，UMAT子程序的计算流程为：（1）读取由ABAQUS传递给UMAT子程序的)1(,)(,)()(,,)(,)(nTnSnpnDnn，和T，作为计算初值；（2）采用迭代法，联立方程（6）-（10）式，求解)1()1()1(,)1(nSnpnDn及，；（3）更新应力及全部状态变量，更新Jacobian矩阵。其中，迭代法的计算流程具体如下：（1）迭代循环开始，针对)1()1()1(,)1(nSnpnDn及，赋计算初值；（2）将方程（6）-（9）式写成形如),(1kkkyxfy的迭代格式，由第k步的)1()1(,)1(npnDn，及)1(nS计算第k+1步的)1()1(,)1(npnDn，及)1(nS；（3）分析比较第k步与第k+1步的)1()1(,)1(npnDn，及)1(nS，若它们之间的差满足精度要求，结束循环；否则，继续循环。若循环次数大于预定最大循环次数时，迭代失败。向后Euler方法具有绝对数值稳定性，误差具有一阶精度。虽然是绝对稳定的，但是迭代步长仍要受到一定限制。-4-3.2中点法+普通迭代中点法计算公式为：))1,1(),((210)0(nynxfnynxfhnynyyxy中点法也是隐式方法，计算1ny时要用到迭代法解隐式方程。先用向前显式Euler方法的计算结果作为初值，再作迭代，同样采用普通迭代方法，计算格式为：))(1,1(),((2)1(1),()0(1knynxfnynxfhnyknynynxhfnyny采用上述算法，（1）-（5）式数值计算的离散格式可以表述如下：)18()())()1((:))()1((32)(*)17())()1((:))()1((32))(*)(1()(*)(10)()16())1()((2)()1()15()(:)(23)(1))()(:)(23sinh())((23)()14())1()((2)()1()13()]0)(()()([::)]0)(()()([)()12())1()((2)()1()11()]0)1(()1()1([:))1(1()1(nnTRQetAnpnpnpnpSnSnpnpnpnpnSnSsignnSnShnhnhnhtnSnSnsnsnsmnSnsnsnTRQAengngngtnpnpTnTInpnCTnTInpnnfnfnftnDnDTnTInpnCnDn根据上述计算格式，UMAT子程序的计算流程为：（1）读取由ABAQUS传递给UMAT子程序的)1(,)(,)()(,,)(,)(nTnSnpnDnn，和T，作为计算初值，并计算)(,)(,)(nhngnf。（2）采用迭代法，联立方程（11）-（18）式，求解)1()1()1(,)1(nSnpnDn及，；（3）更新应力及全部状态变量，更新Jacobian矩阵。其中，迭代法的计算流程具体如下：（1）迭代循环开始，针对)1()1()1(,)1(nSnpnDn及，赋计算初值；（2）将方程（11）（12）（14）（16）式写成形如),(1kkkyxfy的迭代格式，由第k步的-5-)1()1(,)1(npnDn，及)1(nS计算第k+1步的)1()1(,)1(npnDn，及)1(nS；（3）分析比较第k步与第k+1步的)1()1(,)1(npnDn，及)1(nS，若它们之间的差满足精度要求，结束循环；否则，继续循环。若循环次数大于预定最大循环次数时，迭代失败。中点法也具有绝对数值稳定性。它的迭代收敛条件比用向后Euler方法的迭代步长可以大一倍，误差具有二阶精度，比向后Euler高一阶。但中点法每积分一步需要计算两次函数值，精度的提高时以增加计算量为代价的。3.3中点法+弦位法迭代数值计算同样采用中点法，同3.2中的离散格式（11）-（18）。弦位法迭代格式为：)1()1()()(1kxkxkxfkxfkxfkxkx根据弦位法迭代格式，UMAT子程序的迭代格式表示为：)22(111)21(111)20(111)19()kSk(S)kk(σ)kk(σ)nk(σkSkS)kDk(D)kk(σ)kk(σ)nk(σkDkD)kσk(σ)kk(σ)kk(σ)nk(σkσkσkσtrykσ)kp(σUMAT子程序中迭代的计算流程为：（1）迭代循环开始，针对)1()1()1(,)1(nSnpnDn及，赋计算初值；（2）将方程（11）（12）（14）（16）式写成形如式（19）-（22）的迭代格式，由第k步的)1()1(,)1(npnDn，及)1(nS计算)kp(σ，再利用式（20）-（22）计算第k+1步的)1()1(,)1(npnDn，及)1(nS；（3）分析第k+1步的)kp(σ1，若它的值满足精度要求，结束循环；否则，继续循环。若循环次数大于预定最大循环次数时，迭代失败。弦位法比普通迭代收敛得快，但是计算工作量要增多，需要计算前两步的值。具体的优劣要针对不同的模型来定。3.4应力、应变增量的Jacobian矩阵-6-采用隐式数值算法，在UMAT子程序中，需要更新应力、应变增量的Jacobian矩阵。推导后得到一致的应力、应变增量的Jacobian矩阵，如下式：0000000000000000002000200021)000000000000000000340003400034)1()1(231000000100000010000001000