大尺度群桩随机波浪力场的降维模拟

刘章军，董心宇，张文远

（武汉工程大学土木工程与建筑学院，湖北武汉 430074）

引 言

近年来，随着中国海洋战略的提出，大型离岸海工建筑物的荷载设计逐渐成为工程界的研究热点。由于海工结构长期处于复杂海况中，因此，合理地模拟海工基础性结构在服役期间所承受的波浪力是进行结构荷载设计的基础和前提。其中，探讨海工基础性结构大小尺度桩柱波浪力模拟的异同对工程结构设计有着重要的参考意义。

在实际工程中，桩柱的大小尺度由直径与波长之比决定，比值大于0.2 的桩柱被称为大尺度桩柱，反之，则被称为小尺度桩柱［1］。计算不同尺度桩柱的波浪力时所采用的原理和公式不尽相同。对于波浪力作用下的大尺度桩柱，国内外学者多采用绕射理论进行研究。在大尺度单桩方面，MacCamy 等［2］发展了以绕射理论为基础的波浪力计算方法。邱大洪［3］基于一阶近似椭圆余弦波理论研究了浅水区非线性波浪对圆柱桩的绕射作用，改进了作用于大尺度桩柱的波浪力计算公式。王俊杰等［4］应用波浪力计算公式，建立了大尺度桩柱波浪力的Nyström 数值模型。邓莎莎等［5］以大尺度桥墩为研究对象，对该桥墩基础所受的波浪载荷进行了数值模拟。在大尺度群桩方面，Spring 等［6］首先提出了针对群桩结构的矩阵算法。Linton 等［7］对该方法进行了改进，简化了群桩波浪力的计算。陶建华等［8］基于矩阵算法，由谱分析法推导得到群桩系数。季新然等［9］通过物理模型实验模拟了多向不规则波浪与群桩结构的相互作用。小尺度桩柱的波浪力计算则多以Morison 方程作为理论基础［10］，该方程由Morison 于1950年提出。毛鸿飞等［11］基于Morison 方程，对大波幅波浪作用下，非完全淹没水平圆柱的波浪力进行了数值模拟。程永舟等［12］基于波浪力的数值模拟，通过波浪水槽试验，对孤立波作用下垂直桩柱周围的局部冲刷特性进行了分析研究。

然而，上述研究仅对给定尺度的桩柱进行了受力分析，并未对桩柱尺度进行判别。实际上，在工程设计中通常依据设计波长进行大小尺度桩柱的判定，而设计波长则通常根据工程经验来确定，这无疑会对工程应用造成困惑。此外，大小尺度桩柱的判别并非一成不变，其将随着海况变化而发生改变。因此，如何选择波浪力计算公式将非常困难。通常，随机波浪力的模拟多采用传统Monte Carlo 方法。然而，传统Monte Carlo 方法所需的随机变量与样本数量巨大，导致结构动力计算量大，且生成的样本无法在概率密度层面上精细地描述随机过程的概率特性，这为结构精细化的动力反应分析与动力可靠度评价带来挑战。

针对上述研究现状，本文首先提出了随机波浪力作用下大小尺度桩柱的判别方法；其次，引入随机函数降维思想［13⁃14］，实现了仅用一个基本随机变量即可精确地模拟大尺度桩柱波浪力随机场的目的；最后，考察了基于Morison 方程计算大尺度桩柱时所带来的误差，以及群桩效应对桩柱波浪力的影响。值得指出的是，由于本文方法生成的代表性样本集合的概率信息完备，因此可与概率密度演化方法［15⁃16］相结合，进而实现海工结构在波浪力作用下的精细化动力反应分析与动力可靠性评价。

1 基于绕射理论的大尺度桩柱波浪力谱模型

在海工结构中，当桩柱之间的间距较小时，在波浪荷载的作用下，其受力与单桩有很大差别。桩中心距与桩直径的比值不大于4 时，需要考虑周边桩柱对波浪力的影响［17］。同时，大尺度桩柱对波动场的影响不可忽略，即，当波浪接触到桩柱后，将在柱面产生一个向外扩散的波，在入射波与散射波的叠加达到稳态时，会形成一个新的波动场。

现假设水体无黏性且不可压缩，波浪做有势运动，群桩结构由N=p×q根桩柱构成，如图1（a）所示，其中p代表行数，q代表列数。为简便计，本文仅考虑直立圆柱桩的情况。对于任意的两个桩柱i和j，以各自桩柱中心为原点建立局部的极坐标系，桩柱的相对位置及相关参数如图1（b）所示。

图1 群桩波浪力的计算简图Fig.1 Diagram for waves interact with multiple circular cyl⁃inders

假设桩柱i的中心坐标为(xi，yi)，桩柱j的中心坐标为(xj，yj)，根据绕射理论，波动场中任意一点的总速度势空间分量为［7］：

其中：

式中 i 为虚数单位，g为重力加速度，H为波高，z是以海底为坐标原点的竖向坐标，ω为圆频率，d为静水面处水深，β为入射波方向；κ表示波数，由色散关系ω2=gκtanh(κd)求得；Ai，n为待定系数；Hn(κri)表示n阶第一类汉克尔函数；Jn(κri)表示n阶第一类贝塞尔函数；Yn(κri)表示n阶第二类贝塞尔函数；J′n(κai)表示函数Jn(κai)对变量κai的一阶导数；H′n(κai)表示函数Hn(κai)对变量κai的一阶导数，ai为桩柱i的半径。

当该点位于圆柱i的表面时，利用欧拉公式和Graf 加法定理［18］，式（1）可以改写为：

式中m与n均为整数。

在式（5）中，通过对等号左边的无穷项进行截断，如取2M+1 项，则可得到N×(2M+1)个未知数Ai，n的联立方程组，从而求解待定系数Ai，n。根据文献［7］，本文中直接取M=6。

于是，圆柱j周围速度势的空间分量为：

在式（6）中，对贝塞尔函数应用Wronskian 关系式，并对通过速度势求得的桩柱表面压强进行积分和无因次化，得到波浪作用于圆柱j时波浪荷载的传递函数［7］：

式中ρ0为海水密度。

这样，群桩结构中桩柱j任意坐标z处的波浪力功率谱密度函数为：

式中Sη(ω)为波高过程的单边功率谱密度函数。

特别地，当结构为大尺度单桩（即N=1）时［2］，待定系数Aj，n=-in，此时传递函数为：

2 大尺度桩柱波浪力场的降维模拟

根据随机过程模拟的降维方法，作用于大尺度群桩中第j个桩柱上的水平波浪力可表示为［19⁃21］：

式中 Δω为频率离散步长；Nω为频域离散点数；ωu为上限截止频率，且有为一组标准的正交随机变量集，满足以下基本条件（必要条件）：

式中k，l=1，2，…，Nω；δkl为Kronecker 符号。

式中Θ为基本随机变量，服从(0，2π]上的均匀分布；α为任意的确定性参数，在本文中取α=π/4。可以验证，式（12）是满足式（11）的基本条件（必要条件）。

在式（12）中，(=1，2，…，Nω) 与k(k=1，2，…，Nω)是一种确定性的一一映射关系，这种一一映射关系可由MATLAB 工具箱的内置函数，即rand（‘state’，0）和temp = randperm（Nω）来实现。(=1，2，…，Nω)与k(k=1，2，…，Nω)的一一映射关系即为kˉ=temp（k），这恰是降维模拟方法的一个充分条件。

于是，将式（12）代入式（10），即可得到随机波浪力场的降维模拟公式：

式中ϕkˉ(Θ)为基本随机变量Θ的函数，其表达式为：

通过选取基本随机变量Θ的代表性点集，即可得到大尺度桩柱随机波浪力的代表性时程集合，进而结合概率密度演化理论，进行桩柱结构随机波浪力作用下的动力反应分析与可靠度精细化计算。

3 桩柱尺度的判别方法

绕射理论仅适用于直径与波长之比大于0.2 的大尺度桩柱波浪力计算，但对于随机波浪而言，其波长是一个变量，这将导致桩柱尺度的判别较为困难。因此，可以通过确定其特征频率的方式来确定其特征波长。这样，便可通过特征波长来判别其是否属于大尺度桩柱。

对于窄谱波浪，其能量主要集中于波浪谱中心位置处。因此，本文采用波浪谱重心处的频率作为特征频率ˉ，其计算如下：

式中m0和m1分别为0 阶与1 阶谱矩，即：

根据特征频率，利用色散关系，即可得到与特征频率对应的特征波长为：

这是一个超越方程，可以用MATLAB 中fsolve函数进行数值求解。

这样，可以给出判别大尺度桩柱的量化公式如下：

为了进一步说明式（18）的合理性，本文利用式（18）对文献［22⁃23］中桩柱的尺度进行判别，所得结果与相应的文献一致。此外，对于同一桩柱，其尺度并非是一成不变的，也会随着海况的变化而发生改变。图2给出了不同直径桩柱随风速v19.5的尺度变化情况。 其中，波高谱采用Pierson⁃Moskow⁃itz 谱［24］：

图2 桩柱尺度随风速变化情况Fig.2 Cylinder scale varying with wind speed

式中 无因次常数α=8.1×10-3，β=0.74，v19.5为海面上19.5 m 高度处的平均风速。

由图2可知，桩柱直径越大，使桩柱尺度发生改变所需的临界风速也越大，即风速的变化会影响桩柱尺度的判别。因此，在应用中，首先需要根据桩柱直径和风速对桩柱尺度进行判别，在此基础上，利用绕射理论或Morison 方程进行随机波浪力的降维模拟。

4 工程算例

4.1 数值模拟分析

为了说明本文建议的降维模拟方法的有效性，首先考虑群桩的特例，即大尺度单桩的随机波浪力场模拟，其参数取值如表1所示。

表1 模拟参数的取值Tab.1 Calculation parameters and corresponding values

利用本文所提降维模拟方法生成610 条随机波浪力代表性样本集合。图3给出了z=20 m 处随机波浪力的第400 条和第500 条代表性样本，它们具有随机波浪力的典型特征。

图3 随机波浪力代表性样本Fig.3 Representative samples of random wave force

在模拟精度方面，图4为z=12 m 处随机波浪力的自相关函数比较，图5为z=12 m 和z=15 m之间的互相关函数比较。可见，采用本文方法模拟的随机波浪力自（互）相关函数拟合效果要好于Monte Carlo 方法，初步验证了降维模拟方法的有效性。图6，7 为模拟的均值，标准差与目标值比较（以z=12 m 为例），图8为z=12 m 处随机波浪力的功率谱密度函数在有效频率0.9～2.5 rad/s 范围内与目标值的比较。其中，有效频率指提供频谱总能量99.9%的频谱范围。表2为随机波浪力样本均值误差，标准差误差及功率谱误差的具体结果，从表2可以看出，本文方法的均值、标准差及功率谱密度函数均与目标值拟合一致，进一步验证了本文方法的有效性。且其均值误差、标准差误差及功率谱误差均小于Monte Carlo 方法。在模拟效率方面，采用本文方法生成单条样本耗时6.467 s，Monte Carlo 方法为6.623 s。可见，本文方法在模拟效率方面与Monte Carlo 方法大致相当，但本文方法的模拟精度更高，体现了其优越性。

图4 自相关函数比较Fig.4 Comparison of auto⁃correlation functions

图5 互相关函数比较Fig.5 Comparison of cross⁃correlation functions

图6 均值的比较Fig.6 Comparison between simulated mean values and target value

图7 标准差的比较Fig.7 Comparison between simulated standard deviations and target value

图8 功率谱的比较Fig.8 Comparison between simulated power spectrum and target value

表2 模拟精度比较Tab.2 Comparison of simulation accuracy

4.2 大尺度桩柱计算公式对比分析

为了考察大尺度桩柱基于Morison 方程带来的误差。现以4.1 中桩柱模型和波浪力参数为例，分别基于绕射理论和Morison 方程对大尺度桩柱进行计算，并将两者的结果进行对比，分析基于Morison方程计算大尺度桩柱所带来的误差。图9为z=18 m 处基于绕射理论和Morison 方程所生成的第400 条随机波浪力代表性样本比较。

图9 z=18 m 处随机波浪力代表性样本的比较Fig.9 Comparison of representative sample of random wave force at z=18 m

从图9可知，在桩柱z=18 m 处，基于Morison方程计算得到的桩柱波浪力明显偏大。为进一步分析两种计算公式的差异性，通过对610 条随机波浪力代表性样本的统计分析，图10 给出了桩柱不同位置处基于两种计算公式得到的最大波浪力均值对比。

从图10 可以看出，在大尺度桩柱竖直方向上存在某一临界点，在该临界点以下，基于绕射理论模拟的最大波浪力均值大于基于Morison 方程模拟的结果，在该临界点以上，则小于基于Morison 方程模拟的结果。

图10 桩柱不同位置处最大波浪力均值的比较Fig.10 Comparison of maximum wave force at different po⁃sitions of the cylinder

图11 为两种计算公式得到最大波浪力均值的相对误差随海况（即特征波长Lˉ）的变化情况，其中，相对误差ε计算如下：

式中Smax为基于Morison 方程得到的最大波浪力均值，Lmax为基于绕射理论得到的最大波浪力均值。

由图11 可知，当Smax小于Lmax时，随着z的增大，相对误差ε的绝对值将随之减小，当Smax大于Lmax后，随着z的增大，相对误差ε也将随之增大。同时，随着特征波长Lˉ的增大，相对误差ε=0 时所对应的z值也将随之增大。

图11 在3 种海况下桩柱不同位置的差值比Fig.11 The difference ratio of cylinder at different positions under three sea conditions

4.3 群桩效应影响分析

为了说明大尺度群桩效应对随机波浪力场模拟的影响，现基于绕射理论对平行于波向的串列三桩（N=3×1）和其中某一单桩（N=1）分别进行随机波浪力场的降维模拟。群桩的排列及波浪的作用方向如图12 所示，其中桩柱模型和波浪参数均与上文相同，桩的中心距为8 m。

图12 串列三桩位置示意图Fig.12 Schematic diagram of three cylinder positions

图13 给出了z=12 m 处1 号桩按单桩和群桩情况下第400 条随机波浪力代表性样本的比较。从图中可见，1 号桩在群桩情况下的波浪力最大值大于单桩情况。事实上，通过对610 条代表性样本的统计分析，可得1 号桩在群桩情况下的最大波浪力均值也大于单桩情况。因此，对于1 号桩，群桩效应将会使波浪力在总体上偏大。

图13 1 号桩在单桩情况与群桩情况下随机波浪力代表性样本比较Fig.13 Comparison of representative samples of random wave force of 1st cylinder in single pile case and pile group case

同样地，对2 号和3 号桩进行类似的群桩影响分析，可知群桩效应会使2 号桩波浪力在总体上略偏大，而使3 号桩波浪力在总体上偏小。

图14 给出了在群桩情况下，桩柱z=12 m 处1 号桩与3 号桩第400 条随机波浪力代表性样本比较，显然1 号桩受到的波浪力明显大于3 号桩，且通过对610 条代表性样本的统计分析，可知1 号桩受到的最大波浪力均值比3 号桩大10.35%，该结论与文献［8］谱分析法得到的群桩系数结果相似，进一步验证了本文方法的正确性。

图14 1 号桩与3 号桩随机波浪力代表性样本比较Fig.14 Comparison of representative samples of wave force between 1st cylinder and 3rd cylinder

5 结 论

本文首先根据波高过程的功率谱密度函数，建议了不规则波浪作用下桩柱大小尺度的判别方法。然后，基于绕射理论推导了大尺度桩柱随机波浪力的功率谱密度函数，并引入随机函数的降维思想，实现了随机波浪力场的降维模拟。最后，探讨了大尺度桩柱基于Morison 方程所带来的误差，以及群桩效应对桩柱波浪力的影响。通过数值分析，可得以下结论：

（1）本文方法实现了仅用一个基本随机变量即可精细地模拟大尺度群桩随机波浪力场，极大地减少了随机变量的数量。而且，相对于传统Monte Carlo 方法，本文方法具有更高的精度，并能够与概率密度演化方法结合来实现海工结构的动力反应与可靠度精细化分析，因此具有明显的优势。

（2）大尺度桩柱在竖直方向上存在某一临界点，在该临界点以下，采用Morison 方程计算桩柱时会导致波浪力偏小，在该临界点以上，则会导致波浪力偏大，且该临界点随着特征波长Lˉ的增加而向桩顶移动。

（3）当单向不规则波浪作用于群桩结构时，沿着波浪传播方向，群桩效应将使第1 桩柱的波浪力在总体上比单桩情况偏大。后面的桩柱由于受到前面桩柱的遮蔽，其波浪力逐渐减小。