劉叢林,郜 冶,賀 征,康文武
(1.哈爾濱工程大學(xué) 航天與建筑工程學(xué)院,哈爾濱 150001;2.北京航空航天大學(xué) 能源與動力工程學(xué)院,北京 100191)
?
相對流速對非球形鋁滴旋轉(zhuǎn)受力的影響①
劉叢林1,郜 冶1,賀 征1,康文武2
(1.哈爾濱工程大學(xué) 航天與建筑工程學(xué)院,哈爾濱 150001;2.北京航空航天大學(xué) 能源與動力工程學(xué)院,北京 100191)
針對固體火箭發(fā)動機(jī)中非球形鋁滴的旋轉(zhuǎn)受力問題進(jìn)行數(shù)值模擬,詳細(xì)分析流速變化對顆粒的受力影響。計算結(jié)果表明,非球顆粒在流場中的升力系數(shù)與無量綱轉(zhuǎn)速之間已經(jīng)不再滿足經(jīng)典的理論關(guān)系式。雷諾數(shù)Re<1時,周向升力系數(shù)存在著交替出現(xiàn)的2組極大值與極小值;升力系數(shù)均為正值,且周向變化較為平緩。隨來流速度增加,阻力系數(shù)變得均衡,系數(shù)值向0值逼近;升力系數(shù)也急速降低。當(dāng)相對流速由0.5 m/s提高到5 m/s時,顆粒所受阻力銳降89%。顆粒的阻力系數(shù)與雷諾數(shù)之間滿足新的關(guān)系式。
非球形鋁滴;相對流速;阻力系數(shù);升力系數(shù)
在含顆粒的多相流計算中,一般視顆粒為離散相,顆粒在流場中的運動軌跡及空間分布取決于作用于其上的力和力矩,對流場中的離散相進(jìn)行受力分析是獲得準(zhǔn)確計算結(jié)果的基礎(chǔ)。尤其固體火箭發(fā)動機(jī)中,金屬鋁滴在流場中相變?nèi)紵龝r,會因部分氧化鋁凝相回落到顆粒表面而轉(zhuǎn)變?yōu)椴灰?guī)則的非球體[1-2],導(dǎo)致動力學(xué)特性發(fā)生一系列改變,從而影響整個多相流場的流動特性。但經(jīng)典的顆粒受力理論是以規(guī)則球形為基礎(chǔ)進(jìn)行推導(dǎo)的,實際應(yīng)用中將不可避免地產(chǎn)生偏差[3]。另外,實際流場中,顆粒的運動狀態(tài)往往是轉(zhuǎn)動伴隨平動同時存在,但為簡化計算,數(shù)值模擬中,前者往往被忽略,理論上講,顆粒在發(fā)動機(jī)中所受阻力與相對流速度呈平方反比關(guān)系。所以,流速變化對顆粒的受力影響很大。固體火箭發(fā)動機(jī)中,在推進(jìn)劑剛開始燃燒時,燃?xì)馑俣容^小,但幾十毫秒后,燃?xì)鈯A帶著顆粒向噴管加速流動,由于慣性作用,顆粒相速度滯后,與主燃?xì)忾g存在著速度差,再加上自身旋轉(zhuǎn)作用,受力狀態(tài)較為復(fù)雜。在固體火箭發(fā)動機(jī)內(nèi)流場中,金屬顆粒的運動特性對整個工作過程起至關(guān)重要的作用,綜合考慮實際環(huán)境中可能存在的顆粒形態(tài)及運動狀態(tài),針對旋轉(zhuǎn)非球形單顆粒的受力問題進(jìn)行詳細(xì)數(shù)值分析,對進(jìn)一步分析顆粒間凝聚、碰撞等運動特性以及顆粒與流體間相互作用有一定的參考作用,對進(jìn)一步認(rèn)清多相流的本質(zhì)特征、完善多相流模型有重要意義。
金屬鋁極易氧化,鋁顆粒存放期間,其表面就包裹著氧化物,初始進(jìn)入發(fā)動機(jī)時,顆粒形狀已經(jīng)發(fā)生變化,不再是規(guī)則的球形[4],受力狀態(tài)與運動特性亦發(fā)生相應(yīng)改變。Merrill[5]研究了鋁顆粒在固體火箭發(fā)動機(jī)環(huán)境中的相變?nèi)紵^程,對其結(jié)構(gòu)變化過程做了仔細(xì)地推導(dǎo)和計算。結(jié)果表明,顆粒將在90 ms內(nèi)完成燃燒過程,其間主要以Al和Al2O3兩球形相切的方式出現(xiàn),隨燃燒的發(fā)展,外形不斷變化,如圖1所示。數(shù)值模擬中,通常以等體積球體代替非球體進(jìn)行受力分析。對Merrill教授的研究結(jié)果進(jìn)行統(tǒng)計可知,隨鋁滴在燃燒室中狀態(tài)的改變,等體積顆粒在來流方向投影的橫截面面積偏離真實顆粒的水平較大,最大相差49.9%,即便在顆粒初始進(jìn)入流場時,也有1.1%的區(qū)別,這將對顆粒的受力分析產(chǎn)生極大誤差。
(a)t=0 ms (b)t=30 ms
(c)t=60 ms (d)t=90 ms
以Merrill的研究為基礎(chǔ),利用Fluent商用軟件中的滑移網(wǎng)格技術(shù),對典型階段下旋轉(zhuǎn)鋁滴的受力狀態(tài)進(jìn)行分析。鋁滴顆粒初始半徑R=50 μm,計算區(qū)域取得足夠大,以保證顆粒的運動特性不受邊界流動影響:平行于來流方向的長度l=3 mm,垂直于來流方向的寬度d為2 mm。顆粒置于流場中靠近來流方向處,中心與入口邊界相距1.1 mm。采用非正交網(wǎng)格進(jìn)行計算,近壁面處做邊界層處理,顆粒附近區(qū)域進(jìn)行了加密處理,最小網(wǎng)格尺寸為5 μm,最大為170 μm,總網(wǎng)格數(shù)約為35萬,如圖2所示。為便于分析,記來流v的方向為x,顆粒繞z軸順時針旋轉(zhuǎn),軸線與來流方向瞬時針夾角為θ,如圖3所示。
(a)三維模型 (b)中心截面圖
圖3 θ的表示方法Fig.3 Representation of θ
假設(shè)來流是穩(wěn)態(tài)的,其物性在運動過程中不發(fā)生變化,計算的控制方程為
·U=0
(1)
(2)
為簡化模擬條件,不考慮氣固兩相間的傳熱,以空氣代替發(fā)動機(jī)內(nèi)的燃?xì)?,顆粒所受x方向阻力系數(shù)CD由顆粒所受到的阻力FD定義[6]:
(3)
同樣,顆粒所受到的與來流垂直方向的阻力系數(shù),即升力系數(shù)CL,由升力FL定義:
(4)
計算中,顆粒僅有轉(zhuǎn)動,無平動,顆粒雷諾數(shù)Rep小于100,屬層流不可壓縮流動。采用三階精度Quick差分離散格式對方程組離散,利用Simplec算法求解控制方程組。收斂判據(jù)為殘差和小于10-6,迭代時間步長為10-3。
3.1 模型驗證
對于規(guī)則球形顆粒的旋轉(zhuǎn)受力情況,已有學(xué)者進(jìn)行了大量研究,并總結(jié)出相應(yīng)的升力經(jīng)驗公式CL=2Γ[7],Γ為無量綱轉(zhuǎn)速Γ=rω/vr,ω為顆粒旋轉(zhuǎn)角速度。無量綱轉(zhuǎn)速是顆粒角速度與來流速度的比值,可更好地反映顆粒的旋轉(zhuǎn)狀態(tài)。為驗證模型的適用性,針對典型工況進(jìn)行計算,并與經(jīng)驗公式對比。取初始半徑R=50 μm的球形顆粒進(jìn)行驗證計算,無量綱轉(zhuǎn)速Γ分別為0.1、1、2、3。模擬得顆粒的阻力系數(shù)對比結(jié)果如圖4所示。當(dāng)雷諾數(shù)較小時,計算結(jié)果與理論解的最大誤差不超過5%,二者符合良好,證明文中所用模型是正確的。以此模型為計算基礎(chǔ),對其他狀態(tài)下的旋轉(zhuǎn)顆粒進(jìn)行了受力分析。
圖4 升力系數(shù)模擬結(jié)果與理論值的比較Fig.4 Simulation result vs theoretical value
3.2 相對流速變化對顆粒受力的影響
顆粒旋轉(zhuǎn)速度發(fā)生改變時,會影響因旋轉(zhuǎn)而產(chǎn)生的升力,也導(dǎo)致顆粒表面壓力梯度的變化。選取t=90 ms的典型狀態(tài),對比計算轉(zhuǎn)速分別為5、50、100、200 rad/s時顆粒的受力情況。計算結(jié)果如圖5所示。
計算表明,旋轉(zhuǎn)速度對顆粒受力的影響較小,各大轉(zhuǎn)速下,顆粒的升力系數(shù)與阻力系數(shù)相差較小。因此,以鋁滴燃盡狀態(tài)(t=90 ms)為基礎(chǔ),取50 rad/s轉(zhuǎn)速狀態(tài),模擬計算相對流速vr分別為0.5、1.0、1.5、5.0 m/s時顆粒的受力情況,對比分析流場對顆粒受力的影響。4種工況中,顆粒雷諾數(shù)Re分別為0.035、0.07、0.14、0.35,無量綱轉(zhuǎn)速Γ分別為0.007、0.003 5、0.002 3、0.000 7。計算結(jié)果表明,非球顆粒在流場中的升力系數(shù)與無量綱轉(zhuǎn)速之間已經(jīng)不再滿足CL=2Γ的關(guān)系。隨相對流速增加,顆粒的無量綱轉(zhuǎn)速減小,阻力系數(shù)也隨之減小。圖6~圖8詳細(xì)反映了顆粒在一個旋轉(zhuǎn)周期內(nèi)不同方位角處的受力情況。
圖6反映了當(dāng)相對流速由0.5 m/s逐漸變化到5 m/s時,顆粒在y方向所受升力系數(shù)的變化過程。當(dāng)來流速度相對較小時,在一個旋轉(zhuǎn)周期內(nèi),顆粒的升力系數(shù)變化較為明顯。周向不同位置處存在著2組極大值與極小值,它們交替出現(xiàn),每對極值所在位置相差約180°,2組極值的連線交叉近似“十”字形,且阻力系數(shù)在負(fù)向的最大振幅大于正向,這一點在圖6(a)表現(xiàn)得尤為明顯。對比圖6(a) ~(d)可看出,各工況下,阻力系數(shù)極值出現(xiàn)的位置固定,與相對流速無關(guān)。隨來流速度的增加,正負(fù)向極值分別迅速衰減和增加,周向上阻力系數(shù)變得更加均衡,均逐漸逼近于零。由圖6(d)可看出,當(dāng)流速增加到5 m/s時,周向上各位置處的阻力系數(shù)均在0值環(huán)線附近,說明此時流場對顆粒在x向所產(chǎn)生的阻力作用已經(jīng)很弱。
(a) 升力系數(shù)對比
(b) 阻力系數(shù)對比
圖7反映了當(dāng)相對流速由0.5 m/s逐漸變化到5 m/s時,顆粒在x方向所受阻力系數(shù)的變化過程。與x向變化規(guī)律類似,在一個旋轉(zhuǎn)周期內(nèi),周向不同位置處存在著2組交替出現(xiàn)的極值,每對極值所處位置也相差約180°。與x向不同的是顆粒所受到的y向阻力系數(shù)均為正值,且周向變化較為平緩,曲線近似橢圓形。相對而言,水平位置處極值稍大于垂直位置,但二者差別不大。流速對顆粒阻力影響很大,由圖7(a)~(d) 明顯可見,阻力系數(shù)迅速減小,當(dāng)相對流速由0.5 m/s增加至5 m/s時,阻力系數(shù)均值由614.3銳降到64.5。流速越大,對顆粒產(chǎn)生的阻力作用越弱。
圖8(a)、(b)分別直觀反映了顆粒在周向不同方位時所受x方向阻力和y方向升力的波動情況。顆粒在x向的阻力系數(shù)變化接近“正弦”曲線,數(shù)值正負(fù)交替出現(xiàn),說明顆粒在不同位置時所受到x向的阻力方向是變化的,且相對流速越大時,相對應(yīng)的“振幅”越小,曲線越加光滑。當(dāng)相對流速增加到5 m/s時,曲線已接近于數(shù)值為零的直線,比流速為0.5 m/s時的平均值銳減了89%。在顆粒主軸與來流約成90°和270°兩位置處,顆粒所受x向的阻力系數(shù)均為零。顆粒在y向的阻力系數(shù)曲線波動較小,相對流速愈大,曲線愈加平滑,當(dāng)流速增加到5 m/s時,曲線已經(jīng)接近于直線。與x向阻力變化類似,y向阻力均值也銳減了89%。
(a) vr=0.5 m/s
(b) vr=1 m/s
(c) vr=2 m/s
(d) vr=5 m/s
(a) vr=0.5 m/s
(b) vr=1 m/s
(c) vr=2 m/s
(d) vr=5 m/s
3.3 顆粒升力系數(shù)和阻力系數(shù)與顆粒雷諾數(shù)的關(guān)系
計算表明,顆粒的阻力系數(shù)與來流相對速度成正比,當(dāng)物性條件確定時,顆粒的雷諾數(shù)Re只與相對流速有關(guān)。當(dāng)旋轉(zhuǎn)速度確定時,顆粒的無量綱轉(zhuǎn)速Γ也只與相對來流速度有關(guān),即與雷諾數(shù)Re相關(guān),由此可得到阻力系數(shù)與雷諾數(shù)Re的函數(shù)關(guān)系,即CD=f(Re)。為探索二者的關(guān)系,首先統(tǒng)計顆x向無量綱相對阻力系數(shù)與相對雷諾數(shù)的比值,如表1所示。
(a) 升力系數(shù)對比
(b) 阻力系數(shù)對比
表1 顆粒相對阻力系數(shù)與相對雷諾數(shù)的比值統(tǒng)計Table 1 Statistics of the numberless drag coefficient and Reynolds number
由表1可知,當(dāng)旋轉(zhuǎn)角速度恒定時,顆粒的相對阻力系數(shù)與顆粒雷諾數(shù)之比值滿足:
(5)
式中CD1和CD2與Re1和Re2代表2種不同工況下的阻力系數(shù)和雷諾數(shù);α為一常數(shù)。
由此可推斷,顆粒阻力系數(shù)與雷諾數(shù)滿足以下關(guān)系:
CD=βRe-1
(6)
式中β為常數(shù)。
引入球形度(Sphericity)的概念[8-9],ψ=(dv/ds),可得到顆粒的球形度為0.912。Haider和Chhabra[10]曾提出球形度與阻力系數(shù)間存在著一定的關(guān)系:
(7)
其中
A=exp(2.328 8-6.458 1ψ+2.448 6ψ2)
B=0.096 4+0.556 5ψ
C=exp(4.905-13.894 4ψ+18.422 2ψ2-10.259 9ψ3)
D=exp(1.468 1+12.258 4ψ-20.732 2ψ2+15.885 5ψ3)
當(dāng)顆粒的球形度ψ=0.912時,可求得顆粒的阻力系數(shù)與雷諾數(shù)間關(guān)系為
CD=24.48Re-1
(8)
此式與式(6)完全符合,此時β取為24.48,這也再次印證了本文計算所得結(jié)果可信。
統(tǒng)計顆粒y向升力系數(shù)計算結(jié)果,得到表2。與表1的統(tǒng)計結(jié)果相比,表2中旋轉(zhuǎn)周向不同位置處,顆粒的無量綱升力系數(shù)與雷諾數(shù)比值存在著差別,并不是常數(shù)關(guān)系。
表2 顆粒相對升力系數(shù)與相對雷諾數(shù)的比值統(tǒng)計Table2 Statistics of the numberless lift coefficient and Reynolds number
式(6)應(yīng)修正為
(9)
即比例系數(shù)α是方位角θ的函數(shù),α=f(θ)。統(tǒng)計表明,當(dāng)顆粒雷諾數(shù)Re=0.14~0.35時,α值為0.8~1.17。此時,顆粒升力系數(shù)與雷諾數(shù)Re之間的關(guān)系也應(yīng)修正為
CL=C(θ)Re-1
(10)
(1)固體火箭發(fā)動機(jī)燃燒室環(huán)境中,非球形鋁滴的旋轉(zhuǎn)受力與經(jīng)典球形顆粒的理論解有很大區(qū)別,顆粒所受x向阻力系數(shù)與y向升力系數(shù)變化規(guī)律差別較大,阻力系數(shù)明顯大于升力系數(shù)。小雷諾數(shù)條件下,一個旋轉(zhuǎn)周期內(nèi),周向不同位置處的升力系數(shù)方向發(fā)生變化,存在著交替出現(xiàn)的兩組極大值與極小值。其中,負(fù)向振幅大于正向,阻力系數(shù)極值出現(xiàn)的位置固定,與相對流速無關(guān);阻力系數(shù)均為正,但周向變化較為平緩,各方位角處數(shù)值相差不大。隨來流速度增加,升力系數(shù)變得均衡,系數(shù)值逐漸逼近于零;阻力系數(shù)也急速降低。當(dāng)相對流速由0.5 m/s提高到5 m/s時,顆粒所受阻力銳降89%。流速越大,對顆粒產(chǎn)生的阻力作用越弱。
(2)顆粒的相對升力系數(shù)與顆粒雷諾數(shù)之比值滿足:CL1/CL2=α(θ)Re2/Re1,比例系數(shù)α是方位角θ的函數(shù)。當(dāng)顆粒雷諾數(shù)Re=0.14~0.35時,升力系數(shù)中α值為0.8~1.17,阻力系數(shù)中α=1。顆粒阻力系數(shù)與雷諾數(shù)滿足關(guān)系:CD=β(θ)Re-1。其中,x向阻力系數(shù)中,β為常數(shù);y向升力系數(shù)中,β是方位角θ的函數(shù)。
[1] 方丁酉.兩相流體力學(xué)[M]. 長沙:國防科技大學(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] Robert Geisler. A global view of the use of aluminum fuel in solid rocket motors[R].AIAA 2002-3748.
[5] Merrill K King. Aluminum combustion in a solid rocket motor environment[J]. Proceedings of the Combustion Institute, 2009, 32(2):2107-2114.
[6] 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.
[7] Rubinow S I, Keller J B. The transerverse force on a spinning sphere moving in a viscous fluid[J]. Fluid Mechanics, 1961, 11: 447-459.
[8] 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), 2000, 5(1):31-38.
[9] 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.
[10] Chhabra R P, Agarwal L, Sinha N Kl. Drag on non spherical particles: an evaluation of available methods[J].Powder Technology,1999,101:288-295.
(編輯:崔賢彬)
Influence of relative velocity on the force of non-spherical rotating aluminum droplet
LIU Cong-lin1, GAO Ye1, HE Zheng1, KANG Wen-wu2
(1. Aerospace and Civil Engineering, Harbin Engineering University, Harbin 150001, China; 2. School of Energy and Power Engineering, Beihang University, Beijing 100191,China)
The rotation characteristics of non-sphere aluminum in different kinds of fluid velocity were simulated in order to analyze the effect of flow rate on the forces acting on particles. Results show that, the relationship between lift coefficient and dimensionless rotational speed doesn't obey the typical law. When Reynolds number is less than 1, two groups of maximum and minimum values appear alternately on certain circumferential location for drag coefficient. Lift coefficient is always positive and smooth. With the increasing of the inflow velocity, drag coefficient becomes even, and the value is approaching to zero, while the lift coefficient falls quickly. If the relative velocity increase from 0.5 m/s to 5 m/s, drag coefficient of the droplet decreases by 89% sharply. Drag coefficient and Reynolds number obey a new law.
non-spherical aluminum droplet; relative velocity; drag coefficient; lift coefficient
2014-05-06;
:2014-06-18。
國家自然科學(xué)基金(11372079);中央高?;究蒲袠I(yè)務(wù)費專項資金(HEUCF 130203)。
劉叢林(1981—),女,助理研究員,研究方向為多相流內(nèi)流場燃燒。E-mail:liuconglin2006@126.com
V435
A
1006-2793(2015)02-0214-06
10.7673/j.issn.1006-2793.2015.02.012