斷層錯動數值模擬中的幾個關鍵影響因素

張炳焜

(四川省建筑科學研究院，四川成都610081)

針對斷層錯動類問題的有限元數值模擬，有幾個關鍵問題值得注意，全面的考慮這些問題無疑能使計算結果更為準確。本項研究分別考慮了以下幾個影響因素：準靜態法、網格尺寸、應變軟化和材料阻尼。為了能夠說明這四個影響因素的重要性，本文建立簡單均一上覆砂土層二維有限元模型，深度H=20 m，寬度L=80 m，土層內摩擦角φ=35°，剪脹角為ψ=15°，彈性模量E=60 MPa，泊松比ν=0.3，不考慮內聚力的影響。以各影響因素為變量展開計算，模擬逆斷層60°，垂直位錯量d=1 m時上覆土層應變與位移變形情況。

1 準靜態分析

準靜態法是一種用靜力學方法近似解決動力學問題的簡易方法，該方法能在有限程度上反映荷載的動力特性，它物理概念清晰，與全面考慮結構物動力相互作用的分析方法相比，計算方法較為簡單，計算工作量很小、參數易于確定，并積累了豐富的使用經驗，易于設計工程師所接受。在進行數值模擬時，有兩個方面值得注意，一是在對地震導致的地面震動問題進行分析時，動力分析法是必不可少的。而另一方面，對于斷層錯動問題，若是綜合考慮斷層錯動產生的波動和位移這兩種因素的影響，則模擬結果會缺乏針對性。為了能夠有針對性的說明斷層錯動時上覆土層的響應情況，選擇準靜態法進行模擬。ABAQUS/EXPLICIT模塊常常在解決準靜態類型問題是非常有效的。

本文對斷層垂直錯動速度分別為v=1 m/s，v=0.5 m/s，v=0.1 m/s和v=0.05 m/s的情況進行了計算和對比，從上覆土層塑性應變云圖(圖1)可以發現，當v=1 m/s時，塑性應變帶未曾貫通上覆土層，而是較為集中的分布在了上覆土下部，由于高速的下部斷層錯動，導致上覆土層下部，即靠近斷層處的土體受到高速的撞擊而產生局部的應變，該區域的應變集中，消耗了本應傳遞至上覆土上部的應變能量。隨著錯動速度的降低，當v=0.1 m/s時，塑性應變帶較為完整的貫通了上覆土層。根據以上分析得出結論：相對高速的斷層錯動速度會導致上覆土下部應變的集中，并且地面傾斜值低于相對低速的情況。而當斷層錯動速度減小到一定程度后，若是繼續降低斷層錯動速度，覆土的塑性應變和地面位移將不再發生明顯的變化，反而增加了計算所需的時間。因此，本文推薦采用的斷層錯動速度為0.1 m/s。

(a)v=1 m/s

(b)v=0.5 m/s

(c)v=0.1 m/s

(d)v=0.05 m/s圖1 塑性應變云圖

2 網格尺寸

對有限元模擬計算而言，除了材料參數的選取以及模型的建立外，模型的網格劃分也是非常關鍵的，網格劃分的好壞直接影響計算的精度。斷層錯動情況的網格尺寸是個關鍵因素，與其他模擬不同，斷層錯動伴隨著較大的位移和應變的發生，若是模型網格尺寸較大，則可能導致塑性應變區的寬度發生變化，即破裂帶可能變寬，而影響了計算精度，對后續的研究造成誤導。本文有限元模型中部斷層區域是研究的重點，在該區域內的網格尺寸更應該進行嚴格的控制，力求達到均勻和精細。

本文對模型中部網格尺寸分別為L=10 %H，L=7.5 %H，L=5 %H和L=2.5 %H的情況進行了計算和對比，H為覆土深度。從上覆土層塑性應變云圖(圖2)可以發現，隨著網格尺寸的減小，塑性應變帶的寬度逐漸減小，當L=2.5 %H時，塑性應變帶寬度達到合理水平。因此，本文推薦網格尺寸L=2.5%H。

(a)L=10 %H

(b)L=7.5 %H

(c)L=5 %H

(d)L=2.5 %H圖2 塑性應變云圖

3 應變軟化

應變軟化是指材料試件經一次或多次加載和卸載后，進一步變形所需的應力比原來的要小，即應變后材料變軟的現象。在斷層發生錯動后，涉及到較大的上覆土體變形和應變，土體強度的降低是應該考慮的，若是不考慮此因素的影響而進行數值模擬，得出的模擬結果與實際可能會產生一定的差別，之前學者們對此得出了較好的對比結果[1-3](Roth et al.1981,1982;Loukidis 1999)，針對路堤邊坡穩定性而言，亦有類似結論[4-5](Potts et al.1990,1997)。針對摩爾庫倫模型，內摩擦角、剪脹角和內聚力是該模型最主要的強度參數，將上覆土層的塑性剪應變e作為控制參數，考慮內摩擦角、剪脹角和內聚力隨塑性剪應變的增加而數值逐漸降低，當塑性剪應變達到某一個臨界值后，內摩擦角、剪脹角和內聚力將不再發生變化，保持其殘余值。此處，內摩擦角和剪脹角的峰值分別記為φmax和ψmax，殘余值分別記為φres和ψres，塑性剪應變臨界值記為e0。則內摩擦角和剪脹角與塑性剪應變e之間的關系為：

(1)

(2)

本文分別對考慮應變軟化和未考慮應變軟化的情況進行數值模擬，當考慮應變軟化時，內摩擦角和剪脹角的殘余值分別為φres=25°，ψres=0°，塑性剪應變臨界值e0=0.1。

從以上兩種情況的塑性應變云圖可以看出(圖3)，二者的塑性應變帶均貫通了上覆土層，并且分布較為均勻，但是，當未考慮應變軟化時，破裂帶和水平面的夾角明顯較小。從某種意義上降低了斷層的傾角。為了能夠充分的估計斷層錯動對上覆土層的影響，考慮應變軟化模型是合理的，有必要的。

(a)考慮應變軟化

(b)未考慮應變軟化圖3 塑性應變云圖

4 材料阻尼

對于準靜態問題而言，雖然該方法大大的降低了動力波動對模型的影響，但是材料的阻尼仍然是不能忽視的問題，本文考慮Rayleigh阻尼的影響。

在有限元分析中，當采用振型分解法時，即把多自由度結構的復雜振動分解為按各個振型的獨立振動的疊加，則常假定阻尼矩陣[C]滿足振型的正交條件，即:

(3)

瑞利(Rayleigh)曾提出阻尼矩陣采用質量矩陣與剛度矩陣的線性組合，必滿足正交條件，即：

[C]=α[M]+β[K]

(4)

阻尼矩陣用轉置振型矩陣及振型矩陣前后相乘后，將得: 如令α=0，即C=β[K]，此時β=2ξn/ωn

ξn=(α/ωn+βωn)/2

(5)

如令β=0，即C=α[M]，此時α=2ξnωn

可見，當阻尼矩陣僅與質量矩陣有關時，阻尼比與振動頻率成反比，高階振型的阻尼就非常小；同樣，當阻尼矩陣僅與剛度矩陣有關時，阻尼比正比于頻率，高階振型的阻尼非常大，這就使結構的高階振型對結構動力反應的影響減弱甚至消除了。

系數α、β可由給定的第i、第j階振型的阻尼比ξi、ξj反算得到：

(6)

通常結構振動分析中廣泛采用的是簡單的Rayleigh阻尼模型。

在工程中常采用簡化求法，即：

α=ξ1·ω1β=ξ1/ω1

(7)

本文分別對考慮材料阻尼和未考慮材料阻尼的情況進行數值模擬，從塑性應變云圖可以看出(圖4)，當考慮材料阻尼時，塑性應變帶貫通了上覆土層，并且分布較為均勻，而未考慮材料阻尼時，破裂帶除了貫通上覆土層外，在破裂帶左側地面出現了一個局部的塑性應變區域，上覆土層破裂的不規律。因此，雖然采用準靜態法進行數值模擬，但是材料的阻尼仍是不容忽視的關鍵因素之一。

(a)考慮材料阻尼

(b)未考慮材料阻尼圖4 塑性應變云圖

5 結論

通過以上計算分析，得出以下結論：對于地震引發的地面震動問題，動力分析方法是切實有效的。不過對于下覆斷層錯動類問題，為了能夠較有針對性的對問題進行分析，建議采用準靜態法進行研究，并綜合考慮有限元網格尺寸、土體應變軟化特性以及材料阻尼等影響因素。避免對上覆土層及結構的影響估計不足，使分析結果與實際貼切，為日后的工程設計和研究提供幫助。

[1] Roth WH, Scott RF, Austin I (1981). Centrifugemodelling of fault propagation through alluvial soils.Geophys Res Lett 8(6):561-564

[2] Roth WH, Sweet J, Goodman RE (1982) Numerical and physical modelling of flexural Slip phenomena and potential for fault movement. Rock Mech 12:27-46 (Suppl.)

[3] Loukidis D (1999). Active fault propagation through soil. Diploma Thesis, School of Civil Engineering, National Technical University, Athens

[4] Potts DM, Dounias GT, Vaughan PR (1990). Finite element analysis of progressive failure of Carsington Embankment. Géotechnique 40(1):79-101

[5] Potts DM, Kovacevic N, Vaughan PR (1997). Delayed collapse of cut slopes in stiff clay. Géotechnique 47(5):953-982