波浪作用下海床孔壓累積過程離散元數(shù)值模擬

王岳，劉春*，劉曉磊，劉輝，李亞沙

( 1. 南京大學(xué) 地球科學(xué)與工程學(xué)院，江蘇 南京 210023；2. 中國海洋大學(xué) 環(huán)境科學(xué)與工程學(xué)院，山東 青島 266100)

1 引言

黃河口海床沉積物工程地質(zhì)特性復(fù)雜，土體以粉質(zhì)土為主，易液化。在風(fēng)暴潮期間，波浪對土體施加巨大的循環(huán)荷載作用，造成海床土體孔壓累積，并發(fā)生一定程度的液化，進而誘發(fā)海床表層出現(xiàn)滑坡、塌陷凹坑等災(zāi)害，直接威脅到海上工程如石油平臺、海底管線等構(gòu)筑物的安全穩(wěn)定，最終造成巨大的經(jīng)濟損失[1]。

黃河口海床沉積物孔壓累積問題已成為國內(nèi)外海洋工程及地質(zhì)工程領(lǐng)域研究的熱點問題。已有大量學(xué)者對波浪作用下海床孔壓響應(yīng)機制進行了研究，Anderson等[2]基于波浪水槽實驗對沖浪帶沙洲上波浪引起的孔壓梯度和床位響應(yīng)進行了觀測研究；Niu等[3]在寬波浪水槽中進行了一系列實驗，研究了隨機波浪作用下粉質(zhì)海底的孔壓響應(yīng)；常方強和賈永剛[1]在黃河口海床若干深度處埋設(shè)孔壓探頭，研究黃河口原狀和重塑粉質(zhì)土在波浪循環(huán)荷載作用下的液化特征；劉曉磊等[4]在現(xiàn)代黃河水下三角洲潮間帶岸灘進行現(xiàn)場試驗，測定黃河口原狀海床沉積物在循環(huán)荷載作用不同階段的孔壓響應(yīng)與強度變化；張少同等[5]研究了黃河現(xiàn)代三角洲沉積物的堆積、固結(jié)、液化、侵蝕再懸浮、海床變形滑動及后期改造等一系列動態(tài)變化過程；呂豪杰等[6]基于改造的圓筒試驗設(shè)備，通過加載裝置控制模型土體上方的波壓力，實現(xiàn)了波浪荷載作用下的樁周土體孔壓響應(yīng)模擬；宋玉鵬等[7]在黃河口使用自行研發(fā)的孔壓監(jiān)測設(shè)備對海底粉土孔壓力進行了有效監(jiān)測。

在試驗研究的基礎(chǔ)上，數(shù)值模擬是探究海床孔壓響應(yīng)機理的重要方法。例如，吳雷曄等[8]通過與離心模型試驗數(shù)據(jù)的對比，驗證了基于Biot固結(jié)理論及“交變移動”土體本構(gòu)模型的數(shù)值分析方法對波浪作用下海床液化問題模擬的有效性；潘冬子等[9]結(jié)合波浪槽模型試驗研究的結(jié)果，提出了一種粉質(zhì)海床在波浪作用下超靜孔壓增長的有限差分?jǐn)?shù)值分析方法，并揭示了粉質(zhì)海床孔壓發(fā)展的機理；劉濤等[10]采用有限差分法，對波浪作用下孔壓響應(yīng)進行了數(shù)值計算，發(fā)現(xiàn)海床中由波浪引起的超孔壓隨深度分布不均勻，存在極值點，并且孔壓響應(yīng)存在滯后效應(yīng)；周援衡等[11]采用Biot固結(jié)理論分析了海床沙層對波浪荷載的響應(yīng)，用有限元數(shù)值模擬法得到沙床中孔壓分布，數(shù)值計算結(jié)果與實測結(jié)果吻合良好。

前人對波浪作用下海床沉積物孔壓響應(yīng)問題開展了大量的試驗研究，并采用有限差分法、有限元法等數(shù)值方法，較好地模擬了孔壓的累積過程。但是這些方法難以模擬巖土體的離散性以及孔壓累積造成的土體破壞和災(zāi)害效應(yīng)。本文基于非連續(xù)的離散元法和多孔介質(zhì)孔隙網(wǎng)絡(luò)模型，提出離散元孔隙密度流方法，實現(xiàn)滲流和流固耦合模擬?；诂F(xiàn)場試驗裝置和數(shù)據(jù)，建立了相應(yīng)的離散元數(shù)值模型，對波浪作用下海床孔壓累積過程進行了模擬分析，探究了孔壓累積規(guī)律和機理。并進一步模擬了周期性波浪荷載作用下土顆粒液化懸浮過程。

2 離散元法基本原理

在離散元法中，巖土體模型由一系列遵循牛頓運動定律的顆粒堆積組成（圖1）。顆粒間通過彈簧力相互作用，兩顆粒間的法向力（Fn）由下式給出[12]：

圖1 離散顆粒堆積模型（a）、顆粒間法向彈簧力（b）和顆粒間切向彈簧力（c）Fig. 1 Packing model of discrete element (a), normal spring force between particles (b), and shear spring force between particles (c)

式中，Kn為法向剛度；Xn為法向相對位移；Xb為斷裂位移。初始時，顆粒與其相鄰顆粒相互連接，受拉力或壓力作用（式（1a），拉力為正）。當(dāng)兩顆粒之間的Xn超過斷裂位移時，彈簧斷裂，顆粒間拉力消失（式（1c）），僅存在壓力作用（式（1b））。

顆粒間通過切向彈簧來模擬剪切變形和切向力，剪切彈簧的切向力（Fs）通過下式?jīng)Q定：

式中，Ks為切向剛度；Xs為切向相對位移。對于完整的顆粒連接，其最大切向力（Fsmax）由摩爾庫倫準(zhǔn)則確定：

式中，F(xiàn)s0為顆粒內(nèi)部的抗剪力；μp為摩擦系數(shù)。當(dāng)顆粒受切向外力超過式（3）中的Fsmax時，顆粒間連接斷裂，單元連接允許的最大切向力F'smax為

當(dāng)外力超過最大摩擦力時，兩相鄰顆粒開始滑移，滑動摩擦力等于F'smax。

在數(shù)值模擬中，單元所受合力為彈簧力、阻尼力和重力等的矢量和。假定在非常小的時間步內(nèi)，單元的速度和加速度不變，基于時間步迭代算法和牛頓運動方程實現(xiàn)離散單元的動態(tài)模擬[12]。

3 離散元孔隙密度流法

基于傳統(tǒng)的巖體孔隙網(wǎng)絡(luò)模型[13]，本文提出了離散元孔隙密度流法來實現(xiàn)巖土體流固耦合的高效計算，這種方法的基本步驟為：（1）利用離散元堆積顆粒建立孔隙流體域（圖2b）；（2）通過孔隙流體密度和溫度來確定孔壓，在孔隙間采用類似于達(dá)西定律的方法計算滲流；（3）將孔壓作用于固體顆粒，實現(xiàn)流固耦合模擬。

如圖2a所示，經(jīng)過離散元計算獲得離散元堆積顆粒后，根據(jù)顆粒間的接觸關(guān)系，通過連接相鄰顆粒的中心點，在堆積模型中剖分出一系列相互連通的流體域。建立以孔隙流體域單元為節(jié)點，孔喉通道相連接的孔隙流體域（圖2b），最后得到顆粒?孔隙系數(shù)（圖2c）。

圖2 離散元堆積顆粒（a）、孔隙流體域（b）和顆粒?孔隙系統(tǒng)（c）Fig. 2 Discrete element stacked particles (a), pore fluid domain (b), and particle and pore system (c)

數(shù)值模擬假定孔隙均為飽和狀態(tài)，孔隙流體壓力ρ可以由孔隙流體密度ρ和溫度T來確定[14]：

根據(jù)已知的飽和水蒸氣壓、溫度、密度數(shù)據(jù)擬合得到水溫為20℃時，孔隙流體壓力與孔隙流體密度的關(guān)系[15]為

如圖3所示，當(dāng)兩相鄰孔隙間存在壓力差時，流體將通過二者間的孔喉通道發(fā)生滲流，其滲流量可通過類似達(dá)西定律的方法來計算。本文中，在單位時間內(nèi)通過孔喉的流量（q）定義為

圖3 顆粒微觀滲流示意圖Fig. 3 Schematic diagram of particle microscopic seepage

式中，k為孔喉的滲透系數(shù)；A為孔喉通道面積；dP為孔喉間壓力差；l為孔喉長度。兩顆粒半徑分別為R1、R2；兩顆粒中心的距離為L；孔喉長度l定義為

對于二維問題，單位寬度孔喉通道面積可用孔喉的直徑（dw）表示。由于二維顆粒堆積通常會封閉孔喉，從而使孔喉直徑為0，導(dǎo)致無法發(fā)生滲流。因此，在進行二維模型滲流計算時，需要給每個單元定義較小的滲流作用半徑（Rw），默認(rèn)取單元半徑的0.975倍，從而保證一定的孔喉直徑，并使流體能夠通過孔喉運移。孔喉直徑（dw）定義為

式中，Rw1和Rw2分別為單元1和單元2的滲流作用半徑。

以上計算中，流體在孔隙間壓差驅(qū)動下運移，實現(xiàn)了滲流模擬，進一步實現(xiàn)流固耦合作用，將流體壓力作用在單元上，并獲得流體對單元的合力。通過常規(guī)離散元法計算單元運動，實現(xiàn)流體對固體的作用。同時，當(dāng)單元運動后，孔隙的體積變化，需計算孔隙密度，并根據(jù)公式（5）獲得孔隙新的壓力，從而實現(xiàn)固體對流體的作用。最后，根據(jù)單元坐標(biāo)劃分新的網(wǎng)格，并更新網(wǎng)格中的流體屬性，再返回到滲流計算。

4 離散元建模與數(shù)值模擬

4.1 建模與賦材料

本文依據(jù)常方強和賈永剛[1]在黃河三角洲現(xiàn)行的清水溝流路南部的潮間帶上開展的現(xiàn)場試驗，實現(xiàn)了波浪作用下海床孔壓累積過程的離散元建模與數(shù)值模擬。現(xiàn)場試驗裝置主要包括加載系統(tǒng)和孔壓監(jiān)測系統(tǒng)，加載系統(tǒng)采用自制活塞式壓水方式，主要由圓柱形鐵桶、活塞和手柄組成；孔壓監(jiān)測系統(tǒng)由孔壓探頭、多通道數(shù)據(jù)采集系統(tǒng)、通訊電纜和計算機組成。

根據(jù)現(xiàn)場試驗裝置，利用MatDEM進行建模[16]。第一步先建立長、高為0.5 m×1 m的模型箱，向模型箱內(nèi)隨機添加平均半徑為0.005 m的顆粒，總顆粒數(shù)為5629，顆粒經(jīng)過重力沉積和壓實作用，最終形成顆粒堆積模型（圖4）。

圖4 顆粒堆積模型Fig. 4 Particle accumulation model

第二步再給模型賦材料，根據(jù)現(xiàn)場試驗[1]及經(jīng)驗得到海洋土體宏觀力學(xué)性質(zhì)（表1），基于宏微觀轉(zhuǎn)換公式[12]和MatDEM自動訓(xùn)練模塊，軟件自動訓(xùn)練得到符合宏觀要求的離散元材料，具體離散元微觀力學(xué)參數(shù)見表2。

表1 海洋土體宏觀力學(xué)性質(zhì)Table 1 Macro mechanical properties of ocean soil

表2 離散元微觀力學(xué)參數(shù)Table 2 Micro mechanical parameters of the discrete element

第三步將模型箱上部0.2 m范圍內(nèi)的顆粒篩選出并刪除，建立海床表面的流體域。根據(jù)土體滲透率確定微觀滲透系數(shù)，并計算出模擬海床沉積物滲透系數(shù)K為1.32×10?5cm/s，數(shù)值在現(xiàn)場試驗測得的數(shù)據(jù)范圍內(nèi)。

4.2 數(shù)值模擬計算

現(xiàn)場試驗利用人工操作活塞施加水壓循環(huán)荷載以模擬波浪荷載，循環(huán)荷載控制在0～4 kPa，加載周期為6.67 s，循環(huán)加載40次，觀察孔壓累積規(guī)律。同樣地，在數(shù)值模擬中，對海床流體域施加0～4 kPa循環(huán)荷載，模擬波峰、波谷的作用。模擬設(shè)定時間步為6.67×10?4s，每個周期迭代計算10000次，實現(xiàn)加載周期6.67 s。循環(huán)加載40次，模擬現(xiàn)實時間266.8 s，其時間段與現(xiàn)場試驗分析時間段一致。

現(xiàn)場試驗分別將孔壓探頭埋設(shè)在深度為5 cm、8 cm、12 cm、20 cm和30 cm處，因此在離散元模型內(nèi)選取土體深度為0.05 m、0.08 m、0.12 m、0.20 m、0.30 m、0.40 m和0.5 m的位置，對應(yīng)于現(xiàn)場試驗，記錄不同深度土體孔隙水壓力，并進一步用于繪制孔壓分布圖和孔壓累積過程圖。

5 模擬結(jié)果與分析

5.1 單周期波浪荷載下的孔壓變化特征

圖5顯示了施加荷載1個周期內(nèi)，不同深度土體孔隙水壓力隨時間變化曲線。圖中黑色實線代表數(shù)值模擬時，施加于海床表面流體域上即深度為0處用來模擬波浪作用的荷載，增加循環(huán)次數(shù)可模擬正弦荷載作用。

圖5 單周期荷載下孔隙水壓力變化趨勢Fig. 5 Trend of pore water pressure under single-cycle loading

前T/2時間內(nèi)施加最大荷載4 kPa，在荷載作用下，各深度土體顆粒間逐漸產(chǎn)生孔隙水壓力。深度為0.05 m時土體最先產(chǎn)生孔壓，由于波浪直接作用于表層土體，而水下土體處于完全飽和狀態(tài)，當(dāng)受到荷載作用時，表層土顆粒間的水最先受到擠壓，并且不能及時排出，從而導(dǎo)致表層土體最先產(chǎn)生孔壓。在滲流作用下，孔壓逐漸向深層傳遞，由于滲流需要一定的時間，因此土體深度越大，滲流時間所需越長，所以，土體孔壓增大速率隨深度增加而降低。

后T/2時間內(nèi)施加最小荷載0 kPa，觀察到深度為0.05 m時土體顆粒間孔壓急劇減小，緊接著深度為0.08 m、0.12 m的土體顆粒間孔壓開始陸續(xù)減小，且深度越大，減小越緩慢。施加荷載1個周期內(nèi)，深度大于0.2 m的土體，孔壓一直緩慢上升，因為土體深度越大，顆粒間孔壓受荷載變化影響需要的時間越長，在6.67 s周期內(nèi)表現(xiàn)不明顯。

施加荷載1個周期后，觀察到深度為0.12 m的土體顆粒間孔壓最大，表層土體顆粒與波浪直接接觸，受荷載作用最明顯，在0 kPa荷載作用時孔壓減小速率快；荷載減小對深層土體產(chǎn)生影響所需的時間更長，并且上層土體顆粒間孔壓減小過程中仍然發(fā)生滲流作用，因此深層土體顆粒間孔壓減小速率相對較慢。

5.2 循環(huán)波浪荷載下的孔壓累積過程

為了解周期循環(huán)波浪荷載的作用，將施加荷載次數(shù)設(shè)為40次。圖6顯示了施加荷載第1個、第10個、第20個和第40個周期后，數(shù)值模型中的孔壓分布情況，其分別對應(yīng)于6.67 s、66.7 s、133.4 s和266.8 s 4個時刻。從圖中可以看出，隨著荷載循環(huán)次數(shù)增加，孔壓累積范圍逐漸變寬，孔壓最大區(qū)域穩(wěn)定在土體深度0.20 m上下。加載40個周期后，深度大于0.30 m土體孔壓基本穩(wěn)定在1.5 kPa。

圖6 荷載循環(huán)1次（a）、10次（b）、20次（c）、40次（d）孔壓分布Fig. 6 Pore pressure distribution of load cycle 1 (a), 10 (b), 20 (c), 40 (d) times

圖7為不同深度土體在荷載循環(huán)40個周期內(nèi)，孔壓隨時間的累積曲線。從圖中可以看出，深度從0.05 m到0.30 m處的孔壓均呈現(xiàn)出“邊振蕩，邊累積”的趨勢，振蕩幅值隨著深度的增加而減小。隨著孔壓累積時間延長，各個深度土體孔壓都累積增加[17]，增加速率逐漸降低。

圖7 循環(huán)荷載40次孔壓累積Fig. 7 Pore water pressure accumulation process under cyclic loading 40 times

將孔壓累積過程曲線進行擬合，可得到不同土體深度處孔壓P與時間t之間的函數(shù)關(guān)系。土體深度為0.05 m時：

土體深度為0.12 m時：

土體深度為0.2 m時：

土體深度為0.3 m時：

以上公式可歸納為

式中，A為振幅；w為角速度；φ為初相；B為系數(shù)；α為指數(shù)。

可以看到，孔壓累積過程曲線為正弦函數(shù)與冪函數(shù)的疊加。正弦波的周期與加載周期基本吻合，正弦波振幅A、系數(shù)B和土體深度呈負(fù)相關(guān)，冪函數(shù)指數(shù)α和土體深度呈正相關(guān)。

表3顯示了采用數(shù)值模擬與現(xiàn)場試驗兩種方法得到不同深度土體孔壓累積266.7 s時的數(shù)據(jù)對比，此時孔壓數(shù)據(jù)都基本達(dá)到穩(wěn)定。結(jié)果顯示，數(shù)值模擬與現(xiàn)場試驗孔壓數(shù)據(jù)接近；累積時間266.7 s后，兩者孔壓都緩慢上升，呈現(xiàn)出良好的相關(guān)性。分析存在細(xì)小誤差的原因為：數(shù)值模擬基于平面應(yīng)變假設(shè)，建立的離散元模型相對均質(zhì)，而現(xiàn)場試驗海床土體較為復(fù)雜，可能存在分層，并且其滲透性、密實度等會隨深度和孔壓的變化而發(fā)生變化，因此孔壓向深層傳遞過程會有所差異。導(dǎo)致0.12 m深度誤差較大的原因可能是現(xiàn)場試驗數(shù)據(jù)采用人工加載、人工讀數(shù)，存在偶然誤差。

表3 數(shù)值模擬與現(xiàn)場試驗孔壓數(shù)據(jù)對比Table 3 Data comparison of pore pressure between numerical simulation and field test

繼續(xù)增加施加荷載次數(shù)，當(dāng)孔壓累積時間足夠長時（500 s以上）可以發(fā)現(xiàn)0.20 m深度的土體孔壓最終收斂于所施加最大荷載與最小荷載的平均值（2 kPa）。現(xiàn)場試驗測得土體重度為19.8 kN/m3，則土體浮重度r'為10 kN/m3，當(dāng)施加荷載為0 kPa時，0.20 m深度土體的初始有效應(yīng)力σ'=r'×h，結(jié)果為2 kN/m2。傳統(tǒng)土體力學(xué)以初始有效應(yīng)力作為判斷海床液化的臨界值[18]，在離散元法中同樣適用，此時土體孔壓將達(dá)到初始有效應(yīng)力，土體開始發(fā)生液化，液化一方面會降低海床表層沉積物的抗侵蝕性；另一方面，在海床內(nèi)部垂向滲流的驅(qū)動下，內(nèi)部部分細(xì)粒沉積物會被“泵送”輸運至海床表面，進而進入上覆水體成為再懸浮沉積物[5]。

6 海床沉積物液化懸浮模擬

張少同等[5]分析了黃河三角洲沉積物的侵蝕再懸浮過程，波浪循環(huán)荷載作用不僅可以導(dǎo)致海床沉積物發(fā)生孔壓累積、液化破壞、強度弱化甚至喪失，還為沉積物再懸浮提供剪切力，微弱的近底剪切力就能夠?qū)姸热趸蛲耆夯暮５壮练e物侵蝕再懸浮。

以上數(shù)值模擬研究表明，在波浪循環(huán)荷載作用下，海床沉積物內(nèi)部孔隙水壓力發(fā)生累積，并逐漸達(dá)到液化標(biāo)準(zhǔn)，進而土體顆粒發(fā)生液化懸浮。為了進一步分析土體發(fā)生液化后的情況，開展了波浪周期性荷載作用下，海床沉積物液化懸浮過程的模擬。

6.1 建模與加載

模擬過程中離散元建模方法同4.1節(jié)，為了更接近海海床實際情況，設(shè)置模型箱長、高為1 m×0.5 m，隨機添加平均半徑為0.01 m的顆粒，總顆粒數(shù)為1579。加載方式與之前不同，在海床流體域左側(cè)施加4 kPa荷載，相當(dāng)于實際波高0.80 m的作用效果[4]，線性遞減到0 kPa模擬波浪剪切作用。根據(jù)5.2節(jié)模擬結(jié)果，在土體0.20 m深度處施加2 kPa荷載，用來模擬已累積孔壓對后續(xù)土體孔壓累積以及顆粒移動的影響。模擬設(shè)定時間步為1×10?3s，每個周期迭代計算10000次，模擬真實世界10 s，得到土體孔壓分布和顆粒位移場。

6.2 模擬結(jié)果分析

從圖8中可以看出，在波浪剪切作用和已累積孔壓共同作用下左側(cè)表層土體內(nèi)積累的孔壓較大，為3.5 kPa左右，右側(cè)較小，為0.5 kPa左右，深層土體內(nèi)孔壓累積大小介于兩者之間。

圖9為數(shù)值模擬10 s后土體顆粒位移場，可以觀察到深度在0～0.30 m范圍內(nèi)土體顆粒發(fā)生剪切位移，結(jié)合圖8孔壓分布可以發(fā)現(xiàn)，0～0.30 m土體內(nèi)孔壓累積超過2 kPa，達(dá)到液化標(biāo)準(zhǔn)，土體發(fā)生液化，在波浪剪切力的作用下，土顆粒懸浮，淺層土體由于直接受到波浪剪切作用，位移量更大，隨著土體深度增大，位移量逐漸減小。

圖8 孔壓分布Fig. 8 Distribution of pore pressure

圖9 模擬10 s位移場Fig. 9 Simulation of 10 s displacement field

觀察整個模型箱內(nèi)土體顆粒位移情況，發(fā)現(xiàn)海床整體呈現(xiàn)圓弧形移動（圖10）。土體液化后土體顆粒間的黏聚力喪失，進而造成海床抗剪性降低，在波浪荷載作用下海床發(fā)生圓弧形振蕩剪切。由于波浪荷載是循環(huán)荷載，呈周期性變化，波峰和波谷交替出現(xiàn)，在海床內(nèi)引起的循環(huán)應(yīng)力呈現(xiàn)相反的方向變化，為振蕩剪切破壞創(chuàng)造了條件。隨著波浪周期性荷載的持續(xù)作用，土體液化深度增大，振蕩面向更深層延展，最終破壞穩(wěn)定在一定深度[15]。

圖10 圓弧形移動示意圖Fig. 10 Schematic diagram of circular arc movement

7 結(jié)論

本文采用離散元數(shù)值模擬的方法，提出離散元孔隙密度流方法，實現(xiàn)滲流和流固耦合模擬?；诂F(xiàn)場試驗裝置和數(shù)據(jù)，通過自主研發(fā)的離散元模擬軟件MatDEM，探究了波浪荷載作用下海床沉積物的孔壓累積規(guī)律和機理，并分析了顆粒液化懸浮過程，取得了以下結(jié)論：

（1）提出了離散元孔隙密度流方法，并基于非連續(xù)的離散元法和多孔介質(zhì)孔隙網(wǎng)絡(luò)模型，實現(xiàn)了海床沉積物孔壓的累積過程模擬，模擬結(jié)果與現(xiàn)場試驗吻合較好。

（2）對海床沉積物施加波浪荷載后，表層土體中產(chǎn)生較高孔壓，并逐漸向深層傳遞；在循環(huán)波浪荷載作用下，土體顆粒間孔壓累積范圍逐漸增加。

（3）當(dāng)孔壓累積時間足夠長時，土層中孔壓收斂于所施加最大荷載與最小荷載的平均值，此時若孔壓達(dá)到初始有效應(yīng)力，土體將發(fā)生液化，內(nèi)部土體顆粒成為再懸浮沉積物。

（4）土體顆粒液化懸浮后發(fā)生移動，淺層顆粒位移量大，隨著土體深度增大，位移量逐漸減小，土體整體表現(xiàn)為圓弧形移動。

本研究中模擬了土體顆粒間孔壓的累積變化情況，在此基礎(chǔ)上，初步分析了土體液化懸浮過程，驗證了基于孔隙密度流法的流固耦合模擬的有效性，關(guān)于土體再懸浮過程以及振蕩剪切破壞的精細(xì)分析和機理研究有待后續(xù)開展。