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

用MATLAB求解钢在渗硼过程中的渗层厚度摘 要：本文介绍了如何用MATLAB中的PDE工具箱求解钢在渗硼过程中的渗层厚度。利用PDE工具箱，不需编程，直接进入用户图形界面对实际问题求解。对渗层厚度的计算表明，PDE工具箱求解问题快捷、简单、通用性强。关键词：MATLAB；渗层厚度；渗硼；PDE1.概述在材料热处理生产实际中，通常需要确定渗层的厚度、渗入物的浓度分布，因而研究特殊非稳态扩散尤为重要。非稳态扩散问题的控制方程一般为多维非线性方程，不能直接计算出解析解。一般都是将时间、空间坐标划分为许多的网格，然后借助于计算机编程求解。而利用MATLAB中的PDE工具箱，无需编程则可直接对特殊边界条件的非稳态导热问题进行求解。2.PDE工具箱简介MATLAB中的偏微分方程(PDE)工具箱是用有限元法求解偏微分方程得到数值近似解，可以求解线性的椭圆型、抛物线型、双曲线型偏微分方程及本征型方程和简单的非线性偏微分方程，具体可处理数学模型形式如下：椭圆型PDE(1)非线性PDE(2)本征型问题(3)抛物线型PDE(4)双曲线型PDE(5)式中： u域上求解变量t时间变量特征值c、a、f常数或者变量上述几种数学模型多在热传导（扩散）、电磁学和声学的波传导等问题求解中应用。3、扩散问题的MATLAB解法对于大多数的扩散问题，求解浓度场时很难得到解析解，只能利用计算机得到数值解来无限接近精确解。数值解方法又分有限元法、有限差分法、混合微分差分法、离散元法、拉格朗日元法等，其中有限元法是利用部分插值把区域连续求解的微分方程离散成求解线性代数方程组。在使用MATLAB的PDE工具箱进行有限元计算前需要有一些预处理的工作，如对所求解模型的几何形状或者形体进行离散化，即用比较简单的形状和形体来逼近和代替实际的形状和形体，这样可以把比较复杂的曲线和曲面问题转化为相对简单的直线或平面问题。在实际求解扩散问题时，主要有以下五步。第一，建立一个用来描述对应浓度问题的物理模型，第二，根据需要对求解问题赋予边界条件，MATLAB指定了如下3种边界条件：Dirichlet 条件，hu=r; Neumann条件，; 混合边界条件，Dirichlet条件和Neumann条件的组合。式中n为垂直于边界的单位矢量，h、r、q、g为常量或与u有关的变量。第三，确定偏微分方程的类型，并结合已知条件设定方程参数。第四，创建初始三角形网格以及细化网格。第五，设定求解参数并求解偏微分方程。4用PDE工具箱求解浓度场浓度场中用来描述三维非稳态扩散微分方程的一般形式为：(6)其中，D为扩散系数，t为时间。要通过控制方程(6)获得某一具体扩散问题的浓度分布，需要结合定解条件（初始条件、边界条件）来求解，通过下面的扩散模型来体现。对一10cm×10cm×5cm的碳含量为C0.8长方体钢材进行渗硼，渗硼温度为850。已知850时硼在T8中的扩散系数D=1.10×10-13 m2/s，分析渗层厚度随时间（09000s）变化的情况。图 1 Fe-B相图根据图1可知：在1000以下渗硼时，随着硼浓度的升高形成铁硼化合物，当含硼量达到8.84左右时，形成稳定的中间化合物Fe2B；当含硼量达到16.23左右时，形成含硼量更高的稳定化合物FeB。硼在铁中的扩散属于反应扩散。只有当表面形成高含硼量的稳定化合物FeB时，硼才会向内部扩散形成FeB或Fe2B。因而假定开始处理时材料表面即达到wB16.23，wB8.84的面距表面的距离即为渗层厚度。首先建立物理模型。根据实验观察，850时，9000s时的渗硼厚度达不到100m。相对于长宽高的数值，可将模型简化为长宽比为1：10的长方形。设宽为100m，高为10m。单击工具，在窗口拉出一个矩形，双击矩形区域，在Object Dialog对话框输入Left为0，Bottom为0，Width为1e-4,Height为1e-5。由于数值过小（默认值为1左右），图形不可见，所以要调整坐标显示比例。方法：选择Options->Axes Limits,把X，Y轴的自动选项打开，并打开axes equal。物理模型如图2所示：图 2 扩散的物理模型等价的偏微分方程组如式（7）：(7)式（7）中，t为时间，C为t时刻时（x，y）处硼的浓度，D为扩散系数，D1.10×10-13 m2/s。n为垂直于边界的单位矢量。然后，设定边界条件。此种条件下，硼在Fe中主要沿X单向扩散，在Y向可认为没有扩散。单击 ，使边界变红色，然后分别双击每段边界，打开Boundary Conditions对话框，设置边界条件。上下边界：由于上下边界与外界绝缘，因而选择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'%将网格坐标的单位由m化为mu1=u1(:,1);%只取第一列（即x）u1=u1,100*u;%将浓度化为百分数并与x值组合成新矩阵。u1=sortrows(u1,1);%按x升序排列图3是根据MATLAB计算的渗硼进行一小时后硼浓度的分布曲线。根据图1可知，渗层中的硼浓度并不是连续变化的，图4是根据Fe-B相图修正出的图像。图 3 MATLAB计算的渗层中硼浓度渗层深度的变化曲线（850×1h）图 4根据Fe-B相图修正的渗层硼浓度渗层深度的变化曲线（850×1h）为了反映渗层随时间变化的趋势，需求出各个时刻的渗层厚度。如下：rows cols=size(u1); %求出u1的行列数a=0; %临时存储某一列最接近8.48的x值c=0; %保存每分钟的渗层厚度for j=2:cols b=1; for i=1:rows%找出与8.84最接近的x值 if abs(u1(i,j)-8.84)<b a=u1(i,1); b=abs(u1(i,j)-8.84); end end c(j-1)=a;endc=c'c=(0:60:9000)',c;%将渗层厚度加上时间标志图5是根据上述语句计算出的渗层厚度随时间变化的曲线图。从图5可看出，初始阶段扩散速度很快，扩散层的厚度与时间呈直线关系，随扩散时间延长，反应层增厚，扩散层厚度与时间呈抛物线关系，这非常符合反映扩散的特点1，且计算出的渗层厚度与实验测得的数据（表1）也非常吻和。图 5渗层厚度随时间变化的曲线表 1 实验测得的渗层厚度（850）时间（s）实验测得厚度值(m)计算厚度值(m)相对误差36001.60E-051.67E-054.05%54002.00E-052.03E-051.23%72002.50E-052.44E-052.56%90002.75E-052.69E-052.33%6.结论研究浓度场的分布在工程领域中有重要的作用。本文中提及的扩散模型对固体材料中的原子扩散计算具有重要的实际意义。本模型不仅适用于渗硼，对于渗碳、渗氮等问题同样适用。此外，可根据计算来确定达到规定厚度所需的时间，对热处理具有很好的指导意义。从上述使用PDE工具箱求解扩散模型的结果可以看出使用MATLAB的快捷灵活。若自己编程，过程复杂麻烦，当步长选取不当会导致求解结果发散或振荡，因而结合PDE工具箱直接调用其中网格划分。求解方程的命令来求解浓度场可以快速得到理想的结果。在GUI上可以处理复杂几何形状的扩散问题,如多种物性参数材料的叠加或嵌入，这是一般的程序很难用语言描述的，有了网格的精化,使得有限元数值解的精度大大提高，MATLAB使用灵活,易于修改,这是其它语言无法比拟的。因而可见，使用MATLAB的PDE工具箱可以方便快速地解出特殊几何形状和复杂边界条件的导热问题，多种可视化图形使得结果便于理解，也为工程上一些实际扩散问题提供了解决方法。参考文献1 王掌忠.材料科学基础材.北京:机械工业出版社,2005,121-141.2 苏金明,张莲花,刘波.MATLAB工具箱应用M.北京:电子工业出版社,2004.3 何仁斌.MATLAB6工程计算及应用M.重庆:重庆大学出版社,2001.4 苏金明,王永利.MATLAB7.0实用指南M.北京:电子工业出版社,2005.4 陆君安.偏微分方程的MATLAB解法M.武汉:武汉大学出版社,2001.渗硼!渗铬和渗钒都是反应扩散,反应扩散的过程是渗入元素首先溶于被渗金属中,当它超过溶解度时就形成新相层2,5"所谓渗硼!渗铬和渗钒深度都是指这些新相层厚度"与新相层相邻的为过渡区,即扩散元素(渗入元素)在C2Fe中的固溶区"如果试样预先渗碳或N2C共渗,在过渡区内将同时存在碳!氮和渗入元素的浓度梯度"在高温下,碳!氮原子向心部方向扩散的瞬间,将会把与其邻接的铁原子推挤至高能不稳定状态,这无疑有利于铬或钒原子和高能态的铁原子交换位置,使铬或钒原子向心部方向进行置换扩散变得容易,也就是说能与奥氏体形成间隙固溶体且扩散系数较大的碳或氮的扩散,加快了铬和钒在过渡区的扩散速度"硼在C2Fe中既以间隙固溶体存在,也以置换固溶体存在,这是目前已知同一元素可以与铁形成两类固溶体的唯一元素14"但与碳和氮相比,硼的原子半径较大14,即使是进行间隙扩散,其扩散激活能也较碳氮扩散激活能高,所以过渡区碳!氮向心部方向的扩散,同样有利于硼元素的扩散"为固体硼法1.2,然而,固体渗硼剂多以碳化硼为主,一方,碳化硼在渗硼过程中硼势较高,渗硼层脆性大,一方面,碳化硼价格昂贵,因此,开发廉价且能满性能要求的渗硼剂具有一定的现实意义"本文介一种新研制的廉价固体渗硼剂,并对其渗硼工艺行了试验"工件经渗硼后可以显著提高其表面的耐磨性、抗蚀性、红硬性、抗高温氧化性。因而成功地应用在各种工模具及零件上，并取得了明显的经济效果,但常规渗硼需要在较高的温度下进行数小时，能耗很大。章桥新等人2实验证实在微波场下能实现快速渗硼，能有克服常规渗硼的不足，提高效率，