旋转摩擦焊(rotary friction welding, RFW)是一种先进的固态焊接技术，被广泛应用于汽车和航空航天领域，用于制造如发动机阀门和驱动轴等关键部件[1-2]。该工艺通过回体件金属在高压力(通常为10～200 MPa)和高速旋转(200～5 000 mm/s)的条件下相互摩擦来实现连接。在摩擦界面附近，金属材料因受热和剧烈的塑性变形而软化。随着摩擦过程的推进，系统逐渐从初始的加热阶段过渡到稳定状态。在此阶段结束后，通过施加或不施加顶锻压力的制动，最终形成牢固的焊接接头。
在焊接过程温度达到准稳态温度时，界面金属处于黏塑性流动状态，黏滞耗散随之成为重要热源。随着摩擦焊热源研究的深入，研究表明，摩擦焊界面的热源最终演变为与应变速率相关的黏塑性金属层。作为热源本身，该层不仅产热，还通过热传导调节温度场。但考虑到焊接过程涉及复杂的热力耦合，因此，分别描述摩擦焊工艺的热分析或力分析是方便和常规的。
第1 个热分析模型由Rykalin 等[3]于1959 年提出，假设产热恒定并完全由界面摩擦产生，其核心在于摩擦系数的确定，此模型被广泛应用于计算各种摩擦焊接工艺，如旋转摩擦焊[4-5]，线性摩擦焊[6-8]和搅拌摩擦焊[9-11]。然而，摩擦焊接过程中相变和组织演化会改变材料表面性质，从而增加摩擦系数测量的难度[12]。为解决传统模型对摩擦系数和剪切屈服应力等参数的依赖，Maalekian[13]采用逆向热传导(inverse heat conduction,IHC)方法，通过实验数据反推出界面热流量，降低了对材料物性参数的要求，并显著提升了温度场预测的精度。Yang 等[14]进一步引入共轭梯度法和误差原理，利用温度测量值直接反推出界面生热分布，无需预设摩擦系数或热生成函数。研究表明，温度依赖性热导率和比热容的引入显著提高了预测精度，但测量点位置对结果影响较大，越接近界面，逆问题求解的准确性越高。然而，这类热模型的主要局限性在于未考虑塑性流动，而塑性流动是RFW 稳态阶段的核心生热机制，对接头焊合质量起着至关重要的作用。
第2 种方法是引入焊接接头的响应参数，以减少未知参数的数量，因为耦合机制的缺失不可避免地导致欠定方程组。这种方法通常用于基于动力学理论开发的分析模型，因为耦合机制本质上是热力学标准。引入响应参数，如温度分布形式[15]或轴向缩短率[16]，使得求解欠定和动力学方程系统成为可能，但是，由于预先假定的响应，这些动力学模型实际上是回顾性分析。正如Bhamji 等[17]指出的那样，他们无法根据给定的初始焊接参数集进行前瞻性预测分析。摩擦焊接形成的大变形体系是远离平衡态的耗散系统，其热力耦合机理可从Onsager-Ziegler 最大熵产生原理(Onsager-Ziegler maximum entropy production principle, OZ-MEPP) 推导[18]。该理论指出，远离平衡的开放系统的熵产生(单位时间和单位体积的不可逆动力学行为引起的熵增加)具有最大值，并受到其边界条件的限制。这一最大值原理为建立动力学行为耦合机制提供了热力学判据。最近，基于OZ-MEPP，本课题组提出了二维分析模型来描述线性摩擦焊(linear friction welding, LFW)产生的大变形体系的稳态[19]。LFW 过程是由两个矩形试件之间的相对往复运动完成的，其物质流和熵流(熵流是系统与环境之间的能量和质量交换所引起的熵变，相当于稳态时的负熵产)具有二维特征，因此二维模型适合描述LFW。同样的，在旋转摩擦焊中将接头划分为： 高速率剪切变形区(high-rate shear deformation zone,HSDZ)与低速率压缩扭转变形区(low-rate compression-torsion deformation zone, LCDZ)，并假设HSDZ 区域产热主要是塑性变形产热，体系熵产生主要集中于此区[20]。LCDZ 区的剪切应力与应变平行于飞边生长方向。得到了二维模型下两个变形区域的本构方程，随后借鉴实际焊接过程中摩擦焊的飞边是沿着整个界面三维方向生长的特点，结合二维模型大于实际三维状态下飞边拘束度，将两变形区本构模型推广至三维，运用热力学原理，建立了摩擦焊稳态阶段中工艺参数与最终接头响应参量之间的联系。这些方程定量地描述了接头不同区域的应力应变关系，而且在无需引入未知参量的情况下保持了模型的自洽性。但上述摩擦焊热力学模型仅实现了对同种摩擦焊预测性模型建立，对异种摩擦焊相关研究还是较少。
本研究旨在建立异种旋转摩擦焊稳态大变形体系热力耦合模型。通过最大熵产生、热传导和质量守恒理论，并结合摩擦焊大变形区金属流动方式和塑性变形本构方程，推导出异种摩擦焊大变形区稳态阶段解析模型；随后对异种界面焊合过程进行描述，建立异种接头界面焊合判据，并结合上述解析模型，最终建立异种旋转摩擦焊稳态阶段热力耦合解析模型。选用TC4/TC17 异种钛合金作为实验材料，在不同焊接参数下进行焊接，总结焊接过程中稳态阶段黏塑性区宽度、轴向缩短速率和界面温度变化规律，和TC4/SUS321 异种材料稳态阶段黏塑性区平均温度变化规律，验证模型的可靠性。
1 实验方法
1.1 旋转摩擦焊稳态描述
图1为旋转摩擦焊稳态阶段的接头结构示意图，其中x y z 坐标系分别为飞边挤出方向(y 向)、摩擦方向(z 向)和单位长度内焊接界面一侧沿焊接(或摩擦)界面法线方向(x 向)。图中接头可大致分为3个区域，即弹性区Z1、沿界面法向分布的大变形区Z2(线AB为大变形区与弹性区的分界线)以及包覆大变形区的飞边区Z3。其中大变形区(Z2 区)以塑性剪切变形为主，因此，大变形区为主要产热区域。在焊接压力P 的作用下，Z2 区内金属由侧向(垂直于OB)沿y 方向挤出进入Z3 区形成飞边，同时Z1 区的金属从垂直于AB 方向进入Z2 区。Z2 区的黏塑性区宽度OB、Z1 区金属向Z2 区流入速度(即单边轴向缩短速度)-νx 和界面处温度T 3 个参量是随着焊接参数变化的响应参量。在焊接压力作用下，Z2 区受压往Z3 区挤出的过程与Z1 区受压往Z2 区推进的过程均是处于与时间无关的稳态，因此Z2 区尺寸、Z1区金属向Z2 区流入速度-νx、Z2 区金属流出速度(即飞边生长速度)νy、温度和应变速率分布等各动力学参数均保持恒定，以下仅为一侧接头模型，另一侧共享相同模型推导过程，仅存在参数上的变化。
图1 旋转摩擦焊稳态的接头形貌特征
Fig.1 Morphological characteristics of the joint welding in the steady state of RFW
针对建立大变形区较为复杂的热力学解析模型，需要首先对模型进行简化分析，忽略飞边对大变形Z2 区中的热传导与热对流的影响，先构建相对简单的一维解析模型，再考虑飞边对大变形体系的影响，并建立稳态阶段二维解析模型
1.2 异种接头稳态一维OZ-MEPP 解析模型
图2为当摩擦焊过程已经达到稳定状态时垂直于摩擦方向-z 的横截面的示意图，并基于以下假设开发了异种接头一维分析模型。
图2 异种旋转摩擦焊接头的大变形系统一维物理模型示意图：(a)一侧接头；(b)另一侧接头
Fig.2 Schematic diagram of a one-dimensional physical model of a large deformation system for a dissimilar rotary friction welding joint:(a)one side of the joint;(b)the other side of the joint
(1)大变形体系区域为A1B1OC1(A2B2OC2)，由于在实际接头中，该体系被飞边和弹性区包围，因此忽略了大变形区的热辐射。
(2)忽略了飞边方向对热传导和对流的影响。从大变形体系中挤出的金属以-νx1 的速度沿-x 方向流过焊接界面，因此等温线平行于焊接界面。
(3)只存在平行于摩擦方向z 的剪切应力和剪切应变速率(σz 和
)。
(4)塑性变形和产热都集中在大变形体系中(黏塑性区、变形功和产热全部集中在大变形区，而弹性区黏塑性变形及其变形功和产热忽略不计)。
(5)沿焊接界面法向的试样长度为无限长，忽略Z1 对大变形体系的影响，大变形体系在热力学上处于远离平衡的非线性非平衡稳态，参数如焊接压力(P)或摩擦速度(U)、大变形体系在z 方向上的高度(OB1)、界面温度(T)和摩擦方向z 上的剪应力(σz)不随时间变化，保持恒定。
(6)被焊试样的热物理参数，如热导率(λ)、热扩散率(κ)和密度为恒定，与温度无关。
根据傅里叶热传导方程和能量守恒定律，可得到大变形体系中单位体积的热积累φ，且体系处于稳态，体系内任何一确定的空间位置φ 值为0，因此体系中的热积累可通过下式计算：
式中，q为体积产热功率；λ1为单侧的热导率(为常数)；T″x为T 在轴向x 的二阶导数。并根据假设(3)和(4)可得到方程：
式中，νx为温度场的平移速度，其负值为试样一侧的缩短速率；Qs、σe 和
分别为产热功率、等效应力和等效应变率。通过假设大变形体系内塑性变形功以热的形式全部转化，此时体系内产热功率的数值与σe 与
的乘积相同。
根据最大熵产生原理，并结合式(1)～(3)，大变形系统的熵产生率θ(单位时间和体积内的熵产生率)可以表示为：
式中，T′x为界面法向(轴向x)的温度梯度；νx1为温度场的平移速度(与焊接过程中的轴向缩短速度大小相等方向相反)。
基于模型假设(5)可知，大变形体系的稳态阶段对应大变形体系熵产生率θ 最大的状态，可以θ为目标函数，φ为限制函数，基于双曲正弦本构方程，确定函数自变量中的独立变量仅为T、
或T、σz，因此可依据拉格朗日乘子法和变分法，建立下式：
式中，q′x 是q 对x 的导数。
在旋转摩擦焊过程中，大变形体系的产生代表着焊接界面的消失[21]，此时的体系为在物质和动力学行为上具有连续性的塑变整体，即界面两侧的温度梯度T′是连续的。而异种材料接头界面处温度的峰值在靠近黏塑性区宽度较大的一侧。图3为接头处温度分布的示意图，温度峰值假定在x=(B1-B2)/2(B1>B2)处，此处的温度梯度为0，因此大变形系统中的温度分布可写为：
图3 异种旋转摩擦焊接头温度分布示意图
Fig.3 Schematic of the temperature distribution of dissimilar rotary friction welded joints
一维大变形体系的熵产生率θ 可通过以下方法推导，由式(2)～(6)可知：
式中，
和
是T、T′x、T″x 和
的平均值，其位于x=0 处，且σz 不随x 变化，即σz=σza。
由于大变形系统中平行摩擦方向的剪切应变速率
沿x 向的变化很小，因此
可被忽略，且结合式(1)～(3)，式(9)可简化为：
而进一步将等式(7)和(8)代入式(10)中，可得出熵产生率θ为：
1.3 异种摩擦焊接头稳态二维OZ-MEPP 解析模型
在实际焊接过程中，大变形区域的金属是沿着y 向以OC1 界面被挤出形成飞边，这使得实际接头的大变形区相对摩擦方向z 向发生偏转，将大变形区域划分为低速压扭变形区LCDZ 与高速剪切变形区HSDZ 两个区域。图4为由旋转摩擦焊产生的大变形体系的二维模型。在这个模型中，飞边(图4 中的B1H1O，B2H2O)方向对熵流(或熵产生)的影响被认为是在y0 方向上，因此，异种旋转摩擦焊过程所使用的样本被假设为一对矩形块，如图4a 所示。图4b表示两侧垂直于摩擦方向z 1/4 的横截面。LCDZ 与弹性区Z1 相邻，HSDZ 区与飞边区Z3 相邻，LCDZ主要承受扭转变形，此区域为实际接头中所谓的“热力耦合区”(A1B1C1、A2B2C2)；HSDZ 主要承受沿摩擦焊方向的高速剪切变形，此区域对应实际接头中所谓的大变形区(B1OC1、B2OC2)。
图4 异种旋转摩擦焊接头的大变形系统二维物理模型示意图：(a)两侧接头示意图；(b)两侧接头二维模型剖面图
Fig.4 Schematic of the 2D model of the large deformation system of a heterogeneous rotary friction welding head:(a)schematic diagram of the joint on both sides;(b)2D model profile of the two side joints
基于引入飞边对大变形体系的影响所导致的大变形区发生区域分化以及流动的偏转，在垂直于摩擦方向z 的方向上，建立稳态阶段二维模型，并对二维模型进行以下假设开发。
(1)从大变形区挤出的金属以速度νy1、νy2 沿y0(或-y0)方向流动，形成飞边，如图4 所示。
(2)异种大变形体系统中HSDZ 和LCDZ 之间的界面为H-L 界面。体系等温线与H-L 界面平行，但与摩擦界面成一定夹角，实际摩擦界面法向定义为x0 方向，总的热交换方向沿H-L 界面法线方向进行，并继承一维模型相关公式。
(3)HSDZ 仅具有平行于摩擦方向z 的剪切应力和应变率，即σz 和
。变形功与Qs 相当，LCDZ 仅具有扭转和压缩变形，其对应的变形功由焊接压力P和一侧缩短率-νx0 提供。且LCDZ 具有平行于飞边生长方向的剪切应力和应变率(即：沿y0 方向的流速随x0 线性减小，在弹性区与大变形区界面处时减小为零)。
(4)一维模型假设中的(5)和(6)仍然适用于二维模型。
(5)继承一维模型假设(1)，将一维模型的假设(4)更改为：忽略弹性区黏塑性变形和产热。此外，由于大变形区的产热主要集中在HSDZ，LCDZ 并未对大变形区产热做出贡献，因此大变形体系的熵产生主要集中在HSDZ，忽略了LCDZ 的熵产生。
由图4b 所示二维模型，以最快热传导方向建立坐标系xOy (虚线所示)，此坐标系与一维坐标系相同，因此式(11)仍对二维奏效，但二维的熵产生率必须采用虚线坐标系对应的相关参数(黏塑性区宽度B，轴向缩短速度νx)。然而xOy 坐标系中黏塑性区宽度B 与轴向缩短速度νx 不方便测量，因此可以根据xOy 轴与x0Oy0 的位向关系，可得出ν′x=νx1/cosα1、B′1=B1cosα1 关系式，并将其带入式(11)中可得二维模型熵产生率为：
由于大变形体系处于稳定状态，式中νx1 是一个常数且与位置无关。再由一维中的假设(6)和二维中的假设(4)可知，式(12)是与位置无关的常数体系，体系内的熵产生率在各区域均匀相同，不随x 变化。
在整个大变形区域中，黏塑性区较宽的一侧产热较多，这是因为其有较宽的塑性变形产热区域。根据以界面为原点，沿轴向x=(B1-B2)/2 处为温度的峰值(图3)这一假设，建立焊接界面与体系平均熵产生率关系如下式所示：
式中，
和Tx1，Tx2 分别为温度峰值处x1=(B1-B2)/2 与x2=B1/2处σz，
和T 的值，由于摩擦焊试件所受切向剪切应力在不随x 轴变化，因此可得出σz,x1=σz,x2，且根据ν′x=νx1/cosα1 和B′1=B1cosα1，式(13)可变化为：
大变形体系内，切向应变速率
的变化远小于温度的变化，因此
≈1[19]，其次，根据图4b 所示温度分布的映射关系(即等温线方向平行实际弹/黏界面B1C1)可知，Tx1 与Tx2 在实际轴向x0 的位置分别为x0=(B1-B2)/2、x0=B1/2，因此，式(14)进一步推导为：
由于
是接近于0 的数，可根据泰勒展开(cosx 部分保留前两项，sinx 部分等价于x)，并基于式(7)和(8)，上述等式可表达为：
经过进一步的等式变换，可将上述等式变化为：
将式(12)中所表示的熵产生率表达式带入到上述等式中，最终可以得出针对二维大变形体系下的空间边界条件，如下式所示：
根据图4 所示，式(18)中的cosα1 可以表示为：
式中，R为试样的半径。
综上，摩擦焊大变形体系的一、二维模型为分析实际摩擦焊过程提供了重要的理论基础。
1.4 异种接头稳态二维OZ-MEPP 定解模型
建立旋转摩擦焊稳态体系模型，需要结合实际旋转摩擦焊稳态阶段的具体特征，构建二维定解模型。针对模型分析作出如下定解模型假设。
(1)在大变形体系内，热传导最快的方向x 轴发生偏移，偏向轴向一侧，这种偏移主要是由于焊接界面边缘区域具有更高的线速度和更厚的黏塑性变形区域，因此产热更多。
(2)继承二维模型假设(2)中的区域划分，但不同的是，尽管旋转摩擦焊中的LCDZ/HSDZ 界面(H-L界面)仍然是B1C1，但定解模型中此界面具有弯曲特征(图5)，而非二维模型所示的平直特征(图4)。
图5 大变形体系的二维定解模型物理模型异种接头示意图与剖面图：(a,b)一侧接头；(c,d)另一侧接头
Fig.5 Physical modelling of heterogeneous joints with 2D fixed solution models for large deformation systems:(a,b)one side of the joint;(c,d)the other side of the joint
(3)二维模型中假设(3)对LCDZ 与HSDZ 变形方向的假设，以及LCDZ 与HSDZ 变形功分别对应垂直焊接界面C1O、C2O(即轴向)压力做功与扭矩做功的假设对二维定解模型保持不变。
(4)大变形区与弹性区的弹黏界面A1B1(A2B2)上沿飞边的流速(νy1,x0=B1、νy2,x0=B2)为0，且LCDZ 区内且垂直焊接界面各面内向飞边处的平均流速沿x0向呈线性分布。
(5)由于实际焊接过程HSDZ 区的剪切变形具有三维特征，在此区域的热传导、对流所对应的飞边生长或熵流具有不同方向上的变化，这一特征使得HSDZ 内熵产生率θ 的约束条件变弱，因此造成自发向熵值减小的方向上移动，进而引发熵产生率θ减小。在二维模型中，可以通过y0Oz 平面相对某一飞边生长方向的平均共轭速度作为当量摩擦速度，对上述减小的熵产生率θ 适当估计。
(6)继承二维模型的假设(1)、(4)和(5)。
根据假设(3)可知，式(12)和(18)仍对二维定解模型大变形系统有效。由于假设(2)，H-L 界面为曲线，这意味着角度α 随着半径r 而增加(参考图5b 和d)，因此式(19)不适用于旋转摩擦焊稳态阶段二维最大熵产生定解模型。图4 表明，二维模型中理论法向x与实际法向x0 之间的夹角α 实际上是由流体连续性和质量守恒定律决定的，而二维定解模型中的大变形体系α1 的平均值α1，a 仍可根据流体的连续性和质量守恒定律计算出来，即：
由式(19)和(20)可知，α1 也可以表示为速度的函数：
再根据流体连续性和质量守恒定律，轴向缩短的体积等于焊接过程期间沿着y0 向流出的体积(即飞边生长的体积)，即：
式中，R为试样的半径；νy1为二维定解模型中的温度场的平移速度(即轴向缩短速率)；νy1,R为沿y0 向金属流出的速度(即飞边生长的平均速度)。
利用上述方法，旋转摩擦焊二维定解模型的α1也可以描述为速度的函数。α1,a 可以表示为速度的函数，根据式(22)获得以下关系：
因此旋转摩擦焊二维定解模型中的cosα1 表示为：
1.4.1 二维大变形体系中LCDZ 的变形特征
对旋转摩擦焊二维模型中的变形区LCDZ 进行分析。首先基于二维定解模型的假设(3)，在LCDZ(图5 中A1B1C1、A2B2C2 区域)中的变形仅具有沿着飞边生长方向y0 的剪切。LCDZ 中y0 向的剪切变形对应的变形功由焊接界面的法向压力P 和一侧缩短速率νx1 提供，即：
式中，
和σy1 分别为LCDZ 内的应变速率和剪切应力，其方向平行于飞边生长方向y0；∧′1为系统内LCDZ 区的体积。
根据假设(6)和图5 所示，在大变形体系中沿y0方向的流速随x0 方向呈线性减小趋势，并在B1 时减小至0，因此，等式(25)中ε˙y1 可表征为：
式中，νy1,a 是νy1,r 的平均值；νy1,r为在LCDZ 中距离中心轴r 的点沿y0 方向上的流速，νy1,a 的位置被假定为x=B1/2 处。而LCDZ 中νy1,r 的平均值可以估计为νy1,R。根据式(26)可知，νy1,a 通过下式表述：
结合式(26)和(27)，将
表示为：
在进行LCDZ 的体积∧′1 计算时，可先计算HSDZ 的体积∧″1：
由二维定解模型假设(3)可知，HSDZ 中无沿径向方向的剪切变形，因此HSDZ 内距轴线为r 的(图5中的E1F1，E2F2)线上各点具有y0 向相同的流速νy1,r。再根据流体的连续性和物质的质量守恒定律，此时的νy1,r 应满足：
式中，νy1,R为试样半径为R 处各点沿y0 向的流 速。
根据图5 所示，对于半径为r 的圆柱体，轴向x0缩短量与径向y0 挤出量保持守恒，即轴向缩短的体积与飞边生长的体积相同，所以νy1,r 还应满足：
式中，Br 代表半径为r 时焊接界面距离H-L 界面的值。再由式(30)可得：
将式(27)和(32)带入到式(30)中，可得：
此时根据式(29)和(33)可得出HSDZ 区的体积∧′1为：
将计算结果与图5 所示的LCDZ 和HSDZ 体积比较可得出：
联立式(28)和(34)并将其代入到式(25)中可得出LCDZ 的σy1 值为：
此推导过程表明，大变形体系中LCDZ 内的应力σy1 与应变速率ε˙y1 是具有能量与质量守恒定律获得的，其结果不受旋转摩擦焊熵流的影响。
1.4.2 二维大变形体系中HSDZ 的变形特征
根据二维假设(3)可知，在大变形体系中，HSDZ只存在沿着摩擦方向z 向的剪切变形，高速剪切带只存在沿着摩擦方向的剪切，其相应的剪切应力和应变率可以根据高速剪切带的θ 计算，HSDZ 的θ是大变形体系的θ，二维大变形体系的θ 由等式(12)所示，在计算HSDZ 区的θ 时，由于H-L 界面和飞边生长的径向方向导致νx1 不是冷金属流入二维定解HSDZ 的流动速度，如图6 所示。图6a 展示出了实际的二维定解模型的HSDZ，其中焊接界面是曲面，并且飞边在径向方向上生长。二维HSDZ 如图6b 所示，焊接界面为矩形，且飞边仅沿y0 方向生长。由图6可知二维定解模型飞边生长的横截面积(H1B1I1I′1OH′1)与二维飞边生长的横截面积(H′3H′2I′2I′3)相等，二维定解和二维模型下的冷金属流出速度相同，可将二维定解模型HSDZ 流出的冷金属速度νy1,R，用于表示二维定解模型HSDZ 的θ。因此需要将式(12)中的速度发生替换，以适用于二维定解模型的θ 转化计算。
图6 旋转摩擦焊摩擦焊的大变形区模型示意图：(a)二维定解模型；(b)二维模型
Fig.6 Schematic diagram of the large deformation zone model for rotary friction welding friction welding:(a)2D fixed solution model;(b)2D model
根据图4b 可知，存在有νx1/νy1=B1/R 的关系，因此式(12)可以改写为式(37)。式(37)即为可适用于图6 二维模型下的θ 表示。
进一步可对二维定解中的θ 进行表示，由于二维定解和二维模型下的冷金属流出速度相同，因此，可得出二维定解下θ2 的表达式为：
式(38)可应用于图6a 所示二维定解模型熵产生率的计算，然而金属的流出速度在现实测量中较难实现，结合图6 所示，对于二维定解模型中的-νx 却可根据试样实际的缩短率来测量，因此式(38)需进一步优化。依据质量守恒定律，二维定解模型以及二维模型中的高速剪切带中的流动应当满足以下等式：
由于两模型中的冷金属流出速度相同(νy1,R=νy1)，因此式(39)和(40)可表示为νx1,2=2νx1，即二维定解模型下的缩短速率为二维模型的两倍，所以结合此关系并根据式(12)得出HSDZ 的熵产生率表达式为：
式(41)所表示的是HSDZ 的熵流表达式，然而熵产生率也可以用动力学行为表达，因此可将二维模型下的熵产生率θ 近似表示为：
根据二维定解模型假设(5)，在应对因熵流在不同方向上的变化而造成熵产生率θ 减小的方面上，需用焊接界面上相对于某一飞边生长方向(图7)的平均共轭速度U′来计算。在式(42)中，二维模型在平行于焊接界面的截面内，仅考虑了飞边生长方向(y0向)上的熵流，但可认为剪切方向及飞边生长方向的熵各占总熵流的1/2，因此，对HSDZ 区的应变速率的估计可表现为：
图7 焊接界面边缘任意一点C′相对y0 方向的边缘焊接速度U′示意图
Fig.7 Schematic of the conjugated velocity U′at point C′corresponding to the flash direction
根据图7 的示意图所示，可以得出焊接边缘速度U′的表达式为：
式中，ω1为单侧旋转摩擦焊的焊接旋转速度；β为y0向与任意飞边取向OC′的夹角。综上所述，根据式(43)和(44)可建立二维的HSDZ 中熵产生率表达式：
式(45)仅适用于二维模型下的熵产生率，要计算动力学条件下二维定解的熵产生，需要将式(45)发生转化。在将二维模型的熵产生率转化为二维定解模型熵产生率时，根据图6 可以得出，二维定解模型下的大变形体系的体积是二维模型下的一半，因此，θ2 应当为θ2D 的一半，即：
对比式(43)和(46)，可以得到：
结合式(42)和式(46)进一步可推导出：
基于上述推导结果，可以建立出二维定解模型下的旋转摩擦焊热力过程解析方程。虽然旋转摩擦焊的大变形中两个不同区域HSDZ 和LCDZ 具有不同的应变速率，但两个区域的变形可以用相同的本构方程来描述，因此本文采用了适用宽范围的应力和应变速率的反双曲正弦本构方程来描述，其表达式为：
式中，σe 和
分别为等效应力和等效应变率；α、n 和A为材料参数；Z为Zener-Hollomon 参数；Q为激活能；R为气体常数。
在RFW 过程中，LCDZ 和HSDZ 的
值很高，尤其是HSDZ。两区域的塑性流动均可使得Z/A>>1[22]，因此式(49)进而可得以下简化结果：
根据关系式：
，和νx1,3=2νx1，将式(28)和(36)带入到式(50)中即可得到LCDZ 的本构方程表达式：
为了与气体常数R 作区分，将等式中试样半径改为R0。同样根据关系时：
，将式(47)和(48)代入到式(50)中即可得出HSDZ 的本构方程表达式：
此外，替换式(24)到式(12)，旋转摩擦焊二维定解模型的边界条件可以改写为:
综上，已得出某一具体的旋转摩擦焊单侧稳态阶段热力过程解析模型，另一侧的推导过程与上述推导相似。通过相似的推导过程，可得出另一侧的HSDZ 和LCDZ 的本构方程以及边界条件。
另一侧的HSDZ 的本构方程表达式为：
另一侧的LCDZ 本构方程表达式为：
另一侧的二维定解旋转摩擦焊模型的边界条件为：
另一侧的旋转摩擦焊模型的熵产生率为：
结合式(53)与(56)，可得出：
1.5 异种旋转摩擦焊接头界面焊合判据的建立
与同种材料旋转摩擦焊相比，异种材料旋转摩擦焊接头界面更加复杂。在稳态阶段，异种界面产生大量的塑性变形和热量，界面两侧的少部分金属呈黏塑性状态，流动性增强。随着过程进行，远离界面的温度迅速升高，更多金属在高温、高应变速率和高应力下呈现黏塑性状态，温度、应变速率和应力沿轴向形成梯度分布。由于异种材料在比热、热传导和热扩散等物理性质上的差异，导致两侧温度梯度不一致。此外，材料的屈服强度随温度变化不同，进而使得两侧的塑性变形程度不同，产生明显的应力和应变速率梯度。在适当的温度和焊接压力作用下，摩擦焊大变形区两侧的黏塑性区宽度有所区分。
当界面温度和流速保持连续时，塑性流动均匀，材料充分接触，界面温度和应力促使两侧发生原子级互扩散，形成牢固的金属键，最终实现焊合。因此，可以将满足界面流速、温度及熵产生率连续的条件作为异种旋转摩擦焊界面焊合的热力学判据，并依据上述推导建立焊合判据的表达式。
当摩擦界面处两侧的温度相同时，得到界面温度连续的表达式为：
式中，
和
分别代表着界面处两侧不同材料的温度；Ta 代表界面平均温度。
当摩擦界面两侧的流速达到一致，可得出如下关系式：
式中，ε˙e1 和ε˙e2 分别为焊接界面两侧的等效应变速率；ω为焊接转速。并结合HSDZ 的等效应变速率等式，得出界面处流速相同条件为：
在RFW 过程中，当处于稳态阶段的界面处两侧熵产生率相同，可作为焊合的判断依据，由式(41)和(57)可得出表达式为：
综上所述，对于某一具体的异种旋转摩擦焊稳态阶段，在式(51)、(52)、(54)、(55)、(58)、(61)、(62)所构成的方程组中，已知参量包含：被焊试件半径R0，被焊两种母材的本构参量α1、n1、A1、Q1(α2、n2、A2、Q2)，气体常数R，热物性参数κ1(κ2)、λ1(λ2)以及工艺参数ω、P，而νx1、νx2、B1、B2、Ta、ω1、ω2为7 个未知参量，其与方程数相等，即方程可解。
2 实验结果及讨论
表1为TC4 被焊试样的几何参数和热物性参数[23]，表2为TC17 被焊试样的几何参数和热物性参数[24]。
表1 TC4钛合金试样相关参数[23]
Tab.1 Parameters of the TC4 titanium alloy samples[23]
ParameterValueParameterValue Densities,ρ/(kg·m-3)4 440Constitutive parameter,α/MPa-10.009 2 Heat capacity,c/(J·kg-1·K-1)807.05Constitutive parameter,n3.862 6 Thermal diffusion,κ/(m2·s-1)4.73×10-6Constitutive parameter,A/s-11.128 8×1020 Thermal conductivity,λ/(W·m-1·K-1)16.95Deformation activation energy,Q/(kJ·mol-1)456×103 Specimen radius,R0/mm12.50Gas constant,R/(J·mol-1·K-1)8.314
表2 TC17钛合金试样相关参数[24]
Tab.2 Parameters of the TC17 titanium alloy samples[24]
ParametersValueParametersValue Densities,ρ/(kg·m-3)4 650Constitutive parameter,α/MPa-10.012 93 Heat capacity,c/(J·kg-1·K-1)717.50Constitutive parameter,n3.036 48 Thermal diffusion,κ/(m2·s-1)5.32×10-6Constitutive parameter,A/s-11.1×108 Thermal conductivity,λ/(W·m-1·K-1)17.75Deformation activation energy,Q/(kJ·mol-1)214.322×103 Specimen radius,R0/mm12.50Gas constant,R/(J·mol-1·K-1)8.314
2.1 TC4/TC17 异种接头黏塑性区宽度验证
对异种旋转摩擦焊稳态阶段响应参量黏塑性区宽度B 采用变焊接压力实验组和变焊接转速实验组进行验证。两侧黏塑性区宽度代表着接头大变形区域的宽度，是大变形体系熵产生率表达的重要组成参数，也是本文摩擦焊稳态最大熵产生模型中的一个重要参数，对黏塑性区宽度的验证能够有效验证模型的准确性。
2.1.1 不同焊接转速下黏塑性区宽度验证
图8 是恒定焊接压力60 MPa 下，焊接转速分别为600、900、1 200、1 500、1 800、2 100 r/min 时，响应参量黏塑性区宽度B1(TC17)和B2(TC4)随焊接转速的变化关系。图中黑色线条代表模型对两侧黏塑性区宽度的计算值，红点为实验值。整体上，黏塑性区宽度随着转速的增加而增大。这主要是因为转速提高会增强界面的摩擦产热，导致界面温度升高进而增大了B1 和B2，值得注意的特征转速600 r/min条件下，TC17 一侧的B1 值显著偏高，仅次于高转速下的数值。其原因可能在于低转速下升温速率较缓，焊接时间较长，TC17 侧能够充分进行轴向热传导和塑性变形，从而形成较宽的大变形区。然而对比两侧黏塑性区宽度发现，在任意转速条件下B1＞B2。根据表1 与2 可知，这种差异主要源于两种材料的热物性和变形特性。一方面，TC17 的热扩散系数与热传导系数均高于TC4，这使得在焊接过程中TC17 侧沿着轴向传递的热量要高于TC4 侧。另一方面，TC17 具有较低的变形激活能，意味着其更容易在高温条件下发生塑性变形。高热传导性和低变形激活能的共同作用造成TC17 侧的黏塑性区宽度B1 大于TC4 侧B2。根据模型计算同样也表现出相同的数值差异，对比两侧计算值与实验值的相关性，二者计算值与实验值平均偏差为9.3%和8.3%，计算结果具有较好的准确性。
图8 焊接压力60 MPa，不同焊接转速下TC4/TC17 黏塑性区宽度B 实验值与计算值对比：(a)TC17(B1);(b)TC4(B2)
Fig.8 Comparison of the correlations between the experimental and calculated values of the TC4/TC17 viscoplastic zone width B with various welding speeds:(a)TC17(B1);(b)TC4(B2)
2.1.2 不同焊接压力下黏塑性区宽度验证
图9 是恒定转速1 200 r/min 下，焊接压力分别为60、90、120、150、180、210 MPa 时，响应参量黏塑性区宽度B1(TC17)和B2(TC4) 宽度随焊接压力的变化关系。图9a 和b 中黑线和红点随着焊接压力增加，黏塑性区宽度B1，B2 呈递减趋势。焊接压力的增加，使得焊接界面温度降低，从而抑制了黏塑性区的扩展，导致两侧宽度均呈下降趋势。此外，由于母材两侧不同的热物性和变形特性，保持在相同焊接参数下时，由图可知B1＞B2。模型计算结果与实验值的对比显示，TC17 和TC4 两侧的计算值与实验值的平均偏差分别为8.1%和6.5%，表现出较好的准确性。
图9 焊接转速1 200 r/min，不同焊接压力下TC4/TC17 黏塑性区宽度B 实验值与计算值对比：(a)TC17(B1);(b)TC4(B2)
Fig.9 Comparison of the correlations between the experimental and calculated values of TC4/TC17 viscoplastic zone width B with various weld pressures:(a)TC17(B1);(b)TC4(B2)
2.2 TC4/TC17 异种接头轴向缩短速率验证
对异种旋转摩擦焊稳态阶段响应参量轴向缩短速率νx 同样采用变焊接压力实验组和变焊接转速实验组进行验证。两侧轴向缩短速率代表着接头温度场的推进速度，同样也是大变形体系熵产生率表达的重要组成参数，对轴向缩短速率采用高速摄影机进行拍摄测量，验证结果能够有效验证模型的准确性。
2.2.1 不同焊接转速下轴向缩短速率验证
图10 是恒定焊接压力60 MPa 下，焊接转速分别为600、900、1 200、1 500、1 800、2 100 r/min 时，响应参量轴向缩短速度νx1(TC17)和νx2(TC4)随焊接转速的变化关系。图中黑色线条代表模型对两侧轴向缩短速度的计算值，红点为实验值。轴向缩短速度与试样轴向缩短量相关，轴向缩短量与飞边挤出量对应，而飞边挤出量取决于母材在焊接温度下的屈服强度。根据实验结果可知，TC4 的挤出量要远大于TC17，因此在相同焊接参数下，νx1 要大于νx2。从图中可以观察到，TC17 侧的νx1 实验值随焊接转速呈现明显的递增趋势，这与模型计算结果一致。然而，TC4 侧的νx2 表现出先上升后下降的非线性变化规律，可能的形成原因为： 在转速低于1 200 r/min的特征参数范围内，焊接转速的提升有效缩短了焊接时间，增加了界面摩擦产热和轴向缩短量，从而使缩短速度显著提升；但当转速超过1 200 r/min后，随着转速的增加，虽在一定程度上减少了焊接时间，但在实际焊接过程中，高转速的热输入下显著减少了初始摩擦产热阶段的时间，相对延长了稳态塑变产热阶段的时间，然而此时轴向缩短量的增量并没有发生显著变化，因此缩短速度随着转速的升高呈递减趋势。而TC17 一侧由于母材有较好的热传导和热扩散能力以及较低的变形激活能，使得缩短量增量持续上升，因此νx1 随转速变化表现出上升的趋势。模型计算值与实验值的相关性分析显示平均偏差为8.7%和7.2%，具有良好的准确性。
图10 焊接压力60 MPa，不同焊接转速下TC4/TC17 轴向缩短速率νx 实验值与计算值对比：(a)TC17(νx1);(b)TC4(νx2)
Fig.10 Comparison of the correlations between the experimental and calculated values of the TC4/TC17 shortening rate νx with various welding rotation:(a)TC17(νx1);(b)TC4(νx2)
2.2.2 不同焊接压力下轴向缩短速率验证
图11 是恒定转速1 200 r/min 下，焊接压力分别为60、90、120、150、180、210 MPa 时，响应参量轴向缩短速度νx1(TC17)和νx2(TC4)随焊接压力的变化关系。随着焊接压力的增加，黏塑性金属被大量挤出，相比于变转速实验组，变压力实验组轴向缩短速度要更大。值得注意的是，缩短速度的变化并不是随着焊接压力的递增而递增。当焊接压力达到150 MPa为变化特征点，往后随着焊接压力的增加呈平稳趋势，这种先递增后趋于稳定的变化趋势在计算值中同样也有所表现。且模型计算值与实验值的平均相对误差为9.1%和8.0%，具有良好的准确性。
图11 焊接1 200 r/min，不同焊接压力下TC4/TC17 轴向缩短速率νx 实验值与计算值对比：(a)TC17(νx1);(b)TC4(νx2)
Fig.11 Comparison of the correlations between the experimental and calculated values of the TC4/TC17 shortening rate νx with various welding pressures:(a)TC17(νx1);(b)TC4(νx2)
2.3 TC4/TC17 异种接头界面温度验证
2.3.1 不同焊接转速下界面温度验证
通常转速的增加能够增大摩擦产热，进而使得界面温度上升。图12 是恒定焊接压力60 MPa 下，稳态阶段界面温度随转速变化，实验值与计算值相关性对比结果。由图可知，随着转速的增加，通过计算与实验得到的界面温度均呈单调递增的趋势，且二者模型计算的结果与实验结果进行对比，平均偏差为4.3%具有较好的准确性。
图12 焊接60 MPa，不同焊接转速下TC4/TC17 界面温度T 实验值与计算值对比
Fig.12 Comparison of the correlations between the experimental and calculated values of the TC4/TC17 interface temperature T with various welding rotation
2.3.2 不同焊接压力下界面温度验证
图13 是恒定转速1 200 r/min 下，稳态阶段界面温度随焊接压力的变化。焊接压力的增大，可使得黏塑性金属不断被挤出，进而降低界面温度。图中随着焊接压力的增大，实验与计算得出的界面温度单调减小，变化趋势一致，且平均偏差为3.0%，具有较好的准确性。
图13 焊接1 200 r/min，不同焊接压力下TC4/TC17 界面温度T 实验值与计算值对比
Fig.13 Comparison of the correlations between the experimental and calculated values of the TC4/TC17 interfacial temperature T with various welding pressures
2.4 其他材料验证
通过上述分析可知，由异种旋转摩擦稳态热力耦合解析模型预测的基本特征已通过TC4/TC17异种材料的焊接实验得到了验证，验证结果初步证实了解析计算具有一定的准确性。针对模型的可靠性，本文根据已有文献的TC4/SUS321 异种旋转摩擦焊界面温度的实验数据[25]，对模型进一步的验证，表3是SUS321 材料相关的参数[26]。
表3 SUS321不锈钢试样相关参数[26]
Tab.3 Parameters for the SUS321 steel samples[26]
ParametersValueParametersValue Densities,ρ/(kg·m-3)7 930Constitutive parameter,α/MPa-10.00 9 Heat capacity,c/(J·kg-1·K-1)500.12Constitutive parameter,n6.1 Thermal diffusion,κ/(m2·s-1)5.42×10-6Constitutive parameter,A/s-19.76×101 7 Thermal conductivity,λ/(W·m-1·K-1)21.50Deformation activation energy,Q/(kJ·mol-1)465×103 Specimen radius,R0/mm12.50Gas constant,R/(J·mol-1·K-1)8.314
图14 是恒定焊接压力150 MPa 下，稳态阶段界面温度随焊接转速的变化。由图可知，焊接转速的增大，使得界面温度呈先增大后减小的趋势，在600 r/min 下界面温度达到最小值，模型计算结果与实验结果的平均偏差为4.4%，具有较好的可靠性。
图14 TC4/SUS321 界面温度T 实验值与计算值随焊接转速的相关性对比
Fig.14 Comparison of the correlations between the experimental and calculated values of the TC4/SUS321 interface temperature T with the welding rotation
3 结论
(1)建立了一维通用的异种旋转摩擦焊分析模型来描述温度分布和熵产生的基本特征；其次，考虑到飞边方向对热传导和对流的影响，建立了和二维定解旋转摩擦焊的空间约束条件。
(2)通过描述异种旋转摩擦焊界面焊合过程，建立异种旋转摩擦焊焊合判据，结合二维定解模型和金属塑性本构方程，建立了旋转摩擦焊稳态热力耦合模型。该模型只需根据初始条件、试样尺寸、本构参数和热物性参数，就可直接预测大变形体系的轴向缩短速率νx、黏塑性区宽度B 和界面温度T 等相关响应参数，保证了模型的自洽性和可预测性。
(3)通过实验手段(TC4/TC17)与文献调阅(TC4/SUS321)验证了该模型。其中，实验手段以摩擦焊过程中3 个响应参量为验证对象，通过变压力组(60～210 MPa)和变转速组(600～2 100 r/min)焊接参数下，实验结果与计算结果变化趋势相近，两侧验证结果相对平均误差均在10%以内，具有良好的准确性；此外，查阅文献获得TC4/SUS321 异种接头界面温度随转速的变化规律。验证结果表明，实验值与计算值变化趋势相似，相对误差为4.3%，具有可靠性。
[1] NICHOLAS E D,THOMAS W M.A review of friction processes for aerospace applications [J]. International Journal of Materials and Product Technology,1998,13(1-2):45-55.
[2] KIMURA M,INOUE H,KUSAKA M,KAIZU K,FUJI A.Analysis method of friction torque and weld interface temperature during friction process of steel friction welding[J]. Journal of Solid Mechanics and Materials Engineering,2010,4(3):401-413.
[3] RYKALIN N N, PUGIN A I, VASIL’EVA V A. The heating and cooling of rods butt welded by the friction process[J].Svarochnoe Proizvodstvo,1959,6:42-52.
[4] MIDLING O T,GRONG Ø.A process model for friction welding of Al-Mg-Si alloys and Al-SiC metal matrix composites—I. Haz temperature and strain rate distribution[J]. Acta Metallurgica et Materialia,1994,42(5):1595-1609.
[5] SINGH S K, CHATTOPADHYAY K, PHANIKUMAR G, DUTTA P. Experimental and numerical studies on friction welding of thixocast A356 aluminum alloy[J].Acta Materialia,2014,73:177-185.
[6] LI W Y,SHI S X,WANG F F,MA T J,LI J L,GAO D L,VAIRIS A. Heat reflux in flash and its effect on joint temperature history during linear friction welding of steel[J]. International Journal of Thermal Sciences,2013,67:192-199.
[7] GENG P H,QIN G L,ZHOU J,ZOU Z D.Hot deformation behavior and constitutive model of GH4169 superalloy for linear friction welding process[J]. Journal of Manufacturing Processes, 2018, 32:469-481.
[8] LI W Y, WANG F F, SHI S X, MA T J. Numerical simulation of linear friction welding based on ABAQUS environment: Challenges and perspectives[J].Journal of Materials Engineering and Performance,2014,23:384-390.
[9] JAIN R, PAL S K, SINGH S B. A study on the variation of forces and temperature in a friction stir welding process: A finite element approach [J]. Journal of Manufacturing Processes, 2016,23:278-286.
[10] ZHAO Y, WU A P, REN J L, SATO Y S, KOKAWA H, MIYAKE M, YAN D Y. Temperature and force response characteristics of friction stir welding on Invar 36 alloy[J]. Science and Technology of Welding and Joining,2013,18(3):232-238.
[11] YANG X W, FENG W Y, LI W Y, XU Y X, CHU Q, MA T J,WANG W B. Numerical modelling and experimental investigation of thermal and material flow in probeless friction stir spot welding process of Al 2198-T8[J].Science and Technology of Welding and Joining,2018,23(8):704-714.
[12] UDAY M B,AHMAD FAUZI M N,ZUHAILAWATI H,ISMAIL A B. Advances in friction welding process: A review[J].Science and technology of Welding and Joining,2010,15(7):534-558.
[13] MAALEKIAN M. Thermal modeling of friction welding[J]. ISIJ International,2008,48(10):1429-1433.
[14] YANG Y C, CHEN W L, LEE H L. A nonlinear inverse problem in estimating the heat generation in rotary friction welding [J].Numerical Heat Transfer, Part A: Applications, 2011,59(2):130-149.
[15] MENDEZ P F,TELLO K E,LIENERT T J.Scaling of coupled heat transfer and plastic deformation around the pin in friction stir welding[J].Acta Materialia,2010,58(18):6012-6026.
[16] XIONG J T, LI J L, WEI Y N, ZHANG F S, HUANG W D. An analytical model of steady-state continuous drive friction welding[J].Acta Materialia,2013,61(5):1662-1675.
[17] BHAMJI I,PREUSS M,THREADGILL P L,ADDISON A C.Solid state joining of metals by linear friction welding:A literature review[J].Materials Science and Technology,2011,27(1):2-12.
[18] MARTYUSHEVLM,SELEZNEVVD.Maximumentropy production principle in physics,chemistry and biology[J].Physics Reports,2006,426(1):1-45.
[19] XIONG J T, ZHOU W, LI J L, ZHANG F S, HUANG W D. The thermodynamic analytical models for steady-state of linear friction welding based on the maximum entropy production principle[J].Materials&Design,2017,129:53-62.
[20] XIONG J T,YUAN L,DU Y J,SHI J M,LI J L.The thermodynamic analytical models for steady-state of continuous drive friction welding based on the maximum entropy production principle[J].Applied Thermal Engineering,2019,159:113840.
[21] ALEXANDROV S,ALEXANDROVA N.On the maximum friction law in viscoplasticity[J]. Mechanics of Time-Dependent Materials,2000,4:99-104.
[22] CHO S H,KIM S I,YOO Y C.The effects of deformation variables on hot workability of austenitic stainless steel[J]. Metals and Materials,1998,4:732-736.
[23] 徐勇，杨湘杰，何毅，杜丹妮.TC4 钛合金流动软化行为及本构模型研究[J].稀有金属材料与工程，2017,46(5):1321-1326.XU Y, YANG X J, HE Y, DU D N. Flow softening behavior and constitutive equation of TC4 titanium alloy during hot deformation[J].Rare Metal Materials and Engineering,2017,46(5):1321-1326.
[24] LI L, LI M Q. Constitutive model and optimal processing parameters of TC17 alloy with a transformed microstructure via kinetic analysis and processing maps [J]. Materials Science and Engineering:A,2017,698:302-312.
[25] LI X, LI J L, LIAO Z X, JIN F, ZHANG F S, XIONG J T.Microstructure evolution and mechanical properties of rotary friction welded TC4/SUS321 joints at various rotation speeds[J].Materials&Design,2016,99:26-36.
[26] NKHOMA R K C, SIYASIYA C W, STUMPF W E. Hot workability of AISI 321 and AISI 304 austenitic stainless steels[J].Journal of Alloys and Compounds,2014,595:103-112.