Abaqus 焊接残余应力仿真：顺序热力耦合与 DFLUX 移动热源实战

Abaqus 焊接残余应力仿真：顺序热力耦合与 DFLUX 移动热源实战

焊接是制造与结构工程中最常见的连接工艺，但焊接过程中局部高温加热与快速冷却会在接头区域产生显著的残余应力与变形。残余应力直接影响结构的疲劳寿命、抗应力腐蚀能力和尺寸稳定性，是桥梁、压力容器、船舶、轨道交通等关键结构必须评估的指标。借助 Abaqus 数值仿真，可以在物理试验之前预测焊接温度场与残余应力分布，为工艺优化提供依据。本文从工程实战角度，系统梳理 Abaqus 焊接残余应力仿真的完整流程：顺序热力耦合方法、DFLUX 移动热源子程序、温度相关材料设置以及关键收敛技巧。

一、工程背景：为什么要做焊接仿真

焊接过程本质上是一个高度非线性的瞬态热-力耦合问题。电弧或激光在极小区域内输入大量热量，使局部金属熔化，随后热量向周围母材传导散失，熔池凝固收缩。由于加热和冷却在空间上极不均匀，焊缝及热影响区会产生复杂的塑性应变，冷却后被”冻结”为残余应力。

实际工程中，焊接仿真主要解决三类问题：

  • 残余应力与变形预测：评估焊后结构的尺寸精度和装配可行性，指导反变形与夹具设计。
  • 工艺参数优化：对比不同热输入、焊接速度、焊接顺序对接头质量的影响，减少试错成本。
  • 服役性能评估：将残余应力作为初始状态，叠加服役载荷进行疲劳、断裂和应力腐蚀分析。

对于厚板、多道焊等复杂结构，物理试验成本高、周期长，数值仿真几乎是不可替代的前期评估手段。

二、理论基础：顺序热力耦合与 Goldak 热源模型

2.1 顺序热力耦合方法

焊接仿真有两种主流路线：完全耦合（directly coupled）与顺序耦合（sequentially coupled）。完全耦合在一个分析步内同时求解温度场与位移场，计算代价高、收敛困难，工程上较少采用。顺序热力耦合将问题解耦为两步：

  1. 非耦合热传递分析：先独立求解瞬态温度场，忽略应力/变形的影响；
  2. 应力分析：将上一步得到的温度场作为预定义场（热载荷）导入，求解热应力与残余应力。

这一方法的物理依据是：温度场对位移场有强影响（热膨胀、材料软化），但机械变形对温度场的反馈通常较弱，可以忽略。顺序耦合在保证精度的同时大幅降低计算成本，是焊接仿真的事实标准。

顺序热力耦合仿真流程

2.2 Goldak 双椭球移动热源模型

移动热源的描述是焊接温度场仿真的核心。Goldak 等（Goldak J., Chakravarti A., Bibby M., 1984, Metallurgical Transactions B, 15(2): 299-305）提出的双椭球热源模型是目前应用最广的体热源模型。该模型用前后两个半椭球描述熔池形状：焊接方向前方的半椭球较短（参数 c_f），后方较长（参数 c_r），以反映熔池”前陡后缓”的真实形态。

热流密度在椭球内按高斯分布衰减，前半部分与后半部分分别表示为：

前半椭球（热源前方）：

q_f(x,y,z) = (6·√3 · f_f · Q · η) / (a·b·c_f · π·√π) · exp( -3·(x²/a² + y²/b² + z²/c_f²) )

后半椭球（热源后方）：

q_r(x,y,z) = (6·√3 · f_r · Q · η) / (a·b·c_r · π·√π) · exp( -3·(x²/a² + y²/b² + z²/c_r²) )

其中 Q 为焊接热输入（电压×电流），η 为热效率，a、b 为椭球半长轴与半宽，f_f、f_r 为前后热量分配系数（满足 f_f + f_r = 2）。热源中心位置随时间沿焊接方向移动：x_c = x_0 + v·t。

Goldak双椭球移动热源示意

三、官方文档解读：DFLUX 子程序接口

在 Abaqus 中实现移动体热源，标准做法是编写 DFLUX 用户子程序。根据 Abaqus User Subroutines Reference Guide（文档路径：SIMACAESUBRefMap/simasub-c-dflux.htm），DFLUX 用于定义热传递或质量扩散分析中的非均匀分布热流（nonuniform distributed flux），可表示为位置、时间、温度、单元号、积分点号等的函数。子程序在每个热流积分点被调用。

官方给出的子程序接口骨架如下：

      SUBROUTINE DFLUX(FLUX,SOL,KSTEP,KINC,TIME,NOEL,NPT,COORDS,
     1 JLTYP,TEMP,PRESS,SNAME)
C
      INCLUDE 'ABA_PARAM.INC'
C
      DIMENSION FLUX(2), TIME(2), COORDS(3)
      CHARACTER*80 SNAME
C
C     user coding to define FLUX(1) and FLUX(2)
C
      RETURN
      END

关键参数含义（摘自官方文档）：

  • FLUX(1)（输出）：流入模型的热流密度。体热流单位为 JT⁻¹L⁻³（即 W/mm³）。
  • FLUX(2)（输出）：dq/dθ，热流对温度的导数，用于改善非线性收敛；热流不随温度变化时设为 0。
  • COORDS(3)（输入）：当前积分点坐标，是计算积分点与热源中心相对位置的关键。
  • TIME(1)（输入）：当前分析步时间，用于计算热源中心的瞬时位置。
  • JLTYP（输入）：热流类型标识，体热流 BFNU 对应 JLTYP=1。

此外，根据非耦合热传递分析文档（路径：SIMACAEANLRefMap/simaanl-c-heattransfer.htm），热传递分析可包含热传导、边界对流与边界辐射，可为瞬态或稳态、线性或非线性。焊接仿真的非线性主要来源于材料属性随温度变化、潜热效应以及辐射边界（温度越高非线性越强）。

热传递分析流程

四、实战要点：多年经验总结的坑

焊接仿真看似流程清晰，实际操作中有大量细节决定成败。以下是多年实战中反复验证的要点：

4.1 网格与时间增量的匹配

移动热源每步前进的距离必须小于热源特征尺寸，否则温度场会出现”跳跃”，导致结果失真。经验准则：单步热源移动距离 ≤ 最小网格尺寸的 1/2~1/3。例如焊接速度 v=5 mm/s、网格尺寸 1 mm，则时间增量 Δt 应取 0.05~0.1 s 量级。焊缝附近网格要足够细（0.5~1 mm），远离焊缝处可逐步过渡粗化以节省计算量。

4.2 温度相关材料必须覆盖全温度范围

材料属性（热导率、比热、密度、弹性模量、屈服应力、热膨胀系数）必须从室温一直定义到熔点以上。高温区材料强度急剧下降，接近熔点时屈服应力应趋近于一个很小的值（如 10 MPa），否则熔池区域会产生不真实的应力。比热容在相变温度附近可设置一个尖峰来近似考虑潜热。

4.3 收敛控制

热分析中用 DELTMX 限制单步最大温度增量（如 50℃），防止热源刚进入时温度骤升导致发散。应力分析建议开启 NLGEOM=YES 考虑大变形，并使用较小的初始增量步配合自动增量。残余应力分析中刚体位移必须充分约束，但约束方式不能引入附加应力（推荐对称边界或最小约束）。

4.4 单元选择

热分析用一阶热传递单元 DC3D8，应力分析用一阶实体单元 C3D8（或 C3D8R 减缩积分）。两套分析建议使用相同网格，以便温度场通过预定义场准确映射。

五、Abaqus 实现：DFLUX 子程序与 INP 关键字

5.1 Goldak 双椭球 DFLUX 完整实现

下面是基于官方接口实现的 Goldak 双椭球移动热源子程序，焊接方向沿 X 轴正向：

      SUBROUTINE DFLUX(FLUX,SOL,KSTEP,KINC,TIME,NOEL,NPT,COORDS,
     1 JLTYP,TEMP,PRESS,SNAME)
C
      INCLUDE 'ABA_PARAM.INC'
C
      DIMENSION FLUX(2), TIME(2), COORDS(3)
      CHARACTER*80 SNAME
C
      REAL*8 Q, ETA, V, A, B, CF, CR
      REAL*8 X0, Y0, Z0, XC, YC, ZC
      REAL*8 FF, FR, PI
      REAL*8 DX, DY, DZ, R2
C
      PARAMETER (PI=3.14159265358979D0)
C
C     === 焊接工艺参数（按实际工艺修改）===
      Q    = 3000.0D0        ! 热输入(W) = 电压*电流
      ETA  = 0.8D0           ! 热效率(MIG焊约0.8)
      V    = 5.0D0           ! 焊接速度(mm/s)
      A    = 5.0D0           ! 焊接方向半长轴(mm)
      B    = 4.0D0           ! 椭球半宽(mm)
      CF   = 1.0D0           ! 前半椭球深度参数
      CR   = 2.0D0           ! 后半椭球深度参数
C
      X0   = 0.0D0           ! 起始点坐标
      Y0   = 0.0D0
      Z0   = 0.0D0
C
C     === 当前时刻热源中心位置 ===
      XC = X0 + V * TIME(1)
      YC = Y0
      ZC = Z0
C
C     === 积分点相对热源中心的距离 ===
      DX = COORDS(1) - XC
      DY = COORDS(2) - YC
      DZ = COORDS(3) - ZC
C
C     === Goldak 双椭球热流密度 ===
      IF (DX .GE. 0.0D0) THEN
C         前半椭球（热源前方）
          FF = 2.0D0 * CF / (CF + CR)
          R2 = (DX/A)**2 + (DY/B)**2 + (DZ/CF)**2
          FLUX(1) = FF * 6.0D0*SQRT(3.0D0) * Q * ETA
     1              / (A*B*CF * PI*SQRT(PI))
     2              * EXP(-3.0D0 * R2)
      ELSE
C         后半椭球（热源后方）
          FR = 2.0D0 * CR / (CF + CR)
          R2 = (DX/A)**2 + (DY/B)**2 + (DZ/CR)**2
          FLUX(1) = FR * 6.0D0*SQRT(3.0D0) * Q * ETA
     1              / (A*B*CR * PI*SQRT(PI))
     2              * EXP(-3.0D0 * R2)
      END IF
C
      FLUX(2) = 0.0D0        ! dq/dT = 0
C
      RETURN
      END

代码要点：通过 TIME(1) 计算热源中心瞬时位置 XC，再用 COORDS 求积分点相对位置；根据 DX 正负切换前/后椭球公式；FLUX(2)=0 表示热流不随温度变化。

5.2 热传递分析 INP 关键字

** 热传递分析步
*STEP, NAME=HeatTransfer, NLGEOM=NO, INC=10000
*HEAT TRANSFER, TRANSIENT, DELTMX=50.
0.01, 10., 1e-5, 0.01
**
** 施加体热流（调用DFLUX子程序）
*DFLUX
WELD_SET, BFNU, 1.0
**
** 对流边界条件
*FILM
SURF_ALL, F1, 25., 1e-5
**
** 辐射边界条件
*RADIATE
SURF_ALL, R1, 25., 0.8
**
** 温度输出请求
*NODE OUTPUT
NT
*END STEP

说明：WELD_SET 是施加体热流的单元集；BFNU 表示非均匀体热流，会触发 DFLUX 子程序（JLTYP=1）；DELTMX=50. 限制最大温度增量防止发散；时间增量 0.01 s 需根据焊接速度与网格尺寸调整。

5.3 顺序耦合应力分析 INP 关键字

** 应力分析步
*STEP, NAME=Stress, NLGEOM=YES, INC=10000
*STATIC
0.01, 10., 1e-6, 0.1
**
** 导入热传递分析的温度场
*PREDEFINED FIELD, TYPE=TEMPERATURE, FILE=heat_transfer_job
**
** 温度相关弹性属性
*ELASTIC
210000., 0.3, 20.
190000., 0.3, 200.
150000., 0.3, 600.
50000.,  0.3, 1000.
10000.,  0.3, 1500.
**
*EXPANSION
1.2e-5, 20.
1.3e-5, 200.
1.5e-5, 600.
2.0e-5, 1000.
**
*PLASTIC
250., 0.0, 20.
200., 0.0, 200.
100., 0.0, 600.
50.,  0.0, 1000.
10.,  0.0, 1500.
**
** 输出请求
*NODE OUTPUT
U, RF
*ELEMENT OUTPUT
S, E, PEEQ
*END STEP

说明：FILE=heat_transfer_job 对应热传递分析的 job 名；材料属性必须覆盖室温到熔点以上的完整范围；高温时屈服应力趋近于零是物理真实，不能省略；建议 NLGEOM=YES 考虑大变形。

六、结果验证

仿真结果的可信度需要通过验证来保证。常见的验证手段包括：

  • 温度场验证：将仿真得到的熔池形状、热影响区宽度与金相试验对比；用热电偶实测温度-时间曲线与仿真节点温度历史对比。
  • 残余应力验证：采用盲孔法、X 射线衍射法或切条法测量表面/内部残余应力，与仿真结果对比分布趋势与峰值。
  • 变形验证：测量焊后工件的角变形、纵向/横向收缩量，与仿真位移结果对比。

典型的焊接温度场仿真结果如下图所示，熔池区域达到峰值温度，沿焊接方向呈现明显的温度梯度，热源后方为逐渐冷却的尾迹。

焊接温度场仿真云图

残余应力分布的一般规律是：焊缝及近缝区纵向残余拉应力接近材料屈服强度，远离焊缝逐渐衰减并转为压应力，整体满足自平衡条件。这一规律与大量文献报道及工程实测一致。

七、总结

焊接残余应力仿真的核心可归纳为四句话：方法上用顺序热力耦合，热源上用 Goldak 双椭球 DFLUX，材料上覆盖全温度范围，网格时间增量上保证热源平稳移动。掌握这条主线，再结合对流/辐射边界、收敛控制和结果验证，就能搭建起可靠的焊接仿真流程。

焊接仿真涉及工艺参数众多、非线性强、调试周期长，初学者常在热源参数标定、收敛发散、材料数据缺失等环节卡壳。建议从简单的单道平板堆焊入手，先跑通温度场并验证熔池形状，再叠加应力分析，逐步扩展到多道焊与复杂结构。

参考文献

  1. Goldak J., Chakravarti A., Bibby M. A new finite element model for welding heat sources. Metallurgical Transactions B, 1984, 15(2): 299-305.
  2. Dassault Systèmes. Abaqus User Subroutines Reference Guide — DFLUX (SIMACAESUBRefMap/simasub-c-dflux.htm).
  3. Dassault Systèmes. Abaqus Analysis User’s Guide — Uncoupled heat transfer analysis (SIMACAEANLRefMap/simaanl-c-heattransfer.htm).