耦合多螺旋桨滑流影响的低雷诺数机翼设计

王科雷， 周洲,*， 祝小平

1.西北工业大学 航空学院， 西安 710072 2.西北工业大学 无人机特种技术重点实验室， 西安 710065

耦合多螺旋桨滑流影响的低雷诺数机翼设计

王科雷1,2， 周洲1,2,*， 祝小平2

1.西北工业大学 航空学院， 西安 710072 2.西北工业大学 无人机特种技术重点实验室， 西安 710065

以某型手抛式太阳能无人机(UAV)模型为对象进行考虑多螺旋桨滑流影响的低雷诺数机翼平面形状设计研究。首先，基于升力面理论发展了准定常求解多螺旋桨/机翼相互气动干扰问题的涡格法(VLM)程序，并采用建立参考翼型气动特性数据库的形式发展了相关低雷诺数修正(LRC)方法；然后，通过对翼型、低雷诺数机翼及单螺旋桨/机翼算例的数值模拟及与相关实验结果的对比，验证了本文数值方法具备模拟低雷诺数复杂流动问题的可靠性及准确性；最后，对某型手抛式太阳能无人机简化拉力多螺旋桨/机翼模型进行了直接优化设计及反设计，并通过具有较高精度的CFD准定常求解技术对优化结果进行了验证。结果表明：以CFD方法计算结果为参考，本文涡格法程序及低雷诺数修正方法能够准确高效地计算相关低雷诺数复杂流动问题；传统未考虑多螺旋桨滑流影响的设计机翼在实际螺旋桨工作状态下将偏离设计点，机翼气动特性得不到提高；考虑螺旋桨滑流影响的优化设计方法能够有效改善机翼阻力特性，相对应地，在设计状态下优化机翼总阻力能够降低19.52 counts。

手抛式太阳能无人机; 低雷诺数机翼; 多螺旋桨/机翼相互气动干扰; 涡格法; 低雷诺数修正; 优化设计; 反设计

随着近年来无人机(UAV)技术的快速发展，世界各国都在积极拓展民用无人机的应用范围。无人机自身并没有固定的标签，它可以如同“互联网+”一般被广泛地应用到几乎所有行业中，展现出极大的商业价值和应用前景。而中国地域广阔，存在自然灾害频繁发生且种类繁多、边境监控困难、海洋安全管控较弱等诸多问题，亟需依据不同领域或行业的具体需求从各个方面积极提高自身无人机的研发能力。

从边境巡逻的角度出发，在无人机平台安全、可靠的基础上还需要具备2项能力：① 简单、易操作的起降模式；② 足够大的航程和航时。此时手抛式固定翼太阳能无人机[1]就显示出其他种类无人机无法比拟的巨大优势。然而，由于此类无人机尺度较小，且飞行速度较低，其低雷诺数复杂流动特征[2-4]极为显著，数值模拟难度及气动设计难度均相对常规飞行器较大。

作为典型的螺旋桨类飞行器，将螺旋桨滑流影响纳入无人机气动设计阶段将显得尤为重要。自1986年Kroo[5]首次提出在螺旋桨旋转的真实状态下进行机翼设计相比干净机翼设计更具意义以来，Veldhuis[6-8]及Rakshith[9]等相继对螺旋桨/机翼构型各种布局参数进行了气动影响数值研究及气动布局设计研究；乔宇航等[10]基于升力线理论，通过对机翼几何扭转分布的控制达到了降低机翼诱导阻力的目的；徐家宽等[11]采用基于雷诺平均Navier-Stokes(Reynolds Averaged Navier-Stokes, RANS)方程的多重参考坐标系(Multi Reference Frame, MRF)方法[12-13]对某型运输机机翼局部区域剖面翼型进行了考虑螺旋桨滑流影响的气动优化设计，并且采用Kriging代理模型[14-15]进一步减少计算量，提高设计效率。

尽管国内外众多学者已经针对螺旋桨/机体一体化设计问题进行了大量研究，但其气动设计过程主要是在常规雷诺数状态下进行，而且受到计算效率的限制，其所采用的数值计算方法大多是以无黏近似求解方法为主，这对于太阳能无人机气动设计所处的低雷诺数强黏性环境显然行不通。另外，就目前而言，低雷诺数流动本身的数值计算便已经较为困难，而当引入多个螺旋桨滑流气动干扰后其数值计算更将是难上加难，因此常规的太阳能无人机设计过程[16-18]仍未将动力系统的气动影响考虑在内。总的来说，低雷诺数流动数值模拟方法精度和效率之间的矛盾是制约太阳能无人机螺旋桨/机翼一体化气动设计发展的主要因素。

基于以上考虑，本文首先发展了一套准定常求解多(拉力)螺旋桨/机翼相互气动干扰问题的涡格法(Vortex Lattice Method, VLM)程序，并以建立参考翼型气动特性数据库的形式发展了相关低雷诺数修正(Low Reynolds Correction, LRC)方法；然后以CFD方法作为参考，通过对低雷诺数翼型及机翼算例的数值模拟开展了数值方法的验证及评价；最后基于某型手抛式太阳能无人机模型开展了低雷诺数机翼平面形状优化设计研究，并基于CFD方法对设计结果进行了验证及分析。

需要说明的是，本文涡格法及相关低雷诺数修正方法程序均采用FORTRAN语言编写，气动优化设计方面采用iSIGHT商业软件实现，而CFD计算方面则采用Fluent商业软件实现。

1 数值方法

1.1 涡格法

涡格法是基于薄物体的假定简化，在流体无旋、无黏、不可压缩的无界流域内，使用离散的源汇、涡和偶极子等奇点组合来模拟流动现象，其数值计算较为简便，应用极为广泛。

面对多个螺旋桨加机翼的布局形式，采用马蹄涡和源的奇点组合，共5个涡系形成离散数学模型：螺旋桨桨叶涡系、螺旋桨尾迹涡系、机翼翼面涡系、机翼翼面源系及机翼尾迹涡系。其中，随着机翼模型的确定，机翼翼面涡系及尾迹涡系中涡格节点及控制点的空间位置始终保持不变，而各螺旋桨桨叶及尾迹涡系涡格节点及控制点的瞬时空间位置则由输入文件中螺旋桨安装位置、螺旋桨转速及计算时刻所确定。

涡格法离散数学模型中，机翼翼面区域内展向划分50格，并在机翼翼尖和螺旋桨下游区域进行局部加密，而弦向则按余弦函数划分40格，共计2 000个涡格；螺旋桨模型仅包含桨叶而忽略轮毂，各桨叶叶面区域展向均匀划分25格，弦向均匀划分40格，共计1 000个涡格；整个计算模型总涡格量为(1 000×nblade×npro+2 000)个，nblade为单个螺旋桨的桨叶数目，npro为计算模型所采用的螺旋桨数目。

另外，考虑到涡格法旨在进行气动特性快速计算，在其数学建模过程中主要进行了2项简化：① Michael和Brian[19]将机翼和螺旋桨的自由尾迹形式与预定尾迹形式进行了对比，计算表明基于2种尾迹形式的升力、阻力计算结果误差不超过0.1%，因此本文涡格法程序中采用预定尾迹的形式以提高计算效率。其中，机翼尾迹涡系采用10倍弦长的考虑机翼展向两端尾涡卷曲的二次抛物线型固定尾迹涡模型，而螺旋桨尾迹涡系则近似采用5倍弦长的由桨叶后缘各面元节点带出的螺旋面型预定尾迹涡模型。② Rangwalla和Wilson[20]提出，为了符合实际流动问题中机翼表面不可穿透的条件，近似地将与机翼实体相交的螺旋桨尾迹面元进行“消除”(Wake Snipping)处理，同时该处理方式在对复杂外形的非定常势流求解迭代过程中能够有效提高数值计算的效率和稳定性，因此本文涡格法程序中也采用同样的对螺旋桨尾迹面元的“消除”处理。

螺旋桨与机翼之间的相互气动干扰通过各不同涡系之间的速度相互诱导来实现。将螺旋桨旋转周期划分为数个瞬时状态，当某一时刻螺旋桨瞬时相位角即桨叶位置确定时，各涡系之间诱导速度影响系数矩阵也可以通过Biot-Savart法则确定，基于物面不可穿透条件解得当前时刻各涡系内各面元涡强度，进而解得瞬时气动力系数，最终以一个周期内各瞬时气动力系数之和求平均的方法获得准定常气动特性参数。

1.2 低雷诺数修正方法

目前有关低雷诺数复杂流动数值计算及修正的方法发展仍相对不够成熟，大部分研究[21-22]所发展的低雷诺数修正方法主要是依据经验公式对数值结果直接进行修正，缺乏物理内涵，这对简单线性问题可能还较为适用，但当应用于多螺旋桨/机翼相互气动干扰的低雷诺数复杂流动问题时其修正精度及可靠性均无法得到保证。

因此，本文由物理层面出发，对低雷诺数平直机翼复杂绕流问题作出2项简化：① 将低雷诺数机翼绕流问题简化为机翼无黏绕流问题与低雷诺数黏性影响的叠加，也即在涡格法数值计算基础上假设低雷诺数黏性影响是直接作用于各涡格内的涡线段强度；② 借鉴叶素理论[23]基本原理及思想，将机翼型阻的估算问题简化为各展向位置当地迎角下的低雷诺数翼型阻力积分问题。基于以上2项简化，本文发展的低雷诺数修正方法具体步骤可以概括如下：

步骤1建立当前低雷诺数条件下的翼型气动特性参考数据库。分别采用无黏求解方法(求解欧拉方程)和黏性求解方法(求解RANS方程)计算出当前雷诺数条件不同迎角下的机翼剖面翼型气动力参数，将2种情况下翼型计算升力的比值作为低雷诺数条件对翼型环量的影响系数λ，建立当前雷诺数Re条件不同计算迎角α下的黏性影响系数Re-α-λ数据库。另外，基于黏性计算结果建立当前雷诺数条件下的翼型计算升阻力系数Re-cl-cd数据库。

步骤2在涡格法计算过程中，基于机翼不同展向位置的当地来流迎角，参考Re-α-λ数据库选取对应黏性影响系数λ对该区域所有弦向涡格中涡线段强度进行修正，之后继续求解以得出修正后的升力系数CL及诱导阻力系数CDi。另外，考虑到：① 机翼翼尖附近区域气流下洗作用强度极大，造成了该区域翼段所产生升力在机翼整体升力中仅为一小量，对其进行修正与否对机翼整体气动特性计算结果影响不大；② 尽管机翼翼尖仍处于低雷诺数流动中，但其绕流流动特征是以翼尖涡形式为主，且流动形态极为复杂，难以直接建立相关数据库以进行低雷诺数修正。因此，在低雷诺数修正程序编写中对机翼翼尖附近区域内涡格的涡线段强度将不作任何修正处理，且将该区域范围初步定义为由翼尖向内0.5倍弦长的长方形区域。

步骤3基于沿机翼展向分布的升力系数计算值，参考Re-cl-cd数据库对应求出不同展向位置的各剖面翼型阻力系数，则机翼总型阻系数CDp可以按式(1)求解：

(1)

式中：b为机翼展长；c为弦长；y为机翼展向位置；Sw为机翼总面积。最终机翼总阻力系数CD即为机翼诱导阻力系数与机翼型阻系数之和。

另一方面，由于螺旋桨沿径向的扭转角、桨叶宽度以及桨叶剖面翼型变化较大，流动形态极为复杂，且随着螺旋桨的尺寸及转速等改变时其各剖面的当地迎角、当地雷诺数等亦改变较大，若仍采用上述建立参考数据库的方法来进行修正将显得不切实际。考虑到本文主要目的是将电推进系统(螺旋桨)滑流诱导的气动影响纳入机翼优化设计过程中，而并不着重强调对螺旋桨气动及流动特性的准确计算，因此，在低雷诺数修正程序编写中对螺旋桨模型中弧面上涡格的涡线段强度亦将不作任何修正处理。

1.3 CFD方法

采用基于结构-非结构混合网格技术耦合转捩模型准定常求解RANS方程的MRF方法对多螺旋桨/机翼气动干扰问题进行数值模拟。空间离散采用二阶迎风MUSCL(Monotone Upstream-centered Scheme for Conservation Laws)插值的Roe格式，时间离散与推进则采用隐式AF(Approximate Factorization)方法。

1) MRF模型方法

MRF模型方法是一种对螺旋桨滑流进行准定常数值模拟的数学方法，相比于过分耗费计算资源的非定常求解方法，MRF模型方法在更加节省计算资源的同时仍能获得较高的数值模拟精度，在定轴旋转体的气动计算中应用较为广泛[12-13]。

MRF模型方法主要思想是：通过在各螺旋桨周围建立一个规则封闭圆柱流动区域来模拟螺旋桨的旋转运动。建立与螺旋桨具有相同旋转运动方式的旋转坐标系，通过相应的数学转换以及旋转区域与非旋转区域的数据插值传递，实现在静态网格下的包含旋转气流的流场数值模拟。其旋转坐标系下积分形式控制方程为

(2)

式中：t、V和A分别为时间项、流体微元控制体和围成控制体的闭曲面；n为曲面或曲线的单位法向量；Q、H、HV和G分别为守恒变量、无黏通量项、黏性通量项和坐标转换的添加源项，其具体求解公式可参考文献[12]。

2) 结构-非结构混合网格技术

与MRF模型方法的远场静止流动区域及圆柱旋转流动区域相对应，将计算网格划分为静止域网格和运动域网格。针对静止区域划分结构化网格以减小网格总量，节约计算时间；针对运动区域，由于螺旋桨桨叶在径向位置具有不同的叶素安装角，桨叶高度扭转，几何外形比较复杂，划分非结构化网格可以在保证计算精度的同时降低桨叶的网格难度，提高网格生成效率。

图1所示为螺旋桨/机翼构型混合网格结构。机翼所处的非旋转区域结构化网格近壁面网格y+约为0.5，网格量始终保持在600万左右。而螺旋桨所处圆柱形旋转区域非结构化网格近壁面网格y+约为1.5，网格量始终保持为200万左右。

图1 螺旋桨/机翼构型混合网格结构示意图Fig.1 Schematic of hybrid grids structure of propeller/wing configuration

3)k-kL-ω转捩模型

通常，巡航状态下的太阳能无人机绕流问题是处于104～105量级的低雷诺数范畴，在数值计算中必须考虑层流转捩的问题。因此，采用Walters等[24-25]以及Volino[26]提出的基于局部变量“层流动能”的k-kL-ω转捩模型进行求解，其输运方程组可写为

(3)

(4)

(5)

湍流动能及层流动能生成项及近壁面耗散项表达式分别为

(6)

(7)

式中：x为坐标轴系，其下标i、j表示各轴系方向；k为动能；ν为黏性系数；下标T和L分别表示湍流和层流；下标s和l分别表示小尺度和大尺度；ω为湍流频率；αT为湍流标量扩散率；S为张力率梯度；R和RNAT分别为由旁路转捩和自然转捩引起的湍流产生项；其相关变量及相关参数具体值可参考文献[24]。

2 数值方法验证及评价

本文基于涡格法发展低雷诺数修正方法(统称为VLM-LRC方法)的主要目的是为了在低雷诺数飞行器数值优化过程中取代计算效率相对较低的CFD方法，在提供足够计算精度的同时进一步提高设计效率。考虑到本文所发展低雷诺数修正方法的关键环节为相关参考数据库的建立，因此该参考数据库的准确性及可靠性就显得尤为重要，而这些均依赖于所采用数值模拟方法求解低雷诺数问题的计算能力。因此，本节将首先对CFD方法数值模拟低雷诺数翼型绕流问题的计算精度进行验证，之后以CFD方法计算结果为参考来对VLM-LRC方法计算精度及效率进行验证和评价。主要包含以下3部分内容：

1) 通过低雷诺数条件下翼型算例验证CFD方法具备精确描述低雷诺数复杂流动特性的能力。

2) 通过低雷诺数机翼及单螺旋桨/机翼算例验证VLM-LRC方法在低雷诺数复杂流动问题求解方面具备较高的计算精度。

3) 通过对CFD方法与VLM-LRC方法数值计算耗时的对比，证明VLM-LRC方法在获得足够计算精度的前提下能够有效提高计算效率。

2.1 CFD方法验证

[27]中相关实验条件及结果，对NACA633-018翼型进行数值模拟。选取计算状态为：远场来流速度V∞=15 m/s，来流湍流度Tu∞=0.25%，弦长雷诺数Rec=3.0×105。图2给出了不同迎角下翼型计算压力系数Cp分布与实验结果对比。

图2 翼型压力系数分布对比 Fig.2 Comparison of airfoil pressure coefficient distribution

由图2可以看出，本文CFD方法在不同迎角下翼型压力分布计算结果与实验值吻合很好，k-kL-ω转捩模型对翼型表面压力平台区域预测精度较高，也即对翼型绕流流场中层流分离、转捩及湍流再附位置的预测精度较高。这也验证了本文CFD方法精确描述低雷诺数复杂流动的有效性，为下文计算分析提供可靠的气动特性计算基础。

2.2 VLM-LRC方法验证

参考文献[28]中相关实验条件及结果，对展弦比为8的FX 63-137低雷诺数平直机翼模型(Wing)进行数值模拟。选取计算状态为：远场来流速度V∞=30 m/s，飞行高度H=20 km，来流湍流度Tu∞=0.1%，弦长雷诺数Rec=3.0×105。另外，在该低雷诺数机翼的基础上再构造单螺旋桨/机翼模型(Spro)进行数值模拟，其中螺旋桨模型采用直径为1.2 m的某型双叶螺旋桨X1，以0° 安装角安装在机翼中心正前方0.8 m处，旋转方向为顺气流逆时针方向，螺旋桨转速为2 500 r/min。

图3为分别采用CFD方法及VLM-LRC方法计算得到的各构型机翼气动力数值计算结果与实验结果对比。图4为2° 迎角下2种构型机翼展向升力系数分布对比。其中，整个机翼沿展向平均划分为40份，而黑色箭头方向表示螺旋桨诱导下气流上洗、下洗特征。由于本文更加侧重于考虑机翼不同展向位置各等展长翼段升力对机翼整体升力的贡献量，因此各翼段升力系数将始终采用机翼总面积作为参考面积进行计算，而下文机翼沿展向分布的升力计算方式亦与此处相同。

图3 机翼气动特性数值计算与实验结果对比Fig.3 Comparison of wing aerodynamic performance between numerical and experimental results

图4 2° 迎角下机翼展向升力系数分布对比 Fig.4 Comparison of wing span-wise lift coefficient distribution (α=2°)

由图3可以看出，采用CFD方法计算得到的干净机翼气动力系数与实验数据符合良好，仅在大迎角条件下相对误差较大；采用VLM-LRC方法计算得到的2种构型气动力系数均与CFD方法计算结果吻合较好，但VLM-LRC方法计算升、阻力曲线随迎角的增长率相对较小，全迎角范围内2种数值方法计算结果之间相对误差均不超过6.5%；当引入螺旋桨滑流影响时，机翼升、阻力均有所增大，而VLM-LRC方法及CFD方法计算结果基本一致，所得到的气动力系数增加量亦基本相同。

由图4可以看出，VLM-LRC方法计算得到的2种构型机翼展向升力分布曲线均与CFD方法计算结果吻合良好，仅在螺旋桨中心轮毂下游区域以及翼尖附近区域的计算升力差别较大，这可能与涡格法对螺旋桨几何模型的简化以及低雷诺数修正方法对螺旋桨及机翼涡格强度的修正处理等均有关系，但整体上并不影响VLM-LRC方法计算结果与CFD方法计算结果之间的一致性。另外，当前状态下螺旋桨对远场来流的加速效应显著强于旋转效应，螺旋桨下游机翼翼段当地雷诺数相对增大，机翼剖面升力系数均显著增高，同时由于旋转效应始终存在，机翼翼段上洗侧剖面翼型升力相对增量大于下洗侧剖面翼型升力增量。

2.3 VLM-LRC方法评价

本文使用Fluent商业软件建立Rec=3.0×105条件下的FX 63-137低雷诺数翼型气动特性参考数据库，耗时约3 h。表1为分别使用VLM-LRC方法和CFD方法完成该雷诺数条件下干净机翼算例及单螺旋桨/机翼算例所需要的计算模型涡格量或网格量，以及各自对应消耗的机时(所有计算均在同一台4核8线程计算机中进行)。其中，CFD方法计算耗时是在达到与VLM-LRC方法计算结果相同精度时所耗费的时间。可以看出，在达到相同计算精度的前提下，VLM-LRC方法计算效率相对CFD方法显著较高。

表1 气动性能计算消耗时间Table 1 Aerodynamic performance calculated time

3 太阳能无人机机翼平面形状设计

通过2.2节研究可以看出，利用螺旋桨滑流影响来降低机翼诱导阻力具有一定的可能性与可行性。螺旋桨滑流影响下的机翼升力明显增大，也就是说，在提供相同升力的情况下，考虑螺旋桨滑流影响的机翼环量需求更低，继而其诱导阻力也会进一步降低。因此，本节将以机翼展向升力分布接近理论椭圆形曲线(对应诱导阻力最小)为目标进行机翼平面形状设计研究。基于VLM-LRC方法针对某型手抛式太阳能无人机采用直接优化设计与反设计相结合的方法，分别对螺旋桨停机以及螺旋桨工作2种情况下的机翼平面形状进行优化设计，并采用CFD方法对优化设计结果进行验证、比较及分析。

3.1 手抛式太阳能无人机

图5(a)给出了某型手抛式太阳能无人机构型示意图。该无人机采用了双尾撑和倒V尾的布局形式，机翼采用FX 63-137平直机翼，弦长0.3 m，展长4.5 m，2个螺旋桨均为直径0.3 m的X1螺旋桨，以0° 安装角对称安装在机翼两侧正前方0.2 m的位置，螺旋桨间距为1.8 m，左、右螺旋桨分别沿顺气流方向逆时针、顺时针旋转。该无人机巡航(设计)状态选为：远场来流速度V∞=15 m/s，飞行高度H=500 m，机翼弦长雷诺数Rec=3.0×105，螺旋桨工作转速为4 000 r/min。

为提高计算及设计效率(见图5(b))，对该无人机模型进行简化，忽略机身、吊舱以及倒V尾的影响，且仅针对螺旋桨同步旋转时其滑流影响下的机翼平面形状进行优化设计。

图5 手抛式太阳能无人机构型示意图Fig.5 Schematic of hand-throw solar-powered unmanned aerial vehicle

3.2 气动优化设计

图6为相关的低雷诺数机翼平面形状具体优化设计流程，该研究内容共由2轮设计完成：第1轮直接优化设计，以及第2轮反设计。

图6 优化设计流程Fig.6 Flowchart of optimization design process

3.2.1 直接优化设计

在直接优化设计过程中，首先，结合实际工程需求，以太阳能无人机机翼面积，也即太阳能电池铺设面积不变为约束条件，对机翼平面形状进行参数化建模，建模过程中采用三次B样条曲线[29-30]对机翼后缘(Trailing Edge, TE)边界形状进行控制，沿单侧机翼展向共选取6个控制点，各控制点参数设计空间为(-3.0，2.0)，当控制前后机翼面积变化量ΔSw≤0.001 m2时，即可认为满足机翼面积不变的约束条件；然后，采用VLM-LRC方法直接并行求解当前模型在0° 和2° 迎角时的机翼气动力系数，采用线性插值方法计算出CL=0.9时的机翼阻力系数CD，同时计算出2° 迎角时优化机翼展向升力分布曲线与理论椭圆形分布曲线的均方根误差ε；最后，采用NSGA-II遗传算法[31-32]为搜索器进行直接多目标数值优化，其中种群个体数设定为12，遗传代数设定为20，交叉率设定为0.9，变异率设定为0.5。

直接优化目标设定为：在太阳能无人机巡航状态提供升力系数CL=0.9的约束条件下，获得最优CD及ε。该优化设计问题可以表示为

Objective function:minJ=ω1CD+ω2ε

式中：ω1和ω2为权重系数，取值均为0.5。

相比于传统多目标优化算法，NSGA-II遗传算法运行一次即可获得多个Pareto最优解，且其运算效率高，解集亦具有良好的分布性，特别是对目标个数较少的问题有很好的表现。图7给出了

图7 NSGA-II遗传算法优化搜索迭代过程 Fig.7 Iteration history of optimization search by NSGA-II genetic algorithm

本文机翼平面形状寻优迭代历程。可以看出，经过优化设计系统的优化搜索，可以得到阻力特性明显改善的优化结果。

3.2.2 反设计

在直接优化设计结果的基础上进一步对机翼诱导阻力系数进行反设计，将机翼展向升力分布按理论椭圆形曲线分布为目标，以机翼展向升力分布曲线与理论椭圆形曲线的逼近程度，也即均方根误差ε的大小作为衡量标准，对机翼弦长分布进行修正设计。

首先，判断当前模型机翼展向升力分布曲线与理论椭圆形曲线之间的均方根误差ε是否满足设计终止条件(ε≤0.005)，若是，则输出设计结果，若否，则对当前模型机翼各控制点处截面升力系数与最优椭圆形理论值之间的关系进行判断及分析，基于该关系对各控制点参数进行步长为 0.1 的余量修正；然后，更新计算模型，如此反复；最后，设计出满足要求的机翼平面形状。

3.3 优化结果分析

图8为优化前后机翼平面形状对比。为了便于区分，将原始机翼、螺旋桨停机状态优化机翼以及螺旋桨工作状态优化机翼分别命名为Original wing、W-optimal wing及P-optimal wing以示区分。可以看出，当不考虑螺旋桨滑流影响时，优化机翼平面形状呈现出类椭圆形特征，翼尖区域机翼弦长显著减小；而当考虑螺旋桨滑流影响时，整个机翼后缘形状呈波浪形，翼尖区域和螺旋桨下游区域机翼弦长均显著减小。

图9为2° 迎角下3种机翼在螺旋桨停机及工作状态下的展向升力系数分布曲线与理论椭圆形升力系数分布曲线的对比。其中，各机翼沿展向被均分为50份，且螺旋桨靠近机翼对称面的一侧为上洗侧，靠近机翼翼尖的一侧为下洗侧。考虑到本文引入理论椭圆形升力系数分布曲线的主要目的是为了通过与计算升力系数分布曲线之间的对比为后续机翼设计提供指导，并不具备实际意义，因此在椭圆形曲线构建过程中不考虑曲线与x轴所围成的面积，也即椭圆形机翼理论升力与当前机翼计算升力一致，其半短轴取值始终为1.0，半长轴取值则始终为当前机翼中心截面翼型计算升力系数。

图8 优化前后机翼平面形状对比Fig.8 Comparison of wing planform between before and after optimization

图9 2° 迎角下3种机翼展向升力系数分布对比Fig.9 Comparison of span-wise lift coefficient distribution among three wings (α=2°)

由图9可以看出：① 螺旋桨滑流主要改变其对应下游有限区域内的机翼计算升力分布，其对机翼翼尖及翼根处的气动力分布几乎没有任何影响，同时当前状态下螺旋桨对远场来流的旋转效应相对更强，左右螺旋桨下游机翼展向升力分布形态均呈现出典型的气流上洗侧升力增大、下洗侧升力减小特征。② 当螺旋桨停机时，原始机翼沿展向的升力系数分布曲线与椭圆形理论曲线差别显著，而在螺旋桨停机状态下进行设计的优化机翼“W-optimal wing”沿展向的计算升力系数分布曲线与椭圆形理论曲线几乎完全吻合，而优化机翼“P-optimal wing”沿展向升力系数分布曲线则与其呈波浪形分布的机翼弦长相对应，在机翼弦长显著减小的螺旋桨下游区域内计算升力系数相对理论亦较小，此时原始机翼、优化机翼“W-optimal wing”以及优化机翼“P-optimal wing”3种机翼展向升力系数分布曲线与理论椭圆形分布曲线的均方根误差依次为0.024 2、0.003 6和0.013 3。③ 当螺旋桨工作时，螺旋桨滑流的加速作用十分显著，3种机翼在各自螺旋桨下游区域均与螺旋桨旋转方向相对应地呈现出较为显著的气流上洗及下洗特征，原始机翼及优化机翼“W-optimal wing”展向升力分布均方根误差亦显著增大，而优化机翼“P-optimal wing”由于其呈波浪形分布的机翼弦长能够与螺旋桨滑流影响形成有效互补，其沿展向的计算升力系数分布曲线与椭圆形理论曲线吻合相对较好，此时3种机翼展向升力系数分布曲线与理论椭圆形分布曲线的均方根误差依次为0.027 4、0.007 1 和0.004 7。

表2优化前后机翼阻力系数对比(CL=0.9)

Table2Comparisonofwingdragcoefficientsbeforeandafteroptimization(CL=0.9)

PropellerWingΔCDi/countΔCDp/countΔCD/countΔCD/%OffOriginal193．0143．0336．0W⁃optimal172．2148．4320．6-4．58P⁃optimal184．4141．3325．7-3．06OnOriginal190．0154．3344．3W⁃optimal186．4174．3360．7+4．76P⁃optimal176．6150．0326．6-5．14

表2为螺旋桨停机及工作2种状态下在达到设计升力系数CL=0.9时的优化前后机翼计算阻力系数(1 count=10-4)对比。其中，总阻力系数变化百分比均是以相同状态下的原始机翼总阻力为参考进行计算。

可以看出：① 在螺旋桨工作状态下，原始机翼诱导阻力系数相对螺旋桨停机状态稍有减小，这表明与翼尖涡方向相反的螺旋桨滑流可以在一定程度上促使机翼诱导阻力降低，但机翼型阻相对显著增大，整体上机翼总阻力计算值增大达8.3 counts。② 在螺旋桨停机状态下，优化机翼“W-optimal wing”计算诱导阻力相对同状态原始机翼减小达20.8 counts，但计算型阻相对增加约5.4 counts，计算总阻力相对减小4.58%，这表明以机翼展向升力分布形态为衡量标准的设计思路可行有效。但在螺旋桨工作状态下，优化机翼“W-optimal wing”并没有表现出如螺旋桨停机状态一般的优秀特性，其诱导阻力相对同状态原始机翼仅有3.6 counts的减小量，而由于螺旋桨下游机翼弦长的增加，对应其剖面升力系数有所增加，从而导致机翼型阻相对显著增大，机翼总阻力亦显著增大达4.76%。也就是说，在螺旋桨工作状态下，该优化机翼阻力特性相比原始机翼反而有所恶化。③ 在螺旋桨停机状态下，优化机翼“P-optimal wing”计算诱导阻力及型阻均相对同状态原始机翼稍有减小，最终机翼总阻力相对减小3.06%。而在螺旋桨工作状态下，优化机翼“P-optimal wing”计算诱导阻力减小量显著增大至13.4 counts，同时由于螺旋桨下游机翼弦长的减小，其计算型阻亦相对稍有减小，最终优化机翼“P-optimal wing”计算总阻力相对减小达17.7 counts，占同状态原始机翼总阻力的5.14%。

由此可见，对于螺旋桨类飞机而言，传统未考虑螺旋桨滑流影响而仅针对干净机翼进行气动设计的思路在实际应用中存在较严重问题，而考虑实际工作状态螺旋桨滑流影响下的机翼气动优化设计显然更具有意义。

3.4 优化结果验证

考虑到实际工程中，三次曲线形状的机翼后缘边界将为机翼结构的设计、加工与制造，以及机翼舵面的设计、安装与操纵等带来困难。本文在优化机翼实体建模过程中，针对性地对机翼后缘形状进行了修正——利用若干条直线段代替曲线边界。图10为优化机翼“P-optimal wing”后缘修正前后的平面形状对比。图中：左侧为B样条曲线型后缘模型，右侧为修正后的折线型后缘模型。

图11为采用CFD方法计算得到的实际螺旋桨工作状态下的优化前后几种机翼升阻特性计算结果对比。图中：B-spline和Polyline分别代表曲线型后缘和折线型后缘的优化机翼。

由图11可以看出：① 折线型后缘与曲线型后缘对机翼气动力结果影响并不明显，全计算迎角范围内二者之间的升阻力相对误差始终在0.3% 以内。② 相对于原始机翼，2种优化机翼在各迎角下的计算升阻力均有所减小，而在设计状态机翼CL=0.9时，优化机翼“W-optimal wing”阻力系数相对原始机翼增大约20.13 counts，占原始机翼总阻力的6.29%，而优化机翼“P-optimal wing”阻力系数则相对减小约19.52 counts，占原始机翼总阻力的6.10%。

由此可见，CFD方法计算结果及气动力变化趋势与VLM-LRC方法计算结果较为一致，这表明本文基于VLM-LRC方法发展的耦合螺旋桨滑流影响的机翼优化设计思路及方法有效可靠。

图10 优化机翼后缘修正示意图Fig.10 P-optimal wing TE-corrected view

图11 优化前后机翼CFD方法计算结果对比Fig.11 Comparison of CFD numerical results before and after wing optimization

4 结 论

1) 以CFD准定常求解方法为参考，本文发展的涡格法准定常求解程序及低雷诺数修正方法针对低雷诺数条件下的多螺旋桨/机翼相互气动干扰问题具备较高的数值计算精度以及突出的数值计算效率，且其计算结果能够在一定程度上准确反映螺旋桨滑流影响下机翼表面力分布特征。

2) 在本文研究的低雷诺数状态下，螺旋桨滑流对机翼气动特性的影响表现为增升增阻，同时由机翼展向力分布还可以看出螺旋桨对气流的加速效应极为显著，而旋转效应则相对较弱。

3) 本文基于VLM-LRC方法提出的通过控制机翼表面形状来获得最优椭圆形升力分布的设计思路有效可靠，在螺旋桨停机及工作状态下分别进行优化设计的机翼相对同等状态原始机翼诱导阻力可降低达20.8 counts和13.4 counts。而最终采用CFD方法进行验证的优化机翼“P-optimal wing”总阻力可相对降低达19.52 counts，减阻效果显著。

4) 对于多螺旋桨驱动的太阳能无人机而言，螺旋桨滑流影响将为全机带来较为显著的附加气动力，这可能会导致传统仅仅针对主要升力面进行优化设计的无人机飞行性能偏离设计点，甚至产生负面作用。因此，耦合动力系统气动影响的无人机优化设计思路更具实际工程意义。

参 考 文 献

[1] ROMEO G, FRULLA G, CESTINO E, et al. Heliplat: Design, aerodynamic, structural analysis of long-endurance solar-powered stratospheric platform[J]. Journal of Aircraft, 2004, 41(6): 1505-1520.

[2] EHAB A E, COLIN P B. Experimental investigation of the effect of propeller slipstream on boundary layer behavior at low Reynolds number[J]. AIAA-Applied Aerodynamics Conference, 2000, 194(10): 267-276.

[3] 王科雷, 周洲, 甘文彪, 等. 太阳能无人机低雷诺数翼型气动特性研究[J]. 西北工业大学学报, 2014, 32(2): 163-168.

WANG K L, ZHOU Z, GAN W B, et al. Studying aerodynamic performance of the low-Reynolds-number airfoil of solar energy UAV[J]. Journal of Northwestern Polytechnical University, 2014, 32(2): 163-168 (in Chinese).

[4] GAVIN K A, ROBERT W D, MICHAEL S S. Propeller induced flow effects on wings at low Reynolds numbers: AIAA-2013-3193[R]. Reston： AIAA, 2013.

[5] KROO I. Propeller-wing integration for minimum induced loss[J]. Journal of Aircraft, 1986, 23(7): 561-565.

[6] VELDHUIS L L M, HEYMA P M. A simple wing optimization code including propeller effects[C]//21st Congress of the International Council of the Aeronautical Sciences. [S.l.]： ICAS, 1998.

[7] VELDHUIS L L M, HEYMA P M. Aerodynamic optimization of wings in multi-engined tractor propeller arrangements[J]. Aircraft Design, 2000, 3(3): 129-149.

[8] VELDHUIS L L M. Propeller/wing aerodynamic interference[D]. Delft: Delft University of Technology, 2005.

[9] RAKSHITH B R, DESHPANDE S M, RODDAM N, et al. Optimal low-drag wing planforms for tractor-configuration propeller-driven aircraft[J]. Journal of Aircraft, 2015, 52(6): 1-11.

[10] 乔宇航, 马东立, 邓小刚. 基于升力线理论的机翼几何扭转设计方法[J]. 北京航空航天大学学报, 2013, 39(3): 320-324.

QIAO Y H, MA D L, DENG X G. Wing geometric twist design method based on lifting-line theory[J]. Journal of Beijing University of Aeronautics and Astronautics, 2013, 39(3): 320-324 (in Chinese).

[11] 徐家宽, 白俊强, 黄江涛, 等. 考虑螺旋桨滑流影响的机翼气动优化设计研究[J]. 航空学报, 2014, 35(11): 2910-2920.

XU J K, BAI J Q, HUANG J T, et al. Aerodynamic optimization design of wing under the interaction of propeller slipstream[J]. Acta Aeronautica et Astronautica Sinica, 2014, 35(11): 2910-2920 (in Chinese).

[12] 陈广强, 白鹏, 詹慧玲, 等. 高空长航时无人机螺旋桨滑流效应影响研究[J]. 飞机设计, 2014, 34(4): 1-9.

CHEN G Q, BAI P, ZHAN H L, et al. Numerical simulation study on propeller slipstream effect on high altitude long endurance unmanned air vehicle (HALE UAV)[J]. Aircraft Design, 2014, 34(4): 1-9 (in Chinese).

[13] 陈广强, 白鹏, 詹慧玲, 等. 一种推进式螺旋桨无人机滑流效应影响研究[J]. 空气动力学学报, 2015, 33(4): 554-562.

CHEN G Q, BAI P, ZHAN H L, et al. Numerical simulation study on propeller slipstream effect on unmanned air vehicle with propeller engine[J]. Acta Aerodynamica Sinica, 2015, 33(4): 554-562 (in Chinese).

[14] 孙美建, 詹浩. Kriging模型在机翼气动外形优化中的应用[J]. 空气动力学报, 2011, 29(6): 759-764.

SUN M J, ZHAN H. Application of Kriging surrogate model for aerodynamic shape optimization of wing[J]. Acta Aerodynamica Sinca, 2011, 29(6): 759-764 (in Chinese).

[15] 穆雪峰, 姚卫星, 余雄庆, 等. 多学科设计优化中常用代理模型的研究[J]. 计算力学学报, 2005, 22(5): 608-612.

MU X F, YAO W X, YU X Q, et al. A survey of surrogate models used in MDO[J]. Chinese Journal of Computational Mechanics, 2005, 22(5): 608-612 (in Chinese).

[16] ESMAEEL E, MEHRAN T, SAMAN N. Aerodynamic performance of Parastoo UAV[J]. Aircraft Engineering and Aerospace Technology: An International Journal, 2013, 85(2): 97-103.

[17] 甘文彪, 周洲, 许晓平. 仿生全翼式太阳能无人机气动数值模拟[J]. 航空学报, 2015, 36(10): 3284-3294.

GAN W B, ZHOU Z, XU X P. Aerodynamic numerical simulation of bionic full-wing typical solar-powered unmanned aerial vehicle[J]. Acta Aeronautica et Astronautica Sinica, 2015, 36(10): 3284-3294 (in Chinese).

[18] 甘文彪, 周洲, 许晓平. 仿生全翼式太阳能无人机分层协同设计及分析[J]. 航空学报, 2016, 37(1): 163-178.

GAN W B, ZHOU Z, XU X P. Multilevel collaboration design and analysis of bionic full-wing typical solar-powered unmanned aerial vehicle[J]. Acta Aeronautica et Astronautica Sinica, 2016, 37(1): 163-178 (in Chinese).

[19] MICHAEL D P, BRIAN J G. Wing aerodynamic analysis incorporating one-way interaction with distributed propellers： AIAA-2014-2852[R]. Reston： AIAA, 2014.

[20] RANGWALLA A A, WILSON L N. Application of a panel code to unsteady wing-propeller interference[J]. Journal of Aircraft, 1987, 24(8): 568-571.

[21] 迟永一. 基于CFD方法的轻型公务机雷诺数修正研究[J]. 黑龙江科技信息, 2014(17): 122.

CHI Y Y. A study of CFD-based Reynolds correction for light business jets[J]. Heilongjiang Science and Technology Information, 2014(17): 122 (in Chinese).

[22] 张辉, 李杰, 龚志斌. 某运输机全机构型雷诺数效应阻力修正研究[J]. 飞行力学, 2012, 30(1): 20-24.

ZHANG H, LI J, GONG Z B. Drag correction to higher Reynolds number for a transport aircraft[J]. Flight Dynamics, 2012, 30(1): 20-24 (in Chinese).

[23] 王巍, 李东升, 刘春. 基于小扰动理论的桨叶叶素气动载荷计算方法[J]. 北京航空航天大学学报, 2016, 42(2): 294-302.

WANG W, LI D S, LIU C. Calculation method of blade element aerodynamic loads based on small perturbation theory[J]. Journal of Beijing University of Aeronautics and Astronautics, 2016, 42(2): 294-302 (in Chinese).

[24] WALTERS D K, COKLJAT D. A three-equation eddy-viscosity model for Reynolds-averaged Navier-Stokes simulations of transitional flows[J]. Journal of Fluids Engineering, 2008, 130(1): 1-14.

[25] WALTERS D K, LEYLEK J H. Impact of film-cooling jets on turbine aerodynamic losses[J]. Journal of Turbomachinery, 2000, 122(3): 537-545.

[26] VOLINO R J. A new model for free-stream turbulence effects on boundary layers[J]. Journal of Turbomachinery, 1998(120): 613-620.

[27] HSIAO F B, LIU C F, TANG Z. Experimental studies of airfoil performance and flow structures on a low Reynolds number airfoil[C]//19th AIAA, Fluid Dynamics, Plasma Dynamics, and Lasers Conference, Honolulu, 1987.

[28] ABTAHI A A, MARCHMAN J F. Aerodynamics of an aspect ratio 8 wing at low Reynolds number[J]. AIAA Journal, 1985, 22(7): 628-634.

[29] 张家树, 杨恬. 三次B样条曲线的快速生成算法[J]. 西南师范大学学报(自然科学版), 1997, 22(6): 644-647.

ZHANG J S, YANG T. A fast generating algorithm of gubic B-spline curve[J]. Journal of Southwest China Normal University (Natural Science), 1997, 22(6): 644-647 (in Chinese).

[30] 荣焕宗, 张伟荣, 张蔚. 带有面积约束的B样条曲线拟合方法[J]. 应用数学与计算数学学报, 1990, 4(2): 19-25.

RONG H Z, ZHANG W R, ZHANG W. B-spline curve fitting method with an area constraint[J]. Communication on Applied Mathematics and Computation, 1990, 4(2): 19-25 (in Chinese).

[31] 王小平, 曹立明. 遗传算法-理论、应用与软件实现[M]. 西安： 西安交通大学出版社, 2002.

WANG X P, CAO L M. Genetic algorithm-theory, application and realization[M]. Xi’an: Xi’an Jiaotong University Press, 2002 (in Chinese).

[32] 谭艳艳. 几种改进的分解类多目标进化算法及其应用[D]. 西安： 西安电子科技大学, 2013.

TAN Y Y. Several modified decomposition-based multi-objective evolutionary algorithms and their applications[D]. Xi’an: Xidian Univeristy, 2013 (in Chinese).

(责任编辑： 鲍亚平)

Aerodynamic design of low-Reynolds-number wing taking intoaccount the multiple propellers induced effects

WANGKelei1,2,ZHOUZhou1,2,*,ZHUXiaoping2

1.CollegeofAeronautics,NorthwesternPolytechnicalUniversity,Xi’an710072,China2.LaboratoryofScienceandTechnologyonUAV,NorthwesternPolytechnicalUniversity,Xi’an710065,China

Based on a certain hand-throw solar-powered unmanned aerial vehicle (UAV), the optimization design approaches for low-Reynolds-number wing coupled with multiple propellers induced effects are studied. The corresponding quasi-steady procedure based on the vortex lattice method (VLM) of lifting line theory and the low Reynolds correction (LRC) method based on the reference airfoil aerodynamic properties database are developed to simulate the multiple propellers/wing aerodynamic interference at low Reynolds numbers. The reliability and accuracy of the simplified numerical method (VLM procedure and LRC method) are testified with several cases studies and their comparison with experimental results. Both the direct optimization design and inverse design of the simplified hand-throw solar-powered UAV model in tractor configuration are conducted, and the optimization results are examined with high-accuracy CFD technique. It shows that (a) the low-Reynolds-number flow can be simulated by the VLM-LRC method efficiently and accurately; (b) the aerodynamic properties of the optimal wing cannot be improved when the propeller slipstream effect is not taken into consideration in the conventional design approach; (c) the wing drag performance can be greatly improved with the optimization approach that takes into account the multiple propeller slipstream effects, and the optimized wing has a drag reduction of 19.52 counts at the design state.

hand-throw solar-powered unmanned aerial vehicle; low-Reynolds-number wing; multiple propellers/wing aerodynamic interference; vortex lattice method; low Reynolds correction; optimal design; inverse design

2016-09-23；Revised2016-12-03；Accepted2016-12-19；Publishedonline2017-01-041433

URL：www.cnki.net/kcms/detail/11.1929.V.20170104.1433.004.html

s：CivilAircraftProject(MIZ-2015-F-009);ShaanxiProvinceScienceandTechnologyProject(2015KTCQ01-78)

2016-09-23；退修日期2016-12-03；录用日期2016-12-19； < class="emphasis_bold">网络出版时间

时间：2017-01-041433

www.cnki.net/kcms/detail/11.1929.V.20170104.1433.004.html

民机专项 (MIZ-2015-F-009)； 陕西省科技统筹 (2015KTCQ01-78)

*

.E-mailzhouzhou@nwpu.edu.cn

王科雷， 周洲， 祝小平． 耦合多螺旋桨滑流影响的低雷诺数机翼设计J． 航空学报,2017,38(6):120813.WANGKL,ZHOUZ,ZHUXP.Aerodynamicdesignoflow-Reynolds-numberwingtakingintoaccountthemultiplepropellersinducedeffectsJ.ActaAeronauticaetAstronauticaSinica,2017,38(6):120813.

http://hkxb.buaa.edu.cnhkxb@buaa.edu.cn

10.7527/S1000-6893.2016.120813

V211

A

1000-6893(2017)06-120813-13

*Correspondingauthor.E-mailzhouzhou@nwpu.edu.cn