三維圓柱型顆粒堆坍塌問題的全相態數值模擬1)

陳福振 李亞雄 史騰達 嚴紅

(西北工業大學動力與能源學院,西安 710072)

(西北工業大學太倉長三角研究院,江蘇太倉 215400)

引言

顆粒堆坍塌問題是自然界和工業工程中最為常見的一種物理現象,如倉庫堆積物料突然塌方現象、砂粒輸運過程中的傾翻現象、地球物理中的巖石崩塌、雪崩、海底塌方等等,給人類的生命和財產安全造成嚴重的影響.盡管這種現象無處不在,但即使是涉及顆粒介質的靜態問題,例如確定靜止沙堆底部中心的壓力問題等,仍然存在一定的爭議.而運動的顆粒介質則呈現出更加豐富和復雜的現象,如最為典型的顆粒介質在不同體積分數和載荷作用下會表現出類固體行為、類液體行為、類氣體以及完全離散的慣性行為等[1-5].因此,建立可有效描述這些不同行為的顆粒介質運動宏觀預測理論一直是物理學領域的一個重要目標,其中Science也將顆粒介質力學行為的通用理論模型定義為人類尚未解決的125 個科學難題之一.本文就主要針對顆粒堆坍塌過程中所表現出的復雜相態特性,建立一種描述顆粒介質全部相態的物理模型和數值模擬方法,為預測顆粒堆坍塌規律提供一條新的有效的途徑.

在顆粒堆坍塌實驗方面,Lube 等[6]研究了三維軸對稱顆粒堆坍塌過程的三種流動狀態,分析了不同高寬比下的運動規律.在此基礎上又補充開展了二維不同材料下的坍塌過程實驗[7].Lajeunesse 等[8-9]采用矩形通道和半圓形管道兩種不同的裝置,實驗獲得了二維和軸對稱顆粒坍塌運動過程,特別是獲得了顆粒堆內部流動結構.Roche 等[10]開展了細顆粒(約75 μm)和粗顆粒(約330 μm)情況下,最初與空氣流化形成的軸對稱顆粒堆釋放產生重力顆粒流鋪展的實驗結果.Artoni 等[11]研究了濕顆粒堆的坍塌特性,旨在收集濕顆粒介質中的坍塌觸發和休止現象的實驗數據,獲得了不同粒徑、不同液體表面張力和不同液體量的玻璃微珠樣品的最終沉積形態和鋪展動力學過程.以上為典型的在水平表面上開展的顆粒堆坍塌實驗.除此之外,文獻[12-16]開展了在斜面上進行顆粒堆坍塌實驗的工作.雖然實驗可以獲得顆粒堆坍塌過程典型時刻結果,但是費時費力,同時顆粒運動的細節無法全部獲取,對于揭示顆粒運動的機理支撐度有限.因此,基于顆粒介質運動理論,采用數值模擬的方法對顆粒堆坍塌過程進行深入研究,對于全面掌握顆粒堆坍塌運動規律,預測顆粒介質運動的不同現象具有重要意義.

在顆粒堆坍塌數值模擬方面,主要包括兩大類方法,一類是以離散元(DEM)為代表的顆粒軌道追蹤方法,通過計算作用在每個離散單元上的作用力給出其單獨運動的軌跡.文獻[17-19]采用DEM 方法對平面上顆粒堆的坍塌過程進行了數值模擬.孫其誠和王光謙[20]對重力作用下12000 個球心共面的二維等徑顆粒靜態堆積進行了離散動力學模擬,從力鏈角度揭示顆粒靜態和動態性質.成浩等[21]基于DEM 方法對松散體滑動堆積特性及影響因素進行了分析.Zhang 等[22]采用DEM 方法對三維顆粒堆坍塌過程進行了數值模擬,研究了顆粒堆厚度對坍塌過程的影響.DEM 方法對于揭示顆粒介質的運動機理、擬合獲得顆粒介質的宏觀統計特征、小規模的工程應用非常有效,但是該方法的計算消耗對于模擬工業系統中存在較大數量的顆粒來說往往遙不可及;該方法所使用的參數和模型是單顆粒層次的,往往與可開展的宏觀實驗不匹配;計算的時間步長與顆粒的硬度息息相關,往往很受限制.雖然有學者[23-24]開發了DEM 的粗粒化處理技術,但存在著條件苛刻、放大不足、顆粒變硬需要進一步降低時間步長而增大計算時長等問題.

另一類是基于宏觀連續介質力學的數值模擬方法,如求解顆粒彈塑性本構模型的有限元方法(finite element method,FEM)[25]和有限體積方法(finite volume method,FVM)[26],求解基于流變學的黏塑性本構模型的FEM[27-29]和FVM[30-32]等,此類方法的特點是基于網格離散求解數理模型,已成功用于顆粒堆坍塌、剪切造粒、化工流化床、顆粒管道輸送等問題中,但此類方法存在的最大問題是必須依賴網格求解帶來的缺陷,如有限元在計算類固態和小變形的類液態時是合理的,但對于顆粒介質宏觀大變形問題或者類氣態問題往往會出現網格的扭曲和纏繞,計算終止;而采用有限體積法解決顆粒類氣態問題和類液態問題較為合適,但類固態本構模型無法求解,同時該方法無法獲得顆粒的運動軌跡,易產生偽擴散,不易加入顆粒蒸發、燃燒等物理化學模型.

為了克服網格方法在求解顆粒宏觀連續介質力學模型上的不足,有學者嘗試采用粒子類方法進行模擬,如文獻[33-38]均采用物質點法(material point method,MPM)根據土壤的連續介質力學模型計算了類土體的顆粒介質運動問題.此類方法一方面較難處理顆粒的類氣態問題,甚至是體積分數更小的慣性態,這時需要將氣相應力強制置零[38],或者引入其他模型來處理干燥顆粒材料的離散運動;另一方面物質點法必須使用背景網格,依賴于背景網格求解動量方程,不可避免地存在網格重分、大范圍布置背景網格、網格與物質點之間需要不斷反復插值的問題.

光滑粒子流體動力學方法(smoothed particle hydrodynamics,SPH)[39-40]作為另外一種完全拉格朗日粒子方法,對離散顆粒進行模擬表征具有很大優勢.SPH 不僅適合于模擬處于類固態和類液態的顆粒相體積分數幾乎保持恒定的問題[41-46],同時對于類氣態問題也非常適合模擬.Chen 等[47-51]前期就在傳統SPH 方法的基礎上建立了SPH 粒子與類氣態離散顆粒間的一一對應關系,將SPH 改造成適于分散性顆粒相求解的光滑離散顆粒流體動力學方法(smoothed discrete particle hydrodynamics,SDPH),耦合FVM 求解連續氣體相,成功應用于化工流化床[48]、風沙躍移[49]、氣-粒傳熱[50]、空氣燃料拋灑[51]等領域中,解決了氣流載體作用下離散顆粒的類氣態數值模擬問題.Chen 和Yan[52-53]近些年嘗試將SDPH 方法應用于顆粒其他相態的數值模擬,取得了一定的進展.

本文即在此基礎上,進一步探索建立顆粒介質的全部相態物理模型.從顆粒表現出不同相態的物理機理出發,將描述顆粒介質的彈塑性理論、黏塑性理論、顆粒動理學理論以及單顆粒輸運理論有效結合起來,通過確定不同相態之間的過渡關系和轉化準則,建立描述顆粒介質經歷全部相態的耦合模型理論,不同相態之間不僅可以共存,同時可以正向和反向轉化.然后,采用SDPH 方法模擬類固態、類液態和類氣態,同時耦合DEM 模擬顆粒離散慣性運動狀態,建立求解顆粒介質全部相態的數值模擬新方法.新方法針對不同的相態采用不同的與之相匹配的粒子方法,既保證了對各個相態的準確描述和動態再現,又降低了計算量.最后,采用新的理論模型和數值方法,對不同長徑比條件下的三維圓柱型顆粒堆坍塌問題進行了數值模擬,一方面驗證了模型和方法的準確性和適用性,另一方面探明了影響顆粒堆坍塌特性的因素和機理,為顆粒介質運動問題的解決提供理論和數據支撐.

1 顆粒介質全相態物理模型

顆粒介質作為一種由介觀離散顆粒組成的宏觀無定型態物質,根據顆粒的體積分數不同和載荷作用條件不同會表現出一些不同的行為特征,例如顆粒體在堆積狀態時,整體表現出類似固體結構的行為,運動非常緩慢,在屈服的過程中會形成剪切帶,應力狀態也表現出率無關狀態,這種狀態為顆粒介質的類固態,稱之為顆粒固相;當顆粒堆積過程中承受的等效剪應力與等效體應力的比值超過一定閾值時,顆粒介質發生塑性流動,產生類似于液體之間的黏性剪切作用效果,這種狀態為顆粒介質的類液態,稱之為顆粒液相;當顆粒之間運動的速度梯度繼續加大,顆粒體積分數減小,不再等價于不可壓縮狀態時,顆粒與顆粒之間的相互作用力也不再遵循多顆粒接觸假設,而主要以頻繁的二體碰撞為主,碰撞速度相互獨立,應力表現為剪切率的二次函數,這時的顆粒介質表現出類似氣體運動的行為,這種狀態為顆粒介質的類氣態,稱之為顆粒氣相.

以上是目前國內外對于顆粒介質存在的三種相態的描述,作者通過分析發現當顆粒的體積分數繼續減小,顆粒間的二體碰撞假設也不再滿足,顆粒間碰撞的概率已經非常微小,顆粒主要以受到的體力作用和外部流場作用為主,這時從相對尺度上來說顆粒屬于介觀尺度范圍,不再遵循擬流體的連續介質力學定律,因此不適用于以上三種狀態描述,將這種狀態稱為離散慣性相,遵循質點運動定律.由此,顆粒介質從濃密到稀疏、從連續到離散、體積分數從1 至0 的全部相態可全覆蓋,如圖1 所示.

圖1 顆粒介質全相態定義及物理模型示意圖Fig.1 Definition and diagram of the physical model of all phases of granular media

本文在建立顆粒介質全相態理論之前,根據顆粒稀疏程度的不同將全相態劃分為三個區域,將顆粒固相和液相區域統稱為濃密顆粒介質區(φl,min≤φp≤φs,max,φp為顆粒體積分數,φs,max為顆粒固相最大體積分數,φl,min為顆粒液相最小體積分數,通常 φl,min≈φs,max);將顆粒液相和氣相之間的過渡相和氣相區域稱之為稀疏顆粒介質區(φg,min≤φp≤φl,min);將離散慣性相區域稱之為超稀疏顆粒介質區(φp<φg,min).劃分區域的目的是分區域建模,使復雜的多相態建模能夠實現一定程度的簡化.建立顆粒介質的全部相態本構理論就需要根據其所處的不同相態區域進行建模,同時建立不同相態之間的轉變原則:

(1) 對于濃密顆粒介質區域狀態采用彈-黏-塑性本構理論描述,具體針對濃密顆粒介質固相采用彈塑性理論描述,對于濃密顆粒介質液相采用流變學理論描述;

(2)對于稀疏顆粒介質區域狀態采用顆粒動理學理論和摩擦動力學理論描述,具體針對稀疏顆粒介質的液氣過渡相采用顆粒動理學和摩擦動力學相結合的方式描述,對于稀疏顆粒介質氣相采用顆粒動理學描述;

(3)對于超稀疏顆粒介質區域狀態采用質點動力學理論描述,具體就是指針對超稀疏顆粒介質的慣性相采用質點動力學理論進行描述.

下面就分別介紹這些不同的理論,以及這些理論之間的轉變原則.

1.1 描述濃密顆粒介質運動狀態的彈-黏-塑性本構模型(φl,min ≤φp ≤φs,max )

對于濃密顆粒介質來說,首先介紹其運動控制方程,包括質量守恒方程和動量守恒方程兩個,如下

本文應用了愛因斯坦求和約定的指數表示法,α和 β 分別表示笛卡爾坐標系下的1,2 和3 三個分量,ρ 為材料的密度,v為速度,σαβ為材料的全應力張量分量,通常由兩部分組成:各向同性靜水壓力p和應力偏張量s

式中,δαβ是克羅內克函數,當α=β時δαβ=1,當α≠β時δαβ=0.靜水壓力p直接采用本構方程中的應力計算得到

全應力張量的具體形式將根據其所處的狀態不同而發生改變,具體見以下章節介紹.fα為其他外力,如本文中所需要考慮的重力.d/dt為全導數.

有效耦合現有的彈塑性理論和基于流變學的黏塑性理論,獲得基于稠密顆粒運動彈-黏-塑性理論的顆粒固液相計算模型[52].假設準靜態和流動狀態之間的屈服準則用Drucker-Prager 屈服準則表示,則可以建立從完全彈性到彈性-微黏塑性再到完全彈性-黏塑性的過渡過程,如圖2 所示.

圖2 完全彈性到彈性-微小黏塑性再到完全彈性-黏塑性的過渡過程Fig.2 The transition from complete elastic to elastic-micro viscoplastic and then to complete elastic-viscoplastic

1.1.1 完全彈性下的顆粒介質應力-應變關系

1.1.2 彈性-微小黏塑性狀態下的顆粒介質應力-應變關系

隨著顆粒受力的增加,變形逐漸增大.彈性引起的剪應力逐漸增大,首先達到.此時,根據顆粒流變學理論,顆粒介質開始流動,和 μ(I)逐漸增加.雖然發生了流動,但流動速度非常小,尚未達到塑性屈服狀態.在此階段,只有剪應力逐漸增大.

彈性-微小黏塑性狀態下顆粒介質的法向應力仍按線彈性本構模型計算,如式(5),剪應力按流變學模型計算,即

1.1.3 彈性-完全黏塑性狀態下的顆粒介質應力-應變關系

以及其中的塑性乘子變化率公式

1.2 描述稀疏顆粒介質區域的顆粒動理學理論和摩擦動力學理論(φg,min ≤φp<φl,min )

1.2.1 描述顆粒氣相的動理學理論(KTGF)

對于稀疏顆粒流區域中的顆粒氣相,采用顆粒動理學理論(簡稱KTGF)[54-55]進行描述,一共包括質量守恒、動量守恒、擬溫度守恒等三個方程,質量守恒方程和動量守恒方程如式(1)和式(2),擬溫度守恒方程式為顆粒動理學理論所特有,如下

式中,ds為顆粒的直徑,ess為顆粒之間的碰撞恢復系數.g0為顆粒的徑向恢復系數

式中,φs,max為顆粒介質可達到的最大體積分數值.

1.2.2 描述顆粒氣相-液相過渡區的摩擦動力學理論

摩擦動力學主要用來描述處于長程接觸以及存在不可忽略的摩擦應力的顆粒.一些學者研究表明,摩擦動力學可以用于體積分數較高的顆粒[56-57].在稠密顆粒流和稀疏顆粒流之間,必然存在一個從長程接觸到瞬態碰撞接觸的過渡階段.這里,摩擦動力學被用來彌補顆粒氣相和液相之間存在的過渡缺陷問題.Johnson 等[58-59]提出的半經驗模型用來計算由摩擦引起的法向應力,如下所示

式中,Fr,a和b是材料的經驗常數.摩擦黏度的計算最早由Schaeffer[60]在研究筒倉顆粒物料在重力作用下從錐形出口流動的過程中,假定服從理想剛塑性本構理論,獲得了摩擦剪應力與正應力之間存在以下關系

式中,? 為內摩擦角.通過分析可知[53,59],摩擦動力學可以被認為是 μ(I) 流變學[61-62]在慣性數或者說剪切應變率無限大時的極限情況.

1.3 描述超稀疏顆粒流區域的質點動力學理論(φp<φg,min )

離散質點動力學是指針對具有一定質量但幾何尺寸大小可以忽略的物體,采用牛頓運動定律進行描述的動力學過程

式中,右端項Fdrag為顆粒所受到的曳力.Fcol為顆粒之間的碰撞作用力,Fg為重力.本文重點是考慮顆粒間的碰撞相互作用力(見第3 節).

1.4 濃密顆粒流與稀疏顆粒流兩個狀態之間的轉變原則

KTGF 結合摩擦動力學模型[56-57]用于計算顆粒液相和氣相之間的轉換.一些學者通過研究發現,摩擦動力學可以用于具有較高體積分數的顆粒,但當體積分數增加到某個值時存在不確定性.然而,將摩擦動力學與彈-黏-塑性模型相結合可以克服這一缺點.在轉變過渡的界面處應力應保持守恒,彈性-黏塑性模型計算的法向應力和剪應力應與摩擦動力學計算所得的法向應力和剪應力相同.顆粒應力計算公式表示如下

式中,pEVP和 μEVP分別為采用彈-黏-塑性本構模型計算獲得的正應力和剪應力,pKTGF和pfriction分別是顆粒壓力pp中的動理學部分和摩擦部分,μKTGF和μfriction分別是顆粒剪切黏度 μp的動理學部分和摩擦部分.φg,max是摩擦應力開始逐漸增加時的顆粒體積分數.小于 φg,max值時,通過已知資料發現未觀察到顆粒間的摩擦行為,因此假設,當均勻分布的粒子不再接觸時,摩擦相互作用不再發生,主要由碰撞相互作用產生相互間的應力.

顆粒液相和氣相之間的過渡如圖3 所示.從顆粒液相和顆粒氣相之間的過渡相開始,顆粒擬溫度從零開始升高,表明顆粒之間的碰撞逐漸加劇.從過渡相到氣相轉變的界面上,顆粒之間的摩擦完全消失,轉變為完全二體碰撞.從上述計算公式可以看出,用這種方法建立的過渡態可以有效地連接液相和氣相,實現平穩過渡,符合物理定律.

圖3 濃密顆粒流與稀疏顆粒流之間的過渡轉化Fig.3 Transition between dense granular flow and dilute granular flow

1.5 稀疏顆粒流與超稀疏顆粒流兩個狀態之間的轉變原則

對于顆粒體積分數小到一定數值后(φg,min),顆粒的宏觀流動特征時間尺度明顯小于微觀碰撞特征時間尺度,對顆粒類氣態的宏觀連續描述不再成立,這時轉化為顆粒質點進行追蹤.由于在顆粒氣相計算的過程中,進行了宏觀擬流體假設,一個單元表征了在空間該位置處顆粒的統計平均信息,如顆粒的有效密度、顆粒的均值粒徑、顆粒的均值速度等,這時在轉化的過程中需要保證轉化的質量、動量和能量的守恒.后面的章節中會介紹對于稀疏顆粒流采用的是拉格朗日粒子方法SDPH 進行離散求解,而超稀疏顆粒流的離散質點動力學也是采用粒子法DEM 進行模擬,因此轉化的過程較為自然,可以保證物理量的守恒,同時為了更加貼近實際,還可以結合粒子分裂算法,使兩者在空間位置上也對應起來.具體粒子間的轉化方法見3.3 節.

2 顆粒介質全相態數值模擬方法

第2 節闡述了顆粒介質全相態理論,要對顆粒介質全相態進行數值模擬除了理論模型之外,還需要引入合適的數值模擬方法對模型進行離散求解.本文對于顆粒介質全相態中的類固態、類液態和類氣態三種狀態采用SDPH 方法進行模擬,因為這三種狀態的本構模型均是基于宏觀連續介質力學所建立的,與SDPH 的理論基礎一致.當散體顆粒分散到一定程度即顆粒相的體積分數小到一定程度時,采用2.3 節的離散顆粒動力學進行描述,引入描述單一顆粒行為的離散單元法進行計算.

2.1 基于顆粒連續介質力學的光滑離散顆粒流體動力學方法

為明確SPH 粒子與顆粒之間的對應關系,Chen 等[47-51]前期已經從顆粒動理學角度出發,將顆粒相視為擬流體,擬流體區域采用SPH 方法離散求解,同時將傳統SPH 方法改造成了適用于離散顆粒相求解的SDPH 方法,本文在此基礎上繼續拓展SDPH 方法的應用范圍,對于顆粒固相和液相同樣采用SDPH 離散闡述,與顆粒氣相的SDPH 描述保持一致.

2.1.1 光滑離散顆粒流體動力學方法

在SDPH 方法中,實際顆粒的質量、速度、壓力、位置、粒徑分布、體積分數等均在計算粒子上有體現.同時,顆粒氣相的擬溫度值也作為計算粒子的一個參量.計算粒子與實際顆粒之間的參量對應關系如下

其中,p為顆粒下標,φp為顆粒體積分數,ρp為顆粒實際密度.顆粒的有效密度可進一步推導獲得

式中,n為空間內真實顆粒的總數目,Vp為每個實際顆粒所占據的空間體積,mp為每個實際顆粒的質量,V0表征了所有實際顆粒占據的空間總體積.式(23)表示計算粒子的密度等于真實顆粒的有效密度.單個計算粒子的質量和體積表征了空間中某些真實顆粒的總質量和體積.由每個計算粒子表示的真實顆粒數通過每個計算粒子的質量與每個真實顆粒的質量之比獲得.計算粒子的速度、壓力和擬溫度取真實顆粒群相應參數的平均值.

從以上對應關系可以看出,當使用SDPH 方法進行計算時,每個計算粒子可以代表一定數量的真實顆粒群,計算粒子的數量可以根據計算精度的要求確定.對于某些精度要求不太嚴格的問題,計算量將大大減少,計算效率將顯著提高.

2.1.2 顆粒固相和液相的SDPH 離散公式

在SPH 方法中,流體域被一系列粒子離散和求解.函數f(r) 及其空間導數 ? ·f(r) 在粒子i處的估計值為

其中,m,ρ,r分別為物質的質量、密度和空間位置矢量,為光滑函數或核函數,本文采用三次樣條核函數,h為光滑函數W的影響域或支持域的長度.

式(1)和式(2)中顆粒固相和液相的質量和動量守恒方程采用SDPH 方法離散,如下所示

其中,總應力張量 σαβ同樣采用SDPH 離散,得到

為了滿足零階和一階一致性條件,本文采用了一種基于核及其梯度正規化的SPH 方法的修正形式.核近似的修正形式:為了實現C0一致性,函數f(xi)的修正核近似為[63]

為了滿足動量守恒和C1一致性,本文采用混合核和梯度修正公式[64],這種混合校正是一種改進核梯度的梯度校正和內核修正的組合,它改變了內核本身.梯度修正保證了角動量守恒,并且在內力與密度方程變化一致的情況下,精確計算了線性場的梯度,公式如下

2.1.3 顆粒氣相的SDPH 離散公式

采用SDPH 方法對稀顆粒氣相的擬溫度守恒方程(11)進行離散,得到

顆粒擬溫度梯度 ? θ 采用SDPH 離散獲得

不同于顆粒固相和液相,顆粒氣相的擬溫度變化率 d θ 需要額外進行計算.然后,采用式(33a)獲得新的時刻的顆粒擬溫度值 θ .

2.1.4 基于SDPH 方法的濃密顆粒介質與稀疏顆粒介質相態轉變

采用SDPH 方法不論是求解濃密顆粒相還是稀疏顆粒相,計算粒子與實際顆粒之間的對應關系不會發生改變,只不過在求解的本構模型上存在不同.決定顆粒介質是處于濃密態還是稀疏態的參數是體積分數,體積分數的計算公式為.圖4 展示了轉變策略.如果計算粒子M的體積分數 φp大于φl,min,它表示該粒子正處于顆粒固相和液相.當 φp降低至小于 φl,min,它表示顆粒已經轉變為稀疏顆粒狀態或顆粒氣相.將變換后的計算粒子標記為N.粒子M和N之間的變量關系如圖4 所示.

圖4 基于SDPH 方法的濃密顆粒介質與稀疏顆粒介質相態轉變Fig.4 Phase transition between dense granular media and dilute granular media based on SDPH method

首先要保持粒子的物性參量均不變,包括粒子的位置、速度、密度、能量等,主要區別在粒子間相互作用力上.由于濃密顆粒流的正應力采用彈性定律計算,剪應力采用彈性剪應力與塑性剪應力加和的方式計算,在向稀疏顆粒介質算法轉變時,這兩個作用力也應該保持不變,即由濃密顆粒介質計算的彈性正應力轉化為稀疏顆粒介質的摩擦正應力,見式(17),由式(17)反向計算求得Fr作為不變量,而后根據體積分數 φp的變化更新由摩擦產生正應力的值;對于摩擦剪應力則繼續采用 μ (I) 流變學剪應力公式計算 τ=μ(I)p|γ|/γ,只不過這時的p開始由式(17)計算;對于由長時接觸產生的彈性剪應力則置零;以上是對于SDPH 采用摩擦動力學在過渡區產生的正應力與剪應力的數值計算,保證了轉化的動量的守恒;同時從轉化開始,擬溫度的值由零開始計算,從而由碰撞產生的正應力和剪應力逐漸增大,直到過渡區全部轉化為顆粒動理學模型計算.

同樣地,對于顆粒從稀疏態轉化為濃密態時,粒子變量值應該保持不變,包括粒子的密度、速度、能量、坐標等.兩個狀態的粒子存在的差別主要在內力的計算上.從稀疏態到濃密態過程中,顆粒的擬溫度值逐漸降低,顆粒間的碰撞頻次也逐漸降低,而由摩擦作用產生的正應力和剪切力數值則逐漸增加.在轉化的界面上,由摩擦作用產生的正應力和剪切力與顆粒液相計算的正應力和剪切力相同.在轉化為濃密顆粒態之后,顆粒間的塑性剪切力逐漸降低,彈性剪切力逐漸增加,顆粒速度逐漸降低,體積分數逐漸增加,按照塑性流動法則計算顆粒介質材料的卸載過程,直至恢復到靜止狀態.

2.1.5 基于SDPH 方法的濃密顆粒介質與稀疏顆粒介質間相互作用

當在空間中存在處于濃密狀態的顆粒與處于稀疏狀態顆粒共存時,兩種相態的顆粒之間將有機會產生相互作用.圖5 顯示了SDPH 法中濃密介質和稀疏介質之間的相互作用.因為SDPH 粒子之間的相互作用與否取決于臨近粒子搜索,當被搜索粒子位于搜索粒子的支持域時,被搜索粒子將對搜索粒子產生相互作用.假定處于稀疏狀態的粒子為搜索粒子,處于濃密狀態的粒子為被搜索粒子同時位于搜索粒子的支持域內時,兩者產生類似于二體碰撞的應力,處于濃密狀態的粒子對處于稀疏狀態的粒子有力的貢獻;當處于濃密狀態的粒子為主粒子,處于稀疏狀態的粒子為被搜索粒子同時位于搜索粒子的支持域內時,由于濃密顆粒之間的計算依賴于咬合接觸,而稀疏顆粒與濃密顆粒之間的距離超出了咬合接觸的范圍,因此處于稀疏狀態的粒子對濃密狀態的粒子之間不產生內力作用,而轉化為一種外力施加于動量方程右端項.

圖5 基于SDPH 方法的濃密顆粒介質與稀疏顆粒介質間相互作用Fig.5 Interaction between dense granular media and dilute granular media based on SDPH method

2.2 基于顆粒質點顆粒動力學的離散單元法

采用DEM 對方程(20)中描述超稀顆粒流的微分方程進行離散,并給出以下方程

2.3 SDPH 與DEM 之間算法的耦合

2.3.1 SDPH 與DEM 之間算法的轉化

在稀疏顆粒流SDPH 粒子的體積分數降到一定閾值后(φg,min)時,其不再遵循二體碰撞假設的顆粒動理學模型,因此將SDPH 粒子轉化為DEM 粒子進行計算.圖6 顯示了SDPH 和DEM 之間算法的轉化方案.

圖6 SDPH 和DEM 算法之間的轉化方案Fig.6 Conversion between SDPH and DEM

轉化策略是,將一個SDPH 粒子轉化為一個DEM 顆粒,SDPH 粒子的質量、速度、剛度、位置等參數與轉化后的DEM 的粒子相同,DEM 粒子的密度為實際表征的顆粒的密度,因此根據DEM 粒子的質量和密度便可計算出DEM 在轉化后的粒徑大小,即與每個SDPH 粒子表征具有一定分布的顆粒群一樣,采用該方法轉化之后的DEM 粒子也是代表了一個顆粒群,代表的顆粒群的屬性和參量與SDPH 粒子相同.存在的不同是,DEM 粒子中的顆粒體積分數為1,DEM 粒子的密度只有一個值即真實顆粒的密度.因此,與SDPH 粒子相比,DEM 粒子的尺寸明顯小于SDPH 粒子的尺寸.可以看出,該處的DEM 粒子表征了具有相同粒徑分布的顆粒群,等于是粗粒度的DEM 算法,采用文獻中的算法計算[24].

2.3.2 SDPH 粒子與DEM 顆粒之間作用力傳遞

假如在計算顆粒介質過程中,空間中既有SDPH 粒子,又有DEM 粒子,當他們之間處于相互接觸的狀態時,需要計算他們之間的接觸作用力.采用DEM 粒子之前的接觸力計算方法進行計算.因此,需要將處于接觸狀態的SDPH 粒子轉化為DEM 虛粒子進行計算(等效兩個DEM 粒子接觸計算).圖7 展示了SDPH 粒子與DEM 顆粒之間的接觸力計算規則.SDPH 粒子與DEM 顆粒之間相互作用力包括接觸力Fc,ij=Fcn,i j+Fct,ij和法向接觸阻尼力Fd,ij=Fdn,ij+Fdt,i j

圖7 SDPH 粒子與DEM 顆粒之間的相互作用Fig.7 Interaction between SDPH particle and DEM particle

考慮DEM 粒子對SDPH 粒子作用的SDPH 方法動量方程為

考慮SDPH 粒子對DEM 粒子作用的DEM 方法的動量方程

式中,FDEM為DEM 粒子作用于SDPH 粒子上的作用力矢量,FSDPH為DEM 粒子作用于SDPH 粒子上的作用力矢量.

2.4 邊界條件

SDPH 和DEM 粒子與固壁邊界之間的作用力分為法向邊界力fn和切向邊界力fτ,如圖8 所示.

圖8 SDPH 和DEM 粒子邊界力施加方法Fig.8 Boundary force on SDPH and DEM particles

本文采用罰函數方法[67]計算SDPH 粒子與固壁邊界之間的法向接觸力,同時,為了保持算法的魯棒性和參數選擇的方便性,對罰參數進行了改進以保證法向力與SDPH 粒子到邊界的距離成反比

其中,ε 為罰參數,rij為粒子i和粒子j之間的徑矢,vi為粒子i的速度矢量,為邊界粒子j的速度矢量,nj為邊界粒子j的法向向量,如圖8 所示.Wij為粒子i和粒子j之間的核函數,Aj為邊界粒子的面積.DEM 顆粒和邊界之間的法向接觸力采用方程(38)所示的赫茲模型[65]施加.

顆粒與邊界之間的切向作用力fτ,采用摩擦模型計算.該模型認為顆粒與邊界的切向力與法向力成正比關系

其中,t anφ 為顆粒流與邊界的摩擦系數,φ 為摩擦角.

3 三維圓柱型顆粒堆坍塌過程數值模擬

顆粒堆坍塌問題是認識顆粒材料運動規律、檢驗模型方法有效性的基礎案例.很多學者已經對二維單側顆粒堆坍塌過程進行了數值模擬,但由于在厚度方向受限于前后壁面的制約,很多三維的現象無法獲取.同時,之前的坍塌過程數值模擬工況僅限于長徑比較小的案例,顆粒基本上遠離快速流態化狀態,也就是基本是處于類固態和類液態,采用濃密顆粒介質的本構理論便可實現有效計算,但對于長徑比較大的案例未曾涉及.這里選取不同長徑比下的軸對稱圓柱型顆粒堆坍塌算例進行數值模擬,捕捉全三維顆粒堆坍塌過程中一些典型現象,如截錐形結構、圓錐形結構、圓球形頂部結構以及墨西哥帽結構等,檢驗新的理論和方法在描述這種全三維、多相態問題上的有效性,同時深入認識和理解其物理機理.

算例模型示意圖如圖9 所示.三維圓柱型顆粒堆高度為H0,半徑為R0(具體數值根據工況不同而發生改變),長徑比用a來表示,a=H0/R0.在實驗過程中,初始由空心圓形容器固定,在t=0 時刻,將空心圓形容器沿垂直方向提升,并假定容器移除速度非常快,顆粒流不受其影響.而數值模擬實施過程中,顆粒堆從零時刻開始計算,即表征圓形容器已經移除.顆粒介質由初始靜止狀態,開始沿徑向塌陷,重力勢能轉化為運動動能,像流體一樣流動.在流動的過程中,存在于顆粒堆中的能量逐漸耗散,流動逐漸轉變為準靜態狀態,流速降低,流動性減小,形成一個對稱的堆積體.假定最終沉積時的狀態中顆粒介質的高度用Hf,半徑用Rf表示,通過相關實驗研究[14-16]表明,顆粒的材料、粒度、表面粗糙度以及初始顆粒堆的幾何形態均對最終的沉積形態以及顆粒流的運動過程具有顯著的影響,其中,初始長徑比起著至關重要的作用.因此,本文重點通過改變結構的初始參數(堆體的高度、半徑和體積)獲得顆粒流動的標度定律,與實驗進行對比驗證.同時,通過分析濃密顆粒流中顆粒的流動以及最終沉積規律,揭示形成不同流動特征的機理,為實際自然現象規律的揭示提供支撐.

圖9 模型示意圖Fig.9 Geometric model of granular pile

算例中的顆粒材料設定為干燥無黏性沙子,顆粒直徑為0.32 mm,實際密度為2600 kg/m3,初始體積分數為0.6,體積密度為1560 kg/m3,彈性模量為20 GPa,泊松比0.3,內摩擦角30°.采用SDPH 方法進行初始離散,SDPH 粒子的密度為顆粒的有效密度1560 kg/m3,初始體積分數為0.6,SDPH 粒子的直徑為5 mm,粒子總數量根據工況不同而發生改變,光滑長度為6.5 mm.底部邊界與顆粒之間有法向作用力和切向摩擦力存在,法向接觸力按照閾值函數方法計算,切向摩擦力按照摩擦模型計算,摩擦力系數 t anφ=0.6 .該算例中,決定相態轉變的體積分數臨界值來源于不同模型的適用范圍數值,屬于經驗參數,這里取 φl,min=0.55,φg,max=0.5,φg,min=0.02 .

首先,為了驗證本文所選取的SDPH 粒子直徑的合理性,在長徑比a=0.72 條件下對改變不同的SDPH 粒子初始直徑進行了計算,并與文獻實驗結果[6]進行了對比,如圖10 所示.從結果來看,隨著初始粒徑的減小,計算的精度也越高,但是計算量隨之增加,初始粒徑為5 mm 與2.5 mm 的結果相差較少,同時與實驗吻合較好.因此,為了保持足夠的精度,同時計算量又得到控制,本文選取的直徑5 mm 是較為合理的.

圖10 不同SDPH 粒子直徑下計算結果與實驗[6]對比Fig.10 Comparison between calculated results and experiments[6]under different SDPH particle diameters

圖11 為計算獲得的在a=0.55 工況下不同時刻的顆粒運動過程及最終沉積形態圖.可以看到,在150 ms 時刻,處于圓堆體上部外層的顆粒已經開始斜向下、向四周坍塌運動,而處于內層的顆粒則一直處于靜止狀態,它們之間存在一個明顯的間隔面將外部坍塌區域與內部非變性區域分開,這在固體力學領域稱為剪切帶.隨著時間的發展,剪切帶進一步向內部中心移動,外層坍塌的區域進一步擴展.由于底部摩擦力的作用以及顆粒-顆粒間接觸作用和重力作用的共同影響,在最終顆粒停止運動后,在堆體頂部留有一部分未受任何干擾的區域,形成類似“平頭帽”狀形態,該區域保持在初始高度H0,并且其與周圍自由表面坡形成一個約為靜態休止角的狀態.

圖12 為計算獲得的在a=0.9 工況下,不同時刻的顆粒運動過程及最終沉積形態圖.相比a=0.55 工況下的計算結果,該流動更加復雜.它不像a=0.55工況下直接在靜態區和非靜態區之間形成分割面,而是由外層一層層向內坍塌擴散,分割面不清晰,自由表面的坡度從流動前沿到堆體的上頂部逐漸變陡.隨著屈服顆粒一層層的發展侵蝕,最終整個顆粒堆的上表面完全達到屈服狀態,僅在最高峰處留下一個尖尖的錐形角,最終形成的坡度的角度也明顯小于靜態休止角.與圖11 不同的是,由于該工況下最終堆積形態的最高峰不再是初始高度值,因此數值模擬結果采用高度等高線進行云圖顯示.通過圖像可以看出本文的數值模擬很好地捕獲到了該實驗現象,與實驗圖像吻合非常好.

圖11 a=0.55 工況下顆粒堆坍塌過程與實驗[6]對比Fig.11 Comparison of the calculated collapse of the granular column with that of the experiment[6] under the condition of a=0.55

圖12 a=0.9 工況下顆粒堆坍塌過程與實驗[6]對比Fig.12 Comparison of the calculated collapse of the granular column with that of the experiment[6] under the condition of a=0.9

圖13 為a=3.0 工況下顆粒堆在不同時刻的運動過程形態計算結果與實驗對比.可以看出從0 時刻開始整個上表面開始向下運動,位于堆體底部外側的顆粒則開始向外運動,在顆粒堆前沿兩側觀察到較大的速度.由于堆體高度較大,堆體頂部的粒子以較高的速度坍塌,可以看到處于堆體周圍的顆粒已經有少部分處于快速相態化狀態,但位于堆體頂部的顆粒表面形態不發生改變.在堆體向下運動一段距離后,由于堆體周圍顆粒受重力和剪切力作用,堆體頂部變形形成一個凸形圓頂形態,其曲率半徑也隨著時間的發展而逐漸變小.由于水平面上速度較大的粒子可以在水平表面的兩個方向上同時對稱地運動,因此,其逐漸轉變成正弦函數的形狀.通過與Lajeunesse 等[9]通過剖開顆粒堆積體形態可以看到,最終的形態可以分為三個區域,一個是中心未發生任何改變的靜態區域,半徑與初始堆體半徑相同;一個是顆粒流先到達的區域;另一個是外部的擴展流動區域.本文數值模擬同樣觀察到了該現象,與實驗吻合較好.

圖13 a=3.0 工況下顆粒堆坍塌過程與實驗[6]對比Fig.13 Comparison of the calculated collapse of the granular column with that of the experiment[6] under the condition of a=3.0

繼續增大初始長徑比至a=13.8,顆粒堆的表面速度擴展更大,從圖14(b)和14(c)顯示的堆體部分顆粒的分布即可看出,堆體表面的顆粒已經出現很大的速度波動,處于明顯的快速相態化狀態.在484 ms 上部堆體已經完全消失,由于向下運動的顆粒堆速度較高,到達鋪展表面的顆粒仍有垂直向下運動的趨勢,造成對于鋪展表面顆粒堆的沖擊,中心鋪展區域形成一個環形凹狀結構,隨著時間推移,凹狀結構進一步擴大,同時能量進一步耗散,直到514 s 后基本達到穩定狀態.在此基礎上進一步增大長徑比至a=16.7,如圖15 所示,凹狀結構更加明顯,形成明顯的墨西哥帽形態,中心是一個尖型的錐角,四周形成一個脊狀.通過與實驗對比,除了中心錐角的角度稍大于實驗獲得的錐角角度之外,脊狀的位置、脊狀的高度、鋪展的范圍等數值均與實驗[6,10]吻合較好.

圖14 a=13.8 工況下顆粒堆坍塌過程與實驗[6]對比Fig.14 Comparison of the calculated collapse of the granular column with that of the experiment[6] under the condition of a=13.8

圖15 a=16.7 工況下顆粒堆最終鋪展形態與實驗[10]對比Fig.15 Comparison of the calculated final spreading morphology of the granular column with that of experimental results[10] under the condition of a=16.7

在本文模擬算例中,對于較大長徑比的三維顆粒堆來說,其顆粒運動的速度較大,鋪展的范圍也很大,存在一些顆粒從主體顆粒堆中分離出來,不再遵循顆粒的類固態和類液態假設.顆粒的體積分數變化較大,顆粒之間的接觸不再以長時碰撞接觸為主,屬于快速顆粒流動,具有較大的慣性數,剪切速率較高,堆積的圍壓接近于零,這時必須采用類氣態顆粒流模型和離散單元模型進行求解.

圖16 展示了a值為16.7 條件下顆粒堆在坍塌中所經歷的不同相態過程,黑色表示顆粒處于濃密態(顆粒固相和液相),紅色表示顆粒處于稀疏態(顆粒氣相和液氣過渡狀態),白色表示顆粒處于超稀疏態(顆粒離散慣性相).可以看出,由于初始長徑比較大,顆粒到達壁面鋪展時速度較高基本上均處于快速顆粒流狀態,而處于邊緣附近的部分少數顆粒由于速度更大,體積分數較小,基本上脫離顆粒堆的運動,達到超稀疏狀態,計算很好地捕捉到了相態演變過程.

圖16 a=16.7 三維圓柱型顆粒堆坍塌過程不同相態演變(黑色表示顆粒處于濃密態,紅色表示顆粒處于稀疏態,白色表示顆粒處于超稀疏態)Fig.16 The evolution of different phases during the collapse of cylindrical granular pile (Black indicates that the particles are in dense state,red indicates that the particles are in dilute state,and white indicates that the particles are in ultra-dilute state)

圖17 為a<1.7 情況下顆粒鋪展范圍隨高徑比的不同變化曲線,該曲線理論上是一條線性直線,因為在這個長徑比范圍內,顆粒鋪展的范圍有限,中心部分的顆粒堆基本保持不變,只有超出一定半徑范圍之后的顆粒參與鋪展運動,初始顆粒堆的高度直接決定了該鋪展范圍,而與初始的顆粒堆半徑無關,因此該直線等價于鋪展范圍r∞?r0與初始顆粒堆高度h0之間的關系.從量綱角度分析其必定是線性關系,實驗擬合數據比例系數為1.24,從圖17 可以看出數值模擬結果與實驗結果[6]對比較好.

圖17 顆粒堆鋪展范圍隨著a 不同的變化曲線(a <1.7)Fig.17 The runout range of the granular column as a function of a (a <1.7)

圖18 為不同a值條件下顆粒堆高度與鋪展范圍半徑之間的關系曲線,可以看出在a值較小的情況下(a<1.0),由于顆粒堆頂部部分中心顆粒未受到坍塌擾動作用,其高度與鋪展面積之間無相關關系.在1 .0

圖18 不同a 值條件下顆粒堆高度與鋪展范圍半徑之間的關系曲線Fig.18 The height of the granular column as a function of spreading radius under different conditions of a

圖19 為三種不同a值條件下顆粒堆鋪展范圍隨時間變化曲線,隨著a值增加鋪展的時間增長,鋪展的范圍增加,曲線基本上呈S 形,在初始階段鋪展的較為緩慢,主要與速度啟動存在一定的滯后有關,最終在底部摩擦力和內部剪切作用力下能量耗散,緩慢回到靜止狀態.從數值模擬結果和實驗[14]對比可以看出,兩者數據吻合較好.

圖19 顆粒堆鋪展范圍隨時間變化曲線Fig.19 The runout range of the granular column over time

最后,為了展示本文所采用的數值方法在計算效率方面的優勢,對相同工況(a=0.55)和相同粒子數情況下采用本文方法和DEM 方法的計算耗時進行了對比,具體DEM 計算是采用商業軟件EDEM 2018 完成.對比數據如表1 所示.由于DEM 方法時間步長受顆粒的硬度和直徑影響較大,無法采用更大的時間步長.而SDPH 方法雖然也受顆粒屬性的影響,但是該顆粒屬性為宏觀局部的平均值,計算獲得的時間步長更大些.因此,計算獲得同樣的結果,DEM 計算的總步數較多,同時每一時間步對于接觸力的計算耗時相較SDPH 方法來說也較大一些,總耗時較長.

表1 本文方法與DEM 方法計算耗時對比Table 1 Comparison of calculation time between this method and DEM method

4 結論

論文針對三維圓柱型顆粒堆坍塌過程中表現出來的堆積、緩慢流動、快速流動、分離流動等不同相態問題開展了數值模擬研究,成功再現了不同工況下的顆粒堆坍塌過程,獲得的主要結論如下.

(1)建立的全相態物理模型和數值方法對于顆粒介質運動過程中呈現出來的不同相態可以較好地捕捉,同時由于采用的數值方法中除了超稀疏區域外均為基于連續介質力學的粒子類方法,計算所使用的參數均是宏觀實驗可獲得的參量,克服了采用單一單顆粒軌道追蹤依賴于單顆粒屬性的缺陷;同時,采用的粒子方法屬于粗粒度方法,既可以實現對顆粒的拉格朗日運動軌跡再現,又可以降低實際單顆粒建模所帶來的計算量大的問題.

(2)三維圓柱型顆粒堆在不同初始高寬比條件下最終沉積物形態存在較大的差別,可形成平凸形、尖角形、圓頂形、墨西哥帽形等,在高寬比小于1.7 范圍內顆粒堆鋪展范圍與高寬比近似成直線形關系.

(3)雖然本文采用基于連續介質力學模型的數值方法對顆粒介質運動問題模擬計算取得了較好的結果,但不可否認此類方法也存在一定的局限:一方面此類方法必須依賴于連續介質力學模型的發展,而顆粒介質的宏觀理論預測模型目前尚屬于發展階段,還需要不斷改進和完善;另一方面此類方法描述的是物質的連續介質力學行為,每個單元參量表征的是物質局部的統計平均特性,即使本文采用粒子方法獲得了顆粒的運動軌跡,但是這種粒子是屬于一種“偽顆粒”,而非真實的顆粒,無法獲得真實單顆粒的運動特征以及真實單顆粒間的相互作用行為;第三方面由于此類方法計算顆粒運動問題經歷了連續介質力學模型的建立和數值模擬方法的使用兩個關鍵過程,每個過程中均存在一定的假設和簡化,理論上計算的精度相比單顆粒直接模擬的微觀方法來說存在一定的差距,需要不斷發展和提高物理模型和數值方法的準確度.