附体对分层流中潜艇水动力特性的影响

2021-03-16刘双何广华王威高云

兵工学报 2021年1期

刘双，何广华,，王威，高云

(1.哈尔滨工业大学机电工程学院, 黑龙江哈尔滨 150001;2.哈尔滨工业大学(威海) 船舶与海洋工程学院, 山东威海 264209)

0 引言

潜艇是各国发展军事力量不可或缺的一环，研究潜艇的水动力特性对建立潜艇作战系统至关重要。实际海水中存在的温、盐差会使水体出现分层现象，潜艇航行于多层流之间时会产生内波[1]，针对潜艇航行兴起内波尾迹的研究具有重要意义。

各国学者采用理论及数值方法对航行体源致内波进行了研究[2-4]。在兴波尾迹研究方面，Yeung等[5]通过求解Green函数，探讨了运动潜体产生的表面波模式及内波模式对自由液面和界面处兴波的影响；赵先奇等[6]对细长体在三层流体中运动生成内波的问题进行了模型实验研究；Chang等[7]研究了双层流中航行潜艇的尾迹特性；Chomaz等[8]研究了分层流体中球体的近场尾流情况。在分层流中物体受力研究方面，Motygin等[9]基于势流理论探讨了二维物体分别在上层、下层流体中运动所产生的波浪阻力公式；Grue[10]对“Fram”号极地船航行过程进行研究，得出的受力变化规律与Ekman[11]一致；勾莹等[12]采用模型实验方法研究了箱型结构在双层流中拖动时的阻力特性，并与单层流中拖航阻力实验结果进行了对比研究；Wei等[13]、Wang等[14]针对内孤立波与水下潜体相互作用问题进行了实验研究。在潜艇附体研究方面，Bull[15]通过改变湍流模型、求解器及网格形式，系统地验证了SUBOFF潜艇尾流场数值模拟方法；Bensow等[16]采用大涡模拟(LES)方法计算了SUBOFF潜艇的黏性流场；Posa等[17]研究了雷诺数对全附体SUBOFF潜艇尾迹特征的影响；Yao等[18]采用LES及边界元方法对全附体潜艇周围的流场及激流噪声进行了研究；Seil等[19]研究了附体及其加装方式对BB2型潜艇水动力载荷的影响。

综上可知，尾迹研究常常仅关注界面波模式，对于航行体受力，尤其是航行体外形对其自身水动力特性的影响研究较少。针对密度分层流这一特殊海况，分析各个附体对潜艇水动力特征带来的影响显得十分必要。

本文采用黏性计算流体力学(CFD)理论，基于工程物理仿真软件STAR-CCM+平台建立了密度分层流数值仿真模型。首先进行收敛性验证，并将数值结果与实验结果进行对比，验证本方法的准确性；再以SUBOFF潜艇为基础研究其在裸艇、仅有驾驶舱附体、仅有尾翼附体及全附体工况下以不同航速在密度分层流中航行时的兴波及阻力特征，分析各附体工况下航速对潜艇所受阻力的影响以及阻力系数峰值的出现与艇艏、艉兴波干扰的关系；并结合压力分布探讨附体存在带来的影响。本模型具有较好的计算精度，可为潜艇在密度分层流中航行时的水动力性能预报提供新的手段，同时对潜艇设计优化与规避策略提供参考。

1 数值计算模型

1.1 控制方程

本模型基于雷诺平均(RANS)方程，不可压缩牛顿流体的控制方程为

(1)

(2)

(1)式为连续性方程，(2)式为动量方程，二者均为时均计算处理后的形式。

1.2 湍流模型

为使(2)式封闭，需要引入湍流模型。采用Rea-lizablek-ε湍流模型(k为湍动能，ε为耗散率)，该模型稳定性良好，压力梯度求解精度高，工程应用比较广泛，其湍动能及耗散率输运方程为

(3)

(4)

式中：xt表示在x轴方向上对时间t求偏导；μ为分子扩散所造成的动力黏性；μt为湍流黏性系数；σk和σε分别为湍动能k和耗散率ε的湍流普朗特数；Gk为平均速度梯度产生的湍动能；Gb为由浮力产生的湍动能；YM表示湍流脉动膨胀对总耗散率的影响；C1=max[0.43,η/(η+5)]，系数η=Sk/ε，S为平均应变率；C2和C1ε为常数；C3ε为浮力对耗散率影响的函数，C3ε=tanh|up/uv|,up为平行于重力方向的速度分量，uv为垂直于重力方向的速度分量。参数取值为C2=1.9、C1ε=1.44、σk=1.0、σε=1.2，在文献[20]中已证实是可靠的。

1.3 平面进行波理论

根据平面进行波理论分析兴波阻力特性，对水下航行潜艇的快速性研究具有参考价值。宽度为b、长度为2λ的封闭区域中平面进行波波能Eb为

(5)

式中：g为重力加速度；H为波高。

平面进行波理论中兴波阻力系数Cw可表示为

(6)

式中：Rw为兴波阻力；v为航速；mL表示兴波长度，m为系数，与船型和弗劳德数有关，L为艇长；A为湿表面面积；C、D为常数；Fr为弗劳德数。

由于cos(2πmL/λ)值是在1.0和-1.0之间变动，因此兴波阻力系数Cw-Fr曲线存在波动。

1.4 流体体积方法

通过流体体积方法进行界面捕捉，引入变量相体积分数αq，定义每个控制单元中第q相的体积分数为

(7)

式中：Vq为第q相在控制单元中所占的体积；V为控制单元的体积。当αq=0时，表示在控制单元内不存在第q相流体；当αq=1时，表示在控制单元中充满第q相流体；当0<αq<1时，说明在控制单元中存在第q相流体和其他相流体的界面。

2 数值模型与验证

2.1 计算模型可靠性验证

分层流中潜艇阻力及兴波特性的CFD研究较少，研究基于STAR-CCM+软件进行仿真，首先对数值模型进行验证。根据文献[12]中的实验装置参数进行建模如图1所示。其中水深为0.6 m，淡水在上层，密度ρw=997.56 kg/m3；下层是密度ρs=1 024.00 kg/m3的盐水；上层淡水厚hw为0.30 m，下层盐水厚hs为0.30 m. 实验模型为60 cm×45 cm×35 cm的方型浮箱，吃水深度定为0.2 m，速度分别为0.14 m/s、0.16 m/s、0.18 m/s、0.20 m/s、0.24 m/s.

图1 实验装置示意图[12]Fig.1 Schematic diagram of test device[12]

与相关实验结果的对比见图2，其中数值模拟结果为箱体在双层流中拖航时所受的总阻力，用Rtot表示。由图2可以看出，数值模拟结果与文献[12]中的实验数据基本吻合，说明本文所采用的CFD计算模型可靠。

图2 本文数值模拟结果与实验数据[12]对比Fig.2 Comparison of the simulated results in the present paper with the experimental data in Ref.[12]

2.2 数值模拟相关设置

采用对称面条件来减少计算资源的消耗，提高计算效率，故模拟得到的是半艇阻力值，需乘以2来展示全艇阻力。为了防止数值模拟过程中回流的产生，适当延长计算域长度，取艇前3倍艇长，艇后12倍艇长。潜艇表面边界条件类型为壁面，物面条件为无滑移。基于欧拉多相流模型，通过用户自定义场函数(UDF)的方式分别给定空气、淡水、盐水3相的密度及初始体积分数，使流体域分层。入口边界条件类型为速度入口，在此处指定流体速度；出口边界条件类型为压力出口，在此处通过UDF方式给定分层流体的压力分布。

对各流体相分界面及潜艇附体处网格进行加密以准确捕捉潜艇在分层流中航行时的流场特性。此外，使用切割体网格并结合棱柱层网格生成器来处理边界处的网格，其中棱柱层厚度为0.018 m，棱柱层数为6，增长率为1.2. 调节网格增长率等参数控制网格质量，如图3所示，通过计算机(中央处理器i7-8700,主频3.19 GHz)进行计算，一个算例耗时30 h左右。

图3 潜艇周围网格分布Fig.3 Grid distribution around submarine

研究采用的潜艇模型为标准SUBOFF潜艇模型，分别选择4种附体工况(裸艇、仅有驾驶舱、仅有尾翼、全附体)，其中全附体SUBOFF模型见图4，具体尺寸如表1所示，潜艇潜深为h，即潜艇重心到静水面的垂向距离。

图4 SUBOFF全附体潜艇模型Fig.4 Full-appendage SUBOFF model

表1 潜艇尺寸参数

计算域中流体密度沿垂向变化，分层情况如图5所示。SUBOFF潜艇以不同速度航行于上层即淡水层，潜深h=0.12L,潜艇重心距淡水- 盐水交界面的距离同样为0.12L.

图5 分层示意图Fig.5 Schematic diagram of density stratified fluid

2.3 网格及时间步长收敛性验证

采用加装全附体的潜艇以航速Fr=0.5进行收敛性研究，由于网格与航速、波长、波陡等诸多因素有关，参考文献[21]的建议，并经过一系列数值模拟最终找到适合本文的经验方式。数值模型网格划分的宗旨是：保证数值模型能够在潜艇附体周围、交界面等处有足够的网格密度，使单位兴波波长内大于80个网格、单位波高内大于15个网格，用以捕捉潜艇在各交界面处的兴波特征。通过改变网格的基础尺寸来控制网格数量，取时间步长为0.01 s，结果如表2所示。其中误差是指采用某一网格量计算结果与采用上一个较少网格量所得结果之间的差别；表2中Rf代表摩擦阻力,在STAR-CCM+软件中，摩擦阻力通过(8)式计算：

(8)

式中：Tn为面n上的应力张量；an为面网格面积矢量。

表2 网格收敛性分析

由表2可以看出，当网格基础尺寸由网格3降低到网格4时，阻力计算结果几乎没有变化，证明可采用网格3方式进行高效的CFD研究。

为找到快速、稳定的时间步长方案，参考文献[21]中的建议，时间步长Δt按照0.005L/U～0.01L/U，并考虑库朗数的要求来给出时间步长的预估值，再通过大量的数值实验最终确定时间步长。表3为采用网格3网格方案进行的时间步长收敛性研究，可见随着时间步长降低，计算逐渐达到收敛。当时间步长达到0.01 s时，再降低时间步长到0.005 s对计算结果的影响不大，故本研究中均选择时间步长为0.01 s.

表3 时间步长收敛性分析

3 数值结果与分析

3.1 不同附体类型对潜艇快速性的影响

3.1.1 阻力变化

图6 阻力随航速的变化Fig.6 Change of resistance with navigational speed

模拟中计算域如图5所示，淡水层总厚度为0.24L，密度取为997.56 kg/m3，盐水层厚为0.67L，密度为1 020.00 kg/m3. 根据文献[20]，潜艇离交界面越近、兴波情况越明显，故使潜艇位于潜深h=0.12L处，同时潜艇中心距内波面距离同样为0.12L，以便观察潜艇附体改变时其在两个界面处兴波的变化。针对每一种附体工况均模拟了Fr为 0.3、0.4、0.5、0.6、0.7、0.8 6种速度，在研究附体影响同时分析潜艇水动力性能随航速的变化。图6为各附体工况下总阻力及摩擦阻力随航速的变化。

由图6可知，随着航速增加，潜艇的总阻力和摩擦阻力不断增加。低速时，附体对阻力的影响较小；高速时，附体对阻力有增大的作用；尾翼对摩擦阻力的影响要略大于驾驶舱，但二者对总阻力的影响几乎一致。

3.1.2 阻力系数变化

密度分层流中各附体工况下潜艇兴波阻力系数随航速的变化趋势如图7所示。由图7可以看出，曲线随航速的增加出现先上升、后下降的趋势，这是由于在兴波阻力系数计算(6)式中存在一个余弦项，此项会根据潜艇艏、艉兴波的相互干扰情况在-1～1之间变化。低速状况下，仅加装尾翼时兴波阻力系数最低，加装全附体工况时系数最高；高速时各附体对阻力系数的影响不大；整体来看，尾翼的存在会使兴波阻力系数降低。

图7 阻力系数随航速的变化Fig.7 Change of resistance coefficient with navigational speed

各附体工况下：在Fr=0.5附近兴波阻力系数均最高，说明此时艇艏兴波与艇艉兴波发生了不利干扰；Fr=0.5时仅加装尾翼工况的阻力系数最低，而仅加装驾驶舱的潜艇阻力系数最高。

可见，附体的不同会对阻力和阻力系数带来一定影响，且这种影响与航速有关。

3.2 不同附体存在对潜艇兴波特性的影响

潜艇在密度分层流中航行时，会在自由液面及内波面处产生兴波，兴波尾迹特征可以作为水下潜体探测的重要参考，兴波会随着潜艇航速的变化而改变。本节取Fr为0.3、0.5、0.7工况研究不同附体工况下潜艇以不同航速航行时在两个分层界面上的兴波及速度特性。

潜艇在裸艇与全附体工况下航行时在自由液面及内波面处的兴波尾迹随航速变化情况如表4和表5所示。由表4和表5可见：随着航速的增加，裸艇与全附体自由液面均呈Kelvin波形，兴波范围越来越广，波长逐渐增加；此外，相比低速(Fr=0.3)以及高速(Fr=0.7)情况，中速(Fr=0.5)附近的波峰峰值最大。内波波形与表面波存在较大区别，兴波将不呈现Kelvin波形，而是波峰集中在整个兴波区域的中间部分，并呈现峰谷交替情况，波形以狭长的V形向后方延伸，在整体兴波的尾部出现环状的兴波区域。且自由液面处第1个兴波必为波峰，称为伯努利水丘，而内波面的第1个兴波必为波谷。

表4 裸艇与全附体工况下自由液面兴波云图随航速变化

表5 裸艇与全附体工况下内波面兴波云图随航速变化

随着航速的增加，内波面兴波波长也会逐渐增加，并在Fr=0.5时波峰峰值较高。潜艇兴波会消耗能量，波的能量与波长的平方呈正比，故随着速度增加，潜艇所受阻力会出现上升的趋势；在Fr=0.5时峰值较高，说明此时潜艇艏、艉兴波发生了不利干扰，这与兴波阻力系数曲线所展示的情况一致。由于淡水、盐水之间的密度差远小于淡水与空气之间，且仅为其千分之一量级，故很小的扰动即会激发大振幅的内波，且其密度梯度小、恢复力弱，所以内波周期长、波长大，整体波形对速度变化的响应迟缓。

表6和表7分别为裸艇及全附体工况时自由液面和内波面上沿潜艇前进方向的水质点速度分布情况，与表4、表5一致，上侧为裸艇工况、下侧为全附体工况，潜艇航行方向为x轴负向。由表6和表7可见：裸艇及全附体自由液面与内波面速度云图均呈现V形；自由液面处水质点速度要高于内波面，说明内波面上水质点运动对潜艇航速的改变不敏感。

表6 裸艇与全附体工况下自由液面水质点速度云图随航速变化

表7 裸艇与全附体工况下内波水质点速度云图随航速变化

将裸艇与全附体工况对比可知，潜艇航行于密度分层流中时，无论从兴波尾迹、水质点速度的分布形状还是从波形中峰、谷及波长的数值来看，附体存在带来的影响均较小。经过分析研究，驾驶舱及尾翼附体单独存在对潜艇兴波尾迹及水质点速度分布同样影响极小，故不在展示其兴波及速度分布情况。

3.3 不同附体存在对潜艇上、下方波面抬升的影响

表8为潜艇在加装不同附体工况下以航速Fr为0.3、0.5、0.7航行时，其上方自由液面及下方内波面在y=0 m剖面上的波切线。表8中两条虚线之间为潜艇在x轴方向的位置，潜艇沿x轴负向航行。

由表8可以发现，随着航速增加，兴波波长不断加大。相比内波面处，自由液面处的波形清晰，变化比较明显，这与3.2节中阐述的内波面波形变化缓慢原因一致；自由液面上第1个波形为波峰，内波面上第1个波形为波谷；对比不同附体工况下的结果发现，低速时驾驶舱的存在会使艇艏上方波形产生波动(表8中方框)，高速时几乎没有影响；整体来看，附体改变对各个速度下交界面处兴波的波面抬升影响不大。

表8 各附体工况不同航速下兴波波面抬升

图8 不同航速下各附体工况潜艇表面压力分布Fig.8 Surface pressure distribution of submarine with different appendages at various speeds

3.4 不同附体存在对潜艇表面压力分布的影响

在Fr分别为0.3、0.5、0.7时加装不同附体工况下的潜艇表面压力分布情况如图8所示。

由图8可见：潜艇表面压力随着速度的增加而逐渐加大，其中潜艇艏、艉压力变化较大，而主体部分的压力几乎不变；当存在附体时，潜艇表面曲率突变处，如驾驶舱及尾翼附近(图8中方框)的压力会产生较大波动，且加装全附体时的压力要明显高于其他附体工况。分析原因是由于加装附体后潜艇表面出现拐点，改变了局部流体的速度，进而引起压力变动。

下面将结合具体压力数值进行定量分析，Fr=0.5时加装全附体潜艇的表面压力及上方自由液面兴波如图9所示，由3部分构成，最下方是潜艇位置，中间是潜艇表面压力，最上方为潜艇航行时其中纵剖面与自由液面波切线上的波面抬升。

图9 潜艇表面压力分析图Fig.9 Analysis diagram of submarine surface pressure

在水中航行的潜体，其压力组成可分为静压、动压两部分，其中：静压是垂直于潜艇表面以上单位流体所具有的静压能，可以通过公式p=ρgh来计算；动压是由局部速度变化所引起，可由总压减去静压得到。表9为各个分析点处的压力值。

表9中E、J为潜艇艏、艉两点，其压差为3 264.5 Pa，这是形成潜艇压差阻力的主要原因，由图8可知，此数值会随着潜艇航速的增加而不断增加。为便于对比，后续的压力计算值将以此数值为基准，用百分比形式展示。F点静压与G点静压的差值为潜艇艏、艉压差的40.5%，但总压差值却是潜艇艏、艉压差的272.2%，说明动压的变化占较大比例；再对比H、I两点，静压差为16.5%，而总压差为187.9%，H、I两点在潜艇上处于同一水平位置，上方波面抬升带来的静压变化很小，说明主要是由于动压的变化导致了总压的差别。分析以上现象的原因是由于F、I点分别处于驾驶舱及尾翼的拐点处，附体带来的曲率突变导致此处局部流体速度改变，动压急剧变化，进而影响了各点的压力分布。