基于MATLAB和ANSYS的悬臂梁拓扑优化.docx
《基于MATLAB和ANSYS的悬臂梁拓扑优化.docx》由会员分享，可在线阅读，更多相关《基于MATLAB和ANSYS的悬臂梁拓扑优化.docx（29页珍藏版）》请在冰豆网上搜索。
基于MATLAB和ANSYS的悬臂梁拓扑优化
88-51-FB-5B-3A-66
现代设计理论和方法大作业
----基于MATLAB和ANSYS的悬臂梁拓扑优化
指导老师：
志刚
项目组长：
薛亚波
项目成员：
机自66
学院：
机械工程学院
基于MATLAB和ANSYS的悬臂梁拓扑优化
一、计划和任务安排表
1.1学习时间安排表：
任务
时间
查阅整理资料
5月03日——5月08日
学习相关知识和论文
5月09日——5月15日
设计实验方案
5月16日——5月17日
学习优化理论
5月18日——5月20日
撰写论文
5月21日——5月22日
集体交流学习心得
5月23日——5月24日
1.2任务分配表：
学号
任务
06011136
齐寰宇
查阅资料、答辩PPT制作
06011137
王晓庆
查阅资料、答辩PPT制作
06011139
冰
查阅资料、研究问题分析
06011141
康辰龙
查阅资料、ansys拓扑模块研究
06011143
梁伟云
查阅资料、项目背景研究
06011145
卢岗
查阅资料ansys拓扑模块研究
06011153
清伟
查阅资料ansys拓扑模块研究
06011154
薛亚波
查阅资料、后期工作总结
06011155
严宏伟
查阅资料ansys拓扑模块研究
06011156
志亮
查阅资料、项目背景研究
06011157
龙
查阅资料、答辩PPT制作
06012020
乔百杰
查阅资料、Matlab结构优化
06073012
邓磊
查阅资料、项目背景研究
06182047
立峰
查阅资料、研究问题分析
06055164
杜博
查阅资料、研究问题分析
06011146
牛玉辉
查阅资料、研究问题分析
06011147
朴群星
查阅资料、项目背景研究
二、项目背景介绍及问题描述
2.1项目背景及意义：
2.1.1工程背景及基本原理：
通常把结构优化按设计变量的类型划分成三个层次：
结构尺寸优化、形状优化和拓扑优化。
尺寸优化和形状优化已得到充分的发展，但它们存在着不能变更拓扑结构的缺陷。
在这样的背景下，人们开始研究拓扑优化。
拓扑优化的基本思想是将寻求结构的最优拓扑问题转化为在给定的设计区域寻求最优材料的分布问题。
寻求一个最佳的拓扑结构形式有两种基本的原理：
一种是退化原理，另一种是进化原理。
退化原理的基本思想是在优化前将结构所有可能杆单元或所有材料都加上，然后构造适当的优化模型，通过一定的优化方法逐步删减那些不必要的结构元素，直至最终得到一个最优化的拓扑结构形式。
进化原理的基本思想是把适者生存的生物进化论思想引入结构拓扑优化，它通过模拟适者生存、物竞天择、优胜劣汰等自然机理来获得最优的拓扑结构。
2.1.2本文研究意义：
目前，结构优化大部分集中在尺寸设计变量（如板厚、杆的剖面积及管梁的直径）。
拓扑结构优化较尺寸优化复杂,但对于有些问题拓扑结构优化比尺寸优化有效,悬臂梁是其中的例子之一。
本文讨论悬臂梁的拓扑优化问题,围绕这一问题,怎样使结构具有最大刚度的设计占有相当重要的地位;怎样优化结构的形状使材料的分布，更加合理从而达到使结构具有最大刚度的目的是本文要研究的问题。
2.2研究现状
2.2.1理论研究现状：
结构拓扑优化是近20年来从结构优化研究中派生出来的新分支，它在计算结构力学中已经被认为是最富挑战性的一类研究工作。
目前有关结构拓扑优化的工程应用研究还很不成熟，在国外处在发展的初期，尤其在国尚属于起步阶段。
1904年Michell在桁架理论中首次提出了拓扑优化的概念。
自1964年Dorn等人提出基结构法，将数值方法引入拓扑优化领域，拓扑优化研究开始活跃。
20世纪80年代初，程耿东和N.Olhoff在弹性板的最优厚度分布研究中首次将最优拓扑问题转化为尺寸优化问题，他们开创性的工作引起了众多学者的研究兴趣。
1988年Bendsoe和Kikuchi发表的基于均匀化理论的结构拓扑优化设计，开创了连续体结构拓扑优化设计研究的新局面。
1993年XieYM和StevenGP提出了渐进结构优化法。
1999年Bendsoe和Sigmund证实了变密度法物理意义的存在性。
2002年罗鹰等提出三角网格进化法，该方法在优化过程中实现了退化和进化的统一，提高了优化效率。
2.2.2应用研究现状：
在前人提出的重要理论基础上，后人也将其跟其他现代设计的方法相结合，衍生出了其他一些拓扑结构优化方法：
如与可靠性相结合的情况下，MAUTE等应用变密度法并结合可靠性分析对一微机电系统进行了基于可靠性的拓扑优化设计，PAPADRAKAKIS等将遗传算法应用于具有可靠性约束的桁架结构拓扑优化设计中，国学者马洪波也对基于遗传算法的结构可靠性优化问题进行了讨论。
华南理工大学机械工程学院欧阳高飞等对基于水平集方法的结构可靠性拓扑优化进行了研究。
2.3研究目标：
2.3.1设计目标：
设计一悬臂梁在自由端受到竖直向下集中力P的作用，要求保持悬臂梁长度尺寸不变。
2.3.2优化目标：
在满足以上条件情况下对原实体悬臂梁做结构拓扑优化设计，使其结构刚度最大、重量最轻。
（优化的结果应该使原设计区域产生孔洞，使结构拓扑发生变化。
）
优化目标分析：
要使梁满足以上条件情况必须同时兼顾重量和刚度问题。
重量太轻不能满足刚度问题。
刚度太大又必须足够的材料以满足其要求。
将两个因素进行耦合分析，我们可以得出的结论是必须通过结构层次来进行优化。
连续体结构拓扑优化较成熟的方法主要有均匀化方法、变密度方法和渐进结构优化法，通过综合比较我们最终选用变密度法。
2.3.3实现方法：
本文在对题目进行分析后，决定从两个独立的方向来分别研究，以此来互相验证，保证结果的可靠性。
实验中要使用ansys和matlab两种软件分别独立进行优化分析，然后将两种结果来进行比较。
2.3.4实践目标：
通过这次的作业我们期望通过搜索资料，团体讨论，分组作业的方式，以期望达到对拓扑优化的方法和过程有个具体的了解和学习。
另外对分析软件的应用能够达到一个新的高度。
这些不仅能使我们现在的知识体系得到充实和优化，而且也是我们今后人生的财富
2.4研究容
2.4.1问题描述
如图所示，悬臂梁在自由端受到竖直向下集中力P的作用，要求在保持悬臂梁长度尺寸不变的情况下对原实体悬臂梁做结构拓扑优化设计，优化目标是使结构刚度最大、重量最轻。
（优化的结果应该使原设计区域产生孔洞，使结构拓扑发生变化。
）
原实体悬臂梁横截面为矩形，矩形宽度为0.2m，悬臂梁的材料为45钢，密度ρ为7.8×10-6kg/mm3，弹性模量E=2×105MPa，许用切应力[τ]=60MPa。
竖直向下集中力为
。
图悬臂梁受集中力作用（单位：
m）
2.4.2问题转化
本研究中，要求在保持悬臂梁在结构尺寸不变的情况下，使悬臂梁的刚度最大、质量最轻。
一般悬臂梁都为三维立体结构，载荷是在竖直方向上加载的，所以在垂直纸面方向上，并且在悬臂梁的对称平面上加载，根据材料力学知识可知，所加载在和只对加载区周围较小围产生影响，其余部分与加载形式无关。
其受力与如图所示平面的手里相同，故可以将问题抽象为平面问题来处理。
当占有区域S的一个结构上作用有强度为P（S）的载荷时,使结构具有最大刚度的问题和使结构所受外力功W具有最小值的问题是等价的。
结构所受的外力功W与结构变形能C之间的关系为W=2C，所以我们可以将该问题处理为求在一定的约束条件下，该悬臂梁的的最小变形能最小的问题。
2.4.3数学模型：
该问题中，要求同时满足刚度最大，质量最轻，这两个变量若同时改变，则问题复杂度太大，并亲切可能导致问题不可求解。
所以我们采用在确定的质量下，来讨论刚度最大的问题。
由于对特定的材料，其质量和体积有一定的关系，并且我们采用去除法的思想来建立模型的，故我们可以采用给优化后的体积与优化前的体积比赋确定的值，来达到在给定质量条件下满足刚度最大的问题。
其数学模型如下：
注：
其中C（x）为结构变形能，U为结构变形总位移矩阵，K为结构总刚度矩阵，N为划分单元总数，
为单元位移向量，
为单元刚度，（由于划分单元的时候，我们采用等分举行单元，所以个单元的刚度可用一个常量来处理）V（x）是拓扑结构优化过程中变化着的体积，
为未经过优化前悬臂梁的体积。
F为结构所受的总载荷。
x为悬臂梁的相对密度。
2.4.3模型分析求解：
该问题的优化方法有很多种，常用的有如下方法：
OptimalityCriteria（OC）methods,（优化准则方法）
SequentialLinearProgramming（SLP）methods（序列线性规划法）
MethodofMovingAsymptotes（MMAbySvanberg1987）等
为了简化问题的复杂度，此处我们采用standardOC-method.方法来实现。
在处理过程中，关于设计变量相对密度x每一步的更新，我们采用
在1995年提出的如下算法来实现：
注：
其中m是移动限制量，即相对密度变化的最大步长，更新条件中，加上最小边界条件判断，是防止相对密度变为0，使问题出现奇异；加上相对密度最大边界条件判断，是为了防止出现相对密度大于1，是问题是去物理意义。
是数值阻尼系数。
可以由以下条件来确定：
其中
为拉格朗日乘数，由双向切片算法来确定。
目标函数对单元相对密度的变化率，即单元敏感度如下确定：
为了确保该问题的拓扑结构优化解的存在，在处理技术上，我们得对该问题的求解过程加些限制。
此处我们采用过滤技术来进行处理。
我们采用网格独立性滤波器来对单元敏感度的计算来进行改进，其改进算法如下：
其中
为卷积因子（即权重因子），其由以下表达式来确定：
采用该算法对单元敏感度的更新进行改进，是为了在更新单元e的敏感度时，将距其距离不超过
的临域的单元来进行均匀化，防止各乡邻单元的敏感度变化过大。
注：
其中dist（e,f）定义为单元体e的中心到单元体f中心之间的距离，
为过滤大小。
在优化过程中，我们采用改进的敏感度表达式（5）来代替表达式（4）来对（3）是进行计算。
三、实验设计及结果分析：
3.1技术路线：
在本实验中，由于我们提出了两种独立的方法，故我们将小组分为两部分;
小组一：
该组成员采用Ansys软件来对该悬臂梁模型进行建模分析，并进行拓扑结构优化。
小组二：
该组成员采用上面的数学模型，用Matlab原来来实现整个数学模型的算法。
并对处理结果进行动态显示，以便让我们清晰地看到整个拓扑结构优化的过程。
3.2试验结果及分析：
3.2.1Ansys分析的拓扑结构优化结果：
3.2.1.2建立几何模型：
3.2.1.2.1几何参数设置如下：
3.2.1.2.1在Ansys/workbench中建立如下几何模型：
3.2.1.2网格化：
将几何模型直接调入有限元环境。
在workbench里建立的模型可以无缝导入到有限元分析环境里。
采用自动网格划分工具，得到的模型如下：
其中，网格单元和节点数设置如下：
Statistics
Nodes节点数目
3120
Elements单元数目
561
3.2.1.3施加载荷和设置边界条件
选择优化类型为“shapeoptimization”——形状优化
在上截面施加15000N的力
3.2.1.4求解：
3.2.1.4.1设置材料属性：
Structural
Young'sModulus氏模量
2.e+011Pa
Poisson'sRatio
0.3
Density密度
7850.kg/m³
ThermalExpansion
1.2e-0051/°C
TensileYieldStrength
2.5e+008Pa
CompressiveYieldStrength
2.5e+008Pa
TensileUltimateStrength
4.6e+008Pa
CompressiveUltimateStrength
0.Pa
Thermal
ThermalConductivity
60.5W/m·°C
SpecificHeat
434.J/kg·°C
Electromagnetics
RelativePermeability
10000
Resistivity
1.7e-007Ohm·m
3.2.1.4.2定义期望目标值的减少量，即期望减少的质量或体积分数：
Scope
Geometry
AllBodies
Definition
TargetReduction
50.%
3.2.1.4.3拓扑结果如下：
3.2.1.4.4其他两种结果：
定义质量减少率为80%
Scope
Geometry
AllBodies
Definition
TargetReduction
80.%
拓扑结果：
定义质量减少率为60%
Scope
Geometry
AllBodies
Definition
TargetReduction
60.%
60%的拓扑结果如下
3.2.1.4试验结果分析及结论：
从拓扑结构我们可以看到，在不同的材料去除率下，总的趋势都是边界面梁，其中包含着等强度设计的思想，这种拓扑结果较为符合现实情况，在工程上具有很强的实用性。
在不同的材料去除率下，我们发现结果具有不同的拓扑结构。
3.2.2拓扑结构优化：
注：
以下结果是在材料去除率为40%的情况下进行的优化。
拓扑结构优化前，该结构的各项参数如下：
3.2.2.1优化分析：
在上面Ansys图中，我们可以发现，改图朴结构十分不规则。
工程实际中，这样的结构在加工上很难实现。
基于此，我们决定对上面的结果进行优化，得到较为规则的拓扑结构，以便工程上的加工。
3.2.2.1有限元分析：
优化前材料去除率为34%
Box
LengthX
1.6m
LengthY
1.m
LengthZ
0.2m
Properties
Volume
0.21212m³
Mass
1665.1kg
OriginalMass
2512kg
CADParameters
ds_d
520
ds_s
150
ds_D1
320
ds_ss
250
ds_dd
990
总位移变形图
应变图
应力图
优化后的分析结果，
ObjectName
TotalDeformation
EquivalentElasticStrain
EquivalentStress
Minimum
0.m
5.4554e-008m/m
10911Pa
Maximum
5.1147e-006m
3.119e-006m/m
6.2381e+005Pa
3.2.2.1优化设计：
根据拓扑结果，我们设置了5个参数作为优化设计变量，其各参数取值如下：
我们采用多目标优化，将质量和结构最大变形作为目标函数，使其取值最小。
在将多目标转化为单目标时，我们对质量和最大变形分别赋予了不同的权重，具体设置如下表：
优化结果：
通过对比以上三种结果，我们发现结果A较为理想。
Ratingdefinitions:
∙Themostpositiveratingof
meansthatthedesignisexcellentintermsofsatisfyingthegoaldefinedfortheparameter.
∙Themostnegativeratingof
meansthatthedesignispoorintermsofsatisfyingthegoal.
∙Aratingof
isneutralwithrespecttothegoal.
响应面:
3.2.3Matlab语言实现的拓扑结构优化结果：
3.2.2.1重要参数设计：
Nelx=32x方向单元的数目为32
nely=20y方向单元的数目为20
注：
该悬臂梁的长宽比为16：
10
f/volfrac=0.4保留原材料的体积分数为0.4
penal=3.0抑制权值为3.0（该取值是资料建议的典型值）
=1.2过滤大小为1.2（该取值是资料建议的典型值）
3.2.2.2优化实验结果：
3.2.2.3试验结果分析
从该实验结果来看，在我们给定的体积保留率的情况下，每经过一次拓扑结构优化，该优化程序就将悬臂梁的拓扑结构中强度要求不高处材料的密度减小，直到所有无用的材料都将被去除为止。
我们的拓扑结构优化模型是建立在结构变形能最小、体积去除率自己给定的基础上进行的，故我们可以根据实际情况，自行确定体积去除率。
在拓扑优化的过程中，我们可以观察到，我无论体积压缩率如何变化，悬臂梁模型最终都向桁架结构进化。
这说明，在结构件中，在自身材料多少相同的条件下，桁架具有很高的刚度和强度，其实这也就是为什么拓扑结构优化首先在桁架结构领域提出。
故工程上，我们常见工程人员采用桁架结构来作为一些工程的支撑结构，如塔吊等。
在实验中，在验证不同的体积压缩率时，我们发现，在不同的给定体积压缩率下，算法的有效性也不同，但在验证过程中，发现算法一直会收敛。
体积压缩率小的时候，该算法能很快终止；体积压缩率较大的时候，该算法的收敛速度较慢，并且还会出现不同程度的震荡，并且体积压缩率越大，该算法的振动也震荡。
3.3实验结论及总结：
实验中，为了充分发挥我们小组人员的优势，我们从两条独立的途径来分别对该问题进行研究。
Ansys中，我们利用其自带的的模块对该问题进行了建模，拓扑结构优化。
为了工程上的加工方便，我们在拓扑结构优化的基础上，对其形状进行了规则化，与此同时，进行再次优化。
Matlab中，我们将建立的数学模型用数值模拟的形式进行了处理，并将所得到的处理结果进行了图像动态显示，以此来清晰的观察拓扑结构优化的动态过程，给人以直观的印象。
在用Ansys进行拓扑结构优化的时候，我们发现，当材料去除率为60%时，其所得到的拓扑结构与我们用Matlab进行拓扑结构优化是所得到的结果的拓扑结构是一致的。
这验证了我们的数学模型是对的。
但是我们的模型也存在着一定的问题，即在材料去除率很小的时候，预加载荷未变的情况下，该模型的解均存在。
出现这一情况的原因是，我们在计算悬臂梁的位移矩阵，一直认为其应力未超过破坏强度，整个过程均在弹性变形围。
故在实际应用的时候，需要我们结合一定的工程经验来确定个重要参数的取值。
四、参考文献及附录：
4.1主要参考文献：
1.Sigmund,O.1994:
Designofmaterialstructuresusingtopology
optimization.Ph.D.Thesis,DepartmentofSolidMechanics,
TechnicalUniversityofDenmark
2.Sigmund,O.1997:
Onthedesignofcompliantmechanisms
usingtopologyoptimization.Mech.Struct.Mach.25,495–526
3.《ANSYS/WORKBENCH设计、仿真与优化》兵何正嘉雪峰等
清华大学2008年8月第一版
4..topopt.dtu.dk.
5.《基于结构最大刚度的形状优化方法》欧阳高飞建民
机械工程学报2008年10月
4.2相关附件：
附件1：
悬臂梁拓扑结构优化Matlab程序：
%%%%A99LINETOPOLOGYOPTIMIZATIONCODEBYOLESIGMUND,OCTOBER1999%%%
functiontop（nelx,nely,volfrac,penal,rmin）;%主程序，即调用程序。
其中的参数需要用户提供。
依次为x方向单元数、y方向单元数、期望优化的体积或质量百分比、抑制权值、过滤因子
%INITIALIZE初始化
x（1:
nely,1:
nelx）=volfrac;%在整个求解域均分材料
loop=0;
change=1.;%新旧改变之差
%STARTITERATION
whilechange>0.01%第一层，设置图像显示次数，即优化的次数
loop=loop+1;
xold=x;
%FE-ANALYSIS有限元
[U]=FE（nelx,nely,x,penal）;%调用有限元子程序，返回位移矢量
%OBJECTIVEFUNCTIONANDSENSITIVITYANALYSIS
[KE]=lk;%调用单元刚度阵
c=0.;
forely=1:
nely%16-24行一次循环。
外层循环，y方向循环nely次数
forelx=1:
nelx%层循环，x方向循环次数为nelx
n1=（nely+1）*（elx-1）+ely;%在全局坐标系中，左上方的节点数目
n2=（nely+1）*elx+ely;%在全局坐标系中，右上方的节点数目%n1、n2作用是从全局位移矢量提取局部或单元位移矢量
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）;%返回单元位移
c=c+x（ely,elx）^penal*Ue'*KE*Ue;
dc（ely,elx）=-penal*x（ely,elx）^（penal-1）*Ue'*KE*Ue;
end
end
%FILTERINGOFSENSITIVITIES灵敏度检验
[dc]=check（nelx,nely,rmin,x,dc）;%调用网格滤波子函数
%DESIGNUPDATEBYTHEOPTIMALITYCRITERIAMETHOD标准优化
[x]=OC（nelx,nely,x,volfrac,dc）;%调用优化子函数
%PRINTRESULTS拓扑过程显示
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）]）%优化过程指标显示
%PLOTDENSITIES密度分布
colormap（gray）;imagesc（-x）;axisequal;axis
升级会员 