計(jì)算材料學(xué)Fortran程序匯總

INTEGERM(1:10000),NUMBER1(0:364),NUMBER2REALX,YISEED=RTC()DOJ=1,10000NUMBER1=0X=RAN(ISEED)NUMBER1(0)=INT(365*X+1)JJJ=1DOI=1,364Y=RAN(ISEED)NUMBER2=INT(365*Y+1)ETR=COUNT(NUMBER1.EQ.NUMBER2)IF(ETR==1)THENEXITELSEJJJ=JJJ+1M(J)=JJJNUMBER1(I)=NUMBER2ENDIFENDDOENDDODOI=1,10000IF(M(I).LE.23)SUM=SUM+1ENDDOPRINT*,SUM/10000ENDMonteCarloSimulationofOneDimensionalDiffusionINTEGERX,XX(1:1000,1:1000)REALXXM(1:1000)! X:INSTANTANEOUSPOSITIONOFATOM!XX(J,I)：X*X,J:第幾天實(shí)驗(yàn),I:第幾步跳躍! XXM(I):THEMEANOFXXWRITE(*,*)"實(shí)驗(yàn)天數(shù)JMAX，實(shí)驗(yàn)次數(shù)IMAX"READ(*,*)JMAX,IMAXISEED=RTC()DOJ=1,JMAX!第幾天實(shí)驗(yàn)X=0!!!DOI=1,IMAX!第幾步跳躍RN=RAN(ISEED)IF(RN<0.5)THENX=X+1ELSEX=X-1ENDIFXX(J,I)=X*XENDDOENDDOOPEN(1,FILE="C:\DIF1.DAT")DOI=1,IMAXXXM=0.0XXM(I)=1.0*SUM(XX(1:JMAX,I))/JMAX!!WRITE(1,*)I,XXM(I)ENDDOCLOSE(1)ENDMonteCarloSimulationofTwoDimensionalDiffusionINTEGERX,Y,XY(1:1000,1:1000)REALXYM(1:1000)X:INSTANTANEOUSPOSITIONOFATOM! XY(J,I)：X*Y,J:第幾天實(shí)驗(yàn),I:第幾步跳躍XYM(I):THEMEANOFXYWRITE(*,*)"實(shí)驗(yàn)天數(shù)JMAX，實(shí)驗(yàn)次數(shù)IMAX"READ(*,*)JMAX,IMAXISEED=RTC()DOJ=1,JMAX!第幾天實(shí)驗(yàn)X=0!!!Y=0!!!DOI=1,IMAX!第幾步跳躍RN=RAN(ISEED)IF(RN.LT.0.25)THENx=xy=y-1ENDIFIF(RN.LT.0.5.AND.RN.GE.0.25)THENx=xy=y+1ENDIFIF(RN.LT.0.75.AND.RN.GE.0.5)THENx=x-1y=yENDIFIF(RN.GE.0.75)THENx=x+1y=yENDIFXY(J,I)=X*X+Y*YENDDOENDDOOPEN(1,FILE="C:\DIF2.DAT")DOI=1,IMAXXYM=0.0XYM(I)=1.0*SUM(XY(1:JMAX,I))/JMAX!!WRITE(1,*)I,XYM(I)ENDDOCLOSE(1)END! MonteCarloSimulationofOneDimensionalDiffusionINTEGERX,XY(1:1000,1:1000),y,XN(1:4),YN(1:4),RNREALXYM(1:1000)! X:INSTANTANEOUSPOSITIONOFATOM! XY(J,I)：X*Y,J:第幾天實(shí)驗(yàn),I:第幾步跳躍!XYM(I):THEMEANOFXYWRITE(*,*)"實(shí)驗(yàn)天數(shù)JMAX，實(shí)驗(yàn)次數(shù)IMAX"READ(*,*)JMAX,IMAXXN=(/0,0,-1,1/)YN=(/-1,1,0,0/)ISEED=RTC()DOJ=1,JMAX!第幾天實(shí)驗(yàn)X=0!!!Y=0!!!DOI=1,IMAX!第幾步跳躍RN=4*RAN(ISEED)+1X=X+YN(RN)Y=Y+YN(RN)XY(J,I)=X*X+Y*YENDDOENDDOOPEN(1,FILE="C:\DIF2.DAT")DOI=1,IMAXXYM=0.0XYM(I)=1.0*SUM(XY(1:JMAX,I))/JMAX!!WRITE(1,*)I,XYM(I)ENDDOCLOSE(1)END做三維空間隨機(jī)行走？？留作業(yè)! MonteCarloSimulationofOneDimensionalDiffusionINTEGERX,XY(1:1000,1:1000),y,XN(1:6),YN(1:6),ZN(1:6),RNREALXYM(1:1000)! X:INSTANTANEOUSPOSITIONOFATOM! XY(J,I)：X*Y,J:第幾天實(shí)驗(yàn),I:第幾步跳躍!XYM(I):THEMEANOFXYWRITE(*,*)"實(shí)驗(yàn)天數(shù)JMAX，實(shí)驗(yàn)次數(shù)IMAX"READ(*,*)JMAX,IMAXXN=(/0,0,-1,1,0,0/)YN=(/-1,1,0,0,0,0/)ZN=(/0,0,0,0,1,-1/)ISEED=RTC()DOJ=1,JMAX!第幾天實(shí)驗(yàn)X=0!!!Y=0!!!Z=0DOI=1,IMAX!第幾步跳躍RN=6*RAN(ISEED)+1X=X+XN(RN)Y=Y+YN(RN)Z=Z+ZN(RN)XY(J,I)=X*X+Y*Y+Z*ZENDDOENDDOOPEN(1,FILE="C:\DIF2.DAT")DOI=1,IMAXXYM=0.0XYM(I)=1.0*SUM(XY(1:JMAX,I))/JMAX!!WRITE(1,*)I,XYM(I)ENDDOCLOSE(1)END通過該程序了解fortran語言如何畫圖(通過像素畫圖)USEMSFLIBINTEGERXR，YR!在的區(qū)域中畫一個(gè)圓PARAMETERXR=400,YR=400INTEGERR,S(1:XR,1:YR)X0=XR/2!圓心位置X0，YOY0=YR/2R=MIN(X0-10,Y0-10)!圓半徑S=0!像素的初始狀態(tài)(顏色)DOI=1,XRDOJ=1,YRIF((I-X0)**2+(J-Y0)**2<=R**2)S(I,J)=10IER=SETCOLOR(S(I,J))IER=SETPIXEL(I,J)ENDDOENDDOEND畫一個(gè)圓(1、如何選出晶界區(qū)域；2、進(jìn)一步加深對(duì)畫圖的理解)USEMSFLIBINTEGERXR，YR!在的區(qū)域中畫一個(gè)圓PARAMETERXR=400,YR=400INTEGERR,S(0:XR+1,0:YR+1),XN(1:4),YN(1:4),SNSXN=(/0,0,-1,1/)YN=(/-1,1,0,0/)X0=XR/2!圓心位置X0，Y0Y0=YR/2R=MIN(X0-10,Y0-10)!圓半徑S=0!像素的初始狀態(tài)(顏色)DOI=1,XRDOJ=1,YRIF((I-X0)**2+(J-Y0)**2<=R**2)S(I,J)=10IER=SETCOLOR(S(I,J))IER=SETPIXEL(I,J)ENDDOENDDODOI=1,XR!畫晶界DOJ=1,YRNDS=0DOK=1,4IF(S(I,J).NE.S(I+XN(K),J+YN(K)))NDS=NDS+1ENDDOIF(NDS>0)THENIER=SETCOLOR(9)ELSEIER=SETCOLOR(8)ENDIFIER=SETPIXEL(I,J)ENDDOENDDOEND如何畫有一定寬度的晶界？！MC模擬一個(gè)晶粒的縮小(1、 如何定義基體和晶粒；(2、 如何尋找邊界；(3、 如何計(jì)算能量(構(gòu)造或描述問題的概率過程)步驟：1、在基體上畫一個(gè)原始晶粒，并賦予基體和原始晶粒不同的狀態(tài)(取向號(hào)或體積能)2、尋找晶界3、 MC的一個(gè)時(shí)間步(晶粒長(zhǎng)大過程中一種隨機(jī)性)4、 計(jì)算晶界上網(wǎng)格點(diǎn)的相互作用，每個(gè)格點(diǎn)的相互作用導(dǎo)致晶粒長(zhǎng)大隨機(jī)選取一個(gè)初始格點(diǎn)；若此點(diǎn)屬于晶界，那么可以隨機(jī)轉(zhuǎn)變?yōu)樗従拥臓顟B(tài)，若不是晶界，則跳出循環(huán)，不發(fā)生轉(zhuǎn)變；計(jì)算轉(zhuǎn)變前后的能量變化，/E=/E+/E+/Evsq(自由能=體積能+表面能+能量起伏)若/E小于0則新晶相被接受，網(wǎng)格狀態(tài)發(fā)生轉(zhuǎn)變。USEMSFLIBPARAMETERIR=400,JR=400INTEGERIS(0:IR+1,0:JR+1),TMAX,ISN(1:8),NSTATE,T,NR,IX,IYWRITE(*,*)"PLEASEINPUTTHETIMESTEP"READ(*,*)TMAXISEED=RTC()定義圓心和半徑IRC=IR/2JRC=JR/2R=MIN(IRC,JRC)-10!定義基體和圓晶粒分別為狀態(tài)1、狀態(tài)2IS=1DOI=1,IRDOJ=1,JRDISTANCE=SQRT(1.0*(I-IRC)**2+1.0*(J-JRC)**2)IF(DISTANCE.LT.R)IS(I,J)=2ISE=SETCOLOR(IS(I,J))ISE=SETPIXEL(I,J)ENDDOENDDOOPEN(1,FILE="E:\LUKE.DAT")!尋找晶粒邊界，計(jì)算能量，改變狀態(tài)。DOT=1,TMAXDOX=1,IRDOY=1,JRIX=IR*RAN(ISEED)+1JY=JR*RAN(ISEED)+1ISN=(/IS(IX-1,JY-1),IS(IX-1,JY),IS(IX-1,JY+1),IS(IX,JY-1),IS(IX,JY+1),IS(IX+1,JY-1),IS(IX+1,JY),IS(IX+1,JY+1)/)E0=COUNT(ISN.NE.IS(IX,JY))!如果不是圓晶粒邊界，則跳出，重新循環(huán)IF(E0.EQ.0)CYCLE!隨機(jī)尋找一個(gè)相鄰點(diǎn)NR=8*RAN(ISEED)+1NSTATE=ISN(NR)!判斷與相鄰點(diǎn)的能量差，并決定是否改變狀態(tài)E=COUNT(ISN.NE.NSTATE)RD=RAN(ISEED)DE=E-E0+NSTATE-IS(IX,JY)+2.5*RD-1.25IF(DE.LT.0.0)IS(IX,JY)=NSTATEISRE=SETCOLOR(IS(IX,JY))ISRE=SETPIXEL(IX,JY)ENDDOENDDOWRITE(1,*)T,COUNT(IS.EQ.2)ENDDOCLOSE(1)END(1、 演示體積能互換的轉(zhuǎn)變情況；2、演示IX=IR*RAN(ISEED)+1JY=JR*RAN(ISEED)+1的替換；3、能量關(guān)系變化；)作業(yè)3：基體上單晶粒形核長(zhǎng)大(1、 選取初始形核點(diǎn)(2、 定義體積能的變化)USEMSFLIBPARAMETERIR=400,JR=400INTEGERIS(0:IR+1,0:JR+1),TMAX,ISN(1:8),NSTATE,T,NR,IX,IYWRITE(*,*)"PLEASEINPUTTHETIMESTEP"READ(*,*)TMAXISEED=RTC()!定義圓心和半徑IRC=IR/2JRC=IR/2!定義基體和圓晶粒分別為狀態(tài)10、狀態(tài)2IS=10IS(IRC,JRC)=2OPEN(1,FILE="E:\LUKE.DAT")!尋找晶粒邊界，計(jì)算能量，改變狀態(tài)。DOT=1,TMAXDOX=1,IRDOY=1,JRIX=IR*RAN(ISEED)+1JY=JR*RAN(ISEED)+1ISN=(/IS(IX-1,JY-1),IS(IX-1,JY),IS(IX-1,JY+1),IS(IX,JY-1),IS(IX,JY+1),IS(IX+1,JY-1),IS(IX+1,JY),IS(IX+1,JY+1)/)E0=COUNT(ISN.NE.IS(IX,JY))!如果不是圓晶粒邊界，則跳出，重新循環(huán)IF(E0.EQ.0)CYCLE!隨機(jī)尋找一個(gè)相鄰點(diǎn)NR=8*RAN(ISEED)+1NSTATE=ISN(NR)!判斷與相鄰點(diǎn)的能量差，并決定是否改變狀態(tài)基體上多晶粒形核長(zhǎng)大，（1基體上多晶粒形核長(zhǎng)大，（1、 如何生成多晶粒形核核心（2、 如何定義、并區(qū)別多晶粒3、為什么要有邊界條件）1步驟：1、定義變量：體積能變量（一維數(shù)組）；基體各像素點(diǎn)變量（二維數(shù)組），鄰居的取向號(hào)（一維數(shù)組）等2、讀取晶粒生長(zhǎng)步長(zhǎng)；3、尋找晶界4、計(jì)算晶界上網(wǎng)格點(diǎn)的相互作用，每個(gè)格點(diǎn)的相互作用導(dǎo)致晶粒長(zhǎng)大或縮?。孩匐S機(jī)選取一個(gè)初始格點(diǎn)；②若此點(diǎn)屬于晶界，那么可以隨機(jī)轉(zhuǎn)變?yōu)樗従拥臓顟B(tài)，若不是晶界，則跳出循環(huán)，不發(fā)生轉(zhuǎn)變；③計(jì)算轉(zhuǎn)變前后的能量變化，/E=/E+/E+/EvsqE=COUNT(ISN.NE.NSTATE)RD=RAN(ISEED)DE=E-E0+NSTATE-IS(IX,JY)+2.5*RD-1.25IF(DE.LT.0.0)IS(IX,JY)=NSTATEISRE=SETCOLOR(IS(IX,JY))ISRE=SETPIXEL(IX,JY)ENDDOENDDOWRITE(1,*)T,SQRT(1.0*COUNT(IS.EQ.2))ENDDOCLOSE(1)END（自由能=體積能+表面能+能量起伏）④若/E小于0則新晶相被接受，網(wǎng)格狀態(tài)發(fā)生轉(zhuǎn)變。USEMSFLIB

定義常數(shù)名PARAMETERIR=400,JR=400,NMAX=100定義變量：體積能變量(一維數(shù)組)；基體各像素點(diǎn)變量(二維數(shù)組)，鄰居的取向號(hào)(一維數(shù)組)等INTEGERIS(0:IR+1,0:JR+1),TMAX,ISN(1:8),NSTATE,T,NR,IX0,IY0,IX,IYINTEGERIGV(0:NMAX)讀取晶粒生長(zhǎng)步長(zhǎng)；WRITE(*,*)"PLEASEINPUTTHETIMESTEP"READ(*,*)TMAXISEED=RTC()定義基體體積能為10,所有晶粒體積能為1IGV=1IGV(0)=10定義晶粒長(zhǎng)大位置(100個(gè)形核點(diǎn)初始形核位置)，并附已不同的取向號(hào)DOI=1,NMAXIX0=IR*RAN(ISEED)+1JY0=JR*RAN(ISEED)+1IS(IX0,JY0)=IENDDO邊界條件IS(0,1:JMAX)=IS(IMAX,1:JMAX)IS(IMAX+1,1:JMAX)=IS(1,1:JMAX)IS(0:IMAX+1,0)=IS(0:IMAX+1,JMAX)IS(0:IMAX+1,JMAX+1)=IS(0:IMAX+1,1)OPEN(1,FILE="E:\LUKE.DAT")!尋找晶粒邊界。DOT=1,TMAXDOX=1,IRDOY=1,JRIX=IR*RAN(ISEED)+1JY=JR*RAN(ISEED)+1ISN=(/IS(IX-1,JY-1),IS(IX-1,JY),IS(IX-1,JY+1),IS(IX,JY-1)&,IS(IX,JY+1),IS(IX+1,JY-1),IS(IX+1,JY),IS(IX+1,JY+1)/)E0=COUNT(ISN.NE.IS(IX,JY))!如果不是晶粒邊界，則跳出，重新循環(huán)IF(E0.EQ.0)CYCLE!隨機(jī)尋找一個(gè)相鄰點(diǎn)NR=8*RAN(ISEED)+1NSTATE=ISN(NR)!判斷與相鄰點(diǎn)的能量差，并決定是否改變狀態(tài)RD=RAN(ISEED)E=COUNT(ISN.NE.NSTATE)IG=IGV(NSTATE)-IGV(IS(IX,JY))DE=E-E0+IG+2.5*RD-1.25IF(DE.LT.0.0)IS(IX,JY)=NSTATEISRE=SETCOLOR(MOD(IS(IX,JY),15))ISRE=SETPIXEL(IX,JY)ENDDOENDDOWRITE(1,*)T,1.0*COUNT(IS.NE.0)ENDDOCLOSE(1)END! iii=MOD(IS(IX,JY),15)! If(iii==0)iii=iii+1物理意義1晶粒長(zhǎng)大與時(shí)間的關(guān)系2形核數(shù)目與的關(guān)系3、與EBSD比較4、為什么圓弧會(huì)長(zhǎng)成近似的直線Monte-Carlo方法模擬晶粒長(zhǎng)大過程中的新研究進(jìn)展第一個(gè)問題：采用Monte-Carlo方法存在以下幾方面的問題（1）模擬過程只計(jì)及最近鄰格點(diǎn)的相互作用，難于考慮畸變場(chǎng)等物理場(chǎng)以及晶界形態(tài)對(duì)晶粒長(zhǎng)大過程的影響，不利于在唯象層次上深入開展對(duì)晶粒長(zhǎng)大過程的研究。應(yīng)當(dāng)指出，大部分的計(jì)算機(jī)模擬主要集中于二維系統(tǒng)，雖然能夠在一定程度上提供晶粒長(zhǎng)大的有關(guān)動(dòng)力學(xué)和拓?fù)鋵W(xué)的信息，但是，任何多晶體材料中的晶粒都是三維的，所有二維晶粒長(zhǎng)大的模擬都是基于某種特定的假設(shè)或簡(jiǎn)化，因此模擬三維晶粒的長(zhǎng)大更具有實(shí)際意義。同時(shí)，在目前的研究中，對(duì)晶粒生長(zhǎng)過程的形貌演化模擬的逼真度比較好，但在生長(zhǎng)模型的邊界處理等方面還有待完善，特別是晶粒生長(zhǎng)與晶化溫度、時(shí)間、氣氛等參數(shù)密切相關(guān)，這種復(fù)雜工藝條件下的晶粒生長(zhǎng)過程的模擬是當(dāng)前該領(lǐng)域正待解決的難題。另外，國(guó)內(nèi)外用各種方法模擬晶粒長(zhǎng)大，主要集中在對(duì)方法本身的討論、模擬晶粒分布函數(shù)特點(diǎn)及晶粒形貌等方面，而將模擬方法與實(shí)際工藝聯(lián)系起來的研究者比較少。第二個(gè)問題：算法種類：1、標(biāo)準(zhǔn)算法2、改進(jìn)型算法3、其它算法第三個(gè)問題：模擬晶粒長(zhǎng)大的一些應(yīng)用1、 金屬凝固過程的形核；2、 再結(jié)晶及存在織構(gòu)的再結(jié)晶（形核位置時(shí)考慮能量關(guān)系）；3、 第二相粒子的影響（第二相粒子對(duì)晶粒長(zhǎng)大有阻礙作用.這種情況下模擬較復(fù)雜.由于無法用解析解定量描述,用MC方法模擬顯得尤為重要；(a)j--H-ql*6000$j-i-IWOMCSgf—2.5l1ODD0Mc$ (d)/-S-ft—i2&afl-MC$SI<5不同暉粗時(shí)創(chuàng)址二粕粒子時(shí)兄粧就觀胡級(jí)的匪噸儼I5、異常晶粒的長(zhǎng)大；6、三維晶粒的長(zhǎng)大(模擬不同的點(diǎn)陣結(jié)構(gòu)生長(zhǎng)過程)；山東大學(xué)中國(guó)有色金屬學(xué)報(bào)，2008，4252()105MCS圈斗不同楝型欖擬得到閭晶粒辰人指數(shù)Flg+252()105MCS圈斗不同楝型欖擬得到閭晶粒辰人指數(shù)Flg+4Grfliugi'owthexponentfb]Ldififei'enrniodels?一W)3I、0.4舫±6024°—Model2.0391=0.010*—Model\0.501±0.011200 400 600 800 10007、蒙特卡洛方法不是數(shù)學(xué)的方法，而是一種多次模擬實(shí)驗(yàn)的方法。元胞自動(dòng)機(jī)程序(生命永不停止)usemsflibPARAMETERIR=400,JR=400,NMAX=40000!NMAX隨機(jī)產(chǎn)生的生命種子integeris(0:1001,0:1001),is1(0:1001,0:1001),ISN(1:8),Tmax,num!is-基體的二維數(shù)組print*,'pleaseinputloop(Tmax)'read*,TmaxISEED=RTC()is=15！“死”的狀態(tài)，基體為白色！賦予生命的種子，“活”的狀態(tài)1DOI=1,NMAXIX0=IR*RAN(ISEED)+1JY0=JR*RAN(ISEED)+1IS(IX0,JY0)=1ENDDO!executetheruledot=1,Tmaxis1=is！邊界條件IS(0,1:JMAX)=IS(IMAX,1:JMAX)IS(IMAX+1,1:JMAX)=IS(1,1:JMAX)IS(0:IMAX+1,0)=IS(0:IMAX+1,JMAX)IS(0:IMAX+1,JMAX+1)=IS(0:IMAX+1,1)！搜索生命存在的位置doix=1,IRdojy=1,JR！判斷鄰居狀態(tài)ISN=(/IS(IX-1,JY-1),IS(IX-1,JY),IS(IX-1,JY+1),IS(IX,JY-1)&,IS(IX,JY+1),IS(IX+1,JY-1),IS(IX+1,JY),IS(IX+1,JY+1)/)num=COUNT(ISN.eq.1)！賦予生存的條件if((is(ix,jy)==15.and.num==3).or.(is(ix,jy)==1.and.(num==3.or.num==4)))thenis1(ix,jy)=1elseis1(ix,jy)=15endif！畫圖ISRE=SETCOLOR(IS1(IX,JY))ISRE=SETPIXEL(IX,JY)enddoenddois=is1enddoend元胞自動(dòng)機(jī)—產(chǎn)生和發(fā)展.四個(gè)階段:1940s誕生：VonNeumann自我復(fù)制機(jī).1960-70S起步：JH?Conway生命游戲.1980s理論研究：S.WolframCA分類.1980-90s應(yīng)用：HPP-FHP格子氣自動(dòng)機(jī)?C?LangtonN?Packard人工生命元胞自動(dòng)機(jī)（CellularAutomata,簡(jiǎn)稱CA,也有人稱其為細(xì)胞自動(dòng)機(jī)、點(diǎn)格自動(dòng)機(jī)、分子自動(dòng)機(jī)或單元自動(dòng)機(jī)）是一種建立在離散的時(shí)間和空間上的動(dòng)力學(xué)系統(tǒng)。散布在規(guī)則格網(wǎng)（LatticeGrid）中的每一元胞（Cell）取有限的離散狀態(tài)，遵循同樣的作用規(guī)則，依據(jù)確定的局部規(guī)則作同步更新。大量元胞通過簡(jiǎn)單的相互作用而構(gòu)成動(dòng)態(tài)系統(tǒng)的演化。與一般的動(dòng)力學(xué)模型不同，元胞自動(dòng)機(jī)不是由嚴(yán)格定義的物理方程或函數(shù)確定，而是由一系列模型構(gòu)造的規(guī)則構(gòu)成，凡是滿足這些規(guī)則的模型都可以算做是元胞自動(dòng)機(jī)模型。因此，元胞自動(dòng)機(jī)是一類模型的總稱，或者是一個(gè)方法框架。其特點(diǎn)是時(shí)間、空間、狀態(tài)都是離散的，每個(gè)變量只取有限多個(gè)狀態(tài)，且其狀態(tài)改變的規(guī)則在時(shí)間和空間上都是局部的。元胞自動(dòng)機(jī)是一種對(duì)具有局域連通性的格點(diǎn)，應(yīng)用局部（有時(shí)為中等范圍）確定性或概論性的轉(zhuǎn)換規(guī)則來描述在離散空間和時(shí)間上復(fù)雜系統(tǒng)演化規(guī)律的同步算法。元胞自動(dòng)機(jī)可以采取任意的維數(shù)，自動(dòng)機(jī)的格點(diǎn)描述被認(rèn)為與所選模型密切相關(guān)的基本的系統(tǒng)實(shí)體。獨(dú)立的格點(diǎn)可以表示連續(xù)的體積單元，原子顆粒大原子團(tuán)、汽車、人口、晶格缺陷等。通過作用于每個(gè)格點(diǎn)的確定性或概率性轉(zhuǎn)換規(guī)則，自動(dòng)機(jī)的演化得以實(shí)施。這些規(guī)則決定了每個(gè)元胞的狀態(tài)，并將其描述為其前一狀態(tài)和鄰近元胞的函數(shù)。在計(jì)算狀態(tài)轉(zhuǎn)變過程中所考慮的鄰近位置的數(shù)目確定了相互作用的范圍和自動(dòng)機(jī)的局部演化。！轉(zhuǎn)變概率元胞自動(dòng)機(jī)特點(diǎn)：時(shí)間、空間、狀態(tài)都離散的動(dòng)力學(xué)模型，利用確定性或概論性的轉(zhuǎn)換規(guī)則來描述在離散空間和時(shí)間上復(fù)雜系統(tǒng)演化2_4.2蘭頓的螞蟻唱蟻規(guī)則元胞自動(dòng)機(jī)是由ChrisLangton利GregTurk發(fā)現(xiàn)的叫滋.自動(dòng)札過過極簡(jiǎn)單的兩條規(guī)則來模擬螞蟻的行為規(guī)律，嗎蟻在方形網(wǎng)格上遠(yuǎn)動(dòng)，網(wǎng)格分為黑色和白色，演化規(guī)則対’D如果碼蟻當(dāng)前處于黑色網(wǎng)格」則向左轉(zhuǎn)9臚.并將網(wǎng)格涂成H色.2｝細(xì)果嶋蟻芻前處于白儀網(wǎng)格，則向啟轉(zhuǎn)9仇并將網(wǎng)格涂成黒色.結(jié)昊發(fā)現(xiàn)*這只小媽蟻表現(xiàn)出了用當(dāng)復(fù)雜的著為"從H色的元胞窣間開始’經(jīng)過入綃500步厶胡“它就繞回碾位，之后又進(jìn)入-種混沌狀態(tài)’這種無規(guī)則運(yùn)動(dòng)的狀態(tài)-畝持縝到大約1M00步后，螞蟻竇然又出現(xiàn)非常堀則的運(yùn)動(dòng)形式，它會(huì)F辟一條笙自的路線井遷旦前逵一2.6中臨位卜（b）分別表現(xiàn)了處于混沱狀態(tài)和規(guī)則運(yùn)動(dòng)狀態(tài)的螞蟻?zhàn)詣?dòng)機(jī).（b）圖中螞蟻開辟的筆直的路線被稱為公路⑶】,與網(wǎng)格方向成4和角，公路在不遇到障礙物的情捏下會(huì)一直延伸下喪一閨3第391步,,一個(gè)鵬蟻圖5笫10500歸一個(gè)訥蟻圖4第冗00步』一個(gè)螞蟻3.2多個(gè)螞蟻的情況這里対心網(wǎng)格屮加入名個(gè)螞蟻的情況進(jìn)彳亍了模擬。圖6是件網(wǎng)楙屮加入9個(gè)螞蚊的初始情況。圖7是這9個(gè)媽蟻運(yùn)行到250牛時(shí)步時(shí)的情況，表現(xiàn)岀一種混沌狀態(tài)口闈S是這9個(gè)螞蟻運(yùn)彳」到1250個(gè)時(shí)步時(shí)的情況，從圖屮川以看到,左邊口經(jīng)有兩個(gè)螞蟻丿I：辟山『自l2的公路，并不斷遠(yuǎn)離初始位置中多個(gè)賂蟻村互作川，使開辟公路的時(shí)間明顯減少。圖＜5第0步『圖＜5第0步『9個(gè)螞蟻國(guó)8第1250歩，9個(gè)螞蟻Langton螞蟻規(guī)則的一般結(jié)杲是，在螞蟻的運(yùn)動(dòng)過程小，無論初始空間結(jié)構(gòu)怎樣（灰色或口色元胞的構(gòu)形），螞蟻訪問的是無邊界的空間區(qū)域“Buniinovitch和Thioberkoy發(fā)現(xiàn)的

2.43概率元胞自動(dòng)機(jī)上述元胞自動(dòng)機(jī)都屬于確定型元胞自動(dòng)機(jī)「元胞自動(dòng)機(jī)的規(guī)則即元胞狀態(tài)的轉(zhuǎn)化函數(shù)也是確定的，稱為確定型演化規(guī)則，局部映射也是單值映射，然而本小節(jié)所提出的概率元胞自動(dòng)機(jī)的局部演化規(guī)則為加入隨機(jī)性的規(guī)則，局部映射為多值映射，概率元胞自動(dòng)機(jī)也稱為隨機(jī)元胞自動(dòng)機(jī).概率元胞自動(dòng)機(jī)可以描寫粒子的隨機(jī)運(yùn)動(dòng).多應(yīng)用于隨機(jī)擴(kuò)散.多粒子運(yùn)動(dòng),分形生長(zhǎng)等粒子運(yùn)動(dòng)模型，如下例即為使用概率元胞自動(dòng)機(jī)建立的森林火災(zāi)模型.元胞空間采用二維正方形網(wǎng)格自動(dòng)機(jī)，元胞節(jié)點(diǎn)包括三種狀態(tài)匚正在生長(zhǎng)的正在燃燒的樹變成空格位:1)2)樹，正在燃燒的樹和空白狀態(tài).在最初時(shí)刻所有網(wǎng)格用樹和正在燃燒的樹或空白這三種狀態(tài)填充，在每個(gè)時(shí)間步，元胞按以下規(guī)則進(jìn)行狀態(tài)正在燃燒的樹變成空格位:1)2)如杲生長(zhǎng)中的樹格位的最近的鄰居中有一棵樹處于燃燒狀態(tài)，則將它變?yōu)楣逕隣顟B(tài);3） 在空白格位的節(jié)點(diǎn)，穩(wěn)木以概率戸生長(zhǎng)14） 考慮到閃電的作用，在最近的鄰居中沒有正在懶燒的樹的情況下，樹在每個(gè)時(shí)間步以概率f變成燃燒中的樹模型演化情況如圖2.7所示，取參數(shù)曲風(fēng)發(fā)現(xiàn)當(dāng)上t￡）時(shí)，此模塑具有自組織的臨界狀態(tài)*這也說明」若聲維持穩(wěn)定狀態(tài)，表征系統(tǒng)的幾個(gè)物理量將具有乘方律行為.圖2圖2J基于槪率元鮑自動(dòng)機(jī)的森林火衣模型砂=1/500)rig.2.7Forestfiremodelbasedonprobabilitj1cellularautomaumode^在現(xiàn)實(shí)生活中，許許多多尙然現(xiàn)象都和概率有關(guān)F諸如擴(kuò)散、傳播、流行等等很多物理現(xiàn)象都可以使用概率元胞自動(dòng)機(jī)進(jìn)行模擬,概率元胞自動(dòng)機(jī)也搭比較重要的元腕自動(dòng)機(jī)模型之一”245格子氣自動(dòng)機(jī)格子氣自動(dòng)機(jī)（LatticeGasAutomata,LGA又稱格氣機(jī)）是元胞自動(dòng)機(jī)在流體力學(xué)與統(tǒng)計(jì)物瑾中的具體化*也是元胞自動(dòng)機(jī)在科學(xué)研究領(lǐng)域成功應(yīng)用的范例〔叫相對(duì)于竹生命游戲“來說?格子氣自動(dòng)機(jī)更注重于模型加實(shí)問性”它利用元胞自動(dòng)機(jī)的動(dòng)態(tài)特征，來模擬流體粒子的運(yùn)動(dòng)」("由于流體粒予不會(huì)輕易從模型空間中消失，這個(gè)特征需要格子氣自動(dòng)機(jī)是一個(gè)可逆元胞自動(dòng)機(jī)模型.(2)格了氣日動(dòng)機(jī)的鄰居模梨逓常采目Margulos棗型，即它的迎妁呈基于一個(gè)2X2的網(wǎng)格空間的一它的規(guī)則如圖入9所示：0■■Q￡2.9Mai^ulos類塑艷子氣住丈機(jī)鄰居模昭

Fig.29NeighborsnodelofLatticeGasAutomata區(qū)里老色球代董涼體粒子，白色球代曩空初元胞?町以看出,格于氣自動(dòng)機(jī)不像其它的元胞自動(dòng)機(jī)模型那拝，只以一個(gè)元胞(常被稱為中心元胞)為研究對(duì)象*考慮幷狀態(tài)的轉(zhuǎn)換’而是舟慮包含四個(gè)兀胞的一*側(cè)方塊.(3)依陽上述規(guī)則麻鄰馬棋里在計(jì)算完-次后?需耍濤這個(gè)2x2的模板沿跡角方向滑動(dòng).再計(jì)算一次.那么*一個(gè)流體粒子的運(yùn)動(dòng)需要兩個(gè)時(shí)間步才能完成.從時(shí)何和空間的角度看，格F氣口動(dòng)機(jī)扣對(duì)其地的尸胞門動(dòng)機(jī)模冬軾有較為獨(dú)轉(zhuǎn)的持社-格P自動(dòng)機(jī)作為…種特殊類型的元胞H動(dòng)機(jī)已成為流休動(dòng).力學(xué)中的一個(gè)重要領(lǐng)域，幾乎獨(dú)立于元胞自動(dòng)機(jī)研究之外了.元胞自動(dòng)機(jī)模擬單晶長(zhǎng)大USEMSFLIBPARAMETERIR=400,JR=400INTEGERIS(0:IR+1,0:JR+1),TMAX,ISN(1:8),NSTATE,T,NR,IX0,IY0,IX,JY!!根據(jù)過去狀態(tài)IS,定義現(xiàn)在狀態(tài)為IS1；為了突出邊界，特別定義ISN1INTEGERIS1(0:IR+1,0:JR+1),ISN1(1:8)WRITE(*,*)"PLEASEINPUTTHETIMESTEP"READ(*,*)TMAXISEED=RTC()IRC=IR/2JRC=JR/2!!!定義基體體積能為10,晶粒體積能為1IS=10IS(IRC,JRC)=1!!在循環(huán)前定義現(xiàn)在狀態(tài)數(shù)組IS1的初始值IS1=ISOPEN(1,FILE="E:\LUKE.DAT")DOT=1,TMAX!!每次循環(huán)前重新定義過去狀態(tài)數(shù)組ISIS=IS1!邊界條件IS(0,1:JMAX)=IS(IMAX,1:JMAX)IS(IMAX+1,1:JMAX)=IS(1,1:JMAX)IS(0:IMAX+1,0)=IS(0:IMAX+1,JMAX)IS(0:IMAX+1,JMAX+1)=IS(0:IMAX+1,1)!! 整個(gè)基體上，各個(gè)點(diǎn)上的狀態(tài)都要根據(jù)規(guī)則改變，而非隨機(jī)取點(diǎn)改變DOIX=1,IRDOJY=1,JRISN=(/IS(IX-1,JY-1),IS(IX-1,JY),IS(IX-1,JY+1),IS(IX,JY-1)& ,IS(IX,JY+1),IS(IX+1,JY-1),IS(IX+1,JY),IS(IX+1,JY+1)/)E0=COUNT(ISN.NE.IS(IX,JY))!如果不是晶粒邊界，則跳出，重新循環(huán)IF(E0.EQ.0)CYCLE!隨機(jī)尋找一個(gè)相鄰點(diǎn)NR=8*RAN(ISEED)+1NSTATE=ISN(NR)!判斷與相鄰點(diǎn)的能量差，并決定是否改變狀態(tài)E=COUNT(ISN.NE.NSTATE)RD=RAN(ISEED)IG=NSTATE-IS(IX,JY)DE=E-E0+IG+2.5*RD-1.25!! 用現(xiàn)在狀態(tài)數(shù)組IS1記錄邊界狀態(tài)的改變IF(DE.LT.0.0)IS1(IX,JY)=NSTATEENDDOENDDO!! 每循環(huán)20次在顯示屏幕上刷新狀態(tài)(顏色)DOIX=1,IRDOJY=1,JR!IF(MOD(T,20).EQ.0)THENISN1=(/IS1(IX-1,JY-1),IS1(IX-1,JY),IS1(IX-1,JY+1),IS1(IX,JY-1)&,IS1(IX,JY+1),IS1(IX+1,JY-1),IS1(IX+1,JY),IS1(IX+1,JY+1)/)ISRE=SETCOLOR(IS(IX,JY))! 如果是邊界，定義特別顏色15,加以區(qū)分IF(COUNT(ISN1.NE.IS1(IX,JY)).GT.0)ISRE=SETCOLOR(15)ISRE=SETPIXEL(IX,JY)!ENDIFENDDOENDDO!!!記錄循環(huán)次數(shù)和對(duì)應(yīng)的晶粒面積WRITE(1,*)T,SQRT(1.0*COUNT(IS.EQ.1))ENDDOCLOSE(1)END元胞自動(dòng)機(jī)模擬單晶長(zhǎng)大，再附加一個(gè)規(guī)則USEMSFLIBPARAMETERIR=400,JR=400,NMAX=100INTEGERIS(0:IR+1,0:JR+1),TMAX,ISN(1:8),NSTATE,T,NR,IX0,IY0,IX,JYINTEGERIGV(0:NMAX)!!根據(jù)過去狀態(tài)IS,定義現(xiàn)在狀態(tài)為IS1；為了突出邊界，特別定義ISN1INTEGERIS1(0:IR+1,0:JR+1),ISN1(1:8)WRITE(*,*)"PLEASEINPUTTHETIMESTEP"READ(*,*)TMAXISEED=RTC()IRC=IR/2JRC=JR/2! 定義基體體積能為10,晶粒體積能為1IS=10IS(IRC,JRC)=1!!在循環(huán)前定義現(xiàn)在狀態(tài)數(shù)組IS1的初始值IS1=IS!!OPEN(1,FILE="E:\LUKE.DAT")DOT=1,TMAX!!每次循環(huán)前重新定義過去狀態(tài)數(shù)組ISIS=IS1!!!邊界條件IS(0,1:JR)=IS(IR,1:JR)IS(IR,1:JR)=IS(iR+1,1:JR)IS(0:IR+1,0)=IS(0:IR+1,JR)IS(0:IR+1,JR+1)=IS(0:IR+1,1)!! 整個(gè)基體上，各個(gè)點(diǎn)上的狀態(tài)都要根據(jù)規(guī)則改變，而非隨機(jī)取點(diǎn)改變DOIX=1,IRDOJY=1,JR!!ISN=(/IS(IX-1,JY-1),IS(IX-1,JY),IS(IX-1,JY+1),IS(IX,JY-1)& ,IS(IX,JY+1),IS(IX+1,JY-1),IS(IX+1,JY),IS(IX+1,JY+1)/)E0=COUNT(ISN.NE.IS(IX,JY))!如果不是晶粒邊界，則跳出，重新循環(huán)IF(E0.EQ.0)CYCLE!隨機(jī)尋找一個(gè)相鄰點(diǎn)sss=0doiii=1,8NSTATE=ISN(iii)!判斷與相鄰點(diǎn)的能量差，并決定是否改變狀態(tài)E=COUNT(ISN.NE.NSTATE)RD=RAN(ISEED)IG=NSTATE-IS(IX,JY)DE=E-E0+IG+2.5*RD-1.25!! 用現(xiàn)在狀態(tài)數(shù)組IS1記錄邊界狀態(tài)的改變IF(DE.LT.0.0)thensss=sss+1endifenddoif(sss.gt.5)IS1(IX,JY)=11-is(ix,jy)ENDDOenddo!! 每循環(huán)20次在顯示屏幕上刷新狀態(tài)(顏色)DOIX=1,IRDOJY=1,JRIF(MOD(T,20).EQ.0)THENISN1=(/IS1(IX-1,JY-1),IS1(IX-1,JY),IS1(IX-1,JY+1),IS1(IX,JY-1)&,IS1(IX,JY+1),IS1(IX+1,JY-1),IS1(IX+1,JY),IS1(IX+1,JY+1)/)ISRE=SETCOLOR(IS(IX,JY))! 如果是邊界，定義特別顏色15,加以區(qū)分IF(COUNT(ISN1.NE.IS1(IX,JY)).GT.0)ISRE=SETCOLOR(15)ISRE=SETPIXEL(IX,JY)ENDIFENDDOENDDO!!!記錄循環(huán)次數(shù)和對(duì)應(yīng)的晶粒面積WRITE(1,*)T,SQRT(1.0*COUNT(IS.EQ.1))ENDDOCLOSE(1)END3.2基于CA的品粒長(zhǎng)大模擬原理和方法3.2.1元胞選擇及其狀態(tài)確定本文模擬采用四邊形元胞，所研究空間為二錐四邊形網(wǎng)格平面，元胞鄰居為圖3-1所示的Moore型鄰居.a:dELU圖3-1Mooit樂鄰居定文毎一牛元胞在甘1時(shí)刻的狀態(tài)白其相鄰元胞在t時(shí)刻的狀態(tài)所決定，即：旳+1)=FS⑴上⑴，附d⑴岸⑴丿◎薊(啪 <3-1)式中>鞏什1)為元胞e在1+1時(shí)刻的狀態(tài)，口⑴?附)分別為元胞e及其相鄒兀?胞在『時(shí)刻的狀態(tài)”F対元砲狀態(tài)轉(zhuǎn)變規(guī)貝J,322元胞狀態(tài)的轉(zhuǎn)換規(guī)則在本文每一模擬時(shí)間步C1CAS)內(nèi)，網(wǎng)格內(nèi)每個(gè)元胞秋態(tài)都按臥下規(guī)皿判斷其足否轉(zhuǎn)換；(1)如元胞憑圍的八個(gè)元胞■b*s乩仁旨h、i的狀態(tài)均與其相同「則茫下一個(gè)CAS時(shí).元胞亡的狀態(tài)囉持不變?⑵@)如兀胞匕d.i；h中任意三個(gè)同為R狀態(tài)，貝J注下-個(gè)CAS時(shí)，元胞亡的狀態(tài)轉(zhuǎn)變?yōu)?(b)ta.元胞血G臥i中任意三個(gè)同為/狀態(tài)，貝恠下一個(gè)CAS時(shí)，元胞匕的狀態(tài)轉(zhuǎn)變?yōu)閮孩羌僭O(shè)晶界能均質(zhì)分右，且一&元胞必須克服能量壁空以達(dá)錨其新狀態(tài)?USEMSFLIBPARAMETERIR=400,JR=400,NMAX=100INTEGERIS(0:IR+1,0:JR+1),TMAX,ISN(1:8),NSTATE,T,NR,IX0,IY0,IX,JYINTEGERIGV(0:NMAX)!!根據(jù)過去狀態(tài)IS,定義現(xiàn)在狀態(tài)為IS1；為了突出邊界，特別定義ISN1INTEGERIS1(0:IR+1,0:JR+1),ISN1(1:8)WRITE(*,*)"PLEASEINPUTTHETIMESTEP"READ(*,*)TMAXISEED=RTC()!定義基體體積能為10,晶粒體積能為1IGV=1IGV(0)=10!定義晶粒長(zhǎng)大位置，并附能量DOI=1,NMAXIX0=IR*RAN(ISEED)+1JY0=JR*RAN(ISEED)+1IF(IS(IX0,JY0).NE.0)CYCLEIS(IX0,JY0)=IENDDO!! 在循環(huán)前定義現(xiàn)在狀態(tài)數(shù)組IS1的初始值IS1=IS!!OPEN(1,FILE="E:\LUKE.DAT")DOT=1,TMAX!!每次循環(huán)前重新定義過去狀態(tài)數(shù)組ISIS=IS1!!!邊界條件IS(0,1:JMAX)=IS(IMAX,1:JMAX)IS(IMAX+1,1:JMAX)=IS(1,1:JMAX)IS(0:IMAX+1,0)=IS(0:IMAX+1,JMAX)IS(0:IMAX+1,JMAX+1)=IS(0:IMAX+1,1)!尋找晶粒邊界，計(jì)算能量，改變狀態(tài)。!! 整個(gè)基體上，各個(gè)點(diǎn)上的狀態(tài)都要根據(jù)規(guī)則改變，而非隨機(jī)取點(diǎn)改變DOIX=1,IRDOJY=1,JR!!ISN=(/IS(IX-1,JY-1),IS(IX-1,JY),IS(IX-1,JY+1),IS(IX,JY-1)& ,IS(IX,JY+1),IS(IX+1,JY-1),IS(IX+1,JY),IS(IX+1,JY+1)/)E0=COUNT(ISN.NE.IS(IX,JY))!如果不是晶粒邊界，則跳出，重新循環(huán)IF(E0.EQ.0)CYCLE!隨機(jī)尋找一個(gè)相鄰點(diǎn)NR=8*RAN(ISEED)+1NSTATE=ISN(NR)!判斷與相鄰點(diǎn)的能量差，并決定是否改變狀態(tài)E=COUNT(ISN.NE.NSTATE)RD=RAN(ISEED)IG=IGV(NSTATE)-IGV(IS(IX,JY))DE=E-E0+IG+2.5*RD-1.25!! 用現(xiàn)在狀態(tài)數(shù)組IS1記錄邊界狀態(tài)的改變IF(DE.LT.0.0)IS1(IX,JY)=NSTATE!!!! ISRE=SETCOLOR(MOD(IS(IX,JY),15))!! ISRE=SETPIXEL(IX,JY)ENDDOENDDO!!每循環(huán)20次在顯示屏幕上刷新狀態(tài)(顏色)DOIX=1,IRDOJY=1,JRIF(MOD(T,20).EQ.0)THENISN1=(/IS1(IX-1,JY-1),IS1(IX-1,JY),IS1(IX-1,JY+1),IS1(IX,JY-1)& ,IS1(IX,JY+1),IS1(IX+1,JY-1),IS1(IX+1,JY),IS1(IX+1,JY+1)/)ISRE=SETCOLOR(MOD(IS(IX,JY),15))! 如果是邊界，定義特別顏色15,加以區(qū)分IF(COUNT(ISN1.NE.IS1(IX,JY)).GT.0)ISRE=SETCOLOR(15)ISRE=SETPIXEL(IX,JY)ENDIFENDDOENDDO!!!記錄循環(huán)次數(shù)和對(duì)應(yīng)的晶粒面積WRITE(1,*)T,SQRT(1.0*COUNT(IS1.gt.0))ENDDOCLOSE(1)END錯(cuò)誤的多晶長(zhǎng)大USEMSFLIBPARAMETERIR=400,JR=400,NMAX=100INTEGERIS(0:IR+1,0:JR+1),TMAX,ISN(1:8),NSTATE,T,NR,IX0,IY0,IX,JYINTEGERIGV(0:NMAX)!! 根據(jù)過去狀態(tài)IS,定義現(xiàn)在狀態(tài)為IS1；為了突出邊界，特別定義ISN1INTEGERIS1(0:IR+1,0:JR+1),ISN1(1:8)WRITE(*,*)"PLEASEINPUTTHETIMESTEP"READ(*,*)TMAXISEED=RTC()!!!定義基體體積能為10,晶粒體積能為1IGV=1IGV(0)=10!定義晶粒長(zhǎng)大位置，并附能量DOI=1,NMAXIX0=IR*RAN(ISEED)+1JY0=JR*RAN(ISEED)+1IS(IX0,JY0)=IENDDO!! 在循環(huán)前定義現(xiàn)在狀態(tài)數(shù)組IS1的初始值IS1=IS!!OPEN(1,FILE="E:\LUKE.DAT")DOT=1,TMAX!!每次循環(huán)前重新定義過去狀態(tài)數(shù)組ISIS=IS1!!!邊界條件IS(0,1:JMAX)=IS(IMAX,1:JMAX)IS(IMAX+1,1:JMAX)=IS(1,1:JMAX)IS(0:IMAX+1,0)=IS(0:IMAX+1,JMAX)IS(0:IMAX+1,JMAX+1)=IS(0:IMAX+1,1)!! 整個(gè)基體上，各個(gè)點(diǎn)上的狀態(tài)都要根據(jù)規(guī)則改變，而非隨機(jī)取點(diǎn)改變DOIX=1,IRDOJY=1,JRISN=(/IS(IX-1,JY-1),IS(IX-1,JY),IS(IX-1,JY+1),IS(IX,JY-1)& ,IS(IX,JY+1),IS(IX+1,JY-1),IS(IX+1,JY),IS(IX+1,JY+1)/)E0=COUNT(ISN.NE.IS(IX,JY))!如果不是晶粒邊界，則跳出，重新循環(huán)IF(E0.EQ.0)CYCLE!隨機(jī)尋找一個(gè)相鄰點(diǎn)sss=0doiii=1,8! NR=8*RAN(ISEED)+1NSTATE=ISN(iii)!判斷與相鄰點(diǎn)的能量差，并決定是否改變狀態(tài)E=COUNT(ISN.NE.NSTATE)RD=RAN(ISEED)IG=IGV(NSTATE)-IGV(IS(IX,JY))DE=E-E0+IG+2.5*RD-1.25!! 用現(xiàn)在狀態(tài)數(shù)組IS1記錄邊界狀態(tài)的改變! IF(DE.LT.0.0)IS1(IX,JY)=NSTATEIF(DE.LT.0.0)thensss=sss+1nstate1=nstateif(sss.gt.5)IS1(IX,JY)=nstate1endifenddoENDDOENDDO!!每循環(huán)20次在顯示屏幕上刷新狀態(tài)(顏色)DOIX=1,IRDOJY=1,JRIF(MOD(T,40).EQ.0)THENISN1=(/IS1(IX-1,JY-1),IS1(IX-1,JY),IS1(IX-1,JY+1),IS1(IX,JY-1)& ,IS1(IX,JY+1),IS1(IX+1,JY-1),IS1(IX+1,JY),IS1(IX+1,JY+1)/)ISRE=SETCOLOR(MOD(IS(IX,JY),15))! 如果是邊界，定義特別顏色15,加以區(qū)分IF(COUNT(ISN1.NE.IS1(IX,JY)).GT.0)ISRE=SETCOLOR(15)ISRE=SETPIXEL(IX,JY)ENDIFENDDOENDDO!!記錄循環(huán)次數(shù)和對(duì)應(yīng)的晶粒面積WRITE(1,*)T,COUNT(IS1.GT.0)ENDDOCLOSE(1)ENDUSEMSFLIBPARAMETERIR=400,JR=400,NMAX=100INTEGERIS(0:IR+1,0:JR+1),TMAX,ISN(1:8),NSTATE,T,NR,IX0,IY0,IX,JYINTEGERIGV(0:NMAX)!!根據(jù)過去狀態(tài)IS,定義現(xiàn)在狀態(tài)為IS1；為了突出邊界，特別定義ISN1INTEGERIS1(0:IR+1,0:JR+1),ISN1(1:8),ISN2(1:8)WRITE(*,*)"PLEASEINPUTTHETIMESTEP"READ(*,*)TMAXISEED=RTC()!!! 定義基體體積能為10,晶粒體積能為1IGV=1IGV(0)=10!定義晶粒長(zhǎng)大位置，并附能量DOI=1,NMAXIX0=IR*RAN(ISEED)+1JY0=JR*RAN(ISEED)+1IS(IX0,JY0)=IENDDO!!在循環(huán)前定義現(xiàn)在狀態(tài)數(shù)組IS1的初始值IS1=IS!!OPEN(1,FILE="E:\LUKE.DAT")DOT=1,TMAX!!每次循環(huán)前重新定義過去狀態(tài)數(shù)組ISIS=IS1!!!邊界條件IS(0,1:JMAX)=IS(IMAX,1:JMAX)IS(IMAX+1,1:JMAX)=IS(1,1:JMAX)IS(0:IMAX+1,0)=IS(0:IMAX+1,JMAX)IS(0:IMAX+1,JMAX+1)=IS(0:IMAX+1,1)!! 整個(gè)基體上，各個(gè)點(diǎn)上的狀態(tài)都要根據(jù)規(guī)則改變，而非隨機(jī)取點(diǎn)改變DOIX=1,IRDOJY=1,JRISN=(/IS(IX-1,JY-1),IS(IX-1,JY),IS(IX-1,JY+1),IS(IX,JY-1)& ,IS(IX,JY+1),IS(IX+1,JY-1),IS(IX+1,JY),IS(IX+1,JY+1)/)E0=COUNT(ISN.NE.IS(IX,JY))!如果不是晶粒邊界，則跳出，重新循環(huán)IF(E0.EQ.0)CYCLE!隨機(jī)尋找一個(gè)相鄰點(diǎn)iiii=0doiii=1,8! NR=8*RAN(ISEED)+1NSTATE=ISN(iii)!判斷與相鄰點(diǎn)的能量差，并決定是否改變狀態(tài)E=COUNT(ISN.NE.NSTATE)RD=RAN(ISEED)IG=IGV(NSTATE)-IGV(IS(IX,JY))DE=E-E0+IG+2.5*RD-1.25!! 用現(xiàn)在狀態(tài)數(shù)組IS1記錄邊界狀態(tài)的改變! IF(DE.LT.0.0)IS1(IX,JY)=NSTATEIF(DE.LT.0.0)theniiii=iiii+1isn2(iiii)=nstateid1=iiii*ran(iseed)+1if(iiii.gt.5)IS1(IX,JY)=isn2(id1)elseendifenddoENDDOENDDO!!每循環(huán)20次在顯示屏幕上刷新狀態(tài)(顏色)DOIX=1,IRDOJY=1,JRIF(MOD(T,40).EQ.0)THENISN1=(/IS1(IX-1,JY-1),IS1(IX-1,JY),IS1(IX-1,JY+1),IS1(IX,JY-1)& ,IS1(IX,JY+1),IS1(IX+1,JY-1),IS1(IX+1,JY),IS1(IX+1,JY+1)/)ISRE=SETCOLOR(MOD(IS1(IX,JY),15))! 如果是邊界，定義特別顏色15,加以區(qū)分IF(COUNT(ISN1.NE.IS1(IX,JY)).GT.0)ISRE=SETCOLOR(15)ISRE=SETPIXEL(IX,JY)ENDIFENDDOENDDO!!!記錄循環(huán)次數(shù)和對(duì)應(yīng)的晶粒面積WRITE(1,*)T,COUNT(IS1.GT.0)ENDDOCLOSE(1)ENDA23局部規(guī)則的確定局部規(guī)則在建立元胞自動(dòng)機(jī)模型中是最遁要的步驟，溫度是如何傳遞的是根據(jù)這個(gè)規(guī)則而制定的，本文采用的是Wn.Neumann型的鄰居，所以研究中心元胞周圍的四個(gè)元胞對(duì)它溫度的影響?用熱力學(xué)來研究元胞之間溫度的傳遞是局部規(guī)則制定的基礎(chǔ).物體內(nèi)溫度的變化與熱量有關(guān).熱量的增減引起溫度的升降，為了確定溫度場(chǎng)，必須研究熱量在物體內(nèi)傳導(dǎo)的過程。溫度在時(shí)間域和空間域中的分布，稱為溫度場(chǎng).它可表示為公式(4-1》T=T(x,齊為t) (4J)若溫度不隨時(shí)間變化，即#=0,則T-T(k,y,z)稱為穩(wěn)定溫度場(chǎng).若溫度沿z向不變，則稱為平面溫度場(chǎng)。在任一瞬時(shí).連接場(chǎng)內(nèi)相同溫度值的各點(diǎn)，就得到此時(shí)刻的等溫面。沿等溫面切向,溫度不變口而沿垂直等溫面的法向，溫度變化率最大。表示一點(diǎn)最大埴溫率的矢量，稱為溫度梯度.即公式(42)：“st.arjrdTt.dr(4-2)V"%需知云r齊+無瓦其中血是單位矢量，沿等溫面的法線指向增溫方向。(4-2)在單位時(shí)問內(nèi)通過單位面積的熱JL稱為熱漬密度，即公式(4-3)；IdQ(切dH——旦

sdt(切菓■點(diǎn)的最大熱流密度矢量是沿等溫面的法線且指向降溫方向「根韜傳導(dǎo)定律，熱流密度與溫度梯度成正比而方問相反，即公式(M)：(4-4)(4-5)率：忑耍為x方向上的溫度梯度.我們可以得到導(dǎo)熱方程為公式(46)(4-4)(4-5)率：忑耍為x方向上的溫度梯度.我們可以得到導(dǎo)熱方程為公式(46)：OX式中p為材料的密度(kg/M)：c為材料的比熱容(J/(kgK))；q=—AVT=-%久2~dn幾稱為導(dǎo)熱系數(shù)。熱流密度在x方向的分量是公式件刃；q$=-Z—cosC^x)=-A—^x.y.z)dn dx類似，在任意方向的熱流密度，等于導(dǎo)熱系數(shù)乘以溫度在該方向的變化率的負(fù)值&其中，亦為x方向上的熱流密度；九為材料沿K,y,z方向的熱導(dǎo)r為時(shí)間G)；人，九人分別為材料沿x,y,z方向的熱導(dǎo)率(W/(mK))；Q二0(x,Z玷)為材料內(nèi)部的熱源密度(W/kg)。上式中，第一項(xiàng)為體元升溫需要的熱量：右側(cè)第一、二和三項(xiàng)是由x,y和2方向流入體元的熱量：躍后一項(xiàng)體元內(nèi)熱源產(chǎn)生的熱量。微分方程的物理意義，體元升溫所需的熱量應(yīng)該等于流入體元的熱量與體元內(nèi)產(chǎn)生的熱量的總和。由于我們研究的是二維的溫度場(chǎng)模型，所以簡(jiǎn)化后的方程，二維非穩(wěn)態(tài)熱傳導(dǎo)方程為公式(4?〃(4-7)(4-7)在區(qū)域內(nèi)進(jìn)行網(wǎng)格劃分，將連續(xù)的求解域離散為不連續(xù)的點(diǎn)，形成離散網(wǎng)格,網(wǎng)格步長(zhǎng)分別為兀何一九=Ar,格,網(wǎng)格步長(zhǎng)分別為兀何一九=Ar,畑_九=切劃分的單元格步長(zhǎng)可以是均勻的，也可以是不均勻的。此模型我們?nèi)卧癫介L(zhǎng)是均勻的，如3所表示:對(duì)于二維各向同性，無內(nèi)熱源的非穩(wěn)態(tài)熱傳導(dǎo)微分方程，利用中心差分法來研究，差分解法是將函數(shù)用離散點(diǎn)的函數(shù)值表示，將微分用有限差分表示,將導(dǎo)數(shù)用有限差商表示，從而將微分方程化為差分方程一代數(shù)方程，使求解微分方程的問題化為求解代數(shù)方程的問題.中心差分法公式件8)、(4-9).(4-10)如下t圖4-3鄰居模型Fig.4-3Neighborttiodel白m=遲tj一2氏寸加 (Axpdy2 M竺迫m&t di將(4將人(4-9).(4-10)代入(4?7)可得公式(4-11)：叫嘰_arJf(4-8)(4-9)(4-10)戸中碼腫-2山｛j+dW(Ax)2 如,在劃分單元時(shí)，我們令則上式可化為公式(4-12).A/4+4*4/供<Ar)2J(4-11)(4-12)上式即為元胞在絕熱條件下的溫度傳遞公式?即可以通過中心元胞的溫度値和鄰居元胞的溫度值得到下一時(shí)刻該元胞的溫度“我們把公式(?12)定義為元胞自動(dòng)機(jī)的局部規(guī)則口通過對(duì)元胞自動(dòng)機(jī)溫度場(chǎng)模型中元胞，元胞空間，鄰居規(guī)則「局部規(guī)則的建立憲成，就刀片溫度場(chǎng)模型建立完咸，我們把離散的網(wǎng)格定文為元胞，元胞空間為二維p鄰居來用VbmNeumann型，局部規(guī)則為推導(dǎo)公式件1幼!相變byUsingBoundaryTrackingMethord,2008,11.23!graingrowthvelocityis"1.0*(IT-NT(IST))"or