计算机在材料学中综合作业

计算机在材料学中的应用综合作业

目录一、传热学模型建立 -2-（一）【模型一】平板焊接过程的二维温度场的数值模拟 -2-（二）【模型二】立方体淬火过程三维非稳态温度场的有限差分计算 -5-（三）【模型三】平板激光淬火过程的三维温度场数值模拟 -9-（四）【模型四】长方体正火过程三维非稳态温度场的有限差分计算 -14-（五）【模型五】端淬实验的数值模拟 -18-三、计算机模拟计算·结果处理 -26-四、团队工作情况 -29-五、成员投入程度评定排序 -30-六、总结与体会 -31-一、传热学模型建立（一）【模型一】平板焊接过程的二维温度场的数值模拟1.模型的建立二维模型的假设：1.材料表面热物性参数不随温度变化。2.不考虑相变潜热。3.考虑工件的辐射与空气的对流换热。4.焊接热源能量分布为高斯分布。5.材料各向同性。6.工件为二维有限大物体。（Z方向厚度不计）2.热传导方程采用直角坐标系，则材料内部热传导方程为：式中，QUOTE，为材料热扩散率，ρ为密度，QUOTE为比热，λ为热传导率，T为温度，t为时间。材料的热物性参数，比热QUOTE，热传导系数λ和密度ρ均不随温度变化。3.边界条件上表面：QUOTE其它表面：QUOTE式中，T是工件表面的温度，Ta是环境温度。n是其它表面的外法线方向。h是材料表面总的换热系数，包括空气对流和热辐射换热。h=QUOTE+QUOTE，QUOTE为对流换热系数，QUOTE为辐射换热系数。q(x,y,t)为焊接热源能量分布函数P为热源有效功率，QUOTE为热源半径，v为热源运动速度。4.有限差分方程由于对称性，取工件整体区域的一半求解。X方向节点数为N1，Y方向节点数为N2。X方向步长为△x，Y方向步长为△y，时间步长为△t。热源沿X方向运动。内部节点差分方程：令QUOTE，QUOTE，QUOTE，则内部节点差分方程变为：四条线边界节点的差分方程：四个角边界节点的差分方程：（二）【模型二】立方体淬火过程三维非稳态温度场的有限差分计算1.模型的建立钢的淬火是将钢加热到临界温度Ac3（亚共析钢）或Ac1（过共析钢）以上温度，保温一段时间，使之全部或部分奥氏体化，然后以大于临界冷却速度的冷速快冷到Ms以下（或Ms附近等温）进行马氏体（或贝氏体）转变的热处理工艺，，以大幅提高钢的刚性、硬度、耐磨性、疲劳强度以及韧性等，从而满足各种机械零件和工具的不同使用要求。也可以通过淬火满足某些特种钢材的铁磁性、耐蚀性等特殊的物理、化学性能。淬火过程是一个涉及到相变，热传导，热对流，热辐射的三维非稳态传热问题。为计算淬火过程的温度场，首先应建立淬火过程的传热学模型。本模型中采用以下几点假设。材料的热物性参数不随温度变化不考虑相变潜热材料为各向同性考虑工件辐射与空气对流换热。工件为立方体2.热传导方程根据工件形状，采用三维直角坐标系，在工件内部的热传导方程为，式中，，为材料的热扩散率，λ为热传导率，ρ为材料密度，为比热，T为温度，t为时间。λ,ρ,均不随温度变化。3.边界条件式中，T为工件表面温度，为淬火介质温度，n为其他表面的的外法线方向，h为淬火介质的对流换热系数。4.初始条件初始时刻工件整体温度分布均匀，即，为一常数，取淬火加热温度，对35CrMo钢，取860℃。5.有限差分方程首先对求解区域离散化，求解区域为立方体，在X,Y,Z方向用网格来离散求解区域，X方向节点数为N1，Y方向节点数为N2，Z方向节点数为N3。X方向步长为△X，Y方向步长为△Y，Z方向步长为△Z。时间步长为△t。内部某节点编号为（i,j,k).左右相邻两节点为(i,j-1,k),(i,j+1,k),前后相邻两节点为(i.j,k-1),(i,j,k+1).上下相邻两节点为（i,j-1,k),(i,j+1,k)。内部节点差分方程为令，，，。则内部节点差分方程变成了。同理，六个面边界节点的差分方程分别为：十二条棱上边界节点的差分方程：八个角边界节点的差分方程为：（三）【模型三】平板激光淬火过程的三维温度场数值模拟1.建立传热学模型图1是平板类工件激光热处理示意图。在平板类工件激光热处理过程的三维传热学模型中采用以下几点假设：（1）材料表面对激光的吸收系数不随温度变化；（2）材料的热物性参数随温度变化；（3）考虑相变潜热；（4）考虑工件的辐射与空气对流换热；（5）入射激光束能量分布为高斯分布（TEM00）；（6）激光功率恒定；（7）激光束扫描速度不变；（8）材料各向同性；（9）工件为三维有限大物体；（10）工件初始温度恒定。图1平板类工件激光热处理示2.热传导方程P.S.如果激光束扫描的速度很快，则沿激光束扫描方向（Y方向）的传热可以忽略，这样三维传热问题就可简化为X方向和Z方向传热的二维传热问题，这就是平板类工件激光热处理过程的二维传热模型。在直角坐标系中，工件内部的热传导方程为：式中，ρ为材料密度，Cp为材料比热，λ为热传导系数，T为温度，t为时间。3.边界条件边界条件为：上表面：其他表面：Q（x,y,t）是激光光斑能量分布函数。P为激光功率，A为吸收系数，R为激光光斑直径，V为激光光斑运动速度。Ta是环境温度，h是材料表面总的换热系数，包括空气对流和热辐射换热。n是其他表面的外法线方向。4.初始条件初始条件为：初始时刻工件整体温度分布均匀，即，Ta为一常数。5.相变潜热的处理激光热处理加热时发生奥氏体相变，冷却时发生马氏体相变，其他相变不可能发生。目前对于在激光相变硬化等超快速加热条件下的奥氏体转变动力学尚不明白，因此，采用把奥氏体化的相变潜热折算成温度降低值的办法，把温度降低值等分，分步加到计算得到的温度场中。HA为奥氏体的相变潜热，n为一整数，一般取为10，或者其他合适的整数。可以用Koistinen-Marburger公式来计算钢铁材料中马氏体相变的转变量Fm为马氏体转变的体积百分数，Ms为马氏体开始转变温度，T为节点温度，k和n为两个常数，与具体的钢种有关。每一步马氏体相变所放出的相变潜热可折算成温升，然后利用ΔTmn对温度场进行修正：Hm为马氏体转变的相变潜热。ΔTmn是每一步马氏体转变放出的相变潜热所引起的温升。有限差分方程（激光热处理的三维差分方程）由偏微分方程和边界条件，初始条件构成传热的定解问题。用能量平衡法推导了这一问题的非均匀空间网络的三维显式差分方程。直角坐标系下差分方程：X方向节点数为101，y方向节点数为101，z方向节点数为51。X方向步长为△X，y方向步长为△Y，z方向步长为△Z。内部节点差分方程为令，，，。则内部节点差分方程变成了。六个面上节点的差分方程分别为：上表面：下表面：左表面：右表面：前表面：后表面：棱上的节点的差分方程：八个角上节点的差分方程为：（四）【模型四】长方体正火过程三维非稳态温度场的有限差分计算1、模型的建立钢铁材料在炉中加热到超过奥氏体转变点后，发生奥氏体转变，然后放到空气中冷却，就可得到含有珠光体的均匀组织，从而实现钢铁材料的正火处理。长方体是一种在工业界中常用的工件，他的正火过程是一个涉及到相变、热传导、热对流的三维非稳态传热过程，可用三维模型近似的描述长方体正火冷却的传热过程，在这个三维模型中采用以下几个假设：工件为三维长方体材料的热物性参数不随温度变化不考虑相变潜热考虑工件与空气的热对流换热材料各向同性2、直角坐标系的建立及网格的划分根据工件的形状，建立如图1所示的坐标系，将OX，OY分成100份，OZ分成50份，即节点为从（i,j,k)=(1,1,1)到（i,j,k)=(101,101,51)。此模型为长×宽×高等于100×100×50（cm）的钢铁材料，材料的密度为7.8g/,比热为0.5J/(g.k),导热系数为0.3。工件冷却过程为从边界向中心逐渐冷却。3、热传导方程在工件内部的热传导方程为，式中，，为材料的热扩散率，λ为热传导率，ρ为材料密度，为比热，T为温度，t为时间。λ,ρ,均不随温度变化。初始条件初始时刻工件整体温度分布均匀，即=，设为900°C。设环境温度=20°C。5、边界条件外表面：式中，T为工件表面温度，为环境温度，n为外表面的外法线方向，h为材料表面的换热系数。h=+,为对流换热系数，为辐射换热系数。6、有限差分方程X方向节点数为101，y方向节点数为101，z方向节点数为51。X方向步长为△X，y方向步长为△Y，z方向步长为△Z。内部节点差分方程为令，，，。则内部节点差分方程变成了六个面上节点的差分方程分别为：上表面：下表面：左表面：右表面：前表面：后表面：十二条棱上的节点的差分方程：八个顶角上节点的差分方程为：（五）【模型五】端淬实验的数值模拟1、模型的建立端淬实验是将标准尺寸的端淬试样(Φ25mm×100mm)奥氏体化后，在专用设备上对其一端面喷水冷却，后沿轴线方向的测出硬度-距水冷端距离的关系曲线的试验方法，是测定钢的淬透性的方法之一。由于模型轴向对称，计算时可以采用二维模型如图。此次模拟实验的模型为标准尺寸（Φ25mm×100mm）的金属试样。先将式样加热到860℃并保温，使式样完全奥氏体化，然后将式样取出，放在端淬试验机上，在式样下端喷20℃的水120s。模型中采取如下假设：材料的热物性参数不随温度变化忽略相变潜热对温度场的影响不考虑轴向传热材料为各向同性材料考虑材料与空气的对流换热，采取第三类边界条件网格划分取式样延一个直径的剖面的右半部分建立如上右图的直角坐标系，将or方向分成125份，oz方向分成1000份。[即从（i，k)=（1,1）到（i，k）=（N1，N2）=（126,1001）]。3.热传导方程材料内部的热传导方程是，式中为密度，为比热，为传导热，T为温度，t为时间。4.边界条件下边：；左边：0；其他边：式中，T是工件表面温度，Ta是环境温度（水温与环境温度相图）。n是其他表面的法线方向。是式样与冷却水的对流换热系数，是式样与空气的换热系数。材料的热物性参数，比热，热传导系数和密度均不随温度的变化而变化。初始条件端淬开始时，式样温度为860℃且温度均匀，完全奥氏体化。冷却水温度与环境温度始终为20℃。有限差分方程r方向的节点数为31，Z方向的节点数为201，r方向的步长为，Z方向的步长为。内部节点的差分方程为令，，则内部节点的差分方程变为令则：四个边上的节点的差分方程：上边：下边：左边：右边：四个角上节点的差分方程：右上角：右下角：左上角：左下角：二、计算机编程模拟·程序设计程序如下：/****有限长圆柱体端淬****//*****初始温度860℃******//****采用显式差分方程****/#include<stdio.h>#include<math.h>#definePI3.14#defineK0.3//导热系数#defineCP0.5//比热#defineP7.8//密度#defineH12.0//与淬火介质的换热系数#defineH20.0003//与空气的换热系数#defineN151#defineN2401doubleS,DS=0.001,S1;//时间与时间步长intT0=860,Ta=20;//初始温度doubleT1[N1][N2],T2[N1][N2];doubleFR,Fz,F1,F2;//傅里叶数doubleDR=0.025,DZ=0.025;//空间步长voidf_number()//傅里叶数的计算{ FR=K*DS/(P*CP*DR*DR); Fz=K*DS/(P*CP*DZ*DZ); F1=8*N1-12; F1=F1/(4*N1-5); F2=8*(N1-1); F2=F2/(4*N1-5); printf("F1=%f\tF2=%f\n",F1,F2); printf("FR=%f\tFz=%f\n",FR,Fz);}voidinitial_temperature()//初始温度{ inti,j; for(i=0;i<N1;i++) for(j=0;j<N2;j++) { T1[i][j]=T0; }}voidheat_conduction()//热传导方程{ inti,j; /*内部节点*/ for(i=1;i<(N1-1);i++) for(j=1;j<(N2-1);j++) { T2[i][j]=T1[i][j]+FR*(1.0-1.0/(2.0*i))*(T1[i-1][j]-T1[i][j])+FR*(1.0+1.0/(2.0*i))*(T1[i+1][j]-T1[i][j])+Fz*(T1[i][j-1]-2*T1[i][j]+T1[i][j+1]); } /*边上的节点*/ /****上边****/ for(i=1;i<(N1-1);i++) { T2[i][N2-1]=T1[i][N2-1]+FR*(1.0-1.0/(2.0*i))*(T1[i-1][N2-1]-T1[i][N2-1])+FR*(1.0+1.0/(2.0*i))*(T1[i+1][N2-1]-T1[i][N2-1])+2*Fz*(T1[i][N2-2]-T1[i][N2-1])+2*Fz*DZ*H2/K*(Ta-T1[i][N2-1]); } /****下边****/ for(i=1;i<(N1-1);i++) { T2[i][0]=T1[i][0]+FR*(1.0-1.0/(2.0*i))*(T1[i-1][0]-T1[i][0])+FR*(1.0+1.0/(2.0*i))*(T1[i+1][0]-T1[i][0])+2*Fz*(T1[i][1]-T1[i][0])+2*Fz*DZ*H1/K*(Ta-T1[i][0]); } /****左边****/ for(j=1;j<(N2-1);j++) { T2[0][j]=T1[0][j]+Fz*(T1[0][j-1]-2*T1[0][j]+T1[0][j+1])+4*FR*(T1[1][j]-T1[0][j]); } /****右边****/ for(j=1;j<(N2-1);j++) { T2[N1-1][j]=T1[N1-1][j]+Fz*(T1[N1-1][j-1]-2*T1[N1-1][j]+T1[N1-1][j+1])+FR*F1*(T1[N1-2][j]-T1[N1-1][j])+FR*F2*DR*H2/K*(Ta-T1[N1-1][j]); } /*角上的节点*/ /***右上角***/ T2[N1-1][N2-1]=T1[N1-1][N2-1]+FR*F1*(T1[N1-2][N2-1]-T1[N1-1][N2-1])+FR*F2*DR*H2/K*(Ta-T1[N1-1][N2-1])+2*Fz*(T1[N1-1][N2-2]-T1[N1-1][N2-1])+2*Fz*DZ*H2/K*(Ta-T1[N1-1][N2-1]); /***右下角***/ T2[N1-1][0]=T1[N1-1][0]+FR*F1*(T1[N1-2][0]-T1[N1-1][0])+FR*F2*DR*H2/K*(Ta-T1[N1-1][0])+2*Fz*(T1[N1-1][1]-T1[N1-1][0])+2*Fz*DZ*H1/K*(Ta-T1[N1-1][0]); /***左上角***/ T2[0][N2-1]=T1[0][N2-1]+2*Fz*(T1[0][N2-2]-T1[0][N2-1])+2*Fz*DZ*H2/K*(Ta-T1[0][N2-1])+4*FR*(T1[1][N2-1]-T1[0][N2-1]); /***左下角***/ T2[0][0]=T1[0][0]+2*Fz*(T1[0][1]-T1[0][0])+2*Fz*DZ*H1/K*(Ta-T1[0][0])+4*FR*(T1[1][0]-T1[0][0]);}voidoutput_1(){FILE*fp1; FILE*fp2; FILE*fp3; FILE*fp4; FILE*fp5;fp1=fopen("1.dat","a"); fp2=fopen("2.dat","a");fp3=fopen("3.dat","a");fp4=fopen("4.dat","a");fp5=fopen("5.dat","a"); fprintf(fp1,"%f\n",T1[0][0]); fprintf(fp2,"%f\n",T1[0][20]); fprintf(fp3,"%f\n",T1[0][40]); fprintf(fp4,"%f\n",T1[0][60]); fprintf(fp5,"%f\n",T1[0][80]); fclose(fp1); fclose(fp2); fclose(fp3); fclose(fp4); fclose(fp5);}intoutput_2(){inti,j;FILE*fp6;fp6=fopen("6.dat","a");for(j=N2-1;j==0||j>0;j=j-2){for(i=0;i<N1;i=i+2) {fprintf(fp6,"%f\t",T1[i][j]);}fprintf(fp6,"\n");}fclose(fp6);}intmain(){ inti,j,m=0,p;printf("请输入淬火时间:");scanf("%lf",&S1); printf("\n");f_number(); printf("请判断傅里叶数，输入任意数字加回车继续\n"); scanf("%d",&p); printf("结果写入文件中，请稍后\n");initial_temperature(); for(S=0;S<S1;S+=DS){heat_conduction(); if(m%1000==0)output_1(); m=m+1; if(m%10000==0) output_2();for(i=0;i<N1;i++)for(j=0;j<N2;j++){T1[i][j]=T2[i][j]; } } return0; }三、计算机模拟计算·结果处理不同时刻端淬试样的温度场情况50s 100s150s200s250s300s温度标尺从以上图中可以看出，水冷端开始时温度降低很快，降低到100℃后开始减慢；试样上端开始时温度变化很慢。整个的端淬过程中，水冷端温度迅速下降，低温区逐渐向上移动。端淬过程的时间——温度图及分析材料的CCT曲线由cct曲线可知，材料的鼻尖温度大约为550℃，时间为1s。对比上面的端淬过程时间温度曲线可知，距水冷端2mm处1s时温度大约为550℃由此可知，材料可淬成马氏体的深度大约在2