二维射流中翼型气动特性计算与分析

随着固定翼类垂直起降飞机的发展,动力喷流中机翼的气动特性越来越受到人们关注。在零来流或小来流状态下,传统的固定翼类飞机由于升力不足,无法正常起飞。然而,固定翼类垂直起降飞机凭借动力喷流对机翼的作用,可产生额外的升力,从而实现垂直或短距起降[1-2]。因此,清楚地认识喷流中机翼的气动特性对此类飞机的安全起降来说非常重要。

动力喷流与升力面的相互作用是较为复杂的空气动力学问题,国内外学者已经进行了不少研究。龚志斌等[3]采用CFD(computational fluid dynamic)技术对外吹式襟翼的气动特性进行了计算,指出襟翼环量的增加与动力喷流对襟翼下表面的冲击和上表面的吹气加速有关。杨小川等[4]采用激励盘方法来等效螺旋桨,对比了不同组合形式下分布式螺旋桨对机翼的影响,结果表明不同形式的滑流均能使机翼的升阻力增大。

CFD类方法可以获取精细的流场细节,但需要求解整个流域,计算花费较大。因此,不少学者提出了一些简化方法来研究动力喷流对机翼的作用。Cole等[5]采用高阶自由尾迹方法对螺旋桨/机翼干扰进行了计算,结果表明该方法对轻载荷的桨有效。Veldhuis[6]采用滑流管模型来模拟滑流,研究了巡航阶段螺旋桨对机翼气动力的影响。但滑流管模型属于固定尾迹的一种,忽略了滑流的收缩和偏转,同样只能适用于自由来流速度大、桨盘载荷小的情况。还有一些学者采用半经验的喷流速度计算方法[7-8],但应用场景比较局限。

以上研究大多针对巡航阶段动力/机翼间的干扰,Shollenberger[9-10]较早对小来流情况下,二维无黏射流中翼型的气动特性进行了计算。其采用离散涡面来等效射流,基于射流内外总压差不变假设,进行射流边界位置的更新,对飞机起飞阶段喷流/机翼干扰进行了初步研究。与Shollenberger的求解思路类似,Lewis[11]采用自由流线模型,利用半无限长涡面和有限长离散涡面来等效射流,对零自由来流情况下二维射流的偏转进行了模拟。Bontempo等[12]利用自由流线模型,对均匀载荷桨盘的收缩效应进行了计算,进一步验证了自由流线模型的有效性。

真实动力产生的喷流不仅有轴向速度,往往还存在旋转速度,流动状态比较复杂。直接进行实际动力的三维计算,时间花费较大,且不易对比分析各部分速度对机翼影响的差异。另外,三维机翼的气动特性与二维翼型的气动特性密切相关,工程中常基于二维翼型的气动数据进行三维机翼特性的计算,如螺旋桨计算采用的片条理论。因此,为了解三维射流对机翼的影响机理,对二维射流中翼型的气动特性要有较清楚的认识,目前这方面的研究还有待深入。

面元法可看作是一种半解析的方法,能排除黏性的影响,方便变化网格并进行不同影响因素的研究,有利于得到一些规律性的结论。因此,参考文献[9]的求解思路,基于自由流线模型,利用涡面对射流进行建模,结合一阶面元法,本文建立了二维射流中翼型气动特性的计算方法,分析了射流中影响翼型气动力的不同因素,以期为进一步认识动力增升原理提供帮助。

1 射流建模

射流内的流体一般比环境流动具有更大的速度,射流与周围环境流动的边界属于自由剪切层,在此边界上切向速度不连续,这与涡面的特性一致,因此可以用涡面来等效射流流动。

本文建立了二维射流中翼型气动特性的计算方法,主要假设有:

1) 忽略黏性的影响,不考虑射流的卷吸[13];

2) 流动是二维且不可压的;

3) 射流内的总压保持不变。

1.1 无限长涡面

首先对有限长涡面的诱导速度计算进行介绍,设涡面的涡强γ保持不变,规定γ顺时针转为正,逆时针为负，位置在y=0,即x轴上,如图1所示。

图1 有限长涡面

根据文献[14],有限长等强度涡面对平面上一点的诱导速度的积分表达式如下所示：

当涡面长度无限大时,即x1→-∞,x2→+∞,积分可得

(3)

因此,单根强度为γ的无限长涡面诱导的速度将平面分为了两部分。涡面上方诱导速度为涡面下方诱导速度为

再来考虑两根强度相反的无限长涡面,将强度为-γ的无限长涡面置于将强度为γ的无限长涡面置于处,涡面均与x轴平行,则此时涡面诱导的速度场为

(4)

图2 上下强度相反的无限长涡面

可以看出,相距为h,上下强度相反的无限长涡面诱导的速度场等价于宽度为h的射流。因此,可采用上下强度相反的无限长涡面来等效射流。当涡强不变的情况下,射流内部的速度不随宽度h变化。

1.2 右半无限长上下涡面

沿x轴方向对称分布上下2条半无限长涡面, 半无限长上涡面起点为(xR,yR),下涡面起点为(xR,-yR),涡面轴向与x轴平行,向右直到+∞。上涡面涡强为-γ,下涡面涡强为γ,进行积分计算,可得右半无限长上下涡面对给定点(x,y)的诱导速度如下所示：

(5)

(6)

式中

当-yR<y<yR时

当y>yR时

当y<-yR时

同理可得左半无限长上下涡面诱导的速度场。可以验证,无限长上下涡面诱导的速度场就等于左半和右半无限长上下涡面诱导速度场的叠加。

1.3 本文采用的射流模型

实际翼型一般距射流出口为有限值,引入有限长无厚度壁面来刻画射流出口的壁面。另外,翼型存在于射流中,会对射流造成偏转,影响偏转部分涡面的涡强,因此引入有限长上下离散涡面来计算射流的偏转。基于以上考虑,最终采用的射流模型如图3所示。

图3 采用的射流模型

射流模型包括左半无限长上下涡面+上下壁面+离散涡面+右半无限长上下涡面。基于此射流模型,本文将对射流中翼型的气动特性进行计算。

2 计算模型介绍

2.1 射流边界条件

图4 射流偏转示意

对于单独射流,给定射流内部速度Vw,外部自由来流速度V∞,均沿x轴正向。在不受翼型影响时,射流不发生偏转,其内部的压力始终等于环境压力,可求得射流内外总压差

翼型存在的情况下,射流内部的压力和速度将发生变化。基于射流内外总压差不变假设,对于射流边界上的任一点均有

由于自由涡面不能承力,在涡面上任一位置有pin=pout,则总压差可写为

令为涡面上一点的速度,γsheet=Vin-Vout为涡面上一点的涡强。根据总压差不变条件最终可得

ΔH=ργsheetVsheet

(7)

(7)式就是本文射流边界满足的计算条件,可用于离散涡面涡强的更新。

2.2 自由流线模型

自由流线模型认为射流边界是一条流线,对于一个离散涡面,在起点坐标固定的情况下,根据当前涡面上各单元的速度,就可以对涡面节点的位置进行更新[11]。

设当前涡面上有N个单元,N+1个节点,每个单元的长度为Δsj(j=0,…N-1),每个单元中心的速度为(uj,vj),起始节点位置固定,则从第二个节点开始,各个节点的坐标可按(8)式进行更新

(8)

采用(8)式进行迭代,当计算过程中涡面节点的位置不再变化时,可认为计算收敛。

2.3 物面计算模型

参考文献[14],无厚度壁面采用lumped涡模型进行求解,基于法向速度边界条件,得到壁面上各单元的环量。有厚度翼型则采用等强度涡面元进行计算,采用切向速度边界条件,求解得到翼型上各单元的涡强。线性方程组形式如下所示

AX=b

A为影响系数矩阵,只和物面本身有关,本文计算中存在3个物面,矩阵阶数为3个物面单元数之和。X为方程解,包括壁面上各单元的环量和翼型上各单元的涡强。b为源项,指速度沿物面法向或切向的投影。对自由来流中的翼型而言,源项b只包含自由来流的影响。对本文计算的射流而言,源项b除了自由来流的速度外还需包含射流涡面对物面的诱导速度。

具体的系数矩阵组建过程可参考文献[14],整个线性方程组采用Gauss消元法求解。

2.4 广义Kutta-Joukowski定理

经典的Kutta-Joukowski定理(以下简称K-J定理)只适合于自由来流情况下单个翼型升力的计算。本文所研究的流场则包含多个自由涡和附着涡,经典的K-J定理不再适用。文献[15]将传统的K-J定理进行了推广,使之能适用于多个自由涡和附着涡的情况。参考文献[15],本文将采用广义K-J定理对射流中翼型的气动力进行计算。

自由来流V∞沿x轴正向,定义y轴正向为升力方向,x轴正向为阻力方向。设翼型上有N个单元,则翼型的升阻力可按(9)式计算

(9)

式中：Vi,x指除翼型外,流场中其他自由涡(射流边界)和附着涡(无厚度壁面)对翼型第i个面元x方向的诱导速度；Vi,y指这些涡对第i个面元y方向的诱导速度；Γi是翼型上第i个面元的环量。

3 方法验证

3.1 求解流程

确定了射流和翼型的计算模型后,方法的整个求解步骤如下所示:

1) 确定射流内外总压差,给定翼型的位置,确定离散涡面的长度和单元数,初始化涡面的位置和涡强。

若射流内部的速度为Vw,自由来流速度为V∞,则涡面上的涡强初始化为

γ0=Vw-V∞

(10)

2) 根据当前涡面的位置和强度求涡面对物面上各单元的诱导速度。

3) 求解线性方程组,更新物面上的涡强和环量。

4) 求物面对涡面,涡面对涡面的诱导速度,基于公式(8),更新涡面节点的位置。

5) 根据涡面上的速度,利用公式(7),更新涡面单元的涡强。

6) 重复步骤2)～5),直到翼型环量或涡面位置不再变化。

图5 本文方法框架

在计算过程中,左半无限长上下涡面的位置和涡强保持不变,右半无限长上下涡面的涡强不变,起点位置随离散涡面最后一个节点的坐标变化,有限长离散涡面的涡强和位置均会变化,整个计算过程采用松弛迭代进行。

本文计算包括的几何参数有射流的宽度h,无厚度壁面的长度L,翼型的弦长c,离散涡面的长度Lw。

选用NACA0012翼型作为研究对象,翼型的旋转中心定义在四分之一弦长处,在计算坐标系中对应为(xc,yc),翼型的攻角为α。计算坐标系原点定义在射流出口的中心,如图6所示。

图6 射流中翼型位置示意

为减小离散涡面上的网格量,可采用节点非均匀分布的方式,在靠近射流出口和翼型的位置上,布置较密的网格。在翼型后方,涡面上节点的间距可逐渐加大。

射流中翼型气动力系数按(11)式计算

(11)

式中：cl表示升力系数；cd表示阻力系数；cp表示压力系数。ρ∞和p∞分别指无穷远处空气的密度和静压；Vw指给定的射流速度。

本文中取ρ∞=1.225 kg/m3,p∞=101 325 Pa,射流速度Vw=30 m/s,射流的宽度h=0.16 m。

3.2 物面求解方法

选取NACA0012翼型,对自由来流中翼型进行计算,并与二维翼型计算开源软件XFoil的无黏解进行对比, 检验本文程序对无黏情况下翼型气动力计算的有效性。

对于面元法,翼型表面的压力系数根据文献[14],可表示为

(12)

式中,γ表示为翼型单元的涡强。

翼型上的离散面元N=256,攻角范围从-4°～12°,计算得到升力系数曲线如图7所示。

图7 无黏计算升力系数对比

本文计算得到的升力系数与XFoil计算的基本符合,表明本文翼型气动力计算方法是有效的。

3.3 单独射流模型

在没有翼型的情况下,对单独射流模型进行计算,验证本文建立的射流模型的有效性。

取无厚度壁面长度L=0.32 m,单元数96。离散涡面的长度Lw=4 m,单元数分别取为Nw=180,200,250和300,检验离散涡面网格量对射流速度的影响。

自由来流速度取V∞=1 m/s,射流速度保持不变。定义中点在x=0.32 m处,与y轴平行,长度为0.9h的线段,对比线段上的速度随网格量的变化如图8所示。

图8 单元数对射流轴向速度的影响

由图8可知,对涡面进行离散,会引入离散误差,因此得到的射流速度与理论值会存在差异。但随着网格数的增加,误差逐渐减小。当Nw=200时,轴向速度与理论值的误差不超过0.2%,以上结果表明本文建立的射流模型是有效的。

3.4 射流中翼型气动特性计算

由于缺少实验数据,本文将利用CFD计算结果作为参考,对不同偏角下射流中翼型的气动力进行计算,检验本文方法的有效性。

选用对称的NACA0012翼型进行计算,翼型的弦长c=0.2 m,旋转中心取为xc=0.32 m,yc=0 m。离散涡面单元数Nw=300,其余参数与3.3节单独射流计算相同。 翼型攻角从-4°～20°,角度间隔为4°。自由来流速度V∞=1 m/s,射流速度Vw=30 m/s,均沿x轴正向。

根据下离散涡面最后一个节点的位置变化来判断收敛,当前后两次迭代节点位置误差绝对值小于0.000 1时,认为计算收敛。给出4°和12°偏角时,涡面位置误差随迭代次数的变化如图9所示。

图9 位置误差收敛曲线

可以看出,经过多次迭代,涡面的位置最终不再变化,表明本文方法是收敛的。

计算收敛后,给出翼型4°和12°攻角时射流边界的位置如图10所示。可以看出,随着翼型攻角增大,射流方向的偏转也在变大。

图10 不同攻角下的射流偏转对比

CFD计算域远场边界分别采用压力入口和压力出口边界。压力入口给定自由来流的总压,压力出口给定自由来流静压。射流出口采用压力入口边界,根据射流速度计算出口处的总压,上下壁面则采用无滑移边界条件。

图11 CFD计算域

计算域网格数为4万,通过旋转翼型,生成不同攻角下的翼型网格。12°攻角时翼型网格如图12所示。

图12 CFD计算网格

CFD计算自由来流速度V∞=1.0 m/s,射流速度Vw=30 m/s。翼型采用无滑移边界条件,保持xcr=0.32 m,ycr=0 m。湍流模型分别选用S-A和SST k-w模型,计算收敛后对不同方法得到的气动力系数进行对比。由于无黏计算无法得到摩擦阻力,CFD的阻力将只取压差阻力。

本文方法以下简称为fsm(free streamline model),将采用广义K-J定理对翼型的升力和阻力进行计算。 fsm方法与CFD计算的气动力对比如图13所示。

图13 不同方法翼型气动力对比

可以看出,在20°攻角范围内, fsm与CFD计算的气动力符合较好,升力系数基本呈线性变化。

以上结果表明,本文的fsm方法可以反映出射流中翼型的气动特性,可用于射流中翼型气动力的计算。

4 射流中翼型气动特性影响因素研究

在保持射流总压不变的情况下,利用fsm方法,对翼型在射流中不同上下位置,射流宽度与翼型弦长比值等参数进行研究,分析这些参数对射流中翼型气动特性的影响。

4.1 不同上下位置下翼型气动力对比

保持旋转中心xcr=0.32 m不变,改变ycr的取值,研究翼型在射流中的上下位置对其气动力的影响。令ycr=0,±0.03,±0.05 m,取攻角α=0°,4°,8°,12°对应的升力系数和阻力系数进行对比。

图14 不同上下位置翼型气动力对比

总体上看,翼型位置越靠上,升力和阻力系数就越小。值得一提的是,虽然NACA0012是对称翼型,但在射流中,攻角为0°时翼型的升力也不一定为0。

取攻角α=12°时,翼型在不同上下位置时射流上边界的偏转情况对比如图15所示。

图15 翼型不同上下位置射流偏转对比

可以看出,相同攻角下,翼型处于射流中不同位置引起的射流偏角是不同的。翼型越靠下,射流偏转越明显,因此气动力越大。

4.2 不同宽弦比下翼型气动力对比

将射流宽度与翼型弦长的比值定义为宽弦比表示为RA=h/c。基于fsm方法,保持翼型的弦长和位置不变,通过改变射流宽度,研究宽弦比对射流中翼型气动特性的影响。

仍然采用NACA0012翼型,翼型的旋转中心固定为xcr=0.32 m,ycr=0。射流宽度分别取为h=0.16,0.32,3.2 m,对应宽弦比RA=0.8,1.6和16,得到翼型的升力系数对比如图16所示。

图16 不同宽弦比下翼型升力系数对比

其中,自由来流指面元法计算的,在自由来流中,NACA0012翼型对应的无黏升力系数曲线。可以看出,翼型在有限宽射流中的气动特性与无界的自由来流中差异较大。随着宽弦比的增大,翼型的升力线斜率不断增大,二者间的差异逐渐变小。

取攻角α=12°,给出不同宽弦比下,射流上边界的偏转情况对比如图17所示。

图17 不同宽弦比下射流偏转对比

由图17可知,随着宽弦比的增大,射流的偏转程度逐渐变小,这主要是因为射流边界与翼型间的距离变大,翼型对射流边界的影响在逐渐变小。

继续增大射流宽度,给出α=12°时,射流中翼型阻力系数随宽弦比的变化如图18所示。由图可知,当宽弦比较小时,射流中翼型的阻力随射流宽度的增大会先变大。当宽弦比RA≥16后,阻力则随宽弦比增大而不断减小。可以预见,当宽弦比趋于无穷时,翼型周围的流动将接近于自由来流,翼型的阻力也将趋于零。

图18 阻力系数随宽弦比的变化

4.3 无黏射流中翼型气动特性分析

为对射流中翼型气动特性的变化有更深的认识,本文将利用动量理论对射流中翼型的受力特性进行分析。

假设翼型距离射流出口足够远,认为翼型的存在对射流出口的动量没有影响,则二维射流出口动量保持为在自由来流为零的情况下,根据动量理论,可推导出射流中翼型的气动力与射流偏角θ的关系如(13)式所示

L=Tsinθ, D=T(1-cosθ)

(13)

定义参数代入上面升力和阻力的表达式有

可以看出,在零来流情况下,当翼型对射流出口动量没有影响时,参数F≡T,且与射流偏角θ无关。也就是说,翼型处于不同攻角下,计算出的参数F均是相同的。

定义可得到等价的无量纲等式

(14)

图19 不同攻角下的参数f对比

下面对这一等式进行验证。将第3.4节fsm方法计算出的不同攻角下翼型的气动力系数代入(14)式，计算得到参数f，如图19所示。在第3.4节中,cT=1.6,可以看出fsm方法计算的f与cT值的相对误差很小,不超过1%,表明动量理论分析的结果是有效的,射流中翼型气动力的大小与射流的偏转密切相关,本文fsm方法可以计算出这一特征。

再来看宽弦比对升力系数的影响,当射流宽度增加时,虽然射流的偏角θ在变小,但由于h增加,射流出口的动量是变大的,二者的乘积最终决定升力的变化。

对翼型阻力而言,从动量理论得到的阻力公式来看,只要射流偏角非零,翼型就会受到阻力。从广义K-J定理阻力公式来看,射流偏转后,翼型会受到向下的诱导速度(可看作下洗),因而产生阻力,这与三维有限展长机翼诱导阻力的产生类似。当宽弦比足够大时,此时翼型受到的下洗速度接近于零,因此翼型阻力也趋于零。

5 结 论

基于自由流线模型,本文发展了一种无黏射流中翼型气动力的计算方法。在射流内总压不变的情况下,研究了不同因素对翼型气动特性的影响,并根据动量理论对射流中翼型的受力特性进行了分析,得出的主要结论有:

1) 翼型的攻角、翼型在射流中上下位置对其气动力的影响,均可通过射流偏角反映出来,本文的fsm方法能够计算出这一特征。

2) 在二维有限宽射流中,翼型产生升力的同时也会带来阻力。

3) 当射流宽度给定时,射流偏角越大,射流中翼型的升力和阻力就越大。

4) 二维有限宽射流中翼型的升力线斜率小于自由来流中翼型的升力线斜率。

5) 攻角不变时,翼型在射流中位置越靠下,射流偏角越大。

6) 当射流宽弦比无限大时,翼型在射流中的气动特性就接近于自由来流中的情形。无黏情况下,二维翼型的阻力为零。