Chinese Journal of Ship Research
基于体积力模型的潜艇应急上浮运动数值模拟分析
引用格式：魏可可,高霄鹏.基于体积力模型的潜艇应急上浮运动数值模拟分析[J]. 中国舰船研究, 2021, 16(2): 49–56. WEI K K, GAO X P. Numerical simulation of emergency surfacing motion of submarine based on volumetric force model[J]. Chinese Journal of Ship Research, 2021, 16(2): 49–56.
魏可可*，高霄鹏海军工程大学舰船与海洋学院，湖北武汉 430033
摘 要:［目的］为了探索潜艇应急上浮运动的规律，基于体积力的螺旋桨简化模型，开展不同航速下潜艇的应急上浮运动数值模拟。［方法］对艇体水下定深直航、上浮以及上浮出水进行数值模拟，对比分析不同航速下潜艇上浮过程中艇的螺旋桨转速与航速的匹配，以及上浮运动时间和艇体姿态的变化。［结果］结果显示，随着螺旋桨转速的增大，潜艇获得稳定航速的时间越短，在水下定深直航运动的时间也越短；无论是纵倾角还是横倾角，其第1次发生变化的时刻、产生峰值的时刻等随转速的增大呈减小趋势，同时随着转速的增大，纵倾角的来回震荡波动时间会越来越长，而横倾角的来回震荡波动时间则越来越短。［结论］相关计算方法和研究结果对潜艇应急上浮运动研究具有一定的实际参考价值。关键词：体积力；应急上浮运动；定深直航；上浮出水中图分类号: U661.33 文献标志码:A DOI：10.19693/j.issn.1673-3185.01879
Numerical simulation of emergency surfacing motion of submarine based on volumetric force model
WEI Keke*, GAO Xiaopeng
College of Naval Architecture and Ocean Engineering, Naval University of Engineering, Wuhan 430033, China
Abstract: ［Objectives ］ In order to explore the laws of the emergency surfacing motion of submarines, a simplified model of propeller volume force was used to carry out the numerical simulation of submarine emergency surfacing motion at different speeds.［Methods］Numerical simulations of the underwater fixed depth, ascent and surfacing of the hull were carried out, and submarine buoyancy under different speeds was compared and analyzed by matching the screw rotation speed with the speed of the submarine, as well as its ascent time and attitude.［Results］As the results show, the faster the rotation speed, the shorter the time for the submarine to gain steady speed and maintain fixed depth underwater. Regardless of the trim angle or heeling angle, the time when the angle changes for the first time and that when the peak occurs show a trend of reduction although the rotational speed increases. At the same time, with the increase in rotation speed, the fluctuation time of trim angle increases, while the fluctuation time of heeling angle decreases.［Conclusions］The relevant calculation methods and results of this study can provide practical references for research on submarine emergency surfacing motion. Key words: volume force；emergency surfacing motion；depth direct sailing；rise out of the water
0 引 言
潜艇的应急上浮是指当发生舱室破损进水、卡舵等事故时，采取排除部分或者全部主压载水
舱内的水，从而使潜艇迅速浮出水面的一种紧急操纵措施。为保证潜艇应急上浮的安全性和稳定性，准确预报潜艇在应急上浮过程中的运动特性十分有必要。
和正常的上浮相比，潜艇应急上浮中艇的排水方式不同且排水速率快，艇体姿态变化幅度大。常见的潜艇应急上浮运动预报方法是通过拘束模试验[1-4] 或参数辨识[5] 等方法求解水动力系数，然后再基于六自由度运动模型进行预报。戴余良等[6-7] 应用非线性系统运动稳定性与分叉理论，以及同伦延拓数值计算方法，分析了潜艇应急上浮运动的稳态响应及运动稳定性。王晓玢等[8]通过对潜艇垂直面操纵运动进行非线性建模，利用分叉与突变理论方法对潜艇上浮运动中由静态分叉和动态分叉引发的状态突变进行了分析。潜艇水动力系数的求解较复杂且有一定的难度，而上述预报方法主要是针对潜艇某自由度上的运动，并不能真实反映潜艇在六自由度完全释放情况下的应急上浮运动特性，且预报精度还有待提高。随着计算机技术的快速发展，计算流体力学在船舶运动研究中的应用越来越广泛。Carrica等[9-10] 基于自主研发的流体软件 CFDShip-Iowa V4，实现了不同工况下潜艇六自由度运动的数值模拟计算，并对其操纵性进行了较高精度的预报。Chase 等[11] 采用重叠网格技术，对 Suboff 潜艇模型的直航拖曳和自由自航等运动进行了数值模拟计算。周广礼等[12-13] 基于 RANS 方程及流体体积模型，采用整体动网格技术，对Suboff 模型的应急上浮六自由度运动及其出水运动予以了研究。廖欢欢等[14] 采用流体体积法（volume of fluid，VOF）六自由度求解器和重叠网格技术，研究了潜艇在一定潜深下的应急上浮六自由度运动，并对上浮过程中的姿态和水动力进行了分析。以上大多学者对潜艇应急上浮运动的研究都没有考虑螺旋桨以及初始航速等问题的嵌入，为此，本文将基于STAR-CCM+软件平台，采用基于体积力模型的虚拟桨盘，针对潜艇在不同航速下的应急上浮运动进行数值模拟，并对比分析潜艇在应急上浮运动中的姿态和上浮运动时间等，用以为潜艇在实际中的上浮运动操控提供一定的借鉴。
1 计算对象
本文将以美国泰勒研究中心提供的Suboff 潜艇试验模型为研究对象。国内外针对该模型开展了大量的拖曳及流场测量试验，其三维模型图如图 1所示，模型主尺度参数如表1所示。
2 数值模拟2.1 控制方程及湍流模型
应用黏流CFD 方法开展数值计算的基础是纳维−斯托克斯（Navia-Stokes，N-S）方程，然而受限于计算能力，当前还不能通过直接求解N-S 方程（direct numerical simulation，DNS）来获得潜艇周围复杂的黏性流场。将N-S方程中的湍流脉动项进行时均化处理可以极大地减小计算量，所得雷诺平均方程（ RANS ）的具体表达如式（ 1）所示。然而，N-S方程在时均化的过程中增加了含未知量时均的雷诺应力项，通过引入湍流模型，可使方程封闭。目前，该方法已在舰艇流体力学理论研究及工程实践中得到广泛应用。
式中： ρ为流体密度； µ为流体粘性； p为静压； fi为单位质量的质量力；ui，uj分别为i, j点上x方向的速
度分量；t 为时间； xi ， xj分别为i, j点上x方向的位移分量；上加线“-”表示参量时均化； −ρu′ iu′ 为雷j诺应力。开展计算时，控制方程离散均采用二阶迎风格式，并应用分离式解法（segregated flow）对离散后的方程组进行求解。压力速度耦合迭代采用SIMPLE算法，涉及非稳态计算时，设定计算过程中的时间步长为∆t且固定不变，时间迭代采用二阶形式，即
V∆V
分，其中 为体积， ∆V为控制体的体积值， u为x方向的速度分量， ϕ是广义变量； χ(n)为第t = n × ∆t时刻χ的量值；n为时间步数。
采用 RANS 方法开展潜艇水动力计算时，湍流模型的选取是确保计算精度的关键要素，众多学者已针对该问题开展了大量研究，但所得结论存在较大的不同，且湍流模型的适用性也与网格布局及操作者的数值计算经验相关。综合各类文献资料，本文选取 Wilcox k − ω与 SST k − ω湍流模型进行优选分析。Wilcox k − ω湍流模型最早由 Wilcox 于 1988年提出，后经过不断的修正，模拟精度及其适用范围均得到了提升，STAR-CCM+软件平台中就提供有最新的修正模型，具体表达形式如下：
式中：k为湍动能；ω为特殊湍动能耗散；其余参数的具体定义及参数值参见文献[15]。Menter 等[16] 结合标准k − ε与k − ω湍流模型，提出了剪切应力模型SST k − ω，式（4）给出了该模型的最新修正形式，具体推导过程及参数定义详见文献 [15]。
2.2 计算域及网格划分
计算域选取为方形域，除域尾端边界设置为压力出口外，其余边界条件均设置为速度入口，域大小及计算域与自由液面的初始相对位置如图 2所示。进行网格划分时，为能更好地捕捉两相流自由液面，对自由液面可能通过的区域进行了网格加密。此外，为准确模拟潜艇出水后艇体周围的流场，对模型首部及指挥室围壳周围网格进行了细化处理。同时，为满足网格贴体的要求，对近壁面区域进行了局部加密，以确保艇体壁面y+在30~300 之间[12]。当 y+在 30~200之间时，表征第1层网格布置于对数层，可通过假定的壁面函数对近壁面流动进行模拟。该方法不直接求解黏性底层与过渡层内的流动，适用性较强，可极大地减少近壁面网格数量。艇体表面、域剖面及近壁面域的网格如图3所示，网格总数为700 万。
2.3 体积力模型
通过对螺旋桨入流面流速进行加权积分，获得桨的进流速度，然后代入到已知螺旋桨水动力曲线中，预估得到推力及扭矩，接着，通过添加源项的方式，对设定的体积力域内的流体施加等效的力及力矩，该方法即为体积力模型。体积力模型不仅能够将力及力矩施加于潜艇上，而且能更
为准确地模拟由螺旋桨引起的潜艇尾流场变化，相较于直接力模型，其计算量大，但计算效率远优于开展的真实螺旋桨计算。为了重点分析潜艇的运动特征，提高计算效率，螺旋桨用基于体积力模型的虚拟桨代替。该虚拟桨的参数采用DTMB 4 383标准桨的，模型主参数如表2所示。
2.4 数值计算的验证
针对漂角β = 2◦ ∼ 14◦、间隔2◦的 Suboff 模型的无因次纵向水动力X ′进行数值计算，计算结果与试验值[17] 的对比如图4所示。由图4可看出，不同漂角下艇体纵向水动力数值计算结果与试验值吻合良好，证明该方法可行。因此，选取该方法进行后续的数值仿真计算。
2.5 计算工况
为探讨不同航速对潜艇上浮运动的影响规律，选取上浮力为总排水量的5%、螺旋桨转速分别为 300 ，400 ，500 ， 600 r/min的工况进行数值模拟计算，具体计算工况如表3所示。
3 不同航速下数值模拟结果分析3.1 螺旋桨转速与航速的匹配
假定潜艇上浮之前在水下做定深直航运动，为获取模型水下稳定直航状态，设定初始直航航速U=0，在螺旋桨推力的作用下，潜艇沿x方向做直航运动。为避免在初始计算阶段出现计算结果不稳定的现象，在 t=0~4 s时间段内将艇体按静态处理，待螺旋桨推力及流场计算结果稳定后，即在 t1 时刻，将 x方向运动完全释放，潜艇开始做直航运动。不同转速下螺旋桨推力及艇体阻力的时历曲线如图5所示。由图可知，艇体阻力及螺旋桨推力随时间逐渐向自航值逼近，且二者随时间的绝对变化率较为相近，航速随时间逐步增大并在t2时刻后趋于稳定；随着转速的增大，螺旋桨获得的推力以及艇体的阻力值越大，其运动稳定的时间就越短，最终稳定的航速值也越大，这从图6即可看出。表 4给出了不同转速下潜艇做定深直航运动时的计算结果。由表4 可知，4 个转速下 t1 值间的差别较小，在t2 时刻是随转速的增大而减小，转速越大，需要稳定的时间就越短；转速为300~ 600 r/min时，最终稳定的定深下的航速分别为0.97，1.35，1.72 和 2.1 m/s。
3.2 艇体上浮运动
图 7给出了不同转速下艇体垂向位移运动的时历曲线。由图可知，艇体做应急上浮至出水一般会经历水下定深直航运动—上浮运动—水面航行这3种状态的运动。当螺旋桨转速为300 r/min时，其水下定深直航时间为 101 s ，上浮时间 8 s，上浮出水的时刻为 109 s ；当转速为 400 r/min 时，其水下定深直航时间为 92 s，上浮时间 9 s，上浮出水的时刻为 101 s；当转速为 500 r/min 时，其水下定深直航时间为 91 s，上浮时间 9 s，上浮出水的时刻为 100 s；当转速为 600 r/min 时，其水下定深直航时间为71 s，上浮时间 12 s，上浮出水的时
刻为 83 s。由此可知，随着转速的增大，潜艇获得稳定航速的时间越短，艇体在水下做定深直航运动的时间也就越短；不同转速下艇体上浮的时间虽然有所差别，但其差值量级较小，上浮运动的时间主要还是受上浮力的影响；同时，潜艇的上浮运动是一个非线性的复杂的六自由度运动，当上浮力一致时，螺旋桨转速越大，其上浮运动速度不一定越大。
3.3 艇上浮时的姿态分析3.3.1 纵倾角分析
图 8所示为不同转速下艇体纵倾角的时历曲线，表 5 给出了纵倾角第1 次产生变化的时刻、发生峰值的时刻、峰值以及来回振荡波动的时间。由图8 和表5 可知，当转速为 300 r/min 时，纵倾角第1次产生变化的时刻是 100.2 s，第 1次产生峰值的时刻是 109 s ，其峰值为−22°，纵倾角来回振荡波动的时间是29.8 s；当转速为 400 r/min时，纵倾角第1次产生变化的时刻是92 s，第 1次产生峰值的时刻是99 s，其峰值为−23°，纵倾角来回振荡波动的时间是47 s；转速为 500 r/min 时，纵倾角第1次产生变化的时刻是 91 s，第 1 次产生峰值的时刻是96 s，其峰值为−22°，纵倾角来回振荡波动的时间是47 s；当转速为 600 r/min 时，纵倾角第1次产生变化的时刻是 71 s，第 1 次产生峰值的时刻是75 s，其峰值为−27°，纵倾角来回振荡波动的时间是54 s。结果显示，随着转速的增大，纵倾角第1次产生变化及峰值的时刻呈减小趋势，而纵倾角来回振荡波动的时间则越来越大，即转速越大，艇体出水时受自由液面扰动越小。纵倾角的峰值变化主要发生在潜艇上浮运动过程
中，来回振荡波动主要发生在潜艇出水时刻。对比 3.2节有关上浮的时刻可知，无论转速多大，潜艇第1次倾角的峰值都是发生在潜艇上浮运动过程中，由对比分析可知，随着转速的增大，纵倾角的峰值并未呈现出明显的变化规律。
3.3.2 横倾角分析
图 9所示为不同转速下艇体横倾角的时历曲线，表 6 给出了横倾角第1 次产生变化的时刻、发生峰值时刻、峰值大小以及来回振荡波动的时间。由图9 和表6 可知，当转速为 300 r/min 时，横倾角第1 次产生变化的时刻是 101 s，第 1次产生峰值的时刻是 108 s，其峰值为 28.5°，横倾角来回振荡波动的时间是 87 s ；当转速为 400 r/min时，横倾角第1次产生变化的时刻是92 s，第 1次
产生峰值的时刻是 99 s ，其峰值为 22.3°，横倾角来回振荡波动的时间是 79 s ；当转速为 500 r/min时，横倾角第1次产生变化的时刻是88 s，第 1次产生峰值的时刻是 98 s ，其峰值为 18.1°，横倾角来回振荡波动的时间是 73 s ；当转速为 600 r/min时，横倾角第1次产生变化的时刻是70 s，第 1次产生峰值的时刻是 78 s ，其峰值为 22.2°，横倾角来回振荡波动的时间是69 s。结果显示，随着转速的增大，横倾角第1次产生变化时刻和发生峰值的时刻以及来回振荡波动的时间越来越小。横倾角的峰值变化主要发生在潜艇上浮运动过程中，来回振荡波动主要发生在艇出水时刻。对比3.2节有关上浮的时刻可知，无论转速是多少，潜艇第1次横倾角的峰值均发生在潜艇上浮运动过程中，随着转速的增大，横倾角的峰值并未呈现出明显的变化规律。
4 结 语
本文基于体积力模型的虚拟桨盘模型，对不同转速下潜艇的水下定深直航、六自由度上浮、出水以及水面航行进行了数值模拟，并对比分析了不同转速下潜艇螺旋桨与转速的匹配及其姿态等，该方法可为潜艇的应急上浮运动提供一定的参考。通过对比分析不同转速下艇体上浮运动过程中的姿态，可得出如下结论：随着转速的增大，潜
艇获得稳定航速的时间越短，艇体在水下做定深直航运动的时间也越短；不同转速下其上浮运动的时间相差不大；转速越大，无论是纵倾角还是横倾角，其倾角第1次产生变化的时刻和发生峰值的时刻会越来越小；同时随着转速的增大，纵倾角的来回震荡波动时间会越来越长，而横倾角的来回震荡波动时间则越来越短，且倾角的峰值随转速的变化并未呈现出明显的变化规律。
参考文献：
[1] SMALLWOOD D A, WHITCOMB L L. Adaptive identification of dynamically positioned underwater robotic vehicles[J]. IEEE Transactions on Control Systems Technology, 2003, 11(4): 505–515.
[2] WATT G D. Estimating the hydrodynamic characteristics of the dorado remote minehunting vehicle[R]. [S. l.]: Defence Research and Development Canada, 2002.
[3] KIM J, KIM K, CHOI H S, et al. Estimation of hydrodynamic coefficients for an AUV using nonlinear observers[J]. IEEE Journal of Oceanic Engineering, 2002, 27(4): 830–840.
[4] WATSON K P, WEBSTER J S, CRANE J W, et al. Prediction of submersible maneuvering performance at high incidence angles[C]//Proceedings of OCEANS'93. Victoria, BC, Canada: IEEE, 1993: 289–294.
[5] 贾欣乐,杨盐生.船舶运动数学模型——机理建模与辨识建模[M]. 大连:大连海事大学出版社, 1999. JIA X L, YANG Y S. Mathematical model of ship motion[M]. Dalian: Dalian Maritime University Press, 1999 (in Chinese)..
[6] 戴余良,俞科云.潜艇应急上浮运动稳定性与分叉分析[J].上海交通大学学报, 2010, 44(10): 1400–1404, 1413. DAI Y L, YU K Y. Stability and bifurcations analysis of submarine during emergency ascent[J]. Journal of Shanghai Jiao Tong University, 2010, 44(10): 1400–1404, 1413 (in Chinese).
[7] 陈志法, 戴余良, 林雄伟,等.潜艇水下高速运动稳定性分析[J]. 舰船科学技术, 2015, 37(3): 15–20. CHEN Z F, DAI Y L, LIN X W, et al. Motion stability analysis of high-speed submarines under water[J]. Ship Science and Technology, 2015, 37(3): 15–20 (in Chinese).
[8] 王晓玢, 孙尧,莫宏伟.潜艇操纵运动分叉突变特性[J].工程力学, 2009, 26(10): 252–256. WANG X F, SUN Y, MO H W. Bifurcation and catastrophe characteristics in submarine motion[J]. Engineering Mechanics, 2009, 26(10): 252–256 (in Chinese).
[9] CARRICA P M, ISMAIL F, HYMAN M, et al. Turn and zigzag maneuvers of a surface combatant using a URANS approach with dynamic overset grids[J]. Journal of Marine Science and Technology, 2013, 18(2): 166–181.
[10] CARRICA P M, SADAT-HOSSEINI H, STERN F. CFD analysis of broaching for a model surface combatant with explicit simulation of moving rudders and rotating propellers[J]. Computers & Fluids, 2012, 53: 117–132.
[11] CHASE N, MICHAEL T, CARRICA P M. Overset simulation of a submarine and propeller in towed, self-propelled and maneuvering conditions[J]. International Shipbuilding Progress, 2013, 60(1/2/3/4): 171–205.
[12] 周广礼, 欧勇鹏, 高霄鹏.潜艇上浮出水运动及粘性流场直接数值模拟方法[J].海军工程大学学报, 2017, 29(4): 86–92. ZHOU G L, OU Y P, GAO X P. Numerical simulation method of surfacing maneuvers and flow field of submarine[J]. Journal of Naval University of Engineering, 2017, 29(4): 86–92 (in Chinese).
[13] 周广礼, 董文才, 欧勇鹏.潜艇应急上浮六自由度运动及黏性流场数值模拟[J].国防科技大学学报, 2017, 39(2): 199–206. ZHOU G L, DONG W C, OU Y P. Numerical simulation of six degree of freedom motion and viscous flow for submarine's emergency ascent[J]. Journal of National University of Defense Technology, 2017, 39(2): 199–206 (in Chinese).
[14] 廖欢欢, 庞永杰,李宏伟,等.基于重叠网格技术的潜艇应急上浮空间运动的数值模拟[C]//第二十七届全国水动力学研讨会文集.南京, 2015. LIAO H H, PANG Y J, LI H W, et al. Numerical simulation based on the overset grid technique about submarine emergency ascent space movement[C]//Proceedings of the 27th National Hydrodynamics Symposium. Nanjing: [s.n.], 2015 (in Chinese).
[15] WILCOX D C. Formulation of the k-ω turbulence model revised[J]. AIAA Journal, 2008, 46(11): 2823–2838.
[16] MENTER F R, KUNTZ M, LANGTRY R. Ten years of industrial experience with the SST turbulence model[C]// Turbulence, Heat and Mass Transfer. Antalya: Begell House, 2003.
[17] FELDMAN J. Revised standard submarine equations of motion: DTNSRDC report: SPD-0393-09 [R]. Bethesda,Maryland: David Taylor Naval Ship Research and Development Center, 1979.