基于DEM模擬氣固鼓泡床中顆粒碰撞參數對流場間歇性的影響

彭麗，吳迎亞，李佳瑤，高金森，藍興英

（中國石油大學（北京）重質油國家重點實驗室，北京 102249）

引 言

鼓泡床具有良好的傳熱和傳質特點，在化學工程、生物環境工程以及食品加工方面有著廣泛的應用[1]。計算流體力學方法（CFD）[2-3]常用于研究鼓泡床的流動行為[4]，其中基于歐拉-拉格朗日的CFD-DEM 方法常用于模擬追蹤顆粒相的運動，該方法可以從單顆粒層面上研究顆粒間的碰撞過程和氣體-顆粒間的相互作用過程，從而能夠更加準確地模擬出整個床層的氣固流動過程，所以DEM 方法是研究鼓泡床內微觀機制的重要方法。目前，DEM方法主要采用軟球模型處理顆粒間接觸、變形、受力和運動，將顆粒的碰撞視為彈簧-阻尼系統[5-8]，分別采用顆粒彈性系數（k）和顆?；謴拖禂担╡）描述顆粒碰撞作用力和碰撞過程中的能量損失。因此，在氣固鼓泡床流動的CFD-DEM 模擬中這兩個參數是影響模擬結果的關鍵參數。已有的實驗和模擬研究表明，鼓泡床內顆粒運動并非是完全無規則的運動，存在著大量的有序運動，這種有序運動被稱為相干結構，具體表現為顆粒渦團[9]，這種相干結構從微觀上表征顆粒速度脈動的生成與發展。Sun 等[10]和Jiradilok 等[11]指出正是由于相干結構的存在引起了流場間歇性，即顆粒渦團在時空中的不均勻出現或變化。因此，流場相干結構和間歇性均反映了鼓泡床內氣固流動的微觀特性。已有的研究主要是從宏觀角度[12-16]研究顆粒碰撞屬性對鼓泡床氣固流動的影響，很少從微觀角度進行分析。但鼓泡床內氣固運動的微觀特性對傳熱和傳質有著重要的影響，因此本研究重點考察了顆粒碰撞屬性對流場相干結構和間歇性等微觀行為的影響，為找到一個準確模擬氣固鼓泡流動行為的方法提供依據。有研究者在研究湍流相干結構時引入了小波分析方法，發現小波分析能在多個尺度上剖析湍流的相干結構[17]。前期本課題組已應用CFD-DEM 方法對氣固鼓泡床內氣固流動進行了模擬研究，并利用小波分析方法分析了床層內的多尺度相干結構[18]。本研究在前期工作基礎上分析CFD-DEM 方法中的關鍵參數——顆粒彈性系數和恢復系數對鼓泡床氣固流動行為的影響，并利用小波分析方法剖析顆粒彈性系數和恢復系數對鼓泡床內流場間歇性和相干結構的影響。

1 模擬對象

本研究的模擬工況為Goldschmidt 等[19]的擬二維鼓泡床冷態實驗，計算模型如圖1所示。氣體從底部均勻進氣，速度入口，無分布板；氣體出口采用壓力出口。具體的實驗條件及相關模擬參數見表1。

圖1 計算模型的幾何結構Fig.1 Geometry structure of simulation domain

2 數學模型

2.1 CFD-DEM 模型

本研究采用 MFIX 軟件（開源代碼：https://mfix.netl.doe.gov/），基于歐拉-拉格朗日的CFD-DEM 模型對鼓泡床內氣固流動過程進行模擬研究。該模型通過跟蹤流場中每一個顆粒的運動軌跡來模擬整個流場中顆粒運動。采用Gidaspow 曳力模型[20]描述氣-固相兩相作用力。氣相采用層流模型。顆粒相的碰撞采用軟球模型[21-22]描述，將顆粒碰撞視為非彈性碰撞，并考慮摩擦力的存在。關于模型的詳細描述及相關表達式見文獻[19]。

表1 模擬條件Table 1 Simulation conditions

2.2 小波分析方法

函數f(t)在小波基函數φ(t)的連續小波變換定義為f(t)在φ(t)時的卷積[23]

式中，?ab(t) 是 ?(t) 經過平移變換和伸縮變換得到的一系列小波族函數，目前常用的小波基函數有墨西哥帽、法蘭西帽、Morlet 函數和Daubechies函數等[24]。小波變換不僅能表達頻域信息，還能表達時域信息，具有時頻雙局部化的能力。為了定量分析鼓泡床中顆粒渦等相干結構，小波變換是一種很有效的工具，這一點已經在單相湍流上得到了證實[25]。一般分析湍流信息時小波基函數采用Daubechies 函數，所以本研究中采用Daubechies 函數為小波基函數分析鼓泡床中相干結構。

3 結果與討論

3.1 顆粒彈性系數和恢復系數對鼓泡床流場間歇性的影響

Maxwell 曾指出[26]，采用DEM 方法模擬氣固鼓泡流化床時，顆粒彈性系數k一般取400～1500 N·m-1，顆?；謴拖禂礶一般取0.8～1.0。目前，在氣固流化床模擬中k和e一般設為常數[27]，k和e取值的不同會影響床層中顆粒脈動，從而影響整個流場的間歇性。本研究考察了e=0.90 時，k分別為400、800、1000、1200、1500 N·m-1，以及k=1000 N·m-1時，e分別為0.80、0.90、0.95、0.99 時鼓泡床的氣固流動行為。

在鼓泡床中，顆粒速度脈動能和平均床層高度反映了顆粒在碰撞過程中能量的變化，也間接反映出流場的間歇性。鼓泡床的瞬時顆粒平均高度反映了床層隨時間的波動特性。瞬時顆粒平均高度計算公式如下

式中，n是顆粒數目，hk是k顆粒的瞬時高度。

圖2和圖3分別為顆粒在床中心處（X=7.5 cm，Y=7.5 cm）Y方向上的瞬時速度（顆粒速度）脈動能隨彈性系數和恢復系數的變化情況。由圖2可知，隨著彈性系數的增大，在鼓泡床低頻高能區能譜曲線呈上升趨勢；但在高頻耗散區，k為1200、800 N·m-1時的能量高于k為1000 N·m-1時的能量，能譜曲線偏離了Kolmogorov-5/3 的標度律[28]，說明顆粒流場中存在著引起流場間歇性的相干結構。由圖3可知，在高頻耗散區，隨著顆?；謴拖禂档脑龃?，碰撞過程中的能量損耗降低。

圖2 不同彈性系數下顆粒速度脈動能譜圖Fig.2 Fluctuating energy spectrum of particle velocity at different elasticity coefficients (e=0.90)

圖3 不同恢復系數下顆粒速度脈動能譜圖Fig.3 Fluctuating energy spectrum of particle velocity at different restitution coefficients (k=1000 N·m-1)

顆粒碰撞過程中的能量變化會引起平均床層高度的變化，圖4和圖5給出了床層高度隨彈性系數和恢復系數的變化情況。如圖4所示，k為1000 N·m-1時平均床層高度最高，對應圖2中顯示的耗能區能量最低。圖5顯示，隨著恢復系數的增大，高頻區顆粒速度脈動能損耗降低，床層高度增高。

圖4 顆粒彈性系數對平均床層高度的影響Fig.4 Effect of elasticity coefficients on mean bed height (e=0.90)

圖5 顆?；謴拖禂祵ζ骄矊痈叨鹊挠绊慒ig.5 Effect of restitution coefficients on mean bed height (k=1000 N·m-1)

上述顆粒速度脈動能譜分析表明顆粒流場存在較強的間歇性，但是速度脈動能譜圖無法量化流場間歇性的強弱。類似于單相湍流的流場間歇性的 量化指標[18]，前期工作已表明采用平坦因子FF 可以定量描述鼓泡床流場間歇性的強度[29]。對于顆粒的脈動速度，其FF 定義為

式中，v(x) 表示顆粒在Y方向上的脈動速度，表示求平均。當FF＜3，表明顆粒的脈動速度存在較強的周期性；當FF=3，表明顆粒流場不存在間歇結構，顆粒速度脈動是由于隨機脈動造成的；當FF＞3，表明顆粒的脈動速度存在間歇性，而且FF 越大流場的間歇性越強。

為此，本研究進一步通過平坦因子定量分析顆粒彈性系數和恢復系數對流場間歇性的影響。圖6和圖7分別為平坦因子隨彈性系數和恢復系數的變化情況。如圖所示，在不同的k和e時FF 均大于3，說明鼓泡床內流場存在較強的間歇性。隨著k的 增大，平坦因子先降低后增加，當k增大到1000 N·m-1時平坦因子最小，反映了流場間歇性最弱，整體能量耗散最低，對應床層最高，這與圖2、圖4的分析結果一致。隨著e取值的增大，顆粒在動量傳遞過程中的能量損耗降低，導致顆粒速度脈動能降低，流場間歇性也降低，床層高度增加，這也與圖3、圖5的分析結果一致。

圖6 不同顆粒彈性系數對平坦因子FF 的影響Fig.6 Effect of elasticity coefficients on flatness factor (e=0.90)

圖7 不同顆?；謴拖禂祵ζ教挂蜃覨F 的影響Fig.7 Effect of restitution coefficients on flatness factor (k=1000 N·m-1)

3.2 顆粒彈性系數和恢復系數對鼓泡床內相干結構的影響

采用小波變換分析方法[30-31]進一步研究了顆粒彈性系數和恢復系數對鼓泡床流場間歇性以及相干結構的影響。圖8和圖9分別為不同彈性系數和恢復系數下H為18 cm 截面處的平坦因子與頻率之間的關系。由圖可知，在低頻區FF 較低，說明低頻區流場的間歇性較弱，流場的相干結構少。隨著彈性系數的增大，在高頻區k為1000 N·m-1時出現了一個極小值；隨著恢復系數的增大，在高頻區FF 變小，流場的間歇性變化規律與3.1 節的分析一致。根據單相湍流理論，類比于大渦模擬，將鼓泡床顆粒湍動區域分為大渦和小渦，大渦區對應低頻區，小渦區對應高頻區。由于整個顆粒的能量主要集中在大渦區，能量高，大渦在時間和空間上的分布比較均勻，間歇性較弱，不會產生相干性結構；在小渦區，FF 明顯大于3，小渦的不斷擾動存在著大量的相干結構，造成了流場較強的間歇性。

圖8 不同彈性系數的平坦因子與頻率之間的關系Fig.8 Effect of frequency on flatness factor at different elasticity coefficients (e=0.90,H=18 cm)

圖9 不同恢復系數的平坦因子與頻率之間的關系Fig.9 Effect of frequency on flatness factor at different restitution coefficients (k=1000 N·m-1，H=18 cm)

進一步分析了顆粒彈性系數和恢復系數對小波系數的概率密度函數（PDF）分布的影響。Onorato等[32]指出小波系數的PDF 拖尾長度可以表示為相干結構引起的流場間歇性的強弱。如圖10所示，k為1000 N·m-1時主峰較強，拖尾較弱，說明k為1000 N·m-1時流場間歇性較弱。由圖11可以看出，隨著恢復系數的增大，小波系數的主峰增高，拖尾減小，進一步說明恢復系數的增大降低了流場的間歇性。

圖10 不同彈性系數下小波系數的PDF 分布Fig.10 Probability density function of wavelet coefficient at different elasticity coefficients (e=0.90)

圖11 不同恢復系數下小波系數的PDF 分布Fig.11 Probability density function of wavelet coefficient at different restitution coefficients (k=1000 N·m-1)

3.3 顆粒彈性系數和恢復系數對流場相干結構的多尺度的影響

進一步考察了顆粒彈性系數和恢復系數對流場相干結構的多尺度的影響。圖12和圖13為不同彈性系數和恢復系數下不同小波分解尺度下小波系數的相干平面圖，以可視化的方式顯示出了在不同彈性系數和恢復系數下顆粒渦團隨時間的演化過程。由圖可以看出，在渦的中心處相干性最強，相鄰渦團之間相關性最弱的是兩個渦的分界線[33]；鼓泡床中相干結構往往不止一個，大尺度渦和小尺度渦并存，沿著時間的歷程某些尺度的渦會呈現周期性特征，并伴隨著與相鄰尺度渦的合并和分解，說明該時間段內相干渦與鄰近尺度渦存在著強烈的相互作用，導致了渦的聚并和破碎[34]。由圖12可知，在小波尺度為0～20 時，隨著彈性系數的增大，渦中心數目基本一致，這是由于該尺度下相干結構（渦）主要由背景湍流造成。隨著小波尺度的增加，渦團的數目逐漸減少。在60～128 尺度下，渦團的數目急劇減少，但渦的強度大、范圍廣，大渦分解成一部分小渦。在小波尺度為40 左右，隨著彈性系數的增大，渦中心數目先降低后增加，說明流場間歇性先降低后增加，這與圖6平坦因子隨彈性系數先降低后增大的變化趨勢一致。由圖13可知，在小波尺度為0～20 時，隨著恢復系數的增大，渦中心數目基本一致。隨著小波尺度的增加，渦團數目逐漸減少。在60～128 尺度下，渦團的數目也急劇減少，大渦會分解成一部分小渦。在小波尺度為40左右，隨著恢復系數的增大，渦中心數目呈下降的趨勢，說明流場間歇性下降，這與圖7平坦因子隨恢復系數增大而降低的分析一致。

圖12 不同彈性系數下的小波系數自相關平面圖Fig.12 Auto-correlation of wavelet in different elasticity coefficients (e=0.90)

圖13 不同恢復系數下的小波系數自相關平面圖Fig.13 Auto-correlation of wavelet in different restitution coefficients (k=1000 N·m-1)

4 結 論

本研究采用DEM 方法對鼓泡床內氣固流動過程進行了模擬研究，重點分析了顆粒彈性系數和恢復系數對鼓泡床流場間歇性的影響，并利用小波變化方法分析了彈性系數和恢復系數對鼓泡床內相干結構的影響，得到以下結論。

（1）顆粒彈性系數和恢復系數對顆粒速度脈動能、平均床層高度以及平坦因子有一定影響。隨著顆粒彈性系數取值的變大，高頻區能量和平坦因子先降低后增加，床層高度先增加后降低，流場間歇性先減弱后增強。在k為1000 N·m-1時，高頻區能量和平坦因子最低，床層高度最大，對應流場間歇性最弱。顆粒恢復系數越大，高能區能量和平坦因子越低，床層高度越大，流場間歇性越弱。

（2）顆粒彈性系數和恢復系數對小波系數的概率密度函數有一定影響。在k為1000 N·m-1時，主峰較強，拖尾較弱，流場間歇性較弱。隨著恢復系數的增大，小波系數的主峰增強，拖尾減弱，流場間歇性減弱。

（3）顆粒彈性系數和恢復系數對顆粒流場中引起流場間歇性的顆粒渦團有顯著的影響。在高頻區小波尺度為40 附近，隨著彈性系數的增大，渦中心數目先減少后增加，進一步說明流場間歇性先減弱后增強；隨著恢復系數的增大，渦中心數目逐漸減少，流場間歇性逐漸減弱。本研究結果表明，在DEM 方法中顆粒碰撞參數的取值對氣固鼓泡床流場間歇性有較大的影響，針對不同體系選擇合適的顆粒碰撞參數可為更加準確地模擬鼓泡床氣固流動現象提供基礎。

[1]Jin Yong (金涌),Zhu Jingxu (祝京旭),Yu Zhiqing (俞芷青).Fluidization Engineering Principles (流態化工程原理) [M].Beijing:Tsinghua University Press,2001:72-75.

[2]Cundall P A,Strack O D L.A discrete numerical model for granular assemblies [J].Otechnique,1997,29 (1):331-336.

[3]Zhu H P,Zhou Z Y,Yang R Y,et al.Discrete particle simulation of particulate systems:a review of major applications and findings [J].Chemical Engineering Science,2008,63 (23):5728-5770.

[4]Enwald H,Peirano E,Almstedt A E.Eulerian two-phase flow theory applied to fluidization [J].Ⅰnternational Journal of Multiphase Flow,1996,22:21-66.

[5]Tsuji Y,Kawaguchi T,Tanaka T.Discrete particle simulation of two-dimensional fluidized bed [J].Powder Technology,1993,77 (1):79-87.

[6]Kawaguchi T,Tanaka T,Tsuji Y.Numerical simulation of two-dimensional fluidized beds using the discrete element method (comparison between the two- and three-dimensional models) [J].Powder Technology,1998,96 (2):129-138.

[7]Feng Y Q,Yu A B.Assessment of model formulations in the discrete particle simulation of gas-solid flow [J].Ⅰndustrial & Engineering Chemistry Research,2004,43 (26):8378-8390.

[8]Pandit J K,Wang X S,Rhodes M J.Study of Geldart’s Group A behaviour using the discrete element method simulation [J].Powder Technology,2005,160 (1):7-14.

[9]Sun J,Wang J,Yang Y.CFD simulation and wavelet transform analysis of vortex and coherent structure in a gas-solid fluidized bed [J].Chemical Engineering Science,2012,71 (26):507-519.

[10]Sun J,Wang J,Yang Y.CFD investigation of particle fluctuation characteristics of bidisperse mixture in a gas-solid fluidized bed [J].Chemical Engineering Science,2012,82 (12):285-298.

[11]Jiradilok V,Gidaspow D,Breault R W.Computation of gas and solid dispersion coefficients in turbulent risers and bubbling beds [J].Chemical Engineering Science,2007,62 (13):3397-3409.

[12]Fang Mingming (房明明),Luo Kun (羅坤),Yang Shiliang (楊世亮),Zhang Ke (張科),Fan Jianren (樊建人).CFD-DEM investigation and parameter sensitivity analysis of a typical bubbling fluidization process [J].Journal of Engineering Thermophysics(工程熱物理學報),2013,34 (6):1106-1111.

[13]Navarro H A,de Souza Braun M P.Determination of the normal spring stiffness coefficient in the linear spring-dashpot contact model of discrete element method [J].Powder Technology,2013,246:707-722.

[14]Kharaz A H,Gorham D A,Salman A D.An experimental study of the elastic rebound of spheres [J].Powder Technology,2001,120 (3):281-291.

[15]Fu J S,Cheong Y S,Reynolds G K,et al.An experimental study of the variability in the properties and quality of wet granules [J].Powder Technology,2004,140 (3):209-216.

[16]Müller P,Antonyuk S,Tomas J,et al.Investigations of the Restitution Coefficient of Granules—Micro-Macro-Interaction [M].Berlin:Springer,2008:235-241.

[17]Xia Zhenyan (夏振炎),Tian Yan (田硯),Jiang Nan (姜楠).Wavelet spectrum analysis on energy transfer of multi-scale structures in wall turbulence [J].Applied Mathematics and Mechanics(應用數學和力學),2009,30 (4):409-416.

[18]Wu Yingya (吳迎亞),Lan Xingying (藍興英),Gao Jinsen (高金森).Analysis of flow field intermittency and coherent structure of particles based on DEM simulation of gas-solid bubbling bed [J].CⅠESC Journal(化工學報),2014,65 (7):2724-2732.

[19]Goldschmidt M J V,Beetstra R,Kuipers J A M.Hydrodynamic modelling of dense gas-fluidised beds:comparison and validation of 3D discrete particle and continuum models [J].Powder Technology,2004,142 (1):23-47.

[20]Gidaspow D.Multiphase Flow and Fluidization:Continuum and Kinetic Theory Descriptions [M].Boston:Academic Press,1994.

[21]Zhang Yong (張勇),Jin Baoshen (金保升),Zhong Wenqi (鐘文琪).Three-dimensional DEM simulation on particle mixing characteristics of spout-fluid bed [J].Proceedings of the CSEE(中國電機工程學報),2008,28 (2):33-38.

[22]Sommerfeld M,Tsuji Y,Crowe C T.Multiphase Flows with Droplets and Particles [M].Boca Raton:CRC Press,1997:12-14.

[23]Wang Jiajun (王嘉駿),Zhang Wenfeng (張文峰),Feng Lianfang (馮連芳),Gu Xueping (顧雪萍).Wavelets analysis of pressure fluctuation in agitated fluidized bed [J].Journal of Chemical Ⅰndustry and Engineering(China) (化工學報),2006,57 (12):2855-2859.

[24]Mallat S G.A theory for multiresolution signal decomposition:the wavelet representation [J].Pattern Analysis and Machine Ⅰntelligence,ⅠEEE Transactions on,1989,11 (7):674-693.

[25]Giri B K,Mitra C,Panigrahi P K,et al.Multi-scale dynamics of glow discharge plasma through wavelets:self-similar behavior to neutral turbulence and dissipation [OL].[2014-01-13].http://arxiv.org/abs/ 1401.2742.

[26]Zhang Ke (張科).CFD-DEM simulation of complex dense gas-solid flow [D].Hangzhou:Zhejiang University,2012.

[27]Tian F G,Zhang M C,Fan H J,et al.Numerical study on microscopic mixing characteristics in fluidized bedsviaDEM [J].Fuel Process Technology,2007,88 (2):187-198.

[28]Tamburini A,Cipollina A,Micale G,et al.CFD simulations of dense solid-liquid suspensions in baffled stirred tanks:prediction of solid particle distribution [J].Chemical Engineering Journal,2013,223:875-890.

[29]Jiang Nan (姜楠),Tian Yan (田硯).Wavelet analysis of detecting multi-scale coherent structures and intermittency in wall turbulence [J].Journal of Harbin Engineering University(哈爾濱工程大學學報),2005,26 (1):7-12.

[30]Liu Haifeng (劉海峰),Wu Tao (吳韜),Wang Fuchen (王輔臣),Gong Xin (龔欣),Yu Zunhong (于遵宏).Coherent structure of turbulence based on wavelet analysis (Ⅰ):Determining coherent structure with energy maxima criterion [J].Journal of Chemical Ⅰndustry and Engineering(China) (化工學報),2000,51 (6):761-765.

[31]Liu Haifeng (劉海峰),Zhao Tiejun (趙鐵鈞),Wang Fuchen (王輔臣),Gong Xin (龔欣),Yu Zunhong (于遵宏).Coherent structure of turbulence based on wavelet analysis (Ⅱ) :Wave shape and local singularity of coherent structure [J].Journal of Chemical Ⅰndustry and Engineering(China) (化工學報),2000,51 (6):766-770.

[32]Onorato M,Camussi R,Iuso G.Small scale intermittency and bursting in a turbulent channel flow [J].Physical Review E,2000,61 (2):1447.

[33]Zhang Bin (張斌),Wang Tong (王彤),Gu Chuangang (谷傳綱),Dai Zhengyuan (戴正元).Identification of turbulent coherent structures based on wavelet and bi-spectrum analysis [J].Transactions of the Chinese Society for Agricultural Machinery(農業機械學報),2009,40 (11):203-207.

[34]Oberleithner K,Sieber M,Nayeri C N,et al.Three-dimensional coherent structures in a swirling jet undergoing vortex breakdown:stability analysis and empirical mode construction [J].Journal of Fluid Mechanics,2011,679 (1):383-414.