目前在國內船舶與海洋工程工業界，尚無任何一款國產的基于勢流理論的邊界元水動力求解程序得到廣泛應用和認可。長期以來，國內相關專業分析領域被WAMIT、AQWA、以及WADAM所占據。在當前嚴峻的形勢下什么情況下要計算節點坐標，立足自身，開展相關專業軟件的研究與應用實踐的需求十分迫切。
本文針對開源邊界元三維輻射繞射水動力計算程序Nemoh的基本理論，程序組成，使用方法進行介紹。對本人圍繞Nemoh進行的一些開發進行介紹，以供相關研究人員參考。
全文約6000字，圖42個，閱讀時間約為15分鐘。
01
—
Nemoh簡介
Nemoh是基于邊界元法( ,BEM)計算勢流三維輻射繞射問題的水動力計算程序，由Ecole de （南特中央理工大學）開發，距今已有30多年的歷史，2014年1月部分代碼開源，也是目前唯一一個開源的BEM水動力計算程序[1]。
Nemoh自2014年開源以來主要應用于海上波浪能分析領域及其他相關領域，據稱，目前有超過900個用戶在持續使用該軟件并進行相關研究工作[2]。
目前版本的Nemoh能夠實現的基本功能主要包括：
靜水力計算，包括排水量、浮心、漂心以及靜水剛度；
F-K力的計算與輸出；
繞射波浪力的計算與輸出；
輻射阻尼與附加質量的計算與輸出；
（科欽函數）的計算與輸出；
單元波浪壓力的計算與輸出；
流場關注點的波面計算與輸出；
脈沖響應函數（ ，IRF）的計算與輸出。
02
—
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）
此時，科欽函數（ ）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采用的是參數計算查詢表插值的方法。
03
—
Nemoh的運行
Nemoh由三個子程序組成，子程序及其對應功能分別是：
mesh：用于模型檢查和靜水力計算；
：用于模型檢查和F-K力的計算；
：三維輻射繞射計算；
：將計算結果整理輸出。
3.模型文件
Nemoh的水動力模型文件后綴名為.dat，文件組成為節點及對應編號以及單元及對應節點編號。一個典型的Nemoh模型文件格式如下圖所示。模型文件第一行要標記模型對稱性。Nemoh可以考慮模型關于XOZ平面對稱。節點要有編號，節點描述完畢要有結束行。單元可以為四邊形或者三角形單元，單元描述完畢需要有結束行標記。
圖3.1 典型Nemoh模型文件格式（上圖為節點編號及坐標，下圖為單元節點編號）
Fig3.1. A Model File（Node and Value，Top ；Node ， ）
3.2運行Mesh.exe
dat模型文件可以用水動力計算，但如果需要進行靜水力計算（運行Mesh.exe），則使用者需要單獨編寫一個無后綴名的模型文件，該文件包含單元數、節點數、節點坐標值以及組成單元的節點編號，如下圖所示。
圖3.2用于Nemoh Mesh運行的模型文件（上圖為節點數、單元數及節點坐標，下圖為單元節點編號）
Fig3.2A Model File for Mesh run（Node ， and Node
Value，Top ；Node ， ）
在運行Mesh.exe之前需要編寫一個名為Mesh.cal的輸入文件，該文件包括模型文件名、對稱性、重心位置、單元數、水線相對位置、縮比、外部流體密度以及重力加速度，如圖3.3所示。
圖3.3用于Mesh.exe運行的Mesh.cal文件
Fig3.3Mesh. Mesh.exe
通過命令提示符（CMD）運行Mesh.exe，正常運行會在CMD窗口提示相關信息，如下圖所示。
圖3.系統下成功運行Mesh.exe
Fig3.4Run Mesh.exe s
Mesh.exe運行完畢后會在Mesh文件夾中生成以下文件:
.dat根據無后綴名的模型描述文件重新生成的水動力計算模型文件；
.tec可以使用打開查看的模型文件；
_info.dat模型的節點和單元數目；
.dat指定的重心位置；
.dat靜水力計算數據，包括浮心、重心、漂心、排水體積以及水線面面積；
.dat程序計算的模型慣性矩信息；
KH.dat程序計算的模型靜水剛度。
經過作者測試，目前版本的Mesh.exe存在bug。當模型曲面較復雜時，改變單元數，Mesh.exe計算的靜水剛度會發生變化。建議使用者在使用Mesh.exe時盡量采用比較密的網格進行計算，以減少誤差的發生。
3.3使用Nemoh進行水動力計算
在進行水動力計算之前需要編制Nemoh.cal文件。該文件包括以下主要內容：
用戶定義的流體密度、重力加速度、水深以及波浪生成參考點；
計算的體的數目；
對應體的模型文件名、節點數、單元數以及計算自由度；
計算最大、最小頻率以及二者之間的計算頻率數；
計算最大、最小波浪方向以及二者之間的計算浪向數；
后處理設置，包括IRF計算參數、科欽函數計算參數、是否保存單元波浪壓力文件、是否計算自由表面波高以及計算范圍等。
典型Nemoh.cal文件內容如下圖所示。
圖3.5 典型Nemoh.cal文件內容
Fig3.5A Nemoh.cal‘s Input Data
另外還需要準備一個input.txt文件，文件中只要輸入0即可。input.txt保存的是程序數值解法的相關設置。在官方例子中input.txt中有直接求解法（以0標記）和廣義最小殘值法（ ，GRESS）（以1標記）以及GRESS的相關設置。但打開Nemoh的代碼可以發現，程序并沒有提供GRESS解法，換句話說，目前版本的Nemoh只提供直接求解法，因而input.txt中只需輸入0即可。
Nemoh水動力計算的運行順序為：
1. 運行.exe
程序讀入模型文件，隨后在mesh文件夾中生成L10.dat、L12.dat、.dat、.dat、.dat、Mesh.tec文件并在讀入的模型目錄下生成.dat。L10.dat、L12.dat為重新編排的模型文件。L10保留模型的單元信息，包括對稱性、型心、法線方向、面積等。L12.dat對讀入的模型文件進行格式化整理。.dat保留模型法線方向信息和面積信息，用于壓力的積分計算。.dat、.dat用于保存自由表面輸入信息和科欽函數的相關輸入數據。.dat保存對應各個求解工況的模型單元速度向量。
圖3.6 運行.exe
Fig3.6 Run .exe
程序將需要進行計算的自由度、波浪頻率、波浪方向、科欽函數等參數保留在文件夾中的index.dat文件中。
程序會根據波浪速度勢在各個單元型心位置的分布進行入射壓力計算并進行積分，給出F-K力，并將結果保存在文件夾下的.dat和.tec文件中。
2. 運行.exe
程序讀入Nemoh.cal、L12.dat、.dat和.dat。程序分配內存，生成計算工況，讀入單元速度向量，進行輻射速度勢和繞射速度勢的求解。求解完畢后進行壓力計算，讀入.dat進行積分，最終給出對應輻射壓力（實部附加質量，虛部輻射阻尼）和繞射波浪壓力（實部幅值，虛部相位）。所有結果都將保存在文件夾中的.dat文件中。
圖3.7 運行solve.exe
Fig3.7 Run solve.exe
3. 運行.exe
程序讀入.dat和.dat文件，將附加質量和輻射阻尼結果保存在CA.dat、CM.dat以及s.tec中。
圖3.8 運行.exe
Fig3.8 Run .exe
繞射力保存在.tec文件中。總的波浪力保存在.tec中。
如果計算科欽函數，則相關計算結果保存在 XXX.dat文件中，XXX為對應波浪方向。如果輸出單元波浪壓力，則數據會保存在. XXX.dat中，XXX為對應單元號。如果計算了自由表面波面情況，則結果數據會保存在.dat文件中，XXX為對應點。
至此，Nemoh完成全部工作
04
—
Nemoh目前的局限性與優勢
從軟件開發歷史來看，Nemoh的開發源自Géard 的博士論文[4]。論文中對于三維輻射繞射的相關理論和數值求解方法進行了非常詳細的介紹。
在此基礎上，Géard 開發了耐波性計算程序（零航速）和（考慮航速）。Nemoh程序與這兩個軟件有緊密的淵源，其格林函數及數值解法完全脫胎于Géard 的博士論文。Nemoh作為目前唯一一個開源的BEM水動力求解程序，這本身已經是非常有意義的了。
從應用角度來看，目前版本的Nemoh還有很多不完善的地方：
1. 代碼的可讀性和可維護性較差
Nemoh有相當多的代碼是格式，在影響矩陣的求解中有大量的GOTO語句，使得后來者比較難以理解。Nemoh的程序中有很多數值計算函數是作者編寫的，現在已經可以通過庫函數來解決。
2. 程序架構不利于并行處理
Nemoh的核心求解程序完全是線性特性的程序架構，存在大量的GOTO語句以及程序嵌套，必須對源代碼進行重構才能實現多線程并行處理。
3.求解速度較低
Nemoh目前僅提供直接求解法。對1000個單元左右的模型進行20個周期的計算耗時在15分鐘左右，常規的商業軟件對于同樣的問題求解時間在2~5分鐘左右。
4.后處理較弱
實際上目前版本的Nemoh是不存在常規意義上的后處理功能。用戶需要圍繞Nemoh進行大量的后續開發工作來實現諸如運動幅值響應算子（ ， RAO）、遠場波浪力（基于運動RAO和科欽函數）、近場波浪力（基于速度勢、運動RAO、速度勢梯度等）、模型顯示、波面模擬，乃至截面載荷計算和整體彎矩/剪力計算等功能。
相對而言，Nemoh已經將最復雜的問題（一階速度勢以及波浪載荷的計算）解決，這是它最大的優勢。用戶可以在此基礎上結合自身特定要求進行相關的開發活動。另外，由于速度向量文件已知，用戶可以通過重寫該文件來實現更廣泛的應用，諸如浮體的廣義自由度的應用分析等。
05
—
圍繞Nemoh進行的開發工作
: a -based of Nemoh
重新編寫了Nemoh的求解程序，將這部分嵌入中形成了[5]。目前的最新版本為1.2（2020.4）。
嵌入通過重寫的Nemoh格林函數部分。同時將和的科學計算庫進行整合。目前能夠實現的功能包括：附加質量、輻射阻尼、繞射波浪力以及入射波浪力的求解； API接口；對格林函數的編寫求解進行了優化；基于利用多線程加速求解；實現運動與自由波面的三維模擬等。
相比于Nemoh，的代碼可讀性和可維護性更好。將和通過f2py進行跨平臺編程，將耗費資源較高的求解內容通過計算，通過來實現前后處理，有更好的功能分工和跨平臺性能。
未來將會把開發重點放到程序求解效率的提高以及程序通用性。
代碼已開源，可在進行下載。
Wave , WARP[3]
Wave ，WARP，著作權登記號。是作者圍繞Nemoh進行開發的軟件，直白一點的講，是本人給Nemoh套了個殼并對原生Nemoh做了部分修改。
WARP目前有兩個版本：Demo版和開發版，Demo版用于原生Nemoh的簡單前處理與運行文件的編寫，開發版中對原生Nemoh進行了部分代碼優化，并添加一些后處理和數據計算功能。整合了以實現模型顯示與修復以及靜水力計算。
圖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 Mesh
用戶通過輸入參數來定義Nemoh.cal文件，程序通過調用Nemoh來實現波浪力與附加質量和輻射阻尼的計算。
圖5.3 WARP 1.23工況設置
Fig5.3 WARP 1.23 Case
2. RAO與遠場波浪力的求解
用戶輸入浮體質量信息后，程序調用Nemoh計算的數據來進行浮體RAO的計算。
程序調用Nemoh計算的科欽函數和計算完畢的RAO數據進行遠場波浪載荷的求解。作者對原生Nemoh進行了部分修改，現在可以通過WARP計算遠場法計算的艏搖二階定常波浪力。
圖5.4 WARP 1.23 RAO參數輸入界面
Fig5.4 WARP 1.23 RAO
3. 簡單的頻域運動分析
用戶可以通過定義關注點和波浪譜（目前僅支持譜）進行波頻頻域運動分析。對于頻域系泊分析，可以通過來實現。
4. 數據曲線顯示
WARP可以將主要的計算結果（附加質量、輻射阻尼、一階波浪力、遠場波浪力、重心位置運動RAO及相位）進行顯示。對模型能夠進行顯示查看。
圖5.5 WARP 1.23 顯示升沉附加質量和輻射阻尼曲線
Fig5.5 WARP 1.23, Added Mass and of Heave
圖5.6 WARP 1.23 顯示橫搖波浪力矩曲線及相位
Fig5.6 WARP 1.23, Roll Wave Force and Phase Angle Curve
圖5.7 WARP 1.23 顯示橫搖運動RAO幅值曲線及相位
Fig5.7 WARP 1.23, Roll 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的開發重點將集中在以下幾個方面：
對Nemoh源代碼進行重構，實現多線程并行計算；
對Nemoh源代碼進行重構，將運動RAO、遠場算法以及近場算法進行集成；
對后處理進行擴展，添加波面與運動顯示以及面元壓力分布顯示；
評估集成或內核轉為的可行性。
06
—
總結與展望
目前在國內船舶與海洋工程工業界，尚無任何一款國產的基于勢流理論的邊界元水動力求解程序得到廣泛應用和認可。長期以來，國內相關專業分析領域被WAMIT、AQWA、以及WADAM所占據。在當前嚴峻的形勢下，立足自身開展相關專業軟件的研究十分迫切。立足于開源軟件進行開發不失為一種快速解決目前難題的捷徑。
本文針對開源邊界元三維輻射繞射水動力計算程序Nemoh的基本理論，程序組成，使用方法進行介紹。對本人圍繞Nemoh進行的一些開發進行介紹，以供相關研究人員參考。
圍繞Nemoh可以開展的工作有很多，但核心還是在于提高計算效率和擴展計算功能兩大方面，這也將是未來相關開發工作的方向。
限于篇幅，本文不對Nemoh的計算精度進行對比分析什么情況下要計算節點坐標，如有條件將以單獨文章對這部分內容進行介紹。
歡迎相關領域研究人員和對本文所述開發工作感興趣的老板共同交流探討。
本號也將開設專欄對Nemoh的使用和在其基礎上的開發和應用進行介紹，歡迎關注。
參考文獻：
[1]Aurélien , Géard . and of the opend BEM NEMOH[C]. EWTEC 2015.
[2] Aurélien , Géard . The of “ and of the opend BEM NEMOH”.EWTEC 2015.
[3]WARP V1.00 軟件幫助[R].2018.6.
[4]Géard . Les de - et de de : Etude et par la des [D]PhD , école Supé de Mé de , 1987.
[5]: a -based of Nemoh,
Email：