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

北京理工大学车辆工程王路基于99行程序的拓扑优化学习报告（一）背景和前言随着汽车工业的飞速发展以及日益突出的能源问题，汽车工业面临的挑战以及竞争环境也越来越激烈，对汽车产品提出了降低其制造成本及燃油经济性的新要求。在提高汽车安全性、减少汽车排放和解决能源消耗的背景下，提出了汽车轻量化技术。实现汽车轻量化的途径包括三个方面：结构优化技术、新型材料和先进性制造工艺。其中，我们所讨论的是结构优化技术，其中结构优化设计分为三个层次：尺寸优化（SizeOptimization）、形状优化（ShapeOptimization）和拓扑优化（TopologyOptimization）。本文我们基于99行matlab程序初步学习拓扑优化技术中的理论和优化方法。拓扑优化技术指的是在给定的设计空间内寻求最佳的材料布局，同时在满足平衡方程、物理关系、几何关系和边界约束条件下使得结构达到某种性能最优的应用技术。拓扑优化的理论研究最早可以追溯到Michel提出的桁架理论，连续体结构的拓扑优化由于描述和数值计算得困难，发展一直相对缓慢，直到Bendsoe和Kikuchi在1988年提出的均匀化方法之后才得到迅速的发展，其基本思想是在组成拓扑结构中引入微结构，通过微结构的几何参数作为设计变量，通过微结构的增加和删减实现结构的拓扑形状的改变，实现拓扑优化和尺寸优化的统一。在微结构的基础上，我们介绍变密度法的应用，变密度法是在均匀化方法的基础上产生的，把材料引入微结构代之以密度在0~1之间变化的假想材料，把密度作为设计变量，从而实现材料的删减，因其模型简单、计算变量相对较少成为目前广泛采用的方法。根据不同的插值模式，变密度法又有不同的插值模型：SIMP法（SolidIsotropicMaterialwithPenalization）、Hashin-Shtrikman法，以及RAMP法（RationalApproximationofMaterialProperties）。（二）拓扑优化问题的描述（1）材料插值模型变密度理论的材料插值模型通过引入中间密度单元，将离散型问题转化成连续型优化问题，而实际上，中间密度单元是无法存在和制造，因此要尽量避免中间密度单元的产生，减少中间密度单元的数目，此时就需要对设计变量中出现的中间密度值进行惩罚。北京理工大学车辆工程王路目前材料插值模型主要分为两类，SIMP法和RAMP法，基于SIMP格式的材料插值模型为：0,1),()()(0iminpiminixEExExE（1）其中：)E(xi——插值以后的弹性模量0E——实体部分材料的弹性模量minE——孔洞部分材料的弹性模量ix——单元相对密度，取值为1时表示有材料，为0时表示无材料即孔洞p——惩罚因子在本文中，在本篇文章中，讨论的SIMP的表达为：0pi,ji,jExxE)()(（2）其中：i,jx为第i个子域内第j个单元的相对密度。（2）数学模型根据在体积或质量约束下求最小柔度，即最大刚度来建立优化模型，基于变密度理论的SIMP法的周期性拓扑优化问题的数学模型可以表达为：findRxxxxXTji),,,,(,,,,312111mi,,,21nj,,,21minmiminjjiTjipjinjjijiTjiTukuxukuKUUFUXC11101,,,,,,)()(ts.FKUminjjijixfVV110,,10max,minxxxji其中：X——单元相对密度矢量C——结构的柔度F——载荷矢量U——位移矢量北京理工大学车辆工程王路k——结构刚度矩阵jiu,——单元位移矢量jik,——单元刚度矩阵0k——初始段元刚度矩阵ji,——单元体积V——优化后的体积f——保留的体积分数0V——初始体积minx——设计变量的取值下限maxx——设计变量的取值上限n——子域内单元的个数（三）程序解读99行程序主要包括主函数、OC优化准则、网格过滤、有限元分析、灵敏度分析5个部分，其流程图为：北京理工大学车辆工程王路设计域的离散化和程序参数的初始化主程序调用有限元分析求解子程序灵敏度分析子程序网格过滤子程序OC优化准则优化设计变量是否满足精度要求否结束循环是位移Udcdc(new)x(new)图1.99行程序的流程图具体程序分析：总体程序为：%%%%A99LINETOPOLOGYOPTIMIZATIONCODEBYOLESIGMUND,JANUARY2000%%%%%%%CODEMODIFIEDFORINCREASEDSPEED,September2002,BYOLESIGMUND%%%functiontop(nelx,nely,volfrac,penal,rmin);nelx=80;nely=20;volfrac=0.4;penal=3;rmin=2;%INITIALIZE北京理工大学车辆工程王路x(1:nely,1:nelx)=volfrac;loop=0;change=1.;%STARTITERATIONwhilechange0.01loop=loop+1;xold=x;%FE-ANALYSIS[U]=FE(nelx,nely,x,penal);%OBJECTIVEFUNCTIONANDSENSITIVITYANALYSIS[KE]=lk;c=0.;forely=1:nelyforelx=1:nelxn1=(nely+1)*(elx-1)+ely;n2=(nely+1)*elx+ely;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;endend%FILTERINGOFSENSITIVITIES[dc]=check(nelx,nely,rmin,x,dc);%DESIGNUPDATEBYTHEOPTIMALITYCRITERIAMETHOD[x]=OC(nelx,nely,x,volfrac,dc);%PRINTRESULTSchange=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;axistight;axisoff;pause(1e-6);end%%%%%%%%%%OPTIMALITYCRITERIAUPDATE%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%function[xnew]=OC(nelx,nely,x,volfrac,dc)l1=0;l2=100000;move=0.2;while(l2-l11e-4)lmid=0.5*(l2+l1);xnew=max(0.001,max(x-move,min(1.,min(x+move,x.*sqrt(-dc./lmid)))));ifsum(sum(xnew))-volfrac*nelx*nely0;l1=lmid;elsel2=lmid;endend%%%%%%%%%%MESH-INDEPENDENCYFILTER%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%function[dcn]=check(nelx,nely,rmin,x,dc)dcn=zeros(nely,nelx);fori=1:nelxforj=1:nelysum=0.0;fork=max(i-floor(rmin),1):min(i+floor(rmin),nelx)forl=max(j-floor(rmin),1):min(j+floor(rmin),nely)fac=rmin-sqrt((i-k)^2+(j-l)^2);sum=sum+max(0,fac);dcn(j,i)=dcn(j,i)+max(0,fac)*x(l,k)*dc(l,k);endend北京理工大学车辆工程王路dcn(j,i)=dcn(j,i)/(x(j,i)*sum);endend%%%%%%%%%%FE-ANALYSIS%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%function[U]=FE(nelx,nely,x,penal)[KE]=lk;K=sparse(2*(nelx+1)*(nely+1),2*(nelx+1)*(nely+1));F=sparse(2*(nely+1)*(nelx+1),1);U=zeros(2*(nely+1)*(nelx+1),1);forelx=1:nelxforely=1:nelyn1=(nely+1)*(elx-1)+ely;n2=(nely+1)*elx+ely;edof=[2*n1-1;2*n1;2*n2-1;2*n2;2*n2+1;2*n2+2;2*n1+1;2*n1+2];K(edof,edof)=K(edof,edof)+x(ely,elx)^penal*KE;endend%DEFINELOADSANDSUPPORTS(HALFMBB-BEAM)F(2*(nelx/2+1)*(nely+1),1)=1;fixeddofs=[2*(nely/2+1),2*nelx*(nely+1)+2*(nely/2+1)];alldofs=[1:2*(nely+1)*(nelx+1)];freedofs=setdiff(alldofs,fixeddofs);%SOLVINGU(freedofs,：）=K(freedofs,freedofs)\F(freedofs,：）;U(fixeddofs,：）=0;%%%%%%%%%%ELEMENTSTIFFNESSMATRIX%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%function[KE]=lkE=1.;北京理工大学车辆工程王路nu=0.3;k=[1/2-nu/61/8+nu/8-1/4-nu/12-1/8+3*nu/8...-1/4+nu/12-1/8-nu/8nu/61/8-3*nu/8];KE=E/(1-nu^2)*[k(1)k(2)k(3)k(4)k(5)k(6)k(7)k(8)k(2)k(1)k(8)k(7)k(6)k(5)k(4)k(3)k(3)k(8)k(1)k(6)k(7)k(4)k(5)k(2)k(4)k(7)k(6)k(1)k(8)k(3)k(2)k(5)k(5)k(6)k(7)k(8)k(1)k(2)k(3)k(4)k(6)k(5)k(4)k(3)k(2)k(1)k(8)k(7)k(7)k(4)k(5)k(2)k(3)k(8)k(1)k(6)k(8)k(3)k(2)k(5)k(4)k(7)k(6)k(1)];%根据程序的运行的先后顺序，依次详细分析各行程序：functiontop(nelx,nely,volfrac,penal,rmin);nelx=80;x轴方向的单元数目nely=20;y轴方向的单元数目volfrac=0.4;体积比，这个参数是用来确定保留的材料量penal=3;材料插值的惩罚因子，由于变密度理论的材料插值模型引入中间密度单元，所以将离散型问题转化成连续型优化问题，而在实际过程中，是无法存在和制造的，所以要尽量避免中间密度单元的产生，减少中间密度单元的数目，所以需要对中间密度值进行惩罚，penal是惩罚因子。rmin=2;灵敏度过滤的过滤半径关于体积比、惩罚因子和过滤半