0 引　言

海运活动安全性和海上资源开采作业率要求的提高对船舶与海洋结构物在波浪中的生存周期提出了新的要求。三维势流理论因其计算效率高且结果准确的特点被广泛应用到船舶早期适航性设计中^[1]。当自由面格林函数被用于流体速度势的求解时，可以极大地降低网格离散量并提高数值的计算效率。格林函数的数值计算已成为船舶耐波性数值计算的关键。

三维脉动源格林函数被用于规则波中船舶运动响应的求解。1942年，Havelock^[2]推导出了这一格林函数的解析积分表达式。对三维脉动源格林函数进行数值计算，通常会采用数值积分法、级数展开法或多项式逼近法。Newman^[3]采用Romberg积分法获得了10⁻¹⁵精度的频域零航速自由面格林函数的数值结果。尽管如此，受限于计算效率，实际水动力计算中很少使用此数值积分法。Newman^[4]在采用多项式逼近方法计算零航速自由面函数的数值结果时，将整个计算域划分为4个子区域，在不同的子区域将格林函数分解为不同的解析函数和剩余函数，并采用双重Chebyshev多项式对剩余函数进行逼近。这一算法数值效率高，数值精度高达10⁻⁶，已被成功地应用到WAMIT软件中。Chen^[5]采用另一种多项式逼近的方法对Noblesse^[6]提出的剩余函数进行了数值逼近，并将该算法应用到了水动力学软件HydroStar中。此外，Clément^[7]发现了三维脉动源格林函数满足的四阶常微分方程。Wu等^[8]提出了一种全域内计算格林函数的多项式逼近公式，避免了计算域的分区，提高了数值计算的效率。

本文对Wu等^[8]的格林函数全域计算公式进行了改进，将格林函数与该全域计算公式的差值作为剩余函数，采用双重Chebyshev多项式对剩余函数进行逼近，从而提高了全域计算公式的精度。将本文格林函数计算方法应用到浮体波浪中运动响应及二阶平均漂移力的预报中，并与试验结果及其他数值方法进行比较，数值结果验证了本文计算方法在水动力计算中的可靠性。

1 理论基础 1.1 频域零航速自由面格林函数

建立描述浮体运动的笛卡尔坐标系。z轴竖直向上，o-xy平面与自由面重合。对于频域稳态问题，若规则入射波的圆频率为ω，浮体的简谐摇荡运动以及绕射力的频率均为ω。无限水深频域零航速自由面格林函数G(P,Q)e^−iωt表示源点$ Q\left(\tilde{x},\; \tilde{y},\; \tilde{z}\right) $在点P(x, y, z)处生成的速度势，G(P,Q)满足如下控制方程和边界条件^[8]：

$ \left\{ \begin{gathered} {\nabla ^2}G\left( {P,Q} \right) = - 4{\text{π}} \delta \left( {P - Q} \right)，{\text{ }}z < 0，\\ \frac{{\partial G}}{{\partial z}} - \nu G = 0，{\text{ }}z = 0，\\ \mathop {\lim }\limits_{z \to - \infty } \nabla G = 0 ，\\ \mathop {\lim }\limits_{R \to \infty } \sqrt R \left( {\frac{{\partial G}}{{\partial R}} - i\nu G} \right) = 0 。\\ \end{gathered} \right. $

(1)

式中：δ为狄拉克函数；k₀ = ω²/g为波数；g为重力加速度；R为场点和源点的水平距离。由式(1)可推导出频域零航速自由面格林函数的表达式如下：

$\begin{split} &G\left( {P,Q} \right) = \frac{1}{r} + \frac{1}{d} + \\ & \frac{{{k_0}}}{{\text{π}} }\int_0^\infty {{\mathrm{d}}k\int_0^{2{\text{π}} } {\dfrac{{{e^{k\left( {z + \tilde z} \right) + ik\left[ {\left( {x - \tilde x} \right)\cos \theta + \left( {y - \tilde y} \right)\sin \theta } \right]}}}}{{k - {k_0}}}{\rm{d}}\theta } } 。\end{split}$

(2)

式中：r和d分别为场点和源点及源点在水平面上镜像的距离。即：

$ \left\{ \begin{gathered} r = \sqrt {{{\left( {x - \tilde x} \right)}^2} + {{\left( {y - \tilde y} \right)}^2} + {{\left( {z - \tilde z} \right)}^2}} ，\\ d = \sqrt {{{\left( {x - \tilde x} \right)}^2} + {{\left( {y - \tilde y} \right)}^2} + {{\left( {z + \tilde z} \right)}^2}} 。\\ \end{gathered} \right. $

(3)

Newman^[3]经过一系列的推导，将频域零航速自由面格林函数简化为如下积分形式：

$ G\left( {P,Q} \right) = \frac{1}{r} + \frac{1}{d} + {k_0}F\left( {X,Y} \right) + 2{\text{π}} i{k_0}{e^{ - Y}}{J_0}\left( X \right) 。$

(4)

其中，

$ \begin{split} F\left( {X,Y} \right) = \;& - 2{e^{ - Y}}\int_0^Y {{{\left( {{X^2} + {t^2}} \right)}^{{{ - 1} \mathord{\left/ {\vphantom {{ - 1} 2}} \right. } 2}}}{e^t}{\rm{d}}t} - \\ \;&{\text{π}} {e^{ - Y}}\left[ {{Y_0}\left( X \right) + {H_0}\left( X \right)} \right] 。\end{split}$

(5)

J_n(X)为第一类n阶贝塞尔函数（Bessel function）；Y_n(X)为第二类n阶贝塞尔函数；H_n(X)为n阶斯特鲁夫函数（Struve function）；X=k₀R；$ {Y}=-{k}_{0}\left(z+\tilde {z}\right) $。函数F(X, Y)的导数可通过如下公式计算：

$ \begin{split} {F_X}\left( {X,Y} \right) = \;& 2X\int_0^Y {{{\left( {{X^2} + {t^2}} \right)}^{{{ - 3} \mathord{\left/ {\vphantom {{ - 3} 2}} \right. } 2}}}{e^{t - Y}}{\rm{d}}t} + \\ \;&{\text{π}} {e^{ - Y}}\left[ {{Y_1}\left( X \right) + {H_1}\left( X \right)} \right] - 2{e^{ - Y}} ，\end{split}$

(6)

$ {F_Y}\left( {X,Y} \right) = {{ - 2} \mathord{\left/ {\vphantom {{ - 2} {\sqrt {{X^2} + {Y^2}} }}} \right. } {\sqrt {{X^2} + {Y^2}} }} - F\left( {X,Y} \right) 。$

(7)

1.2 格林函数的全域计算公式

Noblesse^[6]将格林函数表达式为：

$ G = - \frac{1}{r} + L + W。$

(8)

式中，波动项W表示为：

$ W\left( {h,v} \right) = 2{\text{π}} \left[ {{{\tilde H}_0}\left( h \right) - i{J_0}\left( h \right)} \right]{e^v}。$

(9)

式中，$ {\tilde H_n}\left( h \right) $为第n类零阶Struve函数。与之对应的，近场项L为：

$ L\left( {h,v} \right) = - \dfrac{1}{d} - \dfrac{4}{{\text{π}} }\int_0^{\frac{{\text{π}} }{2}} {{Re} \left\{ {{e^M}{E_1}\left( M \right)} \right\}{\rm{d}}\theta }。$

(10)

式中：$ h = \sqrt {{{\left( {x - \tilde x} \right)}^2} + {{\left( {y - \tilde y} \right)}^2}} $ ，$ v = z + \tilde z $，M = v + ihcosθ；E₁(·)为复指数积分函数。

格林函数的梯度▽G ≡ (G_x, G_y, G_z)被表示为：

$ \left\{\begin{aligned}&{G}_{z}= \dfrac{z-\tilde{z}}{{r}^{3}}+{L}_{z}+W，{L}_{z}=\dfrac{v}{{d}^{3}}-\frac{1}{d}+L ，\\ &{G}_{h}=\dfrac{h}{{r}^{3}}+{L}_{h}+{W}_{h}，{L}_{h}=\dfrac{h}{{d}^{3}}+{L}_{*}，\\ &{G}_{x}={G}_{h}\dfrac{x-\tilde{x}}{h}，\\ &{G}_{y}={G}_{h}\dfrac{y-\tilde{y}}{h}。&\end{aligned}\right. $

(11)

式(11)中波动项W_h表达式为：

$ {W_h}\left( {h,v} \right) = 2{\text{π}} \left[ {{2 \mathord{\left/ {\vphantom {2 {\text{π}} }} \right. } {\text{π}} } - {{\tilde H}_1}\left( h \right) + i{J_1}\left( h \right)} \right]{e^v}。$

(12)

近场项$ {L}_{*} $表达式为：

$ {L_*}\left( {h,v} \right) = - \frac{4}{{\text{π}} }\int_0^{\frac{{\text{π}} }{2}} {{Im} \left\{ {{e^M}{E_1}\left( M \right) - {1 \mathord{\left/ {\vphantom {1 M}} \right. } M}} \right\}\cos \theta {\rm{d}}\theta }。$

(13)

Wu等^[8]给出了L的实用近似计算公式：

$ L\left( {h,v} \right) \approx - \frac{1}{d} + \frac{{2P}}{{1 + {d^3}}} + 2\rho {\left( {1 - \rho } \right)^3}R。$

(14)

式中：

$ P = {e^v}\left( {\ln \frac{{d - v}}{2} + \gamma - 2{d^2}} \right) + {d^2} - v，$

(15)

$ R = \left( {1 - \beta } \right)A + \beta B - \frac{{\alpha C}}{{1 + 6\alpha \rho \left( {1 - \rho } \right)}} + \beta \left( {1 - \beta } \right)D。$

(16)

式中：$ \gamma =0.577 $为欧拉常数，$ \alpha =-v/d $，$ \mathrm{\beta }={h}/{d} $和$ \rho =d/\left(1+d\right) $。式(16)中关于$ \rho $的多项式$ A\left(\rho \right) $、$ B\left(\rho \right) $、C$ \left(\rho \right) $和D$ \left(\rho \right) $的定义如下：

$ \left\{\begin{split} &\begin{split} {A}=& 1.21-13.328\rho +215.896{\rho }^{2}-1763.96{\rho }^{3}+\\ &8418.94{\rho }^{4}-24314.21{\rho }^{5}+42002.57{\rho}^{6}-\\ &41592.9{\rho}^{7}+21859{\rho }^{8}-4838.6{\rho }^{9} ，\end{split}\\ &\begin{split} {B}=& 0.938-5.373\rho -67.92{\rho }^{2}+796.534{\rho }^{3}-\\ & 4780.77{\rho }^{4}+17137.74{\rho }^{5}-36618.81{\rho }^{6}+\\ & 44894.06{\rho }^{7}- 29030.24{\rho }^{8}+7671.22{\rho }^{9} ，\end{split}\\ &\begin{split} {C}=&1.268-9.747\rho +209.653{\rho }^{2}-1397.89{\rho }^{3}+\\ & 5155.67{\rho }^{4}-9844.35{\rho }^{5}+\\ & 9136.4{\rho }^{6}-3272.62{\rho }^{7} ，\end{split}\\ &\begin{split} {D}= &0.632-40.97\rho +667.16{\rho }^{2}-6072.07{\rho }^{3}+\\ & 31127.39{\rho }^{4}-96293.05{\rho }^{5}+181856.75{\rho }^{6}-\\ & 205690.43{\rho }^{7}+128170.2{\rho }^{8}-33744.6{\rho }^{9}。\end{split}\\ \end{split} \right.$

(17)

L*的实用近似计算公式为^[7]：

$ {L_ * }\left( {h,v} \right) \approx \frac{{2{P_ * }}}{{1 + {d^3}}} - 4{Q_ * } + 2\rho {\left( {1 - \rho } \right)^3}{R_ * } 。$

(18)

其中，$ {P}_{*} $、$ {Q}_{*} $和$ {R}_{*} $的表达式为：

$ {P_*} = \frac{{\beta + h}}{{d - v}} - 2\beta + {e^v}d - h ，$

(19)

$ {Q_*} = {e^{ - d}}\left( {1 - \beta } \right)\left( {1 + \frac{d}{{1 + {d^3}}}} \right) ，$

(20)

$ {R_*} = \beta {A_*} - \left( {1 - \alpha } \right){B_*} + \beta \left( {1 - \beta } \right)\rho \left( {1 - 2\rho } \right){C_*}。$

(21)

式(21)中关于$ \rho $的多项式$ {A}_{\mathrm{*}}\left(\rho \right) $、$ {B}_{\mathrm{*}}\left(\rho \right)\mathrm{和}{C}_{\mathrm{*}}\left(\rho \right) $的定义如下：

$\left\{\begin{split} & \begin{split} {A}_{*}=& 2.948-24.53\rho +249.69{\rho }^{2}-754.85{\rho }^{3}-\\ & 1187.71{\rho }^{4}+ 16370.75{\rho }^{5}-48811.41{\rho }^{6}+\\ & 68220.87{\rho }^{7}- 46688{\rho }^{8}+12622.25{\rho }^{9} ，\end{split}\\& \begin{split} {B}_{*}=& 1.11+2.894\rho -76.765{\rho }^{2}+1565.35{\rho }^{3}-\\ & 11336.19{\rho }^{4}+44270.15{\rho }^{5}-97014.11{\rho }^{6}+\\ & 118879.26{\rho }^{7}-76209.82{\rho }^{8}+19923.28{\rho }^{9} ，\end{split}\\& \begin{split} {C}_{*}=& 14.19-148.24\rho +847.8{\rho }^{2}-2318.58{\rho }^{3}+\\ & 3168.35{\rho }^{4}-1590.27{\rho }^{5} 。\end{split}\\ \end{split}\right. $

(22)

1.3 剩余函数的选取及多项式的逼近

Liang等^[9]指出，式(14)的误差大于10⁻³。为提高格林函数计算精度，将式(8)及式(14)中第3项的一半选为误差函数并采用Chebyshev多项式进行数值逼近，即：

$ \begin{split} 2f =& G(P,Q) + \frac{1}{r} + \frac{1}{d} - \frac{{2P}}{{1 + {d^3}}} - 2\rho {\left( {1 - \rho } \right)^3}R - \\ & 2{\text{π}} {e^{ - Y}}\left[ {{H_0}\left( X \right) - {Y_0}\left( X \right)} \right] 。\end{split}$

(23)

则有剩余函数的表达式为：

$ \begin{aligned} f = \;& - \frac{1}{2}F\left( {X,Y} \right) - \frac{P}{{1 + {d^3}}} - {\text{π}} {H_0}\left( X \right){e^{ - Y}} = \\ \;& {e^{ - Y}}\int_0^Y {{{\left( {{X^2} + {t^2}} \right)}^{{{ - 1} \mathord{\left/ {\vphantom {{ - 1} 2}} \right. } 2}}}{e^t}{\rm{d}}t} - \frac{P}{{1 + {d^3}}} + \\ \;& \frac{{\text{π}} }{2}{e^{ - Y}}\left[ {{Y_0}\left( X \right) - {H_0}\left( X \right)} \right] 。\end{aligned}$

(24)

在使用Chebyshev多项式进行多项式逼近之前，先采用龙贝格积分^[2]获得式(23)中单重积分的值。剩余函数的双重Chebyshev多项式表达式为：

$ f\left( {\xi ,\eta } \right) \approx \sum\limits_{m = 0}^M {\sum\limits_{n = 0}^N {{d_{mn}}{T_m}\left( \xi \right){T_n}\left( \eta \right)} }。$

(25)

Chebyshev多项式系数表达式为：

$ {d_{mn}} = \frac{{{\varepsilon _m}{\varepsilon _n}}}{{MN}}\sum\limits_{m = 0}^{M - 1} {\sum\limits_{n = 0}^{N - 1} {f\left( {{\xi _i},{\eta _j}} \right){T_m}\left( {{\xi _i}} \right){T_n}\left( {{\eta _j}} \right)} }。$

(26)

式中：$ {\xi }_{i}=\mathrm{cos}\dfrac{2i+1}{2M}{\text{π}} $；$ {\eta }_{j}=\mathrm{cos}\dfrac{2j+1}{2N}{\text{π}} $，$ {T}_{m}\left({\xi }_{i}\right)=\mathrm{cos} \dfrac{2i+1}{2M}m{\text{π}} $；$ {T}_{n}\left({\eta }_{j}\right)=\mathrm{cos}\dfrac{2j+1}{2N}n{\text{π}} $；$ {\xi}_{0}=1 $；$ {\xi}_{j}=2\left(j\geqslant 1\right) $。对于在无限流域上的剩余函数，对X、Y作如下变换：

$ \left\{ \begin{gathered} \xi = \frac{{X - 1}}{{X + 1}} ，\\ \eta = \frac{{Y - 1}}{{Y + 1}} 。\\ \end{gathered} \right.$

(27)

此时半无穷流域：$ 0\leqslant {X} < \mathrm{\infty }\text{、}0\leqslant {Y} < \mathrm{\infty } $映射到$ -1\leqslant \xi \leqslant 1 $、$ -1\leqslant \eta \leqslant 1 $的单位空间。

1.4 流体速度势及水动力计算

在线性频率势流理论框架内，可以采用速度势表示流体的速度，速度势满足拉普拉斯方程，可被表示为：

$ \Phi \left( {x, y, z; t} \right) = \left[ {{\phi _I}\left( {x,y,z} \right) + {\phi _D}\left( {x,y,z} \right) + {\phi _R}\left( {x,y,z} \right)} \right]{e^{ - i\omega t}} 。$

(28)

式中，$ {\phi _I} $、$ {\phi _D} $和$ {\phi _R} $分别为入射势、辐射势和绕射势，ω为波浪圆频率。

进一步地，辐射速度势$ {\phi _R} $可被分解为：

$ {\phi _R} = i\omega \sum\limits_{j = 1}^6 {{\xi _j}{\phi _j}}。$

(29)

式中，${\xi _j}$为第j个运动模态的幅值；${\phi _j}$为规范化速度势，控制方程和边界条件如下：

$ \begin{cases} \nabla^2 \phi_j = 0 ，{{流体域内}}，\\ \dfrac{\partial \phi_j}{\partial z} - v \phi_j = 0 ，{{平均水线面上}}，\\ \dfrac{\partial \phi_j}{\partial n} = n_j ，{{船体湿表面上}}，\\ \lim_{z \to -\infty} \phi_j = 0 ，{{底部条件}}，\\ \lim_{R \to \infty} \sqrt{R} \left( \dfrac{\partial \phi_j}{\partial R} - iv \phi_j \right) = 0 ，{{远方辐射条件}}。\end{cases}$

(30)

式中，${n_j}$为浮体湿表面第j个运动模态的法向量。

可采用源分布模型或者源-偶混合分布模型进行速度势的求解，具体如下：

$ \alpha \left( P \right)\phi \left( P \right) + \iint_{{s_B}} {\phi \left( Q \right)}\frac{{\partial G\left( {P,Q} \right)}}{{\partial {n_Q}}}{\rm{d}}s = \iint_{{s_B}} {G\left( {P,Q} \right)}\frac{{\partial \phi \left( Q \right)}}{{\partial {n_Q}}}{\rm{d}}s，$

(31)

$ \phi \left( P \right) = \iint_{{s_B}} {\sigma \left( Q \right)G\left( {P,Q} \right)}{\rm{d}}s。$

(32)

其中，式(31)为源分布模型；式(32)为源-偶混合分布模型，S_B为流体边界。

浮体附加质量和阻尼系数可通过下式计算获得：

$ {\mu _{jk}} + \frac{i}{\omega }{\lambda _{jk}} = \rho \iint_{{s_B}} {{\phi _k}{n_j}\sigma \left( Q \right)G\left( {P,Q} \right)}{\rm{d}}s。$

(33)

式中：$ {\mu _{jk}} $为附加质量系数；$ {\lambda _{jk}} $为阻尼系数。

在频域范围内，船体在波浪中的运动方程可以被表达为：

$ \left\{ { - {\omega ^2}\left[ {{\boldsymbol{M}} + {\mathbf{\mu }}\left( \omega \right)} \right] - i\omega {\mathbf{\lambda }}\left( \omega \right) + {\boldsymbol{C}}} \right\}{\mathbf{\xi }} = {{\boldsymbol{f}}^{F - K}} + {{\boldsymbol{f}}^D}。$

(34)

M为质量矩阵，有：

$ {{\boldsymbol{M}} = \left[ {\begin{array}{*{20}{c}} m&0&0&0&{m{z_g}}&{ - m{y_g}} \\ 0&m&0&{ - m{z_g}}&0&{m{x_g}} \\ 0&0&m&{m{y_g}}&{ - m{x_g}}&0 \\ 0&{ - m{z_g}}&{m{y_g}}&{{I_{11}}}&{{I_{12}}}&{{I_{13}}} \\ {m{z_g}}&0&{ - m{x_g}}&{{I_{21}}}&{{I_{22}}}&{{I_{23}}} \\ { - m{y_g}}&{m{x_g}}&0&{{I_{31}}}&{{I_{32}}}&{{I_{33}}} \end{array}} \right]。} $

(35)

式中：(x_g,y_g,z_g)为船体在动坐标系下的重心坐标；I_ij为浮体转动惯量。μ(ω)为波浪圆频率ω对应的附加质量矩阵，有

$ {\boldsymbol{\mu }}\left( \omega \right) = \left[ {\begin{array}{*{20}{c}} {{\mu _{11}}}&{{\mu _{12}}}&{{\mu _{13}}}&{{\mu _{14}}}&{{\mu _{15}}}&{{\mu _{16}}} \\ {{\mu _{21}}}&{{\mu _{22}}}&{{\mu _{23}}}&{{\mu _{24}}}&{{\mu _{25}}}&{{\mu _{26}}} \\ {{\mu _{31}}}&{{\mu _{32}}}&{{\mu _{33}}}&{{\mu _{34}}}&{{\mu _{35}}}&{{\mu _{36}}} \\ {{\mu _{41}}}&{{\mu _{42}}}&{{\mu _{43}}}&{{\mu _{44}}}&{{\mu _{45}}}&{{\mu _{46}}} \\ {{\mu _{51}}}&{{\mu _{52}}}&{{\mu _{53}}}&{{\mu _{54}}}&{{\mu _{55}}}&{{\mu _{56}}} \\ {{\mu _{61}}}&{{\mu _{62}}}&{{\mu _{63}}}&{{\mu _{64}}}&{{\mu _{65}}}&{{\mu _{66}}} \end{array}} \right]。$

(36)

λ(ω)为波浪圆频率ω对应的阻尼系数矩阵，有

$ {\boldsymbol{\lambda }}\left( \omega \right) = \left[ {\begin{array}{*{20}{c}} {{\lambda _{11}}}&{{\lambda _{12}}}&{{\lambda _{13}}}&{{\lambda _{14}}}&{{\lambda _{15}}}&{{\lambda _{16}}} \\ {{\lambda _{21}}}&{{\lambda _{22}}}&{{\lambda _{23}}}&{{\lambda _{24}}}&{{\lambda _{25}}}&{{\lambda _{26}}} \\ {{\lambda _{31}}}&{{\lambda _{32}}}&{{\lambda _{33}}}&{{\lambda _{34}}}&{{\lambda _{35}}}&{{\lambda _{36}}} \\ {{\lambda _{41}}}&{{\lambda _{42}}}&{{\lambda _{43}}}&{{\lambda _{44}}}&{{\lambda _{45}}}&{{\lambda _{46}}} \\ {{\lambda _{51}}}&{{\lambda _{52}}}&{{\lambda _{53}}}&{{\lambda _{54}}}&{{\lambda _{55}}}&{{\lambda _{56}}} \\ {{\lambda _{61}}}&{{\lambda _{62}}}&{{\lambda _{63}}}&{{\lambda _{64}}}&{{\lambda _{65}}}&{{\lambda _{66}}} \end{array}} \right]。$

(37)

C为回复力矩阵，有

$ {\begin{split}&{\boldsymbol{C}} = \\ & \left[ \begin{array}{*{20}{c}} 0 &0& 0& 0& 0& 0 \\ 0 &0& 0& 0& 0& 0 \\ 0 &0& {\rho g{A_{wp}}}& 0& { - \rho g\iint\limits_{wp} {x{\rm{d}}x{\rm{d}}y}}& 0 \\ 0 &0& 0& \begin{split}\rho & g[ V( {{z_b} - {z_g}} ) + \\ & \iint\limits_{wp} {{y^2}{\rm{d}}x{\rm{d}}y} ]\end{split}& 0& 0 \\ 0 &0& { - \rho g\iint\limits_{wp} {x{\rm{d}}x{\rm{d}}y}}& 0& \begin{split} & \rho g[V( {{z_b} - {z_g}} ) + \\ & \iint\limits_{wp} {{x^2}{\rm{d}}x{\rm{d}}y} ]\end{split}& 0 \\ 0 &0& 0 & 0 & 0 & 0 \end{array} \right] 。\end{split}}$

(38)

式中：(x_b,y_b,z_b)为船体在动坐标系下的浮心坐标；V为船舶排水体积；下标wp的双重积分表示沿船舶水线面的积分；A_wp为船体水线面面积。

计算二阶平均漂移力的远场积分法也叫远场法，Maruo^[10]和Newman^[11]分别根据动量定理和能量守恒定理导出了计算零航速浮体纵荡和横荡方向上二阶平均漂移力及艏摇力矩的远场公式。

$ \bar F_1^{\left( 2 \right)} = \frac{{\rho {\nu ^2}}}{{8{\text{π}} }}\int_0^{2{\text{π}} } {{{\left| {H\left( \theta \right)} \right|}^2}\left( {\cos \beta - \cos \theta } \right){\rm{d}}\theta }，$

(39)

$ \bar F_2^{\left( 2 \right)} = \frac{{\rho {\nu ^2}}}{{8{\text{π}} }}\int_0^{2{\text{π}} } {{{\left| {H\left( \theta \right)} \right|}^2}\left( {sin\beta - sin\theta } \right){\rm{d}}\theta }，$

(40)

$ \bar F_6^{\left( 2 \right)} = - \frac{{\rho \nu }}{{8{\text{π}} }}\int_0^{2{\text{π}} } {{H^ * }\left( \theta \right)H'\left( \theta \right){\rm{d}}\theta } + \frac{1}{{2\nu }}\rho \omega {\zeta _a}{Im} \left\{ {H'\left( \beta \right)} \right\} 。$

(41)

式中：ν为波数；ρ为流体密度；β为浪向角；H(θ) 为柯钦函数，计算式为

$ H\left( \theta \right) = \iint_{{S_b}} {\left( {\frac{{\partial {\phi _B}}}{{\partial n}} - {\phi _B}\frac{\partial }{{\partial n}}} \right)}\exp \left[ {\nu \zeta - i\nu \left( {\xi \cos \theta + \eta \sin \theta } \right)} \right]。$

(42)

式中：${\phi _B}$为扰动速度势，表示绕射势$ {\phi _D} $和辐射势$ {\phi _R} $的和，H*(θ)为H(θ)的共轭复数；H'(θ)为H(θ)关于θ的导数；ζ_a为波浪幅值。

压力积分法也叫近场法，Pinkster^[11]通过直接在船体湿表面上进行压力积分获得船体受到的二阶漂移力：

$ \begin{split} {{\boldsymbol{F}}^{\left( 2 \right)}} =& \int_{WL} {\frac{{\rho g}}{2}\left[ {\zeta _r^{\left( 1 \right)} \cdot \zeta _r^{\left( 1 \right)}} \right]{\boldsymbol{n}}dl} - \frac{\rho }{2}\iint_{{S_b}} \left[ {\nabla {\Phi ^{\left( 1 \right)}} \cdot \nabla {\Phi ^{\left( 1 \right)}}} \right]\\ &{\boldsymbol{n}}dS - \rho \iint_{{S_b}} \left[ {\left( {{{\boldsymbol{X}}_T} + {{\boldsymbol{X}}_R} \times {\boldsymbol{r}}} \right) \cdot \nabla \frac{{\partial {\Phi ^{\left( 1 \right)}}}}{{\partial t}}} \right]\\ &{\boldsymbol{n}}dS - {{\boldsymbol{X}}_R} \times {{\boldsymbol{F}}^{\left( 1 \right)}} - \rho \iint_{{S_b}} {\frac{{\partial {\Phi ^{\left( 2 \right)}}}}{{\partial t}}{\boldsymbol{n}}dS } ，\end{split} $

(43)

$ \begin{split} &{{\boldsymbol{M}}^{\left( 2 \right)}} = \int_{WL} {\frac{{\rho g}}{2}\left[ {\zeta _r^{\left( 1 \right)} \cdot \zeta _r^{\left( 1 \right)}} \right]\left( {{\boldsymbol{r}} \times {\boldsymbol{n}}} \right)dl} -\frac{\rho }{2}\iint_{{S_b}} \left[ \nabla {\Phi ^{\left( 1 \right)}} \cdot \right.\\ &\left. \nabla{\Phi ^{\left( 1 \right)}} \right] \left( {{\boldsymbol{r}} \times {\boldsymbol{n}}} \right)dS - \rho \iint_{{S_b}} \left[ \left( {{{\boldsymbol{X}}_T} + {{\boldsymbol{X}}_R} \times {\boldsymbol{r}}} \right) \cdot \nabla \frac{{\partial {\Phi ^{\left( 1 \right)}}}}{{\partial t}} \right] \\ & \left( {\boldsymbol{r}} \times {\boldsymbol{n}} \right)dS - {{\boldsymbol{X}}_R} \times {{\boldsymbol{M}}^{\left( 1 \right)}} - \rho \iint_{{S_b}} \frac{{\partial {\Phi ^{\left( 2 \right)}}}}{{\partial t}}\left( {{\boldsymbol{r}} \times {\boldsymbol{n}}} \right)dS 。\end{split} $

(44)

式中：F⁽¹⁾和M⁽¹⁾表示船体受到的一阶波浪力和力矩; ζ_r⁽¹⁾为水线处波面相对船体的升高；X_T和X_R分别为浮体的线位移和角位移。对上述式(43)和式(44)取平均值即可获得二阶平均漂移力。

2 数值结果与分析

采用双重Chebyshev多项式系数对上述剩余函数进行Chebyshev多项式逼近的数值误差如图1所示。图中M和N分别为双重Chebyshev多项式的项数。可以看出，M=8、N=10时，剩余函数的误差在10⁻³量级；M=50、N=20时，绝大多数剩余函数的误差在10⁻⁶量级，少数点位误差在10⁻⁵量级。随着M和N的增大，剩余函数的误差逐渐变小。理论上M和N取足够大的数量，剩余函数的误差可以无限接近于0。考虑到格林函数的计算效率，本文截取Chebyshev多项式M=8、N=10，进行剩余函数的拟合。

3 格林函数在船舶水动力计算中的应用 3.1 数值算例

选用半球浮体、Wigley-3数学船型^[12]和驳船^[13]作为算例，采用本文格林函数计算方法计算这些浮体在特定海况下的运动响应以及二阶平均漂移力，验证本文频域零航速自由面格林函数计算方法在浮体水动力计算中的有效性。与Wigley-3 和驳船的船型相关参数如表1所示。

3.2 波浪中的运动响应预报

本节计算了Wigley-3在迎浪海况下以及驳船在迎浪以及浪向角135°海况下的运动RAO。本文格林函数计算方法以及Romberg积分法计算格林函数均被用于Wigley-3和驳船在特定海况下的RAO的数值计算中。除此以外，商业软件AQWA也被用于驳船的RAO计算中。Wigley-3和驳船平均湿表面网格划分见图2和图3。Wigley-3迎浪海况下在X轴方向上运动RAO的数值结果与试验值的对比分别见图4。驳船在135°浪向以及迎浪海况下运动RAO的数值结果与试验值的对比分别见图5和图6。这些图中，“本文方法”为本文所用格林函数计算方法计算所得数值结果；Romberg为Romberg积分计算格林函数获得的数值结果；Wigley-3的试验结果源自Journee的试验报告^[12]。图5和图6驳船的RAO计算结果为商业软件AQWA计算所得数值结果，试验数据源自文献[13]。

从图4～图6可以看出，本文格林函数计算方法计算所得Wigley-3在迎浪工况下的纵向波浪力、垂荡和纵摇运动RAO以及驳船在135°浪向和迎浪下的运动RAO结果与Romberg积分获得的高精度格林函数计算所得水动力结果完全重合。这些数值结果与试验结果以及商业软件AQWA计算所得数值结果并没有显著的区别，这说明本文零航速自由面格林函数的计算方法可以被用于波浪中浮体的运动响应预报。

3.3 二阶平均漂移力预报

分别采用本文格林函数计算方法和高精度格林函数（Romberg）计算了半球浮体在完全固定和自由漂浮情况下在波浪中受到的纵向二阶平均漂移力，并与试验值和其他学者的数值结果进行对比，详见图7。图中二阶力的计算采用的远场公式^{[10 – 11]}；“Analytical”和“Exp”分别为解析解和试验值，均源自文献[14]；“HOBEM”为Choi等^[15]采用高阶面元法计算所得数值结果。

与半球浮体的二阶力计算类似，采用本文格林函数计算方法和高精度的格林函数计算了驳船在135°浪向下受到的纵向、横向以及艏摇方向的二阶平均漂移力和力矩，并与试验结果进行了对比，详见图8。图中“本文方法”和“Romberg”分别表示采用本文格林函数计算方法以及高精度格林函数计算所得二阶力数值结果；“Cal”和“Exp”源自文献[13]的数值结果以及试验数据。鉴于文献[13]数值结果是近场公式计算所得的，而本文自编的程序在计算二阶力时采用的远场公式；计算驳船在这2个浪向下受到的二阶平均漂移力时，增加了在商业软件AQWA中分别使用远场公式和近场公式计算所得数值结果，即图中“Far field”和“Near field”。此外，本文还计算了该驳船在迎浪状态下受到的纵向二阶平均漂移力，具体结果详见图9。

由图7 ～图9可知，采用本文格林函数计算方法计算所得二阶平均漂移力与高精度格林函数计算所得数值结果完全重合。这说明了本文格林函数计算方法可以被用于波浪中浮体受到的二阶平均漂移力的预报。图中商业软件AQWA中使用远场公式计算所得二阶力数值结果与本文自编程序计算所得数值结果存在较小的差异，这可能是由于AQWA和自编程序求解边界积分方程的方式不同引起的。AQWA采用的是分布源模型，而本文采用的是源-偶混合分布模型。

4 结　语

本文提出了一种新的三维脉动源格林函数的数值计算方法，并在浮体运动响应和二阶平均漂移力的计算中进行验证，得出如下结论：

1） 本文方法可以提高格林函数全域计算公式的数值精度，随着Chebyshev多项式项数增加，剩余函数拟合误差理论上可以无限接近于0。

2） 本文格林函数计算方法在浮体运动响应和二阶平均漂移力计算中是可靠且可行的。