李定機，毛成立，鄭 慶，童 悅，馬 超

(1.上海航天動力技術研究所，上海 201109；2.上海市塑料研究所有限公司，上海 201702)

0 引言

柔性噴管是固體火箭發(fā)動機推力矢量控制的結構形式之一，柔性接頭作為其關鍵部件，對柔性噴管的擺動特性有重要影響。柔性接頭一般由金屬增強件、橡膠彈性件、前后法蘭構成，前后法蘭中間包夾著數(shù)層增強件與數(shù)層彈性件交替層疊粘接的結構，從而形成可全軸擺動結構。橡膠材料在不同形式的載荷下表現(xiàn)許多與金屬材料不同的力學特性，如Payne效應、Mullins效應、超彈性以及粘彈性的應力松弛、蠕變、遲滯等。導彈飛行過程中，在燃氣高壓、伺服機構驅(qū)動的聯(lián)合動態(tài)載荷作用下，柔性接頭結構內(nèi)部受力情況較為復雜，對擺動力矩特性具有決定性影響的彈性件橡膠會表現(xiàn)出復雜的幾何、材料雙重非線性。因此，對彈性件材料力學行為的研究是柔性噴管靜態(tài)、動態(tài)力矩特性計算的重要基礎。

SHANI等在計算其設計的柔性接頭擺動力學性能時，在有限元程序中直接采用線彈性的假設對彈性件材料進行建模，彈性件在計算中的剪切模量被設為恒定值。王成林等及史宏斌等采用彈簧單元模擬柔性接頭彈性件的力學行為，彈簧單元的剛度系數(shù)矩陣來自噴管擺動試驗的擺動剛度。文獻[4-12]從彈性件材料本構出發(fā)，通過材料的基礎力學試驗確定了基于唯象理論的超彈性本構模型來描述彈性件的非線性彈性行為。目前，相關文獻資料中對柔性噴管彈性件建模都是基于線彈性或非線性超彈性的理想模型，適用于描述靜態(tài)或準靜態(tài)過程中的彈性行為，由于未考慮材料的粘彈性性質(zhì)，因此無法描述動態(tài)過程中與時間相關的各類力學行為。

本文針對柔性噴管用彈性件橡膠材料，建立了橡膠材料的超-粘彈性模型，進行不同加載速率和不同應變水平的試驗研究，并通過試驗驗證該模型和方法在模擬彈性件材料動態(tài)遲滯行為上的有效性，從而為后續(xù)柔性接頭動態(tài)特性的計算研究提供材料層面的建模方法。

1 超-粘彈性本構模型的建立

1.1 超彈性模型

橡膠通常被認為是不可壓縮的各向同性材料,其應力-應變的本構關系由含應變能函數(shù)W偏導項的等式描述，如式(1)所示。

(1)

式中為主應力；為應變能函數(shù)；為主伸長比；，分別為第一應變、第二應變不變量；為靜水壓強。

目前，有大量關于超彈性本構方程的研究，按機理可分為兩大類，即基于連續(xù)介質(zhì)力學的唯象理論和基于分子結構和構象熵改變的統(tǒng)計理論。在數(shù)值模擬中，常用的模型有Ogden、多項式模型及縮減多項式模型。在不可壓縮的假設下，一些常用的唯象理論模型的應變能函數(shù)如下：

Mooney-Rivlin：

=(-3)+(-3)

(2)

Neo Hooke：

=(-3)

(3)

Yeoh：

(4)

Ogden：

(5)

對于材料的單軸拉伸(UT)、平面拉伸(PT)和雙軸拉伸(ET)三種變形狀態(tài)而言，其Cauchy應力可通過應變能的偏導數(shù)形式表示成式(6)～式(8)，從而能通過各類試驗來確定應變能表達式中的各個材料系數(shù)。

(6)

(7)

(8)

1.2 超-粘彈性模型

將粘彈性引入本構模型中，考慮時間相關性從而計算材料動態(tài)的特性。常用的超-粘彈性模型的建立方式是通過粘壺和彈簧單元的串聯(lián)和并聯(lián)構造一個力學單元，以模擬材料的彈性和粘性。BERGSTR?M等從高分子鏈的蠕動運動出發(fā)，推導建立了非線性超-粘彈性模型。模型網(wǎng)絡如圖1所示，由平行的A 網(wǎng)絡和B網(wǎng)絡部分構成，A網(wǎng)絡能夠描述力學單元平衡狀態(tài)下的響應，而B網(wǎng)絡則描述時間影響下偏離平衡態(tài)的響應。兩個網(wǎng)絡的應變梯度張量相同，即==，二者貢獻的應力之和即為模型的瞬時應力。

圖1 Bergstr?m-Boyce網(wǎng)絡模型示意圖Fig.1 Schematic of Bergstr?m-Boyce network model

基于超彈性應變能函數(shù)理論，網(wǎng)絡中的彈性部分的本構關系形式如式(1)所示，描述與時間無關的應力-應變關系。

而粘性效應由式(9)所示的等效應變率描述：

(9)

另外，還需要定義應力比例因子，其意義是在瞬時載荷作用下，B網(wǎng)絡和A網(wǎng)絡承載的應力大小之比。

2 彈性件材料力學試驗

對于彈性件使用的橡膠材料來說，需要通過材料基礎力學試驗確定所使用的超彈性模型及其參數(shù)，這些試驗包括單軸拉伸和壓縮、雙軸拉伸和壓縮以及平面拉伸和壓縮試驗。考慮到柔性接頭的變形形式主要是剪切變形，并在擬合過程中，除了縮減多項式模型，Ogden模型和Mooney-Rivlin模型均需要兩種及以上的變形類型試驗數(shù)據(jù)，以得到確定的材料參數(shù)。因此，本文采用單軸拉伸試驗和四板剪切試驗研究彈性件材料的力學性能。

單軸拉伸試驗方法參照GB/T 528—2009，選擇總長度75 mm、厚度2 mm、中間寬度4 mm的Ⅱ型啞鈴型試樣；四板剪切試驗方法得到參照GB/T 12830—2008。兩種試驗均在電子拉力機上完成，如圖2所示。

(a)Uniaxial tensile test (b)4-slice shear test圖2 單軸拉伸和四板剪切試樣Fig.2 Uniaxial tensile specimen and 4-slice shear specimen

為了研究Mullins效應對該橡膠材料的影響，需要對試件進行數(shù)次加載循環(huán)，對單軸拉伸試件和四板剪切試件分別進行4次加載循環(huán)得到圖3和圖4的應力-應變曲線。由多次加載的曲線可以看出，對于兩種不同的變形形式來說，材料第一次加載和第二次加載的應力-應變曲線，都有顯著差異。在相同的應變下，第二次加載的應力相較于第一次加載的應力來說顯著減小，而在第二次加載之后，應力-應變曲線則逐漸穩(wěn)定下來。其中，單軸拉伸的第一次加載和第二次加載曲線之間差異較大，第二次加載應力值相比第一次最大減少了50.6%；拉伸應變在0%～150%的范圍內(nèi)，應力穩(wěn)定較快，在該應變區(qū)間內(nèi)可以選擇第三次或第四次的數(shù)據(jù)進行擬合，但應變超過150%時則需要更多次數(shù)的試驗來確定更穩(wěn)定的應力曲線；相比而言，剪切變形的第一次和第二次加載之間的差異相對更小，第二次加載應力值相比第一次最大減少了27.1%，并且應力值在更少的加載次數(shù)下能達到更穩(wěn)定的狀態(tài)。因此，剪切試驗進行3次加載循環(huán)就可確定應力軟化后的曲線。

圖3 單軸拉伸機械調(diào)節(jié)Fig.3 Mechanical adjustment on uniaxial tensile specimen

圖4 四板剪切機械調(diào)節(jié)Fig.4 Mechanical adjustment on 4-slice shear specimen

由于柔性接頭在工作過程中需要連續(xù)、多次擺動，因此應當采用機械調(diào)節(jié)之后的應力-應變數(shù)據(jù)擬合超彈性模型，即對試樣進行4次重復加載，并記錄每次應力-應變，取最后一次加載時的應力-應變，以消除Mullins效應。分別取圖3和圖4中最后一次加載曲線，并取二者同應變范圍內(nèi)數(shù)據(jù)以提高擬合準確度，如圖5所示，采用此數(shù)據(jù)用于擬合本構模型。Abaqus中超彈性參數(shù)的擬合需要平面拉伸的試驗數(shù)據(jù)來描述材料的剪切變形性能，而本文采用的是四板剪切的試驗。二者的區(qū)別在于前者是純剪切的變形形式，后者則是簡單剪切的變形形式，簡單剪切可以看作是一種純剪切加上轉動的變形形式。因此，在擬合之前要對數(shù)據(jù)進行轉換。參考文獻[14]可以推導出簡單剪切和純剪切的應力-應變關系如下：

(10)

式中，分別為四板剪切獲得剪切應變和剪切應力；、分別為平面拉伸的名義應變和名義應力。

圖5 超彈性模型擬合所用數(shù)據(jù)Fig.5 Data for fitting hyperelastic model

柔性噴管在不同的擺角幅度和擺動速率下，表現(xiàn)出的力矩特性也不同。為了研究柔性噴管往復擺動過程中的動態(tài)力矩特性與彈性件材料動態(tài)力學特性的關聯(lián)性，本文基于四板剪切試驗形式，在不同速率和不同最大應變水平的試驗條件下進行了剪切加載回復試驗。按照本文所研究柔性接頭的尺寸估算彈性件在在擺動到較大角度6°和中等角度3°時的最大剪切應變分別對應四板剪切試驗70%和150%的應變水平；而估算柔性接頭在中等擺角3°時，在1、0.5、0.1 Hz的擺動頻率下，剪切應變速率分別對應四板剪切試驗加載速率為420、240、50 mm/min時的剪切應變速率。因此，從上述估算得到的2個應變水平和3個速率水平中選擇并組合，形成四板剪切動態(tài)試驗的加載工況，圖6是四板剪切試件以240 mm/min的加載速率分別加載到應變?yōu)?50%和70%時，再以相同速度回復的遲滯曲線，圖7是四板剪切試件分別以420、240、50 mm/min的加載速率加載到應變150%附近的遲滯曲線。

圖6 不同應變水平下四板剪切遲滯曲線Fig.6 Shear hysteresis curve at differentstrain levels

圖7 不同加載速率下的四板剪切遲滯曲線Fig.7 Shear hysteresis curve at differentrates of loading

遲滯曲線圍成的面積表示材料在一個加載循環(huán)周期內(nèi)的能量損耗，這種損耗和材料的粘彈性密切相關。為了量化能量損耗大小，將遲滯環(huán)包絡面積與加載段曲線和應變軸所圍面積之比值定義為一個循環(huán)加載過程中的能量損耗率。從圖6中70%和150%應變幅度的遲滯曲線面積可以看出，應變幅度對能量損耗的影響很大，隨著最大應變的增加，循環(huán)過程的能量損耗增大，這主要是由于阻尼力做功的行程增加了，70%和150%應變條件下的能量損耗率分別為21.59%和17.59%。圖7反映了應變速率對遲滯曲線的影響，隨著加載速率的增加，加載在試件上的載荷也隨之增加。而隨著變形的增大，加載速率對載荷的影響也越來越顯著，在應變?yōu)?50%時，加載速率為240 mm/min和50 mm/min的曲線應力值分別相比于加載速率為420 mm/min的曲線對應應力值下降了4.24%和6.03%。而在所研究的加載速率范圍內(nèi)，能量損耗的變化不大，三個速率下能量損耗率為21.46%，21.59%和21.44%。這些應變幅度和應變速率對材料力學性能的影響規(guī)律與文獻[7]中所描述的擺動幅度和擺動頻率對柔性接頭力矩特性的影響規(guī)律是一致的。因此，針對彈性件材料本構進行建模,通過數(shù)值計算的方法可預示柔性接頭的靜態(tài)、動態(tài)力矩特性。

3 模型參數(shù)確定

3.1 超彈性模型參數(shù)的確定

根據(jù)圖5中的單軸拉伸和四板剪切的試驗數(shù)據(jù)，分別使用2階Ogden、3階Ogden、Mooney-Rivlin、Neo Hooke和Yeoh模型進行擬合，得到的原始數(shù)據(jù)與擬合曲線的對比圖如圖8和圖9所示。

圖8 單軸拉伸擬合結果Fig.8 Curves fitting to uniaxial tensile test

圖9 剪切擬合結果Fig.9 Curves fitting to simple shear test

對于單軸拉伸，Ogden模型相比多項式及縮減多項式模型更為逼近試驗數(shù)據(jù)；而對于四板剪切，列出的模型擬合效果較好，其中最大的平均誤差不超過20%。為了比較超彈性本構模型擬合單軸拉伸應力-應變和剪切應力-應變的整體效果，利用組試驗-預測數(shù)據(jù)可以計算擬合的歸一化殘差，同時考慮模型對單軸拉伸和四板剪切的整體預測效果，則整體殘差如式(11)所示。

(11)

采用式(11)計算得到表1中各個模型的整體殘差，由此判斷采用3階Ogden模型時整體擬合效果最好，相應的材料參數(shù)在表2中列出。根據(jù)Drucker穩(wěn)定性條件可知，這組本構參數(shù)在所研究的應變范圍內(nèi)能夠保證材料計算的穩(wěn)定。

表1 超彈性模型預測單軸拉伸和剪切的整體殘差

表2 三階Ogden模型材料參數(shù)

3.2 粘彈性模型參數(shù)的確定

粘性部分的材料參數(shù)具有相互耦合的特點，因此文獻[16-17]中給出了該模型下一般彈性體的粘性部分參數(shù)，即=5，=-1，=4，并建議從給出這組一般性參數(shù)出發(fā)，采取試錯法改變參數(shù)，以使得計算結果符合試驗曲線。為節(jié)省時間和成本，本文利用Python對Abaqus進行了二次開發(fā)，實現(xiàn)了參數(shù)化建模和計算結果的自動化處理。采用模擬退火算法，將材料參數(shù)和作為設計變量，將計算得到的應力-應變曲線和試驗的應力-應變曲線的平方誤差作為目標函數(shù)，編寫了相應的腳本，迭代優(yōu)化后得到如表3的材料參數(shù)。

表3 粘彈性模型材料參數(shù)

4 計算與結果分析

4.1 有限元模型的建立

在Abaqus中建立與試驗中所使用試件同尺寸的有限元模型，如圖10和圖11所示，分別為單軸拉伸試件和四板剪切試件的網(wǎng)格劃分以及邊界條件加載情況。橡膠部分均應用C3D8H雜交單元以適應橡膠的不可壓縮性，剛性板則采用C3D8R單元。另外，在四片剛性板中央分別定義參考點，基于這些參考點定義剛性板的剛體約束以節(jié)約計算資源。分析采用Dynamic，Implicit過程求解，以計算與時間相關的力學行為，利用表2和表3的參數(shù)定義橡膠的材料屬性。

圖10 單軸拉伸試件計算網(wǎng)格Fig.10 Computation grids on uniaxialtensile specimen

4.2 模型計算檢驗

采用第2節(jié)中建立的超-粘彈性模型進行計算，超彈性部分采用三階Ogden模型，粘性部分采用Bergstrom-Boyce遲滯模型。試件達到最大位移時的應力云圖如圖12和圖13所示，應力的分布符合實際標準件的變形規(guī)律。由于柔性噴管擺動時彈性件主要應變形式為剪切應變，故下文只針對剪切形式變形進行計算仿真。分別在加載速率為420、240、50 mm/min的條件下進行最大應變水平為100%的四板剪切試驗和對應的仿真計算，得到如圖14所示的應力-應變曲線。試驗和計算結果呈現(xiàn)出一致的規(guī)律，即隨著加載速率減小，應力值也減小，尤其是在應變較大的區(qū)域，加載速率對剪切應力的影響更為明顯，而應變較小時加載速率對剪切應力幾乎沒有影響。三個加載速率下的計算結果與試驗的平均誤差均不超過8%。

圖11 四板剪切試件計算網(wǎng)格Fig.11 Computation grids on 4-slice shear specimen

圖12 單軸拉伸應力云圖Fig.12 Stress contour of uniaxial tensile simulation

圖13 四板剪切應力云圖Fig.13 Stress contour of 4-slice shear simulation

(a)420 mm/min (b)240 mm/min (c)50 mm/min圖14 不同剪切速率下的剪切試驗與仿真對比Fig.14 Comparison between test and simulation of 4-slice shear at different rate of loading

在四板剪切有限元試件上施加速率為240 mm/min位移邊界條件，在位移達到5.6 mm，即70%應變時再以相同速度回復到零位移處，以模擬四板剪切加載回復循環(huán)的試驗。模擬計算的應力-應變數(shù)據(jù)與試驗數(shù)據(jù)對比如圖15所示，計算結果與試驗數(shù)據(jù)吻合較好，遲滯曲線形狀一致，曲線包圍面積相對誤差21.3%。由于該模型沒有考慮材料塑性，無法反映永久形變，導致曲線在原點附近的小應變范圍內(nèi)相對誤差較大，因此除去應變0.05以下的數(shù)據(jù)點，模型計算結果的平均相對誤差為9.80%。說明所使用的超-粘彈性模型能夠較好地計算材料的遲滯特性，從而描述材料力學行為的時間相關性。

因此，數(shù)學模型在超彈性和粘彈性力學行為的模擬上效果較好，能夠較準確地計算出具有非線性性、速率相關性和遲滯特性的應力-應變曲線，結果符試驗規(guī)律，可以進一步用于其他結構件的動態(tài)特性的計算。

圖15 四板剪切循環(huán)加載試驗與仿真結果Fig.15 Result of test and simulation in 4-slicecycle loading process

4.3 數(shù)值研究

4.3.1 加載條件對計算結果影響

針對不同加載速率和不同應變水平的情況進行剪切循環(huán)的計算。

圖16為70%應變水平下以速率240 mm/min和420 mm/min循環(huán)加載的計算結果對比，圖17為以240 mm/min速率分別加載到應變?yōu)?0%和110%的循環(huán)加載計算結果對比。

圖16 應變水平70%不同加載速率的循環(huán)加載計算結果Fig.16 Result of cycle loading simulation with 70%strain at different rates of loading

最大應變水平一定的情況下，隨著加載速度的增加，遲滯曲線的面積略有減小，加載到相同應變處所對應的應力略有增大；加載速率一定的情況下，最大應變水平遲滯環(huán)面積影響較大，即應變水平增大，循環(huán)過程能量損耗顯著增加，主要是由于回復過程中載荷力減小造成的。因此，模型能夠反映橡膠的動態(tài)遲滯行為及其在不同應變速率和不同應變水平下的變化規(guī)律。

圖17 加載速率240 mm/min時不同應變水平的循環(huán)加載計算結果Fig.17 Result of cycle loading simulation at 240 mm/minloading rate with different strain levels

4.3.2 材料參數(shù)對模型的影響

Bergstr?m-Boyce遲滯模型的材料參數(shù)需要通過不同速率下的試驗數(shù)據(jù)進行辨識，但是由于這幾個材料參數(shù)之間的耦合特性，確定所研究橡膠的粘彈性參數(shù)比較困難。參數(shù)和參數(shù)與高聚物分子的類型有關，與填充物無關；而參數(shù)和參數(shù)對模型影響較大。為了研究等效彈性伸長率和蠕變參數(shù)對模型計算遲滯特性的影響，以表3中參數(shù)為基礎，分別改變和進行加載速率240 mm/min，應變70%的加載回復計算，結果如圖18和圖19所示。

圖18 改變參數(shù)A的遲滯計算結果Fig.18 Calculation results of hysteresis processunder changing parameter A

結果顯示，增大參數(shù)，則遲滯曲線的包圍面積增大，并且在應變較大的區(qū)域，更大的值使得應力值更小。這是因為在式(9)中，參數(shù)以比例形式改變等效蠕變速率，在越大的值下計算的等效蠕變速率也越大，分子鏈來不及松弛便發(fā)生了滑移使得滯后效應增強。相比之下，參數(shù)對遲滯計算的影響更為顯著，參數(shù)增大，模型的遲滯特性也隨之增強。參數(shù)是模型網(wǎng)絡粘性部分(Network B)應力和彈性網(wǎng)絡部分(Network A)應力的比例，因此增大時粘性部分的耗散也更大，當減小到圖19中0.1的情況時，幾乎體現(xiàn)不出遲滯效應，所采用的超-粘彈性模型便退化為只有A網(wǎng)絡部分的超彈性模型。

圖19 改變參數(shù)S的遲滯計算結果Fig.19 Calculation results of hysteresis processunder changing parameter S

4 結論

(1)Mullins效應對四板剪切應力-應變曲線的影響比對單軸拉伸應力-應變曲線的影響更小，前者在第一次和第二次加載中應力值最大減少了50.6%，而后者應力值最大僅減少27.1%，另外四板剪切應力-應變曲線穩(wěn)定較快，僅在兩次加載之后的曲線就有較好的一致性。

(2)加載速率和應變水平對材料的遲滯特性都有影響；加載速率越大，應力值越大；應變水平增大，遲滯曲線面積顯著增加，主要是由于回復段的加載力減小造成的。

(3)所使用的超-粘彈性模型，能夠較準確地捕捉材料在單軸拉伸和剪切變形時的非線性、速率相關性和遲滯行為，計算不同速率剪切應力-應變曲線時最大誤差不超過8%，而計算遲滯曲線時應力平均誤差為9.80%，遲滯曲線包圍面積計算誤差為21.3%，并且對于不同的加載條件也能符合試驗規(guī)律，在遲滯模型中，參數(shù)和的增大都能增加粘性的遲滯效應；以此為基礎，可以指導柔性接頭的動態(tài)計算建模。