陳俊宇，佘成學(xué)

(武漢大學(xué) 水資源與水電工程科學(xué)國(guó)家重點(diǎn)實(shí)驗(yàn)室, 湖北 武漢 430072)

在復(fù)雜的地質(zhì)條件下，天然巖體受到各種構(gòu)造作用，其內(nèi)部會(huì)產(chǎn)生大量的貫通節(jié)理。節(jié)理的存在會(huì)影響巖體的力學(xué)行為，同時(shí)降低巖體的強(qiáng)度[1-3]。當(dāng)節(jié)理巖體的應(yīng)力達(dá)到強(qiáng)度極限后，隨著變形的增加，其強(qiáng)度會(huì)逐漸降低到一個(gè)較低的水平，此即節(jié)理巖體破壞后的應(yīng)變軟化現(xiàn)象[4]。其中一些節(jié)理巖體的軟化速率非常大，其軟化過程在應(yīng)力-應(yīng)變曲線中表現(xiàn)為一個(gè)很陡的峰后應(yīng)力下降段。對(duì)于此類節(jié)理巖體，可將其破壞后的軟化現(xiàn)象做近似處理，考慮為脆性破壞。

觀察節(jié)理巖體破壞后的物理狀態(tài)，發(fā)現(xiàn)其破壞模式并不唯一，可能因巖塊破壞導(dǎo)致整體破碎；也可能僅僅沿著節(jié)理面變形破壞[5]。在以往的研究中，鄭宏等[6]通過拓展經(jīng)典塑性理論，提出了脆塑性的巖塊和單結(jié)構(gòu)面的應(yīng)力跌落方式，但對(duì)于包含節(jié)理面的脆性巖體，并未在其變形破壞的模擬方法上展開研究。此外，朱澤奇等[7]建立了改進(jìn)的層狀巖體遍布節(jié)理模型。金華等[8]建立了含多組貫穿節(jié)理的擴(kuò)展遍布節(jié)理巖體模型，這些模型是基于理想彈塑性理論提出，不能模擬節(jié)理巖體的軟化破壞現(xiàn)象。朱道建等[9]引入軟化參數(shù)來(lái)描述節(jié)理巖體的軟化特性曲線，但在軟化曲線的確定方式上缺乏理論依據(jù)。同時(shí)，上述的多種節(jié)理巖體模型對(duì)材料屈服后的應(yīng)力修正方式也過于繁瑣復(fù)雜，其理論有待進(jìn)一步完善。

本文將節(jié)理巖體概化為一種復(fù)合材料[10]，分析節(jié)理巖體的變形破壞模式，基于彈塑性理論提出一種模擬其脆性破壞全過程的有限元方法，并利用ABAQUS編寫計(jì)算子程序。利用該模型模擬青松電站引水隧洞的開挖，證明所建模型可以運(yùn)用于實(shí)際工程，能較好反映脆性節(jié)理巖體的變形與破壞；用該方法可以正確模擬隧洞開挖后圍巖的屈服破壞情況，驗(yàn)證了薄弱節(jié)理面對(duì)節(jié)理巖體的屈服破壞的重大影響。

1 節(jié)理巖體的脆性破壞模型

1.1 對(duì)節(jié)理巖體復(fù)合材料的基本假定

節(jié)理巖體由巖塊和無(wú)厚度的節(jié)理面組成。巖塊為各向同性體，應(yīng)力-應(yīng)變關(guān)系在各個(gè)坐標(biāo)系中都相同；節(jié)理有對(duì)應(yīng)的局部坐標(biāo)系n-r-s，n為節(jié)理面外法線方向，r、s在面內(nèi)相互正交，分別沿節(jié)理的走向和傾向。面上有垂直于面的法向應(yīng)力σn，和在面內(nèi)相互正交的兩個(gè)剪切應(yīng)力τnr、τns。以層狀節(jié)理巖體為例，巖體與節(jié)理面的關(guān)系見圖1。

圖1節(jié)理巖體與節(jié)理面的關(guān)系

假定巖塊和節(jié)理面為串聯(lián)關(guān)系，則二者在整體坐標(biāo)系中承擔(dān)的荷載{σ}相等，節(jié)理巖體的總應(yīng)變{ε}等于巖塊與節(jié)理面的應(yīng)變之和[11]，即

{σ}={σR}={σJ}

(1)

{ε}={εR}={εJ}

(2)

式中：{σR}、{εR}分別為巖塊在整體坐標(biāo)系下的應(yīng)力和應(yīng)變；{σJ}、{εJ}分別為節(jié)理在整體坐標(biāo)系下的應(yīng)力和應(yīng)變。上標(biāo)R表示巖塊，上標(biāo)J表示節(jié)理，大寫表示整體坐標(biāo)系。

節(jié)理面局部坐標(biāo)系下的應(yīng)力{σj}和應(yīng)變{εj}可由整體坐標(biāo)系轉(zhuǎn)換得到，根據(jù)整體坐標(biāo)系與局部坐標(biāo)系的關(guān)系[12]，有

{σj}=[Tσ]{σJ}

(3)

{εj}=[Tε]{εJ}

(4)

式中：[Tσ]、[Tε]分別為應(yīng)力轉(zhuǎn)換矩陣和為應(yīng)變轉(zhuǎn)換矩陣，小寫的上標(biāo)j表示在局部坐標(biāo)系下的節(jié)理。

1.2 節(jié)理巖體的變形破壞模式及過程分析

1.2.1 巖塊與節(jié)理的變形破壞概化模型分析

脆性節(jié)理巖體包含巖塊和節(jié)理兩種材料，二者都可能發(fā)生破壞，其變形破壞過程對(duì)應(yīng)的應(yīng)力-應(yīng)變曲線如圖2所示。

根據(jù)圖2可將巖塊和節(jié)理的變形破壞過程概化為：材料由彈性狀態(tài)加載到峰值強(qiáng)度后，發(fā)生脆性破壞并產(chǎn)生塑性變形，應(yīng)力迅速跌落至殘余強(qiáng)度。

圖2脆性體變形破壞應(yīng)力-應(yīng)變曲線簡(jiǎn)圖

其中巖塊采用D-P屈服準(zhǔn)則進(jìn)行描述：

(5)

節(jié)理面采用去張拉強(qiáng)度的M-C屈服準(zhǔn)則進(jìn)行描述：

(6)

式中：σt為節(jié)理的法向抗拉強(qiáng)度；α、K、φ、c為材料的力學(xué)參數(shù)。

1.2.2 節(jié)理巖體的破壞模式分析

從上述巖塊和節(jié)理的變形破壞概化模型可以看出，巖塊和節(jié)理的峰值強(qiáng)度并不一致，所以在節(jié)理巖體中，兩者不一定同步破壞，由此節(jié)理巖體會(huì)呈現(xiàn)下列幾種破壞模式：

(1) 巖塊破壞而節(jié)理未破壞。由于巖塊破壞后，節(jié)理不可能仍以原來(lái)的完整狀態(tài)存在，必然受到巖塊的影響，因此本文假定巖塊破壞后節(jié)理也相應(yīng)破壞。即巖塊破壞后，節(jié)理隨著巖塊一起進(jìn)入殘余強(qiáng)度階段。

(2) 節(jié)理破壞而巖塊未破壞。節(jié)理破壞后，若巖塊仍保持完好，則巖塊的力學(xué)性能不受節(jié)理影響。即節(jié)理破壞而巖塊未破壞時(shí)，僅節(jié)理進(jìn)入殘余強(qiáng)度階段。

(3) 巖塊與節(jié)理同時(shí)破壞。同理，受巖塊破壞的影響，節(jié)理不可能以原來(lái)的完整狀態(tài)存在。巖塊與節(jié)理都進(jìn)入殘余強(qiáng)度階段。

1.2.3 節(jié)理巖體的變形破壞過程分析

由于在節(jié)理巖體中巖塊與節(jié)理面不一定同時(shí)破壞，且兩者的彈塑性狀態(tài)相互獨(dú)立，所以其中一種材料還處于彈性狀態(tài)時(shí)，另一種材料的應(yīng)力可能已達(dá)到峰值強(qiáng)度而發(fā)生脆性破壞[13]。

現(xiàn)將節(jié)理巖體的變形破壞過程劃為破壞前—應(yīng)力跌落—?dú)堄鄰?qiáng)度三個(gè)階段：

(1) 節(jié)理巖體破壞前，巖塊與節(jié)理都處于彈性狀態(tài)；

(2) 若巖塊或節(jié)理發(fā)生脆性破壞，或兩者同時(shí)破壞，稱該階段為節(jié)理巖體的應(yīng)力跌落階段；

(3) 應(yīng)力跌落階段結(jié)束后，稱之為節(jié)理巖體的殘余強(qiáng)度階段。

1.3 節(jié)理巖體變形破壞全過程的本構(gòu)關(guān)系

1.3.1 破壞前的節(jié)理巖體本構(gòu)關(guān)系

{ΔσR}和{ΔεR}分別為巖塊在整體坐標(biāo)系下的應(yīng)力、應(yīng)變?cè)隽?，{ΔσJ}和{ΔεJ}分別為節(jié)理在整體坐標(biāo)系下的應(yīng)力、應(yīng)變?cè)隽?，?/p>

(7)

(8)

節(jié)理巖體在整體坐標(biāo)系下的應(yīng)力、應(yīng)變?cè)隽糠謩e為{Δσ}和{Δε}，根據(jù)上文對(duì)節(jié)理巖體復(fù)合材料的基本假定，得到

(9)

綜上，記[De]為節(jié)理巖體的總體彈性矩陣，得到節(jié)理巖體在破壞前的本構(gòu)關(guān)系可統(tǒng)一表示為

{Δσ}=[De]{Δε}

(10)

(11)

1.3.2 應(yīng)力跌落階段的節(jié)理巖體本構(gòu)關(guān)系

應(yīng)力跌落計(jì)算引用鄭宏等[6]提出的方法：采用關(guān)聯(lián)流動(dòng)法則，巖塊或節(jié)理在跌落過程中總應(yīng)變保持不變，勢(shì)函數(shù)的非連續(xù)變化導(dǎo)致產(chǎn)生非微分量的塑性應(yīng)變?cè)隽?，其與彈性應(yīng)變?cè)隽恐蜑?。

對(duì)應(yīng)三種破壞模式，在應(yīng)力跌落階段，節(jié)理巖體的本構(gòu)關(guān)系如下：

(1) 巖塊破壞而節(jié)理未破壞。跌落時(shí)的應(yīng)力增量由巖塊產(chǎn)生，節(jié)理巖體的本構(gòu)關(guān)系為

(12)

式中：ΔλR為巖塊的塑性跌落因子，其計(jì)算方法參見鄭宏等[6]的文獻(xiàn)。

(2) 節(jié)理破壞而巖塊未破壞。跌落時(shí)的應(yīng)力增量由節(jié)理產(chǎn)生，先在局部坐標(biāo)系中計(jì)算，再轉(zhuǎn)換到整體坐標(biāo)系中，得到節(jié)理巖體的本構(gòu)關(guān)系為

(13)

若節(jié)理壓剪破壞，有

(14)

若節(jié)理張拉破壞，有

Δλj=σt/(Knd)

(15)

式中：φp和φr分別是節(jié)理的峰值內(nèi)摩擦角和殘余內(nèi)摩擦角；Ks、Kn分別是節(jié)理的切向剛度和法向剛度；σt為節(jié)理法向抗拉強(qiáng)度；d為節(jié)理間距。

(3) 巖塊與節(jié)理同時(shí)破壞。此時(shí)巖塊和節(jié)理都會(huì)應(yīng)力跌落，將跌落時(shí)的應(yīng)力增量疊加，得到節(jié)理巖體的本構(gòu)關(guān)系為

(16)

1.3.3 殘余強(qiáng)度階段的節(jié)理巖體本構(gòu)關(guān)系

巖塊或節(jié)理面破壞后，采用理想彈塑性屈服準(zhǔn)則進(jìn)行描述。若巖塊發(fā)生破壞，則巖塊改用其殘余強(qiáng)度對(duì)應(yīng)的力學(xué)參數(shù)。節(jié)理若發(fā)生壓剪破壞，則改用其殘余強(qiáng)度對(duì)應(yīng)的切向力學(xué)參數(shù)；若發(fā)生張拉破壞，則改用其殘余強(qiáng)度對(duì)應(yīng)的法向力學(xué)參數(shù)。

根據(jù)上文對(duì)節(jié)理巖體復(fù)合材料的基本假定，若巖塊或節(jié)理產(chǎn)生塑性變形，則將對(duì)應(yīng)的彈性矩陣改為彈塑性矩陣，顯然式(9)與式(10)、式(11)仍然成立。

記節(jié)理巖體的總體彈塑性矩陣為[Dep]，則在殘余強(qiáng)度階段，節(jié)理巖體的本構(gòu)關(guān)系可統(tǒng)一表示為：

{Δσ}=[Dep]{Δε}

(17)

(18)

2 節(jié)理巖體脆性破壞全過程的有限元模擬

本文針對(duì)所建立的節(jié)理巖體概化模型，利用有限元軟件ABAQUS進(jìn)行二次開發(fā)，采用FORTRAN語(yǔ)言編寫UMAT計(jì)算子程序，實(shí)現(xiàn)對(duì)節(jié)理巖體變形破壞全過程的模擬。子程序計(jì)算流程如下：

分別計(jì)算巖塊和節(jié)理的試探應(yīng)力，代入各自的屈服條件式判定彈塑性狀態(tài)。若兩種材料都處于彈性狀態(tài)，則應(yīng)力無(wú)需調(diào)整。若某種材料的應(yīng)力已達(dá)到峰值強(qiáng)度，則進(jìn)行應(yīng)力跌落計(jì)算，再結(jié)合殘余強(qiáng)度參數(shù)判斷其彈塑性狀態(tài)，進(jìn)行應(yīng)力修正并計(jì)算彈塑性矩陣，進(jìn)而求出總體彈塑性矩陣并更新應(yīng)力。最后將總體彈塑性矩陣作為雅可比矩陣，與整體應(yīng)力一起返回給主程序，子程序退出。

3 工程實(shí)例

青松引水隧洞[15]通過地區(qū)屬于涼山中高山區(qū)，隧洞全長(zhǎng)近14 km。本文選取樁號(hào)4+900的斷面為代表性斷面進(jìn)行三維有限元計(jì)算。該斷面所對(duì)應(yīng)洞段圍巖為砂巖，屬于硬巖，完整性較好；局部節(jié)理發(fā)育，一組節(jié)理與洞線大角度相交，交角取53°，傾角為41°。洞段最大埋深約1 000 m，存在較高的地應(yīng)力場(chǎng)，垂直應(yīng)力最大為26.5 MPa。

隧洞設(shè)計(jì)斷面為圓形，直徑5.08 m，計(jì)算時(shí)一次開挖。隧洞模型的上下左右各取50 m圍巖，沿隧洞線路方向取10 m長(zhǎng)度。圍巖的左、右側(cè)邊界和沿隧洞線路方向采用水平向單鏈桿約束，模型底面施加垂直向約束，頂部則視為自由邊界，但需考慮其上部巖體的壓力作用。

有限元網(wǎng)格采用8節(jié)點(diǎn)等參單元，共計(jì)18 780個(gè)單元，22 854個(gè)節(jié)點(diǎn)，具體的模型網(wǎng)格見圖3。在整體坐標(biāo)系中，Y軸正向?yàn)閺纳嫌蜗蛳掠窝厮矶捶较?，Z軸正向?yàn)榇怪毕蛏戏较颍琗軸正向以右手法則確定。令巖塊和節(jié)理的內(nèi)摩擦角在硬化過程中保持不變，節(jié)理巖體的力學(xué)參數(shù)見表1和表2。

圖3 隧洞模型有限元網(wǎng)格

表2 隧洞模型節(jié)理力學(xué)參數(shù)

開挖前先計(jì)算初始地應(yīng)力場(chǎng)，并進(jìn)行地應(yīng)力平衡。隧洞開挖后，圍巖受到地應(yīng)力釋放荷載作用，向洞內(nèi)變形。

由于硬巖巖塊的強(qiáng)度較高，開挖完成后圍巖中的完整巖塊基本沒有屈服。但由于巖體中存在大量的薄弱節(jié)理面，因而在隧洞頂部和底部仍然有部分巖體發(fā)生屈服破壞，其破壞模式以沿節(jié)理傾向的滑移為主。圍巖的屈服區(qū)分布見圖4，由于節(jié)理走向與洞線斜交，屈服區(qū)左右分布不對(duì)稱。

圖4圍巖屈服區(qū)分布

圖5給出隧洞開挖后圍巖的合位移增量等值線圖和位移矢量圖。從圖5中可看到，隧洞開挖后圍巖向洞內(nèi)變形，洞頂及洞底的位移增量最大，最大值為2.95 cm。同樣由于節(jié)理走向與洞線斜交，位移增量的左右分布不對(duì)稱。由于隧洞開挖后圍巖中巖塊仍處于彈性狀態(tài)，巖體變形不大，能保持穩(wěn)定。

圖5圍巖的位移分布

計(jì)算中以拉應(yīng)力為正，壓應(yīng)力為負(fù)，隧洞開挖后圍巖中的應(yīng)力以壓力為主，圖6給出圍巖的最小主應(yīng)力等值線分布圖。從圖6中可看出，受薄弱節(jié)理面的影響，隧洞開挖后，洞頂及洞底周圍部分巖體發(fā)生屈服破壞，導(dǎo)致這兩處圍巖的主應(yīng)力值較小。圍巖的最大壓應(yīng)力數(shù)值約63.93 MPa，位于隧洞洞腰處中。同樣由于節(jié)理走向與洞線斜交，應(yīng)力的左右分布不對(duì)稱。

圖6圍巖的最小主應(yīng)力分布

綜合上述計(jì)算結(jié)果，表明對(duì)于強(qiáng)度很高的硬巖，若其中有薄弱節(jié)理面存在，則隧洞開挖時(shí)圍巖的變形破壞主要由節(jié)理所控制。圍巖容易沿節(jié)理發(fā)生剪切錯(cuò)動(dòng)，導(dǎo)致圍巖變形不斷增大，最終發(fā)生屈服破壞[16]。實(shí)際工程中，考慮脆性破壞的節(jié)理巖體能很好反映自身的變形破壞情況。

4 結(jié) 論

(1) 將節(jié)理巖體概化為連續(xù)復(fù)合材料，分析脆性節(jié)理巖體的三種破壞模式。根據(jù)各破壞模式分別建立對(duì)應(yīng)的節(jié)理巖體變形破壞本構(gòu)關(guān)系?；趶椝苄岳碚摰玫焦?jié)理巖體的應(yīng)力修正方法，提高了計(jì)算效率?；贏BAQUS編寫UMAT計(jì)算子程序，實(shí)現(xiàn)對(duì)節(jié)理巖體變形破壞全過程的模擬。

(2) 對(duì)青松引水隧洞的開挖過程進(jìn)行數(shù)值模擬，證明了本文的理論模型可運(yùn)用于實(shí)際工程，用該方法能較好地模擬脆性節(jié)理巖體的變形與破壞，正確反映隧洞開挖后圍巖的屈服破壞情況。

(3) 通過計(jì)算結(jié)果得出，在強(qiáng)度較高的硬巖中，圍巖的破壞主要是由節(jié)理破壞所導(dǎo)致，同時(shí)表明薄弱節(jié)理面對(duì)節(jié)理巖體的位移、應(yīng)力，和屈服破壞都有重大影響。