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

蒙特卡罗模拟在材料科学中应用举例扩散晶粒形核与长大再结晶1、MonteCarloSimulationofDiffusionMechanism:RandWalk方均位移平方max122)(max1max),...,3,2,1(max),...,3,2,1(jjjijiijiiiiiXXjXjjXXXXiiX!MonteCarloSimulationofOneDimensionalDiffusionINTEGERX,XX(1:1000,1:1000)REALXXM(1:1000)!X:INSTANTANEOUSPOSITIONOFATOM!XX(J,I)：X*X,J:第几天实验,I:第几步跳跃!XXM(I):THEMEANOFXXWRITE(*,*)实验天数Jt，实验次数IcREAD(*,*)Jt,IcISEED=RTC()DOJ=1,Jt!第几天实验X=0!每天都是从原点出发DOI=1,Ic!第几步跳跃RN=RAN(ISEED)IF(RN0.5)THENX=X+1ELSEX=X-1ENDIFXX(J,I)=X*X!记录下原子每天每次跳动后离原点的距离ENDDOENDDOOPEN(1,FILE=“f:\DIF1.DAT)DOI=1,IcXXM=0.0XXM(I)=1.0*SUM(XX(1:Jt,I))/Jt!!WRITE(1,*)I,XXM(I)ENDDOCLOSE(1)ENDrnRn2随机行走：原子扩散的平均距离与原子跳动次数的平方根成正比。变量的定义：X——记录原子的位置xx(j,i)——记录原子每次跳动后距离原点的距离j：统计的天数，i：跳动的次数xxm(i)——统计距离与次数的关系实验天数：Jt，实验次数：Ic二维随机行走随机的确定方法一、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=yENDIF方法二XN=(/0,0,-1,1/)YN=(/-1,1,0,0/)RN=4*RAN(ISEED)+1X=X+YN(RN)Y=Y+YN(RN)!MonteCarloSimulationofTwoDimensionalDiffusionINTEGERX,Y,XY(1:1000,1:1000)REALXYM(1:1000)WRITE(*,*)实验天数Jt，实验次数IcREAD(*,*)Jt,IcISEED=RTC()DOJ=1,Jt!第几天实验X=0!!!Y=0!!!DOI=1,Ic!第几步跳跃RN=RAN(ISEED)IF(RN.LT.0.25)THENx=xy=y-1elseIF(RN.LT.0.5.AND.RN.GE.0.25)THENx=xy=y+1elseIF(RN.LT.0.75.AND.RN.GE.0.5)THENx=x-1y=yelseIF(RN.GE.0.75)THENx=x+1y=yENDIFXY(J,I)=X*X+Y*YENDDOENDDOOPEN(1,FILE=f:\DIF2.DAT)DOI=1,IcXYM=0.0XYM(I)=1.0*SUM(XY(1:Jt,I))/Jt!!WRITE(1,*)I,XYM(I)ENDDOCLOSE(1)END!MonteCarloSimulationofTwoDimensionalDiffusionINTEGERX,XY(1:1000,1:1000),y,XN(1:4),YN(1:4),RNREALXYM(1:1000)!X:INSTANTANEOUSPOSITIONOFATOM!XY(J,I)：X*Y,J:第几天实验,I:第几步跳跃!XYM(I):THEMEANOFXYWRITE(*,*)实验天数Jt，实验次数IcREAD(*,*)Jt,IcXN=(/0,0,-1,1/)YN=(/-1,1,0,0/)ISEED=RTC()DOJ=1,Jt!第几天实验X=0!!!Y=0!!!DOI=1,Ic!第几步跳跃RN=4*RAN(ISEED)+1X=X+YN(RN)Y=Y+YN(RN)XY(J,I)=X*X+Y*YENDDOENDDOOPEN(1,FILE=C:\DIF2.DAT)DOI=1,IcXYM=0.0XYM(I)=1.0*SUM(XY(1:Jt,I))/Jt!!WRITE(1,*)I,XYM(I)ENDDOCLOSE(1)！做三维空间随机行走？ENDSimulationofRecrystallisation&GrainGrowthbyMeansofaMonteCarloModelModelMicrostructureThegrainstructureismappedontoadiscretetwoorthreedimensionallattice.ToeachlatticesiteanintegerSbetween1andamaximumvalueQisassigned.Theseintegersrepresentthecrystallographicorientationofthedifferentgrainsandarecalledspinsororientations.Thusbetweentwoadjacentlatticesitesofunlikespinsthereisagrainboundarysegment,whiletwoneighboringsitesoflikespinsbelongtothesamegrainandthereisnograinboundarybetweenthem.2Algorithm(proposedbyAndersonetal.1984)1.PickupalatticesiteiwithspinS_iatrandom2.CalculateitsenergyE_iaccordingtoequation(1)3.PickupanewspinS_jdifferentfromtheoldone4.CalculatetheenergyofthechosensiteiifithadthisnewspinS_j5.FlipsiteitospinS_jwithprobabilityPaccordingtoequation(2)6.Repeatsteps1.-5.EnergyE_iofalatticesite:(1)J:EnergyfactorNN:Numberofnearestneighbours:Kronecker-deltaFlipProbabilityP:(2):Energydifferenceduetotheflip()Kb:BoltzmannkT:SystemtemperatureTimeUnit:Onetimeunit(1MonteCarlostep:MCS)equalsthenumberoflatticesitesofthesystem2.蒙特卡罗模拟双晶粒的长大如何表示晶粒：颜色（状态）相同的像素（元胞）；如何定义一段晶界：不同状态元胞的交线；如何表示一段晶界：形成晶界的两个元胞；如何定义邻居：最近邻，次近邻；画一个圆晶粒USEMSFLIBINTEGERXR，YR!在的区域中画一个圆PARAMETERXR=400,YR=400INTEGERR,S(1:XR,1:YR)X0=XR/2!圆心位置X0，YOY0=YR/2R=MIN(X0-10,Y0-10)!圆半径S=0!像素的初始状态（颜色）DOI=1,XRDOJ=1,YRIF((I-X0)**2+(J-Y0)**2=R**2)S(I,J)=10IER=SETCOLOR(S(I,J))IER=SETPIXEL(I,J)ENDDOENDDOEND画晶界XN=(/-1,0,1,-1,1,-1,0,1/)YN=(/-1,0,-1,0,0,1,1,1/)DOI=1,XR!画晶界DOJ=1,YRNDS=0DOK=1,8IF(S(I,J).NE.S(I+XN(K),J+YN(K)))NDS=NDS+1ENDDOIF(NDS0)THENIER=SETCOLOR(9)ELSEIER=SETCOLOR(S(I,J))ENDIFIER=SETPIXEL(I,J)ENDDOENDDO转变规则——能量最小原理如何实现界面的迁移：单胞状态的转变转变规则：随机选择一个邻居，如果转变后系统能量降低（考虑能量起伏）。如何计算能量：体积能（动力）；界面能（阻力）：体积能EV计算：EV=0EV(0)=ER界面能如何计算一个单胞界面能：异类邻居数之和。XN=(/-1,-1,-1,0,0,1,1,1/)YN=(/-1,0,1,-1,1,-1,0,1/)DOII=1,8ISB(II)=IS(I+XN(II),J+YN(II))ENDDOE0=COUNT(ISB.NE.IS(I,j))任意选邻居再计算能量随机选取一个邻居CELL及能量变化ISTR=ISB(8*RAN(ISEED)+1)IF(ISTR.EQ.IS)CYCLEETR=COUNT(ISB.NE.ISTR)能量判断单个元胞的体积能Ev与元胞的一个面的能量Es之比：Ev=8Es能量变化与能量起伏DEB=ETR-E0DEV=EV(ISTR)-EV(IS)DE=DEB+DEV+2.5*RAN(ISEED))-1.25MC模拟双晶粒的长大如何定义基体和晶粒；如何寻找边界；如何计算能量（构造或描述问题的概率过程）步骤：1、读取晶粒生长步长；2、在基体上画一个原始晶粒，并赋予基体和原始晶粒不同的状态（取向号或体积能）3、寻找晶界4、计算晶界上网格点的相互作用，每个格点的相互作用导致晶粒长大或缩小：①随机选取一个初始格点；②若此点属于晶界，那么可以随机转变为它邻居的状态，若不是晶界，则跳出循环，不发生转变；③计算转变前后的能量变化，⊿E=⊿Ev+⊿Es+⊿Eq(自由能=体积能+表面能+能量起伏)④若⊿E小于0，则新晶相被接受，网格状态发生转变。双晶长大/收缩①USEMSFLIB②PARAMETERIR=400,JR=400③INTEGERIS(0:IR+1,0:JR+1),TMAX,ISB(1:8),ISTR,T,NR,IX,IY,X,Y,RD④INTEGERXN(1:8),YN(1:8)⑤XN=(/-1,0,1,-1,1,-1,0,1/)⑥YN=(/-1,0,-1,0,0,1,1,1/)⑦WRITE(*,*)PLEASEINPUTTHETIMESTEP⑧READ(*,*)TMAX⑨ISEED=RTC()⑩!定义圆心和半径11X0=IR/212Y0=JR/213R=MIN(X0,Y0)-10014!定义基体和圆晶粒分别为状态5、状态215IS=516DOI=1,IR17DOJ=1,JR18IF((I-X0)**2+(J-Y0)**2=R**2)IS(I,J)=219IER=SETCOLOR(IS(I,J))20IER=SETPIXEL(I,J)21ENDDO22ENDDO23OPEN(1,FILE=E:\data.DAT)24!寻找晶粒边界，计算能量，改变状态。25DOT=1,TMAX26DOX=1,IR27DOY=1,JR28IX=IR*RAN(ISEED)+129JY=JR*RAN(ISEED)+130DOII=1,831ISB(II)=IS(IX+XN(II),JY+YN(II))32ENDDO33E0=COUNT(ISB.NE.IS(IX,JY))34!如果不是圆晶粒边界，则跳出，重新循环35IF(E0.EQ.0)CYCLE36!随机寻找一个相邻点37NR=8*RAN(ISEED)+138IS