《结构拓扑优化与材料设计》第四次课99行程序解析.pptx

第3次课回顾》结构拓扑优化设计参数：设计变量、状态变量》结构拓扑优化问题列式：最小柔顺性》结构拓扑优化迭代算法：不动点迭代格式》拓扑优化的数值问题2010-10-92计算结构力学第4次课主要内容99行程序解析2010-10-93计算结构力学Matlab 执行top (nelx, nely, volfrac, penal, rmin)》nelx, nely水平方向、竖直方向单元数》Volfrac体积份数》Penal惩罚指数》Rmin过滤半径》 top(60,20,0.5,3.0,1.5)2010-10-94计算结构力学Matlab 执行》 默认边界：半“MBB-beam”》集中荷载作用在左上角点》对称性边界条件作用在左侧边》右下角点作用有竖直支撑算结构力学 52010-10-9计初始化%%%% A 99 LINE TOPOLOGY OPTIMIZATION CODE BY OLESIGMUND, JANUARY 2000 %%%%%%% CODE MODIFIED FOR INCREASED SPEED, September 2002, BY OLE SIGMUND %%%function top(nelx,nely,volfrac,penal,rmin);% INITIALIZEx(1:nely,1:nelx) = volfrac; loop = 0;change = 1.;变量赋初值，=体积份数 循环计数 自变量的变化量% START ITERATIONwhile change 0.01 loop = loop + 1; xold = x;% FE-ANALYSIS[U]=FE(nelx,nely,x,penal);循环计数 变量存储有限元分析位移% OBJECTIVE FUNCTION AND SENSITIVITY ANALYSIS[KE] = lk; c = 0.;单元刚度矩阵，相同 柔顺性指标定义for ely = 1:nelyfor elx = 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;end end灵敏度分析(2*n1-1;2*n1)(2*n2-1;2*n2)(2*n1+1;2*n1+2)(2*n2+1;2*n2+2)连续型优化列式》将 ρ 看成构成结构的密度，将原“黑白”离散设 计转化为具有灰度的连续变量优化。2010-10-99计算结构力学% 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 RESULTSchange = 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 )])% PLOT DENSITIEScolormap(gray); imagesc(-x); axis equal; axis tight; axis off;pause(1e-6);end top主程序结束%%%%%%%%%% 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 = 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; elsel2 = lmid; endend二分法更新拉氏乘子，