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

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

0 引 言

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

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

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

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

1 考慮溫度影響的海冰材料本構(gòu)模型

1.1 考慮溫度影響的海冰材料本構(gòu)理論

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

即

式（3）可以寫為

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

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

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

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

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

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

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

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

式中，dγ為塑性一致性參數(shù)增量。采用關(guān)聯(lián)流動法則，總應(yīng)變?yōu)閺椥詰?yīng)變與塑性應(yīng)變之和，

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

此外，還必須建立合理的單元失效準(zhǔn)則來模擬海冰材料的脆性破壞特征，海冰存在多種破壞形式，在受壓情況下剪切和擠壓破壞是兩種主要形式。但在實際碰撞問題的處理中，當(dāng)采用剪切破壞準(zhǔn)則時，實際上過高估計了材料的抗剪強度，往往會導(dǎo)致單元提前失效。這是因為冰在經(jīng)歷塑性階段的擠壓過程中，單元靜水應(yīng)力增加，此時單元的強度明顯增加，此種失效標(biāo)準(zhǔn)顯然更難達(dá)到?；谝陨戏治觯疚牟捎靡月薀o關(guān)和累積塑性應(yīng)變?yōu)楸碚鞯膭討B(tài)失效準(zhǔn)則來定義單元的失效，該準(zhǔn)則起源于Jordaan的經(jīng)驗失效方程[19]，如下式所示：

式中，a、p1為常數(shù)，p為靜水應(yīng)力。

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

1.2 子程序算法實現(xiàn)

積分率本構(gòu)方程的數(shù)值算法被稱為應(yīng)力更新算法或本構(gòu)積分算法，對于率無關(guān)塑性，常用的本構(gòu)積分算法有顯式算法、隱式算法以及半隱式算法。一個廣義的率無關(guān)的彈塑性本構(gòu)模型[20]可以寫為

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

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

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

圖3 本構(gòu)算法流程Fig.3 Flow chart of constitutive algorithm

1.3 單元測試和材料驗證

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

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

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

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

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

圖7給出了碰撞力-時間曲線。可以看到，在整個作用過程中，碰撞力呈現(xiàn)典型的非線性特征，剛開始接觸時，由于球體接觸區(qū)域面積變化較快，碰撞力上升明顯且上升速度較快，其后，面積變化速率逐漸減小，碰撞力上升趨勢變緩，總體趨于平穩(wěn)。計算得到的壓力-面積曲線如圖8 所示，可以看到，壓力-面積曲線總體趨勢與ISO 標(biāo)準(zhǔn)[21]以及其他學(xué)者實際測得的變化趨勢相同，初始階段由于接觸面積較小，其值小于標(biāo)準(zhǔn)值，且此時誤差較大，隨著面積增大，平均壓力值逐漸減小，與標(biāo)準(zhǔn)值逐漸接近，最后階段，壓力值趨于平緩，且與標(biāo)準(zhǔn)曲線的一致性較好，實際情況中，有圍壓下海冰強度最大約為30 MPa 左右，因此在使用中應(yīng)以面積較大（≥0.25 m2）時為準(zhǔn)。由于ISO 標(biāo)準(zhǔn)的壓力-面積曲線是基于極限狀態(tài)設(shè)計（ULS）以及偶然極限狀態(tài)（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）對角線節(jié)點的位置在0.2 s、0.4 s、0.6 s、0.8 s四個時間點各個節(jié)點的節(jié)點力?？梢钥吹剑S著作用過程的進(jìn)行，存在著高壓區(qū)與低壓區(qū)的相互轉(zhuǎn)化，越接近中心的節(jié)點力變化幅度較大，沿著四周的逐漸減小。

圖9 節(jié)點位置圖Fig.9 Nodal position

圖10 節(jié)點力Fig.10 Nodal force

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

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

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

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

2 船-冰碰撞場景仿真

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

2.1 場景介紹

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

圖14 碰撞場景Fig.14 Collision scenario

表1 DH32鋼的力學(xué)參數(shù)Tab.1 DH32 steel mechanical parameters

圖15 真實應(yīng)力-應(yīng)變曲線Fig.15 Real curves of stress-strain

2.2 結(jié)果分析

冰體、結(jié)構(gòu)的損傷變形如圖16 所示?？梢钥吹剑鲎策^程中冰體在與船艏接觸區(qū)域發(fā)生較大程度的破壞，艏部結(jié)構(gòu)也發(fā)生較大程度變形。冰體最大應(yīng)力達(dá)到15 MPa 左右，破碎表面不規(guī)則并伴有碎片產(chǎn)生，冰體有明顯的高低應(yīng)力區(qū)且高應(yīng)力區(qū)分布具有顯著的局部性。由于碰撞區(qū)域較小，速度較大，艏部結(jié)構(gòu)最大應(yīng)力達(dá)到764 MPa 左右，外板變形呈現(xiàn)連續(xù)凹陷的形式，內(nèi)部構(gòu)件交叉區(qū)域產(chǎn)生一定程度的應(yīng)力集中現(xiàn)象，應(yīng)力分布區(qū)域較為集中，可見球鼻艏內(nèi)部加強結(jié)構(gòu)在碰撞過程中起到重要作用。

圖16 冰體與結(jié)構(gòu)損傷變形Fig.16 Damage and deformation of the ship structure and ice block

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

圖17 碰撞力-時間關(guān)系曲線Fig.17 Collision force-time curve

3 結(jié) 論

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

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

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

（3）船-冰碰撞會對船體結(jié)構(gòu)產(chǎn)生顯著的變形和破壞，碰撞過程中冰體在與船艏接觸區(qū)域發(fā)生較大程度的破壞時，艏部結(jié)構(gòu)也發(fā)生較大程度變形。冰體的高應(yīng)力分布均有明顯的局部性，內(nèi)部加強結(jié)構(gòu)在碰撞過程中起到重要作用。