开源边界元三维辐射绕射水动力计算程序Nemoh及在此基础上的二次开发

本文摘要（由AI生成）：

本文主要介绍了Nemoh，一个基于边界元法的三维辐射绕射水动力计算程序。Nemoh自2014年开源以来，主要应用于海上波浪能分析领域及其他相关领域，目前有超过900个用户在持续使用该软件并进行相关研究工作。Nemoh能够实现的基本功能包括静水力计算、F-K力的计算与输出、绕射波浪力的计算与输出、辐射阻尼与附加质量的计算与输出、Kochin Function的计算与输出、单元波浪压力的计算与输出等。Nemoh的运行包括三个步骤：预处理、求解和后处理。Nemoh目前存在局限性，如后处理较弱，但最大的优势是已经将最复杂的问题（一阶速度势以及波浪载荷的计算）解决。用户可以在此基础上结合自身特定要求进行相关的开发活动。未来Nemoh的开发重点将集中在提高求解效率和实现更多功能。

目前在国内船舶与海洋工程工业界，尚无任何一款国产的基于势流理论的边界元水动力求解程序得到广泛应用和认可。长期以来，国内相关专业分析领域被WAMIT、AQWA、Hydrostar以及WADAM所占据。在当前严峻的形势下，立足自身，开展相关专业软件的研究与应用实践的需求十分迫切。

本文针对开源边界元三维辐射绕射水动力计算程序Nemoh的基本理论，程序组成，使用方法进行介绍。对本人围绕Nemoh进行的一些开发进行介绍，以供相关研究人员参考。

全文约6000字，图42个，阅读时间约为15分钟。

01

—

Nemoh是基于边界元法(Boundry Element Method,BEM)计算势流三维辐射绕射问题的水动力计算程序，由Ecole Centrale de Nantes（南特中央理工大学）开发，距今已有30多年的历史，2014年1月部分代码开源，也是目前唯一一个开源的BEM水动力计算程序[1]。

Nemoh自2014年开源以来主要应用于海上波浪能分析领域及其他相关领域，据称，目前有超过900个用户在持续使用该软件并进行相关研究工作[2]。

目前版本的Nemoh能够实现的基本功能主要包括：

1. 静水力计算，包括排水量、浮心、漂心以及静水刚度；
   

2. F-K力的计算与输出；

3. 绕射波浪力的计算与输出；

4. 辐射阻尼与附加质量的计算与输出；

5. Kochin Function（科钦函数）的计算与输出；

6. 单元波浪压力的计算与输出；

7. 流场关注点的波面计算与输出；

8. 脉冲响应函数（Impulse Response Function，IRF）的计算��输出。

Nemoh的理论基础[1][3]

对公式不感兴趣可以跳过本节

2.1 边界条件

Nemoh基于传统的势流边界元法。假设流体无旋，则速度势可以用下式表达：   

          （2.1）

考虑到流体不可压缩，则：

             （2.2）

则边界条件可以表达为：

  （2.3）

    (2.4）

      (2.5）

     (2.6）

其中，f(M)为标量函数。m0为波数，满足色散关系（2.7）。

        （2.7）

式子（2.3）~（2.6）分别对应：自由表面边界条件，海底边界条件，物体表面运动边界条件，辐射边界条件。

2.2 基本表达式

入射势的表达式为：

    （2.8）

β为入射波波浪方向。

有限水深条件下，f0表达式为：

          （2.9）

根据格林定理，速度势可以表达为：

    （2.10）

其中G为线性化的格林函数，σ（M’）为源分布，由下式表达：

  （2.11）

有限水深条件下的格林函数：

       （2.12）

其中：

      （2.13）

其中：

其中：

格林函数的梯度为：

           （2.14）

有限水深条件下，速度势可以表达为：

  （2.15）

此时，科钦函数（Kochin Function）H(α)表达式为：

      （2.16）

浸入水中的物体湿表面通过三角形或四边形平面单元进行离散，采用常数面元法进行求解，此时速度势和其法相导数认为在湿表面的网格单元上是定常的。

2.3 常数面元法

以N1代表物体离散的单元数目，对式子（2.10）、（2.11）进行离散进而得到以下线性表达式：

            (2.17)

其中Kmn和Smn称为影响系数，fn为对应面元n型心位置On的法相速度分布。

影响系数由下列式子求出：

        (2.18)

影响系数可以由两部分的和求得：

i.  对于1/MM'的积分，M'为面元单元上的一点，M为流场中的一点，可以通��解析方法进行求解；

ii.  对M'进行θ在-π/2至-π/2范围内进行积分。该部分内容可以通过理论求解，也可以通过二维参数构成的函数查询表插值得来，Nemoh采用的是参数计算查询表插值的方法。

Nemoh的运行

Nemoh由三个子程序组成，子程序及其对应功能分别是：

1. mesh：用于模型检查和静水力计算；

2. preProcessor：用于模型检查和F-K力的计算；

3. Solver：三维辐射绕射计算；

4. postProcessor：将计算结果整理输出。

3.1 Nemoh模型文件

Nemoh的水动力模型文件后缀名为.dat，文件组成为节点及对应编号以及单元及对应节点编号。一个典型的Nemoh模型文件格式如下图所示。模型文件第一行要标记模型对称性。Nemoh可以考虑模型关于XOZ平面对称。节点要有编号，节点描述完毕要有结束行。单元可以为四边形或者三角形单元，单元描述完毕需要有结束行标记。

图3.1 典型Nemoh模型文件格式（上图为节点编号及坐标，下图为单元节点编号）

Fig3.1. A TypicalNemoh Model File（Node Number and Coordinate Value，Top Figure；Node Numbers Formed  Elements，Bottom Figure）

3.2 运行Mesh.exe

dat模型文件可以用水动力计算，但如果需要进行静水力计算（运行Mesh.exe），则使用者需要单独编写一个无后缀名的模型文件，该文件包含单元数、节点数、节点坐标值以及组成单元的节点编号，如下图所示。

图3.2 用于Nemoh Mesh运行的模型文件（上图为节点数、单元数及节点坐标，下图为单元节点编号）

Fig3.2 A TypicalNemoh Model File for Mesh run（Node  Number，Element Number and Node

Coordinate Value，Top Figure；Node Numbers Formed  Elements，Bottom Figure）

在运行Mesh.exe之前需要编写一个名为Mesh.cal的输入文件，该文件包括模型文件名、对称性、重心位置、单元数、水线相对位置、缩比、外部流体密度以及重力加速度，如图3.3所示。

图3.3 用于Mesh.exe运行的Mesh.cal文件

Fig3.3 Mesh.cal for Running Mesh.exe

通过命令提示符（CMD）运行Mesh.exe，正常运行会在CMD窗口提示相关信息，如下图所示。

图3.4 Window系统下成功运行Mesh.exe

Fig3.4 Run Mesh.exe under Microsoft Windows

Mesh.exe运行完毕后会在Mesh文件夹中生成以下文件:

1. .dat根据无后缀名的模型描述文件重新生成的水动力计算模型文件；
   

2. .tec可以使用Tecplot打开查看的模型文件；

3.  _info.dat模型的节点和单元数目；

4.  GC_hull.dat指定的重心位置；

5.  Hydrostatics.dat静水力计算数据，包括浮心、重心、漂心、排水体积以及水线面面积；

6.  Inertia_hull.dat程序计算的模型惯性矩信息；

7.  KH.dat程序计算的模型静水刚度。

经过作者测试，目前版本的Mesh.exe存在bug。当模型曲面较复杂时，改变单元数，Mesh.exe计算的静水刚度会发生变化。建议使用者在使用Mesh.exe时尽量采用比较密的网格进行计算，以减少误差的发生。

3.3 使用Nemoh进行水动力计算

在进行水动力计算之前需要编制Nemoh.cal文件。该文件包括以下主要内容：

1. 用户定义的流体密度、重力加速度、水深以及波浪生成参考点；

2. 计算的体的数目；

3. 对应体的模型文件名、节点数、单元数以及计算自由度；

4. 计算最大、最小频率以及二者之间的计算频率数；

5. 计算最大、最小波浪方向以及二者之间的计算浪向数；

6. 后处理设置，包括IRF计算参数、科钦函数计算参数、是否保存单元波浪压力文件、是否计算自由表面波高以及计算范围等。

典型Nemoh.cal文件内容如下图所示。

图3.5 典型Nemoh.cal文件内容

Fig3.5 A Typical Nemoh.cal‘s Input Data

另外还需要准备一个input.txt文件，文件中只要输入0即可。input.txt保存的是程序数值解法的相关设置。在官方例子中input.txt中有直接求解法（以0标记）和广义最小残值法（Generalized Minimal Residual，GRESS）（以1标记）以及GRESS的相关设置。但打开Nemoh的代码可以发现，程序并没有提供GRESS解法，换句话说，目前版本的Nemoh只提供直接求解法，因而input.txt中只需输入0即可。

Nemoh水动力计算的运行顺序为：

1. 运行preProcesser.exe

程序读入模型文件，随后在mesh文件夹中生成L10.dat、L12.dat、Integration.dat、Freesurface.dat、Kochin.dat、Mesh.tec文件并在读入的模型目录下生成Normalvelocities.dat。L10.dat、L12.dat为重新编排的模型文件。L10保留模型的单元信息，包括对称性、型心、法线方向、面积等。L12.dat对读入的模型文件进行格式化整理。Integration.dat保留模型法线方向信息和面积信息，用于压力的积分计算。Freesurface.dat、Kochin.dat用于保存自由表面输入信息和科钦函数的相关输入数据。Normalvelocities.dat保存对应各个求解工况的模型单元速度向量。

图3.6 运行preProcesser.exe

Fig3.6 Run preProcesser.exe

程序将需要进行计算的自由度、波浪频率、波浪方向、科钦函数等参数保留在results文件夹中的index.dat文件中。

程序会根据波浪速度势在各个单元型心位置的分布进行入射压力计算并进行积分，给出F-K力，并将结果保存在results文件夹下的FKforce.dat和FKforce.tec文件中。

2. 运行solver.exe

程序读入Nemoh.cal、L12.dat、Normalvelocities.dat和Integration.dat。程序分配内存，生成计算工况，读入单元速度向量，进行辐射速度势和绕射速度势的求解。求解完毕后进行压力计算，读入Integration.dat进行积分，最终给出对应辐射压力（实部附加质量，虚部辐射阻尼）和绕射波浪压力（实部幅值，虚部相位）。所有结果都将保存在results文件夹中的Forces.dat文件中。

图3.7 运行solve.exe

Fig3.7 Run solve.exe

3. 运行postProcesser.exe

程序读入Forces.dat和FKforce.dat文件，将附加质量和辐射阻尼结果保存在CA.dat、CM.dat以及RadiationCoefficients.tec中。

图3.8 运行postProcesser.exe

Fig3.8 Run postProcesser.exe

绕射力保存在DiffractionForce.tec文件中。总的波浪力保存在ExcitationForce.tec中。

如果计算科钦函数，则相关计算结果保存在Kochin  XXX.dat文件中，XXX为对应波浪方向。如果输出单元波浪压力，则数据会保存在pressure.  XXX.dat中，XXX为对应单元号。如果计算了自由表面波面情况，则结果数据会保存在freesurfaceXXX.dat文件中，XXX为对应点。

至此，Nemoh完成全部工作

Nemoh目前的局限性与优势

从软件开发历史来看，Nemoh的开发源自Géard Delhommeau的博士论文[4]。论文中对于三维辐射绕射的相关理论和数值求解方法进行了非常详细的介绍。

在此基础上，Géard Delhommeau开发了耐波性计算程序AQUADYN（零航速）和AQUAPLUS（考虑航速）。Nemoh程序与这两个软件有紧密的渊源，其格林函数及数值解法完全脱胎于Géard Delhommeau的博士论文。Nemoh作为目前唯一一个开源的BEM水动力求解程序，这本身已经是非常有意义的了。

从应用角度来看，目前版本的Nemoh还有很多不完善的地方：

1. 代码的可读性和可维护性较差

Nemoh有相当多的代码是Fortran77格式，在影响矩阵的求解中有大量的GOTO语句，使得后来者比较难以理解。Nemoh的程序中有很多数值计算函数是作者编写的，现在已经可以通过库函数来解决。

2. 程序架构不利于并行处理

Nemoh的核心求解程序完全是线性特性的程序架构，存在大量的GOTO语句以及程序嵌套，必须对源代码进行重构才能实现多线程并行处理。

3.求解速度较低

Nemoh目前仅提供直接求解法。对1000个单元左右的模型进行20个周期的计算耗时在15分钟左右，常规的商业软件对于同样的问题求解时间在2~5分钟左右。

4.后处理较弱

实际上目前版本的Nemoh是不存在常规意义上的后处理功能。用户需要围绕Nemoh进行大量的后续开发工作来实现诸如运动幅值响应算子（Response Amplitude Operater， RAO）、远场波浪力（基于运动RAO和科钦函数）、近场波浪力（基于速度势、运动RAO、速度势梯度等）、模型显示、波面模拟，乃至截面载荷计算和整体弯矩/剪力计算等功能。

相对而言，Nemoh已经将最复杂的问题（一阶速度势以及波浪载荷的计算）解决，这是它最大的优势。用户可以在此基础上结合自身特定要求进行相关的开发活动。另外，由于速度向量文件已知，用户可以通过重写该文件来实现更广泛的应用，诸如浮体的广义自由度的应用分析等。

围绕Nemoh进行的开发工作

Capytaine: a Python-based distribution of Nemoh

Matthieu Ancellin重新编写了Nemoh的求解程序，将这部分嵌入python中形成了Capytaine[5]。目前Capytaine的最新版本为1.2（2020.4）。

Capytaine嵌入通过Fortran重写的Nemoh格林函数部分。同时将meshmagick和python的科学计算库进行整合。目前能够实现的功能包括：附加质量、辐射阻尼、绕射波浪力以及入射波浪力的求解；python API接口；对格林函数的编写求解进行了优化；基于OpenMP利用多线程加速求解；实现运动与自由波面的三维模拟等。

相比于Nemoh，Capytaine的代码可读性和可维护性更好。Capytaine将python和fortran通过f2py进行跨平台编程，将耗费资源较高的求解内容通过Fortran计算，通过python来实现前后处理，有更好的功能分工和跨平台性能。

未来Capytaine将会把开发重点放到程序求解效率的提高以及程序通用性。

Capytaine代码已开源，可在Github进行下载。

Wave Analysis Response Program, WARP[3]

Wave Analysis Response Program，WARP，著作权登记号2018R841654。是作者围绕Nemoh进行开发的软件，直白一点的讲，是本人给Nemoh套了个壳并对原生Nemoh做了部分修改。

WARP目前有两个版本：Demo版和开发版，Demo版用于原生Nemoh的简单前处理与运行文件的编写，开发版中对原生Nemoh进行了部分代码优化，并添加一些后处理和数据计算功能。整合了meshmagick以实现模型显示与修复以及静水力计算。

图5.1 WARP Demo界面

Fig5.1 WARP Demo GUI

1. WARP的前处理

基于工作习惯，目前模型通过ANSYS APDL进行建模，将节点和单元输出后，WARP将其组装为.dat模型格式，将模型进行显示与修复。

图5.2 WARP 1.23 显示船体面元模型

Fig5.2 WARP 1.23 A ISO View of MPV Vessel Mesh

用户通过输入参数来定义Nemoh.cal文件，程序通过调用Nemoh来实现波浪力与附加质量和辐射阻尼的计算。

图5.3 WARP 1.23工况设置

Fig5.3 WARP 1.23 Case Setting

2. RAO与远场波浪力的求解

用户输入浮体质量信息后，程序调用Nemoh计算的数据来进行浮体RAO的计算。

程序调用Nemoh计算的科钦函数和计算完毕的RAO数据进行远场波浪载荷的求解。作者对原生Nemoh进行了部分修改，现在可以通过WARP计算远场法计算的艏摇二阶定常波浪力。

图5.4 WARP 1.23 RAO参数输入界面

Fig5.4 WARP 1.23 RAO Calculation Setting

3. 简单的频域运动分析

用户可以通过定义关注点和波浪谱（目前仅支持JOSNWAP谱）进行波频频域运动分析。对于频域系泊分析，可以通过WARP_Fmo来实现。

4. 数据曲线显示

WARP可以将主要的计算结果（附加质量、辐射阻尼、一阶波浪力、远场波浪力、重心位置运动RAO及相位）进行显示。对模型能够进行显示查看。

图5.5 WARP 1.23 显示升沉附加质量和辐射阻尼曲线

Fig5.5 WARP 1.23, Added Mass and Radiation Damping of Heave Motion

图5.6 WARP 1.23 显示横摇波浪力矩曲线及相位

Fig5.6 WARP 1.23, Roll Wave Force Moment and Phase Angle Curve

图5.7 WARP 1.23 显示横摇运动RAO幅值曲线及相位

Fig5.7 WARP 1.23, Roll Motion RAO and Phase Angle Curve

图5.8 WARP 1.23 远场二阶定常纵荡力曲线

Fig5.8 WARP 1.23, Far-Field Wave Drift Force Curve

目前WARP开发版版本号为1.23。未来WARP的开发重点将集中在以下几个方面：

1. 对Nemoh源代码进行重构，实现多线程并行计算；

2. 对Nemoh源代码进行重构，将运动RAO、远场算法以及近场算法进行集成；

3. 对后处理进行扩展，添加波面与运动显示以及面元压力分布显示；

4. 评估集成或内核转为Capytaine的可行性。

总结与展望

目前在国内船舶与海洋工程工业界，尚无任何一款国产的基于势流理论的边界元水动力求解程序得到���泛应用和认可。长期以来，国内相关专业分析领域被WAMIT、AQWA、Hydrostar以及WADAM所占据。在当前严峻的形势下，立足自身开展相关专业软件的研究十分迫切。立足于开源软件进行开发不失为一种快速解决目前难题的捷径。

本文针对开源边界元三维辐射绕射水动力计算程序Nemoh的基本理论，程序组成，使用方法进行介绍。对本人围绕Nemoh进行的一些开发进行介绍，以供相关研究人员参考。

围绕Nemoh可以开展的工作有很多，但核心还是在于提高计算效率和扩展计算功能两大方面，这也将是未来相关开发工作的方向。

限于篇幅，本文不对Nemoh的计算精度进行对比分析，如有条件将以单独文章对这部分内容进行介绍。

欢迎相关领域研究人员和对本文所述开发工作感兴趣的老板共同交流探讨。

本号也将开设专栏对Nemoh的使用和在其基础上的开发和应用进行介绍，欢迎关注。

参考文献：

[1]Aurélien Babarit, Géard Delhommeau. Theoretical and numerical aspects of the opend sourse BEM solver NEMOH[C]. EWTEC 2015.

[2] Aurélien Babarit, Géard Delhommeau. The Presentation of “Theoretical and numerical aspects of the opend sourse BEM solver NEMOH”.EWTEC 2015.

[3]WARP V1.00 软件帮助[R].2018.6.

[4]Géard Delhommeau. Les Problemes de Diffraction-Radiation et de Resistance de Vahues: Etude Theorique et Resolution Numerique par la Methode des Singularites[D]PhD Thesis, École Nationale Supérieure de Mécanique de Nantes, 1987.

[5]Capytaine: a Python-based distribution of Nemoh,

   https://ancell.in/capytaine/latest/index.html