基于MATLAB的平面剛架靜力分析

1、工程計(jì)算實(shí)踐作業(yè)基于MATLAB的平面剛架靜力分析 為了進(jìn)一步理解有限元方法計(jì)算的過程，本文根據(jù)矩陣位移法的基本原理應(yīng)用MATLAB編制計(jì)算程序?qū)σ云矫鎰偧芙Y(jié)構(gòu)進(jìn)行了靜力分析。本文還利用ANSYS大型商用有限元分析軟件對(duì)矩陣位移法的計(jì)算結(jié)果進(jìn)行校核，發(fā)現(xiàn)兩者計(jì)算結(jié)果相當(dāng)吻合，驗(yàn)證了計(jì)算結(jié)果的可靠性。一、 問題描述如圖1所示的平面剛架，各桿件的材料及截面均相同，E=210GPa，截面為0.12×0.2m的實(shí)心矩形，現(xiàn)要求解荷載作用下剛架的位移和內(nèi)力。圖1二、矩陣位移法計(jì)算程序編制為編制程序方便考慮，本文計(jì)算中采用“先處理法”。具體的計(jì)算步驟如下。（1） 對(duì)結(jié)構(gòu)進(jìn)行離散化，對(duì)結(jié)點(diǎn)和單元進(jìn)

2、行編號(hào)，建立結(jié)構(gòu)（整體）坐標(biāo)系和單元（局部）坐標(biāo)系，并對(duì)結(jié)點(diǎn)位移進(jìn)行編號(hào)；（2） 對(duì)結(jié)點(diǎn)位移分量進(jìn)行編碼，形成單元定位向量；（3） 建立按結(jié)構(gòu)整體編碼順序排列的結(jié)點(diǎn)位移列向量，計(jì)算固端力、等效結(jié)點(diǎn)荷載 及綜合結(jié)點(diǎn)荷載列向量；（4） 計(jì)算個(gè)單元局部坐標(biāo)系的剛度矩陣，通過坐標(biāo)變換矩陣形成整體坐標(biāo)系下的單元?jiǎng)偠染仃?；（5） 利用單元定位向量形成結(jié)構(gòu)剛度矩陣；（6） 按式 求解未知結(jié)點(diǎn)位移；（7） 計(jì)算各單元的桿端力。根據(jù)上述步驟編制了平面剛架的分析程序。程序中單元?jiǎng)偠染仃嚢聪率接?jì)算。轉(zhuǎn)換矩陣則按下式計(jì)算。計(jì)算程序框圖如圖2所示，具體的程序代碼見附錄1。圖2 MATLAB矩陣分析法程序框圖三、解題步

3、驟取整體坐標(biāo)系如圖3所示，對(duì)結(jié)構(gòu)進(jìn)行離散化，對(duì)結(jié)點(diǎn)和單元進(jìn)行編號(hào)如圖4所示，局部坐標(biāo)系用單元中箭頭的方向表示，原始數(shù)據(jù)如下： 圖3 圖4剛架結(jié)點(diǎn)輸入矩陣為，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=210GPaa=0.12; %矩形截面長(zhǎng)0.12mb=0.2; %矩形截面寬0.2m輸入整體坐標(biāo)系下各單

4、元結(jié)點(diǎn)荷載列陣，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é)點(diǎn)荷載q=10e3; %10kN/mfe=0.5*q*l(1),0,-q*l(1)2/12,0.5*q*l(1),0,q*l(1)2/12;由此計(jì)算得到平面剛架整體坐標(biāo)系下的結(jié)點(diǎn)位移(m),d= 0.0035 0.0000 -0.0004 0.0030 -0.0005 -0.0004 0.0027 0.0000 0.0016 -0.0051 0.0000

5、-0.0006各個(gè)單元的桿端力如表1所示，表1 各單元桿端力單元12345i端Fx(kN)-17917.04717917.0517917.0517917.05-32083Fy(kN)17507.3731-17507.422492.6322492.6322492.63Mz(kN·m)1897.83076-1897.832092.833-26668.344999.85j端Fx(kN)-32082.953-17917-1791732082.9532082.95Fy(kN)-17507.373-22492.6-22492.6-22492.6-22492.6Mz(kN·m)-3731

6、2.596-2092.8326668.34-44999.851249.01四、計(jì)算結(jié)果校核在ANSYS中使用beam3單元，按照如圖4所示的離散結(jié)構(gòu)建立平面剛架模型施加約束和荷載，得到的有限元模型如圖5所示。計(jì)算分析后得到結(jié)構(gòu)的軸力圖、剪力圖和彎矩圖如圖6、7、8所示，命令流見附錄2。圖5 有限元模型圖6 軸力圖(kN)圖7 剪力圖（kN）圖8 彎矩圖（kN·m）從ANSYS計(jì)算結(jié)果中提取各結(jié)點(diǎn)位移、內(nèi)力，并與矩陣位移法分析的結(jié)果比較，得到表2、3。表2 各節(jié)點(diǎn)位移比較節(jié)點(diǎn)號(hào)項(xiàng)目矩陣位移法ANSYS誤差1Ux(m)000Uy(m)000Rotz(rad)0002Ux(m)0.0034

7、780.00348-2E-06Uy(m)0.00001740.00001740Rotz(rad)-0.00037-0.000365-5E-063Ux(m)0.0030180.00302-2E-06Uy(m)-0.00051-0.0005122E-06Rotz(rad)-0.00038-0.000378-2E-064Ux(m)0.0026870.00269-3E-06Uy(m)0.00003120.00003120Rotz(rad)0.0016240.001624E-065Ux(m)-0.00513-0.0051451.5E-05Uy(m)0.00001340.00001340Rotz(rad)

8、-0.00056-0.000557-3E-066Ux(m)000Uy(m)000Rotz(rad)000表3 各結(jié)點(diǎn)內(nèi)力比較節(jié)點(diǎn)號(hào)項(xiàng)目矩陣位移法ANSYS誤差1Fx(kN)-32.083-32.080-0.003Fy(kN)-17.507-17.503-0.004Mz(kN·m)-37.313-37.307-0.0062Fx(kN)-17.917-17.9200.003Fy(kN)17.50717.5030.004Mz(kN·m)1.8981.909-0.0113Fx(kN)-17.917-17.9200.003Fy(kN)-22.493-22.4970.004Mz(kN&

9、#183;m)-2.093-2.1100.0184Fx(kN)-17.917-17.9200.003Fy(kN)-22.493-22.4970.004Mz(kN·m)26.66826.682-0.0145Fx(kN)-32.083-32.080-0.003Fy(kN)-22.493-22.4970.004Mz(kN·m)-45.000-44.999-0.0016Fx(kN)32.08332.0800.003Fy(kN)-22.493-22.4970.004Mz(kN·m)51.24951.2390.010由表2、表3的結(jié)果對(duì)比可知，兩種方法的計(jì)算結(jié)果十分接近，誤差

10、均可以忽略不計(jì)，從而驗(yàn)證了計(jì)算結(jié)果的可靠性與準(zhǔn)確性。四、結(jié)論通過對(duì)一個(gè)平面剛架靜力分析問題的求解，本文編制的平面剛架靜力分析程序計(jì)算結(jié)果與有限元軟件ANSYS計(jì)算結(jié)果校核后，發(fā)現(xiàn)兩者計(jì)算結(jié)果十分接近，誤差可忽略不計(jì)，說明該程序計(jì)算結(jié)果的可靠性與精確性。附錄1 矩陣位移法計(jì)算程序pmgj.m 計(jì)算主程序clear clc%結(jié)點(diǎn)坐標(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

11、 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=210GPaa=0.12; %矩形截面長(zhǎng)0.12mb=0.2; %矩形截面寬0.2mA=a*b;I=(a*b3)/12;clear a b% 單元長(zhǎng)度計(jì)算for i=1:5 dx=x(i,1)-x(i+1,1); dy=x(i,2)-x(i+1,2); l(i)=(dx2+dy2)0.5;end%生成轉(zhuǎn)換矩陣t1=coortrans(x(2,1),x(2

12、,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é)點(diǎn)荷載（結(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é)

13、點(diǎn)荷載（結(jié)構(gòu)坐標(biāo)系下）q=10e3; %10kN/mfe=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;%生成單元?jiǎng)偠染仃?單元坐標(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)系下的單元?jiǎng)偠染?/p>

14、陣轉(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 IJNonLog=(locvec2=0);ij=find(NonLog)

15、;IJ=locvec2(NonLog);K(IJ,IJ)=K(IJ,IJ)+K2(ij,ij);F(IJ)=F(IJ)+f(2,ij)'clear NonLog ij IJNonLog=(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 IJNonLog=(locvec4=0);ij=find(NonLog);IJ=locvec4(NonLog);K(IJ,IJ)=K(IJ,IJ)+K4(ij,ij);F(IJ

16、)=F(IJ)+f(4,ij)'clear NonLog ij IJNonLog=(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é)點(diǎn)位移d=KF;%單元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-

17、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)系下）

18、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)系-

19、>結(jié)構(gòu)坐標(biāo)系）%y=coortrans(x1,y1,x2,y2),x1,y1,x2,y2為單元兩端結(jié)點(diǎn)在結(jié)構(gòu)坐標(biāo)系下的坐標(biāo)，單位都為mfunction 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;endelemstiffm.m 單元?jiǎng)偠染仃嚿珊瘮?shù)% 單元?jiǎng)偠染仃嚿珊瘮?shù)% y=elemstif

20、fm(E,A,l,I),E,A,l,I均采用國(guó)際單位制 Pa m2 m m4function y=elemstiffm(E,A,l,I)i1=E*A/l;i2=12*E*I/(l3);i3=6*E*I/(l2);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)=