用MATLAB求解钢在渗硼过程中的渗层厚度

用MATLAB求解钢在渗硼过程中的渗层厚度

摘 要：本文介绍了如何用MATLAB中的PDE工具箱求解钢在渗硼过程中的渗层厚度。利用PDE工具箱，不需编程，直接进入用户图形界面对实际问题求解。对渗层厚度的计算表明，PDE工具箱求解问题快捷、简单、通用性强。 关键词：MATLAB；渗层厚度；渗硼；PDE

1.概述

在材料热处理生产实际中，通常需要确定渗层的厚度、渗入物的浓度分布，因而研究特殊非稳态扩散尤为重要。非稳态扩散问题的控制方程一般为多维非线性方程，不能直接计算出解析解。一般都是将时间、空间坐标划分为许多的网格，然后借助于计算机编程求解。而利用MATLAB中的PDE工具箱，无需编程则可直接对特殊边界条件的非稳态导热问题进行求解。

2.PDE工具箱简介

MATLAB中的偏微分方程(PDE)工具箱是用有限元法求解偏微分方程得到数值近似解，可以求解线性的椭圆型、抛物线型、双曲线型偏微分方程及本征型方程和简单的非线性偏微分方程，具体可处理数学模型形式如下： 椭圆型PDE

非线性PDE

本征型问题

抛物线型PDE

双曲线型PDE

式中： u－域Ω上求解变量

t－时间变量 ε－特征值

c、a、f常数或者变量

???(c?u)?au?f (1)

???(c(u)?u)?a(u)u?f(u) (2)

???(c?u)?au???du

(3)

d?u???(c?)?au?f ?t(4)

?2ud2???(c?)?au?f ?t(5)

1

上述几种数学模型多在热传导（扩散）、电磁学和声学的波传导等问题求解中应用。

3、扩散问题的MATLAB解法

对于大多数的扩散问题，求解浓度场时很难得到解析解，只能利用计算机得到数值解来无限接近精确解。数值解方法又分有限元法、有限差分法、混合微分差分法、离散元法、拉格朗日元法等，其中有限元法是利用部分插值把区域连续求解的微分方程离散成求解线性代数方程组。在使用MATLAB的PDE工具箱进行有限元计算前需要有一些预处理的工作，如对所求解模型的几何形状或者形体进行离散化，即用比较简单的形状和形体来逼近和代替实际的形状和形体，这样可以把比较复杂的曲线和曲面问题转化为相对简单的直线或平面问题。

在实际求解扩散问题时，主要有以下五步。第一，建立一个用来描述对应浓度问题的物理模型，第二，根据需要对求解问题赋予边界条件，MATLAB指定了如下3种边界条件：①Dirichlet 条件，hu=r; ②Neumann条件，n?(c?u)?qu?g; ③混合边界条件，Dirichlet条件和Neumann条件的组合。式中n为垂直于边界的单位矢量，h、r、q、g为常量或与u有关的变量。第三，确定偏微分方程的类型，并结合已知条件设定方程参数。第四，创建初始三角形网格以及细化网格。第五，设定求解参数并求解偏微分方程。

4．用PDE工具箱求解浓度场

浓度场中用来描述三维非稳态扩散微分方程的一般形式为：

?C?D?2C ?t(6)

其中，D为扩散系数，t为时间。要通过控制方程(6)获得某一具体扩散问题的浓度分布，

需要结合定解条件（初始条件、边界条件）来求解，通过下面的扩散模型来体现。 对一10cm×10cm×5cm的碳含量为ωC＝0.8长方体钢材进行渗硼，渗硼温度为850℃。已知850℃时硼在T8中的扩散系数D=1.10×10-13 m2/s，分析渗层厚度随时间（0～9000s）变化的情况。

2

图 1 Fe-B相图

根据图1可知：在1000℃以下渗硼时，随着硼浓度的升高形成铁硼化合物，当含硼量达到8.84％左右时，形成稳定的中间化合物Fe2B；当含硼量达到16.23％左右时，形成含硼量更高的稳定化合物FeB。

硼在铁中的扩散属于反应扩散。只有当表面形成高含硼量的稳定化合物FeB时，硼才会向内部扩散形成FeB或Fe2B。因而假定开始处理时材料表面即达到wB＝16.23％，wB＝8.84％的面距表面的距离即为渗层厚度。

首先建立物理模型。根据实验观察，850℃时，9000s时的渗硼厚度达不到100μm。相对于长宽高的数值，可将模型简化为长宽比为1：10的长方形。设宽为100μm，高为10μm。单击

工具，在窗口拉出一个矩形，双击矩形区域，在Object Dialog对话框输入

Left为0，Bottom为0，Width为1e-4,Height为1e-5。由于数值过小（默认值为1左右），图形不可见，所以要调整坐标显示比例。方法：选择Options->Axes Limits,把X，Y轴的自动选项打开，并打开axes equal。

物理模型如图2所示：

图 2 扩散的物理模型

等价的偏微分方程组如式（7）：

3

??C2?D???t??C(x,y,0)?0?C?16.23%,左边界（y轴）??C?0,右边界(x=0.0001) ???C?0,下边界(x轴)??n??C??0,上边界(y=0.00001)??n

(7)

式（7）中，t为时间，C为t时刻时（x，y）处硼的浓度，D为扩散系数，D＝1.10×10

-13

m/s。n为垂直于边界的单位矢量。

然后，设定边界条件。此种条件下，硼在Fe中主要沿X单向扩散，在Y向可认为没有

，使边界变红色，然后分别双击每段边界，打开Boundary Conditions对

2

扩散。单击

话框，设置边界条件。上下边界：由于上下边界与外界绝缘，因而选择Neumann条件，令g为0，q为0。左边界：保持定值16.23％，所以选择Dirichlet条件，令h=1,r=0.1623。右边界：保持定值0，所以选择Dirichlet条件，令h=1,r=0。

接着，确定方程类型并输入相关参数。由式（7）可以看出，方程属于Parabolic（抛物型）。单击

，打开PDE Specification对话框，输入c=1.1e-13，a=0,f=0,d=1。

，网格个数设为5000个左右。

其次创建网格，单击

最后设定求解参数并解方程。单击Solve菜单中Parameters…选项，打开Solve Parameters对话框，输入Time为0:60:9000（每分钟一个观察点），u(t0)为0,其他不变。完成后，单击

，开始解方程。

5.分析处理数据

MATLAB只能显示出浓度的分布情况，不能直接求出渗层厚度，因而有必要对MATLAB计算的结果进一步处理。

点击Mesh->Export Mesh输出网格数据，点击Solve->Export Solution输出按网格序号排列的数值解u。通过以下语句可将u整理成按距表面距离（即x）排列。

u1=[1e6*p']; u1=u1(:,1); u1=[u1,100*u];

图3是根据MATLAB计算的渗硼进行一小时后硼浓度的分布曲线。根据图1可知，渗层中的硼浓度并不是连续变化的，图4是根据Fe-B相图修正出的图像。

%将网格坐标的单位由m化为μm %只取第一列（即x）

%将浓度化为百分数并与x值组合成新矩阵。 %按x升序排列

u1=sortrows(u1,1);

4