基于稀薄氣-粒兩相流的蒙特卡洛顆粒輻射模型研究①

徐振富，李 潔，石于中，劉 英，胡建峰

(國防科技大學(xué)航天與材料工程學(xué)院，長沙 410073)

基于稀薄氣-粒兩相流的蒙特卡洛顆粒輻射模型研究①

徐振富，李 潔，石于中，劉 英，胡建峰

(國防科技大學(xué)航天與材料工程學(xué)院，長沙 410073)

在氣-粒兩相相變模型及液態(tài)和固態(tài)顆粒碰撞、聚合和分離模型的基礎(chǔ)上，發(fā)展稀薄條件下考慮顆粒輻射的蒙特卡洛顆粒輻射模型。通過對高超聲速稀薄環(huán)境中的氣-粒兩相噴流流場的數(shù)值模擬，得到氣-粒兩相流的流場參數(shù)，利用所得流場參數(shù)作為顆粒輻射模型的初始參數(shù)進行顆粒輻射計算，同時考慮了有無探照發(fā)射時的光譜輻射強度。結(jié)果表明，在顆粒濃度較大時計算兩相稀薄流的流場參數(shù)，考慮顆粒輻射是必要的，并且考慮有無探照發(fā)射對光譜輻射強度數(shù)值的影響。

稀薄氣體動力學(xué);顆粒輻射;氣-粒兩相流;DSMC方法

0 引言

固體火箭發(fā)動機羽流中占噴流質(zhì)量流量30%的微尺寸顆粒Al2O3對紅外輻射特性的影響極大，粒子輻射的散射過程中，包含衍射、折射、透射、反射和吸收等現(xiàn)象。粒子的這種光學(xué)散射特性還將對飛行器光學(xué)元件的正常工作產(chǎn)生影響，亟需正確預(yù)測高超聲速飛行器羽流的兩相流效應(yīng)和羽流的輻射特性，以滿足空間目標識別和突防需求。隨著航天飛行器在高空的機動飛行和突防的需求發(fā)展，基于高空氣-粒兩相稀薄羽流的輻射特性研究日益受到重視。目前，Gimelshein［1］課題組在美國空軍經(jīng)費支持下，采用基于DSMC方法的SMILE軟件，研究了120 km高空的固體姿控發(fā)動機的橫向羽流干擾流場，分析了兩相羽流流場的輻射特性。密歇根大學(xué)的Burt等［2］利用Gallis推導(dǎo)的單個球形顆粒熱力學(xué)模型，基于分子動力學(xué)理論，首次實現(xiàn)了顆粒相對氣相作用的DSMC模擬，建立了雙向耦合的氣-粒兩相相互作用模型，同時探討了顆粒的非球形效應(yīng)、顆粒旋轉(zhuǎn)和顆粒相變，數(shù)值模擬了兩相羽流場的輻射特性。上述工作雖然考慮了氣-粒兩相的雙向耦合作用，為兩相羽流的輻射特性研究建立了較好的數(shù)學(xué)描述方法，但這些方法存在某些缺陷，主要是由于氣體分子和顆粒的相間作用處理模式各不相同，導(dǎo)致動量和能量的守恒僅是從時間平均層面上實現(xiàn)，不能保證單個時間步長內(nèi)的氣-粒相間碰撞的動量守恒和能量守恒，誤差將會隨著氣體分子和顆粒的數(shù)密度差距增大而增大［2］。

國內(nèi)對氣-粒兩相稀薄流的輻射特性研究工作還未見公開報道，目前絕大部分工作集中于低空連續(xù)介質(zhì)范疇的噴流輻射特性研究，而對高空稀薄條件下的氣-粒兩相羽流的輻射特性研究還停留在流場建模和氣動特性分析的工作上，如樊菁等［3］圍繞空間氣液、氣-粒多相羽流場展開研究，由于局限于空間真空環(huán)境下，未考慮氣-粒兩相的雙向耦合作用。李潔［4］基于粒子微觀行為的熱力學(xué)模型，建立了氣-粒兩相雙向耦合的動量和能量傳輸機制，構(gòu)造了氣-粒相互作用的熱力學(xué)模型。國內(nèi)的研究工作都未考慮稀薄條件下顆粒的輻射，因而忽略了顆粒輻射對流場的流動特性和輻射特性產(chǎn)生的影響。

本文在已有的工作基礎(chǔ)上，對稀薄環(huán)境下的氣-粒兩相羽流場的輻射特性的建模和數(shù)值模擬進行了完善和深入研究分析。

1 蒙特卡洛顆粒輻射模型［2］

該顆粒輻射模型應(yīng)用蒙特卡洛射線追蹤方法，用拉格朗日法追蹤通過計算網(wǎng)格的大量光子群。模型中所關(guān)心的輻射傳熱光譜部分(波長約在0.5～5 μm)被分成了一系列波段。假設(shè)Nη不同的波段，每一個寬度為Δηi并以波數(shù)ηi為中心，從穿過網(wǎng)格的來源粒子中每隔幾個時間步隨機進行挑選，從而產(chǎn)生大量具有代表性的能量束Nb。每一束都代表一定數(shù)量的能量，并且相等數(shù)量的能量束分配給每一個Nη。新產(chǎn)生的能量束根據(jù)隨機產(chǎn)生的單位矢量u來確定傳播方向。根據(jù)分配的波數(shù)段和來源粒子的性質(zhì)來賦予能量束初始能量為Pb。

以下是Plass［5］基于Mie理論計算得到的關(guān)系，來源粒子在i處的平均光譜發(fā)射率可近似表述為

式中Tp和Rp分別代表顆粒的溫度和顆粒的半徑;k是Al2O3在溫度為Tp、波數(shù)為ηi時吸收系數(shù)的值。

由于實驗數(shù)據(jù)的缺乏和k對固體顆粒不純的極度敏感，本文忽略了k對固體顆粒組成的依賴。將式(1)運用于普朗克的黑體函數(shù)，能計算出被賦值的波數(shù)范圍內(nèi)來源粒子發(fā)射能量Pp，i:

式中c0表示光在真空中的速度;h是普朗克常數(shù);kB是波爾茲曼常數(shù)。

能量束的初始能量被定義為

式中Np是網(wǎng)格中典型顆粒的總數(shù);Wp是權(quán)重因子。

Ωi(Tp)在每一個波數(shù)段模擬啟動時進行計算，在溫度為T，用實驗方法測定得到k(T，ηi)。式(2)中Ωi(Tp)的值通過對顆粒溫度的線性插值求得。

能量束選定后，每個能量束在當(dāng)前的時間步在網(wǎng)格中穿行，直到其通過流入、流出或者吸收的邊界而離開。當(dāng)能量束經(jīng)過了有顆粒布置(或者在前一時間步布置)的網(wǎng)格時，一部分分配的能量將被吸收，能量束也有可能被散射。網(wǎng)格中給定波數(shù)段的吸收和散射性質(zhì)由平均的光譜吸收比αi和散射系數(shù)σi決定。考慮到包括Nspec種不同顆粒組分j的仿真，模型中各個組分按照顆粒半徑Rj來分類。在網(wǎng)格中感興趣的是各組分顆粒都有一個平均的溫度Tj和數(shù)密度nj，此處Tj和nj的平均值通過在時間步長內(nèi)平均而得到。將式(1)和切爾科夫定律用于定義光譜吸收比，可用以下的公式計算αi:

各組分的k(Tj，ηi)值通過查表得到。與此相應(yīng)的散射系數(shù)σi以所有顆粒組分總和的形式給出:

式中Θi，j代表了顆粒組分j在波數(shù)為i時的散射效率因子。

Θi，j的值在模擬啟動時進行計算，利用了接近于Siege 和 Howell［6］的 Mie 第一原理。假設(shè)ni?k，此處ni是顆粒在波數(shù)為ηi時折射率的實部。Θi，j可由無量綱的參數(shù)xi，j=2πηiRj的一個函數(shù)給定:

當(dāng)xi，j?1 時，可假定式(7)對 Al2O3是準確的，而對于更大的xi，j值大大地超過了對于Θi，j的估算。為了使xi，j能在更大范圍內(nèi)運用，加入了限制條件Θi，j≤2。這和Plass應(yīng)用Mie理論的計算結(jié)果能較好地符合，并且避免了為尋找確切的Mie理論方法而做的詳細計算。值得注意的是，本文忽略了ni對顆粒溫度的依賴關(guān)系，資料表明，Al2O3的ni值在較大溫度范圍內(nèi)幾乎保持恒定。實驗表明ni值隨著顆粒尺寸的增加而升高，但由于缺乏有用的實驗數(shù)據(jù)，忽略了這種依賴關(guān)系。

當(dāng)能量束進入了一個給定的網(wǎng)格，沿著初始軌道離開此網(wǎng)格的總距離De由與指定的波數(shù)相應(yīng)的αi和σi值決定。散射的距離Ds由能量束傳播了距離Ds后還未發(fā)生散射的概率Pns決定，即

求解Ds，并設(shè)置Pns屬于任意數(shù) R∈［0，1］，可得到

如果Ds＞De，顆粒會沿著初始軌道離開網(wǎng)格。否則，能量束將會被散射。如果被散射，能量束將會沿著軌道移動Ds的距離。此后，其方向重新進行指定。各向異性的散射過程可用Henyey Greenstein散射相函數(shù)近似描述:

其中，自由參數(shù)g是散射角θ的平均余弦。可重設(shè)相應(yīng)的分配函數(shù)f(θ)=2πφ(θ)sinθ，如果設(shè)定:

其中，任意數(shù)R在0～1之間，在式(10)和式(11)之后，通過以下公式定義θ:

最終傳播方向的單位矢量u*按以下計算:

其中，方位角φ是［0，2π］之間的任意值;u是初始方向;單位矢量t1和t2如下:

為方便起見 被定義為沿x軸方向的單位矢量。

程序?qū)κ?9)、式(12)、式(13)進行賦值循環(huán)，直到能量束離開網(wǎng)格為止。假設(shè)能量束通過此網(wǎng)格傳播的總距離為Dt，分配的能量Pb被減少的部分占顆粒相吸收對穿透輻射強度影響的1-exp(-αiDt)。

隨著能量束通過網(wǎng)格，被半徑為Rp的單個粒子吸收的能量ΔQabs可表示如下:

式中Vcell是網(wǎng)格的體積。

能量束對于相應(yīng)波數(shù)的平均輻射能量流Δqi的貢獻如下:

在視覺上很小的網(wǎng)格內(nèi)，當(dāng)αiDt?1時，式(15)所得的結(jié)果會導(dǎo)致較大的誤差。為了糾正此誤差，當(dāng)αiDt＜10-5時，利用式(15)的線性形式計算Δqi:

值得注意的是，式(16)給出了αi=0的顆粒區(qū)域以外的Δqi的一個精確解。在能量束穿過網(wǎng)格的每一時間步，每一網(wǎng)格的能量流由通過網(wǎng)格的所有能量束的貢獻Δqi之和決定。所得結(jié)果通過在許多時間步長內(nèi)平均來減少統(tǒng)計發(fā)散。一旦流場達到了穩(wěn)定狀態(tài)的條件，將在進行輻射計算的所有時間步長內(nèi)進行平均。

如上所討論的，強烈的雙向耦合存在于流場特征和羽流輻射之間。輻射傳熱對顆粒溫度的影響很大，并會間接地影響其他性質(zhì)，如顆粒相的組成、物質(zhì)的密度以及顆粒和氣體之間的動量和能量轉(zhuǎn)移的比率。為說明輻射發(fā)射和吸收對顆粒溫度的影響，每一個典型顆粒的溫度在每個時間步長內(nèi)用以下公式修正:

式中 Δt為時間間隔;mp為顆粒的質(zhì)量;cp是顆粒的比定壓熱容;Qrad為顆粒的凈輻射傳熱率。

這里假設(shè)顆粒溫度在空間是統(tǒng)一的，基于Al2O3的小尺寸顆粒和較高的熱傳導(dǎo)率而得出的低的比奧(電流單位=10 A)數(shù)。輻射傳熱率可由式(18)計算:

在基爾霍夫定律和式(1)、式(2)之后。上述符號qi為顆粒所在網(wǎng)格中波數(shù)為i時平均時間和平均方向的能量流。

噴口內(nèi)的輻射發(fā)射顯著增加了出口平面附近的輻射率。這種輻射最初由噴管內(nèi)壁產(chǎn)生，在過去的假設(shè)下一般稱為“探照燈發(fā)射”，其主要來源于喉部上游。按照廢氣流的光學(xué)厚度，噴管內(nèi)顆粒的發(fā)射也會有很大貢獻。由于探照燈發(fā)射估計會影響羽流中的溫度和顆粒組成，應(yīng)考慮一種對輻射和流場的耦合方法進行模擬。通過在噴管出口的流入邊界處產(chǎn)生附加的能量束來說明探照燈發(fā)射。

在出口平面的每一流入邊界處描繪成為具有典型特征溫度Tw的黑體墻。沿著位于流入邊界的每個網(wǎng)格面，在進行輻射計算的每個時間步長產(chǎn)生新的能量束Nf。每一能量束都隨機的分配了波數(shù)i，并被賦予了初始的方向u為

式中nf為網(wǎng)格表面的單位內(nèi)法線矢量;θ是u相對于表面的天頂角;R是［0，1］之間的任意數(shù)。

初始能量Pb從式(19)中定義，Af為表面積，Nη/Nf為權(quán)重因子，并對相應(yīng)波數(shù)段的普朗克函數(shù)積分:

為了計算羽流輻射率，一個或更多的模擬輻射計放在了網(wǎng)格區(qū)域以外的某個地方。考慮一個表面積為As，單位外法向為ns，角分辨率由天頂角ω定義的傳感器。當(dāng)每一個能量束離開網(wǎng)格時，可判斷沿著單位矢量為u的能量束的軌跡是否會與傳感器表面相交。沿著由光束方向的單位矢量確定的軌道橫穿傳感器表面。如果相交發(fā)生，并滿足條件-u·ns≥cosω，在當(dāng)前的時間步內(nèi)，分配給光束的能量Pb加入到了相應(yīng)波數(shù)的總的吸收能量ΣPi中。平均波段光譜輻射強度Ii的瞬時值可看作ΣPi與傳感面積、吸收輻射的立體角、波束波長的比。可得到

其中:

每一束的Ii值在較多時間步長內(nèi)進行平均，以減少發(fā)散，取樣僅在流場達到定常態(tài)進行。

2 算例描述

本文采用了從噴管出口平面向下延伸100 m，向外擴展40 m的矩形區(qū)域。噴管出口處的直徑為7.85 cm，顆粒和氣體在出口平面的數(shù)據(jù)都從Anfimov［7］中得到，數(shù)據(jù)是基于Star-27發(fā)動機中顆粒質(zhì)量分數(shù)為30%的管口流動特征的模擬。數(shù)據(jù)中的氣體為N2、H2和CO的混合物，氣相之間的碰撞用可變硬球碰撞模型。在出口平面，氣體的速度達到了3 113 m/s，溫度為1 433 K，密度為0.011 kg/m3，摩爾分數(shù)H2為0.38，N2和CO各占0.31。顆粒相的尺寸分布是離散的，按粒徑分為7個不同的種類，直徑范圍在0.3～6 μm之間。由于噴管出口流場信息的缺乏，由Anfimov［7］給出的顆粒性質(zhì)在整個出口平面內(nèi)認為是相同的。噴管出口處顆粒的性質(zhì)由表1給出。

設(shè)顆粒表面的導(dǎo)熱調(diào)節(jié)系數(shù)為0.9，這樣90%涉及到漫反射的相間碰撞都完全適應(yīng)了顆粒溫度，其余的10%涉及到了鏡面反射。相間的動量和能量轉(zhuǎn)換用雙向耦合法進行計算，液態(tài)Al2O3小滴的結(jié)晶化作用用一個非平衡的相變模型描述。相變模型解決了溫度對結(jié)晶率和相關(guān)放熱的依賴，并忽略了固態(tài)Al2O3的γ-α的轉(zhuǎn)變及不同相密度的差異。

表1 噴流出口顆粒相參數(shù)分布Table 1 Particle properties at the nozzle exit

顆粒吸收系數(shù)k的值從 Konopka、Reed、Calia［8］的實驗數(shù)據(jù)中得到。基于此數(shù)據(jù)，本文使用了10個波數(shù)，相當(dāng)于中紅外的范圍為1.3～4.5 μm 的波長。通過對文獻［9］中2種固體火箭發(fā)動機的廢氣流中的顆粒數(shù)據(jù)收集及研究，發(fā)現(xiàn)第2種(火箭2)所給出的k值更多地與其他實驗數(shù)據(jù)和相關(guān)文獻相一致。這里應(yīng)用由SEM測量法得到的第二種流動計算值。

根據(jù)文獻［10］，本文設(shè)散射角的平均余弦g=0.5，并對噴管出口的探照發(fā)射運用了一個有效溫度Tw=1 300 K。

3 結(jié)果分析

圖1為本文所得及文獻［2］中的顆粒質(zhì)量密度等值線圖的對比。

圖1 顆粒質(zhì)量密度等值線圖Fig.1 Contours of particles＇mass density

如圖1(a)所示，首先注意到顆粒僅在大體上為一半的仿真區(qū)域中出現(xiàn)，在噴管出口鄰近區(qū)域的羽流中，顆粒被擴張的氣體從中心線往外推移，由于最大顆粒相的發(fā)散角被顆粒質(zhì)量所限制及顆粒尺寸分布范圍的存在，所以大顆粒基本上聚集在中心線附近，顆粒質(zhì)量密度隨著與中心線的距離的增大而逐漸減小，較小的顆粒將會在離中心線較遠的區(qū)域找到。由于模擬區(qū)域尺寸規(guī)模比噴管出口半徑大幾個數(shù)量級，在圖1(a)中所列出的顆粒特性僅反映了在較遠區(qū)域的一種趨勢，在較遠區(qū)域的顆粒和氣體之間動量和能量的轉(zhuǎn)換假定為可忽略，顆粒在遠離噴口處幾乎沿著直線前進。

由圖1(b)可見，文獻中的顆粒質(zhì)量密度等值線明顯彎向中心線，而本文的計算結(jié)果較平滑。這是由于本文采用的兩相流模型與文獻采用的模型有所不同，從而使得顆粒質(zhì)量密度分布狀態(tài)有所差異。

圖2為有、無探照發(fā)射的光譜輻射強度等值線圖。由圖2(a)可知，在軸向隨著與噴口距離的增加，有探照發(fā)射的光譜輻射強度逐漸減小。這主要是由于光譜輻射強度與顆粒的質(zhì)量密度及顆粒的尺寸成比例，在遠離噴口處，顆粒的質(zhì)量密度減小，大尺寸顆粒的數(shù)量也減小，故而輻射強度減小。同理，在徑向顆粒的質(zhì)量密度和大尺寸顆粒的數(shù)量也隨著距中心線距離的增加而減小，所以徑向光譜輻射強度也隨著距中心線距離的增加而減小。

圖2 光譜輻射強度等值線圖Fig.2 Contours of spectral radiance

對比有無探照發(fā)射的光譜輻射強度等值線圖可知，在有無探照發(fā)射的情況下，輻射強度等值線的趨勢基本相同，不同的是在有探照發(fā)射的情況下，由于探照發(fā)射這種點聲源特性的作用，在噴口附近徑向的輻射強度增加。因此，有探照發(fā)射的光譜輻射強度等值線圖較平緩些。

圖3為X=0出口截面輻射強度沿Y向變化的對比曲線圖。在中心線與噴口之間，輻射強度先減小后增加。這可能是由于在中心線附近大顆粒較多，導(dǎo)致輻射強度較大，稍向外由于氣體的膨脹，帶動較小顆粒的向兩邊移動，導(dǎo)致輻射強度下降。而在噴口邊緣，較小顆粒被氣體膨脹所帶動，在噴口邊緣積聚了較多顆粒，使得輻射強度增加，而在噴口邊緣以外由于顆粒數(shù)量的急劇減小，所以輻射強度急劇下降。由圖3可知，在有無探照發(fā)射時和顆粒吸收的輻射強度的變化趨勢基本一致，造成原因相同，只是變化幅度有所差異。

圖3 出口截面輻射強度變化對比曲線Fig.3 Spectral radiance change contrast of exit section

4 結(jié)論

(1)如果忽略了各網(wǎng)格中心之間的彎曲變形以及統(tǒng)計發(fā)散的影響，從噴管出口噴出來的顆粒軌跡應(yīng)是直的，而顆粒質(zhì)量密度由于顆粒的膨脹沿下游降低。

(2)光譜輻射強度的強弱與顆粒的質(zhì)量密度成正比，說明了在兩相稀薄流中的顆粒輻射對流場影響較大。在顆粒濃度較大時，計算流場參數(shù)時，考慮顆粒輻射是必要的。

(3)在計算流場的顆粒輻射數(shù)據(jù)時，考慮有無探照發(fā)射時的光譜輻射強度是有區(qū)別的，說明在計算光譜輻射強度時，要針對具體問題具體分析;雖然得到的結(jié)果趨勢大致一致，但數(shù)值上有所差別。

［1］Sergey Gimeishei，GenIlady Marltclov Jean Muyiacrt.Numerical modeling of low thrust solid propellant nozzles at high altitudes［R］.AIAA 2006-3273.

［2］Burt J M，Boyd I D.A Monte Carlo radiation model for simulating rarefied multiphase plume flows［R］.AIAA 2005-4691.

［3］樊菁，劉宏立，蔣建政，等.火箭剩余推進劑排放過程的分析與模擬［J］.力學(xué)學(xué)報，2004，36(2):129-139.

［4］李潔，任兵，陳偉芳.稀薄流過渡區(qū)氣-粒兩相噴流的建模與數(shù)值模擬［J］.空氣動力學(xué)學(xué)報，2005，23(4):484-489.

［5］Plass G N.Temperature dependence of the Mie scattering and absorption cross sections for aluminum oxide［J］.Applied Optics，1965，4(12):1616-1619.

［6］Siegel R，Howell J R.Thermal radiation heat transfer［M］.Washington:Hemisphere Publishing，1981.

［7］Anfimov N A，Karabadjak G F，Khmelinin B A，et al.Analysis of mechanisms and nature of radiation from aluminum oxide in different phase states in solid rocket exhaust plumes［R］.AIAA 93-2818.

［8］Konopka W L，Reed R A，Calia V S.Measurements of infrared optical properties of Al2O3rocket particles［R］.AIAA 83-1568.

［9］Duval R，Soufiani A，Taine J.Coupled radiation and turbulent multiphase flow in an aluminised solid propellant rocket engine［J］.Journal of Quantitative Spectroscopy and Radiative Transfer，2004，84(4):513-526.

［10］Reed R A，Beale K S，Neese D W，et al.The effect of seachlight emission on radiation from solid rocket plumes［R］.AIAA 92-2918.

Study of Monte Carlo particles radiation model based on rarefied gas-particle two-phase plume flows

XU Zhen-fu，LI Jie，SHI Yu-zhong，LIU Ying，HU Jian-feng
(College of Aerospace and Material Engineering，National Univ.of Defense Technology，Changsha 410073，China)

Based on the phase change model of gas-particle two-phase and the collision，consolidation and separation model of liquid and solid particals，the Monte Carlo particle radiation model for rarefied flows was developed.Through simulating the gas-particle two-phase jet flow field of hypersonic rarefied flows，the parameter of two-phase jet flow field was obtained.By using the parameter as the initial coming parameter for the particles radiation model，and considering spectral radiant intensity，particle radiation was calulated.The results show that it＇s essential to consider the particles radiation for the hypersonic rarefied gas-particle two-phase flows，and the searchlight emission will impact the numerical result.

rarefied gas dynamics;particles radiation;gas-particle two-phase flow;DSMC method

V435

A

1006-2793(2012)04-0479-06

2011-11-17;

2011-12-15。

國家自然科學(xué)基金資助項目(NO.10602063)。

徐振富(1984—)，男，碩士生，主要研究計算流體力學(xué)。E-mail:zhenfu.100@163.com

(編輯:呂耀輝)