噴管兩相流耦合數(shù)值模擬與粒子熱增量預測

蘇鵬輝 艾邦成 潘宏祿

(中國航天空氣動力技術(shù)研究院， 北京100074)

現(xiàn)代固體火箭發(fā)動機為了抑制不穩(wěn)定燃燒，提高工作效能，常使用含鋁粉的復合推進劑，其燃燒產(chǎn)物中含有氧化鋁凝聚相粒子，質(zhì)量分數(shù)有時可高達30%，形成典型兩相流動［1］.兩相流中粒子對噴管熱環(huán)境的影響主要體現(xiàn)在兩個方面:一方面粒子在流場中存在速度和溫度滯后，通過與氣相間的動量和能量交換引起氣相流場變化;另一方面高速運動的高溫粒子撞擊噴管壁面產(chǎn)生壁面熱流增加和粒子沖刷燒蝕等現(xiàn)象，對噴管喉襯熱防護提出了更高要求.因此基于燃氣-粒子相互作用的噴管兩相流研究已成為噴管喉襯熱防護研究的重要方面.

噴管兩相流研究主要有數(shù)值模擬和試驗測量這兩種手段.在數(shù)值模擬方面，目前國內(nèi)外學者主要基于歐拉-拉格朗日描述的粒子軌道模型對不同狀態(tài)下噴管兩相流的流場參數(shù)和粒子運動狀態(tài)進行數(shù)值計算分析.文獻［2－4］利用顆粒軌道模型對噴管兩相流進行了數(shù)值模擬研究，得到了豐富的成果，顯示了顆粒軌道模型的應用價值.

在噴管兩相流試驗，特別是噴管粒子熱增量試驗方面，文獻［5］進行了不同攻角平板的粒子熱增量試驗，粒子垂直平板速度所對應動能的70%轉(zhuǎn)化為平板表面熱能.文獻 ［6］進行了超聲速球頭繞流粒子熱增量試驗，在球頭駐點處粒子動能70%轉(zhuǎn)化為球頭駐點熱能.

目前已有的噴管兩相流數(shù)值模擬集中關(guān)注于粒子運動以及流場參數(shù)變化，而粒子熱增量研究則主要集中于試驗方面.通過噴管內(nèi)流場兩相流數(shù)值計算獲得壁面粒子熱增量分布并不多見.本文基于顆粒軌道模型，以含有氧化鋁粒子的發(fā)動機噴管兩相流為研究對象，通過有限體積方法和單元內(nèi)顆粒源 (PSIC，Particle Source in Cell)兩相耦合方法，模擬了不同直徑粒子對噴管流場及噴管壁面熱環(huán)境的影響.

1 控制方程及邊界條件

1.1 氣相控制方程

兩相流中氣相控制方程組采用守恒形式完全氣體三維N-S方程組

其中，F(xiàn)，G和H為對流項通量;Fv，Gv和Hv為粘性及傳熱耗散項通量.為封閉控制方程組需補充理想氣體狀態(tài)方程p=ρRT，氣體動力粘性系數(shù)μ通過薩特蘭公式得到.右端源項S表示粒子對氣相作用，其表達式在下文給出.

1.2 粒子相控制方程

氧化鋁粒子的熔點2345K，沸點3250K［7］，當燃燒室溫度在3000～3100K時，燃燒室內(nèi)及噴管喉部前端的氧化鋁粒子在表面張力作用下形成球形液滴［8］.使用粒子軌道模型對粒子相建模，粒子相與氣體相方程獨立求解，利用源項耦合，通過跟蹤計算大量粒子的實際運動軌跡及粒子狀態(tài)隨軌跡變化，實現(xiàn)對兩相流物理過程的模擬［9］.定常流動中，每個粒子軌道均代表了在流場入口截面上這個粒子所占據(jù)的面積微元內(nèi)所有流經(jīng)此微元的粒子軌跡及狀態(tài)變化.此方法節(jié)省存儲空間，易于并行計算，易于模擬粒子與氣相之間的速度和溫度滑移及粒子相復雜變化過程.

針對固體火箭發(fā)動機噴管兩相流特點，對粒子相做出假設(shè):①粒子均為一定直徑的球形;②粒子與氣相之間無質(zhì)量交換;③粒子與氣相之間存在對流換熱;④粒子與氣相之間只考慮阻力;⑤粒子與粒子之間沒有相互作用;⑥粒子不考慮湍流擴散效應.粒子速度、位移及溫度隨時間變化的拉格朗日型控制方程組如下:

其中，up，xp，Tp，Ap，mp和 cp分別為粒子速度、位置、溫度、表面積、質(zhì)量和熱容;u和T為流體速度和溫度，其余系數(shù)計算可參考文獻［10］.

對有限體積法中每個網(wǎng)格單元，統(tǒng)計穿過此單元的全部粒子在單元內(nèi)交換的動量和熱量，除以該控制體的體積Ω，即得到此控制體內(nèi)粒子源項［11］，S的表達式如下:

2 數(shù)值方法及定解條件

2.1 數(shù)值方法

氣相流場控制方程采用有限體積方法求解，其中對流項通量使用Roe格式求解，耗散項通量采用高斯方法求解，時間離散采用LU-SGS(Lower-Upper Symmeffic Gauss-Seidel)隱式方法.粒子相控制方程采用歐拉顯式格式求解.

兩相耦合過程采用PSIC方法求解［12－13］，過程是:首先獲取氣相場定常解，之后求解粒子軌道方程，統(tǒng)計粒子與氣相之間的動量及能量交換，作為源相加入氣相流場，重新計算氣相流場定常解，反復迭代，直至收斂到預期精度.迭代過程中為減小振蕩，加快收斂速度，在更新粒子源項時對粒子源項進行如下低松弛處理:

其中，Sn為第n次迭代中氣相控制方程的粒子源項;S*為求解粒子軌道方程得到的新源項;Sn+1為第n+1次迭代中氣相控制方程的粒子源項;α為松弛因子，在計算初期取0.5左右.

2.2 定解條件

2.2.1 粒子相入口條件

粒子在噴管入口加入流場，初始速度和初始溫度與當?shù)貧庀鄥?shù)相同，粒子質(zhì)量流量在入口截面上滿足均勻分布，每條粒子軌道所代表粒子質(zhì)量流量為mi=m/Nk，其中m為入口截面上總的粒子質(zhì)量流量，Nk為粒子軌道總數(shù).

2.2.2 粒子相壁面條件

有關(guān)試驗研究［5－6］表明，粒子撞擊壁面時，粒子垂直壁面速度所對應動能的70%轉(zhuǎn)化為壁面熱能，設(shè)粒子撞擊前平行和垂直壁面的速度分別為Vp1、Vn1，撞擊之后反彈速度分別為 Vp2、Vn2，不計粒子內(nèi)能變化，根據(jù)能量守恒有:

粒子撞擊壁面時，粒子攜帶的能量部分轉(zhuǎn)化為壁面熱量，根據(jù)上述試驗結(jié)論，每條粒子軌道在與壁面碰撞時，產(chǎn)生的壁面熱增量為

對全部顆粒軌道的壁面熱增量進行累加，可得到總的壁面熱增量.

3 結(jié)果與分析

選取典型噴管外形和燃燒室工況進行噴管兩相流的計算.圖1是噴管三維網(wǎng)格示意圖，左側(cè)為入口截面，右側(cè)為出口截面.計算所用參數(shù):燃燒室總壓5MPa，總溫3000 K，噴管壁面為等溫壁，壁溫300K.粒子由入口處加入流場，粒子質(zhì)量流量253kg/(m2·s)，粒子熱容為cp=880 J/(kg·K)，粒子物質(zhì)密度為ρp=3960kg/m3.

圖1 噴管計算網(wǎng)格示意圖

3.1 粒子軌道數(shù)量對熱增量均勻性影響

粒子軌道模型中需要模擬的軌道數(shù)量應根據(jù)具體問題的特點而確定，為了確定適合本文算例的軌道數(shù)量，進行了軌道數(shù)量與粒子熱增量均勻性的分析.由于該噴管幾何形狀以及流場邊界條件的軸對稱性，在理想情況下同一軸向位置處的周向各面粒子熱增量應保持一致，但是離散軌道的不均勻性會導致計算結(jié)果軸對稱性下降.模擬軌道數(shù)量越多，計算結(jié)果的均勻性越好，但計算所消耗機時也會成比例增加，因此通過對不同軌道數(shù)的模擬結(jié)果進行熱增量均勻性分析，可以確定本問題需要的軌道數(shù)量，在數(shù)據(jù)均勻性和計算時間之間取得平衡.

圖2粒子軌道數(shù)量N分別為1萬，10萬，40萬和100萬時在每個軸向位置x處的熱增量計算值沿周向分布的標準差σ.

圖2 不同軌道數(shù)量熱增量分布標準差

從圖2可以看到:當N＜40萬時，數(shù)據(jù)均勻性較差，隨著N增加，數(shù)據(jù)均勻性有明顯改善;當N＞40萬時，數(shù)據(jù)均勻性隨N增加改善較小.因此權(quán)衡數(shù)據(jù)均勻性和所需計算時間，本研究中軌道數(shù)量取為40萬條.

3.2 粒子軌道分析

為了考察粒子軌道模型對噴管內(nèi)燃氣-粒子兩相流問題求解的有效性，以及觀察粒子在噴管兩相流中的流動特點，對不同直徑粒子的軌道進行了計算和分析.粒子直徑?jīng)Q定了粒子質(zhì)量和受力，隨著粒子直徑增大，粒子在流場中跟隨性變差，分布范圍逐漸向軸線靠攏，在擴張段壁面附近出現(xiàn)逐漸增大的粒子稀疏區(qū)域.

粒子直徑較小時慣性也較小，在流場中的跟隨性很好，圖3a是直徑1 μm粒子的軌道分布，可以觀察到粒子軌跡在噴管中分布均勻，在噴管收縮段內(nèi)粒子軌跡隨噴管型線轉(zhuǎn)彎，向噴管中心軸線靠攏，而在噴管擴張段內(nèi)粒子軌跡隨噴管型線一起擴展至整個噴管出口截面，粒子在收縮段內(nèi)與噴管喉襯表面沒有明顯碰撞，在擴張段內(nèi)亦無明顯的粒子稀疏區(qū)域.圖3b～圖3d是直徑為10～40 μm粒子的軌跡分布，可以觀察到，隨著粒子直徑增加，粒子慣性愈加明顯，在噴管收縮段內(nèi)粒子由于無法及時隨噴管型線轉(zhuǎn)向而與噴管壁面發(fā)生碰撞.在噴管擴張段，粒子軌道集中分布于軸線周圍的圓形區(qū)域內(nèi)，在噴管壁面附近產(chǎn)生了粒子稀疏區(qū)域.同一軸向位置處，小直徑粒子的分布區(qū)域半徑較大，粒子稀疏區(qū)域較小，而大直徑粒子的分布區(qū)域半徑較小，粒子稀疏區(qū)域較大，大直徑粒子的軌道更加趨向于噴管軸心.由以上結(jié)果可知，使用粒子軌道模型可以得到粒子在兩相流中的典型分布特征，粒子軌道模型對噴管兩相流問題是適用的.

圖3 不同直徑粒子軌道分布

3.3 兩相流流場分析

含有粒子的兩相流流場與純氣相流場相比有較大區(qū)別，為了研究粒子對流場的影響，需要采用兩相耦合方法對流場進行數(shù)值計算.

圖4給出了粒子直徑為20 μm時兩相流與純氣相流流場參數(shù)分布對比.圖4a～圖4c中，上半部分為兩相流流場參數(shù)分布，下半部分為純氣相流場參數(shù)分布.可以看到，兩相流中，氣相速度和溫度受粒子影響較大，而壓力分布受粒子影響較小，在噴管擴張段內(nèi)兩相流中氣相速度較純氣相流動時降低，兩相流中氣相溫度較純氣相流動時升高，這是由于擴張段內(nèi)氣相速度增加溫度降低，而粒子速度和溫度都存在滯后，因此粒子對氣相產(chǎn)生減速加溫作用.圖4a和圖4b中，噴管擴張段軸線附近，與純氣相流動相比，兩相流的速度和溫度等值線出現(xiàn)朝向噴管出口方向的偏移，這是由于在軸線附近位置粒子軌道分布較為密集，氣相受到粒子的減速加溫作用更加明顯.由圖4d～圖4f可見，在擴張段內(nèi)相同軸向位置處，兩相流比純氣相流速度更低、溫度更高.由上述結(jié)果可知，兩相耦合計算可以較為全面地反映流場中粒子與氣體之間的相互作用情況，所得結(jié)果更符合真實情況.

圖4 兩相流與純氣相流流場參數(shù)對比

3.4 粒子相壁面熱增量分析

為了預測粒子撞擊噴管壁面時動能轉(zhuǎn)化為熱能所產(chǎn)生的壁面熱增量，在兩相流計算中對粒子軌道與壁面碰撞情況進行統(tǒng)計分析，得到了噴管壁面熱增量的分布范圍和熱流大小.

圖5給出了對應于不同粒子直徑的壁面熱增量分布.可以觀察到，粒子熱增量集中分布于喉部前端收縮段內(nèi)，收縮段的前后粒子熱增量均迅速減小，噴管擴張段內(nèi)粒子熱增量并不明顯.粒子熱增量的分布范圍隨粒子直徑增大略有增加，而粒子熱增量峰值則隨粒子直徑的增大迅速增加，這是由于小直徑粒子隨流性很好，不會與壁面發(fā)生劇烈撞擊，因此其熱增量相對較小，而大直徑粒子慣性較大，在噴管收縮段內(nèi)無法及時隨流場轉(zhuǎn)向，因此撞擊壁面的現(xiàn)象比小直徑粒子更加明顯，壁面熱增量也就隨之增加.

圖5 不同直徑粒子熱增量分布

4 結(jié)論

本文采用基于粒子軌道模型和有限體積方法的PSIC兩相耦合算法，結(jié)合基于粒子熱增量試驗的粒子壁面反彈模型，對含有氧化鋁粒子的固體火箭發(fā)動機噴管兩相流及其對噴管壁面熱環(huán)境的影響進行了研究，得到以下結(jié)論:①使用粒子軌道模型得到的粒子軌道分布和流場參數(shù)分布能夠反映出噴管兩相流區(qū)別于純氣相流動的特點，說明粒子軌道模型對噴管內(nèi)燃氣-粒子流動是適用的;②兩相流耦合計算所得到的氣相溫度和速度分布與純氣相流動相比有較大區(qū)別，耦合計算所得結(jié)果反映了氣相和顆粒相之間的相互作用，對實際流動現(xiàn)象模擬更加準確;③粒子動能向熱能轉(zhuǎn)化產(chǎn)生的粒子熱增量集中分布于喉部前端收縮段內(nèi)，收縮段前后熱增量迅速減小，噴管擴張段內(nèi)粒子熱增量并不明顯，熱增量峰值隨粒子直徑增加明顯增大.

References)

［1］張宏安，葉定友，侯曉.固體火箭發(fā)動機凝聚相微粒分布研究現(xiàn)狀［J］.固體火箭技術(shù)，2000，23(3):25－28

Zhang Hongan，Ye Dingyou，Hou Xiao.Recent status of research on condensed-phase particles distribution in solid rocket motor［J］.Journal of Solid Rocket Technology，2000，23(3):25－28(in Chinese)

［2］李雅娣，陳林泉，蹇澤群.粒子尺寸分級的噴管兩相流場計算［J］.固體火箭技術(shù)，2003，26(3):32－34

Li Yadi，Chen Linquan，Jian Zequn.Numerical study of twophase nozzle flow with classified particles［J］.Journal of Solid Rocket Technology，2003，26(3):32－34(in Chinese)

［3］李強，李江，劉佩進.模型發(fā)動機內(nèi)凝相顆粒碰撞的數(shù)值模擬［J］.推進技術(shù)，2008，29(1):18－21

Li Qiang，Li Jiang，Liu Peijin.Numerical simulation of particle collision and impingement for a model SRM ［J］.Journal of Propulsion Technology，2008，29(1):18－21(in Chinese)

［4］Sabnis J S，Jong F J，Gibeling H J.Calculation of particle trajectories in solid-rocket motors with arbitrary acceleration［J］.Journal of Propulsion and Power，1992，8(5):961 －967

［5］姜貴慶，王淑華.噴管內(nèi)粒子熱增量機理研究［J］.固體火箭技術(shù)，1990，13(2):18－22

Jiang Guiqing，Wang Shuhua.The investigation of particle impingement heating mechanisms of SRM ［J］.Journal of Solid Rocket Technology，1990，13(2):18－22(in Chinese)

［6］Dunbar L E，Courtney J F，Mcmillen L D.Heating augmentation in particle erosion environments［R］.AIAA-74-607，1974

［7］Patnaik P.Handbook of inorganic chemicals［M］.New York:McGraw-Hill，2002:11 －12

［8］傅德彬，姜毅，袁增鳳.燃氣射流氣固兩相流場研究［J］.北京理工大學學報，2004，24(5):387－391

Fu Debin，Jiang Yi，Yuan Zengfeng.Numerical study on gassolid two phase jet flow ［J］.Journal of Beijing Institute of Technology，2004，24(5):387－391(in Chinese)

［9］劉君，黃琳，夏智勛.潛入式擺動噴管兩相內(nèi)流場數(shù)值模擬 ［J］.固體火箭技術(shù)，2003，26(4):18－21

Liu Jun，Huang Lin，Xia Zhixun.Numerical simulation on two-phase inner flow field of submerged vectored nozzle ［J］.Journal of Solid Rocket Technology，2003，26(4):18 －21(in Chinese)

［10］方丁酉.兩相流動力學［M］.長沙:國防科技大學出版社，1988:46－47

Fang Dingyou.Two-phase flow dynamics［M］.Changsha:National University of Defense Technology Press，1988:46 －47(in Chinese)

［11］尚慶.湍流-顆粒反應相互作用的理論分析 ［D］.北京:清華大學航天航空學院，2005

Shang Qing.Studies and application of theoretical model for turbulence-particle reaction interactions［D］. Beijing:School of Aerospace，Tsinghua University，2005(in Chinese)

［12］傅德彬，姜毅.燃氣射流氣固兩相數(shù)值模擬與顆粒沖刷分析 ［J］.彈箭與制導學報，2004，24(1):63－66

Fu Debin，Jiang Yi.Numerical study on gas-solid two-phase jet flow and solid particle erosion［J］.Journal of Projectiles Rockets Missiles and Guidance，2004，24(1):63 － 66(in Chinese)

［13］張旭.固體火箭發(fā)動機兩相流特性研究 ［D］.哈爾濱:哈爾濱工程大學航天與建筑工程學院，2006

Zhang Xu.The investigation of characteristics of two-phase flow in solid rocket motor［D］.Harbin:College of Aerospace and Civil Engineering，Harbin Engineering University，2006(in Chinese)