船舶阻力是船舶快速性研究的重要内容之一。兴波阻力作为船舶总阻力的重要组成部分，对船型变化十分敏感，适当改变船型能有效降低兴波阻力。以兴波阻力为优化目标，可以对船体型线进行优化。作用在船舶的水动力都是通过船体湿表面计算获得的，较大航行姿态变化会影响水动力计算结果，纵倾较大时还会影响船舶性能。如浅水中航行的船舶，为避免发生触底现象，纵倾和升沉限制了船舶的最高航速。

Michell提出了兴波阻力的理论计算公式，之后国内外许多学者投入到理论计算兴波阻力的研究中。Dawson^[1]提出了用Rankine源代替Kelvin源求解船行波问题，取得了很好的应用效果。但是Dawson忽略了自由面条件中的非线性项，将自由面条件从真实的位置转移到了未扰动的静水面上，未考虑水面以上船体的影响。Jensen等^[2-5]在此基础上均进行了有关研究。Wyatt^[6]提出了一种非线性迭代格式，将波高与源强同时进行求解，从而降低自由面条件中速度势导数的阶数，但同时也增大了未知数数目。王中等^[7-8]采用Wyatt迭代格式，计算了Wigley、S60及考虑升沉和纵倾的DTMB5415，且对非线性兴波阻力计算过程中程序编制，网格配置，阻力计算等诸多环节进行了探讨。此外，Tarafder^[9]对自由面条件进行高阶摄动展开计算了非线性船行波问题，不足之处在于自由面条件仍在静水面满足，He^[10]基于线性理论计算船舶兴波阻力，并考虑了水深对兴波阻力影响。

本文基于三维势流理，采用Rankine高阶面元法，通过定常迭代求解非线性船行波问题。其中自由面条件迭代格式在文献[11]的基础上稍做改进，对当前迭代步仅保留线性项。为了提高计算效率，采用一种基于不完全LU分解的预条件迭代法求解方程。姿态计算采用一种迭代格式。计算结果包括Wigley和S60兴波阻力、姿态、自由面波形等，通过与试验结果对比，验证了本文非线性方法具有较高的准确性，讨论了非线性因素对船行波问题的影响。为进一步探索有限水深船行波问题，本文以Wigley船为对象，计算了在不同水深下兴波阻力和自由面波形，并探讨分析了水深对兴波阻力及自由面波形的影响。

1 数学模型和数值计算方法 1.1 船行波的数学模型

假定船舶以航速U沿x轴正向运动，O-xyz为随船平动坐标系，如图 1。原点O位于船舯静水面上，xoy平面平行于未扰动的静水面，x轴指向船艏，y轴指向左舷，z轴垂直向上。

假定流体无粘不可压，流动无旋，流场中的速度势满足以下方程和边界条件。

在流场中满足Laplace方程如下

$ \left\{ \begin{align} &{{\Delta }^{2}}\varphi =0\text{ } \\ &Z\leqslant\zeta \left( x, y \right) \\ \end{align} \right. $

(1)

在自由面上，满足运动学边界条件:

$ \left\{ \begin{align} &{{\varphi }_{x}}{{\zeta }_{x}}+{{\varphi }_{y}}{{\zeta }_{y}}-{{\varphi }_{z}}=0 \\ &\text{ }Z=\zeta \left( x, y \right) \\ \end{align} \right. $

(2)

以及动力学边界条件:

$ g\zeta +\frac{1}{2}(\nabla \varphi \cdot \nabla \varphi -{{U}^{2}})=0, Z=\zeta \left( x, y \right) $

(3)

船体湿表面，满足不可穿透条件

$ \Delta \varphi \cdot n=0 $

(4)

无穷远前方的辐射条件

$ \varphi =-Ux $

(5)

由伯努利方程计算船体表面动压力

$ p=\frac{1}{2}\rho (\nabla \varphi \cdot \nabla \varphi -{{U}^{2}}) $

(6)

兴波阻力系数由压力沿x方向积分并无因次化

$ {{C}_{w}}=\frac{-\iint\limits_{{{\overline{S}}_{b}}}{p{{n}_{x}}\text{d}s}}{0.5\rho {{U}^{2}}{{\overline{S}}_{b}}} $

(7)

式中：ϕ为速度势，ζ为波高，下标x、y表示偏导数，n=[n_x, n_y, n_z]为物面内法向量，p为计算船体表面的动压力，C_w为兴波阻力系数，S船体实际湿表面，S_b为船体平均湿表面。

1.2 非线性自由面条件的处理

假定速度势ϕ可以分解为均匀来流速度势-Ux和扰动势φ，即φ=-Ux+φ，将其代入由(2)、(3)组成的联合自ϕ=Ux+φ由面运动学和动力学边界条件中：

$ \begin{align} &{{U}^{2}}{{\varphi }_{xx}}+g{{\varphi }_{z}}=2U{{\varphi }_{x}}{{\varphi }_{xx}}+\text{ }2U{{\varphi }_{y}}{{\varphi }_{xy}}+U{{\varphi }_{z}}{{\varphi }_{xz}}- \\ &{{\varphi }_{x}}{{\varphi }_{x}}{{\varphi }_{xx}}-\text{ }2{{\varphi }_{x}}{{\varphi }_{y}}{{\varphi }_{xy}}-{{\varphi }_{y}}{{\varphi }_{y}}{{\varphi }_{yy}}-\text{ }{{\varphi }_{y}}{{\varphi }_{z}}{{\varphi }_{zy}}-{{\varphi }_{x}}{{\varphi }_{z}}{{\varphi }_{xz}} \\ \end{align} $

(8)

求解非线性兴波问题的一个难点就是处理非线性自由面边界条件，本文在求解中将自由面条件改写为如下迭代格式：

$ \begin{align} &{{U}^{2}}{{\varphi }_{xx}^{n}}+g{{\varphi }_{z}^{n}}=2U{{\varphi }_{x}^{n-1}}~{{\varphi }_{xx}^{n-1}}+2U{{\varphi }_{y}^{n-1}}{{\varphi }_{y}^{n-1}}+ \\ &\text{ }U{{\varphi }_{z}^{n-1}}{{\varphi }_{xz}^{n-1}}-\text{ }{{\varphi }_{x}^{n-1}}{{\varphi }_{x}^{n-1}}{{\varphi }_{xx}^{n-1}}-\text{ }2{{\varphi }_{x}^{n-1}}{{\varphi }_{y}^{n-1}}{{\varphi }_{xy}^{n-1}}-\text{ } \\ &{{\varphi }_{y}^{n-1}}{{\varphi }_{y}^{n-1}}{{\varphi }_{yy}^{n-1}}-\text{ }{{\varphi }_{y}^{n-1}}{{\varphi }_{z}^{n-1}}{{\varphi }_{zy}^{n-1}}-\text{ }{{\varphi }_{x}^{n-1}}{{\varphi }_{z}^{n-1}}{{\varphi }_{xz}^{n-1}} \\ \end{align} $

(9)

1.3 纵倾和升沉计算

本文计算船舶纵倾和升沉采用如下迭代方法，最终能得到船舶的真实平衡姿态。首先计算两次迭代间升沉和纵倾的变化值Δh和Δθ，定义下沉和艉倾为正，可表示为

$ \left\{ \begin{align} &\Delta h{{c}_{33}}+\Delta \theta {{c}_{35}}=\text{ }\iint\limits_{{{S}_{b}}}{p{{n}_{z}}\text{d}s}-\text{ }\iint\limits_{{{S}_{b}}}{\rho gz{{n}_{z}}\text{d}s}+\text{ }\iint\limits_{{{\overline{S}}_{b}}}{\rho gz{{n}_{z}}\text{d}s} \\ &\text{ }\Delta \theta {{c}_{55}}+\Delta h{{c}_{53}}=\iint\limits_{{{S}_{b}}}{p(z{{n}_{x}}-x{{n}_{z}})\text{d}s}~-\text{ } \\ &\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \iint\limits_{{{S}_{b}}}{\rho gz(z{{n}_{x}}-x{{n}_{z}})\text{d}s}+\text{ }\iint\limits_{{{\overline{S}}_{b}}}{\rho gz(z{{n}_{x}}-x{{n}_{z}})\text{d}s} \\ \end{align} \right. $

(10)

式中：c_ij为船舶静水回复力系数，具体表达式参考文献[12]。计算得到纵倾和深沉变化后，叠加前一次迭代对应的结果，即可得到更新后的纵倾和升沉:

$ \left\{ \begin{align} &{{\theta }^{n}}={{\theta }^{n-1}}+\Delta \theta \text{ } \\ &{{h}^{n}}={{h}^{n-1}}+\Delta h \\ \end{align} \right. $

(11)

当第n次和n-1次的姿态小于给定误差时认为结果收敛。

1.4 高阶面元法的数值实现

相比于常值面元法，高阶面元法能在较小的面元离散数目下获得更精确的解。本文采用如图 2所示的9节点等参数高阶单元离散船体和自由面。单元内任意点坐标和源强可表示为

$ \left\{ \begin{align} &\left[x, y, z \right]=\underset{j=1}{\overset{9}{\mathop{\sum }}}\, {{N}_{j}}\left( \xi, \eta \right)\text{ }[{{x}_{j}}, {{y}_{j}}, {{z}_{j}}];\text{ } \\ &\sigma =\underset{j=1}{\overset{9}{\mathop{\sum }}}\, {{N}_{j}}\left( \xi, \eta \right){{\sigma }_{j}} \\ \end{align} \right. $

(12)

式中j为单元内节点编号。

采用源分布法求解上述定解问题，场点p处的速度势可以表达为

$ {{\varphi }_{p}}=\sum\limits_{e=1}^{\text{Ne}}{\sum\limits_{j=1}^{9}{{{\sigma }_{j}}\int_{-1}^{1}{\int_{-1}^{1}{{{G}_{p, q(\xi, \eta )}}{{N}_{j}}\left( \xi, \eta \right)\left| \mathit{\boldsymbol{J}}\left( \xi, \eta \right) \right|\text{d}\xi \text{d}\eta }}}} $

(13)

式中：N_e为总面元数，q为源点，σ为源强，G_{p, q(ξ, η)}=1/r_{p, q(ξ, η)}，N_j(ξ, η)为形函数J(ξ, η)为Jacob矩阵，形函数和Jacob矩阵表示方法参考文献[13]。

将式(13)代入自由面条件和物面条件，边界条件中速度势的一阶导数通过对格林函数解析求导得到，二阶数通过对一阶导进行数值差分得到。自由面为贴体网格，数值差分沿x方向采用Dawson的四阶迎风格式，y方向采用中心差分格式，归纳得到按节点离散的边界积分方程:

$ \left[\begin{matrix} {{\mathit{\boldsymbol{H}}}_{bb}}&{{\mathit{\boldsymbol{H}}}_{bf}} \\ {{\mathit{\boldsymbol{H}}}_{fb}}&{{\mathit{\boldsymbol{H}}}_{ff}} \\ \end{matrix} \right]\left[\begin{matrix} {{\sigma }_{b}} \\ ~{{\sigma }_{f}} \\ \end{matrix} \right]=\left[\begin{matrix} {{\mathit{\boldsymbol{C}}}_{b}} \\ {{\mathit{\boldsymbol{C}}}_{f}} \\ \end{matrix} \right] $

(14)

式中：H为系数矩阵，C是右端项。

按节点离散的边界积分方程是一个线性的方程组，其系数矩阵维数较大，而且方程组的系数矩阵是稠密非对称的。采用传统的高斯消去法效率极低，本文采用saad^[14]推荐的GMRES(generalized minimum residual method)方法，该方法适用于大型非对称的线性方程组求解。直接迭代法求解线性方程组稳定性较差，采用预处理方法可以提高计算效率和迭代的稳定性。本文采用单机并行计算的方法组装系数矩阵，采用不完全LU分解的预处理方法求解线性方程组，当计算总结点数为5 584时，求解方程组时间仅为2.25 s，相比于高斯消元法，求解方程组时间为122 s，计算效率大大提高。

1.5 船舶湿表面网格自动生成

非线性计算的迭代过程中，自由面在不断变化，要求数值计算程序能够根据自由面变化重新划分船体湿表面网格。本文在计算中对基础型值数据进行三次B样条插值实现船舶湿表面网格值自动划分。自由面计算域自船艉向下游延伸约2倍船长，自由面自船艏向上游延伸0.5倍船长，船侧为1倍船长，为满足辐射条件，纵向差分采用Dawson四阶迎风格式，同时将源点后置0.3倍纵向网格间距^[15]，与控制点形成错位。自由面源点同时上至1倍纵向网格间距以去除奇异性，同时提高数值稳定性。

如图 3依次为Wigley和S60在Fr为0.328时第8次迭代开始时刻的贴波面船体面元网格划分，以及初始时刻自由面计算域网格划分。计算过程中，对Wigley船沿纵向和垂向离散31×17个节点，自由面沿纵向和横向离散95×25个节点；对S60船沿纵向和垂向离散49×15个节点，自由面沿纵向和横向离散123×31个节点。

2 无限水深算例结果和分析

基于上述理论，采用Fortran语言开发了数值计算程序，以Wigley和S60船为例，其中Wigley船主尺度为L为1, B为0.1, D为0.062 5，排水量为2.78 kg，重心为(0，0，0)S60船的主尺度为L为1, B为0.133，D为0.053，排水量为4.25 kg，重心为(-0.014 6，0，0)，转动惯量取0.25 L_pp，L_pp为垂线间长。计算其在不同航速下的兴波阻力、航行姿态和自由面波形等，并与相关实验数值结果进行对比分析。

图 4(a)、(b)分别为线性和非线性计算Wigley和S60兴波阻力系数与相关试验值^[16-17]的比较。兴波阻力的试验值是由测得的总阻力扣除摩擦阻力的k+1倍得到的，其中k表示形状因子。非线性计算兴波阻力系数包括两种方法：考虑姿态影响和不考虑姿态影响。前者对应自由模状态，后者对应固定模状态。从图中可以看出计算Wigley和S60非线性和线性兴波阻力系数随傅汝德数Fr的变化规律和实验值一致，非线性计算结果与自由模实验值吻合良好。当航速较低时(Fr＜0.25)姿态变化对兴波阻力系数几乎没有影响。随着航速的增加，姿态变化对兴波阻力系数的影响变大，相比于自由面条件的非线性项对兴波阻力系数的影响而言，姿态变化影响更大。

考虑到计算的效率，非线性计算S60船，总结点数为4 551，当Fr为0.2时，迭代6次后收敛，总耗时161 s; 当Fr增加时，迭代步数会相应增大，计算用时将更长，相比之下线性计算效率更高，仅需要10 s左右，其中求解方程组的时间约为3 s。因此，在最初的设计阶段，如果更关心兴波阻力曲线的趋势和拐点位置，选用线性理论效率更高。

图 5和图 6分别为线性和非线性计算Wigley和S60纵倾和升沉与试验值^[16-17]对比。从图中可以看出：当Fr＜0.25时，非线性计算和线性计算纵倾和升沉相差不大；当Fr＞0.25时非线性和线性计算纵倾和深沉有明显的差别。非线性计算纵倾和升沉更加精确。这是由于在计算纵倾和升沉的时候回复力系数需要在船体湿表面积分，而非线性计算相比线性计算考虑到船体湿表面变化。非线性计算纵倾和升沉与实验值吻合良好。

图 7和图 8为线性和非线性计算Wigley和S60船侧波形与文献[18]对比。线性计算波形在船舯以后与试验结果接近，船舯以前波长和波幅均与实验值相比有明显的差别，且波形向艉部移动。非线性计算结果试验值吻合良好，尤其是非线性计算有效地改善船艏波幅降低的现象，但是由于在船艏存在波浪飞溅、破碎现象，无法用势流理论计算，因此非线性计算船艏波形与实验值仍然有些差异。由非线性自由模和固定模计算结果比较可知，考虑姿态和不考虑姿态对船侧波形的影响不大。

图 9为Fr=0.316时非线性和线性计算S60纵切波波形图。从图中可以看出，线性预报结果与试验值^[18]相比有所差异，主要体现在三点：1)低估了船艏波幅；2)高估了船艉后方波高；3)波长略有不同。而非线性预报结果在精度上有明显提高，与试验基本吻合一致。姿态变化对波形影响不大，相比之下姿态变化对阻力有一定影响。

3 有限水深的船舶兴波问题

船舶在浅水水域中航行与在深水水域中航行相比，兴波阻力和自由面波形有显著的差别。本文以Wigley船为例，计算在不同水深下兴波阻力和自由面波形并讨论了水深变化对兴波阻力和船侧波形的影响。为了满足水底不可穿透条件，在数值计算过程中，将船体表面和自由面关于水底镜像。自由面计算域自船艉向下游延伸1.5倍船长，自船艏向上游延伸0.75倍船长，船侧延伸1.25倍船长。船体网格沿纵向和垂向离散为33×13个节点，自由面沿纵向和横向离散为77×29个节点。

图 10为Fr=0.316时线性和非线性计算兴波阻力系数与文献[19]数值计算结果^[19]对比。从图中可以看出本文线性计算兴波阻力与文献结果一致，非线性计算兴波阻力相比于线性计算兴波阻力略大。当水深接近临界值时(Fr_h=1)，兴波阻力急剧增大，Fr_h=Fr√(L/h)为水深傅汝德数。

图 11表示不同水深(h/T=1.6，2，3.2，∞)非线性计算兴波阻力系数随航速变化。图中可以看出，浅水对兴波阻力影响体现在：低航速区浅水计算兴波阻力相比深水计算值大，高航速区浅水计算兴波阻力比深水计算值小。对某一固定水深，兴波阻力系数曲线峰值出现在0.96≤Fr_h≤0.97，接近船舶临界速度区。对比不同水深兴波阻力系数曲线比较可知，随着水深减少兴波阻力系数曲线峰值变大。

图 12为Fr=0.316时非线性计算不同水深(h/T=1.6，2，3.2，∞)的船侧波高。由图可知，波长和波幅随着水深减少而增大，同时船艏艉处的波形变得更加陡峭。

图 13为Fr=0.316时不同水深下(H/T=1.6，2)自由面波形图。图(a)中可以明显得观察到横波和散波。随着H/T减小，凯尔文角逐渐增大，当H/T=1.6时，凯尔文角达到最大值90°，如图(b)，此时Fr_h=1。

4 结论

1) 本文所采用的计算非线性船行波的数值方法能够准确预报船舶兴波阻力、航行姿态和船侧波高，对不同船型、不同航速、不同水深均广泛适用。

2) 兴波阻力预报与实验值变化趋势一致。线性计算相比非线性方法更高效，对于文中全域3 000网格量问题，10 s即可完成单个航速计算。非线性计算与自由模实验值吻合良好。

3) 线性预报船侧波形与纵切波形与试验验值相比船艏波幅较小而船艉增大，且波长略有不同。非线性计算有效改善线性计算与试验值的差异，船侧波形与波切均与实验值吻合。

4) 姿态对阻力的影响与非线性项相比更大，而对波面升高贡献不大。

5) 浅水计算兴波阻力系数随着航速增加增大，在Fr_h接近1(0.96≤Fr_h≤0.97)时达到最大值，之后随航速增加兴波阻力系数逐渐减小。水深越浅，兴波阻力曲线峰值越大。

6) 浅水船行波计算中，随着水深减少，船侧波形的波长和波幅变大，且艏艉处波形更加陡峭。