波流作用下黃河三角洲硬殼層液化滲流形成機制研究*

冷 浩 胡瑞庚 劉紅軍② 王兆耀

(①中國海洋大學(xué)環(huán)境科學(xué)與工程學(xué)院， 青島 266100， 中國) (②山東省海洋環(huán)境地質(zhì)工程重點實驗室， 青島 266100， 中國)

0 引 言

黃河三角洲地區(qū)普遍發(fā)育著表面硬殼層，又稱“鐵板砂”，該表面硬殼層具有密度大、顆粒粗、滲透性大等特點，厚度一般1～3 m左右，在床面分布不連續(xù)。由于其上述特點，黃河三角洲表面硬殼層常作為海底管線和海洋平臺的持力層，所以對硬殼層的形成機制以及其在波浪作用下的性質(zhì)變化研究具有重要意義(賈永剛等， 2007)。

賈永剛等(2007)通過在黃河三角洲海床表面進行人工振動模擬波浪循環(huán)荷載以及室內(nèi)模型試驗，每輪荷載施加期間以及停止作用后的一段時間內(nèi)海床表層一定范圍內(nèi)的土體發(fā)生液化，孔壓消散后進行下一輪荷載，長期循環(huán)荷載作用下海床表層土體呈現(xiàn)重復(fù)液化，研究發(fā)現(xiàn)海床重復(fù)液化的結(jié)果導(dǎo)致了海床強度增強并且一定范圍內(nèi)的細顆粒物質(zhì)向海床表面發(fā)生移動，海床表層土體密度增大。Yang et al. (2021)在黃河三角洲的新灘和廣利港地區(qū)進行了原位靜力觸探試驗，海床表層土體中的錐尖阻力qt和側(cè)摩阻力fs數(shù)值較大(圖 1)，表明該地區(qū)海床表層土體具有較高的土體強度，常方強等(2009)取黃河三角洲硬殼層存在區(qū)域原位1 m深度的土體進行土工試驗，研究了硬殼層土體的粒度分布、孔隙比、超固結(jié)比等，砂粒、粉粒含量較高，砂粒含量約為16%，粉粒含量在44%～87%之間，而且該層土體相較于下伏土層處于明顯的超固結(jié)狀態(tài)，由此表明黃河三角洲硬殼層土體強度高，密度大，顆粒粗等特點。王虎等(2019)通過從原生沉積、沖刷分選和液化滲流3種模式對國內(nèi)外鐵板砂形成機制研究進行了總結(jié)分析，認為黃河三角洲鐵板砂分布不連續(xù)的特點以及現(xiàn)場試驗研究都表明液化滲流模式導(dǎo)致了黃河三角洲硬殼層的形成與發(fā)展。

圖 1 原位CPT測量結(jié)果Fig. 1 In-situ CPT measured results

波浪作用下會在海床內(nèi)部產(chǎn)生附加應(yīng)力，循環(huán)附加應(yīng)力會引起土體單元產(chǎn)生塑性應(yīng)變，這種應(yīng)變使得土體產(chǎn)生不可恢復(fù)的變形體積壓縮，由于土體內(nèi)部孔隙水壓力來不及排出而不斷累積，形成超孔隙水壓力，當(dāng)超孔隙水壓力超過該點處的有效應(yīng)力時土體發(fā)生完全液化(王虎， 2012； Hu et al.，2021)。

根據(jù)超孔隙水壓力發(fā)展模式，波致海床液化可分為兩種機制： (1)瞬態(tài)液化，在波谷位置，土骨架受到向上的滲流作用力，當(dāng)向上的滲流作用力大于豎向有效應(yīng)力時，發(fā)生瞬態(tài)液化(Sakai et al.，1992; Tsai， 1995;Jeng et al.，2015； 李正輝， 2019)； (2)累積液化，波浪作用下會在海床內(nèi)部產(chǎn)生附加應(yīng)力，循環(huán)附加應(yīng)力會引起土體單元產(chǎn)生塑性應(yīng)變，這種應(yīng)變使得土體產(chǎn)生不可恢復(fù)的變形體積壓縮，由于土體內(nèi)部孔隙水壓力來不及排出而不斷累積，形成超孔隙水壓力，當(dāng)超孔隙水壓力超過該點處的有效應(yīng)力時土體發(fā)生完全液化(劉占閣， 2008； Sumer et al.，2012； 王虎， 2012； 廖晨聰， 2016； 段釗等， 2020； 劉曉磊等， 2020; Hu et al.，2021)。在海洋環(huán)境中存在海流作用時，波浪與海流耦合作用會引起海床表面波壓力的變化，進而影響海床內(nèi)部累積孔隙水壓力的大小，順流條件下會導(dǎo)致孔隙水壓力累積幅值增大； 逆流條件則相反(劉小麗等， 2018)。本文以數(shù)值模擬方法研究波流聯(lián)合作用下黃河三角洲海床內(nèi)部孔隙水壓力的累積液化特征。

海床表面為排水界面，由波浪作用而產(chǎn)生的超孔隙水壓力通過排水邊界消散，從而形成穩(wěn)定的垂直于床面的滲流場，海床內(nèi)部土顆粒受到重力、水流拖拽力、上舉力以及滲流力的作用(圖 2)。滲流力會顯著降低泥沙啟動的臨界切應(yīng)力，當(dāng)土顆粒受到的切應(yīng)力大于臨界切應(yīng)力時泥沙開始啟動(Cheng et al.，1999； 程永舟等， 2012； 曹志剛等， 2019)。細顆粒物質(zhì)受上述作用力影響明顯首先被啟動上浮，通過粗顆粒之間的孔隙被帶到海床表面并逐漸發(fā)展成滲流通道，導(dǎo)致海床表層一定范圍內(nèi)的土顆粒粗化(圖 3)。與此同時海床受到波致循環(huán)剪應(yīng)力作用產(chǎn)生的塑性變形使海床表層土體密化，密化和粗化作用兩者聯(lián)合作用于海床從而形成了鐵板砂。王虎等(2014)通過將滲流力考慮到傳統(tǒng)泥沙啟動公式中：推導(dǎo)并驗證了波浪作用下海床臨界沖刷深度的計算方法，有效解釋了波浪對海床的改造作用。

圖 2 泥沙顆粒受力情況Fig. 2 Stress of sediment particles

圖 3 液化滲流模式下硬殼層的形成(王虎等， 2019)Fig. 3 Formation of the upper seabed layer under liquefaction seepage model(Wang et al.，2019)

基于土力學(xué)與泥沙運動學(xué)相結(jié)合，本文以黃河三角洲粉土海床為研究對象，采用數(shù)值模型研究波浪與海流同一時間共同作用下海床內(nèi)部孔隙水壓力累積特征，分析不同流速以及不同硬殼層厚度對累積孔隙水壓力的影響。根據(jù)數(shù)值計算求得的累積孔隙水壓力沿深度的分布值，基于王虎推導(dǎo)建立的海床臨界沖刷深度計算方法，研究考慮滲流力作用的海床泥沙啟動以及臨界沖刷深度，分析不同波流條件對鐵板砂形成深度的影響。

1 數(shù)值模型建立與驗證

1.1 控制方程

利用數(shù)值模擬軟件COMSOL來對波浪與海流同時刻共同作用于海床時的累積孔隙水壓力響應(yīng)進行計算。累積孔隙孔壓響應(yīng)計算模型控制方程和邊界條件如下。

本文基于biot固結(jié)理論和滲流連續(xù)方程，同時考慮土骨架和流體的壓縮性質(zhì)，海床土體假定為多孔彈性介質(zhì)，海床土體在x、z方向上的二維控制方程為

(1)

(2)

(3)

(4)

(5)

式中：u，w分別為土體的水平和豎向位移;εV為土單元體應(yīng)變；pt為波致孔隙水壓力；β為流體壓縮系數(shù)；G為土體的壓縮模量；υ為泊松比；k為滲透系數(shù)；γ為海床土重度；n為孔隙度；τ為土單元剪應(yīng)力；γW為海水重度。

(6)

式中：K為孔隙流體體積模量；pW為絕對靜水壓力；Sr為飽和度。

海床累積孔隙水壓力發(fā)展控制方程為

(7)

式中：pres為累積孔隙水壓力；CV為土體固結(jié)系數(shù)；f為孔隙水壓力累積源項。

采用已建立的適用于黃河三角洲海床的孔隙水壓力累積源項f(胡瑞庚等， 2021)：

(8)

(9)

(10)

式中：N為動荷載循環(huán)次數(shù)；Nl土體達到液化時的循環(huán)次數(shù)；a，b，c為經(jīng)驗系數(shù)，根據(jù)動三軸試驗分別得到a=0.326，b=0.8，c=0.75(劉紅軍等， 2005)。ar，βr為土體類型與密度相關(guān)參數(shù)，根據(jù)試驗分別取0.221， 0.312(Sumer et al.，2012)。σ′0為平均有效應(yīng)力，當(dāng)累積孔隙水壓力超過該值時海床即發(fā)生液化。

(11)

式中：K0為土體側(cè)壓力系數(shù)；γ′為土體浮容重；z為海床深度。

1.2 邊界條件

對于波流共同作用在海床表面生成的波壓力，采用hsu et al.(2009)提出的解析解：

(12)

式中：ρW為海水密度；H為波高；d為水深；λ為波數(shù)；U0為海流流速；ω，ω0，ω2分別為總圓頻率、第一階圓頻率和第三階圓頻率， 表達式為:

(13)

(14)

ω=(λH)2ω2+ω0

(15)

海床表面處波致孔隙壓力為波流共同作用產(chǎn)生的波壓力，海床表面為排水界面，累積孔隙水壓力為0，并且平均有效應(yīng)力和剪應(yīng)力消失：

(16)

海床兩側(cè)為周期性邊界條件，上下兩層海床交界面處為連續(xù)邊界條件，即：

pres左=pres右;pres上=pres下

(17)

海床底部為不可壓縮不可排水邊界，即

(18)

1.3 數(shù)值模擬驗證

通過Sumer et al.(2012)的試驗數(shù)據(jù)來驗證本文孔壓數(shù)值模型的有效性，將試驗的土體和波浪參數(shù)輸入到本文模型中，計算結(jié)果與試驗數(shù)據(jù)對比如圖4所示，可以看出本文孔壓數(shù)值模型計算結(jié)果與Sumer et al.(2012)的試驗數(shù)據(jù)基本吻合，表明了本模型的可靠性。

圖 4 累積孔隙水壓力隨時間的變化(z=-0.085 m)Fig. 4 Variation of residual pure water pressure with time

2 累積孔隙水壓力響應(yīng)分析

基于建立的數(shù)值模型，根據(jù)試驗以及相關(guān)資料獲得了黃河三角洲地區(qū)原始海床以及表面硬殼層的物理力學(xué)指標(biāo)，波浪條件為黃河三角洲地區(qū)一年一遇海況，并選取不同的海流流速，研究波浪共同作用下不同流速對累積孔隙水壓力分布的影響。波流參數(shù)與海床土體參數(shù)見表 1。

表 1 波流參數(shù)與海床土體參數(shù)取值Table 1 Wave-current parameters and seabed soil parameters

2.1 原始單層海床累積孔壓分布特征

考慮未形成表面硬殼層的原始單層海床，觀察不同波浪作用時間累積孔隙水壓力的變化以及液化情況。

由圖 5 表明隨著波浪作用時間的持續(xù)增長，海床累積孔隙水壓力逐漸增加，但累積孔壓每隔30 min的增量隨著時間而減小，這是由于孔壓的消散以及附加應(yīng)力導(dǎo)致的海床孔隙壓縮極限無法再產(chǎn)生塑性變形，最終穩(wěn)定后的孔隙水壓力稱為穩(wěn)態(tài)超孔壓力。由于真實海況波浪持續(xù)時間以及數(shù)值模型計算量的原因，本文只考慮波浪作用120 min內(nèi)的累積孔隙水壓力結(jié)果。由圖 5 可知，累積孔壓沿深度分布存在極大值，該極值隨波浪持續(xù)作用時間逐漸下移，并且由圖 5 還可以看出累積孔隙水壓力沿深度分布曲線與平均有效應(yīng)力線的交點隨時間逐漸下移，海床液化深度范圍逐漸增大，海床的液化深度在波浪持續(xù)作用90 min的時間間隔內(nèi)，液化深度由0.62 m增加至2.14 m。

圖 5 波流作用下累積孔隙水壓力沿深度分布隨時間的變化Fig. 5 Variation of residual pore water distribution along depth with time under wave and current action

海流與波浪聯(lián)合作用會顯著影響海床累積孔隙水壓力的分布，取表中的不同海流流速參數(shù)，分析不同流速對海床累積孔壓分布特征的影響規(guī)律。圖 6為波流持續(xù)作用120 min時，不同流速下海床孔隙水壓力沿深度的分布圖。

由圖 6 可知，隨著流速的增大，海床累積孔壓相應(yīng)增大，流速由0 m · s-1增大至2 m · s-1時，海床累積孔隙水壓力極值由12.15 kPa增加到29.63 kPa。但流速為負值，即流速流向與波浪行進方向相反時，海床孔隙水壓力則相應(yīng)減小，流速為- 2 m · s-1時的累積孔壓與流速為0 m · s-1相比僅下降了約40%。所以，當(dāng)流速流向與波浪行進方向一致時，海流會對海床累積孔隙水壓力的增長起促進作用，顯著增加了海床累積孔隙水壓力，而流向相反時，會對海床累積孔隙水壓力起到抑制作用，這種抑制作用相對于流速一致時累積孔壓的變化，影響較小。

圖 6 波流聯(lián)合作用不同流速下累積孔隙水壓力沿深度分布Fig. 6 Vertical distribution of residual pore water pressure for different velocity under wave and current action

2.2 硬殼層厚度對海床累積孔壓的影響

黃河三角洲地區(qū)分布著厚度不一的硬殼層，選取表 1 中的硬殼層相關(guān)物理力學(xué)指標(biāo)，并分別取硬殼層厚度為1 m， 2 m， 3 m，研究不同硬殼層厚度對海床累積孔壓的影響規(guī)律。圖 7為波流持續(xù)作用120 min時，不同硬殼層厚度下海床累積孔隙水壓力沿深度的分布圖。

圖 7 波流聯(lián)合作用不同硬殼層厚度對累積孔隙水壓力的影響(t=120 min)Fig. 7 Effect of hard shell thickness on cumulative pore water pressure under wave and current actiona. 硬殼層厚度1 m； b. 硬殼層厚度2 m； c. 硬殼層厚度3 m

由圖 7可以看出硬殼層存在時，海床累積孔隙水壓力顯著減小，隨著硬殼層厚度的增加而相應(yīng)減小，累積孔隙水壓力沿深度分布出現(xiàn)轉(zhuǎn)折，在原始海床中沿深度增加迅速。這表明密度大、孔隙大的表面硬殼層有利于海床累積孔壓的消散。海床累積孔隙水壓力沿深度分布的極值均出現(xiàn)在下層原始單層海床中，流速為0 m · s-1時硬殼層厚度由1 m增加到3 m，極值點深度下降了1.38 m。

由圖 7還可以看出累積孔壓沿深度分布極值所處深度均靠近兩層海床界面處，硬殼層厚度越小越接近，而且該極值大小隨著厚度的增加也約接近海床界面處累積孔隙水壓力值。這是由于波流作用于海床影響顯著的區(qū)域大約在1～3 m之間(圖 6)，然而硬殼層的存在有利于海床孔隙水壓力的消散，因此隨著硬殼層厚度的增加，極值的增量也相應(yīng)減小。

3 滲流模式下硬殼層形成機制

波浪及潮流共同作用于海床表面，會使海床表層附近的細粒泥沙顆粒發(fā)生啟動，造成表層顆粒較粗、密度較大的情況。對于厚度較大的硬殼層是由形成于海床內(nèi)部的滲流力導(dǎo)致的，由上文可知在波浪以及潮流作用下，海床內(nèi)部會產(chǎn)生較大的累積孔隙水壓力，這些累積孔隙水壓力會逐漸向外排出，海床相鄰點處的累積孔隙水壓力差值形成的水力梯度產(chǎn)生滲流力，進而引起水流動，因為海床表面為排水界面，從而形成向上的滲流。海床內(nèi)部的細顆粒泥沙受到滲流力的作用打破受力平衡狀態(tài)產(chǎn)生位移，繼而沿著粗顆粒之間的孔隙運移至海床表層。海床內(nèi)部累積孔隙水壓力隨著時間逐漸增大，因此在海床深度較大的位置也可產(chǎn)生指向床面的滲流力，形成較厚的粗顆粒層。累積孔隙水壓力的形成是由于波浪作用于海床，在海床內(nèi)部產(chǎn)生循環(huán)附加應(yīng)力，使海床內(nèi)部顆粒之間的孔隙壓縮產(chǎn)生超孔隙水壓力，來不及排除以至逐漸累積，壓縮孔隙體積使海床變得致密，細顆粒向上運移和孔隙壓縮兩方面同時進行形成了硬殼層海床。

根據(jù)王虎等(2014)的研究， 當(dāng)波流引起的海床表面最大切應(yīng)力超過臨界啟動切應(yīng)力時，土體顆粒啟動，進而形成沖刷。很多學(xué)者在傳統(tǒng)泥沙顆粒力學(xué)平衡方程中考慮了滲流力的作用，并建立了新的力學(xué)平衡方程，并且竇國仁等(2001)給出了該力學(xué)平衡方程中的重力FG、水流拖拽力FD、上舉力FL以及滲流力FS的表達式。根據(jù)該新的力學(xué)平衡方程以及各個受力的表達式，王虎等(2014)推導(dǎo)建立了海床臨界沖刷深度計算公式，研究滲流力對海床沖刷的影響。推導(dǎo)過程如下：

泥沙顆粒受力平衡方程(圖 2) 表達式為：

FD=(FG-FL-FS)tanφ

(19)

式中：FD為水流拖拽力；φ為飽和土靜止內(nèi)摩擦角。其中，F(xiàn)G，F(xiàn)L，F(xiàn)D，F(xiàn)S表達式為(竇國仁等， 2001)：

(20)

(21)

(22)

(23)

式中：τs為臨界啟動切應(yīng)力；d50為中值粒徑； 阻力系數(shù)CD=0.4； 上舉力系數(shù)CL=0.1； ΔP為相距ΔL兩點間的累積孔隙水壓力差值。

根據(jù)式(17)～式(21)，臨界啟動切應(yīng)力τs表達式為：

(24)

通過數(shù)值計算模型求得的海床土體累積孔隙水壓力值，進而可得到海床任意深度處的滲流壓力梯度，然后根據(jù)式(22)求解該深度處土體顆粒的臨界啟動切應(yīng)力。

波流作用引起的海床表面最大切應(yīng)力(錢寧等， 2003)為：

(25)

式中：ρ1為床面附近流體密度；um為海底最大軌跡質(zhì)點速度；fW為瓊森阻力系數(shù)，決定于水流雷諾數(shù)Re和海床相對粗糙度2am/d50，兩參數(shù)取值表達式如下：

(26)

(27)

取得Re和2am/d50值，根據(jù)Sassa et al.(2001)文獻可得fW值，最終，可由式(25)求得海床表面最大切應(yīng)力值。

床面最大切應(yīng)力τm大于臨界啟動切應(yīng)力τs時，泥沙可發(fā)生啟動。細顆粒泥沙首先發(fā)生啟動，通過粗顆粒之間形成的孔隙運移至海床表面。當(dāng)τm小于τs時，則不會發(fā)生沖刷。相關(guān)計算結(jié)果見表 2。

表 2 相關(guān)計算參數(shù)及結(jié)果Table 2 Relevant calculation parameters and results

由表 2 可知，當(dāng)流速U0=1 m·s-1時海床深度為0 m處的臨界啟動切應(yīng)力為負值，流速U0=2 m · s-1時臨界啟動切應(yīng)力負值已達到海床深度0.5 m處，這說明該兩種流速下引起的ΔP/ΔL值較大，以致(ρ-ρW)g-ΔP/ΔL值為負，τs為負值。流速流向與波浪行進方向一致時可形成較大的滲流力作用于泥沙顆粒上，滲流力越大可以將較粗的泥沙顆粒運移至海床表面，形成更粗的硬殼層。并且由表 2 還可以看出隨著流速的增大，泥沙顆粒的臨界啟動切應(yīng)力也相應(yīng)減少，這表明流速越大，泥沙顆粒發(fā)生啟動所需要的海床表面最大切應(yīng)力也越小，泥沙越容易啟動。

圖 8 不同流速下切應(yīng)力沿深度變化Fig. 8 Shear stress changes with depth at different velocitya. 流速U0=0 m·s-1； b. 流速U0=1 m·s-1

由圖 8 可知，τs隨著海床深度的增加急劇下降，τm隨著深度增加而增加，兩值隨深度快速接近，兩條線相交處的地方泥沙啟動停止。從圖 8 也可以看出兩條線并沒有相交，但是這兩個不同流速下的泥沙啟動在1.5 m深度處停止，這是由于波浪引起的累積孔隙水壓力沿深度分布極值在1.5 m深度附近，并且流速存在對于該極值的深度分布貢獻也極小(圖 6)，所以海床泥沙啟動只發(fā)展到1.5 m深度左右。該深度符合黃河三角洲地區(qū)硬殼層厚度普遍為1～3 m的情況。黃河三角洲地區(qū)在遭受5年一遇以及50年一遇的復(fù)雜海況條件時，累積孔隙水壓力沿深度分布極值會分布在海床更深的位置，并且引起較大的ΔP/ΔL值，從而形成較厚以及顆粒較粗的表面硬殼層。

4 結(jié) 論

(1)波浪引起的累積孔隙水壓力與波浪持續(xù)時間成正相關(guān)，累積孔隙水壓力隨時間逐漸增加，但增速逐漸放緩，最終會達到穩(wěn)定狀態(tài)，海床孔隙不可再壓縮。

(2)海流對海床累積孔隙水壓力有重要影響。海流流向與波浪行進方向一致時，對累積孔隙水壓力起促進作用，流速越大累積孔隙水壓力越大，反之對累積孔隙水壓力有抑制作用。

(3)表面硬殼層的存在會顯著促進累積孔壓的消散，累積孔隙水壓力沿深度分布的極值均出現(xiàn)在下層原始海床中，流速U0=0 m · s-1時硬殼層厚度由1 m增加到3 m，極值點深度下降了1.38 m。

(4)累積孔隙水壓力引起的滲流力對于海床泥沙啟動影響顯著，在流速U0=0 m · s-1，U0=1 m · s-1時泥沙啟動深度均為海床1.5 m深度處，并且海流流向與波浪行進方向一致時，會產(chǎn)生較大ΔP/ΔL值帶動較粗的泥沙顆粒至海床表層，但對泥沙啟動的最大深度影響不大。