規(guī)則波中船舶Froude-Krylov 力的解析積分計(jì)算

張 騰，任俊生，張秀鳳，梅天龍

(1. 大連海事大學(xué)航海學(xué)院，大連 116026；2. 上海交通大學(xué)船舶海洋與建筑工程學(xué)院，上海 200240)

船舶航行在大洋中時(shí)時(shí)刻刻遭受波浪的作用，船舶在波浪中的大幅度搖蕩運(yùn)動(dòng)也會(huì)對(duì)船體結(jié)構(gòu)造成破壞，從而大大降低船舶的適航性，更為嚴(yán)重的會(huì)導(dǎo)致海難的發(fā)生．因此，準(zhǔn)確計(jì)算船體所受的波浪作用力對(duì)船舶操作、設(shè)計(jì)及運(yùn)營(yíng)顯得尤為重要．

自20 世紀(jì)50 年代以來，船舶在波浪中所受的Froude-Krylov 力及運(yùn)動(dòng)的研究取得了一定的進(jìn)展，主要體現(xiàn)在從無限水深到有限水深及從數(shù)值計(jì)算到解析計(jì)算方法的轉(zhuǎn)變．在無限水深范疇內(nèi)，Korvin-Kroukovsky 等[1]基于二維理論，采用切片法將船體離散成若干等截面的切片，數(shù)值計(jì)算出船舶各個(gè)橫剖面由于入射波作用引起的F-K 力，沿船長(zhǎng)方向積分可得到總的F-K 力，但其并未考慮到由入射波引起的繞射力．Salvesen 等[2]進(jìn)一步提出了STF 法，該理論有效地考慮由入射波引起的繞射力對(duì)船體的作用，其運(yùn)動(dòng)預(yù)報(bào)精度有了較大提高．基于三維范疇內(nèi)，Beck等[3]及孫葳等[4]將船體平均濕表面離散成一系列四邊形平面面元，F(xiàn)-K 力可由高斯積分定理在平面面元上進(jìn)行直接積分得到．鄒元杰等[5]基于“高頻低速”假定對(duì)零航速格林函數(shù)進(jìn)行航速修正，并分析了波浪中行駛船舶在不同波長(zhǎng)處所受的水動(dòng)力及F-K 力對(duì)船舶運(yùn)動(dòng)的影響．

以上文獻(xiàn)中F-K 力計(jì)算都是基于無限水深及數(shù)值積分進(jìn)行的，實(shí)際上，水深對(duì)船舶運(yùn)動(dòng)具有很大的影響，比如船舶破損事故經(jīng)常在有限水深海域發(fā)生，因此F-K 力的計(jì)算需要考慮水深的影響是必不可少的．在有限水深下，基于二維理論，賀五洲等[6]采用切片法預(yù)報(bào)船舶所受到的F-K 力、輻射力及繞射力．基于三維頻域理論，馮乾棟等[7]基于有限水深三維頻域格林函數(shù)，對(duì)有限水深中的破損船舶波浪載荷進(jìn)行數(shù)值計(jì)算．在三維時(shí)域內(nèi)，劉昌鳳[8]采用數(shù)值積分方法計(jì)算船舶在有限水深及無限水深下所受到的F-K 力，并基于三維時(shí)域Green 函數(shù)[9]，成功地預(yù)報(bào)了船舶在規(guī)則波中的運(yùn)動(dòng)．上述不管有限水深還是無限水深中F-K 力的計(jì)算，都是基于數(shù)值積分方法，對(duì)于具有強(qiáng)烈波動(dòng)性的F-K 力的計(jì)算，即使將船體劃分大量的面元，仍然存在數(shù)值積分誤差．Rodrigues等[10]采用解析壓力積分表達(dá)式對(duì)無限水深以正弦表達(dá)式形式的規(guī)則入射波產(chǎn)生的F-K 力進(jìn)行精確計(jì)算，其船體表面網(wǎng)格劃分與計(jì)算水動(dòng)力相同的粗網(wǎng)格相同，精度卻毫無損失，避免了數(shù)值積分計(jì)算F-K 力產(chǎn)生的數(shù)值誤差．但是該研究進(jìn)行F-K 力解析積分計(jì)算時(shí)僅僅考慮無限水深這種特殊的情形，并沒有將水深因素考慮進(jìn)去，在近海處，由有限水深帶來的淺水效應(yīng)對(duì)F-K 力的計(jì)算有較大的影響．

為避免數(shù)值積分方法計(jì)算F-K 力帶來的誤差及減少船體表面面元?jiǎng)澐謹(jǐn)?shù)量，并將流域水深因素考慮進(jìn)去，基于格林定理，本文詳細(xì)推導(dǎo)了任意水深下FK 力在四邊形平面面元上的解析積分表達(dá)式．研究了流域水深在不同入射波浪向角及波長(zhǎng)下對(duì)作用在WigleyⅠ型船[11]船體上的F-K 力的影響，并進(jìn)一步基于三維時(shí)域Green 函數(shù)建立了船舶任意水深下的波浪激勵(lì)力計(jì)算模型．最后，對(duì)迎浪航行在規(guī)則波中WigleyⅠ船舶受到的波浪激勵(lì)力進(jìn)行預(yù)報(bào)．

1 數(shù)學(xué)模型的建立

如圖1 所示，參考坐標(biāo)系Oxyz 固定在以速度U前進(jìn)的船體上，原點(diǎn)O 為船舶中點(diǎn)，Oxy 與靜水面重合，Ox軸正方向指向船首，Oy 軸正方向指向船舶左舷，Oz軸正方向垂直向上．面元局部坐標(biāo)系O′x′yz′ ′的原點(diǎn)O′為面元幾何中心，O′xy′ ′與四邊形面元平面重合，O z′ ′軸正向垂直于四邊形平面面元，且指向船體內(nèi)部．面元局部坐標(biāo)坐標(biāo)與參考坐標(biāo)系Oxyz 坐標(biāo) r ( x , y ,z)的轉(zhuǎn)換式為

圖1 參考坐標(biāo)系與面元坐標(biāo)系Fig.1 Reference coordinate system and panel coordinate system

1.1 F-K力的解析積分表達(dá)式推導(dǎo)

假定流體不可壓縮、無黏性，且流動(dòng)是無旋的，不考慮流體表面張力效應(yīng)．水深為d的入射波速度勢(shì)及入射波波高在參考坐標(biāo)系下的表達(dá)式分別為

式中：0η為入射波波幅；g為重力加速度；ω為入射波原頻率；k =ω2g 為波數(shù)；β為浪向角，其傳播方向與Ox軸正方向的夾角，迎浪為為遭遇頻率．

式中ρ為流體密度．

將船體濕表面處理為N塊四邊形平面面元，則第i塊面元上的入射波壓力積分為

式中Si為第i 塊面元面積．

經(jīng)過推導(dǎo)，則IF可寫為

式中 fI1及 fI2可分別寫為

第i塊四邊形面元上第j條邊在面元坐標(biāo)系內(nèi)可以參數(shù)化得

則第i塊面元第j條邊上的C、D 及D′ 關(guān)于參數(shù)α的表達(dá)式分別為

fI1可進(jìn)一步寫作

式中 ΔC Eij和 ΔS Eij可表示為

式中SEij及CEij可表示為

同理，可以采用相同方法求取fI2．最終將fI1及代入式(6)，解得入射波壓力積分IF．

同理，可以采用解析積分表達(dá)式求解F-K 力及力矩．值得注意，當(dāng)水深 d→∞時(shí)，此時(shí)IF可以表示為

注意：式(22)為IF為無限水深時(shí)F-K 力的解析積分表達(dá)式，相對(duì)于有限水深F-K 力的解析表達(dá)式(6)，式(22)中不包含水深d 的影響．

1.2 繞射問題求解

在Oxyz 參考坐標(biāo)下的入射波遇到船體反射產(chǎn)生繞射勢(shì)DΦ，其滿足如下條件：

式中：Ω為流域；SF為自由液面；SH為船體平均濕表面；SB為水底面積；n 為單位法向量，指向流域外部．令分別代表船舶縱蕩、橫蕩、垂蕩、橫搖、縱搖及首搖.

時(shí)域格林函數(shù)用于求解滿足繞射勢(shì)DΦ邊界方程

G0與的數(shù)值求解可由文獻(xiàn)[9，12-13]所得，其繞射勢(shì)ΦD滿足的邊界積分方程[8]為

式中：Γ為船體與自由液面的交線；為垂直水線指向船體內(nèi)部的單位法向量．

采用常數(shù)面元法求解式(29)，最終求解出繞射勢(shì)

DΦ，則由繞射勢(shì)產(chǎn)生的壓力 pD為

則i運(yùn)動(dòng)模態(tài)上的繞射力FDi為

i模態(tài)上的頻域波浪激勵(lì)力 Fwi可由 FDi與 FIi合力經(jīng)過傅里葉變換得到[8]．

完成建模工作后，有必要對(duì)Rodrigues 方法[10]與本文方法做一個(gè)比較，如表1 所示．

如表1 所示，本文方法相對(duì)于Rodrigues 方法，適用范圍更為廣泛，主要體現(xiàn)在：①本文采用解析積分方法計(jì)算任意水深下的F-K 力，而Rodrigues 采用解析積分方法只計(jì)算無限水深下的F-K 力；②本文采用直接時(shí)域方法計(jì)算流體擾動(dòng)力，可以有效避免采用間接時(shí)域方法在高頻處的數(shù)值誤差；③本文入射波波高以余弦形式表示，而Rodrigues 方法以正弦形式表示，最終兩者F-K 力解析積分表達(dá)式不僅體現(xiàn)在有無水深項(xiàng)的不同，還體現(xiàn)在式(17)及式(21)括號(hào)項(xiàng)中表達(dá)式的不同．

2 數(shù)值模擬與結(jié)果分析

2.1 面元上F-K力的解析積分計(jì)算驗(yàn)證與分析

首先選取1 塊1 m×1 m 的正方形面元進(jìn)行研究，其4 個(gè)頂點(diǎn)坐標(biāo)在參考坐標(biāo)系中依次為(-0.5，0，-1)、(0.5，0，-1)、(0.5，0，1)和(-0.5，0，0)．該面元前進(jìn)速度為U ＝1 m/s，入射波波長(zhǎng)λ=5 m，波幅0η＝0.1 m，水深d ＝6 m．在t＝3 s 時(shí)刻，對(duì)浪向角β為及在該面元上F-K 力由壓力積分解析表達(dá)式分別求得為446.637 N、-88.462 N(保留到0.001 N)．對(duì)面元上FK 力采用高斯數(shù)值積分進(jìn)行計(jì)算，并進(jìn)一步對(duì)面元進(jìn)行如圖2 的剖分，直到與解析解求得的結(jié)果相同為止．

圖2 面元剖分示意Fig.2 Panel subdivision schematic

由表2 可知，當(dāng)劃分次數(shù)為3 時(shí)，數(shù)值積分解與解析積分解具有相同的精度，此時(shí)面元被劃分為64個(gè)子面元，數(shù)值積分所需要的面元數(shù)目遠(yuǎn)遠(yuǎn)大于解析積分．對(duì)于形狀比較規(guī)則的箱形體船，比如駁船[10]，需要6 塊面元即可求出準(zhǔn)確的F-K 力值，極大地提高了計(jì)算效率．

表2 t＝3 s面元上F-K壓力積分值Tab.2 Values of F-K pressure integration on the panels at t＝3 s N

2.2 船體F-K力的計(jì)算與分析

本文采用WigleyⅠ型船驗(yàn)證數(shù)學(xué)模型的可靠性與準(zhǔn)確性，其主要參數(shù)如表3 所示．并將船體平均吃水一下的船體劃分為80×4 個(gè)面元，如圖3 所示．

表3 WigleyⅠ型船主尺度Tab.3 Main particulars of Wigley Ⅰ hull

圖3 Wigley Ⅰ型船體面元?jiǎng)澐諪ig.3 Panel distribution for Wigley Ⅰ hull

由圖4 和圖5 可得，“數(shù)值積分”與“解析積分”的F-K 力時(shí)間歷程完全重合，從而驗(yàn)證了本文開發(fā)的F-K 力解析積分算法的可靠性．

圖4 無量綱垂蕩F-K力時(shí)間歷程(=0.5，d =0.5λ，F(xiàn)r=0.2)Fig.4 Time history of non-dimensional heave F-K forces(=0.5，d =0.5λ，F(xiàn)r=0.2)

圖5 無量綱縱搖F-K力時(shí)間歷程(=0.5，d =0.5λ，F(xiàn)r=0.2)Fig.5 Time history of non-dimensional pitch F-K forces(=0.5，d =0.5λ，F(xiàn)r=0.2)

WigleyⅠ型船前后左右?guī)缀螌?duì)稱，如圖 6 所示．本文選取第1象限中及3 個(gè)浪向角進(jìn)行研究.在船舶航速下，對(duì)于不同入射波波長(zhǎng)船長(zhǎng)比λL 處，本文對(duì)影響F-K 力的流域水深d做進(jìn)一步分析與探討．其中，無量綱水深為

由圖7～圖12 可得，在不同入射波波長(zhǎng)及浪向角下，當(dāng)無量綱水深較小時(shí)，水深對(duì)F-K 力影響較為顯著，不可忽略．當(dāng)無量綱水深時(shí)，船體所受的 F-K 力逐漸趨于無窮水深所受到的F-K 力.

圖9 無量綱垂蕩F-K力幅值(β =)Fig.9 Amplitude of non-dimensional heave F-K forces(β =)

圖11 無量綱垂蕩F-K力幅值(β =)Fig.11 Amplitude of non-dimensional heave F-K forces(β =)

2.3 繞射問題分析與驗(yàn)證

圖6 浪向角示意Fig.6 Diagram of incident wave directions

圖8 無量綱縱搖F-K力幅值(β =π)Fig.8 Amplitude of non-dimensional pitch F-K forces(β =π)

圖10 無量綱縱搖F-K力幅值()Fig.10 Amplitude of non-dimensional pitch F-K forces()

圖12 無量綱縱搖F-K力幅值(β =)Fig.12 Amplitude of non-dimensional pitch F-K forces(β =)

為與Journée[11]試驗(yàn)結(jié)果做比對(duì)，本文只對(duì)無限水深(即水深d 趨于無窮的情況)的WigleyⅠ型船進(jìn)行垂蕩、縱搖2 個(gè)自由度的波浪激勵(lì)力數(shù)值模擬，船舶前進(jìn)速度分別為Fr=0.2與Fr=0.3，基于本文計(jì)算模型所得結(jié)果表示為“本文方法”，基于Journée 試驗(yàn)所得結(jié)果為“試驗(yàn)結(jié)果”．

垂蕩波浪激勵(lì)力為Fw3，F(xiàn)w3無量綱形式為縱搖波浪激勵(lì)力為其無量綱形式為無量綱波數(shù)為

由圖13～圖16 可知，在大部分頻率范圍內(nèi)，WigleyⅠ型船所受的垂蕩波浪激勵(lì)力與縱搖波浪激勵(lì)力的頻率響應(yīng)函數(shù)均與試驗(yàn)值吻合良好，其中垂蕩波浪激勵(lì)力的預(yù)報(bào)效果優(yōu)于縱搖波浪激勵(lì)力．結(jié)果驗(yàn)證了本文提出算法的可靠性與有效性．

圖13 無因次垂蕩波浪激勵(lì)力幅值(Fr=0.2)Fig.13 Non-dimensional heave amplitude of wave exciting force(Fr=0.2)

圖14 無因次縱搖波浪激勵(lì)力幅值(Fr=0.2)Fig.14 Non-dimensional pitch amplitude of wave exciting force(Fr=0.2)

圖15 無因次垂蕩波浪激勵(lì)力幅值(Fr=0.3)Fig.15 Non-dimensional heave amplitude of wave exciting force(Fr=0.3)

此外，由于本文采用自由面時(shí)域格林函數(shù)計(jì)算波浪激勵(lì)力，該方法基于Neumann-Kelvin 線性化自由面條件，其不能準(zhǔn)確考慮航速效應(yīng)的影響，在個(gè)別頻率處計(jì)算結(jié)果與試驗(yàn)結(jié)果偏差較大．為了進(jìn)一步提高波浪激勵(lì)力的數(shù)值預(yù)報(bào)精度，可以采用簡(jiǎn)單Rankine 法，該方法不僅能準(zhǔn)確地考慮航速效應(yīng)，也能計(jì)入移動(dòng)興波對(duì)船體所受波浪激勵(lì)力的影響，從而進(jìn)一步提高船體受力及運(yùn)動(dòng)預(yù)報(bào)的數(shù)值精度．

圖16 無因次縱搖波浪激勵(lì)力幅值(Fr=0.3)Fig.16 Non-dimensional pitch amplitude of wave exciting force(Fr=0.3)

3 結(jié) 論

(1)本文系統(tǒng)地推導(dǎo)了任意水深下船體所受F-K力的解析積分方法，可以避免由于數(shù)值積分所產(chǎn)生的誤差，并大大減少船舶所需剖分的面元數(shù)量．對(duì)于形狀規(guī)則的船體，比如駁船或者海洋鉆井平臺(tái)，只需要?jiǎng)澐謹(jǐn)?shù)量級(jí)面元即可滿足計(jì)算要求．

(2)在一定入射波長(zhǎng)及浪向角下，當(dāng)流域水深較小時(shí)，水深因素對(duì)WigleyⅠ型船船體所受的F-K 力有較為顯著的影響，不可忽視．隨著水深的增加，特別當(dāng)水深船長(zhǎng)比大于0.5 時(shí)，有限水深中WigleyⅠ型船船體所受到的F-K 力越來越趨近于無限水深中船體所受到F-K 力．

(3)通過對(duì)無限水深中的以不同航速迎浪前進(jìn)的WigleyⅠ型船進(jìn)行計(jì)算仿真，其波浪激勵(lì)力與試驗(yàn)值吻合良好，進(jìn)一步證明了本文任意水深下船體所受F-K 力解析積分計(jì)算模型的精確性與可靠性．