Matlab与有限元
速来！下面是矩形薄板基于三节点三角形单元的有限元分析（matlab程序），请问怎么将其修改成受均布荷载的简支梁（相同的分析模型）
下图为简支梁的半结构模型：
原代码如下：
clc
tic;
Initial_info = [0.09 0.06 18 12];
disp(['该程序计算的是',num2str(Initial_info(3)+1),'×',num2str(Initial_info(4)+1),'=',...
num2str((Initial_info(3)+1)*(Initial_info(4)+1)),'个结点的有限元模型']);
LX = Initial_info(1);
LY = Initial_info(2);
nx = Initial_info(3);
ny = Initial_info(4);
ne = 2 * nx * ny;
np = (nx+1)*(ny+1);
for i = 1:nx+1; j = 1:ny+1;
Np(i,j) = j+(i-1)*(ny+1);
end
% 生成节点编号矩阵Np
for i = 1: nx+1; j = 1:ny+1;
XX(i,j) = (i-1)*LX/nx;
YY(i,j) = (j-1)*LY/ny;
end
XY = [reshape(XX',np,1),reshape(YY',np,1)];
nx2 = nx/2;
Np1 = Np(1:nx2+1,:);
Np2 = Np(nx2+1:end,:);
for i = 1: nx2*ny
if rem(i,nx2) == 0
xp = nx2;
yp = i/nx2;
else
xp = rem(i,nx2);
yp = fix(i/nx2)+1;
end
Dof1(i,:) = [Np1(xp,yp),Np(xp+1,yp),Np(xp,yp+1)];
Dof1(i+nx2 * ny,:) = [Np(xp+1,yp),Np(xp+1,yp+1),Np(xp,yp+1)];
Dof2(i,:) = [Np2(xp,yp),Np(xp+1,yp),Np(xp+1,yp+1)];
Dof2(i+nx2 * ny,:) = [Np2(xp,yp),Np(xp+1,yp+1),Np2(xp,yp+1)];
end
Dof = [Dof1;Dof2];
for i = 1:ne
unit(i,:) = [XY(Dof(i,1),1),XY(Dof(i,2),1),XY(Dof(i,3),1),...
XY(Dof(i,1),2),XY(Dof(i,2),2),XY(Dof(i,3),2)];
end
disp('前处理完成');
% 前处理完成
%%% 单元刚度矩阵
E = 2*10^11;
u = 0; % 平面应力问题
D = E/(1-u^2)*[1 u 0; u 1 0;0 0 (1-u)/2];
for i = 1:ne
xi = unit(i,1);
yi = unit(i,4);
xj = unit(i,2);
yj = unit(i,5);
xm = unit(i,3);
ym = unit(i,6);
ai = xj * ym - xm * yj;
aj = xm * yi - xi * ym;
am = xi * yj - xj * yi;
bi = yj - ym;
bj = ym - yi;
bm = yi - yj;
ci = -xj + xm;
cj = -xm + xi;
cm = -xi + xj;
A = abs((ai+aj+am)/2);
B = [bi 0 bj 0 bm 0
0 ci 0 cj 0 cm
ci bi cj bj cm bm];
Be{i,1} = B/2/A;
ke{i,1} = [Be{i,1}]'*D*Be{i,1}*A;
end
%%% 总刚度矩阵叠加
K = sparse(2*np,2*np);
for ie = 1:ne
a = Dof(ie,1);
b = Dof(ie,2);
c = Dof(ie,3);
DOF(1) = 2*a-1;
DOF(2) = 2*a;
DOF(3) = 2*b-1;
DOF(4) = 2*b;
DOF(5) = 2*c-1;
DOF(6) = 2*c;
for n1 = 1:6
for n2 = 1:6
K(DOF(n1),DOF(n2)) = K(DOF(n1),DOF(n2))+ke{ie,1}(n1,n2);
end
end
end
%%% 单元等效结点荷载
y = (0:LY);
P = (10^7/0.03)*y - 10^7;
Re = sparse(ne,6);
for i = 1:ne
switch i
case num2cell(1:ne/2-nx2)
Pe = [0 0 0 0 0 0];
case num2cell(ne/2-nx2+1:ne/2)
Pe = -LX*P*[0,0,0,1,0,1]/nx/2;
case num2cell(ne/2+1:ne-nx2)
Pe = [0 0 0 0 0 0 ];
otherwise
Pe = -LX*P*[0,0,0,1,0,1]/nx/2;
end
Re(i,:) = Pe;
end
%%% 荷载叠加
Rr = sparse(1,2*np);
for i = 1:ne
a = Dof(i,1);
b = Dof(i,2);
c = Dof(i,3);
DOF(1) = 2*a-1;
DOF(2) = 2*a;
DOF(3) = 2*b-1;
DOF(4) = 2*b;
DOF(5) = 2*c-1;
DOF(6) = 2*c;
for ni = 1:6
Rr(DOF(n1)) = Rr(DOF(n1)) + Re(i,n1);
end
end
%%% 生成需处理的序列
cp = [1:nx+1];
ctype = ones(1,length(cp));
ctype(nx2+1) = 2;
cp_all = (cp-1)*(ny-1)+1;
p_stake = zeros(1,2*length(cp));
for i = 1:length(cp)
switch ctype(i)
case{2}
p_stake(2*i) = 2*cp_all(i);
p_stake(2*i-1) = 2*cp_all(i)-1;
case{1}
p_stake(2*i) = 2*cp_all(i);
p_stake(2*i-1) = [];
otherwise
p_stake(2*i) = [];
p_stake(2*i-1) = 2*cp_all(i)-1;
end
end
[m,j] = find(p_stake == 0);
p_stake(:,j) = [];
%%% 处理对应的行列
K_d = K;
K_f = K;
K_d(p_stake,:) = [];
K_d(:,p_stake) = [];
K_f(:,p_stake) = [];
K_f = K_f(p_stake,:);
Rr_unkown = Rr;
Rr_unkown(:,p_stake) = [];
RR = transpose(Rr_unkown);
[L,U] = lu(K_d);
UU = U\(L\RR);
Rx = K_f*UU;
%%% 数值计算部分
UU_all = UU';
for i = 1:length(p_stake)
UU_all = [UU_all(:,1:p_stake(i)-1),0 ,UU_all(:,p_stake(i):end)];
end
for i = 1:np
UU_info(i,:) = [UU_all(2*i-1),UU_all(2*i)];
end
%%% 出图部分
%%%梁尺寸及荷载图
figure;
set(gcf,'outerposition',get(0,'ScreenSize'));
set(gcf,'name','梁的尺寸及荷载');
line([-0.015,-0.01],[0.06,0.06]),hold on
line([-0.01,-0.01],[0.06,0]),hold on
line([-0.01,-0.005],[0,0]),hold on
line([-0.005,-0.015],[0,0.06]),hold on
line([0.105,0.1],[0.06,0.06]),hold on
line([0.1,0.1],[0.06,0]),hold on
line([0.1,0.095],[0,0]),hold on
line([0.095,0.105],[0,0.06]),hold on
rectangle('position',[0,0,0.09,0.06]),hold on
quiver(-0.01,0.06,-0.0055,0,'LineWidth',2.0,'MaxheadSize',0.8,'color','k'),hold on
quiver(-0.01,0.05,-0.0038,0,'LineWidth',2.0,'MaxheadSize',0.8,'color','k'),hold on
quiver(-0.01,0.04,-0.002,0,'LineWidth',2.0,'MaxheadSize',0.8,'color','k'),hold on
quiver(-0.01,0.02,-0.002,0,'LineWidth',2.0,'MaxheadSize',0.8,'color','k'),hold on
quiver(-0.01,0.01,-0.0038,0,'LineWidth',2.0,'MaxheadSize',0.8,'color','k'),hold on
quiver(-0.01,0,-0.0055,0,'LineWidth',2.0,'MaxheadSize',0.8,'color','k'),hold on
quiver(0.1,0.06,-0.0055,0,'LineWidth',2.0,'MaxheadSize',0.8,'color','k'),hold on
quiver(0.1,0.05,-0.0038,0,'LineWidth',2.0,'MaxheadSize',0.8,'color','k'),hold on
quiver(0.1,0.04,-0.002,0,'LineWidth',2.0,'MaxheadSize',0.8,'color','k'),hold on
quiver(0.1,0.02,-0.002,0,'LineWidth',2.0,'MaxheadSize',0.8,'color','k'),hold on
quiver(0.1,0.01,-0.0038,0,'LineWidth',2.0,'MaxheadSize',0.8,'color','k'),hold on
quiver(0.1,0,-0.0055,0,'LineWidth',2.0,'MaxheadSize',0.8,'color','k'),hold on
text(-0.0155,0.062,'1000N/cm^2')
text(-0.0105,-0.002,'1000N/cm^2')
text(0.1,0.062,'1000N/cm^2')
text(0.095,-0.002,'1000N/cm^2')
text(0,-0.002,'梁宽9cm，高6cm，厚1cm。E=2×10^11N/m^2，μ=0（此图单位：m）')
axis([-0.025 0.12 -0.015 0.07]),hold on;
%网格图
figure;
set(gcf,'outerposition',get(0,'ScreenSize'));
set(gcf,'name','网格划分');
for i = 1:ne
line([unit(i,1:3),unit(i,1)],[unit(i,4:6),unit(i,4)]);
end
axis([-0.025 0.12 -0.015 0.07]),hold on
xytext = num2str([1:np]');
text(XY(:,1)+LX/nx/8,XY(:,2)+LY/ny/4,xytext);
print(gcf,'-dbitmap','model.bmp');
%位移矢量图
figure;
set(gcf,'outerposition',get(0,'ScreenSize'));
set(gcf,'name','结点的位移矢量场');
quiver(XY(:,1),XY(:,2),-UU_info(:,1),UU_info(:,2)),
axis([-0.025 0.12 -0.015 0.07]),hold on;
for i = 1:ne
plot([unit(i,1:3),unit(i,1)],[unit(i,4:6),unit(i,4)],'r:');
end
disp('梁顶面结点');
JD = 13:13:247
disp('梁顶面结点位移');
WY = [-UU_info(13);-UU_info(26);-UU_info(39);-UU_info(52);-UU_info(65);-UU_info(78);-UU_info(91);...
-UU_info(104);-UU_info(117);-UU_info(130);-UU_info(143);-UU_info(156);-UU_info(169);-UU_info(182);...
-UU_info(195);-UU_info(208);-UU_info(221);-UU_info(234);-UU_info(247);]
%结点位移
figure;
set(gcf,'outerposition',get(0,'ScreenSize'));
set(gcf,'name','结点位移');
plot(JD,WY,'r:','LineWidth',2.0),grid on
toc;
- 你可以看下这个问题的回答https://ask.csdn.net/questions/7546294