2017, Vol. 39
近年来，空腔结构的阻力分布及空腔附近流交换问题成为研究的热点问题，空腔在海洋航行器和航空飞行器上广泛应用，如潜艇上的流水孔等。潜艇的流水孔是指在潜艇的非耐压非水密结构上，用于潜艇上浮下潜时供液体自由进出的开口。当流体流经流水孔时，流场的非定常性会产生各种尺度的漩涡及湍流引发的脉动，作用于艇体外表面导致压力波动，破坏表面湍流附面层内的流动进而形成流噪声，其对潜艇航行隐身性能产生一定的影响[1]，很可能会暴露潜艇的位置，严重影响潜艇的作战能力和隐身性能。
空腔流动包含多种复杂的物理现象，如剪切层不稳定性、边界层分离、漩涡运动、激波边界层互相干扰等。研究空腔流动的方法主要有实验方法、数值分析等方法。赖焕新等[2]运用大涡模拟的数值分析方法对空腔流动及气动噪声进行研究；邬明[3]以潜艇流水孔为对象研究了流水孔引起的阻力增量集中在导流板和流水孔内壁上；陈灿等[4]运用大涡模拟研究了不可压缩空腔流的特性及噪声；孟生等[5]以潜艇流水孔为研究对象，通过改变流水孔的模型参数，得出了流水孔的噪声特性；张楠等[6]以潜艇流水孔为对象通过求解RANS方程，结合PISO算法精确仿真了流水孔内外流场的情况。Rossiter等[7]通过实验得出了空腔内存在涡和噪声的反馈循环，从而提出了空腔流动振荡频率的模型，并给出了预估振荡频率的半经验公式；涂海文等[8]运用RNG，k-ε湍流模型对潜艇的阻力分布和流场分布进行了模拟，为进一步优化潜艇的艇型和分析潜艇的流噪声打下了基础；李霖等[9]运用大涡模拟对于圆柱绕流进行了模拟，显示出这种数值模拟方法的优势。
本文采用大涡模拟方法，以Suboff艇型为对象，研究开孔纵向分布对空腔流动特性的影响，包括内外流的交换，频谱特性及其对阻力的影响。
1 几何模型取Suboff主艇体母线建立二维模型如图1所示。
|
图 1 二维模型 Fig. 1 2D model |
从左至右孔依次标记为①②③④⑤，盖板如图所示标记为①②③④，挡板与盖板延长线的夹角为45o，船长为4.356 m。
分别对光体、孔集中分布在平行中体首部、集中在平行中体中间、集中在平行中体尾部和均匀分布在平行中体上5种工况进行研究。各工况具体如表1所示。
|
表 1 工况 Tab.1 Conditions |
大涡模拟（LES）主要利用滤波的方法将涡分为大尺度涡和小尺度涡，然后分别进行数值模拟分析。大涡模拟控制方程为连续性方程以及N-S方程进行过滤之后的方程：
| $\frac{{\partial \rho }}{{\partial t}} + \frac{\partial }{{\partial {x_i}}}\left( {\rho \overline {{\mu _i}} } \right) = 0\text{，}$ | (1) |
| $\frac{\partial }{{\partial t}}\left( {\rho \overline {{\mu _i}} } \right) + \frac{\partial }{{\partial {x_j}}}\left( {\rho \overline {{\mu _i}} \overline {{\mu _j}} } \right) = \frac{\partial }{{\partial {x_j}}}\left( {\mu \frac{{\partial {\sigma _{ij}}}}{{\partial {x_j}}}} \right) - \frac{{\partial \bar \rho }}{{\partial {x_i}}} - \frac{{\partial {\tau _{ij}}}}{{\partial {x_i}}}\text{。}$ | (2) |
方程右边含有不封闭项
本文采用Smagorinsky涡粘模型来模化亚格子应力，其张量的偏量部分为：
| ${\tau _{ij}} = - 2{v_i}{S_{ij}} = - 2C_s^2{\Delta ^2}\left| {\bar S} \right|\overline {{S_{ij}}} \text{，}$ | (3) |
| $\left| {\bar S} \right| = \sqrt {2{S_{ij}}} \text{。}$ | (4) |
式中：vi为涡黏系数Cs为Smagorinsky常数；∆为滤波尺度，通常取0.1～0.23；Sij为经过滤波后的速度变形张量，定义如下：
| $\overline {{S_{ij}}} = \frac{1}{2}\left( {\frac{{\partial \overline {{\mu _i}} }}{{\partial {x_j}}} + \frac{{\partial \overline {{\mu _j}} }}{{\partial \overline {{x_i}} }}} \right)\text{。}$ | (5) |
整个流域采用非均匀的结构化网格，对于靠近艇体和空腔的挡板和盖板处进行了网格加密，边界层厚度约为0.8 mm。5 种工况的网格量约为1 650 000，y+最大值大约为0.6，计算域上游取3倍艇长，下游取9倍艇长，宽度为3倍艇长。
|
图 2 计算域 Fig. 2 Computational domain |
本文对空腔数值模拟的边界条件如下：
1）入口条件：入口处设置为速度入口，速度大小均为4 m/s，速度方向垂直于入口边界，为x轴正方向。
2）出口条件：出口处设置为压力出口。
3）壁面条件：艇体主体、水槽、挡板、盖板均设置为无滑移壁面条件。
4）对称条件：流域顶部和流域下部均设置为对称边界条件。
5）内部条件：5 个孔和内部流场均设置为内部条件。
2.4 数值离散采用大涡模拟计算方法，压力项采用二阶隐式格式求解，时间项采用二阶隐式，动量项采用中心差分求解，时间步长为0.000 25，每一时间步最大迭代次数为15，迭代收敛标准为0.000 001。使用收敛性较好的SIMPLEC求解压力耦合方程。
3 计算结果3.1 各工况下各部位阻力将各工况各部位阻力求解结果如表2所示。
|
|
表 2 各工况各部位阻力情况 Tab.2 The resistance of different parts on each condition |
由表2可知，工况1与开孔工况2～工况5对比发现，工况1的总阻力最小，而且阻力增量主要在艇体上；将4种开孔工况之间对比发现，工况3的总阻力相对最小，且艇体的阻力增量也最小，这跟此种开孔分布下的孔的内外流交换和压力分布都有关，此种工况下5个孔的流量最小，说明内外流交换最小，由内外流引起的主艇体的纵向压差最小，因而阻力增量最小。
3.2 各工况的压力分布与流线图
|
图 3 工况1流线图 Fig. 3 Streamlines of condition 1 |
|
图 4 工况2流线图 Fig. 4 Streamlines of condition 2 |
|
图 5 工况3流线图 Fig. 5 Streamlines of condition 3 |
|
图 6 工况4流线图 Fig. 6 Streamlines of condition 4 |
|
图 7 工况5流线图 Fig. 7 Streamlines of condition 5 |
由图3～图7可看出，工况1艇体附近几乎没有涡，流线沿艇体流过，并且越靠近艇体，流线越密集；工况2流线图在艇体首部出现涡，这些涡分布在开孔附近；工况3流线图上艇体中间部位出现涡，并且这些涡出现在开孔附近；工况4流线图上艇体尾部出现涡，并且这些涡出现在开孔附近；工况5流线图上艇体空腔内各部位都出现涡，这些涡的位置与开孔位置吻合。由于涡分布在孔的附近，即涡的形成位置与孔的位置有关，说明由于孔附近的内外流交换产生了涡，当水从一个孔流入后，在空腔内流动，受到空腔壁的阻挡改变方向，做螺旋运动，最终由其他孔流出，进而形成涡。
|
图 8 各工况压力分布曲线 Fig. 8 Pressure distribution curve of each condition |
各部位压力分布如图8所示。由图可知，各工况在开孔部位压力变化较为剧烈，并引起首尾部部分压力分布发生改变。
3.3 各工况下5个开孔流量流量定义流入为负，流出为正，各工况下5个开孔流量如表3所示。
|
|
表 3 各工况各孔流交换情况 Tab.3 Flow exchange of each hole on each condition |
由表3可知，工况2，水从孔③④⑤流入，从孔①②流出；工况3，水从孔③④⑤流入，从孔①②流出，但相比工况2，此种工况下5个孔的流量都较小，说明孔附近流交换小；工况4，水从孔①②③流入，从孔④⑤流出；工况5，水从孔②③④流入，从①⑤流出。综合这 4 种开孔分布工况可知，工况3的孔的流量最小，即此种工况下孔的内外流交换最小。
3.4 频谱分析
|
图 9 工况1的频谱特性曲线 Fig. 9 Spectral characteristic curve of condition 1 |
|
图 10 工况2的频谱特性曲线 Fig. 10 Spectral characteristic curve of condition 2 |
|
图 11 工况3的频谱特性曲线 Fig. 11 Spectral characteristic curve of condition 3 |
|
图 12 工况4频谱特性曲线 Fig. 12 Spectral characteristic curve of condition 4 |
|
图 13 工况5的频谱特性曲线 Fig. 13 Spectral characteristic curve of condition 5 |
结合图3～图8及表2并将光体工况与开孔工况对比知，开孔后阻力增加，而且阻力增量主要在主艇体上，这与开孔后孔附近流动交换有关。开孔后，水会从一部分孔流入空腔，从其他孔流出，从而影响首尾部压力分布，进而影响主艇体压差阻力，使得开孔后总阻力比光体大。
结合表2表3及流线图可知，4 种开孔工况下，孔集中分布在中部时，总阻力最小，而且相比其他工况主要是艇体上阻力较小。由于水按照顺压流动，孔开在首部时，水从后部孔流入空腔，从靠近首部孔流出，与首部来流交汇阻碍了首部流动速度，从而使首部压力增大，艇体首尾压差增大，艇体压差阻力增大，而摩擦阻力几乎不变，从而主艇体总阻力增大；当孔开在尾部时，距离尾部较近，水从远离尾部的孔流入，从靠近尾部孔流出，从而使得尾部水流动加快，尾部压力相应降低，而首部压力基本不变，故首尾压差增大，主艇体压差阻力增大，总阻力增大；当孔开在中部时，由于距离首尾均较远，对于首尾压力影响较小，所以阻力增量也最小；当孔均匀分布时，水从中间几个孔流入，从靠近首尾的 2 个孔流出，分别阻塞首部流动，加速尾部流动，使得首部压力增大，尾部压力减小，首尾压差明显增大，主艇体压差阻力增大，总阻力增加最为显著。
由图9～图13可知，工况1阻力波动幅度较大频带较窄，约为f=13.6 Hz，阻力幅值峰值为10.1 N；开孔以后，阻力波动幅度较大的频带明显增加，波动幅值相比于光体增加明显。主要因为开孔引起的空腔内外流交换明显；工况2下阻力波动幅值较大频率范围10～25 Hz，阻力幅值峰值为15.6 N；工况3下阻力波动幅值较大频率范围15～25 Hz，阻力幅值峰值为11.3 N；工况4下峰值频率为13～30 Hz，阻力幅值峰值为13.5 N，开孔后峰值处频率向高频率值移动，并且开孔后阻力幅值均比光体阻力幅值大，工况5工况阻力波动幅值较大频率范围10～30 Hz，在这个频率区间内没有明显的峰值频率，最大波动幅值约为15 N，主要因为开孔引起大范围的内外流的交换。4 种开孔分布之间对比，发现开孔集中分布在平行中体中间时阻力幅值最小，而且波动最小，说明此种开孔分布对阻力频谱特性影响较小。
5 结 语本文采用大涡模拟的数值分析方法，以Suboff艇型为对象，研究平行中体上开孔分布对于空腔流动特性的影响，研究了光体、开孔集中分布在平行中体左边、开孔集中分布在平行中体右边、开孔集中分布在平行中体中间、开孔均匀分布在平行中体上 5 种工况下的压力分布、阻力情况、流交换情况以及阻力的频谱特性。得出如下结论：
1）开孔后阻力相较于光体总阻力增大，而且当开孔位于平行中体中间时总阻力的增量最大，而且阻力增量是由于开孔附近流交换引起的首尾压差变化，进而使压差阻力增大，对于真实的潜艇设计时，流水孔设计位置因与光体压力分布结合起来，尽量避免内外流的交换。
2）开孔引起的阻力增加主要是主艇体阻力增加，而非开孔部位的阻力增加，因此潜艇设计时注意开孔与主艇体的配合。
3）开孔分布对艇体水动力噪声有较大影响，不仅增加了噪声强度，而且改变了峰值频率的分布。
| [1] | 曹志敏, 杜阳华. 潜艇武器系统综合仿真系统及关键技术[J].计算机仿真, 2007, 24(6): 31–37. |
| [2] | 赖焕新, 周邵萍, 苏永升, 等. 空腔流动的大涡模拟及气动噪声控制[J].工程热物理学报, 2008(2): 228–232. |
| [3] | 邬明. 带流水孔潜艇流噪声特性分析[J].四川兵工学报, 2015(10): 27–31. |
| [4] | 陈灿, 吴方良, 李环, 等. 不可压缩空腔流振荡模式和声学特性研究[J].水动力学研究与进展(A辑), 2015(3): 272–278. |
| [5] | 孟生, 张宇文, 王艳峰. 潜艇流水孔流噪声特性研究[J].计算机仿真, 2011(7): 57–62. |
| [6] | 张楠, 沈泓萃, 姚惠之, 等. 带流水孔潜体流场数值模拟[J].船舶力学, 2004(1): 1–11. |
| [7] | ROSSITER J E. Wind tunnel experiments on the flow over rectangular cavities at subsonic and transonic speeds[R]. London, UK:Aeronautical Research Council, 1966. http://iopscience.iop.org/article/10.1088/0022-3727/47/50/502002/pdf |
| [8] | 涂海文, 孙江龙. 基于CFD的潜艇阻力及流场数值计算[J].舰船科学技术, 2012, 34(3): 19–25. |
| [9] | 李霖, 张志国, 王先洲, 等. 低雷诺数圆柱绕流的大涡模拟分析[J].舰船科学技术, 2013, 35(1): 22–26. |