结冰对扑动翼气动特性影响的计算分析

2018-08-28谢凌望张兴伟廖智豪

汕头大学学报（自然科学版） 2018年3期

谢凌望，张兴伟，廖智豪

（汕头大学工学院，广东汕头 515063）

0 引言

如今，风力发电作为开发新能源的一种重要技术手段越来越受到世界各国的重视.然而，在高纬度地区工作的风力机由于气候寒冷，风力机叶片时常发生结冰现象[1]，这对风力机的性能及安全性都有很大的影响.20世纪以来，无论是从实验上还是理论上，对于风力机叶片结冰问题的研究都取得了很大的进展.1992年，Bose等人[2]在自然条件下对水平轴风力机叶片进行结冰实验，他们观察到风力机叶片结冰大多发生在前缘和尾缘位置，且吸力面结冰较少.1998年，Jasinski等人[3]对风力机专用翼型NREL-S809在冰风洞内进行了结冰前后翼型性能参数对比实验，指出了结冰会导致翼型的气动性能下降，影响风力机叶片的性能.2007年，Clement等人[4]在不同结冰环境中对风力机叶片进行了风洞实验，得到了风力机叶片在不同结冰条件下的几何外形.2010年，邓晓湖等[5]通过求解N-S方程对风力机叶片表面结冰进行了模拟，得出了叶片结冰会使叶片提前进入失速区，对叶片翼型的气动特性造成恶劣的影响.2011年，朱程香等[6]采用了拉格朗日算法计算了翼型表面结冰的形成过程，指出了明冰对翼型的气动特性影响较大，结冰会使翼型提前发生流动分离，使得升力系数下降，阻力系数上升.然而，大部分的研究只针对一种结冰情况，对于多种结冰情况下气动特性变化的研究并不系统和全面.而且，大部分研究的对象都是静态下的翼型，对于动态情况下的翼型研究较少.虽然风力机叶片在实际工作状态下为旋转状态，并非扑动状态，但由于本文将三维的风力机叶片简化为二维的翼型，并没有办法真实地模拟叶片在实际情况下的旋转状态.而且，风力机叶片在实际工作中也会发生震动，故选择扑动状态下的翼型来模拟工作状态下的风力机叶片.本文将基于风力叶片的结冰形式，对扑动状态下的翼型进行研究，将翼型的结冰情况分为不同结冰类型、不同结冰时间以及不同结冰位置，并探究对称和非对称翼型在相同结冰情况下的气动特性差异.本文通过有限体积法全面地探究了结冰对于扑动翼气动特性的影响.

1 物理模型的建立

翼型前缘的结冰类型分为两类：明冰和霜冰，它们的几何外形如图1所示.两者的几何外形有很大的不同，当翼型前缘结明冰时，在弦长方向并没有很大的增长，结冰的区域是往上下翼面纵向发展的，其结冰形状比较复杂，典型的特征是会在上下翼面形成角状突出物，称为“羊角冰”.而结霜冰时，主要的结冰区域是往弦长的横向发展的，对弦长的增长比较明显，其结冰形状比较规则与光滑.本文采用NACA4412翼型作为主要研究对象，其结冰10 min时的明冰翼型与霜冰翼型如图2所示.另外，在其它结冰时间下的明冰翼型如图3所示.

图1 明冰与霜冰在不同结冰时间下的几何外形[7]

图2 结冰10 min时的NACA4412明冰翼型和霜冰翼型

图3 不同结冰时间下的明冰翼型

翼型真实结冰情况如图4所示，翼型结冰主要分布在前缘、靠近前缘的上下翼面以及在尾缘位置.据此，为探究在不同位置结冰对风力机叶片翼型所带来的影响，将完整结冰翼型进行简化，分为前缘结冰、上翼面结冰、下翼面结冰以及尾缘结冰四个部分.前缘结冰翼型在前面已经给出，其余位置的翼型结冰模型如图5所示.

图4 翼型真实结冰图[8]

图5 在不同结冰位置下的翼型模型

图6 结冰10 min时的NACA0012明冰和霜冰翼型

由于对称翼型NACA0012与非对称翼型NACA4412具有相同的厚度，所以本文选用其作为另一个研究对象，分别在结明冰10 min和结霜冰10 min的两种情况下，探究对称翼型与非对称翼型在相同结冰情况下的气动特性差异.NACA0012翼型在结明冰10 min和结霜冰10 min时的几何模型如图6所示.

2 计算理论基础

二维非定常不可压缩的N-S方程如下所示.

式中ui为流速，p为动压，ρ为流体密度，δij为克罗内克函数，vt是湍流粘度，k是湍流动能.本文采用k-ω湍流模型，雷诺数为135 000.本文采用一阶时间积分并运用SIMPLE算法来求解流场内压力和速度耦合的问题，在对流通量和扩散通量上均采用二阶中心差分格式.

本文的计算域与边界示意图如图7所示.本文采用三角形网格，网格总数大约为60 000个.计算域进口边界条件设定为速度入口，出口边界条件设定为流出边界，上下边界条件设定为对称边界，翼型表面满足壁面无滑移条件.其中，前缘与入口边界的距离为7倍弦长，尾缘距离出口边界为11倍弦长，弦线距离上下边界均为8倍弦长.本文采用动网格技术来实现翼型在扑动情况下的网格变形模拟.在动态的情况下，本文采用的扑动方式为简谐运动：其中A＝0.1c，f＝1.59 rad/s.扑动周期内的时均升力系数和时均阻力系数分别定义为：

本文采用NACA0012翼型进行静态下的翼型算例验证，在雷诺数为135 000的条件下，通过Fluent软件计算其升力系数，并绘制了其升力系数曲线图与实验值[9]相比较，如图8所示.可以看到，尽管模拟值的失速攻角提前了两度且峰值较小，但整体上的拟合程度还是较好的.对于扑动状态下的算例验证，详见文献[10].

图7 计算域与边界示意图

图8 实验值与计算值的比较

3 计算结果与分析

将图2中的物理模型进行网格划分后，通过Fluent软件计算后可以得到NACA4412翼型在不同结冰类型下的时均升阻力系数和时均升阻比的数值，如表1所示.从总体上看，无论是结霜冰还是结明冰都会改变翼型的气动性能，主要体现在时均阻力系数增大和时均升阻比下降上.而且从表1中可以看出翼型结明冰时其时均升阻比相比结霜冰时下降的幅度要更大，对翼型的气动性能影响更加明显.以8°攻角为例，翼型结霜冰时，其时均升阻比下降了大约59.3%，而翼型结明冰时，其时均升阻比下降了大约86.2%.通过Tecplot后处理软件可以得到在14°攻角下NACA4412翼型在不同结冰类型下的流线图和压力系数图，见图9和图10所示.从图9中可以看到，结霜冰和结明冰会改变翼型前缘的流线型，从而增大翼型的阻力.相比起干净翼型，翼型结霜冰和结明冰时都会在其前缘产生涡，而且结明冰时其前缘涡更大，对翼型的气动性能会产生更大的影响.从图10中可看出，无论是结霜冰还是结明冰都会改变翼型前缘的压力系数，导致翼型的气动性能下降.其中，结明冰时其压力系数改变更加明显，导致其气动性能最差.

表1 在不同结冰类型下的NACA4412翼型的气动特性参数

图9 不同结冰类型下的翼型流线图对比

图10 不同结冰类型下的翼型压力系数图对比

将图3中的物理模型进行网格划分后，通过Fluent软件计算后可以得到NACA4412翼型在不同结冰时间下的时均升阻力系数和时均升阻比的数值，如表2所示.从表2可以看到，在结冰类型为明冰时，结冰时间越长，结冰面积越大，对翼型的气动特性影响也比较明显.从表2还可以看出随着结冰时间的增加，翼型的时均阻力系数不断增大，时均升阻比不断减小，所以翼型的气动性能也不断下降.以12°攻角为例，相比起表1中的干净翼型，结冰时间为2 min时，时均阻力系数增加了0.02，时均升阻比减小了3.92；结冰时间6 min时，时均阻力系数增加了0.04，时均升阻比减小了9.23；结冰时间10 min时，时均阻力系数增加了0.15，时均升阻比减小了13.23.

表2 在不同结冰时间下的明冰翼型的气动特性参数

将图5中的物理模型进行网格划分后，通过Fluent软件计算后可以得到NACA4412翼型在不同位置下的时均升阻力系数和时均升阻比的数值，如表3所示.翼型表面的结冰位置包括前缘、尾缘和上下翼面，其中前缘结冰模型采用结明冰时的翼型.通过与表1中干净翼型的数据相比较，可以看到，翼型在不同位置发生结冰对翼型的气动特性产生的影响也不尽相同.从总体上看，无论在什么位置结冰都会造成翼型的时均升阻比变小，影响翼型的气动性能.其中，前缘结冰比起其它位置结冰会对翼型的气动特性造成更大的影响，其时均阻力系数最大、时均升阻比最小.通过Tecplot后处理软件可以得到在10°攻角时干净翼型和尾缘结冰翼型的涡量图，如图11所示.可以看到，当翼型尾缘结冰时，会改变翼型涡脱落的方向，使得其涡脱落的方向更加水平，导致其时均升阻力系数和时均升阻比都较小.

表3 在不同结冰位置下的翼型的气动特性参数

图11 在10°攻角下的涡量对比图

将图6中的物理模型进行网格划分后，通过Fluent软件计算后可以得到NACA0012翼型在不同结冰类型下的时均升阻比的数值，如表4所示.从总体上看，与非对称翼型类似，无论是结霜冰还是结明冰都会使得对称翼型的时均升阻比下降.通过与表1的非对称翼型NACA4412的时均升阻比相比较后可知，当结冰类型为霜冰时，对称翼型NACA0012的时均升阻比的平均下降幅度为26.6%，而非对称翼型NACA4412翼型的时均升阻比的平均下降幅度为42.5%.当结冰类型为明冰时，对称翼型NACA0012的时均升阻比的平均下降幅度为60.8%，而非对称翼型NACA4412翼型的时均升阻比的平均下降幅度为65.1%.由此可见，无论是结霜冰还是结明冰，在扑动的情况下，非对称翼型受到的影响都会更大.图12给出了NACA0012翼型和NACA4412翼型在相同结冰情况下的涡量对比图.可以看到，在扑动状态下，非对称翼型NACA4412在前缘处和尾缘处所产生的脱落涡都要比对称翼型NACA0012的要大，导致其气动性能较差.