sunouhebianbianjie.zip_Dendrite growth_phase field_相场 枝晶_相场模拟 枝晶

%模拟镍的相场计算(无任何耦合、粗网格、变边界条件)----------------------------- clear;clc;clf;tic %参数---------------------------------------------------------------------- m=0.05; eps=0.005; r=0.08; %动力学系数、界面相关参数、各向异性系数 Tm=1728; L=2350; DT=0.155; Cp=5.42; %熔点、潜热、热扩散率、比热 Delta=-0.5; R=0.05; alfa=400; %无量纲过冷度、形核半径、耦合相关参数alfa EPS=eps^2;A=m/EPS;B=(m/EPS)*EPS;C=eps*alfa*(-Delta);D=30/(-Delta); %网格划分------------------------------------------------------------------- N=50; NTimeSteps=10000; Dx=0.02; Dy=Dx; Dt=1e-5; E=(m*Dt)/eps^2; DXX=Dx^2; [x,y]=meshgrid(0:Dx:1); t=100; 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); 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 theta_x=(phy_x.*phy_xy - phy_y.*phy_xx)./(phy_x.^2+phy_y.^2); theta_y=(phy_x.*phy_yy - phy_y.*phy_xy)./(phy_x.^2+phy_y.^2); sin4theta=4.*(phy_x.^3*phy_y-phy_y.^3*phy_x)./(phy_x.^2+phy_y.^2)^2; cos4theta=1-8.*phy_x.^2.*phy_y.^2./(phy_x.^2+phy_y.^2)^2; rcos4th=1+r*cos4theta; r1cos4th=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*(r1cos4th*rcos4th-r2sin4th2).*theta_x.*phy_y; A4=-EPS*(r1cos4th*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 if n==t t=t+100; 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 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 T=T+1 end toc