该 MATLAB 代码展示了如何根据拓扑优化理论 build 一座桥梁。该代码还基于欧拉-伯努利梁理论。
最后一部分(%% 最优性标准更新)由最优性标准方法组成。不过，我对这句话理解有些困难。
xnew=max(0.001,max(x-move,min(1,min(x+move,x.*sqrt(-dc./lmid)))));
在这句话的最后，我们将比较 1、x+move 和 x.*sqrt(-dc./lmid)，并找到其中的最小值。不过，我不明白这个x.*sqrt(-dc./lmid)是什么意思，为什么突然就出来了。
谁能解释一下这个x.*sqrt(-dc./lmid)的含义，以及为什么我必须比较它？ This image might be of help.
%% Ground-structure based Topology Optimization by D.-M. Kim, May 2016
function gstop(Nx,Ny,Vlim,penal,maxL)
% Generate Nodes & Elements
[node,elem]=GS(Nx,Ny,maxL);
Ne=size(elem,1); % Total number of elements
% Initialize
x=zeros(Ne,1);
x(1:Ne)=0.1; % initial density
loop=0;
change=1;
change2=1;
c=0;
cold=0;
diffc=1;
diffc2=1;
xold=x;
filedir = 'images/Plot_iteration_';
% Start Iteration
while ((change>0.001 && change2>0.001) || (diffc>1e-2 && diffc2>1e-2))
loop=loop+1;
xold2=xold;
xold=x;
cold2=cold;
cold=c;
% FE-analysis
[U]=FE(Nx,Ny,node,elem,x,penal);
% Objective Function and Sensitivity Analysis
[KE]=stiffness(elem); % element stiffness matrix
c=0.;
dc=zeros(Ne,1);
for jj=1:Ne
n1=elem(jj,1);
n2=elem(jj,2);
Ue=U([3*n1-2;3*n1-1;3*n1;3*n2-2;3*n2-1;3*n2],1);
c=c+x(jj)^penal*Ue'*KE(:,:,jj)*Ue;
dc(jj)=-penal*x(jj)^(penal-1)*Ue'*KE(:,:,jj)*Ue;
end
% Design Update by the Optimality Criteria Method
[x]=OC(x,elem,Vlim,dc);
% Print Results
change=max(abs(x-xold));
change2=max(abs(x-xold2));
diffc=abs((c-cold)/c);
diffc2=abs((c-cold2)/c);
disp(['It.:' sprintf('%4i',loop) ' / Obj.:' sprintf('%10.4f',c)...
' / Vol.:' sprintf('%6.3f',x'*elem(:,3))...
' / ch.:' sprintf('%6.3f',change)])
% Plot Densities
clf
hold on
plot(node(:,1),node(:,2),'.','color','r','markersize',3)
[xsort,order]=sort(x);
for jj=1:Ne
if xsort(jj)>0.01
n1=elem(order(jj),1);
n2=elem(order(jj),2);
% undeformed
x1=node(n1,1); y1=node(n1,2);
x2=node(n2,1); y2=node(n2,2);
% deformed
x1=node(n1,1)+U(3*n1-2); y1=node(n1,2)+U(3*n1-1);
x2=node(n2,1)+U(3*n2-2); y2=node(n2,2)+U(3*n2-1);
r=3*xsort(jj);
rr=1-xsort(jj);
plot([x1,x2],[y1,y2],'linewidth',r,'color',rr*[0 0 0])
% plotting the mirror image
hold on
plot([-x1, -x2], [y1, y2],'linewidth',r,'color',rr*[1 1 1])
end
end
axis equal; axis tight; axis off; pause(1e-6);
print(strcat(filedir, int2str(loop)), '-dpng');
end
%% Generate Ground-structure (Nodes & Elements)
function [node,elem]=GS(Nx,Ny,maxL)
Nn=Nx*Ny; % Total Number of nodes
node=zeros(Nn,2); % Location of nodes (col1: x-dir / col2: y-dir)
Lx=80; Ly=40;
% Generate Nodes
no=0;
for jy=1:Ny
for jx=1:Nx
no=no+1;
node(no,1)=(Lx/(Nx-1))*(jx-1); % x-dir. position
node(no,2)=(Ly/(Ny-1))*(jy-1); % y-dir. position
end
end
% Generate Elements
elem=zeros(10000,4); % Element data (node 1 / node 2 / elem. length / rotation angle)
no=0;
for n1=1:Nn
no_temp=0;
the_temp=[];
for n2=n1+1:Nn
x1=node(n1,1); y1=node(n1,2);
x2=node(n2,1); y2=node(n2,2);
dx=x2-x1; dy=y2-y1;
E_len=norm([dx,dy]);
E_theta=acos(dx/E_len);
max_length=maxL*(Lx/(Nx-1));
% Save element data only if the length is shorter than maxL
if (E_len<=max_length)&&(no_temp==0||isempty(find(abs(the_temp-E_theta)<1e-3,1)))
no=no+1;
elem(no,:)=[n1,n2,E_len,E_theta];
no_temp=no_temp+1;
the_temp(no_temp,1)=E_theta;
end
end
end
elem=elem(1:no,:);
%% FE-analysis
function [U]=FE(Nx,Ny,node,elem,x,penal)
[KE]=stiffness(elem); % element stiffness matrix
Nn=size(node,1); % Total number of nodes
Ne=size(elem,1); % Total number of elements
K=zeros(3*Nn,3*Nn); % global stiffness matrix
F=zeros(3*Nn,1); % load vector
U=zeros(3*Nn,1); % displacement vector
for jj=1:Ne
n1=elem(jj,1);
n2=elem(jj,2);
edof=[3*n1-2;3*n1-1;3*n1;3*n2-2;3*n2-1;3*n2];
K(edof,edof)=K(edof,edof)+x(jj)^penal*KE(:,:,jj);
end
% Define Loads and Supports
F(2,1)=-10; %kN
fixeddofs=union([1:3*Nx:3*Nx*(Ny-1)+1],[3:3*Nx:3*Nx*(Ny-1)+3]);
fixeddofs=union([3*Nx-1],fixeddofs);
alldofs=[1:3*Nn];
freedofs=setdiff(alldofs,fixeddofs);
% Solving
U(freedofs,:)=K(freedofs,freedofs)\F(freedofs,:);
U(fixeddofs,:)=0;
%% Stiffness matrix
function [KE]=stiffness(elem)
E=10000;
b=0.2;
h=0.4;
inertia=b*h^3/12;
csarea=b*h;
Ne=size(elem,1); % Total number of elements
KE1=zeros(6,6,Ne);
KE2=zeros(6,6,Ne);
KE=zeros(6,6,Ne);
KEL=zeros(6,6,Ne);
for jj=1:Ne
EL=elem(jj,3);
KE1(:,:,jj)=inertia*E/(EL^3)*...
[0 0 0 0 0 0;
0 12 6*EL 0 -12 6*EL;
0 6*EL 4*EL^2 0 -6*EL 2*EL^2;
0 0 0 0 0 0;
0 -12 -6*EL 0 12 -6*EL;
0 6*EL 2*EL^2 0 -6*EL 4*EL^2];
KE2(:,:,jj)=csarea*E/EL*...
[1 0 0 -1 0 0;
0 0 0 0 0 0;
0 0 0 0 0 0;
-1 0 0 1 0 0;
0 0 0 0 0 0;
0 0 0 0 0 0];
cs=cos(elem(jj,4));
si=sin(elem(jj,4));
T=[cs si 0 0 0 0;
-si cs 0 0 0 0;
0 0 1 0 0 0;
0 0 0 cs si 0;
0 0 0 -si cs 0;
0 0 0 0 0 1];
KEL(:,:,jj)=KE1(:,:,jj)+KE2(:,:,jj);
KE(:,:,jj)=T'*KEL(:,:,jj)*T;
end
%% Optimality Criteria Update
function [xnew]=OC(x,elem,Vlim,dc)
l1=0; l2=100000; move=0.1;
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)))));
vol=xnew'*elem(:,3);
if vol-Vlim>0
l1=lmid;
else
l2=lmid;
end
end
最佳答案
TL;DR:直接回答问题
what this x.*sqrt(-dc./lmid) means, and why do I have to compare it?
(1) x.*sqrt(-dc./lmid) 为 x·Beη，更新后的值xenew。
(2) 您必须对其进行比较，以确保它在适当的区间内(大于或等于 xmin 和 xe-m，且小于或等于 1 和 xe+m)。
说明
什么是Beη
根据最优性标准，您希望将之前的 xe 估计更新为更准确的值。您可以通过将xe乘以因子q来实现这一点:
x<sub>e</sub><sup>new</sup> = x<sub>e</sub>·q
- 如果q恰好是1，那么xe已经在a中(可能是 局部)最优。
- 如果q > 1，那么您低估了xe。
- 如果q < 1，那么您高估了xe。
在您的示例中，q 是 Beη。
η的含义
首先要注意的是，η是一个调整参数，即它用于帮助方法收敛。之所以称为“阻尼系数”，是因为它限制了 Beη 的大小(注意 0 < η <= 1)。在您的示例中，根据经验将其设置为 0.5，即平方根。
为什么要进行比较
xe 限制为 xmin <= xe <= 1。因此，我们必须确保 xenew 是否小于 xmin 或大于 1，我们将截断它。作为引用，此操作通常称为 clamping 。此外，我们还希望将 xenew 限制为有点接近 xe(帮助收敛)。因此，我们需要再次限制 xenew 使其比 m 值更接近。两种比较都将确保 xenew 有效(xmin <= x e·Beη <= 1) 并接近之前的值 (xe - m <= xe·Beη <= xe + m)。
什么是 sqrt(-dc./lmid)
它是Beη。查看您添加的图像，我们可以看到 sqrt 来自 η，-dc 是 -∂c/∂x e 和 lmid 来自 λ·∂V/∂xe。在这种情况下，元素的体积始终相同，dV = 1，因此它脱离了方程。因此，lmid 只是λ。它被称为“mid”，因为为了找到它，代码运行 bisection algorithm ，每次迭代将间隔除以一半。
关于matlab - 最优性准则法，我们在Stack Overflow上找到一个类似的问题： https://stackoverflow.com/questions/47714378/