不同邊界對花崗巖三軸試驗影響的三維離散元數(shù)值研究

王治林，鄭明明，夏敏，熊亮，吳祖銳，王凱

（成都理工大學地質(zhì)災害防治與地質(zhì)環(huán)境保護國家重點實驗室，四川 成都 610059）

0 引言

花崗巖是由不同種類較大晶體顆粒與非晶體顆粒等構成的非連續(xù)、非均質(zhì)且力學行為較為復雜的巖石，且在諸多掘進工程中極為常見，楊達等［1］認為使用氣動潛孔錘并且當球齒以30 J沖擊功鉆進花崗巖時，可有效降低能量損耗與提高碎巖效率。湯鳳林等［2］研究了鉆進速度與規(guī)程參數(shù)的優(yōu)化關系，提出了優(yōu)化正常鉆進的措施。張程等［3］通過室內(nèi)實驗研究了超聲波碎巖技術下花崗巖徑向位移與內(nèi)部損傷的關系，并認為超聲波碎巖技術能提高在花崗巖等硬巖中的鉆進效率，以上學者要么從鉆進設備要么從工藝方法等方面研究如何提高鉆進效率，而未考慮到花崗巖本身的力學特性，目前許多學者利用數(shù)值方法研究巖石的力學特性，而在數(shù)值方法中，離散元法對于探究巖石在變形破壞的細觀機理方面具有巨大的優(yōu)勢［4-6］。

許多學者［7-9］使用顆粒流軟件（PFC3D）分析巖石在不同圍壓下、加載速率下的力學特性時其邊界均使用剛性邊界，而常規(guī)室內(nèi)三軸試驗中，試樣基本都采用橡膠膜或乳膠膜密封，剛性邊界不符合實際邊界情況，因此，金磊等［10］根據(jù)細分概念，提出三維組合墻法模擬三軸邊界，其在加載過程中體現(xiàn)出與室內(nèi)試樣相似的橫向膨脹現(xiàn)象，能良好地模擬試樣的變形特征。Xu等［11］在雙軸實驗時使用一定厚度的光滑顆粒集合體來模擬柔性薄膜，發(fā)現(xiàn)邊界條件一定程度上會影響剪切帶的形成和形態(tài)。Cheung等［12］，Thmos等［13］在雙軸實驗中對比了柔性邊界和剛性邊界對試驗模擬結果的影響，發(fā)現(xiàn)剛性邊界下試樣變形與室內(nèi)試樣變形形態(tài)不符，而柔性邊界下試樣能夠體現(xiàn)出與室內(nèi)試樣相同的橫向膨脹現(xiàn)象。Qu等［14］在PFC中提出用膜顆粒微觀參數(shù)表示柔性膜邊界物理屬性的方法，并與大量三軸實驗結果對比分析，指出使用柔性膜邊界是當下較為適合的離散元邊界模擬方法。但對于室內(nèi)試驗所用的橡膠膜或乳膠膜，無論是柔性邊界還是三維組合墻在力學本質(zhì)上都難以模擬出其彈塑性特征，而對于彈塑性物體的模擬，有限元軟件有著極其顯著的優(yōu)勢。

對此，基于最新的有限元-離散元（FDM-EDM）耦合技術，使用PFC分別生成“剛性”、“柔性膜”、“Shell”單元三種邊界下的花崗巖三軸試驗模型，進而分析“剛性”、“柔性膜”、“Shell”單元3種邊界對花崗巖常規(guī)三軸模擬試驗的影響。

1 花崗巖模型

1.1 花崗巖細觀結構

花崗巖中的礦物顆粒是影響花崗巖力學性質(zhì)的重要因素之一，也是巖石非均質(zhì)性的重要體現(xiàn)，花崗巖中礦物成分主要為長石、石英、黑云母以及其它礦物等，根據(jù)張濤等［15］的研究，本文花崗巖內(nèi)部礦物成分含量為：長石（57.6%）、石英（30.4%）、黑云母（10.1%）、其它礦物（2.0%）。本文中圓柱形巖樣半徑為2.1 cm，高度9 cm，巖樣中顆粒數(shù)量為14799，粒徑范圍為0.028～0.056 cm，粒徑分布服從均勻分布。

1.2 花崗巖模型的生成

顆粒流程序中不能直接生成三維等效晶質(zhì)模型（Grain-Based-Model），須借助其軟件中自帶的“塊體（Rblock）”單元或者“幾何體（geometry）”單元來間接生成三維等效晶質(zhì)模型。張濤等［15］是通過3個不同階段來生成三維等效晶質(zhì)模型：（1）構建試樣初始顆粒模型；（2）礦物晶體分組，在初始顆粒模型外部構建一個外接“幾何體”單元，將外接幾何體單元轉化為“塊體”單元，根據(jù)礦物成分對“塊體”單元進行分組，再將分組后的“塊體”單元集合再次轉化“幾何體”單元集合；（3）利用“幾何體”單元進行填充，最終生成三維等效晶質(zhì)模型。

上述過程最終實現(xiàn)了花崗巖三維等效晶質(zhì)模型的建立，但其生成步驟十分繁瑣，且“塊體”單元與“幾何體”單元之間需要多次轉化，降低了計算效率，為此本文在其基礎上通過直接使用“Rblock”對初始顆粒模型進行劃分，從而簡化了生成步驟，提高了計算效率。花崗巖模型生成過程中的3個階段如圖1所示。

階段一：確定試樣尺寸，構建初始顆粒模型，如圖1（a）所示。

階段二：“塊體”單元分組。在初始顆粒模型外部構建一個外接“幾何體單元”，將外接“幾何體”單元再次轉化為“塊體”單元，對“塊體”單元進行分組，如圖1（b）、（c）所示。

階段三：生成三維等效晶質(zhì)模型。通過判斷“塊體”單元與初始模型顆粒位置的關系對初始模型顆粒進行分組，從而生成三維等效晶質(zhì)模型，如圖1（d）、（e）所示。

圖1 花崗巖模型生成示意Fig.1 Schematic diagram of granite model generation

1.3 模型微觀參數(shù)確定

顆粒流程序中自帶多種接觸模型，其中平行粘結模型（Parallel Bond Model）適合用來模擬巖石顆粒之間的接觸［16-19］，因此本文中巖石顆粒之間的接觸均采用平行黏結模型。此外本文將巖石內(nèi)部顆粒之間的接觸稱為晶體內(nèi)接觸和晶體間接觸。根據(jù)Ding等［20］以往的花崗巖常規(guī)三軸離散元模擬研究以及目前效果較好的“試錯法”，本文中微觀接觸參數(shù)如表1所示。

表1 平行粘結模型微觀參數(shù)Table 1 Micro-parameters of the parallel bonding model

2 三軸模擬試驗建立

邊界的建立和伺服是三軸模擬試驗的關鍵步驟，前面已介紹過如何建立花崗巖模型的方法，因此這里主要介紹本文中花崗巖三軸壓縮試驗的3種邊界。

2.1 “剛性”邊界

PFC可以直接生成常規(guī)三軸模擬試驗中的剛性邊界，剛性邊界屬于剛體，受力后不會產(chǎn)生任何變形，但通過剛性伺服函數(shù)可以獲得徑向方向上的一個平移自由度，細觀上表現(xiàn)為剛性邊界頂點的徑向位移，宏觀上表現(xiàn)為剛性邊界的整體放大與縮?。?1］，由于毛海濤等［21］對剛性伺服機制的深入解釋，本文中便不再詳細解釋剛性伺服機制，“剛性”邊界如圖2（a）所示。

2.2 “柔性膜”邊界

蔣成龍等［22］詳細介紹了柔性膜的生成方法和伺服理論。本文利用其方法生成常規(guī)三軸實驗所用的柔性膜邊界，基于柔性膜伺服理論，使柔性膜邊界下的花崗巖模型達到目標圍壓。柔性膜顆粒在空間上具有6個自由度，可以任意變形，但本質(zhì)上屬于離散元邊界，Qu等［14］經(jīng)過理論推導、試驗驗證找到了柔性膜顆粒間的接觸參數(shù)與柔性膜宏觀參數(shù)（彈性模量、泊松比）的關系，基此，本文使用Qu等［14］發(fā)現(xiàn)的關系通過柔性膜宏觀參數(shù)確定柔性膜顆粒的接觸參數(shù)，本文中柔性膜邊界的彈性模量和泊松比分別為7.8 MPa、0.46。“柔性膜”邊界如圖2（b）所示。

圖2 3種不同邊界條件示意圖Fig.2 Schematic diagram of three different boundary conditions

2.3 “Shell”單元邊界

室內(nèi)試驗所用的橡皮膜在力學本質(zhì)上屬于彈塑性物體，而在模擬彈塑性物體時，有限元軟件比離散元有巨大優(yōu)勢［23］，有限元法中的“Shell”單元由帶有3個節(jié)點的三角面組成，其中每個節(jié)點都具有6個自由度（3個平移自由度、3個旋轉自由度），如圖3所示：其變形分為膜變形以及梁變形，對于各向同性的“Shell”單元，其力學控制方程符合式（1）與式（2），其中E和v為“Shell”單元的彈性模量和泊松比，本文中E=7.8 MPa，v=0.46。同時其顯著優(yōu)點是：“Shell”單元邊界的圍壓可以通過Structure Shell Apply命令精準施加且圍壓在加載過程中保持恒定不變，因此引入FLAC中的“Shell”單元作為三軸模擬試驗的邊界。本文中“Shell”單元邊界由2304個三角形“Shell”單元均勻組成，如圖2（c）所示。

圖3 Shell單元示意圖Fig.3 Example of Shell element

2.4 花崗巖三軸加載

為在離散元中研究邊界對花崗巖三軸壓縮試驗的影響，設計9種巖樣在3種不同邊界下的模擬實驗，共計27組模擬試驗，結果見表2。

表2 模擬巖樣標號及試驗參數(shù)Table 2 The number of samples and experimental parameters

3 不同邊界下模擬結果的對比分析

3.1 3種邊界下的巖樣變形差異

為對比不同情況下的巖樣變形破壞情況，同時鑒于文章篇幅，本文僅列出了X61巖樣在3種邊界下的變形破壞情況，如圖4（a）、（b）、（c）所示：巖樣在柔性膜邊界和“Shell”邊界下都發(fā)生了明顯的橫向膨脹現(xiàn)象，這也與圖4（d）中甘霖等［24］所做的花崗巖三軸試驗中花崗巖發(fā)生橫向膨脹的現(xiàn)象基本吻合，然而巖樣在剛性邊界下卻沒有發(fā)生橫向膨脹的現(xiàn)象，這是一定程度上由于剛性伺服機制，為了保持伺服圍壓穩(wěn)定在目標區(qū)間內(nèi)，整個圓柱形剛性邊界進行整體同步運動，因此在圍壓相同、軸應變（ε1）相同的情況下，相對于Shell邊界而言，剛性邊界下巖樣的徑向應變（ε3）大于“Shell”單元邊界下巖樣的徑向應變（ε3），而根據(jù)MARTIN［25］所提出的體積應變公式（3），剛性邊界下的巖樣體積應變（εv）應大于Shell單元邊界下的巖樣體積應變（εv）。

圖4 試樣變形形態(tài)對比Fig.4 Comparison of deformation of different samples

為對上述猜想進行驗證，列出模擬巖樣的體積應變曲線，鑒于文章篇幅，本文僅列出X61巖樣在3種邊界情況下的體積應變曲線，如圖5所示。從圖5中可以看出，在同一軸應變時，X61巖樣在剛性邊界下的體積應變大于在Shell單元邊界下的體積應變，從而驗證了上訴猜想的合理性。

圖5 3種邊界下X61試樣的應力-體積應變曲線Fig.5 Stress vs volume strain curves of X61 specimens under three boundaries

3.2 巖樣應力-應變曲線

三軸試驗中巖樣的應力-應變曲線是反映不同邊界對花崗巖三軸試驗影響的主要形式，圖6為27組巖樣三軸模擬試驗的應力-應變曲線圖。

從圖6中可知：不同邊界下巖樣的三軸壓縮曲線的應力-應變曲線的變化趨勢大體上是一致的，曲線在峰值應力前近似呈線性變化，巖樣進入屈服階段后，曲線斜率逐漸變小，峰后階段巖樣強度均出現(xiàn)明顯的下降。相同加載速率下，巖樣在3種邊界條件下的峰值應力隨圍壓的增加而增長，其中巖樣在柔性膜邊界下的峰值應力增加幅度最為劇烈。相同圍壓條件下，巖樣在剛性邊界下的峰值強度受加載速度的增加影響最大，在柔性膜邊界下的峰值強度受加載速率的增加影響最小。同一加載條件下，巖樣在柔性膜邊界下的殘余強度大于在其它2種邊界下的殘余強度，這主要存在2個原因：一是由于根據(jù)柔性膜伺服理論計算出的柔性膜表面積略大于實際的柔性膜邊界表面積，造成巖樣在柔性膜邊界下所受的實際圍壓大于設定的目標圍壓，尤其是目標圍壓越大的情況下，巖樣所受的實際圍壓與目標圍壓差距越大；二是由柔性膜邊界中的膜顆粒本身具有質(zhì)量，對巖樣形成了端部約束［14］。因此相對比其他兩種邊界下的巖樣，巖樣在柔性膜邊界下的峰值應力易受圍壓影響且殘余強度普遍大于其他兩種邊界下巖樣的殘余強度。同時發(fā)現(xiàn)在同一加載條件下，巖樣在剛性邊界下的峰值應力大于在Shell邊界下的峰值應力，這一定程度上是由于剛性伺服函數(shù)的連續(xù)運行之間存在時間間隔以及剛性伺服函數(shù)確定圍壓值時所存在的誤差造成的［21］。

3.3 不同邊界對巖樣力學強度影響

根據(jù)圖6中的數(shù)據(jù)繪制出3種邊界下巖樣峰值應力、彈性模量隨加載速率變化的分布曲線，如圖7、8所示。

圖6 3種邊界下不同圍壓及加載速率下試樣應力-應變曲線Fig.6 Axial and radial stress vs strain curves under three boundaries at different confining pressure and loading velocity

從圖7中可以看出3種邊界下，隨著圍壓的增加，不同加載速率下的巖樣峰值應力均大幅增加，各力學參數(shù)之間與圍壓之間滿足線性正相關關系。Shell邊界下的巖樣在加載速率為0.5 mm/s，圍壓為6 MPa時峰值應力為222 MPa，當圍壓增加到12、24 MPa時，峰值應力增加到289、415 MPa，增長了約30%、86%。而當加載速率增加到1.0、1.5 mm/s時，峰值應力增加到239、253 MPa，增長了約7.6%、13.9%。柔性膜邊界下的巖樣在速率為0.5 mm/s，圍壓為6 MPa時，其峰值應力為232 MPa，當圍壓增加到12 、24 MPa時，柔性膜邊界下巖樣的峰值應力增加到352、575 MPa，增長了約51%、147.8%，而當加載速率增加到1.0、1.5 mm/s時，峰值應力增加到242、254 MPa，增長了約4.3%、9.5%。剛性邊界下的巖樣在速率為0.5 mm/s，圍壓為6 MPa時，其峰值應力為270 MPa，當圍壓增加到12、24 MPa時，剛性邊界下巖樣的峰值應力增加到331、447 MPa，增長了約22.6%、65.6%，而當加載速率增加到1.0、1.5 mm/s時，峰值應力增加到293、318 MPa，增長了約8.5%、17.8%。對比發(fā)現(xiàn)，對于3種邊界下巖樣的峰值應力而言，由圍壓效應造成的影響大于由加載效應造成的影響，其中柔性膜邊界下巖樣的峰值應力受所受圍壓影響最大，剛性邊界下巖樣的峰值應力受加載速率影響最大。

圖7 三種邊界下巖樣峰值應力與加載速率關系Fig.7 Relationship between peak stress and loading rate of rock samples under three boundaries

從圖8中可以看出，3種邊界下巖樣的彈性模量都隨著圍壓或者加載速率的增加而增加，但增長幅度不大，其中Shell邊界下巖樣的彈性模量受加載速率影響較大，柔性膜邊界和剛性邊界下巖樣的彈性模量受圍壓影響較大。

圖8 三種邊界下巖樣彈性模量與加載速率關系Fig.8 Relationship between elastic modulus and loading rate of rock samples under three boundaries

3.4 巖樣力鏈分布

荷載在試樣通過顆粒間的接觸力鏈傳遞，分析接觸力鏈可從細觀角度上對比分析不同邊界對巖樣三軸壓縮實驗的影響，因此繪制出3種邊界下巖樣在軸應變（7%）的接觸力鏈圖，如圖9所示，其中線條粗細代表接觸力鏈強度大小。

從圖9可知，巖樣破壞后在不同邊界下的力鏈強度、力鏈位置分布均存在差異。對于同一加載條件下的巖樣，剛性邊界下巖樣在邊界接觸處存在一條強度遠大于其他接觸力鏈的接觸力鏈，這說明此時巖樣在剛性邊界下已經(jīng)發(fā)生了應力集中現(xiàn)象，而Shell邊界下巖樣則未出現(xiàn)此情況。而對于Shell邊界下巖樣，其接觸力鏈強度較為均勻，較少出現(xiàn)個別接觸力鏈強度遠大于其他接觸力鏈強度的現(xiàn)象。

4 結論

本文使用顆粒流程序研究了剛性、柔性膜、“Shell”3種邊界對模擬花崗巖常規(guī)三軸試驗的影響。針對3種邊界下花崗巖的宏觀破壞形態(tài)和應力應變關系、力鏈分布等模擬結果進行了分析。通過本文研究，得出結論如下：

（1）邊界影響巖樣的加載效應和圍壓效應。巖樣在不同邊界下，其受到加載速率和圍壓的影響程度不同，從上面的三軸模擬實驗中可以看出，相同圍壓下，剛性邊界下巖樣的峰值應力易受加載速率的影響，相同加載速率下，柔性膜邊界下巖樣的峰值應力受圍壓影響最大。

（2）邊界對巖樣的破壞變形有顯著影響。從圖4、圖9中可知，柔性膜和Shell邊界下的巖樣會發(fā)生橫向不均勻膨脹破壞，與花崗巖室內(nèi)實驗現(xiàn)象相符，而剛性邊界下的巖樣整體上呈現(xiàn)整體壓縮破壞，不太符合花崗巖室內(nèi)實驗現(xiàn)象。

（3）剛性邊界和柔性膜邊界在三軸加載過程中無法對巖樣施加精準且穩(wěn)定的圍壓。其中剛性邊界由于剛性伺服函數(shù)的連續(xù)運行之間存在時間間隔以及剛性伺服函數(shù)確定實際圍壓值時所存在的細小誤差，導致巖樣在剛性邊界下無法獲得穩(wěn)定的圍壓。而柔性膜邊界在三軸加載過程中根據(jù)柔性膜伺服理論計算出的柔性膜表面積略大于實際的柔性膜邊界表面積，導致巖樣在柔性膜邊界下所受的實際圍壓大于目標圍壓，所以巖樣在柔性膜邊界下的峰值應力易受目標圍壓影響且殘余強度普遍大于其他2種邊界下巖樣的殘余強度。

（4）剛性邊界對巖樣在三軸加載過程中應力集中現(xiàn)象的產(chǎn)生具有影響。從圖9中可知，巖樣在剛性邊界下發(fā)生應力集中的概率較大，這是由于對于Shell邊界和柔性膜邊界而言，二者在力學性質(zhì)上具有6個自由度，而剛性邊界由于其伺服機制，其在力學性質(zhì)只具有一個自由度，因此相對于其他2種邊界，巖樣在剛性邊界下發(fā)生應力集中的概率較大。

圖9 3種邊界下巖樣接觸力鏈圖（軸應變?yōu)?%）Fig.9 Contact force chain of rock samples under three boundaries（ε1=7%）

綜合對比剛性邊界、柔性膜邊界和Shell邊界，發(fā)現(xiàn)在Shell邊界下，巖樣圍壓施加簡單精準且在加載過程中其所受圍壓恒定不變，應力集中發(fā)生概率較小，因此在花崗巖離散元常規(guī)三軸模擬時，Shell邊界是合適的選擇。