流化床內(nèi)顆粒旋轉特性的數(shù)值模擬

李 斌，宋小龍

（華北電力大學 能源動力與機械工程學院，保定071003）

流化床作為高效、低污染清潔燃燒技術在電力行業(yè)中得到了較為廣泛的應用，但是其床內(nèi)復雜的氣固兩相流流動行為還難以被人們充分了解.目前，關于流化床內(nèi)顆粒速度、濃度分布的研究已相對成熟.由于受到不平衡力及顆粒間碰撞的作用，顆粒在床內(nèi)平動的同時還會產(chǎn)生旋轉，顆粒的旋轉不僅會對自身的運動產(chǎn)生影響，還會對其所在的流場產(chǎn)生影響［1］，因此顆粒旋轉特性的研究作為顆粒速度和濃度等研究的深入對于全面認識氣固兩相流的流動特性具有重要意義.

隨著研究的深入，流化床內(nèi)顆粒的旋轉逐漸引起人們的重視.目前，對顆粒旋轉特性的研究主要通過試驗研究和數(shù)值模擬方法，試驗研究的關鍵在于如何準確測量顆粒的轉速，有學者［2-4］對顆粒的旋轉特性進行了試驗研究，但是還需要進一步研究更精確地測量顆粒轉速的方法.此外，也有學者采用硬球模型、直接蒙特卡羅法或雙流體模型對流化床內(nèi)顆粒的旋轉特性進行了數(shù)值模擬［5-8］，但是對流化床內(nèi)顆粒旋轉特性采用軟球模型進行數(shù)值模擬的研究并未深入.

數(shù)值模擬的關鍵是建立正確的物理數(shù)學模型，流化床內(nèi)氣固兩相流數(shù)值模擬的模型主要包括基于歐拉方法的雙流體模型［9-10］和基于拉格朗日方法的顆粒軌道模型［11-12］.筆者將基于歐拉-拉格朗日范疇的離散單元法和計算流體力學結合起來，在自行開發(fā)程序的基礎上對流化床內(nèi)顆粒的旋轉特性進行了數(shù)值模擬，充分發(fā)揮離散單元法可以獲得顆粒非常豐富的微觀信息的優(yōu)勢，模擬得到床內(nèi)氣體的速度場、顆粒的速度場和平均體積分數(shù)以及顆粒平均轉速的分布，分析了顆粒旋轉對顆粒自身運動、流場、床層膨脹和床層空隙率的影響以及影響顆粒平均轉速的因素，為全面認識和深入研究流化床內(nèi)的氣固兩相流提供一定基礎.

1 顆粒旋轉機理

流化床中的顆粒在床內(nèi)運動過程中主要受到以下力的作用：（1）顆粒與顆粒碰撞產(chǎn)生的碰撞力；（2）顆粒與壁面碰撞產(chǎn)生的力；（3）氣體對顆粒產(chǎn)生的曳力；（4）顆粒自身重力.床內(nèi)顆粒遇到以下情形時會發(fā)生旋轉［5］：當2個顆粒發(fā)生偏心碰撞時，在接觸點產(chǎn)生的碰撞力可分解成一個垂直于接觸點與顆粒圓心連線方向的分力，該力產(chǎn)生轉矩使顆粒發(fā)生旋轉；當顆粒與壁面發(fā)生非垂直碰撞時，由顆粒與壁面間的摩擦力產(chǎn)生對顆粒圓心的轉矩，使顆粒旋轉；當顆粒在具有一定速度梯度的流場中運動時，在不平衡力的作用下也會發(fā)生旋轉；此外，當一個旋轉顆粒與一個非旋轉顆粒發(fā)生碰撞時，在顆粒間摩擦力的作用下，也會使非旋轉顆粒產(chǎn)生旋轉.

根據(jù)文獻［13］，由顆粒碰撞引起的顆粒轉速比由氣體速度梯度引起的顆粒轉速大得多，顆粒的碰撞成為引起顆粒旋轉的主要因素.筆者在模型處理時只考慮由顆粒碰撞引起的顆粒旋轉.

2 數(shù)學模型

將計算流體力學和基于歐拉-拉格朗日范疇的離散單元法相結合，對流化床內(nèi)顆粒的旋轉特性進行數(shù)值模擬.在歐拉坐標系下考察氣相的運動，在拉格朗日坐標系下計算顆粒的運動，顆粒間的碰撞采用軟球模型，根據(jù)牛頓第二定律建立每一個顆粒的運動方程，然后對每一個顆粒的運動方程進行求解，氣固兩相耦合采用牛頓第三定律.

2.1 顆粒間的碰撞力

軟球模型中顆粒間的碰撞力包括彈性力和黏性阻尼力.根據(jù)物理學定律，顆粒間的碰撞力可分解為法向力fn，ij和切向力ft，ij.法向力包括法向彈性力fcn，ij和法向黏性阻尼力fdn，ij，切向 力ft，ij包括切向彈性力fct，ij和切向黏性阻尼力fdt，ij.

顆粒所受的法向力為

式中：kn為顆粒法向方向的彈性系數(shù)；Δδn為顆粒i和顆粒j之間的法向相對位移.

式中：ur為顆粒i相對顆粒j的速度矢量，ur＝uiuj；ni為顆粒法向方向的單位矢量，方向由顆粒i的圓心指向接觸點.

法向黏性阻尼力為

式中：ηn為顆粒法向方向的阻尼系數(shù)；un為顆粒i和顆粒j在接觸點處的法向相對速度.

顆粒所受的切向彈性力為

式中：kt為顆粒切向方向的彈性系數(shù)；Δδt為顆粒i和顆粒j之間的切向相對位移.

式中：ti為顆粒i切向方向的單位矢量；ωi為顆粒的旋轉角速度，rad／s，規(guī)定逆時針方向為正；Ri為顆粒質(zhì)心到接觸點的矢量.

切向黏性阻尼力為

式中：ηt為顆粒切向方向的阻尼系數(shù)；ut為顆粒i和顆粒j在接觸點的切向相對速度.

式中：μ為顆?；瑒幽Σ烈驍?shù).

當顆粒i與壁面發(fā)生碰撞時，將壁面看成是速度和角速度均為零的顆粒.

2.2 顆粒運動方程

顆粒在床內(nèi)的運動可以分解為平動和轉動.

平動運動方程為

式中：Fy，i為流體對顆粒的曳力；Fp為顆粒間的碰撞力.

轉動運動方程為

式中：Ii為顆粒i的轉動慣量，對于球形顆粒Ii＝2／5miR2i；Tij為顆粒i與顆粒j發(fā)生碰撞時對顆粒i產(chǎn)生的轉矩；k為與顆粒i發(fā)生碰撞的顆粒個數(shù).

當一個顆粒同時與周圍的幾個顆粒碰撞時，可以通過矢量疊加計算顆粒所受的合力及合力矩.

經(jīng)過一個時間步長Δt后，顆粒i的速度、位移和角速度可由下式確定：

式中：ui0、si0和ωi0分別為上一時刻顆粒i的速度、位移和角速度.

2.3 氣相控制方程

氣相模型采用考慮氣固兩相耦合的Navier-Stocks方程，氣相湍流模型采用k-ε兩方程模型，并采用Simpler算法進行求解，詳細求解模型參見文獻［14］～文獻［15］.

2.4 模擬對象

所模擬的對象為150mm（x）×4mm（y）×900 mm（z）的矩形截面準三維流化床，床深（y方向）為顆粒直徑，床層底部中心位置設置一個空氣進口，空氣進口寬度為10mm.計算區(qū)域的網(wǎng)格劃分為17×1×92，水平方向（x方向）和垂直方向（z方向）的網(wǎng)格長度均為10mm.床料采用2 400個球形顆粒，其中顆粒的恢復系數(shù)e、彈性系數(shù)k和摩擦因數(shù)分別為0.9、200 N／m 和0.3，黏 性 系 數(shù) 采 用2×（m為等效質(zhì)量）進行計算.空氣密度為1.205kg／m3，動力黏度為1.8×10-5kg／（m·s）.

2.5 模擬過程

首先確定床內(nèi)顆粒的初始位置，方法如下：在床內(nèi)隨機生成2 400個顆粒位置，在只考慮重力和顆粒間碰撞力的作用下使其自由落體，顆粒最終的靜止位置作為數(shù)值模擬過程中顆粒的初始位置.為了對比考慮顆粒旋轉和未考慮顆粒旋轉的情況，在模擬過程中將顆粒的旋轉角速度隨時清零，即可計算未考慮顆粒旋轉的情況［5］.

3 模擬結果

筆者模擬了空氣進口速度為29m／s、顆粒直徑為4mm 及密度為2 700kg／m3時床內(nèi)顆粒的旋轉特性，并對不同密度（直徑均為4 mm，顆粒密度分別為3 900kg／m3、2 700kg／m3和1 500kg／m3）以及不同直徑（顆粒密度均為2 700kg／m3，顆粒直徑分別為3mm、4mm 和5mm）顆粒的旋轉特性進行了數(shù)值模擬，得到不同密度及不同直徑顆粒的平均轉速分布.

3.1 流化過程圖

首先模擬了空氣進口速度為29m／s時考慮顆粒旋轉和未考慮顆粒旋轉的床內(nèi)顆粒流化過程，得到加入流化氣體后床內(nèi)顆粒的流化過程圖（圖1和圖2）.對比圖1和圖2可以看出，考慮顆粒旋轉后床層的膨脹高度比未考慮顆粒旋轉時有所提高，且在床內(nèi)更容易形成氣泡.

圖1 未考慮顆粒旋轉的床內(nèi)顆粒流化過程圖Fig.1 Diagram of fluidization process without consideration of particle rotation

圖2 考慮顆粒旋轉的床內(nèi)顆粒流化過程圖Fig.2 Diagram of fluidization process with consideration of particle rotation

3.2 顆粒平均速度分布

圖3和圖4分別給出了空氣進口速度為29m／s時顆粒垂直方向和水平方向平均速度的分布，其中h為床層高度（以下簡稱床高）.從圖3和圖4可以看出，考慮顆粒旋轉和不考慮顆粒旋轉時顆粒垂直方向和水平方向的平均速度分布趨勢一致.從圖3還可知，在垂直方向上床層中心位置處的顆粒平均速度較大且為正值，顆粒向上運動；靠近壁面處的顆粒平均速度為負值，顆粒向下運動，從而形成流化床內(nèi)顆粒環(huán)核流動，并且隨床高的增加，顆粒垂直方向的平均速度減小.此外通過對比考慮顆粒旋轉和不考慮顆粒旋轉情況下顆粒垂直方向的平均速度可知，考慮顆粒旋轉后床層中心位置處的顆粒平均速度變大，但是靠近壁面兩側的顆粒平均速度變小.從圖4還可以看出，在水平方向上，靠近床層底部區(qū)域噴口中心位置左側顆粒的平均速度為正值且向右側運動，右側顆粒的平均速度為負值且向左側運動，兩側的顆粒向噴口中心位置聚集并進入主氣流的噴射區(qū)，由于受到主氣流的卷吸作用，該區(qū)域顆粒水平方向的平均速度較大，其他區(qū)域的顆粒分別向兩側擴散，床層表面的顆粒由于受到中間顆粒向上的排擠以及受周圍顆粒的影響較小，其顆粒平均速度也較大.考慮顆粒旋轉后顆粒水平方向平均速度較不考慮顆粒旋轉時有所增大.

圖3 顆粒垂直方向平均速度分布Fig.3 Vertical velocity distribution of particles

圖4 顆粒水平方向平均速度分布Fig.4 Horizontal velocity distribution of particles

3.3 氣體平均速度分布

圖5和圖6分別給出了空氣進口速度為29m／s時氣體垂直方向和水平方向平均速度的分布.從圖5可知，垂直方向上中心位置處的氣體平均速度較大且方向為正，兩側的氣體平均速度較小且方向為負，并且氣體平均速度隨床高的增加而減小.氣體水平方向平均速度以噴口為中心基本呈對稱分布，噴口中心兩側的平均速度方向相反，在床層底部區(qū)域由于受主氣流的卷吸作用，兩側的氣體向噴口中心聚集，其他區(qū)域的氣體向兩側擴散，并且隨床高的增加氣體水平方向平均速度減小.對比考慮顆粒旋轉和不考慮顆粒旋轉的情況可知，顆粒的旋轉對氣體水平方向平均速度的影響略大于對垂直方向平均速度的影響，但其規(guī)律并不明顯.

圖5 氣體垂直方向平均速度分布Fig.5 Vertical velocity distribution of gas

圖6 氣體水平方向平均速度分布Fig.6 Horizontal velocity distribution of gas

3.4 顆粒平均體積分數(shù)分布

圖7給出了考慮顆粒旋轉和未考慮顆粒旋轉時床內(nèi)顆粒平均體積分數(shù)的分布.從圖7可以看出，兩種情況下顆粒平均體積分數(shù)的分布趨勢相似，床內(nèi)底部區(qū)域顆粒的平均體積分數(shù)較大且為濃相，上部區(qū)域顆粒的平均體積分數(shù)較小且為稀相，靠近壁面兩側的顆粒平均體積分數(shù)大于中心區(qū)域，呈現(xiàn)中心稀兩側密的特點.

圖7 顆粒平均體積分數(shù)的分布Fig.7 Distribution of average volume fraction of particles

3.5 床層平均空隙率隨床高的變化趨勢

圖8給出了床層平均空隙率隨床高的變化趨勢.從圖8可以看出，在床層底部考慮顆粒旋轉時床層平均空隙率較未考慮顆粒旋轉時有所增大，這與文獻［5］采用硬球模型得到的結果基本一致，分析認為這是由于考慮顆粒旋轉后床層底部區(qū)域顆粒平均速度增大的幅度較大，顆粒運動范圍擴大，導致床層的平均空隙率增大.

圖8 床層平均空隙率隨床高的變化趨勢Fig.8 Distribution of average bed void fraction along bed height

3.6 顆粒碰撞次數(shù)的分布

圖9給出了0～2s內(nèi)床內(nèi)不同網(wǎng)格區(qū)域內(nèi)顆粒沿床寬發(fā)生碰撞次數(shù)總和的分布趨勢.從圖9可以看出，顆粒的碰撞次數(shù)隨著床高的增加而減小，兩側的碰撞頻率略高于中心區(qū)域，這與顆粒平均體積分數(shù)的分布一致，這是由于顆粒平均體積分數(shù)較大區(qū)域的顆粒碰撞概率較高，顆粒碰撞次數(shù)增加.

圖9 考慮顆粒旋轉時顆粒的碰撞次數(shù)Fig.9 Distribution of collision frequency with consideration of particle rotation

3.7 顆粒平均轉速的分布

圖10給出了不同床高處顆粒平均轉速沿床寬的分布情況.從圖10可以看出，床內(nèi)噴口中心左側顆粒的平均轉速為正值，顆粒逆時針旋轉；右側顆粒的平均轉速為負值，顆粒順時針旋轉，顆粒的平均轉速分布與顆粒水平方向的平均速度分布趨勢基本一致，水平方向平均速度較大顆粒的平均轉速也較大，這是由于水平方向平均速度大的顆粒與其他顆粒和壁面的碰撞概率增加，顆粒經(jīng)歷的流場相對復雜，從而增加了顆粒增大平均轉速的機會［3］.隨著床高的增加，顆粒平均轉速減小，這與顆粒碰撞次數(shù)的分布相吻合，但靠近壁面處顆粒的平均轉速小于兩側中心區(qū)域，分析認為這是由于顆粒的平均轉速不僅與碰撞次數(shù)有關，還與顆粒間的碰撞強度有關，雖然靠近壁面處的顆粒體積分數(shù)較大，碰撞次數(shù)較多，但是靠近壁面處顆粒的平均速度較小，顆粒碰撞強度較弱，因此顆粒的平均轉速相對較小.從圖10還可以看出，顆粒平均轉速與顆粒垂直方向平均速度的關系不大.

圖10 顆粒平均轉速的分布Fig.10 Average rotating speed distribution of particles

3.8 不同密度顆粒的平均轉速隨床寬的變化趨勢

圖11給出了空氣進口速度為29 m／s、顆粒直徑均為4mm 時不同密度顆粒的平均轉速隨床寬的變化.從圖11可以看出，密度較大的顆粒的平均轉速小于密度較小的顆粒，這是由于密度較大的顆粒轉動慣量較大，在受到相同轉矩時，顆粒的角加速度較小，因而顆粒角速度和平均轉速較小.

圖11 不同密度顆粒的平均轉速隨床寬的變化趨勢Fig.11 Distribution of average rotating speed along bed width for particles of different densities

3.9 不同直徑顆粒的平均轉速隨床寬的變化趨勢

圖12給出了空氣進口速度為29 m／s、顆粒密度均為2 700kg／m3時不同直徑顆粒的平均轉速隨床寬的變化.從圖12可知，直徑小的顆粒平均轉速大于直徑大的顆粒，這與文獻［3］的結果一致，這是由于在受到相同轉矩時，直徑較小的顆粒轉動慣量較小，因而可以獲得更大的平均轉速.

圖12 不同直徑顆粒的平均轉速隨床寬的變化Fig.12 Distribution of average rotating speed along bed width for particles of different sizes

4 結 論

（1）考慮顆粒旋轉和未考慮顆粒旋轉時模擬得到的氣體速度場、顆粒速度場和平均體積分數(shù)的分布趨勢基本一致.

（2）考慮顆粒旋轉后，在垂直方向上床層中心位置處的顆粒平均速度增大，靠近壁面兩側的顆粒平均速度減小，顆粒水平方向的平均速度較未考慮顆粒旋轉時有所增大.

（3）與未考慮顆粒旋轉情況相比，考慮顆粒旋轉后床層膨脹高度有所增加，底部區(qū)域的床層平均空隙率增大.

（4）顆粒間的碰撞次數(shù)與床內(nèi)顆粒平均體積分數(shù)分布趨勢一致，顆粒碰撞次數(shù)隨著床高增加而減小，靠近壁面兩側的顆粒碰撞頻率略高于中心區(qū)域.

（5）床內(nèi)噴口中心左側的顆粒逆時針旋轉，右側的顆粒順時針旋轉，且顆粒的平均轉速隨床高的增加而減小，顆粒平均轉速的分布與顆粒水平方向平均速度的分布趨勢基本一致.

（6）顆粒的平均轉速不僅與顆粒間的碰撞次數(shù)有關，還與顆粒間的碰撞強度有關，并且密度大的顆粒的平均轉速相對較小，密度小的顆粒的平均轉速相對較大，直徑小的顆粒的平均轉速大于直徑大的顆粒.

［1］由長福，祁海鷹，徐旭常.煤粉顆粒所受Magnus力的數(shù)值模擬［J］.工程熱物理學報，2001，22（5）：625-628. YOU Changfu，QI Haiying，XU Xuchang.Numerical simulation of Magnus lift on a coal particle［J］.Journal of Engineering Thermophysics，2001，22（5）：625-628.

［2］吳學成，王勤輝，駱仲泱，等.高速數(shù)字攝影應用于流化床內(nèi)顆粒旋轉特性的測試［J］.中國電機工程學報，2005，25（17）：72-77. WU Xuecheng，WANG Qinhui，LUO Zhongyang，et al.Measurement of particle rotation in CFB riser with high-speed digital video camera techique［J］.Proceedings of the CSEE，2005，25（17）：72-77.

［3］駱仲泱，吳學成，王勤輝，等.循環(huán)流化床中顆粒旋轉特性［J］.化工學報，2005，25（10）：1869-1874. LUO Zhongyang，WU Xuecheng，WANG Qinhui，et al.Particle rotation characteristics in CFB riser［J］.Journal of Chemical Industry and Engineering（China），2005，25（10）：1869-1874.

［4］高瓊.循環(huán)流化床內(nèi)部顆粒流動特性的試驗研究［D］.杭州：浙江大學，2005.

［5］袁竹林，馬明，徐益謙，等.離散數(shù)值模擬顆粒自轉對流化床內(nèi)顆粒流態(tài)的影響［J］.燃燒科學與技術，2001，7（3）：238-241. YUAN Zhulin，MA Ming，XU Yiqian，etal.Study on effect of particle rotation on fluidizing behavior using discrete numerical simulation method［J］.Journal of Combustion Science and Technology，2001，7（3）：238-241.

［6］郝振華，陳巨輝，白穎華，等.考慮顆粒旋轉的顆粒動力學模擬提升管氣固兩相流動特性［J］.高?；瘜W工程學報，2010，24（5）：776-782. HAO Zhenhua，CHEN Juhui，BAI Yinghua，etal.Numerical study on gas-solid flow in a riser using dynamic theory with consideration of particle rotation［J］.Journal of Chemical Engineering of Chinese Universities，2010，24（5）：776-782.

［7］SUN J，BATTAGLIA F.Hydrodynamic modeling of particle rotation for segregation in bubbling gas-fluidized beds［J］.Chemical Engineering Science，2006，61（5）：1470-1479.

［8］TRIESCH O，BOHNET M.Measurement and CFD prediction of velocity and concentration profiles in a decelerated gas-solids flow［J］.Powder Technology，2001，115（2）：101-113.

［9］HUILIN L，GIDASPOW D.Hydrodynamics of binary fluidization in a riser：CFD simulation using two granular temperatures［J］.Chemical Engineering Science，2003，58（16）：3777-3792.

［10］HUILIN L，YURONG H，GIDASPOW D.Hydrodynamic modelling of binary mixture in a gas bubbling fluidized bed using the kinetic theory of granular flow［J］.Chemical Engineering Science，2003，58（7）：1197-1205.

［11］XU B H，YU A B.Numerical simulation of the gassolid flow in a fluidized bed by combining discrete particle method with computational fluid dynamics［J］.Chemical Engineering Science，1997，52（16）：2785-2809.

［12］田鳳國，章明川，顧明言，等.流化床微觀混合過程的定量評價［J］.動力工程，2007，27（1）：45-49. TIAN Fengguo，ZHANG Mingchuan，GU Mingyan，etal.Quantitative evaluation of micro-mixing processes in fluidized beds［J］.Journal of Power Engineering，2007，27（1）：45-49.

［13］WU X，WANG Q，LUO Z，etal.Theoretical and experimental investigations on particle rotation speed in CFB riser［J］.Chemical Engineering Science，2008，63（15）：3979-3987.

［14］李斌，紀律.流化床爐內(nèi)顆?；旌系碾x散單元法數(shù)值模擬［J］.中國電機工程學報，2012，32（20）：42-48. LI Bin，JI Lü.Numerical simulation of particle mixing in circulating fluidized bed with discrete element method［J］.Proceedings of the CSEE，2012，32（20）：42-48.

［15］紀律.循環(huán)流化床流動及磨損特性的DEM 數(shù)值模擬［D］.保定：華北電力大學，2011.