《蒙特卡罗模拟在材料科学中应用举例课件》由会员分享，可在线阅读，更多相关《蒙特卡罗模拟在材料科学中应用举例课件（25页珍藏版）》请在金锄头文库上搜索。
1、蒙特卡罗模拟在材料科学中应用举例,扩散 晶粒形核与长大 再结晶,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(*,*) 实验天数
2、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) E
3、ND,随机行走：原子扩散的平均距离与原子跳动次数的平方根成正比。,变量的定义： 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
4、.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 !第几
5、步跳跃 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,
6、 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
7、,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 & Gra
8、in 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 gra
9、ins 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
10、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. R
11、epeat 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 latti
12、ce 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(
13、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=SETPI
14、XEL(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),任意选邻居再计算能量,随机选取一个邻居
15、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、寻找晶
16、界 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 TIM
