基于某MATLAB的平面刚架静力分析报告.docx
- 文档编号:17788967
- 上传时间:2023-08-03
- 格式:DOCX
- 页数:20
- 大小:123.27KB
基于某MATLAB的平面刚架静力分析报告.docx
《基于某MATLAB的平面刚架静力分析报告.docx》由会员分享,可在线阅读,更多相关《基于某MATLAB的平面刚架静力分析报告.docx(20页珍藏版)》请在冰点文库上搜索。
基于某MATLAB的平面刚架静力分析报告
基于MATLAB的平面刚架静力分析
为了进一步理解有限元方法计算的过程,本文根据矩阵位移法的基本原理应用MATLAB编制计算程序对以平面刚架结构进行了静力分析。
本文还利用ANSYS大型商用有限元分析软件对矩阵位移法的计算结果进行校核,发现两者计算结果相当吻合,验证了计算结果的可靠性。
一、问题描述
如图1所示的平面刚架,各杆件的材料及截面均相同,E=210GPa,截面为0.12×0.2m的实心矩形,现要求解荷载作用下刚架的位移和内力。
图1
二、矩阵位移法计算程序编制
为编制程序方便考虑,本文计算中采用“先处理法”。
具体的计算步骤如下。
(1)对结构进行离散化,对结点和单元进行编号,建立结构(整体)坐标系和单元(局部)坐标系,并对结点位移进行编号;
(2)对结点位移分量进行编码,形成单元定位向量
;
(3)建立按结构整体编码顺序排列的结点位移列向量
,计算固端力
、等效结点荷载
及综合结点荷载列向量
;
(4)计算个单元局部坐标系的刚度矩阵,通过坐标变换矩阵
形成整体坐标系下的单元刚度矩阵
;
(5)利用单元定位向量形成结构刚度矩阵
;
(6)按式
求解未知结点位移;
(7)计算各单元的杆端力
。
根据上述步骤编制了平面刚架的分析程序。
程序中单元刚度矩阵按下式计算。
转换矩阵则按下式计算。
计算程序框图如图2所示,具体的程序代码见附录1。
图2MATLAB矩阵分析法程序框图
三、解题步骤
取整体坐标系如图3所示,对结构进行离散化,对结点和单元进行编号如图4所示,局部坐标系用单元中箭头的方向表示,原始数据如下:
图3图4
刚架结点输入矩阵为,
x=[00;0-5;1.63-6.37;4-5;4-1;42];
各单元定位向量为,
locvec1=[123000];
locvec2=[123456];
locvec3=[456789];
locvec4=[789101112];
locvec5=[101112000];
输入截面参数,
E=2.1e11;%E=210GPa
a=0.12;%矩形截面长0.12m
b=0.2;%矩形截面宽0.2m
输入整体坐标系下各单元结点荷载列阵,
f(1,:
)=zeros(1,6);
f(2,:
)=[000040e30];
f(3,:
)=zeros(1,6);
f(4,:
)=[000-50e300];
f(5,:
)=zeros(1,6);
输入整体坐标系下单元1等效节点荷载
q=10e3;%10kN/m
fe=[0.5*q*l
(1),0,-q*l
(1)^2/12,0.5*q*l
(1),0,q*l
(1)^2/12];
由此计算得到平面刚架整体坐标系下的结点位移(m),
d=
0.0035
0.0000
-0.0004
0.0030
-0.0005
-0.0004
0.0027
0.0000
0.0016
-0.0051
0.0000
-0.0006
各个单元的杆端力如表1所示,
表1各单元杆端力
单元
1
2
3
4
5
i端
Fx(kN)
-17917.047
17917.05
17917.05
17917.05
-32083
Fy(kN)
17507.3731
-17507.4
22492.63
22492.63
22492.63
Mz(kN·m)
1897.83076
-1897.83
2092.833
-26668.3
44999.85
j端
Fx(kN)
-32082.953
-17917
-17917
32082.95
32082.95
Fy(kN)
-17507.373
-22492.6
-22492.6
-22492.6
-22492.6
Mz(kN·m)
-37312.596
-2092.83
26668.34
-44999.8
51249.01
四、计算结果校核
在ANSYS中使用beam3单元,按照如图4所示的离散结构建立平面刚架模型施加约束和荷载,得到的有限元模型如图5所示。
计算分析后得到结构的轴力图、剪力图和弯矩图如图6、7、8所示,命令流见附录2。
图5有限元模型
图6轴力图(kN)
图7剪力图(kN)
图8弯矩图(kN·m)
从ANSYS计算结果中提取各结点位移、内力,并与矩阵位移法分析的结果比较,得到表2、3。
表2各节点位移比较
节点号
项目
矩阵位移法
ANSYS
误差
1
Ux(m)
0
0
0
Uy(m)
0
0
0
Rotz(rad)
0
0
0
2
Ux(m)
0.003478
0.00348
-2E-06
Uy(m)
0.0000174
0.0000174
0
Rotz(rad)
-0.00037
-0.000365
-5E-06
3
Ux(m)
0.003018
0.00302
-2E-06
Uy(m)
-0.00051
-0.000512
2E-06
Rotz(rad)
-0.00038
-0.000378
-2E-06
4
Ux(m)
0.002687
0.00269
-3E-06
Uy(m)
0.0000312
0.0000312
0
Rotz(rad)
0.001624
0.00162
4E-06
5
Ux(m)
-0.00513
-0.005145
1.5E-05
Uy(m)
0.0000134
0.0000134
0
Rotz(rad)
-0.00056
-0.000557
-3E-06
6
Ux(m)
0
0
0
Uy(m)
0
0
0
Rotz(rad)
0
0
0
表3各结点内力比较
节点号
项目
矩阵位移法
ANSYS
误差
1
Fx(kN)
-32.083
-32.080
-0.003
Fy(kN)
-17.507
-17.503
-0.004
Mz(kN·m)
-37.313
-37.307
-0.006
2
Fx(kN)
-17.917
-17.920
0.003
Fy(kN)
17.507
17.503
0.004
Mz(kN·m)
1.898
1.909
-0.011
3
Fx(kN)
-17.917
-17.920
0.003
Fy(kN)
-22.493
-22.497
0.004
Mz(kN·m)
-2.093
-2.110
0.018
4
Fx(kN)
-17.917
-17.920
0.003
Fy(kN)
-22.493
-22.497
0.004
Mz(kN·m)
26.668
26.682
-0.014
5
Fx(kN)
-32.083
-32.080
-0.003
Fy(kN)
-22.493
-22.497
0.004
Mz(kN·m)
-45.000
-44.999
-0.001
6
Fx(kN)
32.083
32.080
0.003
Fy(kN)
-22.493
-22.497
0.004
Mz(kN·m)
51.249
51.239
0.010
由表2、表3的结果对比可知,两种方法的计算结果十分接近,误差均可以忽略不计,从而验证了计算结果的可靠性与准确性。
四、结论
通过对一个平面刚架静力分析问题的求解,本文编制的平面刚架静力分析程序计算结果与有限元软件ANSYS计算结果校核后,发现两者计算结果十分接近,误差可忽略不计,说明该程序计算结果的可靠性与精确性。
附录1矩阵位移法计算程序
pmgj.m计算主程序
clear
clc
%结点坐标
x=[00;0-5;1.63-6.37;4-5;4-1;42];
%x1=0;y1=0;
%x2=0;y2=-5;
%x3=1.63;y3=-6.37;
%x4=4;y4=-5;
%x5=4;y5=-1;
%x6=4;y6=2;
%单元定位向量
locvec1=[123000];
locvec2=[123456];
locvec3=[456789];
locvec4=[789101112];
locvec5=[101112000];
%刚架的材料特性截面特性
E=2.1e11;%E=210GPa
a=0.12;%矩形截面长0.12m
b=0.2;%矩形截面宽0.2m
A=a*b;
I=(a*b^3)/12;
clearab
%单元长度计算
fori=1:
5
dx=x(i,1)-x(i+1,1);
dy=x(i,2)-x(i+1,2);
l(i)=(dx^2+dy^2)^0.5;
end
%生成转换矩阵
t1=coortrans(x(2,1),x(2,2),x(1,1),x(1,2));
t2=coortrans(x(2,1),x(2,2),x(3,1),x(3,2));
t3=coortrans(x(3,1),x(3,2),x(4,1),x(4,2));
t4=coortrans(x(4,1),x(4,2),x(5,1),x(5,2));
t5=coortrans(x(5,1),x(5,2),x(6,1),x(6,2));
%结点荷载(结构坐标系下)
f(1,:
)=zeros(1,6);
f(2,:
)=[000040e30];
f(3,:
)=zeros(1,6);
f(4,:
)=[000-50e300];
f(5,:
)=zeros(1,6);
%单元等效节点荷载(结构坐标系下)
q=10e3;%10kN/m
fe=[0.5*q*l
(1),0,-q*l
(1)^2/12,0.5*q*l
(1),0,q*l
(1)^2/12];
%单元坐标系下的值
fe1=[0,0.5*q*l
(1),q*l
(1)^2/12,0,0.5*q*l
(1),-q*l
(1)^2/12];
%生成单元刚度矩阵(单元坐标系)
k1=elemstiffm(E,A,l
(1),I);
k2=elemstiffm(E,A,l
(2),I);
k3=elemstiffm(E,A,l(3),I);
k4=elemstiffm(E,A,l(4),I);
k5=elemstiffm(E,A,l(5),I);
%将单元坐标系下的单元刚度矩阵转换到结构坐标系下
K1=t1'*k1*t1;
K2=t2'*k2*t2;
K3=t3'*k3*t3;
K4=t4'*k4*t4;
K5=t5'*k5*t5;
%组装总体刚度矩阵
K=zeros(12);
F=zeros(12,1);
NonLog=(locvec1~=0);
ij=find(NonLog);
IJ=locvec1(NonLog);
K(IJ,IJ)=K(IJ,IJ)+K1(ij,ij);
F(IJ)=F(IJ)+f(1,ij)'+fe(ij)';
clearNonLogijIJ
NonLog=(locvec2~=0);
ij=find(NonLog);
IJ=locvec2(NonLog);
K(IJ,IJ)=K(IJ,IJ)+K2(ij,ij);
F(IJ)=F(IJ)+f(2,ij)';
clearNonLogijIJ
NonLog=(locvec3~=0);
ij=find(NonLog);
IJ=locvec3(NonLog);
K(IJ,IJ)=K(IJ,IJ)+K3(ij,ij);
F(IJ)=F(IJ)+f(3,ij)';
clearNonLogijIJ
NonLog=(locvec4~=0);
ij=find(NonLog);
IJ=locvec4(NonLog);
K(IJ,IJ)=K(IJ,IJ)+K4(ij,ij);
F(IJ)=F(IJ)+f(4,ij)';
clearNonLogijIJ
NonLog=(locvec5~=0);
ij=find(NonLog);
IJ=locvec5(NonLog);
K(IJ,IJ)=K(IJ,IJ)+K5(ij,ij);
F(IJ)=F(IJ)+f(5,ij)';
clearNonLogijIJ
%节点位移
d=K\F;
%单元1杆端力(结构坐标系下)
d1=zeros(6,1);
NonLog=(locvec1~=0);
ij=find(NonLog);
ij1=find(~NonLog);
IJ=locvec1(NonLog);
d1=d(IJ);
d1(ij1)=0;
F1=K1*d1-fe';
%单元2杆端力(结构坐标系下)
d2=zeros(6,1);
NonLog=(locvec2~=0);
ij=find(NonLog);
ij1=find(~NonLog);
IJ=locvec2(NonLog);
d2=d(IJ);
d2(ij1)=0;
F2=K2*d2-f(2,:
)';
%单元3杆端力(结构坐标系下)
d3=zeros(6,1);
NonLog=(locvec3~=0);
ij=find(NonLog);
ij1=find(~NonLog);
IJ=locvec3(NonLog);
d3=d(IJ);
d3(ij1)=0;
F3=K3*d3-f(3,:
)';
%单元4杆端力(结构坐标系下)
d4=zeros(6,1);
NonLog=(locvec4~=0);
ij=find(NonLog);
ij1=find(~NonLog);
IJ=locvec4(NonLog);
d4=d(IJ);
d4(ij1)=0;
F4=K4*d4-f(4,:
)';
%单元5杆端力(结构坐标系下)
d5=zeros(6,1);
NonLog=(locvec5~=0);
ij=find(NonLog);
ij1=find(~NonLog);
IJ=locvec5(NonLog);
d5=d(IJ);
d5(ij1)=0;
F5=K5*d5-f(5,:
)';
coortrans.m转换矩阵生成函数
%转换矩阵生成函数(单元坐标系->结构坐标系)
%y=coortrans(x1,y1,x2,y2),x1,y1,x2,y2为单元两端结点在结构坐标系下的坐标,单位都为m
functiony=coortrans(x1,y1,x2,y2)
a=atan((y2-y1)/(x2-x1));
c=cos(a);
s=sin(a);
t=zeros(6);
t(1,1)=c;t(1,2)=s;
t(2,1)=-s;t(2,2)=c;
t(3,3)=1;
t(4,4)=c;t(4,5)=s;
t(5,4)=-s;t(5,5)=c;
t(6,6)=1;
y=t;
end
elemstiffm.m单元刚度矩阵生成函数
%单元刚度矩阵生成函数
%y=elemstiffm(E,A,l,I),E,A,l,I均采用国际单位制Pam2mm4
functiony=elemstiffm(E,A,l,I)
i1=E*A/l;
i2=12*E*I/(l^3);
i3=6*E*I/(l^2);
i4=4*E*I/l;
i5=2*E*I/l;
k=zeros(6);
k(1,1)=i1;k(1,4)=-i1;
k(2,2)=i2;k(2,3)=i3;k(2,5)=-i2;k(2,6)=i3;
k(3,3)=i4;k(3,5)=-i3;k(3,6)=i5;
k(4,4)=i1;
k(5,5)=i2;k(5,6)=-i3;
k(6,6)=i4;
k(3,2)=k(2,3);
k(4,1)=k(1,4);
k(5,2)=k(2,5);k(5,3)=k(3,5);
k(6,2)=k(2,6);k(6,3)=k(3,6);k(6,5)=k(5,6);
y=k;
end
附录2ANSYS建模计算命令流
finish
/clear
/prep7
B=0.120$H=0.200$E=210000000
et,1,beam3
mp,ex,1,E
mp,prxy,1,0.3
r,1,B*H,B*H*H*H/12,H
k,1
k,2,0,5
k,3,1.6304,6.3681
k,4,4,5
k,5,4,1
k,6,4,-2
*set,i
*do,i,1,5
l,i,i+1
*enddo
lesize,all,0.5
latt,1,1,1
lmesh,all
dk,1,all
dk,6,all
fk,3,fy,-40
fk,5,fx,-50
lsel,s,,,1
esll,s
sfbeam,all,1,pres,10
allsel,all
dtran
ftran
sftran
/pbc,all,,2
/psf,pres,norm,2,0,1
eplot
/solu
solve
/post1
/pbc,u,,1!
显示支座约束符号,并图形显示变形
pldisp,1!
将当前主要结果列表显示
presol,elem
!
/pnum,sval,1
etable,mi,smisc,6
etable,mj,smisc,12
plls,mi,mj,-1!
弯矩图kN.m
etable,qi,smisc,2
etable,qj,smisc,8
plls,qi,qj,-1!
剪力图kN
etable,ni,smisc,1
etable,nj,smisc,7
plls,ni,nj,1!
轴力图kN
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 基于 MATLAB 平面 刚架 静力 分析 报告