王 壯，段志波，廖新超，劉一鳴

(湖北工業(yè)大學(xué) 土木建筑與環(huán)境學(xué)院，湖北 武漢 430068)

砂土是人類進(jìn)行工程建設(shè)的主要材料之一且自然界中分布廣泛，工程實(shí)踐活動(dòng)都與砂土的力學(xué)特性息息相關(guān)：如建筑物地基承載力、路橋地基穩(wěn)定性和土質(zhì)邊坡穩(wěn)定性等問題[1]。直剪試驗(yàn)是研究土體抗剪強(qiáng)度及力學(xué)性能的重要方法，而室內(nèi)直剪試驗(yàn)通常只能研究砂土宏觀的力學(xué)特性，無法衡量細(xì)觀層面砂土顆粒之間相互作用對(duì)宏觀力學(xué)變形特性的影響。Cundall等[2]基于牛頓第二定律和力位移方程提出了針對(duì)砂土等理想散體材料為研究對(duì)象的離散單元法。隨著國內(nèi)外學(xué)者對(duì)離散元深化研究及計(jì)算機(jī)算力的巨大提升，離散單元法在研究土體細(xì)觀機(jī)理方面得到廣泛運(yùn)用。基于離散單元法，學(xué)者們做了大量模擬試驗(yàn)，并取得了許多成果：朱遙等[3]探究了顆粒形態(tài)對(duì)砂的抗剪強(qiáng)度的影響，發(fā)現(xiàn)不規(guī)則形態(tài)砂顆粒具有較大抗剪強(qiáng)度；潘遠(yuǎn)陽等[4]研究了剪切帶形成的細(xì)觀機(jī)理；孔亮等[5]分析了顆粒大小和形狀對(duì)砂土力學(xué)特性的影響，并構(gòu)建了顆粒外輪廓特征的性狀系數(shù)。文獻(xiàn)[6]表明當(dāng)顆粒間的摩擦較大時(shí)，決定顆粒材料力學(xué)機(jī)制及微觀變形機(jī)制的主要因素為顆粒之間相對(duì)轉(zhuǎn)動(dòng)作用。傳統(tǒng)的離散元模擬將顆粒設(shè)為球形或圓形，忽略了顆粒轉(zhuǎn)動(dòng)的影響。Iwashita等[7-8]均對(duì)離散單元法本構(gòu)模型進(jìn)行了改良，考慮了顆粒轉(zhuǎn)動(dòng)的影響，并引入了轉(zhuǎn)動(dòng)阻抗的概念。Jiang等[9]提出了抗轉(zhuǎn)動(dòng)線性接觸模型，基于三維接觸模型，設(shè)定顆粒間的接觸為面接觸，并考慮了滾動(dòng)阻力和扭轉(zhuǎn)阻力的影響。Zhu等[10]基于離散單元法，通過模擬三軸試驗(yàn)分析了剪切強(qiáng)度隨轉(zhuǎn)動(dòng)阻抗和細(xì)粒含量變化的機(jī)理。Zhao等[11]綜合比較了滾動(dòng)阻力和顆粒形狀對(duì)顆粒材料內(nèi)部剪切誘發(fā)的組構(gòu)和各向異性的影響，認(rèn)為滾動(dòng)阻力模型可有效再現(xiàn)砂土中剪切誘導(dǎo)各向異性主要特征。而對(duì)轉(zhuǎn)動(dòng)阻抗的選取則需綜合考慮顆粒的實(shí)際形狀和顆粒間的相對(duì)間距，總體而言對(duì)于使用純圓顆粒建模的試驗(yàn)，實(shí)際顆粒形狀越接近球形，相對(duì)間距越小，需選取的抗轉(zhuǎn)動(dòng)系數(shù)越小[12-13]。

本文利用離散元軟件PFC3D對(duì)砂樣直剪試驗(yàn)進(jìn)行一系列標(biāo)準(zhǔn)砂直剪數(shù)值模擬試驗(yàn)，考慮不同圍壓下轉(zhuǎn)動(dòng)阻抗對(duì)砂樣直剪過程中的宏觀受力及變形特性影響，細(xì)觀層面則通過接觸力鏈和配位數(shù)分析顆粒接觸狀態(tài)和剪切帶發(fā)育的演化規(guī)律；結(jié)合玫瑰圖、偏組構(gòu)及組構(gòu)張量分析各向異性演化規(guī)律。通過對(duì)砂土宏細(xì)觀力學(xué)結(jié)合分析，表明轉(zhuǎn)動(dòng)阻抗作用是一個(gè)不容忽視的因素，對(duì)揭示砂土剪切的本質(zhì)機(jī)理提供理論參考。

1 離散元建模

1.1 抗轉(zhuǎn)動(dòng)線性接觸模型

抗轉(zhuǎn)動(dòng)線性接觸模型在考慮顆粒轉(zhuǎn)動(dòng)影響時(shí)已經(jīng)取得良好運(yùn)用[10-11]。由于砂顆粒材料自身各向異性不明顯及不考慮砂顆粒與墻體之間的邊界效應(yīng)，綜合選取砂顆粒-墻體接觸類型為線性接觸模型，砂顆粒-砂顆粒接觸類型為抗轉(zhuǎn)動(dòng)線性接觸模型（圖1）。線性模型的接觸點(diǎn)可以想象為一對(duì)彈性彈簧，線性力由具有恒定法向剛度和切向剛度的線性彈簧產(chǎn)生，主要包括法向接觸力和切向接觸力，其力學(xué)原理為：

圖1 接觸模型Fig. 1 Contact model

在線性模型的基礎(chǔ)上，抗轉(zhuǎn)動(dòng)線性接觸模型[7]增加了轉(zhuǎn)動(dòng)阻抗力矩，其定義為：

式中：kr為 顆粒滾動(dòng)剛度；Rˉ為 有效接觸半徑；Ij為j顆粒慣性矩；rj為j顆粒半徑為與j顆粒接觸q顆粒的切向力為與j顆粒接觸q顆粒的扭矩；R1和R2分別為顆粒q、j的半徑。

式中：μr為 抗轉(zhuǎn)動(dòng)系數(shù)；θb為接觸顆粒間的相對(duì)轉(zhuǎn)動(dòng)為抗轉(zhuǎn)轉(zhuǎn)動(dòng)線性接觸模型中的法向接觸力。

1.2 數(shù)值模型及細(xì)觀參數(shù)

離散元建模過程中顆粒數(shù)量會(huì)影響模擬的結(jié)果和計(jì)算速率，在保證結(jié)果可靠前提下提高計(jì)算效率，通常采用擴(kuò)大顆粒粒徑R或者減小試樣尺寸L的方法，文獻(xiàn)[14]中數(shù)值平面應(yīng)變?cè)囼?yàn)結(jié)果表明，當(dāng)L/R=30～40 時(shí)才能保證尺寸效應(yīng)的影響最小。劉海濤等[15]通過一系列三軸試驗(yàn)的顆粒流模擬，認(rèn)為當(dāng)L/R＞40時(shí)，尺寸效應(yīng)對(duì)粗粒土抗剪強(qiáng)度沒有影響。尹小濤等[16]通過模擬不同內(nèi)尺度比巖土試樣單軸壓縮試驗(yàn)，發(fā)現(xiàn)試樣中的最小粒徑dmin與數(shù)值試樣的長度之比小于0.01時(shí)，數(shù)值試驗(yàn)的結(jié)果具有較好穩(wěn)定性。綜上所述，砂土離散元試樣采用尺寸與室內(nèi)大型直剪儀尺寸一致，尺寸為 305 mm×305 mm×200 mm，砂顆粒粒徑相對(duì)實(shí)際值放大10倍并刪除部分小粒徑，最終取最小粒徑dmin=1 mm，最大粒徑dmax=12 mm，這可同時(shí)滿足dmin/L≤0.01，L/R≥40。按初始孔隙率0.39隨機(jī)生成顆粒43 775個(gè)。離散元數(shù)值模型分別采用5塊剛性墻體模擬上下剪切盒，并生成2道水平剛性擋板防止剪切過程中的顆粒溢出（圖2）。

圖2 數(shù)值模型（單位：mm）Fig. 2 Numerical model of sand （unit：mm）

通過對(duì)比研究前人對(duì)砂土建模細(xì)觀參數(shù)的選取，初步確定砂細(xì)觀參數(shù)取值范圍，經(jīng)過反復(fù)調(diào)試最終確定細(xì)觀參數(shù)見表1。

表1 砂細(xì)觀參數(shù)Tab. 1 Sand parameters

砂土直剪數(shù)值模擬主要包括制樣、固結(jié)和剪切3個(gè)階段，本研究首先采用孔隙率制樣隨機(jī)生成顆粒并消除不平衡力影響，然后對(duì)試樣施加不同圍壓固結(jié)，最終剪切階段保持法向應(yīng)力不變，對(duì)下盒施加一個(gè)向右的剪切速度模擬砂土剪切。分析法向應(yīng)力分別為50、100、150和200 kPa時(shí)不同抗轉(zhuǎn)動(dòng)系數(shù)條件下砂土直剪的宏細(xì)觀力學(xué)特性。

2 宏觀力學(xué)特性

2.1 應(yīng)力-應(yīng)變響應(yīng)

由圖3中剪應(yīng)力-剪切位移曲線可知，曲線初始斜率基本一致，隨著抗轉(zhuǎn)動(dòng)系數(shù)的增大各曲線峰值隨之增大，而增幅逐漸降低；豎向位移-剪切位移曲線均表明試樣先微弱剪縮后剪脹的變形特性。相比于抗轉(zhuǎn)動(dòng)系數(shù)為0的試樣，其他試樣均表現(xiàn)出應(yīng)變軟化特征，抗轉(zhuǎn)動(dòng)系數(shù)越大，軟化越明顯，且在剪切后期仍有明顯剪脹發(fā)生。

圖3 剪切應(yīng)力-剪切位移及豎向位移-剪切位移曲線Fig. 3 Curves of shear stress-shear displacement and volume change-shear displacement

以上規(guī)律表明：轉(zhuǎn)動(dòng)阻抗對(duì)砂土試樣的初始彈性模量無明顯影響；對(duì)試樣的峰值強(qiáng)度有顯著增強(qiáng)作用，但隨著μr的增大，轉(zhuǎn)動(dòng)阻抗的影響逐漸趨于飽和；在轉(zhuǎn)動(dòng)阻抗的作用下，砂體更易表現(xiàn)出應(yīng)變軟化現(xiàn)象，且隨著μr的增大，軟化逐漸加劇。分析主要原因?yàn)椋河捎谵D(zhuǎn)動(dòng)阻抗的作用，顆粒間發(fā)生相對(duì)轉(zhuǎn)動(dòng)需要克服的阻力變大，造成宏觀抗剪強(qiáng)度的增大；而在轉(zhuǎn)動(dòng)阻抗作用下試樣表現(xiàn)出的應(yīng)變軟化現(xiàn)象符合中密砂的性狀。

2.2 抗剪強(qiáng)度指標(biāo)

將模擬結(jié)果抗剪強(qiáng)度包線與王軍等[20-23]所做室內(nèi)標(biāo)準(zhǔn)砂直剪試驗(yàn)擬合包線進(jìn)行對(duì)比。由圖4（a）可知，當(dāng)抗轉(zhuǎn)動(dòng)系數(shù)為0時(shí)，抗剪強(qiáng)度包線與前人室內(nèi)試驗(yàn)結(jié)果基本吻合且為直線，符合庫倫強(qiáng)度準(zhǔn)則的一般規(guī)律，驗(yàn)證了模擬試驗(yàn)的可行性與準(zhǔn)確性?？辜魪?qiáng)度包線的斜率隨抗轉(zhuǎn)動(dòng)系數(shù)的增大而增大，即轉(zhuǎn)動(dòng)阻抗的作用可有效增大峰值內(nèi)摩擦角。對(duì)于純圓顆粒，可根據(jù)實(shí)際顆粒的不規(guī)則程度及實(shí)際試樣的宏觀力學(xué)特性選取抗轉(zhuǎn)動(dòng)系數(shù)。

圖4 抗剪強(qiáng)度包線對(duì)比Fig. 4 Comparison of shear strength envelope

2.3 剪切帶

為了探究剪切過程中試樣內(nèi)不同位置的變形，在試樣中均勻設(shè)置了9個(gè)長方體區(qū)域，并將其中砂顆粒設(shè)置為黑色，監(jiān)測直剪過程中各組砂顆粒移動(dòng)情況。對(duì)比圖5（b）（c）可知，相比于初始階段，剪應(yīng)力達(dá)到峰值時(shí)，靠近兩側(cè)的砂顆粒出現(xiàn)相對(duì)位移，表明試樣內(nèi)砂顆粒的相對(duì)位移由左右兩側(cè)開始，變形由外向內(nèi)發(fā)展，最終變形集中于一條帶狀區(qū)域內(nèi)，即剪切帶（如圖5（d）），其厚度為79.2 mm，為平均粒徑的13～15倍，與潘遠(yuǎn)陽等[4]研究結(jié)果相吻合。

圖5 試樣變形演化Fig. 5 Evolution of specimen’s deformation

3 細(xì)觀結(jié)果分析

3.1 接觸力鏈

力鏈網(wǎng)絡(luò)是土骨架體系上粒間接觸力的宏觀形式，可體現(xiàn)顆粒體系受外荷載作用時(shí)受力響應(yīng)機(jī)制的敏感性。圖6為初始狀態(tài)、峰值狀態(tài)和最終狀態(tài)試樣內(nèi)力鏈的分布情況，其中力鏈直徑與法向接觸力的大小成正比。初始狀態(tài)，由于加載尚未進(jìn)行，伺服后試樣內(nèi)接觸力鏈分布均勻，各向異性不明顯。達(dá)到剪應(yīng)力峰值時(shí)，強(qiáng)力鏈向剪切帶區(qū)域靠攏，且沿直剪儀的斜對(duì)角分布。達(dá)到最終狀態(tài)時(shí)，強(qiáng)力鏈更加集中，且主要分布于剪切帶附近。

圖6 接觸力鏈網(wǎng)絡(luò)演化Fig. 6 Contact force chains in specimen

為進(jìn)一步定量分析試樣內(nèi)部接觸力演化規(guī)律，依據(jù)文獻(xiàn)[24]提出的強(qiáng)弱接觸力系統(tǒng)，分別對(duì)不同條件下試樣的法向和切向的強(qiáng)接觸力比例進(jìn)行量化統(tǒng)計(jì)。圖7和8分別為法向應(yīng)力為50 kPa時(shí)法向強(qiáng)接觸力和切向強(qiáng)接觸力的演化規(guī)律。由圖7（a）、8（a）可知，法向與切向強(qiáng)接觸力比例均隨著剪切進(jìn)行先降低，當(dāng)達(dá)到極值后逐漸趨于穩(wěn)定；抗轉(zhuǎn)動(dòng)系數(shù)越大，強(qiáng)接觸力比例越小。這是由于轉(zhuǎn)動(dòng)阻抗作用越大，顆粒間抗轉(zhuǎn)動(dòng)能力越強(qiáng)，顆粒間傳遞力矩越大，從而使強(qiáng)接觸力更大，而強(qiáng)接觸力主要抵抗外部荷載的作用，因此當(dāng)外部荷載相同時(shí)所需的強(qiáng)接觸力數(shù)越少。由圖7（b）、8（b）可明顯看出，隨著圍壓的增大，強(qiáng)接觸力比例極值呈微弱增大趨勢。

圖7 法向接觸力變化規(guī)律Fig. 7 Variation of normal contact force

圖8 切向接觸力變化規(guī)律Fig. 8 Variation of tangential contact force

3.2 配位數(shù)

配位數(shù)是指試樣中某顆粒與周圍顆粒的接觸數(shù)目，描述顆粒材料內(nèi)部顆粒的接觸程度。配位數(shù)又稱平均接觸數(shù)目，本文采用力學(xué)配位數(shù)[25]，其定義如下：

式中：Nc為試樣內(nèi)部總接觸數(shù)：Ns為試樣內(nèi)部總顆粒數(shù)；Ns0為某顆粒周圍無接觸的顆粒數(shù)；Ns1為某顆粒周圍接觸數(shù)目為1的顆粒數(shù)。

圖9 為200 kPa法向應(yīng)力條件下的平均力學(xué)配位數(shù)隨剪切位移的變化曲線，可看出配位數(shù)均隨剪切進(jìn)行而降低，后趨于穩(wěn)定。隨著抗轉(zhuǎn)動(dòng)系數(shù)的增大，配位數(shù)的穩(wěn)定值隨之減小，這是由于轉(zhuǎn)動(dòng)阻抗的作用使顆粒不易發(fā)生相對(duì)轉(zhuǎn)動(dòng)，進(jìn)而使顆粒間形成較大的接觸力，抗轉(zhuǎn)動(dòng)系數(shù)越大，傳遞能力越強(qiáng)，平均有效接觸就越少。

圖9 配位數(shù)隨剪切位移變化曲線Fig. 9 Variation curve of coordination number with shear displacement

3.3 各向異性

接觸組構(gòu)[26]是反映顆粒散料微細(xì)觀各向異性的常用指標(biāo)。剪切過程中對(duì)試樣內(nèi)顆粒間的法向接觸力和切向接觸力的數(shù)量及大小進(jìn)行統(tǒng)計(jì)分析，圖10和11展示了200 kPa法向應(yīng)力峰值強(qiáng)度時(shí)顆粒間法向接觸力和切向接觸力的各向異性分布玫瑰圖，圖中每10°為一個(gè)區(qū)間，每個(gè)區(qū)間內(nèi)接觸力的平均值與扇形面積相對(duì)應(yīng)。可明顯看出，法向接觸力分布玫瑰圖呈花生狀，切向接觸力分布玫瑰圖呈花瓣?duì)?，與有關(guān)文獻(xiàn)[27]研究結(jié)果相吻合。

圖10 法向接觸力分布Fig. 10 Distribution of normal contact force

圖11 切向接觸力分布Fig. 11 Distribution of tangential contact force

對(duì)比不同抗轉(zhuǎn)動(dòng)系數(shù)條件下的玫瑰圖可知，隨著抗轉(zhuǎn)動(dòng)系數(shù)的增大，“花生”變得更加細(xì)長，“花瓣”亦表現(xiàn)出類似變化，這表明轉(zhuǎn)動(dòng)阻抗的作用使法向接觸力和切向接觸力更加集中分布于部分區(qū)間，增大了這部分區(qū)間內(nèi)的接觸力平均值，進(jìn)而增大了各向異性程度。此外，隨著抗轉(zhuǎn)動(dòng)系數(shù)的增大，法向接觸力與切向接觸力的主方向均發(fā)生了順時(shí)針偏轉(zhuǎn)，其中法向接觸力偏轉(zhuǎn)相對(duì)明顯。

為進(jìn)一步定量分析顆粒材料的各向異性，參考組構(gòu)張量的概念[28]

其中：Fi j為組構(gòu)張量在不同方向上的分量表示組構(gòu)。在此基礎(chǔ)上，張鐸[29]根據(jù)廣義偏應(yīng)力概念采用組構(gòu)張量定義“組構(gòu)偏量”，即：

圖12 （a）為法向應(yīng)力是50 kPa時(shí)，偏組構(gòu)Fq隨剪切位移的變化曲線，可以看出隨著抗轉(zhuǎn)動(dòng)系數(shù)的增大，偏組構(gòu)峰值隨之增大，即試樣整體的各向異性隨之增大，但隨著剪切進(jìn)一步進(jìn)行，抗轉(zhuǎn)動(dòng)系數(shù)不為0的試樣的各向異性均逐漸降低，類似于應(yīng)變軟化趨勢，且抗轉(zhuǎn)動(dòng)系數(shù)越大，降低幅度越大。依據(jù)文獻(xiàn)[24]提出的強(qiáng)弱接觸力系統(tǒng)，將顆粒間接觸分為強(qiáng)弱接觸。圖12（b）給出了強(qiáng)接觸偏組構(gòu)隨剪切位移的變化曲線，其規(guī)律與整體偏組構(gòu)類似，但在整個(gè)剪切過程中，其大小均大于整體偏組構(gòu)，這表明轉(zhuǎn)動(dòng)阻抗對(duì)試樣整體各向異性的增強(qiáng)作用主要是增加強(qiáng)接觸的各向異性。分析主要原因：轉(zhuǎn)動(dòng)阻抗增大了強(qiáng)接觸傳遞力的能力和抵抗力鏈破壞的能力，使強(qiáng)接觸處顆?？蓚鬟f更大的力，同時(shí)減小了接觸數(shù)量，即使接觸力的傳遞更加集中，各向異性程度更大。圖12（c）為強(qiáng)接觸偏組構(gòu)峰值隨法向應(yīng)力的變化曲線，可見隨著法向應(yīng)力的增大，不同抗轉(zhuǎn)動(dòng)系數(shù)試樣的各向異性均呈下降趨勢。

Oda[30]提出采用組構(gòu)張量描述顆粒材料內(nèi)的接觸特性，通過分布密度函數(shù)E(Ω)定義組構(gòu)張量，并采用組構(gòu)張量的第二不變量定義其各向異性系數(shù)ac。

式中：φi j為 組構(gòu)張量；ni、nj為接觸法向在笛卡爾坐標(biāo)系中的分量為描述全局接觸法向分布密度的參量；ac為的第二不變量。

與組構(gòu)張量各向異性系數(shù)的描述類似，法向接觸力張量也可定義其各向異性系數(shù)an。圖13和14為各向異性系數(shù)隨軸向應(yīng)變的變化曲線。隨著剪切進(jìn)行，試樣發(fā)生膨脹變形，ac和an的演化規(guī)律十分相似，均表現(xiàn)為先急劇增加，相比于無轉(zhuǎn)動(dòng)阻抗作用的試樣，其他試樣達(dá)到峰值后隨軸向應(yīng)變顯示出軟化趨勢，且抗轉(zhuǎn)動(dòng)摩擦系數(shù)越大，軟化趨勢越明顯，與剪應(yīng)力-剪切位移曲線相對(duì)應(yīng)；對(duì)比圖13和14可知，整個(gè)剪切過程中an大 于ac，這是由法向接觸力大小的非均勻分布導(dǎo)致的，同時(shí)表明法向接觸力各向異性對(duì)試樣抗剪強(qiáng)度影響較大；隨著抗轉(zhuǎn)動(dòng)系數(shù)的增大，ac與an均隨之增大，即在砂試樣中產(chǎn)生更顯著的剪切誘導(dǎo)各向異性，從而提高了砂的穩(wěn)定性及抗剪強(qiáng)度。

圖13 接觸法向各向異性系數(shù)αc隨剪切位移變化曲線Fig. 13 Curve of contact normal anisotropy coefficient αc with axial strain

圖14 法向接觸力各向異性系數(shù)αn隨剪切位移變化曲線Fig. 14 Curve of normal contact force anisotropy coefficient αn with axial strain

4 結(jié) 語

（1）數(shù)值試驗(yàn)宏觀應(yīng)力-應(yīng)變特性、強(qiáng)度指標(biāo)及剪切帶變形規(guī)律與前人研究結(jié)果相類似，表明砂土直剪數(shù)值模擬試驗(yàn)存在良好的可行性；相比于無轉(zhuǎn)動(dòng)阻抗作用的試樣，轉(zhuǎn)動(dòng)阻抗的作用使砂表現(xiàn)出明顯的應(yīng)力軟化及剪脹特性，且抗轉(zhuǎn)動(dòng)系數(shù)越大，剪應(yīng)力峰值越大。

（2）隨著剪切進(jìn)行，力鏈分布區(qū)域化，達(dá)到峰值后集中于剪切帶區(qū)域，且發(fā)生了主方向的偏轉(zhuǎn)。由于轉(zhuǎn)動(dòng)阻抗增大顆粒間傳遞接觸力的能力，外部荷載相同條件下，抗轉(zhuǎn)動(dòng)系數(shù)越大，所需接觸數(shù)越少，即表現(xiàn)為平均力學(xué)配位數(shù)隨抗轉(zhuǎn)動(dòng)系數(shù)的增大而減小，且強(qiáng)弱接觸力比例隨之更加懸殊。

（3）玫瑰圖表明轉(zhuǎn)動(dòng)阻抗的作用增強(qiáng)了法向接觸力和切向接觸力的各向異性，偏組構(gòu)的變化規(guī)律說明轉(zhuǎn)動(dòng)阻抗主要增大了強(qiáng)接觸力的各向異性，而張量的第二不變量演化規(guī)律則表明轉(zhuǎn)動(dòng)阻抗對(duì)法向接觸力各向異性的增強(qiáng)作用對(duì)提高宏觀抗剪強(qiáng)度具有顯著作用。