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

李超凡， 朱仁傳， 周文俊

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

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

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

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

1.1 參考坐標系的建立

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

圖1 船體運動坐標系Fig.1 The coordinate system for ship motion

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

圖2 液艙內液體晃蕩坐標系示意Fig.2 The three ordinates describing sloshing in the liquid tank

1.2 基于脈沖響應函數法的頻-時域轉換

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

(1)

頻域的水動力系數中附加質量μij和阻尼系數λij在時域的轉換關系為：

(2)

(3)

則時延函數的表達式為：

(4)

1.3 非線性Froude-Krylov力和靜恢復力

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

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

(5)

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

(6)

(7)

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

圖3 瞬時濕表面船體網格Fig.3 Panels for transient wet surface

坐標轉換的方式為：

(8)

(9)

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

考慮液艙內部非定常的流體晃蕩問題，本文采用基于三維勢流理論的邊界元法對此進行時域數值分析。為減少計算的復雜度，液艙內部在動坐標系下求解，此時艙壁上的速度勢法向導數均為0，之后通過坐標轉換實現與船體運動的耦合。場內速度勢滿足的定解條件：

(10)

根據液艙尺寸和計算時裝載深度進行網格劃分實現模型的空間離散，場內任意一點的速度勢可根據格林定理得到，選取三維簡單格林函數，則速度勢可表達為：

(11)

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

(12)

(13)

設已知上一時刻的自由面速度勢與波面升高，根據此刻的速度勢φ和波面升高ζ，再經自由面條件計算此刻每個離散點的?φ/?t和?ζ/?t，實現時間步的遞進。

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

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

(14)

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

2 載液船舶與計算模型

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

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

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

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

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

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

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

3 耦合計算結果與分析

3.1 迎浪工況

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

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

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

已知任意變量可表示為：

n=1,2,…,n

(15)

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

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

(16)

(17)

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

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

表2 傅里葉級數擬合參數Table 2 Fourier series fitting parameters

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

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

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

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

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

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

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

3.2 橫浪工況

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

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

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

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

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

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

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

4 結論

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

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

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

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