王怡舒 沈超敏 劉斯宏,2) 陳靜濤

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

?(中國建筑第二工程局有限公司,上海 200135)

引言

顆粒材料是由大量相互接觸的離散顆粒組成的系統(tǒng)[1],其宏觀力學行為與內(nèi)部顆粒間復雜接觸密切相關,而基于連續(xù)介質力學的唯象描述主要建立在顆粒材料宏觀物理試驗的基礎上,不能追蹤顆粒間的相互作用.離散單元法(discrete element method,DEM)[2]作為研究顆粒材料內(nèi)部接觸結構(組構)特性的一種有效模擬方法,被廣泛用于探究顆粒材料細觀與宏觀力學特性的聯(lián)系[3-6].從細觀角度出發(fā),研究顆粒間的傳力機制[7],并建立宏細觀尺度間的聯(lián)系,是研究顆粒材料宏細觀力學性質的必要方法.

顆粒材料內(nèi)部力的傳遞通過顆粒間的接觸實現(xiàn),在外部載荷的作用下,當接觸點處的切向接觸力大于最大靜摩擦時,形成滑動接觸;另一方面,由于顆粒形狀的不規(guī)則性,在顆粒轉動過程中會產(chǎn)生抗轉動力矩,當顆粒點的轉動力矩大于最大抗轉動力矩時,形成轉動接觸.采用滑動摩擦和抗轉動力矩是反映顆粒材料的摩擦性和抗轉動特性的重要途徑,被廣泛應用于顆粒材料的理論模型和數(shù)值模擬中[8-11],如在DEM 模擬時,為了考慮顆粒的粗糙程度,顆粒的接觸通常設置Mohr-Coulomb 摩擦,為了考慮顆粒形狀的影響則在顆粒接觸模型中引入抗轉動力矩[12-13].

顆粒間接觸也可以根據(jù)法向接觸力的大小,采用平均法向接觸力的大小作為分界線,將體系的接觸分為大于平均法向接觸力的強接觸和小于平均法向接觸力的弱接觸,相互連接形成的強弱力鏈網(wǎng)絡抵抗外部載荷[14-15],該分類方法被很多學者采用[3,16-18].目前,一些學者發(fā)現(xiàn)滑動接觸對顆粒材料的宏觀力學行為影響深刻,如Aloso-Marroquin 等[19]認為滑動接觸對顆粒材料的塑性行為非常重要,可以將全局接觸劃分為非滑動-滑動接觸網(wǎng)絡,以便理解顆粒體系的宏觀力學響應.此外,研究發(fā)現(xiàn)強弱力鏈扮演的角色也有所差異,Radjai 等[14]指出強接觸對宏觀力學響應有著決定性的影響,強接觸承擔全部偏載荷,而幾乎所有的摩擦耗散發(fā)生在弱接觸.與滑動接觸形成較鮮明對比的是,盡管轉動接觸對真實顆粒的動力學影響顯著,目前針對轉動接觸力鏈的類似研究較少.

顆粒材料細觀尺度組構各向異性的產(chǎn)生及變化是其宏觀尺度各向異性力學特征的本質原因,組構各向異性分為兩類:原生各向異性和誘發(fā)各向異性,原生各向異性是材料的固有屬性,在巖土材料中典型表現(xiàn)為橫觀各向同性和沉積面鉛直方向的各向異性,由顆粒本身的性質(裂隙或者孔隙)以及沉積過程造成;誘發(fā)各向異性隨著應力狀態(tài)的改變而發(fā)生變化.各向異性可以通過顆粒接觸的組構張量定量表征,Oda[20]較早提出了組構張量的概念定量統(tǒng)計組構各向異性;在此基礎上,Rothenburg 等[21-22]首先推導了二維顆粒體系的誘發(fā)各向異性表達式,并提出了顆粒材料應力-接觸力-組構理論關系;Chantawarangul[23]推導了三維條件下的應力-接觸力-組構關系;Zhao和Guo[24]初步揭示了顆粒材料中組構各向異性與宏觀應變局部化的直接聯(lián)系.各向異性與顆粒內(nèi)部的接觸類型密切相關,Sufian 等[25]將顆粒體系中的接觸網(wǎng)絡分為非滑動-強接觸、非滑動-弱接觸、滑動-強接觸和滑動-弱接觸,發(fā)現(xiàn)非滑動-強接觸的組構各向異性參數(shù)對顆粒材料抗剪能力的貢獻最大,在不同應力路徑下偏應力比與非滑動-強接觸的組構各向異性參數(shù)具有唯一的線性關系.然而,目前關于抗轉動效應對不同子接觸網(wǎng)絡的組構各向異性演變規(guī)律的影響仍缺乏系統(tǒng)的研究.

本文根據(jù)剪切過程中顆粒接觸點的滑動與否、轉動與否和強弱力情況,將全局接觸網(wǎng)絡分為不同的子接觸網(wǎng)絡.基于不同抗轉動系數(shù)下顆粒材料的常規(guī)三軸剪切試驗數(shù)值模擬結果,分析剪切過程中不同子接觸網(wǎng)絡的組構張量的演變規(guī)律,并探究抗轉動系數(shù)對子接觸網(wǎng)絡的各向異性演變規(guī)律及對宏觀應力的貢獻程度的影響.

1 DEM 模擬試驗

1.1 顆粒接觸模型

顆粒間的相互作用可用接觸法向n、枝向量b、法向接觸力fn、切向接觸力ft等矢量進行描述.接觸法向指顆粒接觸點處切平面外法線方向,可用其方向余弦表征;枝向量是指兩接觸顆粒幾何中心的連線,b=(R1+R2)n,其中R1和R2分別為兩接觸顆粒的半徑;顆粒間的法向和切向接觸力分別指兩接觸顆粒沿接觸法向和在接觸切平面上的接觸力,如圖1(a)所示.在進行離散元數(shù)值模擬時,考慮顆?？罐D動效應的線性接觸模型[26]包括3 個部分:法向、切向和轉動接觸模型,如圖1(b)所示.

圖1 抗轉動-線性接觸模型Fig.1 Linear-based rolling resistance contact model

顆粒間的法向接觸力和切向接觸力分別用下式進行計算

式中,fn,kn和un分別為法向接觸力、法向接觸剛度和法向重疊量;Δfs,ks和Δus分別為切向接觸力增量、切向接觸剛度和切向相對位移.顆粒間法向接觸剛度kn和切向接觸剛度ks可根據(jù)下式計算[26]

式中,E為顆粒彈性模量,A和L分別為接觸點處的橫截面積和長度,A=πr2,r=min(R1,R2),L=R1+R2;κ*為剛度比.當計算的切向接觸力fs大于最大靜摩擦μsfn時,顆粒間發(fā)生滑動,此時的顆粒接觸為滑動接觸.滑動接觸點的切向接觸力fs=μsfn(μs為顆粒間滑動摩擦系數(shù)).

轉動接觸模型[27-28]主要是通過在顆粒間加入一個抗轉動力矩mr,其定義為

式中,θr為顆粒間相對轉角;kr為顆粒轉動剛度,與切向接觸剛度ks相關,可按下式計算

當顆粒點的轉動力矩mr大于最大抗轉動力矩時,顆粒間發(fā)生相對轉動,此時的接觸為轉動接觸.轉動接觸點的抗轉動力矩mr=為顆粒間的抗轉動系數(shù)).

1.2 數(shù)值模擬方案

本文離散單元模擬的材料為粗顆粒材料,其初始尺寸為0.3 m×0.3 m×0.6 m(長×寬×高),根據(jù)試樣尺寸與顆粒的最大粒徑比不應小于5 的原則[29],控制最大顆粒粒徑dmax為60 mm.考慮DEM 的計算效率,試樣最小粒徑dmin設定為10 mm.有關研究表明粗粒料的級配可以用分形分布的函數(shù)來描述[30]

式中,d為顆粒的粒徑;dmax為最大粒徑;M(L＜d)為小于粒徑d的土體質量;MT為土體總質量;D為分形維數(shù).對于最大粒徑在60 mm 以內(nèi)的粗粒料,D=1.88～2.63 為級配良好區(qū)間[30],因此本文數(shù)值試樣級配的分形維數(shù)D=1.9.

在制備試樣時,首先在模型尺寸范圍內(nèi)隨機生成粒徑服從圖2(a) 所示分形分布的球體顆粒,數(shù)量為13 584,然后采用等向壓縮的方法固結至目標圍壓,在等向壓縮的過程中可以通過控制顆粒間滑動摩擦系數(shù)μs生成不同初始孔隙比的試樣.圖2(b)給出了在不同滑動摩擦系數(shù)μs下等向固結至200 kPa時試樣孔隙比e的變化規(guī)律,發(fā)現(xiàn)孔隙比e隨著μs的增大而增大,且增大幅度逐漸減小.在試樣從固結狀態(tài)的不同μs值變成μs=0.3 的穩(wěn)定狀態(tài)中,對于固結過程采用μs＞0.3 的試樣其孔隙比會降低,而對于固結過程采用μs＜0.3 的試樣其孔隙比保持不變.這種孔隙比降低的現(xiàn)象是由于減小滑動摩擦系數(shù)導致了顆粒間容許最大切向接觸力μsfn的減小,造成了更多顆粒接觸間的滑移,突然釋放的剪切彈簧中的應變能轉化為動能驅動顆粒重排,最終形成一個新的穩(wěn)定結構.

圖2 DEM 數(shù)值試樣級配與初始孔隙比Fig.2 Particle size distribution and initial void ratio of DEM specimen

本文在等向壓縮至預定圍壓值200 kPa 的過程中設定μs=1.0,μr=0.0,采用周期性邊界的伺服機制,制備得到相對密實度較低、邊界效應較小的試樣.本文的主要目的是研究剪切過程中的宏細觀力學特性而非重現(xiàn)室內(nèi)試驗過程,因此在保證宏觀力學規(guī)律合理的前提下,通過試算得到模擬中細觀參數(shù)取值見表1,其中彈性模量與Guo[31]和Yan 等[32]的取值在同一數(shù)量級,剛度比取為1.0[33-34],在準靜態(tài)加載情況下局部阻尼系數(shù)取為0.0[35].粒間抗轉動系數(shù)μr取為0.0,0.01,0.05,0.1,0.15,0.2,共進行6 組數(shù)值模擬試驗.

表1 剪切過程DEM 計算參數(shù)Table 1 Input DEM simulation parameters

2 顆粒體系的接觸網(wǎng)絡及表征方法

根據(jù)顆粒間的接觸特性,如接觸點的滑動與否、轉動與否和強弱力情況,可以定義不同的子接觸網(wǎng)絡k,其中k可以代表非滑動(ns)、滑動(s) 或非滑動-非滾動(ns-nr) 子接觸網(wǎng)絡等.圖3 為摩擦系數(shù)μs=0.3,μr=0.05 的試樣在軸向應變ε1=25.0%時顆粒體系的全局接觸網(wǎng)絡和不同子接觸網(wǎng)絡的力鏈圖.

圖3 顆粒接觸網(wǎng)絡分類方法Fig.3 Classification of contact networks(CNs)

為了定量表征不同子接觸網(wǎng)絡k的各向異性,一般用接觸法向n和枝向量b的分布特征描述顆粒材料的幾何各向異性,法向接觸力fn和切向接觸力ft的分布特征描述顆粒材料的力學各向異性.Oda[20]建議采用接觸法向來表征全局組構張量φi j,定義為

子接觸網(wǎng)絡k的各向異性參數(shù)可以根據(jù)接觸法向、枝向量、法向接觸力和切向接觸力的各向異性二階張量表示

表征子接觸網(wǎng)絡k中枝向量、法向接觸力和切向接觸力各向異性程度的二階偏張量可分別由下式進行計算

3 組構各向異性演變規(guī)律

3.1 子接觸網(wǎng)絡組構張量演變規(guī)律

圖4 給出了摩擦系數(shù)μs=0.3,μr=0.05 的試樣在剪切過程中全局組構張量φij和滑動、非滑動兩種子接觸網(wǎng)絡的組構張量變化規(guī)律,以及軸向應變ε1=25% 時顆粒接觸點數(shù)按接觸角的空間分布玫瑰圖.不同虛線表示全局組構張量φij在剪切過程中的變化,可以發(fā)現(xiàn)在ε1=0% 時,系統(tǒng)的組構φ11=φ22=φ33=0.33,表明試樣此時處于各向同性階段;隨著軸向應變增大,大主應力方向組構φ11先增大,之后隨著加載過程稍微減小并最終趨于一個相對穩(wěn)定的狀態(tài);相應地,中主應力方向φ22和小主應力方向φ33先減小,之后稍微增大并保持相對穩(wěn)定.對于滑動、非滑動兩種子接觸網(wǎng)絡,非滑動接觸的各向異性明顯,主要沿大主應力方向集中,滑動接觸在殘余狀態(tài)時各向異性較小,沿小主應力方向有輕微的集中.

圖4 滑動、非滑動接觸的組構張量演變規(guī)律Fig.4 Evolution of fabric tensor of non-sliding and sliding contact networks

考慮到顆粒的抗轉動特性,將接觸網(wǎng)絡進一步劃分為非滑動-非轉動、非滑動-轉動(ns-r)、滑動-非轉動(s-nr)和滑動-轉動(s-r)4 種子接觸網(wǎng)絡,圖5 給出了這4 種子接觸網(wǎng)絡的組構張量變化規(guī)律,發(fā)現(xiàn):

圖5 非滑動-滑動-非轉動-轉動接觸的組構張量演變規(guī)律Fig.5 Evolution of fabric tensor of the nonsliding-sliding-nonrolling-rolling contact neworks

(3)與轉動、非轉動接觸的相關組構張量的變化不是獨立的,受到顆粒間是否滑動的影響,剪切過程中系統(tǒng)的應力誘發(fā)各向異性主要是由非滑動接觸引起的.

研究發(fā)現(xiàn),接觸網(wǎng)絡中的強接觸和弱接觸對系統(tǒng)的力學作用機制完全不同,強接觸對宏觀力學貢獻有決定性的影響.圖6 給出了剪切過程中強、弱接觸網(wǎng)絡的組構張量演變規(guī)律,發(fā)現(xiàn)強接觸sf 的數(shù)量小于弱接觸wf 的數(shù)量,強接觸sf 表現(xiàn)出大于整體的各向異性,且沿豎向集中,弱接觸wf 在殘余狀態(tài)時各向異性不明顯,沿大主應力方向有很小的集中.該現(xiàn)象與文獻[3]的結論一致,證實了強接觸形成柱狀結構承擔外部載荷,并由較多的弱接觸沿四周支撐.

圖6 強、弱接觸的組構張量演變規(guī)律Fig.6 Evolution of fabric tensor of strong force and weak force contact networks

進一步地,當接觸網(wǎng)絡劃分為非滑動-強接觸(ns-sf)、非滑動-弱接觸(ns-wf)、滑動-強接觸(ssf)、滑動-弱接觸(s-wf) 時,4 種子接觸網(wǎng)絡的組構變化有較大地差異(圖7),發(fā)現(xiàn)非滑動-強接觸(圖7(a)) 的組構各向異性最明顯,非滑動-弱接觸(圖7(b)) 和滑動-強接觸(圖7(c)) 的組構與系統(tǒng)組構的大小相差不大,該3 種接觸均沿大主應力方向集中,而滑動-弱接觸(圖7(d))沿小主應力方向有輕微的集中.

圖7 非滑動-滑動-強-弱接觸的組構張量演變規(guī)律Fig.7 Evolution of fabric tensor of nonsliding-sliding-strong-weak contact networks

3.2 抗轉動特性對各向異性參數(shù)的影響

基于不同子接觸網(wǎng)絡的各向異性指標可以用來分析顆粒體系的細觀力學機理.圖8 分別給出了摩擦系數(shù)μs=0.3,μr=0.05 的試樣在剪切過程中不同子接觸網(wǎng)絡k的接觸法向各向異性參數(shù)和法向接觸力各向異性參數(shù)的演變規(guī)律,從圖8(a)可以看出非滑動、強接觸和非滑動-強接觸網(wǎng)絡的均大于全局接觸法向各向異性ac,其中強接觸(特別是非滑動-強接觸)是幾何各向異性的主要貢獻者,該規(guī)律與本文第2 節(jié)組構張量的演變規(guī)律一致.圖8(b)顯示全局的接觸法向力各向異性an最大,這是因為子接觸網(wǎng)絡僅可以從有限的接觸法向力中進行選擇,該結果與文獻[25]的模擬結果一致.

圖8 剪切過程中子接觸網(wǎng)絡的各向異性演變規(guī)律Fig.8 Evolution of anisotropy indexes of sub-contact networks during shearing

圖9 選取了對各向異性貢獻較大的3 個子接觸網(wǎng)絡(ns,sf,ns-sf) 和全局接觸網(wǎng)絡,研究抗轉動效應對四種接觸網(wǎng)絡的各向異性參數(shù)的影響規(guī)律,發(fā)現(xiàn):

圖9 抗轉動系數(shù)對不同接觸網(wǎng)絡各向異性指標演變規(guī)律的影響Fig.9 Influence of rolling resistance on the evolution of anisotropy indexes of different contact networks

圖9 抗轉動系數(shù)對不同接觸網(wǎng)絡各向異性指標演變規(guī)律的影響(續(xù))Fig.9 Influence of rolling resistance on the evolution of anisotropy indexes of different contact networks(continued)

4 應力-接觸力-組構關系

Christonffersen 等[37]提出了基于全局接觸網(wǎng)絡的顆粒體系應力張量

式中,V為顆粒體系總體積;為接觸點c處的法向接觸力;為接觸點c處的枝向量;Nc為顆粒體系的總接觸點數(shù).系統(tǒng)的平均有效應力p,廣義剪應力q,不同方向主應變εi(i=1,2,3)和體應變εv分別由下式進行計算

式中,li0是剪切開始時試樣不同方向的初始長度,li是剪切過程中試樣不同方向的長度.

圖10 給出了剪切過程中應力應變變化規(guī)律,可以看出,應力比出現(xiàn)較為明顯的峰值,并隨抗轉動系數(shù)μr的增大而增大,對體應變的影響也較為明顯,剪脹性增大.

圖10 不同抗轉動系數(shù)下宏觀應力-應變關系曲線Fig.10 Curves of macro stress-strain relationship under different rolling friction coefficients

Ouadfel 等[22]與Li 等[38]針對全局接觸推導了應力-接觸力-組構(stress-force-fabric,SFF)關系,當進一步將全局接觸分為不同接觸網(wǎng)絡時,Sufian 等[25]推導了在若干子接觸網(wǎng)絡中的SFF 關系,全局偏應力比q/p*可以表示為子接觸體系中不同各向異性參數(shù)的函數(shù)

式中,ζk=Nk/N是子接觸網(wǎng)絡k中接觸數(shù)目Nk占整個體系接觸數(shù)N的比例,γk=是子接觸網(wǎng)絡的平均法向接觸力與體系平均法向接觸力之比.

圖11 給出了抗轉動效應對不同子接觸網(wǎng)絡ξk的影響規(guī)律,發(fā)現(xiàn)非滑動和非滑動強接觸網(wǎng)絡的ξns和ξns-sf隨著抗轉動系數(shù)μr的增大而減小,說明提高抗轉動效應會減小接觸數(shù)目和平均法向接觸力的貢獻程度;而強接觸網(wǎng)絡的ξsf不隨抗轉動系數(shù)的改變而改變,說明強接觸網(wǎng)絡對宏觀應力的貢獻程度在不同μr情況下均相同.

圖11 抗轉動系數(shù)對ξk 的影響規(guī)律Fig.11 Influence of rolling resistance on ξk

圖11 抗轉動系數(shù)對ξk 的影響規(guī)律(續(xù))Fig.11 Influence of rolling resistance on ξk (continued)

圖12 對比了由SFF 關系得到的ns-s-sf-wf 接觸網(wǎng)絡應力比q/p*與由式(19) 計算的應力比q/p,可以看出,對于采用不同抗轉動系數(shù)的顆粒體系,SFF的理論預測值與實際的應力比非常接近,從而表明了SFF 理論的適用性.

圖12 SFF 理論應力比驗證Fig.12 Validation of SFF relationship

5 結論

本文采用離散單元法模擬了不同抗轉動系數(shù)μr情況下顆粒材料三軸剪切試驗,根據(jù)接觸點的滑動與否、轉動與否和強弱力情況,將顆粒體系分為不同的子接觸網(wǎng)絡,分析了剪切過程中不同子接觸網(wǎng)絡的組構張量的演變規(guī)律,并通過統(tǒng)計不同子接觸網(wǎng)絡的各向異性參數(shù),研究了抗轉動效應對各向異性演變規(guī)律的影響,主要結論有:

(1) 分析了剪切過程中非滑動-滑動-非轉動-轉動4 種子接觸網(wǎng)絡的組構張量變化規(guī)律,發(fā)現(xiàn)轉動、非轉動接觸的組構變化不是獨立的,受到顆粒間是否滑動的影響,剪切過程中系統(tǒng)的應力誘發(fā)各向異性主要是由非滑動接觸引起的.

(2)剪切過程中非滑動接觸的各向異性明顯,主要沿大主應力方向集中,滑動接觸在殘余狀態(tài)時各向異性較小,沿小主應力方向有輕微的集中;強接觸表現(xiàn)出大于整體的各向異性,且沿豎向集中,弱接觸在殘余狀態(tài)時各向異性不明顯,沿大主應力方向有很小的集中.

(3)通過分析抗轉動效應對不同子接觸網(wǎng)絡的各向異性參數(shù)影響,發(fā)現(xiàn)提高抗轉動效應會導致組構各向異性的增大,并提高接觸網(wǎng)絡的穩(wěn)定性.非滑動接觸網(wǎng)絡和全局網(wǎng)絡的法向接觸力各向異性和an均隨抗轉動系數(shù)μr的增大而增大,且法向接觸力各向異性比接觸法向各向異性大;而與強接觸相關的子接觸網(wǎng)絡的法向接觸力各向異性和變化不大,且接觸法向各向異性占主導地位.

(4)非滑動和非滑動-強接觸網(wǎng)絡的ξns和ξns-sf均隨著抗轉動系數(shù)μr的增大而減小,說明提高抗轉動效應會減小接觸數(shù)目和平均法向接觸力對宏觀應力的貢獻程度;而強接觸網(wǎng)絡的ξsf不隨抗轉動系數(shù)的改變而改變,說明強接觸網(wǎng)絡對宏觀應力的貢獻程度在不同μr情況下均相同.