基于Pasternak 海床模型的橢圓余弦波浪荷載作用下埋置管線動力響應解析解

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

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

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

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

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

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

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

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

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

本節針對因波浪荷載對埋設管線產生的滲流力(圖2)進行研究，本文計算模型中假定：

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

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

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

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

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

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

式中： ?為速度勢；n為常數，其中n=0, 1, 2, …；如圖2 所示，以管線中間圓心處為坐標原點，埋置方向為x軸；垂直于x軸方向的為z軸；t為時間。

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

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

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

根據上述條件，可以推出自由水面的非線性二階波動方程：

求解式(4)，可得：

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

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

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

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

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

將 ?用復數表示，有：

式中，i 為虛數單位。

取η的一階量，則有：

不同的模數κ決定了波面的曲線形狀，其與波浪的要素有以下的關系：

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

由V=L/T，有：

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

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

當波浪在海床上傳播時，對埋設管線產生滲流壓力，結合上述假設，根據Biot 固結理論，海床土體的控制方程可以表示為：

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

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

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

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

式中，d為管線埋深。

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

式中，C=，m為波數，f為波頻率。

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

式中，D為管線直徑。

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

2 管土相互作用研究

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

圖3 動力海床模型圖Fig.3 Dynamic seabed foundation models

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

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

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

如圖3 所示，由第1 節假設4)，根據豎向力的平衡條件，得到第一個平衡方程：

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

整理得到：

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

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

對式(29)的兩邊求關于x的偏導，代入式(27)可以得到：

根據梁的初等變形理論，管線動力響應控制方程可寫為：

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

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

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

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

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

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

假設管線在施加動態載荷之前處于靜止狀態，因此，管線位移和管線速度的初始條件由下式得：

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

式中，n=0、1、2、3 時依次對應管線的撓度、轉角、彎矩和剪力。

因為微分控制方程是四階線性偏微分方程，很難直接求解，故本文引入兩種Fourier 變換和Laplace變換來進行簡化，將等式變成代數方程。為方便后續推導，定義F[·]和F-1[·]分別表示Fourier 變換及其逆變換，L[·]和L-1[·]分別表示Laplace 變換及其逆變換。

針對動力Pasternak 海床的基本運動方程，如式(33)，將等式兩邊進行關于空間域x的Fourier變換如下：

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

針對時域t，將等式兩端進行Laplace 變換，可以得到：

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

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

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

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

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

令：

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

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

管線撓度響應的關于空間域和時間域的最終解可以表示為(c2≤4kDρA)：

結合式(44)，當實 際工況下存在c2>4kDρA時，撓度公式將變為如下形式：

其中：

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

3 算例驗證

3.1 數值模型驗證

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

本文工程案例選址在渤海南部淺海水域的埕島油田[32]，其位于現代黃河三角洲平原的水下三角洲平原上；水深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；海土體為粉質黏土，天然重度為1.84×104N/m3，彈性模量為1.0×106N/m2，泊松比μ為0.42，海床滲透系數ks為0.5×10-6m/s，孔隙率np為0.66，黏聚力為11.1 kPa，摩擦角為8.8°，不排水抗剪強度為13 kN/m3。

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

圖4 波壓力計算CFD 模型圖Fig.4 Model diagram of CFD

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

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

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

圖6 管線撓度對比圖Fig.6 Comparison of pipeline displacement

圖7 管線彎矩對比圖Fig.7 Comparison of pipeline bending moment

3.2 模型試驗算例驗證1

為進一步驗證本文解析解的可靠性，將本文解析解與既有的模型試驗結果進行對比。LIU 等[33]研究分析了海洋波浪誘導下海床內孔隙壓力變化特征，其在一個垂直圓柱體內填充了1.8 m 厚的砂質沉積物，上方注入0.2 m 的水。通過包含橡膠氣囊的氣囊罐將空氣壓力轉換成水壓，在諧波動波壓力上施加額外的靜態孔隙壓力，以此來模擬真實的深水深度?？紫秹毫τ嬜院４脖砻嫫鹨来蜗蛳戮鶆虿贾迷谏百|沉積物中以監測實時的孔壓變化。實際工況下的計算參數為：土壤的彈性模量為3.3×1011N/m2；海床的滲透系數為1.8×10-4m/s；土壤密度為1960 kg/m3；泊松比為0.3；孔隙率為0.425；波浪高度為3.5 m；周期9 s；水深5.2 m。

圖9 顯示了在2 個不同深度的超孔隙壓力的時間變化。可以看到，解析值和試驗結果之間匹配良好，其中圖9(b)擬合結果略好于圖9(a)，原因可能是，越接近海床表面，邊界效應對于孔壓的影響較為明顯。這也進一步表明，海床深度越深，波浪的影響越小。本文解析解能恰當地捕捉真實的水動力波浪荷載作用下海床土壤的孔隙壓力響應，驗證了本文解析解的適用性。

圖9 試驗算例驗證Fig.9 Model test example

3.3 模型試驗算例驗證2

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

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

圖10 為海床表面以下0.16 m 和0.19 m 處的滲流壓力監測值與本文解析解對比?？梢钥闯霰疚慕馕鼋廨^準確地預測了海床內不同深度處滲流壓力隨時間變化的趨勢，滲流壓力的大小隨著深度的增加而有所減小，圖10(b)較圖10(a)的擬合精度更高，這是由于P10 測點更接近上方的埋置管線，埋置管線自身會對滲流壓力產生散射，從而影響到P10 測點的監測值，而P11 相較于P10測點則更加遠離上方埋置管線，受其影響較小。

圖10 試驗算例驗證2Fig.10 Model test example 2

4 參數分析

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

4.1 管線埋置深度d 的影響

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

圖11 不同埋置深度下動力響應對比Fig.11 Comparison under different embedding depth

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

4.2 波浪高度H 的影響

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

圖12 不同波浪高度下動力響應對比Fig.12 Comparison under different wave height

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

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

4.3 剪切層模量G 的影響

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

圖13 不同剪切層模量下動力響應對比Fig.13 Comparison under different shear modulus

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

5 結論

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

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

(2) 波浪高度H會影響橢圓余弦波的非線性形態，波高增加直接會導致波浪力增強，管線的動態響應特性更加明顯，尤其是管線的轉角、彎矩和剪力的動態響應更為顯著，但波高變化不影響管線動力響應的相位差值。

(3) 剪切層模量G決定了Pasternak 動力海床模型預測精度，Winkler 動力海床模型(G=0)會高估管線的動態響應特性，而過高的剪切層模量則又會低估管線的動態響應特性，故管土作用海床計算模型的選擇和海床物理參數的取值均需謹慎。

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