提升管內細長顆粒流化運動的氣固多向耦合模型的構建

蔡杰，鐘文琪，袁竹林

(1東南大學能源與環境學院，江蘇 南京 210096；2南京師范大學能源與機械工程學院，江蘇 南京 210042)

引言

細長顆粒的循環流化在工業生產中具有廣泛的應用背景，如生物質秸稈在循環流化床中的燃燒、煙絲在循環流化床中的干燥或加濕、藥丸在流化床內的干燥成型等[1-3]。因此，研究細長顆粒的循環流化特性具有重要的意義。相比球形顆粒的三自由度，細長顆粒具有六自由度，包括位置矢量和姿態矢量。由剛體動力學理論可知，剛體的姿態矢量用一組歐拉角(進動角、章動角和自轉角) 唯一表示。細長顆粒的氣固兩相流模型研究主要著眼于細長顆粒的受力、運動（平動與轉動）模型的構建，細長顆粒間的碰撞模型構建以及細長顆粒與流場之間耦合關系的構建等。而其運動研究則主要著眼于細長顆粒在沉降、上升或運動過程中的取向分布、濃度分布、旋轉速度及位移等特征。細長顆粒的氣固兩相流研究主要采用歐拉-拉格朗日方法，流場作為連續相，采用歐拉方法處理，而細長顆粒則作為離散相，采用拉格朗日方法追蹤其受力、運動及碰撞[4-7]。在早期階段，細長顆粒與流場之間都是單向耦合，即僅考慮流場對細長顆粒的曳力，而不考慮細長顆粒的存在對流場特征的影響[4]。但事實上，即使是稀相流，細長顆粒的存在也會對流場產生較大的影響，因此，構建細長顆粒與流場[尤其是高Reynolds數（Re）湍流]之間的雙向耦合關系具有極其重要的意義[8-10]。拉格朗日時間尺度與湍動能（κ）-湍流耗散率（ε）的關系的建立，為固粒相-流場雙向耦合模型的構建奠定了基礎[11]。近年來，時均Navier(N)- Stokes（S）方程直接模擬方法被應用于球形固粒在湍流流場中的擴散研究[12-13]，但卻幾乎未見被應用于細長顆粒在湍流場中的擴散研究[14-15]。本文在文獻[16]的基礎上，通過引入拉格朗日時間尺度與κ-ε的關系，構建細長顆粒-湍流之間的雙向耦合關系，并改進細長顆粒間碰碰模型，從而構建起高Re條件下細長顆粒-湍流多向耦合模型。并且，采用此模型對某提升管內的細長顆粒-湍流氣固兩相流場進行了數值模擬研究。

1 數學模型的構建

1.1 數學模型建立

本文將剛體動力學理論與經典的球形粒子氣固兩相流雙向耦合模型結合，建立細長顆粒與氣相場間的雙向耦合關系。基于勒讓德-高斯積分法將細長顆粒沿軸向離散成數個離散分段；獨立計算各離散分段在流場中所受的曳力，同時記錄該曳力的負值及該離散分段的體積；矢量合成各離散分段所受的曳力從而得到細長顆粒所受的合曳力，同時計算細長顆粒所受的力矩；根據歐拉動力學方程，計算細長顆粒的平動及轉動。同時，對于每個流場網格，根據所有顆粒離散分段所占的體積份額獲得網格中流場的空隙率，矢量合成網格內所有離散分段對其施加的反作用力作為時均N-S方程的受力源項；再基于κ-ε模型耦合關聯式，分別計算κ方程及ε方程中的源項。從而構建起細長顆粒與流場之間的雙向耦合關系。

1.1.1 流場數學模型連續方程

式中，αg為空氣空隙率；ρg為空氣的密度，kg·m-3；vg為空氣的速度，m·s-1。

時均N-S方程

式中，vsi為細長顆粒離散分段的速度，m·s-1；m為流場網格中細長顆粒離散分段的數量；μ、μt分別為空氣的黏度系數和湍流黏度系數，Pa·s；Ksg為氣相場與細長顆粒離散分段之間的動量交換系數[17]，此處采用Wen等[18]提出的經驗關聯式。

式中，CDs為細長顆粒離散分段的曳力系數，kg·s-1；ds為離散分段i的體積等效直徑，m。由于是非球形顆粒，本文采用Sabine 等[19]提出的非規則顆粒曳力系數計算方法確定細長顆粒離散分段的曳力系數。

式中，dA為不規則顆粒的表面積等效直徑，m；c為不規則顆粒的表面球形度；Res為Reynolds數；ρs為細長顆粒離散分段的密度，kg·m-3。

μt由式（5）計算

湍動能方程

湍流耗散率方程

這里Gκ為湍動能的產生項，其意義與單相流體κ-ε模型中的意義相同；Iκ和Iε代表了固體相（s）對氣體相（g）的影響。Elgobashi[20]提出Iκ和Iε可采用如下形式

這里ksg是氣體相（g）、固體相（s）的速度的協方差，Simonin等[21]提出采用式（10）計算

其中，b由式（11）計算

式中，CV=0.5，為附加質量系數；ηsg則由式（12）計算

式中，τt為載能紊流渦的特征時間，定義如下

這里θ是細長顆粒離散分段平均速度和平均相對速度的夾角。

式中，Lt為紊流渦的長度標尺，定義為

τF,sg定義如下

這里σκ=1，σε=1.3，C1ε=1.44，C2ε=1.92，C3ε=1.2，Cμ=0.09，均為經驗常數。該湍流源項模型能夠很好地滿足流場計算的穩定性要求。

1.1.2 離散相數學模型文獻[16]中對細長顆粒沿長度方向離散采用均分數段的方法，但這樣處理從數學上來講具有不合理性，會導致一定程度的誤差。因此，本文對之前的模型進行了改進，基于勒讓德-高斯積分法，沿長度方向對細長顆粒進行數值離散。

模型中考慮到的力有氣流對細長顆粒的曳力f和重力mg。細長顆粒在流場中所受的力為細長顆粒長度方向上各點受力的積分。細長顆粒在流場中的曳力f及定坐標系力矩Mf為

式中，l為細長顆粒的長度，m；f(x)是細長顆粒長度方向上各點所受的流體曳力，N·m-1；F(s)是(-1,1)區間單位長度所受的流體曳力，N·m-1；x=[-l+l+(l+l)s]/4。根據勒讓德-高斯積分法，式（19）相應的求積公式為

其中，n=10為細長顆粒離散分段的數量；λi為高斯求積系數；Fi(si)為每個離散分段所受的流體曳力，N·m-1；li為每個離散分段的長度，m；r(si)為定坐標系中細長顆粒離散分段相對于細長顆粒質心位置的單位矢量，m。

Fi(si) 的計算公式為

每個離散分段的速度為

其中，vs0為細長顆粒的質心速度，m·s-1；I為由動坐標系向定坐標系轉換的余弦矩陣；ω為細長顆粒在動坐標系中的角速度，rad·s-1。

細長顆粒的位移（平動）及轉動方程見文獻[16]。

式中，s為細長顆粒的位移，m；m為細長顆粒的質量，kg；g為重力加速度，m·s-2。

本文在文獻[16]基礎上，對Nabu碰撞概率模型進行修正，從而建立基于細長顆粒離散分段的細長顆粒碰撞概率模型。其思路如下：對處于同一流體網格內的細長顆粒離散分段進行排序，并根據Nabu碰撞概率模型逐一判斷每個細長顆粒離散分段將會與其他哪個細長顆粒離散分段發生碰撞（同一細長顆粒不同離散分段之間的碰撞被排除），則其所歸屬的兩個細長顆粒將會發生碰撞，并且這兩個細長顆粒離散分段即為細長顆粒間碰撞的碰點。該處理方式無須額外再構建一套網格。在進行網格尺度確定時需要考慮兩方面的因素：一是同一細長顆粒上不同離散分段受力的差異；二是同一網格內不同細長顆粒離散分段碰撞。因此，本文在定義網格尺寸時，以最長細長顆粒長度的1/2作為網格尺寸。

離散分段i和同一網格內其他細長顆粒離散分段的碰撞概率為

式中，n為真實顆粒數；N為取樣顆粒數；Gij為細長顆粒分段i和j的相對速度，m·s-1；Δt為時間步長。

在細長顆粒離散分段i總碰撞概率Pi小于1時，利用隨機數R(0j/N-Pij，則發生細長顆粒間的碰撞。

限于篇幅，對兩細長顆粒間的剛性碰撞模型此處不做詳細介紹，見文獻[16]。

細長顆粒與壁面的處理采用壁面碰撞處理，壁面碰撞處理同樣基于剛體碰撞理論，并將墻壁設為質量無窮大、速度為零的剛體，具體模型可見文獻[22]。

1.1.3 模型計算連續采用基于同位網格的SIMPLE方法計算。先進行單相連續相的計算，獲得收斂解之后，開始細長顆粒受力、運動方程求解，細長顆粒的平動及轉動軌跡的計算。細長顆粒相連續計算一定的時間步數后，在進行細長顆粒受力運動計算的同時，計算細長顆粒離散分段在連續相網格中所占的空隙率、連續相網格所受反作用力源項以及κ-ε模型中的源項，并代入連續相計算從而獲得新的連續相。然后基于新的連續相繼續進行顆粒場的計算，從而實現連續相和細長顆粒相連續交替耦合求解。本文中設定的時間步數間隔為10個時間步，即每10個時間步重新更新連續相。圖1為計算流程。

圖1 程序計算流程 Fig.1 Computational flow chart of program

連續相-細長顆粒相雙向耦合兩相流程序全部由c++實現。采用Dell precision T5610工作站（內存128G 雙E5-2637處理器）進行計算。由于完成連續相的迭代計算和細長顆粒的碰撞判斷及計算均需較長時間，因此完成一個算例需超過10 d時間。

1.2 數值模擬參數及條件

本文數值模擬對象為一矩形提升管內的細長顆粒氣-固兩相流場，提升管的長寬高為0.5 m×0.5 m×3 m，采用的實驗細長顆粒為密度、外形均勻性良好的秸稈（火柴棒）。為最大程度還原加料環境，細長顆粒初始速度為0，角速度在動坐標系3個軸上投影為0～2 rad·s-1，3個歐拉角為0～π，以上均以隨機方式給出。每隔0.01 s，從提升管入口處放入一定數量的細長顆粒，本文中采用了4種長徑比的細長顆粒，并按一定級配由床層底部放入。實際的工作條件如表1所示。典型的計算區域如表2所示。

表1 實際的工作條件 Table 1 Actual conditions

表2 典型的計算區域 Table 2 Typical calculation regions

2 數值模擬結果與分析

2.1 細長顆粒流化數值模擬圖

受目前實驗技術能力與條件的限制，當前細長顆粒的流化運動實驗主要依靠高速攝像抓拍一些細長顆粒在流化運動過程中的濃度分布特性。數值實驗中，當細長顆粒到達提升管出口之后，提升管內細長顆粒的數量基本保持在10000個左右，相較實驗的顆粒量多，因此實驗照片與數值模擬圖略有差異。圖2和圖3分別按一定百分比配比的4種不同長徑比的細長顆粒物料在不同時刻的流化模擬圖和實驗照片。由圖2可以明顯看出，細長顆粒在上升運動過程中伴隨著強烈的旋轉，這與實驗照片結果吻合。盡管實驗顆粒量相對于模擬顆粒量較少，但是從實驗圖中可以明顯看出細長顆粒在流化運動過程中存在絮團現象，在不同時刻不同區域細長顆粒數量濃度明顯差異，這說明即便細長顆粒的數量濃度較低，數值建模也需考慮細長顆粒間的相互作用。在流化運動過程中，會有一些細長顆粒脫離顆粒群，以較高的速度向上運動，這在數值模擬圖與實驗照片中均有明顯的體現。這說明本文所建模型是合理的。

圖2 提升管內不同時刻細長顆粒流化 運動模擬圖（5 m·s-1） Fig.2 Simulation of fluidization of slender particles in riser at different time (5 m·s-1)

圖3 提升管內不同時刻的細長顆粒流化 運動照片（5 m·s-1） Fig.3 Photos of fluidization of slender particles in riser at different time (5 m·s-1)

2.2 細長顆粒徑向濃度分布

圖4為提升管內沿水平中心軸線上由中心四周方向的細長顆粒的數量濃度分布。由圖4可以看出，由中心向四周方向，細長顆粒的數量濃度由低到高逐漸增加，但到近壁處后又有明顯的下降。細長顆粒數量濃度最低區域在中心區域，而數量濃度最高區域在靠近近壁處區域。這個結論與實驗照片的定性結論（肉眼觀測）是一致的，與球形顆粒在流化床提升管內的數量濃度分布規律也是一致的。這說明與球形顆粒流化運動過程一樣，細長顆粒在提升管內上升運動過程中同樣伴隨著明顯的水平遷移，由中心向四周移動。原因在于流場中心區域的徑向速度大于四周區域的徑向速度。并且，由于壁面的碰撞反彈等原因，細長顆粒最高數量濃度區域并非出現在近壁區，而是靠近近壁區區域。

圖4 提升管內水平中心軸線上不同時間時的 細長顆粒濃度分布 Fig.4 Number concentration distribution of slender particles along horizontal center axis in riser at different time

2.3 流場內細長顆粒的取向分布

由于細長顆粒在流化運動過程中的向上運動速度非常快，因此通過實驗手段獲得其流化過程中的取向分布特性是非常困難的。本文采用了一種簡化的辦法，通過測量實驗圖片中的細長顆粒的長度，將實驗圖片中與拍攝平面接近平行的細長顆粒取出來，測量其與z軸的夾角，從而獲取細長顆粒在流化運動過程中的取向分布特性。圖5顯示了整個提升管內細長顆粒的取向分布的數值模擬與實驗結果。從圖5中可以看出，兩者的結論大體是一致的。在流化運動過程中，大多數細長顆粒以其軸與流場主流速度方向呈較小的夾角向上運動。并且，以0°～15°夾角向上運動的細長顆粒數量最大。這是由于細長顆粒以對流場影響最小的形式（能量最小化）在流場內做流化運動。

2.4 流場內空隙率分布

圖5 流場內細長顆粒的取向分布 Fig.5 Orientation distribution of slender particles in riser

受目前細長顆粒兩相流實驗技術能力的限制，通過實驗手段獲取流化過程中流場的速度、壓強、空隙率等流場信息是非常困難的，而數值模擬卻提供了獲取這些流場信息的有效途徑。圖6反映了提升管內1.63 m和2.55 m高度，沿中心軸線方向上，0.48、0.51、0.54和0.57 s時刻的流場空隙率分布。從圖6可以看出，由于是非穩態過程，因此不同時刻的空隙率差異非常明顯。在1.63 m和2.55 m高度，沿中心軸線方向上，空隙率的分布非常不規則，但大體而言，中心區域的空隙率要大于四周其他區域，這與流化床內顆粒濃度分布特性基本對應，中心區域細長顆粒數量濃度低，四周區域的細長顆粒數量濃度高。空隙率最低的區域出現在中心區與近壁區之間的過渡區域，這說明該區域的細長顆粒的濃度最高。并且，在流化運動過程中，由于細長顆粒間的相互碰撞的影響，細長顆粒的濃度分布的隨機性非常強，甚至某時刻某位置處，如圖6所示，1.63 m高度上，0.51 s時刻及2.55 m高度上，0.48 s、0.51 s、0.54 s時刻，當地流場沒有細長顆粒（離散分段）的存在，因而該處的流體空隙率均為100%。這些充分說明細長顆粒的遷移以及流場參數隨時間變化過程中具有非常強的隨機性。因此，采用非穩態模型模擬細長顆粒的流化運動是合理的。

圖6 提升管內水平中心軸線上不同時刻的流體空隙率分布 Fig.6 Volume fraction distribution of local flow field along horizontal center axis in riser at different time

2.5 流場內靜壓強分布

圖7反映了提升管內1.63 m和2.55 m高度，沿中心軸線方向上，0.48、0.51、0.54和0.57 s時刻的壓強分布。從圖7中可以看出，隨著時間增加，流場的靜壓強總體呈逐漸下降趨勢，這是由于隨著時間的增加，流場內的細長顆粒的數量不斷增加，流場需要將更多的能量傳遞給當地的細長顆粒。與圖6 相對應，空隙率下降明顯的區域，當地靜壓強相對其他地區有更加明顯的突降，這也是由于當地流場需要傳遞更多的能量給處于當地流場的細長顆粒（離散分段），使它們獲得足夠的能夠向上運動的能量，細長顆粒（離散分段）量越多，靜壓強突降越明顯。在1.63 m高度上，0.51 s時刻及2.55 m高度，0.48、0.51、0.54 s時刻，由于沒有細長顆粒停留在當地，因此沿軸線方向上靜壓強變化不明顯。此外，由圖7可以看出，沿該軸線方向上的當地流場的靜壓強值并不是一條水平線，而是有一定程度的上下波動。這是由于當地流場受到上下前后左右流場的影響，而軸線方向上的當地流場上下前后左右的流場是不一樣的，因此，軸線上所有流場的靜壓強值并不是一條水平線。

圖7 提升管內水平中心軸線上不同時刻的壓強分布 Fig.7 Pressure distribution of local flow field along horizontal center axis in riser at different time

2.6 湍流速度的分布特征

圖8顯示了提升管內1.63 m和2.55 m高度，沿中心軸線方向上，0.48、0.51、0.54和0.57 s時刻的當地湍流速度分布。從圖8中可以看出，湍流場速度分布大體呈中間水平，近壁處迅速下降趨勢。流場中間區域的速度梯度較小是由于湍流的擾動作用，而近壁處則由于層流黏性支層的影響，速度值迅速下降，直至與壁面速度相同。并且，在有細長顆粒停留的流場區域，流場速度值會有突降，這部分突降的速度值被轉移給細長顆粒（離散分段），以使它們獲得向上運動的能量。當地流場停留的細長顆粒（離散分段）量越多，當地流場的空隙率越小，速度值下降也越明顯。由于每個流場網格的流場參數變化會受到上下前后左右流場的影響，而軸線方向上的當地流場上下前后左右的流場是不一樣的，此圖上沿軸線方向的速度值并不是嚴格的水平線，而是有一定程度的上下波動。

圖8 提升管內水平中心軸線上不同時刻的湍流速度分布 Fig.8 Turbulent velocity distribution of local flow field along horizontal center axis in riser at different time

3 結論

本文針對當前高Reynolds數條件下細長顆粒氣固兩相流模型發展尚不完善的現實情況，基于κ-ε模型的球形粒子-流場氣固間耦合關聯式，構建了細長顆粒-湍流場氣固雙向耦合模型，并采用該氣固雙向耦合模型對細長顆粒氣固兩相流場進行了數值模擬研究。研究發現，在此提升管內，細長顆粒氣固兩相流場有以下幾個特點。

（1）細長顆粒在提升管內向上運動過程中伴隨有較明顯的水平遷移，由中心向四周區域遷移，但細長顆粒數量濃度最大區域并不是近壁區，而是靠近近壁區區域。

（2）細長顆粒在流化運動過程中具有非常明顯的取向選擇性，細長顆粒多數以其軸趨于與主流速度平行的運動姿態向上運動。

（3）沿水平軸線方向上，流場四周區域的空隙率小于中心區域的空隙率，并且空隙率分布具有較強的隨機性。

（4）隨著提升管內顆粒數量的大量增加，當地流場的靜壓強會有明顯的下降，并且，空隙率越小的區域，當地靜壓強會有突降。

（5）湍流場速度沿徑向呈中間均勻，而近壁處則由于黏性支層的原因而急速下降直至與壁面速度相同分布，并且，有細長顆粒（離散分段）停留的區域當地速度會有較為明顯的下降。

[1] Bernstein O, Shapiro M.Direct determination of the orientation distribution function of cylindrical particles immersed in laminar and turbulent shear flows [j].Journal of Aerosol Science, 1994, 25(1): 113-136.

[2] Lin Jianzhong, Lin Jiang, Shi Xing.Research on the effect of cylinder particles on the turbulent properties in particulate flows [J].Applied Mathematics and Mechanics, 2002, 23(5): 483-488．

[3] Lin Jianzhong(林建忠), Wang Changbin(王常斌), Shi Xing(石興)．Effects of wall on the sedimentation of cylindrical particles in the Newtonian fluid [J].Chinese Journal of Applied Mathematics(應用力學學報), 2003, 20(4): 1-6.

[4] Hilton J E, Mason L R, Cleary P W.Dynamics of gas-solid fluidised beds with non-spherical particle geometry [J].Chemical Engineering Science, 2010, 65(5): 1584-1594.

[5] Yin C, Rosendahl L, K?r K, et al.Modelling the motion of cylindrical particles in a nonuniform flow [J].Chemical Engineering Science, 2003, 58(15): 3489-3498.

[6] Fan L, Mao Z S, Wang Y D.Numerical simulation of turbulent solid-liquid two-phase flow and orientation of slender particles in a stirred tank [J].Chemical Engineering Science, 2005, 60(24): 7045-7056.

[7] Wang Yelong.Numerical study on sedimentation of mutual-collision cylindrical particles in the vertical channel [J].Applied Mathematics and Mechanic, 2006, 27(7): 859-866.

[8] Gore R, Crowe C.Effect of particle size on modulating turbulent intensity [J].Int.J.Multiphase Flow, 1989, 15: 279-288.

[9] Paschkewitz J S, Dubief Y, Dimitropoulos C D, et al.Numerical simulation of turbulent drag reduction using rigid fibres [J].J.Fluid Mech., 2004, 518: 281-317.

[10] Zhang Weifeng, Lin Jianzhong.Research on the motion of particles in the turbulent pipe flow of fiber suspensions [J].Applied Mathematics and Mechanics, 2004, 25(7): 677-685.

[11] Shuen J, Chen L, Faeth G.Evaluation of a stochastic model of particle dispersion in a turbulent round jet [J].AIChE J., 1983, 29(1): 167-170.

[12] Matteo C, Mathiesen V.Numerical simulation of particulate flow by the Eulerian-Lagrangian and the Eulerian-Eulerian approach with application to a fluidized bed [J].Comput.Chem.Eng., 2005, 29(2): 291-304.

[13] Xiong Yuanquan(熊源泉), Yuan Zhunlin(袁竹林)，Zhang Mingyao(章名耀).Three-dimensional numerical simulation on conveying properties of gas-solid injector under pressurization [J].Journal of Chemical Industry and Engineering (China) (化工學報), 2004, 55(10): 1638-1643.

[14] Lin Jianzhong, Lin Jiang, Shi Xing.Research on the effect of cylinder particles on the turbulent properties in particulate flows [J].Applied Mathematics and Mechanics, 2002, 23(5): 483-488.

[15] Qi D W.Simulations of fluidization of cylindrical multiparticles in a three-dimensional space [J].Int J.Multiphase Flow, 2001, 27(1): 107-118.

[16] Cai Jie(蔡杰), Fan Fengxian(凡鳳仙), Wu Xuan(吳晅), et al.Effects of inter-particle collisions on the orientation distribution of slender particles in gas-solid flows [J].Proceedings of the CSEE(中國電機工程學報), 2008, 28(17): 66-69.

[17] Crowe C T, Troutt T R, Chung J N.Numerical models for two-phase turbulent flows [J].Annual Review of Fluid Mechanics, 1996, 28: 11-43.

[18] Wen C Y, Yu Y H.Mechanics of fluidization [J].Chem.Eng.Prog.Sump.Serie., 1966, 62: 100-113.

[19] Sabine T C, Michael G, Efstathios E M.Drag coefficients of irregularly shaped particles [J].Powder Technology, 2004, 139: 21-32.

[20] Elgobashi S.Particle-laden turbulent flows: direct simulation and closure models [J].Applied Scientific Research, 1991, 48: 301-314.

[21] Simonin O, Deutsch E, Miniter J P.Eulerian prediction of the fluid/particle correlated motion in turbulent two-phase flows [J].Applied Scientific Research, 1993, 51: 275-283.

[22] Cai Jie(蔡杰), Peng Zhengbiao(彭正標), Wu Xuan(吳晅), et al.Simulation study on the effects of wall on the fluidizing features of slender particles in gas-solid flows [J].Proceedings of the CSEE(中國電機工程學報), 2008, 28(23): 71-74.