FPSO水动力特性的完全非线性数值模拟

孙雷*1,2，陆婷婷1，邓潇潇1，刘昌凤3 1大连理工大学船舶工程学院，辽宁大连 116024 2高新技术船舶与深海开发装备协同创新中心，上海200240 3大连海洋大学海洋与土木工程学院，辽宁大连 116023摘 要:［目的］为提高海洋结构物的安全性能，针对波浪与结构物相互作用的问题开展完全非线性数值模拟研究。［方法］基于三维完全非线性时域势流理论及高阶边界元法（ HOBEM），建立波浪与结构物相互作用的开敞水域模型。采用速度势分离技术将整个问题分解为入射部分和散射部分，入射势由理论解给定。采用混合欧拉—拉格朗日（ MEL）方法追踪瞬时自由水面的流体质点，并采用四阶龙格—库塔法对瞬时自由水面进行更新。引进虚拟函数计算波浪载荷，而非直接求解速度势的时间导数。在自由水面的外侧设置人工阻尼层，防止波浪从远场边界反射。自由水面网格仅在初始时刻生成一次，并采用弹簧近似法在不改变网格节点顺序的情况下对瞬时水面进行网格重构，以避免数值不稳定。［结果］在验证所提出数值模型有效性和精确性的基础上，针对某浮式生产储卸油轮（ FPSO）模型的水动力特性进行数值模拟，发现考虑非线性影响时 FPSO的运动响应在共振区段明显增大，证明了传统线性方法的预报结果趋于危险。［结论］研究成果既可为海洋浮式结构物的设计提供更可靠的预报工具，也可为其实际应用提供理论依据。关键词：时域；完全非线性；高阶边界元；混合欧拉—拉格朗日方法；水动力特性中图分类号: U661.1 文献标志码:A DOI：10.19693/j.issn.1673-3185.01726 Fully nonlinear numerical simulation on hydrodynamic characteristics of FPSO SUN Lei*1,2, LU Tingting1, DENG Xiaoxiao1, LIU Changfeng3 1 School of Naval Architecture Engineering, Dalian University of Technology, Dalian 116024, China 2 Collaborative Innovation Center for Advanced Ship and Deep-Sea Exploration, Shanghai 200240, China 3 School of Marine and Civil Engineering, Dalian Ocean University, Dalian 116023, China Abstract: ［Objectives］In order to improve the safety performance of marine structures, fully nonlinear numerical simulation study is carried out on the problem of wave-structure interaction.［Methods ］ Based on the three-dimensional fully nonlinear time-domain potential flow theory and the Higher-Order Boundary Element Method (HOBEM), an open water model for the interaction between waves and structures is established. The total velocity potential is decomposed into incident potential and scattering potential by the velocity potential separation technique. The Mixed Euler-Lagrange (MEL) method is used to track the particle in the instantaneous free surface, and the fourth-order Runge-Kutta method is used to update the instantaneous free surface. In order to calculate wave loads, a virtual function is introduced instead of directly predicting the time derivative of the velocity potential. To prevent waves from reflecting from the far field boundary, an artificial damping layer is placed on the outside of the free surface. In addition, the free surface mesh is generated only once at the initial time, and then the spring approximation method is used to reconstruct the instantaneous free surface mesh without changing the order of mesh nodes, so as to avoid the possible numerical instability. ［Results］On the basis of validating and verifying the present numerical model, a Floating Production, Storage and Offloading (FPSO) unit is numerically simulated to study its hydrodynamics. It is found that the motion response of FPSO increases obviously in the resonance section when considering the nonlinear effect, and it is proved that the prediction of the motion of FPSO by the traditional linear method is tend to be risk. ［Conclusions］The work of this paper will provide a more reliable prediction tool for the design of the offshore floating structures and also provide the necessary theoretical basis for the practical application of the offshore floating structures. Key words: time domain； fully nonlinear； Higher-Order Boundary Element Method (HOBEM)；Mixed Euler-Lagrange (MEL) method；hydrodynamics

0 引 言

随着海洋开发逐渐向深海领域发展，海洋环境更加恶劣，在大波幅波浪作用下，海洋结构物载荷与运动的非线性特性显著增强。目前，常用的商用软件（例如 AQWA，Seasam，Hydrostar，Wamit等）都是基于线性或弱非线性的假定进行数值模拟，这必然会导致计算精度不足。因此，针对波浪和海洋结构物相互作用的问题，开展完全非线性数值模拟研究，对于提高海洋工程装备的安全性具有重要意义。1976 年， Longuet-Higgins 等[1] 基于势流理论提出了一种完全非线性方法，采用混合欧拉—拉格朗日方法追踪自由水面质点，自由水面和物面边界条件得到实时满足，每一时刻都对网格进行重构。近年来，随着计算机运算能力的不断提高，完全非线性方法得到了快速发展。周斌珍[2]采用完全非线性方法建立了开敞水域模型，并对漂浮圆柱的受迫振动、绕射问题以及Ringing现象进行了深入研究，揭示了非线性作用的机理。Bai 等[3-6] 基于高阶边界元法（HOBEM），运用区域分解技术，建立了波浪水池的完全非线性三维数值模型，并对圆柱的受迫振动以及非线性波浪作用下圆柱的运动特性及波浪爬高特性等进行了分析研究。宁德志[7]采用快速多极子去奇异边界元法（ Fast Multi-pole Boundary Element Method）分别建立了非线性三维开敞水域模型及非线性数值波浪水槽模型，并对圆球的强迫纵荡、垂荡和转动问题以及波浪对圆柱的绕射等问题进行了探讨。Feng 等[8-9] 利用完全非线性数值模型对并列双箱在波浪作用下的水动力共振问题进行了研究。然而，上述研究都仅局限于简单几何形状物体，鲜有考虑船舶这种带有复杂曲面形状的海洋结构物。本文将针对某浮式生产储卸油轮（ Floating Production Storage and Offloading，FPSO）模型进行数值模拟研究。首先，综合考虑自由水面、物面以及入射波浪的非线性影响，基于势流理论，结合高阶边界元法，建立波浪与海洋结构物相互作用的三维完全非线性开敞水域模型。然后，采用混合欧拉—拉格朗日（MEL）方法追踪瞬时自由水面，利用四阶 Runge-Kutta 法预报更新下一时间步的水面、物体位移和速度，每一时刻都对自由水面和物体湿表面网格进行更新重构。最后，通过对自由漂浮的圆柱、无航速的Wigley 船模型进行模拟，验证所提模型的准确性，并模拟某FPSO模型的水动力特性。1 数学模型1.1 初边值问题图 1 所示为建立的完全非线性开敞水域模型。图中：SF 为瞬时自由水面边界；SB为瞬时物体湿表面；SD 为水底边界；d 为水深；n 为物面单元单位法向矢量，以指出流体方向为正方向。定义 2组坐标系，分别描述流体域和物体的运动。一个是固定坐标系O-xyz，其中 z =0 位于静水面上，垂直向上为正；另一个是船体坐标系 O'-x'y'z'，其中原点O'在船体质心上。当物体静止处于平衡位置时，2组坐标系相互重合。假定流体为无黏性、不可压缩、无旋运动的理想流体，则整个流域可用速度势函数描述，总的流场速度势ϕ满足 Laplace 方程：

由物质导数和伯努利方程推导可得瞬时自由水面的完全非线性运动学和动力学边界条件。本文采用混合欧拉—拉格朗日方法追踪瞬时自由水面，其边界条件在拉格朗日系统下可以写为

式中： d/dt = ∂ / ∂ t + ∇ϕ ·∇ ，为物质导数； X(x, y, z)为瞬时自由水面流体任意质点的位置矢量；η为自由水面高程；g为重力加速度。在瞬时物体湿表面SB 上，满足如下边界条件：

式中，V为物面单元速度矢量。针对固定边界，满足V = 0；针对运动的物体，V可以表示为

式中：U为结构物的平动速度；Ω为结构物的转动速度；rb为位置矢量。

1.2 速度势分离本文采用 Ferrant 等[10-11] 提出的入射势和散射势分离技术，将总波浪( ϕ,η )分解为入射波浪(ϕI ,η I)和散射波浪(ϕS ,η S)，即

式中的入射波浪通过高阶Stokes 理论解给定。将式 (6) 代入式 (1)~式 (4)，整理可得关于散射波浪(ϕS ,η S)的初边值问题。

式中：r为流体质点与空间固定坐标系原点O的水平距离；x0 和 y0为初始时刻流体质点的水平坐标；μ(r)为阻尼系数函数，其作用是对流体出口处的散射波浪进行消波，以防止波浪反射对结构物造成二次作用，具体表示为

式中：ω 为波浪圆频率；λ 为入射波波长；r0 为阻尼层起点半径；α0和 β0 为关系数，本文中，α0=β0=1。为了避免初始效应导致数值不稳定，干扰模拟进程，本文在与 (ϕI, ηI) 有关的项上作用一个缓冲函数：

式中，Tm为缓冲时间，一般取入射波周期的整数倍，这里取为2倍。时域分析理论是一种完全模拟实际作用过程的理论，其控制方程中不仅含有空间项，还含有时间项，因此还需要给出初始条件：在整个计算域内应用格林第二定理，将上述关于散射波浪 (ϕS, ηS) 的边界值问题转化成如下的边界积分方程进行求解：

式中：p= (x0, y0, z0)，为源点坐标；q= (x, y, z)，为场点坐标；α(p) 为固角系数；S1 为整个计算域表面； G(p, q)为格林函数。由于整个计算域是关于O-xz 平面对称，且海底是水平的，所以选取简单Rankine 源及其关于对称面（ y =0 ）和海底（ z =−d ）的镜像为格林函数。鉴此，在计算时只需考虑一半计算域，而且海底被排除在外，式（17）则为G ( p , q )格林函数。这里，采用HOBEM 法，将计算域离散为8节点的四边形单元，并将其表示为参数坐标下的等参单元。在每个单元内，则采用二次形函数来描述任意点的几何坐标、速度势和速度势的法向导数，即

，( n)k式中：( ξ,ς)为参数坐标；Xk，ϕk ∂ ϕ/ ∂ ，hk 分别为单元上第k个节点的坐标、速度势、速度势法向导数和形函数；K为单元内节点的个数。将式(19)代入式 (16) 中，积分方程可以离散成如下形式：

式中： J ( ξ,ς)为雅各比行列式；Ne1 和 Ne2 分别为自由水面和物面边界上离散的曲面单元数。上述积分采用标准的 Gauss-Legendre 积分，将8 节点四边形单元离散为4×4个样本点进行求解。最终积分方程转换为如下线性方程组：

式中： A11 ，A12 ，A21 ，A22 均为系数矩阵A 的元素； B1 和 B2 均为向量B的元素。上述积分边界每一时刻都在变化，因此每一时间步的系数矩阵A和向量B都需重新建立再求解。计算过程中，假定当前时刻物面上的速度势导数和自由水面上的速度势均已知，由积分方程(20) 计算当前时刻的物面速度势和自由水面上的速度法向导数后，采用四阶 Runge-Kutta法，基于自由水面边界条件式(8)~式(11)，计算下一时刻水质点位置的和速度势。自由水面流体质点的速度分量可由式(22) 计算：

式中：nx，ny 和 nz 为法向向量n 沿 x，y 和 z 方向的分量； ∂ ϕS / ∂ ξ和∂ ϕS / ∂ ς可直接通过式 (19) 求解。需要注意的是，一个节点可能由几个单元共用，因此需对节点周围单元速度取平均值，得到该节点的速度分量。2 数值实现因为完全非线性方法的自由水面和物面边界条件都均实时满足，所以每个时间步都要重新捕捉及更新自由水面和物体湿表面网格。自由水面和物面的网格重构是完全非线性模型工作的重点及难点。2.1 自由水面网格更新本文采用弹簧近似法[12-13] 对自由水面动网格进行处理。弹簧近似法的本质是把网格节点看作由弹簧连接起来的，整个网格像弹簧系统一样变形。因此，只需要在初始时刻形成一次网格即可。网格节点位移可由式(23) 确定： 2.2 物面网格更新i

式中：δi 和 δj 分别为节点i和 j 的水平位移矢量； Ni 为网格节点总数；kij为节点i和 j之间的弹簧刚度。kij 的计算式为

式中， li j为节点i与 j之间的距离。更新后，网格节点 i的水平位置矢量Xnew可表示为

式中， Xorigin为节点 i在初始网格中的水平坐标。i至此，新的均匀网格节点位置已经完成更新。下面将通过二次形函数插值求解该网格节点上的速度势以及波面高程。为方便描述，将通过上述方法求得的网格记为网格1，将通过混合欧拉—拉格朗日方法（MEL）方法求得的网格记为网格2。为进行插值求解，需首先确定网格1中的节点在网格2中所属的单元，可以通过式(26)来确定：

式中： Sm 为网格2 中第m 个单元的面积；Sj 为三角形子元素的面积；该三角形由网格1中所求节点和第m 个单元的节点构成；M1 为三角形子元素的个数。当式 (26) 成立时，网格1所求节点一定在第m个单元的内部，此时可以通过式(27) 来求解节点的参数坐标( ξ,ς)：

式中：(xi, yi) 为网格1 上所求节点的坐标；(xkm, ykm)为网格2上第m个单元节点的坐标。通过上述处理得到的自由水面的计算域网格相对均匀，可避免因网格变形、扭曲而导致的数值不稳定性。对于船舶这种具有复杂曲面的海洋结构物来说，物面网格的更新较为困难。为了降低网格更新难度及工作量，本文将物面网格分为靠近水线区域和远离水线区域，如图2所示。因为远离水线区域始终处于瞬时水线以下，所以这部分网格可以不用更新。针对靠近水线区域的网格，通过背景网格插值技术予以更新。初始时刻通过读入船体各站的型线数据进行网格划分，包括水线以上的干舷部分和水线以下的部