设计域离散化初始化，有限元分析，灵敏度分析，网格过滤，OC优化准则设计变量
1 %%%% A 99 LINE TOPOLOGY OPTIMIZATION CODE BY OLESIGMUND, OCTOBER 1999 %%% 99行程序代码
2 function top(nelx,nely,volfrac,penal,rmin);水平方向上的离散单元，竖直方向的离散单元，材料体积与设计域体积之比，惩罚因子，灵敏度过滤半径
3 % INITIALIZE初始化
4 x(1:nely,1:nelx) = volfrac; x为设计变量，给设计域内单元一个初始相对密度
5 loop = 0;迭代次数
6 change = 1.; change是储存迭代之后目标函数的改变值，判断是否收敛。
7 % START ITERATION开始迭代
8 while change > 0.01 改变量小于等于0.01时结束迭代
9 loop = loop + 1;迭代次数加1
10 xold = x; 前一次的设计变量赋值给xold
11 % FE-ANALYSIS有限元分析
12 [U]=FE(nelx,nely,x,penal); 对每次迭代都进行有限元分析，计算节点位移，储存在全局位移U中
13 % OBJECTIVE FUNCTION AND SENSITIVITY ANALYSIS目标函数和灵敏度分析
14 [KE] = lk;计算单元刚度矩阵
15 c = 0.; 储存目标函数变量
16 for ely = 1:nely
17 for elx = 1:nelx
18 n1 = (nely+1)*(elx-1)+ely; 左上角节点编号
19 n2 = (nely+1)* elx +ely; 右上角节点编号
20 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); 局部位移数组储存4个节点，8个自由度
21 c = c + x(ely,elx)^penal*Ue’*KE*Ue;目标函数
22 dc(ely,elx) = -penal*x(ely,elx)^(penal-1)*Ue’*KE*Ue;总体结构的灵敏度
23 end
24end
目标函数：
灵敏度推导过程：
25 % FILTERING OF SENSITIVITIES灵敏度过滤
26 [dc] = check(nelx,nely,rmin,x,dc);
27 % DESIGN UPDATE BY THE OPTIMALITY CRITERIA METHOD oc优化设计准则
28 [x] = OC(nelx,nely,x,volfrac,dc);
29 % PRINT RESULTS 打印结果
30 change = max(max(abs(x-xold))); 更新目标函数改变值
31 disp([’ It.: ’ sprintf(’%4i’,loop) ’ Obj.: ’
sprintf(’%10.4f’,c) ...
32 ’ Vol.: ’ sprintf(’%6.3f’,sum(sum(x))/(nelx*nely)) ...
33 ’ ch.: ’ sprintf(’%6.3f’,change )])
34 % PLOT DENSITIES
35 colormap(gray); imagesc(-x); axis equal; axis
tight; axis off;pause(1e-6);
36 end
37 %%%%%%%%%% OPTIMALITY CRITERIA UPDATE %%%%%%%%%
38 function [xnew]=OC(nelx,nely,x,volfrac,dc) 输入小括号中的，输出为xnew
39 l1 = 0; l2 = 100000; move = 0.2; 定义取值区间，利用二分法得到满足体积约束的拉格朗日算子，最大正向位移
40 while (l2-l1 > 1e-4)
41 lmid = 0.5*(l2+l1); 二分法取中间点
42 xnew = max(0.001,max(x-move,min(1.,min(x+move,x.
*sqrt(-dc./lmid)))));
43 if sum(sum(xnew)) - volfrac*nelx*nely > 0;
44 l1 = lmid;
45 else
46 l2 = lmid;
47 end
48 end
49 %%%%%%%%%% MESH-INDEPENDENCY FILTER %%%%%%%%%%%无关网格敏度过率子程序
50 function [dcn]=check(nelx,nely,rmin,x,dc)
51 dcn=zeros(nely,nelx); dcn清零，保存更新的目标函数灵敏度
52 for i = 1:nelx
53 for j = 1:nely
54 sum=0.0;
55 for k = max(i-round(rmin),1):
min(i+round(rmin),nelx)
56 for l = max(j-round(rmin),1):
min(j+round(rmin), nely)
57 fac = rmin-sqrt((i-k)^2+(j-l)^2);
58 sum = sum+max(0,fac);
59 dcn(j,i) = dcn(j,i) + max(0,fac)*x(l,k)
*dc(l,k);
60 end
61 end
62 dcn(j,i) = dcn(j,i)/(x(j,i)*sum);
63 end
64 end
65 %%%%%%%%%% FE-ANALYSIS %%%%%%%%%%%% 有限元求解子程序
66 function [U]=FE(nelx,nely,x,penal) 输入小括号中，输出为全局节点位移
67 [KE] = lk; 单元刚度矩阵
68 K = sparse(2*(nelx+1)*(nely+1), 2*(nelx+1)*(nely+1)); 总体刚度矩阵的稀疏矩阵
69 F = sparse(2*(nely+1)*(nelx+1),1);力矩阵的稀梳矩阵
U =zeros(2*(nely+1)*(nelx+1),1); U清零用来保存更新的全局节点位移
70 for ely = 1:nely
71 for elx = 1:nelx
72 n1 = (nely+1)*(elx-1)+ely;计算节点编号，同上
73 n2 = (nely+1)* elx +ely;
74 edof = [2*n1-1; 2*n1; 2*n2-1; 2*n2; 2*n2+1;2*n2+2;2*n1+1; 2*n1+2];
75 K(edof,edof) = K(edof,edof) +x(ely,elx)^penal*KE; 将单元刚度矩阵组成总体刚度矩阵
76 end
77 end
78 % DEFINE LOADSAND SUPPORTS(HALF MBB-BEAM) 定义载荷支撑
79 F(2,1) = -1; 左上角垂直向下的力
80 fixeddofs = union([1:2:2*(nely+1)], [2*(nelx+1)*(nely+1)]); 施加约束，左边第一列和右下角固定，每个节点自由度有两个，所以乘以2
81 alldofs = [1:2*(nely+1)*(nelx+1)]; 所有自由度
82 freedofs = setdiff(alldofs,fixeddofs); 从所有自由度中除去固定自由度
83 % SOLVING
84 U(freedofs,:) = K(freedofs,freedofs) \F(freedofs,:); 求解方程，将各节点自由度储存到U之中
85 U(fixeddofs,:)= 0; 固定节点自由度为0
86 %%%%%%%%%% ELEMENT STIFFNESS MATRIX %%%%%%% 求解单元刚度矩阵子程序
87 function [KE]=lk
88 E = 1.;
89 nu = 0.3;
90 k=[ 1/2-nu/6 1/8+nu/8 -1/4-nu/12 -1/8+3*nu/8 ...
91 -1/4+nu/12 -1/8-nu/8 nu/6 1/8-3*nu/8];
92 KE = E/(1-nu^2)*
[ k(1) k(2) k(3) k(4) k(5) k(6) k(7) k(8)
93 k(2) k(1) k(8) k(7) k(6) k(5) k(4) k(3)
94 k(3) k(8) k(1) k(6) k(7) k(4) k(5) k(2)
95 k(4) k(7) k(6) k(1) k(8) k(3) k(2) k(5)
96 k(5) k(6) k(7) k(8) k(1) k(2) k(3) k(4)
97 k(6) k(5) k(4) k(3) k(2) k(1) k(8) k(7)
98 k(7) k(4) k(5) k(2) k(3) k(8) k(1) k(6)文章来源：https://www.toymoban.com/news/detail-419737.html
99 k(8) k(3) k(2) k(5) k(4) k(7) k(6) k(1)];文章来源地址https://www.toymoban.com/news/detail-419737.html
到了这里，关于拓扑优化丨99行拓扑优化详细解释（OC优化准则）的文章就介绍完了。如果您还想了解更多内容，请在右上角搜索TOY模板网以前的文章或继续浏览下面的相关文章，希望大家以后多多支持TOY模板网！
