3DEC中紅砂巖節理模型開發及數值模擬

劉清樸，侯 迪，榮 冠

(1. 武漢大學 水資源與水電工程科學國家重點實驗室，武漢 430072; 2. 武漢大學 水工巖石力學教育部重點實驗室，武漢 430072； 3. 貴州省水利水電勘測設計研究院，貴陽 550002)

0 引 言

天然巖體中存在著大量結構面，如節理、裂隙等。節理力學性質較完整巖石更為復雜，顯著影響整個巖體的力學特性。在工程荷載作用下，節理及其力學性質影響著巖體中的應力分布，進而影響巖體工程的穩定性。同時，節理的抗剪強度模型是巖石力學的基礎問題和研究的難點。

目前，研究者提出很多節理峰值抗剪強度模型。其中，Patton[1]根據齒狀節理直剪試驗，提出了反映節理起伏的線性抗剪強度模型。Ladanyi[2]考慮法向應力與節理鋸齒起伏的互鎖作用關系，提出了節理非線性抗剪強度模型。Barton[3,4]通過引入節理粗糙度系數(JRC)，提出了著名的抗剪強度模型。Plesha[5]假定節理為均勻對稱的齒形起伏面，提出了一個在力學上更嚴格的抗剪強度模型。但上述模型都是基于巖石節理二維形貌參數提出的峰值強度模型，而二維模型難以合理和全面地反映出巖石節理各向異性等力學特性。因此，深入研究節理三維形貌的表征方法，建立巖石節理三維強度模型成為克服節理模型難以準確反映其真實力學特性的關鍵。Jing[6]研究了節理面的各向異性摩擦角，并參考Plesha模型中剪脹角演化算法，提出了考慮三維各向異性的節理強度模型。夏才初等[7]根據實驗結果給出了節理抗剪強度與形貌統計參數Z2的經驗關系。Grasselli[8,9]通過大量巖石節理直剪實驗，認為在剪切過程中只有面對剪切方向的微單元面才相互接觸剪切、磨損，而面向方向的會出現分離，通過研究節理面微元有效剪切傾角與其對應接觸面積的統計數學模型，提出了節理三維粗糙度的擬合算法，建立了節理三維抗剪強度模型。唐志成等[10,11]在Grasselli強度模型的基礎上，考慮峰值剪脹角與三維形貌特征關系，提出一個可以反映剪脹效應的峰值強度修改模型。

由Cundall于1971年提出來的離散元方法允許不同塊體之間產生大的位移甚至脫離接觸，并在計算過程中程序能夠自動判斷塊體之間接觸狀態的變化，適用于解決大位移(諸如滑移、張開或閉合等)的非線性問題，在工程巖體數值分析計算中有著獨特的優勢。目前離散元程序中節理模型主要有庫倫滑移節理模型(Coulomb Slip Model)、連續屈服節理模型(Continuously Yielding Joint Model)、BB模型(Barton-Bandis Joint Model)，但在工程巖體計算分析中，現有模型存在著一些不足。其中，庫倫滑移節理模型假定法向變形為線彈性、切向變形為理想彈塑性，所以不能合理的模擬出節理的非線性力學特性。由Cundall和Hart[12]提出的連續屈服節理模型，采用指數方程來表征節理法向變形特征，而由Bandis等[13]提出的法向雙曲線方程應用較為廣泛，且連續屈服節理模型計算參數需依靠經驗確定，無法在工程計算中廣泛應用。BB模型剪切方向采用Barton基于節理二維形貌參數提出的剪切強度模型，Hong等[14]指出二維形貌參數會低估節理面的粗糙度，BB模型難以全面反映三維形貌影響下巖石節理剪切力學特性。高艷華等[15]基于試驗數據擬合，通過在三維離散元軟件3DEC中編寫FISH語言動態改變調用的內置節理模型參數，研究了節理剪切力學特性，但通過FISH進行動態修改的方法將大大增加計算時間，特別是在大型巖體工程計算中將影響計算效率。如何基于三維形貌參數把節理非線性變形特征引入到三維離散元軟件3DEC中成為了大型工程更加合理模擬分析的關鍵。

綜上所述，本文采用人工劈裂的紅砂巖節理進行直剪試驗，提出了反映節理三維形貌的紅砂巖節理峰值抗剪強度模型。并根據提出的節理峰值抗剪強度修正描述了剪切過程，采用C++語言編寫了相應節理模型[16]。在3DEC主程序中模擬了單軸壓縮試驗和直剪試驗，驗證了所提模型可以合理反映出紅砂巖節理力學特性，為后續研究優勢結構面控制下的邊坡穩定性及在開挖、支護等工程作用下圍巖洞室穩定性分析、安全評價提供基礎。

1 節理直剪試驗

1.1 試件制備

由于天然巖石節理采集難度較大，因此采用人工劈裂節理進行直剪試驗研究。通過采集天然紅砂巖制備成尺寸為200 mm×100 mm×100 mm巖石試件(編號為hs-1、hs-2、hs-3、hs-4、hs-5)，并將其劈裂獲得新鮮耦合巖石節理(見圖1)。試驗中使用的紅砂巖試件采自湖北鄂州。試件巖石新鮮，顏色為磚紅色，主要組成礦物為長石和石英，細砂質結構。利用室內巖石取芯機，對采集的紅砂巖鉆取直徑50 mm 、長度為100 mm 的標準巖樣。通過在三軸伺服試驗機上進行單軸壓縮試驗，測量出紅砂巖的抗壓強度、彈性模量及泊松比。通過巖樣用巴西劈裂法測定紅砂巖的抗拉強度。并通過對平滑節理面進行反復直剪試驗確定試樣基本摩擦角。本次試驗采用試件的基本力學參數為：單軸抗壓強度為65 MPa，抗拉強度2.6 MPa，節理面黏聚力為0 MPa，節理基本摩擦角28°，彈性模量為22.4 GPa，泊松比為0.18。

圖1 紅砂巖劈裂節理試件Fig.1 Joint specimens of sandstone

1.2 試驗方法

試驗采用五組紅砂巖試件，在耦合狀態下使用JAW-1000型巖石節理剪切-滲流耦合試驗機(見圖2)進行剪切試驗。首先按照荷載控制方式施加法向荷載，巖石節理的法向荷載按照σn∶σc(法向荷載：單軸抗壓強度) 分別為0.5%、1.0%、1.5%、2.0%、2.5%施加。每次法向荷載按0.1 kN/min加載速率加至預定法向荷載。然后按照位移加載控制方式進行直剪試驗，剪切速率為0.5 mm/min。當剪切位移達到10 mm時停止試驗，在試驗過程中，試驗系統通過控制軟件可自動采集數據結果。

(2)

式中：m為面對剪切方向的單元數；θ*si為第i個面對剪切方向單元的視傾角；Af為面對剪切方向的三角形單元面積總和；At為粗糙節理面實際面積。

圖2 剪切滲流耦合試驗系統整機圖Fig.2 Photograph of coupled shear-flowtest system

圖3 KEYENCE_LK-G5000掃描儀及其測量現場Fig.3 The scanner and its measurement field

1.3 試驗結果

在法向荷載為0.325、0.65、0.975、1.3、1.625 MPa情況下對紅砂巖節理試件進行直剪試驗。直剪試驗結果及節理三維形貌參數見表1和圖4。

表1 巖石節理形貌參數及強度試驗值Tab.1 Morphology parameters and results of shear tests

圖4 紅砂巖試件剪切位移與剪切荷載曲線Fig. 4 Curves of direct shear tests of sandstone

2 節理峰值抗剪強度模型

根據Barton[4]，Grasselli[8,9]，唐志成[10,11]等研究表明，在巖石節理剪切試驗中隨著法向荷載的增加峰值剪脹角由一個最大值減小到一個恒定值，分析剪脹角與法向荷載的關系可以發現，剪脹角與法向荷載之間符合函數y=1/(1+x)的形式。此外，Grasselli等[9]指出巖石抗拉強度對節理抗剪強度有著重要影響。根據表1和表2分析形貌參數與剪脹角度之間關系可知三維平均傾斜角越大，節理粗糙程度越高，峰值剪脹角度越大；而最大接觸面積比越大，說明節理粗糙程度越低，峰值剪脹角越小。峰值剪脹角與法向荷載之間存在如下邊界條件：

(3)

式中：i0為節理初始剪脹角；σn為節理法向應力。

根據峰值剪脹角曲線變化規律，考慮上述邊界條件，并考慮抗拉強度的影響，而且采用σn/σt作為變量可使得等式左右量綱守恒，提出如下公式：

(4)

基于Parton公式，即可得巖石節理抗剪強度公式為：

(5)

根據直剪試驗數據(見表2)擬合，確定該節理試件經驗系數a=0.057 3，b=0.269 9。

為了驗證新模型的正確性，采用峰值剪切強度公式(5)計算紅砂巖節理試件峰值抗剪強度值，同時采用Barton峰值抗剪強度公式計算，節理試件的JRC值采用Barton[17]的經驗公式計算，公式如下：

(6)

式中：L為節理長度；ξ為節理最大凸起程度。

本文每隔1 mm選取沿剪切方向的共9條剖面線計算其ξ值，節理ξ值由組成節理面的9條剖面線的計算值求平均得到。JRC計算結果與節理峰值剪切強度的計算結果列于表2。

表2 巖石節理強度計算值與試驗值對比表Tab.2 Comparison between calculated and observed results

由表2及誤差平方和結果可以看到，采用新模型計算巖石節理的峰值剪切強度與試驗值更為接近，Barton公式在一定程度上低估了節理的峰值抗剪強度，采用新模型計算節理峰值剪切強度明顯優于Barton公式。

3 節理模型二次開發與數值模擬

3.1 節理模型

國內外學者對節理法向變形規律進行過大量研究，主要是根據室內試驗得到節理法向應力-節理閉合量曲線。目前應用廣泛的為Bandis等[13]基于5種不同巖石(石灰巖、粉砂巖、板巖、砂巖、玄武巖)結構面進行了大量室內試驗，指出節理法向應力與法向變形為非線性關系，提出的如下雙曲線本構方程：

(7)

式中：ΔVj為節理法向變形；節理初始法向剛度Kni和節理最大閉合量Vm可由Bandis等[13]根據節理粗糙系數JRC、節理壁抗壓強度JCS、平均節理縫隙寬度ajn提出的經驗公式求出。其中，Bandis提出的Vm經驗公式中給出了不同循環加卸載次數下計算參數，為了便于計算，本文僅采用第一次加載曲線。

文獻[18]指出對于未風化粗糙節理，影響節理剪切強度的節理因素主要包括節理的粗糙度部分i和基本摩擦角φb兩部分。其中，在剪切過程中基本摩擦角可認為不發生變化，其粗糙度呈遞減變化。Barton等[18]提出節理表面的粗糙度直接影響剪切強度，通過剪切位移與峰值對應位移的比值修正了節理粗糙系數JRC，得到整個剪切過程的剪切強度。因此，根據上述Barton提出未風化粗糙節理在剪切過程中修正方法，結合1.2節中提出的峰值抗剪強度經驗公式(5)，可得出剪切過程中節理切向應力τ為：

τ=σntan (i+φb)

(8)

式中：σn為節理法向應力；φb為基本摩擦角；i=Bip為修正的峰值剪脹角，表示當前狀態下節理形貌剪脹效應的大小，ip為峰值剪脹角，B可由表3中A對應計算得出，A是當前的剪切位移與峰值位移的比值，即：A=u/upeak。

表3 系數B的計算對照表[18]Tab.3 values used for calculating the coefficient B[18]

峰值剪切位移upeak由Barton[18]提出的經驗計算公式可得：

(9)

在剪切過程中，法向位移增量和切向位移增量的比值定義為剪脹角的正切值。參照Barton[4]提出的關于未風化節理剪脹模型，計算剪切過程中剪脹角dn，其公式為：

(10)

在剪切過程中，法向閉合量不僅與法向應力有關，還受到剪切位移的影響，根據剪脹角的定義對其進行修正，剪切位移引起的Δun的計算公式為：

Δun=Δustandn

(11)

3.2 節理模型的二次開發

三維離散元軟件3DEC為用戶提供了兩種方法來修改節理模型。一種是通過3DEC內部的FISH語言命令對計算過程中調用3DEC原有的模型進行動態修改；該方法編寫簡單、容易掌握，但其在計算過程中需不斷對原有模型的進行動態修改，速度比內置DLL文件模型要慢數倍，直接影響計算效率。另一種是通過3DEC為用戶提供的user-defined joint constitutive models 接口，利用C++語言編寫相應的節理模型DLL文件(3DEC程序自帶模型均是以DLL文件的形式提供給用戶)，該方法可以高效、穩定和簡潔地將模型嵌入到3DEC程序中。因此，筆者采用第二種方法編寫新模型的DLL文件。

3DEC為用戶提供了節理模型開發的模板，節理模型開發工作主要是在原有模板基礎上，通過Microsoft Visual Studio 2010平臺，運用C++語言修改編寫新節理模型相應的頭文件(.h)和源文件(.cpp)。最終，編譯形成相應的動態鏈接庫DLL文件，并進行調試修改。

圖5 userjbb節理模型開發流程圖Fig.5 Development flowchart of userjbb joint model

3.3 節理力學特性模擬

為了驗證上述編寫的節理模型正確性，本文采用3DEC模擬了單軸壓縮試驗和直剪試驗，與已有成果進行對比分析。為了保持與室內試驗相一致，模型尺寸采用200 mm×100 mm×100 mm，該模型為含有一水平節理的巖體。對節理下部巖塊底部進行固定，節理上部巖塊施加一水平速度，直剪試驗中采用與室內試驗的相同荷載條件，計算模型如圖6所示。

圖6 數值模擬計算模型Fig.6 Computational model in numerical simulation

計算模型由巖塊和節理兩部分組成。BB模型及userjbb模型計算中巖塊、節理的力學及幾何形貌采用室內試驗中hs-3、hs-4兩組參數(見表4)。

表4 巖體力學及形貌參數Tab.4 Mechanical properties and morphology parameters of rock mass

根據離散元計算原理可知，塊體間的整個接觸又分成了若干個子接觸。因此，在數值模擬計算中，需根據節理面各子接觸作用力計算總的應力，并根據節理面子接觸的相對位移計算出節理平均位移，取法向為例，其表示式如下：

(13)

式中：n為子接觸個數；Fni為第i個子接觸法向力；A為整個節理面面積；uni為第i個子接觸法向位移。同理，可以計算出平均切向應力和平均切向位移。

3.4 模擬結果分析

通過編寫試驗相應的3DEC命令代碼和平均應力和平均位移計算代碼，求出數值模擬過程中整個節理面的平均應力和平均位移。采用userjbb節理模型模擬了單軸壓縮試驗，對比數值模擬和Bandis提出的雙曲線方程理論結果，在三維離散元中可以采用雙曲線方程描述節理法向的變形特征，如圖7所示。兩者計算曲線十分接近，變形規律如下：隨著法向應力的增加，法向閉合量增加速率逐漸增大，法向閉合量不斷增加，最終趨近于最大閉合量。在法向應力較小時，其法向閉合量與法向應力近似為線性。但是當應力較大時，節理法向非線性變形特征尤為突出。

圖7 節理法向應力-法向位移數值試驗與理論結果Fig.7 Numerical and theoretical results of joint normal stress vs normal displacement

采用三維離散元程序中userjbb節理模型及二維離散元程序中BB模型模擬了紅砂巖節理的直剪試驗，并與室內直剪試驗結果進行了對比，如圖8所示。圖8表明數值模擬計算結果、試驗結果均表現了節理剪切應力-剪切變形的總體規律為非線性：在常法向應力下節理剪切應力隨著剪切位移的增加而增加至峰值抗剪強度，之后隨著剪切位移的增加而衰減，衰減到殘余剪切強度后保持不變。其中，采用userjbb節理模型的數值結果和試驗結果吻合較好，而BB模型數值結果明顯低于userjbb節理模型的數值結果和試驗結果。究其原因為：巖體節理表面為凸凹不平的三維粗糙面，其表面上隨機分布了一些微凸體。隨著剪切位移的增加，節理面上的微凸體發生咬合、爬坡等現象。當作用在微凸體上的法向應力大于其所能承受的最大應力時，則會導致節理微凸體的磨損甚至破碎。充分考慮節理形貌特征，將其引入到剪切強度計算中成為合理反映剪切強度的關鍵。而在BB模型中僅考慮了節理的二維形貌，計算結果偏于保守，不能反映出節理三維形貌對剪切強度的影響。

圖8 紅砂巖節理數值模擬與直剪試驗結果Fig.8 Numerical and test results of sandstone joint

4 結 語

(2)分別采用提出的峰值抗剪強度模型和Barton模型計算了相應的峰值強度值，并與室內試驗結果對比分析，分析表明：新模型能更好地描述節理峰值剪切強度，Barton模型在一定程度上低估了節理的峰值抗剪強度。

(3)基于3DEC提供的接口，采用C++語言編寫相應的DLL文件，實現了將建立節理非線性本構模型內嵌于3DEC中。通過數值模擬分析驗證，該方法是一種行之有效的方法。

(4)通過數值模擬與試驗結果對比，可以得出：采用新模型的數值結果與理論、試驗結果吻合較好，可以合理反映出紅砂巖節理力學特性。

□

[1] Patton F D. Multiple modes of shear failure in rock[C]∥1st ISRM Congress. International Society for Rock Mechanics, 1966.

[2] Ladanyi B, Archambault G. Simulation of shear behavior of a jointed rock mass[C]∥The 11th US Symposium on Rock Mechanics (USRMS). American Rock Mechanics Association, 1969.

[3] Barton N. Review of a new shear-strength criterion for rock joints[J]. Engineering Geology, 1973,7(4):287-332.

[4] Barton N, Choubey V. The shear strength of rock joints in theory and practice[J]. Rock Mechanics, 1977,10(1-2):1-54.

[5] Plesha M E. Constitutive models for rock discontinuities with dilatancy and surface degradation[J]. International Journal for Numerical and Analytical Methods in Geomechanics, 1987,11(4):345-362.

[6] Jing L. Numerical modelling of jointed rock masses by distinct element method for two-and three-dimensional problems[D]. Lule?: Lule? Tekniska Universitet,1990.

[7] 夏才初, 孫宗頎. 工程巖體節理力學[M]. 上海：同濟大學,2002.

[8] Grasselli G, Wirth J, Egger P. Quantitative three-dimensional description of a rough surface and parameter evolution with shearing[J]. International Journal of Rock Mechanics and Mining Sciences, 2002,39(6):789-800.

[9] Grasselli G, Egger P. Constitutive law for the shear strength of rock joints based on three-dimensional surface parameters[J]. International Journal of Rock Mechanics and Mining Sciences, 2003,40(1):25-40.

[10] 唐志成, 夏才初, 宋英龍, 等. Grasselli節理峰值抗剪強度公式再探[J]. 巖石力學與工程學報, 2012,31(2):356-364.

[11] 唐志成, 夏才初, 宋英龍. 粗糙節理的峰值抗剪強度準則[J]. 巖土工程學報, 2013,35(3):571-577.

[12] Cundall P A, Hart R D. Analysis of block test no. 1 inelastic rock mass behavior: phase 2-a characterization of joint behavior(final report)[R]. Itasca Consulting Group, 1984.

[13] Bandis S C, Lumsden A C, Barton N R. Fundamentals of rock joint deformation[J]. International Journal of Rock Mechanics and Mining Sciences & Geomechanics Abstracts. Pergamon, 1983,20(6):249-268.

[14] Hong E S, Lee J S, Lee I M. Underestimation of roughness in rough rock joints[J]. International Journal for Numerical and Analytical Methods in Geomechanics, 2008,32(11):1 385-1 403.

[15] 高艷華, 吳順川, 王 賀, 等. 基于持續屈服節理模型的節理直剪數值試驗[J]. 中南大學學報 (自然科學版), 2016,(4):23.

[16] Itasca Consulting Group Inc. 3DEC (3 dimensional distinct element code) Constitutive Models[R]. Itasca Consulting Group Inc.，2013.

[17] Barton N, de Quadros E F. Joint aperture and roughness in the prediction of flow and groutability of rock masses[J]. International Journal of Rock Mechanics and Mining Sciences, 1997,34(3-4):252.e1-252.e14.

[18] Barton N. Modelling rock joint behavior from in situ block tests: implications for nuclear waste repository design[R]. USA: Terra Tek, Inc., Salt Lake City, UT, 1982.