基于DSMC方法的高空固發羽流模擬

丁逸夫，趙 瑜，王平陽

(1.上海交通大學 機械與動力工程學院,上海 200240； 2.上海航天化工應用研究所,上海 201109)

基于DSMC方法的高空固發羽流模擬

丁逸夫1，趙 瑜2，王平陽1

(1.上海交通大學 機械與動力工程學院,上海 200240； 2.上海航天化工應用研究所,上海 201109)

為研究稀薄狀態下固體火箭發動機羽流中顆粒的行為，在稀薄氣相流場DSMC程序的基礎上，添加了氣固兩相相互作用模型和顆粒相變模型。氣固兩相相互作用模型中，顆粒所受的作用由單個顆粒受到周圍分子力和熱作用公式求出，分子所受的作用由耦合求解氣固相互作用方法求出，顆粒相變考慮固化、熔解，以及不發生相變的加熱與冷卻。用算例研究了流場的特性參數及不同粒徑顆粒的溫度分布，與文獻數據符合良好，準確度在國內同類研究中具有一定的優勢。在此基礎上，針對實際固發羽流算例進行計算，分析了不同顆粒直徑對流場的影響。結果表明：粒徑越小顆粒擴散越開，X=0.4 m處粒徑1 μm顆粒的擴散較粒徑100 μm從0.03 m增大到0.2 m；粒徑越小顆粒溫度降低越多，X=0.4 m處近軸線位置粒徑1 μm顆粒的溫度較粒徑100 μm降低了44.4%；粒徑越大，對氣相流場阻礙作用越明顯，X=0.02 m處近軸線位置粒徑100 μm顆粒的速度較純氣相流場降低了54.5%。

固體發動機； 高空羽流； 氣固兩相流； 相互作用模型； 顆粒相變模型； 流場特性參數； 溫度分布； 顆粒直徑； DSMC

0 引言

固體火箭發動機有悠久的歷史和廣泛的應用，在運載火箭的助推器和上面級、彈道導彈上面級、彈道導彈防御系統、飛行員緊急逃生裝置中都有使用。高空環境固體火箭發動機羽流中，氣體膨脹導致溫度降低，顆粒溫度明顯高于周圍氣體，是羽流輻射的主要來源，一方面對航天器羽流可見區域產生一定的影響，另一方面羽流紅外可見特性是導彈隱身和目標識別所關注的。同時，顆粒在衛星姿軌控、二級火箭分離時可能與下游航天器壁面發生碰撞，對其正常工作造成一定的危害。因此，關于固體火箭發動機的稀薄氣固兩相羽流場研究越來越受重視，對顆粒行為的模擬成為研究者關心的問題。

基于連續流假設的CFD方法無法準確模擬稀薄流場[1]。由BIRD提出的DSMC方法，經過多年的發展已成為求解稀薄流的重要方法。將DSMC方法直接用于高空固發羽流的模擬有一定的難度，主要表現在以下：顆粒相對分子幾何尺寸很大，兩者間發生的是多體碰撞；顆粒出口常呈熔融態，隨流動發生復雜的相變過程；顆粒與分子間的碰撞傳能計算；顆粒與周圍環境的輻射換熱計算等。國外對高空固發羽流的研究始于2000年，GALLIS推導出單個固體顆粒受周圍氣體分子作用力和熱的模型[2]。在此基礎上密歇根大學的BOYD和BURT等提出了耦合求解氣固兩相相互作用的技術路線，豐富完善了該模型[3-4]。國內的研究起步較晚，最早的相關報道是文獻[5]采用GALLIS提出的模型求解氣固兩相力與熱的相互作用，并初步論述了過程中動量和能量的傳輸機制。文獻[6]用DSMC/EPSM混合算法結合氣粒耦合作用模型，對羽流參數進行了模擬，但未考慮顆粒相變的影響。文獻[7]通過坐標變換計算了含液相顆粒的羽流場，但同樣未考慮相變。文獻[8]計算了高空固發羽流，同樣未考慮顆粒相變的影響。文獻[9]在文獻[5]的基礎上添加了相變模型，模擬結果與BURT結果的趨勢基本一致，但仍存在較明顯的誤差。相變會導致固體顆粒溫度升高，并導致不考慮相變的羽流熱效應偏低。本文在DSMC框架下耦合氣固兩相相互作用模型和顆粒相變模型，對高空固發羽流模擬進行了研究，對氣體數密度和不同直徑的顆粒溫度等流場參數分別與文獻[3，9]的結果進行比較。在此基礎上，建立了一個固發模型，考察不同顆粒直徑對流場參數的影響規律，以及顆粒的溫度與密度分布，為后續耦合輻射模型研究固發羽流對敏感表面的影響提供參考。

1 數學模型

對液發羽流純氣相流場的研究中使用DSMC方法求解已積累了大量經驗[10]。本文重點是與固發羽流模擬特有的顆粒相關的原理和模型。

1.1氣相對固相的作用

根據GALLIS提出的單個固體顆粒受到周圍氣體分子作用力和熱的模型，設同一個網格中，每個氣體仿真分子作用到某個固體顆粒上的力、力矩、熱量分別為FP，MP，QP，有

FP=FPnr+FPot

(1)

(2)

(3)

(4)

QP=π(RP)2×

(5)

式中：Ur為分子對顆粒的相對速度；Rp為等效顆粒半徑；Ng為仿真分子權因子；Vc為網格體積；m為單個真實氣體分子的質量；Cr為氣體分子相對固體顆粒的速度的模；τ為固體顆粒的表面熱適應系數；kB為波爾茲曼常數；Tp為固體顆粒的溫度；ωp為固體顆粒的轉動角速度；θ為氣體分子轉動自由度；er為單個真實氣體分子的平均轉動能[2]。對每個顆粒，令其與對應網格的分子發生碰撞，效果累加且作用于顆粒，得下一時刻固體顆粒的速度、角速度和溫度。

1.2固相對氣相力和熱的作用

考慮固相對氣相的影響，需選取與顆粒發生碰撞的分子。根據BIRD提出的非時間計數器方法(NTC)，先確定所有可能發生碰撞的網格中仿真分子數N，再計算碰撞概率，由顆粒以Cr以對應時間步長掃過的體積與網格體積之比決定。兩者相乘得到與顆粒發生碰撞的分子數Ns，有

Ns=floor(WpNπ(Rp)2(Cr)maxΔt/Vc+R)

(6)

式中：Wp為顆粒權因子；R為0～1隨機數；floor表示向下取整；Δt為時間步長。選取時根據接受拒絕法，生成0～1的隨機數，與每個分子Cr/(Cr)max比較，Cr/(Cr)max較大時認為該分子與顆粒發生碰撞。

與該顆粒發生碰撞的仿真分子，一部分以概率τ發生漫反射，另一部分以概率1-τ發生鏡面反射。鏡面反射不改變分子的溫度，相對速度大小在碰撞中亦不變。漫反射分子與顆粒充分交換熱量，分子以顆粒溫度為熱速度反射出去。漫反射分子反射后相對速度的概率分布可表示為

(7)

(8)

1.3固體顆粒相變模型

顆粒的相變在固發羽流中十分常見，不僅影響自身的溫度，而且對周圍氣體分子狀態產生影響。熔融態Al2O3顆粒在流場中與周圍分子相互作用溫度降低，到達向α相Al2O3轉變的溫度(1 950 K)時發生相變，結晶由表面不斷向內核發展，同時伴隨有強烈的放熱過程，溫度升高至α相Al2O3向熔融態Al2O3轉變的溫度(2 330 K)時達到平衡。此時溫度不變，顆粒對外散熱與相變產生的熱量平衡，結晶鋒面不斷向核心發展。

假設當顆粒溫度顯著低于熔點溫度Tm達到結晶溫度Tf時，顆粒由外向內發生相變，且一旦開始結晶表面不斷向內發展，結晶鋒面位置取決于

(9)

式中：r1為結晶鋒面半徑相對顆粒半徑的值；t為時間；AC=2.7×10-6ms-1·K-1.8；nc=1.8[11]。

r1通過能量平衡方程對顆粒溫度產生影響，可表示為

(10)

式中：CS為顆粒的比熱容，且CS=1 225 J/(kg·K)；mp為顆粒質量；hf為相變潛熱，且hf=1.07×106J/kg。

為描述顆粒相變過程，將顆粒相變分為三種獨立情況：固化、熔解、不發生相變的加熱和冷卻。以下給出三種情況的判斷依據，在此基礎上闡述三種情況下鋒面半徑與傳遞熱量的關系建立。

如Tp

(11)

式中：Δ(r1)3為(r1)3在時間Δt內的變化量。

如Tp≥Tm，r1<1，QP>0，顆粒將發生熔解。此時假設顆粒溫度保持定值，傳入顆粒的能量全部用于顆粒相變。令式(11)中顆粒溫度變化值為零，計算出r1的值。特別地，如r1的值大于1，表明該熱量全部用于相變仍有結余，此時將r1歸1，顆粒溫度變化可表示為

ΔTp=[hf/CS]((r1)3-1)

(12)

顆粒不發生相變，傳入或傳出顆粒的能量令顆粒溫度升高或降低，顆粒溫度變化可表示為

ΔTp=QPΔt/(CSmp)

(13)

2 驗證算例

國內外關于高空固發羽流的實驗很少，主要集中于對其宏觀輻射的測量，關于流場顆粒的測量未見相關報道，計算無法直接與實驗結果進行對比。為驗證模型的正確，與國外文獻[3]和相同算例國內文獻[9]的計算結果進行比較。

2.1邊界條件

邊界條件如圖1所示。驗證算例為軸對稱結構，對稱軸為X軸，對稱面為長80 m(-10 m到70 m)、寬30 m(0 m到30 m)的矩形。左邊界為環境來流面，來流組成為摩爾分數78%的N2和摩爾分數22%的O2，溫度288 K，馬赫數13.5，密度5.79×10-8kg/m3。在局部放大圖中可清晰看到噴管的結構，因原實驗中為給出噴管結構的具體參數，本文對照原文近似了一個噴管結構，長3.1 m，出口直徑0.78 m，圓柱段長1 m，喉部直徑0.2 m。噴管壁面均為表面溫度為300 K的漫反射面。

噴管出口為主要來流面，包含H2(摩爾分數0.38)、N2(摩爾分數0.31)和CO(摩爾分數0.31)三種組分的氣體，溫度1 433 K，速度3 113 m/s，密度0.011 kg/m3。除這三種氣體外，還包括固體顆粒Al2O3，其占噴管出口總質量流量的30%，按粒徑大小分為7組，直徑為0.3～6 μm，顆粒基本參數見表1。

表1 固體顆粒(Al2O3)參數

2.2結果驗證

氣相主要流場參數分布如圖2所示。由圖2可知：分子流出噴管后迅速膨脹，沿+X、+Y向數密度迅速減小，膨脹的最大角度約45°。這條線將計算區域分為兩部分：左側為環境來流區域，右側為高空固發羽流區域。在羽流區域，壓力急劇降低，由于環境分子稀薄，對羽流阻礙較小，分子的內能轉化為動能，靜溫降低，X、Y向速度均有顯著增加。環境來流區域，初始來流速度高，但數密度較低，與羽流區域相互作用，速度降低，壓強升高，在45°線上形成一個高壓區，動能轉為內能，且此處的羽流速度亦有所降低，溫度升高，在該位置出現一個高溫區域。

分子密度計算結果與文獻[3]的比較如圖3所示。其中文獻[3]給出了4個密度的等值線。由圖3可知：計算結果的流場結構與參考文獻較接近，表明本文計算結果較可信。

在計算區域中的(0 m，0.02 m)到(70 m，25 m)的線段上統計固體顆粒的溫度信息，并與文獻[3，9]比較。跟蹤線上氣相溫度變化如圖4所示。由圖4可知：在跟蹤線上氣相溫度迅速降低，然后逐漸趨于平緩。本文結果與文獻[3]的結果較一致，最大誤差小于60 K，此時相對誤差小于10%。跟蹤線上粒徑0.4 μm顆粒溫度變化如圖5所示。粒徑0.4 μm顆粒在噴管出口的溫度1 634 K，隨著流動與分子碰撞傳能溫度降低。由圖5可知：粒徑0.4 μm的顆粒溫度與文獻結果較一致，先迅速下降然后趨于平緩，最大誤差小于100 K，此時相對誤差小于15%。跟蹤線上粒徑4 μm顆粒溫度變化如圖6所示。粒徑4 μm顆粒溫度變化較大，初始時顆粒溫度2 178 K，隨著流動溫度降低，0.4 m位置達到1 950 K，此時熔融Al2O3發生相變結晶。由于相變向外放熱，加之相變開始溫度1 950 K明顯低于熔點2 330 K，顆粒溫度迅速上升，接近2 330 K時相變產生的熱量與顆粒和分子碰撞散發的熱量平衡，顆粒溫度保持穩定。由圖6可知：粒徑4 μm的顆粒溫度與文獻結果較一致，最大溫差約100 K，此時相對誤差約5%。綜上，可認為本文的高空固發羽流計算結果較準確。

3 固發算例

對實際高空固發羽流算例進行計算，重點比較顆粒粒徑對氣相流場分布的影響，以及顆粒的空間與溫度分布。因沒有耦合輻射模型，計算時僅選取距離噴管0.5 m的區域，此時認為耦合輻射對結果影響可忽略。

3.1邊界條件

計算中發動機為錐形噴管(如圖7所示)，噴管入口半徑6 mm，喉部半徑4 mm，R1為半徑1 mm外切圓，噴管出口半徑8.4 mm，喉部至噴管入口的距離6 mm，至噴管出口的距離7.38 mm。設總溫2 803 K，總壓2.4 MPa，喉部靜壓1.341 5 MPa；混合氣體的摩爾分數分別為H2O(0.238 0)、CO(0.131 2)、CO2(0.045 76)、H2(0.292 8)、N2(0.292 1)。取顆粒粒徑分別為1，10，100 μm，質量流量均為0.0 172 kg/s。

先用CFD方法對噴管內部含顆粒的連續流流場進行計算，取出噴管出口流場參數，再用DSMC程序進行計算羽流的稀薄流流場。網格按等比級數劃分，X向100網格，末網格與首網格之比為10；Y方向60網格，末網格與首網格之比為5。

3.2結果分析

計算所得3種粒徑顆粒的數密度分布云圖如圖8所示。由圖8可知：顆粒擴散程度隨粒徑增加而逐漸降低，X=0.5 m處粒徑1，10，100 μm的顆粒分別擴散至0.2，0.1，0.03 m處。

X=0.4 m處顆粒密度分布如圖9所示。由圖9可知：當質量流量一定時，顆粒粒徑越小擴散越開，密度也越低。X=0.4 m處顆粒溫度分布云圖如圖10所示。由圖10可知：顆粒溫度隨顆粒粒徑減小而降低，大粒徑100 μm顆粒溫度降至2 700 K，同比減小3.7%(與燃燒室溫度2 803 K比較，下同)；粒徑10 μm顆粒溫度降至2 400 K，同比減小14.4%；粒徑1 μm顆粒溫度降至1 500 K，同比減小46.5%。綜合圖9、10，分布相同時，顆粒溫度越高對反流區輻射影響就越顯著，此時粒徑越大影響越明顯；溫度相同時，擴散越開，相當于輻射面積越大，對反流區輻射影響就越顯著，此時粒徑越小影響越明顯。綜上，具體對反流區輻射影響效果還需后續進一步分析才能得出。

X=0.02 m處分子速度分布云圖如圖11所示。由圖11可知：顆粒對所在區域分子的運動起阻礙作用，且粒徑越大阻礙越明顯。近軸線處，不添加顆粒時，分子速度2 750 m/s，添加粒徑100 μm顆粒后速度降至1 250 m/s，同比減小54.5%；添加粒徑10 μm顆粒后速度降至1 770 m/s，同比減小35.6%；添加粒徑1 μm顆粒后速度降至2 630 m/s，同比減小4.4%。顆粒越大出口的速度越慢，且越集中在軸線附近，對氣體運動阻礙越明顯。

4 結束語

本文基于DSMC框架下耦合氣固相互作用模型和顆粒相變模型，對文獻[3]中的算例進行了計算，并與文獻[3，9]的結果進行比較。在粒徑1，10，100 μm三種顆粒條件下，對小型錐形噴管固體火箭發動機羽流進行了計算，并分析了計算結果，得到結論如下：第一，將GALLIS單個顆粒受周圍分子作用力和熱模型用于固發羽流模擬，計算結果與文獻符合良好，氣相最大誤差小于60 K，最大誤差出現處相對誤差小于10%，固相最大誤差小于100 K，此時相對誤差均小于15%。第二，追蹤線上的小粒徑(0.4 μm)顆粒，因未發生相變，加之氣體分子碰撞傳遞熱量，溫度一直降低；大粒徑(4 μm)顆粒，發生相變后表面溫度迅速升高，因無耦合輻射，計算范圍內溫度維持在2 330 K(α相Al2O3向熔融態的轉變溫度)。第三，顆粒加入對反流區輻射影響取決于顆粒空間分布和溫度分布，粒徑越小擴散越大，X=0.4 m處粒徑1 μm的擴散較100 μm顆粒從0.03 m增大至0.2 m，對輻射有促進作用(圖9)，但同時溫度降低更顯著。X=0.4 m處近軸線位置粒徑1 μm的溫度較100 μm顆粒降低了44.4%，對輻射有抑制作用(圖10)，但具體作用效果還需后續分析。第四，顆粒粒徑越大，對流場分子運動的阻礙越明顯(圖11)，X=0.03 m處近軸線位置添加粒徑100 μm顆粒流場的分子X向速度較純氣相流場降低了54.5%。國內外很少有關于高空固發羽流的實驗，實驗數據主要為羽流宏觀輻射特性，開展高空固發羽流流場參數測量實驗，為模擬工作提供數據支撐，是后續研究方向之一。為形成對高空固發羽流中顆粒行為的完整描述，還需在本文研究的基礎上耦合輻射模型。與氣體分子相比，顆粒相關模型中存在大量的簡化假設，部分假設對結果有顯著的影響(如Al2O3相變過程晶型轉變)，后續將分析這些簡化假設對計算結果的影響，完善顆粒相關模型。

[1] BIRD G A. Molecular gas dynamics[R]. NASA STI/Recon Technical Report A, 1976.

[2] GALLIS M A, TORCZYNSKI J R, RADER D J. An approach for simulating the transport of spherical particles in a rarefied gas flow via the direct simulation Monte Carlo method[J]. Physics of Fluids, 2001, 13(13): 3482-3492.

[3] BURT J M. Monte Carlo simulation of solid rocket exhaust plumes at high altitude[D]. Michigan: The University of Michigan, 2006.

[4] BURT J M, BOYD I D. High altitude plume simulations for a solid propellant rocket[J]. AIAA Journal, 2007, 45(12): 2872-2884.

[5] 李潔, 尹樂, 顏力. 稀薄流氣粒兩相耦合作用的熱力學模型[J]. 國防科技大學學報, 2009, 31(3): 6-10.

[6] 張斌, 毛根旺, 胡松啟, 等. 固體微推力器氣粒兩相羽流場的數值模擬[J]. 固體火箭技術, 2011, 34(3): 314-318.

[7] 何小英, 賀碧蛟, 蔡國飆. 直接模擬蒙特卡羅羽流模擬的兩相作用模型[J]. 推進技術, 2011, 32(2): 214-219.

[8] 李中華, 李志輝, 彭傲平, 等. 一種稀薄兩相流動的數值模擬方法[J]. 空氣動力學學報, 2015, 33(2): 266-271.

[9] 劉英, 李潔, 胡建峰, 等. 多相稀薄流動中的顆粒相變研究[J]. 固體火箭技術, 2014, 37(3): 330-335.

[10] 劉騰, 王平陽, 包軼穎, 等. 高超聲速稀薄流球面氣動熱工程計算公式的建立[J]. 導彈與航天運載技術, 2012(5): 32-38.

[11] PLASTININ Y, ANFIMOV N, BAULA G, et al. Modeling of aluminum oxide particle radiation in a solid propellant motor exhaust[C]∥ 31st Thermophysics Conference. [S. l.]: [n. l.], 1996: 1879.

SimulationofHighAltitudeSolidRocketPlumeBasedonDSMCMethod

DING Yi-fu1， ZHAO Yu2， WANG Ping-yang1

(1. School of Mechanical and Power Engineering, Shanghai Jiao Tong University, Shanghai 200240, China; 2. Shanghai Space Propulsion Technology Research Institute, Shanghai 201109, China)

In order to study the behavior of particles in the solid plume in the rarefied state, the gas-solid two-phase interaction model and the particle phase transition model were added on the basis of DSMC program. In the gas-solid two-phase interaction model, the force and heat of the particles are obtained by using force and thermal effects of single particle subject to the surrounding molecular. The force and heat of the molecules are obtained by the method in the technical route for solving gas-solid interaction. The particle phase change is based on the model which includes solidification, liquating and heating and cooling without phase change. The characteristic parameters of the flow field and the temperature distribution of particles with different particle sizes were studied by one sample. The results were in good agreement with the literature data. The accuracy of the results in the same research in domestic has certain advantages. The calculation of a real solid-rocket plume was carried out, and the effect of different particle diameters on the flow field was analyzed. The results show that the smaller the particle size, the more the particles diffuse. The diffusion of particles with diameter 1 μm is changed from 0.03 m to 0.2 m atX=0.4 m compared with particles with diameter 100 μm. The smaller the particle size, the more the particle temperature decreased. The temperature of particles with diameter 1 μm is decreased by 44.4% atX=0.4 m near-axis position compared with particles with diameter 100 μm. The larger the particle size, the more obvious obstruction of molecular motion in the flow field. The molecularXvelocity is decreased by 54.5% atX=0.02 m near-axis position compared with the pure gas flow field.

solid rocket engine; plume; gas-solid two-phase flow; interaction model; particle phase transition model; characteristic parameter of flow field; temperature distruction; particle diameter; DSMC

1006-1630(2017)05-0110-07

2017-02-13；

2017-03-31

國家自然科學基金資助(50306013)；航天先進技術聯合研究中心基金資助(USCAST2013-31)

丁逸夫(1992—)，男，碩士生，主要研究方向為固體火箭發動機羽流模擬。

V43

A

10.19328/j.cnki.1006-1630.2017.05.018