顆粒毛細效應影響因素的離散元分析1)

張華騰 凡鳳仙 王志強

(上海理工大學能源與動力工程學院,上海 200093)

(上海理工大學上海市動力工程多相流動與傳熱重點實驗室,上海 200093)

引言

在外加擾動作用下顆粒物質的敏感性和非線性響應等使得顆粒物質表現出復雜而奇特的動力學行為,日益吸弓著不同領域科研工作者的興趣[1-6].近幾年,文獻報道了顆粒的毛細效應[4-6],即將一根細管插入填充有顆粒物質的容器中并對管施加豎直振動時,顆粒物質在管內迅速上升并最終達到一個穩定的高度的現象.該現象為實現顆粒物料的逆重力輸送提供了一種潛在的技術途徑,弓起了研究者的高度關注.

自Liu 等[7]對顆粒毛細效應進行實驗報道以來,研究者迅速開展了針對顆粒毛細效應的實驗、理論和數值模擬研究.首先,在實驗研究方面,獲得了管振動強度、振幅、頻率、管插入深度、管口形狀、管內顆粒填充高度、間隙氣體、管內徑、管截面形狀、顆粒直徑對顆粒上升高度的影響[7-11];對比了管頂開口與封閉兩種情況下顆粒的上升過程,發現管頂封閉對顆粒上升起促進作用,分析了管徑及管長對這種促進作用的影響[12].其次,在理論研究方面,將管內顆粒視為一個整體,建立了顆粒受力模型,對顆粒的上升高度和速度進行了理論分析,并提出了顆粒上升的“空穴填充”機理[7];建立了顆粒上升高度隨時間變化的半經驗關系式,進一步分析了顆粒上升的“空穴填充”機理,并且基于能量平衡假設,推導出顆粒最終毛細上升高度的理論表達式[13];基于一個振動周期內管內顆粒的流進和流出特性,建立了顆粒毛細上升高度的理論模型,該模型能夠在一定程度上解釋振幅、頻率和管徑對顆粒最終毛細上升高度的影響,但無法解釋較大振動頻率下顆粒最終毛細上升高度隨頻率下降的現象,以及較小管徑下顆粒最終毛細上升高度隨管徑增加的現象[14].再者,在數值模擬方面,利用離散元方法,獲得了顆粒毛細效應全過程[15];揭示了顆粒毛細效應的對流機理,并發現了顆粒最終毛細上升高度隨管徑的變化特性[4];展現了顆粒在豎直方向的運動規律[5];獲得了顆粒填充率的變化特性,以及振幅、振動強度、恢復系數對顆粒上升過程的影響規律[16];展示了不同初始填充高度下顆粒毛細效應的過程,探討了顆粒毛細上升過程中管內顆粒柱高度、顆粒速度場和填充率分布隨時間的演變特性,獲得了顆粒毛細效應過程中由容器傳輸到管內的顆粒的占比分布[6].

雖然已有研究在顆粒毛細效應的宏觀效果、動力學過程和內在機理上取得了諸多進展.然而,由于顆粒系統參數的復雜性和多樣性,目前顆粒毛細效應研究尚處于起步階段,宏觀的實驗研究仍不系統,尚缺乏對不同操作參數下管徑影響的探討,而這方面研究對推動顆粒毛細效應在顆粒逆重力輸運上的應用具有重要意義.此外,已有研究中提出的顆粒毛細上升高度的理論模型都是建立在把大量顆粒視為整體的基礎之上,事實上管內不同高度和徑向位置顆粒的運動速度往往不同,因此現有理論模型無法對顆粒的輸運提供有效的定量指導.離散元方法通過跟蹤每一個離散單元(顆粒)的運動以實現對整個顆粒系統的模擬,是從顆粒尺度探究顆粒動力學行為的行之有效的方法[17-19].本文將在課題組前期工作的基礎上[4-6,15],利用離散元方法,開展顆粒毛細效應過程數值模擬,獲得不同容器寬度、管振幅和頻率下顆粒毛細上升高度隨管徑的變化規律,以探明顆粒毛細效應的影響因素,并為顆粒物料逆重力輸運方案的優化提供科學依據.

1 離散元模型與計算方法

1.1 離散元模型

在離散元模型中,顆粒系中任意顆粒i的平動與轉動方程由牛頓第二定律給出,即

式中,mi和Ii分別為顆粒i的質量和轉動慣量;vi和ωi分別為顆粒i的速度和角速度;t為時間;n為與顆粒i相接觸的顆粒的個數;Fij,n和Fi j,t分別為顆粒i和與它相接觸的顆粒j之間的法向作用力和切向作用力;Mi j,t為顆粒j對顆粒i的切向作用力產生的力矩;Mij,r為顆粒j對顆粒i的滾動摩擦力產生的力矩;g為重力加速度.

采用黏彈性接觸模型[20]描述法向作用力,采用修正的Cundall-Strack 模型[21]描述切向作用力,則

式中,ρ 為彈性參數,其是楊氏模量Y、泊松比υ、有效半徑Reff=RiRj/(Ri+Rj)的函數;Ri和Rj分別為顆粒i和j的半徑;δi j,n為顆粒i與j之間的法向重疊量,δi j,n=Ri+Rj-|ri-rj|;ri和rj分別為顆粒i和j的位置矢量;An為法向耗散系數,根據其與恢復系數ε 之間的關系[22-24],可采用Pad′e近似方法進行計算[25];en為法向單位矢量,en=(rj-ri)/|rj-ri|;μs為滑動摩擦系數;G為剪切模量,G=Y/[2(1+υ)];At為切向耗散系數,At=AnY/(1-υ2)[26];vij,t為接觸點上顆粒i與j的切向相對速度;et為切向單位矢量.此外,式(4)中的積分路徑為兩顆粒接觸期間在接觸點的相對位移.

彈性參數ρ 的計算式為

根據定向恒轉矩模型[27],Mi j,r可寫為

式中,μr為滾動摩擦系數;ωi j為顆粒i和j的相對角速度,ωi j=ωi-ωj.

對于顆粒i與管壁面及容器壁面之間的相互作用,將壁面視為半徑無限大的顆粒,利用上述模型進行處理[28-30].

1.2 計算方法

基于DEM 模型,借助開源顆粒系統數值模擬軟件LIGGGHTS[31]對顆粒毛細效應動力學行為開展數值模擬.本文數值模擬采用的顆粒系統示意圖如圖1 所示.其中,容器的橫截面為正方形,其邊長為w;容器底面位于z=-18 mm 處,容器高度足夠高;管豎直插入容器中心,管底面中心位于坐標原點,管內徑為D,管壁厚為0.3 mm,管長度足夠長;顆粒直徑d=0.6 mm.數值模擬采用的顆粒物性參數見表1.數值模擬過程為:首先將一定數目的顆粒隨機加入容器內(含管內區域),顆粒在重力作用下發生沉降,由于重力作用、顆粒與壁面以及顆粒與顆粒之間的相互作用,系統內顆粒總動能先增大后減小,趨于總動能為0 的松弛狀態;待顆粒達到松弛狀態后,對管施加豎直方向的周期性正弦振動,使得管底面位置zb隨管振動時間t的變化規律為

式中,a為振幅,f為頻率.

數值模擬時,時間步長的選擇非常重要,時間步長越小,計算精度越高,但計算時間成本也越高.因而,在選擇時間步長時,應兼顧計算精度和計算成本.在離散元模擬中,通常選擇顆粒碰撞時間的1/50～1/10 作為時間步長[32-33].為確定時間步長,采用無阻尼、無黏性碰撞模型估算顆粒碰撞時間[34]

圖1 數值模擬采用的顆粒系統示意圖Fig.1 Schematic of the granular system used in numerical simulations

表1 數值模擬采用的顆粒物性參數Table 1 Physical properties of the particles used in numerical simulations

式中,meff為顆粒對的有效質量,meff=mimj/(mi+mj);vimp為碰撞速度.取vimp=1 m/s 作為碰撞速度參考值,則有tcol≈6.0 × 10-5s.因此,選擇時間步長Δt=tcol/30 ≈2.0×10-6s.在數值模擬時,發現采用更小的時間步長,數值模擬結果不發生改變,這表明時間步長的選擇是合理的.

2 結果與討論

2.1 不同管徑下顆粒毛細效應過程

在容器寬度與粒徑比w/d=40、管振幅與粒徑比a/d=14.33、管振動頻率f=12 Hz 時,得到的典型管徑與粒徑比(D/d)下,顆粒毛細效應過程快照如圖2 所示.圖2(a)～圖2(c)分別為D/d=3.33,8.33,15時的結果.圖中,T為管振動周期,T=1/f;Hc為毛細上升高度;豎直方向速度接近0(|vz| ＜0.05 m·s-1)的顆粒用藍色表示,其余顆粒中,向上運動的用紅色、向下運動的用綠色表示.由圖2(a)可見,在D/d很小時,由于堵塞嚴重,毛細上升速度較慢,且管內顆粒發生結團現象,使得顆粒柱出現中斷.需要說明的是:圖中給出的是振動周期整數倍時刻的結果,此時管內顆粒處于最密堆積狀態[4,6],為結團和中斷最弱的時刻,非整數倍周期時刻結團和中斷更為嚴重.由圖2(b)可以看出,D/d增大到8.33 時,起初毛細上升高度增大迅速,隨后毛細上升高度的增大逐漸減緩,直至達到一個穩定的最終毛細上升高度;由于D/d仍較小,受堵塞影響,管徑方向幾乎不存在速度梯度;在12T時刻,顆粒柱頂部的少量顆粒因相互作用較弱,自由程較大,運動行為不同于管內其余顆粒而向下運動,管內其余顆粒間頻繁的摩擦與非彈性碰撞不斷消耗顆粒系統的動能,使得顆粒聚集在一起,顆粒自由程較小,顆粒呈現隨管向上運動的狀態;圖中24T～108T各時刻,管內顆粒均向上運動,這與顆粒體系的能量耗散有關,經過更多振動周期,顆粒間的摩擦與非彈性碰撞使得顆粒柱頂部顆粒呈密堆積狀態,管內顆粒柱頂部與下部的顆粒運動協同起來,因而顆粒柱整體向上運動.由圖2(c)可知,在D/d=15 條件下,顆粒毛細上升高度先迅速增大,而后在t≥24T時出現起伏,這與圖2(b)差異顯著.由圖2(c)還可以看出,t≥24T時,管內顆粒柱分離為速度截然不同的兩層:管以最大速度向上運動時,管內上層顆粒向下運動,與管的運動方向相反;管內下層顆粒隨管向上運動,且管壁附近顆粒的速度大于管中心區域顆粒的速度,在管徑方向存在明顯的速度梯度.

圖2 數值模擬得到的顆粒毛細效應過程快照Fig.2 Snapshots of the granular capillarity process obtained by numerical simulations

圖3 不同管徑下顆粒毛細上升高度隨時間的演變Fig.3 Evolution of capillary height of the granular materials with time at different tube diameters

為定量描述顆粒毛細上升高度(即管內顆粒表面高度和容器內顆粒表面高度之差)隨時間的演變,將管內z＞0 區域沿z軸劃分為高度為5d的連續的單元,計算各單元內顆粒填充率φ,并定義滿足φ ＞0.1 的單元的最大高度為管內顆粒表面高度.類似地,將容器內z＞0,r＞D/2 區域沿z軸劃分為高度為d的連續的單元,計算各單元內顆粒填充率,采用相同的方法確定容器內顆粒表面高度.在此基礎上,即可獲得顆粒毛細上升高度.不同D/d條件下顆粒毛細上升高度隨時間的演變如圖3 所示.其中,圖3(a)給出了D/d=5,6.67,8.33 時的結果,圖3(b)給出了D/d=10,15,20,30 時的結果.從圖3 可以清楚地看到,顆粒毛細效應過程中伴隨著顆粒毛細上升高度的波動,且顆粒毛細上升速度和最終毛細上升高度受到管徑的影響.由圖3(a)可知,管徑較小時,隨著管徑的增加,顆粒毛細上升速度更快,所能達到的最終毛細上升高度也更大.總體看來,管徑較大時,隨著管徑的增加,顆粒毛細上升速度減小,所能達到的最終毛細上升高度也更小,如圖3(b)所示.同時,需要指出的是:D/d=15 時,顆粒毛細上升高度呈現出反復波動的非穩定特征,這一特征產生的原因可能是豎直振動過程中顆粒的填充率分布和堆積結構、顆粒與顆粒以及顆粒與壁面之間的碰撞等因素影響了顆粒的自組織行為.復雜體系中顆粒物質的行為模式對其內部結構及擾動具有高敏感性,其物理機制十分復雜,有待進一步深入研究.

2.2 臨界管徑的影響因素分析

Fan 等[4]在僅改變管徑而保持其他參數不變時,發現在毛細效應能夠發生的管徑范圍內,顆粒最終毛細上升高度隨管徑的增大,先增大后減小.為了敘述的方便,本文將顆粒最終毛細上升高度隨管徑增大而增大的管徑區間的上限稱為臨界管徑.基于此,臨界管徑下顆粒最終毛細上升高度達到最大,一旦管徑大于臨界管徑,顆粒最終毛細上升高度下降.臨界管徑影響因素研究對于利用顆粒毛細效應輸運顆粒物料的效果優化具有重要意義.

2.2.1 容器寬度的影響

為探究容器寬度對臨界管徑的影響,在a/d=14.33,f=12 Hz 條件下,對容器寬度與粒徑比w/d=40,60 時顆粒毛細效應過程開展數值模擬,獲得顆粒最終毛細上升高度隨D/d的變化關系,如圖4 所示.為體現顆粒最終毛細上升高度隨管的振動而上下波動的特性,圖中利用點表示顆粒最終毛細上升高度的平均值,利用誤差條表示標準差.由圖可知,w/d=40 情況下,臨界管徑對應的D/d=9.17,在D/d≥10 時,顆粒表現出液體的性質,由于管內顆粒的流化,阻礙了顆粒的上升,使得在臨界管徑附近顆粒最終毛細上升高度急劇下降.確切地說,w/d=40 情況下,在D/d=9.17 時,在D/d=10 時,兩者相比顆粒最終毛細上升高度下降了67.6 mm.在w/d=60 情況下,顆粒最終毛細上升高度隨著D/d的增大,同樣呈現出先增大后減小的特征,并且在管徑超出臨界管徑時顆粒最終毛細上升高度急劇下降.然而,兩種容器寬度下,臨界管徑并不相同,在w/d=60 情況下,臨界管徑對應的D/d=10,相應的顆粒最終毛細上升高度為138.1 mm;在D/d=11.67 時,顆粒最終毛細上升高度為60.7 mm,與D/d=10 時相比下降了77.4 mm.究其原因是:有摩擦的容器側壁邊界的存在對于誘導容器內顆粒的對流起到了不可缺少的作用[4],隨著容器寬度的增大,容器側壁的邊界效應降低,容器內顆粒的對流強度隨之減小,這造成了由顆粒對流傳輸至管口的顆粒的質量通量的減小,從而削弱了管內顆粒的流化,因此,顆粒由堵塞到流化的轉折點延后,表現為臨界管徑隨容器寬度的增大而增大.

圖4 容器寬度對臨界管徑的影響Fig.4 Effect of container width on critical tube diameter

2.2.2 振幅的影響

圖5 振幅對臨界管徑的影響Fig.5 Effect of vibration amplitude on critical tube diameter

為探究振幅對臨界管徑的影響,在w/d=40 的條件下對振幅與粒徑比a/d=14.33 與20 兩種情況下顆粒毛細效應開展數值模擬,得到顆粒最終毛細上升高度H∞c隨管徑與粒徑比D/d的變化情況,如圖5 所示.在a/d=20 情況下,臨界管徑對應的D/d=11.67,對比顯示,該情況下的臨界管徑明顯大于a/d=14.33 時的臨界管徑;同時,顆粒最終毛細上升高度H∞c的最大值可達到203.2 mm,也明顯大于a/d=14.33 時H∞c的最大值133.5 mm.振幅增大弓起臨界管徑增大的原因是:振幅增大,管的振動強度增大,管內壁附近的顆粒受到了更大的壁面剪切力作用,使得更多顆粒能夠跟隨管運動,減小了管內顆粒自由流動的區域,相當于減少了有效管徑,因而增強了管內顆粒的堵塞效應,使得顆粒堵塞能夠在更大的管徑下影響顆粒毛細效應.在D/d=5 時,a/d=20情況下顆粒最終毛細上升高度比a/d=14.33 情況下略低,這同樣能夠通過顆粒堵塞進行解釋.由于振幅較大時,更易出現堵塞,管徑很小時,管內顆粒堵塞嚴重抑制了顆粒的上升,使得最終顆粒毛細上升高度下降,特別是當管徑過小時,顆粒出現完全堵塞,顆粒毛細效應將無法發生[4].此外,D/d＞5 時,a/d=20 情況下的顆粒最終毛細上升高度均大于a/d=14.33 時的值,這與文獻中實驗報道的顆粒最終上升高度隨振幅的增加而增加的現象一致[11],該現象可解釋為振幅增加弓起管的振動強度增大,使得管外壁對容器內的顆粒產生了更大的剪切力作用,促進了容器內顆粒的對流,增強了容器內顆粒向管口的傳輸,有效地提高了顆粒毛細上升高度.

2.2.3 頻率的影響

圖6 頻率對臨界管徑的影響Fig.6 Effect of vibration frequency on critical tube diameter

在w/d=40,a/d=14.33 條件下,考察f=12 Hz與f=15 Hz 時最終穩定毛細上升高度隨D/d的變化情況,結果如圖6 所示.對比發現,D/d≤8.33 時,在較高頻率下,略低;D/d≥9.17 時,在較高頻率下,也較高.特別是當10 ≤D/d≤15 時,兩種頻率下顆粒最終毛細上升高度差異顯著:f=15 Hz 條件下,顆粒毛細效應動力學處于受堵塞影響的區域,保持很高的值;f=12 Hz 條件下,顆粒毛細效應動力學處于不受顆粒堵塞影響的流化區域,較低.與振幅相比,頻率對顆粒毛細效應的影響更為復雜,頻率升高弓起振動強度增大、振動周期減小.由于本文研究中采用的兩種頻率相差不大,振動強度的增大對顆粒毛細效應動力學的影響起主導作用.較高頻率下,振動強度也較大,靠近管內壁的顆粒因受到更大的壁面剪切力作用更難以與管壁分離,導致管內顆粒的有效流動區域變小,使得管內顆粒更容易發生堵塞.頻率較高時,更為嚴重的堵塞效應造成較小管徑(D/d≤8.33)下顆粒更難以上升,從而產生更小的最終毛細上升高度.同時,頻率較高時,堵塞效應能夠在更大的管徑(D/d≤15)下發揮作用,使得臨界管徑增大.此外,由于頻率增大,增加了管的振動強度,因而增強了管外顆粒在對流機制作用下向管口的傳輸,導致D/d＞15 時,較高頻率下,顆粒最終毛細上升高度較大.然而,需要說明的是:若振動頻率在很大的范圍內變化,則管振動周期的改變將可能對顆粒動力學行為起主導作用.來自文獻的實驗結果表明,保持振幅不變時,隨著頻率的增大,顆粒最終上升高度先增加后減小,存在一個使得最大的最佳頻率值[11].更大的頻率變化范圍內,臨界管徑的變化特性,有待進一步研究.

3 結論

基于離散元方法,對顆粒毛細效應動力學行為開展數值模擬研究,給出了不同管徑下顆粒毛細效應過程中顆粒豎直方向速度和毛細上升高度隨時間的演變規律,分析了容器寬度、管的振幅和頻率對顆粒最終毛細上升高度隨管徑變化的影響特性.通過本文研究,得到如下結論:

(1)在容器寬度與粒徑比為40、管振幅與粒徑比為14.33、管振動頻率為12 Hz 情況下,管徑與粒徑比D/d=3.33 時,管內顆粒堵塞嚴重,使得顆粒上升緩慢,并造成顆粒柱中斷;D/d=8.33 時,起初顆粒毛細上升高度迅速增加,隨后顆粒毛細上升高度增速減緩,管內顆粒在管徑方向幾乎不存在速度梯度;D/d=15 時,隨著顆粒毛細上升高度的增大,管內顆粒柱分離為速度截然不同的兩層.

(2)在顆粒毛細效應能夠發生的管徑范圍內,存在一個對應于顆粒最終毛細上升高度最大值的臨界管徑,當管徑小于臨界管徑時,顆粒最終毛細上升高度隨管徑的增加而增加,當管徑大于臨界管徑時,顆粒最終毛細上升高度隨管徑的增大而減小;臨界管徑受到容器寬度和振動參數的影響.

(3)臨界管徑隨著容器寬度的增加而有所增大;增加振幅和適當增大頻率能夠有效促進臨界管徑的增大,從而將更多的顆粒輸運到更高的高度,這對于利用顆粒毛細效應實現顆粒物料的逆重力輸運具有重要意義.