潜艇应急上浮机动不稳定性建模与仿真*

吕帮俊，陈 佳，黄 斌，彭利坤

（海军工程大学动力工程学院，武汉 430033）

0 引言

潜艇在水下遭遇火灾或舱室进水，需要快速上浮至水面状态。一般的操纵方法是，用尾升降舵造成一定的尾倾，增加螺旋桨的转速，同时吹除主压载水舱使潜艇快速上浮到水面。如果潜艇上浮时跃出水面，浮心将会迅速下移，此时上层建筑内的水尚来不及排出导致重心下移滞后，因此，将会造成潜艇的短时失稳，倘若潜艇在浮至水面状态之前已经形成了比较大的横倾，浮至水面后极易形成危险横倾［1］。Itard 通过自航模试验观察到“枯叶”现象，潜艇上浮最大横倾角可达到60°［2］。戴余良等采用奇异性与分叉理论，分析了潜艇在左右舷对称的情况下，应急上浮运动并不限于垂直面内，可能出现横滚，在某种临界条件下，该运动可能失稳而出现分叉现象［3］。Watt 等认为潜艇应急上浮过程中，造成横倾失稳，最关键的因素是指挥台围壳受到非对称流动的影响［4］。如何在潜艇设计阶段针对此强非线性机动过程的潜艇运动状态进行准确预报，是避免出现上浮横摇甚至失稳问题的核心和关键。

从潜艇的设计和安全航行的角度出发，需要能够精确预报潜艇在各种战术机动，或者应急机动时的运动状态变化情况。在应急吹除过程中，需要对极端非稳态机动过程进行计算仿真。目前有多种手段可用于潜艇应急浮起过程的运动状态预报，自航模试验方法是全尺度机动预报最有效的方法［4］。然而，进行这些测试与试验，一方面需要高额的费用，另一方面，相当一部分试验由于其高风险性，在很多国家的海军是被禁止的。此外，由于试验条件限制，很多试验达不到实艇实尺度所能达到的雷诺数，尺度效应的问题不能得到有效解决。

计算机仿真提供了成本相对较低的研究方法，主要包括两类：一类是采用计算流体力学软件，模拟实艇或缩比模型在水中的运动，如CFD 流体计算软件，可以同时输出各种流动特征参数，包括压力，速度，涡流强度，切应力，紊流程度，等［5］。分析这些流动特征可以对该问题的物理过程提供深度视角，而试验方法中，无法得到这些细节，尤其是非稳态机动过程中测量流动区域的相关参数非常困难。另一类是采用潜艇六自由度运动微分方程和水动力模型试验结果，进行潜艇运动仿真求解和操纵性能预报，该方法最早由Gertler 和Hagen 在1967 年提出［6］，并由Feldman 在1979 年进行了修正［7］。运动方程中的水动力和力矩都是潜艇运动状态参数及其微分的函数，并采用一定阶数的泰勒展开进行近似和简化。由于泰勒展开式中，各项的形式都是运动参数乘以定常的水动力系数，因此，该方法也被称为基于水动力系数的仿真方法。Feldman 修正的模型考虑了潜艇在高速大舵角回转运动中横向流的影响及水平面运动过程中的水动力非线性，运动过程中水动力产生的历程效应，但该运动模型，基于来流攻角小于18°的水动力试验数据，因此，该模型仅适用于潜艇常规操纵时的性能预报［8］。Park 等采用高阶多项式，分别对潜艇大攻角拖曳水池和风洞试验结果进行曲面拟合，取得了较高的吻合度，该模型适用于潜艇高速回转以及应急上浮大攻角机动时的运动预报［9］。刘辉等结合大攻角及螺旋桨变负荷拘束船模水动力试验，采用潜艇大攻角应急挽回操纵模型，较好地预报了潜艇舱室进水时成功挽回能力和运动规律［10］。

随着现代潜艇朝着大潜深、高航速的方向发展，尤其是战术训练中往往需要潜艇在垂直面以高航速、大攻角作强机动，上述机动过程中，水动力的非线性、不对称性和耦合特性，使得潜艇的空间机动呈现高维复杂非线性特征，也是目前潜艇操纵控制和运动预报领域研究的重点和难点。区别于传统的基于攻角和漂角定义下的潜艇操纵运动水动力描述形式［11］，本文采用了潜艇大攻角操纵运动数学模型，分别结合高压气吹除指数模型和拉瓦尔喷管模型，对潜艇高压气吹除应急上浮过程的运动状态进行了预报，并重点研究了应急上浮过程横倾角的影响因素，为制定科学安全的潜艇应急上浮策略提供解决方案。

1 潜艇大攻角操纵运动数学模型

图1 描述潜艇空间机动的坐标系

2 稳定性判据

为了给出大攻角强机动工况下潜艇的横摇稳定性判据，基于上述大攻角操纵运动数学模型中的横摇力矩方程式（4），可以得到单自由度准静态稳定性指标［4］

3 高压气吹除主压载水舱模型

高压气吹除主压载水舱数学模型是进行高压气吹除计算最便捷而有效的手段，文献［4］采用了指数模型对高压气吹除过程进行了计算，并进一步对吹除后潜艇的运动状态进行了预报。然而，指数模型没有考虑排水孔面积以及水舱外背压变化等因素的影响，计算结果存在一定误差，因此，本文将采用更为准确的拉瓦尔喷管模型［15］对高压气吹除过程进行计算。在此，分别给出指数模型和拉瓦尔喷管模型，用于计算和比较不同吹除模型对潜艇上浮状态预报结果的差异。

3.1 指数吹除模型

3.2 拉瓦尔喷管模型

瑞典国家哥德堡船模试验池Bystrom 基于拉瓦尔喷管理论，分3 部分给出了高压气吹除主压载水舱数学模型［15］，成为潜艇高压气吹除系统计算的理论基础。

压载水舱的排水模型

压载水舱的压力变化模型

气瓶释放气体流量模型

式中，Ph为压载水舱的舷外瞬时压力，ρ 为海水密度，Cn为损耗系数，Vh为排水速度，Ah为水舱排水孔面积，m˙F为高压气流量，At为喷嘴的喷口面积；Ct为阀流量系数（0≤Ct≤1），k 为等熵常数，R 为气体常数，PF0为气瓶中的初始压力，mF0为气瓶中初始质量；TF0为气瓶初始温度。

软件各模块组成如图3 所示，其中，初始化模块用于设定潜艇的初始航速、深度、纵横倾角、选取的吹除模型以及水舱吹除顺序、吹除延时等，高压气吹除模块分别按照3.1、3.2 中的相关理论计算高压气吹除过程的各参数，潜艇大攻角运动模型解算模块为整个软件的核心，采用四阶龙格库塔数值积分算法求解潜艇操纵运动参数，并将潜艇运动参数及水舱水量、高压气压力变化等状态参数存储为excel 表格，用于Matlab 调用进行作图分析。

图3 潜艇高压气吹除应急上浮机动仿真分析软件组成

4 仿真平台搭建

为了开展潜艇应急上浮过程的仿真研究，基于VC++平台，开发了仿真软件，界面如图2 所示，界面上方采用曲线形式显示潜艇深度、纵倾、横倾、冲角、首尾升降舵及方向舵角等操纵运动参数随时间的变化历程，界面下方中间部分采用数值实时显示仿真过程中潜艇运动参数变化情况，界面其余部分为航速、操舵、均衡、高压气吹除等控制输入。

图2 潜艇高压气吹除应急上浮机动仿真分析软件界面

5 潜艇应急上浮过程运动仿真

5.1 吹除模型的影响

为了研究采用不同吹除模型对潜艇应急上浮横倾角的影响，采用式（1）～式（6）给出的潜艇大攻角操纵运动数学模型进行仿真计算，同时为了验证计算结果的准确性，采用与文献［4］一致的初始条件和挽回控制手段：潜艇均衡好，处于无纵倾定深定常运动状态，初始指令航速与实际航速uc=u=3 m/s，初始深度ζ=100 m。从0 时刻开始吹除，指令航速同时提高到uc=6 m/s，尾升降舵控制潜艇抬首且纵倾保持为尾倾20°。潜艇首部和尾部4 个主压载水舱容积及容积中心纵向坐标如表1 所示。

表1 主压载水舱容积及容积中心纵向坐标

如图4 所示，实线和虚线分别为采用高压气吹除拉瓦尔喷管模型和指数模型计算得到潜艇各状态参数变化情况，虚线所示结果与文献［4］基本吻合，从而验证了本文所采用仿真计算方法的正确性。从仿真结果可以看出，随着吹除过程的进行以及潜艇的上浮，无论是横摇角速度p 还是横倾角φ都在逐渐增大，通过对比可以发现：

图4 分别采用高压气吹除指数模型和拉瓦尔喷管模型计算得到的潜艇应急上浮时运动参数变化情况

1）指数模型计算的首部的1、2 号主压载水舱尽管容积差别很大，几乎同时吹空（1 号略快于2号），尾部的3、4 号主压载水舱也有类似的结果（3号略快于4 号），这主要是因为该模型假设进入每个主压载水舱的高压空气质量与该水舱的体积成正比，而不考虑水舱排水孔面积、本身容积大小，以及外部背压对进入该水舱高压气流量的影响等因素，影响吹除率的最主要因素为水舱的位置布置。

2）拉瓦尔喷管模型计算得到的吹空顺序为：4、2、3、1 号，尾部的4 号主压载水舱尽管由于艇体尾倾导致排水背压大，但此压载水舱容积最小（仅30 m3），所以最先吹空。2、3 号主压载水舱虽然容积相同，但尾倾使得2 号主压载水舱排水背压较小，所以第2排空，3 号主压载水舱第3 排空。1 号主压载水舱虽然在上浮过程中排水背压最小，但体积最大（100 m3），所以1 号主压载水舱最后排空。因此，综合考虑水舱的容积、位置布置、艇的运动状态等因素，拉瓦尔喷管模型的计算结果与实际情况更为符合。

3）在上述拉瓦尔喷管模型计算得到的吹空顺序下，潜艇上浮初期，尽管尾升降操上浮满舵使潜艇抬首，由于4 号主压载水舱容积最小，排水速度快，很快产生较大的首倾力矩，导致潜艇尾倾形成相对较慢。随着吹除过程的继续，首部主压载水舱的排水量逐渐增大，在上浮过程的中后期，吹除主压载水舱形成的合力矩逐渐由首倾变为尾倾力矩，最终的合力矩大小为640 t·m，因此，该模型计算得到上浮接近水面时的最大横倾为6.8°，而指数模型计算结果仅为2.0°。这主要是因为应急上浮过程中尾倾的形成，能将原本与负垂速（-w）方向一致的负浮力转移了一部分到沿Gx 轴方向来使u 加速，一方面使垂速绝对值增长速度放缓，也即负攻角形成的慢且绝对值将减小，另一方面Gx 轴向速度u 的增加，进一步使攻角减小。因此，尾倾形成快有助于减小应急上浮过程的横倾角。

5.2 初始深度的影响

为了进一步分析初始深度对潜艇上浮机动的影响，分别设定仿真初始深度为100 m 和150 m，同时吹除全部主压载水舱，采用拉瓦尔喷管模型计算吹除过程，潜艇状态参数变化情况如图5 所示，可以看出，如果初始深度增加到150 m（图中虚线），上浮时间将会大幅增加，一方面是因为上浮距离增大，另一方面是因为主压载水舱的吹除率明显减缓。式（7）中的静稳性项的变化也减缓，从正值减小到0 的时刻也相较于100 m（图中实线）吹除的情况推后约20 s，而应急上浮过程中横倾角迅速增大主要发生在直至上浮到水面的这一阶段。因此，从图5可以看出，分别在100 m 和150 m 吹除全部主压载水舱，直至上浮到水面状态，产生的最大横倾角较为接近。

图5 分别在100 m 和150 m 吹除主压载水舱时潜艇状态参数变化情况

5.3 水舱吹除顺序的影响

考虑到1 号主压载水舱位于首部，吹除时有利于尾倾形成进而抑制横倾增大，而采取同时吹除全部主压载水舱时，根据4.1 的仿真结果，1 号水舱将会最后吹空，没有很好地发挥其产生尾倾力矩的效用。为了在上浮过程中较快形成尾倾，可以考虑改变水舱的吹除顺序，例如，先吹除1 号主压载水舱，延时5 s 后吹除剩余的2～4 号主压载水舱。相应的仿真计算结果如图6 所示，图中实线为同时吹除，虚线为先吹除1 号主压载水舱的情况。可以看出，潜艇尾倾形成时间大大缩短，且上浮过程的最大横倾角也由6.8°减小到4.1°。结合4.2 可以得出，应急上浮过程的横倾角主要由所持续的时长决定，并受水舱吹除顺序的影响较大，而受初始深度的影响较小。

图6 延时吹除尾部主压载水舱时潜艇状态参数变化情况

6 结论

本文以潜艇大攻角水动力试验结果为基础，对基于水流入射角和方向角定义下的潜艇大攻角操纵运动数学模型进行了数值仿真求解，采用不同初始状态以及不同的高压气吹除模型，预报了潜艇应急上浮过程的运动状态，重点分析了上述因素对潜艇应急上浮过程横倾的影响，仿真结果表明：1）高压气吹除拉瓦尔喷管模型，相较于指数模型加入了水舱排水孔面积、潜艇上浮过程水舱外部背压变化等因素的影响，更能准确计算实艇高压气吹除过程，从而提高潜艇应急上浮机动过程运动状态预报的准确性；2）应急上浮过程中尾倾的形成有助于抑制大的负攻角的形成，从而有助于减小应急上浮过程的横倾角；3）应急上浮过程中的横倾角受水舱吹除顺序的影响较大，而受初始深度的影响较小。