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

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

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

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

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

0 引 言

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

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

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

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

1.1 Rahman 法

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

1.2 CSR 法

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

1.3 算例分析

運用ANSYS 非線性有限元（FEM）、Rahman法、JBP 規(guī)范計算得到了4 種加筋板單元平均應(yīng)力—應(yīng)變曲線，并與Patrick［3］采用的ISUM 程序計算結(jié)果進行了對比。加筋板的尺寸（單位為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īng)驗公式擬合得到的關(guān)系曲線，初始撓度已經(jīng)包含在解析式里面，所以沒有這一項。非線性有限元計算時采用的初始撓度是按ISO 標(biāo)準(zhǔn)選取，以模態(tài)的形式疊加在模型上面［4］。

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

圖1 4 種尺寸加筋板單元采用不同方法計算得到的應(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

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

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

2.1 材料屬性

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

2.2 計算模型邊界條件及模型范圍

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

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

2.3 初始變形

由于焊接和施工等因素，加筋板必然會產(chǎn)生比較復(fù)雜的初始變形。本文采用文獻［4］中施加加筋板的初始撓度的方法，考慮了3 種類型的初始變形（圖3），分別為板上面的初始變形ωopl，筋上面的柱形初始變形ωoc及筋上面的側(cè)傾初始變形ωos，未考慮殘余應(yīng)力。為了得到初始變形形狀，首先對模型進行了特征值屈曲計算，得到了不同階的屈曲模態(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 有限元計算結(jié)果分析

為了建立內(nèi)河船舶加筋板單元平均應(yīng)力—應(yīng)變關(guān)系數(shù)據(jù)庫，首先必須確定影響關(guān)系曲線的主要因素。參照加筋板計算極限強度經(jīng)驗公式，選取板的柔度系數(shù)β 和筋的柔度系數(shù)（r 為慣性半徑）作為主要變量，在對內(nèi)河一些標(biāo)準(zhǔn)散貨船及油船的主尺度和結(jié)構(gòu)尺寸進行統(tǒng)計歸納之后，確定板的柔度系數(shù)β 需控制在1～3 之間，筋的柔度系數(shù)λ 需控制在0.1～1 之間。圖4 所示為幾組計算結(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ù)變大，相同應(yīng)變對應(yīng)的應(yīng)力會不斷減小，說明板變?nèi)踔螅咏畎宓呐R界應(yīng)力在不斷減小，抵抗能力變?nèi)?。由圖4（c）和圖4（d）可以看出，當(dāng)板的柔度系數(shù)一定時，隨著筋柔度系數(shù)的不斷增大，相同應(yīng)變對應(yīng)的應(yīng)力不斷減小，說明筋變?nèi)踔螅咏畎宓呐R界應(yīng)力在不斷減小，抵抗能力變?nèi)?，只是變化幅度沒有前面的大，因此，板的強弱對加筋板單元極限強度的影響更為重要。

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

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

3.1 程序迭代方法

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

3.2 插值方法

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

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

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

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

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

當(dāng)β0和λ0都不在節(jié)點上面，β0的左右節(jié)點分別為β1和β2，λ0的左右節(jié)點分別為λ1和λ2時，下一步需找到β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)河集裝箱船進行船體梁極限強度計算，采用MARC非線性有限元計算船體梁極限強度，并與文獻［5］中用直接法計算的結(jié)果進行比較。內(nèi)河集裝箱船橫剖面如圖5所示。

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

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

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

1）采用文獻［5］中的直接法進行計算。

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

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

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

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

2）采用非線性有限元法進行計算。

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

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

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

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

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

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

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

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

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

5 結(jié) 論

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

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

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

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

［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)河雙殼船極限強度分析［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.