基于断裂力学的船型网箱疲劳损伤评估研究

刘安民，谢凌俊，苏一鸣，丁洪涛，何文涛

（中国海洋大学 工程学院，山东 青岛 266000）

利用简化计算方法进行海洋结构物疲劳寿命评估时，波浪载荷与结构响应都是通过经验公式计算的。这一方面会限制简化计算方法在缺乏数据样本的新型结构物上的应用，另一方面也会带来计算精度上的问题[1-3]。直接计算方法是公认较为准确的海洋结构物疲劳强度计算方法，与简化计算方法相比，直接计算方法中结构的疲劳载荷并不是通过经验公式计算的，而是通过波浪载荷软件，直接运用势流理论和莫里森方程对有限元模型进行计算得到的，并且通过有限元结构分析得到了结构的应力。这不仅能更好地反映全船疲劳载荷分布，也能够更好地反映出结构细节之处的载荷特点[4]。

谱分析方法是直接计算方法的一个重要方法，谱分析方法物理意义明确，过程合理，考虑因素全[5]，被认为是一种准确的方法，也得到了较为广泛的研究。赖明雁[6]认为，谱分析方法是建立在真实的海况与真实装载工况下的直接计算方法，计算精度高。徐帅[7]利用谱分析方法开展了LNG 船疲劳可靠性的分析，对谱分析方法的计算精度表示了肯定。宋宪仓[8]展开了半潜平台关键节点疲劳寿命评估的研究，并对结果进行了宽带修正，研究表明，相比于确定性方法，谱分析方法能够考虑波浪能量的分布，结果更为精确。曾琪[9]利用挪威船级社开发的海洋结构分析软件SESAM 对江海直达船进行了基于S-N曲线法的疲劳谱分析研究，并且将结果与简化计算方法进行了比较，认为谱分析相比于简化计算方法能得到更加准确的结果。白晓东[10]开展了深海网箱结构时域方法和频域方法的疲劳分析，并探讨了包括Rayleigh 分布、Gamma 分布和GEV 分布在内的多种分布模型对疲劳寿命结果的影响，表明对于非高斯过程，利用Rayleigh 分布对应力范围概率密度进行描述可能会产生较大的误差。由于谱分析方法是建立在应力响应窄带高斯过程的基础上的，因此不少学者对非高斯过程[11]、双峰或者三峰高斯过程[12]、高低频载荷共同作用[13]、窄带修正[14]等方面也做了不少工作，以提高谱分析方法的精度。

为了加深对基于断裂力学方法的海洋结构物疲劳裂纹扩展分析与寿命预测研究的认识和了解，形成更加成熟完善的研究体系，学者们从不同层面展开了广泛的研究。Tanaka 等[15-17]对管节点结构开展了疲劳裂纹扩展分析研究，丰富了海洋工程中常见节点结构的疲劳性能研究。王小松[18]以在役FPSO 为对象，根据英国标准BS 7910 进行了基于断裂力学的疲劳研究，并进行了多种规范之间的比较，研究结果对于营运安全具有一定的指导意义。

本文所采用的海洋结构物疲劳寿命评估方法将谱分析方法和断裂力学方法相结合，利用有限元计算的方式得到结构物在真实海况下的载荷条件，从疲劳载荷的角度保证疲劳评估的精度[19]，然后通过断裂力学进行疲劳评估。本文将以船型网箱为例，进行海洋结构物的运动特性以及结构应力分析，并利用子模型技术进行局部应力分析。之后，运用多尺度子模型实现壳单元到实体单元的转变，并在实体单元部分进行裂纹扩展，得到疲劳寿命。

1 疲劳分析方法和有限元模型

1.1 谱分析方法

谱分析方法的主要分析流程：

1）对海洋结构物进行水动力分析与应力响应分析，并提取疲劳热点的应力传递函数Hstress。

2）基于海洋结构物为线性系统的假设，根据海况i的波浪谱，求得应力响应谱密度函数。

3）基于应力响应是窄带高斯过程的假设，求得应力范围的概率密度函数f(ΔS)。

4）根据线性累积损伤理论，可以得到海况i下每年的疲劳损伤度，见式（1）。

式中：pi为海况概率；f0i为海况i对应的平均跨零率；ni为循环次数；Ni为循环破坏的寿命；S表示应力范围；A和m为疲劳实验参数。

5）循环步骤1）、2）、3），在对所有短期海况进行遍历后即可得到1 a 内的总损伤度，见式（2）。

1.2 裂纹扩展增量计算

对于裂纹稳定扩展阶段，最熟悉的扩展模型莫过于Paris 模型[20]：

式中：C和m为材料参数；为裂纹扩展速率；ΔK为应力强度因子范围；ΔKth为应力强度因子范围门槛值。当应力强度因子范围ΔK小于门槛值ΔKth时，裂纹将不能够扩展。

单一短期海况下的裂纹扩展增量E(da)可以根据式（4）–（6）进行计算：

式中：ΔNi为第i个短期海况下的循环次数；ΔTi为对应短期海况下的持续时间；fi为对应海况下的平均过零率；fiΔK为应力强度因子范围的概率密度函数；m0和m2为应力强度因子响应谱0 阶和2 阶谱矩。

由此可得，单一短期海况下的裂纹扩展增量为：

1.3 有限元模型

1.3.1 模型参数和坐标系

船型网箱主要由船艏舱室、船艉浮箱以及船艏船艉之间的箱型梁组成。船型网箱主要参数见表1。

表1 船型网箱主要参数Tab.1 Main parameters of ship shaped cages

坐标系原点在船底中线上，x轴正方向由船艏指向船艉，y轴由左舷指向右舷，z轴由船底指向上，如图1 所示。同时，波浪入射角θ的定义方式也在图1 中进行了展示，波浪入射角θ即波浪入射方向与x轴正方向之间的夹角。

图1 网箱参考坐标系定义Fig.1 Definition of the cage reference coordinate system

1.3.2 网箱有限元模型

由于海洋结构物主尺度较大，结构复杂，同时结构的厚度效应也会在水动力分析中有明显的体现，因此对于海洋结构物，一般利用壳单元以及梁单元进行有限元建模。网格主要为四边形单元，在结构不连续、曲率变化大处等将利用三角形单元进行离散。对于网箱上大量的T 形材、角钢等细长结构将利用梁单元进行建模，如图2a 所示。

图2 船型网箱有限元模型Fig.2 Finite element model of ship shaped cage:a) element type;b) constraint settings

结构应保持自由动态平衡状态，需要对其自由度进行约束[21]。DNV 等船级社做出了相应的规范要求[22]，采用三点约束法对网箱进行约束，如图2b 所示。约束点1 位于船艏底板中线纵舱壁下，约束其x、y、z三个自由度上的平动位移；约束点2 位于船艉右浮箱横纵型材的交点下，约束其y、z两个自由度上的平动位移；约束点3 位于船艉左浮箱横纵型材的交点下，约束其z轴方向上的自由度。

1.3.3 服役环境与载荷

海洋结构物在服役期间，同时遭受风浪流等载荷的共同作用，由于风和流的作用较弱，因此本文只考虑波浪载荷的作用。本文船型网箱服役于北大西洋海域，采用北大西洋海域波浪散布图。进行频域内分析时，由于结构的对称性，单位波幅规则波入射方向设为0°～180°，具体为0°、30°、60°、90°、120°、150°、180°。同时，波浪频率设置为 0.2～2.4 rad/s，间隔0.1 rad/s。

1.4 子模型

针对海洋结构物这类尺寸能够达到百米，并且结构复杂、细节丰富的超大型结构，在建模时保留全部细节将会导致建模工作量的显著增加，同时也会大大降低计算效率，还可能带来应力奇异现象等不利影响。因此，一般利用简化模型、粗网格模型对结构进行整体建模与分析，而采用子模型技术进行局部区域的建模与应力分析。

子模型技术的主要分析流程如图3 所示，主要步骤如下[23-24]：

图3 SESAM 子模型Fig.3 SESAM submodel

1）全局模型的建模与分析。对全局模型进行简化建模，采用较粗的网格进行单元划分，然后进行全局模型的分析。根据全局模型的分析结果，确定需要进行重点分析的局部区域（子模型区域）。

2）子模型的建模。根据实际尺寸、结构以及分析目标等，对子模型进行精细化建模，添加整体模型中未考虑的细节结构，并利用精细网格进行单元划分。

3）子模型边界继承。根据子模型的切割边界，从全局模型的分析结果中读取计算当前子模型的边界条件。

4）子模型分析。利用子模型进行局部结构的应力分析等，得到局部结构更精确的应力分析结果等。

在对大型结构进行有限元分析时，为了兼顾计算效率与计算精度，通常不会采用单一尺度（如壳单元、实体单元、梁单元）的单元进行结构离散，而是利用多尺度子模型技术建立多尺度耦合模型，在保证计算精度的同时提高计算速率，壳-实体耦合模型便是其中的重要形式[25]。在本文海洋结构物裂纹扩展分析中，壳-实体耦合模型的采用主要基于以下2 点考虑：1）载荷计算是在SESAM 中以壳、梁单元的形式完成的，因此壳单元的保留有利于更加灵活便捷地实现载荷的传递；2）实体单元的采用主要为了体现结构细节，得到更加精准的局部应力，实现裂纹在实体中的扩展。

1.5 多尺度耦合模型

壳-实体耦合模型中，一部分单元为壳单元，一部分单元为实体单元。由于壳单元节点具有六自由度，而实体单元节点仅仅具有3 个平动自由度，因此壳单元与实体单元连接处将会出现自由度不一致的情况，需要进行内部分布耦合约束以实现节点之间的协调[26]。

壳-实体耦合在 ABAQUS 中的实现（通过shell-to-solid coupling 命令）如图4 所示。为了减少后续有限元计算的时间，壳-实体耦合模型的范围只是SESAM 中子模型的局部，相当于从子模型中再次提取最为关心的疲劳热点部分作为子模型，这样处理的好处是既能更加突出重点区域的应力分布，也能缩小模型，提高计算效率。从SESAM 中子模型导出分析区域（二级子模型），并将其中实体区域内结构以实体单元重新建模，最后通过shell-to-solid coupling实现壳体单元与实体单元之间的连接，形成壳-实体耦合子模型。

图4 壳-实体耦合子模型Fig.4 Shell-solid coupling submodel

此外，壳-实体耦合子模型的载荷是通过节点对节点的形式从上一级子模型中继承的，并赋予到图4中壳单元最外层节点上。这部分工作主要是将从SESAM 软件中导出的边界上所有节点的六自由度运动逐一导入到ABAQUS 壳-实体耦合子模型中对应的节点上，从而实现载荷的传递，如图5 所示。

图5 载荷传递示意图Fig.5 Schematic diagram of load transfer

2 应力分析结果及疲劳分析

2.1 运动特性分析

在SESAM 软件中，根据三维势流理论，可以求得波浪场的速度势，并可以进一步求解出单位波幅规则波作用下结构的六自由度运动响应。本文船型网箱在0°、30°、60°、90°、120°、150°、180°入射方向下的六自由度（垂荡、纵荡、横荡、横摇、纵摇、艏摇）运动响应传递函数如图6 所示。

图6 运动响应传递函数RAOFig.6 Motion response transfer function RAO:a) heave;b) roll;c) surge;d) pitch;e) sway;f) yaw

由图6 可见，对于船型网箱，六自由度运动中垂荡运动是其中最激烈的运动模式。随着波浪频率逐渐增大，垂荡运动幅值先增大、后快速减小到0 附近。这可能是由于在波浪频率等于1.1 rad/s 附近时，网箱出现共振现象，导致垂荡运动的迅速增大。共振区的左侧区域，频率较小，波周期较大，结构所承受的浮力发生变化，出现了明显的回复力。对于纵荡以及横荡运动，基本都呈现出随着波浪频率的减小、波浪周期的增大，运动幅值不断增大的趋势，这主要是由于水动力分析时未考虑系泊力的影响，无法提供回复力以抵消结构运动导致的。同时可以看出，对于纵荡运动，波浪角度为0°以及180°时，运动幅值最大。因为此时波浪方向与船型网箱的长度方向一致。横荡的运动幅值则在90°时最大，因为此时波浪方向正好与船型网箱长度方向垂直，承受着最大的横向波浪力。

横摇与纵摇的运动幅值随着波浪频率变化趋势基本一致，均是先增大到一定程度后逐渐减小。这是因为横摇与纵摇均会受到回复力的影响，在超过共振区之后，运动幅值受到回复力的作用逐渐减小。同时，纵摇与纵荡一致，运动幅值在0°及180°时最大，而在90°时最小。横摇则相反，运动幅值在90°时最大，而在0°和180°时最小。艏摇运动是六自由度运动最为微弱的运动，并且运动幅值小，基本不发生运动。

2.2 结构应力分析

网箱在服役过程中，始终遭受着复杂的静水压力与波动载荷的作用，并最终反映到结构应力上。对结构进行应力分析有利于了解网箱整体的受力情况，判断结构整体强度，评估最危险工况与结构上最危险的节点等。同时，局部节点疲劳强度分析等工作的进行也需要应力分析的支持，特别是局部应力分析。因此，本节对网箱模型在波浪载荷作用下的应力进行分析研究。

2.2.1 整体应力分析

网箱在静水压力与单位波幅波浪力共同作用下的整体应力分析结果如图7 所示。由于结构的对称性，仅对0°、30°、60°、90°波浪入射方向进行展示。从图7 中可以看出，网箱整体上应力分布并不均匀，呈现出纵向箱型梁与船艏底部应力高，而艉浮箱与横向箱型梁的应力偏低的趋势。同时，尽管网箱整体应力并不高，但是在各波浪入射角下，均出现了明显的局部节点应力集中现象。如船中角隅处、船艏与两侧箱型梁连接处、上下箱型梁与之间斜撑的连接点、横纵箱型梁连接点等局部区域都是应力集中明显的区域，这些区域都是潜在的裂纹源或者疲劳损伤严重的区域。

图7 网箱整体应力响应Fig.7 Overall stress response of cages

2.2.2 局部应力分析

利用子模型技术得到的网箱局部应力分析结果如图8 所示，其中节点D1位于船艏浮箱与横箱型梁相贯处，节点D2和D3位于船中角隅处。从图8 中可以看出，通过子模型技术得到的局部应力与通过整体分析得到的相同位置处的节点应力在应力分布与应力值上大体一致，这说明子模型技术能够比较真实完整地继承整体模型的载荷条件。同时，由于子模型技术采用了更加精细的网格进行单元划分，在各拐角处的应力集中现象均更加显著。由此可见，利用子模型技术进行网箱上局部结构的应力分析能够更好地反映局部结构的真实应力分布，同时能够更好地突出应力集中现象，对疲劳热点的疲劳强度分析、疲劳寿命预测具有重要的意义。

图8 网箱局部应力响应Fig.8 Local stress response of cages

壳-实体耦合子模型与壳模型在相同载荷条件下的位移云图对比如图9 所示。可以看出，壳-实体耦合模型的计算结果与壳模型计算结果的一致性较好，证明了利用多尺度子模型技术建立壳-实体耦合子模型的可行性。

图9 多尺度子模型计算效果Fig.9 Result of multi-scale submodel:a) shell-solid coupling model;b) shell model

2.3 裂纹扩展和疲劳寿命

表面裂纹一般用半椭圆描述，因此本文也利用双参数半椭圆裂纹进行裂纹的描述。裂纹长度为2a，裂纹深度为c，裂纹长深比为a/c。基于此，在疲劳裂纹扩展分析中，仅需要针对裂纹深度与裂纹长度进行应力强度因子传递函数与裂纹增量的计算。本文将初始裂纹定义为a=c=0.5 mm 的圆形裂纹。

基于断裂力学的疲劳裂纹扩展主要流程：

1）初始化裂纹半长、深度、循环次数与时间。

2）在壳-实体耦合模型中，逐次施加某一浪向下所有不同频率的单位波幅规则波对应的边界条件，插入当前尺寸裂纹，并计算对应的SIF，以此构建该裂纹尺寸、该浪向下的应力强度因子传递函数。

3）根据应力强度因子传递函数，结合各短期海况的波浪谱，进行各短期海况下裂纹扩展增量的计算，并进行叠加。

4）当该浪向下所有短期海况均已循环完毕或者达到裂纹更新阈值时，进行裂纹尺寸更新，并重复步骤2）—4）。

D1、D2、D3这3 个节点在裂纹扩展过程中裂纹半长及裂纹深度随时间变化的曲线如图10 所示。可以看出，3 个节点的裂纹扩展速率存在着明显的差别，节点D1的扩展速率最高，其次分别为节点D2、D3。对应地，按照裂纹深度达到80%板厚时认为结构达到疲劳失效的判定标准，节点D1的疲劳寿命最短，而节点D3的疲劳寿命最长，三者疲劳寿命分别为11.3、17.3、33.1 a。

3 结论

本文以船型网箱为例，进行了海洋结构物频域内的响应分析，介绍了网箱的主要参数、参考坐标系和有限元模型建模与边界条件，同时开展了结构的运动特性与应力分析，并对疲劳热点开展疲劳分析。可得到以下主要结论：

1）通过水动力分析得到了该船体网箱的运动特性，为疲劳寿命评估提供了更准确的波浪载荷与结构响应。六自由度的运动幅值均与浪向和波频有着明显的关联，对于有明显共振频率的运动，其幅值往往在共振频率附近达到最大值，而对于没有共振频率的运动，其幅值往往随着波频的减小而增大。垂荡是六自由度运动中最激烈的运动模式，在0°与180°波浪入射角、1.1 rad/s 的频率下运动响应最激烈。纵摇和纵荡的运动幅值也在0°以及180°时最大。

2）网箱整体应力符合强度要求，但在不同浪向下均可观察到明显的应力集中现象，应力集中区域主要为船中角隅处、船艏与两侧箱型梁连接处、上下箱型梁与之间斜撑的连接点、横纵箱型梁连接点等。

3）子模型技术能够完整准确地继承局部区域的载荷条件并展开应力分析，能够更准确地反映局部区域的应力集中现象，为疲劳强度分析提供更准确的载荷条件。

4）实现了基于子模型技术的海洋结构物强度分析与寿命预测。运用子模型进行局部应力计算，并将多尺度子模型技术应用到其中，解决了大尺度模型与小尺寸裂纹之间的耦合、模型与载荷的传递、计算效率与计算精度的平衡等问题。