用有限差分法和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时的三维图：。