砂巖直接剪切試驗三維離散元分析

楊開新，陳有亮，蔣明鏡,2,3,5，王華寧

(1.上海理工大學 土木工程系,上海 200093；2.同濟大學 土木工程防災重點實驗室,上海 200092；3.同濟大學 土木工程學院,上海 200092；4.同濟大學 航空航天與力學學院,上海 200092；5.天津大學 建筑工程學院,天津 300072)

1 研究背景

隨著 “一帶一路”建設的不斷推進，交通、水利水電工程和資源開采等正在向深部地層快速發展，超千米深的隧道、礦井、地下洞室等工程越來越普遍。對巖體強度和變形特性以及破壞機理的研究是巖體穩定性分析的基礎，進而提高巖體工程設計的合理性[1]。直接剪切試驗能夠有效地模擬剪切面的受力狀況，并能充分考慮在不同法向應力下的裂紋擴展貫通情況，是研究巖土類材料剪切力學特性的有效手段之一[2-3]。由于直剪試驗結果對試樣制備技術和邊界加載條件十分敏感，并且試樣破壞時裂紋擴展情況等微觀信息也無法全面地獲取，因而許多研究人員青睞于采用數值方法來研究巖土類材料的剪切特性，其中由Cundall等[4]提出的離散單元法(discrete element method,DEM)是一種研究巖石剪切特性和破壞機理的可行方法。夏才初等[5]、宋英龍等[6]采用PFC2D生成粗糙節理面并模擬其剪切性質。結果表明：宏觀剪切區域的產生主要是由細觀剪裂紋的累積形成的，剪切破壞區域集中在爬坡效應顯著的位置，剪切過程中形成的剪裂紋數量在出現剪切峰值應力后顯著增加，此時粗糙節理面的破壞最為顯著。Jiang等[7-9]將基于室內試驗得到的巖石微觀力學模型植入離散元軟件，模擬了由化學風化導致的共面非貫通與非共面非貫通節理巖體的直剪試驗，分析了預制節理連通率、節理傾角、法向應力和節理間距對巖體力學特性和貫通模式的影響，并揭示了相應的宏微觀機制。結果表明：節理巖體的全剪切過程包括彈性剪切、裂縫擴展、巖橋破壞和貫穿不連續4個階段；非貫通非共面節理巖體試樣的破壞由節理端點處因拉應力導致的翼裂紋開始，隨著剪切位移的增加，翼裂紋朝著最大壓力方向擴展，并最終貫通巖橋。Shang等[10]基于離散單元法研究了邊界條件(恒定法向載荷(constant normal lad,CNL)和恒定法向剛度(constant normal stiffeness,CNS))對直剪初始巖石不連續面破壞機理的影響。結果表明：CNL和CNS邊界條件顯著影響初始巖石不連續面的剪切特性；施加相同的法向應力下，與CNL直剪的結果相比，CNS直剪中的峰值剪切應力顯著增加，與新產生的微裂紋的開放和巖橋內破裂帶的產生有關，導致正應力的急劇增加。Jiang等[11]研究了直接剪切試驗中砂土的臨界狀態特性和主應力旋轉問題，并分析了作者提出的考慮抗滾動和抗扭轉的三維接觸模型的影響。結果表明：在剪切應變達到15%以后，剪切區域達到臨界狀態，而整個試樣并未達到臨界狀態；主應力旋轉趨勢與剪切應力的變化相似；考慮抗滾動和抗扭轉的三維接觸模型可以很好地反映砂土的剪切力學行為及變形特性。

上述研究多集中于節理巖體的剪切力學行為及變形特性，能夠較好地反映節理巖體的剪切特性以及裂紋發展情況，但對于直剪試驗過程的宏、微觀分析仍不夠全面，且多局限于二維數值模擬層面，對于深部高地應力作用下巖石的剪切力學特性的研究明顯不足。針對以上問題，本文基于本團隊的三維砂巖半經驗半理論膠結接觸模型，將該模型植入商業軟件PFC3D中，進行了不同法向應力下的直剪數值試驗，對砂巖的剪切特性進行了全面地宏、微觀分析，試圖探尋巖石在直剪試驗過程中的裂紋擴展機制，為研究深部巖體的剪切力學行為與深部巖體工程設計提供借鑒。

2 接觸模型與微觀參數

為了更加合理有效地模擬出真實巖石地力學行為與變形特性，本文微觀接觸模型的半經驗半理論膠結強度準則[12]曲線如圖1所示，其中圖1(a)、1(b)和1(c)分別描述了膠結在單純法向應力作用下的抗剪、抗彎和抗扭強度。剪-彎-扭耦合強度包面采用Shen等[13]提出的表達式，如圖1(d)所示，其描述了膠結在法向、剪切、彎轉和扭轉共同作用下的強度。由圖1(d)可知，剪-彎-扭耦合強度包面的大小取決于其3個方向的軸長，而這3個方向的軸長又取決于膠結的法向應力，因此剪-彎-扭耦合強度包面的大小取決于膠結的法向應力。

圖1 半經驗半理論膠結強度準則曲線[12]

本團隊在半經驗半理論膠結強度準則的基礎上，將其引入三維顆粒含抗彎轉和抗扭轉接觸模型[14]形成新的膠結接觸模型，并將該新膠結接觸模型命名為半經驗半理論膠結接觸模型。新模型吸收了雙球膠結接觸試驗的經驗，并同時考慮了材料的理論強度，因此與單純的經驗模型和理論模型相比更具合理性也更加可靠。本文參數標定的對象為深部Berea砂巖[15]，通過Bera等[16]的CT掃描圖(圖2)間接獲得其顆粒級配曲線(圖3)。本文選擇深部Berea砂巖為目標巖石主要基于以下幾個方面的考慮：(1)此種巖石的室內三軸試驗資料比較完整，包含9個圍壓的全過程試驗數據，并且圍壓的范圍符合高應力的要求；(2)能夠獲得此種巖石的微觀CT掃描圖；(3)砂巖的成巖過程相對易于模擬，有利于定量模擬的成功實施。為了提高計算效率，在進行參數標定[12]時采用1.27×104顆粒的小尺寸試樣來模擬巖石的常規三軸試驗。試樣長寬高的比例約為1∶1∶2.5，參數標定模擬試樣如圖4所示。

圖2 Berea砂巖CT掃描圖[16] 圖3 Berea砂巖顆粒級配曲線 圖4 DEM參數標定模擬試樣[12]

以Berea砂巖室內三軸壓縮試驗[16]分別在10和40 MPa兩種圍壓下的全過程應力-應變曲線為逼近對象，通過對比分析不斷調整微觀參數，并以其他圍壓下的全過程應力-應變曲線為輔助，最終選定試樣的微觀參數[12]如表1所示。

3 結果與分析

3.1 直接剪切數值試驗模型

巖石的抗剪強度是巖石對剪切破壞的極限抵抗能力，通常用直剪試驗測定。本文進行了高應力下砂巖的直剪數值試驗，采用正方體試樣，包含4×104顆粒，長、寬、高的比例為1∶1∶1，試樣尺寸為5.5 mm×5.5 mm×5.5 mm，制樣時采用與參數標定時一致的顆粒級配，試樣模型如圖5所示。具體直剪試驗的DEM模擬分為以下4個步驟：(1)分層欠壓成樣。采用分層欠壓法[17]成樣，試樣分5層，每層顆粒數為8 000，試樣目標孔隙比為0.8，保證了成樣后顆粒的重疊量足夠小，且考慮了目標巖石為孔隙相對較多的Berea砂巖(實測孔隙比為0.215±0.007)[16]。這里造成制樣孔隙比的差異是由于離散元中膠結部分是虛擬的，沒有計算體積，只考慮了力學響應。(2)預壓固結。本文模擬的對象為深部巖石，假定目標試樣地應力在垂直方向上為30 MPa，應力狀態為K0狀態。在分層欠壓穩定后，固定左右前后四面墻的位置，通過上下墻伺服施加豎向壓力30 MPa并達到穩定。(3)生成膠結。在試樣預壓固結穩定后，在顆粒之間一次性生成膠結，模擬Berea砂巖的膠結形成過程。然后將膠結試樣卸載至等向受壓應力狀態，壓應力取1 kPa，模擬巖石的取樣卸載過程(理論上壓應力應該卸載至0，數值試驗中取1 kPa僅僅是為了便于伺服控制)。(4)加載。令上下墻以0.025%/s的應變速率相向運動，模擬巖石直剪試驗的加載過程。

圖5 直剪數值試驗試樣模型圖

3.2 宏觀結果分析

模擬了法向應力為20、30和40 MPa下的直剪試驗全過程，其剪應力-剪應變曲線如圖6所示。由圖6可知，剪應力-剪應變曲線可以分為3段：(1) 線性剪切階段。剪應力隨著剪應變的增大呈線性增大，在此階段的末端，曲線變為非線性，表明此時剪切面上已經有膠結破壞產生，粗糙的區域開始磨平[18]。(2) 應力衰減階段。此階段內，剪應力隨著剪應變的增加而減小，對應著膠結破壞繼續發生，剪切面上的粗糙區域繼續破壞。(3) 殘余強度階段。進入此階段后，剪應力達到一個相對穩定的值，即殘余強度，表明試樣剪切面上粗糙區域的磨平效應基本結束。

圖6 不同法向應力下直剪試驗應力-應變曲線 圖7 直剪試驗全過程中各能量(功)的演化曲線

離散元模擬膠結散粒體材料受力變形過程中涉及的能量(功)[19]包括外界輸入功、儲存于膠結接觸處的膠結彈性能、儲存于無膠結接觸處的無膠結彈性能、膠結破壞耗散能、黏附型膠結接觸處的膠結耗散能、無膠結耗散能、數值阻尼耗散能和顆粒動能，具體含義及表達方式見參考文獻[19]。基于上述基礎，將能量劃分為顆粒彈性能、顆粒膠結彈性能、顆粒摩擦耗能、顆粒彎轉耗能、顆粒扭轉耗能、阻尼耗能、顆粒動能。圖7給出了直剪試驗全過程中各能量(功)的演化情況。分析圖7可知，在線性剪切階段，外界做功主要用于增加顆粒彈性能與顆粒膠結彈性能。隨著剪應變的增加，進入應力衰減階段，顆粒摩擦耗能迅速增加，阻尼耗能也略微增加，表明此時，外界做功主要被顆粒摩擦所消耗，而顆粒彈性能與顆粒膠結彈性能基本達到穩定。在進入殘余階段后，阻尼耗能也基本達到穩定，顆粒摩擦耗能繼續增加。整個DEM直剪試驗過程中，顆粒彎轉耗能、顆粒扭轉耗能與顆粒動能(準靜態)非常小，可忽略不計，且由能量總和與外界做功近乎相等可知，離散元模擬過程中始終存在著能量守恒。總體而言，整個DEM直剪試驗過程中，外界做功主要轉化為顆粒摩擦耗能，即剪切過程的本質是剪切面上粗糙區域的磨平，磨平的過程中伴隨著大規模的膠結破壞，顆粒間因而產生不平衡力，但整個顆粒體系內部結構處于快速調整狀態。

3.3 微觀結果分析

基于半經驗半理論膠結接觸模型，對微觀膠結的破壞形式采用了系統的分類[20]：(1)根據膠結破壞時的法向應力狀態，將膠結分為壓區破壞和拉區破壞，即當膠結破壞時，如果膠結法向應力為壓應力稱為壓區破壞，如果膠結法向應力為拉應力稱為拉區破壞。(2)根據膠結破壞準則耦合強度包面[13]將膠結破壞形式分為剪破壞、彎破壞和扭破壞，其中剪切項占比最大時稱為剪破壞，彎轉項占比最大時稱為彎破壞，扭轉項占比最大時稱為扭破壞。結合以上兩種分類方法最終可以將膠結破壞分為壓應力狀態下的剪、彎、扭破壞和拉應力狀態下的剪、彎、扭破壞，即壓剪、壓彎、壓扭、拉剪、拉彎和拉扭6種破壞類型。

圖8給出了直剪數值試驗中各破壞類型的膠結破壞數與剪應變的關系曲線。由圖8可知，在線性剪切階段，膠結破壞數目增長迅速，表明剪切面上粗糙區域的磨平會產生大量膠結破壞；在應力衰減階段，膠結破壞增長率較小，膠結破壞數目略微增多；進入殘余階段時，膠結破壞數目隨著應變的增加近似呈線性增加，此時剪切面上粗糙區域的磨平基本完成，但膠結破壞總數目會隨著剪應變的增加而不斷增加。直剪數值試驗全過程中，膠結破壞形式按上述(1)分類時，膠結破壞形式主要是拉區破壞；按(2)分類時，膠結破壞形式主要是彎破壞，其次是剪破壞，扭破壞最小。結合來看，膠結拉區破壞數目遠遠大于膠結壓區破壞數目，無論在拉應力或壓應力狀態下膠結的破壞形式均由彎破壞形式主導。

圖8 各破壞類型的膠結破壞數與剪應變的關系曲線 圖9 直剪數值試驗試樣裂紋擴展狀態(剪應變為7%)

為了研究直剪數值試驗過程中的裂紋擴展機理，本文主要采用格里菲斯強度理論[21]來解釋巖石內部的裂紋擴展現象，說明巖石的剪切破壞機理。格里菲斯認為在外力作用下，微裂紋的存在改變了材料內部的應力狀態，隨著外力的逐漸增大，裂紋將沿著與最大拉應力成直角的方向擴展。

圖9展示了直剪數值試驗剪應變為7%時的試樣裂紋擴展情況。圖9中位于中央的為試樣裂紋的3D顯示狀態，將其按圖中箭頭顯示方式旋轉0°、90°、180°、270°后得到其平面裂紋顯示狀態，左上與右下分別為其上視圖與下視圖，圖中黑色裂紋為壓區破壞裂紋，紅色裂紋為拉區破壞裂紋。從圖9中0°與180°視圖來看，裂紋沿x-z平面呈對稱分布，剪切面兩端分布裂紋較少，中間區域分布裂紋較多，形狀近似一個扁平的橢圓；從90°與270°視圖來看，裂紋沿y-z平面呈對稱分布，裂紋在剪切面上呈矩形分布，存在著明顯的破碎帶，即中間剪切面是壓區破碎帶，剪切面上下方是拉區破碎帶；從上視圖與下視圖來看，整個剪切面上均存在裂紋，這是二維顯示不能達到的效果，也說明直剪數值試驗的裂紋擴展在殘余強度階段分布于整個剪切面，且裂紋沿x-y平面呈對稱分布。綜合來看，試樣直剪數值試驗的裂紋擴展具有明顯的空間對稱性，分布形狀近似呈橢球狀，且對稱分布于剪切面。

將本文直剪數值模擬試樣裂紋擴展狀態與Zhu等[22]得出的砂巖剪切力學行為試驗結果相比較，如圖10所示，圖10(b)中紅色虛線內部為裂紋擴展區域。

圖10 本文直剪數值試驗試樣裂紋擴展狀態與參考文獻[22]試驗結果比較

由圖10可看出，盡管兩者加載的法向應力大小不一致，但處于高法向應力下的砂巖剪切破壞特性基本一致，即最終形成一個截面呈扁平的橢圓狀分布的裂紋域，這與低法向應力下的砂巖剪切破壞特性明顯不同。

為了進一步了解剪切面上裂紋擴展演化情況，對試樣剪切面進行不同剪應變階段的切片觀測，結果如表2所示。由表2可知，在線性剪切階段，裂紋隨著剪應變的增加而增加。在剪應變為0.2%時，即準靜態加載初期，裂紋還未產生；在剪應變為0.4%時，出現少量裂紋分布于剪切面的一端；在剪應變為0.6%時，裂紋進一步發展，壓區裂紋主要分布于剪切面兩端，而此時拉區裂紋主要分布于剪切面內部；在剪應變為0.8%時，裂紋大量產生，分布規律與剪應變為0.6%時呈現的分布規律一致，主要是由于在直剪試驗的應力狀態下，最大拉應力近乎沿剪切面方向。進入應力衰減階段，裂紋數量而繼續增加。進入殘余強度階段(應變為3%～7%)，裂紋數量隨著剪應變的增加逐漸增加，這與圖9顯示的裂紋擴展規律一致，在整個剪切面上無論是兩端還是中間，均分布著大量的壓區和拉區的裂紋，布滿整個剪切面。

表2 直剪數值試驗試樣剪切面裂紋隨剪應變擴展演化狀況

為了更加細致地研究直剪試驗全過程中的裂紋擴展情況，試驗輸出了剪應變為0.4%、0.6%、1%、3%和7%時的巖體破碎、裂縫發展、力鏈及顆粒位移場狀態圖，如表3所示。根據前文直剪試驗應力-應變曲線分析結果，剪應變為0.4%、0.6%時屬于線性剪切階段，剪應變為1%時屬于應力衰減階段，剪應變為3%、7%時屬于殘余強度階段。由表3可見，在剪應變為0.4%時，巖體試樣部分區域發生破碎，但破碎區域分布比較零星，多分布于試樣邊界處；在剪應變為0.6%時，破碎繼續發展，但規模不大；在剪應變為1%時，進入應力衰減階段，此時大規模膠結破壞已經發生，剪切面處的破碎面也已經初步形成；隨著剪應變不斷增加到達殘余強度階段(7%)時，剪切面呈橢球狀的破碎區域基本穩定。

表3 直剪數值試驗不同剪應變時巖體試樣破碎、裂縫發展、力鏈與顆粒位移場狀態

圖9已經展示了殘余階段剪應變為7%的裂紋擴展情況，下面將依據表3分析直剪數值試驗全過程中試樣的空間裂紋擴展情況。在線性剪切階段，剪應變為0.4%時，巖體試樣開始出現裂紋，主要分布于試樣左下與右上區域，剪切面上出現少許裂紋；在剪應變為0.6%時，裂紋繼續擴展，試樣左下與右上區域的裂紋數量增加并不明顯，但是剪切面上的裂紋明顯增加，且拉區裂紋占主導地位。在應力衰減階段，剪應變為1%時，剪切面裂紋區域基本形成，可知直剪數值試驗中裂紋的擴展機制先是從剪切面周圍應力集中區域發生，然后剪切面上由于拉應力的影響也開始產生裂紋，最后剪切面上的裂紋之間相互貫通，形成一片裂紋區域。在殘余強度階段，剪應變為3%時，裂紋從剪切面開始沿垂直于最大拉應力方向且平行于大主應力方向發展，即按照格里菲斯強度理論所示方向發展(圖9)；在剪應變為7%時，裂紋發展基本達到穩定。

表3給出了直剪數值試驗全過程的力鏈和顆粒位移場的發展情況，力鏈中黑色表示壓力鏈，紅色表示拉力鏈。隨著剪應變的增加，壓力鏈由豎直狀態逐漸向著小主應力方向發生偏轉，而拉力鏈向著大主應力方向偏轉，剪應變越大，偏轉現象越明顯。顆粒位移場顯示的是試樣中間沿y-z平面的一個薄長方體。顆粒位移場是反映顆粒位移方向以及大小的一個矢量場，是研究裂紋擴展的一個重要手段，本研究不考慮顆粒位移的大小，只考慮顆粒位移方向，即裂紋發展趨勢，來揭示直剪數值試驗中巖石內部裂紋擴展機理。下面將依據表3分析直剪數值試驗全過程中試樣顆粒位移場的發展情況。在線性剪切階段，剪應變為0.4%時，試樣內部發生以剪切方向為主導的偏轉，試樣上下兩端偏轉程度較小，而剪切面處的偏轉程度較大，偏轉中心為剪切面的中心區域；在剪應變為0.6%時，試樣內部偏轉程度加劇，且偏轉中心發生傾斜，傾向于大主應力方向。在應力衰減階段，剪應變為1%時，大規模膠結破壞已產生，剪切面處的破碎帶已經基本形成，位移場中與之對應的是弧形滑動面的基本形成。在殘余強度階段，剪應變為3%時，貫通的滑動面的形狀十分穩定，但在剪應變達到7%時，貫通的滑動面的寬度變大，滑動面開始從剪切面分別向上、下發展，這與巖體破碎、裂紋擴展所表現的趨勢一致，說明利用顆粒位移場來研究試樣內部位移場也不失為一種良好的手段。

圖11給出了本文直剪數值試驗及參考文獻[22]中砂巖剪切試驗所得出的峰值剪應力與法向應力的關系曲線。由圖11可知，峰值剪應力隨著法向應力的增加而近似呈線性增加，采用莫爾庫倫強度準則對其擬合，發現擬合效果較好。與參考文獻[22]中砂巖的試驗結果擬合曲線相比，兩條強度包線近似平行，本文砂巖內摩擦角擬合結果小4.82%，黏聚力大28.87%。由此可知，本文砂巖直剪數值模擬結果與試驗結果在定性上基本一致，在定量上有較小差值。

圖11 直剪數值試驗峰值剪應力-法向應力變化曲線

基于上文的研究可以得出，高應力下砂巖的剪切力學行為表現出應力軟化的特質，與室內試驗結果[21]一致性較好。而不同應力下的裂紋擴展存在區別，高法向應力下，砂巖在剪切荷載作用下，會形成一個截面呈扁平橢圓狀分布的裂紋域，而低法向應力下，會形成沿剪切面的寬度較小的裂紋域[21]，使得深部巖體的剪切破壞機制明顯區別于淺層巖體。

對于砂巖三維直剪數值試驗的裂紋顯示問題，本文還存在不足，未能更加清晰地給出直剪過程中巖石試樣內部中裂紋在空間的具體擴展模式，將在以后進一步研究。

4 結 論

本文通過離散單元法，采用基于本團隊的半經驗半理論膠結接觸模型，模擬了砂巖的直接剪切試驗全過程，對數值試驗結果進行了詳細的宏、微觀分析。主要結論如下：

(1)砂巖的直剪數值試驗全過程可分為3個階段：線性剪切階段、應力衰減階段和殘余強度階段。

(2)砂巖的直剪數值試驗過程中，外界做功主要轉化為顆粒摩擦耗能，離散元模擬過程中始終存在著能量守恒；膠結破壞數目會隨著剪應變的增加而增加，膠結拉區破壞數目大于膠結壓區破壞數目；在拉或者壓應力狀態下膠結的破壞形式均由彎破壞形式主導。

(3)砂巖直剪數值試驗得出的巖體裂紋擴展機制是：裂紋首先從剪切面兩端產生，拉區裂紋產生于剪切面中間，向剪切面兩端以及上、下側擴展，壓區裂紋從剪切面兩端向中間擴展，最后拉區裂紋與壓區裂紋相互貫通布滿整個剪切面，在空間上最終形成一個截面呈橢球狀的分布域。

(4)砂巖直剪數值試驗的裂紋擴展結果、峰值剪應力結果與砂巖剪切試驗結果一致性較好。