顆粒物質內自旋小球運動行為的數值模擬研究?

程琦 冉憲文 劉蘋 湯文輝 Raphael Blumenfeld2)

1)(國防科學技術大學理學院工程物理研究所,長沙 410073)

2)(劍橋大學卡文迪許實驗室,劍橋 CB30HE,英國)

1 引 言

許多工業中生產過程,如產品分揀、原材料攪拌[1?4]以及產品運輸[5,6]等,都涉及到顆粒物質系統動力學.同時,顆粒物質系統動力學也是模擬不同特性顆粒分離[7?9]以及顆粒斑圖[10]形成的基礎.顆粒物質系統動力學行為非常復雜,涉及多種未知因素,科學研究過程中必須對具體問題進行簡化,研究顆粒物質體系內較大物體的運動及對體系動力學行為的影響則是許多科技工作者選擇的出發點.該動力學過程與許多現實過程密切相關,如巴西果效應問題[7?9]、對顆粒物質的沖擊成坑以及侵徹問題[11?15]、地下管道的抬升[16]以及沙漠生物的運動等[17?19].

近年來,如何解析沙漠蜥蜴在沙漠上奔跑、潛伏以及沙子中的游動問題得到相關科學家的關注[17?19],并對該問題進行了較多的研究,尤其是顆粒物質中“犁地”和顆粒物質內水平驅動物體的上升問題[20,21],以及顆粒物質內深埋物體運動所引起的浮力等問題[1,3,22?26].相關結果表明,顆粒物質內物體所受到的升力和阻力與物體的運動速度以及所處的深度密切相關,這些結果的得出對于解析沙漠蜥蜴潛入沙子內部以及如何在沙子內部的游動提供了相關依據,但這些結果依然無法解釋大耳沙蜥蜴(Phrynocephalus mystaceus)的運動行為.與其他沙蜥屬的蜥蜴不同,大耳沙蜥是通過左右橫向擺動和轉動身體,激發身體附近沙子的運動形成周圍沙子的“液化”,進而實現潛入沙子內部.針對左右擺動問題,我們以球形物體為研究對象,采用離散元仿真進行了具體研究,獲得了影響球體升/降運動的相圖,并給出了橫向振動頻率和振幅相關的理論模型[27].然而,對于轉動問題,尚缺乏針對性的研究.

本文采用美國圣地亞實驗室開發的三維離散元程序LIGGGHTS[28],對顆粒物質內旋轉的球形物體所導致的升/降問題和平動問題進行了研究,具體研究了角速度、表面摩擦系數等對球形物體運動的影響.

2 模擬設置

2.1 顆粒間的相互作用模型

模擬中顆粒之間的相互作用采用經典的Hertz模型,也即軟球接觸模型[29],如圖1所示.Hertz模型將相互接觸的兩個顆粒之間的相互作用簡化為顆粒之間的切向力(即兩個顆粒中心的連線)和法向力(即在接觸點垂直于切線方向),采用彈簧來模擬相互作用的可恢復部分,采用黏壺來模擬相互作用的耗散特性.

圖1 顆粒間相互作用的Hertz模型 (a)法向相互作用;(b)切向相互作用Fig.1.The Hertz model for the interactions of particles:(a)The normal interaction;(b)the tangential interaction.

根據Hertz模型和彈簧阻尼模型[30],顆粒物質中顆粒i和顆粒j相互接觸時,顆粒i的受力可分為法向力Fnij和切向力Ftij兩個部分,具體表達式為:

其中,nij是由j的球心指向i的球心的單位向量,μij是顆粒i和顆粒j之間的摩擦系數,knij和ktij分別為法向和切向的彈性系數,γnij和γtij分別為法向和切向的阻尼系數,δnij和unij分別是法向重疊量和法向相對速度,δtij和utij分別是切向滑移和切向相對速度.需要指出的是,utij不一定與切方向平行.knij和ktij與楊氏模量和泊松系數相關,而γnij和γtij不僅與楊氏模量和泊松系數相關,還與恢復系數相關.需要明確的是,法向重疊量δnij是不考慮顆粒變形時兩個顆粒之間的重疊量,如圖2所示.切向滑移的定義是:一個時間步長內顆粒i和顆粒j的切向相對位移量.具體計算方法是:utij?dt,然后取切方向的投影得到δtij.

圖2 兩個顆粒的法向重疊量Fig.2.The normal overlapped quantity of two particles.

2.2 模擬參數設置

模擬對象為一個深埋入顆粒系統中的球形物體,通過球形物體的自主旋轉運動來模擬大耳沙蜥身體的轉動狀態.整個模擬空間或模擬容器為立方體,具體尺寸為0.5 m×1 m×1.3 m,相關尺寸參數和模擬參數如表1所列.為了減少邊界對整個系統模擬的影響,我們在X和Y方向設置周期性邊界條件,在Z方向設置固壁邊界條件.

表1 模擬中使用的相關參數Table 1.Related args used in the simulation.

對于多數顆粒系統而言,其內顆粒的分布是完全隨機無規則的,因此無法像生成晶體一樣,按照其排列規則生成顆粒系統.在模擬前,首先采用“落雨法”[31]生成基底顆粒,即從容器上方隨機生成顆粒然后落下,在容器內形成基底顆粒,待生成足夠多的基底顆粒后則停止生成;然后整個顆粒系統在重力的作用下弛豫足夠長的時間,以確保整個系統的總動能穩定并在較低的數值范圍內.而后在容器口生成小球,讓小球以固定速度運動到顆粒系統內基底顆粒中央,帶有小球的顆粒系統再在重力的作用下弛豫,直到系統總動能穩定并在較低的數值范圍內,此時,所生成的模擬構型就是我們要進行系統研究的顆粒系統,如圖3所示.現在每個顆粒的位置就是數值模擬顆粒動力學前的初始位置,并且每個顆粒均接近于靜止狀態.

圖3 重力弛豫后的顆粒系統Fig.3.The granular system after gravity-relaxation.

在獲得模擬顆粒系統的初始狀態后,主要通過數值模擬研究小球的自轉角速度ω以及小球與基底顆粒之間摩擦系數μ對小球在顆粒系統內運動的影響.小球的自轉軸為X軸,其自轉方向如圖4所示.模擬中,小球與基底顆粒間的摩擦系數考慮0.4,0.45,0.5,0.6四種情況;自轉角速度取10,20,40,80,160,320,640,1280 rad/s.

圖4 小球自身的旋轉運動Fig.4.The rotation of the embeded sphere.

3 模擬結果

模擬前,小球埋入的深度約為0.6 m.以小球的初始位置為坐標原點建立笛卡爾直角坐標系,其中X軸方向為小球的自轉方向,Z軸方向為重力加速度反方向,Y軸方向滿足右手系.模擬中發現小球除自身的旋轉運動外,還會有Z方向的上升運動以及X和Y方向的平動.

3.1 Z方向的運動

圖5中給出了50 s內小球在顆粒系統內沿Z方向上的運動情況:1)圖5(a)中,小球與基底顆粒間摩擦系數為μ=0.4,小于基底顆粒之間的摩擦系數μ=0.45,此時不論小球自轉的角度度多大,基本上都沒有產生小球的上升運動,僅在角速度為1280 rad/s時,小球才有少許的向上波動;2)圖5(b)中,小球與基底顆粒間摩擦系數為μ=0.45,與基底顆粒間的摩擦系數相同,此時開始出現上升趨勢,但在50 s的時間內上升距離有限,并且不同角速度之間上升的速度開始出現差異,基本上是角速度越大上升越明顯;3)圖5(c)中,小球與基底顆粒間的摩擦系數為,大于基底顆粒間的摩擦系數,此時小球的上升趨勢十分明顯,在50 s的時間內,角速度為160,320,640,1280 rad/s的四種情況中,小球均上升到了顆粒系統頂部,成功地“飄”在顆粒物質上;4)圖5(d)中,小球與基底顆粒間的摩擦系數為μ=0.6,大于基底顆粒間的摩擦系數,此時小球的上升趨勢更為明顯,在50 s的時間內,角速度為40,80,160,320,640,1280 rad/s六種情況中的小球都成功上升到了顆粒系統頂部,比圖5(c)多出了兩種狀態.

圖5(c)和圖5(d)中,當角速度為1280 rad/s時,小球的上升速度反而比其角速度為640 rad/s時慢.其原因可能是小球的角速度超過了臨界點,使得Ftij的計算變為所導致,但具體的原因還需要更深入的探索.

綜合比較圖5的結果可以看出,在小球與基底顆粒的摩擦系數μ小于基底顆粒的摩擦系數時,無論小球自身旋轉速度多少,小球都基本上“懸浮”于顆粒系統中;當μ與基底顆粒間的摩擦系數相同時,小球開始出現上升運動,并且隨著μ的增加,上升趨勢更為明顯;對比μ相同、但小球自身旋轉速度不同的情況,自轉速度越大,小球上升的趨勢就越為明顯.

圖5 Z方向上小球的運動Fig.5.The motion of sphere in the Z direction.

3.2 X方向上的運動

小球在X方向上的運動如圖6所示,可以看出X方向上的運動沒有明顯的規律.當小球與基底顆粒間的摩擦系數小于基底顆粒間的摩擦系數時,除角速度為1280 rad/s的情況外,小球位置在X方向上基本不變,即使是角速度1280 rad/s時,小球在X方向上的位移也很小,遠小于小球的半徑0.06 m;當小球與基底顆粒間的摩擦系數大于基底顆粒間的摩擦系數時,小球在X方向上的位移變得相對明顯,但隨著小球與基底顆粒間摩擦系數的增加,X方向上的位移并沒有明顯的增加,而且X方向上的位移均小于小球半徑0.06 m,這說明小球X方向上的運動與X方向上的自轉關系不大.

3.3 Y方向上的運動

小球在Y方向上的運動如圖7所示,與X方向的運動相比,顯然Y方向的運動要明顯得多.當μ小于基底顆粒間的摩擦系數時,如圖7(a)所示,小球位置在整個50 s的模擬時間內在Y方向上基本保持不變,只有角速度為1280 rad/s時,小球才在44 s后開始產生Y方向運動,但所產生的位移也小于小球的半徑0.06 m.當μ等于基底顆粒間的摩擦系數時,如圖7(b)所示,小球在Y方向上的位移開始變得明顯,并且呈現明顯的分叉狀態,即小于160 rad/s的自轉在Y方向上引起Y軸正向位移;大于320 rad/s的自轉在Y方向上引起Y軸負向位移.當μ大于基底顆粒間的摩擦系數時,如圖7(c)和圖7(d)所示,小球自轉所引起Y方向的位移更為明顯,可達到小球半徑的4—6倍;對比圖7(c)和圖7(d)可以看出,小球自轉角速度以及小球與基底顆粒之間的摩擦系數都會對Y位移產生影響,二者的增加都會帶來Y方向位移的增加;特別是對于角速度160 rad/s的情況,當小球與基底顆粒之間的摩擦系數為μ=0.5時,小球的運動趨勢是Y軸正向,而當μ=0.6時,小球的運動趨勢則為Y軸負向;當小球上升到顆粒系統表面時,由于小球自轉的影響,則在Y軸負向產生明顯的位移,此種情況與沙漠表面車輪的運動類似.

圖6 X方向上小球的運動Fig.6.The motion of sphere in the X direction.

圖7 Y方向上小球的運動Fig.7.The motion of sphere in the Y direction.

4 結 論

根據數值模擬和討論結果,對于深埋入顆粒物質中物體的旋轉運動,可以得出以下結論:

1)深埋入顆粒物質內的物體,水平旋轉運動能夠引起豎直以及水平方向上相對明顯的運動,但在角速度方向上所引起運動位移有限,小于物體自身的尺度;

2)物體與基底顆粒間的摩擦系數對物體在顆粒物質內部的運動軌跡有顯著影響,只有當物體與基底顆粒間的摩擦系數大于顆粒物質自身基底顆粒間的摩擦系數時,物體才會獲得顯著的運動,這也許就是大多數沙漠生物外表披掛鱗片的原因;

3)在物體與基底顆粒間的摩擦系數大于顆粒物質自身基底顆粒間的摩擦系數的前提下,深埋入顆粒物質中的物體自身的旋轉角速度對其運動軌跡也會產生很大的影響,旋轉角速度越快所引起物體的運動就越明顯.

[1]Soller R,Koehler S A 2006Phys.Rev.E.Stat.Nonlin.Soft Matter Phys.74 021305

[2]Guillard F,Forterre Y,Pouliquen O 2013Phys.Rev.Lett.110 138303

[3]Guillard F,Forterre Y,Pouliquen O 2014Phys.Fluids26 043301

[4]Guillard F,Forterre Y,Pouliquen O 2015Phys.Rev.E:Stat.Nonlin.Soft Matter Phys.91 022201

[5]Knight J B,Ehrichs E E,Kuperman V Y,Flint J K,Jaeger H M,Nagel S R 1996Phys.Rev.E54 5726

[6]Eshuis P,van der Meer D,Alam M,van Gerner H J,van der Weele K,Lohse D 2010Phys.Rev.Lett.104 038001

[7]Shinbrot T,Muzzio F J 1998Phys.Rev.Lett.81 4365

[8]M?bius M E,Lauderdale B E,Nagel S R,Jaeger H M 2001Nature414 270

[9]Shinbrot T 2004Nature429 352

[10]Cheuk C Y,White D J,Bolton M D 2008J.Geotech.Geoenviron.Eng.134 154

[11]Metcalfe G,Shattuck M 1996Physica A:Stat.Mech.Appl.233 709

[12]Wood R J K,Wheeler D W 1998Wear220 95

[13]Lohse D,Rauhe R,Bergmann R,van der Meer D 2004Nature432 689

[14]Lohse D,Bergmann R,Mikkelsen R,Zeilstra C,van der Meer D,Versluis M,van der Weele K,van der Hoef M,Kuipers H 2004Phys.Rev.Lett.93 198003

[15]Katsuragi H,Durian D J 2007Nat.Phys.3 420

[16]Clark A H,Petersen A J,Kondic L,Behringer R P 2015Phys.Rev.Lett.114 144502

[17]Dowling K J 1996Ph.D.Dissertation(Pittsburgh:Carnegie Mellon University)

[18]Maladen R D,Ding Y,Li C,Goldman D I 2009Science325 314

[19]Maladen R D,Ding Y,Umbanhowar P B,Kamor A,Goldman D I 2011J.R.Soc.Interface8 1332

[20]Costantino D J,Bartell J,Scheidler K,Schiffer P 2011Phys.Rev.E:Stat.Nonlin.Soft Matter Phys.83 011305

[21]Percier B,Manneville S,Mcelwaine J N,Morris S W,Taberlet N 2011Phys.Rev.E:Stat.Nonlin.Soft Matter Phys.84 051302

[22]Wassgren C R,Cordova J A,Zenit R,Karion A 2003Phys.Fluids15 3318

[23]Chehata D,Zenit R,Wassgren C R 2003Phys.Fluids15 1622

[24]Ding Y,Gravish N,Goldman D I 2011Phys.Rev.Lett.106 028001

[25]Maladen R D,Umbanhowar P B,Ding Y,Masse A,Doldman D I 2011Robotics and Automation(ICRA),2011 IEEE International Conference on.IEEEShanghai,China,May 9–13,2011 p1398

[26]Potiguar F Q,Ding Y 2013Phys.Rev.E88 012204

[27]Huang L,Ran X,Blumenfeld R 2016Phys.Rev.E94 062906

[28]Kloss C,Goniva C,Hager A,Amberger S,Pirker S 2012Prog.Comput.Fluid Dyn.12 140

[29]Gran model hertz model,Christoph K http://www cfdem com/media/DEM/docu/gran_model_hertz html[2017-04-03]

[30]Sun Q C,Wang G Q 2009Introduction to Mechanics of Granular Materials(Beijing:Science Press)p31(in Chinese)[孫其誠,王光謙 2009顆粒物質力學導論(北京:科學出版社)第31頁]

[31]Huang D C,Feng Y D,Xie W M,Lu M,Wu H P,Hu F L,Deng K M 2012Acta Phys.Sin.61 124501(in Chinese)[黃德財,馮耀東,解為梅,陸明,吳海平,胡鳳蘭,鄧開明2012物理學報61 124501]