深隧道圍巖分區(qū)破裂的數(shù)學(xué)模擬

戚承志，錢七虎，王明洋，陳劍杰

（1. 北京建筑工程學(xué)院 工程結(jié)構(gòu)與新材料北京市高校工程研究中心，北京 100044；2. 解放軍理工大學(xué) 工程兵工程學(xué)院，南京 210007；3. 西北核技術(shù)研究院，西安 710024）

1 引 言

自俄羅斯學(xué)者發(fā)現(xiàn)深部隧道圍巖的分區(qū)破裂現(xiàn)象以來[1]30多年過去了，俄羅斯學(xué)者[2－8]、中國學(xué)者[9－13]進(jìn)行了一系列理論和試驗(yàn)研究，這些研究從不同的角度分析了分區(qū)破裂現(xiàn)象的機(jī)制，促進(jìn)了對于這一現(xiàn)象的理解。筆者在前期工作中[9,13－14]，利用連續(xù)相變理論成功地模擬了深部隧道圍巖分區(qū)破裂化現(xiàn)象的時(shí)空結(jié)構(gòu)，也解釋了與這一現(xiàn)象密切相關(guān)的受損巖石試樣臨界后變形增量的逆轉(zhuǎn)現(xiàn)象。但上述工作的不足是，沒有把解析解與巖石的臨界后變形特征參數(shù)（例如巖石臨界后變形的軟化模量等）聯(lián)系起來。20世紀(jì)80年代之后發(fā)展起來的梯度塑性理論[15－22]在解釋固體的變形局部化和有規(guī)律的變形模式空間分布方面取得了很大的進(jìn)展，這些理論采用應(yīng)變梯度或者其他內(nèi)變量梯度作為額外的模型變量，這些梯度變量進(jìn)入到介質(zhì)的平衡方程中或者內(nèi)變量演化方程中，得到了不可逆變形的局部化解，或者空間分布解。在文獻(xiàn)[23－24]中作者利用了內(nèi)變量梯度理論，模擬了深部隧道圍巖的分區(qū)破裂現(xiàn)象。但這一理論的缺點(diǎn)是梯度項(xiàng)沒有進(jìn)入到平衡方程中去，因此，不能夠通過平衡方程把隧道圍巖彈性區(qū)和塑性區(qū)的受力和變形狀態(tài)統(tǒng)一解出來。本文將致力于把梯度項(xiàng)融進(jìn)平衡方程中去，以便建立自洽的深部隧道圍巖的分區(qū)破裂模型。

2 介質(zhì)變形的非局部變量

位移矢量u的分量為ui，應(yīng)變張量ε的分量為εij。

其中

梯度理論可以分成兩類：應(yīng)變梯度和內(nèi)變量梯度理論。這兩種理論的根本區(qū)別在于，應(yīng)變梯度理論中，應(yīng)變梯度被選為狀態(tài)變量，與應(yīng)變梯度共軛的高階應(yīng)力進(jìn)入到平衡方程中去；而內(nèi)變量梯度與某些耗散熱力學(xué)力相對應(yīng)，這些熱力學(xué)力可以進(jìn)入到內(nèi)變量的演化方程中去，但不出現(xiàn)在平衡方程或者動(dòng)量守恒方程中。因此，內(nèi)變量梯度理論只修正本構(gòu)關(guān)系的描述，而幾何方程與平衡方程沒變。從熱力學(xué)的角度來看，內(nèi)變量梯度模型只修正了自由能函數(shù)和耗散勢函數(shù)，而應(yīng)變梯度理論需要擴(kuò)展內(nèi)力功與外力功的表達(dá)式。

在本文中將利用內(nèi)變量模型來研究巖體的變形。

在描述連續(xù)介質(zhì)的非局部或者梯度本構(gòu)關(guān)系時(shí)通常引入下列形式的單位體積的內(nèi)能[21]：

3 虛功原理

對于給定介質(zhì)的一個(gè)子區(qū)域V，作用在V內(nèi)材料上的外力功等于花費(fèi)于V內(nèi)材料上的內(nèi)力功。外部功是由宏觀體力、作用于介質(zhì)表面Γ上的宏觀面力產(chǎn)生：

式中：δ為變分算符；δu為虛位移。

為了描述介質(zhì)的不可逆變形，內(nèi)變量q可以取為介質(zhì)的全應(yīng)變ε、累積有效塑性應(yīng)變=而非局部變量Q(q)取為應(yīng)變梯度?ε。

因?yàn)?/p>

把 εe= ε- εp代入式（6），得

虛功原理可以表示為

將式（5）、（9）代入式（10），得

式（12）為介質(zhì)的平衡方程，其中考慮了梯度效應(yīng)；式（13）為介質(zhì)的宏觀邊界條件；式（15）為高階的應(yīng)力邊界條件；而式（14）定義了介質(zhì)的屈服條件。對于該方程兩邊取絕對值得

式中：Ωij為材料隨動(dòng)硬化參數(shù)，當(dāng)材料各向同性硬化時(shí)，有 Ωij=0，此時(shí)式（16）變?yōu)?/p>

所以B具有材料塑性限的物理意義。

4 本構(gòu)關(guān)系的熱力學(xué)限制

介質(zhì)的不可逆變形的出現(xiàn)引起了介質(zhì)內(nèi)部能量的耗散，所以介質(zhì)的本構(gòu)關(guān)系必須符合熱力學(xué)定律的限制。

熱力學(xué)第一定律可以表示為

如果要求能量守恒方程對于每一點(diǎn)都成立，那么必須在等式的右邊增加一項(xiàng)非局部能量剩余項(xiàng)R，來表示位于該點(diǎn)的介質(zhì)單元由于相互作用得到的來自介質(zhì)中其他單元的能量：

對于非局部能量剩余項(xiàng)R，Polizzotto[21－22]引入下列隔離條件來表示體積V內(nèi)沒有任何非局部能量流出該體積之外，即

非局部能量剩余項(xiàng)R如同Abu-Alrub等[20]所指出的那樣，除了與塑性梯度效應(yīng)外，還可能與表面效應(yīng)和界面效應(yīng)有關(guān)。如果不考慮表面和界面效應(yīng)、在彈性區(qū)域以及不考慮梯度效應(yīng)時(shí)，R=0。在現(xiàn)有的很多梯度塑性模型中忽略了非局部能量剩余項(xiàng)R，此處將忽略非局部能量剩余項(xiàng)R。

引入自由能函數(shù)ψ =U-Ts ，對于等溫過程，即T˙=0，由式（19）可得

自由能密度函數(shù)ψ可以表示為

對于式（22）關(guān)于時(shí)間求導(dǎo)得

因?yàn)橛墒剑?）有

把上式代入到式（23）中得

把式（24）代入到式（21）得

由式（25）得

可以引入下列廣義力

這樣式（27）可以寫為

式（25）可簡化為

所以為了滿足式（31）必須有如下關(guān)系：

為了簡化分析，在此忽略掉隨動(dòng)硬化效應(yīng)，即設(shè)塑性應(yīng)變張量不進(jìn)入到塑性自由能函數(shù)中去。設(shè)塑性自由能函數(shù)為

式中：H為硬化（軟化）模量，當(dāng)H＞0時(shí)材料具有硬化效應(yīng)，當(dāng)H＜0時(shí)材料具有軟化效應(yīng)。

這樣除了得到廣義虎克定律外還會得到下列廣義力：

把式（34）代入到式（16），得

5 考慮應(yīng)變梯度時(shí)問題的支配方程

對于平衡方程式（11），當(dāng)忽略體力時(shí)變?yōu)?/p>

對于深部圓形隧道圍巖，利用極坐標(biāo)分析較為方便。取圓柱坐標(biāo)的徑向、切向坐標(biāo)分別為r和θ。高階應(yīng)力導(dǎo)數(shù)的展開式較復(fù)雜，為了簡化求解，只保留項(xiàng)，更一般的情況將在以后研究。這樣在極坐標(biāo)里式（36）變?yōu)?/p>

5.1 彈性情況下的支配方程

幾何方程為

對于Tijk，在軸對稱情況下可以取為下列簡單的形式

其中c1＞0為系數(shù)。這樣平衡方程（37）變?yōu)?/p>

式（41）的解為

式中：A、B 、C 、 D 為待定常數(shù)； I1、 K1分別為第一類和第二類變型Bessel函數(shù)。

因?yàn)樵趓→∞時(shí)ur應(yīng)該有限，所以B=0、C= 0，所以式（26）最終變?yōu)?/p>

所以徑向變形εr和環(huán)向變形εθ為

而應(yīng)力為

5.2 具有下降段的彈塑性情況下的支配方程

具有軟化段的應(yīng)力-應(yīng)變曲線如圖1所示[26]。

圖1 具有軟化段的應(yīng)力-應(yīng)變曲線Fig.1 Stress-strain curve with softening

下降段用應(yīng)力可以表示為

如果引入側(cè)向變形系數(shù)

那么式（46）變?yōu)?/p>

對于彈性變形如果為了簡單起見取泊松比為μ= 0.5，那么彈性變形可以表示為

而總的變形為

式（51）可以變?yōu)?/p>

其中

形式上式（52）可以變?yōu)?/p>

對于彈塑性問題，通常假定式（39）成立[19]。利用公式（53）和方程（31），用σ1和σ2分別代替（31）中的σr和σθ，最終得到具有圓形橫截面的深部隧道圍巖的平衡方程

在隧道壁 r=r0處由式（13）、（15）可得邊界條件

由于邊界條件展開式復(fù)雜，此處不給出展開式。

式（54）～（56）的解析解較難得到?，F(xiàn)在考慮一種特殊情況 l1- k2= （k1- l2），如果忽略體積變形εr+εθ=0。那么

這樣式（54）變?yōu)?/p>

其中 ξ2= k1c1。

式（58）的解為

圍巖的應(yīng)變?yōu)?/p>

圍巖的應(yīng)力為

式（59）～（61）中：J0和J1分別為第一類 0階和1階Bessel函數(shù)；Y0和Y1分別為第二類0階和1階Bessel函數(shù)。

由式（41）和式（58）的比較可見，由于巖體本構(gòu)關(guān)系中出現(xiàn)了下降段，使得式（41）中的算符(?2- β2)中常數(shù)β前面的“-”符號變?yōu)槭剑?8）中算符 (ξ2+?2)中常數(shù)ξ前面的“+”符號，從而使得方程的解由單調(diào)性的函數(shù)式（42）轉(zhuǎn)變?yōu)槭剑?9）的具有準(zhǔn)周期振蕩的解。因此，可見巖體本構(gòu)關(guān)系中下降段的存在是引起分區(qū)破裂的必要條件。

由位移、應(yīng)力和應(yīng)變表達(dá)式（59）～（61）可見，這3組解的右邊前兩項(xiàng)對應(yīng)著彈塑性理論的解，而后面的Bessel函數(shù)項(xiàng)代表著對于彈塑性解的修正，這種修正隨著不可逆變形的增加而增加。這里的位移、應(yīng)變和應(yīng)力表達(dá)式與筆者在文獻(xiàn)[9, 13]中得到的解的差別在于，此處出現(xiàn)了彈塑性理論解，是在連續(xù)介質(zhì)力學(xué)框架內(nèi)的自洽解。

位移、應(yīng)變和應(yīng)力表達(dá)式具有隨距隧道壁距離增加幅值減小的準(zhǔn)周期性，能夠反映分區(qū)破裂現(xiàn)象的位移、應(yīng)變和應(yīng)力變化情況。文獻(xiàn)[9, 13]的研究表明，選取適當(dāng)?shù)腂essel函數(shù)解的常數(shù)可以模擬試驗(yàn)和現(xiàn)場得到的破裂區(qū)隨距隧道壁距離的分布情況，此處適當(dāng)選取式（59）～（61）中Bessel函數(shù)前面的系數(shù)可以模擬試驗(yàn)和現(xiàn)場得到的破裂區(qū)隨距隧道壁距離的分布情況。將在后續(xù)研究中研究式（59）～（61）的解。

5.3 忽略彈性變形的塑性情況下的支配方程

對于巖體的極限后變形，在此采用具有下降段的Tresca 本構(gòu)模型。為了簡化計(jì)算，在此采用剛塑性模型，即忽略彈性變形。

對于深部圍巖的分區(qū)破裂化，一般認(rèn)為，在地應(yīng)力達(dá)到巖體的單軸抗壓強(qiáng)度后開始出現(xiàn)分區(qū)破裂化現(xiàn)象，此時(shí)按照Tresca 屈服條件，兩個(gè)滑移面上的剪應(yīng)力T和T1達(dá)到極限狀態(tài)：

式中：τc為巖體的抗剪強(qiáng)度。

在極限后區(qū)域，巖體中的剪力將下降，兩個(gè)剪應(yīng)力隨著剪應(yīng)變的下降情況如圖2所示?？梢员硎緸?/p>

圖2 忽略彈性變形時(shí)具有軟化段的應(yīng)力-應(yīng)變曲線Fig.2 Stress-strain curves with softening without consideration of elastic deformations

對于式（63）第2式，σz等于巖體的原始地應(yīng)力：σz=σ∞，而 εz=0，所以式（63）第2式變?yōu)棣襯= σ∞- σc- M εr，dσr=- Mdεr。

式（64）的解為[27]

式中：η=kr， k2=Mc1，其他符號含義詳見文獻(xiàn)[27]。

將式（65）代入到應(yīng)變和應(yīng)力的表達(dá)式中，可以得到應(yīng)變和應(yīng)力的表達(dá)式。

與式（41）中的算符 (?2- β2)中常數(shù)β前面的“-”符號不同，式（58）中的算符 (ξ2+?2)中常數(shù)ξ前面的符號為“+”，從而使得方程（64）的解式（65）具有隨距隧道壁距離增加幅值減小的準(zhǔn)周期性，能夠反映分區(qū)破裂現(xiàn)象。方程（64）的解將在后續(xù)工作中進(jìn)行數(shù)值研究。

6 結(jié) 論

（1）本文利用塑性梯度理論研究了深部隧道圍巖的分區(qū)破裂現(xiàn)象。作為額外的狀態(tài)變量，在此引入應(yīng)變梯度這一新變量。利用虛功原理得到了巖體的平衡方程、邊界條件和流動(dòng)準(zhǔn)則。利用非局部的Clausius-Duhem不等式獲得了巖體的本構(gòu)方程。

（2）對于圓形深部隧道，由上述理論得到了彈性變形情況下、具有下降段的彈塑性變形情況下和不考慮彈性變形的塑性變形情況下深部隧道圍巖的支配方程。

（3）對于隧道圍巖的臨界后課題來講，應(yīng)變梯度的引進(jìn)能夠描述巖體處于臨界后狀態(tài)的自組織行為。

（4）巖體應(yīng)力-應(yīng)變曲線下降段的存在，是引起隧道圍巖出現(xiàn)隨距隧道壁距離增加位移、應(yīng)變和應(yīng)力準(zhǔn)周期性變化的必要條件。

[1]SHEMYAKIN E I, FISENKO G L, KURLENYA M V,et al. Zonal disintegration of rocks around underground workings. Part 1: Data of in situ observations[J]. Journal of Mining Science, 1986, 22(3): 157－168.

[2]SHEMYAKIN E I, FISENKO G L, KURLENYA M V,et al. Zonal disintegration of rocks around underground workings. Part II: Rock fracture simulated in equivalent materials[J]. Journal of Mining Science, 1986, 22(4):223－232.

[3]SHEMYAKIN E I, FISENKO G L, KURLENYA M V,et al. Zonal disintegration of rocks around underground mines. Part III: Theoretical concepts[J]. Journal of Mining Science, 1987, 23(1): 1－5.

[4]SHEMYAKIN E I, FISENKO G L, KURLENYA M V,et al. Zonal disintegration of rocks around underground workings. Part IV: Practical applications[J]. Journal of Mining Science, 1989, 25(4): 297－302.

[5]ODINTSEV V N. Mechanism of the zonal disintegration of a rock mass in the vicinity of deep-level workings[J].Journal of Mining Science, 1994, 30(4): 334－343.

[6]CHANYSHEV A I. To the problem of deformable medium failure. Part I: Basic equations[J]. Journal of Mining Science, 2001, 37(3): 273－288.

[7]CHANYSHEV A I. To the problem of deformable medium failure. Part II: Discussion of the results of analytical solutions[J]. Journal of Mining Science, 2001,37(4): 392－400.

[8]GUZEV M A, POROSHIN A A. Non-Euclidean model of the zonal disintegration of rocks around an underground working[J]. Journal of Applied Mechanics and Technical Physics, 2001, 42(1): 131－139.

[9]QI Cheng-zhi, QIAN Qi-hu, WANG Ming-yang.Evolution of the deformation and fracturing in rock masses near deep-level tunnels[J]. Journal of Mining Science, 2009, 45(2): 112－119.

[10]潘一山, 李英杰, 唐鑫, 等. 巖石分區(qū)破裂化現(xiàn)象研究[J]. 巖石力學(xué)與工程學(xué)報(bào), 2007, 26(增刊1): 3335－3341.PAN Yi-shan, LI Ying-jie, TANG Xin, et al. Study of zonal disintegration of rock[J]. Chinese Journal of Rock Mechanics and Engineering, 2007, 26(Supp.1): 3335－3341.

[11]賀永年, 蔣斌松, 韓立軍, 等. 深部巷道圍巖間隔性區(qū)域斷裂研究[J]. 中國礦業(yè)大學(xué)學(xué)報(bào), 2008, 37(3): 300－304.HE Yong-nian, JIANG Bin-song, HAN Li-jun, et al.Study of intermittent zonal fracturing of surrounding rock in deep roadways[J]. Journal of China University of Mining & Technology, 2008, 37(3): 300－304.

[12]顧金才, 顧雷雨, 陳安敏, 等. 深部開挖洞室圍巖分層斷裂破壞機(jī)制模型試驗(yàn)研究[J]. 巖石力學(xué)與工程學(xué)報(bào),2008, 27(3): 433－438.GU Jin-cai, GU Lei-yu, CHEN An-min, et al. Model test study of mechanism of layered fracture within surrounding rock of tunnels in deep stratum[J]. Chinese Journal of Rock Mechanics and Engineering, 2008,27(3): 433－438.

[13]WANG Ming-yang, QI Cheng-zhi, QIAN Qi-hu, et al.Physical modeling of the deformation increment sign change effect of rock sample under compression[J].Journal of Mining Science, 2010, 46(4): 21－28.

[14]戚承志, 錢七虎. 巖體動(dòng)力變形破壞的基本問題[M].北京: 科學(xué)出版社, 2009.

[15]AIFANTIS E C. On the microstructural origin of certain inelastic models[J]. Journal of Engineering Materials and Technology, 1984, 106(2): 326－330.

[16]VARDOULAKIS I, AIFANTIS E C. A gradient flow theory of plasticity for granular materials[J]. Acta Mechanics, 1991, 87(3－4): 197－217.

[17]AIFANTIS E C. Pattern formation in plasticity[J].International Journal of Engineering Science, 1995,33(15): 2161－2178.

[18]FLECK N A, HUTCHINSON J W. Strain gradient plasticity[C]//Advances in Applied Mechanics. New York:Academic Press, 1997: 295－361.

[19]CHAMBON R, CAILERIE D, HASAN N E.One-dimensional localization studied with a second grade model[J]. European Journal of Mechanics/A: Solids,1998, 17(4): 637－656.

[20]ABU AL-RUB R K, VOYIADJIS G Z. A physically based gradient plasticity theory[J]. International Journal of Plasticity, 2006, 22(3): 654－684.

[21]POLIZZOTTO C. Unified thermodynamic framework for nonlocal/gradient continuum theories[J]. European Journal of Mechanics/A: Solids, 2003, 22(5): 651－668.

[22]POLIZZOTTO C. Interfacial energy effects within the framework of strain gradient plasticity[J]. International Journal of Solids and Structures, 2009, 46(7－8): 1685－1694.

[23]QI Cheng-zhi, QIAN Qi-hu, WANG Ming-yang, et al.Mechanical modeling of zonal disintegration of rock mass near deep level tunnels[C]//Proceedings of the 3rd International Conference on Contemporary Problems in Architecture and Construction. Beijing: [s. n.], 2011: 114－119.

[24]WANG Ming-yang, QI Cheng-zhi, QIAN Qi-hu, et al.One plastic gradient model of zonal disintegration of rock mass near deep level tunnels[J]. Journal of Mining Science, 2012, 48(1): 46－54.

[25]ILYJUSHIN A A. Plasticity[M]. Moscow & Leningrad:OGIZ Press, 1948.

[26]BAKLASHOV I V, KARTOZIA B A. Mechanics of underground buildings and structures[M]. Moscow:Nedra Press, 1984.

[27]ZAITSEV V F, POLYANIN A D. Hand book for ordinary differential equations[M]. Moscow: Physmatlit Press, 2001.