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

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

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

0 引言

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

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

1 流固耦合數學物理模型

1.1 物理模型

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

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

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

1.2 流場計算模型

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

(a)前段網格

(b)后段網格

1.2.1 邊界條件及初始條件

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

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

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

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

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

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

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

(1)

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

控制方程：

(2)

(3)

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

湍流粘性系數計算公式為

(4)

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

1.3 結構場計算模型

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

表1 裝藥結構參數

圖3 藥柱網格示意圖

1.3.1 材料特性

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

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

1.3.2 計算方法

瞬態動力學分析求解的基本運動方程為

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

(5)

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

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

1.4 流固耦合方法

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

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

2 計算結果及分析

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

圖4 流固耦合流程示意圖

(a)0.35 ms

(b)1.15 ms

(c)1.40 ms

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

(a)應力分布

(b)變形分布

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

(a)應力分布

(b)變形分布

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

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

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

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

(6)

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

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

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

(a)應力分布

(b)應變分布

3 結論

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

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

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

(4)在點火壓力峰值階段，最大Mise應變產生在底部外壁邊緣，易造成藥柱結構完整性破壞。因此，對此處的結構設計和力學性能要求較高。

參考文獻：

[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點火瞬間流固耦合數值分析[J].固體火箭技術，2009，26(12).

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

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

[5] 崔小強.激波與固體火箭發動機裝藥裂紋流固耦合基理研究[D].國防科學技術大學,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] 張書俊.固體火箭發動機粘彈性藥柱結構可靠性分析[J].固體火箭技術，2006，29(3).