NUMERICAL STUDY OF DYNAMIC CHARACTERISTICS FOR OFFSHORE WIND TURBINE UNDER COMPLEX ATMOSPHERIC INFLOW
-
摘要: 随着风能技术的不断进步, 风机叶片逐渐向大型化发展, 这使得真实复杂大气入流对风机运行性能的影响愈发显著. 为研究真实复杂大气入流下海上风机的力学特性响应, 利用基于大涡模拟的域前模拟方法生成复杂大气入流, 并结合致动线模型模拟风机叶片, 对中性复杂大气入流下海上固定式风机进行数值模拟, 重点分析风机的气动性能及转子和叶片根部的力学特性, 并与均匀入流计算工况进行对比. 计算结果表明, 中性复杂大气入流中的大尺度低速气流团使得风机气动功率输出值在较长一段时间处于较低水平, 此外, 中性复杂大气入流的高湍流强度特征使得风机气动功率的变化幅值和标准差较均匀入流工况大幅增加; 风机轴向推力的标准差值增加到均匀入流的53倍, 中性复杂大气入流的来流流场扰动引起偏航力矩的最大值、均方根和标准差分别增加到均匀入流的10、4.4和4.3倍; 速度垂向分布的不均匀性以及轮毂高度附近的大尺度低速羽流结构导致摆振剪力和弯矩的标准差响应值分别为均匀入流的2倍和4.6倍.Abstract: With the great development of wind energy technology, the blades of wind turbine have gradually developed to large-scale, which makes the real and complex atmospheric inflow have more and more significant impacts on the operating performance of wind turbines. The numerical simulation of bottom-fixed offshore wind turbine under neutral complex atmospheric inflow is performed to study the dynamic responses of wind turbine under that complex inflow. A precursor simulation method based on large eddy simulations is used to generate the complex atmospheric inflow, and the actuator line model is combined to model the wind turbine blades. The numerical results are compared with the uniform inflow condition, and the results are focusing on the analysis of aerodynamic performance and the dynamic characteristics of rotor and blade root. The numerical results show that the large-scale low-velocity airflow in the neutral and complex atmospheric inflow is responsible for the lower output of wind turbine aerodynamic power in a long period time. In addition, the high turbulence intensity characteristics of the neutral and complex atmospheric inflow lead to the significant increase of varying amplitude and standard deviation of wind turbine aerodynamic power. The standard deviation of rotor thrust increased to 53 times of the uniform inflow condition, and the maximum value, root mean square and standard deviation of yaw moment increased to 10, 4.4 and 4.3 times of uniform inflow condition, because of the disturbance of the neutral and complex atmospheric inflow. The standard deviation values of flapwise shear force and bending moment reach up to 2 and 4.6 times of uniform inflow condition, respectively, caused by the collective effects between the inhomogeneity of the velocity vertical distribution and the large-scale low-velocity plume structures near the hub height.
-
引 言
近年来, 传统化石能源的不可再生性以及对环境造成的污染, 使其难以满足当今社会对能源的巨大需求. 风能由于其无污染、可再生及储量大等特点逐渐成为最具发展潜力的新能源之一[1]. 2021年全球风能报告指出[2], 2020年新增装机容量高达93 GW, 同比增长53%. 风能技术的不断进步, 促使风机尺寸逐渐向大型化发展, 并导致大气边界层(atmospheric boundary layer, ABL)复杂来流对风机的运行性能、疲劳载荷等方面产生恶劣影响. 因此, 研究大气边界层内风机的运行性能和力学特性具有重要意义.
目前, 生成满足大气边界层复杂来流统计特性的方法主要有两种[3]: 人工合成和域前模拟. 人工合成的基本思想是在入口平面处基于严格的数学推导给出入流风速的计算公式. 宁旭[4]采用谐波合成法生成了符合高雷诺数湍流能谱分布的随机风速序列, 并将其应用于风机数值模拟. Huang等[5]提出了“DSRFG”方法, 该方法生成的脉动风速场严格满足连续性条件且能够并行计算生成不同空间位置的风速时程. Castro和 Paz[6]在“DSRFG”方法的基础上, 通过引入时间尺度参数考虑了脉动风速场的时间相关性. 但是, 人工合成法给定的风速不能严格满足NS方程, 需沿下游发展一段距离才能获得符合统计特性的大气湍流风场. Xie等[7]对传统2D涡方法进行了参数化研究, 提出了改进的2D涡方法, 将符合湍流统计特性的下游自适应距离由原来的5h减小到3h (h为计算域的高度).
域前模拟的基本做法是预先生成一个符合大气湍流统计特性的风场, 并将其作为计算区域的入口边界条件. 在域前模拟中, 生成大气湍流风场的方法主要有两种: 利用湍流激发装置以及依靠表面粗糙度. 在实验中经常利用湍流激发装置(如粗糙元、网格板等)来促使大气湍流风场的生成[8-10], 后来有部分学者将其应用于数值模拟中. Phuc 等[11]利用尖劈和粗糙元构建了数值风洞, 其数值模拟出的湍流风场统计特性与实验值符合良好. 周桐等[12]利用尖劈和粗糙元等湍流激发装置数值模拟了大气边界层入口湍流, 与“CDRFG”方法对比, 该方法生成的流场结构更加合理. 但是, 该方法模拟出的湍流风场特性高度依赖于湍流激发装置附近的网格分辨率, 且激发装置与流场相互作用机理机制尚不明确.
另一种域前模拟方法是直接依靠计算域底部的粗糙度, 无需进行复杂的网格划分, 直接长时间大范围数值模拟出大气边界层湍流风场, 如美国国家可再生能源实验室研发的风机气动性能数值计算软件“SOWFA”便基于此方法生成复杂大气入流[13]. Lee等[14]研究了不同表面粗糙度和稳定性的大气入流及风机尾流对下游风机疲劳载荷的影响, 结果表明表面粗糙度和上游风机尾流对风机疲劳载荷影响显著. Ning和Wan[15]基于LES结合致动线模型, 研究了不同大气稳定性下尾流弯曲效应及其对下游风机气动性能的影响. 白鹤鸣等[16]基于“SOWFA”软件, 对对流大气边界层下错列排布三风机进行数值模拟, 发现垂向错列可以有效提升后排风机的功率, 且不会恶化其气动性能. 李德顺等[17]利用LES结合致动线模型的方法对一台处于中性大气入流下的外场试验机组进行了叶根载荷分析, 发现叶根挥舞载荷对大气中的湍流结构响应明显. Johlas等[18-19]通过结合“SOWFA”和“OpenFAST”[20], 对大气复杂入流下浮式风机的尾流场进行了分析, 并对比了不同浮式支撑平台对风机尾流场和平台六自由度运动的影响.
以上研究中, 大都采用域前模拟中的依靠地表粗糙度的方式来生成风机的大气复杂入流, 但其中却少有文献重点关注风机的力学特性响应. 因此, 本文采用域前模拟的方法, 通过定义典型海面粗糙度, 直接长时间大范围数值模拟真实环境下中性大气边界层复杂入流. 并结合LES和致动线模型, 对中性大气复杂入流下的海上固定式风机进行数值模拟, 并以均匀入流计算结果作为参考基准, 定量分析大气复杂入流对风机气动性能和叶根载荷的影响.
1. 数值方法
1.1 致动线模型
风机叶片利用致动线模型进行模拟, 该方法无需求解叶片表面边界层, 能大幅降低计算成本并获得较为精确的结果. 其最早由Sørensen和Shen[21]提出, 基本思想是利用致动线上的致动点代表沿着径向离散化的风机叶片, 然后基于叶素理论计算致动点的体积力, 并将其反作用于流场以体现风机叶片对流场的作用. 致动点上的体积力可采用下式计算
$$ \boldsymbol{f}=(\boldsymbol{L}, \boldsymbol{D})=\frac{1}{2} \rho \boldsymbol{U}_{\rm{r e l}}^{2} c {\rm{d}} r\left(C_{L} {{\boldsymbol{e}}_L}+C_{D} {{\boldsymbol{e}}_D}\right) $$ (1) 式中,
$ \boldsymbol{L} $ 和$ \boldsymbol{D} $ 分别代表叶片半径$ r $ 处叶素的升力和阻力,$ \rho $ 为空气密度,${\boldsymbol{U}}_{\rm{rel}}$ 为相对入流风速,$ c $ 为叶片半径$ r $ 处二维翼型的弦长,${\rm{ d}}r$ 为叶素的宽度,$ {C}_{L} $ 和$ {C}_{D} $ 分别为升力系数和阻力系数, 可根据攻角确定,${{\boldsymbol{e}}_L}$ 和${{\boldsymbol{e}}_D}$ 为升力和阻力方向的单位矢量.图1显示了叶片半径处二维翼型的速度矢量, 其中相对入流风速可由下式确定
$$ \begin{array}{c}{\boldsymbol{U}}_{\rm{rel}}=\sqrt{{\boldsymbol{U}}_{\boldsymbol{z}}^{2} + {\left({{\boldsymbol{\varOmega}} }r-{\boldsymbol{U}}_{\boldsymbol{\theta }}\right)}^{2}}\end{array} $$ (2) 式中,
$ {\boldsymbol{U}}_{\boldsymbol{z}} $ 和$ {\boldsymbol{U}}_{\boldsymbol{\theta }} $ 分别为入流风速的轴向速度和切向速度,$ \boldsymbol{\varOmega } $ 为转子角速度.另外, 若直接将致动点上的体积力作用于流场, 会造成数值奇异性, 在此采用高斯核函数进行体积力光顺投影, 光顺后的体积力表达式如下
$$ \begin{split} {\boldsymbol{f}}_{\varepsilon }=&f\otimes {\boldsymbol{\eta }}_{\varepsilon }=\\ &\sum _{i=1}^{N}{\boldsymbol{f}}_{i}\left({\boldsymbol{x}}_{i},{\boldsymbol{y}}_{i},{\boldsymbol{z}}_{i},t\right)\frac{1}{{\varepsilon }^{3}{{\text{π}} }^{\frac{3}{2}}}{\rm{exp}}\left[-{\left(\frac{{\boldsymbol{d}}_{i}}{\varepsilon }\right)}^{2}\right]\end{split} $$ (3) 式中,
$ N $ 为叶片上致动点的数目,$ \left({\boldsymbol{x}}_{i},{\boldsymbol{y}}_{i},{\boldsymbol{z}}_{i}\right) $ 表示第$ i $ 个致动点的位置,$ {\boldsymbol{d}}_{i} $ 为致动点距离投影点的距离,$ \varepsilon $ 代表投影的宽度, 本文取$ \varepsilon \approx 2\mathrm{\Delta }x $ [22],$ \mathrm{\Delta }x $ 为叶片附近的网格尺寸.1.2 控制方程
本文采用大涡模拟方法, 对中性复杂大气入流下的海上固定式风机进行数值模拟, 空间过滤后的流场的控制方程如下
$$ \dfrac{\partial {\stackrel-{\boldsymbol{u}}}_{i}}{\partial {\boldsymbol{x}}_{i}}={\boldsymbol{0}} $$ (4) $$ \begin{split} & \frac{\partial \overline{\boldsymbol{u}}_{i}}{\partial t}+\frac{\partial}{\partial \boldsymbol{x}_{j}}\left(\overline{\boldsymbol{u}}_{j} \overline{\boldsymbol{u}}_{i}\right)=-\underbrace{\frac{\partial \hat{p}}{\partial \boldsymbol{x}_{i}}}_{{\rm{I}}} -\underbrace{\frac{1}{\rho_{0}} \frac{\partial}{\partial \boldsymbol{x}_{i}} \bar{p}_{0}(\boldsymbol{x}, \boldsymbol{y})}_{{\rm{I I}}} -\\ &\qquad \underbrace{2 \boldsymbol{\varepsilon}_{i 3 k} \boldsymbol{\varOmega}_{3} \overline{\boldsymbol{u}}_{k}}_{{\rm{I I I}}}+\underbrace{g\left(\frac{\bar{\theta}-\theta_{0}}{\theta_{0}}\right) \boldsymbol{\delta}_{i 3}}_{{\rm{I V}}}-\underbrace{\frac{\partial}{\partial \boldsymbol{x}_{j}}\boldsymbol{\tau}_{i j}^{D}}_{{\rm{V}}}+\underbrace{\frac{1}{\rho_{0}} \boldsymbol{f}_{i}^{{\rm{T}}}}_{{\rm{VI}} }\end{split} $$ (5) 式中, 动量方程右端第一项为修正压力梯度项, 第二项为背景压力梯度项, 用以克服计算域底部表面粗糙度并驱动大气湍流的生成, 第三项为考虑地球自转的科氏力项, 第四项代表温度引起的浮力项, 第五项为流体应力张量项, 采用“Smagorinsky”亚格子模型封闭, 第六项为风机叶片体积力源项, 由致动线模型计算得出, 在域前模拟阶段忽略此项. 另外, 还需要求解位温输运方程以获得位温场
$$ \begin{array}{c}\dfrac{\partial \stackrel-{\theta }}{\partial t} + \dfrac{\partial }{\partial {\boldsymbol{x}}_{j}}\left({\boldsymbol{u}}_{j}\stackrel-{\theta }\right)=-\dfrac{\partial }{\partial {\boldsymbol{x}}_{j}}{\boldsymbol{q}}_{j}\end{array} $$ (6) 式中, 热通量
${{\boldsymbol{q}}}_{j}$ 由下式计算$$ \begin{array}{c}{\boldsymbol{q}}_{j}=-\dfrac{{\nu }_{{\rm{SGS}}}}{{Pr}_{t}}\dfrac{\partial \stackrel-{\theta }}{\partial {\boldsymbol{x}}_{j}}\end{array} $$ (7) 式中,
$ {Pr}_{t} $ 为普朗特数, 在中性和对流大气边界层中取$ 1/3 $ , 关于控制方程更详细的说明可查阅参考文献[23]. 对于均匀入流下的风机, 其流场的控制方程未考虑科氏力项和温度项的影响, 也不需要背景压力梯度驱动流体流动, 在此不再重复赘述.1.3 局部坐标系
在本文的结果分析中, 需要对风机转子气动力和力矩以及叶片根部的力和力矩进行结果分析, 两者所采用的坐标系不同, 在此进行简要说明.
风机转子气动力和力矩采用的是轮毂局部坐标系, 如图2(a)所示, 其坐标原点为转子轴与旋转平面的交点, x轴为沿着轮毂中心线指向下游方向, z轴垂直于轮毂中心线并与叶片具有相同的方位角, y轴根据右手螺旋定则确定.
叶片根部的力和力矩基于叶片局部坐标系, 如图2(b)所示, 其中i代表叶片编号. y轴平行于弦线并指向叶片后缘, z轴为沿着桨距轴指向叶尖, x轴根据右手螺旋定则确定. 另外, 两个局部坐标系均会随着叶片旋转.
2. 算例设置
2.1 风力机参数
本文所采用的研究对象为美国国家可再生能源实验室(National Renewable Energy Laboratory, NREL)开发的5 MW风机, 其是一种传统迎风向可变桨距和转速的3叶片风力机. 该风机的额定风速为
$ 11.4\; {\rm{m/s }}$ , 对应的额定功率和额定转速分别为5 MW和12.1 r/min, 该风机的主要参数如下表1所示, 关于该风机更详细的数据可参考文献[24].表 1 NREL 5 MW风机主要参数Table 1. Main parameters of NREL 5 MW wind turbineParameter Value rated power/MW 5 rated wind speed/(m·s−1) 11.4 rated rotor speed/(r·min−1) 12.1 hub height/m 90 number of blades 3 rotor orientation upwind 2.2 大气入流算例设置
考虑真实复杂大气入流的海上固定式风机数值模拟需要分两阶段进行: 域前模拟和主模拟. 在域前模拟阶段, 需要长时间大范围的依靠底面粗糙度来生成真实的大气湍流风场. 如图3所示, 计算域的四周侧壁均采用周期性边界条件, 以达到通过减小计算域尺寸降低计算量的目的, 计算域顶部设置为滑移边界条件, 底面采用“Moeng”壁面模型[25], 表面粗糙度设置为0.001, 对应典型低粗糙度的海面情况. 计算域的长度、宽度和高度分别设置为
$ 3000 \;{\rm{m}}\times $ $ 3000\;{\rm{ m}}\times 1000\; {\rm{m}} $ , 三个方向的网格分辨率为$ 10\; {\rm{m}}\times $ $ 10\; {\rm{m}}\times 10\; {\rm{m}} $ , 对应的网格数量为900万. 轮毂高度处的风速设置为风机的额定风速$ 11.4\; {\rm{m/s}} $ , 入流风向设置为240°(西南方向入流)以避免在周期性边界条件中同一高度处的涡结构几乎同时到达下游边界时又同时传入上游边界时“卡死”现象的出现.![]()
图 3 域前模拟计算域及边界条件(单位: km)Figure 3. Calculation region and boundary conditions of precursor simulation (unit: km)对于初始位温的设置, 从底部至700 m高度设置300 K均匀分布, 700 m至800 m高度设置逆温反转层, 其中的温度线性增长到308 K, 800 m至1000 m高度的温度以
$ 0.003\;{\rm{ K/m}} $ 的速率线性增长. 为了使生成的大气湍流达到准平衡态, 域前模拟阶段的计算时间为19000 s, 计算步长为0.4 s, 保留最后1000 s的计算数据作为主模拟阶段的输入.主模拟阶段的计算域和背景网格的分辨率与域前模拟阶段相同, 但北面和东面侧壁边界修改为零梯度边界条件以防止下游风机尾流循环进入上游边界, 如图4所示. 另外, 为了更加精确的捕捉风机尾流区域的流场细节, 对风机附近的区域进行了二级网格加密, 加密后风机附近的网格尺寸为
$ \rm 2.5\; m\times $ $ \rm2.5 \;m\times 2.5\; m $ , 总网格数量约为1400万. 时间步长的设置需要满足CFL收敛性条件$$ \begin{array}{c}{\rm{max}}\left(\left|\dfrac{U\Delta t}{\Delta x}\right|,\left|\dfrac{\mathrm{\varOmega }R\Delta t}{\Delta x}\right|\right) < 1\end{array} $$ (8) 式中,
$ U=11.4\; {\rm{m/s}} $ , 为轮毂高度入流风速;$ R=63 \;{\rm{m}} $ , 为转子半径;$ \mathrm{\varOmega }=12.1\;{\rm{ r/min}} $ , 为入流风速下的转子转速;$ \mathrm{\Delta }x=2.5\;{\rm{ m}} $ , 为风机叶片附近的网格分辨率. 代入上式, 可得$ \mathrm{\Delta }t < 0.031\;{\rm{ s}} $ , 因此, 本文的时间步长设置为0.02 s. 计算时长为1000 s, 对应于域前模拟阶段的流场数据保存时间, 并利用最后350 s的计算数据进行结果分析. 本课题组在风机数值模拟方向具有扎实的研究基础, 关于风机附近网格划分策略的详细细节可参考文献[26-28].2.3 均匀入流算例设置
如图5所示, 均匀入流算例计算域的长度、宽度和高度分别设置为
$ 10 D\times 3 D\times 3 D $ , D为转子直径126 m. 风机距离上游边界为$ 2 D $ , 同样, 与主模拟做法相一致, 背景网格的尺寸为$ 10\;{\rm{ m}}\times 10 \;{\rm{ m}}\times 10 \;{\rm{ m}} $ , 并对风机叶片附近区域采取二级网格加密策略, 加密后的网格尺寸为$ 2.5 \;{\rm{ m}}\times 2.5 \;{\rm{ m}}\times 2.5 \;{\rm{ m}} $ , 网格数量约为320万. 另外, 上游边界为速度入口边界条件, 对应于11.4 m/s的均匀入流, 入流风向为270°(西侧方向入流), 下游边界为自由流出边界条件, 四周壁面为周期性边界条件. 时间步长为0.02 s, 计算时间为400 s, 并利用最后350 s的计算数据进行结果分析.3. 结果分析
3.1 气动功率
图6显示了计算结果最后350 s中性大气入流和均匀入流下海上固定式风机的功率时历曲线. 可以看出, 均匀入流下的转子气动功率输出较为稳定, 但大气复杂入流下转子功率输出变化十分复杂, 转子功率在
$ t=25 \;{\rm{s}} $ 附近开始存在一个较低值, 且持续时间超过50 s, 这是由大气复杂入流中的大尺度低速气流导致. 另外, 随着时间的增加, 气动功率值迅速增加至6 MW附近, 在$ t=90\;{\rm{s}} $ 附近又存在较大的下降梯度, 随后又逐渐增加, 并在均匀入流风机的气动功率值附近振荡.在数据统计特性方面, 分别统计了气动功率的最大值、最小值、平均值、均方根和标准差, 如表2所示. 可以看出, 两种工况的气动功率均方根值差异不大, 分别为5.22 MW和5.30 MW, 因为气动功率的均方根值主要由平均风速决定, 两种计算工况的平均风速均在11.4 m/s附近. 但是, 大气入流中的复杂湍流会引起气动功率振荡幅值的极大增加, 进而引起标准差的增加, 其值约为均匀入流气动功率标准差的4倍.
表 2 气动功率统计值Table 2. Aerodynamic power statisticsCase Aerodynamic power/MW max min mean rms std ABL 6.04 4.20 5.21 5.22 0.26 uniform 5.33 5.26 5.30 5.30 0.07 3.2 轴向推力和偏航力矩
轴向推力和偏航力矩属于风机转子的气动力和力矩, 其局部坐标系为轮毂坐标系, 其中轴向推力沿轮毂坐标系x轴, 偏航力矩绕z轴. 图7(a)和图7(b)分别显示了大气入流和均匀入流两种工况下风机转子的轴向推力和偏航力矩对比曲线. 与气动功率输出类似, 均匀入流下风机的轴向推力稳定在600 kN附近, 而复杂大气入流下的轴向推力输出无论是变化幅度还是变化梯度均较大, 这同样归因于复杂大气入流下的大尺度湍流来流.
当风机转子的其中两个叶片位于塔架同一侧时会导致盘面受力的横向不对称性, 并由此产生偏航力矩. 因此, 即便在均匀入流下, 风机偏航力矩也会随叶片旋转作周期性变化, 变化周期对应于叶片旋转周期. 然而, 复杂大气入流下的流场扰动会使得风机的偏航力矩随时间呈现明显区别于均匀入流工况的复杂变化, 明显加剧偏航力矩的周期性变化幅度, 并显著增加其均方根值, 如图7(b)所示. 因此, 在风机偏航优化控制中需重点关注复杂大气入流带来的不利影响[29].
表3显示了两种入流工况下轴向推力和偏航力矩的统计值. 虽然轴向推力的均方根值差异不大, 但复杂大气入流下的轴向推力标准差与均匀入流相比增大了约53倍. 另外, 针对偏航力矩, 复杂大气入流的流场扰动极大的增加了其幅值, 偏航力矩的最大值和最小值分别增加到均匀入流的10倍和8倍, 并导致均方根和标准差的比值为均匀入流的4.3倍, 相关文献[4]也表明复杂入流下各类结构载荷的标准差值较均匀入流时均有大幅增加.
表 3 轴向推力和偏航力矩统计值Table 3. Statistics of rotor thrust and yaw momentCase Rotor thrust/kN max min mean rms std ABL 745 559 632 633 35.2 uniform 605 598 601 601 0.67 Case Yaw moment/(kN·m) max min mean rms std ABL 2990 −2409 −111 803 796 uniform 294.4 −294.7 1.85 184 184 3.3 摆振剪力和弯矩
摆振剪力和弯矩的参考坐标系为风机叶片局部坐标系, 其中摆振剪力沿x轴, 摆振弯矩绕y轴. 图8分别显示了摆振剪力和弯矩的对比曲线由图8(a)可以看出, 均匀入流下叶根部位的摆振剪力随叶片旋转作周期性变化, 变化周期对应于叶片旋转周期. 另外, 大气入流下的摆振剪力也近似作周期性变化, 但是, 由于复杂大气入流下的高强度湍流结构及速度沿垂向分布的不均匀性, 导致摆振剪力的变化幅度较为剧烈且随机. 摆振弯矩主要由摆振剪力引起, 两者的变化趋势基本一致, 在此不再重复赘述.
表4给出了摆振剪力和弯矩的统计值, 可以看出, 复杂大气入流和均匀入流工况下摆振剪力和弯矩的均方根值差异不大. 但是, 由于复杂大气入流的流向风速在垂向上存在速度梯度, 当风机叶片旋转至轮毂平面上方时, 摆振剪力和弯矩值较均匀入流增大; 反之, 较均匀入流减小, 这导致了摆振剪力和弯矩标准差值的显著增加.
表 4 摆振剪力和弯矩统计值Table 4. Statistics of flapwise shear force and bending momentCase Flapwise shear force/kN max min mean rms std ABL 323.4 174.8 252.8 253.5 18.4 uniform 258.9 230.4 244.5 244.7 9.11 Case Flapwise bending moment/(kN·m) max min mean rms std ABL 11530 6002 8943 8974 741.7 uniform 8861 8347 8609 8610 160.1 3.4 速度云图
图9显示了复杂大气入流下某时刻的流场速度云图. 由图9(a)可以看出, 轮毂高度处的来流速度同时存在大尺度的高速和低速的气流团结构, 并且在风轮旋转平面处存在明显的不对称性, 使得偏航力矩振荡幅值急剧增加. 由于风机吸收了入流风的能量, 导致风机盘面后方的尾流区域存在明显的速度损失. 本文并未对轮毂进行建模, 因此在轮毂中心处存在一条狭长的高速气流, 该气流随着向下游传播迅速和低速风机尾流相互掺混. 另外, 尾流膨胀和蜿蜒效应也清晰可见.
图9(b)展示的是240°入流风垂直面的瞬时尾流速度云图, 可以看出, 由于海面的阻滞作用, 在靠近计算域底部高度处存在着大尺度低速气流结构. 另外, 在轮毂高度附近存在大尺度低速羽流结构, 其与计算域底部的大尺度低速气流团相互作用会显著增加叶片摆振剪力和弯矩的振荡幅值. 在风轮旋转平面下游, 可以明显看出风机尾流与外部大气之间的掺混作用. 但由于海面低粗糙度的作用, 使得中性复杂大气入流的环境湍流强度较高地表粗糙度时低, 导致风机尾流速度随下游流向距离恢复较慢.
3.5 尾涡结构
第三代涡识别方法能够合理回答涡的六大要素问题, 并定量化表示涡的特性. 因此, 本文采用第三代“Liutex”涡识别方法[30], 对中性复杂大气入流下海上固定式风机的尾涡结构进行可视化识别. 根据图10可以看出, 由于中性大气入流的复杂湍流场, 导致在风机上游就已形成多尺度湍流涡结构. 风机叶片的旋转会诱导出清晰可见的叶尖涡, 该叶尖涡与复杂大气入流涡结构的相互作用引起叶尖涡的迅速破碎, 并与小尺度环境涡结构相互掺混向下游传播. 另外, 风机尾涡的膨胀效应也清晰可见.
![]()
图 10 尾涡结构($ \left|Liutex\right|=0.18 $ )Figure 10. Vortex structures ($ \left|Liutex\right|=0.18 $ )4. 结论
本文采用基于大涡模拟的域前模拟方法, 生成了真实环境下中性大气边界层复杂入流, 并利用致动线模型对风机叶片进行数值模拟, 定量研究了中性复杂大气入流下海上固定式风机的力学特性响应, 重点分析了风机的气动性能以及转子和叶片根部的力和力矩, 通过与均匀入流工况进行对比, 得出了以下结论.
(1)气动功率的均方根值由平均入流风速决定, 但是由于中性大气来流的复杂性, 导致气动功率输出振荡幅值加剧. 另外, 大气入流中的大尺度低速气流团使得气动功率输出值在较长一段时间内过低.
(2)中性大气复杂入流使得风机轴向推力的标准差相较与均匀入流急剧增加, 风机盘面前的复杂来流扰动引起偏航力矩最值、均方根值以及标准差值的较大增加.
(3)中性复杂大气入流流向速度沿垂向的梯度、流场的复杂性以及轮毂高度处大尺度低速羽流结构的共同作用, 加剧了叶片根部摆振剪力和弯矩的响应幅值和标准差值.
随着风机逐渐向大型化发展, 复杂大气入流会进一步恶化风机叶片的入流条件, 增加风机叶片的力学响应幅值, 并引起标准差值的极大增加. 因此, 复杂大气入流下风机叶片的力学响应特性需得到更加广泛关注.
-
![]()
图 3 域前模拟计算域及边界条件(单位: km)
Figure 3. Calculation region and boundary conditions of precursor simulation (unit: km)
![]()
图 10 尾涡结构(
$ \left|Liutex\right|=0.18 $ )Figure 10. Vortex structures (
$ \left|Liutex\right|=0.18 $ )表 1 NREL 5 MW风机主要参数
Table 1 Main parameters of NREL 5 MW wind turbine
Parameter Value rated power/MW 5 rated wind speed/(m·s−1) 11.4 rated rotor speed/(r·min−1) 12.1 hub height/m 90 number of blades 3 rotor orientation upwind 
下载: 导出CSV
表 2 气动功率统计值
Table 2 Aerodynamic power statistics
Case Aerodynamic power/MW max min mean rms std ABL 6.04 4.20 5.21 5.22 0.26 uniform 5.33 5.26 5.30 5.30 0.07
下载: 导出CSV
表 3 轴向推力和偏航力矩统计值
Table 3 Statistics of rotor thrust and yaw moment
Case Rotor thrust/kN max min mean rms std ABL 745 559 632 633 35.2 uniform 605 598 601 601 0.67 Case Yaw moment/(kN·m) max min mean rms std ABL 2990 −2409 −111 803 796 uniform 294.4 −294.7 1.85 184 184
下载: 导出CSV
表 4 摆振剪力和弯矩统计值
Table 4 Statistics of flapwise shear force and bending moment
Case Flapwise shear force/kN max min mean rms std ABL 323.4 174.8 252.8 253.5 18.4 uniform 258.9 230.4 244.5 244.7 9.11 Case Flapwise bending moment/(kN·m) max min mean rms std ABL 11530 6002 8943 8974 741.7 uniform 8861 8347 8609 8610 160.1
下载: 导出CSV
-
[1] Chehouri A, Younes R, Ilinca A, et al. Review of performance optimization techniques applied to wind turbines. Applied Energy, 2015, 142: 361-388 doi: 10.1016/j.apenergy.2014.12.043
[2] Council GWE. GWEC global wind report 2021//Global Wind Energy Council, Brussels, Belgium, 2021
[3] 周桐, 闫渤文, 杨庆山等. 大气边界层大涡模拟入口湍流生成方法研究. 工程力学, 2020, 37(7): 68-76 (Zhou Tong, Yan Bowen, Yang Qingshan, et al. Study of inflow turbulence generation methods with large eddy simulation for atmospheric. Engineering Mechanics, 2020, 37(7): 68-76 (in Chinese) doi: 10.6052/j.issn.1000-4750.2019.07.0381 [4] 宁旭. 大气边界层内风机气动及尾流特性的数值研究. [硕士论文]. 上海: 上海交通大学, 2020 Ning Xu. Numerical study of wind turbine aerodynamics and wake characteristics under atmospheric boundary layer. [Master Thesis]. Shanghai: Shanghai Jiao Tong University, 2020 (in Chinese)
[5] Huang SH, Li QS, Wu JR. A general inflow turbulence generator for large eddy simulation. Journal of Wind Engineering and Industrial Aerodynamics, 2010, 98(10-11): 600-617 doi: 10.1016/j.jweia.2010.06.002
[6] Castro HG, Paz RR. A time and space correlated turbulence synthesis method for large eddy simulations. Journal of Computational Physics, 2013, 235: 742-763 doi: 10.1016/j.jcp.2012.10.035
[7] Xie B, Gao F, Boudet J, et al. Improved vortex method for large-eddy simulation inflow generation. Computers & Fluids, 2018, 168: 87-100
[8] Hlevca D, Degeratu M. Atmospheric boundary layer modeling in a short wind tunnel. European Journal of Mechanics-B/Fluids, 2020, 79: 367-375 doi: 10.1016/j.euromechflu.2019.10.003
[9] Li QA, Murata J, Endo M, et al. Experimental and numerical investigation of the effect of turbulent inflow on a horizontal axis wind turbine (Part I: Power performance). Energy, 2016, 113: 713-722 doi: 10.1016/j.energy.2016.06.138
[10] Murata J, Endo M, Maeda T, et al. Experimental and numerical investigation of the effect of turbulent inflow on a horizontal axis wind turbine (Part II: Wake characteristics). Energy, 2016, 113: 1304-1315 doi: 10.1016/j.energy.2016.08.018
[11] Phuc PV, Nozu T, Kikuchi H, et al. Wind pressure distributions on buildings using the coherent structure Smagorinsky model for LES. Computation, 2018, 6(2): 32 doi: 10.3390/computation6020032
[12] 周桐, 杨庆山, 闫渤文等. 大气边界层大涡模拟入口湍流生成方法综述. 工程力学, 2020, 37(5): 15-25 (Zhou Tong, Yang Qingshan, Yan Bowen, et al. Review of inflow turbulence generation methods with large eddy simulation for atmospheric boundary layer. Engineering Mechanics, 2020, 37(5): 15-25 (in Chinese) [13] Fleming P, Gebraad P, van Wingerden JW, et al. SOWFA super-controller: A high-fidelity tool for evaluating wind plant control approaches//National Renewable Energy Laboratory, Golden, CO, USA, 2013
[14] Lee S, Churchfield M, Moriarty P, et al. Atmospheric and wake turbulence impacts on wind turbine fatigue loadings//Proceedings of the 50th AIAA Aerospace Sciences Meeting including the New Horizons Forum and Aerospace Exposition, 2012: 540
[15] Ning X, Wan D. LES study of wake meandering in different atmospheric stabilities and its effects on wind turbine aerodynamics. Sustainability, 2019, 11(24): 6939 doi: 10.3390/su11246939
[16] 白鹤鸣, 万德成, 王尼娜等. 大气边界层入流下错列排布三风机气动性能数值模拟. 水动力学研究与进展(A辑), 2021, 36(1): 10-19 (Bai Heming, Wan Decheng, Wang Nina. Numerical simulation of aerodynamic performance of three wind turbines with staggered strategies under atmospheric boundary layer flow. Chinese Journal of Hydrodynamics, 2021, 36(1): 10-19 (in Chinese) [17] 李德顺, 郭涛, 李伟等. 中性大气边界层中风力机的湍流演化及叶根载荷分析. 科学通报, 2019, 64(17): 1832-1843 (Li Deshun, Guo Tao, Li Wei, et al. Evolution of turbulence in a wind turbine flow field with a neutral atmospheric boundary layer and an analysis of the blade root load. Chinese Science Bulletin, 2019, 64(17): 1832-1843 (in Chinese) doi: 10.1360/N972019-00213 [18] Johlas HM, Martínez-Tossas LA, Schmidt DP, et al. Large eddy simulations of floating offshore wind turbine wakes with coupled platform motion. Journal of Physics:Conference Series, 2019, 1256(1): 012018 doi: 10.1088/1742-6596/1256/1/012018
[19] Johlas HM, Martínez-Tossas LA, Lackner MA, et al. Large eddy simulations of offshore wind turbine wakes for two floating platform types. Journal of Physics:Conference Series, 2020, 1452(1): 012034 doi: 10.1088/1742-6596/1452/1/012034
[20] Jonkman JM, Buhl ML. FAST User's Guide. Golden, CO, USA: National Renewable Energy Laboratory, 2005
[21] Sørensen JN, Shen WZ. Numerical modeling of wind turbine wakes. Journal of Fluids Engineering, 2002, 124(2): 393-399 doi: 10.1115/1.1471361
[22] Troldborg N. Actuator line modeling of wind turbine wakes. [PhD Thesis]. Denmark: Technical University of Denmark, 2008
[23] Churchfield MJ, Lee S, Michalakes J, et al. A numerical study of the effects of atmospheric and wake turbulence on wind turbine dynamics. Journal of Turbulence, 2012, 13: N14
[24] Jonkman J, Butterfield S, Musial W, et al. Definition of a 5-MW reference wind turbine for offshore system development//National Renewable Energy Laboratory, Golden, CO, USA, 2009
[25] Moeng CH. A large-eddy-simulation model for the study of planetary boundary-layer turbulence. Journal of the Atmospheric Sciences, 1984, 41(13): 2052-2062 doi: 10.1175/1520-0469(1984)041<2052:ALESMF>2.0.CO;2
[26] Cheng P, Huang Y, Wan D. A numerical model for fully coupled aero-hydrodynamic analysis of floating offshore wind turbine. Ocean Engineering, 2019, 173: 183-196 doi: 10.1016/j.oceaneng.2018.12.021
[27] Huang Y, Cheng P, Wan D. Numerical analysis of a floating offshore wind turbine by coupled aero-hydrodynamic simulation. Journal of Marine Science and Application, 2019, 18(1): 82-92 doi: 10.1007/s11804-019-00084-8
[28] Huang Y, Wan D. Investigation of interference effects between wind turbine and spar-type floating platform under combined wind-wave excitation. Sustainability, 2020, 12(1): 246
[29] Yang J, Fang L, Song D, et al. Review of control strategy of large horizontal-axis wind turbines yaw system. Wind Energy, 2021, 24(2): 97-115 doi: 10.1002/we.2564
[30] 王义乾, 桂南. 第三代涡识别方法及其应用综述. 水动力学研究与进展(A辑), 2019, 34(4): 413-429 (Wang Yiqian, Gui Nan. A review of the third-generation vortex identification method and its applications. Chinese Journal of Hydrodynamics, 2019, 34(4): 413-429 (in Chinese)
计量
- 文章访问数: 676
- HTML全文浏览量: 313
- PDF下载量: 119
