潘佳禾， 廖晨聰,b,c， 陳錦劍,b,c

(上海交通大學(xué) a. 土木工程系； b. 海洋工程國(guó)家重點(diǎn)實(shí)驗(yàn)室； c. 高新船舶與深海開(kāi)發(fā)裝備協(xié)同創(chuàng)新中心， 上海 200240)

隨著人類對(duì)海洋的開(kāi)發(fā)，近岸海域內(nèi)工程結(jié)構(gòu)物的數(shù)量日益增加.以斜坡為主的近岸海床是海底油氣管道“登陸”的重要過(guò)渡地帶.在近岸淺水環(huán)境下，海床及管道的受力狀態(tài)易受波浪作用的影響，淺水區(qū)域的波浪呈現(xiàn)出諸如孤立波的特殊形態(tài).孤立波在斜坡海床上爬升和回落時(shí)，會(huì)引起孔隙水壓力變化，可能導(dǎo)致海床的液化與滑坡，進(jìn)而危及海底管道等近岸結(jié)構(gòu)物.

孤立波是一種模擬近岸淺水環(huán)境的常用波形，也被廣泛地用于模擬海嘯等極端波浪[1].當(dāng)孤立波傳播到近岸斜坡海底地帶時(shí)，由于水深變淺，波浪的波幅和非線性均有所增加，而到達(dá)海岸線后波浪還會(huì)繼續(xù)向陸地爬升及回落.對(duì)于孤立波在斜坡上的破碎、爬升及回落過(guò)程已有較多實(shí)驗(yàn)觀測(cè)和數(shù)值模擬研究.Synolakis[2]通過(guò)水槽實(shí)驗(yàn)捕捉斜坡上孤立波的變化波面.宣瑞韜[3]采用考慮造波區(qū)波動(dòng)非穩(wěn)態(tài)特性的造波方法，在水槽實(shí)驗(yàn)中模擬孤立波和雙孤立波波形的傳播及爬高過(guò)程.王賀等[4]基于雷諾平均Navier-Stokes(RANS)方程、流體體積(VOF)方法及修正的Goring造波方法，實(shí)現(xiàn)對(duì)孤立波在斜坡上的破碎、爬升及回落過(guò)程的數(shù)值模擬.實(shí)際上，孤立波在近岸斜坡上的這種破碎、爬升及回落過(guò)程將會(huì)影響海床受力及孔壓響應(yīng).

關(guān)于波浪、海床與結(jié)構(gòu)物之間的相互作用問(wèn)題，已有不少學(xué)者開(kāi)展了研究工作.劉博等[5]給出一個(gè)近似u-p(固相位移-孔壓)動(dòng)力模型的解析解，考慮在波流共同作用下多孔介質(zhì)海床的孔隙水壓力及有效應(yīng)力的變化，并研究了潮流流向?qū)４岔憫?yīng)的影響.Jeng等[6]通過(guò)Navier-Stokes(N-S)方程控制波浪，提出一個(gè)2維波浪-海床-防波堤之間相互作用的單向耦合模型，較精準(zhǔn)地模擬了三者之間的相互作用.Zhang等[7]在此基礎(chǔ)上進(jìn)一步提出單樁-波浪-海床之間相互作用的3維模型，研究海床土體的動(dòng)力響應(yīng)問(wèn)題.張軍等[8]利用COMSOL軟件構(gòu)建了波浪-海床-管線計(jì)算模型，模擬了在1階Stokes波作用下管線周?chē)馏w孔壓及其有效應(yīng)力的分布情況.胡翔等[9]采用準(zhǔn)靜態(tài)3維數(shù)值分析方法，研究了波浪載荷作用下海床的孔壓變化規(guī)律及飽和砂質(zhì)海床中的單樁響應(yīng)問(wèn)題.陳寶清等[10]采用計(jì)算流體力學(xué)的OpenFOAM軟件，建立了波浪-海床-結(jié)構(gòu)物數(shù)值模型，模擬了在結(jié)構(gòu)物影響下的海床響應(yīng).Jeng等[11]提出一種新的數(shù)值模型，分別探討了波浪作用下各項(xiàng)同性與各項(xiàng)異性、均質(zhì)與非均質(zhì)海床中海底管道附近的孔壓、有效應(yīng)力以及管道自身的應(yīng)力響應(yīng)，并分析了不同波浪參數(shù)、海床土體參數(shù)、管道尺寸對(duì)計(jì)算結(jié)果的影響.

上述研究主要集中于水平海床(相當(dāng)于遠(yuǎn)海環(huán)境)在周期性波浪作用下的響應(yīng)，對(duì)孤立波的研究主要聚焦在模擬水動(dòng)力學(xué)問(wèn)題的破碎、爬升及回落過(guò)程.最近，有學(xué)者開(kāi)展了對(duì)于近岸斜坡海床及結(jié)構(gòu)物響應(yīng)的研究.Zhao等[12]建立了2維數(shù)值模型，研究在孤立波作用下1∶6和1∶15兩種傾斜度多孔海床的累積孔隙水壓對(duì)液化范圍的影響.Zhao等[13]采用數(shù)值方法研究在波浪作用下防波堤附近傾斜海床的土體響應(yīng)，針對(duì)海床滲透系數(shù)、波浪周期等參數(shù)進(jìn)行分析，并對(duì)比了防波堤附近不同位置的孔壓.Young等[14]在48 m長(zhǎng)的水槽中模擬了海嘯波(破碎孤立波)對(duì)近岸沙灘的侵蝕和沉積作用.Xiao等[15]利用物理實(shí)驗(yàn)數(shù)據(jù)，對(duì)海嘯作用下的沙質(zhì)海床進(jìn)行數(shù)值模擬，分析不同波高和海床滲透系數(shù)對(duì)液化深度的影響.Gao等[16]設(shè)計(jì)管土相互作用設(shè)備進(jìn)行足尺物理試驗(yàn)，研究波浪作用下鋪設(shè)于傾斜海床表面的海底管道上坡及下坡的穩(wěn)定性.

本文針對(duì)近岸斜坡海床，建立孤立波-斜坡海床-海底管道耦合模型.采用N-S方程和k-ε湍流模型對(duì)孤立波進(jìn)行數(shù)值模擬，并將得到的斜坡表面波壓力作用于海底斜坡模型.對(duì)于斜坡土體，基于Biot多孔彈性固結(jié)理論建立海床控制方程；對(duì)于海底管道，根據(jù)線彈性理論建模，模擬分析埋管斜坡海床和管道在孤立波作用下的響應(yīng).通過(guò)對(duì)比數(shù)值模擬結(jié)果與試驗(yàn)數(shù)據(jù)及其解析結(jié)果，驗(yàn)證所建模型的準(zhǔn)確性.根據(jù)計(jì)算結(jié)果，分析在孤立波作用下斜坡上波壓力的水平分布隨時(shí)間的變化特點(diǎn)及不同位置的孔壓隨時(shí)間的變化規(guī)律.此外，對(duì)斜坡模型進(jìn)行參數(shù)分析，研究滲透系數(shù)對(duì)孔壓縱向分布的影響.

1 模型與計(jì)算方法

孤立波在傳播至近岸斜坡時(shí)會(huì)發(fā)生波浪的破碎，并在水平面以上的斜坡繼續(xù)爬升，至最高點(diǎn)后回落.為了準(zhǔn)確地模擬孤立波在斜坡上的破碎、爬升及回落過(guò)程中波壓力的變化，且以此為基礎(chǔ)研究海床的孔壓響應(yīng)，將模型簡(jiǎn)化為波浪、斜坡海床、海底管道3個(gè)部分進(jìn)行耦合計(jì)算.

首先，通過(guò)Flow-3D軟件中的水槽模擬孤立波在斜坡上的破碎、爬升及回落過(guò)程，獲得斜坡表面的水壓力，以水壓力減去初始靜止水壓得到波浪壓力pb；然后，采用COMSOL軟件對(duì)斜坡海床和海底管道進(jìn)行建模，以邊界條件的方式利用插值函數(shù)將pb施加于斜坡模型的上邊界(水土交界處)，計(jì)算與分析斜坡土體在孤立波作用下的孔壓響應(yīng).孤立波-斜坡海床-海底管道耦合模型如圖1所示，模型坐標(biāo)原點(diǎn)O設(shè)置于水平海床底面及傾斜斜坡底面的交界處，孤立波由左側(cè)產(chǎn)生，x軸為波浪行進(jìn)方向，z軸為波浪高度方向.

圖1 孤立波-斜坡-海底管道耦合模型Fig.1 The solitary wave-sloping seabed submarine-pipeline coupling model

1.1 孤立波模型

在孤立波數(shù)值模型中，液體的運(yùn)動(dòng)控制方程采用N-S方程，選取k-ε湍流模型以更真實(shí)地反映波浪與海底斜坡之間的相互作用.Flow-3D軟件中用FAVOR法劃分網(wǎng)格會(huì)產(chǎn)生流體面積和體積分?jǐn)?shù)項(xiàng)，其連續(xù)性方程及動(dòng)量方程表達(dá)式如下：

(1)

(2)

(3)

式中：Ax,Az分別為x和z方向流體面積分?jǐn)?shù)；φF為流體體積分?jǐn)?shù)；vx,vz分別為x和z方向流體速度；ρ為流體密度；p為孔隙水壓力；μ為動(dòng)力黏度；g為重力加速度；ax,az分別為x和z方向黏性加速度.

連續(xù)性方程和動(dòng)量方程的k-ε湍流模型：

PT+fT+DkT-εT

(4)

(5)

式中：kT為湍流動(dòng)能；εT為湍流動(dòng)能耗散率；PT為湍流動(dòng)能產(chǎn)生項(xiàng)；fT為湍流浮力產(chǎn)生項(xiàng)；DkT為湍流動(dòng)能擴(kuò)散項(xiàng)；DεT為湍流動(dòng)能耗散率的擴(kuò)散項(xiàng)；C1,C2,C3為無(wú)量綱計(jì)算系數(shù).

孤立波模型的左側(cè)為波浪入射邊界，基于孤立波理論造波；孤立波模型的右側(cè)為波浪流出邊界(事實(shí)上，由于斜坡長(zhǎng)度足夠長(zhǎng)，保證波浪爬坡后會(huì)完全回落，并不會(huì)從右側(cè)流出)；孤立波模型的上部為水和空氣的交界面，氣壓等于一個(gè)標(biāo)準(zhǔn)大氣壓(101.3 kPa)；孤立波模型的底部為水土交界面，采用墻面邊界，邊界上的流體法向速度為0.

1.2 埋管海床模型

1.2.1海床模型 采用COMSOL軟件對(duì)斜坡海床建模，基于Biot固結(jié)方程，利用偏微分方程模塊對(duì)海底斜坡設(shè)置本構(gòu)關(guān)系方程.以邊界條件的方式用插值函數(shù)將pb施加于斜坡模型的上邊界(水土交界處)，計(jì)算斜坡土體在孤立波作用下的孔壓響應(yīng).基于Biot多孔彈性固結(jié)理論的海床控制方程為

(6)

(7)

假設(shè)海床為均勻、飽和多孔、各向同性的彈性土體，其幾何非線性可以忽略，并且土體的滲透系數(shù)k在各方向上皆相同，則土體孔隙中的液體壓力可以表述為

(8)

式中：ns為土體孔隙率；γw為水的重度；βs為孔隙中液體的壓縮系數(shù)，

(9)

Kw為水的體積模量(此處取2 GPa)；pw0為絕對(duì)水壓力；S為飽和度.

(10)

海床底部可以視為不透水的剛性底面.該底面海床土體的水平位移及縱向位移為0，海水的法向流量也為0，因此邊界條件可以表示為

us=ws=0

(11)

在縱向邊界上，x方向的位移為零.

1.2.2管道模型 采用COMSOL軟件對(duì)海底管道進(jìn)行建模，該模型基于線彈性理論，采用與海床控制方程類似的偏微分方程控制海底管道響應(yīng)，

(12)

(13)

式中：up,wp分別為海底管道在x，z方向的位移；Gp為管道材料的切變模量；υp為管道材料的泊松比.

海底管道外壁視為不透水邊界，故管道外壁表面法向流量為0.假設(shè)管道外壁與海床土體之間無(wú)相對(duì)位移，即管道外壁在x、z方向的位移與海床土體在x、z方向的位移相等，則管道外壁的邊界條件可以表示為

?p/?n=0，us=up，ws=wp

(14)

式中：n為垂直于海底管道壁切線的方向.

2 模型驗(yàn)證

2.1 孤立波模型驗(yàn)證

Synolakis[2]針對(duì)孤立波在斜坡上的破碎、爬升及回落過(guò)程，在水槽和固定斜坡上開(kāi)展試驗(yàn)，并采用該試驗(yàn)數(shù)據(jù)驗(yàn)證孤立波模型的可靠性.設(shè)定模型參數(shù)：最大水深為0.29 m，波高水深比為0.28，斜坡傾斜度為1∶20.

孤立波模型與Synolakis試驗(yàn)數(shù)據(jù)[2](TS1)的對(duì)比結(jié)果如圖2所示，其中t為波浪運(yùn)動(dòng)時(shí)間.圖2(a)～(c)所示為孤立波在斜坡上的破碎及爬升過(guò)程，圖2(d)所示為波浪從最高點(diǎn)回落的階段.

圖2 孤立波模型與TS1數(shù)值結(jié)果的對(duì)比Fig.2 Comparison of numerical results between solitary wave numerical model and TS1

由圖2可見(jiàn)， TS1基本分布在孤立波模型模擬的波浪表面，雖然孤立波模型波峰高度的模擬結(jié)果比TS1值略大，但其波浪的爬升高度及波浪形狀與TS1的結(jié)果吻合得較好.因此，孤立波模型能較好地模擬孤立波在斜坡上的破碎、爬升及回落過(guò)程.

2.2 埋管海床模型驗(yàn)證

為了驗(yàn)證埋管海床模型的準(zhǔn)確性，將埋管海床模型的數(shù)值計(jì)算結(jié)果與文獻(xiàn)[17]的解析結(jié)果(TS2)和文獻(xiàn)[18]的試驗(yàn)數(shù)據(jù)(TS3)進(jìn)行對(duì)比.設(shè)波浪參數(shù)：波高為 5.24 mm；水深d為0.533 m；周期為0.9 s；波長(zhǎng)為1.25 m.海底管道參數(shù)：彈性模量Ep為 68 GPa；υp為0.32；管道外徑R為 0.042 m；埋深e為0.659 m.斜坡海床參數(shù)：海床厚度H為0.826 m；ns為0.42；υs為0.33；側(cè)向土壓力系數(shù)為0.41；k為1.1 mm/s；S為1；土體密度ρs為 2 053 kg/m3；Gp為0.64 MPa.

在波浪作用下，管道周?chē)馏w孔壓沿管道圓周分布，埋管海床模型的計(jì)算結(jié)果與TS2及TS3的數(shù)據(jù)對(duì)比結(jié)果如圖3所示，其中p0為初始孔隙水壓力.由圖3可知：埋管海床模型與TS2在數(shù)值計(jì)算結(jié)果及變化趨勢(shì)上均吻合得較好；與TS3在數(shù)值計(jì)算結(jié)果上略有出入，但變化趨勢(shì)基本一致.因此，可以認(rèn)為埋管海床模型能較好地模擬實(shí)際工程的應(yīng)用情況.

圖3 埋管海床模型與TS2及TS3的數(shù)據(jù)結(jié)果對(duì)比圖Fig.3 Comparison of numerical results among seabed model with a buried pipeline, TS2 and TS3

3 計(jì)算結(jié)果與分析

基于上述已驗(yàn)證的孤立波模型和埋管海床模型，建立孤立波-斜坡海床-海底管道耦合模型；模擬孤立波在斜坡上一個(gè)完整的破碎、爬升及回落過(guò)程；討論管道埋置于水平海床、斜坡坡腳及斜坡海岸線附近3種典型情況下，海床和海底管道響應(yīng)；討論埋深、土體參數(shù)等對(duì)海床和海底管道響應(yīng)的影響.孤立波參數(shù)：波高為3 m；d為10 m；ρ為 1 000 kg/m3；Kw為2 GPa.海底管道參數(shù)：Ep為68 GPa；υp為0.32；管道密度ρp為 2 700 kg/m3；R為1 m；e為2 m；管壁厚度為0.1 m.斜坡海床參數(shù)：海床厚度為50 m；斜坡高度為20 m；水平海床長(zhǎng)度為100 m；斜坡長(zhǎng)度為200 m；斜坡坡度為1∶10；υs為0.35；ns為0.425；ρs為 2 053 kg/m3；側(cè)向土壓力系數(shù)為0.41；k為0.1 mm/s；S為1；Gs為5 GPa.

3.1 埋管斜坡海床的響應(yīng)特征

3.1.1海底管道位置及埋深的影響 海底管道位于3個(gè)典型位置，當(dāng)波峰、波谷經(jīng)過(guò)不同埋深的海底管道圓心正上方時(shí)，p沿海底管道圓周分布的情況如圖4所示.隨著埋深的增加，海底管道上下方孔壓峰值的差距逐漸減小，而峰谷值之間的差距卻在增加，p的谷值均出現(xiàn)在海底管道兩側(cè).

當(dāng)海底管道位于水平海床(見(jiàn)圖4(a))與斜坡坡腳(見(jiàn)圖4(b))時(shí)，p的峰谷值之間的差距以及兩個(gè)峰值之間的差距均較小，不同埋深對(duì)p分布的影響也較?。划?dāng)海底管道位于斜坡海岸線附近時(shí)(見(jiàn)圖4(c))，p的峰谷值間的差距以及兩峰值間的差距均較大，且不同埋深對(duì)p分布的影響非常明顯；在孤立波回落階段(見(jiàn)圖4(d))，海底管道周?chē)讐撼霈F(xiàn)負(fù)值，管道正下方出現(xiàn)明顯的谷值，且隨著埋深的增加，管道上下方孔壓差距逐漸減小.

綜上所述，海底管道的埋深越大，其圓周附近的p越小，管道上下方兩個(gè)p峰值之間的差距越小，管道正上方p峰值與其側(cè)面谷值之間的差距越大.相比于水平海床和斜坡坡腳位置時(shí)的情況，斜坡海岸線附近的海底管道周?chē)鷓分布更不均勻，且受管道埋深的影響更為顯著.

圖4 波峰、波谷經(jīng)過(guò)不同埋深管道時(shí)p沿管道圓周的分布情況Fig.4 Distribution of p around pipe for different pipe buried depth and location when peak or trough passes

圖5 不同k下p及沿縱向及圓周的分布情況Fig.5 The vertical distribution and around pipe distribution of p and for different k

圖6 不同S下p及沿縱向及圓周的分布情況Fig.6 The vertical distribution and around pipe distribution of p and for different S

3.2 管道的縱向受力與位移特點(diǎn)

3.2.1管道的縱向受力 在孤立波作用下，當(dāng)海底管道位于水平海床(x1=-50 m)，斜坡坡腳(x2=0 m)以及斜坡海岸線附近(x3=86 m)3個(gè)位置，e=2 m時(shí)，管道的縱向受力F隨時(shí)間的變化趨勢(shì)如圖7所示.其中，F(xiàn)向上為正，向下為負(fù).在波峰到達(dá)海底管道前數(shù)秒，海底管道受到較小的浮力F1；當(dāng)波峰經(jīng)過(guò)海底管道時(shí)，管道受到向下的壓力F2.|F1| 與 |F2| 的值在x3(斜波海岸線附近)處大于x1(水平海床)和x2(斜坡坡腳)處.與x1和x2位置不同，當(dāng)海底管道埋置于x3處時(shí)，在孤立波回落階段會(huì)受到持續(xù)約20 s的較大上浮力.這可能是由于該時(shí)間段內(nèi)在海底管道上方位置，即斜坡海岸線附近，出現(xiàn)“水洼”，導(dǎo)致海底管道上方p大幅下降，從而使海底管道受到的上浮力增加.

在孤立波作用下，埋置于x3位置，e=2，4，6，8 m時(shí)的F隨時(shí)間的變化趨勢(shì)如圖8所示.由圖8可知，e越小的管道受到的|F|越大.在孤立波回落階段，e=2 m處的管道受到的最大上浮力遠(yuǎn)大于其他位置，約為e=4 m處的2倍，e=6 m處的3倍，e=4 m處的4倍.因此，當(dāng)海底管道埋置在斜坡較淺的位置時(shí)，通過(guò)增加埋深可以大幅降低在孤立波作用下受到的上浮力.

圖7 不同位置海底管道的F隨t的變化趨勢(shì)Fig.7 The variation of F with t for different buried location

圖8 x3處不同埋深海底管道的F隨t的變化趨勢(shì)Fig.8 The variation of F near coastline with t for different pipeline depth on x3

3.2.2海底管道的位移 選取海底管道頂端的位移代表海底管道在海床中的整體位移.當(dāng)海底管道位于3種典型位置，e=2 m時(shí)，wp隨時(shí)間的變化曲線如圖9所示.當(dāng)波峰經(jīng)過(guò)海底管道時(shí)，海底管道出現(xiàn)最大沉降，且x3處的海底管道的沉降大于x2及x1處的海底管道的沉降，其值分別為0.032, 0.026與0.021 m.在x3處的波浪回落階段(t=40～60 s)，受到海底管道上方出現(xiàn)波谷(水洼)的影響，海底管道出現(xiàn)的上浮位移達(dá)到0.013 m，海底管道上下浮動(dòng)的位移差值達(dá)到0.045 m，明顯大于其他兩個(gè)位置的位移.

當(dāng)海底管道處于不同位置且有波峰(谷)經(jīng)過(guò)時(shí)，e與wp的關(guān)系分布如圖10所示.其中：x3(爬升，波峰)對(duì)應(yīng)斜坡海岸線附近位置，孤立波爬升過(guò)程中的波峰出現(xiàn)時(shí)刻；而x3(回落，波谷)對(duì)應(yīng)斜坡海岸線附近位置，孤立波回落過(guò)程中的波谷(水洼)出現(xiàn)時(shí)刻.由圖10可知，e的大小對(duì)wp的極值有一定的影響，并且在3個(gè)位置作用相近.e越大，|wp|越小，且基本呈線性變化.

圖9 不同位置，wp隨t的變化曲線Fig.9 The variation of wp with t for different buried locations

圖10 不同位置有波峰(谷)經(jīng)過(guò)時(shí)wp隨e的變化Fig.10 The variation of wp with e for different buried locations when peak or trough passes

4 結(jié)論

(1)e越大，管道圓周附近|p|越小.p分布越不均勻，管道周?chē)鷓的峰值與谷值之間的差距越大.相比于水平海床和斜坡坡腳位置，斜坡海岸線附近管道周?chē)鷓的分布更不均勻，且受到e的影響更為顯著.

(5) 海岸線附近位置并且e較淺的海底管道在孤立波回落階段會(huì)受到較大的、持續(xù)的上浮力，可以通過(guò)增加e，大幅降低浮力的最大值，從而緩解不利情況.

(6)wp主要受管道埋置位置的影響.其中，海底管道在斜坡海岸線附近出現(xiàn)最大上浮和下沉，最大位移差為0.045 m，該數(shù)值遠(yuǎn)大于水平海床與海床斜坡坡腳位置的最大位移差.此外，e越大，|wp|越小.