基于DEM模擬液固流化床粒度級配對顆粒流動特性的影響

丁冬峰，程 可，陸曉峰，朱曉磊，朱凌雪

（1.南京工業大學機械與動力工程學院，江蘇南京 211816；2.金陵科技學院，江蘇南京 211816）

液固流化床具有混合均勻、傳熱和傳質性能好等優點，目前已被廣泛應用于化工、食品技術、濕法冶金和水處理等諸多領域[1-2]。目前對于氣固流化床已經有了較為深入和系統的研究，而對液固流化床的研究較少。液固流化床的放大、設計和運行主要取決于對顆粒力學行為和流動特征（如顆粒分布、流動形態和顆粒混合狀態）的準確預測。但是由于流化床內顆粒流體系統具有結構非均勻性、狀態多值性、結構突變和結構的多尺度特征[3]，迄今為止仍然沒有一套完整的理論體系來描述流化床內顆粒流體系統的動態特性，更缺乏液固流化床的放大與設計的通用方法，因此國際上此領域的研究受到學者的廣泛關注[4-5]。

在流化床的研究中，實驗方法可以獲得較為可信的結果。但是，受限于實驗條件和測試水平，實驗測試不但成本高昂，而且很難獲得顆粒在微觀尺度上的運動信息。隨著計算機技術的快速發展，數值模擬方法成為研究流化床系統的一個重要手段[6]。長期以來，液固兩相流的模擬主要采用的是Euler-Euler雙流體模型[7-8]。該方法將顆粒相處理為連續的擬流體，但這一假設從本質上削弱了模擬液固兩相系統中非均勻結構的真實性。而Euler-Lagrange法可以很好地考慮液固兩相系統中非均勻結構特性，并且可以獲得豐富的顆粒微觀信息，目前已被廣泛應用于流化床的模擬領域[9-11]。

在實際工業應用中，單一粒度的顆粒系統是沒有的，幾乎都是由多種粒度組成的混合顆粒體系。顆粒粒徑分布（PSD）對顆粒的流化行為有著很重要的影響。孫國剛等[12]的實驗表明，對于存在粒度級配的顆粒群如果僅僅用單一的平均粒徑，很難全面描述其內部復雜的流體動力學行為，應該要考慮粒徑分布的影響。國內外很多學者對多組分顆粒體系進行了研究[13-15]，但是這些研究成果還是主要集中在氣固流化床領域。鑒于此，以液固流化床為研究對象，采用CFD-DEM模擬方法，對二元混合、窄級配與寬級配3種粒度級配方式進行模擬，分別研究粒度級配對顆粒分布、顆粒速度、顆粒溫度波動以及顆粒混合度的影響。

1 數學模型

1.1 連續相控制方程

連續性方程為

動量方程為

式中：ρl為液相密度；εl為液相體積分數；ul為液相速度；p為壓力；τl為液相剪切應力張量，其定義為

式中：μl是液相的剪切黏度；g是重力加速度。Fint為顆粒對流體的作用力，其計算公式為

式中：βint為固相與液相間的曳力系數，βint（ul-up）、Fvm分別為曳力、虛擬質量力。

1.2 顆粒相控制方程

假設顆粒形狀為球體，顆粒相的控制方程可用牛頓第二定律進行描述。每個顆粒都有平動和轉動2種運動類型，其計算公式為

式中：mp和up是顆粒的質量和速度；式（5）右邊各項分別為顆粒的重力、接觸碰撞力、潤滑力、相間力和浮力；Ip是顆粒的轉動慣量；ω是顆粒的角速度；Tp是由顆粒切向接觸力引起的轉矩。

1.2.1 接觸碰撞力

顆粒接觸模型如圖1所示。顆粒與顆粒、顆粒與壁面之間的接觸力采用Hertz-Mindlin軟球接觸模型。接觸力由法向力Fcn和切向力Fct組成，計算公式為

式中：k和η分別為彈性系數和阻尼系數；δ為重疊量；v為相對速度；μs為靜摩擦系數。

1.2.2 潤滑力

考慮到液固兩相流顆粒流動的特殊性，故在經典的Hertz-Mindlin接觸模型的基礎上，考慮增加一個潤滑力來表征顆粒之間填隙流體對顆粒碰撞的影響，如圖1所示。當顆粒相互遠離時，液體將流入顆粒間的縫隙內；當顆粒相互靠近時，液體將從縫隙內流出。流體的運動將產生一個壓力梯度和黏性應力，它們累積后作用于顆粒上。此力將抑制顆粒間的分離。當2個顆粒相互靠近時，潤滑力演變為排斥力；當相互遠離時，潤滑力演變為吸引力。潤滑力[16]的計算公式為

式中：v為2個顆粒的相對滑動速度；dp，i為i顆粒粒徑。從潤滑力的表達式中可以看出，當2個顆粒距離無限近的時候，潤滑力會阻止2個顆粒的接觸，并且會發散到無窮大。但是實際情況下，顆粒上還有其他作用力，并且顆粒表面也是粗糙的，所以在潤滑力達到理論的極值之前，顆粒就會發生接觸。為了模擬真實情況，所以設置截止距離δn＞0.001dp，來防止潤滑力計算發散。

圖1 顆粒接觸模型Fig.1 Model of contact force between two particles

1.2.3 浮力

在液固兩相流中，由于液體的密度相較于氣體較大，所以液體對顆粒的浮力更加重要，浮力的計算公式為

1.3 相間耦合方程

顆粒流體系統的相間作用力很復雜，包括曳力、虛擬質量力、Magnus力、Saffman力、Basset力等。由于液體有較高的密度和黏度，因此本文中的相間耦合模型中考慮了曳力和虛擬質量力，而且在流化床中其他力與這2種力相比都較小，可以忽略。

1.3.1 虛擬質量力

在液固兩相流中，顆粒相對于流體加速運動時，必將帶動其周圍的部分流體產生擾動，使其加速。對于球形顆粒，虛擬質量力的計算公式為

1.3.2 曳力

采用Lu等[17]所提出的Huilin-Gidaspow曳力模型，該模型避免了Wen&Yu和Ergun模型在孔隙率為0.8處的不連續性，通過引入函數讓曳力模型變成一條光滑的曲線。Huilin-Gidaspow曳力模型中的曳力系數計算公式為

單顆粒曳力系數Cd和Reynolds數Re的一般關系為

1.4 數值模擬參數及邊界條件

三維流化床幾何模型如圖2所示，模型直徑為140 mm，高度為1 000 mm。顆粒平均粒徑和密度分別為3 mm和2 500 kg/m3。初始條件下，顆粒在床層內堆積形成一定高度，床層顆粒總質量為9.7 kg。采用速度入口和壓力出口邊界條件，入口處液相速度假設是均勻分布，液相在壁面處采用無滑移邊界條件。模擬參數如表1所示。采用Gambit進行網格劃分，網格類型選擇Hex或Wedge，網格尺寸為10 mm。采用MFIX-DEM開源代碼進行模擬。

圖2 三維流化床幾何模型Fig.2 3D fluidized bed geometry model

表1 模擬參數Tab.1 Parameters used in simulations

2 模型驗證

Limtrakul等[18]采用計算機輔助放射性示蹤粒子成像技術（CARPT）和射線計算機輔助層析成像技術（CT）測量了液固流化床的軸向速度和固體體積分數，實驗所用的床體直徑為1 40mm，高度為1 500 mm，顆粒粒徑為3 mm，密度為2 500 kg/m3，液體表觀流速為0.07 m/s。為了驗證模型的有效性，將模擬結果與Limtrakul等的實驗結果和Wang等[16]采用二維模型的模擬結果進行對比。

圖3和圖4分別為顆粒軸向速度和顆粒濃度的徑向分布。從圖中可以看出，模擬結果與Limtraku的實驗結果吻合較好。而Wang模擬的是二維流化床，顆粒被簡化為一個圓盤，缺少了一個維度上的信息，對顆粒運動產生了明顯影響，所以模擬結果誤差很大。

圖3 顆粒軸向速度的徑向分布Fig.3 Radial distribution of axial velocity of particles

3 粒度級配顆粒流動特性的模擬

3.1 顆粒粒徑分布

為了研究液固流化床內粒度級配顆粒的流動特性，設計了3種不同粒徑分布的顆粒組合，包括二元混合顆粒、正態分布顆粒（窄級配）以及均勻分布顆粒（寬級配），3種分布類的顆粒平均當量直徑均為3 mm。模擬所用顆粒的粒徑分布如表2所示。

表2 顆粒粒徑分布Tab.2 Particle size distribution

3.2 模擬結果

3.2.1 顆粒分布分析

為了研究3種不同的粒度級配方式對液固流化床內顆粒移動及分布情況的影響，將A、B、C、D、E 5種顆粒分別用不同的顏色表示。在液體流速為0.1 m/s條件下，3種不同粒度級配在0～6 s內的顆粒分布如圖5所示。

圖4 顆粒濃度的徑向分布Fig.4 Radial distribution of particle concentration

圖5 粒度級配對顆粒分布的影響Fig.5 Influence of particle size distribution on particle distribution

從圖可以看出，在t=0 s時刻，顆粒在堆積形成的床層內均勻混合。隨著時間的推移，顆粒在流體的作用下，床層均勻向上膨脹，不同尺寸的顆粒逐漸分層。床層底部有小液泡產生，并逐漸向上運動，運動的過程中液泡不斷融合、破裂和重新產生。在4 s以后，床層高度基本保持不變，由此可以認為流化床在4 s以后達到穩定狀態。

從圖5中還可以看出，二元混合和寬級配流化床的上部區域均有黑色的小顆粒聚集，這是由于這2種粒度級配方式中小顆粒的比例比較大。顆粒所受的曳力與粒徑的平方成正比，顆粒的重力卻與粒徑的3次方成正比。當顆粒粒徑減小的時候，重力減小的速度遠大于曳力減小的速度，小顆粒的曳力會大于重力，所以小顆粒會更加容易的往上部運動。該現象表明二元混合和寬級配這2種粒度級配方式更容易產生顆粒的偏析。

為了定量研究流化床內顆粒的分層情況，將流化床分成3個區域，每個區域高度為0.2 m，研究不同高度處的大顆粒與小顆粒分布情況。3種粒度級配方式中大顆粒與小顆粒在不同高度的分布如圖6—8所示。

從圖6a、7a、8a中可以看出，不同區域的小顆粒質量變化一直處于波動過程，且在床層上部區域小顆粒質量都呈現出明顯的增大趨勢。這說明由于小顆粒具有較好的流動特性，會在整個床層中進行循環流動，即由上部區域進入下部區域，再由下部區域流動到上部區域。但是小顆粒的總體運動趨勢是向床層的上部聚集。相反，從圖6b、7b、8b中可以看出，床層下部區域大顆粒質量呈現出明顯的增大趨勢，且增大趨勢比較平緩。這表明大顆粒的流動特性較差，從床層上部一直往下部運動。圖7中還反應出窄級配中大顆粒與小顆粒在不同高度處的質量波動要比二元混合和寬級配更加劇烈。

圖6 二元混合中大顆粒與小顆粒在不同高度的分布Fig.6 Distribution of large particles and small particles at different heights in binary mixture

圖7 窄級配中大顆粒與小顆粒在不同高度的分布Fig.7 Distribution of large particles and small particles at different heights in narrow size distribution

圖8 寬級配中大顆粒與小顆粒在不同高度的分布Fig.8 Distribution of large particles and small particles at different heights in wide size distribution

3.2.2 顆粒速度分布分析

粒度級配對顆粒軸向速度的影響如圖9所示。從圖可以看出，所有的曲線均表現出相同的趨勢：沿著徑向方向顆粒軸向速度呈現出中心高，邊壁低的不均勻分布，說明顆粒流動形式為典型的環核流動結構形式。產生環核流動的主要原因是液體徑向速度分布不均勻，導致顆粒在不同位置所受到的曳力不同。在中心位置液體的流速大，對顆粒的曳力也大，所以顆粒速度較高；在壁面位置，液體速度降低，并且顆粒還受到邊壁的阻滯作用的影響，使得顆粒速度出現負值。從圖9b中可以看出，窄級配顆粒的軸向速度在不同高度差異不大，這與Roy等[19]對單一粒徑顆粒的實驗結果相似，這說明窄級配表現出單分散顆粒的運動特性。而從圖9a和9c中可以看出，二元混合和寬級配的軸向速度在不同高度上差異較明顯，這是由于粒徑分布較寬，導致在不同高度上顆粒分布不均勻，從而使得顆粒軸向速度在不同高度上出現差異。

圖9 粒度級配對顆粒軸向速度的影響Fig.9 Influence of particle size distribution on particle axial velocity

3.2.3 顆粒溫度波動分析

顆粒溫度波動是基于顆粒隨機波動速度的統計函數，顆粒溫度波動與氣體動力學理論中的溫度類似，是度量顆粒碰撞強度的函數，同時也是測量顆粒系統的動力學行為最重要的參數之一。在三維系統中，可以定義為

式中：uxi、uyi和 uzi分別為 i顆粒 x、y、z 3 個方向的隨機波動速度，是瞬時速度和平均速度的差值。N為顆粒個數。

圖10為粒度級配對顆粒溫度波動的影響。從圖中可以看出，顆粒溫度波動主要集中在顆粒體積分數為0.3～0.5這個區間，而在這個區間的兩側顆粒溫度波動均較低。這是由于當顆粒體積分數較低的時候，顆粒個數較少，所以顆粒間發生碰撞的概率也很小；而顆粒在較高的體積分數下，由于顆粒之間接觸十分緊密，阻礙了顆粒的相互運動，從而使顆粒的脈動程度降低，顆粒溫度波動也減小。從圖中還可以看出，3種粒度級配的溫度波動關系：窄級配＞寬級配＞二元混合。這是由于窄級配中的大顆粒比例較少，能表現出很好的流化性能，導致顆粒的脈動程度也較高，從而使窄級配有較大的顆粒溫度波動。由于窄級配中顆粒有較大的溫度波動，因此在不同高度處，窄級配顆粒質量波動會比二元混合和寬級配更加劇烈。

圖10 粒度級配對顆粒溫度波動的影響Fig.10 Influence of particle size distribution on particle temperature

3.2.4 顆粒混合度分析

為了定量分析不同種顆粒的混合程度，引入Lacey[20]混合指數M對3種粒徑分布下的顆粒混合與分離情況進行分析。M表達式為

式中：S2為2種顆粒的實際混合方差；為完全分離時的方差，S02＝q（1-q），q為其中一種顆粒的體積分數；SR2為完全混合時的方差，SR2＝ q（1-q）/N，N 為一個樣本內的平均顆粒數目。

圖11為粒度級配對顆粒混合指數的影響。圖11a為寬級配的混合指數變化，從圖中可以看出，A顆粒和E顆粒、B顆粒和E顆粒之間的混合指數下降幅度較大，均降低到了0.9以下；A-C、B-D、C-E和A-D之間的混合指數也減小到了0.9～0.95之間；而A-B、B-C、D-E和C-D之間的混合指數從混合指數為1的均勻混合狀態降低到一定水平，不再進一步減小，只是略有波動，且混合指數均在0.95以上。

圖11 粒度級配對顆粒混合指數的影響Fig.11 Influence of particle size distribution on mixing index

這是由于小顆粒較容易往床層上部運動，而大顆粒則有向床層下部聚集的趨勢，故2種顆粒粒徑差距越大，越容易分離，而2種顆粒粒徑越接近，越難分離，混合程度越好。圖11b為窄級配的混合指數變化，從圖中可以看出，顆粒間的混合指數下降幅度均有所減小，A-E、B-E和A-D之間的混合指數與其他顆粒相比下降幅度還是較大，但是均在0.9以上。這說明，由于窄級配顆粒溫度波動較大，顆粒波動較劇烈，因此窄級配與寬級配相比，床層內顆粒整體上混合的較均勻，顆粒偏析情況不明顯。

4 結論

采用CFD-DEM耦合方法模擬了三維液固流化床系統，在模型中考慮了潤滑力對顆粒運動的影響。通過與Limtraku的實驗結果對比，發現本文中所采用的數學模型可以有效模擬液固流化床系統。隨后研究了二元混合、窄級配與寬級配3種粒度級配方式對顆粒流動特性的影響，得到如下結論：

1）大顆粒有向床層下部運動的趨勢，而小顆粒則有向床層上部聚集的趨勢，且小顆粒的流動特性要比大顆粒好。

2）與窄級配相比，二元混合和寬級配的軸向速度在不同高度上差異較明顯，且顆粒溫度波動也較小。

3）采用混合指數定量考察流化床內顆粒的混合狀態，發現2種顆粒粒徑差距越大，越容易分離，而粒徑越接近，越難分離，混合程度越好；窄級配與寬級配相比，床層內顆粒整體上混合的較均勻，顆粒偏析情況不明顯。

參考文獻（References）：

[1]WU C，XIAO M，CHEN S，et al.Application of fluidization technology in productions of cemented carbides[J].International Journal of Refractory Metals&Hard Materials，2014，46（1）:188-194.

[2]BELLO M M，RAMAN A A A，PURUSHOTHAMAN M.Applications of fluidized bed reactor in wastewater treatment A review of the major design and operational parameters[J].Journal of Cleaner Production，2016，141:1492-1514.

[3]李靜海，歐陽潔，高示秋，等.顆粒流體復雜系統的多尺度模擬[M].北京:科學出版社，2005:12-25.

[4]PAN H，CHEN X Z，LIANG X F，et al.CFD simulations of gas liquid solid flow in fluidized bed reactors-a review[J].Powder Technology，2016，299:235-258.

[5]LIU G，YU F，LU H，et al.CFD-DEM simulation of liquidsolid fluidized bed with dynamic restitution coefficient[J].Powder Technology，2016，304:186-197.

[6]任立波.稠密顆粒兩相流的CFD-DEM耦合并行算法及數值模擬[D].山東：山東大學，2015.

[7]PAN H，LIANG X，LUO Z.CFD modeling of the gas solid two-fluid flow in polyethylene FBRs:from traditional operation to super-condensed mode[J].Advanced Powder Technology，2016，27（4）:1494-1505.

[8]FENG X，LI X，CHENG J，et al.Numerical simulation of solid-liquid turbulent flow in a stirred tank with a two-phase explicit algebraic stress model[J].Chemical Engineering Science，2012，82:272-284.

[9]OZEL A，MOTTA J C B D J C，ABBAS M，et al.Particle resolved direct numerical simulation of a liquid-solid fluidized bed:comparison with experimental data[J].International Journal of Multiphase Flow，2017，89:228-240.

[10]朱衛兵，李錦時，王猛，等.提升管內濕顆粒流動特性的離散元模擬[J].中國電機工程學報，2013，33（35）:81-88.

[11]王猛，朱衛兵，孫巧群，等.提升管內氣固流動特性的離散元模擬[J].化工學報，2013，64（7）:2436-2445.

[12]孫國剛，李靜海.粒徑分布對循環床內顆粒速度分布的影響[J].中國粉體技術，1998，4（1）:1-5.

[13]趙永志，鄭津洋.寬粒徑分布流化床的微觀尺度模擬與分析[J].中國電機工程學報，2007，27（35）:55-61.

[14]GAUTHIER D，ZERGUERRAS S，FLAMANT G.Influence of the particle size distribution of powders on the velocities of minimum and complete fluidization[J].Chemical Engineering Journal，1999，74（3）:181-196.

[15]WANG S Y，LU H L，LI X，et al.Discrete particle simulations for flow of binary particle mixture in a bubbling fluidized bed with a transport energy weighted averaging scheme [J]. Chemical Engineering Science，2009，64（8）:1707-1718.

[16]WAMG S，GUO S，GAO J，et al.Simulation of flow behavior of liquid and particles in a liquid solid fluidized bed[J].Powder Technology，2012，224:365-373.

[17]LU H L，GIDASPOW D.Hydrodynamics of binary fluidization in a riser:CFD simulation using two granular temperatures[J].Chemical Engineering Science，2003，58（16）:3777-3792.

[18]LIMTRAKUL S，CHEN J，RAMACHANDRAN P A，et al.Solids motion and holdup profiles in liquid fluidized beds[J].Chemical Engineering Science，2005，60（7）:1889-1900.

[19]ROY S，SAI P S T，JAYANTI S.Numerical simulation of the hydrodynamics of a liquid solid circulating fluidized bed[J].Powder Technology，2014，251:61-70.

[20]LACEY P M C.Developments in the theory of particle mixing[J].Journal of Applied Chemistry，1954，4（5）:257-268.