0 引　言

风机是风电技术发展的核心，目前有很多方法可以提高风机性能的途径，但对于风机本身而言，没有从根本上改进与风直接接触捕获风能部件的效率，其叶片作为风机的关键部件，由不同的翼型截面构成，翼型气动特性的好坏直接影响叶片的性能，是叶片设计的关键^[1]，因此，近年来各国的研究者更加关注的是风机叶片捕获风能效率的问题。而仿生学的兴起为改进风机叶片构型提供一种全新的思路，并且近年来也有学者尝试将仿生技术应用到风机上，胡盈彬等^[2]研究仿生尾缘襟翼设计对轴流风机气动性能及其内部流动的影响，增设仿生尾缘襟翼能有效提高风机的做功能力，并扩宽风机的高效运行区。尾缘襟翼的安装角度为140°时，风机性能提升最显著，在设计流量下，全压效率提升8.51%；宋磊等^[3]利用三维数值模拟的方法对仿鱼脊阻力型风机进行研究分析发现，在叶尖速比为0.5的工况下，随着高径比的增加仿生风机的风能转化效率逐渐提高，在高径比达到3.5时风能利用效率达到峰值；王骥月等^[4]通过对小型水平轴风力发电机的仿生设计，发现仿生叶片的升力系数为标准叶片的2.19倍，输出功率提高了25.77%。Xin等^[5]将海鸥翼型的优良气动性能应用在水平轴风力发电机上，利用数值模拟的方法对仿海鸥翼型叶片的气动性能进行研究，结果表明，在不同风速工况下和标准翼型相比，可以达到提升扭矩的作用。马列等^[6]发展了一种吸力面仿鱼形叶片，通过数值模拟方法研究了叶片仿生设计对多翼离心风机气动性能的影响，与原型相比采用吸力面仿鱼形叶片的多翼离心风机风量增加了8.4%，效率提升了6.73%；廖庚华等^[7]对比试验流量-静压曲线发现，仿生风机气动性能明显好于原型风机；刘燕^[8]在充分研究座头鲸鳍结构特点后设计一种带有前缘结节的仿生叶片，分别对仿生叶片不同结节高度和不同结节个数进行比较，得出结节数越多或者前缘结节的幅值越大时，其在大攻角下的升力系数提升越显著；Serson等^[9]发现波状前缘的翼型能够延缓失速，在大迎角下仍能保持更好的气动性；Martiqua等^[10]对具有正弦曲线和光滑前缘的模型进行实验，实验结果表明前缘正弦曲线翼型能抑制升力系数的突然降低。Wang等^[11]将前缘锯齿状结构应用到了垂直轴风机叶片上，在低雷诺数下升力系数比原型叶片有所提高，提高了风机发电量。

本文基于鸟类生物体态特征进行研究，设计具有非光滑前缘凸起的仿生叶片，通过数值模拟探究仿生叶片气动性能的优劣。首先，基于雷诺时均的方法，对5 MW风机叶片的气动特性进行数值计算。然后，通过对比参考叶片的风洞试验数据和数值结果，验证数值计算方法的可靠性。最后，通过分析前缘凸起结构对风机叶片流动分离与失速的影响，揭示延迟失速的作用机理，阐明失速工况下仿生叶片性能提升的根源。

1 仿生叶片几何模型设计

本文选用NREL5 MW风机叶片为标准叶片，如图1(a)所示，叶片的最大弦长为4652 mm，叶片长度为62900 mm。NREL5 MW风机标准叶片的扭转角参数参照文献[12]，由于叶片的扭角、弦长均沿叶片展向方向发生变化，因此仿生非光滑前缘的波长和振幅也随着叶片的弦长而变化，但是变化范围与提取的鸮翅膀非光滑前缘的变化范围保持一致，仿生叶片设计数据如表1所示。依据文献[13]可知翅膀前缘凸起高度与凸起的中线间距比为0.12～0.19，翅膀前缘凸起的间距与弦长比为0.03～0.46，羽翼前缘凸起高度与弦长比大部分集中在0.025～0.1。由此可得仿生非光滑前缘的波长B和波高A与弦长C的关系分别为B = (0.03～0.46)C，A = 0.5(0.025～0.1)C ，基于以上对叶片模型进行设计，如图1(b)所示；并且给出了仿生叶片在一个前缘正弦结节下的示意图，如图2所示。

通过对仿生叶片的气动特性进行数值计算，并与标准叶片对比，从而得出仿生叶片最优攻角下的气动性能，有利于后续仿生风机的进一步研究。

2 数值模型 2.1 控制方程与湍流模型

假设流体不可压，则流场的连续方程和动量方程表示为：

连续方程

$\frac{{\partial {u_i}}}{{\partial {x_i}}} = 0。$

(1)

动量方程

$\frac{{\partial {u_i}}}{{\partial t}} + \frac{\partial }{{\partial {x_j}}}({u_i}{u_j}) = - \frac{1}{\rho }\frac{{\partial p}}{{\partial {x_i}}} + v\frac{\partial }{{\partial {x_j}}}\left( {\frac{{\partial {u_i}}}{{\partial {x_j}}} + \frac{{\partial {u_j}}}{{\partial {x_i}}}} \right)。$

(2)

式中：${u_i}$、${u_j}$为速度分量；${x_i}$、x_j为位置坐标分量；${u_i}$为运动粘度系数；ρ为流体密度；p为压力。

为实现微分方程的封闭性以完成对方程的有效求解，本文采用${\mathrm{SS}}{{\mathrm{T}}^{}}k - \omega$湍流模型^[14]，其结合了传统$k - \omega$湍流模型和传统$k - \varepsilon$湍流模型的优势，在处理近壁面时，采用传统$k - \omega$湍流模型，在边界层边缘和自由剪切层时，则采用传统$k - \varepsilon$湍流模型，提高了计算的准确性。${\mathrm{SS}}{{\mathrm{T}}^{}}k - \omega$的湍流模型公式如下：

$\rho \frac{{{\mathrm{d}}k}}{{{\mathrm{d}}t}} + \nabla \cdot (\rho k\bar v) = \nabla \cdot \left[ {\left( {\frac{{{\mu _t}}}{{{\sigma _k}}} + \mu } \right)\nabla k} \right] + {G_k} + {S_k} - {Y_k}，$

(3)

$\rho \frac{{{\mathrm{d}}\omega }}{{{\mathrm{d}}t}} + \nabla \cdot \left( {\rho \omega \bar v} \right) = \nabla \cdot \left[ {\left( {\frac{{{\mu _t}}}{{{\sigma _\omega }}} + \mu } \right)\nabla k} \right] + {G_\omega } + {S_\omega } + {D_\omega } - {Y_\omega }。$

(4)

式中：$k$为湍流动能；$\omega$为单位耗散率；$\bar v$为平均流体速度；${\mu _t}$为湍流粘度；${\sigma _k}$和${\sigma _\omega }$为模型系数；${G_k}$为速度梯度引起的湍动能生成项；${G_\omega }$为$\omega$方程生成项；${S_k}$和${S_\omega }$为源项；${D_\omega }$为正交发散项；${Y_k}$和${Y_\omega }$为由于湍流运动引起的$k$和$\omega$的耗散项。

2.2 计算域的选取及网格划分

运用STAR-CCM+软件建立计算模型，选择了一个平行六面体的计算域选取方案，该域的长度、宽度和高度分别为7L、6L和L，其中L为叶片展长。叶片前缘位于横向中间，纵向距离速度入口2L。在进行网格划分时，先给出计算域整体网格如图3(a)所示。图3(b)为Z＝32 m截面网格图。为了让叶片表面附近有精细网格分辨率，需要在小尺度上解决涡结构和流动分离问题，在叶片表面增加了5层网格尺寸逐渐增大的棱柱层，表面网格的增长率设置为1.3，通常在实际流动中，粘性底层的气体速度差异很大，速度梯度大，粘性力较大。要准确求解粘性底层，必须将壁面附近的网格分的足够细，取y+值在1以下，第一层网格高度保持在为5.932×10⁻⁵ m附近。图3(c)为靠近叶片表面的棱柱层网格示意图。

2.3 边界条件

本文计算模型的整体计算域包括：入口边界、出口边界和计算域侧边界（包含顶部和底部共4个面），如图4所示。在入口边界，设定均匀来流进入的速度入口；在叶片表面，由于风机工作环境影响，将叶片设定为无滑移壁面条件；计算域侧边界设为滑移壁面条件；出口边界设定为压力出口条件。

3 计算结果分析 3.1 数值方法验证

因DU25翼型为5 MW风机标准叶片重要组成部分之一，故选取其作为研究对象，叶片弦长为4007 m，展长为弦长的17倍，即68119 m。在Solidworks建模后将模型导入STAR-CCM+中进行数值模拟。图5为DU25翼型数值结果与实验数据^[15]的升力系数与阻力系数比较。升力系数C_L和阻力系数C_D定义为：

${C_L} = \frac{L}{{0.5\rho {V^2}A}}，$

(5)

$C_{D}=\frac{D}{0.5 \rho V^{2} A}。$

(6)

在小攻角（0°~12°）下，数值计算结果与实验结果吻合较好，大攻角误差略有增大，但最大相对误差仍在10%以内，计算方法有效性可以确定。为了在保证精度的同时使得网格数尽可能的少，提高计算的效率，需对网格进行无关性验证，如表2所示。

其中，M1网格数为268万，M2网格数为482万，M3网格数为915万，当网格数达到460万，误差相对较小，计算结果已经稳定，因此选择该网格数计算。

3.2 仿生叶片气动性能比较

为了研究仿生叶片的气动力性能，本文在基于鸟类飞行时间，将风机叶片工作风速设置为8 m/s，雷诺数为3.2×10⁷。通过数值模拟得出标准叶片和仿生叶片在攻角在0°~39°的升阻力系数及升阻比。

图6为标准叶片与仿生叶片气动性能对比图。从图6(a)可知标准叶片随着攻角的增大，升力系数呈增大趋势，随着攻角增大到24°，升力系数则呈下降趋势，最大下降率为18.1%，而仿生叶片在0°～25°攻角时升力系数呈上升趋势。尤其在大攻角下，仿生叶片的升力系数明显大于标准叶片，且在攻角25°以后，仿生叶片升力系数最大下降率为11.5%。相较于标准叶片，前缘凸起结构可以有效延缓失速，对失速后的叶片气动性能也有一定改善作用。从图6(b)中可知随着攻角的增大，标准叶片与仿生叶片阻力系数都呈现逐渐递增的趋势，并且在攻角为27°以后，两者上升趋势更快，并且在攻角为39°时差距最大。从图6(c)可知，仿生叶片与标准叶片的升阻比变化趋势大体相同，仿生叶片与标准叶片在12°攻角时达到最高点，之后开始减小，且仿生叶片最大升阻比低于标准叶片，此后在15°～39°攻角下仿生叶片的升阻比都大于标准叶片，进一步表明前缘凸起结构会提高叶片气动性能。

3.3 叶片表面压强

依据升阻力系数特征曲线变化情况（见图6），分别在小攻角$\alpha = 15^\circ$和大攻角$\alpha = 30^\circ$下选取标准叶片和仿生叶片沿叶片展长方向靠近叶根部位且扭转角变化位置最大位置30%截面处、叶片中间50%、70%截面处和叶片末端90%截面处的特征截面进行分析，叶片表面压强分布如图7和图8所示。

当攻角$\alpha = 15^\circ$时，30%截面处的标准叶片与仿生叶片表面压强接近相等，但是在叶片前缘处，标准叶片吸力面与压力面所形成的压力差略大于仿生叶片；叶片中间50%截面处仿生叶片的表面压强小于标准叶片，该处位于前缘凸起位置，导致仿生叶片前缘受力稍晚于标准叶片；在50%截面处标准叶片所表现出的气动性要好于仿生叶片。在展长方向靠近叶尖部位30%截面处，根据叶片表面压强变化示意图，标准叶片前缘吸力面的叶片表面压力大于仿生叶片。原因是因为前缘凸起的存在导致气流向凹处流动，使得仿生叶片表面流线附着较少从而导致叶片表面压强小于标准叶片。在叶尖处，仿生叶片前缘的压强差大于标准叶片，这是因为90%处截面位于凹处位置，凸起处的气流向凹起处汇聚，导致此截面处压差较大。

当攻角$\alpha = 30^\circ$时，在30%截面标准叶片和仿生叶片压强变化示意图曲线所围成的面积相差不大，但是在叶片前缘处，标准叶片吸力面与压力面所形成的压差略大于仿生叶片，且标准叶片吸力面最低点压强最低为−137 Pa，而仿生叶片则为−121 Pa，这是因为由于失速以后，仿生叶片30%截面处于凹处位置，气流聚集，致使流动变得不稳定，流经30%截面的气流沿展向相互影响。在沿叶片展向50%截面处，仿生叶片的表面压强明显高与标准叶片，并且在50%截面处仿生叶片吸力面与压力面的压差大于标准叶片。但是从50%处表面压强变化示意图中发现，仿生叶片的压力系数曲线出现较大地波动，这是因为气流作用到叶片前缘凸起结构上后，气流出现不稳定的流动，气流会沿着凸起流向波谷，在靠近叶尖部位气流会产生沿叶片展向流动，这也是造成叶片表面压力较大的主要原因。叶片沿展向70%、90%截面处（即靠近叶尖部位）相比于在15°攻角下，70%截面处仿生叶片的最大压强为120 Pa，明显大于标准叶片的最大压强47 Pa，而90%截面处仿生叶片的最大压强为98 Pa，明显大于标准叶片的最大压强50 Pa。

由上述分析可知，大攻角下（$\alpha = 30^\circ$）仿生叶片的气动性表现出明显的优势，在展向50%、70%、90%截面处，仿生叶片由于前缘凸起结构的存在使得叶片获得表面压力大于标准叶片，从而使叶片获得更好的气动性。

3.4 速度云图和压力云图

下面给出各个截面速度云图和压力云图的比较，意在进一步分析2种叶片在流场中变化，从而更加清楚直观地了解其气动性能变化。

由图9可知，2种叶片各截面展向速度的对比，在攻角为15°时30%截面处2个叶片速度变化不大；50%截面处仿生叶片低压区明显大于标准叶片，与图7中截面仿生叶片表面压力围成面积大于标准叶片相对应；仿生叶片70%截面处可以看出在弦长的40%左右处边界层发生分离，叶片吸力面在边界层分离后的区域速度基本为0，形成失速状态，且仿生叶片失速区域明显大于标准叶片，说明在此攻角下仿生叶片的气动性不如标准叶片。90%截面处边界层分离发生在弦长60%处，相比于此时的仿生叶片边界层分离程度则更严重，使得仿生叶片的气动性不如标准叶片。

从图10标准叶片与仿生叶片附近速度云图可看出，在30%和50%截面处标准叶片与仿生叶片速度云图变化不大；70%截面处标准叶片相比与仿生叶片来说低速区范围都增大且尾流流场变化明显，整体流动状态呈现恶化状态。在基线切面处压力面靠近尾缘处速度也随之增大，在90%截面处标准叶片低速区域小于仿生叶片。

图11为风速为8 m/s攻角为15°时标准叶片与仿生叶片附近压力对比图，标准叶片与仿生叶片吸力面的压力分布几乎相等，其中在叶片前缘位置的压力分布较大，尤其越靠近叶尖的区域，而沿叶片弦向，从叶片前缘到尾缘呈现较大顺压梯度变化，且沿叶片展长方向压差变化趋于平缓；图12仿生叶片吸力面靠近上半部分的前缘凹处位置出现明显的低压区，叶片前缘的低压分布面积最大，且30%截面处仿生叶片压差大于标准叶片然而叶片其余截面中，50%截面处标准叶片与仿生叶片的低压分布面积相差不大，而70%、90%截面处标准叶片吸力面低压分布面积大于仿生叶片；从叶片前缘低压区的展向分布范围来看，标准叶片的前缘位置低压分布集中在（30%截面～50%截面）范围，而仿生叶片低压分布在叶展的（30%截面～70%截面）范围，并且仿生叶片前缘低压分布区域明显大于标准叶片。此外，虽然在70%、90%截面处标准叶片沿弦长分布低压区域大于仿生叶片，但整体气动性能表型不如仿生叶片。

3.5 涡量图和流线图分析

下面分析升力系数处于上升阶段的攻角（$\alpha = 15^\circ$）以及完全失速对应的攻角（$\alpha = 30^\circ$）时，2种叶片表面涡结构图和流线图比较，意在更进一步的体现叶片在流场中的流动情况。

由图13可知，当攻角为15°时，标准叶片除去叶根和叶尖部分由流体碰撞而产生的涡，叶片表面均保持着附着流动，并且叶片周围除叶根部流线少许紊乱外其余位置的流线分布也较为均匀；反观仿生叶片由于前缘结节的分流作用，流线由两侧凸起向凹处汇集，同时叶根尾部流线也有紊乱的现象，因此凹处更容易发生分离，其前缘结节凹处的尾缘均率先产生分离涡，叶片前缘为直线的部分流动则比较均匀。则当攻角为15°时，标准叶片的气动性能要优于仿生叶片。

由图14可知，当攻角为30°时，标准叶片的流动分离现象十分明显，从图中可知尾缘涡开始大面积脱落，并且叶根部分的涡脱落现象也十分严重，流线也出现紊乱，标准叶片出现剧烈失速现象；与标准叶片相对比而言，仿生叶片尾缘涡脱落面积小于标准叶片30%左右，并且表面涡结构起伏情况较标准叶片表现较轻，叶根部分的涡脱落现象相较于标准叶片也有改善。

4 结　语

本文采用商业软件STAR-CCM+对5 MW风机叶片进行建模仿真并进行气动性能分析，将最优攻角下的叶片装入风机进行数值模拟，与标准风机的气动性能进行对比后，所得结论如下：

1）大攻角下（$\alpha = 30^\circ$），仿生叶片在升力系数和升阻比方面比标准叶片高，升力系数最大提高为21.97%；且由不同截面压强分布可知，在$\alpha = 30^\circ$下靠近叶尖50%截面以后的升力系数要高于标准叶片，侧面说明在大攻角下仿生叶片的气动性能明显优于标准叶片。后续进一步可将仿生叶片应用于5 MW风机上，旨在提高其捕获风能的效率以及扭矩系数，在绿色能源利用方面具有重要意义。

2）由不同截面速度云图和压力云图的分析结果来看，大攻角下（$\alpha = 30^\circ$）仿生叶片靠近叶尖50%截面以后的失速区小于标准叶片，并且吸力面和压力面压差大于标准叶片。因此，在大攻角下仿生叶片相比标准叶片具有更高的升力系数和升阻比。

3）通过对比不同攻角下2种叶片在流场中的涡量和流线可知，大攻角（$\alpha = 30^\circ$）下仿生叶片与标准叶片相对比而言尾缘涡脱落面积小于标准叶片30%左右，且表面涡结构起伏情况较标准叶片表现较轻，叶根部分的涡脱落现象被有效改善。