怎麼樣編程計算浮體的附加質量和附加阻尼？

01-16

謝邀~

我先默認這裡的「浮體的附加質量和附加阻尼」是指在linear potential theory框架下，由airy wave引起的wave radiation force。

其實我覺得這句話應該已經把怎麼計算說的差不多了=。=

不知道提問這個問題的人是出於什麼目的。

1，不管是做industry project還是自己的research topic, 如果要解決的終極問題不是研究如何更好地計算附加質量和附加阻尼，而是在所做的topic里需要浮體的added mass and damping coefficients作為必需的輸入的話，在手上有商業軟體的前提下，用商業軟體計算。因為linear potential theory是很成熟的理論，也是很成熟的計算方法，無論出於學術目的，還是工業項目上的分析，用成熟商業軟體是最快捷並且精確的方法。自主編程基本上無學術價值可言。

2，手上沒有商業軟體，或者這部分代碼是自己一個更大系統的in-house code裡面的一部分的話，我推薦參考這本書裡面的內容。 船舶在波浪上的運動理論。劉應中，繆國平。這本書對於最基本的potential theory, boundary element method, perturbation theory, radiation and diffraction problem, wave force都有比較明確的講解。

求解added mass and damping coefficients, i.e, wave radiation force，的本質還是要求解浮體運動引起的radiation potential。首先要明確governing equation （Laplace equation)+boundary condition (free surface, body surface, bottom, infinite boundary). 其次，了解求解lapalce equation的數值解法，應該都是boundary element method，當然裡面會涉及很多很具體的網格劃分，單元存儲，矩陣解法，源匯布置等（computational我基本只懂點概念，沒自己上手寫過代碼=。=）。potential求解之後，基本上用pressure integral就可以得到added mass damping coefficients了。

PS:想到如果是想計算multi-body interaction下的added mass damping的話，我不確定多少商業軟體能算以及精確度如何（我自己沒做過這方面的工作，只能確定應該是有商業軟體按能算的），但是問題本質沒有變，想算好added mass damping，就是要把流場的radiation potential算準確。

PPS: master的時候上船舶在波浪上的運動理論這門課，最後的課程project就是自己編程算二維的hydrodynamic coefficients, i.e. added mass, potential damping, first-order wave excitation force, 浮體形狀是free surface上是一個半圓，水下部分還有一個圓。自己從畫網格，布置源匯，解potential開始算。理論感覺很簡單吧，物面形狀也很規則吧，網格也很好畫容易存儲吧，二維potential還是比三維簡單一點的吧。

最後結果就是debug到converge，並且形狀看起來像那麼回事已經耗盡心力了，最後同學之間的結果，只要是自己獨立寫的code不是copy過來的，基本上都不太一樣o(╯□╰)o (當然，課程project花費的時間和精力肯定不能跟有目的性的research code相提並論）

所以，恩，成熟的商業軟體受眾廣泛並且賣價那麼高是有道理的，必定核心solver都是做了很多專業的演算法優化，stability test等等的。

以上。僅適用於floating structure under airy wave.

原諒中英文混合的專業名詞描述，我基本上是第一反應是中文的名詞就用了中文，第一反應是英文的名詞就用了英文。。。（其實應該考慮一下以後專業問題都用英文答，正好鍛煉一下現在的渣英文能力。。。以上答案真實的反應了我現在思考專業問題時候大腦的語言反饋，結果就是跟老闆交流的時候一腦子的名詞，短語，構成整句的時候就開始語言障礙~~o(&>_&<)o ~~）

首先理解什麼是附加質量和阻尼係數，它們是線性系統下才有的概念（控制方程及邊界條件都為線性方程，誠如樓上所言），而且在勢流理論的框架下，為波物作用(wave-body interaction)時空解耦 (decoupling)的產物,這裡不做詳細展開。

具體來說，附加質量表徵慣性力大小，這種力與物體加速度成正比，在求解運動方程時，根據牛二，我們直接將它跟質量放在一起，因此我們形象叫它~附加質量(added mass)；對於阻尼係數，這裡的阻尼不是粘性阻尼，因為物理模型我們已經假設無旋無粘不可壓，這裡的阻尼指的是由於物體在含有自由表面的流體中運動，勢必會興波，而興起波浪必然會帶走能量，這個能量就是物體運動做的功，體現為物體受到的阻力，大小與物體運動的速度相關，嚴格來說，這個阻尼係數應該叫 wave-radiated damping, 也就是興波阻力，這就是為什麼勢流工具可以用來進行船型優化，使(興波)阻力最小. 至此，我們描述了物體在流體中運動的兩種力，一種與物體加速正比，一種與速度正比，此外再加上與物體位移正比的恢復力，一旦物體運動確定，我們就可以方便的求解運動方程(顯然，這是線性系統，顯然運動與受力不耦合，顯然實際情況中我們不太可能把物體受力就分離出這幾個還有許多非線性項, 顯然這兩個概念在頻域中尤為合適)

現在再說計算機編程實現，在勢流理論框架下，對於這種線性問題，最方便的求解方法莫過於面元法 (panel method) 它本質上屬於邊界元方法(boundary element method)，最先將它應用於船舶水動力學的是Hess,Smith 兩位大神，他們事實上是搞飛機的，所以這種方法又叫Hess-Smith 方法，手機做答不便貼文獻，中文文獻可參見《船舶在波浪中運動的勢流理論》哈工程戴遺山段文洋那本紅寶書，上面第二章給出了詳細的推倒，只要將其翻譯為程序語言即可，另外交大劉應中編的《船舶計算流體力學》也有詳細介紹，以上兩本都是詳細到快要寫出演算法的程度。外文書籍有（practical ship hydrodynamics）現在已出第二版。如果是學生，以後做這個方向，建議自己動手編一下，編完程序才說明真正掌握了其中概念，不過，It is a journey of study rather than research. 如果是科研需要，講求短平快，許多商業軟體可以解決，網上也有好多開源代碼，Google一下會搜出很多。

春節趕飛機路上作答，倉促，希望有幫助。

@Emily Ma 姐姐已經簡單說了一下基本思路，我就補（zhao）充（chao）一(bi)下(ji)具體細節吧（默默翻出本科時蹭段老師的水動力學時做的作業...... )

知乎公式和圖片環境各種不會用orz，請大家原諒那碩大無比的公式....

以下討論限定在零航速情形.

我們知道，在勢流理論中，浮體周圍的波浪速度勢可以分解為三部分：

$varphi(X)e^{-iomega t}=[(varphi_l+varphi_d)+sum_6^{j=1}varphi_rj x_j]e^{-jomega t}$

式中， $varphi_l$ 為單位波幅的一階入射勢, $varphi _d$ 為對應的繞射勢， $varphi _{rj}$ 為第j模態單位運動引起的輻射勢.
這三種速度勢分別對應了入射波產生的Froude-Krylov力，繞射波產生的繞射力和浮體自身運動引起的輻射波產生的輻射力.
而之所以引入附加質量和附加阻尼這樣的水動力係數正是為了簡單愉快地刻畫輻射力 (^_^).

浮體作第k模態單位運動引起的輻射波產生的第j模態輻射力可以表達為:

$F_{rjk}= -iomega ho int_{S_0}varphi_{rk}n_jdS$

將輻射速度勢展開為實部和虛部並代入上式，即可得到附加質量和波浪阻尼的表達式如下:

其中 $A_{jk}$ 與 $B_{jk}$ 即為附加質量係數與附加阻尼係數. 可以看出，求解附加質量係數與附加阻尼係數關鍵在於求得流體速度勢。

在理想流體假設下，令流場速度勢中與時間無關的項為 $varphi$ ，那麼流場中波浪的繞射/輻射效應就可以用以Laplace方程進行描述, 其邊界條件為:

（1）零航速線性自由麵條件

（2）物面邊界條件

（3）海底不可穿透條件

（4）遠場條件: 無窮遠總波浪擾動趨於0

求解在上述控制方程下的流體速度勢通常採用格林函數法。引入在頻域上周期格林函數 $G(X,xi ,omega )$ (格林函數解析表達式實在太長orz我就不寫下來了), 根據Green定理波浪的速度勢可以通過第二類Fredholm邊界積分表示:

在平均濕表面上引入分布源函數後，流體速度勢可以進一步表示為:

其中平均濕表面上的源強函數 $sigma (xi )$ 根據前面給出的物面邊界條件來確定:

求解上述方程通常採用Hess-Smith方法，即將平均濕表面劃分為許多個四邊形或三角形單元，並用單元內一點的源強和速度勢代表整個單元面. 這樣，上述積分方程就可以轉化為如下的離散形式：

式中， $N_p$ 是在平均濕表面上劃分的面元總數， $Delta S_m$ 為第m個單元的面積，而 $xi_m$ , $X_k$ 分別為第m個單元和第k個單元幾何中心的坐標。

這樣，原積分方程就轉化為了一個 $N_p$ 階線性方程組, 求解該方程組即可得到各個單元速度勢的分布, 再根據水動力係數和速度勢的關係對各個單元遍歷求和, 即可得到所需的附加質量和附加阻尼了。

Reference

[1] 戴遺山, 段文洋. 船舶在波浪中運動的勢流理論. 國防工業出版社, 2008.

[2] Aqwa Theory Manual

按照 @穆言的建議，我把評論里自己的疑問寫成一個回答貼在這裡。實際上算是開一個有關Haskind Relations的新問題了。

以下是書中有關Haskind Relations的論述

--------------------------------------------------------------------

在一階線性假設下，入射勢 $phi_I$ 和繞射勢 $phi_D$ 對於浮體產生的FK力和繞射力之和（激發力）為

$old X_i(omega)=-iomega hointint_{S_B}(phi_I+phi_D)n_idS$

由於輻射勢 $phi_i$ 的在濕表面 $S_B$ 法嚮導數

$frac{partialphi_i}{partial n}=iomega n_i$

所以激發力又可以寫成

$old X_i(omega)=- hointint_{S_B}(phi_I+phi_D)frac{partialphi_i}{partial n}dS$

格林公式指出，對於任意兩個勢函數 $phi_j,phi_k$ 可寫

$intintint_V(phi_j riangledown^2phi_k-phi_k riangledown^2phi_j)dV=intint_S(phi_jfrac{partialphi_k}{partial n}-phi_kfrac{partialphi_j}{partial n})dS$

積分區域為整個流體邊界。

所以在繞射勢 $phi_D$ 和輻射勢 $phi_i$ 均滿足拉普拉斯方程的條件下

$intint_{S_0}phi_Dfrac{partialphi_i}{partial n}dS=intint_{S_0}phi_ifrac{partialphi_D}{partial n}dS$

註：此時積分面已經過簡化只剩浮體濕表面。

又由於入射勢和繞射勢在濕表面不可穿透的特性

$frac{partialphi_I}{partial n}=-frac{partialphi_D}{partial n}$

所以激發力可以改寫為

$old X_i(omega)=- hointint_{S_B}(phi_Ifrac{partialphi_i}{partial n}-phi_ifrac{partialphi_I}{partial n})dS$

此時式子中已經沒有難以求解的繞射勢，其被輻射勢和入射勢的表達所替代。而這個式子就是所謂的Haskind Relations。

---------------------------------------------------------------------------

我的問題就是，前述內容在使用格林公式時，代入的是繞射勢和輻射勢。如果我同樣的把入射勢和輻射勢代入，會不會同樣得到

$intint_{S_0}phi_Ifrac{partialphi_i}{partial n}dS=intint_{S_0}phi_ifrac{partialphi_I}{partial n}dS$

使得最後激發力積分項內部減號兩邊相等，激發力最終等於0？

實際操作中仔細處理格林函數的主值積分，要參考Abramowitz and Stugun的數學手冊第六章指數積分內容