4 changed files with 534 additions and 0 deletions
Binary file not shown.
@ -0,0 +1,133 @@ |
|||||
|
%------------------------------------------------------- |
||||
|
% MMA求解程序 MMA.m |
||||
|
% 本代码将原本的非线性问题线性近似，非凸问题凸近似，得到一系列 |
||||
|
% 线性的、凸的子问题，再通过subsolv.m迭代求解这些线性的、凸的子问题最终找到最优值点 |
||||
|
function [xmma,ymma,zmma,lam,xsi,eta,mu,zet,s,low,upp] = ... |
||||
|
mmasub(m,n,iter,xval,xmin,xmax,xold1,xold2, ... |
||||
|
f0val,df0dx,fval,dfdx,low,upp,a0,a,c,d) |
||||
|
% This function mmasub performs one MMA-iteration, aimed at |
||||
|
% solving the nonlinear programming problem: |
||||
|
% 该函数进行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 |
||||
|
% 输入参数有： |
||||
|
% m = 约束条件的数量 |
||||
|
% n = 变量的个数 |
||||
|
% iter = 迭代步数 |
||||
|
% xval = 包含所有xj的向量 |
||||
|
% xmin = 包含所有xj取值下限的向量 |
||||
|
% xmax = 包含所有xj取值上限的向量 |
||||
|
% xold1 = 上次迭代产生的x_val |
||||
|
% xold2 = 上上次迭代产生的的xval |
||||
|
% f0val = 目标函数的现有值 |
||||
|
% df0dx = 包含目标函数对所有xj偏导数的列向量 |
||||
|
% fval = 时约束函数fi值的列向量 |
||||
|
% dfdx = x_j.fi对xj偏导数的矩阵 |
||||
|
% low = 前一次迭代的渐近线下界 |
||||
|
% upp = 前一次迭代的渐进线上界 |
||||
|
% 引入下面四个量防止在初始值选择不恰当的形况下出现子问题不适定（没有结）的情况 |
||||
|
% a0 = a_0*z中的常系数a_0 |
||||
|
% a = 用来存储a_i*z中系数 a_i的列向量 |
||||
|
% c = 用来存储c_i*y_i中系数c_i的列向量 |
||||
|
% d =用来存储0.5*d_i*(y_i)^2中系数d_i的列向量，d_i应该是一个足够大的数， |
||||
|
% 可以证明当d_i足够大时y_i趋近于0 |
||||
|
%输出的量: |
||||
|
% xmma = 子问题求解出的xj的列向量 |
||||
|
% ymma = 子问题求解出的yi的列向量 |
||||
|
% zmma = z的标量 |
||||
|
% lam = m个约束对应的拉格朗日乘子 |
||||
|
% xsi = n个α<x_j对应的拉格朗日乘子 |
||||
|
% eta = n个x_j<β对应的拉格朗日乘子 |
||||
|
% mu = m个y_i<0对应的拉格朗日乘子。 |
||||
|
% zet = z>0对应的拉格朗日乘子 |
||||
|
% s = m个约束条件对应的松弛因子 |
||||
|
% low = 本次求解的子问题对应的渐近线的下界 |
||||
|
% upp = 本次求解的子问题对应的渐近线上界 |
||||
|
% |
||||
|
%epsimin = sqrt(m+n)*10^(-9);计算机精度的最小值 |
||||
|
epsimin = 10^(-7); |
||||
|
raa0 = 0.00001; |
||||
|
move = 1.0; |
||||
|
albefa = 0.1; |
||||
|
asyinit = 0.5; |
||||
|
asyincr = 1.2; |
||||
|
asydecr = 0.7; |
||||
|
eeen = ones(n,1); |
||||
|
eeem = ones(m,1); |
||||
|
zeron = zeros(n,1); |
||||
|
|
||||
|
% Calculation of the asymptotes low and upp :计算渐近线的上下界 |
||||
|
if iter < 2.5 |
||||
|
low = xval - asyinit*(xmax-xmin); |
||||
|
upp = xval + asyinit*(xmax-xmin); |
||||
|
else |
||||
|
zzz = (xval-xold1).*(xold1-xold2); |
||||
|
factor = eeen; |
||||
|
factor(find(zzz > 0)) = asyincr; |
||||
|
factor(find(zzz < 0)) = asydecr; |
||||
|
low = xval - factor.*(xold1 - low); |
||||
|
upp = xval + factor.*(upp - xold1); |
||||
|
lowmin = xval - 10*(xmax-xmin); |
||||
|
lowmax = xval - 0.01*(xmax-xmin); |
||||
|
uppmin = xval + 0.01*(xmax-xmin); |
||||
|
uppmax = xval + 10*(xmax-xmin); |
||||
|
low = max(low,lowmin); |
||||
|
low = min(low,lowmax); |
||||
|
upp = min(upp,uppmax); |
||||
|
upp = max(upp,uppmin); |
||||
|
end |
||||
|
|
||||
|
% 计算设计变量的精确界限α和β |
||||
|
|
||||
|
zzz1 = low + albefa*(xval-low); |
||||
|
zzz2 = xval - move*(xmax-xmin); |
||||
|
zzz = max(zzz1,zzz2); |
||||
|
alfa = max(zzz,xmin); |
||||
|
zzz1 = upp - albefa*(upp-xval); |
||||
|
zzz2 = xval + move*(xmax-xmin); |
||||
|
zzz = min(zzz1,zzz2); |
||||
|
beta = min(zzz,xmax); |
||||
|
|
||||
|
% 计算 p0, q0, P, Q 和 b. |
||||
|
|
||||
|
xmami = xmax-xmin; |
||||
|
xmamieps = 0.00001*eeen; |
||||
|
xmami = max(xmami,xmamieps); |
||||
|
xmamiinv = eeen./xmami; |
||||
|
ux1 = upp-xval; |
||||
|
ux2 = ux1.*ux1; |
||||
|
xl1 = xval-low; |
||||
|
xl2 = xl1.*xl1; |
||||
|
uxinv = eeen./ux1; |
||||
|
xlinv = eeen./xl1; |
||||
|
% |
||||
|
p0 = zeron; |
||||
|
q0 = zeron; |
||||
|
p0 = max(df0dx,0); |
||||
|
q0 = max(-df0dx,0); |
||||
|
%p0(find(df0dx > 0)) = df0dx(find(df0dx > 0)); |
||||
|
%q0(find(df0dx < 0)) = -df0dx(find(df0dx < 0)); |
||||
|
pq0 = 0.001*(p0 + q0) + raa0*xmamiinv; |
||||
|
p0 = p0 + pq0; |
||||
|
q0 = q0 + pq0; |
||||
|
p0 = p0.*ux2; |
||||
|
q0 = q0.*xl2; |
||||
|
% |
||||
|
P = sparse(m,n); |
||||
|
Q = sparse(m,n); |
||||
|
P = max(dfdx,0); |
||||
|
Q = max(-dfdx,0); |
||||
|
%P(find(dfdx > 0)) = dfdx(find(dfdx > 0)); |
||||
|
%Q(find(dfdx < 0)) = -dfdx(find(dfdx < 0)); |
||||
|
PQ = 0.001*(P + Q) + raa0*eeem*xmamiinv'; |
||||
|
P = P + PQ; |
||||
|
Q = Q + PQ; |
||||
|
P = P * spdiags(ux2,0,n,n); |
||||
|
Q = Q * spdiags(xl2,0,n,n); |
||||
|
b = P*uxinv + Q*xlinv - fval ; |
||||
|
% |
||||
|
%%% 使用反演原对偶牛顿方法求解子问题 |
||||
|
[xmma,ymma,zmma,lam,xsi,eta,mu,zet,s] = ... |
||||
|
subsolv(m,n,epsimin,low,upp,alfa,beta,p0,q0,P,Q,a0,a,b,c,d); |
||||
@ -0,0 +1,211 @@ |
|||||
|
% 子问题求解程序 subsolv.m |
||||
|
% |
||||
|
function [xmma,ymma,zmma,lamma,xsimma,etamma,mumma,zetmma,smma] = ... |
||||
|
subsolv(m,n,epsimin,low,upp,alfa,beta,p0,q0,P,Q,a0,a,b,c,d); |
||||
|
% |
||||
|
% 使用该函数求解MMA生成的子问题可以写成： |
||||
|
% |
||||
|
% minimize SUM[ p0j/(uppj-xj) + q0j/(xj-lowj) ] + a0*z + |
||||
|
% + SUM[ ci*yi + 0.5*di*(yi)^2 ], |
||||
|
% |
||||
|
% subject to SUM[ pij/(uppj-xj) + qij/(xj-lowj) ] - ai*z - yi <= bi, |
||||
|
% alfaj <= xj <= betaj, yi >= 0, z >= 0. |
||||
|
% 然后使用拉格朗日的方法求解子问题的对偶问题 |
||||
|
% Input: m, n, low, upp, alfa, beta, p0, q0, P, Q, a0, a, b, c, d. |
||||
|
% Output: xmma,ymma,zmma, slack variables and Lagrange multiplers. |
||||
|
% |
||||
|
een = ones(n,1); |
||||
|
eem = ones(m,1); |
||||
|
epsi = 1; |
||||
|
epsvecn = epsi*een; |
||||
|
epsvecm = epsi*eem; |
||||
|
x = 0.5*(alfa+beta); |
||||
|
y = eem; |
||||
|
z = 1; |
||||
|
lam = eem; |
||||
|
xsi = een./(x-alfa); |
||||
|
xsi = max(xsi,een); |
||||
|
eta = een./(beta-x); |
||||
|
eta = max(eta,een); |
||||
|
mu = max(eem,0.5*c); |
||||
|
zet = 1; |
||||
|
s = eem; |
||||
|
itera = 0; |
||||
|
|
||||
|
while epsi > epsimin |
||||
|
epsvecn = epsi*een; |
||||
|
epsvecm = epsi*eem; |
||||
|
ux1 = upp-x; |
||||
|
xl1 = x-low; |
||||
|
ux2 = ux1.*ux1; |
||||
|
xl2 = xl1.*xl1; |
||||
|
uxinv1 = een./ux1; |
||||
|
xlinv1 = een./xl1; |
||||
|
|
||||
|
plam = p0 + P'*lam ; |
||||
|
qlam = q0 + Q'*lam ; |
||||
|
gvec = P*uxinv1 + Q*xlinv1; |
||||
|
dpsidx = plam./ux2 - qlam./xl2 ; |
||||
|
|
||||
|
rex = dpsidx - xsi + eta; |
||||
|
rey = c + d.*y - mu - lam; |
||||
|
rez = a0 - zet - a'*lam; |
||||
|
relam = gvec - a*z - y + s - b; |
||||
|
rexsi = xsi.*(x-alfa) - epsvecn; |
||||
|
reeta = eta.*(beta-x) - epsvecn; |
||||
|
remu = mu.*y - epsvecm; |
||||
|
rezet = zet*z - epsi; |
||||
|
res = lam.*s - epsvecm; |
||||
|
|
||||
|
residu1 = [rex' rey' rez]'; |
||||
|
residu2 = [relam' rexsi' reeta' remu' rezet res']'; |
||||
|
residu = [residu1' residu2']'; |
||||
|
residunorm = sqrt(residu'*residu); |
||||
|
residumax = max(abs(residu)); |
||||
|
|
||||
|
ittt = 0; |
||||
|
while residumax > 0.9*epsi & ittt < 100 |
||||
|
ittt=ittt + 1; |
||||
|
itera=itera + 1; |
||||
|
|
||||
|
ux1 = upp-x; |
||||
|
xl1 = x-low; |
||||
|
ux2 = ux1.*ux1; |
||||
|
xl2 = xl1.*xl1; |
||||
|
ux3 = ux1.*ux2; |
||||
|
xl3 = xl1.*xl2; |
||||
|
uxinv1 = een./ux1; |
||||
|
xlinv1 = een./xl1; |
||||
|
uxinv2 = een./ux2; |
||||
|
xlinv2 = een./xl2; |
||||
|
plam = p0 + P'*lam ; |
||||
|
qlam = q0 + Q'*lam ; |
||||
|
gvec = P*uxinv1 + Q*xlinv1; |
||||
|
GG = P*spdiags(uxinv2,0,n,n) - Q*spdiags(xlinv2,0,n,n); |
||||
|
dpsidx = plam./ux2 - qlam./xl2 ; |
||||
|
delx = dpsidx - epsvecn./(x-alfa) + epsvecn./(beta-x); |
||||
|
dely = c + d.*y - lam - epsvecm./y; |
||||
|
delz = a0 - a'*lam - epsi/z; |
||||
|
dellam = gvec - a*z - y - b + epsvecm./lam; |
||||
|
diagx = plam./ux3 + qlam./xl3; |
||||
|
diagx = 2*diagx + xsi./(x-alfa) + eta./(beta-x); |
||||
|
diagxinv = een./diagx; |
||||
|
diagy = d + mu./y; |
||||
|
diagyinv = eem./diagy; |
||||
|
diaglam = s./lam; |
||||
|
diaglamyi = diaglam+diagyinv; |
||||
|
|
||||
|
if m < n |
||||
|
blam = dellam + dely./diagy - GG*(delx./diagx); |
||||
|
bb = [blam' delz]'; |
||||
|
Alam = spdiags(diaglamyi,0,m,m) + GG*spdiags(diagxinv,0,n,n)*GG'; |
||||
|
AA = [Alam a |
||||
|
a' -zet/z ]; |
||||
|
solut = AA\bb; |
||||
|
dlam = solut(1:m); |
||||
|
dz = solut(m+1); |
||||
|
dx = -delx./diagx - (GG'*dlam)./diagx; |
||||
|
else |
||||
|
diaglamyiinv = eem./diaglamyi; |
||||
|
dellamyi = dellam + dely./diagy; |
||||
|
Axx = spdiags(diagx,0,n,n) + GG'*spdiags(diaglamyiinv,0,m,m)*GG; |
||||
|
azz = zet/z + a'*(a./diaglamyi); |
||||
|
axz = -GG'*(a./diaglamyi); |
||||
|
bx = delx + GG'*(dellamyi./diaglamyi); |
||||
|
bz = delz - a'*(dellamyi./diaglamyi); |
||||
|
AA = [Axx axz |
||||
|
axz' azz ]; |
||||
|
bb = [-bx' -bz]'; |
||||
|
solut = AA\bb; |
||||
|
dx = solut(1:n); |
||||
|
dz = solut(n+1); |
||||
|
dlam = (GG*dx)./diaglamyi - dz*(a./diaglamyi) + dellamyi./diaglamyi; |
||||
|
end |
||||
|
|
||||
|
dy = -dely./diagy + dlam./diagy; |
||||
|
dxsi = -xsi + epsvecn./(x-alfa) - (xsi.*dx)./(x-alfa); |
||||
|
deta = -eta + epsvecn./(beta-x) + (eta.*dx)./(beta-x); |
||||
|
dmu = -mu + epsvecm./y - (mu.*dy)./y; |
||||
|
dzet = -zet + epsi/z - zet*dz/z; |
||||
|
ds = -s + epsvecm./lam - (s.*dlam)./lam; |
||||
|
xx = [ y' z lam' xsi' eta' mu' zet s']'; |
||||
|
dxx = [dy' dz dlam' dxsi' deta' dmu' dzet ds']'; |
||||
|
|
||||
|
stepxx = -1.01*dxx./xx; |
||||
|
stmxx = max(stepxx); |
||||
|
stepalfa = -1.01*dx./(x-alfa); |
||||
|
stmalfa = max(stepalfa); |
||||
|
stepbeta = 1.01*dx./(beta-x); |
||||
|
stmbeta = max(stepbeta); |
||||
|
stmalbe = max(stmalfa,stmbeta); |
||||
|
stmalbexx = max(stmalbe,stmxx); |
||||
|
stminv = max(stmalbexx,1); |
||||
|
steg = 1/stminv; |
||||
|
|
||||
|
xold = x; |
||||
|
yold = y; |
||||
|
zold = z; |
||||
|
lamold = lam; |
||||
|
xsiold = xsi; |
||||
|
etaold = eta; |
||||
|
muold = mu; |
||||
|
zetold = zet; |
||||
|
sold = s; |
||||
|
|
||||
|
itto = 0; |
||||
|
resinew = 2*residunorm; |
||||
|
while resinew > residunorm & itto < 50 |
||||
|
itto = itto+1; |
||||
|
|
||||
|
x = xold + steg*dx; |
||||
|
y = yold + steg*dy; |
||||
|
z = zold + steg*dz; |
||||
|
lam = lamold + steg*dlam; |
||||
|
xsi = xsiold + steg*dxsi; |
||||
|
eta = etaold + steg*deta; |
||||
|
mu = muold + steg*dmu; |
||||
|
zet = zetold + steg*dzet; |
||||
|
s = sold + steg*ds; |
||||
|
ux1 = upp-x; |
||||
|
xl1 = x-low; |
||||
|
ux2 = ux1.*ux1; |
||||
|
xl2 = xl1.*xl1; |
||||
|
uxinv1 = een./ux1; |
||||
|
xlinv1 = een./xl1; |
||||
|
plam = p0 + P'*lam ; |
||||
|
qlam = q0 + Q'*lam ; |
||||
|
gvec = P*uxinv1 + Q*xlinv1; |
||||
|
dpsidx = plam./ux2 - qlam./xl2 ; |
||||
|
|
||||
|
rex = dpsidx - xsi + eta; |
||||
|
rey = c + d.*y - mu - lam; |
||||
|
rez = a0 - zet - a'*lam; |
||||
|
relam = gvec - a*z - y + s - b; |
||||
|
rexsi = xsi.*(x-alfa) - epsvecn; |
||||
|
reeta = eta.*(beta-x) - epsvecn; |
||||
|
remu = mu.*y - epsvecm; |
||||
|
rezet = zet*z - epsi; |
||||
|
res = lam.*s - epsvecm; |
||||
|
|
||||
|
residu1 = [rex' rey' rez]'; |
||||
|
residu2 = [relam' rexsi' reeta' remu' rezet res']'; |
||||
|
residu = [residu1' residu2']'; |
||||
|
resinew = sqrt(residu'*residu); |
||||
|
steg = steg/2; |
||||
|
end |
||||
|
residunorm=resinew; |
||||
|
residumax = max(abs(residu)); |
||||
|
steg = 2*steg; |
||||
|
end |
||||
|
epsi = 0.1*epsi; |
||||
|
end |
||||
|
|
||||
|
xmma = x; |
||||
|
ymma = y; |
||||
|
zmma = z; |
||||
|
lamma = lam; |
||||
|
xsimma = xsi; |
||||
|
etamma = eta; |
||||
|
mumma = mu; |
||||
|
zetmma = zet; |
||||
|
smma = s; |
||||
@ -0,0 +1,190 @@ |
|||||
|
%本文给定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 |
||||
Loading…
Reference in new issue