拓扑优化经典99行程序解读.doc

精品文档/thread-993188-1-1.htmlSigmund教授所编写的top优化经典99行程序，可以说是我们拓扑优化研究的基础；每一个新手入门都会要读懂这个程序，才能去扩展，去创新；99行程序也有好多个版本，用于求解各种问题，如刚度设计、柔顺机构、热耦合问题，但基本思路大同小异；本文拟对其中的一个版本进行解读，愿能对新手有点小小的帮助。不详之处，还请论坛内高手多指点读懂了该程序，只能说是略懂拓扑优化理论了，我手里就有一些水平集源程序是成千上万行，虽然在99行的基础上成熟了很多，但依然还有很多的发展空间。源程序如下：% A 99 LINE TOPOLOGY OPTIMIZATION CODE BY OLE SIGMUND, JANUARY 2000 % CODE MODIFIED FOR INCREASED SPEED, September 2002, BY OLE SIGMUND %function top(nelx,nely,volfrac,penal,rmin);nelx=80;nely=20;volfrac=0.4;penal=3;rmin=2;% INITIALIZEx(1:nely,1:nelx) = volfrac; loop = 0; change = 1.;% START ITERATIONwhile change 0.01loop = loop + 1;xold = x;% FE-ANALYSISU=FE(nelx,nely,x,penal); % OBJECTIVE FUNCTION AND SENSITIVITY ANALYSISKE = 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(ely,elx)penal*Ue*KE*Ue; dc(ely,elx) = -penal*x(ely,elx)(penal-1)*Ue*KE*Ue; endend% FILTERING OF SENSITIVITIESdc = check(nelx,nely,rmin,x,dc); % DESIGN UPDATE BY THE OPTIMALITY CRITERIA METHODx = OC(nelx,nely,x,volfrac,dc); % PRINT RESULTSchange = 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 DENSITIEScolormap(gray); imagesc(-x); axis equal; axis tight; axis off;pause(1e-6);end % OPTIMALITY CRITERIA UPDATE %function xnew=OC(nelx,nely,x,volfrac,dc)l1 = 0; l2 = 100000; move = 0.2;while (l2-l1 1e-4)lmid = 0.5*(l2+l1);xnew = max(0.001,max(x-move,min(1.,min(x+move,x.*sqrt(-dc./lmid);if sum(sum(xnew) - volfrac*nelx*nely 0; l1 = lmid;else l2 = lmid;endend% MESH-INDEPENDENCY FILTER %function dcn=check(nelx,nely,rmin,x,dc)dcn=zeros(nely,nelx);for i = 1:nelxfor 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(0,fac); dcn(j,i) = dcn(j,i) + max(0,fac)*x(l,k)*dc(l,k); end end dcn(j,i) = dcn(j,i)/(x(j,i)*sum);endend% FE-ANALYSIS %function U=FE(nelx,nely,x,penal)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:nelxfor 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(edof,edof) = K(edof,edof) + x(ely,elx)penal*KE;endend% DEFINE LOADS AND SUPPORTS (HALF MBB-BEAM)F(2*(nelx/2+1)*(nely+1),1) = 1;fixeddofs = 2*(nely/2+1),2*nelx*(nely+1)+2*(nely/2+1);alldofs = 1:2*(nely+1)*(nelx+1);freedofs = setdiff(alldofs,fixeddofs);% SOLVINGU(freedofs, ： ）= K(freedofs,freedofs) F(freedofs,： ）; U(fixeddofs,： ）= 0; %这两行复制后应换成英文字符，我这里为了防止生成QQ表 % 情修改了一下格式% ELEMENT STIFFNESS MATRIX %function KE=lkE = 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/8nu/6 1/8-3*nu/8;KE = E/(1-nu2)* 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);%程序执行方法：（just for matlab new users）打开matlab，点开new M-file，将上述源程序复制粘贴到M-文件中，修改蓝色部分的格式，保存。按F5即可执行程序个人解读（会针对大家的提问，高手们的解释，不断补充更新）：主程序部分：包括：数据初始化；有限元分析；敏度分析，OC算法，结果显示function top(nelx,nely,volfrac,penal,rmin);nelx=80; % x轴方向的单元数nely=20; % y轴方向单元数volfrac=0.4;%体积比penal=3; %材料插值的惩罚因子rmin=2; %敏度过滤的半径% INITIALIZEx(1:nely,1:nelx) = volfrac; %x是设计变量 loop = 0; %存放迭代次数的变量change = 1.; %每次迭代目标函数的改变值，用以判断何时收敛% START ITERATIONwhile change 0.01 %当两次连续目标函数迭代的差小于等于0.01时，结束迭代 loop = loop + 1; %迭代次数加1 xold = x; %把前一次的设计变量赋值给xold% FE-ANALYSISU=FE(nelx,nely,x,penal); %进行有限元分析 % OBJECTIVE FUNCTION AND SENSITIVITY ANALYSISKE = 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(ely,elx)penal*Ue*KE*Ue; %计算目标函数的值（即柔度 %是多少） dc(ely,elx) = -penal*x(ely,elx)(penal-1)*Ue*KE*Ue; % 灵敏度分析的结果 这一行 %和上一行可参考论文中的公式endend% FILTERING OF SENSITIVITIESdc = check(nelx,nely,rmin,x,dc); %灵敏度过滤，为了边界光顺一点% DESIGN UPDATE BY THE OPTIMALITY CRITERIA METHODx = OC(nelx,nely,x,volfrac,dc); %优化准则法更新设计变量% PRINT RESULTSchange = 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 DENSITIEScolormap(gray); imagesc(-x); axis equal; axis tight; axis off;pause(1e-6); %优化结果的图形显示（个人认为这种图形显示方法很不好，太简单了。比较方便的图形显示应该是：每一次迭代同时显示优化结果、目标函数曲线，然后自动保存每一次的结果）end OC算法子程序：function xnew=OC(nelx,nely,x,volfrac,dc)l1 = 0; l2 = 100000; move = 0.2;%l1，l2用于体积约束的拉格朗日乘子while (l2-l1 1e-4)lmid = 0.5*(l2+l1); xnew = max(0.001,max(x-move,min(1.,min(x+move,x.*sqrt(-dc./lmid); %这里是OC算法的核心所在，具体含义可参考论文中的公式if sum(sum(xnew) - volfrac*nelx*nely 0; %采用了二乘法更新拉格朗日乘子 l1 = lmid;else l2 = lmid;endend敏度过滤技术子程序：function dcn=check(nelx,nely,rmin,x,dc)dcn=zeros(nely,nelx);for i = 1:nelxfor 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(0,fac); dcn(j,i) = dcn(j,i) + max(0,fac)*x(l,k)*dc(l,k); end end dcn(j,i) = dcn(j,i)/(x(j,i)*sum);endend这一段就不多解释了，只是为了光顺边界的，现在二重敏度过滤技术用得更多一点了。理论部分可参考罗震博士的毕业论文看不懂的代码在matlab命令窗口输入：help XXX（即你看不懂的那个关键词，比如这里的floor）有限元求解子程序function U=FE(nelx,nely,x,penal)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:nelxfor 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; %这里的Y轴是反向的，但是不影响最后的结果，详情请见二楼TYNGOD这位高手的解释，感谢TYNGOD。 K(edof,edof) = K(edof,edof) + x(ely,elx)penal*KE; %将单元刚度矩阵组装成总的刚度矩阵endend% DEFINE LOADS AND SUPPORTS (HALF MBB-BEAM)F(2*(nelx/2+1)*(nely+1),1) = 1; %初始的集中力fixeddofs = 2*(nely/2+1),2*nelx*(nely+1)+2*(nely/2+1);%固定结点alldofs = 1:2*(nely+1)*(nelx+1); %所有结点freedofs = setdiff(alldofs,fixeddofs); %自由节点% SOLVINGU(freedofs,：） = K(freedofs,freedofs) F(freedofs,：）; %有限元求解：位移场U(fixeddofs,：）= 0; %固定节点位移为0单元刚度矩阵的子程序：function KE=lkE = 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/8nu/6 1/8-3*nu/8;KE = E/(1-nu2)* 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);这里不解释，自己可查有限元