叶片作为风电机组风能捕获的核心部件，近年来呈大型化趋势发展，运达能源科技集团股份有限公司研制的陆上WD175-6000机型风轮直径达175 m，明阳智慧能源集团股份公司推出的海上MySE16.0-242机型风轮直径达242 m。叶片扫掠面积增大，使得机组发电功率增大，叶片大尺度变化影响叶尖失速和流动分离状态，导致叶片压强分布情况更为复杂。因此，气动性能分析一直是大型风电叶片设计制造领域关注的热点[1,2]。

叶素动量（blade element momentum,BEM）理论由于计算效率高，被风电科技工作者广泛采用[3,4]。Refan等[5]利用BEM理论对直径2.2 m的风轮进行气动分析，并开展风洞实验验证，结果表明BEM理论计算的压力系数、功率系数与实验结果吻合较好；文献[6]为提高风电机组的运行效率，对叶尖损失、偏斜尾迹和旋转效应等因素进行修正，建立了改进的BEM模型；Sriti[7]基于修正的BEM理论开展叶片气动优化设计研究，叶片优化后拥有更高的功率系数。BEM理论计算简单，但由于未考虑空气沿叶片展向的流动变化，计算精度难以保证。计算流体力学（computational fluid dynamics,CFD）方法采用三维建模方式，利用N-S方程和合适的湍流模型，模拟分析空气流经叶片的运动状态，其计算精度较BEM理论有较大提升[8,9]。张野等[10]对半径为1 m的风轮进行CFD计算和风洞实验，与BEM理论相比，CFD方法计算的气动载荷更接近实验数据。Kim等[11]采用CFD方法对半径为102.88 m的风轮进行气动分析，发现CFD方法与BEM理论计算的最大功率差值为6%。Plaza等[12]采用BEM理论和CFD方法分析了叶片压力和转矩，并与风洞实验数据进行比较，发现高风速时CFD的计算结果优于BEM。CFD方法因计算精度高被更多专家采用，但随着风轮尺度的增加，CFD仿真分析网格数量不断增加，计算效率变低。Younas等[13]以翼型NACA63-615为主建立了10 m长的风电叶片，采用200万有限元网格分析了压强、速度和流场分布。刘强等[14]利用BEM理论和CFD方法，比较了两者在不同风速下NREL 5 MW叶片功率与载荷特性，其中CFD网格为532万；Peric等[15]采用近1000万有限元网格分析不同流动状态下10 MW叶片气动性能，结果表明CFD方法计算结果与实验数据一致性良好。精细化CFD模型保证了叶片气动分析的收敛性与可靠性，但网格数量增加，加大了计算机负荷，延长了仿真时间，导致难以直接采用CFD方法对风电叶片进行优化设计，而缩比模型作为提高计算效率的有效手段被广泛应用[16,17]。Giahi等[18]根据几何相似准则按比例缩小2 MW风轮CFD模型，发现CFD缩比模型和相似理论计算的转矩和轴向力分布趋势相同；黄宸武等[19]以保持原模型和缩比模型叶尖速比相同为原则，保证缩比模型运动相似，比较了功率系数、推力系数和流场分布状态，揭示了叶片表面气体流动的差异性是导致气动性能不同的主要原因；李仁年等[20]以雷诺相似准则建立33 kW的缩比模型，保证了叶片表面流动状态、压力系数等与原模型相似，缩比模型计算的载荷误差小于2%。几何相似、运动相似和动力相似，使缩比模型与原模型在结构、运动、受力等方面均符合相似准则，但在实际开展缩比实验时，已有流体介质不能同时满足3个相似准则。CFD方法采用自定义方式假设虚拟存在某种流体介质，保证了叶片几何、运动和动力相似，合适的有限元网格尺寸提高了CFD计算效率，但过小的缩比模型忽略了叶片小特征对流场的影响，导致计算的转矩和轴向力与原模型存在一定差异性。

针对上述问题，本文以不同风轮尺度下BEM理论和CFD方法的气动性能精度对比分析为出发点，基于几何相似、运动相似和动力相似理论，分析缩比模型与原模型之间的几何、运动、力学参数比例尺关系，建立大尺度风轮的高精度CFD缩比模型，提高了20 MW叶片在多种工况下CFD仿真分析的计算效率。

1、叶片气动分析

1.1 BEM理论与分析

BEM理论假设风轮为轮盘状结构，叶片沿叶展方向分为若干个叶素段，相邻叶素段不计叶展方向的影响，将作用于叶素段上的气动力积分并分解为平行于旋转平面的切向力和垂直于旋转平面的轴向力，切向力产生的力矩驱动叶片围绕风轮轴线旋转，轴向力产生气动推力作用在塔架上，轴向力T和转矩Q分别为：

式中：R——叶片长度，m；ρ——空气密度，kg/m3;W——相对入流速度，m/s;N——叶片数量；c——弦长，m;Cn——轴向力系数；r——叶展方向的长度，m;Ct——切向力系数。相对入流速度W为：

式中：α——轴向诱导因子；V——入流风速，m/s;β——周向诱导因子；ω——风轮转速，r/min。

本文以美国可再生能源实验室（NREL）公布的5 MW风力机叶片为研究对象，叶片长度R为61.5 m，额定风速V为11.4 m/s，风轮额定转速ω为12.1 r/min[21]。经BEM理论计算，17节叶素段的轴向力和转矩如图1所示，轴向力和转矩整体沿叶展方向增大，第14节叶素段轴向力为30.2 kN，达到轴向力峰值，第5节叶素段转矩约160 kN·m，为该工况下转矩最大的叶素段。

1.2 CFD分析

CFD方法通过对叶片流场域网格划分、湍流模型选取和方程组求解等，模拟计算叶片流场域流态变化、压强分布等流体力学性能。本文忽略叶片间相互影响作用，利用多重参考系方法将流场分为外部域与内部域（如图2所示），内部域的旋转带动外部域流体运动。

图1 叶片展向轴向力和转矩分布   

图2 CFD流场域示意图   

针对额定服役工况，采用非结构网格、SST k-ω湍流模型和Simple算法进行CFD仿真分析及其网格无关性验证，第一层边界层厚度10-6m保证了y+≤1，并对叶片附近网格进行加密处理，无关性验证结果如图3所示。当网格数量大于480万时，随着网格数量的增加，转矩基本保持不变，700万网格计算的转矩较480万网格计算的相对改变量仅为0.46%，但计算耗时约为480万的1.5倍。考虑到计算精度和效率，本文采用480万网格开展5 MW叶片CFD仿真分析。

图3 5 MW叶片网格无关性验证   

1.3 风轮尺度效应分析

以NERL 5 MW风电叶片翼型数据为参考，假设额定入流风速V均为11.4 m/s，保持叶尖速比相等，利用式（4）计算不同功率的风轮直径，翼型弦长按照风轮直径进行等比例缩放，建立0.5、1、2、10和20 MW叶片几何模型，不同功率的风轮参数如表1所示。

式中：D——风轮直径，m;P——机组额定发电功率，MW。

表1 不同功率下风轮主要参数  

采用BEM理论和CFD方法分析不同功率风轮额定工况下叶片气动性能，计算的轴向力和转矩如图4所示，CFD分析的网格数量及其计算时间如图5所示。随着风轮尺寸的增大，BEM理论和CFD方法计算的转矩和轴向力增大，发展趋势基本一致，但两者计算的转矩和轴向力相对差值变大，额定功率小于5 MW的风轮相对差值小于3%，而20 MW风轮计算的轴向力相对差值达6.54%，转矩相对差值达8.61%。另外，随着风轮额定功率的增加，CFD仿真分析的网格数量和计算时间迅速增加，5 MW叶片CFD计算时间约为10 h，而20 MW叶片的耗时高达34 h。

图4 风轮尺度效应下BEM与CFD计算精度对比   

图5 风轮尺度效应下CFD的计算时长   

BEM理论根据翼型气动力数据计算叶片轴向力与转矩。翼型气动力数据对BEM理论在气动性能预测中有重要影响，不同功率等级风轮尺寸导致风力机运行过程中翼型雷诺数差距较大。本文针对不同功率等级叶片，对叶片第4节翼型DU40、第6节翼型DU35、第7节翼型DU30、第8节翼型DU25、第10节翼型DU21和第12节翼型NACA64的雷诺数进行计算，结果如表2所示。不同功率等级叶片同一翼型段上雷诺数不同，随着功率等级的增大，受翼型弦长的影响雷诺数增大。叶片不同翼型段雷诺数也不同，受叶片光顺性影响，叶片中后段翼型弦长沿叶展方向线性减小，对应的雷诺数相应减小。而雷诺数会影响叶片表面大气流动状态，进而改变叶片表面压强分布状况，使翼型的升力系数和阻力系数发生改变。本文利用Profili软件和表2数据，计算翼型DU40、DU35、DU30、DU25、DU21和NACA64在不同雷诺数下的升力系数与阻力系数，结果如图6所示。随着雷诺数的增大，各翼型的升力系数不断增大，阻力系数不断减小，对应的轴向力系数与切向力系数增大。BEM理论进行气动分析时，普遍以单一雷诺数值下的翼型升阻力曲线为参考，未考虑翼型弦长和风速变化等因素对升阻力曲线的影响，进而影响叶片的轴向力和转矩的计算精度。而CFD方法进行气动分析时，以流体介质、风速、转速和叶片几何形状等构建叶片三维流场域模型，采用合适的湍流模型，模拟仿真大气经过风电叶片的流动状态，获取了叶片叶展方向的压强与流速的分布特征，与大气流经叶片后的流动状态更为接近。 

表2 风轮尺度效应下不同翼型的雷诺数 

图6 不同雷诺数下的升力和阻力系数  

2、比例模型

2.1 相似理论

流体力学相似理论主要是指在至少满足一个相似性准则的前提下，分析两个流体模型的力学规律。本文通过叶片尺寸缩放、叶尖速比不变和流体介质自定义的方式，使叶片原型与比例模型空间形状几何相似、翼型段运动相似且经过翼型段的黏滞力和惯性力相似。设雷诺数为Re、动力粘度系数为μ、L为模型空间几何尺寸，k为比例模型与原型的几何比例尺，下标o为与原型相关的物理量，下标s为与比例模型相关的物理量，在几何相似下，比例模型与原型的几何比例尺k为：

根据运动相似和动力相似，且为保证经过翼型段的空气流态相似，采用雷诺准则，比例模型与原型的叶尖速比和雷诺数相等。

假设入流速度V和密度ρ不变，由于叶素段翼型形状、攻角和流态保持不变，轴向和周向诱导因子不变，经式（3）得相对入流速度W保持不变，则比例模型与原型的动力黏度系数和风轮转速的比值为：

根据BEM理论，比例模型与原型的推力系数CT、力矩系数CM、功率系数CP相同，则原型与比例模型的轴向力、力矩与功率的比值为：

2.2 比例模型构建及分析

利用相似理论，根据式（5）～式（9），以NREL 5 MW叶片为原型，分别建立几何比例尺为0.010、0.100、0.316、0.447和0.633的CFD比例模型，得到的比例模型主要参数如表3所示。  

表3 5 MW叶片不同比例模型参数  

依据表3参数建立不同尺度的几何模型并开展CFD仿真分析，当比例模型网格尺寸按原型等比例缩放时，比例模型与原型的网格数量保持不变，则CFD方法计算时间与效率保持不变。为提高计算效率，比例模型网格总量以几何比例尺为参考进行有限元网格划分、开展叶片流体仿真分析，并根据相似模型轴向力、力矩与功率的关系计算，结果如表4所示。几何比例尺逐渐减小，比例模型计算的轴向力和转矩对应减小，转换后的轴向力和转矩较原型的相对差值逐渐增大。当几何比例尺为0.316时，比例模型与原型的轴向力和转矩相对差值分别为3.49%和3.59%，比例模型仿真计算时间缩至4.33 h；当几何比例尺为0.100时，比例模型叶片长度为6.15 m，网格总量为565172个，轴向力和转矩相对差值分别为8.91%、15.63%；几何比例尺为0.010时，比例模型叶片长度为0.615 m，网格总量为59680个，而轴向力和转矩相对差值分别为40.73%、67.07%。 

表4 比例模型的计算精度与计算效率  

几何比例尺过小，比例模型的轴向力和转矩较原型差别越大，为进一步厘清比例模型与原型的差异性，开展大气流经叶片流态对比分析，绘制风轮平面相对速度分布图和叶片表面极限流线分布图，如图7、图8所示。比例模型风轮平面相对转速分布趋势基本一致，但几何比例尺为0.100和0.010时，叶尖周边区域流速明显增大，而流速变化影响叶尖及其周边压强分布，且靠近叶尖的叶素段能提供较大的轴向力和转矩，导致比例模型计算的转矩和轴向力差别较大。从图8可发现，各模型均从叶根周边区域开始产生流动分离，在翼型过渡段处流动分离加剧，与原型相比，几何比例尺为0.100和0.010的模型产生的失速旋涡更为剧烈，并在叶根尾缘分离处存在一定程度的翼展方向的横向流动，流动分离与叶根后缘横向流动的扩散范围更广。

图7 比例模型风轮平面速度分布图   

图8 比例模型极限流线分布图  

对原型与比例模型叶根周边区域极限流线分布的不同，绘制大气流经叶根周边区域的流线图，并选取0.08R、0.13R、0.18R位置处翼型截面周边横向流动速度进行分析，如图9所示。根据叶片几何形状，从叶根到叶片翼型过渡区域，翼型截面弦长由3.54 m逐渐增至4.56 m，当大气流经该段叶片时，翼型截面面积逐渐增大，使得大气受叶片表面压强作用产生水平向左的分力，又因为叶片表面大气流动速度很小，使得大气向叶根流动。当叶片以一定转速旋转时，带动叶片周边区域气体作圆周运动，气体受离心力影响产生图示水平向右的分力。同时，从叶片过渡区域到叶尖，翼型截面弦长由4.56 m按照线性变化的规律逐渐减至0.70 m，大气流经该段叶片时翼型截面面积减小，气体受到图示水平向右的分力作用，使叶片周边区域气体具有横向流动速度。从图9可发现靠近叶根附近的翼型横向流动速度分布均匀，沿翼展方向叶根尾缘分离处横向流动速度明显增大。笔者认为发生此现象主要是因为在叶根到叶片翼型过渡段，翼型形状由圆形截面过渡到狭长型DU40翼型截面，DU40翼型尾缘形状与圆形形状相差较大，翼型尾缘处结构梯度变化明显比前缘大，导致翼型尾缘分离处横向流动速度明显大于其他区域。

图9 叶根周边区域流态分布图  

当几何比例尺从1.000缩至0.010时，叶片长度从61.5 m缩至0.615 m，叶片翼型后缘厚度按照几何比例尺发生相应变化，在进行CFD网格划分时，叶片相邻流场网格尺寸由0.1 m缩至0.0095 m，当几何比例尺过小时，翼型后缘厚度较薄的区域被忽略，影响翼型后缘速度与压强分布，造成比例模型与原型在该区域流动分布特征差异性较大，进而影响叶根及其过渡区域轴向力和转矩的计算精度。因此，比例模型过小的几何比例尺因网格尺寸和网格总量的限制，忽略了叶片小特征对流场的影响，影响了叶片及其周边的速度、压强、流动分布情况，增大了叶片轴向力、转矩与原型的差异。

3、验证分析

参考NREL 5 MW叶片翼型、扭角等数据，设机组额定风速为11.4 m/s，按比例建立20 MW叶片几何模型，叶片长度为123 m。综合考虑计算的效率与精度，本文以几何比例尺k=0.316构造CFD缩比模型，对低风速8.0 m/s、额定风速11.4 m/s、高风速14.0 m/s这3种不同工况进行CFD仿真计算的对比分析。20 MW叶片原型在这3种不同工况下的转速分别为4.5、6.15和6.15 r/min，在风速为14 m/s时，桨距角为8.3°。根据相似理论计算得到缩比模型在3种工况下的转速分别为14.2、19.5、19.5 r/min。经CFD仿真计算，轴向力、转矩和计算时间如表5所示，风轮平面内速度和叶片极限流线分布如图10、图11所示。  

表5 不同风速下比例模型轴向力与转矩  

图1 0 不同风速下比例模型速度图   

图1 1 不同风速下比例模型流线图  

针对不同工况，缩比模型CFD网格总量为2.97×106，原型网格总量为1.22×107，缩比模型网格总量仅为原型的0.243倍，而缩比模型与原型计算的转矩和轴向力相对差值不大。低风速时，轴向力与转矩的相对差值分别为5.09%和5.0%；额定风速时，轴向力与转矩的相对差值分别为4.73%和5.94%；高风速时，轴向力与转矩的相对差值分别为5.78%和5.7%。低风速、额定风速和高风速工况下，风轮平面内气体相对转速分布基本保持一致。3种工况下吸力面极限流线分布形式基本相同，均从叶根相同位置开始产生流动分离，在翼型过渡段处流动分离加剧，叶片中段与后段极限流线变化相差不大。缩比模型CFD仿真计算时间约为7 h，原型约为34 h，计算时间缩短了27 h，计算效率提高了4倍。

4、结论

本文以某5 MW风电叶片翼型数据为基础，对比分析不同风轮尺度效应下BEM理论与CFD方法计算精度；利用流体力学相似理论，构建大兆瓦风电叶片CFD缩比模型，对比分析不同几何比例尺下叶片的气动性能；以20 MW风电叶片为原型，比较缩比模型不同工况下CFD模型的计算精度与效率，得到以下主要结论：

1）随着风轮尺度的增大，BEM理论与CFD方法计算的轴向力与转矩的相对差值不断增大，CFD方法计算所需的时间大幅增加。

2）基于几何相似、运动相似和动力相似理论，厘清了叶片缩比模型与原型之间的几何、运动、力学参数的比例关系，并发现过小的几何比例尺叶片气动性能差异大。

3）多种工况下20 MW叶片缩比模型计算的轴向力、转矩、相对速度、极限流线等流体力学参数基本与原型相同，计算效率与原型相比得到大幅提高。

参考文献:

[1]姚叶宸,刘兆方,黄典贵.高效低载的三维风力机叶片外形优化设计方法[J].太阳能学报,2023,44(1):257-264.

[9]李雄威,徐家豪,李庚达,等.一种新修正的风电机组尾流分析模型[J].太阳能学报,2022,43(8):260-265.

[10]张野,刘诗尧,李津津,等.风力机叶片三维速度场PIV测量及数值模拟研究[J].工程热物理学报,2022,43(10):2642-2646.

[14]刘强,杨科,黄宸武,等. 5 MW大型风力机气动特性计算及分析[J].工程热物理学报,2012,33(7):1155-1159.

[19]黄宸武,杨科,刘强,等.相同雷诺数和尖速比下缩比风力机流动特性分析[J].工程热物理学报,2013,34(3):440-444.

[20]李仁年,刘姝君,李德顺,等.基于相似理论的风力机气动特性分析[J].太阳能学报,2015,36(12):2916-2921.

基金资助:国家自然科学基金(51805163;52075164;52205098);湖南省自然科学基金(2019JJ50192);

文章来源:易子欣,高国强,戴巨川等.大尺度风轮高精度CFD缩比模型构建及其分析[J].太阳能学报,2023,44(12):341-348.