海上浮式风机时域耦合程序原理及其验证

陈嘉豪, 刘格梁, 胡志强

(1. 上海交通大学 海洋工程国家重点实验室; 高新船舶与深海开发装备协同创新中心，上海 200240; 2. 纽卡斯尔大学 工程学院, 纽卡斯尔 NE1 7RU，英国;3. 中国能源建设集团广东省电力设计研究院有限公司，广州 510663)

近十年来，由于各国日渐重视可再生能源的转型发展，海上风电技术及其相关产业在全球呈现蓬勃发展的态势[1].而随着近海优质风电资源的日益开发，海上风电产业也如同过去油气产业所经历的变革一样，逐渐从陆地走向海洋，从浅海过渡到中等或更深的水域，其相应的支撑平台形式也从固定式逐渐过渡为浮式[2].但由于海上浮式风机具有复杂的多物理场耦合特性，如何准确、快速地对海上浮式风机的耦合动力特性进行数值预报，成为一直以来的研究难点与热点问题之一[3].

目前，海上浮式风机的数值分析主要分为频域计算和时域计算.早期研究人员主要借鉴离岸油气工业的经验，利用频域分析法对海上浮式风机进行设计与分析[4-5].虽然，频域分析法能够快速简便地计算海上结构物的频域动力响应特性，但由于海上浮式风机涉及空气动力学、水动力学、结构动力学、锚泊动力学及控制系统等技术，其动力响应具有高度的非线性特性.频域分析法无法全面、准确地分析海上浮式风机的耦合动力响应特性.在这之后，国际上一些研究人员改写或调用固定式风机或机械行业的计算软件，用以时域计算并分析海上浮式风机的耦合动力特性，如FAST[6]、 HAWC2[7]、 Bladed[8]、 Adams、 Simpack等软件.在国内，唐友刚等[9]利用水动力计算软件(HydroD)对自主设计的新型海上浮式风机在频域范围内进行动力学仿真.Ma等[10]利用FAST软件对立柱型浮式平台耦合动力响应进行时域计算分析.赵文超等[11]利用计算流体动力学(CFD)滑移网格技术对海上浮式风机的气动性能做了仿真计算.刘利琴等[12]针对垂直轴的海上浮式风机开发了气动-水动力耦合数值仿真程序.当前，国内外关于海上浮式风机的全耦合时域仿真程序尚在研究中，且国内对于海上浮式风机相关理论的研究及数值程序的开发与国外同行尚有一定差距.

本文介绍自主研发的用于海上浮式风机的气动-水动-结构-锚链-控制全耦合时域数值仿真程序(DARwind)的相关理论方法.在气动计算方面，利用非定常叶素-动量理论计算风机的气动载荷；在水动力计算方面，利用Airy线性波浪理论、势流理论以及Morison方程计算相关的水动力载荷；在锚泊模型方面，采用准静态的悬链线锚泊模型计算锚链张力；在多体动力学方面，利用Kane动力学方程构建系统的动力学方程，并求解各自由度的加速度信息；在柔性动力学方面，对于如塔架和桨叶等柔性体，采用高阶变形模型以考虑非线性的刚柔耦合效应；对比测试DARwind程序与FAST软件，以验证DARwind程序在海上浮式风机数值仿真方面的有效性.本文可为海上浮式风机耦合动力理论及其相关的数值仿真程序的研发提供一定的参考依据.

1 理论模型

1.1 水动力模型

DARwind程序中，利用Airy线性波浪理论、势流理论以及Morison方程计算相关水动力载荷[13].部分水动力参数，如：复原力系数、波浪激励力传递函数、波浪辐射力传递函数、二阶波浪力的二次传递函数(QTF)矩阵均由三维频域势流软件，如WAMIT、Sesam等，进行频域计算并生成.这些频域水动力参数将作为DARwind程序的读入数据，再依据实际计算的波幅信息、平台的位移、速度和加速度等参数，实时生成时域水动力载荷.其中，黏性阻力采用Morison经验公式进行修正.

波浪激励力FW包含波浪入射力和绕射力，考虑线性化的一阶波浪力，并采用谐波叠加法进行计算，即

(1)

式中：上标W为波浪激励；fn为第n个单元规则波的频率；εn为第n个单元规则波的相位；An为第n个单元规则波的波幅；HW(ifn)表示频率为fn的波浪力的频率响应函数；t为时间.

考虑自由表面记忆效应，波浪辐射力采用间接时域计算方法，利用频域势流理论计算可以获得辐射水动力参数μjk、λjk，再经过下面的公式，可以得到时域辐射力

(2)

j=1,2,…,6

(3)

沿湿表面对平台在水中受到的静压力积分可以获得浮体的静水力载荷，分为平台正浮时的浮力和平台发生位移时的静水回复力，其中静水回复力的表达式为

(4)

对于平台的细长构件，由流体黏性引起的载荷较为明显，势流理论无法准确地计算该载荷.此时，可根据Morison经验公式对势流载荷进行黏性修正.圆柱体在lx处离散的切片长度为dl的流体黏性力FV表达式为

(5)

式中：上标V为黏性；CD为黏性系数，一般由实验测得；D为圆柱直径；dl为圆柱离散的切片厚度；νw为水质点在切片处垂直于切片轴向的速度分量，其值等于波浪运动导致的水质点的速度与海流引起的速度在该方向的速度分量之和；νs为切片的运动速度在垂直于切片轴向的速度分量.黏性力矩阵包含集中力和力矩，沿着圆柱总长度LC进行积分，表达式为

(6)

因此，总的水动力FH根据线性理论，为上述FW、FR、FS和FV之和，则有

FH=FW+FR+FS+FV

(7)

1.2 锚泊动力模型

锚泊系统对海上浮式风机平台起到基本的定位和约束作用，同时影响着浮式结构物的运动响应特性.DARwind程序采用准静态的悬链线锚泊模型[14]计算锚泊对平台的约束载荷.其中，除锚泊的伸长变形外的一些动态特性，如惯性力、阻尼力、锚链变形弯矩均被忽略不计.

图1 悬链线锚泊示意图Fig.1 The schematic diagram of catenary mooring lines

悬链线锚泊模型如图1所示.其中，锚点受到水平方向张力HA和垂直方向张力VA的作用.假定躺底段(锚链与海床接触的部分)长为lb；锚链任意一点的张力为Ts′，水平方向的张力为Hs′，垂直方向的张力为Vs′，其距离锚点的水平距离为xs′，垂直距离为zs′，从躺底点算起的锚链自然长度为ls′；导缆孔的张力为TF，该张力水平方向的分力为HF，垂直方向的分力为VF；导缆孔到海底的垂直距离为zF，锚链整体的水平投影长度(水平距离)为xF；锚链的总长度为L，湿密度为w, 截面面积为A，轴向刚度为EA.

当锚链与海底接触时，即lb≠0，锚点处此时没有垂向张力作用于锚链.根据导缆孔和锚点之间的水平距离和垂向距离，平台导缆孔处对锚链的水平张力和垂向张力有如下关系式：

(8)

(9)

当lb=0时，锚点对锚链的垂向张力等于导缆孔对锚链的垂向张力减去锚泊质量，关系式为

VA=VF-wL

(10)

此时，导缆孔处的水平张力和垂向张力的关系式为

(11)

(12)

只要获得锚链当前时刻的锚点和导缆孔的位置(一般情况下，锚点的位置保持不变，导缆孔的位置需要根据当前的平台位置信息获得)，通过上述非线性方程组即可求得导缆孔对平台的约束张力.其中，非线性方程的数值求解可以借助迭代技术，如牛顿迭代法等进行数值求解.

1.3 空气动力模型

叶素-动量法[15]结合了一维的动量理论和二维的叶素理论.动量理论是基于理想的无穷桨叶正对来流时的致动盘假定推导而获得的.叶素理论则考虑了桨叶的实际形状，将桨叶划分为许多二维的翼型剖面，称之为“叶素”.根据每一个叶素上实际入流风速计算获得局部入流攻角，然后根据入流攻角从翼型气动系数的数据库中插值获得当前翼型剖面的气动系数，如升力系数、阻力系数和转矩系数等.根据这些系数，结合相对入流风速即可计算获得该剖面的气动载荷.最后，再通过积分(叠加)各个桨叶剖面叶素的气动载荷，计算总的风机气动载荷.令由动量理论得到的气动载荷和由叶素理论得到的气动载荷在环形控制体(见图2(a))内保持相等，即可得到各不同径长的环形控制体内的桨叶叶素的诱导系数，其表达式为

(13)

(14)

式中：a为轴向诱导系数；a′为切向诱导系数；σ=crbBblade/(2πrb)为风轮实度，其中Bblade为桨叶数，rb为叶素的径向长度，crb为叶素局部弦长；φ为入流角；CL和CD分别是叶素升力和阻力系数.

图2 气动模型示意图Fig.2 The schematic diagram of aerodynamic model

在经典的叶素动量理论中，假定叶轮只有转动，且来流是空间均匀和时间定常的.但对于实际的海上浮式风机而言，空气流动是时空变化的，浮式平台的运动也会引起桨叶的运动，加剧其相对入流风速的变化，气动诱导系数也随时空所变化，因此必须考虑非定常的局部气动模型，如图2所示.其中，α为叶素的局部入流攻角.为了简化模型，各叶素的气动性能只考虑小范围控制体内的影响.其研究的单元气动影响区域为一个扇形区域，如图2(c)所示.该扇形区域面积是图2(a)所示的环形控制体区域面积的1/3，即在当前时刻每个桨叶的叶素气动计算只考虑该桨叶附近区域的影响.

当不考虑桨叶柔性变形和平台运动所带来的叶素运动，叶素的相对入流风速如图2(b) 所示，其表达式为

(15)

式中：v0为上游未受干扰的入流风速大小；Ω为风轮的转度；vrel为最终合成的相对入流风速矢量；在非定常的局部叶素-动量模型中，DARwind程序考虑了局部桨叶的叶素因平台运动、塔架变形等引起的运动速度，以及由桨叶本身变形导致的叶素运动速度.叶素-动量法虽然计算速度快，且过去一直作为风能领域常用的可靠计算方法，但由于其并不能完全描述风场及风机的细节信息，且计算精度在特定场合可能有一定的误差，故需要引入一些修正方法[16]，例如：通常情况下，风机桨盘面并不正对来流风，此时称之为“偏航”状态.对于海上浮式风机而言，由于浮式平台的大幅度运动，这种偏航状态更为常见，故需要考虑偏航状态下的气动修正；另外，由于动量理论假定叶片数量为无穷多，但实际的风轮叶片数量是有限的，风机尾流涡系与无穷桨叶的理想风轮的尾流场有差别，此时需要采用Prandtl叶尖损失模型进行修正；当轴向诱导因子约大于0.4时，简单的动量理论不再适用，此时需要使用 Glarert 经验公式[16-17]进行适当的修正，以使数值计算结果与实际结果更为符合.上述所提及的修正理论均被采用并编入DARwind程序的气动模型中，程序中的叶素-动量法是一个反复迭代求解诱导系数直至收敛的过程，其计算速度关系到程序整体的仿真速度，所以此模块还需要采用一些算法的优化技术以提高迭代计算的收敛速度.

1.4 动力学模型

海上浮式风机的多个不同部件拥有不同的特性，可各自根据相对变形量简化为刚性体或柔性体，因此海上浮式风机模型为一个刚柔耦合的多体系统.采用浮动坐标系[18]描述刚柔耦合多体系统的运动关系.该方法采用两套坐标系描述每个部件的位置、方向及变形：① 全局坐标系，用以描述各部件的位置和方向；② 各部件的随体坐标系，用以描述部件的柔性变形场.其中对于柔性体，如塔架和桨叶，其变形场采用模态叠加法进行有限模态的近似离散，变形量U如下所示：

U=ΘQ

(16)

式中：Θ为模态形状函数矩阵；Q为模态坐标(广义坐标)矩阵.研究表明[19-20]：当柔性体进行大范围运动时，尤其是高速旋转运动，用传统的线性化方法将会产生较大的数值误差.因此在DARwind程序中，柔性体(塔架和桨叶)均被近似为三维的Euler-Bernoulli 模型.当发生大范围的运动时，梁的侧向变形会导致轴向的耦合收缩效应，并且随着运动频率的增加，柔性体的刚度将会随之增加，该效应称为刚柔混合多体的“动力刚化”.三维梁模型的高阶非线性耦合变形量可简化为

(17)

式中：ush是由梁结构侧向变形引起的轴向收缩量；uy0和uz0分别是梁侧向变形量关于局部随体坐标系y，z方向的分量.当考虑高阶非线性耦合变形效应时，上述的柔性变形可以表示为

(18)

其中非线性高阶耦合效应矩阵H为

(19)

(20)

式中：下标x为对应物理量的坐标分量.通过上述的浮动坐标系法及模态离散，可以描述系统各部件的各个质点之间的运动学模型.研究中需要构建系统的动力学控制方程，以求解系统在外力和系统约束作用下各自由度的信息.采用Kane动力学方程[21]构建系统动力学控制方程.这种方法兼具矢量力学和分析力学特点，既避免了出现理想约束反力，也避免了对动力学函数进行求导，减少了计算工作量，且易于拓展和推广.

一个质点系统，关于第j个广义速度的Kane方程可表示为

Fj+F*j=0

(21)

广义主动力Fj以及广义惯性力F*j分别计算如下：

(24)

(25)

式中：Bki和Dki分别为部件i的平动和转动的加速度余项，其解析表达式较为复杂，可通过线速度和角速度关于时间的一阶导数推导而得，需要注意各个部件的坐标系之间的转换关系.

对于刚柔伺服多体耦合的海上浮式风机系统，其Kane方程可表示为

(28)

Nr+Nf+Ns=N

将式(28)中的各部件的广义主动力和惯性力分别进行累计计算，则有

(29)

(30)

式中：Fi和Mi分别为部件i的动力学参考点的力和力矩；mi和Ii分别为部件i的质量和转动惯量；式(29)右端第一项为刚体和柔性体的广义主动力，比起刚体，柔性体中存在额外的结构变形能Ei；式(29)右端的第二项为控制器的广义主动力或内力；式(30)右端为刚体和柔性体的广义惯性力和力矩.在海上浮式风机系统中，浮式支撑平台、机舱、转轴和轮毂结构均视为刚体，塔架和桨叶均视为柔性体.需将柔性体划分为许多微小单元，对每个微单元分别建立Kane方程，故柔性体的存在增加了动力学方程的复杂度和计算量.

[IiDi+ωi×(Iiωi)]}=0

(31)

令

(32)

[IiDi+ωi×(Iiωi)]}

(33)

(34)

1.5 控制系统模型

DARwind程序采用发电机转矩控制器和桨距角控制器.当风速小于额定风速时，以优化风能捕获效率为目标，通过调节发电机转矩进而调节风机的最佳叶尖速比；当风速达到或超过额定风速后，则主要通过改变桨距角控制系统调节气动性能，进而达到稳定功率的目的.

图3 发电机转矩控制Fig.3 Generator-torque control

桨距角控制器采用比例-积分统一变桨机制，其表达式为

(35)

式中：θ为桨距角；Δθ为当前需要改变的桨距角度；t为桨距角控制器的执行时间；ΔΩ为当前风轮转速与目标转速之差；KP=KPc/[1+(θ/θc)]为桨距角控制器的比例项，其中θc为当桨距角为此值时，功率对桨距角的导数为桨距角等于0时所对应的导数的2倍，KPc为θc桨距角下的比例项；KI=KIc/ [1+(θ/θc)]为桨距角控制器的积分项，KIc为θc桨距角下的积分项.二者的具体取值需要根据实时的桨距角动态进行调整，其推导过程可参考文献[22].

1.6 程序计算流程

运用上述理论模型构建用于海上浮式风机动力性能预报的DARwind程序，计算流程如下：

(1) 程序开始时，读入标记参数、结构信息及环境信息等，用以初始化计算模型.

(2) 在每一个计算时间步下，程序需要根据平台的位置、速度信息及平台参考点处的波浪场信息生成相应的水动力载荷；根据桨叶的变形、运动及局部桨叶在当前时空的风场风速，利用非稳态叶素-动量法计算气动载荷；根据平台的位置信息计算导缆孔的位置；利用悬链线锚泊模型计算当前锚链对平台的约束力.另外，程序还考虑了重力、浮力及控制器的作用力等参数.

(3) 经过步骤(2)求得所有环境载荷后，通过Kane方程构建系统动力学控制方程，求得当前各自由度的加速度信息，并利用4阶Runge-Kutta时间积分法递推下一时刻各自由度的位置及速度信息.

(4) 程序依次循环直到预设的结束时间.

相比于其他类似程序，DARwind程序有两方面的创新[23]：一方面，该程序可以非常简便地根据不同计算目的，对海上浮式风机系统建立不同的多体系统方案.例如：当需要快速获得海上浮式风机系统的动力特性时，可将整个系统作为单体模型；当需要适中的计算速度和精度时，可以将整个系统作为多刚体结构，忽略柔性结构的影响；当需要较为准确的系统模型时，可将海上浮式风机系统作为刚柔耦合多体模型.另一方面，DARwind程序中创新地采用了高阶的刚柔耦合模型，考虑了刚体运动对柔性体动力性能的影响.

2 对比验证

DARwind程序将与美国国家可再生能源实验室(NREL)的风机仿真FAST 程序[6]进行对比分析，以验证DARwind程序在海上浮式风机数值仿真方面的有效性.采用OC3-Hywind Spar海上浮式风机作为研究对象，如图4所示，关于该海上浮式风机更详细的信息可以参考文献[24].

图4 OC3-Hywind Spar海上浮式风机[24]Fig.4 OC3-Hywind Spar floating offshore wind turbine[24]

2.1 固有频率测试

在相同的模型参数设置下，DARwind和FAST程序在自由衰减运动中测得的固有频率以及无因次阻尼系数如表1所示.由表1可知, 除了纵摇运动具有微小的差异之外，DARwind以及FAST程序的自由衰减运动的特性都是非常接近的.该结果表明DARwind程序能准确地反映海上浮式风机系统的质量、惯性、静水回复力、锚链刚度及水动力阻尼等特性.

表1 DARwind和FAST程序中的运动固有频率与无因次阻尼系数

Tab.1 Natural frequency and damping coefficient of platform motions in DARwind and FAST softwares

运动模态固有频率/HzDARwindFAST无因次阻尼系数DARwindFAST纵荡0.00810.00800.48770.5286垂荡0.03230.03230.24230.2415纵摇0.03140.03340.27590.3300首摇0.12180.12090.27230.2820

2.2 风工况

海上浮式风机是用以捕获风能的离岸设备，对其气动性能及结构变形的预测对海上浮式风机而言同样也是需要重点考虑的因素.对DARwind程序进行稳态风工况测试，测试结果如表2所示.数值计算中，风速空间均匀且时间定常；转速为风轮转速；桨距角设置为0°.当不考虑控制系统的作用且平台固定时，该工况的测试结果如表3所示.由表3可知，DARwind 程序的气动载荷值(气动推力和气动转矩)比FAST程序的气动载荷值略大，这可能是由于两个数值软件在气动理论方法上的差异所导致的.DARwind程序采用的是基于叶素扇形影响区域的局部非稳态气动模型；而FAST程序则采用的是基于风轮整体的稳态气动模型.在塔架和桨叶的变形方面，相比于FAST程序，DARwind程序的计算结果相对较小，但两者差别不大.这主要是由于DARwind程序采用了高阶的刚柔耦合模型(见式(17)～(20))，考虑了刚体运动对柔性体的影响.当柔性体发生大范围运动(如高速转动)时，柔性体的刚度会随之提高，变形也会随之减小.总体而言，无论是气动性能还是结构弹性动力特性，两个程序的计算结果都比较接近，这验证了上述柔性体变形理论及叶素-动量理论在海上漂浮式风机时域计算的有效性.

表2 稳态风工况测试Tab.2 Steady wind cases

表3 气动载荷与叶尖变形量Tab.3 Aerodynamic loads and blade-tip deflection

2.3 风波联合工况

海上浮式风机工作于有风及波浪载荷共同作用的海洋环境下，因此，对浮式平台运动特性的预报显得尤为重要.下文将验证DARwind程序在风波联合工况下的运动特性.测试风速为额定风速(11.4 m/s)；风轮保持额定转速(12.1 r/min)；桨距角统一设置为0°；不考虑控制器的作用；波浪模型参考“北海联合海浪计划”(Jonswap)波浪谱，有义波高Hs为2.0 m；谱峰周期Tp为8.0 s；谱形参数r设置为3.3.DARwind 和 FAST 程序的纵荡及纵摇的功率谱密度对比如图5所示.由图5可知，两种数值程序在风波联合工况下的运动特性吻合得较好，且均能反映出纵荡及纵摇运动中明显的共振响应、波浪能量区的强迫运动响应、纵荡及纵摇的耦合运动效应.

图5 平台运动功率谱密度Fig.5 Power spectral density of platform motion

2.4 控制系统模型

在海上浮式风机正常运转的过程中，控制系统起着重要的作用，不但能够调节海上浮式风机的气动载荷与发电功率，而且对海上浮式风机的运动响应有着重要的影响.在变桨距控制过程中，通过桨叶的桨距角、风机功率等参数的变化情况，验证 DARwind 程序控制系统的有效性.测试中，海上浮式风机的平台受载荷运动的影响，对控制系统的要求比较高，有别于传统的固定式风机.在该控制系统下，桨距角随时间的变化情况如图6(a)所示，发电机功率随时间的变化情况如图6(b)所示.由图6(a)可知，上述测试过程中的桨距角随时间的变化呈现出波动，这是由于平台运动导致桨叶的相对入流风速急剧变化导致的.此时，为了保持发电机功率的稳定，桨距角也需要进行相应的调整.对比这两个程序的计算结果，总体上DARwind程序的桨距角变化波动稍微大一些，但是两者的变化情况基本相同.图6(b)中的发电机功率也有类似的情况.其主要原因是由于在额定工况时，FAST程序的转矩控制依然在进行调整，而DARwind程序在额定工况下选择稳定发电机转矩，以减少轴系疲劳损伤.此时，变桨距控制器成为了DARwind程序唯一的发电机功率调节手段，所以其桨距角的波动会相对大一些.

图6 控制系统下的桨距角及发电机功率随时间的变化情况Fig.6 Blade pitch angle and generator power vary with time

3 结语

本文介绍了海上浮式风机时域耦合仿真程序——DARwind程序的相关理论，并做了相应的数值验证分析.研究结果表明：所采用的理论方法在海上浮式风机时域耦合计算中是有效的，DARwind程序能够有效地仿真并准确地预报海上浮式风机复杂的时域耦合动力响应特性.然而，当前的DARwind程序仍有进一步发展和优化的空间，如：可以进一步提升气动模型精度，尤其对于海上浮式风机大范围摇摆运动导致的复杂尾流场而言；当前的控制系统也可以采用更为高级的控制方式，如神经网络控制、模糊控制等现代智能化控制方式.