基于可破碎離散元法的堆石料應力變形及剪脹特性縮尺效應研究

徐 琨，楊啟貴，周 偉，馬 剛，黃泉水

（1.長江勘測規劃設計研究院，武漢 430010；2.水利部土石壩破壞機理與防控技術重點實驗室，南京 210029；3.武漢大學水資源與水電工程科學國家重點實驗室，武漢 430072）

0 引 言

堆石料作為混凝土面板堆石壩的主要筑壩材料，具有較寬的顆粒級配，最小的顆粒粒徑可小于0.1 mm，最大的顆粒粒徑可達1 000 mm 以上，國內一些面板堆石壩工程中采用的粒徑最大可達800～1 600 mm［1］。室內試驗是當前獲取堆石料力學變形參數的重要手段，《土工試驗規程》（SL237-1999）中規定試樣直徑應不小于試驗土料最大粒徑的5～6 倍，由于筑壩堆石料最大粒徑過大，受限于當前試驗設備和技術水平，只能對縮尺后的堆石料開展力學試驗。縮尺后堆石料與原級配堆石料所呈現出的力學變形特性差異，稱之為堆石料的縮尺效應。

當前眾多學者針對堆石料縮尺效應開展了諸多有益的研究。在力學特性的縮尺效應方面，文獻［2-8］通過研究發現不同尺寸試樣的抗剪強度一般呈現出隨最大粒徑的增大而降低，大尺寸試樣破碎率高于小尺寸試樣破碎率的規律；但也有文獻［9-11］研究發現堆石料力學特性呈現了趨勢相反的縮尺規律。

在變形特性縮尺效應方面，王繼莊［12］和酈能惠等［9］通過研究均發現體積彈性模量隨最大粒徑的增大而增加，而李翀等［3］研究發現體積彈性模量隨最大粒徑的增大而減小。褚福永等［13］和武利強等［14］研究表明，堆石料初始彈模和體積切線模量隨最大粒徑的增加而增大，初始泊松比隨最大粒徑的增大逐漸減小。Wei等［15］研究認為粗粒料試樣的壓縮模量隨最大粒徑的增大而增加。花俊杰等［16］研究認為大尺寸試樣的切線模量更高。王永明等［17］研究表明按相同相對密度制樣時，最大粒徑越大的試樣其體積模量和初始彈模越大；按相同干密度制樣時，試樣的變形參數隨最大粒徑的增加呈現出先減小后增大的非單調關系。武利強等［7］研究表明，同一干密度制樣時，粗粒料的變形模量隨最大粒徑的增大而減小。孔憲京等［18］基于超大型三軸儀研究了堆石料相變點處的體變縮尺規律，結果表明試樣尺寸越大，相變點處的體變也越大，也即剪縮性越強。

以上研究表明，縮尺效應不僅對堆石料力學特性有影響，同時對其變形特性也有顯著影響。當前，在堆石料縮尺規律方面取得了豐富的研究成果，在縮尺機理上有一定程度的認知，但由于堆石料縮尺規律的影響因素較多［4，7，15，19，20］，在多個影響因素的耦合作用下，堆石料力學變形特性的縮尺效應仍有待深入研究，縮尺引起的堆石料應力變形差異產生機理還有待明晰。在當前試驗條件和技術水平下，室內試驗一方面在試樣尺寸上有很大的限制，另一方面難以監測到試樣細觀層面的實時演化情況，無法非常有效的揭示堆石料縮尺效應產生的機理。

離散元方法作為近些年來逐漸興起的一種重要數值模擬研究手段，為從細觀尺度研究堆石料力學響應及其機理提供了有效的途徑。為進一步揭示堆石料縮尺效應產生機理，本文基于離散元方法開展堆石料縮尺效應研究。考慮剪脹特性是土體所具有的特殊性質，也是描述堆石料變形特性的關鍵指標之一，已有研究初步揭示了堆石料剪脹特性的縮尺規律，但縮尺產生的機理尚未有深入研究，仍不清晰。鑒于此，本文通過開展不同尺寸試樣的數值三軸剪切試驗，探討堆石料力學變形縮尺規律，并從宏、細觀層面深入分析試樣剪脹特性的縮尺規律，揭示其產生的細觀機理。

1 數值試驗準備

1.1 顆粒破碎的模擬

研究采用PFC3D［21］模擬常規三軸剪切試驗，模擬中使用“圓球”作為堆石料顆粒，從而消除顆粒形狀對試驗結果的影響。顆粒破碎是影響縮尺效應的主要原因之一［22，5］，研究中采用碎片替代法（Fragment Replacement Method，FRM）［23］模擬顆粒破碎。采用FRM 模擬顆粒破碎需要考慮兩大問題：①顆粒破碎準則；②碎片替換模式。本研究中顆粒破碎準則采用八面體剪應力準則［24］，碎片替換模式采用Apollonian 堆積態碎片替換模式［25］，可以較為真實的反映顆粒材料的力學響應特性，具體細節可參見文獻［8］，本文不再贅述。

1.2 數值試樣級配曲線

如圖1所示，本研究采用相似級配法對原級配縮尺，原級配堆石料最大顆粒粒徑（dmax）為600 mm，分別縮尺得到dmax為60、90和120 mm 的三條級配曲線。圖1中所示級配曲線由Tyler公式［26］［pi= 100(di/dmax)3-D］計算得到，其中，pi為小于粒徑di的所有顆粒的累計質量所占所有顆粒總質量的百分比，D為描述級配曲線分形特征的維數。考慮模擬效率，試樣中顆粒的數目不能過多，同時試驗所選取的顆粒材料的級配寬度也不能太寬；經反復模擬測試，最終取D為2.0，顆粒最小粒徑（dmin）為10 mm，對于粒徑10 mm 以下的顆粒以等質量方式替換為粒徑10 mm的顆粒。

1.3 數值試樣制備

研究采用的試樣徑徑比為5，在圓柱形空間內采用圖1 中所示的3 種縮尺后級配，通過粒徑膨脹法生成三軸試樣。制得試樣的尺寸為300×600、450×900 和600×1 200（直徑×高度），單位均為mm，各尺寸試樣中dmax分別為60、90和120 mm。

圖1 原級配曲線及其縮尺后的級配曲線Fig.1 Prototype gradation curve and the scaled gradation curves

由于3 種尺寸試樣的級配不同，為使不同試樣初始密度狀態達到基本相同的水平，制樣方法和參數參考文獻［27］。最終制得100 kPa 圍壓下3 種尺寸的初始試樣（圖2），初始試樣的孔隙率分別為30.65%（dmax=60 mm）、28.23%（dmax=90 mm）和27.03%（dmax=120 mm）。結果表明，試樣的最大粒徑越大，采用相同初始密度狀態控制制樣得到的試樣密度越大，與文獻［19，28］研究結果一致。

圖2 不同縮尺級配的初始試樣示意圖Fig.2 Illustration of the initial samples with different particle size distribution

研究中采用的接觸模型及離散元參數參考文獻［29，30，8］選取，如表1所示。

表1 離散元模擬參數Tab.1 Input parameters for DEM simulation

研究選取0.4，0.8 和1.2 MPa 三種試驗圍壓，對各試樣以最大顆粒粒徑和試驗圍壓為標識編號，如D60CP400 表示最大顆粒粒徑為60 mm，試驗圍壓為0.4 MPa（400 kPa）的試樣，其中D代表“Diameter”，CP 代表“Confining Pressure”，其他試樣編號均遵循該規則。

2 數值試驗結果與分析

2.1 力學特性的縮尺規律

三軸數值剪切試驗的應力-應變演化曲線如圖3 所示。由圖3 可知，不同尺寸試樣間的應力-應變曲線存在顯著差異，主要表現為尺寸越大試樣的峰值強度越小，且圍壓越高時差異越明顯。提取各試樣峰值摩擦角φpeak［sinφpeak=(σ1-σ3)/(σ1+σ3)，其中σ1為大主應力，本文研究中也為軸向應力，σ3為小主應力，本文研究中也即圍壓］表征試樣抗剪強度，得到峰值摩擦角與dmax的關系曲線如圖4所示。可知，試樣峰值摩擦角隨試樣最大粒徑的增大而減小，也即試樣尺寸越大，試樣抗剪強度越低，且圍壓越大時越明顯。以上規律與室內試驗研究結果［31，32］一致，表明研究中開展的三軸數值剪切試驗是合理有效的。

圖3 應力-應變演化曲線Fig.3 Stress-strain curves

圖4 峰值摩擦角與dmax的關系曲線Fig.4 Relationship between peak friction angle and dmax

根據De Mello［33］對粗粒料提出的冪函數非線性公式［τ=A(σn)b，A和b為擬合參數，τ為粗粒料的剪切強度，σn為正應力］來擬合得到研究中各組試樣的剪切強度包絡線，如圖5 所示。可知，擬合得到的b值均小于1，試樣剪切強度包絡線為一條斜率逐漸變緩的曲線，即隨正應力的增加試樣剪切強度的增加幅度有減小的趨勢，且最大粒徑越大的試樣，其增幅減小越明顯。

圖5 三軸試樣剪切強度包絡線Fig.5 Shear strength envelope of triaxial tests

2.2 變形特性的縮尺規律

三軸數值剪切試驗的體變-應變演化曲線如圖6所示，本研究中取體變剪縮為負，剪脹為正。由于生成的初始試樣是較為密實的試樣，加載過程中試樣呈現出的剪縮量很小。由圖6 可知，不同尺寸試樣間的體變演化存在差異，尺寸越大的試樣在演化過程中呈現出更弱的剪脹特性，且圍壓越高時這一差異越明顯。

圖6 不同圍壓下試樣體變-應變曲線Fig.6 Volumetric strain response under different confining pressures

為更直觀對比不同試驗圍壓下各尺寸試樣剪脹剪縮的差異，提取相變點處和10%剪應變處試樣的體變量的絕對值如圖7所示。其中，相變點［34］為剪脹性材料所特有的性質，即體積由體縮變為體脹的特征點，相變點處體變量也即試樣的最大剪縮量。由圖7（a）可知，10%剪應變處，圍壓越高的試樣，剪脹量越小，且試樣剪脹量隨著最大粒徑的增大而降低；圍壓越高時，隨著最大粒徑的增大，剪脹量降低程度越大，如：dmax由60 mm增大至120 mm，圍壓為0.4 MPa 時，剪脹量降低約15.8%，圍壓為1.2 MPa 時，剪脹量降低約66.2%。由圖7（b）可知，相變點處，試驗圍壓越高，試樣的剪縮量越大；試樣剪縮量隨著最大粒徑的增大而增加。以上結果表明，試樣尺寸越大，剪縮特性越強。

圖7 試樣剪脹特性隨最大粒徑的變化關系Fig.7 Relationship of the shear dilatancy characteristics with dmax

2.3 堆石料剪脹縮尺效應產生的原因分析

現有研究指出顆粒破碎是導致堆石料力學響應縮尺效應的關鍵因素，試樣尺寸越大，顆粒破碎率越高，抗剪強度越弱［5］，而導致這一現象的細觀力學機理則是由于剪切過程中，大尺寸試樣較小尺寸試樣產生了更多的力學不穩定顆粒導致的［8］。

對于堆石料剪脹特性，顆粒破碎同樣也是其縮尺規律的關鍵因素。顆粒材料的剪脹主要來源于顆粒錯動、轉動等重新排列后空隙的增加，而顆粒破碎會產生更為細小顆粒碎片，細小顆粒碎片可填充于大顆粒間的空隙中，從而減弱了試樣的剪脹效應［35］；顆粒破碎率越高，則剪脹效應減弱地越多，從而造成了尺寸越大試樣的剪縮特性越強。但受限于與室內試驗監測手段的限制，這一分析目前尚缺少試驗數據的支撐，細觀機理未得到深入分析，為此借助數值試驗的優勢開展深入研究。

（1）有效配位數演化規律。顆粒配位數是散粒體結構特征的重要表征參數之一，表示與目標顆粒發生接觸的顆粒數目［36，37］。Thornton［38］在研究中將配位數小于2 的顆粒為稱之為懸浮顆粒，并將試樣中顆粒配位數小于2 的懸浮顆粒除去后得到的試樣平均配位數定義為有效配位數CNeff，可表示為：

式中：C為試樣中各顆粒間接觸的總數；Ntot、N0和N1分別為試樣中顆粒的總數、配位數為0 的顆粒數目和配位數為1 的顆粒數目。

縮尺試樣有效配位數演化情況如圖8 所示。由圖8 可知，不同圍壓下，各尺寸試樣的有效配位數均呈現出相似的演化規律，表現為：有效配位數隨著試樣剪切先呈現出陡降，之后緩慢減小，隨著剪切繼續，有效配位數逐步達到基本穩定狀態，這與Duran 等［39］和Gu 等［40］的研究結果相似。同時，也發現各尺寸試樣有效配位數在數值上存在較為明顯的差異，表現為：同一剪切位移時，有效配位數隨試樣最大顆粒粒徑的增大而降低；結果表明，剪切試驗中，試樣尺寸越大，其承力結構越稀疏。

圖8 試樣有效配位數演化曲線Fig.8 Evolution curves of effective coordination number

（2）各粒徑組有效配位數。有效配位數的演化結果展示了不同尺寸試樣間承力結構的整體性差異，但不足以解釋堆石料試樣剪脹特性差異產生的原因。為深入分析，將試樣中顆粒按粒徑分組，然后求出各粒徑組內顆粒的有效配位數，得到各粒徑組有效配位數分布。分組規則為，粒徑小于等于10 mm 的為一組，10 到20 mm 的為第二組，依次類推；所有分組均為左開右閉區間，顆粒最大粒徑為60、90 和120 mm 的試樣分別分出6、9和12組。

初始試樣的各粒徑組有效配位數直方圖如圖9 所示。可知，不同尺寸試樣中各粒徑組有效配位數均呈現隨顆粒粒徑的增大而升高的趨勢，但最大粒徑組的有效配位數處于同一量級，相差不大。對于同一粒徑組，dmax越大的試樣，其粒徑組有效配位數越低，如（50，60］粒徑組時，顆粒最大粒徑為60、90 和120 mm 的試樣其粒徑組有效配位數分別為25.09、11.37和6.48。

將粒徑組和有效配位數歸一化處理，得到初始試樣各粒徑組有效配位數在半對數坐標下的分布，如圖10 所示。由圖可知，雖然標準化后的各粒徑組有效配位數分布有一些差異，但大體呈現出較為相似的分布，表明采用相似級配法縮尺得到的各尺寸試樣在一定程度上保持初始承力結構的相似性。

圖10 歸一化的初始試樣各粒徑組有效配位數分布Fig.10 Normalized effective coordination number distribution of each particle size group of the initial sample

圖11所示為剪切前后各粒徑組有效配位數分布對比情況，各試驗圍壓下分布規律相似，僅列出0.4 和1.2MPa 圍壓下的結果。由圖11可知，不同尺寸試樣的各粒徑組有效配位數分布在剪切后（ε1= 10%）呈現出不同的演化趨勢。剪切后，小尺寸試樣（dmax= 60 mm）各粒徑組有效配位數均有所降低，整體降低明顯；隨試樣尺寸的增大，其各粒徑組有效配位數降低程度減弱，并有部分粒徑組下的有效配位數高于剪切前的值。此外，不同試驗圍壓下，試樣剪切前后的各粒徑組有效配位數變化情況有所差異，表現出圍壓越大且試樣尺寸越大時，剪切后的各粒徑組有效配位數較剪切前的上抬越明顯。

圖11 試樣各粒徑組有效配位數分布Fig.11 Effective coordination number distribution of each particle size group

另開展一組不考慮顆粒破碎的三軸剪切試樣做對比研究，得到試樣剪切前后各粒徑組有效配位數分布如圖12所示，僅展示1.2 MPa 試驗圍壓。對比圖12 與圖11（b）可知，顆粒破碎對試樣剪切后各粒徑組有效配位數分布產生了顯著影響，顆粒破碎使剪切后的各粒徑組有效配位分布發生上抬，且試樣尺寸越大時上抬程度越大。采用Marsal［2］對試樣顆粒破碎率的定義，提取試樣顆粒破碎率如圖13 所示。由圖可知，試樣尺寸越大，顆粒破碎率越大，且試驗圍壓越大時，顆粒破碎率越大。綜上，可以證明各粒徑組有效配位數分布的上抬與顆粒破碎率是成正相關的，即試樣顆粒破碎率越大，各粒徑組有效配位數分布上抬越多。

圖12 試樣各粒徑組有效配位數分布（顆粒不破碎）Fig.12 Effective coordination number distribution of each particle size group（uncrushable particle）

圖13 顆粒破碎率BmFig.13 Particle breakage factor Bm

有效配位數在一定程度上表征試樣的密實程度，配位數越高表明顆粒接觸越多，試樣越密實。當試樣顆粒發生破碎后，碎片填充于顆粒間的空隙中，隨著剪切的進行逐漸充滿密實，試樣的承力結構變得密實，也即有效配位數增大；顆粒破碎越多時，則試樣填充越密實，有效配位數增大越多，由此揭示了不同尺寸試樣各粒徑組有效配位數分布剪切前后演化差異產生的原因。

此外，有效配位數的減小在一定程度上代表試樣的剪脹，反之代表試樣的剪縮。為更清晰展現剪切前后試樣各粒徑組有效配位數的變化情況，定義各粒徑組有效配位數分布差ΔCNeff=CNeff|ε1=10%-CNeff|ε1=0%。其中，CNeff|ε1=0%、CNeff|ε1=10%分別代表剪切前、后的各粒徑組有效配位數分布，ΔCNeff為正值代表剪切后發生了剪縮，負值代表剪切后發生了剪脹。對比顆粒不破碎（圖12）和顆粒破碎［圖11（b）］兩種情況，得到各粒徑組有效配位數分布差如圖14 所示。可知，不同尺寸試樣間ΔCNeff的有明顯差異，同時考慮顆粒破碎的試樣較不考慮顆粒破碎的試樣的ΔCNeff均在較大粒徑區間有明顯增大，試樣尺寸越大則ΔCNeff增大越明顯，也即試樣發生了更多程度的剪縮。同時，統計得到各試驗圍壓下不考慮顆粒破碎試樣剪切后（ε1= 10%）的體變均在+5.00±0.1%范圍內，表明不考慮顆粒破碎時各尺寸試樣的剪脹量基本相同。綜上可知，一方面，顆粒破碎促使ΔCNeff向增大的方向演化；另一方面，尺寸越大試樣的大粒徑區間的ΔCNeff增多越明顯，降低了試樣的剪脹量，使得試樣尺寸越大，剪脹性越弱；由此揭示了試樣剪脹縮尺規律發生的原因。

圖14 各粒徑組有效配位數分布差Fig.14 Difference of effective coordination number distribution of each particle size group

3 結 論

本文基于可破碎離散元法，從宏、細觀層面開展了堆石料應力變形及剪脹特性縮尺效應研究，得到了如下主要結論：

（1）數值模擬試驗較好反映了堆石料應力變形的縮尺規律。試樣尺寸越大，抗剪強度越低，剪縮特性越強，最終的剪縮量越大；試驗圍壓越高，縮尺效應對堆石料抗剪強度、剪脹特性的影響越大，縮尺現象越顯著。

（2）采用相似級配法縮尺得到的不同尺寸試樣一定程度上保持著初始承力結構的相似性，但尺寸越大的試樣，其承力結構越稀疏。試樣尺寸越大，各粒徑組有效配位數分布的降低程度減弱，并有部分粒徑組下的有效配位數高于剪切前的值，圍壓越大且試樣尺寸越大時，這一現象越明顯。

（4）試樣的各粒徑組有效配位數分布的演化與試樣顆粒破碎率成正相關，試樣顆粒破碎率越大，各粒徑組有效配位數分布上抬越多；通過對比研究，揭示了不同尺寸試樣間各粒徑組有效配位數分布差（ΔCNeff）演化的差異是試樣剪脹特性縮尺規律產生的細觀機理之一。 □