計(jì)算力學(xué)課程設(shè)計(jì)

計(jì)算力學(xué)課程設(shè)計(jì)說明書班級姓名學(xué)號日期2015年12月25日計(jì)算力學(xué)課程設(shè)計(jì)桁架的有限元計(jì)算方法引言有限元分析是求取復(fù)雜微分方程近似解的一種非常有效的工具，是現(xiàn)代數(shù)字化科技的一種重要的基礎(chǔ)性原理。將它應(yīng)用于工程技術(shù)中，可成為工程設(shè)計(jì)和分析的可靠工具。而在桁架結(jié)構(gòu)中，運(yùn)用有限元的方法，通過現(xiàn)代有限元分析軟件如MATLAB，ANSYS等可輕易的求得各個(gè)桿件的內(nèi)力等。例如下圖的桁架結(jié)構(gòu)，運(yùn)用有限元法可十分清晰的了解各桿件的受力狀況。1基本力學(xué)模型如下圖1所示圖12有限元計(jì)算原理首先，明確Matlab程序要實(shí)現(xiàn)的5個(gè)重要模塊分別為：單元剛度矩陣的求解、單元組裝、節(jié)點(diǎn)位移的求解、單元應(yīng)力的求解、節(jié)點(diǎn)力的求解。下面給出這5個(gè)模塊的實(shí)現(xiàn)。(1)單元剛度矩陣求解定義函數(shù)Bar2D2Node_Stiffness，該函數(shù)計(jì)算單元的剛度矩陣，輸入彈性模量E，橫截面積A，兩個(gè)節(jié)點(diǎn)坐標(biāo)輸出單元剛度矩陣k(4X4)。具體代碼如下：functionk=Bar2D2Node_Stiffness(E,A,x1,y1,x2,y2)L=sqrt((x2-x1)*(x2-x1)+(y2-y1)*(y2-y1));x=acos((x2-x1)/L);C=cos(x);1安徽理工大學(xué)理學(xué)院計(jì)算力學(xué)課程設(shè)計(jì)S=sin(x);k=E*A/L*[C*CC*S-C*C-C*S;C*SS*S-C*S-S*S;-C*C-C*SC*CC*S;-C*S-S*SC*SS*S];(2)單元組裝定義函數(shù)Bar2D2Node_Assembly,該函數(shù)進(jìn)行單元剛度矩陣的組裝，輸入單元剛度矩陣k,單元的節(jié)點(diǎn)編號i、j。輸出整體剛度矩陣KK，具體代碼如下：functionz=Bar2D2Node_Assembly(KK,k,i,j)DOF(1)=2*i-1;DOF(2)=2*i;DOF(3)=2*j-1;DOF(4)=2*j;forn1=1:4forn2=1:4KK(DOF(n1),DOF(n2))=KK(DOF(n1),DOF(n2))+k(n1,n2);endendz=KK;(3)節(jié)點(diǎn)位移的求解定義函數(shù)Bar2D2Node_Disp(KK,num,p)，該函數(shù)輸入KK為總體剛度矩陣；num為活動自由度編號數(shù)組；p為活動自由度方向上的節(jié)點(diǎn)力；輸出節(jié)點(diǎn)位移列陣。具體代碼如下：functionu=Bar2D2Node_Disp(KK,num,p)k=KK(num,num);u=k\p;(4)單元應(yīng)力的求解定義函數(shù)函數(shù)Bar2D2Node_Stress(E,x1,y1,x2,y2,u)，該函數(shù)計(jì)算單元的應(yīng)力輸入彈性模量E，第一個(gè)節(jié)點(diǎn)坐標(biāo)(x1,y1)，第二個(gè)節(jié)點(diǎn)坐標(biāo)(x2,y2)單元節(jié)點(diǎn)位移矢量u，返回單元應(yīng)力標(biāo)量stress。具體代碼如下：functionstress=Bar2D2Node_Stress(E,x1,y1,x2,y2,u)L=sqrt((x2-x1)*(x2-x1)+(y2-y1)*(y2-y1));x=acos((x2-x1)/L);C=cos(x);S=sin(x);stress=E/L*[-C-SCS]*u;(5)計(jì)算節(jié)點(diǎn)力定義函數(shù)Bar2D2Node_Forces(KK,q)，該函數(shù)用于計(jì)算節(jié)點(diǎn)力，KK為剛度矩陣，q為節(jié)點(diǎn)位移陣列。functionP=Bar2D2Node_Forces(KK,q)2安徽理工大學(xué)理學(xué)院計(jì)算力學(xué)課程設(shè)計(jì)P=KK*q;以上5個(gè)函數(shù)寫入MATLAB的M文件中完成函數(shù)的定義，以備命令窗口的調(diào)用。至此，基于MATLAB的桿單元有限元分析的程序設(shè)計(jì)的準(zhǔn)備工作已經(jīng)完成。3 計(jì)算結(jié)果及分析針對上面的具體模型運(yùn)用MATLAB實(shí)現(xiàn)有限元的計(jì)算：（1）結(jié)構(gòu)的離散化與編號對該結(jié)構(gòu)進(jìn)行自然離散，節(jié)點(diǎn)編號和單元編號如上圖所示（2）計(jì)算各單元的剛度矩陣（基于國際標(biāo)準(zhǔn)單位）輸入彈性模量E、橫截面積A，各點(diǎn)坐標(biāo)。然后分別針對單元1，2,3和4,調(diào)用4次Bar2D2Node_Stiffness，就可以得到單元的剛度矩陣。對應(yīng)的主程序中代碼：E=2.95e11;A=0.0001;x1=0;y1=0;x2=0.4;y2=0;x3=0.4;y3=0.3;x4=0;y4=0.3;k1=Bar2D2Node_Stiffness(E,A,x1,y1,x2,y2)k2=Bar2D2Node_Stiffness(E,A,x2,y2,x3,y3)k3=Bar2D2Node_Stiffness(E,A,x1,y1,x3,y3)k4=Bar2D2Node_Stiffness(E,A,x4,y4,x3,y3)計(jì)算結(jié)果如下kl=737500000-73750000□a00□-73750000073750000□a00□k2=1.0e+0070.00000.0000-o.oooa-o.oooa0.00009.3333-O.OOOQ-5.S333-0.0000-0.0000o.oooao.oooa-0.0000-9.S3330.00009.8333圖2單元1與單元2的剛度矩陣3安徽理工大學(xué)理學(xué)院計(jì)算力學(xué)課程設(shè)計(jì)k3-1.Oe+OOr率3.77602.8320-3.7730-2.83202.S3202.1240-2.S320-2.1240-3.7760-2.83203.77502.8320-2.S320-2.12402.83202.1240k4=73F500000-^3750000a00□a-73F500000^3750000a00□a圖3單元3與單元4的剛度矩陣(3)建立整體剛度方程由于該結(jié)構(gòu)共有4個(gè)節(jié)點(diǎn)，因此，設(shè)置結(jié)構(gòu)總的剛度矩陣為KK(8X8),先對KK清零，然后四次調(diào)用函數(shù)Bar2D2Node.Assembly進(jìn)行剛度矩陣的組裝。相關(guān)主程序代碼為：KK=zeros(8,8);%由于該結(jié)構(gòu)共有4個(gè)節(jié)點(diǎn)，因此，結(jié)構(gòu)總的剛度矩陣為KK(8X8),先對K清零，然后四次調(diào)用函數(shù)Bar2D2Node.Assembly進(jìn)行剛度矩陣的組裝。KK=Bar2D2Node_Assembly(KK,k1,1,2);KK=Bar2D2Node_Assembly(KK,k2,2,3);KK=Bar2D2Node_Assembly(KK,k3,1,3);KK=Bar2D2Node_Assembly(KK,k4,4,3);結(jié)果如下KE二1.0e+008卡1.11510.2S32-0.73750-C.3775-0.2832000.20320.2124□a-0.2S32-0.212400-0.737500.7375o.oooa-o.oooa-o.ooao00000.00000.9S33-o.oooa-0.9S3300-0.3776-0.2332-0.0000-o.oooa1.11510.283273^50-0.2832-0.2124-0.0000-0.9S33C.2S321.135700000a-0.73^500.73^50000aa0004安徽理工大學(xué)理學(xué)院計(jì)算力學(xué)課程設(shè)計(jì)圖4整體剛度矩陣(4)邊界條件的處理及剛度方程的求解由圖可以看出，被約束的自由度有：節(jié)點(diǎn)1的x,y方向自由度，節(jié)點(diǎn)2的y方向自由度，4節(jié)點(diǎn)的x、y方向兩個(gè)自由度。則活動自由度編號為3,5,6.活動自由度對應(yīng)的節(jié)點(diǎn)載荷F3=20000N，F(xiàn)5=0N，F(xiàn)6=25000N，采用高斯消去法進(jìn)行求解，對應(yīng)的代碼為：num=[3,5,6];%可活動的自由度編號p=[20000;0;-25000];u=Bar2D2Node_Disp(KK,num,p)結(jié)果如下u=Cle-003*0.27120.0565-a.2225由此可知u2=0.2712*10"-3,u3=0.0565*10"-3,v3=-0.2225*10"-3o(5)支反力的計(jì)算在得到整個(gè)結(jié)構(gòu)的節(jié)點(diǎn)位移后，由原整體剛度方程就可以計(jì)算出對應(yīng)的支反力。這部分對應(yīng)的主程序的代碼如下：q=zeros(8,1);q(num)=u;%節(jié)點(diǎn)位移陣列P=Bar2D2Node_Forces(KK,q)結(jié)果如下5安徽理工大學(xué)理學(xué)院計(jì)算力學(xué)課程設(shè)計(jì)F=1.Oe+004*-1.50330.31252.00001SF50.0000-2.5000-0.4167由此可知支反力為FRX1=-15833N,FRY1=3125N,FRY2=21879N,FRX4=-4167N,FRY4=0N。(6)單元應(yīng)力的計(jì)算先從整體位移列陣q中提取出單元的位移列陣，然后，調(diào)用計(jì)算單元應(yīng)力的函數(shù)Bar2D2Node_Stress,就可以得到各個(gè)單元的應(yīng)力分量。u1=[q(1);q(2);q(3);q(4)]stress1=Bar2D2Node_Stress(E,x1,y1,x2,y2,u1)u2=[q(3);q(4);q(5);q(6)]stress2=Bar2D2Node_Stress(E,x2,y2,x3,y3,u2)u3=[q(1);q(2);q(5);q(6)]stress3=Bar2D2Node_Stress(E,x1,y1,x3,y3,u3)u4=[q(7);q(8);q(5);q(6)]stress4=Bar2D2Node_Stress(E,x4,y4,x3,y3,u4)結(jié)果如下：stress1=200000000,stress2=-2.1875e+008,stress3=-5.2083e+007，stress4=4.1667e+007參考文獻(xiàn)[1]曾攀.有限元基礎(chǔ)教程.北京：高等教育出版社2009,7[2]王勖成.有限單元法.北京:清華大學(xué)出版社,2003,7.附錄(程序)子程序1：functionk=Bar2D2Node_Stiffness(E,A,x1,y1,x2,y2)L=sqrt((x2-x1)*(x2-x1)+(y2-y1)*(y2-y1));6安徽理工大學(xué)理學(xué)院計(jì)算力學(xué)課程設(shè)計(jì)x=acos((x2-x1)/L);C=cos(x);S=sin(x);k=E*A/L*[C*CC*S-C*C-C*S;C*SS*S-C*S-S*S;-C*C-C*SC*CC*S;-C*S-S*SC*SS*S];子程序2：functionz=Bar2D2Node_Assembly(KK,k,i,j)DOF(1)=2*i-1;DOF(2)=2*i;DOF(3)=2*j-1;DOF(4)=2*j;forn1=1:4forn2=1:4KK(DOF(n1),DOF(n2))=KK(DOF(n1),DOF(n2))+k(n1,n2);endendz=KK;子程序3：functionu=Bar2D2Node_Disp(KK,num,p)k=KK(num,num);u=k\p;子程序4：functionstress=Bar2D2Node_Stress(E,x1,y1,x2,y2,u)L=sqrt((x2-x1)*(x2-x1)+(y2-y1)*(y2-y1));x=acos((x2-x1)/L);C=cos(x);S=sin(x);stress=E/L*[-C-SCS]*u;子程序5：functionP=Bar2D2Node_Forces(KK,q)P=KK*q;主程序：E=2.95e11;A=0.0001;x1=0;y1=0;x2=0.4;y2=0;x3=0.4;y3=0.3;7安徽理工大學(xué)理學(xué)院計(jì)算力學(xué)課程設(shè)計(jì)x4=0;y4=0.3;k1=Bar2D2Node_Stiffness(E,A,x1,y1,x2,y2)k2=Bar2D2Node_Stiffness(E,A,x2,y2,x3,y3)k3=Bar2D2Node_Stiffness(E,A,x1,y1,x3,y3)k4=Bar2D2Node_Stiffness(E,A,x4,y4,x3,y3)KK=zeros(8,8);KK=Bar2D2Node_Assembly(KK,k1,1,2);KK=Bar2D2Node_Assembly(KK,k2,2,3);KK=Bar2D2Node_Assembly(KK,k3,1,3);KK=Bar2D2Node_Assembly(KK,k4,4,3)num=[3,5,6];p=[20000;0;-25000];u=Bar2D2Node_Disp(KK,num,p)q=zeros(8,1);q(num)=u;P=Bar2D2Node_Forces(KK,q)u1=[q(1);