拓扑优化经典99行程序解读
本帖最后由 yygao56 于 2011-7-10 22:02 编辑Sigmund教授所编写的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;
% INITIALIZE
x(1:nely,1:nelx) = volfrac;
loop = 0;
change = 1.;
% START ITERATION
while change > 0.01
loop = loop + 1;
xold = x;
% FE-ANALYSIS
=FE(nelx,nely,x,penal);
% OBJECTIVE FUNCTION AND SENSITIVITY ANALYSIS
= lk;
c = 0.;
for ely = 1:nely
for elx = 1:nelx
n1 = (nely+1)*(elx-1)+ely;
n2 = (nely+1)* elx +ely;
Ue = U(,1);
c = c + x(ely,elx)^penal*Ue'*KE*Ue;
dc(ely,elx) = -penal*x(ely,elx)^(penal-1)*Ue'*KE*Ue;
end
end
% FILTERING OF SENSITIVITIES
= check(nelx,nely,rmin,x,dc);
% DESIGN UPDATE BY THE OPTIMALITY CRITERIA METHOD
= OC(nelx,nely,x,volfrac,dc);
% 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 =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;
end
end
%%%%%%%%%% MESH-INDEPENDENCY FILTER %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function =check(nelx,nely,rmin,x,dc)
dcn=zeros(nely,nelx);
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(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);
end
end
%%%%%%%%%% FE-ANALYSIS %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function =FE(nelx,nely,x,penal)
= 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 = ;
K(edof,edof) = K(edof,edof) + x(ely,elx)^penal*KE;
end
end
% DEFINE LOADS AND SUPPORTS (HALF MBB-BEAM)
F(2*(nelx/2+1)*(nely+1),1) = 1;
fixeddofs = ;
alldofs = ;
freedofs = setdiff(alldofs,fixeddofs);
% SOLVING
U(freedofs, : )= K(freedofs,freedofs) \ F(freedofs,: );
U(fixeddofs,: )= 0; %这两行复制后应换成英文字符,我这里为了防止生成QQ表
% 情修改了一下格式
%%%%%%%%%% ELEMENT STIFFNESS MATRIX %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function =lk
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/8nu/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)];
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
程序执行方法:(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; %敏度过滤的半径
% INITIALIZE
x(1:nely,1:nelx) = volfrac; %x是设计变量
loop = 0; %存放迭代次数的变量
change = 1.; %每次迭代目标函数的改变值,用以判断何时收敛
% START ITERATION
while change > 0.01 %当两次连续目标函数迭代的差小于等于0.01时,结束迭代
loop = loop + 1; %迭代次数加1
xold = x; %把前一次的设计变量赋值给xold
% FE-ANALYSIS
=FE(nelx,nely,x,penal); %进行有限元分析
% OBJECTIVE FUNCTION AND SENSITIVITY ANALYSIS
= lk; %单元刚度矩阵
c = 0.; %用来存放目标函数的变量.这里目标函数是刚度最大,也就是柔
%度最小
for ely = 1:nely
for elx = 1:nelx
n1 = (nely+1)*(elx-1)+ely;
n2 = (nely+1)* elx +ely;
Ue = U(,1);
c = c + x(ely,elx)^penal*Ue'*KE*Ue; %计算目标函数的值(即柔度
%是多少)
dc(ely,elx) = -penal*x(ely,elx)^(penal-1)*Ue'*KE*Ue; % 灵敏度分析的结果 这一行
%和上一行可参考论文中的公式
end
end
% FILTERING OF SENSITIVITIES
= check(nelx,nely,rmin,x,dc); %灵敏度过滤,为了边界光顺一点
% DESIGN UPDATE BY THE OPTIMALITY CRITERIA METHOD
= OC(nelx,nely,x,volfrac,dc); %优化准则法更新设计变量
% 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
OC算法子程序:
function =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;
end
end
敏度过滤技术子程序:
function =check(nelx,nely,rmin,x,dc)
dcn=zeros(nely,nelx);
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(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);
end
end
这一段就不多解释了,只是为了光顺边界的,现在二重敏度过滤技术用得更多一点了。理论部分可参考罗震博士的毕业论文
看不懂的代码在matlab命令窗口输入:help XXX(即你看不懂的那个关键词,比如这里的floor)
有限元求解子程序
function =FE(nelx,nely,x,penal)
= 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 = ; %这里的Y轴是反向的,但是不影响最后的结果,详情请见二楼TYNGOD这位高手的解释,感谢TYNGOD。
K(edof,edof) = K(edof,edof) + x(ely,elx)^penal*KE; %将单元刚度矩阵组装成总的刚度矩阵
end
end
% DEFINE LOADS AND SUPPORTS (HALF MBB-BEAM)
F(2*(nelx/2+1)*(nely+1),1) = 1; %初始的集中力
fixeddofs = ;%固定结点
alldofs = ; %所有结点
freedofs = setdiff(alldofs,fixeddofs); %自由节点
% SOLVING
U(freedofs,:) = K(freedofs,freedofs) \ F(freedofs,:); %有限元求解:位移场
U(fixeddofs,:)= 0; %固定节点位移为0
单元刚度矩阵的子程序:
function =lk
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/8nu/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)];
这里不解释,自己可查有限元理论。个人认为这里的单元刚度求解不是很好。
最后是我的运行结果:
优化结果:
论文:A 99 line topology optimization code written in Matlab
如果不看论文,只看我这个程序,是读不懂的。建议新手先大概了解一下拓扑优化的基本知识 3# 417332551
有空把你的算法程序发给我看看吧!
邮箱:yygao56@smail.hust.edu.cn
多谢啦
yygao56 发表于 2011-7-11 18:58 http://forum.simwe.com/images/common/back.gif
见习版主原来是华科的啊,校友好! 请大家注意Sigmund教授发表的代码中,定义坐标轴的Y轴方向为坚直向下,与通常我们熟悉的Y轴竖直向上为正是相反的。从文中和代码中都可以证实,说明如下:
文中指出:The variables n1 and n2 denote upper left and right element node numbers in global node numbers。也就是说n1和n2两个节点分别是四边形单元的左上角点和右上角点。如下图所示,
n1++++++++++++++n2
| |
| |
| |
| |
n4++++++++++++++n3( n1 -> n2 -> n3 -> n4 )
也就是说,Sigmund指定单元四个节点的顺序为顺时针,分别为:左上,右上,右下,左下 四个节点。(通常我们是从左下节点开始逆时针方向定义节点编号的)所示单元的自由度分别有:
edof =
所以,我们可以证实单元内四个节点的编号(n1, n2, n3, n4)有如下的关系:
n1+1=n4;
n2+1=n3;
n1+(nely+1)=n2; ( nely 为结构网格划分后竖向的单元数目)
所以说明,Sigmund教授定义的节点编号是竖直向下编号的,则Y轴是竖直向下为正,其原点为结构的左上角点。
但对于这种刚度优化问题,这个定义是不会影响柔顺度数值的。
此问题已得到学者证实!! 本帖最后由 417332551 于 2011-7-10 20:54 编辑
哥哥,可以用不同的颜色标注一下。。
你真是个大好人。一定要顶你当OS的斑竹。
乐于助人,授人以渔。。才是大侠风范!!!此人日后必成大器。
帖子发表完后,还可以再编辑编辑。。用不同颜色区分。。
另外有一些QQ表情。。
有时间,我给你一部分OC算法和mma算法的源程序。(暂时没空)
楼主,辛苦啦你!
2# TYNGOD
是的。程序中Y轴方向是反的,但是因为是弹性结构,所以对最后的拓扑没有影响。 3# 417332551
嗯,因为我发帖都是先word编辑,或者记事本编辑的,然后复制进来。所以对论坛内发帖的很多功能不太会。我试试改颜色,谢谢你的建议 威武,:victory:版块里一大人才,大家有啥不懂的赶紧问,莫让他闲着,嘿嘿:lol 哈哈哈哈。你干的不错!有心就好。
有心做公益就好!!
[b] 1# yygao56 3# 417332551
有空把你的算法程序发给我看看吧!
邮箱:yygao56@smail.hust.edu.cn
多谢啦 9# yrunze
话说,学校邮箱以前都不用的~~最近才启用的。 本帖最后由 417332551 于 2011-7-28 10:01 编辑
记得查收
你先把这个注释,详细讲解一下。。然后我再给你传资料。。
因为这个是我很久很久以前的师兄做的硕士论文。。。不知道我上传他的东西,会不会出问题。
10# yygao56 本帖最后由 417332551 于 2011-7-28 10:00 编辑
It.: 1 Obj.: 333.0403 Vol.:0.400 ch.:0.200
It.: 2 Obj.: 208.0136 Vol.:0.400 ch.:0.200
It.: 3 Obj.: 151.8037 Vol.:0.400 ch.:0.200
It.: 4 Obj.: 123.7844 Vol.:0.400 ch.:0.200
这个是你的99行代码。。我正看里面的计算公式。
It.: 5 Obj.: 111.9788 Vol.:0.400 ch.:0.173
It.: 6 Obj.: 103.2981 Vol.:0.400 ch.:0.200
It.: 7 Obj.: 97.1750 Vol.:0.400 ch.:0.167
It.: 8 Obj.: 90.7349 Vol.:0.400 ch.:0.200
It.: 9 Obj.: 85.5328 Vol.:0.400 ch.:0.162
It.: 10 Obj.: 80.2626 Vol.:0.400 ch.:0.179
It.: 11 Obj.: 75.5167 Vol.:0.400 ch.:0.154
It.: 12 Obj.: 70.9125 Vol.:0.400 ch.:0.159
It.: 13 Obj.: 66.6570 Vol.:0.400 ch.:0.124
It.: 14 Obj.: 63.3128 Vol.:0.400 ch.:0.127
It.: 15 Obj.: 60.6792 Vol.:0.400 ch.:0.094
It.: 16 Obj.: 59.0614 Vol.:0.400 ch.:0.094
It.: 17 Obj.: 57.8935 Vol.:0.400 ch.:0.082
It.: 18 Obj.: 57.1503 Vol.:0.400 ch.:0.070
It.: 19 Obj.: 56.5448 Vol.:0.400 ch.:0.074
It.: 20 Obj.: 56.1134 Vol.:0.400 ch.:0.057
It.: 21 Obj.: 55.8116 Vol.:0.400 ch.:0.059
It.: 22 Obj.: 55.5296 Vol.:0.400 ch.:0.047
It.: 23 Obj.: 55.3214 Vol.:0.400 ch.:0.048
It.: 24 Obj.: 55.1301 Vol.:0.400 ch.:0.038
It.: 25 Obj.: 55.0099 Vol.:0.400 ch.:0.038
It.: 26 Obj.: 54.8593 Vol.:0.400 ch.:0.034
It.: 27 Obj.: 54.7776 Vol.:0.400 ch.:0.034
It.: 28 Obj.: 54.6600 Vol.:0.400 ch.:0.031
It.: 29 Obj.: 54.5681 Vol.:0.400 ch.:0.029
It.: 30 Obj.: 54.5063 Vol.:0.400 ch.:0.027
It.: 31 Obj.: 54.4336 Vol.:0.400 ch.:0.027
It.: 32 Obj.: 54.3624 Vol.:0.400 ch.:0.025
It.: 33 Obj.: 54.2916 Vol.:0.400 ch.:0.024
It.: 34 Obj.: 54.2521 Vol.:0.400 ch.:0.023
It.: 35 Obj.: 54.1981 Vol.:0.400 ch.:0.022
It.: 36 Obj.: 54.1318 Vol.:0.400 ch.:0.022
It.: 37 Obj.: 54.0924 Vol.:0.400 ch.:0.019
It.: 38 Obj.: 54.0483 Vol.:0.400 ch.:0.020
It.: 39 Obj.: 54.0072 Vol.:0.400 ch.:0.017
It.: 40 Obj.: 53.9663 Vol.:0.400 ch.:0.018
It.: 41 Obj.: 53.9645 Vol.:0.400 ch.:0.016
It.: 42 Obj.: 53.9046 Vol.:0.400 ch.:0.017
It.: 43 Obj.: 53.8879 Vol.:0.400 ch.:0.016
It.: 44 Obj.: 53.8638 Vol.:0.400 ch.:0.016
It.: 45 Obj.: 53.8307 Vol.:0.400 ch.:0.015
It.: 46 Obj.: 53.8200 Vol.:0.400 ch.:0.015
It.: 47 Obj.: 53.8020 Vol.:0.400 ch.:0.015
It.: 48 Obj.: 53.7849 Vol.:0.400 ch.:0.015
It.: 49 Obj.: 53.7713 Vol.:0.400 ch.:0.015
It.: 50 Obj.: 53.7176 Vol.:0.400 ch.:0.015
It.: 51 Obj.: 53.7128 Vol.:0.400 ch.:0.015
It.: 52 Obj.: 53.6724 Vol.:0.400 ch.:0.015
It.: 53 Obj.: 53.6572 Vol.:0.400 ch.:0.015
It.: 54 Obj.: 53.6178 Vol.:0.400 ch.:0.015
It.: 55 Obj.: 53.5669 Vol.:0.400 ch.:0.016
It.: 56 Obj.: 53.5105 Vol.:0.400 ch.:0.017
It.: 57 Obj.: 53.4867 Vol.:0.400 ch.:0.017
It.: 58 Obj.: 53.4436 Vol.:0.400 ch.:0.017
It.: 59 Obj.: 53.3850 Vol.:0.400 ch.:0.017
It.: 60 Obj.: 53.3274 Vol.:0.400 ch.:0.017
It.: 61 Obj.: 53.2989 Vol.:0.400 ch.:0.017
It.: 62 Obj.: 53.2237 Vol.:0.400 ch.:0.016
It.: 63 Obj.: 53.1788 Vol.:0.400 ch.:0.016
It.: 64 Obj.: 53.1438 Vol.:0.400 ch.:0.016
It.: 65 Obj.: 53.0900 Vol.:0.400 ch.:0.015
It.: 66 Obj.: 53.0506 Vol.:0.400 ch.:0.014
It.: 67 Obj.: 53.0082 Vol.:0.400 ch.:0.014
It.: 68 Obj.: 52.9708 Vol.:0.400 ch.:0.013
It.: 69 Obj.: 52.9630 Vol.:0.400 ch.:0.012
It.: 70 Obj.: 52.9436 Vol.:0.400 ch.:0.011
It.: 71 Obj.: 52.9106 Vol.:0.400 ch.:0.010
It.: 72 Obj.: 52.8900 Vol.:0.400 ch.:0.009
慢慢研究研究里面的公式。。。希望你早点上学,回到论坛。。。。祝愿大家将此帖子深加工。研究式的学习和探讨,不能浮于表面代码的复制 粘贴 运行。 F(2*(nelx/2+1)*(nely+1),1) = 1; %初始的集中力
fixeddofs = ;%固定结点
alldofs = ; %所有结点
freedofs = setdiff(alldofs,fixeddofs); %自由节点 13# 417332551
于是乎,你决定要仔细研究算法了? 有些地方看不懂,有些公式不清楚。要查阅很多资料的。
区别一:这些都是二维结构的的。直接定义坐标轴节点数目。工程上结构很复杂,结构不规则。
二、三维结构的话,好像很难定义节点,设计空间。
有点难度!斑竹这个程序和PDF那个程序做了一些相应的变动。
你看看我上传的几个文档。。。。有一个是自动输入坐标轴上的 节点数目以及其他参数,matlab做的一个原型系统。。
不是我做的啊!!!
14# kmani 在此多谢!请问各位前辈有基于ansys的APDL语言的拓扑优化程序吗,ansys自带的拓扑模块无法优化shell63壳单元 斑竹威武,不解释! 蛮不错的,以前没有任何资料,我花了个把月才看懂一个程序,早知道有你,我相信我就不会过得这么纠结,qq358939206,有缘拓扑者可加我多交流 n1++++++++++++++n2
| |
| |
| |
| |
n4++++++++++++++n3( n1 -> n2 -> n3 -> n4 )
前几天看了代码,发现节点的顺序也是这样的,原来是Y轴方向错误,坛子里高手真多。 想请教楼主,99行代码中的敏度推导是不是这样的。以前一直对敏度前面的一个负号耿耿于怀,现在推了一遍,还对倒数第二步有点疑问。