三维网格重划分技术及其在焊接数值模拟中的应用

黄辉，赵耀，袁华，王波

（华中科技大学船舶与海洋工程学院，湖北武汉 430074）

0 引言

焊接是现代机械工业中的重要连接工艺，是结构产生残余变形和应力的主要来源。从20世纪50年代开始，人们开始焊接变形和应力的数值计算，取得了一系列的研究成果。近10多年来，随着计算机的迅猛发展，使焊接过程的热弹塑性计算成为模拟焊接变形、应力、组织转变等复杂现象的主要手段。焊接过程中材料在高温条件下发生分解、相变，而且不同温度下的材料性能具有很大差异，使焊接局部区域经历高度非线性的过程，要模拟整个结构的焊接变形和残余应力，需要进行热机耦合非线性有限元分析。由于焊缝区温度和应力应变的梯度很大，进行有限元建模分析时，网格划分必须足够精细，并且进行时间增量步的分析。对于大型复杂结构，其焊缝长度和焊接时间都很长，有限元网格规模和计算时间步势必急剧膨胀，进行完整的热弹塑性有限元分析将十分困难。

网格重划分技术可以处理网格畸变及单元变形过大的问题，使热轧、岩土等大变形分析得以进行下去，并保证计算结果的精度和可靠性。S．Brown和H．Song［1］结合网格重划分和动态子结构技术，采用壳体单元分析了平板激光焊接后的变形和应力，所使用的时间与完整的热弹塑性计算相比，仅为后者的七分之一。L．－E．Lindgren［2］基于一种可变节点的体单元，开发了自适应网格技术，并应用三维有限元方法成功模拟了大型铜罐的电子束焊接，大幅缩短了计算的时间。本文开展了三维全六面体网格重划分技术研究并将其应用于焊接数值模拟，对平板和圆筒的焊接模型分别进行了焊接热力学仿真。

1 网格重划分

焊接过程中，热影响区的温度、应力应变梯度较大，具有高度的非线性特征，有限元模拟中单元应当划分得很细;然而，远离焊缝的绝大部分区域仍然处于线弹性状态，单元可以划分得相对较粗。如果不断对焊接模型进行网格重新划分，使细密的网格随着热源一起移动，就能节省大量的计算资源并保证计算结果的精度。在网格重划分前后，必须保证新旧网格的节点占据相同几何空间，以顺利地进行应力应变等变量的映射，使分析求解得以连续进行下去。下面介绍节点迁移和变量映射2个步骤。

1.1 节点迁移

由于网格疏密程度在网格重划分前后变化较大，不得不完全改变节点和单元的编号，重新构建有限元模型的网格。为了匹配新网格节点在旧网格中的位置，必须对网格单元一一寻址，逐一判断节点是否处于某单元内部［3］。节点迁移的基本问题就是已知单元8个节点的坐标和位于该单元中1个节点的坐标，反求这1个节点在该单元对应的等参元中的坐标。

对于六面体单元，其形函数为:

根据等参元的性质，单元中任一点处的整体坐标

将表1中坐标值代入式（1）和式（2），可得如下方程组:

式中:ai，bi，ci为已知量xi，yi，zi的函数，它们的值均为常数。上述高阶齐次方程组的求解方法有牛顿法、割线法、布朗方法等，本研究采用了牛顿法，对于比较规则的网格单元，进行为数不多的几次迭代其结果就能收敛。求得节点在局部坐标系中的坐标后，可以用形函数对节点位移进行插值，从而得到下一步分析时的节点在整体坐标系中的坐标值。

1.2 变量映射

确定新网格的节点坐标后，需要将旧网格中求解的应力应变场映射到新网格上，以保证分析的连续性和准确性。采用了如下映射算法［4］:

1）采用插值技术，将上一步网格中积分点上的变量（应力、应变等）外插到每一个单元的节点上，然后对所有邻接同一节点的单元的变量取平均。

2）新网格的积分点的位置由旧网格积分点来确定，分两步:①先找到该点所属旧网格的单元，这样就可以确定该点在单元中的位置（该步骤假定新网格中所有积分点位于旧网格边界内）。②变量由旧网格的节点上内插到新网格的积分点上。

2 计算实例

分析了平板表面堆焊和圆筒环焊缝的焊接过程，分别比较了完整模型与网格重划分模型的计算结果。平板堆焊模型的尺寸为500 mm×300 mm×6 mm（图1）。圆筒的内径为150 mm，外径为156 mm，长度为400 mm（图2）。平板和圆筒的焊接条件见表2。SS400的材料热物性与热力学性能参数曲线见文献［5］。

2.1 热传导计算

采用顺次耦合的热弹塑性有限元计算方法，分别对平板模型和圆筒模型进行焊接热力学仿真分析。首先进行焊接过程的热传导计算，把得到模型的瞬态温度数据信息存储在1个文件中，然后读入每一时刻的温度作为模型的载荷，进行焊接模型的应力应变计算。由于高斯热源模型能很好地反映大多数焊接过程中的热流分布特征，其在焊接热传导的数值仿真中应用甚广。高斯热流分布示意如图4所示，函数表达式为:

图3 焊缝中线处的圆筒横剖面Fig.3Cross section of pipe model

式中:q（r）为半径r处的表面热流，W/m2;q（O）为热源中心处的热流量最大值;c为热源集中系数;r为距热源中心的距离。

图4 高斯热源热流分布示意图Fig.4Heat flux distribution of Gauss heat source

热传导计算中考虑了相变潜热，取为270 kJ/kg，固相线温度为1 465℃，液相线温度为1 544℃;环境温度设为25℃，冷却条件为空冷，表面换热系数取常数10 W/（m2·℃）;考虑辐射散热，发射率取为0.5。加热过程中时间步长最大取0.5 s，冷却过程取自动时间步长。为了兼顾计算的效率和精度，在焊缝附近区域网格划分得很细，在远离焊缝区划分得较粗。沿焊缝方向根据热源的位置，在网格重划分模型中分4步实现焊接过程模拟，平板焊接模型的有限元网格如图5所示。完整网格模型的有限元网格与网格重划分模型第4步的网格相同，沿纵向均匀划分，沿横向由密变疏。圆筒焊接模型的网格划分方法与重划分步骤与平板焊接模型相似，在此不赘述。为了简化，以FM（Full Model）表示完整模型，RM（Rezoning Model）表示网格重划分模型。平板完整模型、平板网格重划分模型、圆筒完整模型、圆筒重划分模型分别表示为FMⅠ，RMⅠ，FMⅡ，RMⅡ。平板和圆筒焊接过程中的瞬时温度场分布云图如图5和图6所示。

2.2 应力应变计算

应力应变计算所使用的有限元模型与热传导计算模型的节点和单元完全相同，单元类型由热单元改为力学单元，材料特性由热物性参数变为热力学性能参数，同时增加了力学边界。在平板有限元模型中只限制了结构的刚体位移，在圆筒有限元模型中设置了一端简支的约束条件。将热传导计算的节点瞬时温度读入有限元模型，开始应力应变的增量步计算。图7为圆筒焊接后的Mises应力分布云图。

图7 圆筒焊后Mises应力云图Fig.7Mises Stress distribution of pipe model after welding

3 计算结果

3.1 温度场计算结果

平板焊缝长度中心处横截面上3点的瞬时温度结果见图8，测点分布于平板表面，与焊缝的距离分别为0 mm，20 mm和44 mm。从图中不难看出，焊接过程中各点经历了加热和冷却2个阶段。距离焊缝越近的点，温度上升过程越快，达到的最高温度也越高。在冷却阶段，各点的温度逐渐趋于一致。网格重划分模型与完整模型的计算结果完全一致，说明采用网格重划分技术在计算热传导过程时满足精度要求。

3.2 变形场计算结果

平板焊接模型的变形结果见图9，曲线表示了横向收缩、角变形沿焊缝方向的分布情况。从图中可以看出，横向收缩在焊缝长度中间段有最大值，角变形随焊缝长度增加略有变大，但整体上看，2种变形沿焊缝方向都较为均匀。

=圆筒焊接模型中A和B两点（标示于图3）的瞬态径向位移示于图10。从图中可以明显看出，当点所在区域处于加热阶段时，焊缝金属发生膨胀，径向位移变大，在热源远离该处，焊缝因冷却而收缩，这时径向位移减小。由于A，B两点在环焊缝中受热2次，可在位移图中看到2个波峰。

由图9和图10可知，平板和圆筒的完整模型与网格重划分模型计算结果相一致，说明网格重划分模型预测的焊接变形结果没有精度损失。

3.3 应力场计算结果

平板表面堆焊模型的残余应力结果以等值线的方式绘制见图11和图12。由图11可见，纵向残余应力的分布情况，焊缝区域存在较大的拉应力，远离焊缝区则呈现较低的压应力以使构件保持自身的平衡。图12为平板焊后的横向残余应力结果，在焊缝两端存在较大的压应力，焊缝区的横向应力水平较低。图11和图12表明，完整计算模型与网格重划分模型计算的焊接残余应力结果非常吻合。

圆筒焊接模型的0°和180°纵向剖面的残余应力结果见图13和14，图中反映了圆筒焊后内外侧的应力分布情况。焊缝起始位置处（0°）的环向应力峰值较低，而焊缝中间段（180°）的环向应力峰值则超过了材料的屈服强度;同时，圆筒内侧的应力要高于外侧的应力水平。图14中的结果表明，圆筒内侧的轴向应力为拉应力，外侧的轴向应力为压应力，0°位置处的应力峰值低于180°位置处的应力水平。从图中结果可以看出，完整模型和网格重划分模型的横向残余应力和纵向残余应力完全吻合，而且有限元分析预测的环向应力和轴向应力的分布规律与Burdekin［6］关于低碳钢圆筒环焊缝的研究结果相一致。

3.4 有限元计算时间

平板与圆筒焊接模型的有限元网格节点数和计算时间如表3所示，2种模型在使用网格重划分技术后，计算时间均减少了约30%。由于网格重划分模型中的节点数一般少于完整模型的节点数，也即模型中总的自由度数前者少于后者，内存的需求也必将因网格重划分技术而得以减少。对大型结构而言，计算时间和内存大小是决定数值模拟能否实现的主要因素。计算结果表明，网格重划分技术可以解决这方面的问题。本文暂时没有研究热源后端的网格粗化问题，若能选择合理的粗化准则并优化网格设计将可能大幅提高计算效率，同时保证结果的可靠性。

4 结语

1）焊接热源影响下的区域具有高度非线性特征，远离热源区域为弱非线性，可以采用网格重划分技术进行模拟。

2）使焊缝区细密网格随热源不断移动，与完整的有限元模型相比，减少了节点自由度，提高了计算效率并保证了应力和变形结果的精度。

3）通过网格重划分模型与完整模型的计算结果可知，提出的网格重划分方法在充分保证结果精度的前提下，大约减少了30%的计算时间。

［1］BROWN SB，SONGH．Rezoninganddynamic substructuring techniques in FEM simulations of welding processes［J］．ASME Journal of Engineering for Industry，1993，155:415－423．

［2］LINDGREN L E，HAGGBLAD H A，et al．Automatic remeshing for three-dimensional finite element simulation 0f welding［J］．Computer Methods in Applied Mechanics and Engineering，1997，147（3）:401－409．

［3］MURTI V，VALLIAPPANS．Numericalinverseisoparametricmappinginremeshingandnodalquantity contouring［J］．Computers and Structures，1986，22（6）: 1011－1021．

［4］HIBBITT D，KARLSONB，SORENSONP．ABAQUS analysisuser'smanual［M］．ABAQUSVersion6.5 Documentation．

［5］WANG Rui，SHERIF R，et al．Numerical and experimental investigations on welding deformation［J］．Trans．JWRI，2008，37（1）:79－90．

［6］BURDEKIN F M．Local stress relief of circumferential butt welds in cylinders［J］．British Welding Journal，1963，10 （9）:483－490．