周東謨，謝旭源，王瑞民，劉向陽，惠步青

（1.中北大學(xué)機電工程學(xué)院， 山西 太原 030051； 2.北京理工大學(xué)宇航學(xué)院， 北京 100081）

0 引 言

復(fù)合固體推進劑固化是指在一定的溫度條件下，推進劑中的黏合劑體系化學(xué)交聯(lián)形成網(wǎng)狀大分子，液態(tài)的推進劑藥漿固化成型的過程［1］。降溫過程主要是靠發(fā)動機殼體散熱使推進劑由固化溫度降低到使用溫度的過程。對于貼壁澆注的固體火箭發(fā)動機，裝藥在固化成型后通常伴隨殘余應(yīng)力/應(yīng)變載荷的產(chǎn)生［2-3］，這種載荷一直伴隨著發(fā)動機的整個服役過程。固體推進劑耐負荷能力差，這些載荷的存在會引起裝藥力學(xué)性能發(fā)生變化，嚴重時會造成推進劑藥柱及界面的裂紋和脫粘，最終導(dǎo)致固體發(fā)動機裝藥結(jié)構(gòu)完整性破壞［4］。因此準(zhǔn)確分析固化降溫作用下發(fā)動機藥柱的應(yīng)力/應(yīng)變響應(yīng)，對保證發(fā)動機安全可靠的應(yīng)用具有重要意義。

推進劑在固化過程中會發(fā)生不同程度的熱膨脹和固化收縮，從而引起熱應(yīng)力和固化收縮應(yīng)力。貼壁澆鑄式發(fā)動機產(chǎn)生的殘余應(yīng)力/應(yīng)變主要有以下4 點原因［5-6］：（1）黏合劑體系在固化過程中自身的體積收縮引起的固化收縮應(yīng)力/應(yīng)變；（2）推進劑固化溫度與發(fā)動機的使用溫度不同；（3）推進劑和發(fā)動機殼體兩者熱膨脹系數(shù)有差異；（4）藥柱外表面與剛性發(fā)動機殼體的內(nèi)表面相粘結(jié)，藥柱的位移變形受到限制。當(dāng)推進劑從高粘性流體狀物體轉(zhuǎn)變成固態(tài)時，黏合劑體系形成交聯(lián)網(wǎng)絡(luò)結(jié)構(gòu)，推進劑產(chǎn)生固化體積收縮，大多數(shù)推進劑在固化起始點（凝膠點）之前開始收縮，但推進劑在凝膠點之前呈流體性質(zhì)，模量較小，不能承受大的剪應(yīng)力，因此其收縮應(yīng)力較?。?］。推進劑在凝膠點之后仍有一定的體積收縮，稱為凈體積收縮，由于受諸多因素的影響，有時凈體積收縮變得很大［8］。這種體積收縮受發(fā)動機殼體的限制，在固化的發(fā)動機藥柱內(nèi)導(dǎo)致產(chǎn)生殘余應(yīng)力，這種在固化過程中由自身的體積收縮引起的殘余應(yīng)力被稱為固化收縮應(yīng)力，該應(yīng)力的大小直接取決于凈體積收縮率［9］。

目前在分析固體發(fā)動機固化降溫后的殘余應(yīng)力/應(yīng)變時，主要采用數(shù)值模擬方法開展研究，通常沒有考慮固化階段藥柱固化收縮應(yīng)力/應(yīng)變的形成過程，而是在降溫計算中基于經(jīng)驗將固化收縮的影響折算為等效溫度效應(yīng)［10］，如復(fù)合推進劑固化溫度為50 ℃時，認為推進劑裝藥零應(yīng)力溫度為58 ℃，然后采用熱-力耦合計算方法分析藥柱在降溫過程中的應(yīng)力/應(yīng)變，文獻［11-12］中關(guān)于固化降溫載荷下藥柱的應(yīng)力/應(yīng)變分析即采用此種方法。對于推進劑固化收縮過程鮮有考慮，固化降溫過程中推進劑裝藥的殘余應(yīng)力/應(yīng)變形成歷程尚不明確。因此，固化收縮應(yīng)力問題在藥柱固化降溫力學(xué)響應(yīng)分析過程中應(yīng)得到足夠的重視。

為研究硝酸酯增塑聚醚（NEPE）推進劑在固化與降溫2 個階段的應(yīng)力/應(yīng)變的形成機制，開展了NEPE 推進劑在固化與降溫過程中藥柱殘余應(yīng)力/應(yīng)變數(shù)值模擬研究，分析了NEPE 推進劑藥柱在固化與降溫過程中殘余應(yīng)力/應(yīng)變的發(fā)展歷程，獲得了各階段不同因素對推進劑藥柱殘余應(yīng)力/應(yīng)變的影響規(guī)律。研究結(jié)果對于提高復(fù)合固體推進劑藥柱殘余應(yīng)力/應(yīng)變的評估精度、明晰其形成機理可提供一定的參考。

1 數(shù)學(xué)模型

1.1 熱傳導(dǎo)模型

推進劑固化過程中的熱量主要來源于2 個方面，一是經(jīng)殼體傳遞給推進劑的外界熱量，另一部分是推進劑固化過程中化學(xué)反應(yīng)產(chǎn)生的非線性內(nèi)熱源。根據(jù)Fourier 熱傳導(dǎo)定律和能量平衡原理，推進劑固化過程中的三維瞬態(tài)熱傳導(dǎo)控制方程表示為［13］：

式 中，ρ為 材 料 的 密 度，kg·m-3；Cp為 材 料 的 比 熱 容，J·kg-1·K-1；T和t分別為溫度和時間，K 和s；ki為材料3個方向的導(dǎo)熱率，J·s-1·m-1·K-1；Q為推進劑固化反應(yīng)產(chǎn)生的熱量，J，可表示為：

式中，Hr為單位體積材料固化反應(yīng)釋放的總熱量，kJ·mol-1；α為固化度；dα/dt為材料的固化速率。

NEPE 推進劑的固化反應(yīng)動力學(xué)模型及參數(shù)，選用文獻［14］所采用的Hsich 動力學(xué)模型，可表示為：

該模型采用基于流變儀測試的推進動態(tài)儲能模量表征推進劑的固化動力學(xué)特性，其中G′0、G′、G′∞分別對于固化初始時刻、固化任意時刻、以及固化結(jié)束時的推進劑動態(tài)儲能模量，MPa；β和η為擬合參數(shù)；η、η0對應(yīng)于與溫度有關(guān)的化學(xué)反應(yīng)的特征松弛時間，s；E為活化能常數(shù)，kJ·mol-1；R為普適氣體常數(shù)，其值為8.314 J·mol-1·K-1；T為絕對溫度，K。

熱傳導(dǎo)固化分析的邊界條件采用第三類邊界條件：

式中，Ta為外界換熱溫度，K；h為對流換熱系數(shù)，W·m-2·K-1。根據(jù)對流換熱系數(shù)的經(jīng)驗計算方法［15］，在靜止空氣條件下，自然對流換熱系數(shù)大致范圍為1.0～10.0 W·m-2·K-1。

1.2 熱力學(xué)模型

式中，Vsh為某固化度下的體積收縮率；V′sh為總的固化體積收縮率。

式中，αc為熱膨脹系數(shù)；ΔT為溫度增量。

考慮到NEPE 推進劑藥柱在固化過程的力學(xué)效應(yīng)時間尺度較長，故采用線彈性模型對其在固化成型過程中的力學(xué)響應(yīng)進行分析。

NEPE 推進劑固化階段完成后，進入降溫階段，將降溫階段的推進劑視為粘彈性材料，三維線性粘彈性材料積分本構(gòu)方程的一般形式為［16］：

對于熱粘彈性邊值問題，溫度變化將引起熱應(yīng)力和熱應(yīng)變，粘彈性材料的力學(xué)性能不僅與時間有關(guān)，還和溫度有著強烈的依賴關(guān)系，即為時溫等效原理［17］：

式中，αT為轉(zhuǎn)換因子，采用WLF（Williams Landel Ferry）方程［17］描述：

式中，C1和C2為常數(shù)，與材料屬性相關(guān)，可通過力學(xué)性能試驗獲得；Tref為參考溫度，K。

2 有限元建模

2.1 計算流程

基 于ABAQUS 有 限 元 軟 件 平 臺［18］，首 先 采 用熱-化學(xué)模型模擬推進劑藥柱的生熱、傳熱過程，計算固化階段每個單元節(jié)點的溫度和固化度值；然后采用熱-力模型計算推進劑固化降溫的應(yīng)力/應(yīng)變，通過模塊化的順序耦合方法來完成固化降溫過程的數(shù)值計算。數(shù)值仿真流程如Scheme 1 所示。

Scheme 1 Curing and cooling simulation for propellant

數(shù)值仿真中用到的子程序功能如下：

USDFLD：定義固化度場量和固化速率場量，并實現(xiàn)HETVAL、UMAT 和UEXPAN 子程序之間的數(shù)據(jù)傳遞；

FILM：定義邊界對流換熱條件，控制環(huán)境溫度；

HETVAL：定義固化過程中的反應(yīng)放熱，結(jié)合USDFLD 子程序?qū)⑼七M劑固化過程嵌入到傳熱分析中，以實現(xiàn)溫度場和固化度場的耦合求解；

UEXPAN：定義NEPE 推進劑材料的熱應(yīng)變及其化學(xué)反應(yīng)收縮應(yīng)變，結(jié)合HETVAL、USDFLD 子程序求解得到的溫度場和固化度場計算應(yīng)變；

UMAT：定義推進劑固化過程中的彈性本構(gòu)模型和降溫過程中的粘彈性本構(gòu)模型，結(jié)合UEXPAN 子程序求解得到的熱應(yīng)變及固化收縮應(yīng)變計算藥柱應(yīng)力。

2.2 有限元建模

以圓管型NEPE 推進劑裝藥為研究對象，建立1/4三維有限元模型，模型包含殼體、絕熱層、推進劑藥柱，絕熱層兩端含人工脫粘層?？紤]到芯模的熱傳導(dǎo)率較大，在降溫過程對發(fā)動機裝藥的溫度歷程影響較大，故模型還包含芯模部分。模型網(wǎng)格類型采用三維8 節(jié)點實體單元C3D8H，共37744 個單元，45269 個節(jié)點。圓管型NEPE 推進劑裝藥模型及網(wǎng)格劃分如圖1所示。

圖1 發(fā)動機模型與網(wǎng)格（單位： mm）Fig.1 Motor model and mesh（unit： mm）

2.3 邊界條件與計算工況

邊界條件如下：

（1）殼體-絕熱層-推進劑界面均采用綁定（Tie）約束，分析過程中界面不分離，內(nèi)部熱傳遞模式為接觸熱傳導(dǎo)；

（2）殼體外表面固定，發(fā)動機對稱面設(shè)置對稱約束條件；考慮到芯模表面一般涂抹脫模劑以便于脫模，推進劑與芯模表面設(shè)置為自由面；

（3）芯模兩端、殼體外表面與空氣對流換熱，對流換熱系數(shù)取值為10 W·m-2·K-1［19］。

計算工況如下：

Step-1：固化階段，NEPE 推進劑藥漿、絕熱層、殼體與芯模初始溫度為50 ℃，NEPE 推進劑藥漿在50 ℃的環(huán)境溫度下固化7 d。

Step-2：降溫階段，固化完成的發(fā)動機放置20 ℃室溫環(huán)境，自然冷卻至室溫。

2.4 材料參數(shù)

NEPE 推進劑的固化動力學(xué)參數(shù)取自文獻［14］，如表1 所示。在固化階段，根據(jù)國內(nèi)外學(xué)者提出的固態(tài)高聚合物模量動靜態(tài)轉(zhuǎn)換公式［20］，可將固化階段NEPE 推進劑的模量表示為：

表1 NEPE 推進劑固化動力學(xué)參數(shù)Table 1 Curing kinetic parameters of NEPE propellant

E′(t) ≈G′(ω)|ω=2/(πt)（13）由式（13）可知NEPE 推進劑固化階段的彈性模量可以近似采用固定頻率ω測得的儲能模量。文獻［14］中NEPE 推進劑儲能模量隨固化度的變化曲線如圖2 所示。

圖2 G′隨固化度的變化圖［14］Fig.2 G′ varies along the degree of curing［14］

擬合圖2 數(shù)據(jù)，可獲得NEPE 推進劑固化階段的彈性模量為：

式中，E′為彈性模量，MPa。

NEPE 推進劑藥柱在固化過程中，反應(yīng)熱主要源于NEPE 黏合劑的反應(yīng)，則固化過程中推進劑的反應(yīng)熱為：

式中，Vs為反應(yīng)物的體積占比；Qh為NEPE 黏合劑的反應(yīng)熱，參照文獻［21］取值為134.68 kJ·mol-1。NEPE黏合劑在NEPE 推進劑中的質(zhì)量占比為約為8%，考慮到推進劑藥柱各組成組分密度相差不大，可將NEPE黏合劑在NEPE 推進劑中的質(zhì)量占比視為其體積占比。

在降溫階段，將NEPE 推進劑視為粘彈性材料，其松弛模量采用Prony 級數(shù)表示為：

式中，E（t）為松弛模量，MPa。

取 參 考 溫 度 為20 ℃，C1=-7.053，C2=171.513。其余材料參數(shù)如表2 所示。

表2 材料參數(shù)Table 2 Material properties parameters

3 結(jié)果分析

3.1 固化過程藥柱溫度與固化度

利用ABAQUS 軟件中的Visualization 模塊顯示計算結(jié)果云圖。NEPE 推進劑藥柱在50 ℃環(huán)境溫度下，在固化前期（24 h）、固化中期（72 h）和固化后期（168 h）的溫度場云圖如圖3 所示。

圖3 NEPE 推進劑50 ℃下固化的前中后期溫度場云圖Fig.3 Contours of the temperature under the different curing period of the NEPE propellants at 50 ℃

由圖3 可知，在固化過程中，推進劑藥柱內(nèi)部中心溫度較周邊高。固化中期（圖3b）藥柱內(nèi)部溫度高于固化初期（圖3a）的溫度，在固化結(jié)束時（圖3c），藥柱內(nèi)部溫度趨于均勻，接近環(huán)境溫度50 ℃。顯然這是由于在固化初期和固化后期藥柱固化反應(yīng)速率較慢，在固化中期固化反應(yīng)速率較快從而使得放熱速率也較快，且藥柱中心位置離換熱邊界較遠，存在一定的熱聚集引起的。

推進劑藥柱固化24，72 h 和168 h 時的固化度場云圖如圖4 所示。

圖4 NEPE 推進劑固化度場云圖Fig.4 Contours of the curing degree of the NEPE propllants

由圖4 可知，固化度分布與溫度分布呈現(xiàn)出相同趨勢，在固化初期和固化中期藥柱內(nèi)部溫度相對邊緣較高，因此固化初期（圖4a）和固化中期（圖4b）推進劑藥柱內(nèi)部中心固化度較周邊高。固化結(jié)束時（圖4c）藥柱內(nèi)部固化度趨于均勻，達到0.986 以上。

為進一步分析NEPE 推進劑藥柱內(nèi)部固化度和溫度的變化規(guī)律，根據(jù)溫度和固化度的分布特點選取10個節(jié)點進行分析，如圖5 所示。

圖5 分析節(jié)點的選取Fig.5 Schematic diagram of node selection

各節(jié)點溫度及固化度隨時間的變化曲線如圖6所示。

圖6 不同節(jié)點的固化度和溫度隨時間變化圖Fig.6 Variation histories of curing degree and temperature with curing time for different node selection

由圖6a 可知，在固化20 h 內(nèi)NEPE 推進劑固化度增速較慢；在20 h 到40 h 時，固化度增速逐漸加快，在40 h 時，固化度增速達到最大；在固化40 h 之后，固化度增速逐漸減慢，并在160 h 之后，固化度逐漸趨于穩(wěn)定。在固化25 h 到130 h 間，推進劑內(nèi)部F、G、H、I、J點的固化度明顯比周邊其他點高，最大固化度差在65 h 時達到0.10，在其余時間段內(nèi)各點的固化度差異不是很明顯。

由圖6b 可知，在固化過程中NEPE 推進劑藥柱內(nèi)部溫度整體呈現(xiàn)出先增大，后減小的趨勢，在固化30 h 左右內(nèi)部各點的溫度達到最大，各點之間的溫差也達到最大。位于推進劑藥柱內(nèi)部的F、G、H、I、J 點的溫度較高，而位于推進劑藥柱端部內(nèi)孔邊緣上的A 點和絕熱層處的D 點，基本上溫度維持在環(huán)境溫度50 ℃左右。在固化150 h 之后，推進劑內(nèi)部溫度基本上趨于均勻，約為50 ℃。

推進劑固化整體上是一放熱反應(yīng)，在固化前期與中期，固化程度不高，溫度的不均勻引起藥柱內(nèi)部固化速率分布不均，從而導(dǎo)致藥柱在固化前、中期內(nèi)部存在較大溫度梯度與固化度梯度。在固化后期，藥柱基本固化完成，相應(yīng)的固化反應(yīng)速率很低，藥柱內(nèi)部產(chǎn)熱速率與散熱速率基本平衡，其內(nèi)部溫度也趨于一致。

3.2 固化過程藥柱殘余應(yīng)力應(yīng)變

以等效應(yīng)力與總等效應(yīng)變表征藥柱的殘余應(yīng)力與殘余應(yīng)變。NEPE 推進劑藥柱固化完成時，藥柱總等效應(yīng)力σ與總等效應(yīng)變ε結(jié)果如圖7 所示。

圖7 固化結(jié)束時的應(yīng)力應(yīng)變云圖Fig.7 Contours of stress and strain at the end period of curing

圖8 分析路徑選取圖Fig.8 Schematic map of path selection

圖9 應(yīng)力/應(yīng)變隨路徑變化曲線圖Fig.9 Variation histories of the stress and strain with path 1-4

圖10 應(yīng)力/應(yīng)變隨時間變化曲線圖Fig.10 Variation histories of stress and strain with curing time

3.3 降溫過程藥柱殘余應(yīng)力應(yīng)變

在固化的基礎(chǔ)上對NEPE 推進劑藥柱進行降溫過程進行計算，降溫完成時藥柱總等效應(yīng)力σtot與總等效應(yīng)變εtot云圖如圖11 所示。

圖11 固化降溫結(jié)束時的應(yīng)力應(yīng)變云圖Fig.11 Contours of stress and strain at the end of cooling

由圖11 可知，降溫結(jié)束時，藥柱產(chǎn)生的σtot與εtot的最大值均位于藥柱內(nèi)孔中部，最大σtot為0.0357 MPa，最大εtot為0.033。降溫完成后藥柱的σtot與εtot沿路徑1～4 的變化如圖12 所示。

圖12 應(yīng)力/應(yīng)變隨路徑1～4 變化曲線圖Fig.12 Variation histories of the stress and strain with path 1-4

綜合圖10 與圖13 可知，NEPE 推進劑藥柱在固化與降溫2 個階段的殘余應(yīng)力/應(yīng)變基本符合應(yīng)力/應(yīng)變疊加原理。在固化過程中，藥柱內(nèi)部的溫差會引起固化速率不均勻，影響殘余應(yīng)力/應(yīng)變的形成歷程，但最終藥柱的固化度與溫度趨于一致，因此，藥柱內(nèi)部的溫差不會影響藥柱最終的殘余應(yīng)力和殘余應(yīng)變。

圖13 各應(yīng)力/應(yīng)變隨路徑1 變化曲線圖Fig.13 Variation histories of the stress and strain with path 1

圖14 應(yīng)力/應(yīng)變隨路徑1～4 變化曲線圖Fig.14 Stress and strain varied along path 1-4

表3 進一步給出了在藥柱固化降溫過程中路徑1上中心位置各部分應(yīng)力/應(yīng)變的比例。

表3 各應(yīng)力應(yīng)變占比表Table 3 The stress and strain ratio

由表3 可知，NEPE 推進劑藥柱內(nèi)孔處的殘余應(yīng)力與應(yīng)變，主要是由藥柱降溫載荷引起的，即推進劑固化溫度與藥柱使用溫度不同，在固化成型后的降溫過程中，由于殼體的約束，推進劑產(chǎn)生的熱應(yīng)力占比約為80%，熱應(yīng)變占比約為70%。同時推進劑由于固化體積收縮產(chǎn)生的固化收縮應(yīng)力/應(yīng)變也不容忽視，固化收縮應(yīng)力占比約為18%，固化收縮應(yīng)變占比約為30%。

4 結(jié) 論

開展了考慮NEPE 推進劑在固化降溫過程中藥柱殘余應(yīng)力/應(yīng)變數(shù)值分析，得到以下結(jié)論：

（1）NEPE 推進劑藥柱在50 ℃環(huán)境固化過程中，由于藥柱中心位置離換熱邊界較遠，存在一定的熱聚集，導(dǎo)致藥柱內(nèi)部存在溫度梯度，最高溫度產(chǎn)生在藥柱截面中心位置，在固化結(jié)束時藥柱溫度場趨于與環(huán)境溫度一致。內(nèi)部溫差導(dǎo)致固化過程中固化速率不均勻，造成藥柱截面中心位置固化速率較快，但在固化完成時內(nèi)部固化度趨于一致，藥柱內(nèi)部的溫差不會影響藥柱最終的殘余應(yīng)力和殘余應(yīng)變。

（2）NEPE 推進劑藥柱在固化與降溫過程中的殘余應(yīng)力/應(yīng)變主要由固化收縮應(yīng)力/應(yīng)變與降溫過程產(chǎn)生熱應(yīng)力/應(yīng)變構(gòu)成。對于殘余應(yīng)力，二者占比分別約為20% 與80%，對于殘余應(yīng)變，二者占比分別約為30%與70%?？山普J為，在固化與降溫2 個階段中藥柱內(nèi)部的殘余應(yīng)力/應(yīng)變是由固化收縮應(yīng)力/應(yīng)變與熱應(yīng)力/應(yīng)變線性疊加而成。

（3）由于傳統(tǒng)方法在計算時只將推進劑固化體積收縮率折算為等效溫度而未考慮推進劑的裝藥結(jié)構(gòu)因素，故按傳統(tǒng)方法計算得到的殘余應(yīng)力/應(yīng)變相較本研究計算結(jié)果整體偏大，其中，殘余應(yīng)力最大偏差約為20%，殘余應(yīng)變最大偏差約為10%，其余部分基本一致。