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

（1.海軍航空工程學(xué)學(xué)院研究生管理大隊(duì)，山東煙臺(tái)264001；2.92635部隊(duì)，山東青島266041）

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

圖1 發(fā)動(dòng)機(jī)裝藥結(jié)構(gòu)圖

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

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

1 計(jì)算模型及方法

1.1 燃面計(jì)算

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

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

圖2 燃面推移實(shí)現(xiàn)方法

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

1.2 網(wǎng)格劃分

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

圖4 計(jì)算區(qū)域網(wǎng)格劃分與邊界條件

1.3 求解方法

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

1.4 邊界條件

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

2 內(nèi)流場(chǎng)數(shù)值模擬結(jié)果

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

2.1 典型截面上的流速、流跡分布特點(diǎn)

1）流速在不同部位的變化。第1個(gè)流速顯著變化的部位是裝藥懸臂段末端，如圖5 b）～d）所示。由于燃?xì)馔ǖ劳蝗粩U(kuò)大，使燃?xì)獾玫脚蛎?，?nèi)能轉(zhuǎn)化為動(dòng)能，流速在僅十幾cm 的燃?xì)馔ǖ纼?nèi)增大約5 倍，隨后由于周?chē)細(xì)饬魉佥^低，使在該區(qū)域加速的燃?xì)饬魉儆钟兴档汀?/p>

第2個(gè)流速變化的部位是裝藥開(kāi)槽段與懸臂段相連部位，如圖5 b）～d）所示，是由于燃?xì)馔ǖ赖拿娣e變大而產(chǎn)生。該燃?xì)饧铀賲^(qū)域與上一加速區(qū)域正好相反，在助推段裝藥燃燒前期加速并不明顯，在助推段燃燒后期，該區(qū)燃?xì)庖约铀僦林暗? 倍。在經(jīng)歷短暫的加速后，其流速有所下降。

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

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

2）流速在槽內(nèi)流動(dòng)的特點(diǎn)。雖然裝藥兩側(cè)開(kāi)槽部分燃?xì)馔ǖ涝诔跏茧A段較為狹長(zhǎng)，有發(fā)生侵蝕燃燒的可能性，但是在本文模型假設(shè)基礎(chǔ)上，經(jīng)計(jì)算后發(fā)現(xiàn)開(kāi)槽內(nèi)部的燃?xì)饬魉佥^低，認(rèn)為槽內(nèi)發(fā)生侵蝕燃燒效應(yīng)的可能性較低，燃?xì)鈨H在開(kāi)槽出口處加速明顯。旋轉(zhuǎn)對(duì)發(fā)動(dòng)機(jī)開(kāi)槽出口處燃?xì)饬魉儆绊戄^為明顯，由圖6 中流跡分布所示，在XOY截面上流出的燃?xì)赓|(zhì)點(diǎn)在流出槽口時(shí)，已經(jīng)有較大角度的流跡偏移，因而在發(fā)動(dòng)機(jī)整個(gè)工作過(guò)程中，開(kāi)槽出口部位推進(jìn)劑發(fā)生侵蝕燃燒的可能性要大于開(kāi)槽內(nèi)部推進(jìn)劑。

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

流速沿燃燒室軸向變化較大，且呈現(xiàn)出流速沿軸向先升高后降低再升高的過(guò)程。尤其是在通道面積變化明顯的部位。由于裝藥懸臂段尾部產(chǎn)生了較大的流速梯度分布，加之燃?xì)獾恼承栽谒幹捕藝姽苋肟谥暗膮^(qū)域形成了較長(zhǎng)的渦流區(qū)，該渦流區(qū)內(nèi)的流速較低，阻礙了后續(xù)燃燒產(chǎn)生燃?xì)獾牧鞒?，該渦流區(qū)域在裝藥點(diǎn)燃時(shí)刻起產(chǎn)生，并一直到助推段燃燒后期逐漸消失，由圖6 a）～d）所示，這是該種裝藥流場(chǎng)與其他貼壁澆注式裝藥流場(chǎng)的最大不同之處，將對(duì)發(fā)動(dòng)機(jī)的正常工作造成不利的影響，主要體現(xiàn)在以下幾個(gè)方面：

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

2）由于有該渦流區(qū)域的存在，等效于減小了發(fā)動(dòng)機(jī)的燃?xì)馔ǖ烂娣e，在質(zhì)量流量不變的情況下，則勢(shì)必增加了燃?xì)庠跍u流區(qū)兩側(cè)的流速，再加上發(fā)動(dòng)機(jī)旋轉(zhuǎn)造成燃?xì)獾碾x心流動(dòng)，使該型發(fā)動(dòng)機(jī)燃燒室尾段絕熱層的工作環(huán)境更為嚴(yán)酷。

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

2.3 懸臂段燃?xì)饧铀賲^(qū)的影響

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

3 結(jié)論

發(fā)動(dòng)機(jī)工作過(guò)程中，燃?xì)獾牧魉傺厝紵逸S向變化較大，開(kāi)槽出口附近、懸臂段尾部及中部先后經(jīng)歷了增速燃?xì)獾臎_刷，因而該部位殼體絕熱層工作條件較為嚴(yán)酷。為了防止發(fā)動(dòng)機(jī)殼體燒穿，必須提高燃燒室尾部以及開(kāi)槽部位絕熱層耐燒蝕性能，如采用性能更好的耐燒蝕材料或增加其厚度。

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

[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特征造型技術(shù)的固體發(fā)動(dòng)機(jī)裝藥燃面計(jì)算[J]. 固體火箭技術(shù)，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] 熊文波，劉宇，任軍學(xué)，等.基于單元法的三維裝藥通用燃面計(jì)算[J].航空學(xué)報(bào)，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，F(xiàn)RANCOIS VUILLOT，TUGRUL H YINAZTEPE. Viscous internal flow applications for solid propellant rocket motors，AIAA 2006- 5116[R].AIAA，2006.

[8] 何國(guó)強(qiáng)，王國(guó)輝，蔡體敏，等.過(guò)載條件下固體發(fā)動(dòng)機(jī)內(nèi)流場(chǎng)數(shù)值模擬[J].推進(jìn)技術(shù)，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] 張亞林. 高速旋轉(zhuǎn)固體火箭發(fā)動(dòng)機(jī)內(nèi)流場(chǎng)研究[D]. 南京：南京理工大學(xué)，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]嚴(yán)聰，何國(guó)強(qiáng)，劉洋.旋轉(zhuǎn)條件下長(zhǎng)尾噴管發(fā)動(dòng)機(jī)三維兩相流場(chǎng)數(shù)值模擬[J].空軍工程大學(xué)：自然科學(xué)版，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] 陳亮. 旋轉(zhuǎn)固體火箭發(fā)動(dòng)機(jī)內(nèi)流場(chǎng)數(shù)值研究[D]. 哈爾濱：哈爾濱工程大學(xué)，2007.

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

[13]武曉松，王棟，余陵，等.高速旋轉(zhuǎn)固體火箭發(fā)動(dòng)機(jī)的動(dòng)態(tài)燃速特性研究[J].彈道學(xué)報(bào)，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]王革，陳亮，郜冶，等.旋轉(zhuǎn)對(duì)固體火箭發(fā)動(dòng)機(jī)燃燒室燃?xì)饬鲃?dòng)的影響[J]. 空氣動(dòng)力學(xué)學(xué)報(bào)，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）