课程作业课程名称科学计算选讲姓名陈学敏学号1012208003学院材.docVIP

课程作业 课程名称：科学计算选讲 姓名：陈学敏 学号：1012208003 学院：材料学院 专业：材料学 有限元法求解偏微分方程在微波渗硼扩散 1. 引 言 模具及零件经渗硼后可显著提高其表面的耐磨性、抗蚀性、红硬性、抗高温氧化性，已广泛应用在各种工模具及零件的生产实践中。常规渗硼温度高、时间长，不符合高效节能原则，近来采用的微波技术可实现快速渗硼，为传统的渗硼技术开辟了一条新途径。MATLAB中的偏微分方程PDE工具箱是用有限元法求解偏微分方程得到数值近似解，可求解线性的椭圆型、抛物线型、双曲线型偏微分方程及本征型方程和简单的非线性偏微分方程。 为研究微波渗硼的过程和机理，本文采用MATLAB中的PDE工具箱，对微波渗硼非稳态扩散问题进行了研究，采用数值计算方法求解了非稳态扩散微分方程，并且计算出了渗硼层的厚度、硼的浓度分布，实现了用MATLAB对微波渗硼过程的模拟仿真。 2.1 微波渗硼扩散方程 研究表明，虽然微波改变了扩散阿累乌斯公式中的扩散常数，但微波场渗硼的三维非稳态扩散微分方程的一般形式仍可以写为: (1) 其中，D为扩散系数，t为时间。结合定解条件(初始、边界条件)，即可实现对微分方程(1)的求解。 .2 MATLAB解法 对大多数的非稳态扩散问题求解，一般很难得到解析解，只能用数值方法求得数值解。数值解方法又分有限元法、有限差分法等，其中有限元法是利用部分插值把区域连续求解的微分方程离散成求解线性代数方程组。用MATLAB的PDE工具箱求解即是采用有限元法求解。用MATLAB的PDE工具箱进行有限元计算前需进行预处理，例如对所求解模型的几何形状进行离散化，即用比较简单的形状和形体来逼近和代替实际的形状和形体，这样可以把比较复杂的曲线和曲面问题转化为相对简单的直线或平面问题。求解非稳态微波扩散问题时，我们共设计了以下3个步骤。 (1) 建立描述微波渗硼的物理模型 根据Fe-B相图，在1000 ℃以下渗硼时，随着硼浓度的升高形成铁硼化合物，当含硼量为8.84 %左右时，形成稳定的中间化合物Fe2B；当含硼量达16.23 %左右时，形成含硼量更高的稳定化合物FeB。 (2) 边界条件和初始条件确定 MATLAB指定了如下3种边界条件:①Dirichlet条件，hu=r；②Neumann条件，n·(c u)+ qu=g；③混合边界条件，即Dirichlet条件和Neumann条件的组合。式中n为垂直于边界的单位矢量，h、r、q、g为常量或与u有关的变量。微波渗硼材料碳含量为0.8 %，渗硼温度为850 ℃。扩散系数D=1.10×10-13 m2/s，时间变化为0～9000 s。简化后的边界条件如所示。 (3) 网格划分及求解偏微分方程创建网格，将网格个数设为个左右，输入时间参数，C(x,y,0)=0，求解偏微分方程。 850 ℃×9000s微波渗硼后硼浓度的三维分布 2.3 数据分析处理 2.3.1 微波渗硼浓度分布 微分方程求解后只能显示出浓度的分布，为进一步求出渗层厚度，还需对计算的结果进行处理。点击Mesh-Export Mesh输出网格数据，点击Solve-Export Solution输出按网格序号排列的数值解u。通过以下语句可将u整理成按距表面距离(即x)排列。 u1=[1e6*p]; 将网格坐标的单位由m化为μm u1=u1(:,1); 只取第一列(即x) u1=[u1,100*u]; 将浓度化为百分数并与x值组合成新矩阵。 u1=sortrows(u1,1); 按x升序排列 x=u1(:,1);取u1的第一列，即X轴 y=u1(:,3602); 取u1的第3602列，即Y轴 图3：850 ℃×1 h微波渗硼后硼浓度的分布曲线 图是根据处理后数据绘出的经850 ×1 h微波渗硼后硼浓度的分布曲线。可看出，这与反应扩散的特点非常符合 2.3.2 渗硼层厚度 为了反映渗层随时间的变化和进行验证，用计算出的各个时刻的渗层厚度和实验测定的渗层厚度进行对比。计算各个时刻的渗层厚度MATLAB程序段如下: [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; end c=c; c=[(0:60:9000),c]; 将渗层厚度加上时间标志 s=c(:,1); x t=c(:,2); y轴渗层厚度 plot（s, t） 图4