MATLAB code: Explicit finite difference method for phase-field dendritic growth model
Date:
%Main code
time0 = clock();
format long;
Nx = 300;
Ny = 300;
NxNy = Nx * Ny;
dx = 0.03;
dy = 0.03;
nstep = 3500;
nprint = 50;
dtime = 1.0e-4;
tau = 0.0003;
epsilonb = 0.01;
kappa = 1.8;
delta = 0.02;
aniso = 6.0;
alpha = 0.9;
gamma = 10.0;
teq = 1.0;
theta0 = 0.3;
seed = 10.0;
phi = zeros(Nx, Ny);
tempr = zeros(Nx, Ny);
for i = 1:Nx for j = 1:Ny
if ((i - 150)^2 + (j - 150)^2 < seed )
phi(i, j) = 1.0;
end
end end
X = 0:0.03:9 - 0.03;
Y = 0:0.03:9 - 0.03;
Z = zeros(300, 300);
for istep = 1:nstep
%----
% compute the laplacians and epsilon:
%---
for i = 1:Nx
for j = 1:Ny
jp = j + 1;
jm = j - 1;
ip = i + 1;
im = i - 1;
if (im == 0)
im = Nx;
end
if (ip == (Nx + 1))
ip = 1;
end
if (jm == 0)
jm = Ny;
end
if (jp == (Ny + 1))
jp = 1;
end
hne = phi(ip, j);
hnw = phi(im, j);
hns = phi(i, jm);
hnn = phi(i, jp);
hnc = phi(i, j);
lap_phi(i, j) = (hnw + hne + hns + hnn -4.0 * hnc) / (dx * dy);
hne = tempr(ip, j);
hnw = tempr(im, j);
hns = tempr(i, jm);
hnn = tempr(i, jp);
hnc = tempr(i, j);
lap_tempr(i, j) = (hnw + hne + hns + hnn -4.0 * hnc) / (dx * dy);
phidx(i, j) = (phi(ip, j) - phi(im, j)) / (2.0 * dx);
phidy(i, j) = (phi(i, jp) - phi(i, jm)) / (2.0 * dy);
theta = atan2(phidy(i, j), phidx(i, j));
epsilon(i, j) = epsilonb * (1.0 + delta * cos(aniso * (theta - theta0)));
epsilon_deriv(i, j) = -epsilonb * aniso * delta * sin(aniso * (theta - theta0));
end %for j
end %for i
for i = 1:Nx
for j = 1:Ny
jp = j + 1;
jm = j - 1;
ip = i + 1;
im = i - 1;
if (im == 0)
im = Nx;
end
if (ip == (Nx + 1))
ip = 1;
end
if (jm == 0)
jm = Ny;
end
if (jp == (Ny + 1))
jp = 1;
end
phiold = phi(i, j);
%-- first term:
term1 = (epsilon(i, jp) * epsilon_deriv(i, jp) * phidx(i, jp) - ...
epsilon(i, jm) * epsilon_deriv(i, jm) * phidx(i, jm)) / (2.0 * dy);
%-- second term:
term2 = -(epsilon(ip, j) * epsilon_deriv(ip, j) * phidy(ip, j) - ...
epsilon(im, j) * epsilon_deriv(im, j) * phidy(im, j)) / (2.0 * dx);
%-- factor m:
m = alpha / pi * atan(gamma * (teq - tempr(i, j)));
%-- Time integration:
phi(i, j) = phi(i, j) +(dtime / tau) * (term1 + term2 + epsilon(i, j)^2 * lap_phi(i, j) + ...
phiold * (1.0 - phiold) * (phiold -0.5 + m));
if (phi(i, j) < 10^ - 30)
phi(i, j) = 10^ - 30;
end
if (phi(i, j) > 10^30)
phi(i, j) = 10^30;
end
Z(i, j) = phi(i, j);
%-- evolve temperature:
tempr(i, j) = tempr(i, j) +dtime * lap_tempr(i, j) + kappa * (phi(i, j) - phiold);
if (tempr(i, j) < 10^ - 30)
tempr(i, j) = 10^ - 30;
end
if (tempr(i, j) > 10^30)
tempr(i, j) = 10^30;
end
end
end
if (mod(istep, nprint) == 0)
figure(1);
surf(X,Y,Z);
shading interp
colormap jet;
axis image;
view(0,90);
M(istep) = getframe;
end %if
end %istep
compute_time = etime(clock(), time0);
function [phi,tempr] = nucleus(Nx,Ny,seed)
format long;
phi = zeros(Nx,Ny);
tempr = zeros(Nx,Ny);
for i=1:Nx
for j=1:Ny
if ((i-Nx/2)^2+(j-Ny/2)^2 < seed)
phi(i,j) = 1.0;
end
end
end
end %endfunction
function [phi,tempr] = nucleus(Nx,Ny,seed)
format long;
phi = zeros(Nx,Ny);
tempr = zeros(Nx,Ny);
for i=1:Nx
for j=1:Ny
if ((i-Nx/8)^2+(j-Ny/2)^2 < seed)
phi(i,j) = 1.0;
end
if ((i-2*Nx/8)^2+(j-Ny/2)^2 < seed)
phi(i,j) = 1.0;
end
if ((i-3*Nx/8)^2+(j-Ny/2)^2 < seed)
phi(i,j) = 1.0;
end
if ((i-4*Nx/8)^2+(j-Ny/2)^2 < seed)
phi(i,j) = 1.0;
end
if ((i-5*Nx/8)^2+(j-Ny/2)^2 < seed)
phi(i,j) = 1.0;
end
if ((i-6*Nx/8)^2+(j-Ny/2)^2 < seed)
phi(i,j) = 1.0;
end
if ((i-7*Nx/8)^2+(j-Ny/2)^2 < seed)
phi(i,j) = 1.0;
end
if ((i-8*Nx/8)^2+(j-Ny/2)^2 < seed)
phi(i,j) = 1.0;
end
end
end
end %endfunctio