计算材料学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:第几天实验,I:第几步跳跃! XXM(I):THEMEANOFXXWRITE(*,*)"实验天数JMAX，实验次数IMAX"READ(*,*)JMAX,IMAXISEED=RTC()DOJ=1,JMAX!第几天实验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:第几天实验,I:第几步跳跃XYM(I):THEMEANOFXYWRITE(*,*)"实验天数JMAX，实验次数IMAX"READ(*,*)JMAX,IMAXISEED=RTC()DOJ=1,JMAX!第几天实验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:第几天实验,I:第几步跳跃!XYM(I):THEMEANOFXYWRITE(*,*)"实验天数JMAX，实验次数IMAX"READ(*,*)JMAX,IMAXXN=(/0,0,-1,1/)YN=(/-1,1,0,0/)ISEED=RTC()DOJ=1,JMAX!第几天实验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做三维空间随机行走？？留作业! 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:第几天实验,I:第几步跳跃!XYM(I):THEMEANOFXYWRITE(*,*)"实验天数JMAX，实验次数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!第几天实验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!在的区域中画一个圆PARAMETERXR=400,YR=400INTEGERR,S(1:XR,1:YR)X0=XR/2!圆心位置X0，YOY0=YR/2R=MIN(X0-10,Y0-10)!圆半径S=0!像素的初始状态(颜色)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画一个圆(1、如何选出晶界区域；2、进一步加深对画图的理解)USEMSFLIBINTEGERXR，YR!在的区域中画一个圆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!像素的初始状态(颜色)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模拟一个晶粒的缩小(1、 如何定义基体和晶粒；(2、 如何寻找边界；(3、 如何计算能量(构造或描述问题的概率过程)步骤：1、在基体上画一个原始晶粒，并赋予基体和原始晶粒不同的状态(取向号或体积能)2、寻找晶界3、 MC的一个时间步(晶粒长大过程中一种随机性)4、 计算晶界上网格点的相互作用，每个格点的相互作用导致晶粒长大随机选取一个初始格点；若此点属于晶界，那么可以随机转变为它邻居的状态，若不是晶界，则跳出循环，不发生转变；计算转变前后的能量变化，/E=/E+/E+/Evsq(自由能=体积能+表面能+能量起伏)若/E小于0则新晶相被接受，网格状态发生转变。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!定义基体和圆晶粒分别为状态1、状态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")!寻找晶粒边界，计算能量，改变状态。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))!如果不是圆晶粒边界，则跳出，重新循环IF(E0.EQ.0)CYCLE!随机寻找一个相邻点NR=8*RAN(ISEED)+1NSTATE=ISN(NR)!判断与相邻点的能量差，并决定是否改变状态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、 演示体积能互换的转变情况；2、演示IX=IR*RAN(ISEED)+1JY=JR*RAN(ISEED)+1的替换；3、能量关系变化；)作业3：基体上单晶粒形核长大(1、 选取初始形核点(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!定义基体和圆晶粒分别为状态10、状态2IS=10IS(IRC,JRC)=2OPEN(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))!如果不是圆晶粒边界，则跳出，重新循环IF(E0.EQ.0)CYCLE!随机寻找一个相邻点NR=8*RAN(ISEED)+1NSTATE=ISN(NR)!判断与相邻点的能量差，并决定是否改变状态基体上多晶粒形核长大，（1基体上多晶粒形核长大，（1、 如何生成多晶粒形核核心（2、 如何定义、并区别多晶粒3、为什么要有边界条件）1步骤：1、定义变量：体积能变量（一维数组）；基体各像素点变量（二维数组），邻居的取向号（一维数组）等2、读取晶粒生长步长；3、寻找晶界4、计算晶界上网格点的相互作用，每个格点的相互作用导致晶粒长大或缩小：①随机选取一个初始格点；②若此点属于晶界，那么可以随机转变为它邻居的状态，若不是晶界，则跳出循环，不发生转变；③计算转变前后的能量变化，/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则新晶相被接受，网格状态发生转变。USEMSFLIB

定义常数名PARAMETERIR=400,JR=400,NMAX=100定义变量：体积能变量(一维数组)；基体各像素点变量(二维数组)，邻居的取向号(一维数组)等INTEGERIS(0:IR+1,0:JR+1),TMAX,ISN(1:8),NSTATE,T,NR,IX0,IY0,IX,IYINTEGERIGV(0:NMAX)读取晶粒生长步长；WRITE(*,*)"PLEASEINPUTTHETIMESTEP"READ(*,*)TMAXISEED=RTC()定义基体体积能为10,所有晶粒体积能为1IGV=1IGV(0)=10定义晶粒长大位置(100个形核点初始形核位置)，并附已不同的取向号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))!如果不是晶粒边界，则跳出，重新循环IF(E0.EQ.0)CYCLE!随机寻找一个相邻点NR=8*RAN(ISEED)+1NSTATE=ISN(NR)!判断与相邻点的能量差，并决定是否改变状态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晶粒长大与时间的关系2形核数目与的关系3、与EBSD比较4、为什么圆弧会长成近似的直线Monte-Carlo方法模拟晶粒长大过程中的新研究进展第一个问题：采用Monte-Carlo方法存在以下几方面的问题（1）模拟过程只计及最近邻格点的相互作用，难于考虑畸变场等物理场以及晶界形态对晶粒长大过程的影响，不利于在唯象层次上深入开展对晶粒长大过程的研究。应当指出，大部分的计算机模拟主要集中于二维系统，虽然能够在一定程度上提供晶粒长大的有关动力学和拓扑学的信息，但是，任何多晶体材料中的晶粒都是三维的，所有二维晶粒长大的模拟都是基于某种特定的假设或简化，因此模拟三维晶粒的长大更具有实际意义。同时，在目前的研究中，对晶粒生长过程的形貌演化模拟的逼真度比较好，但在生长模型的边界处理等方面还有待完善，特别是晶粒生长与晶化温度、时间、气氛等参数密切相关，这种复杂工艺条件下的晶粒生长过程的模拟是当前该领域正待解决的难题。另外，国内外用各种方法模拟晶粒长大，主要集中在对方法本身的讨论、模拟晶粒分布函数特点及晶粒形貌等方面，而将模拟方法与实际工艺联系起来的研究者比较少。第二个问题：算法种类：1、标准算法2、改进型算法3、其它算法第三个问题：模拟晶粒长大的一些应用1、 金属凝固过程的形核；2、 再结晶及存在织构的再结晶（形核位置时考虑能量关系）；3、 第二相粒子的影响（第二相粒子对晶粒长大有阻碍作用.这种情况下模拟较复杂.由于无法用解析解定量描述,用MC方法模拟显得尤为重要；(a)j--H-ql*6000$j-i-IWOMCSgf—2.5l1ODD0Mc$ (d)/-S-ft—i2&afl-MC$SI<5不同晖粗时创址二粕粒子时兄粧就观胡级的匪吨俨I5、异常晶粒的长大；6、三维晶粒的长大(模拟不同的点阵结构生长过程)；山东大学中国有色金属学报，2008，4252()105MCS圈斗不同楝型榄拟得到闾晶粒辰人指数Flg+252()105MCS圈斗不同楝型榄拟得到闾晶粒辰人指数Flg+4Grfliugi'owthexponentfb]Ldififei'enrniodels・一W)3I、0.4舫±6024°—Model2.0391=0.010*—Model\0.501±0.011200 400 600 800 10007、蒙特卡洛方法不是数学的方法，而是一种多次模拟实验的方法。元胞自动机程序(生命永不停止)usemsflibPARAMETERIR=400,JR=400,NMAX=40000!NMAX随机产生的生命种子integeris(0:1001,0:1001),is1(0:1001,0:1001),ISN(1:8),Tmax,num!is-基体的二维数组print*,'pleaseinputloop(Tmax)'read*,TmaxISEED=RTC()is=15！“死”的状态，基体为白色！赋予生命的种子，“活”的状态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！判断邻居状态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元胞自动机—产生和发展.四个阶段:1940s诞生：VonNeumann自我复制机.1960-70S起步：JH・Conway生命游戏.1980s理论研究：S.WolframCA分类.1980-90s应用：HPP-FHP格子气自动机・C・LangtonN・Packard人工生命元胞自动机（CellularAutomata,简称CA,也有人称其为细胞自动机、点格自动机、分子自动机或单元自动机）是一种建立在离散的时间和空间上的动力学系统。散布在规则格网（LatticeGrid）中的每一元胞（Cell）取有限的离散状态，遵循同样的作用规则，依据确定的局部规则作同步更新。大量元胞通过简单的相互作用而构成动态系统的演化。与一般的动力学模型不同，元胞自动机不是由严格定义的物理方程或函数确定，而是由一系列模型构造的规则构成，凡是满足这些规则的模型都可以算做是元胞自动机模型。因此，元胞自动机是一类模型的总称，或者是一个方法框架。其特点是时间、空间、状态都是离散的，每个变量只取有限多个状态，且其状态改变的规则在时间和空间上都是局部的。元胞自动机是一种对具有局域连通性的格点，应用局部（有时为中等范围）确定性或概论性的转换规则来描述在离散空间和时间上复杂系统演化规律的同步算法。元胞自动机可以采取任意的维数，自动机的格点描述被认为与所选模型密切相关的基本的系统实体。独立的格点可以表示连续的体积单元，原子颗粒大原子团、汽车、人口、晶格缺陷等。通过作用于每个格点的确定性或概率性转换规则，自动机的演化得以实施。这些规则决定了每个元胞的状态，并将其描述为其前一状态和邻近元胞的函数。在计算状态转变过程中所考虑的邻近位置的数目确定了相互作用的范围和自动机的局部演化。！转变概率元胞自动机特点：时间、空间、状态都离散的动力学模型，利用确定性或概论性的转换规则来描述在离散空间和时间上复杂系统演化2_4.2兰顿的蚂蚁唱蚁规则元胞自动机是由ChrisLangton利GregTurk发现的叫滋.自动札過过极简单的两条规则來模拟蚂蚁的行为规律，吗蚁在方形网格上远动，网格分为黑色和白色，演化规则対’D如果码蚁当前处于黑色网格」则向左转9胪.并将网格涂成H色.2｝细果嶋蚁刍前处于白仪网格，则向启转9仇并将网格涂成黒色.结昊发现*这只小妈蚁表现出了用当复杂的着为"从H色的元胞窣间开始’经过入绡500步厶胡“它就绕回碾位，之后又进入-种混沌状态’这种无规则运动的状态-亩持缜到大约1M00步后，蚂蚁窦然又出现非常堀则的运动形式，它会F辟一条笙自的路线井迁旦前逵一2.6中临位卜（b）分别表现了处于混沱状态和规则运动狀态的蚂蚁自动机.（b）图中蚂蚁开辟的笔直的路线被称为公路⑶】,与网格方向成4和角，公路在不遇到障碍物的情捏下会一直延伸下丧一閨3第391步,,一个鹏蚁图5笫10500归一个訥蚁图4第冗00步』一个蚂蚁3.2多个蚂蚁的情况这里対心网格屮加入名个蚂蚁的情况进彳亍了模拟。图6是件网楙屮加入9个蚂蚊的初始情况。图7是这9个妈蚁运行到250牛时步时的情况，表现岀一种混沌状态口闱S是这9个蚂蚁运彳」到1250个时步时的情况，从图屮川以看到,左边口经有两个蚂蚁丿I：辟山『自l2的公路，并不断远离初始位置中多个賂蚁村互作川，使开辟公路的时间明显减少。图＜5第0步『图＜5第0步『9个蚂蚁国8第1250歩，9个蚂蚁Langton蚂蚁规则的一般结杲是，在蚂蚁的运动过程小，无论初始空间结构怎样（灰色或口色元胞的构形），蚂蚁访问的是无边界的空间区域“Buniinovitch和Thioberkoy发现的

2.43概率元胞自动机上述元胞自动机都属于确定型元胞自动机「元胞自动机的规则即元胞狀态的转化函数也是确定的，称为确定型演化规则，局部映射也是单值映射，然而本小节所提出的概率元胞自动机的局部演化规则为加入随机性的规则，局部映射为多值映射，概率元胞自动机也称为随机元胞自动机.概率元胞自动机可以描写粒子的随机运动.多应用于随机扩散.多粒子运动,分形生长等粒子运动模型，如下例即为使用概率元胞自动机建立的森林火灾模型.元胞空间采用二维正方形网格自动机，元胞节点包括三种状态匚正在生长的正在燃烧的树变成空格位:1)2)树，正在燃烧的树和空白状态.在最初时刻所有网格用树和正在燃烧的树或空白这三种状态填充，在每个时间步，元胞按以下规则进行状态正在燃烧的树变成空格位:1)2)如杲生长中的树格位的最近的邻居中有一棵树处于燃烧状态，则将它变为罐烧状态;3） 在空白格位的节点，穩木以概率戸生长14） 考虑到闪电的作用，在最近的邻居中没有正在懒烧的树的情况下，树在每个时间步以概率f变成燃烧中的树模型演化情况如图2.7所示，取参数曲风发现当上t£）时，此模塑具有自组织的临界状态*这也说明」若聲维持稳定状态，表征系统的几个物理量将具有乘方律行为.图2图2J基于槪率元鮑自动机的森林火衣模型砂=1/500)rig.2.7Forestfiremodelbasedonprobabilitj1cellularautomaumode^在现实生活中，许许多多尙然现象都和概率有关F诸如扩散、传播、流行等等很多物理现象都可以使用概率元胞自动机进行模拟,概率元胞自动机也搭比较重要的元腕自动机模型之一”245格子气自动机格子气自动机（LatticeGasAutomata,LGA又称格气机）是元胞自动机在流体力学与统计物瑾中的具体化*也是元胞自动机在科学研究领域成功应用的范例〔叫相对于竹生命游戏“来说•格子气自动机更注重于模型加实问性”它利用元胞自动机的动态特征，来模拟流体粒子的运动」("由于流体粒予不会轻易从模型空间中消失，这个特征需要格子气自动机是一个可逆元胞自动机模型.(2)格了气日动机的邻居模梨逓常采目Margulos棗型，即它的迎妁呈基于一个2X2的网格空间的一它的规则如图入9所示：0■■Q£2.9Mai^ulos类塑艳子气住丈机邻居模昭

Fig.29NeighborsnodelofLatticeGasAutomata区里老色球代董涼体粒子，白色球代曩空初元胞•町以看出,格于气自动机不像其它的元胞自动机模型那拝，只以一个元胞(常被称为中心元胞)为研究对象*考虑幷状态的转换’而是舟虑包含四个兀胞的一*側方块.(3)依阳上述规则麻邻馬棋里在计算完-次后・需耍濤这个2x2的模板沿跡角方向滑动.再计算一次.那么*一个流体粒子的运动需要两个时间步才能完成.从时何和空间的角度看，格F气口动机扣对其地的尸胞门动机模冬轼有较为独转的持社-格P自动机作为…种特殊类型的元胞H动机已成为流休动.力学中的一个重要领域，几乎独立于元胞自动机研究之外了.元胞自动机模拟单晶长大USEMSFLIBPARAMETERIR=400,JR=400INTEGERIS(0:IR+1,0:JR+1),TMAX,ISN(1:8),NSTATE,T,NR,IX0,IY0,IX,JY!!根据过去状态IS,定义现在状态为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!!在循环前定义现在状态数组IS1的初始值IS1=ISOPEN(1,FILE="E:\LUKE.DAT")DOT=1,TMAX!!每次循环前重新定义过去状态数组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)!! 整个基体上，各个点上的状态都要根据规则改变，而非随机取点改变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))!如果不是晶粒边界，则跳出，重新循环IF(E0.EQ.0)CYCLE!随机寻找一个相邻点NR=8*RAN(ISEED)+1NSTATE=ISN(NR)!判断与相邻点的能量差，并决定是否改变状态E=COUNT(ISN.NE.NSTATE)RD=RAN(ISEED)IG=NSTATE-IS(IX,JY)DE=E-E0+IG+2.5*RD-1.25!! 用现在状态数组IS1记录边界状态的改变IF(DE.LT.0.0)IS1(IX,JY)=NSTATEENDDOENDDO!! 每循环20次在显示屏幕上刷新状态(颜色)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,加以区分IF(COUNT(ISN1.NE.IS1(IX,JY)).GT.0)ISRE=SETCOLOR(15)ISRE=SETPIXEL(IX,JY)!ENDIFENDDOENDDO!!!记录循环次数和对应的晶粒面积WRITE(1,*)T,SQRT(1.0*COUNT(IS.EQ.1))ENDDOCLOSE(1)END元胞自动机模拟单晶长大，再附加一个规则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)!!根据过去状态IS,定义现在状态为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!!在循环前定义现在状态数组IS1的初始值IS1=IS!!OPEN(1,FILE="E:\LUKE.DAT")DOT=1,TMAX!!每次循环前重新定义过去状态数组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)!! 整个基体上，各个点上的状态都要根据规则改变，而非随机取点改变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))!如果不是晶粒边界，则跳出，重新循环IF(E0.EQ.0)CYCLE!随机寻找一个相邻点sss=0doiii=1,8NSTATE=ISN(iii)!判断与相邻点的能量差，并决定是否改变状态E=COUNT(ISN.NE.NSTATE)RD=RAN(ISEED)IG=NSTATE-IS(IX,JY)DE=E-E0+IG+2.5*RD-1.25!! 用现在状态数组IS1记录边界状态的改变IF(DE.LT.0.0)thensss=sss+1endifenddoif(sss.gt.5)IS1(IX,JY)=11-is(ix,jy)ENDDOenddo!! 每循环20次在显示屏幕上刷新状态(颜色)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,加以区分IF(COUNT(ISN1.NE.IS1(IX,JY)).GT.0)ISRE=SETCOLOR(15)ISRE=SETPIXEL(IX,JY)ENDIFENDDOENDDO!!!记录循环次数和对应的晶粒面积WRITE(1,*)T,SQRT(1.0*COUNT(IS.EQ.1))ENDDOCLOSE(1)END3.2基于CA的品粒长大模拟原理和方法3.2.1元胞选择及其状态确定本文模拟采用四边形元胞，所研究空间为二錐四边形网格平面，元胞邻居为图3-1所示的Moore型邻居.a:dELU图3-1Mooit樂邻居定文毎一牛元胞在甘1时刻的状态白其相邻元胞在t时刻的狀态所决定，即：旳+1)=FS⑴上⑴，附d⑴岸⑴丿◎蓟(啪 <3-1)式中>巩什1)为元胞e在1+1时刻的状态，口⑴〜附)分别为元胞e及其相邹兀•胞在『时刻的状态”F対元砲状态转变规贝J,322元胞状态的转换规则在本文每一模拟时间步C1CAS)内，网格内每个元胞秋态都按臥下规皿判断其足否转换；(1)如元胞憑围的八个元胞■b*s乩仁旨h、i的状态均与其相同「则茫下一个CAS时.元胞亡的状态囉持不变・⑵@)如兀胞匕d.i；h中任意三个同为R状态，贝J注下-个CAS时，元胞亡的状态转变为4(b)ta.元胞血G卧i中任意三个同为/状态，贝恠下一个CAS时，元胞匕的狀态转变为儿⑶假设晶界能均质分右，且一&元胞必須克服能量壁空以达錨其新状态•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)!!根据过去状态IS,定义现在状态为IS1；为了突出边界，特别定义ISN1INTEGERIS1(0:IR+1,0:JR+1),ISN1(1:8)WRITE(*,*)"PLEASEINPUTTHETIMESTEP"READ(*,*)TMAXISEED=RTC()!定义基体体积能为10,晶粒体积能为1IGV=1IGV(0)=10!定义晶粒长大位置，并附能量DOI=1,NMAXIX0=IR*RAN(ISEED)+1JY0=JR*RAN(ISEED)+1IF(IS(IX0,JY0).NE.0)CYCLEIS(IX0,JY0)=IENDDO!! 在循环前定义现在状态数组IS1的初始值IS1=IS!!OPEN(1,FILE="E:\LUKE.DAT")DOT=1,TMAX!!每次循环前重新定义过去状态数组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)!寻找晶粒边界，计算能量，改变状态。!! 整个基体上，各个点上的状态都要根据规则改变，而非随机取点改变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))!如果不是晶粒边界，则跳出，重新循环IF(E0.EQ.0)CYCLE!随机寻找一个相邻点NR=8*RAN(ISEED)+1NSTATE=ISN(NR)!判断与相邻点的能量差，并决定是否改变状态E=COUNT(ISN.NE.NSTATE)RD=RAN(ISEED)IG=IGV(NSTATE)-IGV(IS(IX,JY))DE=E-E0+IG+2.5*RD-1.25!! 用现在状态数组IS1记录边界状态的改变IF(DE.LT.0.0)IS1(IX,JY)=NSTATE!!!! ISRE=SETCOLOR(MOD(IS(IX,JY),15))!! ISRE=SETPIXEL(IX,JY)ENDDOENDDO!!每循环20次在显示屏幕上刷新状态(颜色)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,加以区分IF(COUNT(ISN1.NE.IS1(IX,JY)).GT.0)ISRE=SETCOLOR(15)ISRE=SETPIXEL(IX,JY)ENDIFENDDOENDDO!!!记录循环次数和对应的晶粒面积WRITE(1,*)T,SQRT(1.0*COUNT(IS1.gt.0))ENDDOCLOSE(1)END错误的多晶长大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)!! 根据过去状态IS,定义现在状态为IS1；为了突出边界，特别定义ISN1INTEGERIS1(0:IR+1,0:JR+1),ISN1(1:8)WRITE(*,*)"PLEASEINPUTTHETIMESTEP"READ(*,*)TMAXISEED=RTC()!!!定义基体体积能为10,晶粒体积能为1IGV=1IGV(0)=10!定义晶粒长大位置，并附能量DOI=1,NMAXIX0=IR*RAN(ISEED)+1JY0=JR*RAN(ISEED)+1IS(IX0,JY0)=IENDDO!! 在循环前定义现在状态数组IS1的初始值IS1=IS!!OPEN(1,FILE="E:\LUKE.DAT")DOT=1,TMAX!!每次循环前重新定义过去状态数组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)!! 整个基体上，各个点上的状态都要根据规则改变，而非随机取点改变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))!如果不是晶粒边界，则跳出，重新循环IF(E0.EQ.0)CYCLE!随机寻找一个相邻点sss=0doiii=1,8! NR=8*RAN(ISEED)+1NSTATE=ISN(iii)!判断与相邻点的能量差，并决定是否改变状态E=COUNT(ISN.NE.NSTATE)RD=RAN(ISEED)IG=IGV(NSTATE)-IGV(IS(IX,JY))DE=E-E0+IG+2.5*RD-1.25!! 用现在状态数组IS1记录边界状态的改变! IF(DE.LT.0.0)IS1(IX,JY)=NSTATEIF(DE.LT.0.0)thensss=sss+1nstate1=nstateif(sss.gt.5)IS1(IX,JY)=nstate1endifenddoENDDOENDDO!!每循环20次在显示屏幕上刷新状态(颜色)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,加以区分IF(COUNT(ISN1.NE.IS1(IX,JY)).GT.0)ISRE=SETCOLOR(15)ISRE=SETPIXEL(IX,JY)ENDIFENDDOENDDO!!记录循环次数和对应的晶粒面积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)!!根据过去状态IS,定义现在状态为IS1；为了突出边界，特别定义ISN1INTEGERIS1(0:IR+1,0:JR+1),ISN1(1:8),ISN2(1:8)WRITE(*,*)"PLEASEINPUTTHETIMESTEP"READ(*,*)TMAXISEED=RTC()!!! 定义基体体积能为10,晶粒体积能为1IGV=1IGV(0)=10!定义晶粒长大位置，并附能量DOI=1,NMAXIX0=IR*RAN(ISEED)+1JY0=JR*RAN(ISEED)+1IS(IX0,JY0)=IENDDO!!在循环前定义现在状态数组IS1的初始值IS1=IS!!OPEN(1,FILE="E:\LUKE.DAT")DOT=1,TMAX!!每次循环前重新定义过去状态数组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)!! 整个基体上，各个点上的状态都要根据规则改变，而非随机取点改变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))!如果不是晶粒边界，则跳出，重新循环IF(E0.EQ.0)CYCLE!随机寻找一个相邻点iiii=0doiii=1,8! NR=8*RAN(ISEED)+1NSTATE=ISN(iii)!判断与相邻点的能量差，并决定是否改变状态E=COUNT(ISN.NE.NSTATE)RD=RAN(ISEED)IG=IGV(NSTATE)-IGV(IS(IX,JY))DE=E-E0+IG+2.5*RD-1.25!! 用现在状态数组IS1记录边界状态的改变! 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!!每循环20次在显示屏幕上刷新状态(颜色)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,加以区分IF(COUNT(ISN1.NE.IS1(IX,JY)).GT.0)ISRE=SETCOLOR(15)ISRE=SETPIXEL(IX,JY)ENDIFENDDOENDDO!!!记录循环次数和对应的晶粒面积WRITE(1,*)T,COUNT(IS1.GT.0)ENDDOCLOSE(1)ENDA23局部规则的确定局部規则在建立元胞自动机模型中是最遁要的步骤，温度是如何传递的是根据这个规则而制定的，本文采用的是Wn.Neumann型的邻居，所以研究中心元胞周围的四个元胞对它温度的影响・用热力学来研究元胞之间温度的传递是局部规则制定的基础.物体内温度的变化与热量有关.热量的增减引起温度的升降，为了确定温度场，必须研究热量在物体内传导的过程。温度在时间域和空间域中的分布，称为温度场.它可表示为公式(4-1》T=T(x,齐為t) (4J)若温度不随时间变化，即#=0,则T-T(k,y,z)称为稳定温度场.若温度沿z向不变，则称为平面温度场。在任一瞬时.连接场内相同温度值的各点，就得到此时刻的等温面。沿等温面切向,温度不变口而沿垂直等温面的法向，温度变化率最大。表示一点最大埴温率的矢量，称为温度梯度.即公式(42)：“st.arjrdTt.dr(4-2)V"%需知云r齐+无瓦其中血是单位矢量，沿等温面的法线指向增温方向。(4-2)在单位时问内通过单位面积的热JL称为热渍密度，即公式(4-3)；IdQ(切dH——旦

sdt(切菓■点的最大热流密度矢量是沿等温面的法线且指向降温方向「根韬传导定律，热流密度与温度梯度成正比而方问相反，即公式(M)：(4-4)(4-5)率：忑耍为x方向上的温度梯度.我们可以得到导热方程为公式(46)(4-4)(4-5)率：忑耍为x方向上的温度梯度.我们可以得到导热方程为公式(46)：OX式中p为材料的密度(kg/M)：c为材料的比热容(J/(kgK))；q=—AVT=-%久2~dn几称为导热系数。热流密度在x方向的分量是公式件刃；q$=-Z—cosC^x)=-A—^x.y.z)dn dx类似，在任意方向的热流密度，等于导热系数乘以温度在该方向的变化率的负值&其中，亦为x方向上的热流密度；九为材料沿K,y,z方向的热导r为时间G)；人，九人分别为材料沿x,y,z方向的热导率(W/(mK))；Q二0(x,Z玷)为材料内部的热源密度(W/kg)。上式中，第一项为体元升温需要的热量：右侧第一、二和三项是由x,y和2方向流入体元的热量：躍后一项体元内热源产生的热量。微分方程的物理意义，体元升温所需的热量应该等于流入体元的热量与体元内产生的热量的总和。由于我们研究的是二维的温度场模型，所以简化后的方程，二维非稳态热传导方程为公式(4・〃(4-7)(4-7)在区域内进行网格划分，将连续的求解域离散为不连续的点，形成离散网格,网格步长分别为兀何一九=Ar,格,网格步长分别为兀何一九=Ar,畑_九=切划分的单元格步长可以是均匀的，也可以是不均匀的。此模型我们取单元格步长是均匀的，如3所表示:对于二维各向同性，无内热源的非稳态热传导微分方程，利用中心差分法来研究，差分解法是将函数用离散点的函数值表示，将微分用有限差分表示,将导数用有限差商表示，从而将微分方程化为差分方程一代数方程，使求解微分方程的问题化为求解代数方程的问题.中心差分法公式件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 如,在划分单元时，我们令则上式可化为公式(4-12).A/4+4*4/供<Ar)2J(4-11)(4-12)上式即为元胞在绝热条件下的温度传递公式・即可以通过中心元胞的温度値和邻居元胞的温度值得到下一时刻该元胞的温度“我们把公式(・12)定义为元胞自动机的局部规则口通过对元胞自动机温度场模型中元胞，元胞空间，邻居规则「局部规则的建立宪成，就刀片温度场模型建立完咸，我们把离散的网格定文为元胞，元胞空间为二维p邻居来用VbmNeumann型，局部规则为推导公式件1幼!相变byUsingBoundaryTrackingMethord,2008,11.23!graingrowthvelocityis"1.0*(IT-NT(IST))"or