基于MATLAB的平面刚架有限元分析,基于MATLAB的平面刚架静力分析

工程计算实践作业

基于MATLAB的平面刚架静力分析

为了进一步理解有限元方法计算的过程，本文根据矩阵位移法的基本原理应用MATLAB编制计算程序对以平面刚架结构进行了静力分析。本文还利用ANSYS大型商用有限元分析软件对矩阵位移法的计算结果进行校核，发现两者计算结果相当吻合，验证了计算结果的可靠性。

一、 问题描述

如图1所示的平面刚架，各杆件的材料及截面均相同，E=210GPa，截面为0.120.2m的实心矩形，现要求解荷载作用下刚架的位移和内力。

图1

二、矩阵位移法计算程序编制

为编制程序方便考虑，本文计算中采用“先处理法”。具体的计算步骤如下。

(1) 对结构进行离散化，对结点和单元进行编号，建立结构(整体)坐标系和单元(局部)坐标系，并对结点位移进行编号；

(2) 对结点位移分量进行编码，形成单元定位向量；

(3) 建立按结构整体编码顺序排列的结点位移列向量，计算固端力、等效结点荷载 及综合结点荷载列向量；

(4) 计算个单元局部坐标系的刚度矩阵，通过坐标变换矩阵形成整体坐标系下的单元刚度矩阵 ；

(5) 利用单元定位向量形成结构刚度矩阵；

(6) 按式 求解未知结点位移；

(7) 计算各单元的杆端力。

根据上述步骤编制了平面刚架的分析程序。程序中单元刚度矩阵按下式计算。

转换矩阵则按下式计算。

计算程序框图如图2所示，具体的程序代码见附录1。

图2 MATLAB矩阵分析法程序框图

三、解题步骤

取整体坐标系如图3所示，对结构进行离散化，对结点和单元进行编号如图4所示，局部坐标系用单元中箭头的方向表示，原始数据如下：

图3 图4

刚架结点输入矩阵为，

x=[0 0;0 -5;1.63 -6.37;4 -5;4 -1;4 2];

各单元定位向量为，

locvec1=[1 2 3 0 0 0];

locvec2=[1 2 3 4 5 6];

locvec3=[4 5 6 7 8 9];

locvec4=[7 8 9 10 11 12];

locvec5=[10 11 12 0 0 0];

输入截面参数，

E=2.1e11;%E=210GPa

a=0.12; %矩形截面长0.12m

b=0.2; %矩形截面宽0.2m

输入整体坐标系下各单元结点荷载列阵，

f(1,:)=zeros(1,6);

f(2,:)=[0 0 0 0 40e3 0];

f(3,:)=zeros(1,6);

f(4,:)=[0 0 0 -50e3 0 0];

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(kNm)

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(kNm)

-37312.596

-2092.83

26668.34

-44999.8

51249.01

四、计算结果校核

在ANSYS中使用beam3单元，按照如图4所示的离散结构建立平面刚架模型施加约束和荷载，得到的有限元模型如图5所示。计算分析后得到结构的轴力图、剪力图和弯矩图如图6、7、8所示，命令流见附录2。

图5 有限元模型

图6 轴力图(kN)

图7 剪力图(kN)

图8 弯矩图(kNm)

从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(kNm)

-37.313

-37.307

-0.006

2

Fx(kN)

-17.917

-17.920

0.003

Fy(kN)

17.507

17.503

0.004

Mz(kNm)

1.898

1.909

-0.011

3

Fx(kN)

-17.917

-17.920

0.003

Fy(kN)

-22.493

-22.497

0.004

Mz(kNm)

-2.093

-2.110

0.018

4

Fx(kN)

-17.917

-17.920

0.003

Fy(kN)

-22.493

-22.497

0.004

Mz(kNm)

26.668

26.682

-0.014

5

Fx(kN)

-32.083

-32.080

-0.003

Fy(kN)

-22.493

-22.497

0.004

Mz(kNm)

-45.000

-44.999

-0.001

6

Fx(kN)

32.083

32.080

0.003

Fy(kN)

-22.493

-22.497

0.004

Mz(kNm)

51.249

51.239

0.010

由表2、表3的结果对比可知，两种方法的计算结果十分接近，误差均可以忽略不计，从而验证了计算结果的可靠性与准确性。

四、结论

通过对一个平面刚架静力分析问题的求解，本文编制的平面刚架静力分析程序计算结果与有限元软件ANSYS计算结果校核后，发现两者计算结果十分接近，误差可忽略不计，说明该程序计算结果的可靠性与精确性。

附录1 矩阵位移法计算程序

pmgj.m 计算主程序

clear

clc

%结点坐标

x=[0 0;0 -5;1.63 -6.37;4 -5;4 -1;4 2];

% 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=[1 2 3 0 0 0];

locvec2=[1 2 3 4 5 6];

locvec3=[4 5 6 7 8 9];

locvec4=[7 8 9 10 11 12];

locvec5=[10 11 12 0 0 0];

%刚架的材料特性 截面特性

E=2.1e11;%E=210GPa

a=0.12; %矩形截面长0.12m

b=0.2; %矩形截面宽0.2m

A=a*b;

I=(a*b^3)/12;

clear a b

% 单元长度计算

for i=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,:)=[0 0 0 0 40e3 0];

f(3,:)=zeros(1,6);

f(4,:)=[0 0 0 -50e3 0 0];

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);

clear NonLog ij IJ

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);

clear NonLog ij IJ

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);

clear NonLog ij IJ

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);

clear NonLog ij IJ

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);

clear NonLog ij IJ

%节点位移

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

function y=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均采用国际单位制 Pa m2 m m4

function y=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

附录2 ANSYS建模计算命令流

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

17

展开阅读全文