0 引　言

由于海洋物产丰富，各国都把海洋资源的勘探和开发作为发展重点^[1]。但水下环境险恶复杂，潜水员水下作业十分危险。为了提高安全性，水下作业装置已经成为水下工作和研究的重要工具^[2]。水下作业装置在工作时，由于工作环境复杂，会在水下以不同的偏航和俯仰姿态工作，导致其自身受力点和周围流场的变化，进而影响其工作时的稳定性。因此，分析水下作业装置在不同偏航、俯仰姿态下的水动力性能非常重要。肖昌润等^[3]用动网格法对比了RNG $ k-\varepsilon $湍流模型和SST $ k-\omega $湍流模型的计算精度。结果表明，两者计算精度相差不大，但是RNG $ k-\varepsilon $湍流模型计算时间较短。易文彬等^[4]采用STAR-CCM+软件预测船模阻力。结果表明，棱柱层的层数和总厚度对阻力预测有很大影响。王兴茹^[5]研究了潜艇在偏航和俯仰状态下的涡旋。结果表明，随着偏航角的增大，流场中的涡流强度增大。Xavier等^[6]研究了潜艇在不受自由面或海底影响情况下的偏航角。潜艇周围的负压随着偏航角的增大而增大。张思聪^[7]在偏航角下对潜水器进行仿真，结果表明，随着偏航角的增大，潜水器迎流面一侧的船体所受压力越来越大，流速越来越慢。Zhang等^[8]使用雷诺平均（RANS）和标准$ k-\varepsilon $模型方法，研究了蝶形机器人的阻力性能问题，并用实验证明了数值模拟的准确性。Abedi等^[9]采用风洞试验和数值模拟相结合的方法对潜水器流场进行分析，得到了潜水器表面压力分布，Seo等^[10]利用PIV测量了多截面KVL-CC2标准船模的二维速度场，分析了其流动特性和速度场。Waston等^[11]利用循环水槽作为试验环境对航母舰面流场进行测量，结果表明，实验方法在评估航母舰面流场时具有较高的准确性，并且与数值模拟结果相吻合程度较高。

本文利用STAR-CCM+软件和循环水槽PIV试验，对不同偏航角和俯仰角的水下作业装置表面压力、涡旋和尾流场进行仿真和试验，探索不同偏航角和俯仰角对水下作业装置水动力特性的影响。

1 数值理论

对于不可压缩粘性流体，在忽略热尾流的情况下，质量守恒和动量守恒被用来描述水下工作装置流场中流体流动的基本规律^[12]。对黏性流场的模拟，采用RANS方程进行求解^[13]。

连续性方程：

$ \frac{{\partial {u_i}}}{{\partial {x_i}}} = 0。$

(1)

动量守恒方程：

$ {f_i} - \frac{1}{\rho }\frac{{\partial p}}{{\partial {x_i}}} + v{\nabla ^2}{u_i} = \frac{{\partial {u_i}}}{{\partial t}} + {u_j}\frac{{\partial {u_i}}}{{\partial {x_j}}}。$

(2)

时均值和脉动值是通过变量分解法从瞬时的物理量中得出的，在RANS方法中被采用，即：

$ \varPhi = \bar \varPhi + \varPhi ' ，$

(3)

$ \bar \varPhi = \frac{1}{{\Delta t}}\int_t^{t + \Delta t} {\varPhi \left( t \right)} {\mathrm{d}}t。$

(4)

式中：$\varPhi $为某一物理量的瞬时值；$\bar \varPhi $为时均值；$\varPhi '$为脉动值。

将式(1)和式(2)的速度和压强值用时均值和瞬时值代替，得到时均连续性方程和Reynolds时均动量方程。

$ \frac{{\partial {U_i}}}{{\partial {x_i}}} = 0。$

(5)

式中：$ {U}_{i} $为平均速度，$ {\mathrm{m}}/{\mathrm{s}} $。

动量方程（雷诺方程）：

$ \frac{{\partial {U_i}}}{{\partial t}} + {U_j}\frac{{\partial {U_i}}}{{\partial {x_j}}} = {F_i} - \frac{1}{\rho }\frac{{\partial P}}{{\partial {x_i}}} + v\frac{{{\partial ^2}{U_i}}}{{\partial {x_j}\partial {x_j}}} + \frac{1}{\rho }\frac{{\partial \left( { - \rho \overline {u_i^{'}u_j^{'}} } \right)}}{{\partial {x_j}}}。$

(6)

由于雷诺应力的存在，导致方程组不封闭，无法直接求解，因此引入RNG $ k-\varepsilon $湍流模型来封闭控制方程并进行计算。

湍流动能$ k $方程：

$ \begin{aligned}[b]\frac{\partial }{{\partial t}}\left( {\rho k} \right) + \frac{\partial }{{\partial {x_i}}}\left( {\rho k{u_i}} \right)=& \frac{\partial }{{\partial {x_j}}}\left( {{\alpha _k}{\mu _{eff}}\frac{{\partial k}}{{\partial {x_j}}}} \right) +\\& {G_k} + {G_b} - \rho \varepsilon - {Y_M} + {S_k}。\end{aligned}$

(7)

湍流耗散率$\varepsilon $方程：

$ \begin{aligned}[b]&\frac{\partial }{{\partial t}}\left( {\rho \varepsilon } \right) + \frac{\partial }{{\partial {x_i}}}\left( {\rho \varepsilon {u_i}} \right)= \frac{\partial }{{\partial {x_j}}}\left( {{\alpha _\varepsilon }{\mu _{eff}}\frac{{\partial \varepsilon }}{{\partial {x_j}}}} \right)+ \\ & {G_{1\varepsilon }}\frac{\varepsilon }{k}\left( {{G_k} + {G_{3\varepsilon }}{G_b}} \right) - {G_{2\varepsilon }}\rho \frac{{{\varepsilon ^2}}}{k} - {R_\varepsilon } + {S_\varepsilon }。\end{aligned}$

(8)

在计算中，假设流体不可压缩，即其密度保持不变，$ {G}_{b} $和$ {Y}_{M} $为0，并且在计算过程中不需要考虑源项的影响，即$ {S}_{k} $和$ {S}_{\varepsilon } $为0，式(7)和式(8)变为：

$ \frac{\partial }{{\partial t}}\left( {\rho k} \right) + \frac{\partial }{{\partial {x_i}}}\left( {\rho k{u_i}} \right) = \frac{\partial }{{\partial {x_j}}}\left( {{\alpha _k}{\mu _{eff}}\frac{{\partial k}}{{\partial {x_j}}}} \right) + {G_k} - \rho \varepsilon ，$

(9)

$ \begin{aligned}[b]&\frac{\partial }{{\partial t}}\left( {\rho \varepsilon } \right) + \frac{\partial }{{\partial {x_i}}}\left( {\rho \varepsilon {u_i}} \right)= \\& \frac{\partial }{{\partial {x_j}}}\left( {{\alpha _\varepsilon }{\mu _{eff}}\frac{{\partial \varepsilon }}{{\partial {x_j}}}} \right) + {G_{1\varepsilon }}\frac{\varepsilon }{k}{G_k} - {G_{2\varepsilon }}\rho \frac{{{\varepsilon ^2}}}{k} - {R_\varepsilon }。\end{aligned}$

(10)

式中：$ {\mu _{eff}} = \mu \;+\; {\mu _t} $；${\mu _t} = \rho {C_\mu }\;\dfrac{{{k^2}}}{\varepsilon }$；${G_k} = - \rho \overline {u_i^{'}u_j^{'}} \;\;\dfrac{{\partial {u_j}}}{{\partial {x_i}}}$；$ {R_\varepsilon } = \dfrac{{{C_\mu }\rho {\eta ^3}\left( {1 - \eta /{\eta _0}} \right)}}{{1 + \beta {\eta ^3}}}\dfrac{{{\varepsilon ^2}}}{k} $；$\eta = Sk/\varepsilon $；$S = \sqrt {2{S_{ij}}{S_{ij}}} $；${S_{ij}} =\dfrac{1}{2}\left( {\dfrac{{\partial {u_i}}}{{\partial {x_j}}} + \dfrac{{\partial {u_j}}}{{\partial {x_i}}}} \right)$。

方程中包含了一些可调的经验参数，可以使计算更加趋于合理。一般情况下，可以取值如下：$ {G_{1\varepsilon }} $=1.42，$ {G_{2\varepsilon }} $=1.68，$ {C}_{\mu } $=0.0845，${\sigma _k}$=1.0，${\sigma _\varepsilon }$=1.3，${\eta _0}$=4.38，$\beta $=0.012。

2 计算方法和网格无关性验证

首先，对全附体1∶24缩小模型下全附体SUBOFF潜艇的阻力值进行仿真，并与美国大卫泰勒船池实验中心公布的实验结果进行对比^[14]，结果如表1所示。

可知，综合6个航速，表中计算结果和实验结果最大误差为4.9%，平均误差为2.95%，这一误差在一般工程要求的范围内。

随后，将潜艇在10 kn时（雷诺数$ Re=2.2\times {10}^{7} $）中纵面上半边缘线与指挥台围壁之间的压力系数与实验值进行比较^[15]，用压力系数验证数值结果的准确性，对比结果基本一致。

如图1所示，SUBOFF潜艇中纵面上半边缘线和指挥台中部的压力系数分布与实验数据高度吻合。同时可以看出，压力系数在艇首、指挥台前端和潜艇尾部前端较高，这是由于该位置流体流速较低造成的。仿真结果与美国大卫泰勒船池实验中心公布的实验结果^[14]基本一致，说明了本文数值计算方法预测水下作业装置流场的可行性和准确性。

最后，验证了水下作业装置的网格无关性。在雷诺数$ Re = 1.56 \times {10^5} $工况下，划分了5组不同大小的网格，对流场敏感区域进行加密后计算出阻力系数，评估网格无关性。在不同网格数下，水下作业装置的阻力见表2，根据阻力系数公式可得到不同总网格下的阻力系数。其中$ {F}_{d} $是水下作业装置所受阻力值。

根据图2可知，水下作业装置在第3套网格后（305万）阻力系数趋于稳定。考虑到网格质量对计算结果影响较大，最终选择第4套网格（480万）作为计算水下作业装置流场的标准。

3 计算模型和网格划分 3.1 几何模型

图3(a)为水下作业装置的3D图，图3(b)为水下作业装置1∶16下的实物模型图，包括主体和机械臂两部分，主要基本参数见表3。

3.2 计算域和流体域

本文旨在研究水下作业装置在偏航和俯仰工况下对表面压力、速度以及涡量的影响，以了解其发展规律，因此，为了减小阻塞比对计算精度的影响，应该增加Y、Z方向的计算域尺寸，如图4(a)所示。在试验中，循环水槽截面为0.6 m（宽）×0.6 m（水深），工作段长度为2.0 m，作业装置模型最大阻塞比小于5%，可忽略阻塞效应。

3.3 网格划分

对设定好的计算域进行网格划分，网格划分的质量会对水下作业装置的精度产生直接影响^[5]。在STAR-CCM+软件中，依次选择表面重建、自动表面修复、切割体网格单元生成器、棱柱层网格生成器来提高流体求解精度^[16]。由于计算水下作业装置的表面压力需要能够精确求解近壁流体的求解器，因此可通过“层选项”设置增加棱柱层，以提高计算结果的精度。根据网格无关性的结果，可以得到水下作业装置表面网格示意图，如图5所示。

4 结果分析

水下作业装置水下运动有6个自由度。通常情况下，绕Z轴旋转运动叫做偏航运动，对应的夹角叫做偏航角$ {\alpha _\gamma } $；绕Y轴旋转运动叫做俯仰运动，对应的夹角叫做俯仰角$ {\alpha _P} $，如图6所示。在$ Re = 1.56 \times {10^5} $工况下，对不同偏航角$ {\alpha _\gamma } $和不同俯仰角$ {\alpha _P} $下的水下作业装置进行数值模拟计算。

4.1 不同偏航、俯仰角对压力的影响

水下作业装置上最大压力点的位置随着来流方向的不同而不同，图7为水下作业装置迎流侧和背流侧表面压力云图。可知，偏航角非0时，水下作业装置迎流侧的压力明显高于背流侧。从水下作业装置背流侧表面的压力图中可以看出，主体首部和机械臂的尾部附近出现明显的低压区，低压区的范围随着偏航角的增大而扩大。相反，从水下作业装置的迎流侧可以看出，主体的尾部和机械臂的前端都有高压区，在相同的偏航角下，由于迎流侧的压力远高于背流侧，水下作业装置会受到一个从迎流侧到背流侧的力。随着角度的增大，迎流侧高压区的范围不断扩大，并逐渐向尾部蔓延。

为了进一步观察水下作业装置压力分布，选取4个偏航角下Z=0平面的压力系数图进行对比，如图8所示。可知，水下作业装置的最大压力系数随着偏航角$ {\alpha _\gamma } $的变化而变化，最大压差随着偏航角的增大而逐渐增大，最大压差左侧的整体压力系数在增大。

水下作业装置上最大压力的位置随着俯仰角的不同而不同。为了观察不同俯仰角对水面压力的影响，对俯仰角$ {\alpha }_{P} $=0°、5°和10°的水下工作装置的流场进行了模拟计算。来流方向与俯仰角$ {\alpha }_{P} $的相对位置见图6。

水下作业装置在俯仰角$ {\alpha }_{P} $为0°、5°、10°下的作业装置表面压力图如图9所示。随着俯仰角的增加，水下作业装置首部的最大压力逐渐向下发展，首部上方的低压区明显增大，并逐渐向中部发展，最后与尾部的低压区相连，机械臂上方的低压区也明显增大。

为了进一步观察水下作业装置压力分布，选取了俯仰角$ {\alpha }_{P} $为0°、5°、10°三种情况下X=1截面的轮廓曲线压力系数图进行对比，如图10所示。可知，随着俯仰角的改变，最大压力差产生的位置和大小基本相同；在−60°～60°之间，随着俯仰角的增大，压力差受机械臂的影响也在增大。

4.2 不同偏航、俯仰角对涡旋的影响

利用刘超群^[17]的第二代涡识别准则$ \Omega $准则提取出不同偏航角$ {\alpha }_{\gamma } $的流场中的涡旋结构。如图11(a)所示，在偏航状态下，水下作业装置上会出现3类涡旋，第1类涡旋是由机械臂顶端引起的梢涡，随着偏航角变化，不断的偏向背流侧偏移；第2类涡旋是由机械臂和机身流体交汇区域的形成的交汇涡；第3类涡旋是当流体在机身表面流动时，部分流动会分离，形成涡旋结构。随着偏航角的增大，水下作业装置机械臂上的梢涡逐渐向背流侧偏移，主体上产生的表面分离涡逐渐增多。如图11(b)所示，在俯仰状态下，随着俯仰角的增加，由机械臂顶端引起的梢涡和交汇区域的形成的交汇涡基本不变，而艇体中间产生的涡旋逐渐增加。

4.3 不同偏航、俯仰角对尾流速度的影响

根据雷诺相似准则，在$ Re = 1.56 \times {10^5} $条件下对不同偏航角$ {\alpha }_{\gamma } $的和不同俯仰角$ {\alpha }_{p} $的水下作业装置进行PIV试验。在偏航状态下，取作业装置左侧尾部中间尾流为拍摄平面，当水下作业装置的偏航角发生变化时，会影响周围流场的速度，如图12所示。可以看出，在偏航运动过程中，水下作业装置尾部后会形成一条低速带，负方向时，随着角度的增加，速度阻塞区域减小，低速带逐渐减小，正方向时，随着角度的增加，速度阻塞区域增大，低速区范围逐渐扩大。

在俯仰状态下，取作业装置尾部中间尾流为拍摄平面，如图13所示，当水下作业装置的俯仰角发生变化时，尾流场变化剧烈，俯仰角在负方向时，随着角度的增大，尾流下部分速度减小，上部分速度增大；俯仰角在正方向时，随着角度的增加，尾流下部分速度增加，上部分速度减小。

5 结　语

通过对水下作业装置的分析可得出如下结论：数值模拟阻力和压力系数与试验结果吻合度较高，6种航速下阻力结果与实验结果最大误差为4.9%，平均误差为2.95%，可以满足工程所需要求。水下作业装置随着偏航和俯仰角的增大，其表面压力、涡旋和尾流速度的变化也越大。而表面压力差、涡旋和速度阻塞区域的增大，则会导致其周围湍流强度和增加，能量耗散加剧，从而对流体运动的稳定性产生较大影响，后续可以针对机械臂和作业装置的形状进行优化，或者通过增加附体等方法来增强水下作业装置的稳定性。