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

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

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

0 引言

含鋁金屬顆粒的氣固多相流在固體火箭發動機工作過程中很普遍。其中，離散相顆粒的動力學特性是不可忽視的重要問題。當金屬顆粒在流場中發生相變燃燒時，部分凝相將回落到顆粒表面，令其轉變為不規則的非球體［1－2］，導致動力學特性發生一系列改變，進而影響整個多相流場的流動特性。經典的顆粒受力理論是以規則球形為基礎進行推導的，在實際應用中不可避免地會產生一定偏差［3］。尤其當來流與非球顆粒主軸呈一定角度時，其受力更加復雜。

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

1 物理模型

鋁顆粒初始進入發動機時，其表面就包裹著一層氧化物，顆粒形狀已經發生變化［4］，受力狀態發生一定改變。Merrill［5］研究了鋁顆粒在固體火箭發動機環境中的燃燒過程，對Al2O3的生成過程進行了推導和計算。結果表明，顆粒在相變燃燒過程中，其外形不斷變化，主要以Al顆粒和Al2O3顆粒兩球相切的形式出現，如圖1所示。

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

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

圖2 物理模型Fig.2 Physical model

2 控制方程

假設來流是穩態的，其物理性質在運動過程中不發生變化，計算的控制方程為

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

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

其中，FD由兩部分組成，即顆粒表面壓強梯度所產生的壓差阻力及流體粘性所產生的粘性阻力。

隨顆粒相變燃燒的進行，尺寸不斷改變，顆粒雷諾數Rep隨之發生變化，以當量直徑計，對應于t=0、30、60、90 ms時，Rep分別為 274.0、116.2、10.4 和 9.6［5］。

3 計算結果與分析

3.1 模型驗證

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

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

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

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

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

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

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

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

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

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

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

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

3.3 旋轉角度變化對鋁顆粒受力的影響

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

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

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

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

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

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

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

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

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

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

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

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

4 結論

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

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

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