馬 永，趙世平，盧丙舉，冷鼎鑫

(1.中國船舶集團有限公司第七一三研究所 河南 鄭州 450015；2.河南省水下智能裝備重點實驗室 河南 鄭州 450015；3.中國海洋大學(xué) 機械電子工程學(xué)院，山東 青島 266003)

0 引 言

橡膠減震器是發(fā)射裝置重要組成部分，布置在發(fā)射筒的周向環(huán)形空間，起到氣密、導(dǎo)向、適配、減震等作用。發(fā)射裝置減震器一般采用橡膠等材質(zhì)，工作過程需要承受長時間預(yù)壓，間歇式瞬態(tài)大能量沖擊載荷作用。橡膠除了具有明顯的超彈性外，還具有蠕變等時間、溫度相關(guān)的粘彈性力學(xué)特性。橡膠長時間壓縮下力學(xué)性能的變化，影響減震器長期貯存使用壽命問題，以及發(fā)射過程的安全性。

蠕變研究的是材料“長期使用性能”，即形變、應(yīng)力隨時間和溫度變化的現(xiàn)象與過程，是材料靜態(tài)粘彈性的基本表現(xiàn)。由于問題的復(fù)雜性和特殊性，工程一般多采用試驗方式分析橡膠減震器的蠕變性能[1-3]，但試驗往往有一定的離散性，且長時間蠕變試驗周期過長、成本大，很難得出規(guī)律性結(jié)論。在仿真應(yīng)用方面，目前常采用Kelvin模型、Burgers模型和標準線性固體模型等線性機械模型描述材料的蠕變行為，但這些本構(gòu)模型往往都有一定的局限性。近些年隨著非線性粘彈性材料的研究急劇增加，非線性粘彈性理論不斷得到完善，對非線性粘彈性本構(gòu)模型的研究取得了大量的進展。Hanyga[4]、高慶等[5]研究了材料的多重積分型本構(gòu)模型，但多重積分型本構(gòu)方程由于應(yīng)力、應(yīng)變關(guān)系式過于復(fù)雜，本構(gòu)模型參數(shù)不能從試驗中直接確定，難以在工程上推廣應(yīng)用；冪律關(guān)系型本構(gòu)關(guān)系函數(shù)數(shù)形式簡單、數(shù)值實現(xiàn)方便，在蠕變變形范圍小的工況中應(yīng)用廣泛[6]；Bemstein，Kearsley 和 Zapas[7-8]依據(jù)Green-Rivlin，Coleman和Noll理論，詳盡地研究了有限變形的應(yīng)力松弛問題，推導(dǎo)了不可壓縮固體和流體的非線性粘彈性本構(gòu)關(guān)系，建立了BKZ理論；Christensen將非線性橡膠彈性理論擴展到了非線性粘彈性的情況，提出了一種適用于不可壓縮橡膠類材料的非線性粘彈性理論。Luo[9-12]根據(jù)工程學(xué)原理，引入損傷函數(shù)，在超彈性本構(gòu)的基礎(chǔ)上，通過蠕變/松弛強度系數(shù)和蠕變/松弛軟化系數(shù)描述橡膠材料的蠕變特性。工程中實測獲得的工程對象隨時間的響應(yīng)（如載荷、變形位移等）隱式地蘊含著材料的本構(gòu)信息，即工程對象的應(yīng)力-應(yīng)變本構(gòu)關(guān)系和模型參數(shù)，學(xué)者們應(yīng)用不同系統(tǒng)識別模擬工具，針對本構(gòu)模型選擇（模型辨識）及其參數(shù)確定等問題進行研究[13]。

本文開展減震器材料靜態(tài)力學(xué)以及48 h變載壓縮蠕變試驗，構(gòu)建基于損傷函數(shù)的時變?nèi)渥儽緲?gòu)模型，利用abaqus/UHYPER用戶子程序接口對本構(gòu)模型進行二次開發(fā)，完成本構(gòu)模型參數(shù)識別。為更好地描述橡膠減震器蠕變損傷的力學(xué)特性，采用Von Mises屈服準則來評估橡膠材料的蠕變殘余變形，建立一種考慮永久變形的修正蠕變本構(gòu)模型。

1 試驗研究

1.1 靜態(tài)力學(xué)試驗

靜態(tài)力學(xué)試驗的目的是獲取橡膠材料在單軸拉伸、等雙軸拉伸及平面拉伸變形下的應(yīng)力-應(yīng)變關(guān)系，并擬合得到相應(yīng)的本構(gòu)模型參數(shù)。

試驗設(shè)備包括電子萬能試驗機（型號WDW-10）、激光引伸儀（型號LE-05）裁刀、試樣基板、游標卡尺等。試驗溫度保持在常溫恒溫（23 ℃ ± 2°環(huán)境溫度）。試驗件與試驗過程、應(yīng)力應(yīng)變曲線，如圖1所示。

1.2 時變?nèi)渥兞W(xué)試驗

圖1 橡膠材料的靜態(tài)力學(xué)性能試驗Fig.1 Static mechanical properties test of rubber materials

時變?nèi)渥兞W(xué)試驗參考國家標準《GB19242-2003硫化橡膠在壓縮或剪切狀態(tài)下蠕變的測定》進行，試驗設(shè)備與橡膠材料靜態(tài)力學(xué)性能試驗設(shè)備相同。具體試驗流程為：將橡膠材料試樣放置于實驗機夾具正中心，調(diào)整試驗機，使夾具剛好與橡膠材料試樣接觸并且處于不受力狀態(tài)，調(diào)節(jié)壓縮試樣變形至中等變形，再返回原位，重復(fù)幾次，完成機械調(diào)節(jié)。機械調(diào)節(jié)的全過程在標準室內(nèi)溫度下進行，機械調(diào)節(jié)在試驗之間間隔最長不超過48 h。機械調(diào)節(jié)結(jié)束后，開始正式試驗，在不超過15 s的時間內(nèi)施加全力于試樣上，并且不明顯過量；達到滿負荷力之后，在不同的時間測定試樣的變形，這個力在整個試驗過程中保持恒定，力保載即48 h。橡膠材料的壓縮蠕變試驗外載荷分別取為1 000 N，1 500 N。對于橡膠材料，僅進行常溫壓縮蠕變試驗。在試驗過程中，記錄橡膠材料在不同加載條件下的蠕變試驗數(shù)據(jù)（時間、變形及載荷）。

橡膠試樣在不同加載水平下的蠕變試驗結(jié)果，如圖2所示。

由圖2可知，隨著時間的增加，橡膠材料的蠕變曲線斜率逐漸減小，表明橡膠的蠕變變形速率隨著時間的增加而逐漸變慢。

圖2 橡膠試件不同加載水平下蠕變試驗結(jié)果Fig.2 Creep test results of rubber specimens under different loading levels

2 材料本構(gòu)模型及參數(shù)辨識

2.1 超彈性本構(gòu)模型

橡膠材料的靜態(tài)力學(xué)特征可采用基于應(yīng)變能函數(shù)的超彈本構(gòu)模型來描述?；趹?yīng)變不變量表示的連續(xù)介質(zhì)唯象模型可寫成：

式中，F(xiàn)為變形梯度張量，計算如下式：

根據(jù)式（1）可推導(dǎo)出不同的超彈本構(gòu)模型，其中多項式模型被廣泛應(yīng)用于預(yù)測橡膠材料在外荷載激勵下的力學(xué)響應(yīng)。多項式模型可寫為：

式中：Cij，Dij為橡膠材料的超彈系數(shù)，可通過橡膠的超彈力學(xué)試驗來確定。

對于式（6）取前兩項，則可得目前有限元計算中應(yīng)用最廣泛的Mooney-Rivlin模型，其表達式為：

2.2 唯象蠕變本構(gòu)模型

超彈本構(gòu)模型可用于分析橡膠材料的靜力學(xué)特性，但其與時間效應(yīng)無關(guān)，無法描述橡膠的粘彈性特性。為了評估橡膠材料在常量載荷作用下的蠕變效應(yīng)，需對傳統(tǒng)超彈本構(gòu)模型進行修正。引入帶有損傷意義的蠕變函數(shù)Wc-r，將橡膠材料的蠕變變形與時間相連，并假設(shè)蠕變變形與應(yīng)變不變量，及時間歷程有關(guān)。由此，傳統(tǒng)的超彈性模型與時間歷程建立了聯(lián)系。修正后的唯象蠕變本構(gòu)模型為：

為描述時變損傷效應(yīng)，構(gòu)建冪函數(shù)形式的蠕變損傷函數(shù)，其具體表達式為[10-13]：

從數(shù)學(xué)角度考慮，式（8）為連續(xù)函數(shù)；從工程學(xué)方面考慮，式（9）提出的蠕變函數(shù)其參數(shù)均具有物理含義。其中T為荷載加載結(jié)束時刻的時間，k表征蠕變的損傷幅值，代表蠕變強度，r為衰減系數(shù)，描述了損傷隨時間衰減的程度。

2.3 本構(gòu)模型參數(shù)識別

從本質(zhì)上來說，式（8）提出的唯象蠕變本構(gòu)模型是修正的超彈性本構(gòu)模型，因此超彈性本構(gòu)模型參數(shù)識別的準確度極為重要。根據(jù)最小二乘原理，利用橡膠材料的單軸拉伸、等雙軸拉伸、平面拉伸試驗數(shù)據(jù)，擬合得到本構(gòu)模型參數(shù)。文獻[14]提出采用擬合優(yōu)度R2來評估擬合誤差，其表達式為：

表1 橡膠材料超彈性本構(gòu)模型參數(shù)Tab.1 Parameters of hyperelastic constitutive model for rubber materials

為進一步驗證模型參數(shù)的準確性，通過Abaqus有限元軟件，對橡膠材料的靜態(tài)力學(xué)特性進行數(shù)值模擬。根據(jù)試樣尺寸，建立橡膠試件的軸對稱幾何模型（Φ29 mm×12.5 mm），模型的上端與解析剛體相連，并在橡膠與剛體之間進行自由度耦合設(shè)置。采用CAX4HT單元對模型進行網(wǎng)格劃分，有限元模型如圖3所示。

圖3 橡膠試件有限元模型Fig.3 Finite element model of rubber specimen

將數(shù)值模擬的靜態(tài)壓縮特性與試驗結(jié)果進行對比，如圖4所示。

圖4 橡膠材料靜態(tài)壓縮試驗與仿真結(jié)果對比Fig.4 Comparison of static compression test and simulation results of rubber materials

由圖4可知，試驗與仿真結(jié)果吻合度較好，這說明超彈本構(gòu)模型及參數(shù)選取合理，能夠表征橡膠在加載階段的非線性彈性特性，數(shù)值模擬精度良好。

基于式（10），利用Abaqus用戶自定義子程序UHYPER，將提出的蠕變本構(gòu)模型通過Fortran語言進行編譯，并將編譯的子程序與有限元軟件進行關(guān)聯(lián)，開展橡膠材料的蠕變數(shù)值模擬。參照壓縮蠕變的試驗過程，仿真分析2種加載水平（1.0 kN和1.5 kN）下48 h的時變?nèi)渥冃袨椤?/p>

為明確蠕變本構(gòu)參數(shù)對蠕變特性的影響程度，開展參數(shù)敏感性分析，敏感性范圍定為±20%。具體分析過程為：1）選擇一組蠕變參數(shù)作為參考參數(shù)，這里選取常溫1.5 kN作用下的橡膠蠕變參數(shù)（k0，r0）作為參考值；2）改變一個蠕變參數(shù)的數(shù)值并保持另一個蠕變參數(shù)數(shù)值不變，分析時變?nèi)渥冏冃蔚淖儎臃秶?，如圖5所示。

不斷調(diào)整蠕變參數(shù)k和r，使仿真結(jié)果與試驗得出時間-蠕變變形誤差變?。环謩e對不同加載水平（1 000 N 和 1 500 N）的時間-蠕變變形試驗數(shù)據(jù)進行擬合，最終得出蠕變參數(shù)擬合結(jié)果（見表2）。

圖5 冪指數(shù)蠕變本構(gòu)模型參數(shù)敏感性分析Fig.5 Parameter sensitivity analysis of power exponent creep constitutive model

表2 橡膠蠕變參數(shù)擬合結(jié)果（23 ℃）Tab.2 Fitting results of rubber creep parameters（23 ℃）

為考察時變?nèi)渥儽緲?gòu)模型及參數(shù)的模擬精度，與橡膠試件不同載荷下的時間-蠕變變形進行對比。結(jié)果表明，試驗與仿真預(yù)測結(jié)果吻合度較好（見圖6）。

圖6 橡膠試件時間-變形仿真與試驗結(jié)果對比Fig.6 Comparison of time deformation simulation and experimental results of rubber specimen

為了詳細地說明時變?nèi)渥冏冃我?guī)律，圖7對比了仿真得出的不同時刻橡膠材料的軸向蠕變變形。在蠕變過程中，橡膠材料的變形隨時間增加呈現(xiàn)非線性的變化特性。

圖7 橡膠試件軸向蠕變變形分布Fig.7 Distribution of axial creep deformation of rubber specimens

由圖7可見，橡膠的軸向變形云圖呈現(xiàn)分層的現(xiàn)象：在材料底部，由于受到邊界條件的限制，變形量近似為0；在材料頂部，由于載荷條件的影響，變形量最大；而在其他區(qū)域，變形是逐漸變化的。這種蠕變變形的分層現(xiàn)象并未因蠕變時間的改變而發(fā)生變化。此外，材料的最大蠕變變形隨著時間的增加而增大，材料自由表面的膨脹效應(yīng)也隨著時間的推移而愈發(fā)明顯。

在超彈性本構(gòu)模型的基礎(chǔ)上引入蠕變損傷函數(shù)，構(gòu)建的時變?nèi)渥儽緲?gòu)模型，很好描述了橡膠材料的壓縮蠕變性能。

3 減震器蠕變性能分析

圖8 橡膠減震器有限元模型Fig.8 Finite element model of rubber shock absorber

建立橡膠減震器的有限元模型，如圖8所示。模型上端與解析剛體相接觸，下端與剛性體連接，采用具有6個自由度的C3D8HT單元對模型進行網(wǎng)格劃分。載荷沿橡膠減震器垂直方向施加于上端解析剛體上，并通過剛體與橡膠器件之間的自由度耦合傳遞至器件。為了防止解析剛體與橡膠材料相互滲透，對接觸區(qū)域進行面與面接觸設(shè)置。

在橡膠減震器上表面加載2.1 t的載荷，橡膠減震器在在室溫條件下48 h蠕變過程中的應(yīng)力分布云圖（見圖9）。

橡膠減震器在蠕變過程中的位移分布云圖，如圖10所示。

圖9 橡膠減震器在不同時刻的應(yīng)力分布云圖Fig.9 Stress distribution nephogram of rubber shock absorber at different times

圖10 橡膠減震器在不同時刻的位移分布云圖Fig.10 Displacement distribution nephogram of rubber shock absorber at different times

由圖9和圖10可知，在蠕變過程中橡膠減震器的應(yīng)力最大值始終出現(xiàn)在器件的外側(cè)位置；橡膠器件與材料的分布云圖相似，均出現(xiàn)了分層現(xiàn)象，但器件的分層變形不如材料的明顯，這是由于器件結(jié)構(gòu)效應(yīng)引起的。隨著蠕變時間的增加，減震器的蠕變變形也在不斷增加，橡膠減震器蠕變行為的時間-蠕變變形曲線如圖11所示。

圖11 橡膠減震器時間-蠕變變形曲線Fig.11 Time creep deformation curve of rubber shock absorber

由圖11可知，橡膠減震器在15 s時的蠕變變形為5.43 mm，48 h 后蠕變變形達到 5.97 mm，在整個時變過程中，蠕變變形較為平緩。

為了描述橡膠減震器在蠕變條件下的損傷永久變形，我們對時變?nèi)渥儽緲?gòu)模型進行進一步擴展，采用Von Mises屈服準則來評估橡膠材料的蠕變殘余變形，提出一種考慮永久變形的修正蠕變本構(gòu)模型，其表達式為：

式中，基于屈服準則的分段本構(gòu)方程為：

式中： σy0， σy1， σy2為第1，2，3階屈服應(yīng)力； εy0，εy1， εy2為各屈服應(yīng)力對應(yīng)的應(yīng)變；為硬化系數(shù)；KP是長期模量。

橡膠器件蠕變前后的剛度變化，如圖12所示。

由圖12可知，當位移為5 mm時，減震器未蠕變前的承載為17 045.4 N，引入蠕變損傷后的器件承載力為13 430.2 N，降低了21.2%。定義等效剛度為：

圖12 橡膠器件蠕變前后剛度對比Fig.12 Comparison of stiffness of rubber components before and after creep

其中， ΔF為載荷變化量， Δx為位移變化量。由圖12可知，橡膠器件在蠕變前與蠕變后的等效剛度分別為3.89×106N/m 和 3.53E6 N/m，退化了 9.1%。

橡膠減震器蠕變前與損傷蠕變后靜態(tài)變形不同位移處的應(yīng)力分布云圖，如圖13所示。

圖13 蠕變前后不同位移處應(yīng)力分布云圖Fig.13 Nephogram of stress distribution at different displacements before and after creep

由圖13可知，減震器在考慮蠕變損傷效應(yīng)后，其靜態(tài)變形的最大應(yīng)力呈現(xiàn)退化趨勢。當靜態(tài)位移為1 mm時，蠕變前與蠕變后的應(yīng)力最大值分別為0.41 MPa和0.17 MPa；當位移為3 mm時，蠕變前/后的應(yīng)力最大值分別為 1.10 MPa 和 0.67 MPa；當位移為 5 mm 時，蠕變前/后的應(yīng)力最大值分別為 1.80 MPa 和 1.33 MPa。

4 結(jié) 語

本文開展了減震器橡膠材料在 1 000 N 和 1 500 N兩種載荷水平下48 h蠕變試驗研究，通過對蠕變試驗數(shù)據(jù)分析，可以得出以下結(jié)論：

1）在超彈性本構(gòu)模型的基礎(chǔ)上，引入損傷意義的蠕變函數(shù)，將荷載引起的橡膠材料的蠕變與時間相連，構(gòu)建了時變?nèi)渥儽緲?gòu)模型；采用Von Mises屈服準則來評估橡膠材料的蠕變殘余變形，提出一種考慮永久變形的修正蠕變本構(gòu)模型。

2）通過Abaqus/UPYPER用戶子程序，將橡膠材料的時變?nèi)渥儽緲?gòu)模型用Fortran語言編譯，構(gòu)建了橡膠材料及減震器蠕變性能仿真分析模型。

3）對不同載荷水平下的橡膠材料的蠕變試驗曲線對比，仿真與實驗結(jié)果吻合度較好，說明時變?nèi)渥儽緲?gòu)模型很好的描述了橡膠材料的蠕變性能，且預(yù)測精度滿足工程應(yīng)用。

4）通過橡膠減震器48 h壓縮蠕變試驗與仿真結(jié)果對比，驗證了方法的可行性。