趙 洲,魏江波
(西安科技大學(xué) 地質(zhì)與環(huán)境學(xué)院,陜西 西安 710054)
基于數(shù)值分析的滑坡破壞過程模擬、機(jī)理分析與強(qiáng)度評價是滑坡災(zāi)害研究的熱點之一。目前,國內(nèi)外對滑坡破壞機(jī)理等方面的數(shù)值模擬研究多采用有限單元法(FEM)和有限差分法(FDM),這2種方法均假設(shè)研究體為連續(xù)體,忽略了巖土體顆粒間的復(fù)雜相互作用及實際巖土體的高度非線性行為,對大變形大位移求解尤其土質(zhì)滑坡分析也不夠理想,對整個滑坡的穩(wěn)定情況評價力度不足。顆粒流方法(particle flow code,PFC)是一種特殊的離散單元法[1-2],通過剛性的圓形或者球體來模擬巖土體顆粒,以力-位移定律和牛頓第二定律為基本理論模擬顆粒的運(yùn)動和相互作用,運(yùn)用特定的粘結(jié)模型將其集合于一體來模擬巖土體的整體性,用滑動模型模擬顆粒之間的剪切運(yùn)動及巖土體宏觀大變形,通過細(xì)觀顆粒的運(yùn)動特征來反映宏觀巖土體的變形破壞特征。同時,顆粒流方法最大優(yōu)點在于不用假設(shè)滑動面位置以及可直觀的觀察巖土體的動態(tài)真實破壞過程[3]。
近年來運(yùn)用顆粒流方法對滑坡以及邊坡的研究也越來越多。周健等、廖彪應(yīng)用顆粒流方法分析了土體細(xì)觀參數(shù)對土質(zhì)邊坡破壞形式的影響,認(rèn)為隨著土體顆粒黏性的增大,其破壞類型從塑性破壞向脆性破壞發(fā)生過渡,并得出了顆粒流方法與傳統(tǒng)方法對均質(zhì)土坡穩(wěn)定性計算結(jié)果比較接近的結(jié)論,為邊坡穩(wěn)定性分析開辟了一條新的途徑[3-5]。張龍以雞尾山巖質(zhì)滑坡為例,對滑坡的運(yùn)動過程、破壞方式及堆積規(guī)律進(jìn)行了顆粒流分析,認(rèn)為摩擦系數(shù)及土體強(qiáng)度對滑體物質(zhì)在堆積區(qū)的分布具有重要影響,但對滑體的最大位移影響不大[6]。馬建全對黑方臺焦家崖頭灌溉影響下的黃土滑坡進(jìn)行了基于PFC2D的顆粒流模擬,展現(xiàn)了滑坡的整個運(yùn)動過程,并分析了滑坡各階段的破壞特征[7]。廖靜薇基于試驗分析與顆粒流理論,建立了粉質(zhì)粘土宏觀力學(xué)參數(shù)與細(xì)觀力學(xué)參數(shù)的對應(yīng)關(guān)系,提出了細(xì)觀力學(xué)下粉質(zhì)粘土邊坡的失穩(wěn)判據(jù)和邊坡安全系數(shù)計算方法,并結(jié)合實例應(yīng)用證明了顆粒流方法在土坡穩(wěn)定性分析中的有效性[8]。N.Li等利用PFC進(jìn)行了荷載作用下土質(zhì)邊坡的變形破壞機(jī)制研究,并與室內(nèi)物理模擬試驗進(jìn)行了對比,取得的模擬結(jié)果十分相近[9]。Valentino等通過顆粒水槽試驗?zāi)M了顆粒流的運(yùn)動過程,并采用離散元軟件PFC2D進(jìn)行數(shù)值模擬水槽試驗,同時獲取顆粒流沖出距離和沖擊力頻譜特性[10]。唐紅梅等借助PFC2D軟件對三峽庫區(qū)龔家方2號斜坡破壞過程進(jìn)行了數(shù)值模擬,得出模擬結(jié)果與實際破壞情況相一致的結(jié)論[11]。吳越等基于功能關(guān)系,通過物理實驗與數(shù)值模擬結(jié)合對比,從滑坡力學(xué)特征的能量角度研究了滑坡運(yùn)動過程中的能量耗散以及滑體對承載體的沖擊作用,建立了承災(zāi)體損毀程度與滑坡沖擊作用的定量函數(shù)關(guān)系[12-15]。
由此可見,目前采用顆粒流方法對滑坡的研究多集中在滑坡顆粒流模型的構(gòu)建、巖土體細(xì)觀參數(shù)對滑坡破壞方式的影響、滑坡運(yùn)動過程的模擬、滑坡穩(wěn)定性計算的對比分析等方面,而對滑坡破壞過程中的力學(xué)機(jī)理分析、滑坡速度、動能及沖擊力等反應(yīng)滑坡強(qiáng)度變化特征指標(biāo)的分析和研究相對較少。
利用顆粒流軟件PFC2D,以府谷縣新府山滑坡為例,建立了滑坡顆粒流模型,模擬了滑坡從變形失穩(wěn)到運(yùn)動停止的動態(tài)破壞全過程,分析了滑坡變形破壞特征及其力學(xué)機(jī)理,研究了滑坡強(qiáng)度的變化特征,給出了以沖擊力為表征的滑坡沖擊破壞強(qiáng)度定量分析方法。
新府山滑坡位于陜西省府谷縣府谷鎮(zhèn)新府山安居小區(qū)北側(cè)斜坡(圖1)。該滑坡體前后緣相對高差25 m,寬約148 m,長約40 m,厚約3~3.5 m,體積約1.1×104m3,坡體中部較緩,坡頂、坡腳部較陡,滑坡坡面總體近直線型,傾向東南?;轮芙绶置?,左右兩側(cè)以沖溝為界?;麦w為第四系全新統(tǒng)褐色風(fēng)積粉土,下伏上古生界二疊系棕紅色泥巖,產(chǎn)狀200°∠35°。由于坡體前緣開挖使得坡體下部變陡,加上降雨下滲作用使得坡體強(qiáng)度降低,降雨匯集于泥巖頂面,軟化泥巖結(jié)構(gòu),降低抗剪強(qiáng)度,從而誘發(fā)坡體沿粉土-泥巖接觸面下滑,屬典型的土—巖接觸性滑坡(圖2)。
圖1 新府山滑坡平面位置及發(fā)育特征示意圖Fig.1 Location and characteristics of Xinfushan landslide
圖2 新府山滑坡地質(zhì)剖面Fig.2 Geological section of Xinfushan landslide
巖土體力學(xué)參數(shù)是滑坡穩(wěn)定性評價的關(guān)鍵和難點[16]。文中根據(jù)土樣試驗及工程經(jīng)驗分析,確定的新府山滑坡相關(guān)巖土體宏觀力學(xué)參數(shù)見表1.
PFC賦予巖土體的是其顆粒的細(xì)觀參數(shù),包括顆粒密度、粒徑、法向剛度、切向剛度、摩擦系數(shù)等,并通過顆粒的細(xì)觀參數(shù)來反映巖土體的宏觀特征。
表1 新府山滑坡巖土體力學(xué)參數(shù)Table 1 Rock and soil mechanical parameters of Xinfushan landslide
為了獲得宏觀參數(shù)相對應(yīng)的細(xì)觀參數(shù),利用PFC2D程序中的雙軸試驗?zāi)P头椒?gòu)建了6 m×12 m的矩形模型(圖3),上下墻(wall)施加壓力,左右墻模擬約束。在構(gòu)建好的矩形墻內(nèi),用半徑擴(kuò)張法結(jié)合generate命令生成一定數(shù)量一定空隙率大小的顆粒[17-20]。為了更好的模擬巖土體的真實特征,顆粒(ball)半徑越小越好,但受計算機(jī)的容量與運(yùn)行速度的限制以及所需精度要求,可適當(dāng)擴(kuò)大顆粒半徑。最后根據(jù)土體顆粒特性,選用顆粒的接觸粘結(jié)模型,并用property命令生成接觸粘結(jié)模型。
在雙軸試驗?zāi)P徒⒑弥?,賦予顆粒一組細(xì)觀參數(shù)及不同圍壓,進(jìn)行模擬并獲得對應(yīng)偏應(yīng)力-應(yīng)變曲線,進(jìn)而獲得對應(yīng)的摩爾應(yīng)力圓。進(jìn)行大量的試算模擬,通過每一組細(xì)觀參數(shù)與不同的圍壓所獲得的摩爾應(yīng)力圓繪制出對應(yīng)的強(qiáng)度包絡(luò)線并反算出宏觀參數(shù)與表1參數(shù)進(jìn)行對比,獲取與表1相匹配的那組細(xì)觀參數(shù)作為實際巖土體的模擬參數(shù)。借鑒文獻(xiàn)12中的宏細(xì)觀參數(shù)的定量關(guān)系式[21],并通過多次試算獲得了新府山滑坡巖土體顆粒的細(xì)觀參數(shù)(表2)。
圖3 雙軸試驗?zāi)P虵ig.3 Biaxial test model
表2 新府山滑坡巖土體顆粒細(xì)觀參數(shù)Table 2 Rock and soil microscopic parameters of Xinfushan landslide
以粉土為例,賦予表2的細(xì)觀參數(shù),對粉土進(jìn)行圍壓分別為50,100和200 kPa的雙軸試驗?zāi)M,并得知50 kPa圍壓下峰值強(qiáng)度為190,100 kPa圍壓下峰值強(qiáng)度為310,200 kPa圍壓下峰值強(qiáng)度為570 kPa,則可繪制Mohr應(yīng)力圓與Mohr-Coulumn破壞包絡(luò)線,計算得粘聚力約為21.3 kPa,內(nèi)摩擦角約為26.8°(圖4)。對比表1可知,表2中的細(xì)觀參數(shù)可用于后文滑坡模型模擬。
圖4 Mohr應(yīng)力圓與Mohr-Coulumn包絡(luò)線(粉土)Fig.4 Mohr stress circle and Mohr-Coulumn line(silt)
為了更真實的模擬并反映滑坡的特性,解決傳統(tǒng)半徑擴(kuò)張法存在的巖土體層間無接觸問題與放大系數(shù)難確定問題,文中結(jié)合“分層下落法”的思路構(gòu)建滑坡模型,該方法不需要建立坡面約束墻和層間分界墻,不用消除浮球,且存在層間接觸[22],與實際情況更加吻合。根據(jù)圖1尺寸與表2參數(shù)構(gòu)建滑坡模型,過程如圖5所示。
在建模過程中,以單元平均不平衡力達(dá)單元平均接觸力的1%,且最大不平衡力達(dá)最大接觸力的1.0×10-5時作為平衡條件,即循環(huán)停止,且平衡后進(jìn)行速度與位移清零。為便于區(qū)分上下層巖土體,將其進(jìn)行分組(group),并設(shè)置不同顏色。保留斜坡邊界無剛度墻便于觀察顆粒運(yùn)動過程中的情況,其只作為參考作用,并無阻礙作用。
初始模型構(gòu)建完成后,刪除右側(cè)豎直墻(wall),并在坡腳處建一水平剛性墻,設(shè)置其摩擦系數(shù)(friction)為1.0,用以模擬坡腳處水平地面。同時為了研究分析滑坡體應(yīng)力變化,在滑動面上中下不同位置設(shè)置測量圓,以監(jiān)測(history)記錄對應(yīng)位置在坡體滑動過程中的應(yīng)力變化狀況。最終模型如圖6所示。
圖5 分層下落法建模過程Fig.5 Modeling process of stratified rain fall
顆粒流強(qiáng)度折減法就是將顆粒的接觸粘結(jié)強(qiáng)度和摩擦系數(shù)同時折減一定的比例,使得滑坡剛好處于臨界位移狀態(tài),以臨界位移作為滑坡破壞判斷依據(jù),此時,將折減前強(qiáng)度參數(shù)(n_b,s_b和fric)與折減后的臨界強(qiáng)度參數(shù)(n_bcr,s_bcr和friccr)的比值定義為滑坡穩(wěn)定性系數(shù)。其公式如下
式中n_b為顆粒的法向接觸粘結(jié)強(qiáng)度;s_b為顆粒的切向接觸粘結(jié)強(qiáng)度;fric為顆粒間的摩擦系數(shù);n_bcr為臨界狀態(tài)時的顆粒法向接觸粘結(jié)強(qiáng)度;s_bcr為臨界狀態(tài)時的顆粒切向接觸粘結(jié)強(qiáng)度;friccr為臨界狀態(tài)時顆粒間摩擦系數(shù)。
圖6 新府山滑坡顆粒流模型Fig.6 PFC2D model of Xinfushan landslide
根據(jù)公式(1)進(jìn)行滑坡強(qiáng)度折減模擬。起初設(shè)定折減系數(shù)為1.0,逐漸增大折減系數(shù),并進(jìn)行模擬觀察顆粒位移場的變化,當(dāng)折減系數(shù)增大至1.115時滑坡剛好處于臨界位移狀態(tài),此時的折減系數(shù)即為滑坡的穩(wěn)定性系數(shù),與采用SLOP/W模塊求得的滑坡穩(wěn)定性系數(shù)為1.116的結(jié)果十分接近[23]??芍摶略谧匀还r下處于基本穩(wěn)定狀態(tài),但在人類工程活動、強(qiáng)降雨等誘發(fā)因素的影響下,不排除其復(fù)活的可能性,一旦滑坡強(qiáng)度達(dá)到臨界狀態(tài)時,該滑坡將繼續(xù)失穩(wěn)破壞。因此,有必要對滑坡的破壞過程及其機(jī)理做進(jìn)一步分析。
3.2.1 從變形特征分析滑坡破壞機(jī)理
選取上述強(qiáng)度折減法所求的臨界位移狀態(tài)時的細(xì)觀參數(shù),即將巖土體細(xì)觀強(qiáng)度參數(shù)折減1.115倍,采用PFC2D程序?qū)缕茐倪^程進(jìn)行仿真模擬。
根據(jù)顆粒流模擬結(jié)果分析,滑坡在初始時段處于變形狀態(tài),應(yīng)力重新分布,坡腳顆粒受上覆壓力作用下變形不斷積累,3.0×104步時坡腳擠壓破壞(圖7(a))。坡腳擠壓破壞后,坡體失去支撐,變形加劇,顆粒根據(jù)自身的接觸力自行調(diào)整位置,并且使得坡體中部抗剪強(qiáng)度最弱處產(chǎn)生剪切破壞。在上覆壓力作用下,變形繼續(xù)加劇,坡體中的剪切裂隙處產(chǎn)生應(yīng)力集中,當(dāng)裂隙尖端應(yīng)力大于其抗剪強(qiáng)度時,裂隙尖端繼續(xù)破壞,裂隙延伸,表現(xiàn)為坡體中剪切裂隙向上下附近的粘結(jié)最弱處延伸。變形過程中,破體向下滑移,坡頂失去支撐,5.0×104步時坡頂產(chǎn)生明顯拉張裂縫(圖7(b))。當(dāng)上部拉張裂縫與中部剪切裂隙相貫通,并與坡腳擠壓破壞處相連接后,坡體將失穩(wěn)破壞,產(chǎn)生整體破壞并滑移。當(dāng)滑坡整體破壞后,上部臺階處坡頂顆粒向下滑落,達(dá)到一定時步后,其臺階之上的原滑坡后壁所形成的小邊坡失去下部支撐,1.8×105步時發(fā)生次生滑坡(圖7(c))。
圖7 不同時步滑坡破壞形態(tài)Fig.7 Failure form of landslide at different time
在滑坡變形—失穩(wěn)—運(yùn)動過程中觀察顆粒間的粘結(jié)(即黑色短線)變化情況(圖8),可看出,在滑坡再次失穩(wěn)破壞過程中,滑體顆粒間的粘結(jié)基本破壞,坡體愈加破碎,并順勢下滑堆積于坡腳。對于以位移分析法為主的顆粒流法,在滑坡運(yùn)動過程中,通過位移矢量圖可以直觀看出顆粒自動搜索的滑動面位置(圖9)。
圖8 t=150e4時步接觸粘結(jié)圖Fig.8 t=150e4 time step contact bond
3.2.2 從應(yīng)力特征分析滑坡破壞機(jī)理
根據(jù)滑坡破壞過程中應(yīng)力變化曲線的分析可知(圖10),坡頂由于產(chǎn)生拉張裂縫和垮塌,所以1號測量圓內(nèi)y方向的應(yīng)力明顯減小,x方向略有減
圖9 t=150e4時步位移矢量Fig.9 t=150e4 time step displacement vector
小,最后當(dāng)滑坡破壞后,1號測量圓所在滑坡后壁處于穩(wěn)定狀態(tài),故該位置應(yīng)力最后均趨于穩(wěn)定(圖10(a))。對于坡中由于坡體持續(xù)性向下滑移,且破碎程度不斷增大,使得2號測量圓所處位置上覆巖土體不斷減少,因而其所受豎向壓力減小,所以2號測量圓內(nèi)y方向應(yīng)力明顯降低,由于上部坡體破壞對下部巖土體的擠推作用與持續(xù)下滑使得x方向應(yīng)力起初略有增大后又稍微減小(圖10(b))。對于坡低,由于坡體下滑擠推作用及堆積作用使得3號測量圓內(nèi)x和y應(yīng)力都顯著增加(圖10(c))。
圖10 應(yīng)力變化模擬曲線Fig.10 Stress change simulation curve
滑坡強(qiáng)度是表達(dá)滑坡所儲備的破壞能力及對承災(zāi)體破壞程度的一項重要指標(biāo),也有學(xué)者認(rèn)為滑坡強(qiáng)度指滑坡具有的能量大小[24],其研究對滑坡風(fēng)險評估具有重要意義?;聫?qiáng)度大小與滑坡規(guī)模、斜坡坡度、坡高、滑坡巖土體類型等滑坡發(fā)育特征均具有一定聯(lián)系,然而目前仍沒有一種確定的滑坡強(qiáng)度定量分析與評價的方法。一般認(rèn)為,滑坡體的位置越高、滑體體積越大、移動速度越快,則滑坡強(qiáng)度也就越高,對承災(zāi)體的破壞能力也就越大。
為了分析新府山滑坡在變形與破壞過程中強(qiáng)度變化特征,文中從速度、動能、沖擊力等指標(biāo)入手,間接分析新府山滑坡強(qiáng)度大小的變化特征。
由于顆粒流模型中難以監(jiān)測出滑坡整體的運(yùn)動速度,因此文中選取滑坡破壞后顆粒最大速度vmax進(jìn)行分析(圖11),以及進(jìn)行實時監(jiān)測滑坡動能的變化(圖12)。雖然顆粒最大速度不能代表滑坡的整體運(yùn)動速度,但就特定滑坡而言,仍然可以反映出滑坡能量大小的變化特征。
從圖11,圖12對比可以看出,滑坡速度特征與動能特征基本一致?;略?.0×105時步破壞時,雖然速度最大,但滑坡破壞體積較小,所以此時的動能并非最大?;缕茐暮螅捎诨P作用使得滑坡破壞體積增大,直至1.5×105時步時,滑坡的破壞體積和速度都處于較大狀態(tài),此時動能達(dá)到最大。隨后,由于滑坡在水平地面的緩慢運(yùn)動使得動能逐漸降低。這一變化過程也符合滑坡破壞的突發(fā)性和堆積的緩慢性特征。
因此,從速度和動能指標(biāo)來看,新府山滑坡的強(qiáng)度在其破壞的初始階段為最高,然后逐漸變小。
圖11 顆粒最大速度變化曲線Fig.11 Maximum velocity variation curve of particles
圖12 滑坡運(yùn)動過程動能變化曲線Fig.12 Landslide kinetic energy variation curve
滑坡強(qiáng)度研究的主要目的之一是評價滑坡對建筑物等承災(zāi)體的沖擊破壞能力大小,強(qiáng)度越高的滑坡,對承災(zāi)體的破壞能力也越大[25-26]。從致災(zāi)后果來看,無論滑坡動能、速度等指標(biāo)如何變化,若不考慮承災(zāi)體自身抗災(zāi)能力特征,就滑坡自身而言,最終反映建筑物等承災(zāi)體損毀程度的變量是滑坡沖擊力大小。因此,為了分析新府山滑坡強(qiáng)度即滑坡對承災(zāi)體的沖擊破壞能力大小變化特征,文中選擇滑坡沖擊力這一指標(biāo)來進(jìn)行滑坡強(qiáng)度大小分析。
從圖11,圖12的分析可以看出,隨著新府山滑坡運(yùn)動的停歇,其速度、動能等呈衰減趨勢,也就是說,在水平地面上距離滑坡坡腳不同位置(文中將其稱為坡腳距)處的承災(zāi)體,其遭受滑坡的沖擊力大小不同,破壞程度也不同。
通過在沿滑面傾向方向上不同坡腳距(用L表示)處構(gòu)建豎直剛性wall,并分別監(jiān)測滑坡破壞過程中坡腳剛性wall上的沖擊力隨時間的變化,以此來分析滑坡對不同坡腳距處建筑物的沖擊作用,從而間接反映出滑坡強(qiáng)度大小及其變化特征。通過上文模擬分析和計算,新府山滑坡的最大滑移距離約為16 m,因此,在坡腳距從0 m開始每間隔1.0 m分別設(shè)置滑坡沖擊力測試的剛性wall,共建立17個沖擊力測試模型(圖13),分別實時監(jiān)測剛性wall上的水平?jīng)_擊力大小(圖14)。
圖13 剛性wall設(shè)置示意圖Fig.13 Rigid wall setup illustration
圖14 滑坡沖擊力隨時間的變化曲線Fig.14 Variation curves of impact force with time
由圖14可知,當(dāng)坡腳距一定時,起初滑動階段,由于滑坡速度與破壞的滑體體積不斷增大使得沖擊力隨滑坡的運(yùn)動而增加,當(dāng)滑體減速堆積并停止運(yùn)動時,沖擊力逐漸趨于穩(wěn)定(此時的沖擊力表現(xiàn)為滑體對剛性墻的靜止土壓力)。同時,當(dāng)坡腳距不同時,沖擊力隨坡腳距的增大明顯減小。然而由于滑坡體特征及運(yùn)動特征等因素影響,導(dǎo)致滑坡沖擊力峰值基本上接近于滑坡體失穩(wěn)運(yùn)動后的穩(wěn)定狀態(tài)下堆積體對擋墻的靜止土壓力。
對不同坡腳距上的最大沖擊力進(jìn)行統(tǒng)計分析可知(圖15),滑坡最大沖擊力Fmax與坡腳距Li之間的關(guān)系可以表達(dá)為
Fmax=761.46e-0.17/Li,(0≤Li≤16)
式中Fmax為滑坡最大水平?jīng)_擊力;Li為坡腳距。
圖15 滑坡沖擊力隨坡腳距的變化曲線Fig.15 Variation curves of impact force with Li
因此,以沖擊力表達(dá)的新府山滑坡強(qiáng)度大小隨著滑坡在水平地面運(yùn)動距離的增大而減小,并在運(yùn)動停止前呈現(xiàn)出指數(shù)衰減的特征。從滑坡破壞及運(yùn)動過程分析,當(dāng)滑坡體在整體發(fā)生破壞并順著滑面運(yùn)移到斜坡坡腳處時,由于具有較大的體積和速度,因而在坡腳處具有較大的沖擊力即破壞強(qiáng)度,但隨著滑坡體與水平地面的碰撞、摩擦運(yùn)動等,滑體內(nèi)部能耗、與地面之間的摩擦能耗增加,導(dǎo)致滑體顆粒運(yùn)動速度變緩,沖擊力減小,滑坡強(qiáng)度也隨之降低,直到滑坡停止運(yùn)動強(qiáng)度為零。
1)通過利用顆粒流模型對新府山滑坡失穩(wěn)破壞過程的仿真模擬及破壞機(jī)理分析,得出該滑坡變形破壞的力學(xué)機(jī)制表現(xiàn)為坡腳擠壓、坡中剪切、坡頂拉張;
2)通過對滑坡體顆粒運(yùn)動速度和動能的變化分析,認(rèn)為新府山滑坡的強(qiáng)度在其破壞的初始階段為最高,然后逐漸變??;
3)通過對不同坡腳距處滑坡沖擊力強(qiáng)度指標(biāo)的監(jiān)測分析,得出該滑坡對承災(zāi)體的沖擊破壞能力隨著坡腳距的增大呈指數(shù)衰減特征。