拓扑优化经典 行程序注释

序%%%%%%%%%%%%%%%%%%%%
143. function [KE]=lk
144.
145. E = 1.;
146. nu = 0.3;
147. k=[ 1/2-nu/6 1/8+nu/8 -1/4-nu/12 -1/8+3*nu/8 ...
148. -1/4+nu/12 -1/8-nu/8 nu/6
157.
k(7) k(4) k(5) k(2) k(3) k(8) k(1) k(6)
158.
k(8) k(3) k(2) k(5) k(4) k(7) k(6) k(1)];
44.
45. % PRINT RESULTS 屏幕上显示迭代信息
46. change = max(max(abs(x-xold))); %计算目标函数的改变量
47. disp([' It.: ' sprintf('%4i',loop) ' Obj.: ' sprintf('%10.4f',c) ...
6);
53. end
54.
55.
56.
57. %%%%%%%%%% OPTIMALITY CRITERIA UPDATE %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
58. function [xnew]=OC(nelx,nely,x,volfrac,dc)
最优化准则
59.
60. l1 = 0; l2 = 100000; %用于体积约束的拉格朗日乘子
%单位刚度矩阵
因为元素刚度矩阵对任何元素都是一样的，所以该子程序仅被调用一次
24. c = 0.;
%用来存放目标函数的变量，这里刚度最大，柔度最小
25.
26. for ely = 1:nely
表示ely从1到nely，+1递增
27. for elx = 1:nelx
28.
n1 = (nely+1)*(elx-1)+ely; %左上角的单元节点
78. function [dcn]=check(nelx,nely,rmin,x,dc)
79.
80. dcn=zeros(nely,nelx);
81.
82. for i = 1:nelx
83. for j = 1:nely
84. sum=0.0;
85.
86. for k = max(i-floor(rmin),1):min(i+floor(rmin),nelx)
29.
n2 = (nely+1)* elx +ely; %右上角的单元节点
30.
31.
%所示单元的自由度，左上，右上，右下，左下
32.
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);
33.
Ue：元素位移矢量； U：整体位移矢量
125. end
126. end
127.
128. % DEFINE LOADS AND SUPPORTS (HALF MBB-BEAM)
129. F(2,1) = -1; % 应用了一个在左上角的垂直单元力。
130. %按着图上来的，最左边和右下角已经固定
131. fixeddofs = union([1:2:2*(nely+1)],[2*(nelx+1)*(nely+1)]); %固定结点
92.
end
93. end
94.
95. dcn(j,i) = dcn(j,i)/(x(j,i)*sum);
96. end
97. end
98.
99.
100. %%%%%%%%%% FE-ANALYSIS 有限元求解子程
序%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
101. function [U]=FE(nelx,nely,x,penal) %自定义函数，最后返回[U]
107. % *2 是因为 x，y 都有一个数
108. K = sparse(2*(nelx+1)*(nely+1), 2*(nelx+1)*(nely+1));
109. %力矩阵的稀疏矩阵
110. F = sparse(2*(nely+1)*(nelx+1),1);
111. U = zeros(2*(nely+1)*(nelx+1),1); %零矩阵
61. move = 0.2;
62.
63. while (l2-l1 > 1e-4)
64. lmid = 0.5*(l2+l1);
65.
66. %即论文公式的综合
67. xnew = max(0.001,max(x-move,min(1.,min(x+move,x.*sqrt(-dc./lmid)))));
12. change = 1.;
%每次迭代，目标函数（柔度）的改变值，用来判断何时
收敛
13.
14. % START ITERATION
15. while change > 0.01 %当两次连续目标函数迭代的差<=0.01 时，迭代结束
16. loop = loop + 1;
17. xold = x; %把前一次的设计变量付给 xold
40. [dc] = check(nelx,nely,rmin,x,dc); %灵敏度过滤，为了边界光顺一点
41. 42. % DESIGN UPDATE BY THE OPTIMALITY CRITERIA METHOD 最优准则法优化更新设计变量
43. [x] = OC(nelx,nely,x,volfrac,dc);
主程序部分包括： 数据初始化；有限元分析；敏度分析，OC算法，结果显示
1. function top(nelx,nely,volfrac,penal,rmin);
2.
3. nelx=80;
% x 轴方向上单元个数
4. nely=20;
% y 轴方向上单元个数
{ }是用于元胞数组，即cell，其 中的元素可以是不同格式的，如 字符和数值，大小也可以不同； [ ] 是用于描述矩阵，初始化或 赋值时使用； ( ) 是用于提取元素，或函数调 用，定义时使用。
102.
103. [KE] = lk; %单元刚度矩阵
104.
105. % sparse 把一个全矩阵转化为一个稀疏矩阵，只存储每一个非零元素的三个值：元素值，元
素的行号和列号 106. %总体刚度矩阵的稀疏矩阵
在矩阵中，若数值为0的元素数目远远多于非0元素的数目， 并且非0元素分布没有规律时，则称该矩阵为稀疏矩阵
68.
69. if sum(sum(xnew)) - volfrac*nelx*nely > 0;%二乘法减半
70. l1 = lmid;
71. else
72. l2 = lmid;
73. end
74. end
75.
76.
77. %%%%%%%%%% MESH-INDEPENDENCY FILTER 敏度过滤技术子程序%%%%%%%%%%%%%%%%%%%
5. volfrac=0.4;
%体积比
6. penal=3;
%材料插值的惩罚因子
7. rmin=2;
%敏度过滤半径
8.
9. % INITIALIZE
10. x(1:nely,1:nelx) = volfrac; % x 是设计变量(单元伪密度) 材料设计区域的均匀分配
11. loop = 0;
%存放迭代次数的变量
153.
k(3) k(8) k(1) k(6) k(7) k(4) k(5) k(2)
154.
k(4) k(7) k(6) k(1) k(8) k(3) k(2) k(5)
155.
k(5) k(6) k(7) k(8) k(1) k(2) k(3) k(4)
156.
k(6) k(5) k(4) k(3) k(2) k(1) k(8) k(7)
87.
for l = max(j-floor(rmin),1):min(j+floor(rmin),nely)
88.
fac = rmin-sqrt((i-k)^2+(j-l)^2);
89.
sum = sum + max(0,fac);
90.
91.
dcn(j,i) = dcn(j,i) + max(0,fac)*x(l,k)*dc(l,k);
120. % 一个点有两个，所以要*2。第一个从 1 开始，所以*2 之后要-1。
121. edof = [2*n1-1;2*n1; 2*n2-1;2*n2; 2*n2+1;2*n2+2; 2*n1+1;2*n1+2];
122.
123. %将单元刚度矩阵组装成总的刚度矩阵
124. K(edof,edof) = K(edof,edof) + x(ห้องสมุดไป่ตู้ly,elx)^penal*KE;
34.
c = c + x(ely,elx)^penal*Ue'*KE*Ue; %计算目标函数的值（柔度）
35.
dc(ely,elx) = -penal*x(ely,elx)^(penal-1)*Ue'*KE*Ue; %目标函数的灵敏度
36. end