船舶舷側(cè)與小型冰山碰撞數(shù)值模擬分析

楊碧野，黃志剛，劉寧，孫哲，張桂勇*,3,4

1 大連理工大學(xué) 船舶工程學(xué)院，遼寧省深海浮動結(jié)構(gòu)工程實驗室，遼寧 大連，116024 2 中汽研汽車檢驗中心（天津）有限公司，天津 300399 3 大連理工大學(xué) 工業(yè)裝備與結(jié)構(gòu)分析國家重點實驗室，遼寧 大連，116024 4 高新船舶與深海開發(fā)裝備協(xié)同創(chuàng)新中心，上海 200240

0 引 言

北極有著豐富的自然資源以及獨特的戰(zhàn)略位置，近年來，北極的開發(fā)逐漸成為國際熱點[1]。在北極，冰區(qū)航路環(huán)境十分復(fù)雜，海面上存在很多尺度不一的浮冰，甚至是體積巨大的冰山和冰脊。得益于現(xiàn)代造船技術(shù)的發(fā)展，目前冰區(qū)航行船舶的性能及結(jié)構(gòu)強度一般能夠滿足其在一定厚度的層冰和浮冰區(qū)內(nèi)自由航行，但冰山對于極區(qū)船舶來說仍然是一種威脅。雖然船舶與冰山碰撞的概率并不大，但是一旦發(fā)生，就會給船舶結(jié)構(gòu)以及船上人員安全帶來危險。據(jù)Hill[2]的統(tǒng)計，在1619～2004 年間，共發(fā)生了約670 起船舶與冰山的碰撞事故，這些事故不但造成了人員傷亡和經(jīng)濟損失，還給周圍環(huán)境帶來了嚴(yán)重破壞。海冰的密度與海水非常接近，規(guī)模龐大的冰山位于海面以上部分的體積很小，隱蔽性很高，即便是借助船用雷達，也很難及時觀測到這些海面上尺度較小的冰山[3]。特別是近年來隨著全球氣候變暖，小規(guī)模的冰山從融化的大型冰山中分離，導(dǎo)致航道內(nèi)小型冰山出現(xiàn)的概率增加，因此，船舶與小型冰山的碰撞問題以及考慮碰撞載荷下船舶結(jié)構(gòu)的設(shè)計問題亟待研究。

在船舶與冰山的碰撞問題中，涉及接觸、大變形、流固耦合等多種復(fù)雜的非線性問題，給理論計算和數(shù)值模擬帶來了挑戰(zhàn)。為方便預(yù)測碰撞力，有學(xué)者提出了相應(yīng)的理論和經(jīng)驗公式。Cammaert等[4]在考慮海冰擠壓破壞時的能量損失的基礎(chǔ)上，提出了浮冰或者冰山與海洋結(jié)構(gòu)物碰撞的簡化公式；Liu 等[5]基于Stronge 的碰撞力學(xué)模型，從外部動力學(xué)的角度提出了用于分析船?船或者船?冰碰撞的三維數(shù)學(xué)模型；Timco[6]通過分析大量的全尺度實測數(shù)據(jù)以及模型試驗數(shù)據(jù)，認為海冰與海洋結(jié)構(gòu)物之間碰撞力的大小與海冰的動能直接相關(guān)。這些理論或經(jīng)驗公式在一定程度上可以快速預(yù)測冰山與船舶的碰撞力，方便研究碰撞過程中的機理，但假設(shè)眾多，過于簡化，無法捕捉到碰撞中一些細節(jié)現(xiàn)象，并且不能滿足現(xiàn)代船舶的需要。為了能夠較為精細地模擬船舶與冰山的碰撞過程，以非線性有限元法為主的數(shù)值模擬方法在近些年逐漸得到應(yīng)用。Gagnon[7]采用泡沫材料模型模擬了冰山與船舶或者海洋結(jié)構(gòu)物的碰撞過程，重點模擬了冰山的壓碎行為。Liu 等[8]按照Tsai-Wu屈服面準(zhǔn)則，提出了一種適用于冰山的材料模型，并將此模型運用到了球鼻艏與冰山相撞的模擬中。Gao 等[9]分別采用彈塑性材料模型和泡沫模型對不同局部形狀冰山與船體的碰撞過程進行了研究。在數(shù)值模擬中，大部分學(xué)者主要關(guān)注冰山的材料模型以及船體的結(jié)構(gòu)變形，而對于碰撞過程中水的影響則考慮得并不全面，采用常數(shù)形式的附加質(zhì)量系數(shù)并不能完全反映水動力的影響，而且冰山與結(jié)構(gòu)物之間會出現(xiàn)例如波浪緩沖效應(yīng)等自由液面變化。為了得到精確的結(jié)果，需要同時進行流體的計算。

本文將主要基于商業(yè)軟件LS-DYNA 對小型冰山與船舶舷側(cè)結(jié)構(gòu)的碰撞進行數(shù)值模擬。首先，采用多物質(zhì)任意拉格朗日?歐拉（ALE）算法進行流固耦合計算；然后，在不同碰撞角度下對船舶舷側(cè)與冰山的碰撞過程進行仿真，同時對碰撞過程中的現(xiàn)象和不同工況下的碰撞力、碰撞速度以及能量吸收情況進行分析，用以指導(dǎo)極地船舶的建造和航行操作。

1 流固耦合算法

本文采用ALE 算法模擬冰山與舷側(cè)結(jié)構(gòu)在水中的碰撞問題。在ALE 算法中，采用歐拉方法對流體的運動和變形進行描述，而對于船體結(jié)構(gòu)和冰山，則采用拉格朗日描述，并在每個時間步內(nèi)進行載荷和邊界的信息交換。

ALE 算法控制方程包括以下基本方程。

質(zhì)量守恒方程：

式中：ρ 為材料密度；v為物質(zhì)質(zhì)點的速度；w為網(wǎng)格的速度；σ 為應(yīng)力張量；g為重力加速度；e為能量；t為時間。

在流固耦合交界面處，流體與固體之間的信息傳遞采用罰函數(shù)法，即在每個時間步內(nèi)，一旦流體節(jié)點穿過結(jié)構(gòu)表面，就會產(chǎn)生一個耦合力Fc將貫穿的流體節(jié)點推回至耦合交界面。該耦合力的大小正比于貫穿深度d，圖1 所示為該耦合算法示意圖（圖中，tn為第n個時間步），其中數(shù)值剛度k為：

圖1 耦合算法示意圖Fig. 1 Schematic diagram of coupling algorithm

式中：K為流體單元的體積模量；Vc為流體單元的體積；A為位于耦合區(qū)域結(jié)構(gòu)單元的平均面積，為了避免數(shù)值不穩(wěn)定，引入了罰函數(shù)系數(shù)pf（0≤pf≤ 1）。

2 海冰的材料模型

冰的本構(gòu)關(guān)系十分復(fù)雜，其力學(xué)性質(zhì)與應(yīng)變率、溫度等外在和內(nèi)在物理因素有關(guān)[10-12]。在船舶與冰山的碰撞研究方面，國內(nèi)外學(xué)者在冰山材料的選擇上主要有3 種：雙線性彈塑性模型(MAT_PLASTIC_KINEMATIC)[13]、根 據(jù)Tsai-Wu 屈 服 準(zhǔn)則定義的彈塑性模型[8,14]，以及由Gagnon[7]提出的可壓碎泡沫材料(MAT_CRUSHABLE_FOAM)。由于海冰復(fù)雜的力學(xué)性質(zhì)，目前還沒有一種材料模型能夠完全描述海冰的力學(xué)行為，即便是上述3 種比較常見的冰山材料模型，也存在一些缺點，例如可壓碎泡沫材料雖然可以模擬船舶與冰山碰撞過程中的壓碎行為，但在理論上缺乏對海冰材料性質(zhì)的系統(tǒng)描述，無法模擬海冰裂紋等問題。本文研究關(guān)注的重點是碰撞力和船體舷側(cè)的應(yīng)力變化情況。因冰山的體積較小，破壞并不是十分明顯，因此將采用雙線性彈塑性材料模型來模擬冰山。其中，冰山的材料參數(shù)參考文獻[15]的北極冰山實驗數(shù)據(jù)，如表1 所示。

表1 冰山的主要材料參數(shù)Table 1 Main material parameters of iceberg

本文將以一個球形冰與剛性板的碰撞試驗為例來驗證所用材料模型的可行性，圖2 所示為計算工況示意圖（圖中虛線表示對這條線以下進行約束）。球形冰半徑為2 m，球形冰遠離碰撞位置的一側(cè)剛性約束，剛性板以1 m/s 的速度與球形冰發(fā)生碰撞。海冰的壓力?面積曲線常被用來記錄海冰的力學(xué)性質(zhì)。Kim 等[17]基于大量的試驗室內(nèi)海冰以及自然海冰的相關(guān)力學(xué)試驗，總結(jié)得到壓力?面積公式P=0.35A?0.5，如圖3 中紅線所示。同時，為檢驗網(wǎng)格大小對計算結(jié)果的影響，分別采 用0.1 m×0.1 m×0.1 m，0.15 m×0.15 m×0.15 m和0.2 m×0.2 m×0.2 m 這3 種不同的網(wǎng)格尺寸進行了數(shù)值模擬，并將記錄得到的壓力?面積曲線與Kim 等[17]試驗所得的曲線進行了對比，如圖3 所示。由圖可見，3 種網(wǎng)格尺寸下的壓力?面積關(guān)系與Kim 等[17]試驗所得曲線的規(guī)律相似，并且數(shù)值模擬得到的關(guān)系曲線也基本上在經(jīng)驗曲線周圍波動，說明了所采用材料模型的可行性。

圖2 球形冰與剛性板的碰撞示意圖Fig. 2 Schematic diagram of the collision between spherical ice and rigid plate

圖3 壓力?面積曲線Fig. 3 Curves of pressure with nominal contact area

3 舷側(cè)?冰山碰撞仿真模型

3.1 計算模型

在上述流固耦合模型以及冰材料模型的基礎(chǔ)上，本文對船舶舷側(cè)與小型冰山的碰撞過程進行了數(shù)值模擬。研究對象為一艘抗撞擊能力較強的油船，選擇船舶平行中體的部分雙層舷側(cè)結(jié)構(gòu)作為計算模型。圖4 所示為數(shù)值模擬中的舷側(cè)模型，其結(jié)構(gòu)包含了雙層舷側(cè)板及其之間的肋板、縱桁等構(gòu)件，雙層舷側(cè)板之間的距離為1 m。整個舷側(cè)構(gòu)件使用LS-DYNA 軟件中的Shell 單元進行建模，板厚按照表2 中參數(shù)。船體板的材料選擇Q235 鋼，材料模型同樣選擇MAT_PLASTIC_KINEMATIC 模型，主要材料參數(shù)如表3 所示。

圖4 舷側(cè)模型Fig. 4 Model of the side structure

表2 舷側(cè)板的厚度Table 2 Thickness of the side shell

表3 舷側(cè)板的主要材料參數(shù)Table 3 Main material parameters of the side shell

基于Gao 等[9]的研究成果，發(fā)現(xiàn)船舶與冰山碰撞時，尖銳形狀的冰山在碰撞過程中容易發(fā)生破碎現(xiàn)象，同等條件下鈍形冰的接觸面將有更大的碰撞力，會對船舶造成更大的傷害，因此數(shù)值模擬中選擇的是具有代表性的球形冰山。模型中，球形冰山的直徑為4 m，材料模型參數(shù)如表1 所示，冰山的密度為900 kg/m3，水的密度為1 025 kg/m3，對應(yīng)的冰山露出高度約為0.81 m。在建立水域和空氣域時，保證兩者的交界面與船舶吃水線位置相同。船舶吃水為9 m，流體計算域尺度大于雙側(cè)舷側(cè)和球形冰山的尺度，且流體域要在空間內(nèi)包括整個模型，整個計算域如圖5 所示。

圖5 整體計算域模型Fig. 5 The model of whole computational domain

對于流體部分，采用NULL 材料模型，搭配線性多項式狀態(tài)方程描述流體變形與壓力之間的關(guān)系，其中水的黏性系數(shù)設(shè)置為8.684×10?4(N·S)/m2，整個計算模型受到的重力加速度為9.8 m/s2。

采用侵蝕算法，即CONTACT_ERODING_SURFACE_TO_SURFACE 關(guān)鍵字處理船體結(jié)構(gòu)與冰山之間的碰撞，設(shè)置CONSTRAINED_LAGRANGE_IN_SOLID 關(guān)鍵字用于處理結(jié)構(gòu)和冰山與水之間的流固耦合。同時，為了避免冰山的結(jié)構(gòu)網(wǎng)格與水域的網(wǎng)格之間重疊所帶來的求解錯誤，INITIAL_VOLUME_FRACTION_GEOMETRY功能被用來進行流固耦合界面的初始化。

3.2 網(wǎng)格尺寸

根據(jù)LS-DYNA 軟件理論指導(dǎo)手冊的建議[19]，當(dāng)流體域的網(wǎng)格尺寸與結(jié)構(gòu)物的網(wǎng)格尺寸相近時，其流固耦合ALE 算法將有更理想的計算結(jié)果。因此，在對舷側(cè)、冰山以及水域劃分網(wǎng)格的過程中，同時考慮了計算精度以及計算效率的問題，船體結(jié)構(gòu)采用Shell 單元建模，將舷側(cè)肋板的網(wǎng)格尺寸定義為0.25 m，舷側(cè)模型其他區(qū)域的網(wǎng)格尺寸定義為0.3 m。球形冰山采用Solid 單元建模，根據(jù)圖3 所示的網(wǎng)格無關(guān)性結(jié)果，球形冰山的網(wǎng)格尺寸選擇為0.15 m，水域和空氣域的網(wǎng)格單元尺寸定義為0.3 m。計算模型各部分的網(wǎng)格類型以及對應(yīng)的網(wǎng)格數(shù)量如表4 所示。

表4 模型網(wǎng)格參數(shù)Table 4 Parameters of the model mesh

3.3 模擬工況

在數(shù)值模擬中對舷側(cè)邊緣兩側(cè)進行剛性約束，使冰山以一定的速度和角度與舷側(cè)碰撞。為了增加冰山運動過程中的穩(wěn)定性，將浮冰z方向的自由度進行位移約束，其余方向自由運動。碰撞前，冰山距離舷側(cè)法向的初始距離為1.5 m，以初速度V=5 m/s 向船舶運動。冰山與舷側(cè)碰撞角度的俯視圖如圖6 所示。本文分別考慮了冰山與舷側(cè)的碰撞角度α = 30°，45°，60°，75°和90°時的工況，模擬的物理時間為2 s。需要說明的是，本文所模擬的冰山撞擊船舶與船舶撞擊冰山在物理過程上略有不同，且未考慮船舶運動的影響，作為一種簡化的處理手段[9,20-22]，模擬冰山撞擊船舶可以更好地捕捉該過程中的流體阻滯效應(yīng)和結(jié)構(gòu)變形。

圖6 碰撞角度俯視圖Fig. 6 Top view of collision

4 結(jié)果及分析

選取舷側(cè)與冰山的碰撞角度α = 60°的情況進行分析。圖7 所示分別為數(shù)值模擬過程中在初始狀態(tài)、0.45 s 以及1.21 s 不同時刻（T）舷側(cè)與冰山的相對位置以及碰撞過程中的自由液面變化情況。由圖可以看出，隨著冰山的運動，自由液面變化劇烈，在冰山與舷側(cè)結(jié)構(gòu)之間逐漸興起波浪，甚至導(dǎo)致舷側(cè)上發(fā)生了波浪砰擊，這與Song等[22]在模型試驗中所觀察到的現(xiàn)象類似。

圖7 α = 60°時不同時刻舷側(cè)?冰山碰撞現(xiàn)象Fig. 7 Collision phenomenon between side structure and iceberg with α = 60° at different times

4.1 碰撞速度分析

圖8 所示為碰撞過程中冰山在3 個方向的速度改變曲線圖。由圖可以看出，垂直于舷側(cè)方向的速度Vy由于受水的阻力以及波浪緩沖作用的影響，在碰撞發(fā)生前速度逐漸減?。辉诩s0.51 s時，由于舷側(cè)與冰山發(fā)生碰撞，Vy急劇減小，直至發(fā)生回彈；而平行于舷側(cè)方向的速度Vx變化則比較緩慢，僅在碰撞時刻發(fā)生了小幅度的跳躍變化。與Vy相比，Vx在碰撞的瞬間更多地是受到摩擦力的影響，而Vy則是受到正面碰撞的影響，因此Vx的速度跳躍變化相對較小。此外，在碰撞前的速度衰減方面，Vy也更加明顯。從圖8 所示的自由液面變化中可以發(fā)現(xiàn)，舷側(cè)法向的興波更加劇烈，波浪緩沖影響更大，而Vx更多地是受到流體阻力的影響，變化不明顯。

圖8 α = 60°時的冰山速度變化情況Fig. 8 Velocity of iceberg with α = 60°

圖9 所示為不同碰撞角度下，初始速度相同的冰山在碰撞過程中速度幅值隨時間變化的情況。由圖可以看到，隨著碰撞角度的增大，冰山在接觸法相方向上的速度分量隨之增大，速度幅值的衰減越發(fā)明顯。如上所述，在此過程（圖7）中，碰撞角度越大，冰山和船體結(jié)構(gòu)之間的法向興波越劇烈，冰山因水動力的影響所受到的斥力作用也越明顯，在斥力和阻力的共同作用下，速度也就會發(fā)生更快速的衰減。

圖9 不同碰撞角度下的冰山速度幅值變化Fig. 9 Velocity amplitude of iceberg with different collision angles

4.2 應(yīng)力分析

圖10 所示為碰撞角度α = 90°時舷側(cè)在不同時刻的等效應(yīng)力圖。在碰撞過程中，應(yīng)力云圖的面積逐漸增大，其中碰撞區(qū)域中心的應(yīng)力最大，然后向周圍逐漸變??；在碰撞接觸最充分的時刻，應(yīng)力云圖的區(qū)域面積以及中心應(yīng)力達到最大值，碰撞結(jié)束后，高應(yīng)力區(qū)域面積隨之減小。對應(yīng)過程中的冰山應(yīng)力云圖如圖11 所示，其與碰撞過程中舷側(cè)的應(yīng)力特征相似，即冰山碰撞區(qū)域中心的應(yīng)力最大，周圍的應(yīng)力較小。同時，從圖11(b)中可以看到，在垂直碰撞過程中，冰山發(fā)生了輕微的破碎現(xiàn)象。

圖10 α = 90°時舷側(cè)不同時刻的等效應(yīng)力圖Fig. 10 Equivalent stress contours of the side structure with α = 90° at different times

圖11 垂直碰撞過程中不同時刻的冰山應(yīng)力圖Fig. 11 Stress contours of the iceberg at different times during vertical collision

4.3 不同碰撞角度下的碰撞力峰值

圖12 所示為舷側(cè)與冰山在不同碰撞角度下的碰撞力曲線。從圖中可以看出，隨著碰撞角度的增大，碰撞時刻不斷提前。當(dāng)α = 30°，45°，60°，75°時，碰撞力曲線的變化趨勢相似，即在某一時刻碰撞力達到最大值后開始變小，碰撞力曲線光滑；且在這4 種碰撞角度下，均發(fā)生了冰山與舷側(cè)碰撞的反彈過程，冰山未發(fā)生破壞。當(dāng)α = 90°時，碰撞力發(fā)生了波動，其變化趨勢與其他碰撞角度下的碰撞力曲線不同，原因是冰山在這種情況下發(fā)生了破碎，導(dǎo)致碰撞力出現(xiàn)波動現(xiàn)象。

圖12 不同工況下的碰撞力曲線Fig. 12 Curves of collision force in different cases

在圖12 所示各種碰撞角度下的碰撞力曲線中，提取碰撞力峰值并繪制碰撞力峰值與碰撞角度的關(guān)系曲線圖如圖13 所示。由圖中的曲線趨勢可以看出，隨著碰撞角度的增大，碰撞力峰值呈明顯上升的趨勢，在α = 30°～75°的階段，碰撞力峰值上升明顯，且增長幅度基本相同；在α = 75°～90°范圍內(nèi)，碰撞力稍有增長，增長幅度較小。其原因可能是當(dāng)α = 30°，45°，60°，75°時，在碰撞的過程中冰山未發(fā)生破碎，而在α = 90°情況下冰山發(fā)生了局部破碎；另外，在α = 90°工況下，除碰撞過程中舷側(cè)的塑性變形消耗能量外，一部分能量還被冰山破碎所消耗。

圖13 不同碰撞角度下的碰撞力峰值Fig. 13 Peak collision force with different collision angles

4.4 船體結(jié)構(gòu)的能量吸收

圖14 給出了不同碰撞角度下船體結(jié)構(gòu)的能量吸收情況。碰撞過程中，很大一部分能量被船殼的塑性拉伸所吸收，但在當(dāng)前算例中，由于冰山的體積較小，在與船體結(jié)構(gòu)發(fā)生碰撞時，船體結(jié)構(gòu)并未發(fā)生明顯破壞，碰撞結(jié)束后，船體板在發(fā)生一定的塑性變形后回彈，能量吸收曲線表現(xiàn)為迅速增加后略微下降。由圖14（a）～圖14（d）可以看到，舷側(cè)外板為船體結(jié)構(gòu)的主要能量吸收部件，且隨著碰撞角度的增加，分配到舷側(cè)外板上的能量越多。而當(dāng)α = 90°時（圖14（e）），由于在碰撞過程中冰山破碎帶走了一部分能量，分配到船體結(jié)構(gòu)上的碰撞能量明顯下降，且除舷側(cè)外板之外的其他構(gòu)件也在能量吸收中占據(jù)有相當(dāng)?shù)谋壤?/p>

圖14 不同碰撞角度下船體結(jié)構(gòu)的能量吸收Fig. 14 Energy absorption of ship structure at different collision angles

5 結(jié) 語

本文基于非線性有限元方法和ALE 算法建立了船舶舷側(cè)結(jié)構(gòu)與小型冰山相互作用的數(shù)值模型，并在模型中考慮了船舶結(jié)構(gòu)變形、海冰破壞以及兩者之間水動力的影響。其中，由海冰材料模型所得的壓力?面積曲線與試驗所得結(jié)果吻合較好，且整個數(shù)值模擬結(jié)果在現(xiàn)象上也與模型試驗結(jié)果相似。然后，在上述基礎(chǔ)上，研究了小型冰山以不同角度與舷側(cè)結(jié)構(gòu)發(fā)生碰撞時的過程中速度、碰撞力和船體結(jié)構(gòu)能量吸收的變化。經(jīng)分析發(fā)現(xiàn)，在冰山向舷側(cè)運動的過程中，由于興波消耗能量以及水的阻力作用，冰山的速度逐漸減小，冰山在以初速度運動撞擊舷側(cè)后發(fā)生了回彈；船舶舷側(cè)結(jié)構(gòu)的應(yīng)力區(qū)域面積在碰撞過程中逐漸增加，應(yīng)力由中心向四周逐漸減小，在碰撞接觸最充分時產(chǎn)生了最大的中心應(yīng)力和應(yīng)力面積；隨著舷側(cè)與冰山碰撞角度的逐漸增大，碰撞力和船體吸收能量相應(yīng)增加，但增幅不同，并且冰山是否發(fā)生破碎會影響到碰撞力和船體能量吸收的大小。

在下一步的工作中，計劃引入船體的運動模型并開展船體防撞結(jié)構(gòu)的優(yōu)化設(shè)計。