拉/剪載荷下橡膠材料的Mullins誘導各向異性研究

李飛，丁智平，卜繼玲，黃友劍

（1.湖南工業大學 機械工程學院，湖南 株洲412007；2. 株洲時代新材料科技股份有限公司，湖南 株洲412007）

拉/剪載荷下橡膠材料的Mullins誘導各向異性研究

李飛1，2，丁智平1，卜繼玲2，黃友劍2

（1.湖南工業大學 機械工程學院，湖南 株洲412007；2. 株洲時代新材料科技股份有限公司，湖南 株洲412007）

使用矩形橡膠試樣依次進行不同位移載荷水平下純剪切和軸向拉伸/純剪切試驗，研究橡膠材料在拉/剪載荷作用下的Mullins軟化誘導各向異性。試驗結果表明，橡膠純剪切試樣在經過軸向拉伸后，由于Mullins軟化誘導各向異性的影響，再進行純剪切試驗會使試樣的橫向剛度增大。利用試樣的軸向拉伸/純剪切試驗數據，擬合得到相應的Mullins軟化控制參數，反映Mullins軟化誘導各向異性對剪切試樣橫向剛度的影響。分別對純剪切試樣在2種試驗載荷作用下的橫向剛度曲線進行有限元數值模擬，仿真估計值與試驗值的殘差平方和均值的對比結果表明，考慮Mullins軟化誘導各向異性可以顯著提高拉-剪載荷作用下橡膠元件橫向剪切剛度的預測精度。

橡膠；純剪切；Mullins軟化；誘導各向異性；有限元；剛度

0 引言

橡膠作為一種獨特的彈性阻尼材料而被廣泛地應用于各種隔振系統中[1]，其在承載過程中不僅承受很大的軸向載荷，同時還受到橫向剪切載荷作用，產生較大的剪切變形，例如軸箱彈簧、旁承和錐形簧等。為了確保橡膠彈性元件的使用性能，對其在剪切載荷作用下的橫向剛度的設計和預測，就顯得非常重要。

橡膠受載變形后會發生Mullins軟化效應，且這種軟化效應具有誘導各向異性。這是因為橡膠材料在拉伸過程中各個方向發生了不均勻的軟化，所以產生了各向異性[2]，且橡膠承受的載荷愈大，Mullins軟化誘導各向異性愈加顯著。利用有限元軟件對橡膠彈性元件的剛度進行數值模擬時，存在較大誤差，特別是剪切剛度預測與試驗結果存在很大差異。

對于Mullins軟化誘導各向異性的研究，國外學者從材料本構關系的角度進行了探索。R. Dargazany等人認為分子網絡結構的變化是橡膠材料在發生Mullins軟化時出現誘導各向異性的微觀機理，并提出模擬橡膠材料Mullins軟化各向異性的本構模型[3]。J. Diani等人通過試驗驗證了Mullins軟化誘導各向異性現象，并基于統計熱力學理論提出了相應的本構理論來描述這一現象[4]。M. Rebouaha等人基于唯象學理論，提出了計算硅膠材料各向異性軟化的本構理論[5]。以上學者所提出的本構理論雖然考慮了Mullins軟化的誘導各向異性，能夠提高計算精度，但是本構模型非常復雜，本構參數擬合難度很大，難以用于工程實際。本文通過對橡膠純剪切試樣進行多向載荷拉伸試驗，研究在剪切載荷作用下橡膠材料的Mullins誘導各向異性，基于偽彈性理論，提出一種Mullins效應軟化參數擬合方法，能夠方便地應用于工程有限元分析，為提高橡膠彈性元件受剪切載荷的橫向剛度預測精度提供設計依據和技術支持。

1 偽彈性理論本構模型

1.1 Mullins軟化效應

橡膠材料的Mullins軟化效應是指橡膠材料的加載路徑與卸載路徑不同，在相同位移載荷下，橡膠材料的卸載應力明顯低于加載應力。當再次施加小于初始位移水平的載荷時，橡膠材料的應力應變會沿著卸載曲線變化。當施加超過初始位移水平的載荷時，橡膠材料的應力應變會沿著新的加載曲線進行變化。圖1為發生Mullins軟化效應的加載-卸載應力應變曲線。[6]

圖1 理想的Mullins力學行為Fig.1 The ideal Mullins mechanical behavior

在工程實際中，橡膠彈性元件產品檢驗都要經過多次的加載卸載循環，以消除Mullins效應，然后再進行剛度檢測。而現有的工程有限元分析軟件中，經典的橡膠材料超彈本構模型沒有考慮Mullins效應，分析的結果往往與試驗檢測的結果相差很大。為了考慮Mullins效應的影響，提高計算精度，人們發展了偽彈性本構理論[7]。

1.2 偽彈性理論本構參數

基于唯象學的超彈材料本構模型，通過建立應變或拉伸比等加載變量與應變能之間的關系來描述橡膠材料的應力應變行為。為使本構方程能夠方便描述Mullins效應，將應變能密度分成畸變應變能密度和體積變化應變能密度2部分[8]，即

式中：U是應變能；Udev是材料初始加載至某一應變下的應變能密度；Uvol是材料體積變化所產生的應變能密度。在超彈模型中增加損傷變量來反映材料的Mullins軟化效應，并且將該變量與應變能函數聯系起來，偽彈性本構模型可以表示為[8]

式中：η為損傷變量，其物理意義為橡膠材料在某一變形下，其卸載應力與初始加載應力的比值，如圖2所示[8]，圖中是加載至相同位移下的應力值，是卸載至相同位移下的應力值；為損傷應變能密度函數；為主拉伸比。

圖2 損傷變量的物理意義Fig.2 The physical significance of damage variable

利用應變能計算損傷變量，其解析式為[9]

由式（3）可知，3個軟化控制參數都與材料軟化損傷的大小有著緊密的聯系。其中，r值決定了橡膠材料的最大損傷程度，r值越小，橡膠材料所能發生的軟化程度就越大；的值與損傷的大小成反比。值越大，軟化損傷越小；m為材料軟化所需最大應變能的閥值，m值越大；材料發生軟化所需的最大應變能的閥值越大，m值的大小對低應變水平加載的軟化影響非常大。

偽彈性本構理論沒有考慮橡膠的Mullins軟化誘導各向異性，將Mullins軟化損傷視為是各向同性的。對橡膠彈性元件在單向受載下同向剛度進行仿真時，偽彈性本構理論具有較高的計算精度，但是在預測受多向載荷歷程的橡膠產品剛度時，仿真結果與試驗值有很大差別。

2 拉/剪載荷下Mullins誘導各向異性試驗

2.1 試驗原理

大多數橡膠減振元件的使用工況，都是在承擔較大軸向載荷的同時還承受橫向載荷。某型空氣彈簧輔助彈簧如圖3所示，該元件在實際使用過程中除承受96.9 kN軸向載荷外，還在頂部有10 mm的橫向位移。在這種工況下，橡膠材料同時受到軸向載荷和橫向剪切載荷作用。

圖3 輔助彈簧承載后變形情況Fig.3 The deformation of auxiliary spring under load

課題組設計了2組試驗來進行對比，第1組為試驗A，第2組為試驗B，試驗流程如圖4所示。

試驗A為純剪切試驗，由于試樣的寬度遠大于其拉伸方向的長度，試樣在拉伸過程中寬度方向的拉伸比近似等于1，因而處于純剪切狀態。

試驗B是先將試樣進行軸向拉伸，然后進行純剪切試驗。橡膠元件通常承受壓縮載荷，但是由于橡膠材料的不可壓縮性，在受到軸向的壓縮載荷時，會在橫向發生膨脹，產生很大的拉應力。因此使用軸向拉伸/純剪切試驗（試驗B，下文中簡稱拉/剪試驗）可以模擬橡膠彈性元件承受軸向壓縮和橫向剪切載荷歷程的工況，研究橡膠材料在拉/剪載荷下的Mullins誘導各向異性。

圖4 試驗示意圖Fig.4 Schematic of experiments

2.2 試驗方案

試驗設備為SANS微機控制萬能拉伸試驗機，試驗溫度為24℃，加載速率為5 mm/min。橡膠試樣尺寸為60 mm×26 mm×2 mm，剪切試驗標距為60 mm× 6 mm×2 mm，試驗方案見表1。

表1 Mullins誘導各向異性試驗方案Table 1Experiment scheme of Mullins induced anisotropy

2.3 試驗結果與分析

試驗結果如圖5～7所示。由試驗結果可以看出，剪切橡膠試樣在經過軸向拉伸以后，橫向剛度出現了明顯的差異（圖中橡膠純剪切試樣的橫向剛度計算是用試樣在橫向所受載荷除以橫向位移），而且這種差異會隨著軸向拉伸位移載荷水平的增大而增大。由于偽彈性本構理論將Mullins軟化損傷看成各向同性的，所以難以準確計算橡膠剪切試樣經軸向拉伸變形后的橫向剛度。線和試驗B剛度曲線的擬合值。

圖5 6 和 18 mm 位移載荷下剪切試樣橫向剛度對比Fig.5 The lateral stiffness contrast under the displacement load of 6 mm and 18 mm

圖6 12和 24 mm拉伸位移下剪切試樣橫向剛度對比Fig.6 The lateral stiffness contrast under the tensiondisplacement of 12 mm and 24 mm

圖7 21 mm拉伸位移下各向異性對比Fig.7 Anisotropy contrast under tension displacement of 21 mm

圖8 各級位移載荷水平下試驗A和試驗B軟化曲線擬合Fig.8 The soften curves fitting of experiment A and experiment B under different displacement loading levels

為定量描述Mullins軟化誘導各向異性的大小隨位移載荷水平變化的規律，計算各個位移載荷水平下的試驗A剛度曲線和試驗B剛度曲線的殘差平方和均值，計算式如下[11]，即

式中：yi為試驗A中每個位移載荷水平下載荷；為試驗B中每個位移水平下加載至相同位移下的載荷；n為所有的數據點的個數。

表2所示為計算結果，可以看出，隨著位移水平的增大，Mullins軟化誘導各向異性對于橫向剛度的影響越發顯著。

表2 殘差平方和均值計算結果Table 2The calculated results of mean sum of squared residuals

3 Mullins軟化誘導各向異性的數值模擬

3.1 軟化控制參數擬合

在Abaqus軟件中擬合軟化控制參數，首先將最大位移載荷的純剪切試驗數據輸入到超彈材料本構模型中，確定試樣純剪切的加載曲線；然后，將每個位移載荷水平純剪切穩定循環試驗數據輸入到偽彈性本構模型中。為了分析Mullins軟化誘導各向異性的影響，應分別擬合純剪切試樣的剪切軟化控制參數與拉/剪軟化控制參數。在擬合剪切軟化控制參數時，輸入試驗A每個位移載荷水平的穩定循環數據；而在擬合拉/剪軟化控制參數時，輸入試驗B的每個位移載荷的穩定循環數據。

在Abaqus軟件中建立試樣的有限元模型如圖9所示。

圖9 有限元仿真邊界條件Fig.9 The boundary conditions for finite element simulation

按照試驗情況設置邊界條件，并選擇Data Check命令，在Abaqus軟件中自動生成的DAT文件中輸出軟化控制參數。利用Abaqus軟件擬合得出的軟化控制參數見表3。

表3 有限元擬合參數Table 3The finite element fitting parameters

3.2 Mullins軟化有限元仿真與試驗數據對比

由于實際的Mullins軟化曲線與理想化的軟化曲線存在一些差別，即在循環加載試驗中，試樣加載曲線與卸載曲線并不會完全重合。為方便仿真結果與試驗數據對比，本文利用最小二乘法將試驗中的載荷-位移穩定循環加載試驗數據擬合成一條五次多項式曲線，用以檢驗仿真效果。擬合結果為

其中：y剪和y軸分別為純剪切試驗和拉/剪試驗的載荷值；x為位移。

在24 mm位移載荷水平下，純剪切和拉/剪試驗數據擬合曲線如圖10～11所示。純剪切試驗沒有考慮Mullins軟化誘導各向異性效應，拉/剪試驗則考慮了Mullins軟化誘導各向異性效應。

從圖10～11中可以看出，不考慮Mullins軟化效應的有限元分析值與試驗結果相差較大，而考慮了軟化效應的分析值與試驗結果吻合較好。

圖12所示為考慮誘導各向異性與不考慮誘導各向異性仿真曲線對比。可以看出，考慮Mullins軟化誘導各向異性擬合出的軟化控制參數，能夠更好的模擬橡膠試樣在拉/剪載荷作用下橫向Mullins軟化效應。殘差分析表明，沒有考慮誘導各向異性的仿真軟化曲線與試驗值的殘差平方和均值為932.39，而考慮誘導各向異性模擬出的軟化曲線與試驗值的殘差平方和均值為253.47，計算精度顯著提高。

圖10 24 mm純剪切試驗數據擬合曲線與仿真曲線對比Fig.10 The contrast of fitting curve and simulation curve ofpure shear experiment data of 24 mm

圖11 24 mm拉/剪試驗數據擬合曲線與仿真曲線對比Fig.11 The contrast of fitting curve and simulation curve of tension/shear experiment data of 24 mm

圖12 考慮誘導各向異性與不考慮誘導各向異性仿真曲線對比Fig.12 The simulation curves contrast of considering induced anisotropy and without considering induced anisotropy

4 結論

1）使用矩形試樣進行5個位移載荷水平的純剪切和軸向拉伸/純剪切循環加載試驗，結果表明，橡膠材料在發生Mullins軟化的同時，還具有Mullins軟化誘導各向異性，并且這種軟化誘導各向異性會隨軸向位移載荷水平的增加而增大。

2）利用橡膠試樣的純剪切和軸向拉伸/純剪切循環加載試驗數據，分別擬合得出橡膠試樣在兩種載荷歷程下的橫向Mullins軟化控制參數，提出了橡膠材料受拉/剪載荷的橫向剛度預測方法。

3）利用擬合軟化控制參數對橡膠試樣在2種載荷歷程下的橫向剛度進行有限元仿真，并對仿真估計值與試驗值的殘差平方和均值進行對比，結果表明，考慮和沒有考慮Mullins 軟化誘導各向異性的殘差平方和均值分別為253.47和932.39，仿真精度明顯得到提高。

[1]方建輝，丁智平，卜繼玲，等. 橡膠減振支座動態性能仿真分析研究[J]. 湖南工業大學學報，2014，28(2)：16-22. Fang Jianhui，Ding Zhiping，Bu Jiling. Simulation Research on Dynamic Performance of Rubber Damping Brackets of Wind Turbine[J]. Journal of Hunan University of Technology，2014，28(2)：16-22.

[2]Diani J，Fayolle B，Gilormini P. A Review on the Mullins Effect[J]. European Polymer Journal，2009，45(3)：601-612.

[3]Dargazany R，Itskov M. A Network Evolution Model for the Anisotropic Mullins Effect in Carbon Black Filled Rubbers[J]. International Journal of Solids and Structures，2009，46(16)：2967-2977.

[4]Diani J，Brieu M，Vacherand J M. A Damage Directional Constitutive Model for Mullins Effect with Permanent Set and Induced Anisotropy[J]. European Journal of Mechanics，2006，25(3)：483-496.

[5]Rebouah M，Machad G，Chagnon G，et al. Anisotropic Mullins Stress Softening of a Deformed Silicone Holey Plate [J]. Mechanics Research Communications，2013，49：36-43.

[6]Cantournet S， Desmorat R， Besson J. Mullins Effect and Cyclic Stress Softening of Filled Elastomers by Internal Sliding and Friction Thermodynamics Model[J]. International Journal of Solids and Structures，2009，46(11/12)：2255-2264.

[7]Ogden R W，Roxburgh D G. A Pseudo-Elastic Model for the Mullins Effect in Filled Rubber[J].Proceedings of the Royal Society，1999，455：2861-2877.

[8]張友南，黃友劍. 多軸承載下基于Mullins效應的非線性剛度的分析方法[J]. 計算機輔助工程，2013，22(增刊2)：360-363. Zhang Younan，Huang Youjian. Analysis Method for Non-Linear Stiffness Under Multi-Axis Loading Condition Basedon Mullins Effect[J]. Computer Aided Engineering，2013，22(S2)：360-363.

[9]Ryan E P，Will V M. Implications of the Mullins Effect on the Stiffness of a Pre-Loaded Rubber Component ABAQUS Users’ Conference[EB/OL]. 2014-07-15. http：//www. axelproducts.com/downloads/PaigeMars.Pdf.

[10]Hoffman J D. Numerical Methods for Engineers and Scientists [M]. New York：McGraw Hill Inc，1992：378.

[11]王巖，隋思漣，王愛青. 數理統計與MATLAB工程數據分析[M]. 北京：清華大學出版社，2006：165. Wang Yan，Sui Silian，Wang Aiqing. Mathematical Statistics and MATLAB Engineering Data Analysis[M]. Beijing：Tsinghua University Press，2006：165.

（責任編輯：申劍）

Study on Mullins Softening Induced Anisotropy of Rubber Materials Under Tension/Shear Loads

Li Fei1,2，Ding Zhiping1，Bu Jiling2，Huang Youjian2
（1. School of Mechanical Engineering,Hunan University of Technology,Zhuzhou Hunan 412007,China；2. Zhuzhou Time New Material Technology Co.,Ltd.,Zhuzhou Hunan 412007,China）

Conducts pure shear and axial tension/ pure shear experiments on rectangular rubber specimens under different displacement load levels,and investigates Mullins softening induced anisotropy of rubber materials under tension/ shear load. The experimental results show that because of the effect of Mullins softening induced anisotropy,the pure shear test on the specimens after axially tensions causes the specimen lateral stiffness increasing. Uses the axial tension/ pure shear experiment data of rubber specimens to fit and obtain the relevant Mullins softening control parameters,and reflects the effect of Mullins softening induced anisotropy on the transverse stiffness. Makes finite element simulation on he transverse stiffness curves of the specimens under two kinds of load tests,and compares the mean sum of squared residuals between the simulation estimate value and the experimental value. It is considered that Mullins softening induced anisotropy will significantly increase the prediction accuracy of transverse shear stiffness of rubber elements under tension and shear load.

rubber；pure shear；Mullins softening；induced anisotropy；finite element；stiffness

Q334.9

A

1673-9833(2015)01-0038-06

2014-11-02

中央財政創新團隊專項基金資助項目（0420036017），湖南省十二五機械工程重點學科基金資助項目（湘教發[2011]-76號）

李飛（1986-），男，遼寧丹東人，湖南工業大學碩士生，主要研究方向為機械結構強度理論及應用，

E-mail：woshitiancailifei@163.com

丁智平（1956-），男，湖南桃江人，湖南工業大學教授，博士，主要研究方向為機械結構強度理論及應用，E-mail：472447034@qq.com

10.3969/j.issn.1673-9833.2015.01.007