波浪上船體與液艙晃蕩耦合的非線性時域模擬

李超凡， 朱仁傳， 周文俊

(1.上海交通大學(xué) 海洋工程國家重點實驗室， 上海 200240； 2.中國船舶及海洋工程設(shè)計研究院， 上海 200011)

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

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

1 載液船舶運動的基本理論

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

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

圖1 船體運動坐標(biāo)系Fig.1 The coordinate system for ship motion

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

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

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

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

(1)

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

(2)

(3)

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

(4)

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

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

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

(5)

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

(6)

(7)

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

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

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

(8)

(9)

1.4 邊界元法求解時域液艙晃蕩

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

(10)

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

(11)

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

(12)

(13)

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

1.5 基于時域混雜法的船舶時域運動方程

考慮弱散射的時域混雜法船舶時域運動的方程為：

(14)

圖4 非線性時域耦合運動計算流程Fig.4 Flowchart of non-linear prediction for coupled motion in time domain

2 載液船舶與計算模型

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

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

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

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

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

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

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

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

3.1 迎浪工況

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

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

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

已知任意變量可表示為：

n=1,2,…,n

(15)

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

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

(16)

(17)

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

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

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

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

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

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

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

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

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

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

3.2 橫浪工況

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

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

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

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

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

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

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

4 結(jié)論

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

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

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

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