有限元鋼架結(jié)構(gòu)分析~手算+matlab+ansys模擬

1、有限元大作業(yè)鋼架結(jié)構(gòu)分析選題人：日 期：2016年6月2日目錄：第一章：問題重述1一、題目內(nèi)容：1二、題目要求：1第二章：有限元法手工求解2一、平面兩單元離散化2二、單元分析2三、單元組裝5四、邊界條件引入及組裝總體方程5五、求解整體剛度方程，計算節(jié)點2的位移和轉(zhuǎn)角6六、求節(jié)點1、3支撐反力6七、設(shè)定數(shù)據(jù)，求解結(jié)果7八、繪制軸力圖、彎矩圖、剪力圖8第三章、matlab編程求解：9一、總體流程圖繪制：9二、輸入數(shù)據(jù)：9三、計算單元剛度矩陣：10四、建立總體剛度矩陣：10五、計算未約束點位移：10六、計算支反力：10七、輸出數(shù)據(jù)：10八、編程：10第四章 有限元求解11一、預(yù)處理11二、模型建立：

2、12二、分析計算14三、求解結(jié)果15四、繪制圖像16第五章 結(jié)果比較19第六章 心得體會19一、王小燦：19二、孫明哲：20三、張國威20第七章 附錄22一、matlab程序22第一章：問題重述一、題目內(nèi)容：圖示平面鋼架結(jié)構(gòu)圖1.1 題目內(nèi)容二、題目要求：（1）采用平面梁單元進行有限元法手工求解，要求寫出完整的求解步驟，包括： a）離散化：單元編號、節(jié)點編號；b）單元分析：單元剛度矩陣，單元節(jié)點等效載荷向量；c）單元組長：總體剛度矩陣，總體位移向量，總體節(jié)點等效載荷；d）邊界條件的引入及總體剛度方程的求解；e）b點的位移，a、c處支撐反力，并繪制該結(jié)構(gòu)的彎矩圖、剪力圖和軸力圖。（2）編制通用平

3、面鋼架分析有限元matlab程序，并計算蓋提，與手工結(jié)果進行比較；（3）利用ansys求解，表格列出b點的位移，a、c處支反力，繪制彎矩圖、剪力圖和軸力圖，并與手算和matlab程序計算結(jié)果比較。（4）攥寫報告，利用a4紙打??；（5）心得體會，并簡要說明各成員主要負責(zé)完成的工作。第二章：有限元法手工求解一、平面兩單元離散化將平面梁離散為兩個單元，單元編號分別為和，節(jié)點號分別為1、2、3；如圖2-1所示：圖2-1 單元離散化示意圖二、單元分析首先建立整體坐標系與局部坐標系如圖所示；1、求單元剛度矩陣對于單元，求局部坐標系的單元剛度矩陣：由于單元局部坐標系與整體坐標系的夾角為：，則單元的局部坐標變

4、換矩陣為：可以得到在總體坐標系下的單元的剛度矩陣：對于單元，求局部坐標系的單元剛度矩陣：由于單元局部坐標系與整體坐標系的夾角為,則。2、求單元節(jié)點等效載荷向量將p等效在單元兩側(cè)節(jié)點1,2上：將均布載荷等效在單元兩側(cè)的節(jié)點2，3上：與作用在節(jié)點上的力疊加為整體坐標系下的節(jié)點載荷：三、單元組裝將兩個整體坐標系下的單元剛度矩陣組裝為整體剛度矩陣：四、邊界條件引入及組裝總體方程由于節(jié)點1、3為固定約束，所以節(jié)點1和3的x、y方向的位移以及轉(zhuǎn)角均為0，節(jié)點2無位移約束，不存在支反力，所以力約束即為外力約束。 五、求解整體剛度方程，計算節(jié)點2的位移和轉(zhuǎn)角提取節(jié)點2位移的相關(guān)要素：求得：六、求節(jié)點1、3支撐

5、反力根據(jù)總體方程，提取求解節(jié)點1支撐反力所需方程：根據(jù)總體方程，提取求解節(jié)點2支撐反力所需方程：七、設(shè)定數(shù)據(jù)，求解結(jié)果設(shè)定各個數(shù)據(jù)：楊氏模量：泊松比：力：截面面積：慣性矩：將數(shù)據(jù)代入結(jié)果。節(jié)點2的位移和轉(zhuǎn)角：節(jié)點1支撐反力：節(jié)點3支撐反力：八、繪制軸力圖、彎矩圖、剪力圖應(yīng)用材料力學(xué)的分析方法，對梁單元進行分析。軸力圖：圖2-2 軸力圖剪力圖：圖2-3 剪力圖彎矩圖圖2-4 彎矩圖第三章、matlab編程求解：一、總體流程圖繪制：圖 3.1 總體流程圖二、輸入數(shù)據(jù)：考慮到后續(xù)計算和以下參數(shù)相關(guān)：節(jié)點個數(shù)，單元數(shù)，楊氏模量，慣性矩，單元長度，單元截面積，單元的旋轉(zhuǎn)角度，節(jié)點與單元的對應(yīng)關(guān)系，力與轉(zhuǎn)

6、矩的約束以及結(jié)構(gòu)約束。考慮到鋼架結(jié)構(gòu)，每個單元的楊氏模量，慣性矩，單元長度，單元截面積以及單元的旋轉(zhuǎn)角度都可能不一樣，所以采用矩陣的形式進行輸入。（注：由于本題除長度外一樣，故將其余幾項改為常量進行計算）單元與節(jié)點對應(yīng)關(guān)系為：一個單元對應(yīng)2個節(jié)點，且按順序連接。力與轉(zhuǎn)矩的約束以及結(jié)構(gòu)約束：應(yīng)包括約束值，作用節(jié)點，作用類型，3種，并以作用節(jié)點與作用類型來反推此約束在完整的約束矩陣中的位置。三、計算單元剛度矩陣：圖3.2 單元剛度矩陣生成流程圖考慮到每個單元的剛度矩陣與坐標變換的矩陣形式相同，只是數(shù)據(jù)不同，故采取建立模板，利用eval()，函數(shù)來帶入不同單元的值，生成一系列單元剛度矩陣，并用一個三

7、維數(shù)組存儲這些矩陣。四、建立總體剛度矩陣：考慮到每個單元剛度矩陣都是66的形式，表述了2個節(jié)點間的相互關(guān)系；故建立元胞數(shù)組，并使元胞數(shù)組的階數(shù)與節(jié)點個數(shù)相同，利用元胞數(shù)組存儲節(jié)點間關(guān)系。首先建立與節(jié)點個數(shù)相同階數(shù)的空元胞數(shù)組，之后檢索每個單元剛度矩陣對應(yīng)的2個節(jié)點間的關(guān)系，將其分離成4個33的矩陣，按節(jié)點與單元對應(yīng)關(guān)系，存儲到元胞數(shù)組中。最后將元胞數(shù)組展開形成的大矩陣即為總體剛度矩陣。五、計算未約束點位移：利用總體位移與外力間的關(guān)系，采用矩陣求解，求取非約束點的位移。并針對結(jié)果進行對應(yīng)處理，使結(jié)果與作用點、作用形式對應(yīng)。六、計算支反力：利用約束點位移皆零的特點，簡化總體剛度矩陣，同時由于部分節(jié)

8、點的部分方向上為內(nèi)力而非支反力，再度簡化總體剛度矩陣。利用兩次簡化后的剛度矩陣與計算出的位移結(jié)果相乘，求得不計直接作用在節(jié)點約束方向上時的支反力，將結(jié)果加上由于直接作用在節(jié)點約束方向上時產(chǎn)生的支反力，即為最后的支反力結(jié)果。七、輸出數(shù)據(jù)：將計算所得的未約束點位移與支反力，采用與輸入方式相似的方式進行處理并進行輸出。八、編程：見附錄一第四章 有限元求解一、預(yù)處理1、選擇單元類型：ansys main menu: preprocessor element typeadd/edit/delete add beam：2d elastic 3ok (返回到element types 窗口) close圖4

9、.1 選擇單元類型2、定義材料參數(shù)：ansys main menu: preprocessor material props material modelsstructural linear elasticisotropic: ex:3e10 (彈性模量)，prxy:0.3(泊松比)ok 圖4.2定義材料參數(shù)3、定義單元截面積和慣性矩：ansys main menu: preprocessor real constant add type beam 3 ok cross-sectional area area:0.05(橫截面積) area moment of inteia izz:1(慣性矩

10、) ok圖4.3定義單元截面積和慣性矩二、模型建立：1、畫出關(guān)鍵點：ansys main menu: preprocessor modeling creatkeypoint in active csnode number 1 x:0,y:0,z:0 apply node number 2 x:0,y:1,z: apply node number 3 x:2,y:1,z:0ok 2、構(gòu)造連線：ansys main menu: preprocessor modeling creatline linesstraight line 依次連接特征點ok圖4.4模型建立3、劃分網(wǎng)格：ansys main

11、menu: preprocessor meshingmeshtool set 選擇1，2節(jié)點之間部分apply選擇2，3節(jié)點之間部分單元長度分別為0.1和0.2okmeshingmeshtool mesh分別選擇1和2，2和3節(jié)點之間部分ok圖4.4 劃分網(wǎng)格4、添加約束和載荷： 左下角和右上角添加約束：ansys main menu: preprocessor solutiondefine loads apply structural displacement on nodes 選擇1節(jié)點all dofapplyon nodes 選擇1節(jié)點all dofok添加頂部均布載荷：ansys ma

12、in menu: preprocessor solutiondefine loads apply structural pressure on beams 選擇頂部所有的單元vali pressure value node i :1000 valj pressure value node j :1000 ok添加力矩和力：ansys main menu: preprocessor solutiondefine loads apply structural force/ monmenton nodes 選擇2節(jié)點apply lab mzvalue 100 . (輸入力矩)on nodes 選擇8

13、節(jié)點apply lab fx value 1000 （輸入力）圖4.5 添加約束和載荷二、分析計算ansys main menu: solution solve current ls ok should the solve command be executed? y close (solution is done! ) 關(guān)閉窗口圖4.6 求解模型三、求解結(jié)果1、位移ansys main menu: general postproc list resultnodal solutiondof solutionx-component of displacement apply y-componen

14、t of displacementok圖4.7 x方向位移解圖4.8 y方向位移解2、支反力：ansys main menu: general postproc list resultreaction soluall itemsok圖4.9 支反力結(jié)果四、繪制圖像1、設(shè)置參數(shù)ansys main menu: general postproc element table difine table add 在user label for item 中輸入fx-i, 在results data item中選擇by sequence num并輸入 smisc,1apply在user label for

15、 item 中輸入fx-j, 在results data item中選擇by sequence num,并輸入 smisc,7apply在user label for item 中輸入fy-i, 在results data item中選擇by sequence num,并輸入 smisc,2apply在user label for item 中輸入fy-j, 在results data item中選擇by sequence num,并輸入 smisc,8apply在user label for item 中輸入mz-i, 在results data item中選擇by sequence num

16、,并輸入 smisc,6apply在user label for item 中輸入mz-j, 在results data item中選擇by sequence num,并輸入 smisc,12ok圖4.10 圖表參數(shù)設(shè)置2、圖像輸出a）軸力圖：ansys main menu: general postproc plot resul tcontour plot line elem res 選擇 fx_i fx_j apply圖4.11 軸力圖b）剪力圖：ansys main menu: general postproc plot resul tcontour plot line elem res

17、 選擇 fy_i fy_j apply圖4.12 剪力圖c）彎矩圖：ansys main menu: general postproc plot resul tcontour plot line elem res 選擇 mz_i z_j ok圖4.13 彎矩圖第五章 結(jié)果比較結(jié)果手算matlabansys-0.4142210-8-0.4142210-8-0.4142210-8-0.3302310-7-0.3302310-7-0.3302310-7-1033.1-1033.1-1033.149.53549.53549.53564.50164.50164.5013.10673.10673.1067

18、1950.51950.51950.5-1462.3-1462.3-1462.3通過對比知道，三種方式的結(jié)果完全一樣，顯示了結(jié)果的正確性。第六章 心得體會在有限元課程結(jié)束之際，通過這份大作業(yè)，我整理總結(jié)了這個學(xué)期學(xué)到的有限元思想。深化我對有限元計算流程的理解，提升了我的matlab的編程水平。編程時首先要對整體的流程有一個清晰地構(gòu)想，從數(shù)據(jù)的獲取開始到計算輸出進行分步處理；通過手算例題確定計算過程中使用的參數(shù)。再通過模板獲得每個單元的矩陣。難點在于如何將這些單元正確的拼接到一起，一開始用的是biadiag對角拼接指令，但由于每個單元對應(yīng)的節(jié)點并非按順序排列，不能實現(xiàn)所需效果；之后我重新計算例題，

19、探究到組裝的本質(zhì)是將每個單元的矩陣按照矩陣所對應(yīng)節(jié)點的信息進行分塊，并累加到總體矩陣中，由于矩陣實現(xiàn)較為麻煩，所以我想利用元胞數(shù)組來存儲節(jié)點信息，并成功實現(xiàn)了所需功能，又利用cat拼接指令，將元胞數(shù)組展開成一個大型矩陣完成拼接。再獲得總體剛度矩陣后，只要找到對應(yīng)的位移與外力的參數(shù)行，再進行乘除計算就相對簡單了。我認為在學(xué)習(xí)過程中，提高編程能力是很重要的一個方面，有著良好的編程能力，可以讓很多工程問題得以用計算機解決，更容易獲得結(jié)果。而想要提高編程能力，首先要有一定的編程思想與數(shù)學(xué)建模能力。我認為編程思想指的是，要對所做的事情有一個過程性與結(jié)構(gòu)性的認識，并根據(jù)使用對象進行相關(guān)的調(diào)整。過程性指的是

20、在編程要有一個完整的流程圖，將復(fù)雜問題轉(zhuǎn)化成若干簡單的小問題，針對小問題進行求解；結(jié)構(gòu)性指的是在編程中要利用條件if，循環(huán)for等結(jié)構(gòu)簡化程序，同時也要對每部分程序的輸入、輸出以及執(zhí)行作用有一個明確的認知；最后根據(jù)使用對象，對輸入、輸出，進行調(diào)整，也可以用gui編制界面方便使用。以上是我對這次大作業(yè)的心得體會，希望老師輔導(dǎo)校正第七章 附錄一、matlab程序clcclearformat compactformat shortgjd=input(請輸入節(jié)點數(shù)：);dy=input(請輸入單元數(shù)：);e=input(請輸入楊氏模量e：);i=input(請輸入慣性矩i：);l=input(請輸入單

21、元長度l：);a=input(請輸入單元截面積：);fai=input(請輸入單元相對旋轉(zhuǎn)角度：);%輸入對應(yīng)關(guān)系時，小節(jié)點放前面單元 節(jié)點1 節(jié)點2dy_jd=input(請輸入單元與節(jié)點對應(yīng)關(guān)系：);%輸入力與扭矩約束值 作用節(jié)點 作用類型（轉(zhuǎn)矩為3 x方向為1 y方向為2）lys=input(力與轉(zhuǎn)矩約束矩陣：);%輸入結(jié)構(gòu)約束作用節(jié)點 作用類型（轉(zhuǎn)角為3 x方向為1 y方向為2）wys=input(結(jié)構(gòu)約束矩陣：);%原始數(shù)據(jù)% l=1;% e=3*1010;% p=1000;% a=0.05;% dy=2;jd=3;ll=l 2*l;i=20*a;% dy_jd=1 1 2;2 2

22、3;% fai=pi/2 0;% q=p/l;m=p*l/10;% lys=44/125*p 1 1;-12*p*l/125 1 3;81/125*p 2 1;-p 2 2;-67/750*p*l 2 3;-p 3 2;p*l/3 3 3;% wys=1 1;1 2;1 3;3 1;3 2;3 3;%對力約束與位移約束式子分別進行編號處理wys(:,3)=(wys(:,1)-1)*3+wys(:,2);lys(:,4)=(lys(:,2)-1)*3+lys(:,3);%對力約束與位移約束式子進行排序lys=sortrows(lys,4);wys=sortrows(wys,3);%單元剛度矩陣s

23、yms fai e a i l realk=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;t= cos(fai), sin(fai), 0; -sin(fai), cos(fai), 0; 0, 0, 1;%坐標變換矩陣

24、t=blkdiag(t,t);%總體坐標系下的單元剛度矩陣k=t*k*t;%帶入每個單元的數(shù)，生成單元剛度矩陣kk，其每一頁對應(yīng)相應(yīng)頁數(shù)的單元的剛度矩陣for j=1:dy; e=e; i=i; l=ll(j); a=a; fai=fai(j); kk(:,:,j)=eval(k);end%生成總體剛度矩陣kk%采用元胞數(shù)組的方式對各項進行保存%生成空元胞數(shù)組，元胞數(shù)組的行列大小與節(jié)點數(shù)相同for j=1:jd; for jj=1:jd; ling1j,jj=zeros(3); endendling2=ling1;%將對單元剛度矩陣部分分成4分加入元胞數(shù)組中for j=1:dy; kk1=kk(1:3,1:3,j); kk2=kk(1:3,4:6,j); kk3=kk(4:6,1:3,j); kk4=kk(4:6,4:6,j); ling2dy_jd(j,2),dy_jd(j,2)=kk1+ling2dy_jd(j,2),dy_jd(j,2); ling2dy_jd(j,2),dy_jd(j,3)=kk2+ling2dy_jd(j,2),dy_jd(j,3); ling2dy_jd(j,3),dy_jd(j,2)=kk3+ling2dy_jd(j,3),dy_jd(j,2); ling2dy_jd(j,3),dy_jd(j,3)=kk4+