分享：SUS316不锈钢马鞍形管-管接头的残余应力数值模拟及高效计算方法开发

骆文泽, 胡龙, 邓德安,

奥氏体不锈钢具有优越的综合力学性能与抗腐蚀性能,因此成为核电装备中关键部件的最主要结构材料之一。在制造过程中,奥氏体不锈钢部件通常采用熔焊方法进行连接。近年来,国际原子能机构IAEA (International Atomic Energy Agency)的调查结果表明,在现有运行的400余套核电机组中,为数不少的焊接接头(包括奥氏体不锈钢接头)中检测出了应力腐蚀裂纹[1,2]。已有研究[3~5]表明,除了材料敏化和腐蚀环境外,拉伸应力被认为是导致奥氏体不锈钢接头产生应力腐蚀裂纹的重要因素之一。应力腐蚀裂纹的产生会给核电设备的安全稳定运行带来严重的威胁。焊接制造过程局部的快速加热和冷却会产生很高的残余应力,其峰值在焊缝及热影响区往往会达到甚至超过材料的屈服极限,这会导致这些区域更易产生应力腐蚀裂纹。到目前为止,尽管有10多种实验方法(包括无损检测方法、半有损法和有损检测法)可以用于获得焊接接头与结构的残余应力,但是由于实验手段自身的局限性,采用这些方法只能获得有限位置的应力分布情况,很难全面地获得整个焊接接头或结构的残余应力分布。1971年,Ueda和Yamakawa[6]首次提出了热-弹-塑性有限元方法来计算焊接残余应力,基于热-弹-塑性有限元方法的数值模拟技术弥补了实验手段的不足,这种方法可以获得整个焊接接头或结构的残余应力分布与大小。进入21世纪后,该方法已经广泛应用于实际工程结构中焊接接头和结构的残余应力预测[7~9]。但是,由于焊接过程具有非稳态、非线性以及多物理场交互作用的特点,若再加上焊接结构的接头数量多、结构复杂以及尺寸厚大等特点,在采用热-弹-塑性有限元方法计算焊接残余应力时需要耗费更多的计算资源和更长的计算时间,很多情况下甚至会超出计算硬件能力范围而不能实施计算。

针对如何以现有商用有限元软件为平台来高效、高精度地计算焊接残余应力的工程需求,一些学者对此展开了研究。戴培元等[10]采用三维全尺寸模型和二维轴对称模型分别计算了管-管对接接头的焊接残余应力,通过对比发现,二维轴对称模型也能较高精度地预测焊接残余应力,且能大幅节约计算时间和存储空间。Pu等[11]分别采用瞬间热源和移动热源对厚大平板对接接头进行了焊接残余应力的数值模拟,计算结果表明,瞬间热源除了具备耗时少的优点外,其残余应力的计算结果与移动热源的计算结果也有较好的吻合度。胡兴等[12]在对厚大板材进行多层多道焊时,采用合并焊道的方法来计算接头的焊接残余应力,该方法不仅可以节约计算时间,而且得到的残余应力与实验结果也吻合较好。

在实际工程中如果遇到厚大复杂结构,且不能进行降维(如：由三维简化为二维)简化时,往往需要在具体问题具体分析的基础上提出能兼顾计算时间与精度的方法来获取焊接结构的残余应力。本工作的研究对象是SUS316马鞍形管-管接头,这种结构是压水堆核电站中第二水循环系统中的典型结构。该接头具有尺寸大、管壁厚和焊道多的特点,且马鞍形焊缝的坡口大小及局部拘束条件随位置的不同而发生变化。理论上,这类马鞍形接头不适合采用降维简化法来计算焊接残余应力。到目前为止,关于马鞍形管-管复杂接头的残余应力数值模拟的研究并不多见,以往研究者们研究的马鞍形管-管焊接接头尺寸都相对较小,且焊道数量较少。张敏等[13]采用生死单元方法对6道焊的马鞍形管-管接头进行了数值计算,但是该研究对结果的报道相对简略,甚至未给出整个焊接接头的应力分布云图。Katsuyama等[14]采用移动热源对尺寸较小的马鞍形管-管接头进行了数值模拟,但在该研究中,作者并未考虑材料的加工硬化与退火软化效应,也没有将计算结果与实测值进行比较。从所能查阅的公开文献来看,对于厚大的马鞍形管-管接头焊接残余应力的数值模拟和实验测量的比较未见报道。

针对SUS316马鞍形管-管厚大焊接接头,为了提高计算精度和计算效率,本工作以商用有限元软件MSC. Marc为平台,开发了2种焊接残余应力的计算方法。第1种方法建立与实际接头尺寸一致的尺寸全相同的模型,采用移动热源与瞬间热源混合使用的方法,兼顾效率与精度——即采用移动热源模拟打底与盖面焊道的热输入,采用瞬间热源模拟填充焊道热输入；第2种方法则利用接头模型几何形状的对称性,简化模型为1/4局部模型,并采用瞬间热源模拟全部焊道的热输入。本工作在得到SUS316马鞍形管-管接头整体残余应力分布的同时,对比了这2种方法计算得到的典型位置的焊接热循环和残余应力分布,同时,将计算结果与实验测量结果进行比较,并验证了所开发计算方法的可靠性。此外,还分析了这2种计算方法各自的优缺点及它们的工程适用性。

1 实验方法

1.1 焊接实验

SUS316马鞍形管-管垂直相贯接头(下称马鞍形管-管接头)由台管(图1a)和母管(图1b) 2个部分组成。两构件间的坡口为U形坡口,该U形坡口形状会随着构件位置变化发生变化,但是坡口的间距相差不大,在打底完成后,上表面坡口宽度范围在33.6~35.2 mm之间。马鞍形管-管接头的照片如图1c所示,其整体管尺寸如图1d所示。马鞍形管-管接头的特征尺寸如下：母管的长度为850.0 mm,内径为609.6 mm,壁厚为38.9 mm,台管上端面与母管的距离为296.1 mm,台管的外径为184.1 mm,壁厚为46.0 mm,鞍型焊道最外端部距母管端面距离为284.2 mm。台管为锻造结构件,它与母管的材料均为SUS316奥氏体不锈钢,焊接填充材料为Y316L,它们的化学成分如表1所示。图2是U形坡口实际焊接时的焊道布置,焊道布置为16层51道,焊接方法为钨极惰性气体保护焊,焊接参数如表2所示。焊接时,层间温度控制在100℃以下。在图2中,A、B、C和D 4个点是用于测量焊接过程温度循环的热电偶位置,其中A、B位置在母管的内表面上,它们分别与第一道焊缝中心线的距离为10和15 mm,C、D位置在母管的外表面上,它们与第46道焊趾的距离分别为5和10 mm。

图1

图1   SUS316台管、母管、管-管接头的照片及特征尺寸示意图

Fig.1   Photographs of a vice-tube (a), a circular pipe (b), and a saddle tube-pipe joint (c), the dimension diagram of the tube-pipe welded joint (unit: mm) (d)

表1   SUS316母材及Y316L焊材的化学成分 (mass fraction / %)

Table 1  Chemical compositions of SUS316 base metal and Y316L welding wire

新窗口打开| 下载CSV

图2

图2   焊道布置及热电偶位置示意图

Fig.2   Schematic of an arrangement of weld passes and locations of thermo-couples

表2   各焊道焊接速率、热输入及电弧效率

Table 2  Welding speed, heat input, and arc efficiency of each bead

新窗口打开| 下载CSV

1.2 应力测量方法及位置

焊接完成后,采用切片法[15]来获取马鞍形管-管接头90°与180°截面上、下表面的焊接残余应力。采用切片法测量残余应力时,需先将应变片贴于工件表面,然后采用线切割方法将贴片处附近的材料切成小片(本实验为4 mm × 4 mm × 4 mm),来获取释放的弹性应变,最后再通过理论计算得到残余应力[15]。应变片位置布置示意图如图3所示,图3a为俯视示意图,图3b为90°与180°横截面应变片布置示意图。为了方便后续讨论,90°上下表面测量位置与180°上下表面测量位置,分别用line 1、line 2、line 3和line 4表示,如图3b所示。

图3

图3   应变片位置俯视图及各截面应变片位置示意图

Fig.3   Schematics of strain gauge locations in top view (a), and strain gauge locations in 90° and 180° cross-section (b)

2 有限元计算方法

基于MSC. Marc软件平台,本工作建立了2种计算马鞍形管-管接头残余应力的有限元方法,为了方便后续讨论,这2种有限元计算模型分别用Model 1和Model 2表示。第1种方法(Model 1)建立了与实际焊接工件尺寸完全一致的三维有限元网格模型,如图4a所示。在模拟焊接热输入时,坡口下部打底焊道(第1条焊道)和紧接着的第2、3焊道及上部的第46~51盖面焊道采用移动热源模型,其余填充焊道均采用瞬间热源模型。第2种方法(Model 2)为了更进一步提高计算效率与节省计算储存空间,根据该接头几何形状的对称性简化为1/4局部模型,如图4b所示,并在对称面上设置绝热边界条件与对称位移拘束,所有焊道均采用瞬间热源进行温度场计算。Model 1与Model 2的焊道布置如图5所示。焊道采用“生死单元”的方式进行激活。采用移动热源模型时,当热源移动到当前材料位置时,对应的单元被依次激活；当采用瞬间热源时,当前焊道对应的所有单元同时被激活。

图4

图4   全尺寸(Model 1)和1/4尺寸(Model 2)有限元网格模型

Fig.4   Schematic finite element models of full size model (Model 1) (a) and a quarter model (Model 2) (b)

图5

图5   有限元模型中焊道布置

Fig.5   Arrangement of weld beads in the finite element models

2.1 温度场计算

在Model 1中采用移动热源来计算打底焊道和盖面焊道的温度场,本工作采用了Goldak双椭球热源模型[16]来模拟焊接热输入。在该模型中,前、后半部分的热流密度(qf和qr)分布分别定义为[16]：

式中,Q = ηUI (其中,Q为电弧功率(W),η为电弧热效率,U为电弧电压(V),I为焊接电流(A))；af、ar、b和c为椭球形状参数；ff和fr为前、后椭球热量分配系数,且ff + fr = 2,在本研究中,取ff = 0.6,fr = 1.4。

理论上而言,对于多层多道的焊接构件,所有焊道采用移动热源模型来模拟热输入可以获得更接近于实际情况的温度场、残余应力和变形。但是,本工作中的马鞍形管-管接头具有焊接道数多和焊道轨迹复杂的特点,如果全部采用移动热源需要消耗大量的计算资源和冗长的计算时间。为了解决上述矛盾,本研究在可变长度热源模型[17]的基础上,提出了“移动热源+瞬间热源”混合计算法,即采用移动热源模型来模拟打底焊道和盖面焊道的热输入,采用瞬间热源模型来模拟中间填充焊道的热输入。这样的处理方式,既可以大幅节省计算时间和存储空间,也可以兼顾计算精度。

本工作采用的“瞬间热源模型”是一个等密度体积热源[18],其热源长度为整条焊道的长度,横截面积与对应的焊道截面面积相等,加热方式为把整条焊道定义为内部热源并同时加热。该热源模型是在保证每条焊道的总热量与实际的总热量不变的情况下,通过调整加热时间th (s)和热流密度qv (J/(mm3·s2)),使得焊接熔池的最高温度控制在一个合理的范围以及熔池截面控制在一定面积内[19]。th为整条焊缝的假定加热时间,一般设定在数秒的时间范围内[20],qv由 式(3)[19]计算：

式中,L为整条焊缝的长度(mm),v为移动热源焊接速率(mm/s),V为该整条焊道的体积(mm3)。

2.2 应力场计算

应力与变形计算采用热-弹-塑性有限元计算法,即计算各个节点温度场后,再将其以热载荷加载到结构分析有限元模型中,从而求解整个焊接过程的位移、应变和应力,本工作采用的SUS316材料高温热物理性能和高温力学性能参数见文献[10]。

金属材料在焊接过程中热-冶金-力学的耦合行为非常复杂,在焊接热循环作用下,理论上材料的应变包括弹性应变、塑性应变、热应变、相变应变(固态组织转变引起的应变)和蠕变应变等成分。由于本工作采用的母材和焊接填充材料均为奥氏体不锈钢,不需要考虑固态相变问题。此外,焊接过程中加热和冷却速率都较快,这样在高温停留的累积时间相对较短,所以蠕变应变的影响也可以忽略不计[21]。对于SUS316钢,材料的总应变增量(Δεtotal)可以表达为：

式中,Δεe为弹性应变增量,Δεth为热应变增量,Δεp为塑性应变增量。

材料模型中,假定弹性行为遵循各向同性Hook定律[22],塑性行为遵循Von Mises屈服准则[22]。以往的研究[23]表明,SUS316奥氏体不锈钢是一种加工硬化较显著的材料,若不考虑加工硬化效应,得到的残余应力将远低于实际水平；若采用各向同性硬化准则考虑不锈钢的加工硬化,但不考虑退火软化效应,焊接数值模拟得到的应力结果会高于实验结果。因此本工作中,同时采用各向同性加工硬化模型和“阶跃式”退火软化模型考虑SUS316不锈钢的加工硬化和退火软化计算焊接残余应力。各向同性硬化法则中,流变应力与塑性应变(εp)有关,在模拟残余应力时,本工作采用的SUS316奥氏体不锈钢不同温度与不同塑性应变对应的流变应力如表3[24]所示。采用“阶跃式”模型来考虑退火软化效应时,退火温度设置为900℃[25]。

表3   SUS316不锈钢各向同性硬化准则模型参数[24]

Table 3  Parameters in the isotropic strain-hardening model for SUS316 stainless steel[24]

新窗口打开| 下载CSV

3 结果与分析

3.1 焊接热循环

为了比较分别由数值模拟和实验手段获得的焊接热循环,选取母管上4个代表位置(图2中A、B、C和D点)来输出热循环曲线,这4个位置均在母管的90°截面上。A、B位置在第1道焊接过程中的温度循环曲线如图6a所示,C、D位置在第46道焊接过程中的温度循环如图6b所示。从图6可以看出,尽管计算值与实验值在局部有些许差异,但整体而言,不论是各点的峰值温度还是升温和降温过程都与实验结果比较吻合。从数值模拟结果与实测值的比较可知,不管是移动热源模型+瞬间热源模型的Model 1,还是采用纯瞬间热源的Model 2计算得到的温度循环与实验结果均吻合较好。因此,可以推测,采用瞬间热源也能较高精度预测焊接残余应力,而且能大幅节约计算时间。

图6

图6   焊接温度循环结果和实验结果的比较

Fig.6   Comparison between simulated and experimental results for the welding temperature as a function of time

(a) the 1st pass (the measurement locations at points A and B (inset))

(b) the 46th pass (the measurement locations at points C and D (inset))

3.2 焊接残余应力

为了与实验结果进行定量比较和便于分析残余应力的分布特征,对残余应力计算结果进行了坐标变换,即将直角坐标系(x, y, z)转换为圆柱坐标(r, θ, z1) (其中,r方向为径向方向(应力实验中横向应力方向),θ方向为周向方向(应力实验中纵向应力方向),z1方向为轴向方向(为径向与周向共同决定平面的法向方向))。坐标变换及柱坐标原点位置如图7所示。柱坐标转化后,沿着r方向的应力为径向残余应力,沿着θ方向的残余应力为周向残余应力,z1方向的残余应力为轴向残余应力。

图7

图7   残余应力结果坐标系变换(直角坐标转柱坐标)

Fig.7   Schematic showing the transformation from the coordinate to cylindrical coordinate for the residual stress (r—direction of transverse residual stress, θ—direction of longitudinal residual stress, z1—the axial direction)

3.2.1 周向残余应力

图8为Model 1与Model 2的整体周向残余应力(类似平板接头的纵向残余应力)分布云图。总体上而言,这2个模型计算得到的周向残余应力的分布形态十分相似。在Model 1 (图8a)的外表面焊缝始终端位置(0°)可以清晰地看到该处的应力分布与其他位置有较显著的差异。而在Model 2 (图8b)中,由于全部焊道采用瞬间热源模型及对称的热边界和力学边界条件,因此沿着焊缝周向方向的残余应力,不会像Model 1在0°处出现不连续现象,故而呈现出光滑的连续分布特征。图9为Model 1在0°、90°、180°和270°截面上焊缝及其附近区域的周向残余应力分布云图。从4个截面上应力分布云图可以看到,焊缝及其附近的周向残余应力都为拉应力,最后一道焊缝位置及其附近区域的残余应力最大,该处峰值应力达到甚至超过600 MPa。由于盖面焊道采用了移动热源,加之各个截面处的内拘束不同,4个截面上的应力分布也有所不同,不过分布形态基本相似,峰值应力的大小和所在位置也基本相同。

图8

图8   焊接完成后Model 1与Model 2整体周向残余应力分布云图

Fig.8   Overall hoop residual stress distributions after welding for Model 1 (a) and Model 2 (b)

图9

图9   Model 1典型截面上的周向残余应力云图

Fig.9   Hoop residual stress distributions in the characteristic cross-sections for Model 1

(a) 0° (b) 90° (c) 180° (d) 270°

图10中的Model 2在90°和180°截面上的周向残余应力分布与图9b (Model 1的90°截面)和图9c (Model 1的180°截面)的周向残余应力的分布基本一致,这说明采用1/4部分模型也能较高精度地获得周向残余应力。

图10

图10   Model 2典型截面上的周向残余应力云图

Fig.10   Hoop residual stress distributions in the characteristic cross-sections for Model 2

(a) 90° (b) 180°

图11为管-管接头在90°截面的上表面(line 1)和下表面(line 2)的周向残余应力分布,该图比较了数值模拟结果与实验测量结果。在line 1上的应力分布如图11a所示,在-20~0 mm范围内,2种模型的计算结果与实验测量结果吻合很好。在-10 mm附近,周向残余应力达到最大值,2种模型计算得到的最大值和实验测量得到的最大值分别为534、547和506 MPa。Model 1的计算结果在0~15 mm范围内低于实验值,平均低了约40 MPa。Model 2的计算结果在15~25 mm范围内高于实验测量结果,平均高了约50 MPa；应力分布形态与测量结果基本一致。总体上来看,在line 1上数值模拟结果得到的应力分布与实验值一致。line 2上的周向残余应力分布如图11b所示,从此图可知,整体上而言实验值与计算值具有较高的吻合度。

图11

图11   90°处上(line 1)下(line 2)表面周向残余应力实验值与计算值

Fig.11   Experimental and numerical results of the hoop residual stress in 90° cross-section, located at line 1 (a) and line 2 (b) (Insets show the measuring directions)

图12为180°截面的上表面(line 3)和下表面(line 4)的周向残余应力分布,此图比较了数值模拟结果与实验测量结果。位于line 3上的应力分布如图12a所示,在-15~0 mm的区间,数值模拟结果与实验测量值吻合良好,而在0~15 mm的范围内实验值要明显高于数值模拟结果,而且Model 2的预测结果与实验值更接近一些。在line 3上,实验、Model 1和Model 2的最大值分别为529、543和530 MPa,3者的峰值非常接近。位于line 4上的周向残余应力分布如图12b所示,Model 1与Model 2的计算结果与测量结果的分布形态及数值的大小都吻合良好。

图12

图12   180°处上(line 3)下(line 4)表面周向残余应力实验值与计算值

Fig.12   Experimental and numerical results of the hoop residual stress in 180° cross-section, located at line 3 (a) and line 4 (b) (Insets show the measuring directions)

图11和12中的周向残余应力实验测量结果与数值计算结果存在一定程度的差异,但在应力峰值、最小值及分布形态上基本一致,整体上而言,实验测量结果验证了Model 1与Model 2周向残余应力计算的准确性。上述的比较表明,采用Model 2的方法,也能较高精度地模拟厚大管-管接头的周向残余应力。

图13为焊缝中心线line 5上周向残余应力的分布。图中5条曲线分别表示了第2、6、9、12和16层焊接完成时的周向残余应力分布。从此图可以看到一个有趣的周向应力的演化过程,即随着焊接进行,焊道的增加,沿着焊缝中线的周向应力的大小和分布在不断发生变化。也就是说,后焊焊道对先焊焊道产生的应力分布有显著的影响。对于下表面焊缝中心位置的周向应力而言,随着焊道数增加,其值有不断减小的趋势。产生这种变化趋势主要是因为接头本身板厚较厚,加之焊道数和焊接层数较多,周向应力分布除了在与焊缝中线垂直方向(横向方向)上产生平衡外,在厚度方向也产生平衡。

图13

图13   周向残余应力沿中心线分布(90°截面上不同焊层完成后)

Fig.13   Hoop residual stress distributions along the welding centerline in 90° section after welding at different layers (Inset shows the measuring direction)

3.2.2 径向残余应力结果

图14a和b分别为Model 1与Model 2整体径向残余应力(类似平板接头的横向残余应力)分布云图。整体上而言,Model 1与Model 2残余应力分布形态非常相似,但在Model 1的外表面焊缝始终端位置(0°)可以清晰地看到径向残余应力分布与周向残余应力呈现相似的不光滑的局部区域。而Model 2中,由于全部焊道采用瞬间热源模型及对称的热边界和力学边界条件,故径向方向的残余应力呈现出光滑连续的分布特征。

图14

图14   焊后Model 1与Model 2整体径向残余应力分布云图

Fig.14   Overall radial residual stress distributions for Model 1 (a) and Model 2 (b)

图15为Model 1在0°、90°、180°和270°截面上焊缝及其附近区域的径向应力分布云图。径向残余应力在焊缝上表面为拉应力,且呈现出离焊缝越近拉应力越高的特征。根据Dong[26]的研究,管-管对接环焊缝在轴向方向(与子母线平行的方向)应力(在本工作中为径向应力)的分布,可分为弯曲型(bending-type)和自平衡型(equilibrating-type) 2种类型。本工作虽然是马鞍型接头,但它也有类似于环焊缝的一些特征,其沿厚度方向上的径向残余应力呈现拉-压-拉的分布特征,类似于厚大尺寸环焊缝自平衡类型[10]。对0°、90°、180°和270°不同截面位置的径向应力分布比较可知,90°与270°截面的径向残余应力大小与分布均非常相似,同时0°与180°截面的径向残余应力大小与分布也非常相似；而90°与180°截面位置径向残余应力有显著的差异。因此可以推测,对于该马鞍形管-管接头而言,采用从90°至180°截面区域的1/4模型即可描述整体结构的径向残余应力分布。

图15

图15   Model 1典型截面上的径向残余应力云图

Fig.15   Radial residual stress distributions in the characteristic cross-sections for Model 1

(a) 0° (b) 90° (c) 180° (d) 270°

图16a和b分别为Model 2在90°和180°截面上焊缝及其附近区域的径向残余应力分布云图。与图15b中的90°截面和图15c中的180°截面上的径向残余应力分布比较可以看出,Model 1与Model 2的径向残余应力分布及大小基本一致。

图16

图16   Model 2典型截面上的径向残余应力云图

Fig.16   Radial residual stress distributions in 90° (a) and 180° (b) cross-sections for Model 2

图17为90°截面的上表面(line 1)和下表面(line 2)上径向残余应力分布。图17a为line 1上的径向残余应力分布,整体上看2者的分布形态基本一致,2者的峰值拉伸应力位置也很接近。实验测量的最大值为448 MPa,相同位置Model 1的值为451 MPa,Model 2的值为484 MPa；在-20 mm处实验测量的最小值为128 MPa,相同位置上Model 1与Model 2的值分别为41和25 MPa。在90°截面位置的下表面line 2上的径向残余应力分布如图17b所示,可以看到,无论是应力分布形态还是大小,数值模拟结果与实验结果吻合良好。整体上而言,实验值略高于计算值。在line 2上,实验测量最大值为393 MPa,而Model 1与Model 2对应的计算值分别为367和348 MPa,实验与计算得到的残余应力峰值和分布吻合较好。

图17

图17   90°处上(line 1)下(line 2)表面径向残余应力实验值与计算值

Fig.17   Experimental and numerical results of radial residual stress in 90° cross-section, located at line 1 (a) and line 2 (b) (Insets show the measuring directions)

图18为180°截面的上表面(line 3)和下表面(line 4)上的径向残余应力分布图。line 3上的径向应力分布如图18a所示,在-30~0 mm范围内,数值模拟结果与实验测量值之间存在较大的差异,但在0~30 mm范围内,数值模拟结果与实验值在分布形态和大小上均吻合较好。line 4上的径向残余应力分布如图18b所示,可以看到,计算结果与测量结果吻合很好,整体上而言,实验值略高于计算值。

图18

图18   180°处上(line 3)下(line 4)表面径向残余应力实验值与计算值

Fig.18   Experimental and numerical results of radial residual stress in 180° cross-section, located at line 3 (a) and line 4 (b) (Insets show the measuring directions)

图17和18中径向残余应力的实验结果与数值模拟结果尽管存在一定程度的差异,但是从应力分布形态和峰值应力大小上看,2者吻合较好,实验测量结果验证了Model 1与Model 2径向残余应力计算的准确性。上述比较表明,即使采用Model 2也能较准确地获得管-管接头的径向残余应力大小与分布。

图19为焊缝中心线line 5上径向残余应力的分布。图中5条曲线分别表示了第2、6、9、12和16层焊接完成时的径向残余应力分布。随着焊接进行,焊层数的增加,沿着焊缝中线的径向应力的大小和分布不断发生变化,这说明后焊焊道对先焊焊道产生的径向应力分布有显著的影响。随着焊道数的增加,下表面焊缝中心位置的径向应力呈现出先减小随后又增加的趋势,沿着焊缝中线的径向应力在焊接过程中的演化过程与周向应力的演化过程有明显的不同,前者先减小后增大,而后者随着焊道增加又有一直减小的趋势。出现这种现象是因为径向应力分布除了在与焊缝中线垂直方向(横向方向)上产生平衡外,在厚度方向也产生平衡。

图19

图19   径向残余应力沿中心线分布(90°截面上不同焊层完成后)

Fig.19   Radial residual stress distributions along the welding centerline in 90° section after welding different layers (Inset shows the measuring direction)

3.3 计算时间与储存空间

下面对2种方法的计算效率进行定量比较。2个模型的计算均在i9-9900K (CPU参数为8核16线程,主频为4.7 GHz) 的计算机上完成,采用移动热源模型计算第1道焊的应力时所花费的时间为9.65 h,第2道焊的时间花费9.68 h,第3道焊用时为9.69 h。从数值上看,3个焊道计算时间的差异较小,平均每个焊道的计算时间为9.67 h。根据前3道的所用计算时间可以得到这样一个推测,全模型+全移动热源模型,51道全部采用移动热源模型来计算焊接应力,所需计算时间约为493 h。当采用Model 1来计算焊接应力时,实际花费时间为137 h,而采用Model 2,计算时间仅仅为8.4 h。单纯从计算时间上看,采用全模型及全焊道移动热源模型花费的时间最长；采用Model 1的计算时间居中,所用时间仅为前者的27.7%；而采用Model 2的计算时间最少,以全模型+全移动热源模型为参照,所用的计算时间仅为它的1.9%。

从所需的存储空间而言,Model 1与Model 2所需的储存空间分别为1604和191 GB。比较而言,采取Model 2可大大节约储存空间。基于计算时间和储存空间的比较,Model 2具备更强的工程应用潜力。

3.4 讨论

本工作采用Model 1和Model 2计算了厚大马鞍形管-管焊接接头的温度场和残余应力,并与实验结果进行了比较。不论是典型位置的焊接热循环还是焊缝及其附件区域的残余应力,数值模拟结果与实验测量值总体上吻合良好。

3.4.1 计算误差

数值模拟得到的残余应力与实验测量结果尽管在峰值应力和分布形态上基本一致,但依然存在一定程度的差异。通过分析,主要的原因如下：

(1) 初始材料状态难考虑。在实际制备母管和台管时,材料经过固溶处理、矫正变形及机械加工等焊前的工序,这些工序可能会产生不同程度的塑性变形以及初期残余应力,在有限元模型中这些焊前材料初始状态没有被考虑,从而产生误差。

(2) 实验测量残余应力带来的误差。实验测量采用应力释放法进行残余应力测量时,需要在测量位置进行应变片的粘贴。在应变片粘贴前,需要将测量位置打磨平整,由于在180°截面的位置存在一定弧度,故而需要在基体上打磨掉的材料较多,才能达到平整的效果。在这个过程中,打磨较多会对焊缝位置引入额外的塑性应变,使测量结果出现误差。从位置看,在180°处因打磨造成的误差可能会更大。

(3) 有限元计算模型误差。在材料模型中,本工作采用“各向同性硬化模型”和“阶跃式退火模型”来分别考虑加工硬化和退火软化效应,而实际的材料会出现包辛格(Bauschinger)效应,同时软化效应与时间和温度都有密切关系。在本工作采用的材料模型中,上述的问题没有完整地考虑,在今后的研究中,需要开发更高精度的材料模型。

3.4.2 应力峰值

实验测量结果和计算结果均显示,在焊接后接头中存在大于550 MPa的高拉应力区,而SUS316不锈钢的常温屈服强度仅为220 MPa,应力峰值远远超过了材料的常温屈服强度。出现如此高的拉伸残余应力的主要原因如下：

(1) SUS316奥氏体不锈钢的加工硬化效应显著。在焊接过程中,焊缝金属及其附近的母材都要经历峰值温度较高的热循环,由于焊缝及其附近的材料在升温和降温过程中产生的热应变要受到周围温度相对较低金属的拘束作用而发生塑性变形,塑性变形的产生致使材料发生加工硬化而提高流变应力。

(2) 由于本工作中的焊接接头板厚较厚,而且形状相对复杂,因此焊缝处的三向内拘束较大。在内拘束大的情况下,焊接过程中容易产生多轴(二轴或三轴)拉伸应力,这样会导致材料不易屈服。

实验结果表明,金属塑性屈服准则中最接近实际条件的是Von Mises屈服准则[27]。在本工作的数值模拟中材料模型也是采用的Von Mises屈服准则。Von Mises屈服准则中规定金属屈服条件是在第一主应力(σ1)、第二主应力(σ2)和第三主应力(σ3)共同作用下,材料会发生屈服,其遵循 式(5)[27]：

式中,σs为屈服强度。

式(5)表示的屈服面如图20a所示, 式(5)左侧的值小于σs时,材料将不发生塑性变形,处于弹性状态。从纯数学的角度考察 式(5)可知,若当满足σ1 = σ2 = σ3条件时,无论σ1、σ2和σ3多大,材料仍不发生屈服,所以空间中3个应力共同作用时存在某个方向的应力值会大于σs的情况。

图20

图20   SUS316不锈钢空间Von Mises屈服面与平面Von Mises屈服面及其扩大示意图

Fig.20   Schematics of SUS316 Von Mises yield surface in the triaxial stress state (a) and the Von Mises yield surface expansion in the biaxial stress state (b) (σ1—major principal stress, σ2—second principal stress, σ3—third principal stress, σs—yield strength, σ${}_{\mathrm{s}}^{\text{'}}$—yield strength after work hardening)

若只考虑2个应力,可假设σ3 = 0时, 式(5)可简化为一个椭圆方程式：

式(6)椭圆方程的几何形状如图20b所示。此时材料的屈服面则为椭圆,由于采用各向同性硬化准则,当材料发生屈服后,随着塑性应变的增加,屈服面将增大,此时材料将发生应变强化即屈服强度增大。从图21中90°截面Von Mises应力分布云图可知,在焊缝上半部及其附近区域的Von Mises应力明显高于材料常温的屈服强度(220 MPa),其峰值应力达到550 MPa。高应力的产生,既有多轴应力的效应,也有加工硬化的贡献。

图21

图21   Model 1模型90°截面等效Von Mises应力

Fig.21   Equivalent Von Mises stress distributions in 90° cross-section for Model 1

采用的加工硬化及退火软化模型计算得到残余应力,与实验得出的残余应力接近,验证了材料模型准确性。从计算结果和实验结果来看,由于对模型热源的简化和对材料模型的近似处理,使该模型占用的计算资源较少,同时计算时间也大幅缩短。基于应力分布形态及峰值应力实验结果有较高吻合度的事实,从面向工程应用的角度出发,本工作开发的2种计算手法具有较高的工程应用价值。

4 结论

(1) 基于通用有限元软件,针对结构较大、焊道布置复杂的SUS316马鞍形管-管接头温度场和残余应力的计算,开发了全尺寸有限元网格模型+“移动热源+瞬间热源”混合计算模型(Model 1)以及1/4部分网格有限元模型+全焊道瞬间热源模型(Model 2)。计算结果与实验结果吻合较好,验证了2种模型的可行性和妥当性。

(2) Model 1能计算出整体结构焊接时的残余应力,并且也能抓住始终端位置的应力分布特征；Model 2虽然不能模拟始终端位置的应力,但可以较准确地获得除始终端位置及其附近以外区域的残余应力,而且计算时间可大幅缩短,存储空间也能大幅减少。在实际的工程应用中,对于厚大焊接接头或结构的焊接残余应力,若不重点考虑始终端位置的应力特征,可以利用结构的对称性简化模型,同时采用瞬间热源模型来模拟焊接热输入,这种方法可以兼顾计算精度并能大幅提高计算效率。

(3) 从数值模拟结果和实验测量结果来看,焊缝及其附近的残余应力呈现出多轴拉伸的特征,局部位置的单向(周向)应力达到了600 MPa,远远高于材料常温时的屈服极限。产生如此高的残余应力的原因主要有两点,一是SUS316材料有较显著的加工硬化倾向,另一个是在厚大接头焊接时内拘束强,而引起的多轴应力效应。