相對流速對非球形鋁滴旋轉受力的影響①

劉叢林，郜 冶，賀 征，康文武

(1.哈爾濱工程大學 航天與建筑工程學院，哈爾濱 150001;2.北京航空航天大學 能源與動力工程學院，北京 100191)

?

相對流速對非球形鋁滴旋轉受力的影響①

劉叢林1，郜 冶1，賀 征1，康文武2

(1.哈爾濱工程大學 航天與建筑工程學院，哈爾濱 150001;2.北京航空航天大學 能源與動力工程學院，北京 100191)

針對固體火箭發動機中非球形鋁滴的旋轉受力問題進行數值模擬，詳細分析流速變化對顆粒的受力影響。計算結果表明，非球顆粒在流場中的升力系數與無量綱轉速之間已經不再滿足經典的理論關系式。雷諾數Re<1時，周向升力系數存在著交替出現的2組極大值與極小值；升力系數均為正值，且周向變化較為平緩。隨來流速度增加，阻力系數變得均衡，系數值向0值逼近；升力系數也急速降低。當相對流速由0.5 m/s提高到5 m/s時，顆粒所受阻力銳降89%。顆粒的阻力系數與雷諾數之間滿足新的關系式。

非球形鋁滴；相對流速；阻力系數；升力系數

0 引言

在含顆粒的多相流計算中，一般視顆粒為離散相，顆粒在流場中的運動軌跡及空間分布取決于作用于其上的力和力矩，對流場中的離散相進行受力分析是獲得準確計算結果的基礎。尤其固體火箭發動機中，金屬鋁滴在流場中相變燃燒時，會因部分氧化鋁凝相回落到顆粒表面而轉變為不規則的非球體[1-2]，導致動力學特性發生一系列改變，從而影響整個多相流場的流動特性。但經典的顆粒受力理論是以規則球形為基礎進行推導的，實際應用中將不可避免地產生偏差[3]。另外，實際流場中，顆粒的運動狀態往往是轉動伴隨平動同時存在，但為簡化計算，數值模擬中,前者往往被忽略，理論上講，顆粒在發動機中所受阻力與相對流速度呈平方反比關系。所以，流速變化對顆粒的受力影響很大。固體火箭發動機中，在推進劑剛開始燃燒時，燃氣速度較小，但幾十毫秒后，燃氣夾帶著顆粒向噴管加速流動，由于慣性作用，顆粒相速度滯后，與主燃氣間存在著速度差，再加上自身旋轉作用，受力狀態較為復雜。在固體火箭發動機內流場中，金屬顆粒的運動特性對整個工作過程起至關重要的作用，綜合考慮實際環境中可能存在的顆粒形態及運動狀態，針對旋轉非球形單顆粒的受力問題進行詳細數值分析，對進一步分析顆粒間凝聚、碰撞等運動特性以及顆粒與流體間相互作用有一定的參考作用，對進一步認清多相流的本質特征、完善多相流模型有重要意義。

1 物理模型

金屬鋁極易氧化，鋁顆粒存放期間，其表面就包裹著氧化物，初始進入發動機時，顆粒形狀已經發生變化，不再是規則的球形[4]，受力狀態與運動特性亦發生相應改變。Merrill[5]研究了鋁顆粒在固體火箭發動機環境中的相變燃燒過程，對其結構變化過程做了仔細地推導和計算。結果表明，顆粒將在90 ms內完成燃燒過程，其間主要以Al和Al2O3兩球形相切的方式出現，隨燃燒的發展，外形不斷變化，如圖1所示。數值模擬中，通常以等體積球體代替非球體進行受力分析。對Merrill教授的研究結果進行統計可知，隨鋁滴在燃燒室中狀態的改變，等體積顆粒在來流方向投影的橫截面面積偏離真實顆粒的水平較大，最大相差49.9%，即便在顆粒初始進入流場時，也有1.1%的區別，這將對顆粒的受力分析產生極大誤差。

(a)t=0 ms (b)t=30 ms

(c)t=60 ms (d)t=90 ms

以Merrill的研究為基礎，利用Fluent商用軟件中的滑移網格技術，對典型階段下旋轉鋁滴的受力狀態進行分析。鋁滴顆粒初始半徑R=50 μm，計算區域取得足夠大，以保證顆粒的運動特性不受邊界流動影響：平行于來流方向的長度l=3 mm，垂直于來流方向的寬度d為2 mm。顆粒置于流場中靠近來流方向處，中心與入口邊界相距1.1 mm。采用非正交網格進行計算，近壁面處做邊界層處理，顆粒附近區域進行了加密處理，最小網格尺寸為5 μm，最大為170 μm，總網格數約為35萬，如圖2所示。為便于分析，記來流v的方向為x，顆粒繞z軸順時針旋轉，軸線與來流方向瞬時針夾角為θ，如圖3所示。

(a)三維模型 (b)中心截面圖

圖3 θ的表示方法Fig.3 Representation of θ

2 控制方程

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

·U=0

(1)

(2)

為簡化模擬條件，不考慮氣固兩相間的傳熱，以空氣代替發動機內的燃氣，顆粒所受x方向阻力系數CD由顆粒所受到的阻力FD定義[6]：

(3)

同樣，顆粒所受到的與來流垂直方向的阻力系數，即升力系數CL，由升力FL定義：

(4)

計算中，顆粒僅有轉動，無平動，顆粒雷諾數Rep小于100，屬層流不可壓縮流動。采用三階精度Quick差分離散格式對方程組離散，利用Simplec算法求解控制方程組。收斂判據為殘差和小于10-6，迭代時間步長為10-3。

3 結果與分析

3.1 模型驗證

對于規則球形顆粒的旋轉受力情況，已有學者進行了大量研究，并總結出相應的升力經驗公式CL=2Γ[7]，Γ為無量綱轉速Γ=rω/vr，ω為顆粒旋轉角速度。無量綱轉速是顆粒角速度與來流速度的比值，可更好地反映顆粒的旋轉狀態。為驗證模型的適用性，針對典型工況進行計算，并與經驗公式對比。取初始半徑R=50 μm的球形顆粒進行驗證計算，無量綱轉速Γ分別為0.1、1、2、3。模擬得顆粒的阻力系數對比結果如圖4所示。當雷諾數較小時，計算結果與理論解的最大誤差不超過5%，二者符合良好，證明文中所用模型是正確的。以此模型為計算基礎，對其他狀態下的旋轉顆粒進行了受力分析。

圖4 升力系數模擬結果與理論值的比較Fig.4 Simulation result vs theoretical value

3.2 相對流速變化對顆粒受力的影響

顆粒旋轉速度發生改變時，會影響因旋轉而產生的升力，也導致顆粒表面壓力梯度的變化。選取t=90 ms的典型狀態，對比計算轉速分別為5、50、100、200 rad/s時顆粒的受力情況。計算結果如圖5所示。

計算表明，旋轉速度對顆粒受力的影響較小，各大轉速下，顆粒的升力系數與阻力系數相差較小。因此，以鋁滴燃盡狀態(t=90 ms)為基礎，取50 rad/s轉速狀態，模擬計算相對流速vr分別為0.5、1.0、1.5、5.0 m/s時顆粒的受力情況，對比分析流場對顆粒受力的影響。4種工況中，顆粒雷諾數Re分別為0.035、0.07、0.14、0.35，無量綱轉速Γ分別為0.007、0.003 5、0.002 3、0.000 7。計算結果表明，非球顆粒在流場中的升力系數與無量綱轉速之間已經不再滿足CL=2Γ的關系。隨相對流速增加，顆粒的無量綱轉速減小，阻力系數也隨之減小。圖6～圖8詳細反映了顆粒在一個旋轉周期內不同方位角處的受力情況。

圖6反映了當相對流速由0.5 m/s逐漸變化到5 m/s時，顆粒在y方向所受升力系數的變化過程。當來流速度相對較小時，在一個旋轉周期內，顆粒的升力系數變化較為明顯。周向不同位置處存在著2組極大值與極小值，它們交替出現，每對極值所在位置相差約180°，2組極值的連線交叉近似“十”字形，且阻力系數在負向的最大振幅大于正向，這一點在圖6(a)表現得尤為明顯。對比圖6(a) ～(d)可看出，各工況下，阻力系數極值出現的位置固定，與相對流速無關。隨來流速度的增加，正負向極值分別迅速衰減和增加，周向上阻力系數變得更加均衡，均逐漸逼近于零。由圖6(d)可看出，當流速增加到5 m/s時，周向上各位置處的阻力系數均在0值環線附近，說明此時流場對顆粒在x向所產生的阻力作用已經很弱。

(a) 升力系數對比

(b) 阻力系數對比

圖7反映了當相對流速由0.5 m/s逐漸變化到5 m/s時，顆粒在x方向所受阻力系數的變化過程。與x向變化規律類似，在一個旋轉周期內，周向不同位置處存在著2組交替出現的極值，每對極值所處位置也相差約180°。與x向不同的是顆粒所受到的y向阻力系數均為正值，且周向變化較為平緩，曲線近似橢圓形。相對而言，水平位置處極值稍大于垂直位置，但二者差別不大。流速對顆粒阻力影響很大，由圖7(a)～(d) 明顯可見，阻力系數迅速減小，當相對流速由0.5 m/s增加至5 m/s時，阻力系數均值由614.3銳降到64.5。流速越大，對顆粒產生的阻力作用越弱。

圖8(a)、(b)分別直觀反映了顆粒在周向不同方位時所受x方向阻力和y方向升力的波動情況。顆粒在x向的阻力系數變化接近“正弦”曲線，數值正負交替出現，說明顆粒在不同位置時所受到x向的阻力方向是變化的，且相對流速越大時，相對應的“振幅”越小，曲線越加光滑。當相對流速增加到5 m/s時，曲線已接近于數值為零的直線，比流速為0.5 m/s時的平均值銳減了89%。在顆粒主軸與來流約成90°和270°兩位置處，顆粒所受x向的阻力系數均為零。顆粒在y向的阻力系數曲線波動較小，相對流速愈大，曲線愈加平滑，當流速增加到5 m/s時，曲線已經接近于直線。與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 顆粒升力系數和阻力系數與顆粒雷諾數的關系

計算表明，顆粒的阻力系數與來流相對速度成正比，當物性條件確定時，顆粒的雷諾數Re只與相對流速有關。當旋轉速度確定時，顆粒的無量綱轉速Γ也只與相對來流速度有關，即與雷諾數Re相關，由此可得到阻力系數與雷諾數Re的函數關系，即CD=f(Re)。為探索二者的關系，首先統計顆x向無量綱相對阻力系數與相對雷諾數的比值，如表1所示。

(a) 升力系數對比

(b) 阻力系數對比

表1 顆粒相對阻力系數與相對雷諾數的比值統計Table 1 Statistics of the numberless drag coefficient and Reynolds number

由表1可知，當旋轉角速度恒定時，顆粒的相對阻力系數與顆粒雷諾數之比值滿足：

(5)

式中CD1和CD2與Re1和Re2代表2種不同工況下的阻力系數和雷諾數；α為一常數。

由此可推斷，顆粒阻力系數與雷諾數滿足以下關系：

CD=βRe-1

(6)

式中β為常數。

引入球形度(Sphericity)的概念[8-9]，ψ=(dv/ds)，可得到顆粒的球形度為0.912。Haider和Chhabra[10]曾提出球形度與阻力系數間存在著一定的關系：

(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)

當顆粒的球形度ψ=0.912時，可求得顆粒的阻力系數與雷諾數間關系為

CD=24.48Re-1

(8)

此式與式(6)完全符合，此時β取為24.48，這也再次印證了本文計算所得結果可信。

統計顆粒y向升力系數計算結果，得到表2。與表1的統計結果相比，表2中旋轉周向不同位置處，顆粒的無量綱升力系數與雷諾數比值存在著差別，并不是常數關系。

表2 顆粒相對升力系數與相對雷諾數的比值統計Table2 Statistics of the numberless lift coefficient and Reynolds number

式(6)應修正為

(9)

即比例系數α是方位角θ的函數，α=f(θ)。統計表明，當顆粒雷諾數Re=0.14～0.35時，α值為0.8～1.17。此時，顆粒升力系數與雷諾數Re之間的關系也應修正為

CL=C(θ)Re-1

(10)

4 結論

(1)固體火箭發動機燃燒室環境中，非球形鋁滴的旋轉受力與經典球形顆粒的理論解有很大區別，顆粒所受x向阻力系數與y向升力系數變化規律差別較大，阻力系數明顯大于升力系數。小雷諾數條件下，一個旋轉周期內，周向不同位置處的升力系數方向發生變化，存在著交替出現的兩組極大值與極小值。其中，負向振幅大于正向，阻力系數極值出現的位置固定，與相對流速無關；阻力系數均為正，但周向變化較為平緩，各方位角處數值相差不大。隨來流速度增加，升力系數變得均衡，系數值逐漸逼近于零；阻力系數也急速降低。當相對流速由0.5 m/s提高到5 m/s時，顆粒所受阻力銳降89%。流速越大，對顆粒產生的阻力作用越弱。

(2)顆粒的相對升力系數與顆粒雷諾數之比值滿足：CL1/CL2=α(θ)Re2/Re1，比例系數α是方位角θ的函數。當顆粒雷諾數Re=0.14～0.35時，升力系數中α值為0.8～1.17，阻力系數中α=1。顆粒阻力系數與雷諾數滿足關系：CD=β(θ)Re-1。其中，x向阻力系數中，β為常數；y向升力系數中，β是方位角θ的函數。

[1] 方丁酉．兩相流體力學[M]. 長沙：國防科技大學出版社，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] 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。

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

劉叢林(1981—)，女，助理研究員，研究方向為多相流內流場燃燒。E-mail:liuconglin2006@126.com

V435

A

1006-2793(2015)02-0214-06

10.7673/j.issn.1006-2793.2015.02.012