高前进比旋翼的超常规叶素气动特性分析
文/丁志伟
不负责任的声明：本文是一篇偏学术派的硬核专业笔记，如因强行阅读本文诱发头晕目眩、昏昏欲睡等症状，丁某概不负责。
正经来说，本文本来确实是一份笔记稿，不过成稿已经有一些时间了，最近在回顾高前进比旋翼[注1]相关内容，顺手整理出来与大家交流。主要记录了两个部分的内容，第一个部分阐述了一种关于叶素段如何对应360°翼型气动系数C81表的小技巧，这个小技巧对于刚开始进行旋翼气动分析或者从低阶气动模型逐渐转向高阶气动模型的读者朋友也许有一定的参考价值；第二部分内容对高速直升机旋翼在高前进比状态下桨叶诱导的一些涡运动的特性，其中涉及两篇参考文献，如果没有做过高速直升机旋翼气动方面的研究而想要较好地理解这部分内容的话，建议还是先大致略读以下这两篇参考文献最佳。
作为一篇笔记稿，本文更多的是一些我自己在进行一些问题的探究过程中的思考，所以其置信度低于通过同行评议的学术论文，大致上与预印本相当，望有深度阅读需求的读者朋友知悉。
[注1] 所谓前进比即是$\mu = Vf/Vtip$，其中$\mu$就是前进比，$Vf$是指飞行器的飞行速度，而$Vtip$是指旋翼的桨尖线速度，对于常规直升机来说，其前进比的极限值一般决不可能超过0.4。高前进比就是指的$\mu$值比较高的情况，对于到底多高能被称为高前进比状态，不同的研究人员有不同的说法，不过$\mu = 0.6$（Graham Bowen-Davies；2018）是一个认可度比较高的前进比边界，本文同样认为0.6以下为低前进比，以上则为高前进比状态。
1 叶素升阻力系数与风洞试验数据的对应关系策略
1.1 叶素升阻力求解中的疑难点概述
在旋翼的气动分析中，翼型的升阻力特性获取可谓是相当重要的一环，准确的翼型升阻力特性可以提高计算结果的置信度和精度，而不准确的翼型升阻力特性则可能导致配平发散甚至得到错误的结果，在高速高前进比的状态下，气流的强不对称性愈发明显，旋翼气动环境更为复杂，要对旋翼气动特性进行准确计算，首当其冲的就是获取准确的叶素升阻力系数。
为了准确获得叶素升阻力系数，最常用的方法就是先获取对应翼型的360°吹风数据，然后通过插值的方式来得到所需要的结果。但是在代码编写的过程中，我发现常规的叶素气动力系数插值方法，常常存在针对某种来流状态或某种独特翼型的简化假设，导致这些方法无法用于流动复杂、迎角变化幅度大的高前进比状况，为了解决这一问题，我采取了一种全矢量求解迎角和气动力方向的计算策略。
1.2 风洞试验中的迎角与升阻力方向
在风洞实验中，常规的做法是保持来流不变，以翼型前缘正对来流为零度迎角，绕气动中心抬头旋转，将翼型从前缘至后缘连成一条矢量线，该矢量与来流方向的夹角就是翼型与来流的迎角 α ，每隔一个给定的角度变化小量，记录一次稳态的升阻力系数，最后统计成表就是常见的翼型升阻力系数C81表。
在这种试验方式下，阻力的正方向很容易确定，因为阻力矢量与来流矢量始终保持一致，这是毋庸置疑的，而升力的正方向则由以下公式
$$\vec{V}_a \times \vec\alpha$$
来决定，其中 $\vec{V}_a$ 表示来流速度矢量，$\vec\alpha$表示迎角方向矢量（由零升迎角和右手法则确定），其物理意义实则可以是从茹科夫斯基升力定律中来的：
$$T = \rho_a \vec{v}_a \times \Gamma$$
$$\Gamma = \oint \vec{v} ds$$
上式中，Γ为绕翼型环量，通过它的计算公式可以看出，环量方向实则与翼型表面速度流动方向一致，由伯努利定律，风洞中产生正升力的翼型上表面流动要大于下表面，因而Γ的方向实则与迎角转动的方向一致，由此便可以得到升力的正方向，而在风洞实验中，这个方向一般取为竖直向上。
1.3 理论计算中的气动迎角与升阻力方向
在实际的旋翼叶素气动力计算中，来流的方向不像风洞那样是保持不变的，旋翼的角度也不是某一个确定的静态值。处于飞行来流、旋转来流、诱导来流等多重流动作用中的叶素，其气动迎角和升阻力方向的确定远比风洞要复杂。但是，为了计算方便，在常规直升机的叶素气动力计算中，迎角（如下图所示）一般由安装角 θ 与来流角φ 之差给出，式中uT , uP 分别为切向（周向）来流和轴向来流。
$$ \alpha = \theta - \phi $$
$$\phi = tan^{-1} \frac{u_p}{u_T}$$
▲叶素来流迎角和升阻力方向示意图
对于常规直升机而言，除靠近桨根的少部分桨叶段之外（这部分叶素气动力在常规旋翼计算中往往会进行简化处理），绝大部分叶素周向来流都远大于诱导速度，总体而言气动迎角都会落在-15~20°区间，也就是说在大部分状况下，升力系数与迎角是线性关系，升力的方向也比较容易确定，因而这样的简化求解方式是合理的；但是对于高前进比旋翼而言，反流区较大，后行侧大部分桨叶都将进入反流区，周向来流都是从后缘吹向前缘，而反流区内，靠近桨尖的大部分区域，轴向诱导速度相比周向来流，其量级不再有显著差异，因而在旋翼周期运动中，其叶素来流迎角变化幅度将会较为巨大，这种情况下，常规的迎角计算方法将不再适用，比如下图所示的状况。
▲翼型后缘来流示意图
观察上图，从直观可以判断出其气动迎角应当是落在180°附近的区域，但是直接用安装角与来流角之差难以计算准确迎角，而且其升力方向与风洞升力系数之间的关系也并不明显，无法用常规方法中的某一解析公式来进行计算。
1.4 理论计算与风洞数据一一对应的策略
针对上述问题，一种解决方案是采用分段迎角求解方法，考虑来流的不同方向，将迎角计算域分为多个区域（如分为四个——前缘上方、前缘下方、后缘上方、后缘下方），以此来计算正确的迎角和升阻力系数。这种方法针对特定的翼型能取得较好的效果，实则上却存在通用性的问题，不同的翼型将要选择不同的分段方式，比如说，对于存在零升迎角的非对称翼型，其分段区域的划分显然与对称翼型不同，而在高前进比旋翼的外形优化设计中，将应用多种翼型分段布置，如果用这种分段方式，其通用性和效率都不佳，也不利于研究工作的进一步展开。
为解决上述问题，采用了一种全矢量求解的策略来处理迎角求解和升阻力方向计算的问题，经验证，该方法能做到将实际情况与风洞数据一一对应起来，通用性较好。
▲惯性坐标系组示意图
该方法的基本思路是，首先建立一整套惯性坐标系组（如上图所示），包括桨毂坐标系（下标hub），桨叶旋转坐标系（下标ro），桨叶挥舞坐标系（下标β），叶素当地坐标系（下标be，翼型气动中心为原点，前缘指向后缘为正方向），并计算各坐标系之间的转换矩阵，将速度、控制点等量都作为矢量的形式表示在各个坐标系中，而最后在求解旋翼迎角的时候，将叶素来流从其他坐标系中转换到叶素当地坐标系中，并且表示为矢量的形式，通过计算其与叶素翼型前缘到后缘连线这一矢量之间的夹角来确定准确的来流迎角（为服务于矢量计算，保证角度变化的连续性，本方法中迎角取0°到360°，而非传统的-180°到180°[注2]），在这种坐标系下，升力的矢量方向将由下述方式确定
迎角的确定
v_be = rotobe(v_air, β, θ) # 将来流速度矢量转化到叶素坐标系下
v_be_2d = [v_be[2],v_be[3]] # 在叶素坐标系下，得出速度在叶素剖面二维矢量
α = aoaget(v_be_2d) # 来流矢量与翼型前后缘连线矢量夹角即为实际气动迎角（0-360°）
cl,cd = clcdget(α)
升力方向的确定
if fcl(α) >= 0.0# 判断升力系数方向
lbe = [0, cos(α+π/2), sin(α+π/2)] # 若为正，则升力方向为迎角正转90°方向
else
lbe = [0, cos(α-π/2), sin(α-π/2)] # 若为负，则升力方向为迎角逆转90°方向
end
上述代码中，fcl 为升力系数插值函数，fcl(α) 得到的就是对应迎角的升力系数值，lbe 表示升力正方向的矢量，通过该方式，实际升力的方向能自动与风洞数据匹配，因而可采用下述计算公式计算叶素微段单位升力矢量
$$L_e = \vec{lbe} \times |0.5 \rho v^2_r c \times fcl(\alpha)|$$
式中，Le 为叶素微段单位升力，vr 为叶素微段总的来流速度，c为弦长。
采用上述策略，不仅保证了迎角计算的准确性，同时保证了叶素升阻力系数的插值与风洞数据一一对应起来，且具有良好的通用性，能应用到任意翼型的数据插值中。
[注2] 在我以前的代码中（Fortran），0360°方位角变化对于矢量求解代码编写比较友好，在我现在的代码（Julia）中，不管是按照0360°还是-180~180°来求解，代码编写都一样简单，哈哈，推荐一波Julia语言（虽然这本质上其实不是代码语法的作用，主要是对问题的理解更深入了）
2 叶素的超常规气动问题研究
2.1 高前进比流场的物理特性概述
在常规直升机旋翼计算中，由于旋翼反流区较小，在工程计算中往往直接被忽略，而在高前进比计算中（如下图所示），反流区扩散到整个后行侧，旋翼有一半的时间在其中转动，处于其中的旋翼将不产生升力或者产生负升力，这种超常规的流动特性对旋翼的气动载荷和旋翼配平都产生了直接的影响，须谨慎考虑和处理。
▲旋翼前进比增大，反流区逐渐扩大示意图
在高前进比旋翼流场复杂物理现象中，其中最具有代表性的就是反流区中来流自尖锐后缘吹向前缘时，在后缘背风处形成的强涡结构以及由于前飞高速来流导致的桨叶径向强偏流效应（Yawed Flow)。强后缘分离涡生成之后，受到高速来流以及径向偏流的影响，与桨叶一起快速移动，导致大部分后行侧桨叶段都受到其影响，而受其影响的叶素微段气动力特性与常规状态下显然存在种种差异。就目前而言，这方面的相关研究尚且较少，方法体系尚不完善。基于从物理现象本身出发进行探索的考虑，本文采用了基于格子-玻尔兹曼（Lattice Boltzmann Method, LBM）的CFD方法进行流场可视化，以此对高前进比旋翼叶素的超常规物理特性进行了一些研究。
2.2 反流区后缘涡的物理特征
观察旋翼高速高前进比三维流场，可以看到大约在220°方位角的时刻，受高前进比来流的影响，后行侧桨叶的尾随涡无法再卷起为集中的桨尖涡，而靠近桨根的叶素微段处，来流从后缘吹向前缘，生成尖锐后缘涡（Sharp Edge Vortex, SEV），此时受到较强的朝向桨根的径向偏流影响，但是前期后缘涡黏着在后缘表面，拉伸效应不显著；而随着桨叶装过270°，后缘涡又受到较强的朝向桨尖的径向偏流影响，此时后缘涡已经不再是黏附状态，所以拉伸效应显著，加速了后缘涡的扩散。
▲涡量云图（观察点为桨盘下方）
▲桨盘下方平面内涡量演化图
从涡量演化图中还可以发现，高速来流中桨叶柄的后缘拖出了较强的集中涡，这也从一个方面解释了高前进比旋翼气动性能理论计算中桨盘后向力（阻力）比试验数据要小的现象，高速来流中的桨叶柄显然产生了较强的后向阻力，因而在高前进比气动分析中计入这个现象将提升计算的置信度，此外，该现象的出现也要求我们在进行高前进比旋翼优化设计的过程中，桨叶柄的优化设计也要纳入考量。
▲后缘二维流场
结合三维涡量云图和二维流场图来观察可以发现，尖锐后缘涡自桨叶后缘生成后首先会紧密附着在后缘表面。此后随着桨叶沿方位角转动，后缘涡逐渐变强，涡核逐渐变大，随后与桨叶表面轻微分离，但仍然随着桨叶一起运动，大约300°方位角开始，后缘涡涡核逐渐扩散，汇入自由尾流中。
在普遍认可的尖锐后缘涡流动现象中，后缘涡生成后会从后缘尖角脱离，但就上述流场图所展示的内容来说，后缘涡会吸附在桨叶表面并跟随其移动一段时间之后才开始脱落，该现象直接关系到了后缘动态失速是否存在这一问题。
是否存在后缘失速问题说明
a. NACA 0012 翼型升力系数图
b. OA212翼型升力系数图
▲翼型升力系数图
上图分别为NACA 0012 和OA212翼型360°吹风数据图，从图中我们明显可以看到，静态吹风数据中，来流从后缘吹向前缘时，大约在170°迎角左右，存在由后缘分离导致的失速现象，从曲线中同样可以看出，翼型前缘流动分离也会导致失速，而在旋翼周期变距运动中，这一现象表现为动态失速（Dynamic Stall），对于前缘分离导致的动态失速，Leishman-Beddoes 的旋翼翼型动态失速模型计算简单、效率高、精度好因而应用广泛，取得了普遍的认可。
那么，后缘分离导致的失速现象会不会在高前进比旋翼后行侧桨叶的变距运动中演化为动态失速呢？Potsdam 等人[1]在2016年应用RCAS对UH-60A的旋翼分析中，认为，后缘动态失速是存在的，并且产生了后缘动态失速涡（如下图a所示，动态失速涡受偏流影响从后缘移向前缘，导致了图中出现的吸力峰），但是对其物理特性却无法进行准确分析，而在本文的CFD模拟中发现，从大约220°方位角产生尖锐后缘涡开始，由于旋翼的转动和径向流动的影响，后缘涡经历了从黏附后缘表面到分离的过程，但其始终在桨叶后缘表面附近，在前飞来流的偏流影响下，随桨叶一起运动到约300°扩散为止，也就是没有出现动态失速这一过程，Potsdam 等人观察到的“动态失速涡”很有可能就是径向运动的尖锐后缘涡，Hiremath 等人[2]本年度于乔治亚理工学院低速风洞中做的高前进比试验也佐证了这一观点。
a. Potsdam 观察到的“动态失速涡”在桨叶下表面导致的压力变化
b. Hiremath 的实验结果
▲理论计算和试验结果的对比
综上所述，在后续对后行侧的叶素升阻力系数计算中，必须谨慎对待该现象，可通过计入尖锐后缘涡的影响来对升阻力插值表进行修正，拟采用的做法是对静态数据中后缘分离失速段做平滑处理，延缓后缘失速迎角的插值区间。
2.3 强偏流效应的修正
▲偏流示意图
上图中Λ为偏流角，ut 为周向分量，V∞ 为前飞来流，则偏流角计算公式如下，其中ψ为方位角。
$$\Lambda = cos^{-1} (u_t/V_\infty)$$
$$u_t = V_\infty sin(\psi) $$
在高速飞行的时候，来流速度V∞ 径向分量和周向分量随方位角变化曲线如下图所示，其中径向分量指向桨尖为正，周向分量从前缘吹向后缘为正。
▲前飞来流桨叶周向和径向分量随方位角变化示意图
由此可见，径向分量Ur 在高速高前进比飞行过程中，相比于旋翼转速而言已经不可忽略，因此在高速高前进比旋翼的叶素气动参数计算过程中不计入偏流效应，很可能得到错误的结果，影响配平和性能计算的准确性。前行侧旋翼的偏流效应下叶素的气动力系数的求解中可以通过由Smith等人[3]提出的偏流效应修正公式来计入片偏流效应的影响：
$$c_l(\alpha) = c_{l_2d}(\alpha cos^2\Lambda)/cos^2\Lambda$$
$$c_d(\alpha) = c_{d_2d} ( \alpha cos\Lambda)/cos\Lambda$$
$$c_m(\alpha) = c_{m_2d}(\alpha cos^2\Lambda)/cos^2 \Lambda$$
式中下标2d表示未考虑偏流效应状态下的升阻力系数值。
而本文中为了计入高速状态下后行侧反流区中偏流效应带来的影响，主要针对典型翼型在后行侧反流区中存在不同偏流速度时的升阻力系数进行了计算，而后通过数据拟合，形成偏流效应修正的半经验修正策略，其基本思想就是通过拟合受偏流影响下的升阻力系数随方位角变化函数以及偏流角随方位角变化函数（下图所示就是一种拟合方式），最后通过函数间的对应关系，计算出修正的升力系数值。
▲升力系数随方位角变化示意图（反流区）
▲偏流角随方位角变化示意图
2.4 小结
基于格子-波尔兹曼的CFD方法对流体介质的运动模拟和计算效率都要优于传统基于N-S方程的CFD方法，但是在物面网格划分的细致程度和求解精度方面略逊一筹。
无论哪种CFD方法，从目前来看，受限于湍流模型和数值耗散，还是只能对高速高前进比旋翼的强不对称性周期流动进行定性分析或者精度要求较低的定量分析，尚且还无法达到常规旋翼分析那样的精度。
尽管如此，CFD的流场可视化结果仍然有助于我们对高前进比旋翼叶素的超常规物理特性进行更为深入的了解，而由此得到的结论足以为叶素的气动参数的求解和修正提供有益的方向，从而进一步推进高前进比旋翼的性能分析和优化设计。
参考文献
- Potsdam M, Datta A, Jayaraman B. Computational investigation and fundamental understanding of a slowed uh-60a rotor at high advance ratios[J]. Journal of the American Helicopter Society, 2016, 61(2): 1-17.
- Hiremath N, Shukla D, Komerath N. Physics of Reverse Flow on Rotors at High Advance Ratios[J]. arXiv preprint arXiv:1805.00136, 2018.
- Smith M J, Koukol B C G, Quackenbush T, et al. Reverse-and cross-flow aerodynamics for high-advance-ratio flight[C]. 2009.