王澤平， 胡志強， 劉 昆， 陳 剛

(1.上海交通大學 海洋工程國家重點實驗室，上海 200240；2.哈爾濱工程大學 船舶工程學院，哈爾濱 150001；3.英國紐卡斯爾大學 工程學院，英國紐卡斯爾 NE1 7RU；4.江蘇科技大學 船舶與海洋工程學院，鎮(zhèn)江 212003；5.中國船舶及海洋工程設計研究院，上海 200011)

隨著海上航行船舶數量的不斷增加，船舶發(fā)生碰撞事故的概率也逐漸增加。盡管很多先進的避撞設備已經應用到船舶航行中，但船舶發(fā)生碰撞事故的風險仍然存在；船舶碰撞事故會導致環(huán)境污染、經濟損失甚至人員傷亡。有統計數據表明，船舶發(fā)生斜撞事故的概率要高于船舶正撞事故[1]。然而目前對于船舶碰撞事故的研究主要集中在船舶正撞場景，對于船舶斜撞場景開展的研究相對較少，因此需要針對船舶斜撞場景開展相關結構耐撞性研究。

船舶碰撞的研究方法可以分為四種，包括經驗公式法，試驗法，有限元數值仿真法和解析計算法。Minorsky[2]作為經驗公式法的開創(chuàng)者，通過對已發(fā)生的船舶碰撞事故結構損傷進行統計分析，進而提出內能吸收與結構損傷體積的經驗公式。試驗法分析船舶碰撞問題的準確性較好，不過研究成本較高。Pedersen等[3]通過試驗法研究船舶結構的耐撞性；Amdahl等[5]通過開展模型試驗研究船舶碰撞與擱淺時的結構變形機理；Liu等[6-7]通過落錘沖擊試驗研究了應變率的影響。隨著計算機技術的發(fā)展，有限元數值仿真法被廣泛應用于船舶碰撞問題的研究中。王自力等[8]通過數值仿真法研究了船舶碰撞動力學問題。劉敬喜等[9]通過數值仿真法研究了雙殼油船舷側結構的耐撞性。劉昆等[10]采用全耦合技術分析了船體結構的碰撞性能。此外，解析計算法是一種成本低且可以提供理想精度的方法，其通過提出結構的理論變形模型，運用塑性力學等理論，推導得到解析計算公式。Alexander[11]創(chuàng)造性地將解析計算法應用到薄板結構在軸向加載情況下的分析計算中。高振國等[12]通過解析計算法分析了強肋框結構在面內加載時的變形機理。Hong等[13]研究了船舶碰撞與擱淺過程中桁材結構的變形機理，并提出了解析公式。

近年來，隨著船舶噸位的增加，船舶球鼻艏的安裝也更加普遍，例如補給船，穿梭油輪等。船舶球鼻艏結構通常比較堅硬，發(fā)生船舶碰撞事故時更容易引起被撞船舷側結構的損傷，進而導致艙室進水、漏油等事故的發(fā)生。為了減少因為球鼻艏撞擊導致的損失，研究舷側結構在球鼻艏撞擊場景下的結構響應具有重要意義。然而大多數已有的研究都集中在球鼻艏正撞舷側結構的場景，比如，Lee等[14-15]研究了舷側外板在球鼻艏正撞場景下的結構響應。對于球鼻艏斜撞場景下舷側結構耐撞性的研究較少。在斜撞場景下，舷側結構的變形情況與撞擊角度相關性很大，舷側結構的主要構件，例如舷側外板、橫向肋板等構件在斜撞場景下的變形模式與正撞場景是不同的，因此需要針對舷側結構在球鼻艏斜撞場景下的變形機理開展研究。

1 解析計算基本理論

解析計算方法可以快速評估被撞船的耐撞性，通過應用塑性力學理論，求解舷側結構主要構件在球鼻艏斜撞場景下的抵抗力。根據上限定理，外力對結構的做功功率等于結構內部構件的能量耗散率[16]

(1)

(2)

(3)

(4)

N0=σ0t

(5)

式中：σ0為流動應力，t為舷側外板的厚度。

2 舷側外板的變形機理研究

2.1 舷側外板幾何變形模型

橢球形撞頭斜撞一塊邊長分別為2a和2b的矩形板，如圖1所示。該橢球形撞頭可以用如下方程的形式表示

(6)

圖1 橢球形撞頭斜撞矩形板示意圖Fig.1 A rectangular plate struck by an ellipsoid indenter

假設截面為橢圓的撞頭撞擊矩形板，如圖2所示。該橢圓截面是橢球形撞頭在xoy平面內的截面，那么橢圓形截面可以用如下方程的形式表示

(7)

圖2 截面為橢球形的撞頭斜撞矩形板示意圖Fig.2 A rectangular plate struck by an indenter (b1+b2=2b)

矩形板的變形模式如圖3所示。在矩形板受壓變形的過程中，左側部分矩形板的應變可以用如下方程表示

圖3 矩形板的變形模式示意圖Fig.3 Deformation model of the rectangular plate

因此，變形矩形板左側部分的應變率是

(8)

式中:l1是矩形板左側部分未發(fā)生彎曲部位的瞬時長度;α1是矩形板左側部分的瞬時旋轉角度。

(9)

同樣，右側部分矩形板的應變可以表示為

(10)

變形矩形板右側部分的應變率是

(11)

將式(9)和(11)代入式(3)中，可以得到膜拉伸能量耗散率為

(12)

N0_sideplate=σ0tp

(13)

式中：N0_sideplate是塑性膜拉伸力；tp是矩形板的厚度。

橢球形撞頭的撞擊深度可以用下式表示

(14)

其中，b1,b2的長度可以分別表示為

b1=x1tanα1+l1cosα1

(15)

b2=x2tanα2+l2cosα2

(16)

式中，(x1,y1)和(x2,y2)是矩形板未彎曲部分與彎曲部分的交點。x1,x2分別為變形矩形板未彎曲部分與彎曲部分的交點在x軸方向上的坐標值，如圖3所示。x1和x2可以表示為

(17)

(18)

由此可以得到橢球形撞頭的撞擊深度

(19)

橢球形撞頭的撞擊深度對瞬時旋轉角度求導，可以得到撞擊速率為

(20)

由于矩形板的彎曲能量耗散占總能量耗散的比例很小，所以通常忽略彎曲能量耗散。因此，矩形板在橢球形撞頭撞擊時的瞬時阻力可以表示為

(21)

同樣，假設截面為橢圓的撞頭撞擊矩形板，該橢圓截面是橢球形撞頭在xoz平面內的截面，那么橢圓形截面可以用如下方程的形式表示

(22)

矩形板受到撞擊時的瞬時阻力計算方法與上文提出的方法相同，可以表示為如下形式

(23)

Haris等的研究發(fā)現橢球形撞頭撞擊矩形板的瞬時阻力與上文提出的矩形板受到兩個方向上的橢圓形截面撞頭撞擊時的瞬時阻力關系可以表示為

(24)

在斜撞過程中，橢球形撞頭和矩形板之間的摩擦力可以表示為

FS=μ·Fp·sinβ

(25)

由于摩擦產生的能量耗散可以表示為

(26)

2.2 舷側外板破裂的預報

在碰撞過程中，如果外板發(fā)生破裂，那么外板的抵抗力會出現明顯下降，因此能夠較為準確地預測出外板的破裂十分有必要。通常臨界斷裂應變被當做板受到面外載荷作用時發(fā)生膜拉伸的最大應變，即認為當板發(fā)生膜拉伸的應變值達到臨界應變值時，板出現破裂。這一方法被廣泛使用，因此本文也采用這種方法預測外板的破裂。低碳鋼的塑性拉伸應變值通常在0.2～0.35，在船舶的初步設計階段，設計師可以根據具體情況確定這一數值的大小。將式(8)代入式(19)中，可以得到舷側外板最大膜拉伸應變和球鼻艏撞深之間的關系

(27)

式中，εm為外板的臨界斷裂應變值。根據式(27)可以得到舷側外板發(fā)生破裂時球鼻艏的撞擊深度。

2.3 舷側外板破裂后的解析計算公式

隨著碰撞過程的進行，球鼻艏撞深增加，舷側外板將會發(fā)生破裂。Wang等分析了板的變形過程，提出了預測板破裂后的阻力計算公式

(28)

式中：n為板破裂后的裂紋條數，n>2；l為裂縫的長度；φ為球鼻艏頂角的一半；μ為摩擦因數。

3 桁材變形模型

為了研究桁材結構在不同斜撞角度下的變形機理，首先對桁材的斜撞過程進行了數值模擬。選取三種典型的斜撞角度，數值模擬了球鼻艏結構在30°、45°和60°三種場景下斜撞桁材的過程。桁材在30°、45°和60°場景下的有限元模型、變形模型以及橫截面的變形過程如圖4～6所示。

(a) 有限元模型

(b) 變形模型

(c) 橫截面變形過程圖4 桁材在數值模擬中的變形(β=30°)Fig.4 Deformation of the web girder in numerical simulation(β=30°)

(a) 有限元模型

(b) 變形模型

(c) 橫截面變形過程圖5 桁材在數值模擬中的變形(β=45°)Fig.5 Deformation of the web girder in numerical simulation(β=45°)

(a) 有限元模型

(b) 變形模型

(c) 橫截面變形過程圖6 桁材在數值模擬中的變形 (β=60°)Fig.6 Deformation of the web girder in numerical simulation

從圖4～6中可以看出，在不同的斜撞角度下，桁材腹板的折疊高度比是相似的，因此根據數值模擬的結果提出了桁材的理論變形模型，如圖7和圖8所示。其中，b1為受壓桁材結構的一側長度，另一側長度為b2=b-b1。

圖7 桁材撞擊后的理論變形模型Fig.7 The deformation model of the web girder after collision

圖8 橫截面受壓折疊變形過程(黑點代表塑性鉸所在位置)Fig.8 Folding deformation process of the cross section (the black dots represent the location of the plastic hinges)

根據所提出的桁材理論變形模型，通過推導可以得到適用于不同斜撞角度的桁材變形阻力計算公式。限于篇幅，詳細推導過程請參閱Wang等[17]9-10。推導后得到的桁材結構在形成第1個褶皺時的瞬時阻力公式為

(29)

在形成第1個褶皺時的平均阻力公式為

(30)

根據上限定理，當平均阻力最小時，能量耗散也最小，即：

(31)

將式(30)代入式(31)中，可以得到：

(32)

同樣，當桁材結構形成第n個褶皺時(n≥2)，瞬時阻力公式為

(33)

式中：δ′為球鼻艏在斜撞方向上的位移距離；t為腹板的厚度。

桁材上帶板的阻力公式可以采用梁理論進行計算

(34)

式中：a是帶板寬度的一半；b1+b2是帶板的長度；tf是帶板的厚度。

4 橫向肋板變形模型

球鼻艏斜撞橫向肋板時，假設橫向肋板受到長度

為d的線荷載作用，切向碰撞力使橫向肋板上邊緣產生水平位移，垂向碰撞力使橫向肋板產生二維折疊變形，所提出的橫向肋板理論變形模型，分別如圖9和圖10所示。

圖9 橫向肋板碰撞前后的橫截面Fig.9 The cross section of transverse frame before and after collision

圖10 橫向肋板受壓變形模式Fig.10 Deformation mode of the crushing transverse frame

根據所提出的橫向肋板理論變形模型，通過推導可以得到適用于不同斜撞角度的切向力和垂向力的計算公式。限于篇幅，詳細推導過程請參閱Wang等[17]10-12。推導后得到的切向力和垂向力的計算公式如下

(35)

(36)

式中：δ′為球鼻艏在斜撞方向上的位移距離；t為肋板的厚度。

5 數值仿真

本文通過開展數值仿真分析來驗證所提出解析公式的合理性與準確性。

5.1 有限元模型

在數值仿真中，選取一艘雙殼油輪作為被撞船，一艘60 000 t船舶的橢球形球鼻艏作為撞擊船艏。撞擊船艏和被撞船舷側結構的有限元模型如圖11所示，表1列出了被撞船的主尺度。

圖11 油輪舷側結構和撞擊船艏的有限元模型Fig.11 Finite element model of the side structure of the tanker and the bow

表1 被撞船主尺度表Tab.1 Dimentions of the struck ship

使用有限元軟件LS_DYNA進行數值仿真，使用Johnson-Cook GTN材料模型[18]。撞擊船艏的材料屬性定義為剛體，被撞船舷側結構的材料參數如表2和表3所示。撞擊速度定義為3 m/s，摩擦因數定義為0.3[19]。有限元模型中采用的是四節(jié)點四邊形的Belytschko-Tsay單元。為了在保證數值模擬精度的同時有合理的計算時間，有限元模型的網格尺寸設置為200 mm。被撞船舷側兩端采用固定邊界條件；本文研究的是船舶斜撞過程中相對危險的情況，即撞擊船在碰撞過程中不發(fā)生偏轉，被撞船舷側結構在較大撞深情況下的損傷變形機理。在數值模擬中，被撞船和撞擊船之間的接觸設置為面面接觸。

表2 GTN模型參數Tab.2 GTN parameters

表3 Johnson-Cook模型參數Tab.3 Johnson-Cook parameters

5.2 碰撞場景定義

船舶碰撞時斜撞角度的定義如圖12所示，圖13為三種典型的碰撞位置，兩種碰撞場景的定義分別如表4和圖14所示。工況1中的碰撞位置為位置1，工況2中的碰撞位置為位置2，隨著碰撞的進行，撞擊船艏會經過位置3。

圖12 斜撞角度β的定義Fig.12 Definition of the collision angle β

圖13 三種典型碰撞位置Fig.13 Three typical collision positions

表4 不同碰撞場景中的碰撞位置Tab.4 Definition of collision positions in different collision scenarios

(a) 工況1

(b) 工況2圖14 兩種斜撞場景的定義Fig.14 Definition of two oblique collision cases

5.3 碰撞過程中參與構件的確定

在5.2節(jié)中確定了碰撞場景，可以結合撞擊球鼻艏結構的輪廓圖和被撞船舷側的結構圖確定參與斜撞過程中的船體構件。圖15(a)是工況1，撞擊過程通過將撞擊球鼻艏結構的輪廓圖和被撞船舷側的結構圖結合起來進行模擬。在碰撞開始時，與球鼻艏接觸的構件有舷側外板和桁材，所以舷側結構的抵抗力由這兩個構件產生，舷側外板的變形區(qū)域為CDHI；隨著碰撞的進行，橫向肋板也參與到舷側結構抵抗力的計算中，舷側外板的變形區(qū)域為CEGI，如圖15(b)所示。圖15(c)和(d)是工況2，在碰撞開始時，與球鼻艏接觸的構件有舷側外板和橫向肋板，舷側結構的抵抗力由這兩個構件產生，舷側外板的變形區(qū)域為AHJK，如圖15(c)所示；隨著碰撞的進行，舷側外板的變形區(qū)域為FGJK，如圖15(d)所示。

(a) 工況1

(b) 工況1

(c) 工況2

(d) 工況2圖15 碰撞過程中變形區(qū)域及參與構件的確定Fig.15 Determination of deformation area and participating components during collision

6 數值驗證與討論

使用有限元軟件LS_DYNA對兩種斜撞場景進行數值仿真，將得到碰撞過程中球鼻艏撞擊深度與碰撞力和能量耗散之間的關系曲線；解析方法中的能量耗散是通過對膜拉伸能量耗散率積分得到，將數值仿真結果與解析計算結果進行對比，如圖16～圖19所示。

圖16 碰撞場景1中碰撞力-撞深曲線對比結果Fig.16 The comparison results of resistance-depth curves in case 1

圖17 碰撞場景1中能量-撞深曲線對比結果Fig.17 The comparison results of energy-depth curves in case 1

圖18 碰撞場景2中碰撞力-撞深曲線對比結果Fig.18 The comparison results of resistance-depth curves in case 2

圖19 碰撞場景2中能量-撞深曲線對比結果Fig.19 The comparison results of energy-depth curves in case 2

從圖16～圖19的對比結果中可以看出，數值仿真方法和解析方法得出的結果吻合良好。隨著碰撞過程的進行，碰撞力和能量耗散曲線呈上升趨勢。當碰撞場景1和2中球鼻艏的撞擊深度達到3.2 m時，舷側外板發(fā)生破裂；舷側外板發(fā)生破裂后抵抗力會下降，能量耗散曲線的增長速度也出現下降，解析方法和數值仿真方法都能捕捉到這一特點。在舷側外板破裂初期，解析方法和數值仿真方法得到的碰撞力-撞深曲線存在一些偏差，這是因為在舷側外板破裂初期，膜拉伸能量耗散仍然是舷側外板最主要的能量耗散形式，舷側外板的變形模式不會在很短的時間內轉化為撕裂變形模式，這是解析方法和數值仿真方法得到的結果之間存在偏差的主要原因?？傮w來說，解析方法和數值仿真方法得到的結果對比良好，這也驗證了本文所提出的解析公式的合理性。

7 結 論

本文研究了船舶球鼻艏結構斜撞船舶舷側時，舷側結構主要構件的變形機理，提出了舷側外板、桁材和橫向肋板三種構件的理論變形模型，運用塑性力學理論推導得到了三種構件在不同斜撞角度下變形過程中的解析計算公式。定義了兩種典型的斜撞場景，通過有限元軟件LS_DYNA進行了數值仿真，數值仿真結果驗證了解析計算結果的準確性，解析計算公式可以很好地捕捉到舷側外板發(fā)生破裂后碰撞力和能量耗散增長速度下降這一特點。

球鼻艏也是船舶常見的設計形式，船舶斜撞事故發(fā)生的概率要高于船舶正撞事故，而目前針對船舶碰撞的研究主要集中在船舶正撞場景，對于船舶斜撞場景的研究是不充分的，因此，研究舷側結構在船舶球鼻艏結構斜撞場景下的變形機理是很有意義的。本文對舷側結構在球鼻艏斜撞場景下的變形機理進行了研究，所得到的解析公式對于船舶舷側結構耐撞性的評估具有較好的工程應用價值。