关于拓扑优化99行的问题？

我看了Sigmund的99行代码，想改变的一下受力但是改了几次出来的都是奇怪的图像，下面是我摘录的注释（来自tiiaan）

需要大佬带带

我的想法就是改变205-209行的代码，我想把他改成左右两端下方固定，然后力施加在中间（Y方向向下的力）最后形成如下的结果（大致相似）

就是我狂加感叹号那几行！！！！！！！！！！！！！！！！！！！！！！

%%%%%%%% A 99 LINE TOPOLOGY OPTIMIZATION CODE BY OLE SIGMUND, JANUARY 2000 %%%%%%%%
%%%%%%%% COMMENTED - OUT 1.0 BY HAOTIAN_W, JULY 2020 %%%%%%%%
function top(nelx,nely,volfrac,penal,rmin)
nelx=31;
nely=10;
volfrac=0.5;
penal=3.0;
rmin=1.5;
%===================================================================================
% nelx : 水平方向上的离散单元数;
% nely : 竖直方向上的离散单元数;
%
% volfrac : 容积率，材料体积与设计域体积之比，对应的工程问题就是"将结构减重到百分之多少";
%
% penal : 惩罚因子，SIMP方法是在0-1离散模型中引入连续变量x、系数p及中间密度单元，从而将离
% 散型优化问题转换成连续型优化问题，并且令0≤x≤1，p为惩罚因子，通过设定p>1对中间密
% 度单元进行有限度的惩罚，尽量减少中间密度单元数目，使单元密度尽可能趋于0或1;
%
% 合理选择惩罚因子的取值，可以消除多孔材料，从而得到理想的拓扑优化结果:
% 当penal<=2时 存在大量多孔材料，计算结果没有可制造性;
% 当penal>=3.5时 最终拓扑结果没有大的改变;
% 当penal>=4时 结构总体柔度的变化非常缓慢，迭代步数增加，计算时间延长;
%
% rmin : 敏度过滤半径，防止出现棋盘格现象;
%===================================================================================
% 结构优化的数学模型常用名词:
% 1) 设计变量 在设计中可调整的、变化的基本参数，本算例是单元密度;
% 2) 目标函数 设计变量的函数，优化设计的目标，本算例是柔度最小;
% 3) 约束条件 几何、应力、位移等，难点在建立约束方程，本算例是几何体积约束;
% 4) 终止准则 结束迭代的条件，本算例是目标函数变化量<=0.010;
% 5) 载荷工况 定义结构所有可能的受力情况，本算例是单一载荷;
%===================================================================================
% 基于SIMP理论的优化准则法迭代分析流程:
% 1) 定义设计域，选择合适的设计变量、目标函数以及约束函数等其他边界条件;
% 2) 结构离散为有限元网格，计算优化前的单元刚度矩阵;
% 3) 初始化单元设计变量，即给定设计域内的每个单元一个初始单元相对密度;
% 4) 计算各离散单元的材料特性参数，计算单元刚度矩阵，组装结构总刚度矩阵，计算结点位移;
% 5) 计算总体结构的柔度值及其敏度值，求解拉格朗日乘子;
% 6) 用优化准则方法进行设计变量更新;
% 7) 检查结果的收敛性，如未收敛则转4)，循环迭代，如收敛则转8);
% 8) 输出目标函数值及设计变量值，结束计算。
%===================================================================================

% x是设计变量，给设计域内的每个单元一个初始相对密度，值为volfrac
x(1:nely,1:nelx)=volfrac;

% loop储存迭代次数
loop=0;

% change储存每次迭代之后目标函数的改变值，用以判断是否收敛
change=1.;

% 当目标函数改变量<=0.01时说明收敛，结束迭代
while change > 0.01
loop=loop + 1;

% 将前一次的设计变量赋值给xold，x用来储存这一次的结果，之后还要比较它们以判断是否收敛
xold=x;

% 每次迭代都进行一次有限元分析，计算结点位移，并储存在全局位移数组U中
[U]=FE(nelx,nely,x,penal);

% 计算单元刚度矩阵
[KE]=lk;

% c是用来储存目标函数的变量
c=0.;

% 遍历设计域矩阵元素，从左上角第一个元素开始，一列一列
for ely=1:nely
for elx=1:nelx

% 节点位移储存在U中，如果想获得节点位移必须先知道节点编号，进行索引;
% 节点编号可以根据当前单元在设计域中的位置算出;
% n1是左上角节点编号，n2是右上角节点编号;
n1=(nely+1)*(elx-1)+ely;
n2=(nely+1)*elx+ely;

% 局部位移数组Ue储存4个节点共8个自由度位移，每个节点分别有x、y两个自由度;
% 因为是矩形单元，所以根据n1、n2两个节点的编号可以推演出单元所有节点的自由度编号;
% 顺序是：[左上x；左上y；右上x；右上y；右下x；右下y；左下x；左下y];
% 只适用于矩形单元划分网格;
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);

% SIMP模型，将设计变量x从离散型变成指函数型，指数就是惩罚因子：x(ely,elx)^penal
% 计算总体结构的柔度值，这里目标函数是柔度最小，参见论文中公式（1）
c=c + x(ely,elx)^penal*Ue'*KE*Ue;

% 计算总体结构的敏度值，实际上dc就是c对Xe的梯度，参见论文中公式（4）
dc(ely,elx)=-penal*x(ely,elx)^(penal-1)*Ue'*KE*Ue;
end
end

% 无关网格敏度过滤
[dc]=check(nelx,nely,rmin,x,dc);

% 采用优化准则法（OC）求解当前模型，得出满足体积约束的结果，更新设计变量
[x]=OC(nelx,nely,x,volfrac,dc);

% 更新目标函数改变值
change=max(max(abs(x-xold)));

% 打印迭代信息: It.迭代次数，Obj.目标函数，Vol.材料体积比，ch.迭代改变量
disp([' It.: ' sprintf('%4i',loop) ' Obj.: ' sprintf('%10.4f',c) ...
' Vol.: ' sprintf('%6.3f',sum(sum(x))/(nelx*nely)) ...
' ch.: ' sprintf('%6.3f',change )])

% 当前迭代结果图形显示
colormap(gray); imagesc(-x); axis equal; axis tight; axis off;pause(1e-6);
end


% 采用优化准则法（OC）迭代子程序 -------------------------------------------------------
% 数学模型主要求解算法有：优化准则法(OC)、序列线性规划法(SLP)、移动渐进线法(MMA);
% OC适用于单约束优化问题求解，比如这里的"体积约束下的柔度最小化"问题，当求解复杂的多约束拓扑
% 优化问题时，采用SLP和MMA通常更方便;
% 参见论文中公式（2）（3）
function[xnew]=OC(nelx,nely,x,volfrac,dc)
% Input： 水平单元数nelx, 竖直单元数nely, 设计变量x, 材料体积比volfrac, 目标函数灵敏度dc;
% Output： 更新后的设计变量xnew;

% 定义一个取值区间，二分法，得到满足体积约束的拉格朗日算子
l1=0; l2=100000;

% 正向最大位移
move=0.2;

while (l2-l1 > 1e-4)

% 二分法，取区间中点
lmid=0.5*(l2+l1);

% 参见论文中的公式（2）
% sqrt(-dchttps://www.zhihu.com/question/lmid)对应公式中Be^eta（eta=1/2），eta阻尼系数是为了确保计算的收敛性
xnew=max(0.001,max(x-move,min(1.,min(x+move,x.*sqrt(-dchttps://www.zhihu.com/question/lmid)))));

% sum(sum(xnew))是更新后的材料体积, volfrac*nelx*nely是优化目标，用它们的差值判断是否收敛
% sum()可以去看看help，注意一下行、列问题，不看也行，反正知道sum(sum())是所有元素求和就行
if sum(sum(xnew)) - volfrac*nelx*nely > 0
l1=lmid;
else
l2=lmid;
end
end


% 无关网格敏度过滤子程序 --------------------------------------------------------------
% 参见论文中公式（5）（6）
function[dcn]=check(nelx,nely,rmin,x,dc)
% Input： 水平单元数nelx, 竖直单元数nely, 敏度过滤半径rmin, 设计变量x, 总体结构敏度dc;
% Output： 过滤后的目标函数敏度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)
for l=max(j-floor(rmin),1):min(j+floor(rmin),nely)

% 参见论文中公式（6）,fac即公式中卷积算子Hf
% qrt((i-k)^2+(j-l)^2) 是计算此单元与相邻单元的距离，即公式中dist(e,f)
fac=rmin-sqrt((i-k)^2+(j-l)^2);
sum=sum+max(0,fac);

% 参见论文中公式（5）
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


% 有限元求解子程序 --------------------------------------------------------------------
function[U]=FE(nelx,nely,x,penal)
% Input： 水平单元数nelx, 竖直单元数nely, 设计变量x, 惩罚因子penal;
% Output： 全局节点位移U;

% 计算单元刚度矩阵
[KE]=lk;

% 总体刚度矩阵的稀疏矩阵
K=sparse(2*(nelx+1)*(nely+1), 2*(nelx+1)*(nely+1));

% 力矩阵的稀疏矩阵
F=sparse(2*(nely+1)*(nelx+1),1);

% U清零，用来保存更新的全局节点位移
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;

% 同上主程序，计算单元4个节点8个自由度
edof=[2*n1-1; 2*n1; 2*n2-1; 2*n2; 2*n2+1; 2*n2+2; 2*n1+1; 2*n1+2];

% 将单元刚度矩阵KE 组装成 总体刚度矩阵K
K(edof,edof)=K(edof,edof) + x(ely,elx)^penal*KE;
end
end

% 施加载荷，本算例应用了一个在左上角的垂直单元力!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
F(2,1)=-1;%每个节点有两个自由度 (DOF)，一个在 x 方向，另一个在 y 方向。 自由度被标记为
%x方向：(2·nn)-1
%y方向：(2·nn)
%其中 nn 是节点号。
%由于负载施加在节点 1 的 y 方向上，对应的自由度为 2*1=2。
%因此，回答您的问题，代码中的第 79 行：F(2,1)=-1;，是一个力 的 -1 应用于节点 1 的 y 方向。


% 施加约束，消除线性方程中的固定自由度来实现支承结构，本算例左边第一列和右下角固定!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
fixeddofs=union([1:2:2*(nely+1)],[2*(nelx+1)*(nely+1)]);

% 剩下的不加约束的节点自由度，setdiff()从..中除去..
alldofs=[1:2*(nely+1)*(nelx+1)];
freedofs=setdiff(alldofs,fixeddofs);

% 求解线性方程组，得到各节点自由度的位移值储存在U中
U(freedofs,:)=K(freedofs,freedofs) \\ F(freedofs,:);

% 受约束节点固定自由度位移值为0
U(fixeddofs,:)=0;


% 求解单元刚度矩阵子程序 --------------------------------------------------------------
% 有限元方法计算的一个重要的系数矩阵，表征单元体的受力与变形关系;
% 特点：对称性、奇异性、主对角元素恒正、奇数（偶数）行和为0;
% 矩形单元4节点 8*8矩阵;
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)];