李超凡， 朱仁傳， 周文俊

(1.上海交通大學(xué) 海洋工程國(guó)家重點(diǎn)實(shí)驗(yàn)室， 上海 200240； 2.中國(guó)船舶及海洋工程設(shè)計(jì)研究院， 上海 200011)

船舶行業(yè)的發(fā)展和研究與工程需要密切相關(guān)。如今LNG、LPG和一些特種液貨船舶等的使用比例日益增大，但該類船舶在航行狀態(tài)時(shí)劇烈的艙內(nèi)流體晃蕩卻對(duì)營(yíng)運(yùn)安全性提出考驗(yàn)，尤其當(dāng)船舶沒(méi)有裝滿時(shí)，船舶穩(wěn)性發(fā)生改變，艙內(nèi)液體受船體運(yùn)動(dòng)影響產(chǎn)生劇烈晃蕩，由此產(chǎn)生的晃蕩力與船體外部所受的波浪力相互耦合，給運(yùn)動(dòng)預(yù)報(bào)帶來(lái)困難。此外，對(duì)某些特種載液船型來(lái)說(shuō)，人們有時(shí)需要知道其在波浪中工作時(shí)某一時(shí)刻以后若干秒內(nèi)船體的運(yùn)動(dòng)幅值、位置、姿態(tài)等，對(duì)此展開更加準(zhǔn)確有效的運(yùn)動(dòng)時(shí)歷計(jì)算研究十分必要。Vassalos[1]、Zaraphonitiset[2]等較早地基于勢(shì)流理論模擬破損艙進(jìn)水問(wèn)題，但其忽略內(nèi)部液艙晃蕩的非線性影響。Santos等[3-4]研究了隨機(jī)波浪中破損船舶時(shí)域運(yùn)動(dòng)的數(shù)值方法，Newman[5]開發(fā)了WAMIT用于計(jì)算線性化的波浪中液艙晃蕩與船體耦合運(yùn)動(dòng)。Kim等[6]采用脈沖響應(yīng)函數(shù)法求解線性的船舶時(shí)域運(yùn)動(dòng)，對(duì)非線性液艙晃蕩采用有限差分法進(jìn)行模擬。Tang等[7]計(jì)算了時(shí)延函數(shù)對(duì)運(yùn)動(dòng)的影響。Spanos和Papanikolaou[8]，Liu和Papanikolaou[9]應(yīng)用頻域零航速自由面格林函數(shù)考慮航速修正，在瞬時(shí)濕表面上進(jìn)行壓力積分獲得非線性FK力和回復(fù)力進(jìn)行模擬。Ahmed、Hudson等[10]也采用了類似的方法。Zhao[11-12]和Hu等[13]基于勢(shì)流理論處理內(nèi)部晃蕩問(wèn)題，證實(shí)結(jié)果的可行性。Zou[14]等用時(shí)域勢(shì)流求解船體運(yùn)動(dòng)，商用CFD軟件處理內(nèi)部晃蕩。Huang[15]提出了一種基于線性理論的能量耗散條件，在非線性模型中考慮了晃蕩流的流動(dòng)粘性效應(yīng)。洪亮等[16]對(duì)船體內(nèi)外流場(chǎng)均采用時(shí)域勢(shì)流理論求解，將船體運(yùn)動(dòng)與液艙流體晃蕩構(gòu)建耦合運(yùn)動(dòng)方程。唐愷等[17]結(jié)合輻射和繞射問(wèn)題的線性解建立船舶非線性時(shí)域預(yù)報(bào)方法。李裕龍等[18-19]基于全時(shí)域勢(shì)流理論研究有航速時(shí)的液艙晃蕩耦合問(wèn)題。

本文基于勢(shì)流理論提出了一種計(jì)算在波浪中船體與液艙晃蕩耦合運(yùn)動(dòng)的非線性時(shí)域混雜法，重點(diǎn)考慮外流場(chǎng)非線性對(duì)液艙晃蕩耦合運(yùn)動(dòng)影響。文中內(nèi)部液艙晃蕩采用時(shí)域Rankine元法求解，將艙內(nèi)液體流動(dòng)引起的力矩進(jìn)行修正；對(duì)時(shí)域運(yùn)動(dòng)方程，輻射力與繞射力采用完整船舶計(jì)算的結(jié)果，對(duì)入射波力和靜恢復(fù)力考慮瞬時(shí)濕表面的非線性影響；將液艙晃蕩力矩和非線性入射波力、靜恢復(fù)力與繞射波浪激勵(lì)力放在時(shí)域運(yùn)動(dòng)方程的右端建立耦合運(yùn)動(dòng)方程，實(shí)現(xiàn)對(duì)規(guī)則波或不同工況下載液船舶運(yùn)動(dòng)時(shí)歷的求解。本文中有航速的工況采用移動(dòng)脈動(dòng)源方法計(jì)算頻域結(jié)果，計(jì)算分析對(duì)象為一載液S175船，對(duì)其進(jìn)行有無(wú)液艙的線性時(shí)域模擬與非線性時(shí)域模擬，給出運(yùn)動(dòng)幅值響應(yīng)算子RAO曲線，并與實(shí)驗(yàn)值進(jìn)行對(duì)比驗(yàn)證。

1 載液船舶運(yùn)動(dòng)的基本理論

1.1 參考坐標(biāo)系的建立

大地坐標(biāo)系為o0x0y0z0，隨船舶一起以相同速度運(yùn)動(dòng)而構(gòu)成表征船舶搖蕩位移和姿態(tài)的基準(zhǔn)的坐標(biāo)系為參考坐標(biāo)系oxyz，該坐標(biāo)系的oxy平面與靜水面重合。與船體固結(jié)，隨船體一起搖蕩的坐標(biāo)系為動(dòng)坐標(biāo)系o′x′y′z′。各坐標(biāo)系如圖1所示。

圖1 船體運(yùn)動(dòng)坐標(biāo)系Fig.1 The coordinate system for ship motion

在液艙內(nèi)部由于考慮流體晃蕩問(wèn)題，對(duì)其采用三維時(shí)域勢(shì)流理論求解，為此需提前確定時(shí)域勢(shì)流計(jì)算中的邊界條件。若在大地坐標(biāo)系或參考坐標(biāo)系中計(jì)算，則自由面條件的表達(dá)較為簡(jiǎn)單，但艙壁面邊界條件十分復(fù)雜。因此，本文中液艙先在oxy面與艙內(nèi)平均液面重合的動(dòng)坐標(biāo)系中求解，此時(shí)艙壁上速度勢(shì)的法向偏導(dǎo)數(shù)均為0，使外界激勵(lì)的影響體現(xiàn)在網(wǎng)格數(shù)更少的自由面的條件中，將液艙晃蕩在動(dòng)坐標(biāo)系中產(chǎn)生的力矩模擬計(jì)算后，再通過(guò)運(yùn)動(dòng)和力在坐標(biāo)系中的轉(zhuǎn)換實(shí)現(xiàn)船體運(yùn)動(dòng)與液艙晃蕩的耦合計(jì)算。液艙計(jì)算所在坐標(biāo)系如圖2所示。

圖2 液艙內(nèi)液體晃蕩坐標(biāo)系示意Fig.2 The three ordinates describing sloshing in the liquid tank

1.2 基于脈沖響應(yīng)函數(shù)法的頻-時(shí)域轉(zhuǎn)換

本文輻射力部分采用頻域轉(zhuǎn)時(shí)域方式求得[20]。在時(shí)域計(jì)算時(shí)采用卡明斯的脈沖響應(yīng)法思想，把其時(shí)間歷程看成一系列瞬時(shí)的脈沖運(yùn)動(dòng)組成，由此建立的船體運(yùn)動(dòng)微分方程可表達(dá)為：

(1)

頻域的水動(dòng)力系數(shù)中附加質(zhì)量μij和阻尼系數(shù)λij在時(shí)域的轉(zhuǎn)換關(guān)系為：

(2)

(3)

則時(shí)延函數(shù)的表達(dá)式為：

(4)

1.3 非線性Froude-Krylov力和靜恢復(fù)力

在時(shí)域理論角度，其相對(duì)頻域的優(yōu)勢(shì)還體現(xiàn)在能夠反映瞬變或非線性結(jié)果的影響，由于船舶存在內(nèi)部液艙晃蕩問(wèn)題，在計(jì)算時(shí)考慮一定的非線性影響十分必要。本文在計(jì)算時(shí)通過(guò)瞬時(shí)濕表面上的流體壓力積分得到非線性入射波力和靜回復(fù)力，波幅為ζa的簡(jiǎn)諧規(guī)則波的入射勢(shì)為：

i(ωt-k(xcosβ+ysinβ))]

(5)

非線性回復(fù)力根據(jù)其物理含義直接計(jì)算，船舶靜水回復(fù)力等于實(shí)際濕表面的流體靜水力與船舶正浮平均濕表面流體靜力之差：

(6)

(7)

以上計(jì)算需求取瞬時(shí)濕表面，為此需要先對(duì)船體網(wǎng)格進(jìn)行坐標(biāo)轉(zhuǎn)換。由于運(yùn)動(dòng)是動(dòng)坐標(biāo)系相對(duì)參考坐標(biāo)的，而船體網(wǎng)格在動(dòng)坐標(biāo)系下表示，所以先將網(wǎng)格點(diǎn)坐標(biāo)位置轉(zhuǎn)移到參考坐標(biāo)系下便于用波面來(lái)截船體表面。圖3為瞬時(shí)濕表面船體網(wǎng)格示意圖。

圖3 瞬時(shí)濕表面船體網(wǎng)格Fig.3 Panels for transient wet surface

坐標(biāo)轉(zhuǎn)換的方式為：

(8)

(9)

1.4 邊界元法求解時(shí)域液艙晃蕩

考慮液艙內(nèi)部非定常的流體晃蕩問(wèn)題，本文采用基于三維勢(shì)流理論的邊界元法對(duì)此進(jìn)行時(shí)域數(shù)值分析。為減少計(jì)算的復(fù)雜度，液艙內(nèi)部在動(dòng)坐標(biāo)系下求解，此時(shí)艙壁上的速度勢(shì)法向?qū)?shù)均為0，之后通過(guò)坐標(biāo)轉(zhuǎn)換實(shí)現(xiàn)與船體運(yùn)動(dòng)的耦合。場(chǎng)內(nèi)速度勢(shì)滿足的定解條件：

(10)

根據(jù)液艙尺寸和計(jì)算時(shí)裝載深度進(jìn)行網(wǎng)格劃分實(shí)現(xiàn)模型的空間離散，場(chǎng)內(nèi)任意一點(diǎn)的速度勢(shì)可根據(jù)格林定理得到，選取三維簡(jiǎn)單格林函數(shù)，則速度勢(shì)可表達(dá)為：

(11)

時(shí)間離散方式采用中心差分格式，計(jì)算表達(dá)，

(12)

(13)

設(shè)已知上一時(shí)刻的自由面速度勢(shì)與波面升高，根據(jù)此刻的速度勢(shì)φ和波面升高ζ，再經(jīng)自由面條件計(jì)算此刻每個(gè)離散點(diǎn)的?φ/?t和?ζ/?t，實(shí)現(xiàn)時(shí)間步的遞進(jìn)。

1.5 基于時(shí)域混雜法的船舶時(shí)域運(yùn)動(dòng)方程

考慮弱散射的時(shí)域混雜法船舶時(shí)域運(yùn)動(dòng)的方程為：

(14)

圖4 非線性時(shí)域耦合運(yùn)動(dòng)計(jì)算流程Fig.4 Flowchart of non-linear prediction for coupled motion in time domain

2 載液船舶與計(jì)算模型

本文選取一加載方形液艙的S175高速集裝箱船作為算例模型[14]，船體橫剖線圖如圖5所示。該模型試驗(yàn)由鄒康等在中國(guó)船舶科學(xué)研究中心耐波性水池中完成，實(shí)船與模型的主尺度對(duì)比如表1所示，縮尺比為1∶55。

圖5 S175船型線圖Fig.5 Body plan of S175

表1 S175實(shí)船與模型主尺度Table 1 Principal particulars of S175 and its model

液艙位于第9站至第13站之間，關(guān)于船舶中縱剖面對(duì)稱，重心位置與船模重心重合，如圖6所示。該液艙大小為0.6 m×0.3 m×0.25 m，此時(shí)艙內(nèi)裝載液體深度為0.125 m。

圖6 液艙大小尺寸與位置Fig.6 The size and location of tank

由于考慮了非線性入射波力和非線性靜恢復(fù)力的影響，船體網(wǎng)格需劃分到水線以上一段距離，在gambit中建立四節(jié)點(diǎn)船體網(wǎng)格，液艙網(wǎng)格采用三角形面元?jiǎng)澐?，如圖7所示。

圖7 船體瞬時(shí)濕表面網(wǎng)格與液艙網(wǎng)格示意Fig.7 Panels for ship and tank

3 耦合計(jì)算結(jié)果與分析

3.1 迎浪工況

圖8為零航速迎浪不同頻率下有無(wú)液艙的縱搖和垂蕩運(yùn)動(dòng)時(shí)歷曲線，從考慮非線性因素的計(jì)算結(jié)果來(lái)看，液艙晃蕩力矩對(duì)縱搖運(yùn)動(dòng)的影響較小，隨著頻率的減小這種差異也逐漸變小。有無(wú)液艙的垂蕩運(yùn)動(dòng)時(shí)歷峰值存在一個(gè)小的差值，但整體零航速迎浪的運(yùn)動(dòng)所受影響不大。

圖8 零航速迎浪有無(wú)液艙縱搖和垂蕩運(yùn)動(dòng)時(shí)歷結(jié)果Fig.8 The pitch and heave motions in head sea when Fn=0

對(duì)于數(shù)值計(jì)算得到的船體運(yùn)動(dòng)時(shí)歷，需要進(jìn)行頻譜分析結(jié)合傅里葉級(jí)數(shù)展開法對(duì)各時(shí)歷曲線進(jìn)行擬合，由此將時(shí)域計(jì)算結(jié)果轉(zhuǎn)換到頻域進(jìn)行試驗(yàn)值對(duì)比。

已知任意變量可表示為：

n=1,2,…,n

(15)

以迎浪工況波長(zhǎng)1.5 m，波幅0.1 m有液艙時(shí)船舶縱搖運(yùn)動(dòng)時(shí)歷為例。待運(yùn)動(dòng)穩(wěn)定后，舍去前段選取之后的若干個(gè)周期結(jié)果進(jìn)行傅里葉級(jí)數(shù)擬合，X5可展開為：

X5=X5t0+X5t1cos(ωet+γ1)+X5t2cos(ωet+γ2)…

(16)

(17)

式中：X5為縱搖運(yùn)動(dòng)結(jié)果；ωe為遭遇頻率；X5t0為傅里葉級(jí)數(shù)展開的定常項(xiàng)即平均值；X5t1為一階運(yùn)動(dòng)響應(yīng)幅值；γ1為一階運(yùn)動(dòng)響應(yīng)相位，以此類推。

將時(shí)域運(yùn)動(dòng)通過(guò)傅里葉變換取到3階精度，相關(guān)參數(shù)如表2所示。從擬合結(jié)果可以分析認(rèn)為除主頻率外其余噪聲受非線性計(jì)算影響，其結(jié)果相對(duì)線性得到修正。擬合曲線如圖9所示。

表2 傅里葉級(jí)數(shù)擬合參數(shù)Table 2 Fourier series fitting parameters

圖9 傅里葉級(jí)數(shù)展開示例Fig.9 Example of Fourier series expansion method

圖10將考慮非線性影響的有無(wú)液艙運(yùn)動(dòng)與線性時(shí)域程序計(jì)算結(jié)果和船模試驗(yàn)值進(jìn)行比較，可以發(fā)現(xiàn)線性時(shí)域計(jì)算結(jié)果在中高頻率不及本文所用方法的結(jié)果，再結(jié)合試驗(yàn)值可以看出考慮非線性影響后對(duì)縱搖的準(zhǔn)確預(yù)報(bào)有了較為明顯的改善。

圖10 零航速迎浪下非線性和線性數(shù)值計(jì)算結(jié)果與實(shí)驗(yàn)值的縱搖RAO對(duì)比Fig.10 Comparison of pitch RAO in head sea by experiment and linear and nonlinear calculation when Fn=0

圖11～12中給出傅汝德數(shù)為0.275時(shí)迎浪不同頻率下經(jīng)傅里葉變換得到的有無(wú)液艙的縱搖和垂蕩運(yùn)動(dòng)RAO曲線，試驗(yàn)值為ITTC在1978年給出的S175在Fn=0.275時(shí)無(wú)液艙的RAO，線性數(shù)值曲線為不加載液艙的計(jì)算結(jié)果。

圖11 Fn=0.275迎浪下非線性和線性數(shù)值計(jì)算結(jié)果與實(shí)驗(yàn)值的縱搖RAO對(duì)比Fig.11 Comparison of pitch RAO in head sea by experiment and linear and nonlinear calculation when Fn=0.275

從圖11、12中可以看出，在考慮回復(fù)力和入射波力的非線性情況下，與線性處理結(jié)果對(duì)比均能得到比較滿意的縱搖運(yùn)動(dòng)預(yù)報(bào)，并且由于液艙內(nèi)流體晃蕩和波浪誘導(dǎo)力的相互作用，在某些頻率范圍內(nèi)相互抵消減小了船體的運(yùn)動(dòng)幅值。垂蕩情況下線性時(shí)域方法所得的峰值結(jié)果較大，低頻時(shí)裝載后的非線性方法計(jì)算所得結(jié)果較未加載狀態(tài)時(shí)小，當(dāng)頻率增大到一定范圍時(shí)，裝載液艙后的運(yùn)動(dòng)幅值又較無(wú)載液時(shí)大，其原因在于液艙誘導(dǎo)力矩與波浪誘導(dǎo)力矩疊加造成。

圖12 Fn=0.275迎浪下非線性和線性數(shù)值計(jì)算結(jié)果與實(shí)驗(yàn)值的垂蕩RAO對(duì)比Fig.12 Comparison of heave RAO in head sea by experiment and linear and nonlinear calculation when Fn=0.275

3.2 橫浪工況

對(duì)零航速橫浪時(shí)有無(wú)液艙的S175在線性與非線性的情況下進(jìn)行模擬，圖13中給出零航速橫浪入射波長(zhǎng)為9.1、20、40.1 m時(shí)有無(wú)液艙的橫搖和垂蕩運(yùn)動(dòng)時(shí)歷曲線。從橫搖運(yùn)動(dòng)可以看出，在高頻時(shí)，運(yùn)動(dòng)幅值變化不大，非線性的模擬結(jié)果較線性有明顯改善；在接近共振頻率附近時(shí)運(yùn)動(dòng)幅值明顯減小，原因與液艙晃蕩力與波浪力在相位和量級(jí)上的變化有關(guān)；當(dāng)頻率繼續(xù)逐漸降低時(shí)，有液艙的橫搖幅值開始增大，船舶的共振頻率發(fā)生改變。對(duì)垂蕩運(yùn)動(dòng)來(lái)說(shuō)其運(yùn)動(dòng)幅值的變化不大，可見液艙晃蕩力矩對(duì)垂直方向的影響較小。

圖13 零航速橫浪有無(wú)液艙橫搖和垂蕩時(shí)歷Fig.13 The roll and heave motions in beam sea when Fn=0

圖14分別給出了不同波長(zhǎng)的入射波激勵(lì)下，在液艙晃蕩與船體耦合運(yùn)動(dòng)過(guò)程中，液艙誘導(dǎo)力矩與波浪力的時(shí)歷對(duì)比。其對(duì)應(yīng)的無(wú)因次化頻率分別為1.5、1、0.7。

在短波時(shí)液艙晃蕩力與波浪激勵(lì)力量級(jí)相差較大，此時(shí)液艙晃蕩對(duì)船體運(yùn)動(dòng)造成的影響很小，非線性與線性結(jié)果差異不明顯。當(dāng)頻率位于不加載液艙船舶的共振頻率附近時(shí)，液艙晃蕩力約在波浪力幅值的1/10～1/3內(nèi)，此時(shí)對(duì)船體運(yùn)動(dòng)影響較大，當(dāng)波浪為20 m時(shí)，可以看出液艙晃蕩力的相位與波浪激勵(lì)力存在約180°的相位差，此時(shí)二者疊加船舶所受合外力的幅值減小，而船舶橫搖運(yùn)動(dòng)的幅值也有了明顯的減小。當(dāng)波長(zhǎng)逐漸增加，到達(dá)載液船舶的共振頻率附近時(shí)，如圖14(c)，可以看出此時(shí)波浪液艙晃蕩力與波浪激勵(lì)力仍處于一個(gè)量級(jí)，且相位發(fā)生改變，線性與非線性外流場(chǎng)運(yùn)動(dòng)傳遞后造成的液艙晃蕩力基本一致，平穩(wěn)后與波浪激勵(lì)力可見峰峰疊加的現(xiàn)象，船體所受總的力增加，此時(shí)船舶橫搖運(yùn)動(dòng)更為劇烈，對(duì)應(yīng)頻率下船舶橫搖的 RAO 也有所增加，船舶容易處于危險(xiǎn)狀態(tài)。

圖14 橫浪工況橫搖波浪力與液艙晃蕩力比較Fig.14 Comparison between sloshing force and wave force in beam sea

圖15給出計(jì)算多個(gè)入射波頻率后經(jīng)傅里葉變換得到的有無(wú)液艙船舶在線性和非線性下的橫搖運(yùn)動(dòng)RAO曲線。液艙流體晃蕩使得S175 船模的橫搖共振頻率區(qū)間發(fā)生了偏移，船舶加載液艙后橫搖慣性半徑變大，重心位置和初穩(wěn)性高也發(fā)生改變，橫搖固有周期增大，相應(yīng)的共振頻率減小，這與模擬結(jié)果一致。線性計(jì)算與非線性差距不大，其橫搖運(yùn)動(dòng)峰值主要與橫搖的一次和二次阻尼系數(shù)相關(guān)。船舶進(jìn)水后，非線性計(jì)算結(jié)果相較線性結(jié)果有一定的改善，與圓形散點(diǎn)的實(shí)驗(yàn)值更為接近，但在無(wú)液艙的共振頻率附近的模擬仍不是很理想。

圖15 橫浪下非線性和線性有無(wú)液艙數(shù)值計(jì)算結(jié)果與實(shí)驗(yàn)值的橫搖RAO對(duì)比Fig.15 Comparison of roll RAO in beam sea by experiment and linear and nonlinear calculation

4 結(jié)論

1)模型慮及瞬時(shí)濕表面影響下的入射力與回復(fù)力的非線性，內(nèi)部流動(dòng)在液艙坐標(biāo)系下求解，提高了液艙晃蕩耦合運(yùn)動(dòng)模擬的精度，證實(shí)考慮非線性后的結(jié)果能得到更為準(zhǔn)確的模擬。

2)迎浪有無(wú)液艙的縱搖和垂蕩運(yùn)動(dòng)隨航速的變換趨勢(shì)大致相同，隨頻率的增大運(yùn)動(dòng)幅值差距較小。零航速橫浪時(shí)的橫搖運(yùn)動(dòng)在液艙作用下，船體橫搖幅值減小，峰值位置前移。

3)耦合運(yùn)動(dòng)RAO的變化與液艙晃蕩力和波浪誘導(dǎo)力的相位、量級(jí)改變有關(guān)，當(dāng)兩者相差不大，相位相差180°左右時(shí)，橫搖運(yùn)動(dòng)的劇烈程度得到改善。

4)波浪上船舶液艙晃蕩的時(shí)域運(yùn)動(dòng)預(yù)報(bào)在線性條件下有一定的局限性，但完全非線性的模擬計(jì)算代價(jià)巨大，在工程應(yīng)用上還有一定的距離，而本文所用方法計(jì)算快速高效，在幫助處理載液船舶的檢驗(yàn)或前期設(shè)計(jì)中能提供有力幫助。