周建雄， 魏志剛，毛 歡， 韋雅寧， 劉迎松

(安徽工業(yè)大學(xué) 機(jī)械工程學(xué)院，馬鞍山 243032)

1 引 言

橡膠材料由于其獨(dú)特的物理和化學(xué)性質(zhì)，作為減震和密封等部件廣泛應(yīng)用于工業(yè)工程中，建立真實(shí)可靠的本構(gòu)關(guān)系來(lái)描述橡膠材料的力學(xué)特性變得尤為重要[1]。部分模型已經(jīng)寫(xiě)入有限元材料庫(kù)中用于仿真設(shè)計(jì)，如Arruda-Boyce[2]、Mooney-Rivlin[3]、Yeoh[4]和Ogden[5]等模型在A(yíng)BAQUS中只需給出材料參數(shù)就可用于仿真模擬，而本構(gòu)模型的選取也會(huì)影響仿真結(jié)果的準(zhǔn)確性。近些年涌現(xiàn)較多擬合性和可靠性較好的本構(gòu)模型，尚未開(kāi)發(fā)相應(yīng)模型本構(gòu)子程序用于有限元做仿真模擬以更好地滿(mǎn)足工程應(yīng)用和科學(xué)研究[6]。

橡膠材料的彈性本構(gòu)方程可表示為應(yīng)變不變量的函數(shù)，如Mooney-Rivlin模型和Yeoh模型，也可表示為主伸長(zhǎng)率 的函數(shù)形式，如Ogden模型和Valanis-Landel[7]模型。在A(yíng)BAQUS中，不變量表示的超彈性模型，可用UHYPER編寫(xiě)材料子程序，F(xiàn)ranus等[8]對(duì)MCIZ模型編寫(xiě)UHYPER子程序用于有限元實(shí)現(xiàn)對(duì)微可壓縮超彈性薄殼的膨脹分析，Jemiolo等[9]編寫(xiě)Mooney-Rivlin和Blatz-Ko模型的UHYPER子程序，研究材料模型泊松比的取值對(duì)殼體結(jié)構(gòu)的穩(wěn)定性影響。由伸長(zhǎng)率表示的超彈性本構(gòu)模型需要用更具一般性的UMAT子程序編寫(xiě)。Asim等[10]開(kāi)發(fā)了用于表征彈性材料線(xiàn)粘彈性和非線(xiàn)性超彈性的UMAT子程序，Suchocki等[11]開(kāi)發(fā)了冪指數(shù)模型不可壓縮超彈性UMAT子程序。

本文開(kāi)發(fā)了幾個(gè)新提出的橡膠彈模型的UHYPER和UMAT子程序?；陂_(kāi)發(fā)的UHYPER子程序進(jìn)行了多孔橡膠板單軸拉伸有限元仿真，驗(yàn)證了子程序的正確性以及本構(gòu)模型對(duì)復(fù)雜應(yīng)變場(chǎng)預(yù)測(cè)的準(zhǔn)確性；基于開(kāi)發(fā)的UMAT子程序進(jìn)行單軸、雙軸和剪切變形實(shí)驗(yàn)?zāi)M，驗(yàn)證了子程序的可靠性。

2 本構(gòu)關(guān)系的數(shù)值計(jì)算

2.1 相關(guān)物理量的計(jì)算方法

在UMAT中編寫(xiě)主伸長(zhǎng)率形式的彈性模型，需要定義應(yīng)力張量和彈性張量，下面簡(jiǎn)要介紹這些量的得出方法。

2.1.1 變形張量

在連續(xù)介質(zhì)力學(xué)變形理論中[12]，對(duì)于任何detF≠0的張量F，可極分解為F=RU=VR。式中R為正交張量，代表純轉(zhuǎn)動(dòng)；右和左伸長(zhǎng)張量U和V為對(duì)稱(chēng)正定張量，代表純粹的變形。在有限應(yīng)變狀況下，常用右和左Cauchy-Green變形張量C和B來(lái)描述物體的變形。

C=FTF=(URT)(RU)=U2

B=FFT=(VR)(RTV)=V2

(1)

體積變化率J=V/V0=|F|，通過(guò)求解變形梯度F的行列式獲得。

(2)

在數(shù)值計(jì)算中，所需特征值和特征向量可用迭代雅可比方法顯式計(jì)算。Kopp[13]利用Jacobi算法、QL算法和Cuppen算法完成了三階矩陣對(duì)角化的數(shù)值計(jì)算，并提供了C語(yǔ)言和Fortran語(yǔ)言開(kāi)源代碼，其中Jacobi算法最為精確，適合求解本文所需的特征值和特征向量，可將該代碼寫(xiě)成子程序編入U(xiǎn)MAT中。

2.1.2 應(yīng)力張量和彈性張量

Connolly等[14]給出了主伸長(zhǎng)率表示的超彈性本構(gòu)模型數(shù)值實(shí)現(xiàn)方法，并基于Fortran77編寫(xiě)UMAT子程序用于A(yíng)BAQUS，但其存在一定的局限性，其編寫(xiě)的UMAT子程序中未加入體積應(yīng)變能部分。本文給出了應(yīng)力張量和彈性張量詳細(xì)推導(dǎo)過(guò)程，并將體積應(yīng)變能函數(shù)代入公式推導(dǎo)，基于Fortran77編寫(xiě)UMAT子程序用于有限元模擬。

(1) 應(yīng)力張量

材料的變形應(yīng)變能可分解為等容變形和體積變形兩部分，即

(3)

對(duì)應(yīng)變能函數(shù)求導(dǎo)可得第二類(lèi)Piola-Kirchhoff應(yīng)力

(4)

式中 利用鏈?zhǔn)角髮?dǎo)得

(a=1,2,3) (5)

橡膠材料變形一般為大變形，有限元仿真中使用Cauchy應(yīng)力σ，其可由第二類(lèi)Piola-Kirchhoff應(yīng)力S獲得

(6)

第一類(lèi)Piola-Kirchhoff主應(yīng)力Pa(名義應(yīng)力)在工程中普遍應(yīng)用，Cauchy應(yīng)力與其變換關(guān)系為

(7)

可根據(jù)本構(gòu)模型的函數(shù)形式選擇應(yīng)力變換，本文用于編寫(xiě)UMAT的Wei模型是用名義應(yīng)力P表示的柯西彈性模型，則需將P應(yīng)力變換為σ應(yīng)力。

(2) 彈性張量

由于超彈性具有率無(wú)關(guān)性，物質(zhì)彈性張量可以通過(guò)第二類(lèi)Piola-Kirchhoff應(yīng)力和變形張量等價(jià)定義為

Q=2?S/?C

(8)

物質(zhì)彈性張量分解為等容部分和體積變化部分,Q=Qiso+Qvol。

等容彈性張量的譜形式為[15]

Na?Nb?Nb?Na)

(9)

式(9)第二項(xiàng)，當(dāng)λa和λb數(shù)值相近時(shí)，分子分母趨于零，由洛必達(dá)法則得

(10)

(11)

體積變化部分為

(12)

上述本構(gòu)模型中，材料是假定為不可壓縮的，體積應(yīng)變能為零。但完全不可壓縮會(huì)在有限元仿真變形中導(dǎo)致數(shù)值奇異，因此有限元仿真實(shí)際上認(rèn)為材料是近似不可壓縮的。Abaqus中采用的體積應(yīng)變能表達(dá)式為

(13)

式中Da取值決定了材料的可壓縮性，Da取零則認(rèn)為材料是完全不可壓縮的。初始體積模量為k0=1/D1，體積模量過(guò)大會(huì)造成模型體積自鎖，根據(jù)有限元模擬結(jié)果對(duì)Da取值進(jìn)行不斷調(diào)整，Da≈10-3時(shí),有限元結(jié)果與模型理論值一致，且求解較為迅速。

在A(yíng)BAQUS中，空間彈性張量用Cauchy應(yīng)力的Jaumann率來(lái)定義[16]，表示為QJ。

QJ=J-1FFQFTFT=

J-1(Qiso+Qvol)+[(σ⊙I)+(I⊙σ)]=

[(na?nb)?(na?nb+nb?na)]+

[(na?nb)?(na?nb+nb?na)]+

[(σ⊙I)+(I⊙σ)]

(14)

最終的表達(dá)式中要求的量很多，且都與第二類(lèi)Piola-Kirchhoff應(yīng)力S有關(guān)，在應(yīng)力張量部分有提到，Cauchy應(yīng)力σ、第二類(lèi)Piola-Kirchhoff應(yīng)力S和第一類(lèi)Piola-Kirchhoff應(yīng)力P(名義應(yīng)力)存在變換關(guān)系，則不管本構(gòu)模型由何種應(yīng)力描述，都可變換為第二Piola-Kirchhoff應(yīng)力S進(jìn)行推導(dǎo)。本文基于Fortran77完成對(duì)Wei模型和Ogden模型UMAT子程序開(kāi)發(fā)，并將其用于有限元模擬來(lái)驗(yàn)證子程序的正確性。

2.2 使用模型

2.2.1 伸長(zhǎng)率形式模型

(1) Wei模型[17]

在不考慮材料不可壓縮特性時(shí)，Wei模型名義主應(yīng)力為

(15)

式中i,j,k=1,2,3,且取值互不相等，a,c,d,e和m為材料常數(shù)，擬合橡膠三類(lèi)簡(jiǎn)單試驗(yàn)數(shù)據(jù)[18]得,a=0.001，c=2.04，d=0.42，e=0.17,m=0.85。

對(duì)于不可壓縮材料，該模型的柯西主應(yīng)力偏量表達(dá)式為

(16)

將柯西主應(yīng)力轉(zhuǎn)變?yōu)榫幊趟杳x應(yīng)力偏量

(17)

名義主應(yīng)力偏量關(guān)于修正伸長(zhǎng)率的偏導(dǎo)為

(18)

(2) Ogden模型

(19)

式中μi和αi為材料參數(shù)，本文選擇擬合性較好的三階Ogden模型，擬合橡膠三類(lèi)簡(jiǎn)單試驗(yàn)數(shù)據(jù)得[19],μ1=0.006，α1=4.17，μ2=0.709，α2=1.108，μ3=-0.006，α3=-2.125。

名義主應(yīng)力偏量可由等容應(yīng)變能對(duì)修正伸長(zhǎng)率求偏導(dǎo)得

(20)

對(duì)于所需的名義主應(yīng)力偏量及其關(guān)于修正伸長(zhǎng)率偏導(dǎo)可參考式(17,18)求解。

2.2.2 不變量形式模型

使用UHYPER開(kāi)發(fā)子程序遵循固定格式，需定義應(yīng)變能函數(shù)，并計(jì)算應(yīng)變能關(guān)于不變量和體積變化率的一階、二階和三階偏導(dǎo)數(shù)，將其存儲(chǔ)到相應(yīng)數(shù)組中，常用變量有

U(1)=W

(21)

選取三個(gè)近來(lái)提出的超彈性模型編寫(xiě)UHYPER子程序來(lái)探究其預(yù)測(cè)復(fù)雜應(yīng)變場(chǎng)的能力。模型如下。

(1) Exp-Ln[20]模型

(22)

式中A，a和b為材料常數(shù)，擬合未填充多孔硅橡膠板拉伸試驗(yàn)數(shù)據(jù)[21]得A=0.166，a=0.202，b=0.386。

(2) Mansouri[22]模型

W=A1{exp[m1(I1-3)]-1}+

B1{exp[n1(I1-3)]-1}

(23)

式中A1,B1,m1和n1為材料常數(shù)，擬合未填充多孔硅橡膠板拉伸試驗(yàn)數(shù)據(jù)得A1=1.54，B1=5.57，m1=0.089，n1=0.004。

(3) New neo -Hookean[23]模型

(24)

3 數(shù)值模擬

3.1 非均勻拉伸變形有限元模擬

Meunier對(duì)未填充多孔硅橡膠板進(jìn)行拉伸試驗(yàn)，描述橡膠材料在復(fù)雜應(yīng)變狀態(tài)下的力學(xué)性質(zhì)，且拉伸試驗(yàn)數(shù)據(jù)可用于評(píng)估超彈性本構(gòu)模型預(yù)測(cè)復(fù)雜應(yīng)變場(chǎng)的準(zhǔn)確性。本文使用UHYPER子程序編寫(xiě)上述本構(gòu)模型的子程序進(jìn)行多孔橡膠板的拉伸仿真，評(píng)估這些模型預(yù)測(cè)結(jié)果的精度。

圖1為采用的多孔橡膠板的形狀尺寸，孔中心C3和C5之間沿圓心線(xiàn)做了切割。仿真中將底端固定，限制橫縱方向上的自由度，在上頂面沿豎直方向施加57.3 mm的位移載荷，限制水平自由度。

圖1 多孔橡膠板幾何尺寸

有限元仿真得出的變形圖輪廓(紅色輪廓)與Meunier試驗(yàn)變形圖(灰色圖片)重合比較如圖2所示。

圖2 試驗(yàn)與仿真變形對(duì)比

在單軸拉伸變形狀態(tài)下，伸長(zhǎng)率為1.7時(shí)，圖2中紅色輪廓線(xiàn)與灰色圖片重合度非常好，表明三個(gè)模型仿真結(jié)果和試驗(yàn)結(jié)果十分契合，所編子程序能夠按照本構(gòu)方程模擬計(jì)算，且所建橡膠本構(gòu)模型可以準(zhǔn)確預(yù)測(cè)橡膠在復(fù)雜應(yīng)變場(chǎng)下的力學(xué)性質(zhì)。由于此拉伸試驗(yàn)伸長(zhǎng)率僅為1.7，是小變形試驗(yàn)，很難分辨出三個(gè)模型預(yù)測(cè)復(fù)雜應(yīng)變場(chǎng)的優(yōu)劣，可進(jìn)一步增加伸長(zhǎng)率，圖3為伸長(zhǎng)率為2.5時(shí)各模型預(yù)測(cè)多孔橡膠板的應(yīng)力分布及變形仿真。

圖3 應(yīng)力分布及變形輪廓對(duì)比

可以看出，Exp-Ln模型和New neo-Hookean模型預(yù)測(cè)多孔橡膠板在伸長(zhǎng)率2.5時(shí)的應(yīng)力分布及大小較為相近，且變形輪廓重合度較好，相比之下，Mansouri模型預(yù)測(cè)的應(yīng)力要小近30%，且孔洞附近的變形較前兩個(gè)模型要稍大些。

3.2 單軸、等雙軸和純剪變形有限元模擬

在A(yíng)BAQUS中使用由Wei模型及Ogden模型編寫(xiě)的UMAT子程序?qū)θ?lèi)簡(jiǎn)單拉伸實(shí)驗(yàn)做有限元模擬，獲得柯西應(yīng)力σ關(guān)于名義應(yīng)變?chǔ)诺姆抡娼Y(jié)果，由P=σ/(1+ε)將其轉(zhuǎn)換為名義應(yīng)力關(guān)于伸長(zhǎng)率的關(guān)系，并將仿真結(jié)果與本構(gòu)模型求解所得主變形方向下名義應(yīng)力解析解進(jìn)行對(duì)比。

3.2.1 Wei模型

(25)

圖4 單軸拉伸狀態(tài)下理論值與應(yīng)用UMAT子程序獲得有限元仿真結(jié)果對(duì)比

對(duì)比上述單軸、等雙軸和純剪變形有限元仿真結(jié)果和理論結(jié)果，直觀(guān)上發(fā)現(xiàn)兩者的曲線(xiàn)較為吻合，進(jìn)一步繪制三類(lèi)簡(jiǎn)單變形下仿真結(jié)果與理論結(jié)果間的相對(duì)誤差，如圖7所示。

圖5 等雙軸拉伸狀態(tài)下理論值與應(yīng)用UMAT子程序獲得有限元仿真結(jié)果對(duì)比

圖6 剪切變形狀態(tài)下理論值與應(yīng)用UMAT子程序獲得的有限元仿真結(jié)果對(duì)比

圖7 有限元模擬與Wei本構(gòu)方程理論結(jié)果相對(duì)誤差

由圖7可知，對(duì)Wei模型所編UMAT用于橡膠單軸與雙軸變形模擬的相對(duì)誤差隨伸長(zhǎng)率逐漸增大，但最大誤差不超過(guò)3%，而UMAT對(duì)剪切變形計(jì)算的相對(duì)誤差隨伸長(zhǎng)率逐漸減小，最大相對(duì)誤差出現(xiàn)在初始變形階段，數(shù)值接近7%。表明所編UMAT與本構(gòu)方程解析解存在偏差，該偏差可通過(guò)在A(yíng)BAQUS中減小 取值來(lái)減小相對(duì)誤差。

3.2.2 Ogden模型

與Wei模型的處理方法相同，結(jié)合式(20,25)得Ogden模型名義主應(yīng)力理論值

(26)

圖8 單軸拉伸狀態(tài)下理論值與應(yīng)用UMAT子程序獲得的有限元仿真結(jié)果對(duì)比

圖9 等雙軸拉伸狀態(tài)下理論值與應(yīng)用UMAT子程序獲得的有限元仿真結(jié)果對(duì)比

圖10 剪切變形狀態(tài)下理論值與應(yīng)用UMAT子程序獲得的有限元仿真結(jié)果對(duì)比

圖11 有限元模擬與Ogden本構(gòu)方程理論結(jié)果相對(duì)誤差

對(duì)Ogden模型所編UMAT用于橡膠三類(lèi)簡(jiǎn)單變形模擬的相對(duì)誤差隨伸長(zhǎng)率逐漸增大，其相對(duì)誤差較小，最大誤差不超過(guò)2.5%。結(jié)合上文由Wei模型編寫(xiě)UMAT對(duì)三類(lèi)簡(jiǎn)單變形實(shí)驗(yàn)?zāi)M結(jié)果可知，所編子程序能夠按照本構(gòu)模型給出的力學(xué)表達(dá)式正確模擬計(jì)算，其相對(duì)誤差來(lái)源于體積應(yīng)變能函數(shù)Da的取值無(wú)法趨近于0，在模擬計(jì)算中無(wú)法避免。

4 結(jié) 論

針對(duì)各向同性超彈性本構(gòu)模型數(shù)值計(jì)算，本文主要推導(dǎo)了主伸長(zhǎng)率表示的超彈性本構(gòu)模型數(shù)值實(shí)現(xiàn)方法，并基于Fortran77編寫(xiě)UMAT用于有限元模擬，所得結(jié)果如下。

(1) 結(jié)合Exp-Ln模型、Mansouri模型和New neo-Hookean模型開(kāi)發(fā)UHYPER子程序用于多孔橡膠板拉伸仿真，模擬結(jié)果與實(shí)驗(yàn)數(shù)據(jù)十分契合，其中Exp-Ln模型和New neo-Hookean模型能更準(zhǔn)確地預(yù)測(cè)較高伸長(zhǎng)率下的復(fù)雜應(yīng)變場(chǎng)。

(2) 由Wei模型和Ogden模型開(kāi)發(fā)UMAT子程序用于單軸、雙軸和剪切變形有限元仿真，仿真結(jié)果與本構(gòu)方程解析解十分接近，表明所編子程序能夠按照本構(gòu)模型給出的力學(xué)表達(dá)式正確模擬計(jì)算，其中的偏差來(lái)源于A(yíng)BAQUS中為避免出現(xiàn)體積自鎖，體積應(yīng)變能函數(shù)中Di的取值不能完全趨于0。

本文所提方法可用于主伸長(zhǎng)率表示的超彈性本構(gòu)模型數(shù)值計(jì)算，并編寫(xiě)UMAT子程序，可供后續(xù)用戶(hù)開(kāi)發(fā)使用。改進(jìn)之處在于可編寫(xiě)相關(guān)求導(dǎo)子程序來(lái)替代用戶(hù)求解主應(yīng)力對(duì)伸長(zhǎng)率的偏導(dǎo)。