matlab_矩阵位移法编程_结构力学

矩阵位移法编程大作业

一、编制原理

本程序的原理是基于结构力学矩阵位移法原理，以结构结点位移作基本未知量，将要分析的结构拆成已知节点力—结点力位移关系的单跨梁集合，通过强令结构发生待定的基本未知位移，在各个单跨梁受力分析结果的基础上通过保证结构平衡建立位移法的线性方程组，从而求得基本未知量。

二、程序说明

本程序是计算3层11跨框架右侧结点的位移和弯矩的程序，编译过程是按照矩阵位移法的先处理法进行的。首先将结构杆件的交汇点作为结点，共有36个结点和108个位移编号，然后根据梁、柱、斜杆的不同分别建立单元刚度矩阵，然后转换为整体坐标系下的刚度矩阵，然后将所有杆件的单元刚度矩阵整合成为总体刚度矩阵，在进行整合时连续运用for函数，最终形成108阶的总体刚度矩阵。然后通过对荷载的分析自己确定出荷载矩阵，直接写进程序。这样就可以把36个结点的108个位移求得，然后再利用各个单元的单元刚度矩阵和所得的位移求得单元杆件的内力。

离散化编号如下图：

三、算法流程

结构离散化编号单元分析建立梁、柱、斜杆的局部

坐标系下的单元刚度矩阵

确定梁、柱、斜杆在整体坐标系下的刚度矩阵先处理法把各个杆件的单元刚度矩阵整合成总体刚度矩阵

确定综合结点荷载矩阵建立方程，

求解位移

利用杆件单元刚度矩

阵和所求位移求内力

四、源代码

%结构力学大作业 3层11跨框架矩阵位移法编程王贝 091210211

h=input('输入单层高h：');

L=input('输入单跨度L：');

EIc=input('输入柱子的抗弯刚度EIc：');

EAc=input('输入柱子的抗压刚度EAc：');

EIb=input('输入梁的抗弯刚度EIb：');

EAb=input('输入梁的抗压刚度EAb：');

EIo=input('输入斜杆的抗弯刚度EIo：');

EAo=input('输入斜杆的抗压刚度EAo：');

q=input('输入侧向均布荷载集度q:');

T1=[1,0,0,0,0,0;

0,1,0,0,0,0;

0,0,1,0,0,0;

0,0,0,1,0,0;

0,0,0,0,1,0;

0,0,0,0,0,1];%角度为0°的转换矩阵

T2=[0,1,0,0,0,0;

-1,0,0,0,0,0;

0,0,1,0,0,0;

0,0,0,0,1,0;

0,0,0,-1,0,0;

0,0,0,0,0,1];%角度为90°的转换矩阵

x=atan(h/L);

T=[cos(x),sin(x),0,0,0,0;

-sin(x),cos(x),0,0,0,0;

0,0,1,0,0,0;

0,0,0,cos(x),sin(x),0;

0,0,0,-sin(x),cos(x),0;

0,0,0,0,0,1];%斜杆的转换矩阵

T3=T;

%梁的单元刚度矩阵

kb0=[EAb/L 0 0 -EAb/L 0 0;

0 12*EIb/(L*L*L) 6*EIb/(L*L) 0 -12*EIb/(L*L*L) 6*EIb/(L*L);

0 6*EIb/(L*L) 4*EIb/L 0 -6*EIb/(L*L) 2*EIb/L;

-EAb/L 0 0 EAb/L 0 0;

0 -12*EIb/(L*L*L) -6*EIb/(L*L) 0 12*EIb/(L*L*L) -6*EIb/(L*L);

0 6*EIb/(L*L) 2*EIb/L 0 -6*EIb/(L*L) 4*EIb/L];

%柱子的单元刚度矩阵

kc0=[EAc/h 0 0 -EAc/h 0 0;

0 12*EIc/(h*h*h) 6*EIc/(h*h) 0 -12*EIc/(h*h*h) 6*EIc/(h*h);

0 6*EIc/(h*h) 4*EIc/h 0 -6*EIc/(h*h) 2*EIc/h;

-EAc/h 0 0 EAc/h 0 0;

0 -12*EIc/(h*h*h) -6*EIc/(h*h) 0 12*EIc/(h*h*h) -6*EIc/(h*h);

0 6*EIc/(h*h) 2*EIc/h 0 -6*EIc/(h*h) 4*EIc/h;];

%斜杆的单元刚度矩阵

H=sqrt(h*h+L*L);

ko0=[EAo/H 0 0 -EAo/H 0 0;

0 12*EIo/(H*H*H) 6*EIo/(H*H) 0 -12*EIo/(H*H*H) 6*EIo/(H*H);

0 6*EIo/(H*H) 4*EIo/H 0 -6*EIo/(H*H) 2*EIo/H;

-EAo/H 0 0 EAo/H 0 0;

0 -12*EIo/(H*H*H) -6*EIo/(H*H) 0 12*EIo/(H*H*H) -6*EIo/(H*H);

0 6*EIo/(H*H) 2*EIo/H 0 -6*EIo/(H*H) 4*EIo/H];

kb=T1'*kb0*T1;%总体坐标下梁的单元刚度矩阵

kc=T2'*kc0*T2;%总体坐标下柱子的单元刚度矩阵

ko=T3'*ko0*T3;%总体坐标斜杆的单元刚度矩阵

X=zeros(108,108);Y=zeros(108,108);Z=zeros(108,108);%定义108阶0矩阵K1=zeros(108,108);K2=zeros(108,108);K3=zeros(108,108);

K4=zeros(108,108);K5=zeros(108,108);K6=zeros(108,108);

K7=zeros(108,108);K8=zeros(108,108);K9=zeros(108,108);

%把梁杆单元矩阵整合到总体刚度矩阵的循环语句

for ii=1:11

X(3*ii-2:3*ii+3,3*ii-2:3*ii+3)=kb;

K1=K1+X;X=zeros(108,108);

end

for ii=13:23

Y(3*ii-2:3*ii+3,3*ii-2:3*ii+3)=kb;

K1=K1+Y;Y=zeros(108,108);

end

for ii=25:35

Z(3*ii-2:3*ii+3,3*ii-2:3*ii+3)=kb;

K1=K1+Z;Z=zeros(108,108);

end

%把柱杆单元矩阵整合到总体刚度矩阵的循环语句

for jj=1:36

K2(3*jj-2:3*jj,3*jj-2:3*jj)=kc(4:6,4:6);

end

for jj=1:24

K3(3*jj-2:3*jj,3*jj-2:3*jj)=kc(1:3,1:3);

end

for jj=1:24

K4(3*jj-2:3*jj,3*jj+34:3*jj+36)=kc(1:3,4:6);

end