陳福振 李亞雄 史騰達(dá) 嚴(yán)紅

(西北工業(yè)大學(xué)動(dòng)力與能源學(xué)院,西安 710072)

(西北工業(yè)大學(xué)太倉(cāng)長(zhǎng)三角研究院,江蘇太倉(cāng) 215400)

引言

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

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

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

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

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

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

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

1 顆粒介質(zhì)全相態(tài)物理模型

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

以上是目前國(guó)內(nèi)外對(duì)于顆粒介質(zhì)存在的三種相態(tài)的描述,作者通過(guò)分析發(fā)現(xiàn)當(dāng)顆粒的體積分?jǐn)?shù)繼續(xù)減小,顆粒間的二體碰撞假設(shè)也不再滿足,顆粒間碰撞的概率已經(jīng)非常微小,顆粒主要以受到的體力作用和外部流場(chǎng)作用為主,這時(shí)從相對(duì)尺度上來(lái)說(shuō)顆粒屬于介觀尺度范圍,不再遵循擬流體的連續(xù)介質(zhì)力學(xué)定律,因此不適用于以上三種狀態(tài)描述,將這種狀態(tài)稱為離散慣性相,遵循質(zhì)點(diǎn)運(yùn)動(dòng)定律.由此,顆粒介質(zhì)從濃密到稀疏、從連續(xù)到離散、體積分?jǐn)?shù)從1 至0 的全部相態(tài)可全覆蓋,如圖1 所示.

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

本文在建立顆粒介質(zhì)全相態(tài)理論之前,根據(jù)顆粒稀疏程度的不同將全相態(tài)劃分為三個(gè)區(qū)域,將顆粒固相和液相區(qū)域統(tǒng)稱為濃密顆粒介質(zhì)區(qū)(φl(shuí),min≤φp≤φs,max,φp為顆粒體積分?jǐn)?shù),φs,max為顆粒固相最大體積分?jǐn)?shù),φl(shuí),min為顆粒液相最小體積分?jǐn)?shù),通常 φl(shuí),min≈φs,max);將顆粒液相和氣相之間的過(guò)渡相和氣相區(qū)域稱之為稀疏顆粒介質(zhì)區(qū)(φg,min≤φp≤φl(shuí),min);將離散慣性相區(qū)域稱之為超稀疏顆粒介質(zhì)區(qū)(φp<φg,min).劃分區(qū)域的目的是分區(qū)域建模,使復(fù)雜的多相態(tài)建模能夠?qū)崿F(xiàn)一定程度的簡(jiǎn)化.建立顆粒介質(zhì)的全部相態(tài)本構(gòu)理論就需要根據(jù)其所處的不同相態(tài)區(qū)域進(jìn)行建模,同時(shí)建立不同相態(tài)之間的轉(zhuǎn)變?cè)瓌t:

(1) 對(duì)于濃密顆粒介質(zhì)區(qū)域狀態(tài)采用彈-黏-塑性本構(gòu)理論描述,具體針對(duì)濃密顆粒介質(zhì)固相采用彈塑性理論描述,對(duì)于濃密顆粒介質(zhì)液相采用流變學(xué)理論描述;

(2)對(duì)于稀疏顆粒介質(zhì)區(qū)域狀態(tài)采用顆粒動(dòng)理學(xué)理論和摩擦動(dòng)力學(xué)理論描述,具體針對(duì)稀疏顆粒介質(zhì)的液氣過(guò)渡相采用顆粒動(dòng)理學(xué)和摩擦動(dòng)力學(xué)相結(jié)合的方式描述,對(duì)于稀疏顆粒介質(zhì)氣相采用顆粒動(dòng)理學(xué)描述;

(3)對(duì)于超稀疏顆粒介質(zhì)區(qū)域狀態(tài)采用質(zhì)點(diǎn)動(dòng)力學(xué)理論描述,具體就是指針對(duì)超稀疏顆粒介質(zhì)的慣性相采用質(zhì)點(diǎn)動(dòng)力學(xué)理論進(jìn)行描述.

下面就分別介紹這些不同的理論,以及這些理論之間的轉(zhuǎn)變?cè)瓌t.

1.1 描述濃密顆粒介質(zhì)運(yùn)動(dòng)狀態(tài)的彈-黏-塑性本構(gòu)模型(φl(shuí),min ≤φp ≤φs,max )

對(duì)于濃密顆粒介質(zhì)來(lái)說(shuō),首先介紹其運(yùn)動(dòng)控制方程,包括質(zhì)量守恒方程和動(dòng)量守恒方程兩個(gè),如下

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

式中,δαβ是克羅內(nèi)克函數(shù),當(dāng)α=β時(shí)δαβ=1,當(dāng)α≠β時(shí)δαβ=0.靜水壓力p直接采用本構(gòu)方程中的應(yīng)力計(jì)算得到

全應(yīng)力張量的具體形式將根據(jù)其所處的狀態(tài)不同而發(fā)生改變,具體見(jiàn)以下章節(jié)介紹.fα為其他外力,如本文中所需要考慮的重力.d/dt為全導(dǎo)數(shù).

有效耦合現(xiàn)有的彈塑性理論和基于流變學(xué)的黏塑性理論,獲得基于稠密顆粒運(yùn)動(dòng)彈-黏-塑性理論的顆粒固液相計(jì)算模型[52].假設(shè)準(zhǔn)靜態(tài)和流動(dòng)狀態(tài)之間的屈服準(zhǔn)則用Drucker-Prager 屈服準(zhǔn)則表示,則可以建立從完全彈性到彈性-微黏塑性再到完全彈性-黏塑性的過(guò)渡過(guò)程,如圖2 所示.

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

1.1.1 完全彈性下的顆粒介質(zhì)應(yīng)力-應(yīng)變關(guān)系

1.1.2 彈性-微小黏塑性狀態(tài)下的顆粒介質(zhì)應(yīng)力-應(yīng)變關(guān)系

隨著顆粒受力的增加,變形逐漸增大.彈性引起的剪應(yīng)力逐漸增大,首先達(dá)到.此時(shí),根據(jù)顆粒流變學(xué)理論,顆粒介質(zhì)開(kāi)始流動(dòng),和 μ(I)逐漸增加.雖然發(fā)生了流動(dòng),但流動(dòng)速度非常小,尚未達(dá)到塑性屈服狀態(tài).在此階段,只有剪應(yīng)力逐漸增大.

彈性-微小黏塑性狀態(tài)下顆粒介質(zhì)的法向應(yīng)力仍按線彈性本構(gòu)模型計(jì)算,如式(5),剪應(yīng)力按流變學(xué)模型計(jì)算,即

1.1.3 彈性-完全黏塑性狀態(tài)下的顆粒介質(zhì)應(yīng)力-應(yīng)變關(guān)系

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

1.2 描述稀疏顆粒介質(zhì)區(qū)域的顆粒動(dòng)理學(xué)理論和摩擦動(dòng)力學(xué)理論(φg,min ≤φp<φl(shuí),min )

1.2.1 描述顆粒氣相的動(dòng)理學(xué)理論(KTGF)

對(duì)于稀疏顆粒流區(qū)域中的顆粒氣相,采用顆粒動(dòng)理學(xué)理論(簡(jiǎn)稱KTGF)[54-55]進(jìn)行描述,一共包括質(zhì)量守恒、動(dòng)量守恒、擬溫度守恒等三個(gè)方程,質(zhì)量守恒方程和動(dòng)量守恒方程如式(1)和式(2),擬溫度守恒方程式為顆粒動(dòng)理學(xué)理論所特有,如下

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

式中,φs,max為顆粒介質(zhì)可達(dá)到的最大體積分?jǐn)?shù)值.

1.2.2 描述顆粒氣相-液相過(guò)渡區(qū)的摩擦動(dòng)力學(xué)理論

摩擦動(dòng)力學(xué)主要用來(lái)描述處于長(zhǎng)程接觸以及存在不可忽略的摩擦應(yīng)力的顆粒.一些學(xué)者研究表明,摩擦動(dòng)力學(xué)可以用于體積分?jǐn)?shù)較高的顆粒[56-57].在稠密顆粒流和稀疏顆粒流之間,必然存在一個(gè)從長(zhǎng)程接觸到瞬態(tài)碰撞接觸的過(guò)渡階段.這里,摩擦動(dòng)力學(xué)被用來(lái)彌補(bǔ)顆粒氣相和液相之間存在的過(guò)渡缺陷問(wèn)題.Johnson 等[58-59]提出的半經(jīng)驗(yàn)?zāi)Ｐ陀脕?lái)計(jì)算由摩擦引起的法向應(yīng)力,如下所示

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

式中,? 為內(nèi)摩擦角.通過(guò)分析可知[53,59],摩擦動(dòng)力學(xué)可以被認(rèn)為是 μ(I) 流變學(xué)[61-62]在慣性數(shù)或者說(shuō)剪切應(yīng)變率無(wú)限大時(shí)的極限情況.

1.3 描述超稀疏顆粒流區(qū)域的質(zhì)點(diǎn)動(dòng)力學(xué)理論(φp<φg,min )

離散質(zhì)點(diǎn)動(dòng)力學(xué)是指針對(duì)具有一定質(zhì)量但幾何尺寸大小可以忽略的物體,采用牛頓運(yùn)動(dòng)定律進(jìn)行描述的動(dòng)力學(xué)過(guò)程

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

1.4 濃密顆粒流與稀疏顆粒流兩個(gè)狀態(tài)之間的轉(zhuǎn)變?cè)瓌t

KTGF 結(jié)合摩擦動(dòng)力學(xué)模型[56-57]用于計(jì)算顆粒液相和氣相之間的轉(zhuǎn)換.一些學(xué)者通過(guò)研究發(fā)現(xiàn),摩擦動(dòng)力學(xué)可以用于具有較高體積分?jǐn)?shù)的顆粒,但當(dāng)體積分?jǐn)?shù)增加到某個(gè)值時(shí)存在不確定性.然而,將摩擦動(dòng)力學(xué)與彈-黏-塑性模型相結(jié)合可以克服這一缺點(diǎn).在轉(zhuǎn)變過(guò)渡的界面處應(yīng)力應(yīng)保持守恒,彈性-黏塑性模型計(jì)算的法向應(yīng)力和剪應(yīng)力應(yīng)與摩擦動(dòng)力學(xué)計(jì)算所得的法向應(yīng)力和剪應(yīng)力相同.顆粒應(yīng)力計(jì)算公式表示如下

式中,pEVP和 μEVP分別為采用彈-黏-塑性本構(gòu)模型計(jì)算獲得的正應(yīng)力和剪應(yīng)力,pKTGF和pfriction分別是顆粒壓力pp中的動(dòng)理學(xué)部分和摩擦部分,μKTGF和μfriction分別是顆粒剪切黏度 μp的動(dòng)理學(xué)部分和摩擦部分.φg,max是摩擦應(yīng)力開(kāi)始逐漸增加時(shí)的顆粒體積分?jǐn)?shù).小于 φg,max值時(shí),通過(guò)已知資料發(fā)現(xiàn)未觀察到顆粒間的摩擦行為,因此假設(shè),當(dāng)均勻分布的粒子不再接觸時(shí),摩擦相互作用不再發(fā)生,主要由碰撞相互作用產(chǎn)生相互間的應(yīng)力.

顆粒液相和氣相之間的過(guò)渡如圖3 所示.從顆粒液相和顆粒氣相之間的過(guò)渡相開(kāi)始,顆粒擬溫度從零開(kāi)始升高,表明顆粒之間的碰撞逐漸加劇.從過(guò)渡相到氣相轉(zhuǎn)變的界面上,顆粒之間的摩擦完全消失,轉(zhuǎn)變?yōu)橥耆w碰撞.從上述計(jì)算公式可以看出,用這種方法建立的過(guò)渡態(tài)可以有效地連接液相和氣相,實(shí)現(xiàn)平穩(wěn)過(guò)渡,符合物理定律.

圖3 濃密顆粒流與稀疏顆粒流之間的過(guò)渡轉(zhuǎn)化Fig.3 Transition between dense granular flow and dilute granular flow

1.5 稀疏顆粒流與超稀疏顆粒流兩個(gè)狀態(tài)之間的轉(zhuǎn)變?cè)瓌t

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

2 顆粒介質(zhì)全相態(tài)數(shù)值模擬方法

第2 節(jié)闡述了顆粒介質(zhì)全相態(tài)理論,要對(duì)顆粒介質(zhì)全相態(tài)進(jìn)行數(shù)值模擬除了理論模型之外,還需要引入合適的數(shù)值模擬方法對(duì)模型進(jìn)行離散求解.本文對(duì)于顆粒介質(zhì)全相態(tài)中的類固態(tài)、類液態(tài)和類氣態(tài)三種狀態(tài)采用SDPH 方法進(jìn)行模擬,因?yàn)檫@三種狀態(tài)的本構(gòu)模型均是基于宏觀連續(xù)介質(zhì)力學(xué)所建立的,與SDPH 的理論基礎(chǔ)一致.當(dāng)散體顆粒分散到一定程度即顆粒相的體積分?jǐn)?shù)小到一定程度時(shí),采用2.3 節(jié)的離散顆粒動(dòng)力學(xué)進(jìn)行描述,引入描述單一顆粒行為的離散單元法進(jìn)行計(jì)算.

2.1 基于顆粒連續(xù)介質(zhì)力學(xué)的光滑離散顆粒流體動(dòng)力學(xué)方法

為明確SPH 粒子與顆粒之間的對(duì)應(yīng)關(guān)系,Chen 等[47-51]前期已經(jīng)從顆粒動(dòng)理學(xué)角度出發(fā),將顆粒相視為擬流體,擬流體區(qū)域采用SPH 方法離散求解,同時(shí)將傳統(tǒng)SPH 方法改造成了適用于離散顆粒相求解的SDPH 方法,本文在此基礎(chǔ)上繼續(xù)拓展SDPH 方法的應(yīng)用范圍,對(duì)于顆粒固相和液相同樣采用SDPH 離散闡述,與顆粒氣相的SDPH 描述保持一致.

2.1.1 光滑離散顆粒流體動(dòng)力學(xué)方法

在SDPH 方法中,實(shí)際顆粒的質(zhì)量、速度、壓力、位置、粒徑分布、體積分?jǐn)?shù)等均在計(jì)算粒子上有體現(xiàn).同時(shí),顆粒氣相的擬溫度值也作為計(jì)算粒子的一個(gè)參量.計(jì)算粒子與實(shí)際顆粒之間的參量對(duì)應(yīng)關(guān)系如下

其中,p為顆粒下標(biāo),φp為顆粒體積分?jǐn)?shù),ρp為顆粒實(shí)際密度.顆粒的有效密度可進(jìn)一步推導(dǎo)獲得

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

從以上對(duì)應(yīng)關(guān)系可以看出,當(dāng)使用SDPH 方法進(jìn)行計(jì)算時(shí),每個(gè)計(jì)算粒子可以代表一定數(shù)量的真實(shí)顆粒群,計(jì)算粒子的數(shù)量可以根據(jù)計(jì)算精度的要求確定.對(duì)于某些精度要求不太嚴(yán)格的問(wèn)題,計(jì)算量將大大減少,計(jì)算效率將顯著提高.

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

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

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

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

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

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

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

2.1.3 顆粒氣相的SDPH 離散公式

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

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

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

2.1.4 基于SDPH 方法的濃密顆粒介質(zhì)與稀疏顆粒介質(zhì)相態(tài)轉(zhuǎn)變

采用SDPH 方法不論是求解濃密顆粒相還是稀疏顆粒相,計(jì)算粒子與實(shí)際顆粒之間的對(duì)應(yīng)關(guān)系不會(huì)發(fā)生改變,只不過(guò)在求解的本構(gòu)模型上存在不同.決定顆粒介質(zhì)是處于濃密態(tài)還是稀疏態(tài)的參數(shù)是體積分?jǐn)?shù),體積分?jǐn)?shù)的計(jì)算公式為.圖4 展示了轉(zhuǎn)變策略.如果計(jì)算粒子M的體積分?jǐn)?shù) φp大于φl(shuí),min,它表示該粒子正處于顆粒固相和液相.當(dāng) φp降低至小于 φl(shuí),min,它表示顆粒已經(jīng)轉(zhuǎn)變?yōu)橄∈桀w粒狀態(tài)或顆粒氣相.將變換后的計(jì)算粒子標(biāo)記為N.粒子M和N之間的變量關(guān)系如圖4 所示.

圖4 基于SDPH 方法的濃密顆粒介質(zhì)與稀疏顆粒介質(zhì)相態(tài)轉(zhuǎn)變Fig.4 Phase transition between dense granular media and dilute granular media based on SDPH method

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

同樣地,對(duì)于顆粒從稀疏態(tài)轉(zhuǎn)化為濃密態(tài)時(shí),粒子變量值應(yīng)該保持不變,包括粒子的密度、速度、能量、坐標(biāo)等.兩個(gè)狀態(tài)的粒子存在的差別主要在內(nèi)力的計(jì)算上.從稀疏態(tài)到濃密態(tài)過(guò)程中,顆粒的擬溫度值逐漸降低,顆粒間的碰撞頻次也逐漸降低,而由摩擦作用產(chǎn)生的正應(yīng)力和剪切力數(shù)值則逐漸增加.在轉(zhuǎn)化的界面上,由摩擦作用產(chǎn)生的正應(yīng)力和剪切力與顆粒液相計(jì)算的正應(yīng)力和剪切力相同.在轉(zhuǎn)化為濃密顆粒態(tài)之后,顆粒間的塑性剪切力逐漸降低,彈性剪切力逐漸增加,顆粒速度逐漸降低,體積分?jǐn)?shù)逐漸增加,按照塑性流動(dòng)法則計(jì)算顆粒介質(zhì)材料的卸載過(guò)程,直至恢復(fù)到靜止?fàn)顟B(tài).

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

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

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

2.2 基于顆粒質(zhì)點(diǎn)顆粒動(dòng)力學(xué)的離散單元法

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

2.3 SDPH 與DEM 之間算法的耦合

2.3.1 SDPH 與DEM 之間算法的轉(zhuǎn)化

在稀疏顆粒流SDPH 粒子的體積分?jǐn)?shù)降到一定閾值后(φg,min)時(shí),其不再遵循二體碰撞假設(shè)的顆粒動(dòng)理學(xué)模型,因此將SDPH 粒子轉(zhuǎn)化為DEM 粒子進(jìn)行計(jì)算.圖6 顯示了SDPH 和DEM 之間算法的轉(zhuǎn)化方案.

圖6 SDPH 和DEM 算法之間的轉(zhuǎn)化方案Fig.6 Conversion between SDPH and DEM

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

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

假如在計(jì)算顆粒介質(zhì)過(guò)程中,空間中既有SDPH 粒子,又有DEM 粒子,當(dāng)他們之間處于相互接觸的狀態(tài)時(shí),需要計(jì)算他們之間的接觸作用力.采用DEM 粒子之前的接觸力計(jì)算方法進(jìn)行計(jì)算.因此,需要將處于接觸狀態(tài)的SDPH 粒子轉(zhuǎn)化為DEM 虛粒子進(jìn)行計(jì)算(等效兩個(gè)DEM 粒子接觸計(jì)算).圖7 展示了SDPH 粒子與DEM 顆粒之間的接觸力計(jì)算規(guī)則.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 粒子對(duì)SDPH 粒子作用的SDPH 方法動(dòng)量方程為

考慮SDPH 粒子對(duì)DEM 粒子作用的DEM 方法的動(dòng)量方程

式中,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

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

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

顆粒與邊界之間的切向作用力fτ,采用摩擦模型計(jì)算.該模型認(rèn)為顆粒與邊界的切向力與法向力成正比關(guān)系

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

3 三維圓柱型顆粒堆坍塌過(guò)程數(shù)值模擬

顆粒堆坍塌問(wèn)題是認(rèn)識(shí)顆粒材料運(yùn)動(dòng)規(guī)律、檢驗(yàn)?zāi)Ｐ头椒ㄓ行缘幕A(chǔ)案例.很多學(xué)者已經(jīng)對(duì)二維單側(cè)顆粒堆坍塌過(guò)程進(jìn)行了數(shù)值模擬,但由于在厚度方向受限于前后壁面的制約,很多三維的現(xiàn)象無(wú)法獲取.同時(shí),之前的坍塌過(guò)程數(shù)值模擬工況僅限于長(zhǎng)徑比較小的案例,顆?；旧线h(yuǎn)離快速流態(tài)化狀態(tài),也就是基本是處于類固態(tài)和類液態(tài),采用濃密顆粒介質(zhì)的本構(gòu)理論便可實(shí)現(xiàn)有效計(jì)算,但對(duì)于長(zhǎng)徑比較大的案例未曾涉及.這里選取不同長(zhǎng)徑比下的軸對(duì)稱圓柱型顆粒堆坍塌算例進(jìn)行數(shù)值模擬,捕捉全三維顆粒堆坍塌過(guò)程中一些典型現(xiàn)象,如截錐形結(jié)構(gòu)、圓錐形結(jié)構(gòu)、圓球形頂部結(jié)構(gòu)以及墨西哥帽結(jié)構(gòu)等,檢驗(yàn)新的理論和方法在描述這種全三維、多相態(tài)問(wèn)題上的有效性,同時(shí)深入認(rèn)識(shí)和理解其物理機(jī)理.

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

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

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

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

圖10 不同SDPH 粒子直徑下計(jì)算結(jié)果與實(shí)驗(yàn)[6]對(duì)比Fig.10 Comparison between calculated results and experiments[6]under different SDPH particle diameters

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

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

圖11 a=0.55 工況下顆粒堆坍塌過(guò)程與實(shí)驗(yàn)[6]對(duì)比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 工況下顆粒堆坍塌過(guò)程與實(shí)驗(yàn)[6]對(duì)比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 工況下顆粒堆在不同時(shí)刻的運(yùn)動(dòng)過(guò)程形態(tài)計(jì)算結(jié)果與實(shí)驗(yàn)對(duì)比.可以看出從0 時(shí)刻開(kāi)始整個(gè)上表面開(kāi)始向下運(yùn)動(dòng),位于堆體底部外側(cè)的顆粒則開(kāi)始向外運(yùn)動(dòng),在顆粒堆前沿兩側(cè)觀察到較大的速度.由于堆體高度較大,堆體頂部的粒子以較高的速度坍塌,可以看到處于堆體周圍的顆粒已經(jīng)有少部分處于快速相態(tài)化狀態(tài),但位于堆體頂部的顆粒表面形態(tài)不發(fā)生改變.在堆體向下運(yùn)動(dòng)一段距離后,由于堆體周圍顆粒受重力和剪切力作用,堆體頂部變形形成一個(gè)凸形圓頂形態(tài),其曲率半徑也隨著時(shí)間的發(fā)展而逐漸變小.由于水平面上速度較大的粒子可以在水平表面的兩個(gè)方向上同時(shí)對(duì)稱地運(yùn)動(dòng),因此,其逐漸轉(zhuǎn)變成正弦函數(shù)的形狀.通過(guò)與Lajeunesse 等[9]通過(guò)剖開(kāi)顆粒堆積體形態(tài)可以看到,最終的形態(tài)可以分為三個(gè)區(qū)域,一個(gè)是中心未發(fā)生任何改變的靜態(tài)區(qū)域,半徑與初始堆體半徑相同;一個(gè)是顆粒流先到達(dá)的區(qū)域;另一個(gè)是外部的擴(kuò)展流動(dòng)區(qū)域.本文數(shù)值模擬同樣觀察到了該現(xiàn)象,與實(shí)驗(yàn)吻合較好.

圖13 a=3.0 工況下顆粒堆坍塌過(guò)程與實(shí)驗(yàn)[6]對(duì)比Fig.13 Comparison of the calculated collapse of the granular column with that of the experiment[6] under the condition of a=3.0

繼續(xù)增大初始長(zhǎng)徑比至a=13.8,顆粒堆的表面速度擴(kuò)展更大,從圖14(b)和14(c)顯示的堆體部分顆粒的分布即可看出,堆體表面的顆粒已經(jīng)出現(xiàn)很大的速度波動(dòng),處于明顯的快速相態(tài)化狀態(tài).在484 ms 上部堆體已經(jīng)完全消失,由于向下運(yùn)動(dòng)的顆粒堆速度較高,到達(dá)鋪展表面的顆粒仍有垂直向下運(yùn)動(dòng)的趨勢(shì),造成對(duì)于鋪展表面顆粒堆的沖擊,中心鋪展區(qū)域形成一個(gè)環(huán)形凹狀結(jié)構(gòu),隨著時(shí)間推移,凹狀結(jié)構(gòu)進(jìn)一步擴(kuò)大,同時(shí)能量進(jìn)一步耗散,直到514 s 后基本達(dá)到穩(wěn)定狀態(tài).在此基礎(chǔ)上進(jìn)一步增大長(zhǎng)徑比至a=16.7,如圖15 所示,凹狀結(jié)構(gòu)更加明顯,形成明顯的墨西哥帽形態(tài),中心是一個(gè)尖型的錐角,四周形成一個(gè)脊?fàn)?通過(guò)與實(shí)驗(yàn)對(duì)比,除了中心錐角的角度稍大于實(shí)驗(yàn)獲得的錐角角度之外,脊?fàn)畹奈恢?、脊?fàn)畹母叨取佌沟姆秶葦?shù)值均與實(shí)驗(yàn)[6,10]吻合較好.

圖14 a=13.8 工況下顆粒堆坍塌過(guò)程與實(shí)驗(yàn)[6]對(duì)比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 工況下顆粒堆最終鋪展形態(tài)與實(shí)驗(yàn)[10]對(duì)比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

在本文模擬算例中,對(duì)于較大長(zhǎng)徑比的三維顆粒堆來(lái)說(shuō),其顆粒運(yùn)動(dòng)的速度較大,鋪展的范圍也很大,存在一些顆粒從主體顆粒堆中分離出來(lái),不再遵循顆粒的類固態(tài)和類液態(tài)假設(shè).顆粒的體積分?jǐn)?shù)變化較大,顆粒之間的接觸不再以長(zhǎng)時(shí)碰撞接觸為主,屬于快速顆粒流動(dòng),具有較大的慣性數(shù),剪切速率較高,堆積的圍壓接近于零,這時(shí)必須采用類氣態(tài)顆粒流模型和離散單元模型進(jìn)行求解.

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

圖16 a=16.7 三維圓柱型顆粒堆坍塌過(guò)程不同相態(tài)演變(黑色表示顆粒處于濃密態(tài),紅色表示顆粒處于稀疏態(tài),白色表示顆粒處于超稀疏態(tài))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 情況下顆粒鋪展范圍隨高徑比的不同變化曲線,該曲線理論上是一條線性直線,因?yàn)樵谶@個(gè)長(zhǎng)徑比范圍內(nèi),顆粒鋪展的范圍有限,中心部分的顆粒堆基本保持不變,只有超出一定半徑范圍之后的顆粒參與鋪展運(yùn)動(dòng),初始顆粒堆的高度直接決定了該鋪展范圍,而與初始的顆粒堆半徑無(wú)關(guān),因此該直線等價(jià)于鋪展范圍r∞?r0與初始顆粒堆高度h0之間的關(guān)系.從量綱角度分析其必定是線性關(guān)系,實(shí)驗(yàn)擬合數(shù)據(jù)比例系數(shù)為1.24,從圖17 可以看出數(shù)值模擬結(jié)果與實(shí)驗(yàn)結(jié)果[6]對(duì)比較好.

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

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

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

圖19 為三種不同a值條件下顆粒堆鋪展范圍隨時(shí)間變化曲線,隨著a值增加鋪展的時(shí)間增長(zhǎng),鋪展的范圍增加,曲線基本上呈S 形,在初始階段鋪展的較為緩慢,主要與速度啟動(dòng)存在一定的滯后有關(guān),最終在底部摩擦力和內(nèi)部剪切作用力下能量耗散,緩慢回到靜止?fàn)顟B(tài).從數(shù)值模擬結(jié)果和實(shí)驗(yàn)[14]對(duì)比可以看出,兩者數(shù)據(jù)吻合較好.

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

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

表1 本文方法與DEM 方法計(jì)算耗時(shí)對(duì)比Table 1 Comparison of calculation time between this method and DEM method

4 結(jié)論

論文針對(duì)三維圓柱型顆粒堆坍塌過(guò)程中表現(xiàn)出來(lái)的堆積、緩慢流動(dòng)、快速流動(dòng)、分離流動(dòng)等不同相態(tài)問(wèn)題開(kāi)展了數(shù)值模擬研究,成功再現(xiàn)了不同工況下的顆粒堆坍塌過(guò)程,獲得的主要結(jié)論如下.

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

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

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