王銀龍, 李 釗, 李子然, 汪 洋

(1. 中國科學技術大學 近代力學系,合肥 230027; 2. 中國科學技術大學 中國科學院材料力學行為和設計重點實驗室,合肥 230027; 3. 武漢第二船舶設計研究所, 武漢 430205)

橡膠制品的生產(chǎn)工藝一般包括原料混煉、半成品部件擠出/壓延、成型、硫化定型4個過程.硫化定型工藝中,橡膠內(nèi)部高分子鏈與硫原子之間發(fā)生化學反應,形成硫化交聯(lián)網(wǎng)絡,橡膠的彈性水平提高,黏塑性降低,耐老化性能增強[1],因此硫化橡膠得到較為廣泛的應用,對其力學性能的研究也受到廣大學者的重視.相比于硫化橡膠,未硫化橡膠對產(chǎn)品性能的影響卻常被忽視,目前對其力學行為的研究依然較少.然而,未硫化橡膠的力學性能在很大程度上決定了前期的擠出、壓延及成型工藝效果.事實上,上述過程中工藝參數(shù)如果不能與未硫化橡膠性能相匹配,往往會嚴重降低產(chǎn)品最終使用性能,甚至導致生產(chǎn)缺陷[2-3].近年來,隨著各國對橡膠工業(yè)制造工藝研究的愈加深入,未硫化橡膠的力學性能也受到越來越多的關注.

文獻[4-5]中最早較為系統(tǒng)地測試了未硫化炭黑填充橡膠在常溫環(huán)境中不同變形模式(拉伸、壓縮、剪切)下的非線性力學行為.結果表明,未硫化橡膠的力學響應具有明顯的應變率相關性,以及遲滯效應、Mullins損傷、殘余變形等非線性黏彈塑性特征.根據(jù)實驗結果,文獻[4-5]中提出一個基于Yeoh模型和廣義Maxwell模型的非線性黏彈性本構模型,但此模型無法較好刻畫未硫化橡膠應變率相關的塑性力學行為以及損傷軟化等變形特征.為了改善表征效果,文獻[6-8]中基于微球模型(Micro-sphere Model)提出一個新的內(nèi)變量型黏塑性本構模型,此模型采用微球模型運動學框架對非彈性響應進行描述,無需對變形梯度進行乘法分解,但其表征效果并未得到明顯提升.隨后,上述研究基于黏彈性分數(shù)Zener模型[8]對未硫化橡膠的力學行為進行表征,該模型可以較好反映材料的彈性水平,但對黏塑性的表征效果仍較差.文獻[9-11]中將1個 Marlow 彈簧模型和3個Maxwell模型并聯(lián),獲得了未硫化橡膠Prony級數(shù)形式的黏彈性本構模型.此模型與Kaliske等[4]的本構模型并無本質(zhì)區(qū)別,模型對未硫化橡膠單一工況下的松弛實驗取得了較好的表征效果,尚無法表征未硫化橡膠復雜變形下的非線性黏彈塑性力學行為.需要說明的是,以上研究獲得的本構模型雖然都包含應變率效應,但在本構模型參數(shù)擬合過程中均未同時包含多種應變率工況,本構模型能否較好地描述未硫化橡膠不同應變率下的非線性黏彈塑性響應仍需進一步驗證.Kim等[12]通過多種拉伸實驗探究未硫化橡膠應變率相關的非線性拉伸力學行為,并采用一個基于廣義Maxwell模型和修正Kelvin模型的黏彈-黏塑性本構模型對其非線性力學行為進行表征,但該模型在較高應變率下(0.1 s-1)的表征效果略差,也未能較好反映出未硫化橡膠應變率相關的塑性變形特征.文獻[13-14]中對未硫化橡膠開展不同應變率和不同溫度下的拉伸實驗,并基于八鏈模型和Bergstr?m-Boyce流動法則提出了等溫黏彈塑性本構模型和溫度相關的黏彈塑性本構模型,較好地表征了未硫化橡膠較為復雜的非線性力學行為特征,如應變率效應、遲滯效應、殘余變形等.值得一提的是,上述表征未硫化橡膠力學性能的本構模型[4-14]形式都較為復雜,參數(shù)眾多,只有文獻[6,12]中對提出的基于微球的本構模型給出了具體有限元實現(xiàn)過程,并進行了必要的數(shù)值驗證,而其他本構模型能否用于未硫化橡膠復雜變形條件下的數(shù)值仿真仍有待研究.

本文對未硫化橡膠開展了不同應變率下的單調(diào)拉伸及循環(huán)拉伸加卸載實驗,獲得了未硫化橡膠較為復雜的非線性黏彈塑性力學行為.根據(jù)實驗結果,提出了一個三網(wǎng)絡黏彈塑性本構模型(TN模型)對未硫化橡膠的力學行為進行表征,模型充分考慮了膠料的超彈應力響應、非線性黏塑性流動、Mullins損傷軟化等力學特征,能夠較好表征未硫化橡膠非線性的黏彈塑性拉伸力學行為.最后,依托Abaqus有限元仿真軟件用戶材料子程序接口(UMAT),完成了上述TN模型的有限元實現(xiàn),并對未硫化橡膠多步松弛加卸載和輪胎胎面膠壓入模具的過程進行了有限元仿真,驗證了TN模型的有效性和數(shù)值穩(wěn)定性.

1 實驗及結果討論

1.1 實驗材料及試件準備

實驗采用的膠料均來自安徽佳通輪胎有限公司,為避免試件制備過程中橡膠內(nèi)部發(fā)生硫化反應,膠料均為不含硫化體系的母煉膠.

實驗試件采用啞鈴狀,具體尺寸如圖1(a)所示.未硫化橡膠首先經(jīng)過反復開煉,再置于平板硫化機的高溫高壓環(huán)境下獲得厚度為2 mm的橡膠板,最后通過裁切獲得啞鈴狀試件.實驗過程中采用自動網(wǎng)格法[15-16]記錄試件全場的應變信息,實驗前需在試件表面制作網(wǎng)格點陣,如圖1(b)所示.更為具體的試件制備方法以及網(wǎng)格點陣制備方式已在此前的研究[13]中發(fā)表,在此不再贅述.

1.2 橡膠實驗及結果分析

單調(diào)拉伸工況下未硫化橡膠的拉伸實驗結果如圖2(a)所示.圖中:εeng為工程應變;σeng為工程應力.從圖中可以看出,不同應變率下未硫化橡膠具有相似的非線性力學特征,在拉伸初始時具有很小的線性段,隨后發(fā)生軟化行為,最終應力趨于平緩.膠料拉伸應力水平表現(xiàn)出明顯的應變率相關性,隨著應變率增加,初始剛度顯著增加,應力水平提高.相比于硫化橡膠,未硫化橡膠在100%的應變范圍內(nèi)未出現(xiàn)反“S”現(xiàn)象,可能由其內(nèi)部未形成交聯(lián)網(wǎng)絡所致.

圖2 未硫化橡膠不同應變率下拉伸實驗結果Fig.2 Mechanical behaviors of uncured rubber at different strain rates

循環(huán)拉伸加卸載實驗中,未硫化橡膠同樣表現(xiàn)出明顯的應變率相關性,且在加載段的應變率相關性略大于卸載段的應變率相關性,如圖2(b)所示.此外,還可以看出,未硫化橡膠存在較為明顯的Mullins軟化現(xiàn)象,以及顯著的遲滯效應和殘余變形.Mullins軟化效應反映了試件在首次拉伸時的損傷,為了更清晰地觀察未硫化橡膠的軟化現(xiàn)象,圖3給出了其在不同應變率下隨拉伸循環(huán)次數(shù)(C)增加,拉伸上升段在前20%應變范圍內(nèi)的切變模量(G).可見,隨著循環(huán)次數(shù)增加,膠料切變模量明顯降低,且在較高應變率時降低速度更快.遲滯效應反映了循環(huán)加卸載實驗中的能量耗散,未硫化橡膠遲滯效應同樣受到應變率的影響,隨著應變率增加,遲滯效應有所增強.由于未硫化橡膠黏塑性較強,每次卸載后均存在一定的塑性殘余變形,且殘余變形量隨著循環(huán)次數(shù)的增加而累計增加,隨著應變率的增加有所降低.

圖3 不同應變率下4次上升段切變模量Fig.3 Shear modulus in four loading steps at different strain rates

2 未硫化橡膠黏彈塑性本構模型

基于對未硫化橡膠實驗結果的分析,提出了TN模型對未硫化橡膠的力學行為進行表征,模型的黏塑性響應基于對變形梯度的乘法分解,具體推導過程介紹如下.

2.1 變形梯度分解

對于不可壓縮材料,變形梯度F通過乘法分解為形變分量(Distortional Part)與體積分量(Dilatational Part)[17]:

F=FdisFdil

(1)

形變分量與體積分量分別具有以下形式[18]:

(2)

式中:J為變形梯度的行列式;I為單位張量;det(Fdis)=1.

此時左Cauchy-Green變形張量可表示為

B=FFT=J2/3B*

(3)

(4)

式中:B*為B的形變分量.

對于近似不可壓縮橡膠類材料黏塑性變形,形變分量也可分解為彈性變形分量和黏塑性分量[19]:

Fdis=FeFv

(5)

式中:Fe為變形梯度的彈性分量;Fv為變形梯度的黏塑性分量.

此時可獲得速度梯度Ldis:

(6)

(7)

(8)

2.2 TN模型

本文提出的TN模型包含3個并聯(lián)網(wǎng)絡,如圖4所示.其中一個網(wǎng)絡(網(wǎng)絡A)為超彈性網(wǎng)絡,用于表征未硫化橡膠的超彈性應力響應;另外兩個網(wǎng)絡為黏塑性網(wǎng)絡(網(wǎng)絡B、C),用于表征其黏塑性響應.不同于通常所采用的廣義Maxwell模型,黏塑性網(wǎng)絡B、C采用了基于高分子鏈微觀運動學機制的Bergstr?m-Boyce流動模型,能更好地表征未硫化橡膠的非線性黏塑性流動特征.同時,模型引入了基于分子鏈運動演化機制的損傷模型來表征未硫化橡膠的Mullins損傷軟化效應.

圖4 TN模型框架Fig.4 Frame of TN model

模型Cauchy應力σtol為3個網(wǎng)絡Cauchy應力的和:

σtol=σA+σB+σC

(9)

2.2.1TN模型應力響應 橡膠類材料的本構行為常通過應變能密度函數(shù)進行表征[1, 21],目前針對橡膠類高分子材料已經(jīng)發(fā)展出多種不同形式的應變能密度函數(shù),包括唯象應變能密度函數(shù)和基于分子網(wǎng)絡的應變能密度函數(shù)[22-24].其中,Arruda等[25]基于非高斯鏈分子網(wǎng)絡所提出的八鏈模型參數(shù)簡單,表征效果優(yōu)異,應用較為廣泛.考慮到未硫化橡膠具有相似的分子網(wǎng)絡構型,在此采用八鏈模型應變能密度函數(shù)來表征未硫化橡膠的超彈應力響應,應變能密度表達式為

(10)

(11)

對于超彈性力學響應,Cauchy應力可用應變能密度函數(shù)[18]表示為

(12)

式中:C為右Cauchy-Green變形張量.將應變能密度函數(shù)表達式(10)和(11)代入上式,并將C以不變量的形式寫入,即可獲得八鏈模型的Cauchy應力表達式[25]:

K(J-1)I

(13)

對于本文中三網(wǎng)絡框架,3個網(wǎng)絡的Cauchy應力可通過八鏈模型分別表示為

K(J-1)I

(14)

K(J-1)I

(15)

K(J-1)I

(16)

超彈性網(wǎng)絡(網(wǎng)絡A)的Cauchy應力由整體變形梯度決定;黏塑性網(wǎng)絡(網(wǎng)絡B、C)的變形梯度可進行乘法分解為彈性分量和黏塑性分量:F=FeBFvB,F=FeCFvC,其Cauchy應力僅由彈性分量提供,即

(17)

(18)

式中:FeB和FeC分別為B、C網(wǎng)絡的彈性變形梯度;FvB和FvC為B、C網(wǎng)絡的黏塑性變形梯度.

(19)

(20)

(21)

2.2.2黏塑性流動模型 要獲得黏塑性網(wǎng)絡B、C的彈性變形梯度,首先需要求解其黏塑性變形梯度,故需要合適的模型來對黏塑性流動進行表征.由1.2節(jié)實驗結果可知,未硫化橡膠力學行為較為復雜,通常的Maxwell模型不足以較好地表征其非線性黏塑性,故本文采用了基于高分子鏈微觀動力學機制的Bergstr?m-Boyce流動模型.該流動模型認為橡膠類聚合物材料的黏塑性變形是分子鏈的布朗運動(Brownian Motion)和“爬行”運動(Reptation Motion)聯(lián)合作用的結果[18],基于此假設,模型具有以下的表達形式:

(22)

(23)

(24)

將式(22)和(23)代入黏塑性變形速度梯度表達式(8),即可獲得黏塑性變形梯度的時間導數(shù):

(25)

網(wǎng)絡B、C中黏塑性變形梯度的時間導數(shù)分別為

(26)

(27)

對式(26)和(27)進行時間積分即可獲得B、C網(wǎng)絡的黏塑性變形梯度.

2.2.3Mullins損傷效應 此TN模型從橡膠高分子鏈的分子演化理論[26-27]出發(fā)來考慮未硫化橡膠的Mullins軟化效應.對于顆粒填充聚合物材料,填充顆粒和高分子鏈之間的聯(lián)接包括強化學聯(lián)接(共價鍵)及弱物理聯(lián)接[28].在拉伸過程中,物理聯(lián)接(填充顆粒與高分子鏈之間),以及部分化學聯(lián)接(高分子鏈之間或填充顆粒與高分子鏈之間)達到極限伸長而發(fā)生斷裂,引起軟化現(xiàn)象.對于未硫化炭黑填充橡膠,物理聯(lián)接的強度遠遠小于化學共價鍵聯(lián)接,故可認為拉伸過程中發(fā)生斷裂的化學共價鍵數(shù)量相對于物理聯(lián)接點很少,可以忽略不計[27].

TN模型中,A網(wǎng)絡代表了由物理聯(lián)接構成的炭黑顆粒填充高分子鏈網(wǎng)絡,用于表征未硫化橡膠的超彈性響應,而B、C黏塑性網(wǎng)絡代表了由化學共價鍵聯(lián)接構成的自由分子鏈網(wǎng)絡,用于表征未硫化橡膠的黏塑性響應.Mullins損傷主要來自于A網(wǎng)絡中填充顆粒與高分子鏈之間物理聯(lián)接的破壞,該破壞一方面減少了填充顆粒與高分子鏈之間的聯(lián)接點數(shù)目,引起分子網(wǎng)絡中N的增加,另一方面導致n減少[26].物理聯(lián)接的破壞程度與λcha有關,故n的減少由λcha決定,其變化情況可用A網(wǎng)絡中切變模量GA=nAkBT的整體變化來表示[27]:

(28)

N的變化同樣與λcha有關:

(29)

式中:C4～C9均為材料常數(shù).

將式(28)和(29)代入Cauchy應力表達式(14),可得包含應力軟化效應的炭黑填充分子鏈網(wǎng)絡的Cauchy應力表達式:

dev(B*)+K(J-1)I

(30)

化學共價鍵聯(lián)接構成的自由鏈網(wǎng)絡(網(wǎng)絡B、C)對Mullins效應的影響主要來自于分子鏈的運動演化回復,并表現(xiàn)為膠料模量的增大.在此采用微觀泡沫模型(MFM)[29]的形式來計及分子鏈運動的影響,并將原模型中的縮減密度改寫為應力回復因子R(取值應大于1,本文取為1.2)以體現(xiàn)分子鏈的運動回復對應力的增強效果,從而獲得了切變模量與體積模量的表達式:

(31)

(32)

式中:G0、K0為網(wǎng)絡初始切變模量和體積模量,G0=nkBT;α和h為模量縮放因子.

將式(31)和(32)代入Cauchy應力表達式(15)和(16),可獲得包含分子運動演化效應的B、C網(wǎng)絡Cauchy應力表達式:

(33)

(34)

式中:G0B、G0C和K0B、K0C分別為B、C網(wǎng)絡的初始切變模量和體積模量.

2.3 本構模型表征效果

曲線擬合即實驗曲線與表征曲線之間誤差最小化的過程,本文采用最小二乘法進行誤差計算:

(35)

表1 三網(wǎng)絡模型擬合參數(shù)Tab.1 Ultimate optimization parameters of NT model

模型參數(shù)擬合過程采用梯度下降法,擬合結果如圖5所示.圖中:εtrue為真應變;σtrue為真應力.可以看出,不同應變率下的擬合曲線與實驗曲線吻合良 好, 未硫化橡膠的復雜力學行為特征, 包括遲滯效應、殘余變形及軟化現(xiàn)象均得到了較好的表征.4種不同應變率下的表征結果在同一套擬合參數(shù)下完成,擬合參數(shù)如表1所示.

圖5 TN模型表征結果Fig.5 Fitting result using TN model

3 TN模型的有限元實現(xiàn)及驗證

獲得的本構模型若要應用于工程實踐,則需要將其進行數(shù)值實現(xiàn).本文依托于有限元軟件Abaqus,完成了TN模型用戶材料子程序(UMAT)的開發(fā),并進行了仿真計算以驗證模型的有效性和數(shù)值穩(wěn)定性.

模型有限元實現(xiàn)過程中采用的應力更新算法為二階精度的顯式龍格庫塔算法(Explicit Midpoint Rule)[32],該算法流程清晰,靈活性強,能夠較好地平衡計算精度、數(shù)值穩(wěn)定性以及易數(shù)值植入性.采用t代表前一個增量步,t+1代表后一個增量步,則該應力更新算法步驟如下.

(3) 對t時刻和t+1時刻總變形梯度平均獲得t+1/2時刻總變形梯度Ft+1/2.

(6) 檢驗B、C網(wǎng)絡的黏塑性變形梯度是否滿足誤差準則,若滿足則進入第(7)步,否則更新時間步長后返回第(2)步.

(7) 計算t+1時刻B、C網(wǎng)絡的彈性變形梯度并根據(jù)式(14)～(16)計算t+1時刻Cauchy應力.

(8) 儲存自定義狀態(tài)變量并計算一致切線剛度矩陣(Jacobian矩陣).

獲得TN模型用戶材料子程序后,為驗證其有效性,首先采用此子程序?qū)σ粋€單元在不同工況下(拉伸、壓縮、剪切)的力學行為進行了數(shù)值計算,各工況下單元的Mises應力(S)分布情況如圖6所示.圖7給出了在不同應變率下拉伸工況的仿真結果與測試結果的對比情況.可以看出,仿真結果能準確再現(xiàn)不同應變率下的拉伸應力應變曲線.

圖6 一個單元的仿真結果Fig.6 Simulation results of a single element

圖7 拉伸仿真結果與單調(diào)拉伸實驗結果對比Fig.7 Comparison of simulation results with experiments under uniaxial tensile conditions

圖8 多段松弛加卸載實驗與仿真結果對比Fig.8 Comparison of multistep tensile relaxation test and simulation result

另一方面,未硫化橡膠多用于成型工藝,此過程中膠料需經(jīng)歷較為復雜的變形,故本構模型在復雜變形工況下的計算穩(wěn)定性尤為重要.為了驗證TN模型的數(shù)值穩(wěn)定性,對輪胎胎面膠壓入模具的成型過程進行了數(shù)值仿真,如圖9所示.此模型實際上是對輪胎胎面膠縱向花紋溝成型工藝的簡化,如圖9(a)所示,模型軸對稱分布,且僅包含一條花紋溝.采用軸對稱有限元模型對此過程進行仿真,如圖9(b)所示,將胎面膠分為5層以便于觀察膠料流動,總厚度為20 mm,寬度為80 mm,采用的單元類型為CAX4H,節(jié)點數(shù)為 2 721 個,單元數(shù)為 2 594 個.模具采用離散剛體,深度為15 mm,中間的拱形突起高8.75 mm,寬7.5 mm.對胎面膠上層節(jié)點施加徑向向下的位移邊界條件,速度為0.5 mm/s,橡膠與模具之間的接觸摩擦因數(shù)設為0.5.胎面膠壓入模具后的各層膠料分布如圖9(c)所示,由于未硫化橡膠具有近似不可壓縮性,膠料充滿模具后由兩側溢出.圖9(d)和9(e)給出了此時刻胎面膠的應力和應變分布情況,可以看出,中間的拱形突起處膠料流動較為劇烈,存在一定的應力應變集中,靠近模具邊緣處的應力應變分布高于底部平坦部位.膠料流動情況以及應力應變分布情況與文獻[6, 8]中實驗及仿真結果基本一致.綜上可見,TN模型在復雜變形工況下仍能給出合理的計算結果,表明TN模型具有良好的數(shù)值穩(wěn)定性.

圖9 胎面膠壓入模具過程仿真Fig.9 Simplified molding process of a tire tread

4 結論

(2) 根據(jù)實驗結果,提出一個用于表征未硫化橡膠黏彈塑性力學行為的TN本構模型,模型由一個超彈性網(wǎng)絡和兩個非線性黏塑性網(wǎng)絡構成.3個網(wǎng)絡的應力響應均基于超彈八鏈模型,同時在黏塑性網(wǎng)絡中采用了Bergstr?m-Boyce流動法則描述膠料的非線性黏塑性流動.此外,模型從橡膠高分子鏈的分子演化理論出發(fā)考慮了未硫化橡膠的 Mullins 損傷軟化效應.采用TN模型對未硫化橡膠循環(huán)拉伸加卸載實驗曲線進行表征,表征結果與實驗結果吻合良好,膠料的各種拉伸力學特征均得到了較好的體現(xiàn).

(3) 依托有限元仿真軟件Abaqus,完成了TN模型用戶材料子程序(UMAT)的開發(fā),并對試件單調(diào)拉伸和多段松弛拉伸加卸載實驗進行了仿真,仿真結果與實驗結果吻合較好,驗證了此TN模型的有效性.最后,采用TN模型對輪胎胎面膠壓入模具的過程進行數(shù)值計算,進一步驗證了模型在復雜變形工況下良好的數(shù)值穩(wěn)定性.