基于PFC2d的平行粘結模型宏-細觀參數標定研究

汪 洋， 孟祥林， 高 進， 郭治岳

(1.成都建工路橋建設有限公司，四川成都，610000；2.西南交通大學土木工程學院，四川成都，610031)

在利用離散元模型來進行數值模擬時，顆粒間的每個接觸模型都有多個細觀參數，在選用顆粒間的接觸模型時，往往需要定義較多的細觀參數，為了使顆粒集合體所呈現出的宏觀力學性質與其模擬的土石材料相接近，對細觀參數對宏觀參數的影響進行研究，并利用影響規律進行宏-細觀參數的標定是非常重要的。

國內有許多學者針對離散元軟件的宏-細觀參數影響及參數標定進行了研究。趙國彥等[1]采用理論分析與數值模擬結合的手段，對平行粘結模型中細觀參數對宏觀特性進行了系統性研究，定量給出了部分細觀參數與宏觀參數的關系式；劉暢等[2]通過PFC2d建立巖石材料的單元試驗數值模擬，總結了部分細觀參數與宏觀參數的關系；阿比爾的等[3]采用控制變量法，較為全面的分析了各細觀參數與宏觀參數間的定量關系，并考慮參數相互影響，擬合了經驗公式，對細觀參數選擇進行了優化；郝保欽等[4]采用多因素方差分析法確定了平行黏結模型中細觀參數對宏觀力學參數的影響程度并建立了參數關系式，提出了巖石彈性模量、泊松比單軸抗壓強度等宏觀力學參數的細觀參數確定方法。綜上所述，如何對巖體材料的細觀參數與宏觀力學參數進行標定值得深入探討，為了更好地利用給定巖體宏觀參數匹配出響應細觀參數，本文利用PFC2d離散元軟件，基于單軸壓縮試驗和雙軸壓縮試驗，通過敏感新分析得到了優化的宏-細觀參數標定流程，并進行標定驗證，對準確、高效地進行參數標定具有指導作用。

1 平行黏結模型

利用離散元軟件進行數值模擬時，不需要設置巖體的宏觀力學參數，可針對顆粒選擇合適的接觸模型即反映巖體的宏觀力學特性。本文觸模型選用平行黏結模型，因其具有能夠同時傳遞力和力矩的特點，可以更好地體現完整巖體的力學特性。

基于平行黏結模型，選取emod、krat、pb_emod、pb_krat、pb_ten、pb_coh、pb_fa和fric 8個細觀參數，建立室內單軸壓縮試驗和雙軸壓縮試驗進行數值模擬，考慮巖體的彈性模量、泊松比、內摩擦角和粘聚力，對巖體材料宏-細觀參數的影響規律進行研究并進行參數標定。

2 抗剪強度參數計算

針對多組雙軸試驗數據，抗剪強度參數計算目前最常用的回歸方法為“p-q”法和“σ1-σ3”法[5]。

“p-q”法，是一種間接現行回歸方法，在p-q坐標內進行先行回歸(即縱坐標和橫坐標分別為p和q)，利用最小二乘法進行線性回歸，求得擬合方程的斜率和截距，然后根據式(1)～式(4)可求得抗剪強度參數：

(1)

(2)

f1=tanφ=tan(arcsin(tanα))

(3)

c1=a/cos(arcsin(tanα))

(4)

式中：σ1為峰值應力；σ3為圍壓；a為“p-q”法擬合方程截距；α為“p-q”法擬合方程夾角。

“σ1-σ3”法，是一種直接回歸方法，由于試驗中直接測試得到的為σ1和σ3，因此用該方法回歸更為直接，相關計算公式如式(5)～式(7)所示：

σ1=2ctan(φ/2+π/4)=σ3tan2(φ/2+π/4)

(5)

f2=tanφ=tan(2arctanB0.5-π/2)

(6)

c2=A/(tan(φ/2+π/4))

(7)

式中：A-σ1-σ3法截距，取2ctan(φ/2+π/4)；B-回σ1-σ3法斜率，取tan2(φ/2+π/4)。

利用上述2種方法求得抗剪強度參數后，取平均值作為最終強度參數，公式如式(8)～式(9)所示：

f=f1+f2

(8)

c=c1+c2

(9)

式中：f為最終摩擦系數；c為最終粘聚力。

3 數值模擬

3.1 單軸壓縮試驗

單軸壓縮試驗又被成為無側限壓力試驗，是對圓柱形土樣不加側向壓力，而只在中心軸線上逐步加垂直壓力，直到土樣破壞為止的試驗。模型尺度為50 mm×100 mm，最小粒徑取模型直徑的1/80，粒徑比設置為1.66。為保證單軸壓縮試驗的真實性，數值模擬從取土、制作土樣、土體物理力學性質等方面縮小與實際試驗的差距，對試樣進行預壓模擬土體原來的儲存狀態，考慮巖體埋深區間，按照巖體處于地下50 m進行預壓，預壓值為1 MPa，對土體進行卸載模擬取土后的儲存狀態。單軸試驗數值模型及模擬過程見圖 1。

圖1 單軸壓縮試驗

3.2 三軸壓縮試驗

雙軸壓縮試驗是目前研究土樣抗剪強度較為完善的方法，在進行試驗時，先加液壓即通過把壓力水通入盛土樣的壓力室，使土樣在橫向受到相等的壓力，即為施加“圍壓”，然后維持液壓不變，在試樣軸向施加垂直壓力并不斷增大，直到土樣被壓壞。模型尺度為50 mm×100 mm，最小粒徑取模型直徑的1/80，粒徑比設置為1.66。為保證單軸壓縮試驗的真實性，數值模擬從取土、制作土樣、土體物理力學性質等方面縮小與實際試驗的差距，對試樣進行預壓模擬土體原來的儲存狀態，預壓值設為1 MPa。由于雙軸壓縮試驗至少需要進行3組才可得到有效的土樣抗剪強度指標，選取圍壓為2 MPa、3 MPa和4 MPa進行試驗。雙軸試驗數值模型及模擬過程見圖 2。

圖2 雙軸壓縮試驗

4 敏感性分析

PB設計試驗能通過各因素的2個水平的差異與整體的差異進行比較來確定各影響因素的顯著性，并可以對各因素進行敏感性排序，篩選出具有顯著影響的因素，避免試驗資源的浪費，從而減少不必要的工作量，有助于標定工作的準確性。

本次PB設計試驗中平行粘結模型中各細觀參數設置為兩水平，由于有8組細觀參數，在設置設計矩陣時，選取N=16的設計矩陣開展試驗，試驗設計矩陣及試驗結果見表 1。

表1 設計矩陣及試驗結果

標準化效應圖可反映影響因素的敏感性大小，原理為對各影響因素進行檢驗，并將獲得的效應值的絕對值作為比較指標，按其大小進行排序，并與顯著性水水平所計算得出的顯著性臨界值進行比較，該大小可作為其敏感性排序的依據，效應值越大，越敏感。將上述結果代入計算軟件中，對結果進行因素設計分析，基于16組試驗結果的各宏觀參數的標準化效應圖見圖3。

圖3 各宏觀參數的標準化效應圖

由圖3可知：對宏觀參數彈性模量影響顯著的因素為pb_emod、emod和pb_krat；對宏觀參數泊松比影響顯著的因素為pb_krat、krat和pb_emod；對宏觀參數粘聚力影響顯著的因素為pb_ten和pb_coh；對宏觀參數內摩擦角影響顯著的因素為fric和pb_fa。

5 實例驗證

敏感性分析為宏-細觀參數標定指明了基本流程。通過敏感性分析結果可以在進行參數標定時對細觀參數調整有更為定量的把控。基于“粗調”與“微調”相結合的原則，針對與目標宏觀參數的差距合理調整細觀參數。標定流程見圖4。

圖4 標定流程

本文選取四川成都龍泉山隧道附近圍巖參數，見表 2。考慮到細觀參數較多且試驗過程中的誤差性，標定后的細觀參數反映的宏觀參數在標定值的90%～110%范圍內即可認定為符合標定要求。

表2 現場圍巖參數

按照標定流程，對圍巖參數進行標定，單軸壓縮試驗與雙軸壓縮試驗的模型尺寸與模擬流程見本文第3節，標定結果見表 3。

表3 標定結果

選用標定結果對試驗進行單軸壓縮試驗可得到單軸試驗的應力應變曲線見圖 5。從應力應變曲線中可直觀的得到單軸抗壓強度為1.32 MPa；彈性模量為曲線中直線段終點的應力差與對應的應變差的比值[6]，該組數據得到的應力應變曲線上升段較為平穩，因此用應力峰值與開始加載時的應力差值比應變差值即為彈性模量為2.69 GPa。

圖5 單軸壓縮試驗-應力應變曲線

利用雙軸壓縮試驗可得到不同圍壓下的應力應變曲線見圖 6，試樣在2 MPa、3 MPa和4 MPa下的峰值應力分別為8.91 MPa、12.03 MPa和15.37 MPa。

圖6 雙軸壓縮試驗-應力應變曲線

提取試驗的橫向應變-軸向應變曲線見圖7，可發現不同圍壓下的曲線的變化規律相似，均經歷一段線性增長后加速增長，圍壓越大增長速率越慢，通過線性階段的數據可求得試樣的泊松比為0.285。

圖7 橫向應變-軸向應變曲線

通過“p-q”法和“σ1-σ3法”法計算得到的巖石抗剪強度參數見表4，取2種方法的平均值，求得巖石粘聚力為0.670 MPa，內摩擦角為31.85°。

表4 抗剪強度參數計算結果

6 結論

本文基于PFC2d離散元軟件中的平行黏結模型，選取emod、krat、pb_emod、pb_krat、pb_ten、pb_coh、pb_fa和fric 8個細觀參數，考慮巖體的彈性模量、泊松比、內摩擦角和粘聚力，研究了巖體材料宏-細觀參數的相互影響關系，得到結論：

(1)總結分析了雙軸壓縮試驗中抗剪強度參數的計算方法，確定了以“p-q”法和“σ1-σ3法”法計算得到的抗剪強度參數取平均即為最終抗剪強度參數。

(2)建立考慮原狀土、取土和制作土樣不同儲存條件的單軸壓縮試驗和雙軸壓縮試驗進行數值模擬，最大程度的縮小與實際試驗的差距。

(3)選取N=16的設計矩陣開展PB設計試驗，對8個細觀參數在兩水平條件下進行敏感性分析，得到彈性模量的顯著因素為pb_emod、emod和pb_krat；泊松比為pb_krat、krat和pb_emod；粘聚力為pb_ten和pb_coh；內摩擦角為fric和pb_fa。

(4)確定了細觀參數的標定流程并采用該流程對龍泉山圍巖進行參數標定，基于應力應變曲線、橫向應變-軸向應變曲線完整計算了宏觀參數，驗證了參數標定流程的可行性。