《基于MATLAB的平面剛架靜力分析.doc》由會員分享，可在線閱讀，更多相關(guān)《基于MATLAB的平面剛架靜力分析.doc（17頁珍藏版）》請在裝配圖網(wǎng)上搜索。

基于MATLAB的平面剛架靜力分析 為了進一步理解有限元方法計算的過程，本文根據(jù)矩陣位移法的基本原理應(yīng)用MATLAB編制計算程序?qū)σ云矫鎰偧芙Y(jié)構(gòu)進行了靜力分析。本文還利用ANSYS大型商用有限元分析軟件對矩陣位移法的計算結(jié)果進行校核，發(fā)現(xiàn)兩者計算結(jié)果相當(dāng)吻合，驗證了計算結(jié)果的可靠性。 一、 問題描述 如圖1所示的平面剛架，各桿件的材料及截面均相同，E=210GPa，截面為0.120.2m的實心矩形，現(xiàn)要求解荷載作用下剛架的位移和內(nèi)力。 圖1 二、矩陣位移法計算程序編制 為編制程序方便考慮，本文計算中采用“先處理法”。具體的計算步驟如下。 （1） 對結(jié)構(gòu)進行離散化，對結(jié)點和單元進行編號，建立結(jié)構(gòu)（整體）坐標(biāo)系和單元（局部）坐標(biāo)系，并對結(jié)點位移進行編號； （2） 對結(jié)點位移分量進行編碼，形成單元定位向量； （3） 建立按結(jié)構(gòu)整體編碼順序排列的結(jié)點位移列向量，計算固端力、等效結(jié)點荷載 及綜合結(jié)點荷載列向量； （4） 計算個單元局部坐標(biāo)系的剛度矩陣，通過坐標(biāo)變換矩陣形成整體坐標(biāo)系下的單元剛度矩陣 ； （5） 利用單元定位向量形成結(jié)構(gòu)剛度矩陣； （6） 按式 求解未知結(jié)點位移； （7） 計算各單元的桿端力。 根據(jù)上述步驟編制了平面剛架的分析程序。程序中單元剛度矩陣按下式計算。 轉(zhuǎn)換矩陣則按下式計算。 計算程序框圖如圖2所示，具體的程序代碼見附錄1。 圖2 MATLAB矩陣分析法程序框圖 三、解題步驟 取整體坐標(biāo)系如圖3所示，對結(jié)構(gòu)進行離散化，對結(jié)點和單元進行編號如圖4所示，局部坐標(biāo)系用單元中箭頭的方向表示，原始數(shù)據(jù)如下： 圖3 圖4 剛架結(jié)點輸入矩陣為， x=[0 0;0 -5;1.63 -6.37;4 -5;4 -1;4 2]; 各單元定位向量為， locvec1=[1 2 3 0 0 0]; locvec2=[1 2 3 4 5 6]; locvec3=[4 5 6 7 8 9]; locvec4=[7 8 9 10 11 12]; locvec5=[10 11 12 0 0 0]; 輸入截面參數(shù)， E=2.1e11;%E=210GPa a=0.12; %矩形截面長0.12m b=0.2; %矩形截面寬0.2m 輸入整體坐標(biāo)系下各單元結(jié)點荷載列陣， f(1,:)=zeros(1,6); f(2,:)=[0 0 0 0 40e3 0]; f(3,:)=zeros(1,6); f(4,:)=[0 0 0 -50e3 0 0]; f(5,:)=zeros(1,6); 輸入整體坐標(biāo)系下單元1等效節(jié)點荷載 q=10e3; %10kN/m fe=[0.5*q*l(1),0,-q*l(1)^2/12,0.5*q*l(1),0,q*l(1)^2/12]; 由此計算得到平面剛架整體坐標(biāo)系下的結(jié)點位移(m), d= 0.0035 0.0000 -0.0004 0.0030 -0.0005 -0.0004 0.0027 0.0000 0.0016 -0.0051 0.0000 -0.0006 各個單元的桿端力如表1所示， 表1 各單元桿端力 　 單元 1 2 3 4 5 i端 Fx(kN) -17917.047 17917.05 17917.05 17917.05 -32083 Fy(kN) 17507.3731 -17507.4 22492.63 22492.63 22492.63 Mz(kNm) 1897.83076 -1897.83 2092.833 -26668.3 44999.85 j端 Fx(kN) -32082.953 -17917 -17917 32082.95 32082.95 Fy(kN) -17507.373 -22492.6 -22492.6 -22492.6 -22492.6 Mz(kNm) -37312.596 -2092.83 26668.34 -44999.8 51249.01 四、計算結(jié)果校核 在ANSYS中使用beam3單元，按照如圖4所示的離散結(jié)構(gòu)建立平面剛架模型施加約束和荷載，得到的有限元模型如圖5所示。計算分析后得到結(jié)構(gòu)的軸力圖、剪力圖和彎矩圖如圖6、7、8所示，命令流見附錄2。 圖5 有限元模型 圖6 軸力圖(kN) 圖7 剪力圖（kN） 圖8 彎矩圖（kNm） 從ANSYS計算結(jié)果中提取各結(jié)點位移、內(nèi)力，并與矩陣位移法分析的結(jié)果比較，得到表2、3。 表2 各節(jié)點位移比較 節(jié)點號 項目 矩陣位移法 ANSYS 誤差 1 Ux(m) 0 0 0 Uy(m) 0 0 0 Rotz(rad) 0 0 0 2 Ux(m) 0.003478 0.00348 -2E-06 Uy(m) 0.0000174 0.0000174 0 Rotz(rad) -0.00037 -0.000365 -5E-06 3 Ux(m) 0.003018 0.00302 -2E-06 Uy(m) -0.00051 -0.000512 2E-06 Rotz(rad) -0.00038 -0.000378 -2E-06 4 Ux(m) 0.002687 0.00269 -3E-06 Uy(m) 0.0000312 0.0000312 0 Rotz(rad) 0.001624 0.00162 4E-06 5 Ux(m) -0.00513 -0.005145 1.5E-05 Uy(m) 0.0000134 0.0000134 0 Rotz(rad) -0.00056 -0.000557 -3E-06 6 Ux(m) 0 0 0 Uy(m) 0 0 0 Rotz(rad) 0 0 0 表3 各結(jié)點內(nèi)力比較 節(jié)點號 項目 矩陣位移法 ANSYS 誤差 1 Fx(kN) -32.083 -32.080 -0.003 Fy(kN) -17.507 -17.503 -0.004 Mz(kNm) -37.313 -37.307 -0.006 2 Fx(kN) -17.917 -17.920 0.003 Fy(kN) 17.507 17.503 0.004 Mz(kNm) 1.898 1.909 -0.011 3 Fx(kN) -17.917 -17.920 0.003 Fy(kN) -22.493 -22.497 0.004 Mz(kNm) -2.093 -2.110 0.018 4 Fx(kN) -17.917 -17.920 0.003 Fy(kN) -22.493 -22.497 0.004 Mz(kNm) 26.668 26.682 -0.014 5 Fx(kN) -32.083 -32.080 -0.003 Fy(kN) -22.493 -22.497 0.004 Mz(kNm) -45.000 -44.999 -0.001 6 Fx(kN) 32.083 32.080 0.003 Fy(kN) -22.493 -22.497 0.004 Mz(kNm) 51.249 51.239 0.010 由表2、表3的結(jié)果對比可知，兩種方法的計算結(jié)果十分接近，誤差均可以忽略不計，從而驗證了計算結(jié)果的可靠性與準(zhǔn)確性。  四、結(jié)論 通過對一個平面剛架靜力分析問題的求解，本文編制的平面剛架靜力分析程序計算結(jié)果與有限元軟件ANSYS計算結(jié)果校核后，發(fā)現(xiàn)兩者計算結(jié)果十分接近，誤差可忽略不計，說明該程序計算結(jié)果的可靠性與精確性。 附錄1 矩陣位移法計算程序 pmgj.m 計算主程序 clear clc %結(jié)點坐標(biāo) x=[0 0;0 -5;1.63 -6.37;4 -5;4 -1;4 2]; % x1=0;y1=0; % x2=0;y2=-5; % x3=1.63;y3=-6.37; % x4=4;y4=-5; % x5=4;y5=-1; % x6=4;y6=2; %單元定位向量 locvec1=[1 2 3 0 0 0]; locvec2=[1 2 3 4 5 6]; locvec3=[4 5 6 7 8 9]; locvec4=[7 8 9 10 11 12]; locvec5=[10 11 12 0 0 0]; %剛架的材料特性 截面特性 E=2.1e11;%E=210GPa a=0.12; %矩形截面長0.12m b=0.2; %矩形截面寬0.2m A=a*b; I=(a*b^3)/12; clear a b % 單元長度計算 for i=1:5 dx=x(i,1)-x(i+1,1); dy=x(i,2)-x(i+1,2); l(i)=(dx^2+dy^2)^0.5; end %生成轉(zhuǎn)換矩陣 t1=coortrans(x(2,1),x(2,2),x(1,1),x(1,2)); t2=coortrans(x(2,1),x(2,2),x(3,1),x(3,2)); t3=coortrans(x(3,1),x(3,2),x(4,1),x(4,2)); t4=coortrans(x(4,1),x(4,2),x(5,1),x(5,2)); t5=coortrans(x(5,1),x(5,2),x(6,1),x(6,2)); %結(jié)點荷載（結(jié)構(gòu)坐標(biāo)系下） f(1,:)=zeros(1,6); f(2,:)=[0 0 0 0 40e3 0]; f(3,:)=zeros(1,6); f(4,:)=[0 0 0 -50e3 0 0]; f(5,:)=zeros(1,6); %單元等效節(jié)點荷載（結(jié)構(gòu)坐標(biāo)系下） q=10e3; %10kN/m fe=[0.5*q*l(1),0,-q*l(1)^2/12,0.5*q*l(1),0,q*l(1)^2/12]; %單元坐標(biāo)系下的值 fe1=[0,0.5*q*l(1),q*l(1)^2/12,0,0.5*q*l(1),-q*l(1)^2/12]; %生成單元剛度矩陣(單元坐標(biāo)系） k1=elemstiffm(E,A,l(1),I); k2=elemstiffm(E,A,l(2),I); k3=elemstiffm(E,A,l(3),I); k4=elemstiffm(E,A,l(4),I); k5=elemstiffm(E,A,l(5),I); %將單元坐標(biāo)系下的單元剛度矩陣轉(zhuǎn)換到結(jié)構(gòu)坐標(biāo)系下 K1=t1*k1*t1; K2=t2*k2*t2; K3=t3*k3*t3; K4=t4*k4*t4; K5=t5*k5*t5; %組裝總體剛度矩陣 K=zeros(12); F=zeros(12,1); NonLog=(locvec1~=0); ij=find(NonLog); IJ=locvec1(NonLog); K(IJ,IJ)=K(IJ,IJ)+K1(ij,ij); F(IJ)=F(IJ)+f(1,ij)+fe(ij); clear NonLog ij IJ NonLog=(locvec2~=0); ij=find(NonLog); IJ=locvec2(NonLog); K(IJ,IJ)=K(IJ,IJ)+K2(ij,ij); F(IJ)=F(IJ)+f(2,ij); clear NonLog ij IJ NonLog=(locvec3~=0); ij=find(NonLog); IJ=locvec3(NonLog); K(IJ,IJ)=K(IJ,IJ)+K3(ij,ij); F(IJ)=F(IJ)+f(3,ij); clear NonLog ij IJ NonLog=(locvec4~=0); ij=find(NonLog); IJ=locvec4(NonLog); K(IJ,IJ)=K(IJ,IJ)+K4(ij,ij); F(IJ)=F(IJ)+f(4,ij); clear NonLog ij IJ NonLog=(locvec5~=0); ij=find(NonLog); IJ=locvec5(NonLog); K(IJ,IJ)=K(IJ,IJ)+K5(ij,ij); F(IJ)=F(IJ)+f(5,ij); clear NonLog ij IJ %節(jié)點位移 d=K\F; %單元1桿端力（結(jié)構(gòu)坐標(biāo)系下） d1=zeros(6,1); NonLog=(locvec1~=0); ij=find(NonLog); ij1=find(~NonLog); IJ=locvec1(NonLog); d1=d(IJ); d1(ij1)=0; F1=K1*d1-fe; %單元2桿端力（結(jié)構(gòu)坐標(biāo)系下） d2=zeros(6,1); NonLog=(locvec2~=0); ij=find(NonLog); ij1=find(~NonLog); IJ=locvec2(NonLog); d2=d(IJ); d2(ij1)=0; F2=K2*d2-f(2,:); %單元3桿端力（結(jié)構(gòu)坐標(biāo)系下） d3=zeros(6,1); NonLog=(locvec3~=0); ij=find(NonLog); ij1=find(~NonLog); IJ=locvec3(NonLog); d3=d(IJ); d3(ij1)=0; F3=K3*d3-f(3,:); %單元4桿端力（結(jié)構(gòu)坐標(biāo)系下） d4=zeros(6,1); NonLog=(locvec4~=0); ij=find(NonLog); ij1=find(~NonLog); IJ=locvec4(NonLog); d4=d(IJ); d4(ij1)=0; F4=K4*d4-f(4,:); %單元5桿端力（結(jié)構(gòu)坐標(biāo)系下） d5=zeros(6,1); NonLog=(locvec5~=0); ij=find(NonLog); ij1=find(~NonLog); IJ=locvec5(NonLog); d5=d(IJ); d5(ij1)=0; F5=K5*d5-f(5,:); coortrans.m 轉(zhuǎn)換矩陣生成函數(shù) % 轉(zhuǎn)換矩陣生成函數(shù)（單元坐標(biāo)系->結(jié)構(gòu)坐標(biāo)系） %y=coortrans(x1,y1,x2,y2),x1,y1,x2,y2為單元兩端結(jié)點在結(jié)構(gòu)坐標(biāo)系下的坐標(biāo)，單位都為m function y=coortrans(x1,y1,x2,y2) a=atan((y2-y1)/(x2-x1)); c=cos(a); s=sin(a); t=zeros(6); t(1,1)=c;t(1,2)=s; t(2,1)=-s;t(2,2)=c; t(3,3)=1; t(4,4)=c;t(4,5)=s; t(5,4)=-s;t(5,5)=c; t(6,6)=1; y=t; end elemstiffm.m 單元剛度矩陣生成函數(shù) % 單元剛度矩陣生成函數(shù) % y=elemstiffm(E,A,l,I),E,A,l,I均采用國際單位制 Pa m2 m m4 function y=elemstiffm(E,A,l,I) i1=E*A/l; i2=12*E*I/(l^3); i3=6*E*I/(l^2); i4=4*E*I/l; i5=2*E*I/l; k=zeros(6); k(1,1)=i1;k(1,4)=-i1; k(2,2)=i2;k(2,3)=i3;k(2,5)=-i2;k(2,6)=i3; k(3,3)=i4;k(3,5)=-i3;k(3,6)=i5; k(4,4)=i1; k(5,5)=i2;k(5,6)=-i3; k(6,6)=i4; k(3,2)=k(2,3); k(4,1)=k(1,4); k(5,2)=k(2,5);k(5,3)=k(3,5); k(6,2)=k(2,6);k(6,3)=k(3,6);k(6,5)=k(5,6); y=k; end 附錄2 ANSYS建模計算命令流 finish /clear /prep7 B=0.120$H=0.200$E=210000000 et,1,beam3 mp,ex,1,E mp,prxy,1,0.3 r,1,B*H,B*H*H*H/12,H k,1 k,2,0,5 k,3,1.6304,6.3681 k,4,4,5 k,5,4,1 k,6,4,-2 *set,i *do,i,1,5 l,i,i+1 *enddo lesize,all,0.5 latt,1,1,1 lmesh,all dk,1,all dk,6,all fk,3,fy,-40 fk,5,fx,-50 lsel,s,,,1 esll,s sfbeam,all,1,pres,10 allsel,all dtran ftran sftran /pbc,all,,2 /psf,pres,norm,2,0,1 eplot /solu solve /post1 /pbc,u,,1 !顯示支座約束符號，并圖形顯示變形 pldisp,1 !將當(dāng)前主要結(jié)果列表顯示 presol,elem !/pnum,sval,1 etable,mi,smisc,6 etable,mj,smisc,12 plls,mi,mj,-1 !彎矩圖 kN.m etable,qi,smisc,2 etable,qj,smisc,8 plls,qi,qj,-1 !剪力圖 kN etable,ni,smisc,1 etable,nj,smisc,7 plls,ni,nj,1 !軸力圖 kN