節(jié)理網(wǎng)絡(luò)模擬的程序?qū)崿F(xiàn)及其在地下洞室中的動態(tài)響應(yīng)分析

賈 超，王志鵬，朱維申

（山東大學(xué) 土建與水利學(xué)院，濟(jì)南 250061）

1 引 言

隨著國內(nèi)外地下水電站、地下隧道、大型礦山巷道等的大規(guī)模興建，巖石力學(xué)的研究重點(diǎn)日益轉(zhuǎn)入地下。在地下工程的建設(shè)中，不可避免地要遇到巖體的開挖和處理問題，其中節(jié)理巖體是各種巖體工程和環(huán)境工程中經(jīng)常遇到的一種復(fù)雜巖體，節(jié)理裂隙巖體的穩(wěn)定性對整個地下工程巖體結(jié)構(gòu)的安全性具有至關(guān)重要的作用[1－4]。因此，正確地模擬節(jié)理巖體裂隙網(wǎng)絡(luò)對研究巖體節(jié)理裂隙的特征和分布規(guī)律具有重要的意義，同時為地下工程的安全穩(wěn)定性分析提供重要依據(jù)[5]。

節(jié)理巖體地下洞室在地應(yīng)力和地震載荷的作用下，失穩(wěn)破壞的主要形式有[6]：邊墻或拱頂?shù)拇竺娣e垮塌、大巖塊的塌落過度的塑性區(qū)產(chǎn)生和相應(yīng)的邊墻大位移、脆性巖體中高邊墻出現(xiàn)大深度的裂縫帶、突然的巖爆等。這些宏觀現(xiàn)象會對工程安全產(chǎn)生極其重要的影響，加之工程開挖所誘發(fā)的卸荷裂紋，具有平行于臨空面的定向性特征，并經(jīng)常會滯后發(fā)生。目前，國內(nèi)外對地下洞室的抗震研究還處于探索階段，在這方面還缺乏系統(tǒng)的研究。因此，研究地下洞室的抗震安全性十分重要。國內(nèi)學(xué)者隋斌等[6]研究了大型地下洞室群在地震荷載作用下的動態(tài)響應(yīng)模擬分析，并采用新的劈裂破壞判據(jù)對震后可能出現(xiàn)的劈裂損傷區(qū)范圍進(jìn)行了預(yù)測。張麗華等[7]采用動力離散元法和模型試驗的方法分析了大型巖體地下洞室群的動力響應(yīng)，認(rèn)為節(jié)理巖體開挖后在地震動力的作用下將加劇塊體的運(yùn)動，進(jìn)而影響洞室的整體穩(wěn)定性。李海波等[8]分析了地震荷載作用下埋深、洞室形狀等特征對地下巖體洞室位移特征的影響，并得出對工程設(shè)計有指導(dǎo)意義的結(jié)論。馬行東、劉華北等[9－10]研究了地震波入射方向?qū)Φ叵聨r體洞室動態(tài)響應(yīng)的影響。

隨著計算機(jī)技術(shù)的發(fā)展和普及以及各種數(shù)值分析方法在工程地質(zhì)領(lǐng)域的廣泛應(yīng)用，建立在統(tǒng)計學(xué)和概率論基礎(chǔ)上的結(jié)構(gòu)面網(wǎng)絡(luò)模擬技術(shù)逐步得到了深入發(fā)展。而基于結(jié)構(gòu)面網(wǎng)絡(luò)模擬結(jié)果，開展進(jìn)一步的數(shù)值計算分析，在巖體工程中具有良好的發(fā)展前景。本文首先編制了結(jié)構(gòu)面網(wǎng)絡(luò)模擬分析程序，并將程序生成的可執(zhí)行數(shù)據(jù)文件與離散元軟件相融合，最后以含節(jié)理巖體的地下洞室為例，分別計算了靜力和地震荷載作用下，地下洞室圍巖在不同時段的應(yīng)力及位移變化情況。計算方法與結(jié)果對含節(jié)理巖體地下洞室圍巖變形和破壞機(jī)制以及抗震工程設(shè)計具有指導(dǎo)意義，同時為更進(jìn)一步研究巖體結(jié)構(gòu)的性質(zhì)提供幫助。

2 相關(guān)基本原理

2.1 節(jié)理巖體裂隙網(wǎng)絡(luò)模擬程序原理

節(jié)理巖體裂隙網(wǎng)絡(luò)模擬的生成過程與現(xiàn)場實測統(tǒng)計過程恰好相反。現(xiàn)場實測統(tǒng)計過程是根據(jù)巖體結(jié)構(gòu)形式推求結(jié)構(gòu)面幾何參數(shù)的分布函數(shù)形式及建立結(jié)構(gòu)面的概率模型過程。而計算機(jī)模擬形式是上述過程的逆過程，即根據(jù)實測統(tǒng)計確立的結(jié)構(gòu)面幾何參數(shù)的概率模型，求服從這一模型的幾何圖形。

廣泛應(yīng)用于節(jié)理巖體裂隙網(wǎng)絡(luò)模擬中的蒙特卡羅(Monte-Carlo)方法是一種基于“隨機(jī)數(shù)”的計算方法。它常用于模擬過程中隨機(jī)數(shù)的產(chǎn)生，通過對隨機(jī)變量的隨機(jī)模擬和統(tǒng)計試驗來反推實際問題的近似解。研究表明，節(jié)理巖體結(jié)構(gòu)面的發(fā)育具有隨機(jī)性，各幾何特征參數(shù)（方位、規(guī)模、間距、隙寬等）服從一定的概率分布形式，因此，通過對節(jié)理巖體結(jié)構(gòu)面測量數(shù)據(jù)的分析可以建立各參數(shù)的統(tǒng)計概率模型[3]。模擬產(chǎn)生與這個概率模型相似或平行的隨機(jī)數(shù)，再通過抽樣統(tǒng)計求出服從這個分布規(guī)律的統(tǒng)計估算值，作為分析問題的數(shù)值近似解，從而模擬了與實際結(jié)構(gòu)面在統(tǒng)計特征上相同的結(jié)構(gòu)面網(wǎng)絡(luò)，且在計算機(jī)上再現(xiàn)節(jié)理巖體結(jié)構(gòu)面的分布形態(tài)。

2.2 動力離散元計算原理

在離散元中塊體的運(yùn)動一般采用牛頓第2定律來表示[11－12]，即：

進(jìn)行中心差分得：

則速度方程轉(zhuǎn)化為

引進(jìn)阻尼后，動力學(xué)基本方程為

在動力分析中阻尼需要再現(xiàn)對象在承受動力荷載作用下系統(tǒng)能量的損失。對于節(jié)理巖體這種材料，自然阻尼一般情況下是滯后的，與路徑相關(guān)，因此，也就很難用數(shù)值方法再現(xiàn)這種阻尼，所以阻尼參數(shù)的確定通常采用瑞利阻尼[6,11－12]。

3 程序?qū)崿F(xiàn)及與離散元的融合

3.1 程序的實現(xiàn)

本文應(yīng)用巖體結(jié)構(gòu)統(tǒng)計理論及研究成果，在前人研究的基礎(chǔ)上，編制了巖體節(jié)理裂隙網(wǎng)絡(luò)模擬分析程序，該程序主要功能及優(yōu)點(diǎn)如下：

①能夠根據(jù)需要實時選取參數(shù)的分布函數(shù),生成任意概率分布的節(jié)理裂隙形態(tài)。

②該程序可輸出多種格式的數(shù)據(jù)文件，便于用戶將數(shù)據(jù)結(jié)果結(jié)合其他軟件進(jìn)行下一步的分析計算。

③在程序調(diào)試通過的基礎(chǔ)上，運(yùn)用可視化的VB語言將網(wǎng)絡(luò)模擬程序進(jìn)行了擴(kuò)充與發(fā)展，得到了新的可視模擬程序，并且增加了友好的用戶界面。

④該程序集成了多種功能，成為一個比較完善的計算程序，并且用戶能夠容易地使用其操作界面，同時無需再去編譯數(shù)據(jù)文件可有效地減少錯誤。

該程序的技術(shù)路線如圖1所示：

圖1 節(jié)理網(wǎng)絡(luò)模擬流程圖Fig.1 Flowchart of joint network simulation

表1所列為某巖體三組結(jié)構(gòu)面的統(tǒng)計參數(shù)，應(yīng)用本文程序生成的結(jié)構(gòu)面網(wǎng)絡(luò)圖如圖2所示，生成的結(jié)構(gòu)面數(shù)據(jù)信息如圖3所示。

表1 巖體結(jié)構(gòu)面幾何參數(shù)Table 1 Rock structure surface geometric parameters

圖2 節(jié)理網(wǎng)絡(luò)模擬圖Fig.2 Diagram of joint network simulation

圖3 節(jié)理數(shù)據(jù)信息Fig.3 The data of joints

3.2 程序與離散元軟件UDEC的融合

UDEC是基于離散單元法理論的一款計算分析程序。離散單元法最早由Peter Cundall在1971年提出理論雛形，Cundall等在1980年開始又把這一方法思想拓展到研究顆粒狀物質(zhì)的微破裂、破裂擴(kuò)展、和顆粒流動問題。與有限元技術(shù)等通用連續(xù)力學(xué)方法相比較，屬于非連續(xù)力學(xué)方法范疇的 UDEC程序基于離散的角度來對待物理介質(zhì)，以最為樸素的思想分別描述介質(zhì)內(nèi)的連續(xù)性元素和非連續(xù)性元素。

該軟件特別適合于模擬節(jié)理巖石系統(tǒng)或者不連續(xù)塊體集合體系在靜力或動力荷載條件下的響應(yīng)問題[12]。同時UDEC軟件預(yù)留了生成節(jié)理的外部接口平臺，按照UDEC計算程序的要求，將隨機(jī)節(jié)理生成程序結(jié)果與UDEC程序進(jìn)行有機(jī)融合，可以進(jìn)行有效的應(yīng)用計算，使之更有利于如節(jié)理巖體地基、地震、地下結(jié)構(gòu)等問題的分析研究。本文基于結(jié)構(gòu)面網(wǎng)絡(luò)模擬結(jié)果，建立了含節(jié)理巖體的地下洞室?guī)r體力學(xué)分析模型，然后分別計算了靜力和地震荷載作用下，地下洞室圍巖在不同時段的應(yīng)力及位移變化情況。該方法為計算機(jī)節(jié)理網(wǎng)絡(luò)模擬結(jié)果與工程模型的融合計算提供了參考。

4 算例分析

4.1 計算模型

本文中計算模型為某含節(jié)理巖體的地下洞室，洞室平面尺寸為6 m×7.5 m。為體現(xiàn)結(jié)構(gòu)面模擬結(jié)果與工程模型的融合計算效果，由本文所編制的程序隨機(jī)生成了兩組節(jié)理并融合UDEC軟件，生成的模型圖以及數(shù)值分析中的驗算點(diǎn)位置如圖4所示。巖體及結(jié)構(gòu)面參數(shù)如表2、3所示。

4.2 可變形塊體及結(jié)構(gòu)面模型

在數(shù)值模擬計算中，地下洞室?guī)r體模型及結(jié)構(gòu)面模型均采用理想彈塑性模型，屈服準(zhǔn)則采用Mohr-Coulomb強(qiáng)度準(zhǔn)則，屈服函數(shù)如下：

圖4 計算模型Fig.4 Computational model

表2 巖塊變形和強(qiáng)度參數(shù)Table 2 Deformation and strength parameters of rock

表3 節(jié)理變形和強(qiáng)度參數(shù)Table 3 Deformation and strength parameters of joints

式中： N?=1 + sin? 1 - sin ?；其中，?為內(nèi)摩擦角；σ1、σ3分別為最大和最小主應(yīng)力；c為黏聚力；σt為抗拉強(qiáng)度。

當(dāng)巖體內(nèi)部某一點(diǎn)應(yīng)力滿足 fs＜0，就發(fā)生剪切屈服；當(dāng)巖體內(nèi)部某一點(diǎn)應(yīng)力滿足 ft＞0，發(fā)生張拉屈服。

4.3 邊界條件和地震荷載的確定

由于在動力分析中，計算模型的區(qū)域邊界有可能造成波的反射，給數(shù)值分析結(jié)果的準(zhǔn)確性造成一定影響[13]。因此，在本文模型問題域的外周邊界引入黏滯邊界用以消除波的反射，從而模擬有限的巖體；考慮到剪切波對洞室響應(yīng)影響較大，模擬的地震波為施加在模型底面的正弦剪切應(yīng)力波[8]，其速度時程表征為

最大振幅取1.25 m/s，卓越周期T為1 s，頻率為5 Hz，持續(xù)時間為2 s，峰值應(yīng)力為1.25 MPa。由于本文模型邊界為黏滯邊界，所以速度與加速度不能直接作用在模型邊界上，而要轉(zhuǎn)換成應(yīng)力作用在模型邊界上[14－15]，轉(zhuǎn)換公式為

式中：σn為施加的正應(yīng)力；σs為剪應(yīng)力；ρ為密度；Cp為傳入的p波波速；Cs為傳入的s波波速，vn為振動點(diǎn)的垂直速度分量；vs為振動點(diǎn)切向速度分量。

Cp與Cs計算公式分別為

對加速度時程進(jìn)行積分轉(zhuǎn)化為相應(yīng)的速度時程，然后由上式轉(zhuǎn)化為相應(yīng)的應(yīng)力時程：

本文阻尼系數(shù)采用瑞利阻尼[6,16]，其中質(zhì)量阻尼系數(shù)α=1.57，剛度阻尼系數(shù)β=1.59×10-3。

5 計算過程及結(jié)果分析

整個數(shù)值模擬分為3步。

①在開挖之前，對模型進(jìn)行固結(jié)，形成初始地應(yīng)力場并且達(dá)到平衡狀態(tài)，初始應(yīng)力狀態(tài)按各向同性估計為24 MPa（假定垂直荷載由覆蓋深度大約為800 m的巖層產(chǎn)生），固結(jié)后，初始地應(yīng)力場的主應(yīng)力分布圖如圖 5所示，初始最小主應(yīng)力為-26.50 MPa，最大主應(yīng)力為-19.64 MPa，圖中密集分布的部位不是應(yīng)力集中，而是差分網(wǎng)格密集的部位[7]。

圖5 洞室開挖前主應(yīng)力分布圖Fig.5 Principal stress distribution before excavation

②模擬在靜力條件下地下洞室的開挖過程。洞室開挖后位移主要發(fā)生在拱頂部位，位移矢量圖如圖6所示，拱頂A點(diǎn)在垂直方向的位移隨時間歷程曲線圖如圖7所示，最大位移為10.34 mm，模型最終循環(huán)至平衡狀態(tài)。洞室開挖引起了圍巖應(yīng)力重新分布，主應(yīng)力分布如圖8所示，最小主應(yīng)力為-53.1 MPa，最大主應(yīng)力為3.11 MPa，由計算結(jié)果可知洞室在開挖后處于穩(wěn)定狀態(tài)。

圖6 洞室開挖后位移矢量圖Fig.6 Displacement vector diagram after the cavern excavation

圖7 開挖后洞頂A在Y方向上位移隨時間歷程圖Fig.7 Time history responses of displacements of the point A in Y direction after excavation

圖8 洞室開挖后主應(yīng)力分布圖Fig.8 Principal stress distribution after excavation of cavern

③最后施加地震波，模擬洞室開挖完畢后在動力荷載作用下圍巖的位移變化及應(yīng)力分布情況。由于在動力分析中，計算模型的區(qū)域邊界有可能造成波的反射，給數(shù)值分析結(jié)果的準(zhǔn)確性造成一定影響[13]。因此，在本文模型問題域的外周邊界引入黏滯邊界用以消除波的反射；考慮到剪切波對洞室響應(yīng)影響較大，模擬的地震波為施加在模型底面的正弦剪切應(yīng)力波。施加地震作用0.02 s后洞室的位移矢量圖如圖9所示，洞頂A點(diǎn)在垂直方向上的最大位移為 9.12 mm；根據(jù)公路隧道施工技術(shù)規(guī)范(2006-12-25)[17]中規(guī)定：當(dāng)覆蓋層厚度大于300 m時隧道周邊允許相對位移值為1%～3%，相對位移為實測位移值與兩點(diǎn)間距離之比或洞頂實測值與隧道寬度之比。所以在施加地震作用0.02 s后洞室圍巖相對位移值為 1.52%，在允許范圍之內(nèi)，洞室此時是穩(wěn)定的。施加地震荷載后，計算結(jié)果顯示，洞室的圍巖應(yīng)力發(fā)生了重分布，地下洞室的主應(yīng)力圖如圖10所示，洞室最小主應(yīng)力為-53.16 MPa，最大主應(yīng)力為3.13 MPa。

圖9 地震荷載作用下0.02 s后洞室位移矢量圖Fig.9 The displacement vector diagram of cavern under the earthquake loads after 0.02 s

圖10 地震荷載作用下0.02 s后洞室主應(yīng)力分布圖Fig.10 Principal stress distribution under the earthquake loads after 0.02 s

施加地震作用0.04 s后，洞室位移矢量圖及拱頂 A點(diǎn)在垂直方向的位移隨時間歷程曲線圖如圖11、12所示，此時最大位移為242.2 mm，相對位移值為4%，大于允許相對位移值1%～3%，洞室在地震荷載的作用下已經(jīng)開始破壞。洞室的主應(yīng)力圖如圖13所示，洞室最小主應(yīng)力為-53.24 MPa，最大主應(yīng)力為3.31 MPa。

計算結(jié)果表明，在地震荷載作用下，巖體的結(jié)構(gòu)破壞更加突出。本文所示算例，地震時洞室頂部破壞的后果要遠(yuǎn)較靜力時嚴(yán)重，因此，進(jìn)行地下洞室的抗震設(shè)計時應(yīng)予以重視。

圖11 地震作用下洞頂A在Y方向上位移隨時間歷程圖Fig.11 Time history responses of displacements of the point A in Y direction under earthquake loads

圖12 地震荷載作用下0.04 s后洞室位移矢量圖Fig.12 The displacement vector diagram of cavern under earthquake loads after 0.04 s

圖13 地震荷載作用下0.04 s后洞室主應(yīng)力分布圖Fig.13 Principal stress distribution under earthquake loads after 0.04 s

6 結(jié) 語

基于結(jié)構(gòu)面網(wǎng)絡(luò)模擬結(jié)果，建立巖體力學(xué)分析模型，并開展進(jìn)一步的數(shù)值計算分析，在巖體工程中具有良好的發(fā)展前景。本文以含節(jié)理巖體的地下洞室為例，分別計算了靜力和地震荷載作用下，地下洞室圍巖在不同時段的應(yīng)力及位移變化情況。

計算結(jié)果表明：在地震荷載作用下，巖體的結(jié)構(gòu)破壞更加突出。本文所示算例，地震時洞室頂部破壞的后果要遠(yuǎn)較靜力時嚴(yán)重，因此，進(jìn)行地下洞室的抗震設(shè)計時應(yīng)予以重視。但地震動荷載作用下洞室的位移響應(yīng)分析是非常復(fù)雜的問題，本文基于簡單的動力荷載形式和計算模型來分析含節(jié)理巖體地下洞室的破壞問題，僅僅是一個初步的探討，但采用的研究思路及結(jié)果對含節(jié)理巖體地下洞室圍巖變形和破壞機(jī)制以及抗震工程設(shè)計具有指導(dǎo)意義。

[1]陶振寧. 巖石力學(xué)理論與實踐[M]. 北京: 水利電力出版社, 1989.

[2]GOODMAN R E. Introduction to rock mechanics[M].New York: John Wiley Sons, 1989.

[3]周維垣. 高等巖石力學(xué)[M]. 北京: 水利電力出版社,1990.

[4]黃國明. 節(jié)理巖體描述及其工程應(yīng)用[D]. 成都: 成都理工大學(xué), 1999.

[5]陳劍平. 巖體隨機(jī)不連續(xù)面三維網(wǎng)絡(luò)數(shù)值模擬技術(shù)[J].巖土工程學(xué)報, 2001, 23(4): 397－402.CHEN Jian-ping. 3-D net work numerical modeling technique for random discontinuities of rock mass[J].Chinese Journal of Geotechnical Engineering, 2001,23(4): 397－402.

[6]隋斌, 朱維申, 李曉靜. 地震荷載作用下大型地下洞室群的動態(tài)響應(yīng)模擬[J]. 巖土工程學(xué)報, 2008, 30(12):1877－1882.SUI Bin, ZHU Wei-shen, LI Xiao-jing. Simulation of dynamic response of large underground opening complex under seismic loads[J]. Chinese Journal of Geotechnical Engineering, 2008, 30(12): 1877－1882.

[7]張麗華, 陶連金. 節(jié)理巖體地下洞室群的地震動力響應(yīng)分析[J]. 世界地震工程, 2002, 18(2): 158－162.ZHANG LI-hua, TAO Lian-jin. Dynamic response analysis of large underground excavations in jointed rock[J]. World Earthquake Engineering, 2002, 18(2):158－162.

[8]李海波, 馬行東, 李俊如, 等. 地震荷載作用下地下巖體洞室位移特征的影響因素分析[J]. 巖土工程學(xué)報,2006, 28(3): 358－362.LI Hai-bo, MA Xing-dong, LI Jun-ru, et al. Study on influence factors of rock cavern displacement under earthquake[J]. Chinese Journal of Geotechnical Engineering, 2006, 28(3): 358－362.

[9]馬行東, 李海波. 地震波入射方向?qū)Φ叵聨r體洞室動態(tài)響應(yīng)的初步分析[J]. 水力發(fā)電, 2007, 33(1): 23－25.MA Xing-dong, LI Hai-bo. Primary analysis of different incidence propagation on dynamic response of underground rock cavern under earthquake[J]. Water Power, 2007, 33(1): 23－25.

[10]LIU Hua-bei, Song Er-xiang. Seismic response of large underground structures in liquefiable soils subjected to horizontal and vertical earthquake excitations[J].Computers and Geotechnics, 2005, 32(4): 223－244.

[11]LIN M, KICKER D, DAMJANC B, et al. Mechanical degradation of emplacement drifts at Yucca mountain—A modeling case study—Part I: Nonlithophysal rock[J].International Journal of Rock Mechanics & Mining Sciences, 2007, (44): 351－367.

[12]王泳嘉, 邢紀(jì)波. 離散元法及其在巖土力學(xué)中的應(yīng)用[M]. 遼寧: 東北工學(xué)院出版社, 1991.

[13]CHEN Jian-yun, JING Li, ZHOU Jing. Seismic analysis of large 3-D underground caverns based on high performance recursive method[J]. International Journal of Solids and Structures, 2004, 41(11－12): 3081－3094.

[14]夏祥, 李俊如, 李海波, 等. 爆破荷載作用下巖體振動特征的數(shù)值模擬[J]. 巖土力學(xué), 2005, 26(1): 50－56.XIA Xiang, LI Jun-ru, LI Hai-bo, et al. Blasting vibration characteristics under the load of the numerical simulation of rock[J]. Rock and Soil Mechanics, 2005, 26(1): 50－56.

[15]呂濤. 地下水與地震作用下節(jié)理巖體地下洞室圍巖力學(xué)響應(yīng)研究[D]. 北京: 北京工業(yè)大學(xué), 2005.

[16]王濤, 熊將, 郭武祥, 等. 地震荷載作用下地下廠房圍巖穩(wěn)定的離散元計算方法研究[J]. 長江科學(xué)院院報,2005, 26(12): 58－62.WANG Tao, XIONG Jiang, GUO Wu-xiang, et al.Research on stability of rock mass around underground powerhouse cavern by discrete element method under earthquake[J]. Journal of Yangtze River Scientific Research Institute, 2005, 26(12): 58－62.

[17]交通部重慶公路科學(xué)研究所. JTJ042－94 公路隧道施工技術(shù)規(guī)范[S]. 北京: 人民交通出版社, 2006.