潜艇在大国的海军中有着极其重要战略地位，潜艇下潜后再上浮时的稳定性对潜艇姿态和艇中人员安全是极其重要的．周广礼等[1-3]以SUBOFF模型为研究对象，基于RANS方程及体积力模型，并采用整体动网格技术，对潜艇的上浮及其出水运动进行了数值计算；黄祥兵等[4-5]设计了潜艇大型集体逃生舱安全上浮优化方案，对其上浮状态进行数值模拟，计算了上浮过程中的最大速度；吕帮俊等[6]采用了一种新形式的大攻角非稳态潜艇操纵运动仿真数学模型，对高压吹除主压载水舱后潜艇应急上浮过程的运动状态进行了预报．以上研究都是基于潜艇在静水中的上浮运动，然而潜艇在实际的海面上上浮，遭遇的波浪是极其不规则的，每一个波的波高、波长和周期都是随机变化的，即潜艇一般在海面上浮遭遇的波浪是不规则波．关于潜艇在不规则波中上浮运动的研究很少，因此有必要探索研究潜艇在不规则波中的上浮运动规律．随着计算技术的快速发展，建立数值波浪水池应用于数值模拟的条件越来越成熟，与试验相比，数值波浪水池[7-9]具有成本低、效率高、且能获得较好的流场信息等优点．这里基于波浪理论，构造了数值波浪水池，对不规则波的波浪模拟进行了分析，以SUBOFF为研究对象，对其在不规则波中的上浮运动开展研究．1 计算方法1.1　控制方程潜艇的上浮运动是从单介质向多介质变化的复杂流体运动问题，应用CFD数值计算的控制方程主要包含连续性方程和动量方程，方程的形式为∂∂xi(ρui)=0;    ∂∂t(ρui)+∂∂xj(ρuiuj)=-∂p∂xi+∂∂xj(μ∂ui∂xj-ρui'uj'¯)+ρfi, (1)式中：ρ为流体密度；xi，xj为i，j方向的坐标；ui，uj为流体速度沿i，j方向的分量；ui'uj'¯为无因次速度的时均量；fi为单位质量力．1.2　湍流模型潜艇上浮运动是一个须要考虑湍流边界层的运动，针对潜艇上浮运动的数值计算主要选取k-ω Turbulence[10]，综合相关文献和资料[11]，选取SST k-ω湍流模对潜艇上浮运动进行计算，具体为    ∂(ρk)∂t+∂(pkuj)∂xj=∂∂xjμ+μtσk∂k∂xj+Pk-β*ρωk;    ∂(ρω)∂t+∂(pωuj)∂xj=∂∂xjμ+σωμt∂ω∂xj-γυtp-βρω2+2(1-F1)ρσω2ω∂k∂xj∂ω∂xj,(2)式中：k为湍动能；p为压力项；μ为流体动力黏性系数；μt为湍动黏度；ω为特殊湍动能耗散率；σk和β*为模型系数，分别取1-0.15F1和0.09；Pk=min(P,20β*ρkω)为结果项；σω为ω对应的普朗特数；γ为模型转换系数；υt=k/ω为涡黏系数；β为热膨胀系数；F1为混合函数．1.3　自由液面方程自由液面的波动采用VOF (volume of fluid)方法[12]进行模拟，具体为∂Cq∂t+u∂Cq∂x+v∂Cq∂y+w∂Cq∂z=0;∑q=12aq=1, (3)式中：Cq为自由液面网格单元内流体体积与总体积的比值；q=1表示空气项，q=2表示水项；aq为体积分数，这里默认aq=0.5处为自由液面．2 计算对象及网格划分2.1　计算对象SUBOFF为美国泰勒研究中心提供的潜艇试验模型，国内外各水池针对该模型开展了大量的拖曳及流场测量试验，这里在标准的SUBOFF模型基础上，在其指挥台围壳上加有舵，模型的艇长为4.356 m，艇宽为0.508 m，水下排水体积为0.708 m3，重心位置为(2.017，0，-0.01)m，三维模型图如图1所示．10.13245/j.hust.220320.F001图1艇体三维模型(m)2.2　坐标系为了研究潜艇的上浮运动规律，确定潜艇运动中的位置和姿态，采用E-ξηζ为大地坐标系．潜艇的上浮运动是一个六自由度多维运动，主要可分为水平面运动和垂直面运动，因此其随艇坐标系可分解为水平面运动坐标系和垂直平面运动坐标系，如图2所示．10.13245/j.hust.220320.F002图2潜艇的运动坐标系潜艇重心G处相对于大地的速度为V，航速V的方向与Gx轴之间的夹角β为漂角，α为来流的攻角；V在随艇坐标系上的投影为(u,v,w)，潜艇所受的外力F在随艇坐标系上的投影为(Fx,Fy,Fz)，潜艇的位置坐标和姿态角，可用随艇坐标系原点的地面坐标值(ξ,η,ζ)和相对于大地坐标系的3个姿态角(φ,θ,ψ)来确定，其中：φ为横倾角；θ为纵倾角；ψ为航向角．2.3　网格划分对于潜艇上浮的数值模拟，采用重叠网格方法进行，将计算域分为背景域及重叠网格域，重叠网格域内嵌于背景域中，重叠网格域随艇体运动，如图3所示．10.13245/j.hust.220320.F003图3潜艇上浮运动重叠网格应用示意图计算域采用方形域，尺寸大小为9L×4L×4.8L，为准确地表达波形，在计算域中波面附近设定了多层加密块，网格总数量为1.056×107．3 数值不确定度分析采用STERN等[13]提出的不确定分析方法，评估潜艇上浮运动数值计算结果的不确定度，验证该数值计算方法的可靠性，不确定度分析分为验证和确认两个过程[14]．3.1　验证验证过程主要针对网格收敛性和时间步长收敛性．进行网格收敛性分析时，时间步长Δt=0.005 s，根据网格加细比rG=2可知网格基本尺寸可分为1，1.414 2，2，简称为细、中、粗3种网格，网格的具体信息如表1所示．进行时间步长收敛性分析时，选取细网格，根据时间步长的加细比rT=2，选取3个时间步长Δt1=0.005 s，Δt2=0.007 s，Δt2=0.010 s．10.13245/j.hust.220320.T001表 1不同网格类型网格类型无量纲基本尺寸网格数量/104背景域网格重叠域网格网格总数细1215.3743.2958.5中1.414 2113.4435.5548.9粗276.5261.3337.8表2和表3分别给出了潜艇在静水中3个正上浮力与艇重的占比分别为6%，9%和12%时，艇体零速自由上浮中垂向速度最大峰值的网格和时间步长的收敛性结果．表2中：ΔF为正上浮力占比，SG1，SG2，SG3分别为网格收敛性研究中细、中、粗网格的垂向速度最大峰值的数值计算结果；UG为网格不确定度；UGC为修正值的不确定度；D为对应的试验值，基于模型上浮试验可求得．表3中：ST1，ST2，ST2分别为时间步长收敛性研究中Δt1=0.005 s，Δt2=0.007 s，Δt2=0.010 s的垂向速度最大峰值的数值计算结果；UT，UTC的含义与表2类似，下标T表示对时间步长收敛性的研究．10.13245/j.hust.220320.T002表 2网格收敛性验证结果ΔF/%SG1SG2SG3UGUGCD60.7820.7660.7430.0780.0630.85690.9500.9310.8980.0820.0631.060121.2611.2361.1940.1140.0891.32410.13245/j.hust.220320.T003表 3时间步长收敛性验证结果ΔF/%ST1ST2ST3UTUTCD60.7820.7710.7480.0430.0320.85690.9500.9340.9010.0640.0481.040121.2611.2381.1870.0860.0631.3243.2　确认表4分别给出了潜艇静水中3个不同正浮力下零速自由上浮时艇体垂向速度最大峰值的数值计算的确认结果．数值确认不确定度UV=USN2+UD2，由网格收敛性和时间步长收敛性的数值计算结果可求得该数值计算不确定度USN=UG2+UT2，UD为试验不确定度，在此取2%D，由表4可知：数值计算结果与试验值的对比误差E均小于确认不确定度UV，由此可知该数值计算结果均得到确认，从而证明了上浮运动数值计算的可靠性．10.13245/j.hust.220320.T004表 4数值计算的确认结果ΔF/%USN/DUV/DE/D610.4310.628.64910.0110.218.651210.8110.994.764 数值波浪水池4.1　数值造波依据线性叠加原理，假定不规则波是有许多不同波长、不同波幅和随机相位的单元波叠加而成的，则不规则波的波面函数和速度场可表示为ζ(x,t)=∑i=1Naicos(kix-ωit+φi)；(5)u(x,z,t)=∑i=1Naiωiekizcos(kix-ωit+φi),v=0,w(x,z,t)=∑i=1Naiωiekizsin(kix-ωit+φi), (6)式中：ai，ωi，ki，φi分别为第i个组成波的波幅、圆频率、波数和初相位，φi在(0,2π)范围内均布随机取值；u(x,z,t)，v和w(x,z,t)分别为波浪在纵向、横向和垂向上的速度分量．不规则波的模拟是通过频谱的模拟来实现的，设模拟的对象谱Sζ(ω)的能量大多数分布在ωL∼ωH频率范围内，把波频划分为N个区间，则每个区间Δω=(ωH-ωL)/N，并定义波频ωi=ωL+(i-0.5)Δω，进而可求解波幅ai=2Sζ(ωi)Δω．在求得每个子波的波幅、相位以及波长后，调用STAR-CCM+软件中的Java宏命令，完成对子波参数的赋值，通过线性叠加的方式，将子波的参数添加到叠加的波形中．4.2　波浪谱为了研究潜艇在不规则波中上浮出水问题，选取了ITTC双参数谱作为海浪谱密度的估算，ITTC双参数谱的具体形式为Sζ(ω)=173H1/32T14ω5exp-691T14ω4，(7)式中：H1/3为有义波高；T1为谱形心周期；ω=2π/T1为波频．4.3　不规则波的波浪模拟分析潜艇在不规则波中开展上浮运动之前，首先对造取的不规则波的可行性进行分析，选取有义波高为0.080，0.128和0.180 m，对应的谱形心周期为1.25，1.50和1.60 s等3个不同的有义波高和不同谱形心周期的不规则波开展数值计算．10.13245/j.hust.220320.F004图4H1/3=0.08 m，T1=1.25 s不规则波的波面时历曲线10.13245/j.hust.220320.F005图5H1/3=0.128 m，T1=1.5 s不规则波的波面时历曲线图4~6分别给出了有义波高H1/3=0.08 m、谱形心周期T1=1.25 s，H1/3=0.128 m、谱形心周期T1=1.5 s，H1/3=0.18 m、谱形心周期T1=1.6 s在不同纵向位置处(艇艏为坐标原点，分别取原点处和艇艏后方3L处)不规则波的波高时历曲线．可以看出：无论是哪种有义波高和谱形心周期，在x=0(上游)和x=3L(下游)处波高的变化规律随时间的推进而并未呈现一定的衰减，由此证明构造的不规则波可用于研究潜艇上浮运动．表5给出不规则波的有义波高和谱形心周期的计算值与目标值的对比，可见有义波高和谱形心周期的误差绝对值都在4%～6%之内．图7给出了不同有义波高和谱形心周期下波浪谱的对比，从图中可知：数值计算的波浪谱与目标波浪谱的变化规律高度一致，且较为相符．综上所述：该不规则波的模拟有较高的精度，具有很好的可行性，可为潜艇在不规则波开展上浮运动的研究奠定基础．10.13245/j.hust.220320.F006图6H1/3=0.180 m，T1=1.6 s不规则波的波面时历曲线10.13245/j.hust.220320.T005表 5不规则波的有义波高和谱形心周期序号H1/3T1计算值/m目标值/m误差/%计算值/m目标值/m误差/%10.075 60.080 0-5.501.321.255.6020.121 00.128 0-5.471.591.506.0030.172 50.180 0-4.171.671.604.3810.13245/j.hust.220320.F007图7不规则波波浪谱模拟的对比5 不规则波中的上浮运动在开展不规则波的上浮运动计算时，先初始化波形，此时潜艇是处于水下的静止状态，不释放任何自由度，待波形完全形成后再释放潜艇6个方向自由度使其上浮．选定3个不同的有义波高H1/3=0.080 0，0.128 0和0.180 0 m，对应的谱形心周期T1=1.25，1.50和1.60 s，浪向为迎浪，正上浮力占比为9%，上浮的潜深为3.5 m．为了更容易表征潜艇上浮运动的规律，使其更具普适性，将潜艇上浮及出水时艇体的受力F、速度u、位移ξ等参数进行无量纲化，L为艇长，B为艇宽，其无量纲公式为：F'=2F/(ρLBV2);u'=u/V;ξ'=ξ/L;V=u2+v2+w2.图8和图9分别给出了在不同有义波高的不规则波中上浮到出水时潜艇的纵向、横向和垂向等3个方向的力、速度、位移和姿态角的时历变化结果，图中的虚线表示潜艇上浮和出水的分界线，即虚线以左为潜艇的上浮过程，虚线以右为潜艇出水后的状态．从图中可知：在上浮出水以前，由于潜艇不受波面的影响，不同有义波高下艇体受到的力差别较小，在潜艇出水以后，受不规则波波浪的影响，发生一定的摇摆运动，当不规则波中有义波高的波高值越大，潜艇发生摇摆的幅度越大，因此潜艇受到的力和力矩越大，即潜艇在纵向方向受到的阻力也就越大，导致潜艇在出水后再纵向方向的速度值和位移就越小；潜艇在上浮未出水时，艇体受波浪的影响很小，因此艇体的横倾角和纵倾角的峰值基本一致，待潜艇在上浮出水后，有义波高的波高值越大，潜艇的垂向运动响应越大，因此潜艇垂向位移的波动越大，艇体的横倾角和纵倾角的波动值越大．10.13245/j.hust.220320.F008图8不同有义波高下艇体上浮受力和速度的时历曲线10.13245/j.hust.220320.F009图9不同有义波高下艇体上浮位移和姿态角的时历曲线采用CFD软件中Java宏命令功能，基于ITTC双参数谱，实现了不同有义波高下不规则波的模拟，同时通过添加人工阻尼项，实现了消波功能，构造了不规则波的数值波浪水池．基于构造的波浪水池计算了不规则波有义波高和谱形心周期，并与目标值进行了对比，发现误差的绝对值都在4%～6%之内；同时，对数值计算的波浪谱与目标谱也进行了对比，发现其变化规律高度一致，且较为相符．综上可见：该不规则波的模拟有较高的精度，具有很好的可行性，该数值计算方法可为潜艇在波浪中的上浮运动研究提供一定的借鉴．基于构造的数值波浪水池，对迎浪中不同有义波高下潜艇的上浮运动开展了研究，发现潜艇在不规则波浪中上浮出水后，有义波高的波高值越大，潜艇的垂向运动响应越大，艇体的姿态角越大，可为实艇在波浪中的上浮提供一定的理论指导．