張 珍,陳明輝,邵 菲,謝永和?

(1.浙江海洋大學(xué) 船舶與海運(yùn)學(xué)院,浙江 舟山316022;2.中國艦船研究設(shè)計(jì)中心,湖北 武漢430064)

載液貨船在航行過程中引起的液艙晃蕩問題一直以來都備受關(guān)注。強(qiáng)非線性以及隨機(jī)性使得液艙晃蕩一直是流體力學(xué)領(lǐng)域極具挑戰(zhàn)的研究課題之一。對液艙晃蕩引起的晃蕩載荷、結(jié)構(gòu)響應(yīng)以及兩者耦合效應(yīng)的研究是保證船舶航行的重要依據(jù),也為載液船液艙設(shè)計(jì)提供了參考價(jià)值。

目前液艙晃蕩的研究主要集中在理論分析,試驗(yàn)研究以及數(shù)值模擬3 種途徑。最早源于ABRAMSON[1]對線性勢流理論方面的總結(jié),隨后FALTINSEN[2]采用非線性勢流理論對水平激勵(lì)作用下的矩形液艙晃蕩進(jìn)行了探究。相比理論分析,試驗(yàn)探究可以為液艙晃蕩提供直觀借鑒,FABRIZIO PISTANI[3]開展晃蕩模型試驗(yàn),通過對比不同箱體晃蕩沖擊壓力,發(fā)現(xiàn)棱形箱體斜面設(shè)計(jì)可有效減緩沖擊壓力。BRAEUNIG 等[4]對橫蕩激勵(lì)下LNG 液艙內(nèi)的晃蕩現(xiàn)象進(jìn)行了試驗(yàn)和分析。雖然模型試驗(yàn)可以直觀展現(xiàn)晃蕩劇烈程度,但在部分條件限制下很難完成,且花費(fèi)的成本較高。隨著社會的進(jìn)步和計(jì)算機(jī)的廣泛運(yùn)用,數(shù)值模擬方法備受青睞。其中有基于自由表面跟蹤技術(shù)的探究,如:BATTAGLIA[5]通過水平(Level Set)法追蹤自由液面對三維矩形液艙晃蕩進(jìn)行探究。SHAMSODDINI[6]基于ISPH 法分析了擋板對于減輕液艙晃蕩的影響。還有根據(jù)計(jì)算網(wǎng)格的研究,如:JIANG[7]等基于Open FOAM 模擬棱形液艙晃蕩,ZHAO[8]等基于重疊網(wǎng)格法對部分裝載的LNG 液艙進(jìn)行了研究。

近年來,液艙晃蕩問題研究越來越廣泛和趨于成熟,然而,大部分學(xué)者是將液艙剛性化,但在實(shí)際問題中,液艙其實(shí)是彈性結(jié)構(gòu),導(dǎo)致晃蕩引起的沖擊壓力會受到艙壁彈性的影響。隨著水彈性力學(xué)的發(fā)展,陸續(xù)有學(xué)者開展了彈性液艙晃蕩的研究。朱仁慶和劉艷敏[9]將有限體積法、顯示積分法、ALE 法相結(jié)合對彈性液艙晃蕩進(jìn)行數(shù)值模擬探究,不僅考慮了外部激勵(lì),還考慮了裝載率變化對彈性液艙內(nèi)晃蕩載荷與響應(yīng)的影響。ZHANG[10]采用MPSFEM 耦合方法對晃蕩沖擊載荷引起的彈性結(jié)構(gòu)響應(yīng)進(jìn)行數(shù)值分析,改變材料的楊氏模量探究彈性艙內(nèi)晃蕩現(xiàn)象,分析沖擊壓力變化和結(jié)構(gòu)在時(shí)域頻域的動(dòng)態(tài)響應(yīng)特征。蔣梅榮等[11]基于振動(dòng)臺試驗(yàn)對彈性壁液艙的共振情況做了探究,將液艙置于簡諧激勵(lì)振動(dòng)下,研究液艙內(nèi)液體晃動(dòng)現(xiàn)象。

本工作基于數(shù)值模擬方法對剛性和彈性液艙晃蕩問題進(jìn)行充分探究。通過改變外部激勵(lì)頻率參數(shù),模擬剛性和彈性液艙晃蕩過程,對比分析因晃蕩而產(chǎn)生的對艙壁沖擊壓力以及現(xiàn)象圖,探究外部激勵(lì)對剛性與彈性液艙晃蕩的影響。

1 理論方法

1.1 控制方程

本文主要是求解N-S方程,計(jì)算隱式不定常以及不可壓縮的粘性流體在流場中的一系列問題。N-S方程如下所示:

1.2 自由液面捕捉方法

本工作應(yīng)用流體體積(VOF)方法[12],通過追蹤每一個(gè)網(wǎng)格中水和體積的體積分?jǐn)?shù)來捕捉兩相流的相界面。該方法假設(shè)網(wǎng)格分辨率足以求解相之間的交界面位置和形狀。

交界面的相分布和位置由相體積分?jǐn)?shù)的場來描述,根據(jù)相體積分?jǐn)?shù)的值,可以區(qū)分網(wǎng)格單元中不同相或流體是否存在,相φ的體積分?jǐn)?shù)定義如下:

其中,w和a分別表示水和空氣,本文采用φ＝0.5作為水和空氣間的相界面。應(yīng)用VOF可以實(shí)現(xiàn)良好的質(zhì)量守恒以及提高計(jì)算效率。

1.3 FVM 方 法

本工作采用的數(shù)值模擬核心方法是有限體積法(FVM)[13],FVM 法需要離散化求解區(qū)域,并將其劃分為有限大小的離散網(wǎng)格。每個(gè)網(wǎng)格節(jié)點(diǎn)以特定方法形成圍繞該節(jié)點(diǎn)的控制體積。因此,其核心是將控制微分方程集成到控制體積中:

其中,等式左邊第一項(xiàng)表示控制變量V中特征變量φ的總量隨時(shí)間變化。n·(ρφu)是沿n方向的流率,n是由于對流流動(dòng)而在控制體積表面外部的法線方向。n·(Γ·gradφ)表示φ向外部控制量的擴(kuò)散率。

1.4 流固耦合(FSI)

流固耦合[14]用來描述移動(dòng)或變形的固體與流體內(nèi)流場或外流場之間相互作用、相互影響的問題。當(dāng)流動(dòng)的流體與固體結(jié)構(gòu)接觸時(shí),固體會受到應(yīng)力與應(yīng)變的作用,這些力會使結(jié)構(gòu)產(chǎn)生變形。應(yīng)力使作用于表面的單位面積力的度量。

Von Mises應(yīng)力[15]是很重要的韌性失效度量,表示偏應(yīng)力的大小,與平均應(yīng)力無關(guān):

如果材料處于單軸應(yīng)力狀態(tài),則Von Mises應(yīng)力等于應(yīng)力的非零分量。通常,通過加載大致處于單軸應(yīng)力狀態(tài)的試樣來研究材料屈服。式(9)中,S1,S2,S3分別為第一,二,三主應(yīng)力。其大概的含義是在一定的變形條件下,當(dāng)材料的單位體積形狀改變的彈性位能(又稱彈性形變能)達(dá)到某一常數(shù)時(shí),材料就屈服;即當(dāng)單元體的形狀改變比能達(dá)到一定程度,材料開始屈服。

2 模型驗(yàn)證

為驗(yàn)證數(shù)值模擬方法的可靠性,采用西班牙造船廠Navantia建造的138 000 m2的LNG 液艙縱剖面做橫搖模擬,該剖面按照1∶50的比例模型化(如圖1),液艙的長度為L＝0.90 m,液艙的高度為H＝0.58 m,液艙內(nèi)部的液體高度h＝0.093 m。在圖(1)中,位于液艙左側(cè)固壁上位于自由液面處的一點(diǎn)P點(diǎn)為壓力監(jiān)測點(diǎn)。該液艙模型的運(yùn)動(dòng)方程表示為

圖1 液倉物理模型Fig.1 Physical model of tank

其中,θmax＝4°,為角位移振幅,周期T＝2.112 s,此為液艙的自然共振周期。

液艙橫搖運(yùn)動(dòng)的中心選在液艙底部的中心位置,檢測位于液艙左側(cè)固壁上自由液面處一點(diǎn)P的壓力隨液艙橫搖運(yùn)動(dòng)的變化。為了便于同實(shí)驗(yàn)值進(jìn)行比較,對時(shí)間和壓力進(jìn)行了無量綱化處理,壓力無量綱化處理公式為P＝P/ρgh,其中,P為P點(diǎn)處的壓力值,ρ為液艙內(nèi)液體的密度,h為液艙內(nèi)的液體高度;g為重力加速度。

由于本工作采用的數(shù)值方法需要進(jìn)行網(wǎng)格劃分,因此在計(jì)算前進(jìn)行了網(wǎng)格獨(dú)立性驗(yàn)證,從而來確保精度。圖2列出了15 s內(nèi)最大壓力隨網(wǎng)格數(shù)變化的情況。采用了4種網(wǎng)格數(shù)來進(jìn)行驗(yàn)證,從圖中可以看出,最大誤差在3%左右,因而可以驗(yàn)證網(wǎng)格的無關(guān)性。

圖2 網(wǎng)格無關(guān)性驗(yàn)證Fig.2 Grid independence verification

為驗(yàn)證時(shí)間步長的獨(dú)立性,對不同時(shí)間步長的壓力歷時(shí)曲線進(jìn)行了對比探究。圖3分別顯示了壓力在步長為0.001,0.001 5和0.002 s的歷時(shí)曲線。從圖中可見,最大誤差在0.28%左右,因此,驗(yàn)證步長的獨(dú)立性。而本文采取0.001 s的時(shí)間步長來更精準(zhǔn)探究液艙晃蕩的影響。

圖3 P 點(diǎn)壓力歷時(shí)曲線Fig.3 Time history of pressure at the point P

圖4展現(xiàn)了液艙縱剖面橫搖晃蕩前十五個(gè)周期的角度和液艙左側(cè)壁面上P點(diǎn)處的壓力時(shí)歷值,包括實(shí)驗(yàn)結(jié)果與數(shù)值模擬結(jié)果。從圖中可以明顯看出液艙左側(cè)固壁上P點(diǎn)的壓力在自由表面波拍擊壁面時(shí)刻發(fā)生徒增現(xiàn)象。壓力最大值約為1 kPa,壓力值誤差較小,且數(shù)值模擬曲線波動(dòng)基本與試驗(yàn)曲線一致,結(jié)果吻合的很好。從而驗(yàn)證了本工作采用的數(shù)值模擬方法在解決液艙晃蕩問題方面是相當(dāng)有效的。

圖4 P 點(diǎn)的壓力歷時(shí)曲線Fig.4 Time history of pressure at the point P

3 數(shù)值模擬

本工作采用液艙幾何尺寸如圖5所示,液艙長為1.2 m,高為0.6 m,艙內(nèi)水深0.12 m,對應(yīng)充水率為20%,右側(cè)和底部艙壁共設(shè)置3個(gè)測壓點(diǎn),位置如圖1所示。液艙沿x方向做如下運(yùn)動(dòng):其中,A＝0.06 m,為液艙運(yùn)動(dòng)幅值,由解析理論得到液艙晃蕩一階共振頻率ω＝2.975 s－1,為驗(yàn)證外部激勵(lì)頻率對液艙晃蕩的影響,外部激勵(lì)頻率設(shè)置如表1所示的3種工況。f為此裝載下液艙晃蕩的共振頻率。

圖5 液艙幾何模型Fig.5 Geometric model of the tank

表1 計(jì)算工況Table 1 Cases for simulation

3.1 剛性液艙晃蕩

為了計(jì)算方便,本節(jié)將液艙簡化為二維結(jié)構(gòu)。將艙壁視為剛性。網(wǎng)格采用四面體棱柱層網(wǎng)格。圖6給出了不同激勵(lì)頻率下晃蕩流體的運(yùn)動(dòng)情況。從圖6可以看出,激勵(lì)頻率對晃蕩流體晃蕩特征具有顯著的影響。隨著激勵(lì)頻率變大,液體沖擊左側(cè)艙壁的高度越高,出現(xiàn)碎波現(xiàn)象越來越嚴(yán)重。從圖6(h)可見,液面卷入空氣比另兩種工況嚴(yán)重,而到圖6(i)所示時(shí)刻時(shí),出現(xiàn)了明顯不同,且伴有液體飛濺現(xiàn)象,這種飛濺現(xiàn)象因激勵(lì)頻率變大,而更嚴(yán)重。

圖6 晃蕩過程現(xiàn)象圖Fig.6 Impact process of the tank sloshing

圖7～9為3個(gè)壓力點(diǎn)所監(jiān)測的壓力變化情況。圖7可見,P1壓力點(diǎn)所監(jiān)測壓力差異較小。在7 s前,3種工況曲線基本重合,7 s后,呈現(xiàn)的雙峰特性的第二個(gè)峰值附近有些許差異。在12.5 s左右,0.4f狀態(tài)下,曲線出現(xiàn)明顯波動(dòng),而波動(dòng)未在另兩種工況下發(fā)生。在14 s時(shí),2倍的激勵(lì)頻率作用下,產(chǎn)生的壓力比另兩種工況產(chǎn)生的壓力值要大,約為1 k Pa,另兩種工況值則低于1 kPa。0.4f工況下,壓力最大值約為2 000 Pa,而2f工況下,最大壓力達(dá)到近2 250 Pa。

圖7 P 1 壓力點(diǎn)壓力變化歷史曲線Fig.7 Time history of pressure at the point P 1

圖8為P2壓力點(diǎn)所監(jiān)測壓力的對比圖。0.4f工況下,最大壓力值約1 750 Pa,而2f工況下,最大值約為1 950 Pa。14 s時(shí),曲線有明顯差異,2f工況下壓力大于其他兩種工況。壓力曲線仍然呈雙峰特性。

圖8 P 2 壓力點(diǎn)壓力變化歷史曲線Fig.8 Time history of pressure at the point P 2

圖9為P3壓力所監(jiān)測壓力對比圖,因?yàn)楸狙芯可婕皟上嗔?所以壓力值出現(xiàn)了負(fù)值的情況。從圖中可見,2f工況下,液面卷入空氣情況較嚴(yán)重。作用P3壓力點(diǎn)的最大壓力為f工況下,達(dá)到1 220 Pa,而2f工況下最大壓力低于1 200 Pa。在11.3 s左右,f工況下壓力比另兩種工況下壓力要大,但在13.6 s左右,2f工況下最大壓力大于其余兩種工況。此外,2f工況下,曲線波動(dòng)沒有另兩種工況平滑。

圖9 P 3 壓力點(diǎn)壓力變化歷史曲線Fig.9 Time history of pressure at the point P 3

綜合圖7～9可見,外部激勵(lì)頻率對液艙晃蕩有一定的影響,頻率越大,晃蕩越明顯,自由液面出現(xiàn)碎波,液體飛濺現(xiàn)象越嚴(yán)重。在7 s前,3種工況壓力歷時(shí)曲線基本吻合,但在7 s后,曲線開始發(fā)生變化,2f工況下,曲線的波動(dòng)次數(shù)比其他兩種工況要多,導(dǎo)致曲線不如另兩種工況平滑。在14 s附近,P1,P2點(diǎn)壓力差異較明顯。13.5 s左右,P3點(diǎn)的峰值差異較明顯。由于兩相流的作用,P3值出現(xiàn)了負(fù)值情況,這是由于空氣對其影響而產(chǎn)生的。

3.2 彈性液艙晃蕩

為探究外部激勵(lì)頻率對彈性液艙晃蕩的影響,殼體單元采用鋁材質(zhì),因?yàn)殇X的特性相對較柔,可以達(dá)到模擬的目的。液艙艙壁厚度為10 mm,具體材料參數(shù)如表2所示。圖10為20%充液率狀態(tài)下三維液艙模型圖。因?yàn)榱鞴恬詈鲜欠浅?fù)雜的問題,所以在對模型進(jìn)行網(wǎng)格劃分時(shí),固體域和流體域采用不同網(wǎng)格:固體域部分采用四面體網(wǎng)格,流體域部分采用多面體網(wǎng)格。圖11為模型的網(wǎng)格劃分圖。

圖10 液艙計(jì)算模型Fig.10 Tanks calculation model

圖11 模型網(wǎng)格劃分圖Fig.11 Mesh used in the simulation

表2 材料參數(shù)表Table 2 Material parameters

圖12為考慮艙壁的彈性影響,在不同激勵(lì)頻率作用下,P1點(diǎn)和P2的壓力歷時(shí)曲線的比較。從圖12中可以看出,不同頻率作用下,壓力時(shí)歷曲線復(fù)雜程度均有所不同。對于彈性液艙,當(dāng)外部激勵(lì)頻率恰為共振頻率時(shí),液艙晃蕩更劇烈,產(chǎn)生的抨擊壓力也比另兩種工況要大。當(dāng)所加的外部激勵(lì)頻率遠(yuǎn)離共振頻率時(shí),不論是增大還是減小,對P1,P2點(diǎn)作用的抨擊壓力會減小,且前后兩個(gè)脈沖數(shù)值相差不明顯。當(dāng)外部激勵(lì)頻率為0.4f時(shí),壓強(qiáng)時(shí)歷曲線比2f作用下產(chǎn)生的時(shí)歷曲線更平緩。

圖12 P 1 和P 2 點(diǎn)的壓力時(shí)歷曲線對比Fig.12 Comparison for time history of pressure at the point P 1 and P 2

圖13是當(dāng)外部激勵(lì)頻率等于共振頻率時(shí),剛性液艙與彈性液艙上的3個(gè)壓力監(jiān)測點(diǎn)所監(jiān)測的壓力曲線對比圖。圖13(a)、(b)分別是P1點(diǎn)和P2的壓力歷時(shí)曲線對比圖,P1點(diǎn)在剛性液艙上的最大壓力約為2 180 Pa,而在彈性液艙上最大壓力約為1 980 Pa。最大壓力值減小了9.17%。P2點(diǎn)在剛性液艙上的最大壓力約為1 900 Pa,而在彈性液艙上最大壓力約為1 750 Pa,壓力減小了7.9%?？梢娕摫诘膹椥詫问帥_擊載荷起到了一定的緩和作用,且艙底彈性的緩和作用比艙兩壁的緩和作用更明顯。液艙晃蕩引起的沖擊載荷會產(chǎn)生雙峰值,彈性液艙中,第二次產(chǎn)生的峰值相較于第一次峰值減小的更明顯。

圖13(c)是P3壓力點(diǎn)的對比圖,在前9 s,彈性液艙中的沖擊載荷比剛性液艙中的值要小,但在9 s后,彈性液艙中P3壓力點(diǎn)比剛性液艙的峰值大,在11.3 s左右,晃蕩壓力增大約27.5%。15 s內(nèi),彈性液艙中最大壓力峰值約為1 380 Pa,剛性液艙中最大峰值約為1 230 Pa,增大約10.9%。但在彈性液艙中,晃蕩引起的空氣卷入現(xiàn)象要改善很多,從圖中可見,彈性艙中,P3壓力點(diǎn)上無負(fù)值出現(xiàn)。且在彈性液艙中,壓力歷時(shí)曲線表現(xiàn)的更光滑,歷時(shí)曲線形式比起剛性液艙更簡單。

圖13 當(dāng)外部激勵(lì)頻率等于共振頻率時(shí),剛彈性液艙上不同壓力點(diǎn)時(shí)歷曲線對比Fig.13 Comparison for different pressure points in elastic and rigid tanks when the external excitation frequency is equal to resonance frequency

綜合分析圖13的3幅對比圖可見,液艙的彈性對沖擊載荷有一定的緩和作用,且艙底的彈性緩和作用比艙兩壁緩和作用要更大。通過對比平均最大壓力峰值可見,彈性液艙中沖擊載荷比剛性液艙中沖擊載荷平均減小了約8.5%。此外,在彈性液艙中,峰值產(chǎn)生的相位相較于剛性液艙中有一定的延時(shí)作用。

圖14(a)是當(dāng)外部激勵(lì)頻率等于共振頻率作用下,自由液面達(dá)到左艙壁最高處時(shí)的現(xiàn)象圖。對比圖6(b)中的第一幅圖可見,在彈性液艙中,晃蕩引起的自由液面爬壁情況比剛性液艙中液面爬壁高度要低。圖6(b)表示在自由液面達(dá)到左艙壁最高處時(shí),作用在液艙上的范式等效應(yīng)力。從圖中可見,最大應(yīng)力為0.106 12 Pa,存在于左艙壁上。因此,在設(shè)計(jì)液艙結(jié)構(gòu)的過程中,需要對此部分進(jìn)行加強(qiáng),防止疲勞引起艙壁結(jié)構(gòu)破壞。

圖14 自由液面達(dá)到左艙壁最高處現(xiàn)象圖以及瞬時(shí)Stress Von Mises圖Fig.14 Impact process of the free surface arriving to the left wall highest place and the instantaneous Stress Von Mises

4 結(jié) 論

建立了剛性與彈性液艙的數(shù)值模型,模擬了不同外部激勵(lì)頻率作用下液艙晃蕩過程,主要結(jié)論如下:

1)采取的數(shù)值方法是可靠且準(zhǔn)確的,并且該方法可以很好的模擬液艙晃蕩過程中的碎波,液體粒子飛濺現(xiàn)象。

2)外部激勵(lì)頻率的變化對液艙晃蕩有一定的影響。在剛性液艙中,激勵(lì)頻率越大,產(chǎn)生的碎波,液體飛濺現(xiàn)象越嚴(yán)重,但因晃蕩引起的沖擊載荷卻相差不大。

3)在彈性液艙中,當(dāng)外部激勵(lì)頻率等于液艙晃蕩的共振頻率時(shí),作用在艙壁上的平均壓力最大。當(dāng)外部激勵(lì)頻率遠(yuǎn)離共振頻率時(shí),無論是增大還是減小,沖擊載荷會減小。

4)液艙的彈性作用對沖擊載荷有一定的減緩作用,且艙底的緩和作用比艙兩壁更明顯。在設(shè)計(jì)液艙結(jié)構(gòu)時(shí),因?yàn)榕搩杀诘膽?yīng)力比其他地方的應(yīng)力要大,所以需要注意艙兩壁的加厚或材料的選取。