自由裝填藥柱冷流沖擊過程中流固耦合數(shù)值分析①

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

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

0 引言

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

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

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

1.1 物理模型

冷流沖擊模擬固體火箭發(fā)動(dòng)機(jī)點(diǎn)火過程實(shí)驗(yàn)臺(tái)主要由高壓氣體暫存式儲(chǔ)氣罐、過渡體及實(shí)驗(yàn)發(fā)動(dòng)機(jī)等部分組成，根據(jù)實(shí)驗(yàn)裝置建立的物理模型如圖1所示，采用三維1/4對(duì)稱模型。

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

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

1.2 流場(chǎng)計(jì)算模型

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

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

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

1.2.1 邊界條件及初始條件

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

(2)噴管出口壓強(qiáng)為1個(gè)大氣壓；

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

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

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

1.2.2 流場(chǎng)計(jì)算控制方程和計(jì)算方法

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

(1)

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

控制方程：

(2)

(3)

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

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

(4)

流場(chǎng)計(jì)算選用SIMPLE算法，RNGk-ε模型，時(shí)間步長(zhǎng)為10-5s。

1.3 結(jié)構(gòu)場(chǎng)計(jì)算模型

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

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

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

1.3.1 材料特性

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

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

1.3.2 計(jì)算方法

瞬態(tài)動(dòng)力學(xué)分析求解的基本運(yùn)動(dòng)方程為

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

(5)

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

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

1.4 流固耦合方法

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

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

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

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

圖4 流固耦合流程示意圖

(a)0.35 ms

(b)1.15 ms

(c)1.40 ms

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

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

(b)變形分布

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

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

(b)變形分布

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

圖8 最大變形量隨時(shí)間變化曲線

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

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

(6)

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

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

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

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

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

3 結(jié)論

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

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

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

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

參考文獻(xiàn)：

[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] 曹琪，李進(jìn)賢，唐金蘭.帶徑向翼槽SRM點(diǎn)火瞬間流固耦合數(shù)值分析[J].固體火箭技術(shù)，2009，26(12).

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

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

[5] 崔小強(qiáng).激波與固體火箭發(fā)動(dòng)機(jī)裝藥裂紋流固耦合基理研究[D].國(guó)防科學(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ā)動(dòng)機(jī)粘彈性藥柱結(jié)構(gòu)可靠性分析[J].固體火箭技術(shù)，2006，29(3).