旋转作用下水平轴风力机气动性能分析

2023-10-09熊宇轩叶祖洋

科学技术与工程 2023年26期

熊宇轩, 叶祖洋*

(1.武汉科技大学资源与环境工程学院, 武汉 430081; 2.冶金矿产资源高效利用与造块湖北省重点实验室, 武汉 430081)

随着能源消耗问题日益严峻,能源结构调整迫在眉睫,风能作为一种清洁能源,为推动能源可持续转型,在能源中比重不断增加[1]。风力机作为一种通过能量转换将风能转化为机械能的发电设备,其生产发电效益受风电场中环境风变化规律的影响。风力机大部分时间在非稳定流动环境中运行,受风剪切、旋转效应、大气湍流、倾斜和偏航失调等非定常因素影响的各种流动条件影响,导致风轮的不稳定负载[2]。风力机故障、服役寿命缩短、输出功率降低和运行维护增加都与此直接相关[3]。同时,因为叶片流场是不稳定的三维分离流,会导致风力机动态失速,并伴随着巨大的负载变化,致使空气动力学建模和非定常载荷的精确预测面临许多挑战。

风力机叶片在复杂环境下气动性能分析是研究整个风电机组发电效率的重要方法。对风力机非定常来流作用下气动特性的研究具有重要意义。李欢等[4]对风力机进行了气动力和离心载荷模拟,并得出两者对风力机结构特性的影响程度。张建平等[5]通过函数拟合不同桨距角下叶片位移和应力,给出模拟计算下叶片动力特性变化规律。代元军等[6]分析对比了叶尖结构对风力机流场计算流体力学(computational fluid dynamics,CFD)数值模拟结果和风洞实验规律,证明了两者的一致性。李德源等[7]采用升力线理论模型和叶素动量理论对风力机叶片气动性能进行分析对比,验证了计算模型的可靠性。朱呈勇等[8]建立了基于Leishman-Beddoes模型和Bak模型的动态失速和旋转效应气动响应预测模型,有效预测了耦合作用对非定常气动力的影响程度。王占洋等[9]通过考虑流固耦合作用探讨风力机应力应变分布特征,证明了耦合作用下风力机输出功率更贴近试验值。Cao等[10]研究了不同条件下叶素动量理论计算作用在叶片上的气动荷载动态响应分析,验证了理论结果和实验结果的正确性。Li等[11]采用多体动力学代码集成CFD模拟方法,研究了在大气非定常情况下尾流风湍流和转子相互作用,表明气动载荷和叶尖挠度在时域和频域上具有良好的一致性。

针对风力机非定常影响的气动性能研究,大多只考虑正常来流风条件下的影响,现考虑均匀来流和三维旋转效应共同作用下非定常条件,结合ANSYS平台流固耦合分析,通过结构场接收流场稳态计算结果,传递耦合面压力、速度等数据研究风轮固体结构在流体作用下气动耦合特性。

1 流固耦合方程

流固耦合是通过数值模拟软件求解模块采用弱耦合分离式求解方法按照设定顺序将所建立的流体域计算流体动力学(computational fluid dyna-mics,CFD)方程和固体域结构动力学(computational structural dynamics,CSD)方程计算数据进行交互传输处理,从而实现各个时间步上的流固耦合求解。流体域模拟采用ANSYS FLUENT,该模块代码使用有限体积法求解流体的控制方程。在本文研究中,对整个流动区域求解不可压缩、非定常雷诺平均Navier-Stokes(URANS)方程。选择基于耦合压力的求解器采用二阶隐式公式,以提高精度,其中解变量采用二阶迎风离散化格式求解。

1.1 流体控制方程

质量守恒方程为

(1)

动量守恒方程为

(2)

式中:t为时间;ff为体积力矢量;ρf为流体密度;v为流体速度矢量;τf为剪切力张量。

剪切力张量表达式为

(3)

1.2 固体控制方程

固体部分系统守恒方程由牛顿第二定律推导为

(4)

1.3 流固耦合基本守恒方程

计算耦合交界处存在关于应力、位移等变量的守恒,满足

(5)

式(5)中:τ为应力;d为位移;n为方向余弦;下标f表示流体;下标s表示固体。

2 计算模型

2.1 湍流模型

水平轴风力机气动耦合分析采用混合k-ω/k-ε剪切应力传输(shear stress transfer,SST)湍流模型,k-ω的SST模型可以考虑湍流剪切的传输并准确预测不同压力梯度条件下流量传输情况,k-ε模型则用于外部区域和自由剪切流。利用指标残差监测来确定计算结果的收敛性,将流量连续性、xyz方向速度、湍流动能和比耗散率的残差收敛精度设置为1×10-6、1×10-3。压力、动量、湍流动能、比耗散率和比能量的空间离散均近似于二阶。SST模型为

(6)

(7)

混合函数F1定义为

(8)

2.2 风力机几何模型

叶片截面叶素坐标使用Profili软件中NACA44XX翼型坐标,通过多线段样条曲线导入翼型坐标,在同一轴线上按翼型攻角缩放、扭转叶素曲线,形成叶片基础骨架,使用铺层表面完成叶片建模。几何结构参考1.5 MW NREL水平轴风力机模型,风轮直径为6.5 m,塔筒高度为8.75 m,塔筒顶端截面直径为0.25 m,底部截面直径为0.45 m,机舱通过简化设计为0.9 m×0.7 m×0.525 m立方体结构,如图1所示。

图1 风力机三维几何模型

2.3 网格划分和边界条件配置

为了模拟旋转状态下的风力机,将封闭整个风力机转子的移动圆柱形旋转内域与包含机舱和塔筒的固定长方体外域相结合。圆柱形内部区域的直径和厚度分别为6.8 m和0.6 m。叶片与内部旋转域通过布尔运算完成内流域的创建。模拟均匀来流风速的非定常影响并减少边界条件的影响,采取将外流域中塔筒和内流域挖除的方法,整个计算区域的尺寸为40 m×30 m×21 m,如图2所示。两个相对运动的区域共享一个滑动界面,作为数据传输耦合面。

图2 风力机计算域示意图

2.3.1 CFD模型处理

流场设置由ANSYS Workbench中FLUENT模块实现,流场域计算抑制整体计算域中结构场部分。叶片表面是一个无滑移剪切条件的旋转壁面,也是两个求解器之间传递数据的流体域和结构域的边界。用非结构化网格对区域进行离散化,并在叶片表面添加更精细的网格和额外的膨胀层,以提高分辨率。在Meshing中设置对称壁面、地面单元尺寸为2.719 m,内旋转域设置体网格大小为0.1 m,采用曲率捕捉划分方式,定义曲率法向角为18°。流场网格划分如图3所示。

图3 流体域边界层细节网格

上游速度入口边界,距离整机10 m,设置为自由流风速,湍流动能强度为0.05,湍流黏度比为10。下游压力出口边界,距离整机30 m,压强为一个标准大气压。外域固定对称壁面、地面和内部旋转域叶片表面设置不可滑移边界,应用Frame motion创建旋转域相对于静止域绕z轴旋转,旋转轴初始坐标为(-4.956 7×10-3,2.541,0.753 5),定义转速为200 r/min。由于模拟过程中叶片与旋转域处于相对运动的状态,对叶片交界面采用同于旋转域的壁面运动方式,选择运动壁面,运动类型为旋转,相对转速为0。

2.3.2 有限元结构模型处理

结构场的计算设置考虑受到流场数值计算压力数据,忽略流场体积部分,网格划分前,为风轮和塔筒分别赋予材料属性,如表1所示。

表1 风力机材料特性参数

划分叶片网格时,为了避免叶片前后翼缘转角尖端在划分网格时出现单元扭曲,嵌入表面网格划分单元,分别设置叶片和塔筒单元尺寸为0.02 m和0.05 m,如图4(a)所示。因为叶片表面是耦合边界,交换气动力和结构位移,采用大挠度假设,离心力考虑了叶片结构的转速,定义风轮旋转方式。求解设置添加风轮旋转轴和幅度值,选择风力机支架的底部,施加固定约束。在流固耦合和结构求解时以叶片表面为耦合界面,导入流场压力数据,最大压力值为2 448 Pa,最小压力值为-3 559.1 Pa,如图4(b)所示。

图4 静力结构分析

3 结果分析

为了评估CFD分析的收敛性,监测残差、x速度、y速度、z速度、湍流动能k和比耗散率ω,设置迭代步数2 000次,采用标准初始化方法,并从进口边界计算初始值,通过判断残差值小于1×10-6和净质量百分比小于0.1%证明计算收敛[12],残差曲线如图5所示。

图5 残差

3.1 流场计算结果分析

针对风场处于大气边界层底部,获取近地面风能资源,受到非定常来流复杂的影响问题,采用三维数值仿真还原风场模型,剖取流场典型截面的速度数值云图,分析研究风力机来流下气动特性。

如图6所示,为风力机z=1 m剖面和风轮叶片表面速度分布云图。从图6分析得:风轮剖面区域速度变化范围为0～89.25 m/s。在旋转域交界面轮廓处存在风速最小的速度梯度,形成了旋转环形流场分布,从叶尖扩展至远离交界面的位置。

图6 流场速度云图

图6(a)显示在叶片前后缘有明显的速度扩张,一直延伸到叶尖,考虑到叶片结构对气流的阻滞作用,叶片前缘梯度范围较小,跨越叶片表面气流补充过程中速度增大,即在后缘有较好的梯度过渡。图6(b)显示轮毂区域处于相对静止状态,发展到叶根开始出现速度线性增长趋势,并在叶片跨中处逐渐明显发展至叶尖。

如图7所示,为风场在x=0 m和y=2.5 m剖面处风力机速度分布云图。从图7(a)中可以看出旋转域在来流风速下产生阻塞,近风轮区域速度变化比较集中,致使风力机下游区域出现尾流拖滞,速度变化范围扩大,沿着入流方向风力机远风区风速急剧降低,最终到达初始流速。

图7 x、y剖面速度云图

图7(b)风力机周围流场诱导区风速开始下降,气流扰动的影响范围越来越大,压力梯度随风速降低而增大,至近尾流区由于叶片几何外形、塔筒和机舱的分布位置产生速度亏损,这是由于塔筒效应使得其附近流场发生变化,引起塔架后方风速降低。

图8为计算域z=1 m剖面处风力机旋转域压力分布云图和叶片轮毂表面压力云图。从图8(a)中可以看出风轮旋转剪切空气时,顺着叶轮旋转方向,叶片下表面开始产生负压力,叶片跨中至叶尖压力梯度开始增大,考虑到沿叶片相对展向旋转速度增大,叶片上下表面压力差导致气流补充交换速度增大。在旋转域内存在明显的压力变化,旋转域周围随着径向距离的增大气流扰动作用减弱并趋于稳定。表明正面入流风对风轮结构体周围压力影响较小。

图8 流场压力分布云图

图8(b)中叶片迎风面负压范围为-3 160～-217 Pa,叶片跨中部分开始有明显的负压区产生,同时在叶片表面有负压集中而不连续,因风轮受到垂直入流,叶片前缘最先接收到风速扰动,掠过叶片上下表面后逐渐恢复至稳定。另外叶片受到法向合力驱动,迎风面和背风面均匀分布的压力由于离心力作用不再完全附着于叶片表面,而是通过能量转换最终有所耗散,符合能量守恒定律,也证明了风力机气动特性的合理性。

3.2 风力机结构响应分析

图9(a)为流固耦合作用下风力机结构应力云图。分析可得,固定约束下整体变形以悬挑梁式分布,应力主要集中在塔筒两侧,从塔筒底部至顶部应力逐渐增大,且在与机舱连接处产生应力集中,塔顶最大应力为塔底应力的8倍。随着气流的来临叶片在定桨距角条件下,叶片前缘应力沿翼型弦长方向顺梯度发展。

图9 风力机结构响应图

风剪切作用主要是改变叶片表面应力梯度,顺时针转矩作用下,迎风侧表面应力大于背风侧,背风侧最小应力为70.058 Pa,沿展向应力先减小后增大再减小,在叶尖处突变至0。风轮应力集中区位于叶根同轮毂连接区域,入流风正向剪切造成应力主要集中区域在叶片跨中表面,且应力峰值为5.923 9 MPa,远小于叶片材料设计强度。

风力机整机应变云图如图9(b)所示。分析可知,应力应变云图分布情况基本一致,展现出应力应变呈现线性相关关系。应变主要集中区域位于叶片迎风侧,且应变峰值为8.611 6×10-5m/m,说明叶片是风机运行过程中维护的重点。

风力机停机工况和旋转工况下叶片变形曲线如图10所示。从停机位置右上方叶片开始逆时针编号叶片1、叶片2、叶片3。对比分析可得,上风区叶片1和叶片2在不同工况下总位移曲线以正指数形式增长。位移曲线起点表示叶片在叶根处位移且均不为零。在停机工况下,叶尖处变形达到峰值,上风区叶片沿展向位移变化规律基本相同,在旋转工况下也有类似的表现。下风区由于重力和塔影效应作用,叶片3整体位移沿展向先减小后增大,叶根至跨中的变形并不明显,甚至接近于零,脱离塔影重叠后位移量突增至正常水平。

图10 叶片变形曲线图

旋转工况下,叶片结构响应频率变大,总体位移较停机时有较大增加,上风区峰值增量分别达到2.143 7、0.867 4 mm。说明风轮气动力和离心力作用下显著影响运行过程中叶片变形量。同样,下风区叶片3位移先减小后增大,与停机时不同的是位移突变出现在后半跨,而后恢复至初始值。当叶片与塔筒位置重叠时,气流扰动相互干涉导致叶片变形波动,表明风轮方位角对结构动态响应有较大影响。

3.3 风轮叶片扭矩及输出功率分析

为了研究均匀入流条件下不同风速风力机输出功率的变化,设置入流工况为5、8、10、12、15、18、20 m/s。在FLUENT中分别模拟计算上述来流工况,绘制风轮扭矩收敛曲线,通过数据处理获得不同工况下风轮扭矩值M,利用下列公式得出风轮输出功率P,进而分析输出数据。

P=Mω

(9)

式(9)中:ω为角速度。

风轮在不同来流工况下风力机输出功率变化曲线如图11所示。从图11可以总结出,风场不同均匀入流下风轮输出功率呈现近似线性变化。入流风速在5～20 m/s增加过程中,风轮输出功率由68.33 kW增大到84.33 kW,在输出功率不断增长的过程中,考虑到风轮扭矩是随迎风速度的增大而增大的,风力机叶片在较大的入流风速下所捕获的气动压力也更多,由于叶片的结构特征,在流场作用下叶片表面受载使叶片形态参数发生变化,风轮升阻比增大,致使叶片的性能得到更好的体现,从而输出功率增大。