堆石料宏細觀力學特性離散元分析

摘要：堆石料的力學特性直接影響堆石壩工程的設計、施工和運行，對壩體變形及安全有重要影響。基于大石峽混凝土面板堆石壩工程，采用離散單元法（DEM）建立了接近于3BA區堆石料真實顆粒形狀的橢球顆粒試樣，開展了一系列常規三軸數值試驗，研究了3BA區堆石料的宏、細觀力學特性。通過不同圍壓下的宏觀應力應變演化曲線，獲得堆石料的宏觀應力、變形指標。采用平均接觸力將接觸網絡劃分為強、弱2個子網絡，分析強、弱接觸網絡的配位數分布與拓撲演化規律。結果表明：割線模量、泊松比、峰值內摩擦角與圍壓呈線性變化，剪脹角與圍壓呈非線性變化；顆粒體系內接觸的拓撲演化主要發生在弱接觸網絡，進而影響堆石料的宏觀體積變化，強接觸網絡內的接觸保持不變，為堆石料提供穩定的傳力通道。研究結果可用于指導實際工程中不同圍壓下宏觀力學指標的選取。

關 鍵 詞：離散元； 橢球顆粒試樣； 宏觀力學指標； 強弱接觸網絡； 拓撲； 宏細觀特性

中圖法分類號： TU431 文獻標志碼： A DOI：10.16232/j.cnki.1001-4179.2024.02.026

0 引 言

堆石料由離散的摩擦性顆粒材料組成，如碎石、砂礫石、塊石等，由于其對地形和地質條件適應性好、可就地取材，在大型水電、港口、交通等基礎設施的建設中被廣泛應用。以堆石壩為例，堆石料的力學特性直接影響堆石壩工程的設計、施工和運行，其力學性能的演化對壩體變形及安全有重要影響［1］。堆石料本質上是離散顆粒材料，連續變形方法（如有限元法）只能從宏觀上分析堆石體的應力變形特征；對堆石料進行宏觀唯象建模不足以反映其所有力學特性。從細觀角度出發，研究堆石料宏細觀力學特性及影響關系，從而深入理解其基本的物理力學特性，是目前認識、應用堆石料的重要立腳點。

從細觀層次揭示堆石料力學特性的方法，目前運用最廣泛的是離散單元法（discrete element method，DEM）［2-4］。該方法通過細觀數值模擬方法可以實時監測堆石料體系內顆粒的位置、速度、角速度、接觸力等幾何和力學特性的細觀演化［5-8］，研究其細觀組構等力學參數對宏觀力學響應的影響，還可以高效開展大量敏感性分析。邵磊等［9］采用離散單元法模擬不同圍壓下堆石料的大三軸排水剪切試驗，得到的應力、應變曲線與室內試驗一致，證明離散單元法可以較好地模擬堆石料的大三軸試驗過程。

堆石料一般來源于巖體的爆破和破碎，顆粒形狀是不規則的，呈現角狀或次角狀［10］。在常規的DEM模擬中，為了優化接觸檢測的計算效率，將堆石料理想化為圓形剛性圓盤或球體的集合體，但這種理想化的處理方法無法再現真實顆粒的一些幾何行為，如顆粒間的互鎖和滾動阻力。有些學者提出可以用抗轉動模型表征復雜的顆粒形狀，同時提高計算效率，但抗轉動模型模擬可以在多大程度上考慮顆粒形狀的影響尚不清楚。在研究顆粒材料微觀組構和宏觀力學行為的聯系時，Zhao等［11］發現抗轉動阻力會產生不切實際的各向異性，建議謹慎使用抗轉動模型。鄒宇雄等［12］指出轉動摩擦和復雜形狀顆粒分別通過限制顆粒轉動、增大咬合作用而增加顆粒間的承載能力，兩者的作用機制有本質差別。還有一些學者采用顆粒簇［13］、多面體［14］、橢球［15-16］等表示不規則的顆粒形狀進行數值模擬研究，其中橢球顆粒是描述各向異性堆石料的最簡單形狀，而且還可以排除棱角、表面粗糙度等因素的干擾。因此，本文采用橢球顆粒進行堆石料的宏、細觀力學性能研究。

本文基于大石峽混凝土面板堆石壩工程，通過統計分析3BA區堆石料的物理、力學特性，建立接近堆石料真實形狀且無滾動摩擦的橢球顆粒數值模型，并進行細觀參數率定；在橢球顆粒試樣的宏觀應力變形與室內三軸實驗一致的前提下，進行堆石料的數值剪切試驗研究，分析圍壓與堆石料宏觀應力、變形指標的關系，探討3BA區堆石料的細觀力學特性。

1 DEM數值試驗

1.1 大石峽堆石料顆粒形狀分析

大石峽水電站的擋水建筑物為混凝土面板砂礫石壩，最大壩高247 m。將現場取樣的3BA區堆石料放入篩分儀進行篩分處理，顆粒的粒徑分為10～20 mm，20～30 mm，30～40 mm，40～50 mm，50～60 mm 5組，每組顆粒取40個。采用先臨三維Einscan-S掃描儀對每組顆粒進行掃描（見圖1）。考慮到計算量和試驗精度的要求，采用MESHLAB軟件，在保留顆粒典型特征信息的基礎上對顆粒的三維信息進行簡化。通過提取、分析5組顆粒試樣（共400個顆粒）的短、中、長軸長度（S、M、L）等形態特征（見表1），發現采用顆粒形狀為S∶M∶L=1.0∶1.5∶2.0的橢球顆粒能夠很大程度上代替3BA區砂礫料的真實顆粒形狀。

1.2 制樣與加載

采用LIGGGHTS生成橢球顆粒的數值試樣［17］，試樣中橢球顆粒形狀為S∶M∶L=1.0∶1.5∶2.0，如圖2所示。制樣時，在尺寸為300 mm×300 mm×300 mm的立方體試樣內，按照3BA區堆石料的顆粒級配形成相應粒徑大小的顆粒個數，總數70 000個，最大粒徑60 mm，試樣級配曲線如圖3所示。其中，橢球顆粒大小采用與圓球顆粒相等體積球體的等效半徑表示。固結時采用各向同性壓縮的方法，6面剛性墻體緩慢向中心移動直至達到目標圍壓。固結時，粒間摩擦系數設為0.03，形成較為密實的顆粒試樣，初始孔隙比為0.20，圍壓為0.3 MPa，如圖4所示。

采用Modified Hertz-Mindlin接觸模型［18］計算顆粒間的相互作用力。在常規三軸壓縮過程中，頂部和底部的剛性邊界以恒定速度分別向下和向上移動，其余4個側向邊界單獨移動。通過控制慣性數I小于0.001，確定大主應變方向的加載速率，保證試樣為準靜態加載［19］；同時保持小主應力σ3不變，恒定在目標圍壓；控制試樣的中主應力σ2=σ3。

1.3 細觀參數的率定

試樣細觀參數的選取對DEM數值試驗的準確性和可靠度至關重要。數值模擬時，試樣的細觀參數參照南京水利科學研究院進行的3BA區堆石料室內三軸剪切試驗的數據確定，密度為2 380 kg/m3，泊松比為0.15，滑動摩擦系數為0.4。恢復系數與顆粒間接觸的法向和切向阻尼系數相關，取值為0.95［20-21］，比較符合堆石料碰撞與回彈的相關特性［22］，具體細觀參數如表2所列。實驗模擬得到的模型偏應力-軸應變曲線與室內試驗結果吻合度較高，且能準確模擬出室內試驗剪縮和剪脹狀態轉換對應的軸應變，如圖5所示。

2 宏觀力學響應及力學指標

2.1 加載過程中的應力應變

開展不同圍壓下的三軸壓縮試驗，圍壓分別為0.3，0.6，0.9，1.2 MPa。圖6為不同圍壓下，試樣的偏應力、體應變隨軸應變的演化過程。隨著圍壓的增大，加載伊始，試樣偏應力的增大速率增加，達到的偏應力峰值增大，且偏應力達到峰值時對應的軸向應變增大。宏觀應變表現為：圍壓越高，試樣剪脹程度越小，且試樣剪縮轉變為剪脹所對應的軸向應變越大。在高圍壓下，試樣主要表現為剪縮，是因為較大的圍壓限制了顆粒爬升、翻越相鄰顆粒的能力，導致試樣體積由剪縮變為剪脹狀態對應的軸向應變較大。

2.2 不同圍壓下的強度、變形指標

堆石料常用的宏觀力學指標有割線模量E50、泊松比ν、峰值內摩擦角φ和剪脹角ψ等，這些強度和變形指標可以由常規三軸加載試驗的宏觀力學響應曲線計算得到［23］，具體定義如圖7所示。

依據不同圍壓下試樣在常規三軸加載試驗的宏觀應力變形曲線，得到3BA區堆石料相應的強度和變形指標，如表3所列；依據其變化規律擬合曲線，如圖8所示。割線模量、泊松比、峰值內摩擦角與圍壓呈線性變化，剪脹角與圍壓呈非線性變化。隨著圍壓的增加，割線模量增大，兩者之間的線性關系可以表示為E50=kσ3Pa+ΔE（1）式中：k、ΔE均為模型參數，k=16.32，ΔE=192.69 MPa；σ3為圍壓；Pa為標準大氣壓。

峰值內摩擦角φ、泊松比ν與圍壓也可以用類似帶截距的線性方程表達，分別為φ=kφσ3Pa+Δφ（2）式中：kφ =-0.16；Δφ=50.65°。ν=kνσ3Pa+Δν（3）式中：kν =-0.023；Δν=0.47。

剪脹角ψ隨著圍壓σ3的增加而減小，兩者之間呈非線性關系，可表示為ψ=Δψσ3Pan（4）式中：Δψ、n為模型參數，表征非線性程度，本次研究Δψ=37.15°，n=0.42。在高圍壓條件下，顆粒間的運動受限會影響堆石料宏觀強度及變形特性，泊松比、峰值內摩擦角和剪脹角隨之減小。對圍壓與強度、變形指標的關系進行分析，有助于選取不同圍壓下的宏觀力學指標，更好地指導實際工程。

3 細觀組構張量與拓撲演化

3.1 組構張量與強接觸網絡

堆石料的宏觀力學特性與加載過程中細觀組構的演化規律密切相關［24］。本文采用Oda ［25］提出的組構張量表達式量化接觸法向的方向分布。Φij=1NcNcc=1ninj（5）式中：n是沿著顆粒間接觸面法線方向的單位矢量；Nc為顆粒體系中的接觸總數；Φii=1。

Thornton和Antony［26］在三維顆粒體系中定義傳遞大于平均法向接觸力的接觸為強接觸網絡。強接觸網絡的組構張量定義為Φsij=1NsNss=1nsinsj（6）式中：nis是單位法向量在第s個強接觸處的第i個分量；Ns是顆粒體系中的強接觸總數。

根據圍壓為0.3 MPa的常規三軸試驗結果分析加載過程中堆石料細觀組構的演化，如圖9所示。強接觸網絡具有較高的各向異性，隨著軸應變的增加，總接觸網絡的Φ1逐漸增大，與偏應力（σ1-σ3）演化較為一致，并在軸應變為2.5%時達到峰值狀態；而Φ1s的峰值所對應的軸應變早于Φ1，說明強接觸網絡在受外荷載作用時反應迅速。Huang等［27］的研究表明張量的主值可以表征顆粒在主方向的聚集度，即試樣內部大主應力方向的接觸聚集度高于小主應力方向。

3.2 接觸網絡的拓撲演化

堆石料體系內部是由顆粒間接觸連接形成的復雜網絡結構，為了更深入地研究堆石料體系內部的連通性，圖10給出了圍壓為0.3 MPa時常規三軸加載試驗中試樣在不同加載狀態下（軸應變為0.5%，1.5%，2.5%，5.0%，7.0%）的配位數分布情況，圖中Pn表示具有n個接觸的顆粒個數，n表示顆粒與周邊顆粒的接觸數。Bratberg等［28］和Liu等［29］分別對二維和三維狀態下的試樣結構進行表征，認為配位數越少的顆粒各向異性程度高，越不穩定，反之配位數越多的顆粒各向異性程度低，越穩定。初始狀態，試樣內P7的顆粒個數最多，隨著加載的進行，整個試樣的配位數分布逐漸左移，說明試樣內顆粒的接觸不斷丟失。在達到峰值狀態以后，配位數的分布穩定在P5的顆粒個數最多，并且隨著軸應變的增大，P5的個數逐漸減小，同時也反映出試樣細觀平均接觸數隨軸應變增大而逐漸減小的變化規律，及宏觀體應變由剪縮變為剪脹的特征。

顆粒體系內部的接觸連接形成的復雜空間網狀結構，稱之為接觸網絡，接觸網絡結構的變化與宏觀應力、變形密切相關。在二維顆粒體系中，許多學者借助拓撲分析，如力環、簇、力鏈結構［29-30］的概念深入了解局部或者整體結構演化以及其對顆粒材料的影響。一般而言，顆粒體系內部的顆粒在外荷載作用下會發生重排，在加載過程中的每一步都存在接觸的新生和消失。這里將顆粒體系內的接觸分為3類，消失（loss）的接觸、新生（gain）的接觸以及不變（constant）的接觸。在加載過程中，對于任意加載步，比如第N步的顆粒體系來說，存在于第N步但不存在于第N-1步的接觸為新生的接觸，存在于第N步但是不存在于第N+1步的接觸為消失的接觸，剩下沒有改變的接觸為不變的接觸。實際加載試驗中，每0.02 s的時間間隔，軸應變加載幅度約為0.04%。圖11給出試樣在加載過程中總接觸網絡新生、消失和不變的接觸所占的比例，分別用Rg，all、Rl，all和Rc，all表示。在試樣中幾乎有98%的接觸保持不變，消失和新增的接觸占比僅為2%。

新生、消失和不變的接觸個數在強接觸網絡中分別表示為Ng，s、Nl，s和Nc，s，在弱接觸網絡中分別表示為Ng，w、Nl，w和Nc，w表示。由圖11可知，試樣中新生和消失的接觸占比非常小，強、弱接觸網絡中大部分為不變的接觸。通過對比強、弱接觸網絡中新生和消失的接觸（見圖12），發現強接觸網絡中幾乎沒有接觸的新生和消失，絕大部分的拓撲變化都發生在弱接觸網絡。在軸應變2.5%～5.0%之間，弱接觸網絡中新生和消失的接觸處于此消彼長的恒定狀態，說明達到峰值狀態后的這一加載階段試樣處于彈性階段。而在軸應變2.5%（亦即應力峰值狀態）之前，消失的接觸數大于新增的接觸數，試樣發生剪縮變化，同時增大體系內的各向異性以抵抗外荷載的施加；在軸應變5.0%之后，新增的接觸數大于消失的接觸數。因此，在顆粒體系內，弱接觸網絡內的接觸在加載過程中不斷發生新生和消失的更替，影響宏觀體積的變化。強接觸網絡內的接觸保持不變，為試樣提供穩定的傳力通道，以抵抗外荷載的施加。

4 結 論

本文基于大石峽混凝土面板堆石壩工程，采用離散單元法建立了接近3BA區堆石料真實形狀的橢球試樣并開展常規三軸剪切試驗，研究了堆石料的宏、細觀力學特性，探討了顆粒體系內細觀拓撲的演化機制。主要結論為：

（1） 在加載過程中，隨著圍壓的增大，宏觀堆石料的偏應力峰值增大，進入剪脹狀態更緩慢，且后期剪脹程度更不明顯。

（2） 宏觀力學指標中，割線模量、泊松比、峰值內摩擦角與圍壓呈線性變化，剪脹角與圍壓呈非線性變化。

（3） 加載過程中，弱接觸網絡內的接觸不斷發生新生和消失的更替，影響顆粒體系宏觀體積變化，強接觸網絡內的接觸保持不變，為顆粒體系提供穩定的傳力通道，抵抗外荷載的施加。

本文的分析仍然處于一個初步階段，后續可以進一步探討不同圍壓對顆粒體系細觀拓撲演化的影響，并嘗試尋找不同圍壓條件下，峰值應力剪應變、脹縮轉化點剪應變、峰值摩擦角、剪脹角等關鍵工程參數與堆石料接觸特性的宏細觀聯系，進一步為本構模型的構建提供參考。

參考文獻：

［1］陳生水.特高土石壩建設與安全保障的關鍵問題及對策［J］.人民長江，2018，49（5）：74-78.

［2］劉波，劉鵬，馬剛，等.顆粒形狀對粗粒土密實度的影響研究［J］.人民長江，2022，53（4）：155-162.

［3］孫晨，韓文喜，王昊.不同粒徑對粗粒土力學參數影響的研究［J］.人民長江，2018，49（10）：97-103.

［4］劉嘉英，馬剛，周偉，等.基于離散元的顆粒材料三維臨界狀態與剪脹特性研究［J］.水利學報，2017，48（9）：1107-1117.

［5］楊舒涵，周偉，馬剛，等.粒間摩擦對巖土顆粒材料三維力學行為的影響機制［J］.巖土工程學報，2020，42（10）：1885-1893.

［6］MIRGHASEMI A A，MATYAS E L，ROTHENBURG L.Influence of particle shape on engineering properties of assemblies of two-dimensional polygon-shaped particles［J］.Geotechnique，2002，52（3）：209-217.

［7］DELUZARCHE R，CAMBOU B.Discrete numerical modelling of rockfill dams［J］.International Journal for Numerical and Analytical Methods in Geomechanics，2006，30（11）：1075-1096.

［8］董啟朋，姚海林，詹永祥.基于顆粒流模型的顆粒材料宏-細觀力學研究［J］.人民長江，2016，47（4）：68-73.

［9］邵磊，遲世春，賈宇峰.堆石料大三軸試驗的細觀模擬［J］.巖土力學，2009，30（增）：239-243.

［10］ZHOU W，MA G，CHANG X，et al.Influence of particle shape on behavior of rockfill using a three-dimensional deformable DEM［J］.Journal of Engineering Mechanics，2013，139（12）：1868-1873.

［11］ZHAO S，EVANS T M，ZHOU X.Shear-induced anisotropy of granular materials with rolling resistance and particle shape effects［J］.International Journal of Solids and Structures，2018，150：268-281.

［12］鄒宇雄，陳遠，馬剛，等.抗轉動對顆粒材料組構特性的影響研究［J］.巖土力學，2020，41（8）：1-11.

［13］SHI D D，WANG W，XUE J F，et al.EM modelling of hollow cylinder shear test on granular materials under combined compression-shear stress using non-spherical particles［J］.Journal of Hydraulic Engineering，2019，50（9）：1052-1062.

［14］ZHAO S W，ZHOU X W，LIU W H.Discrete element simulations of direct shear tests with particle angularity effect［J］.Granular Matter，2015，17（6）：793-806.

［15］NG T.Particle shape effect on macro-and micro-behaviors of monodisperse ellipsoids［J］.International Journal for Numerical and Analytical Methods in Geomechanics，2009，33（4）：511-527.

［16］NG T，GE L.Packing void ratios of very dense ternary mixtures of similar ellipsoids［J］.Granular Matter，2020，22（2）：53-67.

［17］KLOSS C，GONIVA C.LIGGGHTS—open source discrete element simulations of granular materials based on Lammps［J］.Supplemental Proceedings：Materials Fabrication，Properties，Characterization，and Modeling，2011，2：781-788.

［18］PODLOZHNYUK A，PIRKER S，KLOSS C.Efficient implementation of superquadric particles in discrete element method within an open-source framework［J］.Computational Particle Mechanics，2017，4（1）：101-118.

［19］ZHAO S W，ZHOU X W.Effects of particle asphericity on the macro-and micro-mechanical behaviors of granular assemblies［J］.Granular Matter，2017，19（3）：1-18.

［20］劉嘉英，周偉，馬剛，等.顆粒材料三維應力路徑下的接觸組構特性［J］.力學學報，2019，51（1）：26-35.

［21］ZHOU W，LIU J，MA G，et al.Three-dimensional DEM investigation of critical state and dilatancy behaviors of granular materials［J］.Acta Geotechnica，2017，12：527-540.

［22］IMRE B，RBSAMEN S，SPRINGMAN S M.A coefficient of restitution of rock materials［J］.Computers and Geosciences，2008，4（34）：339-350.

［23］馬剛.堆石體的連續-離散耦合分析方法與宏細觀力學特性［D］.武漢：武漢大學，2014.

［24］史旦達，周健，劉文白，等.砂土直剪力學形狀的非圓顆粒模擬與宏細觀機理研究［J］.巖土工程學報，2010，32（10）：1557-1565.

［25］ODA M.Fabric tensor for discontinuous geological materials［J］.Soils and Foundations，1982，22（4）：96-108.

［26］THORNTON C，ANTONY S J.Quasi-static deformation of particulate media［J］.Philosophical Transactions of the Royal Society of London.Series A：Mathematical，Physical and Engineering Sciences，1998，356（1747）：2763-2782.

［27］HUANG X，HANLEY K J，O′SULLIVAN C，et al.DEM analysis of the influence of the intermediate stress ratio on the critical-state behaviour of granular materials［J］.Granular Matter，2014，16（5）：641-655.

［28］BRATBERG I，RADJAI F，HANSEN A.Dynamic rearrangements and packing regimes in randomly deposited two-dimensional granular beds［J］.Physical Review E，2002，66（3）：31303.

［29］LIU J，ZHOU W，MA G，et al.Strong contacts，connectivity and fabric anisotropy in granular materials：a 3D perspective［J］.Powder technology，2020，366：747-760.

［30］ZHU H，NGUYEN H，NICOT F，et al.On a common critical state in localized and diffuse failure modes［J］.Journal of the Mechanics and Physics of Solids，2016，95：112-131.

（編輯：郭甜甜）

Discrete element analysis on macro-microscopic mechanical properties of rockfill materialYANG Shuhan QI Tianqi LIU Jiaying3，4

（1.Changjiang Survey，Planning，Design and Research Co.，Ltd.，Wuhan 430010，China； 2.National Dam Safety Engineering Technology Research Center，Wuhan 430010，China； 3.Department of Civil Engineering，Hangzhou City University，Hangzhou 310015，China； 4.Zhejiang Engineering Research Center of Intelligent Urban Infrastructure，Hangzhou 310015，China）

Abstract： Mechanical properties of rockfill directly affect the design，construction and operation of the rockfill dam，and have an important influence on the deformation and safety of the dam.In this study，based on the Dashixia concrete face rockfill dam project，ellipsoidal particle samples close to the real particle shape of the rockfill in 3BA zone were established by discrete element method （DEM）.A series of conventional triaxial numerical tests were carried out to study the macro-microscopic mechanical properties of the rockfill in 3BA zone.The macro stress and deformation indexes of rockfill were obtained by the evolution curves of macro stress and strain under different confining pressures.The average contact force was used to divide the contact network into strong and weak subnetworks，and the coordination number distribution and topological evolution of the strong and weak contact networks were analyzed.The results showed that the secant modulus，Poisson′s ratio，peak internal friction angle change linearly with confining pressures，while the dilatancy angle changes nonlinearly with confining pressures.The topological evolution of the contact behaviors in the particle system mainly occurs in the weak contact network，which affects the macroscopic volume change of the rockfill，while the contacts in the strong contact network remains unchanged，providing a stable force transfer way for the rockfill.The research results can be used to guide the selection of macroscopic mechanical indexes under different confining pressures in practical engineering.

Key words： DEM；ellipsoidal particle sample；macroscopic mechanical index；strong and weak contact network；topology；macro-microscopic properties

收稿日期：2023-07-03；接受日期：2023-10-18

基金項目：國家自然科學基金青年基金項目（51909194）；第八屆中國科協青年人才托舉工程全額資助項目（2022QNRC001）；博士后創新崗位基金項目（BSH2021G03）

作者簡介：楊舒涵，女，工程師，博士，主要從事高壩結構設計及數值仿真研究工作。E-mail：yangshuhan@cjwsjy.com.cn

通信作者：漆天奇，男，工程師，博士，主要從事高壩結構設計及數值仿真研究工作。E-mail：qitianqi@cjwsjy.com.cn