基于流固耦合的柔性翼型动态失速特性研究

戴 昱， 陈德龙， 信志强

（河海大学力学与材料学院，南京 211100）

随着时代的进步，微型飞行器的机翼、飞机的机翼以及风机叶片都在飞速发展，这对能源的利用和信息的收集都有着巨大的帮助.因此对翼型进行深入研究和改良创新是十分重要的.

如何提高翼型的气动性能是研究的核心问题.目前已有众多研究采用主动控制的方法使翼型发生变形来提升翼型的气动性能.张敬斌等［1］通过主动控制技术研究翼型尾缘部分的变形，发现翼型尾缘部分发生适当变形能够有效地提高翼型的升力系数和升阻比.而Zhuang等［2］的研究发现，在一定范围内，随着翼型尾缘部分主动变形的范围增大，翼型的气动性能也会随之提高.对于俯仰振荡翼型，张俊伟等［3］研究的结果表明，在翼型上仰过程中翼型主动变形量的增大可以抑制涡的分离.Kan等［4］则对翼型后缘的主动变形进行了进一步的研究，发现变形的角度如果发生周期性的变化能够增大升力系数并且可以延迟失速.尾缘主动变形的翼型和普通的襟翼相比存在一些差异，在这之前，普通襟翼的研究已有一些，例如贾亚雷等［5-6］的研究，而由于无缝光滑的特征，翼型尾缘的变形具有更优越的气动性能［7］.通过这些研究可以看出，对翼型施加主动控制使翼型变形从而提高气动性能的方法具有较为优越的效果.虽然主动变形控制的方法对气动性能的提升已经呈现出较为优越的效果，但它需要主动控制系统，这往往导致其结构较为复杂，被动控制则相对简单.被动控制的方式主要是对翼型采用柔性的结构，使其可以产生自适应变形从而提升翼型的气动性能.自适应变形的潜能是巨大的［8］.Macphee和Beyene［9］采用被动控制，利用流固耦合的方法对柔性翼型的气动性能展开了研究.研究表明，在垂直轴风机中柔性翼型的功率系数要大于刚性翼型，即使是低弹性模量的柔性翼型，相比于刚性翼型依然可以提高翼型的气动特性［10］.柔性翼型在低于风机的额定风速下所表现出来的性能要优于刚性翼型［11］.前人对柔性翼型的研究中还未考虑到柔性翼型的动态失速特性，而动态失速现象在风力机和各种飞行器中又是普遍存在的［12-16］.因此对于柔性翼型在动态失速现象中所表现出来的气动特性进行研究是十分有必要的.

本文的研究是在前人对柔性翼型研究的基础上，进一步考虑柔性翼型在动态失速现象中是否还具有一定的价值.因此本文基于流固耦合的方法，研究柔性翼型在不同减缩频率下的动态失速特性.通过对比刚柔性翼型的气动特性、流场特点以及分析柔性翼型的结构响应，得到柔性翼型的动态失速特性及其优势.

1 问题描述

本文采用的柔性翼型的模型如图1所示，以NACA0012为基础翼型，弦长c=1 m，并在距离翼型前端0.25 m处设置了一个0.02 m×0.06 m的矩形孔洞作为旋转轴，整个翼型除矩形孔洞外，其余部分均为柔性材料.

图1 柔性翼型模型Fig.1 Flexible airfoil model

来流介质为空气，密度取为1.225 kg/m3，来流速度为4.5 m/s，动力黏度取1.837 5×10-5Pa·s，雷诺数为30 万.参考前人对于柔性翼型材料参数的选择，本文中的柔性材料密度为100 kg/m3，泊松比为0.37，这种材料可用聚氨酯泡沫实现，具有低密度和广泛的弹性模量的特征［9］.在本次动态失速研究中，通过翼型绕矩形旋转轴做俯仰运动来模拟攻角的非定常变化，攻角的变化运动规律如式（1）所示：

其中：α0为初始攻角；αt为翼型俯仰过程中最大攻角；T为俯仰运动周期；n为俯仰运动周期的个数.前3 s 翼型固定不动以获取稳定流场，3 s后做角速度大小不变的俯仰运动.

减缩频率是翼型动态失速的重要影响因素［17］，其表达式为k=，其中ω表示振荡角速度，c表示弦长，u为来流速度.本文针对刚柔两种翼型分别设置了三种减缩频率，具体工况设置如表1所示.

表1 工况设置Tab.1 Parameters of various working conditions

2 数值方法及验证

2.1 控制方程

2.1.1 流体控制方程 对于不可压缩流体而言，流体密度为常数，质量守恒方程可写为

式中：i为下标，表示坐标轴的方向；u代表流体速度.

对于湍流情况，雷诺时均方程可表示为

式中：j与前者i同为下标；uˉ为时均流速；u′为脉动流速；t代表时间；f表示质量力；ρ为流体密度；p代表压强；μ为流体的动力黏度.而雷诺时均方程相比于N-S方程多出了一项雷诺应力，该项为对称的二阶张量，与速度梯度有关.这导致了雷诺时均方程无法封闭，需要建立湍流模型.本文将采用SST K-omega湍流模型［9，18］，该模型结合了k-omega模型和k-epsilon模型，能够精确地捕捉近壁面的流动分离.

2.1.2 固体控制方程 柔性翼型采用的是齐次和各向同性结构，可使用以下方程［19］：

其中：u为位移场；ρs为固体的密度；P为皮奥拉-基尔霍夫第一应力张量；S表示皮奥拉-基尔霍夫第二应力张量；F代表变形梯度；E代表格林应变张量；I表示特征张量；μ和λ为拉梅参数.

2.1.3 流固耦合 在耦合面上，流体模型和固体模型须满足运动学条件和动力学条件进行耦合.运动学条件是指耦合面上的速度和位移必须连续，可用以下表达式进行描述：

式中：Vf、uf、Vs、us分别为液体的速度、位移，固体的速度、位移.

动力学条件指的是耦合面上的力必须处于平衡状态，可用以下表达式描述：

式中：n为法线方向；σf、σs分别为液体、固体的应力.

流体侧和固体侧通过在耦合面上进行变量交换来完成耦合.在每一个时间步内，都有流体部分提供压力和黏性力给固体，固体部分提供位移增量和速度给流体［20］.一个时间步内的耦合求解流程如图2所示.

图2 一个时间步内耦合求解流程Fig.2 The process of fluid-structure interaction solution within one time step

2.2 数值方法验证

为了验证本文中流固耦合的流场及固体变形计算的准确性，采用Krawczyk等［11］的研究，定量对比不同来流速度下柔性翼型变形的大小和升阻力系数.计算的模型、参数等均与参考文献保持一致.即采用C型计算区域，柔性翼型以弦长1 m的NACA4412为基础翼型.本文取文献中弹性模量为2.5 MPa下的五种来流速度，分别为3.79、5.04、6.309、7.59、10.22 m/s，对比其变形大小和升阻力系数.

计算的结果与文献中的数据对比如图3所示.图中给出的变形大小及升阻力系数都是经过归一化处理的.以柔性翼型的尾缘部分在来流速度为7.59 m/s下的变形大小为0 m，其余工况下的变形大小为与其的差值，其中变形为正代表尾缘部分相对变形方向朝上，变形为负代表尾缘部分相对变形方向朝下；而升阻力系数则为与刚性翼型在来流速度为7.59 m/s工况下所得的升阻力系数的比值.可以看出本文的计算结果与文献相比整体吻合较好，可以证明本文对柔性翼型的流固耦合计算是准确的.

图3 验证算例结果对比Fig.3 Comparison of calculation results of verification cases

2.3 计算域网格划分及边界条件

计算域的左半边为半圆形，其半径为3.5 m，右半边为长8.5 m、宽7 m 的矩形.翼型的尾部与半圆形的圆心重合，整个计算域的C形网格划分如图4所示.左端设定为速度入流边界条件：u=4.5 m/s，v=0 m/s；上下设定为对称边界条件：∂u/∂y=v=0；右端出口设定为纽曼边界条件：∂u/∂x=∂v/∂x=0；翼型则设定为流固耦合边界条件.为了使边界层区域内的流场求解更加精确，壁面第一层网格y+设置为0.7，增长率取1.15，同时，为了验证网格的无关性，对计算域划分了三种不同数量的网格，网格数量分别为15 700、27 890和43 560并进行了计算.结果表明，后两者的升阻力系数最大仅相差0.5%，固体最大变形相差0.7%，而前两者升阻力系数相差23%，固体最大变形相差18%，考虑到计算资源和时间成本，本文采用网格数量为27 890的网格进行计算.对于时间步长，取库朗数为0.1时计算得出的时间步长为1.5×10-5s，满足计算所需的精度要求.

图4 计算域网格Fig.4 Grid of the computational domain

3 计算结果与分析

3.1 不同弹性模量下的固定攻角柔性翼型气动特性对比

为了便于研究减缩频率对翼型动态失速特性的影响，首先对四组不同弹性模量下的固定攻角柔性翼型进行模拟，分析不同弹性模量下的升阻比，选取使翼型具有最佳气动性能的弹性模量作为动态失速时柔性翼型的弹性模量.

如图5 所示为四种弹性模量下的升阻比.可以看出，在2.5°～17.5°攻角范围内，四种弹性模量的柔性翼型的升阻比均高于刚性翼型，这表明柔性翼型相对于刚性翼型普遍具有较优的气动性能.在本次模拟的四种弹性模量中，弹性模量为3×107Pa 的柔性翼型在5°～17.5°攻角范围内所具有的气动特性都是最优的，因此，对于俯仰振荡中的柔性翼型，弹性模量均选用3×107Pa.

图5 柔性翼型在不同弹性模量下的升阻比Fig.5 Lift-to-drag ratios of flexible airfoil under different elastic moduli

3.2 不同减缩频率下刚柔翼型的气动特性分析

柔性翼型在不同减缩频率下的气动特性具有明显的区别，且同一减缩频率下，刚性翼型和柔性翼型也存在区别.

图6 所示为刚性翼型和柔性翼型在不同减缩频率下的升阻力系数迟滞回线图.柔性翼型的升阻力系数迟滞回线除了局部具有波动外，整体变化趋势与刚性翼型保持一致，均有以下规律：当k=0.08 时，升力系数在上仰过程中一直保持较低的状态，在下俯过程中一直保持较高的状态，阻力系数在整个过程中没有出现大幅度的改变；当k=0.024时，在攻角大于7°后升力系数一直大于k=0.08时的值，但在攻角变小的过程中却突然骤降，直到18.5°时才回升，阻力系数在20°之后明显提升；当k=0.009 7时，上仰过程中升力系数在22°时发生了大幅度下降，在24°时出现回升，下俯过程中由25°至24°时有短暂的上升，随后一直处于下降趋势，阻力系数的变化趋势则与升力系数相同.

图6 刚性翼型和柔性翼型升阻力系数迟滞回线Fig.6 Hysteresis loops of the lift-drag coefficients for rigid airfoil and flexible airfoil

再通过比较同一减缩频率下，刚性翼型和柔性翼型的升阻力系数迟滞回线可以发现，在k较小时，二者无明显区别；而k较大时，即k=0.08时，柔性翼型的升力系数整体略高于刚性翼型，阻力系数整体小于刚性翼型，且阻力系数的最大值明显小于刚性翼型.由此可见，在高减缩频率下，柔性翼型的气动特性要显著优于刚性翼型，具有高升阻比和低阻力的特点.

3.3 柔性翼型的结构响应分析

在柔性翼型的俯仰过程中，柔性翼型会受到流体的作用力而发生自适应变形，柔性翼型的变形具有一定的规律.

图7 展示了柔性翼型在k=0.08 时整个俯仰过程中尾缘变形量的变化，其中上仰至25°攻角时的变形及应力在图8 中进行展示.通过图8 中展示的应力可以看出高应力主要集中在矩形固定孔洞附近及离孔洞最近的上下表面处，从孔洞至尾缘应力逐渐减小.观察图8中的变形分布规律可以明显看到，在整个过程中柔性翼型的变形主要集中在尾缘部分，其余部位的变形量接近零.变形量从翼型的中间至尾缘逐渐增大，在尾部达到最大值.从定性的角度来看，变形的分布规律与Macphee和Beyene［9］的研究结果不谋而合，但从定量的角度来分析，动态失速中的柔性翼型的变形在数值及变化规律上又有其独有的特征.通过图7可以看出在整个俯仰过程中，变形的大小随着攻角的增大（减小）快速地上下波动，这表明柔性翼型的变形十分灵活，对气动荷载较为敏感.在上仰过程中，随着攻角的增大，变形的数值整体趋势是略有降低的，而在下俯过程中却相反，随着攻角的减小，其变形的整体趋势也随之减小.上仰阶段和下俯阶段的变形差距存在一定的规律，在小攻角范围内（0°～12.5°），上仰阶段和下俯阶段变形大小的差距并不明显，而在大攻角范围内（12.5°～25°），下俯阶段的变形在数值上大约是上仰阶段的1.5倍至2倍，差距十分明显.

图7 柔性翼型俯仰过程中尾部变形量的变化Fig.7 Deformation variations of flexible airfoil trailing edges during pitching

图8 柔性翼型上仰至25°攻角的变形及应力Fig.8 Deformations and stresses of flexible airfoil pitching up to 25°angle of attack

3.4 不同减缩频率下刚柔翼型的流场分析

本小节将通过分析流场来对柔性翼型动态失速的机理以及刚柔翼型之间气动特性存在差异的原因进行研究.

图9为柔性翼型在k=0.08时上仰过程的涡量和流线.图10为柔性翼型在k=0.009 7时上仰过程的涡量和对应的流线.通过对比图9和图10中的流场可以发现，上仰至15°攻角时，二者尾流均出现了波动，尾缘上表面出现回流现象，且k=0.08的流场尾流部分已经出现涡街；上仰至25°时，二者出现了明显区别，k=0.08的流场尾流部分的涡街增强，翼型尾缘上表面出现流动分离现象，而k=0.009 7的流场已经出现涡旋从翼型表面大规模脱落的情况，翼型已经处于失速状态，这导致升力系数在上仰过程中突然下降.低减缩频率下的翼型在上仰过程出现失速现象，表明了高减缩频率下的翼型相比低减缩频率具有失速延迟的特性.

图9 柔性翼型上仰过程涡量、流线、变形及应力，k=0.08Fig.9 Vorticity contours，streamlines，deformations and stresses of the flexible airfoil during pitching up，k=0.08

图10 柔性翼型上仰过程涡量及流线，k=0.009 7Fig.10 Vorticity contours and streamlines of the flexible airfoil during pitching up，k=0.009 7

图11所示为刚性翼型在k=0.08时上仰至25°攻角的涡量，与图9（b）中相同条件下的柔性翼型相比，二者的区别十分明显，刚性翼型此时已经在上表面产生了集中强涡旋，而柔性翼型上表面仅仅是尾缘部分出现了流动分离，并未出现涡旋.图8展示了此时柔性翼型因为受到流体的作用力所产生的自适应变形，该变形会反作用于流场而改变流场的特性，此时柔性翼型的变形抑制了翼型表面涡旋的生成，对流场起到了一定的调控作用，使柔性翼型表面仅出现了流动分离.通过对流场差异性的分析可以得出，刚柔翼型的气动特性差异是由变形对流场的调控导致的.再通过对应图6（a）和（b）中上仰至25°时的升阻力系数发现柔性翼型的升力系数仅有微小提升，而阻力系数出现了大幅度降低.

图11 刚性翼型上仰至25°攻角涡量及流线，k=0.08Fig.11 Vorticity contours and streamlines of the rigid airfoil pitching up to 25°of attack，k=0.08

4 结论

本文采用流固耦合数值模拟方法，研究了柔性翼型在不同减缩频率下其升阻力系数和流场的变化规律，通过结构响应分析揭示了刚柔翼型气动特性具有明显差异的物理机理.研究发现即使在动态失速中柔性翼型依然具有较优的气动性能，这对未来翼型的发展具有一定的指导意义.具体的研究结论如下：

1）柔性翼型相比刚性翼型在气动特性上所展现的优势是大幅度减阻和小幅度增升.2）柔性翼型的被动变形主要产生在尾缘部分，变形在动态失速下的分布规律与前人研究的静态中的变形分布规律相同，但区别在于动态失速中变形的大小会随攻角的变化上下波动.3）柔性翼型良好的气动特性取决于柔性翼型变形对流场调控的结果.相比于刚性翼型，该变形能够抑制表面涡旋的产生.