固體火箭發動機噴管內氣粒兩相流動的CFD-DSMC模擬①

吳限德，張 斌，陳衛東，夏廣慶，陳茂林

(1.哈爾濱工程大學航天與建筑工程學院，哈爾濱 150080;2.西北工業大學燃燒、流動和熱結構國家級重點實驗室，西安 710072;3.大連理工大學工業裝備結構分析國家重點實驗室航空航天學院，大連 116024)

0 引言

從20世紀后期起，美軍就開始了空天項目的研究與開發，其中包括作戰保障衛星系統、反導武器系統、反衛星武器、空天飛機、軌道轟炸機等。在眾多的空天項目中，各類固體火箭發動機作為一種重要動力裝置，扮演著重要角色。如超高聲速攔截彈中采用的側噴直接力控制技術，在美國的PAC-3、俄羅斯的S-400和歐共體的ASTER中已獲得成功運用［1－3］。在許多情況下，由于從固體火箭發動機噴管噴出的高溫高速顆粒的沖刷效應會給空天平臺的熱防護系統造成嚴重破壞，此時氣粒兩相稀薄羽流對空天平臺的影響不容忽視。因此，急需開展高空環境中固體火箭發動機氣粒兩相流動的研究工作。

為了能更好地模擬高空環境中噴管內外的氣粒兩相流動，需將二者結合進行一體化計算。實現一體化計算的難點并不在于CFD方法和DSMC方法在模擬氣相流動時界面處的處理，而是在于界面處顆粒相的處理。在連續流動中，顆粒相的描述通常采用“顆粒軌道模型”或“擬流體模型”［4］。

以此兩種方法模擬得到界面處顆粒相的參數，很難在稀薄氣粒兩相流動的模擬中再現顆粒運動的隨機性和真實性，同時不良的界面處理方法也會引入較大的人為誤差。采用DSMC方法描述連續流中顆粒相的運動，則在界面處無需過多處理(僅需更換相互作用模型)，即可實現對其在稀薄流中的描述。

目前，國內僅有國防科大的陳偉芳［5］和黃琳［6］開展了超聲速管道和氣粒兩相自由射流的CFD-DSMC研究，但未見對噴管這一典型的壓-跨-超流動進行CFDDSMC的氣粒兩相模擬。文中采用CFD-DSMC方法，模擬了JPL噴管中的連續氣粒兩相流動，研究不同顆粒含量和顆粒粒徑下氣相和顆粒相的流動規律，并與“顆粒確定性軌道模型”和“擬流體模型”得到的結果進行了對比。

1 氣相數值方法

1.1 控制方程

氣相控制方程為不計體積力和無內熱源的三維非定常N-S方程組。在笛卡爾坐標系下，其積分形式如式(1)所示。

控制方程中其余各項的物理意義和具體表達式見文獻［1］。

1.2 計算方法

采用中心有限體積法求解。對流項采用中心差分格式，并添加各向異性人工粘性求解，擴散項采用中心差分格式求解，時間項采用顯式四步Runge-Kutta法求解，同時采用隱式殘值光順法加速收斂，湍流粘性系數采用B-L模型求解。

1.3 邊界條件

氣相入口給定來流總溫、總壓、來流馬赫數和方向角;氣相出口處參數由內場按二階外插獲得;氣相壁面采用絕熱、無滑移固壁邊界條件。

模擬工況:噴管入口處總壓1.034 2 MPa，總溫555.0 K，來流馬赫數0.1，來流方向角0°。

2 顆粒相數值方法

2.1 氣相對顆粒相作用的數學模型

顆粒相的位置、速度和溫度的計算公式分別如式(2)～式(4)所示。

式中 下標k表示第k束顆粒;Dk為第k束顆粒的直徑;為第k束顆粒的滑移速度;CDk為第k束顆粒的滯止系數;Cp和C分別為顆粒和氣體的比熱容;Nuk為第k束顆粒的努謝爾特數;Pr為Prandtl數。

2.2 顆粒相對氣相作用源項的計算

顆粒相對氣相作用源項的計算公式如式(5)所示。

式中 Vi，j，k為網格單元體積;Nk為第k束顆粒中的真實顆粒數。

2.3 邊界條件

顆粒相入口給定顆粒的質量流率、溫度、速度和相關物性參數(如密度、比熱容)。在每個時刻步下，入口邊界均有仿真顆粒進入計算區域，仿真顆粒的位置按隨機方式給出。

顆粒相出口處按吸收邊界處理，同時注銷顆粒的編號和相應屬性。

顆粒相壁面邊界采用 Tabakoff［7］的經驗公式，確定顆粒碰撞后的速度方向和大小。

模擬工況:顆粒的速度和溫度與氣相相同，顆粒密度為 4 004.62 kg/m3，顆粒比熱容為 1 380 J/(kg·K)。計算中，選擇的顆粒質量分數分別為0、0.10、0.30，選擇的顆粒粒徑分別為 1、10、20 μm。

3 計算結果與分析

JPL噴管的幾何結構參見文獻［7］。計算中，氣相的計算區域為實體噴管的1/4，采用結構化網格計算;顆粒相的計算區域為氣相計算區域的一個扇面，采用軸對稱模型描述，每一時刻下對顆粒的速度和位置進行軸對稱修正，并在徑向對顆粒做復制-刪除操作;兩相耦合算法采用PSIC(Particle Source In Cell)算法，將顆粒相對氣相的作用源項按位置施加給氣相計算區域中的所有扇面。

3.1 氣相和顆粒相參數分布規律

圖1和圖2給出了在顆粒質量分數0.30、顆粒粒徑1 μm計算條件下，得到的噴管中氣相和顆粒相的主要參數分布情況。由圖1可看出，顆粒相對氣相流動的阻滯作用非常明顯，能顯著影響氣相流場中的速度和馬赫數分布規律。同時，顆粒相對氣相的加熱作用也十分顯著。

圖1 氣相參數分布Fig.1 Contours of gas-phase flow parameters

從圖2可看出，顆粒在氣相作用下加速運動，同時與氣相間的傳熱，導致自身溫度明顯降低。由顆粒相濃度云圖(c)可看出，壁面附近顆粒濃度最大。同時，在噴管喉部處，顆粒濃度梯度明顯，且在軸線附近顆粒濃度較小。在相同工況下，由顆粒確定軌道模型計算得到的無粒子區中，仍有少量顆粒存在，且其分布具有一定的非連續性。

圖2 顆粒相參數分布Fig.2 Contours of particle-phase parameters

3.2 顆粒質量分數對氣相流動參數的影響

圖3給出了不同顆粒質量分數下，氣相在軸線上馬赫數的分布。由圖3可知，隨顆粒質量分數的增加，氣相在軸線上的馬赫數逐漸降低。這是由于顆粒質量分數越高，單位體積內的顆粒數密度越大，對氣相流動的阻礙作用越大。

圖4給出了不同顆粒質量分數下，氣相在軸線上溫度的分布。由圖4可知，隨顆粒質量分數的增加，氣相在軸線上的溫度逐漸上升。這是由于顆粒質量分數越大，單位體積內顆粒對氣相的熱交換更為強烈，從而導致氣相升溫。以上結果與曾卓雄［8］采用擬流體模型獲得的結果的分布趨勢相同。

圖3 不同顆粒質量分數下氣相在軸線上的馬赫數分布Fig.3 Distributions of Ma number along the axis at different concentrations of particles

3.3 顆粒粒徑對氣相流動參數的影響

圖5和圖6分別為不同顆粒直徑下，氣相在軸線上的馬赫數和溫度分布。顆粒直徑越小，其表面積越大，與氣相動量和能量的交換越顯著。因此，隨顆粒粒徑的減小，氣相在軸線上的馬赫數逐漸降低，氣相在軸線上的溫度逐漸升高。

以上結果與于勇等［7］采用顆粒確定性軌道模型獲得的結果分布趨勢相同。

圖5 不同顆粒粒徑下氣相在軸線上的馬赫數分布Fig.5 Distributions of Ma number along the axis at different diameters of particles

圖6 不同顆粒粒徑下氣相在軸線上的溫度分布Fig.6 Distributions of temperature along the axis at different diameters of particles

4 結論

采用CFD-DSMC方法，模擬了JPL噴管中的氣粒兩相流動，得到了氣相和顆粒相的主要參數分布規律，并研究了顆粒含量和顆粒粒徑對氣相流場的影響。研究表明:

(1)顆粒含量越高，對氣相的減速作用越大，同時對氣相的加熱效應也越明顯。

(2)顆粒的粒徑越小，對氣相的減速作用越大，同時對氣相的加熱效應也越明顯。

(3)與采用顆粒確定性軌道模型計算得到的無粒子區對比發現，按照文中計算方法，該區域仍有少量且非連續的顆粒存在。

［1］徐敏，陳士櫓.側噴干擾流場計算及分析［J］.西北工業大學學報，2003，21(3).

［2］Roger R P.The aerodynamics of jet thruster control for supersonic/hypersonic endo-interceptors:lessons learned［R］.AIAA 1999-0804.

［3］Graham M J，Weinacht P.Numerical simulation of lateral control jets［R］.AIAA 1999-0510.

［4］方丁酉，張為華，楊濤.固體火箭發動機內彈道學［M］.長沙:國防科技大學出版社，1997.

［5］陳偉芳，常雨.三維管道超聲速氣固兩相流動的CFD/DSMC 仿真［J］.推進技術，2005，26(3).

［6］黃琳，劉君.氣-固兩相自由射流的粒子仿真方法［J］.國防科技大學學報，2002，24(3).

［7］于勇，劉淑艷，張世軍，等.固體火箭發動機噴管氣固兩相流動的數值研究［J］.航空動力學報，2009，24(4).

［8］曾卓雄.密相可壓氣固兩相湍流流動的理論分析及其數值模擬［D］.西安:西安交通大學，2001.