賈樂凡，史宏斌，沙寶林，許 楊

(西安航天動力技術(shù)研究所，西安 710025)

0 引言

HTPB推進(jìn)劑被廣泛用于各型航天器、空間飛行器固體發(fā)動機(jī)上[1]。其力學(xué)響應(yīng)行為與溫度、加載方式、加載歷程等諸多因素相關(guān)[2]。國內(nèi)外針對HTPB推進(jìn)劑本構(gòu)方程進(jìn)行了大量研究。HTPB推進(jìn)劑基于本身的高填充比和固化彈性體結(jié)構(gòu)表現(xiàn)為典型的非線性粘彈性復(fù)合材料，此類材料在應(yīng)變較大時呈現(xiàn)明顯的應(yīng)力-應(yīng)變非線性特征[3]。針對HTPB推進(jìn)劑本構(gòu)方程的建立方法主要包括單獨(dú)使用彈性/粘彈性/超彈性模型[4-5]、彈性模型耦合蠕變模型[6-7]、粘彈性模型耦合超彈性模型[8-9]。此外，還有HTPB推進(jìn)劑本構(gòu)方程耦合塑性模型[10]、老化模型[11]、損傷模型[12]等其他模型的研究。目前，大部分本構(gòu)方程建立的核心是表征HTPB材料的粘彈性響應(yīng)特征。而Prony級數(shù)方程作為一種較為簡單的線性粘彈性本構(gòu)方程得到了廣泛應(yīng)用。在基于Prony級數(shù)的模型建立過程中，首先需要確定Prony級數(shù)階數(shù)，隨著Prony級數(shù)階數(shù)增高，其占用的計算資源呈指數(shù)增加。因此，在實(shí)際應(yīng)用過程中Prony級數(shù)階數(shù)并非越高越好，而要結(jié)合數(shù)值模擬目的與計算資源現(xiàn)狀選取合適的Prony級數(shù)的階數(shù)。一般工程計算中這一取值由經(jīng)驗(yàn)給出，未見針對不同階數(shù)的Prony級數(shù)對數(shù)值模擬結(jié)果影響的專門研究。

本文根據(jù)試驗(yàn)獲得HTPB松弛模量主曲線，建立其5～10階Prony級數(shù)模型并進(jìn)行數(shù)值模擬，探究不同階數(shù)的Prony級數(shù)對HTPB推進(jìn)劑數(shù)值模擬研究產(chǎn)生的影響，并嘗試給出不同情況下適用的級數(shù)階數(shù)。所得結(jié)論可應(yīng)用于HTPB推進(jìn)劑相關(guān)數(shù)值模擬研究，為其相關(guān)工程應(yīng)用提供借鑒。

1 粘彈性本構(gòu)方程

對于簡單的粘彈性材料，其應(yīng)力-應(yīng)變曲線與時間具有相關(guān)性。在此基礎(chǔ)上，如果其本構(gòu)關(guān)系符合疊加性原理：

σ(ε1+ε2)=σ(ε1)+σ(ε2)

(1)

以及齊次性原理：

σ(βε)=βσ(ε1)，(β=const)

(2)

則可將其視為線性粘彈性材料。

對于簡單的線性粘彈性材料，其應(yīng)力與應(yīng)變歷史有關(guān)：

(3)

式中E(t)為松弛模量，通過定應(yīng)變試驗(yàn)測得。

對于與溫度相關(guān)的線性粘彈性材料，需要在積分中將溫度依賴項(xiàng)與時間依賴項(xiàng)分離：

(4)

式中ξ為縮減時間。

(5)

式中αT為時間-溫度平移因子，是關(guān)于溫度的函數(shù)，其值通過定應(yīng)變試驗(yàn)擬合而出。

目前，常用的線性粘彈性模型主要包括Maxwell模型、Voigt模型等[13]。其中，Prony級數(shù)模型，作為一種廣義Maxwell模型在工程上應(yīng)用最廣。該模型可視為一個胡克彈簧與若干個Maxwell模型并聯(lián)，其階數(shù)可等效視為并聯(lián)的Maxwell模型個數(shù)。常用的Maxwell模型可視作1階Prony級數(shù)模型。無量綱形式的Prony級數(shù)為

(6)

式中g(shù)R(t)為無量綱松弛模量；N為方程階數(shù)；τi為材料的松弛時間。

對基于Prony級數(shù)的粘彈性本構(gòu)方程，其數(shù)值計算結(jié)果與其階數(shù)N有關(guān)，在實(shí)際應(yīng)用過程中，需結(jié)合數(shù)值模擬目的與計算資源現(xiàn)狀選取合適的級數(shù)階數(shù)。本文針對HTPB推進(jìn)劑粘彈性本構(gòu)方程，研究不同階數(shù)對其數(shù)值模擬研究產(chǎn)生的影響。

2 數(shù)值模擬

2.1 物理模型

通過單向拉伸試驗(yàn)得到的HTPB推進(jìn)劑應(yīng)力-應(yīng)變曲線，獲得其松弛模量主曲線，如圖1所示。

圖1 THPB應(yīng)力松弛模量主曲線(Ts=293.15K)

根據(jù)松弛模量主曲線及時間-溫度平移因子擬合得到無量綱Prony級數(shù)各階的松弛模量系數(shù)，如表1所示。以表1數(shù)據(jù)作為材料模型，根據(jù)單軸拉伸試驗(yàn)標(biāo)準(zhǔn)建立數(shù)值仿真模型。圖2為試件幾何尺寸，圖3為其模型示意圖。

圖2 試件模型(單位：mm)

圖3 啞鈴形試件模型示意圖

表1 各階Prony級數(shù)參數(shù)

分別將5～10階Prony級數(shù)方程作為試件材料參數(shù)代入，對模型一端施加固定邊界條件，另一端施加60 mm位移邊界條件。分別設(shè)置分析步時間間隔為3600、360、36、12、7.2 s，以模擬試件在試驗(yàn)時所經(jīng)歷各拉伸速率下的力學(xué)響應(yīng)得到其應(yīng)力-應(yīng)變曲線，并與試驗(yàn)結(jié)果進(jìn)行比對。

2.2 仿真結(jié)果

圖4為Prony級數(shù)取5～10階時模型在1、10、100、300、500 mm/min拉伸速率下的應(yīng)力-應(yīng)變曲線。

(a)N=5 (b)N=6

由圖4可以看出，各階級數(shù)模型仿真結(jié)果均隨拉伸速率增高而增大。在5階時，模型在500 mm/min拉伸速率下應(yīng)力最大值約0.79 MPa，在1 mm/min拉伸速率下應(yīng)力最大值僅有0.15 MPa；當(dāng)拉伸速率較低時，其曲線斜率變化不大，而當(dāng)拉伸速率較高時，其整體拉伸曲線呈現(xiàn)典型的雙線性特征。這是因?yàn)镻rony級數(shù)模型在較短時間歷程的拉伸過程中主要表現(xiàn)彈性體特征，而在較長時間歷程內(nèi)表現(xiàn)蠕變特性。

此外，綜合分析各階級數(shù)仿真結(jié)果可以看出，拉伸速率較高時，各階Prony級數(shù)仿真結(jié)果較為接近；在拉伸速率較低時，Prony級數(shù)階數(shù)越高，低拉伸速率下仿真結(jié)果越大。這是因?yàn)樵诶焖俾瘦^低(1 mm/min或10 mm/min)時，模型的力學(xué)響應(yīng)受蠕變特性影響較大，高階Prony級數(shù)能夠更好地模擬材料的蠕變行為。因此，其結(jié)果隨階數(shù)升高而增大。

3 試驗(yàn)結(jié)果與討論

3.1 試驗(yàn)儀器與方案

試驗(yàn)儀器為深圳三思縱橫UTM520HB電子萬能試驗(yàn)機(jī)，試驗(yàn)對象為依照QJ 924—1985《復(fù)合固體推進(jìn)劑單向拉伸試驗(yàn)方法》[14]制成的某配方HTPB推進(jìn)劑試件，其配方如表2所示。

表2 HTPB推進(jìn)劑配方

在20 ℃條件下進(jìn)行拉伸速率為1、10、100、300、500 mm/min的HTPB推進(jìn)劑單軸拉伸試驗(yàn)，得到試件的應(yīng)力-應(yīng)變曲線。

3.2 試驗(yàn)結(jié)果

試驗(yàn)在每個拉伸速率下均進(jìn)行3次。其中，除一組試件在300 mm/min拉伸速率下提前斷裂(分析排除試驗(yàn)設(shè)計原因)外，其余各組在相同拉伸速率下的結(jié)果均十分接近。在兩組完整數(shù)據(jù)中，選取一組結(jié)果作為后續(xù)數(shù)值模擬研究參考。試驗(yàn)獲得的HTPB推進(jìn)劑不同拉伸速率下應(yīng)力-應(yīng)變曲線如圖5所示?？梢钥闯?，HTPB推進(jìn)劑在不同拉伸速率下表現(xiàn)出明顯的粘彈性特征，其應(yīng)力應(yīng)變與加載歷程有關(guān)；在同樣拉伸應(yīng)變下，應(yīng)力隨應(yīng)變率的增大而增大。

圖5 不同拉伸速率下推進(jìn)劑應(yīng)力-應(yīng)變曲線

結(jié)合試驗(yàn)結(jié)果可以看出，粘彈性材料拉伸過程中，應(yīng)力-應(yīng)變曲線與拉伸速率有關(guān)，需要分析不同拉伸速率對各階Prony級數(shù)方程數(shù)值模擬結(jié)果的影響。

圖6為不同拉伸速率下的各階Prony級數(shù)模型仿真結(jié)果與試驗(yàn)結(jié)果對比。比較不同拉伸速率下的試驗(yàn)與數(shù)值模擬結(jié)果看出，不同拉伸速率下Prony級數(shù)模型仿真結(jié)果均隨階數(shù)的增高而向試驗(yàn)結(jié)果逼近，這是因?yàn)楦唠A的Prony級數(shù)模型能夠更好地反映材料模型的實(shí)際力學(xué)特性。

(a)1 mm/min (b)10 mm/min

由圖6可見，當(dāng)拉伸速率較高時，各階Prony級數(shù)模型結(jié)果差別不大；而當(dāng)拉伸速率較低時，仿真結(jié)果誤差隨著階數(shù)降低而快速增大。這是因?yàn)镻rony級數(shù)模型瞬時響應(yīng)呈現(xiàn)彈性體特征，而時間效應(yīng)呈現(xiàn)粘性體特征。當(dāng)拉伸速率較高時，材料主要表現(xiàn)彈性效應(yīng)，因此不同階模型結(jié)果差異較?。欢?dāng)拉伸速率較低時，高階級數(shù)模型能夠更好地表征材料的粘性響應(yīng)，因此高階模型相比低階模型誤差明顯減小。此外，可以看到試驗(yàn)中試件相較仿真總是更早進(jìn)入損傷階段。這是因?yàn)橄啾扔趯?shí)際情況，仿真模型無法反映試件在成型過程中的微觀缺陷或損傷，因此其結(jié)果更加理想化。

3.3 不同拉伸速率下各階Prony級數(shù)模型仿真誤差比較

在單軸拉伸試驗(yàn)中，其應(yīng)力通過測力裝置讀數(shù)除以截面面積計算得到，為試驗(yàn)原始數(shù)據(jù)。因此，以試驗(yàn)所測得應(yīng)力值作為基準(zhǔn)分析模型仿真誤差，取每一有效試驗(yàn)應(yīng)力數(shù)據(jù)點(diǎn)與各個仿真結(jié)果應(yīng)力數(shù)據(jù)點(diǎn)距離最小值為其誤差，由此得出不同拉伸速率下的各階Prony級數(shù)模型仿真誤差。

圖7反映了不同拉伸速率下各階模型仿真誤差之間的關(guān)系?？梢钥闯?，不同拉伸速率下Prony級數(shù)模型仿真誤差總體上均隨階數(shù)的增高而減小。當(dāng)拉伸速率較低時，階數(shù)對仿真誤差的影響尤為明顯，階數(shù)過低時，其平均相對誤差甚至大于50%。隨著拉伸速率的增大，仿真誤差對階數(shù)取值的敏感性逐漸降低。當(dāng)拉伸速率在100 mm/min以上時，其各階平均相對誤差均能保持在15%以內(nèi)。因此，當(dāng)需要計算的物理情形下應(yīng)變率較高(拉伸速率較高)時，可以適當(dāng)降低模型的階數(shù)來提高計算效率，而當(dāng)應(yīng)變率較低(拉伸速率較低)時，需要適當(dāng)提高模型階數(shù)來保證計算結(jié)果精度。

圖7 不同拉伸速率下各階級數(shù)仿真誤差比較

進(jìn)一步觀察數(shù)據(jù)可以看出，拉伸速率較低時，低階方程誤差極大，但其誤差隨著階數(shù)增高迅速減??；當(dāng)階數(shù)達(dá)到10階時，模型仿真誤差受拉伸速率變化影響較小。各階方程誤差大體上隨拉伸速率增大而減小，同一拉伸速率下的方程誤差大致上隨階數(shù)增高而減小。當(dāng)拉伸速率較低時，小于9階的方程誤差較大，不能符合實(shí)際情況；當(dāng)拉伸速率較大時，6階以上的模型誤差均相對較小，但考慮到實(shí)際仿真時的情況，即使拉伸速率較大時，基于Prony級數(shù)的粘彈性模型階數(shù)也應(yīng)取不小于7階。

4 結(jié)論

(1)利用Prony級數(shù)計算HTPB推進(jìn)劑粘彈性力學(xué)響應(yīng)時，其階數(shù)最好不小于7階；在小于7階時，模型數(shù)值仿真計算結(jié)果誤差大多偏離試驗(yàn)數(shù)據(jù)10%以上，與材料實(shí)際力學(xué)行為有較大差異；當(dāng)計算貯存、運(yùn)輸?shù)染徛虞d工況時，其階數(shù)最好不小于9階。

(2)當(dāng)拉伸速率較高時，各階級數(shù)仿真結(jié)果與試驗(yàn)相差較小，但隨著拉伸速率下降，低階Prony級數(shù)的仿真誤差迅速增大。

(3)當(dāng)拉伸速率較高時，可采用較低階數(shù)的Prony級數(shù)模型以提高計算效率；而當(dāng)拉伸速率較低時，需要采用高階Prony級數(shù)方程以更好地模擬材料蠕變特性。

本文結(jié)論適用范圍為未進(jìn)入損傷階段的HTPB推進(jìn)劑。其結(jié)論可用于HTPB推進(jìn)劑的簡單工程計算及強(qiáng)度校核，并為其安全裕度設(shè)計及失效判定提供理論依據(jù)。