拓扑优化丨99行拓扑优化详细解释（OC优化准则）

设计域离散化初始化，有限元分析，灵敏度分析，网格过滤，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)

99 k(8) k(3) k(2) k(5) k(4) k(7) k(6) k(1)];