MATLAB在处理平面有限元问题和梁弯曲问题上有很强的能力,主要体现在以下几个方面:
建模与网格划分
MATLAB内置了方便的图形界面工具(pdetoolbox等),可以快速对几何模型进行二维三维网格划分,生成有限元分析需要的网格。
求解器
MATLAB内置了多种求解偏微分方程的有限元求解器,如:适用于线性和非线性问题的有限元求解器assembler,适用于时间依赖问题的时间自适应求解器timestepper等。这些求解器可以自动组装刚度矩阵并求解。
后处理
MATLAB支持多种后处理,可以方便的进行解的可视化,绘制云图、矢量图、剪切面;也可以计算量化解的误差与收敛性,绘制相关曲线。
应用程序接口
MATLAB提供了应用程序接口,允许自定义复杂几何模型,加载网格,实现自定义的有限元运算与后处理。
梁弯曲分析
对于梁的弯曲问题,MATLAB提供了专门的梁挠度计算与绘制工具beamCrippling,可以快速获得梁的弯矩与切力图,并自动检查其强度。
总的来说,利用MATLAB强大的工具箱与API,可以高效的对平面及梁挠问题进行预处理、求解与后处理,得到定量结果及直观的可视化。
主程序:
%% ---------------四边形八节点等参元 matlab计算程序----------------------------
clear all;clc; close all;
format short e ;
%%读入控制数据
E=1E5; %弹性模量
v=0.25; % 泊松比
h=1; %厚度
NumberElement=4; %单元数
NumberNode=23; % 总结点数
ElementNode=8; %单元节点数
ForcePoint=1; %受力结点数
NumberConstraint=3; %约束结点个数
%% 节点坐标 x y
% 结点号 x,y坐标(整体坐标下)
gNdt = [0.0000000e+000 -0.5000000e+000;
1.0000000e+000 -0.5000000e+000;
2.0000000e+000 -0.5000000e+000;
3.0000000e+000 -0.5000000e+000;
4.0000000e+000 -0.5000000e+000;
5.0000000e+000 -0.5000000e+000;
6.0000000e+000 -0.5000000e+000;
7.0000000e+000 -0.5000000e+000;
8.0000000e+000 -0.5000000e+000;
0.0000000e+000 0.0000000e+000;
2.0000000e+000 0.0000000e+000;
4.0000000e+000 0.0000000e+000;
6.0000000e+000 0.0000000e+000;
8.0000000e+000 0.0000000e+000;
0.0000000e+000 0.5000000e+000;
1.0000000e+000 0.5000000e+000;
2.0000000e+000 0.5000000e+000;
3.0000000e+000 0.5000000e+000;
4.0000000e+000 0.5000000e+000;
5.0000000e+000 0.5000000e+000;
6.0000000e+000 0.5000000e+000;
7.0000000e+000 0.5000000e+000;
8.0000000e+000 0.5000000e+000];
%% 单元节点
gElt = [1 2 3 11 17 16 15 10;
3 4 5 12 19 18 17 11;
5 6 7 13 21 20 19 12;
7 8 9 14 23 22 21 13]; % 单元定义： 单元结点号（逆时针）
FPOIN=[9 0 -1];% 节点力：结点号、X方向力(向右正)，Y方向力(向上正)
FIXED=[1 1 1;
10 1 1;
15 1 1];
%约束信息数组(n,3) n:受约束节点数目, (n,1):约束点号
%(n,2)与(n,3)分别为约束点x方向和y方向的约束情况,受约束为1否则为0
%========平面应力问题的求解==============
% 刚度矩阵的生成
%计算刚度矩阵，并对约束条件进行处理
Ke=zeros(2*ElementNode,2*ElementNode); % 单元刚度矩阵并清零
HK=zeros(2*NumberNode,2*NumberNode); % 总刚矩阵并清零
%调用子程序 生成单元刚度矩阵
for m=1:NumberElement %m为单元号
Ke=K(E,v,h,...
gNdt(gElt(m,1),1),gNdt(gElt(m,1),2),...
gNdt(gElt(m,3),1),gNdt(gElt(m,3),2),...
gNdt(gElt(m,5),1),gNdt(gElt(m,5),2),...
gNdt(gElt(m,7),1),gNdt(gElt(m,7),2)); %调用单元刚度矩阵
a=gElt(m,:); %临时向量,用来记录当前单元的节点编号
% 对总刚度矩阵的处理
for j=1:8
for k=1:8
HK((a(j)*2-1):a(j)*2,(a(k)*2-1):a(k)*2)=HK((a(j)*2-1):a(j)*2,(a(k)*2-1):a(k)*2)+...
Ke(j*2-1:j*2,k*2-1:k*2);
end
end
end
%—————————————————————————————————
% 对荷载向量进行处理
FORCE=zeros(2*NumberNode,1); % 张成总荷载向量并清零
for i=1:ForcePoint
b1=FPOIN(i,1)*2-1;b2=FPOIN(i,1)*2; %FPION(i,1)为作用点
FORCE(b1)=FPOIN(i,2); %FPION(i,2)为x方向的节点力
FORCE(b2)=FPOIN(i,3); %FPION(i,3)为y方向的节点力
end
%—————————————————————————————————
%将约束信息加入总刚,总荷载
for i=1:NumberConstraint
if FIXED(i,2)==1
c1=2*FIXED(i,1)-1;
HK(c1,:)=0; %将一约束序号处的总刚列向量清0
HK(:,c1)=0; %将一约束序号处的总刚行向量清0
HK(c1,c1)=1; %将行列交叉处的元素置为1
FORCE(c1)=0;
end
if FIXED(i,3)==1
c2=2*FIXED(i,1);
HK(c2,:)=0;
HK(:,c2)=0;
HK(c2,c2)=1;
FORCE(c2)=0;
end
end
%—————————————————————————————————
HK
Displacement=HK\FORCE;%计算节点位移向量
%% 转换
Displacement2=reshape(Displacement,2,length(Displacement)/2)';
gNdt2=gNdt+Displacement2*100;
figure;
gNdt;
for i=1:size(gElt,1)
for j=1:8
%画点
plot(gNdt(gElt(i,:),1),gNdt(gElt(i,:),2),'b-');
hold on;
%画线
plot(gNdt2(gElt(i,:),1),gNdt2(gElt(i,:),2),'r-');
hold on;
end
end
axis equal;
xlabel('x(m)');
ylabel('y(m)');
title('八节点四边形等参单元:梁位移(放大100倍)');
%———————————求解单元应力————————————————
stress=zeros(3,NumberElement);
for m=1:NumberElement
u(1:16)=0;
d=gElt(m,:); %临时向量,用来记录当前单元的节点编号
for i=1:ElementNode
u(i*2-1:i*2)=Displacement(d(i)*2-1:d(i)*2);
%从总位移向量中取出当前单元的节点位移
end
D=(E/(1-v*v))*[1 v 0;v 1 0;0 0 (1-v)/2];%弹性矩阵
%形成应变矩阵BM
BM=zeros(3,16);
for i=1:ElementNode
J=Jacobi(gNdt(gElt(m,1),1),gNdt(gElt(m,1),2),...
gNdt(gElt(m,3),1),gNdt(gElt(m,3),2),...
gNdt(gElt(m,5),1),gNdt(gElt(m,5),2),...
gNdt(gElt(m,7),1),gNdt(gElt(m,7),2),0,0);
[N_s,N_t]=DHS(0,0);
B1i=J(2,2)*N_s(i)-J(1,2)*N_t(i);
B2i=-J(2,1)*N_s(i)+J(1,1)*N_t(i);
BM(1:3,2*i-1:2*i)=[B1i 0;0 B2i;B2i B1i]/det(J);
end
stressm=D*BM*u';
stress(:,m)=stressm;
end
stress %输出应力
程序结果:
