規(guī)則波浪中艦船搖蕩耦合切片計算方法

朱軍，夏齊強，王智宇，朱韜，黃昆侖，葛義軍

1海軍工程大學(xué)艦船與海洋學(xué)院，湖北武漢430033

2海軍研究院，上海200235

0 引 言

切片理論方法首次發(fā)表于20世紀50年代［1］，之后衍生出了各種不同的切片方法，并被應(yīng)用到某些船舶上用于確定速度的效應(yīng)。Salvesen等［2］提出了線性切片理論方法，該方法可用于描述滿足交叉耦合項的對稱定理。在線性和搖蕩諧振假設(shè)下，線性切片勢流理論方法滿足物面和自由表面的邊界條件，盡管該方法采用勢流理論計算了水動力并提出了六自由度的計算模型，但由于物面邊界條件是在平均位置上得到滿足的，因此垂直面和水平面實際上是分離的。北約耐波性一般準則和共同程序?qū)＜医M［3］的研究表明，采用各種計算機代碼計算出的驅(qū)逐艦垂直平面的運動結(jié)果非常一致。Ogilvie和 Tuck［4］的研究表明，當(dāng)波長在船寬量級的高頻范圍內(nèi)時，切片方法是嚴格有效的。實踐經(jīng)驗表明，對于船長與波長等量級的情況，切片方法能給出良好的結(jié)果。Timman等［5］采用攝動展開方法，將原本在瞬時物面上滿足的邊界條件變換成了在物面平均位置上滿足的邊界條件，這雖然給求解帶來了極大的方便，但卻使得垂蕩、縱搖與橫搖實際上在非線性條件下被解耦了。

為了計算甲板上浪和艏底砰擊等非線性波浪載荷問題，需采用線性切片理論方法預(yù)報船體的搖蕩運動。Kaplan等［6-8］］將線性切片方法預(yù)報的運動響應(yīng)作為輸入，在頻域中有效解決了運動問題，而載荷預(yù)報則在時域中執(zhí)行。Cong等［9］采用脈沖響應(yīng)函數(shù)方法計算了輻射波浪力和衍射波浪力，同時修正了入射波浪壓力分布的Smith效應(yīng)。Gu等［10］采用線性切片理論方法預(yù)報了垂直面的搖蕩（特指垂蕩和縱搖）運動，并直接在時域中計算了運動、甲板上浪與砰擊載荷。林超友和朱軍等［11-12］為了預(yù)報波浪中的操縱運動，采用勢流理論方法計算了波浪力。

線性切片理論方法本質(zhì)上沒有耦合船體的垂直面（垂蕩和縱搖）和水平面（橫搖）運動，這是因為在計算搖蕩速度勢和入射衍射波浪速度勢時，船體切片只取了其靜水的平均位置，導(dǎo)致垂蕩、縱搖與橫搖的相互影響被忽略。盡管速度勢攝動方法為非線性方法，但因其船體切片只取了靜水平均位置，故在本質(zhì)上并沒有改善垂蕩、縱搖與橫搖的耦合性。

為了能計入船體垂蕩、縱搖與橫搖的運動耦合效應(yīng)，本文擬基于廣義縱傾角和廣義吃水增量參數(shù)［13］，在船體瞬時運動狀態(tài)（ΔT，θ和φ）和瞬時波浪狀態(tài)下提取船體切片，以滿足瞬時波面壓力為零的條件修正壓力分布的計算公式，采用文獻［13］所應(yīng)用的圖形面域技術(shù)計算切片的靜水力與波浪擾動力（傅汝德—克雷洛夫力）之和，采用經(jīng)驗方法估算瞬時狀態(tài)下的船體切片慣性和阻尼水動力，并由所建立的搖蕩耦合運動方程預(yù)報規(guī)則波浪中的船體搖蕩運動。

1 坐標(biāo)系及波面方程

1.1 坐標(biāo)系及廣義參數(shù)

坐標(biāo)系如圖1所示，圖中χ為浪向角，t為時間。慣性坐標(biāo)系為E-ξηζ，坐標(biāo)平面E-ξη與靜水面重合；船體運動坐標(biāo)系為O-xyz，坐標(biāo)平面O-xy與靜水面重合，沿慣性坐標(biāo)縱向軸以速度U等速運動；船體固定坐標(biāo)系O′-xO′yO′zO′與船體運動坐標(biāo)保持等速前進運動。

將船體固定坐標(biāo)系繞Ox軸轉(zhuǎn)動φ角，然后沿Oz軸平移ΔT（吃水增量），再繞Oy軸（靜水面軸線）旋轉(zhuǎn)θ角，運動坐標(biāo)系（x，y，z）和固定坐標(biāo)系（xO′，yO′，zO′）的轉(zhuǎn)換關(guān)系（圖2）為：

值得注意的是，這里的吃水T和橫傾角φ與通常的定義是一致的，但吃水增量ΔT為垂直于靜水面的變化量，縱傾角θ為繞船舯剖面與靜水面交線的轉(zhuǎn)動角度，文獻［13］分別稱其為廣義吃水增量ΔT和廣義縱傾角θ。

1.2 波面方程

令波長為λ，波高為Hw，浪向角為χ，慣性坐標(biāo)系下的波面方程為：

將式（2）代入到式（1）的第3項中，得到船體固定坐標(biāo)系下的波面方程為

這里的x，y由式（1）的第1項和第2項確定。顯然，這是一個隱式方程，方程右端前4項是靜水面方程在船體坐標(biāo)系下的表達式，第5項是波面相對靜水面的升高量。

2 船體搖蕩耦合動力學(xué)方程及水動力

2.1 船體搖蕩耦合運動方程

令瞬時波面下船體靜水力與入射波浪力之和的水動力為FFK+B，即將靜水力與傅汝德—克雷洛夫波浪力作為一個整體計算，不難寫出船體橫搖的二階運動方程式為

在船體縱向平面內(nèi)，僅考慮縱搖與垂蕩運動的相互耦合項，按類似方法可得到簡化后的垂蕩與縱搖耦合運動方程式為

式中，P為船體重力。水動力FFK+B作用點（xFK+B,yFK+B,zFK+B）和船體重心（xG,yG,zG）系指船體運動坐標(biāo)，由式（1）確定與船體坐標(biāo)系下的轉(zhuǎn)換關(guān)系。

式（4）和式（5）構(gòu)成了船體垂蕩與縱搖、橫搖耦合的搖蕩運動方程。

2.2 瞬時波面下船體靜水力與波浪擾動力之和

普通線性切片方法是將搖蕩復(fù)原力與波浪擾動力分開處理，且船體切片只取平均水線位置，這導(dǎo)致波峰和波谷通過切片時無差異化。本文基于切片處的瞬時波面方程（3），根據(jù)切片瞬時波面下壓力分布積分計算水動力，并將靜水力與波浪擾動力之和作為整體處理。由勢流理論可知，小振幅規(guī)則波的壓力分布為：

式中，ρ為水的質(zhì)量密度。該壓力分布式忽略了速度平方的二階小量，瞬時波表面壓力等于零的條件并不滿足，波峰處的壓力估計偏大，波谷處的壓力估計偏小［14］，即史密斯效應(yīng)。為此，令

該轉(zhuǎn)換就是將波面作為壓力計算的起點。將瞬時波面下壓力分布改寫為

式中，A=Hw/2。式（9）滿足瞬時波面（z1=0時）壓力p=0，同時還滿足z1→-∞時p→-ρgz。在t時刻x處的船體切片，對式（9）的指數(shù)項進行線性化處理并用y=0近似，則近似的壓力分布為

于是，由式（10）對瞬時波面下的切片進行積分，然后再沿船體縱向積分，便可得到船體水動力FFk+B。

顯然，式（10）的壓力修正體現(xiàn)了波峰和波谷位于切片時的差異，同時，還包含了船體運動姿態(tài)（ΔT，θ和φ）的效應(yīng)，即垂蕩、縱搖與橫搖的耦合效應(yīng)。

2.3 船體水動力系數(shù)

式（4）和式（5）中的附加慣性水動力系數(shù)及阻尼水動力系數(shù)，將由所估算的切片水動力系數(shù)沿船體縱向積分計算得到。

1）船體橫搖附加慣性矩和阻尼系數(shù)。

橫搖附加慣性水動力和阻尼系數(shù)取切片靜水平均位置，按通常的經(jīng)驗公式估算：

式中：ρφ=CB，為橫搖慣性半徑，其中B為船寬，經(jīng)驗系數(shù)C取0.40；h為初穩(wěn)性高；Fn為傅汝德數(shù)。船體橫搖運動的非線性效應(yīng)明顯，Himeno［15］針對橫搖阻尼預(yù)測進行了綜述，更加有針對性的估算公式可參看文獻［15］。

2）船體垂蕩、縱搖附加慣性矩和阻尼系數(shù)。

式（5）中的船體水動力系數(shù)由瞬時波面下切片的水動力系數(shù)積分得到，瞬時波面下切片的垂蕩附加質(zhì)量和阻尼系數(shù)Δb估算為：

式中，系數(shù)C1和由遭遇頻率ωe、切片寬—吃水比Bn/dn及 切 片 面 積An的 系 數(shù)βn（βn=An/(Bndn)）按二維浮體試驗圖譜計算［16］。

值得注意的是，式（12）所估算的是瞬時波面下的船體切片，此時，切片狀態(tài)與船體縱傾角θ、橫傾角φ、吃水增量ΔT以及波浪瞬時相位有關(guān)，故估算所得到的船體水動力系數(shù)是時變的，同時也包含了垂蕩、縱搖與橫搖的耦合影響。

由式（12）沿船體縱向積分，不難得到式（5）等號左端的慣性和阻尼水動力系數(shù)。

2.4 船體搖蕩耦合運動方程特征

式（4）和式（5）構(gòu)成的非線性時變耦合的船體搖蕩耦合運動方程組具有時變和耦合的非線性特征。

1）水動力系數(shù)的時變性。

式（5）中的垂蕩和縱搖的附加慣性矩、阻尼系數(shù)及附加慣性靜矩共6個水動力系數(shù)均是針對瞬時波面下切片經(jīng)積分而得到，瞬時波面下切片的大小和形狀隨船體運動狀態(tài)（ΔT，θ和φ）和波浪相位的不同而變化。所以，式（5）具有系數(shù)時變的特征，而水動力系數(shù)的時變性則意味著運動方程包含了船體運動狀態(tài)和波浪相位的效應(yīng)。

2）垂蕩、縱搖和橫搖的運動耦合性。

水動力FFK+B為靜水力與波浪擾動力（傅汝德—克雷洛夫波浪力）之和，其值取決于瞬時船體運動狀態(tài)（ΔT，θ，φ）和波浪相位。水動力FFK+B既是構(gòu)成船體橫搖的激勵力，也是垂蕩、縱搖運動的激勵力，可見水動力FFK+B耦合了船體垂蕩、縱搖和橫搖的運動。

3 數(shù)值計算及分析

3.1 船型主尺度及波浪參數(shù)

計算艦船的水線長L=142 m，靜水吃水T=4.535 m，航速為0。計算的波浪參數(shù)為：

1）浪向角χ=0°，45°，90°，135°，180°，即頂浪、首斜浪、正橫浪、尾斜浪和順浪的典型浪向角。

2）波長船長之比λ/L=0.2～5.0，波浪頻率從低頻到高頻，涵蓋共振頻率。

3）波陡Hw/λ=0.01，0.04，對應(yīng)的最大波傾角α0=1.8°，7.2°。

3.2 數(shù)值計算及分析

數(shù)值計算初始運動參數(shù)：T=4.535 m，ΔT0=0.0 m，θ0=0.0°，φ0=-2.0°；其他運動參數(shù)的值取為0。其中T是靜水平衡狀態(tài)值，為一常數(shù)，φ0的初值不為0，相當(dāng)于一個初始的外部擾動作用。

圖3～圖7所示為所計算的各個浪向角的搖蕩幅頻曲線，其計算結(jié)果的變化規(guī)律與理論分析一致。圖3～圖5的結(jié)果顯示，垂蕩和縱搖對橫搖的耦合呈現(xiàn)出明顯的非線性效應(yīng)，圖8和圖9則表明，在橫搖共振區(qū)內(nèi)耦合還呈現(xiàn)出偏搖現(xiàn)象。

1）橫搖幅頻曲線計算分析。

由圖3～圖5中的橫搖幅頻曲線可以看出，波陡較高（Hw/λ=0.04）的共振區(qū)的幅頻曲線明顯低于波陡較低（Hw/λ=0.01）共振區(qū)的幅頻曲線，這表明了橫搖共振區(qū)幅頻曲線對波陡的非線性，即波高對橫搖幅頻曲線具有非線性效應(yīng)。由于計算中的阻尼系數(shù)按式（11）取常數(shù)值，因此可以認為該效應(yīng)是因垂蕩、縱搖耦合所導(dǎo)致，即垂蕩、縱搖對橫搖運動的耦合非線性效應(yīng)。

圖4和圖5顯示，高頻區(qū)（ωe≈1.0～1.2）的橫搖幅頻曲線有一個較小的峰值，該頻率范圍大致對應(yīng)垂蕩、縱搖的共振頻率，因此高頻小峰值表明了橫搖對垂蕩和縱搖耦合的響應(yīng)。

2）垂蕩、縱搖幅頻曲線計算分析。

圖6和圖7所示分別為頂浪和順浪時的垂蕩、縱搖幅頻曲線。圖4～圖7所示的垂蕩、縱搖幅頻曲線的計算結(jié)果顯示，不同波陡的垂蕩、縱搖幅頻曲線在所有頻率下基本重合，這說明波高對垂蕩、縱搖幅頻曲線不具有非線性效應(yīng)。但在正橫浪時，在橫搖共振區(qū)內(nèi)，不同波陡情況下的垂蕩、縱搖幅頻曲線具有一定程度的差異，該差異被認為是橫搖對垂蕩、縱搖運動耦合的非線性效應(yīng)。

圖4～圖7所示的縱搖幅頻曲線具有相同的響應(yīng)規(guī)律，即在共振頻率（ωe≈1.0～1.2）范圍內(nèi)有一個小的峰值響應(yīng)。斜浪（首、尾斜浪）的縱搖幅頻曲線與縱向波浪（頂浪和順浪）的縱搖幅頻曲線具有相似的響應(yīng)規(guī)律，表明橫搖對縱搖不存在耦合效應(yīng)。

垂蕩幅頻曲線有明顯的差異說明橫搖對垂蕩的耦合效應(yīng)顯著，該耦合效應(yīng)出現(xiàn)在頻率約大于橫搖共振頻率之后。縱向波浪中的穩(wěn)態(tài)橫搖趨近于0，也即圖6和圖7所示的垂蕩幅頻曲線不包含橫搖的效應(yīng)，對比圖4和圖5斜浪中響應(yīng)曲線所呈現(xiàn)出的差異，說明橫搖對垂蕩具有抑制的耦合作用。

3）橫搖偏搖現(xiàn)象。

圖8所示為波陡Hw/λ=0.04時的正橫浪搖蕩運動時歷曲線。橫搖角時歷曲線存在著明顯的左、右搖蕩幅值間的偏差，即在波浪傳播方向，搖蕩幅值大于其反向幅值，本文稱其為橫搖偏搖現(xiàn)象。橫搖偏搖現(xiàn)象也是垂蕩、縱搖耦合的一種效應(yīng)，計算表明，當(dāng)波高較大且在橫搖共振區(qū)內(nèi)時，會呈現(xiàn)出顯著的橫搖偏搖。圖9所示為Hw/λ=0.04時的橫搖偏搖幅值φa0幅頻曲線，其中正橫浪的橫搖偏搖幅值最大。

4）搖蕩耦合的非線性效應(yīng)機理初步分析。

上述搖蕩幅頻曲線計算分析表明，垂蕩、縱搖與橫搖的耦合會導(dǎo)致運動響應(yīng)的非線性效應(yīng)，其中，最為明顯的是垂蕩、縱搖對橫搖幅頻曲線的非線性效應(yīng)，以及橫搖偏搖現(xiàn)象。由于船體線型具有上、下和前、后不對稱的特點，所以搖蕩耦合的非線性效應(yīng)即來自于該船體線型的不對稱性。

就切片而言，波陡越大，平均水線以上和水線以下線型的不對稱性也越大，船體平均水線上、下線型的不對稱性是共振區(qū)橫搖幅頻曲線非線性的主要根源。

船體前、后線型的不對稱性是產(chǎn)生橫搖偏搖現(xiàn)象的原因。圖8所示為垂蕩（ΔT）、縱搖（θ）和橫搖（φ）的時歷計算曲線。由縱搖（θ）時歷曲線圖可以看出，縱搖運動呈現(xiàn)不均勻的搖蕩，埋艏（θ＜0）滯留時間較長，抬艏（θ＞0）滯留時間較短，此現(xiàn)象源自船體艏、艉線型的不對稱性，符合預(yù)期。對比橫搖與縱搖時歷曲線表面發(fā)現(xiàn)，二者基本是反向同步，縱搖的不均勻性傳遞給橫搖就會導(dǎo)致其左、右橫搖不對稱運動的出現(xiàn)，從而引起偏搖現(xiàn)象——左、右橫搖幅值不等。可以預(yù)見，前、后線型對稱的幾何船體將不會出現(xiàn)橫搖的偏搖現(xiàn)象。

4 結(jié) 論

本文基于廣義縱傾角θ和廣義吃水增量ΔT，以及船體坐標(biāo)系下的波面方程，采用切片理論方法建立了垂蕩、縱搖與橫搖耦合的動力學(xué)方程組；基于瞬時波面方程，為滿足波面壓力處為0的條件，合理修正了瞬時波面下壓力分布的計算公式，將靜水力與波浪擾動力之和FFK+B作為整體積分予以了計算。規(guī)則波浪中，艦船耦合搖蕩運動響應(yīng)數(shù)值計算結(jié)果表明：

1）采用耦合搖蕩切片計算方法可預(yù)測出垂蕩、縱搖對橫搖的耦合非線性效應(yīng)；

2）采用耦合搖蕩切片計算方法可預(yù)測出共振區(qū)內(nèi)的橫搖偏搖現(xiàn)象；

3）船體線型的上、下不對稱性是引起共振區(qū)橫搖非線性的主要原因；

4）船體線型的前、后不對稱性是引起共振區(qū)橫搖偏搖現(xiàn)象的主要原因。

本文提出了實用性的耦合搖蕩運動計算方法，開發(fā)了相應(yīng)的計算軟件，較好地計入了垂蕩、縱搖與橫搖之間的相互耦合效應(yīng)，對高海況下耐波性預(yù)報有積極的作用。后續(xù)，還將深入研究并完善計算方法。