橡膠隔振器靜態(tài)力－位移關(guān)系計算方法的研究

何小靜，上官文斌

(華南理工大學(xué) 機械與汽車工程學(xué)院，廣州 510641)

橡膠材料作為一種高性能的超彈性材料，廣泛地應(yīng)用于汽車隔振器，如汽車動力總成懸置、車身懸置、排氣管吊耳、底盤襯套等。橡膠材料是一種大變形的非線性超彈性材料，對其進行力－位移分析時，其本構(gòu)模型的選擇和本構(gòu)模型常數(shù)的獲取是一項重要的工作。目前可用于橡膠隔振器靜態(tài)力－位移關(guān)系計算的本構(gòu)模型的較多，如Mooney-Rivlin、Ogden等。不同本構(gòu)模型的常數(shù)不一樣，如何合理的選擇本構(gòu)模型和開展哪些實驗以獲取本構(gòu)模型的常數(shù)是橡膠隔振器靜態(tài)力－位移分析的重要工作。

Seibert等［1］對 Arruda-Boyce 模型、Van Der Waals模型、Yeoh模型三種本構(gòu)模型進行對比研究。黃建龍等［2］利用Mooney-Rivlin模型和Yeoh模型對橡膠材料進行有限元分析，介紹兩種本構(gòu)模型材料常數(shù)的獲取方法。王麗榮等［3］利用Mooney-Rivlin模型和Ogden模型對橡膠隔振器靜態(tài)力－位移特性進行有限元分析，介紹了單元特性和網(wǎng)格質(zhì)量對橡膠隔振器靜態(tài)特性的影響。Gracia等［4］研究了一種應(yīng)用于工業(yè)橡膠材料彈塑性特性分析的overlay模型，并介紹其常數(shù)擬合過程。Stefan Hartmann等［5］研究了一種獲取橡膠本構(gòu)模型材料常數(shù)的新方法－光學(xué)測量法，確定采用光學(xué)測量法獲取材料常數(shù)的準(zhǔn)確性。

本文選取了目前廣泛應(yīng)用的幾種本構(gòu)模型來研究橡膠隔振器力－位移關(guān)系計算結(jié)果的影響因素。為獲得各本構(gòu)模型的材料常數(shù)，開展了橡膠材料在不同應(yīng)力－應(yīng)變狀態(tài)(單軸拉伸、等雙軸拉伸和平面剪切)和不同最大應(yīng)變下的應(yīng)力－應(yīng)變關(guān)系的測試。利用所測得的應(yīng)力－應(yīng)變關(guān)系，和不同本構(gòu)模型得到的理論應(yīng)力－應(yīng)變關(guān)系，根據(jù)最小二乘法獲得不同本構(gòu)模型中的材料常數(shù)。利用獲得的材料常數(shù)，對兩種橡膠隔振器(橡膠懸置1，2)的靜態(tài)力－位移特性進行了計算并和實測值進行了對比分析。

分析結(jié)果表明，利用同一本構(gòu)模型計算隔振器的靜剛度時，對懸置1，在拉壓變形的方向，采用三種應(yīng)力－應(yīng)變狀態(tài)擬合得到的本構(gòu)模型常數(shù)計算得到靜剛度值，比單軸拉伸應(yīng)力－應(yīng)變狀態(tài)下的計算精度高;對只有剪切變形的方向，用平面剪切應(yīng)力－應(yīng)變狀態(tài)擬合得到的本構(gòu)模型常數(shù)計算得到的靜剛度，比三種應(yīng)力－應(yīng)變狀態(tài)下的計算精度高。對懸置2，在同一本構(gòu)模型下，利用單軸拉伸、等雙軸拉伸和平面剪切三種應(yīng)力－應(yīng)變狀態(tài)擬合得到的本構(gòu)模型常數(shù)計算出的靜剛度值，比在兩種應(yīng)力－應(yīng)變狀態(tài)下計算出的靜剛度值相對誤差更小。

以橡膠懸置2為研究對象，分析了不同本構(gòu)模型對橡膠隔振器靜態(tài)力－位移特性計算結(jié)果的影響。結(jié)果表明，Mooney-Rivlin模型的計算精度最高，其相對誤差均小于10%，其次是Arruda-Boyce模型，計算相對誤差最大的是Marlow模型。

1 橡膠隔振器力－位移分析的本構(gòu)模型

對橡膠隔振器進行靜態(tài)特性計算時，假定橡膠材料為各向同性的超彈性材料［1］。通過對其應(yīng)變能密度函數(shù)對應(yīng)變不變量求導(dǎo)，可以得到材料的工程應(yīng)力與工程應(yīng)變之間的本構(gòu)關(guān)系［6］。工程應(yīng)力和工程應(yīng)變是忽略因施加載荷的增加或者減少而引起的橫截面面積的變化所得到的應(yīng)力、應(yīng)變［15］。

超彈性材料的應(yīng)變能密度函數(shù)有多種形式，如Mooney-Rivlin 模型［2］、Van Der Waals 模型［9］、Marlow模型［9］、Ogden 模型［10］、Yeoh 模型［2］、Arruda-Boyce 模型［3］、Neo-Hookean 模型、Ploynomial模型［7］等。本文只討論目前廣泛應(yīng)用的前六種本構(gòu)模型，及其常數(shù)的獲取。

應(yīng)變能密度函數(shù)的一般表達式為［8］:

其中，I1、I2、I3分別為一階、二階、三階應(yīng)變不變量，它們?yōu)槿齻€主拉伸比的函數(shù);Ci(i=1，2，…，m)為m個表示超彈性材料剪切特性的常數(shù)，dj(j=1，2，…，n)為n個表示超彈性材料壓縮特性的常數(shù)。I1、I2、I3與超彈性材料的三個主拉伸比λ1、λ2、λ3的關(guān)系為:

由文獻［2－3］、文獻［9－10］中對應(yīng)的各種應(yīng)變能密度函數(shù)，對三個主拉伸比 λ1、λ2、λ3分別求導(dǎo)數(shù)，就可以計算出材料在不同變形狀態(tài)時的工程應(yīng)力［8］:

其中λU、λB和λP分別為測試得到的單軸拉伸、等雙軸拉伸和平面剪切收縮率。三種不同的變形方式的應(yīng)變不變量(I1，I2)與主拉伸比之間的關(guān)系，可由式(2)求得。

采用最小二乘法擬合實測的應(yīng)力－應(yīng)變和由本構(gòu)模型計算得到的應(yīng)力－應(yīng)變，可以求得各本構(gòu)模型的材料常數(shù):

2 橡膠材料的應(yīng)力－應(yīng)變關(guān)系實驗測試

根據(jù)上節(jié)所述的不同本構(gòu)模型的材料常數(shù)獲取方法，在對橡膠材料進行應(yīng)力－應(yīng)變測試時，對橡膠試件進行單軸拉伸、等雙軸拉伸以及平面剪切實驗［11］，以獲得橡膠材料在不同應(yīng)力－應(yīng)變狀態(tài)和不同最大應(yīng)變下的應(yīng)力－應(yīng)變關(guān)系。橡膠材料的應(yīng)力－應(yīng)變實驗是委托美國Axel Products，Inc測試的［12］。測試用的橡膠材料的邵氏硬度為50。

圖2為在不同應(yīng)力－應(yīng)變狀態(tài)和不同最大應(yīng)變下，實測的工程應(yīng)力－工程應(yīng)變關(guān)系。在同一應(yīng)力－應(yīng)變狀態(tài)，不同的最大工程應(yīng)變時，工程應(yīng)力－工程應(yīng)變曲線是不重合的，在加載和卸載過程中的應(yīng)力－應(yīng)變曲線也是不重合的。這是由于橡膠材料存在Mullins效應(yīng)［13］的作用:即在對橡膠試件進行實驗過程中，在加載和卸載中應(yīng)力會出現(xiàn)變小的現(xiàn)象。在進行橡膠材料的應(yīng)力－應(yīng)變測試時，為保證計算精度以及減小Mullins效應(yīng)的影響，取循環(huán)加載(循環(huán)次數(shù)一般為5次)的最后一次應(yīng)力－應(yīng)變數(shù)據(jù)為實測的工程應(yīng)力與工程應(yīng)變的關(guān)系［14］。

圖1 橡膠材料的應(yīng)力－應(yīng)變測試圖Fig.1 The stress versus strain experimental diagram of rubber material

圖2 實測工程應(yīng)力－工程應(yīng)變關(guān)系圖Fig.2 The experimental engineering stress versus engineering strain relations3本構(gòu)模型材料常數(shù)

3 本構(gòu)模型材料常數(shù)

3.1 應(yīng)力－應(yīng)變組合

由于等雙軸拉伸等效單軸壓縮［1］，且本文計算的橡膠隔振器的總體變形的應(yīng)變水平小于0.5，故等雙軸拉伸應(yīng)變的最大值僅取1。為了研究應(yīng)變狀態(tài)和最大應(yīng)變對橡膠隔振器的力－位移計算結(jié)果的影響，選取其中不同應(yīng)力－應(yīng)變狀態(tài)和最大應(yīng)變進行組合，以得到同一本構(gòu)模型在不同應(yīng)力－應(yīng)變狀態(tài)和最大應(yīng)變下的本構(gòu)模型常數(shù)。不同應(yīng)力－應(yīng)變狀態(tài)和最大應(yīng)變組合如表1所示。

表1 應(yīng)力－應(yīng)變狀態(tài)和最大應(yīng)變組合Tab.1 Stress versus strain conditions and maximum strain combination

3.2 應(yīng)力－應(yīng)變擬合值和實測值的對比及材料常數(shù)

圖3所示為邵氏硬度為50的橡膠材料實測以及擬合的工程應(yīng)力－工程應(yīng)變曲線。選取1號、3號組合計算結(jié)果進行分析。

由圖3可見:

(1)在單軸拉伸應(yīng)力－應(yīng)變狀態(tài)下，1號、3號組合時，Mooney-Rivlin模型擬合效果最好，Marlow模型擬合測試曲線趨勢較好，Van Der Waals模型擬合效果相對較差，其他三種本構(gòu)模型擬合材料常數(shù)效果較好。

(2)在等雙軸拉伸應(yīng)力－應(yīng)變狀態(tài)下Mooney-Rivlin模型擬合效果最好，其次是Van Der Waals模型，Ogden模型和Yeoh模型相對于Mooney-Rivlin模型來說，擬合效果相對較差，擬合效果最差的是Marlow模型，特別在3號組合時，Marlow模型擬合曲線基本偏離實測應(yīng)力－應(yīng)變曲線。

(3)在平面剪切應(yīng)力－應(yīng)變狀態(tài)下，Mooney-Rivlin模型擬合誤差最小，其次是Arruda-Boyce模型，擬合效果相對較差的是Marlow模型。

綜合以上分析結(jié)可知，Mooney-Rivlin模型在單軸拉伸、等雙軸拉伸和平面剪切三種應(yīng)力－應(yīng)變狀態(tài)下的擬合效果最好，Arruda-Boyce模型和Van Der Waals模型次之，擬合效果最差的是Marlow模型。表2列出了在表1中的九種應(yīng)力－應(yīng)變狀態(tài)和最大應(yīng)變組合下，各本構(gòu)模型的材料常數(shù)。

4 本構(gòu)模型常數(shù)對橡膠隔振器力－位移關(guān)系計算結(jié)果的影響

本文選取兩種不同類型的懸置做為研究對象。懸置1主要承受拉壓或者剪切變形，而懸置2則同時承受拉壓和剪切變形。

圖3 1、3號組合下，實測與擬合的工程應(yīng)力－工程應(yīng)變曲線Fig.3 The experiment and estimated curves of engineering stress versus strain under combination 1，3

4.1 懸置1

圖4所示的橡膠懸置在X向只承受拉壓、Y向只承受剪切。在對此橡膠懸置進行有限元分析時，在X、Y、Z 三向均加載 5 mm的位移，計算所需的外力。

用Mooney－Rivlin本構(gòu)模型，利用表1中的1號、2號、3號、4號、5號和7號、8號、9號應(yīng)力－應(yīng)變狀態(tài)組合和不同的最大應(yīng)變狀態(tài)得到的本構(gòu)模型材料常數(shù)，對懸置1進行三向靜剛度的有限元計算。計算得到的懸置1在X、Y、Z三個方向的剛度和實測值的對比見表3～表5。

由表3可見，在懸置只承受拉壓變形的X向，利用1號、2號、3號組合擬合得到的本構(gòu)模型常數(shù)，計算得到的靜剛度，比只在單軸拉伸應(yīng)力－應(yīng)變狀態(tài)(7號)下計算的靜剛度相對誤差小。

由表4可見，在懸置只受剪切變形的Y向，利用三種應(yīng)力－應(yīng)變狀態(tài)擬合得到的材料常數(shù)，計算得到的靜剛度，比只在平面剪切應(yīng)力－應(yīng)變狀態(tài)(9號)下計算得到的靜剛度相對誤差大。

圖4 懸置1有限元模型Fig.4 FEA mode of mount I

由表5可見，在同時承受壓縮和剪切應(yīng)力的Z向，在單軸、等雙軸拉伸兩種應(yīng)力－應(yīng)變狀態(tài)下計算的靜剛度相對誤差最小，而在實測等雙軸拉伸應(yīng)力－應(yīng)變狀態(tài)(8號組合)下計算得到的靜剛度相對誤差最大。

4.2 懸置2

圖5所示的懸置2橡膠主簧同時承受拉壓變形和剪切變形。在進行力－位移關(guān)系計算時，X向加載力的范圍從 －3 500～3 500 N;Y向加載力的范圍從－500～500 N;Z向加載力的范圍從 －1 000～1 000 N。

(1)利用Mooney-Rivlin本構(gòu)模型，和不同實測應(yīng)力－應(yīng)變組合下擬合得到的本構(gòu)模型材料常數(shù)，計算其三向力－位移曲線，計算結(jié)果見圖6。在X向和Z向，Mooney-Rivlin模型在1號應(yīng)力－應(yīng)變組合下計算得到靜剛度值與實測值的相對誤差較小;在Y向，在3號應(yīng)力－應(yīng)變組合下計算的靜剛度與實測值的誤差較大。

(2)利用Arruda-Boyce模型、Van der Waals模型、Ogden模型和Yeoh模型，選取4組應(yīng)力－應(yīng)變狀態(tài)組合，對懸置2進行三向靜剛度計算。計算結(jié)果分別見表6、表7、表8和表9。由表6～表9可見，由三種應(yīng)力－應(yīng)變狀態(tài)下計算的靜剛度值，比用兩種應(yīng)力－應(yīng)變狀態(tài)下計算得到的靜剛度值更準(zhǔn)確。

圖5 懸置2有限元模型Fig.5 FEA mode of mount II

表2 各本構(gòu)模型的材料常數(shù)Tab.2 The constitute constants in different constitute models

表3 X向靜剛度計算值、實測值及其相對誤差Tab.3 The estimated and experimental stiffness and their relative error in direction X

表4 Y向靜剛度計算值、實測值及其相對誤差Tab.4 The estimated and experimental stiffness and their relative error in direction Y

表5 Z向靜剛度計算值、實測值及其相對誤差Tab.5 The estimated and experimental stiffness and their relative error in direction Z

表6 利用Arruda-Boyce模型計算得到的懸置2的三向靜剛度與實測值的相對誤差Tab.6 The estimated and experimental stiffness and their relative error of mount II in three vertical directions using Arruda-Boyce mode

表7 利用Van der Waals模型計算得到的懸置2的三向靜剛度與實測值的相對誤差Tab.7 The estimated and experimental stiffness and their relative error of mount II in three vertical directions using Van der Waals mode

表8 利用Ogden模型計算得到的懸置2的三向靜剛度與實測值的相對誤差Tab.8 The estimated and experimental stiffness and their relative error of mount II in three vertical directions using Ogden mode

表9 利用Yeoh模型計算得到的懸置2的三向靜剛度與實測值的相對誤差Tab.9 The estimated and experimental stiffness and their relative error of mount II in three vertical directions using Yeoh mode

5 本構(gòu)模型對力－位移特性計算結(jié)果的影響

在同種應(yīng)力－應(yīng)變狀態(tài)組合下，用不同本構(gòu)模型對懸置2進行靜剛度計算并與實測值對比。選取1號、3號應(yīng)力－應(yīng)變組合下的計算結(jié)果進行分析，如圖6所示。

由圖6可見，用Mooney-Rivlin模型計算的懸置在X、Y和Z三向靜剛度計算值與測試值最接近，其相對誤差小于10%。用Van der Waals模型和Marlow模型計算得到的靜剛度與實測值相對誤差較大;利用Arruda-Boyce計算得到的靜剛度與實測值吻合的較好，計算精度較高，大部分的靜剛度計算值與實測值的相對誤差小于10%。在3號應(yīng)力－應(yīng)變組合下，Ogden模型和Yeoh模型計算得到的靜剛度與實測值相對誤差小于10%。

由此可見，對于研究的幾種本構(gòu)模型，利用Mooney-Rivlin模型計算得到的懸置的靜剛度計算值與實測值的相對誤差最小。其次是Arruda-Boyce模型，誤差最大的是Marlow模型和Van der Waals模型。

6 結(jié)論

根據(jù)以上計算結(jié)果，可得到如下結(jié)論:

(1)對于只受拉壓變形的橡膠隔振器，用三種應(yīng)力－應(yīng)變狀態(tài)比只用單軸拉伸應(yīng)力－應(yīng)變狀態(tài)下擬合得到的本構(gòu)模型常數(shù)計算得到的靜剛度計算值精度高。對于只受剪切變形的橡膠隔振器，用平面剪切應(yīng)力－應(yīng)變狀態(tài)比用三種應(yīng)力－應(yīng)變狀態(tài)擬合得到的本構(gòu)模型材料常數(shù)計算得到的靜剛度精度高。

(2)對橡膠隔振器進行力－位移特性分析時，需選取合適的超彈性橡膠本構(gòu)模型。本文的計算結(jié)果表明Mooney－Rivlin模型的計算精度最高，其相對誤差均小于10%，其次是Arruda－Boyce模型，計算相對誤差最大的是Marlow模型。在實際工程應(yīng)用時，需要根據(jù)橡膠隔振器的結(jié)構(gòu)特點以及受載情況，選取最適合的本構(gòu)模型。

圖6 利用三種應(yīng)力－應(yīng)變應(yīng)變擬合得到的本構(gòu)模型材料常數(shù)計算得到的懸置2三向力－位移關(guān)系(圖中:Kt－測試剛度;Kmr－利用Mooney-Rivlin模型計算得到的剛度;Kv－利用Van der Waals模型計算得到的剛度;Ko－利用Ogden模型計算得到剛度;Ky－利用Yeoh模型計算得到的剛度;Ka－利用Arruda－Boyce模型計算得到的剛度;Kma－利用Marlow模型計算得到的剛度)Fig.6 The relation of force versus displacement for mount II calculated by the constitute constants obtained by the stress versus strain conditions(In this diagram，Kt stand for experiment stiffness，Kmr stand for estimated stiffness using Mooney-Rivlin model，Kv stand for estimated stiffness using Van der Waals model，Ko stand for estimated stiffness using Ogden model，Ky stand for estimated stiffness using Yeoh model，Ka stand for estimated stiffness using Arruda-Boyce model，Kma stand for estimated stiffness using Marlow model)

(3)利用不同的最大應(yīng)變和不同的應(yīng)力－應(yīng)變狀態(tài)組合得到的本構(gòu)模型常數(shù)，對橡膠隔振器力－位移特性分析有很大影響。選取更接近實際變形狀態(tài)的應(yīng)力－應(yīng)變狀態(tài)，所計算出的橡膠隔振器的力－位移特性更準(zhǔn)確。

［1］ Seibert D J，Schoche N.Direct comparison of some recent rubber elasticity models［J］. RubberChemistry and Technology，2000，73:366384.

［2］黃建龍，解廣娟，劉正偉.基于Mooney-Rivlin模型和Yeoh模型的超彈性橡膠材料有限元分析［J］.橡膠工業(yè)，2008，55(8):467－470.

［3］Wang L R，Lu Z H.Modeling methods of constitutive law of rubber hyperelasticity based on finite elemen simulations［J］.Rubber Chemistry and Technology，2003，76(1):271 －285.

［4］Gracia L A，Liarte E，Pelegay J L，et al.Finite element simulation of the hysteretic behavior of an industrial rubber［J］.Application to Design of Rubber Components.Finite Elements in Analysis and Design，2010，46(4):357 －368.

［5］ Hartmann S， Tschope T， SchreiberL， etal. Finite deformations of a carbon black-filled rubber.Experiment，optical measurement and material parameter identification using finite elements［J］.European Journal of Mechanics，2003，22(3):309 －324.

［6］郭仲衡.非線性彈性理論［J］.北京:科學(xué)出版社，1980.

［7］ Morman K N，Nagtegall J C.Finite element analysis of viscoelastic elastomeric structures vibrating about non-linear statically stresses configurations［C］.SAE 811309.

［8］Shangguan W B，Lu Z H，Shi J J.Finite element analysis of static elastic characteristics ofthe rubberisolators in automotive dynamic systems［R］.SAE Technical Paper Series，2003，2003－01－0240.

［9］ ABAQUS ，Inc.Abaqus Analysis User's Manual［M］.2006.

［10］ Yeoh O H.On the Ogden strain-energy function［J］.Rubber Chemistry and Technology，1997，70(2):175 －182.

［11］楊曉翔.非線性橡膠材料的有限元單元法［M］.北京:石油工業(yè)出版社，1999:4.

［12］ Olsson A K.Finite element procedures in modelling the dynamic properties of rubber［D］. Lund:Lund University，2007.

［13］ Mullins L.Softening of rubber by deformation［J］.Rubber Chemical Technology，1969，42:339 －362.

［14］ Kim W D，Kim D J，Kim W S，et al，A study on the equibiaxial tension test of rubber material［J］.Korean Society of Mechanical Engineering，2003，11(5):95 －104.

［15］崔忠圻.金屬學(xué)與熱處理［M］.北京:機械工業(yè)出版，2001.

［16］潘孝勇.橡膠隔振器動態(tài)特性計算與建模方法的研究［D］.杭州:浙江工業(yè)大學(xué)，2009.