topmma.zip

  • 小拓扑
    了解作者
  • matlab
    开发工具
  • 6KB
    文件大小
  • zip
    文件格式
  • 0
    收藏次数
  • 1 积分
    下载积分
  • 0
    下载次数
  • 2021-04-19 21:11
    上传日期
利用移动渐近线法即MMA法优化设计变量,适合梯度类、大规模求解
topmma.zip
  • topmma.m
    18.6KB
内容介绍
% a 99 line topology optimization code by Ole Sigmund,October 199 %整块板均匀加热,左边界给定温度T=0,其余边界为绝热边界 function topmma(nelx,nely,volfrac,penal,rmin) % initialize x(1:nely,1:nelx)=volfrac; % x是设计变量 loop=0; % 存放迭代次数的变量 change=1.; % 每次迭代,x的改变值,用来判断何时收敛 % start ineration xval=reshape(x,nely*nelx,1); xold1=xval; xold2=xold1; low=[]; upp=[]; while change>0.04 %当两次连续目标函数迭代的差<=0.01时,且迭代次数小于50次,迭代结束 loop=loop+1; xold=x; %把前一次的设计变量付给xold % FE analysis [U]=FE(nelx,nely,x,penal); %有限元分析,得到位移矢量U % objective function and sensitivity analysis 目标函数和灵敏度分析 [KE]=lk; % 单位刚度矩阵 obj=0.; % 用来存放目标函数的变量柔度,使得刚度最大,柔度最小 for ely=1:nely for elx=1:nelx n1=(nely+1)*(elx-1)+ely; %左上角的单元节点 n2=(nely+1)*elx+ely; %右上角的单元节点 Ue=U([n1;n2;n2+1;n1+1],1); %所示单元的自由度,左上,右上,右下,左下 obj=obj+(0.001+0.999*penal*x(ely,elx)^penal)*Ue'*KE*Ue; % 计算目标函数的值(柔度) dc(ely,elx)=-0.999*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 最佳优化方式设计更新 % DEFINE INITIAL CONDITIONS m=1; n=nelx*nely; iter=loop; xmin=0.001*ones(n,1); xmax=ones(n,1); f0val=obj; df0dx=reshape(dc,n,1); df0dx2=0*df0dx; % 定义约束条件 fval=sum(sum(x))/(volfrac*n)-1; dfdx=ones(1,n)/(volfrac*n); dfdx2=0*dfdx; a0=1; a=0*ones(m,1); c=10000*ones(m,1); d=0*ones(m,1); % RUN MMA [xmma,ymma,zmma,lam,xsi,eta,mu,zet,s,low,upp] = ... mmasub(m,n,iter,xval,xmin,xmax,xold1,xold2, ... f0val,df0dx,df0dx2,fval,dfdx,dfdx2,low,upp,a0,a,c,d); % END UPDATE BY MMA METHOD xold2=xold1; xold1=xval; xval=xmma; x=reshape(xval,nely,nelx); % PRINT RESULTS change=max(max(abs(x-xold))); disp(['It.:' sprintf('%4i',loop) 'Obj.:' sprintf('%10.4f',obj) '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 % axis坐标轴控制函数 axis off 坐标系设为不可见,(1e-6)10的-6次方 end % imagesc 函数显示灰度图像 下面的代码是具有两个输入参数的 imagesc 函数显示一副灰度图像 imagesc(1,[0,1]); colormap(gray); imagesc 函数中的第二个参数确定灰度范围。 % FE analysis 有限元求解子程序 function [U]=FE(nelx,nely,x,penal) %自定义函数,最后返回[U] [KE]=lk; % sparse 把一个全矩阵转化为一个稀疏矩阵,只存储每一个非零元素的三个值:元素值,元素的行号和列号 % 总体刚度矩阵的稀疏矩阵 % *2是因为x,y都有一个数 K=sparse((nelx+1)*(nely+1),(nelx+1)*(nely+1)); %k表示整体刚度矩阵 F=sparse((nely+1)*(nelx+1),1); %F表示力的向量 U=sparse((nely+1)*(nelx+1),1); %U表示整体变形 for elx=1:nelx for ely=1:nely %一列列的排序,y*x n1=(nely+1)*(elx-1)+ely; %左上 n2=(nely+1)*elx +ely; %右上 % 左上,右上,右下,左下 自由度 % 一个点有两个,所以要*2。第一个从1开始,所以*2之后要-1。 edof=[n1;n2;n2+1;n1+1]; K(edof,edof)=K(edof,edof)+(0.001+0.999*x(ely,elx)^penal)*KE; % 将单元刚度矩阵组装成总的刚度矩阵 end end % define loads and supports 定义荷载和支撑(半MBB梁) F(:,1)=0.01; %施加载荷 垂直单元力 %施加位移约束 fixeddofs =[nely/2+1-((nely/20)):2:nely/2+1+(nely/20)]; %固定节点 alldofs =[1:(nely+1)*(nelx+1)]; %固定所有节点 % 因无约束自由度与固定自由度的不同来找到无约束自由度 freedofs =setdiff(alldofs,fixeddofs); % solving 求解线性方程组,得到节点自由度的位移值 U(freedofs,:)=K(freedofs,freedofs)\F(freedofs,:); U(fixeddofs,:)=0; %受约束的节点自由度的位移值设为0 end % mesh-independency filter 独立网格过滤 敏度过滤技术子程序 function [dcn]=check(nelx,nely,rmin,x,dc) %dcn 清零,dcn用来保存更新的目标函数灵敏度 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) %floor()函数时取整的意思 for l=max(j-floor(rmin),1):min(j+floor(rmin),nely) fac=rmin-sqrt((i-k)^2+(j-l)^2); %卷积算子(重量系数)Hf=过滤半径-控制函数(两元素中心距离) 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 end % Element stiffness matrix 单元刚度矩阵 function [KE]=lk KE=[2/3 -1/6 -1/3 -1/6 -1/6 2/3 -1/6 -1/3 -1/3 -1/6 2/3 -1/6 -1/6 -1/3 -1/6 2/3]; end % This is the file mmasub.m % function [xmma,ymma,zmma,lam,xsi,eta,mu,zet,s,low,upp] = ... mmasub(m,n,iter,xval,xmin,xmax,xold1,xold2, ... f0val,df0dx,df0dx2,fval,dfdx,dfdx2,low,upp,a0,a,c,d); % % Written in May 1999 by % Krister Svanberg <krille@math.kth.se> % Department of Mathematics % SE-10044 Stockholm, Sweden. % % Modified ("spdiags" instead of "diag") April 2002 % % % This function mmasub performs one MMA-iteration, aimed at % solving the nonlinear programming problem:(这个函数mmasub执行一个mma迭代,目标是求解非线性规划问题:) % % Minimize f_0(x) + a_0*z + sum( c_i*y_i + 0.5*d_i*(y_i)^2 ) % subject to f_i(x) - a_i*z - y_i <= 0, i = 1,...,m % xmin_j <= x_j <= xmax_j, j = 1,...,n % z >= 0, y_i >= 0, i = 1,...,m %*** INPUT: % % m = The number of general constraints.(一般约束的数量) % n = The number of variables x_j.(变量的个数xj) % iter = Current iteration number ( =1 the first time mmasub is called).(当前迭代数(第一次调用mmasub时=1)。) % xval = Column vector with the current values of the variables x_j. ( xval =带有变量x_j的当前值的列向量。) % xmin = Column vector with the lower bounds for the variables x_j. (xmin =具有变量x_j下界的列向量。) % xmax = Column vector with the upper bounds for the variables x_j. ( xmax =变量x_j的上界列向量。) % xold1 = xval, one iteration ago (provided that iter>1). ( xold1 = xval,一次迭代前(假设iter>1)。) % xold2 = xval, two iterations ago (provided that iter>2). (xold2 = xval,两次迭代前(假设iter>2)。) % f0val = The value of the objective function f_0 at xval. (f0val =目标函数在xval处的值f_0。) % df0dx = Column vector with the derivatives of the objective function (目标函数导数的列向量f_0对于变量x_j,在xval处计算。) % f_0 with respect to the variables x_j, calculated at xval. % df0dx2 = Column vector with the non-mixed second derivatives of the (非混合二阶导数的列向量,目标函数f_0对于变量x_j按xval计算。df0dx2(j) =二阶导数 % objective function f_0 with respect to the variables x_j,f_0相对于x_j的%(2次)。重要提示:如果有二阶导数,)让df0dx2 = 0*df0dx。) % calculated at xval. df0dx2(j) = the second derivative % of f_0 with respect to x_j (twice). % Important note: If second derivatives are not available, % simply let df0dx2 = 0*df0dx. % fval = Column vector with the values of the constraint functions f_i, (约束函数f_i的值的列向量,按xval计算。) % calculated at xval. % dfdx = (m x n)-matrix
评论
    相关推荐
    • 88行matlab拓扑优化代码-TOSSE:高效的51行Matlab代码用于拓扑优化
      88行matlab拓扑优化代码托斯 高效的51行Matlab代码,用于拓扑优化。 TOSSE(相同尺寸元素的拓扑优化)是用于2D和3D拓扑设计问题的Matlab代码。 该代码使用称为TOP88的经典88行代码作为基础,以开发一种硬0-1进化算法...
    • matlab拓扑优化代码-UNVARTOP:使用UNsmoothVARiationalTopologyOPtimization
      使用UNVARTOP方法进行2D拓扑优化的Matlab代码(用于教育目的) 目录 关于该项目 这是该文件的代码:D。Yago,J。Cante,O。Lloberas-Valls和J. Oliver。 使用非平滑变分拓扑优化(UNVARTOP)方法进行拓扑优化:...
    • 88行matlab拓扑优化代码-TopOpt:使用Matlab进行拓扑优化设计
      88行matlab拓扑优化代码TopOpt 用Matlab编写的99行拓扑优化代码 使用88行代码在MATLAB中进行高效的拓扑优化本页上的Matlab代码旨在用于工程教育。 拓扑优化领域的学生和新手可以在此处找到代码并下载。 可以在结构...
    • matlab拓扑优化代码-topopt:拓扑优化
      matlab模拟优化代码拓扑优化 这是Ole Sigmund中描述的实现。 以或开头。 在下可以找到具有更多安装要求的更有效的版本。
    • 拓扑优化过滤器.rar
      拓扑优化三角形四边形单元过滤器,可以完成敏度和密度过滤
    • 88行matlab拓扑优化代码-topopt_jl:使用Julia进行2D拓扑优化
      88行matlab拓扑优化代码使用Julia进行拓扑优化 使用Julia()进行二维拓扑优化的代码。 该代码在Julia中转换了88行的Matlab代码(),以进行二维拓扑优化。 本文讨论了Matlab实现的详细信息: “使用88行代码在...
    • 88行matlab拓扑优化代码-Topopt:Julia的拓扑优化
      88行matlab拓扑优化代码拓普 Julia的拓扑优化 我想测试Julia在拓扑优化中的性能。 我使用经典的88行matlab代码作为比较的基础。 要运行matlabe代码,只需下载top88.m代码并在MATLAB中使用以下命令运行: top88(60,...
    • 拓扑优化matlab实例,
      对于拓扑优化的matlab实例,可以应用于分析,快捷便利。
    • matlab拓扑优化代码
      matlab拓扑优化资源共享: 1)SIMP/BESO/LSM/ESO/ICM/HM等拓扑优化程序; 2)柔度拓扑、频率拓扑、应力拓扑、疲劳拓扑、解耦拓扑、流体拓扑、电磁拓扑、压电拓扑、多材料拓扑、多尺度拓扑、跨尺度拓扑、多目标拓扑、...
    • 拓扑优化matlab程序
      基于matlab结构拓扑优化的71行程序代码,迭代步骤缩小,收敛更快,优化图像更清晰 基于matlab结构拓扑优化的71行程序代码,迭代步骤缩小,收敛更快,优化图像更清晰