賀 征, 劉叢林，李 卓，顧 璇，郜 冶

(哈爾濱工程大學(xué) 航天與建筑工程學(xué)院，哈爾濱 150001)

?

金屬顆粒燃燒過程表面曳力變化的數(shù)值研究

賀 征, 劉叢林，李 卓，顧 璇，郜 冶

(哈爾濱工程大學(xué) 航天與建筑工程學(xué)院，哈爾濱 150001)

為了考察固體火箭發(fā)動機(jī)內(nèi)金屬顆粒在生長過程中所產(chǎn)生的非球顆粒受力問題，采用數(shù)值模擬方法，在通過和文獻(xiàn)及實(shí)驗(yàn)數(shù)據(jù)對比驗(yàn)證確定最佳計(jì)算模型的基礎(chǔ)上，對初始半徑為100 μm的鋁顆粒處于不同相變?nèi)紵A段時(shí)所受到的曳力進(jìn)行對比分析。結(jié)果表明，在多相流場中，非球顆粒表面可能存在兩處壓力為零的點(diǎn)，曳力系數(shù)普遍大于其當(dāng)量直徑顆粒的計(jì)算結(jié)果，當(dāng)顆粒外形嚴(yán)重偏離球體時(shí)，經(jīng)驗(yàn)公式計(jì)算所得的曳力系數(shù)有失精準(zhǔn)，需進(jìn)行相應(yīng)修正。

金屬顆粒；非球顆粒；燃燒；當(dāng)量直徑；曳力系數(shù)

0 引言

含金屬顆粒的氣固多相流是一種十分典型的多相流動。其中，分散相顆粒的動力學(xué)特性是工程應(yīng)用中不可忽視的重要問題。當(dāng)金屬顆粒在流場中發(fā)生相變?nèi)紵龝r(shí)，部分凝相將回落到顆粒表面，令其轉(zhuǎn)變?yōu)椴灰?guī)則的非球體[1-2]，導(dǎo)致動力學(xué)特性發(fā)生一系列改變，進(jìn)而影響整個(gè)多相流場的流動特性。而經(jīng)典的顆粒受力理論是以規(guī)則球形為基礎(chǔ)進(jìn)行推導(dǎo)的，在實(shí)際應(yīng)用中，不可避免地會產(chǎn)生一定偏差[3]。采用數(shù)值模擬方法，針對固體火箭發(fā)動機(jī)內(nèi)，金屬顆粒在生長過程中，因形狀變化而引起受力改變的現(xiàn)象進(jìn)行研究，對工程應(yīng)用有一定參考作用。

早期的研究集中于球形液滴，并已經(jīng)形成了成熟的經(jīng)驗(yàn)關(guān)系式[4-5]。相對而言，對于非球形顆粒的研究開展的較少，Komar[6]、Loth[7]開展了圓柱體的阻力系數(shù)計(jì)算，并推導(dǎo)其阻力系數(shù)的變化。Sabine[8]等研究了若干顆粒結(jié)合時(shí)的阻力系數(shù)。Youngho and Changhoon[9]模擬計(jì)算了變形顆粒的的阻力和升力，發(fā)現(xiàn)其數(shù)值與顆粒的形變有直接關(guān)系。

目前，對于非球形不規(guī)則顆粒的數(shù)值計(jì)算，大多采用近似球形的處理方法[10]。這種近似過程的實(shí)現(xiàn)主要通過2個(gè)途徑進(jìn)行——引入等效直徑或引入球形度。其中，前者應(yīng)用的較為廣泛，最常用的等效直徑是等體積直徑dv[11]。引入dv后，認(rèn)為非球顆粒與其等效直徑相同的球形顆粒受力相同，文中對這2種計(jì)算結(jié)果進(jìn)行了對比分析。

1 物理模型

鋁顆粒初始進(jìn)入發(fā)動機(jī)時(shí)，其表面將迅速被氧化物所包裹，顆粒形狀已經(jīng)發(fā)生變化[12]，顆粒表面積以及同來流方向相垂直的橫截面面積與球形顆粒不再相同，受力狀態(tài)發(fā)生一定改變。

Merrill[13]研究了鋁顆粒在固體火箭發(fā)動機(jī)環(huán)境中的燃燒過程，對Al2O3的生成過程做了仔細(xì)的推導(dǎo)和計(jì)算。相比于傳統(tǒng)的R2定律，Merrill所得結(jié)果更加符合實(shí)際情況。以Merrill的研究結(jié)果為基礎(chǔ)，對初始半徑為100 μm的鋁顆粒在熱燃?xì)饬髦刑幱诓煌嘧內(nèi)紵A段時(shí)的受力問題進(jìn)行了分析，計(jì)算工況列于表1中。

表1 計(jì)算工況匯總Table1 Summary of computation cases

其中，t表示顆粒進(jìn)入發(fā)動機(jī)燃燒室中的時(shí)間，RAl、Rox分別表示顆粒中鋁的半徑和與其相連的球形氧化物半徑，Req為其相應(yīng)的等體積半徑，Ap、Apeq分別代表真實(shí)顆粒與等體積顆粒垂直于來流方向的橫截面面積。統(tǒng)計(jì)表明，隨顆粒燃燒狀態(tài)的改變，等體積顆粒的橫截面面積偏離真實(shí)顆粒的水平(Δ%)較大，最大相差49.9%，即便在顆粒初始進(jìn)入流場時(shí)，也有1.1%的差別。

利用Fluent軟件，對流場中顆粒的受力狀態(tài)進(jìn)行分析。顆粒半徑r取為100 μm，為準(zhǔn)確計(jì)算顆粒附近流場的流動狀態(tài)，同時(shí)兼顧計(jì)算網(wǎng)格的數(shù)量，計(jì)算區(qū)域中，取平行于來流方向的長l為30r，即3 mm，垂直于來流方向的寬d為20r，即2 mm。顆粒置于流場中，靠近來流方向，中心與入口邊界相距l(xiāng)1為10r，即1 mm。采用非正交網(wǎng)格進(jìn)行劃分，顆粒的近壁面處做加密處理，最小網(wǎng)格尺寸為5 μm，最大為100 μm，總網(wǎng)格數(shù)約為22萬，如圖1所示。

2 控制方程

假設(shè)來流是穩(wěn)態(tài)的，其物理性質(zhì)在運(yùn)動過程中不發(fā)生變化，計(jì)算的控制方程為

▽·U=0

(1)

▽·(ρUU)=-▽p+μ▽2U

(2)

為簡化模擬條件，不考慮氣固兩相間的傳熱，取空氣為介質(zhì)，密度為1.225kg/m3，粘性系數(shù)為1.789×10-5kg/(m·s)，流場入口速度取發(fā)動機(jī)燃燒室環(huán)境下的均值，為20m/s，流場出口為自由邊界。

顆粒曳力系數(shù)Cd由顆粒所受到的阻力FD來定義[14]:

(3)

顆粒所受的阻力FD由2部分組成：顆粒表面壓強(qiáng)梯度所產(chǎn)生的壓力阻力，以及流體粘性所產(chǎn)生的粘性阻力。因此，顆粒的曳力系數(shù)也分為2部分——壓力曳力系數(shù)和粘性曳力系數(shù)。

隨顆粒相變?nèi)紵倪M(jìn)行，尺寸不斷改變，顆粒雷諾數(shù)Rep隨之發(fā)生變化，以當(dāng)量直徑計(jì)，對應(yīng)于t=0、30、60、90 ms時(shí)，Rep分別為274.0、116.2、10.4、9.6。

(a)三維計(jì)算模型

(b)模型網(wǎng)格分布

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

3.1 模型驗(yàn)證

為驗(yàn)證各種湍流模型對計(jì)算結(jié)果的影響，首先對顆粒在1≤Rep≤1 000范圍內(nèi)的受力問題進(jìn)行計(jì)算，以尋求最佳計(jì)算模型。采用的湍流模型包括realizablek-ε模型、RNGk-ε模型、standardk-ε模型和Spalart-Allmadas模型。

驗(yàn)證工作共包括32種工況，顆粒雷諾數(shù)Rep分別取為1、10、30、50、100、300、500及1 000，流場工質(zhì)均為空氣。文獻(xiàn)[15]與[16]分別對Rep在18.7～87.3以及10～809.7范圍內(nèi)的顆粒曳力系數(shù)做了實(shí)驗(yàn)，文獻(xiàn)[17]采用直接模型方法對Rep在10～300范圍內(nèi)的顆粒曳力系數(shù)進(jìn)行了計(jì)算，以上結(jié)果均可作為參考。對于Rep=1和1 000的工況，尚無實(shí)驗(yàn)可詢，取標(biāo)準(zhǔn)阻力曲線值相對比。

圖2直觀地反映了各模型的計(jì)算結(jié)果與標(biāo)準(zhǔn)阻力曲線和文獻(xiàn)實(shí)驗(yàn)值的對比。當(dāng)顆粒雷諾數(shù)Rep<100時(shí)，除Spalart-Allmadas模型外，其余3種模型均能很好地與標(biāo)準(zhǔn)阻力曲線相符。其中，realizablek-ε模型下各種情況的計(jì)算值與標(biāo)準(zhǔn)阻力曲線值的均方差為3.63，RNGk-ε模型的均方差為3.77、標(biāo)準(zhǔn)k-ε模型為4.05，而Spalart-Allmadas模型的均方差則高達(dá)14.28。當(dāng)顆粒雷諾數(shù)Rep較高(>100)時(shí)，4種模型的計(jì)算值與標(biāo)準(zhǔn)曲線相比，差別都很大，但可較好符合文獻(xiàn)[2-3]的實(shí)驗(yàn)結(jié)果。為兼顧較大范圍內(nèi)的計(jì)算，選擇Realizablek-ε模型進(jìn)行顆粒受力分析的數(shù)值模擬更為適宜。

圖2 不同湍流模型的計(jì)算結(jié)果與文獻(xiàn)結(jié)果的對比Fig.2 Computed results of different models vsreference results

3.2 顆粒周圍流場壓力與速度分布對比分析

在多相流場中，顆粒會隨著流動發(fā)生位置偏轉(zhuǎn)，為分析真實(shí)顆粒與球形當(dāng)量徑顆粒的差別，取典型狀態(tài)——顆粒在垂直來流方向投影面最大時(shí)進(jìn)行計(jì)算。不同時(shí)刻下，真實(shí)顆粒與當(dāng)量直徑顆粒流場附近的壓力與速度分布如圖3、圖4所示。在迎風(fēng)方向上，顆粒表面附近的氣流速度由流場入口沿顆粒中心軸線方向迅速降低，在顆粒表面處出現(xiàn)駐點(diǎn)，同時(shí)當(dāng)?shù)貕毫_(dá)到最高值。顆粒背風(fēng)面處也存在著一個(gè)速度為零的區(qū)域，同時(shí)當(dāng)?shù)貕毫档阶畹椭?。流場中顆粒背風(fēng)域形成長長的流動尾跡表明，迎風(fēng)面的氣流速度衰減梯度明顯大于背風(fēng)面的速度增長梯度。

觀察顆粒表面的氣流變化可發(fā)現(xiàn)，氣流由駐點(diǎn)的零速度向滯止點(diǎn)的零速度轉(zhuǎn)變過程是氣流速度沿顆粒表面先增加，至顆粒與來流相垂直的軸線中心處，速度達(dá)到最大值，且此值大于流場的平均速度；之后，因顆粒背風(fēng)面存在渦流區(qū)，氣流速度逐漸降低，至顆粒平行于來流的軸線中心處，再次降為0。在所計(jì)算的4種雷諾數(shù)下，顆粒迎風(fēng)面與背風(fēng)面的壓差達(dá)550 Pa左右，顆粒表面處的最大速度值比流場平均速度高10%。隨著顆粒在流場中相變?nèi)紵倪M(jìn)行，當(dāng)顆粒尺寸減小時(shí)，顆粒背風(fēng)域所形成的流動尾跡隨之變短，同時(shí)顆粒附近的壓力變化區(qū)(以-50～50 Pa計(jì))半徑減小。

(a)當(dāng)量直徑顆粒

(b)真實(shí)顆粒

對比圖3、圖4中當(dāng)量直徑顆粒流場與真實(shí)顆粒流場流動情況可知，當(dāng)t=0 ms，顆粒剛剛開始燃燒時(shí)，真實(shí)顆粒軸向長度大于當(dāng)量直徑下的顆粒，顆粒背風(fēng)域所形成的流動尾跡較長。但此后，因當(dāng)量直徑下顆粒背風(fēng)域所引起的渦流較強(qiáng)，氣流需要較長的路徑來恢復(fù)原來的速度。所以，其余3種工況下，真實(shí)顆粒背風(fēng)域所形成的流動尾跡相對較短。因?yàn)轭w粒在相變過程中呈現(xiàn)出不規(guī)則的幾何形體，真實(shí)顆粒表面附近的壓力分布出現(xiàn)較大波動，尤其在t=60 ms和t=90 ms時(shí)，顆粒外形明顯分成兩個(gè)區(qū)域。因此，在顆粒的表面處存在兩處壓力為0的點(diǎn)，這與當(dāng)量直徑下顆粒表面處的壓力分布有很大區(qū)別。

(a)當(dāng)量直徑顆粒

(b)真實(shí)顆粒

3.3 顆粒曳力對比分析

圖5反映了鋁顆粒在燃燒過程中其曳力系數(shù)的變化規(guī)律，以Cd表示；同時(shí)，也計(jì)算了當(dāng)量直徑下顆粒曳力系數(shù)變化，以Cdeq表示。由理論分析知，隨著燃燒的進(jìn)行，顆粒直徑越來越小，將導(dǎo)致顆粒雷諾數(shù)不斷降低，由此也會使得曳力系數(shù)不斷升高。應(yīng)用當(dāng)量直徑計(jì)算所得的結(jié)果也符合這一規(guī)律，但以真實(shí)顆粒計(jì)算時(shí)，二者卻存在一定差別，整體而言，真實(shí)顆粒的曳力系數(shù)普遍大于當(dāng)量直徑下的顆粒。圖5(a)表明，隨顆粒相變?nèi)紵倪M(jìn)行，真實(shí)顆粒表面的壓力曳力系數(shù)發(fā)生很大變化，且無規(guī)律可循。因?yàn)轭w粒的形狀在燃燒過程中不斷發(fā)生不規(guī)則的改變，所以其表面壓力分布狀況也表現(xiàn)出很大的波動。圖5(b)表明，真實(shí)顆粒表面的粘性曳力系數(shù)一直大于當(dāng)量直徑下的顆粒，尤其當(dāng)燃燒進(jìn)行了60 ms后，以當(dāng)量直徑計(jì)算所得到的粘性曳力系數(shù)僅為真實(shí)顆粒的68.9%，當(dāng)顆粒終止燃燒時(shí)，前者仍比后者高18%。圖5(c)為顆粒的總曳力系數(shù)對比，其變化趨勢與圖5(b)相似。

從所計(jì)算的4種工況而言，以當(dāng)量直徑計(jì)算所得的壓力曳力系數(shù)略高于真實(shí)顆粒，而粘性曳力系數(shù)則低于真實(shí)顆粒。所以，總曳力系數(shù)的偏差并不明顯。但在顆粒燃燒進(jìn)行到60 ms時(shí)，因外形已嚴(yán)重偏離球體，顆粒表面壓力曳力系數(shù)與粘性曳力系數(shù)均大于當(dāng)量直徑下的計(jì)算值，其總曳力系數(shù)與經(jīng)驗(yàn)公式計(jì)算所得的結(jié)果大不相同。此時(shí)，標(biāo)準(zhǔn)阻力曲線已不再適用，需對其進(jìn)行相應(yīng)的修正。

(a)壓力曳力系數(shù)對比

(b)粘性曳力系數(shù)對比

(c)總曳力系數(shù)對比

4 結(jié)論

因?yàn)榻饘兕w粒在相變過程中呈現(xiàn)出不規(guī)則的幾何形體，真實(shí)顆粒表面附近的壓力分布與周圍流場速度分布同球形顆粒有較大區(qū)別，其曳力系數(shù)亦有較大的變化。

(1)隨著顆粒在流場中相變?nèi)紵倪M(jìn)行，當(dāng)顆粒尺寸減小時(shí)，顆粒背風(fēng)域所形成的流動尾跡隨之變短；同時(shí)，顆粒附近的壓力變化區(qū)半徑減小。

(2)在真實(shí)非球體顆粒的表面處，可能存在兩處壓力為0的點(diǎn)。顆粒直徑的減小將導(dǎo)致其雷諾數(shù)降低，曳力系數(shù)不斷升高。整體而言，真實(shí)顆粒的曳力系數(shù)普遍大于其當(dāng)量直徑顆粒。其中，壓力曳力系數(shù)變化很大，但無規(guī)律可循。

(3)粘性曳力系數(shù)普遍大于當(dāng)量直徑下的顆粒，當(dāng)燃燒進(jìn)行60 ms后，其值為以當(dāng)量直徑計(jì)算所得結(jié)果的1.46倍。當(dāng)顆粒終止燃燒時(shí)，前者仍比后者高出18%。

(4)當(dāng)顆粒外形嚴(yán)重偏離球體時(shí)，經(jīng)驗(yàn)公式計(jì)算所得的阻力系數(shù)有失精準(zhǔn)。此時(shí)，標(biāo)準(zhǔn)阻力曲線已不再適用，可考慮對當(dāng)量直徑的計(jì)算結(jié)果進(jìn)行系數(shù)修正。

[1] 方丁酉.兩相流體力學(xué)[D].長沙:國防科技大學(xué)出版社，1988:55-58.

[2] Olsen S E,Beckstead M W.Burn time measurements of single aluminum particles in steam and CO2mixtures[J].Journal of Propulsion and Power,1996,12(4):662-671.

[3] 林建中.超常顆粒多相流體動力學(xué)——圓柱狀顆粒兩相流[M].北京:科學(xué)出版社,2008:85-86.

[4] Khan A R,Richardson J F.The resistance to motion of a solid sphere in a fluid[J].Chem.Eng.Commun,1987,62:135-150.

[5] Hartman M,Yates J G.Free-fall of solid particles through fluids[J].Collet.Czechoslov.Chem.Commun,1993,58(5):961-982.

[6] Komar P D.Settling velocities of circular cylinders at low Reynolds number[J].J.Geol.,1980,48:327-328.

[7] Loth E.Drag of non-spherical solid particles of regular and irregular shape[J].Powder Technology,2008,182:342-353.

[8] Sabine Tran-cng,Michael Gay,efastathios E Michaelides.Drag coefficients of irregular shaped particles[J].Powder Technology,2004,139:21-32.

[9] Youngho Suh,Changhoon Lee.A numerical method for the calculation of drag and lift of a deformable droplet in shear flow[J].Journal of Computational Physics,2013,241:35-57.

[10] Lu Hui-lin,Liu Wen-tie,Zhao Guang-bo.Computational modeling of dense gas particle flow in a pipe:kinetic theory approach of granular flow[J].Journal of Chemical Industry and Engineering (China),2005,5(1):31-38.

[11] Gan Lin,Xu Mao-sheng,zhu Bing-chen.Heat transfer parameters of packed bed with ring pellet[J].Journal of Chemical Industry and Engineering (China),2000,5(6):778-783.

[12] Robert Geisler.A global view of the use of aluminum fuel in solid rocket motors[R].AIAA 2002-3748.

[13] Merrill K King.Aluminum combustion in a solid rocket motor environment[J].Proceedings of The Combustion Institute,2009,32(2):2107-2114.

[14] Ounis H,Ahmadi G,McLaughlin J B.Brownian diffusion of submicrometer particles in the viscous sublayer[J].Journal of Colloid and Interface Science,1991,143(1):266-277.

[15] Unnikrisham A,Chhabra R P.An experimental study of motion of cylinders in Newtonian fluids:wall effects and drag coefficient[J].Can.J.Chem.Eng.,1991,9(9):729-735.

[16] Warnica W D,Renksizbulut M,Strong A B.Drag coefficients of spherical liquid droplets,Part 2:Turbulent gaseous fields[J].Experiments in Fluids.,1995,18:265-276.

[17] By Prosenjit Bagchi,Balachandar S.Inertial and viscous forces on a rigid sphere in straining flows at moderate Reynolds numbers[J].J.Fluid Mech.,2003,481:105-148.

(編輯：崔賢彬)

Numerical study on the change of surface drag force during metal particle combustion

HE Zheng，LIU Cong-lin，LI Zhuo，GU Xuan，GAO Ye
(College of Aerospace and Civil Engineering, Heilongjiang, Harbin 150001, China)

To study the force acted on the non-spherical particle which is produced by burning metal particle in solid ramjet motor,numerical simulation was employed to compare and analyze the drag of aluminum particle whose initial radius is 100 μm in variable phases of the burning process. Before that,the best computational model was determined and validated by comparing with the reference and experimental data.The results show that,in multiphase flow there are two points where the pressure is zero at the non-spherical particle's surface.Compared to the equivalent diameter particle, the real particle's drag coefficient is generally bigger.If the particle is not close to the spherical any more, the drag coefficient computed by the empirical formula is of great roughness,and the formula should be consequentially amended.

metal particle;non-spherical particle;combustion;equivalent diameter;drag coefficient

2014-07-04；

：2014-08-24。

國家自然科學(xué)基金(11372079)；中央高?；究蒲袠I(yè)務(wù)費(fèi)專項(xiàng)基金(HEUCF130203)。

賀征(1978—)，男，博士/副教授，研究方向?yàn)榘l(fā)動機(jī)內(nèi)燃燒與多相流動。E-mail:hezheng1978@163.com

劉叢林(1981—)，女，博士。E-mail:383523445@qq.com

V435

A

1006-2793(2015)04-0492-05

10.7673/j.issn.1006-2793.2015.04.008