Appendix: MATLAB Code
%Nanolobes Project%
clear all
clc
clf
%Define constant parameters
g=9e8;%bulk free energy%
gamma=0.0447;%constant surface energy
L=5e-17;%interface mobility%
M=1e-32;%atoms mobility%
E=5e6;%elastic modulus
sigma0=5e7;%thermal stress
r=3e6;%proportionality constant
%define initial cosine wave surface
y0=500e-9;%initial film thickness
dX=1000e-9;%area of interest
dY=200e-9;%amplitude of waves
lambda=500e-9;%wavelength of perturbations
n=50;%starting number of nodes
for i=1:n
X(i)=dX*i/n;
Y(i)=y0-dY*cos(2*pi*X(i)/lambda);
end
subplot (3,1,1), scatter(X,Y)
%decide on number and size of time steps
dt=0.002; %seconds
T=1000;%number of steps
%calculate initial element sizes
l=zeros(n-1,1);
for i=1:(n-1)
l(i)=((X(i+1)-X(i))^2+(Y(i+1)-Y(i))^2)^0.5;
end
%keep track of number of elements with each iteration
N=zeros(T,1);
for s=1:T
% recalculate element sizes
n=length(X);
l=zeros(n-1,1);
for i=1:(n-1)
l(i)=((X(i+1)-X(i))^2+(Y(i+1)-Y(i))^2)^0.5;
end
% define upper and lower bounds for element size
lmin=mean(l)/10;
lmax=3*mean(l);
% check if any elements are too large or too small
Xc=zeros(n-1);
Yc=zeros(n-1);
p=0; %number of nodes to be added
q=0; %number of nodes to be deleted
for i=1:(n-1)
% if too large, calculate position of new node at midpoint
if l(i)>=lmax
Xc(i)=0.5*(X(i)+X(i+1));
Yc(i)=0.5*(Y(i)+Y(i+1));
p=p+1;
end
% if too small, mark node to be deleted
if l(i)<=lmin
Xc(i+1)=-1;
Yc(i+1)=-1;
q=q+1;
end
end
% adjust arrays and input new nodal positions
XX=zeros(n+p-q,1);
YY=zeros(n+p-q,1);
k=1;
for i=1:n
XX(k)=X(i);
YY(k)=Y(i);
% add new nodes at midpoints
if Xc(i)>0
k=k+1;
XX(k)=Xc(i);
YY(k)=Yc(i);
end
% remove nodes for undersized elements
if Xc(i)==-1;
k=k-1;
end
k=k+1;
end
X=XX;
Y=YY;
% recalculate number and size of elements
n=length(X);
N(s,1)=n;
l=zeros(n,1);
for a=1:(n-1)
l(a)=((X(a+1)-X(a))^2+(Y(a+1)-Y(a))^2)^0.5;
end
m=4*n-1; %size of global H and f matrices
%Defining theta matrix%
sintheta=zeros(n-1,1);
costheta=zeros(n-1,1);
for i=1:(n-1)
sintheta(i)=(Y(i+1)-Y(i))/l(i);
costheta(i)=(X(i+1)-X(i))/l(i);
end
%Define f matrix (constant gamma, with thermal stress)
f=zeros(m,1);
Yavg=mean(Y);
for i=1:(n-1)
%Yavg=0.5*(Y(i)+Y(i+1));
W1=((sigma0*(1-r*(Y(i)-Yavg)))^2)/E;
W2=((sigma0*(1-r*(Y(i+1)-Yavg)))^2)/E;
A=(2*W1+W2)/3;
B=(W1+2*W2)/3;
Q=[costheta(i); sintheta(i); 0; 0; -costheta(i); -sintheta(i); 0];
R=[-sintheta(i); costheta(i); 0; 0; -sintheta(i); costheta(i); 0];
S=[A*sintheta(i);-A*costheta(i);0;0;B*sintheta(i);-B*costheta(i); 0];
%reset and redefine local f matrix for element i
F=zeros(7,1);
F=gamma*Q+0.5*l(i)*(g*R+S);
%add local f matrix elements into global f matrix
k=4*(i-1);
for a=1:7
f(k+a)=f(k+a)+F(a);
end
end
%Define H matrix%
H=zeros(m,m,1);
mu=zeros(n-1);
for i=1:(n-1)
%reset and redefine local H matrix for element i
h=zeros(7,7);
mu(i)=(L*l(i)^2)/(5*M);
%diagonal elements
h(1,1)=2*(l(i)*sintheta(i))^2;
h(5,5)=h(1,1);
h(2,2)=2*(l(i)*costheta(i))^2;
h(6,6)=h(2,2);
h(3,3)=14+4*mu(i);
h(7,7)=h(3,3);
%node i elements
h(1,2)=-2*(l(i)^2)*sintheta(i)*costheta(i);
h(2,1)=h(1,2);
h(1,3)=5*l(i)*sintheta(i);
h(3,1)=h(1,3);
h(2,3)=-5*l(i)*costheta(i);
h(3,2)=h(2,3);
%node i+1 elements
h(5,6)=-2*(l(i)^2)*sintheta(i)*costheta(i);
h(6,5)=h(5,6);
h(5,7)=-5*l(i)*sintheta(i);
h(7,5)=h(5,7);
h(6,7)=5*l(i)*costheta(i);
h(7,6)=h(6,7);
%node i and node i+1 cross elements
h(1,5)=(l(i)*sintheta(i))^2;
h(5,1)=h(1,5);
h(1,6)=-(l(i)^2)*sintheta(i)*costheta(i);
h(6,1)=h(1,6);
h(1,7)=-l(i)*sintheta(i);
h(7,1)=h(1,7);
h(2,5)=-(l(i)^2)*sintheta(i)*costheta(i);
h(5,2)=h(2,5);
h(2,6)=(l(i)*costheta(i))^2;
h(6,2)=h(2,6);
h(2,7)=l(i)*costheta(i);
h(7,2)=h(2,7);
h(3,5)=l(i)*sintheta(i);
h(5,3)=h(3,5);
h(3,6)=-l(i)*costheta(i);
h(6,3)=h(3,6);
h(3,7)=2-mu(i);
h(7,3)=h(3,7);
%jm elements
h(1,4)=-4*l(i)*sintheta(i);
h(4,1)=h(1,4);
h(2,4)=4*l(i)*costheta(i);
h(4,2)=h(2,4);
h(3,4)=-16+2*mu(i);
h(4,3)=h(3,4);
h(4,5)=4*l(i)*sintheta(i);
h(5,4)=h(4,5);
h(4,6)=-4*l(i)*costheta(i);
h(6,4)=h(4,6);
h(4,7)=-16+2*mu(i);
h(7,4)=h(4,7);
h(4,4)=32+16*mu(i);
h=h/(6*L*l(i));
%add local h matrix elements into global H matrix
k=4*(i-1);
for a=1:7
for b=1:7
H(k+a,k+b)=H(k+a,k+b)+h(a,b);
end
end
end
%force first and last node to be stationary along X-axis
H(1,1)=-1e90;
H(3,3)=-1e90;
H(m-2,m-2)=-1e90;
H(m-1,m-1)=-1e90;
%Calculate nodal velocities
V=H\f;
%Extract nodal X,Y velocities and prevent runaway velocities
VX=zeros(n,1);
VY=zeros(n,1);
Vmax=1e-6;
Vc=5e-8;
for i=1:n
k=4*(i-1);
VX(i)=V(k+1);
VY(i)=V(k+2);
if abs(VX(i))>=Vmax
VX(i)=VX(i)*Vc;
end
if abs(VY(i))>=Vmax
VY(i)=VY(i)*Vc;
end
end
%calculate new nodal positions
for i=1:n
X(i)=X(i)+dt*VX(i);
Y(i)=Y(i)+dt*VY(i);
end
%delete anomolous nodes
for i=n:-1:1
% if nodal position exceeds bounds, delete node
if X(i)<0 || X(i)>dX || Y(i)<0
X(i)=[];
Y(i)=[];
end
end
n=length(X);
subplot (3,1,2), plot(X,Y)
hold all
disp(s);
end