99行拓扑优化代码学习，转载于另外两篇博客

‘’‘Matlab
%%%% A 99 LINE TOPOLOGY OPTIMIZATION CODE BY OLE SIGMUND, JANUARY 2000 %%%
function top99_notation(nelx,nely,volfrac,penal,rmin)
% INITIALIZE
x(1:nely,1:nelx) = volfrac;
dc(1:nely,1:nelx) = 0.0;
loop = 0;
change = 1.;
% 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.;      % 保存目标函数(柔度)的变量，%  遍历所有单元% 总的单元数量：nely * nelx;% 网格划分数量：（nely+1）*(nelx+1)for ely=1:nely for elx=1:nelx % 根据单元在设计域中位置 计算单元节点编号n1=(nely+1)*(elx-1)+ely;  % 左上角节点编号  n2=(nely+1)*elx    +ely;  % 右上角节点编号  % 根据以上两个节点的编号 推演出 单元所有节点的自由度编号% 这就限制了此程序只能使用四边形单元剖分四边设计域所得到的有限元模型%% 从数组 U 中根据节点自由度编号 提取该单元的 位移向量 % 左上，右上，右下，左下 自由度% 一个点有两个，所以要*2。第一个从1开始，所以*2之后要-1。%所示单元的自由度，左上，右上，右下，左下，注意顺序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
%%%%%%%%%% 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 > 0l1 = lmid;elsel2 = lmid;end
end
%%%%%%%%%% MESH-INDEPENDENCY FILTER %%%%%%%%%%%%%%%%%%%
function [dcn]=check(nelx,nely,rmin,x,dc)
% dcn 清零，dcn 用来保存 更新的 目标函数灵敏度
dcn=zeros(nely,nelx);
% 遍历所有单元, 加权平均的一个过程
for i = 1:nelxfor j = 1:nelysum=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);endenddcn(j,i) = dcn(j,i)/(x(j,i)*sum);end
end%%%%%%%%%% FE-ANALYSIS %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% FE analysis
function [U]=FE(nelx,nely,x,penal) %自定义函数，最后返回[U]
[KE]=lk; %计算单元刚度矩阵
% sparse 把一个全矩阵转化为一个稀疏矩阵，只存储每一个非零元素的三个值：元素值，元素的行号和列号
% 总体刚度矩阵的稀疏矩阵
% *2是因为x，y都有一个数
% F = KU
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); % 要求的结果u
%组装整体刚度矩阵
for elx=1:nelx for ely=1:nely %一列列的排序，y*xn1=(nely+1)*(elx-1)+ely; % 左上n2=(nely+1)*elx    +ely; % 右上% 左上，右上，右下，左下 自由度% 一个点有两个，所以要*2。第一个从1开始，所以*2之后要-1。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 求解线性方程组, 得到节点自由度的位移值
% 注意在此仅对 freedofs做计算，原理见！！！那里
U(freedofs,:)=K(freedofs,freedofs)\F(freedofs,:);
U(fixeddofs,:)=0;  % 受约束的节点自由度的位移值 设为  0 %%%%%%%%%% ELEMENT STIFFNESS MATRIX %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [KE]=lk
E = 1.; % 杨氏模量，直观可理解为这个材料的硬度，越大则越硬
nu = 0.3; % 泊松比，当这个材料受到压缩时，nu>0则会膨胀（常见），nu<0则会收缩
% 预计算了一个矩形的单元刚度矩阵，具体求解可以看曾攀的《有限元分析》
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)];

原始的博文地址：

Sigmund的99行Matlab拓扑优化程序简析_cocoonyang的专栏-CSDN博客_99行拓扑优化

[更新] 99行拓扑优化代码解析_我要上火星（停止更新，转为私人笔记）-CSDN博客_拓扑优化代码

99行拓扑优化代码学习，转载于另外两篇博客相关推荐

[更新] 99行拓扑优化代码解析
background 1. 什么是拓扑优化所谓拓扑优化,即优化材料的分布,使得最终的结果能够满足结构势能最小,即柔顺度(compliance) min_x c = 1/2 * u^T * K(x) ...
TopOpt | 99行拓扑优化程序完全注释
博客搬家到自己搭建的主页啦q(≧▽≦q),大家快来逛逛鸭! The Corresponding Files (Click to Save): Code: top.m Refe ...
TopOpt | 针对99行改进的88行拓扑优化程序完全注释
博客搬家到自己搭建的主页啦q(≧▽≦q),大家快来逛逛鸭! The Corresponding Files (Click to Save): Code top88.m R ...
Matlab拓扑优化99行步骤,[TopOpt] 针对99行改进的88行拓扑优化程序完全注释
The program can be promoted by line: top88(120,40,0.5,3.0,3.5,1) 代码注释 88行程序为了提高效率,逻辑.流程没有99行程序那么清楚,建 ...
python博客项目评论_Python 爬虫入门——小项目实战（自动私信博客园某篇博客下的评论人，随机发送一条笑话，完整代码在博文最后）...
之前写的都是针对爬虫过程中遇到问题的解决方案,没怎么涉及到实际案例.这次,就以博客园为主题,写一个自动私信博客下的评论人员(在本篇留下的评论的同学也会被自动私信,如果不想被私信,同时又有问题,请私信我 ...
python代码封装供第三方使用_python发博客
广告关闭腾讯云11.11云上盛惠 ,精选热门产品助力上云,云服务器首年88元起,买的越多返的越多,最高返5000元! python生成csdn博客分享图一.前言我们分享博客的方式有很多种,最常见的无 ...
第一篇博客，开启我的编程学习生涯
目录: 一.自我介绍二.编程的目标三.打算怎么学编程四.准备花费多长时间在编程上五.目标的公司六.结束正文: 一.自我介绍大家好,我是大三的学生,现在开始全力准备秋招,为了能在秋招时找到 ...
《软技能——代码之外的生存指南》之博客篇
昨晚拜读软技能-代码之外的生存指南讲真,收益不少,感同身受的太多.而让我决定重新开始写博客,也是在拜读了第21章创建大获成功的博客之后博客的重要性早已烂熟于心,除了总结.记录外,更多的是分 ...
android 智能指针的学习先看邓凡平的书扫盲再看前面两片博客提升
android 智能指针的学习先看邓凡平的书扫盲再看前面两片博客提升转载于:https://www.cnblogs.com/jeanschen/p/3507512.html

99行拓扑优化代码学习，转载于另外两篇博客

99行拓扑优化代码学习，转载于另外两篇博客相关推荐

最新文章

热门文章