非線性粘彈性本構模型在固體火箭發動機裝藥結構完整性分析中應用

李 毅,田維平,董新剛,姚 東

(中國航天科技集團公司四院四十一所，西安 710025)

0 引言

復合固體推進劑是一種顆粒增強的高聚型含能材料，廣泛用于運載火箭、導彈等裝備的動力系統，其性能直接影響裝備的運載能力、射程、貯存性等關鍵性能。復合固體推進劑力學特性的粘彈屬性與其自身的配方特點有關，也與服役環境如溫度、加載應變率、應力應變狀態等因素高度相關。在變形較小時，可認為其本構模型為線粘彈性；隨著變形量增大，推進劑內部微裂紋、微孔洞等損傷逐漸演化，出現明顯的脫濕現象，此時線性理論不再適用，應選取合適的非線性粘彈性本構模型。

早期基于復合推進劑的線粘彈性本構模型開展了工程領域的結構分析，結合產品研制的經驗修正取得了良好的應用效果[1]。Scahpery[2]及其合作者一直致力于含損傷非線性粘彈性本構模型的研究。Ulrich Mandel等[3]對纖維增強復合材料進行研究，提出一種三維的非線性粘彈性本構模型，并指出此模型可適用于任意載荷情況下。Muliana[4]從細觀角度入手，提出一種非線性粘彈性本構模型。張曉等[5]從應變能出發推導出HTPB推進劑非線性粘彈性本構模型。常新龍等[6]研究了HTPB高應變率下的非線性粘彈性本構模型。王哲君等[7]對HTPB壓應力狀態下的本構模型進行了研究。姚東等[8]從損傷效應的擬合出發，提出了考慮應力狀態的非線性本構模型。唐志平、朱兆祥等[9]提出了被后來的學者稱為朱-王-唐非線性本構模型的成果，均為復合固體推進劑非線性粘彈性本構模型的研究作出了貢獻。

本文基于朱-王-唐非線性本構模型，針對其積分項繁多、本構參數獲取難度較大的特點，提出一種形式簡單、便于工程應用的本構模型。

1 推進劑本構模型的提出

朱-王-唐非線性本構模型如式(1)所示：

σ=E0ε+αε2+βε3+

(1)

式中E0、E1、E2、α、β、θ1、θ2均為材料常數，前3項表示非線性彈性部分，后2項積分項表示粘性部分。

式(1)表明，非線性粘彈性本構模型可寫成彈性部分和粘性部分之和的形式。參考該結論，本文將推進劑本構模型表述為如下形式：

σ(t,ε)=fE(t,ε)+fv(t,ε)

(2)

式中fE(t,ε)為彈性項；fv(t,ε)為粘性項。

Staverman和Schwarzl[10]在大量實驗觀測基礎上指出，采用式(3)：

σ(t)=E(t)f(ε)

(3)

代替σ(t)=E(t)ε可較好地描述實驗結果。趙榮國[11]由此假設過去不同時刻的應變對現時應力的影響可忽略不計，并對其進行了驗證。聯想到應力松弛實驗，本文認為推進劑的受力分為應變變化引起的應力及松弛現象導致的應力耗散。對于應變從0～ε的變化，都可認為是以上2種應力疊加產生的效果，這樣應力的表達形式就可簡單地寫成：

σ(t,ε)=fE(ε)-σ松弛

(4)

式中fE(ε)、σ松弛分別為彈性、粘性應力項，其函數形式及其中的材料參數還有待通過實驗數據進行確定。

2 基于松弛及拉伸實驗的本構模型建模

推進劑在變形較小時，服從線粘彈性規律，當變形達到一定值后，損傷的逐漸演化導致其非線性顯著增加，本文假設應力耗散項不受損傷變量的影響。因此，式(4)中的σ松弛就可用式(5)表示：

σ松弛=[E(0)-E(t)]ε

(5)

式中E(0)為初始時刻的松弛模量值；E(t)為t時刻松弛模量值；ε為t時刻應變。

粘性應力耗散項是狀態量，僅與當前狀態有關，如應變、時間等，不受以往應力-應變狀態的影響。

彈性項的應力-應變關系包括大應變狀態的數據，因此，還需進行一定的拉伸實驗。

2.1 松弛實驗及其分析

推進劑松弛實驗采用啞鈴模型，如圖1所示。

試驗件拉伸到應變為5%，隨后測其應力值，并將數據轉化為松弛模量。為得到完備的結果，需進行多組溫度下的松弛實驗。

根據時溫等效原理E(t,T)=E(t/αT,T0)，對不同溫度下的松弛模量進行處理，使所有數據基本重合為一條曲線，即松弛模量主曲線。隨后，將松弛模量主曲線擬合為6階Prony級數形式：

(6)

式中各參數的值如表1所示。在對不同溫度下松弛模量進行處理時，得到不同溫度對應的時溫等效因子的值，時溫等效因子lgαT僅和溫度有關。根據WLF方程(7)，擬合得到的C1和C2的值分別為8.719 0、144.695 8。擬合后的松弛模量主曲線與實驗數據的對比如圖2所示。

(7)

2.2 拉伸實驗及其分析

拉伸試件同圖1所示模型，實驗在Instron 4210型試驗機上完成，測試溫度 25 ℃，應變率 0.023 8 s-1，測得數據結果如圖3所示。

式(4)中的應力彈性項不包括粘性應力。因此，根據拉伸實驗得到的松弛模量及式(5)，對拉伸數據進行修正，得到彈性應力項的表達式：

(8)

fE(ε)=σ實驗+σ松弛

修正后的應力用式(9)形式的函數進行擬合效果如圖4所示。其中，E0和Es均為材料參數，不同溫度下的E0和Es值如表2所示。

(9)

將E0和Es表示為溫度T的函數關系式，擬合的函數關系式如式(10)、式(11)所示，擬合曲線如圖5所示。

(10)

(11)

3 本構模型及其參數的數值驗證

3.1 本構模型一維MATLAB驗證

本文采用MATLAB編寫代碼對本構模型進行一維數值驗證，MATLAB中僅模擬邊界條件的變化，無需定義單元類型。在應變率0.023 8 s-1、25 ℃的情況下，計算出的應力-應變走勢與實驗數據的對比如圖6所示。可看出，應力隨應變變化的趨勢基本和實驗數據吻合，脫濕點對應的應變約為20%，脫濕后的誤差明顯大于線性段，相對偏差達8%。本文假設應力耗散與損傷無關，損傷的增加有可能會導致粘性應力耗散加大，具體規律需進一步研究。

粘彈性材料的力學性能具有溫度相關性和載荷率相關性，需進行多組不同邊界條件的數值仿真，現將不同條件下的數值仿真結果同實驗的對比匯總如圖7所示。從圖7可看出，初始階段應力應變曲線差別不大，符合線粘彈性規律；在相同拉伸速率條件下，隨著溫度的升高，應力值有一個減小的趨勢，溫度越低，應力水平越高；當溫度不變時，拉伸速率越大，應力值越大。數值擬合結果與實驗、理論都具有相同的規律，相對偏差不超過10%。

用MATLAB對松弛實驗進行數值仿真，結果如圖8所示。結果表明，從拉伸實驗數據推導出的本構模型同樣適用于非拉伸情況，且相對偏差不超過10%。

3.2 本構模型模型三維ABAQUS驗證

利用ABAQUS用戶材料子程序UMAT對本構模型進行二次開發，將本文提出的本構模型編寫成UMAT子程序，在ABAQUS中進行本構模型的三維算例驗證。

將本構模型從一維擴展到三維[12-13]如式(12)。

Aijkl(E(0)-E(t))εkl

(12)

其中，Aijkl為一個四階張量，此處可表示為矩陣形式：

(13)

(14)

粘性項的應力計算可按照線彈性部分進行三維計算，需注意一點：粘性項是狀態量，不可按照增量形式進行計算。

應力更新過程分為2步：(1)計算增量步結束時的彈性應力。t時刻傳入的應力σ(t)為真實應力，欲得到增量步結束時的彈性應力，需先得到增量步開始時的彈性應力值，加上上一步減去的粘性應力項，才可得到增量步結束時的彈性應力。(2)計算增量步結束時的真實應力。用步驟(1)中得到的彈性應力減去增量步結束時的粘性應力，便可得到增量結束時的真實應力。增量步結束時真實應力為

(15)

則Δt時間段內應力增量為

Δσij(t+Δt)=σij(t+Δt)-σij(t)

=[Dij-E(0)+E(t+Δt)]AijklΔεkl(t+Δt)+

[E(t+Δt)-E(t)]Aijklεkl(t)

(16)

令：

(17)

依據上述有限元增量形式，將其編寫為UMAT在ABAQ中進行數值仿真計算。

數值仿真模型尺寸與圖1所示試件尺寸相同，在ABAQUS中建立模型如圖9所示。模型包括試驗件和夾具兩部分。試驗件網格尺寸1 mm，網格數量23 230，單元類型為C3D20H。夾具網格尺寸2 mm，與圖9所示垂直紙面方向尺寸為1 mm，網格數量5410，材質為鋼。計算過程中固定一個夾具，另一夾具勻速向另一方向移動。

推進劑泊松比取為0.496，單元類型選擇雜交單元(Hybrid Formulation)以消除沙漏效應，應力更新過程有UMAT實現。

對不同溫度、不同加載速率條件進行數值仿真，取試件中心節點為研究對象，查看其應力-應變曲線并與實驗數據進行對比，結果如圖10所示。三維數值仿真結果符合推進劑拉伸實驗結果，不同溫度、拉伸速率下的變化規律也符合理論預測。

從圖10可看出，三維結果較一維結果誤差有所增大，相對偏差最大為14.8%，這是因為三維計算在ABAQUS下進行，所用應力應變為真實應力應變，較名義應力應變誤差增大。

誤差原因分析：

(1)單向拉伸試驗得到的應力-應變均為名義值，轉換為真實值，進而擬合本構參數時存在一定偏差；

(2)泊松比按經驗取為0.496，與真實值存在一定偏差；

(3)粘性應力耗散項的計算由松弛模量確定，文中假設這一項與損傷無關，而損傷對推進劑本構模型影響極大。因此，在損傷較大的情況下，松弛模量的形式不變必然導致誤差的存在。

4 模擬發動機熱-結構分析

根據本文提出的非線性粘彈性本構模型，對貼壁澆注固體火箭發動機硫化降溫過程進行分析。計算模型如圖11所示，主要結構參數采取無量綱表示。

在ABAQUS/CAE中建立模型。藥柱、絕熱層、殼體之間采取布爾加運算進行合并。為控制網格對分析結果的影響，對關鍵部位采取了不同的種子控制策略(尺寸及過渡要求)，并對藥柱、絕熱結構等對應材料泊松比較大的部位采用雜交單元技術以消除沙漏效應，如圖12所示。在藥柱兩側施加對稱邊界條件，在后裙端面施加固支條件以消除剛體位移。

推進劑、絕熱層、殼體各部件的材料屬性如表3所示。推進劑模量部分本文提出的本構模型在UMAT進行計算。

熱邊界：藥柱硫化降溫過程是放在一定的溫度環境中，該環境溫度為10 ℃；藥柱、絕熱層、殼體初始溫度均為60 ℃，推進劑和絕熱層與空氣的對流換熱系數為1 K/(m2·W)，殼體與空氣的對流換熱系數為10 K/(m2·W)，時間歷程設置為30 d。

表3 推進劑、絕熱層、殼體的材料屬性

藥漿澆注后的硫化過程涉及粘合劑基體高分子網絡的交聯、粘合劑基體與固體顆粒組分的粘接等一系列化學-熱力學耦合過程，反應機理復雜、影響因素眾多。本文所述硫化降溫過程，以藥漿硫化反應達到“零應力溫度”為起點，并假設：

(1)“零應力溫度”直至目標環境溫度的過程中，藥柱熱物理屬性及力學性能均可由推進劑試件的測試結果進行有效表征；

(2)在“零應力溫度”直至目標環境溫度范圍內，推進劑熱物理屬性不隨溫度變化。

受限于實驗測試條件，在進行硫化降溫實驗中，僅可測試藥柱中孔處的位移變化，因此選取此作為本構模型的驗證標準。

取模型上3個典型位置如圖13所示，其輸出溫度隨時間變化的數據如圖14所示。外壁面與目標環境的對流換熱系數最大，故溫度下降快；藥柱內孔處與目標環境的對流換熱系數較小，因此溫度下降較為緩慢；肉厚中間處熱損失最少，故溫度下降最為緩慢。

硫化降溫過程結構的變化主要是因為溫度場的變化引起的，初始階段溫度變化劇烈，所以中孔位移變化也會較快，當溫度場平衡后，位移也不會發生變化。從圖14可看出，約200 h后，溫度場趨于平衡。因此，中孔位移必然也有類似的規律。將中孔位移隨時間變化的數據繪制在圖15中，數值仿真結果與分析一致。

數值仿真計算得到的中孔位移最大為0.014 48R(R為殼體半徑)。采用相同結構尺寸的發動機進行了10 ℃環境下的溫度平衡試驗，8 d后測得藥柱內孔處溫度與環境溫度一致，表明藥柱溫度達到了平衡；此時進行藥柱內徑測量，并與設計值對比，發現藥柱內孔徑向位移約為0.013 8R。數值仿真與實驗數據的相對偏差僅為4.9%，驗證了本構模型的準確性。

5 結論

(1)本文構建的考慮損傷、溫度的非線性粘彈性本構模型，由彈性項和粘性項組成；彈性項僅為應變的函數，結合松弛實驗和拉伸實驗對參數進行了擬合；粘性項為狀態量，且不受損傷的影響，可由松弛模量和應變計算得到。

(2)通過MATLAB和QBAQUS對本構模型進行了數值仿真驗證，仿真結果與實驗數據吻合較好，一維仿真結果與實驗數據相對偏差不超過10%，三維仿真結果相對偏差略有增大，最大為14.8%。

(3)基于本文所建立的本構模型對發動機硫化降溫過程進行了數值仿真，得到中孔徑向位移為0.144 8R，與實驗測得的數據(0.013 8R)相對偏差僅為4.9%，說明本構模型在發動機熱-結構分析具有較好的精度。

[1] 劉中兵,利鳳祥,李越森,等.高過載條件下固體推進劑藥柱結構完整性分析計算[J].固體火箭技術,2003,26(2):12-16.

LIU Zhongbing,LI Fengxiang,LI Yuesen,et al.A calculation analysis for structural integrity of solid propellant grains under high overload[J].Journal of Solid Rocket Technology,2003,26(2):12-16.

[2] Ha K,Schapery R A.A three dimensional viscoelastic constitutive model for particulate composites with growing damage and its experimental validation[J].International Journal of Solid Structures,1998,35(26-27):3497-3517.

[3] Ulrich Mandel,Robin Taubert,Roland Hinterh?lzl.Three-dimensional nonlinear constitutive model for composites[J].Composite Structures,2016,142:78-86.

[4] Muliana A,Rajagopal K R,scharnuterD T,et al.A nonlinear viscoelastic constitutive model for polymeric solids based on multiple natural configuration theory[J].International Journal of Solid Structures,2016,100-101:95-100.

[5] 張曉,鄭堅,彭威,等.HTPB復合固體推進劑粘彈性應變能及非線性本構模型[J].固體火箭技術,2015,38(6):827-832.

ZHANG Xiao,ZHENG Jian,PENG Wei,et al.Viscoelastic strain energy and nonlinear constitutive model for HTPB composite solid propellant[J].Journal of Solid Rocket Technology,2015,38(6):827-832.

[6] 常新龍,賴建偉,張曉軍,等.HTPB推進劑高應變率粘彈性本構模型研究[J].推進技術,2014,35(1):123-127.

CHANG Xinlong,LAI Jianwei,ZHANG Xiaojun,et al.High strain-rate viscoelastic constitutive model for HTPB propellant[J].Journal of Propulsion Technology,2014,35(1):123-127

[7] 王哲君,強洪夫,王廣,等.中應變率下HTPB推進劑壓縮力學性能和本構模型研究[J].推進技術,2016,37(4):776-782.

WANG Zhejun,QIANG Hongfu,WNAG Guang,et al.Mechanical properties and constitutive model for HTPB propellant under intermediate strain rate compression[J].Journal of Propulsion Technology,2016,37(4):776-782.

[8] 姚東,張光喜,高波.考慮應力狀態的HTPB/AP推進劑含損傷熱-粘彈性本構方程[J].固體火箭技術,2014,34(4):496-499.

YAO Dong,ZHANG Guangxi,GAO Bo.Constitutive equations involving damage for HTPB/AP propellant considering stress state[J].Journal of Solid Rocket Technology,2014,34(4):496-499.

[9] 唐志平,田蘭橋,朱兆祥.高應變率下環氧樹脂的力學性能[C]//全國第二屆爆炸力學會議文集.揚州,1981.

TANG Zhiping,TIAN Lanqiao,ZHU Zhaoxiang.Mechanical properties of epoxy resin under high strain rate[C]//Anthology of the Second National Conference on Mechanics of Explosion.Yangzhou,1981.

[10] Staverman A J,Schwarzl F.Nonlinear deformation behavior of high polymers[M].Stuart H A,Ed.Die Phyaile der Hochpolymeren,1956.

[11] 趙榮國,張淳源.一個考慮損傷的非線性粘彈性本構關系[J].湘潭大學自然學科學報,2001,23(3):28-33.

ZHAO Rongguo,ZHANG Chunyuan.A nonlinear viscoelastic constitutive relation with damage[J].National Science Journal of Xiangtan University,2001,23(3):28-33.

[12] 馮震宙,王新軍,王富生,等.朱-王-唐非線性粘彈性本構模型在有限元分析中的實現及其應用[J].材料科學與工程學報,2007,25(2):269-272.

FENG Zhenzhou,WANG Xinjun,WANG Fusheng,et al.Implementation and its application in finite element analysis of constitutive model for ZWT nonlinear viscoelastic material[J].Joural of Materials Sciense and Engineering,2007,25(2):269-272.

[13] 彭云.基于 ABAQUS 的非線性粘彈性本構模型二次開發[J].南昌航空大學學報,2011,25(2):38-41.

PENG Yun.Developing of nonlinear viscoelastic constitutive model based on ABAQUS[J].Joural of Nanchang Hangkong University (Natural Sciences),2011,25(2):38-41.