拓扑优化学习报告-北理工-王路

1、北京理工大学车辆工程王路基于99行程序的拓扑优化学习报告(-)背景和前言随着汽车工业的飞速发展以及日益突出的能源问题，汽车工业面临的挑战以及竞争环境也越來越激烈，对汽车产品提出了降低其制造成本及燃油经济性的新要求。在提高汽车安全性、减少汽车排放和解决能源消耗的背景下，提出了汽车轻量化技术。实现汽车轻量化的途径包括三个方面：结构优化技术、新型材料和先进性制造工艺。其中，我们所讨论的是结构优化技术，其中结构优化设计分为三个层次：尺寸优化(SizeOptimization)、形状优化(ShapeOptimization)和拓扑优化(TopologyOptimization)。本文我们基于99彳亍ma

2、tlab程序初步学习拓扑优化技术中的理论和优化方法。拓扑优化技术指的是在给定的设计空间内寻求最佳的材料布局，同时在满足平衡方程、物理关系、几何关系和边界约朿条件下使得结构达到某种性能最优的应用技术。拓扑优化的理论研究最早可以追溯到Michel提出的桁架理论，连续体结构的拓扑优化由于描述和数值计算得困难，发展一直相对缓慢，直到Bendsoe和Kikuclii在1988年提出的均匀化方法之后才得到迅速的发展，其基本思想是在组成拓扑结构中引入微结构，通过微结构的儿何参数作为设计变量，通过微结构的增加和删减实现结构的拓扑形状的改变，实现拓扑优化和尺寸优化的统一。在微结构的基础上，我们介绍变密度法的应用

3、，变密度法是在均匀化方法的基础上产生的，把材料引入微结构代之以密度在01之间变化的假想材料，把密度作为设计变量，从而实现材料的删减，因其模型简单、计算变量相对较少成为目前广泛采用的方法。根据不同的插值模式，变密度法乂有不同的插值模型：SIMP法(SolidIsotropicMaterialwithPenalization)、Hasliin-Shtrikinan法，以及RAMP法(RationalApproximationofMaterialProperties)o(二)拓扑优化问题的描述材料插值模型变密度理论的材料插值模型通过引入中间密度单元,将离散型问题转化成连续型优化问题，而实际上，中间密

4、度单元是无法存在和制造，因此要尽量避免中间密度单元的产生，减少中间密度单元的数目，此时就需要对设计变量中出现的中间密度值进行惩罚。目前材料插值模型主要分为两类，SIMP法和RAMP法，基于SIMP格式的材料插值模型为：EO0=+(X,)P(Eg-EQ&W0,1(1)其中：Eg)插值以后的弹性模量民一一实体部分材料的弹性模量孔洞部分材料的弹性模量耳一一单元相对密度，取值为1时表示有材料，为0时表示无材料即孔洞P一一惩罚因子在本文中，在本篇文章中，讨论的SIMP的表达为：E(勒)二随)卩耳(2)其中：先为第i个子域内第j个单元的相对密度。数学模型根据在体积或质量约束下求最小柔度，即最大刚度来建立优

5、化模型，基于变密度理论的SIMP法的周期性拓扑优化问题的数学模型可以表达为：findX=(X11,2,3,-.,j)TeRi=1,2,-、mJ=1,2,-、门mnmnminC(X)=FU=UtKU二工工代占,旳二丫(切川旅叫1=1J=11=1j=lS.tKU=Fmni=lj=l0.01loop=loop+1;xold=x;%FE-ANALYSISU=FE(nelx,nely,x5peiial);%OBJECTIVEFUNCTIONANDSENSITIVITYANALYSISKE=lk;c=0.;forely=l:nelyforelx=l:nelxnl=(nely+1)*(elx-1)+ely;

6、112=(nely+1)*elx+ely;Ue=U(2*nl-l;2*nl;2*n2-l;2*n2;2*n2+l;2*ii2+2;2*nl+l;2*nl+2,l);c=c4-x(ely,elx)Apenal*Uer*KE*Ue;dc(ely,elx)=-penal*x(ely,elx)A(penal-l)*Uef*KE*Ue;endend%FILTERINGOFSENSITIVITIESde=check(nelxelyjmin兀de);%DESIGNUPDATEBYTHEOPTIMALITYCRITERIAMETHODx=OC(nelx4iely,x,volfiac,dc);%PRINTRES

7、ULTSchange=max(max(abs(xxold);disp(It.:sprintfC%4i；loop)Obj.:sprii】tf(%10.4f,c).1Vol.:rspriiitf(F%6.3f,siim(siim(x)/(nelx*nely)1ch.:rsprintf(r%6.3f.change)%PLOTDENSITIEScolonnap(gray);imagesc(-x);axisequal;axistight;axisofif;pause(le-6);end%OPTIMALITYCRITERIAUPDATE%fiinctionxiiew=OC(nelx,nely,x,volf

8、iac,dc)11=0;12=100000;move=0.2;wliile(12-11le-4)lmid=0.5*(12+11);xiiew=max(0.001,max(xmove,niin(1.4iiin(x+move.*sqrt(-dc./lmid);ifsiim(siim(xiiew)volfrac*i】elx*nely0;=lmid;else=lmid;endend%MESH-INDEPENDENCYFILTER%fiinctiondcn=check(nelxj】elyjmin,x,dc)dcn=zeros(nelyjielx);fori=1:nelxforj=l:nelysiim=0

9、.0;fork=max(ifloor(rniin)J):niin(i+floor(rniin),nelx)fori=max(j-floor(rmin),1):niin(j+floor(rmin)4iely)fac=nnin-sqit(i-k)A2+(j-l)A2);sum=siim+max(0ac);dcn(j,i)=dcn(j,i)+max(0,fac)*x(lJc)*dc(l,k);endenddcn(jj)=dcii(jj)/(x(jendend%FEANALYSIS%fiinctionU=FE(nelx,nely,x,penal)KE=lk;K=sparse(2*(nelx+1)*(n

10、ely+1),2*(nelx+1)*(nely+1);F=sparse(2*(nely+1)*(nelx+1),1);U=zeros(2*(nely+1)*(nelx+1),1);forelx=l:nelxforely=l:nelynl=(nely+1)*(elx-1)+ely;112=(nely+1)*elx+ely;edof=2H：nl-l;2+nl;2*112-1;2*112;2*n2+l;2*n2+2;2*nl+l;2*111+2;K(edof,edof)=K(edof,edof)+x(ely,elx)ApenalM：KE;endend%DEFINELOADSANDSUPPORTS(H

11、ALFMBB-BEAM)F(2*(nelx/2+1)*(nely+1),1)=1;fixeddofs=2*(nely/2+1),2nelx*(nely+1)+2*(nely/2+1);alldofs=1:2*(nely+1)*(nelx+1);fieedofs=setdiff(alldofs,fixeddofs);%SOLVINGU(fieedofs,:)=K(fieedofs,freedofs)F(freedofs5：);U(fixeddofs,：)=0;%ELEMENTSTIFFNESSMATRIX%fiinctionKE=lkE=l.;nu=0.3;k=l/2-nu/6I/8+1111/

12、8-l/4-ini/12-l/8+35,:nu/8.-1/4+11U/12-l/8-mi/8nu/6l/8-3*nu/8;KE=E/(l-nuA2)*k(l)k(2)k(3)k(4)k(5)k(6)k(7)k(8)k(2)k(l)k(8)kkk(5)k(4)k(3)k(3)k(8)k(l)kkk(4)k(5)k(2)k(4)k(7)k(6)k(l)k(8)k(3)k(2)k(5)k(5)k(6)k(7)k(8)k(l)kkkk(6)k(5)k(4)k(3)k(2)k(l)k(8)k(7)k(7)k(4)k(5)k(2)k(3)k(8)k(l)k(6)k(8)k(3)k(2)k(5)k(4)kk

13、(6)k(l);%根据程序的运行的先后顺序，依次详细分析各行程序：fiinctiontop(nelx,nely,volfiac,peiialjmin);nelx=80;x轴方向的单元数目nely=20;y轴方向的单元数目volfiac=0.4;体积比，这个参数是用來确定保留的材料量penal=3;材料插值的惩罚因子，由于变密度理论的材料插值模型引入中间密度单元，所以将离散型问题转化成连续型优化问题，而在实际过程中，是无法存在和制造的，所以要尽量避免中间密度单元的产生，减少中间密度单元的数目，所以需要对中间密度值进行惩罚，penal是惩罚因子。niiin=2;灵敏度过滤的过滤半径关于体积比、惩罚

14、因子和过滤半径不同引起结果的不同，我们将在之后的讨论中进行。以上是初始化参数的程序段x(l:nely,l:nelx)=volfiac;x是设计变量，为某一个单元的相对密度，在此给出整体网格编排的顺序表：北京理工大学车辆工程王路北京理工大学车辆工程王路x(1tnelx)x(nelyj)xiney.nex：图2.整体网格编排顺序表故对于整个网格，在垂直方向上，有nely个单元，在水平方向上，有nelx个单元。loop=0;存放迭代次数的变量change=1.;每次迭代之后，目标函数的改变值，可以用來判断何时收敛wliilechange0.01当两次连续目标函数迭代的差小于0.01时，结束迭代loo

15、p=loop+1；迭代次数加1xold=x;将前一次的设计变量赋值给xoldU=FE(nelx,nely,penal);进彳亍有限元分析,需要输入的参数为nelx、nely、x、penal进入有限元分析子程序fiinctionU=FE(nelx,nely,x,penal)KE=lk;单元刚度矩阵分析进入单元刚度矩阵子程序fiinctionKE=lkE=l.;k=l/2-nu/61/8+11U/8-l/4-nii/12-l/8+3*nu/8.-1/4+11U/12-l/8-ini/8mi/6l/8-3*nu/8;KE=E/(l-nuA2)*k(l)k(2)k(3)k(4)k(5)k(6)k(7)

16、k(8)k(2)k(l)k(8)k(7)k(6)k(5)k(4)k(3)k(3)k(8)k(l)kkkk(5)k(2)k(4)k(7)k(6)k(l)k(8)k(3)k(2)k(5)k(5)k(6)k(7)k(8)k(l)k(2)kk(4)k(6)k(5)k(4)kk(2)k(l)k(8)k(7)k(7)k(4)k(5)k(2)k(3)k(8)k(l)k(6)k(8)kkk(5)kkkk(l);此处单独解释有限元理论，通过平面问题的4节点矩形单元描述（1）单元的几何和节点描述平面4节点矩形单元如图3所示,单元的节点位移共有8个自由度（DOF）,节点的编号为1、2、3、4,各自的位置坐标为Yi）

17、1=1,2,3,4,各个节点的位移（分别沿X和y方向），为（,vj.i=1,2,3,4u2u3节点4(x2,y2)节点3(x3,y3)v2ulv4u4节点4(x4,v4)节点1(xl,vl)北京理工大学车辆工程王路北京理工大学车辆工程王路图3平面4节点矩形单元将所有节点上的位移组成一个列矩阵，记作0,将所有节点上的力组成一个列矩阵，记作F,贝IJ：U)=比人U2v2u3v3比vj(3)北京理工大学车辆工程王路北京理工大学车辆工程王路北京理工大学车辆工程王路北京理工大学车辆工程王路(4)F(e)=ki%Fx2FygFx3F岁FX1Fy1（2）单元场的位移描述从图3可以看出，节点条件共右8个，*方

18、向上4个（Ui，U2，Us，u.i）y方向上4个（匕少2，仏，匕），故在乂和y方向的位移场可以各有4个待定系数，取以下多项式作为位移场模式：u(x,y)=a0+axx+a2y+a3xyy)=b+b】x+b2y+b3?y要考虑到的特点是当X或y不变时，沿y或X方向位移函数是线性变化(6)的，因此不选X2和才项。由于节点条件，在x=xi,y=yi处，有:=1,2,3,4将四个节点的坐标（,yx）,y2）,ys）,（x,y4）代入（5）方程中，综合（6）北京理工大学车辆工程王路北京理工大学车辆工程王路北京理工大学车辆工程王路北京理工大学车辆工程王路方程得:北京理工大学车辆工程王路北京理工大学车辆工程

19、王路4a_U1-U2-U3+U!d少24bux-u2+u3-u4L4ab经过整理，令(5)式的格式为:4abU1+u2+u3+u4TOC o 1-5 h z304a_U14-U2-U3-U1b5b_%+匕一冬一旳14ab叫-冬+耳北京理工大学车辆工程王路北京理工大学车辆工程王路北京理工大学车辆工程王路北京理工大学车辆工程王路(7)lzKo)zy-by-by-by-bu(x,y)=NJxyX+N2(x,y)u2+N3(x,y)u3+N,(x,y)u4i/(x,y)=N(x,y)+N2(x,y)v2+N3(x,y)v3+N(x,y)匕可以整理得：Nx(x,y)=j(l+-)(l+4axN2(x,y

20、)=-(1-)(1+4aN3(x,y)=l(l-)(1-4aN4(x,y)=|(l+-)(1-4a将（7）式写成矩阵形式:u(x,y)(x,y)北京理工大学车辆工程王路北京理工大学车辆工程王路N00NN20N300N20N3U1匕%匕U3=N-U(e)(2x8)(8x1)(9)北京理工大学车辆工程王路北京理工大学车辆工程王路北京理工大学车辆工程王路北京理工大学车辆工程王路其中y）为该单元的形状函数矩阵。（3）单元应变场的描述由弹性力学中节点的儿何方程:dvOudii%得到单元应变的表达:北京理工大学车辆工程王路北京理工大学车辆工程王路北京理工大学车辆工程王路北京理工大学车辆工程王路=du=dN

21、U(e)(3x2)(2x1)(3x2)(2x8)(8x1)BU(e)(3x8)(8x1)(10)北京理工大学车辆工程王路北京理工大学车辆工程王路北京理工大学车辆工程王路北京理工大学车辆工程王路O0一5yo一axO0一內-引X2)3VZ其中:北京理工大学车辆工程王路北京理工大学车辆工程王路北京理工大学车辆工程王路北京理工大学车辆工程王路单元的形状函数矩阵:O凡No3oN3NOon2n2ooN1N=(2x8)北京理工大学车辆工程王路北京理工大学车辆工程王路北京理工大学车辆工程王路北京理工大学车辆工程王路定义单元的儿何矩阵为：don4NooN3n3oon2n2ooNdx.B(x,y)=dN=0(3x

22、8)(3x2)d二耳D耳B(3x2)(3x2)(3x2)(3x2)j二123,4(12)式(11)中的子矩阵Bi为:B1二(3x2)（4）单元应力场的表达由弹性力学中平面问题的物理方程:北京理工大学车辆工程王路北京理工大学车辆工程王路北京理工大学车辆工程王路北京理工大学车辆工程王路心二吉(6-勺)E夠二云一6)y_2(1+“)心-rxy将(13)式换成应力和应变的关系表达为:q二(邑+夠)叨占7(%+心E故有:E=D(3x3)(3xl)其中定义D为弹性矩阵:(3x3)把(10)式代入到(15)其中应力函数矩阵为:D-己(3x3)1-2式中得:b=D=DEU(e)=SU(e)(3x1)(3x3)

23、(3x1)(3x3)(3x8)(5*1)(3x8)41S=DB(获8)(荻3)(48)(13)(14)(15)(16)(17)(18)北京理工大学车辆工程王路北京理工大学车辆工程王路北京理工大学车辆工程王路北京理工大学车辆工程王路(5)单元势能的表达至此，已经通过(9)(10)(17)式将单元的三个基本变量6&U用节点位移矩阵U來进行表达，将这三个式代入到单元的势能表达式中有：T7e=丄U(e)TK(e)U(e)-F(e)TU(e)(19)其中K是四节点矩形单元的刚度矩阵:北京理工大学车辆工程王路北京理工大学车辆工程王路88K(e)=fA()BtDBdAt=(8x8)JA(8x3)(3x3)(

24、3x8)k21k22syin313233kk】2kdk(20)其中t为平面单元的厚度，按照（11）式，可以对（20）式中的各个部分进行分块，各个子块矩阵为：虬二h畸為為t闊，r，s*2,3,4（2x2）（2x3）2754(21)己经求得了儿何矩阵B和弹性矩阵D,故可以对每一个刚度分量进行（2）积分求解，其中的积分过程省略，最终可以得到刚度矩阵：(3x3)EtK(e)Qi)ab(l-/z2)KlS*222333(22)其中:北京理工大学车辆工程王路北京理工大学车辆工程王路比2=亠一(1一)bO%3=M(1yz)a2o对于程序中的问题，给定a=b=0.5,则：K1=%3=丘5=h产&2=K严=%=

25、-台=K厶O虬=匕=监=-+=hkl5=%7=%=広8=-扌+=h16=&5=丘1=%=-=%OO以上定义的匕,匕,垢七七七,即为程序中定义的刚度矩阵。回到程序中，走完刚度矩阵的定义，回到有限元分析子程序中:K=sparse(2*(nelx+1)*(nely+1),2*(nelx+1)*(nely+1);此处开始涉及到单元内所有节点在x和y两个方向上的计算分析，综合图2和图3,观察所有节点的节点定义分析：图4.节点定义分析在所绘制的图中：红色区域是单兀编号，从.x(nely,l).x(nely,nelx);蓝色的是节点的定义和顺序，从节点1,节点2节点(nely+lXnelx+1)；绿色的是每

26、一个节点上X和y方向上的坐标，从2x1-1,2x1,到2x(nely+l)(nelx+1)-12x(nely+l)(nelx+1)。所以定义刚度的时候，考虑到每一个坐标上的应力会对其他任何一个位置上的坐标产生相应应变，所以这一句定义刚度矩阵的行列数都为2(nely+l)(nelx+1)的稀疏矩阵。F=sparse(2*(nely+1)*(nelx+l),1);U=zerosQaiely+laielx+l),!.);对于每一个节点上的位移只有一个，所以定义位移矩阵为行为2(nely+l)(nelx4-1),列为1的零矩阵。由于有：F=KU所以分析行列数，可知力矩阵的行为2(nely+lXnelx

27、+l),列为1,同样定义一个稀疏矩阵。forelx=1:nelxforely=l:nelynl=(nely+1)*(elx-l)+ely;112=(nely+l)*elx+ely;edof=2H：nl-l;2+nl;2*112-1;2*112;2*n2+l;2*n2+2;2*nl+l;2*nl+2;北京理工大学车辆工程王路北京理工大学车辆工程王路K(edof,edof)=K(edof,edof)+x(ely,elx)Apenal*KE;endend解释这一段循环程序，定义ely和elx为单元编号的X和y坐标，即图4中的红色坐标。节点节点ni=(nely+l)(elxT)+elyn2=(nely

28、+l)elx+ely2ni2n222爲1x(ely，elx)nelxl)x(l,nelx)/x(j-l,i)、工1)xnelx)x(jri)x(ji+1)x(nely-l,1)亠_上/x(nely-l,nelx)x(j,i+1)zx(nel，,2)K(nely,/x(nely,1)aelx-1)x(nely,nelx)图7敏度过滤图解分析fori=1:nelxforj=l:nely总循环，以每一个单元的坐标进行循环，其中j是列坐标，从1到nely,i是横坐标，从1到nelx。所以每一次分析都是以每一个单元作为分析，其顺序和图4中蓝色箭头方向一致。sum=0.0;这个变量是用來计算算子fac的变

29、量。fork=max(i-floor(rmin)J):niin(i+floor(rmin),nelx)for1=max(jfloor(miin),l):mii】(j+floor(rmii】)qely)主循环己经确定了某一个单元坐标这段程序是遍历附近的单元，如果满足以单元双川)中心为圆心以rg作为半径的圆，如图7中的黄色区域所示，在圆内(包括边线)上如果包含其他单元的中心坐标的话，就要添加进入加权叠加灵敏度计算。对于语句：fork=max(i-floor(rmin),1):min(i+floor(rmin),nelx)floor是matlab中的向下取整函数。程序语句代表如果在X坐标上，初始状态

30、：max(i-floor(rmin),l),代表圆的北京理工大学车辆工程王路北京理工大学车辆工程王路范围超过了边界，如图7中单元x(l,l)的横坐标(左侧圆黄色部分)所示，则在边界上k的过滤边界(绿色部分)选择较大值1,此时它的最大值为min(i+floor(nnin),nelx)=i+floor(rmin)=l+floor(rmin),故灵敏度过滤的选择单元在X方向上的范围为为1:l+floor(rmin),可以通过图4中x(l,l)附近绿色区域理解。最终状态：min(i+floor(rmin)jielx),同理，如果圆所包含的边界超过了设计范围，如单元x(l,nelx),右侧圆(黄色部分)

31、超过了设计范围，则过滤边界(绿色区域)选择较小值nelx,此时灵敏度过滤在X方向上的范围为nelx-floor(niiin):nelxo可以通过图4中单元x(l,nelx)附近的绿色过滤区域理解。同理，如果在F坐标上，如果是类似于x(l,l)和x(nely,l)单元一样，圆部分(黄色区域)超过了y方向的设计范围，则同理按照以上解释去理解。fac=rmin-sqi-t(i-k)A2+(j-l)A2);sum=sum+max(0曲c);此处解释数学模型：Sigmund提出的敏度过滤技术，是一种局部意义上的约束方法，这种敏度过滤技术通过引入卷积算子(或权重因子)对过滤半径以内的其他单元的敏度进行加权平均，修正单元的灵敏度，从而避免棋盘格现象，通常能够减少棋盘格现象的方法都能减少网格依赖性。所谓棋盘格现象，就是在计算区域内，材料密度为1和0的单元呈现周期性的分布状态，如下图所示，棋盘格的优化结果在工程上的可制造性很差，没有实际意义。所谓的网格依赖性就是对于拓扑优化的计算结果，经常与区域的网格划分密度有关，选择不同的网格划分密度会产生不同的优化计算结果，通常情况下,当网格划分密度变大时，会出现很多的细小分支结构，网格依赖性也使得计算结果的可制造性下降。