立式上行管內(nèi)液固兩相流中橢球形顆粒分布和運(yùn)動(dòng)特性

彭德其， 王依然， 侯家鑫， 俞天蘭， 吳淑英， 王志平

(1. 湘潭大學(xué)機(jī)械工程學(xué)院， 湖南湘潭411105； 2. 湖南工業(yè)大學(xué)機(jī)械工程學(xué)院, 湖南株洲412007；3. 湘潭鍋爐有限責(zé)任公司, 湖南湘潭411202)

液固兩相流換熱器具有強(qiáng)化傳熱及防垢除垢特性，被廣泛應(yīng)用于化工、 海水淡化等工業(yè)領(lǐng)域。由于目前實(shí)驗(yàn)獲取顆粒運(yùn)動(dòng)特性的方法對(duì)實(shí)驗(yàn)設(shè)備要求較高， 實(shí)驗(yàn)成本受到限制,因此，越來(lái)越多學(xué)者選擇建立各種數(shù)值模型解決該問(wèn)題，其中CFD-DEM模型所需要的經(jīng)驗(yàn)參數(shù)較少，對(duì)顆粒尺度微觀信息獲取較簡(jiǎn)便，有利于研究液固兩相流動(dòng)的微觀機(jī)理。

由于球形顆粒簡(jiǎn)單，過(guò)去以球形顆粒為研究對(duì)象，使用CFD方法進(jìn)行研究[1]，如Kerst等[2]通過(guò)對(duì)液固流化床中兩相流實(shí)驗(yàn)與CFD-DEM耦合模擬對(duì)比， 獲得了固體顆粒在流化床中分布形態(tài)和運(yùn)動(dòng)信息；Wang等[3]的研究表明，豎直液固流化床中流體進(jìn)口速度增加會(huì)增大顆粒在軸向方向的速度；Razzak等[4]發(fā)現(xiàn)液固兩相中直徑較大顆粒有更好的均勻性。在實(shí)際應(yīng)用中顆粒形狀多為非球形[5]，由于非球形顆粒建模難度較大，對(duì)異形顆粒研究較少。針對(duì)這種情況，提出多種模型用于描述非球形顆粒，包括超橢球模型[6-7]、 復(fù)合球模型[8]、 多面體模型[9]、 真實(shí)形狀模型[10]等。借助于這些顆粒模型，近年來(lái)對(duì)非球形顆粒系統(tǒng)的CFD-DEM數(shù)值模擬研究取得了較大進(jìn)展[11-12]。Ma等[13]建立非結(jié)構(gòu)化CFD網(wǎng)格和棒狀顆粒與復(fù)雜幾何接觸檢測(cè)的CFD-DEM模型,對(duì)復(fù)雜幾何形狀的氣-固流化床中棒狀顆粒流態(tài)化進(jìn)行模擬，發(fā)現(xiàn)具有較高的效率和精度。Campbell等[14]研究表明：與普通球形顆粒相比，橢球形顆粒的長(zhǎng)徑比對(duì)剪切流的相變變化有重要影響，橢球形顆粒間更容易形成影響顆粒周邊流場(chǎng)變化的“力鏈”。

上述文獻(xiàn)研究對(duì)象大部分是微米級(jí)粒徑球形顆粒及其在流化床流場(chǎng)中運(yùn)動(dòng)和分布特性，對(duì)大粒徑顆粒研究涉及不多。大粒徑顆粒相較于微米級(jí)顆粒，體積變大， 顆粒運(yùn)動(dòng)受阻，顆粒之間容易聚團(tuán)，而聚團(tuán)形成與破裂對(duì)液固兩相運(yùn)動(dòng)有較大影響。部分學(xué)者針對(duì)較大空間池式容器內(nèi)異形顆粒形狀對(duì)周圍流場(chǎng)影響進(jìn)行了仿真模擬[15-16]。雖然顆粒形狀是引起流化床內(nèi)流動(dòng)特性顯著變化的重要參數(shù)，但對(duì)立式換熱管內(nèi)大粒徑異形顆粒研究較少。

圖1 模擬流程圖Fig.1 Simulation flow chart

本文中以液固兩相流換熱設(shè)備中立式上行換熱管為研究對(duì)象，采用CFD-DEM數(shù)值模擬方法，對(duì)從圓形到長(zhǎng)圓形粒子形狀的橢球形顆粒的長(zhǎng)徑比β和顆粒體積分?jǐn)?shù)α進(jìn)行數(shù)值模擬，并深入分析管內(nèi)橢球形顆粒的分布和運(yùn)動(dòng)行為，對(duì)進(jìn)一步闡述立式列管式換熱器中液固兩相流流動(dòng)機(jī)理和橢球形顆粒在工業(yè)換熱器中實(shí)際應(yīng)用具有重要意義。

1 研究方法

1.1 模型描述及方案設(shè)定

選用常見(jiàn)立式換熱管的結(jié)構(gòu)如圖1所示。以管底部圓心(0,0,0)位置為參考點(diǎn)，換熱管內(nèi)徑Din=32 mm， 換熱管高度L=2 000 mm。 不考慮壁厚影響，流體從豎管底部流入，從頂部流出，顆粒工廠位于管底部，流體域體積為1.608 5×10-3m3。橢球形顆粒長(zhǎng)徑比如圖2所示。橢球形顆粒是一個(gè)類球體，其中半長(zhǎng)徑為a，其他2個(gè)半短徑b、c相等，即b=c，橢球形長(zhǎng)徑比為半長(zhǎng)軸與半短軸之比，即β=a/b。選取體積與半徑1 mm的圓球形顆粒相等的不同長(zhǎng)徑比β的橢球形顆粒進(jìn)行研究。

圖2 橢球形顆粒長(zhǎng)徑比示意圖Fig.2 Schematic diagram of aspect ratio of ellipsoidal particles

圖3為ICEM軟件生成的立式直管網(wǎng)格示意圖。為兼顧計(jì)算精度與工作效率，對(duì)管壁面進(jìn)行邊界層加密處理，設(shè)置近壁面第1層網(wǎng)格尺寸0.05 mm，增長(zhǎng)率1.1，從外往里劃分15層，對(duì)流體區(qū)域使用六面體網(wǎng)格進(jìn)行劃分。同時(shí)，DEM和CFD中立式上行管模型網(wǎng)格文件相同。

圖3 立式上行管網(wǎng)格Fig.3 Vertical upline pipe grid

管內(nèi)使用水作為流體相，流體物性參數(shù)如表1所示。流體進(jìn)口流速選擇實(shí)際中廣泛應(yīng)用1.5 m/s。通過(guò)改變橢球形顆粒長(zhǎng)徑比β和顆粒體積分?jǐn)?shù)α，構(gòu)建了25組不同工況，如表2所示。選取5種不同長(zhǎng)徑比β橢球形顆粒，采用CFD-DEM的歐拉-拉格朗日耦合方法，顆粒體積分?jǐn)?shù)α要低于10%，考慮到過(guò)大顆粒入口濃度會(huì)產(chǎn)生較大流動(dòng)阻力，且會(huì)造成管道進(jìn)口堵塞及管內(nèi)壁嚴(yán)重磨蝕[17]，所以顆粒體積分?jǐn)?shù)α應(yīng)不大于5%。

表1 流體物性參數(shù)

表2 模擬工況

1.2 數(shù)值模擬設(shè)置

使用Fluent軟件對(duì)流體相進(jìn)行數(shù)值模擬，選用基于壓力瞬態(tài)求解器下標(biāo)準(zhǔn)k-ε湍流模型，選用近壁面增強(qiáng)壁面處理法處理。選用速度入口邊界條件，壓力出口邊界條件定義出口靜壓為0，參考?jí)毫?01 325 Pa，管壁處設(shè)置恒定壁溫為350 K，選擇非滑移壁面邊界條件。計(jì)算過(guò)程首先用一階迎風(fēng)格式的SIMPLE算法進(jìn)行初始流場(chǎng)的穩(wěn)態(tài)計(jì)算，然后在此基礎(chǔ)上使用二階迎風(fēng)格式的PISO算法進(jìn)行非穩(wěn)態(tài)計(jì)算。

壁面與顆粒之間的參數(shù)參照文獻(xiàn)[18-19]設(shè)置，見(jiàn)表3，顆粒材料為陶瓷，立式管材料為鋼。顆粒相在EDEM軟件中與顆粒、 壁面間均采用內(nèi)置無(wú)滑移Hertz-Mindlin接觸模型，重力方向?yàn)榱黧w流動(dòng)的相反方向。

表3 壁面與顆粒參數(shù)

1.3 模型驗(yàn)證

為了驗(yàn)證CFD-DEM耦合方法可靠性，采用Shokri等[20]在相同雷諾數(shù)和濃度下不同粒徑玻璃球顆粒實(shí)驗(yàn)結(jié)果進(jìn)行對(duì)應(yīng)建模及數(shù)值模擬。圖4為顆粒流體在管道中心和管壁處的平均速度分布，其中，r表示管內(nèi)任意位置到管中心的徑向距離，R表示管的半徑。模擬與實(shí)驗(yàn)的誤差范圍在7%～10%之間，認(rèn)為誤差可接受，說(shuō)明在立式上行換熱管內(nèi)采用CFD-DEM耦合方法模擬橢球形顆粒運(yùn)動(dòng)是可行的。

圖4 數(shù)值模擬與實(shí)驗(yàn)結(jié)果對(duì)比圖Fig.4 Comparison of numerical simulation and experimental results

2 結(jié)果與分析

2.1 顆粒分布位置周期性波動(dòng)分析

為了研究在立式上行換熱管內(nèi)顆粒自下而上運(yùn)動(dòng)中顆粒位置分布，沿流體和顆粒運(yùn)動(dòng)方向，考慮顆粒自身尺寸，選取Z=(400±5)、 (800±5)、 (1 200±5)、 (1 600±5) mm等4個(gè)區(qū)域的顆粒位置信息。因?yàn)橐汗虄上噙\(yùn)動(dòng)具有隨機(jī)性，一定時(shí)間后管內(nèi)會(huì)達(dá)到動(dòng)態(tài)穩(wěn)定，動(dòng)態(tài)穩(wěn)定的表現(xiàn)可以通過(guò)顆粒質(zhì)量流量來(lái)觀測(cè)。選取顆粒體積分?jǐn)?shù)α=2%，長(zhǎng)徑比β=1.5的橢球形顆粒在監(jiān)測(cè)段(Z=1 800～2 000 mm)內(nèi)的質(zhì)量流量隨時(shí)間變化進(jìn)行分析，結(jié)果如圖5表示。由圖可以觀察到，t≥1.5 s后達(dá)到動(dòng)態(tài)穩(wěn)定，顆粒質(zhì)量流量在一定數(shù)值范圍內(nèi)波動(dòng)，選取波動(dòng)曲線中的中軸線、波峰和波谷t1、t2和t33個(gè)時(shí)刻的顆粒質(zhì)量流量進(jìn)行研究。

圖5 質(zhì)量流量變化圖Fig.5 Change diagram of mass flow rate

圖6—8分別為上述3個(gè)時(shí)刻的顆粒在立式換熱管橫截面(即XY平面， mm)上投影， 圖中每一個(gè)圓點(diǎn)表示一個(gè)顆粒， 在Z=(400±5) mm區(qū)域， 在XY平面上顆粒分布較均勻， 顆粒進(jìn)入Z=(800±5) mm區(qū)域，觀察到近壁面區(qū)域顆粒增加，在Z=(1 200～2 000) mm顆粒運(yùn)動(dòng)后期，顆粒幾乎堆積在管壁附近，管中心有少量零星顆粒，說(shuō)明在顆粒隨流體自下而上運(yùn)動(dòng)過(guò)程中，均勻進(jìn)入管內(nèi)顆粒逐漸受到流體作用向近壁面附近移動(dòng)。選取周期不同時(shí)刻t1、t2和t3，顆粒在相同軸向位置變化趨勢(shì)相同。動(dòng)態(tài)穩(wěn)定后的其他工況橢球形顆粒分布與圖6—8相似，在此不一一闡述。

a)Z=(400±5) mmb)Z=(800±5) mmc)Z=(1 200±5) mmd)Z=(1 600±5) mm圖6 t1時(shí)刻立式管內(nèi)不同高度的顆粒分布圖Fig.6 Particle distribution at different heights in vertical pipes at t1

a)Z=(400±5) mmb)Z=(800±5) mmc)Z=(1 200±5) mmd)Z=(1 600±5) mm圖7 t2時(shí)刻立式管內(nèi)不同高度的顆粒分布圖Fig.7 Particle distribution at different heights in vertical pipes at t2

a)Z=(400±5) mmb)Z=(800±5) mmc)Z=(1 200±5) mmd)Z=(1 600±5) mm圖8 t3時(shí)刻立式管內(nèi)不同高度的顆粒分布圖Fig.8 Particle distribution at different heights in vertical pipes at t3

2.2 顆粒徑向相對(duì)體積分?jǐn)?shù)分布分析

將立式換熱管沿Z軸方向，以管圓心為中心選擇7個(gè)不同半徑(r=2、 4、 6、 8、 10、 12、 14 mm)的圓柱形區(qū)域，以此來(lái)研究經(jīng)充分發(fā)展后換熱管徑向方向上顆粒相對(duì)體積分?jǐn)?shù)變化規(guī)律(見(jiàn)圖9)。對(duì)近壁面區(qū)域顆粒相對(duì)體積分?jǐn)?shù)分布分析，由于顆粒自身體積較大，因此選擇r=14 mm截面看作近壁面區(qū)域。

圖9 區(qū)域劃分示意圖Fig.9 Diagram of regional division

圖10為不同顆粒體積分?jǐn)?shù)下的顆粒徑向相對(duì)體積分?jǐn)?shù)分布圖。如圖分析得到，1)經(jīng)充分發(fā)展階段后，徑向方向上顆粒相對(duì)體積分?jǐn)?shù)主要集中在近壁區(qū)域，管中心區(qū)域顆粒相對(duì)體積分?jǐn)?shù)較低，說(shuō)明立式管內(nèi)顆粒受流體徑向作用力影響發(fā)生遷移；2)在r為2～8 mm區(qū)域，顆粒相對(duì)體積分?jǐn)?shù)變化不明顯，而在r為8～14 mm區(qū)域，越靠近壁面區(qū)域，顆粒相對(duì)體積分?jǐn)?shù)增加的幅度越來(lái)越大，這是由于絕大部分管中心顆粒受到流體的徑向作用力，逐漸向管壁面做遷移運(yùn)動(dòng)，最終導(dǎo)致近壁面區(qū)域顆粒相對(duì)體積分?jǐn)?shù)最大；3)在低顆粒體積分?jǐn)?shù)(α=1%、 2%、 3%)時(shí)，顆粒相對(duì)體積分?jǐn)?shù)在r為10～14 mm區(qū)域有較大的增長(zhǎng)變化，而高顆粒體積分?jǐn)?shù)(α=4%、 5%)中有較大增長(zhǎng)的顆粒相對(duì)體積分?jǐn)?shù)則提前至r為8～14 mm。觀察圖10中各曲線斜率可知，隨著顆粒體積分?jǐn)?shù)α由1%增長(zhǎng)到5%，在r=10 mm處的曲線斜率越來(lái)越大，這是由于隨管內(nèi)顆粒數(shù)量越來(lái)越多，近壁面區(qū)域(r=14 mm)顆粒數(shù)量逐漸達(dá)到飽和，其他受到流體作用力的顆粒慢慢向r=10～12 mm區(qū)域堆積。

不同顆粒長(zhǎng)徑比β條件下，近壁區(qū)域(r=14 mm)顆粒相對(duì)體積分?jǐn)?shù)隨顆粒體積分?jǐn)?shù)α變化如圖11所示。由圖可知，隨著顆粒體積分?jǐn)?shù)α增大，顆粒在近壁區(qū)域顆粒相對(duì)體積分?jǐn)?shù)均呈上升趨勢(shì)，但是增長(zhǎng)速率逐漸變小。這是因?yàn)殡S著顆粒體積分?jǐn)?shù)增加，進(jìn)入管內(nèi)顆粒數(shù)量變多，顆粒在流體徑向力作用下向近壁面區(qū)域遷移，導(dǎo)致近壁面區(qū)域顆粒相對(duì)體積分?jǐn)?shù)增大，然而，當(dāng)近壁面區(qū)域處顆粒數(shù)量達(dá)到飽和時(shí)，顆粒相對(duì)體積分?jǐn)?shù)值變化較小。

a)α=1%b)α=2%c)α=3%d)α=4%e)α=5%圖10 顆粒徑向相對(duì)體積分?jǐn)?shù)分布圖Fig.10 Radial relative volume fraction distribution of particles

圖11 截面r=14 mm處顆粒相對(duì)體積分?jǐn)?shù)對(duì)比圖Fig.11 Comparison chart of relative volume fraction of particles at section r=14 mm

長(zhǎng)徑比β=1.5的橢球形顆粒近壁區(qū)域相對(duì)體積分?jǐn)?shù)最大，分別比β=1.0、β=2.0、β=2.5和β=3.0顆粒的提高16.43%、4.82%、11.72%和19.47%，說(shuō)明近壁區(qū)域顆粒相對(duì)體積分?jǐn)?shù)與橢球形顆粒長(zhǎng)徑比β有關(guān)。在低顆粒體積分?jǐn)?shù)(1%～2%)，長(zhǎng)徑比β=1.0的球形顆粒與β=3.0的橢球形顆粒在近壁區(qū)域的顆粒相對(duì)體積分?jǐn)?shù)值相差不大，隨著顆粒體積分?jǐn)?shù)增大(3%～5%)，球形顆粒相對(duì)體積分?jǐn)?shù)大于β=3.0橢球形顆粒，說(shuō)明橢球形顆粒在近壁區(qū)域內(nèi)顆粒相對(duì)體積分?jǐn)?shù)值并非恒大于球形顆粒。

2.3 顆粒跟隨性分析

顆粒與流體在立式管內(nèi)運(yùn)動(dòng)過(guò)程分初始段、過(guò)渡段和充分發(fā)展段，在充分發(fā)展段，顆粒與流體運(yùn)動(dòng)較為穩(wěn)定?；谏鲜鲱w粒徑向濃度分布，在充分發(fā)展階段(Z=1 500～2 000 mm)，分別提取同一位置顆粒相和流體相瞬時(shí)軸向速度， 觀察橢球形顆粒長(zhǎng)徑比對(duì)液固兩相速度差的分布規(guī)律。 由于考慮顆粒具體尺寸， 所以顆粒參數(shù)在流體相應(yīng)位置Y方向(Y=0 mm)增加一定厚度(Y=-3～3 mm)， 其他方向(X=-16～16 mm、Z=1 500～2 000 mm)的兩相選取參數(shù)均相同。 采用r/R來(lái)表征液固兩相位置信息。 在r/R=±1處(管內(nèi)壁處)流體的速度為0， 而且無(wú)法得到該位置顆粒的速度信息， 為了研究同一時(shí)刻流體和顆粒的速度信息選擇r/R=-0.95～0.95的徑向范圍。

根據(jù)圖10顆粒徑向濃度分布圖可知，長(zhǎng)徑比β的顆粒體積分?jǐn)?shù)變化趨勢(shì)大致相同，因此選取β對(duì)徑向濃度分布影響差異最小的α=3%工況進(jìn)行分析。圖12為橢球形顆粒長(zhǎng)徑比β對(duì)液固兩相瞬時(shí)軸向速度分布的影響。

a)β=1.0b)β=1.5c)β=2.0d)β=2.5e)β=3.0圖12 顆粒與流場(chǎng)速度對(duì)比圖Fig.12 Velocity comparison between particles and flow field

由圖可知，1)長(zhǎng)徑比β顆粒對(duì)流體沿徑向方向的軸向速度分布不同，加入球形顆粒的流體軸向速度呈倒U形并呈管中心對(duì)稱分布，加入β>1.0橢球形顆粒流體軸向速度分布并不嚴(yán)格對(duì)稱，這也表明橢球形顆粒在流場(chǎng)中運(yùn)動(dòng)隨機(jī)性更大；2)橢球形顆粒與周圍流體速度差增大，說(shuō)明橢球形顆粒在流體內(nèi)速度波動(dòng)程度更大，即橢球形顆粒在流體中有更活躍運(yùn)動(dòng)狀態(tài)，從而降低了橢球形顆粒跟隨流體同步運(yùn)動(dòng)能力。

為了更好表征橢球形顆粒長(zhǎng)徑比對(duì)顆粒跟隨性影響，參照文獻(xiàn)[21]定義滑移系數(shù)η進(jìn)行量化評(píng)估。

η=Uf/Up，

(1)

式中：Uf為流體平均速度；Up為顆粒平均速度。

滑移系數(shù)η越大，表示周圍流場(chǎng)與顆粒速度差越大。在近壁面處，流體滿足無(wú)滑移邊界條件，而顆粒與壁面碰撞有較大動(dòng)量，導(dǎo)致在近壁面區(qū)域(r/R=-0.95～-0.7和r/R=0.7～0.95)顆粒速度大于流體速度，因此選擇管中心主要區(qū)域(r/R=-0.7～0.7)對(duì)液固兩相滑移系數(shù)分布規(guī)律進(jìn)行研究。

將立式管自下而上按高度為200 mm劃分為A—J等10個(gè)子流體區(qū)域，如圖13所示。圖14為上述工況下各子區(qū)域滑移系數(shù)變化，在立式管軸向中心區(qū)域，周圍流體速度恒大于顆粒速度，即滑移系數(shù)η>1。 進(jìn)入顆粒穩(wěn)定運(yùn)動(dòng)階段，β>1.0橢球形顆?；葡禂?shù)始終大于球形顆粒，因?yàn)闄E球形顆粒比球形顆粒在流體中的波動(dòng)范圍更大，其中β=1.5橢球形顆粒有最大的滑移系數(shù)，說(shuō)明在流體運(yùn)動(dòng)中長(zhǎng)徑比β=1.5橢球形顆粒比其他顆粒更加活躍。

圖13 區(qū)域劃分圖Fig.13 Regional division diagram圖14 各子區(qū)域滑移系數(shù)變化圖Fig.14 Variation diagram of slip coefficient in each sub-region

3 結(jié)論

采用CFD-DEM耦合方法研究了顆粒在立式換熱管內(nèi)分布與運(yùn)動(dòng)特性，得到如下結(jié)論。

1)在動(dòng)態(tài)穩(wěn)定周期內(nèi)，橫截面上顆粒徑向分布在任意時(shí)刻均相似，顆粒在立式換熱管下半部分橫截面上較均勻，顆粒隨著流體向上運(yùn)動(dòng)，管中心處顆粒數(shù)量逐漸減少；在換熱管上半部分，管中心處幾乎不存在顆粒，顆粒主要分布在管壁附近。

2)經(jīng)充分發(fā)展階段后，顆粒相對(duì)體積分?jǐn)?shù)沿徑向方向逐漸增加，管中心顆粒數(shù)量較少且變化不大，隨著顆粒體積分?jǐn)?shù)增加，近壁區(qū)域顆粒相對(duì)體積分?jǐn)?shù)先逐漸增大，達(dá)到飽和后趨于穩(wěn)定。顆粒體積分?jǐn)?shù)α為1%～3%時(shí)，顆粒相對(duì)體積分?jǐn)?shù)明顯增長(zhǎng)位置在r=10～14 mm區(qū)域，顆粒體積分?jǐn)?shù)α為4%～5%時(shí)有明顯增長(zhǎng)的顆粒相對(duì)體積分?jǐn)?shù)位置則提前至離管中心較近的r=8～14 mm區(qū)域。長(zhǎng)徑比β=1.5橢球形顆粒在近壁區(qū)域的相對(duì)體積分?jǐn)?shù)最大，但橢球形顆粒在近壁區(qū)域內(nèi)顆粒相對(duì)體積分?jǐn)?shù)值并非恒大于球形顆粒。

3)截面上顆粒和流體速度均呈現(xiàn)倒U型近似軸對(duì)稱分布，在r/R=-0.7～0.7的管內(nèi)截面區(qū)域周圍流體速度均大于顆粒速度；在靠近壁面區(qū)域流體速度小于顆粒速度；橢球形顆?；葡禂?shù)均大于球形顆粒的，其中β=1.5橢球形顆粒的滑移系數(shù)最大。