2016基本平面剛架各種荷載MATLAB程序(共5頁)

1、精選優(yōu)質(zhì)文檔-傾情為你奉上% 平面剛架MATLAB程序% 2003.9.16 2007.2.28 2008.4.1 2009.10 2011.10 2013.9 2014.09 2016.03%*% 變量說明% NPOIN NELEM NVFIX NFPOIN NFPRES% 總結(jié)點數(shù),單元數(shù), 約束個數(shù), 受力結(jié)點數(shù), 非結(jié)點力數(shù)% COORD LNODS YOUNG% 結(jié)構(gòu)節(jié)點坐標數(shù)組, 單元定義數(shù)組, 彈性模量% FPOIN FPRES FORCE FIXED% 結(jié)點力數(shù)組，非結(jié)點力數(shù)組，總體荷載向量, 約束信息數(shù)組% HK DISP% 總體剛度矩陣，結(jié)點位移向量%*format sho

2、rt e %設定輸出類型clear %清除內(nèi)存變量FP1=fopen('6-6.txt','rt') %打開初始數(shù)據(jù)文件%讀入控制數(shù)據(jù) NELEM=fscanf(FP1,'%d',1); %單元數(shù)NPOIN=fscanf(FP1,'%d',1); %結(jié)點數(shù)NVFIX=fscanf(FP1,'%d',1); %約束數(shù)NFPOIN=fscanf(FP1,'%d',1); %作用荷載的結(jié)點個數(shù)NFPRES=fscanf(FP1,'%d',1); %非結(jié)點荷載數(shù)YOUNG=fscanf(FP

3、1,'%f',1); %彈性模量% 讀取結(jié)構(gòu)信息LNODS=fscanf(FP1,'%f',6,NELEM)' % 單元定義： 左、右結(jié)點號，面積，慣性矩，線膨脹系數(shù)，截面高度（共計 NELEM組）COORD=fscanf(FP1,'%f',2,NPOIN)' % 坐標： x,y坐標（共計 NPOIN 組）FPOIN=fscanf(FP1,'%f',4,NFPOIN)' % 節(jié)點力（共計 NFPOIN 組）：受力結(jié)點號、X方向力（向右正），% Y方向力（向上正），M力偶（逆時針正）FPRES=fscanf(

4、FP1,'%f',7,NFPRES)' % 均布力（共計% NFPRES 組）：單元號、荷載類型、荷載大小、距離左端長度，溫差=（下端-上端）梯形上邊。下邊(改)% 荷載類型 1-均布荷載 2-橫向集中力 3-縱向集中力 4-三角形荷載 5-溫度荷載 6-梯形荷載FIXED=fscan f(FP1,'%f',NVFIX)' % 約束信息：約束對應的位移編碼（共計 NVFIX 組）%-HK=zeros(3*NPOIN,3*NPOIN); % 張成總剛矩陣并清零FORCE=zeros(3*NPOIN,1); % 張成總荷載向量并清零%形成總剛for

5、i=1:NELEM % 對單元個數(shù)循環(huán) % 生成局部單剛(局部坐標) 右手坐標系 EK=ele_EK(i,LNODS,COORD,YOUNG); T=zbzh(i,LNODS,COORD); % 坐標轉(zhuǎn)換矩陣 EKT=T*EK*T； % 生成整體單剛(整體坐標系) % 組成總剛 按3*3子塊加入總剛中（共計4塊) for j=1:2 %對行進行循環(huán)-按結(jié)點號循環(huán) N1=LNODS(i,j)*3; % j結(jié)點第3個位移的整體編碼 for k=1:2 %對列進行循環(huán)-按結(jié)點號循環(huán) N2=LNODS(i,k)*3; % k結(jié)點第3個位移的整體編碼 HK(N1-2):N1,(N2-2):N2)=HK(

6、N1-2):N1,(N2-2):N2).+EKT(j*3-2:j*3,k*3-2:k*3); % 單剛3 x 3子塊疊加到總剛中 endend% 由結(jié)點力與非結(jié)點力生成總荷載向量列陣for i=1:NFPOIN % 對結(jié)點荷載個數(shù)進行循環(huán) N1=FPOIN(i,1); % 作用荷載的結(jié)點號 N1=N1*3-3; % 該結(jié)點號對應第一個位移編碼 - 1 for j=1:3 FORCE(N1+j)=FORCE(N1+j)+FPOIN(i,j+1);%取結(jié)點荷載 end end% 計算由非結(jié)點荷載引起的等效結(jié)點荷載for i=1:NFPRES % 對非結(jié)點荷載個數(shù)進行循環(huán) F0=ele_FPRES(

7、i,FPRES,LNODS,COORD,YOUNG); %計算單元固端力% 對單元局部桿端力要進行坐標轉(zhuǎn)換 ele=FPRES(i,1); % 取荷載所在的單元號T=zbzh(ele,LNODS,COORD); % 坐標轉(zhuǎn)換矩陣F0=T*F0; NL=LNODS(ele,1); NR=LNODS(ele,2); %單元的左右結(jié)點號 % 將單元固端力變成等效結(jié)點荷載（注意固端力與等效結(jié)點荷載符號相反）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); e

8、nd% 總剛、總荷載進行邊界條件處理for j=1:NVFIX % 對約束個數(shù)進行循環(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%-DISP=HKFORCE % 方程求解，HK先求逆再與力向量左乘% 求結(jié)構(gòu)各個單元內(nèi)力EDISP=zeros(6,1); % 單元位移列向量清零for i=1:NELEM % 對單元個數(shù)進行循環(huán) for j=1:2 %對桿端循環(huán) % i單元左右端結(jié)點號*3 = 該結(jié)點的最后一個位移編碼 N1=LNO

9、DS(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); % 坐標轉(zhuǎn)換矩陣 FE=EK*T*EDISP; %計算局部坐標桿端力（由結(jié)點位移產(chǎn)生） for j=1:NFPRES if FPRES(j,1) = i %成立時，當前單元上有非結(jié)點荷載 F0=ele_FPRES(j,FPRES,LNODS,COORD,YOUNG);%單元固端力 FE=FE+F0; % 考慮由非結(jié)點荷載引起

10、的桿端力 end end FE % 打印桿端力endend%- ele_FPRES.m %計算單元固端力函數(shù)（ 正方向：X向右 Y向上 M逆時針）% 入口參數(shù)：荷載序號，荷載信息，單元信息，結(jié)點坐標% 出口參數(shù)：單元固端力左右兩端的軸力、剪力、彎矩function F0=ele_FPRES(iFPRES,FPRES,LNODS,COORD, E) ele=FPRES(iFPRES,1); %取荷載所在的單元號 G=FPRES(iFPRES,3); %單元荷載大小C=FPRES(iFPRES,4); %單元荷載與左端距離W= FPRES(iFPRES,5); %單元下上端溫差S= FPRES(i

11、FPRES,6); %梯形長邊荷載X= FPRES(iFPRES,6); %梯形短邊荷載H= LNODS(ele,6); %單元截面高度P = LNODS(ele,5); %單元截面線膨脹系數(shù) NL=LNODS(ele,1); NR=LNODS(ele,2); %單元的左右結(jié)點號 dx=COORD(NR,1)-COORD(NL,1); % x 坐標差 dy=COORD(NR,2)-COORD(NL,2); % y 坐標差 L=sqrt(dx2+dy2); %單元長度% 計算公式中一些常出現(xiàn)的項 D=L-C; C1=C/L; C2=C1*C1; C3=C1*C2; B1=D/L; B2=B1/L

12、; 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)=-G

13、*C1;case 4 %三角形荷載 F0(2)=-7*G*L/20;F0(3)=G*L2/20; F0(5)=3*G*L/20; F0(6)=G*L*L/30;case 5 %溫度荷載 F0(3)=E*I*W*P/H; F0(6)=-E*I*W*P/H;case 6 %梯形荷載（改） F1(2)=-X*C*(2-2*C2+C3)/2.0; F1(3)=-X*C*C*(6-8*C1+3*C2)/12.0; F2(5)=-X*C-F0(2); F2(6)=X*C*C*C1*(4-3*C1)/12.0; F3(2)=-7*(S-G)*L/20;F3(3)=(S-G)*L2/20; F4(5)=3*(

14、S-G)*L/20; F4(6)=(S-G)*L*L/30;F0(2)= F1(2)+ F3(2);F0(3)= F1(3)+ F3(3);F0(5)= F2(5)+ F4(5);F0(6)= F2(6)+ F4(6); endreturnele_EK.m % 計算單元剛度矩陣函數(shù) EK% 入口參數(shù)：單元號、單元信息數(shù)組、結(jié)點坐標、彈性模量% 出口參數(shù)：局部單元剛度矩陣EKfunction EK=ele_EK(i,LNODS,COORD,E) NL=LNODS(i,1); NR=LNODS(i,2); %左右結(jié)點號 dx=COORD(NR,1)-COORD(NL,1); % x 坐標差 dy=

15、COORD(NR,2)-COORD(NL,2); % y 坐標差 L=sqrt(dx2+dy2); %單元長度 A=LNODS(i,3); I=LNODS(i,4); %面積；慣性矩% 生成單剛(局部坐標) 右手坐標系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單元的坐標轉(zhuǎn)換矩陣函數(shù) T(6,6)% 入口參數(shù)：單元號，單元信息，結(jié)點坐標% 出口參數(shù)：坐標轉(zhuǎn)換矩陣（整體向局部投影）function T=zbzh(i,LNODS,COORD) NL=LNODS(i,1); %左結(jié)點號 NR=LNODS(i,2); %右