邊坡軟弱夾層宏-細觀參數率定關系研究

閻 飛，梁 兵，孟軻荊

(北京京能地質工程公司，北京 102300)

0 引言

目前，人們對巖土介質的研究正處于由傳統宏觀尺度分別向細觀尺度延伸的過渡時期，顆粒流方法(PFC2D/PFC3D)的興起，為人們提供了一種宏細觀研究的有效手段。PFC是基于Cundall[1]程序發展而來的顆粒流分析程序。Bock等[2]考慮土顆粒形態、結合水、自由水和固結后結構特征構建了新的模型。Utili、Nova[3]用PFC2D系統探討了宏觀強度參數與細觀參數之間的關系。Park、Song[4]用PFC3D實現了節理直剪試驗的三維模擬。張曉平等[5]用PFC2D分析了軟弱夾層幾何參數對試樣力學行為的影響，認為當夾層厚度很小時，實質問題為巖土介質的接觸問題，而夾層厚度較大時，實質問題為夾層的強度問題。周喻等[6]用PFC2D模擬了巖石節理直剪試驗，從宏細觀的角度探討了節理的力學演化特征和破壞機制。顏敬等[7]運用正交設計和人工神經網絡的方法，探索了宏細觀參數之間的互演計算。叢宇等[8]分析了細觀參數取值對彈性模量、泊松比、裂紋分布等宏觀參數和形態的影響。徐金明等[9]結合花崗巖圖像，建立了考慮細觀組分實際分布的顆粒流模型，通過正交設計研究了宏細觀力學性質關系和相應的調參方法。夏磊等[10]將層理面力學參數分為黏結強度、摩擦作用和咬合作用3個方面，研究了不同組合對層狀巖體強度的影響。總結上述成果，可以發現通過PFC2D進行的宏細觀參數研究主要集中在二維接觸剛度模型、滑動模型、平行黏結模型和定義群組所涉及參數的標定上，雖然獲得了一些定性的規律，但尚無法實現應用中的快速調參，針對顆粒間應力鏈和裂紋分布的演化研究明顯滯后。

近幾年來，國內的一些學者也從別的角度對細觀尺度下的巖土介質開展了相關工 作。潘鵬志等[11]運用彈塑性理論、細胞自動機自組織演化理論、統計原理以及巖石力 學理論等，建立了模擬巖石三維破裂過程的張量型細胞自動機模型。羅榮等[12]基于巖石礦物種類及其含量提出了巖石礦物細胞元隨機性參數賦值方法，探討了巖石材料礦物細胞元彈性模量與宏觀模量之間的函數關系。周翠英等[13]采用重整化群的方法概化出紅層軟巖的骨架構件和膠結物構件，推導計算了逾滲閾值。王益棟等[14]給出了顆粒破碎粒徑的分布圖，基于壓縮破碎試驗得到了分形維數。李學豐等[15]把組構張量和加載應力張量合成一個組合張量，定義了一個新的各向異性狀態變量，并建立起各向異性破壞準則。這些成果豐富了人們的研究思路，一定程度上彌補了顆粒流方法的不足，具有很好的借鑒啟示作用。但以DEM為核心，分析顆粒體系的細觀結構及其運動規律，構建巖土介質宏細觀參數關系依然是當前極富學術價值的基礎研究課題，尚處于起步階段。

本文利用顆粒流PFC軟件和均勻設計的思想，通過平行黏結條件下的初步模擬、優化模擬，分自然休止角、無側限抗壓強度、直剪試驗三種模式，對邊坡軟弱夾層的宏—細觀參數率定關系進行了研究，并依據模擬中出現的演化現象，劃分破壞模式，進而得到細觀參數率定建議。

1 均勻設計原理

均勻設計是由方開泰和王元兩位數學家于1978年創立的一種試驗方法[16]，屬于數論中的“偽蒙特卡羅方法”，其思想是要使試驗點在維度空間中“均勻分散”，不同于正交設計的“整齊可比”，每個水平在試驗中都出現且僅出現一次，試驗次數大為減少，結果的處理多采用回歸分析，范圍由大至小，逐步對模型進行修正，提高精度。

1.1 均勻設計表

均勻設計表是試驗設計的依據，具有三條特征：

1)各因素的水平試驗次數僅為1次；

2)任意兩因素的試驗點在平面上，每行每列有且僅有1個；

3)任意兩列之間的均勻性不一定相等，設計時應選擇均勻性較好的列項。類似于正交設計表，每張均勻設計表都有一個代號Ua(bc)，其中U表示均勻設計，a為行數，表示試驗總次數，b表示各因素的水平數，c為列數，代表最多能安排的因素數量。

1.2 均勻設計基本步驟

均勻設計的主要包括選評價指標、定影響因素、定試驗水平、定設計表、回歸分析、優化模型六步。具體在使用時可以采用多步線性回歸分析實現。即先進行初步均勻設計分析，在初步分析的基礎上做進一步的優化均勻設計分析。

由于均勻設計充分利用了試驗點均勻分散的特性，可以大幅減少工作量，在確定自變量與因變量之間響應規律、篩選主控因素、參數率定等方面廣泛適用。本文數值模擬所用到的軟件PFC是利用宏觀試驗結果反演細觀參數組(通常涉及5個參數以上的調參)，從而通過計算再現研究對象的力學響應過程，所以采用均勻設計手段對相關因素進行分析是必要而且合理的。

2 宏-細觀參數演化關系

考慮到常規室內試驗主要包括無側限抗壓和直剪試驗，所以本文的宏細觀參數研究以此兩者為主，同時因為夾層體體現黏性土特性，將黏結模型統一選用為平行黏結，并增加僅存在摩擦參數fri作用時自然休止角試驗的內容。

2.1 自然休止角模擬

對于自然休止角的模擬，采用散落模式，即在一定范圍內隨機生成顆粒，在其左側與下方設定足夠長墻體，施加重力自由下落，顆粒與下墻碰撞反彈后充分迭代平衡，最終形成斜坡角，此角度即為散落模式時的自然休止角，記為φ散，分初步模擬和優化模擬兩步進行方案優化設計。

初步模擬。模型中設置set-disk-on，需要考慮的因素有總顆粒數N、顆粒密度ρ(因為有一定的彈性撞擊)、墻體法向與切向剛度(Kwn、Kws)、顆粒法向與切向剛度(Kbn、Kbs)、摩擦系數fri共計7個。試驗的次數取為因素數的4倍。選用表U28(288)進行各因素水平相等試驗，第8個因素為評判指標φ散。為節省機時，模型中顆粒的粒徑R處于單位cm級別，整個模型的范圍取為長20 m，高12 m的矩形區域，粒子生成的橫向區域為x向0～8 m，y向0.2～12 m，初始設定R為0.1 m，N為1000，顆粒粒徑R隨N的變化而調整，同時應滿足墻體剛度大于顆粒剛度這一要求，故令N范圍為1000～10000，ρ范圍為1000～2700 kg/m3，Kbn與Kbs范圍為107N/m～108N/m，Kwn與Kws范圍為108N/m～109N/m，fri的范圍為0.1～3.0。得出的散落模式(圖1)。

以φ散為因變量，以N、ρ、Kbn、Kbs、Kwn、Kws、fri共計7個參數為自變量，對計算結果進行回歸分析可以得到回歸方程：

φ散=16+1.92e-3N-8.93e-3ρ+1.13Kbn+0.469Kbs+

433e-2Kwn+2.33e-2Kws+7.43fri

(1)

按照方開泰等人的計算方法[16]，各因素貢獻率依次為24.8%、19.2%、8.55%、1.48%、1.27%、0.368%、38.7%。結果說明所試驗的7個因素中，起主要因素的為摩擦參數fri、顆粒數N(粒徑R)和密度ρ，顆粒體的法向剛度Kbn需要進一步分析，其余3個因素為非相關項。

圖1 散落模式示意圖Fig.1 Schematic diagram of scattering mode

優化模擬。由初步的分析，進一步優化均勻設計方案，設置球的兩向剛度相同，墻體剛度取為其10倍，以分析顆粒法向剛度Kbn對φ散的影響。試驗次數取為因素數4的5倍即20水平，選用表U20(204)進行各因素水平相等均勻試驗。粒子生成區間為x向0～8 m，y向0.3～15.3 m，下落距離統一為0.3 m，顆粒數N不再作為直接因素，而以粒徑R代替，模擬粒組為粗砂0.5～2.0 mm，計算時為節約機時，統一將R擴大100倍，即為0.05～0.20 m，N由等面積換算得到(初始面積為0.1 m，N為1000時的面積)，ρ范圍為1000～2700 kg/m3，Kbn范圍為107N/m～108N/m，fri范圍為0.1～1.5，評判指標取φ散。

以φ散為因變量，以R、ρ、Kbn、fri共計4個參數為自變量，對計算結果進行回歸分析可以得到回歸方程：

φ散=14.6-0.147R-8.16e-4ρ+3.36e-3Kbn+

16.4fri

(2)

各因素貢獻率依次為0.707%、0.225%、0.0132%、61.8%，從中可以看出起主要因素的只有參數fri。

令生成區間為x向0～8 m，y向0.3 ～15.3 m，顆粒數N為1500，顆粒半徑為0.1 m，球體兩向剛度設置為106N/m，墻體兩向剛度設置為107N/m，計算fri在0.1～1.5之間變化時，fri與自然休止角φ散的對應關系見表1。可以看到，φ當在0°～42°范圍內時，φ與φ散只差Δφ基本處于12°之內，建議回歸方程為

φ散=5.38+28.618fri

(3)

2.2 無側限抗壓試驗模擬

對于無側限抗壓試驗的模擬，參考叢宇等人的工作[8]，設定數值模型長100 mm，寬50 mm，初始顆粒半徑最小為0.25 mm，最大為1.0 mm，采用半徑擴大法生成指定孔隙率n的試樣。試驗前先利用四面墻體令顆粒充分平衡，再刪去兩側約束，隨后給上部墻體一個向下的速度進行加壓，獲得評判指標彈性模量E和極限抗壓強度Rc。

表1 自然休止角與fri關系Table 1 The relation between natural angle of repose and fri

初步模擬。由于黏結模型為平行黏結，故初步選定平行黏結切向剛度Kps、平行黏結剛度比Kpn/Kps、平行黏結法向強度Ppn、平行黏結切向強度Pps、平行半徑乘子λ和參數fri六項為作用因素，研究平行黏結模型對E與Rc的影響。根據文獻中正交試驗所得結果，令試樣n=0.20，顆粒法向剛度為23 GPa，切向剛度為8.8 GPa，密度ρ=2000 kg/m3，墻體兩向剛度為200 GPa，同時Kps為10 GPa～40 GPa，Kpn/Kps值為0.2～5.0，Ppn為20 MPa～100 MPa，Pps為20 MPa～100 MPa，λ為0.1～1.0，fri為0.3～1.5。試驗次數取為因素數6的4倍即24水平，選用表U24(246)進行各因素水平相等均勻試驗。

模擬過程中，根據《土工試驗方法標準》(GB/T 50123-2019)中相關條文，對有峰值試樣取峰值點作為無側限抗壓強度Rc，而對無峰值試樣取應變為0.15處的應力值作為Rc，彈性模量E通過近直線段的斜率計算得到。首先以E為評判指標進行回歸分析，方程為

E=216+10.9Kps+61.1Kpn/Kps-2.83Ppn-

0.788Pps-56.6λ-77.5fri

(4)

各因素貢獻率依次為23.3%、15.3%、11.2%、0.691%、0.452%、1.54%。結果表明，對于E而言Kps、Kpn/Kps值和Ppn顯然為主控因素，其余三項居于從屬地位，可不計。

然后以Rc為評判指標進行回歸分析，可得線性方程如下：

Rc=-2.97+0.949Kps+2.31Kpn/Kps+0.126Ppn-

0.18Pps+19.5λ+3.05fri

(5)

各因素貢獻率依次為16.6%、7.60%、7.69%、12.5%、18.7%、0.828%。由結果可知，對Rc而言，Kps、Pps、λ為主控因素，Kpn/Kps值、Ppn為次要因素，fri作用可忽略。

優化模擬。根據初步模擬的結果以及在自然休止角試驗中獲得的認識，分別以E和Rc為評判指標進行均勻設計方案優化。

指標E除平行黏結中的Kps、Kpn/Kps值、Ppn三個主控因素外，考慮顆粒及試樣本身性質影響，將顆粒切向剛度Kbs、顆粒兩向剛度比Kbn/Kbs和孔隙率n列為考察對象。墻體為體現加載板的剛性，依然取200 GPa，其余非分析項參數的取值分別為ρ=2000 kg/m3、Pps=80 MPa、λ=0.5、fri=0.8。而分析項參數Kbs范圍為2 GPa～20 GPa、Kbn/Kbs值范圍為0.2～5.0、n范圍為0.05～0.50、Kps范圍為2 GPa～20 GPa、Kpn/Kps值范圍為0.2～5.0、Ppn范圍為10 MPa～300 MPa，試驗次數取為因素數6的4倍即24水平，選用表U24(246)進行各因素水平相等均勻試驗。

對試驗結果進行回歸后可得

E=-339+0.366Kbs+2.49Kbn/Kbs+363n+

8.55Kps+8.55Kpn/Kps-0.0260Ppn

(6)

各因素貢獻率依次為0.0953%、0.256%、58.6%、41.5%、6.69%、0.102%。上述結果說明，主控因素為試樣孔隙率n和平行黏結切向剛度Kps，次要因素為Kpn/Kps值，剩余3項為可忽略項。

指標Rc類似處理，這里將平行黏結中的λ、Kps、Pps和試樣顆粒的Kbs、Kbn/Kbs值、n共6項選為考察因素，取值如下：λ范圍為0.1～1、Kps范圍為2 GPa～20 GPa、Pps范圍為10 MPa～300 MPa、Kbs范圍為2 GPa～20 GPa、Kbn/Kbs值范圍為0.2～5.0，n范圍為0.05～0.50。其余非分析項參數取值依次為墻體兩向剛度 200 GPa、ρ=2000 kg/m3、Ppn=80 GPa、Kpn/Kps=1.0、fri=0.8。試驗次數取為因素數6的4倍即24水平，選用表U24(246)進行各因素水平相等均勻試驗。回歸可得

Rc=-32-0.00723Kbs+0.912Kbn/Kbs+36.8n+

11.6λ+0.372Kps+0.00952Pps

(7)

各因素貢獻率依次為0.00419%、3.87%、67.9%、21.7%、8.87%、1.54%。結果表明，主控因素為試樣孔隙率n和平行黏結半徑乘子λ，次要因素為Kps、Kbn/Kbs值、Pps，Kbs為可忽略項。

在優化模擬過程中，試樣的破壞形態和應力應變曲線見圖2。

圖2 試樣破壞模式及其應力應變曲線Fig.2 Failure mode and stress-strain curve(a) 擠出式破壞 (b) 擠出應力應變 (c) 鼓脹式破壞 (d) 鼓脹應力應變 (e) 剪切式破壞 (f) 剪切應力應變

破壞形態分為三種，分別是擠出式、鼓脹式和剪切式，其中對于以E為指標的24個模擬試樣，擠出式占37.5%(9個)，鼓脹式占58.3%(14個)，剪切式占4.2%(1個)，而對于以Rc為指標的24次試驗，擠出式占33.3%(8個)，鼓脹式占66.7%(9個)，剪切式沒有出現。結合均勻設計的結果，可以得到以下幾點認識：① 試樣的孔隙率n對無側限抗壓試驗起重要作用，當n處于0.3～0.5之間時，試樣破壞以擠出式為主。這是因為生成的顆粒之間空隙較大，自由生成時會出現部分顆粒脫空懸浮，從而使得試樣整體性較差，當加載進行時，受壓部分無法形成力鏈的有效傳導，從而出現擠出效果，此時其應力應變曲線呈現寬幅穩定發展的態勢，而當荷載作用范圍傳遞至底層剛性板時，應力逐步產生集中，從而令曲線轉而變陡。整體而言，擠出式破壞對于模擬高孔隙比軟土試樣較為合適。② 當n<0.3時，顆粒之間基本上都可以形成有效接觸，平行黏結在荷載的作用下充分發揮作用，變形呈現出整體性，比較初步模擬和優化模擬結果可知，試樣無側限抗壓強度試驗的破壞模式主要為鼓脹式和剪切式兩種，反映了中硬土、硬土和巖質試樣的變形特征，且由其各自所占比例的變化推測，各參數之間存在某一界限式組合決定了試樣的最終破壞形態。③ 觀察圖2e，試樣的剪切破壞較為明顯，殘存體呈啞鈴狀，兩端大中間窄，通過觀察變形歷史發現，試樣起初為鼓脹式變形，到達一定程度之后，中心處顆粒突然張裂，隨即形成一組近共軛的斜向裂隙，貫穿后試樣逐步解體，并在一定程度上形成了新的應力傳導途徑進而再次發生剪切破壞，與其對應的是應力應變曲線上峰后發生的二次應力攀升。④ 初步模擬和優化模擬所篩選出主控因素的不同，部分源于破壞模式的差異，由于試驗本身不存在拉力源，故其本質上均屬于壓剪，剪切面是否貫通直接影響了平行黏結切向參數的作用比例。

上述結果也進一步說明，利用顆粒流對巖土介質力學響應的試驗模擬應首先明確其破壞模式，并以此為前提進行均勻方案設計和正交試驗方案設計，從而做到有的放矢。

對所有模式而言，影響E的主控因素為孔隙率n和平行黏結Kps、Kbn/Kbs值，影響Rc的主控因素為孔隙率n和平行黏結半徑乘子λ，其余參數項則需要結合破壞模式分析是否列入主控因素范疇。

2.3 直剪試驗模擬

對于直剪試驗的模擬，在周喻[6]等人工作的基礎上進行，設定試樣模型長100 mm，高40 mm，初始顆粒半徑最小為0.3 mm，最大為0.525 mm，采用半徑擴大法生成指定孔隙率n的試樣。試驗前先利用四面墻體令顆粒充分平衡，并編制伺服程序對試樣上下墻體進行加壓控制，保證壓力達到模擬要求后，對試樣進行橫向剪切直至破壞，該過程可以獲得典型的剪切曲線。重復該過程，即可得到不同壓力下試樣的破壞點，從而統計計算出試樣的黏聚力c和內摩擦角φ。

初步模擬。由于黏結模型為平行黏結，故初步選定平行黏結切向剛度Kps、平行黏結剛度比Kpn/Kps、平行黏結法向強度Ppn，平行黏結切向強度Pps、平行半徑乘子λ和參數fri六項為作用因素，分別模擬壓力為0.25 MPa、0.50 MPa、0.75 MPa、1.00 MPa、1.25 MPa五種情況時試樣的直剪試驗，從而得到相應的c和φ值。

根據文獻中提供的參數，令試樣n=0.20，顆粒法向剛度15 GPa，切向剛度2.5 GPa，密度2750 kg/m3，墻體兩向剛度為50 GPa，同時Kps為0.8～8.0 GPa，Kpn/Kps值為0.2～5.0，Ppn為1 MPa～40 MPa，Pps為1 MPa～40 MPa，λ為0.1～1.0，fri為0.3～1.5。試驗次數取為因素數6的4倍即24水平，選用表U24(246)進行各因素水平相等均勻試驗。

在24組試樣模擬過程中，剪切曲線總體表現為直線型、硬化型和峰值型三類(圖3)，其中硬化型又包含有硬化緩變型與硬化直線型兩亞類。根據《土工試驗方法標準》(GB/T 50123-2019)中相關條文，對峰值型取峰值剪應力作為抗剪強度，對其余試樣取剪切位移4 mm所對應的剪應力作為抗剪強度，進而得到各組試樣的c和φ值。

首先以c為評判指標進行回歸，方程為

c=-411+59.3Kps+66.6Kpn/Kps-0.847Ppn+

0.471Pps+574λ-61.3fri

(8)

各因素貢獻率依次為26.3%、12.1%、0.158%、0.039%、30.8%、0.638%。結果表明，對于c而言，λ、Kps、Kpn/Kps值是主控因素，其余三項居于從屬地位，可忽略。

以φ為評判指標進行回歸分析，可得線性回歸方程：

φ=9.51+0.907Kps+0.308Kpn/Kps+0.0826Ppn+

0.0607Pps+9.87λ+2.87fri

(9)

各因素貢獻率依次39.2%、1.65%、9.52%、4.12%、57.3%、8.91%。即對于φ而言，λ、Kps是主控因素，Ppn、fri是次要因素，而Pps、Kpn/Kps值作用可忽略。

24組試樣中，上述四種類型所占比例分別為直線型∶硬化緩變型∶硬化直線型∶峰值型為9∶10∶2∶3，前兩類占絕對多數。其中直線型主要表現為剪切時曲線寬幅振蕩，低圍壓(0.25 MPa)條件下顆粒間力鏈發育不完全，而高圍壓下顆粒間聯系則較好，同時伴隨剪切的進程，試樣有壓密現象，曲線前部較陡后部呈一較緩直線形態。由動態過程記錄可知，剪切的顯著影響范圍為2～3倍墻體長度，該現象普遍存在于所有模擬試樣當中，剪切面及顆粒破壞形態由于細觀參數的差異而不同，對于直線型表現為前部拖拽和后部下傾。硬化緩變型主要表現為，剪切曲線斜率逐漸變緩，振幅趨小，顆粒間力鏈發展較完全，剪切面也較為平順。與之相比，硬化直線型由于顆粒間作用更強，剪切曲線振幅更小，趨于直線，力鏈發展更為完全，但剪切面不明顯，前后部呈現出拖拽、翻轉現象。峰值型的剪切曲線與巖土體的力學響應行為最為相似，表現為峰值后的軟化現象，力鏈發展完全且剪切面較為平順。

圖3 試樣剪切曲線類型

顯然峰值型是下一步優化模擬的重點，故在選取λ和Kps為分析對象的基礎上，加入球體細觀參數的影響，分析其主控因素。同時根據已有試樣的破壞形態，對模型進行縮減長度和增加墻體剛度的調整。

優化模擬。優化模擬中選定試樣顆粒的Kbs、Kbn/Kbs值、n和平行黏結中的Kps、λ五項因素作為考察對象，分別模擬壓力為0.25 MPa、0.50 MPa、0.75 MPa、1.00 MPa、1.25 MPa五種情況時試樣的直剪試驗，從而得到相應的c和φ值。根據初步模擬結果，令Kbs為2 GPa～20 GPa，Kbn/Kbs值為0.2～5.0，n為0.05～0.35，Kps為0.8 GPa～24 GPa，λ為0.1～1.0，同時Ppn為6.09 MPa，Pps為28.13 MPa，fri為1.187，Kpn/Kps值固定為1.0，密度ρ為2750 kg/m3，調整模型長度為75 mm，墻體剛度為200 GPa。試驗次數取為因素數5的4倍即20水平，選用表U20(205)進行各因素水平相等均勻試驗。

通過優化后20組試樣的剪切計算，所得曲線總體隸屬于峰值型和硬化緩變型兩類，出現比例為10∶10，并且計算時步內均具有峰值或極大值，故取其對應點剪應力作為抗剪強度，進而由5個圍壓下的數據點，繪圖得到各組試樣的c和φ值。

對c進行回歸方程為

c=-1060-16Kbs+114Kbn/Kbs+778n+

57.3Kps+1910λ

(10)

各因素貢獻率依次為2.73%、9.80%、1.44%、46.5%、75.7%。結果表明，對于c而言，λ和Kps是主控因素，其余三項居于從屬地位，可忽略。然后以φ為評判指標進行回歸分析，方程

φ=31.8+0.273Kbs-2.87Kbn/Kbs+14n-

0.241Kps-4.99λ

(11)

各因素貢獻率為9.71%、76.6%、5.68%、10%、6.35%，即對于φ而言，Kbn/Kbs值是主控因素，而Kbs、n、Kps和λ所占影響權重基本處于同一量級，起到局部調節作用。

模擬過程中出現了一些細觀現象。首先是剪切破裂面的形態，如圖4所示有三種情況。第一種為斜剪，存在于個別試樣低圍壓0.25 MPa時，由于顆粒間隨之擠壓錯動，當上表面集聚的應力明顯大于初始應力時，試樣局部發生抬動，剪切面隨之向上發展，直至剪破試樣，其應力位移曲線表現為近線性增大。第二種表現為中后部拱起，這種形態的形成主要是因為顆粒間本身的作用較強，而力鏈的傳導效應有一定的過程和范圍，中后部段顆粒體容易形成較難剪破的團粒，如同凸起一樣阻礙了變形的發展，使得兩側無法很好協調，從而局部孔隙度增大，出現擴容，當剪切面爬升越過該部分后，試樣完全破壞。這種現象在低圍壓向高圍壓轉變的過程中較為常見，應力位移曲線表現為硬化緩變型。第三種形態，剪切面較為平順，基本沿水平向發生破裂，圖中由于上部顆粒的拖拽以及局部錯動，形成較大的空隙，并且在模擬過程中，已破壞顆粒會隨荷載板的運動重新發生黏合，力鏈再次形成發揮傳導作用，也構成了部分的殘余強度，試樣的應力位移曲線是典型的峰值型。

其次是力鏈的發展演化，基本影響范圍仍然是2～3倍墻體長度，發育密度低壓力下稀疏，高壓力下密集，試樣隨剪切的進行，左下角和右上角存在一定的應力遮蔽現象(圖5)。此外，試樣孔隙度的不同會對模擬所得強度的穩定性造成較大影響，在圖6中列舉了兩種典型情況。其中密實型試樣孔隙率低，由計算的20組結果統計，n<0.24，顆粒間相互接觸較為容易，圍壓加載時細觀結構相對保持穩定，庫侖強度線符合線性規律。而松散型試樣孔隙率較大，本次計算中主要是n處于0.24～0.35的剪切組，因為顆粒間聯系并不充分，圍壓加載時試樣便發生了一定的壓密，從而改變了初始的細觀結構，這就不難解釋不同圍壓下抗剪強度點的離散特征，這樣的情況有6組，但總體而言依然遵循壓力越高，線性關系越顯著的規律。上述現象也說明巖土體的結構性，包括宏觀與細觀兩個方面，都會對最終的抗剪強度及其演化關系造成較大影響。

基于上述初步模擬和優化模擬的結果，可以得到直剪試驗強度參數c和φ值與顆粒細觀參數的基本關系如下：

1)對于黏聚力c而言，起主要作用的細觀參數為平行半徑乘子λ、平行黏結切向剛度Kps和平行黏結剛度比Kpn/Kps，影響權重為λ>Kps>Kpn/Kps，且均與c成正相關性，其他參數影響可以忽略不計。

2)對于摩擦角φ而言，顆粒本身的剛度比Kbn/Kbs

圖4 剪切破壞面形態Fig.4 Shape of the shear failure surface(a)斜剪 (b) 中后部拱起 (c) 平順

圖5 剪應力遮蔽Fig.5 Shear stress shielding(a)初始段 (b) 剪切段 (c) 破壞段

圖6 孔隙率對結果影響

是主控因素，其與φ成負相關性，其他參數λ、Kps、Kbs、n屬于次要因素，且各自影響權重相當。

3)孔隙率n的取值對于模擬過程中試樣的力學響應行為影響較大，對于此次的模擬計算組，當n≤0.24時，顆粒間相互接觸較為容易，圍壓加載時細觀結構相對保持穩定，庫侖強度線符合線性規律，而當n>0.24時，顆粒間聯系并不充分，圍壓加載時試樣便發生了一定的壓密，從而改變了初始的細觀結構，令抗剪強度統計點離散性增大。這一現象也從細觀角度說明結構性是控制巖土體力學性能的重要因素。

4)細觀參數平行黏結法向強度Ppn、平行黏結切向強度Pps和fri對于試樣剪切行為的影響可以忽略，這三者都是決定顆粒與顆粒之間局部是否破壞和相互之間摩擦關系的，但對于宏觀上力鏈的發展方向、發育程度以及試樣的破壞模式等，并非主控因素。

5)24組試樣的剪切結果表明，低圍壓向高圍壓轉變的過程中，試樣的破壞模式可能發生變化，具體表現斜剪、中后部拱起、平順三種，并且高圍壓條件下更有助于提升數據點的穩定性，獲得理想的剪切模擬效果。

6)由于剪切模擬過程中力鏈的動態演化范圍有限，約為剪切板長度的2～3倍，故試樣的長度不宜過長，留足剪切位移長度即可，否則可能會出現不合理的結果和現象，這也是尺寸效應的一種表現。

3 宏細觀參數率定建議

通過對自然休止角、無側限抗壓強度和直剪試驗的均勻設計，可以發現不同細觀參數的組合對于模擬結果有著重要影響，主要體現在局部顆粒的響應行為、整體試樣的破壞模式、強度參數的離散性等方面。結合已進行的工作和規律認識，對細觀參數的率定進行討論和建議。

1)fri反映的是顆粒之間的摩擦作用，但對于宏觀上庫侖定律中摩擦角的影響很小，可以忽略不計，該參數主要用于加速模型生成時不平衡力的迭代平衡以及調節顆粒細觀上的響應行為。如對非黏性土體的模擬，當fri<0.1時，試樣會出現局部擠出效應，這一特點可以用于模擬軟弱夾層的擠出-滑移型破壞。

2)無側限抗壓強度試驗模擬中，可以根據E和Rc的實測值對細觀參數進行率定。對于E而言，主控因素為孔隙率n和平行黏結Kps，次要因素為Kpn/Kps值，三者均與E成正相關。對于無側限抗壓強度Rc而言，主控因素為孔隙率n和λ，次要因素為Kps、Kbn/Kbs值和Pps。

3)直剪試驗模擬中，可以根據c和φ的實測值對細觀參數進行率定。對于黏聚力c而言，平行半徑乘子λ、平行黏結切向剛度Kps和平行黏結剛度比Kpn/Kps為主控因素，三者均與c成正相關。對于摩擦角φ而言，顆粒本身的Kbn/Kbs值為主控因素，與φ成負相關，其他參數λ、Kps、Kbs、n屬于次要因素，且各自影響權重相當。

4)孔隙率n對于模擬結果的影響較大，主要是n的改變使得試樣的細觀結構發生了變化。體現在無側限抗壓強度試驗模擬中，表現為當n處于0.3～0.5時，試樣的破壞模式為擠出式，當n<0.3時，試樣的破壞模式為鼓張式、剪切式；體現在直剪試驗模擬中，表現為當n≤0.24時，顆粒間相互接觸較為容易，圍壓加載時細觀結構相對保持穩定，庫侖強度線符合線性規律，當n>0.24時，顆粒間聯系并不充分，圍壓加載時試樣便發生了一定的壓密，從而改變了初始的細觀結構，令抗剪強度統計點離散性增大。上述兩者均從細觀角度說明結構性是控制巖土體力學性能的重要因素。

5)試樣模型的構建中，縱向長度宜大于施壓板長度的2～3倍，孔隙率n宜設置小一些，并且圍壓宜設置大一些，提升模擬結果的穩定性。

仔細分析上述幾條建議，特別是細觀參數與宏觀參數之間的關系，不難發現他們具有很好的對應性。如c、Rc等強度參數的主控因素中，λ和n都起到了很大的作用，前者代表了顆粒間平行黏結的作用范圍，而n則是通過細觀結構影響破壞模式，最終改變試樣的強度。同樣在變形參數E的主控因素中，可以看到影響因素包括Kps和Kpn/Kps，其說明連接的側向變形抵抗能力決定著試樣縱向變形的響應情況。φ主要受到顆粒剛度比Kbn/Kbs的影響，fri作用可以忽略，其他因素對φ的影響權重基本處于同一量級。總結起來就是，n決定了初始試樣的細觀結構差異，λ決定了顆粒間黏結作用的范圍，剛度和剛度比(包括球體和平行黏結)決定了力鏈的發育方向。

上述的建議和討論可以為PFC的有效模擬與細觀參數快速率定提供思路。即可以根據不同的試驗項目，結合實測試驗數據分別對主控細觀參數進行分析，確定范圍。例如直剪試驗，可以先由φ獲得Kbn/Kbs值的初步取值，再根據c確定λ、Kps和Kpn/Kps等參數也可以為無側限抗壓強度試驗的細觀參數確定提供借鑒；反之，對于同一種巖土體，無側限抗壓強度試驗的參數組合也可以與直剪試驗組對比校核。這樣所得的細觀參數會與宏觀實際具有更好的一致性。在率定的過程中，應注意模型試樣的長度設定，并且宜采用高圍壓條件，提升模擬結果的穩定性。

4 結論

利用顆粒流PFC軟件和均勻設計的思想，通過平行黏結條件下的初步模擬、優化模擬，分自然休止角、無側限抗壓強度、直剪試驗三種模式對巖土體的宏細觀參數關系進行了研究，進而得到細觀參數率定建議，所得主要結論如下：

1)對于非黏性顆粒，自然休止角與細觀參數fri關系密切。

2)無側限抗壓強度模擬中，主要破壞呈現出擠出式、鼓脹式和剪切式三種類型。宏觀參數彈性模量主要受控于孔隙率和平行黏結切向剛度。無側限抗壓強度主要受控于孔隙率和平行半徑乘子。

3)直剪試驗模擬中，破壞面有斜剪、中后部拱起和平順三種典型形態。對于黏聚力而言，起主要作用的細觀參數為平行半徑乘子、平行黏結切向剛度和平行黏結剛度比。摩擦角主要受控于顆粒本身的剛度比。

4)孔隙率n對試樣的力學響應有較大影響。試樣模型的構建中，縱向長度宜大于施壓板長度的2～3倍，孔隙率n宜設置小一些，并且圍壓宜設置大一些，提升模擬結果的穩定性。