基于連續(xù)-非連續(xù)耦合方法的降雨滑坡數(shù)值模擬研究

陳聞瀟， 石 崇， 李汪洋， 張一平

（1.河海大學(xué)巖石力學(xué)與堤壩工程教育部重點實驗室，南京 210098； 2.河海大學(xué)巖土工程科學(xué)研究所，南京 210098）

滑坡過程分析是巖土工程評價與災(zāi)害評估的重要內(nèi)容［1］，工程師通過分析滑坡過程中出現(xiàn)的裂隙位置、分布、孕育發(fā)展過程判斷滑坡破壞模式，進(jìn)一步評估潛在滑坡災(zāi)害的規(guī)模和范圍，探討合理的防災(zāi)減災(zāi)措施［2-5］.

滑坡數(shù)值仿真是進(jìn)行滑坡過程分析的重要手段. 對于降雨滑坡的模擬分析，李寧等［6］提出了改進(jìn)的Mein-Larson 降雨入滲模型；蔣中明等［7］建立了基于FLAC有限差分的非飽和滲流計算方法；謝強等［8］分析了降雨影響下邊坡的安全系數(shù)變化. 但是這些學(xué)者提出的降雨計算模型都是基于連續(xù)介質(zhì)力學(xué)，不能分析降雨滑坡造成的大變形破壞、滑坡運動過程、堆積形態(tài)與災(zāi)害影響范圍等.

離散元方法被廣泛運用到計算邊坡的滑坡過程研究中［9-10］，與連續(xù)力學(xué)方法相比，其具有能模擬大變形的優(yōu)勢. 在離散元方法中，考慮降雨影響的主要方法是參數(shù)折減或計算等效入滲深度［11］，或采用流固耦合［12］的方法，但是這些方法中假定土體飽和，不滿足土的飽和-非飽和特性. 離散元中采用的簡化方法勢必會造成結(jié)果的偏差. 另外，滑坡往往僅限于一定范圍，如果整個模型都采用離散單元方法，很容易造成顆粒或塊體數(shù)目過多，導(dǎo)致計算工作量超出電腦承受能力.

因此部分學(xué)者建議采用連續(xù)-非連續(xù)耦合計算方法分析滑坡等大變形問題［13-14］，該方法將可能發(fā)生破壞的土體用離散元求解，而變形量小的周邊巖土體用連續(xù)介質(zhì)方法計算，能極大地提高計算速度與精度［15］. 由于離散元進(jìn)行飽和-非飽和降雨入滲模擬相對困難，然而在連續(xù)數(shù)值模擬方法中很容易做到，如果將兩者結(jié)合，進(jìn)行連續(xù)-非連續(xù)耦合計算，那么不僅能模擬降雨入滲，同時可以模擬降雨誘導(dǎo)的滑坡大變形破壞.

本文基于連續(xù)-非連續(xù)耦合分析方法，利用邊界墻耦合傳遞方法，提出了降雨條件下滑坡的模擬分析方法，借助案例分析了該方法的可行性，并探討了耦合滑坡數(shù)值仿真過程中的控制因素及影響規(guī)律，研究結(jié)果可為連續(xù)-非連續(xù)數(shù)值仿真在巖土工程中的應(yīng)用提供依據(jù).

1 連續(xù)-非連續(xù)耦合分析實現(xiàn)方法

目前連續(xù)-非連續(xù)數(shù)值模擬實現(xiàn)方法主要有兩種：基于邊界控制顆粒［16］和基于邊界控制墻體［17］的耦合方法. 本文采用的連續(xù)-非連續(xù)耦合分析算法以有限差分算法（FLAC）和顆粒離散元方法（PFC）原理為基礎(chǔ)，通過基于邊界控制墻體的方法來進(jìn)行信息交互與耦合計算，屬于邊界控制墻體方法. 邊界墻由三角形面組成，其頂點速度和位置為時間的函數(shù). 耦合邏輯的工作原理是：在每一個耦合分析步中，獲取非連續(xù)模型與邊界墻接觸的力和力矩，通過等效力系統(tǒng)，將其施加于連續(xù)模型的節(jié)點上，在連續(xù)模型中計算后，將節(jié)點的速度和位置傳遞給非連續(xù)顆粒，通過非連續(xù)模型中力-位移法則，計算得出新時間步的力和力矩，如此就完成了一個耦合循環(huán).

在耦合模型中當(dāng)離散顆粒與三角形的邊界墻相接觸時，在墻面上的接觸點為CP. 如圖1所示，三角形邊界墻頂點位置為xi（i=1，2，3），將三角形的頂點與點CP連接，得到三個區(qū)域的面積為Ai（i=1，2，3），三角形的總面積為A. 定義每個頂點的權(quán)重因子wi（i=1，2，3）為一個頂點對邊的三角形面積除以三角形的總面積，即wi=Ai/A. 定義ri（i=1，2，3）為從點CP指向各自三角形頂點的向量，有ri=xi-CP.

非連續(xù)模型每一時步將力和力矩傳遞至連續(xù)模型，其實現(xiàn)需要通過等效力方法. 作用于接觸點CP處的接觸力為F，黏接產(chǎn)生的接觸力矩為M. 通過等效力方法計算后的作用力位于邊界墻的每個頂點xi（i=1，2，3），其力的大小為Fi（i=1，2，3），通過下式計算：

圖1 連續(xù)-非連續(xù)耦合邊界墻Fig.1 Continuous-discontinuous coupling boundary wall

定義n為指向三角形法線方向的單位向量，則沿三角形邊界墻表面的切向力矢量為：

切向單位矢量為：

力的求解通過建立局部坐標(biāo)系，使得局部坐標(biāo)系中x分量與n方向一致，y分量與s方向一致，此時每個ri的x分量為0，即ri,x=0 . 此外，將重心加權(quán)wi（i=1，2，3）施加在三角形平面上最大接觸力的方向. 這種簡化可以直接求解局部坐標(biāo)系統(tǒng)中頂點的力和力矩.

連續(xù)模型每一時步將節(jié)點的速度傳遞至非連續(xù)模型. 節(jié)點傳遞至邊界墻頂點xi（i=1，2，3）的速度分別為Vi（i=1，2，3），假定速度場在三角形中線性變化，則可以通過重心插值方式得到接觸點CP處的速度值，并將該點的速度值作為非連續(xù)顆粒的速度值，其計算公式如下：

2 降雨荷載條件考慮與計算

2.1 降雨入滲的飽和-非飽和計算

本文提出的降雨滑坡的連續(xù)-非連續(xù)耦合分析流程如圖2所示，模型的模擬分為降雨過程的模擬和滑坡過程的模擬. 降雨模擬通過連續(xù)力學(xué)方法計算，由于降雨是長時間尺度的過程，其持續(xù)時間通常以小時和天為單位，在分析時通常不考慮變形［7，18］，故固定變形以節(jié)省效率，采用基于飽和-非飽和滲流理論對降雨進(jìn)行計算，并對孔壓場和滲流場進(jìn)行記錄. 滑坡過程的模擬將允許模型自由變形，將可能發(fā)生破壞的區(qū)域進(jìn)行離散化，基于耦合邊界墻的方法生成連續(xù)-非連續(xù)模型，并將孔壓場和滲流場通過等效力和土體弱化的方法疊加進(jìn)模型，進(jìn)而對降雨滑坡進(jìn)行連續(xù)-非連續(xù)耦合計算分析.

圖2 耦合方法流程圖Fig.2 Flowchart of coupling method

降雨過程模擬采用飽和-非飽和滲流理論. 非飽和土滲流分析中考慮土-水特征線和水力傳導(dǎo)方程，其中土-水特征線用來描述含水率與吸力的關(guān)系，水力傳導(dǎo)方程用來描述滲透系數(shù)與吸力的關(guān)系. 本文采用常用的Van Genuchten（VG）模型［19］來模擬非飽和土的水力特性. VG模型中土-水特征曲線的表達(dá)式為：

式中：θ為土的體積含水率；θr為殘余體積含水率；θs為飽和體積含水率；P 為土體的孔隙水壓力；ɑ、n′、m′為擬合參數(shù). 由于體積含水率與飽和度s 關(guān)系為θ=n·s，通過上式可得到土體飽和度與負(fù)孔隙水壓力的關(guān)系式為：

式中：Sr為殘余飽和度，非飽和土的滲透系數(shù)k( s) 與飽和滲透系數(shù)K的關(guān)系為：

為了實現(xiàn)降雨邊坡的模擬，需要對邊坡表面的降雨邊界條件進(jìn)行實時調(diào)整修正，即動態(tài)邊界條件. 當(dāng)降雨強度小于土體的入滲率時，邊坡表面不會產(chǎn)生積水，入滲量的值取降雨強度值. 當(dāng)降雨強度大于土體的入滲率時，邊坡表面會產(chǎn)生積水，入滲量的值取土體的最大入滲率. 在計算的每一個時間步都會通過自定義FISH函數(shù)，將降雨強度與土體的入滲率的大小進(jìn)行對比，從而實現(xiàn)動態(tài)邊界條件的模擬.

2.2 降雨影響的連續(xù)-非連續(xù)等效力方法

在通過飽和-非飽和滲流計算，得到邊坡的滲流場與孔壓場之后，需要通過等效簡化的方式，將降雨影響施加到非連續(xù)力學(xué)的顆粒中. 首先必須考慮由于雨水入滲而引起的土體重度的增加，受降雨影響的土體重度ρ 可按下式計算：

式中：ρd為土的干重度；θ為體積含水率；ρw為水的重度. 在飽和土體中由于水頭壓力差的存在會對土體產(chǎn)生滲流力，滲流力的方向與滲流方向一致［20］，且忽略非飽和土體內(nèi)滲流力的影響. 在飽和區(qū)域內(nèi)施加的等效滲透作用力Ff為：

式中：J為滲透力；V為單元網(wǎng)格或顆粒單元的體積；γW為水的重度（N/m3），i為水力梯度，其大小值為土體中兩點水勢之差與其滲透距離的比值. 顆粒在飽和作用下受到的浮力的大小可以由以下公式計算得出：

此外，土體在飽和狀態(tài)下，由于軟化作用使其物理力學(xué)參數(shù)存在一定的折減，本文根據(jù)工程經(jīng)驗將飽和區(qū)內(nèi)的土體參數(shù)折減15%. 在非連續(xù)模型中，降雨荷載的影響考慮由公式（8）～（10）計算的土重增加、滲流力和浮力及參數(shù)折減. 在連續(xù)模型中，由于不是滑坡發(fā)生的主要區(qū)域，僅考慮公式（8）計算的土重增加和參數(shù)折減，以保證合理的計算效率.

3 降雨滑坡的耦合計算分析

本文通過一個案例來分析該連續(xù)-非連續(xù)耦合降雨滑坡模擬方法的可行性，為避免關(guān)注點在于誘發(fā)滑坡的降雨閾值問題，本文選取的邊坡安全系數(shù)為1.1，在給定降雨影響后條分法安全系數(shù)小于1.0，即邊坡將在降雨影響下產(chǎn)生滑動. 如圖3所示建立連續(xù)邊坡模型，在此基礎(chǔ)上開展降雨條件下連續(xù)-非連續(xù)耦合滑坡分析.為了兼顧問題的三維性，使用假三維模型來進(jìn)行平面問題的模擬. 模型共有14 636個網(wǎng)格單元和29 796個網(wǎng)格節(jié)點，假三維厚度為3 m，邊坡土體飽和滲透系數(shù)為1.0×10-4cm/s，公式（5）～（7）所需的VG 模型參數(shù)采用文獻(xiàn)［7］中給出值，以反映滲透系數(shù)、飽和度與基質(zhì)吸力之間的關(guān)系.土體內(nèi)非飽和區(qū)的滲透系數(shù)由式（5）～（7）計算.如圖3所示，初始孔壓場通過在模型兩側(cè)邊界設(shè)定孔壓邊界條件，然后采用飽和-非飽和計算穩(wěn)定生成，初始地下水位線（黑線）為孔隙水壓力為0的等值線，孔隙水壓力為負(fù)數(shù)的區(qū)域為非飽和區(qū).

圖3 初始計算模型及孔壓分布（單位：kPa）Fig.3 Initial calculation model and pore pressure distribution

給定降雨強度為1.0×10-6m/s，連續(xù)降雨10 d，通過上文描述的動態(tài)邊界條件方法，將降雨荷載施加于模型的上部邊界. 圖4 為邊坡在降雨10 d 后的孔隙水壓力分布及飽和度等值圖，可見此時邊坡的坡腳，坡面以及坡頂已經(jīng)由負(fù)孔隙水壓力升至正孔壓，這表明此區(qū)域內(nèi)的土體已經(jīng)由非飽和轉(zhuǎn)化為暫態(tài)飽和區(qū)，且在坡腳處與坡體下部的正孔壓區(qū)連通，暫態(tài)飽和區(qū)內(nèi)孔隙水壓力最大不超過50 kPa. 由于邊坡坡面及坡頂?shù)乃衷谥亓ψ饔孟孪蚱履_處下滲，所以坡腳處的飽和區(qū)范圍更大，地下水位上升幅度也較為明顯.

圖4 降雨計算結(jié)果Fig.4 Rainfall calculation results

在降雨影響結(jié)果計算完成后，為了繼續(xù)對降雨導(dǎo)致滑坡過程進(jìn)行模擬，需通過邊界墻耦合方法，建立連續(xù)-非連續(xù)耦合模型. 將降雨計算的結(jié)果保存至單元內(nèi)另開辟的存儲節(jié)點內(nèi)，目的是在耦合模型生成后，通過上文的等效力方法計算降雨荷載并分別施加于連續(xù)模型和非連續(xù)模型中. 然后通過連續(xù)-非連續(xù)耦合方法，將滑坡區(qū)域的潛在危險部分用離散介質(zhì)替換. 本案例中認(rèn)為整個邊坡均屬于危險區(qū)域，在工程應(yīng)用中，可結(jié)合實際滑坡體及其鄰近區(qū)域設(shè)為危險區(qū)域［2］. 為了保證潛在滑坡體與連續(xù)介質(zhì)的連續(xù)性，防止平滑的薄弱交界面出現(xiàn)，連續(xù)-非連續(xù)的接觸面采用鋸齒狀分布.

耦合模型如圖5 所示，共生成18 153 個半徑處于0.15～0.30 m 之間的顆粒. 通過宏觀-細(xì)觀力學(xué)參數(shù)標(biāo)定，得到模型宏-細(xì)觀力學(xué)參數(shù)如表1、表2所示，局部阻尼值為0.3. 由圖5（c）、（d）、（e）可見離散顆粒和連續(xù)單元之間通過邊界墻來形成連接，兩者緊密接觸，這是典型的基于邊界控制墻的離散-連續(xù)耦合方法.圖5（a）為連續(xù)非連續(xù)初始模型的總應(yīng)力圖，可見在將連續(xù)介質(zhì)替換成離散介質(zhì)后，傳遞至下方的力滿足自然應(yīng)力分布，模型的總應(yīng)力符合條件，但由于離散顆粒具有一定的隨機(jī)性，所以在容許情況下應(yīng)力存在細(xì)微的波動. 圖5（b）為模型初始平衡時產(chǎn)生的位移云圖，可見連續(xù)-非連續(xù)模型之間的位移是連續(xù)的，這也證明了耦合的可行性. 但是與單一介質(zhì)模型平衡時產(chǎn)生的位移不同，耦合時因為初始狀態(tài)生成時，離散顆粒與單元網(wǎng)格之間的接觸并不緊密，接觸力也并不均勻，所以最大位移產(chǎn)生在豎直耦合交界面的附近，直到模型逐漸平衡，傳遞均勻的接觸力，位移才不再發(fā)展，形成初始模型，因此在之后的模型計算中，需先對此位移清零.

表1 計算模型細(xì)觀力學(xué)參數(shù)Tab.1 The mesomechanical parameters of calculation model

表2 計算模型宏觀力學(xué)參數(shù)Tab.2 The macromechanical parameters of calculation model

圖5 連續(xù)-非連續(xù)耦合模型Fig.5 Continuous-discontinuous coupling model

圖6 滑坡過程Fig.6 Landslide process

降雨滑坡耦合模型計算結(jié)果如圖6所示. 由圖6（a）、（b）可見，在降雨入滲形成飽和區(qū)的影響下，邊坡滑動且形成了明顯的圓弧滑動面，滑坡類型主要為淺層滑坡，滑動區(qū)域以飽和區(qū)為主，圓弧最大深度處超過了暫態(tài)飽和區(qū)分界線，這說明了非飽和區(qū)的土發(fā)生了滑坡. 在降雨影響下，滑坡體迅速沿著坡面滑落，并由于局部阻尼的影響，滑坡體最終在坡腳堆積并停止進(jìn)一步運動. 滑坡的最大位移為29.3 m，滑坡結(jié)束后堆積角約為31.6°. 由圖6（c）、（d）可見在滑坡啟動階段，各測點的速度迅速增大，而測點3的速度在7 s左右才明顯增大，這表明了該滑坡為牽引式滑坡，前緣由于在降雨作用下更容易飽和，飽和區(qū)范圍較大，土體強度弱化而率先滑動，而處于后緣的土體破壞較遲，其破壞原因受降雨影響和前緣破壞失穩(wěn)的雙重影響. 滑坡體速度在10～25 s左右達(dá)到峰值，隨后速度逐漸減小，在175 s左右滑坡結(jié)束. 圖6（e）為滑坡過程中的裂隙數(shù)變化曲線，由于降雨作用的影響，在滑坡的初期就產(chǎn)生了一定的裂隙，并在滑坡過程中，出現(xiàn)了幾個裂隙陡增的時段，這說明可能發(fā)生了小的崩塌. 計算結(jié)果表明，采用連續(xù)-非連續(xù)耦合計算方式，對降雨-滑坡進(jìn)行耦合分析，能夠較好地進(jìn)行降雨作用下滑坡過程和破壞特征分析.

4 降雨滑坡中的阻尼影響

在滑坡分析中，通常通過設(shè)置阻尼來控制能量衰減，防止動能累積過快［21］. 故在降雨滑坡耦合分析過程中阻尼的參數(shù)尤為重要，在不同阻尼條件下的滑坡過程分析圖7所示. 由圖7（a）可見，在滑坡進(jìn)行的0～30 s內(nèi)，飽和區(qū)內(nèi)土體的速度和阻尼呈負(fù)相關(guān)的趨勢，阻尼越小，平均速度峰值越大，達(dá)到峰值所需的時間也越短. 在滑坡發(fā)生30 s之后，平均速度在不同阻尼影響下呈現(xiàn)相同的趨勢，即隨著滑坡進(jìn)行不斷減小，滑坡體逐漸穩(wěn)定. 總體上，在不同阻尼下，滑坡趨勢均保持一致，即速度先迅速增大至峰值，然后逐漸減小至穩(wěn)定，位移表現(xiàn)為先迅速增大而后逐步穩(wěn)定. 由圖7（b）可見，暫態(tài)飽和區(qū)的平均位移受局部阻尼的影響，具體體現(xiàn)在距坡腳的堆積距離上，阻尼越大堆積距離越小，平均位移也就越小.

從圖7（c）、（d）可以看出，當(dāng)阻尼小于0.3時，滑坡體積與堆積距離明顯增加，這是因為阻尼較小，能量耗散緩慢，在滑坡中產(chǎn)生的沖擊碰撞以及堆積擠壓等具有更大的動能，將會使得滑坡造成更進(jìn)一步的破壞，故滑坡影響區(qū)域也越大. 而當(dāng)阻尼大于0.3時，滑動區(qū)域變化不明顯，堆積距離也基本一致. 合理的獲取阻尼參數(shù)，在降雨滑坡耦合分析過程中尤為重要，通過將不同阻尼的計算結(jié)果，與真實滑坡中的滑坡征兆、運動過程和堆積狀態(tài)對比，來確定合理的阻尼參數(shù)，進(jìn)而實現(xiàn)真實情況下的降雨滑坡耦合模擬，并進(jìn)行更進(jìn)一步的分析.

圖7 局部阻尼影響因素分析Fig.7 Analysis of influencing factors of local damping

5 結(jié)論

本文基于連續(xù)-非連續(xù)數(shù)值模擬耦合分析方法，采用邊界墻耦合理論和飽和-非飽和降雨入滲分析，通過案例對降雨滑坡過程進(jìn)行了模擬，分析了該耦合計算方法的可行性，并研究了局部阻尼在降雨滑坡耦合分析中的影響規(guī)律，得到主要結(jié)論如下：

1）利用連續(xù)模型對飽和-非飽和降雨入滲分析，通過邊界墻耦合方法生成連續(xù)-非連續(xù)耦合模型，并將滲流場通過等效力施加于耦合模型中，實現(xiàn)了降雨影響下的滑坡耦合分析，從而建立了基于連續(xù)-非連續(xù)耦合分析的降雨滑坡分析方法.

2）通過案例分析驗證了降雨滑坡耦合分析的可行性. 計算結(jié)果表明該方法可以較好地模擬降雨滑坡，并進(jìn)行滑坡過程和破壞特征的分析，為降雨滑坡的分析提供了新思路.

3）分析了降雨滑坡耦合分析中，局部阻尼影響因素對計算結(jié)果的影響. 結(jié)果表明局部阻尼對滑坡結(jié)果的影響較大，通過實際滑坡的滑坡征兆、運動過程和堆積狀態(tài)，合理地獲取阻尼參數(shù)，在降雨滑坡耦合分析過程中尤為重要.