第23卷第2期2007年10月 山西大同大学学报(自然科学版)Jou rnal of Shanxi Dat ong University(Nat u ral Science) Vol .23.No.2O ct .2007凝固过程数值计算的模拟许芝卉(山西大同大学数学与计算机科学学院,山西大同037009)摘 要:通过对铸件凝固过程中各换热边界条件的研究,建立了凝固过程的二维非稳态温度场计算数学模型;并运用了有限差分方法对模型进行离散,得到大型方程组,并利用超松驰迭代法(即SOR 法)解该方程组,据此,利用Turbo C 编制了计算机程序.上机运行结果表明,可较满意地模拟铸件凝固过程温度场的分布.关键词:铸件 凝固过程 有限差分法 温度场 数值模拟中图分类号:O242.1 文献标识码:A 文章编号:167420874(2007)022******* 随着材料加工科学和技术的发展,凝固技术不断进步,并成为一门科学.凝固科学与凝固技术相互促进,并进一步促进先进材料加工科学和技术的发展,在国民经济和国防建设中发挥着重要作用.凝固科学与技术体系的发展是建立在现代科学的基础上,不断地以数学、物理、化学及工程科学的新成就充实自己,同时又不断从冶金、晶体生长、材料科学、空间科学、化工、电子、信息、计算科学等领域汲取营养,迄今已初步构筑成一个凝固科学与材料凝固加工技术的应用与研究体系,其应用目标是以控制材料组织结构为核心,进而控制构件形状并获得所需要的性能.随着社会需求与科学技术进步的牵引,特别是新材料与制备加工技术的需要,推动凝固科技向更深、更高、更精细和开发新的、综合性更强、超常规的方向发展;并直接推动新材料的研究开发.凝固模拟是采用计算机模拟铸件的凝固过程,并以直观可视的形式把铸件的凝固过程和缺陷位置等形象地显示出来,它是虚拟制造技术在铸造行业的应用和体现.计算机模拟铸件的凝固过程不仅为铸造工艺设计提供了一种科学依据,而且为彻底改变铸造生产劳动密集型传统模式,提高铸造技术水平,进一步发展铸造工艺计算机设计(CA D )、计算机辅助工程(CAE )、计算机辅助制造(CA M )技术奠定了基础.计算机数值模拟包括几何模型的建立,网格点划分,求解条件(初始条件和边界条件)的确定,数值计算,计算结果的处理及图形显示.其所用的数值模拟的基本方法主要是有限差分法,有限元法和边界元法.目前铸造中应用的较多的领域是:(1)凝固过程数值模拟,主要进行铸造过程的传热分析.包括数值计算方法的选择,潜热处理、缩孔缩松预测判别,铸件、铸型界面传热问题处理.(2)流动场数值模拟,涉及动量、能量与质量传递,其难度较大.(3)铸造应力模拟,此项研究开展较晚,主要进行弹塑性状态应力分析.(4)组织模拟,目前尚处起步阶段.分宏观、中观和微观模拟.能计算形核数,分析初晶类型,枝晶生长速度,模拟组织转变,预测机械性能.1 数学模型1.1 凝固过程的导热方程为揭示温度随空间与时间的变化规律,须求助于傅立叶导热微分方程.三维非稳态导热偏微分方程如下:收稿日期5作者简介许芝卉(66),女,山西应县人,硕士,副教授,研究方向计算数学:2007-04-2:19-:.ρc ρ9T 9t =λ(92T 9x 2+92T 9y 2+92T 9z2)+q (1.1)q =ρL 9f s9t(1.2)式中ρ为材料密度(kg/m 3),c ρ为材料的定压比热(J /kg k ),λ为材料的导热系数(W /m K ),T 为温度(K),t 为时间,q 为单位体积的金属在单位时间内释放的潜热,f s 为固相率,L 为凝固潜热(J /kg ).(1.1)和(1.2)联立整理得:ρ(c ρ-L9f s 9T )9T9t =λ(92T 9x 2+92T 9y 2+92T 9z2)(1.3)令c ′=c ρ-L 9f s9T则有:ρc ′9T 9t =λ(92T 9x 2+92T 9y 2+92T 9z2)(1.4)方程(1.4)的形式更适合于利用差分方法求解.剩下的问题就是看如何确定c ′了,图1为合金状态图的一角.图1 具有凝固温度范围的合金状态图由于k 0=c s c L,m =-T L -T c L -c 0,f s =c L -c 0c L -c s上面3式中k 0为分配平衡系数,c s 为固相中溶质浓度,m 为液相线斜率,c L —液相中溶质浓度,T L —液相线温度,T s —固相线温度.由上面3个方程联立可得:f =T -T L()(+T T L )(5)9f 9T =()(+T T L )因此有c ′=c ρ-m c 0L(1-k 0)(c 0m +T -T L )(1.6)另一方面由(1.1)可变成另一种形式:9T 9t =a (92T 9x 2+92T 9y 2+92T9z 2)+q ′(1.7)其中a =λ/ρc ρ,a 为导温系数q ′=q /ρc ρ.上述为非稳态导热条件下的傅立叶导热微分方程.当为稳态导热条件时,有9T /9t =0.1.2 初始条件和边界条件1)初始条件初始条件是指求解问题的初始温度场,也就是在零时刻温度场的分布.它可以是均匀的,此时有T |t =0=T 0,也可以是不均匀的,各点的温度值已知或者遵从某一函数关系,即T |t=0=T 0(x,y,z ).2)边界条件热传导问题的边界条件常以三种形式给出,分别称之为第一类边界条件、第二类边界条件和第三类边界条件.第一类边界条件(又称D irichle t 边界条件),也称已知温度的边界.第二类边界条件(又称为Neum ann 边界条件),也称已知热流密度的边界.第三类边界条件(又称为Robin 边界条件),也称为对流换热边界.2 数学模型的有限差分方法 文中所用几何体网格为图2所示,所建立二维热传导方程初边值问题为(2.1),文中采用交替方向的隐式差分格式中的Douglas -Rachf ord 格式对方程进行离散.Douglas 和R achford (1956)年提出了另一个交替方向隐式差分格式,即D -R 格式,它是第一个能推广到三维情形的保持无条件稳定的交替方向隐式格式.图 几何体网格41 山西大同大学学报(自然科学版) 2007年s 1-k 0c 0m -1.s m c 01-k 0c 0m -22考察如下二维热传导方程初边值问题9T 9t =a (92T 9x 2+92T9y2)+q T (x,y,0)=7000C 0≤x,y ≤1T (0,y,t )=T (1,y,t )=T (x,0,t )=T (x,1,t )=T 00C, 0≤x,y ≤1,t >0(2.1)逼近(2.1)D -R 格式为:T 3j,k -T nj,kτ=a (δ2x T 3j ,kh 2x+δ2y T nj ,kh 2y )+qT n+1j,k-T 3j ,kτ=a (δ2y T n+1j ,k h 2y -δ2y T nj,kh 2y)上式可写为(I -r x a δ2x )T 3j ,k =(I +r y a δ2y )T nj,k +q τ(I -r y a δ2y )T n +1j ,k=T 3j ,k -r y a δ2y T nj ,k(2.2)由上式第二个方程中解出T 3j ,k ,并代入上式中的第一个方程,得:(I -r x a δ2x )(I -r y a δ2y )T n +1j,k=(I +r x r y a 2δ2x δ2y )Tn j ,k +qτ利用初始条件及边界条件可得如下方程组ATn +1=B T n+b n(2.3)其中A =A 1A 2A 2A 1A 2ωωω A 2A 1A 2 A 2A 1B =B 1B 2B 2B 1B 2ωωωB 2B 1B 2B 2B 1b =(b 1,b 2,…,b M -1,b M )T3 用迭代法解大型稀疏矩阵3.1 基本迭代法的加速(超松驰法,即S OR 法)超松驰迭代法的计算公式为:x k +1=(ω)x (k)+ω(∑j =x(+)j∑j=+x ()j)=x ()+ωa ii(b i -∑i-1j =1aiix(k )j-∑nj=ia(k)ii)i =1,2,…,n (3.1)(3.1)可以写成x (k+1)i=x (k )i+△x i△x i =ωa ii[b i -∑i-1j =1a ii x (k+1)j-∑nj=ia ii x (k )j]k =0,1,2,…;i =1,2,…,n(3.2)(3.1)或(3.2)称为逐次超松驰法,简称S OR 法,ω为松驰因子.S OR 法的收敛速度与ω的取值有关,当ω=1时,它就是G -S 法.因此,可选取ω值使(3.2)的收敛速度较G -S 法快,从而起到了加速作用.为了讨论ω的取值与收敛性的关系,特将(3.1)改写成矩阵形式.由(3.1)可得:a ii x (k +1)i =(1-ω)a ii x (k)i +ω(b i -∑i -1j=1a ii x(k +1)j-∑nj=i+1a ii x (k)j)则上式可写成矩阵的形式:D x (k +1)=(1-ω)D x (k)+ω(b +L x (k+1)+U x (k))整理后得:(D -ωL )x (k+1)=[(1-ω)D +ωU ]x(k)+ωb,x(k +1)=(D -ωL )-1[(1-ω)D +ωU ]x(k)+ω(D -ωL )-1b .令Lω=(D -ωL )-1[(1-ω)D +ωU ]f =ω(D -ωL )-1b则x(k +1)=L ωx(k )+f(3.3)其中L ω为S OR 法的迭代矩阵.显然,(3.3)收敛的充要条件为ρ(L ω)<1,可以证明ρ(L ω)≥|ω-1|,故若(3.3)收敛,则|ω-1|≤ρ(L ω)<1,即0<ω<2是(3.3)收敛的必要条件.下面再给出一个充分条件.定理1 若系数矩阵A 是严格对角占优的或不可约对角占优的,且松驰因子ω∈(0,1),则S OR 收敛.定理2 若系数矩阵A 是实对称的正定矩阵,则当0<ω<2时S OR 收敛.3.2 数值计算与分析采用数值模拟的方法,可以对拟实的铸件和铸型进行模拟计算,方便对于不同热物性影响效果的观察铸件凝固过程中热物性参数主要包括铸件和512007年 许芝卉:凝固过程数值计算的模拟 i1-ia iib i -i -11a i i k 1-n i 1a i i k k i.铸型的比热c ρ,潜热L,导热系数λ和密度ρ,它们皆随温度而变化,这种变化反映到数值模拟方法中,每一时间步长都可能包含很多个热物性值,对模拟精度有影响.本论文中所用几何体网格为图2所示,离散方程和边界条件为(2.1),用超松驰迭代法解的矩阵方程组为(2.3),在程序实现中松驰因子取ω=1.3,初始温度取T =700℃,边界温度取T 0=20℃,导热系数取λ=35.97W /m K,材料密度ρ=7200kg/m 3,比热容取c ρ=630J /kg .k,计算结果由图3和图4表示.图3 节点与温度的关系图 图4 时间与温度的关系参考文献[1]柳百成.铸件充型凝固过程数值模拟国内外研究进展[J ].铸造,1999(8):40245.[2]赵健,张毅.铸件凝固电子计算机数值模拟发展概况[J ].铸造,1985(5):126.[3]康进武.铸钢件凝固过程热应力数值模拟研究[M ].北京:清华大学出版社,1998.[4]徐树方,高立,张平文.数值线性代数[M ].北京:北京大学出版社,2000.[5]余德浩,汤华中.微分方程数值解法[M ].北京:科学出版社,2003.[6]Radovi cZ,Lal ovic M.Nume ri ca l si m ulation of stee l ingot solidifi ca ti on p r ocess[J ].Jou rna l of M a teria ls P rocessing Technology ,2005,160:1562159.S i m u l a t i o n for Num er ica l Eva lua tion of the Solid ifica tion ProcessX U Zhi 2hui(School of M athem atic s and Co mpute r Science,Shanxi Dat ong University,Da t ong Shanxi,037009)Ab stra ct:We hav e constructed a m athe m ati ca l mode l of t w o 2di mensi ona l unsteady te mpe rature f i e ld by ana lyzi ng casting s olidifi 2cati on proce ss fro m exchangi ng v a ri ou s hea t boundary conditi ons .The system of equati ons is obta ined by using F inite diffe rence me t hod t o dis pe rs e the model and by us e of S OR to re s olve this syste m of equati on,writ e the co rresponding co mputati on program by Turbo C .B y analyzing the re s ult,it sh ows that we have mo re approvingl y si m ul a ted the s olidificati on p rocess temperature fie l d .Key wor ds:casting;solidification p r ocess ;finite difference me t hod;temperature fie l d;stre ss field;nu me ri ca l si mulati on61 山西大同大学学报(自然科学版) 2007年。
![一种基于冷态实验模拟凝固过程的装置[发明专利]](https://img.taocdn.com/s1/m/d256b8cc3169a4517723a3f4.png)
![一种金属凝固过程的模拟方法[发明专利]](https://img.taocdn.com/s1/m/25f9ee194693daef5ff73d38.png)
