【物理应用】基于相场法模拟金属镍的晶粒的长大过程附matlab完整代码

1 内容介绍

金属镍是支撑信息技术产业和光伏产业发展的基础材料,它的性能直接影响和决定着相关高技术产业的发展.随着直拉硅单晶的直径不断增大以及集成电路特性线宽持续减小,金属镍的缺陷控制及其质量保证尤为重要,特别是新型的基于快速热处理(RTP)的”内吸杂”结构在热处理过程中硅片内部演变行为的研究.近些年来,计算机数值模拟技术越来越多地被用于材料微缺陷演变的研究中.其中,以吉布斯-朗道自由能理论为基础的相场法,已经成为研究材料介-微观缺陷的一种重要数值模拟方法.

2 完整代码

%模拟镍的相场计算(无任何耦合、粗网格、变边界条件)—————————–
clear;clc;clf;tic
close all
%参数———————————————————————-
m=0.05; eps=0.005;   r=0.05;        %动力学系数、界面相关参数、各向异性系数
Tm=1728; L=2350; DT=0.155; Cp=5.42;  %熔点、潜热、热扩散率、比热
Delta=-0.6; R=0.05; alfa=400;        %无量纲过冷度、形核半径、耦合相关参数alfa
EPS=eps^2;A=m/EPS;B=(m/EPS)*EPS;C=eps*alfa*(-Delta);D=30/(-Delta);
%网格划分——————————————————————-
N=100;  NTimeSteps=5000;   Dx=0.002; Dy=Dx;  Dt=1e-5; E=(m*Dt)/eps^2; DXX=Dx^2;                 
[x,y]=meshgrid(0:Dx:2);     t=20; T=0; h=Dt/Dx/Dx;
%初始条件－—————————————————————–
nn = length(x);
for i=1:nn
    for j=1:nn
        %phy_n(i,j)=0.5*(tanh((sqrt(x(i,j).^2+y(i,j)^2)-R)/(2*sqrt(2)*eps))+1);
        phy_n(i,j)=0.5*(tanh((sqrt((x(i,j)-x(N/2,N/2)).^2+(y(i,j)-y(N/2,N/2))^2)-R)/(2*sqrt(2)*eps))+1);
            U_n(i,j)=Delta* phy_n(i,j);
        end
end
%主程序———————————————————————
for n=1:NTimeSteps  
    %边界条件－—————————————————————–
    phy_n(:,1)=phy_n(:,2);  phy_n(:,N+1)=phy_n(:,N);
    U_n(:,1)=U_n(:,2);      U_n(:,N+1)=U_n(:,N);
    phy_n(1,:)=phy_n(2,:);  phy_n(N+1,:)=phy_n(N,:);
    U_n(1,:)=U_n(2,:);      U_n(N+1,:)=U_n(N,:);
%————————————————————————–
 for i=2:nn-1
        for j=2:nn-1
            phy_x=(phy_n(i+1,j)-phy_n(i-1,j))./2./Dx;
            phy_y=(phy_n(i,j+1)-phy_n(i,j-1))./2./Dy;
            phy_xy=(phy_n(i+1,j+1)-phy_n(i+1,j-1)-phy_n(i-1,j+1)+phy_n(i-1,j-1))./4./Dy.^2;
            phy_xx=(phy_n(i+1,j)-2*phy_n(i,j)+phy_n(i-1,j))./Dx.^2;
            phy_yy=(phy_n(i,j+1)-2*phy_n(i,j)+phy_n(i,j-1))./Dy.^2;
             
            D_phy=phy_n(i,j-1)+phy_n(i,j+1)+phy_n(i+1,j)+…
                  phy_n(i-1,j)+0.5.*(phy_n(i+1,j-1)+phy_n(i-1,j-1)+phy_n(i+1,j+1)+phy_n(i-1,j+1))-6.*phy_n(i,j);
            D_U=U_n(i,j-1)+U_n(i,j+1)+U_n(i+1,j)+U_n(i-1,j)-4.*U_n(i,j);
            D_U_PHY=phy_n(i,j).*(1-phy_n(i,j)).*(phy_n(i,j)-0.5+30*C*U_n(i,j).*phy_n(i,j).*(1-phy_n(i,j)));
             
            if  phy_x==0 & phy_y==0
                 phy_n1(i,j)=phy_n(i,j)+Dt*A*D_U_PHY+B*h.*D_phy;
                 U_n1(i,j)=U_n(i,j)+h*D_U-D*phy_n(i,j).^2.*(1-phy_n(i,j)).^2.*(phy_n1(i,j)-phy_n(i,j));             
            else
              
            F=phy_x.^2+phy_y.^2;
            theta_x=(phy_x.*phy_xy – phy_y.*phy_xx)./F;
            theta_y=(phy_x.*phy_yy – phy_y.*phy_xy)./F;
            sin4theta=4.*(phy_x.^3*phy_y-phy_y.^3*phy_x)./F^2;
            cos4theta=1-8.*phy_x.^2.*phy_y.^2./F^2;
            rcos4th=1+r*cos4theta;    r16cos4th=16.*r*cos4theta;  r2sin4th2=16.*r^2.*sin4theta.^2;
                         
            A1=EPS/DXX*rcos4th.^2.*D_phy;
            A2=-EPS*8.*r*sin4theta*rcos4th*(theta_x.*phy_x+theta_y.*phy_y);
            A3=EPS*(r16cos4th*rcos4th-r2sin4th2).*theta_x.*phy_y;
            A4=-EPS*(r16cos4th*rcos4th-r2sin4th2).*theta_y.*phy_x;
                                            
           phy_n1(i,j)=phy_n(i,j)+E*(D_U_PHY+A1+A2+A3+A4);                
           U_n1(i,j)=U_n(i,j)+h*D_U-D*phy_n(i,j).^2.*(1-phy_n(i,j)).^2.*(phy_n1(i,j)-phy_n(i,j));  
            end
        end
    end  
%更新   
 for i=2:nn-1
        for j=2:nn-1
       phy_n(i,j)=phy_n1(i,j); U_n(i,j)=U_n1(i,j);
        end
 end
%画图
 if n==t
        t=t+50;
        figure(1);meshc(x,y,phy_n);colormap(cool);
        hidden off;
        shading flat;
        colorbar;view(0,90);
        figure(2);meshc(x,y,U_n);colormap(cool);
        hidden off;
        shading flat;
        colorbar;view(0,90);
    end
 T=T+1
end
toc

3 运行结果

4 参考文献

博主简介：擅长智能优化算法、神经网络预测、信号处理、元胞自动机、图像处理、路径规划、无人机、雷达通信、无线传感器等多种领域的Matlab仿真，相关matlab代码问题可私信交流。

部分理论引用网络文献，若有侵权联系博主删除。

来源：Matlab科研工作室