應用非連續模型研究斷層對地殼應力的影響

(中南大學 地球科學與信息物理學院 計算地學研究中心，湖南 長沙，410083)

地殼應力是地震預測和地震危險性評價的一個重要依據，也是地質環境、地質工程、巖土工程設計和施工的重要基礎參數[1?2]。地質構造運動會引起地殼形變和應力，地應力局部集中會造成地震活動和巖體破壞。若巖體內存在斷裂等不連續面, 則其附近的應力狀態會發生復雜的變化[3?5]。地應力是一個地區的區域穩定性和巖體穩定性的重要因素，它對建筑物、構筑物和地下工程的安全有直接的影響。圍巖發生失穩破壞大多沿巖體不連續面即節理、斷層或者附近進行。在巖體高應力區，邊坡開挖、隧洞穿鑿、巷道采礦、地下空間工程等巖體開挖會引起一系列與應力釋放相關的巖體變形和破壞現象，惡化巖體的工程地質條件，直接威脅到施工工程的安全和穩定[6?8]。因此，對地應力的研究受到工程地質界的廣泛重視。在工程力學模擬計算中，往往采用連續變形算法，這對斷裂構造的模擬有很大的不足，在理論計算和工程分析中都難免出現相當大的誤差甚至錯誤[9]，特別是當存在一些大斷裂和斷裂兩側地殼變形躍變劇烈時，仍然采用連續變形算法是不夠的。計算機技術和計算方法的發展，不連續體概念和不連續變形理論逐漸成熟，為工程力學數值模擬和計算地球科學提供了新的方法[10?13]。在此，本文作者采用非連續接觸單元模擬斷層，根據庫侖摩擦準則判斷剪切滑動和計算剪切應力，通過有限元數值模擬方法來探討斷層各要素對地殼應力的影響。

1 理論和方程

地殼巖石往往被斷層分割成眾多離散的塊體，塊體內部的變形可以采用固體力學方程來描述，而不同塊體彼此之間則表現為斷層接觸關系。模型遵循的物理力學方程如下[14]。

(1) 平衡方程。物體平衡的力學條件為：

(2) 幾何方程。物體在應力作用下產生形變，應變和位移滿足如下幾何關系：

(3) 本構方程。對于各向同性線彈性材料，應力與應變對載荷的響應符合胡克定律：

其中：σij為應力；εij為應變；fi為體力；xi為坐標；ui為位移；θ為體應變；λ和μ為拉梅系數，λ=Eν/[(1+ν)(1?2ν)]，μ=E/[2(1+ν)]；E和ν分別為彈性模量和泊松比；δij為 Kronecker變量，若i=j，則δij=1；若i≠j，則δij=0。

(4) 接觸應力方程。斷層往往不只是一個簡單的錯動面，而是相對圍巖稍微軟弱的斷層帶。斷層帶和圍巖自身的變形是連續的，但二者之間是不連續的接觸關系。由于斷層帶和圍巖的接觸區域通常并不完全確定，所以，本文模型采用點/面接觸方式來處理二者的接觸關系。將斷層帶與圍巖的接觸邊界相互視為對方的接觸面和目標面，從而組成接觸對，然后，通過接觸力學方法分析斷層的受力和運動狀態[15]。接觸分析需要計算垂直于目標面的法向接觸應力和平行于目標面的切向接觸應力，采用Pinball算法進行計算。

接觸面和目標面之間的間隙(穿透量)用g來表示，當接觸面穿過目標面時，就發生接觸穿透，規定此時g為負值。采用罰函數法計算法向接觸應力：

其中：fn為法向接觸應力；Kn為法向接觸剛度。

切向接觸應力是由接觸面在目標面上移動所產生的摩擦力引起的，若接觸面沿目標面的切向位移的彈性分量為，則切向接觸應力為：

其中：ft為切向接觸應力；Kt為切向接觸剛度；F為靜態/動態摩擦因子；τ為庫侖滑動摩擦力：τ=κfn；κ為摩擦因數。

2 模型和邊界條件

2.1 基本假設

對于發育斷層的二維平面地殼模型，進行以下假設：

(1) 斷層帶與圍巖都為各向同性的線彈性材料；

(2) 斷層帶相對圍巖較軟弱；

(3) 斷層帶與圍巖為非連續的接觸關系；

(4) 斷層帶與圍巖之間的相對滑動符合庫侖摩擦定律；

(5) 頂板和底板的垂向剛度極大而水平剪切剛度極小。

2.2 模型的建立

采用二維平面模型，考慮水平面上的正方形巖體和其中的不同幾何形態的非貫通型斷層。巖體長×寬為1 km×1 km，斷層帶寬度為10 m。對于具有一定寬度的斷層帶，非連續模型將其外邊界處理為與巖體緊挨的接觸面。

2.3 模型力學參數

真實材料參數本應從實驗中獲取，但不同時代、不同地區、不同深度、不同性質巖石的參數都不盡相同。本文進行一般性討論，巖體和斷層帶的密度、彈性模量、泊松比和摩擦因數在合理范圍內取值[16]。對于接觸力學參數，如法向接觸剛度因子、切向接觸剛度因子和靜態/動態摩擦因子，本文均取 1.0，這適合于ANSYS求解大多數接觸問題[15]。巖體和斷層接觸面的摩擦因數為0.4。除了在分析斷層帶力學參數對地殼應力的影響時需要改變斷層帶的彈性模量和泊松比外，所有模型均采用表1所示的參數。

表1 有限元模型力學參數Table 1 Mechanical parameters of finite element model

2.4 邊界條件

模型受到區域最大主應力σ1和最小主應力σ3的擠壓作用，σ1=10 MPa，σ3=5 MPa。模型左邊界和前部邊界在法向上被約束，右邊界施加最大主應力σ1，后部邊界施加最小主應力σ3。所有邊界剪切自由。

巖體中發育的1條平直斷層模型的有限元分網及邊界條件如圖1所示。其他模型有與此相同的邊界條件，只是斷層幾何形態有所不同。

圖1 有限元模型及邊界條件Fig.1 Finite element model with boundary conditions

3 計算結果及分析

3.1 連續模型和非連續模型地殼應力狀態的比較

連續模型與非連續模型的區別在于后者考慮了斷層帶與圍巖之間不連續運動的接觸關系，在結合處通過接觸單元把二者聯系起來，二者在接觸位置可以發生錯動；而對于連續模型，斷層帶與圍巖在結合處具有公共節點，因此，二者之間不會發生錯動(圖2)。

圖2 連續模型與非連續模型示意圖Fig.2 Schematic diagrams of continuous and discontinuous model

采用相同的模型參數和邊界條件分別計算連續模型和非連續模型的地殼應力，結果表明無論是連續模型還是非連續模型，斷層及圍巖的應力狀態表現出 3個相同的特點：(1) 弱化了斷層帶兩端局部較小區域的最大主應力；(2) 斷層帶具有較高的剪切應力；(3) 最大主應力和剪切應力的高值都位于斷層端部。但是，非連續模型在更大程度上明顯地影響了斷層端部的最大主應力和剪切應力分布狀態。沿著斷層的最大主應力和剪切應力進行對比，發現非連續模型比連續模型更加突出了斷層對其附近地殼應力的影響作用(圖3)。上述特征表明：非連續模型能夠更好地模擬斷層，有效克服了單純使用軟弱帶模擬斷層而不考慮其不連續運動所導致的應力削弱的問題。

3.2 連續模型和非連續模型斷層走滑變形特征的比較

圖3 連續模型和非連續模型地殼應力沿斷層的變化Fig.3 Changes of crustal stress along fault of continuous and discontinuous model

斷層走滑量可以通過計算斷層兩側與斷層帶接觸的巖體單元節點在斷層走滑方向上的位移之差得到，而且可以根據位移關系推測斷層的走滑性質。圖4所示為連續模型和非連續模型斷層走滑方向上的位移沿 1條垂直于斷層的跨斷層測線的變化(測線見圖1，方向為左上—右下)。從圖4可以看出：在給定的區域構造應力作用下，連續模型和非連續模型都反映出斷層發生了左旋走滑作用，但二者給出的走滑量相差很大：非連續模型斷層走滑量為1.73 cm，而連續模型的走滑量僅為0.32 cm；非連續模型斷層走滑變形由3段構成，其中第1和第3段為斷層帶與圍巖之間的不連續錯動，第2段為斷層帶內部的連續變形，且走滑量主要為不連續錯動所貢獻；而連續模型斷層僅有斷層帶的連續變形這一部分，由于斷層帶與圍巖在結合處具有公共節點，故二者并沒有發生錯動。野外地質觀察到的大量擦痕和階步是斷層不連續錯動的具體體現，而非連續模型可以通過接觸單元描述斷層帶與圍巖之間的這種不連續變形，進一步表明非連續模型比連續模型能夠更好地模擬斷層。

圖4 連續模型和非連續模型斷層走滑方向上的位移沿跨斷層測線的變化(測線見圖1)Fig.4 Changes of displacement in fault strike-slipping direction along test line A-B striding over fault of continuous and discontinuous kinematic model(test line shown in Fig.1)

由于斷層端部是地殼應力的高值分布區即應力集中區，該處應力狀態能夠較好地反映斷層對地殼應力的影響。

3.3 斷層力學參數對地殼應力的影響

3.3.1 彈性模量的影響

采用相同的非連續模型計算只改變斷層彈性模量時地殼應力的變化，結果如圖5所示。從圖5可以看出：隨著彈性模量的增大，斷層端部的最大主應力和剪切應力明顯減小，盡管最大主應力并沒有像剪切應力那樣一直遞減，但數值的變化是顯而易見的；斷層中部的最大主應力變化很大，但剪切應力的變化較小；斷層遠處的最大主應力和剪切應力的變化都很小，表明斷層端部地殼應力對斷層的彈性模量非常敏感。確切地說，單獨改變斷層彈性模量實際上是改變了圍巖與斷層彈性模量的差別。當斷層彈性模量取較小值時，斷層與圍巖的彈性模量的差別較大，斷層端部具有較大的最大主應力和剪切應力。因此推斷：在相同的區域主應力作用下，圍巖與斷層的彈性模量相差越大，斷層端部的地殼應力就越大。

3.3.2 泊松比的影響

采用相同的非連續模型計算只改變斷層泊松比時地殼應力的變化，如圖6所示。從圖6可以看出：斷層中部和遠處的最大主應力和剪切應力幾乎都不隨泊松比的改變而變化。盡管斷層端部應力有隨泊松比增大而減小的微弱趨勢，但變化量非常小，表明斷層泊松比對地殼應力的影響是微弱的。

圖5 非連續模型斷層附近地殼應力與斷層彈性模量的關系Fig.5 Relationship between crustal stress near faults and elastic modulus of faults for discontinuous model

圖6 非連續模型斷層附近地殼應力與斷層泊松比的關系Fig.6 Relationship between crustal stress near faults and poisson ratio of faults for discontinuous model

3.4 斷層與區域主應力具有不同夾角時斷層端部地殼應力的變化

采用相同的力學參數計算改變斷層與區域最大主應力的夾角(θ)時斷層端部地殼應力的變化，如圖7所示。從圖7可以看出：斷層端部最大主應力和剪切應力都隨θ的改變而變化。最大主應力在θ約為45°時具有極大值，剪切應力在θ約為25°和70°時具有極大值。盡管應力極大值對應的具體角度可能有誤差，但在一定程度上反映了斷層端部地殼應力與斷層和區域主應力之間夾角的變化關系。連續模型與非連續模型有相似的變化特征，只是前者在數量上比后者的小。

3.5 復合斷層對地殼應力的影響

在地質構造活動強烈的地區，復合斷層是很常見的，根據斷層的發育關系大致可歸納為2類：近似平行斷層和交叉斷層。平行斷層和交叉斷層對地殼應力的影響采用非連續模型進行計算和討論。

3.5.1 平行斷層的影響

考慮2條長、寬相同的平行斷層，改變其中一條斷層與另一條固定斷層的距離，計算固定斷層端部地殼應力的變化，如圖8所示。從圖8可以看出：最大主應力和剪切應力都隨斷層間距減小而逐漸增大，但是，當斷層間距較大時，應力的變化非常小，此時平行斷層對地殼應力的影響幾乎可以忽略。僅在斷層間距很近時，應力隨距離減小而快速增大。但是，無論是最大主應力還是剪切應力，應力的絕對改變量都并不是很大，表明平行斷層是地殼應力比較弱的影響因素。

圖7 斷層端部地殼應力與斷層和區域最大主應力之間夾角θ的關系Fig. 7 Relationship between crustal stress at end of faults and angle θ between fault-strike and regional maximum principal stress

圖8 非連續模型斷層端部地殼應力與平行斷層間距d的關系Fig.8 Relationship between crustal stress at end of faults and distance d between two parallel faults for discontinuous model

3.5.2 交叉斷層的影響

交叉斷層的宏觀形態多為“X”形，計算“X”形交叉斷層具有不同夾角(β)時的地殼應力，結果表明：除了斷層端部，斷層交叉處也是應力的集中區，這與孤立斷層的地殼應力狀態明顯不同；應力擾動范圍隨交叉斷層夾角的變化而發生改變；在斷層交叉處附近，位于區域最大主應力(σ1)方向上的 2個三角區的最大主應力表現為應力強化，而位于區域最小主應力(σ3)方向上的2個三角區的最大主應力表現為應力弱化；斷層交叉處及附近的剪切應力表現為高值且呈中心對稱形態。圖9顯示了地殼應力與交叉斷層夾角(β)的關系。從圖 9可以看出：交叉斷層夾角對地殼應力的影響很大。

圖9 非連續模型斷層交叉處地殼應力與斷層夾角β的關系Fig.9 Relationship between crustal stress at faults crossing point and intersecting angle β of two faults for discontinuous model

4 結論

(1) 非連續模型比連續模型能夠更好地模擬斷層，突出了斷層對地殼應力的影響作用，克服了單純使用軟弱帶模擬斷層所導致的應力削弱現象。

(2) 非連續模型通過接觸單元可以描述連續模型所不能反映的斷層帶與圍巖之間的不連續錯動。

(3) 斷層端部和復合斷層交叉處附近是應力集中區；斷層帶具有較高的剪切應力，但其自身及附近局部區域的最大主應力受到了一定的弱化。

(4) 斷層的彈性模量、斷層與區域主應力的夾角和交叉斷層是地殼應力的重要影響因素，平行斷層是弱影響因素，斷層泊松比對地殼應力的影響不大。

[1] 謝富仁, 陳群策, 崔效鋒, 等. 中國大陸地殼應力環境基礎數據庫[J]. 地球物理學進展, 2007, 22(1): 131?136.XIE Fu-ren, CHEN Qun-ce, CUI Xiao-feng, et al. Fundamental database of crustal stress environment in continental China[J].Progress in Geophysics, 2007, 22(1): 131?136.

[2] 張咸恭, 王思敬, 張倬元. 中國工程地質學[M]. 北京: 科學出版社, 2000: 97?103.ZHANG Xian-gong, WANG Si-jing, ZHANG Zhuo-yuan.Chinese engineering geology[M]. Beijing: Science Press, 2000:97?103.

[3] 蘇生瑞, 朱合華, 王士天, 等. 斷裂物理力學性質對其附近地應力場的影響[J]. 西北大學學報: 自然科學版, 2002, 32(6):655?658.SU Sheng-rui, ZHU He-hua, WANG Shi-tian, et al. The effect of fracture properties on stress field in the vicinity of a fracture[J].Journal of Northwest University: Natural Science Edition, 2002,32(6): 655?658.

[4] 朱維申, 阮彥晟, 李曉靜, 等. 斷層附近應力分布的異常和對隧洞穩定的影響[J]. 地下空間與工程學報, 2008, 4(4):685?689.ZHU Wei-shen, RUAN Yan-sheng, LI Xiao-jing, et al. Abnormal stress distribution adjacent to a fault and it’s influence on stability of tunnel[J]. Chinese Journal of Underground Space and Engineering, 2008, 4(4): 685?689.

[5] 孫禮健, 朱元清, 楊光亮, 等. 斷層端部及附近地應力場的數值模擬[J]. 大地測量與地球動力學, 2009, 29(2): 7?12.SUN Li-jian, ZHU Yuan-qing, YANG Guang-liang, et al.Numerical simulation of ground stress field at ends and vicinity of a fault[J]. Journal of Geodesy and Geodynamics, 2009, 29(2):7?12.

[6] 伍佑倫, 王元漢. 不連續面激活與巷道圍巖破壞區關系的探討[J]. 巖土力學, 2007, 28(6): 1197?1200.WU You-lun, WANG Yuan-han. Discussion on relationship between activation of discontinuities and disturbed zone in surrounding rock mass of tunnels[J]. Rock and Soil Mechanics,2007, 28(6): 1197?1200.

[7] 榮冠, 朱煥春, 王思敬. 錦屏一級水電站左岸邊坡深部裂縫成因初探[J]. 巖石力學與工程學報, 2008, 27(S1): 2855?2863.RONG Guan, ZHU Huan-chun, WANG Si-jing. Primary research on mechanism of deep fracture formation in left bank of Jinping first hydropower station[J]. Chinese Journal of Rock Mechanics and Engineering, 2008, 27(S1): 2855?2863.

[8] 晏長根, 劉彤, 伍法權. 復雜條件下大型地下洞室群的變形穩定性分析[J]. 工程地質學報, 2008, 16(1): 84?88.YAN Chang-gen, LIU Tong, WU Fa-quan. Deformation and stability analysis of large-scale underground cavity group under complicated geological conditions[J]. Journal of Engineering Geology, 2008, 16(1): 84?88.

[9] 周維垣, 楊強. 巖石力學數值計算方法[M]. 北京: 中國電力出版社, 2005: 282?284.ZHOU Wei-yuan, YANG Qiang. Numerical methods for rock mechanics[M]. Beijing: Chinese Electric Press, 2005: 282?284.

[10] 石根華. 塊體系統不連續變形數值分析新方法[M]. 北京: 科學出版社, 1993: 50?80.SHI Gen-hua. New numerical analysis methods for discontinuous deformation of block system[M]. Beijing: Science Press, 1993: 50?80.

[11] 黨亞民, 陳俊勇, 晁定波. 中國及其鄰區地殼應力場數值特征的研究[J]. 測繪科學, 2001, 26(2): 11?14.DANG Ya-min, CHEN Jun-yong, CHAO Ding-bo. Numerical study on stress characteristic in China[J]. Science of Surveying and Mapping, 2001, 26(2): 11?14.

[12] 鄭勇, 陳颙, 傅容珊, 等. 應用非連續性模型模擬斷層活動對青藏高原應力應變場的影響[J]. 地球物理學報, 2007, 50(5):1398?1408.ZHENG Yong, CHEN Yong, FU Rong-shan, et al. Simulation of the effects of faults movement on stress and deformation fields of Tibetan Plateau by discontinuous movement models[J].Chinese Journal of Geophysics, 2007, 50(5): 1398?1408.

[13] Bird P, Kong X H. Computer simulations of California tectonics confirm very low strength of major faults[J]. Geol Soc Am Bull,1994, 106: 159?174.

[14] 傅容珊, 黃建華. 地球動力學[M]. 北京: 高等教育出版社,2001: 80?90.FU Rong-shan, HUANG Jian-hua. Geodynamics[M]. Beijing:Higher Education Press, 2001: 80?90.

[15] ANSYS Inc. Release 11.0 documentation for ANSYS: Contact technology guide[R]. Canonsburg PA:ANSYS Inc, 2007: 1?10.

[16] Turcotte D L, Schubert G. Geodynamics[M]. Cambridge:Cambridge University Press, 2002: 433?436.