-
Notifications
You must be signed in to change notification settings - Fork 4
/
homo3D.m
156 lines (155 loc) · 7.03 KB
/
homo3D.m
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
function CH = homo3D(lx,ly,lz,lambda,mu,voxel)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% lx = Unit cell length in x-direction.
% ly = Unit cell length in y-direction.
% lz = Unit cell length in z-direction.
% lambda = Lame's first parameter for solid materials.
% mu = Lame's second parameter for solid materials.
% voxel = Material indicator matrix. Used to determine nelx/nely/nelz
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% INITIALIZE
[nelx, nely, nelz] = size(voxel); %size of voxel model along x,y and z axis
% Stiffness matrix
dx = lx/nelx; dy = ly/nely; dz = lz/nelz;
nel = nelx*nely*nelz;
[keLamda, keMu, feLamda, feMu] = hexahedron(dx/2,dy/2,dz/2);
% Node numbers and element degrees of freedom for full (not periodic) mesh
nodenrs = reshape(1:(1+nelx)*(1+nely)*(1+nelz),1+nelx,1+nely,1+nelz);
edofVec = reshape(3*nodenrs(1:end-1,1:end-1,1:end-1)+1,nel,1);
addx = [0 1 2 3*nelx+[3 4 5 0 1 2] -3 -2 -1];
addxy = 3*(nely+1)*(nelx+1)+addx;
edof = repmat(edofVec,1,24) + repmat([addx addxy],nel,1);
%% IMPOOSE PERIODIC BOUNDARY CONDITIONS
% Use original edofMat to index into list with the periodic dofs
nn = (nelx+1)*(nely+1)*(nelz+1); % Total number of nodes
nnP = (nelx)*(nely)*(nelz); % Total number of unique nodes
nnPArray = reshape(1:nnP, nelx, nely, nelz);
% Extend with a mirror of the back border
nnPArray(end+1,:,:) = nnPArray(1,:,:);
% Extend with a mirror of the left border
nnPArray(:, end+1, :) = nnPArray(:,1,:);
% Extend with a mirror of the top border
nnPArray(:, :, end+1) = nnPArray(:,:,1);
% Make a vector into which we can index using edofMat:
dofVector = zeros(3*nn, 1);
dofVector(1:3:end) = 3*nnPArray(:)-2;
dofVector(2:3:end) = 3*nnPArray(:)-1;
dofVector(3:3:end) = 3*nnPArray(:);
edof = dofVector(edof);
ndof = 3*nnP;
%% ASSEMBLE GLOBAL STIFFNESS MATRIX AND LOAD VECTORS
% Indexing vectors
iK = kron(edof,ones(24,1))';
jK = kron(edof,ones(1,24))';
% Material properties assigned to voxels with materials
lambda = lambda(1)*(voxel==1) + lambda(2)*(voxel==0);
mu = mu(1)*(voxel==1) + mu(2)*(voxel==0);
% The corresponding stiffness matrix entries
sK = keLamda(:)*lambda(:).' + keMu(:)*mu(:).';
K = sparse(iK(:), jK(:), sK(:), ndof, ndof);
K = 1/2*(K+K');
% Assembly load cases corresponding to the strain cases
iF = repmat(edof',6,1);
jF = [ones(24,nel); 2*ones(24,nel); 3*ones(24,nel);...
4*ones(24,nel); 5*ones(24,nel); 6*ones(24,nel);];
sF = feLamda(:)*lambda(:).'+feMu(:)*mu(:).';
F = sparse(iF(:), jF(:), sF(:), ndof, 6);
%% SOLUTION
% solve by PCG method, remember to constrain one node
activedofs = edof(voxel==1,:); activedofs = sort(unique(activedofs(:)));
X = zeros(ndof,6);
try
%Incomplete Cholesky for the K matrix: K(activedofs,activedofs) = L*L'
L = ichol(K(activedofs(4:end),activedofs(4:end))); % L=Lower triangular matrix
for i = 1:6
X(activedofs(4:end),i) = pcg(K(activedofs(4:end),...
activedofs(4:end)),F(activedofs(4:end),i),1e-10,300,L,L');
end
catch
% Solve directly for the displacements U=inv(K(activedofs,activedofs))*F
X(activedofs(4:end),:) = K(activedofs(4:end),activedofs(4:end))...
\F(activedofs(4:end),:);
end
%% HOMOGENIZATION
% The displacement vectors corresponding to the unit strain cases
X0 = zeros(nel, 24, 6);
% The element displacements for the six unit strains
X0_e = zeros(24, 6);
%fix degrees of nodes [1 2 3 5 6 12];
ke = keMu + keLamda; % Here the exact ratio does not matter, because
fe = feMu + feLamda; % it is reflected in the load vector
X0_e([4 7:11 13:24],:) = ke([4 7:11 13:24],[4 7:11 13:24])...
\fe([4 7:11 13:24],:);
X0(:,:,1) = kron(X0_e(:,1)', ones(nel,1)); % epsilon0_11 = (1,0,0,0,0,0)
X0(:,:,2) = kron(X0_e(:,2)', ones(nel,1)); % epsilon0_22 = (0,1,0,0,0,0)
X0(:,:,3) = kron(X0_e(:,3)', ones(nel,1)); % epsilon0_33 = (0,0,1,0,0,0)
X0(:,:,4) = kron(X0_e(:,4)', ones(nel,1)); % epsilon0_12 = (0,0,0,1,0,0)
X0(:,:,5) = kron(X0_e(:,5)', ones(nel,1)); % epsilon0_23 = (0,0,0,0,1,0)
X0(:,:,6) = kron(X0_e(:,6)', ones(nel,1)); % epsilon0_13 = (0,0,0,0,0,1)
CH = zeros(6);
volume = lx*ly*lz;
for i = 1:6
for j = 1:6
sum_L = ((X0(:,:,i) - X(edof+(i-1)*ndof))*keLamda).*...
(X0(:,:,j) - X(edof+(j-1)*ndof));
sum_M = ((X0(:,:,i) - X(edof+(i-1)*ndof))*keMu).*...
(X0(:,:,j) - X(edof+(j-1)*ndof));
sum_L = reshape(sum(sum_L,2), nelx, nely, nelz);
sum_M = reshape(sum(sum_M,2), nelx, nely, nelz);
% Homogenized elasticity tensor
CH(i,j) = 1/volume*sum(sum(sum(lambda.*sum_L + mu.*sum_M)));
end
end
end
%% COMPUTE ELEMENT STIFFNESS MATRIX AND LOAD VECTOR
function [keLambda, keMu, feLambda, feMu] = hexahedron(a, b, c)
% Constitutive matrix contributions
CMu = diag([2 2 2 1 1 1]); CLambda = zeros(6); CLambda(1:3,1:3) = 1;
% Three Gauss points in both directions
xx = [-sqrt(3/5), 0, sqrt(3/5)]; yy = xx; zz = xx;
ww = [5/9, 8/9, 5/9];
% Initialize
keLambda = zeros(24,24); keMu = zeros(24,24);
feLambda = zeros(24,6); feMu = zeros(24,6);
for ii = 1:length(xx)
for jj = 1:length(yy)
for kk = 1:length(zz)
%integration point
x = xx(ii); y = yy(jj); z = zz(kk);
%stress strain displacement matrix
qx = [ -((y-1)*(z-1))/8, ((y-1)*(z-1))/8, -((y+1)*(z-1))/8,...
((y+1)*(z-1))/8, ((y-1)*(z+1))/8, -((y-1)*(z+1))/8,...
((y+1)*(z+1))/8, -((y+1)*(z+1))/8];
qy = [ -((x-1)*(z-1))/8, ((x+1)*(z-1))/8, -((x+1)*(z-1))/8,...
((x-1)*(z-1))/8, ((x-1)*(z+1))/8, -((x+1)*(z+1))/8,...
((x+1)*(z+1))/8, -((x-1)*(z+1))/8];
qz = [ -((x-1)*(y-1))/8, ((x+1)*(y-1))/8, -((x+1)*(y+1))/8,...
((x-1)*(y+1))/8, ((x-1)*(y-1))/8, -((x+1)*(y-1))/8,...
((x+1)*(y+1))/8, -((x-1)*(y+1))/8];
% Jacobian
J = [qx; qy; qz]*[-a a a -a -a a a -a; -b -b b b -b -b b b;...
-c -c -c -c c c c c]';
qxyz = J\[qx;qy;qz];
B_e = zeros(6,3,8);
for i_B = 1:8
B_e(:,:,i_B) = [qxyz(1,i_B) 0 0;
0 qxyz(2,i_B) 0;
0 0 qxyz(3,i_B);
qxyz(2,i_B) qxyz(1,i_B) 0;
0 qxyz(3,i_B) qxyz(2,i_B);
qxyz(3,i_B) 0 qxyz(1,i_B)];
end
B = [B_e(:,:,1) B_e(:,:,2) B_e(:,:,3) B_e(:,:,4) B_e(:,:,5)...
B_e(:,:,6) B_e(:,:,7) B_e(:,:,8)];
% Weight factor at this point
weight = det(J)*ww(ii) * ww(jj) * ww(kk);
% Element matrices
keLambda = keLambda + weight * B' * CLambda * B;
keMu = keMu + weight * B' * CMu * B;
% Element loads
feLambda = feLambda + weight * B' * CLambda;
feMu = feMu + weight * B' * CMu;
end
end
end
end