引言
Sigmund在2001年在Structural and Multidisciplinary Optimization 发表一篇名为 "A 99 line topology optimization code written in Matlab"论文。该论文附带一个使用Matlab编写的拓扑优化程序。这个只有99行代码的程序基于Matlab环境构建了一个完整的拓扑优化流程:前处理(构建有限元仿真模型), 有限元模型分析计算,拓扑优化迭代和后处理(分析结果显示)。Sigmund在论文中对这个程序的拓扑优化流程做了详细解释。本文仅从程序设计角度解析这段代码。
Sigmund的99行Matlab拓扑优化程序使用模块化方法设计,主要包含以下几个模块:
- 程序主流程
- 有限元模型求解模块
- Filter模块
- 单元刚度阵计算模块:计算平面四边形单元的刚度矩阵
- 优化模块: 使用优化准则法更新设计变量
- 目标函数计算和灵敏度分析模块
下面几个小节首先列出Sigmund99行拓扑优化程序完整代码,接着解释程序中5个全局变量在优化流程中的功用,然后对程序的各个功能模块做详细解析,最后总结程序的局限性并列出扩展程序功能时所需要修改的相应模块。
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
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-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)];
% 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;
else
l2=lmid;
end
end
变量初始化
nelx=60; % x轴方向单元数量
nely=40; % y轴方向单元数量
volfrac=0.5; % 单元的材料体积比
penal=3.; % 惩罚因子,将设计变量从[0,1]转化为指数函数,
% 主要目的是使单元材料密度更加"黑白分明",
% 惩罚因子就是这个指数函数的指数, 其值越大,效果越明显
rmin=1.5; % 过滤半径,目的是预防出现棋盘格现象
设计变量初始化
x(1:nely,1:nelx)=volfrac; % 所有设计变量的赋初始值
主程序
主程序迭代流程
% start ineration
while change>0.01
loop=loop+1; % loop是迭代计数
xold=x; % 将设计变量保存在 变量xold中
% FE analysis
[U]=FE(nelx,nely,x,penal); % 有限元模型分析,将计算得到的各个节点位移值保存在数组U中
% 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; % 右上角节点编号
% 根据以上两个节点的编号 推演出 单元所有节点的自由度编号
% 这就限制了此程序只能使用四边形单元剖分四边设计域所得到的有限元模型
%
% 从数组 U 中根据节点自由度编号 提取该单元的 位移向量
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);
% 计算单元柔度,叠加到目标函数变量中
% 使用SIMP(Solid Isotropic Material with Penalization)材料插值模型
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
优化模型为:
{ m i n c ( x ) = U T K U = ∑ e = 1 n ( x e ) p u e T K e u e s . t . V ρ V = f K U = F x m i n < x < 1 \left\{ \begin{array}{c c c l } min & c(x) & = & U^{T}KU = \sum_{e = 1}^{n} (x_{e})^{p} u_{e}^{T}K_{e} u_{e} \\ s.t. & & \\ & \frac{V_{\rho}}{V} & = & f \\ & KU & = & F \\ & x_{min} & < & x < 1 \\ \end{array} \right. ⎩ ⎨ ⎧mins.t.c(x)VVρKUxmin===<UTKU=∑e=1n(xe)pueTKeuefFx<1
代码 x(ely,elx)^penal*Ue’KEUe 就是计算 $(x_{e})^{p} u_{e}^{T}K_{e} u_{e} $
对 x e x_{e} xe 的梯度是
d c d x = − p ( x e ) p − 1 u e T K e u e \frac{dc}{dx} = - p (x_{e})^{p-1} u_{e}^{T}K_{e} u_{e} dxdc=−p(xe)p−1ueTKeue
此程序使用的SIMP(Solid Isotropic Material with Penalization)材料插值模型,具体代码为:
x(ely,elx)^penal ;
网格划分为
1-------(nely+1)+1------------------------------------nely+1)*nelx
| (1,1) |
2-------(nely+1)+2
|
... -----------(nely+1)*(elx-1)+ely----(nely+1)*elx+ely
| | (elx,ely) |
...
|
nely+1---(nely+1)*2----------
单元节点编号顺序为:
Node_1 Node_2
(nely+1)*(elx-1)+ely----(nely+1)*elx+ely
| (elx,ely) |
(nely+1)*(elx-1)+ely+1----(nely+1)*elx+ely+1
Node_4 Node_3
过滤器模块
%%%%%%%%%% MESH-INDEPENDENCY FILTER %%%%%%%%%%%%%%%%%%%
function [dcn]=check(nelx,nely,rmin,x,dc)
% dcn 清零,dcn 用来保存 更新的 目标函数灵敏度
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)
% sqrt((i-k)^2+(j-l)^2) 是计算此单元与相邻单元的距离
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
Sigmund的论文中对此函数有详细解释,论文中符号与代码中变量名的对应关系如下表所示:
变量 | 论文中符号 |
---|---|
rmin |