FLAC3D中錨桿剪切斷裂的數值實現與應用

張志強，劉 銀，趙梓彤，李化云，2

(1.西南交通大學 交通隧道工程教育部重點實驗室, 四川 成都 610031；2.西華大學 土木建筑與環境學院，四川 成都 610039)

作為新奧法隧道支護體系中唯一能深入圍巖內部的支護結構，錨桿在地下工程圍巖穩定性控制方面有著廣泛的應用[1-4]。對處于大埋深、高地應力等復雜地質條件下的地下工程(交通隧道、煤礦巷道、水工隧道等)而言，當圍巖中存在軟弱結構面時，結構面兩側巖體的相互錯動會使穿過結構面的錨桿承受較大的剪切作用，此時，錨桿容易產生剪切斷裂而失去對圍巖形變的控制能力，甚至導致整個支護體系失效，引發圍巖大變形災害[5-7]。

國內外學者針對錨固體中錨桿的抗剪機制以及地下工程中錨桿的剪斷機理進行了許多研究。張偉等[8]基于錨桿的切向和軸向形變特征，建立了節理巖體錨固體中錨桿的剪切變形分析模型；Li等[9]通過室內試驗，對比分析了玻璃纖維錨桿與螺紋鋼錨桿錨固體在剪切荷載作用下的力學行為；陳文強等[10]建立了加錨結構面抗剪強度計算模型，對剪切過程中錨桿的軸向和橫向作用進行了分析；Forbes等[11]通過新的光纖測試技術，對錨固體剪切過程中錨桿的力學行為進行了分析；肖同強等[6]通過建立肩角錨桿力學模型，揭示了深部構造應力區煤巷肩角錨桿剪斷機制；繆協興等[12]分析了綜放沿空巷道頂部全長黏結式錨桿承受剪切變形的機理；劉金海等[13]研究發現，深埋巷道肩角、肩窩處錨桿最容易在剪應力作用下產生斷裂。同時，數值計算也常用于對錨桿的支護效果進行分析。但目前的數值計算中[14-16]，均未考慮錨桿在剪切荷載作用下的斷裂效應。因此，導致計算結果與實際支護效果不符，夸大錨桿的作用。

FLAC3D軟件在巖土工程大變形計算方面有著廣泛的應用。根據特定研究目的，研究者采用FISH語言對軟件自帶的結構單元進行了針對性的二次開發。如李為騰等[17]將拉斷準則引入FLAC3D中的CABLE單元，實現了錨桿的局部拉伸斷裂；劉少虹[18]提出了爆破震動下錨桿單元破壞損傷判據，建立了基于FLAC3D平臺的錨桿軸力動力損失修正模型；言志信[19]等將界面破壞時砂漿的應力軟化效應引入錨桿結構單元，分析了地震作用下巖質邊坡錨固界面的剪切失效過程。但這些二次開發模型均未考慮錨桿的剪斷效應。

因此，本文基于FLAC3D軟件平臺，對PILE結構單元進行修正，使其能準確地反映錨固體中錨桿剪切形變的全過程；然后將修正PILE單元用于地下工程圍巖穩定性分析中，對錨桿的剪斷與拉斷失效機制進行研究。

1 FLAC3D中PILE單元力學特性

1.1 PILE單元基本特性

FLAC3D中，單個PILE單元由2個節點以及節點間的結構單元構成，結構單元具有相同的截面形狀及物理力學特性。當PILE單元被錨固劑包裹時，其與錨固劑之間通過法向與切向彈簧—滑塊體系來實現形變的耦合與內力的傳遞，如圖1所示。

圖1 錨桿單元結構組成

1.2 PILE單元剪力學特性

剪切荷載作用下，PILE單元承受的剪力隨著剪切位移的增加而線性增長。剪力與剪切應變之間的關系為

Fs=γGA

(1)

式中：Fs為錨桿桿體所承受的剪力；γ為剪切應變；G為桿體的剪切剛度；A為桿體的截面積。

1.3 PILE單元拉伸力學特性

PILE單元中，引入錨桿的拉斷準則，即

εt≥εtmax

(2)

式中：εt為桿體的應變；εtmax為桿體的極限應變，其取值與材料的形變能力有關。

桿體的軸向本構關系如圖2所示，其拉伸表達式為

(3)

式中：Ft為桿體的拉力；Et為桿體的彈性模量；Ftmax為桿體的極限抗拉承載力。

可見，錨桿桿體存在以下3種受拉力學狀態。

狀態1：桿體所受到的外力小于其極限抗拉能力，且桿體的應變小于其極限應變。此時桿體處于彈性階段，桿體軸力隨著形變值的增加而線性增長。

狀態2：桿體承受的外力達到其極限承載能力，而桿體的應變小于其極限應變。此時桿體軸力保持不變，應變繼續增加。

狀態3：桿體的應變達到極限應變。此時桿體產生斷裂，失去承載能力。

通過上述對PILE單元的分析可以發現，PILE單元可以較好地模擬錨桿的拉伸力學效應，但不能真實地反映錨桿的剪力學效應。

圖2 拉伸力學特性

2 PILE單元剪斷力學模型的實現

2.1 錨桿剪力學模型

Li等[7]通過對錨桿進行剪切試驗，發現錨桿的剪力與剪切位移之間存在著圖3(a)所示的關系。可以看出，在桿體剪力達到極限承載能力之前，剪力隨著剪切位移的增加而近似線性增大(OA段)；當桿體剪力達到剪切極限力時，桿體的剪力會突降至其最大殘余承載能力(AB段)；在BC段，隨著剪切位移的增加，桿體承擔的剪力出現波動，但整體上呈現出減小的趨勢；當桿體的剪切位移超過其極限剪切位移時，錨桿產生斷裂，失去承載能力。

根據文獻[6]的試驗結果，可得到簡化的錨桿剪力學模型如圖3(b)所示，其表達式為

(4)

式中：Fs,max為桿體所能承擔的極限剪切荷載；γmax為桿體極限剪切應變；Fs1和Fs2為峰后特征點B和C的剪力；γs1為剪力峰值點對應的剪切應變。

圖3 錨桿剪力學模型

2.2 PILE單元剪斷力學模型的數值實現

為了使PILE單元能夠真實地反映錨桿的剪力學效應，采用FISH語言，以式(4)為基礎對PILE單元的剪力學行為進行修正，修正流程如圖4所示。具體步驟如下。

(1)采用主程序計算得到每一個桿體單元的剪力與剪切應變。

(2)判斷桿體單元受到的剪力是否超過其極限承載能力，若超過，則根據桿體單元的剪切應變通過插值計算得到此時桿體所能承擔的最大剪力；同時，判斷桿體的剪切應變是否超過其極限剪切應變，若超過，則將桿體的剪切極限力設為0；錨桿單元的軸力按式(3)進行計算。

(3)判斷計算是否收斂，若收斂，則計算結束；若不收斂，則重復步驟(1)和步驟(2)。

圖4 PILE單元剪斷開發流程

3 錨桿剪斷力學模型的驗證

以文獻[7]的錨桿剪切試驗結果為基礎，對本文錨桿剪斷力學模型的有效性進行驗證。建立的數值模型如圖5所示，混凝土試樣采用實體單元進行模擬，尺寸為0.4 m×0.4 m×0.4 m。混凝土單元的下端固定，上端為自由面，其余側面施加相應方向的法向約束。錨桿直徑為28 mm，所能承受的最大剪切荷載Fsmax=279 kN，峰后剪力特征值Fs1=150 kN，Fs2=80 kN，極限抗拉承載力Ftmax=558 kN，極限剪切應變γmax=83 mm。錨桿長度L=0.55 m，劃分成15個PILE桿體單元，桿體的CID編號為1～15。桿體與錨固劑參數見表1和表2。

表1 錨桿桿體力學參數

表2 錨固劑力學參數

圖5 剪切數值模型

計算過程中，在A點施加沿x方向的恒定速度，大小為1×10-6m·計算步-1，對編號為11的PILE單元，記錄其不同剪切位移的剪力計算值，結果如圖6所示；同時，將試驗值也繪于圖6中。

圖6 剪力與剪切位移曲線室內試驗與數值計算結果對比

從圖6可以看出：修正模型可以較好地反映出錨桿剪切過程中的彈性階段、塑性變形階段以及破壞階段；彈性階段內，當桿體單元產生同樣的剪切位移時，數值計算得到的剪力小于室內試驗結果，且兩者之間的差值隨著剪切位移的增加而減小；塑性變形階段，室內試驗的波動幅度較數值計算大，剪力整體上均呈現出減小的特征；破壞階段，錨桿的剪力變為0。可以看出，雖然本文的修正模型認為錨桿的殘余抗剪強度與剪切位移之間為線性關系(圖3(b)中BC段)，但由于錨桿在塑性剪切變形階段內形變值會持續波動，導致數值計算結果得到的剪力在該階段內也持續波動，但整體上體現出線性減小的趨勢。

原始PILE單元、修正PILE單元最終形變與室內試驗實物的結果對比如圖7所示。可以看出：原始PILE單元不能反映錨桿在剪力作用下的非連續大變形特征，而修正PILE單元可以較好地反映這一形變特征；修正PILE單元得到的計算結果與室內試驗結果吻合較好，可以較好地反映錨桿剪切形變全過程的力學特征。

圖7 最終形變室內試驗與數值計算結果對比

4 錨桿剪斷力學模型的應用

單弱面層狀巖體是隧道工程建設中常見的一類巖體，采用本文建立的修正PILE單元，分析該隧道中錨桿剪斷及拉斷效應對圍巖穩定性的影響。

4.1 數值模型

建立的數值模型如圖8所示。模型左右寬度均為40 m(約為4倍開挖跨度)，仰拱以下高度為40m，模型厚度為2 m。前、后、左、右邊界以及上邊界施加相應方向的法向應力邊界，下邊界約束豎向位移，地應力取值為σx=σy=σz=10 MPa。層理弱面與水平方向的夾角θ=0°，45°，90°，層理面的強度與剛度取值為圍巖基質體的1/2。圍巖基質體參數見表3。

圖8 層狀巖體隧道計算模型

表3 復合體及圍巖基質體的參數

噴射混凝土采用實體單元模擬，鋼拱架的作用按其彈性模量折算給噴射混凝土[20]，即

(5)

式中：E為噴射混凝土—鋼拱架復合體彈性模量；E0為噴射混凝土彈性模量；Ag為鋼拱架截面積；Eg為鋼材彈性模量；Sc為噴射混凝土截面積。

等效后的復合體力學參數見表3-復合體初期支護的厚度為20 cm。錨桿長度為3.7 m，其布置方式如圖8(b)所示，錨桿的環向間距為1 m，縱向間距為1.2 m。為了后續描述方便，對錨桿從左下角沿逆時針方向依次進行編號，編號為C1—C24。錨桿為全長黏結式錨桿，桿體和錨固劑的力學參數分別見表1和表2。采用PILE 結構單元模擬錨桿，為了提高模擬精度，每根錨桿劃分為37 段，每段長0.1 m，桿體所能承擔的極限剪切荷載Fsmax=140 kN，抗拉載荷為255 kN，拉伸延伸率為16%。

4.2 計算結果分析

4.2.1 錨桿斷裂力學特性

以θ=45°為例，通過對錨桿剪斷單元(錨桿編號為C3，CID編號為77)與錨桿拉斷單元(錨桿編號為C2，CID編號為43)的拉力與剪力隨計算步的變化情況進行記錄，分析錨桿的斷裂力學特性，結果如圖9所示，可以看出，剪斷錨桿單元與拉斷錨桿單元的力學行為有著較大的差異。對于拉伸斷裂單元，其計算過程中的拉力最大值為其抗拉能力極限值；對于剪切斷裂單元，計算過程中的拉力最大值小于抗拉極限值。

1)剪斷錨桿單元

單元的剪力變化過程(圖9(a))為：初始變形階段，剪力隨著計算步的增加而近似線性增加；當計算至1 120步時，剪力為140 kN，達到剪切破壞臨界值；此后，單元在下一計算步內斷裂，剪力驟降至殘余抗剪能力最大值(60 kN)；當計算至1 310步時，剪力降至40 kN，此時單元的應變達到其剪切極限應變；隨著計算步的繼續增加，錨桿出現剪切斷裂，剪力從40 kN降至0 kN，之后隨著計算步的增加保持不變。

單元的拉力變化過程(圖9(b))為：初始形變階段(計算步小于400步)，拉力近似線性增加；當計算步大于400步后，拉力出現劇烈的波動上升；當錨桿剪切斷裂時，單元的拉力從190 kN降至0 kN，之后隨著計算步的增加保持不變。

可見，當錨桿單元剪斷時，單元的拉力并未達到抗拉能力極限值(225 kN)。因此，在數值計算中，對于受剪切荷載作用較大的錨桿而言，若不考慮錨桿剪斷效應，則會同時夸大錨桿的抗拉與抗剪支護能力。

2)錨桿拉斷單元

單元的拉力變化過程(圖9(c))為：拉力先近似線性增加(計算步小于400步)，后波動增加(400步到890步)；當計算至890步時，單元所受到的拉力達到了抗拉能力最大值，即255 kN，此時，由于單元的拉伸應變未達到其極限拉伸應變，因此隨著計算步的進一步增加，單元的拉力保持在255 kN 不變；當計算至1 390步時，單元的應變達到極限應變，此時錨桿拉伸斷裂，拉力降至0 kN，后隨著計算步的增加拉力保持不變。

單元的剪力變化過程(圖9(d))為：在錨桿拉伸斷裂之前，單元所受到的剪力均較小，最大值僅為4.8 kN，且隨著計算步的增加波動變化。當錨桿拉伸斷裂時，單元的剪力也降為0 kN，后隨著計算步的增加剪力保持不變。可見，對于以抗拉作用為主的錨桿單元而言，當剪切荷載小于極限抗剪能力時，不考慮錨桿的剪斷效應不會對錨桿的支護效果產生影響。

3)內力沿桿體的分布規律

剪切斷裂錨桿(C3)與拉伸斷裂錨桿(C2)各單元內力沿桿體的分布規律如圖10所示。從圖中可以看出：當錨桿發生破斷后，錨桿的剪力沿長度方向基本全變為0 kN，說明此時錨桿在全長范圍內失去了抗剪切的能力；而錨桿的拉力沿長度方向的分布規律為中間大兩端小，說明此時錨桿還能對圍巖的局部拉伸變形起到一定的抑制作用。

圖9 錨桿單元剪斷和拉斷力學特性

圖10 錨桿內力分布

4)錨桿支護體系斷裂過程

計算至平衡狀態時，所有產生斷裂的錨桿其斷裂位置處單元的拉力隨時間變化曲線如圖11所示。從圖可以看出：錨桿的斷裂順序與斷裂模式為C3→C22→C2→C4→C20→C5→C21→C23→C9→C16→C18→C8(粗體為拉伸斷裂，下同)，即首先是左右邊墻處的錨桿被剪斷，之后剪斷錨桿相鄰位置處的錨桿被拉斷，然后左右拱腰處錨桿被拉斷，最后左右拱肩處錨桿被剪斷。可見，該工況下，邊墻和拱腰附近錨桿容易產生拉伸或剪切破壞，拱肩附近錨桿容易產生剪切破壞。

4.2.2 不同傾角條件下錨桿的斷裂過程

θ=0°，90°時，錨桿的斷裂過程如圖12所示。從圖12(a)可以看出：當θ=0°時，圍巖未產生大的相對剪切滑動，錨桿只產生了拉伸斷裂，斷裂順序為C3(C21)→C5(C18)→C8(C16)→C9(C15，C10，C14)，即邊墻附近錨桿首先發生斷裂，然后，錨桿的斷裂部位集中分布在拱腰與拱頂之間。當θ=90°時，錨桿的斷裂順序為C10→C15→C11→C9→C16→C19(C18)→C6(C7)→C14，即拱肩與拱頂之間的錨桿首先產生剪切斷裂，然后拱腰與拱肩之間的錨桿產生拉伸斷裂。

4.2.3 與原始PILE單元結果的比較

1)錨桿受力

以θ=45°為例，采用原始PILE單元的錨桿受力如圖13所示。由圖13可知：由于未考慮錨桿剪斷，C3錨桿編號為77的桿體單元所承受的剪力隨著計算步的增加持續增加，承受的拉力先持續增加至抗拉能力極限值，之后保持拉力極限值不變；從錨桿斷裂過程來看，錨桿的拉伸斷裂順序為C2→C4→C6(C22,C21)→C20→C17，集中在邊墻與拱腰附近。

圖11 錨桿破斷過程

圖12 不同傾角下錨桿破斷過程

圖13 原始模型錨桿受力特性(θ=45°)

將圖13與采用改進PILE單元的圖9和圖10對比可知：錨桿斷裂數量少于改進PILE單元計算結果，且斷裂開始的時間要滯后于改進PILE單元的計算結果。

2)圍巖特征位移

各種工況下圍巖的特征位移見表4，可以看出，對于不同層理面傾角，修正PILE單元與原始PILE單元計算結果之間的差異具有不同的規律。當θ=0°時，由于錨桿未產生剪切斷裂，修正模型與原始模型的計算結果相等。當θ=45°時，修正模型的拱頂沉降與洞周收斂比原始模型分別小17和44 mm，洞周收斂之間的差值更大；這是由于錨桿的剪切斷裂主要發生在左右邊墻附近，錨桿剪斷會更多地減弱支護系統對圍巖水平形變的約束能力。當θ=90°，修正模型的拱頂沉降與洞周收斂比原始模型分別小27和17 mm，拱頂沉降之間的差值更大；這是由于錨桿的剪切斷裂主要發生在拱肩與拱頂附近，錨桿剪斷會更多地減弱支護系統對圍巖拱頂豎向形變的約束能力。

表4 圍巖特征位移值

5 結 論

本文基于FLAC3D平臺，將錨桿剪切斷裂模型引入軟件自帶的PILE結構單元，對PILE單元的剪力學行為進行修正。采用修正的PILE單元對單弱面層狀巖體中錨桿的剪斷與拉斷失效機制進行了研究，得到以下結論。

(1)提出了PILE單元在剪切荷載作用下的全過程修正力學模型。該模型增加了錨桿剪切變形的塑性階段與破裂階段，即當剪力達到抗剪能力極限值時，錨桿進入塑性殘余變形階段；當剪切應變達到極限剪切應變時，錨桿剪切破斷，桿體的軸力和剪力突降至零并保持恒定。

(2)錨桿剪斷力學模型的驗證結果表明，修正PILE單元得到的計算結果與室內試驗結果吻合較好，可以較好地反映錨桿剪切形變全過程的力學特征以及剪切荷載作用處錨桿的非連續大變形。

(3)單弱面層狀巖體隧道錨桿支護應用的對比分析結果表明，對于受剪切荷載作用較大的錨桿而言，不考慮錨桿的剪斷效應會同時夸大錨桿的抗拉與抗剪支護能力。與不考慮錨桿剪斷相比，錨桿剪斷后，支護系統對圍巖形變約束能力會有所減弱。但由于錨桿剪斷位置的差異，不同層理面傾角條件下，支護系統約束能力減弱的程度與層理面的傾角有關。