0 引　言

在船舶及海上工程机械领域，静水和波浪条件下的兴波运动是一种常见运动工况，此时船舶的运动是一种非线性运动，研究该运动过程的动力、阻力特性对于指导船舶设计有非常重要的意义。传统的船舶水动力性能研究以模型试验为主，这种研究方式存在多种缺点，比如模型设计和生产的周期长、试验成本高等。近年来，计算机和仿真技术获得了迅速的发展，基于该技术进行工程领域的试验和仿真成为了一种趋势。本文研究的重点是船舶兴波运动过程的非线性现象建模与仿真，介绍船舶动力学分析的基础理论，采用面元法对船舶的表面结构进行有限元划分，利用NAPA有限元求解器和Fluent仿真软件进行了船体的特性仿真。

1 船舶兴波运动非线性特性研究的理论基础

本文使用计算流体动力学理论进行船舶兴波运动的非线性研究，首先建立船舶兴波运动的坐标系如图1所示。

船舶兴波运动坐标系o-xyz中，oz方向垂直于海水平面竖直向上，海水的流速为 $\bar V = \nabla \phi $ ，且流体满足拉普拉斯方程，如下：

$ {\nabla ^2}\phi = 0，$

式中 $\phi $ 为海水的速度势。

船舶兴波运动分析的边界条件^[1]包括：

1）海水动力学边界条件

海水作为不可压缩液体，表面上方的气压为标准大气压且处处相同，可得海水的动力学边界条件如下式：

$ \frac{{\partial \phi }}{{\partial t}} + \frac{1}{2}\nabla {\phi ^{^2}} + g\delta = 0。$

式中： $g$ 为重力加速度； $\delta $ 为海水粘度系数。

2）船舶运动学边界条件

船舶的兴波运动是船舶在自由液面上速度和加速度从0开始不断提升，直至达到额定速度的过程，建立其运动方程为：

$ F = f\left[ {x(t),y(t),t} \right] - z\left( t \right) = 0 。$

其中： $x = x(t) ，y = y(t) ，z = z(t) ，$ 为船舶的航迹，求导可得船舶的运动学边界条件如下式：

$ \frac{{{\rm{d}}F}}{{{\rm{d}}t}} = \frac{{\partial F}}{{\partial t}}\frac{{{\rm{d}}z}}{{{\rm{d}}t}} + \frac{{\partial F}}{{\partial t}}\frac{{{\rm{d}}y}}{{{\rm{d}}t}} + \frac{{\partial F}}{{\partial t}}\frac{{{\rm{d}}x}}{{{\rm{d}}t}} = 0。$

3）扰流边界条件

扰流边界条件是指液体湍流特性对船舶造成的扰动力条件，由于船舶周围的液体流速随着距离的增加而不断发生衰减，建立扰流边界条件如下：

$ \mathop {\lim }\limits_{n \to \infty } \sqrt R \left( {\frac{{{\rm{d}}\phi }}{{{\rm{d}}t}} + \frac{1}{B}\frac{{{\rm{d}}\phi }}{{{\rm{d}}r}}} \right) = 0。$

式中：R为液体与船舶表面的距离；B为湍流强度。

4）RANS方程

船舶的兴波运动采用雷诺平均理论进行湍流求解，将湍流、速度量、位移量等进行不同坐标系下的分解，一方面能够降低复杂方程的计算量，提高流体动力学的求解速度；另一方面也能提高求解的精度，保证计算的准确性。

首先将船舶的速度量分解为 $u = {u_i} + {u_j}$ ，海水的速度势分解为 $\phi = {\phi _i} + {\phi _j}$ ，建立船舶兴波运动的RANS方程为：

$ \left( {\begin{array}{*{20}{c}} {\dfrac{{\partial {u_j}}}{{\partial t}} + \dfrac{\partial }{{\partial t}}(\rho {u_i}) = 0}，\\ {\dfrac{{\partial (\rho {u_i})}}{{\partial t}} + \dfrac{\partial }{{\partial t}}(\rho {u_i}{u_j}) = \dfrac{{\partial \rho }}{{\partial t}} + \mu \left[\dfrac{{\partial {u_i}}}{{\partial t}} + \dfrac{{{\partial ^2}{u_i}}}{{\partial {t^2}}}\right]} 。\end{array}} \right. $

式中 $ \rho $ 为海水的密度。

海水湍流模型的雷诺应力为：

$ \Delta {F_L} = \delta \left[\frac{{\partial {u_i}}}{{\partial {x_i}}} + \frac{{\partial {u_j}}}{{\partial {x_j}}}\right] - \frac{1}{3}pk{L_{ij}}。$

式中： ${L_{ij}}$ 为雷诺应力常数，其计算公式为

$ {L_{ij}} = \frac{{{C_e}\rho {\eta ^4}\left(1 - \dfrac{\eta }{{{\eta _0}}}\right)}}{{1 + b{\eta ^4}}}\frac{{{\varepsilon ^3}}}{k}。$

式中： ${C_e}$ 为湍流模型的动能系数； $\eta $ 为平均形变模数； ${\eta _0}$ 为初始模数； $\varepsilon $ 为流体对船舶的浮力系数。

2 基于面元法的船舶兴波运动非线性现象分析与仿真 2.1 面元法基本理论

不论是结构力学还是计算流体力学领域，有限元计算方法是最常用的一种方法，其核心是将复杂的力学模型和力学方程进行有限元分解，将复杂模型划分为一个个小单元，通过求解单元模型及节点方程，实现整个复杂系统的求解^[2]。

面元法是指将物体表面的光滑曲面等进行有限元化，将曲面离散为多个有限面积的面元，并建立该面元的流体特性方程，最后通过求解面元方程，获取船舶与海洋结构物的水动力特性。

在分析船舶兴波运动时，考虑到船体表面的复杂特性，首先将曲面按照船体结构模块进行划分，分为船首、船体、尾部等，船体又根据尺寸长度划分为N个，示意图如图2所示。

面元法的基础数值理论是格林函数法^[3]，对于任意的三维空间域 $\Omega $ ，存在

$ \iiint_\Omega \nabla \cdot \vec A{\rm{d}}V = \iint_S {\vec n} \cdot \vec A{\rm{d}}S，$

假设 $ \vec A = \phi \nabla G $ ，可得：

$ \iint_S G \cdot \frac{{\partial \phi }}{{\partial n}}{\rm{d}}S = \iiint \nabla \phi \cdot \nabla G{\rm{d}}V + \iiint_\Omega G{\nabla ^2}\phi {\rm{d}}V，$

两式相减可得：

$ \iint_S {\left( {\phi \frac{{\partial G}}{{\partial n}} - G\frac{{\partial \phi }}{{\partial n}}} \right)}{\rm{d}}S = \iiint_\Omega {\left( {{\phi ^2}G - G{\nabla ^2}\phi } \right)}{\rm{d}}V。$

该式为格林公式，当空间区域 $\varOmega$ 中满足：

$ {\nabla ^2}G = \phi {\nabla ^2} = 0 ，$

则可将格林公式化为：

$ \iint_\Omega {\left( {\phi \frac{{\partial G}}{{\partial n}} - G\frac{{\partial \phi }}{{\partial n}}} \right)}{\rm{d}}S = 0。$

在计算船舶兴波运动特性时，将 $ \phi $ 定义为流体的速度势，选择恰当的格林函数G，可得：

$ \iint_\Omega {\phi \left( Q \right){\nabla ^2}G\left( {P,Q} \right)}{\rm{d}}V = \phi \left( P \right) ，$

式中： $ P $ 为流域内的离散点^[4]； $ Q $ 为流域内的源点。

可知，速度势 $ \phi $ 在任一点的数值可以使用边界上的函数值和法向导数确定，定义S为空间区域 $\varOmega $ 的边界， $ \overrightarrow n $ 是S曲面的单位法线方向，可得：

$ \left\{ {\begin{array}{*{20}{l}} {\nabla \cdot \vec A = \phi \cdot \nabla G + {\nabla ^2}G} ，\\ {\vec n \cdot \vec A = \phi \cdot \dfrac{{\partial G}}{{\hat n}}} ，\end{array}} \right. $

船体的面元函数可得：

$ \iint_S \phi \cdot \frac{{\partial G}}{{\partial n}}{\rm{d}}S = \iiint_\Omega \nabla \phi \cdot \nabla G{\rm{d}}V + \iiint_\Omega {{\phi ^2}}G{\rm{d}}V 。$

图3为面元法的示意图。

2.2 船舶兴波运动的阻力建模和仿真求解流程

在进行船舶非线性兴波运动的数值分析时，使用NAPA求解器和有限元分析软件Fluent，并配合面元法函数模型进行流场仿真的简化^[5]。

船舶曲面上的点 $ \left( {{x_i},{y_i}} \right) $ ，其近似函数可用下式表示：

$ {h_{ij}}\left( {{x_i},{y_i}} \right) = \sum\limits_i^{} {\sum\limits_j^{} {{y_{ij}}} } h\left( x \right) 。$

利用有限元积分算法，可得船体曲面的阻力函数如下式：

$ {F_w} = \frac{{4{g^2}\rho }}{{\text{π} {v^2}}}\int_1^\infty {\frac{{{\lambda ^2}}}{{\left( {{\lambda ^2} - 1} \right)}}} \cdot \left[ {{I^2}(\lambda ) + {J^2}(\lambda )} \right]{\rm{d}}\lambda ，$

其中： $ g $ 为重力加速度； $ \rho $ 为海水的密度； $ v $ 为表面与海水的相对航行速度； $ \lambda $ 为海水波长； $ I(\lambda ) $ 和 $ J(\lambda ) $ 分别为波浪水平和垂向的扭矩。分别用下式计算：

$ I(\lambda ) = \int_0^T {{e^{\frac{g}{{{\nu ^2}}}{\lambda ^2}}}} \left[ {\cos \left( {\frac{g}{{{v^2}}}\lambda } \right)} \right]{\rm{d}}s ，$

$ J(\lambda ) = \int_0^T {{e^{{\lambda ^2}}}} \left[ {\sin \left( {\frac{g}{{{v^2}}} \cdot \lambda } \right)} \right]{\rm{d}}s ，$

简化可得：

$ {F_w} = \frac{{8\rho g}}{\pi } \cdot \frac{{{B^2}{T^2}}}{L} 。$

其中：B为船宽；L为船的长度。

船舶兴波运动的流体动力学仿真流程如图4所示。

2.3 船舶兴波运动非线性现象仿真研究

在进行船舶兴波运动仿真研究时，采用的船舶对象为标准的2 000 t集装箱运输船，该船的主尺度参数见表1。

为了提高船体模型的表面网格质量，采用数据插值和二次曲线拟合的方法进行有限元模型细化，具体如下：

1）插值计算

采用插值算法提高有限元的精度，在坐标系内给定n个数据点 $ pi\left( {{x_i},{y_i}} \right),i = 1,2,\cdots,n $ ，给定m个插值点： $ {x_1},{x_2},\cdots,{x_m} $ ，插值点对应的函数值为 $ {y_1},{y_2},\cdots,{y_m} $ ，在区间 $ \left[ {{x_k},{x_{k + 1}}} \right] $ 上可得插值函数^[6]为：

$ \left\{ \begin{gathered} {y_k} = f\left( {{x_k}} \right)，\hfill \\ {y_{k + 1}} = f\left( {{x_{k + 1}}} \right)，\hfill \\ \frac{{{\rm{d}}{y_k}}}{{{\rm{d}}t}} = {g_k}，\hfill \\ \frac{{{\rm{d}}{y_{k + 1}}}}{{{\rm{d}}t}} = {g_{k + 1}}。\hfill \\ \end{gathered} \right. $

式中 $ {g_k} $ 为曲线的曲率。

2）根据插值点确定有限元的型线

$ {f_0}\left( x \right) = a + b\left( {x - {x_k}} \right) + c{\left( {x - {x_k}} \right)^2} 。$

3）生成有限元模型

图5为船首位置的有限元模型示意图。

在Fluent中导入有限元模型和边界条件，得到图6船舶兴波运动的船首形变量仿真结果。

3 结　语

船舶兴波运动是一种非线性运动，该过程的兴波阻力是船舶设计时需要重点考虑的影响因素，本文结合计算流体力学和面元法，借助Fluent软件和NAPA求解器，进行了船舶兴波运动的数学建模和有限元仿真。