船-冰碰撞下冰材料有限元數值方法研究進展

蔡 偉，朱 凌，3，畢璐澤

（1.高性能船舶技術教育部重點實驗室（武漢理工大學），武漢 430063；2.武漢理工大學交通學院，武漢 430063；3.高新船舶與深海開發裝備協同創新中心，武漢 430063）

0 引 言

極地資源的開發、航道運輸以及科考成為當今世界各國關注的熱點話題。隨著全球氣候變暖，北極冰川加速融化，極地生產運輸活動逐漸增多，致使浮冰，特別是冰山逐漸增多，這給航行于冰區的船舶的結構安全帶來了巨大的挑戰。極地地區海面上漂浮的冰體主要分為淡水冰和咸水冰，其中淡水冰主要是冰山冰，或者是一些從冰川分離而出的小型冰山。極地船舶在冰區航行時難免會遭受到浮冰的碰撞作用，在極端冰碰載荷作用下船體結構難免會出現塑性變形，甚至會出現破壞斷裂，如圖1和圖2所示，這將會導致人員傷亡、貨物泄露以及環境污染等問題，甚至會造成沉船事故。例如，2015年，俄羅斯大型拖網漁船“遠東”號與浮冰相撞，導致船舶沉沒，56人死亡以及13人失蹤；2019年1月，“雪龍”號破冰船在南極阿蒙森海與冰山相撞，船艏桅桿被撞倒，部分舷墻受損，對考察任務造成了嚴重影響。因此，亟需開展冰碰載荷下船舶結構動力響應評估方法研究，以探討冰體碰撞下結構的塑性變形損傷。

船舶與海洋工程結構與浮冰碰撞是結構在短時間內遭受巨大沖擊載荷作用下的一種復雜非線性動態響應過程，存在著材料非線性、幾何非線性、接觸非線性和運動非線性等問題。針對船-冰碰撞結構動力響應問題，目前尚未有較為成熟全面的理論方法，尤其是在冰材料力學模型方面，目前國內外學者主要是采用經驗公式法、試驗方法和數值方法來研究冰體碰撞下結構動力響應問題。隨著計算機技術的不斷提高，計算力學方法在不斷發展，數值模擬方法已成為一種最為經濟高效的方法[1]，主要包括有限元方法（FEM）、離散元方法（DEM）[2]、光滑粒子動力學方法（SPH）[3]以及近場動力學方法（PD）[4]等。有限元法在模擬介質間的間斷和離散性方面能力有限，相比之下，離散元方法在模擬冰體離散特性方面具有一定的優勢[5]。除此之外，有限元方法在模擬冰體斷裂以及裂紋擴展方面有一定的局限性，而光滑粒子動力學方法和近場動力學方法的粒子離散形式使其在海冰破碎以及裂紋擴展方面具有更突出的優勢[4]。然而，在船舶結構動力響應研究方面，尤其是針對結構塑性大變形以及斷裂失效問題，有限元方法比離散元、光滑粒子動力學以及近場動力學方法顯得更有優勢，更加成熟。

由于冰體是一種力學性質較為復雜的天然復合材料，其材料力學特性的影響因素眾多，尤其是對于海冰材料，主要受鹽分、雜質等成分含量、溫度、加載速率和形成年份等影響[6]。船-冰碰撞過程是一種瞬態動力響應問題，由于冰材料的固有特性，冰隨著應變率的增加會由韌性破壞到脆性破壞[7-8]。由于冰體材料的復雜力學特性，不同條件下冰的力學性質差異較大，其材料數值模型仍然處于研究探索階段。目前，國內外學者提出冰體理想化材料模型，主要為各向同性彈性失效模型[9-14]、各向同性彈塑性失效應變率相關材料模型[15-18]、可壓碎泡沫型材料模型[19-25]、損傷材料模型[26-27]以及一些其它自定義彈塑性材料數值模型[28-35]等。在船-冰碰撞過程中，冰體的破壞包括了擠壓、彎曲破壞、堆積等動態過程，其破壞力學特性不僅與應變率、溫度等因素有關，還與其自身形成條件有密切關系，如年份、鹽分、雜質等成分含量，利用數值模型方法來研究冰體破壞失效過程具有較大的挑戰性。除此之外，通過一些冰體壓縮、彎曲等試驗測得冰體材料參數，如楊氏模量、屈服應力、剪切模量和泊松比等，這些結果具有較大的離散性，開發一種涵蓋所有條件下的冰體材料模型是不現實的。因此，目前尚未有良好的數值模型來描述所有條件下的冰體材料力學特性，上述所提出的材料模型都有各自的局限性。

因此，本文主要總結了現有幾種常見的冰體有限元數值材料模型，并與現有的一些試驗結果、ISO 推薦的壓力-面積曲線以及單位體積能量吸收經驗值進行對比分析，討論現有一些冰體材料模型的適用性和優缺點。

1 冰材料有限元數值模型研究進展

有限元數值計算方法是船-冰碰撞下結構動力響應的重要評估手段，但是碰撞過程中冰材料的數值模擬是其最為關鍵的難點問題，下面將介紹幾種常見的冰材料數值模型。

在大多數船-冰碰撞情況過程中，船體與冰在接觸過程中高應變率使冰體較明顯地體現出脆性性質，故國內外許多學者將冰體材料模型簡化為各向同性彈性失效材料模型。最早在針對冰雹與復合材料結構碰撞響應進行研究時，Kim 等[9-10]以及Anghileri 等[11]利用ANSYS/LS-DYNA 軟件中的各向同性彈性斷裂模型來模擬在碰撞過程中冰體的力學特性。該模型采用Von-mises 失效準則，將最大塑性應變定義為其破壞模式，把恒定最小壓力模式定義為其分離模式，來模擬冰體失效行為。雖然冰體材料模型比較簡單，可以節省計算時間，但是僅僅依據最大塑性應變和截斷壓力兩個恒定值來判斷冰體失效模式有所不足。結合相關冰體力學參數，一些學者[12-14]逐漸將此模型運用于船-冰碰撞分析之中。考慮到冰體材料力學特性受應變率影響較大，Carney 等[15]采用一種各向同性彈塑性失效應變率相關材料模型來模擬冰體失效，此模型廣泛運用于冰雹的高速沖擊數值模擬[16-17]。胡志寬[18]將各向同性彈塑性失效應變率相關材料模型運用于冰與螺旋槳的碰撞數值模擬研究中。此種材料模型雖然考慮到了應變率變化對冰體力學特性的影響，但是在船-冰碰撞結構動態響應研究中，碰撞速度有限，其應變率變化范圍并不明顯，沒有冰雹高速沖擊的應變率大，故該材料模型的適用性具有一定的局限性。

由于冰體材料參數試驗測量數據離散性較大，具體的數值大小難以確定，一種可壓縮泡沫型材料模型[19-21]被廣泛運用于船-冰碰撞研究之中，通過運用自定義應力-體積應變關系來定義可壓碎泡沫模型的力學行為。考慮到在碰撞過程中冰體與結構物接觸區域存在著高壓區，在其接觸區域周圍存在著低壓區，如圖3所示，但是上述的可壓縮泡沫型材料模型還不能體現出這一區別。為此，Gagnon[22]對此材料模型所定義應力-體積應變關系進行了改進，通過對冰體模型進行分層，在每層賦予不同材料特性的鋸齒形的應力-體積應變關系曲線，如圖4所示，來對冰體的高壓區域和低壓區域進行模擬。類似于Gagnon 的分層方法，Kim 等[23-25]通過對冰體壓縮試驗現象與結果的觀察和分析，定義了兩條鋸齒形的應力-體積應變關系曲線，并賦予冰體不同區域處，以表達冰體與結構物接觸時高壓區域和低壓區域。此外，Kim 還運用“最大主應力”失效準則用來定義冰體破壞，即當冰體網格單元大于最大失效應力時，單元自動失效刪除，采用此種材料模型的數值計算結果與試驗結果吻合較好。相比于其它冰材料模型，可壓碎泡沫材料模型在描述冰體力學特性時在一定程度上有較高的準確性，它可以宏觀地體現出冰體與結構物相互接觸時冰體的變化，可以比較真實地反映實驗過程中的冰體與結構物的接觸力大小，但是其缺乏強有力的物理解釋，無法模擬出冰體脆性失效特性以及冰體在碰撞過程中出現的內部缺陷，該模型只反映特定條件下的冰體材料力學行為。

冰體在碰撞過程中主要體現出脆性失效性質，考慮到巖土、混凝土等材料本構模型與冰體材料在某種程度上有一定的相似性，而混凝土材料模型研究相對于冰體材料數值模型研究較為成熟，一些學者逐漸將一些成熟的混凝土損傷材料模型用于冰體材料模型數值模擬當中。Polach 和Ehlers[26]使用ANSYS/LS-DYNA 軟件中損傷材料模型（*MAT_DAMAGE_3）來模擬冰體材料，使用Lemaitre 損傷模型來模擬冰體失效，當冰體損傷值大于臨界損傷值時，冰材料就會失效。該模型的數值計算結果與實驗結果吻合較好[37]。此外，Cai 等[27]基于混凝土材料本構模型，結合冰體材料力學模型，建立了表述應力-應變關系的本構方程的屈服面公式，以及失效應變與截斷壓力之間關系等式，提出了基于混泥土材料本構模型的冰體數值材料模型，數值計算結果與模型試驗結果吻合較好。由于冰體材料特性的影響因素眾多，如溫度、應變率和鹽度，開發出一種涵蓋所有條件下的冰體材料模型是不現實的。相對于海冰來說，冰山冰為淡水冰，材料屬性較為簡單，可將其視作各向同性材料[38]。Liu 等[28-29]通過開發LS-DYNA 軟件中的自定義各向同性的彈塑性材料模型來模擬冰山的力學行為，在所開發的彈塑性模型中，使用屈服面公式和流動法則來表征冰材料的應力-應變關系，建立失效應變與壓力之間的關系曲線來定義冰體的失效準則。基于Liu 等人的自定義冰體彈塑性材料模型，Gao 等[31]提出了一種新的失效準則來模擬冰體失效，Song 等[33]利用Liu 的冰材料模型開展了考慮附連水質量的船-冰碰撞有限元計算研究。此外，Shi 等[34]提出了一種與溫度梯度相關的冰體彈塑性材料模型，用于研究溫度對船舶-冰山碰撞沖擊力的影響。上述的自定義材料模型雖然在冰體失效準則和冰體三向應力狀態模擬上更加準確，但是依然很難反映冰體真實裂紋延伸和斷裂形態，相比于可壓縮材料模型，無法準確模擬計算出冰體在碰撞過程中接觸面壓力大小分布特征。除此之外，Ince等[35]基于KOSORI冰模型的本構關系[36]，考慮應變速率、溫度和鹽度等參數的影響，對Johnson-Cook 材料模型進行了修正，通過ABAQUS軟件用戶自定義材料模型（UMAT）來實現冰材料模型的數值模擬。

在自然界中，無論是海冰還是淡水冰，不同環境下冰體的形成機理各不相同[6,36,38-39]，如在極地地區有1年冰和多年冰，不同年份的冰層的材料力學特性相差較大。冰的應力-應變特性依賴于應變率和溫度，與冰的形成成分如鹽分等也有直接關系，不同條件下測得的冰體材料參數也各不相同，如冰體密度、屈服應力、楊氏模量、剪切強度、拉伸強度、壓縮強度以及斷裂韌度等參數。開發一種涵蓋所有條件下的冰體材料模型是難以實現的，每種冰體數值材料模型都具有其局限性。

2 模型試驗與數值模擬結果對比分析

船-冰相互作用試驗研究方法主要包括壓縮和拉伸試驗、彎曲試驗、貫穿與剪切實驗、壓痕試驗以及碰撞試驗[5]。其中，壓縮與拉伸試驗、彎曲試驗、貫穿與剪切實驗以及壓痕試驗主要用于冰體力學特性研究。在船-冰相互作用過程中，有些工況是低速大質量碰撞情形，通常采用準靜態加載實驗方法來研究結構的彈塑性動力響應；對于碰撞速度較大時，便采用動態沖擊實驗方法來研究結構的彈塑性動力響應。下面針對三種模型試驗工況來進行數值模擬分析，并對不同冰材料模型的計算結果進行分析。值得說明的是，由于所分析的模型實驗工況有各自的特點，所以對下面所選取的冰材料模型參數與所參考的原文中的冰材料模型進行了微小的改動，以提高數值計算結果與模型實驗結果的可比性。

2.1 冰體壓縮實驗數值模擬

Kim 等[24]開展了不同加載速率圓錐冰體壓縮實驗，實驗裝置如圖5 所示，還采用可壓縮泡沫型冰體材料模型來進行相應的數值仿真研究。在數值模擬中通過定義兩條鋸齒形的應力-體積應變關系曲線，并賦予冰體不同區域處，采用“最大主應力”失效準則來定義冰體破壞失效準則，數值結果與實驗結果吻合較好。將Kim 等人所開展的冰體壓縮實驗其中一組實驗工況來作為有限元數值模擬對象，其主要實驗參數如表1所示，有限元數值模型如圖6所示。圖7對比了模型實驗和6組數值模型計算結果，可以發現：各向同性彈性失效模型（MAT_13）和混凝土損傷模型（MAT_78）的計算結果會出現多個峰值，并出現多次加卸載過程。這是由于這兩種材料模型在加載過程中冰模型會出現網格失效和刪除，導致壓塊與冰體存在多次接觸的過程。因此，在模擬冰體壓縮實驗過程中，采用網格單元失效刪除的方式來模擬冰體失效是不準確的，很難對冰體壓縮力-位移變化機制進行預測，會導致最終所預測的載荷值偏小。

表1 冰體壓縮實驗和數值模擬主要參數Tab.1 Main parameters of ice compression test and numerical simulation

圖7 圓錐冰壓縮實驗與數值的力-位移曲線結果對比Fig.7 Comparison of experimental and numerical force-displacement curves for conical ice compression

2.2 冰體水平碰撞實驗數值模擬

Cai 等[27]利用水平沖擊試驗機開展了楔形冰與船體板模型碰撞實驗，實驗裝置如圖8 所示。本文選取的碰撞實驗工況主要參數如表2 所示，對應的有限元數值模型如圖9 所示，冰體前端長度為200 mm，前端夾角為90°。下面選取了5種材料模型來進行對比，模型實驗與有限元數值計算結果對比如圖10所示。圖10對比了各向同性彈性失效模型（MAT_13）、可壓碎泡沫型材料模型（MAT_63）以及損傷模型（|MAT_78，MAT_153），從對比結果可以發現：采用低壓區的低應力-體積應變曲線（M2）的可壓碎泡沫型材料模型的計算結果與實驗結果相比誤差比較大；在碰撞力預測方面，采用各向同性彈性失效模型和損傷材料模型（MAT_78）的計算結果與實驗結果比較接近；在船體板結構變形預測方面，采用損傷材料中的MAT_78 和MAT_153 單元材料模型的計算結果與實驗結果比較接近。在船-冰碰撞數值模擬過程中，由于可壓碎泡沫型材料模型很難模擬出冰體破壞失效特性，在模擬冰體破壞幾何特征方面存在很大的不足之處。相比較而言，船-冰碰撞過程中冰體與船體局部結構相互作用特性研究方面，由于冰體在碰撞擠壓過程中會出現以顆粒形式剝落的現象，采用如各向同性彈性失效模型（MAT_13）和損傷材料模型（MAT_78、MAT_153）等具有單元失效刪除功能的冰材料模型來進行數值模擬具有較大的優勢。

表2 冰體水平碰撞實驗主要參數Tab.2 Main parameters of ice horizontal impact test

圖10 楔形冰碰撞實驗和數值的力-位移曲線結果對比[27]Fig.10 Comparison of experimental and numerical force-displacement curves for wedge ice collision

2.3 冰體落體碰撞實驗數值模擬

Ince等[35]開展了圓錐冰體落體碰撞實驗，實驗裝置如圖11所示，此外還進行相對應的有限元數值模擬，基于KOSORI冰模型本構關系，采用了改進的Johnson-cook材料模型來模擬冰體力學行為，數值計算結果與實驗結果吻合較好。本文選取的碰撞實驗工況主要參數如表3 所示，對應的有限元數值模型如圖12 所示。從圖13 的對比結果可以發現：不同的冰材料模型數值計算結果差距較大，在板的變形為60 mm之前時，MAT_63（M2）材料單元模型、MAT_78單元模型與KOSORI冰模型的計算結果比較接近，與實驗結果吻合較好：當板的變形超過60 mm時，上述的冰材料模型計算結果差異較大。

表3 冰體落體碰撞實驗主要參數Tab.3 Main parameters of ice drop impact test

圖13 圓錐冰落體碰撞實驗和數值模擬的力-位移曲線結果對比Fig.13 Comparison of the experimental and numerical force-displacement curves for conical ice drop impact

3 冰材料數值模型可行性評估方法

在船-冰碰撞下結構動力響應數值模擬中，冰材料數值模型的準確性與可行性評估是船-冰碰撞數值研究中的一個重要環節。目前國內外學者主要采用的評估方法有模型實驗驗證方法[20,24,26-27,35]、ISO 規范等推薦的壓力面積曲線（P-A）對比方法[13,27-28,31,40]以及比吸能（SEA，Specific Energy Absorption）判定方法[41-43]。模型實驗方法是評估數值模擬方法可行性的最好方法之一，可以通過實驗直接觀察和分析冰體失效過程，并獲得直觀的實驗數據，但是模型實驗成本較高，耗時較長，且其提供的驗證數據有實驗條件限制，具有較大的不確定性。為此，許多學者基于大量實驗測量得到的數據，對冰體碰撞過程中接觸壓力和接觸面積進行分析，通過擬合得到壓力-面積關系，從而為數值模型驗證提供參考標準。比如，Masterson 等[44]總結了一系列的試驗數據，提出了一個由許多研究人員根據現場測量數據擬合得到的P-A關系式，并將其納入ISO 規范[45]。Cai等[27]開展了冰體沖擊剛性墻的有限元數值計算，參考一些現有的有限元計算結果，總結了不同冰體材料模型下壓力-面積值，并與ISO/CD 19906[43]以及API/CSA[46]推薦的壓力-面積曲線進行對比，如圖14 所示。利用現有的一些基于實驗數據得到的P-A曲線，來對所提出的冰材料數值計算模型進行標定和校對，是目前比較常用的方法。

圖14 不同冰體材料模型P-A值與ISO規范推薦曲線對比[27]Fig.14 Comparison between P-A values calculated by different ice material models and P-A curve recommended by ISO rule[27]

對于一種破壞失效吸能材料，單位的質量材料破壞失效所吸收的能量值（SEA）常用來衡量材料的吸能特性。船-冰碰撞過程中冰體體現出一種較強的擠壓破壞失效特性，因此，SEA也用來表征冰體材料的吸能特性。圖15展示了不同縮尺比實驗測定得到的SEA 值，根據冰體破壞體積Vc，主要劃分為三個區域：區域1，當1.0E-8 m3≤Vc≤1.0E-7 m3時，SEA值的范圍大致在1.0E4 J/kg～1.0E5J/kg之間；區域2，當1.0E-7 m3≤Vc≤1.0E-5m3時，SEA值的范圍大致在1.0E3 J/kg～1.0E4J/kg之間；區域3，當1.0E-5m3≤Vc≤1.0E-3m3時，SEA值主要集中在1.0E4 J/kg 左右。在各自區域內，隨著冰體破碎體積的增大，SEA 值在逐漸減小。由于圖15 中總結的實驗數據所對應的模型實驗冰體尺度不同，所以展現出了SEA 值的三個劃分區域[42]。根據Cai等[27]文中的楔形冰與板模型碰撞實驗工況，主要碰撞參數如表2 工況1 所示，開展不同冰材料模型下有限元數值計算，得到不同破碎體積下的SEA 值，如圖16所示。由于受到模型實驗尺寸的限制，隨著碰撞時間的變化，冰體的破壞體積范圍主要集中在區域3。圖16 給出了不同冰材料模型下不同冰體破壞體積所對應的SEA數值計算值，可以發現：圖中的冰材料模型所計算得到的SEA值都小于實驗參考平均值10 000 J/kg，各向同性彈性失效模型（MAT_13）[13]和混凝土損傷模型（MAT_78）[27]計算得到的SEA 值隨著破壞體積的增大有下降的趨勢，而模擬高壓區的可壓碎泡沫型材料模型（MAT_63，M1）[22]計算得到的SEA 值隨著破壞體積的增大有著明顯的上升趨勢，對于模擬低壓區的可壓碎泡沫型材料 模 型（MAT_63，M2）[22]和 損 傷 模 型（MAT_153）[26]計算得到的SEA 值沒有很明顯的變化趨勢，相對比較穩定。除此之外，從表4 可以看出，各向同性彈性失效模型（MAT_13）和損傷模型（MAT_153）所計算得到的SEA 平均值相對于其它模型的SEA 平均值和實驗參考值較小。每種材料模型所定義的材料參數都有著各自對應的適用條件，與圖15 中的實驗所對應的條件有所區別，因此模型實驗的SEA值和數值模擬的SEA值有所差別。

圖15 不同縮尺比實驗測定的SEA值[42]Fig.15 SEA values measured by model tests with different scaled ratios[42]

圖16 不同冰材料模型下不同冰體破壞體積所對應的SEA數值計算值Fig.16 Numerical calculation values of SEA under different ice material models versus different ice crushed volumes

表4 不同冰材料模型下的SEA數值計算值Tab.4 Numerical calculation values of SEA under different ice material models(Experimental reference value:10 000 J/kg)

4 結 語

在自然界中，無論是海冰還是淡水冰，不同環境下冰體的形成機理各不相同，如在極地地區有1年冰和多年冰，不同年份的冰層的材料力學特性相差較大，此外，冰體一些主要的力學參數和力學性質還受應變率、溫度以及鹽度、雜質等冰的形成成分因素的影響。因此，對冰船作用有限元數值模擬方法來說，開發一種涵蓋所有條件下的冰體材料數值模型是目前冰工程領域里一項極具挑戰性的研究工作，其中建立準確的冰材料本構關系是其關鍵技術和難點問題，仍處于研究探索階段；除此之外，在船-冰碰撞數值模擬過程中，冰體變形破壞失效模式和裂紋擴展的真實模擬也是今后亟需解決的有限元數值模擬難點問題。

本文總結了一些現有冰材料有限元數值計算模型，分析了各自的特點，并選取了幾種模型實驗工況，開展了有限元仿真計算，對比并分析了不同冰材料模型計算結果的差異性；此外，還總結了一些現有的冰材料數值模型計算可靠性驗證方法，針對一些算例，討論并分析了用于冰材料模型驗證的P-A曲線驗證方法和比吸能（SEA）判定方法。從實驗結果與數值模型計算結果對比可知，每種材料模型都有各自的適用范圍及前提條件，目前所開發的冰材料數值模型很難適用于各種實驗工況，利用模型實驗驗證方法來對冰材料模型準確性和可行性進行判斷具有一定的局限性。為此，考慮溫度、應變率以及冰體形成條件等影響因素，建立不同實驗工況下的實驗結果數據庫，如SEA 值，為冰材料數值模型提供驗證數據，是目前船-冰碰撞數值方法研究的關鍵之處。