不同覆冰形态下风力机叶片翼型气动特性分析.pdf

第第 5 5 卷卷 第第 2 2 期期2 20 02 24 4 年年 4 4 月月新新能能源源科科技技N Ne ew w E En ne er rg gy y S Sc ci ie en nc ce e a an nd d T Te ec ch hn no ol lo og gy yV Vo ol l.5 5,N No o.2 2A Ap pr ri il l,2 20 02 24 4DOI:10.20145/j.32.1894.20240205 作者简介:徐琛苑(1996),女,硕士研究生;研究方向:风力机叶片气动特性,新能源电力系统规划。xuchenyuan_grid 不同覆冰形态下风力机叶片翼型气动特性分析徐琛苑1,汤斯琦2(1.国网江苏省电力有限公司句容市供电分公司,江苏 镇江 212400;2.国网江苏省电力有限公司镇江供电分公司,江苏 镇江 212000)摘要:低温环境下风力机叶片常面临覆冰的危险,研究覆冰对翼型气动特性的影响对覆冰叶片气动特性分析具有重要意义。文章采用基于有限体积法的数值模拟算法对 S809 二维翼型气动特性及覆冰影响进行了模拟分析。通过对比分析明冰与霜冰两种覆冰形态对翼型气动特性的影响程度,发现霜冰对翼型气动特性影响不大,而明冰则会严重恶化翼型的气动特性,甚至可能导致负阻力。关键词:风力机叶片;气动特性;覆冰影响;数值模拟中图分类号:TK89 文献标志码:A0 引言 随着全球变暖,极端气候发生越来越频繁,风力机冰冻灾害问题也越来越突出。叶片表面覆冰是风力机冰冻灾害最为突出的问题之一。叶片表面覆冰不仅给叶片增加了一个附加重力,还会改变叶片的几何外形,从而影响其气动特性1-3。对于 MW 级大型风力机超长柔性叶片,该问题尤为突出。叶片覆冰机理十分复杂,不仅与温度、湿度、风速及风向变化等周围环境有关,还与风力机叶片本身气动特性及其工作状态有关,目前尚无原型风力机三维叶片覆冰形成机理的理论解释,二维翼型覆冰形成机制也在研究中。根据叶片表面覆冰的几何形状,可以将叶片覆冰形态分为两类:明冰4-5和霜冰6-7。当叶片周围温度低且湿度不是很高时,大气中的小尺寸液滴在与叶片表面碰撞的过程中,瞬间冻结,形成不透明霜状,此时的覆冰形态称为霜冰6-7。霜冰形成过程中,液滴与叶片碰撞冻结时间极短,液滴碰撞冻结位置主要受叶片周围流场及叶片表面边界层流动影响,使得覆冰形态与叶片几何外形类似,对叶片气动特性影响并不是很大。当环境温度相对较高但低于冰点,且湿度很高时,大气中液滴的尺寸较大,大尺寸液滴与结构碰撞后冻结形成的透明冰型称为明冰4-5。因为温度不是很低且液滴尺寸较大,液滴与结构发生碰撞后并不能马上完全冻结,只有一部分在碰撞处发生冻结,而另外一部分则沿着结构表面发生流动逐渐冻结,其在叶片表面的流动不仅与叶片表面边界层流动有关,还与叶片运动状态,如转速、桨距角、方位角等有关,这导致明冰的几何外形非常复杂且很难预测,通常会在液滴撞击点附近形成一些尖角。明冰通常具有非流线型气动外形,它将严重恶化结构的气动特性,导致大幅度的流动分离从而使阻力大幅增加,升力大幅减小。为了准确模拟覆冰对叶片气动特性的影响,需要准确获得叶片表面的压力分布。一种方法是采用大涡模拟方法或者直接数值模拟方法,对叶片表面边界层流动进行准确解析,但解析边界层流动需要大量网格,计算成本过高。另一种方法则是采用湍流模型结合壁面函数的方法对边界层流动进行模化。不同的湍流模型对应于不同的壁面函数,例如k-湍流模型通常和标准壁面函数结合使用,增强型壁面函数则通常和 k-湍流模型结合使用。采用何种湍流模型与壁面函数则需要一定的经验。例如,文献8通过 S809 二维翼型的数值模拟发现,精确的预测转捩分离点的位置是准确模拟翼型压力分布曲线的关键。文献9认为,湍流边界层对过渡区非常敏感,针对涡粘模型对边界层流动分离预测通常有所提前的弊端,提出了一种人工减小过渡14新新能能源源科科技技区涡粘性的方法来延缓分离的发生,但该方法需要人为设定部分参数,且不同工况下参数并不相同,参数设定需要一定的经验。风力机叶片一般都是通过一系列翼型扭转堆叠而成,研究二维翼型的气动特性对三维叶片的气动研究具有重要意义。此外,叶素理论10、制动线11、制动面12等叶片气动特性计算方法都是基于二维翼型气动特性附加一些修正所建立的,因此,研究二维翼型气动特性覆冰影响分析具有重要意义。本文选取S809 翼型为研究对象,首先通过对比不同的湍流模型与边界层处理方法的模拟效果,验证数值模拟算法并确定最优湍流模型与边界层处理方法,进一步分析翼型表面不同覆冰形态对翼型气动特性的影响规律。1 数值模拟参数1 1.1 1 几几何何模模型型 本研究选用经典翼型 S809 翼型作为研究对象,其几何外形如图 1 所示。模拟时取弦长为 C=600 mm,计算域设置为矩形计算域,翼型前缘距离入口设为 10C,距离出口为 20C,翼型中心线距离计算域上下边界都为 5C。边界条件设置为速度入口、压力出口,上下边界为周期性边界条件。来流选为均匀层流来流,来流速度设为 U=51.7 m/s,对应的雷诺数为Re=UC/=2106,其中运动粘度=1.5510-6 m2/s。分别计算翼型在来流风攻角(Angle of attack,AOA)为0、1.02、5.13、9.22、14.24 和 20.15 时的气动特性。图 1 S809 翼型几何示意1 1.2 2 网网格格划划分分 采用 C 形拓扑结构的结构化网格划分方案进行网格划分。为了模化边界层流动,边界层网格划分需要满足一定要求,且对于不同的壁面模型,边界层网格划分要求并不相同,即第一层网格的高度与边界层内网格节点的数量需要满足一定要求。边界层网格划分时需要保证边界层内有足够数量的网格节点以模拟边界层的发展,一般需要在边界层内至少布置 15 个网格节点。根据壁面定律,第一层网格的位置可以通过其与边界层流动的关系来估算:y1=y+u(1)式(1)中 y1为第一层网格高度。u为摩擦速度,其与壁面剪切应力有关。在计算前并没有壁面剪切应力信息,因此,需要采用经验公式对其进行估算,经验公式为:u=0.5CfU2(2)式(2)中 Cf为摩擦系数,根据 F.M.White 算法,Cf=0.026/Re1/7。对于不同的湍流模型及壁面模型,边界层网格第一层网格高度不同,即 y+取值不同。采用高雷诺数湍流模型时,结合标准壁面函数,要求第一层网格节点布置在湍流充分发展区,即要求 y+30;采用低雷诺数湍流模型时,结合增强壁面函数,要求第一层网格节点布置在层流底层,即要求 y+5。针对本文模拟工况,当取 y+=30 时,由式(2)计算可得第一层网格高度 y1=0.2 mm;当取 y+=5,由式(2)计算可得第一层网格高度 y1=0.03 mm。需要注意的是,这里计算所得的第一层网格高度是根据经验公式估算所得,可能并不满足假设条件,一般需要根据模拟所得 y+的真实值进行调整,可以说,第一层网格高度的确定是一个重复调整的过程。来流风攻角为 5.13时的整体网格划分结果及局部网格示意如图 2 所示。图 2 网格划分示意24新新能能源源科科技技2 数值验证 通过对非覆冰状态下的 S809 翼型进行模拟分析,验证数值模拟方案的可靠性,并确定最优的湍流模型及壁面函数组合。2 2.1 1 湍湍流流模模型型对对比比分分析析 根据前人实验结果可知,S809 翼型在攻角为14.24和 20.15时处于失速区,边界层流动涉及转捩及流动分离,此时边界层流动状态的数值模拟对湍流模型、壁面函数的选择及网格的划分最为敏感。因此,通过模拟翼型在来流风攻角为 14.24和 20.15时的气动特性与流场特征,可以最为可靠的对比不同湍流模型与壁面模型组合及网格划分对模拟结果的影响。本研究共对比了 Realize k-湍流模型(以下简称“RKE 湍流模型”)、SST k-湍流模型和 Spalart-Allmaras 湍流模型(以下简称“S-A 湍流模型”)3 种湍流模型。其中 RKE 湍流模型为高雷诺数模型,可采用标准壁面函数,y+取值为 y+=30;SST k-湍流模型为低雷诺数模型,y+可取值为 y+=5;S-A 湍流模型可采用增强壁面函数,y+取值为 y+=5。通过模拟所得 y1值对网格进行验证。经过反复模拟与调整发现,当第一层网格高度取值为 y1=0.6 mm 时,对应的第一层网格无量纲高度 y+=30;而当第一层网格的无量纲高度为 y+=5 时,第一层网格厚度取值应为 y1=0.09 mm。y1=0.09 mm 时翼型上下表面不同位置处模拟计算所得 y+的结果如图 3 所示,可以看出此时的 y1取值满足要求。通常采用无量纲化的翼型表面压力系数分布曲线来描述翼型的气动特性。翼型表面压力系数计算公式如下:Cp=P-P0.5U2(3)式(3)中 P 为翼型表面压力,P为环境压力,U为来流风速,为大气密度。.图 3 模拟计算所得 y+与预测值对比不同湍流模型模拟所得翼型压力系数曲线与实验结果的对比图如图 4 所示。从图 4 中可以看出,在来流风攻角为 14.24时,3 种湍流模型模拟所得的压力系数分布曲线与实验所得压力系数分布曲线都符合的很好;而在来流风攻角为 20.15时,RKE 湍流模型模拟结果与实验结果符合最好,SST k-湍流模型模拟所得压力曲线在翼型前缘上侧有所偏低,即对边界层分离预测有所提前,而 S-A 湍流模型模型模拟所得压力曲线在翼型前缘上侧有所偏大,在翼型前缘下侧则有所偏低,即对边界层分离的预测有所延迟。综上所述,RKE 湍流模型在 3 种湍流模型中对分离点的预测效果最好。不同湍流模型对边界层分离的预测也可以从速度流场云图中看出,如图 5 所示。图 4 不同湍流模型计算所得翼型压力系数分布曲线34新新能能源源科科技技图 5 S809 翼型速度流场云图注:左侧图攻角为 14.24,右侧图攻角为 20.15 综合上述分析可得:在较小攻角情况下,3 种湍流模型对翼型气动特性模拟结果都很准确;但在攻角较大时,RKE 湍流模型相对于 SST k-湍流模型和 S-A 湍流模型对翼型的气动特性模拟相对更为 准确。2 2.2 2 数数值值模模拟拟结结果果验验证证 通过对比实验与模拟所得的不同攻角下 S809 翼型的升阻力系数,对本文数值模拟方案进行验证。基于上述对比分析,其他攻角下翼型气动特性的模拟都采用 RKE 湍流模型进行模拟。升力系数及阻力系数计算公式如下:CL=L0.5U2C(4)CD=D0.5U2C(5)式中 L,D 分别为翼型升力和阻力。模拟结果与实验结果符合的很好,验证了本文所用数值方案的可靠性,如图 6 所示。此外,从图 6 中还可以看出,在来流风攻角为 9.2220.15时,由于翼型上侧边界层流动分离的发生,随着来流攻角的增大,翼型升力系数增长变缓,在较大攻角时,翼型升力系数减小。对于阻力系数,随着攻角的增大,阻力系数大幅增加。这与翼型进入失速区升阻力随来流攻角的变化规律一致。图 6 升力系数与阻力系数随来流风攻角变化曲线3 覆冰影响分析3 3.1 1 覆覆冰冰状状态态 本节对比分析明冰和霜冰两种覆冰形态对翼型气动特性的影响。基于大量叶片覆冰现场观测数据对叶片覆冰形态进行测绘并进行一定简化,确定翼型覆冰形态,其中霜冰与叶片几何外形相似,其厚度分布为类抛物线分布,如图 7 所示,在迎风点处覆冰厚度最大,本研究中霜冰最大厚度选取为 0.05C,覆冰区域长度大致为 00.3C;明冰覆冰形态较为复杂,如图 8 所示,其尖角突出长度为 0.1C,覆冰区域大致为00.2C。图 7 S809 翼型霜冰覆冰形态图 8 S809 翼型明冰覆冰形态44新新能能源源科科技技3 3.2 2 覆覆冰冰气气动动影影响响分分析析 覆冰状态下 S809 二维翼型气动特性的模拟方案与非覆冰状态下 S809 二维翼型气动特性的模拟方案一致,即相同的计算域、相似的网格划、同样采用 RKE湍流及标准壁面函数,具体模拟参数设置此处不再赘述。霜冰覆冰形态下,模拟所得覆冰翼型压力系数分布曲线与非覆冰状态下实验所得压力系数分布曲线对比,如图 9 所示。从图 9 可以看出,覆冰状态下翼型上下侧压差略有减小。在来流攻角较小时(0、1.02和 5.13),覆冰翼型非覆冰区域(0.3C C)的压力系数分布与非覆冰翼型相同。而在来流攻角为9.22及 14.24 时,由于覆冰的影响,在翼型上侧0.3C 0.5C 区域,压力系数有所降低,在翼型上侧0.5CC 区域压力系数与非覆冰翼形一致,翼型下侧非覆冰区域则不受覆冰的影响。值得一提的是,在这两个攻角下,边界层流动分离点发生的位置也没有改变。在来流风攻角为 20.15时,非覆冰区域的压力分布系数也不受覆冰的影响,但此时流动分离点由于覆冰的存在有所提前。图 9 覆冰状态为霜冰时翼型压力系分布曲线,图中圆圈符号为实验所得非覆冰状态下压力系数分布 不同攻角下霜冰覆冰翼型绕流场速度云图如图10 所示,可以看出,覆冰形态为霜冰时,速度流场与非覆冰翼型绕流场类似,这与霜冰覆冰翼型压力系数分布曲线结果一致。此外,霜冰覆冰翼型表面边界层流动分离点的位置也可以从速度流场云图看出。从速度流场云图所得分离点位置与通过压力系数分布曲线所得分离点位置一致,且从速度流场云图可以看出,在来流风攻角小于 9.22时,没有流动分离发生,54新新能能源源科科技技 图 10 覆冰状态为霜冰时翼型速度流场云图 在来流风攻角为 9.22、14.24及 20.15时,翼型上侧发生了稳定的流动分离,但没有漩涡脱落产生,与非覆冰翼型类似。明冰覆冰形态下模拟所得翼型压力系数分布曲线与实验所得非覆冰翼型压力系数分布曲线对比,如图 11所示。从图11 中可以看出,由于明冰复杂几何外形的影响,翼型压力系数发生了剧烈变化,翼型上下侧压力差急剧减小,即明冰使得翼型的气动特性发生了严重的恶化。在来流风攻角5.13时,翼型上侧流动发生了完全分离。值得一提的是,在来流风攻角为20.15时,由于非覆冰翼型分离点位于距离前缘 0.2C 处,非常接近迎风点,所以在翼型非覆冰区域(0.2CC),明冰覆冰翼型压力系数分布曲线与非覆冰翼型压力系数分布曲线基本重合,而在覆冰区域(00.2C),明冰复杂的几何外形导致流场十分复杂,压力系数变化十分剧烈。图 11 覆冰状态为明冰时翼型压力系分布曲线,图中圆圈符号为实验所得非覆冰状态下压力系数分布64新新能能源源科科技技 明冰覆冰状态下模拟所得翼型周围绕流场速度云图如图 12 所示。从图 12 中可以看出,在来流风攻角为 0和 1.02时,明冰的复杂几何外形在明冰尖角尾流区产生了两个稳定的漩涡。当来流风攻角5.13时,由于来流攻角的增大,使得翼型下侧明冰尖角正对流线方向,翼型下侧明冰尖角后部漩涡消失,而翼型上侧则发生完全流动分离。图 12 覆冰状态为明冰时翼型速度流场云图 明覆冰翼型绕流场涡量云图如图 13 所示,从图13 中可以看出,在来流攻角为 14.24和 20.15时,在尾流区产生周期性的漩涡脱落。模拟所得明冰覆冰形态及霜冰覆冰形态下 S809翼型升阻力系数随来流攻角的变化曲线如图 14 所示。从图 14 中可以看出,霜冰覆冰形态下,虽然覆冰使得翼型前缘上下侧压差有所降低,但覆冰增大了翼型的弦长,从而部分抵消了覆冰带来的影响,使得霜冰翼型升阻力系数与非覆冰翼型升阻力系数差别不大,且在来流风攻角较小时,霜冰的存在使得升力系数有所提高,在来流攻角为 0时,升力系数提高了36.24%。不同覆冰状态对翼型的升力系数与阻力系数影响程度如表 1 所示。同样,由于霜冰在某种意义上增加了叶片的几何长度,使得翼型表面摩擦阻力增大,从而使得霜冰覆冰状态下翼型的阻力系数有所增大。相反,对于明冰覆冰翼型,由于明冰的非流线几何外形,翼型表面边界层流动分离大大提前,流动分离也大大加剧,使得翼型的升力系数大幅减小,最大降幅达 83.38%,如表 1 所示;且由于明冰的非流线几何外形,在大攻角下,覆冰翼型表面压差阻力大大增加,导致明冰翼型阻力系数大幅增加;值得一提的是,在来流风攻角为 0、1.02和 5.12时,由于明冰尖角回流区的影响,覆冰的存在导致了负阻力的产生。图 13 覆冰状态为明冰时翼型涡量流场云图图 14 覆冰状态时的翼型升阻力曲线74新新能能源源科科技技表 1 不同覆冰状态对翼型升阻力系数影响程度来流风攻角霜冰状态升力系数变化程度霜冰状态阻力系数变化程度明冰状态升力系数变化程度明冰状态阻力系数变化程度036.24%15.62%-37.72%-466.88%1.0225.74%15.77%-42.03%-850.17%5.1318.50%25.96%-81.57%-464.01%9.22-12.78%85.94%-83.38%274.15%14.24-4.23%14.42%-73.26%448.25%20.151.79%8.94%-39.99%122.41%4 结语 本文通过 S809 二维翼型的数值模拟,对比分析了不同湍流模型及边界层处理方法对模拟结果的影响。通过和实验数据对比,发现 RKE 湍流模型结合标准壁面函数可以很好地预测翼型表面流动的转捩与分离,得到的压力系数分布曲线与实验结果符合很好;而当流动分离较大时,SST k-湍流模型对流动分离点预测有所延迟,S-A 湍流模型对流动分离点预测则大幅提前。在此基础上,本文对比分析了霜冰与明冰两种覆冰形态对翼型的气动特性的影响。主要结论如下:(1)覆冰对翼型的气动特性的影响程度不仅与覆冰形态有关,还与来流攻角有关。(2)霜冰对翼型的气动外形影响并不是很大,在小攻角时,霜冰对分离点几乎无影响,此时翼型升力系数有所增大;大攻角时,分离点有所提前,此时升力系数有所减小。所有攻角下,阻力系数都有所增大。(3)明冰复杂的非流线几何外形导致流动分离的提前与加剧,严重恶化翼型的气动特性,使得翼型的升力系数大幅减小,阻力系数大幅增加。(4)在来流攻角为 0、1.02和 5.12时,受明冰尖角后回流区的影响,明冰的存在将导致负阻力的产生。参考文献1陈彦.冷面结霜研究及风力机叶片覆冰的数值模拟D.北京:清华大学,2012.2马强,吴晓敏,陈彦.风力机叶片覆冰的数值模拟 J.工 程 热 物 理 学 报,2014(4):770-773.3郝艳捧,刘国特,阳林,等.风力机组叶片覆冰数值模拟及其气动载荷特性研究J.电工技术学报,2015(10):292-300.4ADDY H E.Ice accretions and icing effects for modern airfoilsM.Washington D.C.:Nation-al Aeronautics Administration,Glenn Research Center,2000.5张强,曹义华,胡利.翼型表面明冰的数值模拟J.航空动力学报,2009(1):91-97.6MYERS T G,CHARPIN J PF,THOMPSON C P.Slowly accreting ice due to supercooled water impacting on a cold surfaceJ.Physics of Flu-ids,2002(1):240-256.7杨胜华,林贵平.霜冰生长过程的数值模拟J.计算机工程与设计,2010(1):191-194.8SANEI M,RAZAGHI R.Numerical investiga-tion of three turbulence simulation models for S809 wind turbine air-foilJ.Proceedings of the Institution of Mechanical Engineers,Part A:Journal of Power and Energy,2018(8):1037-1048.9ROCHA P C,ROCHA H B,CARNEIRO F M,et al.A case study on the calibration of the k SST(shear stress transport)turbulence model for small scale wind turbines designed with cam-bered and symmetrical airfoils J.Energy,2016(97):144-150.10 BURTON T,JENKINS N,SHARPE D,et al.Wind Energy HandbookM.Chichester:Wi-ley,2011.11 SORENSEN J N,SHEN W Z.Numerical model-ing of wind turbine wakesJ.Journal of Fluids Engineering,2002(2):393-399.12 YANG X,SOTIROPOULOS F.A new class of actuator surface models for wind turbinesJ.Wind Energy,2018(5):285-302.(编辑 何琳)84新新能能源源科科技技Investigation on the aerodynamic characteristic of wind turbine blade airfoils under different ice cover patternsXu Chenyuan1,Tang Siqi2(1.State Grid Jurong County Electric Power Supply Company,Zhenjiang 212400,China;2.State Grid Zhenjiang Power Supply Company,Zhenjiang 212000,China)Abstract:The blades of wind turbine often face the danger of icing under low temperature environment,making it crucial to study the impact of icing on the aerodynamic characteristics of the wind turbine airfoil.This paper employs a numerical simulation based on the fi-nite volume method to investigate the aerodynamic characteristics of the S809 airfoil and the effect of icing.Specifically,this study ana-lyzes and compares the influence of two types of ice coverings,glaze ice and rime ice,on the aerodynamic characteristics of the airfoil.The finding reveal that the rime icing has negligible influence on the aerodynamic characteristics of the airfoil,whereas the glaze icing will significantly deteriorate the aerodynamic characteristics of the airfoil and may resulting negative drag under certain conditions.Key words:wind turbine blade;aerodynamic characteristics;icing effects;numerical simulation(上接第 40 页)30 唐欣薇,张楚汉.基于均匀化理论的混凝土宏细观力学特性研究J.计算力学学报,2009(6):876-881.31 张容国,盛冬发,李忠君,等.周期性复合材料力学性能的多尺度分析J.科学技术与工程,2022(36):15994-16000.32 LEONG M,OVERGAARD L C T,THOMSEN O T,et al.Investigation of failure mechanisms in GFRP sandwich structures with face sheet wrin-kle defects used for wind turbine blades J.Composite Structures,2012(2):768-778.(编辑 姚鑫)Prediction method for mechanical properties of composites materialsbased on mesoscopic model Wei Xin1,Li Hui1,Zhou Bo1,Zheng Haocheng2(1.School of Architecture and Civil Engineering,Shenyang University of Technology,Shenyang 110870,China;2.School of Mechanical Engineering,Shenyang University of Technology,Shenyang 110870,China)Abstract:Based on the mesomechanics properties,in this paper,the mecroscopic Representative volume element(RVE)model has been established and studied.By designing the single cell size and arrangement of the RVE model,and applying periodic conditions,and a prediction method of mechanical properties of composites has been presented based on the micro-model.Firstly,for the glass fi-ber/epoxy resin unidirectional composites,the transverse and longitudinal elasticity modulus and Poissons ratio has been predicted,and the influence of uniform fiber arrangement on the mechanical properties of fiber resin matrix composites was investigated.And the accu-racy of prediction method has been verified by the prediction results of this paper are compared with the experimental values in the litera-ture.Then,for the carbon fiber/epoxy resin unidirectional composites,by changing the effective diameter of the fiber,and then predic-ting the difference in mechanical properties of composites of the quadrilateral RVE model and the hexagon RVE model under the influ-ence of different fiber volume contents has been observed.The results show the more accurate prediction result of the effective modulus of the quadrilateral RVE model under the effect of lateral stretching loads,and presenting a new thought for the prediction and applica-tion of this model.The results show the more accurate prediction result of the effective modulus of the quadrilateral RVE model under the effect of lateral stretching loads,and presenting a new thought for the prediction and application of this model.Key words:representative volume element;fiber volume fraction;periodic boundary conditions;mechanical property94