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

1、矩阵位移法编程大作业姓名:子一、编程原理本程序的原理是基于结构力学矩阵位移法原理，以结构结点位移作基本未知量，将要分析的结构拆成已知节点力一结点力位移关系的单跨 梁集合，通过强令结构发生待定的基本未知位移，在各个单跨梁受力分 析结果的基础上通过保证结构平衡建立位移法的线性方程组，从而求得基本未知量。二、程序说明本程序是计算10个节间距的悬索-拱组合体系主塔顶节点水平位 移、主塔底截面弯矩、拱顶节点竖向位移、拱顶截面弯矩和轴力的程序。 首先将各杆件的交汇点作为结点，共有 20个结点，51个位移，然后根 据不同结构单元分别建立单元刚度矩阵，然后转换为整体坐标系下的刚度矩阵，然后将所有杆件的单元刚度

2、矩阵整合成为总体刚度矩阵，在进 行整合时连续运用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:）;

3、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,

4、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=

5、(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(11)=5*L/2;y=zeros(10,1);for i=1:10y(i)=atan(f2(i+1)-f2(i)/L);endT3=zeros(4,40);for i=1:10T3(1:4,4*i-3:4*i)

6、=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:10s2(i)=sqrt(f2(i+1)-f2(i)A2+LA2);endKbO2=zeros(4,40);for i=1:10KbO2(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/

7、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:10z(i)=atan(f1(i+1)-f1(i)/L);endT4=zeros(6,60);for i=1:10T4(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

8、,-sin(z(i),cos(z(i),0;0,0,0,0,0,1;end %拱的转换矩阵s3=zeros(10,1);for i=1:10s3(i)=sqrt(f1(i+1)-f1(i)A2+LA2);endKaO=zeros(6,60);for i=1:10KaO(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

9、/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) 012*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(

10、9,1);for i=1:9s4(i)=f2(i+1)-f1(i+1);endKtO=zeros(4,36);for i=1:9KtO(1:4,4*i-3:4*i)=(EAt/s4(i)*1 0 -1 0;0 0 0 0;-1 0 1 0;0 0 0 0;end %吊杆的单元刚度矩阵Kc=T1*KcO*T1;%总体坐标下主塔的单元刚度矩阵Kb1=T2*KbO1*T2;Kb11=T21*KbO1*T21;%总体坐标下斜索的单元刚度矩阵Kb2=zeros(4,40);for i=1:10T3O=T3(1:4,4*i-3:4*i);Kb2(1:4,4*i-3:4*i)=T3O*KbO2(1:4,4*

11、i-3:4*i)*T3O;end %总体坐标下悬索的单元刚度矩阵Ka=zeros(6,60);for i=1:10T4O=T4(6*i-5:6*i,6*i-5:6*i);Ka(1:6,6*i-5:6*i)=T4O*KaO(1:6,6*i-5:6*i)*T4O;end %总体坐标下拱的单元刚度矩阵Kt=zeros(4,36);for i=1:9KtOO=KtO(1:4,4*i-3:4*i);Kt(1:4,4*i-3:4*i)=T5*KtOO*T5;end %总体坐标下吊杆的单元刚度矩阵%定义51阶0矩阵K1=zeros(51,51);K2=zeros(51,51);K3=zeros(51,51)

12、;K4=zeros(51,51);K5=zero s(51,51);X=zeros(51,51);Y=zeros(51,51);Z=zeros(51,51);咐巴主塔整合到整体刚度矩阵中：K1(1:3,1:3)=KcO(4:6,4:6);K1(22:24,22:24)=KcO(4:6,4:6);引巴斜索整合到整体刚度矩阵中：K2(1:2,1:2)=Kb1(3:4,3:4);K2(22:23,22:23)=Kb11(1:2,1:2);引巴悬索整合到整彳刚度矩阵中：K3(1:2,1:2)=KbO2(1:2,1:2);K3(1:2,4:5)=KbO2(1:2,3:4);for i=2:10X(2*i

13、:2*i+3,2*i:2*i+3)=KbO2(1:4,4*i-3:4*i);K3=K3+X;end引巴拱整合到整体刚度矩阵中：K4(25:27,25:27)=KaO(4:6,4:6);K4(49:51,49:51)=KaO(1:3,55:57);for i=2:9Y(3*i+19:3*i+24,3*i+19:3*i+24)=KaO(1:6,6*i-5:6*i);K4=K4+Y;end咐巴吊杆整合到整体刚度矩阵中：for i=1:9Z(2*i+2:2*i+3,2*i+2:2*i+3)=KtO(1:2,1:2);Z(2*i+2:2*i+3,3*i+22:3*i+23)=KtO(1:2,3:4);Z

14、(3*i+22:3*i+23,2*i+2:2*i+3)=KtO(3:4,1:2);Z(3*i+22:3*i+23,3*i+22:3*i+23)=KtO(3:4,3:4);K5=K5+Z;endK=K1+K2+K3+K4+K5;%荷载矩阵：P=zeros(51,1);P(26,1)=-q*L/(2*cos(s3(1);P(27,1)=q*L*L/(12*cos(s3(1);P(50,1)=-q*L/(2*cos(s3(10);P(51,1)=-q*L*L/(12*cos(s3(10);for i=2:9P0=zeros(51,1);P0(3*i+20,1)=-q*L/(2*cos(s3(i);P

15、0(3*i+21,1)=-q*L*L/(12*cos(s3(i);P0(3*i+23,1)=-q*L/(2*cos(s3(i);P0(3*i+24,1)=q*L*L/(12*cos(s3(i);P=P+P0;endA=KP; %吉构的位移%主塔底截面的弯矩：Ac(4:6,1)=A(1:3,1);Bc=KcO*Ac;Mc=Bc(3,1);%拱顶截面的弯矩和轴力：Aa=A(34:39,1);KaO17=KaO(1:6,25:30);Ba=KaO17*Aa;Ma=Ba(6,1);Fa=Ba(4,1);%输出结果fprintf( 主塔顶结点的水平位移 %fn ,A(1,1);fprintf(主塔底截面的弯矩 %fn ,Mc);fprintf(拱顶结点的竖向位