基于显式算法的船体梁极限强度非线性有限元分析

李政杰，赵 南，司海龙，虞 昊

(中国船舶科学研究中心，江苏无锡214082)

0 引言

随着现代船舶的大型化发展，船体结构强度问题不断凸显，有关船体结构真实强度的极限强度和破损后的极限承载能力，近几年来已受到越来越多的关注。多家船级社［1－3］的设计标准中已经规定了关于船体结构总纵极限强度的要求，未来向其他船型普及和贯彻极限强度要求的趋势亦愈来愈明显。

船体梁极限强度作为评价结构极限承载能力的重要指标，目前，其研究手段主要有直接计算法、逐步破坏法、理想结构单元法、实船事故调查和模型试验、非线性有限元法6种。表1分别对各种计算方法进行了对比，通过对比可以看出非线性有限元法是计算和分析船体极限强度的一种强有力工具，能够合理地处理结构的边界条件、载荷形式、模型细节、单元类型等影响因素，可以综合考虑结构的初始缺陷和材料的非线性特性，比较直观地模拟结构屈曲崩溃的过程并清晰地反映其后屈曲特性。

船体梁极限强度非线性有限元计算中，隐式算法较为常见，但其计算过程需要多次迭代，存在收敛问题，增量步长较大，总的增量步数较少。隐式算法对增量步的大小没有限制，增量步大小主要取决于求解精度和收敛条件，每个增量步上需要进行多次迭代，求解复杂非线性问题的计算代价较大。因此，有必要采用显式算法，改变求解策略，提高计算效率。

1 显示算法

显式算法和隐式算法都是求解非线性问题强有力的工具。所谓显式算法和隐式算法，是指求解策略的不同，即数学上的出发点不一样。显式算法中一般适用于求解应力波效果显著或模拟时间非常短暂 (一般小于1 s)的分析，如高度不连续的高速的动力学分析、复杂的接触分析、复杂的后屈曲分析、高度非线性分析的准静态分析、边界条件极度不连续的问题、材料磨损和材料失效等，其计算过程不需要迭代，不存在收敛问题，增量步很小，总的增量步很多，增量步的大小与模型的最大固有频率有关，与载荷类型和载荷持续时间无关。

显式算法一般采用中心差分法显示地对运动方程在时间域上进行积分，利用上一步增量步的平衡方程动态地计算下一步增量步的状态。显示求解器的中心差分算子在每一增量步的开始时刻t满足动态平衡方程，通过计算t时刻的加速度，将速度的计算结果推进到t+，将位移的计算结果推到t+Δt。每一增量步结束时刻的状态完全由增量步开始时刻的位移、速度和加速度决定。

隐式算法是将一个分析步分解为多个增量步，载荷则以增量的形式逐步施加，在每个载荷增量步中都要进行一系列迭代，在每次迭代中使用上一次的迭代得到的修正刚度矩阵进行计算。这样，对计算结果逐步修正，直到满足平衡方程，达到收敛。在每一增量步内都需要对平衡方程组 (隐式表达式)进行迭代求解，每次都需要对非线性刚度矩阵求逆，计算成本相对较大，但是隐式算法的计算比较稳定，大多数时候是无条件的稳定，对步长的要求不高。

在计算复杂的船体梁极限强度时发现由于结构复杂，单元数量过多，或者材料的非线性、结构几何非线性、边界非线性等的影响运用隐式算法常常不收敛，得不到可靠地数值解。因此，可考虑使用动态分析近似模拟静力加载过程，即准静态分析，采用显式算法进行求解。

假设结构的非线性静态平衡方程为:

式中:P为载荷列阵;I为内力列阵。

显示有限元法运用中心差分法对在时间上对运动方程进行积分，应用一个增量步的动力学条件计算下一个增量步的动力学条件。在增量步开始时，程序求解动力学平衡方程，表示为节点质量矩阵M和加速度u·的乘积等于节点的内力，即施加的外力P与内力I的差值:

在当前增量步开始时，假设是t时刻，计算加速度为:

因为显式算法总是采用对角的或者集中的质量矩阵，所以求解加速度并不复杂，不需要同时求解两联立方程。任何节点的加速度完全取决于节点质量和作用在节点的合力，使得节点计算成本非常低。

采用中心差分法对加速度在时间上进行积分，在计算速度的变化时假设加速度为常数。应用这个速度的变化值加上前一个增量步中点的速度确定当前增量步中点的速度:

速度对时间的积分加上增量步开始时的位移来确定在增量步结束时的位移:

这样，在增量开始时提供了满足动力学平衡条件的加速度。得到了加速度，在时间上“显式地”前推速度和位移。所谓“显式”是指在增量步结束时的状态仅依靠于该增量步开始前的位移、速度和加速度。为了得到精确的结果，需要假定在载荷随时间增长非常慢，这样在增量步中的加速度几乎为零，这就是所谓的准静态。下面给出显示动力方法的总结:

1)节点计算

①动力平衡方程

②对时间进行显示积分

2)单元计算

②根据本构关系计算应力σ

③ 集成节点内力I(t+Δt)，

3)设置时间t为t+Δt，返回到步骤1)。

2 显示算法与隐式算法算例对比

利用大型通用软件Abaqus，分别采用显式算法和隐式算法对Nishihara箱型梁模型进行对比计算，并将计算结果与试验值［4］比较，验证和比较隐式算法和显示准静态算法在计算船体梁极限强度中的实用性和可靠性。

Nishihara箱型梁的结构尺寸如表2，采用理想弹塑性材料，模型两端采用多点约束，两主节点简支，施加垂向角位移。假设焊接残余拉压应力影响相互抵消，忽略残余应力影响;初始变形采用一阶屈曲模态来模拟，幅值为1.48 mm。箱型梁的有限元模型如图1所示。

图1 Nishihara箱型梁有限元模型Fig.1 The FEM model of Nishihara box girder

表2 箱型梁模型参数Tab.2 The parameter of box girder

图2和图3分别给出了隐式算法和显式算法极限状态下箱型梁应力云图，对比两图可看出，2种应力分布云图基本一致。

图2 极限状态应力分布云图 (隐式算法)Fig.2 The stress tensor of ultimate state(implicit algorithm)

图3 极限状态应力分布云图 (显式算法)Fig.3 The stress tensor of ultimate state(explicit algorithm)

图4给出了2种算法的载荷－位移曲线，2条弯矩－转角曲线基本一致，在达到峰值点之前，结构处于线性状态，两曲线重合;达到峰值后，进入塑性的部件逐渐增加，显式算法存在着能量损耗，与静态分析存在细微差别，由此表现出来结构的后屈曲性能存在稍微差别。总体而言，两者计算的结果比较一致。

图4 不同算法弯矩－转角曲线Fig.4 Bending moment-angle curve of different algorithms

表3给出了2种算法的计算结果及试验值。隐式算法与试验结果比较的相对误差为3.06%，显式算法与试验结果比较的相对误差为5.44%;显式算法相对隐式算法计算结果稍大，这是由于显式算法属于准静态求解，其求解步长稍大，以致结果存在稍许误差;而显式算法求解的时间约为隐式求解时间的一半，其计算效率较高。当结构较为复杂时，隐式算法的计算时间将大大增加。

表3 箱型梁计算结果对比 (108 N·mm)Tab.3 The calculate result of box girder comparison(108N·mm)

对比分析2种计算的结果:隐式算法的计算精度较高，但其计算效率较低;显式算法计算精度相对隐式算法稍低，但也满足工程设计需要，其计算效率较高。当结构较为复杂时，隐式算法的收敛性较差，往往不能得到收敛的结果，造成计算失败;显式算法其收敛性较好，选择适当的计算步长，能够得到满足要求的收敛解。综合考虑计算的精度和效率，基于显式算法的船体梁极限强度计算具有较高的优势。

3 目标实船计算

Dow's 1/3 test hull模型［5］具有较高的代表性，许多学者对其进行了对比计算，其结果常被用来检验计算方法的精度和适用性。Dow试验模型尺寸如图5所示。

图5 Dow's test hull剖面Fig.5 The section of Dow's test hull

图6 有限元模型Fig.6 The FEM model

运用显式算法，采用单跨模型计算其极限承载能力。模型初始变形形状采用其一阶模态，其变形幅值为试验测得变形量;忽略残余应力的影响，其有限元模型如图6所示。采用理想弹塑性材料，模型两端采用多点约束，两主节点简支，施加垂向转角位移。

图7和图8分别给出了Dow模型中拱和中垂2种极限状态下的应力分布云图，图9给出了模型的端面弯矩－转角曲线。

图7 中拱极限状态下应力分布云图Fig.7 The stress tensor of hog ultimate state

图8 中垂极限状态下应力分布云图Fig.8 The stress tensor of sag ultimate state

图9 Dow模型 (1/2)弯矩－转角曲线Fig.9 Bending moment-angle curve of Dow model(1/2)

表4给出了2种极限状态下模型的极限强度值，并与其他学者的计算结果进行比较。显示计算的结果稍微偏大，这是因为采用1/2对称模型，约束强化结构，并且未考虑残余应力的影响。

表4 Dow's test hull计算结果比较Tab.4 The comparison for the result of Dow's test hull

4 结语

本文对有限元计算中的显示算法进行了介绍，并采用显式算法和隐式算法对Nishihara模型的极限强度的计算结果进行对比分析，发现显式算法具有较好的精度且计算效率较高;基于显式算法对Dow船体梁模型的极限强度进行计算，并与相关研究结果比较，结果表明该方法具有良好的计算精度，可以用于实船的极限强度分析，为船体结构的设计提供参考。

［1］ IACS．Common Structural Rules for Bulk Carriers［S］．2008．

［2］ ABS Rules for Building and Classing Steel Vessel Rules［S］．Houston，2013．

［3］ DNV Rules for Classification of Ships［S］．Norway，2013．

［4］ Nishihara S．Ultimate longitudinal strength of midship cross section［J］．Naval Arch．Ocean Eng，1984，22:200 －214．

［5］ ISSC report of special task committee VI．2，ultimate hull girder strength．In:Ohtsubo H， Sumi Y， editors．Proceedings of the 14th international ship and offshore structures congress，Nagasaki，Japan;October 2000．