鄒宇雄 馬 剛,2) 李易奧 王 頔 邱煥峰 周 偉
* (武漢大學水資源與水電工程科學國家重點實驗室,武漢 430072)
? (水工巖石力學教育部重點實驗室,武漢 430072)
** (中國電建集團貴陽勘測設(shè)計研究院有限公司,貴陽 550081)
顆粒物質(zhì)是由大量離散固體相互作用而形成的復(fù)雜體系,在自然界和工業(yè)界廣泛存在.以水利和巖土工程中的堆石體為例,由于取材方便,且對地形、地質(zhì)條件有較強的適應(yīng)性,在堆石壩、路堤、機場高填方地基等工程建設(shè)中得到廣泛應(yīng)用.這些工程一旦出現(xiàn)問題和隱患,將會嚴重威脅人民生命財產(chǎn)安全,因此迫切需要深入研究顆粒材料的物理力學特性.
由于顆粒材料具有結(jié)構(gòu)異質(zhì)性、各向異性和多尺度結(jié)構(gòu)等特點,其表現(xiàn)出復(fù)雜的物理力學特性,如剪脹[1]、應(yīng)變局部化[2],固?液轉(zhuǎn)變[3-4]等.對巖土類顆粒材料,自1963 年Roscoe 等[5]提出劍橋模型以來,許多學者根據(jù)顆粒材料的宏觀應(yīng)力變形特性發(fā)展了各種各樣的本構(gòu)理論[6-8].但目前還沒有統(tǒng)一的理論能夠解釋顆粒材料所有的應(yīng)力變形特性.蔣明鏡[9]指出這些基于唯象的土力學本構(gòu)理論雖然在巖土工程實踐中發(fā)揮了重要作用,但是難以從機理上刻畫巖土顆粒系統(tǒng)的非連續(xù)性、大變形、破壞等問題.因此,越來越多學者嘗試從細觀尺度研究顆粒材料復(fù)雜力學行為的起源,并試圖將顆粒形狀、表面摩擦、粒徑分布等顆粒尺度的特性與宏觀力學響應(yīng)建立聯(lián)系.
巖土顆粒材料的形狀各異,大量研究表明顆粒形狀對顆粒材料的剪切強度[10-13]、塑性行為[14]以及堆積特性[15-16]等有顯著影響.深入研究顆粒形狀對顆粒材料物理力學特性的影響,剖析這些影響的細觀機理,對顆粒材料的理論研究和工程應(yīng)用有重要意義.受試驗技術(shù)限制,目前只能通過光彈試驗技術(shù)[17]、X-ray 斷層掃描技術(shù)[18]等研究顆粒材料宏觀力學響應(yīng)的細觀機理.與此同時,細觀數(shù)值模擬方法在顆粒材料的宏細觀力學特性研究方面得到廣泛應(yīng)用,如離散單元法(discrete element method,DEM)[19]和連續(xù)離散耦合分析方法(combined finite and discrete element method,FDEM)[20-21].其中FDEM 中將顆粒離散為有限元網(wǎng)格,理論上可以考慮任意顆粒形狀,在模擬復(fù)雜顆粒形狀方面具有顯著優(yōu)勢.
雖然顆粒材料的結(jié)構(gòu)特性與其物理力學性質(zhì)的相關(guān)性已經(jīng)得到大量研究證明[22-26],但是具有復(fù)雜顆粒形狀和寬粒徑分布的顆粒材料的結(jié)構(gòu)表征,以及建立宏觀力學響應(yīng)和細觀結(jié)構(gòu)的聯(lián)系仍是一個重大挑戰(zhàn).基于該現(xiàn)狀,顆粒統(tǒng)計力學得到了快速發(fā)展,例如Edwards 等[27-28]建立了顆粒系統(tǒng)的體積系綜理論(Edwards' volume ensemble),用統(tǒng)計物理的框架來研究非平衡態(tài)的顆粒堆積問題.該理論強調(diào)了局部自由體積(孔隙)在描述靜態(tài)顆粒材料阻塞行為中的重要作用.在工程中,常采用孔隙率、孔隙比、體積分數(shù)等表征顆粒材料的密實程度,但是這些表征量僅能提供宏觀尺度的自由體積信息,忽視了顆粒尺度自由體積的形態(tài)、分布和演化.
已有研究表明,顆粒尺度的自由體積與顆粒材料的力學性能、變形特征等密切相關(guān)[29-30],所以有必要探索顆粒材料的自由體積的統(tǒng)計分布特性及其在剪切過程中的演化規(guī)律.將顆粒與其周圍的自由體積組成一個局部元胞(以Voronoi 元胞較為典型)用以定量分析顆粒材料自由體積特征是一種常見思路.不少學者研究了顆粒材料在特定狀態(tài)下的局部自由體積[31-32].但是顆粒材料在剪切過程中會產(chǎn)生顆粒重排列現(xiàn)象,顆粒材料的細觀結(jié)構(gòu)和自由體積會不斷演化.因此,探索剪切過程非球顆粒材料的局部自由體積的演變規(guī)律,可以為研究顆粒材料的宏細觀多尺度力學特性提供新的思路.
本文將選擇具有不同非球度的橢球顆粒體系為研究對象,采用連續(xù)離散耦合分析方法進行三軸剪切數(shù)值模擬,分析不同非球度的橢球顆粒試樣在剪切過程中的宏細觀力學響應(yīng),并基于Voronoi 元胞對橢球顆粒體系的局部自由體積進行定量分析,探討局部自由體積的分布特征和演化規(guī)律以及受顆粒形狀的影響.
FDEM 最早由Munjiza 等[20]提出,充分融合了有限元和離散元方法的優(yōu)勢.顆粒的應(yīng)力變形計算采用有限單元法,顆粒運動和接觸模型采用離散單元法.FDEM 中顆粒均由有限元網(wǎng)格構(gòu)成,相比其他不連續(xù)數(shù)值模擬方法而言,最吸引人的優(yōu)點是能夠考慮各種復(fù)雜形狀[33-36].
有限元網(wǎng)格不僅定義了顆粒的形狀和邊界,而且可以很方便進行接觸檢索以及引入接觸模型.在FDEM 中,罰函數(shù)法被用于顆粒間的接觸分析,假定接觸力偶相互侵入,根據(jù)重疊量在節(jié)點間產(chǎn)生分布式接觸力.切向力的最大值根據(jù)庫侖摩擦定律控制.本文采用線性接觸模型計算顆粒間的接觸力,基于兩個接觸顆粒的相對速度,通過阻尼力定律確定阻尼力.接觸力的計算公式如下
式中,和 分別為節(jié)點法向與切向接觸剛度,定義knkt為kn=pnAnode,kt=ptAnode,pn和pt分別為法向和切向罰剛度,Anode為接觸面上的節(jié)點控制區(qū)的面積,為滑移率,Δt為時間增量,δn為節(jié)點侵入量,βn和βt為法向與切向臨界阻尼系數(shù),ζn和 ζt為法向和切向 恢復(fù)系數(shù),m為節(jié)點質(zhì)量.
Barrett[37]將顆粒形狀量化分為形態(tài)、圓度和表面紋理3 個層面.本文重點關(guān)注形態(tài),故選取橢球顆粒體系作為研究對象.共生成7 種主軸長度不同的橢球顆粒,顆粒的各種形態(tài)參數(shù)如圖1 所示.其中L,I,S分別代表橢球長軸、中軸和短軸的長度.α 為Domokos系數(shù)[38],定義為α=(1/S+1/I+1/L)表征顆粒主軸各向異性,在橢球形顆粒中可以反映顆粒的非球度.ψ為球度,定義為表征顆粒形狀偏離圓球的程度,式中A和V分別代表顆粒的表面積和體積.后文分析中選取 α 作為橢球顆粒的形狀量化指標.
圖1 顆粒形狀描述Fig.1 Particle shape descriptors
在立方體容器中生成等效半徑(定義為等體積球體的半徑)在6.8~ 10.2 mm 范圍內(nèi)均勻分布的6394 個顆粒,顆粒位置和傾向均隨機以避免產(chǎn)生初始各向異性.采用各向同性壓縮制樣,制樣過程中顆粒間摩擦系數(shù)設(shè)為0.1,六面無摩擦的剛性墻體緩慢壓縮試樣直至達到目標圍壓0.5 MPa.圖2 為α=3.394的數(shù)值試樣,試樣尺寸為0.42 m × 0.42 m ×0.42 m.采用常規(guī)三軸應(yīng)力路徑加載,圍壓保持σ2=σ3=0.5 MPa,控制頂部和底部兩塊剛性墻體以恒定速度相向移動.
圖2 數(shù)值試樣Fig.2 Numerical sample
數(shù)值試驗的可靠度很大程度取決于參數(shù)取值.本試驗中顆粒密度、彈性模量和泊松比參考天然砂礫石的自然屬性.根據(jù)巖石材料的回彈性質(zhì)[39],臨界阻尼系數(shù)取0.03,對應(yīng)的恢復(fù)系數(shù)為0.9.Tatone 和Grasselli[40]驗證了當接觸罰剛度設(shè)為10 倍彈性模量時所得到的彈性響應(yīng)與室內(nèi)物理試驗匹配較好,故本試驗接觸罰剛度設(shè)為彈性模量的10 倍.通過顆粒柱坍塌物理試驗和數(shù)值試驗對比分析來率定顆粒摩擦系數(shù).具體步驟如圖3 所示:選取粒徑在6~ 10 mm均勻分布的砂礫石裝入無底板圓筒,隨后緩慢抽離圓筒,重復(fù)10 次試驗,測得休止角均值為29.19°.對物理試驗中的砂礫石進行掃描,通過主成分分析得到主軸長度,在此基礎(chǔ)上重構(gòu)出形狀相近的橢球顆粒.調(diào)整顆粒摩擦系數(shù)進行顆粒柱坍塌數(shù)值試驗,使得休止角逼近物理試驗的結(jié)果,最終摩擦系數(shù)設(shè)為0.5 時休止角為29.5°.FDEM 數(shù)值模擬所需細觀參數(shù)匯總于表1.
表1 FDEM 數(shù)值模擬細觀參數(shù)Table 1 Parameters used in the FDEM simulation
圖3 顆粒柱坍塌試驗Fig.3 Granular column collapse tests
圖4 和圖5 分別為7 組不同非球度的橢球顆粒試樣的偏應(yīng)力?軸向應(yīng)變和體積應(yīng)變?軸向應(yīng)變關(guān)系曲線.不同試樣的宏觀響應(yīng)呈現(xiàn)相似的規(guī)律,偏應(yīng)力在加載初期迅速上升達到峰值,隨后開始下降,呈現(xiàn)出明顯的應(yīng)變軟化現(xiàn)象.試樣在加載初期有少量體積壓縮,隨后發(fā)生明顯的體脹.在較大軸向應(yīng)變時偏應(yīng)力和體積應(yīng)變都趨于穩(wěn)定或呈緩慢變化趨勢.這種應(yīng)力和體積響應(yīng)在密實的無黏性顆粒材料中非常典型[6].隨形狀參數(shù) α 增大,即顆粒非球度增大,試樣的剪切強度明顯增大,并伴隨著更加顯著的剪脹.目前普遍認為這種影響源自非球顆粒的互鎖效應(yīng),其增強了顆粒抵抗轉(zhuǎn)動的能力[12,41].
圖4 不同非球度的橢球顆粒偏應(yīng)力?軸向應(yīng)變Fig.4 Curves of deviatoric stress versus axial strain for ellipsoidal particles with different asphericity
圖5 不同非球度的橢球顆粒體積應(yīng)變?軸向應(yīng)變Fig.5 Curves of volumetric strain versus axial strain for ellipsoidal particles with different asphericity
Voronoi 剖分作為一種基本幾何結(jié)構(gòu)劃分方法,將離散的顆粒在空間上劃分,利用所得到的空間結(jié)構(gòu)來表征材料的細觀結(jié)構(gòu),目前已成為結(jié)構(gòu)表征的一種常用手段.其中圓球顆粒的Voronoi 剖分被廣泛應(yīng)用,計算軟件也較為成熟(如Voro++[42]).相比而言,非球顆粒的Voronoi 剖分由于實現(xiàn)困難,在顆粒材料的計算分析中鮮有報道.為實現(xiàn)非球顆粒的Voronoi 剖分,Schaller 等[43]提出了Set Voronoi 剖分方法,現(xiàn)已被成功運用于非球顆粒體系的細觀結(jié)構(gòu)分析[44-45].具體步驟如圖6 所示:(1)在顆粒表面均勻分布足夠多的點,得到每個顆粒表面的離散點集;(2)計算所有顆粒表面離散點的Voronoi 元胞;(3)將屬于同一顆粒的Voronoi 元胞合并,形成該顆粒的Voronoi 元胞.
圖6 Set Voronoi 剖分二維示意圖Fig.6 Two-dimensional illustration of Set Voronoi tessellation
上述步驟適用于完全分離的非球顆粒體系.然而在FDEM 中,由于顆粒間會相互侵入,會造成局部計算失準.為了避免這種情況,在進行離散點的Voronoi 剖分之前,將每個顆粒表面離散點沿其法向方向收縮一定距離,可在不影響Voronoi 剖分結(jié)果的同時避免局部誤差[43].圖7 為圓球顆粒和橢球顆粒的Voronoi 元胞示意圖,可以看出圓球的Voronoi元胞是凸多面體,而橢球的Voronoi 元胞表面會存在凹面,形態(tài)更加復(fù)雜.
圖7 圓球和橢球顆粒的Voronoi 元胞Fig.7 Voronoi cells of spherical and ellipsoidal particles
對于顆粒材料而言,即使是圓球顆粒體系,在剪切過程中仍然會產(chǎn)生各向異性[46].由于各向異性與顆粒材料的力學性質(zhì)密切相關(guān),相關(guān)研究受到廣泛關(guān)注.例如Ouadfel 和Rothenburg[24]基于組構(gòu)和接觸力的各向異性張量,推導(dǎo)了(stress-force-fabric,SFF)關(guān)系,建立了組構(gòu)各向異性、接觸力各向異性與宏觀應(yīng)力比之間的關(guān)系.組構(gòu)和接觸力是對接觸狀態(tài)的描述,相比而言,Voronoi 元胞是對顆粒局部空間的描述,其形狀的不規(guī)則程度可以用來量化局部空間各向異性.由于Voronoi 元胞本質(zhì)上是一個多面體,在量化其形狀特性上可以借鑒于顆粒形狀的描述指標.引入最常見的球度,描述Voronoi 元胞形狀偏離圓球(各向同性)的程度,定義其為與元胞等體積球體的表面積與元胞表面積的比值
式中,Vc和Sc分別代表Voronoi 元胞的體積與表面積.
為了避免邊界效應(yīng),本文分析排除了在邊界2 倍平均粒徑范圍內(nèi)的顆粒.圖8 為不同非球度的橢球顆粒試樣在不同加載階段(εa為軸向應(yīng)變)的Voronoi 元胞球度 ψc的概率密度分布函數(shù)(probability density function,PDF).可以發(fā)現(xiàn) ψc呈現(xiàn)高斯分布,圖中實線為高斯分布函數(shù)擬合曲線.在加載過程中 ψc分布逐漸向左移動(均值降低),且在加載到較大應(yīng)變時趨于穩(wěn)定.這說明顆粒材料在剪切過程中,局部空間的各向異性會逐漸增強直至系統(tǒng)達到臨界狀態(tài),且隨顆粒非球度增大各向異性增強越明顯,這與組構(gòu)各向異性和接觸力各向異性的演化類似[41,47].在剪切過程中,ψc分布的標準差也逐漸增大.
圖8 不同非球度的橢球顆粒在不同加載階段的Voronoi 元胞球度概率密度分布Fig.8 Probability density distributions of sphericity of Voronoi cells at the different loading states for ellipsoidal particles with different asphericity
當顆粒材料充分剪切后,會達到一個與初始狀態(tài)無關(guān)的力學穩(wěn)定態(tài),即臨界狀態(tài).顆粒材料在該狀態(tài)時的力學性質(zhì)受到廣泛關(guān)注[7,11].本文試樣雖未完全達到理想的臨界狀態(tài),但是在加載后期應(yīng)力狀態(tài)和體積已經(jīng)趨近于穩(wěn)定,故選取軸向應(yīng)變40%時為試樣的“臨界狀態(tài)”,代表充分剪切后的狀態(tài).圖9 為不同非球度的橢球顆粒試樣在臨界狀態(tài)的Voronoi元胞球度概率密度分布.隨形狀參數(shù) α 增大 ψc分布向左平移,分布更廣.圖10 展示了臨界狀態(tài)下Voronoi元胞球度與顆粒形狀參數(shù) α 的線性關(guān)系.這說明局部各向異性的分布強烈依賴顆粒形狀.相比圓球,橢球顆??梢援a(chǎn)生更為復(fù)雜的局部結(jié)構(gòu).
圖9 不同非球度的橢球顆粒在臨界狀態(tài)的Voronoi 元胞球度概率密度分布Fig.9 Probability density distributions of sphericity of Voronoi cells at critical state for ellipsoidal particles with different asphericity
為了量化顆粒尺度自由體積的大小,引入局部孔隙 比eloc=(Vc?Vp)/Vp,其中Vp和Vc分別為顆粒體積和其對應(yīng)的Voronoi 元胞體積.圖8 展示了不同非球度的橢球顆粒試樣在剪切過程中局部孔隙比eloc的概率密度分布.局部孔隙比eloc呈現(xiàn)類高斯分布.針對局部自由體積的統(tǒng)計分布特性已有一定報道,Schaller 等[45]發(fā)現(xiàn)扁橢球顆粒局部體積分數(shù)1/(1+eloc)近似服從高斯分布,Dong 等[32]發(fā)現(xiàn)橢球顆粒和柱狀顆粒局部孔隙比服從對數(shù)正態(tài)分布,Zhao 等[30]認為局部孔隙率eloc/(1+eloc)服從一個修改的對數(shù)正態(tài)分布.此外,還有一些研究發(fā)現(xiàn)圓球顆粒堆積體的局部體積分數(shù)的倒數(shù)(1+eloc) 服從k?Γ 分布[48-49].上述研究中以k?Γ 分布的擬合參數(shù)最少,其函數(shù)為
式中,k為表示分布函數(shù)形狀的參數(shù),Γ 為伽馬函數(shù),emin為當前堆積狀態(tài)下最小局部孔隙比,loc為局部孔隙比的均值.
圖11 中實線為k?Γ 函數(shù)擬合曲線.雖然k?Γ 分布是針對圓球顆粒體系提出的,但是對于橢球顆粒體系仍然適用.在加載過程中eloc分布逐漸向右移動,峰值減小并且分布范圍更寬,這意味著顆粒在剪切過程中局部自由體積變大且分布更廣.在加載后期這種演變逐漸趨于穩(wěn)定,與試樣先剪脹后達到臨界狀態(tài)的現(xiàn)象相吻合.隨非球度增大這種演變越明顯,這與圖5 中剪脹隨 α 增強相對應(yīng).
圖11 不同非球度的橢球顆粒在不同加載階段的局部孔隙比概率密度分布Fig.11 Probability density distributions of local void ratio at the different loading states for ellipsoidal particles with different asphericity
圖12 為全局孔隙比eglo=0.638時不同非球度的橢球顆粒試樣的局部孔隙比的概率密度分布,此時7 組試樣剪切狀態(tài)不同,軸向應(yīng)變分別為9.0%,12.7%,15.2%,21.0%,27.1%,31.8%和37.7%.可以發(fā)現(xiàn)不同非球度試樣的局部孔隙比分布基本一致,這說明顆粒局部孔隙比的分布僅與全局孔隙比相關(guān),不受顆粒形態(tài)和剪切狀態(tài)的影響.圖13 為不同非球度的橢球顆粒試樣的局部孔隙比分布標準差σ(eloc) 與全局孔隙比eglo的關(guān)系.σ (eloc)與eglo呈現(xiàn)明顯的線性關(guān)系,且這個線性關(guān)系不受顆粒形態(tài)影響,進一步證明了顆粒局部孔隙比的分布僅受全局孔隙比控制.
圖12 不同非球度的橢球顆粒在全局孔隙比相同時局部孔隙比概率密度分布Fig.12 Probability density distributions of local void ratio at the states as same global void ratio for ellipsoidal particles with different asphericity
圖13 局部孔隙比的標準差與全局孔隙比的關(guān)系Fig.13 Relationship between standard deviation of local void ratio and global void ratio
從局部孔隙比分布的演化可知剪切過程中局部自由體積的變化是非均勻的,探索局部孔隙比波動可以幫助了解顆粒體系剪脹的細觀機理.選取了5 個加載階段(εa=4%,10%,20%,30%和40%),分析不同非球度的橢球顆粒試樣在應(yīng)變窗口 Δ εa≈0.8%的局部孔隙比波動.圖14 為不同非球度的橢球顆粒試樣在不同加載階段的局部孔隙比波動的概率密度分布.局部孔隙比的波動呈現(xiàn)非對稱拉普拉斯分布(asymmetric laplace distribution,ALD),其函數(shù)為
圖14 不同非球度的橢球顆粒試樣在不同加載階段的局部孔隙比波動概率密度分布Fig.14 Probability density distributions of local void ratio fluctuations at the different loading states for particles with different shapes
式中,m為分界點,λ 為尺度參數(shù),κ 為非對稱參數(shù).
Δeloc的分布可以分為左右兩個非對稱的指數(shù)函數(shù)(在單對數(shù)坐標系下表現(xiàn)為斜率不一致的兩段線性分布),分界點在 Δeloc=0附近.左右兩側(cè)分布的非對稱性刻畫了局部自由體積收縮和膨脹的博弈.在加載初期(εa=4%)體積膨脹明顯處于優(yōu)勢地位(圖中圖14(a)斜率絕對值小于圖14(b)),試樣表現(xiàn)出明顯的剪脹行為.隨形狀參數(shù) α 增大,兩側(cè)斜率差異增大,這意味著有更顯著的體積膨脹,與宏觀體積響應(yīng)相符合(圖5).隨著剪切進行,代表體積膨脹的圖14(a)體積波動分布變化較小,而圖14(b)分布的斜率有明顯下降,兩側(cè)趨于對稱,即剪脹會隨剪切過程逐漸停止.
ALD 中非對稱參數(shù) κ 為分布兩側(cè)斜率比值的平方,描述了ALD 的非對稱程度,也刻畫了體積膨脹的強度.κ 越小體積膨脹越明顯,κ=1時體積無變化.如圖15 所示,κ 與全局孔隙比eglo呈現(xiàn)明顯的線性關(guān)系,并在 κ=1附近截斷.這意味著體積膨脹的強度很大程度上取決于試樣當前全局孔隙與臨界全局孔隙比的差值,故隨 α 增大體脹更顯著可以歸因于橢球顆粒在相同固結(jié)條件下可以形成更密實的堆積.
圖15 非對稱參數(shù) κ 與全局孔隙比的關(guān)系Fig.15 Relationship between asymmetric parameter κ and global void ratio
本文采用連續(xù)離散耦合分析方法,對具有不同非球度的橢球顆粒試樣進行了三軸剪切數(shù)值模擬.基于Set Voronoi 剖分技術(shù)對剪切過程中的顆粒試樣進行Voronoi 元胞分割,研究了顆粒試樣在剪切過程中局部自由體積的統(tǒng)計分布特性和演化規(guī)律,以及顆粒形態(tài)的影響.主要結(jié)論如下:
(1)不同橢球顆粒試樣的Voronoi 元胞的球度在剪切過程中均服從高斯分布,但其均值降低,標準差略有上升,表明顆粒材料在剪切過程局部各向異性的逐漸增強.局部各向異性程度與顆粒形態(tài)顯著相關(guān),臨界狀態(tài)時Voronoi 元胞球度隨顆粒形態(tài)參數(shù) α線性減小.
(2) 不同橢球顆粒試樣的局部孔隙比均服從k?Γ 分布,且這個分布僅與全局孔隙比相關(guān),不受顆粒形態(tài)和剪切狀態(tài)的影響.非球度更大的顆粒在剪切過程中會產(chǎn)生更強烈的重排列,導(dǎo)致局部孔隙比和Voronoi 元胞球度在剪切過程中的變化隨顆粒形狀參數(shù) α 的增大而增大.
(3)不同非球度的橢球顆粒試樣的局部孔隙比波動均服從非對稱拉普拉斯分布,這種不對稱性刻畫了局部自由體積收縮和膨脹的博弈.非對稱參數(shù)與全局孔隙比呈現(xiàn)明顯的線性關(guān)系,表明體積膨脹的強度很大程度上取決于試樣當前全局孔隙此與臨界全局孔隙比的差值.