波浪与浮体相互作用的时域模拟具有重要意义^[1]。在势流理论框架下，根据所选取格林函数的 不同，主要分为时域Rankine源法和时域自由面格 林函数法。对于时域自由面格林函数法，只需要在 船体表面进行积分，国内外许多学者都在此方面做了研究^[2-3]。此方法对于不同浮体进行模拟时，都 需要重新计算时域格林函数在每一时刻的值，而时域格林函数的计算又需要占用较长的CPU时间，所 以该方法计算效率较低。此外，格林函数在近水面单元的高频振荡特性^[4]，导致外飘型船舶的计算，结果时常发散。对于时域Rankine源法，其形式简单，且易与高阶边界元法相结合，以进行非线性问题 的模拟，所以现在仍然有很多学者对此进行研究^[5-6]。但是，时域Rankine法需要对计算域进行截断，且需要布置相应的人工远方辐射条件。

通过对时域Rankine源法和时域自由面格林函 数法各自优缺点的分析，有学者提出了匹配方 法^[7-8]，此种方法不仅避免了不同浮体都需对格林 函数进行计算的问题，同时适用于外飘型船舶的计 算。但是此方法的准确性和有效性仍需进一步验 证。无论是时域Rankine源法还是匹配方法，在计 算与时间相关的记忆速度势时，时间卷积项的存在，常常导致长时间的模拟无法进行。

针对上述问题，本文首先建立了三维内外场时 域快速匹配边界元模型，通过与解析解的对比，验证 了数值模型的有效性。其次利用该模型研究了流场 中不同点的记忆速度势随时间的变化特点。最后根 据记忆速度势随时间的变化特性，提出了利用 Laguerre函数^[9]对记忆速度势时间变量进行级数逼 近，并通过引人一放缩系数加快了级数的收敛速度。

1 基础理论 1.1 内外场匹配模型

根据文献[10]，假定流场为不可压缩的理想流体，则流场中的物理量可以用流体速度势φ来进行 描述。为建立三维时域匹配模型，引人一虚拟控制面，如图 1所示。

流体速度势在内域满足以Rankine源所建立的边界积分方程：

$ 2\pi {\phi ^I} + \int {\int\limits_S {{\phi ^I}} \frac{{\partial {G_r}}}{{\partial {n_q}}}d{S_q}} = \int {\int\limits_S {{G_r}\frac{{\partial {\phi ^I}}}{{\partial {n_q}}}d{S_q}} } $

(1)

式中：G_r=1/|p-q|为简单格林函数。

流体速度势在外域满足以时域自由面格林函数 所建立的边界积分方程：

$ - 2\pi {\phi ^{II}} + \int {\int\limits_c {[{\phi ^{II}}\frac{{\partial {G_0}}}{{\partial {n_q}}} - {G_0}\frac{{\partial {\phi ^{II}}}}{{\partial {n_q}}}]d{S_q}} } = \int\limits_0^t {dT\int {\int\limits_c {[\tilde G\frac{{\partial {\phi ^{11}}}}{{\partial {n_q}}} - {\phi ^{II}}\frac{{\partial \tilde G}}{{\partial {n_q}}}]d{S_q}} } } $

(2)

式中：G₀和$ \tilde G $分别为时域自由面格林函数的瞬时项和波动项。根据文献[10]，时域自由面格林函数可以写为

$ \left\{ \begin{array}{l} G\left( {p,q,t} \right) = \delta \left( t \right) \times {G_0} + H\left( t \right) \times \tilde G\\ \\ {G_0}\left( {p,q} \right) = \left( {\frac{1}{r} - \frac{1}{{r'}}} \right)\\ \\ \bar G\left( {p,q,t} \right) = 2\int\limits_0^\infty {\sqrt {gk} } \times {e^{k\left( {z + \zeta } \right)}}{J_0}\left( {kR} \right)\sin \left( {\sqrt {gkt} } \right)dk \end{array} \right. $

(3)

式中：r表示流场中场点p(x，y，z)和源点q(ξ，η，ζ)之间的距离，r'表示场点p(x，y，z)和源 点关于自由面的镜像点q(ξ，η，-ζ)之间的距离。

速度势φ¹(p,t)在内域满足的边界积分方程(１)可以离散为如下矩阵形式：

$ \left[\begin{array}{l} {A_{11}}\;\;\;{A_{12}}\;\;\;{A_{13}}\\ {A_{21}}\;\;\;{A_{22}}\;\;\;{A_{23}}\\ {A_{31}}\;\;\;{A_{32}}\;\;\;{A_{33}} \end{array} \right]\left[\begin{array}{l} \phi _C^I\\ \phi _F^I\\ \phi _B^I \end{array} \right] = \left[\begin{array}{l} {B_{11}}\;\;\;{B_{12}}\;\;\;{B_{13}}\\ {B_{21}}\;\;\;{B_{22}}\;\;\;{B_{23}}\\ {B_{31}}\;\;\;{B_{32}}\;\;\;{B_{33}} \end{array} \right]\left[\begin{array}{l} \tilde \phi _C^I\\ \tilde \phi _F^I\\ \tilde \phi _B^I \end{array} \right] $

(4)

式中：$ \tilde \phi $表示速度势对边界面的法向导数，下标1、2、3分别表示边界面S_C、S_F和S_B。

同理，速度势φ^II在外域满足的边界积分方程(2)也可以离散为相应的矩阵形式

$ C\phi _C^{II} = D\tilde \phi _C^{II} - {r_C} $

(5)

式中：r_C表示外域所满足的边界积分方程的时间卷积项。

在虚拟控制面S_C上，流体速度势及其法向导数应满足如下连续条件：

$ \phi _C^I = \phi _C^{II},\tilde \phi _C^I = \tilde \phi _C^{II} $

(6)

通过速度势在内域满足的边界积分方程式(4)和外域满足的边界积分方程式(5)，便可以得到内外域匹配求解方程组：

$ \left[\begin{array}{l} {A_{11}}\;\;\;{A_{12}}\;\;\;{A_{13}}\\ {A_{21}}\;\;\;{A_{22}}\;\;\;{A_{23}}\\ {A_{31}}\;\;\;{A_{32}}\;\;\;{A_{33}} \end{array} \right]\left[\begin{array}{l} {[C]^{ - 1}}\left( {[D]\left\{ {\tilde \phi _C^{II}} \right\} - \left\{ {{r_C}} \right\}} \right)\\ \;\;\;\;\;\;\;\;\;\;\;\;\;\phi _F^I\\ \;\;\;\;\;\;\;\;\;\;\;\;\;\phi _B^I \end{array} \right] = \left[\begin{array}{l} {B_{11}}\;\;\;{B_{12}}\;\;\;{B_{13}}\\ {B_{21}}\;\;\;{B_{22}}\;\;\;{B_{23}}\\ {B_{31}}\;\;\;{B_{32}}\;\;\;{B_{33}} \end{array} \right]\left[\begin{array}{l} \tilde \phi _C^I\\ \tilde \phi _F^I\\ \tilde \phi _B^I \end{array} \right] $

(7)

式中矩阵的系数如下

${A_{ij}} = \left\{ \begin{array}{l} 2\pi ,\;\;\;\;\;\;i = j\\ \iint \limits_{{S_C} + {S_F} + {S_B}}{\frac{\partial }{{\partial {n_j}}}\left( {\frac{1}{{{r_{ij}}}}} \right)ds,\;\;\;i \ne j} \end{array} \right.$

(8)

$ {B_{ij}} = \iint\limits_{{S_C} + {S_F} + {S_B}} {\frac{1}{{{r_{ij}}}}ds} $

(9)

$ {{C}_{ij}}=\left\{ \begin{align} & -2\pi -\iint\limits_{{{S}_{c}}}{\frac{\partial }{\partial {{n}_{j}}}\left( \frac{1}{r{{'}_{ij}}} \right)ds},\ \ \ \ \ i=j \\ & \iint\limits_{{{S}_{C}}}{\frac{\partial }{\partial {{n}_{j}}}}\left( \frac{1}{{{r}_{ij}}}-\frac{1}{r{{'}_{ij}}} \right)ds,\ \ i\ne j \\ \end{align} \right. $

(10)

$ {{D}_{ij}}=\iint\limits_{{{S}_{C}}}{\left( \frac{1}{{{r}_{ij}}}-\frac{1}{r{{'}_{ij}}} \right)}ds $

(11)

由式(7)，在求解流场速度势时，只有虚拟控制面S_C上涉及了时域自由面格林函数波动项的计算。

1.2 流场速度势分解

由于流体速度势所满足的场方程、边界条件和初始条件均是线性化的，故可以对流体速度势进行线性化分解^[11]

$ \phi \left( {p,t} \right) = {\phi _I}\left( {p,t} \right) + {\phi _D}\left( {p,t} \right) + \sum\limits_{k = 1}^6 {{\phi _k}\left( {p,t} \right)} $

(12)

式中：φ_I、φ_D和φ_k分别为入射波速度势、绕射波速度势和辐射波速度势，其中φ_I为已知量，φ_D和φ_k为待求解的未知量。

对于辐射波速度势φ_k，根据其所代表的物理意义，可以将其分解为与时间无关的瞬时速度势Ψ_k和与时间相关的记忆速度势X_k:

$ \phi \left( {p,t} \right) = {\psi _k}\left( p \right){v_k}\left( t \right) + \int\limits_0^t {{X_k}\left( {p,t - \tau } \right){v_k}\left( \tau \right)d\tau } $

(13)

根据瞬时速度势和记忆速度势在内外域满足的边界条件和初始条件^[3]，利用式(7)便可求出相应的值。根据线性化的Bernoulli方程^[10]便可求出对应的辐射波力

$ {F_{jk}}\left( t \right) = - {μ_{jk}}{v_k}\left( t \right) - \int_0^t {{K_{jk}}\left( {t - \tau } \right){v_k}\left( \tau \right)d\tau } $

(14)

式中：μ_jk和K_jk分别为与频率无关的附加质量和时域延迟函数。其可以相应的转化为频域附加质量A_jk和阻尼系数B_jk:

$ {A_{jk}}\left( \omega \right) = {μ_{jk}} - \frac{1}{\omega }\int_0^\infty {{K_{jk}}\left( \tau \right)\sin \left( {\omega \tau } \right)d\tau } $

(15)

$ {B_{jk}}\left( \omega \right) = \int_0^\infty {{K_{jk}}\left( \tau \right)\cos \left( {\omega \tau } \right)d\tau } $

(16)

1.3 时间卷积积分项

根据流体速度势在外域所满足的边界积分方程(2)可知，随着时间的增长，时间卷积余项{r_C}所占用的内存将逐渐增大，最终可能导致计算失败。

对于内外场匹配边界元，虚拟控制面S_C可以设为规则的半球，流体速度势在外域上所满足的边界积分方程便转化为了在确定积分区域上的第二类Fredholm型积分方程。对于任何函数，根据其自变量的变化区间，可以选择相应的正交多项式进行级数逼近。速度势记忆项的时间变量变化区间为（0,∞），故可用Laguerre函数^[9]进行正交级数展开:

$ {X_k}\left( {p,t} \right) = \sum\limits_{n = 0}^\infty {{a_n}} \left( p \right){l_n}\left( t \right) $

(17)

$ {a_n}\left( p \right) = \int_0^\infty {{X_k}} \left( {p,t} \right){l_n}\left( t \right)dt $

(18)

式中：l_n（x）=L_n（x）e^-x/2，为归一化的Laguerre函数。

根据记忆速度势随时间的变化特性，为加快级数的收敛速度，引入一放缩系数a:

$ {X_k}\left( {p,t} \right) = \sum\limits_{n = 0}^\infty {a{'_n}} \left( p \right)L{'_n}\left( t \right) $

(19)

$ a{'_n}\left( p \right) = \int_0^\infty {{X_k}} \left( {p,t} \right)L{'_n}\left( t \right)dt $

(20)

其中，L'_n（t）为修正的Laguerre函数：

$ L{'_n}\left( t \right) = \frac{1}{{\sqrt a }}{e^{ - \frac{t}{{2a}}}}{L_n}\left( {\frac{t}{a}} \right),0 < a \le 1 $

(21)

2 数值结果 2.1 匹配边界元有效性验证

利用时域快速匹配边界元法求解浮体和波浪相互作用问题时，影响计算精度的主要因素有时间步长、边界网格数目、自由面网格形式、虚拟控制面和物体之间的距离。

为讨论自由面网格形式的影响，针对相同工况，分别选择了矩形网格、辐射状网格和矩形三角形混合网格进行了计算，通过与解析解^[12]的对比验证发现辐射状的自由面网格能够最好的描述辐射波和绕射波的传播，如图 2所示，图中半球半径为1 m，控制面半球半径为5 m。此与辐射和绕射波代表外传柱面波的物理意义相一致。

为讨论时间步长对计算结果的影响，针对相同工况，分别选取了△t=0.02、0.01、0.005 s进行了计算。计算得到的无量纲时域延迟函数如图 3所示，3条曲线几乎重合。此外，针对边界网格数目和虚拟控制面与物体之间的距离也做了同样的对比分析。最后，考虑数值计算精度和计算效率，在文章后续计算中，自由面划分为辐射式网格，计算时间步长选为0.01 s，虚拟控制面和物面间距为5倍的物体特征长度，网格数目为物面130、自由面200、控制面150。

对于半球在静水中做强迫垂荡运动，分别计算了其附加质量系数和阻尼系数，如图 4、5所示。由图可知，利用本文建立的快速时域匹配边界元法计算的结果与文献[12]提供的解析解符合良好，说明本文提供的数值模型能够准确模拟波浪和浮体相互作用的流场。

2.2 Laguerre函数有效性验证

为验证Laguerre函数的有效性，首先对记忆速度势随时间的变化特性进行了研究。记忆速度势代表浮体干扰形成表面波浪，开始波浪传播之后的那部分速度势，并且随时间逐渐衰减，如图 6所示，图中点1位于虚拟控制面和自由表面相交处(5,0,0)，点2位于虚拟控制面的中部(4.33,0.00,2.50)，点3位于虚拟控制面的底部(0,0,5)。图 6表明随着场点远离自由面，记忆速度势的变化逐渐趋于平缓，此与记忆速度势代表物体兴波部分速度势的物理意义相一致；且记忆速度势随时间出现一个较大的脉冲峰以后，便开始迅速衰减，此与修正Laguerre函数中的指数函数部分相符合，初步说明了利用Laguerre函数进行记忆速度势时间变量级数逼近的可行性。

根据图 6中的3组数据，本文对随时间变化最为剧烈的点1进行了Laguerre函数逼近研究。在不引入放缩系数的情况下，函数逼近随级数项数的变化如图 7所示，绝对误差随项数的变化如图 8所示。由图可知，随着逼近级数项数的增加，级数解越趋近于原数据，但是收敛效果随项数的增加改善较为缓慢，例如当级数项数增加一倍时，逼近绝对误差缩小不到一倍，而计算求解矩阵的维数却大大增加。

为了增加Laguerre级数逼近的收敛速度，同时又不增加计算机求解负担，引入了放缩系数的概念，如式(21)所示。由图 10可知，随着放缩系数的减小，相对误差逐渐减小，当放缩系数等于0.2时，级数项取为15的结果要明显优于不取放缩系数的级数项取为20的结果，如图 11所示。由此可以确定，通过引入放缩系数，可以利用较少的项数来逼近函数的真解，从而确保了利用Laguerre函数进行记忆速度势时间变量进行级数逼近的经济性和可行性。

3 结论

本文通过引入虚拟控制面，将流场划分为內域和外域，并根据速度势及其法向导数在控制面上连续的条件建立了时域快速匹配边界元模型。此外，提出了利用Laguerre函数对记忆速度势的时间轴进行级数逼近，以解析计算时间卷积积分的思路，并通过引入一放缩系数，大大加快了级数的收敛速度。通过对计算结果的分析，得到以下结论：

1) 通过与解析解的对比表明本文建立的内外场时域快速匹配边界元法具有较高的数值计算精度，且由求解矩阵可知，时域快速匹配边界元法与直接边界元法相比具有较高的计算效率，为研究浮体和波浪相互作用下流场的特性提供了一种精确、高效的方法；

2) 自由面划分为辐射状网格，控制面距离物面为5倍的物体的特征尺度，时间步长取为0.01 s时，时域快速匹配边界元法便能给出令人满意的计算结果；

3) 代表浮体兴波之后流场特性的记忆速度势在出现一个很大的峰之后，随时间迅速衰减，当t=2.5 s后，基本趋于平缓；

4) 利用Laguerre函数对记忆速度势的时间变量做级数逼近时，随着级数项数的增加，级数逐渐收敛，但级数收敛速度较为缓慢；

5) 分析Laguerre函数及其正交权函数的性质，引入一放缩系数，数值结果显示，放缩系数可大大加快级数的收敛速度，N=15，a=0.2的逼近结果要明显的优于N=20，a=1时的结果,保证了时域快速解析匹配边界元的可行性与高效性。