蒙特卡罗模拟在材料科学中应用举例

蒙特卡罗模拟在材料科学中应用举例,扩散 晶粒形核与长大 再结晶,1、Monte Carlo Simulation of Diffusion,Mechanism : Rand Walk 方均位移平方,!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(*,*) “实验天数Jt，实验次数 Ic“ READ(*,*) Jt, Ic ISEED=RTC() DO J=1,Jt !第几天实验 X=0 ! 每天都是从原点出发 DO I=1,Ic !第几步跳跃 RN=RAN(ISEED) IF(RN0.5)THEN X=X+1 ELSE X=X-1 END IF XX(J,I)=X*X !记录下原子每天每次跳动后离原点的距离 END DO END DO OPEN(1,FILE=“f:DIF1.DAT“) DO I=1,Ic XXM=0.0 XXM(I)=1.0*SUM(XX(1:Jt,I)/Jt ! WRITE(1,*) I, XXM(I) END DO CLOSE(1) END,随机行走：原子扩散的平均距离与原子跳动次数的平方根成正比。,变量的定义： X记录原子的位置 xx( j, i)记录原子每次跳动后距离原点的距离 j ：统计的天数，i：跳动的次数 xxm(i)统计距离与次数的关系 实验天数：Jt，实验次数：Ic,二维随机行走随机的确定,方法一、 RN=RAN(ISEED) IF(RN.LT.0.25)THEN x=x y=y-1 END IF IF(RN.LT.0.5.AND.RN.GE.0.25)THEN x=x y=y+1 END IF IF(RN.LT.0.75.AND.RN.GE.0.5)THEN x=x-1 y=y END IF IF(RN.GE.0.75)THEN x=x+1 y=y END IF,方法二 XN=(/0,0,-1,1/) YN=(/-1,1,0,0/) RN=4*RAN(ISEED)+1 X=X+YN(RN) Y=Y+YN(RN),! Monte Carlo Simulation of Two Dimensional Diffusion INTEGER X,Y,XY(1:1000,1:1000) REAL XYM(1:1000) WRITE(*,*) “实验天数Jt，实验次数 Ic“ READ(*,*) Jt,Ic ISEED=RTC() DO J=1,Jt !第几天实验 X=0 ! Y=0 ! DO I=1,Ic !第几步跳跃 RN=RAN(ISEED) IF(RN.LT.0.25)THEN x=x y=y-1 else IF(RN.LT.0.5.AND.RN.GE.0.25)THEN x=x y=y+1 else IF(RN.LT.0.75.AND.RN.GE.0.5)THEN x=x-1 y=y else IF(RN.GE.0.75)THEN x=x+1 y=y END IF XY(J,I)=X*X+Y*Y END DO END DO OPEN(1,FILE=“f:DIF2.DAT“) DO I=1,Ic XYM=0.0 XYM(I)=1.0*SUM(XY(1:Jt,I)/Jt ! WRITE(1,*) I, XYM(I) END DO CLOSE(1) END,! Monte Carlo Simulation of Two Dimensional Diffusion INTEGER X,XY(1:1000,1:1000),y,XN(1:4),YN(1:4),RN REAL XYM(1:1000) ! X:INSTANTANEOUS POSITION OF ATOM ! XY(J,I)：X*Y ,J:第几天实验,I:第几步跳跃 ! XYM(I): THE MEAN OF XY WRITE(*,*) “实验天数Jt，实验次数 Ic“ READ(*,*) Jt,Ic XN=(/0,0,-1,1/) YN=(/-1,1,0,0/) ISEED=RTC() DO J=1,Jt !第几天实验 X=0 ! Y=0 ! DO I=1,Ic !第几步跳跃 RN=4*RAN(ISEED)+1 X=X+YN(RN) Y=Y+YN(RN) XY(J,I)=X*X+Y*Y END DO END DO OPEN(1,FILE=“C:DIF2.DAT“) DO I=1,Ic XYM=0.0 XYM(I)=1.0*SUM(XY(1:Jt,I)/Jt ! WRITE(1,*) I, XYM(I) END DO CLOSE(1) ！做三维空间随机行走？ END,Simulation of Recrystallisation & Grain Growth by Means of a Monte Carlo Model,Model Microstructure The grain structure is mapped onto a discrete two or three dimensional lattice. To each lattice site an integer S between 1 and a maximum value Q is assigned. These integers represent the crystallographic orientation of the different grains and are called “spins“ or “orientations“. Thus between two adjacent lattice sites of unlike spins there is a grain boundary segment, while two neighboring sites of like spins belong to the same grain and there is no grain boundary between them.,2 Algorithm (proposed by Anderson et al. 1984),1. Pick up a lattice site i with spin S_i at random 2. Calculate its energy E_i according to equation (1) 3. Pick up a new spin S_j different from the old one 4. Calculate the energy of the chosen site i if it had this new spin S_j 5. Flip site i to spin S_j with probability P according to equation (2) 6. Repeat steps 1.-5.,Energy E_i of a lattice site: (1) J: Energy factor NN: Number of nearest neighbours : Kronecker-delta,Flip Probability P: (2) : Energy difference due to the flip ( ) Kb : Boltzmann k T: System temperature Time Unit: One time unit (1 Monte Carlo step: MCS) equals the number of lattice sites of the system,2.蒙特卡罗模拟双晶粒的长大,如何表示晶粒：颜色（状态）相同的像素（元胞）； 如何定义一段晶界：不同状态元胞的交线； 如何表示一段晶界：形成晶界的两个元胞； 如何定义邻居：最近邻，次近邻；,画一个圆晶粒,USE MSFLIB INTEGER XR，YR !在的区域中画一个圆 PARAMETER XR=400,YR=400 INTEGER R,S(1:XR,1:YR) X0=XR/2 ! 圆心位置X0，YO Y0=YR/2 R=MIN(X0-10,Y0-10) !圆半径 S=0 !像素的初始状态（颜色） DO I=1,XR DO J=1,YR IF(I-X0)*2+(J-Y0)*2=R*2) S(I,J)=10 IER=SETCOLOR(S(I,J) IER=SETPIXEL(I,J) END DO END DO END,画晶界,XN=(/-1,0,1,-1,1,-1,0,1/) YN=(/-1,0,-1,0,0,1,1,1/) DO I=1,XR !画晶界 DO J=1,YR NDS=0 DO K=1,8 IF(S(I,J).NE.S(I+XN(K),J+YN(K)NDS=NDS+1 END DO IF(NDS0)THEN IER=SETCOLOR(9) ELSE IER=SETCOLOR(S(I,J) END IF IER=SETPIXEL(I,J) END DO END DO,转变规则能量最小原理,如何实现界面的迁移：单胞状态的转变 转变规则：随机选择一个邻居，如果转变后系统能量降低（考虑能量起伏）。 如何计算能量：体积能（动力）；界面能（阻力）： 体积能EV计算： EV=0 EV(0)=ER,界面能,如何计算一个单胞界面能：异类邻居数之和。 XN=(/-1,-1,-1,0,0,1,1,1/) YN=(/-1,0,1,-1,1,-1,0,1/) DO II=1,8 ISB(II)=IS(I+XN(II),J+YN(II) END DO E0=COUNT(ISB.NE.IS(I,j),任意选邻居再计算能量,随机选取一个邻居CELL 及能量变化 ISTR=ISB(8*RAN(ISEED)+1) IF(ISTR.EQ.IS)CYCLE ETR=COUNT(ISB.NE.ISTR),能量判断,单个元胞的体积能Ev与元胞的一个面的能量Es之比： Ev=8Es 能量变化与能量起伏 DEB=ETR-E0 DEV=EV(ISTR)-EV(IS) DE=DEB+DEV+2.5*RAN(ISEED)-1.25,MC模拟双晶粒的长大,如何定义基体和晶粒； 如何寻找边界； 如何计算能量（构造或描述问题的概率过程） 步骤： 1、读取晶粒生长步长； 2、在基体上画一个原始晶粒，并赋予基体和原始晶粒不同的状态（取向号或体积能） 3、寻找晶界 4、计算晶界上网格点的相互作用，每个格点的相互作用导致晶粒长大或缩小： 随机选取一个初始格点； 若此点属于晶界，那么可以随机转变为它邻居的状态，若不是晶界，则跳出循环，不发生转变； 计算转变前后的能量变化，E=Ev+Es +Eq (自由能=体积能+表面能+能量起伏) 若E小于0，则新晶相被接受，网格状态发生转变。,双晶长大/收缩,USE MSFLIB PARAMETER IR=400,JR=400 INTEGER IS(0:IR+1,0:JR+1),TMAX,ISB(1:8),ISTR,T,NR,IX,IY,X,Y,RD INTEGER XN(1:8),YN(1:8) XN=(/-1,0,1,-1,1,-1,0,1/) YN=(/-1,0,-1,0,0,1,1,1/) WRITE(*,*)“PLEASE INPUT THE TIME STEP “ READ(*,*)TMAX ISEED=RTC() ! 定义圆心和半径 X0=IR/2 Y0=JR/2 R=MIN(X0,Y0)-100 ! 定义基体和圆晶粒分别为状态5、状态2 IS=5 DO I=1,IR DO J=1,JR IF(I-X0)*2+(J-Y0)*2=R*2) IS(I,J)=2 IER=SETCOLOR(IS(I,J) IER=SETPIXEL(I,J) END DO END DO OPEN(1,FILE=“E:data.DAT“) ! 寻找晶粒边界，计算能量，改变状态。 DO T=1,TMAX DO X=1,IR DO Y=1,JR IX=IR*RAN(ISEED)+1 JY=JR*RAN(ISEED)+1 DO II=1,8 ISB(II)=IS(IX+XN(II),JY+YN(II) END DO E0=COUNT(ISB.NE.IS(IX,JY) ! 如果不是圆晶粒边界，则跳出，重新循环 IF(E0.EQ.0)CYCLE ! 随机寻找一个相邻点 NR=8*RAN(ISEED)+1 ISTR=ISB(NR) ! 判断与相邻点的能量差，并决定是否改变状态 ETR=COUNT(ISB.NE.ISTR) RD=RAN(ISEED) DE=ETR-E0+ISTR-IS(IX,JY)+2.5*RD-1.25 IF(DE.LT.0.0)IS(IX,JY)=ISTR ISRE=SETCOLOR(IS(IX,JY) ISRE=SETPIXEL(IX,JY) ENDDO ENDDO WRITE(1,*)T,SQRT(1.0*COUNT(IS.EQ.5) ENDDO CLOSE(1) END,多晶的长大,与单晶相比较 多个晶粒的定义 多个晶粒体积能的定义 周期性边界条件,! 定义晶粒位置，并附能量 DO I=1, NMAX IX0=IR*RAN(ISEED)+1 JY0=JR*RAN(ISEED)+1 IS(IX0,JY0)=I END DO,!定义基体体积能为8,晶粒体积能为1。 IGV=1 IGV(0)=8,IS(0,1:JR)=IS(Ic,1:JR) IS(IR+1,1:JR)=IS(1,1:JR) IS(0:IR+1,0)=IS(0:IR+1,JR) IS(0:IR+1,JR+1)=IS(0:IR+1,1),USE MSFLIB PARAMETER IR=400,JR=400,NMAX=1000 INTEGER IS(0:IR+1,0:JR+1),TMAX,ISB(1:8),ISTR,T,NR,IX,IY,X,Y,RD INTEGER XN(1:8),YN(1:8), IGV(0:NMAX) XN=(/-1,0,1,-1,1,-1,0,1/) YN=(/-1,0,-1,0,0,1,1,1/) WRITE(*,*)“PLEASE INPUT THE TIME STEP “ READ(*,*)TMAX ISEED=RTC() IGV=1! 定义基体体积能为8,晶粒体积能为1 IGV(0)=8 DO I=1, NMAX ! 定义晶粒长大位置，并附能量 IX0=IR*RAN(ISEED)+1 JY0=JR*RAN(ISEED)+1 IS(IX0,JY0)=I END DO ! OPEN(1,FILE=“E:MULTI.DAT“) ! 寻找晶粒边界，计算能量，改变状态。 DO T=1,TMAX ! 边界条件 IS(0,1:JR)=IS(Ic,1:JR) IS(IR+1,1:JR)=IS(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 X=1,IR DO Y=1,JR IX=IR*RAN(ISEED)+1 JY=JR*RAN(ISEED)+1 DO II=1,8 ISB(II)=IS(IX+XN(II),JY+YN(II) END DO E0=COUNT(ISB.NE.IS(IX,JY) ! 如果不是圆晶粒边界，则跳出，重新循环