DOI:10.11990/jheu.201510073
, , , 2. 中海石油(中国)有限公司 天津分公司, 天津 300452;
3. 日本九州大学 应用力学研究所, 福冈县 春日 816-8580
, , , 2. CNOOC China Ltd.Tianjin 300452, China;
3. Research Institute for Applied Mechanics, Kyushu University, Fukuoka 816-8580, Japan
近年来，能源开发与环境保护的需求强力推动着风能科学与技术的研究，产业规模和新技术不断发展，2009~2014年全球风电装机容量的年均增长率高达18.47%[1]。
海上风能因分布广泛、发电稳定等优势而受到重视，海上风场建设和装机容量年均增长率大于风电总装机增长率，发展十分迅速[1-3]。海上风力机按照支撑结构可分为固定式和浮式两种，其中浮式风力机(floating offshore wind turbine, FOWT)可作业于深水海域，是海上风力机的主要发展方向[3-5]。根据FOWT基础的形式，可进一步分为Spar式、张力腿式以及半潜式等。其中，Spar平台是浮式风力机支撑结构中应用较广泛的形式，其结构设计可使整个系统重心降低[6]，在各种工况下[7-8]仍能保持良好稳性；另外，Spar式风力机在风浪环境下载荷与运动较小，经济性比张力腿式有很大的优势[7]。目前对于Spar型FOWTs的研究已经取得了显著成果，在运动响应、气动性能等方面发展了一些理论模型，在线性频域内进行系统性能预报[7-8]，并且开展相关的试验技术研究。2009年，世界上第一座漂浮式风力机Hywind在挪威附近的北海正式启用。然而，FOWTs在风浪中的运动是耦合的、非线性的。就Spar式平台本身已发展了非线性理论预报方法[9-11]。Hanson等[16-18]把HAWC2和SIMO/RIFLEX集成在一起，计算了FOWT平台在各种风浪工况下的性能，并与模型试验结果对比。Shim等用WAMIT计算NREL张力腿式风力机系统的水动力性能，用有限元方法分析锚泊系统载荷和强度，对风力机、浮体和锚泊系统耦合运动进行频域、时域分析。但对于平台上安装大型风力机的FOWTs，在风浪流载荷下的耦合运动研究主要限于线性简化模型[12]，其耦合特性规律和机理并不明确。本文对Spar型海上FOWTs建立时域非线性运动模型，应用时域耦合以及时频转换的方法研究定常风、规则波对系统运动响应的影响以及该系统在随机风浪下的运动响应特性。
1 计算方法考虑Spar型FOWTs(含风力机组、浮式平台及其锚系)在风浪流作用下的运动。假设平台为刚体，平台瞬时姿态由平台欧拉角Θ=(α, β, γ)T、重心瞬时位置X0=(x, y, z)T、重心运动线速度U0=(ux, uy, uz)T及角速度Ω=(ωx, ωy, ωz)T来描述。矢量形式的系统非线性动力学方程[10]组写为
| $\left\{ \begin{array}{l} \boldsymbol{M} \cdot \frac{{{\rm{d}}{\boldsymbol{U}_0}}}{{{\rm{d}}t}} = \boldsymbol{F}\\ \boldsymbol{I} \cdot \frac{{{\rm{d}}\mathit{\pmb{\Omega}} }}{{{\rm{d}}t}} + \mathit{\pmb{\Omega}} \times \left( {\boldsymbol{I} \cdot \mathit{\pmb{\Omega}} } \right) = N\\ \frac{{{\rm{d}}{\boldsymbol{X}_0}}}{{{\rm{d}}t}} = {\boldsymbol{U}_0}\\ \boldsymbol{B} \cdot \frac{{{\rm{d}}\mathit{\pmb{\Theta}} }}{{{\rm{d}}t}} = \mathit{\pmb{\Omega}} \end{array} \right.$ | (1) |
式中：M、I分别为系统的质量和转动惯量矩阵；F、N分别为作用于重心处的合外力和力矩，表示为
| $\left\{ \begin{array}{l} \boldsymbol{F} = {\boldsymbol{F}_{{\rm{hydrostatic}}}} + {\boldsymbol{F}_{{\rm{wave}}}} + {\boldsymbol{F}_{{\rm{current}}}} + {\boldsymbol{F}_{{\rm{wind}}}} + {\boldsymbol{F}_{{\rm{mooring}}}}\\ \boldsymbol{N} = {\boldsymbol{N}_{{\rm{hydrostatic}}}} + {\boldsymbol{N}_{{\rm{wave}}}} + {\boldsymbol{N}_{{\rm{current}}}} + {\boldsymbol{N}_{{\rm{wind}}}} + {\boldsymbol{N}_{{\rm{mooring}}}} \end{array} \right.$ | (2) |
式中：下标分别表示作用于平台重心处的静水回复力、波浪力、流力、风力(含风轮)和系泊力及其它们的力矩。式(1)中B是平台绕x、y、z轴的瞬时角位移和角速度间的关系矩阵，由平台欧拉角表示为[10]
| $B = \left[ {\begin{array}{*{20}{c}} {\cos \beta }&{\sin \gamma }&0\\ { - \cos \beta \sin \gamma }&{\cos \gamma } &0\\ {\sin \beta }&0 &1 \end{array}} \right]$ | (3) |
关于风力机气动载荷数值计算模型，分别基于叶素理论、广义动态尾流法[14]计算叶元诱导速度，基于B-L(Beddoes-Leishman)动态失速理论[15]修正叶元气动力系数，在此基础上计算风轮气动载荷。广义动态尾流法模拟偏航和动态入流效应精度高，而B-L动态失速理论可修正动态入流引起的动态失速，以进一步提高风力机气动载荷的计算精度，是较适合FOWT复杂气动载荷的数值模型。考虑到Spar平台的结构特征，作用于平台的水动力载荷由细长体理论计算，系泊线载荷采用准静态悬链线法计算。此外，考虑了风力机的变桨控制，其控制模式与Hywind OC3相同。
在方程组(1)中，各物理量是时间的函数，各方程的参数之间相互关联，而且待求物理量存在乘积项，因此该方程组为非线性的耦合运动方程，需要在时域中迭代求解。为方便数值计算，将方程(1)转化为
| $\left\{ \begin{array}{l} \frac{{{\rm{d}}{\boldsymbol{U}_0}}}{{{\rm{d}}t}} = {\boldsymbol{M}^{ - 1}} \cdot \boldsymbol{F}\\ \frac{{{\rm{d}}\mathit{\pmb{\Omega}} }}{{{\rm{d}}t}} = {\boldsymbol{I}^{ - 1}} \cdot \left[ {\boldsymbol{N} - \mathit{\pmb{\Omega}} \times \left( {\boldsymbol{I} \cdot \mathit{\pmb{\Omega}} } \right)} \right]\\ \frac{{{\rm{d}}{\boldsymbol{X}_0}}}{{{\rm{d}}t}} = {\boldsymbol{U}_0}\\ \frac{{{\rm{d}}{\mathit{\pmb{\Theta}} _0}}}{{{\rm{d}}t}} = {\boldsymbol{B}^{ - 1}} \cdot \mathit{\Omega} \end{array} \right.$ | (4) |
方程(4)为二阶常微分方程。方程的求解采用四阶龙格-库塔法[12]，取时间步长Δt，在时域中步进迭代，直至收敛至设定精度。基于上述给出的FOWTs方程组、流体动力数值模型及其计算方法，编制了Fortran时域耦合程序。数值计算流程如图 1所示。
|
| 图1 数值计算流程 Figure 1 Flowchart of the numerical model |
美国NREL设计的OC3-Hywind[13]浮式风力机系统的支撑平台属典型的Spar型。将该FOWTs作为计算模型，验证本文时域耦合程序的正确性。取工况1(环境参数见表 1)运动响应计算结果，与文献[14]所给数据进行对比分析，如图 2所示。
| 风浪环境参数 | 数值 |
| 定常风风速/(m·s-1) | 8.0 |
| 规则波波高/m | 6.0 |
| 规则波周期/s | 10.0 |
|
| 图2 工况1 FOWTs响应时历曲线 Figure 2 Response history of the FOWTs in case I |
文献[14]中用HAWC2、Accoina SESAM等5种不同计算工具进行模拟分析，不同计算工具得出的计算结果有一定的差异。这里取文献中计算结果峰值最大值“上限”、谷值最小值“下限”以及一组有代表性的曲线“文献值”与本文计算结果“计算值”对比。
从图 2可以看出，本文耦合模型计算出的纵荡、升沉、艏摇运动响应都在文献计算值的上下限范围内，且运动响应曲线趋势与文献中有代表性的运动响应曲线变化趋势一致；在风轮气动响应方面，本文风轮转速计算值峰值比文献风轮转速峰值上限大0.5 r·m-1，使风轮功率峰值也超出文献峰值。总体来说，虽然本文数值模型计算的结果与文献值略有差异，但得出的风力机系统各项性能曲线与文献中性能曲线走势相一致，在FOWTs耦合运动分析和选型设计上具有应用和参考价值。
3 定常风对FOWTs运动响应的影响在工作工况下，风力机的气动推力会使FOWTs产生纵荡和纵摇运动，风力机的扭矩会使系统产生横摇运动；同时，在横摇和纵摇运动的作用下，上部风力机会偏航和俯仰，产生偏航力矩和俯仰力矩，从而引起系统艏摇，因此，定常风会对FOWTs带来多方面影响。选取3种工况(表 2)对定常风作用下FOWTs的水/气动力响应进行计算和分析。
| 工况编号 | 定常风风速/(m·s-1) | 波浪 |
| 1(较低风速) | 8.0 | 无浪 |
| 2(额定风速) | 11.2 | 无浪 |
| 3(较高风速) | 16.0 | 无浪 |
除艏摇的5自由度运动响应平均值及峰峰值列于图 3，风轮推力、扭矩和功率的平均值和峰峰值列于图 4。在该FOWTs设计时，为避免艏摇运动过大影响风力机性能，加入了艏摇弹簧抑制艏摇运动，致系统艏摇运动很小(约10-3度数量级)，因此，不对该系统的艏摇运动作详细分析。
|
| 图3 工况1~3 FOWTs运动响应均值和峰峰值对比 Figure 3 Comparison of motion in case 1~3 |
|
| 图4 工况1~3风轮气动载荷随风速变化 Figure 4 Aerodynamic loading of the turbine rotor in case 1~3 |
对于峰峰值，从图 3(b)、(d)、(f)看出，在8 m/s(较低)和16 m/s(较高)风速下，系统的纵向(纵荡、升沉和纵摇)运动均可平衡到某个运行点，运动响应的峰峰值几乎为零；而在额定风速11.2 m/s下，系统的纵向运动呈周期性变化，峰峰值分别为6.77 m、2.06 m、14°。此外，11.2 m/s风速下系统横向(横荡和横摇)运动响应峰峰值也明显比8 m/s及16 m/s风速下的峰峰值大很多(图 3(h)、(j))。这是因为在额定风速附近风力机推力变化较为剧烈(图 4(b))，从风力机的稳态推力曲线来说，在低于额定风速时推力随风速的增加而增大，而在高于额定风速时推力随风速的增加而减小，因此，在额定风速附近时，风力机系统的整体响应很难平衡在某一点，而是呈现出随时间周期性变化的特点。
对于平均值，风力机系统在11.2 m/s风速下纵向运动平均值较大(图 3(a)、(c)、(e))。这主要归因于风力机的气动推力(图 4(a))，从风力机稳态推力曲线也可进行解释，风力机在额定风速11.2 m/s附近的推力明显大于8 m/s和16 m/s风速附近的推力，推力大导致系统的纵向运动响应值变大。另外，横荡和横摇的平均值随着风速的增加而不断增大(图 3(g)、(i))，这归因于风力机的气动扭矩(图 4(c))，且从稳态扭矩曲线来讲，低于额定风速时扭矩随风速的增加而增大，高于额定风速时变化较小，因此，风力机的扭矩平均值在8 m/s附近时最小、16 m/s附近时最大，从而导致在16 m/s风速时系统的横荡和横摇运动均值最大。
4 规则波对FOWTs运动响应的影响从FOWTs风轮的推力、扭矩以及功率稳态曲线可知，在低风速、额定风速和高风速等不同风速区域，风力机气动性能有明显的差别。着重分析同一规则波中不同风速下风力机气动性能以及系统运动性能的影响，三种计算工况环境参数见表 3。
| 工况编号 | 定常风 | 规则波 | |
| 风速/(m·s-1) | 波高/m | 周期/s | |
| 4 | 8.0 | 3.0 | 8.0 |
| 5 | 11.2 | 3.0 | 8.0 |
| 6 | 16.0 | 3.0 | 8.0 |
工况4~6与工况1~3下风力机系统的5自由度运动响应结果见图 5，风力机气动载荷结果见图 6。从图 6看出，规则波对FOWT平台的纵荡、纵摇、横荡、横摇运动响应平均值的影响较小，这是因为这4个自由度的运动响应主要归因于风轮的气动推力和扭矩(图 6(a)、(b))。
从图 5、图 6看出，在定常风和规则波作用下，浮式风力机系统的各自由度运动响应、叶轮气动载荷响应峰峰值较无浪工况都变大。
|
| 图5 工况4~6FOWTs运动响应均值和峰峰值对比 Figure 5 Comparison of FOWT's motion in case 4~6 |
|
| 图6 工况4~6风轮气动载荷随风速变化 Figure 6 Aerodynamic loading of the turbine rotor in case 4~6 |
图 7和图 8分别给出了FOWTs风轮偏航角、俯仰角和偏航力矩、俯仰力矩在各风速下的平均值和峰峰值计算结果。
|
| 图7 工况4~6FOWTs风轮偏航和俯仰角均值与峰峰值 Figure 7 Yaw and tilt angle of the FOWTs in case 4~6 |
|
| 图8 工况4~6FOWTs风轮偏航和俯仰力矩均值与峰峰值 Figure 8 Yaw and tilt moment of the FOWTs in case 4~6 |
从图 7(a)、(b)看出，在8 m/s和16 m/s风速下，风轮的偏航角平均值几乎为零，偏航角的峰峰值在0.1°~0.2°范围内，可认为风力机风轮不偏航；而在11.2 m/s风速下，偏航角平均值为-0.02°、峰峰值为1°左右，虽然比其他两个风速下的平均值和峰峰值大很多，但该角度仍很小，也认为不偏航。
从图 7(c)、(d)看出，在8、16 m/s风速下，风轮的俯仰角平均值基本在3°左右，而11.2 m/s风速下为4.5°，约为前者的1.5倍；11.2 m/s风速下风轮俯仰角的峰峰值高达15°，而其它两个风速下均低于2°。
综上可见，该FOWTs风轮的俯仰明显比偏航严重很多，这从图 8所示的力矩进一步说明。在各风速下风力机偏航力矩平均值的绝对值均小于4 kN·m，峰峰值不超过50 kN·m；而俯仰力矩平均值最大可达180 kN·m，峰峰值高达500 kN·m。
风力机俯仰会明显影响风轮的气动性能，同时俯仰角峰峰值较大，还会明显影响风力机的疲劳性能，因此，应该针对风轮的俯仰性能，对FOWTs进行仔细的优化设计。
5 湍流风不规则波下FOWTs运动响应特性分析从3、4节可知，在低于和高于额定风速区域，定常风和规则波对Spar型FOWTs运动响应特性的影响规律基本一致。因此，对湍流风、不规则波下FOWTs运动响应特性的分析，只取额定风速11.2 m/s(工况7)和高风速区域的代表性风速16.0 m/s(工况8)，参数见表 4。
| 工况 编号 |
湍流风 (改进Von Karman谱) |
不规则波 (ITTC谱) |
||
| 平均风速(m·s-1) | 波高/m | 周期/s | ||
| 7 | 11.2 | 3.0 | 10.0 | |
| 8 | 16.0 | 3.0 | 10.0 | |
图 9为工况7、8的归一化后的湍流风谱和波浪谱。可以看出，本文湍流风数值模型模拟出的风谱能与目标风谱基本吻合(图 9(a)、(b))，不规则波浪数值模型模拟出的波浪谱也很好地与目标波浪谱吻合(图 9(c))，表明了本文湍流风以及不规则波数值模型的正确性。
|
| 图9 工况7、8环境风浪谱 Figure 9 Wind and wave spectra for case 7、8 |
图 10为工况7、工况8下Spar式FOWTs纵荡、纵摇及升沉运动归一化后的响应谱。从图 10(a)、(b)看出，在额定及高风速下FOWTs纵荡运动均呈现出明显的低频特性，这说明湍流风的低频特性对FOWTs的纵荡运动影响明显。从图 10(c)、(d)看出，在额定和高风速区域FOWTs的纵摇运动响应谱特性有明显不同，在额定风速处，系统的纵摇运动呈现出明显的低频性，而在高风速区域，系统的纵摇运动呈现出明显的波频特性。从图 10(e)、(f)看出，与额定风速处相比，高风速区域FOWTs升沉运动响应波频特性明显增强。这主要是因为在额定风速附近，风轮的气动推力变化较为剧烈，对FOWTs运动响应影响较大，而湍流风呈现出显著的低频特性，故致使额定风速处的运动响应低频特性要比高风速区域的运动响应低频特性明显。
|
| 图10 工况7、8FOWTs纵荡(摇)、升沉运动响应谱 Figure 10 Surge, pitch and heave spectra for case 7、8 |
1)考虑浮式风力机系统在风浪流联合作用下的瞬时受力与运动，建立了FOWTs的时域非线性耦合动力学方程组、各部件的流体动力数学模型、以及在时域中迭代求解耦合方程的方法，编写了Fortran程序源代码。针对OC3-Hywind系统，与已有结果对比验证了数学模型和程序代码的有效性，能够应用于典型Spar式FOWTs耦合性能分析。
2)定常风影响FOWTs运动响应的程度与风速范围有关。在额定风速处，定常风对5自由度运动响应的平均值和峰峰值都有显著影响。而在低于和高于额定风速的范围，仅对5自由度运动响应的均值和横荡(摇)运动的峰峰值产生一定影响。
3)规则波对FOWTs5自由度运动响应的平均值影响很小，但峰峰值明显增大。风浪环境对风轮推力、扭矩和功率的平均值和峰峰值的影响特性分别与其对纵荡(摇)、升沉的影响特性一致；从稳态推力曲线的不同区域对比，风浪环境在额定风速处对FOWTs运动响应和风轮气动性能影响最大；此外，在定常风、规则波作用下，风轮俯仰运动剧烈，影响风力机气动及疲劳性能，需针对俯仰性能对FOWTs进行优化设计。
4)随机风浪环境对FOWTs的纵荡运动响应呈现出明显的低频性；系统的纵摇运动响应在额定风速处低频性占主导地位，而在高于额定风速区波频性占主导地位；在额定及高风速区，系统升沉运动响应的低频性和波频性都较为明显，在高风速区的波频性比额定风速处明显增强。
| [1] | 徐涛. 2014年世界风能产业概况[J]. 风能产业, 2015(10): 22–28. |
| [2] | SAHIN A D. A review of research and development of wind energy in Turkey[J]. Clean-soil, air, water, 2008, 36(9): 734–742. DOI:10.1002/clen.v36:9 |
| [3] | TWIDELL J, GAUDIOSI G. Offshore wind power[M]. (s.l.): Multi-science Publishing Co. Ltd, 2009: 2 -14. |
| [4] | 张亮, 吴海涛. 海上浮式风力机技术现状及难点[J]. 风能产业, 2013(7): 22–28. |
| [5] | WANG Yi, DUAN Menglan, SHANG Jinghong. Application of an abandoned jacket for an offshore structure base of wind turbine in Bohai heavy ice conditions[C]//Proceedings of the 19th International Offshore and Polar Engineering Conference. Osaka, Japan, 2009:384-389. |
| [6] | WITHEE J E. Fully coupled dynamic analysis of a floating wind turbine system[D]. Cambridge, MA, USA:Massachusetts Institute of Technology, 2004. |
| [7] | KARIMIRAD M. Dynamic response of floating wind turbines[J]. Scientia Iranica, 2010, 17(2): 146–156. |
| [8] |
张亮, 赵玉娜, 马勇, 等. Spar式风力机平台设计及水动力影响因素研究[J].
哈尔滨工程大学学报, 2015, 36(1): 19–23.
ZHANG Liang, ZHAO Yu'na, MA Yong, et al. Design of a Spar-type wind turbine platform and analysis of the hydrodynamic influencing factors[J]. Journal of Harbin engineering university, 2015, 36(1): 19–23. |
| [9] | MA Q W, PATEL M H. On the non-linear forces acting on a floating spar platform in ocean waves[J]. Applied ocean research, 2001, 23(1): 29–40. DOI:10.1016/S0141-1187(00)00025-0 |
| [10] | MA Q W. Numerical simulation of nonlinear interaction between structures and steep waves[D]. London:University of London, 1998. |
| [11] | CHITRAPU A S, SAHA S, SALPEKAR V Y. Time-domain simulation of spar platform response in random waves and current[C]//Proceedings of the 17th International Conference on Offshore Mechanics and Arctic Engineering. Lisbon Portugal, 1998. |
| [12] | SHI Wei, PARK H C, CHUNG C W, et al. Time domain and frequency domain characterization of floating offshore wind turbine[C]//The Proceedings of the 9th ISOPE Pacific/Asia Offshore Mechanics Symposium, PACOMS-2010. Busan, Korea, 2010:91-98. |
| [13] | JONKMAN J. Definition of the floating system for phase IV of OC3[R]. Technical Report NREL/TP-500-47535. National Renewable Energy Laboratory, 2010. |
| [14] | JONKMAN J, MUSIAL W. Offshore code comparison collaboration (OC3) for IEA task 23 offshore wind technology and deployment[R]. (s.l.):Technical Report NREL/TP-500-48191. National Renewable Energy Laboratory, 2010. |
| [15] | HANSEN A C, BUTTERFIELD C P. Aerodynamics of horizontal-axis wind turbines[J]. Annual review of fluid mechanics, 1993, 25: 115–149. DOI:10.1146/annurev.fl.25.010193.000555 |
| [16] | NIELSEN F G, HANSON T D, SKAARE B. Integrated dynamic analysis of floating offshore wind turbines[C]//Proceedings of 25th International Conference on Offshore Mechanics and Arctic Engineering. Hamburg, Germany, 2006. http://cn.bing.com/academic/profile?id=2044212094&encoded=0&v=paper_preview&mkt=zh-cn |
| [17] | LARSEN T J, HANSON T D. A method to avoid negative damped low frequent tower vibrations for a floating, pitch controlled wind turbine[J]. Journal of physics:conference series, 2007, 75(1): 012073. |
| [18] | SKAARE B, HANSON T D, NIELSON F G. Importance of control strategies on fatigue life of floating wind turbines[C]//Proceedings of the 26th International Conference on Offshore Mechanics and Arctic Engineering. San Diego, CA, United States, 2007. |