基于反應譜分析的邊坡動力響應與地震穩定性數值模擬

梁健偉， 陳海嘯， 戴旭明， 張潤明， 曾偉雄

(1.中國能源建設集團廣東省電力設計研究院有限公司， 廣州 510663； 2.中國地質大學(武漢)工程學院， 武漢 430074； 3.中國科學院武漢巖土力學研究所, 武漢 430071)

中國是多山的國家，在降雨或工程開挖擾動作用下，土質或巖質山體極易出現坡體滑動從而導致滑坡工程地質災害頻發。分析滑坡災害演化過程以及不同作用環境下坡體的穩定性，對于工程施工環境以及生命財產的合理保護具有重要意義[1]。特別是中國東南沿海及西南山地處于板塊構造邊緣，受地殼運動影響，該區域地震頻發，因此，亟需開展工程擾動及地震共同作用下邊坡穩定性分析研究。

針對邊坡穩定性分析問題，目前主要研究的是靜力作用環境下的分析方法：林澤雨等[2]通過室內試驗研究了降雨及坡體組成對邊坡穩定性的影響規律；李龍起等[3]結合離散元數值分析軟件PFC2D對降雨條件下邊坡穩定性進行模擬分析；閆超等[4]則提出了一種基于強度折減法的滑動面分析方法，系統研究了次級滑動面的潛滑因素；劉興等[5]利用有限差分數值軟件FLAC3D對庫水升降影響條件下下邊坡穩定性進行了分析研究；宋彥琦等[6]基于強度折減理論，采用大型有限元軟件ABAQUS，以特征部位位移突變作為破壞判據，研究了邊坡坡度和凍融循環次數對邊坡穩定性的影響；史俊濤等[7]分析了容重增加法和強度折減法在邊坡進行穩定性分析中的應用，結果表明后者比前者更加合理。胡衛東等[8]研究了不同的非線性準則條件下的邊坡穩定性問題。上述研究通過理論、試驗及數值模擬方法系統討論了降雨環境、庫水位升降、凍融循環及分析方法等對邊坡穩定性的影響，但并未考慮動力響應對安全系數的弱化作用。

對于動力響應下邊坡穩定性研究：孔宇陽等[9]利用擬動力方法推導了考慮地震荷載及巖體放大系數等不同工況下巖石邊坡抗傾覆性的安全系數；張迎賓等[10]利用有限差分軟件FLAC3D，開展了考慮坡高、坡角及坡面形狀等地形因素的地震動力響應分析；段現花等[11]聯合有限元方法、極限平衡理論與Newmark法，分析了地震作用下錨碇邊坡的動力穩定性；胡訓健等[12]利用PFC2D研究了地震作用下斷續節理的巖橋長度及間距組合形式對層狀巖質邊坡的破壞模式、動力響應規律以及應力演化特征的影響規律。上述研究對不同因素下邊坡穩定性分析，均取得了較好研究成果。但是，大部分研究工作采用商業軟件自帶的分析模塊，難以表征地震動力特性對邊坡峰值加速度及強度的影響。

反應譜法作為主要應用于結構抗震分析領域一種動態分析方法[13]，目前較少應用到邊坡工程領域。基于反應譜理論，提出將反應譜法運用到邊坡穩定性分析和地震穩定性分析中。同時考慮了地震動力特性的反應譜法動態分析，通過地面地震動峰值加速度，研究邊坡地震動峰值加速度，進而進行強度折減，得到滑動面的安全系數。

1 反應譜分析

1.1 基本原理

反應譜法相較靜力法，能夠考慮在地震時結構和地面所表現的動力特性，標志著對地震分析的方法正式進入了動力研究的階段。該方法主要思路是：將結構自由度進行拆解，使得多自由度結構能夠簡化成多個單自由度結構組合，針對這些單自由度結構體，可以分別計算其在地震作用下的最大響應，然后，將上述最大響應按照合理的組合方法計算最終反映多自由的結構地震響應。

起初工程界常采用的模態組合方法是平方和平方根法(square root of the sum of square，SRSS)，此法假定了結構的各階振型反應之間互不影響，當結構有彼此接近的固有頻率時，該方法計算得到的最大響應值可能會出現偏大或是偏小的情況。后來，Wilson等[14]根據結構的隨機振動理論在 SRSS 法的基礎上研究給出了完全二次組合法(complete quadratic combination，CQC)，即常規的模態分解反應譜法。經過中外大量的學者研究，反應譜的振型組合方法有多種，但就目前而言，最常用方法的是SRSS法和CQC法。

(1)SRSS法計算方程為

(1)

式(1)中：F為結構的地震作用效應；Si為結構第i階振型的地震作用效應；n為振型階數。

SRSS法對于計算主振型的周期差別較大結構的抗震具有較高的精度，適用于地震作用效應中單一方向問題。

(2)CQC法計算方程為

(2)

(3)

(4)

式中：式(2)為 CQC 法適用的條件方程；ρT為周期比；Tj和Ti(Tj≥Ti)為相鄰兩階振型的周期；ξ為阻尼比；Si和Sj分別為第i階和第j階的地震作用效應；rij為振型組合相關系數。

1.2 基于反應譜法動態分析計算流程

采用反應譜法分析邊坡的地震影響系數，通過以地震安全性評價所獲得的地震作用下動峰值加速度為基本參數，研究上述參數不同放大系數分布特征,從而獲得邊坡的振型，在此基礎上，計算不同位置邊坡的地震影響系數與地震慣性力，然后將地震慣性力疊加到各個邊坡的單元體中進行穩定性分析,并可計算邊坡在地震動力下的模態和應力,根據強度折減法并可獲得邊坡穩定系數,這是一種更為先進的準動態分析方法。

計算流程分為兩個部分：①基于ANSYS的彈性有限元的模態分析和反應譜分析，求得邊坡各巖土體單元的節點地震響應加速度譜值，計算單元的地震影響系數；②基于FLAC3D的疊加ANSYS反應譜計算成果的彈塑性非線性大變形和破壞分析，導入ANSYS反應譜計算的各單元地震慣性力，考慮開挖+重力+水荷載，進行強度折減計算，搜索滑動面并計算最小安全系數。具體計算流程如下。

(1)建立既有點單元方法的模態分析模型,計算目標邊坡的第1～10階振型的特性。

(2)開展反應譜分析研究，反應譜分析如下。①對于邊坡的每一個模態，計算在不同地震激勵方向的模態參與因子與模態系數;②根據模態系數，計算不同模態最大響應下的位移;③采用均方根法，計算多自由度結構的整體地震響應。

大寫字母為核電廠設備代碼圖1 邊坡分布示意圖Fig.1 Schematic diagram of slope distribution

根據核電廠抗震設計規范要求，同時考慮水平和垂直地震加速度的共同作用，輸出水平向和豎直向的模態參與因子和模態系數，計算水平方向和豎直方向的邊坡各單元的節點最大地震響應加速度值，換算為各單元的地震慣性力，繪制地震慣性力云圖，分析其最大動力響應特征及變化規律。

(3)將ANSYS軟件計算的各單元最大地震慣性力輸入FLAC3D數值模型中，作為SL-2的地震荷載力，采用彈塑性模型，考慮開挖+自重應力場+水荷載，計算邊坡穩定系數和位移。

2 邊坡動力響應與地震穩定性分析

2.1 工程地質概況

某核電工程廠址位于廣東省惠州市，初步計劃場地平整工作分兩階段實施：第一階段的場地平整范圍包括一期、二期的工程場地(以下統稱為一期工程場地)，第二階段的場地平整范圍為三期的工程場地。某自然邊坡位于一期核電廠區西側，總長550 m，最大高度106 m，坡腳與核島區最近的直線距離86 m，小1.4倍的邊坡高度(圖1)，根據《核電廠抗震設計規范》(GB 50267—2019)，該邊坡屬于核安全相關的自然斜坡。

該邊坡總體呈北西走向，傾向北東，植被較發育，坡度為20°～30°，覆蓋層厚度相對較厚,一般5～10 m,局部地段15 m，以強風化巖為主；邊坡中段局部坡度達40°，覆蓋層一般以基巖為主，厚度一般2～5 m。其可能發生的破壞模式為上部土體沿巖土界面發生的滑移破壞。根據巖土層室內物理力學試驗：粉質黏土的天然黏聚力和內摩擦角分別為17.8 kPa和27.9°，強風化熔結凝灰巖的天然黏聚力和內摩擦角分別為22.3 kPa和30.8°，中等風化熔結凝灰巖的顆粒密度為2.72 g/cm3，含水率為0.59%，微風化熔結凝灰巖的顆粒密度為2.74 g/cm3，含水率為0.16%，中等風化熔結凝灰巖的單軸抗壓強度為59.8 MPa，微風化熔結凝灰巖的單軸抗壓強度為132.5 MPa。根據勘察結果，選取典型剖面(26-26)，進行邊坡的穩定性分析。

廠區內發育斷裂構造7條，主要走向以北東為主，這些斷裂零散出露，對該邊坡的直接影響不大。如圖2所示，在該自然邊坡坡腳處發育一規模較小不穩定邊坡(編號HB01)，坡體主要巖土層為坡積粉質黏土，混部分碎石，滑面下部為強風化熔結凝灰巖。

圖2 HB01滑坡體Fig.2 HB01slope

2.2 邊坡動力響應計算模型與結果

反應譜法采用彈性有限元法計算，計算參數采用動態參數，根據一期邊坡勘察結果，考慮到廠址區雨量充沛，覆蓋土層及強風化巖體物理力學參數選用飽和狀態下的推薦值，中等風化、微風化巖石的物理力學參數參考巖石試驗結果，巖土體動彈性參數如表1所示。

表1 采用動態參數的本構模型的計算參數

采用大型商業有限元計算軟件ANSYS進行SL-2 地震作用下的邊坡模態分析和反應譜分析，獲得各節點最大地震響應。

采用ANSYS建立的邊坡分析模型有限元單元共6 959個，節點共7 184個，主要采用四邊形單元,網格質量較好,網格單元、節點和材料分組信息文件如圖3所示。

圖3 研究剖面網格劃分示意圖Fig.3 Schematic diagram of grid division of research section

邊坡前10階振動特性(固有頻率、模態參與因子)及譜值如表2所示。邊坡在SL-2地震作用下的第1、5、10階振型圖如圖4所示。由表2、圖4可知，當振型由第1階逐漸增加到第10階時，震動頻率逐漸由3.89 Hz增加到9.56 Hz，振型也發生了顯著變化。在第1階時，振型矢量主要呈現水平位移特征，形式比較平緩，沒有出現較大激震現象；在第5階時，振型矢量發生了一定偏轉，出現了垂直坡面振型情況，最大振型出現在了坡體中部；在第10階振型時，坡頂和坡腳均出現了較大激震，坡頂矢量方向基本沿坡面分布，在坡腳則在震動作用下出現了大幅的垂直坡面的振型。

表2 邊坡前 10 階振動特性

圖4 邊坡前10階典型振型矢量圖Fig.4 Typical vector figures of the first 10 vibration modes of slope

在有限單元方法中，所謂模態分析，實際上是求單元或者剛度矩陣的特征值，特征值對應的數量即該方法中的階數，階次則是表示上述特征值由大到小的排序。一般來說，實際工程問題模態階數無窮，但考慮到對研究對象起主導作用的階數情況，一般根據經驗取前面部分階數，計算過程中，只計算了前10階。然后采用SRSS法疊加各振型后的邊坡最大加速度及剪應力如圖5所示。

圖5 SRSS疊加后的邊坡加速度及應力圖Fig.5 Acceleration and stress diagram of slope after SRSS superposition

通過分析可知，邊坡坡面在SL-2地震作用下產生了放大效應，并以坡頂和坡中表層的強風化層受到的激勵最大，最大合加速度矢量方向順坡面向上，分析剪應力云圖可知，強風化層頂面及底面為應力集中區，為地震潛在破壞區域。根據反應譜分析方法的計算過程可以計算出地震影響系數的最大值為1.12。剖面地震響應統計如表3所示。

表3 26-26剖面地震響應分析極值

2.3 地震穩定性分析數值模型與結果

考慮采用的動態有限元計算方法(反應譜法)僅考慮了彈性計算，沒有考慮塑性變形和破壞準則，因此，需補充采用FLAC3D軟件計算彈塑性變形下，疊加地震荷載加速度后的邊坡穩定性的計算。

基于工程地質概況中巖體力學參數復核結果，采用如表4所示參數進行數值模擬研究。

表4 計算模型與參數

2.3.1 自然邊坡的穩定

由圖6可以看出：最小主應力的分布區間為-1.518～0.032 MPa,該區間的上限為正值，且正值是分布在滑坡體表面的，意味著在滑坡體表面存在拉應力；在滑坡體表面以下最小主應力均為負值，說明該滑坡體內部受壓。

由圖7可以看出：最大主應力的分布區間為-5.44～-0.016 MPa,該區間的上限和下限均為負值，分布在整個滑坡體，意味著整個滑坡體受壓應力。

針對研究剖面，施加自重及水位，利用表4提供的巖體力學參數計算至平衡狀態，得到圖6、圖7所示最大、最小主應力分布，結果表明自重情況下研究剖面可能在淺地表附近產生30 kPa左右的拉應力。在此基礎上采用FLAC3D自帶的強度折減法，計算得自然工況安全系數為1.85，滑動面分布如圖8所示。

圖6 最小主應力分布云圖Fig.6 Cloud diagram minimum principal stress distribution

圖7 最大主應力分布云圖Fig.7 Cloud diagram maximum principal stress distribution

圖8 自然工況下邊坡滑動面分布示意圖Fig.8 Distribution diagram of slope sliding surface distribution under natural conditions

2.3.2 地震作用下邊坡穩定性

在自重+水位平衡后，分別施加反應譜計算結果轉化而來的地震力，然后分析不同剖面的變形、塑性區、穩定性。計算得到地震工況下位移、最小主應力和最大主應力分布云圖如圖9～圖11所示。針對26-26研究剖面，圖9結果表明地震作用下該剖面可能產生約7.2 cm的塑性變形。

圖9 地震工況下位移分布云圖Fig.9 Cloud diagram of displacement distribution under earthquake condition

由圖10可以看出：最小主應力的分布區間為-2.08～1.12 MPa,該區間的上限為正值，且正值是分布在滑坡體內部的，意味著在滑坡體內部存在拉應力；在滑坡體其他位置最小主應力均為負值，說明該滑坡體大部分受壓。

圖10 最小主應力分布云圖Fig.10 Cloud diagram of minimum principal stress distribution

由圖11可以看出：最大主應力的分布區間為-8.09～-0.007 5 MPa, 該區間的上限和下限均為負值，分布在整個滑坡體，意味著整個滑坡體受壓應力。

圖11 最大主應力分布云圖Fig.11 Cloud diagram of maximum principal stress distribution

在此基礎上采用FLAC3D自帶的強度折減法，計算得自然工況安全系數為1.21，折減后滑動面和塑性區分布如圖12所示，塑性區主要位于強風化巖體、坡度30°～40°范圍，應引起注意。

地震工況條件下26-26剖面塑性區的分布范圍均比自然工況條件下的范圍分布更廣，地震工況條件下26-26剖面的塑性區主要分布在坡面。

針對不同工況進行計算，得出剖面的自然工況穩定性、地震工況穩定性安全系數，分析可知自然和地震工況條件下研究剖面的安全系數分別為1.85和1.21，與估算的邊坡安全系數相比較，所以該剖面安全。根據研究中有關主要巖土層物理力學試驗結果，對邊坡穩定性進行了試算，按照地質提供參數估算的邊坡安全系數自然狀態下，邊坡的安全系數在1.426～2.056，按照核安全邊坡控制剖面的地震安全系數為0.907，計算結果偏小。

none表示未破壞；shear-n表示現在剪切破壞；shear-p 表示現 在彈性狀態，但是過去剪切破壞；tension-n表示現在張拉破 壞；tension-p表示現在彈性狀態，但是過去張拉破壞圖12 地震工況下邊坡滑動面及塑性區分布Fig.12 Distribution of slip surface and plastic zone of slope under earthquake condition

3 結論

基于反應譜理論，提出把反應譜法運用到邊坡穩定性分析和地震穩定性分析當中去。該方法同時考慮了地震動力特性的反應譜法動態分析，通過地面地震動峰值加速度，研究邊坡地震動峰值加速度，進而進行強度折減，得到滑動面的安全系數。通過研究，得出如下結論。

(1)采用反應譜方法可考慮邊坡不同位置的地震影響系數的變化規律，符合邊坡地震動響應的實際規律。

(2)在自然工況、降雨工況能夠保持整體穩定，滿足規范安全性要求，剛體極限平衡法地震工況按核安全相關邊坡與SL-2級地震峰值驗算研究剖面穩定達到規范1.50要求，動力法地震工況滿足1.20要求。

(3)局部穩定按照核安全相關邊坡考慮地震，安全系數為1.507，滿足規范要求(靜力核安全地震工況1.50)。按照SL-2級地震峰值考慮地震(水平0.25，垂直0.20)考慮，研究剖面安全系數為1.592滿足規范要求(靜力核安全地震工況1.50)。