有限元钢架结构分析手算+matlab+ansys模拟

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%采用元胞数组