劉 仔,李艷臣

(1 哈爾濱工程大學(xué)航天與建筑工程學(xué)院,哈爾濱 150001;2 上海衛(wèi)星裝備研究所,上海 200240)

?

五種顆粒平均直徑計(jì)算方法模擬研究

劉 仔1,李艷臣2

(1 哈爾濱工程大學(xué)航天與建筑工程學(xué)院,哈爾濱 150001;2 上海衛(wèi)星裝備研究所,上海 200240)

為了確定更加合理的顆粒平均直徑計(jì)算方法,文中分別模擬研究了多尺寸分布與單尺寸分布的發(fā)動(dòng)機(jī)內(nèi)流場(chǎng)與性能參數(shù)的差異。采用高雷諾數(shù)條件下的k-ε湍流模型與顆粒軌道模型進(jìn)行數(shù)值模擬。模擬結(jié)果表明,5種顆粒平均直徑計(jì)算方法得到的單尺寸分布的計(jì)算結(jié)果與多尺寸分布的計(jì)算結(jié)果均存在偏差。采用質(zhì)量平均法與體積表面積平均法求解顆粒平均直徑較好,而線性平均法不宜采用。

顆粒軌道模型;質(zhì)量平均法;體積表面積平均法;線性平均法

0 引言

兩相流動(dòng)在固體火箭發(fā)動(dòng)機(jī)內(nèi)普遍存在,研究發(fā)動(dòng)機(jī)內(nèi)的兩相流場(chǎng)對(duì)發(fā)動(dòng)機(jī)的合理設(shè)計(jì)具有重大的意義,因此國(guó)內(nèi)外通過數(shù)值模擬對(duì)固體火箭發(fā)動(dòng)機(jī)內(nèi)的兩相流場(chǎng)與發(fā)動(dòng)機(jī)性能參數(shù)進(jìn)行了大量的研究工作[1-6]。發(fā)動(dòng)機(jī)中的顆粒存在粒徑的顆粒數(shù)分布和粒徑的顆粒質(zhì)量分布兩種多尺寸分布[7]。雖然多尺寸分布計(jì)算結(jié)果優(yōu)于單尺寸分布計(jì)算結(jié)果[8],但是當(dāng)采用多尺寸分布進(jìn)行數(shù)值計(jì)算時(shí),會(huì)導(dǎo)致計(jì)算機(jī)的計(jì)算量成倍增加,導(dǎo)致計(jì)算效率會(huì)降低,因此為了縮短計(jì)算機(jī)的計(jì)算時(shí)間,一般數(shù)值模擬中都采用單尺寸分布代替多尺寸分布。目前能夠代替多顆粒的單顆粒直徑計(jì)算方法有5種,研究5種方法中哪種方法計(jì)算得到的發(fā)動(dòng)機(jī)性能與兩相流場(chǎng)最接近多尺寸分布具有較大的工程實(shí)際意義。

1 數(shù)值計(jì)算方法

1.1 氣相控制方程與湍流模型選擇

氣相控制方程的通用形式為:

(1)

其中:φ是任意獨(dú)立變量;Γφ是輸運(yùn)系數(shù);Sφ是氣相源項(xiàng);Spφ是氣相與顆粒之間相互作用源項(xiàng)。在進(jìn)行數(shù)值模擬過程中,文中選擇RNGk-ε模型兩方程湍流模型,采用標(biāo)準(zhǔn)壁面函數(shù)法對(duì)發(fā)動(dòng)機(jī)壁面進(jìn)行處理。

1.2 顆粒相控制方程

連續(xù)方程:

(2)

其中:A是發(fā)動(dòng)機(jī)的通道面積;vkn是垂直于通道截面的顆粒分速度。

在i方向上的動(dòng)量方程:

(3)

其中:τrk是顆粒相的動(dòng)量弛豫時(shí)間;gi是重力加速度在i方向上的分量。

能量方程:

(4)

其中,在方程的右邊分子中的第一項(xiàng)是顆粒與氣相之間對(duì)流換熱,而后一項(xiàng)是顆粒與氣相之間的輻射換熱。TR是輻射溫度,εp是發(fā)射率,σ是玻爾茲曼常數(shù),ck是顆粒的比熱。

1.3 物理模型與邊界條件

計(jì)算模型如圖1所示。其中兩相流中的顆粒相為Al2O3,氣相為推進(jìn)劑燃燒產(chǎn)生的燃?xì)?。Al2O3質(zhì)量分?jǐn)?shù)為25.0%,密度為3 960.0 kg/m3,比熱為1 408.0 J/(kg·K);燃?xì)獾亩▔罕葻釣? 811.0 J/(kg·K),相對(duì)原子質(zhì)量為20.174,比熱比k為1.226。其中入口a設(shè)定為質(zhì)量流量入口,燃?xì)馀c顆粒的質(zhì)量流量為4.0 kg/s,顆粒與燃?xì)獾臏囟仍O(shè)定為3 200.0 K,壁面b設(shè)置為無滑移壁面,顆粒與壁面之間的作用設(shè)置為反彈,對(duì)稱軸c設(shè)置為對(duì)稱邊界,出口d設(shè)置為壓力出口。

圖1 物理模型

2 算例驗(yàn)證

文中對(duì)文獻(xiàn)[9]的數(shù)據(jù)進(jìn)行模擬對(duì)比。圖2給出了文中與文獻(xiàn)中顆粒的分布對(duì)比,驗(yàn)證了文中計(jì)算方法的合理性與正確性。

圖2 顆粒直徑為3.90 μm時(shí)噴管中的無顆粒區(qū)對(duì)比

3 不同顆粒平均直徑計(jì)算方法

其中對(duì)多顆粒進(jìn)行求解平均直徑計(jì)算方法有5種[7],包括線性平均法、表面積平均法、體積平均法、體積表面積平均法和質(zhì)量平均法。選取平均直徑計(jì)算方法的原則是按照顆粒平均直徑計(jì)算所得的結(jié)果應(yīng)該與按多尺寸分布計(jì)算所得的結(jié)果相近,因此需要根據(jù)具體問題合理的選擇顆粒平均直徑,下面給出5種顆粒平均半徑計(jì)算方法的計(jì)算公式。

(5)

(6)

(7)

(8)

(9)

文中采用文獻(xiàn)[7]中給出的某種多尺寸分布,其中多尺寸分布規(guī)律如表1。利用5種顆粒平均直徑計(jì)算方法對(duì)該多尺寸分布進(jìn)行平均直徑的計(jì)算,其中得到的5種顆粒平均直徑計(jì)算結(jié)果如表2。

表1 多尺寸分布

表2 不同計(jì)算方法的計(jì)算結(jié)果

4 計(jì)算結(jié)果分析

圖3與圖4給出了5種不同的單尺寸分布與多顆粒分布下發(fā)動(dòng)機(jī)軸線上的氣相馬赫數(shù)與溫度分布。由圖3可知,發(fā)動(dòng)機(jī)軸線上的氣相馬赫數(shù)在噴管收斂段之前5種方法得到的結(jié)果與多尺寸分布的結(jié)果基本一致,但在噴管的擴(kuò)張段內(nèi)氣相馬赫數(shù)存在偏差。因?yàn)樵诎l(fā)動(dòng)機(jī)擴(kuò)張段內(nèi)顆粒在軸線上的積聚數(shù)量相差較大,則顆粒對(duì)氣相的膨脹抑制作用很強(qiáng),同時(shí)顆粒對(duì)氣相的對(duì)流換熱作用也很強(qiáng)。5種平均直徑計(jì)算方法得到的結(jié)果均偏大,其中線性平均法、表面積平均法、體積平均法、體積表面積平均法和質(zhì)量平均法的計(jì)算結(jié)果偏差依次減小。由圖4可知,在發(fā)動(dòng)機(jī)燃燒室和噴管收斂段內(nèi)5種計(jì)算方法得到的軸線氣相溫度基本一致,偏差較大的位置在噴管的擴(kuò)張段內(nèi),其中5種方法的計(jì)算結(jié)果在噴管的擴(kuò)張段均偏小。因?yàn)轭w粒直徑較大時(shí)顆粒的隨流性差,向發(fā)動(dòng)機(jī)軸線積聚的數(shù)量越多,對(duì)氣相的傳熱量也越大。質(zhì)量平均法計(jì)算結(jié)果與多尺寸分布結(jié)果最接近,其次是體積表面積平均法和體積平均法,線性平均計(jì)算結(jié)果偏差最大。通過氣相馬赫數(shù)與溫度的比較分析得,5種方法中質(zhì)量平均法計(jì)算效果更好,其次是體積表面積平均法,模擬效果最差的是線性平均法。

推力F是固體火箭發(fā)動(dòng)機(jī)的主要性能參數(shù)[10]。表3給出了多尺寸分布與5種顆粒平均直徑計(jì)算方法得到的發(fā)動(dòng)機(jī)推力。由表3可知,多尺寸分布計(jì)算得到的推力介于體積表面積平均法與質(zhì)量平均法之間,其中質(zhì)量平均法的計(jì)算結(jié)果偏小,而體積表面積平均法的計(jì)算結(jié)果偏大,因此更好的顆粒平均直徑值介于以上兩種方法計(jì)算結(jié)果之間。可以對(duì)以上兩種方法得到的平均直徑進(jìn)行求算術(shù)平均值能夠使計(jì)算結(jié)果更加接近多尺寸分布結(jié)果。5種計(jì)算方法中線性平均法的計(jì)算結(jié)果偏差最大,因此對(duì)于計(jì)算發(fā)動(dòng)機(jī)推力時(shí)一般不宜采用。

圖4 發(fā)動(dòng)機(jī)軸線氣相溫度

表3 單尺寸與多尺寸分布計(jì)算推力比較

推力系數(shù)Cf是表征噴管性能的參數(shù)[10],其中Cf=F/(PcAt)。表4給出了多尺寸分布與5種顆粒平均直徑計(jì)算方法得到的發(fā)動(dòng)機(jī)推力系數(shù)。由表4可知,多尺寸分布計(jì)算得到的推力系數(shù)介于體積表面積平均法與質(zhì)量平均法之間,其中質(zhì)量平均法的計(jì)算結(jié)果偏小,而體積表面積平均法的計(jì)算結(jié)果偏大,因此更好的顆粒平均直徑大小介于以上兩種方法計(jì)算結(jié)果之間?？梢詫?duì)以上兩種方法得到的平均直徑進(jìn)行求算術(shù)平均值能夠使計(jì)算結(jié)果更加接近多尺寸分布結(jié)果。5種計(jì)算方法中線性平均法的計(jì)算結(jié)果偏差最大,因此對(duì)于計(jì)算發(fā)動(dòng)機(jī)推力系數(shù)時(shí)一般不宜采用。

表4 單尺寸與多尺寸分布計(jì)算推力系數(shù)比較

5 結(jié)論

通過采用5種不同的單尺寸分布與多尺寸分布對(duì)固體火箭發(fā)動(dòng)機(jī)內(nèi)的兩相流場(chǎng)進(jìn)行模擬研究,得到了5種不同的單尺寸分布與多尺寸分布計(jì)算得到的發(fā)動(dòng)機(jī)的內(nèi)流場(chǎng)、推力與推力系數(shù)。模擬結(jié)果表明,發(fā)動(dòng)機(jī)軸線上的氣相馬赫數(shù)與溫度質(zhì)量平均法的計(jì)算結(jié)果最接近多尺寸分布,其次是體積表面積平均法,偏差最大是線性平均法。多尺寸分布計(jì)算得到的發(fā)動(dòng)機(jī)推力與推力系數(shù)均介于體積表面積平均法與質(zhì)量平均法得到的單尺寸分布計(jì)算結(jié)果之間,質(zhì)量平均法得到的結(jié)果均偏小,體積表面積平均法得到的結(jié)果均偏大,最好的平均直徑值介于以上兩種平均直徑之間。其中線性平均法計(jì)算結(jié)果偏差最大,因此一般在計(jì)算發(fā)動(dòng)機(jī)性能參數(shù)時(shí)不宜采用。

[1] Jun Sung Park, Seung, Wook Baek. Interaction of a moving shock wave with a two-phase reacting medium [J]. Heat and Mass Transfer, 2003, 46: 4717-4732.

[2] Sachdev J S, Groth C P T, Gottlieb J J. A parallel solution-adaptive scheme for predicting multi-phase core flows in solid propellant rocket motors [C]∥St. John’s Newfoundland, 13th Annual Conference of CFD Society of Canada, 2005.

[3] SH AFQ AT Wahab, XIE Kan, LIU Yu. Numerical simulation of multi-phase combustion flow in solid rocket motors with metalized propellant [J]. 航空動(dòng)力學(xué)報(bào), 2009, 24(7): 1654-1660.

[4] 賀征, 郜冶, 顧璇. 錐柱型裝藥固體火箭發(fā)動(dòng)機(jī)兩相內(nèi)流場(chǎng)中顆粒運(yùn)動(dòng)的數(shù)值模擬 [J]. 宇航學(xué)報(bào), 2005, 26(3): 354-357.

[5] 張志峰, 劉洋, 蔡體敏. 隨機(jī)顆粒軌道模型在長(zhǎng)尾噴管發(fā)動(dòng)機(jī)流場(chǎng)計(jì)算中的應(yīng)用 [J]. 固體火箭技術(shù), 2007, 30(5): 376-380.

[6] 劉靜, 徐旭. 隨機(jī)軌道模型在噴管兩相流計(jì)算中的應(yīng)用 [J]. 固體火箭技術(shù), 2006, 29(5): 333-353.

[7] 方丁酉. 兩相流體力學(xué) [M]. 長(zhǎng)沙: 國(guó)防科技大學(xué)出版社, 1988: 1-68.

[8] 王增輝, 陳軍濤, 陳林泉. Al2O3尺寸分布對(duì)噴管性能預(yù)示的影響研究 [J]. 彈箭與制導(dǎo)學(xué)報(bào), 2008, 28(5): 142-144.

[9] 劉曉俐. 高金屬含量固體火箭發(fā)動(dòng)機(jī)的相關(guān)問題研究 [D]. 哈爾濱: 哈爾濱工程大學(xué), 2006.

[10] 董師顏, 張兆良. 固體火箭發(fā)動(dòng)機(jī)原理 [M]. 北京: 北京理工大學(xué)出版社, 1996.

LIU Zai1,LI Yanchen2

(1 College of Aerospace and Civil Engineering, Harbin Engineering University, Harbin 150001, China; 2 Shanghai Institute of Spacecraft Equipment, Shanghai 200240, China)

To determine more reasonable calculation method for average particle diameter, the differences of internal flow field and performance of solid rocket motor in multi-size and sin gle-size distribution conditions were separately simulated. Numerical simulation usin gk-εturbulence under high Reynolds number and particle trajectory model was introduced. It shows that the results calculated by multi-size distribution are different from that of sin gle size distribution obtained through five calculation methods for average particle diameter. Mass average method and averaging method in volume and surface area are better to calculate particle average diameter, while linear average method should not be used.

particle trajectory model; mass average method; averaging method in volume and surface area; linear average method

2014-07-11

劉仔(1991-),男,四川鄰水人,學(xué)士,研究方向:固體火箭發(fā)動(dòng)機(jī)內(nèi)的流動(dòng)、傳熱與燃燒。

V435

A

Simulation Research of Five Average Particle Diameter Calculation Methods