不同細料含量土石混合料塑性行為離散元模擬1)

王濤 朱俊高,2) 劉斯宏

* (河海大學巖土力學與堤壩工程教育部重點實驗室,南京 210098)

? (河海大學水利水電學院,南京 210098)

引言

土石混合料是指由大粒徑塊石和作為填充成分的細粒土組成的混合料[1-2].土石混合料在自然界中分布廣泛,可作為基礎填筑材料應用于水利、鐵路、公路等工程建設之中.作為典型的二元混合顆粒材料,土石混合料的力學特性與細料含量(fc)密切相關.細料含量較低時,土石混合料由石料承擔骨架,其力學特性主要由石顆粒間的摩擦特性決定;細料含量超過閾值細料含量時,土石混合料由土料承擔骨架,石顆粒懸浮在土料骨架基質中,其力學特性主要由土顆粒間的摩擦特性決定[3].

目前,國內外學者針對土石混合料已開展了大量的室內試驗[4-6]、現場試驗[7-8]與數值模擬[9-11]研究,發現土石混合料的抗剪強度、剪脹特性、臨界狀態線與破壞模式等均與細料含量密切相關.然而,有關細料含量對土石混合料塑性行為的影響卻鮮有研究,從細觀層面闡明細料含量對土石混合料塑性行為的影響機制,對建立合理的彈塑性本構模型具有重要意義.

為了模擬土石混合料的應力應變關系,近年來眾多唯象類本構模型相繼被提出[12-13],這些本構模型是在總結試驗規律基礎上提出的,能夠較好地擬合試驗結果,但是缺乏物理機制的支撐.土石混合料是離散顆粒組成的集合體,其宏觀應力應變由集合體內部接觸的幾何分布與接觸力大小決定[14-15],但是在唯象類本構模型中其離散的本質并未被顯性表達,在此層面上,離散單元法由于易獲得代表體積單元中顆粒位置和接觸信息,是研究細觀參數如何影響宏觀特征的有效方法.顆粒集合體由無數個細觀結構單元組成,這些細觀結構的幾何與力學特性直接決定了集合體宏觀的強度變形特性[16].近年來,針對細觀結構的研究成果有效解釋了顆粒材料的破壞模式[17]、失穩特性[18]與剪切帶的形成過程[19].因此,探究加載過程中不同細料含量土石混合料內細觀結構的演化規律可揭示細料含量對土石混合料塑性行為的影響機制.

本文目的在于探究細料含量對土石混合料塑性行為的影響規律及其細觀機制.需要說明的是,細料含量對石料骨架結構與土料骨架結構塑性行為的影響機制不同,因此本文重點討論石料骨架結構土石混合料.本文基于二維離散元數值模擬,研究細料含量對石料骨架結構土石混合料(fc=0,0.05,0.1)屈服面、塑性流動方向與失穩破壞等特性的影響規律,并從細觀結構角度對數值模擬結果進行解釋.

1 DEM 數值試驗與分析方法

1.1 雙軸壓縮數值模擬

離散單元法最早是由Cundall 和Strack[20]提出,之后成功應用于巖土不連續介質的非線性力學分析中.離散單元法的基本思想是將不連續介質離散化為獨立的單元,通過對各個獨立結構之間的相互作用進行分析,進而研究整體的力學性質.離散單元法一般采用顯式時步算法,假設單個時步內顆粒速度和加速度不變,求解顆粒的運動方程和力-位移方程,并在下一步中更新顆粒、墻體和接觸的信息,如此不斷更新交替,直至整個計算結束.

本研究采用開源軟件Yade[21]進行計算,Yade主要應用于巖土工程領域,對土力學細觀機理、山體滑坡、地震等工程問題有很好的支持,具有計算效率高的優點.

法向與切向接觸力與接觸位移可由下式確定

式中,kn為法向剛度,kt為切向剛度,δn為法向重疊量,δt為切向重疊量,?為顆粒間摩擦角.

法向剛度kn可由接觸顆粒的尺寸和材料模量E決定,并假定切向剛度與法向剛度的比值取定值r,kn與kt可由下式確定

式中,Rp和Rq分別為接觸顆粒p和q的半徑.

本次研究所采用的模型參數見表1,表1 的參數參照文獻[17-19]進行選取,其中設置顆粒-墻體摩擦角為0°可確保墻體的應力為主應力(上下墻體為大主應力,左右墻體為小主應力).

表1 DEM 模型參數取值Table 1 Parameters for DEM tests

為確保數值試樣為石料骨架結構,制備低細料含量試樣(fc=0,0.05,0.1),采用間斷級配,土顆粒在0.4～0.8 mm 粒徑范圍內均勻分布,石顆粒在4～8 mm 粒徑范圍內均勻分布.試樣制備分三個步驟.(1)生成試樣:首先在長方形墻體內生成相互無接觸的顆粒群(見圖1),將顆粒球心的z方向坐標固定以形成二維情況.(2)固結:首先采用內部加載法,通過增大顆粒粒徑使墻體處達到90 kPa 的圍壓,隨后停止增大顆粒粒徑,利用墻體的移動來施加設計的圍壓(100 kPa 或200 kPa).(3)剪切:保持側面墻體的圍壓不變,以足夠小恒定的速率移動頂部和底部墻體以剪切試樣,以保證剪切過程為準穩態[22].

圖1 DEM 雙軸壓縮模型Fig.1 DEM biaxial model

本研究采用相對密實度指標,保持不同細料含量的試樣初始相對密實度相同(Dr=10%),相對密實度的計算公式為

式中emax和emin分別為試樣的最大與最小孔隙比.

在DEM 模擬中,通常通過設置固結過程中不同顆粒-顆粒摩擦角控制試樣的密度.本研究將顆粒間摩擦角設為0°獲得試樣的最小孔隙比,這與文獻[23-24]中所采用的方法相同.對于最大孔隙比,文獻[25-26]研究發現,設置顆粒間摩擦角為35°可獲得試樣最大孔隙比,因此本研究通過設置顆粒間摩擦角為35°制備最大孔隙比試樣.獲取不同細粒含量試樣的最大、最小孔隙比后,可在固結過程中不斷試算顆粒間摩擦角直至其相對密實度達到10%.不同細料含量土石混合料試樣的最大孔隙比、最小孔隙比與制樣孔隙比列于表2.

表2 不同細料含量試樣最大最小孔隙比、制樣孔隙比與相對密實度Table 2 Minimum and maximum void ratios,void ratio for sample preparing and relative density for samples with different fines contents

圖2 給出了fc=0 試樣剪切過程中應力比η(η=q/p)與體積應變εv的變化曲線.可見,η與εv均隨軸向應變增大而增大,隨后達到某一恒定值,試樣達到臨界狀態,體現出典型的松散試樣的特性.加載初期,應力比急速增長,隨后波動上升,在波動上升階段,準穩態假設在某些加載步中可能不再成立.為了研究試樣在不同應力狀態下的穩定性,在雙軸壓縮過程中保存處于不同應力比的試樣(如圖2(a)中藍色實心點),這些應力狀態對應的η∈{0.15,0.22,0.33,0.43,0.46,0.49,0.52}.

圖2 雙軸壓縮過程中應力比η 與體積應變εv 隨軸向應變變化曲線(fc=0)Fig.2 Evoluitons of stress ratio η and volumetric strain εv with axial strain εyy (fc=0)

1.2 分岔失穩點分析

作為典型的非關聯流動材料,巖土材料在達到塑性破壞極限之前便可能發生失穩破壞,通常采用二階功失穩準則來描述失穩破壞.二階功失穩準則最早由Hill[27]提出,并被廣泛應用于巖土工程計算中.Darve等[28-29]詳細描述了顆粒材料分岔失穩特性,并建立了顆粒材料失穩分析的框架.

二階功失穩準則是指,若存在一個加載路徑,在此應力路徑下二階功W2=dσ:dε 為負,在該應力路徑下土體失穩.dσ 和 dε 分別為應力增量矢量和應變增量矢量.

利用二階功失穩準則,采用經典的應力控制的方向加載法研究不同應力狀態下試樣的穩定性.方向加載分為兩步:預穩定和方向加壓.

(1)預穩定:1.1 節中試樣是在加載過程中保存的,需要讓其在恒定的應力下達到穩定,穩定的判別標準是試樣的不平衡力Funb＜ 10-5.不平衡力是判斷試樣是否平衡的指標,其計算公式為

式中,Nc為顆粒Np的接觸總數;Fcp為顆粒Np所受到的合外力,Fc為顆粒Np所受所有接觸力的平均值.

1.1 節中保存的試樣在預穩定過程中產生的體積應變如圖2(b)所示,藍色實心點為預穩定前的體積應變,藍色空心點為預穩定后的體積應變.可見,當η ≤0.46 時,預穩定過程中的體積應變很小,可忽略;當 η=0.49,0.52,由于顆粒重排布預穩定過程中產生了不可忽略的體積應變.

(2)方向加壓:預穩定結束后,施加不同方向的應力增量 dσ,如圖3(a) 所示,應力增量的數值為,加載方向角為α.可根據式(5)計算施加在側向墻和豎直墻上的應力 dσxx與 dσyy.當圍壓為100 kPa 時,dσ 值設置為5 kPa,當圍壓為200 kPa 時,由于在固結過程中試樣儲存了更多的彈性能,需要采用更大的 dσ 來釋放這些彈性能,因此當圍壓為200 kPa 時,dσ 值設置為10 kPa

圖3 應力增量dσ 及其應變響應增量dεFig.3 Stress probe dσ and incremental strain response dε

對于施加的應力增量dσ,可以得到相應的應變增量響應 dε,dε 的方向角為β,如圖3(b)所示.可由dσ和 dε 計算歸一化的二階功[30]

圖4為fc=0 的試樣在不同初始應力比下的歸一化二階功,圖4 中的紅色圓圈代表為0,曲線在紅圈內時,為負,反之,為正.可以看出對于所有的應力比,主要在第一和第三象限降低,且最小值出現在第三象限.圖4(a)是η≤0.46 的情況,即預穩定過程中的體積應變可以忽略的情況,可以看出當應力比較小時,對于所有的加載方向均不存在負的,即試樣在低應力比時是穩定的.隨著應力比增大,最小值逐漸降低,當η=0.43,0.46 時開始出現負的,產生失穩錐(所有產生負的的加載方向的集合).圖4(b)是η＞ 0.46 的情況,可見在η=0.49 時無負值,η=0.52 時存在負值,即隨著η=0.46 增大至0.49 時,試樣由不穩定狀態轉變為穩定狀態,出現這一現象的原因是η=0.49 的試樣在預穩定過程中產生了較大的體積變形,試樣的細觀結構因此產生了較大改變,即進行方向加載的試樣與前期加載過程中保存的試樣并不相同,因此研究失穩特性時,應重點關注預穩定過程中不產生較大體積變形的試樣.

圖4 在不同初始應力比下的歸一化二階功 (fc=0)Fig.4 Circular diagrams of the normalized second-order work at different stress ratio η (fc=0)

2 細料含量對土石混合料塑性行為的影響

2.1 失穩特性

圖5 給出了不同細料含量的土石混合料在相同應力比下(η=0.43)的歸一化二階功.可以看出,隨著細料含量增大,試樣的最小值逐漸增大.當fc為0 時,存在一個較寬的失穩錐;當fc為0.05 時,失穩錐變窄;當fc增大至0.1 時,失穩錐消失,試樣在所有加載方向上的均大于0.這說明對于石料骨架結構,增大細料含量有利于試樣的穩定.

圖5 不同細料含量的試樣(fc=0,0.05,0.1)的歸一化二階功(η=0.43)Fig.5 Circular diagrams of the normalized second-order work computed for the specimen with different fines content (fc=0,0.05,0.1)under the same stress ratio

如前所述,巖土材料在達到莫爾庫倫破壞極限之前就有可能發生破壞,即應力空間中出現分岔失穩區域[28-30].分岔失穩區域僅存在于非關聯流動材料,對于關聯流動材料,莫爾-庫倫極限線與分岔失穩區域的下邊界線重合,無分岔失穩區域.分岔區域越大,表示材料的非關聯流動性越強.

為了進一步探究細料含量對土石混合料分岔失穩區域的影響,圖6 給出了不同細料含量試樣莫爾-庫倫極限線與分岔失穩區域的下邊界線的夾角θ,兩條直線之間的區域即為分岔失穩區,可見隨著細料含量增加,θ逐漸降低,說明細料含量的增加會降低材料的非關聯流動特性.

關于細料含量增大會降低材料的非關聯流動特性這一發現,將在2.2 節中重點討論.

2.2 非關聯流動性

在彈塑性力學框架中,應變增量可寫成彈性應變增量(dεe)和塑性應變增量(dεp)之和

為了獲得彈性應變,可在加載過程中設置顆粒間摩擦角為90°,此時顆粒間不會發生滑移,產生的應變全部為彈性應變.獲得彈性應變增量之后,可根據式(7)計算塑性應變增量.

圖7 給出了fc=0.1 的試樣在應力比為0.43 時不同加載方向下的應變增量.可見彈性應變增量數據點形成了一個中心位于坐標原點的橢圓,表現出純彈性行為[31].塑性應變增量數據點形成一條直線,表明塑性應變的方向與加載方向無關,證明了流動法則的存在.

圖8 給出了不同應力加載方向角α下的塑性應變增量方向角βp及其數值大小,可見βp保持不變,根據流動法則,塑性應變方向與加載方向無關且與該應力狀態下的塑性勢面法向一致,因此βp即為塑性勢面法方向.在彈塑性力學框架中,當應力在屈服面以內變化,只產生彈性應變,應力若超過屈服面,將產生塑性應變.當α ∈[60,240]時,試樣產生塑性應變,可知α=60和α=240為屈服面法方向的兩個切方向,因此屈服面的法方向為=150.

圖9 比較了不同細料含量試樣塑性應變大小,可見在彈性區(0≤α≤60與240≤α≤360),不產生塑性變形;在塑性區(60≤α≤240),產生塑性變形.由于W2=dσ:dε=dσ:dεe+dσ:dεp,且dσ:dεe≥0,因此塑性應變很大程度上決定了二階功的正負.由圖9 可見增大細料含量可以減小塑性應變大小,這也解釋了為什么增大細料含量有利于試樣的穩定.

圖9 不同細料含量試樣(fc=0,0.05,0.1)在相同應力比下(η=0.43)的塑性應變大小Fig.9 Norms of plastic strain for specimens with fc=0,0.05,0.1 at the same η

圖10 給出了屈服面f與塑性勢面g的示意圖,及其法方向n和m隨細料含量的變化情況.可見,隨著細料含量增加,屈服面法方向n保持不變,而塑性勢面法方向m逐漸靠近n,m與n之間夾角減小,試樣的非關聯流動特性降低,這與2.1 節中的結果一致.圖11 為有細粒和無細粒情況下塑性變形流動方向的示意圖.集合體變形主要是來自力鏈在載荷作用下崩塌引起的顆粒重新排布,圖11(a)所示,初始狀態下,力鏈的方向與加載方向(y方向)一致,一旦力鏈發生崩塌豎向會產生較大的豎向變形,同時水平向會產生一定的膨脹,需要注意的是豎直方向的變形遠大于水平向的變形.當試樣中含有細顆粒時,如圖11(b)所示,由于細顆粒會側向支撐力鏈,因此試樣豎向的變形相較于無細顆粒時更小,所以加入細顆粒后會使塑性流動方向發生逆時針的旋轉,這與圖10 所展示的m方向的變化一致.

圖10 屈服面f 與塑性勢面g 的示意圖及其法方向隨細料含量的變化情況(η=0.43)Fig.10 Schematic diagram of yield surface f and plastic potential surface g (on the left) and their normal directions m and n (on the right)for specimens with different fc at the same η=0.43

圖11 有細粒和無細粒情況下流動方向示意圖Fig.11 Illustration of the plastic flow direction with and without fine grain

3 細觀機理

3.1 基本概念回顧

(1)各向異性

各向異性是顆粒材料的一個重要特性[32],可分為幾何各向異性和力學各向異性[33].幾何各向異性可用標量ac定量描述

(2)力鏈

力鏈是顆粒集合體最基本的細觀結構,力鏈在顆粒結合體中負責傳遞較大的接觸力并具有準直線的幾何特征.Peters等[36]對產生力鏈的三個條件歸納如下:(1) 力鏈上的顆粒的最大主應力大于所有顆粒的平均最大主應力;(2) 力鏈顆粒與其相鄰顆粒的圓心連線和最大主應力方向夾角小于45°;(3) 一條力鏈至少含有3 個顆粒.

力鏈的屈曲會使集合體承載力降低,甚至導致試樣發生崩塌[37].如圖12 所示,考慮3 個顆粒組成的力鏈,力鏈的偏角θ定義為中心顆粒與兩端顆粒圓心連線的夾角.假設在t時刻力鏈偏角為 θt,經過時間 dt后,力鏈偏角為θt+dt,則力鏈屈曲角θb為θb=θt-θt+dt.需要定義一個閾值θc,當θb＜θc時,認為力鏈未發生屈曲,反之,力鏈發生屈曲[37].本研究中,作者經過反復試算,發現閾值取為1 度時試驗規律較好,因此選取θc=1°.

圖12 力鏈偏角示意圖Fig.12 Illustration of force chain bending angle

3.2 細料含量對內部幾何結構與力學狀態的影響

為探究細料含量對試樣幾何與力學各向異性的影響,圖13 繪制了不同細料含量試樣在雙軸壓縮過程中ac(幾何各向異性)和an(力學各向異性)隨應力比的變化過程.可見,不同細料含量試樣在相同應力比時an相同,但ac不同,這意味著在相同應力比時不同細粒含量試樣具有相似的力學狀態,但內部幾何結構不同.

圖13 不同細料含量試樣 ac (幾何各向異性)和 an (力學各向異性)隨應力比的變化過程Fig.13 Evolution of ac (filled symbols) and an (empty symbols) with respect to η in specimens with different fc during biaxial tests

在外載荷作用下,試樣的塑性變形主要取決于顆粒間的滑移、滾動、原有接觸的打開與新接觸的生成,這些都與試樣內部的幾何結構相關,而不同細料含量試樣具有不同的內部幾何結構,因此細料含量會影響試樣的塑性流動方向(塑性勢面法向).這一結果與一些學者發現的剪脹與試樣內部幾何結構密切相關這一結論一致[38-39].然而,屈服面的方向與集合體的應力狀態相關,隨著細料含量增加,集合體力學各向異性不發生改變,因此細料含量不改變屈服面的方向.

材料的屈服對應著承載力的降低,細觀上表現為集合體內部的力鏈難以抵抗外載荷,發生屈曲變形,并伴隨著大量的接觸滑動.以下從力鏈屈曲、接觸滑動角度進一步解釋為何細料含量不改變屈服面方向.圖14 為不同細料含量試樣在雙軸壓縮過程中力鏈屈曲角的概率密度函數.可見,不同細料含量的試樣屈曲角的概率密度函數幾乎重合,這說明試樣的內部力學狀態基本相同,不受細料含量的影響.

圖14 不同細料含量試樣(fc=0,0.05,0.1)在雙軸壓縮過程中力鏈屈曲角的概率密度函數(應力比為0.43)Fig.14 Buckling angle probability density functions of specimens with fc=0,0.05,0.1(at stress ratio η=0.43)

為了進一步探究細料含量對試樣力學狀態的影響規律,定義接觸滑動因子指標Ip

式中,? 為顆粒-顆粒間摩擦角;Fn和Ft分別為法向和切向接觸力.Ip介于0 至1 之間,反映了某一接觸是否接近滑動,Ip=1 說明接觸滑動.

圖15 為不同細料含量試樣的所有接觸數Ip概率密度函數,可見,不同細料含量的試樣Ip具有相似的概率密度函數,這一結果進一步佐證了細顆粒的加入并不會很大程度影響顆粒集合體的內部力學狀態,因此,細顆粒的加入不影響屈服面的方向.

圖15 不同細料含量試樣的所有接觸數 Ip 概率密度函數(應力比為0.43)Fig.15 Sliding index probability density functions of specimens with fc=0,0.05,0.1 at the same stress ratio η=0.43

4 結論與展望

通過二維離散元數值模擬,研究了細料含量對石料骨架結構整體穩定性和塑性行為的影響規律,主要結論如下.

(1)細顆粒有利于顆粒集合體的整體穩定,細顆粒能夠限制集合體的塑性變形,提高試樣的穩定性.

(2)細顆?？刂祁w粒集合體塑性變形的方向,隨著細料含量增加,塑性勢面法方向和屈服面法方向之間的夾角減小,材料的非關聯流動性降低,材料分岔失穩區域變窄.

(3)盡管加入到石顆粒中的部分細顆粒與石顆粒共同承擔骨架作用,但是細顆粒的加入不影響顆粒集合體的力學狀態,不改變試樣的屈服面方向.

細觀分析結果表明不同細料含量的試樣在經歷相同加載路徑至同一應力比時,具有相似的力學狀態和不同的內部幾何結構.從實用性角度來說,這一發現對構建考慮細料含量影響的土石混合料彈塑性本構模型具有指導意義,例如對于石料骨架結構土石混合料,屈服函數應與細粒含量無關,塑性勢函數應為細粒含量相關,且保證材料非關聯流動性隨細粒含量增大而減弱.

本文主要研究目的在于揭示細粒含量影響土石混合料塑性行為的細觀機制,模擬時僅采用了常規二維球形顆粒,擬在后續研究中拓展至三維情況并考慮顆粒形狀的影響