基于高阶Rankine源法的波物相互作用完全非线性模拟

2018-11-01，，，，，

中国海洋平台 2018年5期

，，，，，

(1. 大连理工大学船舶工程学院，辽宁大连， 116024； 2. 大连海洋大学海洋与土木工程学院，辽宁大连， 116023)

0 引言

海洋结构物是海洋资源开发长期作业的依托平台，开发过程中往往伴随着恶劣的海洋环境条件，引起结构物的不利荷载和剧烈运动[1]，在极端情况下，甚至会引起结构物整体倾覆，导致重大的工程事故。但是，传统水动力软件如AQWA[2]、SESAM等都采用线性波浪假定，忽略了波浪的非线性影响，无法满足极限波浪情况下海洋结构的设计需要。因此，有必要对非线性波浪与结构物的作用问题做进一步研究，在设计阶段准确地预测出海洋结构物在服役期间可能受到海洋环境的最大冲击，避免海洋工程事故的发生。

目前，关于完全非线性波浪与结构物相互作用的研究有2种不同的研究方向：一种方法统称数值波浪水槽技术[3]，模拟物理试验水槽技术建立数值波浪水槽，结构物放置于数值水槽中央，两端使用侧壁，前后分别设置造波与消波边界，模拟整个造波、波物作用与消波的过程。另一种方法则基于入散射势分离技术开发的开敞水域技术[4]，将速度势分离成入射势和散射势，对散射势建立积分方程，在一个开敞圆域内建立波浪与结构物相互作用的模型。由于开敞水域技术更符合工程实际情况，且具有便于造波、消波以及计算量少、计算效率高等优点，逐渐被学者研究和采用。2013年，ZHOU等[4]采用开敞水域完全非线性技术对漂浮圆柱体的绕射、受迫振动以及Ringing现象进行深入的研究，揭示非线性作用的机理。2014年，BAI等[5]对完全浸没水中的自由及系泊的柱状结构物进行完全非线性模拟。2015年，FENG等[6]利用完全非线性数值模型对并列双箱在波浪作用下的水动力共振问题进行研究。2016年，HANNAN等[7]利用完全非线性水波理论对完全浸没水中的起重船的有效工作载荷进行数值模拟分析。但是，上述研究仅限于考虑波浪与圆柱以及方箱等简单几何形状的物体作用。这主要是由于完全非线性波物作用理论要求自由水面与物面条件均在瞬时满足，计算中需要在每一时刻对自由水面以及物面进行更新。当考虑波浪对复杂曲面结构的作用时，复杂几何形状的数值计算极易导致数值不稳定的情况发生。因此，建立精确自由水面及物面网格的更新方法，准确模拟复杂曲面波物运动过程，是完全非线性波物作用模型应用于工业设计所亟需解决的关键性问题。

本文采用开敞水域数值模拟技术以及高阶边界元法[8]，在开敞圆域内建立一个完全非线性波浪与海洋结构物作用的三维数值分析模型。通过完全非线性波浪对直立圆柱作用的高阶波浪力进行模拟分析，验证本文模型的正确性。为验证本文模型对复杂曲面问题的适用性，对非线性波浪与Wigley船的相互作用问题进行模拟，并通过与JOURNEE[9]的物理模型试验对比，验证模拟方法的可靠性。

1 控制方程与边界条件

图1 计算模型示例

波浪与结构物作用问题的计算模型如图1所示(图中：SD为水底固面；SB为物面边界；SF为自由水面)。定义2组坐标系O-xyz和O-x′y′z′分别描述整个流体域和物体的运动，法向量定义为出流体方向为正。

在势流理论的假定下，以三维拉普拉斯方程[10]为基本控制方程，在时域内建立完全非线性波浪与浮体间作用的数学模型，并且瞬时自由水面条件和物面条件均是非线性的。根据入散射势分离技术，将总速度势和波面(φ，η)分解成入射势(φI，ηI)和散射势(φS，ηS)，代入控制方程和边界条件并整理，得

2φS=0 (在计算域D中) (1)

式中：U为结构物的平动速度；Ω为结构物的转动速度；t为时间；n为法向量；rb为任一点到转动中心的距离。

采用五阶Stokes波理论计算入射势的有关项，并乘以缓冲函数来避免初始效应。为保证辐射波浪向外传播而不反射回来，在计算域外围布置一层人工阻尼层。

2 模型建立

2.1 边界积分方程建立

采用Rankine源格林函数[11]，在整个计算域内应用格林第二定理，将散射势φS的边值问题转化为边界积分方程问题：

(7)

式中：c为固角系数；G为格林函数。

通过高阶边界元法将积分方程进行离散，可得

(ξ,ζ)|(8)

式中：Hk为权系数；N为离散节点总数；m为任一离散单元；M为离散单元总数；K为任一单元上的节点总数；j为离散单元上节点；J为雅克比函数；ξ、ζ分别为任一节点在参数坐标系下的横、纵坐标值。

其中，Rankine源格林函数为

2.2 复杂船体曲面的网格更新

如何准确地追踪自由水面上流体质点和复杂船体曲面上湿表面的变化情况，并根据相应的变化对水面和物面上的网格节点进行更新与重构，是完全非线性波浪与三维船体相互作用模拟的重点也是难点。

本文采用混合欧拉-拉格朗日法(Mixed Euler-Lagrange method，MEL)追踪自由水面流体质点，采用跟随式网格拓扑技术保证网格单元的高质量。自由水面节点位置可以用水平坐标(x,y)和垂向坐标z(η)来表示。在每一时间步内，首先通过MEL对各节点位置坐标进行更新。为保持在水平投影面上各节点与运动物体的水平相对位置不变，通过初始时刻的网格节点水平坐标和刚体水平位移得到新的网格节点水平坐标

Δxi(i=1,2)(10)

最后，通过插值技术，根据新的节点水平坐标，在最开始运用MEL得到的水面节点坐标中，插值得到新节点处的波面高程和速度势。首先找到距离目标点最近的节点，在其所在的每个单元中有

(11)

式中：h为形函数。

因为h只与参数坐标(ξ,η)有关，将上式改写为

(12)

采用平面二分法插值可以求得参数坐标(ξ,η)，进而通过型函数得到新节点的波面高程和速度势，得到新的相对均匀的自由水面单元节点参数。处于自由水面和物面相交的水线上的节点具有特殊性(同属于水面和物面)。为了保证经过自由水面条件更新后的水线节点仍处于物面上，需要对其进行特殊处理。

图2 水线节点处理示例

图3 船体表面网格重构示例

图2为空间坐标系和随体坐标系下的船体横剖面图。点P0(x0,y0,z0)位于水线上，经过自由水面条件更新后，与船体表面分离，需要将其移至船体表面上的目标点P(x,y,z)处。首先，通过坐标转换矩阵得到点P0在随体坐标系中对应的坐标点P0′(x0′,y0′,z0′)；然后，根据船体横剖曲线关系(Wigley船体函数关系或者船体横剖线图)将P0′强行移动到船体表面上，得到点P1′(x0′,y1′,z0′)；最后，运用转换矩阵，得到空间坐标系中船体表面上的点P1(x1,y1,z1)。由图2可以看出：点P1与目标点P的位置非常接近。因为模拟中船体在每一时间步的位移值极小，可以近似认为点P1即为交线上的点。

船体湿表面形状不断变化，物面网格也需要重构。船体曲面形状的准确捕捉和还原对船体在非线性波浪作用下的受力和运动响应的模拟准确性非常重要。在随体动坐标系下，假设物面网格节点的横坐标x在整个模拟过程中保持不变，改变各节点的y轴和z轴坐标值。如图3所示，首先根据水线节点的改变量，通过弹簧近似法[4]对物面各节点的z坐标进行处理：

Δzn+1=∑kijΔzn/∑kij(13)

式中：Δz为各节点z坐标的改变量；kij为弹簧刚度系数；n为迭代次数。

最后，由节点的x坐标以及移动后的z坐标，根据船体表面的曲面函数或者通过横剖线图进行插值，即可得到移动后节点的y坐标值。上述物面网格重构方法对实际曲面具有较高的还原性和通用性，适用于各类复杂曲面海洋结构物。

综上所述，经过自由水面网格的更新、水线节点的处理以及船体表面网格的重构等3个处理过程，就可以实现捕捉瞬时液面和瞬时复杂船体湿表面的目标，进而实现对整个完全非线性过程的时域模拟。整个处理过程具有较强的适用性，理论上可以处理实际海洋工程中大多数具有复杂曲面形状的结构物。

2.3 ∂φ/∂n奇异性处理

当计算域边界表面(如非线性自由水面)变化较大时，格林函数的导数∂φ/∂n会产生柯西主值积分奇异性问题，对完全非线性的时域模拟造成困难。∂φ/∂n奇异性的处理方法有2种：一种是BAI等[5]采用的间接计算法，方法简单但是需要增加额外的计算控制面，增加了数值模拟的计算量与时间；另一种是NING等[12]采用的直接计算法，可以得到一定精确度的奇异积分值，但是当源点位于单元边界中点时效果不是很好。

本文在第2种方法的基础上进行适当改进，通过单元细分、积分域分解、极坐标变换以及级数展开等过程将奇异积分项逐步转化为常规积分，实现奇异项的消除。包含格林函数导数部分的积分在转换后的局部坐标系中可以表示为

式中：φ(x)为速度势；Tp为参数坐标系下的实际单元的映射；h为形函数，代表单元节点对单元内任一点参数值的贡献。

图4 参数空间中的极坐标系

当源点与j节点重合时，式(14)中第j项的积分即为奇异积分。如图4所示，根据与源点重合的节点位置不同，可以将单元节点分为角点和边中点等2种情况。以第1种情况为例，采用直接法求解∂φ/∂n奇异积分。首先将参数单元分为4个小正方形单元，在小单元2、3、4中无奇异性影响，仍采用常规高斯积分。在1号小单元中，可以看作奇异性存在于α-β之间的区域内，在这一区域内积分可以表示为

ξdζ(15)

在剩下的区域中积分为

ξdζ(16)

式中：Tp2为单元区域减去奇异性区域后剩下的区域。

将实际区域α-β看作平面，在源点附近可以对空间向量进行泰勒展开：

代入式(16)，可以得到

(18)

式中：θ为极坐标系下角度的参数；f(θ)为表达式的函数；R为单元边界上点距源点的距离；α为奇异区域外缘距源点的距离；ρ为极坐标系下径向的参数。其中：

积分IS2中包含2个同样的奇异积分可以相互抵消，采用常规高斯积分即可。

将式(17)代入式(15)，可以得到

式中：ε为奇异区域边界半径；A为定义的偏导数参量Ai(θ)=cos(θ)∂xi/∂ζ+sin(θ)∂xi/∂ζ。

式(20)等号右侧第1项为无奇异项，可以通过常规高斯积分计算。式(20)等号右侧第2项积分随ε趋于0而趋于无穷，但当源点周围所有对应积分相加后总和必然等于0，故可以不考虑。

通过上述单元细分、积分方程分解、极坐标变换以及空间向量级数展开等处理过程，可以将包含格林函数导数积分中的奇异积分部分消除。在保证计算精度的前提下，只需要额外增加很小的计算量，就可以得到较高的数值模拟效率，有助于对非线性波浪与复杂曲面船体相互作用的模拟模型的建立。

3 数值结果

3.1 完全非线性波浪对直立圆柱的绕射问题

图5 直立圆柱和自由水面初始网格划分示例

本文通过与HUSEBY等[13]的试验结果及ZHOU等[4]的数值结果进行对比，验证计算模型中波浪力计算结果的准确性。采用与试验相同的模型进行分析，如图5所示，直立贯底圆柱半径r=0.03 m，计算域水深为d=0.6 m，水域半径取2倍波长。根据对称性，在1/2物面和自由水面上划分网格，其中物面上剖分12×10个单元，在水面上剖分12×20个单元。入射势采用5阶Stokes波理论解，波数k=8.1，与试验值相同，仅考虑波幅A的变化对圆柱受力的影响情况。将数值试验结果通过傅里叶变换后与试验结果进行比较。各阶波浪力对比如图6所示。

图6 各阶波浪力对比

从图6可以看出：本文数值模型计算结果大体上与HUSEBY等[13]的试验结果及ZHOU等[4]的数值结果吻合，证明了数值模型的准确性。在未损失太多计算精度的前提下可以模拟更大波幅波浪的入射情况。其中波浪一阶力的计算结果略小于其他2种结果，但相差很小，在可接受的误差范围内。随着振幅增加波浪力变化平缓，变化趋势与试验结果相同。波浪二阶力作为最主要的高阶力，在一定波高范围内，本文的计算结果与试验吻合良好。本文波浪三阶力和四阶力的结果随着波幅A变化较大；本文五阶波浪力的计算结果小于其他2种结果，但差距不大。

3.2 完全非线性波浪对Wigley船的作用问题

为了验证本文数值模型对凸出水线型复杂曲面海洋结构物在完全非线性波浪作用下的模拟能力，对零航速Wigley船在非线性波浪作用下受到的波浪力和运动响应进行模拟，并将结果与JOURNEE[9]的试验结果进行对比。本文数值模拟所采用的Wigley船型的数学表达式为

图7 Wigley船及自由水面网格划分

船长L=3 m，船宽B=0.3 m，吃水T=0.187 5 m，船排水量0.078 m3，重心距船底0.17 m，纵向惯性半径Ryy=0.75。入射波采用五阶Stokes波，入射角180°，波幅A=0.01 m，波长λ/L取值0.5～2.5。根据对称性，在1/2 Wigley壳体上划分8×16个单元。计算域半径取为R=5L，在1/2水面上划分16×30个单元。Wigley船及自由水面网格划分如图7所示。

图8为Wigley船受到的波浪力和力矩随无因次波长λ/L的变化情况。由图8 a)可以看出：本文关于Wigley船的纵向波浪力的数值结果与试验值大部分吻合良好，只有在波长较小时，两者存在一些偏差。由图8 b)可以看出：本文关于Wigley船垂向波浪力的模拟计算结果与试验值存在一定误差，波长较大时数值结果略大于试验值。由图8 c)可以看出：绕y轴纵摇力矩的计算结果在波长较小时与试验结果吻合良好，在大波长时略小于试验值。

图8 Wigley船波浪力和波浪力矩对比

图9为Wigley船在波浪作用下的位移响应与试验结果对比，可以看出：Wigley船的垂向位移在波长较小时比试验值偏小且没有反应出试验值的波动性，可能是模型试验精度误差或者表面张力黏性引起的；而在波长λ/L大于1.0时，数值计算结果偏大一些，但是误差不大。本文的纵摇计算结果与试验结果基本吻合，除个别点外偏差不大，具有一定的可靠性。因此，本模型对具有非规则表面的Wigley船的时域模拟基本与试验值吻合良好，进而证明了本模型的处理方法对瞬时自由水面和船体复杂湿表面的跟踪与还原的准确性，具有较高的实用价值。