偏心率對(duì)顆粒介質(zhì)次生各向異性的影響

周小文 許衍彬 趙仕威 陳昊 張昌輝

（華南理工大學(xué)亞熱帶建筑科學(xué)國(guó)家重點(diǎn)實(shí)驗(yàn)室∥華南巖土研究院，廣東 廣州 510640）

顆粒介質(zhì)于巖土水利工程中應(yīng)用較為廣泛，常見于尾礦砂蓄積、大砂袋吹填、路基碎石鋪設(shè)等［1-3］。然而，顆粒介質(zhì)在細(xì)觀上呈現(xiàn)出較宏觀尺度更為明顯的不連續(xù)性與組構(gòu)各向異性，而其宏觀層面物理力學(xué)特性的各向異性與其細(xì)觀組構(gòu)各向異性密不可分［4-5］，因此對(duì)各向異性的細(xì)觀尺度研究具有重要意義。Casagrande［6］依據(jù)各向異性成因?qū)⑵鋭澐譃樵飨虍愋院痛紊飨虍愋裕飨虍愋栽从陬w粒集合體在天然沉積過程中因顆粒級(jí)配及顆粒形狀導(dǎo)致的顆粒排列及接觸的空間差異，而次生各向異性則源于顆粒在外部荷載擾動(dòng)下的原生組構(gòu)重新排列。Oda等［7］進(jìn)一步地基于顆粒集合體的接觸法向與接觸力等細(xì)觀組構(gòu)張量給出了各向異性的定量表征，為從細(xì)觀尺度解析顆粒介質(zhì)諸多典型宏觀力學(xué)現(xiàn)象（例如，擠壓膨脹效應(yīng)［8］、糧倉(cāng)效應(yīng)［9］、拱效應(yīng)［10］、剪切帶［11-12］及液化［13-14］等）的機(jī)理提供了可能；Rothenburg等［15］提出的應(yīng)力-力-組構(gòu)（SFF）關(guān)系即是從細(xì)觀尺度出發(fā)研究顆粒介質(zhì)宏觀力學(xué)特性的典型代表，相關(guān)研究表明了顆粒介質(zhì)宏觀抗剪強(qiáng)度與細(xì)觀組構(gòu)各向異性之間密切的跨尺度聯(lián)系。因此，為進(jìn)一步揭示顆粒介質(zhì)宏觀力學(xué)現(xiàn)象的細(xì)觀機(jī)理，基于組構(gòu)各向異性的SFF關(guān)系進(jìn)行跨尺度研究即為一條重要途徑。

對(duì)顆粒形狀特征進(jìn)行準(zhǔn)確描述是進(jìn)行細(xì)觀尺度研究的基礎(chǔ)。然而由于建模方便及計(jì)算效率限制，既有顆粒材料組構(gòu)各向異性的研究多基于二維模擬或形狀規(guī)則顆粒（如理想的圓盤或球顆粒［16-19］）。但是，理想球顆粒與自然界中真實(shí)復(fù)雜顆粒形狀相差甚遠(yuǎn)，前者不僅顆粒間接觸形式單一，且由于法向接觸力指向顆粒質(zhì)心而無(wú)法考慮真實(shí)顆粒間的轉(zhuǎn)動(dòng)阻抗與咬合互鎖等作用。導(dǎo)致外部荷載驅(qū)動(dòng)下顆粒旋轉(zhuǎn)自由度增大，進(jìn)而顆粒定向排列現(xiàn)象減弱且組構(gòu)空間分布趨向均勻，顆粒分布隨機(jī)性增強(qiáng)，所以細(xì)觀組構(gòu)各向異性及其影響下的力學(xué)特性無(wú)法被真實(shí)地揭示［20］。基于此，為考慮顆粒形狀特征的影響并保證計(jì)算效率，有學(xué)者［21-22］嘗試在球顆粒中引入轉(zhuǎn)動(dòng)阻抗模型來(lái)間接考慮顆粒形狀特征的影響，但此類折中表征方式在模擬組構(gòu)各向異性演化方面尚有不足［23-24］。因此，直接采用趨真非球顆粒來(lái)反映顆粒形狀特征的影響并結(jié)合更為豐富的非球顆粒組構(gòu)各向異性表征方式進(jìn)行研究是一條更具科學(xué)意義的途徑。非球顆粒往往具有顆粒長(zhǎng)軸并存在長(zhǎng)軸排列優(yōu)先性，因此有學(xué)者［25］通過顆粒長(zhǎng)軸分布方向性來(lái)額外表征組構(gòu)各向異性。此外，由于非球顆粒間接觸法向與支向量不相重合，且二者夾角隨顆粒不規(guī)則性增大而增大［26］（此種不重合也產(chǎn)生了額外力矩［27］），因此非球顆粒集合體中支向量對(duì)幾何各向異性奉獻(xiàn)較大，Ouadfel等［28］針對(duì)橢球顆粒額外引進(jìn)支向量表征組構(gòu)各向異性。

鑒于上文所述，顆粒形狀特征對(duì)顆粒介質(zhì)組構(gòu)各向異性的演化起著重要作用，趨真非球顆粒組構(gòu)各向異性的研究也逐漸增多。比如，由細(xì)觀組構(gòu)各向異性引起的典型宏觀力學(xué)現(xiàn)象受顆粒形狀的影響逐漸被關(guān)注，諸如拱效應(yīng)［29］、剪切帶［30］、液化［31-32］等；應(yīng)力-力-組構(gòu)（SFF）關(guān)系理論的適用性也先后被多種趨真非球顆粒所驗(yàn)證，如橢球［28］、超級(jí)球［33］、超級(jí)橢球［24］、多面體［31］及球簇［34］等。同時(shí)較多學(xué)者基于上述組構(gòu)各向異性的表征方式，致力于顆粒形狀特征對(duì)次生細(xì)觀各向異性的影響研究。譬如，有學(xué)者采用準(zhǔn)球多面體［35-36］、準(zhǔn)球球簇［37］以及超級(jí)球［33］來(lái)探究棱角度、球度對(duì)次生細(xì)觀各向異性演化的影響。亦有學(xué)者探究長(zhǎng)細(xì)比對(duì)次生細(xì)觀各向異性演化的影響，此方面不乏橢球［26，38］、超級(jí)橢球［24，39］、延長(zhǎng)多面體［40］、RCR顆粒［41］、球簇［42-43］等非球顆粒。但作者注意到，以上研究將長(zhǎng)細(xì)比與棱角度等指標(biāo)糅合起來(lái)，以致于難以精準(zhǔn)探究單一形態(tài)指標(biāo)對(duì)組構(gòu)各向異性的影響。并且，在已有的相關(guān)文獻(xiàn)中，未能找到有學(xué)者對(duì)顆粒偏心不對(duì)稱進(jìn)行量化與研究。而實(shí)際上偏心率越大，法向接觸力偏離質(zhì)心愈發(fā)顯著，產(chǎn)生了額外的繞質(zhì)心轉(zhuǎn)動(dòng)阻抗。繼而，顆粒之間互鎖咬合效應(yīng)隨之增大，顆粒旋轉(zhuǎn)自由度也隨之減小。因此，組構(gòu)各向異性和宏觀抗剪強(qiáng)度受到影響。對(duì)此，本研究提出偏心率量化顆粒主軸不對(duì)稱的程度，一定程度上完善了顆粒形態(tài)特征的量化指標(biāo)體系。并繼而嘗試基于非對(duì)稱擴(kuò)展超級(jí)橢球進(jìn)行了不同偏心率的顆粒介質(zhì)趨真建模，以期探究偏心率對(duì)顆粒介質(zhì)次生細(xì)觀各向異性演化的影響，為進(jìn)一步揭示眾多力學(xué)現(xiàn)象的細(xì)觀機(jī)理并建立考慮顆粒形狀特征的顆粒材料各向異性細(xì)觀本構(gòu)提供參考。

1 離散元模型建立

1.1 趨真顆粒形狀模型

超橢球在局部笛卡爾坐標(biāo)系下方程［44］為

式中：rx、ry、rz分別代表x、y、z軸的半主軸長(zhǎng)；ε1、ε2衡量顆粒表面棱角度。

超橢球雖可刻畫自然界中80%的顆粒形態(tài)特征，但仍不能定量刻畫現(xiàn)實(shí)顆粒的偏心現(xiàn)象。為對(duì)顆粒偏心率的變化進(jìn)行定量描述，筆者所在團(tuán)隊(duì)在超級(jí)橢球模型的基礎(chǔ)上開發(fā)了擴(kuò)展超級(jí)橢球［45］趨真顆粒模型，如圖1所示。其由8個(gè)不同形狀超級(jí)橢球的八分體在保證表面連續(xù)、光滑的前提下組成，在表征顆粒形狀特征上相比其他非球顆粒具有更高廣度，且其接觸檢測(cè)算法的魯棒性和效率性也已通過堆積試驗(yàn)（動(dòng)態(tài)）和三軸試驗(yàn)（準(zhǔn)靜態(tài)）得以驗(yàn)證［45］。

圖1 擴(kuò)展超級(jí)橢球的構(gòu)建原理Fig.1 building principle of poly-superellipsoid

擴(kuò)展超級(jí)橢球表面方程及相關(guān)參數(shù)為

式中：rix、riy、riz代表第i個(gè)八分體在x、y、z主軸的半主軸長(zhǎng)；八分體的序號(hào)i如圖1所示；r+x、r-x分別代表擴(kuò)展超級(jí)橢球在x主軸正負(fù)方向上的半主軸長(zhǎng)，其他同理；εi1、εi2衡量第i個(gè)八分體的表面棱角度。

偏心率具體定義如下

式中：i遍歷x、y、z軸；ζi為i軸的偏心率；r+i、r-i分別為i軸在正負(fù)方向上的長(zhǎng)度。

特別地，本研究獨(dú)立于顆粒形態(tài)棱角度與長(zhǎng)細(xì)比展開研究，即令ε1=ε2=1.0以及l(fā)x=ly=lz。同時(shí)將x、y、z軸三個(gè)方向的偏心率統(tǒng)一起來(lái)，即令ζ=ζx=ζy=ζz并在保證體積一致的情況下將統(tǒng)一后的偏心率ζ設(shè)置為0、0.2、0.4、0.6、0.8，以S1-S5分別代之。

數(shù)值模擬基于筆者所在團(tuán)隊(duì)開發(fā)的非球顆粒開源離散元程序SudoDEM［25，46］。顆粒接觸模型采用線性彈簧模型及庫(kù)倫滑移模型。擴(kuò)展超橢球接觸檢測(cè)及計(jì)算的詳細(xì)信息參見https：//sudodem.github.io。

1.2 數(shù)值模擬流程

1.2.1 試樣制備

參考Zhou等［47］采用的半徑擴(kuò)大法生成如圖2所示的試樣，其中顆粒等效直徑即等體積球直徑。試樣尺寸為11.36 mm×11.36 mm×11.36 mm，其與最大顆粒直徑之比為13.36，滿足Jamiolkowski等［48］為減小尺寸效應(yīng)和應(yīng)變局部化而提出的標(biāo)準(zhǔn)。顆粒數(shù)目為6000個(gè)（預(yù)試驗(yàn)證明此顆粒數(shù)目已達(dá)基本要求）。模擬不考慮重力以避免原生各向異性影響。

圖2 離散元試樣生成Fig.2 Generations of DEM specimens

模擬細(xì)觀參數(shù)如表1所示。保證模擬穩(wěn)定的前提下采取密度放大法將顆粒材料密度放大106倍［18，49］后可將時(shí)步提升至5×10-5s。接觸剛度采用100MPa×r［18，50］，其中r為顆粒等效半徑均值。進(jìn)行準(zhǔn)靜態(tài)模擬時(shí)，系統(tǒng)慣性系數(shù)Iinertia不超過0.001［51-52］，即

表1 DEM數(shù)值模擬細(xì)觀參數(shù)Table 1 Meso-parameters used in the DEM simulations

式中：ε為軸向加載應(yīng)變率，控制為0.01/s；d為顆粒等效直徑的均值（0.6 mm）；ρ為顆粒材料密度（2 650×106kg/m3）；σ3為圍壓（100 kPa）。經(jīng)計(jì)算Iinertia為0.000 976，滿足準(zhǔn)靜態(tài)模擬標(biāo)準(zhǔn)要求。

1.2.2 模擬方案

數(shù)值模擬分固結(jié)與剪切兩個(gè)階段，如圖3所示。固結(jié)過程生成初始各向同性試樣，而剪切過程研究次生各向異性的演化。固結(jié)開始前，固定邊界墻位置不動(dòng)的同時(shí)周期性施加人工阻尼使系統(tǒng)快速穩(wěn)定，然后按應(yīng)力控制伺服機(jī)制，讓顆粒在100kPa圍壓下固結(jié)，此過程調(diào)控摩擦系數(shù)［18］來(lái)保證不同顆粒形狀試樣的初始孔隙比一致為0.626［24，33，53］。固結(jié)完成后開展真三軸剪切，即頂?shù)變擅鎵σ?.01/s的恒定應(yīng)變加載速率相向移動(dòng)，側(cè)壁四面墻按應(yīng)力控制伺服機(jī)制獨(dú)立移動(dòng)，維持圍壓100 kPa。為將顆粒間疊合量控制于平均粒徑（d50）的0.7%內(nèi)，最終剪切應(yīng)變?yōu)?0%。

圖3 不同偏心率試樣的真三軸剪切Fig.3 True triaxial shear of varying eccentricity samples

2 結(jié)果分析

2.1 宏觀力學(xué)響應(yīng)

參考Christoffersen等［54］基于細(xì)觀層面顆粒間接觸力與支向量來(lái)計(jì)算應(yīng)力張量σij。

式中：V為試樣體積（包括孔隙體積）；Nc為總接觸數(shù)；f和l為顆粒間接觸力（包含法向與切向）和支向量；i、j∈{ }1，2，3，代表笛卡爾坐標(biāo)系下x、y、z方向。由于邊界光滑，因此試樣邊界剪應(yīng)力均為零，即σij為對(duì)角矩陣。

平均主應(yīng)力p和偏應(yīng)力張量σ′ij采用下式計(jì)算：

式中，tr(σij)為σij的跡，δij為科羅內(nèi)多張量。

采用廣義剪應(yīng)力q量化偏應(yīng)力張量［18］：

式中，σ′ijσ′ij根據(jù)愛因斯坦求和約定計(jì)算。

由于本研究采用剛性邊界墻進(jìn)行真三軸剪切，可從宏觀層面根據(jù)邊界墻的位移來(lái)計(jì)算對(duì)數(shù)軸向應(yīng)變［18］。

式中，H0、H分別為剪切前和剪切過程的試樣高度。

不同偏心率試樣的偏應(yīng)力比q/p隨軸向應(yīng)變的演化曲線如圖4所示?？捎^察到偏心率的小偏差導(dǎo)致了顯著不同的宏觀抗剪強(qiáng)度演化形態(tài)。顯然，試樣臨界強(qiáng)度隨著偏心率增大而單調(diào)增大，偏心顆粒有著顯著增強(qiáng)的臨界抗剪強(qiáng)度。這是因?yàn)椋弘S偏心率增大，顆粒間法向接觸力偏離質(zhì)心愈明顯，產(chǎn)生了更大的繞質(zhì)心額外力矩［27］，增大了轉(zhuǎn)動(dòng)阻抗，進(jìn)而增強(qiáng)了顆粒互鎖和咬合牢固性，使得試樣可承擔(dān)更大的剪應(yīng)力。Kozicki等［55］、Zhao等［36］也有過類似報(bào)道。另外，雖試樣的初始孔隙比一致，但試樣的相對(duì)密實(shí)度隨顆粒偏心率改變而不同。具體地，偏心率較大的試樣相對(duì)密實(shí)度較小，偏心率較小的試樣相對(duì)密實(shí)度較大。因此，隨偏心率增大試樣逐漸由應(yīng)變軟化轉(zhuǎn)向應(yīng)變硬化。

圖4 偏應(yīng)力比-軸向應(yīng)變曲線Fig.4 Deviatoric stress ratio-axial strain

2.2 細(xì)觀組構(gòu)演化

2.2.1 配位數(shù)

配位數(shù)為表征顆粒材料內(nèi)部組構(gòu)的細(xì)觀標(biāo)量，某顆粒的配位數(shù)即與此顆粒相接觸的顆粒數(shù)。鑒于配位數(shù)為0或者1的顆粒對(duì)組構(gòu)穩(wěn)定無(wú)貢獻(xiàn)，本文通過力學(xué)平均配位數(shù)M描述內(nèi)部組構(gòu)穩(wěn)定性［57］，即

式中：N1表示配位數(shù)為1的顆粒數(shù)目；N0表示配位數(shù)為0的顆粒數(shù)目；Nc表示總接觸數(shù)目（包括顆粒—墻接觸）；N表示總的顆粒數(shù)目。

力學(xué)平均配位數(shù)M隨軸向應(yīng)變的演化曲線如圖5所示?？梢姼髟嚇恿W(xué)平均配位數(shù)隨應(yīng)變發(fā)展呈指數(shù)衰減并最終達(dá)臨界值，即力學(xué)穩(wěn)定jamming態(tài)。根據(jù)Rothenburg等［58］和Zhao等［59］的相關(guān)研究，配位數(shù)的減小主要為密實(shí)的試樣在剪切過程中體積剪脹所帶來(lái)的塑性變形所致。另外，隨偏心率的增大，配位數(shù)下降幅度減小。所以，臨界配位數(shù)隨偏心率單調(diào)增加，對(duì)應(yīng)前文所述較高的臨界抗剪強(qiáng)度?？赡茉?yàn)椋浩穆视螅w粒間接觸面積更大，內(nèi)嵌愈加縝密，因此試樣配位數(shù)愈大，內(nèi)部組構(gòu)愈穩(wěn)定，為其擁有更高強(qiáng)度作出細(xì)觀解釋。Santamarina等［60］也指出非球顆粒需要更大的配位數(shù)來(lái)維持穩(wěn)定。

圖5 力學(xué)平均配位數(shù)-軸向應(yīng)變曲線Fig.5 Mechanical average coordination number-axial strain

2.2.2 接觸力傳遞

接觸間承擔(dān)著或大或小的接觸力，當(dāng)給定接觸的法向接觸力小于試樣的平均法向接觸力，此接觸即定義為弱接觸，反之為強(qiáng)接觸。弱接觸比例ξw即為弱接觸與總接觸的數(shù)量之比。

顆粒材料在變形過程中往往伴隨著顆粒相互滑動(dòng)。Alonso等［61］認(rèn)為滑動(dòng)接觸的各向異性為顆粒材料塑性變形的主要原因。滑動(dòng)接觸比例ξs即為滑動(dòng)接觸與總接觸的數(shù)量之比。

摩擦發(fā)揮系數(shù)I用于衡量接觸間摩擦發(fā)揮程度，其定義如式（18）［62］所示。由于切向接觸力的計(jì)算采用了庫(kù)倫滑移模型，摩擦發(fā)揮系數(shù)I不大于1。當(dāng)I為1時(shí)表明接觸發(fā)生相互滑動(dòng)，顆粒集合體重新排列。

式中：fn和ft分別為法、切向接觸力；μ為摩擦系數(shù)。

弱接觸比例ξw、滑動(dòng)接觸比例ξs隨軸向應(yīng)變演化以及臨界狀態(tài)（εa=50%）的摩擦發(fā)揮系數(shù)概率分布（PDF）曲線分別如圖6、圖7所示。

由圖6可見，弱接觸比例超過60%并始終大于強(qiáng)接觸比例，說明弱接觸網(wǎng)格始終比強(qiáng)接觸網(wǎng)格承擔(dān)著更大的作用。Zhao等［33］也曾有過相關(guān)報(bào)道；弱接觸比例及滑動(dòng)接觸比例均隨偏心率增大而增大，但后者較前者對(duì)偏心率更為敏感。

圖6 弱接觸、滑動(dòng)接觸比例-軸向應(yīng)變曲線Fig.6 Weak and slide contacts proportion-axial strain

由圖7可見，增大的偏心率減小了低摩擦發(fā)揮接觸的比例并增大了高摩擦發(fā)揮接觸的比例，當(dāng)I小于0.75時(shí)隨偏心率的增大I的PDF逐漸減小，而當(dāng)I大于0.80時(shí)隨偏心率增大I的PDF逐漸增大，Nie等［34］也曾得到類似結(jié)論。對(duì)此，Estrada等［63］曾作相關(guān)解釋：顆粒的偏心不規(guī)則性加劇了接觸互鎖現(xiàn)象，提高了顆?？剐D(zhuǎn)能力，進(jìn)而導(dǎo)致顆粒間接觸需通過相互滑動(dòng)來(lái)抵抗外荷載。最終，滑動(dòng)接觸比例和高摩擦發(fā)揮接觸的比例均增大。而從另一角度，根據(jù)Zhao等人［24，33］以及Radjai等人［64］的相關(guān)研究，絕大部分滑動(dòng)接觸發(fā)生在弱接觸。換言之，較小法向接觸力致最大靜摩擦力較易被打破。所以弱接觸比例的提升間接導(dǎo)致了滑動(dòng)接觸比例和高摩擦發(fā)揮比例的增大。

圖7 臨界狀態(tài)摩擦發(fā)揮系數(shù)PDF分布Fig.7 PDF of friction mobilized index at critical states

通過力鏈圖可將顆粒間法向接觸力可視化，進(jìn)而直觀展現(xiàn)各向異性，如圖8所示。法向接觸力方向由顆粒質(zhì)心連線（即支向量）代表，而大小則由連線粗細(xì)與顏色代表。

圖8 初始及臨界狀態(tài)接觸力鏈圖Fig.8 Contact force chains at initial and critical states

各試樣在剪切前的力鏈分布相近，因此僅針對(duì)S5試樣進(jìn)行展示，如圖8（a）所示。隨剪切的進(jìn)行，稀疏的強(qiáng)力鏈似柱子般承擔(dān)著外部豎向荷載并形成有效的傳力路徑，而密集且始終占據(jù)主導(dǎo)地位的弱力鏈則維持著強(qiáng)力鏈的穩(wěn)定。值得注意的是，弱力鏈幾乎呈各向同性分布，而強(qiáng)力鏈卻僅沿豎直加載方向分布，呈現(xiàn)出顯著的各向異性，表明強(qiáng)力鏈主導(dǎo)了法向接觸力的組構(gòu)各向異性，F(xiàn)oroutan等［65］有著類似結(jié)論。隨偏心率增大，強(qiáng)力鏈數(shù)量上雖無(wú)明顯變化，但其管徑膨脹，加劇了豎直向與水平向法向接觸力的各向異性。雖高偏心率試樣對(duì)應(yīng)更高的強(qiáng)度（見圖4），但其剪應(yīng)力的承擔(dān)卻更加不均勻。此外如前文所述，隨偏心率增大，接觸數(shù)增大的同時(shí)弱接觸比例也增大，因此弱力鏈更密集，保證了強(qiáng)力鏈的穩(wěn)定，進(jìn)而使其可承受更大的外部荷載。

2.2.3 組構(gòu)各向異性

諸多學(xué)者［18，24，28，47，65］通過接觸法向、接觸力和支向量的組構(gòu)張量對(duì)組構(gòu)各向異性進(jìn)行定量描述。

接觸法向c指接觸的外法線方向，其各向異性可用組構(gòu)張量φij定量描述［7］，

式中：φij為3×3二階對(duì)稱張量；Nc為總接觸數(shù)；n為接觸法向單位向量，ni（i∈{1，2，3}）代表n在笛卡爾坐標(biāo)系下x、y、z3個(gè)方向的分量；上式連續(xù)空間積分向離散求和形式的轉(zhuǎn)化詳見文獻(xiàn)［66］。

法向接觸力fn、切向接觸力ft的組構(gòu)張量為

式中：acij為接觸法向偏組構(gòu)張量，詳見后文；根據(jù)愛因斯坦求和約定計(jì)算；|fn|為法向接觸力大小；|ft|為切向接觸力大??；t為接觸切向單位向量，t（ii∈{1，2，3}）代表t在笛卡爾坐標(biāo)系下x、y、z3個(gè)方向的分量。

支向量指顆粒間質(zhì)心連線。由于顆粒形狀的影響，非球顆粒的支向量與接觸法向不重合。對(duì)此，以接觸法向?yàn)閰⒖枷祵⒅蛄糠纸鉃榉ㄏ蚍至亢颓邢蚍至浚ㄏ蛑蛄縟n與切向支向量dt的組構(gòu)張量分別為

式中：|dn|為支向量在接觸法向方向的分量；|dt|為支向量在接觸法向的垂直分量。

基于上述認(rèn)識(shí)，5類組構(gòu)的偏組構(gòu)張量a*ij分別為

求得偏組構(gòu)張量后便可通過第二不變量A*來(lái)量化各向異性［18］。

式中：a（*ij*代表c、fn、ft、dn、d）t為偏組構(gòu)張量；σ′ij為偏應(yīng)力張量；當(dāng)二者共軸時(shí)A*為正，反之為負(fù)；5類組構(gòu)各向異性的第二不變量A*隨應(yīng)變?chǔ)臿的演化曲線如圖9所示，隨偏心率ζ的演化曲線（在εa為50%的狀態(tài)下）如圖10所示。

圖9 組構(gòu)各向異性-軸向應(yīng)變曲線Fig.9 Fabric anisotropy-axial strain

圖10 組構(gòu)各向異性-偏心率曲線Fig.10 Fabric anisotropy-eccentricity

接觸法向在單位球面下的概率密度函數(shù)E(Θ)為

式中：a1∈［0，2π］，a2∈［0，π］。

5類組構(gòu)在單位球面下的概率密度函數(shù)可以用傅里葉級(jí)數(shù)展開?？紤]到方向量的對(duì)稱性，奇數(shù)階張量對(duì)級(jí)數(shù)解無(wú)貢獻(xiàn)［67］，可采用二階傅里葉級(jí)數(shù)展開［15，65］。

式中：E(Θ)、Fn(Θ)、Fit(Θ)、Dn(Θ)、Dti(Θ)分別為接觸法向c、法向接觸力fn、切向接觸力ft、法向支向量dn、切向支向量dt的概率密度函數(shù)。

為將組構(gòu)各向異性可視化，可基于式（31）-（35）對(duì)5類組構(gòu)的拓?fù)浞植紙D進(jìn)行傅里葉擬合。

由于剪切前試樣處于各向同性應(yīng)力狀態(tài)，5類組構(gòu)在剪切前也均處于初始各向同性狀態(tài)。圖8（a）中各向同性的弱力鏈以及圖9中5類組構(gòu)各向異性第二不變量的初始值均為零便作了證明，因此鑒于篇幅有限不對(duì)初始狀態(tài)的組構(gòu)分布圖進(jìn)行展示。

圖9中5類組構(gòu)各向異性第二不變量A*在剪切過程中的演化各有不同，但均從各向同性狀態(tài)隨剪切發(fā)展出次生各向異性，并在不同的應(yīng)變水平達(dá)到不同的各向異性臨界狀態(tài)。其中Ac演化形態(tài)最接近圖4的強(qiáng)度演化形態(tài)。

在試樣的臨界狀態(tài)，5類組構(gòu)的拓?fù)浞植几饔胁煌?。從圖11可觀察到，接觸法向c和法向接觸力fn的拓?fù)浞植汲省昂J”狀（其中，法向接觸力fn的“葫蘆”更高且更細(xì)）。這是因?yàn)椋佑|法向c和法向接觸力fn均在豎直方向（大主應(yīng)力方向）上一定程度地增長(zhǎng)。然而，水平向則不然，甚至出現(xiàn)了收縮現(xiàn)象。這導(dǎo)致了c與fn次生各向異性的產(chǎn)生。另一方面，切向接觸力ft與切向支向量dt的拓?fù)浞植汲省半p葉草”狀。這是因?yàn)?，ft與dt均在45°+n×90°(n=0，1，2，3)的方向上出現(xiàn)較大幅度的增長(zhǎng)。但同時(shí)，二者在0°+n×90°(n=0，1，2，3)的方向上卻基本保持為零。特別地，ft與dt“雙葉草”狀的拓?fù)浞植寂c現(xiàn)實(shí)中砂土試樣常見的剪切帶形態(tài)較為相像。這表明，切向接觸力ft與切向支向量dt的各向異性與砂土的剪切帶有著密不可分的聯(lián)系。值得注意的是，法向支向量dn的拓?fù)浞植紴閳A形。這表明，法向支向量dn為各向同性分布或其各向異性可忽略不計(jì)。這也是圖9（d）中Adn在整個(gè)剪切過程中始終很小且其演化曲線崎嶇不平的原因所在。

在試樣的臨界狀態(tài)，偏心率對(duì)組構(gòu)各向異性的影響也各不相同。偏心率對(duì)接觸法向c的各向異性影響很小。圖9（a）、圖11（a）和Zhao等［33］均證明了此結(jié)論。另一方面，在圖9（b）當(dāng)中，Afn隨偏心率的增大而小幅度地增大；在圖11（b）當(dāng)中，隨偏心率的增大，法向接觸力的拓?fù)浞植甲兊酶印案呤荨?；在圖8中，隨偏心率增大，強(qiáng)力鏈的管徑出現(xiàn)了膨脹。這三個(gè)現(xiàn)象均表明，偏心率增大了法向接觸力fn的各向異性。切向接觸力ft的各向異性隨偏心率增大而大幅度增長(zhǎng)。此觀點(diǎn)可從圖9（c）中隨偏心率增大而大幅增大的Aft、圖11（c）中隨偏心率增大而整體變大的“雙葉草”得出。最后，圖9（e）和圖11（e）均可表明切向支向量dt各向異性則隨偏心率增大而不等幅度增大。這是因?yàn)椋山佑|法向與支向量的不相重合產(chǎn)生的切向支向量隨顆粒形態(tài)特征偏離球體而逐漸增大［26］。

圖11 臨界狀態(tài)組構(gòu)分布圖Fig.11 Distribution of fabric at critical state

從圖10可見，切向接觸力ft、切向支向量dt的各向異性與偏心率的關(guān)系線的斜率較高，表明偏心率對(duì)二者的影響較大?？梢哉J(rèn)為，在偏心顆粒試樣中，切向接觸力ft、切向支向量dt對(duì)幾何各向異性奉獻(xiàn)較大，不應(yīng)忽視，這與其他學(xué)者的結(jié)論有所偏差［33］。

5類組構(gòu)各向異性對(duì)剪應(yīng)力比（宏觀抗剪強(qiáng)度）的貢獻(xiàn)各不相同。圖12展示了5類組構(gòu)各向異性的貢獻(xiàn)比重隨軸向應(yīng)變的演化。從中可見，Afn權(quán)重始終最大。這表明，法向接觸力fn的各向異性對(duì)抗剪強(qiáng)度的貢獻(xiàn)最大。同時(shí)，在圖9（b）中，Afn的演化形態(tài)與宏觀強(qiáng)度的演化形態(tài)始終最為接近，此現(xiàn)象也對(duì)上述觀點(diǎn)作出了佐證。文獻(xiàn)［18，24，31］也曾報(bào)道過類似的觀點(diǎn)。法向接觸力fn及接觸法向c各向異性總占比約80%，而其余三者權(quán)重均較小。這表明接觸法向c與法向接觸力fn各向異性對(duì)強(qiáng)度的貢獻(xiàn)較大。

圖12 組構(gòu)各向異性比重-軸向應(yīng)變曲線Fig.12 Weights of fabric anisotropy-axial strain

隨軸向應(yīng)變的發(fā)展，Ac的貢獻(xiàn)比重呈現(xiàn)出先增加后減小的趨勢(shì)，Afn的貢獻(xiàn)比重則呈現(xiàn)出先減小后增大的趨勢(shì)。它們分別在約15%的軸向應(yīng)變處達(dá)到峰值與谷值。這說明，剪切過程中Ac的權(quán)重和Afn相互轉(zhuǎn)移。這意味著，剪切過程中幾何各向異性與力學(xué)各向異性相互轉(zhuǎn)化。其原因正如Zhao等［33］所述：剪切過程中因剪脹而不斷減小的平均配位數(shù)使得接觸分散性降低。

在偏心率的影響方面，從圖12得知，隨偏心率增大，Ac、Afn的相對(duì)權(quán)重減小，Aft的相對(duì)權(quán)重增大。這表明，偏心率的增大會(huì)讓法向接觸力的各向異性向切向接觸力的各向異性轉(zhuǎn)化。此觀點(diǎn)解釋如下：一方面，如前文所述，弱力鏈呈各向同性，強(qiáng)力鏈才為法向接觸力各向異性的主導(dǎo)。Guo等［18］也曾指出，強(qiáng)接觸網(wǎng)格中法向接觸力各向異性比弱接觸網(wǎng)格大得多。因此，隨偏心率增大，弱力鏈比例的提高、強(qiáng)力鏈比例的減小導(dǎo)致了法向接觸力各向異性比重的減?。涣硪环矫?，Sufian等［68］認(rèn)為：滑動(dòng)接觸網(wǎng)格中切向接觸力的各向異性相比法向接觸力各向異性更大。也即，較高偏心率帶來(lái)的滑動(dòng)接觸比例的增大會(huì)導(dǎo)致切向接觸力各向異性比重增大。綜合兩方面原因，偏心率的增大會(huì)讓法向接觸力各向異性向切向接觸力各向異性轉(zhuǎn)化。

3 結(jié)論

基于筆者所在團(tuán)隊(duì)開發(fā)的非球顆粒開源離散元程序SudoDEM，開展了擴(kuò)展超級(jí)橢球的顆粒趨真建模，進(jìn)行了不同偏心率試樣的真三軸剪切。基于組構(gòu)張量表征各向異性，并通過二階傅里葉拓?fù)鋽M合圖和力鏈圖將各向異性可視化，著重分析了顆粒偏心率對(duì)顆粒介質(zhì)細(xì)觀次生各向異性演化的影響，得到以下主要結(jié)論：

（1）隨偏心率增大，細(xì)觀組構(gòu)各向異性不同程度地發(fā)展，共同承擔(dān)著宏觀抗剪強(qiáng)度的提升。同時(shí)，臨界配位數(shù)更大，滑動(dòng)接觸比例更高，接觸力網(wǎng)格的空間非均勻性也愈加顯著。這主要因?yàn)轭w粒偏心增強(qiáng)了顆粒間的互鎖與咬合。

（2）三軸剪切發(fā)展了次生各向異性。五種組構(gòu)各向異性表現(xiàn)各不相同，并在不同的應(yīng)變水平達(dá)到各向異性臨界狀態(tài)。其中，法向接觸力與接觸法向的空間分布呈“葫蘆”狀，前者的各向異性演化形態(tài)最接近強(qiáng)度演化形態(tài)，并始終占據(jù)最大權(quán)重；后者的各向異性權(quán)重次之（二者總占比約80%）；法向支向量空間分布為圓形，其各向異性可忽略不計(jì)；切向接觸力及切向支向量的空間分布呈頗似剪切帶的“雙葉草”狀。特別地，二者受偏心率影響較為敏感，在偏心顆粒試樣中二者對(duì)強(qiáng)度的貢獻(xiàn)不可忽略。

（3）剪切過程中幾何各向異性和力學(xué)各向異性相互轉(zhuǎn)化；隨偏心率增大，法向接觸力各向異性逐漸向切向接觸力各向異性轉(zhuǎn)化。