危巖體斷裂特性的網格重構數值仿真研究

程 靜

(浙江同濟科技職業學院水利工程學院，浙江 杭州 311231)

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

針對危巖體的斷裂力學模型分析方面，國內外學者進行了大量有益的研究工作。如王林峰等構建了危巖體的震動理論模型，對地震作用下危巖體的穩定性進行了分析；李部等對滑塌式危巖體的主控結構面斷裂擴展方向和穩定性進行了理論分析；王景宏對鏈子崖危巖體的地質特性及其破壞機制進行了探討；李克森等基于斷裂力學與水力學理論對危巖體主控面斷裂的滲流特性進行了分析研究。同時，相關規范也運用赤平極射投影及模糊數學建立了危巖體穩定性的評價標準。

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

鑒于以上研究的不足之處，本文基于Abaqus軟件平臺，利用Python語言開發了基于斷裂力學理論的危巖體任意裂紋自由擴展的網絡重構數值模擬方法，該方法利用圍道積分計算每一步裂紋前緣的應力強度因子(SIF)，根據最大主應力法則判斷危巖體裂紋的擴展方向，并利用網格重剖分技術實現了危巖體裂紋的擴展?；谝陨暇W絡重剖分數值模擬新方法對降雨條件下危巖體的斷裂過程及應力變形規律進行評價，并對危巖體在不同裂紋長度及裂紋角度的斷裂過程和應力變形規律進行敏感性分析，研究結果可為降雨條件下危巖體穩定性分析提供參考，同時可為降雨條件下危巖體任意裂紋擴展的數值仿真技術研究提供新的思路。

1 計算理論與數值手段

1. 1 圍道積分計算危巖體裂紋尖端的應力強度因子

在準靜態分析的背景下，

J

積分在二維中定義為

(1)

式中：

Γ

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

Γ

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

q

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

n

為

Γ

的法線方向；

H

可以表達為下式：

(2)

其中：對于彈性材料，

W

定義為彈性應變能，對于彈塑性或彈黏塑性材料，

W

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

I

為剛度矩陣；

σ

為裂紋尖端的應力場；

u

為裂紋尖端的位移場；

x

為裂紋尖端的位置坐標。

J

積分可以表達為

(3)

裂紋尖端應力強度因子在線性彈性斷裂力學中起著重要的作用，它們描述了荷載或變形對裂紋尖端應力和應變場大小的影響，并測量了裂紋擴展傾向或裂紋驅動力。此外，裂紋尖端應力強度因子與能量釋放速率(

J

積分)有關：

(4)

式中：

K

為裂紋尖端應力強度因子，

K

=[

K

,

K

,

K

]；

B

為前對數能量因子矩陣，對于均質、各向同性材料，

B

是對角陣。

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

(5)

因此，公式(5)建立了

J

積分與裂紋尖端應力強度因子之間的關系，可以根據該式計算裂紋尖端的應力強度因子。

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

1.2 基于網格重剖分的危巖體裂紋擴展數值方法

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

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

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

(2) 對初始裂紋模型進行靜力有限元計算，判斷裂紋尖端的應力強度因子是否達到材料的斷裂韌度，若達到材料的斷裂韌度，則裂紋發生擴展，刪除原始模型網格，裂紋往前擴展一個增量，重新劃分網格。

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

1. 3 最大周向應力準則

危巖體的裂紋擴展需要解決裂紋擴展方向的判別問題。均質、各向同性線彈性材料的近裂紋尖端應力場可以表達為

(6)

(7)

式中：

σ

為環向應力；

τ

為切向應力；

r

和

θ

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

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

(8)

2 重慶市萬州太白危巖體穩定性分析

2. 1 工程概況

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

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

2. 2 計算模型及模型概況

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

L

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

α

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

α

、不同裂隙長度

L

以及不同降雨強度下危巖體的斷裂力學特性，將裂隙傾角

α

分別取為0°、30°、60°，裂隙長度

L

分別取為2 m和4 m，降雨分為天然工況與暴雨工況，為了簡化計算，根據文獻[16]的研究結果，天然工況下裂隙內部水壓力取為裂隙深度的1/3，而暴雨工況下裂隙內部水壓力則取為裂隙深度的2/3。

2. 3 計算結果分析

2.3.1 不同工況下危巖體的應力及位移分析

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

L

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

α

為0°,裂隙長度

L

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

α

為0°,裂隙長度

L

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

α

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

L

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

對圖4至圖6中不同工況下危巖體的最大主應力、剪應力及水平位移進行統計，其結果列于表1。

表1 不同工況下危巖體的應力及變形統計

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

L

=4 m工況較之

L

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

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

2.3.2 危巖體裂紋擴展過程分析

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

L

=2 m和4 m)工況下危巖體裂隙擴展的全過程，見圖7。

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

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

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

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

為了探究危巖體裂紋擴展過程中裂紋尖端應力強度因子的變化規律，對圖7中兩種工況下危巖體裂紋尖端的應力強度因子隨裂紋擴展長度的變化進行了追蹤，其結果見圖8。

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

由圖8可見，隨著裂紋長度的不斷擴展，危巖體裂紋尖端的Ⅰ型應力強度因子和Ⅱ型應力強度因子的絕對值也在不斷地增大，同時危巖體裂紋尖端的Ⅰ型應力強度因子與Ⅱ型應力強度因子在危巖體裂紋擴展初期的量級一致，但是隨著危巖體裂紋長度不斷擴展后，危巖體裂紋尖端的Ⅰ型應力強度因子的絕對值也不斷增大，而其Ⅱ型應力強度因子的絕對值則變化不大，表明危巖體裂紋擴展過程中Ⅰ型擴展占主導地位，而值得注意的是，初始裂紋長度越長，危巖體裂紋尖端的Ⅰ型應力強度因子和Ⅱ型應力強度因子的絕對值越大，表明此時危巖體越容易失穩。

3 結 論

本文基于Abaqus軟件平臺，利用Python語言開發編制了基于斷裂力學理論的危巖體裂紋任意擴展的網絡重構數值模擬方法，并利用該方法對重慶市萬州太白編號為W15的危巖體穩定性進行了數值仿真分析，得到了以下結論：

(1) 基于圍道積分及斷裂力學理論，對Abaqus軟件進行了二次開發，利用Python語言實現了危巖體裂紋擴展、方向判斷及應力強度因子計算的全過程，為危巖體斷裂變形特性模擬提供了一種新的數值模擬方法。

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

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

(4) 危巖體裂紋朝著危巖體突出部位進行擴展，主要受到Ⅰ型擴展主導，隨著裂紋的不斷擴展，危巖體裂紋尖端的Ⅰ型應力強度因子呈現指數型上升，Ⅱ型應力強度因子變化不大，初始裂紋長度越長，危巖體裂紋尖端的Ⅰ型應力強度因子和Ⅱ型應力強度因子的絕對值越大。