內(nèi)河船舶極限強(qiáng)度計(jì)算的逐步破壞法程序設(shè)計(jì)

駱文剛，楊 平，崔虎威，白小溪

1 中船重工建筑工程設(shè)計(jì)研究院有限責(zé)任公司，北京100121

2 武漢理工大學(xué)交通學(xué)院，湖北武漢430063

3 中國船級社大連分社，遼寧大連116013

0 引 言

在船舶設(shè)計(jì)與強(qiáng)度評估中，如果不能明確結(jié)構(gòu)的極限承載能力，依舊單純采用基于線彈性理論的許用工作應(yīng)力方法，便無法確定結(jié)構(gòu)的真實(shí)安全極限，也不再適合當(dāng)前的設(shè)計(jì)理念和發(fā)展要求。ISO，IMO，IACS 等國際組織已于2006年開始正式將船體梁極限強(qiáng)度評估內(nèi)容列入規(guī)范，將之用于指導(dǎo)大型船舶結(jié)構(gòu)的設(shè)計(jì)。目前，極限強(qiáng)度在內(nèi)河船舶設(shè)計(jì)中還未采用。但隨著內(nèi)河船舶標(biāo)準(zhǔn)化、大型化發(fā)展趨勢的日趨明顯，小噸位船舶因其運(yùn)能及經(jīng)濟(jì)性劣勢將被快速淘汰，上萬噸的船舶將越來越常見。由此，大型船舶的安全性就須引起足夠重視，有必要在極限強(qiáng)度方面開展相關(guān)研究。

就現(xiàn)階段的發(fā)展來說，計(jì)算船體梁極限載荷的方法主要有非線性有限元法、理想結(jié)構(gòu)單元法（ISUM）、逐步破壞法和直接法。其中，逐步破壞法由于其計(jì)算效率高，結(jié)果也比較可靠。該方法的核心為加筋板單元平均應(yīng)力—應(yīng)變關(guān)系，本文將主要在這方面開展一些研究工作，并將其運(yùn)用于實(shí)際的計(jì)算過程之中。

1 加筋板單元平均應(yīng)力—應(yīng)變關(guān)系

采用簡化計(jì)算方法——逐步破壞法來計(jì)算船體梁的總縱極限彎矩，最核心的方面就是要合理處理加筋板單元的載荷—端縮曲線，即單元的平均應(yīng)力—應(yīng)變關(guān)系。

1.1 Rahman 法

在Rahman 法［1］中，加筋板單元的壓縮行為分3 個(gè)階段，即穩(wěn)定階段、非卸載階段和卸載階段。采用梁柱理論，將加筋板單元理想化為梁柱模型，并基于邊緣屈服準(zhǔn)則，假設(shè)加筋板帶板上邊緣或者面板下邊緣這其中之一發(fā)生了屈服，那么整個(gè)結(jié)構(gòu)就發(fā)生失效。Rahman 法由于其公式相對簡單，計(jì)算結(jié)果也比較可靠，應(yīng)用廣泛。但其不足之處在于，其將加筋板單元簡單地等價(jià)為了梁柱模型，而且只能計(jì)算加筋板單元梁柱屈曲模式，同時(shí)在處理初始撓度及殘余應(yīng)力等方面也不是很合理。

1.2 CSR 法

在CSR 計(jì)算方法［2］中，將單元分為了硬角單元、普通加強(qiáng)筋單元及橫骨架板格單元。對于硬角單元、受拉的普通加強(qiáng)筋單元和橫骨架板格單元，將發(fā)生彈塑性崩潰破壞。對于受壓單元，根據(jù)其不同的失效模式，采用了不同的計(jì)算公式，但這樣將各種失效模式人為地獨(dú)立分開，與實(shí)際情況也有一定差別。而且CSR 法中的應(yīng)力—應(yīng)變關(guān)系表達(dá)式主要是在一些經(jīng)驗(yàn)公式的基礎(chǔ)上總結(jié)得到，對于按CSR 規(guī)范設(shè)計(jì)的海船比較實(shí)用，如果用到內(nèi)河船舶，則不一定實(shí)用。例如，計(jì)算光板極限強(qiáng)度的Faulkner 經(jīng)驗(yàn)公式在計(jì)算比較厚的板時(shí)較為準(zhǔn)確，但在計(jì)算薄板時(shí)誤差較大。同時(shí)，內(nèi)河船舶結(jié)構(gòu)尺寸與海船也相差較大，也需要一種更為合理的計(jì)算方法。

1.3 算例分析

運(yùn)用ANSYS 非線性有限元（FEM）、Rahman法、JBP 規(guī)范計(jì)算得到了4 種加筋板單元平均應(yīng)力—應(yīng)變曲線，并與Patrick［3］采用的ISUM 程序計(jì)算結(jié)果進(jìn)行了對比。加筋板的尺寸（單位為mm）為：

I板：2 400×800×15 角鋼250×90×10/15

II板：2 400×800×15 扁鋼220×20

III板：4 000×800×15 角鋼250×90×10/15

IV 板：4 000×800×15 型材400×11/150×19

在初始撓度處理方面，Rahman 法取為Δ =a/750，a 為加筋板跨距長度；JBP 規(guī)范計(jì)算方法因是基于經(jīng)驗(yàn)公式擬合得到的關(guān)系曲線，初始撓度已經(jīng)包含在解析式里面，所以沒有這一項(xiàng)。非線性有限元計(jì)算時(shí)采用的初始撓度是按ISO 標(biāo)準(zhǔn)選取，以模態(tài)的形式疊加在模型上面［4］。

4 種加筋板計(jì)算模型采用不同的計(jì)算方法得到的關(guān)系曲線對比如圖1 所示。其中ISUM 代表理想結(jié)構(gòu)單元法計(jì)算結(jié)果。

圖1 4 種尺寸加筋板單元采用不同方法計(jì)算得到的應(yīng)力—應(yīng)變關(guān)系曲線對比圖Fig.1 The stress-strain curves of four different dimensions of the stiffened plate element calculated using different methods

從圖中曲線的極值點(diǎn)可以得到加筋板在單向受壓下的極限強(qiáng)度，而且可以看出，Rahman 法很明顯有一段平行的階段，即非卸載階段，是一種理想化的模型；圖中ISUM 法是在沒有考慮塑性變形情況下計(jì)算的結(jié)果，所以曲線變化相對比較緩慢。在圖1（c）和圖1（d）的算例中，有限元法在經(jīng)歷峰值之后，應(yīng)力突然下降很快，這是由于部分單元進(jìn)入塑性階段，發(fā)生了塑性應(yīng)變，因而結(jié)構(gòu)抵抗能力也很快降低。從極限值上可以看出，有限元法與其它方法也比較接近。總體來說，采用非線性有限元法對加筋板單元進(jìn)行大變形分析，而且也考慮了材料的非線性及結(jié)構(gòu)由于焊接產(chǎn)生的初始變形，計(jì)算結(jié)果比較可靠。

2 加筋板單元有限元分析方法

2.1 材料屬性

本節(jié)計(jì)算的加筋板單元，假設(shè)其材料是連續(xù)的，各向同性且彈性理想塑性（elastic-perfectly plastic），屈服準(zhǔn)則為Von-Mises 屈服準(zhǔn)則。材料的屈服極限σy=235 MPa，楊氏模量E=205 800 MPa，泊松比v=0.3。

2.2 計(jì)算模型邊界條件及模型范圍

計(jì)算模型取1 個(gè)強(qiáng)橫梁間距，為了充分考慮板格發(fā)生的屈曲失效模式，將加筋板單元橫向各擴(kuò)展b/2 范圍，延伸至相鄰加強(qiáng)筋的位置。采用文獻(xiàn)［4］中提出的邊界條件，即在A-D，B-C 邊施加Ry，Rz轉(zhuǎn)動(dòng)約束，且在加筋板的縱向邊界處施加z向約束，兩邊中點(diǎn)施加x 向約束；在A-B，C-D 邊施加Rx，Rz轉(zhuǎn)動(dòng)約束，且?guī)О迳系墓?jié)點(diǎn)施加z 向約束，兩邊中點(diǎn)施加y 向約束；由于計(jì)算模型取自連續(xù)加筋板中的一部分，因而各邊應(yīng)施加直邊界條件，如圖中A-B，C-D 施加x 方向直邊界條件，該截面所有點(diǎn)在受力時(shí)x 向位移相同；A-D，B-C施加y 方向直邊界條件，該截面所有點(diǎn)在受力時(shí)y向位移相同。邊界條件如圖5 所示（其中x 方向指筋的長度方向）。

圖2 加筋板單元有限元模型示意圖Fig.2 The FE model of stiffened plate element

2.3 初始變形

由于焊接和施工等因素，加筋板必然會(huì)產(chǎn)生比較復(fù)雜的初始變形。本文采用文獻(xiàn)［4］中施加加筋板的初始撓度的方法，考慮了3 種類型的初始變形（圖3），分別為板上面的初始變形ωopl，筋上面的柱形初始變形ωoc及筋上面的側(cè)傾初始變形ωos，未考慮殘余應(yīng)力。為了得到初始變形形狀，首先對模型進(jìn)行了特征值屈曲計(jì)算，得到了不同階的屈曲模態(tài)。然后根據(jù)板和筋不同的屈曲模式選擇相應(yīng)的模態(tài)，分別放大相應(yīng)的倍數(shù)至一定的幅值，幅值大小按ISO 18072 取值，為Ao=0.1β2t，Bo=Co=0.001 5a，a為跨距，β=(b/t)為板的柔度系數(shù)。

圖3 加筋板的初始撓度形式示意圖Fig.3 Stiffened plate with initial deflection sketch map

2.4 有限元計(jì)算結(jié)果分析

為了建立內(nèi)河船舶加筋板單元平均應(yīng)力—應(yīng)變關(guān)系數(shù)據(jù)庫，首先必須確定影響關(guān)系曲線的主要因素。參照加筋板計(jì)算極限強(qiáng)度經(jīng)驗(yàn)公式，選取板的柔度系數(shù)β 和筋的柔度系數(shù)（r 為慣性半徑）作為主要變量，在對內(nèi)河一些標(biāo)準(zhǔn)散貨船及油船的主尺度和結(jié)構(gòu)尺寸進(jìn)行統(tǒng)計(jì)歸納之后，確定板的柔度系數(shù)β 需控制在1～3 之間，筋的柔度系數(shù)λ 需控制在0.1～1 之間。圖4 所示為幾組計(jì)算結(jié)果。

圖4 內(nèi)河船舶加筋板單元平均應(yīng)力—應(yīng)變關(guān)系Fig.4 Unit average stress-strain curves of inland ship stiffened plate

由圖4（a）和圖4（b）可以看出，隨著加筋板應(yīng)變的不斷增加，應(yīng)力也隨之不斷變大；當(dāng)經(jīng)過峰值之后，應(yīng)力開始緩慢減小。當(dāng)筋的柔度系數(shù)一定時(shí)，板的柔度系數(shù)變大，相同應(yīng)變對應(yīng)的應(yīng)力會(huì)不斷減小，說明板變?nèi)踔螅咏畎宓呐R界應(yīng)力在不斷減小，抵抗能力變?nèi)酢Ｓ蓤D4（c）和圖4（d）可以看出，當(dāng)板的柔度系數(shù)一定時(shí)，隨著筋柔度系數(shù)的不斷增大，相同應(yīng)變對應(yīng)的應(yīng)力不斷減小，說明筋變?nèi)踔螅咏畎宓呐R界應(yīng)力在不斷減小，抵抗能力變?nèi)酰皇亲兓葲]有前面的大，因此，板的強(qiáng)弱對加筋板單元極限強(qiáng)度的影響更為重要。

另外，曲線對應(yīng)的峰值就是加筋板單元在這種單向受壓情況下的極限強(qiáng)度值，隨著板柔度系數(shù)和筋柔度系數(shù)的不斷變大，極限強(qiáng)度值會(huì)變小。本文為了得到不同尺寸加筋板單元的平均應(yīng)力—應(yīng)變關(guān)系，將板柔度系數(shù)控制在了1.0～3.0 之間，以0.2 的大小遞增，筋的柔度系數(shù)控制在了0.1～1.0 之間，以0.1 的大小遞增，一一對應(yīng)，并建立每條曲線的數(shù)據(jù)庫。如果某個(gè)加筋板單元對應(yīng)的柔度系數(shù)超出了數(shù)據(jù)庫范圍，將通過外插的方法實(shí)現(xiàn)。

3 逐步破壞法程序設(shè)計(jì)

3.1 程序迭代方法

逐步破壞法主要包含2 個(gè)迭代計(jì)算。第1 個(gè)為中和軸迭代，即在某一曲率下，對瞬時(shí)中和軸位置進(jìn)行調(diào)整直至滿足剖面拉壓平衡，然后計(jì)算所有單元的彎矩貢獻(xiàn)，獲得該指定曲率下的彎矩值。第2 個(gè)為曲率迭代，即依次增加曲率，計(jì)算得到各曲率下的上述迭代過程彎矩值，得到彎矩—曲率（M-χ）曲線，并取曲線峰值作為計(jì)算模型的極限垂向彎矩。本文采用FORTRAN 語言編寫了逐步破壞法計(jì)算程序。

3.2 插值方法

當(dāng)單元類型為普通加筋板單元時(shí)，程序會(huì)自行判定其受拉或者受壓，如受拉，則選擇彈塑性破壞曲線；如受壓，會(huì)根據(jù)輸入的加筋板單元尺寸大小計(jì)算得到帶板柔度系數(shù)β 和筋柔度系數(shù)λ，然后在數(shù)據(jù)庫里自動(dòng)選擇合適的曲線。

當(dāng) β=m×0.2 且 5 ≤m ≤15，λ=n×0.1 且1 ≤n ≤10（其中m，n 為正整數(shù)）時(shí)，數(shù)據(jù)里會(huì)存在與之對應(yīng)的曲線。

當(dāng) β=m×0.2 且 5 ≤m ≤15，λ ≠n×0.1 且1 ≤n ≤10（其中m，n 為正整數(shù)）時(shí)，則按如下插值方法：

假設(shè)需要計(jì)算應(yīng)變?yōu)棣臿ve時(shí)的應(yīng)力σave，已知該加筋板單元對應(yīng)的柔度系數(shù)分別為β0和λ0，β0在節(jié)點(diǎn)上面，而λ0則在節(jié)點(diǎn)之間，且左右節(jié)點(diǎn)分別為λ1和λ2，下一步需找到β0，λ1對應(yīng)的曲線F1(x)及β0，λ2對應(yīng)的曲線F2(x)，則

另外，如果λ0在節(jié)點(diǎn)上面，而β0在節(jié)點(diǎn)之間，插值方法類似。

當(dāng)β0和λ0都不在節(jié)點(diǎn)上面，β0的左右節(jié)點(diǎn)分別為β1和β2，λ0的左右節(jié)點(diǎn)分別為λ1和λ2時(shí)，下一步需找到β1和λ1對應(yīng)的曲線F1(x)，β1和λ2對應(yīng)的曲線F2(x)，β2和λ1對應(yīng)的曲線F3(x)以及β2和λ2對應(yīng)的曲線F4(x)，則

最后，

4 算例及結(jié)果分析

應(yīng)用本文的程序?qū)δ骋粌?nèi)河集裝箱船進(jìn)行船體梁極限強(qiáng)度計(jì)算，采用MARC非線性有限元計(jì)算船體梁極限強(qiáng)度，并與文獻(xiàn)［5］中用直接法計(jì)算的結(jié)果進(jìn)行比較。內(nèi)河集裝箱船橫剖面如圖5所示。

圖5 某內(nèi)河集裝箱船橫剖面圖Fig.5 Cross section of an inland river container ship

船體強(qiáng)框架間距為2 000 mm，材料特性為：彈性模量E = 205 800 MPa，泊松比v = 0.3，屈服強(qiáng)度σy=235 MPa。

主尺度為：總長Loa= 69.8 m，垂線間長Lpp=66.6 m，型寬B = 16.0 m，型深D = 5.0 m，設(shè)計(jì)吃水d =4.3 m，肋距0.5 m。

1）采用文獻(xiàn)［5］中的直接法進(jìn)行計(jì)算。

首先，將船體剖面離散成7 塊加筋板，即甲板、兩外舷側(cè)、兩內(nèi)舷側(cè)、外底和內(nèi)底。然后，根據(jù)Paik 和Lee 經(jīng)過試驗(yàn)，并考慮構(gòu)件初始變形得到受壓加筋板極限強(qiáng)度的經(jīng)驗(yàn)公式

式中：β 和λ 分別為板和筋的柔度系數(shù)；σy為材料屈服應(yīng)力。

運(yùn)用上式計(jì)算每個(gè)加筋板的極限屈曲強(qiáng)度，再計(jì)算各構(gòu)件的相當(dāng)面積，然后，再假設(shè)船體崩潰時(shí)中拱和中垂兩種工況下的剖面應(yīng)力分布情況（圖6）。根據(jù)平衡條件（截面壓應(yīng)力等于拉應(yīng)力），計(jì)算得到中和軸位置，將整個(gè)剖面對中和軸取矩并求和，即可得到此狀況下的計(jì)算彎矩值。

圖6 直接法假設(shè)剖面應(yīng)力分布狀況Fig.6 Section stress distribution hypothesis by direct method

2）采用非線性有限元法進(jìn)行計(jì)算。

假設(shè)強(qiáng)框架足夠強(qiáng)，不會(huì)發(fā)生整個(gè)板架失穩(wěn)的情況，計(jì)算模型取兩強(qiáng)框架之間，板和骨材全部采用板單元模擬。縱骨腹板上劃6 個(gè)網(wǎng)格，面板上劃2 個(gè)網(wǎng)格，縱骨間距劃10 個(gè)網(wǎng)格，單元總數(shù)為6 萬多個(gè)，具體模型如圖7 所示。在左、右兩端面施加剛體約束，主節(jié)點(diǎn)設(shè)置在型心位置。左端面主節(jié)點(diǎn)施加沿x，y，z 軸的線位移約束Ux，Uy，Uz及繞x，z 軸的轉(zhuǎn)動(dòng)約束Rx，Rz；右端面主節(jié)點(diǎn)施加沿y，z 軸的線位移約束Uy，Uz及繞x，z 軸的轉(zhuǎn)動(dòng)約束Rx，Rz。然后，選擇合適的模態(tài)放大一定的比例模擬初始變形。計(jì)算模型不考慮殘余應(yīng)力的影響。圖8所示為有限元法計(jì)算得到的彎矩—曲率曲線。

圖7 69.8 m 內(nèi)河集裝箱船有限元計(jì)算模型Fig.7 FE calculation model of 69.8 m inland river container ship

圖8 有限元計(jì)算得到的彎矩—曲率曲線Fig.8 Moment-curvature curve by FEM

3）采用本文編寫的逐步破壞法程序進(jìn)行計(jì)算。

將船體離散為加筋板單元，其中有82 個(gè)一般加筋板單元，28 個(gè)硬角單元。在輸入文件中，必須寫入每一個(gè)單元帶板和筋的尺寸及其垂向坐標(biāo)，以及材料的屈服應(yīng)力等已知參數(shù)。最后，計(jì)算得到彎矩—曲率關(guān)系圖（圖9）。圖中，Mp為塑性彎矩值（橫剖面所有構(gòu)件都達(dá)到屈服應(yīng)力時(shí)對應(yīng)的彎矩值），用本文程序計(jì)算的值為292 146 kN·m ；Mo為瞬時(shí)曲率對應(yīng)的彎矩值。

圖9 本文程序計(jì)算得到的彎矩—曲率曲線Fig.9 Moment-curvature curve by the developed program in this paper

3 種方法的計(jì)算結(jié)果對比如表1 所示。

表1 某內(nèi)河集裝箱船極限強(qiáng)度計(jì)算結(jié)果Tab.1 The limit strength calculation results of an inland river container ship

由表1 中可以看出，3 種方法計(jì)算得到的極限彎矩值還是比較接近的。與直接法計(jì)算結(jié)果相比，誤差大小比文獻(xiàn)［5］中的結(jié)果要相對偏小，主要是因?yàn)樵谟?jì)算加筋板屈曲強(qiáng)度時(shí)采用的是Paik的經(jīng)驗(yàn)公式［6-8］，Paik 在文獻(xiàn)中也提到該計(jì)算公式比較保守，偏于安全。

5 結(jié) 論

1）算例結(jié)果表明，本文設(shè)計(jì)的逐步破壞法計(jì)算程序與非線性有限元法（MARC）和直接法的計(jì)算結(jié)果相近，具有良好的計(jì)算精度。

2）與常規(guī)的海船極限強(qiáng)度計(jì)算方法相比，內(nèi)河船在船舶結(jié)構(gòu)單元尺寸、船用鋼材的性能及初始變形的大小方面都有很大差別，而且，很多基于大型海船結(jié)構(gòu)強(qiáng)度計(jì)算的經(jīng)驗(yàn)公式在這里也已不太適用。與傳統(tǒng)的海船采用一些理論及經(jīng)驗(yàn)公式簡化計(jì)算方法相比，本文采用非有限元法直接計(jì)算得到單元的應(yīng)力—應(yīng)變曲線，更加符合實(shí)際情況。

3）本文的加筋板單元應(yīng)力—應(yīng)變關(guān)系數(shù)據(jù)庫由有限元的系列計(jì)算建立，因此，對于加筋板的失效模式，其比Rahman 的梁—柱模型的計(jì)算方法更加符合實(shí)際，避免了將各種失效模式人為獨(dú)立分開所可能產(chǎn)生的問題。

4）本文建立的加筋板單元應(yīng)力—應(yīng)變關(guān)系數(shù)據(jù)庫已包含了大多數(shù)內(nèi)河船舶結(jié)構(gòu)尺寸，對于內(nèi)河不同類型的船舶極限強(qiáng)度計(jì)算既方便，又實(shí)用。

［1］RAHMAN M K，CHOWDHURY M. Estimation of ulti?mate longitudinal bending moment of ships and box girders［J］. Journal of Ship Research，1996，40（3）：244-257.

［2］IACS. Common structural rules for bulk carriers［S］.2008.

［3］PATRICK K. Development of new idealized structural unit method for the collapse analysis of stiffened plate structures［D］.Japan：Hiroshima University，2001.

［4］PAIK J K，KIM B J，SEO J K.Methods for ultimate lim?it state assessment of ships and ship-shaped offshore structures：Part II stiffened panels［J］. Ocean Engi?neering，2008，35（2）：271-280.

［5］鄧林峰，楊平.內(nèi)河雙殼船極限強(qiáng)度分析［J］.船海工程，2008，37（6）：33-35.DENG Linfeng，YANG Ping. Ultimate strength analysis for inland double-hull ships［J］. Ship and Ocean Engi?neering，2008，37（6）：33-35.

［6］PAIK J K，SEO J K. Nonlinear Finite Element Method models for ultimate strength analysis of steel stiff?ened-plate structures under combined biaxial compres?sion and lateral pressure actions—Part. I Plate ele?ments［J］. Thin-Walled Structures，2009，47（8/9）：1008-1017.

［7］PAIK J K，KIM B J，SEO J K. Methods for ultimate limit state assessment of ships and ship shaped off?shore structures：Part III hull girders［J］. Ocean Engi?neering，2008，35（2）：281-286.

［8］PAIK J K，THAYAMBALLI A K，YANG S H. Residual strength assessment of ships after collision and ground?ing［J］.Marine Technology，1998，35（1）：38-54.