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

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

0 引言

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

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

1 物理模型

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

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

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

圖2 物理模型Fig.2 Physical model

2 控制方程

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

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

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

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

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

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

3.1 模型驗(yàn)證

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

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

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

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

3.2 形狀變化對(duì)鋁顆粒受力的影響

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

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

顆粒在主軸垂直于來流方向的橫截面變化，必然會(huì)引起受力的改變。圖4反映了鋁顆粒在燃燒過程中所受到的阻力系數(shù)變化，以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)量直徑下的顆粒。

表1 計(jì)算工況匯總Table 1 Summary of computation case

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

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

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

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

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

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

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

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

3.4 顆粒受力隨燃燒時(shí)間的變化

因?yàn)榉乔蛐晤w粒不具有各向同性的性質(zhì)，所以當(dāng)顆粒主軸線與來流呈一定角度時(shí)，顆粒外表面的受力不再具有對(duì)稱性，除受平行于來流方向的x向阻力外，顆粒還會(huì)受到垂直于來流方向的y向升力。圖6和圖7分別反映了顆粒角度變化時(shí)，所受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

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

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

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

4 結(jié)論

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

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

［1］方丁酉.兩相流體力學(xué)［D］.長(zhǎng)沙:國(guó)防科技大學(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］林建中.超常顆粒多相流體動(dòng)力學(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.