桂曉波，余 陵，蔡文祥，連 立

(1.南京理工大學(xué) 機械學(xué)院，南京 210094；2.淮海工業(yè)集團第三研究所，長治 046000)

0 引言

固體火箭發(fā)動機點火過程包含一系列復(fù)雜的物理化學(xué)反應(yīng)，點火過程時間較短，但發(fā)生故障的概率卻相對較高。自由裝填藥柱能夠方便地更換裝藥，不存在粘結(jié)脫粘問題，大大延長了武器系統(tǒng)的存貯時間，因此得到廣泛應(yīng)用。發(fā)動機點火過程中自由裝填藥柱受到瞬態(tài)高壓燃氣壓力作用發(fā)生變形，藥柱結(jié)構(gòu)完整性容易受到破壞，這是典型的流固耦合現(xiàn)象。近10年來，國內(nèi)外逐步采用計算模擬手段對流固耦合現(xiàn)象進行分析。例如，John Montesano[1]提出了一種簡化的單向流固耦合計算模型；曹琪等人[2-3]進行了貼壁澆注藥柱的流固耦合模擬。目前，國內(nèi)很少有人對自由裝填藥柱點火過程中的應(yīng)力應(yīng)變狀態(tài)進行分析，自由裝填藥柱內(nèi)孔和外環(huán)同時受到壓力作用，且壓力變化和分布各不相同，對此進行研究，有助于了解發(fā)動機故障發(fā)生的內(nèi)在原因，并為發(fā)動機結(jié)構(gòu)設(shè)計、裝藥設(shè)計和材料選擇提供依據(jù)。冷流沖擊實驗主要利用高壓氣體模擬固體火箭發(fā)動機點火壓強峰值，作為極限載荷作用在藥柱上，對其進行流動沖擊，為沖擊載荷條件下藥柱結(jié)構(gòu)完整性的研究提供實驗條件，為數(shù)值分析提供驗證數(shù)據(jù)。

本文利用System Coupling耦合器作為流體計算軟件FLUENT和結(jié)構(gòu)分析軟件ANSYS的數(shù)據(jù)交換平臺，針對大長徑比的自由裝填藥柱固體火箭發(fā)動機點火瞬態(tài)過程展開研究，對冷流沖擊過程高壓氣體流動和推進劑藥柱變形之間的相互作用進行了數(shù)值仿真，得到高壓氣體在火箭發(fā)動機內(nèi)的流動狀況，以及藥柱在壓力沖擊下的應(yīng)力與變形的變化規(guī)律，從而分析出容易失效的位置。

1 流固耦合數(shù)學(xué)物理模型

1.1 物理模型

冷流沖擊模擬固體火箭發(fā)動機點火過程實驗臺主要由高壓氣體暫存式儲氣罐、過渡體及實驗發(fā)動機等部分組成，根據(jù)實驗裝置建立的物理模型如圖1所示，采用三維1/4對稱模型。

圖1 流固耦合物理模型示意圖

為簡化模型，不考慮殼體、襯層和絕熱層。發(fā)動機總長914 mm，直徑90 mm，長徑比>10，屬于大長徑比發(fā)動機范疇。儲氣罐內(nèi)徑96 mm，計算時選取長度200 mm。

1.2 流場計算模型

流場區(qū)域進行分區(qū)建模，將氣源、過渡體、燃燒室、噴管等分成12個分區(qū)，使用Workbench Mesh前處理軟件對其進行網(wǎng)格劃分，采用六面體占主體，生成對象如圖2所示，共有單元219 041個，節(jié)點190 065個。

(a)前段網(wǎng)格

(b)后段網(wǎng)格

1.2.1 邊界條件及初始條件

(1)儲氣罐設(shè)置為5 MPa壓力氣源；

(2)噴管出口壓強為1個大氣壓；

(3)1/4對稱面設(shè)置為周期性對稱邊界；

(4)裝藥頭部、外壁、內(nèi)壁為耦合邊界。

計算模型的初始條件：p=0.101 325 MPa，T=295 K，速度為0，湍流動能k和動能耗散率ε取小值。

1.2.2 流場計算控制方程和計算方法

三維氣相湍流流動守恒方程組包括連續(xù)性方程、動量方程和RNGk-ε方程，可表達成通用形式：

(1)

式中Φ為氣相通用變量，在連續(xù)性方程中取1，在動量方程中表示u、v、w3個方向的速度，在RNGk-ε湍流模型中表示湍動能k和耗散率ε；S表示源項；Γ為擴散系數(shù)。

控制方程：

(2)

(3)

式中Gk表示由于平均速度梯度引起的湍動能產(chǎn)生；Gb表示由于浮力影響引起的湍動能產(chǎn)生；YM表示可壓縮湍流脈動膨脹對總的耗散率的影響；αk和αε分別為湍動能k和耗散率ε的有效湍流普朗特數(shù)的倒數(shù)；C1ε、C2ε、C3ε為常數(shù)。

湍流粘性系數(shù)計算公式為

(4)

流場計算選用SIMPLE算法，RNGk-ε模型，時間步長為10-5s。

1.3 結(jié)構(gòu)場計算模型

固體推進劑藥柱為前后不對稱的單孔管狀裝藥，藥柱后段外通道沿裝藥軸線方向減小，尺寸見表1。結(jié)構(gòu)場模型網(wǎng)格劃分采用六面體結(jié)構(gòu)化網(wǎng)格，整個藥柱生成單元19 800個，節(jié)點92 778個，見圖3。

表1 裝藥結(jié)構(gòu)參數(shù)

圖3 藥柱網(wǎng)格示意圖

1.3.1 材料特性

固體推進劑藥柱是一種典型的時間溫度相關(guān)的粘彈性材料，在較低應(yīng)力作用下近似線性粘彈性，其泊松比接近0.5，與不可壓縮材料較為相似。這里采用線性粘彈性本構(gòu)模型，基本可滿足工程的需要，其拉壓松弛模量可用Prony級數(shù)表示[4]：

其中，泊松比為0.495；瞬態(tài)彈性模量為2 671.2 MPa；密度為1 775 kg/m3。

1.3.2 計算方法

瞬態(tài)動力學(xué)分析求解的基本運動方程為

[M]{u″}+[C]{u′}+[k]{u}={F(t)}

(5)

式中 [M]表示質(zhì)量矩陣；[C]表示阻尼矩陣；[k]表示剛度矩陣； {u″}表示節(jié)點加速度向量； {u′}表示節(jié)點速度向量； {u}表示節(jié)點位移向量； {F(t)}表示t時刻的載荷向量。

計算時，將藥柱后端面約束固定，前端面、外壁面和內(nèi)壁面為流固耦合面，兩個對稱面；采用Newton-Raphson法求解，時間步長與流場迭代步相同；每次迭代步計算后，利用System Coupling耦合器在FLUENT和ANSYS之間進行數(shù)據(jù)交換。

1.4 流固耦合方法

本文數(shù)值模型，推進劑具有較大的長徑比，在發(fā)動機點火流場建立過程中，推進劑裝藥結(jié)構(gòu)場變形很小，從整個發(fā)動機燃燒室尺度出發(fā)，可認為固體域的變形對流體域的影響可忽略不計，而流場的建立對結(jié)構(gòu)分析有重大影響，故在研究中，流固耦合分析采用單向耦合分析法，僅考慮冷流沖擊流場對裝藥結(jié)構(gòu)的影響。在沖擊瞬間，冷流氣體升溫幅度并不大，可認為藥柱處于常溫下，不涉及粘彈性材料隨溫度變化關(guān)系，只研究藥柱在峰值壓力作用下隨時間變化的受力情況。

System Coupling作為ANSYS workbench解決流固耦合問題的求解器，可按照設(shè)定的順序在FLUENT和ANSYS中分別進行流體計算和固體計算，通過流固耦合面(FSI Interface)把流體域和固體域的計算結(jié)果相互交換傳遞。待某一時刻的流固計算均達到收斂要求，進行下一時刻的計算，依次而行求解得到隨時間變化的流場及結(jié)構(gòu)場的變化過程。具體計算流程[5]如圖4所示。

2 計算結(jié)果及分析

通過仿真計算可知，在0.35 ms時，流場壓力前鋒抵達藥柱前端面，流場開始對藥柱產(chǎn)生作用力。圖5(a)為0.35 ms時刻的流場速度馬赫數(shù)云圖，高壓氣體進入燃燒室，在藥柱前端面擴展區(qū)形成高度欠膨脹超聲速射流[6]，進入內(nèi)通道形成高速流區(qū)。圖5(b)為1.15 ms時刻內(nèi)通道氣流進入后腔形成約束管內(nèi)的射流波系結(jié)構(gòu)，內(nèi)通道氣流壓力波傳遞快于外通道，當(dāng)內(nèi)通道壓力波前鋒到達后端面時，外通道壓力前鋒還沒有到達，此時的內(nèi)外通道壓力差導(dǎo)致裝藥呈內(nèi)壁膨脹狀態(tài)。圖5(c)為1.40 ms時刻內(nèi)外通道氣流在后腔相遇，此時波系混亂，后腔內(nèi)流場狀況非常復(fù)雜。

圖4 流固耦合流程示意圖

(a)0.35 ms

(b)1.15 ms

(c)1.40 ms

經(jīng)過耦合傳遞，流場壓強作用在藥柱推進劑上，產(chǎn)生應(yīng)力和變形。圖6為0.35 ms時刻藥柱應(yīng)力與變形分布圖。此時流場壓力前鋒剛進入藥柱內(nèi)通道，正面軸向沖擊使前端面內(nèi)孔邊緣出現(xiàn)最大變形，變形量為0.005 7 mm，這一階段變形較小。最大應(yīng)力出現(xiàn)在內(nèi)孔離前端面距離5 mm左右的內(nèi)壁面，這主要是因為軸向壓力對前端面的擠壓，應(yīng)力值為0.657 3 MPa。

(a)應(yīng)力分布

(b)變形分布

圖7為1 ms時刻藥柱應(yīng)力與變形云圖。這時，左端藥柱頭部內(nèi)壁應(yīng)力4.119 9 MPa，底部外壁邊緣，也就是后擋藥板爪卡固定處，應(yīng)力達4.119 6 MPa，如圖7(a)所示。由于壓力波的傳遞及內(nèi)外壓力差作用，藥柱最大應(yīng)力位置隨之轉(zhuǎn)移至藥柱底部外壁邊緣。這時的藥柱頭部最大變形量為0.390 49 mm。從圖7(b)可知，藥柱變形主要為頭部冷流沖擊所導(dǎo)致的軸向變形。

(a)應(yīng)力分布

(b)變形分布

圖8為藥柱最大變形量隨時間變化曲線，呈震蕩變化趨勢，在1.73 ms時，變形量上升到最大值，為1.078 5 mm，位置在藥柱頭部端面。由于粘彈性材料彈性，最大變形量呈震蕩變化，彈性模量決定了最大變形的大小。隨正面沖擊的冷流氣體壓力下降，震蕩幅度也逐漸減小，變形量也隨之減小。

圖8 最大變形量隨時間變化曲線

圖9給出藥柱最大應(yīng)力和底部外壁邊緣應(yīng)力隨時間變化曲線。底部外壁邊緣應(yīng)力呈震蕩變化，同樣是粘彈性材料彈性特性的體現(xiàn)。兩者重合部分即最大應(yīng)力由頭部轉(zhuǎn)移至底部時刻。1.91 ms時，藥柱底部外壁邊緣的應(yīng)力最大，其值為13.32 MPa，最大應(yīng)變0.005 115 5(見圖10，左端為底部)。因此，整個沖擊過程中，底部外壁邊緣和頭部內(nèi)壁最易發(fā)生結(jié)構(gòu)完整性破壞，在結(jié)構(gòu)設(shè)計和材料選擇上應(yīng)加以重視。

實驗過程中，在藥柱頭部外壁、頭部內(nèi)壁、頭部正面及底部外壁等5處粘貼電阻應(yīng)變片，進行多次測量，得到的應(yīng)變數(shù)據(jù)也反映了在頭部內(nèi)壁和底部外壁處的應(yīng)變要比其他所有監(jiān)測點大，從而驗證了數(shù)值分析的正確性。藥柱的失效準則一般采用八面體剪切理論[7]，其失效函數(shù)為

(6)

式中γ8m為八面體剪應(yīng)變的極限值；μ為泊松比；εe為最大等效應(yīng)變；εM為藥柱單向拉伸的最大延伸率,可由實驗得到。

文中使用的藥柱延伸率為35%，遠大于最大等效應(yīng)變，Z′>0。因此，藥柱沒有發(fā)生失效。

圖9 藥柱最大應(yīng)力和底部外壁邊緣應(yīng)力隨時間變化曲線

(a)應(yīng)力分布

(b)應(yīng)變分布

3 結(jié)論

(1)仿真結(jié)果表明System Coupling耦合器能夠很好的將FLUENT和ANSYS結(jié)合進行耦合計算，同時同步計算流場和結(jié)構(gòu)場，更能準確模擬冷流沖擊實驗過程的實際工作狀況。

(2)藥柱變形沿軸向逐漸減小，材料的彈性特性使最大變形量隨時間呈震蕩變化，瞬態(tài)彈性模量決定最大變形的大小。

(3)在壓力上升及傳播過程中，自由裝填藥柱頭部內(nèi)孔壁及底部外壁邊緣Mise應(yīng)力較大，最大應(yīng)力在這2個位置交替變化，這與實驗測得的數(shù)據(jù)相一致。

(4)在點火壓力峰值階段，最大Mise應(yīng)變產(chǎn)生在底部外壁邊緣，易造成藥柱結(jié)構(gòu)完整性破壞。因此，對此處的結(jié)構(gòu)設(shè)計和力學(xué)性能要求較高。

參考文獻：

[1] John Montesano,Kamran Behdinan,David R Greatrix,et al.Internal chamber modeling of a solid rocket motor: effects of coupled structural and acoustic oscillations on combustion[J].Journal of Sound and Vibration,2008,311: 20-38.

[2] 曹琪，李進賢，唐金蘭.帶徑向翼槽SRM點火瞬間流固耦合數(shù)值分析[J].固體火箭技術(shù)，2009，26(12).

[3] 高雙武.含缺陷固體火箭發(fā)動機點火瞬態(tài)流固耦合數(shù)值模擬[D].西安高科技研究所,2007.

[4] 曹杰.自由裝填固體火箭發(fā)動機裝藥點火沖擊特性研究[D].南京理工大學(xué),2013.

[5] 崔小強.激波與固體火箭發(fā)動機裝藥裂紋流固耦合基理研究[D].國防科學(xué)技術(shù)大學(xué),2007.

[6] Purcell Steven P,Daines Weldon L,Christensen Leon W.Simulation of the pressure gradient in a segmented solid rocket motor during the ignition transient[C]//Joint Propulsion Conference and Exhibit,1992: 24.

[7] 張書俊.固體火箭發(fā)動機粘彈性藥柱結(jié)構(gòu)可靠性分析[J].固體火箭技術(shù)，2006，29(3).