找回密码
 注册
Simdroid-非首页
查看: 3128|回复: 63

拓扑优化经典99行程序解读

[复制链接]
发表于 2011-7-10 18:18:11 | 显示全部楼层 |阅读模式 来自 湖北武汉
本帖最后由 yygao56 于 2011-7-10 22:02 编辑

Sigmund教授所编写的top优化经典99行程序,可以说是我们拓扑优化研究的基础;
每一个新手入门都会要读懂这个程序,才能去扩展,去创新;
99行程序也有好多个版本,用于求解各种问题,如刚度设计、柔顺机构、热耦合问题,但基本思路大同小异;
本文拟对其中的一个版本进行解读,愿能对新手有点小小的帮助。
不详之处,还请论坛内高手多指点
读懂了该程序,只能说是略懂拓扑优化理论了,
我手里就有一些水平集源程序是成千上万行,虽然在99行的基础上成熟了很多,但依然还有很多的发展空间。

源程序如下:
%%%% A 99 LINE TOPOLOGY OPTIMIZATION CODE BY OLE SIGMUND, JANUARY 2000 %%%
%%%% CODE MODIFIED FOR INCREASED SPEED, September 2002, BY OLE SIGMUND %%%
function top(nelx,nely,volfrac,penal,rmin);
nelx=80;
nely=20;
volfrac=0.4;
penal=3;
rmin=2;
% INITIALIZE
x(1:nely,1:nelx) = volfrac;
loop = 0;
change = 1.;
% START ITERATION
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 RESULTS
  change = max(max(abs(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 > 0;
    l1 = lmid;
  else
    l2 = lmid;
  end
end
%%%%%%%%%% 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
%%%%%%%%%% 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 = zeros(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 (HALF MBB-BEAM)
F(2*(nelx/2+1)*(nely+1),1) = 1;
fixeddofs   = [2*(nely/2+1),2*nelx*(nely+1)+2*(nely/2+1)];
alldofs     = [1:2*(nely+1)*(nelx+1)];
freedofs    = setdiff(alldofs,fixeddofs);
% SOLVING
U(freedofs, : )= K(freedofs,freedofs) \ F(freedofs,: );      
U(fixeddofs,: )= 0;   %  这两行复制后应换成英文字符,我这里为了防止生成QQ表

                                       %   情修改了一下格式
%%%%%%%%%% 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)];
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
程序执行方法:(just for matlab new users)
打开matlab,点开new M-file,将上述源程序复制粘贴到M-文件中,修改蓝色部分的格式,保存。按F5即可执行~~




程序个人解读(会针对大家的提问,高手们的解释,不断补充更新):
主程序部分:
包括:  数据初始化;有限元分析;敏度分析,OC算法,结果显示

function top(nelx,nely,volfrac,penal,rmin);
nelx=80;   % x轴方向的单元数
nely=20;   % y轴方向单元数
volfrac=0.4;  %体积比
penal=3;      %材料插值的惩罚因子
rmin=2;       %敏度过滤的半径
% INITIALIZE
x(1:nely,1:nelx) = volfrac;   %x是设计变量
loop = 0;                     %存放迭代次数的变量
change = 1.;                  %每次迭代目标函数的改变值,用以判断何时收敛
% START ITERATION
while change > 0.01           %当两次连续目标函数迭代的差小于等于0.01时,结束迭代   
  loop = loop + 1;            %迭代次数加1
  xold = x;                   %把前一次的设计变量赋值给xold
% 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 RESULTS
  change = max(max(abs(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
OC算法子程序:
function [xnew]=OC(nelx,nely,x,volfrac,dc)  
l1 = 0; l2 = 100000; move = 0.2;  %l1,l2用于体积约束的拉格朗日乘子
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))))); %这里是OC算法的核心所在,具体含义可参考论文中的公式
  if sum(sum(xnew)) - volfrac*nelx*nely > 0;   %采用了二乘法更新拉格朗日乘子  
    l1 = lmid;
  else
    l2 = lmid;
  end
end

敏度过滤技术子程序:
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
这一段就不多解释了,只是为了光顺边界的,现在二重敏度过滤技术用得更多一点了。理论部分可参考罗震博士的毕业论文
看不懂的代码在matlab命令窗口输入:help XXX(即你看不懂的那个关键词,比如这里的floor)

有限元求解子程序
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 = zeros(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];   %这里的Y轴是反向的,但是不影响最后的结果,详情请见二楼TYNGOD这位高手的解释,感谢TYNGOD。
    K(edof,edof) = K(edof,edof) + x(ely,elx)^penal*KE;   %将单元刚度矩阵组装成总的刚度矩阵
  end
end
% DEFINE LOADS AND SUPPORTS (HALF MBB-BEAM)
F(2*(nelx/2+1)*(nely+1),1) = 1;                         %初始的集中力
fixeddofs   = [2*(nely/2+1),2*nelx*(nely+1)+2*(nely/2+1)];  %固定结点
alldofs     = [1:2*(nely+1)*(nelx+1)];                      %所有结点  
freedofs    = setdiff(alldofs,fixeddofs);                   %自由节点
% SOLVING
U(freedofs,:) = K(freedofs,freedofs) \ F(freedofs,:);       %有限元求解:位移场
U(fixeddofs,:)= 0;                                          %固定节点位移为0

单元刚度矩阵的子程序:
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)];
这里不解释,自己可查有限元理论。个人认为这里的单元刚度求解不是很好。

最后是我的运行结果:
优化结果:

论文:A 99 line topology optimization code written in Matlab

如果不看论文,只看我这个程序,是读不懂的。建议新手先大概了解一下拓扑优化的基本知识

本帖子中包含更多资源

您需要 登录 才可以下载或查看,没有账号?注册

×

评分

1

查看全部评分

发表于 2011-7-21 15:27:56 | 显示全部楼层 来自 湖北十堰
Simdroid开发平台
3# 417332551  
有空把你的算法程序发给我看看吧!
邮箱:yygao56@smail.hust.edu.cn
多谢啦
yygao56 发表于 2011-7-11 18:58

见习版主原来是华科的啊,校友好!
回复 1 不支持 0

使用道具 举报

发表于 2011-7-10 20:47:08 | 显示全部楼层 来自 上海浦东新区
请大家注意Sigmund教授发表的代码中,定义坐标轴的Y轴方向为坚直向下,与通常我们熟悉的Y轴竖直向上为正是相反的。从文中和代码中都可以证实,说明如下:
文中指出:The variables n1 and n2 denote upper left and right element node numbers in global node numbers。也就是说n1和n2两个节点分别是四边形单元的左上角点和右上角点。如下图所示,

n1++++++++++++++n2
|                                |
|                                |
|                                |
|                                |
n4++++++++++++++n3  ( n1 -> n2 -> n3 -> n4 )
也就是说,Sigmund指定单元四个节点的顺序为顺时针,分别为:左上,右上,右下,左下 四个节点。(通常我们是从左下节点开始逆时针方向定义节点编号的)所示单元的自由度分别有:
edof = [2*n1-1; 2*n1; 2*n2-1; 2*n2; 2*n2+1; 2*n2+2; 2*n1+1; 2*n1+2]

所以,我们可以证实单元内四个节点的编号(n1, n2, n3, n4)有如下的关系:
n1+1=n4;
n2+1=n3;
n1+(nely+1)=n2; ( nely 为结构网格划分后竖向的单元数目)
所以说明,Sigmund教授定义的节点编号是竖直向下编号的,则Y轴是竖直向下为正,其原点为结构的左上角点。

但对于这种刚度优化问题,这个定义是不会影响柔顺度数值的。
此问题已得到学者证实!!

评分

1

查看全部评分

回复 1 不支持 0

使用道具 举报

发表于 2011-7-10 20:50:03 | 显示全部楼层 来自 江苏镇江
本帖最后由 417332551 于 2011-7-10 20:54 编辑

哥哥,可以用不同的颜色标注一下。。

你真是个大好人。一定要顶你当OS的斑竹。

乐于助人,授人以渔。。才是大侠风范!!!此人日后必成大器。

帖子发表完后,还可以再编辑编辑。。用不同颜色区分。。


另外有一些QQ表情。。


有时间,我给你一部分OC算法和mma算法的源程序。(暂时没空)


楼主,辛苦啦你!
回复 不支持

使用道具 举报

 楼主| 发表于 2011-7-10 21:08:55 | 显示全部楼层 来自 湖北武汉
2# TYNGOD
是的。程序中Y轴方向是反的,但是因为是弹性结构,所以对最后的拓扑没有影响。
回复 不支持

使用道具 举报

 楼主| 发表于 2011-7-10 21:11:41 | 显示全部楼层 来自 湖北武汉
3# 417332551
嗯,因为我发帖都是先word编辑,或者记事本编辑的,然后复制进来。所以对论坛内发帖的很多功能不太会。我试试改颜色,谢谢你的建议
回复 不支持

使用道具 举报

发表于 2011-7-10 21:14:07 | 显示全部楼层 来自 湖北武汉
威武,版块里一大人才,大家有啥不懂的赶紧问,莫让他闲着,嘿嘿
回复 不支持

使用道具 举报

发表于 2011-7-11 15:49:52 | 显示全部楼层 来自 江苏镇江
哈哈哈哈。你干的不错!有心就好。
有心做公益就好!!

[b] 1# yygao56
回复 不支持

使用道具 举报

 楼主| 发表于 2011-7-11 18:58:41 | 显示全部楼层 来自 湖北武汉
3# 417332551
有空把你的算法程序发给我看看吧!
邮箱:yygao56@smail.hust.edu.cn
多谢啦
回复 不支持

使用道具 举报

 楼主| 发表于 2011-7-21 17:19:12 | 显示全部楼层 来自 湖北武汉
9# yrunze
话说,学校邮箱以前都不用的~~最近才启用的。
回复 不支持

使用道具 举报

发表于 2011-7-28 09:23:08 | 显示全部楼层 来自 江苏镇江
本帖最后由 417332551 于 2011-7-28 10:01 编辑

记得查收




你先把这个注释,详细讲解一下。。然后我再给你传资料。。


因为这个是我很久很久以前的师兄做的硕士论文。。。不知道我上传他的东西,会不会出问题。

10# yygao56

本帖子中包含更多资源

您需要 登录 才可以下载或查看,没有账号?注册

×
回复 不支持

使用道具 举报

发表于 2011-7-28 09:34:48 | 显示全部楼层 来自 江苏镇江
本帖最后由 417332551 于 2011-7-28 10:00 编辑

It.:    1 Obj.:   333.0403 Vol.:  0.400 ch.:  0.200
It.:    2 Obj.:   208.0136 Vol.:  0.400 ch.:  0.200
It.:    3 Obj.:   151.8037 Vol.:  0.400 ch.:  0.200
It.:    4 Obj.:   123.7844 Vol.:  0.400 ch.:  0.200

这个是你的99行代码。。我正看里面的计算公式。
It.:    5 Obj.:   111.9788 Vol.:  0.400 ch.:  0.173
It.:    6 Obj.:   103.2981 Vol.:  0.400 ch.:  0.200
It.:    7 Obj.:    97.1750 Vol.:  0.400 ch.:  0.167
It.:    8 Obj.:    90.7349 Vol.:  0.400 ch.:  0.200
It.:    9 Obj.:    85.5328 Vol.:  0.400 ch.:  0.162
It.:   10 Obj.:    80.2626 Vol.:  0.400 ch.:  0.179
It.:   11 Obj.:    75.5167 Vol.:  0.400 ch.:  0.154
It.:   12 Obj.:    70.9125 Vol.:  0.400 ch.:  0.159
It.:   13 Obj.:    66.6570 Vol.:  0.400 ch.:  0.124
It.:   14 Obj.:    63.3128 Vol.:  0.400 ch.:  0.127
It.:   15 Obj.:    60.6792 Vol.:  0.400 ch.:  0.094
It.:   16 Obj.:    59.0614 Vol.:  0.400 ch.:  0.094
It.:   17 Obj.:    57.8935 Vol.:  0.400 ch.:  0.082
It.:   18 Obj.:    57.1503 Vol.:  0.400 ch.:  0.070
It.:   19 Obj.:    56.5448 Vol.:  0.400 ch.:  0.074
It.:   20 Obj.:    56.1134 Vol.:  0.400 ch.:  0.057
It.:   21 Obj.:    55.8116 Vol.:  0.400 ch.:  0.059
It.:   22 Obj.:    55.5296 Vol.:  0.400 ch.:  0.047
It.:   23 Obj.:    55.3214 Vol.:  0.400 ch.:  0.048
It.:   24 Obj.:    55.1301 Vol.:  0.400 ch.:  0.038
It.:   25 Obj.:    55.0099 Vol.:  0.400 ch.:  0.038
It.:   26 Obj.:    54.8593 Vol.:  0.400 ch.:  0.034
It.:   27 Obj.:    54.7776 Vol.:  0.400 ch.:  0.034
It.:   28 Obj.:    54.6600 Vol.:  0.400 ch.:  0.031
It.:   29 Obj.:    54.5681 Vol.:  0.400 ch.:  0.029
It.:   30 Obj.:    54.5063 Vol.:  0.400 ch.:  0.027
It.:   31 Obj.:    54.4336 Vol.:  0.400 ch.:  0.027
It.:   32 Obj.:    54.3624 Vol.:  0.400 ch.:  0.025
It.:   33 Obj.:    54.2916 Vol.:  0.400 ch.:  0.024
It.:   34 Obj.:    54.2521 Vol.:  0.400 ch.:  0.023
It.:   35 Obj.:    54.1981 Vol.:  0.400 ch.:  0.022
It.:   36 Obj.:    54.1318 Vol.:  0.400 ch.:  0.022
It.:   37 Obj.:    54.0924 Vol.:  0.400 ch.:  0.019
It.:   38 Obj.:    54.0483 Vol.:  0.400 ch.:  0.020
It.:   39 Obj.:    54.0072 Vol.:  0.400 ch.:  0.017
It.:   40 Obj.:    53.9663 Vol.:  0.400 ch.:  0.018
It.:   41 Obj.:    53.9645 Vol.:  0.400 ch.:  0.016
It.:   42 Obj.:    53.9046 Vol.:  0.400 ch.:  0.017
It.:   43 Obj.:    53.8879 Vol.:  0.400 ch.:  0.016
It.:   44 Obj.:    53.8638 Vol.:  0.400 ch.:  0.016
It.:   45 Obj.:    53.8307 Vol.:  0.400 ch.:  0.015
It.:   46 Obj.:    53.8200 Vol.:  0.400 ch.:  0.015
It.:   47 Obj.:    53.8020 Vol.:  0.400 ch.:  0.015
It.:   48 Obj.:    53.7849 Vol.:  0.400 ch.:  0.015
It.:   49 Obj.:    53.7713 Vol.:  0.400 ch.:  0.015
It.:   50 Obj.:    53.7176 Vol.:  0.400 ch.:  0.015
It.:   51 Obj.:    53.7128 Vol.:  0.400 ch.:  0.015
It.:   52 Obj.:    53.6724 Vol.:  0.400 ch.:  0.015
It.:   53 Obj.:    53.6572 Vol.:  0.400 ch.:  0.015
It.:   54 Obj.:    53.6178 Vol.:  0.400 ch.:  0.015
It.:   55 Obj.:    53.5669 Vol.:  0.400 ch.:  0.016
It.:   56 Obj.:    53.5105 Vol.:  0.400 ch.:  0.017
It.:   57 Obj.:    53.4867 Vol.:  0.400 ch.:  0.017
It.:   58 Obj.:    53.4436 Vol.:  0.400 ch.:  0.017
It.:   59 Obj.:    53.3850 Vol.:  0.400 ch.:  0.017
It.:   60 Obj.:    53.3274 Vol.:  0.400 ch.:  0.017
It.:   61 Obj.:    53.2989 Vol.:  0.400 ch.:  0.017
It.:   62 Obj.:    53.2237 Vol.:  0.400 ch.:  0.016
It.:   63 Obj.:    53.1788 Vol.:  0.400 ch.:  0.016
It.:   64 Obj.:    53.1438 Vol.:  0.400 ch.:  0.016
It.:   65 Obj.:    53.0900 Vol.:  0.400 ch.:  0.015
It.:   66 Obj.:    53.0506 Vol.:  0.400 ch.:  0.014
It.:   67 Obj.:    53.0082 Vol.:  0.400 ch.:  0.014
It.:   68 Obj.:    52.9708 Vol.:  0.400 ch.:  0.013
It.:   69 Obj.:    52.9630 Vol.:  0.400 ch.:  0.012
It.:   70 Obj.:    52.9436 Vol.:  0.400 ch.:  0.011
It.:   71 Obj.:    52.9106 Vol.:  0.400 ch.:  0.010
It.:   72 Obj.:    52.8900 Vol.:  0.400 ch.:  0.009


慢慢研究研究里面的公式。。。希望你早点上学,回到论坛。。。。祝愿大家将此帖子深加工。研究式的学习和探讨,不能浮于表面代码的复制 粘贴 运行。

本帖子中包含更多资源

您需要 登录 才可以下载或查看,没有账号?注册

×
回复 不支持

使用道具 举报

发表于 2011-7-28 10:47:58 | 显示全部楼层 来自 江苏镇江
F(2*(nelx/2+1)*(nely+1),1) = 1;                         %初始的集中力
fixeddofs   = [2*(nely/2+1),2*nelx*(nely+1)+2*(nely/2+1)];  %固定结点
alldofs     = [1:2*(nely+1)*(nelx+1)];                      %所有结点  
freedofs    = setdiff(alldofs,fixeddofs);                   %自由节点
回复 不支持

使用道具 举报

发表于 2011-7-28 23:09:38 | 显示全部楼层 来自 湖北武汉
13# 417332551
于是乎,你决定要仔细研究算法了?
回复 不支持

使用道具 举报

发表于 2011-7-29 09:19:47 | 显示全部楼层 来自 江苏镇江
有些地方看不懂,有些公式不清楚。要查阅很多资料的。


区别一:这些都是二维结构的的。直接定义坐标轴节点数目。工程上结构很复杂,结构不规则。
二、三维结构的话,好像很难定义节点,设计空间。

有点难度!斑竹这个程序和PDF那个程序做了一些相应的变动。


你看看我上传的几个文档。。。。有一个是自动输入  坐标轴上的   节点数目以及其他参数,matlab做的一个原型系统。。


不是我做的啊!!!
14# kmani
回复 不支持

使用道具 举报

发表于 2011-8-7 23:49:46 | 显示全部楼层 来自 湖南长沙
在此多谢!请问各位前辈有基于ansys的APDL语言的拓扑优化程序吗,ansys自带的拓扑模块无法优化shell63壳单元
回复 不支持

使用道具 举报

发表于 2011-8-12 11:08:50 | 显示全部楼层 来自 湖北武汉
斑竹威武,不解释!
回复 不支持

使用道具 举报

发表于 2011-8-31 10:56:24 | 显示全部楼层 来自 广西柳州
蛮不错的,以前没有任何资料,我花了个把月才看懂一个程序,早知道有你,我相信我就不会过得这么纠结,qq358939206,有缘拓扑者可加我多交流
回复 不支持

使用道具 举报

发表于 2011-9-2 10:29:34 | 显示全部楼层 来自 湖北宜昌
n1++++++++++++++n2
|                                |
|                                |
|                                |
|                                |
n4++++++++++++++n3  ( n1 -> n2 -> n3 -> n4 )
前几天看了代码,发现节点的顺序也是这样的,原来是Y轴方向错误,坛子里高手真多。
回复 不支持

使用道具 举报

发表于 2011-9-2 11:01:12 | 显示全部楼层 来自 湖北宜昌
想请教楼主,99行代码中的敏度推导是不是这样的。以前一直对敏度前面的一个负号耿耿于怀,现在推了一遍,还对倒数第二步有点疑问。

本帖子中包含更多资源

您需要 登录 才可以下载或查看,没有账号?注册

×
回复 不支持

使用道具 举报

您需要登录后才可以回帖 登录 | 注册

本版积分规则

Archiver|小黑屋|联系我们|仿真互动网 ( 京ICP备15048925号-7 )

GMT+8, 2024-3-29 08:38 , Processed in 0.076464 second(s), 18 queries , Gzip On, MemCache On.

Powered by Discuz! X3.5 Licensed

© 2001-2024 Discuz! Team.

快速回复 返回顶部 返回列表