張治國，葉 銅，張成平，PAN Yu-tao，沈安鑫，吳鐘騰

(1.上海理工大學(xué)環(huán)境與建筑學(xué)院，上海 200093；2.自然資源部丘陵山地地質(zhì)災(zāi)害防治重點(diǎn)實(shí)驗(yàn)室 福建省地質(zhì)災(zāi)害重點(diǎn)實(shí)驗(yàn)室，福州 350002；3.國家海洋局北海預(yù)報(bào)中心，山東省海洋生態(tài)環(huán)境與防災(zāi)減災(zāi)重點(diǎn)實(shí)驗(yàn)室，青島 266061；4.北京交通大學(xué)城市地下工程教育部重點(diǎn)實(shí)驗(yàn)室，北京 100044；5.Department of Civil and Environmental Engineering, National University of Singapore, Singapore 119077)

為解決日益嚴(yán)峻的“能源危機(jī)”，開發(fā)海洋石油天然氣的步伐正逐漸加快，而海底管線作為經(jīng)濟(jì)高效的工程設(shè)施，逐漸成為了能源運(yùn)輸?shù)氖滓x擇?？紤]到施工難度以及后期復(fù)雜海洋環(huán)境的影響，海洋管線大多數(shù)情況下采用淺埋的方式(圖1)。由于海床的可滲透特性，波浪傳播過程中在海床表面產(chǎn)生的隨時(shí)間周期性變化的波壓力，將會(huì)進(jìn)一步向海床中傳遞，改變海床中有效應(yīng)力的分布，對(duì)埋置管線產(chǎn)生滲流壓力。與此同時(shí)，管線使用過程中管土相互作用以及海床土體的液化影響，容易導(dǎo)致管線自身出現(xiàn)穩(wěn)定性問題，尤其在淺水區(qū)，管線受波浪荷載的影響較大，循環(huán)的波浪荷載會(huì)對(duì)管線的結(jié)構(gòu)強(qiáng)度與疲勞壽命[1-3]產(chǎn)生較大的影響。

圖1 波浪荷載下埋置管線示意圖Fig.1 Buried pipe under wave loading

目前，針對(duì)海床、隧道或管線在動(dòng)荷載作用下的響應(yīng)研究方法主要有：理論解析[4-6]、模型試驗(yàn)[7-9]、數(shù)值模擬分析[10-13]等。在理論解析研究中，文獻(xiàn)[4 - 5]針對(duì)陸地地震波為輸入對(duì)象，沒有考慮海洋波浪荷載的動(dòng)力響應(yīng)，且文獻(xiàn)[4 - 6]均以隧道或管線橫截面為分析對(duì)象，而沒有考慮結(jié)構(gòu)物的縱向動(dòng)力響應(yīng)。淺水區(qū)由于極端環(huán)境造成管線破壞而帶來的經(jīng)濟(jì)損失以及原油泄漏導(dǎo)致的環(huán)境污染的現(xiàn)象一直頻繁出現(xiàn)[14-16]。因此，在管線的前期設(shè)計(jì)中，波浪力的分析是不可或缺的一個(gè)步驟。

既有研究中，學(xué)者采用了不同的波浪理論來描述波浪場(chǎng)。XU 等[17]在研究規(guī)則波作用下，海底管線周圍海床演變行為和管線水動(dòng)力的實(shí)驗(yàn)研究時(shí)，采用Stokes 的二階波動(dòng)理論用于評(píng)估近床水平和垂直流速，并與實(shí)驗(yàn)結(jié)果進(jìn)行了對(duì)比分析。王小雯和張建民[18]基于JONSWAP 頻譜模擬隨機(jī)波浪，采用砂土振動(dòng)液化大變形本構(gòu)結(jié)合Biot 動(dòng)力固結(jié)理論，給出了波浪作用下飽和砂質(zhì)海床土體中超靜孔壓瞬態(tài)變化與液化過程彈塑性動(dòng)力變化規(guī)律；MCDOUGAL 等[19]根據(jù)線性波浪方程在淺水區(qū)建立了分析模型，用于估算土壤中孔隙壓力和埋置管線上產(chǎn)生的壓力，這往往會(huì)引起較大的預(yù)測(cè)誤差，故波浪的非線性對(duì)淺水區(qū)滲流壓力的影響是不可忽視的。

現(xiàn)有的有關(guān)地基梁動(dòng)荷載作用下振動(dòng)特性的研究中一般將海床土體看作動(dòng)力Winkler 海床模型。YU 等[20]基于動(dòng)力Winkler 地基模型，推導(dǎo)出無限長梁在任意動(dòng)荷載作用下的動(dòng)力響應(yīng)的解析解，通過考慮脈沖載荷和時(shí)滯載荷這一類的特殊動(dòng)態(tài)載荷，將退化解與已有研究成果進(jìn)行比較，驗(yàn)證了所提出解的有效性。ZHEN 等[21]研究了在簡(jiǎn)諧運(yùn)動(dòng)載荷作用下，位于非線性Winkler 地基上的無限長Euler-Bernoulli 梁的穩(wěn)態(tài)響應(yīng)，發(fā)現(xiàn)基礎(chǔ)剛度的非線性部分不僅影響無限Euler-Bernoulli 梁在簡(jiǎn)諧移動(dòng)荷載作用下的定性分析結(jié)果，而且影響定量分析結(jié)果；YANG 等[22]引入動(dòng)力Winkler地基模型，研究了瑞利地震波作用下長隧道的縱向地震響應(yīng)，發(fā)現(xiàn)土-結(jié)構(gòu)物的相對(duì)剛度比和波頻對(duì)隧道縱向動(dòng)力響應(yīng)的影響；馮浩等[23]得到了動(dòng)力Winkler 地基模型地基上受恒定軸向壓力和橫向行波作用下無限長Euler-Bernoulli 梁的解析解，并與數(shù)值模擬結(jié)果進(jìn)行了對(duì)比驗(yàn)證。上述文獻(xiàn)中將地基看作動(dòng)力Winkler 模型，土體單元各自獨(dú)立，未能體現(xiàn)管土作用導(dǎo)致土體變形連續(xù)性[24]，與實(shí)際土體受力狀態(tài)不符。

本文采用兩階段解析分析法進(jìn)行波浪荷載作用下埋置管線的動(dòng)力響應(yīng)研究。首先，在既有的研究成果的基礎(chǔ)上，考慮土骨架的變形以及孔隙水的壓縮性，基于橢圓余弦波理論和經(jīng)典Biot 固結(jié)理論，推導(dǎo)出了埋置管線所受關(guān)于時(shí)間因素的縱向水平分布滲流壓力；然后，引入三參數(shù)的動(dòng)力Pasternak 海床模型，在充分考慮海床變形連續(xù)性的條件下，考慮海床土體流變阻尼帶來的黏彈性特性，基于Euler–Bernoulli 梁理論，結(jié)合第一階段推導(dǎo)得到的滲流壓力，得到梁的動(dòng)態(tài)微分平衡方程，并依據(jù)Fourier 變換和Lapalce 變換將高階微分控制方程簡(jiǎn)化為代數(shù)方程進(jìn)行求解。根據(jù)卷積定理，得到梁在橢圓余弦波浪荷載作用下?lián)隙?、速度、轉(zhuǎn)角、彎矩和剪力的動(dòng)態(tài)響應(yīng)解。

1 滲流力F(x, t)推導(dǎo)

本節(jié)針對(duì)因波浪荷載對(duì)埋設(shè)管線產(chǎn)生的滲流力(圖2)進(jìn)行研究，本文計(jì)算模型中假定：

圖2 滲流示意圖Fig.2 Seepage schematic diagram

1) 無限長管線埋置于均質(zhì)、剛度和滲透性均各向同性的無限海床中，考慮土骨架的變形和孔隙水的壓縮性，壓縮系數(shù)為常數(shù)；海床內(nèi)的滲流滿足Darcy 定律，海床土壤的滲透系數(shù)為常數(shù)。

2) 采用適用于淺水區(qū)(水深與波長之比小于0.04)的非線性橢圓余弦波[25]，不考慮波浪力傳播過程中的能量損耗。

3) 主要研究埋置管線循環(huán)波浪荷載下土體的豎向相互作用，不考慮管土橫向作用。

4) 假設(shè)埋置管線所受的滲流壓力為自由海床時(shí)由海床表面的波壓力傳遞至管壁時(shí)的周期孔隙水壓，并直接參與管土相互作用。

在淺水區(qū)域，波浪的非線性特性較為明顯，故引入橢圓余弦波理論，取冪級(jí)數(shù)作為勢(shì)函數(shù)[26]：

式中： ?為速度勢(shì)；n為常數(shù)，其中n=0, 1, 2, …；如圖2 所示，以管線中間圓心處為坐標(biāo)原點(diǎn)，埋置方向?yàn)閤軸；垂直于x軸方向的為z軸；t為時(shí)間。

假定水平向的無窮遠(yuǎn)處為靜止水平面，圖2中的波浪的自由水面z=η+h處的邊界條件為：

式中：g為重力加速度； η為靜止水平面以上的波面高度；vx和vz分別為水質(zhì)點(diǎn)的水平流速和豎向流速；h為水深。

淺水區(qū)波浪的豎向流速遠(yuǎn)低于水平流速，忽略vz，將vx用線性化的水平分速取代，式(2)可轉(zhuǎn)換為：

根據(jù)上述條件，可以推出自由水面的非線性二階波動(dòng)方程：

求解式(4)，可得：

式中：(sn)、(cn)、(dn)分別為模數(shù)為κ的橢圓正弦，余弦及模弦函數(shù)：

F(κ)和E(κ)分別為模數(shù)為κ的第一類和第二類完全橢圓積分：

F(κ′)和E(κ)分別為模數(shù)為κ′的第一類和第二類完全橢圓積分，并有κ′2+κ2=1；H為波浪高度；m為波數(shù)[27]；f=，f為波浪頻率，T為波浪周期。

由于忽略掉了豎向流速vz的影響，僅取一階近似橢圓余弦波的近似，只對(duì)式(6)首項(xiàng)進(jìn)行積分，有：

式中，ε=e-πF(κ′)/F(κ)。

將 ?用復(fù)數(shù)表示，有：

式中，i 為虛數(shù)單位。

取η的一階量，則有：

不同的模數(shù)κ決定了波面的曲線形狀，其與波浪的要素有以下的關(guān)系：

橢圓余弦波的波速V可由下式計(jì)算得到[28]：

由V=L/T，有：

模數(shù)κ可由式(13)迭代求得，則此時(shí)橢圓余弦波的波面形狀即可確定。一般將海床表面的孔隙水壓近似等于波浪在海底產(chǎn)生的超靜波壓強(qiáng)Pw，有：

式中： γw為海水重度； ρw為海水密度。

當(dāng)波浪在海床上傳播時(shí)，對(duì)埋設(shè)管線產(chǎn)生滲流壓力，結(jié)合上述假設(shè)，根據(jù)Biot 固結(jié)理論，海床土體的控制方程可以表示為：

可變形介質(zhì)中的可壓縮性流體控制方程可以表述為：

式中：np為孔隙率；K為海水的體積彈性模量。

將式(17)和式(18)中的位移項(xiàng)消除，可得到關(guān)于孔隙水壓p的控制方程：

計(jì)算自由海床時(shí)波浪引起的滲流壓力p1的控制方程以及邊界條件為：

式中，d為管線埋深。

由式(20)可知方程的解p1=p1(x,z,t)，通過分離變量法即可求得自由海床下引起的滲流壓力為：

式中，C=，m為波數(shù)，f為波頻率。

基于已求得的滲流壓力p1，可以得到管線上方自由海床的孔隙水壓F1為：

式中，D為管線直徑。

作用在管線上的滲流壓力F可以表示為：

2 管土相互作用研究

在海洋環(huán)境中，海洋軟土經(jīng)常受到波浪等循環(huán)荷載作用，產(chǎn)生累積變形(沉降)，累積變形主要是累積剪切變形。因此，考慮海床剪切特性的Pasternak 海床模型更適合于海洋環(huán)境。而海洋軟土的流變特性也尤其顯著，本文在Pasternak 模型的彈簧元件k一側(cè)并聯(lián)增加一個(gè)黏壺c，形成一個(gè)Kelvin 模型，并與剪切層G形成了一個(gè)三參數(shù)的黏彈性Pasternak 海床模型，如圖3(a)所示。

圖3 動(dòng)力海床模型圖Fig.3 Dynamic seabed foundation models

黏彈性Pasternak 海床模型，是在黏彈性Winkler 海床模型的彈簧單元上增加一個(gè)只產(chǎn)生剪切變形而不可壓縮的剪切層，將原本獨(dú)立的土彈簧單元連接起來，以此體現(xiàn)相鄰彈簧間的變形連續(xù)特性，還原土體的剪切作用；與此同時(shí)，這兩者也均是梁的瞬態(tài)響應(yīng)、譜關(guān)系以及振動(dòng)特性問題中常見的黏彈性海床模型。當(dāng)因波浪荷載產(chǎn)生的循環(huán)滲流力荷載在黏彈性的海床梁中傳播時(shí)，不可避免地會(huì)涉及到海床的剪切特性帶來的變形連續(xù)響應(yīng)，考慮了海床土連續(xù)性的Pasternak 海床模型相較于Winkler 海床模型顯得更為貼近實(shí)際工況。管線假定為Euler–Bernoulli 梁，其擱置于動(dòng)力Pasternak 海床模型上。

取圖3(b)中Euler-Bernoulli 梁上任意截面x處的微元段dx作為隔離體進(jìn)行受力分析，其截面上存在彎矩Mz、剪力Qz、外部滲流力動(dòng)荷載F(x,t)、地基反力q(x,t)、軸向拉力N以及動(dòng)荷載帶來的附加效應(yīng)慣性力m(x,t)，在這里慣性力可以表示為[20]：

式中：m為管線單位長度質(zhì)量； ρ為管線的材料密度；A為管線的截面面積；w(x,t)定義為管線的垂直撓度，其中荷載沿x軸分布，t為時(shí)間。

如圖3 所示，由第1 節(jié)假設(shè)4)，根據(jù)豎向力的平衡條件，得到第一個(gè)平衡方程：

式中，F(xiàn)(x,t)為式(24)求得的海浪引起的滲流力荷載。

整理得到：

由力矩平衡條件，對(duì)圖3(b)中右下角點(diǎn)O取彎矩，得到第二個(gè)平衡方程：

忽略掉式(28)的微分平方項(xiàng)，整理可得：

對(duì)式(29)的兩邊求關(guān)于x的偏導(dǎo)，代入式(27)可以得到：

根據(jù)梁的初等變形理論，管線動(dòng)力響應(yīng)控制方程可寫為：

式中：EI為管線等效抗彎剛度；E為管線彈性模量；I為截面慣性矩。

如圖3(a)所示，本文動(dòng)力Pasternak 海床模型海床反力q(x,t)為：

式中：k為海床彈簧剛度；c為海床黏滯系數(shù)；G為Pasternak 海床模型剪切層的模量。

將式(32)代入式(31)，可得波浪荷載作用下埋置管線基于動(dòng)力Pasternak 海床模型的控制方程：

當(dāng)剪切層參數(shù)G=0 時(shí)，動(dòng)力Pasternak 海床模型則退化為動(dòng)力Winkler 海床模型(廣義Kelvin 模型)，海床彈簧系數(shù)k、剪切模量G分別按ATTEWELL等[29]、TANAHASHI[30]的代表性公式進(jìn)行取值：

式中：Es為土體彈性模量；μ為土體泊松比；t為剪切層厚度，根據(jù)文獻(xiàn)[31]，t取2.5 倍的管徑，即t=2.5D。

假設(shè)管線在施加動(dòng)態(tài)載荷之前處于靜止?fàn)顟B(tài)，因此，管線位移和管線速度的初始條件由下式得：

由于管線是無限的，兩端的邊界條件為：

式中，n=0、1、2、3 時(shí)依次對(duì)應(yīng)管線的撓度、轉(zhuǎn)角、彎矩和剪力。

因?yàn)槲⒎挚刂品匠淌撬碾A線性偏微分方程，很難直接求解，故本文引入兩種Fourier 變換和Laplace變換來進(jìn)行簡(jiǎn)化，將等式變成代數(shù)方程。為方便后續(xù)推導(dǎo)，定義F[·]和F-1[·]分別表示Fourier 變換及其逆變換，L[·]和L-1[·]分別表示Laplace 變換及其逆變換。

針對(duì)動(dòng)力Pasternak 海床的基本運(yùn)動(dòng)方程，如式(33)，將等式兩邊進(jìn)行關(guān)于空間域x的Fourier變換如下：

考慮到Fourier 變換的微分性質(zhì)以及式(36)的邊界條件，可以得到：

針對(duì)時(shí)域t，將等式兩端進(jìn)行Laplace 變換，可以得到：

考慮Laplace 變換的微分性質(zhì)和式(35)的初始條件，得到：

因此，該問題在頻域中的解可以表示為：

為了獲得時(shí)域中的解，需要將等式(42)進(jìn)行Laplace 逆變換：

式中，δ=c/(2ρA)。

如果ψ(ξ)≥0(c2≤4kDρA)，則有：

令：

并且設(shè)f(t)、h(t)分別表示F(ω)和H(ω)的Laplace逆變換，有：

利用卷積定理，式(43)求解得：

管線撓度響應(yīng)的關(guān)于空間域和時(shí)間域的最終解可以表示為(c2≤4kDρA)：

結(jié)合式(44)，當(dāng)實(shí) 際工況下存在c2>4kDρA時(shí)，撓度公式將變?yōu)槿缦滦问剑?/p>

其中：

需要說明的是，管線的速度響應(yīng)v(x,t)和加速度響應(yīng)a(x,t)可以通過獲取撓度w(x,t)響應(yīng)相對(duì)于時(shí)域t的一階和二階導(dǎo)數(shù)獲得；管線的轉(zhuǎn)角響應(yīng)θ(x,t)、彎矩響應(yīng)M(x,t)和剪力響應(yīng)Q(x,t)也可通過取撓度w(x,t)響應(yīng)相對(duì)于空間域x的一階、二階和三階偏導(dǎo)數(shù)，各乘以抗彎剛度系數(shù)EI得到。

3 算例驗(yàn)證

3.1 數(shù)值模型驗(yàn)證

為驗(yàn)證橢圓余弦波作用下埋置管線動(dòng)力響應(yīng)解析解的正確性，將其與三維數(shù)值模型的運(yùn)行結(jié)果進(jìn)行對(duì)比。采用計(jì)算流體力學(xué)CFD 軟件(Flow3D)和有限元分析FEM 軟件(ABAQUS)相結(jié)合的方式進(jìn)行數(shù)值分析。首先，通過CFD 軟件建立三維數(shù)值波浪水池來模擬淺水區(qū)的波浪場(chǎng)及作用在海床上的波壓力；其次，將得到的海床表面的周期波壓函數(shù)輸入到FEM 軟件中，作為波浪作用下埋置管線動(dòng)力響應(yīng)研究的邊界條件，基于管線與海床的三維有限元模型，對(duì)管土相互作用的動(dòng)力響應(yīng)問題進(jìn)行數(shù)值驗(yàn)證分析。

本文工程案例選址在渤海南部淺海水域的埕島油田[32]，其位于現(xiàn)代黃河三角洲平原的水下三角洲平原上；水深2 m～15 m，面積約1710 km2。管線埋置深度d為0.5 m，管線截面直徑D為0.324 m，壁厚tw為0.0205 m，管線的材料密度ρ 為7850 kg/m3，彈性模量為2.07×1011N/m2，泊松比為0.28；波浪的有效波高2 m，波浪周期為8.6 s；海水密度為1025 kg/m3，體積彈性模量K為2.18×109N/m2；海土體為粉質(zhì)黏土，天然重度為1.84×104N/m3，彈性模量為1.0×106N/m2，泊松比μ為0.42，海床滲透系數(shù)ks為0.5×10-6m/s，孔隙率np為0.66，黏聚力為11.1 kPa，摩擦角為8.8°，不排水抗剪強(qiáng)度為13 kN/m3。

圖4 為CFD 軟件建立的三維數(shù)值波浪水池模型，下部三維海床實(shí)體尺寸為120 m×40 m×10 m，網(wǎng)格單元數(shù)為720 292；上部為流體子域，其網(wǎng)格單元數(shù)為718 676。自由表面追蹤類型為尖銳接觸面(sharp interface)，流體的流動(dòng)類型為不可壓縮；由于本工程背景中流體為海水，接觸的是空氣，故選用單流體法，忽略氣相的作用；物理模型選用泥沙沖刷模型，波浪類型為斯托克斯波和橢圓余弦波(Stokes wave and Cnoidal wave)，這與本文解析解中橢圓余弦波理論相契合；波面的豎向邊界條件為壁面，造波口(wave)及出流口(outflow)為模型的x軸方向(見圖4)，正交y軸方向不設(shè)置波浪。

圖4 波壓力計(jì)算CFD 模型圖Fig.4 Model diagram of CFD

圖5 為FEM 軟件建立的三維有限元模型。三維海床實(shí)體尺寸為60 m×40 m×10 m，實(shí)體海床單元數(shù)為68 280，節(jié)點(diǎn)數(shù)為80 412；管線長度為60 m，實(shí)體單元數(shù)為1320，節(jié)點(diǎn)數(shù)為1544；共計(jì)69 600 個(gè)實(shí)體單元。管線采用線彈性本構(gòu)模型，海床土體采用Prony 級(jí)數(shù)定義的黏彈性本構(gòu)模型進(jìn)行模擬；模型上邊界條件為無約束透水邊界，四周邊界條件為上下自由滑動(dòng)的不透水邊界，底面為固定的不透水邊界；將CFD 模型得到的海床表面波壓力導(dǎo)入FEM 軟件的Amplitude 模塊并施加于海床實(shí)體模型的上表面進(jìn)行數(shù)值計(jì)算。

圖5 管土作用FEM 模型圖Fig.5 Model diagram of FEM

從圖6、圖7、圖8 可以看出，隨著水深的增大，管線的撓度、彎矩、加速度響應(yīng)特性都有所衰減。需要說明的是，根據(jù)式(12)～式(15)及案例實(shí)際波浪參數(shù)，獲得的模數(shù)κ較小，因此，波浪力荷載使管線在撓度、速度、加速度響應(yīng)趨勢(shì)呈現(xiàn)出較強(qiáng)的周期線性分布。水深為3 m、6 m 和9 m對(duì)應(yīng)的管線峰值撓度分別為3.0×10-3m、2.1×10-3m和1.2×10-3m，較水深3 m 的工況，后兩者衰減幅度分別為30%、60%；水深為3 m、6 m 和9 m 對(duì)應(yīng)的管線峰值彎矩分別為1.55 kN·m、1.01 kN·m和0.54 kN·m，較水深3 m 的工況，后兩者衰減幅度分別為35%和65%；水深3 m、6 m、9 m 對(duì)應(yīng)管線峰值加速度分別為6.2×10-4m/s2、3.9×10-4m/s2和2.1×10-4m/s2，較水深3 m 工況，后兩者的衰減幅度分別為37%和66%。總體上，本文解析解與數(shù)值模擬結(jié)果吻合度較高，進(jìn)一步驗(yàn)證了本文解析解的正確性與適用性。

圖6 管線撓度對(duì)比圖Fig.6 Comparison of pipeline displacement

圖7 管線彎矩對(duì)比圖Fig.7 Comparison of pipeline bending moment

3.2 模型試驗(yàn)算例驗(yàn)證1

為進(jìn)一步驗(yàn)證本文解析解的可靠性，將本文解析解與既有的模型試驗(yàn)結(jié)果進(jìn)行對(duì)比。LIU 等[33]研究分析了海洋波浪誘導(dǎo)下海床內(nèi)孔隙壓力變化特征，其在一個(gè)垂直圓柱體內(nèi)填充了1.8 m 厚的砂質(zhì)沉積物，上方注入0.2 m 的水。通過包含橡膠氣囊的氣囊罐將空氣壓力轉(zhuǎn)換成水壓，在諧波動(dòng)波壓力上施加額外的靜態(tài)孔隙壓力，以此來模擬真實(shí)的深水深度?？紫秹毫τ?jì)自海床表面起依次向下均勻布置在砂質(zhì)沉積物中以監(jiān)測(cè)實(shí)時(shí)的孔壓變化。實(shí)際工況下的計(jì)算參數(shù)為：土壤的彈性模量為3.3×1011N/m2；海床的滲透系數(shù)為1.8×10-4m/s；土壤密度為1960 kg/m3；泊松比為0.3；孔隙率為0.425；波浪高度為3.5 m；周期9 s；水深5.2 m。

圖9 顯示了在2 個(gè)不同深度的超孔隙壓力的時(shí)間變化?？梢钥吹?，解析值和試驗(yàn)結(jié)果之間匹配良好，其中圖9(b)擬合結(jié)果略好于圖9(a)，原因可能是，越接近海床表面，邊界效應(yīng)對(duì)于孔壓的影響較為明顯。這也進(jìn)一步表明，海床深度越深，波浪的影響越小。本文解析解能恰當(dāng)?shù)夭蹲秸鎸?shí)的水動(dòng)力波浪荷載作用下海床土壤的孔隙壓力響應(yīng)，驗(yàn)證了本文解析解的適用性。

圖9 試驗(yàn)算例驗(yàn)證Fig.9 Model test example

3.3 模型試驗(yàn)算例驗(yàn)證2

SUN 等[34]利用波浪水槽進(jìn)行了一系列綜合的實(shí)驗(yàn)室試驗(yàn)，以研究海溝層中部分嵌入管道周圍的波浪引起的孔隙壓力。實(shí)驗(yàn)在河海大學(xué)一個(gè)長55 m、高1.3 m、寬1.0 m 的波浪水槽中進(jìn)行。波浪水槽在上游裝有液壓活塞式造波機(jī)，在下游端裝有一個(gè)海綿式消波器，以耗散來波能量，減小波浪反射效應(yīng)。造波機(jī)能夠產(chǎn)生周期為0.6 s～2.5 s規(guī)則波，最大波高為0.2 m。在實(shí)驗(yàn)中，通過使用孔隙壓力傳感器和波高計(jì)同時(shí)測(cè)量回填溝槽中管道周圍的波浪引起的孔隙水壓力變化和水面高度。海床土的粒徑為0.173 mm，單位重度為26.5 kN/m3，海床滲透系數(shù)為3.56×10-5m/s，泊松比為0.32，孔隙率為0.564。

選擇SUN 等[34]全埋管道的測(cè)試工況10 的P10 和P11 測(cè)點(diǎn)的實(shí)測(cè)滲流壓力值與本文解析解進(jìn)行對(duì)比。其溝槽深度0.15 m，回填高度0.15 m，波浪高度為0.14 m，周期為1.6 s。

圖10 為海床表面以下0.16 m 和0.19 m 處的滲流壓力監(jiān)測(cè)值與本文解析解對(duì)比。可以看出本文解析解較準(zhǔn)確地預(yù)測(cè)了海床內(nèi)不同深度處滲流壓力隨時(shí)間變化的趨勢(shì)，滲流壓力的大小隨著深度的增加而有所減小，圖10(b)較圖10(a)的擬合精度更高，這是由于P10 測(cè)點(diǎn)更接近上方的埋置管線，埋置管線自身會(huì)對(duì)滲流壓力產(chǎn)生散射，從而影響到P10 測(cè)點(diǎn)的監(jiān)測(cè)值，而P11 相較于P10測(cè)點(diǎn)則更加遠(yuǎn)離上方埋置管線，受其影響較小。

圖10 試驗(yàn)算例驗(yàn)證2Fig.10 Model test example 2

4 參數(shù)分析

針對(duì)管線埋置深度d、波浪有效波高H、剪切層模量G進(jìn)行敏感參數(shù)分析。相關(guān)計(jì)算參數(shù)設(shè)置參考埕島油田[32]，基本參數(shù)為管線埋置深度d為0.5 m，波高H為2m，剪切層模量G為100 kPa，水深h為3 m，管線外徑D為0.3 m，管線壁厚tw為0.02 m ，管線材料密度ρ 為7850 kg/m3，黏滯系數(shù)c為2000 Pa?s。基于基本工況中求得的峰值浪高z0、海床表面孔壓p0、峰值撓度w0、峰值轉(zhuǎn)角θ0、峰值彎矩M0、峰值剪力Q0、峰值速度v0、峰值加速度a0、周期T=8.5 s 和水深h=3 m 進(jìn)行參數(shù)歸一化處理。

4.1 管線埋置深度d 的影響

分別選取管線埋置深度d為0.5 m、1 m、2 m、4 m 和8 m，在5 種不同工況下的波面方程以及的梁的撓度、轉(zhuǎn)角、彎矩、剪力、速度和加速度的動(dòng)力響應(yīng)如圖11 所示，可以發(fā)現(xiàn)，波面方程與自由海床下滲流壓力的衰減規(guī)律并不受管線埋置深度的影響(見圖11(a)和圖11(b))，但管線的動(dòng)態(tài)響應(yīng)特征則受管線埋置深度的影響尤為顯著，響應(yīng)幅值隨著埋置深度的增加有所減弱，且不同埋置深度會(huì)導(dǎo)致滲流壓力到達(dá)管線的時(shí)間不一致，導(dǎo)致管線動(dòng)力響應(yīng)的周期函數(shù)存在較為明顯的相位差。

圖11 不同埋置深度下動(dòng)力響應(yīng)對(duì)比Fig.11 Comparison under different embedding depth

圖11(c)～圖11(h)中，管線的撓度、轉(zhuǎn)角、彎矩、剪力、速度和加速度響應(yīng)特性均因埋置深度的增加而減小，以轉(zhuǎn)角響應(yīng)為例，埋置深度為0.5 m、1 m、2 m、4 m 和8 m 時(shí)對(duì)應(yīng)的歸一化峰值轉(zhuǎn)角分別為1.0、0.64、0.21、0.072 和0.015，其余四種工況較基本工況減幅分別為36%、79%、92.8%和98.5%，管線埋置深度越深，受波浪滲流壓力影響越小，故在淺水區(qū)域，將管線進(jìn)行合理深度的埋設(shè)有利于管線安全。圖11(c)中可以明顯地看出，這五種工況之間存在大小不等的相位差，隨著埋置深度的增大，管線動(dòng)力響應(yīng)函數(shù)會(huì)在不同工況之間產(chǎn)生相位差。

4.2 波浪高度H 的影響

波浪高度H取為1 m、2 m、4 m、8 m 和16 m，并以波高2 m 為基本工況，五種工況下的波面方程以及梁的撓度、轉(zhuǎn)角、彎矩、剪力、速度和加速度的動(dòng)力響應(yīng)如圖12 所示，可以發(fā)現(xiàn)，波浪高度的變化對(duì)于波面方程以及管線的動(dòng)態(tài)響應(yīng)的影響都顯得較為敏感，尤其是管線的轉(zhuǎn)角、彎矩和剪力的動(dòng)態(tài)響應(yīng)特性。此時(shí)管線的埋置深度以及外徑尺寸保持不變，五種工況下滲流壓力到達(dá)管線的時(shí)間不變，不存在相位差值。

圖12 不同波浪高度下動(dòng)力響應(yīng)對(duì)比Fig.12 Comparison under different wave height

圖12(a)中，波面高度隨著波浪高度的增大而增大，波浪高度為1 m、2 m、4 m、8 m 和16 m時(shí)對(duì)應(yīng)的歸一化峰值高度分別為0.7、1、1.4、2 和2.8，其余四種工況較基本工況的增幅分別為43%、100%、186%和300%。圖12(b)中，自海床表面向下滲流壓力隨著深度的加深有所衰減，波浪高度顯著地影響著海床表面的波壓力，波浪高度為1 m、2 m、4 m、8 m 和16 m 時(shí)對(duì)應(yīng)的歸一化峰值高度分別為0.8、1、1.15、1.4 和1.8，其余四種工況較基本工況的增幅分別為25%、44%、75%和125%。

圖12(c)～圖12(h)中，管線的撓度、轉(zhuǎn)角、彎矩、剪力、速度和加速度響應(yīng)特性均因波高H的增加而增加。以剪力響應(yīng)為例，波浪高度為1 m、2 m、4 m、8 m 和16 m 時(shí)對(duì)應(yīng)的歸一化峰值彎矩分別為0.25、1、3.12、12.56 和48.9，其余四種工況較基本工況依次增加了0.75%、2.87%、12.31%、48.65%。

4.3 剪切層模量G 的影響

在保留其他參數(shù)不改變的情況下，剪切層模量G取值為0 Pa、1.0×104Pa、1.0×105Pa、1.0×106Pa 和1.0×107Pa。五種工況下的波面方程以及梁的撓度、轉(zhuǎn)角、彎矩、剪力、速度和加速度的動(dòng)力響應(yīng)如圖13 所示。剪切層參數(shù)不影響波面方程，管線的動(dòng)態(tài)響應(yīng)特性隨著剪切層參數(shù)的增大而有所衰減。

圖13 不同剪切層模量下動(dòng)力響應(yīng)對(duì)比Fig.13 Comparison under different shear modulus

圖13(a)～圖13(f)中，Winkler 動(dòng)力海床模型(G=0)會(huì)高估管線的動(dòng)態(tài)響應(yīng)特性，而過高的剪切層模量則又會(huì)低估管線的動(dòng)態(tài)響應(yīng)特性，故剪切層模量的取值大小是影響Pasternak 動(dòng)力海床模型預(yù)測(cè)精度的決定性因素。以管線的速度響應(yīng)為例，剪切層模量G取為0 Pa、1.0×104Pa、1.0×105Pa、1.0×106Pa 和1.0×107Pa 對(duì)應(yīng)的歸一化峰值速度分別為為1.3、1.25、1、0.55 和0.15，其余四種工況較基本工況的減小幅度分別為3.8%、23%、58%和88%，故合理的剪切層厚度顯得尤其重要。

5 結(jié)論

考慮到既有波浪-海床-管線相互作用機(jī)制研究中一般將海床土體看作動(dòng)力Winkler 海床模型，不能考慮管土界面接觸的變形連續(xù)性，本文首先基于橢圓余弦波理論和Biot 固結(jié)理論推導(dǎo)得到了埋置管線的滲流力荷載，接著將管線看作Euler-Bernoulli 梁，引入3 參數(shù)的Pasternak 海床模型，結(jié)合已求得的滲流力動(dòng)荷載推導(dǎo)得到管線的微分動(dòng)平衡方程，并最終由Fourier 和Laplace 積分變換求得管線的動(dòng)力響應(yīng)解。通過多個(gè)數(shù)值工況和2 個(gè)模型試驗(yàn)驗(yàn)證可以看出，本文解析方法具備較好的可靠性和合理性，可用于淺海區(qū)海洋管線前期設(shè)計(jì)參考，具有一定的工程實(shí)用意義。通過參數(shù)分析，得出主要結(jié)論如下：

(1) 管線埋置深度d對(duì)管線的動(dòng)態(tài)響應(yīng)特征的影響較為顯著，管線的撓度、轉(zhuǎn)角、彎矩、剪力、速度和加速度響應(yīng)幅值隨著埋置深度的增加有所減弱；隨著埋置深度的增大，海床表面的波壓力向下傳播至管線埋深處的時(shí)間隨之增加，管線動(dòng)力響應(yīng)會(huì)產(chǎn)生相位差。

(2) 波浪高度H會(huì)影響橢圓余弦波的非線性形態(tài)，波高增加直接會(huì)導(dǎo)致波浪力增強(qiáng)，管線的動(dòng)態(tài)響應(yīng)特性更加明顯，尤其是管線的轉(zhuǎn)角、彎矩和剪力的動(dòng)態(tài)響應(yīng)更為顯著，但波高變化不影響管線動(dòng)力響應(yīng)的相位差值。

(3) 剪切層模量G決定了Pasternak 動(dòng)力海床模型預(yù)測(cè)精度，Winkler 動(dòng)力海床模型(G=0)會(huì)高估管線的動(dòng)態(tài)響應(yīng)特性，而過高的剪切層模量則又會(huì)低估管線的動(dòng)態(tài)響應(yīng)特性，故管土作用海床計(jì)算模型的選擇和海床物理參數(shù)的取值均需謹(jǐn)慎。

(4) 淺水區(qū)海底埋置管線受波浪力影響較為顯著，基于Pasternak 地基模型進(jìn)行海床-管線相互作用分析時(shí)，對(duì)部分土體參數(shù)進(jìn)行了簡(jiǎn)化；而深水區(qū)靜水壓力較大造成管線屈服及管線埋置段和非埋段交替出現(xiàn)等復(fù)雜工況。