matlab拓扑优化流程图,Sigmund的99行Matlab拓扑优化程序简析

引言

Sigmund在2001年在Structural and Multidisciplinary Optimization 上发表一篇名为 “A 99 line topology optimization code written in Matlab”论文。该论文后附带了一个Matlab拓扑优化程序。这个只有99行代码程序基于Matlab环境构建了一个完整的拓扑优化流程，其中包含：前处理(构建有限元仿真模型), 有限元模型分析计算,拓扑优化迭代和后处理(分析结果显示)。Sigmund在论文从拓扑优化角度对这个程序做了详细解释。本文仅从程序设计角度解析这段代码。

Sigmund的99行Matlab拓扑优化程序使用模块化方法设计，主要包含以下几个模块：

- 程序主流程

- 有限元模型求解模块

- Filter模块

- 单元刚度阵计算模块：计算平面四边形单元的刚度矩阵

- 优化模块： 使用优化准则法更新设计变量

- 目标函数计算和灵敏度分析模块

99行拓扑优化代码

Sigmund的99行Matlab拓扑优化程序如下所示

% a 99 line topology optimization code by Ole Sigmund,October 1999

clear

nelx=60;

nely=40;

volfrac=0.5;

penal=3.;

rmin=1.5;

% initialize

x(1:nely,1:nelx)=volfrac;

loop=0;

change=1;

% start ineration

while 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: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)*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 result

change=max(max(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 densities

colormap(gray);imagesc(-x);axis equal;axis tight; axis off;pause(1e-6);

end

% 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=sparse(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,edof)+x(ely,elx)^penal*KE;

end

end

% define loads and supports

ip=(nelx+1)*(nely+1);

F(2*ip,1)=-1;

fixeddofs=[1:2*(nely+1)];

alldofs=[1:2*(nely+1)*(nelx+1)];

freedofs=setdiff(alldofs,fixeddofs);

% solving

U(freedofs,:)=K(freedofs,freedofs)\F(freedofs,:);

U(fixeddofs,:)=0;

% 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(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

% Element stiffness matrix

function [KE]=lk