地震作用下隧道洞口段巖質邊坡損傷演化分析

趙杰　徐劍　王桂萱　尹訓強

摘要：從巖石材料細觀非均勻性特點出發，運用連續介質損傷力學理論，推導了巖石材料細觀損傷本構關系式。以ANSYS二次開發平臺為主體，基于細觀損傷力學，提出基于宏觀均質假定的等效損傷單元來模擬隧道洞口段巖質邊坡的損傷開裂過程。計算結果表明：在細觀方面邊坡損傷破壞是個不斷損傷的過程，在地震荷載作用下，邊坡結構首先會出現微裂紋，隨著地震動時間增加，已出現的微裂紋逐漸擴展形成可見的裂縫，最后形成連續的損傷區域；尤其在隧道拱頂上方形成連續貫通的損傷區域，可看出拱頂已完全破壞，給隧道安全運行帶來不利影響，因此在設計施工時要對拱頂上方進行加固處理。該計算方法和分析成果可用于同類復雜結構的損傷破壞分析。 關鍵詞：巖質邊坡；非均質性；損傷力學；等效損傷單元；損傷破壞；ANSYS

中圖分類號：U45 文獻標識碼：A 文章編號：1000-0666（2016）01-0053-07

0 引言

近年來發生的重大巖石動力災害事件表明：礦山巷道、地下洞室、隧道工程、邊坡工程垮塌和破壞，多數是由于巖體內的節理、裂隙及其擴展造成（羅國煌等，2004；李明田等，2003）。巖體內的節理、裂隙就是其損傷的實質性表現。所以，從細觀損傷理論出發，研究圍巖體細觀損傷擴展機理和演化特征對巖體工程穩定性的評價與支護設計是非常有意義的。

Kachanov在1958年提出損傷力學思想，至今已經過50多年的發展，一些理論已日益成熟并在各學科和工程技術領域得以應用。周維垣等（1998）應用連續介質損傷力學理論，從巖體內部微裂紋產生和擴展的損傷機理出發，建立了堅硬巖體的彈性一損傷耦合的各向異性彈脆性損傷本構模型。朱萬成等（2003）從巖石的細觀非均勻性特點出發，應用彈性損傷力學對巖石材料細觀的損傷演化過程進行了描述。趙寶友等（2003）基于損傷力學原理，通過引入損傷變量的方法詳細推導了ABAQUS中的損傷塑性本構模型屈服階段及其卸載后的損傷機制和力學行為，并以某水電站地下廠房洞室結構為例，進行二維動力時程非線性分析，得出強震下采用此損傷塑性模型是合理和必要的。隋斌等（2008）利用FLAC3D軟件模擬了地震作用下某地下洞室群圍巖的動態響應，采用新的劈裂破壞判據對震后可能出現的劈裂損傷區范圍進行了預測。左雙英和肖明（2009）基于深埋大型地下洞室群圍巖損傷機理，對映秀灣水電站洞室群在三向地震荷載耦合作用下的動力響應進行數值模擬，得出地震作用下圍巖損傷分布的特征。張志國等（2010）采用動力時程分析法研究地下洞室地震響應特性，結合三維彈塑性損傷有限元理論，編寫了動力時程計算程序，分析了漁子溪水電站地下廠房洞室群的圍巖穩定。

筆者在上述研究的基礎上，基于巖石材料的宏觀均質假定，運用連續介質彈性損傷力學模型來描述單元的材料特性，采用最大拉應力與摩爾一庫倫破壞準則作為損傷閾值，以通用有限元軟件ANSYS開發平臺為主體，運用UPFs的二次開發特點，基于用戶自定義單元的功能將損傷單元法引入到ANSYS中，用來分析地震作用下隧道洞口段巖質邊坡從微裂紋產生、擴展到形成可見裂縫，直至最終的損傷過程。

1 損傷演化本構模型

1.1 巖石材料非均質性特征

巖石作為一種自然地質體，由于各種因素使巖石內部存在大量不同的微缺陷，所以巖石介質的構成非常復雜，通常對其進行數學描述也是非常困難的。唐春安（1993）認為可以利用統計的方法對離散后的巖石介質進行近似描述。建立材料細觀非均質特性隨機模型時，所謂非均質性是指材料的力學參數，如強度、彈性模量、泊松比等。在給定的空間上對模型進行網格劃分，為了描述材料參數的這種空間變化，需要解決兩個問題：第一，如何使材料參數在數值上滿足非均勻分布規律；第二，如何將這些不同材料參數隨機地分布到劃分的單元中去。本文將巖質邊坡看成一個大的樣本空間，采用等效損傷單元對巖石材料進行離散，使得離散后的單元形式具有類似于普通有限單元的特性，并且認為其本身為均勻、連續的介質，然后在樣本空間內引入概率統計的方法，假定各單元的材料力學性質存在差異且服從Weibull概率統計分布，從而可表現出一定的非均勻型和離散性。楊天鴻等（2004）認為采用Weibull隨機分布可以較好地模擬材料性質的非均勻性，該分布可以用分布密度函數來定義：式中，μ為巖石材料隨機分布力學性質參數（強度、彈性模量、泊松比等）；μ⁰為巖石材料基元體力學性質平均值；m為巖石力學性質分布函數的密度形狀函數，其物理意義反應了巖石材料的均質性；f（μ）是巖石力學性質μ的統計分布函數。

式（1）很好地反映了巖石材料的細觀力學性質非均勻性分布情況。隨著均勻性系數m的增加，基元體的力學性質將集中于一個狹窄的范圍之內，表明單元強度趨于均勻；當m值減小時，基元體的力學性質分布范圍變寬，表明巖石介質的性質趨于均勻?？梢妋值反映了統計模型中材料結構的離散程度，能夠較為真實地反映巖石材料的非均質性。

1.2 巖體本構關系與屈服準則

余大慶（1993）認為巖石應力-應變曲線是由于其受力后的不斷損傷引起微裂紋萌生和擴展而造成的，而不是由于其塑性變形。尤其是在拉伸應力作用下，其脆性更加明顯。因此，用彈性損傷力學本構模型來描述巖石材料細觀單元力學性質是合適的。在本文的研究方法中，等效單元應力一應變本構關系亦采用了彈性損傷模型。在施加地震荷載前，用彈性模型和泊松比來定義每一個等效損傷單元，并且將其假定為線彈性、均值、無損傷的。當單元進入損傷階段后，依據Le-maitre應變等價原理，認為應力作用在受損材料上引起的應變與有效應力作用在無損材料引起的應變等價。據這一原理，受損材料單元的本構關系可通過無損材料中的名義應力得到：式中，E⁰為材料初始彈性模量；E為材料損傷后的彈性模量，D為損傷變量。D=0表示材料處于無損傷狀態，D=1表示處于完全損傷（斷裂或者破壞）狀態，00

在巖石材料損傷演化中，隨著等效損傷單元應力會不斷增加，當單元的應力或應變狀態滿足某個給定的損傷閾值（準則）時，單元開始損傷。這里選擇兩個損傷閾值準則，第一個采用最大拉應變準則，當細觀單元的最大拉伸主應變達到其給定的應變閥值時，該單元開始發生拉伸損傷；第二個采用摩爾一庫侖準則，當細觀單元的應力狀態滿足摩爾一庫侖準則時，該單元發生剪切損傷。需要說明的是，巖石材料在外部荷載作用下，其抗拉強度遠低于抗壓強度的特點，且單元的拉伸損傷和剪切損傷不會同時發生于同一荷載步下。換句話說拉伸準則具有優先權，若細觀單元滿足拉伸準則，則不需再判斷其是否滿足摩爾一庫侖準則。只有在損傷單元未滿足拉伸準則之后才判定其是否滿足摩爾一庫侖則。圖1、2分別給出了單元材料拉伸本構關系曲線與壓縮本構關系曲線圖。單元滿足拉損傷模型以及摩爾一庫侖準則對應的損傷變量表達為式中，f_tr為單元參與強度；ε_t0為拉損傷應變閾值；ε_tu為單元材料極限拉伸應變；ε_c0為單元應力達到單元材料單軸抗壓強度時的主要應變；ε_cu為單元材料的極限應變；ε_c為單元的單軸壓縮應變；N為與下降段形狀有關的參數。

2 等效損傷單元在ANSYS中的實現

本文等效損傷單元采用的本構關系與屈服準則是其與普通有限單元最大的區別之處?；诘刃卧厥庑?，需要以ANSYS為開發平臺，利用UPFs為開發工具進行二次開發。圖3給出了基于ANSYS為二次開發平臺的等效損傷單元建立流程，其中虛線方框部分由UPFs的功能實現。詳細的程序流程如下：

首先，建立隧道洞口段巖質邊坡有限元模型，對用戶進行單元編譯和鏈接，施加位移邊界條件，控制求解流程以及結構的輸出（此部分利用GUI和APDL完成）；然后，進入地震動力分析求解階段，在每一荷載步施加地震荷載向量；最后，在ANSYS數據庫中將建立單元的所有信息輸入到接口子程序UserElem.f中，并完成以下步驟（此部分利用子程序接口UserElem.f完成）：

①檢驗是否第一次進入等效損傷單元子程序，若“是”，則利用蒙特卡羅方法生成服從Weibull分布的材料參數隨機數，然后生成最終的單元材料力學參數，將此數據存儲在數組中并使用FOR-TRAN語言中的COMMON將其定義為全局變量，以便以后的平衡迭代步和荷載步調用，進入②；若“否”，則將直接獲得上一迭代步的材料力學參數，進入②。

②通過節點真實位移{u}_t+{△u}^it計算單元的應變s；由損傷變量Dⁱ_l（Dⁱ_lt表示拉損傷；Dⁱ_lc表示壓損傷，其中i表示當前迭代步，l表示當前荷載步）判斷該單元是否發生過損傷，即Dⁱ_l是否大于D⁰_l，其中，0是上一荷載步的最后迭代步，若“否”，則說明該單元此前未發生損傷，進入③；若“是”，則說明該單元此前發生過損傷，進入④。

③由②所求應變s檢驗是否發生拉損傷，如發生了則進入⑤，若未發生拉損傷，則檢驗是否為壓損傷，如發生了則進入⑥，未發生ANSYS數據庫中，計算整體結構響應。

④判斷發生損傷的類型：若為拉損傷，計算壓損傷變量Dⁱ_lc。

⑤計算相應的拉損傷變量Dⁱ_lt以及弱化后彈性模量E。

⑥計算相應的壓損傷變量Dⁱ_lc以及弱化后彈性模量E。

⑦計算單元矩陣，包括單元剛度陣、阻尼陣、質量陣以及單元等效荷載向量，然后將上面得到數據返回到ANSYS數據庫中，計算整體結構響應。如果計算結果不滿足收斂條件，則進入下一平衡迭代步，若滿足則進入下一荷載步，直至地震分析結束。

3 巖質邊坡損傷演化模擬

3.1 邊坡計算參數及有限元模型

計算參數選?。簭椥阅Ａ縀=3000GPa，泊松比v=0.23，質量密度ρ=2300kg/m³，動抗拉強度f_t=1.5×10⁵。選取材料的隨機力學參數為：彈性模量和抗拉強度的均質度m=2，泊松比的均質度m=100。選取瑞利阻尼進行計算，計算時取第1階與第10階的固有頻率，陣型阻尼比為0.05。本文所建立的邊坡模型坡度為45°，坡高為20m，坡前計算長度取20m，坡后計算長度取44m，坡腳以下取20m。為研究隧道洞口段巖質邊坡在地震作用下的損傷破壞過程，在平面應變狀態下，采用等效損傷單元（幾何形狀如四節點平面應變單元）對邊坡模型進行有限元網格離散，離散后的單元尺寸約為0.4m×0.4m，模型共分18589個單元、18941個節點。模型左右兩邊施加水平約束，底部施加水平和豎向約束，邊坡有限元模型見圖4。分水平與豎直2個方向輸入地震波，其中水平峰值加速度為0.474g，豎向加速度峰值為0.312g，地震動總的持續時間為10.56s，加速度時程曲線見圖5。

3.2 邊坡地震損傷演化分析

本文從地震荷載作用整個過程中，選取6個典型的時間點（2.64s、2.92s、5.78s、6.96s、8.98s、9.42s）來分析隧道洞口段巖質邊坡從開始出現微裂紋到形成清晰可見的裂縫，裂縫不斷擴展至形成最終損傷區域總過程，其損傷演化過程及其主拉應力的分布見圖6。隧道洞口段巖質邊坡損傷演化分析過程分析中：在靜力荷載（本文僅施加了重力荷載）作用的基礎上施加地震荷載。地震荷載激勵的初期階段，邊坡絕大部分區域處于受壓狀態下，雖然結構沒有出現嚴重的裂縫，但由于單元材料隨機分布的影響，少數材料強度較弱的單元依舊出現了不同程度的損傷，坡面開始出現微裂紋，但邊坡結構整體上處于完好的狀態。隨著地震激勵時間的增加，已出現的微裂紋不斷擴展，新的微裂紋也不斷出現，這些微裂紋逐漸擴展、延伸，最后在坡后形成數條清晰可見的裂縫，如圖6a、b所示；隨著地震動的繼續震蕩，坡后裂縫范圍繼續擴展、延伸，坡面發生破壞區域逐漸擴大；洞口拱頂上方和坡面中間開始出現裂紋并擴展成較大裂縫，形成多處貫通的損傷區域，同時造成了主拉應力的重新分布，如圖6c、d所示；最后，隨著地震動的繼續震蕩，拱頂上方裂縫繼續擴展，最終形成了貫穿整個拱頂上方的損傷區域，此時可以認為拱頂上方已完全破壞，會給隧道的安全運行帶來影響；坡面中間形成的自上而下裂縫會繼續擴展、延伸，最終形成了一條清晰可見的裂縫，可以近似將該裂縫看成邊坡在強震作用下產生的滑坡趨勢，如圖6e、f所示。

4 結論

本文通過引入服從Weibull函數的細觀損傷微元體，運用連續損傷介質理論，理論上推導出了可反映巖石非均質性和宏觀破裂與失穩過程的巖石細觀損傷本構關系式。

地震作用下巖質邊坡的破壞是材料細觀損傷不斷產生、累積的結果。在地震荷載激勵的初期階段，邊坡絕大部分區域處于受壓狀態下，雖然結構沒有出現嚴重的裂縫，但由于單元材料隨機分布的影響，少數材料強度較弱的單元依舊出現了不同程度的損傷，坡后開始出現微裂紋，這些微裂紋逐漸貫通并形成數條清晰可見的裂縫。隨著地震動的繼續震蕩，在洞口拱頂上方和坡面中間開始出現新的裂縫。拱頂上方裂縫最終形成一個連續貫通的損傷區域，此時可以認為拱頂上方已完全破壞，這會給隧道的安全運行帶來影響，因此在洞口設計、施工時要對拱頂上方進行必要的加固處理；坡面中間形成的自上而下的裂縫，可以看成是邊坡結構在強震作用下產生的滑坡趨勢，設計、施工也要給與重視。