蔡杰，鐘文琪，袁竹林

(1東南大學(xué)能源與環(huán)境學(xué)院，江蘇 南京 210096；2南京師范大學(xué)能源與機(jī)械工程學(xué)院，江蘇 南京 210042)

引言

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

1 數(shù)學(xué)模型的構(gòu)建

1.1 數(shù)學(xué)模型建立

本文將剛體動(dòng)力學(xué)理論與經(jīng)典的球形粒子氣固兩相流雙向耦合模型結(jié)合，建立細(xì)長(zhǎng)顆粒與氣相場(chǎng)間的雙向耦合關(guān)系。基于勒讓德-高斯積分法將細(xì)長(zhǎng)顆粒沿軸向離散成數(shù)個(gè)離散分段；獨(dú)立計(jì)算各離散分段在流場(chǎng)中所受的曳力，同時(shí)記錄該曳力的負(fù)值及該離散分段的體積；矢量合成各離散分段所受的曳力從而得到細(xì)長(zhǎng)顆粒所受的合曳力，同時(shí)計(jì)算細(xì)長(zhǎng)顆粒所受的力矩；根據(jù)歐拉動(dòng)力學(xué)方程，計(jì)算細(xì)長(zhǎng)顆粒的平動(dòng)及轉(zhuǎn)動(dòng)。同時(shí)，對(duì)于每個(gè)流場(chǎng)網(wǎng)格，根據(jù)所有顆粒離散分段所占的體積份額獲得網(wǎng)格中流場(chǎng)的空隙率，矢量合成網(wǎng)格內(nèi)所有離散分段對(duì)其施加的反作用力作為時(shí)均N-S方程的受力源項(xiàng)；再基于κ-ε模型耦合關(guān)聯(lián)式，分別計(jì)算κ方程及ε方程中的源項(xiàng)。從而構(gòu)建起細(xì)長(zhǎng)顆粒與流場(chǎng)之間的雙向耦合關(guān)系。

1.1.1 流場(chǎng)數(shù)學(xué)模型連續(xù)方程

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

時(shí)均N-S方程

式中，vsi為細(xì)長(zhǎng)顆粒離散分段的速度，m·s-1；m為流場(chǎng)網(wǎng)格中細(xì)長(zhǎng)顆粒離散分段的數(shù)量；μ、μt分別為空氣的黏度系數(shù)和湍流黏度系數(shù)，Pa·s；Ksg為氣相場(chǎng)與細(xì)長(zhǎng)顆粒離散分段之間的動(dòng)量交換系數(shù)[17]，此處采用Wen等[18]提出的經(jīng)驗(yàn)關(guān)聯(lián)式。

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

式中，dA為不規(guī)則顆粒的表面積等效直徑，m；c為不規(guī)則顆粒的表面球形度；Res為Reynolds數(shù)；ρs為細(xì)長(zhǎng)顆粒離散分段的密度，kg·m-3。

μt由式（5）計(jì)算

湍動(dòng)能方程

湍流耗散率方程

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

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

其中，b由式（11）計(jì)算

式中，CV=0.5，為附加質(zhì)量系數(shù)；ηsg則由式（12）計(jì)算

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

這里θ是細(xì)長(zhǎng)顆粒離散分段平均速度和平均相對(duì)速度的夾角。

式中，Lt為紊流渦的長(zhǎng)度標(biāo)尺，定義為

τF,sg定義如下

這里σκ=1，σε=1.3，C1ε=1.44，C2ε=1.92，C3ε=1.2，Cμ=0.09，均為經(jīng)驗(yàn)常數(shù)。該湍流源項(xiàng)模型能夠很好地滿足流場(chǎng)計(jì)算的穩(wěn)定性要求。

1.1.2 離散相數(shù)學(xué)模型文獻(xiàn)[16]中對(duì)細(xì)長(zhǎng)顆粒沿長(zhǎng)度方向離散采用均分?jǐn)?shù)段的方法，但這樣處理從數(shù)學(xué)上來(lái)講具有不合理性，會(huì)導(dǎo)致一定程度的誤差。因此，本文對(duì)之前的模型進(jìn)行了改進(jìn)，基于勒讓德-高斯積分法，沿長(zhǎng)度方向?qū)?xì)長(zhǎng)顆粒進(jìn)行數(shù)值離散。

模型中考慮到的力有氣流對(duì)細(xì)長(zhǎng)顆粒的曳力f和重力mg。細(xì)長(zhǎng)顆粒在流場(chǎng)中所受的力為細(xì)長(zhǎng)顆粒長(zhǎng)度方向上各點(diǎn)受力的積分。細(xì)長(zhǎng)顆粒在流場(chǎng)中的曳力f及定坐標(biāo)系力矩Mf為

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

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

Fi(si) 的計(jì)算公式為

每個(gè)離散分段的速度為

其中，vs0為細(xì)長(zhǎng)顆粒的質(zhì)心速度，m·s-1；I為由動(dòng)坐標(biāo)系向定坐標(biāo)系轉(zhuǎn)換的余弦矩陣；ω為細(xì)長(zhǎng)顆粒在動(dòng)坐標(biāo)系中的角速度，rad·s-1。

細(xì)長(zhǎng)顆粒的位移（平動(dòng)）及轉(zhuǎn)動(dòng)方程見(jiàn)文獻(xiàn)[16]。

式中，s為細(xì)長(zhǎng)顆粒的位移，m；m為細(xì)長(zhǎng)顆粒的質(zhì)量，kg；g為重力加速度，m·s-2。

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

離散分段i和同一網(wǎng)格內(nèi)其他細(xì)長(zhǎng)顆粒離散分段的碰撞概率為

式中，n為真實(shí)顆粒數(shù)；N為取樣顆粒數(shù)；Gij為細(xì)長(zhǎng)顆粒分段i和j的相對(duì)速度，m·s-1；Δt為時(shí)間步長(zhǎng)。

在細(xì)長(zhǎng)顆粒離散分段i總碰撞概率Pi小于1時(shí)，利用隨機(jī)數(shù)R(0j/N-Pij，則發(fā)生細(xì)長(zhǎng)顆粒間的碰撞。

限于篇幅，對(duì)兩細(xì)長(zhǎng)顆粒間的剛性碰撞模型此處不做詳細(xì)介紹，見(jiàn)文獻(xiàn)[16]。

細(xì)長(zhǎng)顆粒與壁面的處理采用壁面碰撞處理，壁面碰撞處理同樣基于剛體碰撞理論，并將墻壁設(shè)為質(zhì)量無(wú)窮大、速度為零的剛體，具體模型可見(jiàn)文獻(xiàn)[22]。

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

圖1 程序計(jì)算流程 Fig.1 Computational flow chart of program

連續(xù)相-細(xì)長(zhǎng)顆粒相雙向耦合兩相流程序全部由c++實(shí)現(xiàn)。采用Dell precision T5610工作站（內(nèi)存128G 雙E5-2637處理器）進(jìn)行計(jì)算。由于完成連續(xù)相的迭代計(jì)算和細(xì)長(zhǎng)顆粒的碰撞判斷及計(jì)算均需較長(zhǎng)時(shí)間，因此完成一個(gè)算例需超過(guò)10 d時(shí)間。

1.2 數(shù)值模擬參數(shù)及條件

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

表1 實(shí)際的工作條件 Table 1 Actual conditions

表2 典型的計(jì)算區(qū)域 Table 2 Typical calculation regions

2 數(shù)值模擬結(jié)果與分析

2.1 細(xì)長(zhǎng)顆粒流化數(shù)值模擬圖

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

圖2 提升管內(nèi)不同時(shí)刻細(xì)長(zhǎng)顆粒流化 運(yùn)動(dòng)模擬圖（5 m·s-1） Fig.2 Simulation of fluidization of slender particles in riser at different time (5 m·s-1)

圖3 提升管內(nèi)不同時(shí)刻的細(xì)長(zhǎng)顆粒流化 運(yùn)動(dòng)照片（5 m·s-1） Fig.3 Photos of fluidization of slender particles in riser at different time (5 m·s-1)

2.2 細(xì)長(zhǎng)顆粒徑向濃度分布

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

圖4 提升管內(nèi)水平中心軸線上不同時(shí)間時(shí)的 細(xì)長(zhǎng)顆粒濃度分布 Fig.4 Number concentration distribution of slender particles along horizontal center axis in riser at different time

2.3 流場(chǎng)內(nèi)細(xì)長(zhǎng)顆粒的取向分布

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

2.4 流場(chǎng)內(nèi)空隙率分布

圖5 流場(chǎng)內(nèi)細(xì)長(zhǎng)顆粒的取向分布 Fig.5 Orientation distribution of slender particles in riser

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

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

2.5 流場(chǎng)內(nèi)靜壓強(qiáng)分布

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

圖7 提升管內(nèi)水平中心軸線上不同時(shí)刻的壓強(qiáng)分布 Fig.7 Pressure distribution of local flow field along horizontal center axis in riser at different time

2.6 湍流速度的分布特征

圖8顯示了提升管內(nèi)1.63 m和2.55 m高度，沿中心軸線方向上，0.48、0.51、0.54和0.57 s時(shí)刻的當(dāng)?shù)赝牧魉俣确植肌膱D8中可以看出，湍流場(chǎng)速度分布大體呈中間水平，近壁處迅速下降趨勢(shì)。流場(chǎng)中間區(qū)域的速度梯度較小是由于湍流的擾動(dòng)作用，而近壁處則由于層流黏性支層的影響，速度值迅速下降，直至與壁面速度相同。并且，在有細(xì)長(zhǎng)顆粒停留的流場(chǎng)區(qū)域，流場(chǎng)速度值會(huì)有突降，這部分突降的速度值被轉(zhuǎn)移給細(xì)長(zhǎng)顆粒（離散分段），以使它們獲得向上運(yùn)動(dòng)的能量。當(dāng)?shù)亓鲌?chǎng)停留的細(xì)長(zhǎng)顆粒（離散分段）量越多，當(dāng)?shù)亓鲌?chǎng)的空隙率越小，速度值下降也越明顯。由于每個(gè)流場(chǎng)網(wǎng)格的流場(chǎng)參數(shù)變化會(huì)受到上下前后左右流場(chǎng)的影響，而軸線方向上的當(dāng)?shù)亓鲌?chǎng)上下前后左右的流場(chǎng)是不一樣的，此圖上沿軸線方向的速度值并不是嚴(yán)格的水平線，而是有一定程度的上下波動(dòng)。

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

3 結(jié)論

本文針對(duì)當(dāng)前高Reynolds數(shù)條件下細(xì)長(zhǎng)顆粒氣固兩相流模型發(fā)展尚不完善的現(xiàn)實(shí)情況，基于κ-ε模型的球形粒子-流場(chǎng)氣固間耦合關(guān)聯(lián)式，構(gòu)建了細(xì)長(zhǎng)顆粒-湍流場(chǎng)氣固雙向耦合模型，并采用該氣固雙向耦合模型對(duì)細(xì)長(zhǎng)顆粒氣固兩相流場(chǎng)進(jìn)行了數(shù)值模擬研究。研究發(fā)現(xiàn)，在此提升管內(nèi)，細(xì)長(zhǎng)顆粒氣固兩相流場(chǎng)有以下幾個(gè)特點(diǎn)。

（1）細(xì)長(zhǎng)顆粒在提升管內(nèi)向上運(yùn)動(dòng)過(guò)程中伴隨有較明顯的水平遷移，由中心向四周區(qū)域遷移，但細(xì)長(zhǎng)顆粒數(shù)量濃度最大區(qū)域并不是近壁區(qū)，而是靠近近壁區(qū)區(qū)域。

（2）細(xì)長(zhǎng)顆粒在流化運(yùn)動(dòng)過(guò)程中具有非常明顯的取向選擇性，細(xì)長(zhǎng)顆粒多數(shù)以其軸趨于與主流速度平行的運(yùn)動(dòng)姿態(tài)向上運(yùn)動(dòng)。

（3）沿水平軸線方向上，流場(chǎng)四周區(qū)域的空隙率小于中心區(qū)域的空隙率，并且空隙率分布具有較強(qiáng)的隨機(jī)性。

（4）隨著提升管內(nèi)顆粒數(shù)量的大量增加，當(dāng)?shù)亓鲌?chǎng)的靜壓強(qiáng)會(huì)有明顯的下降，并且，空隙率越小的區(qū)域，當(dāng)?shù)仂o壓強(qiáng)會(huì)有突降。

（5）湍流場(chǎng)速度沿徑向呈中間均勻，而近壁處則由于黏性支層的原因而急速下降直至與壁面速度相同分布，并且，有細(xì)長(zhǎng)顆粒（離散分段）停留的區(qū)域當(dāng)?shù)厮俣葧?huì)有較為明顯的下降。

[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(應(yīng)用力學(xué)學(xué)報(bào)), 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) (化工學(xué)報(bào)), 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(中國(guó)電機(jī)工程學(xué)報(bào)), 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(彭正標(biāo)), 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(中國(guó)電機(jī)工程學(xué)報(bào)), 2008, 28(23): 71-74.