0 引　言

水下滑翔机^[1]作为新型水下机器人，因其长航程、高续航和低噪声等优点^[2]，在军事和民用领域受到广泛重视，是军事侦察和资源开发的主要设备^[3]。为了更好的探究水下滑翔机的滑翔运动特性，需要对其滑翔运动过程进行深入研究。张华等^[4]利用CFX水动力计算软件，对垂直面内运动的水下滑翔机运动参数进行探究，在获取运动参数的基础上，结合LQR控制方法完成垂直面内的水下滑翔机运动特性的研究。顾建农等^[5]从柯基霍夫方程出发，建立水下滑翔机的非线性仿真模型，对其滑翔能力及回转特性进行深入研究。

随着计算方法和计算机技术的不断创新和发展，计算流体力学（CFD）迎来了新的飞跃，逐渐被应用于水下航行器的运动仿真^[6]。房萍萍^[7]采用动网格技术和6DOF（6 Degrees of Freedom）方法，对不同结构特征的AUV在洋流干扰下的运动进行了仿真分析，研究AUV在复杂海况下的抗干扰能能力。艾晓锋^[8]将MRF方法^[9]与动网格技术结合起来，采用多重坐标系对螺旋桨进行数值模拟，并将螺旋桨与AUV结合起来进行里数值模拟，获得了不同AUV速度下的速度场与压力场。刘金夫^[10]基于重叠网格技术^[11]和6DOF方法，分析了水下滑翔机在复杂洋流影响下的滑翔运动性能。李健欣^[12]结合动态网格技术，对水下滑翔机的静水阻力、直航以及滑翔运动进行直接模拟。

现有研究多采用升阻比作为经济性指标，未关联之字形运动的水平-垂直位移协同关系,且大攻角滑翔的内在机理尚未明确,亟需开展针对性研究。本文通过CFD仿真与机理分析，明确了重浮力差和重浮心沿机身轴线距离对水下滑翔机滑翔性能的协同影响规律，揭示了大攻角现象的力矩失衡机理，最终确定最佳下潜工况为重力增加1.2%G、重心前移0.2%L。研究结果为水下滑翔机的参数设计、运动控制及性能优化提供了定量理论依据。

1 研究对象和理论基础 1.1 滑翔机模型及其基本参数

本文研究的水下滑翔机为课题组在研的仿生型可变翼水下滑翔机，内部耐压壳采用传统圆柱耐压壳设计，为了降低阻力非耐压主体外形仿照鲭鱼设计，收放翼机构采用涡轮蜗杆。其三视图如图1所示。其基本参数和表1所示。

1.2 理论基础 1.2.1 基本控制方程和湍流模型

流体三大守恒方程分别是质量守恒方程、动量守恒方程和能量守恒方程^[13]。对于所研究的不可压缩流体，主要遵循质量守恒方程和动量守恒方程。质量守恒方程也叫连续性方程，其形式为：

$ \frac{\partial u}{\partial x}+\frac{\partial v}{\partial y}+\frac{\partial w}{\partial z}=0 。$

(1)

式中：u、v、w是速度矢量$ \mathbf{\nu } $在x、y、z方向的分量。

动量守恒方程方程，也叫Navier-Stokes方程其形式如式(2)所示：

$ \frac{\text{d}\mathbf{\nu }}{\text{d}t}={F}_{b}-\frac{1}{\rho }\nabla p+v\Delta \mathbf{\nu }+\frac{1}{3}\nu \nabla \left(\nabla \mathbf{\nu }\right) 。$

(2)

式中：$ v $为粘性运动系数。

本文采用标准k-ε两方程湍流模型，适用于高雷诺数流动，具有计算稳定性强、工程适用性广的特点，其控制方程如下：

湍动能k方程：

$ \begin{split}&\frac{\partial (\rho k)}{\partial t}+\frac{\partial (\rho {u}_{i}k)}{\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}+{G}_{b}-\rho \varepsilon 。\end{split} $

(3)

耗散率ε方程：

$ \begin{split}&\frac{\partial (\rho \varepsilon )}{\partial t}+\frac{\partial (\rho {u}_{i}\varepsilon )}{\partial {x}_{i}}=\frac{\partial }{\partial {x}_{j}}\left[\left(\mu +\frac{{\mu }_{t}}{{\sigma }_{\varepsilon }}\right)\frac{\partial \varepsilon }{\partial {x}_{j}}\right]+\\ & {C}_{1\varepsilon }\frac{\varepsilon }{k}({G}_{k}+{C}_{3\varepsilon }{G}_{b})-{C}_{2\varepsilon }\rho \frac{{\varepsilon }^{2}}{k}。\end{split} $

(4)

式中：$ {\mu }_{t}=\rho {C}_{\mu }\dfrac{{k}^{2}}{\varepsilon } $为湍流黏度；$ {G}_{k} $为平均速度梯度产生的湍动能；$ {G}_{b} $为由浮力产生的湍动能；模型常数$ {C}_{1\varepsilon } $=1.44，$ {C}_{2\varepsilon } $=1.92，$ {C}_{\mu } $=0.09，$ {\sigma }_{k} $=1.0，$ {\sigma }_{\varepsilon } $=1.3，采用Launder-Spalding标准取值。

1.2.2 受力分析

水下滑翔机在水中运动时受到重力、浮力、水动力及流体阻力的综合作用，在对其进行受力分析时，需要首先在大地坐标系下进行求解之后再转换到体坐标系。其中体坐标系oxyz与滑翔机固连，为动坐标系，惯性坐标系OXYZ固定在地面上，为静坐标系，如图2所示。受力分析见图3。

机身轴线与OX轴所成的夹角称作俯仰角$ \theta $，与速度方向所成的夹角称作攻角$ \alpha $，速度与OX轴所成的夹角称作姿态角$ \beta $。当滑翔机处于最终稳定状态滑翔时，其受力平衡方程为：

$ L\sin \alpha +(G-B)\sin \theta -D\cos \alpha =0 ，$

(5)

$ L\cos \alpha +D\sin \alpha -(G-B)\cos \theta =0，$

(6)

$ G({x}_{G}-{x}_{B})-{M}_{z}=0 。$

(7)

式中：L和D分别为滑翔机所受到的升力和阻力；M_z为俯仰力矩；B为浮力；G为重力；x_G和x_B分别为重心和浮心的x坐标。

为了研究其运动过程，采用6DOF运动方程描述水下滑翔机在水中的滑翔运动过程，其数学模型包含平动方程和旋转方程。

平动方程（质心运动）：

$ m\frac{{\mathrm{d}}\boldsymbol{V}}{{\mathrm{d}}t}=\sum\boldsymbol{F} 。$

(8)

式中：m为滑翔机质量；$ \boldsymbol{V} $为速度向量；$ \boldsymbol{F} $为所受合力。

旋转方程（绕质心转动）：

$ {\boldsymbol{I}}\frac{{\mathrm{d}}\mathbf{\omega }}{{\mathrm{d}}t}+\mathbf{\omega }\times ({\boldsymbol{I}}\mathbf{\omega })=\sum\boldsymbol{M}。$

(9)

式中：$ {\boldsymbol{I}} $为转动惯量矩阵；$ \mathbf{\omega } $为角速度向量；$ \boldsymbol{M} $为滑翔机受力矩。

1.2.3 之字形运动经济性指标定义

单次之字形运动由下潜和上浮组成，假设运动对称，那么单次之字形运动的总水平位移为：

$ {S}_{X}=\left|2\times\frac{H}{{V}_{Y}}\times{V}_{X}\right|。$

(10)

式中：$ {S}_{X} $为水平位移；H为目标深度；$ {V}_{Y} $为垂直速度；$ {V}_{X} $为水平速度。

由于目标深度H固定，$ {S}_{X} $与$ |{V}_{X}/{V}_{Y}| $成正比。因此，定义水平速度与垂直速度的比值$ |{V}_{X}/{V}_{Y}| $为滑翔经济性指标，该值越大，单位垂直位移的水平推进效率越高。

2 滑翔运动仿真分析 2.1 仿真基础 2.1.1 计算域选取及网格划分

DEMIRDŽIĆ等^[14]对中计算域的选择方法，进行多次数值计算验证，将计算域设定为70 m×40 m×8 m。考虑到水下滑翔机具有对称性结构，本文仅研究其在X、Y方向的平动以及绕Z轴的旋转运动，将计算域划分为7个子区域，具体划分方式如图4所示。其中几何对象被包裹在旋转域Fluid_r中，旋转域跟随滑翔机进行主动运动，其余的Fluid_C、Fluid_L、Fluid_R、Fluid_U、Fluid_D通过对其施加求解出的滑翔机的速度进行平动，剩下的Fluid_s为静止域。这些流体域的运动方式见表2。

网格划分采用混合网格（非结构网格加结构网格），针对不同区域特性选择适配的网格类型。

Fluid_r域（含滑翔机主体）：因几何结构复杂，采用非结构化四面体网格离散，同时，为准确捕捉滑翔机周围的流场特性，对滑翔机近壁区域进行了局部网格加密，以提高数值模拟的精度。

各区域交界面通过interface边界条件耦合，Fluid_C、Fluid_L及Fluid_U域的相邻交界面定义为interior类型，确保流场变量连续性。为验证网格无关性，生成388万、549万、776万3套网格，以1 m/s下的阻力为考核指标进行验证，结果如表3所示。采用“相对误差=|当前网格阻力-较稀疏网格阻力|/较稀疏网格阻力×100%”计算，776万网格与549万网格的相对误差仅1.29%，表明549万网格已满足计算精度要求，最终采用该网格方案，网格划分结果如图5所示。

2.1.2 计算参数设置

仿真采用非稳态数值模拟，结合标准k-ε湍流模型与VOF多相流模型（空气为第一相，水为第二相）；启用Implicit Body Force模型添加浮力，结合Open Channel Flow模型模拟近水面下潜过程。边界条件设置：入口为Pressure Inlet（流速0 m/s，模拟静水环境），出口为Pressure Outlet（相对压力0 Pa），滑翔机表面和计算域边界设为无滑移壁面。

引入两类用户自定义函数（UDF）：

1）DEFINE_SDOF_PROPERTIES：定义滑翔机6DOF被动运动特性，包括质量、转动惯量、重心位置等参数；

2）DEFINE_CG_MOTION：实时获取滑翔机质心瞬时速度，传递至Fluid_C、Fluid_R等平动域，驱动其同步运动。

数值求解采用SIMPLE算法，压力项采用PRESTO!离散格式，动量、k、ε方程均采用二阶迎风格式提高精度。时间步长设为0.005 s，总计算步数14000步（对应仿真时长70 s）。通过监测速度及俯仰角的时程曲线，确认50 s后各参数波动幅度小于5%，视为进入稳定滑翔状态，后续分析采用50～70 s的稳态平均数据。

2.2 船池拖曳试验验证

为验证网格划分方法与仿真模型的有效性，在某船池开展1∶1比例滑翔机模型的拖曳试验，测量不同速度下的阻力值，并与仿真结果进行对比。拖曳速度范围为1～6 m/s，每级速度稳定30 s后记录阻力数据，每个速度工况重复3次，取平均值作为试验结果。

仿真与试验的阻力对比结果如表4所示。由表可知：

速度≤2 m/s时，阻力误差小于1.6%（1 m/s时误差1.39%，2 m/s时误差1.54%），与试验结果高度吻合；

速度>2 m/s时，误差逐渐增大（最大为6.41%，5 m/s时），主要原因是高速下试验模型表面出现轻微流场分离，而仿真中采用的标准k-ε模型对强分离流的预测精度存在一定偏差，但误差仍在工程允许范围内（≤7%）。

综上，仿真模型与网格划分方法具有较高的可靠性，可用于后续滑翔运动特性分析。

3 滑翔运动特性分析 3.1 不同重浮力差和重浮心距离对滑翔性能的影响

为探究重浮力差（$ \Delta G=G-B $，以$ \mathrm{\% }G $表示）和重心前移距离（$ \Delta x={x}_{G}-{x}_{B} $，以$ \mathrm{\% }L $表示）对滑翔性能的影响，设计3种重浮力差（0.5%G、1.0%G、1.5%G）与3种重心前移距离（0.1%L、0.15%L、0.2%L）的正交组合，共9种工况，具体如表5所示。

9种工况的稳态性能参数（50～70 s平均值）如表6所示，图6和图7分别给出了重心前移0.2%L和重浮力差1.0%G时的运动时历曲线。结合数据与曲线分析，得出以下规律：

当重浮心距离固定时，重浮力差的影响呈现差异化特征：小重浮心距离（0.1%L）下，增大重浮力差使滑翔角、攻角显著增大，水平速度和垂直速度均上升，但滑翔经济性（1.21～1.55）较低；中等重浮心距离（0.15%L）下，重浮力差从0.5%G增至1.5%G时，攻角从4.789°跃升至34.796°，经济性从2.58降至1.39；大重浮心距离（0.2%L）下，中等重浮力差（1.0%G）对应最优经济性（2.62），重浮力差超过1.2%G后，攻角急剧增大至20°以上，经济性显著下降。

当重浮力差固定时，重浮心距离的影响同样与重浮力差水平相关：小重浮力差（0.5%G）下，增大重心前移距离使攻角从31.350°降至3.430°，经济性从1.55提升至2.04；中等重浮力差（1%G）下，大重心前移距离（0.2%L）表现出最优协同效应，经济性达2.62；大重浮力差（1.5%G）下，重心前移距离对经济性的调节作用较弱（1.21～1.77）。

值得注意的是，工况1、工况2、工况3、工况5、工况6、工况9出现大攻角现象（|α|>20°），且伴随滑翔经济性下降，其内在机理需进一步探究。

3.2 大攻角滑翔成因机理

定义攻角绝对值超过20°为大攻角工况（工况1、工况2、工况3、工况5、工况6、工况9），小于20°为小攻角工况（工况4、工况7、工况8）。为揭示大攻角形成机理，通过拟合建立俯仰力矩与速度、攻角的定量关系，结合重浮力差产生的自身力矩进行对比分析。

首先，基于滑翔机动力学特性，选取典型工况下的仿真数据（涵盖速度0.5～2.5m/s与攻角范围0～90°），采用最小二乘法拟合俯仰力矩M_z与速度U、攻角$ \alpha $的关系，拟合函数如下：

$ {{M}_{z}=\begin{cases} {U}^{2}(1.127\alpha +0.1415) , \quad\alpha \leqslant 14{^{\circ}}，\\ {U}^{2}(-0.1185{\alpha }^{2}+5.3277\alpha -31.3356), \quad \alpha \gt 14{^{\circ}}。\end{cases} }$

(11)

拟合优度$ {R}^{2}=0.968 $，表明该模型能准确描述俯仰力矩的变化规律，曲面拟合结果如图8所示。从拟合结果可以看出，在小攻角阶段俯仰力矩和攻角呈现线性关系，并且增长缓慢，当攻角过大时，俯仰力矩和攻角呈现出二次函数的关系，具体表现为俯仰力矩随之攻角的增大先增大后减小。

根据式(7)计算重浮力差产生的自身力矩M_self，将各工况的稳定速度与攻角代入式(11)预测流体动力俯仰力矩M_z^pred，结果如表7所示。分析表明：

小攻角工况：自身力矩略大于预测流体动力矩，差值在−0.35～−0.10 N·m之间，滑翔机通过姿态微调实现力矩平衡，维持小攻角稳定滑翔；

大攻角工况：预测流体动力矩大于自身力矩，差值在0.14～3.44 N·m之间，多余力矩阻止滑翔机向小攻角姿态调整，最终形成大攻角稳态。

这一机制揭示了大攻角滑翔的本质是俯仰力矩失衡，而非主动控制策略导致。

3.3 最佳滑翔工况确定

最佳滑翔工况定义为：小攻角（|α|<20°）范围内，滑翔经济性最大的工况。该定义兼顾水平推进效率（经济性指标）与运动稳定性（小攻角避免流场分离），符合水下滑翔机长航程、低能耗的核心需求。

基于3.1节结果，在重心前移0.2%L（表现出良好协同效应）条件下，将重浮力差细化为0.5%G～1.5%G（步长0.1%G），共11个细化工况，仿真时长40 s（确保达到稳定状态），重点监测攻角和滑翔经济性的变化。细化工况的仿真结果如表8所示。随着重浮力差增大：

攻角呈持续上升趋势：重浮力差≤1.2%G时，攻角均小于20°（小攻角范围）；重浮力差从1.2%G增至1.3%G时，攻角从7.095°跃增至23.188°，进入大攻角阶段；

滑翔经济性呈“先升后降”趋势：在重浮力差1.2%G时达到最大值2.42（对应攻角5.852°，水平速度0.984 m/s，垂直速度−0.406 m/s）；重浮力差超过1.2%G后，经济性从2.49降至2.07，与攻角跃变同步，验证了大攻角对经济性的恶化作用。

对该最优工况（重浮力差1.2%G、重心前移0.2%L）进行力矩验证：预测流体动力矩（8.25 N·m）与自身力矩（9.61 N·m）差值−1.36 N·m，差值较小，力矩基本平衡，确保姿态稳定。图9给出了该工况下的速度分布云图，可见滑翔机周围流场平稳，无明显分离现象，进一步验证了该工况的合理性。

需说明的是，“最佳滑翔工况”是针对本文研究对象（特定几何参数和运动范围）的优化结果。若滑翔机几何参数发生显著变化，最佳工况需通过相同方法重新标定。

4 结　语

为研究重浮力差和重浮心距离对水下滑翔机滑翔运动的影响，本文通过CFD仿真软件对多种滑翔工况进行了数值仿真分析，得出以下主要结论：

1）重浮力差与重心前移距离对滑翔机姿态和经济性具有显著协同影响：重浮力差增大使滑翔攻角单调增大、经济性降低；重心前移距离增大使滑翔攻角单调减小、经济性提高；小攻角范围内，两者的协同作用可显著提升水平推进效率。

2）大攻角滑翔的本质是俯仰力矩失衡：滑翔机俯仰力矩与攻角的关系在小攻角阶段（α<14°）近似线性，大攻角阶段（α>14°）呈二次函数关系；当流体动力俯仰力矩大于重浮力差产生的自身力矩时，滑翔机无法维持小攻角姿态，最终形成大攻角稳态，导致经济性显著恶化。

3）本文研究对象的最佳下潜工况为：重浮力差1.2%G、重心前移0.2%L。该工况下滑翔经济性达2.49，攻角7.095°（小攻角范围），流场平稳，兼顾了运动效率与稳定性，可为水下滑翔机的浮力调节系统设计、重心配置及运动控制策略优化提供定量理论依据。