船舶舷側與小型冰山碰撞數值模擬分析

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

1 大連理工大學 船舶工程學院，遼寧省深海浮動結構工程實驗室，遼寧 大連，116024 2 中汽研汽車檢驗中心（天津）有限公司，天津 300399 3 大連理工大學 工業裝備與結構分析國家重點實驗室，遼寧 大連，116024 4 高新船舶與深海開發裝備協同創新中心，上海 200240

0 引 言

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

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

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

1 流固耦合算法

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

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

質量守恒方程：

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

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

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

式中：K為流體單元的體積模量；Vc為流體單元的體積；A為位于耦合區域結構單元的平均面積，為了避免數值不穩定，引入了罰函數系數pf（0≤pf≤ 1）。

2 海冰的材料模型

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

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

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

圖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 舷側?冰山碰撞仿真模型

3.1 計算模型

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

圖4 舷側模型Fig. 4 Model of the side structure

表2 舷側板的厚度Table 2 Thickness of the side shell

表3 舷側板的主要材料參數Table 3 Main material parameters of the side shell

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

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

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

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

3.2 網格尺寸

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

表4 模型網格參數Table 4 Parameters of the model mesh

3.3 模擬工況

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

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

4 結果及分析

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

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

4.1 碰撞速度分析

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

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

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

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

4.2 應力分析

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

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

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

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

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

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

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

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

4.4 船體結構的能量吸收

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

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

5 結 語

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

在下一步的工作中，計劃引入船體的運動模型并開展船體防撞結構的優化設計。