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

2、2二、分析計算14三、求解結果15四、繪制圖像16第五章結果比較19第六章心得體會19一、王小燦：錯誤!未定義書簽。二、孫明哲：錯誤!未定義書簽。三、張國威錯誤!未定義書簽。第七章附錄21一、matlab程序21、題目內容:圖示平面鋼架結構第一章：問題重述Eh20EA 1=100cm0 41圖i.i題目內容、題目要求:1采用平面梁單元進行有限元法手工求解，要求寫出完整的求解步驟，包括：a離散化：單元編號、節點編號；b單元分析：單元剛度矩陣，單元節點等效載荷向量；c單元組長：總體剛度矩陣，總體位移向量，總體節點等效載荷；d邊界條件的引入及總體剛度方程的求解；eB點的位移，A、C處支撐反力，并繪制

3、該結構的彎矩圖、剪力圖和軸力圖。2編制通用平面鋼架分析有限元Matlab程序，并計算蓋提，與手工結果進行比較；3利用Ansys求解，表格列出B點的位移，A、C處支反力，繪制彎矩圖、剪力圖和軸力圖，并與手算和Matlab程序計算結果比較。4攥寫報告，利用A4紙打印；5心得體會，并簡要說明各成員主要負責完成的工作。3第二章：有限元法手工求解、平面兩單元離散化將平面梁離散為兩個單元，單元編號分別為和，節點號分別為1、2、3;如圖2-1所示:圖2-1單元離散化示意圖、單元分析首先建立整體坐標系與局部坐標系如下圖;1、求單元剛度矩陣對于單元，求局部坐標系的單元剛度矩陣:K1EAl00EAl00012EI

4、丁6EI了012EIl36EI06EIl24EIl06EI7T2EIlEAl00EAl00012EI6EIF012EIl36EIF06EIl22EIl06EI7r4EIl則單元的局部坐標變換矩陣為：cossin0000sincos0000T1e001000000cossin0000sincos0000001由于單元局部坐標系與整體坐標系的夾角為:90可以得到在總體坐標系下的單元的剛度矩陣:01000010000000100000001000010000000112EIK1T1eTK；T1e6EIl212EI不6EI7T06EI12EI06EIl2l3l2EAl00EAl0104EI6EI103

5、EIlll206EI12EI06EIl2l3l2EAl00EAl0103EI6EI104EI2對于單元，求局部坐標系的單元剛度矩陣:EA2lEA2l3EI3EI2l33EI2l22EI3EI2l3K22l23EI2l2EA2lEA2l3EI3EI2l33EI2l2EI3EI2l22l33EI2l203EI2l2EIl03EI2l22EIl由于單兀局部坐標系與整體坐標系的夾角為0,則K2K2。2、求單元節點等效載荷向量將P等效在單元兩側節點1,2上：11Fx1Fx2018F；F；244 P12581P125Fb23abL3Fa2a3bL3M11Fab212plL2725m2將均布載荷等效在單元兩

6、側的節點2, 3 上:2Fa2b18pl2-L125Fx2,FX30xvxF；2Fy23m1 m2pi212Pl3與作用在節點上的力疊加為整體坐標系下的節點載荷:x144 p125Fy10M112pl125Fx281P125)Fy2PM218125110Pl67Pl750(Fx3Fy1PplM33三、單元組裝12EI06EI12EI06EI000l3l2l3l2EAEA0l00l00006EI104EI6EI103EIl000l2ll212EI06EI12EIEA06EIEA00l3l23l2-l32l2l0EAl00EA3E1l2l33EI2l203EI-2l33EI2l26EI02EI6E

7、I3EI6EIl03EIEIl2ll22l2-2l2l000EA00EA00-2l2l00003EI3EI03EI3EI-2l32l22l32l200003EIEI03EI2EI2l2l2l2l將兩個整體坐標系下的單元剛度矩陣組裝為整體剛度矩陣:K四、邊界條件引入及組裝總體方程由于節點1、3為固定約束，所以節點不存在支反力，所以力約束即為外力約束。1和3的x、y方向的位移以及轉角均為0,節點2無位移約束,U1ViU2V2M1Qx2Qy2M2Fx2Fy2M281p125P67PlU3750U3根據總體方程，提取求解節點1支撐反力所1希方程:12EI i 3-06EI , 2一'Qx11E

8、A11U2Fx1 一'Qy100V2F y1M16EI1203EI21M1根據總體方程，EA提取求解節點02支撐反力所需方程:0,Qx3213EI53EI -Fx3Qy30213了V2212Fy3M303EI212EI21M3六、求節點1、3支撐反力五、求解整體剛度方程，計算節點2的位移和轉角提取節點2位移的相關要素：12EIEA06EI一81p132112U2Qx2F，212pEA3EI3EI125'0，CI32V2Qy2Fy2P6EI1213EI216EI2MM'67Pl122121Ivl2Ivl2750求得：1413P(316A12273I)U.375E(8A21

9、4105AI1272I2)232713P(419A124080I)2/24227375E(8A1105AI172I)212P(134A21412831AI1231680I2)1125EI(8A214105AI1272I2)2P(394A21420643AIl213500I2)3758A214105AI1272I2)7A12P(419A124080I)3758A214105AI1272I2),24222P1(298A132655AI11170I)11258A214105AI1272I2)7A12P(316A12273I)3758A214105AI1272I2)P(3067A21450190AI1

10、25400012)3758A214105AI1272I2),2422P1(6268A1130809AI1247680I)22508A214105AI127212)七、設定數據，求解結果設定各個數據：楊氏模量：E3l010Pa泊松比：0.3力：P1KN2截面面積：A0.05m慣性矩：I1m4將數據代入結果。節點2的位移和轉角:1413P(316A12273I)U2375E(8A214105AI1272I2)713P(419A124080I)V2/242Z2T375E(8A1105AI172I)212P(134A21412831AI1231680I2)1125EI(8A214105AI1272I2

11、)節點1支撐反力：2R394A21420643AI1213500I2)Q375(8A214105AI1272I2)Qx1Q7A12R419A124080I)y1375(8A214105AI127212)M12P1(298A21432655AIl21170I2)1125(8A214105AI1272I2)節點3支撐反力：0.0414-7100.33020.11901003.149.564.5Qx3Qy3M37A12P(316A12273I)124227375(8A1105AI172I)P(3067A21450190AIl254000I2)375(8A214105AI1272I2)Pl(6268A

12、214130809AI12247680I2)2250(8A214105AI1272I2)3.11950.51462.3八、繪制軸力圖、彎矩圖、剪力圖應用材料力學的分析方法，對梁單元進行分析。軸力圖:剪力圖:彎矩圖圖2-2軸力圖圖2-3剪力圖圖2-4彎矩圖第三章、matlab編程求解:、總體流程圖繪制:開始輸入數據計算單元喇度矩陣求氈體剛度矩陣計算未約束點位移計算支反力輸出數據結束圖3.1總體流程圖二、輸入數據：考慮到后續計算和以下參數相關：節點個數，單元數，楊氏模量，慣性矩，單元長度，單元截面積，單元的旋轉角度，節點與單元的對應關系，力與轉矩的約束以及結構約束。考慮到鋼架結構，每個單元的楊氏模

13、量，慣性矩，單元長度，單元截面積以及單元的旋轉角度都可能不一樣，所以采用矩陣的形式進行輸入。注：由于此題除長度外一樣，故將其余幾項改為常量進行計算單元與節點對應關系為：一個單元對應2個節點，且按順序連接。力與轉矩的約束以及結構約束：應包括約束值，作用節點，作用類型，3種，并以作用節點與作用類型來反推此約束在完整的約束矩陣中的位置。三、計算單元剛度矩陣:Jiten建立修植K里標年授U輸入京據摩3在曜疵的13結天圖3.2單元剛度矩陣生成流程圖考慮到每個單元的剛度矩陣與坐標變換的矩陣形式相同，只是數據不同，故采取建立模板，利用eval(),函數來帶入不同單元的值，生成一系列單元剛度矩陣，并用一個三維

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

15、節點的部分方向上為內力而非支反力，再度簡化總體剛度矩陣。利用兩次簡化后的剛度矩陣與計算出的位移結果相乘，求得不計直接作用在節點約束方向上時的支反力，將結果加上由于直接作用在節點約束方向上時產生的支反力，即為最后的支反力結果。七、輸出數據：將計算所得的未約束點位移與支反力，采用與輸入方式相似的方式進行處理并進行輸出。八、編程：見附錄一第四章有限元求解、預處理1、選擇單元類型：fElement Type f Add/Edit/Delete-Add-2D-elbstacr13一 OK （ANSYSMainMenu:Preprocessor回至ijElementTypes窗口）一Close圖4.1選擇

16、單元類型2、定義材料參數:ANSYSMainMenu:Preprocessor-MaterialProps-MaterialModels-Structural-LinearfElasticIsotropic:EX:3e10（彈性模量）,PRXY:0.3（泊松比）-OK出卡24ANSY, d JLM 刖Et M-W工 iK4repirMwM. feWa*!-.E bum mJ FLk fir * UHIW. |F Amc-rECE 抵 NETfcTM- * C*w . Tlwui LHHftca_Z_* FrlrUjmnii d田叫STI GT! tef Up3J MidWMd hLLpE田 P

17、hywRHI *E If-filMlI AlUlFMncbfcat FdrD= IXUMtiH iMMftaCPFAMMd Kil llttfiMAi 用"111回里匈皿國國J1JOJ國日用國包回7nltJeJTI圖4.2定義材料參數3、定義單元截面積和慣性矩:ANSYSMainMenu:PreprocessorfRealconstantfAddfTypebeam3OkfCross-sectionalareaAREA:0.05（橫截面積）AreamomentofinteiaIZZ:1（慣性矩）一OK圖4.3定義單元截面積和慣性矩、模型建立:1、畫出關鍵點:ANSYSMainMenu:

18、PreprocessorfModelingCreatKeypointInActiveCSfNodenumber1一X:0,Y:0,Z:0一ApplyNodenumber2一X:0,Y:1,Z:一ApplyNodenumber3一X:2,Y:1,Z:0一OK2、構造連線：ANSYSMainMenu:PreprocessorModelingCreat-Lineflinesfstraightline一依次連接特征點一Ok圖4.4模型建立ANSYSMainMenu:PreprocessorMeshingMeshtoolSet一選擇1,2下點之間部分一Apply一選擇2,3節點之間部分一單元長度分別為和

19、一OKMeshingMeshtoolMesh一分別選擇1和2,2和3節點之間部分一OKiArt.DU Rf y_DH d 甘 | a&'fllr&F PH巫nz JRn|u ，后閉二| a|£用噩巧岫tin Hg網即tftenEE 圖 9p*國三 網B*圖4.4劃分網格3«lo|4、添加約束和載荷：左下角和右上角添加約束：ANSYSMainMenu:Preprocessor-Solution-Defineloads-Apply-Structural-Displacement-Onnodes一選擇1節點一ALLDOFfApplyOnnodes一選擇1節

20、點一ALLDOF一OK添加頂部均布載荷：ANSYSMainMenu:Preprocessor-Solution-Defineloads-Apply-Structural-Pressure-Onbeams-選擇頂部所有的單元fVALIpressurevaluenodeI:1000VALJpressurevaluenodeJ:1000fOK添加力矩和力：ANSYSMainMenu:PreprocessorfSolutionfDefineloadsApplyStructuralfForce/MonmentfOnnodes一選擇2節點一ApplyLABMZVALUE100.（輸入力矩）一Onnodes

21、一選擇8節點一ApplyLABFXVALUE1000輸入力J| D_ _:，一 .1. 一。卜包堡1回創凰劇工AMSVS raolb»圖4.5添加約束和載荷、分析計算ANSYSMainMenu:SolutionSolveCurrentLS一OK一ShouldtheSolveCommandbeExecuted?丫一Close(Solutionisdone!)一關閉窗口32圖4.6求解模型三、求解結果1、位移solution f X-component ofANSYSMainMenu:GeneralPostprocListresultNodalsolution一DOFdisplaceme

22、ntfApplyfY-componentofdisplacementfOKA具峭WedqHRF“VL叫“也"u.amay國,APQLg&N爐EM*5mNSlDE-fl.Z-WS-lF-W4lsl%4lfL11.d,匿療14ye.UTVSLtt¥-*.*1*11-*-e.i-fl.tijminfl.MilHFMHffl*B-fli/FMF-BRfl/mw.由NubenI“I*FftDD*u&IftX*EFilhinStlM”一IrRy,OGMUiSuma5Ql«riHQ>*M4ilSiBlM<.SHinr-r，r!，!"二L&#

23、171;tU4ryiMUlllQDptenitorOwtpDMwgwnwPWVMtEHiiidddl"圖4.7x方向位移解小AIMW£0LUrtyIlerIFil*MmIm，舊PltCykWoikllirfP*1iEmMtmMiruClrKWes畫回國創囹工國AhSYEiDQd-書制S3E-R5>dIuLi>nSokiSio一mi9E5-W.2*721L-4i>7NodalilLNdiM2-B-M3MM-97RLNRIIIF印7f6MlM歸一迪占帕處射JNMEflr7T.El&E-SIHft；lMbIJTsow項GENU-二一maii一二一二二三E

24、lFrefereritEiEHlPr(prMtWprfgluliwi臼GeneraliMtpirDc50ata&FIe:Opts3fcMuht(urnmiry國融*Mu3FxihireCFitaria3fQLMfovuHiatera-bonEumniqP4rt|4lErT$F*aX3LiU；咐iS口jticVDctorDa>taiPMhiMtrtilie-T；HflH-3QuwyResults-5&0i中SfttiuluVllWff5WFfttKAFUQhodnd亡alaJ233711F-W7-B；1T251E-0T75HmE-的n圖4.8y方向位移解2、支反力：ANSY

25、SMainMenu:GeneralPostprocListresultReactionSolufAllitems一OK圖4.9支反力結果四、繪制圖像1、設置參數ANSYSMainMenu:GeneralPostprocElementTableDifineTable一Add在userlabelforitem中輸入FX-I,在Resultsdataitem中選擇Bysequencenum并輸入smisc,1-Apply在userlabelforitem中輸入FX-J,在Resultsdataitem中選擇Bysequencenum,并輸入smisc,7-Apply在userlabelforitem

26、中輸入FY-I,在Resultsdataitem中選擇Bysequencenum,并輸入_smisc,2-Apply在userlabelforitem中輸入FY-J,在Resultsdataitem中選擇Bysequencenum,并輸入smisc,8-Apply在userlabelforitem中輸入MZ-I,在Resultsdataitem中選擇Bysequencenum,并輸入smisc,6-Apply在userlabelforitem中輸入MZ-J,在Resultsdataitem中選擇Bysequencenum,并輸入smisc,12fOKTirwwapip Tlrw iJWQI n

27、(w= uaMfl 靠一 2JfiMa TimTirwr 1.MC0SMjs 帽gg 了 sum) CwmnC 求kmw / 1arrihBUJ印百ANSYSffif 1 徹* LlhlEiBI""fZndli irirvuirT3 Id 聃IrhUft D RriptnvrHiiflw E !WFQ DMdjim ihwiTiVk三F力 Cpermms圖4.10圖表參數設置2、圖像輸出a軸力圖:ANSYSMainMenu:GeneralPostprocPlotresultContourplotfLineElemRes一選擇FX_IFX_J一Apply圖4.11軸力圖b男力

28、圖:ANSYSMainMenu:GeneralPostprocPlotresultContourplotfLineElemRes一選擇FY_IFY_J一ApplyA靦05ft加如期UPUkiPkqQcrhWQfUFWnt修KeroMtyBls電一Rl更魚到倒倒Jt回UX理颯中sMix絲巴帆弘JMJMuuri|ta由 Eprw.iwS WM1MB 3 %l Pwlw S ML1Q Mt O h 母 trMtt !kWn«np t田目 Vlitari Crlwii El W：mwitoS Dflfw mril lh«pi C4flWJr Hm3 h* dii M 加 iS Ih

29、miRt Mm S Hi hi f由 Ymhh rkit QPWWMmB Cmu«1a Fl At STMunimt Um Niidi. S flHduhi Ei .SHS lor H gwtr H做訕r« rat 匕 »Md Cilu E fltmtni. T j E full /*NU Iuf4i»E ia*d C«m 出 Chwi. ttaiHm *i TD£-£HK -49CD rim。圖4.12男力圖c彎矩圖:ANSYSMainMenu:GeneralPostprocPlotresultContourplotfLi

30、neElemRes一選擇MZ_IZ_J一OK圖4.13彎矩圖第五章結果比較結果手算MATLABANSYSU2XX10-8X10-8X10-8U2YX10-7X10-7X10-7FAXFAYMaFCXFCYMe通過比照知道，三種方式的結果完全一樣，顯示了結果的正確性。第六章心得體會在有限元課程結束之際，通過這份大作業，我整理總結了這個學期學到的有限元思想。深化我對有限元計算流程的理解，提升了我的matlab的編程水平。編程時首先要對整體的流程有一個清晰地設想，從數據的獲取開始到計算輸出進行分步處理；通過手算例題確定計算過程中使用的參數。再通過模板獲得每個單元的矩陣。難點在于如何將這些單元正確的拼

31、接到一起，一開始用的是biadiag對角拼接指令，但由于每個單元對應的節點并非按順序排列，不能實現所需效果；之后我重新計算例題，探究到組裝的本質是將每個單元的矩陣按照矩陣所對應節點的信息進行分塊，并累加到總體矩陣中，由于矩陣實現較為麻煩，所以我想利用元胞數組來存儲節點信息，并成功實現了所需功能，又利用cat拼接指令，將元胞數組展開成一個大型矩陣完成拼接。再獲得總體剛度矩陣后，只要找到對應的位移與外力的參數行，再進行乘除計算就相對簡單了。我認為在學習過程中，提高編程能力是很重要的一個方面，有著良好的編程能力，可以讓很多工程問題得以用電腦解決，更容易獲得結果。而想要提高編程能力，首先要有一定的編程

32、思想與數學建模能力。我認為編程思想指的是，要對所做的事情有一個過程性與結構性的認識，并根據使用對象進行相關的調整。過程性指的是在編程要有一個完整的流程圖，將復雜問題轉化成假設干簡單的小問題，針對小問題進行求解；結構性指的是在編程中要利用條件if,循環for等結構簡化程序，同時也要對每部分程序的輸入、輸出以及執行作用有一個明確的認知；最后根據使用對象，對輸入、輸出，進行調整，也可以用GUI編制界面方便使用。以上是我對這次大作業的心得體會，希望老師輔導校正第七章附錄一、matlab程序clcclearformatcompactformatshortGjd=input（'請輸入節點數：

33、9;）;dy=input（'請輸入單元數：'）;E=input（'請輸入楊氏模量E：'）;I=input（'請輸入慣性矩I：'）;L=input（'請輸入單元長度L：'）;A=input（'請輸入單元截面積：'）;FAI=input（'請輸入單元相對旋轉角度：'）;%輸入對應關系時，小節點放前面單元節點1節點2dy_jd=input（'請輸入單元與節點對應關系：'）;%輸入力與扭矩約束值作用節點作用類型轉矩為3x方向為1y方向為2lys=input（'力與轉矩約束矩陣：

34、9;）;%輸入結構約束作用節點作用類型轉角為3x方向為1y方向為2wys=input（'結構約束矩陣：'）;%原始數據%L=1;%E=3*10A10;%P=1000;%A=0.05;%dy=2;jd=3;LL=L2*L;I=20*A;%dy_jd=112;223;%FAI=pi/20;%q=P/L;M=P*L/10;%lys=44/125*P11;-12*P*L/12513;81/125*P21;-P22;-67/750*P*L23;-P32;P*L/333;%wys=11;12;13;31;32;33;%對力約束與位移約束式子分別進行編號處理wys（:,3）=（wys（:,1

35、）-1）*3+wys（:,2）;lys（:,4）=（lys（:,2）-1）*3+lys（:,3）;%對力約束與位移約束式子進行排序lys=sortrows(lys,4);wys=sortrows(wys,3);%單元剛度矩陣symsfaieailrealk=e*a/l00-e*a/l00;012*e*i/F36*e*i/F20-12*e*i/F36*e*i/F2;06*e*i/lA24*e*i/l0-6*e*i/lA22*e*i/l;-e*a/l00e*a/l00;0-12*e*i/lA3-6*e*i/lA2012*e*i/lA3-6*e*i/lA2;06*e*i/lA22*e*i/l0-6*e*i/lA24*e*i/l;t=cos(fai),sin(fai),0;-sin(fai),cos(fai),0;0,0,1;%坐標變換矩陣T=blkdiag(t,t);%總體坐標系下的單元剛度矩陣K=T'*k*T;%帶入每個單元的數，生成單元剛度矩陣kk，其每一頁對應相應頁數的單元的剛度矩陣forj=1:dy;e=E;i=I;l=LL(j);a=A;fai=FAI(j);kk(:,:,j)=eval(K);end%生成總體剛度矩陣KK%采用元胞數組