漂浮式风－浪能混合利用系统运动模型数值分析

韩端锋，丁松，马庆位

(1.哈尔滨工程大学船舶工程学院，黑龙江哈尔滨150001;2.伦敦城市大学工程及数学科学学院，伦敦英国EC1V 0HB)



漂浮式风－浪能混合利用系统运动模型数值分析

韩端锋1，丁松1，马庆位2

(1.哈尔滨工程大学船舶工程学院，黑龙江哈尔滨150001;2.伦敦城市大学工程及数学科学学院，伦敦英国EC1V 0HB)

摘要:风浪能混合发电技术具备可提高平台利用效率、增加并稳定发电量和减小发电成本的特点。本文以STC型漂浮式风浪能混合利用系统为例，建立并求解风浪能混合利用系统的非线性联合运动方程组，分析风浪能混合利用系统的波浪载荷、风载荷。计算风浪能混合利用系统的运动响应，并将数值计算的结果与已有的试验数据进行比较，验证了风浪能混合利用系统数学运动模型的准确性和适用性。应用该数学运动模型分析PTO(power takeoff)参数BPTO对波能装置发电性能的影响。在一定范围内，增加BPTO使俘获宽度比曲线带宽变小，且俘获宽度峰值对应入射波频率向低频区转移。

关键词:浮式风机;波能浪装置;风浪能混合发电;数值模拟;联合运动方程;能量吸收装置;运动模型

研究漂浮式风电技术对于开发我国深水海域风电有重要的理论价值和长远的战略意义。世界上很多地方已开始对漂浮式风电技术进行研究［1－2］。但是，深海风电技术还面临着若干问题。一是像陆上风电一样，海上风电波动性、间歇性、不规则性强、稳定性差。二是海上风电开发成本较高，企业的回报率较低。

为提高海上风电系统的稳定性、降低发电成本，可利用漂浮式风电系统中的平台结构，将波浪能装置与风机结合，形成风－浪能混合发电系统。目前已有一些关于这种风浪能混合利用系统的设计方案。Peiffer将WindFloat浮式风机分别与振荡浮子式波能装置、振荡水柱式波能装置和摆式波能装置结合［3－5］。Muliawan等开发了STC型风浪能混合利用系统(spar－torus－combination)，将Spar类型的浮式风机与振荡浮子式波能装置结合［6］。由于波能装置和浮式平台的多样性，关于风浪能混合利用系统尚无明确的公认分类方式，Pérez－Collazo对风浪能混合利用系统的发展现状进行了综述［7］。

本文以STC型风浪能混合利用系统为例，建立该风浪能混合利用系统的数学运动模型。该运动模型考虑了Spar平台六自由度运动以及浮子沿Spar平台轴向的升沉运动。计算其作业过程中的风载荷、波浪载荷。编写MATLAB程序对该风浪能混合利用系统的运动进行数值模拟。在两种不同操作模式下，将数值计算结果与已有的试验数据进行对比，验证了数学计算模型的准确性。最后，应用本文的计算方法，初步研究PTO参数的选取对发电性能的影响。

1　联合运动方程

在船舶与海洋工程计算中，通常将线性化的运动方程应用于小幅度的运动［8］。然而，当物体面临陡波，运动响应较大，此时线性化的运动方程并不能满足计算精度的需求。本文在计算风浪能混合利用系统的运动时，将其视为刚体，不考虑受力后的变形。因而，系统的运动可用刚体的完全非线性运动方程来描述。在推导中，为了书写简洁，常以WEC代表垂荡浮子式波能装置。

1.1建立坐标系

在描述风浪能混合利用系统运动的过程中，共建立3个坐标系:

绝对坐标系OXYZ，其中O为绝对坐标系的原点，可取空间中的任意一固定点，OZ轴垂直向上，OX轴指向船艏(纸面的右侧)，OX、OY、OZ符合右手定则;随体平动坐标系O'X'Y'Z'，它是一个辅助坐标系，只随SPAR坐平移运动而不作旋转运动;固体坐标系oxyz，整个坐标系被固定在Spar平台上，其坐标原点o取在Spar平台的重心，oz轴沿Spar平台轴向中心线垂直向上。

在初始时刻，3个坐标系重合，三者关系如图1所示。

图1　描述风浪能混合利用系统运动的坐标系Fig.1 Three coordinate systems used in motion simulation

在此风浪能混合利用系统的运动中，浮子仅能沿Spar平台轴向运动，也就是说浮子与Spar平台具有相同的旋转运动。绝对坐标系OXYZ与固体坐标系oxyz之间的坐标转换关系为

式中:T为由欧拉角(α，β，γ)组成的随体平动坐标系O'X'Y'Z'与固体坐标系oxyz之间的坐标转换关系矩阵。

在本文所建立坐标系中，对于任意一个向量L的时间导数计算式为:

设浮子在固体坐标系oxyz下的位移xb、速度ub、加速度，固体坐标系oxyz原点o在绝对坐标系OXYZ中的位移X0、速度U0、加速度，固体坐标系oxyz旋转角速度Ω、角加速度。

根据式可得到浮子在绝对坐标系OXYZ下的位移XW、速度UW、加速度:

式中:Ω×xb为牵连速度，固体坐标系oxyz绕其原点o以角速度Ω转动并带动浮子一起转动而引起的速度;为牵连加速度，其中是由于角速度Ω大小发生改变所产生的，第二项Ω×(Ω×xb)是由于固体坐标系oxyz以角速度Ω转动所产生的;2Ω×ub为科氏加速度，是由牵连运动与相对运动相互影响所产生的。

1.2联合运动方程

Spar平台和浮子在计算中被假设为刚体，Spar平台的运动姿态可根据其质心的位移X0和速度U0、随体平动坐标系O'X'Y'Z'相对于固体坐标系oxyz之间的欧拉角(α，β，γ)以及角速度Ω来描述;由于浮子仅能沿Spar平台轴线上下运动，两者具有相同的旋转运动，因此，浮子的运动姿态仅需相对坐标系oxyz下浮子的位移xb、速度ub便可以确定。基于1.1节的推导，给出该风浪能混合利用系统的运动联合方程:

式中:FS、NS为Spar平台的合外力、合外力矩;FW为浮子的合外力;MS、IS为Spar平台的质量矩阵、转动惯量矩阵;MW为浮子的质量矩阵;B是由欧拉角(α，β，γ)组成的矩阵:

为便于计算，Spar平台的平动运动方程基于绝对坐标系OXYZ，而浮子的平动运动方程、Spar平台的转动运动方程基于固体坐标系oxyz。因此MS和IS中所有非对角线元素均为零。另外，如系统无初始角速度，基于细长体理论的假设条件，可以得Ω× ISΩ=0，此时Spar平台的运动仅有5个自由度。

2　风浪能混合利用系统载荷分析

STC型风浪能混合利用系统在风、波浪及锚泊等外部载荷的作用下运动，本文按照各载荷的来源对其进行分类，分别计算各种载荷，再将各载荷矢量相加可得到作用在Spar平台、浮子上的合外力与合外力矩，最后代入联合运动方程的右侧并进行数值求解。因此，将方程中合外力FS、FW分解为

式中:GS、GW为Spar平台、浮子的重力;FWAVE－S、FWAVE－W为Spar平台、浮子的波浪载荷;FMOOR为锚泊力;FWIND为风载荷;FPTO－S、FPTO－W为Spar平台、浮子之间由于波能吸收装置产生的相互作用力。

2.1波浪载荷分析

本文基于细长体理论，将水动力分为非粘性力和粘性力两部分，非粘性力应用Rainey提出的完全细长体公式进行计算［9－10］，将非粘性水动力分成3部分:单位浸没长度力F1、底端力F2和自由表面力F3;粘性力Fdrag则应用莫里森公式中的阻力项进行计算。由于计算原理一致，在本节对波浪载荷分析中，以Spar平台为例进行说明。

2.1.1单位浸没长度力F1

单位浸没长度力表达式如式(11)，将单位长度浸没力沿Spar平台水下垂直方向积分得到该项波浪力作用于Spar平台的合力:

2.1.2底端力F2

底端力F2作用点位于Spar平台底部，表达式为式(13)，其中压力项p=ph+ps，为水动压力与静水压力之和，静水压力项psSl与式中ρS{ g }n一起，提供Spar轴向的浮力和静水回复力。式中后两项与Munk力矩相关，与波幅的平方成比例。

2.1.3自由表面力F3

自由表面力F3作用于Spar平台与波浪瞬时湿表面的相交处，具体计算式如下

式中:α为Spar平台轴线与波浪自由表面法向量之间的夹角，它的取值由Spar平台的运动参数和波浪瞬时自由表面的空间位置共同决定;t为垂直于Spar平台轴线的单位切向量，方向指向流体外侧，从式中不难看出，自由表面力与波幅的三次方成比例。

2.1.4粘性水动力Fdrag

单位长度Spar平台的粘性力采用莫里森公式中的拖曳力项计算:

式中:CD为拖曳力系数，wn为相对速度w在与SPAR轴线垂直的平面内分量。

2.1.5附加质量系数和拖曳系数的取值

应用Sesam软件，对Spar平台、浮子的水下部分建模，如图2所示。在得到附加质量曲线后，根据入射波波浪周期，即可得到该工况下的各附加质量系数。例如，在之后的对比计算中，入射波波浪周期为12 s，计算得到Spar平台水平方向附加质量系数为0.8，竖直方向附加质量系数为0.06;浮子水平方向附加质量系数为0.8，竖直方向附加质量系数为1.5。拖曳力系数的选取参考文献［8］中的取值，Spar平台的拖曳系数为0.6，波能装置WEC的拖曳系数为1。

图2　Spar平台、浮子的网格模型Fig.2 Mesh model for spar platform and WEC

2.2风载荷FWIND的计算

关于风浪能混合利用系统中风载荷，本文根据Statoil针对浮式风机风载荷的处理方法［11］，将风载荷简化为水平推力，作用于风机转子中心处。根据推力系数曲线，将该推力表示为与相对风速urel相关的函数进行实时计算，具体计算式为:

式中:ρa为空气密度;A为转子扫过的总面积;urel为风机转子中心处的相对风速;CT(urel)与相对风速相关的推力系数函数。

2.3锚泊力FMOOR的计算

对于锚泊力的计算，将锚泊线简化为不同刚度系数K的非线性弹簧，在计算锚泊力时考虑了导缆孔的位置锚泊力计算式如下［13］

式中:K为刚度系数;ΔXf为锚泊线的长度变化量。

2.4能量吸收阻尼力FPTO的计算

波能装置的工作原理是将其沿Spar平台轴线垂荡运动的动能通过能量吸收装置(power take－off)最终转化为可传输的电能，FPTO－W被视为能量吸收装置在浮子工作期间对其运动产生的‘阻力’，计算式为

式中:BPTO、KPTO分别为能量吸收装置的阻尼系数和回复刚度系数，xb、ub分别为浮子距平衡位置位移、运动速度。显然FPTO－W与FPTO－S是一对相互作用力，即两个等值反向的矢量。

2.5Spar平台合力矩NS的计算

至此已将式中FS、FW的各项分力表达式给出，下面给出Spar平台转动时受到的合力矩的表达式:

根据2.1～2.4中对单位浸没长度力F1、底端力F2、自由表面力F3、粘性力Fdrag、风载荷FWIND、锚泊力FMOOR的描述以及各分力的作用点到Spar平台重心处(即固体坐标系原点)的位置矢量，可得到上述各分力对Spar平台转动运动形成的力矩:

式中:rb为Spar平台底部相对于固体坐标系原点的位置矢量;ru为自由表面与Spar平台轴线相交处相对于固体坐标系原点的位置矢量;rWIND为风机转子中心相对于固体坐标系原点的位置矢量;rMOOR为锚链导缆孔相对于固体坐标系原点的位置矢量。

3　联合运动方程的数值解法及验证

经过之前联合运动方程的推导以及方程组右侧合外力、合外力矩的分析，对整个风浪能混合利用系统的联合运动方程已经有了清晰明确的表达。然而，整个风浪能混合利用系统的运动方程是非线性方程，方程中各运动参数是相互耦合的，如欧拉角矩阵B，很难单独求解其中的某个运动参数。此外，联合方程右侧的各项载荷是随Spar平台、浮子的位移、速度等运动参数而实时变化的。因此，在通常情况下，常利用数值计算的方法求解风浪能混合利用系统的运动响应。

3.1数值计算步骤

本文的数值计算程序在MATLAB中编写，数值计算步骤如下:

1)给出系统初始时刻的运动参数即Spar平台位移、速度、欧拉角、角速度;浮子在固体坐标系oxyz中位移、速度。本文计算中假定它们在初始时刻为静止状态。

2)由Spar平台、浮子的运动姿态，计算它们的瞬时浸没长度、相对流速、相对风速以及导缆孔在绝对坐标系中的位置。

3)根据Spar平台、浮子的瞬时浸没长度、相对流速、相对风速以及导缆孔在绝对坐标系中的位置计算Spar平台、浮子各自的合外力、合外力矩。

4)根据之前的推导，将计算好的各合外力、合外力矩代入联合运动方程式(7)中，并整理为一阶非线性常系数微分方程组，采用自适应阶的四阶龙格库塔法求解Spar平台、浮子的联合运动方程，在每一个时间阶内，系统的各项运动参数可被求解，通过迭代法达到所要求的误差控制限，输出并记录Spar平台、浮子的运动参数。

5)返回步骤2)，进行下一时刻的计算。计算结束于最终计算时刻。

3.2数学运动模型验证

对于本文算例中的风浪能混合利用系统STC，共有3种工作模式:Released模式、MWL模式和SUB模式。Released模式中浮子可沿Spar平台轴线方向自由运动，如处于极限工况下，浮子停止吸收、转化波浪能;MWL模式中浮子与Spar平台锁死，无相对运动，浮子仍处于自由水面位置;SUB模式中浮子与Spar平台锁死，无相对运动。系统整体下沉，浮子处于水面之下。在本文中，分别选取Released模式和MWL模式作为对比算例，进行数值计算，并与已有的试验数据进行对比验证［14］。需要指出的是，两种模式的对比算例中不考虑风载荷，且波浪输入为微幅波。

3.2.1Released模式

Released模式的对比试验中波浪输入为规则波，波浪周期1.697 s，波高0.08 m，试验的Froude比例为1∶50。图3中分别为Spar平台横荡、垂荡、纵摇运动的位移时历对比曲线以及浮子相对升沉运动的位移时历对比曲线。

图3　Release模式下数值计算结果与试验数据对比Fig.3 The comparison between numerical calculation and experimental result in Release Mode

3.2.2MWL模式

MWL模式的对比试验中波浪输入为规则波，波浪周期1.273 s，波高0.04 m，试验的Froude比例为1∶50。由前所述，MWL模式下浮子与Spar平台锁死，无相对运动，在该模式的数值计算中，将Spar平台与浮子简化成一个刚体，图4为系统纵摇运动对比图。从图可以看出，Release模式和MWL模式下数值计算结果均与试验数据相一致，初步验证了本文关于风浪能混合利用系统的模型建立、数值计算方法的准确性和适用性。

图4　MWL模式下纵摇运动数值计算结果与试验数据对比Fig.4　The pitch motion comparison between numerical calculation and experimental result in MWL Mode

4　PTO阻尼参数分析

由FPTO－W和ub可以得到一段时间内的浮子平均输出功率，也就是反应了瞬时输出功率曲线的均值，根据式(19)，平均输出功率的计算公式为

式中:第2项的积分几乎为零，浮子输出功率的来源是第一项，因此，阻尼系数BPTO的取值对波能输出功率具有很大的影响。

应用本文的数学运动模型，计算不同PTO阻尼系数BPTO下吸波浮子的平均输出功率，再根据单位宽度内的波能，得到浮子俘获宽度比。如图5所示，入射波为波高2 m的一组规则波。

图5　描述风浪能混合利用系统运动的坐标系Fig.5 WEC absorbed power width under different PTO damps

经过数值计算得到的各工况下浮子吸波效率曲线如图5所示，通过分析这些图表，可以得出:当BPTO增大时，波能吸收效率曲线的峰值增大，例如，波高2 m时，吸波效率曲线的峰值由BPTO为1 000 kN·s/m时的21.6％增长到PTO阻尼系数为8 000 kN·s/m时的57.6％;而随着BPTO的增大，峰值所对应的波浪频率减小，带宽也随之减小，例如，当BPTO为1 000 kN·s/m时，吸波效率曲线峰值对应的入射波频率为1.3 rad/s，设备对0.9～1.7 rad/s波浪频率区间的入射波都有良好的吸收效率(吸波效率大于15％)，而当BPTO增长到8 000 kN·s/m时，吸波效率曲线峰值对应的入射波频率减小到0.7 rad/s，此时设备仅对0.7～0.9 rad/s波浪频率区间的入射波有良好的吸收效率。

5　结论

本文就风浪能混合利用系统STC的运动响应进行了数值分析，得出以下结论:

1)建立了适用于STC风浪能混合利用系统的数学运动模型:基于刚体假设条件，建立包含Spar平台、浮子在内的风浪能混合利用系统的非线性联合运动方程。编写MATLAB程序，运用自适应阶的四阶龙格库塔法求解联合运动方程组，模拟风浪能混合利用系统的运动响应。在Release模式和MWL模式下，分别将系统运动的数值计算结果与已有的试验数据进行验证，数值计算结果与试验数据相符合，验证了本文关于风浪能混合利用系统的模型建立、数值计算方法的准确性和适用性。

2)应用该数值方法探讨了PTO阻尼系数B_ PTO对吸波浮子发电性能的影响:在一定范围内，增加B_PTO将使得波能吸收效率曲线峰值"变尖"，且峰值对应的入射波频率向低频区转移。通过分析阻尼系数B_PTO的取值对波能输出功率的影响，为之后波能吸收装置(PTO)系统的设计提供依据。

参考文献:

［1］张亮，吴海涛，荆丰梅，等.海上漂浮式风力机研究进展及发展趋势［J］.海洋技术，2010，29(4):122－126.ZHANG Liang，WU Haitao，JING Fengmei，et al.Research progress and development trend of offshore floating wind turbine［J］.Ocean technology，2010，29(4):122－126.

［2］SHIM S，KIM M H.Rotor－floater－tether coupled dynamic analysis of offshore floating wind turbines［C］//The Eighteenth International Offshore and Polar Engineering Conference.International Society of Offshore and Polar Engineers，2008.

［3］PEIFFER A，RODDIER D，AUBAULT A.Design of a point absorber inside the WindFloat structure［C］//ASME 2011 30th International Conference on Ocean，Offshore and Arctic Engineering.American Society of Mechanical Engineers，2011:247－255.

［4］AUBAULT A，ALVES M，SARMONTO A，et al.Modeling of an oscillating water column on the floating foundation WindFloat［C］//ASME 2011 30th International Conference on Ocean，Offshore and Arctic Engineering.American Society of Mechanical Engineers，2011:235－246.

［5］PEIFFER A，RODDIER D.Design of an oscillating wave surge converter on the windfloat structure［C］//4th International Conference on Ocean Energy.Dublin:ICOE2012.2012.

［6］MULIAWAN M J，KARIMIRAD M，MOAN T.Dynamic response and power performance of a combined spar－type floating wind turbine and coaxial floating wave energy converter［J］.Renewable Energy，2013，50:47－57.

［7］PEREZ－COLLAZO C，GREAVES D，IGLESIAS G.A review of combined wave and offshore wind energy［J］.Renewable and Sustainable Energy Reviews，2015，42:141－153.

［8］MEI C C.The applied dynamics of ocean surface waves ［M］.New York:World Scientific，1989:285－290.

［9］RAINEY R C T.Slender－body expressions for the wave load on offshore structures［C］//Proceedings of the Royal Society of London A:Mathematical，Physical and Engineering Sciences.The Royal Society，1995，450(1939):391－416.

［10］MA Q W，PATEL M H.On the non－linear forces acting on a floating spar platform in ocean waves［J］.Applied O－cean Research，2001，23(1):29－40.

［11］NIELSEN F G，HANSON T D，SKAARE B.Integrated dynamic analysis of floating offshore wind turbines［C］//25th International Conference on Offshore Mechanics and Arctic Engineering.American Society of Mechanical Engineers，2006:671－679.

［12］HU Zhou，WAN Decheng.Numerical investigations on the aerodynamic performance of wind turbine:downwind versus upwind configuration［J］.Journal of Marine Science and Application，2015，14(1):61－68.

［13］XU Zhengqiang，HUANG Shan.Numerical investigate of mooring line damping and the drag coefficients of studless chain links［J］.Journal of Marine Science and Application，2014，13(1):76－84.

［14］WAN L，GAO Z，MOAN T.Model test of the STC concept in survival modes［C］//ASME 2014 33rd International Conference on Ocean，Offshore and Arctic Engineering.A－merican Society of Mechanical Engineers，2014:V09AT09A010－V09AT09A010.

Numerical modeling of floating wind-wave energy hybrid utilization system

HAN Duanfeng1，DING Song1，MA Qingwei2
(1.College of Shipbuilding Engineering，Harbin Engineering University，Harbin 150001，China;2.School of Engineering and Mathematical Sciences，City University London，UK EC1V 0HB)

Abstract:Combining wind turbine system and wave energy converter(WEC)is a new trend which will be cost effective and lead to an overall improvement on power performance.In this paper，the Spar－Torus Combination(STC)type floating wind－wave energy hybrid utilization system is studied.A MATLAB program was built to analyze the wave load and wind load of the hybrid system.The motion response of the hybrid system was calculated and the result was compared with the experimental data that had been gotten.The accuracy and adaptability of the numerical model for the hybrid system were validated.Then this numerical method was used to study the effect of WEC's power take－off system on its power performance.As the PTO coefficient increased，the peak frequency of capture width ratio curve decreased，and the corresponding incident wave frequency turned to low－frequency area.

Keywords:floating wind turbine;wave energy converter;hybrid wind－wave energy system;numerical simulation;motion equation;power take－off;motion model

通信作者:丁松，E－mail:dingsong2014@ hotmail.com.

作者简介:韩端锋(1966－)，男，教授，博士生导师;丁松(1989－)，男，博士研究生.

基金项目:国家自然科学基金资助项目(51379051).

收稿日期:2015－10－27.网络出版时间:2016－01－04.

中图分类号:O35

文献标志码:A

文章编号:1006－7043(2016)01－0138－07

doi:10.11990/jheu.201510064

网络出版地址:http://www.cnki.net/kcms/detail/23.1390.u.20160104.1648.014.html