考虑固相变的钛合金焊接残余应力数值仿真方法

郦 羽，万正权，唐文勇

（1.中国船舶科学研究中心，江苏无锡 214082；2.上海交通大学船舶海洋与建筑工程学院，上海 200011）

0 引 言

焊接温度场、残余应力和变形的实验测量既耗时又昂贵，而且还可能存在与测量相关的不确定性[1]。因此，焊接计算模拟经常被用来预测不同类型材料和接头在焊接过程中和焊接后的瞬态残余应力场和变形的分布和大小。由于焊接过程中的温度梯度非常大，焊缝在小体积内是由复杂的非均匀组织组成。在冷却过程中焊缝可能会发生固态相变，相关的体积变化和相的具体材料特性会对焊接残余应力和变形产生显著影响。许多焊接模型忽略了相变的影响，Ferro等[2]认为这种影响是不可忽略的。

少数研究人员研究了Ti-6Al-4V 的相变动力学，其中大多数研究是针对Ti-6Al-4V 焊接的具体的案例，忽略了复杂相变过程影响，因为Ti-6Al-4V 随化学成分、温度历史和相变前的初始相形态而变化。Malinov 等[3]利用电阻率技术研究了等温条件下的相变动力学，并利用差示扫描热法研究了不同冷却速率下、连续冷却条件下的相变动力学，利用JMAK 方程[4]和叠加性规则确定了β相转化为α相的分数；Ahmed 和Rack[5]研究了在β相转变点以上的冷却速度对冷却过程中相变的影响，并确定在高于410 ℃/s 的高冷却速度下形成马氏体；而Gil Mur 研究了回火热处理过程中马氏体分解为α和β的过程，并在更低的冷却速率下观察到完全马氏体组织；Kelly Kampe 等[6]模拟了Ti-6Al-4V 在激光金属沉积过程中的微观结构演变；Fan 使用Malinov 等[3]的JMAK 等温转变动力学参数，数值研究了相变对Ti-6Al-4V 激光成形的影响；Elmer 等[7]利用同步x 射线衍射技术实验测量了Ti-6Al-4V 气体钨极弧焊（GTAW）过程中的转化动力学；Longuet 等[8]开发了Ti-6Al-4V的通用多相模型，并应用于直接激光制造和激光焊接工艺。

本文分别采用JMAK 方程和K-M 方程表达钛合金固态相变过程中的扩散转变和非扩散转变建立Ti-6Al-4V 焊接过程中固态相变的数学模型，确定数学模型中的参数求解方法，并基于ABAQUS 软件开发考虑固态相变的热-冶金-力学耦合的有限元计算方法，该方法可以考虑温度与组织变化对焊接残余应力的影响。通过一个对接接头的激光焊接的有限元算例，验证所建立的固态相变的Ti-6Al-4V 焊接有限元计算模型的可行性。该工作对于定量分析Ti-6Al-4V 焊接残余应力的影响因素具有参考意义。

1 焊接残余应力数值计算方法

焊接是一个涉及许多复杂物理现象的过程，焊接过程中焊件的冶金、热和力学变化之间存在多场相互作用，包括：焊接热循环引起的热应力以及伴随的熔化、凝固或固态相变，局部膨胀引起的转换应力，非弹性变形产生的热，由应力或应变加速的相变、相变潜热。对于解决工程实际问题的计算方法而言，为了兼顾计算效率与计算精度，在有限元模型中往往重点考虑主要因素，有时要忽略次要因素。

1.1 温度场的计算

由于焊接过程的物理复杂性，采用一种简化方法仅基于热方程的求解，其中热源是模拟的，而不是规定焊接模型的体积热流输入。由傅里叶定律和能量守恒推导出的热方程如式（1）所示，用于求解瞬态温度场（T）在时间（t）和空间（x，y，z）上的温度场：、，

式中，T是温度，ρ为密度，Cp为比热，k为导热系数，Q̇V为体积热通量（W/m3）。为了求解热方程，必须指定热导率、密度和比热容。对于经历相变的材料，相变潜热和微观结构演变也必须考虑。

采用三维锥形高斯热源模型对激光焊接过程进行热分析。图1 所示为三维锥形模型，其径向功率密度分布为高斯分布，轴向为线性分布，可以得到更精确的深熔激光束焊接结果[9]。圆锥热源的上表面热流最大，下表面热流最小，如式（2）所示。

图1 对接焊缝采用的三维高斯型热源模型Fig.1 Three dimensional Gaussian heat source model for butt weld

式中，r0为特定高度z处热源的半径，q0为最大体积功率密度，r为当前内部点的半径，

分布参数r0从圆锥区域的顶部到底部线性减小，如式（4）所示。

式中，ze是锥体区域的上表面，zi是下表面。通过对物体上的体积热流Q进行积分，并将其重新作为总热量输入，可以得到最大的体积功率密度，如式（5）所示。

式中，H=ze-zi，h=z-zi，e 为自然对数的底数，η为热源效率。热源效率或能量传递效率定义为被辐照材料吸收的能量与激光输出功率的比值。

焊接过程中，工件与周围空气之间通过对流和辐射的方式进行热交换。式（6）和式（7）分别将表面对流和辐射造成的热损失定义为边界条件。

式中，T为当前温度，T0为环境温度，Tabs为绝对零度温度，ε是发射率或发射热辐射的能力，而σ为Stefan-Boltzmann 常数（5.68×10-8J/K4⋅m2⋅s）。辐射的影响在较高的温度下（靠近焊缝）更大，而在较低的温度下（远离焊缝）不显著。Yang等[10]指出：Ti-6Al-4V 的发射率在760 ℃以下的温度下几乎不受温度的影响，其发射率在0.3 以下；而在760 ℃以上，它开始氧化，其发射率随着温度的升高而增加；在1000 ℃以上，其发射率几乎达到1.0。相比之下，对流效应在温度较低时更显著，而在温度较高时则不显著。用于传热系数和辐射常数的校准值为hconv=10，ε=0.4。

1.2 组织计算

钛合金在中高温区的相变可以被看成是两个均析转变（一个发生在加热过程，另一个发生在冷却过程），通过过往研究，钛合金α→β和β→α相变被认为是两个互逆的扩散型相变，扩散转变可以用如公式（8）所示的Johnson-Mehl-Avrami方程（简称JMAK方程）[4]来描述。

式中，f为时刻t时产物相的体积分数，feq是产物相T温度下的平衡体积分数，与温度相关的k( )T和与温度无关的n是定义相转变动力学的JMAK方程参数。

为了将等温模型的应用扩展到焊接过程中的非等温相变，使用了Scheil加法规则[11]。Scheil规则将瞬时小的等温时间增量步长叠加的总和近似表示任意连续的温度变化。通过设置一个较小的计算时间步，就可以采用阶梯状的等温变化曲线来近似焊接过程的连续变化曲线。

对于每一个计算时间步，引入虚拟时间t*0，虚拟时间的表达式为

式中，t*0是在温度T1的等温转变过程中达到f(t0,T0)所需的虚拟时间，feq(T1)是产物相在温度T1下的相平衡分数。

在快速冷却过程中，β相转变为马氏体相的过程为非扩散型转变，采用Koisten-Marburger 方程（简称K-M方程）进行描述。

式中，fM是马氏体相的体积分数，kM是依赖材料的参数，MS是马氏体开始温度。

1.3 应力场计算

计算焊接应力场时，将温度场的计算结果以热载荷的方式加载到应力计算模型中，并将组织（相）的计算结果代入材料的参数模型中，采用热-弹塑性有限元方法进行力学分析。假设材料的弹塑性行为具有各向同性硬化规律，钛合金焊接过程中的总应变包括以下几部分：

式中，εTotal是总应变，εe是弹性应变，εp是塑性应变，εth是热应变，εc是蠕变应变，εvp是粘塑性应变，εtp是相变塑性应变。总应变由上述应变分量组成。假定蠕变、粘塑性和突变分量引起的应变可忽略不计，因此不包括在计算中。由于忽略了粘塑性效应，假定屈服应力与应变率无关，且取决于塑性应变和温度。根据Lindgren 的研究[12]，焊接模拟必须考虑弹性、塑性和热应变来计算残余应力。热应变对焊接诱导残余应力的影响最大，因此，热膨胀系数是应力计算中最重要的材料参数之一。Lindgren[13]进一步指出，焊接冷却阶段是导致残余应力和变形的重要阶段。在材料经历相变的情况下，所产生的体积变化可以包括在热应变分量中。

2 Ti-6Al-4V固态相变模型的建立

2.1 加热过程相变数学模型建立

在加热阶段，Ti-6Al-4V 经历了从α相到β相的快速扩散控制转变，其中包括β的形核和生长以及α/β界面的运动导致的α的溶解，这是由β稳定剂在界面上的传输引起的。β的体积分数在β转变温度（Tβ）以上达到100%，并在熔化温度下保持稳定。假设β生长增加量是通过相等量的α溶解，该过程为扩散转变，采用JMAK 方程来描述。为了避免球状、魏氏组织、网篮组织和晶界α相等相形态的差异，假设了一个简化的单α相，如公式（13）所示，计算加热过程中β相的体积分数。

式中，t*0是在温度T1的等温转变过程中达到fβ(t0,T0)所需的虚拟时间，feqβ是β相的平衡分数。

阅读不是让学生记住多少知识，而是在阅读的过程中最终学会阅读。因此，家长在阅读绘本的过程中可以渗透相关绘本阅读的策略，比如，看图获取信息并进行推断，对后来发生的事进行预测，阅读中的自我提问，结合学生的实际等。家长们可以根据自己的优势和特长自主设计大胆尝试，让孩子们充分体验阅读带给他们的快乐，进一步提升他们的阅读能力，培养他们热爱阅读的浓厚兴趣和良好习惯，从而让孩子们逐渐学会阅读绘本，构建阅读力。

2.2 冷却过程相变数学模型建立

在冷却阶段，冷却速率决定了转变动力学和形成的相。在低于β相变点的温度下，当冷却速率小于20 K/s时，β通过成核和扩散反应分解成α[5]。采用JMAK方程描述。

式中，t*0是在温度T1的等温转变过程中达到fα(t0,T0)所需的虚拟时间，feq

α是α相的平衡分数。

在β相变点以上温度快速冷却，将β转变成另外一种形式的α相，即细针状马氏体相α'，其不同于平衡α相，由于相邻片层和界面的结晶取向差异较大，且重复率较低，因此具有较高残余应力[14]。

根据Ahmed和Rack的研究[5]，对于大于410 ℃/s的快速冷却速率，将形成完全马氏体组织，并抑制向α相的扩散转变。

对于20 ℃/s到410 ℃/s之间的中等冷却速率，从β相变点以上温度开始冷却，将导致β向晶界的块状α'部分转变，和靠近原β晶界的晶体内马氏体αm的部分转变。将α'纳入αm中，考虑到缺乏有关该问题的精确数据，可进行另一种选择，并从可用于转变的β相的量中减去β相的当前平衡分数feqβ，以防止完全转变为马氏体α'。

式中，fαm是αm的分数，fβeq是β的平衡分数。在这种情况下，从可用于转化的β相中减去β的当前平衡分数，以防止完全转化为α'。

2.3 相变数学模型中的参数确定

为了求解JMAK 方程的产物相，需要确定方程中的参数k( )T和n，以及产物相的平衡相，图2 显示了JMatPro 模型确定的α相和β相的平衡分数。由于实际的母材组织中存在β相，在低温条件下，β相体积分数预测不能为0。

图2 α相和β相平衡相分数随温度变化的模型Fig.2 Equilibrium fraction of α and β phases varying with temperature

本文基于Szkliniarz 和Smotka[15]的试验数据，在不同加热和冷却速度下材料处于不同温度时的α、β两相体积分数，对其两组试验数据，对方程（8）进行改写，采用最小二乘法拟合求解。

通过拟合的数据发现，n0基本不随温度变化，k0（Ti）随温度而变化。为了使JMAK 方程计算的值更准确，将拟合的k0（Ti）和n0作为初值，将该温度下的JMAK 方程计算值fJMAK(t,Ti)与该温度下的试验测得值fmeasure(t,Ti)之差的平方（fJMAK（t,Ti）-fmeasure（t,Ti））2，与试验测得值fmeasure(t,Ti)的比值记为误差值Eerr，通过搜索求解k(Ti)和n的值，使得Eerr值最小而作为该温度下的一组最优解（k（Ti）,n）。

为了较为准确地表达出k（T）随温度的变化，Elmer 等将随温度变化的参数k（T）表示为式（21）所示的形式，p和Q均为不随温度变化的系数，R为气体常数。将以上优化的一组不同温度下的k（T）值，通过公式（21）拟合求得p、Q值。

对于K-M 方程中的参数求解，通过GLEEBLE 热模拟试验机对试样进行热膨胀试验得到热膨胀曲线。在较大冷却速度范围内冷却组织转变为马氏体相，通过双切线法确定马氏体相变的开始温度MS。由K-M方程的逆运算得到马氏体相变参数kM的值，计算公式如下：

图3 显示了β的平衡相分数以及α'的平衡相分数，它们使用K-M 方程的kM或γ和MS函数确定。选择γ等于0.015，而MS等于650°C。

图3 用Koistenen-Marburger模型计算β和α'的相分数，作为Ms和γ的函数Fig.3 Calculation of phase fractions of β and α'with Koistenen-Marburger model as the functions of Ms and γ

3 Ti-6Al-4V的结果与讨论

3.1 焊接过程有限元模型的建立

为了验证所建立的钛合金焊接过程中的固态相变过程，本文采用一个50 mm×50 mm×1.5 mm的板进行对接焊计算，激光焊接速度为50 mm/s。计算模型的网格尺寸为0.5×0.5×0.5 mm3，单元数量为60 000。在热分析中，采用2 阶三维20 节点二次扩散传热六面体单元（DC3D20）。在力学分析中，采用了一阶三维8结点线性六面体简化积分单元（C3D8R）。力学分析中位移的有限元形状函数通常要比热分析高一个维度。简化积分单元，C3D8R 使用低阶积分来形成单元刚度，并提供更好的近似真实行为。Ti-6Al-4V随温度变化的热物理材料性能参数如图4所示。

图4 Ti-6Al-4V随温度变化的热物理材料特性Fig.4 Thermo-physical material properties of Ti-6Al-4V varying with temperature

3.2 焊接过程中相的成分变化

由第2 章所建立的Ti6Al4V 相变模型可知，Ti6Al4V 的相变转变类型主要与焊接过程中的三个变量相关：当前温度、冷却/加热速率和峰值温度。通过对焊接过程中这三个变量的判断，可以根据相变方程计算出相应的相产物体积分数。为了实现焊接过程的相变模拟，在热分析过程中，本文使用ABAQUS 的USDFLD 子程序的内部状态变量作为冷却/加热速率和峰值温度的函数，计算不同微观结构相的相变体积分数。在力学分析中，通过UEXPAN 子程序调用USDFLD 子程序的相体积分数场变量，将不同相的热膨胀特性参数与相变体积分数叠加运算，从而考虑焊接相变对热应变影响，最终得到残余应力的影响。

在模型计算过程中，当前温度为场变量NT11（℃），冷却/加热速率定义为状态变量SDV3（℃/s），峰值温度定义为状态变量SDV4（℃），场变量FV1 定义为β相，场变量FV2 定义为α相，场变量FV3 定义为马氏体相αm。图5～7 为焊接过程中，0.25 s、1.0 s、4.0 s 三个时刻仿真计算场变量变化；图8～10 为焊接过程中，0.25 s、1.0 s、4.0 s 三个时刻钛合金板的组成成分体积分数变化。Ti6Al4V 试样板初始组成成分：β相体积分数为0.075 26，α相体积分数为0.924 74，马氏体αm相体积分数为0。

图5 焊接过程中0.25 s时刻仿真计算场变量变化Fig.5 Simulation calculation of field variable change at 0.25 s during welding

图8 焊接过程中0.25 s时刻组织成分体积分数变化Fig.8 Volume fraction changes of different micro-phases at 0.25 s during welding

由峰值温度SDV4 可知，焊接过程中，只在焊缝及热影响区域的一个较小范围内峰值温度超过650 ℃，也即焊接过程中的相变主要发生在该区域。由当前温度NT11 可知，随着热源的移动，温度场由熔池（温度超过1650 ℃的区域）向外逐渐降低，呈椭圆分布。温度场的分布特征决定了冷却/加热速率SDV3 的分布，由图5 和图6 可知，随着热源移动，由内到外呈现冷却速率到加热速率的过渡，冷却速率由内到外逐渐减小，最内部超过-410 ℃/s 的区域呈长椭圆分布，逐渐过渡到-410～-20 ℃/s 的区域。

图6 焊接过程中1.0 s时刻仿真计算场变量变化Fig.6 Simulation calculation of field variable change at 1.0 s during welding

对图5进行分析，在0.25 s时刻，焊缝及热影响区域当前温度都在980 ℃以上，所以都转化为β相，与图8 计算的相变组成相吻合。在1.0 s时刻，整个焊缝及热影响区域的温度在650 ℃以上，所以马氏体αm相不会产生，冷却速率在-410～-20 ℃/s 区域时，生成α相，与图9 的相变组成相吻合。在4.0 s 时刻，也即冷却过程，焊缝及热影响区域周围有相当大部分区域的温度已经降到300～650 ℃之间，焊缝及热影响区域的冷却速率在-410～-20 ℃/s，因此该阶段焊缝及热影响区会产生马氏体αm相，由图10可以发现此时焊缝及热影响处主要是马氏体αm相和α相，含有少量β相，与图7也吻合得较好。因此，仿真计算过程基本实现了第2章中Ti6Al4V固态相变的相变模型过程。

图7 焊接过程中4.0 s时刻仿真计算场变量变化Fig.7 Simulation calculation of field variable change at 4.0 s during welding

图9 焊接过程中1.0 s时刻组织成分体积分数变化Fig.9 Volume fraction changes of different micro-phases at 1.0 s during welding

图10 焊接过程中4.0 s时刻组织成分体积分数变化Fig.10 Volume fraction changes of different micro-phases at 4.0 s during welding

图11 显示了焊接热循环过程中，焊缝中心距起始位置20 mm 处积分点位置α相、β相和α'相体积含量变化的数值模拟结果。初始组织由约0.924 74α和0.075 26β组成。在加热过程中，随着温度的迅速升高，α 相分数迅速下降到0，而β相增加到1，实际上，在熔融温度以上也应该变成0，但为了简单起见，在β转变点温度以上，它仍然维持在1。当材料开始冷却到β相变点温度以下时，根据冷却速率的不同，一个不同的微观结构由β相的扩散或非扩散转化产生。由图11 可知，冷却过程中β相先转变为α相，冷却到一定温度后，α相的体积分数保持不变，剩下的β相转变为马氏体αm相，达到β相在该温度下的平衡相后，各相成分保持不变。在冷却至室温过程中，最终相主要由0.472 7α，0.488 64αm以及0.038 66β相组成。

图11 焊接中心处α、β和α'相分数随温度变化的时间-温度历史Fig.11 α,β and α'phase fractions at the weld center line chang⁃ing with temperature in the time-temperature history

3.3 焊后残余应力和平面变形

从图12 和图13 可以看出，在数值焊接模型中引入相变效应会对焊接过程结束时焊缝周围残余应力的最终分布产生一定改变。相变模型的y向应力分量S22 峰值比无相变模型大100 MPa，且焊缝宽度上的应力梯度更陡，残余应力范围更大。同样，通过相变计算得到的x向应力S11 和z向应力S33 值也会增大，且改变了拉压应力的分布情况，但只是增大了很小的一部分。具有相变效应的变形结果比没有相变的变形结果略有增加，弯曲变形和角变形的幅度在0.1 mm以内，由于模型尺寸较小，变形不是很明显。

图12 中横剖面处三个主应力方向的残余应力分布Fig.12 Residual stress distribution in three principal stress directions at the middle cross section

图13 工件焊后残余应力和变形的分布Fig.13 Distribution of residual stress and deformation of workpiece after welding

基于以上仿真结果分析，由于α和β单个晶体的尺寸匹配良好，这意味着在新相的形核和溶解过程中，没有任何原因改变初始晶向，但是钛合金材料在不同相之间的体积比存在较大差异。因此，Ti-6Al-4V 中的相变会引起焊缝体积变化，从而导致焊接残余应力的一定变化，但是不会显著改变残余应力的分布情况。

4 结 语

本文建立了Ti-6Al-4V 焊接过程中固态相变的数学模型，确定了数学模型中的参数求解方法，并基于ABAQUS 软件开发了考虑固态相变的热-冶金-力学耦合的有限元计算方法。基于该计算模型，对Ti-6Al-4V 对接板激光焊接进行了计算验证，得出了焊接接头的组织包含α相、马氏体αm相和少量β相。并基于以上结果，进行了焊接残余应力计算。通过对计算结果分析可以发现，考虑相变后残余应力y向分量S22增加明显，x向分量S11和z向分量S33的分布规律有较明显改变，但是幅值较小。钛合金材料在不同相之间的体积比存在较大差异，因此相变会引起焊接残余应力的一定变化，但是不会显著改变残余应力的分布情况。