结构力学大作业(矩阵位移法)

矩阵位移法编程大作业

姓名：

学号：

一、编程原理

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

二、程序说明

本程序是计算10个节间距的悬索-拱组合体系主塔顶节点水平位移、主塔底截面弯矩、拱顶节点竖向位移、拱顶截面弯矩和轴力的程序。首先将各杆件的交汇点作为结点，共有20个结点，51个位移，然后根据不同结构单元分别建立单元刚度矩阵，然后转换为整体坐标系下的刚度矩阵，然后将所有杆件的单元刚度矩阵整合成为总体刚度矩阵，在进行整合时连续运用for函数，最终形成51阶的总体刚度矩阵。然后通过对荷载的分析确定出荷载矩阵，直接写进程序。这样就可以把20个结点的51个位移求得，然后再利用各个单元的单元刚度矩阵和所得的位移求得单元杆件的内力。

三、算法流程

建立各单位在局部结构离散化编号进行单元分析坐标系下的单位刚度方程

确定各单位在总体将单元刚度矩阵集合确定综合结点坐标系下的

单元矩阵方程成总体刚度矩阵点荷载矩阵

建立方程利用杆件单元刚度矩阵输出结果

求解位移和所求位移求内力

结束

四、源代码

L=input('输入单节间L：');

EIc=input('主塔的抗弯刚度EIc:');

EAc=input('主塔的抗压刚度EAc:');

EAb=input('悬索和斜索的抗拉刚度EAb:');

EAt=input('吊杆的抗拉刚度EAt:');

EIa=input('拱的抗弯刚度EIa:');

EAa=input('拱的抗压刚度EAa:');

q=input('拱上沿轴向均布荷载集度q:');

T1=[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;];%主塔的转换矩阵

h=(5*L)/2;

KcO=[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;];%主塔的单元刚度矩阵

x=atan(2*L/h);

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

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

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

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

y=-atan(2*L/h);

T21=[cos(y),sin(y),0,0;

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

0,0,cos(y),sin(y);

0,0,-sin(y),cos(y);];%斜索的转换矩阵

s1=sqrt(2*L*2*L+h*h);

KbO1=(EAb/s1)*[1 0 -1 0;

0 0 0 0;

-1 0 1 0;

0 0 0 0;];%斜索的单元刚度矩阵

f2(1)=5*L/2;f2(2)=58*L/25;f2(3)=109*L/50;f(4)=52*L/25;f2(5)=101*L/50;f2 (6)=2*L;f2(7)=101*L/50;f2(8)=52*L/25;f2(9)=109*L/50;f2(10)=58*L/25;f2(1 1)=5*L/2;

y=zeros(10,1);

for i=1:10

y(i)=atan((f2(i+1)-f2(i))/L);

end

T3=zeros(4,40);

for i=1:10

T3(1:4,4*i-3:4*i)=[cos(y(i)),sin(y(i)),0,0;

-sin(y(i)),cos(y(i)),0,0;

0,0,cos(y(i)),sin(y(i));

0,0,-sin(y(i)),cos(y(i));];

end%悬索的转换矩阵

s2=zeros(10,1);

for i=1:10

s2(i)=sqrt((f2(i+1)-f2(i))^2+L^2);

end

KbO2=zeros(4,40);

KbO2(1:4,4*i-3:4*i)=(EAb/s2(i))*[1 0 -1 0;

0 0 0 0;

-1 0 1 0;

0 0 0 0;];

end%悬索的单元刚度矩阵

f1(1)=0;f1(2)=9*L/20;f1(3)=4*L/5;f1(4)=21*L/20;f1(5)=6*L/5;f1(6)=5*L/4; f1(7)=6*L/5;f1(8)=21*L/20;f1(9)=4*L/5;f1(10)=9*L/20;f1(11)=0;

z=zeros(10,1);

for i=1:10

z(i)=atan((f1(i+1)-f1(i))/L);

end

T4=zeros(6,60);

for i=1:10

T4(6*i-5:6*i,6*i-5:6*i)=[cos(z(i)),sin(z(i)),0,0,0,0;

-sin(z(i)),cos(z(i)),0,0,0,0;

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

0,0,0,cos(z(i)),sin(z(i)),0;

0,0,0,-sin(z(i)),cos(z(i)),0;

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

end%拱的转换矩阵

s3=zeros(10,1);

for i=1:10

s3(i)=sqrt((f1(i+1)-f1(i))^2+L^2);

end

KaO=zeros(6,60);

for i=1:10

KaO(1:6,6*i-5:6*i)=[EAa/s3(i) 0 0 -EAa/s3(i) 0 0;

0 12*EIa/(s3(i)*s3(i)*s3(i)) 6*EIa/(s3(i)*s3(i)) 0

-12*EIa/(s3(i)*s3(i)*s3(i)) 6*EIa/(s3(i)*s3(i));

0 6*EIa/(s3(i)*s3(i)) 4*EIa/s3(i) 0 -6*EIa/(s3(i)*s3(i)) 2*EIa/s3(i);

-EAa/s3(i) 0 0 EAa/s3(i) 0 0;

0 -12*EIa/(s3(i)*s3(i)*s3(i)) -6*EIa/(s3(i)*s3(i)) 0

12*EIa/(s3(i)*s3(i)*s3(i)) -6*EIa/(s3(i)*s3(i));

0 6*EIa/(s3(i)*s3(i)) 2*EIa/s3(i) 0 -6*EIa/(s3(i)*s3(i)) 4*EIa/s3(i);]; end%拱的单元刚度矩阵

T5=[0 1 0 0;

-1 0 0 0;

0 0 0 1;

0 0 -1 0;];%吊杆的转换矩阵

s4=zeros(9,1);