基于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。基于基本工况中求得的峰值浪高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 地基模型进行海床-管线相互作用分析时，对部分土体参数进行了简化；而深水区静水压力较大造成管线屈服及管线埋置段和非埋段交替出现等复杂工况。