沈王刚，郑尧坤，林志良

(上海交通大学 船舶海洋与建筑工程学院 海洋工程国家重点实验室，上海 200240)

0 引 言

研究潜体近水面航行时的兴波现象及其力学特性，一直是一项具有重要工程意义的课题。尤其近几年各类附加水翼新概念船型的大量涌现，使得潜体兴波研究显得更为重要。由于计算效率的限制，以往对潜体兴波研究主要采用定常方法。如Farell改进了Havelock源格林函数，推导得到潜航回转椭球体在Neumann-Kelvin问题下的半解析理论解[1]；Doctors和Beck研究了潜体在Neumann-Kelvin问题下，采用Havelock源计算的收敛性[2]。随着时代的发展，时域方法以其可以计算非定常兴波，得到运动物体实时兴波信息，便于仿真模拟和后处理等优点，吸引了越来越多学者的注意。

自从Hess和Smith[3]第一次应用面元法求解三维物体绕流问题以来，面元法在船舶与海洋工程领域得到了广泛的应用。其中Rankine源形式简单，具有较高灵活性，且易于推广到非线性计算，已经有不少学者尝试使用Rankine源处理兴波问题[4–5]。

然而，传统的面元法在求解问题时，需要计算每个面元与配置点之间的相互作用，生成一个N阶方阵（N为离散面元个数）。一方面储存这个N阶稠密矩阵需要大量的内存空间，另一方面求解该矩阵也将耗去大量时间。例如采用高斯消去法直接求解，时间复杂度将达到O（N3），即使使用迭代法求解，时间复杂度也有O（N2）。计算效率上的缺陷将面元法的应用限制在了小规模流场，更遑论时域计算时需要进行多次时间步进迭代，计算规模将被进一步压缩。为克服传统面元法这一局限性，引入快速多极子法，可同时减小存储量与时间复杂度。

快速多极子法最早由Rokhlin[6]于1985年提出，用来加速二维势流问题的求解。此后Liu[7]，Nishimura[8]及Yoshida[9]等将快速多极子法和边界元法结合，并对其进行深入研究。多极子法与边界元法的结合已被证实可以大幅度提高其计算效率，适用于求解大规模势流问题[10–11]。

本文结合多极子法与面元法，对三维潜体的兴波特性进行时域研究，将所得结果与传统方法结果进行比较，证明多极子面元法高精度与高效率的特性。

1 潜体兴波问题描述

考虑如图1所示的计算流场，O-xyz为三维笛卡尔直角坐标系，xy平面与自由面平行，潜体沿x轴方向以速度U运动，z轴垂直向上，整个流场边界由潜体物面SB和自由面SF组成。潜体的特征长度为l，特征宽度为d，潜体中心距自由表面深度为h。对于时域兴波问题，自由面条件中对流项起到主要作用，足够远处将自由水面截断，截断边界上的反射对势流解的影响并不敏感[12]。因此本文中采用截断有限自由表面代替无限自由表面，其长为L，宽为D。假设流体无压、无粘并且无旋，流场中存在一速度势，满足拉普拉斯方程：

在潜体表面SB上满足不可穿透条件，即

在本文时域兴波问题中，假设潜体由静止快速启动至匀速运动状态，且在自由面SF上初始条件（t=0）为：

图1 计算流场示意图Fig. 1 The schematic of computational flow field

2 传统面元法求解

首先将整个流场的边界划分成N个离散面元，其中前NB个为三角形常数面元，用于表征潜体表面；后NF个为四边形常数面元，用于表征自由水面，则有N=NB+NF。第i个面元的的面积为Si，其上均匀分布强度为σi的面源。此时，可得到流场中任一点的诱导速度势

诱导速度

此外，选取每一个常数面元的中点为配置点，在每一个配置点上满足已知的边界条件，例如在t=0初始时刻，潜体物面上的配置点应满足式（2），而自由表面上的配置点应满足式（5）中的据此可得到如下线性方程组，即

常数面元法的优点在于公式推导简单，易于应用；但为了提高精度，需要增加离散面元数目，使得计算规模增大、计算效率降低。此外，不少学者尝试使用高阶面元法求解势流问题，如Liu Y H[15]，Romate JE[16]等。尽管计算精度和效率有所提升，高阶面元法仍无法有效克服传统面元法求解大规模流场问题的缺陷。因此，本文采用快速多极子方法加速传统常数面元法计算，求解大规模流场时域兴波问题。

3 多极子面元法

快速多极子法与传统面元法的结合过程和其与传统边界元法结合过程相同，详细内容可参考文献[7 – 9]。本节将简要介绍快速多极子方法的主要思想、基本公式和与面元法结合后的计算步骤。

3.1 多极子法主要思想

以粒子在粒子群中受力为例，两两直接计算再求和策略的计算量级是O（N2），其中N表示粒子数量。快速多极子法将所有粒子按其空间位置分成不同的集合（如图2中的A和B），首先将粒子对外的作用力汇聚到集合内一点（如图2中的①过程），其次计算集合与集合之间的相互作用力（如图2中的②过程），最后将集合所受到的外界作用力分配到集合内每一个粒子（如图2中的③过程）。对于相距很近的粒子，相互间作用力直接计算得到；而对于相距较远的粒子，相互间的作用力则采用上述多极子方法展开得到，从而将计算量级减少至O（N·LogN)[7–9]。其中区分粒子间距关系，可采用树状结构划分策略[17]加以实现。在传统面元法中，主要过程在于计算边界离散面元间的相互作用，因此可以采用快速多极子方法来提高计算效率。

图2 多极子计算示意图Fig. 2 The schematic of multipole calculation

图3 多极子展开关键点Fig. 3 Related points in multipole expansion

3.2 多极子法基本公式

关于多极子法公式具体推导请参见文献[7-9]，这里以面元法中三维Rankine源诱导速度势核函数多极子展开为例，其展开形式如下：

得到核函数的展开形式后，即可得到积分方程的展开形式如下：

对于核函数的局部展开系数则定义如下：

3.3 多极子面元法计算步骤

在此只简要介绍多极子面元法的计算步骤，具体请参见文献[18]。

步骤1与传统面元法一样对边界进行单元离散，布置面源，选取配置点。

步骤2划分树状结构。以二维问题为例，利用逐层嵌套的正方形来划分空间树状结构。首先利用一个足够大的正方形包围整个边界S，并定义其层级数为0级。考察第L层级的正方形单元C，将其均分为4个小正方形单元，层级为L+1，正方形单元C是它们的父单元。若某离散边界单元的中点位于小正方形内，则认为其属于该小正方形单元。对于4个小正方形单元，分成3种情况。若该小正方形单元内离散边界单元数为0，则将它舍弃，不计入树状结构内；若该小正方形单元内离散边界单元数大于0，且不大于事先给定的常数p时，将其定义为叶单元，不再往下划分；若该小正方形单元内离散边界单元数大于常数p时，则继续将其往下划分。如此嵌套割划，直到最后一层均为叶单元，或达到规定最大层数则停止。对于正方形单元之间的空间关系，定义如下（见图4）：若2个正方形单元共用一个角点，则定义它们为相邻单元；若2个正方形单元不相邻，但其父单元相邻，则定义它们为相隔单元；若2个正方形单元的父单元也不相邻，则定义它们为远距单元。

图4 正方形单元相互关系Fig. 4 Relations between square cells

步骤3向上传递。首先通过式（24）或式（25）计算各个叶单元的多极子动量矩，其次利用式（26）将叶单元的多极子动量矩转移到其父单元中。以此类推，多极子动量矩不断上聚，直至第2层。

步骤4向下传递。对于相隔单元，通过式（28）计算其作用；对于远距单元，则利用式（29）将父单元通过得到的局部展开系数继承到子单元，以此来计算其作用。如此类推，直到得到每个叶单元的局部展开系数。向上传递和向下传递过程参见图5。

图5 向上和向下传递过程Fig. 5 Upward and downward passes

步骤5近场远场作用叠加。对于某叶单元C中的配置点，利用传统面元法的直接作用计算该叶单元及相邻单元中其他离散面元对其的作用，利用式（27）将局部展开系数从单元中心转移到配置点以计算相隔或远距单元中离散面元对其的作用。将近场和远场的作用叠加，就得到了整个边界所有离散面元对配置点的作用。

步骤6利用GMRES求解式（10），若迭代误差小于设定值，则终止运算，否则返回第3步进行下一步迭代。本文中迭代误差设定值均取为10–4。

4 计算结果与分析

4.1 近水面回转椭球体匀速直航兴波问题

本节将以近水面回转椭球体匀速直航产生的兴波波形和兴波升力、阻力为切入点，阐述如何选取合适的自由面计算域和自由面及物面网格密度，验证使用多极子面元法计算潜体兴波问题的可行性，并将多极子面元法与传统面元法进行计算效率对比。

计算模型为回转椭球体，长轴长l=1.0 m、短轴长d=0.2 m的回转椭球在自由水面下（球心距静水面距离为h=0.16 m）从静止状态突然启动，达到常速U=1.566 m/s并沿轴正向匀速直航，Froude数为浸深比为椭球表面网格划分如图6所示，其中图6（a）椭球面元数为528，图6（b）椭球面元数为1 304，网格划分时均为两端密，中间疏；图6（c）椭球面元数为2 640，网格划分时两端与中间均加密。

图6 椭球网格划分Fig. 6 The mesh generation of elipsoid

自由面的网格密度和时间步长会直接影响到计算量、波面的光滑程度以及兴波水动力学参数的准确性，网格和时间步长的选取参照的是Dommermuth等人通过大量数值试验得到的经验公式[19]：

4.1.1 自由面静网格

自由面为一足够大的长L=20 m，宽D=6 m的矩形平面，对称于坐标系原点。椭球球心初始坐标为x=–7.5，y=0，沿轴正向匀速直航。

椭球网格划分如图6（b）所示。自由面网格密度有3套，其中粗网格选取为100×30，网格间距为中网格选取为150×45，网格间距为细网格选取为200×60，网格间距时间步长间隔为0.04 s，步进次数200次，总模拟时间8 s。

图7为t=8 s时波面图，经过对比可发现，粗网格与细网格所描绘出的波面波形图及等高线图基本一致，但粗网格仍有一些波面锯齿，而细网格的波面则光滑许多。对比图8亦可发现，3种网格描绘出的自由面中剖线二维波形图在第一个波谷处有比较明显的差别，除此之外中网格和细网格的结果完全贴合，而粗网格则与其他2种网格差别较大。

图7 t=8 s时波面图Fig. 7 The wave surface at t=8 s

其次考察兴波导致的椭球体升力、阻力。图9（a）为升力系数随时间变化情况，图9（b）为兴波阻力系数随时间变化情况。图中圆点直线为Farell解[1]，3种网格所计算得到的力学系数都能收敛到Farell解附近。然而粗网格结果，特别是其升力系数的波动很大，毛刺明显；中网格和细网格的毛刺较粗网格有所减少，但仍然比较明显。此外，粗网格计算得到的升力系数收敛值和其他2种网格比较吻合，但兴波阻力收敛值有比较明显的差别，显示粗网格密度并不够，计算精度低；同时随着自由面网格不断加密，其计算结果也将逐渐收敛。

图8 t=8 s时自由面二维波形图（y=0）Fig. 8 The 2D waveform of free surface at t=8 s (y=0)

图9 兴波水动力参数随时间变化情况Fig. 9 The change of wave making hydrodynamic parameters over time

使用静网格固然能描述潜体的整个尾后兴波波面，但其取的自由面计算域过大，网格数量过多，时间效率低下；此外有大片自由面区域在长时间内兴波影响很小，是一种计算资源浪费。为了能提高计算效率，减少计算资源浪费，使用能随潜体一起前进的动网格来进行兴波时域计算。

4.1.2 自由面动网格与网格密度

自由面为一个长L=8 m，宽D=5 m的矩形平面，初始时刻对称于坐标系原点。椭球球心初始坐标为x=–7.5，y=0，沿轴正向匀速直航。

椭球网格划分有3套（见图6）。自由面网格密度选取为80×50，网格间距时间步长间隔为0.04 s，步进次数200次，总模拟时间8 s。在每一次步进后，若椭球中心与网格中心的距离大于，则将网格向轴正向移动一个网格间距，相当于消去最后一排网格，并在第一排网格前方再加上一排网格，其初始条件如式（5）所示。

从图10中可以看到，虽然动网格的网格密度与静网格中细网格密度相同，但其计算结果中的毛刺现象却有大幅度改善，究其原因可能是静网格所取的自由面计算域过大，潜体相对网格的位置变化明显，在网格密度不够绝对密的情况下毛刺就会比较多，从图9中也可看出，当自由面静网格加密后，毛刺现象有所改善。此外，1304面元的椭球网格2和2640面元的椭球网格3所计算得到的升力曲线完全一致，兴波阻力曲线有略微差别，但只在1.8%左右；而528面元的椭球网格1则与其他2种网格计算结果差别较大，因此可认为随着椭球网格数量的增加，计算结果将逐渐收敛，同时1304面元的椭球网格2精度已经足够。

此外，考察自由面横向网格密度变化率q的影响，q=1代表横向网格均匀分布，q＞1代表中纵线附近网格密，横向两端网格疏。表1显示了不同横向网格密度变化率下计算所得结果，其中2种网格布置的横向网格数相同，均为50个。由表中可看出横向网格变化率对计算结果的影响很小，几乎可以忽略不计。

4.1.3 计算效率对比

图10 不同物面网格下兴波水动力参数随时间变化情况Fig. 10 The change of wave making hydrodynamic parameters over time in different object surface mesh

表1 不同横向网格密度变化率下兴波水动力学系数Tab. 1 The wave making hydrodynamic parameters in different changing rate of transverse mesh density

表2 多极子面元法和传统面元法单步计算效率比较Tab. 2 Comparison of single-step computational efficiency between FMBEM&BEM

多极子法对传统面元法的改进，首先体现在内存使用量上，传统面元法需要占用O（N2）的内存，普通台式计算机无法驾驭面元数量过万的问题，但多极子面元法只需要占用O（N）的内存，使得普通台式计算机也可以计算面元数量达到数百万甚至数千万的大规模问题。表2和图11是分别采用多极子面元法和传统面元法进行时域单步计算的效率对比，计算环境为1.9 GB内存和2.50 GHz*2处理器。从图表中可以看出，当面元数较小时，使用多极子面元法和传统面元法所需要的计算量相近，甚至当面元数只有数百个时，传统面元法还要略优于多极子面元法；然而随着面元数的增长，多极子面元法所需的计算量近似线性增长，而传统面元法所需的计算量却是以三次方的速度增长，当面元数上万时，使用传统面元法求解不仅在内存上受到限制，在计算时间上也显得不现实。此外，由于时域问题较定常问题要多一个时间维度，因此利用多极子面元法来计算时域问题的可行性和高效性不言而喻，具有巨大的应用潜力。

图11 多极子面元法和传统面元法单步计算效率比较Fig. 11 Comparison of single-step computational efficiency between FMBEM&BEM

4.2 不同航速及潜深下椭球定常兴波

本小节将以近水面椭球匀速直航达到稳定状态后所产生的升力和兴波阻力为切入点，说明不同航速和下潜深度对其定常兴波的影响。计算模型与上节一样为回转椭球体，长轴长l=1.0 m、短轴长d=0.2 m，改变航速和潜深以研究其定常兴波力学系数的变化情况。

兴波水动力系数如图12所示，其中图12（a）是升力系数随Froude数变化情况，图12（b）是兴波阻力系数随Froude数变化情况。与Doctor解[2]相比，本文方法所得结果在时吻合较好；在时，高Froude数情况下的升力系数和阻力系数，以及Fr=0.45～0.5段的阻力系数会有较为明显的差别，其余部分吻合较好。由图中可以看出，浸深比下升力系数和阻力系数的值总小于浸深比下的值，升力系数在左右达到最大，在Fr=0.55～0.6后变为负升力；阻力系数在左右达到最大，在低Froude数时接近零。此外，下升力系数和阻力系数最大值所对应的Froude数较略微后移。

5 结 语

从计算结果来看，本文所开发的快速多极子面元法可以较准确地预报近水面潜体时域兴波特性，定常兴波计算结果与已有研究结果吻合良好，证明了其合理性与可靠性；且相较于传统面元法，有效地克服了计算规模和效率的局限性，引进了运动网格提高计算精度和效率。本文方法将进一步改进和发展，以计算完全非线性兴波问题、船体兴波问题等，具有较大的应用前景。

图12 兴波水动力参数Fig. 12 The wave making hydrodynamic parameters