形狀及旋轉(zhuǎn)角度對非球鋁顆粒受力的影響分析①

賀 征，劉叢林，顧 璇，郜 冶

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

0 引言

含鋁金屬顆粒的氣固多相流在固體火箭發(fā)動機工作過程中很普遍。其中，離散相顆粒的動力學(xué)特性是不可忽視的重要問題。當(dāng)金屬顆粒在流場中發(fā)生相變?nèi)紵龝r，部分凝相將回落到顆粒表面，令其轉(zhuǎn)變?yōu)椴灰?guī)則的非球體［1－2］，導(dǎo)致動力學(xué)特性發(fā)生一系列改變，進而影響整個多相流場的流動特性。經(jīng)典的顆粒受力理論是以規(guī)則球形為基礎(chǔ)進行推導(dǎo)的，在實際應(yīng)用中不可避免地會產(chǎn)生一定偏差［3］。尤其當(dāng)來流與非球顆粒主軸呈一定角度時，其受力更加復(fù)雜。

本文采用數(shù)值模擬方法，針對固體火箭發(fā)動機內(nèi)鋁金屬顆粒在燃燒過程中因形狀變化而引起受力改變的現(xiàn)象進行研究，對工程應(yīng)用有一定的參考作用。

1 物理模型

鋁顆粒初始進入發(fā)動機時，其表面就包裹著一層氧化物，顆粒形狀已經(jīng)發(fā)生變化［4］，受力狀態(tài)發(fā)生一定改變。Merrill［5］研究了鋁顆粒在固體火箭發(fā)動機環(huán)境中的燃燒過程，對Al2O3的生成過程進行了推導(dǎo)和計算。結(jié)果表明，顆粒在相變?nèi)紵^程中，其外形不斷變化，主要以Al顆粒和Al2O3顆粒兩球相切的形式出現(xiàn)，如圖1所示。

以Merrill的研究結(jié)果為基礎(chǔ)，針對初始半徑為100 μm的顆粒，分析其在燃燒過程中不同時刻下形狀發(fā)生變化，而且相對于來流旋轉(zhuǎn)到一定角度時的受力情況。記來流的方向為x，顆粒主軸與來流方向所夾銳角為θ。如圖2所示，當(dāng)顆粒主軸平行于來流方向時，θ=0°;顆粒主軸垂直于來流方向時，θ=90°。為詳細反映顆粒的受力情況，對應(yīng)于t=0、30、60、90 ms時，分別取 θ=0°、30°、60°和 90°的工況進行研究。

利用Fluent軟件對流場中顆粒的受力狀態(tài)進行分析，為準(zhǔn)確計算顆粒附近流場的流動狀態(tài)，同時兼顧計算網(wǎng)格的數(shù)量，計算區(qū)域中取平行于來流方向的長度為顆粒半徑的30倍，垂直于來流方向的寬度為顆粒半徑的20倍。采用非正交網(wǎng)格進行劃分，對顆粒表面附件區(qū)域網(wǎng)格做了加密處理。

圖2 物理模型Fig.2 Physical model

2 控制方程

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

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

顆粒阻力系數(shù)CD由其所受的阻力FD來定義［6］:

其中，F(xiàn)D由兩部分組成，即顆粒表面壓強梯度所產(chǎn)生的壓差阻力及流體粘性所產(chǎn)生的粘性阻力。

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

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

3.1 模型驗證

為驗證各種湍流模型對計算結(jié)果的影響和計算方法的可信性，首先對球形顆粒在1≤Red≤1 000范圍內(nèi)的受力問題進行計算，以尋求最佳計算模型。采用的湍流模型包括realizable k-ε模型、RNG k-ε模型、sdandard k-ε模型和Spalart-Allmadas模型。

文獻［7－8］對不同顆粒雷諾數(shù)情況下的顆粒阻力系數(shù)進行了實驗，文獻［9］采用直接模型方法對Rep在10～300范圍內(nèi)的顆粒阻力系數(shù)進行了計算，以上結(jié)果均可作為參考。分別采用上述文獻中實驗和計算條件，對不同顆粒雷諾數(shù)的情況下的球形顆粒進行了計算驗證。圖3直觀地反映了各模型的計算結(jié)果、標(biāo)準(zhǔn)阻力曲線和文獻實驗值的對比。

圖3 不同湍流模型的計算結(jié)果與文獻結(jié)果的對比Fig.3 Computed results of different model vs reference results

由圖3可見，當(dāng) Rep＜100時，除 Spalart-Allmadas模型外，其余3種模型均能很好地與標(biāo)準(zhǔn)阻力曲線相符。其中，realizable k-ε模型下，各種情況的計算值與標(biāo)準(zhǔn)阻力曲線值的均方差最小，為3.63。當(dāng)Rep＞100時，4種模型的計算值與標(biāo)準(zhǔn)曲線相比，差別都很大，但可較好地符合文獻［10－11］的實驗結(jié)果，這說明所使用方法的可行性和可信性。為兼顧較大范圍內(nèi)的計算，選擇realizable k-ε模型進行顆粒受力分析的數(shù)值模擬更為適宜。

3.2 形狀變化對鋁顆粒受力的影響

經(jīng)典顆粒受力的推導(dǎo)過程均是將流場中的顆粒當(dāng)作球形進行的，而對于燃燒的鋁顆粒，原有的計算模型在實際應(yīng)用中產(chǎn)生一定偏差。

在實際的工程應(yīng)用中，對于非球形不規(guī)則顆粒的數(shù)值計算大多以近似球形顆粒的方法進行處理。在多種處理方法中，以等效直徑應(yīng)用較為廣泛。其中，又以等體積直徑dv應(yīng)用較多。顯然，真實鋁顆粒的受力與等效直徑顆粒受力間將存在很大的差別，為深入分析形狀變化對鋁顆粒受力的影響，暫取顆粒主軸平行于來流方向時，二者之間的差別進行計算，計算工況列于表1中。表1中，rAl、rox分別表示顆粒中鋁的半徑及與其相連的圓形氧化物半徑;req為其相應(yīng)的當(dāng)量半徑;Ap、Apeq分別為顆粒主軸垂直于來流方向的橫截面面積。統(tǒng)計結(jié)果表明，隨顆粒燃燒狀態(tài)的改變，二者相差很大，最大為49.9%，即便在顆粒初始進入流場時，也有1.1%的差別。

顆粒在主軸垂直于來流方向的橫截面變化，必然會引起受力的改變。圖4反映了鋁顆粒在燃燒過程中所受到的阻力系數(shù)變化，以Cd表示，同時也計算了當(dāng)量直徑下顆粒所受到的阻力系數(shù)變化，以Cdeq表示。由理論分析知，隨著燃燒的進行，顆粒直徑越來越小，導(dǎo)致顆粒雷諾數(shù)不斷降低，由此也使得阻力系數(shù)不斷升高。應(yīng)用當(dāng)量直徑計算所得的結(jié)果也符合這一規(guī)律，但以真實顆粒計算時，二者卻存在一定差別。整體而言，真實顆粒所受到的阻力系數(shù)普遍大于當(dāng)量直徑下的顆粒。

表1 計算工況匯總Table 1 Summary of computation case

圖4 非圓顆粒與其當(dāng)量直徑所受流體阻力系數(shù)對比Fig.4 Viscous coefficient of non-spherical particles vs the equavilent radical particals

圖4(a)表明，隨顆粒相變?nèi)紵倪M行，真實顆粒表面所受到的壓差阻力系數(shù)發(fā)生很大變化，且無規(guī)律可循。因為顆粒的形狀在燃燒過程中不斷發(fā)生不規(guī)則的改變，所以其表面壓力分布狀況也表現(xiàn)出很大的波動。圖4(b)表明，真實顆粒表面受到的粘性阻力系數(shù)一直大于當(dāng)量直徑下的顆粒，尤其在燃燒進行至60 ms時，以當(dāng)量直徑計算所得到的粘性阻力系數(shù)僅為真實顆粒的68.9%，當(dāng)顆粒終止燃燒時，前者仍比后者高18%。圖4(c)為顆粒所受到的總阻力系數(shù)對比，其變化趨勢與圖4(b)相似。

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

3.3 旋轉(zhuǎn)角度變化對鋁顆粒受力的影響

在不同時刻、不同來流角度下，顆粒表面周圍的流線分布如圖5所示。隨著顆粒主軸與來流間夾角的增大，流場在顆粒背風(fēng)面所產(chǎn)生的旋渦附著點位置有由頂端逐漸向下移動的趨勢，當(dāng)顆粒軸線與來流呈90°直角時，背風(fēng)面的旋渦則分成兩部分，分別附著于顆粒表面流速和壓力最低處。

隨著顆粒主軸與來流間夾角的增大，流場在顆粒背風(fēng)面所產(chǎn)生的旋渦附著點位置有由頂端逐漸向下移動的趨勢，當(dāng)顆粒主軸與來流呈90°直角時，背風(fēng)面的旋渦則分成兩部分，分別附著于顆粒表面流速和壓力最低處。各時刻下，顆粒表面附近的壓力分布沒有明顯區(qū)別，由顆粒所引起的流場壓力變化區(qū)域在迎風(fēng)面不大于3倍的顆粒半徑，背風(fēng)面不超過1.5倍的顆粒半徑，且當(dāng)顆粒與來流夾角θ＜30°時，背風(fēng)面的壓力變化范圍更小。

各時刻下，顆粒尾流的寬度與長度因顆粒主軸與來流夾角不同而有所改變。當(dāng)顆粒主軸不平行或垂直于來流時，即θ介于0°～90°之間時，流場在顆粒背風(fēng)面所形成尾流的寬度有所增加，這種現(xiàn)象在顆粒開始燃燒一段時間后，顆粒外形嚴(yán)重偏離球形時更加明顯。尤其當(dāng)顆粒燃燒60 ms后，θ=60°時，因顆粒主軸與來流相垂直的面積相對較大，其周圍所形成的尾流最寬，但相比于其他工況，此時尾流的長度有所減小，氣流流經(jīng)顆粒后，在背風(fēng)面較短的一段距離處，便恢復(fù)到流場的平均水平。

圖5(a)表明，顆粒初始燃燒時，其形狀沒有發(fā)生明顯改變，顆粒主軸與來流間角度的變化對流場的壓力分布和速度分布沒有產(chǎn)生很大影響，但很大程度上改變了流場漩渦所發(fā)生的位置。圖5(b)表明，顆粒開始燃燒30 ms后，其氧化物所形成的凝相約為金屬顆粒體積的18%，顆粒外形已偏離球體。因此，當(dāng)顆粒主軸與來流間夾角發(fā)生改變時，漩渦位置發(fā)生明顯變化。圖5(c)表明，顆粒開始燃燒60 ms后，其氧化物所形成的凝相體積與原金屬顆粒體積相當(dāng)，顆粒外形接近橢球體，流場分布因顆粒與來流間夾角不同而發(fā)生明顯改變。當(dāng)θ由30°增加到90°時，流場漩渦位置發(fā)生明顯改變，且當(dāng)θ=90°時，漩渦由1個變成2個，對稱分布于顆粒背風(fēng)面。圖5(d)反映出的流場變化規(guī)律與圖5(c)類似，但比較可知，流線分布狀態(tài)的改變程度弱于t=60 ms時的工況。

圖5 顆粒周圍流線分布Fig.5 Stream lines around the particle

3.4 顆粒受力隨燃燒時間的變化

因為非球形顆粒不具有各向同性的性質(zhì)，所以當(dāng)顆粒主軸線與來流呈一定角度時，顆粒外表面的受力不再具有對稱性，除受平行于來流方向的x向阻力外，顆粒還會受到垂直于來流方向的y向升力。圖6和圖7分別反映了顆粒角度變化時，所受x方向阻力系數(shù)和y方向升力系數(shù)的變化規(guī)律。

圖6 x方向顆粒阻力系數(shù)隨來流角度的變化Fig.6 Particle＇s drag coefficient in x direction

圖7 y方向上顆粒升力系數(shù)隨來流角度變化Fig.7 Particle＇s lift coefficient in y direction

隨燃燒的進行，顆粒體積越來越小，使得顆粒雷諾數(shù)逐漸降低。因此，顆粒所受到的阻力系數(shù)整體上是隨時間增大的。但同一時刻下，若顆粒在流場中的擺放位置發(fā)生改變，那么它所受到的阻力系數(shù)也將有所不同。由圖5可見，隨顆粒與來流間角度的增加，其受到的壓差阻力系數(shù)呈不斷增大的趨勢，見圖6(a)，而粘性阻力系數(shù)則不斷減小，見圖6(b)。當(dāng)θ=30°時，隨燃燒的進行，顆粒所受的粘性阻力系數(shù)變化較大，至燃燒終止時，顆粒所受到的粘性阻力系數(shù)增長了26%;但顆粒受到的壓差阻力系數(shù)基本保持不變。當(dāng)θ=60°和90°時，顆粒的受力情況與之相反。

在所研究范圍內(nèi)，顆粒在流場中受到的壓差阻力系數(shù)占據(jù)更重要的地位。因此，作用于顆粒的總阻力系數(shù)在整體上是隨夾角的增加而不斷增加的，見圖6(c)。隨顆粒燃燒活動的進行，這種差別愈加明顯，當(dāng)顆粒燃燒60 ms后，若顆粒主軸垂直于來流方向(θ=90°)，其受到的阻力系數(shù)比 θ=30°時增加6.6%，而當(dāng)顆粒終止燃燒時，在上述2種情形下，顆粒所受到的阻力系數(shù)也仍有5.5%的差別。

當(dāng)顆粒不再平行于來流時，顆粒主軸將受到垂直于來流方向上的力。圖7統(tǒng)計了不同工況下，顆粒在y方向上所受到的升力系數(shù)變化規(guī)律。除θ=90°外，顆粒在y方向上所受的升力系數(shù)一般為負值，說明顆粒受到指向y軸負方向的力。由圖7可見，此方向上顆粒受力的變化情況較復(fù)雜，即便在同一時刻下，其值也不是隨顆粒與來流間角度的變化呈現(xiàn)出單調(diào)增加或減少的趨勢。當(dāng)θ=30°和60°時，隨燃燒的進行，顆粒在y方向上所受的升力系數(shù)表現(xiàn)出先增加、后降低的趨勢，但均保持為負值。當(dāng)θ=90°時，顆粒在y方向基本不受粘性力作用，各時刻下的粘性力系數(shù)接近0，但隨燃燒的進行，顆粒所受到的壓差阻力系數(shù)與總阻力系數(shù)呈不斷升高的趨勢。燃燒初始時，顆粒的阻力系數(shù)為負值，即顆粒受到指向y負方向的力的作用，但當(dāng)t=60 ms后，由于凝相的體積超過了原顆粒，所以整體顆粒的受力方向發(fā)生了變化，由y的負向轉(zhuǎn)變?yōu)檎?，升力系?shù)變?yōu)檎怠?/p>

4 結(jié)論

(1)隨著非球顆粒主軸與來流間夾角的增大，流場在顆粒背風(fēng)面所產(chǎn)生的旋渦附著點位置有由頂端逐漸向下移動的趨勢。當(dāng)顆粒主軸與來流呈90°直角時，背風(fēng)面的旋渦則分成兩部分，分別附著于顆粒表面流速和壓力最低處。當(dāng)顆粒主軸與來流呈一定角度時，顆粒外表面的受力不再具有對稱性，將同時受到x方向阻力和y方向升力的作用。

(2)在金屬顆粒的燃燒初期，顆粒雷諾數(shù)較大，顆粒在流場中的旋轉(zhuǎn)角度對其受力情況影響不大。隨著燃燒的不斷進行，顆粒雷諾數(shù)逐漸減小，當(dāng)顆粒主軸與來流呈不同角度時，其受力情況有很大改變。尤其在x方向上，顆粒所受阻力系數(shù)差值最大為5%。其中，粘性阻力系數(shù)差值最大為18.9%。相對于阻力系數(shù)而言，顆粒受到的升力系數(shù)量級很小，最大值不足阻力系數(shù)的20%。

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

［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.

［4］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.

［5］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.

［6］Robert Geisler.A global view of the use of aluminum fuel in solid rocket motors［R］.AIAA 2002-3748

［7］Merrill K King.Aluminum combustion in a solid rocket motor environment［J］.Proceedings of the combustion institute，2009，32(2):2107-2114.

［8］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.

［9］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.

［10］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.

［11］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.