基本平面剛架MATLAB程序

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

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

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

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

5、 % 組成總剛 按3*3子塊加入總剛中（共計4塊) for j=1:2 %對行進行循環(huán)-按結點號循環(huán) N1=LNODS(i,j)*3; % j結點第3個位移的整體編碼 for k=1:2 %對列進行循環(huán)-按結點號循環(huán) N2=LNODS(i,k)*3; % k結點第3個位移的整體編碼 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% 由結點力與非結點力生成總荷載向量列陣for i=1:NFPOIN % 對結點荷載個數進行循環(huán) N1=FPOIN(i,1); % 作用荷載的結點號

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

7、,2); %單元的左右結點號 % 將單元固端力變成等效結點荷載（注意固端力與等效結點荷載符號相反）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% 總剛、總荷載進行邊界條件處理for j=1:NVFIX % 對約束個數進行循環(huán) N1=FIXED(j);HK(1:3*NPOIN,N1)=0; HK(N1,1:3*NPOIN)=0; HK(N1,N1)=1; % 將零位移約束對應的行、列變成零，主元變成1 FORCE(N1)=0;end%-DI

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

9、K*T*EDISP; %計算局部坐標桿端力（由結點位移產生） for j=1:NFPRES if FPRES(j,1) = i %成立時，當前單元上有非結點荷載 F0=ele_FPRES(j,FPRES,LNODS,COORD);%單元固端力 FE=FE+F0; % 考慮由非結點荷載引起的桿端力 end end FE % 打印桿端力end%- ele_FPRES.m %計算單元固端力函數（ 正方向：X向右 Y向上 M逆時針）% 入口參數：荷載序號，荷載信息，單元信息，結點坐標% 出口參數：單元固端力左右兩端的軸力、剪力、彎矩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); %單元的左右結點號 dx=COORD(NR,1)-COORD(NL,1); % x 坐標差 dy=COORD(NR,2)-COORD(NL,2); % y 坐標差 L=sqrt(dx2+dy2); %單元長度% 計算公式中一些常出現的項 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 % 計算單元剛度矩陣函數 EK% 入口參數：單元號、單元信息數組、結點坐標、彈性模量% 出口參數：局部單元剛度矩陣EKfunction EK=ele_EK(i,LNODS,COORD,E) NL=LNODS(i,1); NR=LNODS(i,2); %左右結點號 dx=COORD(NR,1)-COORD(NL,1); % x 坐標差 dy=COORD(NR,2)-COORD(NL,2); % y 坐標差 L=sqrt(dx2+dy2); %單元長度 A=LNODS(i,3); I=LNODS(i,4); %面積；慣性矩% 生成單剛(局部坐標)

13、 右手坐標系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單元的坐標轉換矩陣函數 T(6,6)% 入口參數：單元號，單元信息，結點坐標% 出口參數：坐標轉換矩陣（整體向局部投影）function T=zbzh(i,LNODS,COORD) NL=LNODS(i,1); %左結點號 NR=LNODS(i,2); %右結點號 dx=COORD(NR,1)-COORD(NL,