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

（華南理工大學(xué)亞熱帶建筑科學(xué)國家重點(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］，因此對各向異性的細(xì)觀尺度研究具有重要意義。Casagrande［6］依據(jù)各向異性成因?qū)⑵鋭澐譃樵飨虍愋院痛紊飨虍愋?，原生各向異性源于顆粒集合體在天然沉積過程中因顆粒級配及顆粒形狀導(dǎo)致的顆粒排列及接觸的空間差異，而次生各向異性則源于顆粒在外部荷載擾動(dòng)下的原生組構(gòu)重新排列。Oda等［7］進(jìn)一步地基于顆粒集合體的接觸法向與接觸力等細(xì)觀組構(gòu)張量給出了各向異性的定量表征，為從細(xì)觀尺度解析顆粒介質(zhì)諸多典型宏觀力學(xué)現(xiàn)象（例如，擠壓膨脹效應(yīng)［8］、糧倉效應(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)行跨尺度研究即為一條重要途徑。

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

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

1 離散元模型建立

1.1 趨真顆粒形狀模型

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

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

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

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

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

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

偏心率具體定義如下

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

特別地，本研究獨(dú)立于顆粒形態(tài)棱角度與長細(xì)比展開研究，即令ε1=ε2=1.0以及l(fā)x=ly=lz。同時(shí)將x、y、z軸三個(gè)方向的偏心率統(tǒng)一起來，即令ζ=ζ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］。顆粒接觸模型采用線性彈簧模型及庫倫滑移模型。擴(kuò)展超橢球接觸檢測及計(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］來保證不同顆粒形狀試樣的初始孔隙比一致為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ì)觀層面顆粒間接觸力與支向量來計(jì)算應(yīng)力張量σij。

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

平均主應(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ù)邊界墻的位移來計(jì)算對數(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)了顆?；ユi和咬合牢固性，使得試樣可承擔(dān)更大的剪應(yīng)力。Kozicki等［55］、Zhao等［36］也有過類似報(bào)道。另外，雖試樣的初始孔隙比一致，但試樣的相對密實(shí)度隨顆粒偏心率改變而不同。具體地，偏心率較大的試樣相對密實(shí)度較小，偏心率較小的試樣相對密實(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的顆粒對組構(gòu)穩(wěn)定無貢獻(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所示。可見各試樣力學(xué)平均配位數(shù)隨應(yīng)變發(fā)展呈指數(shù)衰減并最終達(dá)臨界值，即力學(xué)穩(wěn)定jamming態(tài)。根據(jù)Rothenburg等［58］和Zhao等［59］的相關(guān)研究，配位數(shù)的減小主要為密實(shí)的試樣在剪切過程中體積剪脹所帶來的塑性變形所致。另外，隨偏心率的增大，配位數(shù)下降幅度減小。所以，臨界配位數(shù)隨偏心率單調(diào)增加，對應(yīng)前文所述較高的臨界抗剪強(qiáng)度。可能原因?yàn)椋浩穆视?，顆粒間接觸面積更大，內(nèi)嵌愈加縝密，因此試樣配位數(shù)愈大，內(nèi)部組構(gòu)愈穩(wěn)定，為其擁有更高強(qiáng)度作出細(xì)觀解釋。Santamarina等［60］也指出非球顆粒需要更大的配位數(shù)來維持穩(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ì)算采用了庫倫滑移模型，摩擦發(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)接觸比例均隨偏心率增大而增大，但后者較前者對偏心率更為敏感。

圖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é)論。對此，Estrada等［63］曾作相關(guān)解釋：顆粒的偏心不規(guī)則性加劇了接觸互鎖現(xiàn)象，提高了顆?？剐D(zhuǎn)能力，進(jìn)而導(dǎo)致顆粒間接觸需通過相互滑動(dòng)來抵抗外荷載。最終，滑動(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

各試樣在剪切前的力鏈分布相近，因此僅針對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ù)量上雖無明顯變化，但其管徑膨脹，加劇了豎直向與水平向法向接觸力的各向異性。雖高偏心率試樣對應(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)張量對組構(gòu)各向異性進(jìn)行定量描述。

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

式中：φij為3×3二階對稱張量；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ì)心連線。由于顆粒形狀的影響，非球顆粒的支向量與接觸法向不重合。對此，以接觸法向?yàn)閰⒖枷祵⒅蛄糠纸鉃榉ㄏ蚍至亢颓邢蚍至?，法向支向量dn與切向支向量dt的組構(gòu)張量分別為

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

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

求得偏組構(gòu)張量后便可通過第二不變量A*來量化各向異性［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ù)可以用傅里葉級數(shù)展開。考慮到方向量的對稱性，奇數(shù)階張量對級數(shù)解無貢獻(xiàn)［67］，可采用二階傅里葉級數(shù)展開［15，65］。

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

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

由于剪切前試樣處于各向同性應(yīng)力狀態(tài)，5類組構(gòu)在剪切前也均處于初始各向同性狀態(tài)。圖8（a）中各向同性的弱力鏈以及圖9中5類組構(gòu)各向異性第二不變量的初始值均為零便作了證明，因此鑒于篇幅有限不對初始狀態(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)力方向）上一定程度地增長。然而，水平向則不然，甚至出現(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)較大幅度的增長。但同時(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)，偏心率對組構(gòu)各向異性的影響也各不相同。偏心率對接觸法向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的各向異性隨偏心率增大而大幅度增長。此觀點(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)系線的斜率較高，表明偏心率對二者的影響較大?？梢哉J(rèn)為，在偏心顆粒試樣中，切向接觸力ft、切向支向量dt對幾何各向異性奉獻(xiàn)較大，不應(yīng)忽視，這與其他學(xué)者的結(jié)論有所偏差［33］。

5類組構(gòu)各向異性對剪應(yīng)力比（宏觀抗剪強(qiáng)度）的貢獻(xiàn)各不相同。圖12展示了5類組構(gòu)各向異性的貢獻(xiàn)比重隨軸向應(yīng)變的演化。從中可見，Afn權(quán)重始終最大。這表明，法向接觸力fn的各向異性對抗剪強(qiáng)度的貢獻(xiàn)最大。同時(shí)，在圖9（b）中，Afn的演化形態(tài)與宏觀強(qiáng)度的演化形態(tài)始終最為接近，此現(xiàn)象也對上述觀點(diǎn)作出了佐證。文獻(xiàn)［18，24，31］也曾報(bào)道過類似的觀點(diǎn)。法向接觸力fn及接觸法向c各向異性總占比約80%，而其余三者權(quán)重均較小。這表明接觸法向c與法向接觸力fn各向異性對強(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)出先增加后減小的趨勢，Afn的貢獻(xiàn)比重則呈現(xiàn)出先減小后增大的趨勢。它們分別在約15%的軸向應(yīng)變處達(dá)到峰值與谷值。這說明，剪切過程中Ac的權(quán)重和Afn相互轉(zhuǎn)移。這意味著，剪切過程中幾何各向異性與力學(xué)各向異性相互轉(zhuǎn)化。其原因正如Zhao等［33］所述：剪切過程中因剪脹而不斷減小的平均配位數(shù)使得接觸分散性降低。

在偏心率的影響方面，從圖12得知，隨偏心率增大，Ac、Afn的相對權(quán)重減小，Aft的相對權(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)格中切向接觸力的各向異性相比法向接觸力各向異性更大。也即，較高偏心率帶來的滑動(dòng)接觸比例的增大會(huì)導(dǎo)致切向接觸力各向異性比重增大。綜合兩方面原因，偏心率的增大會(huì)讓法向接觸力各向異性向切向接觸力各向異性轉(zhuǎn)化。

3 結(jié)論

基于筆者所在團(tuán)隊(duì)開發(fā)的非球顆粒開源離散元程序SudoDEM，開展了擴(kuò)展超級橢球的顆粒趨真建模，進(jìn)行了不同偏心率試樣的真三軸剪切。基于組構(gòu)張量表征各向異性，并通過二階傅里葉拓?fù)鋽M合圖和力鏈圖將各向異性可視化，著重分析了顆粒偏心率對顆粒介質(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ì)；切向接觸力及切向支向量的空間分布呈頗似剪切帶的“雙葉草”狀。特別地，二者受偏心率影響較為敏感，在偏心顆粒試樣中二者對強(qiáng)度的貢獻(xiàn)不可忽略。

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