计算材料学Fortran程序汇总.doc

INTEGER M(1:10000), NUMBER1(0:364), NUMBER2 REAL X,Y ISEED=RTC()DO J=1, 10000NUMBER1=0X=RAN(ISEED)NUMBER1(0)=INT(365*X+1)JJJ=1 DO I=1,364 Y=RAN(ISEED) NUMBER2=INT(365*Y+1) ETR=COUNT(NUMBER1.EQ.NUMBER2) IF (ETR= =1) THEN EXIT ELSE JJJ=JJJ+1 M(J)=JJJ NUMBER1(I)=NUMBER2 END IF END DOEND DO DO I=1,10000 IF(M(I).LE.23) SUM=SUM+1END DO PRINT *,SUM/10000 END =Monte Carlo Simulation of One Dimensional Diffusion INTEGER X,XX(1:1000,1:1000)REAL XXM(1:1000)!X:INSTANTANEOUS POSITION OF ATOM!XX(J,I)：X*X ,J:第几天实验,I:第几步跳跃!XXM(I): THE MEAN OF XX WRITE(*,*) 实验天数JMAX，实验次数 IMAXREAD(*,*) JMAX,IMAXISEED=RTC()DO J=1,JMAX !第几天实验X=0 !DO I=1,IMAX !第几步跳跃RN=RAN(ISEED)IF(RN0.5)THEN X=X+1ELSE X=X-1 END IFXX(J,I)=X*XEND DOEND DOOPEN(1,FILE=C:DIF1.DAT)DO I=1,IMAXXXM=0.0XXM(I)=1.0*SUM(XX(1:JMAX,I)/JMAX !WRITE(1,*) I, XXM(I)END DOCLOSE(1)END=!Monte Carlo Simulation of Two Dimensional Diffusion INTEGER X,Y,XY(1:1000,1:1000)REAL XYM(1:1000)!X:INSTANTANEOUS POSITION OF ATOM!XY(J,I)：X*Y ,J:第几天实验,I:第几步跳跃!XYM(I): THE MEAN OF XY WRITE(*,*) 实验天数JMAX，实验次数 IMAXREAD(*,*) JMAX,IMAXISEED=RTC()DO J=1,JMAX !第几天实验X=0 !Y=0 !DO I=1,IMAX !第几步跳跃RN=RAN(ISEED)IF(RN.LT.0.25)THEN x=xy=y-1 END IFIF(RN.LT.0.5.AND.RN.GE.0.25)THEN x=xy=y+1 END IFIF(RN.LT.0.75.AND.RN.GE.0.5)THEN x=x-1y=y END IFIF(RN.GE.0.75)THEN x=x+1y=y END IFXY(J,I)=X*X+Y*YEND DOEND DOOPEN(1,FILE=C:DIF2.DAT)DO I=1,IMAXXYM=0.0XYM(I)=1.0*SUM(XY(1:JMAX,I)/JMAX !WRITE(1,*) I, XYM(I)END DOCLOSE(1)END=!Monte Carlo Simulation of One Dimensional Diffusion INTEGER X,XY(1:1000,1:1000),y,XN(1:4),YN(1:4),RNREAL XYM(1:1000)!X:INSTANTANEOUS POSITION OF ATOM!XY(J,I)：X*Y ,J:第几天实验,I:第几步跳跃!XYM(I): THE MEAN OF XY WRITE(*,*) 实验天数JMAX，实验次数 IMAXREAD(*,*) JMAX,IMAXXN=(/0,0,-1,1/)YN=(/-1,1,0,0/)ISEED=RTC()DO J=1,JMAX !第几天实验X=0 !Y=0 !DO I=1,IMAX !第几步跳跃RN=4*RAN(ISEED)+1X=X+YN(RN) Y=Y+YN(RN)XY(J,I)=X*X+Y*YEND DOEND DOOPEN(1,FILE=C:DIF2.DAT)DO I=1,IMAXXYM=0.0XYM(I)=1.0*SUM(XY(1:JMAX,I)/JMAX !WRITE(1,*) I, XYM(I)END DOCLOSE(1)END做三维空间随机行走？留作业 !Monte Carlo Simulation of One Dimensional Diffusion INTEGER X,XY(1:1000,1:1000),y,XN(1:6),YN(1:6),ZN(1:6),RNREAL XYM(1:1000)!X:INSTANTANEOUS POSITION OF ATOM!XY(J,I)：X*Y ,J:第几天实验,I:第几步跳跃!XYM(I): THE MEAN OF XY WRITE(*,*) 实验天数JMAX，实验次数 IMAXREAD(*,*) JMAX,IMAXXN=(/0,0,-1,1,0,0/)YN=(/-1,1,0,0,0,0/)ZN=(/0,0,0,0,1,-1/)ISEED=RTC()DO J=1,JMAX !第几天实验X=0 !Y=0 !Z=0DO I=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*ZEND DOEND DOOPEN(1,FILE=C:DIF2.DAT)DO I=1,IMAXXYM=0.0XYM(I)=1.0*SUM(XY(1:JMAX,I)/JMAX !WRITE(1,*) I, XYM(I)END DOCLOSE(1)END=通过该程序了解fortran语言如何画图（通过像素画图）USE MSFLIB INTEGER XR，YR !在的区域中画一个圆PARAMETER XR=400,YR=400 INTEGER R, S(1:XR,1:YR) X0=XR/2 ! 圆心位置X0，YOY0=YR/2R=MIN(X0-10,Y0-10) !圆半径S=0 !像素的初始状态（颜色）DO I=1,XRDO J=1,YRIF(I-X0)*2+(J-Y0)*2=R*2)S(I,J)=10IER=SETCOLOR(S(I,J) IER=SETPIXEL(I,J)END DOEND DOEND=!画一个圆（1、如何选出晶界区域；2、进一步加深对画图的理解）USE MSFLIB INTEGER XR，YR !在的区域中画一个圆PARAMETER XR=400,YR=400 INTEGER R, S(0:XR+1,0:YR+1), XN(1:4), YN(1:4), SNS XN=(/0,0,-1,1/)YN=(/-1,1,0,0/)X0=XR/2 ! 圆心位置X0，Y0Y0=YR/2R=MIN(X0-10,Y0-10) !圆半径S=0 !像素的初始状态（颜色）DO I=1,XRDO J=1,YRIF(I-X0)*2+(J-Y0)*20)THENIER=SETCOLOR(9)ELSEIER=SETCOLOR(8)END IF IER=SETPIXEL(I,J)END DOEND DOEND如何画有一定宽度的晶界？=！MC模拟一个晶粒的缩小（1、 如何定义基体和晶粒；（2、 如何寻找边界；（3、 如何计算能量（构造或描述问题的概率过程） 步骤：1、在基体上画一个原始晶粒，并赋予基体和原始晶粒不同的状态（取向号或体积能） 2、寻找晶界 3、MC的一个时间步（晶粒长大过程中一种随机性） 4、计算晶界上网格点的相互作用，每个格点的相互作用导致晶粒长大 随机选取一个初始格点； 若此点属于晶界，那么可以随机转变为它邻居的状态，若不是晶界，则跳出循环，不发生转变； 计算转变前后的能量变化，E=Ev+Es +Eq(自由能=体积能+表面能+能量起伏) 若E小于0，则新晶相被接受，网格状态发生转变。USE MSFLIBPARAMETER IR=400,JR=400 INTEGER IS(0:IR+1,0:JR+1),TMAX,ISN(1:8),NSTATE,T,NR,IX,IYWRITE(*,*)PLEASE INPUT THE TIME STEP READ(*,*)TMAXISEED=RTC()!定义圆心和半径IRC=IR/2JRC=JR/2R=MIN(IRC,JRC)-10!定义基体和圆晶粒分别为状态1、状态2IS=1 DO I=1,IR DO J=1,JR DISTANCE=SQRT(1.0*(I-IRC)*2+1.0*(J-JRC)*2) IF(DISTANCE.LT.R)IS(I,J)=2 ISE=SETCOLOR(IS(I,J) ISE=SETPIXEL(I,J) END DO END DOOPEN(1,FILE=E:LUKE.DAT)!寻找晶粒边界，计算能量，改变状态。 DO T=1,TMAX DO X=1,IR DO Y=1,JR IX=IR*RAN(ISEED)+1 JY=JR*RAN(ISEED)+1 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)+1 NSTATE=ISN(NR)!判断与相邻点的能量差，并决定是否改变状态 E=COUNT(ISN.NE.NSTATE) RD=RAN(ISEED) DE=E-E0+NSTATE-IS(IX,JY)+2.5*RD-1.25 IF(DE.LT.0.0)IS(IX,JY)=NSTATE ISRE=SETCOLOR(IS(IX,JY) ISRE=SETPIXEL(IX,JY) ENDDO ENDDO WRITE(1,*)T,COUNT(IS.EQ.2)ENDDO CLOSE(1)END （1、 演示体积能互换的转变情况；2、 演示IX=IR*RAN(ISEED)+1 JY=JR*RAN(ISEED)+1的替换；3、能量关系变化；）=作业3：基体上单晶粒形核长大（1、 选取初始形核点（2、 定义体积能的变化） USE MSFLIBPARAMETER IR=400,JR=400 INTEGER IS(0:IR+1,0:JR+1),TMAX,ISN(1:8),NSTATE,T,NR,IX,IYWRITE(*,*)PLEASE INPUT THE TIME STEP READ(*,*)TMAXISEED=RTC()!定义圆心和半径IRC=IR/2JRC=IR/2 !定义基体和圆晶粒分别为状态10、状态2IS=10IS(IRC,JRC)=2OPEN(1,FILE=E:LUKE.DAT)!寻找晶粒边界，计算能量，改变状态。 DO T=1,TMAX DO X=1,IR DO Y=1,JR IX=IR*RAN(ISEED)+1 JY=JR*RAN(ISEED)+1 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)+1 NSTATE=ISN(NR)!判断与相邻点的能量差，并决定是否改变状态 E=COUNT(ISN.NE.NSTATE) RD=RAN(ISEED) DE=E-E0+NSTATE-IS(IX,JY)+2.5*RD-1.25 IF(DE.LT.0.0)IS(IX,JY)=NSTATE ISRE=SETCOLOR(IS(IX,JY) ISRE=SETPIXEL(IX,JY) ENDDO ENDDO WRITE(1,*)T,SQRT(1.0*COUNT(IS.EQ.2)ENDDOCLOSE(1)END 基体上多晶粒形核长大，（1、 如何生成多晶粒形核核心（2、 如何定义、并区别多晶粒3、为什么要有边界条件） 步骤：1、定义变量：体积能变量（一维数组）；基体各像素点变量（二维数组），邻居的取向号（一维数组）等；2、读取晶粒生长步长； 3、寻找晶界 4、计算晶界上网格点的相互作用，每个格点的相互作用导致晶粒长大或缩小： 随机选取一个初始格点； 若此点属于晶界，那么可以随机转变为它邻居的状态，若不是晶界，则跳出循环，不发生转变； 计算转变前后的能量变化，E=Ev+Es +Eq(自由能=体积能+表面能+能量起伏) 若E小于0，则新晶相被接受，网格状态发生转变。= USE MSFLIB! 定义常数名PARAMETER IR=400,JR=400,NMAX=100! 定义变量：体积能变量（一维数组）；基体各像素点变量（二维数组），邻居的取向号（一维数组）等； INTEGER IS(0:IR+1,0:JR+1),TMAX,ISN(1:8),NSTATE,T,NR,IX0,IY0,IX,IYINTEGER IGV(0:NMAX)! 读取晶粒生长步长；WRITE(*,*)PLEASE INPUT THE TIME STEP READ(*,*)TMAXISEED=RTC()!定义基体体积能为10,所有晶粒体积能为1IGV=1IGV(0)=10!定义晶粒长大位置（100个形核点初始形核位置），并附已不同的取向号DO I=1, NMAXIX0=IR*RAN(ISEED)+1JY0=JR*RAN(ISEED)+1 IS(IX0,JY0)=IEND DO!边界条件 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)!寻找晶粒边界。 DO T=1,TMAX DO X=1,IR DO Y=1,JR IX=IR*RAN(ISEED)+1 JY=JR*RAN(ISEED)+1 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)+1 NSTATE=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.25 IF(DE.LT.0.0)IS(IX,JY)=NSTATE ISRE=SETCOLOR(MOD(IS(IX,JY),15) ISRE=SETPIXEL(IX,JY) ENDDO ENDDO WRITE(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 方法模拟显得尤为重要；5、异常晶粒的长大；6、三维晶粒的长大（模拟不同的点阵结构生长过程）；山东大学 中国有色金属学报，2008，47、蒙特卡洛方法不是数学的方法，而是一种多次模拟实验的方法。=元胞自动机程序（生命永不停止）use msflibPARAMETER IR=400,JR=400,NMAX=40000 !NMAX-随机产生的生命种子integer is(0:1001,0:1001),is1(0:1001,0:1001),ISN(1:8), Tmax, num ！is-基体的二维数组print*,please input loop(Tmax)read*, TmaxISEED=RTC()is=15 ！“死”的状态，基体为白色 ！赋予生命的种子，“活”的状态1DO I=1, NMAXIX0=IR*RAN(ISEED)+1JY0=JR*RAN(ISEED)+1 IS(IX0,JY0)=1END DO!execute the ruledo t=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)！搜索生命存在的位置do ix=1,IRdo jy=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) then is1(ix,jy)=1 else is1(ix,jy)=15 end if ！画图 ISRE=SETCOLOR(IS1(IX,JY) ISRE=SETPIXEL(IX,JY) end do end dois=is1end doend=元胞自动机产生和发展. 四个阶段: 1940s 诞生：Von Neumann 自我复制机. 1960-70s起步：JH.Conway 生命游戏. 1980s 理论研究：S.Wolfram CA分类. 1980-90s 应用：HPP-FHP格子气自动机. C.Langton N.Packard 人工生命元胞自动机（Cellular Automata，简称CA，也有人称其为细胞自动机、点格自动机、分子自动机或单元自动机）是一种建立在离散的时间和空间上的动力学系统。散布在规则格网（Lattice Grid）中的每一元胞（Cell）取有限的离散状态，遵循同样的作用规则，依据确定的局部规则作同步更新。大量元胞通过简单的相互作用而构成动态系统的演化。与一般的动力学模型不同，元胞自动机不是由严格定义的物理方程或函数确定，而是由一系列模型构造的规则构成，凡是满足这些规则的模型都可以算做是元胞自动机模型。因此，元胞自动机是一类模型的总称，或者是一个方法框架。其特点是时间、空间、状态都是离散的，每个变量只取有限多个状态，且其状态改变的规则在时间和空间上都是局部的。 元胞自动机是一种对具有局域连通性的格点，应用局部（有时为中等范围）确定性或概论性的转换规则来描述在离散空间和时间上复杂系统演化规律的同步算法。元胞自动机可以采取任意的维数，自动机的格点描述被认为与所选模型密切相关的基本的系统实体。独立的格点可以表示连续的体积单元，原子颗粒大原子团、汽车、人口、晶格缺陷等。通过作用于每个格点的确定性或概率性转换规则，自动机的演化得以实施。这些规则决定了每个元胞的状态，并将其描述为其前一状态和邻近元胞的函数。在计算状态转变过程中所考虑的邻近位置的数目确定了相互作用的范围和自动机的局部演化。！转变概率元胞自动机特点：时间、空间、状态都离散的动力学模型，利用确定性或概论性的转换规则来描述在离散空间和时间上复杂系统演化 =元胞自动机模拟单晶长大 USE MSFLIBPARAMETER IR=400,JR=400 INTEGER IS(0:IR+1,0:JR+1),TMAX,ISN(1:8),NSTATE,T,NR,IX0,IY0,IX,JY! 根据过去状态IS，定义现在状态为IS1；为了突出边界，特别定义ISN1 INTEGER IS1(0:IR+1,0:JR+1),ISN1(1:8) WRITE(*,*)PLEASE INPUT THE TIME STEP READ(*,*)TMAXISEED=RTC() IRC=IR/2JRC=JR/2!定义基体体积能为10,晶粒体积能为1IS=10IS(IRC,JRC)=1! 在循环前定义现在状态数组IS1的初始值 IS1=IS OPEN(1,FILE=E:LUKE.DAT) DO T=1,TMAX! 每次循环前重新定义过去状态数组IS IS=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)! 整个基体上，各个点上的状态都要根据规则改变，而非随机取点改变 DO IX=1,IR DO JY=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)+1 NSTATE=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)=NSTATE ENDDO END DO! 每循环20次在显示屏幕上刷新状态（颜色） DO IX=1,IR DO JY=1,JR! IF(MOD(T,20).EQ.0)THEN ISN1=(/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) ! END IF ENDDO ENDDO! 记录循环次数和对应的晶粒面积 WRITE(1,*)T, SQRT(1.0*COUNT(IS.EQ.1)ENDDOCLOSE(1)END =元胞自动机模拟单晶长大，再附加一个规则 USE MSFLIBPARAMETER IR=400,JR=400,NMAX=100 INTEGER IS(0:IR+1,0:JR+1),TMAX,ISN(1:8),NSTATE,T,NR,IX0,IY0,IX,JYINTEGER IGV(0:NMAX)! 根据过去状态IS，定义现在状态为IS1；为了突出边界，特别定义ISN1 INTEGER IS1(0:IR+1,0:JR+1),ISN1(1:8) WRITE(*,*)PLEASE INPUT THE TIME STEP READ(*,*)TMAXISEED=RTC() IRC=IR/2JRC=JR/2!定义基体体积能为10,晶粒体积能为1IS=10IS(IRC,JRC)=1! 在循环前定义现在状态数组IS1的初始值 IS1=IS! OPEN(1,FILE=E:LUKE.DAT) DO T=1,TMAX! 每次循环前重新定义过去状态数组IS IS=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)! 整个基体上，各个点上的状态都要根据规则改变，而非随机取点改变 DO IX=1,IR DO JY=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=0do iii=1,8 NSTATE=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+1 end ifend do if(sss.gt.5) IS1(IX,JY)=11-is(ix,jy) ENDDOend do! 每循环20次在显示屏幕上刷新状态（颜色） DO IX=1,IR DO JY=1,JR IF(MOD(T,20).EQ.0)THEN ISN1=(/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) END IF ENDDO ENDDO! 记录循环次数和对应的晶粒面积 WRITE(1,*)T, SQRT(1.0*COUNT(IS.EQ.1)ENDDOCLOSE(1)END USE MSFLIBPARAMETER IR=400,JR=400,NMAX=100 INTEGER IS(0:IR+1,0:JR+1),TMAX,ISN(1:8),NSTATE,T,NR,IX0,IY0,IX,JYINTEGER IGV(0:NMAX)! 根据过去状态IS，定义现在状态为IS1；为了突出边界，特别定义ISN1 INTEGER IS1(0:IR+1,0:JR+1),ISN1(1:8) WRITE(*,*)PLEASE INPUT THE TIME STEP READ(*,*)TMAXISEED=RTC() !定义基体体积能为10,晶粒体积能为1IGV=1IGV(0)=10!定义晶粒长大位置，并附能量DO I=1,NMAX IX0=IR*RAN(ISEED)+1JY0=JR*RAN(ISEED)+1IF(IS(IX0,JY0).NE.0)CYCLE IS(IX0,JY0)=IEND DO! 在循环前定义现在状态数组IS1的初始值 IS1=IS! OPEN(1,FILE=E:LUKE.DAT) DO T=1,TMAX! 每次循环前重新定义过去状态数组IS IS=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)!寻找晶粒边界，计算能量，改变状态。 ! 整个基体上，各个点上的状态都要根据规则改变，而非随机取点改变 DO IX=1,IR DO JY=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)+1 NSTATE=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) ENDDO ENDDO! 每循环20次在显示屏幕上刷新状态（颜色） DO IX=1,IR DO JY=1,JR IF(MOD(T,20).EQ.0)THEN ISN1=(/IS1