用有限差分法和Matlab计算二维热加工温度场分析

用有限差分法和Matlab 计算二维热加工温度场分析

Eg1 薄板焊接过程温度场分析。

取焊件的一半为模型进行离散化，起始点为o 点，以后以v 速度沿x 轴运动。

根据题意为二维不稳态导热，二维不稳态导热方程为：

k Q

y T x T T 12222-

+∂∂+∂∂=∂∂τα

y

x 左

边

界

下边界

图 二维焊接离散化

题目可化为以下微分方程组：(以y 轴正方向为上，x 轴正方向为右)

(,,0)0,(),(),(),e e e T

C k Q

T x y T

T

k x

T k T T x T

k T T y T

k T T x ρτβββ∂⎧=∆+⎪∂⎪=⎪⎪∂=⎪∂⎪⎪∂⎨=-⎪∂⎪∂⎪=-⎪∂⎪∂⎪=-⎪∂⎩左边界（y 轴）

右边界

下边界(x 轴)

上边界

需要的参数（均已cm,cal,g 为单位，所以不必换算）：

参数

数值 备注 ρ

7.82 g/cm 3 密度 v

0.4 cm/s 焊接速度 h

1cm 板厚度 Qm 4000 cal/cm 3 热源分布密度 β

0.0008cal/cm 2·s ·℃ 换热系数 Te ，T 0

20℃ 周边介质温度，初始温度 k 0.1 cal/c m ·s ·℃ 导热系数

用PDE Tool 解题步骤如下：

1． 区域设置

单击工具，在窗口拉出一个矩形，双击矩形区域，在Object Dialog 对话框输入Left 为0，Bottom 为0，Width 为2,Height 为2。

与默认的坐标相比，图形小的看不见，所以要调整坐标显示比例。方法：选择Options->Axes Limits,把X ，Y 轴的自动选项打开。

设置Options->Application 为Heat Transfer （设置程序应用热传输模型）

2． 设置边界条件

单击

，使边界变红色，然后分别双击每段边界，打开Boundary Conditions 对话框，设置边界条件（根据边界条件）。所有的边界都为Neumann 条件。输入值见下表：

边界

g (Heat Flux) q (Heat Transfer Coefficient) 左边界

0 0 右边界

0.0008*20 0.0008 下边界（x 轴）

-0.0008*20 -0.0008 上边界

0.0008*20 0.0008

3． 设置方程类型

单击，打开PDE Specification 对话框，设置方程类型为Parabolic （抛物型）， rho(密度）为7.82，C （比热）为0.16，k （导热系数）为0.1，Q （热源）为4000*exp(-3*(x.^2+(y-0.4*t).^2)/0.49)，其他参数为0

4． 网格划分

单击，或者加密网格，单击。

5． 初值和误差的设置

222

2

222340000exp(3())()/exp(3()13(())/exp(.)0.49

44000m r Q h r

x y cal cm cm r

x y cal c Q m ττν=•-+-=•--+-=•则：

单击Solve菜单中Parameters…选项，打开Solve Parameters对话框，输入Time为0:0.5:5，u(t0)为20,其他不变。

6．解方程

单击，开始解方程。

7．整理数据

单击Mesh->Export Mesh…输出p e t 的数值，单击Solve->Export Solution…输出u.

回到Matlab主窗口执行下面两条命令：

u1=[p',u(:,7)] %将节点坐标和其在3s时的温度组成新矩阵

u2=sortrows(u1,3) %将u1按温度值大小升序排列。

u1=[p',u(:,4)] %将节点坐标和其在1.5s时的温度组成新矩阵

u2=sortrows(u1,3) %将u1按温度值大小升序排列。

8．温度场分布：

1.5s时

1.5s时的三维图

3s 时的温度场：

3s时的三维图：