賈 超，王志鵬，朱維申

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

1 引 言

隨著國內(nèi)外地下水電站、地下隧道、大型礦山巷道等的大規(guī)模興建，巖石力學(xué)的研究重點(diǎn)日益轉(zhuǎn)入地下。在地下工程的建設(shè)中，不可避免地要遇到巖體的開挖和處理問題，其中節(jié)理巖體是各種巖體工程和環(huán)境工程中經(jīng)常遇到的一種復(fù)雜巖體，節(jié)理裂隙巖體的穩(wěn)定性對整個(gè)地下工程巖體結(jié)構(gòu)的安全性具有至關(guān)重要的作用[1－4]。因此，正確地模擬節(jié)理巖體裂隙網(wǎng)絡(luò)對研究巖體節(jié)理裂隙的特征和分布規(guī)律具有重要的意義，同時(shí)為地下工程的安全穩(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àn)的方法分析了大型巖體地下洞室群的動力響應(yīng)，認(rèn)為節(jié)理巖體開挖后在地震動力的作用下將加劇塊體的運(yùn)動，進(jìn)而影響洞室的整體穩(wěn)定性。李海波等[8]分析了地震荷載作用下埋深、洞室形狀等特征對地下巖體洞室位移特征的影響，并得出對工程設(shè)計(jì)有指導(dǎo)意義的結(jié)論。馬行東、劉華北等[9－10]研究了地震波入射方向?qū)Φ叵聨r體洞室動態(tài)響應(yīng)的影響。

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

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

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

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

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

2.2 動力離散元計(jì)算原理

在離散元中塊體的運(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 程序的實(shí)現(xiàn)

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

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

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

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

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

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

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

表1所列為某巖體三組結(jié)構(gòu)面的統(tǒng)計(jì)參數(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是基于離散單元法理論的一款計(jì)算分析程序。離散單元法最早由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]。同時(shí)UDEC軟件預(yù)留了生成節(jié)理的外部接口平臺，按照UDEC計(jì)算程序的要求，將隨機(jī)節(jié)理生成程序結(jié)果與UDEC程序進(jìn)行有機(jī)融合，可以進(jìn)行有效的應(yīng)用計(jì)算，使之更有利于如節(jié)理巖體地基、地震、地下結(jié)構(gòu)等問題的分析研究。本文基于結(jié)構(gòu)面網(wǎng)絡(luò)模擬結(jié)果，建立了含節(jié)理巖體的地下洞室?guī)r體力學(xué)分析模型，然后分別計(jì)算了靜力和地震荷載作用下，地下洞室圍巖在不同時(shí)段的應(yīng)力及位移變化情況。該方法為計(jì)算機(jī)節(jié)理網(wǎng)絡(luò)模擬結(jié)果與工程模型的融合計(jì)算提供了參考。

4 算例分析

4.1 計(jì)算模型

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

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

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

圖4 計(jì)算模型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 邊界條件和地震荷載的確定

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

最大振幅取1.25 m/s，卓越周期T為1 s，頻率為5 Hz，持續(xù)時(shí)間為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ì)算公式分別為

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

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

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

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

①在開挖之前，對模型進(jìn)行固結(jié)，形成初始地應(yīng)力場并且達(dá)到平衡狀態(tài)，初始應(yīng)力狀態(tài)按各向同性估計(jì)為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)在垂直方向的位移隨時(shí)間歷程曲線圖如圖7所示，最大位移為10.34 mm，模型最終循環(huán)至平衡狀態(tài)。洞室開挖引起了圍巖應(yīng)力重新分布，主應(yīng)力分布如圖8所示，最小主應(yīng)力為-53.1 MPa，最大主應(yīng)力為3.11 MPa，由計(jì)算結(jié)果可知洞室在開挖后處于穩(wěn)定狀態(tài)。

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

圖7 開挖后洞頂A在Y方向上位移隨時(shí)間歷程圖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)力分布情況。由于在動力分析中，計(jì)算模型的區(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時(shí)隧道周邊允許相對位移值為1%～3%，相對位移為實(shí)測位移值與兩點(diǎn)間距離之比或洞頂實(shí)測值與隧道寬度之比。所以在施加地震作用0.02 s后洞室圍巖相對位移值為 1.52%，在允許范圍之內(nèi)，洞室此時(shí)是穩(wěn)定的。施加地震荷載后，計(jì)算結(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)在垂直方向的位移隨時(shí)間歷程曲線圖如圖11、12所示，此時(shí)最大位移為242.2 mm，相對位移值為4%，大于允許相對位移值1%～3%，洞室在地震荷載的作用下已經(jīng)開始破壞。洞室的主應(yīng)力圖如圖13所示，洞室最小主應(yīng)力為-53.24 MPa，最大主應(yīng)力為3.31 MPa。

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

圖11 地震作用下洞頂A在Y方向上位移隨時(shí)間歷程圖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ù)值計(jì)算分析，在巖體工程中具有良好的發(fā)展前景。本文以含節(jié)理巖體的地下洞室為例，分別計(jì)算了靜力和地震荷載作用下，地下洞室圍巖在不同時(shí)段的應(yīng)力及位移變化情況。

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

[1]陶振寧. 巖石力學(xué)理論與實(shí)踐[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é)報(bào), 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é)報(bào), 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é)報(bào),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ì)算方法研究[J]. 長江科學(xué)院院報(bào),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.