劉成林，陳煒昀，陳國興，俞 縉

(1. 河海大學(xué) 港口海岸與近海工程學(xué)院，江蘇 南京 210098; 2. 南京工業(yè)大學(xué) 巖土工程研究所，江蘇 南京 210009; 3. 華僑大學(xué) 土木工程學(xué)院，福建 廈門 361021)

近海工程埋置型結(jié)構(gòu)物(如沉管隧道、海底管道等)在建設(shè)過程中，通常需要預(yù)先在海床表面開挖一段基槽，便于建筑物埋置。這種處理方式使得建筑物能減少周圍環(huán)境作用力，也避免了在不穩(wěn)定的新近沉積土海床表面施工。在自然海洋環(huán)境中，海洋波浪時常與海流(比如洋流、渦流、潮汐流等)同時存在，這種波流共同作用下的荷載將進一步使海床中的孔壓和有效應(yīng)力產(chǎn)生變化[1]。在一定條件下，波流荷載下的超孔隙水壓力超過土體自身的初始有效應(yīng)力，相應(yīng)土體會發(fā)生失穩(wěn)并產(chǎn)生液化[2]。同時，由于海床基槽的開挖，原有海床土體在長期自然環(huán)境影響下形成的初始應(yīng)力和固結(jié)狀態(tài)會發(fā)生改變。因此，波流作用下的海床響應(yīng)和液化評估應(yīng)該基于工程擾動后的海床狀態(tài)。

通常情況下，動荷載作用下的海床響應(yīng)機理可以根據(jù)孔隙水壓力的產(chǎn)生機理分為兩種形式。其一是由于波浪在海床中產(chǎn)生孔隙水壓力的相位滯后性而導(dǎo)致的瞬態(tài)液化，通常發(fā)生在波谷之下較為密實的海床中[3]。另一種是較為松散海床在波浪循環(huán)荷載下密實，土顆粒進行重新分布，殘余孔壓累積最終可能誘發(fā)液化[4]。關(guān)于兩種液化發(fā)生形式和適用范圍的討論在文獻[5]中有詳細論述。然而，不論何種液化發(fā)生，對建設(shè)在海床上的近海工程結(jié)構(gòu)物穩(wěn)定性都會構(gòu)成很大威脅，目前已有大量研究關(guān)注于防波堤、管線周圍海床在波浪荷載下的動態(tài)響應(yīng)[6-8]。由于文中主要研究的是基槽開挖后的臨時海床形態(tài)受波流荷載下的動態(tài)響應(yīng)，因此著重對海床的瞬態(tài)響應(yīng)及瞬態(tài)液化情況進行討論。

眾所周知，海床的失穩(wěn)破壞是在一定的靜荷載和動荷載共同作用下發(fā)生的。其中，靜荷載指海床土體自身的重力，取決于海床的幾何形狀和地質(zhì)情況。動荷載則來自于波流作用對床面產(chǎn)生的動壓力，與海床所處波流環(huán)境有關(guān)。近年來開展了很多關(guān)于海床動態(tài)響應(yīng)的研究，包括解析解形式[5]、試驗研究形式[9-10]以及數(shù)值模擬研究形式[11-12]。最近，Lin等[13]采用水平集方法[14-15]捕捉兩相流中水和空氣的表面，基于有限元模型分析了波浪荷載下幾種特定管道埋置形式下的動態(tài)海床響應(yīng)。年廷凱等[16]基于極限分析上限方法，求解了極端條件下不同時刻的斜坡海床穩(wěn)定性系數(shù)。Duan等[17]采用有限元方法實現(xiàn)了波流-海床-管道耦合計算，計算了在基槽內(nèi)部分回填管道周圍無黏性土的液化情況，得到了管道回填厚度選取與波浪參數(shù)的經(jīng)驗公式。Liao等[18]采用兩階段耦合數(shù)值框架計算了傾斜防波堤端部周圍土體的動態(tài)響應(yīng)。然而，以上研究中大部分都是對未經(jīng)擾動的平坦海床進行動態(tài)響應(yīng)計算，對海底人工開挖基槽形成的不規(guī)則海床在波流共同作用下的響應(yīng)少有研究。

為了研究開挖之后的基槽周圍海床在波流荷載作用下的動態(tài)響應(yīng)，建立了一個二維波流-海床耦合模型。通過雷諾時均納維—斯托克斯(RANS)方程結(jié)合湍流模型模擬波浪場，在模型兩側(cè)邊界設(shè)定穩(wěn)定流速邊界實現(xiàn)流場的模擬，海床的動態(tài)響應(yīng)通過Biot多孔彈性理論計算得出。分析了基槽開挖后的海床初始應(yīng)力和應(yīng)變的變化情況，討論了波浪參數(shù)及海流參數(shù)對海床響應(yīng)結(jié)果的影響(包括孔壓變化和液化情況變化)，以及波流耦合下的波浪形態(tài)變化。最后，提出了流速與海床液化深度的經(jīng)驗關(guān)系，對于近海水下基槽開挖的設(shè)計和施工提供了建議和參考。

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

考慮波流共同作用下開挖后的基槽海床動態(tài)響應(yīng)問題，計算模型如圖1所示。模型主要分為兩部分，波流子模型和海床子模型。波浪部分通過求解RANS方程和k-ε湍流閉合方程實現(xiàn)波浪場模擬，其中通過求解包含動量源項的方程實現(xiàn)動量源造波。海床部分通過Biot理論得到波-流底部壓力作用下的孔壓和位移。在波浪傳遞區(qū)域的上邊界定義定常流入流速，下邊界設(shè)定定常流出流速和壓力輸出邊界實現(xiàn)穩(wěn)定的流場模擬，從而實現(xiàn)波流同向相互作用(如圖1所示)。反之，則可以實現(xiàn)波流反向相互作用。

圖1 模型計算示意Fig. 1 Sketch of the analysis model

1.1 波流模型

1.1.1 波浪模型

模型采用動量源函數(shù)造波，通過將動量源函數(shù)代入動量守恒公式中的源項得以實現(xiàn)。另外，通過水平集(LSM)方法結(jié)合移動網(wǎng)格方法捕捉空氣和波浪兩相流體的接觸表面，即隨時間而改變的波浪自由表面。

模型中采用求解實際應(yīng)用中最為常用的標(biāo)準(zhǔn)k-ε湍流模型結(jié)合RANS方程來實現(xiàn)波浪場的模擬。眾所周知的是，湍流是一種復(fù)雜、非穩(wěn)態(tài)、旋轉(zhuǎn)的不規(guī)則流動形式，由黏性力引起。k-ε湍流模型主要包括兩個輸運方程和兩個獨立變量：湍流動能k和湍流耗散率ε。湍流黏度表示為：

(1)

其中，ρ為流體密度，Cμ表示一個湍流模型參數(shù)。

湍流動能k的輸運方程可表示為：

(2)

耗散率ε的輸運方程表示為：

(3)

上述方程(1)～(3)中的模型常數(shù)Cμ，Cε1，Cε2，σk，σε分別為0.09，1.44，1.92，1.0，1.3[19]。方程(2)～(3)中Pk表示湍流壓力分布項。

由于模型考慮的為二維問題，故質(zhì)量守恒方程和動量守恒方程可以表示為：

(4)

(5)

(6)

(7)

其中，k為波數(shù)；ω為波浪角頻率；β=80/(δ2·L2)為造波源區(qū)的寬度，其中L為波長并且δ為表征造波寬度的一個參數(shù)；Ds為源區(qū)的振動幅度。更多動量源造波的介紹可參考Wei等[20]，Choi和Yoon[21]的研究。

模型中采用水平集方法(LSM)捕捉空氣和波浪兩相流的接觸面，即波浪表面。相關(guān)公式可以表示為：

(8)

其中，φ為水平集函數(shù)；εt是控制兩相流界面厚度的參數(shù)，其默認值為跨越兩相界面處網(wǎng)格最大尺寸的一半；γ是重新初始化的參數(shù)；u為流速場的速度分量，采用的流速通過對流的方法實現(xiàn)水平集參數(shù)的讀取。關(guān)于水平集方法在Desjardins等[22]的研究中有詳細介紹。

1.1.2 海流模型

在數(shù)值波浪傳遞區(qū)域的上側(cè)邊界設(shè)定穩(wěn)定的流速入口，在下側(cè)邊界設(shè)定穩(wěn)定的流速出口，可以實現(xiàn)波流同向傳遞。反之則可以實現(xiàn)波流反向傳遞。法向流速設(shè)置為u=-nU0，其中n表示流入/流出邊界的法向向量，U0為穩(wěn)定海流的速度。

1.2 海床模型

由于模型主要研究的工況為基槽開挖完畢的臨時工況，波浪作用于海床斷面的時間較短，因此文中主要關(guān)注于土體的瞬態(tài)響應(yīng)。通常，沉管隧道的天然海床基礎(chǔ)可視為多相介質(zhì)材料，包括固體相(土骨架)、流體相(水和氣體)。二維海床模型的建立采用Biot固結(jié)方程(“Q-S方程”)，由Biot提出的準(zhǔn)靜態(tài)方程描述土體，模型忽略孔隙流體和土顆粒的慣性加速項，同時不考慮孔隙流體對土顆粒的相對變形[23]。

基于準(zhǔn)靜態(tài)方程，將海床考慮為各向同性，各方向上滲透系數(shù)相同?？紫读黧w的質(zhì)量守恒方程可以表示為：

(9)

(10)

其中，Kf表示孔隙流體的實測體積模量；Kw為水的真體積模量；p0為絕對靜水壓力，在模型中取為p0=γwd，d是水深；Sr為土體飽和度。

孔隙水壓的控制方程可表述為：

(11)

(12)

式(11)和(12)中的有效應(yīng)力可以從多孔彈性理論得出：

(13)

(14)

(15)

其中，Gs表示土體剪切模量，通過楊氏模量Es和泊松比υs可以表示為Gs=Es/[2(1+υs)]。

1.3 邊界條件

數(shù)值模型通過波流模塊模擬波浪場，隨著波浪的傳播，波壓力傳遞到有限厚度的多孔介質(zhì)海床上，有限厚度的海床底部視為剛性不可滲透。靜止水面在z=d+h的位置，波浪傳播沿著x軸正方向，豎向z軸以海床底部開始向上為正(如圖1)。波浪場中波高為H，波長為L，靜止水面距離海床表面的距離為d，有限海床的厚度為h。通過式(11)～(15)可以求解出海床土體中瞬態(tài)孔壓和土體的變形情況?？刂品匠淌降那蠼庑枰Y(jié)合相應(yīng)的邊界條件，包含海床各界面處邊界條件、自由水面的邊界條件。

在海床表面，垂直有效應(yīng)力和剪應(yīng)力變?yōu)?，孔隙水壓力等同于海床表面處波浪壓力。

(16)

其中，pb為波浪場在海床表面的波浪壓力，ps為孔隙水壓力。

海床底部視為不可滲透基巖，海床底部沒有位移和垂向滲流，表示為：

(17)

海床兩側(cè)邊界視為不可滲透邊界，同時沒有水平向位移，表示為：

(18)

在波浪自由表面，由于兩相流模型中的壓力連續(xù)條件，波浪自由表面的壓力等于大氣壓。

1.4 模型耦合

模型耦合流程首先要模擬海床基槽周圍的流場與壓力場，再將壓力作為邊界條件施加于海床之上以實現(xiàn)海床的孔壓和位移計算。最后，將兩部分子模塊在模型中實現(xiàn)物理場耦合，從而在單位時間步長內(nèi)實現(xiàn)波流場壓力下的土體響應(yīng)單向耦合計算。

模型的網(wǎng)格主要分為固定網(wǎng)格和移動網(wǎng)格兩種，在波流場自動剖分為三角形網(wǎng)格，同時在空氣和波浪表面附近采用自動細分方法實現(xiàn)網(wǎng)格較好地與物理場匹配。在開挖基槽部分進行網(wǎng)格細化，以實現(xiàn)更為細致地對此部分流場計算，使計算結(jié)果更可靠。移動網(wǎng)格則被用于捕捉自由水面和空氣的接觸面，自由變形的移動網(wǎng)格采用拉格朗日平滑類型。

1.5 模型驗證

首先，將模型計算結(jié)果與前人的試驗值或解析解進行比較。主要的驗證內(nèi)容分為兩部分，波流場驗證和海床動態(tài)響應(yīng)驗證。為了驗證造出波浪的準(zhǔn)確性，保證海床受到穩(wěn)定波浪荷載，將波浪部分的自由波面形態(tài)和相應(yīng)計算區(qū)域特征點處水質(zhì)點運動與解析公式解進行對比。波浪采用Liu 等[25]研究中的波浪參數(shù)，波浪參數(shù)為波高H=0.143 m，周期T=1.75 s，水深d=0.533 m。選取了造波源區(qū)之外1 m位置處的波面振動情況與解析解進行對比，如圖2所示，模型計算結(jié)果基本上和解析結(jié)果相吻合。同時在波浪部分提取波浪自由表面形態(tài)與解析解對比驗證，結(jié)果如圖3所示。模型計算結(jié)果與解析公式結(jié)果基本趨于一致，僅有起始位置處存在少許相位偏差，可能是造波初始階段未達到穩(wěn)定導(dǎo)致。整體而言，波浪部分的數(shù)值波浪結(jié)果和預(yù)期解析公式保持一致。

圖2 特征點處波面振動情況驗證Fig. 2 Validation of the free wave surface vibration at characteristic point

圖3 波浪自由表面驗證Fig. 3 Validation of the free wave surface

為了驗證波流相互作用下波面水質(zhì)點振動情況，將模型與Chen等[26]試驗結(jié)果進行對比，模擬波流同向(U=0.2 m/s)和波流反向(U=-0.2 m/s)的特征水質(zhì)點處波面振動情況。Chen等[26]所采用波浪參數(shù)為：波高H=0.05 m，水深d=0.6 m，波周期T=1 s。將模型結(jié)果與Chen等[26]試驗中G3波浪儀測量數(shù)據(jù)對比，結(jié)果如圖4所示，結(jié)果顯示模型可以有效反映波流共同作用下的波浪特征。

圖4 波流相互作用下波面驗證Fig. 4 Wave surface validation under wave-current interaction

在海床土體動態(tài)響應(yīng)驗證部分，將模型結(jié)果與Qi和Gao[27]試驗數(shù)據(jù)進行對比，試驗的主要參數(shù)為：波高H=0.102 m，波周期T=1.2 s，流速U=0.23 m/s；土體滲透系數(shù)K=1.88×10-4m/s，海床相對密實度Dr=0.352，土體浮重度γ′=9.03 kN/m3。試驗中裝土的特制土箱固定在水槽中段，采用孔壓計測量土體表面和下部孔壓的變化，具體試驗布置可在Qi和Gao[27]研究中得到。模型結(jié)果與試驗中的孔壓計PPT1和PPT2對比，結(jié)果如圖5所示，展示了兩個波浪周期內(nèi)的土層表面和表面之下的孔壓變化。數(shù)值結(jié)果與試驗數(shù)據(jù)具有良好一致性，即模型對波流壓力下的海床孔壓模擬具有良好應(yīng)用性。

圖5 數(shù)值模型與Qi和Gao試驗數(shù)據(jù)對比結(jié)果Fig. 5 Comparison between numerical and experimental data of Qi and Gao

2 數(shù)值模型應(yīng)用

2.1 海床初始固結(jié)

在工程活動如水下基槽開挖影響下，海床土體會進入新的固結(jié)狀態(tài)?？紫端畨毫ο?，有效應(yīng)力增加，基礎(chǔ)會產(chǎn)生一定沉降[28]。首先確定受擾動之后海床的初始固結(jié)狀態(tài)和有效應(yīng)力分布情況。文中選用的波浪、土體參數(shù)和邊坡參數(shù)見表1。波浪和土體參數(shù)參考渤海區(qū)域極端情況下記錄的參數(shù)[29]，計算沉管隧道開挖形成的基槽在波流共同作用下的動態(tài)響應(yīng)。后文無特別說明情況下均采用表中參數(shù)?；坶_挖后的有效應(yīng)力分布情況如圖6所示，作為初始有效應(yīng)力應(yīng)用在后續(xù)分析中。

表1 研究中的參數(shù)取值Tab. 1 Parameters in this analysis

圖6 基槽開挖完成后的初始有效應(yīng)力分布Fig. 6 Initial effective stress distribution after trench excavation

2.2 瞬態(tài)液化

波浪作用下的海床液化對于近海工程結(jié)構(gòu)物的設(shè)計施工影響重大。目前已有多種波浪作用下的海床液化判別依據(jù)[30-32]，文中采用Zen和Yamazaki[30]提出的二維液化判別準(zhǔn)則：

(19)

圖7給出了計算沉管隧道基槽開挖之后海床土體液化深度為最大值的情況，上半部分為波浪作用于海床表面的動壓力，下半部分深色區(qū)域展示液化區(qū)域。模型采用了表1的參數(shù)，計算結(jié)果和前人得出的結(jié)論保持一致，瞬態(tài)液化主要發(fā)生在海床表層區(qū)域，由于孔壓產(chǎn)生的相位滯后性而周期性出現(xiàn)在波谷之下[33-34]。液化區(qū)域分布和最大液化深度的計算對海底人工基槽設(shè)計和施工具有一定指導(dǎo)意義。

圖7 基槽開挖完成后海床土體液化情況Fig. 7 Seabed soil liquefaction after trench excavation

2.3 波浪參數(shù)的影響

對波浪參數(shù)對海床的液化情況和孔隙水壓力分布影響進行參數(shù)分析。分別取不同波浪高度H=3 m，4 m，5 m，采用造波模型產(chǎn)生相應(yīng)的波浪計算基槽周圍土體的最大液化深度，結(jié)果如圖8(a)所示，隨著波高的增大，基槽部分的最大液化深度分別為0.98 m，1.30 m，1.63 m。取不同的波浪周期T= 6 s，7 s，8 s，分析波浪周期對最大液化深度的影響。從圖8(b)可以看出，最大液化深度隨著波浪周期的增大呈指數(shù)形式增大，海床最大液化深度分別為1.35 m，1.52 m，2.10 m。即長周期的波浪對海床土體液化破壞顯著增大。

圖8 不同波高和波周期條件下海床土體最大液化深度變化Fig. 8 Variation of maximum liquefaction depth of seabed soil under different wave heights and wave periods

同時，計算得出基槽內(nèi)最大孔隙水壓力隨不同波高和不同波周期的垂向分布情況。結(jié)果如圖9(a)所示，最大孔壓|ps|max/0.5γwH隨著波高的增大而增大。但3種不同波高情況下當(dāng)z/hs=-1時(即海床底部)，孔隙水壓力趨近于0。如圖9(b)所示，隨著波浪周期的增大，最大孔壓也相應(yīng)增大。這些規(guī)律和圖8中液化深度分析呈現(xiàn)的規(guī)律相一致，正是因為海床中的超孔隙水壓力增大超過海床自身的初始有效自重，導(dǎo)致了液化的發(fā)生。

圖9 海床土體的孔壓分布隨著波浪高度和周期的變化Fig. 9 Pore pressure distribution of seabed soil varies with wave height and period

2.4 海流參數(shù)的影響

如前所述，波流共同作用下的波浪形態(tài)和僅存在波浪時有很大不同。針對不同流速(范圍-2～2 m/s，以0.5 m/s為間隔)對波浪形態(tài)影響進行參數(shù)分析，旨在得到波-流共同作用對波浪主要特征參數(shù)的影響。結(jié)果如圖10所示，對計算得到的波高和波長采用無量綱化處理，即同輸入波高H和波長L(見表1)的比值。波浪在同向海流作用下波長增大波高減小，波陡度減小。相反，在反向海流的作用下，波長變短而波高增大，波陡度增大。

圖10 波浪在不同流速下的波高和波長變化Fig. 10 Change of wave height and wave length with different current velocities

圖11 不同流速波流共同作用下海床最大孔壓分布情況Fig. 11 Distribution of maximum pore pressure in the seabed with different current velocities

圖12 不同流速下最大液化深度 Fig. 12 Maximum liquefied depth with different current velocities

海流對海床動態(tài)響應(yīng)和破壞情況影響重大，在實際海洋環(huán)境中，波浪和海流(潮流)常常同時存在。在不同的波流組合方式作用下，海床中的孔壓和有效應(yīng)力分布會產(chǎn)生顯著變化。分析不同海流流速與波浪的公共作用對海床中孔壓分布的影響，對海床受波流共同作用影響的分析有重要意義。選取三種不同的海流速度，-1 m/s，0 m/s，1 m/s，分別計算波流反向作用、僅波浪作用和波流同向作用下的海床動態(tài)響應(yīng)結(jié)果。如圖11所示，海床不同深度處的最大孔壓在波流反向的情況下達到最大，而在波流同向的情況下最小。

在上述分析的基礎(chǔ)上，模型計算了海床最大液化深度在相同波浪參數(shù)(見表1)伴隨不同流速下的變化，流速變化范圍為-2～2 m/s，以0.5 m/s為間隔。結(jié)果如圖12所示，在反向海流作用下，海床的最大液化深度同比有所增大。而在同向海流作用下，海床的最大液化深度相應(yīng)減小。值得注意的是，反向海流對液化深度的影響(1.65 m相對于1.32 m，同比增大25%)比同向海流的影響(1.12 m相對于1.32 m，同比減小15%)更大。圖12中基于線性回歸分析給出了不同流速下最大液化深度ZL的擬合曲線為：

ZL=0.013 6u03+0.012 1u02-0.186 1u0+1.347 5

(20)

式中：u0為流速，擬合曲線的相關(guān)參數(shù)R2=0.990 9。

3 結(jié) 語

研究波流荷載下開挖基槽周圍海床的動態(tài)響應(yīng)，提出基于RANS方程和Biot固結(jié)“Q-S”理論的波浪-海床耦合計算數(shù)值模型。將模型詳細地與前人試驗解和解析解對比驗證后，應(yīng)用于計算波流共同作用下的波浪形態(tài)變化，以及基槽海床動態(tài)響應(yīng)計算，包括孔壓分布和液化深度變化?；跀?shù)值計算結(jié)果，結(jié)論可以總結(jié)如下：

1) 基槽開挖之后的初始應(yīng)力場和應(yīng)變場會有很大改變，后續(xù)動態(tài)響應(yīng)分析應(yīng)當(dāng)以開挖后新的固結(jié)狀態(tài)作為初始狀態(tài)考慮。

2) 開挖之后基槽內(nèi)土體的最大孔壓|ps|max/0.5γwH和液化深度ZL隨著波浪周期T和波高H的增大而增大。

3) 波浪在同向海流作用下波長增大，波高減小，波陡度減小；在反向海流作用下波長減小，波高增大，波陡度增大。同樣條件下，反向海流對波浪形態(tài)的影響更大。

4) 反向海流作用下，海床內(nèi)的最大孔壓和液化深度都相應(yīng)增大。同時，反向海流對最大液化深度的影響更大?；诰唧w波浪參數(shù)提出的流速與最大液化深度關(guān)系對相應(yīng)研究和工程應(yīng)用具有參考價值。