埋管流化床內(nèi)濕顆粒流動及混合特性的CFD-DEM數(shù)值模擬

劉道銀，宋誠驍，王錚，馬吉亮，陳曉平

?

埋管流化床內(nèi)濕顆粒流動及混合特性的CFD-DEM數(shù)值模擬

劉道銀，宋誠驍，王錚，馬吉亮，陳曉平

（東南大學能源與環(huán)境學院，能源熱轉換及其過程測控教育部重點實驗室，江蘇南京 210096）

摻液滴的氣固流化床，越來越多地應用于顆粒材料合成和改性工程，但是液滴使顆粒具有一定的濕度。濕顆粒具有黏性，與干燥顆粒流化特性有很大差異。本文應用遲滯碰撞模型，分別在顆粒碰撞的加載和卸載階段考慮液體黏彈阻力（lubrication force）和液橋力引起的動能耗損，建立濕顆粒流化床CFD-DEM模型。對含有埋管的濕顆粒流化床進行數(shù)值模擬研究，比較了不同黏性條件下顆粒的流動及混合特性。發(fā)現(xiàn)隨著顆粒黏性的增大，顆粒聚團加劇，氣泡邊界變得粗糙并逐漸轉變?yōu)闅怏w溝流，最終導致流化失效，埋管在床中起骨架作用，促進了溝流的形成。隨著顆粒黏性的增大，顆粒的相對運動被阻礙，顆粒速度降低，床層壓降逐漸減小，顆?；旌纤俣葴p慢。

濕顆粒流化；遲滯碰撞模型；顆粒聚團；顆?；旌?/p>

在眾多的應用氣固流化床的物理和化學工藝中，經(jīng)常伴有噴霧過程，例如催化劑顆粒制備、藥品顆粒制備、食品顆粒加工、廢液焚燒、烯烴聚合反應等。在此類摻液滴氣固流化床中，液滴使顆粒具有一定的濕度，濕顆粒具有黏性，可能導致顆粒團聚、流化不良甚至失流化現(xiàn)象。

濕顆粒流化特性已經(jīng)引起了較多的關注。在實驗研究方面，ZHOU等[1]結合圖像法和聲頻測量技術研究了不同含液量氣固流化床的流體力學特性。XU等[2]利應用圖像法研究了濕顆粒條件下噴動流化床的流型、最小噴動速度、最大噴射高度等。WANG等[3]運用PIV技術研究了鼓泡流化床中濕顆粒流化形態(tài)，并與模型預測結果進行對比。MA 等[4]研究了黏性顆粒流化床中氣泡特性（氣泡形狀、數(shù)量、直徑等），研究氣泡向氣體溝流轉變的規(guī)律。

在數(shù)值模擬研究方面，應用較多的是在DEM模型（discrete element method）[5]基礎上考慮顆粒黏性，主要有三類思路：①黏性力與接觸力線性疊加；②黏性力歸結到碰撞恢復系數(shù)；③應用遲滯碰撞模型考慮黏性作用。

第一類思路將液橋力和碰撞力疊加。MIKAMI等[6]歸納出靜態(tài)液橋力的回歸形式，在彈簧阻尼碰撞模型的基礎上，沿顆粒間法向疊加靜態(tài)液橋力，模型可以預測聚團流化現(xiàn)象。在MIKAMI濕顆粒碰撞模型基礎上，許多學者對不同對象濕顆粒體系，如滾筒、混合器、料倉、流化床等，進行了數(shù)值模擬研究。LIM等[7]研究了不同流化風速下濕顆粒流化床內(nèi)的顆?；旌咸匦浴ANG等[8-9]模擬流化床提升管中干濕顆粒的流化形態(tài)、聚團特性。HE等[10]對流化床中干濕顆粒的流化特性進行了研究。此類模型不能考慮液體黏度對顆粒碰撞黏結特性，因為靜態(tài)液橋力由液體表面張力系數(shù)和液體含量決定，與液體黏性無關，動態(tài)液橋力雖然與液體黏性關系較大，但是在流化床條件下通?？梢院雎?。實驗和理論研究均表明，顆粒相互靠近，排擠顆粒間的液體時產(chǎn)生黏彈阻力（lubrication force），對于高黏性液體而言尤為如此，因此直接忽略液體黏性的影響存在一定的局限性[11]。

第二類思路將濕顆粒碰撞動能耗損歸結到碰撞恢復系數(shù)，然后用濕恢復系數(shù)（wet restitution coefficient）替代彈簧阻尼碰撞模型中的碰撞恢復系數(shù)[12]。VAN BUIJTENEN等[13]研究了濕恢復系數(shù)對噴霧流化床造粒過程中顆粒相互作用的影響。DARABI等[14]提出了一種濕恢復系數(shù)模型并對流化床進行數(shù)值模擬研究。該模型的優(yōu)點在于同時考慮了液體黏性和液體表面張力對碰撞的影響，但是濕恢復系數(shù)很小時，彈簧阻尼碰撞模型中阻尼力過大，導致在加載階段接觸力可能減小的現(xiàn)象，偏離了實際物理過程；而當恢復系數(shù)接近0時，模型的求解存在奇點，導致該模型無法復現(xiàn)顆粒聚團 現(xiàn)象。

第三類思路運用遲滯碰撞模型[15]分階段考慮黏彈阻力和液橋力對碰撞的影響。DAVIS等[16]研究了黏彈阻力對顆粒碰撞的影響，并提出了基于Stokes量綱數(shù)為1的碰撞恢復系數(shù)表達式。MA等[17]研究了濕顆粒的碰撞過程，區(qū)分了顆?？拷A段的黏彈阻力耗損和分離階段的液橋力耗損。應用遲滯碰撞模型描述濕顆粒碰撞過程，加載過程考慮黏彈阻力耗損，卸載過程考慮液橋力耗損，既可以預測碰撞恢復系數(shù)為0時的情況（即顆粒發(fā)生黏結），同時又能考慮液體黏性和表面張力對動能耗損的影響。

常規(guī)顆粒在含埋管流化床中的運動特性已有很多研究[18-21]，但是針對濕顆粒的含埋管流化床的研究目前較為少見。本文作者將遲滯碰撞模型與CFD-DEM模型耦合，對含埋管濕顆粒流化床進行數(shù)值模擬研究，研究不同黏性和流化風速條件下的濕顆粒流動及混合特性（流動形態(tài)、顆粒速度、床層壓降、混合行為）。限于計算量，本文作者針對含約20萬個顆粒的二維流化床開展CFD-DEM數(shù)值模擬研究，主要目的在于從定性的角度考察分析濕顆粒流化特性以及模型適用性，同時為將來開展?jié)耦w粒三維流化床CFD-DEM定量研究和分析奠定基礎。

1 模型介紹

1.1 數(shù)學模型

在CFD-DEM模型[22-24]中，氣相用Navier-Stokes方程描述，顆粒相用牛頓第二定律跟蹤體系內(nèi)每個顆粒的受力和運動。對于每個顆粒，考慮顆粒自身重力、周圍顆粒對它的接觸力以及氣體對顆粒的曳力和壓力梯度力。氣固兩相之間采用雙向耦合，曳力應用Eugun和Wen&Yu關聯(lián)式。CFD-DEM模型的接觸力通常采用線性彈簧阻尼碰撞模型描述，而本文作者采用的接觸力模型為遲滯碰撞模型，下面詳細介紹。

1.2 濕顆粒碰撞模型

結合一對濕顆粒碰撞過程的接觸力與形變量關系曲線（圖1），對模型介紹如下[25]。圖1所示中軸表示接觸力，大于0表示顆粒間作用力為排斥力，小于0為吸引力。整個碰撞過程分為靠近、加載、卸載和分離。

階段1：加載過程從圖1中的點開始。當顆粒接觸并開始加載時，顆粒持續(xù)減速直到速度為0（點）。段的斜率即加載彈性模量（kload）。BC段顆粒的受力情況見式(1)。

F= kload+p(1)

階段2：在顆粒速度減小至0以后，顆粒在彈性力作用下開始反彈，卸載過程開始。卸載過程進行至點時，塑性形變量為0。段的斜率即卸載彈性模量（kunload）。段顆粒的受力情況見 式(2)。

F= kunload(–0)p(2)

k,unload與k,load的關系用黏彈恢復系數(shù)（wet）表示，wet=(k,load/k,unload)1/2。k,unload和wet是輸 入?yún)?shù)。wet由數(shù)（）的函數(shù)表示[26]，本文取用常見形式wet=exp(–35/)。=pp/(6π2)中，p為顆粒質量；為液體黏度；為顆粒 半徑[27]。

階段3：卸載過程沿線繼續(xù)進行，直到點。點以后彈性力消失，只有液橋力繼續(xù)作用于顆粒之間。液橋力可以由液體表面張力系數(shù)和液體含量函數(shù)確定[7]。在本文中，液橋力的最大值由數(shù)（）給定。數(shù)定義為液橋力與顆粒重力的比值，如式(3)。

式中，為液體表面張力系數(shù)；p為顆粒 密度。

階段4：分離階段從點到點，在此過程中，顆粒沒有直接接觸，但是液橋力仍然存在，直到顆粒間距離大于極限液橋距離c（點）[7]。

除了上述兩種由液體產(chǎn)生的能量耗損外，部分能量耗損以波等形式，耗散于空氣，采用干碰撞恢復系數(shù)（dry）表示。將顆粒發(fā)生碰撞黏結的臨界點記為臨界碰撞速度（p,cr）。p,cr越大，表示顆粒黏性越大，則越容易發(fā)生黏結。

2 結果和討論

2.1 模擬條件及工況介紹

模擬對象為準二維流化床，床體尺寸、顆粒物性等參數(shù)按照MA等[28]的實驗條件設置。床體尺寸30cm×90cm，床體深度為1倍顆粒直徑，顆粒數(shù)236148個，參數(shù)詳見表1。流化床含18根埋管，埋管半徑1cm，埋管布置高度分別為10cm、16cm、22cm、28cm和34cm，具體位置如圖2所示。模擬物理時間12s，為了避免顆粒初始堆積狀態(tài)對流化的影響，3s至12s的數(shù)據(jù)用于統(tǒng)計和分析。模型輸入?yún)?shù)包括量綱為1液橋力數(shù)、液體黏度和彈性模量kunload。數(shù)由液體表面張力和含液量確定。在本文模擬中，固定和kunload，改變數(shù)和流化風速g，共12個計算工況，見表2。

2.2 流動形態(tài)及壓力特征

圖2為g=2.5mf時3種黏性條件下的瞬時流動形態(tài)模擬結果。隨數(shù)增大，顆粒間黏性力作用越來越明顯，顆粒團聚尺寸逐漸增大，氣泡由較規(guī)則的圓形和帽形逐漸發(fā)展為具有分叉、邊界粗糙的不規(guī)則形狀，高黏性時氣泡消失，形成穩(wěn)定的氣體溝流。MIKAMI[7]、HE[10]等的模擬結果也反映類似的氣泡形態(tài)特性。此外，由顆粒顏色看到，黏性增加導致顆粒整體運動速度降低。氣體溝流的出現(xiàn)和顆粒運動速度的降低，均表明流化床內(nèi)氣固流動受到抑制。埋管起到了骨架的作用，使得高黏性條件下顆粒團聚可能阻塞于埋管間、埋管與壁面間，從而促進了流化失效。在MA等[28]的實驗中，黏性力的引入由多聚物涂層法實現(xiàn)，黏性力的大小與床溫呈近似線性關系，無多聚物涂層即無黏性，有多聚物涂層、床溫36℃即低黏性，有多聚物涂層、床溫44℃即高黏性。將本文的模擬結果與MA等[28]的實驗結果（圖3）定性對照，黏性顆粒引起的聚團、氣體溝流等特征基本與實驗相符。

表1 數(shù)值模擬物理參數(shù)表

表2 數(shù)值模擬工況表

圖4為g=2.5mf時3種顆粒黏性條件床層壓降動態(tài)曲線。按照床料重量計算的床層壓降理論值為3635Pa，模型計算得到的各顆粒黏性條件的床層壓降與理論值接近，說明了模型合理性。圖5為時均床層壓降圖，其中圖5(a)對比不同黏性條件，圖5(b)對比不同流化風速條件。隨著顆粒黏性的逐漸增大，床層壓降逐漸減小，并且顆粒黏性導致床層壓降減小這一現(xiàn)象，在不同流化風速條件下都可以觀察到。這是由于黏性增大導致顆粒聚團，并引起流化不良和氣體溝流現(xiàn)象，導致床層壓降下降。

圖3 Ug=2.5Umf時MA[28]等實驗拍攝不同顆粒黏性的瞬時流動形態(tài)圖

2.3 顆粒及氣體速度

圖6為g=2.5mf時3種顆粒黏性條件下顆粒時均速度（p）矢量圖，直觀反映顆粒在床內(nèi)的循環(huán)情況。干顆粒的循環(huán)狀態(tài)清晰明顯，隨著黏性的增加，顆粒速度降低，顆粒循環(huán)逐漸變?nèi)酰?56時接近“死床”。由于埋管的存在，顆粒在床內(nèi)大循環(huán)的同時有繞埋管小循環(huán)趨勢，低黏性時較明顯。

圖7為g=2.5mf時3種顆粒黏性條件下時均氣體速度（gl）矢量圖。由于氣泡集中在密相區(qū)床層中部，干顆粒條件下表現(xiàn)為床層中部gl較大，壁面附近較小。隨著黏性增大，床層中部的時均氣速減小，表明氣固流動受到一定抑制，高黏性條件下溝流穩(wěn)定存在。

圖8為床層20cm高度處，時均顆粒水平分速度（px）和軸向分速度（py）的分布情況，其中圖8(a)和8(b)對比不同顆粒黏性情況。隨數(shù)增大，環(huán)核流動趨勢越來越弱，在低黏性條件下（=6）環(huán)核流動尚可分辨，隨著黏性加強（=17以上工況）環(huán)核流動基本消失。此外，由于采樣點位置接近高度22cm處的埋管，可觀察到該高度的4處埋管（=6cm，12cm，18cm，24cm）對顆粒整體環(huán)核流動的分解作用，其中=12cm和=18cm的埋管對顆粒速度方向改變效果最為明顯。圖8(c)和8(d)對比不同流化風速情況，可以看到提高風速可一定程度改善黏性顆粒流化，促進床內(nèi)顆粒循環(huán)運動。

圖6g=2.5mf時不同顆粒黏性條件下時均顆粒速度矢量分布（為了便于觀察，每4個網(wǎng)格顯示顆粒速度）

圖7 為g=2.5mf時不同顆粒黏性條件下時均氣體速度矢量分布（為了便于觀察，每3個網(wǎng)格顯示顆粒速度）

2.4 顆?；旌?/p>

圖9為g=2.5mf時不同黏性條件的顆?；旌蠒r序圖，選取某一時刻為起始時刻0，時序間隔為2s，其中藍色顆粒為示蹤顆粒。為研究顆粒在水平方向的混合情況，將流化床從左到右等分10個區(qū)域，依次命名為L1～L10，示蹤顆粒初始布置在L4和L5區(qū)域。同理，為研究軸向混合情況，將床體從下至上每3cm共劃分10個區(qū)域，依次命名為V1～V10，示蹤顆粒初始布置在V6和V7 區(qū)域。

圖10為示蹤顆粒在床內(nèi)各區(qū)域的濃度分布曲線圖，(a)、(b)為干顆粒，(c)、(d)為=17的濕顆粒。結合混合時序圖與示蹤顆粒濃度分布圖可以發(fā)現(xiàn)，干顆粒在兩個方向的混合明顯快于濕顆粒，且隨黏度增加，混合速率進一步降低。當顆粒黏度較高時，顆粒的混合進度幾乎停滯，表明流化嚴重不良。此外，干顆粒水平方向達到混合均勻，大約需要5.5s，而豎直方向達到混合均勻，需要約2.5s，所以顆粒水平混合慢于軸向混合，該規(guī)律在黏性顆粒流化時也能觀察到。這是由于流化床顆粒豎直方向運動比豎直方向運動更劇烈。

3 結論

應用CFD-DEM數(shù)值模擬耦合遲滯碰撞模型，研究埋管流化床中濕顆粒的流動及混合特性，主要結論如下。

（1）顆粒間的黏性力阻礙顆粒的相對運動，隨著黏性力增大，更多顆粒形成聚團，并最終導致氣泡轉變?yōu)闅怏w溝流，與實驗結果對照定性合理。埋管起骨架作用，加速流化失效。時均床層壓降隨著顆粒黏性增大而減小。

（2）干顆粒床層中部的時均氣體速度較壁面處大，高黏性條件下時均氣速分布可反映穩(wěn)定溝流的產(chǎn)生。埋管導致顆粒運動發(fā)生較大幅度的轉向，并將顆粒環(huán)核流動分解為若干小環(huán)。顆粒速度隨著顆粒黏性增大逐漸降低，環(huán)核流動逐漸消失，床內(nèi)顆粒循環(huán)被破壞。

（3）顆粒混合速度隨顆粒黏性增大而減小，高黏性時混合難以均勻。顆粒沿水平混合的速率通常大于豎直方向。

[1] ZHOU Y F，SHI Q，HUANG Z L，et al. Effects of liquid action mechanisms on hydrodynamics in liquid-containing gas-solid fluidized bed reactor[J]. Chemical Engineering Journal，2016，285：121-127.

[2] XU H B，ZHONG W Q，YU A B，et al. Spouting characteristics of wet particles in a conical-cylindrical spouted bed[J]. Industrial & Engineering Chemistry Research，2015，54（40）：9894-9902.

[3] WANG T Y，HE Y R，TANG T Q，et al. Experimental and numerical study on a bubbling fluidized bed with wet particles[J]. AIChE Journal，2016，62：1970-1985.

[4] MA J L，LIU D Y，CHEN X P. Bubble behaviors of large cohesive particles in a 2D fluidized bed[J]. Industrial & Engineering Chemistry Research，2016，55（3）：624-634.

[5] TSUJI Y，KAWAGUCHI T，TANAKA T. Discrete particles simulation of 2-dimensional fluidized bed[J]. Powder Technology，1993，77（1）：79-87.

[6] MIKAMI T，KAMIYA H，HORIO M. Numerical simulation of cohesive powder behavior in a fluidized bed[J]. Chemical Engineering Science，1998，53（10）：1927-1940.

[7] LIM C，WEE E，BENG R，et al. Mixing behaviors of wet granular materials in gas fluidized bed systems[J]. AIChE Journal，2013，59（11）：4058-4067.

[8] WANG M，ZHU W B，SUN Q Q，et al. A DEM simulation of dry and wet particle flow behaviors in riser[J]. Powder Technology，2014，267：221-233.

[9] 王猛. 循環(huán)流化床提升管內(nèi)氣固兩相流的DEM數(shù)值模擬[D]. 哈爾濱：哈爾濱工程大學，2013.

WANG M. A DEM numerical simulation of gas-solid two-phase flow in circulating fluidized bed riser[D]. Harbin：Harbin Engineering University，2013.

[10] HE Y R，PENG W G，WANG T Y，et al. DEM study of wet cohesive particles in the presence of liquid bridges in a gas fluidized bed[J]. Mathematical Problems in Engineering，2014，1：1-14.

[11] SEVILLE J P K，CLIFT R. The effect of thin liquid layers on fluidisation characteristics[J]. Powder Technology，1984，37（1）：117-129.

[12] SUTKAR V S，DEEN N G，PATIL A V，et al. CFD-DEM model for coupled heat and mass transfer in a spout fluidized bed with liquid injection[J]. Chemical Engineering Journal，2016，288：185-197.

[13] VAN BUIJTENEN M S，DEEN N G，HEINRICH S，et al. A discrete element study of wet particle-particle interaction during granulation in a spout fluidized bed[J]. Canadian Journal of Chemical Engineering，2009，87（2）：308-317.

[14] DARABI P，POUGATCH K，SALCUDEAN M，et al. DEM investigations of fluidized beds in the presence of liquid coating[J]. Powder Technology，2011，214（3）：365-374.

[15] LUDING S. Cohesive，frictional powders: contact models for tension[J]Granular Matter，2008，10（4）：235-246.

[16] DAVIS R H，RAGER D A，GOOD B T. Elastohydrodynamic rebound of spheres from coated surfaces[J]. Journal of Fluid Mechanics，2002，468：107-119.

[17] MA J L，LIU D Y，CHEN X P. Experimental study of oblique impact between dry spheres and liquid layers[J]. Phys. Rev. E. Stat. Nolin Soft Matter Phys.，2013，88（3）：033018.

[18] YATES J G，RUIZ-MARTINEZ R S，CHEESMAN D J. Prediction of bubble size in a fluidized bed containing horizontal tubes[J]. Chemical Engineering Science，1990，45（4）：1105-1111.

[19] HULL A S，CHEN Z，F(xiàn)RITZ J W，et al. Influence of horizontal tube banks on the behavior of bubbling fluidized beds：1. Bubble hydrodynamics[J]. Powder Technology，1999，103（3）：230-242.

[20] ASEGEHEGN T W，SCHREIBER M，KRAUTZ H J. Investigation of bubble behavior in fluidized beds with and without immersed horizontal tubes using a digital image analysis technique[J]. Powder Technology，2011，210（3）：248-260.

[21] ASEGEHEGN T W，SCHREIBER M，KRAUTZ H J. Numerical simulation and experimental validation of bubble behavior in 2D gas–solid fluidized beds with immersed horizontal tubes[J]. Chemical Engineering Science，2011，66（21）：5410-5427.

[22] VAN WACHEM B G M，VAN DER SCHAAF J，SCHOUTEN J C，et al. Experimental validation of Lagrangian–Eulerian simulations of fluidized beds[J]. Powder Technology，2001，116（2）：155-165.

[23] LIU D Y，BU C S，CHEN X P. Development and test of CFD-DEM model for complex geometry：a coupling algorithm for fluent and DEM[J]. Computers & Chemical Engineering，2013，58：260-268.

[24] LUO K，YANG S L，F(xiàn)ANG M M，et al. LES-DEM investigation of the solid transportation mechanism in a 3-D bubbling fluidized bed. part I: Hydrodynamics[J]. Powder Technology，2014，256：385-394.

[25] LIU D Y，VAN WACHEM B G M，MUDDE R F，et al. An adhesive CFD-DEM model for simulating nanoparticle agglomerate fluidization[J]. AIChE Journal，2016，62（7）：2259-2270.

[26] LEGENDRE D，ZENIT R，DANIEL C，et al. A note on the modelling of the bouncing of spherical drops or solid spheres on a wall in viscous fluid[J]. Chemical Engineering Science，2006，61（11）：3543-3549.

[27] DAVIS R H，SERAVSSOL J M，HINCH E J. The elastohydrodynamic collision of 2 spheres[J]. Journal of Fluid Mechanics，1986，163：479-497.

[28] MA J L，LIU D Y，CHEN X P. Bubbling behaviors of cohesive particles in a two-dimensional fluidized bed with immersed tubes[J]. Particuology，2017，31：152-160. http:// dx.doi.org/10.1016/ j.partic. 20106.005.

CFD-DEM simulation of wet particles flow and mixing behavior in fluidized bed with immersed tubes

LIU Daoyin，SONG Chengxiao，WANG Zheng，MA Jiliang，CHEN Xiaoping

（Key Laboratory of Energy Thermal Conversion and Control of Ministry of Education，School of Energy and Environment，Southeast University，Nanjing 210096，Jiangsu，China）

Gas-solid fluidized beds with the presence of some liquid droplets are applied widely in industry. Due to the existence of droplets，particles become wet and viscous，which makes fluidization significantly different with dry particles. In this paper，a hysteresis contact model is applied to consider energy dissipation of lubrication force and liquid bridge force during wet particle collision，which is coupled with CFD-DEM model. The fluidization and mixing characteristics of wet particles under different conditions in a 2D fluidized bed with immersed tubes are compared. Generally，it is found that with the presence of liquid，the particles begin to agglomerate and the bubbles become gas channels. The immersed tubes play the part of frameworks，promoting the formation of gas channels. With increasing cohesiveness，the particle velocity decreases generally and the internal recirculation of wet particles is restrained. It is also found that with an increase in the cohesiveness，the pressure drop of the fluidized bed and the mixing speed of particles decrease.

wet particle；hysteresis contact model；agglomerate；solid mixing

TQ021

A

1000–6613（2017）06–2070–08

10.16085/j.issn.1000-6613.2017.06.015

2016-09-20；

2016-12-08。

國家自然科學基金項目（51306035，51676042）。

劉道銀（1982—），男，副研究員，研究方向為多相流及數(shù)值模擬。E-mail：dyliu@seu.edu.cn。