拓扑优化方面的经典99行程序

xinyuxf

老师说这个是拓扑优化的入门程序，要看明白的，个人感觉不错，

1. 
   
2. %%%% A 99 LINE TOPOLOGY OPTIMIZATION CODE BY OLE SIGMUND, JANUARY 2000 %%%
   
3. %%%% CODE MODIFIED FOR INCREASED SPEED, September 2002, BY OLE SIGMUND %%%
   
4. %%%% 一个由 OLE SIGMUND编写的99行拓扑优化代码，2000年1月 %%%
   
5. %%%% 为加速而修改的代码，2002年9月,由OLE SIGMUND编写 %%%
   
6. function top(nelx,nely,volfrac,penal,rmin);
   
7. % INITIALIZE
   
8. x(1:nely,1:nelx) = volfrac;
   
9. loop = 0;
   
10. change = 1.;
    
11. % START ITERATION
    
12. while change > 0.01  
    
13.   loop = loop + 1;
    
14.   xold = x;
    
15. % FE-ANALYSIS
    
16.   [U]=FE(nelx,nely,x,penal);         
    
17. % OBJECTIVE FUNCTION AND SENSITIVITY ANALYSIS
    
18.   [KE] = lk;
    
19.   c = 0.;
    
20.   for ely = 1:nely
    
21.     for elx = 1:nelx
    
22.       n1 = (nely+1)*(elx-1)+ely;
    
23.       n2 = (nely+1)* elx   +ely;
    
24.       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);
    
25.       c = c + x(ely,elx)^penal*Ue'*KE*Ue;
    
26.       dc(ely,elx) = -penal*x(ely,elx)^(penal-1)*Ue'*KE*Ue;
    
27.     end
    
28.   end
    
29. % FILTERING OF SENSITIVITIES
    
30.   [dc]   = check(nelx,nely,rmin,x,dc);   
    
31. % DESIGN UPDATE BY THE OPTIMALITY CRITERIA METHOD
    
32.   [x]    = OC(nelx,nely,x,volfrac,dc);
    
33. % PRINT RESULTS
    
34.   change = max(max(abs(x-xold)));
    
35.   disp([' It.: ' sprintf('%4i',loop) ' Obj.: ' sprintf('%10.4f',c) ...
    
36.        ' Vol.: ' sprintf('%6.3f',sum(sum(x))/(nelx*nely)) ...
    
37.         ' ch.: ' sprintf('%6.3f',change )])
    
38. % PLOT DENSITIES  
    
39.   colormap(gray); imagesc(-x); axis equal; axis tight; axis off;pause(1e-6);
    
40. end
    
41. %%%%%%%%%% OPTIMALITY CRITERIA UPDATE %%%%%%%%%%%%%%%%%%%%%%%%%%%%
    
42. function [xnew]=OC(nelx,nely,x,volfrac,dc)  
    
43. l1 = 0; l2 = 100000; move = 0.2;
    
44. while (l2-l1 > 1e-4)
    
45.   lmid = 0.5*(l2+l1);
    
46.   xnew = max(0.001,max(x-move,min(1.,min(x+move,x.*sqrt(-dc./lmid)))));
    
47.   if sum(sum(xnew)) - volfrac*nelx*nely > 0;
    
48.     l1 = lmid;
    
49.   else
    
50.     l2 = lmid;
    
51.   end
    
52. end
    
53. %%%%%%%%%% MESH-INDEPENDENCY FILTER %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    
54. function [dcn]=check(nelx,nely,rmin,x,dc)
    
55. dcn=zeros(nely,nelx);
    
56. for i = 1:nelx
    
57.   for j = 1:nely
    
58.     sum=0.0;
    
59.     for k = max(i-floor(rmin),1):min(i+floor(rmin),nelx)
    
60.       for l = max(j-floor(rmin),1):min(j+floor(rmin),nely)
    
61.         fac = rmin-sqrt((i-k)^2+(j-l)^2);
    
62.         sum = sum+max(0,fac);
    
63.         dcn(j,i) = dcn(j,i) + max(0,fac)*x(l,k)*dc(l,k);
    
64.       end
    
65.     end
    
66.     dcn(j,i) = dcn(j,i)/(x(j,i)*sum);
    
67.   end
    
68. end
    
69. %%%%%%%%%% FE-ANALYSIS %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    
70. function [U]=FE(nelx,nely,x,penal)
    
71. [KE] = lk;
    
72. K = sparse(2*(nelx+1)*(nely+1), 2*(nelx+1)*(nely+1));
    
73. F = sparse(2*(nely+1)*(nelx+1),1); U = zeros(2*(nely+1)*(nelx+1),1);
    
74. for elx = 1:nelx
    
75.   for ely = 1:nely
    
76.     n1 = (nely+1)*(elx-1)+ely;
    
77.     n2 = (nely+1)* elx   +ely;
    
78.     edof = [2*n1-1; 2*n1; 2*n2-1; 2*n2; 2*n2+1; 2*n2+2; 2*n1+1; 2*n1+2];
    
79.     K(edof,edof) = K(edof,edof) + x(ely,elx)^penal*KE;
    
80.   end
    
81. end
    
82. % DEFINE LOADS AND SUPPORTS (HALF MBB-BEAM)
    
83. F(2,1) = -1;
    
84. fixeddofs   = union([1:2:2*(nely+1)],[2*(nelx+1)*(nely+1)]);
    
85. alldofs     = [1:2*(nely+1)*(nelx+1)];
    
86. freedofs    = setdiff(alldofs,fixeddofs);
    
87. % SOLVING
    
88. U(freedofs,:) = K(freedofs,freedofs) \ F(freedofs,:);      
    
89. U(fixeddofs,:)= 0;
    
90. %%%%%%%%%% ELEMENT STIFFNESS MATRIX %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    
91. function [KE]=lk
    
92. E = 1.;
    
93. nu = 0.3;
    
94. k=[ 1/2-nu/6   1/8+nu/8 -1/4-nu/12 -1/8+3*nu/8 ...
    
95.    -1/4+nu/12 -1/8-nu/8  nu/6       1/8-3*nu/8];
    
96. KE = E/(1-nu^2)*[ k(1) k(2) k(3) k(4) k(5) k(6) k(7) k(8)
    
97.                   k(2) k(1) k(8) k(7) k(6) k(5) k(4) k(3)
    
98.                   k(3) k(8) k(1) k(6) k(7) k(4) k(5) k(2)
    
99.                   k(4) k(7) k(6) k(1) k(8) k(3) k(2) k(5)
    
100.                   k(5) k(6) k(7) k(8) k(1) k(2) k(3) k(4)
     
101.                   k(6) k(5) k(4) k(3) k(2) k(1) k(8) k(7)
     
102.                   k(7) k(4) k(5) k(2) k(3) k(8) k(1) k(6)
     
103.                   k(8) k(3) k(2) k(5) k(4) k(7) k(6) k(1)];
     
104. %
     
105. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     
106. % This Matlab code was written by Ole Sigmund, Department of Solid         %
     
107. % Mechanics, Technical University of Denmark, DK-2800 Lyngby, Denmark.     %
     
108. % Please sent your comments to the author: [email]sigmund@fam.dtu.dk[/email]%
     
109. %                                                                          %
     
110. % The code is intended for educational purposes and theoretical details    %
     
111. % are discussed in the paper                                               %
     
112. % "A 99 line topology optimization code written in Matlab"                 %
     
113. % by Ole Sigmund (2001), Structural and Multidisciplinary Optimization,    %
     
114. % Vol 21, pp. 120--127.                                                    %
     
115. %                                                                          %
     
116. % The code as well as a postscript version of the paper can be             %
     
117. % downloaded from the web-site: [url]http://www.topopt.dtu.dk[/url]        %
     
118. %                                                                          %
     
119. % Disclaimer:                                                              %
     
120. % The author reserves all rights but does not guaranty that the code is    %
     
121. % free from errors. Furthermore, he shall not be liable in any event       %
     
122. % caused by the use of the program.                                        %
     
123. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     

复制代码

[ 本帖最后由 ChaChing 于 2010-1-30 23:35 编辑 ]

macle888

请问这个程序用的是什么语言？？

suffer

原帖由 macle888 于 2006-11-3 17:31 发表
请问这个程序用的是什么语言？？

在这个版的显然是matlab的

liweidlut

是拓扑大家Sigmund的程序，经典

bitcader

可惜有限元分析错了，y方向的位移反了，大家可以改改看

86029

拓扑优化，matlab程序，好东西！ 是很经典啊，可是我水平低，看不懂啊！~

拓扑优化研究方法综述
收稿日期：2005-07-05
陈敏志张旭明徐冯君
摘要：
根据两个基本原理把拓扑优化方法分为退化法和进化法两大类，分别介绍了各主要的拓扑优化设计方法，对这些方法进行了比较，并指出各方法适用的范围，介绍了当前拓扑优化的研究热点，并对结构拓扑优化的未来作了展望。
关键词：结构拓扑优化，退化法，进化法，优化算法
中图分类号：TU311.4 文献标识码：A
1 历史及发展概况
结构拓扑优化〔1，2〕是近20年来从结构优化研究中派生出来的新分支，它在计算结构力学中已经被认为是最富挑战性的一类研究工作。目前有关结构拓扑优化的工程应用研究还很不成熟，在国外处在发展的初期，尤其在国内尚属于起步阶段。1904年Michell在桁架理论中首次提出了拓扑优化的概念。自1964年Dorn等人提出基结构法，将数值方法引入拓扑优化领域，拓扑优化研究开始活跃。20世纪80年代初，程耿东和N.Olhoff在弹性板的最优厚度分布研究中首次将最优拓扑问题转化为尺寸优化问题，他们开创性的工作引起了众多学者的研究兴趣。1988年Bendsoe和Kikuchi发表的基于均匀化理论的结构拓扑优化设计，开创了连续体结构拓扑优化设计研究的新局面。1993年XieYM和StevenGP提出了渐进结构优化法。1999年Bendsoe和Sigmund证实了变密度法物理意义的存在性。2002年罗鹰等提出三角网格进化法，该方法在优化过程中实现了退化和进化的统一，提高了优化效率。
2 工程背景及基本原理
通常把结构优化按设计变量的类型划分成三个层次：结构尺寸优化、形状优化和拓扑优化。尺寸优化和形状优化已得到充分的发展，但它们存在着不能变更结构拓扑的缺陷。在这样的背景下，人们开始研究拓扑优化。拓扑优化的基本思想是将寻求结构的最优拓扑问题转化为在给定的设计区域内寻求最优材料的分布问题。寻求一个最佳的拓扑结构形式有两种基本的原理：一种是退化原理，另一种是进化原理。退化原理的基本思想是在优化前将结构所有可能杆单元或所有材料都加上，然后构造适当的优化模型，通过一定的优化方法逐步删减那些不必要的结构元素，直至最终得到一个最优化的拓扑结构形式。进化原理的基本思想是把适者生存的生物进化论思想引入结构拓扑优化，它通过模拟适者生存、物竞天择、优胜劣汰等自然机理来获得最优的拓扑结构。
3 结构拓扑优化设计方法
目前常使用的拓扑优化设计方法可以分为两大类：退化法和进化法。退化法即传统的拓扑优化方法，一般通过求目标函数导数的零点或一系列迭代计算过程求最优的拓扑结构。目前常用于拓扑优化的退化法有基结构方法、均匀化方法、变密度法、变厚度法等。基结构方法（GSA）的思路是假定对于给定的桁架节点，在每
两个节点之间用杆件连结起来，得到的结构称为基结构。按照某种规则或约束，将一些不必要的杆件从基本结构中删除，认为最终剩下的构件决定了结构的最佳拓扑。基结构方法更适合于桁架和框架结构的拓扑优化。基结构法是在有限的子空间内寻优，容易丢失最优解，另外还存在组合爆炸、解的奇异性等问题。均匀化方法（HM）引入微结构的单胞，通过优化计算确定其材料密度分布，并由此得出最优的拓扑结构。均匀化方法主要应用于连续体的拓扑优化设计，它不仅能用于应力约束和位移约束，也能用于频率约束。目前用均匀化方法来进行拓扑优化设计的有一般弹性问题、热传导问题、周期渐进可展曲面问题、非线性热弹性问题、振动问题和骨改造问题等。变密度法〔3〕是一种比较流行的力学建模方式，与采用尺寸变量相比，它更能反映拓扑优化的本质特征。通常，单元密度与弹性模量之间的关系采用人为给出的幂函数规律，亦有采用p的有理分式形式或者E1和E0的组合形式。
目前已有研究将变密度法应用于卧式千斤顶〔4〕以及磁场等拓扑优化设计中，取得了一定的成效。另外，人们对变密度法的不足加以改进，研究出SIMP法以及能量法则等方法，提高了计算能力及应用水平。变厚度法采用满应力法，进行有限元分析得到的各单元在节点处的Misses相当于应力，将围绕每一节点i的所有单元在节点i处应力的加权平均值作为接点i的应力σi。通过迭代不断改变各节点处的厚度使其应力趋近最大的允许值，达到满应力的设计目的。当节点处的厚度低于事先设定值时，节点即被删除。此外，两相法、内力法〔5〕、泡泡法等退化方法也开始被应用于拓扑优化上。
进化法是一类全局寻优方法，目前常用于拓扑优化的进化法主要有遗传算法、模拟退火算法和渐进结构优化法等。遗传算法（GA）是一种新型的基于遗传进化机理的寻优技术，它通过选择、交叉、变异等过程可使群体性能趋于最佳，从而获得全局最优解。作为拓扑优化设计方法之一，GA主要应用于建筑结构优化方面，如桁架结构优化设计、抗震结构智能优化设计等。近几年出现的改进GA法，如自适应GA的研究和复合型GA，使GA法得到进一步发展。
模拟退火算法（SA）〔6〕的思想源于固体退火过程：固体在加温时，其内部粒子随温升变为无序状，内能增大，而徐徐冷却时粒子渐趋有序，在每个温度都达到平衡状态，最后在常温时达到基态，内能减为最小。SA法是最近十年来才应用于拓扑优化的，应用得比较广泛的是桁架的拓扑优化，最近几年将非线性问题转化为
多目标SA法的研究也取得一定的进展〔7〕，SA法适合解决比较复杂的问题。
渐进结构优化法（ESO法）的基本概念很简单，即通过将无效的或低效的材料一步步去掉，结构也将逐步趋于优化。可用ESO法优化的结构为桁架、刚架、板壳或三维连续体等，优化的约束条件包括应力、刚度、位移、频率、临界压力及动响应力等。目前ESO法已被成功应用于热应力优化、热传导优化、触应力优化、可
变荷载结构优化以及机床床身的优化〔8〕等。另外，双向拓扑优化法（BESO）、神经元网络法、自然生长方
法、节点增加方法以及极大熵原理法、软杀法和硬杀法等进化方法也开始被应用于拓扑优化上，并取得了瞩目的进展。此外，为了克服退化法和进化法各种算法自身的不足，一种有效的方法就是采用混合策略，即将两种或两种以上优化方法相结合，如三角网格进化法、复合形遗传算法等等，这里就不再一一介绍了。
4 拓扑优化设计方法比较
1）用于拓扑优化的优化算法比较
目前拓扑优化设计方法常采用的优化算法主要有两类：优化准则法（OC）和数学规划法（MP）。拓扑优化中常用的OC法是满应力准则法，其最大的特点是收敛速度快，要求重分析的次数一般跟设计变量的数目关系不大。拓扑优化中常用MP法是移动渐进算法（MMA法）。对单约束拓扑优化问题，算法一般选用优化准则法；对模型中包含多约束情况或多工况的拓扑优化问题，优化准则法的推导过程复杂，计算和更新拉格朗日乘子比较麻烦，这时应优先选用MP法。
2）几种主要退化法比较
基结构方法中的应力的不连续性会引起奇异最优解现象。基结构方法更适合于桁架和框架结构的拓扑优化。均匀化方法不仅能用于应力约束、位移约束，也能用于频率约束，但其设计变量多，敏度计算复杂，优化后的结构常常含有多孔质材料。变密度法基于各向同性材料，程序实现简单，计算效率高。但是有一点需要说明：它是人为假设的，甚至是基于经验的。
3）几种主要进化法比较
遗传算法不受初始值的影响，其缺点是搜索时间过长、易发生早熟收敛等；模拟退火算法具有较强的全局搜索能力，但它存在着最后结果可能比中间结果差的问题。渐进结构优化法简单通用，但收敛性较差，且优化过程中误删除单元后不能再恢复，双向拓扑优化法则可以弥补此不足。
4）退化法和进化法之间的比较
退化法一般通过求目标函数导数的零点或一系列迭代计算过程求最优解，易陷入局部最优解，且要求目标函数有较好的连续性和可微性。进化法既不要求连续又不要求可微，有较强的全局寻优能力，但需要花费较长的时间，而且它们没有固定的理论背景，其收敛性也未被充分证明。
5）拓扑优化方法面临的问题
拓扑优化的主要困难在于其可行域的奇异性，导致了拓扑优化的全局最优解与局部最优解之间存在很大差异。另外，基于有限元法求解拓扑优化问题，在优化结构中会经常出现棋盘格式、网格依赖性、中间密度材料等数值计算问题。
5 应用领域及发展前景
结构拓扑优化设计研究，已被广泛应用于建筑、航天航空、机械、海洋工程、生物医学及船舶制造等领域。目前结构拓扑优化设计技术的研究热点主要有以下几方面：
1）基于无网格数值技术的拓扑优化设计；并行结构拓扑优化设计技术；双向拓扑优化技术；复合形遗传算法等混合拓扑优化设计技术；
2）多目标拓扑优化设计；结构动力学的拓扑优化设计；非线形的拓扑优化设计〔7〕；可变荷载拓扑优化设计；多工况下的拓扑优化设计；
3）路径规划非线形控制柔性机构的拓扑优化设计；压电智能柔性机构以及柔性力学结构的拓扑优化设计；磁场的拓扑优化设计〔4〕；抗震结构智能优化设计。结构拓扑优化设计虽然已在许多方面取得成果，但拓扑优化设计在实际工程中的应用还需要进一步的较为深入的研究。
参考文献：
〔1〕程耿东.结构优化新方法及其计算机实现〔J〕.力学与实践，1992，14（1）：1-6.
〔2〕许素强，夏人伟.结构优化方法研究方法综述〔J〕.航空学报，1995，16（4）：385-396.
〔3〕邓扬晨.力学建模在变密度法中密度与刚度关系中的应用〔J〕.机械科学与技术，2002（22）：95-98.
〔4〕孟明，汤文成.卧式千斤顶墙板的拓扑优化设计〔J〕.起重运输机械，2002（5）：10-12.
〔5〕汪树玉，刘国华，包志仁.结构优化设计的现状与进展〔J〕.基建优化，1999，20（4）：3-14.
〔6〕蔡文学，程耿东.桁架结构拓扑优化设计的模拟退火算法〔J〕.华南理工大学学报，1998，29（9）：79-84.
〔7〕张长林，余建星.非线性约束最优化问题的多目标模拟退火算法〔J〕.复旦学报，2003，42（1）：93-97.
〔8〕倪晓宇，易红.机床床身结构的有限元分析与优化〔J〕.制造技术与机床，2005（2）：47-50.

[ 本帖最后由 ChaChing 于 2010-1-31 14:55 编辑 ]

蝴蝶兰_2009

请问楼主这个怎么运行啊    直接保存成M文件再运行 怎么运行不了呢  我水平很低  还望不吝赐教

86029

原帖由 bitcader 于 2008-1-21 11:01 发表

                               
登录/注册后可看大图

可惜有限元分析错了，y方向的位移反了，大家可以改改看

Y的方向是不是跟坐标系选取有关系。他取得是竖直向上，而力是竖直向下的，不知道理解对不？

86029

75行左右的组集总刚度矩阵
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

这是在组集总刚度矩阵，怎么个对号入坐。等号后第一项是在累加，有点没看出来是在“对号入座”

等号右边的 K(edof,edof)是指0矩阵中的EDOF行和EDOF列的0元素换成x(ely,elx)^penal*KE，就是这样的对号入座，不知道对不？

[ 本帖最后由 ChaChing 于 2010-5-14 20:54 编辑 ]

ajiang

回复 9 # 86029 的帖子

那个程序是 谁发的呀 ，小弟正困惑程序的 难懂，可以交流下吗，我研究的 就是 优化算法编程，愚昧我 都看不懂这程序，哪位大哥精通拓扑优化程序

freebeita

优化理论与拓扑优化编程还是有点差别的，楼上的会多目标优化不~~

ajiang

谢谢楼上的freebeita ，现在才看到回帖，呵呵，最近看懂了 不少，我 q：358939206,有做拓扑优化研究的同学可以多联系哦 ，大家交流下一起进步

86029

ajiang 发表于 2010-11-14 22:04

                               
登录/注册后可看大图

回复 9 # 86029 的帖子

那个程序是 谁发的呀 ，小弟正困惑程序的 难懂，可以交流下吗，我研究的 就是 优化 ...

我的也是算法，但是才看见你的留言，球球：361191113

骨头19901119

运行说第三行就错了，v那个没定义

wangzhong45

这个程序还是有点不明白，对于复杂一点的结构，怎么弄啊？