封培元，馬寧，顧解忡

（1.中國船舶及海洋工程設(shè)計研究院，上海200011；2.上海交通大學(xué)船舶海洋與建筑工程學(xué)院，上海200240）

基于MFS無網(wǎng)格方法的船舶耐波性預(yù)報

封培元1，馬寧2，顧解忡2

（1.中國船舶及海洋工程設(shè)計研究院，上海200011；2.上海交通大學(xué)船舶海洋與建筑工程學(xué)院，上海200240）

文章提出了一種新的基于MFS法（Method of Fundamental Solutions，即基本解法）的求解船舶勢流問題的無網(wǎng)格數(shù)值方法。該方法采用勢流控制方程的基本解表達(dá)流場中的速度勢，可以得到任意階連續(xù)的速度勢導(dǎo)數(shù)；另外，通過引入流場虛邊界的概念避免了傳統(tǒng)邊界元方法中存在奇點的問題。相比傳統(tǒng)的邊界元方法，MFS法具有無奇性、不依賴網(wǎng)格、不涉及數(shù)值積分、數(shù)學(xué)原理簡單和算法易于編程實現(xiàn)的優(yōu)勢。文中對于MFS法在求解淺水二維剖面水動力系數(shù)中的應(yīng)用進(jìn)行了探索，提出了一套實用的數(shù)值算法，并進(jìn)一步結(jié)合頻域切片法計算了Wigley船體的垂蕩和縱搖響應(yīng)，通過和頻域面元法及試驗結(jié)果的對比驗證了方法的正確性和實用性；最后應(yīng)用開發(fā)的程序分析了不同水深和航速對于耐波性的影響。該研究開創(chuàng)性地將MFS方法應(yīng)用到了船舶水動力研究領(lǐng)域，為之后更深入廣泛地應(yīng)用這一特性優(yōu)良的數(shù)值方法提供了有益的參考。

MFS；無網(wǎng)格；切片法；勢流；耐波性

0 引言

忽略流體粘性的勢流模型在船舶水動力學(xué)研究中有著廣泛的應(yīng)用，其核心是求解不同邊界條件下Laplace方程的解。由于船舶水動力研究中自由面、物面和無窮遠(yuǎn)輻射邊界條件的復(fù)雜性，使得尋求解析解變得困難重重，一般情況下需要借助數(shù)值方法來進(jìn)行求解。目前，以Hess和Smith[1]在1964年提出的面元法應(yīng)用最為廣泛，其基本原理是將物面離散為源強為常數(shù)的面元，再選取滿足自由面和無窮遠(yuǎn)輻射條件的復(fù)雜格林函數(shù)根據(jù)格林公式進(jìn)行數(shù)值積分得到物面上的源強分布，最后計算出物面壓力分布從而得到附加質(zhì)量等水動力系數(shù)。該方法的優(yōu)點是魯棒性強，可以針對任意形狀的物體進(jìn)行計算；但也存在一些限制，比如：（1）三維計算中涉及繁瑣的物面網(wǎng)格生成；（2）假定物面速度勢為常數(shù)后速度勢的導(dǎo)數(shù)不連續(xù)，因而在計算二階力時容易導(dǎo)致發(fā)散；（3）面元分割改變了物面幾何形狀，因此物面邊界條件只能得到近似滿足；（4）對格林公式進(jìn)行積分時存在奇點問題，需要特別處理，因而增加了編程的難度。

相比而言，若采用無網(wǎng)格的數(shù)值方法則不涉及對于流域的離散，可以有效化解以上面元法中存在的問題，因而近年來受到許多學(xué)者的關(guān)注。目前，常用的無網(wǎng)格法有十多種，它們基本特征是將連續(xù)體離散為有限數(shù)量的質(zhì)點集合，提供相應(yīng)的近似函數(shù)，從而避免了繁瑣的單元網(wǎng)格劃分[2]。MFS法就是其中一種較新發(fā)展出的無網(wǎng)格的數(shù)值方法，該方法對于存在基本解的邊值問題特別有效，并且物體幾何形狀越是復(fù)雜就越能體現(xiàn)出MFS法的優(yōu)越性[3]。MFS法用流域控制方程基本解的線性組合作為流場的近似解，可以得到任意階連續(xù)的導(dǎo)數(shù)，并且通過引入虛邊界的概念完全避免了原先邊界積分中存在奇點的問題。另外，研究指出MFS法具有收斂速度快和精度高的優(yōu)點[4-5]；對于存在自由面或者移動邊界的問題，MFS更能發(fā)揮自身的特點[6]，已被成功應(yīng)用于三維造波模擬[7]以及斜向行進(jìn)波遭遇障礙物后的演變模擬[8]等研究中，取得了很好的效果。

本次研究中進(jìn)一步將MFS法引入船舶耐波性領(lǐng)域：利用MFS法計算有限水深二維剖面的水動力系數(shù)，并以此為基礎(chǔ)，結(jié)合頻域切片法理論預(yù)報了船舶在淺水中的耐波特性。以一艘Wigley船為例，證明了所提出算法的正確性和實用性。本文首先詳細(xì)論述了MFS法的原理和計算方法；其次給出了MFS法求解有限水深二維剖面水動力系數(shù)的具體流程；接著簡要回顧了切片理論的相關(guān)計算公式；然后給出了Wigley船在不同水深和航速下的垂蕩和縱搖響應(yīng)計算結(jié)果，并對結(jié)果進(jìn)行了驗證和分析；最后總結(jié)了本文的研究并對之后的研究工作進(jìn)行了展望。

1 MFS法基本原理

勢流問題由Laplace方程控制：

其中：Φ為流場的速度勢。MFS法的核心思想是用基本解的線性組合擬合流場中的物理量，再讓邊界上的物理量滿足給定的邊界條件從而求出擬合物理量時未知的形狀系數(shù)。因此，將速度勢表達(dá)為如下形式：

其中：Pi為屬于流域Ω的任意場點；c=［c1,c2,…,cN］T為擬合速度勢的未知形狀系數(shù)向量；S為包含源點Sj坐標(biāo)的N維向量。

KD為Laplace方程的基本解，對于二維和三維問題分別有：

用MFS法進(jìn)行計算時，首先要分別布置好場點和源點。場點為配置點，位于實際的流域邊界上，需要滿足對應(yīng)的邊界條件；源點則位于真實流場外的虛邊界上（如圖1所示），這樣可以避免場點和源點重合時產(chǎn)生的奇點問題。需要指出的是，虛邊界的選?。丛袋c的布置）對于MFS法的精度至關(guān)重要。若虛邊界十分靠近真實邊界，則形狀系數(shù)求解矩陣的對角線元素值將趨近于無窮大，影響求解精度；而當(dāng)虛邊界遠(yuǎn)離真實邊界時又將導(dǎo)致形狀系數(shù)求解矩陣的條件數(shù)快速增長，使矩陣變得病態(tài)[9]。文獻(xiàn)[10]中就源點的最優(yōu)布置方法進(jìn)行了探討并指出虛邊界的選取問題目前仍是制約MFS法發(fā)展的關(guān)鍵因素。本次研究通過數(shù)值試驗進(jìn)行不斷嘗試后發(fā)現(xiàn)：可以取物面場點間的平均間距作為虛邊界位移量，由此根據(jù)每個場點的位置沿真實邊界的外法向進(jìn)行平移得到源點的布置位置。完成場點和源點的布置后，只要根據(jù)（2）式表達(dá)速度勢并讓每個場點處的速度勢均滿足流域邊界條件即可。一般的邊界條件分為Dirichlet型和Neumann型兩種，因此有：

其中：fD和fN為Dirichlet和Neumann邊界條件的一般表達(dá)；χD和χN為兩種邊界條件的組合系數(shù)。最后，對每個場點建立方程，得到如下形式的方程組：

求解后得到形狀系數(shù)向量c，再由（2）式可求得流場內(nèi)任意點的速度勢及其任意階導(dǎo)數(shù)。

圖1 MFS法求解二維剖面水動力系數(shù)Fig.1 Solving 2D section hydrodynamic coefficients based on MFS

2 有限水深二維剖面水動力系數(shù)計算

二維剖面水動力的求解包括橫蕩、垂蕩和橫搖三個模態(tài)，分別用下標(biāo)2，3，4表示。有限水深情況下，計算各個模態(tài)速度勢的定解問題（如圖1所示）為：

其中：φj為規(guī)范化速度勢，j=2，3，4；nj為物面處單位外法向量的j模態(tài)分量；ω0為波浪圓頻率；k為波數(shù)與ω0間滿足色散關(guān)系：

其中：H為水深。對于無窮遠(yuǎn)邊界，只能在距浮體有限遠(yuǎn)處進(jìn)行截斷，并用兩根垂線代替。這一處理無窮遠(yuǎn)邊界的方法在采用簡單格林函數(shù)的面元法中有過成功的應(yīng)用[11-12]。

將（7）式至（10）式按（4）式的形式進(jìn)行表達(dá)，并對每個場點建立相應(yīng)的方程后進(jìn)行求解，得到形狀系數(shù)向量后可由（2）式得到物面上規(guī)范速度勢的數(shù)值。由此可按（12）式計算出二維剖面的附加質(zhì)量和阻尼系數(shù)：

3 頻域切片理論

本文將二維剖面水動力計算與Salvesen等[13]提出的頻域切片理論相結(jié)合，用于船舶垂蕩和縱搖響應(yīng)的預(yù)報。根據(jù)文獻(xiàn)[13]建立坐標(biāo)系，得到船體垂蕩和縱搖的運動方程為：

其中：η3和η5分別為垂蕩和縱搖；M為船舶排水量；I5為船體縱搖慣量；ω為遭遇頻率，與航速U和浪向β間滿足關(guān)系：

Aij，Bij，Cij，F(xiàn)3和F5分別為對各個二維船體剖面積分后全船的附加質(zhì)量、阻尼系數(shù)、恢復(fù)系數(shù)，一階垂蕩波浪力和一階縱搖波浪力矩。詳細(xì)的計算公式可參考文獻(xiàn)[13]。

4 Wigley船計算結(jié)果

通過對一艘Wigley船的垂蕩和縱搖響應(yīng)計算結(jié)果驗證MFS法結(jié)合切片理論預(yù)報船舶耐波性的有效性和實用性。Wigley船為數(shù)值船型，主要參數(shù)如表1所列。計算時先將船體沿垂線間長等間距分解為21個剖面，分別應(yīng)用MFS法求解附加質(zhì)量和阻尼系數(shù)，最后再通過切片理論積分得到全船的水動力系數(shù)并根據(jù)（13）式和（14）式求解運動響應(yīng)幅值算子（RAO）。

表1 Wigley船主要參數(shù)Tab.1 Main particulars of Wigley Hull

首先計算了傅汝德數(shù)Fn為0.3的情況，并與Delft水池試驗結(jié)果[14]和法國船級社三維頻域面元法軟件HydroSTAR的計算結(jié)果進(jìn)行了對比。計算中，取水深為5倍的吃水；自由面寬度為6倍的船寬；物面上按0.02 m等間距布置場點，同時以0.02 m作為虛邊界位移量布置源點。單側(cè)自由面、水底和左/右無窮遠(yuǎn)邊界分別布置60、60和50個場點。計算得到的垂蕩和縱搖RAO如圖2和3所示。

圖2 垂蕩響應(yīng)計算結(jié)果驗證Fig.2 Validation of heave response

圖3 縱搖響應(yīng)計算結(jié)果驗證Fig.3 Validation of pitch response

從圖中的對比可以發(fā)現(xiàn)，本方法的計算結(jié)果與試驗符合較好，與頻域面元法間存在一定差異，原因可能是由于頻域面元法考慮了船體的三維特性并且在處理航速時與切片法有所不同造成的。

圖4 水深對于垂蕩響應(yīng)的影響Fig.4 Effect of water depth on heave response

圖5 水深對于縱搖響應(yīng)的影響Fig.5 Effect of water depth on pitch response

圖6 航速對于垂蕩響應(yīng)的影響Fig.6 Effect of forward speed on heave response

圖7 航速對于縱搖響應(yīng)的影響Fig.7 Effect of forward speed on pitch response

利用這套算法進(jìn)一步研究了水深與航速對于該Wigley船耐波性的影響。先是計算了水深為3d，4d和5d水深下的垂蕩和縱搖RAO曲線，d為船舶吃水。計算中根據(jù)水深的變化等量調(diào)整了左/右無窮遠(yuǎn)邊界上的場點個數(shù)，分別為30，40和50，其他計算參數(shù)不變。從圖4和5所示的結(jié)果中可以看出，水深對于垂蕩和縱搖的峰值和長波區(qū)域的響應(yīng)有一定影響，隨著水深的增加響應(yīng)變得更劇烈，三維面元法的計算結(jié)果也證實了這一變化規(guī)律。航速對于Wigley船垂蕩和縱搖的影響如圖6和7所示，分別計算了Fn=0.0，0.1，0.2和0.3四種情況下的RAO曲線。從圖中可見，航速對于該船的影響十分顯著，垂蕩和縱搖的峰值隨著航速的上升而顯著增加。

5 結(jié)語

本文采用了一種新的無網(wǎng)格方法MFS研究了有限水深中二維剖面的水動力系數(shù)求解問題；提出了一套實用的數(shù)值方法，在求解過程中體現(xiàn)出了MFS法無奇性、不依賴網(wǎng)格、不涉及數(shù)值積分、數(shù)學(xué)原理簡單和算法易于編程實現(xiàn)的優(yōu)勢。另外，將二維剖面水動力系數(shù)求解與頻域切片法相結(jié)合，用于船舶耐波性的預(yù)報。以一艘Wigley船型為例，驗證了方法的正確性和可行性，并借助該工具研究了水深和航速對于Wigley船耐波性的影響，結(jié)果表明這兩個參數(shù)對于船舶耐波性有著很大的影響。本文研究為更深入地研究MFS法打下了基礎(chǔ)，之后擬開展的工作包括虛邊界的最優(yōu)確定方法，MFS法在三維水動力計算中的應(yīng)用以及MFS法在時域耐波性研究中的應(yīng)用。這些研究將更能體現(xiàn)MFS法自身的優(yōu)勢，具有很好的應(yīng)用前景。

[1]Hess J L,Smith A M O.Calculation of nonlifting potential flow about arbitrary three-dimensional bodies[J].J Ship Res, 1964,8(3)：22-44.

[2]王宇新,孫明,張建臣.無網(wǎng)格MPM法三維前處理系統(tǒng)設(shè)計[J].計算力學(xué)學(xué)報,2008,3：392-396. Wang Yuxin,Sun Ming,Zhang Jianchen.Implementation of 3D preprocessing program for material point method[J].Chinese Journal of Computational Mechanics,2008,3：392-396.

[3]Fairwather G,Karageorghis A.The method of fundamental solutions for elliptic boundary value problems[J].Advances in Computational Mathematics,1998,9：69-95.

[4]Peter M,Youssef F R.Convergence and stability of the method of meshless fundamental solutions using an array of randomly distributed sources[J].Eng Anal Bound Elem,2004,28：143-153.

[5]Li X.On convergence of the method of fundamental solutions for solving the Dirichlet problems of Poisson’s equation[J]. Advances in Computational Mathematics,2005,23：265-277.

[6]Sarler B.Solution of potential flow problems by the modified method of fundamental solutions：formulations with the single layer and the double layer fundamental solutions[J].Eng Anal Bound Elem,2009,33(12)：1374-82.

[7]Wu N J,Tsay T K.Applicability of the method of fundamental solutions to 3-D wave-body interaction with fully nonlinear free surface[J].J Eng Math,2009,63：61-78.

[8]Chen K H,Lu M C,Hsu H M.Regularized meshless method analysis of the problem of obliquely incident water wave[J]. Eng Anal Bound Elem,2011,35：355-362.

[9]Wang H,Qin Q H.Some problems with the method of fundamental solution using radial basis functions[J].Acta Mechanica Solida Sinica,2007,20(1)：21-29.

[10]Karageorghis A.A practical algorithm for determining the optimal pseudo-boundary in the method of fundamental solutions [J].Adv.Appl.Math.Mech.,2009,1(4)：510-528.

[11]Hasanadil M,Duan W Y,Wang Y.Hydrodynamic coefficients of ships with forward speed[J].Journal of Ship Mechanics, 2004,8(3)：46-54

[12]Milgram J H.Strip theory for underwater vehicles in water of finite depth[J].J Eng Math,2007,58：31-50.

[13]Salvesen N,Tuck E O,Faltinsen O M.Ship motions and sea loads[J].Trans.SNAME,1970,78：250-287.

[14]Journée J.Experiments and calculations on 4 Wigley hull forms in head waves[R].Report of Delft University of Technology,Ship Hydromechanics Laboratory,The Netherlands,1992.

Ship seakeeping prediction based on MFS meshfree scheme

FENG Pei-yuan1,MA Ning2,GU Xie-chong2
(1.Marine Design&Research Institute of China,Shanghai 200011,China;2.School of Naval Architecture, Ocean and Civil Engineering,Shanghai Jiaotong University,Shanghai 200240,China)

This paper presents a novel meshfree numerical scheme based on MFS(Method of Fundamental Solutions)for solving ship potential flow problems.This method adopts the fundamental solution of the potential flow control equation to represent the velocity potentials and can obtain any degree of continous velocity potential derivative.Moreover,through the introduction of pseudo boundary,the singularity in conventional Boundary Element Method(BEM)is circumvented.Compared with BEM,MFS has the merits of being desingularized,meshfree,integration-free,mathematically simple and easy for programming.This study explores the application of MFS in computing the hydrodynamic coefficients of two dimensional section in water of finite depth and proposes a practical numerical scheme;furthermore,MFS is combined with frequency domain strip theory to predict the heave and pitch response of a Wigley hull.The validity of the proposed numerical scheme is proved through comparison with panel method and experimental result.Finally,the influence of water depth and forward speed on seakeeping is studied using the current scheme.This study innovatively applies MFS into the filed of ship hydrodynamics and can provide some guidance for the future development of MFS.

MFS;meshfree;strip theory;potential flow;seakeeping

U661.1

A

10.3969/j.issn.1007-7294.2017.04.002

1007-7294（2017）04-0390-06

2016-10-07

國家自然科學(xué)基金（51579144）；工信部高技術(shù)船舶科研項目（K24352）

封培元（1987-），男，博士，E-mail：pyfeng23@163.com；馬寧（1961-），男，教授，長江學(xué)者。