水下固體火箭發動機推力脈動特征研究

唐云龍，李世鵬，劉　筑，隋　欣，王寧飛

(北京理工大學 宇航學院，北京　100081)

?

水下固體火箭發動機推力脈動特征研究

唐云龍，李世鵬，劉筑，隋欣，王寧飛

(北京理工大學 宇航學院，北京100081)

對水下固體火箭發動機的動態特性進行了研究，根據動網格及VOF(Volume of Fluid)基本原理，構建了用于水下固體火箭發動機動態射流響應的軸對稱模型；根據壁面積分基本原理及被動方式求解水下工作發動機推力的計算方法，成功模擬了發動機的受力運動；對不同工況下發動機工作情況進行仿真，研究了水下發動機動態模型下的射流結構，發現在來流作用下，射流氣體會受到來流邊界層的擠壓，產生頸縮現象，之后逐漸向外膨脹；研究結果顯示，在來流邊界層作用下，射流所產生的激波在向下游發展之前，將經歷多次折射，該研究結果與文獻中相關實驗結果較吻合；對發動機不同工況下的推力脈沖峰值進行分析，找出了動態模型下發動機點火初期推力脈沖峰值特征隨工作環境的變化規律。最后，對有/無軸向浮力時的發動機動態特征進行了對比分析。

動網格；水下射流；邊界層；推力峰值

0　引言

利用動網格方法，研究水下固體火箭發動機運動過程的射流結構及推力問題，與水下工作發動機的真實情況更加接近，且有利于掌握水下火箭的工作規律，能為水下發動機設計及工程應用提供幫助。國內外在水下射流流場結構與流動特性研究領域發展均較緩慢[1]。唐嘉寧等[2]將球形氣泡模型與CFD計算結果進行比較，發現球形氣泡模型結構簡單、計算速度快,并能有效反映發動機點火瞬間的推力峰值特性。但該模型基于氣體在整個射流過程中以單一氣泡結構存在，并沿徑向等速無限膨脹假設，并未考慮水-氣間相互作用，因此對推力的解釋并不是很全面。Lindau等[3]對考慮空泡作用下高速超空化水下推進過程進行了數值研究，提出了燃氣、水和水蒸氣的三相流計算模型，但未對水下推進過程中的流動特性及發動機工作性能做進一步分析。唐嘉寧等[4]建立了水下超音速氣體射流軸對稱數學模型，并將數值計算結果與實驗結果進行對比，驗證了該模型的合理性和準確性[5]。

為研究固體火箭發動機在水下工作時水的壓差阻力和來流速度等對其穩態推力的影響，本文建立了發動機水下流場模型,并在發動機前端加設彈體段，通過數值方法，分析了固體火箭發動機水下運動達到平衡速度后的流場結構及發射初期的推力振蕩峰問題。研究結果對導彈水下發射的流體動力特性分析及發動機設計具有指導意義。

1　基本模型及方程

1.1VOF原理

VOF是通過研究網格單元中流體和網格體積比函數F對自由面進行確定，并追蹤流體的變化，并非追蹤自由液面上質點的運動。VOF方法可處理自由面重入等強非線性現象,其有計算時間短、存儲量少等特點，但在處理F變化時稍顯繁瑣，其有一定人為因素。若F=1，則該單元為指定相流體單元；若F=0，則該單元為無指定相流體單元；當0

1.2動網格技術

動網格技術方便研究彈體與水間相對運動，又能保證網格塊之間流動變量的插值精度。在使用動網格技術時，一般將流動域劃分為2個或2個以上的計算域，相鄰的計算域由1個或多個交界面連接，并沿交界面滑動。計算過程中，運動域沿著交界面滑動，每一時間步后，重新計算網格位置，確定交界面上網格點的相對位置，并在交界面處進行通量傳遞[7]。

1.3層變網格方法

在四邊形網格區域內，動態層變方法可根據與運動相鄰的網格層高度增加或減少網格層數。在計算中，動網格模型可指定一理想的高度。如圖1所示，鄰近邊界的網格單元層根據層j的單元層高度h來分裂或與鄰層合并成新層。若層j中單元體積處于膨脹狀態，其膨脹極限如下[8]：

(1)

式中hideal為理想單元高度；ah為全局單元層分裂或合并因子。

當h>(1+ah)hideal時，層j中的單元將分裂成高度hideal的單元層和高度為h-hideal的單元層。若層j中單元體積處于被壓縮狀態，則當壓縮到h

圖1　動態網格層

1.4動網格控制方程

φ代表通用標量，靜止區域和運動區域內的方程可寫為如下形式積分方程[9]：

(2)

式中ρ為密度；ug為網格運動速度；Γ為湍流擴散率；Sφ為源項；φ為運動區域內的任意控制體；?V代表控制體的邊界。

方程(2)的時間導數項可利用一階向后差分寫成如下形式：

(3)

其中，n及n+1表示時間。第n+1時間的層網格體積可通過式(4)計算：

(4)

其中，dV/dt為控制體的體積時間導數。為滿足網格守恒率，控制體的體積時間導數通過式(5)進行計算：

(5)

(6)

其中，?Vj為整個時間步Δt上因控制容積面j膨脹而引起的體積變化。

2　計算方法

2.1推力

在計算過程中，由于動網格邊界運動規律是未知的，故利用被動方式來求解水下工作發動機的推力。其中，Fi和Fe分別表示作用在發動機外表面和內表面上壓強的合力(如圖2所示)，則推力表達式如下[10]：

(6)

圖2　發動機計算方法

2.2速度

本文研究的水下發動機涉及被迫運動問題，其受力非常復雜。研究過程中，利用軸對稱模型，即發動機只僅在軸線方向運動。根據牛頓定理可知軸線上力F可表示為：

F=ma

(7)

故某一時刻速度表達式為：

(8)

式中v(NΔt)為NΔt時刻的速度；Δt為計算時間步長；Fk為kΔt時刻所對應的推力；M為發動機總質量；Δm為推進劑在單個時間步長內的質量損失。

3　仿真模型

在結合水下發動機的實際工作情況的基礎上[9]，建立了相應的二維軸對稱模型。計算域如圖3(a)所示，包含噴管內流區域和外部射流區域。發動機噴管內流區域及鄰近發動機外流區域使用非結構網格作為過渡。在遠離發動機區域使用結構網格作為滑動區域。

(a)計算網格

(b)發動機周邊網格　　　　(c)前端不動網格

圖3(a)中，上下的滑動網格與靜態網格采用接觸邊界進行接觸。在運動區域的起始與結束段采用內部網格作為動態和靜態網格的分界(如圖3(c)所示)。整個區域的外邊界均為壓力出口，而發動機噴管入口則采用壓力入口作為邊界條件。因為采用了滑動網格技術，網格的產生、消失及變形都會對發動機進流場產生較大的影響。所以，為了減少動網格本身對流場的影響，動網格與靜網格的接觸區域選在流場遠端，且滿足長時間非穩態計算的要求。

4　仿真結果及分析

4.1流場

4.1.1流場結構特征分析

圖4為利用圖3中的水下發動機單自由度運動模型達到平衡速度時的計算結果及發動機水洞實驗結果。兩者都是采用燃燒室內壓恒定的射流模型，且在來流作用下，軸線推力達到平衡值，都趨近于零。

(a)水下發動機射流體積分數仿真結果

(b)水下發動機射流體積分數實驗結果[11]

(c)水下發動機射流密度分布仿真結果

(d)水下發動機射流激波分布仿真結果

從來流與水下發動機射流耦合作用下的射流及激波結構上分析，可知圖4 (a)、(b) 中射流仿真的體積分數分布與文獻[11]中仿真結果(如圖4(c)所示)及文獻[12]中實驗結果較一致。由于受到來流慣性作用，噴管出口處會產生邊界層，射流形成的激波(如圖4(d))在邊界層內不斷折射導致通道變窄(如圖中A、B、E、F所示)氣體膨脹受阻隨著射流的方向發展。來流對射流的作用減弱，射流開始向外膨脹(如圖中C、D、G、H所示)，其中A、B、E、F為燃氣射流起始段(如圖(a)中L1所示)，C、D、G、H之后為燃氣射流主體段(如圖(a)中L3所示)，在兩者之間的為過渡段(如圖(a)中L2所示)。通過模擬水下射流流場的運動過程，發現在有來流情況下，由于外界水流場作用，燃氣射流不再斷裂，并未發生回擊現象，與文獻[11]中結論一致，證明該模型能真實反映動態過程。

圖5為射流流動特性與軸線上壓強和速度無因次化后的對應關系。按圖4(a)中的方法對軸線上的壓強和速度進行劃分，研究發現起始段、過段段和主體段特征較明顯，頸縮對射流進行了二次加速。

4.1.2水下發動機運動平衡后軸線特征

何淼生[13]、王春旭[14]、曹嘉怡[15]等學者已經對水下發動機的開機、關機及出水過程進行了詳盡研究。所以，本文將不再對其進行重復性研究，只對在水下發動機運動平衡后的特征進行非穩態分析。

圖6為水下發動機在達到平衡速度后，在不同時刻發動機近流場的分布特征。如圖6(a)、(b)所示，水下發動機射流流場在邊界層與激波共同作用下而產生的頸縮現象。

圖5　某一時刻射流流動特性與軸線上壓強和速度無因次化后的對應關系

(a)壓強分布

(b)馬赫數分布圖

從圖6(a)、(b)可知，當發動機速度達到平衡后，尾流場結構及發動機前段流場均無明顯改變，射流起始段長度基本固定。采用圖5中方法結合流場變化，可描述不同工況下發動機尾流場中頸縮點與射流起始段長度，因為衡后射流結構及發動機軸線特征不再變化，所以不再對發動機速度平衡后時刻的流場分布進行特殊說明。

4.2推力特征分析

圖7(a)為固體火箭發動機推力-時間曲線仿真結果，如圖7所示，其動態模型的推力無大幅振蕩，與靜態模型有明顯差別。在射流初期，射流噴出噴管形成氣泡，并產生強勁的推力脈沖，本文仿真模型與文獻[13]中實驗模型均采用端燃藥柱，認為燃燒室內壓不變，且對徑向位移添加了相應約束，只在軸向上單一自由度運動。

(a)航行器推力-時間曲線

(b)航行器粘性阻力-時間曲線

(c)航行器速度-時間曲線

(d)航行器位移時間曲線

通過對水下發動機軸向推力的監測數據比較，發現本文的仿真結果與文獻[13]中實驗結果一致；隨著速度和推力的增加，發動機在壓差阻力作用下，開始減小直至到零，并在零值附近做小幅值振蕩。如圖7(b)所示，發動機的黏性阻力相對其壓差阻力可忽略不計。如圖7(c)所示，發動機運動速度在推力脈沖的作用下不斷增加，當來流阻力與發動機推力平衡時，速度也趨于平衡。如圖7(d)所示，通過速度對時間積分，可得發動機位移變化曲線，當發動機達到平衡速度后，位移曲線基本呈線性改變。

為了研究水下固體發動機單自由度運動模型在不同工況下發射初期推力的振蕩問題，本文將不同總壓、總溫、水深、擴張比及質量按照單調變化趨勢分成5種工況進行比較(見圖8)。從圖8(a)可知，水下固體火箭發動機在點火初期會出現一個非常大的峰值，該峰值對水下發動機的運動特征具有很大影響。

從圖8(a)～(e)可知，在計算工況下，總壓越大，推力脈沖峰值越大；水深越大，推力脈沖峰值越小；擴張比越小，推力脈沖峰值越大；總溫越高，推力脈沖峰值越大；質量越大，發動機推力脈沖峰值越大，且推力脈沖與發動機推力積分及平衡速度之間具有相關性。

不同質量下的速度-時間曲線見圖9。由圖9可知，水下發動機質量變化不會改變其運動的平衡速度，但在運動相同時間下，水下固體發動機的位移卻隨質量的增加而不斷減小，這說明水下發動機運動速度的爬升過程不同。速度爬升段的面積增加與質量減少(如圖10所示)基本呈線性關系，且是否考慮浮力對其增量變化率基本沒有影響。

(a)燃燒室壓強對發動推力峰值的影響　　　　　　　(b)水深對發動機推力峰值的影響

(e)擴張比對推力峰值的影響

圖9　不同質量下的速度-時間曲線　

圖10　質量減少量與位移增量的關系

5　結論

對發動機在不同工況下的運動過程進行了仿真研究，發現水下發動機在運動過程中來流會與發動機尾部燃氣射流發生耦合；在發動機運動模型中，尾流不會產生斷裂/回擊；來流速度與射流激波的相互作用會在發動機尾部形成邊界層，使射流形成頸縮，且射流進行二次加速；根據射流特征對燃氣射流的起始段、過渡段和主段進行劃分，并與軸線上的速度和壓強特征進行對應；通過對不同時刻軸線及速度特性分析，發現當發動機達到平衡速度后，其運動規律及射流特性基本保持不變；動態模型發動機推力在初期會形成一個脈沖，且粘性阻力相對推力可忽略不計；通過對不同工況下發動機推力特征分析，得出以下結論：水深越大，發動機推力脈沖越小；燃燒室壓強越大，發動機推力脈沖越大；溫度越大，發動機推力脈沖越小；擴張比越大，發動機推力脈沖越大。改變發動機的質量不會改變發動機的平衡速度，但發動機質量越小，加速越快，且發動機軸向浮力對其速度增量影響很小。

[1]甘曉松，賈有軍.水下燃氣射流流場數值研究[J].固體火箭技術，2009，31(1):23-26.

[2]Tang J N,Li S P,Wang N F,et a1.Flow structures of gaseous jet injected into liquid for underwater propulsion[C]//46th AIAA/ASME/ SAE/ASEE Joint Propulsion Conference and Exhibit,Nashville,IN.2010,AIAA 2010-6911.

[3]Lindau J W,Venkateswaran S,Kunz R F,et a1.Multiphase computations for underwater propulsive flows[C]//16th AIAA Computational Fluid Dynamic Conference,Orland,Florida,AIAA 2003-4105.

[4]Tang Jia-ning,Wang Ning-fei,Wei Shyy.Flow structures of gaseous jets injected into water for underwater propulsion[J].Acta Mechanica Sinica,2011,27(4):461-472.

[5]唐嘉寧，李世鵬，王寧飛.水下固體火箭發動機的負推力現象研究[J].固體火箭技術，2012，35(3)：325-344.

[6]康玉東，孫冰．燃氣非平衡流再生冷卻流動傳熱數值模擬[J].推進技術，2011，32(1):119-124.

[7]粱偉,宋文萍.基于滑動網格的螺旋槳復雜繞流數值模擬研究[J].航空計算技術,2011,42(3):45-49.

[8]楊猛.計算流體力學兩相流流動的模擬及兩相流模型的研究[D].天津大學,2005.

[9]楊曉光,陳煥龍,劉華坪,等.導彈水下運動及出水過程的三維流場仿真[J].彈道學報,2010,22(1):107-110.

[10]唐嘉寧,劉向陽.水下固體火箭發動機推力特性研究[J].導彈與航行運載技術,2012(5):15-21.

[11]何曉.噴氣推進導彈空泡流的尺度效應研究[D].上海交通大學,2009.

[12]Ivan N Kirschner,Neal E Fine,James S Uhlman,et a1.Supercavitation research and development[J].Undersea Defense Technologies,Waikiki,HI,2001:1-10.

[13]何淼生,覃粒子,何佳磊,等.水下超聲速噴嘴起動/關機過程射流的動態不穩定性研究[J].推進技術,2014 ,35(4):523-529.

[14]王春旭,張濤,侯國祥.基于渦聲理論的水下射流噪聲預報[J].船舶力學,2010，14(6):670-677.

[15]曹嘉怡,魯傳敬,李杰,等.潛射導彈水下垂直自拋發射過程研究[J].水動力學研究與進展,2006,21(6):752-759.

(編輯：崔賢彬)

Research on thrust fluctuation characteristics of underwater solid rocket motor

TANG Yun-long，LI Shi-peng，LIU Zhu，SUI Xin，WANG Ning-fei

(School of Aerospace Engineering,Beijing Institute of Technology,Beijing100081,China)

The dynamic property of underwater solid rocket motor was studied in this article.Utilizing dynamic grid and the VOF principle,an axisymmetric model describing the dynamic response of the fluidic thrust of underwater solid rocket was established.Movement of the motor was successfully simulated based on the principle of surface integral and the computational method of underwater motor's thrust.By simulating motors under different cases using the dynamic model,the jet structure of underwater motor was studied.It turn out that the jet gas is extruded by the flow in the boundary layer,and generated necking phenomenon,and then expanded outward gradually.The results indicate that the shock wave produced by the jet will refract several times before going downstream,which is consistent with experimental results of other literature.The thrust's pulse peaks in different cases,were analyzed,and the change law of the features of underwater motors with the variation of environment was obtained.Finally,the dynamic characteristics of motors in conditions with/without buoyancy were compared and analyzed.

dynamic grid;water jet;boundary layer;thrust peak

2015-06-04；

2016-01-18。

國家自然科學基金(11072032)。

唐云龍(1986—)，男，博士生，從事固體火箭發動機仿真及測試研究。E-mail：ks669@126.com

V438

A

1006-2793(2016)04-0476-06

10.7673/j.issn.1006-2793.2016.04.005