基于Jourdain原理和有限元离散的浮式风机动力响应分析

刘利琴，肖昌水，郭颖，李焱

(1.天津大学 水利工程仿真与安全国家重点实验室，天津 300072； 2.天津大学 天津市港口与海洋工程重点实验室，天津 300072)

近年来，化石能源带来的环境问题日趋严重，风能作为可再生清洁能源受到人们的关注。陆上风机技术发展较为成熟，海上风机也因其优势得到迅猛发展，风机逐渐由陆地向海洋发展。出于经济效益考虑，水深超过50 m后宜采用浮式基础，包括半潜型、张力腿型(tension leg platform，TLP)、单柱式(SPAR)等。

针对浮式风机已有大量研究。Koo等[1]针对3种基础型式的风机进行模型试验。Shin等[2]针对OC3-Hywind浮式风机进行了模型试验并分析了风机系统运动特性。Jonkman等[3]开发了气动-水动-伺服-弹性全耦合浮式风机仿真程序，并分析了5 MW风机动力特性，Robertson等[4]运用该仿真程序对几种不同基础型式的浮式风机进行环境载荷分析。Ren等[5]基于计算流体动力学(computational fluid dynamics, CFD)方法计算并分析了风浪联合作用下时域运动响应。Li等[6]考虑下沉运动研究了二阶波浪力和气动载荷对TLP型浮式风机动力响应的影响。

目前大部分研究将浮式风机系统考虑为单刚体。本文将浮式风机系统处理为刚-柔耦合的多体系统，其中基础、轮毂为刚体，塔柱和叶片为柔性体，并考虑了叶片轴向变形对弯曲变形(挥舞和摆振)的影响。采用Jourdain速度变分原理，采用有限单元法离散柔性结构，推导浮式风机刚-柔耦合动力学方程，采用自动变步长的向后差分法求解微分-代数方程组，并依此编制Matlab仿真程序。建立了刚-柔耦合多体系统动力学方程，分析了浮式风力机系统的动力响应。

1 浮式风机系统动力学方程

1.1 坐标系的建立

为了分析海上浮式风机系统耦合动力响应，本文简化整机系统，将机舱考虑为塔柱末端质量模型，整机系统分为浮式基础刚体(B1)、塔柱柔体(B2)、轮毂刚体(B3)、叶片柔体(B4、B5、B6)。其中刚体部分考虑6个自由度，引入如图1所示坐标系：1)大地坐标系(O0,e(0))，惯性坐标系；2)浮式基础随体坐标系(O1,e(1))，固结于浮式基础重心；3)塔柱浮动坐标系(O2,e(2))，固结塔柱底端，描述塔柱变形；4)轮毂随体坐标系(O3,e(3))，固结于轮毂中心，随轮毂转动；5)叶片浮动坐标系(O4,e(4))、(O5,e(5))、(O6,e(6))，固结于叶片根部，描述叶片变形。

图1 浮式风机示意Fig.1 Schematic diagram of floating wind turbine

风机系统各体之间的转动位移采用卡尔丹角描述，即绕x轴转动角α、绕y轴转动角β、绕z轴转动角γ，表征转动关系的方向余弦矩阵为[7]：

(1)

式中：c表示cos；s表示sin。

角速度ω表示在随体坐标下坐标列阵，可用相对惯性基转动角的导数表示，即：

(2)

(3)

1.2 柔性梁运动学关系描述

设r0为柔性梁浮动基到惯性基的矢径，ρ为柔性梁上任意点k相对浮动基的矢径，ρ0为未变形时k点关于浮动基的矢径，u为柔性梁关于浮动基的弹性位移，有：

r=r0+ρ

(4)

ρ=ρ0+u

(5)

因此，柔性梁上k点关于惯性基的速度为:

(6)

式中：0()表示惯性基；b()表示浮动基；°表示对浮动基的局部导数。

柔性梁上k点关于惯性基的加速度为:

(7)

1.3 柔性梁变形位移描述

设ui(i=1,2,3)为k点的变形位移在浮动基的坐标分量，根据平断面假设ui(i=1,2,3)可以用中线上对应点k0的变形位移wi(i=1,2,3)表示:

(8)

式中，在小变形假设下，考虑k点的纵向变形位移的二次耦合项[8]：

对于作大范围运动的刚-柔耦合动力学系统，二次耦合项对刚度项的贡献不可不计。就风力机旋转风轮而言，不考虑二次耦合项可能导致数值仿真结果发散。

根据连续介质力学理论，k点处的纵向正应变为:

(9)

不计剪切及扭转效应的柔性梁弹性虚功率为:

(10)

1.4 空间柔性梁有限元离散

目前，变形体的离散技术大致可分为3类：假设模态法、有限元法和部件模态综合法[9]。为了反映旋转叶片的真实振型和频率，本文采用有限元法离散塔柱和叶片。

(11)

其中：

因此，中线上对应点k0的变形位移wi表示为:

(12)

每个单元的节点变形位移列阵pj可通过布尔矩阵Bj及边界条件转化矩阵R从总体变形位移列阵p中“拾取”，即：

pj=BjRp

(13)

式中：布尔矩阵Bj和边界条件转化矩阵R分别为：

本文柔性梁均为一端固支的悬臂梁，因此节点1的位移全为0，总体变形位移列阵p表示为：

p=[w12w22θ22w32θ32…w1n+1w2n+1

θ2n+1w3n+1θ3n+1]T

(14)

式中：元素下标第1位表示方向，第2位表示节点数；n为单元数。

将式(11)～(13)代入式(8)，并令Ni=Nj,iBjR，则k点相对浮动基的变形位移可表示为

(15)

式中耦合形函数为[10]：

1.5 空间梁刚-柔耦合动力学方程

本文通过Jourdain速度变分原理推导空间动力学方程，其形式为：

(16)

式中δWe为外力虚功率。

设广义坐标列阵为：

(17)

现代大型风力机都安装伺服控制系统，对于变桨调节的风力机在未达到额定风速之前，以输出功率最大的转速运行，当达到额定风速之后便以额定转速运行。因此需要在轮毂与塔柱之间施加约束条件形如下，其他连接处设为刚性固定：

(18)

式中：下标t代表塔柱；h代表轮毂；At为塔柱浮动基相对惯性基的余弦矩阵；ρt为考虑塔柱变形的塔柱浮动基到轮毂浮动基的矢径；ω为风轮转速。

因此，约束方程可表示为:

(19)

将式(6)、(7)、(11)、(15)代入动力学方程式(16)，对于刚体方程中的变形位移和弹性虚功率为零，并考虑约束方程得到形如下式的第一类拉格朗日方程：

(20)

式中：()q表示Jacobian矩阵；()t表示对t求偏导数；λ为拉格朗日乘子；M为广义质量阵；Q为广义力阵。式(20)是指标为1的微分-代数方程组，采用传统增广法求解，消去拉格朗日乘子，考虑位置约束和速度约束的违约问题，并采用向后差分法(backward differentiation formulas, BDFs)求解增广后的常微分方程组。本文数值求解为变步长积分，积分精度为10-6。

2 系统载荷计算

2.1 风机叶片非定常气动载荷

本文将经典叶素动量理论(blade element momentum, BEM)扩展到局部，计算每个叶素的非定常气动载荷，考虑风剪切、塔影效应、叶尖轮毂损失因子及Glauert修正，在此基础上引入尾流倾斜模型、B-L动态失速模型，经过修正的BEM能计算得到较为准确的气动载荷，满足本文计算要求。

由图2叶素速度三角关系可得[11]：

图2 叶素速度及作用力三角关系Fig.2 The trigonometric relationship of velocity and force of blade element

(21)

式中：Vrel,x、Vrel,y分别为来流风速叠加叶片运动速度的轴向和周向速度分量。

针对每个叶素，轴向诱导因子和切向诱导因子可表示为：

(22)

(23)

式中：B为叶片数；Cl、Cd为升阻系数；φ为入流角；F为叶尖/轮毂损失因子：

(24)

式中：R为风轮半径；Rhub为轮毂半径。

Buhl基于Glauert修正模型，提出当a≥0.4时，推力系数及轴向诱导因子的表达式[12]为：

(25)

(26)

以上推导过程基于来流风速正对风轮平面，当风轮发生偏航时，尾流将发生倾斜。对于海上浮式风机，当浮式基础发生艏摇，风轮平面发生偏航。因此，必须引入尾流倾斜模型对轴向诱导因子修正。

轴向诱导因子偏航修正可表示为：

(27)

式中：γ为偏航角；θ为风轮方位角；θ0为叶片指向尾流最深处时的方位角。

尾流倾角可用下式近似计算：

χ=(0.6a+1)γ

(28)

考虑风剪切及塔影效应，计算入流速度。风剪切模型可表示为：

(29)

式中：H为轮毂高度；vH为轮毂处风速；η为风切变指数，通常取0.1～0.25。

塔柱的存在影响了尾流分布，如图3所示，对于上风型风力机可采用潜流模型[13]：

图3 风剪切和塔影效应示意Fig.3 Schematic diagram of wind shear and tower shadow

(30)

式中：DT为叶素所在高度塔柱直径；x、y为叶素与塔柱中心的距离。式(30)在叶片处于轮毂中心以下，相对塔柱中心±60°区域内成立。当叶片处于轮毂中心之上，相对塔柱中心±60°区域内无塔影效应，其他区域修正为A(0.5-cosθ)+(0.5+cosθ)，使得各区域平滑过渡。

因此，考虑风剪切和塔影效应后的来流速度可表示为：

U∞=Av(z)

(31)

本文使用B-L动态失速模型进行动态失速的模拟，分为非定常附着流、尾缘流动分离及动态失速3个过程[14]。编制Matlab程序，模拟S809翼型动态失速过程，翼型攻角以正弦规律振动α=15+10sin(ωt)，衰减频率k=0.051 ，马赫数M=0.11，并与静态风洞试验数据对比[15],如图4所示。

图4 动态失速对升阻系数的影响Fig.4 The effect of dynamic stall on lift and drag coefficient

从图4可以发现，动态失速体现在升力失速延迟，在失速区动态气动系数与静态数据有明显差异。攻角增加变化时，升力持续增加并且要大于静态升力数据；攻角回程阶段，升力小于静态升力数据，阻力系数也有类似的现象，攻角持续变化气动系数形成迟滞环。气动仿真时，基础的运动以及叶片自身变形使得翼型攻角持续变化，因此有必要考虑动态失速的影响。

2.2 浮式基础系泊力/矩

本文以外载荷的方式考虑系泊系统对浮式风机系统动力学方面的贡献。以浮式风力机OC4-DeepCwind为例，该风机为悬链线式系泊，如图5。悬链线式系泊根据其形状，忽略水动力的影响，可通过悬链线方程快速估计系泊张力及系泊恢复力(矩)[16]。浮式基础六自由度系泊恢复力(矩)随纵荡位移变化趋势如图6所示。

图6 系泊六自由度恢复力和力矩Fig.6 Restoring force and moment of mooring system

图5 系泊系统示意Fig.5 Schematic diagram of mooring system

2.3 浮式基础水动力

在海洋工程中，小尺度构件用Morison方程计算，大尺度结构物通过三维势流理论计算。本文采用SESAM/Wadam模块进行水动力方面参数的计算，包括一阶波浪力传递函数、附连水质量、势流阻尼、静水恢复刚度。势流理论不计及粘性阻尼，工程中认为粘性阻尼可取临界阻尼β0的5%～8%，即：

(32)

式中：M0为浮式基础质量；Ma为浮式基础附连水质量；Ci为浮式基础第i自由度的静水恢复刚度。

浮式基础时域控制方程为：

(33)

3 半潜式浮式风机算例分析

本文以OC4-DeepCwind浮式风机系统为计算模型，相关参数可参考文献[17-18]。综合上述建模过程，可得到如图7所示程序框图，并编制Matlab仿真程序。

图7 仿真程序框图Fig.7 Simulation block diagram

3.1 柔性叶片固有振动特性

本文采用一次近似耦合模型，考虑二次耦合变形量，体现了叶片离心力和轴向惯性力对挥舞和摆振方向变形的影响，在动力刚度项中产生的附加耦合刚度项。考虑式(16)中δWe=0，且大范围运动已知，则方程退化为非惯性系下自由振动方程。图8为单根叶片在转速为12.1 r/min时，划分不同单元数目计算一阶挥舞固有频率的结果。

图8 验证本文方法单元收敛性Fig.8 Verification of the convergence of element

图8结果显示，叶片划分至少9个单元就能得到收敛的计算结果。采用本文程序计算零转速下叶片自由振动，叶片划分10个单元，分析其固有频率，并与FAST和ADAMS计算结果对比[17]，如表1所示，本文程序结果与二者吻合较好。

表1 本文程序结果验证Table 1 Verification of results by computer program Hz

为了揭示叶片大范围运动产生的“动力刚化”现象，将浮式基础固定并且不考虑塔柱柔性变形，采用本文程序给定转速计算叶片自由振动，固有频率随转速变化如图9所示。

图9 叶片固有频率随转速变化关系Fig.9 The frequencies of blade vary with rotor speed

由图9可以发现：一次近似耦合模型计算的挥舞和摆振方向固有频率均随着转速的增大而增加，并且转速越高，固有频率增大越快。固有频率的增大，反映出结构刚度的增大。而传统零次近似模型的摆振方向固有频率随转速的增大呈减小的趋势，挥舞方向的刚度不随转速变化。这是因为挥舞方向在旋转平面外，摆振方向在旋转平面内，传统零次近似模型忽略始终为正的附加耦合刚度项。一次近似耦合模型考虑了附加耦合刚度项，保证了旋转叶片各方向的刚度不随转速增大而减小，旋转叶片出现了“动力刚化”现象。

3.2 浮式风机系统动力响应

在额定工况下(转速12.1 r/min，风速11.4 m/s)，选定海况规则波波高6 m、周期10 s，波浪和风向均为0°，采用本文程序计算浮式风机系统动力响应。

图10为叶尖挥舞和摆振方向变形，从图中可以看出，一次近似耦合模型的变形量明显小于传统零次近似模型，结果与上述分析的旋转叶片特性相对应，考虑了叶片的“动力刚化”效应。如果采用传统零次近似模型，随着转速不断增大，叶片摆振方向的刚度不断减小，有可能导致叶片变形的仿真失败，与实际情况不符。

图10 叶尖变形时域曲线Fig.10 Blade tip deformation in time domain

通过频谱分析可以发现，叶片轴向变形二次耦合项只对叶片变形幅值有所影响，并不改变其运动规律。叶片挥舞方向振动频率以波浪和1P(风轮转速)为主，这是由于气动载荷耦合了基础的运动，受波浪影响较大。此外，挥舞变形还存在较为明显的波浪频率与1P的耦合频率、2P以及3P等频率成分。与叶片挥舞方向不同，摆振方向受到的切向气动载荷较小，因此波浪频率响应较小，而摆振方向在风轮旋转过程中受交变重力影响较大，因此摆振方向变形主要以重力引起的1P响应为主导。

图11 叶尖变形频谱分析Fig.11 Spectral analysis of blade tip deformation

由于篇幅关系，本文只列出风浪作用平面的3个基础自由度响应，即纵荡、垂荡、纵摇。基础纵荡、纵摇和垂荡的运动响应，如图12所示。

从图12中可以发现，柔性叶片的2种近似模型对浮式基础的三自由度运动影响不大，这是因为考虑叶片轴向变形二次耦合项对叶片的振动速度影响不大，如图13所示，因此2种模型计算的气动载荷相差无几。相比较而言，零次近似模型纵摇平衡位置为2.70°，而一次近似模型纵摇的平衡位置为2.71°，一次近似耦合模型的纵摇平衡位置略有增大。

图13 叶尖挥舞速度时历曲线Fig.13 Blade tip flap velocity in time domain

图12 浮式基础动力响应Fig.12 The dynamic response of floating foundation

此外，对三自由度作频谱分析，如图14所示，可以发现，浮式基础运动响应频率基本上以波频为主，气动载荷对基础运动响应频率成分贡献非常小。这是由于水平轴3根叶片在同一旋转平面内，浮式基础所受气动推力相当于叶片沿展向积分求和，叶片之间相位120°，求和之后气动载荷幅值相互抵消大大减小，与波浪载荷波动幅值相比几乎可忽略不计。

图14 浮式基础运动响应频谱分析Fig.14 Spectral analysis of response of foundation

4 结论

1) 通过计算额定转速下叶片固有频率，验证了本文建模方法的单元收敛性。针对5 MW风机叶片至少划分9个单元可得到收敛的计算结果。仅考虑叶片运动，计算零转速下叶片固有频率，与FAST和ADAMS结果对比吻合较好。

2) 仅考虑叶片运动，计算不同转速下叶片振动，分析其固有频率。传统零次近似模型计算的叶片摆振方向固有频率随转速增大而减小，挥舞方向不发生变化。一次近似耦合模型考虑旋转叶片的“动力刚化”效应，刚度随转速增大而增大，与实际情况相符。

3) 风浪联合作用下，由于一次近似耦合模型考虑了旋转叶片的“动力刚化”效应，叶片振动幅值较小，对叶片振动速度影响较小。叶片挥舞主要以波浪频率和气动1P响应为主，此外还存在二者的耦合频率以及2P、3P响应，摆振方向变形主要受交变重力影响。

4) 在定常风与规则波工况下，气动载荷对基础运动响应频率成分贡献非常小，基础运动响应以波浪频率为主。事实上在某些复杂海况下，气动载荷对基础运动的影响明显增大，因而考虑风机与浮体的耦合效应，对评估风机系统的运动性能至关重要。