%% Project5.m
%%
%% List of nodes with coordinates (Node #, x-coord,
y-coord, z-coord)%
%node = [1 ? ? ?;
% 2 ? ? ?;
%
3 ? ? ?]; % etc...
%% List of elements with its nodes (Element #,
Node #1, Node #2, Node #3)
%element = [1 ? ? ?;
% 2
? ? ?;
%
3 ? ? ?]; % etc...
%% List of nodes at grain boundaries (Node #,
Node # ..etc)
%grain1node = [? ? ? ?]; % etc... Must not include
centergrain node or clamped node
%grain2node = [? ? ? ?]; % etc...
%grain3node = [? ? ? ?]; % etc...
%centergrainnode = ?;
%
%% List of nodes at clamped boundaries (Node
#, Node #, Node # ...etc)
%clampednode = [? ? ? ? ? ? ? ? ?]; % etc...
Perimeter nodes
clear
loadnodes2; % All
of the above is done by executing the file 'loadnodes2.m'
totalnodes = length(node); %
# of nodes ('nodes' must contain more than 3 nodes)
ELnum = length(element); %
# of elements ('element' must contain more than 3 elements)
%Hglobal = zeros(totalnodes*3); % Create global
mobility matrix
%Fglobal = zeros(totalnodes*3,1); % Create global
force vector
m = .001; % Define
mobility
Gamma = 1.2E-10; %
Define inplane force constant (mobility)
GammaB = 1.5E-10; %
Define grain force term
GammaC = 2.0E-10; %
Define triple point force term
DeltaT = 3; % Time
interval per iteration
Iteration = 10000; %
Iterations for each simulation
for it = 1:Iteration, % Start the whole thing
Hglobal = zeros(totalnodes*3); %
Create global mobility matrix
Fglobal = zeros(totalnodes*3,1); %
Create global force vector
for EL = 1:ELnum, %
element has been selected element #dont use i,j or k!
l1 = 0; l2 = 0; l3 = 0; %
Set local nodes to zero
l1node = 0; l2node = 0; l3node = 0; %
Reset the global node numbers of the local nodes
l1node = element(EL,2); %
Connect local node # to global node #
l2node = element(EL,3);
l3node = element(EL,4);
l1 = [node(l1node,2); %
Set local node 1 to global coordinates
node(l1node,3);
node(l1node,4)];
l2 = [node(l2node,2); %
Same for node 2
node(l2node,3);
node(l2node,4)];
l3 = [node(l3node,2); %
And for node 3
node(l3node,3);
node(l3node,4)];
v1 = l2-l1; % Get
two vectors pointing away from local node 1 towards local node 2 and 3
v2 = l3-l1;
n = 0; %reset normal
vector
n = cross(v1,v2); %
Get normal vector to element triangle
if n(3) < 0 %
Check that z component of n is positive since Nz negative cannot happen
n = -n;
end
nx = n(1);
ny = n(2);
nz = n(3);
phi = 0; % Reset
angles
theta = 0;
if nx >= 0 & ny > 0 %
Determines the angle Phi based on the signs of Nx and Ny
phi = -atan(nx/ny);
end
if nx < 0 & ny > 0
phi = -atan(nx/ny);
end
if nx >= 0 & ny < 0
phi = pi - atan(nx/ny);
end
if nx < 0 & ny < 0
phi = pi - atan(nx/ny);
end
if nx == 0 & ny == 0
phi = 0;
end
if nx > 0 & ny == 0
phi = 3*pi/2;
end
if nx < 0 & ny == 0
phi = pi/2;
end
nlength = sqrt(nx^2 + ny^2 + nz^2);
theta = asin(nz/nlength)-pi/2; %MAY BE WRONG!!! Ok, due to order of rotation.
T1 = 0; % Reset rotation
matrices
T2 = 0;
T1 = [cos(phi) sin(phi) 0; -sin(phi) cos(phi)
0; 0 0 1];
T2 = [1 0 0; 0 cos(theta) sin(theta); 0 -sin(theta)
cos(theta)];
l1p = T2*T1*l1; %Transform
global nodes into local, same z-plane coordinates)
l2p = T2*T1*l2;
l3p = T2*T1*l3;
% Determine the determinant of the Jacobian
DetJ = (l1p(1)-l3p(1))*(l2p(2)-l3p(2))-(l2p(1)-l3p(1))*(l1p(2)-l3p(2));
% Determine the local mobility matrix H
T1m = [cos(-phi) sin(-phi) 0; -sin(-phi) cos(-phi)
0; 0 0 1]; % Backward rot matrices
T2m = [1 0 0; 0 cos(-theta) sin(-theta); 0 -sin(-theta)
cos(-theta)]; % Minus angle
t1t2m = T1m*T2m;
t2t1 = T2*T1;
for j = 1:3, % Create
giant Backward and Forward matrices
for k= 1:3,
T1T2m(j,k) = t1t2m(j,k);
T1T2m(j+3,k+3) = t1t2m(j,k);
T1T2m(j+6,k+6) = t1t2m(j,k);
T2T1(j,k) = t2t1(j,k);
T2T1(j+3,k+3) = t2t1(j,k);
T2T1(j+6,k+6)
= t2t1(j,k);
end
end
NTN = [1/12 1/12 0 1/24 1/24 0 1/24 1/24 0; %
Integral calculated using calculator
1/12 1/12 0 1/24 1/24 0 1/24
1/24 0;
0 0 0 0 0 0 0
0 0;
1/24 1/24 0 1/12
1/12 0 1/24 1/24 0;
1/24 1/24 0 1/12
1/12 0 1/24 1/24 0;
0 0 0 0 0 0 0
0 0;
1/24 1/24 0 1/24
1/24 0 1/12 1/12 0;
1/24 1/24 0 1/24
1/24 0 1/12 1/12 0;
0 0 0 0 0 0 0
0 0];
H = abs(DetJ/m)*T1T2m*NTN*T2T1;
for i = 1:9, % Input
local mobility matrix H into global mobility matrix Hglobal
for j = 1:9,
if i == 1 | i
== 2 | i == 3
nop1 = (l1node-1)*3 + i;
end
if j == 1 | j
== 2 | j == 3
nop2 = (l1node-1)*3 + j;
end
if i == 4 | i
== 5 | i == 6
nop1 = (l2node-1)*3 + (i-3);
end
if j == 4 | j
== 5 | j == 6
nop2 = (l2node-1)*3 + (j-3);
end
if i == 7 | i
== 8 | i == 9
nop1 = (l3node-1)*3 + (i-6);
end
if j == 7 | j
== 8 | j == 9
nop2 = (l3node-1)*3 + (j-6);
end
Hglobal(nop1,nop2) = Hglobal(nop1,nop2) + H(i,j); %
Assemble global mobility matrix
end
end
nodCO = [l1 l2 l3];
ite = 1;
for i = 1:3,
for j = 1:3,
if i ~= j
nftemp = [nodCO(1,j)-nodCO(1,i);
nodCO(2,j)-
nodCO(2,i);
nodCO(3,j)-nodCO(3,i)]/sqrt((nodCO(1,j)-nodCO(1,i))^2+(nodCO(2,j)- nodCO(2,i))^2+(nodCO(3,j)-nodCO(3,i))^2);
nf(1,ite)
= nftemp(1);
nf(2,ite)
= nftemp(2);
nf(3,ite)
= nftemp(3);
end
ite = ite + 1;
end
end
for k = 1:3,
n12(k,1) = nf(k,2);
n13(k,1) = nf(k,3);
n21(k,1) = nf(k,4);
n23(k,1) = nf(k,6);
n31(k,1) = nf(k,7);
n32(k,1) = nf(k,8);
end
F = Gamma*[n12 + n13; n21 + n23; n31 + n32]; % Create local force matrix
for i = 1:9, % Input
local force matrix F into global force matrix Fglobal
if i == 1 | i
== 2 | i == 3
nof = (l1node-1)*3 + i;
end
if i == 4 | i
== 5 | i == 6
nof = (l2node-1)*3 + (i-3);
end
if i == 7 | i
== 8 | i == 9
nof = (l3node-1)*3 + (i-6);
end
if i == 3 | i
== 6 | i == 9
p=10000;
end
if i == 1 | i
== 2 | i == 4 | i == 5 | i == 7 | i == 8
p=1;
end
Fglobal(nof) =
Fglobal(nof) + F(i)*p; % Assemble global force
matrix
end
allgrains = [grain1node grain2node grain3node];
% Node #'s of all grain nodes
grainlength = length(allgrains);
for i = 1:grainlength, %
Add grain force to global force vector
num = (allgrains(i)-1)*3 + 3;
Fglobal(num) = Fglobal(num) - GammaB;
end
center = (centergrainnode-1)*3+3; %
Add Center grain force to global force vector
Fglobal(center) = Fglobal(center) - GammaC;
end % End setting up Fglobal and Hglobal for ALL elements
clampedlength = length(clampednode);
for i = 1:clampedlength, %
Clamp nodes by adding large number to diagonal elements
for k = 1:3,
clamp = (clampednode(i)-1)*3
+ k;
Hglobal(clamp,clamp)
= Hglobal(clamp,clamp) + 10^6;
end
end
for i = 1:(totalnodes*3), %
add small number to diagonal elements to avoid singularity
Hglobal(i,i) = Hglobal(i,i) + 0.001;
end
velocity = inv(Hglobal)*Fglobal; %
Get velocity vector
for i=1:totalnodes, % Update node positions
node(i,2) = node(i,2) + velocity((i-1)*3+1)*DeltaT;
node(i,3) = node(i,3) + velocity((i-1)*3+2)*DeltaT;
node(i,4) = node(i,4) + velocity((i-1)*3+3)*DeltaT;
end
it
%end % End the whole thing
time = DeltaT*it/60; %
Total simulation time
def(it,1) = node(3,4); %
Get Center point z-coordinate
def(it,2) = time;
if it == 100 | it == 200 | it == 400 | it ==
600 | it == 800 | it == 1000 | it == 2000 | it == 4000 | it == 6000
| it == 8000 | it == 10000 %Draw model at
these time steps
% DRAW the new surface!!!!
figure
view(3)
for i = 1:ELnum,
DRnode1=element(i,2);
DRnode2=element(i,3);
DRnode3=element(i,4);
DR1 = [node(DRnode1,2); %
Set drawing node 1,2,3 to global coordinates
node(DRnode1,3);
node(DRnode1,4)];
DR2 = [node(DRnode2,2);
node(DRnode2,3);
node(DRnode2,4)];
DR3 = [node(DRnode3,2);
node(DRnode3,3);
node(DRnode3,4)];
Xdr = [DR1(1) DR2(1) DR3(1)];
Ydr = [DR1(2) DR2(2) DR3(2)];
Zdr = [DR1(3) DR2(3) DR3(3)];
patch(Xdr,Ydr,Zdr,'g');
end
xlabel('X-Axis Position of Nodes');
ylabel('Y-Axis Position of Nodes');
zlabel('Z-Axis Position of Nodes');
title(['3D Interface Migration of Ice Block after
',num2str(time),' minutes'])
grid
axis([-2 2 -2 2 -.2 0]);
end
end % End the whole
thing