拓扑优化学习报告_北理工_王路.pdf

《拓扑优化学习报告_北理工_王路.pdf》由会员分享，可在线阅读，更多相关《拓扑优化学习报告_北理工_王路.pdf（26页珍藏版）》请在三一文库上搜索。

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

2、）和拓 扑优化（Topology Optimization） 。本文我们基于 99 行 matlab 程序初步学习拓扑 优化技术中的理论和优化方法。 拓扑优化技术指的是在给定的设计空间内寻求最佳的材料布局， 同时在满足 平衡方程、物理关系、几何关系和边界约束条件下使得结构达到某种性能最优的 应用技术。 拓扑优化的理论研究最早可以追溯到 Michel 提出的桁架理论，连续体结构 的拓扑优化由于描述和数值计算得困难，发展一直相对缓慢，直到 Bendsoe 和 Kikuchi 在 1988 年提出的均匀化方法之后才得到迅速的发展， 其基本思想是在组 成拓扑结构中引入微结构，通过微结构的几何参数作为设

3、计变量，通过微结构的 增加和删减实现结构的拓扑形状的改变，实现拓扑优化和尺寸优化的统一。 在微结构的基础上，我们介绍变密度法的应用，变密度法是在均匀化方法的 基础上产生的， 把材料引入微结构代之以密度在 01 之间变化的假想材料， 把密 度作为设计变量，从而实现材料的删减，因其模型简单、计算变量相对较少成为 目前广泛采用的方法。根据不同的插值模式，变密度法又有不同的插值模型： SIMP 法（Solid Isotropic Material with Penalization） 、Hashin-Shtrikman 法，以及 RAMP 法（Rational Approximation of Mat

4、erial Properties） 。 （二）拓扑优化问题的描述（二）拓扑优化问题的描述 （1）材料插值模型）材料插值模型 变密度理论的材料插值模型通过引入中间密度单元， 将离散型问题转化成连续型优化问 题，而实际上，中间密度单元是无法存在和制造，因此要尽量避免中间密度单元的产生，减 少中间密度单元的数目，此时就需要对设计变量中出现的中间密度值进行惩罚。 北京理工大学 车辆工程 王路 目前材料插值模型主要分为两类，SIMP 法和 RAMP 法，基于 SIMP 格式的 材料插值模型为： 0,1),()()( 0 imin p imini xEExExE （1） 其中： )E(xi 插值以后的弹性

5、模量 0 E实体部分材料的弹性模量 min E 孔洞部分材料的弹性模量 i x单元相对密度，取值为 1 时表示有材料，为 0 时表示无材料即孔洞 p惩罚因子 在本文中，在本篇文章中，讨论的 SIMP 的表达为： 0 p i,ji,j ExxE)()( （2） 其中： i,j x为第i个子域内第j个单元的相对密度。 （2）数学数学模型模型 根据在体积或质量约束下求最小柔度，即最大刚度来建立优化模型，基于 变密度理论的 SIMP 法的周期性拓扑优化问题的数学模型可以表达为： find RxxxxX T ji ),( ,312111 mi, 21 nj, 21 min m i m i n j ji

6、T ji p ji n j jiji T ji T ukuxukuKUUFUXC 111 0 1 , )()( t s. FKU m i n j jiji xfVV 11 0, 10 max,min xxx ji 其中： X单元相对密度矢量 C结构的柔度 F载荷矢量 U位移矢量 北京理工大学 车辆工程 王路 k结构刚度矩阵 ji u , 单元位移矢量 ji k , 单元刚度矩阵 0 k初始段元刚度矩阵 ji, 单元体积 V优化后的体积 f保留的体积分数 0 V初始体积 min x设计变量的取值下限 max x 设计变量的取值上限 n子域内单元的个数 （三（三）程序程序解读解读 99 行程序主要

7、包括主函数、OC 优化准则、网格过滤、有限元分析、灵敏 度分析 5 个部分，其流程图为： 北京理工大学 车辆工程 王路 设计域的离散化和 程序参数的初始化 主程序 调用有限元分析求 解子程序 灵敏度分析子程序 网格过滤子程序 OC优化准则优化设 计变量 是否满足精度要求 否 结束循环 是 位移U dc dc(new) x(new) 图 1. 99 行程序的流程图 具体程序分析： 总体程序为： % A 99 LINE TOPOLOGY OPTIMIZATION CODE BY OLE SIGMUND, JANUARY 2000 % % CODE MODIFIED FOR INCREASED SP

8、EED, September 2002, BY OLE SIGMUND % function top(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.; % START ITERATION while change 0.01 loop = loop + 1; xold = x; % FE-ANALYSIS U=FE(nelx,nely

9、,x,penal); % OBJECTIVE FUNCTION AND SENSITIVITY ANALYSIS KE = lk; c = 0.; for ely = 1:nely for elx = 1:nelx n1 = (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)

10、*Ue*KE*Ue; end end % FILTERING OF SENSITIVITIES dc = check(nelx,nely,rmin,x,dc); % DESIGN UPDATE BY THE OPTIMALITY CRITERIA METHOD x = OC(nelx,nely,x,volfrac,dc); % PRINT RESULTS change = max(max(abs(x-xold); disp( It.: sprintf(%4i,loop) Obj.: sprintf(%10.4f,c) . Vol.: sprintf(%6.3f,sum(sum(x)/(nelx

11、*nely) . ch.: sprintf(%6.3f,change ) % PLOT DENSITIES 北京理工大学 车辆工程 王路 colormap(gray); imagesc(-x); axis equal; axis tight; axis off;pause(1e-6); end % OPTIMALITY CRITERIA UPDATE % % function xnew=OC(nelx,nely,x,volfrac,dc) l1 = 0; l2 = 100000; move = 0.2; while (l2-l1 1e-4) lmid = 0.5*(l2+l1); xnew =

12、 max(0.001,max(x-move,min(1.,min(x+move,x.*sqrt(-dc./lmid); if sum(sum(xnew) - volfrac*nelx*nely 0; l1 = lmid; else l2 = lmid; end end % MESH-INDEPENDENCY FILTER % % function dcn=check(nelx,nely,rmin,x,dc) dcn=zeros(nely,nelx); for i = 1:nelx for j = 1:nely sum=0.0; for k = max(i-floor(rmin),1):min(

13、i+floor(rmin),nelx) for l = 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); end end 北京理工大学 车辆工程 王路 dcn(j,i) = dcn(j,i)/(x(j,i)*sum); end end % FE- ANALYSIS % % function U=FE(nelx,nely,x,penal) KE = lk;

14、 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); for elx = 1:nelx for ely = 1:nely n1 = (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,ed

15、of) + x(ely,elx)penal*KE; end end % DEFINE LOADS AND SUPPORTS (HALF MBB-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); % SOLVING U(freedofs, ： ）= K(freedofs,freedofs) F(freedofs,： ）;

16、 U(fixeddofs,： ）= 0; % ELEMENT STIFFNESS MATRIX % % function KE=lk E = 1.; 北京理工大学 车辆工程 王路 nu = 0.3; k= 1/2-nu/6 1/8+nu/8 -1/4-nu/12 -1/8+3*nu/8 . -1/4+nu/12 -1/8-nu/8 nu/6 1/8-3*nu/8; KE = E/(1-nu2)* 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)

17、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); % 根据程序的运行的先后顺序，依次详细分析各行程序： function top(nelx,nely,volfrac,penal,rmin); nelx=80;

18、 x 轴方向的单元数目 nely=20; y 轴方向的单元数目 volfrac=0.4; 体积比，这个参数是用来确定保留的材料量 penal=3; 材料插值的惩罚因子，由于变密度理论的材料插值模型引入中 间密度单元，所以将离散型问题转化成连续型优化问题，而在实际过程中，是 无法存在和制造的，所以要尽量避免中间密度单元的产生，减少中间密度单元 的数目，所以需要对中间密度值进行惩罚，penal 是惩罚因子。 rmin=2; 灵敏度过滤的过滤半径 关于体积比、惩罚因子和过滤半径不同引起结果的不同，我们将在之后的讨论 中进行。 以上是初始化参数的程序段 x(1:nely,1:nelx) = volfr

19、ac; x 是设计变量， 为某一个单元的相对密度， 在此 给出整体网格编排的顺序表： 北京理工大学 车辆工程 王路 图 2. 整体网格编排顺序表 故对于整个网格，在垂直方向上，有 nely 个单元，在水平方向上，有 nelx 个单元。 loop = 0; 存放迭代次数的变量 change = 1.; 每次迭代之后，目标函数的改变值，可以用来判断何时收敛 while change 0.01 当两次连续目标函数迭代的差小于 0.01 时，结束迭代 loop = loop + 1; 迭代次数加 1 xold = x; 将前一次的设计变量赋值给 xold U=FE(nelx,nely,x,penal)

20、; 进行有限元分析， 需要输入的参数为 nelx、 nely、 x、penal 进入有限元分析子程序 function U=FE(nelx,nely,x,penal) KE = lk; 单元刚度矩阵分析 进入单元刚度矩阵子程序 function KE=lk E = 1.; k= 1/2-nu/6 1/8+nu/8 -1/4-nu/12 -1/8+3*nu/8 . -1/4+nu/12 -1/8-nu/8 nu/6 1/8-3*nu/8; KE = E/(1-nu2)* 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

21、) 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); 此处单独解释有限元理论，通过平面问题的 4 节点矩形单元描述 （1） 单元的几何和节

22、点描述单元的几何和节点描述 平面4节点矩形单元如图3所示， 单元的节点位移共有8个自由度 （DOF） ， 节点的编号为 1、2、3、4，各自的位置坐标为 ),( ii yx4321,i ，各个节点的位 移（分别沿x和y方向） ，为 4321,),(iu ii 。 图 3.平面 4 节点矩形单元 将所有节点上的位移组成一个列矩阵， 记作 )(e U， 将所有节点上的力组成一 个列矩阵，记作 )(e F，则： T e uuuuU 44332211 )( （3） T yxyxyxyx e FFFFFFFFF 44332211 )( （4） （2） 单元场的位移描述单元场的位移描述 从图 3 可以看出

23、，节点条件共有 8 个，x方向上 4 个),( 4321 uuuu，y 方向上 4 个),( 4321 ，故在x和y方向的位移场可以各有 4 个待定系 数，取以下多项式作为位移场模式： xybybxbbyx xyayaxaayxu 3210 3210 ),( ),( （5） 要考虑到的特点是当x或y不变时，沿y或x方向位移函数是线性变化 的，因此不选 2 x和 2 y项。由于节点条件，在 ii yyxx,处，有： 4321, ),( ),( i yx uyxu iii iii （6） 将四个节点的坐标),(),(),(),( 44332211 yxyxyxyx代入（5）方程中，综合（6） 方程

24、得： ab uuuu a b uuuu a a uuuu a uuuu a 4 4 4 4 4321 3 4321 2 4321 1 4321 0 ab b b b a b b 4 4 4 4 4321 3 4321 2 4321 1 4321 0 经过整理，令（5）式的格式为： 44332211 44332211 ),(),(),(),(),( ),(),(),(),(),( yxNyxNyxNyxNyx uyxNuyxNuyxNuyxNyxu （7） 可以整理得： )(),( )(),( )(),( )(),( b y a x yxN b y a x yxN b y a x yxN b y

25、 a x yxN 11 4 1 11 4 1 11 4 1 11 4 1 4 3 2 1 （8） 将（7）式写成矩阵形式： )( )( )( )( ),( ),( ),( 1882 4 4 3 3 2 2 1 1 4321 4321 12 0000 0000 e UN u u u u NNNN NNNN yx yxu yxu （9） 其中),(yxN为该单元的形状函数矩阵。 （3） 单元应变场的描述单元应变场的描述 由弹性力学中节点的几何方程： x u xx ， y yy ， y u x xy 得到单元应变的表达： 北京理工大学 车辆工程 王路 )( )( )()( )( )( )( )( )

26、( )( ),( 18831882 23 12 23 13 ee xy yy xx UBUNuyx （10） 其中： xy y x 0 0 23)( 单元的形状函数矩阵： 4321 4321 82 0000 0000 NNNN NNNN N )( 定义单元的几何矩阵为： )()()()( )( )( )( ),( 23 4 23 3 23 2 23 1 4321 4321 82 23 83 0000 0000 0 0 BBBB NNNN NNNN xy y x NyxB （11） 式（11）中的子矩阵 i B为： 43210 0 23 , )( i x N y N y N x N B ii i

27、 i i （12） （4） 单元应力场的表达单元应力场的表达 由弹性力学中平面问题的物理方程： xyxy xyyy yxxx E E E )( )( )( 12 1 1 （13） 将（13）式换成应力和应变的关系表达为： xyxy xxyyy yyxxx E E E )( )( )( 12 1 1 2 2 （14） 故有： )()()(1333 2 13 2 1 00 01 01 1 D E xy yy xx xy y x （15） 其中定义 )(33 D为弹性矩阵： 2 1 00 01 01 1 2 33 E D )( （16） 把（10）式代入到（15）式中得： 1883188333133

28、313 )( )()( )( )()()()()( ee USUBDD （17） 其中应力函数矩阵为： )()()(833383 BDS （18） （5） 单元势能的表达单元势能的表达 至此，已经通过（9） （10） （17）式将单元的三个基本变量u,用节点 位移矩阵 )(e U来进行表达，将这三个式代入到单元的势能表达式中有： )()()()()(eTeeeTee UFUKU 2 1 （19） 其中 )(e K是四节点矩形单元的刚度矩阵： 44434241 333231 2221 11 83333888 kkkk kkk symkk k tdABDBK T A e e )()()()( )(

29、 )( （20） 其中t为平面单元的厚度，按照（11）式，可以对（20）式中的各个部分 进行分块，各个子块矩阵为： dxdytBDBk T r A rs e )()( )()( )( 2333 3222 ，4321,sr （21） 已经求得了几何矩阵 )(23 B和弹性矩阵 )(33 D，故可以对每一个刚度分量进行 积分求解，其中的积分过程省略，最终可以得到刚度矩阵： 88 7877 686766 58575655 4847464544 383736353433 28272625242322 1817161514131211 2 88 1 k kk kkk kkkk kkkkk kkkkkk

30、kkkkkkk kkkkkkkk ab Et K e )( )( )( （22） 其中： )( 22 77553311 2 1 3 1 abkkkk )( 22 88664422 2 1 3 1 bakkkk )(31 8 67184523 ab kkkk )(31 8 58362714 ab kkkk )(1 8 38475612 ab kkkk )(1 8 78342516 ab kkkk )( 22 3715 2 1 6 1 abkk )( 22 4826 2 1 6 1 bakk )( 22 4628 4 1 3 1 bakk )( 22 5713 4 1 3 1 abkk 22 68

31、42 1 6 1 bakk)( 22 1753 1 6 1 abkk)( 对于程序中的问题，给定50.ba，则： 18866442277553311 62 1 kkkkkkkkk 238475612 1 8 1 kkkkk)( 346285713 124 1 kkkkk 458362714 8 3 8 1 kkkkk 548263715 124 1 kkkkk 678342516 88 1 kkkkk 768425317 6 kkkkk 867184523 8 3 8 1 kkkkk 以上定义的 87654321 kkkkkkkk,即为程序中定义的刚度矩阵。 回到程序中，走完刚度矩阵的定义，回

32、到有限元分析子程序中： K = sparse(2*(nelx+1)*(nely+1), 2*(nelx+1)*(nely+1); 此处开始涉及到单元内所有节点在x和y两个方向上的计算分析，综合 图 2 和图 3，观察所有节点的节点定义分析： x(1,1) x(nely,1) x(1,nelx) x(nely,nelx) x(i,1) . . . 节点节点1 节点节点2 节点节点i 节点节点i+1 节点节点nely 节点节点nely+1 节点节点 nely+2 节点节点 2（nely+1） 节点编号节点编号 顺序顺序 节点节点 (nelx-1)(nely+1)+1 节点节点 nelx(nely+

33、1) 节点节点 nelx(nely+1)+1 节点节点 (nelx+1)(nely+1) 2x1-1 x(i,nelx) 2x1 2x2-1 2x2 2i-1 2i 2(i+1)-1 2(i+1) 2nely-1 2nely 2(nely+1)-1 2(nely+1) 2(nely+2)-1 2(nely+2) 2x2(nely+1)-1 2x2(nely+1) 2(nely+1)(nelx+1)-1 2(nely+1)(nelx+1) 2nelx(nely+1)+1-1 2nelx(nely+1)+1 x y 单元编号单元编号 每一个节点有每一个节点有 x x和和y y两个方两个方 向上坐标

34、向上坐标 图 4.节点定义分析 在所绘制的图中： 红色区域是单元编号，从),().,().,(),(nelxnelyxnelyxxx12111； 蓝色的是节点的定义和顺序，从节点 1，节点 2节点)(11nelxnely； 绿色的是每一个节点上x和y方向上的坐标，从112，12，到 1112)(nelxnely，)(112nelxnely。 所以定义刚度的时候， 考虑到每一个坐标上的应力会对其他任何一个位置 上 的 坐 标 产 生 相 应 应 变 ， 所 以 这 一 句 定 义 刚 度 矩 阵 的 行 列 数 都 为 )(112nelxnely的稀疏矩阵。 F = sparse(2*(nely

35、+1)*(nelx+1),1); U = zeros(2*(nely+1)*(nelx+1),1); 对 于 每 一 个 节 点 上 的 位 移 只 有 一 个 ， 所 以 定 义 位 移 矩 阵 为 行 为 )(112nelxnely，列为 1 的零矩阵。 由于有：KUF 所以分析行列数， 可知力矩阵的行为)(112nelxnely， 列为 1，同样定义一个稀疏矩阵。 for elx = 1:nelx for ely = 1:nely n1 = (nely+1)*(elx-1)+ely; n2 = (nely+1)* elx +ely; edof = 2*n1-1; 2*n1; 2*n2-1

36、; 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; end end 解释这一段循环程序，定义ely和elx为单元编号的x和y坐标，即图 4 中的红色 坐标。 x(ely,elx) x y 节点节点 n n1 1=(=(nelynely+ +1 1)()(elxelx- -1 1)+)+elyely 节点节点 n n2 2=(=(nelynely+ +1 1) )elxelx+ +elyely 2 2n n1 1- -1 1 2 2n n1 1 2 2n n2 2- -

37、1 1 2 2n n2 2 2 2n n2 2+ +1 1 2 2n n2 2+ +2 2 2 2n n1 1+ +1 1 2 2n n1 1+ +2 2 1 12 2 3 3 4 4 图 5.任意单元附近的节点信息 蓝色坐标即为定义的节点编号，在程序中为了避免繁琐，在一个单元上只 定义两个节点： elyelxnelyn)(11 1 elyelxnelyn)(1 2 然后定义一个单元上四个节点 8 个坐标， 其中的坐标顺序如图 5 中橙色的 1,2,3,4 和箭头所示。 故定义一个单元上的坐标为： edof = 2*n1-1; 2*n1; 2*n2-1; 2*n2; 2*n2+1; 2*n2+

38、2; 2*n1+1; 2*n1+2; 接下来，将单元刚度矩阵组装成总体刚度矩阵： K(edof,edof) = K(edof,edof) + x(ely,elx)penal*KE; 其中是根据材料密度模型中的： 0 kxk p jiji )( , （23） 通过循环程序组装成总体刚度矩阵。程序循环结束之后： F(2*(nelx/2+1)*(nely+1),1) = 1; 此处我们通过这个例子来讲解程序中给定坐标的位置确定和给定位置坐 标的写法： 力矩阵 F(2*(nelx/2+1)*(nely+1),1)，是一个列矩阵，很明显是确定y方向上 的加载力，则可知其节点数为(nelx/2+1)*(n

39、ely+1)，对于节点数，每一列有 (nely+1)个， 可以参考图 4 的表示， 则 nelx/2+1 表示列数， 显然对于一个有 nelx 列的单元，nelx/2+1 表示最中间的那一列节点（在此处，我们一般给定 nelx 为 偶数） 。 fixeddofs = 2*(nely/2+1),2*nelx*(nely+1)+2*(nely/2+1); 此处为定义固定节点的坐标数，很容易知道固定的都是y轴上的坐标， nely/2+1 为第一列节点中间位置的节点数，2*(nely/2+1)为该节点y方向上的坐 标；对于坐标为 2*nelx*(nely+1)+2*(nely/2+1)，我们可以知道，

40、它是节点为 nelx*(nely+1)+ (nely/2+1)在y方向上的约束。 此处要注意的是， 在 matlab 中括号的定义是不同的， 中括号用于形成一个 向量或矩阵，小括号用于表达一般的算术预算、还用于表达下标等。 定义载荷是指下标为),)(/(11122nelynelx的力，是用小括号定义的，它是 指下标，而在固定节点中，是用中括号给定的，该程序中定义固定坐标矩阵是 一个21的矩阵，其中包含固定坐标的的两个坐标数。 在此，我们给出固定节点和施加载荷的位置示意图： F 图 6.边界约束情况 alldofs=1:2*(nely+1)*(nelx+1); 定义所有节点的坐标值是从 1 到

41、2*(nely+1)*(nelx+1)。 用根据之前的理解， 中括号定义的是一个)(1121nelxnely的矩阵。 freedofs=setdiff(alldofs,fixeddofs); 其中利用了 matlab 中的 setdiff 函数，它的作用是返回在 alldofs 中有，但 是在 fixeddofs 中没有的量， 相当于集合论中的BAC， 即),(BAs e t d i f fC 。 此处得到除固定节点坐标外的其他坐标，把它们设为自由坐标。 U(freedofs,： ） = K(freedofs,freedofs) F(freedofs,： ） 之前已经组装了总体刚度矩阵 K，

42、并且在稀疏力矩阵中确定了施加载荷的 大小和方向，此处通过公式 K F U e)( 计算自由坐标上的位移值，从而得到位移 场。 U(fixeddofs,： ）= 0; 同时定义固定节点坐标上的位移为 0，固定的坐标就是在侧边中点y方向 上的位移。 至此，已经解释完有限元分析的程序段。 回到主程序中： KE = lk; 定义单元刚度矩阵，不再重复。 c = 0.; 此处的 c 是存放目标函数为柔度的变量， 我们的优化目标是整个结构的刚 度最大，也就是柔度最小。理解为在一定载荷下，优化结构，使结构的变形最 小。 接下来： for ely = 1:nely for elx = 1:nelx n1 =

43、(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; end end 1 n和 2 n已经解释过了， 作为节点的定义， 其中 )(e U为某一个单元上 4 个节 点 8 个方向上的位移， 其中这个循环中节点坐标的位移矩阵已经通过有限元分 析子程序求解出来了，

44、此处计算每一个单元的柔度，此处利用目标函数的表达 式： c = c + x(ely,elx)penal*Ue*KE*Ue; min m i m i n j ji T ji p ji n j jiji T ji T ukuxukuKUUFUXC 111 0 1 , )()( 0 k为单元刚度矩阵。 dc(ely,elx) = -penal*x(ely,elx)(penal-1)*Ue*KE*Ue; 这句是计算每一个单元的灵敏度的，目标函数的灵敏度表达： 在每一个单元上目标函数c对设计变量x的求导： e T e x U F x c （24） 根据方程KUF ，同时两边对设计变量进行求导： ee x U KU x K 0 U x K K x U ee 1 （25） 将式（25）代入式（24）中： e T e p e e T e ukuxpU x K U x c 0 1 )( （26） 对每一个节点进行计算之后，完成循环，程序来到网格过滤子程序： dc=check(nelx,nely,rmin,x,dc); 进入子程序： dcn=zeros(