危巖體斷裂特性的網(wǎng)格重構(gòu)數(shù)值仿真研究

程 靜

(浙江同濟(jì)科技職業(yè)學(xué)院水利工程學(xué)院，浙江 杭州 311231)

危巖體廣泛存在于我國(guó)西部地區(qū)，其主控裂隙往往是影響其穩(wěn)定性的重要因素。危巖致災(zāi)因素較多，如降雨、地震和自重作用等。危巖失穩(wěn)災(zāi)害屢見(jiàn)不鮮，如2008年我國(guó)汶川地震，發(fā)生了鄰近高速公路的危巖失穩(wěn)下墜，導(dǎo)致多輛汽車被砸毀，造成了人員重大傷亡事故；2009年7月我國(guó)四川省汶川縣國(guó)道213 線徹底關(guān)大橋發(fā)生了危巖失穩(wěn)墜落，導(dǎo)致大橋斷落，造成6人死亡、交通癱瘓，經(jīng)濟(jì)損失過(guò)億元。因此，對(duì)危巖體在復(fù)雜工況下的穩(wěn)定性進(jìn)行分析是防治危巖失穩(wěn)的關(guān)鍵所在。

針對(duì)危巖體的斷裂力學(xué)模型分析方面，國(guó)內(nèi)外學(xué)者進(jìn)行了大量有益的研究工作。如王林峰等構(gòu)建了危巖體的震動(dòng)理論模型，對(duì)地震作用下危巖體的穩(wěn)定性進(jìn)行了分析；李部等對(duì)滑塌式危巖體的主控結(jié)構(gòu)面斷裂擴(kuò)展方向和穩(wěn)定性進(jìn)行了理論分析；王景宏對(duì)鏈子崖危巖體的地質(zhì)特性及其破壞機(jī)制進(jìn)行了探討；李克森等基于斷裂力學(xué)與水力學(xué)理論對(duì)危巖體主控面斷裂的滲流特性進(jìn)行了分析研究。同時(shí)，相關(guān)規(guī)范也運(yùn)用赤平極射投影及模糊數(shù)學(xué)建立了危巖體穩(wěn)定性的評(píng)價(jià)標(biāo)準(zhǔn)。

但是以上研究不能實(shí)現(xiàn)危巖體裂紋擴(kuò)展全過(guò)程的評(píng)價(jià)。隨后人們利用各種方法對(duì)危巖體裂紋擴(kuò)展過(guò)程進(jìn)行了數(shù)值仿真研究，如陳洪凱等利用斷裂力學(xué)理論對(duì)危巖體裂紋擴(kuò)展過(guò)程進(jìn)行了數(shù)值仿真研究，得到了危巖體裂紋擴(kuò)展及應(yīng)力強(qiáng)度因子的變化規(guī)律；鄭安興等利用擴(kuò)展有限元對(duì)降雨及自重作用下危巖體斷裂擴(kuò)展過(guò)程進(jìn)行了數(shù)值分析；陳鵬宇等利用離散元軟件PFC對(duì)危巖體破壞過(guò)程進(jìn)行了數(shù)值仿真研究。另外，還有諸如近場(chǎng)動(dòng)力學(xué)法(PD)、數(shù)值流形法等，以上方法解決了部分問(wèn)題，但是有些方法根本不是基于斷裂力學(xué)理論，還有些方法使用起來(lái)具有諸多限制，如需要預(yù)設(shè)裂紋擴(kuò)展路徑，且編程復(fù)雜，難于推廣使用。

鑒于以上研究的不足之處，本文基于Abaqus軟件平臺(tái)，利用Python語(yǔ)言開(kāi)發(fā)了基于斷裂力學(xué)理論的危巖體任意裂紋自由擴(kuò)展的網(wǎng)絡(luò)重構(gòu)數(shù)值模擬方法，該方法利用圍道積分計(jì)算每一步裂紋前緣的應(yīng)力強(qiáng)度因子(SIF)，根據(jù)最大主應(yīng)力法則判斷危巖體裂紋的擴(kuò)展方向，并利用網(wǎng)格重剖分技術(shù)實(shí)現(xiàn)了危巖體裂紋的擴(kuò)展?；谝陨暇W(wǎng)絡(luò)重剖分?jǐn)?shù)值模擬新方法對(duì)降雨條件下危巖體的斷裂過(guò)程及應(yīng)力變形規(guī)律進(jìn)行評(píng)價(jià)，并對(duì)危巖體在不同裂紋長(zhǎng)度及裂紋角度的斷裂過(guò)程和應(yīng)力變形規(guī)律進(jìn)行敏感性分析，研究結(jié)果可為降雨條件下危巖體穩(wěn)定性分析提供參考，同時(shí)可為降雨條件下危巖體任意裂紋擴(kuò)展的數(shù)值仿真技術(shù)研究提供新的思路。

1 計(jì)算理論與數(shù)值手段

1. 1 圍道積分計(jì)算危巖體裂紋尖端的應(yīng)力強(qiáng)度因子

在準(zhǔn)靜態(tài)分析的背景下，

J

積分在二維中定義為

(1)

式中：

Γ

為從底部裂紋表面開(kāi)始到頂部表面結(jié)束的輪廓，如圖1(a)所示，

Γ

→0表明輪廓收縮到裂紋尖端;

q

為在虛擬裂紋擴(kuò)展方向上是否存在單位向量；

n

為

Γ

的法線方向；

H

可以表達(dá)為下式：

(2)

其中：對(duì)于彈性材料，

W

定義為彈性應(yīng)變能，對(duì)于彈塑性或彈黏塑性材料，

W

定義為彈性應(yīng)變能密度加上塑性耗散，從而表示“等效彈性材料”中的應(yīng)變能；

I

為剛度矩陣；

σ

為裂紋尖端的應(yīng)力場(chǎng)；

u

為裂紋尖端的位移場(chǎng)；

x

為裂紋尖端的位置坐標(biāo)。

J

積分可以表達(dá)為

(3)

裂紋尖端應(yīng)力強(qiáng)度因子在線性彈性斷裂力學(xué)中起著重要的作用，它們描述了荷載或變形對(duì)裂紋尖端應(yīng)力和應(yīng)變場(chǎng)大小的影響，并測(cè)量了裂紋擴(kuò)展傾向或裂紋驅(qū)動(dòng)力。此外，裂紋尖端應(yīng)力強(qiáng)度因子與能量釋放速率(

J

積分)有關(guān)：

(4)

式中：

K

為裂紋尖端應(yīng)力強(qiáng)度因子，

K

=[

K

,

K

,

K

]；

B

為前對(duì)數(shù)能量因子矩陣，對(duì)于均質(zhì)、各向同性材料，

B

是對(duì)角陣。

因此公式(4)可以簡(jiǎn)化為

(5)

因此，公式(5)建立了

J

積分與裂紋尖端應(yīng)力強(qiáng)度因子之間的關(guān)系，可以根據(jù)該式計(jì)算裂紋尖端的應(yīng)力強(qiáng)度因子。

圖1 裂紋尖端圍線示意圖Fig.1 Schematic diagram of the crack tip girth

1.2 基于網(wǎng)格重剖分的危巖體裂紋擴(kuò)展數(shù)值方法

危巖體的裂紋擴(kuò)展數(shù)值模擬方法在各種商業(yè)軟件中均有內(nèi)置，但是對(duì)裂紋擴(kuò)展具有諸多限制，如需要預(yù)設(shè)裂紋擴(kuò)展路徑，無(wú)法真正意義上實(shí)現(xiàn)危巖體裂紋的任意擴(kuò)展。同時(shí)，以往數(shù)值模擬方法多采用巖體損傷理論或者彈塑性理論，無(wú)法體現(xiàn)危巖體裂紋尖端的斷裂力學(xué)特征。因此，本文基于圍線積分及Abaqus平臺(tái)實(shí)現(xiàn)了基于斷裂力學(xué)的危巖體任意裂紋擴(kuò)展，并利用Python語(yǔ)言實(shí)現(xiàn)了危巖體裂紋尖端應(yīng)力強(qiáng)度因子提取、裂紋擴(kuò)展準(zhǔn)則判斷、裂紋重構(gòu)及網(wǎng)格重剖分等操作，可達(dá)到危巖體任意裂紋自由擴(kuò)展的目的，同時(shí)可以真實(shí)地反映危巖體裂紋尖端的斷裂力學(xué)特性。基于網(wǎng)絡(luò)重剖分的危巖體裂紋擴(kuò)展數(shù)值模擬方法流程圖如圖2，具體步驟如下：

圖2 基于網(wǎng)格重剖分的裂紋擴(kuò)展數(shù)值方法流程圖Fig.2 Flow chart of numerical method of crack propagation of dangerous rock mass based on grid reconstruction

(1) 建立初始裂紋模型，分配材料屬性，建立Crack的歷程輸出，建立邊界條件及荷載，劃分網(wǎng)格。

(2) 對(duì)初始裂紋模型進(jìn)行靜力有限元計(jì)算，判斷裂紋尖端的應(yīng)力強(qiáng)度因子是否達(dá)到材料的斷裂韌度，若達(dá)到材料的斷裂韌度，則裂紋發(fā)生擴(kuò)展，刪除原始模型網(wǎng)格，裂紋往前擴(kuò)展一個(gè)增量，重新劃分網(wǎng)格。

(3) 利用Abaqus中的*map solution功能將前一步的模型映射到新一步模型的網(wǎng)格上，施加剩余荷載，進(jìn)行下一步裂紋擴(kuò)展。

1. 3 最大周向應(yīng)力準(zhǔn)則

危巖體的裂紋擴(kuò)展需要解決裂紋擴(kuò)展方向的判別問(wèn)題。均質(zhì)、各向同性線彈性材料的近裂紋尖端應(yīng)力場(chǎng)可以表達(dá)為

(6)

(7)

式中：

σ

為環(huán)向應(yīng)力；

τ

為切向應(yīng)力；

r

和

θ

為與裂紋前沿正交的平面上以裂紋尖端為中心的極坐標(biāo)。

利用上述這兩種條件都可以得到裂紋的擴(kuò)展方向，即：

(8)

2 重慶市萬(wàn)州太白危巖體穩(wěn)定性分析

2. 1 工程概況

本文以重慶市萬(wàn)州太白編號(hào)為W15的危巖體(見(jiàn)圖3)為例，利用上述建立的危巖體任意裂紋擴(kuò)展的網(wǎng)絡(luò)重構(gòu)數(shù)值模擬方法對(duì)該危巖體的穩(wěn)定性進(jìn)行數(shù)值仿真分析，以驗(yàn)證該方法的有效性。該地區(qū)多年平均降雨量為1 000 mm，對(duì)危巖體治理過(guò)程中發(fā)現(xiàn)，在2004年5月份由于暴雨導(dǎo)致了危巖體多處地方崩落，因此研究暴雨工況與不同裂隙賦存情況下危巖體的斷裂特性對(duì)于防治危巖失穩(wěn)災(zāi)害十分關(guān)鍵。

圖3 危巖體示意圖Fig.3 Schematic diagram of dangerous rock mass

2. 2 計(jì)算模型及模型概況

該危巖體的高度為20.5 m，長(zhǎng)度為9.5 m，厚度為4 m，其中裂隙位于危巖體中部，裂隙長(zhǎng)度為

L

(m)，裂紋方向與豎直方向的夾角為

α

(°)。危巖體由長(zhǎng)石石英砂組成，其重度為24.8 kN/m，巖體的抗拉強(qiáng)度為480 kPa，巖體的彈性模量為8 000 kPa，其泊松比為0.25，巖體的Ⅰ型斷裂韌度為1.796×10N/m。為了研究在不同裂隙傾角

α

、不同裂隙長(zhǎng)度

L

以及不同降雨強(qiáng)度下危巖體的斷裂力學(xué)特性，將裂隙傾角

α

分別取為0°、30°、60°，裂隙長(zhǎng)度

L

分別取為2 m和4 m，降雨分為天然工況與暴雨工況，為了簡(jiǎn)化計(jì)算，根據(jù)文獻(xiàn)[16]的研究結(jié)果，天然工況下裂隙內(nèi)部水壓力取為裂隙深度的1/3，而暴雨工況下裂隙內(nèi)部水壓力則取為裂隙深度的2/3。

2. 3 計(jì)算結(jié)果分析

2.3.1 不同工況下危巖體的應(yīng)力及位移分析

對(duì)初始條件下的不同裂隙傾角、裂紋長(zhǎng)度以及降雨工況下危巖體裂紋尖端的應(yīng)力云圖進(jìn)行分析，如圖4至圖6所示，為了簡(jiǎn)化篇幅及便于閱讀，圖4給出了裂隙傾角α分別為0°、30°、60°，裂隙長(zhǎng)度

L

為2 m，以及暴雨工況下危巖體的最大主應(yīng)力、剪應(yīng)力和水平位移云圖；圖5給出了裂隙傾角

α

為0°,裂隙長(zhǎng)度

L

分別為2 m和4 m以及暴雨工況下危巖體的最大主應(yīng)力、剪應(yīng)力和水平位移云圖；而圖6則給出了裂隙傾角

α

為0°,裂隙長(zhǎng)度

L

為2 m，以及天然工況與暴雨工況下危巖體的最大主應(yīng)力、剪應(yīng)力和豎向位移云圖。由圖4至圖6可見(jiàn)，無(wú)論何種工況，危巖體的應(yīng)力在裂紋尖端形成應(yīng)力集中，同時(shí)其水平位移在裂紋處呈現(xiàn)不連續(xù)特征，且危巖體主控面的左側(cè)水平位移較小(幾乎為0)，而其右側(cè)水平位移較大，表明在自重及裂隙內(nèi)水壓的作用下，危巖體主控面發(fā)生了張開(kāi)效應(yīng)；對(duì)于不同裂隙角度

α

而言，隨著裂隙角度的增大，裂隙對(duì)危巖體變形的影響逐漸減小，由危巖體水平位移云圖可見(jiàn)其水平位移主要集中于危巖體的突出部位；對(duì)于不同裂隙長(zhǎng)度

L

而言，隨著裂隙長(zhǎng)度的增加，危巖體的應(yīng)力集中程度越大，同時(shí)其最大水平位移主要發(fā)生在裂隙部位；而對(duì)于天然工況與暴雨工況而言，暴雨工況下應(yīng)力集中程度更大，變形更劇烈。

對(duì)圖4至圖6中不同工況下危巖體的最大主應(yīng)力、剪應(yīng)力及水平位移進(jìn)行統(tǒng)計(jì)，其結(jié)果列于表1。

表1 不同工況下危巖體的應(yīng)力及變形統(tǒng)計(jì)

由表1可知：對(duì)于不同裂隙角度而言，危巖體的最大主應(yīng)力值和剪應(yīng)力值隨著裂隙角度的增大呈現(xiàn)先增后減的趨勢(shì)，表明相同條件下在裂隙角度接近30°左右時(shí)危巖體的應(yīng)力集中程度最大，而當(dāng)裂隙角度進(jìn)一步增大后危巖體的應(yīng)力集中程度明顯降低，但是其水平位移值隨著裂隙傾角的增大呈現(xiàn)持續(xù)減小的趨勢(shì)，表明裂隙傾角越小對(duì)危巖體整體穩(wěn)定性的影響越大；對(duì)于不同裂隙長(zhǎng)度而言，裂隙長(zhǎng)度的增大對(duì)危巖體應(yīng)力及變形的影響最大，本文

L

=4 m工況較之

L

=2 m工況危巖體的最大主應(yīng)力值和剪應(yīng)力值分別增長(zhǎng)了88.49%和102.91%，而其水平位移值則增加了104.65%，可見(jiàn)裂隙長(zhǎng)度的增大對(duì)危巖體穩(wěn)定性的影響最大；對(duì)于不同荷載工況而言，暴雨工況增大了危巖體的應(yīng)力集中及變形程度，較之天然工況，暴雨工況下危巖體的最大主應(yīng)力和剪應(yīng)力值增長(zhǎng)了16.25%和30.89%，其水平位移值增加了2.38%，可見(jiàn)暴雨一定程度上增加了危巖體的失穩(wěn)程度，但是其影響要小于危巖體的裂隙長(zhǎng)度。

圖4 不同裂隙傾角下危巖體的應(yīng)力云圖(單位：kPa)Fig.4 Stress cloud map of dangerous rock under different fracture dip angles(unit:kPa)

2.3.2 危巖體裂紋擴(kuò)展過(guò)程分析

限于篇幅，本文即給出危巖體較為危險(xiǎn)的工況：即不同裂隙長(zhǎng)度(

L

=2 m和4 m)工況下危巖體裂隙擴(kuò)展的全過(guò)程，見(jiàn)圖7。

圖5 不同裂隙長(zhǎng)度下危巖體的應(yīng)力云圖(單位：kPa)Fig.5 Stress cloud map of dangerous rock under different fracture lengths(unit:kPa)

圖6 天然工況與暴雨工況下危巖體的應(yīng)力云圖(單位：kPa)Fig.6 Stress cloud map of dangerous rock under natural and rainstorm conditions(unit:kPa)

由圖7可見(jiàn)，危巖體的主控面(即圖7中危巖體的裂隙)在自重及暴雨的作用下，逐漸向右下方擴(kuò)展，且靠近危巖體的尖點(diǎn)處，裂隙長(zhǎng)度越大，危巖體掉落塊體的體積也越大。

圖7 危巖體裂紋擴(kuò)展的全過(guò)程Fig.7 Whole process of crack propagation in dangerous rock mass

為了探究危巖體裂紋擴(kuò)展過(guò)程中裂紋尖端應(yīng)力強(qiáng)度因子的變化規(guī)律，對(duì)圖7中兩種工況下危巖體裂紋尖端的應(yīng)力強(qiáng)度因子隨裂紋擴(kuò)展長(zhǎng)度的變化進(jìn)行了追蹤，其結(jié)果見(jiàn)圖8。

圖8 兩種工況下危巖體裂紋尖端的應(yīng)力強(qiáng)度因子 隨裂紋擴(kuò)展長(zhǎng)度的變化曲線Fig.8 Variation curves of stress intensity factor at crack tip of dangerous rock mass with crack propagation length under two working conditions

由圖8可見(jiàn)，隨著裂紋長(zhǎng)度的不斷擴(kuò)展，危巖體裂紋尖端的Ⅰ型應(yīng)力強(qiáng)度因子和Ⅱ型應(yīng)力強(qiáng)度因子的絕對(duì)值也在不斷地增大，同時(shí)危巖體裂紋尖端的Ⅰ型應(yīng)力強(qiáng)度因子與Ⅱ型應(yīng)力強(qiáng)度因子在危巖體裂紋擴(kuò)展初期的量級(jí)一致，但是隨著危巖體裂紋長(zhǎng)度不斷擴(kuò)展后，危巖體裂紋尖端的Ⅰ型應(yīng)力強(qiáng)度因子的絕對(duì)值也不斷增大，而其Ⅱ型應(yīng)力強(qiáng)度因子的絕對(duì)值則變化不大，表明危巖體裂紋擴(kuò)展過(guò)程中Ⅰ型擴(kuò)展占主導(dǎo)地位，而值得注意的是，初始裂紋長(zhǎng)度越長(zhǎng)，危巖體裂紋尖端的Ⅰ型應(yīng)力強(qiáng)度因子和Ⅱ型應(yīng)力強(qiáng)度因子的絕對(duì)值越大，表明此時(shí)危巖體越容易失穩(wěn)。

3 結(jié) 論

本文基于Abaqus軟件平臺(tái)，利用Python語(yǔ)言開(kāi)發(fā)編制了基于斷裂力學(xué)理論的危巖體裂紋任意擴(kuò)展的網(wǎng)絡(luò)重構(gòu)數(shù)值模擬方法，并利用該方法對(duì)重慶市萬(wàn)州太白編號(hào)為W15的危巖體穩(wěn)定性進(jìn)行了數(shù)值仿真分析，得到了以下結(jié)論：

(1) 基于圍道積分及斷裂力學(xué)理論，對(duì)Abaqus軟件進(jìn)行了二次開(kāi)發(fā)，利用Python語(yǔ)言實(shí)現(xiàn)了危巖體裂紋擴(kuò)展、方向判斷及應(yīng)力強(qiáng)度因子計(jì)算的全過(guò)程，為危巖體斷裂變形特性模擬提供了一種新的數(shù)值模擬方法。

(2) 危巖體的最大主應(yīng)力值和剪應(yīng)力值隨著裂隙角度的增大呈現(xiàn)先增后減的趨勢(shì)，危巖體的水平位移值隨著裂隙傾角的增大持續(xù)減小；隨著裂隙角度的增大，裂隙對(duì)危巖體變形的影響逐漸減小，其水平位移主要集中于危巖體的突出部位。

(3) 危巖體裂隙長(zhǎng)度的增大對(duì)危巖體應(yīng)力及變形的影響最大，使得變形由危巖體的突出部位向裂隙處集中；暴雨工況下危巖體的變形及位移較之天然工況要大，是影響危巖體穩(wěn)定的重要因素之一。

(4) 危巖體裂紋朝著危巖體突出部位進(jìn)行擴(kuò)展，主要受到Ⅰ型擴(kuò)展主導(dǎo)，隨著裂紋的不斷擴(kuò)展，危巖體裂紋尖端的Ⅰ型應(yīng)力強(qiáng)度因子呈現(xiàn)指數(shù)型上升，Ⅱ型應(yīng)力強(qiáng)度因子變化不大，初始裂紋長(zhǎng)度越長(zhǎng)，危巖體裂紋尖端的Ⅰ型應(yīng)力強(qiáng)度因子和Ⅱ型應(yīng)力強(qiáng)度因子的絕對(duì)值越大。