matlab矩陣位移法編程結(jié)構(gòu)力學(xué)

1、矩陣位移法編程大作業(yè)(091210211)、編制原理本程序的原理是基于結(jié)構(gòu)力學(xué)矩陣位移法原理，以結(jié)構(gòu)結(jié)點(diǎn)位移作基本未知 量，將要分析的結(jié)構(gòu)拆成已知節(jié)點(diǎn)力一結(jié)點(diǎn)力位移關(guān)系的單跨梁集合，通過強(qiáng)令結(jié)構(gòu)發(fā)生待定的基本未知位移，在各個單跨梁受力分析結(jié)果的基礎(chǔ)上通過保證結(jié) 構(gòu)平衡建立位移法的線性方程組，從而求得基本未知量。二、程序說明本程序是計(jì)算3層11跨框架右側(cè)結(jié)點(diǎn)的位移和彎矩的程序，編譯過程是按照矩陣位移法的先處理法進(jìn)行的。首先將結(jié)構(gòu)桿件的交匯點(diǎn)作為結(jié)點(diǎn)，共有36個結(jié)點(diǎn)和108個位移編號，然后根據(jù)梁、柱、斜桿的不同分別建立單元剛度矩陣， 然后轉(zhuǎn)換為整體坐標(biāo)系下的剛度矩陣，然后將所有桿件的單元剛度矩陣整

2、合成為 總體剛度矩陣，在進(jìn)行整合時連續(xù)運(yùn)用 for函數(shù)，最終形成108階的總體剛度矩 陣。然后通過對荷載的分析自己確定出荷載矩陣，直接寫進(jìn)程序。這樣就可以把36個結(jié)點(diǎn)的108個位移求得，然后再利用各個單元的單元剛度矩陣和所得的位 移求得單元桿件的內(nèi)力。離散化編號如下圖：*2 LEO黒LJ6I胡LE245US3範(fàn)L64打歷5陽L8637L6738 LJEfl38LE9L47148L50151L53L57L533C厲31L3833L403*L413ELC25L44話L152?L*38 L?l124125L2JL23L29L3005LU19L15藥LIE21L1722L182L132* 120勺12

3、1r U213L2315 L9L112LSLiL5L6L7L9L'OL11L1217二、算法流程輸出結(jié)果a結(jié)束四、源代碼%吉構(gòu)力學(xué)大作業(yè)3層11跨框架矩陣位移法編程 王貝091210211 h=input('輸入單層高h(yuǎn)：');L=input('輸入單跨度L：');Elc=input('輸入柱子的抗彎剛度Elc :');EAc=input('輸入柱子的抗壓剛度EAc ');Elb=input('輸入梁的抗彎剛度Elb :');EAb=input('輸入梁的抗壓剛度EAb ');Elo=inp

4、ut('輸入斜桿的抗彎剛度Elo :');EAo=input('輸入斜桿的抗壓剛度EAo ');q=input('輸入側(cè)向均布荷載集度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°的轉(zhuǎn)換矩陣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°的轉(zhuǎn)換矩陣x=ata n( h/L);T=cos

5、(x),si n(x),0,0,0,0;-si n(x),cos(x),0,0,0,0;0,0,1,0,0,0;0,0,0,cos(x),si n( x),0;0,0,0,-s in (x),cos(x),0;0,0,0,0,0,1;%斜桿的轉(zhuǎn)換矩陣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/(

6、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*

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

8、);0 6*EIo/(H*H) 2*EIo/H 0 -6*EIo/(H*H) 4*EIo/H; kb=T1'*kb0*T1; %總體坐標(biāo)下梁的單元剛度矩陣 kc=T2'*kc0*T2; %總體坐標(biāo)下柱子的單元剛度矩陣 ko=T3'*ko0*T3; %總體坐標(biāo)斜桿的單元剛度矩陣 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,1

9、08);K6=zeros(108,108); K7=zeros(108,108);K8=zeros(108,108);K9=zeros(108,108); %把梁桿單元矩陣整合到總體剛度矩陣的循環(huán)語句 for ii=1:11X(3*ii-2:3*ii+3,3*ii-2:3*ii+3)=kb;K1=K1+X;X=zeros(108,108);end for ii=13:23Y(3*ii-2:3*ii+3,3*ii-2:3*ii+3)=kb; K1=K1+Y;Y=zeros(108,108);end for ii=25:35Z(3*ii-2:3*ii+3,3*ii-2:3*ii+3)=kb; K1

10、=K1+Z;Z=zeros(108,108);end %把柱桿單元矩陣整合到總體剛度矩陣的循環(huán)語句 for jj=1:36K2(3*jj-2:3*jj,3*jj-2:3*jj)=kc(4:6,4:6);end for jj=1:24K3(3*jj-2:3*jj,3*jj-2:3*jj)=kc(1:3,1:3);end for jj=1:24K4(3*jj-2:3*jj,3*jj+34:3*jj+36)=kc(1:3,4:6);endfor jj=1:24 K5(3*jj+34:3*jj+36,3*jj-2:3*jj)=kc(4:6,1:3);end %把斜桿桿單元矩陣整合到總體剛度矩陣的循環(huán)語

11、句 for gg=3:12:27K6(3*gg-2:3*gg,3*gg-2:3*gg)=ko(4:6,4:6);endfor gg=2:12:14K7(3*gg-2:3*gg,3*gg-2:3*gg)=ko(1:3,1:3);endfor gg=2:12:14K8(3*gg-2:3*gg,3*gg+37:3*gg+39)=ko(1:3,4:6);endfor gg=2:12:14K9(3*gg+37:3*gg+39,3*gg-2:3*gg)=ko(4:6,1:3);endK=K1+K2+K3+K4+K5+K6+K7+K8+%總體岡 H度矩陣P=zeros(108,1);P(1,1)=h*q;P

12、(37,1)=h*q;P(73,1)=h*q/2;P(75,1)=q*h*h/12; A=KP; %結(jié)構(gòu)位移B1=kb*A(103:108,1);B2=kb*A(67:72,1);B3=kb*A(31:36,1); D1=zeros(6,1);D1(1:3,1)=A(70:72,1);D1(4:6,1)=A(106:108,1);D2=zeros(6,1);D2(1:3,1)=A(34:36,1);D2(4:6,1)=A(70:72,1); C1=kc*D1;C2=kc*D2;C3=kc(4:6,4:6)*A(34:36,1);M1(1,1)=B1(6,1);M2(1,1)=C1(6,1);M

13、3(1,1)=B2(6,1);M4(1,1)=C2(6,1);M5(1,1)=C1(3,1);M6(1,1)=B3(6,1);M7(1,1)=C3(3,1);M8(1,1)=C2(3,1);for i=1:3m=36*i-2;fprintf('第層最右側(cè)節(jié)點(diǎn)的位移是dn',i,A(m,1)endfprintf( ' 第1層最右側(cè)節(jié)點(diǎn)左側(cè)桿的彎矩是 %fn' ,M6)fprintf( ' 第1層最右側(cè)節(jié)點(diǎn)下側(cè)桿的彎矩是 %fn' ,M7)fprintf( ' 第1層最右側(cè)節(jié)點(diǎn)上側(cè)桿的彎矩是 %fn' ,M8)fprintf( 

14、9; 第2層最右側(cè)節(jié)點(diǎn)左側(cè)桿的彎矩是 %fn' ,M3)fprintf( ' 第2層最右側(cè)節(jié)點(diǎn)下側(cè)桿的彎矩是 %fn' ,M4)fprintf( ' 第2層最右側(cè)節(jié)點(diǎn)上側(cè)桿的彎矩是 %fn' ,M5)fprintf( ' 第3層最右側(cè)節(jié)點(diǎn)左側(cè)桿的彎矩是 %fn' ,M1)fprintf( ' 第3層最右側(cè)節(jié)點(diǎn)下側(cè)桿的彎矩是 %fn' ,M2)五、試算算例 輸入數(shù)據(jù)：輸入單層高 h：1 輸入單跨度 L ：1 輸入柱子的抗彎剛度 EIc：1 輸入柱子的抗壓剛度 EAc： 1 輸入梁的抗彎剛度 EIb：1 輸入梁的抗壓剛度 EAb：1 輸入斜桿的抗彎剛度 EIo：1 輸入斜桿的抗壓剛度 EAo： 1 輸入側(cè)向均布荷載集度 q:1 計(jì)算結(jié)果：第 1 層最右側(cè)節(jié)點(diǎn)的位移是 -6.219850e-003 第 2 層最右側(cè)節(jié)點(diǎn)的位移是 -2.152659e-002 第 3 層最右側(cè)節(jié)點(diǎn)的位移是 -4.131873e-002 第 1 層最右側(cè)節(jié)點(diǎn)左