胡志強，高 巖，姚 琪

（1. 上海交通大學(xué)海洋工程國家重點實驗室，上海 200240；2. 上海交通大學(xué)高新船舶與深海開發(fā)裝備協(xié)同創(chuàng)新中心，上海 200240）

一種理想彈塑性模擬的冰材料本構(gòu)模型

胡志強1，2，高 巖1，姚 琪1

（1. 上海交通大學(xué)海洋工程國家重點實驗室，上海 200240；
2. 上海交通大學(xué)高新船舶與深海開發(fā)裝備協(xié)同創(chuàng)新中心，上海 200240）

船冰碰撞是船舶碰撞研究領(lǐng)域的熱點之一，對冰材料的模擬是船冰碰撞的研究重點。提出一種利用理想彈塑性模型模擬的冰材料本構(gòu)模型，利用半隱式圖形算法計算單元塑性階段的應(yīng)力，利用Tsai-Wu屈服準則和經(jīng)驗失效公式用來描述冰的力學(xué)行為。利用二次開發(fā)功能，將冰材料模型嵌入LS_DYNA程序，并驗證該模型的準確性和適用性。研究中針對不同局部形狀的冰塊與船側(cè)碰撞場景，通過比較分析碰撞力、能量耗散等，探討冰塊的局部形狀對碰撞場景的影響。研究結(jié)果表明：冰材料模型在大接觸面的條件下壓力與已有標準吻合較好；在不同的冰塊局部形狀條件下，船冰碰撞的相互作用過程不同；較鈍形狀的冰塊表現(xiàn)近乎剛體，較尖銳形狀的冰塊較易破碎。

船冰碰撞；彈塑性材料；本構(gòu)模型；冰塊形狀影響；數(shù)值仿真

0 引 言

近年來，隨著在極地地區(qū)航行的船舶日益增多，船舶與冰山發(fā)生碰撞事故的可能性也隨之增加。船舶一旦與冰山碰撞造成破損，很可能發(fā)生漏油事故，影響極地地區(qū)的生態(tài)環(huán)境。因此，為了保證航線的安全性，提高船舶結(jié)構(gòu)的可靠性，需要在船舶結(jié)構(gòu)設(shè)計階段考慮冰載荷的影響。因此，準確模擬船冰碰撞場景，預(yù)測碰撞結(jié)果，為船舶設(shè)計階段提供數(shù)據(jù)參考非常重要。

冰的力學(xué)性質(zhì)十分復(fù)雜。冰的形成過程高度依賴于環(huán)境條件，如海水鹽度、溫度等，因此不同區(qū)域的海冰的性質(zhì)不同。如在俄羅斯北部，因有大量的河流入海，海冰的種類和性質(zhì)與北美北極圈的Beaufort海就有很大的不同［1］。從宏觀組成上，海冰由固體冰、鹽囊、氣體和孔隙等組成，其物理性質(zhì)相當(dāng)復(fù)雜。在微觀層面上，冰的性質(zhì)決定于晶體大小及其結(jié)構(gòu)。粒狀冰、柱狀冰和不連續(xù)冰為三種基本冰晶結(jié)構(gòu)［2］。如果組成海冰的晶體結(jié)構(gòu)不同，或者排列方式不同，海冰的力學(xué)性質(zhì)也會不同。比如多年冰晶體排列無明顯方向性，可近似為各向同性材料。一年冰由于其柱狀冰的方向性，只能簡化為正交異性材料。在宏觀方面，冰的厚度、孔隙率等是冰失效條件的重要影響因素，由于冰在其環(huán)境中較高的同系溫度（Th＞0.8），它的力學(xué)行為與應(yīng)變率有關(guān)。當(dāng)應(yīng)變率較低時，冰主要表現(xiàn)為延性破壞，軟化行為是此破壞模式的最明顯的特征。然而，冰在高應(yīng)變率條件下會表現(xiàn)出明顯的脆性破壞，在載荷峰值之后破壞，徹底失去承受載荷的能力。冰強度的最大值出現(xiàn)在兩者的過渡區(qū)域，應(yīng)變率≈0.001s-1［3］。所以，在低應(yīng)變率條件下，冰可近似為非線性的粘彈塑性材料，同時伴有延性失效模式的材料［4］。在高應(yīng)變率條件下，冰可看作線彈性并脆性失效模式材料。在脆性失效模式下，冰的壓縮強度遠大于拉伸強度［5］。鑒于冰材料的復(fù)雜性，完全模擬其在碰撞過程中的物理變化過程是相當(dāng)困難的。目前對冰材料的模擬主要集中在冰載荷模型上，重點模擬其壓力-面積曲線特性。

許多學(xué)者提出了不同的冰材料本構(gòu)模型。Gagnon［6］根據(jù)冰塊與船舶模型的實驗現(xiàn)象，提出了泡沫模型模擬冰材料本構(gòu)模型，可有效解釋在試驗中觀測到的冰在碰撞區(qū)域的融化現(xiàn)象。Jordaan［7］則從微觀層面描述了冰的破碎過程，利用隨機概率的方法提出了兩種壓力-面積曲線描述冰載荷特性。Jebaraj［8］利用有限元的方法模擬船冰碰撞過程，但是其冰材料模型比較簡單，需要進一步完善。季順迎、岳前進［9］則總體介紹了冰的數(shù)值模型，針對性的研究分析了渤海海冰的相關(guān)性質(zhì)。本文提出了一種基于彈塑性模擬的冰材料本構(gòu)模型，并嵌入到LS_DYNA程序，模擬海冰材料特性，應(yīng)用于船冰碰撞數(shù)值仿真計算，進而探討冰的局部形狀對冰的力學(xué)行為的影響。

1 冰的彈塑性材料特性模擬

1.1彈塑性材料理論

冰材料模型是影響有限元模擬結(jié)果準確性的重要因素之一。但是，由于冰材料的復(fù)雜性，難以模擬船冰碰撞過程中冰材料的所有行為特性。因此，目前的研究主要集中在模擬碰撞過程中的冰載荷特性。Ralston［10］研究了冰破碎的屈服條件和塑性變形，為用塑性理論研究冰材料證明了其可行性。本文利用理想彈塑性材料模擬冰的載荷特性，進而研究船冰碰撞中兩者的相互作用過程。

理想彈塑性是塑性本構(gòu)關(guān)系的特殊形式之一，其在塑性階段不產(chǎn)生硬化，即硬化函數(shù)為零。在彈性階段應(yīng)力應(yīng)變滿足廣義胡克定律。圖 1為一維條件下的理想彈塑性本構(gòu)關(guān)系，冰材料的本構(gòu)關(guān)系為一維形式在三維空間的泛化。

在彈性階段，應(yīng)力應(yīng)變滿足廣義胡克定律，如公式（1）所示。

圖1　理想彈塑性本構(gòu)關(guān)系示意

式中：E——彈性模量；μ——Poisson比；G——剪切彈性模量，且）——正應(yīng)變；）——工程切應(yīng)變，且

其彈性本構(gòu)關(guān)系增量形式為：

屈服方程是判斷彈性狀態(tài)和塑性狀態(tài)的準則。因此，一個較為準確的屈服方程對準確模擬冰材料尤為重要。一般情況下，在碰撞過程中，碰撞面積處的冰單元受周圍單元的高度限制，處于典型的三向應(yīng)力狀態(tài)。Jones［11］和Rist［12］的實驗結(jié)果表明，對于冰材料，靜水壓力對屈服狀態(tài)的影響是不可忽略的。冰的壓縮強度隨著靜水壓力的增加而增加。因此采用Tsai-Wu［13］屈服準則，如公式（4）所示，即第二偏應(yīng)力不變量為靜水應(yīng)力的二次函數(shù)。Tsai-Wu 屈服是以300個實驗數(shù)據(jù)為基礎(chǔ)提出的，在船冰碰撞研究中被廣泛采用。其函數(shù)為：

式中：J2——第二偏應(yīng)力不變量；P——靜水應(yīng)力。a0，a1，a2——分別為常數(shù)，表1中所示為常數(shù)的取值。本文采用Derradji-Aouattu推薦的取值。

表1　常數(shù)取值

提出的理想彈塑性材料在塑性階段采用增量理論本構(gòu)關(guān)系。根據(jù)Drucker公設(shè)中的假設(shè)，應(yīng)變分為彈性應(yīng)變與塑性應(yīng)變，則應(yīng)變增量也分為彈性應(yīng)變增量與塑性應(yīng)變增量，如公式（5）所示，其中——彈性應(yīng)變，——塑性應(yīng)變。

彈性變形規(guī)律不因塑性變形而改變，應(yīng)力仍然符合廣義Hooke定律，即

塑性應(yīng)變增量方向與加載面的外法線方向重合。

式中：dγ——一非負的塑性一致性參數(shù)，塑性加載時dγ＞0，中性變載和卸載時dγ=0。

如果采用關(guān)聯(lián)流動法則，φ則為屈服函數(shù)，即Tsai-Wu 屈服函數(shù)。

所以，塑性階段的本構(gòu)關(guān)系可表示為式（7）。彈性階段的應(yīng)力可用Hook定律求出解析解，塑性階段的應(yīng)力由于塑性一致性參數(shù)的不確定性，需利用算法求取其數(shù)值解。

理想彈塑性本構(gòu)關(guān)系的關(guān)鍵問題之一為失效條件的確定。本文采用Zhenhui Liu［14］提出的以等效塑性應(yīng)變和靜水應(yīng)力為基礎(chǔ)的經(jīng)驗失效準則，如公式（8）：

1.2塑性流動理論算法實現(xiàn)

對于塑性單元或彈塑性過渡單元，其應(yīng)力的計算十分復(fù)雜。常用的應(yīng)力更新積分算法有：顯式和隱式圖形返回法。表2為不同應(yīng)力更新積分算法的優(yōu)缺點比較，本文采取半隱式圖形返回法［15］進行迭代求解。

表2　不同本構(gòu)積分算法比較

根據(jù)所提到的基本理論，第n+1步的應(yīng)力可以表示成公式（9）的形式

因此可以把半隱式圖形返回法分成兩步：初始的彈性預(yù)測步，以及由于應(yīng)力對屈服表面產(chǎn)生偏離所需的塑性調(diào)整步，如式（10）所示。

同時，此種算法在步驟結(jié)束時強化屈服條件，避免了屈服面的漂移，對塑性一致性參數(shù)采用隱式，而對塑性流動方向和塑性模量采用顯示。積分方法可歸結(jié)為：

因此程序的第一部分為求解試驗應(yīng)力σtrial，其應(yīng)力應(yīng)變關(guān)系滿足 Hooke定律，隨之代入屈服函數(shù)進行判別：

如在彈性階段則σtrial=σ，即試驗應(yīng)力為真實應(yīng)力。如在塑性計算，則進行塑性修正，使試驗應(yīng)力回到屈服面上，見圖2。

圖2　圖形返回算法幾何示意

式中：Eσ——彈性空間； ?Eσ——彈性空間界面，即屈服曲面；——由εn計算出的試驗應(yīng)力；σn+1——真實應(yīng)力。

1.3編程過程

步驟1：初始化

步驟2：計算應(yīng)力和屈服函數(shù)

步驟4：計算塑性一致性參數(shù)增量

步驟5：更新狀態(tài)變量和一致性參數(shù)

步驟6：令K=K+1利用更新后的塑性應(yīng)變和一致性參數(shù)代入步驟2進行計算。

2 冰材料本構(gòu)模型驗證

首先利用單元測試驗證冰單元的材料性質(zhì)。采用共8個節(jié)點1個積分點的單元，其一側(cè)的4個節(jié)點固定，另一側(cè)的4個節(jié)點施加0.05m/s的速度以壓縮單元。圖3 （a）、3 （b）為利用單元實驗得到的屈服曲線和失效應(yīng)力。冰單元的屈服方程表現(xiàn)為拋物線型，與“1”中表達一致。屈服應(yīng)變最小值為0.1，發(fā)生在靜水應(yīng)力為50MPa處。一旦應(yīng)變達到失效應(yīng)變，單元失效并被刪除。圖3（c）為改變速度方向進行的拉伸試驗得到的截斷應(yīng)力，當(dāng)拉伸壓力達到2MPa之后則不再增加，從而確定了材料的拉伸性能。為與圖3（a）、圖3（b）中的壓縮方向的靜水應(yīng)力區(qū)分開，圖3（c）中的拉伸方向的靜水應(yīng)力用負值表示。通過與試驗結(jié)果的比較分析，驗證了本文提出的冰材料本構(gòu)模型的準確性。

圖3　試驗曲線分析

以球形冰-鋼板計算工況為例，研究冰本構(gòu)模型失效破碎的相關(guān)性質(zhì)，并驗證冰材料本構(gòu)模型的適用性。球形冰的半徑為1m，碰撞速度為1m/s，碰撞時間設(shè)定為0.5s。冰的壓力-面積曲線是冰力學(xué)的基本特性之一，被廣泛應(yīng)用于計算船冰相互作用過程中的冰載荷分析。Masterson［16］提出了這一壓力-面積公式，并且被ISO（19906）［17］采用。測定冰在相互作用過程的壓力-面積曲線，并與ISO規(guī)范進行對照驗證，結(jié)果如圖4（b）所示。由該圖可知，兩者在碰撞面積＞0.3m2的條件下吻合較好，而船冰碰撞的接觸面積一般大于此值，因此本文提出的冰本構(gòu)模型可較為準確的應(yīng)用到船冰碰撞數(shù)值仿真中。圖4（a）為在結(jié)束時刻的冰壓力云圖。在碰撞區(qū)域存在高壓區(qū)和低壓區(qū)，并且隨著時間推進，兩者相互轉(zhuǎn)換，高壓區(qū)的靜水應(yīng)力可達到77.3MPa，低壓區(qū)的應(yīng)力在2MPa左右，此靜水應(yīng)力范圍與Gagnon的實驗數(shù)據(jù)較為一致。圖4（c）為碰撞力-時間曲線，隨著碰撞的深入，碰撞力成鋸齒形增加。如果冰單元失效，則不再對碰撞貢獻力，碰撞力下降，直到一個新的較大的碰撞面積形成，更多單元參與碰撞過程，碰撞力才進一步增加。此鋸齒形狀為船冰碰撞力曲線的典型形狀。

圖4　球形冰-鋼板碰撞分析

3 利用冰材料本構(gòu)模型作船冰碰撞數(shù)值仿真分析

3.1碰撞場景模型及參數(shù)設(shè)置

不同的冰局部形狀可能導(dǎo)致不同的碰撞場景，構(gòu)建了5種形狀的海冰模型，模型單元與幾何尺寸見圖5。冰單元采用一個積分點的體單元，大小控制在50mm×50mm×50mm左右。在冰塊模型后面加一層與之共享節(jié)點的剛體單元，用來模擬冰塊其他部分對此局部的作用力，同時給以剛體節(jié)點V=2m/s的速度推進冰模型與船側(cè)相撞［18］。碰撞位置選擇為船外板的中心位置處。為了節(jié)省計算時間，只選取一段船側(cè)模型見圖5（6），船舶主尺度見表4。艙段模型長35m，高26m，船側(cè)縱桁布置間隔為900mm，船側(cè)肋板厚21mm，每隔5000mm布置一根。不作針對性的冰區(qū)加強。在船側(cè)模型兩側(cè)固定位移和轉(zhuǎn)動，用來模擬船舶其他部分對此艙段的邊界條件作用。船體結(jié)構(gòu)單元采用Belystchko-Tsay 殼單元，平均大小為225mm。

圖5　船側(cè)與冰塊局部形狀模型

表3　冰塊模型幾何參數(shù) 單位：m

表4　船舶主尺度 單位：m

表5　冰與船體的材料參數(shù)

3.2計算結(jié)果

限于篇幅，選取代表性結(jié)果顯示碰撞仿真結(jié)果。圖6和圖7分別為三棱柱冰塊與球狀冰塊在船冰碰撞過程中的破碎情況；圖8為球體冰塊碰撞過程船側(cè)的變形情況及應(yīng)力云圖；圖9為不同冰塊形狀的碰撞力和能量曲線。由圖6和圖7可知，隨著碰撞的持續(xù)，三棱柱冰塊破碎狀況非常明顯，大量冰單元失效。相反，球體冰仍保持原始形狀，較少冰單元失效刪除?？梢钥闯觯簩τ谳^鈍的冰塊，船側(cè)為能量耗散的主體，冰內(nèi)能/總內(nèi)能在較低水平；對于較為尖銳的冰塊，冰塊為能量耗散的主體，冰內(nèi)能/總內(nèi)能維持在一較高水平（見圖9）。相比于其他形狀，由于擁有較大的碰撞面積，方形冰的碰撞力和總內(nèi)能最大，最大碰撞力可達到58.9MN，最大總內(nèi)能接近70MJ；相對于其他形狀，球體冰表現(xiàn)得最為堅硬。船側(cè)外板在碰撞過程中未被穿透，但發(fā)生較大的塑性變形，船側(cè)縱桁和船側(cè)縱骨變形也較大，并且在碰撞區(qū)域發(fā)生斷裂現(xiàn)象。

圖6　三棱柱冰塊t=0.5s和t=1s破碎情況

圖7　球狀冰塊在t=0.5s和t=1s的破碎情況

圖8　球體冰碰撞中船側(cè)變形情況與應(yīng)力分布

圖9　球狀冰與三棱柱狀冰的碰撞力及能量曲線

4 結(jié) 語

提出一個利用理想彈塑性本構(gòu)關(guān)系的冰材料模型，并嵌入到LS_DYNA程序中，應(yīng)用到船冰碰撞數(shù)值仿真中，研究冰塊在不同局部形狀條件下與船側(cè)碰撞的特性，著重分析在碰撞過程中的碰撞力和能量耗散情況，為船舶在冰載荷作用下的結(jié)構(gòu)設(shè)計提供參考。研究得到以下主要結(jié)論：

1） 提出的理想彈塑性材料本構(gòu)模型，可較合理地模擬冰的受擠壓失效的行為特性；該模型可以較方便地嵌入到數(shù)值仿真程序中；

2） 不同局部形狀的冰塊在碰撞過程中表現(xiàn)出的力學(xué)性質(zhì)不同；較鈍形狀的冰塊表現(xiàn)近乎剛體，而較尖銳形狀的冰塊較易破碎，為能量耗散的主體；正方體冰塊的碰撞力和總內(nèi)能最大，因此可以作為最危險的局部形狀進行進一步研究；

3） 船側(cè)結(jié)構(gòu)發(fā)生較明顯的塑性變形，如航行在冰區(qū)，需采取一定程度的冰區(qū)加強措施。

［1］ Timco G.W.， Weeks W.F.. A review of the engineering properties of sea ice［J］. Cold Regions Science and Technology， 2010. 60 （2）： 107-129.

［2］ 丁德文. 工程海冰學(xué)概論［M］. 北京：海洋出版社，1999.

［3］ Schulson E.M.. Brittle failure of ice［J］. Engineering Fracture Mechanics， 2001. 68 （17）： 1839-1887.

［4］ Sodhi D.S.. Crushing failure during ice-structure interaction［J］. Engineering Fracture Mechanics， 2001. 68 （17）： 1889-1921.

［5］ Schwarz， J.， Weeks W.F.. Engineering properties of sea ice［J］. Journal of Glaciology， 1977. 19 （81）.

［6］ Gagnon R.E.. A numerical model of ice crushing using a foam analogue［J］. Cold Regions Science and Technology， 2011. 65（3）：335-350.

［7］ Jordaan I.J.. Mechanics of ice-structure interaction［J］. Engineering Fracture Mechanics， 2001. 68（17-18）： 1923-1960.

［8］ Jebaraj C.， et al.. Finite element analysis of ship/ice interaction［J］. Computers & Structures， 1992. 43（2）： 205-221.

［9］ 季順迎，岳前進. 工程海冰數(shù)值模型及應(yīng)用［M］. 北京：科學(xué)出版社， 2011.

［10］ Ralston T.D.. Yield and plastic deformation in ice crushing failure［C］//ICSI/AIDJEX Symposium on Sea Ice--Processes and Models， Seattle， Washington. 1977.

［11］ Jones S.J.. The confined compressive strength of polycrystalline ice［J］. Journal of Glaciology， 1982. 28： 171-177.

［12］ Rist M.A.， Murrell S. Ice triaxial deforntation and fracture［J］. Journal of Glaciology， 1994. 40（135）.

［13］ Derradji-Aouat A.. A unified failure envelope for isotropic fresh water ice and iceberg ice［C］//Proceedings of ETCE/OMAE Joint Conference Energy for the New Millenium， 2000.

［14］ Liu Z.， Amdahl J.， L?set S. Plasticity based material modelling of ice and its application to ship-iceberg impacts［J］. Cold Regions Science and Technology， 2011. 65（3）： 326-334.

［15］ Simof J.C.， Hughes T. Computational inelasticity［M］. 2008.

［16］ Masterson D.M.， et al. A revised ice pressure-area curve. in Proceedings［C］//Nineteenth International Conf. on Port and Ocean Engineering under Arctic Conditions， Dalian， China. 2007.

［17］ International S.O.. Petroleum and natural gas industries — arctic offshore structures［Z］， in ISO TC 67/SC 7/WG 8. 2010：Geneva， Switzerland. 434.

［18］ Storheim M.， et al. Iceberg shape sensitivity in ship impact assessment in view of existing material models［C］//ASME 2012 31stInternational Conference on Ocean， Offshore and Arctic Engineering. Rio de Janeiro， Brazil. 2012.

A New Constitutive Model of Ice Material for Ship-ice Interaction based on Ideal Elasto-plastic Property

HU Zhi-qiang1，2，GAO Yan1，YAO Qi1
（1. State Key Laboratory of Ocean Engineering， Shanghai Jiao Tong University， Shanghai 200240；2. Collaborative Innovation Center for Advanced Ship and Deep-Sea Exploration，Shanghai Jiao Tong University， Shanghai 200240）

Ship-ice interaction is currently one of the hot topics of the ship collision research， and the modeling of ice material is one of the key issues in ship-ice collision research. A constitutive model of ice material for ship-ice interaction based on the ideal elastic-plastic property is proposed in this paper， which calculate the element stress in plastic stage with the semi-implicit return mapping algorithm and describes the ice mechanical behavior with the Tsai-Wu yield criterion and empirical failure formula. Secondary development function is used to incorporate the ice material model into LS_DYNA program to validate the accuracy and applicability of the model. In accordance with the scenarios of collisions between a ship and ice blocks of different shapes， the study compares and analyzes the collision force， energy dissipation and etc. to discuss the influence of the ice shape on collision scenarios. The result shows that the pressure of the ice material model agrees well with existing standards in case the contact surface is large. The interaction processes of the ship-ice collision are different with different ice shapes. The blunt-shaped ice performs like a rigid body， while sharp-shaped ice is easy to be broken.

ship-ice collision； elastic-plastic material； constitutive model； influence of iceberg shape； numerical simulation

U661.7

A

2095-4069 （2016） 01-0065-09

10.14056/j.cnki.naoe.2016.01.013

2015-08-31

國家自然科學(xué)基金重點項目（51239007）。

胡志強，男，副教授，1975年生。2008年上海交通大學(xué)船舶及海洋結(jié)構(gòu)物設(shè)計制造專業(yè)獲博士學(xué)位，現(xiàn)于上海交通大學(xué)船舶海洋與建筑工程學(xué)院工作。