0 引　言

水下滑翔机是一种具有长航时、大航程特点的无人水下潜航器，其通过调节净浮力大小及方向可周期性的在上滑、下滑状态中转换，以锯齿状的运动轨迹在较小的垂向空间范围内实现了极大的水平位移，在海洋资源调查、海水团垂直剖面特征提取等领域具有较大的实用价值。出于权衡滑翔阻力以及航行效率的目的，水下滑翔机最大航速一般不超过1 kn，如华盛顿应用物理研究所（UWAPL）的Seaglider、Scripp海洋研究所的Spray、Webb研究所的Slocum^{[1 - 2]}，这也导致当水下滑翔机以岸基或船基为部署方式时，其机动至远场作业区的时间过长，不利于对于突发性任务的快速响应。

与水下滑翔机航速过低的特点相反，由于空气密度约为水密度的1/800，为了维持滑翔器在空中稳态滑翔时受力平衡，滑翔器在空中进行滑翔时水平速度要远大于水下滑翔时的水平速度。黄天鹏等^[3]基于能量管理理论，对某飞翼飞行器的无动力迫降过程进行研究，结果显示，该飞行器在迫降的主要过程中表速均高于75 m/s，仅在初始阶段存在一段由65～75 m/s的速度增长区间，相较前文所述不超过1 kn的水下滑翔航速有明显优势。因此，可兼容空中、水下滑翔功能的跨介质滑翔器可以有效地弥补水下滑翔机投送至远端任务区耗时过长的缺点。

长期以来，受介质过渡过程中的入水冲击问题的限制，跨介质滑翔器的工程化进度缓慢，需得到入水冲击机理、入水策略、减速降载手段的共同推动。在入水冲击机理以及入水策略领域，云忠等^[4]从仿生学的角度出发，结合翠鸟入水过程，利用CFD软件对其设计的翠鸟水空跨介质航行器进行了入水仿真分析，得出结论：入水角度为45°～60°，入水冲击对航行器的结构要求最低。吴正阳^[5]进行了基于翠鸟入水策略的跨介质飞行器构型仿生设计及入水性能研究，揭示了翠鸟入水策略控制机理、翠鸟头颈入水缓冲机理以及入水初始条件对入水特性的影响规律。其指出，模型轴向冲击加速度峰值与初始入水速度呈明显二次关系，与初始入水角度呈线性关系。径向冲击加速度极小值与初始入水速度的关系并不显著，但随初始入水角速的增加显著增大。这也在理论上说明，通过降低入水速度和选择合适的入水姿态，可以有效减小航行器在介质过渡时所受的冲击载荷。

在减速降载手段研究中，已经形成多种基于不同技术思路的工程化方案，既有采用尾挂减速伞在介质过渡前减速的被动减速方案^{[6 - 7]}，也有使用头部喷气反冲的主动减速方案^{[8 - 10]}，还有在头部安装气囊^[11]或单层、多层缓冲头帽^{[12 - 14]}的分流载荷方案。上述方案在实践中已能有效降低航行器主体在介质过渡时所受的冲击载荷，也就使得跨介质滑翔器能克服入水冲击问题对其工程化带来的限制。

综上所述，将水下滑翔机进一步拓展功能发展而来的跨介质滑翔机，在水下滑翔以及介质过渡领域已经有了较为完备的理论研究以及工程实践。因此，出于使跨介质滑翔器在原有水下滑翔能力的基础上兼容空中滑翔能力的目的，本文提出一种使用亚音速翼型的矩形可断机翼，随后采用计算流体力学（CFD）数值模拟方法，对在0.8 Ma初始条件下采用该机翼的一种跨介质滑翔器在空中滑翔时的气动特性进行研究，为跨介质滑翔器后续的具体工程化设计与跨介质性能优化提供基础。

1 机翼设计与计算网格划分 1.1 跨介质滑翔器机翼设计

受到无尾翼且低速舵效不足的影响，飞翼布局飞行器的低速操纵特性不佳^{[15 - 16]}，因此对于滑翔速度不高且无推力的跨介质滑翔器而言，采用带翼旋转体布局比飞翼布局更有利。

在飞行器设计领域的大量理论研究以及工程实践表明，采用带翼旋转体型布局的航行器在实际飞行时的升力绝大部分由机翼产生，而机翼的气动特性又主要受翼型影响，徐世勋等^[17]研究了翼型对于水下滑翔机性能的影响。其认为，在正攻角范围内，正弯度翼型升阻比随弯度增加而变大，可以提高水下滑翔机的滑翔经济性，但在负攻角范围则相反；采用非对称翼型会降低水下滑翔机的稳定性裕度，在水下滑翔机周期性的上滑、下滑过程中可能出现失稳现象。因此，出于不过度损失跨介质滑翔器在水下滑翔时的经济性，同时降低水下滑翔时出现失稳现象风险的目的，跨介质滑翔器的机翼采取对称翼型更为有利。

跨介质滑翔器需要借助母载具赋予一定高度与速度才能开始空中滑翔，为此暂选取常见大型飞机的0.8 Ma、10000 m高度巡航状态为跨介质滑翔器空中滑翔的初始状态，此时因滑翔速度较高，需要选取最小阻力系数较小、临界马赫数较大的翼型，而NACA系列翼型中的NACA64A-010翼型（见图1）可以满足上述要求。

机翼的设计要素除了翼型外，还包括翼面积、机翼形状等，在初始设计阶段，机翼各设计要素、升力等指标间关系可表示为^[18] ：

$ \frac{W}{S} = \frac{1}{2}\rho V_{_\infty }^2{C_{{L}}}。$

(1)

式中：$ W $为航行器完全浸润于介质时所对应的净浮力，当介质为空气时，$ W $可直接取重力；$ S $为机翼面积，即机翼在水平基准面上的投影面积；$ \rho $为工作背景介质密度；$ {V_\infty } $为相对远场的航行器速度；$ {C_{{L}}} $为三维机翼升力系数，可被定义为

$ {C_{{L}}}{\text{ = }}k{C_1}{\cos ^2}{\chi _0} 。$

(2)

式中：$k$为三维修正因子，大小受机翼展弦比和平面形状影响，初次设计可选$k = 1$；${C_1}$为二维翼型升力系数；${\chi _0}$为机翼前缘后掠角，由于选取的初始条件为0.8 Ma、10000 m高度，而滑翔时的速度总体上呈缓慢降低态势，因此暂不对可能出现的短时间机翼局部跨音速现象做针对性的处理，机翼形状选择矩形，即后掠角为0°。

由于空气与水的密度相差近3个数量级，跨介质滑翔器在空中滑翔时几乎不考虑浮力带来的影响，动力来源于自身重力，而水下滑翔时的动力则完全来源于净浮力，因此，同一滑翔器在水下滑翔时相较空中滑翔时所需的翼面积要小得多，通过断点的选择，以控制展长的形式保留水下滑翔所需的机翼。

基于上述分析，本文所设计的机翼主要由2片可向后收折的机翼组成。单片翼面长2.5 m，翼型弦长0.391 m，最大厚度约0.04 m，断点位置待水下滑翔功能具体参数确定后，再进行具体设计。本文所采用的滑翔器主体为水滴型回转体线型，其由瑞典的NYSTYROM首先提出，由可调指数的椭圆以及抛物线组成^[19]。滑翔器主体直径为530 mm，长径比约为8.6，全质量800 kg。空中滑翔时，质心位于滑翔器主体相对位置0.5处，距头部约2.3 m。通过机翼的收折、展开、断开等动作，分别满足空中携行、空中滑翔、介质过渡以及水下滑翔4种状态的需求，如图2所示。

1.2 计算网格划分

数值模拟实验采用Ansys CFD系列软件，主要包括前处理器ICEM、求解器Fluent等，版本为2019R2。目前采用非结构网格建立流场计算域，并对跨介质滑翔器表面局部细节进行网格加密。网格详细信息见表1。

如图3所示，设置远场区域距离跨介质滑翔器模型的尺度约为模型长度的10倍。其中，物面边界条件设置为粘性、绝热物面，同时选择物面函数进行物面积分，以更有效地模拟边界层；远场选择基于黎曼不变量的进出流条件，并设置吸收边界避免激波反射问题。

跨介质滑翔器在空中滑翔时总体上处于减速状态，在不做额外控制的情况下，当滑翔器飞行速度低于失速速度后，就无法维持稳定的滑翔状态。此临界值速度可以通过将式(1)中的C_L替换为C_LMAX，从而反求得该滑翔器失速速度V_S=85.1 m/s。对于民机，从冗余性考量，适航条例规定在实际使用中采取平飞失速速度V_S-1g=0.94V_S作为基准速度^[20]，此时有V_S-1g=90.5 m/s，即该跨介质滑翔器空中滑翔时最小速度不应低于90.5 m/s，相较水下滑翔时一般不超过2 kn（1 m/s）体现出了明显优势。

基于此，本次数值模拟实验需计算0.8、0.6、0.4、0.2 Ma，攻角分别为−16°、−14°、−12°、−10°、−8°、−6°、−4°、−2°、0°、+2°、+4°、+6°、+8°、+10°、+12°、+14°、+16°、+18°、+20°、+22°、+24°状态下跨介质滑翔器的气动特性数据。

1.3 网格无关性验证

在来流马赫数0.8，雷诺数2.6×10⁶（基于参考长度0.391 m）条件下，分别使用610万、904万、1200万、1523万网格数量对跨介质滑翔器以+8°、+24°攻角进行空中滑翔的升力系数进行数值计算，以验证网格无关性，结果如表2所示。

可知，随网格数量的增加，+8°、+24°攻角的升力系数有单调变化的趋势，网格收敛性较好。特别是+8°攻角时，升力系数随网格数量增加向5.27接近，且变化率逐渐减小；在+24°大攻角时，升力系数变化趋势与8°对应的升力系数变化趋势接近，这说明当网格增加至一定数量后，其对滑翔器气动数据计算结果的影响已经较小。

进而将+8°、+24°升力系数随网格数量增加的变化率进行对比。可见，尽管+24°攻角升力系数变化趋势与8°对应的升力系数变化趋势接近，但其变化率相对更大。这主要是由于大攻角时，更为精细的网格有利于再现机翼尾部三维分离等复杂流动特性^[21]。因此，综合考量计算资源以及计算精度，后续均采用1200万数量的网格进行数值计算。

2 计算方法 2.1 流动控制方程

由于速度初始条件为0.8 Ma，因此在采用CFD方法对跨介质滑翔器气动特性进行计算时，需要考虑空气被压缩后所造成的影响，对应的积分形式三维可压缩非定常流场Navier-Strokes(N-S)方程^[22]为：

$ \frac{\partial }{{\partial t}}\int\limits_V {W{\text{d}}V + } \oint_{\partial V} {\left[ {{F_{{e}}} - {F_v}} \right]{\mathrm{d}}A{\text{ = }}} 0 。$

(3)

式中：$V$为控制体的体积；$\partial V$为包围控制体的封闭表面；$A$为控制体表面的面积；$W$、${F_{{e}}}$、${F_{{v}}}$定义为：

$ W = \left\{ {\begin{array}{*{20}{c}} \rho \\ {\rho u} \\ {\rho v} \\ {\rho w} \\ {\rho E} \end{array}} \right\} ，$

(4)

$ {F_{{e}}} = \left\{ {\begin{array}{*{20}{c}} {\rho V} \\ {\rho V + {n_x}p} \\ {\rho V + {n_y}p} \\ {\rho V + {n_z}p} \\ {\rho HV} \end{array}} \right\} ，$

(5)

$ {F_v} = \left\{ {\begin{array}{*{20}{c}} 0 \\ {{n_x}{\tau _{xx}} + {n_y}{\tau _{xy}} + {n_z}{\tau _{xz}}} \\ {{n_x}{\tau _{yx}} + {n_y}{\tau _{yy}} + {n_z}{\tau _{yz}}} \\ {{n_x}{\tau _{zx}} + {n_y}{\tau _{zy}} + {n_z}{\tau _{zz}}} \\ {{n_x}{\Theta _x} + {n_y}{\Theta _y} + {n_z}{\Theta _z}} \end{array}} \right\} 。$

(6)

式中：$\rho $为密度；$v = \left( {u,v,w} \right)$为三维空间速度；$E = e + {{{{\left| v \right|}^2}} / 2} = e + {{( {{u^2} + {v^2} + {w^2}} )} / 2}$为单位质量总能量；$e$为单位质量内能；$( {{n_x},{n_y},{n_z}} )$为控制体面的单位法向量；$V{\text{ = }}v \cdot n = {{n}_x}u + {n_y}v + {n_z}w$为速度在控制面单位法向投影；$p$为压力；$H$为总焓。

2.2 湍流模型选择

流场求解器采用雷诺平均N-S方程为流动控制方程，并使用有限体积法对控制方程进行离散，采用k-epsilon湍流模型，对应的输运方程为：

$ \begin{split}\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}}}{{{\delta _k}}}} \right)\frac{{\partial k}}{{\partial {x_j}}}} \right]+ \\ &{G_k} + {G_b} - \rho \varepsilon - {Y_M} + {S_k}，\end{split}$

(7)

$ \begin{split}\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}}}{{{\delta _\varepsilon }}}} \right)\frac{{\partial \varepsilon }}{{\partial {x_j}}}} \right]+\\ &{C_{1\varepsilon }}\frac{\varepsilon }{k}\left( {{G_k} + {C_{3\varepsilon }}{G_{{b}}}} \right) - {C_{2\varepsilon }}\rho \frac{{{\varepsilon ^2}}}{k} + {S_\varepsilon } 。\end{split}$

(8)

式中：${G_k}$为由于平均速度梯度引起的湍动能k的产生项；${G_{{b}}}$为由浮力引起的湍动能k的产生项；${Y_{M}}$为可压缩湍流中的脉动扩张项；$ {C_{1\varepsilon }} $、$ {C_{2\varepsilon }} $、$ {C_{3\varepsilon }} $为经验常数；$ {\delta _k} $、$ {\delta _\varepsilon } $分别为与湍动能k和耗散率$\varepsilon $对应的Prandtl数；$ {S_k} $、$ {S_\varepsilon } $为源项。

3 计算结果与分析

本次数值模拟实验均采用戴尔T7920工作站6248 R版进行计算，为24核心48线程128 G内存配置，CPU主频3.0 GHz。在1200万网格条件下，单工况平均用时15 h。

3.1 典型流场状态

典型流场状态图如下，其中坐标系的原点位于跨介质滑翔器头部顶点，x方向沿轴线指向机尾，y方向在对称平面内垂直于x方向指向上，z方向遵循右手定则。

图4和图5分别为马赫数0.8、攻角0°时的流场压力分布云图和流场马赫数分布云图。

图6为马赫数0.8、攻角0°时跨介质滑翔器表面的压力分布云图。可以看到跨介质滑翔器头部和机翼、尾舵前缘均存在较大的压力，上表面凸起处压力略大，机翼上下表面均存在压力较小的区域，且上下表面压力分布相似。

图7为马赫数0.8 Ma、攻角24°时跨介质滑翔器表面的压力分布云图。图中滑翔器上表面在头部、机翼及机翼后方段明显压力较小；滑翔器下表面在头部、机翼前缘及尾舵前缘压力较大；滑翔器段在下表面的压力普遍高于上表面；总体来说，滑翔器下表面压力明显大于上表面，迎角对滑翔器上下表面压力影响显著。

3.2 升阻特性分析

对升阻特性等气动数据分析时，采用原点位于跨介质滑翔器顶点的滑翔器坐标系，即x方向沿滑翔器纵向指向前，y方向在滑翔器对称平面内垂直于x方向指向上，z方向遵循右手定则。

图8和图9分别为跨介质滑翔器在不同马赫数下的升力、阻力特性曲线。可知，跨介质滑翔器升力系数在–5°以下的负攻角阶段变化较小，在–5°～+5°之间呈线性增加趋势，在5°～15°之间进入平台期无明显变化，其后继续增加；阻力系数呈先减小后增加的趋势，其最小值位于0°攻角处，这主要是因为–5°～+5°攻角时，来流绕机翼运动较为稳定，此时来流对翼面的作用力随攻角变化较为平稳，使得升力可以随攻角增加线性变化，而摩擦阻力、压差阻力以及诱导阻力变化较小。当攻角继续增加至5°～15°时，来流流经滑翔翼尾部后开始出现三维分离现象，并且随着攻角增加，三维分离现象更加明显，此时滑翔器升力变化不大，但压差阻力迅速上升，就造成了阻力系数迅速增长。

由图还可知，升力、阻力特性在马赫数0.8的结果相对于其他结果曲线有一定程度的偏离，这是因为在滑翔器飞行马赫数达到0.8时，机翼上表面已经出现超音速区，此时滑翔器整体处于跨音速状态，其气动特性随攻角变化趋势更明显。

由图10可知，滑翔器升阻比曲线关于攻角0°近似呈中心对称，在攻角4°～6°时取得最大升阻比5.5～7.5，且在马赫数为0.8时接近6，马赫数小于0.8时均大于6，马赫数0.6时得到升阻比的最大值；最大负升阻比对应攻角约在–5°左右，马赫数0.6时得到负升阻比的最大值。

3.3 气动力与力矩系数分析

图11和图12分别为跨介质滑翔器轴向力、法向力系数。其中，轴向力始终为负，由机头指向机尾。

由图11可知，轴向力系数曲线相对0°攻角具有一定的对称性，总体呈“M”形。主要因为该跨介质滑翔器采用带翼旋转体构型，关于水平面对称性较好，因此轴向力系数在攻角绝对值相同时数值接近。当攻角为0°时，法向力系数与升力系数接近，轴向力系数与阻力系数相近。当攻角在0°～7°，因升力系数相较阻力系数增长更快，使得合力方向不断向升力方向偏转，导致虽然合力大小有所增长，但其在滑翔器轴向上的分力，即轴向力的大小不增反降，在图像上表现为各个速度背景下轴向力系数曲线在6°～7°附近到达极小值点。当攻角由7°进一步增大后，此时攻角增加对于升力的增益减小，而对阻力的增益变大，此时滑翔器所受合力不仅增大，且向阻力方向偏转，使得轴向力迅速增大，在图像上表现为，轴向力系数的绝对值随攻角增加而迅速增长。

由图12可知，法向力系数曲线也关于攻角0°呈近似中心对称，在攻角正半轴与负半轴范围内，法向力系数绝对值随攻角绝对值增加而增加，且在攻角小于–10°以及约5°～12°范围内存在变化速率较小的平台区，这也与一般的认知规律相符合。

由图13可知，俯仰力矩系数随攻角的增加逐渐减小，在攻角约–10°～10°的范围内减小速率较大。这主要是因为随着攻角由–10°增长至10°，各翼面相对质心产生的俯仰力臂先不断减小，而后反向增加，且在–10°～10°范围内变化最为明显。

3.4 静稳定度分析

由攻角所引起的升力在纵轴上的作用点为焦点。由图14可知，在攻角为–11°和+11°附近分别达到焦点的极小值，其中马赫数为0.8时焦点系数的变化值最小，焦点最大值约为0.6。在正负8°～12°的攻角范围内，焦点系数变化较大，这主要在于在攻角附近时，来流在滑翔翼尾部发生了分离，导致滑翔器整体气动特性发生非线性变化。

由于跨介质滑翔器不具备动力推进系统，空中滑翔速度不高，使其气动控制效能不高。并由于空中滑翔工况的需求，滑翔器在空中滑翔时也不需做出各种机动动作，这就使得滑翔器应具备也可以具备较好的静稳定度，特别应具备较好的纵向稳定性。一般而言，当滑翔器焦点位于质心（相对位置0.5）之后时，即认为其具备一定的纵向稳定性，由图易知，在攻角–8°～+8°范围内，焦点系数均大于0.5，并且由于滑翔经济性要求，滑翔器一般倾向于将稳态滑翔姿态相关参数维持在升阻比最大区间内，对应在本文中即攻角4°～6°处，此时滑翔器飞行焦点系数变化较小，基本维持在0.58附近，可知本设计满足稳定性需求。

4 结　语

针对现阶段水下滑翔机水下航速慢，以岸基或船基部署时转场至预定作业区耗时过长的问题，提出一种采用亚音速对称翼型的矩形可断机翼，通过断翼的方式，使水下滑翔机兼容了水下滑翔与空中滑翔功能，具备了跨介质滑翔能力。使用CFD方法对相应条件下的使用该机翼的带翼旋转体型跨介质滑翔器气动特性进行了数值仿真计算，主要有以下结论：

使用亚音速对称翼型、矩形可断机翼的跨介质滑翔器失速速度为85.1 m/s，从冗余性考量，实际滑翔速度不应低于平飞失速速度90.5 m/s。

总体气动特性方面，使用该机翼系统的跨介质滑翔器升力特性、阻力特性、俯仰特性以及稳定性特性随攻角曲线的变化趋势受马赫数影响小。

升阻比特性方面，滑翔器升阻比随攻角增加呈先减后增趋势，且升力、升阻比特性曲线关于零升阻比、零攻角点近似呈中心对称，正升阻比极值点位于6°附近，约为7.5；负升阻比极值点位于–5°附近，约为–8。

俯仰特性方面，滑翔器的俯仰力矩在–10°～10°攻角范围内变化最为明显，且与攻角近似呈线性变化关系。

静稳定度方面，虽然焦点系数总体上受马赫数、攻角影响较大，但在攻角–8°～+8°范围内，焦点系数均大于0.5，且在稳态滑翔攻角4°～6°范围内约为0.58，说明本文所设计机翼可使得跨介质滑翔器在典型空中滑翔状态下具备较好的纵向稳定性。

综上，采用亚音速对称翼型、矩形可断机翼可以使得跨介质滑翔器在不干扰水下滑翔功能的情况下具备空中滑翔能力，可以为跨介质滑翔器的工程化以及气动、水动特性综合优化提供基础。