考虑多场耦合的碰撞载荷下半潜式风机动力响应

2022-01-05王禹凯王加夏

中国海洋平台 2021年6期

王禹凯，宋娜，王加夏，刘昆

(江苏科技大学船舶与海洋工程学院，江苏镇江 212003)

0 引言

随着经济的快速发展，清洁的可再生能源逐渐成为各个国家开发的重点[1]，其中，风能以其无污染且储量大的特点成为近年来发展较快的新能源之一[2-3]。开发风能的方式包含陆上风机和海上风机。海上风机可分为固定式风机和海上浮式风机两种类型。海上浮式风机通常靠近航道，在恶劣天气下，船舶与其发生碰撞事故的概率会大幅提高。此外，风电场需要运维船定期维护，且风电场所在的海域常有较大的风浪，这也会增加碰撞风险。因此，研究半潜式风机在撞击载荷下的运动规律及结构响应情况对其耐撞性的优化及使用寿命的延长有重要的理论价值，对保证海上风电场的安全可靠性具有至关重要的意义。

半潜式风机的动态响应包含风、浪、流和系泊系统等复杂的多场载荷耦合问题。WAYMAN[4]采用线性频域水动力法开发一套在线性频域下计算水动力、气动力和结构耦合响应的分析程序。WANG等[5]提出一种基于AQWA水动力软件的浮式风机耦合响应的简单频域法，研究风、浪、流和系泊系统阻尼对浮式基础响应的影响，但频域法的局限性在于其不能捕捉非线性动态特性和瞬态效应。JONKMAN等[6]开发一套能全耦合考虑空气动力学、水动力学、控制系统、结构分析的程序以分析Barge型浮式风机的结构动力响应。MA等[7]基于此程序分析Hywind Spar浮式风机系统在不同海况下的运动和系统载荷响应特性，讨论风、浪对风机系统的影响，但忽略系泊系统的影响。为了更精确地模拟正常作业条件下浮式风机系统的性能，对浮式风机的风-浪-系泊进行全耦合数值模拟是必要的。

本文以DeepCwind浮式风机模型为研究对象，提出一种风-浪-系泊耦合场景下船舶-浮式风机碰撞动态响应数值分析方法。研究风-浪-系泊耦合场景下半潜式风机在船舶撞击载荷作用下的动力响应，有效地模拟半潜式风机在冲击载荷作用下的动力过程。

1 船舶-浮式风机碰撞耦合算法

1.1 三维势流理论

对于大尺寸构件(D/λ>0.2，其中：D为横截面直径；λ为波长)，入射波经过该结构时会受到严重的辐射和绕射，需要选取基于辐射和绕射理论的势流理论预报浮体结构的水动力系数和波浪力。当入射波与浮体结构互相作用时，波浪的速度势可分解为3个部分：入射波速度势ΦI，绕射波速度势ΦD和辐射波速度势ΦR。其中：绕射波速度势指的是由于浮体结构(海上风机平台)的存在，入射波流过浮体结构时产生绕射波；辐射波速度势是由海上风机平台自身的运动而产生的波浪，这是有关于置于初始静止流场产生微幅振荡的结构与这个流场之间的互相作用问题。

单色波入射情况下的一阶速度势为

Φ=Κ(ΦI+ΦD+ΦR)=Re{[ΦI,1(x,y,z)+ΦD,1(x,y,z)+ΦR,1(x,y,z)]eiω t}

(1)

式中：ΦI,1为一阶线性入射波速度势；ΦD,1为一阶线性绕射波速度势；ΦR,1为一阶线性辐射波速度势；ω为波浪频率；t为时间；Κ为波浪斜率且定义为

(2)

式中：A为波幅;L为波长；k为波数。

式(1)的一阶线性入射波速度势可由线性波浪理论计算得出：

(3)

式中：g为重力加速度，取9.81 m/s2；z为波浪高度；d为水深；k=kexcosθ+keysinθ，其中θ为波浪入射角，ex和ey分别为沿x轴和y轴的单位向量；x=xex+xey。

(1)绕射速度势需满足如下边界条件：

(4)

式中：n为物体表面的单位法向向量，方向朝外；r为物体表面的位置。

(2)辐射速度势需满足如下边界条件：

(5)

式中：ξ1为一阶平移运动；α1为一阶转向运动。

通常利用格林函数计算包含边值条件的拉普拉斯方程，然后建立速度势与格林函数之间的积分方程。需要利用边界单元法计算该积分方程，将物体离散为一定数目的单元，之后假设每个单元内的速度势都可利用单元节点势的多项式表达，基于伽辽金法或配点法建立节点势的线性方程组，得出各节点处的速度势。在得到绕射势和辐射势后，一阶线性水动压强和自由表面高程[8]为

(6)

(7)

式(6)和式(7)中：p1为一阶线性水动压强；ρ为海水密度，取1 025 kg/m3；Φ1为一阶线性速度势；η1为一阶自由表面高程，是时间和空间的函数。

沿物体的湿表面进行积分就能得到总波浪力(F1、F2、F3)和力矩(M1、M2、M3)为

(8)

(9)

式(8)和式(9)中：S为物体的湿表面。另外，一阶波浪力F1可定义为3个部分的和：

F1=FHS,1+FR,1+FEX,1

(10)

式中：FHS,1为一阶刚性静水回复力；FR,1为一阶辐射力；FEX,1为一阶波浪力(由绕射速度势和入射波速度势产生)。

1.2 多场耦合分析方法

图1为DeepCwind浮式风机的频域分析流程。在有限元软件ANSYS Workbench和Hypermesh中分别对DeepCwind浮式风机进行建模与网格划分。利用AQWA的Hydrodynamic Diffraction模块对浮式风机进行静水分析，得到附加质量惯性矩和静水回复刚度矩阵。由于水动力软件AQWA采用势流理论，无法计及黏性阻尼的影响，与实际情况会有较大差异，因此通过黏性阻尼经验公式计算临界阻尼和黏性阻尼并将其添加至软件中以达到阻尼修正的目的。关于半潜式风机基础黏性阻尼的计算，目前还缺乏成熟的方法。垂荡、横摇和纵摇方向上等效黏性阻尼可取临界阻尼的10%。临界阻尼计算式[9]为

图1 DeepCwind浮式风机频域分析流程

(11)

式中：M为质量；Kz为对应自由度的静水回复刚度。式(11)变为

(12)

式中：I为质量惯性矩；ΔI为附加质量惯性矩。

之后利用频域分析研究浮式结构在水中的运动响应。应用AQWA对DeepCwind浮式风机进行水动力性能数值计算，采用经典势流理论的计算方法，利用频域计算的快速性和线性化特点求解单位规则波下的半潜式风机平台的水动力参数，主要包括：附加质量系数、辐射阻尼系数以及不同频率、不同浪向下的一阶波浪力等。

在频域分析中没有考虑系泊系统的回复力和风载荷的影响，因此计算结果不能完全准确地反映半潜式风机的运动响应，而时域分析可用于船舶与海洋结构物在特定的海洋环境条件下的一定运动范围内的运动情况仿真，能够较真实地模拟船舶与海洋结构物的运动响应，并且可以考虑在各时间步内半潜式风机系统各部分之间的耦合作用，充分考虑前一时刻的运动响应对后一时刻的影响[10]。对于一些特定的问题，如结构物之间的撞击、海底搁浅都能清楚地反映其具体发生的时刻和全过程的时间历程。图2为DeepCwind浮式风机时域分析流程。将AQWA计算出来的水动力参数导入美国可再生能源实验室(National Renewable Energy Laboratory，NREL)开发的FAST程序[11]中，使用FAST程序考虑风、浪以及系泊系统与半潜式风机之间的耦合作用，计算半潜式风机的运动响应、气动推力、水动力和锚链力。

图2 DeepCwind浮式风机时域分析流程

在LS-DYNA中无法完成与结构动力响应相关的加载模式，即载荷无法随着结构响应变化而变化，无法达到载荷与结构耦合的效果，因此通过子程序Loadud LS-DYNA加载。

图3为载荷与结构响应相互耦合的实现流程，先将FAST程序中计算出的气动推力载荷和水动力载荷导入LS-DYNA进行结构动力学计算。由于在FAST程序中模型的动力学参考点为水线面中心点，水动力载荷也作用在水线面中心点，因此在LS-DYNA中，将水动力载荷通过子程序Subroutine Loadud作用在水线面中心点，即原点处。气动推力则进行适当简化，时域载荷作为一个集中力直接施加在风机结构上，具体作用点如图4所示。将锚链力进行适当简化后作用于半潜式风机平台系泊点处。当T

图3 耦合算法说明

根据半潜式风机的实际情况，将系泊系统作为边界条件。在LS-DYNA中，如果准确地建立锚链模型并考虑其与水的相互作用，将会额外带来较大的计算负担。因此，为提高计算效率并保证一定的计算精度，对半潜式风机的系泊系统进行适当简化。当半潜式风机小角度运动且位移较小时，可将系泊缆绳简化成线性弹性系统。在FAST中，利用中心差分方法对偏导数进行计算，得到系泊缆绳的弹性刚度矩阵为

(13)

在LS-DYNA中，在半潜式风机的3个系泊点处均采用3个转动和3个平移弹簧单元对系泊系统进行建模，如图4和图5所示。

图4 载荷添加位置示例

图5 弹簧单元有限元模型

2 数值计算模型

2.1 半潜式风机模型描述

以NREL提出的DeepCwind浮式风机模型为研究对象，该风机发电功率为5 MW，额定风速为11.4 m/s，基本参数[12]如表1所示。

表1 DeepCwind浮式风机模型主要参数

采用有限元软件ANSYS Workbench对Deep-Cwind浮式风机进行建模，风机基础主要由3个浮筒组成，并且有横支撑杆和斜支撑杆。此外，对风机模型进行简化，机舱及叶片采用集中质量代替。半潜式风机会产生6个自由度的平移和旋转运动(见图6)，平移运动包括横荡、纵荡和垂荡，旋转运动包括纵摇、横摇和艏摇。纵坐标原点是在水线处由3根圆柱组成的等边三角形的重心，z轴垂直向上，xy平面与海平面重合。

图6 半潜式风机运动示例

2.2 撞击船模型描述

考虑到深海风电场海域内实际的过往船舶状况，选用风电运维船作为撞击船。撞击船主要参数如表2所示。在计算中主要考虑DeepCwind浮式风机在碰撞中的结构动力响应，为缩短计算时长，撞击船采用刚性体且建模时对船体内部结构进行简化，并且省略其上层建筑。船舶有限元模型及网格如图7所示。

表2 撞击船主尺度

图7 碰撞船有限元模型及网格划分

将碰撞船有限元模型与DeepCwind浮式风机装配在一起，为模拟浮式风机在运动状态下遭受船舶撞击的效果，撞击船与浮式风机之间在发生碰撞前有一定的距离，撞击船与风机的撞击点位于风机吃水线附近。

此外，选取主-从接触面方法分别定义在半潜式风机浮筒和撞击船上，选用LS-DYNA有限元软件提供的自动面面接触(*CONTACT_AUTOMATIC_SURFACE_TO_SUFACE)，动、静摩擦因数均取0.2[13]。

在撞击船作用下，被撞半潜式风机平台的立柱和支撑构件可能会发生塑性变形，对于被撞半潜式风机采用常用的塑性动态模型(*MAT_PLASTIC_KINEMATIC)，预先假定一个最大塑性失效应变，若所计算的单元达到最大塑性失效应变，则该单元即刻失效，所对应的构件就认为已损坏而不再承受任何外力。为节约计算时间，撞击船舶采用刚性体(*MAT_RIGID)，半潜式风机和撞击船的材料分别为船用高强钢(Q345)和低碳钢(Q235)，其材料参数[14]如表3所示。

表3 钢材性能参数

撞击船在碰撞前的运动中，需要考虑船舶结构与水之间的相互作用，一般通过附加质量模型和流固耦合模型解决。考虑到流固耦合模型会使计算时间过长，且主要研究对象为半潜式风机，而附加质量模型方法已得到过许多学者[15]的使用与验证，因此选用附加质量模型处理船体结构与水之间的相互作用。此外，附加质量为撞击船质量的0.05倍。

2.3 碰撞场景

选取典型碰撞场景，在天气恶劣的条件下，海上运维船失控后与半潜式风机发生碰撞。根据DNV GL-ST-0126[16]和DNV-ST-0119[17]规范，当海上风机与船舶发生碰撞时，船舶速度不低于0.5 m/s，因此碰撞船舶的撞击速度均取1.0 m/s，撞击船与半潜式风机的撞击点位于风机吃水线附近。海洋环境则取半潜式风机额定风速11.4 m/s，波浪高度为1.5 m，风向和浪向均取y轴负方向，如图8所示。时间步长取0.05 s，模拟总时长为150 s。

图8 简化的碰撞场景示例

3 多场耦合下碰撞仿真结果及分析

3.1 碰撞力

图9为撞击船首部撞击DeepCwind浮式风机碰撞力时程曲线，该曲线总体呈明显的非线性状态，碰撞是在非常短的时间内发生的动力冲击问题，整个碰撞过程持续时间约0.4 s，在船首与DeepCwind浮式风机的浮筒接触开始，从0.05～0.20 s碰撞力急剧增加，这是由于随着风机浮筒参与的撞击范围不断增加，撞击船首部所受抵抗力不断增加，并于0.25 s时达峰值3.723×106N。此后，在0.25～0.45 s，碰撞力出现卸载情况迅速减小，这是由于在撞击船撞击风机浮筒的过程中，船速不断降低，动能下降，直至碰撞过程结束，碰撞力也减小为零。

图9 碰撞力时程曲线

3.2 碰撞能量

系统总体能量转化如图10所示。随着撞击船与DeepCwind浮式风机的浮筒接触，撞击船的动能主要转化为浮筒碰撞部位的变形内能和浮式风机整体动能两部分，还会产生少量的沙漏能、阻尼能和摩擦能。阻尼能和摩擦能的数量级极低，在图10中略去。由图10可知，碰撞始于25.95 s，由于浮式风机靠系泊系统进行定位，浮式基础有较大的灵活性，因此撞击船的动能大部分转换为浮式风机的动能，动能最大值达1.37×108J，少量转化为风机的内能和沙漏能，这与船-船、船-固定式风机等碰撞问题不同。这说明在船-浮式风机的碰撞运动中，撞击对浮式风机的损伤变形影响较小，但对浮式风机的运动响应影响较大。

图10 系统能量转化曲线

3.3 结构应力应变

图11为25.95 s时船舶与DeepCwind浮式风机发生正面碰撞的整体结构应力云图。由图11可知，在船舶碰撞浮式风机的过程中，所产生的应力主要分布于浮筒的撞击区域，存在明显的局部性。

图11 DeepCwind浮式风机与船首碰撞整体结构应力云图

图12为DeepCwind浮式风机遭受撞击的浮筒在不同时刻的结构等效应力云图。由图12可知，在25.95～26.05 s，风机浮筒的结构应力峰值逐渐增大，并于26.05 s达最大值3.626×108Pa，随后随着碰撞力的卸载而减小，且结构应力分布有明显的局部性。

图12 撞击过程中不同时刻浮筒结构应力云图

图13为DeepCwind浮式风机遭受撞击的浮筒在不同时刻的结构塑性应变云图。由图13可知，结构塑性应变在25.95～26.05 s呈逐渐增大的趋势，结构塑性应变区域集中于撞击部位且损伤面积较小，呈现明显的局部性，说明撞击载荷对DeepCwind浮式风机的结构损伤较小，撞击船的动能只有一小部分转化为内能。

图13 撞击过程中不同时刻浮筒结构塑性应变云图

3.4 运动响应

图14为 DeepCwind浮式风机在有无撞击载荷下的横荡、横摇、艏摇和垂荡运动响应对比图。风、浪和船舶撞击方向皆为y轴负方向，纵摇和纵荡因量级较低而忽略。由图14可知，各曲线都呈现明显的非线性状态，垂荡自由度的运动响应与其他不同，在有撞击载荷下垂荡自由度运动幅度低于无撞击载荷下的运动幅度，但是有撞击载荷下的横荡、横摇和艏摇运动幅度明显大于无撞击载荷下的运动幅度，说明撞击载荷对DeepCwind浮式风机的运动响应影响较大，这是由于撞击船所产的动能大部分转化为风机的动能，这些动能使DeepCwind浮式风机的横荡、横摇和艏摇运动幅度加大。