【2019年整理】計(jì)算材料學(xué)論文

1、“L ”形AZ91鑄造鎂合金凝固過程數(shù)值模擬李賓112161366摘 要：以“L”形鎂合金鑄件為例，對(duì)其凝固過程進(jìn)行了合理的假設(shè)和簡化， 建立了凝固過程的二維非穩(wěn)態(tài)溫度場(chǎng)的計(jì)算數(shù)學(xué)模型；并運(yùn)用有限差分方法對(duì)模型進(jìn)行離散，得到顯示差分格式。據(jù)此，利用Matlab編制了計(jì)算機(jī)程序，獲得了鑄件在凝固過程中的溫度分布規(guī)律，數(shù)值模擬能夠反映鑄件冷卻過程溫度場(chǎng)的 動(dòng)態(tài)變化，表明凝固潛熱的影響，也為提高鑄件質(zhì)量，確定凝固時(shí)間和優(yōu)化工藝 參數(shù)提供了參考。關(guān)鍵詞:AZ91鎂合金；鑄件；潛熱處理；凝固過程；溫度場(chǎng)；有限差分；數(shù)值 模擬鑄造是獲得機(jī)械產(chǎn)品毛坯的主要方法之一， 對(duì)國民經(jīng)濟(jì)的發(fā)展起著十分重要 的作用。據(jù)

2、有關(guān)資料表明，鑄件在機(jī)床、內(nèi)燃機(jī)、重型機(jī)器中占70%-90%，農(nóng)業(yè)機(jī)械中占40%-70%。此外，鑄造在冶金工業(yè)、航空、航天、船舶等重要產(chǎn)業(yè) 部門都發(fā)揮著重要作用。然而，傳統(tǒng)的鑄件生產(chǎn)不僅周期較長、浪費(fèi)大，而且 產(chǎn)品質(zhì)量也難以保證。而在實(shí)際生產(chǎn)之前對(duì)鑄件在澆注、 凝固過程中可能產(chǎn)生的 缺陷采用計(jì)算機(jī)數(shù)值模擬的方法進(jìn)行有效的預(yù)測(cè)，從而在澆注前就采取有效對(duì) 策，以確保鑄件的質(zhì)量。由于鎂合金具有較高的比強(qiáng)度、良好的減震性和切削加工性及尺寸的穩(wěn)定性 等，這些優(yōu)良的特性使它成為非常重要的現(xiàn)代工業(yè)材料。目前鎂合金壓鑄件被廣 泛應(yīng)用于汽車、航空航天和計(jì)算機(jī)制造業(yè)等各個(gè)領(lǐng)域。而鎂合金應(yīng)用的關(guān)鍵的一個(gè)環(huán)節(jié)是在材料

3、制備過程中的金屬凝固過程。 因此 可以通過數(shù)值模擬技術(shù)，凝固理論直接與鎂合金鑄造工藝聯(lián)系起來， 推動(dòng)鑄造工 藝技術(shù)的發(fā)展，具有著重要的意義。本文以“L”形鎂合金鑄件為例，進(jìn)行溫度場(chǎng)的數(shù)值模擬，模擬結(jié)果反映了 鑄件冷卻過程溫度場(chǎng)的動(dòng)態(tài)變化(1)1 模型的建立1.1研究對(duì)象研究對(duì)象為L形鎂合金鑄件及其砂型，如圖 1 (a)。砂型壁厚為120mm鑄 件厚度為50mm,橫截面如圖1 (b)所示。本文研究邊界溫度為定值的情況，鎂(a)鎂合金鑄件及其砂型(b)鑄件橫截面圖圖1鎂合金鑄件及其砂型示意圖1.2鑄件溫度場(chǎng)的控制方程人們根據(jù)計(jì)算傳熱學(xué)、金屬凝固理論以及能量守恒原理等理論，建立了鑄件凝固過程三維溫度

4、場(chǎng)數(shù)值模擬的數(shù)學(xué)模型。這個(gè)模型的核心內(nèi)容可以簡單歸納為 以下兩個(gè)方面：(1) 三維不穩(wěn)定導(dǎo)熱偏微分方程的建立和求解；(2) 單值條件的確定。鑄件單元三維不穩(wěn)定導(dǎo)熱偏微分方程的直角坐標(biāo)描述為式中x,y,z坐標(biāo)軸方向的長度,m式中x,y,z坐標(biāo)軸方向的長度,mfs固相率，汀二 CP :tcX22一一 2cy也丿也丿(2)T溫度,Kt時(shí)間,s?-密度,kg/m-熱導(dǎo)率,W/(m K)Cp 比熱，kJ/(kg K)L 合金的結(jié)晶潛熱，KJ/Kg考慮到AZ91鎂合金的結(jié)晶溫度范圍較寬，本文采用等價(jià)比熱容法對(duì)潛熱作 合理的處理3。其推導(dǎo)如下。丙rf fT固相率與溫度密切相關(guān)，則寸可表示為-，式可表示成:其

5、中:CpT z TLIrfCp=Cp L TL T r TsP】P 戲CpT TS式中：兀為液相線溫度；TS為固相線溫度；等效比熱法是將潛熱項(xiàng)折合成比熱容， 以等效比熱容Cp代替Cp。在實(shí)際計(jì)算中假設(shè)潛熱均勻釋放,則L。二 L(TL-TS),Cp二CpL。二CpL(TL- TS )1.3基本假設(shè)和處理影響鑄件凝固過程的因素非常多，在求解中若要把所有復(fù)雜因素都考慮進(jìn)去是不現(xiàn)實(shí)的。因此在鑄件凝固過程復(fù)雜的實(shí)際條件下，需做一些合理的簡化和假 設(shè)：(1)認(rèn)為液體金屬在瞬時(shí)充滿鑄型后開始凝固 熱。(3) 由于鑄模在鋼液凝固過程中的溫度變化不是很大 ，因而將它的熱物性參數(shù) (2)忽略鋼水在凝固過程中的相對(duì)

6、流動(dòng)，即只考慮凝固過程中鋼液的對(duì)流傳7看成是不隨溫度變化的常數(shù)。(4) 不考慮合金的過冷，假定凝固是從給出的液相線溫度開始 ，固相線溫度結(jié) 束。1.4熱物性參數(shù)鑄造模擬中涉及到的最重要的參數(shù)就是熱物性參數(shù)，因此準(zhǔn)確的熱物性參數(shù)是獲得準(zhǔn)確模擬結(jié)果的必要條件。在通常情況下，鑄件與鑄型的密度、比熱、導(dǎo) 熱系數(shù)、換熱系數(shù)等熱物性值不是一恒定的值，而隨溫度不同有所變化，但數(shù)值變化不是很大，在模擬計(jì)算時(shí)均將系統(tǒng)的熱物性值作為常數(shù)處理。因?yàn)檫@樣不僅 降低對(duì)計(jì)算機(jī)內(nèi)存的要求，大大提高模擬運(yùn)算速度，而且對(duì)鑄件凝固推移進(jìn)程、 最后凝固區(qū)域位置等相對(duì)值幾乎沒有影響。(1) 鑄件AZ91的熱物性參數(shù)：表1 AZ91鎂

7、合金熱物性值5,6液相線溫度596 C固相線溫度468 C液態(tài)合金密度1650Kg/m3固態(tài)合金密度1750Kg/m3液態(tài)導(dǎo)熱系數(shù)80W/(K m)固態(tài)導(dǎo)熱系數(shù)105W/(K m)液態(tài)比熱1.35KJ/(Kg K)固態(tài)比熱1.20KJ/(Kg K)凝固潛熱370KJ/Kg 呋喃樹脂砂型的熱物性參數(shù):5呋喃樹脂砂的導(dǎo)熱系數(shù)為1.01W/(m k)，密度為1550kg/m3，比熱為1.08kJ/(kg K)。本文中為方便計(jì)算，假設(shè)砂型保持恒溫。1.5初始條件及邊界條件的確立1)初始條件j = 1 j = 41 i = 1 i(i 二 2,3,80; j 二 2,3,二 C；C）y）2，則顯式差分格

8、式為:初始條件就是在過程開始的時(shí)刻，即t=0時(shí)整個(gè)溫度場(chǎng)內(nèi)部的溫度分布。它可 以是均勻的，此時(shí)有到 T|t=o=To,也可以是不均勻的，各點(diǎn)的溫度值已知或者遵 從某一函數(shù)關(guān)系，即T|t=o=To（x，y，z）.本文中假設(shè)模具型腔始終保持初始溫度為 室溫25 C ,鋼水充滿型腔后其初始溫度為 568 C。2）邊界溫度熱傳導(dǎo)問題的邊界條件常以三種形式給出，分別稱之為第一類邊界條件、第 二類邊界條件和第三類邊界條件.第一類邊界條件（又稱Dirichlet邊界條件），也稱已知溫度的邊界.第二類邊界條件（又稱為Neumann邊界條件），也稱已知熱流密度的邊界.第三類邊界條件（又稱為Robin邊界條件），

9、也稱為對(duì)流換熱邊界.1.6數(shù)學(xué)模型的有限差分方法文中所建立二維熱傳導(dǎo)方程初邊值問題為（3），文中采用顯式差分方法對(duì)方程進(jìn)行離散。若令 N = J. 7X = 40 0. 5 = 80, N = J 勺=20 0. 5 二 40, Lx 和Ly分別為x和y方向的尺寸。邊界上的點(diǎn)屬于點(diǎn)集 B :=81 (i = 11,12,，81,j = 11,12, ,41)鑄件上的點(diǎn)集I:,，0) U(i = 2,3,，10 j = 12,13,，40)再令 f x2 , 1 y必cCp(Ax)2凸=GTiL +T=,j Wil +T)+(1 2f 也2f 曲 Tnj (=2,3，8Q j = 2,3，,40

10、 )T = 25C(j ,j e B)T0= 596C(J , j | )2模擬結(jié)果及分析：模擬開始時(shí)，如圖2所示，鑄件的溫度基本保持在澆注時(shí)的溫度 596C，鑄 件主要部位幾乎沒有發(fā)生溫差，這是一個(gè)極端條件。但鑄件周邊變化較明顯，由 于邊界溫度的作用，鑄件周邊出現(xiàn)溫度梯度，但等溫線較密集，溫度梯度比較大， 與鑄件主體形成強(qiáng)烈對(duì)比。隨著凝固的進(jìn)行，由圖 3、圖4和圖5可知，鑄件主體已經(jīng)出現(xiàn)了溫差，L形較短的一邊度梯度較大，這是短邊較長邊受到比較大的 冷卻效應(yīng)的緣故；由鑄件中心到邊緣，溫度逐漸降低，這也與實(shí)際相符；L形拐角的地方溫度梯度較小，從圖中可以看出，拐角的地方凸出一個(gè)峰，表明此處溫 度較

11、其他區(qū)域溫度高，應(yīng)是最后達(dá)到固相線溫度的，此處可預(yù)測(cè)鑄件最可能會(huì)出 現(xiàn)缺陷，從圖4的等溫線也可以看出以上規(guī)律。隨著時(shí)間的推移，潛熱開始釋放， 圖4的等溫線比圖3的稀疏，鑄件溫度又開始變得均勻一些；由圖中等溫線可以 看到低溫鑄件向鑄件中心推移。圖6是同一點(diǎn)在不同澆注溫度下的冷卻曲線， 此 點(diǎn)溫度一直在下降，開始時(shí)，周圍溫差相差不大，溫度下降速度較平緩，然后溫 度下降速度增大，由于開始時(shí)，液相較多，下降較快，但隨著凝固的進(jìn)行，固相 越來越多，由于凝固潛熱的作用，溫度下降速度變得較以前時(shí)刻??；四條曲線的 趨勢(shì)是相同的，但澆注溫度越高，需要更長的時(shí)間才能冷卻到同一溫度，這與現(xiàn)實(shí)結(jié)果相同，596C澆注溫

12、度下需要大約200s點(diǎn)（6,6）冷卻到固相線。由圖7可 以看出，因?yàn)楣滔酀摕岬尼尫?，使鑄件冷卻曲線變得緩慢。圖3t=60s時(shí)溫度場(chǎng)d!5吋何祖璋握覘曲圖2初始t=0.5s時(shí)溫度場(chǎng)帕5町富觀蜒梔MEE圖5岡M氐于固相線時(shí)溫度場(chǎng)Mflca15OS對(duì)唱眾璋漁址！1i M不同澆注溫度時(shí)，點(diǎn)的冷卻曲線弓0）1OT1 知2DO時(shí)間丿容250500M M M7 7 6 6 6 6 5 5 5 5 4 4 4 4PS3帕5町富觀蜒梔MEEIp-圖4t=150s時(shí)溫度場(chǎng)低十固相線時(shí)溫度場(chǎng)模擬圖圖6不同澆注溫度時(shí)點(diǎn)（6,6）處冷卻曲線圖7澆注溫度為596 C時(shí)，潛熱對(duì)冷卻過程的影響3結(jié)論1） 通過對(duì)凝固過程相關(guān)因

13、素作出合理假設(shè)，建立了“ L”形AZ91鎂合金鑄 件凝固過程的數(shù)學(xué)模型。模擬結(jié)果與實(shí)際過程基本吻合，可以作為實(shí)際鑄造過程 的分析依據(jù)。2） 利用溫度場(chǎng)的模擬結(jié)果可以預(yù)測(cè)可能出現(xiàn)的缺陷及缺陷出現(xiàn)的部位，為優(yōu) 化鑄造工藝提供了科學(xué)依據(jù)。3）由于凝固潛熱的釋放，會(huì)降低合金的冷卻速度。4）澆注溫度不同，鑄件某點(diǎn)溫度的下降走勢(shì)相同，澆注溫度越高，需要更 長的時(shí)間才能達(dá)到固相線溫度。參考文獻(xiàn)1 李青松.AZ91鑄造鎂合金凝固過程數(shù)值模擬D.華東交通大學(xué)碩士學(xué)位論文.20092 于彥東，馬秋，等.鎂合金壓鑄件凝固過程計(jì)算機(jī)模擬J.熱加工工藝，2005, No.123 周建興,劉瑞祥,等.凝固過程數(shù)值模擬中的

14、潛熱處理方法J.鑄造， 2001,V0l.50,No.74 林首位.鑄件凝固過程三維溫度場(chǎng)數(shù)值模擬研究D.太原:華北工學(xué)院工學(xué)碩士 學(xué)位論文,20015 嚴(yán)力.AZ91鎂合金調(diào)壓鑄造充型能力研究D.西北工業(yè)大學(xué)碩士學(xué)位論文，2006趙彥民定向凝固AZ91鎂合金工藝與組織研究D.太原科技大學(xué)碩士學(xué)位論文，20107許芝卉.凝固過程的溫度場(chǎng)數(shù)值模擬J.數(shù)學(xué)的實(shí)踐與認(rèn)識(shí)，2010,Vol.40,No.13附錄：源程序代碼%程序一：取初始一張、凝固過程二張，結(jié)束一張溫度場(chǎng)分布圖clear,clc,clfL1=40;L2=20;N=80;M=40;% 長為 40cm,80等分，寬為 20cm, 40 等

15、分。T0=596;%鑄件初始溫度a=0.01;%a= Cdt=0.5;%時(shí)間步長0.5sdx=L1/N;dy=L2/M;M仁 a*dt/(dxA2);M2=a*dt/(dyA2);T=T0*o nes(N+1,M+1);T1=T0*o nes(N+1,M+1);t=0;l=0;k=0;Tc=;%采集點(diǎn)(6,6)處溫度值for i=1:81for j=1:41if j=1|j=41|i=1|i=81T(i,j)=25;%邊界點(diǎn)溫度為25TelseT(i,j)=596; % 液相線 T0=596Tendendendfor i=11:81for j=11:41T(i,j)=25;endendif (

16、2*M1+2*M2468)%剛低于固相線時(shí)模擬結(jié)束t=t+dt;k=k+1;for i=2:80for j=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);endendfor i=2:80for j=2:40T(i,j)=T1(i,j);endendif t=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 時(shí)溫度場(chǎng)axis tight;ylabel(y,FontSize

17、,14);xlabel(x,FontSize,14);zlabel(T/C ,FontSize,14)title(0.5s 時(shí)溫度場(chǎng)模擬圖,FontSize,18)subplot(1,2,2);C,H=co ntour(x,y,T(i,j);clabel(C,H);axis square;xlabel(x,Fo ntSize,14);ylabel(y,Fo ntSize,14);title(0.5s 時(shí)模擬等溫線圖,FontSize,18)endif t=60i=1:81;j=1:41;x,y=meshgrid(j,i);figure (2);subplot(1,2,1);mesh(x,y,T

18、(i,j)% 畫出 60s 時(shí)溫度場(chǎng)axis tight;ylabel(y,FontSize,14);xlabel(x,FontSize,14);zlabel(T/C ,FontSize,14) title(60s 時(shí)溫度場(chǎng)模擬圖,FontSize,18)subplot(1,2,2);C,H=co ntour(x,y,T(i,j);clabel(C,H);axis square;xlabel(x,Fo ntSize,14);ylabel(y,Fo ntSize,14);title(60s 時(shí)模擬等溫線圖,FontSize,18)endif t=150i=1:81;j=1:41;x,y=mesh

19、grid(j,i);figure);subplot(1,2,1);mesh(x,y,T(i,j)% 畫出 150s 時(shí)溫度場(chǎng)axis tight;ylabel(y,FontSize,14);xlabel(x,FontSize,14);zlabel(T/C ,FontSize,14)title(150s 時(shí)溫度場(chǎng)模擬圖,FontSize,18)subplot(1,2,2);C,H=co ntour(x,y,T(i,j);clabel(C,H);axis square;xlabel(x,Fo ntSize,14);ylabel(y,Fo ntSize,14);title(150s 時(shí)模擬等溫線圖,

20、FontSize,18)endif(k=2)l=l+1;Tc(l)=T(6,6);%每秒采集一次溫度值k=0;endendi=1:81;j=1:41;x,y=meshgrid(j,i);figure(4);subplot(1,2,1);mesh(x,y,T(i,j)%畫出低于固相線時(shí)溫度場(chǎng)axis tight;ylabel(y,FontSize,14);xlabel(x,FontSize,14);zlabel(T/C ,FontSize,14)title(低于固相線時(shí)溫度場(chǎng)模擬圖,FontSize,18)subplot(1,2,2);C,H=co ntour(x,y,T(i,j);clabel

21、(C,H);axis square;xlabel(x,Fo ntSize,14);ylabel(y,Fo ntSize,14);title(低于固相線時(shí)模擬等溫線圖,FontSize,18)figure(5);xx=1: nu mel(Tc);plot(xx,Tc,k-,li newidth,2)xlabel(時(shí)間 /s,FontSize,14);ylabel(溫度 /C ,FontSize,14);title(點(diǎn)(6,6)的冷卻 曲線,Fo ntSize,18);elsedisp(Error!) %如果不滿足穩(wěn)定性條件，顯示“ Error!;end%程序二：某一在不同澆注溫度下的冷卻曲線cl

22、ear,clc,clffor T0=596:20:656%不同的澆注溫度L1=40;L2=20;N=80;M=40;% 長為 40cm,80 等分，寬為 20cm, 40 等分。 a=0.01;dt=0.5; %時(shí)間步長0.5sdx=L1/ N;dy=L2/ M;M仁 a*dt/(dxA2);M2=a*dt/(dyA2);T=T0*o nes(N+1,M+1);T1=T0*o nes(N+1,M+1);t=0;l=0;k=0;Tc=;%每一秒采集一個(gè)點(diǎn)for i=1:81for j=1:41if j=1|j=41|i=1|i=81T(i,j)=25;%邊界點(diǎn)溫度為25TelseT(i,j)=T

23、0;endendendfor i=11:81for j=11:41T(i,j)=25;endendif (2*M1+2*M2468)t=t+dt;k=k+1;for i=2:80for j=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);endendfor i=2:80for j=2:40T(i,j)=T1(i,j);endendif(k=2)%每秒采集一次(6,6)溫度值1=1+1;Tc(l)=T(6,6);k=0;endendxx=1: nu mel(Tc);plot(xx,Tc,k-,li newidth,2)xlabel(時(shí)間 /s,FontSize,14);ylabel(溫度/C ,FontSize,14);title(不同澆注溫度 時(shí)，點(diǎn)(6,6)的冷卻曲線,FontSize,18);hold onelsedisp(Error!) %如果不滿足穩(wěn)定性條件，顯示“ Error！