基于高階半解析移動脈動源法的船舶波浪增阻數(shù)值計算研究

摘 要： 在三維勢流理論框架下，以移動脈動源格林函數(shù)為核函數(shù)構建邊界積分方程，采用九節(jié)點高階曲面單元進行離散求解速度勢，建立了波浪中航行船舶運動響應和波浪增阻的計算模型.為了避免數(shù)值計算的不穩(wěn)定，積分方程中的影響系數(shù)采用了半解析公式進行計算.求解波浪增阻時，采用將輻射能量法、Salvensen法和短波半經(jīng)驗公式相結合的混合法，以便能夠在全頻率范圍內(nèi)得到令人滿意的結果.在此基礎上，自主開發(fā)數(shù)值程序?qū)Σ煌驮诓煌l率、航速下的運動響應和波浪增阻進行預報.通過將預報值與試驗及其它數(shù)值方法的結果進行對比發(fā)現(xiàn)，文中方法的數(shù)值結果與試驗值吻合良好，對于不同船型和工況均有著良好的適用性；相比傳統(tǒng)基于數(shù)值求積公式的移動脈動源法，高階半解析移動脈動源法計算更為穩(wěn)定（尤其對于存在外飄的船型），且由于采用了高階格式，在共振頻率附近的預報精度比常值元法更高.

關鍵詞： 波浪增阻；移動脈動源；高階；半解析；混合法

中圖分類號：U661.32"" 文獻標志碼：A"""" 文章編號：1673-4807（2024）02-015-08

Numerical study on added resistance of ships advancing in waves byhigher-order semi-analytical translating-pulsating source method

Abstract：Within the framework of three-dimensional potential flow theory， a computational model based on translating-pulsating source Green′s function and 9-node higher-order discretization is established for the motion responses and added resistance of ships advancing in waves. In order to avoid the numerical instability， semi-analytical quadrature formula is derived and employed to evaluate the influence coefficients of the boundary integral equation. And the added resistance is calculated by a combined formula of radiation energy method， Salvensen′s method and semi-empirical method in short waves. Based on the above method， a computer code is developed to investigate the motions and added resistance of various ships at different frequencies and forward speeds. By comparing the predicted results with those of experiment and other numerical methods， it is found that the higher-order semi-analytical translating-pulsating source method can work well for different ship types and working conditions. Besides， the present method appears to be more stable than the translating-pulsating source method based on traditional numerical quadrature formula （especially for the non-wall-sided ship） and more accurate than the constant panel method.

Key words：added resistance， translating-pulsating source， higher order， semi-analytical quadrature， combined method

船舶在海上航行時，不可避免地會遭遇波浪.在波浪的作用下，船舶除了會產(chǎn)生各種搖蕩運動之外，受到的阻力相比于靜水也會有所增加，這無疑會大大增加船舶的燃料消耗、影響其保持航速的能力等.因此，尋求合理而實用的數(shù)值方法對航行船舶的波浪增阻進行預報研究，對船舶的設計、使用等具有重要意義.勢流方法因其高效并能保有足夠的精度，在波浪增阻問題中應用廣泛.最簡單的勢流方法是切片法，早期大多數(shù)波浪增阻計算方法都是基于此來表達的[1-2].但是該方法在計算中采用的高頻、低速、船體細長等假定，限制了其進一步的發(fā)展與應用.近年來出現(xiàn)了大量關于三維勢流方法的研究工作，這些方法根據(jù)構建積分方程時采用的核函數(shù)的不同，通常可以分為Rankine源法和自由面格林函數(shù)法兩類.Rankine源法的優(yōu)點在于可以適用于不同形式的邊界條件，考慮非線性、淺水效應等，但在有航速水動力分析中，它存在著離散量大、難以處理Brard數(shù)較小時的輻射條件等缺陷限.

相對而言，自由面格林函數(shù)法由于選擇的格林函數(shù)自動滿足自由面和輻射條件，只需在物面布源，離散量較小.利用該方法研究船-波相互作用問題最為關鍵的一點在于格林函數(shù)的計算.針對無航速自由面格林函數(shù)（又稱脈動源格林函數(shù)），文獻[3]做了大量的研究工作，提出了快速高效的數(shù)值算法，基于此開發(fā)了比較實用的無航速波浪與結構物相互作用問題的求解軟件；在此基礎上，文獻[4-5]基于脈動源格林函數(shù)，通過航速修正考慮航速效應，分別采用不同的遠場公式對船舶的波浪增阻和橫向二階力進行了計算.然而，這類方法只適用于航速較低的工況，對中高航速并不能得到令人滿意的數(shù)值結果.為了充分考慮航速效應，應該采用有航速自由面格林函數(shù)（又稱移動脈動源格林函數(shù)）.近年來，在有航速自由面格林函數(shù)及其面積分數(shù)值計算方法方面的研究工作[6-10]，在一定程度上促進了移動脈動源法的發(fā)展.但從實際工程意義角度來說，該方法距離真正應用于有航速船舶水動力分析仍然任重而道遠.尤其是對存在外飄的中高速船，格林函數(shù)在近自由面附近的劇烈振蕩特性[11]，使得移動脈動源法很難獲得穩(wěn)定的數(shù)值結果.鑒于此，目前很少有學者基于該方法對船舶波浪增阻問題進行研究.文獻[12]曾嘗試基于移動脈動源法的水動力結果，分別采用輻射能量法和Salvesen法[13]對由船舶搖蕩運動以及繞射作用引起的波浪增阻進行了計算.但受限于格林函數(shù)積分在自由面附近的不穩(wěn)定性，數(shù)值格式采用的是基于常值單元離散的低階格式，這是少有的基于移動脈動源法對波浪增阻的研究.

文中針對波浪中航行的船舶，基于移動脈動源格林函數(shù)，建立了一種計算波浪增阻的高階半解析方法.該方法利用九節(jié)點高階曲面單元離散以移動脈動源格林函數(shù)為核函數(shù)的邊界積分方程，分別采用解析公式和Gauss-Legendre求積公式計算影響系數(shù)沿水平和垂直方向的積分，以避免格林函數(shù)沿水平方向的劇烈振蕩所引起的數(shù)值計算不穩(wěn)定.求解波浪增阻時，為了在全頻率范圍內(nèi)得到令人滿意的結果，采用了將輻射能量法（用于計算輻射增阻）、Salvensen法（用于計算繞射增阻）和短波半經(jīng)驗公式相結合的混合法.在此基礎上，采用Fortran自主開發(fā)數(shù)值程序，對不同船型（包括數(shù)學船型modified Wigley和集裝箱船S175）在規(guī)則波中以不同航速航行時的運動和波浪增阻進行數(shù)值計算，并將結果與試驗和其它數(shù)值方法的結果進行對比分析.

1 求解速度勢的高階半解析移動脈動源法

1.1 速度勢與邊界積分方程

在勢流理論框架下，波浪中航行船舶的水動力力問題可以通過引入速度勢函數(shù)來研究.為了便于求解速度勢，建立以船舶平均速度U運動的參考坐標系o-xyz（圖1）.

該坐標系的原點o位于船舶重心正上方的平均靜水面，x軸與船舶前進方向相同，z軸豎直向上.當入射波是微幅波，船舶搖蕩運動幅值不大時，流場的非定常擾動勢可以分解為：

式中：px，y，z代表場點；qξ，η，ζ代表源點；σp和σq分別為p和q點處的源強；SB和CW分別船體的平均濕表面和水線；n1為單位外法線矢量n沿x軸的分量；G（p，q）為滿足自由面條件和輻射條件的移動脈動源格林函數(shù)，Havelock型單積分為[15]：

文中將移動脈動源格林函數(shù)Havelock型表達式中的前兩項記為GR（p，q），最后一項記為GF（p，q）.從表達式可以看出，GR只取決于場、源點的相對位置，而GF則還與航速等物理量相關.圖2為GR和GF沿水平方向的典型空間分布.可以看出，GR是一個平滑變化的函數(shù)，而航速相關部分GF沿水平方向則呈現(xiàn)出劇烈的振蕩特性.

1.2 高階面元離散及影響系數(shù)計算

邊界積分方程式（3）需要采用面元對邊界進行離散轉化為代數(shù)方程組后，才可以進行數(shù)值求解.最常用的離散面元是三角形或者四邊形常值單元.在常值元中，各未知物理量被假定為常值，因而比較簡單易于實現(xiàn).但是它存在離散量大、單元間變量不連續(xù)等固有缺陷.為了克服這些問題，文中采用九節(jié)點二次曲面單元對船體表面進行離散（圖3中的點代表控制點）.該單元中任意一點處的變量值可以利用控制點處的值插值獲得[16].

理論上，標準正方形單元上的積分可以直接采用Gauss-Legendre求積公式進行計算.然而，格林函數(shù)航速相關部分GF沿水平方向的劇烈振蕩特性（圖2）使得很難采用簡單的數(shù)值積分方法獲得關于其的穩(wěn)定積分結果[8].為了克服該問題，文中在計算與GF法向?qū)?shù)相關的面積分時，首先將圖3中的九節(jié)點單元分割為4個四節(jié)點單元（圖4）.然后采用不同的積分策略對四節(jié)點單元上的積分進行計算：垂向仍采用數(shù)值積分方法，而水平方向則推導了解析積分公式進行計算，以避免GF沿水平方向劇烈振蕩所導致的積分計算不穩(wěn)定問題.具體計算為：

系數(shù)的表達式見文獻[10].

圖5為采用上述半解析方法計算單位正方形（其形心位于（0， 0， -0.51）、法向沿y軸）上GF偏導數(shù)面積分的結果與傳統(tǒng)Gauss-Legendre求積公式計算結果的對比.可以看出，當z遠小于0（即場點遠離自由面）時，只要單元中布置的高斯點數(shù)目足夠多，傳統(tǒng)的Gauss-Legendre求積公式也可以獲得令人滿意的結果；但是當z趨于0時，即使在單元中布置了16×16個高斯點，Gauss-Legendre求積公式的結果仍然與理論解[18]存在偏差，而文中的半解析方法則可以在離散水平線段數(shù)目達到8時獲得足夠的精度.

在采用九節(jié)點高階面元對船體表面進行離散，并利用文中建立的公式對邊界積分方程式（3）中與GR和GF相關的面元積分進行計算后，便可求出邊界上的源強，進而得到流場速度勢.

2 運動響應與波浪增阻計算方法

將入射勢I和求得的輻射勢j（j=1， 2，...， 6）、繞射勢7代入水動力計算公式[19]，便可計算得到作用在船體上的波浪激勵力和水動力系數(shù).在此基礎上，求解頻域運動方程式（10）[20]，便可得到船舶的運動響應Xj.

式中：mij為船舶質(zhì)量慣性力系數(shù)；cij為恢復力系數(shù)；μij和λij分別為附加質(zhì)量和阻尼系數(shù)；fIi和fDi分別為入射力和繞射力的復數(shù)幅值.

波浪中航行船舶除了會產(chǎn)生各種搖蕩運動之外，受到的阻力相比于靜水中也會有所增加.這些附加阻力會增加船舶的燃料消耗，影響船舶保持航速的能力.因此，準確預報波浪增阻，對船舶的設計和使用至關重要.為了便于計算和分析，波浪增阻通?？梢员环譃檩椛湓鲎韬屠@射增阻兩部分.

輻射增阻的計算公式可以根據(jù)能量守恒原理得到[21]：

式中：RRAW代表輻射增阻；XiR和XiI分別代表Xi的實部和虛部，XjR和XjI分別代表Xj的實部和虛部.Xi（或Xj）原則上可通過求解頻域運動方程式（10）獲得，但是這個值是相對于船舶在靜水中的平衡位置而言的.為了考慮船體周圍波面變化的影響，文中參照文獻[21]中的做法對船舶垂蕩運動幅值加以修正：

式中：X′3為修正后的垂蕩運動幅值；xcg和ycg分別代表重心處的x和y坐標；T*為船舶平均吃水，近似等于Aw（和Aw分別為船舶的排水體積和水線面面積）.

至于繞射增阻，通過移除Salvesen平均二階力計算公式[13]中與輻射勢相關項，保留其中與繞射勢相關的項來計算，具體表達式為：

式中：RDAW代表繞射增阻；φ*I為入射勢I的共軛復數(shù).

將由式（11、13）分別計算得到的輻射增阻和繞射增阻相加，便可得到船舶在波浪中航行時所受到的總的波浪增阻.然而需要注意的是，上述計算公式均是在線性勢流理論框架內(nèi)建立的，當入射波長較短時，并不適用.為了彌補它們的不足，日本海上技術安全研究所（NMRI）提出了一種短波增阻半經(jīng)驗計算公式為[22]：

式中：RSAW代表短波時的波浪增阻；B為船寬；系數(shù)αd，αU和Bf的表達式見文獻[22].

前面的輻射增阻、繞射增阻公式（11、13）只適用于波長較長的工況，而半經(jīng)驗公式（14）只適用于波長較短的工況.為了克服其局限性，文中參考文獻[23]的做法，提出將長波、短波增阻公式相結合的混合計算方法，以便在全頻率范圍均能夠獲得準確的波浪增阻預報結果.

RAW=（1－αd）（RRAW+RDAW）+αdRSAW（15）

3 數(shù)值結果與分析

基于上述理論和方法，采用Fortran開發(fā)相應的數(shù)值計算程序.為了對程序進行驗證，首先以文獻[23]為研究對象，計算船舶運動響應和波浪增阻，并與試驗和其它數(shù)值方法的預報結果進行對比分析.Modified Wigley是一種數(shù)學船型，船長L為2 m，船寬B為0.3 m，吃水D為0.125 m.圖6為modified Wigley船的型線圖，從圖中可以看出，modified Wigley船是一種直壁型船，橫剖線在水線處與自由面垂直.

圖7為modified Wigley以速度Fr=0.2在不同波長（波長λ船長L比范圍為0.2～2.0，間隔為0.1）的規(guī)則波中迎浪航行時的垂蕩和縱搖運動RAO（response amplitude operator）.為了進行對比分析，圖中除了文中的計算結果，還包括了試驗[24]、基于Gauss-Legendre求積公式的常值元法和高階面元法的結果.從圖中的結果可以看出，常值元法和文中建立的高階半解析移動脈動源法的計算結果差異不大（文中高階方法得到的垂蕩運動響應在峰值附近的精度略高），并且總的來說均與試驗值吻合良好（這主要是因為modified Wigley船是一種簡單的數(shù)值船型，常值元離散也可以獲得比較好的精度）.但是由于船體在采用高階面元進行離散時，最上面一層的控制點位于自由面，而傳統(tǒng)數(shù)值積分法在場、源點接近自由面時無法獲得精確的格林函數(shù)面積分值（圖5），基于Gauss-Legendre求積公式高階面元法計算得到的運動響應呈現(xiàn)出明顯的振蕩特性.

圖8為modified Wigley在Fr=0.2時的波浪增阻.可以看出，類似于前面運動響應的結果，基于常值元法和文中高階半解析移動脈動源法計算得到波浪增阻總的來說也與試驗值吻合良好（其中高階方法的精度在共振頻率附近更高）.而傳統(tǒng)基于Gauss-Legendre求積公式的高階面元法由于無法獲得穩(wěn)定的運動響應結果，計算出的波浪增阻也呈現(xiàn)出明顯的振蕩特性.

以上結果是在處理器為AMD EPYC 7H12（主頻為2.6 GHz）工作站上采用40個線程并行計算獲得的，平均每個頻率的計算時間在17.6 s左右.可以看出，文中程序的計算效率較高，能夠滿足工程需求.

前面對modified Wigley船的計算分析表明，文中建立的高階半解析移動脈動源法以及基于此的波浪增阻方法對于簡單數(shù)學船型可以獲得令人滿意的運動響應和波浪增阻預報結果.為了進一步驗證該方法在實際船舶中的適用性，將進一步對S175集裝箱船的波浪增阻進行計算.S175是ITTC標準船模，船長L為175 m，船寬B為25.4 m，吃水D為9.5 m，排水量Δ為24 742 t，重心位于船尾部分且距離船中2.48 m，重心高度KG為9.52 m，縱搖方向轉動慣量為43 644 888 t·m2.圖9為S175船的九節(jié)點單元離散網(wǎng)格，可以看出，不同于modified Wigley，S175在自由面附近存在一定的外飄.

圖10為S175以航速Fr=0.25在規(guī)則波中迎浪航行時，波浪增阻的計算值和試驗數(shù)據(jù)[25-26]的對比.計算的波長船長比范圍為0.5～2.5，間隔為0.1.除了文中高階半解析方法的結果，圖中數(shù)值結果還包括基于半解析積分公式的常值元法[12]、基于Gauss-Legendre求積公式的常值元法和高階面元法的計算值.從圖中結果可以看出，基于Gauss-Legendre求積公式的常值元法和高階面元法計算得到的波浪增阻比其他兩種數(shù)值方法的結果更為振蕩.這主要是因為S175船在水線附近存在外飄，在外飄處，即使采用常值元進行離散，它的控制點距離自由面也非常接近，從而導致很難采用數(shù)值積分公式獲得準確的格林函數(shù)面積分值.對比基于半解析積分公式的常值元法和文中方法計算得到的波浪增阻可以發(fā)現(xiàn)，雖然前者也可以獲得穩(wěn)定并且總體令人滿意的結果，但是由于常值元離散的固有缺陷，它的計算精度在共振頻率附近明顯不如文中的高階方法.

基于文中建立的高階半解析移動脈動源法，進一步對S175船以速度Fr=0.15， 0.20， 0.30迎浪航行時的波浪增阻進行計算，結果如圖11.從圖中可以看出，數(shù)值結果與試驗值吻合良好，表明本文方法適用于不同航速下船舶波浪增阻的預報.此外，對比圖中不同航速下的結果可以發(fā)現(xiàn)，S175船受到的波浪增阻存在明顯的航速效應，它的峰值及其對應的無因次波長均會隨著航速的增大而增大.

為了進一步分析文中方法的可行性和可靠性，在前面規(guī)則波中波浪增阻計算的基礎上，下面進一步采用譜分析方法對S175在實際海況中航行時的波浪增阻進行預報.數(shù)值計算采用與Nakamura 和文獻[26]相同的波浪譜，有義波高H1/3分別取10.78、9.99、10.56、10.04 cm，相應的平均周期T0分別取1.159、1.413、1.562、1.694 s.圖12為S175以速度Fr=0.15、0.20、0.30分別在4種海況中迎浪航行時的波浪增阻平均值隨平均周期T0的變化.從圖中可以看出，雖然部分航速下的數(shù)值結果在平均周期較短時存在一定的誤差（這可能是試驗本身的誤差或者是線性疊加的譜分析方法與實際物理情況之間存在差別導致），但總的來說預報值與試驗值吻合良好，在工程可以接受的范圍內(nèi).

4 結論

（1） 與傳統(tǒng)Gauss-Legendre求積公式相比，本文建立的半解析積分方法由于在計算與移動脈動源格林函數(shù)航速相關部分GF相關的面積分時，對沿水平方向的積分采用了解析計算公式，從而避免了GF的劇烈振蕩特性所導致的自由面附近邊界積分方程影響系數(shù)計算不準確問題；

（2） 在波浪中航行船舶的運動響應和波浪增阻數(shù)值計算中，文中高階半解析移動脈動源法比相應的基于Gauss-Legendre求積公式的面元法更為穩(wěn)定（尤其是對于存在外飄的船型），比常值元法的預報精度更高；

（3） 無論是簡單的數(shù)學船型Modified Wigley，還是集裝箱船S175，基于高階半解析移動脈動源法計算得到的不同航速下的波浪增阻結果均與試驗值吻合良好，表明本文建立的數(shù)值方法在波浪引起的船舶阻力增加預報中，對于不同船型、不同航速均有廣泛的適用性.

參考文獻（References）

［1］ GERRITSMA J， BEUKELMAN W. Analysis of the resistance increase in waves of a fast cargo ship [J]. International Shipbuilding Progress， 1972， 19（217）： 285-293.

[2] SALVESEN N. Added resistance of ships in waves[J]. Journal of Hydronautics， 2012， 12（1）： 24-34.

[3] NEWMAN J N. Algorithms for the free-surface Green function [J]. Journal of Engineering Mathematics， 1985， 19（1）： 57-67.

[4] PAPANIKOLAOU A， ZARAPHONITIS G. On an improved method for the evaluation of second-order motions and loads on 3D floating bodies in waves [J]. Ship Technology Research， 1987， 34： 170-211.

[5] FANG M C， CHEN G R. On the nonlinear hydrodynamic forces for a ship advancing in waves [J]. Ocean Engineering， 2006， 33（16）： 2119-2134.

[6] IWASHITA H， OHKUSU M. Hydrodynamic forces on a ship moving with forward speed in waves [J]. Journal of the Society of Naval Architects of Japan， 1989，166： 187-205.

[7] MAURY C， DELHOMMEAU G， BA M， et al. Comparison between numerical computations and experiments for seakeeping on ship models with forward speed[J]. Journal of Ship Research， 2003， 47（4）： 347-364.

[8] YAO C B， DONG W C. Study on fast integration method for Bessho form translating pulsating source Green′s function distributing on a panel [J]. Ocean Engineering， 2014， 89： 10-20.

[9] HONG L， ZHU R C， MIAO G P， et al. Study on Havelock form translating pulsating source Green′s function distributing on horizontal line segments and its applications [J]. Ocean Engineering， 2016， 124： 306-323.

[10] YANG Y T， ZHU R C， HONG L， et al. A semi-analytical high-order translating-pulsating source method for forward-speed ship motions [J]. Ocean Engineering， 2019， 182： 627-644.

[11] CHEN X B. Highly oscillatory properties of unsteady ship waves [J]. Journal of Mechanical Engineering Science， 2000， 214（6）：813-823.

[12] HONG L， ZHU R C， MIAO G P， et al. An investigation into added resistance of vessels advancing in waves [J]. Ocean Engineering， 2016， 123： 238-248.

[13] SALVESEN N. Second-order steady-state forces and moments on surface ships in oblique regular waves[C]∥International Symposium on Dynamics of Marine Vehicles and Structures in Waves. University College， London：[s.n.]， 1974.

[14] 朱仁傳， 繆國平. 船舶在波浪上的運動理論[M]. 上海：上海交通大學出版社， 2019.

[15] WU G X， EATOCK T R. A Green′s function form for ship motions at forward speed [J]. International Shipbuilding Progress， 1987， 34（398）：189-196.

[16] 邢瑞山， 卿光輝. 九節(jié)點 Hamilton 等參元列式[J]. 廣東工業(yè)大學學報， 2012， 29（1）： 27-31.

[17] CHEN X， ZHU R C， MA C， et al. Computations of linear and nonlinear ship waves by higher-order boundary element method [J]. Ocean Engineering， 2016， 114： 142-153.

[18] MAURY C. Etude du problème de tenue à la mer avec vitesse d′avance quelconque par une méthode de singularité s de Kelvin [D]. cole France：cole Centrale de Nantes， 2000.

[19] YANG Y T， ZHU R C， Hong L. A frequency-domain hybrid HOBEM for motion responses and added resistance of ships sailing in head and oblique waves [J]. Ocean Engineering， 2019， 194： 106637.

[20] 楊云濤. 航行船舶運動的三維頻域高階面元法數(shù)值計算研究[D]. 上海：上海交通大學， 2020.

[21] 洪亮. 移動脈動源格林函數(shù)法及航行船舶在波浪中阻力增加的研究[D]. 上海：上海交通大學， 2017.

[22] TSUJIMOTO M， SHIBATA K， KURODA M， et al. A practical correction method for added resistance in waves [J]. Journal of the Japan Society of Naval Architects and Ocean Engineers， 2008， 8： 177-184.

[23] GUO B J， STEEN S. Evaluation of added resistance of KVLCC2 in short waves [J]. Journal of Hydrodynamics， 2011， 23（6）： 709-722.

[24] KASHIWAGI M. Hydrodynamic study on added resistance using unsteady wave analysis [J]. Journal of Ship Research， 2013， 57（4）： 220-240.

[25] FUJII H. Experimental study on the resistance increase of a ship in regular oblique waves [C]∥Proceedings of the 14th ITTC. Ottawa，Canda：[s.n.]，1975.

[26] NAKAMURA S， NAITO S. Propulsive performance of a container ship in waves [J]. Journal of the Society of Naval Architects of Japan， 1977，15：24-48.