結(jié)構(gòu)力學(xué)大作業(yè)矩陣位移法

1、矩陣位移法編程大作業(yè)姓名:子一、編程原理本程序的原理是基于結(jié)構(gòu)力學(xué)矩陣位移法原理，以結(jié)構(gòu)結(jié)點位移作基本未知量，將要分析的結(jié)構(gòu)拆成已知節(jié)點力一結(jié)點力位移關(guān)系的單跨 梁集合，通過強令結(jié)構(gòu)發(fā)生待定的基本未知位移，在各個單跨梁受力分 析結(jié)果的基礎(chǔ)上通過保證結(jié)構(gòu)平衡建立位移法的線性方程組，從而求得基本未知量。二、程序說明本程序是計算10個節(jié)間距的懸索-拱組合體系主塔頂節(jié)點水平位 移、主塔底截面彎矩、拱頂節(jié)點豎向位移、拱頂截面彎矩和軸力的程序。 首先將各桿件的交匯點作為結(jié)點，共有 20個結(jié)點，51個位移，然后根 據(jù)不同結(jié)構(gòu)單元分別建立單元剛度矩陣，然后轉(zhuǎn)換為整體坐標(biāo)系下的剛度矩陣，然后將所有桿件的單元剛度

2、矩陣整合成為總體剛度矩陣，在進 行整合時連續(xù)運用for函數(shù)，最終形成51階的總體剛度矩陣。然后通過 對荷載的分析確定出荷載矩陣，直接寫進程序。這樣就可以把 20個結(jié) 點的51個位移求得，然后再利用各個單元的單元剛度矩陣和所得的位 移求得單元桿件的內(nèi)力。三、算法流程結(jié)束四、源代碼L=input（輸入單節(jié)間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;吹塔的轉(zhuǎn)換矩陣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);%斗索的轉(zhuǎn)換矩陣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 %懸索的轉(zhuǎn)換矩陣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 %拱的轉(zhuǎn)換矩陣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;%吊桿的轉(zhuǎn)換矩陣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;%總體坐標(biāo)下主塔的單元剛度矩陣Kb1=T2*KbO1*T2;Kb11=T21*KbO1*T21;%總體坐標(biāo)下斜索的單元剛度矩陣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 %總體坐標(biāo)下懸索的單元剛度矩陣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 %總體坐標(biāo)下拱的單元剛度矩陣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 %總體坐標(biāo)下吊桿的單元剛度矩陣%定義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; %吉構(gòu)的位移%主塔底截面的彎矩：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);%輸出結(jié)果fprintf( 主塔頂結(jié)點的水平位移 %fn ,A(1,1);fprintf(主塔底截面的彎矩 %fn ,Mc);fprintf(拱頂結(jié)點的豎向位