《拓扑优化99行代码翻译(共13页).docx》由会员分享，可在线阅读，更多相关《拓扑优化99行代码翻译(共13页).docx（13页珍藏版）》请在taowenge.com淘文阁网|工程机械CAD图纸|机械工程制图|CAD装配图下载|SolidWorks_CaTia_CAD_UG_PROE_设计图分享下载上搜索。
1、精选优质文档-倾情为你奉上拓扑优化中的99行matlab代码 o.sigmund摘要 这篇文章描述了用matlab语言来简洁的实现在静态负载下符合最小化原理的拓扑优化。总共只需要输入99行代码，包括优化程序和有限元分析子程序。这99行代码中，其中36行为主程序，12行为基于最优控制器的优化程序，16行为敏度过滤分析，其余35行代码作为有限元分析。实际上，除去注释行以及输出行、有限元分析行，仅有49行matlab代码输入用于解决一个适定的拓扑优化问题。再加上3行补充程序代码，这个程序就可以解决多种负载工况问题。这个代码主要是以教育指导为用，完整的matlab代码在附录中给出，同时也可以在网页上下
2、载。关键词 拓扑优化 教育 最优准则 万维网 matlab代码1 简介文中展示的matlab代码主要是为工程教育所用。学生和在拓扑领域的新手可以在网页上下载。这个代码可以用于结构最优化课程学习，学生们可以在多重负载工况、独立网格选择策略、无源场进行扩展应用。另一种可能就是用来激发学生们的直觉来进行最优化设计。研究生可以推测探究在给定边界条件和容量的情况下的拓扑优化并、比较得出最优策略。在文献中，你可以找到很多处理拓扑优化问题的方法。在一篇Bendse and Kikuchi (1988)的原创论文中，基于对现存结论的学习，所谓微观结构或均化作用的方法被使用。均化作用方法在很多文章中都被采用，但
3、它也存在一些缺点，比如对微观结构最优化方法果断的评估与决策很麻烦的，而且结果很难获得如果没有对微观结构进行确定的长度衡量。然而，在这个意义上来说均化作用方法对拓扑优化也是很重要的，它可以在结构的理论分析上提供一定的界限。拓扑优化的另一种方法叫做“幂律法则”或者SIMP法（(Solid Isotropic Material with Penalization) (Bendse1989; Zhou and Rozvany 1991; Mlejnek 1992)）。这里，假设物质性能使恒定不变的同时每个元素是设计区域离散化，变量是元素的相对密度。物质属性在相对物质密度增加到固体材料的物质属性的很多倍
4、时被模板化。这种方法曾一度引起争议因为它认为没有任何物理材料的属性特征能被幂律法则所描述。然而，Bendse and Sigmund (1999)最近发表的文章证实只要在单一条件下能量能够满足，幂律法则在物理上就是可行的。为了确认这个结论的存在性，幂律法则必须与周界约束、梯度约束或者过滤技术相结合。这个拓扑优化的幂律法则方法已经应用于多重约束，多属性，多材料的问题中了。然而，上面提及到的解决方法是基于数学规划法和连续设计变量法，一系列的文章都有涉及到解决拓扑优化的整数问题。Beckers (1999)通过双重途径成功的解决了大规模服从最小化问题，但是其他方法大多基于遗传算法或者其他为了几个元素
5、需要成千上万的函数求值半随机方法，当然这很可能是不切实际的。除了上述提到的方法外，他们都能解决目标明确的问题，一些减少服从或者目标函数的启发式或直觉的方法也已经被提出来。这些理论都被统称为进化设计理论。除了很容易的理解和运用外，进化性分析方法主要的动机似乎是基于数学或连续变量法“涉及一些微积分应用和数学分析”，并且他们包含“一些复杂问题的数学理论”，反之进化性的方法能很好的应用强大的计算处理技术和自然进化过程中的直觉理念。两件事可以反对他。第一，曾经由于更多的比服从最小化理论复杂的目标被考虑进去。进化论方法自身已经变得非常复杂。第二，正如文中所说，以数学理论为基础的方法解决服从最小化问题很容易
6、实现并且计算处理同样很有效率。不仅如此，基于数学理论的方法很容易扩展解决像非共轭和多物理量的无服从问题以及多约束问题。而用扩展进化方法来处理这些问题似乎更加不可行。完整的matlab代码在附录中给出，文章剩余部分包括对优化问题的定义和讨论（第2部分），matlab实现的详细解释（第3部分），接着是扩展问题的讨论（第4部分）和最后总结（第5部分）。2 拓扑优化问题有许多简化方法都是用来简化matlab的代码。第一，假设设计区域是矩形且被平方有限元离散化。这种情况下，元素的数目和结点就很容易表示（一列一列从左上角开始）并且结构的纵横比通过水平（nelx）和垂直（nely）方向元素比率来确定。一个拓
7、扑优化问题基于指数定律法，目标是实现最小化，可以如下表示（1） 式中U和F分别表示整体变形和力的向量，K表示整体刚度矩阵，和分别表示元素的位移矢量和刚度矩阵，x是设计变量的向量，是相对密度的最小向量（非零以避免奇点）N（=nelx * nely）是设计区域离散化的元素数目，p是惩罚因子（通常p=3），和是材料体积和设计区域体积，f (volfrac)为规定的容积率。 该拓扑优化问题（1）可以用几种不同的方法来解决，如Optimality Criteria（OC算法），Sequential Linear Programming (SLP)，Method of Moving Asymptotes
8、(MMA by Svanberg 1987)或者其他的方法。为简单起见，我们这里用标准的OC算法。根据Bendse (1995)对设计变量应用启发式的调整法，公式表示为其中m（移动量）是一个正的移动界限， (= 1/2)是数值阻尼系数，在最优化条件中表示为其中是拉格朗日乘子，在bi-sectioning algorithm中获得。目标函数的敏感度被表示为更多的详细信息将通过最优化准则理论的推导和实现，在文献(e.g. Bendse 1995).中提到。为了确保对拓扑优化问题（1）解决方法的存在性，将会介绍一些作为结果的限制条件。这里我们使用过滤技术，需要强调的是这个过滤技术并不能证明方法的可行
9、性，但是通过作者大量应用证明过滤技术在实际中能产生独立性网格。这个网格独立性滤波器是来修改元素敏度，如下:这个卷积算子（重量系数）表示为控制函数dist(e, f)定义为元素中心e到元素中心f的距离。在过滤面积以外卷积算子为0.卷积算子随着元素f的距离成线性减少。除了原始的敏度（4），修改的元素敏度（5）被用于OC算法（3）。3 matlab算法实现Matlab代码(见附录),建立了作为一个标准拓扑优化的代码。他的主程序为其中nelx和nely分别是横轴和纵轴方向上元素的个数，volfrac是容积率，penal是惩罚因子，rmin是过滤半径。其他的变量如边界条件在matlab代码中定义且可根据
10、实际需要进行调整。在拓扑优化每一次循环迭代中，程序都会生成一幅当前密度分布图来。图1显示了附录代码在给定以下输入条件下的结果密度分布默认的边界条件对应于“黑梁”的一半 (图1)。垂直方向的负载承受在左上角，水平方向结构左边和右边为对称边界。Matlab代码中重要的细节将在下面各小节讨论。3.1 主程序（1-36行）主程序的开始是材料设计区域的均匀分配（第4行）。一些基本初始化后，主要的循环开始调用有限元分析子程序（12行）并返回位移矢量U。因为固体材料的元素刚度矩阵对任何元素来说都是一样的，元素刚度矩阵子程序仅被调用一次（14行）。接下来，一个循环遍历所有元素(第16- 24行)来确定目标函数
11、和敏感性(4)。变量n1和n2表示在所有节点数中左上角和右上角的元素节点数目，并用于从整体位移矢量U中提取元素位移矢量Ue。灵敏度分析将伴随着调用网格独立性过滤器(第26行)和最优准则优化器(28行)。当前的合规以及其它参数在30-33行中体现出来，并绘制出最后结果的密度分布图(第35行)。主循环会在改变值（改变值在30行确定）小于1%是终止，否则以上步骤将重复执行。3.2 基于优化程序的最优化准则（37-48行）更新后的设计变量是通过优化程序找到的（37-48行）。由于材料体积(sum(sum(xnew)是一个拉格朗日乘子的单调递减函数，满足体积约束的拉格朗日乘子数值可以通过bi-secti
12、oning algorithm算法得到（40-48行）。bi-sectioning algorithm的初始化时通过假设一个小于l1和大于 l2结合的拉格朗日乘数（39行）。限定拉格朗日乘数的间隔大小会一直减半直到它的大小小于收敛判别准则（40行）。3.3 网格独立性过滤技术（49-64行）Matlab中第49-64行是为了实现条件（5）。注意，并不是在设计域里所有的元素都要被搜索来寻找那些在处于过滤半径之内的元素，而是在那些被考虑元素周围，2倍过滤半径大小的区域。通过选择过滤半径rmin小于调用子程序后的半径，被过滤后的敏度将会和初始状态一样，并使得过滤器不起作用。3.4有限元分析代码（65
13、-90行）有限元分析代码位于65-90行中。注意解决者是利用matlab中的稀疏选项。整体刚度矩阵是由所有元素的循环建立的（70-77行）。正如主程序中，变量n1和n2指示在所有节点数目中左上部分和右边节点数目并插入到整体刚度矩阵合适的位置中。正如前面所说，节点和元素都是从左到右纵向排布。而且，每个节点都有2个自由度（水平和垂直的），因此指令F(2,1)=-1.第79行应用了一个在左上角的垂直单元力。消除线性方程中的固定自由度来实现支承结构。Matlab可以很好的处理，由这一行其中freedofs表示不受约束的自由度。通常来说，定义一个固定的自由度要简单的多，其后会被自动找到，通过使用matl
14、ab算子setdiff因无约束自由度与固定自由度的不同来找到无约束自由度（82行）。元素的刚度矩阵在86-99行中进行计算。这个针对4个节点的8*8矩阵可通过使用人工智能软件进行处理分析。杨氏模量E和泊松比可以在88和89行进行改变。4 拓展附录给出的matlab代码解决如图1中材料分配的优化问题以达到最小化原则。而算法中许多的拓展和改变都值得思考，其中一部分将在下面进行阐述。4.1 其他边界条件为解决其他优化问题，改变边界条件和支承条件非常容易。为了解决如图2所示的短悬臂梁问题，只需要要79和80作相应改变，根据这些改变，输入行变为4.2 多重负载问题扩展算法来处理多重负载问题也很容易。实际
15、上，仅需额外增加3行，并对其他4行做微小调整即可。就两个负载工况而言，力和适量位移必须定义为两方向向量，这就意味着第69行改变为现在目标函数为二维之和因此第20-22行即被以下行所替代为解决如图3的二重负载问题，一个右上方的负载要添加到底79行，如下4.3 无源元件 在某些情况下,一些元素可能需要采取最低密度值(例如一个管上的洞)。一个nelx*nely的数组中任何元素都可为0和某些元素固定为0，都可以定义到主程序中或者转移到OC算法中。添加行为在OC子程序中寻找无源元素并设置他们的密度为最小值（0.001）。图4显示的是输入以下行生成的结构结果当这紧接着的10行添加到主程序（第4行后），便可
16、在半径为nely/3的圆中找到无源元素，中心为（nely/2，nelx/3）4.4 可选择的优化器不可否认,基于优化器实现的最优准则只适合单一约束并且它基于启发式定点更新方案。为了建立一个更好的优化方案，我们可以得到matlab版的MMA-algorithm (Svanberg 1987) 来自于 Krister Svanberg,KTH, Sweden.其中总共的变量有20个，包括目标函数，约束条件，新旧的密度等等。实现MMA-algorithm相对简单，但是需要定义几个辅助变量。然而，它却可以解决更多不只一个约束条件的复杂设计问题。在少量增加每次迭代过程CPU响应时间的基础上，Matlab
17、优化器将会通过更少的迭代次数来解决标准的拓扑优化问题。4.5 其他拓展扩展到三维空间应该是直截了当的，然而更加复杂问题比如机械设计需要MMA优化器的运用以及定义更多地变量。Matlab指令的简单化可以考虑一些图形输出，交互式输入的简单拓展。5 结论这篇文章展示了基于拓扑优化算法的简单的数学理论实现。这串代码仅有99行输入实现，包括优化准则，独立网格过滤器和有限元分析。这串matlab代码可以在网页上下载以作教育使用。同时这个代码可以很容易的拓展到多重负载问题和无源区域的定义问题。在Matlab运行代码相当缓慢相比于在网站上用Fortran语言实现相同的代码。然而，matlab上一些硬件设备可支
18、持更加有效的代码生成并优化运行时间。值得强调的是运行速度可以通过调整代码来实现，也以代码的简单性为代价。附录 99行代码1 % A 99 LINE TOPOLOGY OPTIMIZATION CODE BY OLESIGMUND, OCTOBER 1999 %2 function top(nelx,nely,volfrac,penal,rmin);3 % INITIALIZE4 x(1:nely,1:nelx) = volfrac;5 loop = 0;6 change = 1.;7 % START ITERATION8 while change 0.019 loop = loop + 1;1
19、0 xold = x;11 % FE-ANALYSIS12 U=FE(nelx,nely,x,penal);13 % OBJECTIVE FUNCTION AND SENSITIVITY ANALYSIS14 KE = lk;15 c = 0.;16 for ely = 1:nely17 for elx = 1:nelx18 n1 = (nely+1)*(elx-1)+ely;19 n2 = (nely+1)* elx +ely;20 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);21 c = c + x(el
20、y,elx)penal*Ue*KE*Ue;22 dc(ely,elx) = -penal*x(ely,elx)(penal-1)*Ue*KE*Ue;23 end24 end25 % FILTERING OF SENSITIVITIES26 dc = check(nelx,nely,rmin,x,dc);27 % DESIGN UPDATE BY THE OPTIMALITY CRITERIA METHOD28 x = OC(nelx,nely,x,volfrac,dc);29 % PRINT RESULTS30 change = max(max(abs(x-xold);31 disp( It.
21、: sprintf(%4i,loop) Obj.: sprintf(%10.4f,c) .32 Vol.: sprintf(%6.3f,sum(sum(x)/(nelx*nely) .33 ch.: sprintf(%6.3f,change )34 % PLOT DENSITIES35 colormap(gray); imagesc(-x); axis equal; axistight; axis off;pause(1e-6);36 end37 % OPTIMALITY CRITERIA UPDATE %38 function xnew=OC(nelx,nely,x,volfrac,dc)3
22、9 l1 = 0; l2 = ; move = 0.2;40 while (l2-l1 1e-4)41 lmid = 0.5*(l2+l1);42 xnew = max(0.001,max(x-move,min(1.,min(x+move,x.*sqrt(-dc./lmid);43 if sum(sum(xnew) - volfrac*nelx*nely 0;44 l1 = lmid;45 else46 l2 = lmid;47 end48 end49 % MESH-INDEPENDENCY FILTER %50 function dcn=check(nelx,nely,rmin,x,dc)5
23、1 dcn=zeros(nely,nelx);52 for i = 1:nelx53 for j = 1:nely54 sum=0.0;55 for k = max(i-round(rmin),1):min(i+round(rmin),nelx)56 for l = max(j-round(rmin),1):min(j+round(rmin), nely)57 fac = rmin-sqrt(i-k)2+(j-l)2);58 sum = sum+max(0,fac);59 dcn(j,i) = dcn(j,i) + max(0,fac)*x(l,k)*dc(l,k);60 end61 end6
24、2 dcn(j,i) = dcn(j,i)/(x(j,i)*sum);63 end64 end65 % FE-ANALYSIS %66 function U=FE(nelx,nely,x,penal)67 KE = lk;68 K = sparse(2*(nelx+1)*(nely+1), 2*(nelx+1)*(nely+1);69 F = sparse(2*(nely+1)*(nelx+1),1); U =sparse(2*(nely+1)*(nelx+1),1);70 for ely = 1:nely71 for elx = 1:nelx72 n1 = (nely+1)*(elx-1)+
25、ely;73 n2 = (nely+1)* elx +ely;74 edof = 2*n1-1; 2*n1; 2*n2-1; 2*n2; 2*n2+1;2*n2+2;2*n1+1; 2*n1+2;75 K(edof,edof) = K(edof,edof) +x(ely,elx)penal*KE;76 end77 end78 % DEFINE LOADSAND SUPPORTS(HALF MBB-BEAM)79 F(2,1) = -1;80 fixeddofs = union(1:2:2*(nely+1),2*(nelx+1)*(nely+1);81 alldofs = 1:2*(nely+1
26、)*(nelx+1);82 freedofs = setdiff(alldofs,fixeddofs);83 % SOLVING84 U(freedofs,:) = K(freedofs,freedofs) F(freedofs,:);85 U(fixeddofs,:)= 0;86 % ELEMENT STIFFNESS MATRIX %87 function KE=lk88 E = 1.;89 nu = 0.3;90 k= 1/2-nu/6 1/8+nu/8 -1/4-nu/12 -1/8+3*nu/8 .91 -1/4+nu/12 -1/8-nu/8 nu/6 1/8-3*nu/8;92
27、KE = E/(1-nu2)* k(1) k(2) k(3) k(4) k(5) k(6) k(7) k(8)93 k(2) k(1) k(8) k(7) k(6) k(5) k(4) k(3)94 k(3) k(8) k(1) k(6) k(7) k(4) k(5) k(2)95 k(4) k(7) k(6) k(1) k(8) k(3) k(2) k(5)96 k(5) k(6) k(7) k(8) k(1) k(2) k(3) k(4)97 k(6) k(5) k(4) k(3) k(2) k(1) k(8) k(7)98 k(7) k(4) k(5) k(2) k(3) k(8) k(1) k(6)99 k(8) k(3) k(2) k(5) k(4) k(7) k(6) k(1);专心-专注-专业