焊接过程伴随着焊接应力与变形，直接影响焊接结构的质量^[1]。焊接过程本身是温度-相变-应力耦合问题，焊接残余应力与变形是温度、相变和热应力耦合作用的结果。应力、相变对温度场的影响很弱，因此，焊接残余应力与变形的有限元计算主要是分析不均匀温度场造成的热应力。此外，国内外很多学者也在研究焊接过程中组织转变对应力的影响^[1-3]。

有限元分析方法中，对耦合问题的分析有2种方法：顺序耦合法和直接耦合法。顺序耦合是将温度场分析的结果作为载荷用于应力分析。直接耦合则是温度和应力同时进行计算。对于焊接问题，由于应力对温度场的影响很弱，计算焊接残余应力时可以不考虑，因此一般采用顺序耦合法。即首先进行非线性瞬态热分析获得整个焊接过程的温度场，再进行非线性静态应力分析，将瞬态热分析的节点温度作为应力分析的载荷，获得应力场分布。在焊接过程中，温度场是焊接过程的主要驱动力，它会引起相变、变形和残余应力^[4]。因此，焊接过程瞬态温度场的准确计算，对焊接残余应力及变形的分析至关重要。

焊接温度场主要取决于热源模型、散热条件以及材料特性。其中，热源模型代表焊接热能在工件上的分布，是温度场分析的基础。常用的热源为体热源模型，如半球、椭球、双椭球热源等^[5]，也有学者在这些热源模型基础上进行了改进^[6]。

在仿真过程中，一般根据焊接方法建立相应的热源模型类型，再根据实际焊缝形状，对热源参数进行调整，以匹配实际焊缝的边界。即使热源功率大小相同，热流密度分布不同时，得到的计算结果差别也会很大，因此热源参数的获取是焊接温度场分析中非常重要的环节。

热源形状参数的重要标准就是获得与实际熔池一致的形状，一般主要指熔池宽度和熔池深度，通过测量试验焊件的截面焊缝形状获得。传统的分析方法中，热源参数通过试算获得。研究人员先设定一组初始热源参数，进行温度场仿真分析后，对比仿真结果与实际熔池的边界，然后再次进行调整，直至仿真结果与实际熔池边界吻合。整个过程较为繁琐，手动调整及结果后处理的工作量较大，增加了仿真的建模成本，效率和精度受到研究人员的经验影响很大。为了能够提高热源参数获取的效率和精度，国内外学者都做了不少相关的研究。

王煜等^[7]基于材料与温度无关的模型，使用解析法获得了镍基合金对接焊件熔宽和熔深等数据与热源参数之间的关系。郭晓凯^[8]以焊接试验温度测点数据误差最小为目标，使用模式搜索法反求了多丝埋弧焊的焊接热源模型参数。李培麟^[9]使用多种智能优化算法对三丝焊接热源参数进行了估计。Jia等^[10]基于准稳态的假设，在二维模型上对双椭球热源参数进行了估计。Ficquet等^[11]使用热源建模工具(HSMT)对热源参数进行了估计，该工具假设焊缝在焊接开始时全部出现。Shan等^[12]使用单元生死的方法对比了材料逐步出现和同时出现的温度场，结果表明未使用单元生死时，温度场的峰值小于使用单元生死的模型。这是由于未考虑单元逐步出现时，热源的部分热量会传到热源前面的材料，导致峰值温度降低，而这些材料实际是还未出现的。因此，焊缝一直存在对于自熔焊是合理的，但是对于填充焊则会导致热量分布与实际有较大差别。

综上所述，试算法获得热源参数过程繁琐，手动调整工作量较大，且效率和精度主要受到研究人员经验的影响。基于二维模型或准稳态假设时，未考虑填充焊接过程中热源前面的材料绝热效应，会导致温度场的峰值温度偏低。针对现有研究中存在的问题，本文提出热源形状参数优化模型，在热源参数获取过程中，使用多软件协同计算的方法，实现参数自动求解。最后以实际焊接试验为案例，使用提出的优化方法获取热源参数进行温度场模拟，并与试验结果进行对比验证。

1 热源参数优化模型的建立

在所有体热源模型中，最具代表性的是Goldak提出的双椭球热源模型^[13]。双椭球热源模型广泛应用于各种熔化焊接的温度场模拟中。沿着焊接方向，双椭球热源分为两部分，如图 1所示。热源内部任意一点(x, y, z)热流密度的表达式如下：

$ \begin{aligned} &q_{\mathrm{f}}=\frac{6 \sqrt{3} f_{\mathrm{f}} Q}{a b c_{\mathrm{f}} {\rm{ \mathsf{ π} }} \sqrt{{\rm{ \mathsf{ π} }}}} \mathrm{e}^{-3\left(\frac{x^{2}}{a^{2}}+\frac{y^{2}}{b^{2}}+\frac{z^{2}}{c_{\mathrm{f}}^{2}}\right)}, \\ &q_{\mathrm{r}}=\frac{6 \sqrt{3} f_{\mathrm{r}} Q}{a b c_{\mathrm{r}} {\rm{ \mathsf{ π} }} \sqrt{{\rm{ \mathsf{ π} }}}} \mathrm{e}^{-3\left(\frac{x^{2}}{a^{2}}+\frac{y^{2}}{b^{2}}+\frac{z^{2}}{c_{\mathrm{r}}^{2}}\right)}, \\ &f_{\mathrm{r}}+f_{\mathrm{f}}=2, \quad Q=\eta U I . \end{aligned} $

(1)

其中：f_f和f_r分别是热源前半部分和后半部分的能量分配系数；Q是焊接热量，η是焊接热效率；a、b、c_f、c_r为双椭球热源形状参数。

双椭球热源模型内部能量分布为Gauss分布，在体热源模型中具有代表性。当c_f=c_r时，该模型为椭球热源模型，当a、b、c_f、c_r均相等时，该模型为半球热源模型。双椭球热源模型形状参数较多，使用传统的试算法确定模型形状参数时，工作量较大，效率和精度受到研究人员主观因素影响较大。其他热源参数较多的组合热源模型也存在类似问题。针对该问题，本文以双椭球热源为例提出一种热源形状参数优化模型，该优化模型也可推广应用于其他类型热源的形状参数获取。

1.1 优化设计变量

以热源形状参数为优化设计变量，对于双椭球热源，有3个方向的尺寸参数：a为焊缝宽度方向上椭球的半轴长度，用于控制能量在焊缝宽度方向上的分布区域；b为焊缝深度方向上椭球的半轴长度，用于控制能量在焊缝深度方向上的能量分布区域；c_f为焊缝长度方向上前半个椭球的半轴长度，c_r为焊缝长度方向上后半个椭球的半轴长度。

热效率决定热源的总能量，形状参数决定能量在热源内部的分布。如果计算模型有对应的温度试验，实际计算时热效率根据距离焊缝较远的热电偶标定得到。当没有温度的试验数据时，热效率η也可作为一个变量。

1.2 目标函数与约束条件

目标函数：以焊缝横截面实际形状作为误差计算的基础。焊缝熔池横截面示意图见图 2。将仿真结果的焊缝横截面熔宽W_S、熔深D_S与实际焊缝横截面熔宽W_M、熔深D_M的误差的平方值进行加权求和，作为目标函数。

约束条件：热源形状参数越小，热源能量在该方向上分布越集中，热源形状参数越大能量分布越平缓。取值过小或过大都会使仿真结果与实际尺寸误差变大。Goldak^[5]指出，热源形状参数小于实际熔池尺寸约10%时，仿真计算得到的结果与实际熔池的尺寸结果吻合会比较好。综合考虑计算量及单元尺寸，热源半轴取值下限推荐在最小单元尺寸与该方向上熔池尺寸的1/2中选择较大值。最大值参考该方向上的实际熔池尺寸。形状参数在取值范围内的变化增量可参考该方向上最小焊缝单元尺寸。

综上，根据设计变量、目标函数及约束条件，提出热源参数优化问题数学模型如式(2)所示。其中：各变量下标1、2分别为取值范围上、下限。

$\min \ \ \ E=w_{1}\left(W_{\mathrm{S}}-W_{\mathrm{M}}\right)^{2}+w_{2}\left(D_{\mathrm{S}}-D_{\mathrm{M}}\right)^{2},\\ \ \ \ \ \ \ \ \ \ \ \ \ \ \text { s. t. }\left\{\begin{array}{l} a \in\left[a_{1}, a_{2}\right] , \\ b \in\left[b_{1}, b_{2}\right] , \\ c_{\mathrm{f}} \in\left[c_{1}, c_{2}\right] , \\ c_{\mathrm{r}}=k c_{\mathrm{f}}, \\ \eta=\left[\eta_{1}, \eta_{2}\right] , \\ a, b, c_{\mathrm{f}}, c_{\mathrm{r}}, k \in \mathbb{R}. \end{array}\right. $

(2)

其中：a、b、c_f、c_r为热源形状参数，η为热效率，不同的焊接方法热效率取值范围有所不同，一般电弧焊的热效率在0.65~0.9之间。k为后半部分半轴与前半部分半轴的比例，一般大于1，以模拟热源移动过程中前半部分和后半部分不同的能量分布。w₁, w₂为误差权重系数。下标S表示有限元仿真的结果，下标M表示实际焊缝的尺寸。

2 优化问题的求解

针对式(2)的优化问题，在Isight中建立优化模型，以最小化熔池边界误差为目标进行优化，使用多软件协同计算的方法进行求解。Abaqus用于对温度场的计算，热源通过Abaqus与Fortran关联的用户自定义子程序施加。Python用于对温度场结果进行后处理，得到有限元模型各个节点的温度结果，并传给MATLAB。MATLAB用于对节点温度进行数据处理，获得熔池的宽度和深度，并将仿真与实际熔池的误差结果传递给Isight。此外，整个计算和后处理软件调用也通过MATLAB编程实现。通过多软件协同自动计算，获得最优热源参数的同时可减少人工操作。完整的求解流程如图 3所示。

Isight将一组热源参数传递给MATLAB，在MATLAB中对用户自定义子程序中的热源参数进行修改，并施加在Abaqus中进行瞬态温度场分析。获得瞬态温度结果后，使用Python进行结果的后处理，将所有节点经历的峰值温度保存至文本文件中。然后在MATLAB中编写自动处理文件的程序，读取分析模型中所有节点对应的坐标，并将坐标信息与峰值温度信息相对应。按照节点坐标排序后，提取热源稳定的一个截面上所有节点的温度，通过线性插值的方式获得熔池宽度方向上峰值温度为1 400℃对应的坐标位置，该位置与焊缝中心的距离为仿真所得熔池半宽；同理获得熔池深度方向上峰值温度为1 400℃对应的坐标位置，该位置与焊缝表面的距离为仿真所得熔池深度。计算与目标熔池宽度和深度的误差平方和并进行加权，得到该组参数对应的误差。在MATLAB获取误差后传递给Isight进行保存。Isight将下一组热源形状参数传递给MATLAB，进入下一轮计算，直至求得最优解。

3 优化算例

设计Q490C对接板件焊接试验，试验过程中使用8个热电偶对板件的温度场进行测量，热电偶具体位置如图 4所示。为了验证试验的重复性，焊接4个试样，其中3个用于测量温度场分布，1个用于切割后获得焊缝熔池形状。板件尺寸为200 mm×200 mm×12 mm。中间有Ⅴ形坡口，角度22°，两端点焊。焊枪速度为0.28 m/min，焊接电压228.6 V，焊接电流312 A。对焊缝横截面进行打破抛光腐蚀，可得到熔池形状。测量焊缝熔池半宽为6.5 mm，熔池深度为8 mm。

按照板件尺寸及焊后焊缝的尺寸建立几何模型。由于对接板件左右对称，因此建立对称模型。如图 5所示，板件大小为200 mm×100 mm×12 mm，x方向垂直于焊缝，y方向沿厚度方向，z方向沿焊缝方向。根据焊缝尺寸建立有限元分析模型。焊缝及热影响区横截面上单元边长为0.8 mm。沿焊接方向单元边长为2 mm。单元总数56 200，节点数62 216。

仿真过程中，将材料与场变量相关联，实现材料填充效果的模拟，以考虑焊接热源前面的隔热现象。在温度场分析中，热源前面的绝热现象通过将热源前面的材料热导率减小2个数量级来实现。在Abaqus材料定义中设置与场变量(FV)相关，材料FV为2时，热导率为常温热导率的1%，模拟未激活的焊缝。FV为1时，为正常的焊缝材料参数。场变量在1和2之间时，通过线性插值来实现过渡。FV为0时，为母材的材料参数。场变量控制通过用户自定义子程序实现。具体实现过程见图 6。

本案例中，焊缝材料认为与母材一致。材料熔化潜热为260 kJ/kg，密度为7.84×10³ kg/m³，Poisson比为0.3。其余材料参数与温度相关，通过材料试验获得，具体材料参数如图 7所示。当温度超过熔点时，热导率设置为室温的4倍，来补偿熔池内部由于熔池内部液体流动导致的对流换热^[14-15]。样件初始温度30℃，在焊件表面设置综合对流换热系数20 W/(m²·℃)。

使用双椭球热源进行温度场模拟，如式(1)所示，通过用户自定义子程序实现热源的移动与加载。

热效率决定热源的总能量，形状参数决定能量在热源内部的分布。在进行温度场分析时，应先对热效率进行标定，确定热源总能量，再进行热源形状参数的选取。热效率标定时，假设其他条件不变，得到不同热效率下温度场的仿真结果，并与试验结果进行对比，选择仿真与试验最为接近的热效率。标定所用的温度测点应距离焊缝较远，测点距离焊缝越远，热源形状参数对测点温度的结果影响越小。此外，测点离焊缝太近时，测点温度结果对位置的偏差更敏感，对应的试验结果受到位置误差的影响较大。本例中，基于5号温度测点的试验结果对热效率进行标定，得到热效率为0.75。

对温度场分析结果进行后处理，提取所有节点在热源作用过程中经历过的峰值温度，绘制云图如图 8所示。图中灰色区域为节点经历的峰值温度大于1 400℃的区域，代表仿真所得熔池区域，黑色区域为节点峰值温度小于100℃的区域。从图中可以看出，对于本案例，在距离板件起点所在边缘24 mm处，熔池已达到稳定状态。峰值温度大于100℃的区域集中在距离焊缝中线50 mm范围内。

在使用热源优化方法时，可对温度场分析模型进行截取，以减少计算量，截取长度应保证计算过程中热源能够达到稳定。本案例中，沿焊接方向截取长度44 mm，横截面方向截取50 mm。

双椭球热源的后半部分半轴长c_r应大于前半部分半轴长c_f，以模拟热源移动过程中前端和后端不同的能量分布(前端梯度较大，后端梯度较小)。本案例中参考其他文献^{[5, 16-17]}的取值方法，令c_r=2c_f。熔池宽度和深度尺寸误差的权重系数取1，则对应的热源参数优化模型如式(3)所示。

$ \begin{aligned} \min E &=\left(W_{\mathrm{S}}-W_{\mathrm{M}}\right)^{2}+\left(D_{\mathrm{S}}-D_{\mathrm{M}}\right)^{2}, \\ \text { s. t. } &\left\{\begin{array}{l} a \in[2, 6], \\ b \in[4, 9], \\ c_{\mathrm{f}} \in[2, 5], \\ c_{\mathrm{r}}=2 c_{\mathrm{f}}, \\ \eta=0.75, \\ a, b, c_{\mathrm{f}}, c_{\mathrm{r}}, k \in \mathbb{N} . \end{array}\right. \end{aligned} $

(3)

计算过程中，对热源参数的取值进行离散化处理。针对确定的有限元模型，变化间隔小到一定程度时，所得结果相差很小，但是计算量会成倍增加。因此，可先设置较大的变化间隔，初步确定热源形状参数。在所得热源形状参数附近，将变化间隔变小，进一步细化形状参数的求解。本例中，先将热源形状参数的变化间隔设置为1。获得初步解对应的结果为a=5, b=7, c_f=2, c_r= 4。

在设计变量a、b的取值范围内，目标函数与设计变量的关系如图 9所示。可以看出，参数过小或者过大都会导致目标函数值(熔池尺寸综合误差)增大。

将设计变量a, b的取值范围缩小，令a∈[4.5, 5.5]，b∈[6.5, 7.5]，c∈[2, 2.5]。热源形状参数的变化间隔设置为0.25，得到最优解对应结果为a=5.5, b=7.25, c_f=2, c_r=4。

在新的取值范围内，目标函数与设计变量的关系如图 10所示。可以看出，设计变量的值相差0.25时，得到的目标函数的数值波动相比于设计变量相差1时要小。

热源参数最优解对应的熔池形状与试验结果对比如图 11所示。

由初步解和最优解计算得到的熔池尺寸与试验结果对比见表 1。可以看出，最优解得到的熔池尺寸与试验结果基本吻合。初步解与最优解得到的熔池尺寸相差也较小。

将热源参数最优解对应的热源加载至完整模型，进行温度场模拟。计算得到各个温度测点的仿真结果，并与试验结果进行对比，如图 12所示。图中S1、S2、S3分别代表3次试验的结果，FEA代表有限元仿真的结果。

3次温度试验中各个测点的峰值温度平均值如表 2所示，可以看出，距离焊缝中心越远，测点峰值温度越低。以试验测点峰值温度的平均值为基准，计算有限元温度场分析结果的相对误差。可以看出，除了测点6/8外，有限元仿真结果的相对误差均小于5%。距离焊缝越近时，测点峰值温度对位置越敏感，测点6/8误差较大的原因可能是由于实际测量测点位置的偏差导致。

通过案例的仿真与试验的结果对比可以看出，热源参数标定的结果较为准确，温度场模拟结果与试验吻合较好，验证了热源模型参数的有效性。

4 结论

本文针对热源形状参数优化问题，以双椭球热源模型为例，建立了热源形状参数优化的数学模型，并使用多软件协同计算的方法对问题进行求解。该方法也可用于其他类型热源模型参数优化。通过算例分析并与试验结果进行对比，表明本文提出的热源参数优化方法能够较为准确地获得热源参数。与传统的试错法相比，使用软件自动调整参数进行协同计算的方法可减少试错法调整热源参数的工作量，同时降低了焊接温度场模拟的效率与精度对研究人员经验的依赖性。