0 引　言

随着船舶大型化快速的发展，人工河道、天然河流以及船舶港口的建设速度没能跟上船舶大型化的发展速度，尤其是船舶航道的深度变化无法跟上船舶吃水深度的变化^[1−2]。和海洋中接近无限深的水域相比，内行河道的深度比较浅，船舶容易受到浅水效应的影响，浅水效应会对船舶的水动力特性以及航行的姿态产生一定影响，因此会极大增加船舶的控制难度^[3]。船舶在浅水区域航行的过程中，由于船体的底部和航道的底部之间的距离比较小，因此会导致航道内部的水流由三维流动转变成二维流动^[4]。航道内部水流运动模式的转变使得船体底部的水流速度加快，进而降低了船体底部的压力，最终使得船舶航行的阻力变大，更加容易下沉，并且船体底部压力的变化朝着船体尾部延伸，则会加剧船舶的纵倾^[5]。因此船舶的浅水效应除了会降低船舶的可操纵性，还会使得船舶出现擦伤、搁浅等事故，极大地威胁到船舶航行的安全^[6]。

1 数值计算方法 1.1 网格划分及控制方程

FLUENTTM求解器可生成输入模型所对应的计算网格，并且针对复杂的模型以及流体可对相应部位的网格进行加密处理，进而得到细化程度更高的网格，同时在对模型进行网格划分的过程中可采用自动适应处理，这样能降低模型网格生成所需要的时间^[7]。其细化的数学模型为：

$ {2^N} = \frac{{\Delta {X_{init}}}}{{\Delta {X_{targ }}}}\text{，} $

(1)

$ N = \log \left( {\frac{{\Delta {X_{init}}}}{{\Delta {X_{targ }}}}} \right) \text{。} $

(2)

式中：N为网格细化的次数；∆X_init为常规网格尺寸；∆X_targ为细化之后的最小网格尺寸。在对模型边界上的网格进行细化加密过程中，FLUENTTM求解器会根据模型的尺寸以及相应的流场速度计算出一个估计值，利用这个估计值对模型边界上的网格尺寸进行解算，该估计值的计算方法如下式：

$ y = \max \left( {{y_{\min }},\min \left( {30 + \frac{{\left( {{Re} - 1{e^6}} \right) \times 270}}{{1{e^9}}}} \right),{y_{\max }}} \right)\text{。} $

(3)

式中：y_max=300，y_min=30，Re为雷诺常数。式（3）一般用于近壁面，并且式（3）中的y值确定之后，则可以对模型边界的厚度进行解算，计算方法如下式：

$ {Y_{wall}} = 6{\left( {\frac{{{V_{ref}}}}{V}} \right)^{ - \frac{7}{8}}}{\left( {\frac{{{L_{ref}}}}{2}} \right)^{\frac{1}{8}}}{y^ + }\text{。} $

(4)

式中：V_ref为特征速度；L_ref为特征长度。在对模型进行网格划分的时候，需设置边界层，计算因子为影响边界层的最主要因素，通常设置计算因子为1.2。对流动问题进行解算的时候需满足质量守恒定律，即：

$ \frac{{{\mathrm{D}}\rho }}{{{\mathrm{D}}t}} + \rho div\left( V \right) = 0\text{。} $

(5)

当需解算的流体是不可压缩的，那么该流体的密度为一个常值，并且满足式（6），则式（5）可化为式（7）。

$ \frac{{{\mathrm{D}}\rho }}{{{\mathrm{D}}t}} = 0\text{，} $

(6)

$ \frac{{\partial u}}{{\partial x}} + \frac{{\partial v}}{{\partial y}} + \frac{{\partial w}}{{\partial z}} = 0\text{。} $

(7)

粘性流体的N-S方程如式（8）和式（9）所示。当流体是粘性不可压缩时，则S-N方程可转变为式（10）。

$ \nabla \cdot v = 0\text{，} $

(8)

$ \frac{{{\mathrm{D}}v}}{{{\mathrm{D}}t}} = {F_b} - \frac{1}{\rho }\nabla p + \upsilon \Delta v + \frac{1}{3}\upsilon \nabla \left( {\nabla \cdot v} \right)\text{，} $

(9)

$ \frac{{\partial {u_i}}}{{\partial t}} + \frac{\partial }{{\partial {x_j}}}\left( {{u_i}{u_j}} \right) = - \frac{1}{\rho }\frac{{\partial p}}{{\partial {x_i}}} + \upsilon \frac{\partial }{{\partial {x_j}}}\left( {\frac{{\partial {u_i}}}{{\partial {x_j}}} + \frac{{\partial {u_j}}}{{\partial {x_i}}}} \right)\text{。} $

(10)

1.2 湍流数值模型

船舶的船尾形状以及雷诺数等因素会对船舶的航行造成影响，并且会呈现出高复杂性的三维流动特性。因此在对流体进行模拟的过程中需将雷诺应力和湍流联系起来，如下式：

$ - \overline {\rho {{u'}_i}{{u'}_j}} = \left| {{\mu _t}\left( {\frac{{\partial {u_i}}}{{\partial {x_j}}} + \frac{{\partial {u_j}}}{{\partial {x_i}}}} \right) - \frac{2}{3}} \right.\left( {\rho k + \frac{{\partial {u_i}}}{{\partial {x_i}}}} \right){\delta _{ij}}\text{。} $

(11)

式中，δ_ij为Kronecker符号，如式（12）所示；k为湍动能，其计算方法如下式：

$ {\delta _{ij}} = \left\{ {\begin{array}{*{20}{c}} {1}，{i = j}，\\ {0}，{i \ne j}。\end{array}} \right. $

(12)

$ k = \frac{{\overline {{{u'}_i}{{u'}_j}} }}{2}\text{。} $

(13)

从式（11）中可看出，涡粘模型将雷诺应力的解算问题转换成湍动粘度的问题。通常采用两方程模型对湍动粘度进行求解，在船舶水动力特性模拟分析中最常用的两方程模型为k-ε模型以及k-ω模型。k-ε模型中引入了湍流耗散率和湍动能2个量对湍动粘度进行求解，其中湍流耗散率的计算方法如式（14）所示。

$ \varepsilon = \frac{\mu }{\rho }\overline {\left( {\frac{{\partial {{u'}_i}}}{{\partial {x_k}}}} \right)\left( {\frac{{{{u'}_j}}}{{\partial {x_k}}}} \right)} \text{。} $

(14)

湍动能k、湍流耗散率ε和湍动粘度μ_t之间的关系可用下式表示，其中，C_μ为经验常数。

$ {\mu _t} = \rho {C_\mu }\frac{{{k^2}}}{\varepsilon }\text{。} $

(15)

k-ω模型中引入了比耗散率和湍动能对湍动粘度进行解算。k-ω模型不但弥补了k-ε模型在解算低雷诺数流动过程中的失真问题，还可对逆压梯度的流动分离现象进行预报。湍流耗散率、湍动能以及比耗散率三者之间的数学关系如下式所示，式中，C为流体特征系数。

$ \omega = \frac{{\varepsilon C}}{k}\text{，} $

(16)

因为流体存在一定粘性，在湍流的影响下，壁面上会出现一层边界层，并且壁面上湍流切向速度的变化比较快，因此湍流模型的解算精确取决于壁面流动的捕捉精度。通常情况下，近壁面流动采用壁面函数法计算，壁面边界层内的位置以及速度的计算式为：

$ {y^ + } = \frac{{\Delta y\rho {u_t}}}{\mu } = \frac{{\Delta y}}{v}\sqrt {\frac{{{\tau _w}}}{\rho }} \text{。} $

(17)

$ {u^ + } = \frac{U}{{{\mu _t}}}\text{。} $

(18)

图1所示为湍流边界层不同位置的流动速度变化情况。可知，随着位置的增大，湍流速度先缓慢增加，然后快速增加，最后线性增大。

2 螺旋桨特性分析

本文使用重叠网格法对船舶螺旋浆的敞水性能进行仿真分析，在仿真过程中，将解算区域设置成圆柱形，并且将船舶螺旋桨的推进系数设置为0.1～0.7，并且螺旋桨的旋转速度为恒定的9.9 r/min，为了能提升解算精度，本文将螺旋桨单位时间内的旋转角度设置成2°，螺旋桨推力和效率的计算式为：

$ {K_T} = \frac{T}{{\rho {n^2}{D^4}}}\text{，} $

(19)

$ {\eta _0} = \frac{{{K_T}}}{{{K_Q}}} \cdot \frac{J}{{2{\text{π}} }}\text{。} $

(20)

图2所示为仿真得到的船舶螺旋桨效率随船舶螺旋桨推进系数曲线。可知，船舶螺旋桨效率先增大后减小，并且在船舶螺旋桨推进系数为60%的时候，达到最大值。

船舶螺旋桨推力随船舶螺旋桨推进系数的变化曲线如图3所示。可以看出，随着船舶螺旋桨推进系数的增大，船舶螺旋桨的推力逐渐降低。

在开放水域中对船舶的拖航实验进行数值仿真分析，船舶螺旋桨的摩擦阻力系数、总阻力系数以及压阻系数的计算式为：

$ {C_f} = \frac{{{R_f}}}{{0.5\rho {V^2}S}}\text{，} $

(21)

$ {C_t} = \frac{{{R_t}}}{{0.5\rho {V^2}S}}\text{，} $

(22)

$ {C_p} = \frac{{{R_p}}}{{0.5\rho {V^2}S}}\text{。} $

(23)

仿真得到的船舶螺旋桨升沉值随佛汝德数的变化曲线如图4所示。可以看出，随着佛汝德数的增大，船舶螺旋桨的升沉值也逐渐变大。

3 船体螺旋桨耦合水动力分析

本文选择KVLCC2船舶的桨-舵系统作为研究对象，并且构建的船体模型的缩小比例为46.426，解算区域的水流速度设置成和KVLCC2船舶航行速度一致，同时傅汝德数取值0.142，船舶螺旋桨的转速设置为9.9 r/min，在该条件下对船舶螺旋桨的水动力性能进行探究，并分析其变化规律。为了能模拟水流的斜向流动，本文在流场的入口处设置了轴向和纵向速度，其计算式为：

$ {V_x} = - V \times \cos \beta \text{，} $

(24)

$ {V_y} = V \times \sin \beta \text{。} $

(25)

在舵角和斜流角均为0^o的工况下，等待解算结束之后，得到的不同频率下扭矩和压强的曲线分别如图5和图6所示。可以看出，船舶螺旋桨的脉动频率是一样的，其扭矩和压强均在频率为39.76 Hz的地方出现最大值。

在吃水深比为23.5工况下，仿真得到的压力随时间变化曲线以及功率谱密度随频率的变化曲线分别如图7和图8所示。可知，船舶螺旋桨受到压力在2000～3000 Pa之间；随着频率的增大，功率谱密度先快速增大然后又迅速的降低。

4 结　语

由于船舶属于高耗能的运输工具，因此研究船舶螺旋桨耦合系统的水动力性能对降低船舶能源消耗有着重要意义。船舶螺旋桨在旋转过程中，其后方会出现水流旋转的现象，为了能提升船舶螺旋桨的效率并收回一定能量，本文基于FLUENTTM求解器对船舶螺旋桨耦合系统进行分析，这有助于推升螺旋桨的推进性能，并且对船舶周围流场起到一定改善作用。