飞翼布局太阳能无人机低雷诺数反弯翼型设计

王科雷, 周洲, 郭佳豪, 徐德

(西北工业大学 航空学院, 陕西 西安 710072)

太阳能无人机采用太阳光辐射能作为动力,理论上具有“绿色永久飞行”的能力,在侦察监视、通信中继等方面具有极为广阔的应用前景。近年来,随着国内外应用需求不断增长,轻小型低空太阳能无人机正在由新概念探索、关键技术攻关,迅速向实用化、工程化迈进[1-2]。然而,由于飞机尺度小、飞行速度低,其低雷诺数翼型气动设计问题日益严峻。

一般在低雷诺数状态下,翼型附面层内易于发生分离、转捩和再附等复杂流动现象,进而导致翼型气动特性急剧恶化[3]。因此,常规低雷诺数翼型设计主要是以增升减阻为目的,如“在翼面上保持附着层流,使不产生分离流”的层流翼型设计思想[4-5],或“保持翼型前缘流动附着,对翼面流动分离和转捩进行控制和设计”的高升力翼型设计思想[6-10]等,这些在常规布局太阳能无人机设计应用中均已取得较好的效果。然而,面对太阳能无人机对更高气动效率和更高光伏组件铺设率的要求,采用常规布局已经很难再挖掘气动和能源系统的潜力,开展飞翼、翼身融合等非常规布局的低雷诺数翼型设计研究已经成为必要。然而,由于飞翼布局自身气动力设计较为独特,在翼型设计时即需对俯仰力矩特性进行综合考虑,但就目前而言,国内外对于低雷诺数翼型俯仰力矩特性的关注相对较为欠缺,仅甘文彪等[11-12]针对某全翼式布局太阳能无人机低雷诺数反弯翼型进行了设计研究和数值分析。

此外,对于轻小型低空太阳能无人机而言,追求长航时性能是其气动布局研究的最终目标,因此该类型无人机往往需要在航时因子较高的高升力系数状态或低需用功率状态下飞行,甚至在靠近最大升力系数附近,以更低的速度、更小的功率飞行,从而获得最长留空时间[12]。这就意味着,采用飞翼布局形式的轻小型低空太阳能无人机将面临高升力、失速和缓、俯仰力矩特性等多重约束和要求,其低雷诺数翼型设计问题亦将更加复杂和困难。

围绕上述问题,本文基于某手抛式小型太阳能无人机的应用需求,结合优化方法开展了飞翼布局约束下的低雷诺数反弯翼型设计研究,对在低雷诺数条件下获取高升力长航时飞行性能的反弯翼型设计思路和设计方法进行探索和分析。

1 数值计算方法

在低雷诺数条件下,翼型附面层内层流分离泡结构(laminar separation bubble,LSB)的产生对压力分布形态影响十分显著,是决定低雷诺数翼型俯仰力矩特性的关键因素。因此,能够准确预测层流分离泡结构的产生和发展对于低雷诺数翼型设计研究至关重要。目前直接数值模拟(direct numerical simulation,DNS)和大涡模拟方法(large eddy simulation,LES)在流动转捩和层流分离泡预测方面具有明显优势[13],但受计算时间、计算资源等限制,并不适用于工程设计。

因此,为提高设计效率,本文基于商业软件FLUENT,采用耦合k-kL-ω转捩模型[14-15]准定常求解雷诺平均N-S(Reyolds-averaged Navier-Stokes,RANS)方程的计算方法来对低雷诺数翼型转捩流动问题进行高效数值模拟和优化设计。此外,为准确反映设计翼型大迎角状态失速特性得到改善的流动机理,采用耦合k-ω剪切应力输运(shear-stress transport,SST)湍流模型[16]的分离涡模拟(detached-eddy simulation,DES)非定常求解方法[17-19]对其大迎角失速分离流动问题进行高精度数值模拟和分析。数值计算过程中,采用二阶迎风MUSCL插值的Roe格式进行空间离散,采用隐式AF方法进行时间离散和推进。

为验证本文数值模拟方法对低雷诺数转捩流动问题和大迎角失速分离流动问题求解的准确性及可靠性,选取拥有丰富气动力试验数据及压强测量数据的Eppler 387翼型低雷诺数绕流算例[20-21]来进行分析研究。参考NASA兰利研究中心低湍流度风洞试验条件,选取数值模拟状态为:Ma=0.09,迎角α=-2°～16°,Δα=2°,翼型弦长雷诺数Rec=3.0×105。如图1所示,分别采用360×150和640×330的O型网格进行划分,近壁面第一层网格高度分别为y+≈0.5和y+≈0.1,网格高度由壁面向计算远场呈指数型增长,增长率分别设为1.1和1.04。在计算过程中,为提高DES模型方法计算效率,采用在RANS准定常数值模拟方法计算结果收敛的基础上继续迭代的方式进行非定常计算,时间步长取为0.000 4 s,子迭代步数设为20,在10 000个时间步之后对计算结果进行采样平均,共采样500个时间步,最终得到时均化气动力结果进行分析研究。

图1 Eppler 387翼型二维网格示意

图2所示为不同计算迎角下Eppler 387翼型准定常/非定常气动力计算结果与试验数据的对比。图3所示为Eppler 387翼型在α=0°～16°,Δα=4°的时均化计算压力分布与试验结果之间的对比。

图2 Eppler 387翼型气动力特性曲线对比

图3 Eppler 387翼型压力分布对比

可以看出,在α=-2°～8°的升力线性范围内,本文准定常/非定常数值模拟方法对Eppler 387翼型的气动特性及压强分布的预测结果均与试验值吻合良好,而相比较于DES非定常数值模拟方法,RANS准定常数值模拟方法对翼型俯仰力矩特性的预测更加准确,这表明k-kL-ω转捩模型对小迎角状态下翼型附面层内层流分离泡结构的捕捉以及对流动转捩位置的预测均较为精准。

随着计算迎角不断增大,RANS准定常数值模拟方法气动力计算误差逐渐增大,尽管其对翼型气动力随迎角变化趋势的整体预测精度尚可,但由图3a)压力分布对比可以看出,该准定常数值模拟方法对大迎角状态下翼型上表面流动分离的预估不足,对翼型失速特性的反映不够准确。而DES非定常数值模拟方法计算精度则相对较高,尤其对翼型大迎角失速状态下的气动力预测与试验值相对吻合更好。然而需要注意的是,由图3b)压力分布对比可以看出,在计算迎角α=16°状态下,试验得到翼型表面压力分布沿弦向变化平缓均匀,但DES非定常时均化结果则波动较大,且总体上对翼型前半部的分离涡强度预估过强,而对翼型后半部的分离涡强度预估稍有不足,这也导致DES非定常数值模拟计算得到的翼型抬头力矩相对试验值较大。

总的来说, 本文基于k-kL-ω转捩模型的RANS准定常数值模拟方法和 DES 非定常数值模拟方法基本能够满足飞翼布局约束下的低雷诺数反弯翼型设计和大迎角失速特性分析的求解精度需求。

2 优化设计问题及方法

2.1 翼型参数化方法

本文采用在翼型设计中应用极为广泛的CST方法[22-23]对低雷诺数反弯翼型进行参数化建模。CST参数化的翼型可以表示为

式中:x,y分别为翼型横坐标和纵坐标;下标u,d分别为上(upper)下(down)的首字母缩写;下标t为翼型后缘(trailing edge,TE)的首字母缩写。C(x)为类函数,取N1=0.5,N2=1.0。S(x)为型函数,其中Aui,Adi为上下翼面待定系数,即翼型的控制变量。Si(x)为Bernstein多项式,其中N取为6,设计翼型上下表面共计14个控制变量。

2.2 优化设计问题描述

针对飞翼布局太阳能无人机这一研究对象,考虑其手抛式起飞方式及追求高升力长航时飞行性能的应用需求,对低雷诺数反弯翼型提出了以下设计要求:一方面,需满足巡航迎角(αcruise)下的高升阻比、高航时因子要求,另一方面,需满足大迎角(αlarge)下的高升力、高航时因子需求。同时,在不同计算迎角下,优化翼型俯仰力矩系数(cm-optimal)需始终不小于基准翼型俯仰力矩系数(cm-baseline)。此外,结合实际工程中结构强度需求及载荷空间需求,将翼型最大相对厚度(t/cmax)及后缘厚度(tTE)均约束为保持不变。

至此,本文低雷诺数反弯翼型优化设计问题可以表达为

(7)

式中:cl与cd为低雷诺数反弯翼型计算升力系数和阻力系数;ω1,ω2为巡航迎角下升阻比和航时因子的权重系数;ω3,ω4为大迎角下升力和航时因子的权重系数;ω5,ω6为巡航迎角特性和大迎角特性的权重系数。

2.3 基于代理模型的优化设计框架

如图4所示,采用多岛遗传算法(multi-island genetic algorithm,MIGA)[24]为搜索器进行单目标寻优,其中子群规模设定为12,岛屿数为12,遗传代数设定为20,交叉率设定为0.7,变异率设定为0.3,岛间迁移率设定为0.5,迁移间隔代数设定为4。为了进一步提高设计效率,在优化过程中,采用均匀试验设计方法[25]与最小化代理模型预测(minimize surrogate prediction,MSP)法则[26]相配合来构建和更新Kriging代理模型[27-28]的方式,实现对实际数值计算过程的逼近和代替。

图4 低雷诺数反弯翼型优化设计流程

具体设计步骤可以概括为:①根据试验设计方法进行初始种群抽样,将初始样本点数目设定为120,采用RANS准定常数值模拟方法并行计算出各样本点响应值;②基于样本点数据建立代理模型,并进行子优化求解;③根据加点法则获取有效的新样本点,并将新的样本点加入构建代理模型的样本点,完成代理模型更新;④判断优化是否结束,未结束则跳转至步骤②,结束则输出优化结果,进入优化翼型分析验证。

3 设计结果及分析

结合某手抛式小型飞翼布局太阳能无人机实际工作状态,将本文低雷诺数反弯翼型设计状态确定为:飞行高度H=6 000 m,来流速度V∞=15 m/s,弦长c=0.4 m,弦长雷诺数Rec=2.5×105,巡航迎角αcruise=8°,计算大迎角αlarge=14°。选取NACA 8-H-12反弯翼型作为基准,在兼顾巡航状态翼型升阻特性的基础上重点关注大迎角失速状态分离流动特性,取ω1=0.4,ω2=0.6,ω3=0.4,ω4=0.6,ω5=0.3,ω6=0.7的权重系数组合作为算例进行优化设计研究。

图5为优化前后翼型轮廓曲线对比,可以看出,优化后翼型前缘半径相对有所增大,翼型头部弯度亦有所增大,而翼型中部正弯度及后缘反弯趋势则相对有所减弱,最大厚度位置亦稍有前移。

图5 优化设计前后翼型轮廓曲线对比

分别采用准定常/非定常数值模拟方法对设计状态下基准翼型和优化翼型的低雷诺数绕流问题进行计算和分析,其中翼型计算网格划分方式及数值计算设置等均与前文低雷诺数翼型算例保持一致,计算迎角范围取为α=-4°～20°,但准定常计算迎角间隔取为Δα=2°,非定常计算迎角间隔则取为Δα=4°。首先,针对采用2种数值模拟方法计算得到的优化前后翼型气动力特性进行对比分析,其次,以RANS准定常数值模拟方法结果为主,对优化前后翼型时均化绕流流场特性变化趋势进行分析,最后,以DES非定常数值模拟方法结果为主,对大迎角失速状态下优化前后翼型非定常流动发展进行研究。

3.1 气动力特性分析

图6为优化前后翼型气动力特性曲线对比。由图可以看出,准定常/非定常2种数值模拟方法在α=-4°～8°迎角范围内计算得到的优化前后翼型气动特性变化趋势一致,而在α=8°～20°迎角范围内计算结果差异较大,但整体上均表明优化翼型在巡航状态及大迎角失速状态下的升阻特性和俯仰力矩特性均相比基准翼型得到明显改善。

图6 优化设计前后翼型气动力特性曲线对比

准定常数值计算结果表明:优化翼型在满足俯仰力矩约束的前提下,全计算迎角范围内升力系数与基准翼型相比较均有所提高,升力线斜率则基本不变,最大升力系数相对提高近12.6%,且最大升力系数对应迎角则由14°增大到16°,可以有效改善小型飞翼布局太阳能无人机的手抛式起飞性能,此外,优化翼型阻力曲线相对基准翼型整体右移,在正迎角状态下阻力特性得到明显改善,同时在较高升力系数下,阻力系数随升力系数的增长相对较为缓慢,这使得优化翼型在设计迎角状态下的升阻比及航时因子均有所提高,其中,在巡航迎角α=8°状态下优化翼型计算升阻比相对提高约3.26%,航时因子相对提高约3.23%,而在计算大迎角α=14°状态下优化翼型计算升阻比相对提高约44.06%,航时因子相对提高约42.70%。

非定常数值计算结果表明:在全计算迎角范围内,优化翼型相比基准翼型的升阻特性始终更优,在巡航迎角α=8°状态下,优化翼型计算升阻比相对提高约6.79%,航时因子相对提高约4.19%。尽管优化前后翼型均在α=12°迎角状态下开始失速,但优化翼型失速特性相对和缓,升力随迎角维持较好,且阻力随迎角增加相对缓慢,如在α=12°,α=16°,α=20°迎角失速状态下,优化翼型计算升阻比相对分别提高约39.77%,55.57%,55.62%,航时因子相对分别提高约44.45%,58.70%,78.21%,这有利于手抛式小型飞翼布局太阳能无人机的高升力长航时飞行。

3.2 时均化绕流流场特性分析

图7给出采用准定常数值模拟方法计算得到的α=8°,12°,16°迎角下优化前后翼型压力分布对比。

图7 优化设计前后翼型压力分布对比

由图可知,在计算迎角α=8°状态下,基准翼型前缘呈陡峭型压力分布,且在计算迎角增大为α=12°时,基准翼型上表面代表层流分离泡结构的低压平台区域迅速前移,且低压平台长度显著减小,基准翼型附面层内流动特征分布发生显著改变,之后随着计算迎角进一步增大,基准翼型头部压力分布特征基本保持不变;相应地,在计算迎角α=8°状态下,优化翼型前缘为圆顶型压力分布,气流加速区内压力变化和缓,且前缘吸力峰值相对基准翼型较低,同时其上表面低压平台相对靠前,低压平台长度相对较小,这在一定程度上抵消了由后缘反弯减弱带来的低头力矩增加,之后随着迎角不断增大,优化翼型上表面低压平台逐渐前移,但前移幅度始终相对较小,压力分布特征变化相对更加缓和,未出现如基准翼型一般的突跃式改变,同时优化翼型前缘压力分布形态逐渐由圆顶型转变为陡峭型。

图8～9分别给出采用准定常数值模拟方法计算得到的优化前后翼型在α=6°～16°迎角范围内的时均化绕流流场特性分布对比,主要包括湍流强度分布和流线分布, 同时也给出翼型上表面层流分离泡结构和二次湍流分离等流场细节的局部放大图,其中标记“S”、“T”、“R”、“ST”分别表示层流分离、流动转捩、湍流再附着、湍流分离。表1～2分别给出各计算迎角下优化前后翼型上表面附面层内流动特性参数,包括图8～9所示层流分离、流动转捩、湍流再附着、湍流分离的弦向相对位置,以及所形成的层流分离泡长度,其中翼型附面层内流动分离及再附着位置主要依靠Vx矢量方向进行判断,流动转捩位置则通过湍流强度It≥0.1判断。可以看出,随着计算迎角变化,翼型附面层内流动呈现2种典型结构:一种是低雷诺数典型层流分离泡结构;另一种是在翼型前缘形成层流分离泡的同时在翼型后部发生湍流分离。

图8 基准翼型绕流流场结构示意

图9 优化翼型时均化绕流流场结构示意

表1 基准翼型附面层内流动特征位置分布

表2 优化翼型附面层内流动特征位置分布

不同的是,对于基准翼型而言,在计算迎角α由8°增大为10°时,翼型附面层内层流分离泡结构显著前移,且层流分离泡长度明显减小,翼型绕流流场结构特征已然发生改变;当计算迎角α由10°增大为12°时,翼型前部层流分离泡位置及长度改变较小,但翼型中后部在后缘反弯轮廓特征的作用下,出现后缘湍流分离现象,翼型绕流流场结构特征又一次发生改变;之后随着计算迎角进一步增大,翼型前部附近层流分离泡位置及长度基本保持不变,翼型中部湍流分离点则逐渐前移,后缘湍流分离区域亦不断增加。

而对于优化翼型而言,其后缘反弯轮廓特征相对较弱,因此在计算迎角α由8°增大为12°的过程中,翼型附面层内流动特征位置基本没有显著改变,流场结构特征始终较为一致;之后随着计算迎角进一步增大,翼型中后部出现湍流分离现象,翼型绕流流场结构特征发生本质改变,但相比于基准翼型,优化翼型附面层内湍流分离位置始终靠后,后缘分离强度亦相对较弱,因此优化翼型在大迎角状态下升阻特性相对更优。

总的来说,在小迎角状态下,基准翼型前缘压力分布呈陡峭型,而优化翼型前缘压力分布则呈圆顶型,这使得优化翼型前部载荷分布变化相对更加缓和,有利于推迟流动分离点的前移,进而拓展升力线性段范围。此外,在大迎角状态下,优化翼型后缘反弯减弱则有利于抑制翼型后部湍流分离的发生,以及降低翼型后部分离流动的强度,是优化翼型失速迎角和最大升力系数均相对基准翼型有所增大的主要原因。

3.3 大迎角失速状态非定常流动特征分析

图10～12分别给出优化设计前后翼型在失速发生前后α=12°,α=16°,α=20°迎角状态下不同时刻(t1=44.20 s,t2=44.28 s,t3=44.36 s,t4=44.44 s,t5=44.52 s,t6=44.60 s)的瞬时压力分布和时均(time averaged)压力分布对比。可以看出,在计算迎角α=12°状态下,优化翼型上表面时均化的低压区域范围(0≤x/c≤0.5)相比基准翼型(0≤x/c≤0.25)明显更大,同时优化翼型后缘附近的压力恢复特征相对基准翼型更优,压差阻力特性相对更佳。此外,基准翼型上表面瞬时压力分布随时间变化规律性相对较差,除翼型前缘吸力峰之外可观测存在3～4处较为明显的低压峰,而优化翼型上表面瞬时压力分布规律性相对较强,除翼型前缘吸力峰之外可观测存在2处明显低压峰,同时在翼型后半段0.6≤x/c≤1.0范围内,优化翼型上表面分布的低压峰值相对基准翼型明显较低,这表明优化翼型瞬态流场内的分离涡流动强度相对基准翼型较弱。

图10 优化设计前后翼型压力 图11 优化设计前后翼型压力 图12 优化设计前后翼型压力分布对比(α=12°) 分布对比(α=16°) 分布对比(α=20°)

当计算迎角增大为α=16°时,基准翼型和优化翼型的瞬时压力分布特征与时均化压力分布特征在宏观上均与α=12°状态较为一致。随着迎角增大,2种翼型前缘吸力峰值均有所增大,而上表面时均化低压区域范围均相对减小,但优化翼型上表面时均化低压区域范围(0≤x/c≤0.35)相比基准翼型(0≤x/c≤0.2)始终较大,这与准定常数值模拟方法计算得到的优化前后翼型压力分布特征变化趋势较为一致(见图7)。

当计算迎角进一步增大为α=20°时,基准翼型和优化翼型的压力分布特征相比α=12°,α=16°状态发生明显改变。对于基准翼型而言,随着迎角增大,其前缘吸力峰值相对降低,同时其上表面时均化压力分布出现范围较广(0.3≤x/c≤0.7)、幅值稳定的低压平台,也即有大范围翼面始终处于分离涡流动的影响之中。而对于优化翼型而言,随着迎角增大,其前缘吸力峰值不断增大,且相比较基准翼型显著较大,这也是优化翼型在当前迎角下的计算升力相比基准翼型更大(见图6a))的主要原因。同时优化翼型上表面瞬时压力分布的波峰、波谷基本上仅出现在0≤x/c≤0.5范围内,这与基准翼型的分离流动特征存在显著差异,表明优化翼型的分离涡流动强度相对基准翼型明显较弱。

为进一步分析计算迎角α=20°状态下翼型绕流流场内分离涡结构的发展,如图13所示,采用Q准则[29]对α=20°迎角下各时刻流场内的涡核进行判断。可以看出,基准翼型在大迎角失速状态下的分离流动是以前缘失速涡的产生、向下游发展和向翼型斜上方脱落为主导,其上翼面中段会伴随产生有一个由剪切层分离形成的分离泡,这与图12a)中0.3≤x/c≤0.7弦长范围内时均化低压平台相对应。相比较下,优化翼型在大迎角失速状态下,其后缘会形成逆时针旋转的尾涡(起动涡),同时其前缘也会形成前缘失速涡,而最显著的特点是,前缘失速涡在向下游发展过程中不脱落,始终滞留在翼型近壁面附近,而顺时针旋转的失速涡在发展至翼型后缘处时会与逆时针旋转的尾涡配对,二者相互作用使两涡心都存在向下的速度,该速度与主流速度叠加后,尾涡向翼型下游泄出,而失速涡则在翼面和新生尾涡的作用下逐渐消弭,显然,这种特殊而复杂的分离涡形态有效抑制了优化翼型近壁面区域内大尺度分离涡的形成和发展。

图13 优化设计前后翼型绕流流场内分离涡结构随时间发展(α=20°)

4 结 论

1) 本文所发展的基于k-kL-ω转捩模型的RANS准定常数值模拟方法和耦合k-ωSST湍流模型的DES非定常数值模拟方法对于低雷诺数转捩流动问题和翼型大迎角失速分离流动问题具有较高的求解精度,其中RANS准定常数值模拟方法对于翼型附面层内流动分离及转捩位置预测较为精准,而DES非定常数值模拟方法则可以较好地捕捉大迎角状态下翼型绕流流场内分离涡的产生和发展,可以满足低雷诺数反弯翼型设计和失速特性研究需求。

2) 本文基于手抛式飞翼布局太阳能无人机的应用需求所发展的低雷诺数反弯翼型优化设计方法切实可行。优化翼型最大升力系数以及大迎角失速状态下的升阻比、航时因子均相对显著提高,同时优化翼型绕流流场内相对基准翼型较多较弱的分离涡结构有效抑制了距离翼型较近区域内的大尺度分离涡的形成和发展,避免了失速状态升力损失,达到了失速和缓的设计目标,有利于小型飞翼布局太阳能无人机的手抛式起飞和高升力长航时飞行。

3) RANS准定常数值模拟结果表明:在巡航状态下,翼型前缘圆顶型压力分布形态使得其前部载荷分布变化更加缓和,有利于推迟流动分离点的前移,进而拓展升力线性段范围;而在大迎角状态下,翼型后缘反弯减弱有利于抑制翼型后部湍流分离的发生,以及降低翼型后部分离流动的强度,可以有效改善翼型大迎角失速特性,提高失速迎角和最大升力系数。

4) DES非定常数值模拟结果表明:在计算迎角α=12°～16°范围内,优化翼型上表面时均化的低压区域范围相比基准翼型更大,同时优化翼型后缘区域压力恢复形态相对更优;而在计算迎角α增大为20°时,基准翼型前缘吸力峰值显著降低,且上表面存在大范围低压平台,而优化翼型前缘吸力峰值有所增大,且翼型上表面瞬时压力分布的波峰、波谷出现的区域保持稳定。这是因为基准翼型在大迎角失速状态下的分离流动是以前缘失速涡的产生、向下游发展和向翼型斜上方脱落为主导,而优化翼型在大迎角失速状态下的绕流流场结构则较为特殊,通过失速涡和后缘尾涡的相互作用,有效避免了近壁面区域内大尺度分离涡的形成和发展。