一種考慮鋼材硬化各向同性的失效準則在船舶碰撞數值仿真中的應用

趙百惠， 胡志強， 陳 剛

(1. 上海交通大學 海洋工程國家重點實驗室， 上海 200240； 2. 英國紐卡斯爾大學 海洋科學與技術學院， 紐卡斯爾， 英國 NE1 7RU； 3. 高新船舶與深海開發裝備協同創新中心， 上海交通大學， 上海 200240；4. 中國船舶及海洋工程設計研究院， 上海 200011)

盡管全球各國均努力防止船舶事故，船舶泄漏污染仍時有發生。發生泄漏污染的主要原因之一是船舶在碰撞和擱淺事故中船體結構斷裂破損。船舶事故不僅會引起漏油和永久性的船舶結構損壞，而且可能導致人員損失和船舶交通堵塞，并對海洋環境造成惡劣影響[1]。除了規范全球船舶運輸行為，增強船上人員安全意識之外，在設計階段預報分析船體結構抗撞能力，并對船體結構進行優化，也是有效減少石油泄漏事故發生的重要措施。

船舶碰撞與擱淺研究方法一般分為三類，即試驗法，簡化解析法和有限元數值模擬法。其中，試驗一般可以提供較為精確和可靠的預測，然而因其需要大量人力和時間，而且價格昂貴，實船試驗很少進行。自20世紀90年代后期[2]，計算機技術的迅速發展使得大尺寸復雜結構船體模型的有限元數值仿真變得可行；相對而言，解析方法的發展較為緩慢。為了滿足造船工業對可靠性和降低成本的效益日益增長的需求，有限元仿真技術在船舶防撞性的直接定量預報以及擱淺場景的模擬應用更加頻繁。此外，數值仿真技術還經常用于對簡化解析方法的驗證應用[3]。

在船舶碰撞與擱淺的有限元數值仿真過程中，判斷模擬鋼板材料的有限元單元的斷裂失效，是一個十分重要的問題。例如，在球鼻艏撞擊船側結構時，船體外板的膜拉伸，是最主要的結構吸能形式。而一旦外板破裂，其吸收能量的能力會大幅度下降。外板破裂早了，低估了外板的吸能能力；外板破裂遲了，則高估了外板的吸能能力。因此，準確預報板材破裂，合理分析船舶的防撞性，就顯得十分關鍵。但是，不同的學者，當利用有限元方法進行數值仿真時，即使是研究同一個對象，也會得出不同的計算結果。造成這一情況的原因主要包括失效準則不同，有限元模型網格尺寸不同，使用的計算程序不同等原因[4]。其中，板材單元失效準則以及網格尺寸是兩個主要原因，并且二者相互影響。

當前國際學術界對船舶碰撞擱淺有限元應用中的鋼材斷裂失效準則已進行了大量研究。Hill[5]基于局部頸縮分析提出了判斷塑性開始的簡化方法，為BWH失效準則的基礎。T?rnqvist[6]依據系列拉伸試驗、壓痕試驗和三點彎曲試驗數據，改進并驗證了RTCL準則[7]。Stoughton等[8]提出了基于應力狀態的成形極限圖，避免了金屬對于應變路徑的依賴問題。在船舶碰撞擱淺的仿真分析中，目前采用較為簡化的臨界等效塑性應變法，而上述鋼材斷裂失效準則因其復雜性尚未得到廣泛應用。本文主要研究目的是基于國際學術界現有的失效準則，推導出一個改進的鋼板板材斷裂失效準則。通過驗證其合理性，對推廣更具準確性的新型失效準則具有重要意義。

1 有限元應用中的鋼材斷裂失效準則

目前，在船舶碰撞擱淺的有限元研究領域，國際學術界已經提出了多種鋼板材料失效準則。它們各有優缺點和不同的模擬適用場景，且均在進一步完善發展中。

1.1 臨界等效塑性應變法

對于一種材料，對應一個經驗塑性應變臨界值εcr，當厚度方向應變達到此臨界值，則認為斷裂發生。因此，臨界等效塑性應變法是在碰撞模擬中最簡便常用的方法。但此失效準則基于較為簡化，沒有考慮應力狀態的影響，結果并不十分準確[9]。

1.2 基于應變狀態的成形極限圖(FLD/FLC)

基于應變狀態的成形極限圖最開始由Keeler等[10]提出。此方法主要是通過材料在單向或雙向應力作用下，給出其在不同的應變路徑中得到的局部失穩極限應變ε1和ε2，以此為橫縱坐標繪制出條帶形的區域或曲線，得到此種材料的成形極限圖。在判斷材料失效時，將塑性失效初期的主應變(ε1,ε2)在成形極限圖中描繪出來，如果落在成形極限圖的成形極限曲線上部，板材變形時就會產生破裂，反之則是安全的。基于應變的成形極限圖的一般形式，如圖1所示。

圖1 基于應變狀態的成形極限圖一般形式

1.3 基于應力狀態的成形極限圖FLD/FLC

在發現基于應變的成形極限曲線的局限性后，Arrieux等[12]提出對于基于應變的成型極限曲線的改進。之后，Stoughton等[13]通過轉換應力應變狀態提出了基于應力狀態的成形極限圖。Stoughton通過對同一金屬樣本不同預應變情況下進行拉伸試驗，發現得到的成型極限曲線差別較大，證實了基于應變狀態得到的成形極限圖對于應變路徑的依賴。隨后，將曲線轉化成應力狀態下成形極限圖，可以發現非線性的應變變化路徑對基于應力狀態的成形極限圖影響不大。

Jie等[14]對此方法進行了改進，并將其應用于有限元法，證明了此方法在數值模擬中輸入的參數簡單，與實驗結果的比對也較合理。

1.4 RTCL失效準則

(1)

當失效累積到一個臨界值時，認為斷裂開始。累積的失效值可以表示為

(2)

在有限元模擬中，一旦積分點積累的失效值達到臨界值Dcr，則單元失效，從有限元結構中移除，表示斷裂失效開始。

1.5 BWH失效準則

金屬在斷裂之前通常是過度的塑性變形，在板和殼結構中，這經常被看作結構在受到外力后在窄帶中的過度塑性流動，其特征在于局部頸縮，材料在這個階段雖然沒有斷裂，但已達到臨界斷裂狀態。從微觀上講，這是因為在材料中已經出現顯著的弱點，其隨后將迅速結合并導致金屬的斷裂。

BWH準則將局部頸縮的開始視為失效，這在對大尺寸板殼結構的碰撞研究中是十分可取的[15]。一方面，結構在發生局部頸縮之后，其應力應變狀態均會發生明顯改變[16]，在有限元數值模擬中，對于較粗的網格劃分無法準確捕捉到此時的狀態變化，而較細的網格劃分會大大增加計算時間和資源，使用BWH準則可以有效減少網格尺寸的影響。此外，失穩的開始可以基于材料的應力應變曲線給出，簡化了此失效準則的檢驗過程[17]。

2 基于power law材料模型的BWH失效準則理論推導

本文提出一種基于power law材料模型的BWH失效準則，該準則適用于船舶碰撞與擱淺的數值仿真分析。

2.1 材料模型

材料模型的選擇是鋼材斷裂失效準則的基礎。考慮各向同性的塑性硬化材料，材料應力-應變曲線由Ludwik’ s power law公式表示，即

σ=Kεn

(3)

式中：K和n為材料的特定常數。

2.2 失效模型

Hill于1952年初步提出以局部頸縮開始判斷鋼材失效的失效準則。Hill假定頸縮方向與最大主應力成φ角度，一旦頸縮形成，沿頸縮帶方向的應變增量為0，則有

(4)

頸縮形成時，應變硬化效果與厚度方向上的尺寸減少并互相平衡，得到

(5)

一般認為主應力增量比例與主應力之間比例一致，即

(6)

(7)

式中：dλ為一非負的比例系數。

在線彈性階段，應力應變滿足Hook定律，存在比例對應的關系；而在塑性變形階段，應力與應變之間不再存在比例對應關系，應變不僅和應力狀態有關，還和變形歷史有關，此時則需要研究應力與應變增量之間的關系，這種增量形式表示的塑形本構關系稱為增量理論。由塑性力學增量理論，塑性應變增量可表示為

(8)

屈服函數f為勢函數，考慮屈服函數變化量的兩種表示方法

(9)

又由塑性流動理論，可得

(10)

(11)

若用最大主應力表示，則有

(12)

需注意的是，上文提出的失效準則只適用于β為負值的情況。而在船舶碰撞擱淺事故中，船體結構受外力導致的應變狀態是多樣的，如圖2所示。考慮平面應力狀態，會發生剪切、雙向拉伸等情況。

(a)

(b)

為使此失效準則適用于不同的應變狀態，現基于Bressan與Williams提出的三個假設，推導出適用于β不為負值時的失效模型。Bressan等[18]提出了預測金屬成型中局部頸縮的失效準則，它基于以下三個假定:① 剪切失穩發生在這樣一個截面，即材料單元沿此截面厚度t方向是沒有應變變化的，即dεt=0 ；② 失穩的發生是由于局部剪應力達到臨界值，即存在τcr，此為材料的一項基本性質；③ 忽略彈性應變。

圖3 應變莫爾圓

由假定①，考慮圖3中應變莫爾圓，圖3中，A點代表發生剪切失穩的材料單元，其厚度方向應變增量為

0，即

(13)

從而得到

(14)

另外由體積不可壓縮條件，即dε1+dε2+dε3=0，有：

dε3=-dε1(1+β)

(15)

可以得到

(16)

由假定②，考慮圖4平面應力狀態莫爾圓，找到對應τ=τcr的點A，可以得到

(17)

由此，可以得到

(18)

圖4 平面應力狀態莫爾圓

考慮β=0時Hill的公式，即由

(19)

得到

(20)

綜上，可以得出適用于power law 材料模型的BWH失效模型，如下

σ1=

(21)

3 材料失效準則的數值應用和驗證

LS-DYNA是著名的通用顯示非線性動力分析程序，可以模擬各種復雜的接觸非線性、材料非線性和幾何非線性問題。對于大部分非線性動力學問題，尤其是大變形結構的瞬時高度非線性問題，如船舶碰撞，對結構整體受力引起的整體應變計算過程當中一般均采用顯式求解方法。與隱式求解方法相比，顯式求解法不會涉及收斂性問題，也不用聯立方程組，避免了因聯立方程組而出現病態的無確定解的問題。本文利用LS-DYNA程序開展數值應用和驗證工作。

LS-DYNA提供自定義材料的二次開發功能，即編寫材料子程序并嵌入LS-DYNA主程序中。LS-DYNA用戶材料子程序從主程序得到應變值，且返回相關應力值，因此二次開發的關鍵問題為編寫材料的自定義本構關系。下面給出通過非線性有限元方法中應力更新算法中的半隱式向后Euler算法，依據上述失效準則作為判定，編寫此材料二次開發材料子程序的主要過程。

3.1 非線性有限元方法——圖形返回算法

在有限元計算時，單元的特性用單元剛度矩陣表示，計算剛度矩陣時，需要進行數值積分。在LS-DYNA中采用的數值積分方法為高斯積分法，因此積分點也稱高斯積分點，單元的應力值由返回給高斯積分點的應力值表示，而計算此應力值的過程是較為復雜的，常用的算法即為非線性有限元之中的應力更新算法。

積分率本構方程的數值算法稱為應力更新算法。應力更新算法可分為向前和向后Euler算法，前者公式簡單易于編程，但應力和內變量的更新值并不滿足屈服條件，常常導致不精確的結果；而向后Euler算法，它強化在時間步結束時屈服條件的一致性，即保證fn+1=0(f為屈服函數)，以避免離開屈服表面的漂移。有許多不同的積分本構方程算法，本文主要采用廣泛應用于實際中的圖形返回算法。其一般包括一個初始的彈性預測步，包含(在應力空間)對屈服表面的偏離，以及一個塑性調整步使應力返回到更新后的屈服表面fn+1=0。該方法可以用兩個組成部分來概括：一個積分算法，它將一組本構方程轉換為一組非線性代數方程；另一部分則是一個對非線性代數方程的求解算法[19]。

積分算法可以寫為

εn+1=εn+Δε，

qn+1=qn+Δλn+1hn,

fn+1=f(σn+1,qn+1)=0

(22)

可以注意到塑性應變增量

(23)

則第n+1步的應力可以表示成公式的形式

(24)

(25)

塑性參數增量為更新其余各變量的基礎。本文的材料塑性流動方向與屈服面法線成正比，則認為塑性流動是關聯的，即采用關聯流動法則，因此r=?fσ。從而得到塑性參數增量δλ

(26)

3.2 程序流程圖

上述應用于此失效準則的二次開發的圖形返回算法可由下文流程圖表示。從主程序得到初始值后，首先帶入屈服函數檢查屈服條件與零關系，若小于零，則退出迭代；若大于零，則計算塑性參數增量，更新各變量值，直至返回屈服面，退出迭代。之后根據此時應力值判斷單元是否失效，并返回應力值。

4 數值仿真應用與分析

研究中采用嵌入所提出材料失效準則的LS-DYNA二次開發程序，以一組NSWCCD船模擱淺試驗數據為依據[20]，驗證所提出的材料失效準則的適用性和合理性。該實驗對四種不同雙層底結構且縮尺比均為1：5的油船，模擬了在同一礁石下發生擱淺的場景，實驗中測出了擱淺過程中的船舶結構受力以及能量耗散情況。本文以其中的一種雙層底結構模型為對象，開展有限元建模以及數值計算，將有限元模擬結果與此模型的實驗結果進行比較，從而驗證斷裂失效準則的可靠性。用到的實驗基本數據如表1。

實驗裝置，如圖5所示。根據實驗要求，礁石頂端從內底板以下約0.05 m處開始進入；同時，船體模型應與水平方向成一定角度，滿足當礁石到達后艙壁時，礁石頂端與船模底部距離為兩倍雙層底間距[21]。

根據實驗模型尺寸，由PATRAN所建模型，如圖6所示，網格尺寸為30 mm。

表1 船模擱淺實驗基本數據

參數數值材料屈服強度σy/MPa283材料極限強度σu/MPa345材料硬化指數n0.22船長/m7.32船寬/m2.54船傾斜角度/(°)3.38礁石頂部半徑/m0.17礁石半頂角/(°)45相對速度14 kn(7.202 m/s)

之后將PATRAN中模型以key文件格式輸出，在LS-DYNA前處理程序中設置參數，并通過二次開發生成的求解器進行計算。因為本文為自定義的材料模型，在前處理設置參數時需要首先選擇相應自定義材料子程序序號；之后，按照子程序中各參數順序輸入對應材料參數，如材料屈服強度、泊松比、硬化參數等。有限元仿真計算時間約為24 h。擱淺場景數值仿真結構損傷截圖，如圖7所示。

圖5 實驗裝置圖(Rodd, 1996)

圖6 NSWCCD擱淺實驗船模有限元模型

圖7 擱淺場景數值仿真結果

圖8 結構變形阻力對比結果

圖9 結構變形能對比結果

結構變形阻力和結構變形能的數值仿真結果與實驗結果對比，如圖8、圖9所示。從兩圖中可知，本文推導出的失效模型數值分析得到的結果與傳統的臨界等效塑性應變方法模擬得到的結果相比，對于結構變形阻力以及結構變形能的模擬結果，均與實驗結果更為接近。因此，此種基于BWH失效準則的失效模型具有較高的可靠性，并具有比傳統失效應變方法精度更高的優點。

5 結 論

本文依據BWH鋼材斷裂失效準則的理論推導了基于power law材料模型的失效準則公式，并利用非線性有限元方法中的圖形返回算法，結合LS-DYNA材料二次開發理論，將此失效準則編程后嵌入了LS-DYNA程序中。利用更新后的LS-DYNA程序，以NSWCCD進行的一項實驗模型為對象，進行了有限元數值模擬，并與實驗結果和經驗等效塑性應變模擬結果對比，證實了此失效準則的可靠性以及更高的準確性。

此失效準則以局部頸縮開始作為判斷失效依據，避免了后頸縮區域的分析，因此對網格尺寸不敏感。也因此，此失效準則的使用尚有局限性和不確定性。此外，在前處理中輸入參數較多，流程較復雜，使用起來不如傳統的臨界塑性失效應變方法簡便。鋼材斷裂的失效準則仍需后續發展和改進，從而達到更廣泛的應用。