孤立波作用下埋管斜坡海床及海底管道的響應分析

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

(上海交通大學 a. 土木工程系； b. 海洋工程國家重點實驗室； c. 高新船舶與深海開發裝備協同創新中心， 上海 200240)

隨著人類對海洋的開發，近岸海域內工程結構物的數量日益增加.以斜坡為主的近岸海床是海底油氣管道“登陸”的重要過渡地帶.在近岸淺水環境下，海床及管道的受力狀態易受波浪作用的影響，淺水區域的波浪呈現出諸如孤立波的特殊形態.孤立波在斜坡海床上爬升和回落時，會引起孔隙水壓力變化，可能導致海床的液化與滑坡，進而危及海底管道等近岸結構物.

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

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

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

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

1 模型與計算方法

孤立波在傳播至近岸斜坡時會發生波浪的破碎，并在水平面以上的斜坡繼續爬升，至最高點后回落.為了準確地模擬孤立波在斜坡上的破碎、爬升及回落過程中波壓力的變化，且以此為基礎研究海床的孔壓響應，將模型簡化為波浪、斜坡海床、海底管道3個部分進行耦合計算.

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

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

1.1 孤立波模型

在孤立波數值模型中，液體的運動控制方程采用N-S方程，選取k-ε湍流模型以更真實地反映波浪與海底斜坡之間的相互作用.Flow-3D軟件中用FAVOR法劃分網格會產生流體面積和體積分數項，其連續性方程及動量方程表達式如下：

(1)

(2)

(3)

式中：Ax,Az分別為x和z方向流體面積分數；φF為流體體積分數；vx,vz分別為x和z方向流體速度；ρ為流體密度；p為孔隙水壓力；μ為動力黏度；g為重力加速度；ax,az分別為x和z方向黏性加速度.

連續性方程和動量方程的k-ε湍流模型：

PT+fT+DkT-εT

(4)

(5)

式中：kT為湍流動能；εT為湍流動能耗散率；PT為湍流動能產生項；fT為湍流浮力產生項；DkT為湍流動能擴散項；DεT為湍流動能耗散率的擴散項；C1,C2,C3為無量綱計算系數.

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

1.2 埋管海床模型

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

(6)

(7)

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

(8)

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

(9)

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

(10)

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

us=ws=0

(11)

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

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

(12)

(13)

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

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

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

(14)

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

2 模型驗證

2.1 孤立波模型驗證

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

孤立波模型與Synolakis試驗數據[2](TS1)的對比結果如圖2所示，其中t為波浪運動時間.圖2(a)～(c)所示為孤立波在斜坡上的破碎及爬升過程，圖2(d)所示為波浪從最高點回落的階段.

圖2 孤立波模型與TS1數值結果的對比Fig.2 Comparison of numerical results between solitary wave numerical model and TS1

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

2.2 埋管海床模型驗證

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

在波浪作用下，管道周圍土體孔壓沿管道圓周分布，埋管海床模型的計算結果與TS2及TS3的數據對比結果如圖3所示，其中p0為初始孔隙水壓力.由圖3可知：埋管海床模型與TS2在數值計算結果及變化趨勢上均吻合得較好；與TS3在數值計算結果上略有出入，但變化趨勢基本一致.因此，可以認為埋管海床模型能較好地模擬實際工程的應用情況.

圖3 埋管海床模型與TS2及TS3的數據結果對比圖Fig.3 Comparison of numerical results among seabed model with a buried pipeline, TS2 and TS3

3 計算結果與分析

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

3.1 埋管斜坡海床的響應特征

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

當海底管道位于水平海床(見圖4(a))與斜坡坡腳(見圖4(b))時，p的峰谷值之間的差距以及兩個峰值之間的差距均較小，不同埋深對p分布的影響也較?。划敽５坠艿牢挥谛逼潞０毒€附近時(見圖4(c))，p的峰谷值間的差距以及兩峰值間的差距均較大，且不同埋深對p分布的影響非常明顯；在孤立波回落階段(見圖4(d))，海底管道周圍孔壓出現負值，管道正下方出現明顯的谷值，且隨著埋深的增加，管道上下方孔壓差距逐漸減小.

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

圖4 波峰、波谷經過不同埋深管道時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 管道的縱向受力與位移特點

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

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

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

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

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

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

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

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

4 結論

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

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

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