You can not select more than 25 topics
Topics must start with a letter or number, can include dashes ('-') and can be up to 35 characters long.
190 lines
6.7 KiB
190 lines
6.7 KiB
%本文给定top3d(30,20,10,0.4,3,1.5)
|
|
function top3d(nelx,nely,nelz,volfrac,penal,rmin)
|
|
% 定义循环参数
|
|
maxloop = 200; % 最大迭代次数
|
|
tolx = 0.01; % 终止条件（残差）
|
|
displayflag = 1; % 显示结构表示
|
|
% 材料的属性
|
|
E0 = 1; % 固体区域的杨氏模量
|
|
Emin = 1e-9; % 非固体区域的杨氏模量，尽可能小但是为了避免奇异性一般不取0
|
|
nu = 0.3; % 泊松比
|
|
% USER-DEFINED LOAD DOFs
|
|
[il,jl,kl] = meshgrid(nelx, 0, 0:nelz); % Coordinates
|
|
loadnid = kl*(nelx+1)*(nely+1)+il*(nely+1)+(nely+1-jl); % Node IDs
|
|
loaddof = 3*loadnid(:) - 1; % DOFs
|
|
% USER-DEFINED SUPPORT FIXED DOFs
|
|
[iif,jf,kf] = meshgrid(0,0:nely,0:nelz); % Coordinates
|
|
fixednid = kf*(nelx+1)*(nely+1)+iif*(nely+1)+(nely+1-jf); % Node IDs
|
|
fixeddof = [3*fixednid(:); 3*fixednid(:)-1; 3*fixednid(:)-2]; % DOFs
|
|
% 有限元分析程序
|
|
nele = nelx*nely*nelz;
|
|
ndof = 3*(nelx+1)*(nely+1)*(nelz+1);
|
|
F = sparse(loaddof,1,-1,ndof,1);
|
|
U = zeros(ndof,1);
|
|
freedofs = setdiff(1:ndof,fixeddof);
|
|
KE = lk_H8(nu);
|
|
nodegrd = reshape(1:(nely+1)*(nelx+1),nely+1,nelx+1);
|
|
nodeids = reshape(nodegrd(1:end-1,1:end-1),nely*nelx,1);
|
|
nodeidz = 0:(nely+1)*(nelx+1):(nelz-1)*(nely+1)*(nelx+1);
|
|
nodeids = repmat(nodeids,size(nodeidz))+repmat(nodeidz,size(nodeids));
|
|
edofVec = 3*nodeids(:)+1;
|
|
edofMat = repmat(edofVec,1,24)+ ...
|
|
repmat([0 1 2 3*nely + [3 4 5 0 1 2] -3 -2 -1 ...
|
|
3*(nely+1)*(nelx+1)+[0 1 2 3*nely + [3 4 5 0 1 2] -3 -2 -1]],nele,1);
|
|
iK = reshape(kron(edofMat,ones(24,1))',24*24*nele,1);
|
|
jK = reshape(kron(edofMat,ones(1,24))',24*24*nele,1);
|
|
% 过滤器
|
|
iH = ones(nele*(2*(ceil(rmin)-1)+1)^2,1);
|
|
jH = ones(size(iH));
|
|
sH = zeros(size(iH));
|
|
k = 0;
|
|
for k1 = 1:nelz
|
|
for i1 = 1:nelx
|
|
for j1 = 1:nely
|
|
e1 = (k1-1)*nelx*nely + (i1-1)*nely+j1;
|
|
for k2 = max(k1-(ceil(rmin)-1),1):min(k1+(ceil(rmin)-1),nelz)
|
|
for i2 = max(i1-(ceil(rmin)-1),1):min(i1+(ceil(rmin)-1),nelx)
|
|
for j2 = max(j1-(ceil(rmin)-1),1):min(j1+(ceil(rmin)-1),nely)
|
|
e2 = (k2-1)*nelx*nely + (i2-1)*nely+j2;
|
|
k = k+1;
|
|
iH(k) = e1;
|
|
jH(k) = e2;
|
|
sH(k) = max(0,rmin-sqrt((i1-i2)^2+(j1-j2)^2+(k1-k2)^2));
|
|
end
|
|
end
|
|
end
|
|
end
|
|
end
|
|
end
|
|
H = sparse(iH,jH,sH);
|
|
Hs = sum(H,2);
|
|
% 迭代初始化
|
|
x = repmat(volfrac,[nely,nelx,nelz]);
|
|
xPhys = x;
|
|
loop = 0;
|
|
change = 1;
|
|
|
|
m=1;
|
|
n=nele;
|
|
xmin=zeros (n, 1) ;
|
|
xmax=ones (n, 1) ;
|
|
xold1=x(:);
|
|
xold2=x(:);
|
|
low=ones(n, 1) ;
|
|
upp=ones(n, 1) ;
|
|
a0=1;
|
|
a=zeros (m, 1) ;
|
|
c_MMA=10000*ones (m, 1);
|
|
d=zeros (m, 1) ;
|
|
% 开始迭代
|
|
while change > tolx && loop < maxloop
|
|
loop = loop+1;
|
|
% 有限元分析
|
|
sK = reshape(KE(:)*(Emin+xPhys(:)'.^penal*(E0-Emin)),24*24*nele,1);
|
|
K = sparse(iK,jK,sK); K = (K+K')/2;
|
|
U(freedofs,:) = K(freedofs,freedofs)\F(freedofs,:);
|
|
% 定义目标函数与灵敏度分析
|
|
ce = reshape(sum((U(edofMat)*KE).*U(edofMat),2),[nely,nelx,nelz]);
|
|
c = sum(sum(sum((Emin+xPhys.^penal*(E0-Emin)).*ce)));
|
|
dc = -penal*(E0-Emin)*xPhys.^(penal-1).*ce;
|
|
dv = ones(nely,nelx,nelz);
|
|
% 对灵敏度进行过滤与投影，需要注意的是，并不能数学上证明对灵敏度进行过滤的
|
|
% 稳定性，如果出现不适定的情况可以尝试更改过滤为密度过滤
|
|
dc(:) = H*(dc(:)./Hs);
|
|
dv(:) = H*(dv(:)./Hs);
|
|
% 移动渐近线法求解
|
|
xval = x(:);
|
|
f0val = c;
|
|
df0dx = dc(:);
|
|
fval = sum(xPhys(:))/(volfrac*nele) - 1;
|
|
dfdx = dv(:)' / (volfrac*nele);
|
|
[xmma, ~, ~, ~, ~, ~, ~, ~, ~, low,upp] = ...
|
|
mmasub(m, n, loop, xval, xmin, xmax, xold1, xold2, ...
|
|
f0val,df0dx,fval,dfdx,low,upp,a0,a,c_MMA,d);
|
|
% 更新迭代变量
|
|
xnew = reshape(xmma,nely,nelx,nelz);
|
|
xPhys(:) = (H*xnew(:))./Hs;
|
|
xold2 = xold1(:);
|
|
xold1 = x(:);
|
|
change = max(abs(xnew(:)-x(:)));
|
|
x = xnew;
|
|
% 输出结果
|
|
fprintf(' It.:%5i Obj.:%11.4f Vol.:%7.3f ch.:%7.3f\n',loop,c,mean(xPhys(:)),change);
|
|
% 输出密度
|
|
if displayflag, clf; display_3D(xPhys); end
|
|
end
|
|
clf; display_3D(xPhys);
|
|
end
|
|
|
|
|
|
% === 生成单元刚度矩阵 ===
|
|
function [KE] = lk_H8(nu)
|
|
A = [32 6 -8 6 -6 4 3 -6 -10 3 -3 -3 -4 -8;
|
|
-48 0 0 -24 24 0 0 0 12 -12 0 12 12 12];
|
|
k = 1/144*A'*[1; nu];
|
|
|
|
K1 = [k(1) k(2) k(2) k(3) k(5) k(5);
|
|
k(2) k(1) k(2) k(4) k(6) k(7);
|
|
k(2) k(2) k(1) k(4) k(7) k(6);
|
|
k(3) k(4) k(4) k(1) k(8) k(8);
|
|
k(5) k(6) k(7) k(8) k(1) k(2);
|
|
k(5) k(7) k(6) k(8) k(2) k(1)];
|
|
K2 = [k(9) k(8) k(12) k(6) k(4) k(7);
|
|
k(8) k(9) k(12) k(5) k(3) k(5);
|
|
k(10) k(10) k(13) k(7) k(4) k(6);
|
|
k(6) k(5) k(11) k(9) k(2) k(10);
|
|
k(4) k(3) k(5) k(2) k(9) k(12)
|
|
k(11) k(4) k(6) k(12) k(10) k(13)];
|
|
K3 = [k(6) k(7) k(4) k(9) k(12) k(8);
|
|
k(7) k(6) k(4) k(10) k(13) k(10);
|
|
k(5) k(5) k(3) k(8) k(12) k(9);
|
|
k(9) k(10) k(2) k(6) k(11) k(5);
|
|
k(12) k(13) k(10) k(11) k(6) k(4);
|
|
k(2) k(12) k(9) k(4) k(5) k(3)];
|
|
K4 = [k(14) k(11) k(11) k(13) k(10) k(10);
|
|
k(11) k(14) k(11) k(12) k(9) k(8);
|
|
k(11) k(11) k(14) k(12) k(8) k(9);
|
|
k(13) k(12) k(12) k(14) k(7) k(7);
|
|
k(10) k(9) k(8) k(7) k(14) k(11);
|
|
k(10) k(8) k(9) k(7) k(11) k(14)];
|
|
K5 = [k(1) k(2) k(8) k(3) k(5) k(4);
|
|
k(2) k(1) k(8) k(4) k(6) k(11);
|
|
k(8) k(8) k(1) k(5) k(11) k(6);
|
|
k(3) k(4) k(5) k(1) k(8) k(2);
|
|
k(5) k(6) k(11) k(8) k(1) k(8);
|
|
k(4) k(11) k(6) k(2) k(8) k(1)];
|
|
K6 = [k(14) k(11) k(7) k(13) k(10) k(12);
|
|
k(11) k(14) k(7) k(12) k(9) k(2);
|
|
k(7) k(7) k(14) k(10) k(2) k(9);
|
|
k(13) k(12) k(10) k(14) k(7) k(11);
|
|
k(10) k(9) k(2) k(7) k(14) k(7);
|
|
k(12) k(2) k(9) k(11) k(7) k(14)];
|
|
KE = 1/((nu+1)*(1-2*nu))*...
|
|
[ K1 K2 K3 K4;
|
|
K2' K5 K6 K3';
|
|
K3' K6 K5' K2';
|
|
K4 K3 K2 K1'];
|
|
end
|
|
% === 展示3D效果图 (ISO-VIEW) ===
|
|
function display_3D(rho)
|
|
[nely,nelx,nelz] = size(rho);
|
|
hx = 1; hy = 1; hz = 1; % 定义效果图的单元大小
|
|
face = [1 2 3 4; 2 6 7 3; 4 3 7 8; 1 5 8 4; 1 2 6 5; 5 6 7 8];
|
|
set(gcf,'Name','ISO display','NumberTitle','off');
|
|
for k = 1:nelz
|
|
z = (k-1)*hz;
|
|
for i = 1:nelx
|
|
x = (i-1)*hx;
|
|
for j = 1:nely
|
|
y = nely*hy - (j-1)*hy;
|
|
if (rho(j,i,k) > 0.5) % 定义展示出的密度阈值
|
|
vert = [x y z; x y-hx z; x+hx y-hx z; x+hx y z; x y z+hx;x y-hx z+hx; x+hx y-hx z+hx;x+hx y z+hx];
|
|
vert(:,[2 3]) = vert(:,[3 2]); vert(:,2,:) = -vert(:,2,:);
|
|
patch('Faces',face,'Vertices',vert,'FaceColor',[0.2+0.8*(1-rho(j,i,k)),0.2+0.8*(1-rho(j,i,k)),0.2+0.8*(1-rho(j,i,k))]);
|
|
hold on;
|
|
end
|
|
end
|
|
end
|
|
end
|
|
axis equal; axis tight; axis off; box on; view([30,30]); pause(1e-6);
|
|
end
|
|
|