%固定兩邊中間的點，在最上方中點施加壓力
function?top（nelxnelyvolfracpenalrmin）;
nelx=80;
nely=50;
volfrac=0.5;
penal=3;
rmin=3;
%?INITIALIZE
x（1:nely1:nelx）?=?volfrac;?
loop?=?0;?
change?=?1.;
%?START?ITERATION
while?change?>?0.01??
??loop?=?loop?+?1;
??xold?=?x;
%?FE-ANALYSIS
??[U]=FE（nelxnelyxpenal）;??
%?objectIVE?FUNCTION?AND?SENSITIVITY?ANALYSIS
??[KE]?=?lk;
??c?=?0.;
??for?ely?=?1:nely
????for?elx?=?1:nelx
??????n1?=?（nely+1）*（elx-1）+ely;?
??????n2?=?（nely+1）*?elx???+ely;
??????Ue?=?U（[2*n1-1;2*n1;?2*n2-1;2*n2;?2*n2+1;2*n2+2;?2*n1+1;2*n1+2]1）;
??????c?=?c?+?x（elyelx）^penal*Ue‘*KE*Ue;
??????dc（elyelx）?=?-penal*x（elyelx）^（penal-1）*Ue‘*KE*Ue;
????end
??end
%?FILTERING?OF?SENSITIVITIES
??[dc]???=?check（nelxnelyrminxdc）;????
%?DESIGN?UPDATE?BY?THE?OPTIMALITY?CRITERIA?METHOD
??[x]????=?OC（nelxnelyxvolfracdc）;?
%?PRINT?RESULTS
??change?=?max（max（abs（x-xold）））;
??disp（[‘?It.:?‘?sprintf（‘%4i‘loop）?‘?Obj.:?‘?sprintf（‘%10.4f‘c）?...
???????‘?Vol.:?‘?sprintf（‘%6.3f‘sum（sum（x））/（nelx*nely））?...
????????‘?ch.:?‘?sprintf（‘%6.3f‘change?）]）
%?PLOT?DENSITIES??
??colormap（gray）;?imagesc（-x）;?axis?equal;?axis?tight;?axis?off;pause（1e-6）;
end?
%%%%%%%%%%?OPTIMALITY?CRITERIA?UPDATE?%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function?[xnew]=OC（nelxnelyxvolfracdc）??
l1?=?0;?l2?=?100000;?move?=?0.2;
while?（l2-l1?>?1e-4）
??lmid?=?0.5*（l2+l1）;
??xnew?=?max（0.001max（x-movemin（1.min（x+movex.*sqrt（-dc./lmid）））））;
??if?sum（sum（xnew））?-?volfrac*nelx*nely?>?0;
????l1?=?lmid;
??else
????l2?=?lmid;
??end
end
%%%%%%%%%%?MESH-INDEPENDENCY?FILTER?%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function?[dcn]=check（nelxnelyrminxdc）
dcn=zeros（nelynelx）;
for?i?=?1:nelx
??for?j?=?1:nely
????sum=0.0;?
????for?k?=?max（i-floor（rmin）1）:min（i+floor（rmin）nelx）
??????for?l?=?max（j-floor（rmin）1）:min（j+floor（rmin）nely）
????????fac?=?rmin-sqrt（（i-k）^2+（j-l）^2）;
????????sum?=?sum+max（0fac）;
????????dcn（ji）?=?dcn（ji）?+?max（0fac）*x（lk）*dc（lk）;
??????end
????end
????dcn（ji）?=?dcn（ji）/（x（ji）*sum）;
??end
end
%%%%%%%%%%?FE-ANALYSIS?%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function?[U]=FE（nelxnelyxpenal）
[KE]?=?lk;?
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）;
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];
????K（edofedof）?=?K（edofedof）?+?x（elyelx）^penal*KE;
??end
end
%?DEFINE?LOADS?AND?SUPPORTS?（HALF?MBB-BEAM）
%固定兩邊中間的點，在最上方中點施加壓力
F（2*（nelx/2+1）*（nely+1）-2*nely1）?=?1;
fixeddofs???=?[2*（nely/2+1）2*nelx*（nely+1）+2*（nely/2+1）];
alldofs?????=?[1:2*（nely+1）*（nelx+1）];
freedofs????=?setdiff（alldofsfixeddofs）;
U（freedofs:）?=?K（freedofsfreedofs）?\?F（freedofs:）;??????
U（fixeddofs:）=?0;
%?SOLVING
%%U（freedofs?:）=?K（freedofsfreedofs）?\?F（freedofs:??;??????
%U（fixeddofs:??=?

?屬性????????????大小?????日期????時間???名稱
-----------?---------??----------?-----??----

?????文件???????3865??2015-10-15?11:21??top99____1.m

-----------?---------??----------?-----??----

?????????????????3865????????????????????1