一种适用于陡波条件下研究半潜平台周围波浪非线性爬升效应的计算方法

单铁兵，杨建民，李 欣，肖龙飞

（上海交通大学 海洋工程国家重点实验室，上海200240）

1 引 言

随着深海油气田的不断开发，浅海开发技术已难以满足要求，传统的重力式平台和导管架平台因其自重和造价随水深增加而大幅度增加，已经不适应深海油气发展的要求，取而代之的是浮式生产系统。其中半潜平台具有性能优良，运动响应较小，工作水深范围广，抗风浪能力强，甲板面积大及装载能力强等优点，在海洋石油勘探和开采中得到了广泛的应用[1]。随着石油勘探开采朝深海以及超深海发展，复杂的深海环境条件以及近年来超强台风等极具破坏力的极端海况的频繁出现，对平台的安全性提出了严峻挑战。

诸多的安全性评估指标中，气隙预报（海洋平台下层甲板底部至波面间垂直距离的预报）一直是半潜平台设计过程中一个关键问题[2]。在波浪与平台的相互作用过程中，有一种现象特别引起关注，即波浪遇平台立柱后将沿其表面往上攀升，有时其高度明显大于其它区域的波浪，甚至能爬升至平台下甲板，以致严重影响平台的气隙性能，该非线性现象称为波浪爬升。波浪的爬升过程从理论上可解释为当入射波浪冲击浮体时，水波的动能转化为势能同时波高发生放大的过程。

导管架平台由于相应管架直径小，对入射波扰动小，故波浪作用机理简单，采用线性理论即可对其周围的波面分布进行预报。对于大型的浮式结构物如半潜平台而言，因体型大，且下浮体与甲板之间通过若干立柱支撑，当波陡较大时，波浪非线性效应较强，柱体周围波浪爬升现象尤为明显，波高易发生放大，并表现出较强的非线性特征，甚至超过平台下甲板发生越浪，此过程中产生的砰击载荷将对平台的安全性构成威胁。陡波条件下的波浪爬升问题近年来已成为水动力学研究的热点问题之一，是平台气隙预报的一个重要方面。

很早以来就有学者对波浪的绕射和爬升效应进行了试探性的研究。McCamy等[3]1954年给出了单根圆柱周围绕射波浪的经典线性解。Niedzwecki等[4]1992年针对该圆柱上的波浪爬升问题开展了模型试验，试验结果表明，波陡（H/λ，其中H为波高，λ为波长）较大时，采用线性散射理论大大低估了波浪的爬升高度。Huston等[5]对一座四立柱的张力腿平台（TLP）在单色规则波下进行了一系列的试验测量，旨在研究波陡、立柱间的间隙等对立柱上的波浪爬坡和甲板下部波浪放大的影响。Lee和Newman等[6]1994年采用基于二阶边界元的方法对平台周围的波浪效应进行了模拟。尽管二阶理论能够提高波浪爬升预报的精确度，但仍不能满足平台详细设计的要求，且无法捕捉到波浪冲击立柱时的高阶非线性现象[7-8]。

计算机的发展使平台立柱周围的波浪散射、海浪沿立柱爬升、甚至砰击平台甲板后波浪发生变形和翻滚等强非线性现象的直接数值模拟成为可能。与二阶理论不同的是，CFD理论基于N-S方程直接求解时域内的流场。近年来有学者针对一些自由液面大变形的强非线性问题提出了较为精确的追踪重构自由面的数值方法。其中常见的有Volume-of-Fluid（VOF）、Smooth-Particle-Hydrodynamics（SPH）、Level-Set（LS）以及CIP方法等。Donald[9]采用两种方法求解重力式平台周围的波面分布，第一种为二阶散射程序（WIMIT），基于势流和小波陡假设,采用截断的摄动展开并结合边界元理论求解波动场，考虑二阶非线性特征。第二种为基于VOF理论的CFD方法，该方法直接求解不可压缩N-S方程，采用的是全非线性自由面边界条件。通过对比发现，线性波条件下，二阶散射理论可较为准确地模拟平台上的波浪爬升，但当波陡较大时，二阶摄动理论限制了其预报的有效性。相比而言，VOF理论即使在入射波浪波陡较大时也能够捕捉到波物相互作用时的高阶非线性效应。Roberto等[10]基于雷诺平均的N-S方程，使用仅求解水相的LS方法模拟自由液面，研究粘性效应对立柱周围波面分布的影响。

本文以连续性方程和动量方程为控制方程，按照与深水试验池相同的摇板造波方式，采用VOF方法捕捉自由液面，建立了三维数值波浪水槽，为了防止水池末端的波浪发生反射而对工作区域造成影响，在动量方程中添加源项进行消波。采用六面体结构化网格对数值模型进行离散，保证了网格的质量，提高了计算的精度。文中还对多种网格划分形式进行了详细研究，同时与试验测量值进行对比，确定了最优的网格划分方案，保证了数值方法的可行性。实际的模型试验中，因浪高仪测量设备数量有限，模型试验前确定平台哪些位置安装浪高仪对研究平台周围的波浪分布更具价值就显得十分重要。此外，入射波陡较大时，波浪的非线性特征较强，与平台之间的水动力作用加剧，将导致平台立柱周围的波浪爬升效应变得更为复杂，甚至产生砰击、越浪等对平台结构和安全性构成一定威胁的强非线性现象，相比平缓的波浪，研究波陡较大的波浪更有实际意义。介此，在该波浪水槽中，对陡波条件下半潜平台周围的波浪分布进行数值模拟。计算时，在平台周围设定多个波浪监测点，用以详细预报不同位置处波浪的绕射和爬升效应，确定波浪非线性较强的区域，为下一步模型试验中确定浪高仪的安装位置提供重要依据。此外，文中还从物理现象上模拟并分析了波浪沿立柱的绕射和爬升过程，为今后进一步研究波物之间的非线性相互作用，如波浪参数对立柱周围波浪爬升的影响、极限波浪下半潜平台的强非线性砰击和越浪、半潜平台的气隙响应等提供参考。

2 数值波浪

2.1 控制方程

半潜平台波浪绕射和爬升的数值模拟是在一个类似于海洋深水试验池的数值水池中进行的。数值模拟包括摇板造波的动态区域，过渡区域，工作区域以及消波区域。

设定流体不可压，则流场的连续性方程为：

相应的动量方程：

式中：x1，x2，x3分别为x，y，z方向的坐标；u1，u2，u3分别为x，y，z方向的速度；f1=0，f2=0，f3=-g；S1，S2，S3分别为x，y，z方向上的附加源项；μ为动力学粘性系数；ρ为流体的密度。

自由液面采用VOF（Volume of Fluid）方法进行求解。该方法的基本原理是通过网格单元中的各流体和网格体积比函数aq来求解自由面的位置。aq=1为该网格单元全部由指定相流体所占据；aq=0说明该网格单元内无指定相流体；0＜aq＜1说明该网格单元内有部分指定相流体。气液两相流的自由面波动的求解方程可写为：

其中：q=1表示空气相，q=2表示水相。

2.2 摇板造波

数值造波和消波的技术是数值波浪水池模拟成功的关键，一直是国内外学者关注的问题。造波的方式有很多种，Boo等[11]通过在入口边界给定速度或速度势，然后基于时间步进方式造出二阶和三阶stokes波。Brorsen和Larsen等[12]提出了一种适合积分方程法的域内源造波方法，即在计算区域内沿水深方向设定一造波源，在源项两边同时产生方向相反的两列波。周勤俊和王本龙等[13]从忽略粘性的欧拉方程出发，以N-S方程为控制方程，将源项添加到动量方程进行造波和消波。董志等[14]采用仿物理的推板造波方式，通过设定造波板边界按简谐运动实现规则波的模拟，该方法造波原理简单，物理造波机理论能较容易地应用到数值造波中。

实际的深水试验池大多采用摇板式造波方法生成波浪，即通过机械驱动使摇板绕固定轴摆动，强迫水质点运动，达到模拟相应波动的目的。波浪波高通过摇板的摆动幅度控制，其摆动频率则决定波浪的周期大小。实验室对规则波的模拟中，摇板驱动信号通过线性波浪理论计算得出。根据波浪理论，若摇板按固定角度做简谐运动，在离开摇板一段距离后波浪将呈现二维规则波特性，即波浪周期与摇板简谐运动相同，波高与摆幅有关。本文基于CFD理论，采用与上海交通大学深水试验池相同的摇板造波方式生成波浪。

对于规则波浪，造波板自由液面处的运动方程可写为：

其中：S0为造波板自由液面处的运动幅值，k为波数。

通过势流理论和造波板的运动方程，可推导出造波板运动幅值和波高H之间的转换关系[15-16]，如（5）式。已知规则波浪的相关参数，为了确定造波板的运动规律，需要知道摇板造波方式下的水力传递函数T（ω,d），该值同样可由势流理论得出，如（6）式所示。

其中：d为水深，k为波数，ω为波浪圆频率。而k可由下面的色散关系求得。

造波板的运动方程可写为：

根据线性造波机理论，在水深为d的水池中，波面方程为：

本文利用动网格模块中的铺层（Laying）功能中的网格更新技术，指定造波板为动边界，并绕某一固定点转动来实现数值造波。按照深水试验池中所采用的线性造波机原理，推导出摇板摆动的角度信号，如（10）式所示。其中R为初始时刻造波板的水下长度。

2.3 源项消波

对于数值波浪水池，消波功能的实现有很多。如主动式消波、阻尼层消波、Sommerfeld辐射边界消波、动量源项消波和多孔介质消波等。主动式消波[17]是指在造波的同时，通过实时的修正造波板的运动，在生成目标波浪的同时产生额外的波动抵消反射波浪。而被动式消波方法如动量源项消波，通过在动量方程中添加源项以消除反射波浪。海绵层阻尼消波[17]主要是通过在特定的计算区域设定1-2倍波长的人工阻尼区，对波浪垂向速度进行强迫衰减，同时开边界处满足Sommerfeld辐射条件。多孔介质消波方法[19]是一种仿物理消波法，通过在N-S方程中增加动量衰减的源项达到消波的目的。

模拟计算过程中，采用王本龙等[13]提出的消波技术，在消波区域内添加源项来消除水池尾部池壁的波浪反射，确保工作区域波浪不受干扰。消波段内波动场可写为：

i=1方向，已加源项的动量方程按照忽略粘性的欧拉方程出发可离散为：

未添加源项的动量方程可离散成：

将（11）式代入方程联立求解可得源项S1的表达式：

同理i=3方向的源项S3的表达式：

其中下标C代表计算值；上标N+1、N分别代表N+1、N时刻的值；λ为与空间位置有关的光滑过渡加权函数，避免速度的剧烈变化影响流场的求解精度。λ=λ（x）在数值水池消波段的表达式为：

式中：xb为消波区的起始位置，xe为消波区的末端，因此可知：

2.4 几何模型的建立

本文采用ICEMCFD软件建立了波浪模拟所需的几何模型，其示意图如图1所示，相应的三维立体模型见图8和图9。设定坐标原点位于平台重心在初始自由液面的投影处。x轴正向与波浪传播方向一致，y轴正方向与平台横向方向一致，z轴沿平台的吃水方向。该数值波浪水池长24 m，宽5 m，水深8 m，静水面以上的高度为1 m，造波板位于X轴-10 m位置处，各区域长分别为摇板及前端过渡区5 m，工作区10 m，尾部过渡区5 m，尾部消波区4 m，在水池试验区域的中央位置设置一虚拟浪高仪来获取波浪的时历。

图1 数值波浪水池的示意图Fig.1 Schematic show of numerical wave tank

图2 电容式浪高仪测量系统Fig.2 Measurement system of capacitive wave probes

2.5 波浪时历的试验测量

为了验证数值波浪水池的精确度，在上海交通大学海洋深水池内对波浪时历进行了实际测量。用于该试验研究的测量仪器主要为电容式精密浪高仪，该仪器由探头、弓形支架、电缆以及信号输出系统等组成，如图2所示，其通常安装在水池拖车下方，通过测量因水位引起的电容变化来获取波浪时历。试验前，首先需对浪高仪进行标定，以确定电压转换系数。待水体充分平静后，通过数据采集系统对浪高仪进行采零处理，正式采集得到的波浪时历均减去该采零值。随后，在控制程序中输入波浪相关参数并转化为摇板信号来驱动造波板实现物理造波。待波形稳定之后进行数据采集，该次试验的采样频率选定为40 Hz，采样时间3分钟，采样点数大约为7 200点。

2.6 网格尺度的测试

CFD计算中，网格的质量、尺度和布置方式等将直接影响数值计算结果的精度。网格质量较差势必影响数值结果的收敛性，甚至可能造成解的发散。六面体结构网格具有存储空间小、数据结构简单、迭代运算速度快、网格质量高的优点，且对于VOF模型来说，排列整齐，形状规则的六面体网格不仅可以提高求解精度，同时可以捕捉到更为光顺的自由液面。因此，本文全部采用H-H型六面体结构化网格对数值波浪水槽进行离散，以确保网格质量良好且方便布置网格节点。此外，过密的网格将大大增加迭代计算所需的时间，过于稀疏的网格不仅无法捕捉到流场的相关细节，而且将引起严重的数值耗散，以致计算结果失真。一般而言，在进行网格的生成时需遵循的原则是：在不影响计算精度的情况下，尽量减少网格数量，提高计算的速度。因此本文在网格划分时对一些流场梯度变化较大的区域进行了加密，而对物理量不敏感的区域采用了粗网格，如自由液面附近采用较细的网格，远离这些区域网格间距按指数分布逐渐拉大，以起到疏密有致的效果。但即便如此，计算精度也不能完全得到保证，因此必须同试验测量值进行对比，才能确定最终的网格划分方式。

结合试验室的造波能力，本文选取了陡度较大的规则波浪作为研究的对象，即波陡（H/λ，其中λ为入射波波长）为1/15，周期T为1.0 s，波高H为0.104 m，采用三种密度由疏到密的网格划分方式对该波浪进行了数值模拟。考虑到AB段网格对数值造波是否能够成功有重要的影响，为了减少波浪传播过程中的数值耗散，对其等距离布置较细密的网格，并按一个波长内100个、50个和25个节点三种方式来检验网格的不同间距对计算结果带来的影响。此外，考虑到自由液面附近，相应的物理量梯度变化大，在EF段同样采用等距离的网格布置方式。且三种方案中，单个波高内的网格数分别为40、20和10个。BC段为模型区域，该段的网格参数无论对计算迭代的收敛时间和精度，还是对模型周围的波动场细节的模拟等都起着至关重要的作用，因此该处的网格应更为细密。其余段的网格划分均按指数分布的形式，尺度由梯度变化大的区域往梯度小的区域逐渐拉大，具体如表1所示。

表1 网格尺度的测试Tab.1 Grid size test

从图4可以看出，网格C的计算结果与试验值偏差严重，且随着模拟时间的增加，波形因数值的耗散逐渐衰减，因此不符合高品质波浪的要求。而网格A和网格B所得到的结果均与试验值较为接近，说明网格B的划分方案可以使模拟结果达到要求。相比网格B，划分最为精细的网格A并没有带来更好的精度，且由于网格数量巨大，将大大地消耗数值运算的时间和资源。因此，从计算的精确度和时间消耗两方面权衡考虑，最终我们选取网格B作为数值波浪水池网格划分的基础。

图5为网格B方案下，数值波浪水池中央位置处未经扰动的波形时历曲线，造波板周围的波形分布如图7。通过与线性Airy波理论值对比发现，计算得到的波峰高而尖，波谷较为平坦，波浪出现了非线性的特征，已经不再满足线性波理论。

图3 数值波浪水池网格划分方案的示意图Fig.3 Schematic show of mesh strategy of numerical wave tank

图4 不同网格划分方案的波浪时历分布Fig.4 Time series of wave elevation model in different mesh scheme

图5 试验值与线性波理论值对比Fig.5 Comparison of wave elevation between test and linear wave theory

图6 计算值与线性波理论值的波浪谱密度对比Fig.6 Comparison of wave spectrum between calculation and linear theory

图7 网格B方案下摇板附近的瞬时波形图Fig.7 Instantaneous wave shape around flap wavemaker in MESH-B

采用FFT理论可将一段时历的波列转化为基于波频的谱密度分布形式，从而得到波浪的各阶谐频成份以及相应的贡献[22]。从图6中可以看出，airy波为线性规则波，其仅包含一阶谐频成份，由于入射波周期T为1.0 s，故谱峰值出现在频率ω约为6.28 rad/s处。此外，我们对网格B划分方案计算得到的时历曲线进行频谱分析发现，谱密度曲线中出现两个峰值：第一个峰值出现在频率ω约6.28 rad/s位置，即主频位置，表示一阶谐频成份（First-harmonic component）；第二个谱峰对应的频率约为12.56 rad/s，其值的大小恰好为主频的2倍，对应的谱峰值称为二阶谐频成份（Second-harmonic component）。故认为该入射波具有较明显的二阶非线性特征。尽管波浪的一阶线性成份所占的比例较大，其高阶谐频成份却是影响波浪非线性爬升效应的重要因素，在以后的研究中将对此重点考察。

3 波浪绕射和爬升的数值计算

3.1 模型建立

本文在上述理论建立的数值波浪水池内，以中海油3 000 m深水半潜式钻井平台为研究对象，对其周围的波浪绕射和波浪爬升现象进行了数值模拟。文章重点研究波浪散射效应（包括绕射、反射等效应）对平台周围波浪爬升的影响，因此设定平台为刚性固定，以消除运动引起的辐射效应的干扰。

表2 中海油半潜平台主要参数Tab.2 Main pariticulars of CNOOC’semisubmersible platform

首先建立了计算所需的几何模型，如图8和9所示。该半潜平台包含水下浮箱、四根直立柱、横撑和平台甲板等，模型缩尺比为1:60，平台实型和模型的具体尺度见表2。半潜平台置于数值波浪水池工作区域的中央位置，平台吃水为34.5 cm，初始气隙（静水面至平台下甲板的垂直距离）为15.6 cm。

为了测量平台周围波面分布情况，在平台前立柱周围设置5排浪高仪，各排浪高仪间的角度为45°，后立柱周围设有1排浪高仪。为了揭示波浪沿立柱爬升的过程，每排均设置四根浪高仪，其总长（R）为立柱的等效半径。考虑到波浪越靠近立柱，波浪爬升效应越复杂，故设置浪高仪沿径向的布置间距（r）由外到内逐渐变小，分别为0.4r/R，0.4r/R，0.193 75r/R和0.006 25r/R，且单排内浪高仪的编号依次为iA，iB，iC和iD（i为浪高仪的排号），如图10所示。

3.2 边界条件

图8 半潜平台三维数值模型Fig.8 3D numerical model of semi-submersible

图9 计算区域的几何模型Fig.9 Geometrical model of calculated domains

如图9所示，水池顶部设置为压力入口边界条件，半潜平台和造波板为无滑移壁面。因数值计算中，所选的波浪为迎浪状态（与平台艏部方向夹角为180°），且平台沿中纵剖面对称，因此水池内流场亦关于中纵剖面对称，为了节省网格的数量，设定水池中纵剖面为对称面，其余均默认为壁面。波动场采用层流模式（Laminar）求解，通过VOF方法追踪自由面；压力—速度耦合项采用SIMPLEC算法进行迭代计算，动量方程中的对流项和扩散项均采用二阶迎风格式；VOF方程采用几何重构法[20]（Geo-Reconstruct）求解。在数值模拟时，设定左边界上部的摇板为运动边界，使摇板绕其末端转动，采用铺层技术控制动态造波区域的网格变形、溃灭和再生来产生波浪。

图10 虚拟浪高仪的布置位置Fig.10 Layout of virtual wave probes

3.3 网格划分

图11 计算区域网格分布Fig.11 Meshes of calculated domains

图12 造波板和自由液面附近网格划分Fig.12 Meshes around wavemaker and free surface

图13 平台周围的网格划分Fig.13 Meshes around semi-submersible

图14 立柱周围的网格划分Fig.14 Meshes around columns

CFD计算中，网格离散化的质量将对数值计算结果的精度产生重要影响，本文采用分块耦合网格划分技术，将整个流场划分成若干区域，各流域均采用六面体结构化网格进行离散。多块结构化网格是单块结构网格和非结构网格的一个折中方案，它简单而又能处理较复杂的形状。按照上面阐述的方案生成流场网格，部分区域相应的网格分布如图11-14所示，整个模型网格单元在360万左右。

3.4 计算结果与分析

图15为Mavrakos等[21]的模型试验中，规则波浪条件下平台立柱周围的波浪爬升情形。本文也得到了类似的波面分布，如图16所示，且波浪沿立柱往上爬升过程中形成的“水舌”也极为相似，基于VOF方法求解得到的自由液面分布也较光顺，因此文中采用的数值计算方法能够较好地模拟这一非线性特征。

图15 模型试验中波浪爬升情形Fig.15 Wave run-up in model tests

图16 数值模拟中立柱周围的波浪爬升Fig.16 Wave run-up in numerical simulation

图17为前端立柱、后端立柱迎浪面中央位置，波浪沿柱体爬升时最大波面分布情况。从图中可以看出，由于正面迎浪，两处的波浪爬升都较大，且随着波浪向前传播，最大波面升高逐渐增大，在靠近立柱壁面时波浪迅速爬升至最高位置，这与Morris-Thomas[22]以及Nielsen等[23]研究结论相一致。

图18给出的是前后立柱壁面附近的波浪爬升时历分布。从图中可以看出，相比平台前立柱壁面附近的3D位置，6D处（后立柱的壁面附近）波浪时历曲线显示出更强的非线性特征，且波浪爬升幅值更大，部分波浪已经爬升至下层甲板，此时上端水体仍具有一定能量，将对甲板局部产生载荷冲击。同时结合图24（b）可知，文中采用的计算方法可以捕捉到波浪砰击后立柱表面时所产生的水体分离和破碎等强非线性现象。由于该类砰击载荷可能造成平台甲板的结构破坏，因此已成为国际上研究的热点问题之一，今后将对该强非线性现象展开详细分析。

图17 最大波面升高随r/R的变化曲线Fig.17 Maximum wave elevation along wave probe row

图18 立柱壁面附近的波浪爬升时历Fig.18 Time history of wave elevation near column wall

图19显示的是虚拟浪高仪组2（虚拟浪高仪组3顺时针方向45°位置）和4处（虚拟浪高仪组3逆时针方向45°位置）的最大波面分布情况。尽管波浪爬升幅度比柱体的迎浪面中央位置小，但仍较为明显，并且表现出较强的非线性特性。虚拟浪高仪组2处最大的波面升高规律与虚拟浪高仪组3与6处基本一致，越靠近立柱，波浪爬坡现象越明显。然而，虚拟浪高仪组4处，最大的波浪爬升并未发生在立柱壁面附近，而是离立柱较远的4B位置，造成这种现象的原因可能是由于波浪遇前立柱后发生绕射和反射，部分反射的水体与入射波在远离柱体位置发生叠加所致，该现象同时由Mavrakos等[21]在模型试验中得以证实。

此外，通过波浪时历对比发现，尽管4B与2D两位置处最大波面升高相接近，但在4B处，由于立柱内侧区域受到的干扰因素多，使得波谷位置单个周期内出现了明显的二次波峰，具有比2D处更强的高阶谐频非线性特征。

当波浪绕过虚拟浪高仪组2和4处继续往后传播时，其爬升幅值迅速减小。且一个立柱半径R范围内，最大波面升高基本无变化，但其非线性特征却越明显,如图21、22所示。

图19 最大波面升高随r/R的变化曲线Fig.19 Maximum wave elevation along wave probe row

图20 立柱壁面附近的波浪爬升时历Fig.20 Time history of wave elevation near column wall

图21 最大波面升高随r/R的变化曲线Fig.21 Maximum wave elevation along wave probe row

图22 立柱壁面附近的波浪爬升时历Fig.22 Time history of wave elevation near column wall

对比沿平台内侧布置的浪高仪组3、4、5发现，3处的波浪爬升最为显著，随着波浪绕立柱继续向前传播，波高逐渐减小，波浪爬升效应明显减弱，这与Nielsen等[23]在对半潜平台的单立柱展开模型试验时所得出的结论相同。值得注意的是，3处波浪爬升幅度虽相对较大，但波形较为规整，波谷的形状较为对称，随着波浪绕立柱往前传播至4处，波谷位置出现了二次波峰，在波浪传至5处时，二次波峰幅值不断的增大，同时主波峰（一个波浪周期内第一次出现的波峰）的形状更尖瘦，该过程将使波浪非线性效应增强，与平台间的相互作用变得更为复杂。

为了更为清晰地理解波物之间的相互作用，文中对沿平台立柱的绕射和爬升现象进行了模拟。图23为平台前立柱周围波浪的爬升过程。该情景可描述为：波浪遇立柱时，部分水流绕过立柱继续朝前传播，而另一部分水体因立柱的阻碍开始沿立柱爬升，此时波浪的能量很大。随着波浪继续往上爬升，其能量不断耗散，兴起的波峰不断增高，波峰周围的速度逐渐减小。待波峰增至一定位置后停止升高，在重力的作用下开始回落，但此时紧贴立柱表面的水体却继续往上爬升，如图（c）所示，直至其速度接近0时，波浪开始下滑，此时可以观察到，回落的水体沿立柱绕射至侧面时将兴起波浪，且其波峰周围的速度显著增大，如图（d）所示。

图24显示的是后立柱迎浪面波浪的爬升现象。从图（a）可以看出，后立柱周围的波浪已上升到平台下甲板并与之发生非线性砰击。砰击结束后，下甲板上的部分水流由于重力作用直接脱落，剩余部分沿立柱回落，同时有小部分水体与尾部水舌发生分离，如图（b），这种现象可能是因水流下滑过程中，两部分水体速度不同引起的。相比前立柱，该规则波浪下，后立柱周围的波浪爬升情况更为严重，非线性效应更强，原因可能因后立柱位于前立柱的尾流中，其周围的波面更易受扰动影响。

图23 前立柱迎浪面波浪的爬升过程Fig.23 Wave run-up process on front side of forward columns

图24 后立柱迎浪面波浪爬升现象Fig.24 Wave run-up process on front side of aft columns

4 结 论

文中基于N-S方程和VOF方法建立了三维数值波浪水槽，运用类似于深水试验池的摇板造波方式产生波浪。采用H-H型六面体结构化网格对含平台模型的数值波浪水槽进行离散，保证了网格的质量，提高了计算的精度。此外，文中还对多种网格划分形式进行了详细研究，同时与试验测量值进行对比，确定了最优的网格划分方案，并保证了数值方法的可行性。随后，在该水池中对陡波条件下半潜平台周围的波浪绕射和爬升效应进行了数值模拟。结果表明，本文采用的计算方法能够对平台周围的波面分布给出详细的数值模拟，同时可预报出波浪沿平台立柱的绕射以及从波浪开始爬升至回落的完整过程，为下一步有关平台气隙的模型试验前确定浪高仪的安装位置提供了重要参考。计算获得的波浪爬升情形和模型试验相似，且文中得到的重要结论与相关文献一致，进一步证明该方法模拟波浪绕射和爬升等现象是可行的。通过计算和分析可以得到以下结论：

（1）立柱正面迎浪位置，波浪爬升较大，且越靠近立柱，波峰越尖瘦，波浪爬升越为严重。因此在下一步的模型试验中，需在该位置附近安装浪高仪进行进一步研究。

（2）一般而言，越靠近立柱波浪爬升效应越明显，但也有例外，如偏离前立柱中纵剖面45°的内侧位置附近（浪高仪组4附近），可能的原因是由于内侧位置发生绕射和反射，部分反射的水体与入射波在远离柱体位置发生叠加所致。

（3）平台立柱的内侧，波面受干扰的因素多，波浪爬升现象比相应的外侧更为明显，爬升幅度更大，成因复杂，应加以详细研究。

（4）由于波浪沿前立柱的绕射、遇后立柱的反射以及波浪叠加等效应的相互作用，平台前立柱的背面附近以及前后立柱之间区域的波面分布将变得较为复杂，波浪的非线性效应增强，应给予足够重视。

[1]王世圣,谢 彬,曾恒一,冯 玮,李晓平,张海滨.3 000米深水半潜式钻井平台运动性能研究[J].中国海上油气,2007,19(4):277-280.

[2]Stansberg C T,Baarholm R,Kristiansen T.Extreme wave amplification and impact loads on offshore structures[C].Offshore Technology Conference,OTC 17487,2004.

[3]McCamy R,Fuchs R.Wave forces on piles:A diffraction theory[R].Tech.Memo No.69,Beach Erosion Board,US Army Corps of Engineers,Washington DC,USA,1954.

[4]Niedzwecki J N,Duggal A S.Wave run-up and forces on cylinders in regular and random waves[J].J Waterways,Port,Coastal,Ocean Eng,1992,118:615-634.

[5]Niedzwecki J M,Huston J R.Wave interaction with tension leg platforms[J].Ocean Engineering,1992,19(1):21-37.

[6]Lee C H,Newman J N.Second-order wave effects on offshore structures[C].Proceedings of the Seventh International Conference on Behaviour of Offshore structures,1994,2.

[7]Donald G Danmeier,Robert K M Seah,Timothy Finnigan,Dominique Roddier.Validation of wave run-up calculation methods for a gravity based structure[C].Proceedings of the 27th OMAE,2008.

[8]Teigen P,Trulsen K.Numerical investigation of non-linear wave effects around multiple cylinders[C].11th International Offshore and Polar Engineering Conference,2001.

[9]Danmeier G D,Robert K M S,Finnigan T,et al.Timothy Finnigan,Dominique Roddier.Validation of wave run-up calculation methods for a gravity based structure[C].Proceedings of the 27th OMAE,2008.

[10]Roberto M,Andrea D M,Riccardo B.Numerical simulation of the flow around an array of free-surface piercing cylinders in waves[C].Proceedings of the 26th OMAE,2007.

[11]Boo S Y,Kim C H.Simultation of fully nonlinear irregular waves in a 3-D numerical wave tank[C].Proceedings of the 17th ISOPE,1994.

[12]Brorsen M,Larsen.Source generation of nonlinear gravity waves with the boundary integral equation method[J].Coastal Engineering,1987,11(2):93-113.

[13]周勤俊,王本龙,兰雅梅,刘 桦.海堤越浪的数值模拟[J].力学季刊,2005,26(4):629-633.

[14]董 志,詹杰民.基于VOF方法的数值波浪水槽以及造波、消波方法研究[J].水动力研究与进展A辑,2009,24(1):15-21.

[15]孙昭晨.推摇混合式造波机理论曲线[J].港口工程,1988(4):30-32.

[16]陶建华.水波的数值模拟[M].天津:天津大学出版社,2004.

[17]李雪临,任 冰,王永学.VOF方法中主动吸收式无反射数值造波研究[C].第二十届全国水动力学研讨会文集,2006.

[18]韩 朋,任 冰,李雪临,王永学.基于VOF方法的不规则波数值波浪水槽的阻尼消波研究[J].水道港口,2009,30(1).

[19]詹杰民,董 志.黏性数值波浪水槽的多孔介质消波方法[C].第十四届中国海洋（岸）工程学术讨论会论文集（上册）,2009.

[20]Youngs D L.Time-dependent multi-material flow with large fluid distortion[C]//Proceedings of Numerical Methods for Fluid Dynamics.New York,1982.

[21]Mavrakos S A,Chatjigeorgiou I K,Grigoropoulos G,Maron A.Scale experiment of motions and wave run-up on a TLP model,subjected to monochromatic waves[C].Proceedings of the 23rd OMAE,2004.

[22]Morris-Thomas M T,Thiagarajan K P.The run-up on a cylinder in progressive surface gravity waves:Harmonic components[J].Applied Ocean Research,2004,26(3-4):98-113.

[23]Nielsen F G.Comparative study on airgap under floating platforms and run-up along platform columns[J].Marine Structures,2003,16(2):97-134.