由于在军事及民用上的特殊作用, 近年来太阳能无人机受到了极大关注^[1]。包括我国在内各个大国, 均投入了对太阳能无人机的研制。从目前公开的太阳能无人机型号来看, 为了追求长时间留空的目标, 太阳能无人机一方面尽可能地增加太阳能电池板的铺设面积; 而另一方面尽可能地减轻机体的结构重量, 普遍具有大展弦比、小翼载荷的特点。其往往采用超轻面密度结构, 具有较大的弹性变形, 这也导致了太阳能无人机的气动弹性问题突出, 对突风敏感的问题。美国“太阳神号”太阳能无人机, 就是在进行燃料电池测试时在太平洋上空遭遇阵风紊流, 进而导致了机体的解体。美国NASA在对联结翼传感器飞机进行测试后指出, 飞机结构设计由阵风载荷主导, 如果飞机的阵风载荷可以主动减少到飞机定常操纵设计水平, 那么机体可以减重20%。由此可见, 阵风载荷对太阳能无人机的影响是其在结构设计过程中所必须考虑的重要一环^[2]。阵风载荷的准确计算对飞机的结构设计以及结构重量有着至关重要的影响。

随着CFD技术的发展, 目前相关学者采用CFD方法对阵风响应现象进行了广泛的研究。Dillsaver等^[3]对柔性飞机的结构刚度对阵风响应特性的影响进行了研究, 研究表明“Helios”这种大柔性飞机在阵风的影响下, 其很容易进入动力学不稳定状态。Patil^[4]分别采用了时域与频域的方法对大柔性飞机的阵风响应进行了研究。Zaide^[5]将网格速度法引入来模拟阵风影响, 并将非定常气动力降阶方法与传统CFD方法进行比较, 结果表明在亚声速与跨声速阶段, 降阶方法具有较好的计算精度。在国内, 詹浩等^[6]采用模态叠加法对弹性飞机的阵风响应进行了研究。许晓平等^[7]通过设计舵面与阵风的同期运动, 对阵风减缓方法进行了研究。综上所述, 国内外相关学者已经对阵风响应进行了比较系统的研究, 也取得了大量的成果。但是, 他们对阵风响应的研究重点主要集中在非线性结构模型的准确建立上, 所建立的气动力模型通常采取解析方法或偶极子法等简化的快速计算理论。这些简化的气动模型无法对低雷诺数下由于层流分离泡而出现的气动特性恶化、气动力非线性等现象做出模拟。在低雷诺数条件下, 翼型的阵风响应由于低雷诺数效应^[8-9]而存在载荷卸载的现象, 其升力系数的响应峰值要大大低于常规雷诺数翼型。对结构设计具有重要的指导作用。与常规雷诺数下的气动响应有很大的不同。因此很有必要采用精确的CFD方法, 对低雷诺数条件下的翼型阵风响应特性进行研究。

本文将基于此, 通过引入网格速度法模拟阵风, 对不同雷诺数以及不同迎角条件下的阵风响应展开研究, 分析雷诺数效应以及大迎角所带来的气动力非线性对阵风响应的影响。

1 数值求解方法介绍及验证 1.1 数值方法

本文采用ICEM软件对低雷诺数翼型FX63-137进行几何建模。整体采用C型结构网格, 在近壁面进行附面层加密。由于要进行低雷诺数计算, 为了精确捕捉附面层流动特征, 消除网格带来的影响, 第一层网格距离以Y⁺为0.5进行选取。翼型前后远场距离均取30倍弦长。近壁面网格如图 1所示:

本文采取有限体积法对二维可压缩积分形式的雷诺平均Navier-Stokes(RANS)方程组进行求解。在控制方程中, 采用二阶精度的Roe迎风离散格式对流通量项进行离散, 采用二阶精度的中心差分格式对黏性通量项进行离散。采用隐式LU-SGS(lower-upper symmetric Gauss-Seidel)进行时间离散。

1.2 湍流模型

在低雷诺数模拟中, 对附面层近壁面流动细节的准确捕捉是模拟准确的关键。γ-Re_θt转捩模型是Langtry和Menter基于流场当地变量提出的。其在两方程SST模型的基础上通过增加了间歇因子γ和转捩动量厚度雷诺数2个变量, 以此构造了2个新的运输方程。其中间歇γ的取值范围为0≤γ≤1, 代表了空间某点的流态为湍流的概率。用来控制γ运输方程中的间歇因子的生成, 保证方程的封闭性。主要的运输方程如下:

γ的输运方程为

(1)

Re_θt的运输方程为

(2)

与SST模型的耦合方程为

(3)

上述运输方程中参数具体取值请参考文献[10-11]。

1.3 网格速度法介绍及验证

本文通过引入网格速度法^[12]来模拟阵风的影响。与动网格方法通过物体边界的运动来获取边界速度不同, 网格速度法通过赋予网格速度, 根据相对运动的原理, 相当于给物面边界一个反方向的速度来模拟阵风的影响。且在采用网格速度法模拟阵风时, 可以采用结构网格进行计算且网格不需要实时更新。相比于动网格方法, 网格速度法具有较高的计算精度与计算效率。在进行CFD计算时, 流域中的速度可以用(4)式表示

(4)

式中，u, v, w是直角坐标系下的速度分量, 而x_t, y_t, z_t指的是网格速度在3个方向上的分量。当流域受到一个向上的阵风干扰时, 相当于在z方向有一个速度wⁱ的突增, 此时, 流场的速度可以表示为

(5)

即

(6)

从上式可以看出, 在z方向上速度的突增相当于给网格一个向下的速度。

为了验证网格速度法的准确性, 本文采用NACA0006翼型进行验证。将网格速度法与动网格方法进行相互验证, 并与文献参考值进行比较。设定翼型在Ma=0.3, 迎角为0°的情况下, 突然受到向上的w_g=0.08V_∞的阵风干扰, 使翼型的有效迎角突增0.08 rad(约4.57°)。定义无量纲时间为:S=2V_∞t/c, 其中c指翼型的气动弦长。翼型的阵风响应曲线如图 2所示。由图 2可知, 通过网格速度法的形式模拟阵风的影响与采用动网格方法的数值模拟结果完全贴合, 且与文献参考值也吻合较好, 最大误差始终保持在4%以内。表明本文采用的网格速度法能够精确地模拟阵风的影响。

1.4 低雷诺数数值求解验证

在低雷诺数条件下, 流动十分复杂, 会发生转捩-分离-再附等现象, 形成长层流分离泡。分离泡的形态, 位置等均会对翼型的气动特性产生较大的影响。对转捩位置的准确判断是模拟的关键所在。本文接下来对γ-Re_θt转捩模型在低雷诺数条件下的数值模拟能力进行验证。以低雷诺数翼型FX63-137为算例, 参考美国伊利诺伊大学(UIUC)^[13]低湍流度亚声速风洞实验数据进行对比。选取的计算状态为:基于翼型弦长雷诺数为Re=3.5×10⁵, 计算迎角范围为α=-6°~12°, 来流湍流度Tu_∞=0.1%。图 3给出了翼型升阻力计算结果与实验值的对比结果。

从图 3可以看出, 计算所得的升、阻力系数与实验结果吻合较好, 特别是升力系数模拟结果, 与实验值误差始终保持在5%以内。总体来说, γ-Re_θt转捩模型对低雷诺数翼型绕流流场具有较好的模拟能力。表明本文采用的数值模拟方法能够准确捕捉低雷诺数近壁面的流场细节, 且具有较高的模拟精度。

1.5 阵风计算模型介绍

本文采用的阵风使翼型在铅垂方向上受到一个竖直向上的1-cos型阵风w_g, 采用的1-cos型阵风表达式参考文献[14], 具体为

(7)

式中, W₀为阵风的幅值, 本文取0.1倍自由来流速度, 使翼型的有效迎角增加5.71°。H指阵风的离散尺度, 本文参考文献[14], 取12.5倍翼型弦长。具体的1-cos阵风速度型如图 4所示。

2 计算结果分析 2.1 雷诺数对阵风响应影响研究

针对不同雷诺数条件下, 低雷诺数翼型FX63-137的阵风响应特性展开研究。选取了雷诺数分别为1.0×10⁵, 3.0×10⁵, 5.0×10⁵, 1.0×10⁶的4种计算状态。图 5所示即为采用前文所述的阵风模型, 在0°迎角下, 4种计算状态下升力系数增量随无量纲时间的变化图。其中, 纵坐标表述的为升力系数在受到阵风扰动之后的增加量。

从图 5可见, 翼型升力系数阵风响应的整体趋势在采用经验算法ONERA以及CFD算法的各个雷诺数下相一致, 且与所采用的阵风速度型相吻合。但是翼型升力系数随阵风响应的幅值在不同计算状况下显示出了较大的不同。采用ONERA方法计算的升力响应幅值要明显大于CFD在各个雷诺数下的计算结果。而在CFD计算结果中, 随着雷诺数的增大, 其升力系数响应峰值随之增大。且在雷诺数从10万增加到30万的过程中, 这种变化更为明显。在雷诺数从30万增长到50万及100万的过程中, 升力系数响应的变化幅度随之减小。这主要跟在雷诺数Re=1.0×10⁵时, 低雷诺数效应显著有关。在雷诺数在10万量级时, 低雷诺数所造成的气动力非线性现象显著, 气动特性恶化, 升力线斜率也随之减小。而随着雷诺数的增大, 气动特性得到改善, 翼型升力线斜率随之增大。故随着雷诺数的增大, 升力系数响应幅值也在同时增大。而采用ONERA算法所得出的升力系数响应由于无法考虑雷诺数对翼型气动特性的影响, 故其计算的升力系数较为理想, 整个响应幅值要远大于采用CFD计算所得到响应值。

2.2 翼型阵风响应流场特征分析

低雷诺数翼型FX63-137在不同雷诺数下的阵风响应特性表现出了一定的区别。这跟其在不同雷诺数下气动特性不同相关。本文选取了比较具有代表性的雷诺数分别为1.0×10⁵及1.0×10⁶的2种算例分析其流场特征。图 6所示为上述2个雷诺数计算状态下翼型阵风响应周期内不同时刻的流场图。其中K_T表示湍动能, 其表征流动的湍流强度大小。

从图 6和图 7中可以看出, 在不同的雷诺数下, 翼型的流场结构有较大的区别。在雷诺数为1.0×10⁵时, 翼型上下表面存在饱满的层流分离泡结构。在翼型上表面, 层流分离泡的影响范围达到了0.5倍翼型弦长左右, 且分离泡一直延续到了翼型后缘。而当雷诺数增加到1.0×10⁶时, 虽然翼型的上下表面依然存在层流分离泡, 但是可以明显看出, 其分离泡结构很短薄, 只存在于翼型弦长中部很小的范围之内。其翼型上表面分离泡影响区域只占到整个翼型弦长的0.05倍。

而在相同雷诺数下, 在阵风响应周期的不同时刻内, 其流场结构也表现出了一定的区别。在2个不同雷诺数状态下的计算结果均表现出随着阵风周期分离泡位置移动的特点。可以看到, 在雷诺数为1.0×10⁵时, 随着阵风扰动幅值的增加, 翼型有效迎角逐步增加, 翼型上表面的层流分离点逐渐前移, 在阵风扰动强度达到最大时, 其上表面流动分离点也随之前移到最前端。在这之后, 阵风扰动幅值开始逐渐减小直至消失, 其流动分离点也随之后移。翼型下表面的流动分离点随阵风扰动而移动的趋势也与上表面相一致。除此之外, 翼型下表面分离泡的形态也随阵风扰动强度的变化而变化, 可以看到随着扰动强度的增大, 下表面分离泡变得越来越扁平。在雷诺数为1.0×10⁶时, 其分离泡随阵风扰动的变化规律也表现出了相同的变化趋势。但是由于其雷诺数相对较大, 低雷诺数效应较弱。上下表面的分离泡都表现出了短薄的特点, 故其上下表面分离泡的形态随阵风扰动并无明显变化。

随着阵风扰动幅度的变化, 翼型的有效迎角随之改变。在低雷诺数状态下, 层流分离泡的位置及形态随翼型迎角的变化而变化。层流分离泡形态及位置的变化导致了翼型的有效弯度及有效厚度发生改变。导致了翼型升力线斜率变化以及翼型的气动特性非线性。总体而言, 随着计算雷诺数的降低, 翼型的升力线斜率减小, 从一定程度上表明, 低雷诺数下分离泡的存在延缓了翼型气动特性的变化。在雷诺数为1.0×10⁵时, 低雷诺数效应显著, 分离泡饱满。故其分离泡位置的移动及形态的变化对翼型的有效弯度及有效厚度的影响较大, 按照上述分离泡延缓气动特性变化的角度而言, 翼型升力线斜率因此较小。而当雷诺数为1.0×10⁶时, 低雷诺数效应已经相对较弱, 虽然分离泡也随阵风扰动位置及形态在随之改变, 但是其翼型表面的分离泡结构已经相对短薄, 对翼型有效弯度及厚度的影响相对较小, 对气动特性变化的延缓作用相对较小, 翼型升力线斜率相对较大。这解释了为何随着雷诺数的增大, 翼型升力系数阵风响应的增量幅值增大。同时由于分离泡的存在对翼型气动特性的变化有一个延缓作用, 可以从图 7中可以看出, 当扰动消失的1T时刻, 其流动分离点的位置并未回到扰动发生前0T时刻的位置。而是有一个相对滞后。这也与图 6中, 在阵风扰动消失后, 升力系数的增量并没有直接回到0, 而是有一个滞后相吻合。

图 8所示为2个计算雷诺数下, 阵风响应周期内不同时刻的压力系数分布图。由图 8可知, 由于雷诺数不同, 翼型表面压力系数分布也有较大的区别。在雷诺数为1.0×10⁵时, 翼型上表面存在较长的压力平台。并随着阵风扰动幅度的增强, 压力平台前移, 导致翼型后缘高压区前移, 压力系数随有效迎角增加而增加的幅值减小, 气动特性变化得到缓和。而在雷诺数为1.0×10⁶时, 可以看到翼型上表面的压力平台相对较短, 虽然其也随阵风扰动强度的变化在移动, 但是其影响区域以及影响范围较小。故从图 8可以看出其翼型表面压力系数随阵风扰动强度变化的幅值要大于雷诺数为1.0×10⁵时的变化。与前文流场结构所表现出的特性相一致。

2.3 大迎角下不同雷诺数阵风响应特性分析

在大迎角下, 翼型进入失速状态, 后缘会发生较大的分离, 气动特性表现出强非线性。而雷诺数对翼型的失速特性也有着相当大的影响。本文接下来将研究大迎角下不同雷诺数的低雷诺数翼型阵风响应特性。图 9给出了大迎角时各个雷诺数下的翼型升力系数增量的阵风响应示意图。

图 9中虚线所示为达到阵风扰动幅值最大的时刻点。从图 9可知, 迎角为10°时, 在阵风扰动还未达到峰值时, 升力系数增量就开始下降, 这种趋势一致延续至扰动结束, 且在阵风扰动结束时, 升力系数的增量是个负值。在这之后, 流场逐渐恢复, 升力系数增量的值逐渐从负变零。此外, 随着雷诺数的减小, 升力系数增量达到峰值的时刻在逐渐前移, 峰值在逐渐减小。且在阵风扰动消失时刻, 升力系数增量负值的幅值随着雷诺数的减小而增大。图 10所示的为迎角为10°在阵风扰动达到峰值时, 不同雷诺数下的流场结构图。由图 10可知, 在T=0.5时刻, 4个雷诺数下的翼型均表现出前缘有一个小分离泡。在这之后, 翼型后缘处发生了较大范围的分离。且随着雷诺数的减小, 翼型上表面的分离区明显增大。这表明随着雷诺数的减小, 翼型的失速迎角在前移。随着阵风扰动强度的增大, 翼型的有效迎角逐渐增大, 小雷诺数下的翼型率先进入失速, 故升力系数增量率先达到峰值并开始逐渐降低。

从图 9还可知, 升力系数增量曲线并不像速度型一样表现出对称性, 在阵风扰动强度相等时刻, 升力系数的增量有着明显的差异。并在扰动结束时, 升力系数要明显低于扰动施加时刻。且雷诺数越低, 这种现象越明显。本文选取雷诺数为1.0×10⁵的算例进行研究, 图 11~图 13所示的分别为雷诺数为1.0×10⁵时, 在一个完整的阵风周期内不同时刻的流场结构图、翼型表面压力系数分布图以及阵风周期内升力系数增量随翼型有效迎角的变化图。

由图 11可知, 随着阵风扰动强度的增大, 翼型的有效迎角增加, 翼型前缘处的层流分离点逐渐前移, 后缘的分离区逐渐扩大。在阵风扰动强度逐渐减弱时, 前缘处的层流分离点逐渐后移, 后缘处的分离区逐渐减小。但是从图 10可以看出, 在阵风扰动强度相等的0.25T与0.75T时刻, 以及阵风扰动为0的0T以及1T时刻, 其流场结构存在一定的差别。扰动强度上行时刻相比下行时刻, 在阵风强度相等的情况下, 其层流分离点相对靠后, 后缘的分离区范围较小。图 12也很好地反映出了上述现象, 由于0.75T时刻相对0.25T时刻, 翼型上表面层流分离点相对靠前, 导致压力平台前移, 后缘高压区前移, 前缘吸力峰值有较大的下降。而由于后缘分离区范围相对较大, 0.75T时刻相比0.25T时刻, 翼型后缘的压力系数相对更大。在这两者的共同作用下, 扰动强度下行时刻在有效迎角相等的情况下相对于扰动上行时刻, 升力系数有较大的减小。造成了如图 13所示的升力系数迟滞现象。图 13中的箭头表示有效迎角的变化方向。

从图 13可以看到, 在有效迎角相等的情况下, 在迎角上行时刻的升力系数随阵风扰动的增量要明显小于有效迎角下行时刻, 并在扰动消失时, 升力系数增量为负值, 并未回到扰动施加点, 形成一个不封闭的迟滞环。同时, 由于低雷诺数效应, 翼型的气动特性变化相对缓和, 更加导致了流场结构相对于有效迎角的滞后。这也导致了图 12中表现出的迟滞环相对普通翼型运动所形成的迟滞环更加饱满。

3 结论

1) 随着雷诺数的减小, FX63-137翼型表面的层流分离泡影响范围越来越大, 分离泡形态越来越饱满, 气动特性非线性增强, 延缓翼型气动特性的变化, 导致翼型的升力系数曲线的斜率减小, 在阵风扰动下的升力系数增量幅值减小。

2) 在迎角α=10°时, 由于进入翼型失速区, 在阵风扰动幅度还未达到峰值时, 升力系数增量就开始减小。在阵风扰动结束时, 由于流场结构相对有效迎角的滞后, 在翼型前缘吸力峰值随有效迎角降低时, 后缘分离区由于之后还未恢复, 导致扰动结束时, 翼型升力系数增量是个负值。并且随着雷诺数的减小, 翼型的失速迎角前移, 升力系数增量达到峰值的时刻越来越前移, 在扰动消失时, 升力系数增量负值的幅值越来越大。

3) 在雷诺数为1.0×10⁵时, 在阵风扰动强度相等的扰动强度上行时刻与扰动强度下行时刻, 翼型流场结构表现出较大差异, 下行时刻翼型前缘层流分离点相对上行时刻相对靠前, 后缘分离区范围相对较大。这两者导致下行时刻翼型阵风响应的升力系数增量相对上行时刻小, 形成一个不封闭的迟滞环。

低雷诺数条件下, 翼型的阵风响应幅值相对于常规雷诺数有所减小, 气动响应特性相对于阵风变化有一个相对滞后, 这对太阳能飞机这类低雷诺数飞机的结构设计与飞行控制提供了相应的指导。