深水SPAR风机系统全耦合动力响应分析研究

闫发锁，门骥远，彭 成

（1.哈尔滨工程大学 船舶工程学院，哈尔滨 150001；2.美国德州农工大学 土木工程系，大学城 77843）

深水SPAR风机系统全耦合动力响应分析研究

闫发锁1，门骥远1，彭 成2

（1.哈尔滨工程大学 船舶工程学院，哈尔滨 150001；2.美国德州农工大学 土木工程系，大学城 77843）

文章采用联合开发的计算程序对深水SPAR风机的浮体、锚泊和风机各子系统进行了水—气动力的全耦合数值分析，研究了深水浮式风机系统的动力响应特点。浮体水动力计算采用基于二阶精度的混合波浪模型（Hybrid Wave Model）的MORISON公式，锚泊系统采用细长杆理论通过非线性有限元方法实现，风机系统的空气动力分析采用基于多体气动弹性理论的FAST模块。以浮体控制方程为主体，通过模块间的载荷与位移传递在每个时间步上迭代求解，形成完全耦合的时域分析方法。通过对NREL的5MW SPAR风机系统在随机海况下的水动力响应分析，验证了该方法的有效性，并分析了浮式风机子系统间的混合动力作用。

浮式风机；动力响应；混合波浪模型；系统耦合

0 引 言

海洋风能储量丰富，在倡导绿色环保、节能减排的可持续发展的全球背景下，海上风电的开发利用是资源环境与社会可持续发展模式的典型示范，已成为国际可再生能源开发的主力方向。由于技术和成本的限制，直至上世纪90年代海上风电场的建设才取得了长足的发展。海上漂浮风机系统由叶轮、塔架、浮体和锚泊系统构成，是深水风电开发的基础装备。目前在深水浮式风机技术方面，欧美和日本等发达国家已积累了丰富的理论成果和实践经验，完成了多个概念设计和试验样机，其中影响较大的有SWAY、HYWIND、WindSea AS、Blue H和WindFloat等（图1）。

图1 几种海上浮式风机系统（由左至右分别为：SWAY-2011、HYWIND-2009、WindSea AS-2008、Blue H-2008和WindFloat-2008）Fig.1 Some types of offshore floatingwind turbine systems

与固定式风力机相比，漂浮式风力机浮式基础在海洋环境载荷作用下有一定幅度的运动。风力机位于距海面近百米的高空，浮式基础的微幅运动即可造成风力机姿态较显著的变化，对风力机气动载荷和功率输出会产生显著影响。因此，研究浮式风机系统的运动响应对发挥海上风力机的性能特点，进行合理的设计优化具有重要的指导意义。

风机叶轮的气动载荷是浮式风机系统重要的环境激励。风机气动力性能的分析目前主要有两种方法：传统的动量叶素理论（BEM）和基于CFD的涡面元方法[1]。BEM方法基于二维假设，理论模型简洁直观、计算快速，但叶素单元的受力计算需要合理的试验或经验系数和端部修正。但从实用意义考虑，针对其局限性已发展了丰富的修正处理方法，BEM仍然是风机叶片设计载荷的主要计算手段。Gutierrez等[2]采用基于BEM的FAST软件开发了浮式风力机的优化设计工具FASTLognoter，具备风力机的参数设计能力。涡面元法考虑了叶片流场旋涡尾流特性，理论模型精确，但数值实现技术难度高，计算量大，目前还主要处于研究阶段。此外，气动力分析方法还要结合风机控制系统的控制模式和叶轮、塔架等的结构特征，才能全面考虑系统的动力性质。基于涡面元法，Sebastian等[3]发展了自由涡放尾流（FVM）方法进行了风力机的气动响应分析。该方法相对BEM在叶片尾流演变和气动载荷方面体现了较高的精度，但计算资源耗费较大。

借助于大型海洋平台的设计分析技术，海上风机系统中浮体的选型与优化、锚泊系统设计以及二者的耦合水动力响应分析已具备一定的基础。Jonkman等[4-5]基于FAST对已有的半潜式、张力腿式和SPAR式三种类型风机浮体的水动力进行了分析。唐友刚等[6]针对一座5MW海上风机设计了SPAR浮体和三立柱—撑杆形式的半潜浮体，并进行了性能和结构分析。Sethuraman等[7]对一个阶梯型SPAR平台风机的水动力响应进行了模型试验和数值预报。Karimirad等[8]采用两种分析代码针对一个单柱张力筋系泊的浮式风机系统进行了平台运动和系泊张力的水弹性时域模拟，并对系统耦合进行了模型试验研究。

综上所述，海上风机各组成系统综合了水、空气结构之间的多域动力耦合分析技术。以浮式平台的相关理论为基础，利用现有浮体水动力分析方法结合风力机械的气动力载荷模块，进行浮式风机系统的动力响应分析成为合理可行的研究方式。目前已有的分析工具包括业界具有较大影响的DeepC、OrcaFlex、Charm3D等水动力分析程序和FAST、GH Bladed、TurbSim等风机气动分析模块[9-11]。针对风机—浮体—系泊三者之间的气动力与水动力的多域耦合，目前的做法需要采用各自独立的分析模块在时域上迭代耦合，分析结果和计算效率也有待检验和提高。本文将水动力耦合分析程序与现有的风机气动力程序结合，完成了适用于深水浮式风机系统的全耦合时域动力响应分析。浮体水动力的分析采用具有二阶精度的混合波浪模型（Hybrid Wave Model）与基于细长体假设的MORISON方法；锚泊系统采用细长杆理论通过非线性有限元方法实现；风机系统的空气动力分析采用基于BEM和多体气动弹性理论的FAST模块。通过对NREL 5MW Hywind OC3风机动力响应的数值模拟与现有结果的比较，分析了浮式风机子系统间的混合动力作用。

1 数值模型

1.1 风机气动载荷分析

系统中风力机的气动载荷采用了FAST[12]（Fatigue,Aerodynamics,Structures,Turbulence）。FAST是由美国国家可再生能源实验室（NREL）开发的一个开放资源程序，用以模拟传统两叶或三叶水平轴风机的气动弹性载荷。FAST中计算风载荷的主要模块是AeroDyn，其主要采用了经典的动量叶素理论[13]。

基于动量叶素理论，风轮气动性能的计算

式中：Q和T分别为作用在风轮上的扭矩和盘面推力；ρ为空气密度；CL和CD分别为升力和阻力系数；φ为入流角；RH和R分别为轮毂及叶片半径；B为叶片数目；c为叶片单元弦长；r为叶片单元的展向半径；Vr为叶片单元处的气流相对速度，其中U为远方来流风速，a和a′分别

r为轴向和周向诱导因子；F为叶片端部尾流的修正系数，源于Prandtl的气动理论方法[13]，其中

1.2 浮体动力分析

采用图2（a）所示坐标系，对某一刚体其上任意一点的6自由度非线性运动方程可以用（3）式和（4）式表示。

图2 分析所采用的坐标系统Fig.2 Coordinate systems for theorymodels

式中：上标‘t’表示转换矩阵，m为结构质量，ξ=（ξ1,ξ2,ξ3）t为坐标系（即大地坐标系）o点坐标值，=d2ξ/ d t2为线性加速度，ω= （ω1,ω2,ω3）t为角速度，rg=（xg,yg,zg）t为中心坐标值，I0为刚体的惯性矩。此外，u、rg和I在随体坐标系即坐标系中定义。要注意到和M分别是作用于刚体上的总力和总力矩，后者是oxyz

00在坐标系中定义的，并且T是（oxyz）坐标系到（OXYZ）坐标系的转换矩阵，其为正交矩阵，正交变换与自身相同 （Tt=T-1）。

当柱体直径相对波长较小时，选取合理的阻力系数和附加质量系数，莫里森方程可较准确地用于计算柱体上的波浪载荷。沿着与柱体轴向垂直的方向n作用于柱体单元上的受力（惯性力和阻力）可以用下式表示：

式中：下标‘n’表示法向分量，Dj为等效柱体直径，zj为第j个单元的湿长，ρf为环境流体的密度，Cm为附加质量的系数，CD为阻力系数。vr和ar分别为环境流体相对于平台结构的速度和加速度，且an为柱体加速度。作用于某一穿过自由液面的柱体的轴向水动力，包括作用于单元底部面积SBj上的动压力，阻力和附加质量力。轴向力可以表示为：

式中：下表‘t’表示轴向分量，φ（1）和φ（2）为一阶和二阶入射波势，CDt和Cmt为垂向阻力系数及附加质量系数，且vrt和art是柱体相对于靠近柱体底部附近的环境流体相对速度和加速度的轴向分量。

1.3 系泊系统分析

当没有额外转矩和力矩作用于系泊线/立管，在图2（b）坐标系下运动的方程可以表示为

式中：r（s,t）为系泊线/立管的位移矢量，是弧长s和时间t的函数；B表示弯曲刚度；q表示单位长度外力，其可通过莫里森方程算得；r为单位长度质量；λ=T-Bκ2为一个标量变量，其中κ为系泊线/立管的曲率，而T为局部张力。“’”表示对s的导数，“·”表示对时间t的导数=∂2t/∂t2为局部加速度。此外，对于不计轴向拉伸的缆索，r必须满足几何约束方程

作用于杆上的外力包括重力、流体静水力和流体动水力。作用于杆上的流体动水力包括附加质量力、阻力和Froude-Krylov力。方程（7）可以通过非线性FEM法数值求解。使用Spar及其系泊系统间的配置条件，柱体及系泊系统的动力方程可以通过Newmark-β方法在时域上同时求解。

2 风机系统参数

为了验证浮体与系泊系统间的耦合作用，本文选取NREL 5MW OC3 SPAR[14]风机进行了系统动力响应分析，并将计算结果与现有的三维势流理论软件分析结果进行了比较。

风机系统包括上部传统的三叶片变螺距风力机及塔架与下部Hywind-OC3 Spar平台及锚泊两部分。Hywind-OC3 Spar风机设计水深320m，系泊系统采用3根系泊线，它们在静水面投影之间的夹角为120°。导缆器位置在静水面以下70m处，导缆口距平台柱体中心线之间的半径为5.2m。海底锚点分布在系统静平衡时距SPAR中心线半径为853.9 m的圆周上。每根系泊线长度为902.2 m，直径0.09m。对系泊系统的模拟做了以下简化：系泊线与平台的三个连接点用铰链连接代替，艏摇方向加入弹簧以模拟原结构刚度；实际上多成分的系泊线近似为均匀湿重（单位长度）和刚度的锚缆。风机结构的主体尺度与系泊系统的布置如图3所示，系统总体的结构参数见表1。

图3 Hywind-OC3 5MW浮式风机系统组成Fig.3 Sketch of Hywind-OC3 5MW floating wind turbine

表1 Hyw ind-OC3 5MW浮式风机系统结构参数Tab.1 Characteristics of Hyw ind-OC3 5MW w ind turbine system

3 系统频域响应分析

为例验证各子系统理论方法及数值实现，首先对风机浮体进行频域水动力分析，得到自由浮体的线性频响函数；之后，考虑锚泊系统的影响，得到系统的水动力耦合频响函数。二者在分析时假设风机塔架和叶片是静止的，未计及上部风机的气动力载荷，但考虑了风机结构的质量及分布。图4为浮体系统纵荡、垂荡、纵摇和首摇4个自由度的波激运动频响函数。浮体在自由状态下，垂荡运动的固有周期为31.4 s.，纵摇/横摇为31 s。

图4 浮体系统的运动频响函数Fig.4 Motions frequency response function of floating body including turbine structure

SPAR平台系统通常使用半张紧式的锚泊系统，锚泊系统的影响主要体现在平动自由度上。基于自由衰减法进行风机系统的时域耦合分析，表2列出了自由衰减模拟的前五个周期的平均值。由表2可见，有了锚泊的约束，系统平面运动固有周期明显减少，垂向及相关运动变化不显著。为了便于验证比较，表中同时给出了文献[9]水动力耦合计算的结果。二者在横荡、纵荡和垂荡运动的周期几乎一致，而在纵摇、横摇和艏摇有较小差别。本文方法计算时较准确地考虑了锚泊系统的影响，与自由浮体相比，垂荡、纵摇固有周期有较小降低，纵摇、横摇因锚泊布置的差别稍有不同，艏摇周期相对有稍大差别。

表2 平台系统的运动固有周期Tab.2 Natural frequencies of platform motions

4 全耦合数值分析与比较

应用数值技术将浮式风机系统的水动力模块与气动弹性模块耦合，得到风机系统结构中心处的运动响应及谱密度，如图5所示。风机作业海域水深320m，分析环境条件选用JONSWAP表示的随机不规则海况，有义波高5.0m，峰值周期8.69 s，参数γ取值3.6，环境风速13m/s。风、浪的传播方向沿x轴正向。因为计算步长小（0.02 s）,且锚泊系统应用细长杆理论需要大量的迭代计算，计算耗时较长。为了合理控制时间和充分利用计算资源，分析采用了分段进行的方式。本例对于3 600 s的计算分三段进行，每段计算1 300 s。为了避免每个分段起始部分重复，各分段之间计算预留了150 s的共用重合区，该区间内应用适当方法进行了处理。为了清晰显示，图中的运动时历的时间为1 300 s，对应的功率谱密度（PSD）的统计是基于3 600 s的计算结果。

图5 风机系统全耦合运动响应与能量谱密度Fig.5 Motion responses and PSD by the full coupledmethod

由于全耦合分析计入了风力机叶片的气动弹性载荷，系统在风浪方向的运动（纵荡）具有较大的运动幅度和平均偏移。根据六自由度运动的时域统计，纵荡运动平均值为9.369 m，最大值为17.708 m，最小值为4.125m，标准差为1.747m。因为在耦合分析中，由于系统上部风力机的叶片的旋转，再加上SPAR主体纵摇速度影响叶片迎面风速，系统纵向受力幅值相比忽略风力作用（图5）时变化大，响应的运动幅值也较大。由此表明，风力机纵向载荷对水动力载荷在频域方面的影响是较大的。本文方法要考虑圆柱体绕流旋涡脱落导致的横向力，而且当风机工作时也会产生一个相当小的横向力。因此，尽管风和浪在x轴方向，主体仍具有一定的横荡和横摇。耦合分析得到的艏摇运动比较小。这实际上当主体迫使风机做艏摇偏转时，风机便可以通过监控机械来实现控制，这种作用于风机上主体首摇运动的反馈在分析时应该被考虑在内。

采用细长杆理论与浮体运动的耦合求解，得到的锚泊系统的受力如图6所示。图6（a）为三根系泊线的顶部张力，图6（b）为底部水平张力。每根系泊线的顶部张力平均比底端张力大约200 kN。由于对称性，理论上#1系泊线的张力应该与#2的张力相同。实际数值计算时，两根系泊线是独立迭代计算的，并且由于系统横向载荷和运动的影响，二者的张力在时域上有较小的局部差别，但总体的吻合度很好。#3的布置与#1和#2相反，其张力大小变化与后两者呈相反趋势，结果也得到了很好的体现。

图6 锚泊线张力Fig.6 Tensions inmooring lines

全耦合在上部风力机气动载荷计算时，考虑了浮体SPAR的运动状态。已有的部分研究在计算风力时假设浮体SPAR是固定的，即认为浮体的姿态和运动对风力机的影响较小。通过将本文全耦合方法与后者假设情况下，在上述随机海况下进行计算分析对比表明，系统的垂荡响应具有较明显的差别（图7）。计算气动力载荷时不考虑浮体的运动状态，系统的最大垂荡响应接近全耦合计算时的2倍。因为系统的纵摇的自然周期（30.7 s）非常接近于升沉运动固有周期（30.9 s），风机所产生的外部纵摇力矩与垂向波浪力的耦合加剧了垂向运动。而在全耦合分析时考虑了SPAR浮体的每时刻的运动和空间姿态，垂荡和纵摇运动的固有频率产生了偏移，二者的谐振效应明显减弱。

图7 全耦合与近似耦合垂荡响应的比较Fig.7 Comparison of heavemotion between full coupling and simple coupling

5 结 论

本文研究了海上浮式风机系统全时域耦合分析方法，将水动力载荷模块与空气动力载荷模块进行联合求解，并对NRELHywind-OC3 5MW风机系统的运动响应进行了全耦合分析。计算结果表明，全耦合方法有效地对风机系统的自由浮体和具有张紧式系泊状态下的运动RAO进行分析。在随机海况下，由于全耦合分析计入了风力机叶片的气动弹性载荷，系统在风浪方向的运动（纵荡）具有较大的运动幅度和平均偏移。由于全耦合考虑了风力机迎风速度的变化，全耦合时系统的动力响应与不考虑浮体运动状态的近似耦合相比，某些自由度的运动可能会有较大差别。针对本文算例所采用SPAR类型平台，因为系统的垂荡和纵摇固有频率相近，两种方法得到的垂荡运动在时域和频域结果有显著的差别。本文的研究表明，针对浮式风机系统发展全耦合分析方法是必要的。

[1]Hansen M O L,Sorensen JN,Voutsinasb S,Sorensen N,Madsen H Aa.State of the art in wind turbine aerodynamics and aeroelasticity[J].Progress in Aerospace Sciences,2006,42:285-330.

[2]Gutierrez JE,Zamora B,García J,Peyrau M R.Tool development based on FAST for performing design optimization of offshorewind turbines:FASTLognoter[J].Renewable Energy,2013,55:69-78.

[3]Sebastian T,Lackner M A.Developmentof a free vortex wakemethod code for offshore floating wind Turbines[J].Renewable Energy,2012,46:269-275.

[4]Jonkman JM,Butterfield S,MusialW,Scott G.Definition of a 5-MW reference wind turbine for offshore system development[R].NREL/TP-500-38060,National Renewable Energy Laboratory.Colorado,USA,2009.

[5]Jonkman JM,Matha D.Dynamics of offshore floating wind turbines-analysis of three concepts[J].Wind Energy,2011, 14:557-569.

[6]唐友刚,王 涵,陶海成,刘中柏.海上风机半潜型浮式基础结构设计及整体强度分析[J].中国造船,2013,54(3):85-93.Tang Yougang,Wang Han,Tao Haicheng,Liu Zhongbai.Structure design and global strength analysis for semi-submersible floating foundation of offshorewind turbine[J].Shipbuilding of China,2013,54(3):85-93.

[7]Sethuraman L,Venugopal V.Hydrodynamic response of a stepped-spar floating wind turbine:Numericalmodelling and tank testing[J].Renewable Energy,2013,52:160-174.

[8]Karimirad M,Meissonnier Q,Gao Zhen,Moan T.Hydroelastic code-to-code comparison for a tension leg spar-type floating wind turbine[J].Marine Structures,2011,24:412-435.

[9]Bae Y H,Kim M H,Im SW,et al.Aero-Elastic-Control-Floater-Mooring coupled dynamic analysis of floating offshore wind turbines[C]//International Offshore and Polar Engineering Conference,June 19-24,2011.Maui,Hawaii,USA,2011.

[10]Masciola M,Robertson A,Jonkman JM,Driscoll F.Investigation of a FAST-OrcaFlex couplingmodule for integrating turbine andmooring dynamics of offshore floatingwind turbines[C]//International Conference on OffshoreWind Energy and Ocean Energy.Beijing,China,2011.

[11]Karimirad M.Dynamic response of floatingwind turbine[J].Mechanical Engineering,2010,17(2):146-156.

[12]Jonkman JM.,Jr,BuhlM L.FASTUser’s Guid[R].NREL/EL-500-38230,2005.

[13]Moriarty P J,Colorado G,Craig Hansen A.AeroDyn Theory Manual[R].NREL/EL-500-36881.National Renewable Energy Laboratory.Colorado,USA,2005.

[14]Jonkman J.Definition of the floating system for phase IV of OC3[R].NREL/TP-500-47535.National Renewable Energy Laboratory.Colorado,USA,2010.

Study on full coup ling dynam ic responses of a deepwater SPAR w ind turbine system

YAN Fa-suo1,MEN Ji-yuan1,PENG Cheng2
(1.Harbin Engineering University,Harbin 150001,China;2.Texas A&M University,College Station 77843,USA)

Themethod to perform coupling analysis for offshore floating wind turbines(OFWT)is developted through integrating hydrodynamic and aerodynamicmodules in time domain.Morisonmethod is used for hydrodynamic computation of floating body and itsmooring system,in which relative velocity between structure elements and waves is implemented by the Hybrid Wave Modelwith second order accuracy.Slender rods theories are applied on themooring systems,and the aerodynamics load induced by wind turbine is simulated by NREL’s code-FAST.Loads and displacements are transferred between the submodules based mainly on floating body control equations in every time step by Newmark-βmethod.Motion response of a 5MW 3 blades spar type OFWT is predicted with and without FAST to validate the combined program.A comparison with results from available 3D linear potential flow method in a random sea condition shows that the code is capable of hydro-aero dynamic analysis for OFWT.

OFWT;dynamic responses;hybrid wavemodel;subsystems coupling analysis

P752 TK89

：Adoi:10.3969/j.issn.1007-7294.2017.02.005

2016-09-30

国家自然科学基金项目（51379051，50909022）；上海交通大学海洋工程国家重点实验室开放基金（1208）；国家国际科技合作专向项目（2013DFE73060）

闫发锁（1977-），男，博士，副教授，E-mail:yanfasuo@hrbeu.edu.cn；

门骥远（1990-），男，博士研究生。

1007-7294（2017）02-159-09