基于数值方法研究带小插片的翼型的流动机理

吴亚东，欧阳华，许 坤，滕金芳，杜朝辉，

（1.上海交通大学航空航天学院，上海 200240；2.上海交通大学机械与动力工程学院，上海 200240）

0 引言

随着石油价格不断上涨及全球气温变暖，考虑到未来能源短缺和环境影响，开发和使用新的能源成为世界各国重点关注的领域。从世界范围看来，风能是最具发展潜力的新能源，风能发电在全球范围内呈快速发展态势。世界各国风机制造商均瞄准了利润丰厚的风机市场，风机制造商生产的风机越来越大。然而，大型系统意味着增加了载荷和系统质量，这也意味着增加了硬件和工程成本。因而，降低风力机发电成本变得尤为重要。提高风力机气动性，是大幅降低风力发电成本的最有效措施之一。流动控制技术可以使叶片工作效率显著提高。流动控制技术，是通过局部区域输入较少的能量，获得非局部或全局性的流动变化，进而改变飞行器的性能。其优点是结构小，重量轻，控制方便。在叶片尖端安装微型小插片，通过增大叶片前后面的压差，提高风力机发电输出功率，目前在风力机叶片上已有初步应用。

小插片的应用，van Dam 的一系列研究奠定了深厚的基础。2000年，Yen［1－2］以机械设计简单、低能耗和有效为出发点，提出微型小插片概念。微型小插片设备的提出，是气动载荷控制一个切实可行和有效的应用。插片的出现改变了部分弦长的局部流动以及尾迹流动，进而影响翼型的空气动力特性。Stand－ish［3］继续这项工作，进行非常全面的二维计算研究，研究了GU25－5（11）8翼型和S809翼型上下表面的小插片高度和位置。同时性能与Gurney机翼结构得到结果大致相同。同样地，插片最优高度与边界层厚度相当，但不超过边界。Chow［4］和Baker［5］对带微型小插片的翼型进行了二维计算，并与风洞实验结果进行对比，得到良好的结果。Mayda［6］进行半无限机翼模拟有限宽度微型小插片三维数值计算，结果证明插片在间距增加的情况下效能降低。国内，仅有西北工业大学的郝礼书［7］，对带微型小插片的NREL S814风力机叶片翼型进行数值模拟计算，结果证明微型小插片明显增加了升力，但是也会导致阻力有不同程度的增加。

以上的研究对带小插片的翼型在不同攻角下的流动特性分析较少，本文针对风力机叶片翼型S809加载小插片前后进行详细的数值模拟研究，并对加载小插片后的翼型在不同攻角下的流场进行了详细分析，总结小插片对翼型气动特性的影响，进而提出小插片影响翼型流动、提升翼型升力的机理。

1 数值方法

1.1 计算网格

计算用翼型采用美国可再生能源实验室（NREL）设计的风力机叶片翼型S809，弦长取50cm。数值计算使用的网格在ICEM CFD 中生成，整个计算区域分为两个部分，包含翼型的内部计算区域以及外部计算区域，如图1所示。外部计算区域的半径为20倍弦长。内部计算区域随着攻角的改变可以旋转，这样便于调整不同攻角的计算，提高效率。

内部计算区域的计算网格采用结构化网格和非结构化网格混合的网格结构，如图1（b）所示，翼型采用C型结构化网格，采用非结构化网格过渡到外部计算区域，可以使边界层网格不会延伸到远场，造成很大的长宽比。外部计算区域采用结构化网格。整个计算网格数目为74198。图1（c）带微型小插片的S809翼型的局部网格示意图。微型小插片位于翼型吸力面95%弦长位置，与翼型表面垂直，其高度为1%弦长。带小插片翼型的计算区域与base翼型相同，整体计算网格数目为124440。

图1 计算区域以及计算网格Fig.1 Computation domain and mesh

1.2 边界条件设置

进口给定速度入口边界，数值为10m／s，出口给定大气压力边界，翼型表面为无滑移绝热边界。参考压力、温度设定为101325Pa和288.15K。收敛残差设定为最大残差1×10－5。

计算基于CFX 软件平台，对流项采用二阶精度格式，湍流模型采用剪切应力输运（Shear Stress Transport，SST）模型。SST 模型是为了弥补 模型的不足而发展起来的。实质上SST 模型是一种分区模型：在边界层的不同层和自由剪切层中采用不同的二方程模型；相对于传统的涡粘性系数的公式，SST模型还对其进行了修正以改进该模型在逆压力梯度流动总的预测性能，因此SST 模型在预测剪切流、压力梯度较大的流动、分离流中表现出很好的适应性。

数值计算时对N－S方程、湍流方程等进行耦合计算求解，方程在迭代算过程中采用隐式格式。每次迭代计算后记录x、y方向速度等参数的差值，同时利用升力系数、阻力系数和扭矩系数来监测解的收敛性，要迭代计算到这些系数均收敛为止。

在小攻角情况下，定常计算能获得收敛结果，随着攻角的增大，计算出现振荡，这时需要进行非定常计算。计算时通过编译CCL 语句，监控翼型的升力和阻力，定常计算时，监控的升阻力是一定值，非定常计算时，监控的升阻力呈周期性变化，升阻力采取几个稳定计算周期结构的平均值。翼型攻角的计算范围为0～30°，每隔2°计算一个工况。

2 结果分析

2.1 升阻力特征

图2和图3给出了原始S809翼型和95%弦长位置处带小插片翼型在不同攻角时的升力和阻力结果。通过对比可以看出，带微型小插片的翼型，在不同攻角下，升力系数均有所提高。在－4°至10°攻角范围内，升力系数增加的幅度逐渐增大，而后基本稳定增加。至16°攻角时，升力开始下降，在18°攻角至24°攻角时，两种翼型的升力系数和阻力系数均出现较大的波动，24°攻角以后，变化平缓，呈略微上升趋势，带小插片的翼型比base翼型的升力要低，而阻力系数在小于16°攻角时，基本上保持不变。由流体力学知识可知，形成升力主要靠翼型吸力面和压力面的压差所造成，摩擦作用很小，所以这时的升力也是逐渐增大的。而阻力由压差阻力和摩擦阻力两部分组成，这时翼型上下面的压差还很小，总的阻力主要是摩擦阻力，所以在小于16°攻角时阻力的变化不大，攻角进一步增大至24°攻角，变化较剧烈，之后变化平缓。翼型的升阻比反应了翼型的气动效率，翼型升阻比越大，翼型的气动效率越高。图4给出了两种翼型的升阻比，从图中可以明显的看出，带小插片翼型的升阻比在小于20°攻角时，始终高于base翼型，在大于20°攻角后，升阻力基本上不会发生变化，这说明小插片起的提升升力作用是在翼型正常工作范围内，超过正常工作范围，升力反而会降低。

图2 升力系数CLFig.2 Lift coefficient

图3 阻力系数CDFig.3 Drag coefficient

图4 升阻比CL／CDFig.4 Lift drag ratio

2.2 表面压力系数

为了进一步说明升力发生的变化，需要提取翼型的表面压力系数进行分析。图5给出了不同攻角下翼型表面压力系数Cp的分布，其中攻角小于16°时，采用定常计算结果，大于或者等于16°时，采用非定常时均结果。

从图5观察可发现，2°攻角时吸力面负压绝对值较低，和压力面压力差距很小，上下表面压力分别呈狭长的细条，包含的面积很小，所以升力不高。当翼型绕流处于小攻角（0°＜α＜6°）时，压力面与升力面压力相差很小，且在翼型表面的前半部分基本是顺压分布，后半部分是逆压分布，这样的压力分布不利用边界层分离，动压能够克服后半部分的逆压静压，不产生分离。当攻角增大时，压力面正压拐点基本都处于同一位置而没有明显的变化，而在翼型的吸力面上，负压力绝对值快速升高，吸力面与压力面的压差在翼型的前缘部分明显增大，因而能够给翼型提供较高的升力。不过，与此同时吸力面上的逆压分布逐渐向翼型前缘移动，使整个吸力面基本都是逆压分布，且逆压强度也逐渐增大，这样气体流过翼型吸力面时所克服的逆压增大，当动压不能克服逆向静压时，流体不再附在壁面上，边界层发生分离，出现回流，形成漩涡。

从翼型表面压力分布还可以发现，随着攻角的增大，翼型吸力面和压力面上的压力系数之间的差距逐渐增大，即两者所包含的面积逐渐增大，所包围的面积越大，则表明做功能力越强。

从图5中可以看出，0°～16°攻角，带小插片翼型和base翼型在0～95%弦长范围内，压力系数的分布趋势接近，尤其是靠近前缘位置处，分布差别不大。从小插片位置至翼型尾缘，这段的压力分布与原始翼型有较大差别，带小插片翼型的压力分布有两段封闭的环面，在尾部单独形成一个环面，这说明小插片的存在使翼型尾缘形成一个压力的二次分布，小插片造成翼型尾缘的压力重新分布，从而造成翼型中部弦长至尾缘部分的压力出现较大的增加，使压力分布的包络面面积增大，致使产生更大的升力。

在18°攻角和20°攻角时，带小插片翼型的变化趋势与原始翼型有较大差别，这和前面分析的在18°～24°攻角波动较大相吻合。24°攻角的压力系数分布，两种翼型又呈接近趋势，而且数值上也较为接近，只是尾缘部分有小插片带来的差别，说明这时候的升力相接近，这和上面分析的升力分布相吻合。

2.3 速度流线图

图6和图7为0°～20°攻角时，原始S809翼型和95%弦长处带微型小插片的翼型的流线图，其中0°～14°为定常计算结果，18°、20°时给出了瞬态结果。从流线分布上分析两种翼型流动结构上的区别。

从原始S809翼型的流线图可以看出，当0°＜α＜10°时，绕翼型的流体紧贴于壁面，没有出现漩涡；当攻角为10°时，在翼型的尾部吸力面出现了一个小的分离涡。随着来流攻角的增大，翼型吸力面上来流脱离现象越来越严重，涡强度逐渐增大，且“滑”向后缘，最后拖向尾流，翼型后出现明显的“涡街”。继续增大攻角，在吸力面的大部分区域，绕翼型的流动不再附体，而是从表面分离出去。翼型吸力面流动分离区域集中在叶型的中后部和尾部，头部开始30%弦长以内区域的流动仍然保持较高的升力系数。当攻角增大到18°左右时，尾缘处较小的脱离涡尺寸增大，变得与主涡尺寸相当了，尾涡扩大到一定程度之后，向前缘扩散，流动分离点到达翼型的头部位置，分离的流动在翼型的吸力面附近形成了一个大的分离涡。

由带微型小插片的翼型可以看出，从0°攻角开始，翼型下表面微型小插片后部形成稳定的逆时针的分离涡，流动经过小插片后，速度降低，类似于台阶流动，在小插片后部形成低速低压区，在翼型的吸力面靠近尾缘处，速度较大，对于base翼型，其压力是小于压力面的，由于小插片的存在，降低了压力，使得在靠近尾缘处，吸力面上的压力大于压力面的压力，迫使一部分流体向压力面流动，使得翼型的后驻点位置向压力面偏移，这说明小插片改变了翼型环量的分布，从而改变了升力。微型小插片后部的压力较低，在尾缘区产生了一个向下的力，增加了翼型的力矩。随着攻角的增大，微型小插片产生的涡变得像一个分离的气泡。当微型小插片产生的诱导涡变大到能影响尾缘的时候，将开始对升力和力矩起作用。

图5 原始S809翼型和带小插片翼型在不同攻角下的压力分布系数Fig.5 Surface pressure coefficient of base airfoil and airfoil with microtab versus different angles of attack

图6 原始S809翼型速度流线图Fig.6 Velocity streamline of base S809

图7 带小插片翼型的速度流线图Fig.7 Velocity Streamline of S809with micro tab

一旦逆时针的涡及其低压区超过尾缘，将发生一个有趣的现象。一些离开尾缘的吸力面的气流，将夹带进入到压力面微型小插片产生的涡当中。而且，这个流动在吸力面和压力面的压力立即开始减缓微型小插片后部低压的集结，并且增加升力和扭矩。受翼型吸力面分离涡影响的气流，由吸力面沿后翼尾缘被拉下，在压力面产生了向上的推力。受到这个涡的影响，这股气流继续沿着压力面向上，沿着微型小插片的吸力面向下流动到微型小插片最低端。这股气流继续沿着微型小插片最低面流动，直到最低面气流流动到微型小插片的压力面。这样，在这个新的停滞点，这两股气流离开翼型表面，从α＝6°的流线图可以看出。分离点出现在后缘到分离点出现在微型小插片末端，这种转变改变了翼型库塔条件。因而，有效地增加了气动曲面轮廓和环量。这样，绕翼型的环量是增加的，所以升力也是增加的。随着攻角的进一步增大，在20°以后，小插片带来的效果就不明显，这时候从速度流线图可以看出，小插片引入的分离涡加强了翼型吸力面上涡的分离，使得升力降低，阻力增加，没有起到有效的作用。

3 结论

本文采用数值模拟的方法，对原型S809翼型和95%弦长位置处带小插片的S809翼型进行了详细的流场分析，通过对比升力系数、阻力系数、升阻比、表面压力系数、速度流线图等，详细分析了带小插片翼型和原始翼型在不同攻角情况下的流动特征，主要结论如下：

（1）采用内部计算区域和外部计算区域相结合的网格拓扑结构，方便调整翼型计算的攻角，提高计算效率，采用SST 湍流模型能够很好模拟翼型流动以及翼型在大攻角时的分离流动；

（2）带小插片翼型提升升力作用是在翼型正常工作范围内，超过正常工作范围，升力反而会降低；

（3）从表面压力系数和速度流线图等分析，带小插片翼型提升升力的机理在于，改变了翼型的库塔条件，使后驻点位置出现在微型小插片末端，有效地增加了气动曲面轮廓和环量。这样，绕翼型的环量是增加的，所以升力也是增加的，由于小插片的尺寸较小，对翼型阻力的影响较小。

［1］YEN D T，van DAM C P，BRAEUCHLE F，SMITH R L，COLLINS S D.Active load control and lift enhancement using MEM translational tabs［R］.AIAA 2000－2422.

［2］YEN D T，van DAM C P，SMITH R L，COLLINS S D.Active load control for airfoils using microtabs［J］.JournalofSolarEnergyEngineering，2001，123：282－289.

［3］STANDISH K J.Aerodynamic analysis of blunt trailing edge airfoils and a microtab－based load control system［D］.［MS Thesis］.University of California，Davis，2003.

［4］CHOW R，van DAM C P.Computational investigations of deploying load control microtabs on a wind turbine airfoil［R］／／2007ASME Wind Energy Symposium／45th AIAA Aerospace Sciences Meeting and Exhibit［C］.AIAA 2007－1018.

［5］BAKER J P，STANDISH K J，van DAM C P.Two－dimensional wind tunnel and computational investigations of a microtab modified airfoil［J］.JournalofAircraft，2007，44（2）：563－572.

［6］MAYDA E A，van DAM C P，YEN－NAKAFUJI D.Computational investigation of finite width microtabs for aerodynamic load control［R］／／2005ASME Wind Energy Symposium／43th AIAA Aerospace Sciences Meeting and Exhibit［C］.AIAA 2005－1185.

［7］郝礼书，乔志德，宋科，等.Microtab对风力机叶片翼型气动特性的影响研究［J］.航空计算技术.2010，40（2）：24－27.