【2019年整理】计算材料学论文

“L”形AZ91铸造镁合金凝固过程数值模拟李宾112161366摘要：以“L”形镁合金铸件为例，对其凝固过程进行了合理的假设和简化，建立了凝固过程的二维非稳态温度场的计算数学模型；并运用有限差分方法对模型进行离散，得到显示差分格式。据此，利用Matlab编制了计算机程序，获得了铸件在凝固过程中的温度分布规律，数值模拟能够反映铸件冷却过程温度场的动态变化，表明凝固潜热的影响，也为提高铸件质量，确定凝固时间和优化工艺参数提供了参考。关键词：AZ91镁合金；铸件；潜热处理；凝固过程；温度场；有限差分；数值模拟铸造是获得机械产品毛坯的主要方法之一，对国民经济的发展起着十分重要的作用。据有关资料表明，铸件在机床、内燃机、重型机器中占70%-90%，农业机械中占40%-70%。此外，铸造在冶金工业、航空、航天、船舶等重要产业部门都发挥着重要作用[1]。然而，传统的铸件生产不仅周期较长、浪费大，而且产品质量也难以保证。而在实际生产之前对铸件在浇注、凝固过程中可能产生的缺陷采用计算机数值模拟的方法进行有效的预测，从而在浇注前就采取有效对策，以确保铸件的质量。由于镁合金具有较高的比强度、良好的减震性和切削加工性及尺寸的稳定性等,这些优良的特性使它成为非常重要的现代工业材料。目前镁合金压铸件被广泛应用于汽车、航空航天和计算机制造业等各个领域[2]。而镁合金应用的关键的一个环节是在材料制备过程中的金属凝固过程。因此可以通过数值模拟技术，凝固理论直接与镁合金铸造工艺联系起来，推动铸造工艺技术的发展，具有着重要的意义。本文以“L”形镁合金铸件为例，进行温度场的数值模拟，模拟结果反映了铸件冷却过程温度场的动态变化。1模型的建立1.1研究对象研究对象为L形镁合金铸件及其砂型，如图1（a）。砂型壁厚为120mm铸件厚度为50mm，横截面如图1（b）所示。本文研究边界温度为定值的情况，镁合金初始温度为其液相线温度596℃，砂型恒为25℃。（a）镁合金铸件及其砂型（b）铸件横截面图图1镁合金铸件及其砂型示意图1.2铸件温度场的控制方程人们根据计算传热学、金属凝固理论以及能量守恒原理等理论,建立了铸件凝固过程三维温度场数值模拟的数学模型。这个模型的核心内容可以简单归纳为以下两个方面[4]:(1)三维不稳定导热偏微分方程的建立和求解;(2)单值条件的确定。铸件单元三维不稳定导热偏微分方程的直角坐标描述为:tfLzTyTxTtTCsp222222(1)式中x,y,z——坐标轴方向的长度,mT——温度,Kt——时间,s——密度,kg/m3——热导率,W/(mK)pC——比热,kJ/(kgK)sf——固相率,％L——合金的结晶潜热，KJ/Kg考虑到AZ91镁合金的结晶温度范围较宽,本文采用等价比热容法对潜热作合理的处理[3]。其推导如下。固相率与温度密切相关，则tfs可表示为tTtfs，式（1）可表示成：222222zTyTxTtTCp（2）其中：SpSLspLppTTCTTTtfLCTTCC。式中：LT为液相线温度；ST为固相线温度；等效比热法是将潜热项折合成比热容，以等效比热容pC代替pC。在实际计算中假设潜热均匀释放,则)(0SLTTLL,)(0SLpppTTLCLCC1.3基本假设和处理影响铸件凝固过程的因素非常多,在求解中若要把所有复杂因素都考虑进去是不现实的。因此在铸件凝固过程复杂的实际条件下,需做一些合理的简化和假设:(1)认为液体金属在瞬时充满铸型后开始凝固。(2)忽略钢水在凝固过程中的相对流动,即只考虑凝固过程中钢液的对流传热。(3)由于铸模在钢液凝固过程中的温度变化不是很大,因而将它的热物性参数看成是不随温度变化的常数。(4)不考虑合金的过冷,假定凝固是从给出的液相线温度开始,固相线温度结束。1.4热物性参数铸造模拟中涉及到的最重要的参数就是热物性参数，因此准确的热物性参数是获得准确模拟结果的必要条件。在通常情况下，铸件与铸型的密度、比热、导热系数、换热系数等热物性值不是一恒定的值，而随温度不同有所变化，但数值变化不是很大，在模拟计算时均将系统的热物性值作为常数处理。因为这样不仅降低对计算机内存的要求，大大提高模拟运算速度，而且对铸件凝固推移进程、最后凝固区域位置等相对值几乎没有影响。(1)铸件AZ91的热物性参数:表1AZ91镁合金热物性值[5,6]液相线温度596℃固相线温度468℃液态合金密度1650Kg/m3固态合金密度1750Kg/m3液态导热系数80W/(Km)固态导热系数105W/(Km)液态比热1.35KJ/(KgK)固态比热1.20KJ/(KgK)凝固潜热370KJ/Kg(2)呋喃树脂砂型的热物性参数:[5]呋喃树脂砂的导热系数为1.01W/(mk)，密度为1550kg/m3，比热为1.08kJ/(kgK)。本文中为方便计算，假设砂型保持恒温。1.5初始条件及边界条件的确立[7]1）初始条件初始条件就是在过程开始的时刻,即t=0时整个温度场内部的温度分布。它可以是均匀的，此时有到T|t=0=T0，也可以是不均匀的，各点的温度值已知或者遵从某一函数关系，即T|t=0=T0(x，y，z).本文中假设模具型腔始终保持初始温度为室温25℃,钢水充满型腔后其初始温度为568℃。2）边界温度热传导问题的边界条件常以三种形式给出，分别称之为第一类边界条件、第二类边界条件和第三类边界条件.第一类边界条件(又称Dirichlet边界条件)，也称已知温度的边界.第二类边界条件(又称为Neumann边界条件)，也称已知热流密度的边界.第三类边界条件(又称为Robin边界条件)，也称为对流换热边界.1.6数学模型的有限差分方法文中所建立二维热传导方程初边值问题为（3），文中采用显式差分方法对方程进行离散。若令,405.020,805.040yLNxLNyxxL和yL分别为x和y方向的尺寸。边界上的点属于点集B：,41)，12,11=j81,，11,12,=i(81=i1=i41=j1=j铸件上的点集I：)40,,13,12;10,,3,2()10，,3,2；80，,3,2(jiji再令22)(,)(yCtfxCtfpypx，则显式差分格式为：IjiTBjiTjiTffTTfTTfTjinjinjiyxnjinjiynjinjixnji,℃596,℃2540,,3,2;80,,3,22210,,,1,1,,1,11,（3）2模拟结果及分析：模拟开始时，如图2所示，铸件的温度基本保持在浇注时的温度596℃，铸件主要部位几乎没有发生温差，这是一个极端条件。但铸件周边变化较明显，由于边界温度的作用，铸件周边出现温度梯度，但等温线较密集，温度梯度比较大，与铸件主体形成强烈对比。随着凝固的进行，由图3、图4和图5可知，铸件主体已经出现了温差，L形较短的一边度梯度较大，这是短边较长边受到比较大的冷却效应的缘故；由铸件中心到边缘，温度逐渐降低，这也与实际相符；L形拐角的地方温度梯度较小，从图中可以看出，拐角的地方凸出一个峰，表明此处温度较其他区域温度高，应是最后达到固相线温度的，此处可预测铸件最可能会出现缺陷，从图4的等温线也可以看出以上规律。随着时间的推移，潜热开始释放，图4的等温线比图3的稀疏，铸件温度又开始变得均匀一些；由图中等温线可以看到低温铸件向铸件中心推移。图6是同一点在不同浇注温度下的冷却曲线，此点温度一直在下降，开始时，周围温差相差不大，温度下降速度较平缓，然后温度下降速度增大，由于开始时，液相较多，下降较快，但随着凝固的进行，固相越来越多，由于凝固潜热的作用，温度下降速度变得较以前时刻小；四条曲线的趋势是相同的，但浇注温度越高，需要更长的时间才能冷却到同一温度，这与现实结果相同,596℃浇注温度下需要大约200s点（6,6）冷却到固相线。由图7可以看出，因为固相潜热的释放，使铸件冷却曲线变得缓慢。图2初始t=0.5s时温度场图3t=60s时温度场图4t=150s时温度场图5刚低于固相线时温度场图6不同浇注温度时点（6,6）处冷却曲线图7浇注温度为596℃时，潜热对冷却过程的影响3结论1）通过对凝固过程相关因素作出合理假设,建立了“L”形AZ91镁合金铸件凝固过程的数学模型。模拟结果与实际过程基本吻合,可以作为实际铸造过程的分析依据。2）利用温度场的模拟结果可以预测可能出现的缺陷及缺陷出现的部位,为优化铸造工艺提供了科学依据。3）由于凝固潜热的释放，会降低合金的冷却速度。4）浇注温度不同，铸件某点温度的下降走势相同，浇注温度越高，需要更长的时间才能达到固相线温度。参考文献[1]李青松.AZ91铸造镁合金凝固过程数值模拟[D].华东交通大学硕士学位论文.2009[2]于彦东,马秋,等.镁合金压铸件凝固过程计算机模拟[J].热加工工艺,2005，No.12[3]周建兴,刘瑞祥,等.凝固过程数值模拟中的潜热处理方法[J].铸造，2001,Vol.50,No.7[4]林首位.铸件凝固过程三维温度场数值模拟研究[D].太原:华北工学院工学硕士学位论文,2001[5]严力.AZ91镁合金调压铸造充型能力研究[D].西北工业大学硕士学位论文，2006[6]赵彦民.定向凝固AZ91镁合金工艺与组织研究[D].太原科技大学硕士学位论文，2010[7]许芝卉.凝固过程的温度场数值模拟[J].数学的实践与认识，2010,Vol.40,No.13附录：源程序代码%程序一：取初始一张、凝固过程二张，结束一张温度场分布图clear,clc,clfL1=40;L2=20;N=80;M=40;%长为40cm,80等分，宽为20cm，40等分。T0=596;%铸件初始温度a=0.01;%a=λ/ρCdt=0.5;%时间步长0.5sdx=L1/N;dy=L2/M;M1=a*dt/(dx^2);M2=a*dt/(dy^2);T=T0*ones(N+1,M+1);T1=T0*ones(N+1,M+1);t=0;l=0;k=0;Tc=[];%采集点（6,6）处温度值fori=1:81forj=1:41ifj==1|j==41|i==1|i==81T(i,j)=25;%边界点温度为25℃elseT(i,j)=596;%液相线T0=596℃endendendfori=11:81forj=11:41T(i,j)=25;endendif(2*M1+2*M2=1)%判断是否满足稳定性条件while(max(max(T))468)%刚低于固相线时模拟结束t=t+dt;k=k+1;fori=2:80forj=2:40T1(i,j)=M1*(T(i-1,j)+T(i+1,j))+M2*(T(i,j-1)+T(i,j+1))+(1-2*M1-2*M2)*T(i,j);endendfori=2:80forj=2:40T(i,j)=T1(i,j);endendift==0.5i=1:81;j=1:41;[x,y]=meshgrid(j,i);figure(1);subplot(1,2,1);mesh(x,y,T(i,j))%画出0.5s时温度场axistight;ylabel('y','FontSize',14);xlabel('x','FontSize',14);zlabel('T/℃','FontSize',14)title('0.5s时温度场模拟图','FontSize',18)subplot(1,2,2);[C,H]=contour(x,y,T(i,j));clabel(C,H);axissquare;xlabel('x','FontSize',14);ylabel('y','FontSize',14);title('0.5s时模拟等温线图','Fon