英语原文共 13 页，剩余内容已隐藏，支付完成后下载完整资料
基于Matlab的利用拓扑优化和能量均化的材料设计方法
摘要：本文针对具有极端性质的材料，提出了一种用于拓扑优化设计的Matlab代码。出于对代码简洁程度的考虑，本文采用基于能量均匀化的渐近逼近法，利用单元互能求出有效参数，求出相应的具有周期边界条件的解。通过约束材料的体积分数，程序可以取均质刚度张量目标函数的最大或最小值，如体积模量、剪切模量和泊松比。整个Matlab代码共88行，于附录给出（出自于Andreassen等人和结构优化的进展43(1):1–16, 2011）
关键词：拓扑优化、微观结构、均匀化、周期性边界条件、MATLAB
符号命名
宏观与微观尺度的纵横比
数值阻尼系数
拉格朗日乘子
单元密度设计变量
更新单元密度变量
周期性波动应变场
单元测试应变场
规定的应变场
附加应变场
体积分数的上界
从最优条件得到的项
目标函数
空间维度
在指数表示法中均化弹性张量
矩阵表示法中的均化弹性张量
固体材料的杨氏模量
一单元的杨氏模量
指标符号的弹性张量
孔隙材料杨氏模量(人造材料)
m 设计变量的位移限制
N 若干有限单元
p 惩罚因子
单元间的相互能量e
单元相互能量总和
U 微尺度位移场
微尺度周期波动场
取决于的宏观尺度位移场
微型Y-周期容许位移场
单位体积
对向节点规定的周期位移
X 宏观尺度笛卡尔坐标
Y 基本单元域
y 微尺度笛卡尔坐标
在j方向上基底细胞大小
W 单位元内的周期性移位
F 外力向量
K 整体刚度矩阵
单元刚度矩阵模量
单元刚度矩阵
U 全局位移矢量
荷载箱i j的单元位移矢量
1.前言
拓扑优化（Bendsoslash;e and Kikuchi 1988）最初是被Sigmund（1994）利用逆均化法应用于材料设计。接着是一系列系统工程（如Sigmund and Torquato 1997; Sigmund 2000; Gibiansky and Sigmund 2000）。这个问题后来被基于密度的方法（如Neves等人. 2000; Guest and Preacute;vost 2007; Zhang等人. 2007）、水平集法（如Challis 等人. 2008; Wang等人. 2014）、拓扑导数法（如Amstutz等人. 2010）和ESO型法（如Huang 等人. 2011）相继解决了。图1，2，3显示了一些典型拓扑优化设计中的微观结构。梯度材料和结构设计已由Paulino（2009）和Almeida（2010）等人给出。另一个密切相关的研究领域是并行材料和结构设计（如Rodrigues等人2002；Zhang和Sun2006；Xia和Breitkopf 2014和2015）。
继Sigmund（2001）在文章中开创性地给出99行MATLAB代码之后，一系列通过MATLAB实现的科研论文极大地促进了拓扑优化的普及和发展。这些包括通过利用耦合水平集方法并使用FEMLAB软件的Liu等人（2005），使用ESO法的Huang和Xie（2010）、使用离散水平集方法的Challis（2010），利用拓扑导数求帕累托最优跟踪提出199行代码的Suresh（2010），由Andreassen（2011）等人提出的88行MATLAB代码。针对多边形网格生成的matlab代码（polymesher）和基于拓扑优化框架（PolyTop）的Talischi（2012a，b），还有实现并行计算的（Mahdavi等人。2006）。
当今作者也从这些科研论文中受益匪浅，例如，多组分结构体系（Xia 等人. 2012, 2013）就是在99行代码（Sigmund 2001）的框架下给出的。简化的多尺度拓扑优化（Xia 和 Breitkopf 2014b）就是采用离散水平集方法（Challis 2010）。此外，在并行材料和结构设计（Xia and Breitkopf 2014a, 2015）的研究上，作者们都基于88行代码框架（Andreassen 等人2011）和ESO的优化器（Huang and Xie 2010）。
本文将88行代码运用到具有极端性质材料的优化设计上。我们遵循Sigmund（1994）提出的设计原则，其中均质材料的本构参数是根据单元相互能量来计算的。其是对有效材料性能的一种预测，我们采用基于平均应力和应变定理（Hashin 1983）的能量等效均匀化方法，而不是使用传统的渐近展开法（Guedes和Kikuchi 1990）。它将在第6节中展示，其应用程序设计算法与MATLAB程序（见附录）可以产生具有极端性质的微观结构和具有类似性的拓扑结构，如图1, 2和3。
论文的其余部分安排如下：第2节简要介绍了均匀化理论。第3节介绍了周期边界条件的实现方法。第4节给出了优化模型。第5节解释了MATLAB的实现方法。第6节列举出了几个使用本文所提出代码的例子。第7节得出结论。MATLAB的安装和使用在附录中给出。
图1最大体积模量的微观结构：前两个来自Sigmund（2000），第三张来自Zhang等人.（2007），最后一张来自于Amstutz等人.（2010）（从左到右）
图2具有最大剪切模量的组织：前两个来自于Neves等人.（2000），第三个来自于Huang等人.（2011），最后一张来自于Amstutz等人.（2010）（从左到右）
图3最小负泊松比的组织：前两个来自于Amstutz等人.（2010），最后两个来自于Wang等人.（2014）（从左到右）
2.均匀化理论
在线性弹性范围内。利用均匀化方法可以评价微观结构周期性的本构表现(图4)。对于在中的单元Y
(1)
,,是在基本单元尺寸上的三个方向。
图4周期性模式下微观结构的物质组成
在渐进均匀化进程中，宏观位移场取决于宏观和微观之间的比列，表示为
(2)
其中所涉及的函数依赖于宏观变量x和局部微观变量y。对于这个关系，一个量在点X附近的一个很小的邻域内变化，它可以被看作是“延伸”。当时，对于固定点X，Y可视为具有周期性的，此时微观尺度就相当于宏观。
当仅仅考虑（2）渐进展开的一阶项，均匀化刚度张量是通过在基元y上平均积分而给出的，可表示为
这里利用了爱因斯坦指数，是下式Y周期的解
其中V是Y-周期容许位移场，相当于三（二维）或六（三维）线性独立的单元测试应变场。
由于本文提出了一个紧凑的matlab代码，因此采用了基于能量的平均应力和应变定理，而不是渐近方法。基于能量的方法直接规定单元测试菌群的基底边界，诱导对应的应变场（见3）。根据Hashin（1983），这些都是对材料有效性能预测的两个等价方法。第3节给出了周期边界条件的详细实现
为了支持在拓扑优化中使用有效的现有算法，（3）依照单元互能原理（Sigmund:1994）以相同的形式重新整理如下
在有限元分析中，基本单元被离散化为N个有限单元，（5）可以近似于
其中与单元测试应变场是对应的单元位移解，是单元刚度矩阵。在二维情况下，我们注意到11→1，22→2，和12→3，允许以扩展的形式写成（6）
其中可以表示为
是的总和
3周期性边界条件(PBC)
通过求解基本单元的平衡问题，对（5）中的应变场进行求值，在周期性的假设下，在给定的应变下，基本单元的位移场可以写成宏观位移场和微观周期波动场的总和（Michel等人. 1999）。
在实践中，（10）不能直接强加在边界上，因为周期性波动项是未知的。这种通用表达式需要转换成约束，该约束在基本单元对应曲面上的节点之间（Xia等人.2003）。如图5所示的二维基底单元，一对相反边界上的位移是
上面的“k ”和“k”表示两个相反的平行边界面，它们的方向垂直于k方向（k=1,2,3）。周期项可以通过位移之间的差来消除
对于任意基于基本单元给出的平行六面体，是不变的，在图5中，我们有，，，。因此由于的确定，等号右边成为确定的值
.
因为, 这种形式的边界条件可以通过约束相应的节点位移来直接施加在有限元模型中。同时，这种形式的边界条件满足了周期性和连续性要求,以及在使用基于位移的有限元分析时的应力要求(Xia等人. 2006)。
图5. 二维矩形基单元模型
4最优模型
4.1改进的SIMP方法
将基本单元离散为N个有限单元，并使用改进的SIMP方法(Sigmund2007)，相应地定义了相同数量的密度设计变量。单元杨氏模量E被定义为
为杨氏固体材料的弹性模量，为Ersatz材料的杨氏模量，这是一种利用柔性材料对孔隙材料的近似方法(Allaire.等人 2004)，从而防止刚度矩阵的奇异性。值在0和1之间，与这些限制分别对应的是替代的固体材料。p是一种惩罚因数，它使密度分布更接近所谓的黑白解决方案。
优化问题的数学公式如下
其中K为全局刚度矩阵，和分别为全局位移矢量和测试用(kl)外力矢量。d是空间维度，表示单元单元体积，thetasym;代表体积分数上限。表示均质刚度张量的目标函数。例如，在二维情况下，材料体积弹性模量的最大化对应于下式的最小值，可表示为
材料剪切模量最大化对应于下式最小值
4.2均化方程的数值解
当几何和加载都显示对称时，情况如第3节中提出的周期边界条件可以简化为常规边界条件(Hassani和Hinton 1998b)。为了保持推导过程的一般性，这种简化在目前的工作中是不适用的。相反，直接施加周期性边界条件(见第3节)。关于(13)的有限单元解，这里采用了直接消除冗余未知数的解决方案。需要注意的是，除了直接解决方案之外，还有另外两种使用惩罚方法和拉格朗日乘子方案(Michel 等人. 1999)。
全局位移分离向量U分为四个部分:表示规定的位移值,表示未知数相应的内部节点,和表示未知节点对应位于边界对面的基础单元，满足 = W ，W的值是根据(13)通过规定计算给定的。(15)的平衡方程可以扩展为
F1是一个未知的向量，其大小等于指定位移节点上的反力。由于周期性假设，F2=0，F3 F4=0。同时，在（18）中，K是对称的，i，e，。消除第一行，加上第三和第四行，利用关系 = W，（18）可以化简为
允许系统的解决方案。
4.3优化准则法
一旦得到位移解，优化问题(15)便可采用标准最优准则法求解。下面提出了一种启发式的修正方案（Bendsoslash;e and Sigmund 2003）
其中m是主动位移限制，是阻尼系数，通过如下最优近似条件得到(Bendsoslash;e and Sigmund 2003)
拉格朗日乘子lambda;的是通过最小二分法，通过限制材料体积分数来选择的。目标函数的灵敏度是用伴随性矩阵计算的(Bendsoslash;e and Sigmund 2003)
按照客观的定义，是一个杨氏模量对应的单位刚度矩阵。当使用均匀网络时，单位体积被设置为1，且=1。
为了确定（15）优化方案的存在性，灵敏度和密度的过滤方案将被用来避免棋盘格的形成和对网格的依赖性。
5 MATLAB的实现
在这一节中，将解释Matlab代码(参见附录)。当前代码是在88行代码基础之上构建的(Andreassen 等人. 2011)。保持前38条线不变。在第4行到第6行中定义材料属性。在第8至17行定义单元
全文共11314字，剩余内容已隐藏，支付完成后下载完整资料
资料编号：[11176]，资料为PDF文档或Word文档，PDF文档可免费转换为Word
以上是毕业论文外文翻译，课题毕业论文、任务书、文献综述、开题报告、程序设计、图纸设计等资料可联系客服协助查找。
