自旋翼无人机因其飞行性能好、不受地形限制等优点, 被广泛用于地质探测、搜救侦查以及科学研究等领域。但是, 当在飞行过程中出现故障时, 可能会发生严重的安全事故甚至坠毁。目前, 已经广泛采用降落伞对固定翼无人机进行减速回收^[1], 而对旋翼机回收的研究工作较少。当失控状态的旋翼机处于垂直下降时, 旋翼产生的旋转流场^[2]加剧了前体尾流的非定常变化, 导致流向降落伞的流场变化紊乱.为探究旋翼扰动下物伞系统的流场特性以及伞衣气动性能变化, 需要对旋翼机/降落伞的非定常复合流场展开深入的研究。

目前, 旋翼流场的研究主要通过风洞试验和数值模拟, 风洞试验能够真实反映旋翼的工作状态^[3-4], 通过控制旋翼的状态参数获得准确的试验数据, 但存在洞壁效应^[5]、实验数据可重复性差^[6]、成本代价高等问题。随着计算流体力学(computational fluid dynamics, CFD)的发展, 数值模拟逐渐成为研究旋翼流场的重要手段, 该方法通过直接求解流场控制方程(Euler方程或Navier-Stokes方程)^[7-8], 能真实还原流场的细节变化, 具备较高的求解精度。如文献[9-11]采用CFD方法对悬停状态下的旋翼流场进行了数值模拟, 得到了桨叶附近真实的流场细节变化; 文献[12-13]观察到前飞状态下的旋翼桨尖涡, 分析了旋翼尾迹涡和下洗流的非定常流场特征变化; 文献[14]发现旋翼垂直下降流场会出现涡环状态, 随下降速度增加旋翼下洗流逐渐上移直至被来流抵消。

众多学者也广泛开展了旋翼复合流场的数值研究, 如文献[15-17]对倾转旋翼机/机翼的复合流场进行了CFD数值模拟, 精确捕捉到了流场的旋涡结构, 观察到旋翼/机身产生的“喷泉效应”^[15-16]现象; 文献[18-20]对旋翼与舰船的复合流场进行数值模拟, 发现旋翼绕流与舰船尾流发生叠加^[18], 会产生旋涡、回流、分离^[19]等多种复杂的流动形式, 在侧风的掺混下进一步形成涡流紊乱区, 涡流区范围和气流下洗明显增加, 甲板流场变化极不稳定^[20]。目前, 对旋翼复合流场的数值研究还有很多^[21-22], 但旋翼机/降落伞非定常复合流场的研究工作却很少。

本文基于PISO(pressure implicit split operator)算法对旋翼机/降落伞复合流场展开非定常数值计算, 为提高网格更新效率, 采用了不同变形尺度网格分类处理的思路。随后, 研究了旋翼转动对物伞系统流场分布和旋涡演变的非定常影响, 分析了不同旋翼转速下流场特性参数和降落伞气动特性的变化情况, 为进一步探究旋翼机与降落伞的气动干扰提供一定参考。

1 计算方法 1.1 数值方法及湍流模型

本文基于有限体积法采用MUSCL(monotone Upstream-centered schemes for conservation laws)三阶格式进行流场方程离散, 为提高瞬态计算精度并加速收敛采用了PISO算法, 压力插值选择Standard格式, 为有效处理瞬态计算伪扩散问题选择Green-gauss Node-based格式梯度插值。边界条件采用速度进口和压力出口, 远场为对称边界条件, 伞衣面以及旋翼机(包括前体和旋翼)为无滑移边界条件, 旋转区域Ω_r和自由流区域Ω_f的交界面Γ_interface为Ineterface边界条件。

为获得精细的流场结构、有效求解固壁绕流问题, 本文计算采用了标准壁面函数法(standard wall functions), 湍流模型选择Realizable k-ε两方程模型, 其控制方程为

$ \left\{ \begin{array}{l} \frac{{\partial \left( {\rho k} \right)}}{{\partial t}} + \frac{{\partial \left( {\rho k{u_i}} \right)}}{{\partial {x_i}}} = \frac{\partial }{{\partial {x_j}}}\left[ {\left( {\mu + \frac{{{\mu _t}}}{{{\sigma _k}}}} \right)\frac{{\partial k}}{{\partial {x_j}}}} \right] + \\ {G_k}-\rho \varepsilon \\ \frac{{\partial \left( {\rho \varepsilon } \right)}}{{\partial t}} + \frac{{\partial \left( {\rho \varepsilon {u_i}} \right)}}{{\partial {x_i}}} = \frac{\partial }{{\partial {x_j}}}\left[ {\left( {\mu + \frac{{{\mu _t}}}{{{\sigma _\varepsilon }}}} \right)\frac{{\partial \varepsilon }}{{\partial {x_j}}}} \right] + \\ \rho {C_1}E\varepsilon -\rho {C_2}\frac{{{\varepsilon ^2}}}{{k + \sqrt {\nu \varepsilon } }} \end{array} \right. $

(1)

式中：经验系数${\sigma _k} = 1.0，{\sigma _\varepsilon } = 1.2，{C_2} = 1.9；\eta = {\left( {2{E_{ij}}\cdot{E_{ij}}} \right)^{\frac{1}{2}}}\frac{k}{\varepsilon }；{E_{ij}} = \frac{1}{2}\left( {\frac{{\partial {u_i}}}{{\partial {x_j}}} + \frac{{\partial {u_j}}}{{\partial {x_i}}}} \right)；{u_t} = \rho {C_\mu }\frac{{{k^2}}}{\varepsilon }；{C_\mu } = \frac{1}{{{A_0} + {A_s}U\cdot\left( {k/\varepsilon } \right)}}；{C_1} = {\rm{max}}\left( {0.43, \frac{\eta }{{\eta + 5}}} \right)$。

1.2 网格更新模型

迭代格式的时间推进方法在每个时间步内都要对流场网格进行更新, 以确保网格变形后能够继续满足计算要求。为了提高网格更新效率, 本文将Diffusion Smoothing和Remeshing两种网格更新方法相结合, 对不同变形尺度的网格进行分类处理：对小幅运动的网格采用Diffusion Smoothing^[23]方法通过求解扩散方程得到位移后网格节点的位置, 但该方法未改变网格拓扑结构, 在大变形区域容易产生负体积, 此时引入Remeshing^[23-24]更新方法进行弥补, 对超过设定的歪斜率或尺寸阈值的网格进行标记, 并局部重新划分这些网格单元, 技术途径如图 1所示。

Diffusion Smoothing方法通过求解以下控制方程来获得网格节点的位置：

$ \left\{ \begin{array}{l} \nabla \cdot \left( {\gamma \;\nabla \mathop u\limits^ \to } \right) = 0\\ \mathop {{x_2}}\limits^ \to = \mathop {{x_1}}\limits^ \to + \mathop {u}\limits^ \to \Delta t \end{array} \right. $

(2)

式中：u为网格运动速度, γ=(1/d^α)为扩散系数(α取1.5), d为正则化后的网格节点与边界的距离, x₂为下一个时间步网格节点的位置。

2 研究对象及数值模型 2.1 研究对象及计算工况

本文的研究对象由三桨叶圆锥形旋翼机(桨叶角为0°)和半球形降落伞组成, 其外形及几何参数关系如图 2所示。不考虑伞衣的织物透气性和伞绳对流场的影响, 物伞系统结构参数、来流及计算工况见表 1、2。

2.2 数值模型

本文建立直径6D_t、高9D_t(D_t为伞衣投影直径)的圆柱形流场模型, 如图 3所示。流场计算域分为两个部分：旋转流域Ω_r和自由流域Ω_f。旋翼位于旋转流域, 其尺寸为直径2.2l、高0.15l的圆柱形流场, 其他部分均处于自由流域, 两个流域交界面为Γ_interface。

旋转流域在每个时间步内均涉及网格的移动和更新, 分别对旋转流域和自由流域划分网格, 并将两流域交界面Γ_interface上的网格节点(如Q₁、Q₂和Q₃)进行连接, 保证了两个流域边界信息传递的连续性, 如图 4所示。计算采用了适应性良好的三角形网格, 为了增加计算精度和捕捉流场细节, 对旋翼机和降落伞大梯度复杂流域的网格进行加密, 网格总数为118万。

2.3 网格适应性验证

本文首先对旋翼转速ω=0 rad/s时物伞系统的绕流流场进行了数值计算, 稳定时伞衣的平均阻力系数为0.71, 与文献[25]中半球形伞阻力系数0.75的相对误差仅为5.33%, 说明本文计算方法具有较好的准确性。

随后, 对旋翼转速ω=50 rad/s时采用了4种不同网格密度的流场进行数值计算, 其稳定后的平均阻力系数C_D见表 3。从表 3可以看出：随网格数增加, 降落伞平均阻力系数C_D增加, 当网格数达到118万时, 平均阻力系数增加延缓, 但是计算时长严重增加, 因此综合考虑, 本文采用118万的网格模型进行数值计算。

3 物伞系统非定常尾流特征 3.1 流场分布规律

本文选取旋翼转速ω=0 rad/s和50 rad/s两种工况下的流场分布和旋涡演变进行分析。图 5为计算稳定(t=15 s)旋翼机- 降落伞系统的流场分布, 可以看出：高速转动的桨叶直接改变了旋翼机的尾流结构, 将其分割成两个扰动区域, 旋涡中心由Q分裂为Q₁和Q₂, 尾流低压区P₁也以桨叶平面P为分界面形成两个低压区; 桨叶运动迫使旋翼机尾流向伞衣入口延伸, 旋涡中心Q₁较Q向上移动; 旋翼机尾流低压区P₁长度增加, 伞衣入口正压区P₂长度减小且呈不对称分布; 旋翼机尾流的低动压区V完全移至旋翼上方, 缩短了到伞衣入口的距离。

3.2 尾涡演变分析

图 6为计算稳定阶段旋翼机- 降落伞系统的尾涡变化情况, 可以看出：旋翼机和伞衣表面时刻出现旋涡的生成与脱离, 旋翼转动导致前体尾流区的涡量呈不对称分布, 干扰了前体尾涡的形成, 前体尾流中的旋涡数量明显减少。同时, 前体尾流的负涡量区在旋翼扰动下逐渐后移, 与伞衣入口的负涡量区相连通, 促进了伞衣尾涡的脱离, 伞衣尾流中的旋涡数量明显增多。

图 7为物伞系统涡量Q的等值面分布, 从图 7(a)可以看出, 由于旋翼机前体和伞衣的钝体效应^[20], 来流绕流物伞表面后形成表面涡流区VS-1和VS-2, 旋翼转动扰乱了旋翼机表面涡流区VS-1的涡量分布, 形成旋转涡流区VS-1R和表面涡流区VS-1S。旋翼机表面涡流区VS-1继续向后发展, 生成脱落涡流区VB进入流场, 如图 7(b)所示, 而旋翼转动时的脱落涡流区VB明显小于无旋翼转动, 旋翼扰乱了前体表面涡流分布后形成旋转涡流区, 将前体尾涡“打碎”, 导致前体尾流中的脱落涡流区范围变小, 进入伞衣的尾涡涡流随之减少, 伞衣表面涡流区VS-2范围变小。

4 旋翼机- 降落伞气动干扰分析 4.1 流场特性参数分析

图 8为不同旋翼转速下伞衣入口处旋涡强度的径向分布, 可以看出：由于分离流动, 旋涡强度在伞衣裙边(R=±2.5 m)处的变化最为剧烈, 此时曲线峰值较高且变化梯度较大, 在伞衣内部变化则较为平缓, 随旋翼转速增加, 旋涡强度曲线的峰值减小, 径向方向的旋涡强度逐渐减弱。

图 9为伞衣入口处旋涡强度随时间的变化, 可以看出：由于旋翼机尾涡的周期性影响, 伞衣入口处的旋涡强度随时间发生着波动, 随旋翼转速增加, 旋翼对前体尾涡的扰乱作用加强, 加速了涡流区的黏性耗散, 进入流场的尾涡涡量减少, 伞衣入口处的旋涡强度逐渐减弱。

为探究旋翼转速对旋翼机和伞衣周围流场压力的影响情况, 在对称轴旋翼机后方0.5 m处(A点)和伞衣前方0.5 m处(B点)进行流场压力监测, 得到两特征位置的压力变化情况(如图 10所示), 从图 10中可以看出：旋翼运动扰乱了旋翼机的尾流负压区, 随旋翼转速增加, 旋翼后方的流场压力逐渐增加, 而伞衣入口的流场压力逐渐减小, 当ω≥100 rad/s时, 旋翼机后方和伞衣入口前方的流场压力变化逐渐延缓。

4.2 降落伞气动特性分析

图 11为伞衣沿子午线上的内外压力及压强系数分布, 可以看出：旋翼转动对伞衣面压强的影响主要在伞衣内侧, 随着旋翼转速增加, 伞衣面外侧压力几乎不变, 而内侧压力和压强系数均逐渐减小, 内外压差减小, 伞衣气动阻力减小。

图 12为伞衣阻力系数C_D随时间的变化, 图 13为平均阻力系数C_D随旋翼转速的变化, 可以看出：旋翼转动会导致伞衣阻力损失, 随着旋翼转速增加, 阻力系数逐渐减小, 但流场达到稳定的时间的更短, 当ω≥100 rad/s时, 平均阻力系数减小趋势渐缓。

5 结论

1) 本文建立的数值方法对旋翼机/降落伞非定常复合流场的计算具有较好的准确性, 同时能精确捕捉流场的尾涡细节变化。

2) 旋翼转动将前体尾流分割成两个扰动区域, 前体尾流区长度增加, 对降落伞的影响逐渐增强, 导致伞衣入口的流场结构呈不对称分布。

3) 旋翼机尾部负涡量区受旋翼影响逐渐后移, 与伞衣入口处的负涡量区相连通, 促进了降落伞尾涡的脱离, 伞衣尾流中的旋涡数量明显增加。

4) 旋翼转动扰乱了前体表面的涡流分布, 在旋翼周围形成旋转涡流区, 前体尾流脱落涡流区范围变小, 进入伞衣的旋涡强度逐渐减弱.旋涡强度在伞衣裙边处变化最为剧烈, 在伞衣入口内侧相对平缓。

5) 随旋翼转速增加, 伞衣入口正压区长度减小, 伞衣外侧压力几乎不变, 内侧压力逐渐减小, 降落伞的阻力系数随之减小, 当ω≥100 rad/s时, 平均阻力系数和压强系数减小趋势逐渐延缓。