張衛(wèi)杰 余瑞華 宋 健 姬 建

(①河海大學(xué)巖土力學(xué)與堤壩工程教育部重點(diǎn)實(shí)驗(yàn)室，南京 210024，中國)(②河海大學(xué)江蘇省巖土工程技術(shù)工程研究中心，南京 210024，中國)

0 引 言

滑坡是地震次生災(zāi)害中非常重要的一部分，其較遠(yuǎn)的滑動距離和較大的沖擊力，給人民群眾的生命、財產(chǎn)安全帶來了巨大威脅，如觸發(fā)滑坡最多的2008年汶川地震，有197 481處滑坡記錄(許沖等，2019)，其中大光包滑坡面積達(dá)7.2 km2、體積約1.2×108m3，最大堆積厚度達(dá)550 m(黃潤秋等，2014)；另一處紅石溝滑坡最大水平距離達(dá)2700 m、最大高差920 m，估算滑坡體積為1.5×107m3(崔圣華等，2019)。可見，由地震引發(fā)的滑坡往往具有速度快、滑程遠(yuǎn)、沖擊力大的特點(diǎn)(Bao et al.，2020a，2020b)，表現(xiàn)出大變形破壞的特征。從已有研究中分析可知，地震誘發(fā)滑坡的體積(Li et al.，2020)和位移(Rathje et al.，2011)是評估土坡抗震性能和滑坡災(zāi)害后果的重要指標(biāo)。另外，滑坡體的致災(zāi)范圍是判斷滑坡能否對已有建構(gòu)筑物造成損失及確定預(yù)警疏散范圍的重要依據(jù)，因此對地震邊坡破壞后的滑坡體大小和致災(zāi)范圍進(jìn)行研究具有重要意義。

地震作用下邊坡的安全評價一直是巖土工程領(lǐng)域的研究熱點(diǎn)之一(Hazari et al.，2020)。對于不同地震動力作用下的邊坡穩(wěn)定性問題，國內(nèi)外學(xué)者的研究主要集中在：(1)不同地震動參數(shù)下邊坡動力響應(yīng)的規(guī)律(張江偉等，2017，2018；羅洋等，2019)；(2)不同地震動幅值作用下邊坡的動力響應(yīng)規(guī)律(Shamy et al.，2019；Song et al.，2019；韓雪等，2020)；(3)不同方向地震作用下邊坡的動力響應(yīng)(侯紅娟等，2013)等。這些研究多數(shù)利用數(shù)值模擬或振動臺試驗(yàn)研究不同地震作用下邊坡的動力響應(yīng)規(guī)律，但由于邊坡大變形動力分析的困難性，地震邊坡失穩(wěn)后滑坡體大小和致災(zāi)范圍研究依然少見。

近年來，越來越多的研究將彈塑性本構(gòu)模型引入光滑粒子流體動力學(xué)(SPH)方法(胡嫚等，2019；Peng et al.，2019)，以更好地模擬土體在大變形狀態(tài)下的應(yīng)力-應(yīng)變行為，揭示更為準(zhǔn)確的土體大變形破壞規(guī)律。在地震邊坡分析方面，Hiraoka et al.(2013)采用SPH方法模擬了邊坡在地震作用下的動力響應(yīng)；Bao et al.(2020a，2020b)提出了基于統(tǒng)一模型的SPH方法，模擬了地震滑坡的啟動過程和滑動過程。盡管如此，地震邊坡失穩(wěn)后滑坡體大小和致災(zāi)范圍方面的研究，特別是邊坡坡角和地震動加速度幅值的影響規(guī)律研究，仍然存在不足。

本文基于SPH分析方法，采用基于Drucker-Prager破壞準(zhǔn)則的彈塑性本構(gòu)模型描述邊坡土體的應(yīng)力-應(yīng)變行為，設(shè)置不同類型的粒子實(shí)現(xiàn)動力荷載的施加和自由場邊界的模擬。利用SPH動力分析方法對已有文獻(xiàn)中的振動臺土坡試驗(yàn)進(jìn)行了模擬，驗(yàn)證了該模型的準(zhǔn)確性和精度。最后，通過對比4種坡度(30°、40°、50°及60°)土坡在不同地震動幅值作用下的SPH動力分析結(jié)果，揭示了坡角和地震動加速度幅值對邊坡滑坡體大小和致災(zāi)范圍的影響規(guī)律。

1 SPH動力分析方法原理

SPH的基本思想是將空間中連續(xù)的實(shí)體離散成一系列粒子，質(zhì)量、速度、應(yīng)力、變形等所有信息由這些粒子承載，粒子之間沒有任何鏈接。在整個分析過程中，SPH方法需要追蹤每個粒子在每個時刻的運(yùn)動信息。其無網(wǎng)格性質(zhì)及粒子間相互作用的特點(diǎn)使其更易處理大變形問題，消除了傳統(tǒng)拉格朗日方法中網(wǎng)格畸變和扭曲等問題(Huang et al.，2014)。

SPH方法的核心問題包括光滑核函數(shù)的光滑近似和粒子近似。其中光滑近似可表示為如下的形式：

(1)

式中：W為光滑核函數(shù)；h為光滑長度；x為粒子的坐標(biāo)信息。

函數(shù)及其導(dǎo)數(shù)的光滑粒子近似可以表示為：

(2)

(3)

式中：m為質(zhì)量；ρ為密度；j為粒子指標(biāo)。

本研究選擇3次樣條B函數(shù)作為光滑核函數(shù)(Huang et al.，2014)，其表達(dá)式為：

(4)

式中：R=|x-xj|/h，αd=15/(7πh3)。

SPH方法的控制方程包括連續(xù)性方程、動量守恒方程和能量守恒方程。在不考慮滑坡體溫度變化的情況下，可以忽略能量守恒方程，只考慮連續(xù)性方程和動量守恒方程。根據(jù)質(zhì)量守恒定律，連續(xù)方程的SPH形式如下，

(5)

式中：Wij為土粒子j在土粒子i處的光滑核函數(shù)值；v為速度；m表示方向；i和j為粒子指標(biāo)。為了獲得更加穩(wěn)定的密度計算，可以對式(5)進(jìn)行正則化修正，得到：

(6)

另一個控制方程是動量方程，由牛頓第二定律推導(dǎo)得出，如下所示：

(7)

式中：σ是土粒子的應(yīng)力張量；m、n是坐標(biāo)方向指標(biāo)；δmn是狄拉克函數(shù)；人工黏度項(xiàng)Πij被用來防止粒子的非物理穿透，其計算方法在文獻(xiàn)(Zhang et al.，2016)中列出。

土體的應(yīng)力-應(yīng)變行為需要通過具體的本構(gòu)模型進(jìn)行描述，在本研究中采用基于Drucker-Prager破壞準(zhǔn)則的彈塑性本構(gòu)模型描述邊坡土體的應(yīng)力-應(yīng)變行為。具有非關(guān)聯(lián)流動規(guī)則的基于Drucker-Prager破壞準(zhǔn)則的彈塑性本構(gòu)模型在Bui et al.(2007，2013)和McDougall et al.(2005)的工作中有詳細(xì)介紹，根據(jù)他們的工作，該模型的增量形式為：

(8)

式中：λe和μe是拉梅常量，可以通過彈性模量E和泊松比υ計算得到；dεij是總應(yīng)變增量，其他組成部分的計算方法如下：

(9)

(10)

C=λefkkgll+2μefklgkl+gkk

(11)

f是可以在Bui et al.(2007)的工作中找到的屈服函數(shù)：

(12)

(13)

g是塑性勢函數(shù)：

(14)

式中：ψ是剪脹角；常量αψ由下式計算：

(15)

如果摩擦角φ和剪脹角ψ不同，為非聯(lián)合流動法則。如果摩擦角φ和剪脹角ψ相同，則是聯(lián)合流動法則。由以上可知，fij和gij可按下式計算：

(16)

(17)

研究方法中，邊界是由若干層虛擬粒子組成的，與運(yùn)動粒子(土粒子)具有不同的粒子類型。在處理邊界的過程中，假設(shè)邊界粒子具有一個虛擬速度，能夠用于運(yùn)動粒子速度梯度的加權(quán)平均。邊界粒子對運(yùn)動粒子的影響是根據(jù)兩者之間的相對距離來確定的。在地震動荷載施加方面，通過對底部固定邊界粒子施加速度，從而在粒子層面上實(shí)現(xiàn)了動力荷載的施加。本研究底部采用Hiraoka et al.(2013)的具有地震荷載的無滑移邊界，來實(shí)現(xiàn)動力荷載的施加；具體可描述為:

VAB=VA-VB=β(VA-Vseismic)

(18)

VB=(1-β)VA+βVseismic

(19)

式中：VA、VB和Vseismic分別為土粒子速度、邊界粒子速度和地震波速度。其中：β=min(βmax，1+dB/dA)，是與運(yùn)動粒子和邊界粒子之間距離有關(guān)的參數(shù)。

基于Bao et al.(2020a)的研究，本文在左右邊界處設(shè)置自由場邊界粒子，以減少地震波的反射。在數(shù)值模擬中，自由場粒子被迫移動，通過邊坡粒子和自由場粒子在地震和重力作用下的同步計算，將自由場粒子的不平衡力施加于邊坡粒子的橫向邊界，使得計算區(qū)域中真實(shí)粒子產(chǎn)生的向外波被適當(dāng)?shù)匚?，以滿足橫向邊界上的位移和應(yīng)力條件，從而實(shí)現(xiàn)自由場邊界的粒子化建模。

2 SPH動力分析方法的驗(yàn)證

為了驗(yàn)證本文SPH動力分析方法在地震土坡模擬上的準(zhǔn)確性，本研究將其應(yīng)用在黃棟等(2017)的土坡振動臺試驗(yàn)?zāi)M中，通過對比地震邊坡的最終變形形態(tài)來驗(yàn)證本文方法的有效性。該土坡的尺寸如圖1所示。計算模型中粒子總數(shù)目為6977，其中邊界粒子有2118個，土坡土粒子有4857個，粒子初始間距為0.01 m。左邊界設(shè)置為自由場邊界，在荷載作用下發(fā)生運(yùn)動，計算區(qū)域內(nèi)運(yùn)動粒子產(chǎn)生的向外波被適當(dāng)?shù)匚找詼p少地震波的反射。初始狀態(tài)下，土坡粒子靜止不動，計算開始后可以在重力作用下運(yùn)動。在初始的計算步中，利用彈性模型對土坡的自重應(yīng)力進(jìn)行初始化，此后在地震動施加過程中，采用彈塑性模型描述土體的應(yīng)力-應(yīng)變行為。SPH模擬的單位時間步長為2.5×10-5s，計算總步數(shù)為1.21×107步，故穩(wěn)定分析的總時間為302.5 s，其中模型自重平衡的時間設(shè)置為1 s，地震動施加時間為300 s。另外，在地震動施加完畢之后，邊坡模型仍然具有速度，需要1.5 s的時間讓邊坡恢復(fù)靜止?fàn)顟B(tài)。SPH中使用的具體計算參數(shù)可見表1，本構(gòu)模型參數(shù)如表2所示。

圖1 振動臺土坡模擬尺寸圖

表1 振動臺土坡模型SPH模擬的參數(shù)

表2 振動臺土坡模型SPH模擬中的本構(gòu)參數(shù)

本案例施加的地震波與黃棟等(2017)論文中振動臺施加的地震波相同，即分別按照100 Gal(2 Hz、5 Hz、10 Hz)，300 Gal(2 Hz、5 Hz、10 Hz)，500 Gal(2 Hz、5 Hz、10 Hz)，700 Gal(10 Hz)逐級加大加速度和頻率進(jìn)行激勵。每級振動時間為30 s，每級振動荷載施加完畢后暫停一段時間。

圖2展示了輸入地震波強(qiáng)度為 500～700 Gal 時SPH模擬得到的土坡變形破壞形態(tài)。黃棟等(2017)將振動臺土坡變形破壞過程分為4個階段：(1)土坡振動固結(jié)和前緣堆積階段；(2)裂縫擴(kuò)展貫通階段；(3)整體滑動階段；(4)堆積體穩(wěn)定階段。本文利用SPH動力分析方法得到的土坡破壞呈現(xiàn)出與黃棟等(2017)類似的變形破壞過程，如圖2所示。土坡在運(yùn)行至地震波500 Gal 2 Hz施加完畢(即844萬步)時，土坡有明顯粒子滑落，于土坡前緣堆積；在地震波500 Gal 5 Hz施加完畢(即964萬步)時，土坡出現(xiàn)貫穿坡頂至坡腳的滑動面；在地震波500 Gal 10 Hz施加完畢(即1084萬步)時，土坡繼續(xù)滑動，但滑動距離不大；在地震波700 Gal 10 Hz施加完畢(即1204萬步)時，土坡的滑動區(qū)域明顯增大。圖3展示了SPH數(shù)值模擬的最終破壞形態(tài)與振動臺試驗(yàn)的最終破壞形態(tài)對比，兩者最終坡形基本吻合，且土體的最大滑動距離在數(shù)值上非常接近。由此可見，本文中提出的SPH動力分析方法在地震土坡變形分析上具有較高的精度和準(zhǔn)確性。

圖2 SPH模擬的試驗(yàn)土坡變形發(fā)展過程圖(單位：m)

圖3 土坡最終形態(tài)對比圖(單位：m)

另外，本研究中在自由場邊界外設(shè)置了振動邊界(圖1)，從振動臺試驗(yàn)邊坡的模擬結(jié)果來看，左側(cè)振動邊界粒子對邊坡動力響應(yīng)的影響很微小，證明了本研究振動邊界及自由場邊界設(shè)置的可行性。

3 不同地震動加速度幅值的影響分析

3.1 模擬案例

為了分析不同坡度土坡在不同地震動加速度幅值作用下的滑坡體大小和致災(zāi)范圍，研究建立了4種坡度的計算模型，分別是30°、40°、50°及60°。具體的尺寸如圖4所示，振動邊界粒子與自由場邊界粒子的設(shè)置也展示在圖4中。

圖4 不同土坡尺寸圖

計算模型的粒子總數(shù)目隨土坡坡度而改變，粒子初始間距統(tǒng)一設(shè)置為0.20 m。在初始狀態(tài)下，土坡粒子靜止不動，計算開始后可以在重力作用及地震作用下運(yùn)動。同樣地，在初始的計算步中，利用彈性模型對土坡的自重應(yīng)力進(jìn)行初始化，此后在地震動施加的過程中，采用彈塑性模型描述土體的應(yīng)力-應(yīng)變行為。SPH模擬的單位時間步長為5.0×10-5s，計算總步數(shù)為6.4×105步，分析的總時間是32 s，其中地震動施加時間為30 s。SPH中使用的具體計算參數(shù)見表3，本構(gòu)模型參數(shù)如表4所示。

表3 土坡模型SPH模擬的參數(shù)

表4 土坡模型SPH模擬中的本構(gòu)模型參數(shù)

本文的數(shù)值模擬首先需要在重力作用下實(shí)現(xiàn)初始地應(yīng)力平衡，因此要求重力作用下的邊坡是穩(wěn)定的，即安全系數(shù)≥1.0。若邊坡在自重應(yīng)力作用下的安全系數(shù)<1.0，自重應(yīng)力下邊坡無法保持穩(wěn)定，后續(xù)就無法考察地震動作用對邊坡穩(wěn)定及失穩(wěn)滑動的影響。另外，本研究通過FLAC2D計算得到的各邊坡安全系數(shù)均>1.0。

圖5 KOBE波加速度時程圖

本案例一共模擬了12個工況，具體設(shè)置如表5所示。基于OpenMP并行計算框架，在配備雙路Intel Xeon E5,2520V4處理器、64 GB DDR4內(nèi)存和Ubuntu Linux 20.04操作系統(tǒng)的工作站上設(shè)置了16線程，每個工況消耗時間約100 min。

表5 SPH模擬土坡工況表

3.2 結(jié)果分析

由于篇幅限制，這里選定40°和60°土坡在不同地震動幅值作用下的總位移云圖進(jìn)行比較，分別如圖6和圖7所示?？梢钥闯觯谕槐壤呦?，40°土坡在1.0倍地震動加速度幅值作用下只在坡頂處有較小變形；而在地震動幅值2.0倍和3.0倍作用下土坡整體有較大變形，出現(xiàn)了清晰且連續(xù)的位移較大區(qū)域，代表著該區(qū)域可能發(fā)生滑動，即滑坡體的范圍，其面積可以用來衡量土坡失穩(wěn)后滑坡體的大小。60°土坡在1.0倍地震動加速度幅值作用下就出現(xiàn)了明顯的位移較大區(qū)域，表示出現(xiàn)了可能的滑坡體。另外，同一土坡在不同地震動幅值作用下存在相似的滑動面，且其潛在的滑動區(qū)域隨著地震動幅值的增大而增加。從圖例中也可以明顯地看出，60°土坡存在更大的位移峰值。由此可知，對于具有相同參數(shù)的土坡，坡度越大穩(wěn)定性越差。

圖6 邊坡角為40°的土坡在不同地震動幅值作用下的總位移云圖(單位：m)

圖7 邊坡角為60°的土坡在不同地震動幅值作用下的總位移云圖(單位：m)

取土坡各工況下滑坡體區(qū)域的最大水平位移值作為地震土坡的滑動距離，可以表征地震土坡失穩(wěn)后的致災(zāi)范圍。圖8所示的是土坡在不同地震動幅值作用下最大水平位移的變化柱狀圖?？梢?，隨著土坡坡度增大，不同地震動幅值作用下土坡的最大水平位移均呈非線性增加趨勢。

圖8 不同地震動幅值作用下土坡最大水平位移變化柱狀圖

目前常用剪應(yīng)變來確定土坡的潛在滑移面，Li et al.(2020)驗(yàn)證了利用位移閾值確定土坡滑動面的準(zhǔn)確性。結(jié)合本案例，比較了40°和60°土坡在2.0倍地震動幅值作用下由剪應(yīng)變和總位移確定的潛在滑動面，如圖9所示?？梢?，最大剪應(yīng)變確定的土坡潛在滑動面與總位移確定的潛在滑動面具有一定的吻合度，滑動面均通過坡腳，且兩者界定的滑體大小也較為相近。因此，利用位移閾值確定土坡滑動面的方法具有一定的準(zhǔn)確性。

圖9所示的坡體在邊界處有些異常。原因是自由場粒子在吸收了地震波后產(chǎn)生速度和位移，其位移與土體位移存在微小差異。邊坡邊界處的異常與邊坡剪應(yīng)變和位移相比較為微小，對模擬結(jié)果的影響也較微小。這在振動臺試驗(yàn)邊坡的模擬結(jié)果中得到了體現(xiàn)。

圖9 在2.0倍地震動幅值作用下40°和60°土坡剪應(yīng)變和總位移云圖(單位：m)

結(jié)合本案例土坡選取的參數(shù)及土坡的最大總位移值，研究選取了不同的位移閾值來確定土坡的滑動面。之后，通過將總位移云圖簡化，得到了更為清晰的滑坡體分布圖，進(jìn)而為計算各土坡的滑坡體大小(以滑動區(qū)域的面積代替)提供了便利。由于1.0倍地震動幅值作用下各土坡的位移值較小，未出現(xiàn)較為清晰的滑動面，這里只計算和展示了2.0倍和3.0倍地震動幅值作用下的數(shù)據(jù)結(jié)果。圖10是以0.60 m為位移閾值的各土坡在2.0倍地震動幅值作用下的滑坡體分布圖。圖11是以0.60 m為位移閾值的各土坡在3.0倍地震動幅值作用下的滑坡體分布圖。

圖10 以0.60 m為位移閾值各土坡在2.0倍地震動幅值作用下滑坡體分布圖

圖11 以0.60 m為位移閾值各土坡在3.0倍地震動幅值作用下滑坡體分布圖

通過觀察滑坡體分布圖，可知，土坡的滑坡體大小與土坡坡度并不是簡單增加的關(guān)系，甚至土坡坡度越大，其滑體區(qū)域可能會越小。之后通過具體計算土坡在不同位移閾值之下的滑坡體區(qū)域，得到了具體的滑坡體數(shù)據(jù)，由此繪出了如圖12和圖13所示的土坡在2.0倍地震動幅值作用和3.0倍地震動幅值作用下的滑坡體大小隨坡度變化曲線?？梢?，同一位移閾值下各土坡滑坡體大小隨坡度的增大近似線性減小，且2.0倍地震動幅值作用和3.0倍地震動幅值作用下的各土坡滑坡體大小的變化規(guī)律相同。

圖12 2.0倍地震動幅值作用下不同位移閾值土坡滑坡體大小變化曲線

圖13 3.0倍地震動幅值作用下不同位移閾值土坡滑坡體大小變化曲線

綜上所述，隨坡度的增加，各土坡在地震作用下的穩(wěn)定性下降，致災(zāi)范圍變大；但滑坡體的區(qū)域面積會減小。

4 結(jié) 論

本文基于SPH動力分析方法，采用基于Drucker-Prager破壞準(zhǔn)則的彈塑性本構(gòu)模型作為土體的彈塑性本構(gòu)模型，設(shè)置不同類型的粒子實(shí)現(xiàn)了動力荷載的施加和自由場邊界的模擬。利用SPH動力分析方法對已有文獻(xiàn)中的土坡破壞振動臺試驗(yàn)進(jìn)行了模擬，對比兩者的最終形態(tài)，驗(yàn)證了該模型的準(zhǔn)確性和精度。通過對比4種坡度土坡在不同地震動幅值作用下的SPH動力分析結(jié)果，揭示了坡度和地震動加速度幅值對地震土坡滑坡體大小和致災(zāi)范圍的影響規(guī)律。本研究得到的主要研究結(jié)論可歸納如下：

(1)本文利用SPH動力分析方法對已有土坡破壞振動臺試驗(yàn)進(jìn)行了模擬，得到了與試驗(yàn)類似的破壞過程、基本吻合的坡形、及數(shù)值上非常接近的最大滑動距離，驗(yàn)證了SPH動力分析方法在地震土坡變形分析中的精度和準(zhǔn)確性。

(2)通過對比同一坡度土坡在不同地震動加速度幅值作用下的SPH模擬結(jié)果，可知同一土坡在不同地震動幅值作用下存在相似的滑動面，且潛在滑動區(qū)域隨著地震動幅值的增大而增加。

(3)通過對比不同坡度土坡在地震動作用下相應(yīng)位移閾值確定的土潛在滑動面與最大剪應(yīng)變確定的潛在滑動面，證實(shí)了利用位移閾值確定土坡滑動面的方法具有一定的準(zhǔn)確性。

(4)通過對比不同坡度土坡在同一地震動加速度幅值作用下的SPH模擬結(jié)果，可知隨坡度的增加，各土坡在地震作用下的穩(wěn)定性下降，致災(zāi)范圍變大；但滑坡體的區(qū)域面積會減小。