旋轉條件下大長徑比固體火箭發(fā)動機三維內流場數值模擬

邱 欣，李新元，丁 彪，張春龍

（1.海軍航空工程學學院研究生管理大隊，山東煙臺264001；2.92635部隊，山東青島266041）

某小型固體火箭發(fā)動機為滿足導彈的射程，采用了大長徑比自由裝填式裝藥、旋轉發(fā)射方式，提高導彈飛行的穩(wěn)定性；通過采用2 種不同燃速的推進劑和雙側開槽半包覆裝藥實現2級推力。其裝藥結構如圖1 所示。在實現上述性能的同時，也增加了發(fā)動機內部流場的復雜程度，因而有必要研究發(fā)動機在工作過程中的流場特點及其對發(fā)動機正常工作可能產生的不利影響，為發(fā)動機的使用與改進提供參考。

圖1 發(fā)動機裝藥結構圖

發(fā)動機燃面計算是流場計算的重要步驟之一，目前較為常用的方法是直接根據燃面平行推移原理計算燃面的變化規(guī)律[1-4]。若考慮燃氣壓強分布對裝藥不同部位燃速影響，則計算結果更為精確，但實現起來較為復雜[5-6]。對于高過載條件下發(fā)動機內部流場的特征，國內外學者進行了大量的研究[7-9]。研究發(fā)現發(fā)動機的運動將對流場造成影響，容易導致燃氣對殼體燒蝕的不確定性。旋轉載荷同樣對發(fā)動機的內部工作過程有明顯影響[10-12]，武曉松等[13]通過大量實驗發(fā)現旋轉對增程火箭發(fā)動機推力影響明顯，根據不同旋轉速度下發(fā)動機的時間推力曲線，采用零維內彈道計算方法，獲得了噴管喉部等效面積和動態(tài)燃速隨轉速變化的經驗公式。王革等[14]采用有限體積法與RNGk-ε湍流模型進行計算，發(fā)現隨著發(fā)動機旋轉速度的增大，發(fā)動機燃氣切向速度峰值升高，燃燒室內的旋流與前封頭相互作用變得更為強烈。

本文采用有限體積法計算了某大長徑比固體火箭發(fā)動機不同燃去肉厚條件下的三維流場，考慮了軸向加速度與旋轉的影響，進行了裝藥燃去不同肉厚條件下的流場計算，獲得了該型發(fā)動機的內流場特征。

1 計算模型及方法

1.1 燃面計算

假設藥柱未包覆段燃面各處燃速相同。首先，以發(fā)動機長軸為軸，以助推段未包覆推進劑邊界為母線旋轉，與原始藥柱相減，可得到助推段未包覆裝藥燃去不同肉厚時的藥型，如圖2 a）所示。然后，以槽的初始輪廓為基礎，建立三維槽體，與原始藥柱相減，實現開槽部位燃面的推移。由于開槽段底部燃面為由球面、柱面和平面相切構成的復雜曲面，因而對開槽底部進行了圓弧導角過渡處理，如圖2 b）所示。

發(fā)動機裝藥內部沿軸向埋入6 根銀絲，如圖2 b）所示。假設銀絲的穩(wěn)定燃面近似為錐形燃面，則銀絲產生的燃面模擬可由一系列錐頂角一定的圓錐體與藥柱相減得到。錐面以燃速r推移，則錐高以速度而增加，θ為錐角。銀絲的燃速放大倍數為5，則得到錐角θ為23°。根據圖2 b）中點A、B、C的位置與間距可確定銀絲燃面出現的位置和時間，然后計算得到錐形燃面的深度。通過以上方法，可以得到藥柱燃去任意肉厚時刻下的外形，如圖3所示。

圖2 燃面推移實現方法

圖3 含銀絲錐形燃面的藥柱變化過程

1.2 網格劃分

由于燃面復雜，計算區(qū)域采用了四面體網格劃分。流場區(qū)域存在銀絲燃面形成的尖銳錐角，導致網格劃分出現負體積，解決辦法是將尖角處削平，使其尺寸與最小網格尺寸接近。為節(jié)約計算資源，并保證一定計算精度，在錐角上劃分較密的網格，然后在錐底劃分較粗的網格，形成從錐角至錐底由細到粗的過渡網格，如圖4所示。

圖4 計算區(qū)域網格劃分與邊界條件

1.3 求解方法

假設推進劑燃燒產物為理想氣體，采用有限體積法、標準k-ε兩方程湍流模型以及二階迎風格式離散方法對流場進行計算。選擇了三維單精度的可耦合求解能量和動量方程的隱式耦合解法器（Coupled Implicit），求解了可壓縮流的三維N-S 方程。由于發(fā)動機工作時以一定角速度旋轉，采用了Moving Reference Frame模型進行模擬。根據裝藥燃去不同肉厚時的燃面及彈體質量，計算發(fā)動機推力及其加速度，并加載到流場計算中。

1.4 邊界條件

入口條件設為質量入口邊界，為了與燃面推移中假設一致，假設同一時刻藥柱各處燃速相同，入口邊界的單位面積質量通量ψ可由ψ=r?ρp求得。式中：ψ為質量通量，單位為kg/(s ?m2)；r為推進劑燃速，單位為m/s，具體參數由內彈道計算求得；ρp為推進劑密度，取1 780 kg/m3；燃面溫度近似取推進劑絕熱燃燒溫度3 000 K。采用壓力出口邊界條件，初始給定出口的靜壓為一個標準大氣壓進行計算，當噴管出口速度為超聲速時，給定的出口壓力將不再起作用，此時出口壓力根據燃燒室內部流動外推得到。采用無滑移壁面條件，利用標準壁面函數法處理邊界湍流。

2 內流場數值模擬結果

經計算得到了燃去不同肉厚時刻的準靜態(tài)流場分布，噴管喉部馬赫數為1，符合實際情況。由于發(fā)動機三維流場分析結果數據量較大，流速矢量、流跡等結果不易于直接觀察分析，因而在分析結果的后處理過程中，建立了沿發(fā)動機徑向相互垂直的2個截面，如圖5 a）所示，可以通過分析截面上的流場參數，獲知三維流場情況。

2.1 典型截面上的流速、流跡分布特點

1）流速在不同部位的變化。第1個流速顯著變化的部位是裝藥懸臂段末端，如圖5 b）～d）所示。由于燃氣通道突然擴大，使燃氣得到膨脹，內能轉化為動能，流速在僅十幾cm 的燃氣通道內增大約5 倍，隨后由于周圍燃氣流速較低，使在該區(qū)域加速的燃氣流速又有所降低。

第2個流速變化的部位是裝藥開槽段與懸臂段相連部位，如圖5 b）～d）所示，是由于燃氣通道的面積變大而產生。該燃氣加速區(qū)域與上一加速區(qū)域正好相反，在助推段裝藥燃燒前期加速并不明顯，在助推段燃燒后期，該區(qū)燃氣以加速至之前的3 倍。在經歷短暫的加速后，其流速有所下降。

第3 個流速變化區(qū)域位于銀絲產生的錐形燃面處，如圖5 e）、f）所示。

圖5 分析截面以及不同燃去肉厚條件下XOY截面的最大流速分布

2）流速在槽內流動的特點。雖然裝藥兩側開槽部分燃氣通道在初始階段較為狹長，有發(fā)生侵蝕燃燒的可能性，但是在本文模型假設基礎上，經計算后發(fā)現開槽內部的燃氣流速較低，認為槽內發(fā)生侵蝕燃燒效應的可能性較低，燃氣僅在開槽出口處加速明顯。旋轉對發(fā)動機開槽出口處燃氣流速影響較為明顯，由圖6 中流跡分布所示，在XOY截面上流出的燃氣質點在流出槽口時，已經有較大角度的流跡偏移，因而在發(fā)動機整個工作過程中，開槽出口部位推進劑發(fā)生侵蝕燃燒的可能性要大于開槽內部推進劑。

2.2 燃燒室尾端渦流區(qū)的影響

流速沿燃燒室軸向變化較大，且呈現出流速沿軸向先升高后降低再升高的過程。尤其是在通道面積變化明顯的部位。由于裝藥懸臂段尾部產生了較大的流速梯度分布，加之燃氣的粘性在藥柱尾端噴管入口之前的區(qū)域形成了較長的渦流區(qū)，該渦流區(qū)內的流速較低，阻礙了后續(xù)燃燒產生燃氣的流出，該渦流區(qū)域在裝藥點燃時刻起產生，并一直到助推段燃燒后期逐漸消失，由圖6 a）～d）所示，這是該種裝藥流場與其他貼壁澆注式裝藥流場的最大不同之處，將對發(fā)動機的正常工作造成不利的影響，主要體現在以下幾個方面：

1）渦流區(qū)的存在，可能會造成噴管效率的下降，更嚴重的情況是造成噴管壅塞，燃燒室壓力過高。另一方面，因為該渦流區(qū)域易受裝藥的擺動和彈體側向加速度的影響，是動態(tài)變化的，有可能造成發(fā)動機工作的不穩(wěn)定性，當達到燃燒室的固有頻率后有可能產生振蕩燃燒，使燃燒室壓強瞬時升高到極限值。

2）由于有該渦流區(qū)域的存在，等效于減小了發(fā)動機的燃氣通道面積，在質量流量不變的情況下，則勢必增加了燃氣在渦流區(qū)兩側的流速，再加上發(fā)動機旋轉造成燃氣的離心流動，使該型發(fā)動機燃燒室尾段絕熱層的工作環(huán)境更為嚴酷。

圖6 不同燃去肉厚條件下的流速（左）、流跡（右）分布

2.3 懸臂段燃氣加速區(qū)的影響

在助推段裝藥燃燒后期，發(fā)動機燃燒室尾部的渦流區(qū)已經消失，但是另一種不利于發(fā)動機安全工作的流場形態(tài)隨之出現，由圖5 e）、f）可知。之前在裝藥尾端的燃氣加速區(qū)域已逐漸前移至懸臂段的中部，其流速約為開槽內部和燃燒室尾部燃氣流速的2 倍，而且在懸臂段燃盡前一直存在。在這個階段，導彈正處于推力最大，加速度最大的時刻，而此時懸臂段裝藥直徑也迅速下降，再加上發(fā)動機的旋轉可能造成的不對稱離心力，極有可能使懸臂段斷裂或因侵蝕燃燒而過早燃斷。若懸臂段斷裂則流至噴管僅需約10 ms，根據懸臂段裝藥初始直徑和推進劑燃速，斷裂殘藥很難順利垂直通過噴管喉部而不對燃燒室壓強產生影響。

3 結論

發(fā)動機工作過程中，燃氣的流速沿燃燒室軸向變化較大，開槽出口附近、懸臂段尾部及中部先后經歷了增速燃氣的沖刷，因而該部位殼體絕熱層工作條件較為嚴酷。為了防止發(fā)動機殼體燒穿，必須提高燃燒室尾部以及開槽部位絕熱層耐燒蝕性能，如采用性能更好的耐燒蝕材料或增加其厚度。

藥柱懸臂段隨著裝藥的燃燒逐漸變細，且由于發(fā)動機工作室以較高的速度旋轉，因而該懸臂段藥柱有斷裂的危險。為了防止懸臂段藥柱斷裂，應采用硬度和模量較高的固體推進劑，并提高藥柱加工精度，避免細小缺陷的出現，或者在保證發(fā)動機內彈道性能的前提下，適當減小懸臂段的長度。

[1] SAINTOUT E，RIBEREAU D，PERRIN P. A tool for 3d surface regression analysis in propellant grains. AIAA 1989-2782[R].AIAA，1989.

[2] FEI QIN. Algorithm study on burning surface calculation of solid rocket motor with complicated grain based on level set methods，AIAA 2006-4774[R].AIAA，2006.

[3] 于勝春，趙汝巖，周紅梅.基于Pro/E特征造型技術的固體發(fā)動機裝藥燃面計算[J]. 固體火箭技術，2005，28（2）：108-111.

YU SHENGCHUN，ZHAO RUYAN，ZHOU HONGMEI.Burning area regression simulation and calculation of solid rocket motor grains based on Pro/E feature modeling technology[J].Journal of Solid Rocket Technology，2005，28（2）：108-111.（in Chinese）

[4] 熊文波，劉宇，任軍學，等.基于單元法的三維裝藥通用燃面計算[J].航空學報，2009，30（7）：1176-1180.

XIONG WENBO，LIU YU，REN JUNXUE，et al. Generalized burning surface calculation of three dimensional propellant based on element method[J]. Acta Aeronautica Et Astronautica Sinica，2009，30（7）：1176-1180.（in Chinese）

[5] JUNG J CHOI，SURESH MENON. Simulation of composite solid propellant combustion with and without internal burning，AIAA 2010-617[R].AIAA，2010.

[6] JUNG J CHOI，SURESH MENON.Solid propellant combustion with surface regression using coupled gas-solid approach，AIAA 2009-4896[R].AIAA，2009.

[7] MINE YUMUSAK，FRANCOIS VUILLOT，TUGRUL H YINAZTEPE. Viscous internal flow applications for solid propellant rocket motors，AIAA 2006- 5116[R].AIAA，2006.

[8] 何國強，王國輝，蔡體敏，等.過載條件下固體發(fā)動機內流場數值模擬[J].推進技術，2002，23（3）：182-185.

HE GUOQIANG，WANG GUOHUI，CAI TIMIN，et al.Numerical simulation on 3-d two-phase flow field in SRM with acceleration load[J]. Journal of Propulsion Technology，2002，23（3）：182-185.（in Chinese）

[9] LANGHENRYM T，MARTIN M. Acceleration effects in solid propellant rocket motors，AIAA 1998-8621577[R].AIAA，1998.

[10] 張亞林. 高速旋轉固體火箭發(fā)動機內流場研究[D]. 南京：南京理工大學，2005.

ZHANG YALIN. Investigation of inner flow field in solid rocket motor with high speed rotation[D]. Nanjing：Nanjing University of Science and Technology，2005.（in Chinese）

[11]嚴聰，何國強，劉洋.旋轉條件下長尾噴管發(fā)動機三維兩相流場數值模擬[J].空軍工程大學：自然科學版，2009，10（5）：10-15.

YAN CONG，HE GUOQIANG，LIU YANG. Three dimension two phase flow field numerical simulation on tail-pipe nozzle srm under spinning condition[J]. Journal of Air Force Engineering University：Natural Science Edition，2009，10（5）：10-15.（in Chinese）

[12] 陳亮. 旋轉固體火箭發(fā)動機內流場數值研究[D]. 哈爾濱：哈爾濱工程大學，2007.

CHEN LIANG. Numerical investigation of flow field in spinning solid rocket motor[D].Harbin：Harbin Engineering University，2007.（in Chinese）

[13]武曉松，王棟，余陵，等.高速旋轉固體火箭發(fā)動機的動態(tài)燃速特性研究[J].彈道學報，2005，17（2）：1-7.

WU XIAOSONG，WANG DONG，YU LING，et al. The dynamic combustion rate characteristics of solid propellant rocket motor with high speed rotating[J]. Journal of Ballistics，2005，17（2）：1-7.（in Chinese）

[14]王革，陳亮，郜冶，等.旋轉對固體火箭發(fā)動機燃燒室燃氣流動的影響[J]. 空氣動力學學報，2008，26（2）：208-211.

WANG GE，CHEN LIANG，GAO YE，et al. Spinning effect on the flow of exhaust gas in combustion chamber of a solid rocket motor[J].ACTA Aerodynamic Sinica，2008,26（2）：208-211.（in Chinese）