《拓扑优化99行代码翻译》由会员分享，可在线阅读，更多相关《拓扑优化99行代码翻译（13页珍藏版）》请在金锄头文库上搜索。
1、拓扑优化中的拓扑优化中的 99 行行 matlab 代码代码 o.sigmund摘要摘要 这篇文章描述了用 matlab 语言来简洁的实现在静态负载下符合最小化原理的拓扑优化。总共只需要输入 99 行代码，包括优化程序和有限元分析子程序。这 99 行代码中，其中 36 行为主程序，12 行为基于最优控制器的优化程序，16 行为敏度过滤分析，其余 35 行代码作为有限元分析。实际上，除去注释行以及输出行、有限元分析行，仅有 49 行 matlab 代码输入用于解决一个适定的拓扑优化问题。再加上 3 行补充程序代码，这个程序就可以解决多种负载工况问题。这个代码主要是以教育指导为用，完整的 matl
2、ab 代码在附录中给出，同时也可以在网页 http:/www.topopt.dtu.dk 上下载。关键词 拓扑优化 教育 最优准则 万维网 matlab 代码1 简介简介文中展示的 matlab 代码主要是为工程教育所用。学生和在拓扑领域的新手可以在网页 http:/www.topopt.dtu.dk 上下载。这个代码可以用于结构最优化课程学习，学生们可以在多重负载工况、独立网格选择策略、无源场进行扩展应用。另一种可能就是用来激发学生们的直觉来进行最优化设计。研究生可以推测探究在给定边界条件和容量的情况下的拓扑优化并、比较得出最优策略。在文献中，你可以找到很多处理拓扑优化问题的方法。在一篇 B
3、endse and Kikuchi (1988)的原创论文中，基于对现存结论的学习，所谓微观结构或均化作用的方法被使用。均化作用方法在很多文章中都被采用，但它也存在一些缺点，比如对微观结构最优化方法果断的评估与决策很麻烦的，而且结果很难获得如果没有对微观结构进行确定的长度衡量。然而，在这个意义上来说均化作用方法对拓扑优化也是很重要的，它可以在结构的理论分析上提供一定的界限。拓扑优化的另一种方法叫做“幂律法则”或者 SIMP 法（(Solid Isotropic Material with Penalization) (Bendse1989; Zhou and Rozvany 1991; Mle
4、jnek 1992)） 。这里，假设物质性能使恒定不变的同时每个元素是设计区域离散化，变量是元素的相对密度。物质属性在相对物质密度增加到固体材料的物质属性的很多倍时被模板化。这种方法曾一度引起争议因为它认为没有任何物理材料的属性特征能被幂律法则所描述。然而，Bendse and Sigmund (1999)最近发表的文章证实只要在单一条件下能量能够满足，幂律法则在物理上就是可行的。为了确认这个结论的存在性，幂律法则必须与周界约束、梯度约束或者过滤技术相结合。这个拓扑优化的幂律法则方法已经应用于多重约束，多属性，多材料的问题中了。然而，上面提及到的解决方法是基于数学规划法和连续设计变量法，一系列
5、的文章都有涉及到解决拓扑优化的整数问题。Beckers (1999)通过双重途径成功的解决了大规模服从最小化问题，但是其他方法大多基于遗传算法或者其他为了几个元素需要成千上万的函数求值半随机方法，当然这很可能是不切实际的。除了上述提到的方法外，他们都能解决目标明确的问题，一些减少服从或者目标函数的启发式或直觉的方法也已经被提出来。这些理论都被统称为进化设计理论。除了很容易的理解和运用外，进化性分析方法主要的动机似乎是基于数学或连续变量法“涉及一些微积分应用和数学分析” ，并且他们包含“一些复杂问题的数学理论” ，反之进化性的方法能很好的应用强大的计算处理技术和自然进化过程中的直觉理念。两件事可
6、以反对他。第一，曾经由于更多的比服从最小化理论复杂的目标被考虑进去。进化论方法自身已经变得非常复杂。第二，正如文中所说，以数学理论为基础的方法解决服从最小化问题很容易实现并且计算处理同样很有效率。不仅如此，基于数学理论的方法很容易扩展解决像非共轭和多物理量的无服从问题以及多约束问题。而用扩展进化方法来处理这些问题似乎更加不可行。完整的 matlab 代码在附录中给出，文章剩余部分包括对优化问题的定义和讨论（第 2 部分） ，matlab 实现的详细解释（第 3 部分） ，接着是扩展问题的讨论（第 4 部分）和最后总结（第 5 部分） 。2 拓扑优化问题拓扑优化问题有许多简化方法都是用来简化 m
7、atlab 的代码。第一，假设设计区域是矩形且被平方有限元离散化。这种情况下，元素的数目和结点就很容易表示（一列一列从左上角开始）并且结构的纵横比通过水平（nelx）和垂直（nely）方向元素比率来确定。一个拓扑优化问题基于指数定律法，目标是实现最小化，可以如下表示（1）式中 U 和 F 分别表示整体变形和力的向量，K 表示整体刚度矩阵，和分别euek表示元素的位移矢量和刚度矩阵，x 是设计变量的向量，是相对密度的最小minx向量（非零以避免奇点）N（=nelx * nely）是设计区域离散化的元素数目，p 是惩罚因子（通常 p=3） ，和是材料体积和设计区域体积，f (volfrac)为规定
8、(x)V0V的容积率。该拓扑优化问题（1）可以用几种不同的方法来解决，如 Optimality Criteria（OC 算法） ，Sequential Linear Programming (SLP)，Method of Moving Asymptotes (MMA by Svanberg 1987)或者其他的方法。为简单起见，我们这里用标准的 OC 算法。根据 Bendse (1995)对设计变量应用启发式的调整法，公式表示为其中 m（移动量）是一个正的移动界限， (= 1/2)是数值阻尼系数，在最优eB化条件中表示为其中 是拉格朗日乘子，在 bi-sectioning algorithm
9、中获得。目标函数的敏感度被表示为更多的详细信息将通过最优化准则理论的推导和实现，在文献(e.g. Bendse 1995).中提到。为了确保对拓扑优化问题（1）解决方法的存在性，将会介绍一些作为结果的限制条件。这里我们使用过滤技术，需要强调的是这个过滤技术并不能证明方法的可行性，但是通过作者大量应用证明过滤技术在实际中能产生独立性网格。这个网格独立性滤波器是来修改元素敏度，如下:这个卷积算子（重量系数）表示为fH控制函数 dist(e, f)定义为元素中心 e 到元素中心 f 的距离。在过滤面积以外卷积算子为 0.卷积算子随着元素 f 的距离成线性减少。除了原始的敏度（4） ，修fH改的元素敏
10、度（5）被用于 OC 算法（3） 。3 matlab 算法实现算法实现Matlab 代码(见附录),建立了作为一个标准拓扑优化的代码。他的主程序为其中 nelx 和 nely 分别是横轴和纵轴方向上元素的个数，volfrac 是容积率，penal 是惩罚因子，rmin 是过滤半径。其他的变量如边界条件在 matlab 代码中定义且可根据实际需要进行调整。在拓扑优化每一次循环迭代中，程序都会生成一幅当前密度分布图来。图 1 显示了附录代码在给定以下输入条件下的结果密度分布默认的边界条件对应于“黑梁”的一半 (图 1)。垂直方向的负载承受在左上角，水平方向结构左边和右边为对称边界。Matlab 代
11、码中重要的细节将在下面各小节讨论。3.1 主程序（主程序（1-36 行）行）主程序的开始是材料设计区域的均匀分配（第 4 行） 。一些基本初始化后，主要的循环开始调用有限元分析子程序（12 行）并返回位移矢量 U。因为固体材料的元素刚度矩阵对任何元素来说都是一样的，元素刚度矩阵子程序仅被调用一次（14 行） 。接下来，一个循环遍历所有元素(第 16- 24 行)来确定目标函数和敏感性(4)。变量 n1 和 n2 表示在所有节点数中左上角和右上角的元素节点数目，并用于从整体位移矢量 U 中提取元素位移矢量 Ue。灵敏度分析将伴随着调用网格独立性过滤器(第 26 行)和最优准则优化器(28 行)。
12、当前的合规以及其它参数在 30-33 行中体现出来，并绘制出最后结果的密度分布图(第 35 行)。主循环会在改变值（改变值在 30 行确定）小于 1%是终止，否则以上步骤将重复执行。3.2 基于优化程序的最优化准则（基于优化程序的最优化准则（37-48 行）行）更新后的设计变量是通过优化程序找到的（37-48 行） 。由于材料体积(sum(sum(xnew)是一个拉格朗日乘子的单调递减函数，满足体积约束的拉格朗日乘子数值可以通过 bi-sectioning algorithm 算法得到（40-48 行） 。bi-sectioning algorithm 的初始化时通过假设一个小于 l1 和大于
13、 l2 结合的拉格朗日乘数（39行） 。限定拉格朗日乘数的间隔大小会一直减半直到它的大小小于收敛判别准则（40 行） 。3.3 网格独立性过滤技术（网格独立性过滤技术（49-64 行）行）Matlab 中第 49-64 行是为了实现条件（5） 。注意，并不是在设计域里所有的元素都要被搜索来寻找那些在处于过滤半径之内的元素，而是在那些被考虑元素周围，2 倍过滤半径大小的区域。通过选择过滤半径 rmin 小于调用子程序后的半径，被过滤后的敏度将会和初始状态一样，并使得过滤器不起作用。3.4 有限元分析代码（有限元分析代码（65-90 行）行）有限元分析代码位于 65-90 行中。注意解决者是利用
14、matlab 中的稀疏选项。整体刚度矩阵是由所有元素的循环建立的（70-77 行） 。正如主程序中，变量 n1和 n2 指示在所有节点数目中左上部分和右边节点数目并插入到整体刚度矩阵合适的位置中。正如前面所说，节点和元素都是从左到右纵向排布。而且，每个节点都有2 个自由度（水平和垂直的） ，因此指令 F(2,1)=-1.第 79 行应用了一个在左上角的垂直单元力。消除线性方程中的固定自由度来实现支承结构。Matlab 可以很好的处理，由这一行其中 freedofs 表示不受约束的自由度。通常来说，定义一个固定的自由度要简单的多，其后会被自动找到，通过使用 matlab 算子 setdiff 因
15、无约束自由度与固定自由度的不同来找到无约束自由度（82 行） 。元素的刚度矩阵在 86-99 行中进行计算。这个针对 4 个节点的 8*8 矩阵可通过使用人工智能软件进行处理分析。杨氏模量 E 和泊松比可以在 88 和 89 行进行改变。4 拓展拓展附录给出的 matlab 代码解决如图 1 中材料分配的优化问题以达到最小化原则。而算法中许多的拓展和改变都值得思考，其中一部分将在下面进行阐述。4.1 其他边界条件其他边界条件为解决其他优化问题，改变边界条件和支承条件非常容易。为了解决如图2 所示的短悬臂梁问题，只需要要 79 和 80 作相应改变，根据这些改变，输入行变为4.2 多重负载问题多
16、重负载问题扩展算法来处理多重负载问题也很容易。实际上，仅需额外增加 3 行，并对其他 4 行做微小调整即可。就两个负载工况而言，力和适量位移必须定义为两方向向量，这就意味着第 69 行改变为现在目标函数为二维之和因此第 20-22 行即被以下行所替代为解决如图 3 的二重负载问题，一个右上方的负载要添加到底 79 行，如下4.3 无源元件无源元件 在某些情况下,一些元素可能需要采取最低密度值(例如一个管上的洞)。一个 nelx*nely 的数组中任何元素都可为 0 和某些元素固定为 0，都可以定义到主程序中或者转移到 OC 算法中。添加行为在 OC 子程序中寻找无源元素并设置他们的密度为最小值（0.001） 。图 4 显示的是输入以下行生成的结构结果当这紧接着的 10 行添加到主程序（第 4 行后） ，便可在半径为 nely/3 的圆中找到无源元素，中心为（nely/2，nelx/3）4.4 可选择的优化器可选择的优化器不可否认,基于优化器实现的最优准则只适合单一约束并且它基于启发式定点更新方案。为了建立一个更好的优化方案，我们可以得到
