0 引　言

焊接技术因其良好的连接工艺性能，被广泛应用于钢质船舶建造。焊接过程中的冷热不均，将使得焊后船体构件完全冷却后，不可避免地产生焊接残余应力与变形^[1]。船体构件焊缝及其相邻区域存在的残余应力与变形不仅会影响船体结构的装配精度，并且对其承载性能有重要影响。目前，极限强度评估理念已成为船体结构强度评估的重要方面。诸多学者采用理论、数值模拟等技术手段对船体结构的极限强度开展了全面研究。Paik等^[2]及Cui等^[3]采用弹性大挠度理论与刚塑性分析相结合的方法对船舶板极限强度进行了研究，提出船舶板极限强度计算的简化公式。Hu等^[4]提出考虑焊接残余应力在单轴压缩载荷作用下，矩形板的极限强度近似方法。Ueda等^[5]研究了复杂初始挠度模式对矩形板极限强度的影响，提出分别用于预测薄板、厚板压缩极限强度的挠度法和曲率法。Gordo等^[6]推导了评估矩形板在双轴压缩载荷下的极限强度公式。Paik等^[7]提出在纵向载荷、边缘剪切和侧向压力组合载荷作用下钢板的极限强度公式。Sadovsk等^[8]研究了初始挠度幅值对矩形板压缩极限强度的影响。Paik等^[9]使用非线性有限元方法，研究了焊接残余应力对单轴压缩载荷作用下高强度钢板极限强度的影响。冯亮等^[10]研究了板格材料、初始缺陷、网格密度、边界条件等几种因素的敏感性，发现理想应力-应变关系会使得结果偏于危险。网格形状、网格密度和边界条件对于有限元结果均有影响。需强调的是，虽意识到焊接残余应力会随船舶运营中的循环载荷作用得到一定程度释放。但从不断追求船体结构极限强度评估准确性的角度看，在船体结构极限强度评估中更为细致地考虑焊接残余应力有必要。因为对于船体结构极限强度评估而言，焊接导致的残余应力与初始变形始终是确定各种评估模型的先决条件。

目前船体结构极限强度研究者对构件初始缺陷的处理，普遍采用的是基于测试数据回归公式，细节性描述可参阅Yao与Fujukubo的专著^[11]。随着计算机软、硬件水平的进步，结构焊接过程数值模拟技术取得了长足进步。为了能动态追踪焊接过程中的应力、应变和焊后的残余应力和变形，Ueda等^[12]提出了热弹塑性有限元法。随后，众多学者基于热弹塑性有限元法对钢结构焊接过程开展了广泛研究。Mollicone等^[13]研究了对接焊矩形板焊缝和夹持条件对面外变形的影响。Ranjbarnode等^[14]考虑热输入的影响，模拟了碳钢与不锈钢对接焊矩形板的焊接过程。为了对焊接过程进行更精确的数值模拟，Okano等^[15]考虑到电弧等离子熔深，研究了焊接工艺条件对焊缝残余应力和变形的影响，开发了基于电弧物理的热源模型。对于中厚板全熔透焊接工艺，Singh等^[16]提出一种自适应体积热源用于焊接过程模拟。Attarha等^[17]使用Abaqus软件进行了三维有限元模拟，预测整个板的温度分布，发现了焊接接头中没有填充材料有助于观察峰值温度和冷却斜率差异与材料性能差异的关系。Sun等^[18]对激光焊接和CO₂气体弧焊薄板接头的焊接引起的应力和变形进行了对比研究，发现激光焊接技术可获得较小的应力和变形。Choobi等^[19]研究了夹紧和夹紧时间对不锈钢板对接焊中残余应力和变形的影响，发现夹紧和松开时间对残余应力分布和最终角变形有较大影响。Chen等^[20]确定了船舶铝板热影响区尺寸、温度历史、变形和残余应力分布。崔虎威等^[21]研究了焊接电流、电弧有效加热半径、焊接速度等焊接参数对平板对接焊温度场和残余应力场的影响。

本文将热弹塑性有限元焊接仿真技术与非线性有限元极限强度数值模拟技术联合，基于热弹塑性有限元焊接模拟所获焊接残余应力与变形，对接焊矩形板极限强度进行分析。所有数值仿真工作，均在有限元软件Abaqus平台上开展。即首先开展对接焊矩形板热弹塑性有限元数值模拟，直接生成矩形板焊接残余应力和变形，并将本文所选厚度的矩形板纵向残余应力和垂向挠度进行比较。随后，对含有焊接残余应力与变形的对接焊矩形板模型，开展压缩载荷作用下的极限强度非线性有限元数值模拟。最后，将数值模拟结果与Faulkner公式计算结果进行对比。

1 对接焊船体矩形板数值模拟 1.1 船体矩形板几何尺寸及材料属性

计算模型几何尺寸选用ISSC所给出的船体板尺寸，具体为2 550 mm×850 mm×9.5 mm，几何模型如图1所示。图中，a为板长，b为板宽，a_c为焊缝宽度。将图中几何模型以焊缝为基准划分为焊缝区1、过渡区（2和3）及远离焊缝区（4和5）。材料采用EH36高强度钢^[22]，材料属性包括热传导系数、比热容、弹性模量、热膨胀系数、屈服应力等，各参数随温度变化情况如图2所示。

1.2 单元属性、网格尺寸与边界条件

在热弹塑性有限元分析中，分别采用壳单元与实体单元离散对接焊船体矩形板获得有限元模型。壳单元模型截面属性厚度为2.375 mm，沿板厚方向有5个截面点，与实体单元模型节点位置相对应。表1为2种模型在不同分析中的单元类型及网格尺寸。为了兼顾计算精度与效率，焊缝附近区域有限元网格划分较为密集，远离焊缝区域的网格相对稀疏。为避免模型发生刚体位移，且允许结构在焊接过程中自由收缩，在有限元模型4个角点（A、B、C、D）处施加位移约束。热弹塑性有限元模型如图3所示。

1.3 热源模型

采用双椭球移动热源子程序进行焊接数值模拟，双椭球热源将总的热能输入分配前半椭球部分和后半椭球部分如图4所示^[23]。

前半部分椭球内热源分布函数：

${{{q}}_1}(x,y,z) =\displaystyle \frac{{6\sqrt 3 {{{f}}_1}\eta UI}}{{{C_f}{{ab}}\text{π} \sqrt {\text{π}} }} \cdot \exp \left( { - \displaystyle\frac{{3{x^2}}}{{{C_f}^2}} - \displaystyle \frac{{3{y^2}}}{{{a^2}}} -\displaystyle \frac{{3{z^2}}}{{{b^2}}}} \right)，$

(1)

后半部分椭球内热源分布函数：

${{{q}}_2}(x,y,z) =\displaystyle \frac{{6\sqrt 3 {{{f}}_2}\eta UI}}{{{C_r}{{ab}}\text{π} \sqrt {\text{π}} }} \cdot \exp \left( { - \displaystyle \frac{{3{x^2}}}{{{C_r}^2}} -\displaystyle \frac{{3{y^2}}}{{{a^2}}} -\displaystyle \frac{{3{z^2}}}{{{b^2}}}} \right)。$

(2)

式中：q₁(x, y, z)和q₂(x, y, z)分别为前半椭球和后半椭球所在空间内的任意位置的热流密度；a₁为热源前端椭球半轴；a₂为热源后端椭球半轴；b、c分别为沿厚度方向的半轴即垂直于焊缝方向的半轴；f₁、f₂分别为前半部分和后半部分的分配系数，应满足f₁+f₂=2；η为焊枪热输入效率；U为焊接电压；I为焊接电流；$v$为焊接速度。在焊接时，焊接方法选用氩弧焊焊接，热源沿焊缝中线施加，焊接工艺参数^[24]见表2。

1.4 焊接过程热弹塑性有限元数值模拟 1.4.1 焊接温度场

采用热弹塑性有限元法，分别对壳单元、实体单元模型开展焊接过程数值模拟。提取了焊接温度、纵向焊接残余应力与垂向挠度信息。图5和图6分别为壳单元、实体单元模型的焊缝上、下表面的温度历程。图7和图8分别为2种模型在稳态阶段时的温度分布图。可以发现，焊缝处2个模型各层的温度分布基本一致,壳、实体单元模型的最高温度均超过材料熔点1500 ℃。

1.4.2 纵向焊接残余应力

横向焊接残余应力对纵向受压矩形板极限强度影响较小，本文主要关注矩形板的纵向焊接残余应力。图9为壳单元模型的纵向焊接残余应力云图，图10为壳、实体单元模型沿中线E-F路径上的纵向焊接残余应力对比，表明2种模型焊后纵向焊接残余应力一致。

1.4.3 垂向挠度

图11为壳、实体单元模型沿中线E-F路径垂向挠度，表明2种模型的垂向挠度变化规律一致，壳单元模型垂向挠度比实体单元模型大。

1.4.4 计算效率

表2为同一软硬件环境下，壳、实体单元模型的单元与节点数与热分析与应力分析的计算时间。具体的计算机配置为：64位Windows操作系统，R7 4800H 处理器，CPU主频为2.90 GHz，内存为16 GB。从表3可知：壳单元模型能够有效地减少单元和节点数目，计算效率较实体单元模型显著提高。

2 极限强度非线性有限元数值模拟与Faulkner简化公式计算

前文通过热弹塑性有限元法，获得了对接焊矩形板的焊接残余应力和变形，可作为极限强度非线性有限元数值模拟的初始缺陷。同时，矩形板材料属性选用EH36高强度钢在室温下的材料参数。为便于分析本文所采用计算路线的可行性，计算矩形板极限强度的有限元模型均采用壳单元模型，一并给出Faulkner简化公式计算结果。Faulkner简化为：

$\frac{{{\sigma _u}}}{{{\sigma _Y}}} = \frac{2}{\beta } - \frac{1}{{{\beta ^2}}}。$

(3)

式中：${\sigma _u}$为极限强度；${\sigma _Y}$为材料屈服强度；$\beta$为板的柔度系数，且$\beta = b/t\sqrt {{\sigma _Y}/E}$。

2.1 极限强度非线性有限元模型

本文采用2块不同厚度的矩形板作为极限强度分析模型，如图12所示。在极限强度分析中，矩形板两端部受纵向压缩荷载作用。矩形板极限强度分析模型几何尺寸及边界条件见表4。

2.2 初始状态时刻纵向残余应力和垂向挠度

图13为初始状态时刻，厚度为9.5 mm和11 mm的对接焊矩形板沿中线E-F路径上的纵向焊接残余应力，对比发现2块矩形板的纵向焊接残余应力变化规律一致。图14为初始状态时刻，厚度为9.5 mm和11 mm的矩形板沿中线E-F路径垂向挠度，由图可知，2种厚度矩形板的垂向挠度变化规律一致，板厚为9.5 mm的矩形板垂向挠度峰值较大。

2.3 极限状态Von-Mises应力云图

图15和图16分别为厚度为9.5 mm和11 mm的矩形板极限状态Von-Mises应力云图，可知：Von-Mises应力沿焊缝对称分布，在焊缝及其热影响区附近，其Von-Mises最大应力值接近钢材屈服强度。

2.4 应力-应变曲线

图17为厚度为9.5 mm、11 mm的对接焊矩形板应力-应变曲线，曲线中的峰值点所对应的应力即为矩形板的极限强度。可知：板厚为9.5 mm矩形板A₁的最大峰值点对应的应力值为173.3 MPa, Faulkner公式计算的极限强度值为166.1 MPa；板厚为11 mm矩形板B₁的最大峰值点对应的应力值为197.0 MPa，Faulkner公式计算的极限强度值为187.6 MPa。

3 结　语

本文首先对对接焊船体矩形板进行热弹塑性有限元分析，对比了相同几何尺寸的壳、实体单元模型的温度场、应力场、变形等结果，并对比了2种模型的计算效率。继而对带有焊接残余应力与变形的壳单元模型，开展了压缩载荷下的极限强度非线性有限元数值模拟。最后，分别比较了2块不同厚度对接焊矩形板初始时刻纵向残余应力和垂向挠度分布规律、极限强度数值模拟结果与Faulkner公式计算结果，得到以下结论：

1）在保证一定计算精度的前提条件下，在热弹塑性有限元数值模拟中采用壳单元模型的计算效率较实体单元模型有显著提高。

2）在同一焊接工艺参数情况下，板厚对纵向焊接残余应力分布影响较小。但厚度为9.5 mm的矩形板垂向挠度峰值比厚度为11 mm的矩形板的垂向挠度峰值大。

3）不同厚度矩形板通过数值模拟得到的极限强度值，与Faulkner简化公式计算的极限强度值误差分别为4.2％和5.1％，因此数值模拟值与公式计算值较吻合。

4）基于热弹塑性有限元法，获得的焊接残余应力和变形，可作为对接焊矩形板极限强度非线性有限元数值模拟的初始缺陷。