内河船舶极限强度计算的逐步破坏法程序设计

骆文刚，杨 平，崔虎威，白小溪

1 中船重工建筑工程设计研究院有限责任公司，北京100121

2 武汉理工大学交通学院，湖北武汉430063

3 中国船级社大连分社，辽宁大连116013

0 引 言

在船舶设计与强度评估中，如果不能明确结构的极限承载能力，依旧单纯采用基于线弹性理论的许用工作应力方法，便无法确定结构的真实安全极限，也不再适合当前的设计理念和发展要求。ISO，IMO，IACS 等国际组织已于2006年开始正式将船体梁极限强度评估内容列入规范，将之用于指导大型船舶结构的设计。目前，极限强度在内河船舶设计中还未采用。但随着内河船舶标准化、大型化发展趋势的日趋明显，小吨位船舶因其运能及经济性劣势将被快速淘汰，上万吨的船舶将越来越常见。由此，大型船舶的安全性就须引起足够重视，有必要在极限强度方面开展相关研究。

就现阶段的发展来说，计算船体梁极限载荷的方法主要有非线性有限元法、理想结构单元法（ISUM）、逐步破坏法和直接法。其中，逐步破坏法由于其计算效率高，结果也比较可靠。该方法的核心为加筋板单元平均应力—应变关系，本文将主要在这方面开展一些研究工作，并将其运用于实际的计算过程之中。

1 加筋板单元平均应力—应变关系

采用简化计算方法——逐步破坏法来计算船体梁的总纵极限弯矩，最核心的方面就是要合理处理加筋板单元的载荷—端缩曲线，即单元的平均应力—应变关系。

1.1 Rahman 法

在Rahman 法［1］中，加筋板单元的压缩行为分3 个阶段，即稳定阶段、非卸载阶段和卸载阶段。采用梁柱理论，将加筋板单元理想化为梁柱模型，并基于边缘屈服准则，假设加筋板带板上边缘或者面板下边缘这其中之一发生了屈服，那么整个结构就发生失效。Rahman 法由于其公式相对简单，计算结果也比较可靠，应用广泛。但其不足之处在于，其将加筋板单元简单地等价为了梁柱模型，而且只能计算加筋板单元梁柱屈曲模式，同时在处理初始挠度及残余应力等方面也不是很合理。

1.2 CSR 法

在CSR 计算方法［2］中，将单元分为了硬角单元、普通加强筋单元及横骨架板格单元。对于硬角单元、受拉的普通加强筋单元和横骨架板格单元，将发生弹塑性崩溃破坏。对于受压单元，根据其不同的失效模式，采用了不同的计算公式，但这样将各种失效模式人为地独立分开，与实际情况也有一定差别。而且CSR 法中的应力—应变关系表达式主要是在一些经验公式的基础上总结得到，对于按CSR 规范设计的海船比较实用，如果用到内河船舶，则不一定实用。例如，计算光板极限强度的Faulkner 经验公式在计算比较厚的板时较为准确，但在计算薄板时误差较大。同时，内河船舶结构尺寸与海船也相差较大，也需要一种更为合理的计算方法。

1.3 算例分析

运用ANSYS 非线性有限元（FEM）、Rahman法、JBP 规范计算得到了4 种加筋板单元平均应力—应变曲线，并与Patrick［3］采用的ISUM 程序计算结果进行了对比。加筋板的尺寸（单位为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 规范计算方法因是基于经验公式拟合得到的关系曲线，初始挠度已经包含在解析式里面，所以没有这一项。非线性有限元计算时采用的初始挠度是按ISO 标准选取，以模态的形式叠加在模型上面［4］。

4 种加筋板计算模型采用不同的计算方法得到的关系曲线对比如图1 所示。其中ISUM 代表理想结构单元法计算结果。

图1 4 种尺寸加筋板单元采用不同方法计算得到的应力—应变关系曲线对比图Fig.1 The stress-strain curves of four different dimensions of the stiffened plate element calculated using different methods

从图中曲线的极值点可以得到加筋板在单向受压下的极限强度，而且可以看出，Rahman 法很明显有一段平行的阶段，即非卸载阶段，是一种理想化的模型；图中ISUM 法是在没有考虑塑性变形情况下计算的结果，所以曲线变化相对比较缓慢。在图1（c）和图1（d）的算例中，有限元法在经历峰值之后，应力突然下降很快，这是由于部分单元进入塑性阶段，发生了塑性应变，因而结构抵抗能力也很快降低。从极限值上可以看出，有限元法与其它方法也比较接近。总体来说，采用非线性有限元法对加筋板单元进行大变形分析，而且也考虑了材料的非线性及结构由于焊接产生的初始变形，计算结果比较可靠。

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

2.1 材料属性

本节计算的加筋板单元，假设其材料是连续的，各向同性且弹性理想塑性（elastic-perfectly plastic），屈服准则为Von-Mises 屈服准则。材料的屈服极限σy=235 MPa，杨氏模量E=205 800 MPa，泊松比v=0.3。

2.2 计算模型边界条件及模型范围

计算模型取1 个强横梁间距，为了充分考虑板格发生的屈曲失效模式，将加筋板单元横向各扩展b/2 范围，延伸至相邻加强筋的位置。采用文献［4］中提出的边界条件，即在A-D，B-C 边施加Ry，Rz转动约束，且在加筋板的纵向边界处施加z向约束，两边中点施加x 向约束；在A-B，C-D 边施加Rx，Rz转动约束，且带板上的节点施加z 向约束，两边中点施加y 向约束；由于计算模型取自连续加筋板中的一部分，因而各边应施加直边界条件，如图中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 初始变形

由于焊接和施工等因素，加筋板必然会产生比较复杂的初始变形。本文采用文献［4］中施加加筋板的初始挠度的方法，考虑了3 种类型的初始变形（图3），分别为板上面的初始变形ωopl，筋上面的柱形初始变形ωoc及筋上面的侧倾初始变形ωos，未考虑残余应力。为了得到初始变形形状，首先对模型进行了特征值屈曲计算，得到了不同阶的屈曲模态。然后根据板和筋不同的屈曲模式选择相应的模态，分别放大相应的倍数至一定的幅值，幅值大小按ISO 18072 取值，为Ao=0.1β2t，Bo=Co=0.001 5a，a为跨距，β=(b/t)为板的柔度系数。

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

2.4 有限元计算结果分析

为了建立内河船舶加筋板单元平均应力—应变关系数据库，首先必须确定影响关系曲线的主要因素。参照加筋板计算极限强度经验公式，选取板的柔度系数β 和筋的柔度系数（r 为惯性半径）作为主要变量，在对内河一些标准散货船及油船的主尺度和结构尺寸进行统计归纳之后，确定板的柔度系数β 需控制在1～3 之间，筋的柔度系数λ 需控制在0.1～1 之间。图4 所示为几组计算结果。

图4 内河船舶加筋板单元平均应力—应变关系Fig.4 Unit average stress-strain curves of inland ship stiffened plate

由图4（a）和图4（b）可以看出，随着加筋板应变的不断增加，应力也随之不断变大；当经过峰值之后，应力开始缓慢减小。当筋的柔度系数一定时，板的柔度系数变大，相同应变对应的应力会不断减小，说明板变弱之后，加筋板的临界应力在不断减小，抵抗能力变弱。由图4（c）和图4（d）可以看出，当板的柔度系数一定时，随着筋柔度系数的不断增大，相同应变对应的应力不断减小，说明筋变弱之后，加筋板的临界应力在不断减小，抵抗能力变弱，只是变化幅度没有前面的大，因此，板的强弱对加筋板单元极限强度的影响更为重要。

另外，曲线对应的峰值就是加筋板单元在这种单向受压情况下的极限强度值，随着板柔度系数和筋柔度系数的不断变大，极限强度值会变小。本文为了得到不同尺寸加筋板单元的平均应力—应变关系，将板柔度系数控制在了1.0～3.0 之间，以0.2 的大小递增，筋的柔度系数控制在了0.1～1.0 之间，以0.1 的大小递增，一一对应，并建立每条曲线的数据库。如果某个加筋板单元对应的柔度系数超出了数据库范围，将通过外插的方法实现。

3 逐步破坏法程序设计

3.1 程序迭代方法

逐步破坏法主要包含2 个迭代计算。第1 个为中和轴迭代，即在某一曲率下，对瞬时中和轴位置进行调整直至满足剖面拉压平衡，然后计算所有单元的弯矩贡献，获得该指定曲率下的弯矩值。第2 个为曲率迭代，即依次增加曲率，计算得到各曲率下的上述迭代过程弯矩值，得到弯矩—曲率（M-χ）曲线，并取曲线峰值作为计算模型的极限垂向弯矩。本文采用FORTRAN 语言编写了逐步破坏法计算程序。

3.2 插值方法

当单元类型为普通加筋板单元时，程序会自行判定其受拉或者受压，如受拉，则选择弹塑性破坏曲线；如受压，会根据输入的加筋板单元尺寸大小计算得到带板柔度系数β 和筋柔度系数λ，然后在数据库里自动选择合适的曲线。

当 β=m×0.2 且 5 ≤m ≤15，λ=n×0.1 且1 ≤n ≤10（其中m，n 为正整数）时，数据里会存在与之对应的曲线。

当 β=m×0.2 且 5 ≤m ≤15，λ ≠n×0.1 且1 ≤n ≤10（其中m，n 为正整数）时，则按如下插值方法：

假设需要计算应变为εave时的应力σave，已知该加筋板单元对应的柔度系数分别为β0和λ0，β0在节点上面，而λ0则在节点之间，且左右节点分别为λ1和λ2，下一步需找到β0，λ1对应的曲线F1(x)及β0，λ2对应的曲线F2(x)，则

另外，如果λ0在节点上面，而β0在节点之间，插值方法类似。

当β0和λ0都不在节点上面，β0的左右节点分别为β1和β2，λ0的左右节点分别为λ1和λ2时，下一步需找到β1和λ1对应的曲线F1(x)，β1和λ2对应的曲线F2(x)，β2和λ1对应的曲线F3(x)以及β2和λ2对应的曲线F4(x)，则

最后，

4 算例及结果分析

应用本文的程序对某一内河集装箱船进行船体梁极限强度计算，采用MARC非线性有限元计算船体梁极限强度，并与文献［5］中用直接法计算的结果进行比较。内河集装箱船横剖面如图5所示。

图5 某内河集装箱船横剖面图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，设计吃水d =4.3 m，肋距0.5 m。

1）采用文献［5］中的直接法进行计算。

首先，将船体剖面离散成7 块加筋板，即甲板、两外舷侧、两内舷侧、外底和内底。然后，根据Paik 和Lee 经过试验，并考虑构件初始变形得到受压加筋板极限强度的经验公式

式中：β 和λ 分别为板和筋的柔度系数；σy为材料屈服应力。

运用上式计算每个加筋板的极限屈曲强度，再计算各构件的相当面积，然后，再假设船体崩溃时中拱和中垂两种工况下的剖面应力分布情况（图6）。根据平衡条件（截面压应力等于拉应力），计算得到中和轴位置，将整个剖面对中和轴取矩并求和，即可得到此状况下的计算弯矩值。

图6 直接法假设剖面应力分布状况Fig.6 Section stress distribution hypothesis by direct method

2）采用非线性有限元法进行计算。

假设强框架足够强，不会发生整个板架失稳的情况，计算模型取两强框架之间，板和骨材全部采用板单元模拟。纵骨腹板上划6 个网格，面板上划2 个网格，纵骨间距划10 个网格，单元总数为6 万多个，具体模型如图7 所示。在左、右两端面施加刚体约束，主节点设置在型心位置。左端面主节点施加沿x，y，z 轴的线位移约束Ux，Uy，Uz及绕x，z 轴的转动约束Rx，Rz；右端面主节点施加沿y，z 轴的线位移约束Uy，Uz及绕x，z 轴的转动约束Rx，Rz。然后，选择合适的模态放大一定的比例模拟初始变形。计算模型不考虑残余应力的影响。图8所示为有限元法计算得到的弯矩—曲率曲线。

图7 69.8 m 内河集装箱船有限元计算模型Fig.7 FE calculation model of 69.8 m inland river container ship

图8 有限元计算得到的弯矩—曲率曲线Fig.8 Moment-curvature curve by FEM

3）采用本文编写的逐步破坏法程序进行计算。

将船体离散为加筋板单元，其中有82 个一般加筋板单元，28 个硬角单元。在输入文件中，必须写入每一个单元带板和筋的尺寸及其垂向坐标，以及材料的屈服应力等已知参数。最后，计算得到弯矩—曲率关系图（图9）。图中，Mp为塑性弯矩值（横剖面所有构件都达到屈服应力时对应的弯矩值），用本文程序计算的值为292 146 kN·m ；Mo为瞬时曲率对应的弯矩值。

图9 本文程序计算得到的弯矩—曲率曲线Fig.9 Moment-curvature curve by the developed program in this paper

3 种方法的计算结果对比如表1 所示。

表1 某内河集装箱船极限强度计算结果Tab.1 The limit strength calculation results of an inland river container ship

由表1 中可以看出，3 种方法计算得到的极限弯矩值还是比较接近的。与直接法计算结果相比，误差大小比文献［5］中的结果要相对偏小，主要是因为在计算加筋板屈曲强度时采用的是Paik的经验公式［6-8］，Paik 在文献中也提到该计算公式比较保守，偏于安全。

5 结 论

1）算例结果表明，本文设计的逐步破坏法计算程序与非线性有限元法（MARC）和直接法的计算结果相近，具有良好的计算精度。

2）与常规的海船极限强度计算方法相比，内河船在船舶结构单元尺寸、船用钢材的性能及初始变形的大小方面都有很大差别，而且，很多基于大型海船结构强度计算的经验公式在这里也已不太适用。与传统的海船采用一些理论及经验公式简化计算方法相比，本文采用非有限元法直接计算得到单元的应力—应变曲线，更加符合实际情况。

3）本文的加筋板单元应力—应变关系数据库由有限元的系列计算建立，因此，对于加筋板的失效模式，其比Rahman 的梁—柱模型的计算方法更加符合实际，避免了将各种失效模式人为独立分开所可能产生的问题。

4）本文建立的加筋板单元应力—应变关系数据库已包含了大多数内河船舶结构尺寸，对于内河不同类型的船舶极限强度计算既方便，又实用。

［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］邓林峰，杨平.内河双壳船极限强度分析［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.