基本平面剛架MATLAB程序

1、% 平面剛架MATLAB程序% 2003.9.16 2008.4.1 2009.10 2011.10%*% 變量說明% NPOIN NELEM NVFIX NFPOIN NFPRES% 總結(jié)點(diǎn)數(shù),單元數(shù), 約束個(gè)數(shù), 受力結(jié)點(diǎn)數(shù), 非結(jié)點(diǎn)力數(shù)% COORD LNODS YOUNG% 結(jié)構(gòu)節(jié)點(diǎn)坐標(biāo)數(shù)組, 單元定義數(shù)組, 彈性模量% FPOIN FPRES FORCE FIXED% 結(jié)點(diǎn)力數(shù)組，非結(jié)點(diǎn)力數(shù)組，總體荷載向量, 約束信息數(shù)組% HK DISP% 總體剛度矩陣,結(jié)點(diǎn)位移向量%*format short e %設(shè)定輸出類型clear %清除內(nèi)存變量FP1=fopen('6-6.tx

2、t','rt') %打開初始數(shù)據(jù)文件%讀入控制數(shù)據(jù) NELEM=fscanf(FP1,'%d',1); %單元數(shù)NPOIN=fscanf(FP1,'%d',1); %結(jié)點(diǎn)數(shù)NVFIX=fscanf(FP1,'%d',1); %約束數(shù)NFPOIN=fscanf(FP1,'%d',1); %作用荷載的結(jié)點(diǎn)個(gè)數(shù)NFPRES=fscanf(FP1,'%d',1); %非結(jié)點(diǎn)荷載數(shù)YOUNG=fscanf(FP1,'%f',1); %彈性模量% 讀取結(jié)構(gòu)信息LNODS=fscanf(F

3、P1,'%f',4,NELEM)' % 單元定義： 左、右結(jié)點(diǎn)號，面積，慣性矩（共計(jì) NELEM組）COORD=fscanf(FP1,'%f',2,NPOIN)' % 坐標(biāo)： x,y坐標(biāo)（共計(jì) NPOIN 組）FPOIN=fscanf(FP1,'%f',4,NFPOIN)' % 節(jié)點(diǎn)力（共計(jì) NFPOIN 組）：結(jié)點(diǎn)號、X方向力（向右正），% Y方向力（向上正），M力偶（逆時(shí)針正）FPRES=fscanf(FP1,'%f',4,NFPRES)' % 均布力（共計(jì)% NFPRES 組）：單元號、荷載類

4、型、荷載大小、距離左端長度FIXED=fscanf(FP1,'%f',NVFIX)' % 約束信息：約束對應(yīng)的位移編碼（共計(jì) NVFIX 組）%-HK=zeros(3*NPOIN,3*NPOIN); % 張成總剛矩陣并清零FORCE=zeros(3*NPOIN,1); % 張成總荷載向量并清零%形成總剛for i=1:NELEM % 對單元個(gè)數(shù)循環(huán) % 生成局部單剛(局部坐標(biāo)) 右手坐標(biāo)系 EK=ele_EK(i,LNODS,COORD,YOUNG); T=zbzh(i,LNODS,COORD); % 坐標(biāo)轉(zhuǎn)換矩陣 EKT=T*EK*T； % 生成整體單剛(整體坐標(biāo)系)

5、 % 組成總剛 按3*3子塊加入總剛中（共計(jì)4塊) for j=1:2 %對行進(jìn)行循環(huán)-按結(jié)點(diǎn)號循環(huán) N1=LNODS(i,j)*3; % j結(jié)點(diǎn)第3個(gè)位移的整體編碼 for k=1:2 %對列進(jìn)行循環(huán)-按結(jié)點(diǎn)號循環(huán) N2=LNODS(i,k)*3; % k結(jié)點(diǎn)第3個(gè)位移的整體編碼 HK(N1-2):N1,(N2-2):N2)=HK(N1-2):N1,(N2-2):N2).+EKT(j*3-2:j*3,k*3-2:k*3); end endend% 由結(jié)點(diǎn)力與非結(jié)點(diǎn)力生成總荷載向量列陣for i=1:NFPOIN % 對結(jié)點(diǎn)荷載個(gè)數(shù)進(jìn)行循環(huán) N1=FPOIN(i,1); % 作用荷載的結(jié)點(diǎn)號

6、N1=N1*3-3; % 該結(jié)點(diǎn)號對應(yīng)第一個(gè)位移編碼 - 1 for j=1:3 FORCE(N1+j)=FORCE(N1+j)+FPOIN(i,j+1);%取結(jié)點(diǎn)荷載 end end% 計(jì)算由非結(jié)點(diǎn)荷載引起的等效結(jié)點(diǎn)荷載for i=1:NFPRES % 對非結(jié)點(diǎn)荷載個(gè)數(shù)進(jìn)行循環(huán) F0=ele_FPRES(i,FPRES,LNODS,COORD); %計(jì)算單元固端力% 對單元局部桿端力要進(jìn)行坐標(biāo)轉(zhuǎn)換 ele=FPRES(i,1); % 取荷載所在的單元號T=zbzh(ele,LNODS,COORD); % 坐標(biāo)轉(zhuǎn)換矩陣F0=T*F0; NL=LNODS(ele,1); NR=LNODS(ele

7、,2); %單元的左右結(jié)點(diǎn)號 % 將單元固端力變成等效結(jié)點(diǎn)荷載（注意固端力與等效結(jié)點(diǎn)荷載符號相反）FORCE(3*NL-2):3*NL)=FORCE(3*NL-2):3*NL)-F0(1:3); FORCE(3*NR-2):3*NR)=FORCE(3*NR-2):3*NR)-F0(4:6); end% 總剛、總荷載進(jìn)行邊界條件處理for j=1:NVFIX % 對約束個(gè)數(shù)進(jìn)行循環(huán) N1=FIXED(j);HK(1:3*NPOIN,N1)=0; HK(N1,1:3*NPOIN)=0; HK(N1,N1)=1; % 將零位移約束對應(yīng)的行、列變成零，主元變成1 FORCE(N1)=0;end%-DI

8、SP=HKFORCE % 方程求解，HK先求逆再與力向量左乘%-% 求結(jié)構(gòu)各個(gè)單元內(nèi)力EDISP=zeros(6,1); % 單元位移列向量清零for i=1:NELEM % 對單元個(gè)數(shù)進(jìn)行循環(huán) for j=1:2 %對桿端循環(huán) % i單元左右端結(jié)點(diǎn)號*3 = 該結(jié)點(diǎn)的最后一個(gè)位移編碼 N1=LNODS(i,j)*3; % 取一端的單元位移列向量 EDISP(3*j-2:3*j)=DISP(N1-2:N1); end % 生成局部單剛(局部坐標(biāo)) 右手坐標(biāo)系 EK=ele_EK(i,LNODS,COORD,YOUNG); T=zbzh(i,LNODS,COORD); % 坐標(biāo)轉(zhuǎn)換矩陣 FE=E

9、K*T*EDISP; %計(jì)算局部坐標(biāo)桿端力（由結(jié)點(diǎn)位移產(chǎn)生） for j=1:NFPRES if FPRES(j,1) = i %成立時(shí)，當(dāng)前單元上有非結(jié)點(diǎn)荷載 F0=ele_FPRES(j,FPRES,LNODS,COORD);%單元固端力 FE=FE+F0; % 考慮由非結(jié)點(diǎn)荷載引起的桿端力 end end FE % 打印桿端力end%- ele_FPRES.m %計(jì)算單元固端力函數(shù)（ 正方向：X向右 Y向上 M逆時(shí)針）% 入口參數(shù)：荷載序號，荷載信息，單元信息，結(jié)點(diǎn)坐標(biāo)% 出口參數(shù)：單元固端力左右兩端的軸力、剪力、彎矩function F0=ele_FPRES(iFPRES,FPRES,

10、LNODS,COORD) ele=FPRES(iFPRES,1); %取荷載所在的單元號 G=FPRES(iFPRES,3); %單元荷載大小C=FPRES(iFPRES,4); %單元荷載與左端距離 NL=LNODS(ele,1); NR=LNODS(ele,2); %單元的左右結(jié)點(diǎn)號 dx=COORD(NR,1)-COORD(NL,1); % x 坐標(biāo)差 dy=COORD(NR,2)-COORD(NL,2); % y 坐標(biāo)差 L=sqrt(dx2+dy2); %單元長度% 計(jì)算公式中一些常出現(xiàn)的項(xiàng) D=L-C; C1=C/L; C2=C1*C1; C3=C1*C2; B1=D/L; B2=

11、B1/L; F0=0;0;0;0;0;0; %單元固端力清零 switch FPRES(iFPRES,2) case 1 %均布荷載 F0(2)=-G*C*(2-2*C2+C3)/2.0; F0(3)=-G*C*C*(6-8*C1+3*C2)/12.0; F0(5)=-G*C-F0(2); F0(6)=G*C*C*C1*(4-3*C1)/12.0; case 2 %橫向集中力 F0(2)=-G*B1*B2*(L+2*C); F0(3)=-G*C*B1*B1; F0(5)=-G*C2*(L+2*D)/L; F0(6)=G*D*C2; case 3 %縱向集中力 F0(1)=-G*B1; F0(4

12、)=-G*C1;endreturnele_EK.m % 計(jì)算單元?jiǎng)偠染仃嚭瘮?shù) EK% 入口參數(shù)：單元號、單元信息數(shù)組、結(jié)點(diǎn)坐標(biāo)、彈性模量% 出口參數(shù)：局部單元?jiǎng)偠染仃嘐Kfunction EK=ele_EK(i,LNODS,COORD,E) NL=LNODS(i,1); NR=LNODS(i,2); %左右結(jié)點(diǎn)號 dx=COORD(NR,1)-COORD(NL,1); % x 坐標(biāo)差 dy=COORD(NR,2)-COORD(NL,2); % y 坐標(biāo)差 L=sqrt(dx2+dy2); %單元長度 A=LNODS(i,3); I=LNODS(i,4); %面積；慣性矩% 生成單剛(局部坐標(biāo))

13、 右手坐標(biāo)系EK =E*A/L 0 0 -E*A/L 0 0;. 0 12*E*I/L3 6*E*I/L2 0 -12*E*I/L3 6*E*I/L2;. 0 6*E*I/L2 4*E*I/L 0 -6*E*I/L2 2*E*I/L;. -E*A/L 0 0 E*A/L 0 0;. 0 -12*E*I/L3 -6*E*I/L2 0 12*E*I/L3 -6*E*I/L2;. 0 6*E*I/L2 2*E*I/L 0 -6*E*I/L2 4*E*I/L;return%-zbzh.m % 形成第i單元的坐標(biāo)轉(zhuǎn)換矩陣函數(shù) T(6,6)% 入口參數(shù)：單元號，單元信息，結(jié)點(diǎn)坐標(biāo)% 出口參數(shù)：坐標(biāo)轉(zhuǎn)換矩陣（整體向局部投影）function T=zbzh(i,LNODS,COORD) NL=LNODS(i,1); %左結(jié)點(diǎn)號 NR=LNODS(i,2); %右結(jié)點(diǎn)號 dx=COORD(NR,1)-COORD(NL,