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

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

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

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

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

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

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

1 研究方法

1.1 模型描述及方案設定

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

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

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

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

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

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

表2 模擬工況

1.2 數(shù)值模擬設置

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

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

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

1.3 模型驗證

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

圖4 數(shù)值模擬與實驗結果對比圖Fig.4 Comparison of numerical simulation and experimental results

2 結果與分析

2.1 顆粒分布位置周期性波動分析

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

圖5 質量流量變化圖Fig.5 Change diagram of mass flow rate

圖6—8分別為上述3個時刻的顆粒在立式換熱管橫截面(即XY平面， mm)上投影， 圖中每一個圓點表示一個顆粒， 在Z=(400±5) mm區(qū)域， 在XY平面上顆粒分布較均勻， 顆粒進入Z=(800±5) mm區(qū)域，觀察到近壁面區(qū)域顆粒增加，在Z=(1 200～2 000) mm顆粒運動后期，顆粒幾乎堆積在管壁附近，管中心有少量零星顆粒，說明在顆粒隨流體自下而上運動過程中，均勻進入管內(nèi)顆粒逐漸受到流體作用向近壁面附近移動。選取周期不同時刻t1、t2和t3，顆粒在相同軸向位置變化趨勢相同。動態(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時刻立式管內(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時刻立式管內(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時刻立式管內(nèi)不同高度的顆粒分布圖Fig.8 Particle distribution at different heights in vertical pipes at t3

2.2 顆粒徑向相對體積分數(shù)分布分析

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

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

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

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

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

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

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

2.3 顆粒跟隨性分析

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

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

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

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

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

η=Uf/Up，

(1)

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

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

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

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

3 結論

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

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

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

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