單治鋼 廖哲賢 董友扣 王 棟 崔 嵐

(①中國電建集團(tuán)華東勘測設(shè)計研究院有限公司， 杭州 311122， 中國) (②海洋巖土工程勘察技術(shù)與裝備浙江省工程研究中心， 杭州 311122， 中國) (③中國地質(zhì)大學(xué)(武漢)， 武漢 430074， 中國) (④中國海洋大學(xué)， 青島 266100， 中國) (⑤中國科學(xué)院武漢巖土力學(xué)研究所， 武漢 430071， 中國)

0 引 言

海底滑坡作為陸相沉積物向海相運移的主要方式之一，可以搬運大量由砂、黏土和碎石組成的沉積物(Hampton et al.，1996； 朱超祁等， 2016； 霍沿東等， 2019)。目前已知的最大規(guī)模海底滑坡Storegga滑坡體積超過5600 km3(Kvalstad et al.，2005)。與陸地滑坡相比，海底滑坡的滑動速度更快，可以高達(dá)35 m · s-1?；麦w材料呈非牛頓流體特性，高速運動的滑坡體可以對下游構(gòu)筑物(管線、纜線、防沉板和立柱等)造成嚴(yán)重破壞，并可能引發(fā)海嘯等次生災(zāi)害(Coyne et al.，2005； Nodine et al.，2006)。解秋紅等(2016)、修宗祥等(2016)和Dong et al.(2017a)采用數(shù)值方法反分析了海底滑坡的滑動過程。部分學(xué)者研究了海底滑坡對海底構(gòu)筑物的沖擊(王立忠等， 2008；Zakeri， 2009； Liu et al.，2015； Dong et al.，2017b， 2020a； 馮斌等， 2019； 王忠濤， 2019； Fan et al.，2021)，并基于此對構(gòu)筑物的外形進(jìn)行了優(yōu)化設(shè)計(范寧等， 2018)。

海底滑坡的滑動距離遠(yuǎn)遠(yuǎn)大于陸上滑坡：比如 Storegga滑坡最遠(yuǎn)滑動距離達(dá)到810 km(Kvalstad et al.，2005)，而陸上滑坡的最大滑動距離僅為數(shù)十公里(朱超祁等， 2016)。滑水被認(rèn)為是海底滑坡高速和長距離滑動的主要原因。滑水是指環(huán)境水侵入滑坡體下方形成水墊層，使滑坡體下的基床阻力遠(yuǎn)低于滑坡體的抗剪強(qiáng)度，滑坡體得以更大的速度向前滑動。Ilstad et al. (2004a)對大量的挪威海底滑坡案例進(jìn)行反分析發(fā)現(xiàn)：滑水產(chǎn)生的水墊層阻力比滑坡體抗剪強(qiáng)度低一個數(shù)量級?；畷鸹麦w前端脫離主體，形成一系列獨立的滑塊。此類滑塊整體性好、底部阻力小，因而滑動距離最遠(yuǎn)可達(dá)到數(shù)百公里。挪威FinneidFjord滑坡區(qū)域內(nèi)發(fā)現(xiàn)了大量的獨立滑塊，證明了真實海底滑坡中滑水的存在(Ilstad et al.，2004a)。De Blasio et al. (2004)通過室內(nèi)試驗發(fā)現(xiàn)滑水發(fā)生后滑坡體的最大滑動速度可達(dá)到35 m · s-1，并推算出滑動距離可達(dá)到125 km。Elverh?i et al. (2005)采用數(shù)值模擬方法對Bear Land和Grand Banks兩處滑坡分別進(jìn)行反分析，發(fā)現(xiàn)滑水導(dǎo)致滑坡體的滑動距離增加約30%。

目前對海底滑坡有效的研究手段還比較有限。由于海底環(huán)境復(fù)雜，且海底滑坡具有不可預(yù)測性和瞬時性，尚無法實時追蹤海底滑坡和滑水現(xiàn)象的全過程?，F(xiàn)有研究主要通過物探技術(shù)獲得海底地層剖面后對海底滑坡案例進(jìn)行事件后分析(來向華等， 2000； 王磊等， 2016)，有相當(dāng)?shù)牟淮_定性。室內(nèi)試驗可以借助高清攝像頭和粒子圖像測速(PIV)技術(shù)捕捉滑水的觸發(fā)過程(De Blasio et al.，2004； Ilstad et al.，2004a)，但是嚴(yán)重的尺度效應(yīng)問題限制了其應(yīng)用范圍。相比之下，數(shù)值分析為研究海底滑坡的滑水提供了更好的選擇，已在現(xiàn)有研究中得到了一定的應(yīng)用(Mohrig et al.，1998； 王立忠等， 2008)。

現(xiàn)有考慮滑水的海底滑坡模擬主要采用深度平均法(Imran et al.，2001)，將滑坡體簡化為插塞層(Plug layer)和剪切層(Shear layer)，對插塞層內(nèi)滑坡體的滑動速度沿深度進(jìn)行平均，可將三維或二維滑坡問題簡化為二維或一維問題。因而該方法計算效率較高，特別適用于大規(guī)模海底滑坡問題。然而，該方法僅能考慮滑坡體的連續(xù)塑性變形，不能模擬滑水造成的滑坡體前端脫離主體并形成獨立滑塊的過程(Dong et al.，2020b)，因而不適用于模擬滑水作用下海底滑坡的滑動全過程。

本文采用一種新型數(shù)值方法物質(zhì)點法(MPM)模擬海底滑坡的滑動過程，考慮滑水對基床阻力的影響，分別對室內(nèi)試驗和真實的海底滑坡案例進(jìn)行反分析，評估滑水對海底滑坡的危害范圍和滑動形態(tài)的影響。

1 滑水的作用機(jī)制及其數(shù)值實現(xiàn)

1.1 滑水的作用機(jī)制

滑水的作用機(jī)制是一個復(fù)雜的多學(xué)科交叉問題，其困難性導(dǎo)致目前對于滑水的針對性研究還比較少(Huang et al.，1997)，國內(nèi)更是未見報道(余斌， 2007； Liu et al.，2020)。國際上De Blasio et al.(2004)和Harbitz et al.(2003)采用改進(jìn)的潤滑理論分析了滑水觸發(fā)的力學(xué)機(jī)制，為后續(xù)研究奠定了理論基礎(chǔ)。

對海底滑坡的滑動過程進(jìn)行簡化分析，將滑坡體考慮為連續(xù)的黏彈性材料，假設(shè)一個平面滑坡體沿剛性海床滑動，滑坡體主要受到自身浮重和水體拖曳力的影響，其中自身浮重垂直于海床方向的分量控制滑坡體的沉積過程，自身浮重沿順坡方向的分量控制滑坡體向前滑動的動力過程?；麦w所受阻力主要包括基床阻力τB、環(huán)境水施加的拖曳阻力τD和前端動水壓力fp。其中基床阻力τB主要由基床土體抗剪強(qiáng)度提供，如圖1所示?；麦w頂部所受到拖曳力計算為：

(1)

式中：ρw是環(huán)境水的密度(ρw=1000 kg · m-3)；CF是一個常數(shù)，文獻(xiàn)(De Blasio et al.，2004)建議取值為0.003～0.005；u是滑坡體的滑動速度。前端動水壓力為：

(2)

式中：CP是一個常數(shù)，其取值可由文獻(xiàn)(Newman， 1977)查得。前端動水壓力部分主要影響滑坡體的形態(tài)，對滑坡體的動力特性影響較小，因此本文中不予考慮。

運動過程中，滑坡體前端趾部滯點處受到環(huán)境水的動水壓力最大，可計算為：

(3)

當(dāng)該動水壓力大于滑坡體前端的浮重時，滑坡體前端上翹，部分環(huán)境水侵入滑坡體底部，逐漸形成水墊層：

(4)

式中：ρ是滑坡體的密度；D是滑坡體厚度；g是重力加速度；θ是坡度。式(4)可變形為：

(5)

式中：Fc是一個表示滑水觸發(fā)時機(jī)的常數(shù)，文獻(xiàn)(De Blasio et al.，2004)根據(jù)一系列小尺度室內(nèi)試驗測得該常數(shù)通常為0.3，而對于現(xiàn)場真實滑坡案例可取值為1。式(4)表明，滑水發(fā)生時滑坡體的滑動速度必須大于某一個臨界速度。當(dāng)滑水發(fā)生后滑坡體底部的阻力由水墊層提供，因而其數(shù)值為：

(6)

式中：μw是水墊層的黏性，文獻(xiàn)(De Blasio et al.，2004)中取值為0.01；Dw是水墊層的厚度；fw是水墊層的消散量，用于表征水墊層維持的時間，本文中不予考慮。

圖 1 滑坡體受力示意圖Fig. 1 Schematic for sliding mass

1.2 滑水的物質(zhì)點法實現(xiàn)

物質(zhì)點法是任意拉格朗日-歐拉法的一種，采用拉格朗日物質(zhì)點描述材料的歷史信息(質(zhì)量、體積、密度、速度、變形梯度和應(yīng)力等)，并在歐拉網(wǎng)格上求解控制方程(de Vaucorbeil et al.，2020； Nguyen et al.，2020；Sun et al.，2020；Yuan et al.，2021)。物質(zhì)點法的主要計算步驟包括： (1)物質(zhì)點向網(wǎng)格節(jié)點投影； (2)計算網(wǎng)格節(jié)點速度和加速度； (3)更新物質(zhì)點的應(yīng)力和其他狀態(tài)信息； (4)更新物質(zhì)點位置。計算前，將滑坡體材料離散為一系列物質(zhì)點，將計算區(qū)域劃分成網(wǎng)格單元。計算時，首先將物質(zhì)點上的速度、質(zhì)量、體積和應(yīng)力等信息投影到網(wǎng)格節(jié)點上，在節(jié)點上建立平衡方程求解出下一時步的速度場，之后將新的速度和加速度投影給物質(zhì)點，用于更新其運動狀態(tài)和變形、應(yīng)力等信息。計算過程中，網(wǎng)格節(jié)點上的信息都為臨時信息，不保存任何歷史信息。網(wǎng)格在空間中的位置固定，不會像傳統(tǒng)有限元法單元那樣隨著材料的運動而發(fā)生變形，避免了網(wǎng)格畸變的問題，因此非常適合處理海底滑坡在滑動過程中的大變形問題(Soga et al.，2016；Yerro et al.，2016； 厲成陽等， 2018)。

本文中的分析基于自主開發(fā)的物質(zhì)點法程序MPM-GeoFluidFlow，基于顯式積分采用二次形函數(shù)進(jìn)行插值求解，采用柯西應(yīng)力描述應(yīng)力狀態(tài)，應(yīng)力增量采用Jaumann率表述以遵從有限應(yīng)變原理。該程序的特色功能有3項： (1)對程序進(jìn)行大規(guī)模GPU并行(Dong et al.，2015， 2018)，計算效率較CPU單核模擬最大提高約900倍； (2)采用先進(jìn)的接觸算法Geo-contact(Ma et al.，2014)，在土體-結(jié)構(gòu)物接觸的模擬中可以得到平滑的接觸反力； (3)對變形過大出現(xiàn)空隙或團(tuán)塊的物質(zhì)點重新生成(Dong， 2020)。

與傳統(tǒng)的CPU相比，GPU擁有更多的計算核心和更強(qiáng)的批處理能力，因此具有更好的并行加速潛力。然而，GPU并行的架構(gòu)與CPU并行差異較大，已有CPU并行的程序并不能直接移植到GPU平臺，因此，需要對程序架構(gòu)進(jìn)行重新設(shè)計，以最大程度發(fā)揮GPU計算核心的效力。并行過程中，需要特別注意保持各個計算線程之間的相互獨立性，防止出現(xiàn)顯存訪問競爭問題。GPU顯存到主機(jī)內(nèi)存之間數(shù)據(jù)傳輸速度較慢，因此，并行計算時，首先將所有變量傳輸?shù)紾PU的顯存中，以避免頻繁的變量傳輸。對于計算步(2)和(3)，每個節(jié)點的信息更新均采用一個核心進(jìn)行計算，所有核心的計算同時進(jìn)行，由于相互之間有較好的獨立性，因此具有良好的可并行性。類似的，對于計算步(4)，每個物質(zhì)點的信息更新采用一個核心進(jìn)行計算，并將所有核心同時進(jìn)行并行計算。但對于計算步(1)，傳統(tǒng)的線行計算中一般將不同的物質(zhì)點信息依次投影給周圍的節(jié)點，直接用到并行計算中時會產(chǎn)生線程間的相互干擾。本文程序采用的解決方案是首先建立每個網(wǎng)格節(jié)點對應(yīng)的物質(zhì)點列表，在并行計算中將每個節(jié)點的投影計算進(jìn)行線行計算，對不同的節(jié)點進(jìn)行并行計算(Dong et al.，2015， 2018)。

(7)

式中：Ai和mi分別是節(jié)點i所表征的面積和質(zhì)量； Δt是時間步長:

(8)

式中：G和λ分別是剪切模量和拉梅常數(shù)；d為網(wǎng)格尺寸；α為系數(shù)。

2 室內(nèi)水槽試驗

小尺度的室內(nèi)水槽試驗可簡化海底滑坡滑動過程中復(fù)雜的海底環(huán)境，對主要控制因素進(jìn)行精細(xì)化還原，因而可用于復(fù)現(xiàn)海底滑坡滑水效應(yīng)的發(fā)生過程。本文中模擬的室內(nèi)水槽試驗包含陸上滑坡(無水環(huán)境)和海底滑坡(有水環(huán)境)兩組對照試驗(De Blasio et al.，2004)。水槽在長度5.7 m內(nèi)坡度為5°； 在此距離(坡折點)以外坡度減小至1°?；麦w材料的初始形狀為高0.11 m，長2 m，截面積0.15 m2的橢圓形?；麦w材料表觀密度為1600 kg · m-3，抗剪強(qiáng)度與剪切應(yīng)變率關(guān)聯(lián)并考慮強(qiáng)度軟化的影響(范寧等， 2018):

(9)

陸上滑坡在滑坡體的自重(ρgcosD)驅(qū)動下向前滑動，并受到基底的摩擦阻力，最終滑動距離約為5 m，在坡折點(5.7 m)前停止，最終滑動形態(tài)如圖 2 所示。5 s時滑坡體的速度云圖如圖 3 所示。

圖 2 小尺度水槽試驗滑動距離Fig. 2 Small-scale flume test results

圖 3 陸上滑坡體速度云圖(時刻5 s)Fig. 3 Velocity contour for subaerial sliding mass(time 5 s)

對于水下滑坡，滑坡體在環(huán)境水的浮力作用下所受的自重驅(qū)動力((ρ-ρw)gcosD)較小。若不考慮滑水的作用，其滑動距離僅為約0.7 m (圖 2)。而在試驗中水下滑坡的滑動距離長達(dá)7.5 m，其主要原因是發(fā)生了滑水現(xiàn)象?？紤]不同的基底阻力τB，進(jìn)行參數(shù)分析，得到了不同的最終滑動距離S：τB=36 Pa，S=1.2 m;τB=20 Pa，S=6.4 m;τB=18 Pa，S=7 m。其中：當(dāng)τB=18 Pa時，滑坡體最終滑動距離與試驗中最為接近。由此可見，試驗中滑坡體所受到的基底阻力平均值約為18 Pa，該值遠(yuǎn)小于滑坡體的抗剪強(qiáng)度初始值。采用式(5)考慮滑水的觸發(fā)時機(jī)，其中，F(xiàn)c取值為0.3。當(dāng)滑坡體發(fā)生滑水時滑坡體底部的阻力由式(6)計算得到，最終預(yù)測的滑動距離與試驗較接近，為8.4 m。

通過上述對比可發(fā)現(xiàn)，物質(zhì)點法可較好地模擬海底滑坡的大變形滑動過程，能夠有效地復(fù)現(xiàn)滑水的觸發(fā)和發(fā)展過程，在預(yù)測滑坡體的滑動距離方面有較高的精度。相比之下，DAM方法由于對速度進(jìn)行深度平均的先天不足，無法準(zhǔn)確預(yù)測滑動過程中滑坡體的速度分布，在對具體案例進(jìn)行分析時限制性較大?？紤]海底滑坡過程中的滑水時，滑坡體前端滑動速度遠(yuǎn)大于滑坡體尾部，因而物質(zhì)點法模擬結(jié)果中滑坡體前后出現(xiàn)分離現(xiàn)象。該現(xiàn)象在試驗中也被觀測到(圖 2)，在真實滑坡中分離的滑坡體前端部分將會形成若干離散的塊體分布在滑坡的運動軌跡上。文獻(xiàn)(De Blasio et al.，2004)的DAM分析中僅考慮滑坡體為連續(xù)體，因而未能反映滑坡體的分離現(xiàn)象，由此可見，在海底滑坡的滑動模擬中，物質(zhì)點法比DAM方法更具優(yōu)勢，而DAM方法主要用于海底管線路由設(shè)計中粗略估計潛在海底滑坡影響范圍，且不能對海底滑坡和管線的相互作用進(jìn)行分析(王立忠等， 2008； Zakeri， 2009； Liu et al.，2015； Dong et al.，2017b， 2020a； 馮斌等， 2019； 王忠濤， 2019)?；瑒舆^程中(5 s)滑坡體的速度分布如圖 4 所示，滑動全過程中滑坡體的速度與位移時程變化趨勢如圖 5 所示。

圖 4 水下滑坡體速度云圖(時刻5 s)Fig. 4 Velocity contour for submarine sliding mass(time 5 s)

圖 5 滑坡體速度位移時程曲線Fig. 5 History of velocity and displacement of sliding mass

3 真實滑坡案例分析

一處真實的岸坡垮塌造成的海底滑坡發(fā)生在北挪威的Finneidfjord灣區(qū)，超過百萬立方米的土體發(fā)生了位移。關(guān)于事件后的沉積物性狀和滑動區(qū)環(huán)境特征等詳見文獻(xiàn)(Ilstad et al.，2004b)，調(diào)查結(jié)果顯示該滑坡主要由淺層氣體造成。整個滑坡區(qū)域可以劃分成4小塊，如圖 6 所示。其中沉積物主體分布在A區(qū)域，該區(qū)域在滑坡體原始位置的1 km范圍內(nèi)，坡度約2.86°； 在A區(qū)域以外，散布著不同尺寸的滑出塊體，其中最大的塊體尺寸110 m×60 m×2 m，體積約1.2萬立方米，最終在D區(qū)停止。區(qū)域B和C的坡度分別為0.94°和0.46°。對現(xiàn)場海床土體進(jìn)行鉆孔取樣并分析得出，土體的軟化系數(shù)為5～35，原狀土不排水抗剪強(qiáng)度為10～15 kPa。根據(jù)地質(zhì)條件綜合分析表明，滑水是造成滑坡體超遠(yuǎn)滑動并分離成離散滑塊的主要原因。

圖 6 Finneidfjord滑坡地形圖和沉積物抗剪強(qiáng)度(Ilstad et al., 2004a)Fig. 6 Morphology and sediment strength in Finneidfjord slide(Ilstad et al., 2004a)a. 滑坡區(qū)分為4個區(qū)域; b. 滑出塊體不排水抗剪強(qiáng)度

物質(zhì)點法分析中，原始滑坡體簡化為高度25 m，上底30 m，下底110 m的梯形平面體，單元尺寸0.25 m。原狀土強(qiáng)度15 kPa，軟化系數(shù)為10?；麦w材料浮密度為500 kg · m-3。模型底邊界模擬海床表面，初始為無滑動邊界，發(fā)生滑水時為部分滑動邊界。模型兩側(cè)和上邊界皆為自由邊界?；麦w的楊氏模量為500su0，泊松比為0.49，重力加速度g=9.81 m · s-2。由于本構(gòu)模型存在分叉不連續(xù)問題，因此對于具有強(qiáng)烈應(yīng)變軟化特性的土的計算往往會受到網(wǎng)格依賴性的困擾(圖 7)。本文的計算中發(fā)現(xiàn)，采用網(wǎng)格尺寸0.25 m和0.125 m時滑坡體前端滑出塊體的運動機(jī)制一致，對計算結(jié)果不會產(chǎn)生過大影響，如圖 7 所示。當(dāng)不考慮滑水效應(yīng)時，滑出塊體長約70 m，高約9 m?；鰤K體在A區(qū)(0～1 km)內(nèi)驅(qū)動力約為2.2 kPa，高于基底阻力(1.5 kPa)，因此滑塊始終保持加速，但速度較小，僅為0.2 m · s-1左右(圖 8a)； 當(dāng)滑塊進(jìn)入B區(qū)后(1～1.2 km)，其內(nèi)驅(qū)動力約為0.7 kPa，低于基底阻力(1.5 kPa)，因此滑塊始終保持減速，很快在B區(qū)停止?？紤]滑水效應(yīng)后，滑出塊體所受基底阻力主要由水墊層提供，其數(shù)值約為350 Pa，低于滑塊在A、B、C區(qū)中的驅(qū)動力，因此滑塊始終保持加速，最終在D區(qū)由于地形變化而停止?？紤]滑水效應(yīng)后，滑坡體最大滑動速度可達(dá)17 m · s-1(圖 8b)，并最終在C區(qū)中部停止(圖 9)。在滑出塊體滑動沿程，可見滑坡體變形產(chǎn)生的明顯滑動痕跡。

圖 7 不同網(wǎng)格尺寸對滑動距離和滑動形態(tài)的影響Fig. 7 Influence of mesh size on runout distance and morphologya. 網(wǎng)格尺寸0.25 m(1000 s); b. 網(wǎng)格尺寸0.125 m(1000 s)

圖 8 滑水對滑坡體滑動距離和滑動形態(tài)的影響Fig. 8 Influence of hydroplaning on runout distance and morphologya. 無滑水； b. 有滑水

圖 9 考慮滑水效應(yīng)時滑坡體的最終形態(tài)Fig. 9 Final morphology of sliding mass with hydroplaning

對比物質(zhì)點法計算結(jié)果和真實滑坡遺跡可發(fā)現(xiàn)，物質(zhì)點法模擬中考慮滑水效應(yīng)時滑坡體的最遠(yuǎn)滑動距離與真實距離較接近，進(jìn)一步表明物質(zhì)點法具有較高的準(zhǔn)確性和可靠性。盡管應(yīng)變軟化本構(gòu)模型具有網(wǎng)格依賴性的不足(其他數(shù)值方法中也同樣存在此問題)，不能準(zhǔn)確捕捉滑坡體的最終形態(tài)，但是前端滑塊的滑動特性較接近，最終預(yù)測的滑塊最遠(yuǎn)距離有較高的可信度。對比分別考慮和不考慮滑水效應(yīng)的模擬結(jié)果發(fā)現(xiàn)，滑水效應(yīng)是Finneidfjord滑坡中滑坡體具有高動力性能的主要原因，進(jìn)一步驗證了前人研究中關(guān)于滑水效應(yīng)的重要性和危害性的分析。

4 結(jié) 論

物質(zhì)點法是一種新型數(shù)值方法，通過追蹤物質(zhì)點的移動描述材料的變形過程，可以有效地模擬海底滑坡等大變形問題。本文采用自主研發(fā)的物質(zhì)點法程序MPM-GeoFluidFlow模擬海底滑坡的滑動過程，考慮滑水對基床阻力的影響，分別對室內(nèi)試驗和真實的海底滑坡案例進(jìn)行了反分析。主要結(jié)論如下：

(1)通過對比有滑水和無滑水條件的計算結(jié)果，驗證了滑水是海底滑坡超遠(yuǎn)距離和超快速度滑動的主要原因，并據(jù)此評估了滑水對海底滑坡危害范圍和滑動形態(tài)的影響，為實際工程中的應(yīng)用提供了便利。

(2)計算結(jié)果進(jìn)一步驗證了物質(zhì)點法在模擬海底滑坡的滑動過程時的可靠性和精確度，為相關(guān)類似研究提供了參考。