船-冰碰撞中考慮溫度影響的冰體材料本構模型研究

俞同強，劉 昆，劉俊杰，王自力

（1.江蘇科技大學船舶與海洋工程學院，江蘇鎮江 212003；2.中國船舶科學研究中心，江蘇無錫 214082）

0 引 言

隨著全球氣候變暖，極地地區的冰雪消融日益加劇，北極海冰尤其是夏季海冰正以每十年10%的速度減少[1]，據估計，北冰洋十年內將出現夏季無冰年，北極航道的全面開通成為可能。對于我國來說，利用北極東北航道，從上海以北港口到歐洲西部港口的航程比傳統航程縮短25%～55%[2]，不僅可以大大減輕貿易成本，而且可以擺脫對馬六甲海峽的依賴，同時對于我國漁業、旅游業以及油氣開采等行業具有深遠影響，其重要性不言而喻。

北極航道的通航對極地船舶的安全性提出了更高的要求，船-冰碰撞所帶來的結構強度與冰載荷問題日益受到人們的重視。由于極地海域復雜的環境特征，即使在夏天也有浮冰存在，可能導致船舶破損、進水甚至沉沒，造成重大的生命財產損失。即使輔以現代設備，船-冰碰撞問題仍然難以避免。2019年“雪龍”號在執行我國第35次南極考察任務期間，在阿蒙森海密集冰區航行中，因受濃霧影響，在南緯69°59.9'、西經94°04.2'位置與冰山碰撞，船艏桅桿及部分舷墻受損，所幸無人員受傷?？梢?，開展極地船舶船-冰碰撞性能研究，對于確定冰載荷大小、評估船體結構性能、減少船-冰碰撞的事故損失以及對未來極地船舶的結構設計等均具有非常重要的意義。

研究船-冰碰撞問題首先需要確定冰載荷的大小，目前主要有理論計算、試驗以及仿真三種方法。理論研究方法主要是在一定尺度下，通過能量法等力學理論結合單軸、壓痕試驗數據，得到不同工況下冰載荷的經驗計算公式[3-5]，方法簡單但是精度較低，適用范圍不夠廣泛。隨著計算機技術快速發展，數值算法不斷進步，商業有限元軟件被越來越多地用于解決海冰與結構作用時冰材料的模擬和冰載荷的計算問題，此時船-冰碰撞問題的研究已發展到一個新的階段。但是，由于海冰材料復雜的物理和力學特性，目前對海冰以及冰載荷的研究仍然面臨諸多困難，其中對其本構關系的合理準確描述一直是影響相關問題計算分析的關鍵因素。國內外學者相繼提出粘塑性、粘彈塑性、彈塑性、彈脆性和泡沫材料等多種材料模型[6-10]，具有一定的實用性但也有各自的局限性。目前由于商業有限元軟件中缺乏海冰材料相關定義，使用自定義材料本構來模擬海冰力學特性和計算冰載荷成為一種切實、有效的手段[11-12]。

本文建立一種帶脆性失效特征的彈塑性冰材料本構模型，考慮溫度對屈服函數及失效準則的影響，編寫相應材料子程序并通過二次開發嵌入到LS-DYNA 材料庫，以此為基礎開展船-冰碰撞研究，分析船體與海冰相互作用時冰體和結構的損傷和變形特點。本文研究成果可為極地船舶碰撞損傷評估和結構設計制造提供參考。

1 考慮溫度影響的海冰材料本構模型

1.1 考慮溫度影響的海冰材料本構理論

自然界中的海冰是由結晶冰、鹵水、空氣和雜質等組成的復雜混合物，其成分比例隨外界條件、時間和空間的變化而變化，其中最重要的控制因素是溫度，所以從材料學角度來看，海冰可以視為溫度敏感性復合材料[13]。隨著海冰力學試驗尤其是多軸試驗的廣泛實施，海冰三維應力狀態和材料力學特點得到全面直觀反映[14-17]，為理論模型建立和數值仿真提供了有效的數據支撐。其中，Derradji-Aouat[18]在總結歸納三軸試驗結果的基礎上給出了各向同性淡水冰、冰山冰以及柱狀冰的多曲面橢圓屈服方程，該屈服方程得到了廣泛接受和使用。

即

式（3）可以寫為

該屈服方程即為Tsai-Wu 屈服準則，其中，f( )p,J2為屈服函數，J2為第二偏應力不變量，a0、a1、a2為系數。Tsai-Wu 準則能夠體現出復雜應力狀態下靜水壓力能影響材料屈服的性質，同時十分適用于描述對于拉壓強度相差很大的海冰材料。該屈服方程在主應力空間下的屈服曲面及在坐標軸上的投影如圖1所示，可以看到在三維空間下，屈服面呈橢球型，在坐標面上的投影為橢圓形，且由于未考慮各向異性，各坐標面下屈服面投影形狀保持一致，對于同一類型冰，屈服面在溫度等因素影響下膨脹或收縮，同時保持橢球焦心不變，即λ、pc、η的數值保持不變。

圖1 主應力空間下的屈服曲面及其在坐標軸上的投影Fig.1 Yield surface in principal stress space and its projection on coordinate axis

由以上分析可知，參數a0、a1、a2的值僅與參數qmax有關，在三軸試驗中，qmax=σ1-σ3，因此可以較為系統地總結出qmax與試驗中各變量的關系，Derradji-Aoua 給出式（5）～（6）來表征qmax與應變率和溫度的關系。

船-冰碰撞下冰體通常處于高應變率（10-3以上）范圍，在此種場景下海冰表現為彈脆性，應變率的影響可以不計。根據Gagnon 等[15-17]的三軸試驗數據，針對該應變率附近不同溫度下的試驗數據進行統計、擬合，得到的屈服曲線如圖2所示。

因此，參數a0、a1、a2可以表示為溫度的函數ai(T)，將該屈服方程發展為受溫度T影響的多重屈服面準則作為本文海冰本構的屈服準則：

圖2 不同溫度下屈服方程曲線Fig.2 Yield equation curves at different temperatures

該本構模型在彈性階段，滿足各向同性廣義胡克定律，其增量形式為

式中，剪切彈性模量G=E/2( 1 +μ， )E為彈性模量，μ為泊松比，γj為工程切應變，γj= 2εj(j = xy,yz,zx )，Δε0為靜水應變增量，Δσ0為靜水應力增量。在塑性階段，由于屈服面的外凸性，塑性應變增量為

式中，dγ為塑性一致性參數增量。采用關聯流動法則，總應變為彈性應變與塑性應變之和，

式中：dσij為應力增量；δij為Kronecker符號，其值為0（i≠j）或1（i=j）。

此外，還必須建立合理的單元失效準則來模擬海冰材料的脆性破壞特征，海冰存在多種破壞形式，在受壓情況下剪切和擠壓破壞是兩種主要形式。但在實際碰撞問題的處理中，當采用剪切破壞準則時，實際上過高估計了材料的抗剪強度，往往會導致單元提前失效。這是因為冰在經歷塑性階段的擠壓過程中，單元靜水應力增加，此時單元的強度明顯增加，此種失效標準顯然更難達到?；谝陨戏治?，本文采用以率無關和累積塑性應變為表征的動態失效準則來定義單元的失效，該準則起源于Jordaan的經驗失效方程[19]，如下式所示：

式中，a、p1為常數，p為靜水應力。

本文進一步發展該經驗失效準則，具體為：考慮溫度的影響，將p1設置為屈服方程較大的一個根，同 時，給 定 初 始 失 效 應 變ε0= 0.01，當＞εf時 材 料 失 效，其 中為 等 效 塑 性 應 變=，εf為失效應變，εf=ε0+(p/p1- 0.5)2。此外，考慮到拉壓強度的不同，給定拉伸強度極限，當p＜pc時，單元刪除。pc為截斷應力，設為拉伸極限2 MPa。該失效準則特點是注重由靜水應力帶來的高應力區和低應力區的承載能力差別，同時抑制摩擦滑移，在碰撞問題中弱化剪切破壞的影響。

1.2 子程序算法實現

積分率本構方程的數值算法被稱為應力更新算法或本構積分算法，對于率無關塑性，常用的本構積分算法有顯式算法、隱式算法以及半隱式算法。一個廣義的率無關的彈塑性本構模型[20]可以寫為

式中：εij分別為總應變、彈性應變以及塑性應變；q為內變量；hα為塑性模量；rij為屈服函數的偏導數為塑性應變增量；γ為塑性參數，由加卸載準則確定。本文采用最近點投影法來求解，其為一種半隱式的返回映射算法。計算中對塑性參數采用隱式而對塑性流動方向和塑性模量采用顯式算法，在步驟結束時計算塑性參數的增量并同時強化屈服條件。在某時刻t的下一個時間步長Δt中，試驗應力由胡克定律求得，若在彈性階段內，則直接更新應力，若超過屈服極限，則進入迭代步，其各參數計算公式為

式中，Δγn+1即為塑性參數γ在n+1 步的增量，C為彈性系數矩陣。在每一個積分步內，首先計算方程塑性參數值與應力分量，

每個迭代步內檢查收斂條件，若收斂，退出迭代算法。由此可得n+1 步的塑性應變增量為

圖3 本構算法流程Fig.3 Flow chart of constitutive algorithm

1.3 單元測試和材料驗證

建立單元測試模型分別施加恒定壓縮速度，測試溫度為-10 ℃，考察單元變形中相關參數的變化情況，在子程序計算中，輸出第二偏應力不變量J2、靜水應力p以及失效應變的值，得到J2-p、εf-p的曲線分別如圖4～5 所示?？梢钥吹?，偏應力不變量J2與靜水應力P呈二次曲線關系，其中標準曲線為理論曲線。圖4 中在a～b范圍內，輸出值低于標準值，此時單元在彈性范圍內，即屈服應力f(p,J2,T)＜0。當到達b點之后，兩線重合，單元達到屈服狀態，輸出值與理論值一致，初始屈服點在10 MPa 左右，與理論解符合。失效應變與靜水應力也呈二次曲線關系，靜水應力輸出值與理論值一致，當靜水應力為55 MPa左右時，失效應變達到最小值為0.01，與理論解相符，且能合理表現冰材料的脆性特征，子程序算法的執行過程是正確的。

圖4 J2-p曲線圖Fig.4 Relationship curve of J2-p

圖5 εf-p曲線圖Fig.5 Relationship curve of εf-p

同時，建立剛性板-球形冰體擠壓模型，通過壓力-面積曲線驗證材料的準確性。球體直徑為1 m，施加1 m/s 的恒定速度，0.5 s 時球體的損傷變形如圖6 所示?？梢钥吹角蛐伪w表面發生破壞失效，失效面不平整，接觸邊緣有冰體脫落，冰體接觸面內應力分布不均，有明顯的高應力區域和低應力區域，最大應力為23.5 MPa左右，表明在三維應力下，材料的強度得到了一定程度的增強。

圖6 0.5 s時球形冰體的破壞情況Fig.6 Damage of ice blocks at 0.5 s

圖7給出了碰撞力-時間曲線?？梢钥吹?，在整個作用過程中，碰撞力呈現典型的非線性特征，剛開始接觸時，由于球體接觸區域面積變化較快，碰撞力上升明顯且上升速度較快，其后，面積變化速率逐漸減小，碰撞力上升趨勢變緩，總體趨于平穩。計算得到的壓力-面積曲線如圖8 所示，可以看到，壓力-面積曲線總體趨勢與ISO 標準[21]以及其他學者實際測得的變化趨勢相同，初始階段由于接觸面積較小，其值小于標準值，且此時誤差較大，隨著面積增大，平均壓力值逐漸減小，與標準值逐漸接近，最后階段，壓力值趨于平緩，且與標準曲線的一致性較好，實際情況中，有圍壓下海冰強度最大約為30 MPa 左右，因此在使用中應以面積較大（≥0.25 m2）時為準。由于ISO 標準的壓力-面積曲線是基于極限狀態設計（ULS）以及偶然極限狀態（ALS）方法，壓力測量取極限值，有一定的富裕度。本文給出Timco 和Molikpaq[22-23]通過試驗測得的壓力-面積曲線（試驗溫度-5 ℃～10 ℃），相比之下，本文所得壓力-面積曲線在面積較小時與試驗得到的曲線吻合得更好。因此，該材料模型能夠反映海冰強度的真實變化情況。

圖7 碰撞力-時程曲線Fig.7 Force-time curve

圖8 壓力面積曲線Fig.8 Pressure-area curve

圖9 和圖10 給出了碰撞過程中剛性板局部（1 m×1 m）對角線節點的位置在0.2 s、0.4 s、0.6 s、0.8 s四個時間點各個節點的節點力?？梢钥吹?，隨著作用過程的進行，存在著高壓區與低壓區的相互轉化，越接近中心的節點力變化幅度較大，沿著四周的逐漸減小。

圖9 節點位置圖Fig.9 Nodal position

圖10 節點力Fig.10 Nodal force

考慮溫度的影響，約束單元底端的四個節點，頂端施加恒定的壓縮速度，直至單元失效，輸出不同溫度下失效時單元的等效應力和等效應變。不同溫度下單元失效時的等效應力和應變以及計算所得壓力-面積曲線如圖11～13所示?？梢钥吹?，隨著溫度的不斷降低，單元失效時的應力增加，表明材料強度不斷增加，材料逐漸“變硬”，同時失效應變逐漸減小，材料逐漸“變脆”，與實際情況相符。不同溫度下的壓力-面積曲線變化趨勢相同，壓力趨于平穩時，在較低溫度下的壓力值較高。

圖11 單元失效時應力-溫度曲線Fig.11 Failure stress-temperature curve in single unit test

圖12 單元失效時應變-溫度曲線圖Fig.12 Failure stain-temperature curve in single unit test

通過以上驗證，可以表明本文提出的自定義冰材料模型子程序執行正確且能夠較好地模擬海冰材料特性和冰載荷，合理體現其溫度敏感性特性，可運用于實際的碰撞場景中。

2 船-冰碰撞場景仿真

圖13 不同溫度的壓力-面積曲線Fig.13 Pressure-area curve at different temperatures

2.1 場景介紹

通過調研分析確定典型計算工況，對于船艏部位的撞擊，選取船舶排水區域內船艏典型結構進行研究。本文選取艏部結構中的球鼻艏結構與塊狀海冰的碰撞場景，如圖14 所示。船體艏部選取球鼻艏至艏部艙壁向后6.5 m處，尺寸為20 m×48 m×26 m，建模采用4 積分點、5 節點殼單元，網格尺寸為200 mm，船體結構所用材料為DH32型低溫鋼，材料屬性由低溫準靜態拉伸試驗所得，參數如表1所示，其真實應力-應變曲線如圖15 所示。塊狀冰體尺寸為10 m×10 m×5 m，采用實體單元，網格尺寸為100 mm。目前極地商船的航速一般不超過16 kn[24]，本文中碰撞相對速度為8 m/s，溫度為-10 ℃。

圖14 碰撞場景Fig.14 Collision scenario

表1 DH32鋼的力學參數Tab.1 DH32 steel mechanical parameters

圖15 真實應力-應變曲線Fig.15 Real curves of stress-strain

2.2 結果分析

冰體、結構的損傷變形如圖16 所示?？梢钥吹?，碰撞過程中冰體在與船艏接觸區域發生較大程度的破壞，艏部結構也發生較大程度變形。冰體最大應力達到15 MPa 左右，破碎表面不規則并伴有碎片產生，冰體有明顯的高低應力區且高應力區分布具有顯著的局部性。由于碰撞區域較小，速度較大，艏部結構最大應力達到764 MPa 左右，外板變形呈現連續凹陷的形式，內部構件交叉區域產生一定程度的應力集中現象，應力分布區域較為集中，可見球鼻艏內部加強結構在碰撞過程中起到重要作用。

圖16 冰體與結構損傷變形Fig.16 Damage and deformation of the ship structure and ice block

圖17 給出了碰撞過程中碰撞力時程曲線?？梢钥吹?，總體上碰撞力呈先迅速上升后逐漸趨緩的特征，在整個作用過程中碰撞力曲線表現出顯著的非線性特征。在碰撞初始階段由于接觸面積迅速增大，碰撞力上升速度較快，隨著碰撞的進行，冰體單元存在同時大量失效的情況，碰撞力有瞬間減小的情況，此后接觸面積增長速度逐漸變緩，碰撞力增長速度減慢，在最后階段，面積趨于最大值，碰撞力變化趨于平穩略有上升。

圖17 碰撞力-時間關系曲線Fig.17 Collision force-time curve

3 結 論

本文通過LS-DYNA 材料子程序二次開發，建立并驗證了基于多曲面屈服準則的考慮溫度影響的海冰材料本構模型，并開展了船-冰碰撞相關研究，主要結論如下：

（1）基于多曲面屈服準則的自定義海冰材料本構模型能夠較好地展現冰材料的材料特性，合理地反映其力學特性，得到的壓力-面積曲線與ISO 標準及相關實驗結果相符，該本構模型可運用于實際的碰撞場景中。

（2）溫度對于海冰材料性能有直接影響，隨著溫度的降低，海冰材料強度不斷增加，材料逐漸“變硬”，同時失效應變逐漸減小，材料“變脆”，與實際情況相符，海冰材料本構模型合理地體現其溫度敏感性特點。

（3）船-冰碰撞會對船體結構產生顯著的變形和破壞，碰撞過程中冰體在與船艏接觸區域發生較大程度的破壞時，艏部結構也發生較大程度變形。冰體的高應力分布均有明顯的局部性，內部加強結構在碰撞過程中起到重要作用。