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
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
| function ParetoOptimalTracing(nelx,nely,desVolFrac,problem,volDecrMax,JIncMax,filterRadius)
% Generate pareto-optimal topologies via fixed point iteration
% Author: Krishnan Suresh; [email]suresh@engr.wisc.edu[/email]
% Reference: "A 199-line Matlab Code for Pareto-Optimal Tracing in Topology Optimization",
% K. Suresh, Vol X, pp xy, Structural and Multidisciplinary Optimization,
% Acknowledgements: "A 99 line topology optimization code written in Matlab"
% by Ole Sigmund (2001), Structural and Multidisciplinary Optimization,
% Vol 21, pp. 120--127.
if (nargin == 0) % default values
nelx = 60;nely = 30; % The grid size for topology optimization
desVolFrac = 0.5; % The final volume fraction desired
problem = 1; % 1 or 2 for cantilevered beam problems
volDecrMax = 0.05*randn(1) % step-size for pareto-tracing
JIncMax = 3; % For steep change in pareto-curve, use additional constraint
filterRadius = 0.8; % Use for smoothening the topological sensitivity field
end
voidEps = 1e-4; % Relative Young's Modulus of void region
filter = fspecial('gaussian', [3 3],filterRadius); % smoothen topological sensitivity field
totalIter = 0;
c = 0;
elemsIn(1:nely,1:nelx) = 1; % intialize the domain
U = FE(nelx,nely,elemsIn,voidEps,problem); % Solve FEA problem
T = ComputeT(U,elemsIn,voidEps); % Compute topological sensitivity
T = filter2(filter,T); % smoothen the field
J(1) = computeCompliance(nelx,nely,elemsIn,voidEps,U); % compute & store compliance
volIndex = 1;volFractions(1) = 1; volfrac = 1; % initialization
volDecrement = volDecrMax; % current decrement of volume fraction
while (volfrac > desVolFrac)
volfrac = volfrac-volDecrement; % move to the next volume fraction
volIndex = volIndex+1;
volFractions(volIndex) = volfrac; % store the volume fraction
iter = 0;
while (iter < 10) % to avoid cycles; typically 10 iterations is sufficient
totalIter= totalIter+1;
[isValid,isParetoOptimal] = analyzeTopology(T,elemsIn);
if ((iter > 0)&&(isValid)&&(isParetoOptimal)) % done with current vol
break
end
% Find the level-set value such that the contour has given vol fraction
value = findContourValueWithVolumeFraction(T,volfrac);
[index] = find(T < value); % eliminate all elements less than this value
elemsIn(1:nely,1:nelx) = 1; % start with the full domain
elemsIn(ind2sub(size(T),index)) = 0; % remove elements
U = FE(nelx,nely,elemsIn,voidEps,problem); % FEA
T = ComputeT(U,elemsIn,voidEps); % Topological Sensitivity
T = filter2(filter,T); % Smoothen the field
iter= iter+1;
end
J(volIndex) = computeCompliance(nelx,nely,elemsIn,voidEps,U);
%mas=repmat(J(volIndex),1);
%mas1=repmat(J(volIndex),1);
value = findContourValueWithVolumeFraction(T,volfrac); % as above
plotContour(T,value,figure(1));
%c = repmat(T,val);
title(['v=' num2str(volfrac) '; J = ' num2str(J(volIndex)) '; #FEA = ' num2str(totalIter)]);
pause(0.001);
dJ = J(volIndex)- J(volIndex-1);
volDecrement = max(volDecrement/5,min(volDecrement,JIncMax*volDecrement/dJ));
end
figure(2); plot(volFractions,J,volFractions,J,'*');
xlabel('Volume'); ylabel('Compliance'); grid on;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [U]=FE(nelx,nely,elemsIn,voidEps,problem)
if (problem == 1) % Cantilevered beam;
fixeddofs = 1:2*(nely+1); % left edge
forcedDof = 2*(nelx+1)*(nely+1)-nely; % y force
elseif (problem == 2) % MBB beam
fixeddofs = 1:2*(nely+1); % left edge
forcedDof = 2*(nelx+1)*(nely+1)-2*nely; % y force
end
K = sparse(2*(nelx+1)*(nely+1), 2*(nelx+1)*(nely+1));
F = sparse(2*(nely+1)*(nelx+1),1); U = zeros(2*(nely+1)*(nelx+1),1);
[KE] = lk;
for elx = 1:nelx
for ely = 1:nely
n1 = (nely+1)*(elx-1)+ely;
n2 = (nely+1)* elx +ely;
edof = [2*n1-1 2*n1 2*n2-1 2*n2 2*n2+1 2*n2+2 2*n1+1 2*n1+2]';
alpha = (1-elemsIn(ely,elx))*voidEps + elemsIn(ely,elx);
K(edof,edof) = K(edof,edof) + alpha*KE;
end
end
F(forcedDof,1) = -1;
alldofs = 1:2*(nely+1)*(nelx+1);
freedofs = setdiff(alldofs,fixeddofs);
U(freedofs,:) = K(freedofs,freedofs) \ F(freedofs,:);
U(fixeddofs,:)= 0;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [KE]=lk % element stiffness
E = 1.; nu = 0.3;
k=[ 1/2-nu/6 1/8+nu/8 -1/4-nu/12 -1/8+3*nu/8 ...
-1/4+nu/12 -1/8-nu/8 nu/6 1/8-3*nu/8];
KE = E/(1-nu^2)*[ k(1) k(2) k(3) k(4) k(5) k(6) k(7) k(8)
k(2) k(1) k(8) k(7) k(6) k(5) k(4) k(3)
k(3) k(8) k(1) k(6) k(7) k(4) k(5) k(2)
k(4) k(7) k(6) k(1) k(8) k(3) k(2) k(5)
k(5) k(6) k(7) k(8) k(1) k(2) k(3) k(4)
k(6) k(5) k(4) k(3) k(2) k(1) k(8) k(7)
k(7) k(4) k(5) k(2) k(3) k(8) k(1) k(6)
k(8) k(3) k(2) k(5) k(4) k(7) k(6) k(1)];
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [T] = ComputeT(U,elemsIn,voidEps)
% Compute the topological sensitivity at the center of each element
[nely,nelx] = size(elemsIn);
gradN =0.5*[-1 1 1 -1;-1 -1 1 1]; % at center
E0 = 1;nu = 0.3;
D0 = 1/(1-nu^2)*[1 nu 0; nu 1 0;0 0 (1-nu)/2]; % plane stress
T(1:nely,1:nelx) = 0; % initialize to 0
for elx = 1:nelx
for ely = 1:nely
n1 = (nely+1)*(elx-1)+ely;
n2 = (nely+1)* elx +ely;
edof = [2*n1-1 2*n1 2*n2-1 2*n2 2*n2+1 2*n2+2 2*n1+1 2*n1+2]';
uGrad = gradN*U(edof(1:2:end));
vGrad = gradN*U(edof(2:2:end));
strains = [uGrad(1); vGrad(2); (uGrad(2)+vGrad(1)) ];
alpha = (1-elemsIn(ely,elx))*voidEps + elemsIn(ely,elx);
E = E0*alpha;
stresses = D0*E*strains;
stressTensor = [stresses(1) stresses(3); stresses(3) stresses(2)];
strainTensor = [strains(1) strains(3)/2; strains(3)/2 strains(2)];
if (elemsIn(ely,elx))
T(ely,elx) = 4/(1+nu)*sum(sum(stressTensor.*strainTensor))- ...
(1-3*nu)/(1-nu^2)*trace(stressTensor)*trace(strainTensor);
else
T(ely,elx) = 4/(3-nu)*sum(sum(stressTensor.*strainTensor))+...
(1-3*nu)/((1+nu)*(3-nu))*trace(stressTensor)*trace(strainTensor);
end
end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [J]=computeCompliance(nelx,nely,elemsIn,voidEps,U)
% Compute the compliance of the entire mesh
[KE] = lk;J = 0;
for elx = 1:nelx
for ely = 1:nely
n1 = (nely+1)*(elx-1)+ely;
n2 = (nely+1)* elx +ely;
edof = [2*n1-1 2*n1 2*n2-1 2*n2 2*n2+1 2*n2+2 2*n1+1 2*n1+2]';
alpha = (1-elemsIn(ely,elx))*voidEps + elemsIn(ely,elx);
Ue = U(edof);
J = J + alpha*Ue'*KE*Ue;
end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function value = findContourValueWithVolumeFraction(field,volfrac)
% Find the level-set value such that the contour has given vol fraction
% The code computes the level-set value with desired external volume
[nely,nelx] = size(field);
externalVolumeDesired = nelx*nely*(1-volfrac);
field = -field; % reverse the sign to compute external volume
valMax = 0; valMin = min(field(:));
bufferedField = valMin*ones(nely+2,nelx+2);% Add buffer to get closed contours
bufferedField(2:end-1,2:end-1) = field;
iterMax = 50;iter = 1;
while (1) % A binary search is used to find the optimal level-set value
valMid = (valMax+valMin)/2;
extVol = computeAreaInContour(bufferedField,valMid);
err = abs(extVol-externalVolumeDesired)/(extVol);
if (err < 0.001) || (iter > iterMax)
value = -valMid; % change the sign before return
return;
end
if (extVol > externalVolumeDesired)
valMin = valMid;
else
valMax = valMid;
end
iter = iter+1;
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function area = computeAreaInContour(bufferedField,value)
% For a given level-set value, compute the area enclosed
% It is assumed that the field has been buffered; see code above
[cntr,h] = contours(bufferedField,[value value]);
indices = find(cntr(1,:) == value);area = 0;
for i = 1:numel(indices)
startCol = indices(i)+1;
endCol = startCol+ cntr(2,indices(i))-1;
xPoly = cntr(1,startCol:endCol);
yPoly = cntr(2,startCol:endCol);
area = area + polyarea(xPoly,yPoly);
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [isValid,isParetoOptimal] = analyzeTopology(T,elemsIn)
% Check if the Topological field is valid and/or pareto-optimal
T_SMin = min(T(elemsIn==1)); % Min of topological field inside the domain
T_AMax = max(T(elemsIn==0)); % Max of topological field outside the domain
T_AMin = min(T(elemsIn==0)); % Min of topological field outside the domain
isValid = 0; isParetoOptimal = 0;
if (T_SMin > 0.8*T_AMax) % See paper
isValid = 1;
end
if (T_AMin+T_SMin >= 0) % See paper
isParetoOptimal = 1;
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function plotContour(T,value,fig)
%Use Matlab's built-in contour command to draw the optimal topology.
[nely,nelx] = size(T);
figure(fig);clf;fill([1 nelx nelx 1],[1 1 nely nely],'b'); hold on;
[cntr,h] =contourf(-T,[-value -value]); % the second argument is essential
%disp(T);
axis('equal'); axis tight;axis off; |
(
