楊海波，劉 楓，李凡珠，盧詠來，張立群

（北京化工大學 有機無機復合材料國家重點實驗室，北京 100029）

采用有限元方法研究橡膠材料/制品的力學響應行為已成為一種方便有效的方法[1-5]。有限元分析中，精確的計算結果與可靠的材料模型和準確的材料參數(shù)密切相關，即本構方程的選取非常重要。在有限元分析中[6-10]較少提及本構方程選取的過程、依據(jù)以及準確性。因此系統(tǒng)地研究和總結不同應力-應變試驗數(shù)據(jù)及不同本構方程對橡膠材料/制品有限元分析結果的影響很有必要。

本工作以炭黑填充天然橡膠（NR）制作的圓柱形橡膠試樣的壓縮變形過程為分析對象，分別采用不同的本構方程對三類應力-應變試驗曲線進行擬合，并用有限元方法預測其壓縮變形。通過對壓縮變形試驗數(shù)據(jù)的對比，判斷不同應力-應變曲線及本構方程的計算精度。

1 實驗

1.1 試樣制備

膠料配方如下：NR 100，炭黑N234 40，氧化鋅 3，硬脂酸 1，硫黃 1，促進劑CBS 1.5。使用密煉機將各組分按一定加料順序在NR中混合和分散均勻，然后將混煉膠停放一定時間，使用圓盤式硫化儀測定硫化曲線，得到t90，再經(jīng)平板硫化機硫化和氣動裁壓機沖壓等步驟分別得到單軸拉伸（UT）試樣、平面拉伸（PT）試樣、等雙軸拉伸（ET）試樣[11]以及圓柱形單軸壓縮試樣。其中UT試樣尺寸為50 mm×5.6 mm×2 mm，PT試樣尺寸為60 mm×10 mm×2 mm，圓柱形單軸壓縮試樣[12]直徑為17.8 mm，高為25 mm，如圖1所示。

圖1 各類試樣示意

1.2 本構（應力-應變）曲線測試

使用電子萬能試驗機對橡膠試樣做多次拉伸回復試驗，以消除橡膠材料的應力軟化效應（Mullins效應），得到穩(wěn)定的力-位移數(shù)據(jù)。再將力-位移數(shù)據(jù)轉為應力-應變數(shù)據(jù)。在進行有限元分析時，選取以下3類應力-應變試驗曲線來擬合橡膠材料的超彈性本構模型：（1）UT應力-應變曲線；（2）單軸壓縮（UC）應力-應變曲線；（3）UT，PT，ET組合的應力-應變曲線。對于PT試驗，由于試樣垂直于拉伸方向上的橫截面積比單軸拉伸大得多，導致試樣與夾具之間易產(chǎn)生滑脫，尤其在大變形加載時。為解決這一問題，采用平面拉伸試驗結合光學測試技術，在試樣表面涂上銀粉，即可輕松地追蹤試樣的真實變形情況，從而修正由于試樣的滑脫而導致力學數(shù)據(jù)不準確的問題。該方法類似于激光測距儀的功能。ET試驗數(shù)據(jù)測定使用易瑞博科技（北京）有限公司的精密測試設備[11]。

1.3 試驗驗證數(shù)據(jù)測試

在有限元分析中，計算的是圓柱形橡膠試樣單軸下壓4.45 mm過程中對應的反作用力。為驗證應力-應變曲線以及本構方程的計算精度，使用拉力機得到圓柱形橡膠試樣的力-位移曲線。為保證試驗數(shù)據(jù)的可靠性，單軸壓縮試驗進行4次，測試結果曲線高度重合。

2 有限元模型

2.1 幾何模型

以炭黑填充NR材料制作的圓柱形單軸壓縮試樣為建模對象，采用軸對稱模型建模，如圖2所示，其中虛線為對稱軸。

圖2 圓柱形試樣軸對稱幾何模型

2.2 材料模型

在材料模型的選取上，對3類應力-應變數(shù)據(jù)采用不同超彈性本構方程擬合，超彈性本構模型包括：一階和二階多項式模型（P_N1，P_N2）、一階至六階減縮多項式模型（RP_N1，RP_N2，RP_N3，RP_N4，RP_N5，RP_N6）、一階至六階Ogden模型（O_N1，O_N2，O_N3，O_N4，O_N5，O_N6）以及Arruda Boyce（AB），Marlow（Mar），Van Der Waals（VdW）模型共6類17種[7，13]，共有51個組合計算結果。

2.3 邊界條件

在邊界條件的施加上，有限元模型中載荷的施加與圓柱形橡膠試樣的壓縮變形試驗完全相同，試樣的上下表面建立兩個剛性面來模擬試驗設備的上下兩個夾持面，載荷施加在兩個剛性面上。剛性面與試樣間建立接觸。其中下端剛性面固定，上端剛性面受到Y軸負方向的4.45 mm的位移加載。整個模型共包含1 150個（CAX8RH）單元、1 224個節(jié)點。

3 結果與討論

3.1 基于UT應力-應變數(shù)據(jù)的超彈性本構方程

3.1.1 方程擬合

不同本構方程對UT應力-應變試驗數(shù)據(jù)的擬合結果見圖3。由圖3可見：P_N2，Mar以及二階以上的O_Ni方程對于UT試驗數(shù)據(jù)的擬合精度較高；P_N1，VdW，AB，O_N1以及三階以下的RP_Ni方程的擬合數(shù)據(jù)小應變下擬合值低于試驗值，中等應變下擬合值高于試驗值，而大應變下擬合值遠低于試驗值；三階及以上的RP_Ni方程擬合效果較好，但擬合曲線不平滑，存在波動，穩(wěn)定性不佳。

圖3 不同本構方程對UT應力-應變試驗數(shù)據(jù)的擬合結果

3.1.2 有限元分析結果

UT試驗數(shù)據(jù)擬合的超彈性本構方程對應的圓柱形試樣的變形計算結果見圖4。由圖4可見：在使用UT應力-應變試驗數(shù)據(jù)擬合本構方程時，擬合精度極高的O_Ni方程（二階及以上）和P_N2方程的預測結果與試驗數(shù)據(jù)差別較大，尤其在大位移變形情況下；而擬合精度極高的Mar方程在整個變形范圍內(nèi)的計算值與試驗值接近；P_N1，VdW和O_N1方程的預測結果類似，在較小位移下（3 mm以內(nèi)）的計算值與試驗值接近，在較大位移下的計算值高于試驗值；RP_N1和AB方程的預測結果類似，在小位移下的計算值低于試驗值，而當位移增大時（5～7 mm），計算值與試驗值接近；RP_N2方程在小位移下（3 mm以內(nèi)）的預測值略低于試驗值，隨著加載位移的增大，其計算值與試驗值吻合；三階及以上的RP_Ni方程的整體預測結果與試驗數(shù)據(jù)接近，但隨著方程階數(shù)的提高，材料模型的穩(wěn)定性下降，計算出現(xiàn)收斂問題。

圖4 UT試驗數(shù)據(jù)擬合的超彈性本構方程對應的圓柱形試樣的變形計算結果

3.2 基于UC應力-應變數(shù)據(jù)的超彈性本構方程

3.2.1 方程擬合

不同本構方程對UC應力-應變試驗數(shù)據(jù)的擬合結果見圖5。由圖5可見：P_N1，P_N2，VdW，Mar和二階及以上的O_Ni方程對UC試驗數(shù)據(jù)的擬合精度較高；AB，O_N1和RP_N1方程在較小應變下（壓縮應變小于0.3）擬合值低于試驗值，較大應變下擬合值高于試驗值；二階及以上的RP_Ni方程擬合效果較好，但擬合曲線同樣出現(xiàn)不平滑的現(xiàn)象，擬合曲線存在波動，穩(wěn)定性不佳。

圖5 不同本構方程對UC應力-應變試驗數(shù)據(jù)的擬合結果

3.2.2 有限元分析結果

UC試驗數(shù)據(jù)擬合的超彈性本構方程對應的圓柱形試樣的變形計算結果見圖6。由圖6可見：使用UC應力-應變數(shù)據(jù)擬合本構方程時，擬合精度較高的P_N1，P_N2，VdW，Mar和二階及以上的O_Ni方程的預測結果與試驗數(shù)據(jù)接近，尤其是Mar方程；AB，O_N1和RP_N1方程在加載位移范圍內(nèi)計算值均低于試驗值；二階及以上的RP_Ni方程預測精度較高，尤其是RP_N5方程；RP_N2，RP_N4和RP_N6方程由于材料模型的不穩(wěn)定性導致計算出現(xiàn)收斂問題。由此可見，相對于UT試驗數(shù)據(jù)，采用UC試驗數(shù)據(jù)擬合的超彈性本構方程預測圓柱形橡膠試樣的單軸壓縮變形行為計算精度大幅度提高。

圖6 UC試驗數(shù)據(jù)擬合的超彈性本構方程對應的圓柱形試樣的變形計算結果

3.3 基于3種應力-應變數(shù)據(jù)的超彈性本構方程

3.3.1 方程擬合

不同本構方程對UT，PT，ET組合的應力-應變試驗數(shù)據(jù)的擬合結果見圖7—9。由圖7可知：P_N1和AB方程與RP_N1方程的擬合結果類似；P_N2方程對ET試驗數(shù)據(jù)的擬合效果較佳，但卻不能很好地描述UT和PT數(shù)據(jù)；VdW方程在應變小于0.7時的整體擬合精度較高；Mar方程在應變小于0.4時可以很好地描述3類應力-應變曲線，當應變大于0.4時，ET的擬合數(shù)據(jù)開始高于相應的試驗數(shù)據(jù)，但對于UT和PT數(shù)據(jù)，擬合數(shù)據(jù)與試驗數(shù)據(jù)完全重合。由圖8可知，RP_N1和RP_N2方程的整體擬合精度不理想，三階以上的RP_Ni方程在較大應變時出現(xiàn)不穩(wěn)定問題，RP_N3方程的擬合精度相對較高，值得一提的是不同階數(shù)的RP_Ni方程在較小應變下（應變小于0.2）的擬合精度均比較理想。由圖9可見，三階及以下的O_Ni方程對UT，PT及ET組合的應力-應變試驗數(shù)據(jù)的擬合精度較差，三階以上的O_Ni方程可以很好地描述UT和ET試驗數(shù)據(jù)，但當應變大于0.6時，PT的擬合數(shù)據(jù)與UT擬合數(shù)據(jù)重合，總體而言，O_N5和O_N6方程的擬合精度較高。

圖7 多項式和AB，Mar，VdW模型本構方程對UT，PT及ET組合的應力-應變試驗數(shù)據(jù)的擬合結果

圖8 減縮多項式模型本構方程對UT，PT及ET組合的應力-應變試驗數(shù)據(jù)的擬合結果

圖9 Ogden模型本構方程對UT，PT及ET組合的應力-應變試驗數(shù)據(jù)的擬合結果

3.3.2 有限元分析結果

采用UT，PT及ET組合的應力-應變數(shù)據(jù)擬合本構方程進行有限元分析時，整體計算結果的精度進一步提高。UC試驗數(shù)據(jù)擬合的超彈性本構方程對應的圓柱形試樣的變形計算結果見圖10。圖10（a）表明，AB，P_N1和P_N2方程的計算結果低于試驗數(shù)據(jù)，VdW方程除了在小位移（小于1 mm）之外計算結果與試驗數(shù)據(jù)相吻合，Mar方程在加載位移范圍內(nèi)的計算精度均較高。由圖10（b）可知，RP_N1和RP_N2方程的計算結果低于試驗數(shù)據(jù)，三階及以上的RP_Ni方程計算精度較高，但仍存在小位移（小于1 mm）下計算精度較低的問題。由圖10（c）可知，O_N1，O_N2和O_N3方程的計算結果低于試驗數(shù)據(jù)，而O_N5和O_N6方程的計算精度高于較低階數(shù)的O_Ni方程，且其計算結果與VdW方程類似，即小位移（小于1 mm）下的計算精度較低。

圖10 UT，PT，EC組合的試驗數(shù)據(jù)擬合的超彈性本構方程對應的圓柱形試樣的變形計算結果

3.4 不同試驗數(shù)據(jù)及超彈性本構方程的有限元分析結果對比

通過上述3類試驗數(shù)據(jù)擬合橡膠材料的超彈性本構方程計算圓柱形橡膠試樣的壓縮變形結果可知，僅提供UT應力-應變數(shù)據(jù)擬合本構方程時，不能根據(jù)擬合精度判斷各個方程的優(yōu)劣，如擬合精度極高的O_Ni方程（二階及以上）和P_N2方程的計算結果與試驗數(shù)據(jù)差別極大。同樣，高擬合精度的Mar方程在整個變形范圍內(nèi)的計算值與試驗值接近，而擬合精度較差的RP_N2方程對應的計算結果卻與試驗數(shù)據(jù)接近。通過UC應力-應變數(shù)據(jù)擬合本構方程計算圓柱形橡膠試樣的壓縮變形行為時，其整體計算精度遠高于通過UT數(shù)據(jù)擬合本構方程，而且擬合精度較高的本構方程對應的計算結果的精度亦較高。

相對于僅采用UT試驗數(shù)據(jù)或UC試驗數(shù)據(jù)，同時使用UT，PT以及ET試驗數(shù)據(jù)擬合本構方程預測圓柱形橡膠試樣的單軸壓縮變形的整體計算精度進一步大大提高，如Mar，VdW以及高階的O_Ni和RP_Ni方程的計算結果與試驗數(shù)據(jù)吻合度極高。這主要是因為UT對應純拉伸狀態(tài)，PT對應純剪切狀態(tài)，ET對應純壓縮狀態(tài)。本工作中，雖然圓柱形橡膠試樣處于單軸壓縮的加載工況，但其并非處在純壓縮狀態(tài)，如圖11所示。由圖11可見，除了縱向壓縮變形外[圖11（b）]，該試樣亦存在橫向拉伸[圖11（a）]和周向拉伸變形[圖11（c）]，而圖11（d）表明由于橡膠試樣上下表面與上下加載面之間存在約束，區(qū)域A處的剪切變形較為明顯。

圖11 圓柱形橡膠試樣的名義應變各分量

4 結論

對于應力-應變試驗數(shù)據(jù)，僅采用單軸拉伸試驗數(shù)據(jù)擬合本構方程時，不能單純依賴擬合精度判斷本構方程的優(yōu)劣；而采用單軸壓縮試驗數(shù)據(jù)擬合本構方程時，其整體計算精度大大提升；同時采用單軸拉伸、平面拉伸和等雙軸拉伸試驗數(shù)據(jù)擬合本構方程時，計算精度進一步提升，且可依據(jù)擬合精度判斷相應超彈性本構方程的優(yōu)劣。

若僅提供一種變形狀態(tài)下的應力-應變數(shù)據(jù)擬合本構方程，不能只根據(jù)擬合精度判斷各方程的優(yōu)劣，如采用單軸壓縮試驗數(shù)據(jù)擬合本構方程預測橡膠試樣單軸壓縮力學行為的整體表現(xiàn)均較好；若提供單軸拉伸、平面拉伸和等雙軸拉伸試驗數(shù)據(jù)擬合本構方程時，可以選用擬合精度高的O_Ni，P_Ni，VdW，Mar等方程作為橡膠材料的超彈性本構方程；而在僅有單軸拉伸數(shù)據(jù)的情況下，選用RP_Ni，AB，Mar方程較好，不要選用擬合精度高的O_Ni和P_N2方程。應變加載條件不同，同種本構方程在不同應變區(qū)間內(nèi)的計算精度亦有所差別。

要準確預測處于復雜變形條件下的橡膠材料及制品的變形行為，需提供相應橡膠材料單軸拉伸、平面拉伸和等雙軸拉伸試驗數(shù)據(jù)，并選用擬合精度高的本構方程作為橡膠材料模型。