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

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

?

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

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

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

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

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

0 引言

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

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

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

1 物理模型

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

Merrill[13]研究了鋁顆粒在固體火箭發(fā)動機環(huán)境中的燃燒過程，對Al2O3的生成過程做了仔細的推導和計算。相比于傳統(tǒng)的R2定律，Merrill所得結果更加符合實際情況。以Merrill的研究結果為基礎，對初始半徑為100 μm的鋁顆粒在熱燃氣流中處于不同相變?nèi)紵A段時的受力問題進行了分析，計算工況列于表1中。

表1 計算工況匯總Table1 Summary of computation cases

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

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

2 控制方程

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

▽·U=0

(1)

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

(2)

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

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

(3)

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

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

(a)三維計算模型

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

3 計算結果與分析

3.1 模型驗證

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

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

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

圖2 不同湍流模型的計算結果與文獻結果的對比Fig.2 Computed results of different models vsreference results

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

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

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

(a)當量直徑顆粒

(b)真實顆粒

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

(a)當量直徑顆粒

(b)真實顆粒

3.3 顆粒曳力對比分析

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

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

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

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

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

4 結論

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

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

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

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

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

[1] 方丁酉.兩相流體力學[D].長沙:國防科技大學出版社，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] 林建中.超常顆粒多相流體動力學——圓柱狀顆粒兩相流[M].北京:科學出版社,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。

國家自然科學基金(11372079)；中央高?；究蒲袠I(yè)務費專項基金(HEUCF130203)。

賀征(1978—)，男，博士/副教授，研究方向為發(fā)動機內(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