基于DSMC方法的高空固發(fā)羽流模擬

丁逸夫，趙 瑜，王平陽

(1.上海交通大學(xué) 機(jī)械與動(dòng)力工程學(xué)院,上海 200240； 2.上海航天化工應(yīng)用研究所,上海 201109)

基于DSMC方法的高空固發(fā)羽流模擬

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

(1.上海交通大學(xué) 機(jī)械與動(dòng)力工程學(xué)院,上海 200240； 2.上海航天化工應(yīng)用研究所,上海 201109)

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

固體發(fā)動(dòng)機(jī)； 高空羽流； 氣固兩相流； 相互作用模型； 顆粒相變模型； 流場特性參數(shù)； 溫度分布； 顆粒直徑； DSMC

0 引言

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

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

1 數(shù)學(xué)模型

對(duì)液發(fā)羽流純氣相流場的研究中使用DSMC方法求解已積累了大量經(jīng)驗(yàn)[10]。本文重點(diǎn)是與固發(fā)羽流模擬特有的顆粒相關(guān)的原理和模型。

1.1氣相對(duì)固相的作用

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

FP=FPnr+FPot

(1)

(2)

(3)

(4)

QP=π(RP)2×

(5)

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

1.2固相對(duì)氣相力和熱的作用

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

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

(6)

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

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

(7)

(8)

1.3固體顆粒相變模型

顆粒的相變?cè)诠贪l(fā)羽流中十分常見，不僅影響自身的溫度，而且對(duì)周圍氣體分子狀態(tài)產(chǎn)生影響。熔融態(tài)Al2O3顆粒在流場中與周圍分子相互作用溫度降低，到達(dá)向α相Al2O3轉(zhuǎn)變的溫度(1 950 K)時(shí)發(fā)生相變，結(jié)晶由表面不斷向內(nèi)核發(fā)展，同時(shí)伴隨有強(qiáng)烈的放熱過程，溫度升高至α相Al2O3向熔融態(tài)Al2O3轉(zhuǎn)變的溫度(2 330 K)時(shí)達(dá)到平衡。此時(shí)溫度不變，顆粒對(duì)外散熱與相變產(chǎn)生的熱量平衡，結(jié)晶鋒面不斷向核心發(fā)展。

假設(shè)當(dāng)顆粒溫度顯著低于熔點(diǎn)溫度Tm達(dá)到結(jié)晶溫度Tf時(shí)，顆粒由外向內(nèi)發(fā)生相變，且一旦開始結(jié)晶表面不斷向內(nèi)發(fā)展，結(jié)晶鋒面位置取決于

(9)

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

r1通過能量平衡方程對(duì)顆粒溫度產(chǎn)生影響，可表示為

(10)

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

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

如Tp

(11)

式中：Δ(r1)3為(r1)3在時(shí)間Δt內(nèi)的變化量。

如Tp≥Tm，r1<1，QP>0，顆粒將發(fā)生熔解。此時(shí)假設(shè)顆粒溫度保持定值，傳入顆粒的能量全部用于顆粒相變。令式(11)中顆粒溫度變化值為零，計(jì)算出r1的值。特別地，如r1的值大于1，表明該熱量全部用于相變?nèi)杂薪Y(jié)余，此時(shí)將r1歸1，顆粒溫度變化可表示為

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

(12)

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

ΔTp=QPΔt/(CSmp)

(13)

2 驗(yàn)證算例

國內(nèi)外關(guān)于高空固發(fā)羽流的實(shí)驗(yàn)很少，主要集中于對(duì)其宏觀輻射的測(cè)量，關(guān)于流場顆粒的測(cè)量未見相關(guān)報(bào)道，計(jì)算無法直接與實(shí)驗(yàn)結(jié)果進(jìn)行對(duì)比。為驗(yàn)證模型的正確，與國外文獻(xiàn)[3]和相同算例國內(nèi)文獻(xiàn)[9]的計(jì)算結(jié)果進(jìn)行比較。

2.1邊界條件

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

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

表1 固體顆粒(Al2O3)參數(shù)

2.2結(jié)果驗(yàn)證

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

分子密度計(jì)算結(jié)果與文獻(xiàn)[3]的比較如圖3所示。其中文獻(xiàn)[3]給出了4個(gè)密度的等值線。由圖3可知：計(jì)算結(jié)果的流場結(jié)構(gòu)與參考文獻(xiàn)較接近，表明本文計(jì)算結(jié)果較可信。

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

3 固發(fā)算例

對(duì)實(shí)際高空固發(fā)羽流算例進(jìn)行計(jì)算，重點(diǎn)比較顆粒粒徑對(duì)氣相流場分布的影響，以及顆粒的空間與溫度分布。因沒有耦合輻射模型，計(jì)算時(shí)僅選取距離噴管0.5 m的區(qū)域，此時(shí)認(rèn)為耦合輻射對(duì)結(jié)果影響可忽略。

3.1邊界條件

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

先用CFD方法對(duì)噴管內(nèi)部含顆粒的連續(xù)流流場進(jìn)行計(jì)算，取出噴管出口流場參數(shù)，再用DSMC程序進(jìn)行計(jì)算羽流的稀薄流流場。網(wǎng)格按等比級(jí)數(shù)劃分，X向100網(wǎng)格，末網(wǎng)格與首網(wǎng)格之比為10；Y方向60網(wǎng)格，末網(wǎng)格與首網(wǎng)格之比為5。

3.2結(jié)果分析

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

X=0.4 m處顆粒密度分布如圖9所示。由圖9可知：當(dāng)質(zhì)量流量一定時(shí)，顆粒粒徑越小擴(kuò)散越開，密度也越低。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，分布相同時(shí)，顆粒溫度越高對(duì)反流區(qū)輻射影響就越顯著，此時(shí)粒徑越大影響越明顯；溫度相同時(shí)，擴(kuò)散越開，相當(dāng)于輻射面積越大，對(duì)反流區(qū)輻射影響就越顯著，此時(shí)粒徑越小影響越明顯。綜上，具體對(duì)反流區(qū)輻射影響效果還需后續(xù)進(jìn)一步分析才能得出。

X=0.02 m處分子速度分布云圖如圖11所示。由圖11可知：顆粒對(duì)所在區(qū)域分子的運(yùn)動(dòng)起阻礙作用，且粒徑越大阻礙越明顯。近軸線處，不添加顆粒時(shí)，分子速度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%。顆粒越大出口的速度越慢，且越集中在軸線附近，對(duì)氣體運(yùn)動(dòng)阻礙越明顯。

4 結(jié)束語

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

[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] 李潔, 尹樂, 顏力. 稀薄流氣粒兩相耦合作用的熱力學(xué)模型[J]. 國防科技大學(xué)學(xué)報(bào), 2009, 31(3): 6-10.

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

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

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

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

[10] 劉騰, 王平陽, 包軼穎, 等. 高超聲速稀薄流球面氣動(dòng)熱工程計(jì)算公式的建立[J]. 導(dǎo)彈與航天運(yùn)載技術(shù), 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

國家自然科學(xué)基金資助(50306013)；航天先進(jìn)技術(shù)聯(lián)合研究中心基金資助(USCAST2013-31)

丁逸夫(1992—)，男，碩士生，主要研究方向?yàn)楣腆w火箭發(fā)動(dòng)機(jī)羽流模擬。

V43

A

10.19328/j.cnki.1006-1630.2017.05.018