大型柱状结构物波浪爬坡的数值波浪水池模拟

， ，

(天津大学 a.水利工程仿真与安全国家重点实验室， b.建筑工程学院，天津 300072)

0 引 言

为保证海洋结构物的安全性，须研究海浪与结构物的相互作用问题。圆柱形结构物为最简单也最常见的海上结构物，国内外学者对其波浪爬坡现象进行试验研究。KRIEBEL[1]进行圆柱形结构波浪爬坡的基准试验。NIEDZWECKI等[2]采用小比例模型试验，对固定圆柱结构的波浪爬坡作用进行测量，结果发现：除了波陡很小的情况之外，线性绕射理论计算出的波浪爬坡结果小于实际值。MARTIN等[3]进行关于圆柱型结构波浪爬坡的试验，并将试验结果与不同理论的计算结果进行对比，结果发现：试验得到的爬坡结果大于理论计算值。De VOS等[4]进行关于规则波与不规则波爬坡的试验，并通过试验结果拟合出关于不规则波爬坡高度的计算公式。除此之外，国内外学者也对波浪爬坡高度进行理论计算研究。KRIEBEL[5]、ISAACSON等[6]、BUCHMANN等[7]使用二阶频域程序计算波浪爬坡高度。BUCHMANN等[8]、ISAACSON等[9]使用二阶时域计算方法计算波浪爬坡高度。但这些研究方法计算出的爬坡结果与KRIEBEL的试验结果相比吻合度不高。

近年来，随着计算流体力学(Computational Fluid Dynamics，CFD)技术的高速发展，数值波浪水池技术的研究成为当今的热点。数值波浪模拟可分为仿物理造波法(推板/摇板造波)和纯数值造波法(源项造波)。在仿物理造波方面，顾挺锋[10]以哈尔滨工程大学的海洋工程水池为原型，采用摇板造波原理模拟稳定的线性规则波、不规则波，成功模拟出消波滩上的波浪爬高以及破碎现象。封星等[11]采用推板造波原理模拟二维线性规则波。在数值造波方法方面，李宏伟[12]将摇板造波方法与源函数造波方法的二维造波进行对比，并对三维模型的数值造波进行初步探讨。季红叶[13]构建三维数值波浪水池，模拟单色波及双色波，并对MARINTEK波浪爬坡试验进行模拟。

目前，国内外学者对波浪爬坡现象数以及数值波浪水池均进行了相关研究，但对波浪模拟精度、双色波模拟等方面还需要更深入的研究。本文采用CFD方法，将数值造波理论应用于波浪爬坡的模拟中，使用数值波浪水池对MARINTEK的波浪爬坡试验进行模拟，并与试验结果以及其他方法进行对比。

1 数值波浪水池模拟基本理论

1.1 Navier-Stokes方程

数值水池模拟的基本方程为Navier-Stokes方程(简称N-S方程)。N-S方程包括连续性方程和动量方程2个部分。连续性方程为

(1)

式中：ux、uy、uz分别为x、y、z3个方向的速度。

动量方程为

(2)

式中：t为时间；ρ为流体密度；g为重力加速度；υ为动力黏性系数；p为压力。

1.2 数值造波方法

采用质量源法进行波浪模拟。质量源法是在连续性方程中引入质量源函数。源项范围取水池左侧的第1层网格，源项顶端在静水面下距静水面1个波高处。对于仅沿x方向传播的三维波浪问题，使用质量源函数进行造波，连续性方程应修改为

(3)

式中：q(x,z,t)为质量源函数，计算公式为

(4)

式中：xs为造波源在x方向的位置；qs(z,t)为造波强度。

式(4)表明：在造波源区域x=xs处造波强度为qs(z,t)，在其他位置造波强度为0。加入质量源后动量方程为

(5)

式中：q为质量源。

造波强度为

(6)

式中：u(xs,z,t)为在造波源处z高度的水质点在t时刻的水平速度；Δx为网格在x方向的长度。

1.3 数值消波方法

采用人工黏性消波法对水池进行消波处理。人工黏性消波法是仿照现实波浪水池中的消波层进行的，其主要思想是在消波区域中加入阻尼项消除反射波浪。在消波区域加入阻尼后，此处的动量方程为

(7)

式中：μ(x)为消波系数，是起点为0的单调递增函数。

本文采用线性消波系数消波，即在消波区域内，消波系数沿水池水平方向呈线性分布，如图1所示。

图1 阻尼消波示意图

消波系数消波的表达式为

(8)

式中：x0为消波区域前端起始位置；x1为消波区末端结束位置；x1-x0为消波区长度，一般取1个波长；α为消波经验系数，选取合适的值可使水池获得最好的消波效果。

1.4 数值水池模拟方法

本文采用有限体积法(Finite Volume Method，FVM)对偏微分方程进行离散，离散格式采用二阶迎风格式，压力速度耦合方法选用压力隐式分离算法(Pressure Implicit Split Operator，PISO)。对于自由液面的捕捉采用较为常用的流体体积法(Volume of Fluid Method，VOF)。

2 MARINTEK波浪爬坡试验模型

本文采用MARINTEK波浪爬坡模型试验测量圆柱周围各个点的波浪爬坡结果，并将数值水池模拟结果与其他公开发表的模拟结果进行对比分析。

模型试验的原型为半潜式平台的大尺度圆柱，服役海域水深为489 m，圆柱吃水24 m，圆柱直径为16 m。模型试验在挪威Trondheim的MARINTEK拖曳水池中进行，水池宽10 m，深10 m。根据原型和试验水深比值作为模型与实物的比例，取1∶48.93，试验圆柱直径为0.327 m，圆柱吃水0.49 m，圆柱中心在水池中距造波区10 m的位置。波面监测点位置如图2所示，图中：A1、A2、A3、A4方向分别为270°、225°、202.5°、180°；点1、2、3、4距圆柱中心的距离分别为0.17 m、0.19 m、0.26 m、0.33 m。试验模拟的波浪条件包括3组线性单色波和3组双色波，本文分别选取其中1组线性规则波和1组双色波进行模拟，模型试验的波浪状况见表1。

图2 波面监测点位置示意图

表1 模型试验波浪状况

3 波浪生成模拟

3.1 模型网格设置

根据表1的波浪参数模拟该深水波浪，在模拟时缩小原来拖曳水池的尺寸，选择数值水池的尺寸为长26 m，宽5 m，深6 m，水深为5 m。消波范围为沿水池长度20 ～26 m的范围，造波源范围为顶端距离静水面1个波高，x方向宽度为1个网格。波浪载荷计算模型如图3所示。

图3 波浪载荷计算模型

为保证模拟的准确性，本文先进行波浪模拟，然后进行波浪与结构物作用的模拟，圆柱体尺寸和位置与模型试验相同。水池左侧即造波源一端采用Symmetry边界条件，顶端采用Pressure-outlet边界条件，其余2个面采用Wall边界条件。由于采用源项造波法，模型网格中没有动网格存在，对时间步长要求较低，时间步长取0.01 s，计算3 000步，即模拟时间为30 s。

在波浪模拟时，为保证波面的捕捉，在静水面上、下分别取1个波高的范围进行网格加密，z方向其他部分采用逐渐变稀疏的方式进行划分；x方向在造波区以及波浪工作范围内采用0.06 m的网格，在消波区域采用逐渐变稀疏的网格进行划分；y方向平均分为50份，即采用0.1 m的网格进行划分。

3.2 波面升高模拟结果

选取10 m处的结果进行对比，在M1单色波工况下的模拟结果如图4所示。使用滤波程序进行统计，可以得到模拟波浪的周期为1.286 2 s，与理论值误差为0.031%；模拟计算得到的波高为0.086 8 m，与理论值误差为0.696%，M1工况的波浪模拟结果较准确。在B1双色波工况下的模拟结果如图5所示。将计算得到的波浪进行统计分析，与理论对比的结果见表2。

图4 M1工况x=10 m处波面升高结果对比 图5 B1工况x=10 m处波面升高结果对比

表2 B1工况10 m处双色波统计结果对比

由表2可以看出：跨零周期计算结果准确度较高，但是波高结果存在一定的差异。可能有以下原因：(1)B1工况网格的划分密度，对于B11单色波分量来说较疏，导致B11分量模拟结果与理论值相比有一定误差，而在双色波中这种误差被扩大了；(2)源项区域顶端距离静水面的高度应该选取1倍波高的距离，此距离偏大时会使得波高偏小，本文中双色波模拟时采用的是1/3波高统计值，可能会有一定影响。

4 波浪爬坡模拟

4.1 模型网格设置

为了更精确地计算结构物周围的压力和波面，在结构物周围进行网格划分加密，自由液面处加密情况与3.1节相同。网格划分如图6所示。

图6 波浪载荷模型结构加密网格俯视图

4.2 模拟结果对比

爬坡模拟时，分别监测A1、A2、A3、A4方向上对应的1、2、3、4点，获得4个点的液面高度以判断4个方向上的波浪爬坡情况。将模拟结果与公开发表的论文的试验值以及A、B、C、D 4种不同数值模拟方法的计算结果进行对比[14]。方法A～D 的模拟参数见表3，表中：BEF为边界元技术，自由表面使用格林方程；BER为边界元技术，使用兰金-格林方程；FD表示使用频域计算；TD表示使用时域计算；2nd表示波浪的扰动采用二阶算法计算，1st表示波高采用线性算法计算。

表3 不同模拟方法参数对比

线性规则波的爬坡采用M1波浪状况模拟计算，对于M1工况本文的数值模拟计算值与方法 A，方法 B以及试验值的对比结果如图7所示，其中横坐标表示测量点距圆柱中的距离与圆柱半径的比值r/R，纵坐标表示爬坡高度与1/2波高的比值A/(H/2)。

图7 M1工况不同方向的波浪爬高计算结果对比

由图7可以看出：4个方向中方法A与方法B的计算结果较吻合；在A1、A2方向上，本文的计算结果与方法A、方法B和试验值的结果均吻合较好；在A3方向上，试验值比数值计算值大出很多，而本文计算结果与其他数值模拟结果较为接近；在A4方向上本文的计算结果比试验值稍大，而其他数值模拟结果比试验值稍小。

观察每个方向上几个点的平均高度并分别将每种计算结果进行比较，对比结果如图8所示，可以看出：本文计算值和方法B的波面增高按照A1～A4的顺序逐渐增大；方法A中A4方向平均波面增高比A3方向接近但稍微小一些；试验值A3方向的波面增高明显大于其他几个方向，其余3个方向的波面增高按照A1、A2、A4的顺序逐渐增大，根据实际经验，A3方向这种突然增大是不合理的，因而A3方向的试验数据忽略不使用。根据检测位置的布置情况，A1与波浪传播方向相互垂直，A4与波浪传播方向相互平行，A2、A3都与波浪传播方向呈一定角度。波浪速度在4个方向的分量按照A1～A4的顺序逐渐增大，根据动量定理，波浪与结构物的相互作用效果按照A1～A4的顺序也逐渐明显，因此可以推测，波面升高平均值应该按照A1～A4的顺序逐渐增大。可见，单色波计算结果与试验及其他模拟结果相比基本吻合。

图8 M1工况不同方向平均波面高度对比结果

双色波选取B1波浪状况进行模拟计算，对于B1工况本文的数值模拟计算值与方法B、C、D和试验值的对比结果如图9所示，其中横坐标表示测量点距圆柱中心的距离与圆柱半径的比值r/R，纵坐标表示爬坡高度A。

图9 B1工况不同方向的波浪爬高结果对比

由图9可以看出：在A1方向上，本文的计算结果与其他模拟和试验的结果基本吻合；在A2和A4方向上，本文计算值大于试验值，而其他数值模拟结果要小于试验值，但本文计算结果更接近试验值；在A3方向上，没有给出试验值，本文计算结果与其他模拟结果相差较大，因此A3方向上结果应舍弃，不予考虑。

观察每个方向上几个点的平均高度并分别将每种计算结果进行比较，结果如图10所示，可以看出：除A3方向上试验值未给出外，其余几种计算结果的波面升高均按照A1～A4的顺序逐渐增大，与实际情况相符。由此可见，双色波平均波面升高计算结果略大于试验值，而B、C、D 3种模拟方法的计算结果小于试验值，但整体变化趋势基本一致。

图10 B1工况不同方向平均波面高度对比结果

5 结 论

根据MARINTEK的波浪爬坡模型试验结果，本文展开数值模拟对比工作。采用数值波浪水池技术，针对单色波和双色波2种工况，模拟了波浪和圆柱周围不同位置的波浪爬坡结果。波浪结果与理论计算值进行对比，而波浪爬坡结果与试验值和其他数值方法进行对比。得出以下结论：

(1) 本文计算得到的线性规则波的模拟周期及波高结果与理论值相比的误差均小于1%；双色波模拟的周期与理论值基本相符，波高存在一定的误差，但基本满足计算要求。

(2) 与其他已发表文献的数值模拟结果相比，本文计算得到的波浪爬坡结果略大于试验值，除了A3方向外，其他方向的计算结果比其他几种模拟方法更接近试验值。这证明本文所述数值方法的计算精度要高于上述进行对比的其他几种数值方法。

[ 1 ] KRIEBEL D. Nonlinear Wave Diffraction by Vertical Circular Cylinder Part 1:Diffraction Theory[J].Ocean Engineering,1990, 17(04): 345-377.

[ 2 ] NIEDZWECKI J M, DUGGAL S D. Wave Run-up and Forces on Cylinders in Regular and Random Waves[J]. Journal of Waterway Port Coastal and Ocean Engineering,2004, 118(06): 615-634.

[ 3 ] MARTIN A J,EASSON W J,BRUCE T. Run-up on Columns in Steep, Deep Water Regular Waves [J]. Journal of Waterway Port Coastal and Ocean Engineering, 2001, 127(01): 26-32.

[ 4 ] De VOS L, FRIGAARD P, De ROUCK J. Wave Run-up on Cylindrical and Cone Shaped Foundations for Offshore Wind Turbines[J]. Coastal Engineering, 2007, 54(01): 17-29.

[ 5 ] KRIEBEL D. Non-Linear Wave Interaction with a Vertical Circular Cylinder Part II: Wave Run-up[J]. Ocean Engineering, 1992, 19 (01): 75-99.

[ 6 ] ISAACSON M D S Q, CHEUNG K. Correction Factors for Non-Linear Run-up and Wave Force on a Large Cylinder[J]. Canadian Journal of Civil Engineering, 1994, 21(05): 762-769.

[ 7 ] BUCHMANN B, FERRANT P, SKOURUP J. Run-up on a Body in Waves and Current, Fully Non-Linear and Finite Order Calculations[C]//13th International Workshop on Water Waves and Floating Bodies, 1998.

[ 8 ] BUCHMANN B, SKOURUP J,CHEUNG K. Run-up on a Structure Due to Waves and Current[C]//Proceeding of the Seventh International Offshore and Polar Engineering Conference, 1997.

[ 9 ] ISAACSON M D S Q, CHEUNG K. Time Domain Solution for Wave Current Interaction with a Two Dimensional Body[J].Applied Ocean Research,1993, 15 (01): 39-52.

[10] 顾挺锋. 海洋工程水池波浪生成的数值模拟[D]. 哈尔滨: 哈尔滨工程大学,2009.

[11] 封星, 吴宛青, 吴文峰, 等.二维数值波浪水槽在FLUENT中的实现[J]. 大连海事大学学报, 2010, 36(03):94-101.

[12] 李宏伟. 数值水池造波方法研究[D]. 哈尔滨: 哈尔滨工程大学, 2009.

[13] 季红叶. 波浪对大尺度固定柱状结构物作用的数值水池模拟[D]. 天津: 天津大学, 2016.

[14] NIELSEN F G.Comparative Study on Airgap Under Floating Platforms and Run-up Along Platform Columns[J]. Journal of Marine Structures, 2003(16): 97-134.