基于有限差分法的大线能量焊接热模拟试件的温度分布

李矗东,2，魏 强，李玉中，黄宏虎

(1.新兴铸管股份有限公司，北京 100026；2.北京科技大学 东凌经济管理学院，北京 100083)

近年来，造船业不断发展，传统的船板钢已不能满足实际工业需求，在焊接方面主要表现为传统的小线能量焊接施工效率低，而大线能量焊接用钢性能则难以保证[1-3]。其中钢板焊接烧穿问题是影响线能量提高的重要原因。烧穿的本质是未熔化的金属不能承受其所受到的应力。例如，在管线钢服役期，维护工作不可避免。在焊接进行时，管道内部存在流体会产生压力，而焊接过程的高温会降低管道的承压能力，当承压能力低于压力时，烧穿问题不可避免。就大线能量焊接船板钢而言，虽然不存在流体压力的影响，但是线能量的增加会降低钢板的承压能力，导致钢板焊穿。Akbari等[4]研究了脉冲激光焊接条件下的熔池温度分布温度，建立了三维导热微分方程，预测了熔池的深度和宽度。Yadaiah和Bag[5]研究了焊接过程的热源问题，利用有限元法建立了氩弧焊的热源模型，并说明了热传递过程中的温度分布情况。然而，焊接热模拟技术是改善焊接工艺的重要研究手段，文中将从温度的角度，分析大线能量焊接热模拟试件的温度分布问题，为焊接工艺提供一定的理论基础。

1 数值模型建立

本模型以焊接热模拟为例，探究线能量为100 kJ·cm-1的热模拟试件热量传递问题。温度场分布可分为三个阶段进行描述，一是升温阶段，热量的总和增加。二是保温阶段，三是冷却阶段，在以下模型中主要讨论试件的升温以及保温阶段。

1.1 焊接热模拟参数

针对试验钢(厚度为25 mm)，进行线能量为100 kJ·cm-1的焊接热模拟实验。将试验钢垂直于轧向取样，加工成Φ6 mm×80 mm试样，在Gleeble3500试验机上进行焊接热模拟[6]。从室温以150 ℃/s的加热速度将试样加热至峰值温度,高温保温3 s，再以一定的速度进行冷却。设定热模拟参数，根据软件所生成的不同t8/5程序，绘制焊接热循环曲线如图1所示。热模拟的各项参数如表1所示。

图1 焊接热循环曲线

表1 热模拟程序参数

用砂纸打磨试样表面直至光亮，避免试样在进行热模拟的过程中发生断裂。用游标卡尺找到试样的中点，并做好记号。在试样中点处点焊一对K型热电偶(用于测温，测温范围在8 ～18 mm)，并保证热电偶的牢固性。将焊有热电偶的试样安装在相应规格的卡具上，保证装配的稳固性。在工作室抽成真空状态之后，填充高纯氩气作为保护气体。利用Gleeble3500试验机专用的编程软件进行程序的编制，然后运行。

1.2 焊接热模拟传热模型

在焊接过程中，试样从室温25 ℃加热1 350 ℃并保温3 s的过程中，试样处于吸热状态。在冷却过程中，试样整体属于放热状态，整个过程中不考虑辐射传热。

根据能量守恒原理以及傅里叶定律，对试件进行微分处理，可建立导热微分方程[7]。

(1)

式中：ρ为试件的密度，7.8 g·cm-3；c为试件金属的比热容，0.7 J·g-1·oC-1；λ为试件金属的导热率，0.5 J·cm-1·oC-1；qv为单位时间，单位体积的物体生成的热量，J·cm-3·s-1。

将整个过程分为2个时间段考虑:①从室温以150 ℃·s-1的加热速度将试样加热至峰值温度(1 350 ℃);②在峰值温度高温保温3 s。根据式(1)建立控制方程和定解条件，在各边界上设定第一类边界条件。然后进行空间离散化，最后建立代数方程。

2 数值模型求解

2.1 推导有限差分法传热方程

在有限元分析中，单元内的任意点(x,y,z)的场变量ψ(x,y,z) 需要通过选定插值形式来确定。根据单元节点值进行插值的求解。以节点温度作为参考量，表示单元内的某一点温度，由此构建插值函数[8-11]。

将温度函数看作三维方向的传导函数，即温度仅在(x,y,z) 三个坐标轴方向传导，即单元中某一点 (x,y,z)的温度函数为

t(x,y,z)=(tx,ty,tz)

(2)

对于空间内的任意点而言，热量可双向传递，但在两个传递方向上，热量的传递情况不同，如图2所示。

图2 热量传递示意

(3)

同理

(4)

(5)

为获得离散化的时间和空间，将空间按立方体分割，取试件上的某一微元体，如图3所示。

图3 空间离散示意

设定网格元的某个顶点坐标为(x,y,z)，可记为

(x,y,z)=(iΔx,jΔy,kΔz)

(6)

式中：Δx,Δy,Δz分别为x轴，y轴，z轴上的步长，i，j，k为整数。在时间上，取时间步长为Δτ，则n时刻为nΔτ，据此划分原则，任意一个空间和时间上的温度函数可表示为

tn(i,j,k)=tn(iΔx,jΔy,kΔz)

(7)

根据微分学原理以及一维稳态传热可知，当试件内的温度是一个连续函数，那么在某一点b处的导数可用极限表示[12-16]。

(8)

根据式(8)建立三维非稳态导热的差分方程，利用向前差商近似代替一阶偏导数式(9)，利用中心差商近似代替二阶偏导数式(10)。

(9)

(10)

将式(9)和式(10)带入式(1)中，可得式(11)。

(11)

若采用正方体网格，可得到式(12)方程。

(12)

2.2 边界条件设立

热模拟的加热原理是依靠电阻自身加热，因此，假设热源在整个试件的中心。加热速度为150 ℃·s-1，在0～8.8 s的时间范围内，加热温度T=150τ+25。设定热电偶的测定范围为8 mm，在保温阶段，热电偶检测范围内的温度均达到1 350 ℃。

在物体的传热问题上存在着三种边界条件，具体如下[17]：

第一类边界条件：物体的边界温度随时间的变化规律一定。

t=f(x,y,z,τ)

(13)

在较为简单的传热过程中，物体的边界温度可为常数。

第二类边界条件：在物体的边界上，各个位置的热通量与时间的函数关系一定。

(14)

式中：qw可为常数或是时间函数。当物体边界处于绝热状态时，qw=0。

第三类边界条件：物体周围介质的温度tf和物体与边界之间的对流传热系数h一定。

(15)

式中：Δt为物体温度与物体周围介质的温度tf的差值,℃。

根据热模拟试验机的加热原理，假设以下边界条件：

(1)在加热升温过程中，假设试件中心3 mm厚的单元体为发热中心，如图4所示。

在加热过程中，以yOz面为界，试件左侧与试件右

图4 试样热影响区示意

侧的传热规律相同。因此以试件左半部分为例进行分析，在升温过程(即0～8.8 s)中，热量变化Q服从式(16)。

Q=cmΔt=cρV(150τ)

(16)

式中：V为试件左半部分的加热单元体体积。

(2)假设试件其余位置的温度为室温25 ℃，不考虑辐射传热。由于工作室抽成真空状态之后，填充高纯氩气作为保护气体，在此条件下，可忽略对流传热。

2.3 偏微分方程的Matlab解法

在寻找某一物体的传热规律时，通常需要求解相应物理条件下所满足的传热方程，通常所遇到的是椭圆型偏微分方程(17)和抛物型偏微分方程(18)[17]。

-▽(cu)+au=f,inΩ

(17)

(18)

式中：Ω为有界区域，c；a,f以及未知函数u是定义在Ω上的函数；d是定义在Ω上的复函数。在抛物型方程中，c,a,f,d均可为时间的函数。

利用Matlab中的PDE工具箱求解偏微分方程，主要存在两个边界条件，一是Dirichlet条件，二是Neumann条件[18-21]。

Dirichlet条件：

hu=t,in∂Ω

(19)

Neumann条件：

(20)

根据式(19)与式(20)中,g,q和h是定义在∂Ω上的函数。对抛物型方程而言，系数g,q和h可以是关于时间t的函数。

3 结果与分析

由于加热区域为3 mm厚的单元体如图3所示。因此，在xOz平面上传热分布有较大区别，观察面定义在xOz平面。试件尺寸为10 mm × 10 mm × 60 mm,利用Matlab中的PDE工具箱对偏微分方程进行求解，求解结果如图5所示。

图5 温度分布

图5(a)和图5(b)分别为升温1 s和2 s时的温度分布图。图5(c)和图5(d)分别为升温3 s和4 s时的温度分布图。图5(e)和图5(f)分别为升温9 s，以及升温9 s且保温3 s时的温度分布图。通过对比发现，在升温2 s后的等温区域比升温1 s后的等温区域小，中心温度也比较高。在升温过程中，加热单元体处于升温状态，因此，中心区域(即观察图左侧)的温度不断升高。在加热初期，由于传热规律不稳定，在等温区域上存在较大的差别。在升温3、4 s后，虽然中心温度及边缘温度均处于升高阶段，但传热规律稳定，等温区域面积基本相同。在4～8 s时间段内，传热规律大致相同。在9 s时，升温停止，处于保温状态，与3～8 s时间范围内的传热规律相比，等温区域稍有扩大。在12 s时，保温结束，等温区域减小，但仍比3～8 s时的等温区域大。

4 结 论

(1)基于不同时间、不同边界条件下的传热规律，建立了大线能量焊接热模拟试件温度分布的数学模型，得出了温度分布图。

(2)在焊接热模拟升温初期，传热规律变化较大，升温中期比较稳定，在升温后期以及保温阶段，传热规律有轻微变化。