Anand粘塑性模型的UMAT子程序及驗(yàn)證

1、- - -Anand粘塑性模型的UMAT子程序及驗(yàn)證高軍1.引言電子封裝及其組件在工藝或者服役過程中,由于功率耗散和環(huán)境溫度的周期變化,會因?yàn)殡娮佑≈齐娐钒?、芯片和焊點(diǎn)的熱膨脹失配，在合金釬焊焊點(diǎn)處產(chǎn)生交變的應(yīng)力應(yīng)變,導(dǎo)致焊點(diǎn)的電、熱或者機(jī)械失效。焊點(diǎn)的熱循環(huán)失效（可靠性）是電子封裝及組裝技術(shù)中的關(guān)鍵問題之一,受到了人們的普遍關(guān)注。焊點(diǎn)體積細(xì)小,應(yīng)力應(yīng)變很復(fù)雜。為了準(zhǔn)確模擬焊點(diǎn)在服役條件下的應(yīng)力應(yīng)變響應(yīng),對可靠性進(jìn)行評估,必須建立合理有效的描述釬焊合金材料力學(xué)響應(yīng)的本構(gòu)方程。SnPb基焊錫釬料廣泛應(yīng)用于電子封裝領(lǐng)域，作為電的連接和機(jī)械的連接。對于釬料的力學(xué)性能的試驗(yàn)和本構(gòu)模型，許多學(xué)者都進(jìn)行了

2、研究。通常SnPb基焊錫釬料具有很強(qiáng)的溫度和加載速率的相關(guān)性，應(yīng)該采用統(tǒng)一型粘塑性本構(gòu)模型描述SnPb釬料的變形行為。在統(tǒng)一型粘塑性本構(gòu)模型中，應(yīng)用最廣泛的是Anand模型。具有形式簡單，模型參數(shù)少等特點(diǎn)，在電子焊點(diǎn)的壽命預(yù)測中廣泛應(yīng)用。它采用與位錯密度、固溶體強(qiáng)化以及晶粒尺寸效應(yīng)等相關(guān)的單一內(nèi)部變量S描述材料內(nèi)部狀態(tài)對塑性流動的宏觀阻抗，可以反映粘塑性材料與應(yīng)變速度、溫度相關(guān)的變形行為，以及應(yīng)變率的歷史效應(yīng)、應(yīng)變硬化和動態(tài)回復(fù)等特征。目前，很多大型商用有限元軟件，如ANSYS、MARC等都把Anand本構(gòu)模型嵌入到通用材料模型庫中供用戶使用，但是，ABAQUS的通用材料模型庫中缺少Anand

3、模型。因此，本報告目的在于通過ABAQUS的用戶子程序接口UMAT,選擇合適的算法，將Anand粘塑性本構(gòu)模型引入ABAQUS中，以便后續(xù)的研究。2Anand本構(gòu)方程統(tǒng)一型粘塑性Anand本構(gòu)模型有兩個基本特征：（1）在應(yīng)力空間沒有明確的屈服面，故在變形過程中不需要加載/卸載準(zhǔn)則,塑性變形在所有非零應(yīng)力條件下產(chǎn)生。（2）采用單一內(nèi)部變量描述材料內(nèi)部狀態(tài)對塑性流動的宏觀阻抗。內(nèi)部變量（或稱變形阻抗）用S標(biāo)記，具有應(yīng)力量綱。粘塑性Anand模型的流動方程采用雙曲蠕變規(guī)律對材料的率相關(guān)性與溫度相關(guān)性進(jìn)行預(yù)測，如下式：Q=C:-p-Ia（TT0）i3s:s28p=3Ae（-詈）p2sign（1-2S*

4、QeRT2.l8p:8pS*=53PPA式中，q為Cauchy應(yīng)力，s為偏應(yīng)力，C為彈性張量，乞?yàn)榭倯?yīng)變，為熱膨脹系數(shù),gp為非彈性應(yīng)變速率，A為常數(shù)，Q為激活能，m為應(yīng)變敏感指數(shù)，g為應(yīng)力乘子，R為氣體常數(shù)，T為絕對溫度，T0為參考溫度，h0為形變硬化-軟化常數(shù)，a為與硬化-軟化相關(guān)的應(yīng)變敏感系數(shù)，S*為變量飽和值，S為系數(shù)，n為指數(shù)。粘塑性Anand本構(gòu)方程中，共有9個材料參數(shù):A,Q,g,m,n,h0，S,a以及初始形變阻抗S0。為真實(shí)模擬釬焊材料內(nèi)部損傷變化，引入損傷，演化率如下式：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)式，主要有三種數(shù)值算法，向前顯式Euler方法，向后隱式Euler方法，中點(diǎn)法。向前顯式Euler方法是條件穩(wěn)定，具有一階精度；向后隱式Euler方法和中點(diǎn)法是絕對穩(wěn)定，分別具有一階和二階精度，它們均為隱式方法，需要利用迭代解隱式方程。迭代主要有四種方法：普通迭代，牛頓法，弦位法，拋物線法。本報告中主要采用數(shù)值絕對穩(wěn)定

6、的向后Euler方法和中點(diǎn)法兩種數(shù)值算法，迭代采用普通迭代和弦位法，進(jìn)行試算比較方法的優(yōu)劣。3.1向后隱式Euler方法+普通迭代向后隱式Euler計算公式為：y(x0)=y0yn+1=yn+hf(xn+Tyn+1)向后Euler方法是隱式方法，計算yn+1時要解隱式方程，通常用到迭代法。例如，先用向前顯式Euler方法的計算結(jié)果作為初值，再作迭代，計算格式為：(0)n+1=yn+hf(xn，yn)(k+1)n+1yn+hf(xn+1,yn(k+)1普通迭代的格式為：xk+1r(xk)，判別迭代過程收斂的條件為：xnj)-xn%“采用上述算法，(1)-(5)式數(shù)值計算的離散格式可以表述如下：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)嚴(yán))2(8(n+1)-即)：(8(n+1)-8聘)*(n+1)人S=SA+Dy即皿+1)-8即)eRTAtA6)(7)(8)(9)10)根據(jù)上

8、述計算格式，UMAT子程序的計算流程為：讀取由ABAQUS傳遞給UMAT子程序的q(n),E(n),A_,D(n),&(n),S(n),T(n+1)p和AT，作為計算初值；采用迭代法，聯(lián)立方程(6)-(10)式，求解Q(n+1),D(n+1)，首+1)及S(n+1)；更新應(yīng)力及全部狀態(tài)變量，更新Jacobian矩陣。其中，迭代法的計算流程具體如下：(1)迭代循環(huán)開始，針對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)，若它們之間的差滿足精度要求，結(jié)束循環(huán)；否則，繼續(xù)循環(huán)。若循環(huán)次數(shù)大于預(yù)定最大循環(huán)次數(shù)時，迭代失敗。向后Euler方法具有絕對數(shù)值穩(wěn)定性，誤差具有一階精度。雖然是絕對穩(wěn)定的，但是迭代步長仍要受到一定限制。3.2中點(diǎn)法+普通迭代中點(diǎn)法計算公式為：y(x0)=y0hyn+1二yn+2(f(Xn，yn)+f(Xn+1,yn+1)中點(diǎn)法也是隱式方法計算yn+1時要用到迭代法解隱式方程。先用向前顯式Euler方法的計算結(jié)果作為初值，再作迭代，同樣采用普通迭代方法，計算格式為

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)式數(shù)值計算的離散格式可以表述如下：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護(hù))-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根據(jù)上述計算格式，UMAT子程序的計算流程為：讀取由ABAQUS傳遞給UMAT子程序的(n),衛(wèi)(n),Aa,D(n)，a(n),S(n),T(n+1)和人卩，作為計算初值，并計算f(n)g(n)h(n)。采用迭代法，聯(lián)立方程(11)-(18)式，求解(n+1),D(n+1)，aP+D及S(n+D；更新應(yīng)力及全部狀態(tài)變量，更新

12、Jacobian矩陣。其中，迭代法的計算流程具體如下：迭代循環(huán)開始，針對(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)，若它們之間的差滿足精度要求，結(jié)束循環(huán)；否則，繼續(xù)循環(huán)。若循環(huán)次數(shù)大于預(yù)定最大循環(huán)次數(shù)時，迭代失敗。中點(diǎn)法也具有絕對數(shù)值穩(wěn)定性。它的迭代收斂條件比用向后

13、Euler方法的迭代步長可以大一倍，誤差具有二階精度，比向后Euler高一階。但中點(diǎn)法每積分一步需要計算兩次函數(shù)值,精度的提高時以增加計算量為代價的。3.3中點(diǎn)法+弦位法迭代數(shù)值計算同樣采用中點(diǎn)法，同3.2中的離散格式(11)-(18)。弦位法迭代格式為：/(xk)xk+1=xkf(xk)-f(xk_1)j)根據(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)迭代循環(huán)開始，針對

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)，若它的值滿足精度要求，結(jié)束循環(huán)；否則，繼續(xù)循環(huán)。若循環(huán)次數(shù)大于預(yù)定最大循環(huán)次數(shù)時，迭代失敗。弦位法比普通迭代收斂得快，但是計算工作量要增多，需要計算前兩步的值。具體的優(yōu)劣要針對不同的模型來定。3.4應(yīng)力、應(yīng)變增量的Jacobian矩陣- - #-采用隱式數(shù)值算法

15、，在UMAT子程序中，需要更新應(yīng)力、應(yīng)變增量的Jacobian矩陣。推導(dǎo)后得到一致的應(yīng)力、應(yīng)變增量的Jacobian矩陣，如下式：-(1-D)C:念1000001000J=(1-D00100000100000100000(1-D)C-3At2b(n+1)eq00000+3At21eq如3九九000九業(yè)3九000九九業(yè)3000000卩00000000)00卩00卩000000卩00卩0000000卩九+2卩九九九九+2卩九九九九+2卩0000000004.單元測試：根據(jù)上述算法和流程，編寫了Anand粘塑性本構(gòu)模型的UMAT子程序，分別為：KANAND.FOR（向后隱式Euler方法+普通迭代）

16、、KMIDDLE.FOR（中點(diǎn)法+普通迭代）、KMIDDLE-XUE.FOR（中點(diǎn)法+弦位法迭代）。三個程序具有相同的狀態(tài)變量，共8個，即非彈性應(yīng)變張量的6個分量、非彈性形變阻抗S和損傷值D。程序中涉及15個相同的材料參數(shù)，它們依次為：Youngs模量，Possion比，熱膨脹系數(shù)，參考溫度，初始損傷，損傷閾值以及式（2）-（5）中的參數(shù)A,Q,g,m,n,h0,j,a以及初始形變阻抗S0。為了測試程序的正確性，對實(shí)際模型進(jìn)行有限元分析。分析采用三維實(shí)體單元，單元類型C3D8：An8-nodelinearbrick。單位取毫米，噸，秒制。材料參數(shù)如下：表1：Pb90Sn10焊料的粘塑性Anan

17、d方程的材料參數(shù)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（中點(diǎn)法+弦位法迭代）暫有收斂問題未解決,故本報告只進(jìn)行前兩程序的測試。分別進(jìn)行正方體模型和長方體模型測試。4.1正方體測試正方體模型邊長取0.02，位移加載2E-007，時間5s。經(jīng)測試，所編程序在拉伸、剪切和熱加載測試中具有統(tǒng)一的正確性，即程序如果可以進(jìn)行一種載荷，就可以進(jìn)行

18、其他載荷加- - -載，反之不能加載。故本報告主要利用拉伸載荷來測試程序在不同網(wǎng)格數(shù)下的正確性。圖1是網(wǎng)格邊13個種子的正方體模型單向拉伸時得到的變形情況。圖1網(wǎng)格邊13個種子的正方體模型(a)向后Euler方法(b)中點(diǎn)法圖2應(yīng)力-應(yīng)變曲線a曲5.C*4.W3.的TiirneDmax.BACK=8.000126912659419E-002Dmax.MIDDLE=8.000126912659419E-002圖3損傷-時間曲線及損傷云圖圖2(a)是采用向后Euler方法計算單向拉伸時的應(yīng)力應(yīng)變曲線，圖2(b)是采用中點(diǎn)法得出的曲線，兩者都正確反映了粘塑性材料的單拉變形特征。圖3是采用兩種方法得出

19、的損傷-時間曲線和損傷云圖，因完全相同只列出一個，從輸出數(shù)據(jù)讀出損傷最大值如圖中所示，數(shù)值相同。從以上分析可以看出采用兩種方法在圖1所示網(wǎng)格數(shù)下可以計算，并得出完全相同結(jié)果。圖4是網(wǎng)格邊15個種子的正方體模型單向拉伸時得到的變形情況。圖4網(wǎng)格邊15個種子的正方體模型(b)中點(diǎn)法(a)向后Euler方法02(KW1WIM.WMOTtow圖5應(yīng)力-應(yīng)變曲線Dmax.BACK=8.000126811290798E-002圖6向后Euler方法計算的損傷-時間曲線及損傷云圖圖5(a)是采用向后Euler方法計算單向拉伸時的應(yīng)力應(yīng)變曲線，圖2(b)是采用中點(diǎn)法得出的曲線，后者在計算中出現(xiàn)錯誤停止，圖中為

20、計算完成了的曲線。對兩者已加載的應(yīng)變部分比較，得出的曲線相近。圖6是采用Euler方法得出的損傷-時間曲線和損傷云圖，從輸出數(shù)據(jù)讀出損傷最大值如圖中所示，與圖3中得到的數(shù)值相近。為找出中點(diǎn)法在計算中出現(xiàn)錯誤的原因和出現(xiàn)錯誤時模型的狀態(tài)，繪出以下圖示。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機(jī)侗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中點(diǎn)法計算結(jié)果從圖7(b)和圖7(d)看出在計算的最后最大主應(yīng)力和最大主應(yīng)變突變?yōu)榱悖瑢?dǎo)致出現(xiàn)圖7(e)的情況，損傷值突變?yōu)闊o窮大。經(jīng)過計算，單元邊種子數(shù)13以下向后歐拉和中點(diǎn)法都可以計算，15以上只有向后歐拉可以繼續(xù)計算，中點(diǎn)法計算中某些變量出現(xiàn)突變，導(dǎo)致無法繼續(xù)計算。具體原因需要進(jìn)一步研究，可能是和算法的具體格式有關(guān)。4.2長方體測試正方體模型邊長取0.02*0.02*0.072,位移加載單向拉伸，時間5s,網(wǎng)格種子9*9*13,采用中點(diǎn)法程序。本部分測試在相同的網(wǎng)格劃分條件下，單向拉伸位移加載的極限。采用位移加載2.88E-5，相應(yīng)長度0.072為0.04%,圖8