計(jì)算機(jī)在材料學(xué)中綜合作業(yè)

計(jì)算機(jī)在材料學(xué)中的應(yīng)用綜合作業(yè)

目錄一、傳熱學(xué)模型建立 -2-（一）【模型一】平板焊接過程的二維溫度場的數(shù)值模擬 -2-（二）【模型二】立方體淬火過程三維非穩(wěn)態(tài)溫度場的有限差分計(jì)算 -5-（三）【模型三】平板激光淬火過程的三維溫度場數(shù)值模擬 -9-（四）【模型四】長方體正火過程三維非穩(wěn)態(tài)溫度場的有限差分計(jì)算 -14-（五）【模型五】端淬實(shí)驗(yàn)的數(shù)值模擬 -18-三、計(jì)算機(jī)模擬計(jì)算·結(jié)果處理 -26-四、團(tuán)隊(duì)工作情況 -29-五、成員投入程度評(píng)定排序 -30-六、總結(jié)與體會(huì) -31-一、傳熱學(xué)模型建立（一）【模型一】平板焊接過程的二維溫度場的數(shù)值模擬1.模型的建立二維模型的假設(shè)：1.材料表面熱物性參數(shù)不隨溫度變化。2.不考慮相變潛熱。3.考慮工件的輻射與空氣的對(duì)流換熱。4.焊接熱源能量分布為高斯分布。5.材料各向同性。6.工件為二維有限大物體。（Z方向厚度不計(jì)）2.熱傳導(dǎo)方程采用直角坐標(biāo)系，則材料內(nèi)部熱傳導(dǎo)方程為：式中，QUOTE，為材料熱擴(kuò)散率，ρ為密度，QUOTE為比熱，λ為熱傳導(dǎo)率，T為溫度，t為時(shí)間。材料的熱物性參數(shù)，比熱QUOTE，熱傳導(dǎo)系數(shù)λ和密度ρ均不隨溫度變化。3.邊界條件上表面：QUOTE其它表面：QUOTE式中，T是工件表面的溫度，Ta是環(huán)境溫度。n是其它表面的外法線方向。h是材料表面總的換熱系數(shù)，包括空氣對(duì)流和熱輻射換熱。h=QUOTE+QUOTE，QUOTE為對(duì)流換熱系數(shù)，QUOTE為輻射換熱系數(shù)。q(x,y,t)為焊接熱源能量分布函數(shù)P為熱源有效功率，QUOTE為熱源半徑，v為熱源運(yùn)動(dòng)速度。4.有限差分方程由于對(duì)稱性，取工件整體區(qū)域的一半求解。X方向節(jié)點(diǎn)數(shù)為N1，Y方向節(jié)點(diǎn)數(shù)為N2。X方向步長為△x，Y方向步長為△y，時(shí)間步長為△t。熱源沿X方向運(yùn)動(dòng)。內(nèi)部節(jié)點(diǎn)差分方程：令QUOTE，QUOTE，QUOTE，則內(nèi)部節(jié)點(diǎn)差分方程變?yōu)椋核臈l線邊界節(jié)點(diǎn)的差分方程：四個(gè)角邊界節(jié)點(diǎn)的差分方程：（二）【模型二】立方體淬火過程三維非穩(wěn)態(tài)溫度場的有限差分計(jì)算1.模型的建立鋼的淬火是將鋼加熱到臨界溫度Ac3（亞共析鋼）或Ac1（過共析鋼）以上溫度，保溫一段時(shí)間，使之全部或部分奧氏體化，然后以大于臨界冷卻速度的冷速快冷到Ms以下（或Ms附近等溫）進(jìn)行馬氏體（或貝氏體）轉(zhuǎn)變的熱處理工藝，，以大幅提高鋼的剛性、硬度、耐磨性、疲勞強(qiáng)度以及韌性等，從而滿足各種機(jī)械零件和工具的不同使用要求。也可以通過淬火滿足某些特種鋼材的鐵磁性、耐蝕性等特殊的物理、化學(xué)性能。淬火過程是一個(gè)涉及到相變，熱傳導(dǎo)，熱對(duì)流，熱輻射的三維非穩(wěn)態(tài)傳熱問題。為計(jì)算淬火過程的溫度場，首先應(yīng)建立淬火過程的傳熱學(xué)模型。本模型中采用以下幾點(diǎn)假設(shè)。材料的熱物性參數(shù)不隨溫度變化不考慮相變潛熱材料為各向同性考慮工件輻射與空氣對(duì)流換熱。工件為立方體2.熱傳導(dǎo)方程根據(jù)工件形狀，采用三維直角坐標(biāo)系，在工件內(nèi)部的熱傳導(dǎo)方程為，式中，，為材料的熱擴(kuò)散率，λ為熱傳導(dǎo)率，ρ為材料密度，為比熱，T為溫度，t為時(shí)間。λ,ρ,均不隨溫度變化。3.邊界條件式中，T為工件表面溫度，為淬火介質(zhì)溫度，n為其他表面的的外法線方向，h為淬火介質(zhì)的對(duì)流換熱系數(shù)。4.初始條件初始時(shí)刻工件整體溫度分布均勻，即，為一常數(shù)，取淬火加熱溫度，對(duì)35CrMo鋼，取860℃。5.有限差分方程首先對(duì)求解區(qū)域離散化，求解區(qū)域?yàn)榱⒎襟w，在X,Y,Z方向用網(wǎng)格來離散求解區(qū)域，X方向節(jié)點(diǎn)數(shù)為N1，Y方向節(jié)點(diǎn)數(shù)為N2，Z方向節(jié)點(diǎn)數(shù)為N3。X方向步長為△X，Y方向步長為△Y，Z方向步長為△Z。時(shí)間步長為△t。內(nèi)部某節(jié)點(diǎn)編號(hào)為（i,j,k).左右相鄰兩節(jié)點(diǎn)為(i,j-1,k),(i,j+1,k),前后相鄰兩節(jié)點(diǎn)為(i.j,k-1),(i,j,k+1).上下相鄰兩節(jié)點(diǎn)為（i,j-1,k),(i,j+1,k)。內(nèi)部節(jié)點(diǎn)差分方程為令，，，。則內(nèi)部節(jié)點(diǎn)差分方程變成了。同理，六個(gè)面邊界節(jié)點(diǎn)的差分方程分別為：十二條棱上邊界節(jié)點(diǎn)的差分方程：八個(gè)角邊界節(jié)點(diǎn)的差分方程為：（三）【模型三】平板激光淬火過程的三維溫度場數(shù)值模擬1.建立傳熱學(xué)模型圖1是平板類工件激光熱處理示意圖。在平板類工件激光熱處理過程的三維傳熱學(xué)模型中采用以下幾點(diǎn)假設(shè)：（1）材料表面對(duì)激光的吸收系數(shù)不隨溫度變化；（2）材料的熱物性參數(shù)隨溫度變化；（3）考慮相變潛熱；（4）考慮工件的輻射與空氣對(duì)流換熱；（5）入射激光束能量分布為高斯分布（TEM00）；（6）激光功率恒定；（7）激光束掃描速度不變；（8）材料各向同性；（9）工件為三維有限大物體；（10）工件初始溫度恒定。圖1平板類工件激光熱處理示2.熱傳導(dǎo)方程P.S.如果激光束掃描的速度很快，則沿激光束掃描方向（Y方向）的傳熱可以忽略，這樣三維傳熱問題就可簡化為X方向和Z方向傳熱的二維傳熱問題，這就是平板類工件激光熱處理過程的二維傳熱模型。在直角坐標(biāo)系中，工件內(nèi)部的熱傳導(dǎo)方程為：式中，ρ為材料密度，Cp為材料比熱，λ為熱傳導(dǎo)系數(shù)，T為溫度，t為時(shí)間。3.邊界條件邊界條件為：上表面：其他表面：Q（x,y,t）是激光光斑能量分布函數(shù)。P為激光功率，A為吸收系數(shù)，R為激光光斑直徑，V為激光光斑運(yùn)動(dòng)速度。Ta是環(huán)境溫度，h是材料表面總的換熱系數(shù)，包括空氣對(duì)流和熱輻射換熱。n是其他表面的外法線方向。4.初始條件初始條件為：初始時(shí)刻工件整體溫度分布均勻，即，Ta為一常數(shù)。5.相變潛熱的處理激光熱處理加熱時(shí)發(fā)生奧氏體相變，冷卻時(shí)發(fā)生馬氏體相變，其他相變不可能發(fā)生。目前對(duì)于在激光相變硬化等超快速加熱條件下的奧氏體轉(zhuǎn)變動(dòng)力學(xué)尚不明白，因此，采用把奧氏體化的相變潛熱折算成溫度降低值的辦法，把溫度降低值等分，分步加到計(jì)算得到的溫度場中。HA為奧氏體的相變潛熱，n為一整數(shù)，一般取為10，或者其他合適的整數(shù)。可以用Koistinen-Marburger公式來計(jì)算鋼鐵材料中馬氏體相變的轉(zhuǎn)變量Fm為馬氏體轉(zhuǎn)變的體積百分?jǐn)?shù)，Ms為馬氏體開始轉(zhuǎn)變溫度，T為節(jié)點(diǎn)溫度，k和n為兩個(gè)常數(shù)，與具體的鋼種有關(guān)。每一步馬氏體相變所放出的相變潛熱可折算成溫升，然后利用ΔTmn對(duì)溫度場進(jìn)行修正：Hm為馬氏體轉(zhuǎn)變的相變潛熱。ΔTmn是每一步馬氏體轉(zhuǎn)變放出的相變潛熱所引起的溫升。有限差分方程（激光熱處理的三維差分方程）由偏微分方程和邊界條件，初始條件構(gòu)成傳熱的定解問題。用能量平衡法推導(dǎo)了這一問題的非均勻空間網(wǎng)絡(luò)的三維顯式差分方程。直角坐標(biāo)系下差分方程：X方向節(jié)點(diǎn)數(shù)為101，y方向節(jié)點(diǎn)數(shù)為101，z方向節(jié)點(diǎn)數(shù)為51。X方向步長為△X，y方向步長為△Y，z方向步長為△Z。內(nèi)部節(jié)點(diǎn)差分方程為令，，，。則內(nèi)部節(jié)點(diǎn)差分方程變成了。六個(gè)面上節(jié)點(diǎn)的差分方程分別為：上表面：下表面：左表面：右表面：前表面：后表面：棱上的節(jié)點(diǎn)的差分方程：八個(gè)角上節(jié)點(diǎn)的差分方程為：（四）【模型四】長方體正火過程三維非穩(wěn)態(tài)溫度場的有限差分計(jì)算1、模型的建立鋼鐵材料在爐中加熱到超過奧氏體轉(zhuǎn)變點(diǎn)后，發(fā)生奧氏體轉(zhuǎn)變，然后放到空氣中冷卻，就可得到含有珠光體的均勻組織，從而實(shí)現(xiàn)鋼鐵材料的正火處理。長方體是一種在工業(yè)界中常用的工件，他的正火過程是一個(gè)涉及到相變、熱傳導(dǎo)、熱對(duì)流的三維非穩(wěn)態(tài)傳熱過程，可用三維模型近似的描述長方體正火冷卻的傳熱過程，在這個(gè)三維模型中采用以下幾個(gè)假設(shè)：工件為三維長方體材料的熱物性參數(shù)不隨溫度變化不考慮相變潛熱考慮工件與空氣的熱對(duì)流換熱材料各向同性2、直角坐標(biāo)系的建立及網(wǎng)格的劃分根據(jù)工件的形狀，建立如圖1所示的坐標(biāo)系，將OX，OY分成100份，OZ分成50份，即節(jié)點(diǎn)為從（i,j,k)=(1,1,1)到（i,j,k)=(101,101,51)。此模型為長×寬×高等于100×100×50（cm）的鋼鐵材料，材料的密度為7.8g/,比熱為0.5J/(g.k),導(dǎo)熱系數(shù)為0.3。工件冷卻過程為從邊界向中心逐漸冷卻。3、熱傳導(dǎo)方程在工件內(nèi)部的熱傳導(dǎo)方程為，式中，，為材料的熱擴(kuò)散率，λ為熱傳導(dǎo)率，ρ為材料密度，為比熱，T為溫度，t為時(shí)間。λ,ρ,均不隨溫度變化。初始條件初始時(shí)刻工件整體溫度分布均勻，即=，設(shè)為900°C。設(shè)環(huán)境溫度=20°C。5、邊界條件外表面：式中，T為工件表面溫度，為環(huán)境溫度，n為外表面的外法線方向，h為材料表面的換熱系數(shù)。h=+,為對(duì)流換熱系數(shù)，為輻射換熱系數(shù)。6、有限差分方程X方向節(jié)點(diǎn)數(shù)為101，y方向節(jié)點(diǎn)數(shù)為101，z方向節(jié)點(diǎn)數(shù)為51。X方向步長為△X，y方向步長為△Y，z方向步長為△Z。內(nèi)部節(jié)點(diǎn)差分方程為令，，，。則內(nèi)部節(jié)點(diǎn)差分方程變成了六個(gè)面上節(jié)點(diǎn)的差分方程分別為：上表面：下表面：左表面：右表面：前表面：后表面：十二條棱上的節(jié)點(diǎn)的差分方程：八個(gè)頂角上節(jié)點(diǎn)的差分方程為：（五）【模型五】端淬實(shí)驗(yàn)的數(shù)值模擬1、模型的建立端淬實(shí)驗(yàn)是將標(biāo)準(zhǔn)尺寸的端淬試樣(Φ25mm×100mm)奧氏體化后，在專用設(shè)備上對(duì)其一端面噴水冷卻，后沿軸線方向的測出硬度-距水冷端距離的關(guān)系曲線的試驗(yàn)方法，是測定鋼的淬透性的方法之一。由于模型軸向?qū)ΨQ，計(jì)算時(shí)可以采用二維模型如圖。此次模擬實(shí)驗(yàn)的模型為標(biāo)準(zhǔn)尺寸（Φ25mm×100mm）的金屬試樣。先將式樣加熱到860℃并保溫，使式樣完全奧氏體化，然后將式樣取出，放在端淬試驗(yàn)機(jī)上，在式樣下端噴20℃的水120s。模型中采取如下假設(shè)：材料的熱物性參數(shù)不隨溫度變化忽略相變潛熱對(duì)溫度場的影響不考慮軸向傳熱材料為各向同性材料考慮材料與空氣的對(duì)流換熱，采取第三類邊界條件網(wǎng)格劃分取式樣延一個(gè)直徑的剖面的右半部分建立如上右圖的直角坐標(biāo)系，將or方向分成125份，oz方向分成1000份。[即從（i，k)=（1,1）到（i，k）=（N1，N2）=（126,1001）]。3.熱傳導(dǎo)方程材料內(nèi)部的熱傳導(dǎo)方程是，式中為密度，為比熱，為傳導(dǎo)熱，T為溫度，t為時(shí)間。4.邊界條件下邊：；左邊：0；其他邊：式中，T是工件表面溫度，Ta是環(huán)境溫度（水溫與環(huán)境溫度相圖）。n是其他表面的法線方向。是式樣與冷卻水的對(duì)流換熱系數(shù)，是式樣與空氣的換熱系數(shù)。材料的熱物性參數(shù)，比熱，熱傳導(dǎo)系數(shù)和密度均不隨溫度的變化而變化。初始條件端淬開始時(shí)，式樣溫度為860℃且溫度均勻，完全奧氏體化。冷卻水溫度與環(huán)境溫度始終為20℃。有限差分方程r方向的節(jié)點(diǎn)數(shù)為31，Z方向的節(jié)點(diǎn)數(shù)為201，r方向的步長為，Z方向的步長為。內(nèi)部節(jié)點(diǎn)的差分方程為令，，則內(nèi)部節(jié)點(diǎn)的差分方程變?yōu)榱顒t：四個(gè)邊上的節(jié)點(diǎn)的差分方程：上邊：下邊：左邊：右邊：四個(gè)角上節(jié)點(diǎn)的差分方程：右上角：右下角：左上角：左下角：二、計(jì)算機(jī)編程模擬·程序設(shè)計(jì)程序如下：/****有限長圓柱體端淬****//*****初始溫度860℃******//****采用顯式差分方程****/#include<stdio.h>#include<math.h>#definePI3.14#defineK0.3//導(dǎo)熱系數(shù)#defineCP0.5//比熱#defineP7.8//密度#defineH12.0//與淬火介質(zhì)的換熱系數(shù)#defineH20.0003//與空氣的換熱系數(shù)#defineN151#defineN2401doubleS,DS=0.001,S1;//時(shí)間與時(shí)間步長intT0=860,Ta=20;//初始溫度doubleT1[N1][N2],T2[N1][N2];doubleFR,Fz,F1,F2;//傅里葉數(shù)doubleDR=0.025,DZ=0.025;//空間步長voidf_number()//傅里葉數(shù)的計(jì)算{ 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()//熱傳導(dǎo)方程{ inti,j; /*內(nèi)部節(jié)點(diǎn)*/ 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]); } /*邊上的節(jié)點(diǎn)*/ /****上邊****/ 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]); } /*角上的節(jié)點(diǎn)*/ /***右上角***/ 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("請(qǐng)輸入淬火時(shí)間:");scanf("%lf",&S1); printf("\n");f_number(); printf("請(qǐng)判斷傅里葉數(shù)，輸入任意數(shù)字加回車?yán)^續(xù)\n"); scanf("%d",&p); printf("結(jié)果寫入文件中，請(qǐng)稍后\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; }三、計(jì)算機(jī)模擬計(jì)算·結(jié)果處理不同時(shí)刻端淬試樣的溫度場情況50s 100s150s200s250s300s溫度標(biāo)尺從以上圖中可以看出，水冷端開始時(shí)溫度降低很快，降低到100℃后開始減慢；試樣上端開始時(shí)溫度變化很慢。整個(gè)的端淬過程中，水冷端溫度迅速下降，低溫區(qū)逐漸向上移動(dòng)。端淬過程的時(shí)間——溫度圖及分析材料的CCT曲線由cct曲線可知，材料的鼻尖溫度大約為550℃，時(shí)間為1s。對(duì)比上面的端淬過程時(shí)間溫度曲線可知，距水冷端2mm處1s時(shí)溫度大約為550℃由此可知，材料可淬成馬氏體的深度大約在2