三維圓柱型顆粒堆坍塌問題的全相態(tài)數(shù)值模擬1)

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

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

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

引言

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

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

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

另一類是基于宏觀連續(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]等,此類方法的特點是基于網(wǎng)格離散求解數(shù)理模型,已成功用于顆粒堆坍塌、剪切造粒、化工流化床、顆粒管道輸送等問題中,但此類方法存在的最大問題是必須依賴網(wǎng)格求解帶來的缺陷,如有限元在計算類固態(tài)和小變形的類液態(tài)時是合理的,但對于顆粒介質(zhì)宏觀大變形問題或者類氣態(tài)問題往往會出現(xiàn)網(wǎng)格的扭曲和纏繞,計算終止;而采用有限體積法解決顆粒類氣態(tài)問題和類液態(tài)問題較為合適,但類固態(tài)本構(gòu)模型無法求解,同時該方法無法獲得顆粒的運動軌跡,易產(chǎn)生偽擴散,不易加入顆粒蒸發(fā)、燃燒等物理化學(xué)模型.

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

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

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

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

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

以上是目前國內(nèi)外對于顆粒介質(zhì)存在的三種相態(tài)的描述,作者通過分析發(fā)現(xiàn)當(dāng)顆粒的體積分?jǐn)?shù)繼續(xù)減小,顆粒間的二體碰撞假設(shè)也不再滿足,顆粒間碰撞的概率已經(jīng)非常微小,顆粒主要以受到的體力作用和外部流場作用為主,這時從相對尺度上來說顆粒屬于介觀尺度范圍,不再遵循擬流體的連續(xù)介質(zhì)力學(xué)定律,因此不適用于以上三種狀態(tài)描述,將這種狀態(tài)稱為離散慣性相,遵循質(zhì)點運動定律.由此,顆粒介質(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)劃分為三個區(qū)域,將顆粒固相和液相區(qū)域統(tǒng)稱為濃密顆粒介質(zhì)區(qū)(φl,min≤φp≤φs,max,φp為顆粒體積分?jǐn)?shù),φs,max為顆粒固相最大體積分?jǐn)?shù),φl,min為顆粒液相最小體積分?jǐn)?shù),通常 φl,min≈φs,max);將顆粒液相和氣相之間的過渡相和氣相區(qū)域稱之為稀疏顆粒介質(zhì)區(qū)(φg,min≤φp≤φl,min);將離散慣性相區(qū)域稱之為超稀疏顆粒介質(zhì)區(qū)(φp<φg,min).劃分區(qū)域的目的是分區(qū)域建模,使復(fù)雜的多相態(tài)建模能夠?qū)崿F(xiàn)一定程度的簡化.建立顆粒介質(zhì)的全部相態(tài)本構(gòu)理論就需要根據(jù)其所處的不同相態(tài)區(qū)域進行建模,同時建立不同相態(tài)之間的轉(zhuǎn)變原則:

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

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

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

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

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

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

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

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

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

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

圖2 完全彈性到彈性-微小黏塑性再到完全彈性-黏塑性的過渡過程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)力逐漸增大,首先達到.此時,根據(jù)顆粒流變學(xué)理論,顆粒介質(zhì)開始流動,和 μ(I)逐漸增加.雖然發(fā)生了流動,但流動速度非常小,尚未達到塑性屈服狀態(tài).在此階段,只有剪應(yīng)力逐漸增大.

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

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

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

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

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

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

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

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

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

摩擦動力學(xué)主要用來描述處于長程接觸以及存在不可忽略的摩擦應(yīng)力的顆粒.一些學(xué)者研究表明,摩擦動力學(xué)可以用于體積分?jǐn)?shù)較高的顆粒[56-57].在稠密顆粒流和稀疏顆粒流之間,必然存在一個從長程接觸到瞬態(tài)碰撞接觸的過渡階段.這里,摩擦動力學(xué)被用來彌補顆粒氣相和液相之間存在的過渡缺陷問題.Johnson 等[58-59]提出的半經(jīng)驗?zāi)Ｐ陀脕碛嬎阌赡Σ烈鸬姆ㄏ驊?yīng)力,如下所示

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

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

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

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

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

1.4 濃密顆粒流與稀疏顆粒流兩個狀態(tài)之間的轉(zhuǎn)變原則

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

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

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

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

1.5 稀疏顆粒流與超稀疏顆粒流兩個狀態(tài)之間的轉(zhuǎn)變原則

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

2.1.3 顆粒氣相的SDPH 離散公式

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

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

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

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

采用SDPH 方法不論是求解濃密顆粒相還是稀疏顆粒相,計算粒子與實際顆粒之間的對應(yīng)關(guān)系不會發(fā)生改變,只不過在求解的本構(gòu)模型上存在不同.決定顆粒介質(zhì)是處于濃密態(tài)還是稀疏態(tài)的參數(shù)是體積分?jǐn)?shù),體積分?jǐn)?shù)的計算公式為.圖4 展示了轉(zhuǎn)變策略.如果計算粒子M的體積分?jǐn)?shù) φp大于φl,min,它表示該粒子正處于顆粒固相和液相.當(dāng) φp降低至小于 φl,min,它表示顆粒已經(jīng)轉(zhuǎn)變?yōu)橄∈桀w粒狀態(tài)或顆粒氣相.將變換后的計算粒子標(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)力采用彈性定律計算,剪應(yīng)力采用彈性剪應(yīng)力與塑性剪應(yīng)力加和的方式計算,在向稀疏顆粒介質(zhì)算法轉(zhuǎn)變時,這兩個作用力也應(yīng)該保持不變,即由濃密顆粒介質(zhì)計算的彈性正應(yīng)力轉(zhuǎn)化為稀疏顆粒介質(zhì)的摩擦正應(yīng)力,見式(17),由式(17)反向計算求得Fr作為不變量,而后根據(jù)體積分?jǐn)?shù) φp的變化更新由摩擦產(chǎn)生正應(yīng)力的值;對于摩擦剪應(yīng)力則繼續(xù)采用 μ (I) 流變學(xué)剪應(yīng)力公式計算 τ=μ(I)p|γ|/γ,只不過這時的p開始由式(17)計算;對于由長時接觸產(chǎn)生的彈性剪應(yīng)力則置零;以上是對于SDPH 采用摩擦動力學(xué)在過渡區(qū)產(chǎn)生的正應(yīng)力與剪應(yīng)力的數(shù)值計算,保證了轉(zhuǎn)化的動量的守恒;同時從轉(zhuǎn)化開始,擬溫度的值由零開始計算,從而由碰撞產(chǎn)生的正應(yīng)力和剪應(yīng)力逐漸增大,直到過渡區(qū)全部轉(zhuǎn)化為顆粒動理學(xué)模型計算.

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

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

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

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

2.2 基于顆粒質(zhì)點顆粒動力學(xué)的離散單元法

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

2.3 SDPH 與DEM 之間算法的耦合

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

圖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

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

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

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

圖16 a=16.7 三維圓柱型顆粒堆坍塌過程不同相態(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 情況下顆粒鋪展范圍隨高徑比的不同變化曲線,該曲線理論上是一條線性直線,因為在這個長徑比范圍內(nèi),顆粒鋪展的范圍有限,中心部分的顆粒堆基本保持不變,只有超出一定半徑范圍之后的顆粒參與鋪展運動,初始顆粒堆的高度直接決定了該鋪展范圍,而與初始的顆粒堆半徑無關(guān),因此該直線等價于鋪展范圍r∞?r0與初始顆粒堆高度h0之間的關(guān)系.從量綱角度分析其必定是線性關(guān)系,實驗擬合數(shù)據(jù)比例系數(shù)為1.24,從圖17 可以看出數(shù)值模擬結(jié)果與實驗結(jié)果[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值條件下顆粒堆高度與鋪展范圍半徑之間的關(guān)系曲線,可以看出在a值較小的情況下(a<1.0),由于顆粒堆頂部部分中心顆粒未受到坍塌擾動作用,其高度與鋪展面積之間無相關(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值條件下顆粒堆鋪展范圍隨時間變化曲線,隨著a值增加鋪展的時間增長,鋪展的范圍增加,曲線基本上呈S 形,在初始階段鋪展的較為緩慢,主要與速度啟動存在一定的滯后有關(guān),最終在底部摩擦力和內(nèi)部剪切作用力下能量耗散,緩慢回到靜止?fàn)顟B(tài).從數(shù)值模擬結(jié)果和實驗[14]對比可以看出,兩者數(shù)據(jù)吻合較好.

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

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

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

4 結(jié)論

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

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

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

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