張富翁，王 立，劉傳平

（1.北京科技大學 機械工程學院，北京100083；2.北京科技大學 北京高校節(jié)能與環(huán)保工程研究中心，北京100083；3.攀枝花市經(jīng)濟和信息化委員會綜合處，四川 攀枝花617000）

顆粒物質廣泛存在于自然界中，與人們的生產(chǎn)生活密切相關.作為原材料或催化劑，顆粒物質被廣泛應用于化工、冶金和材料生產(chǎn)行業(yè)中，諸如礦物分選、氣力輸送以及材料處理等［1-3］.目前，顆粒物質已成為力學、凝聚態(tài)物理和軟物質等學科的研究熱點.顆粒物質的生產(chǎn)、加工及儲運約消耗地球當年所用能量的10%，因此，對顆粒物質的深入研究可起到節(jié)能降耗的作用［4］.

面對冶金化工領域日益嚴苛的工藝要求，如何對生產(chǎn)過程中顆粒材料的運動，特別是組分顆粒的分布進行精準控制已成為目前亟待解決的問題.一些圍繞氣固兩相流（流化床或振動流化床）中顆粒運動的研究成果［5-8］為能夠有效預測顆粒的運動和分布奠定了一定的基礎.然而，在振動作用下，顆粒物質卻表現(xiàn)出許多不同于傳統(tǒng)流體或固體的特殊行為，諸如顆粒對流［9］、成拱［10］、表 面波［11］及分離［12］等.對雙組分顆粒而言，在振動激勵下顆粒體系將發(fā)生小顆粒在下大顆粒在上的“巴西果”分離［13-14］或大顆粒在下小顆粒在上的“反巴西果”分離［15］.分離形式不僅取決于顆粒物性，同時也與振動參數(shù)有關.近年來，研究者從顆粒的受力、幾何結構以及表觀現(xiàn)象等方面進行了不懈的探索［16-19］.然而，由于顆粒系統(tǒng)復雜的非線性特性，有關分離本質機理的解釋仍存在分歧.

離散單元（discrete element method，DEM）法［20-22］能夠準確得到振動床內單個顆粒的受力與運動，為分析能量在床層中的傳遞和耗散規(guī)律提供了足夠信息.本文采用離散元模擬，從床層底部輸入振動能量，分析顆粒床層內能量的傳遞和耗散規(guī)律，并討論振動參數(shù)對顆粒系統(tǒng)能量分布的影響規(guī)律.床層采用雙組分顆?；旌衔?，進一步研究了雙組分顆粒的分布對能量傳遞與耗散及系統(tǒng)能量分布的影響，從能量的角度去認識并討論顆粒的振動分離行為.

1 數(shù)學模型

式中：m為顆粒質量；，其中vg為顆粒的速度，ˉvg為顆粒平均速度.而顆粒床層內的能量傳遞則采用傅里葉定律表示：

式中：k為導熱系數(shù).k與床層顆粒數(shù)量N成指數(shù)方關系，k～N y，其中y為指數(shù)系數(shù)，表征床層內顆粒處于不同的狀態(tài)（有序或無序），y的取值范圍為0～0.5［24-25］.根據(jù)床層耗散特性建立床層能量平衡方程

類比于氣體，將顆粒的平均脈動動能定義為顆粒溫度（T）［23］.基于這一擬氣體假設，可以得到如下顆粒溫度的表達式：

式中：顆粒碰撞所引起的床層能量耗散γ與顆粒粒徑dm、彈性恢復系數(shù)e以及顆粒自身溫度T等參數(shù)相關.

模型計算的物理區(qū)域為長方體，長、寬、高為80、80、1 000 mm.計算區(qū)域底部為振動床底，作振幅為A和頻率為f的正弦振動，床體的位移為s=Asin（2πft），t為運動時間.無量綱振動強度（Γ）的表達式為

式中：a為振動加速度，g為重力加速度.床底振動前，將一定數(shù)量的顆粒隨機初始化于計算區(qū)域內.計算區(qū)域四周為周期邊界，顆粒離開側面邊界后會立刻從對面邊界以相同的速度進入.

床體振動過程中，床層內顆粒隨之發(fā)生運動，床內顆粒碰撞時的運動方程為

式中：m為顆粒質量，F(xiàn)n和Ft分別為顆粒碰撞時受到的法向作用力和切向作用力.顆粒與床底之間、顆粒與顆粒之間均為彈性碰撞，碰撞過程中的碰撞力采用Hertz-Mindlin無滑移接觸模型［26］獲得，計算的時間步長取1×10-6s，碰撞示意圖如圖1所示.法向力和切向力分別表示為

式中：α為法向重疊量，vij是顆粒i相對于顆粒j的速度，n為從顆粒i球心指向顆粒j球心的單位矢量，kn和cn分別為顆粒i的法向彈性系數(shù)和法向阻力系數(shù)；δ是接觸點切向位移，vct是接觸點的滑移速度，kt和ct分別為切向彈性系數(shù)和切向阻尼系數(shù).當時，顆粒i發(fā)生滑動，切向力為

式中：μs為摩擦因數(shù).彈性系數(shù)kn和cn的計算公式為

圖1 顆粒間碰撞模型Fig.1 Collision model between grains

下標分別表示碰撞顆粒；E為彈性模量；G為剪切模量；e為彈性恢復系數(shù).以上碰撞模型中，當i或j代表壁面（i，j=w）時，Rw趨于無窮大，R*=1，m*=1.顆粒與壁面的碰撞為完全彈性碰撞，egw=0且μgw=0.顆粒材料的物性參數(shù)為：泊松比γ=0.25，彈性模量E=108，剪切模量G=108，顆粒彈性恢復系數(shù)e=0.25，顆粒間摩擦系數(shù)μgg=0.9.

2 結果與討論

2.1 單組分顆粒體系

將1 440個直徑為12 mm，密度為1 000 kg／m3的圓形顆粒隨機填充于計算區(qū)域內.對床底施加振動后，床內的顆粒將隨之運動.振動一段時間至顆粒運動達到穩(wěn)定狀態(tài)，對顆粒的運動、動能分布以及能量耗散規(guī)律進行分析.高頻振動時可能會產(chǎn)生湍流，如高頻時會出現(xiàn)振動波和倍周期分叉現(xiàn)象.這些現(xiàn)象可能由間隙空氣引起，但是具體機理目前還不明確.由于采用的顆粒直徑較大（粒徑為12 mm），空氣曳力遠小于顆粒慣性力，本研究忽略間隙空氣的影響.

圖2～5分別展示了振幅A為0.003 m，振動頻率f為10 Hz和100 Hz時處于床內不同高度處內顆粒的運動.其中，v z（空心點曲線）和vχ（實心點曲線）分別為顆粒豎直方向和水平方向上的運動速度.在較低頻率區(qū)，當f=10 Hz時，在豎直方向床層顆粒隨著床底振動將周期性地上拋和下落，如圖2和圖3所示，vχ＜v z，床層整體仍然處于“有序”運動狀態(tài).而在高頻區(qū)，當振動頻率f=100 Hz時，如圖4和圖5所示，顆粒的運動規(guī)律性變差，vχ與v z大致相同.此外，沿床高方向上，低頻和高頻振動情況下顆粒的運動劇烈程度也存在明顯差異.

圖2 f=10 Hz且A=0.003 m時床內顆粒豎直方向的運動Fig.2 Vertical motion of grains in vibrated bed where vibration frequency is 10 Hz and amplitude is 0.003 m

圖3 f=10 Hz且A=0.003 m時床內顆粒水平方向的運動Fig.3 Horizontal motion of grains in vibrated bed where vibration frequency is 10 Hz and amplitude is 0.003 m

由圖可知，在低頻率（f=10 Hz）區(qū)，床層上部顆粒運動幅度最大，顆粒運動的劇烈程度隨著顆粒所處高度的增加而增加；而在高頻（f=100 Hz）區(qū)，位于床底處顆粒的運動幅度大、運動最劇烈，沿著床層高度增大的方向顆粒的運動幅度逐漸減小.

圖4 f=100 Hz且A=0.003 m時床內顆粒豎直方向的運動Fig.4 Vertical motion of grains in vibrated bed where vibration frequency is 100 Hz and amplitude is 0.003 m

圖5 f=100 Hz且A=0.003 m時床內顆粒水平方向的運動Fig.5 Horizontal motion of grains in vibrated bed where vibration frequency is 100 Hz and amplitude is 0.003 m

當床層穩(wěn)定后，統(tǒng)計一定時間（6 s）內顆粒的平均動能（Ek），即可得到床內顆粒的動能分布情況.圖6為不同振動條件下顆粒在振動床內的動能分布規(guī)律.當f=10 Hz時，如圖6（a）所示，在振幅較小的情況下，顆粒動能（Ek）隨著床內高度（h）的增加而增加；而隨著振幅的進一步增大，顆粒動能隨著高度不再是單調變化，即在床底和床層上部顆粒的動能較高，而在床層中部顆粒的動能相對較低.當f=100 Hz時，如圖6（b）所示，在不同的振幅條件下，顆粒動能呈現(xiàn)出相同的變化趨勢.顆粒動能隨著高度的增加而單調減小，且隨著振幅的增加，顆粒動能越來越大.

圖6 不同振動參數(shù)下振動床內顆粒動能的分布規(guī)律Fig.6 Kinetic energy distribution in vibrated bed with different vibration parameters

顆粒間的彈性恢復系數(shù)e是影響床層能量分布規(guī)律的重要因素.圖7為f=10 Hz且A=0.010 m時不同彈性恢復系數(shù)下顆粒床層的動能分布.隨著e的增大，顆粒體系的總動能逐漸增大.當e＜0.75時，顆粒體系動能沿著床高（h）方向有著幾乎一致的分布規(guī)律：顆粒動能在床層中部較低，而在床底和床層上部較高.當e=0.95時，顆粒體系動能分布有所改變，沿著床高增大的方向顆粒動能逐漸增大.

當床層振動時，顆粒之間產(chǎn)生頻繁的相互碰撞，顆粒碰撞前后存在動能差這一動能差即為每次碰撞引起的動能耗散.將單位時間單位體積內床層顆粒碰撞的總能量耗散定義為床層的能量耗散速率，即.圖8為不同振動參數(shù)下顆粒的能量耗散分布規(guī)律，其中γ*=γ／γav，γav為整個床層的平均能量耗散速率.

圖7 不同彈性恢復系數(shù)下振動床內顆粒動能的分布規(guī)律Fig.7 Kinetic energy distribution in vibrated bed with different coefficients of restitution

當Γ＜1時，振動過程中床內顆粒未完全脫離床底，即顆粒與床底一起作上下運動，床層內各個高度動能耗散基本一致，如圖8中實心圓點所示；而當無量綱振動強度增加（頻率增加或振幅增加）后，系統(tǒng)中能量耗散梯度呈負梯度，即床底能量耗散遠多于床層上部.高頻率和大振幅床層中能量耗散均沿著床高的增加而減小，然而在量值上卻存在差異.高頻率時床底與床層上部能量耗散差異大，而高振幅時差異小.當為大振幅振動時，由于單位時間內的碰撞次數(shù)較高頻振動時少，整個床層由于碰撞引起的能量耗散差異不大；而振動為高頻率時，床底顆粒劇烈運動、床層上部顆粒運動幅度相對較小，床底能量耗散速率遠大于上部床層.

如圖8所示，隨著f或A的增加，顆粒體系的平均能量耗散速率（γav）亦隨之增加.這與文獻［27］中的模擬和文獻［28］中的實驗研究結果有著一致的變化趨勢：當無量綱振動強度一定，在f＜100 Hz時，隨著f的增大，床層的平均能量耗散功率（耗散能量）逐漸增大.彭政等［28］根據(jù)能量守恒，將顆粒系統(tǒng)的平均能量耗散功率視作振動臺向顆粒床進行能量傳遞的平均功率W（忽略空氣的影響）.當振幅一定而頻率增大時，振動變得更加劇烈，床底向床層輸入的總能量增加.相應地，床層內顆粒間的非彈性碰撞次數(shù)增加，體系的能量耗散隨之增大.

2.2 雙組分顆粒體系

將2種密度相同尺寸不同的顆粒（大顆粒：直徑為12 mm，密度為1 000 kg／m3，數(shù)量為720顆；小顆粒：直徑為8 mm，密度為1 000 kg／m3，數(shù)量為2 430顆，大小顆?？傮w積之比1：1）先后初始化于床層中，得到2種不同的分布形式.“巴西果”分布（BN distribution），床層下部為小顆粒、上部為大顆粒；“反巴西果”分布（RBN distribution），床層下部為大顆粒，上部為小顆粒.對顆粒床層施加豎直方向振動，在不同的振動條件下，統(tǒng)計一定時間（6 s）內顆粒的平均動能Ek.在統(tǒng)計時間段內，顆粒體系基本保持初始分布，統(tǒng)計過程中并不涉及顆粒分布狀態(tài)的改變.圖9為雙組分顆粒處于不同分布時的動能分布規(guī)律.通過比較，可以得到以下結論.

圖8 床層內能量耗散速率分布Fig.8 Distribution of energy dissipation rate in granular bed

圖9 雙組分顆粒振動床內動能分布Fig.9 Kinetic energy distribution in binary granular mixture

1）相同位置大顆粒動能大于小顆粒動能；

2）當為低頻率（f=10 Hz）振動時，同一組分顆粒動能沿床高方向上為正梯度分布，即顆粒動能隨高度的增加而增加.當振動為高頻率（f=100 Hz）時，同一組分顆粒動能沿床高方向上為負梯度，即顆粒動能隨高度的增加而降低.

3）雙組分顆粒分布不同時，床層內顆粒的動能分布形式存在顯著差異，床層顆?？倓幽芤膊幌嗤?當?shù)皖l振動時，“反巴西果”分布時床層顆?？倓幽芨停划敻哳l振動時，“巴西果”分布時床層顆粒總動能更低.

改變振幅與頻率，比較“巴西果”和“反巴西果”2種不同分布時系統(tǒng)的動能（分別為BN分布總能量EBN和RBN分布總能量ERBN），可得到圖10所示的相圖.在A區(qū)域，Γ＜1，顆粒運動不會離開床層，在B區(qū)域（低頻率、小振幅），當顆粒處于反巴西果分布時，顆粒床層的動能更低，ERBN＜EBN；在C區(qū)域（高頻率、大振幅），當顆粒處于巴西果分布時床層動能更低，EBN＜ERBN.模擬得到2種分布情況下的系統(tǒng)總能量，總能量越低其分布越穩(wěn)定，按照模擬結果進行分區(qū).A區(qū)和B區(qū)的劃分條件為Γ=1，因此，A區(qū)和B區(qū)分界線（見圖10中虛線）的數(shù)學表達式即為A=g／（4π2f2）.通過非線性多項式擬合，可以得到B區(qū)和C區(qū)的分界線（見圖8中實線）的表達式為

圖10 顆粒分布對床層動能的影響Fig.10 Effect of grain size distribution on kinetic energy

顆粒直徑越小，單位體積內顆粒數(shù)量越多，相同條件下碰撞次數(shù)越多，因此，相比于大顆粒床層而言，小顆粒床層對外部輸入能量具有更強的耗散能力.沿床層床高方向顆粒動能分布不同，存在顆粒動能相對更大的“高動能區(qū)”.當耗散能量強的小顆粒處于床層“高動能區(qū)”時，整個床層對外界輸入能量具有更強的耗散能力，使得整個床層自身動能更低.

對于低頻率及低振幅振動，床層“高動能區(qū)”處于床層上部，因此雙組分顆粒處于“反巴西果”分布時床層總動能更低；而對于高頻率振動，床層“高動能區(qū)”處于床底，雙組分顆粒處于“巴西果”分布時床層總動能更低.

3 結 論

（1）作為顆粒體系微觀結構研究的有效工具，離散單元法可以有效地用于顆粒體系的能量傳遞與耗散的模擬研究中.結果表明：當振動參數(shù)（振幅和頻率）發(fā)生改變時，沿著床層的高度方向顆粒的能量分布、傳遞與耗散隨之發(fā)生改變.低頻振動時，碰撞引起的能量耗散沿床高方向逐漸增加；高頻振動時，沿床高方向顆粒運動引起的能量耗散迅速減小.

（2）雙組分顆粒放入振動床層中，系統(tǒng)的總動能隨著振動頻率及振幅的增加而增大.同時，雙組分顆粒在床內的分布也會影響床層總動能：在低振幅、低頻率時，“反巴西果”分布時床層總動能較低，反之則“巴西果”分布時床層總動能較低.

（3）顆粒系統(tǒng)為自身動能分布滿足某一條件極值（例如：總能量最低），而調整床層內雙組分顆粒的分布，最終處于穩(wěn)定的分離狀態(tài).“巴西果”分布和“反巴西果”分布時床層具有的總動能差異是否為顆粒分離的原始驅動力，或與顆粒分離的驅動力之間存在何種聯(lián)系，有待進一步研究.

（

）:

［1］JAEGER H M，NAGEL S R，BEHRINGER R P.Granular solids，liquids，and gases［J］.Reviews of Modern Physics，1996，68：1259- 1273.

［2］KUANG S B，YU A B，ZOU Z S.Computational study of flow regimes in vertical pneumatic conveying［J］.Industrial and Engineering Chemistry Research，2009，48：6846- 6858.

［3］鮑德松，張訓生.顆粒物質與顆粒流［J］.浙江大學學報：理學版，2003，30（5）：514- 517.BAO De-song，ZHANG Xun-sheng.Granular matter and granular flow［J］.Journal of Zhejiang University：Science Edition，2003，30（5）：514- 517.

［4］陸坤權，劉寄星.顆粒物質［J］.物理，2004，33（9）：629- 635.LU Kun-quan，LIU Ji-xing.Granular matter［J］.Physics，2004，33（9）：629- 635.

［5］劉偉偉，盧春喜，范怡平，等.氣固流化床中雙組分混合顆粒的流態(tài)化特性［J］.化工學報，2008，59（8）：1971- 1978.LIU Wei-wei，LU Chun-xi，F(xiàn)AN Yi-ping，et al.Flow behavior of binary mixture particles in gas-solid fluidized beds［J］.Journal of Chemical Industry and Engineering，2008，59（8）：1971- 1978.

［6］張樹青，盧春喜，時銘顯，等.氣固流化床中大差異雙組份顆粒分級特性的實驗研究［J］.高?；瘜W工程學報，2004，55（10）：1581- 1592.ZHANG Shu-qing，LU Chun-xi，SHI Ming-xian，et al.Segregation of binary particle with significant size difference in gas-solid fluidized beds［J］.Journal of Chemical Engineering of Chinese Universities，2004，55（10）：1581- 1592.

［7］江茂強，趙永志，鄭津洋.非等密度顆粒氣固流化床的微觀尺度模擬與分析［J］.浙江大學學報：工學版，2009，43（9）：1703- 1708.JIANG Mao-qiang，ZHAO Yong-zhi，ZHENG Jinyang.Micro-scale simulation and analysis of gas-solid fluidized bed with multi-density distribution of particles［J］.Journal of Zhejiang University：Engineering Science，2009，43（9）：1703- 1708.

［8］DEEN N G，WILLEM G，SANDER G，et al.Numerical analysis of solids mixing in pressurized fluidized beds［J］.Industrial and Engineering Chemistry Research，2010，49：5246- 5253.

［9］KHAIN E，MEERSON B.Onset of thermal convection in a horizontal layer of granular gas［J］.Physical Review E，2003，67：021306.

［10］HSIAU S S，WU M H，CHEN C H.Arching phenomena in a vibrated granular bed ［J］.Powder Technology，1998，99：185- 193.

［11］AOKI K M，AKIYAMA T.Spontaneous wave pattern formation in vibrated granular materials［J］.Physical Review Letters，1996，77：4166- 4419.

［12］YUAN X，ZHENG N，SHI Q，et al.Segregation in mixtures of granular chains and spherical grains under vertical vibration［J］.Physical Review E，2013，87：042203.

［13］ROSATO A D，STRANDBURG K J，PRINZ F，et al.Why the Brazil nuts are on top：Size segregation of particulate matter by shaking ［J］.Physical Review Letters，1987，58：1038- 1040.

［14］HSIAU SS，CHEN W.Density effect of binary mixtures on the segregation process in a vertical shaker［J］.Advanced Powder Technology，2002，13：301- 315.

［15］SHINBORT T，MUZZIO F J.Reverse buoyancy in shaken granular beds［J］.Physical Review Letters，1998，81：4365- 4368.

［16］LIU C P，WANG L，WU P，et al.Effects of gas flow on granular size separation ［J］.Physical Review Letters，2010，104：188001.

［17］ROSATO A D，BLACKMORE D L，ZHANG N，et al.A perspective on vibration-induced size segregation of granular materials ［J］.Chemical Engineering Science，2002，57：265- 275.

［18］LIU C P，WANG L，WU P，et al.Size distribution in gas vibration bed and its application on grain drying［J］.Powder Technology，2012，221：192- 198.

［19］林建忠，游振江，石興.混合層中柱狀粒子運動的研究［J］.工程熱物理學報，2003，24（6）：972- 975.LIN Jian-zhong，YOU Zhen-jiang，SHI Xing.Research on the motion of cylindrical particle in mixing layer［J］.Journal of Engineering Thermophysics，2003，24（6）：972- 975.

［20］張國華，孫其誠，黃芳芳，等.摩擦顆粒體系各向同性壓縮過程中的堵塞行為［J］.物理學報，2011，60（12）：124502.Zh ANG Guo-hua，SUN Qi-cheng，HUANG Fangfang，et al.Jamming phenomena of a two-dimensional frictional granular system under isotropic confining［J］.Acta Physica Sinica，2011，60（12）：124502.

［21］吳錦坤，羅坤，胡桂林，等.鼓泡流化床流動特性的直接數(shù)值模擬［J］.浙江大學學報：工學版，2007，41（3）：504- 508.WU Jin-kun，LUO Kun，HU Gui-lin，et al.Direct particle simulation of flow characteristics in bubbling fluidized bed［J］.Journal of Zhejiang University：Engineering Science，2007，41（3）：504- 508.

［22］趙永志，程易，金涌.提升管與下行床顆粒團聚行為的離散顆粒模擬［J］.化工學報，2007，58（1）：44- 53.ZHAO Yong-zhi，CHENG Yi，JIN Yong.CFD-DEM simulation of clustering phenomena in riser and downer［J］.Journal of Chemical Industry and Engineering，2007，58（1）：44- 53.

［23］WILDMAN R D，HUNTLEY J M.Scaling exponents for energy transport and dissipation in binary vibro-fluidized granular beds［J］.Physics of Fluids，2003，15：3090- 3098.

［24］NARAYAN O，RAMASWAMY S.Anomalous heat conduction in one-dimensional momentum-conserving systems［J］.Physical Review Letters，2002，89：200601.

［25］DHAR A，SAITO K.Heat conduction in the disordered Fermi-Pasta-Ulam chain［J］.Physical Review E，2008，78：061136.

［26］MINDLIN R D.Compliance of elastic bodies in contact［J］.Journal of applied mechanics，1949，16：259- 268.

［27］YANAGIDA T，MATCHETT A J，COULTHARD J M.Dissipation energy of powder beds subject to vibration［J］.Chemical Engineering Research and Design，2001，79：655- 662.

［28］彭政，蔣亦民，劉銳，等.垂直振動激發(fā)下顆粒物質的能量耗散［J］.物理學報，2013，62（2）：024502.PENG Zheng，JIANG Yi-ming，LIU Rui，et al.Energy dissipation of a granular system under vertical vibration［J］.Acta Physica Sinica，2013，62（2）：024502.