西交大有限元原理及应用-大作业
有限元分析大作业报告
有限元分析大作业报告一、引言有限元分析是工程领域中常用的数值模拟方法，通过将连续的物理问题离散为有限个子区域，然后利用数学方法求解，最终得到数值解。
有限元分析的快速发展和广泛应用，为工程领域提供了一种强大的工具。
本报告将介绍在大作业中所进行的有限元分析工作及结果。
二、有限元模型建立本次大作业的研究对象是工程结构的应力分析。
首先，通过对结构进行几何建模，确定了结构的尺寸和形状。
然后，将结构离散为有限个单元，每个单元又可以看作一个小的子区域。
接下来，为了求解结构的应力分布，需要为每个单元确定适当的单元类型和单元属性。
最后，根据结构的边界条件，建立整个有限元模型。
三、材料属性和加载条件在建立有限元模型的过程中，需要为材料和加载条件确定适当的参数。
本次大作业中，通过实验获得了结构材料的弹性模量、泊松比等参数，并将其输入到有限元模型中。
对于加载条件，我们选取了其中一种常见的加载方式，并将其施加到有限元模型中。
四、数值计算和结果分析为了求解结构的应力分布，需要进行数值计算。
在本次大作业中，我们选用了一种常见的有限元求解器进行计算。
通过输入模型的几何形状、材料属性和加载条件，求解器可以根据有限元方法进行计算，并得到结构的应力分布。
最后，我们通过对计算结果进行分析，得出了结论。
五、结果讨论和改进方法根据计算结果，我们可以对结构的应力分布进行分析和讨论。
根据分析结果，我们可以得出结论是否满足设计要求以及结构的强度情况。
同时，根据分析结果，我们还可以提出改进方法，针对结构的特点和问题进行相应的优化设计。
六、结论通过对工程结构进行有限元分析，我们得到了结构的应力分布，并根据分析结果进行了讨论和改进方法的提出。
有限元分析为工程领域提供了一种有效的数值模拟方法，可以帮助工程师进行结构设计和分析工作，提高设计效率和设计质量。
【1】XXX，XXXX。
【2】XXX，XXXX。
以上是本次大作业的有限元分析报告，总结了在建立有限元模型、确定材料属性和加载条件、数值计算和结果分析等方面的工作，并对计算结果进行讨论和改进方法的提出。
有限元分析大作业
有限元大作业一题目要求：图1所示为一悬臂梁，在端部承受载荷，材料弹性模量为E，泊松比为1/3，悬臂梁的厚度（板厚）为t，若该粱被划分为两个单元，单元和节点编号如图所示，试按平面应力问题计算各个节点位移计支反力。
一、单元划分1.计算简图及单元划分如下所示：2.进行节点及单元编号节点i j m单元① 2 3 4② 3 2 13.节点坐标值节点号1 2 3 4坐标值X 2 2 0 0Y 1 0 1 0二、计算单元刚度矩阵1、计算每个单元面积△以及i b ，i c （m j i i ,,=） ①②单元的面积相等，即12121=⨯⨯=∆ 单元①的i b ，i c⎩⎨⎧=--==-=0)(1m j i m j i y x c y y b ⎩⎨⎧=--==-=2)(0i m ji m j x x c y y b ⎩⎨⎧-=--=-=-=2)(1j i mj i m y x c y y b 对平面应力问题，其表达式为[]⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡-+-+-+-+∆-=s r s r sr s r s r s r s r s r b b uc c cb u b uc b c u c ub c c u b b u Et Krs 21212121)1(42 然后对单元①求解单元刚度子矩阵2==i r 2==i s []⎥⎦⎤⎢⎣⎡=3/1001329)1(22Et K 2==i r 3==j s []⎥⎦⎤⎢⎣⎡=03/23/20329)1(23Et K2==i r 4==m s []⎥⎦⎤⎢⎣⎡----=3/13/23/21329)1(24Et K 3==j r 3==j s []⎥⎦⎤⎢⎣⎡=4003/4329)1(33Et K 3==j r 2==i s []⎥⎦⎤⎢⎣⎡=03/23/20329)1(32Et K 3==j r 4==m s []⎥⎦⎤⎢⎣⎡----=43/23/23/4329)1(34Et K 4==m r 4==m s []⎥⎦⎤⎢⎣⎡=3/133/43/43/7329)1(44Et K 4==m r 2==i s []⎥⎦⎤⎢⎣⎡----=3/13/23/21329)1(42Et K 4==m r 3==j s []⎥⎦⎤⎢⎣⎡----=43/23/23/4329)1(43Et K由子矩阵[]e rs K 合成单元刚度矩阵[]⎥⎥⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎢⎢⎣⎡----------------=3/133/443/23/13/23/43/73/23/43/2143/24003/23/23/403/43/203/13/203/23/103/213/2001329)1(Et K将单元①的单元刚度矩阵补零升阶变为单元刚度矩阵，其在总体刚度矩阵中的位置为：节点号→单元②的i b ，i c⎩⎨⎧=--=-=-=0)(1m j im j i y x c y y b ⎩⎨⎧-=--==-=2)(0i m ji m j x x c y y b ⎩⎨⎧=--==-=2)(1j i mj i m y x c y y b 然后对单元 求解单元刚度子矩阵：3==i r 3==i s []⎥⎦⎤⎢⎣⎡=3/1001329)2(33Et K 3==i r 2==j s []⎥⎦⎤⎢⎣⎡=03/23/20329)2(32Et K 3==i r 1==m s []⎥⎦⎤⎢⎣⎡----=3/13/23/21329)2(31Et K 1 2 3 412[])1(22K[])1(23K[])1(24K3[])1(32K[])1(33K[])1(34K4[])1(42K[])1(43K[])1(44K2==j r 2==j s []⎥⎦⎤⎢⎣⎡=4003/4329)2(22Et K 2==j r 3==i s []⎥⎦⎤⎢⎣⎡=03/23/20329)2(23Et K 2==j r 1==m s []⎥⎦⎤⎢⎣⎡----=43/23/23/4329)2(21Et K 1==m r 1==m s []⎥⎦⎤⎢⎣⎡=3/133/43/43/7329)2(11Et K 1==m r 3==i s []⎥⎦⎤⎢⎣⎡----=3/13/23/21329)2(13Et K 1==m r 2==j s []⎥⎦⎤⎢⎣⎡----=43/23/23/4329)2(12Et K 由子矩阵[]e rs K 合成单元刚度矩阵[]⎥⎥⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎢⎢⎣⎡----------------=3/133/443/23/13/23/43/73/23/43/2143/24003/23/23/403/43/203/13/203/23/103/213/2001329)2(Et K将单元②的单元刚度矩阵补零升阶变为单元贡献矩阵，其在总体刚度矩阵中的位置为：节点号→1 2 3 41 [])2(11K[])2(12K[])2(13K2 [])2(21K[])2(22K[])2(23K3 [])2(31K [])2(32K [])2(33K 4三、计算总体刚度矩阵总体刚度矩阵是由各单元的贡献矩阵迭加而成)2()1(][][][][K K K K e +==∑四、进行节点约束处理根据节点约束情况，在总刚矩阵中可采用划行划列处理约束的方法，由题目易知，节点3和4的已知水平位移和垂直位移都为零，划去其相对应的行和列，则总刚矩阵由8阶变为4阶，矩阵如下：⎪⎪⎭⎪⎪⎬⎫⎪⎪⎩⎪⎪⎨⎧--=⎪⎪⎭⎪⎪⎬⎫⎪⎪⎩⎪⎪⎨⎧⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡--------2/02/03/13043/203/73/23/443/23/133/43/23/43/43/73292211p p v u v u Et329][Et K =1 2 3 413/133/43/43/743/23/23/4----3/13/23/21----000243/23/23/4----3/13003/73/43/403/13/23/21----33/13/23/21----3/43/403/13003/743/23/23/4----40003/13/23/21----43/23/23/4----3/133/43/43/7化简⎪⎪⎭⎪⎪⎬⎫⎪⎪⎩⎪⎪⎨⎧--=⎪⎪⎭⎪⎪⎬⎫⎪⎪⎩⎪⎪⎨⎧⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡--------Et p Et p v u v u 3/1603/160130122072412213424472211 五、求解线性方程组方法：采用LU 分解法 1.求解矩阵[]U 各元素⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡-----⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡------⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡--------75/10775/640075/6475/353007/767/27/7502447~7/877/87/7607/87/337/207/767/27/7502447~13012207241221342447⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡----353/44900075/6475/353007/767/27/7502447~ 得到的[]U 矩阵如下:[]⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡----=353/44900075/6475/353007/767/27/7502447U 2.求解矩阵[]L 各元素⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡--------=⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡----⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡----13012207241221342447353/44900075/6475/353007/767/27/75024471353/6475/767/20175/27/40017/40001 得到的[]L 矩阵如下：[]⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡--------=13012207241221342447L3.进行求解⎪⎪⎭⎪⎪⎬⎫⎪⎪⎩⎪⎪⎨⎧---=⇒⎪⎪⎭⎪⎪⎬⎫⎪⎪⎩⎪⎪⎨⎧--=Et p Et p Et p y Et p Et p Ly 79425/850800225/323/1603/1603/160⎪⎪⎭⎪⎪⎬⎫⎪⎪⎩⎪⎪⎨⎧---=⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡----⇒=⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡Et p Et p Et p v u v u y v u v u U 79425/850800225/323/160353/44900075/6475/353007/7675/27/750244722112211 解得Et p v /422.82-= Et p u /497.12-= Et p v /028.91-= Et p u /897.11=于是求得各节点的位移为：⎩⎨⎧-==Etp v Etp u /028.9/897.111 ⎩⎨⎧-=-=Etp v Etp u /422.8/497.122 ⎩⎨⎧==033v u ⎩⎨⎧==044v u 六、求解相应的支反力（运用静力学的平衡方程进行求解）3号节点和4号节点的支反力如下图所示：。
(完整word版)有限元分析大作业报告要点
有限元分析大作业报告试题1：一、问题描述及数学建模图示无限长刚性地基上的三角形大坝，受齐顶的水压力作用，试用三节点常应变单元和六节点三角形单元对坝体进行有限元分析，并对以下几种计算方案进行比较：（1）分别采用相同单元数目的三节点常应变单元和六节点三角形单元计算；（2）分别采用不同数量的三节点常应变单元计算；（3）当选常应变三角单元时，分别采用不同划分方案计算。
该问题属于平面应变问题，大坝所受的载荷为面载荷，分布情况及方向如图所示。
二、采用相同单元数目的三节点常应变单元和六节点三角形单元计算1、有限元建模（1）设置计算类型：两者因几何条件和载荷条件均满足平面应变问题，故均取Preferences 为Structural（2）选择单元类型：三节点常应变单元选择的类型是Solid Quad 4 node182；六节点三角形单元选择的类型是Solid Quad 8 node183。
因研究的问题为平面应变问题，故对Element behavior（K3）设置为plane strain。
（3）定义材料参数：弹性模量E=2.1e11，泊松比σ=0.3（4）建几何模型：生成特征点；生成坝体截面（5）网格化分：划分网格时，拾取lineAB和lineBC，设定input NDIV 为15；拾取lineAC，设定input NDIV 为20，选择网格划分方式为Tri+Mapped，最后得到600个单元。
（6)模型施加约束：约束采用的是对底面BC 全约束。
大坝所受载荷形式为Pressure ，作用在AB 面上，分析时施加在L AB 上，方向水平向右，载荷大小沿L AB 由小到大均匀分布。
以B 为坐标原点，BA 方向为纵轴y ，则沿着y 方向的受力大小可表示为：}{*980098000)10(Y y g gh P -=-==ρρ2、 计算结果及结果分析 （1） 三节点常应变单元三节点常应变单元的位移分布图三节点常应变单元的应力分布图（2）六节点三角形单元六节点三角形单元的变形分布图六节点三角形单元的应力分布图①最大位移都发生在A点，即大坝顶端，最大应力发生在B点附近，即坝底和水的交界处，且整体应力和位移变化分布趋势相似，符合实际情况；②结果显示三节点和六节点单元分析出来的最大应力值相差较大，原因可能是B点产生了虚假应力，造成了最大应力值的不准确性。
有限元分析大作业
《有限元分析及应用》大作业——齿根弯曲应力计算报告班级：无可奉告姓名：无可奉告学号：无可奉告指导老师：无可奉告目录目录 (2)1.概述 (3)1.1工程问题描述 (3)1.2问题分析 (3)2.建模过程 (4)2.1几何建模 (4)2.2CAE网格划分与计算 (5)2.3后处理 (8)3.多方案比较与结果分析 (9)3.1多方案比较 (9)3.2结果分析 (11)1.概述1.1工程问题描述我在本次作业中的选题为齿根弯曲应力的计算与校核。
通过对机械设计的学习，我们可以知道，齿轮的失效形式主要是齿面接触疲劳和齿根弯曲断裂，而闭式传动硬齿面齿轮的失效形式以齿根弯曲断裂，这个时候进行齿根弯曲应力的校核才比较有意义，在设计问题的时候应当选取这种类型的算例。
设计计算的另一个主要思路是将有限元计算的结果与传统机械设计的结算结果进行对比，以从多方面验证计算结果的准确性。
综上，我们最终选取了《机械原理》（第三版）P50例3-1中的问题进行校核计算。
已知起重机械用的一对闭式直齿圆柱齿轮，传动，输入转速n1=730r/min,输入功率P1=35kW，每天工作16小时，使用寿命5年，齿轮为非对称布置，轴的刚性较大，原动机为电动机，工作机载荷为中等冲击。
z1=29，z2=129，m=2.5mm，b1=48mm，b2=42mm，大、小齿轮均为20CrMnTi，渗碳淬火，齿面硬度为58~62HRC，齿轮精度为7级，试验算齿轮强度。
齿面为硬齿面，传动方式为闭式传动。
根据设计手册查出的许用接触应力为1363.6Mpa，计算结果为1260Mpa，强度合格。
根据设计手册查出的许用弯曲应力为613.3MPa，计算结果为619Mpa，强度略显不够。
1.2问题分析大小齿轮啮合，小齿轮受载荷情况较为严峻，故分析对象应当为小齿轮。
可以看出，由于齿轮单侧受载荷，传动过程中每个齿上载荷的变化过程是相同的，故问题可被简化为反对称问题，仅需研究单个齿。
有限元法理论及应用参考答案(推荐文档)
有限元法理论及应用大作业1、试简要阐述有限元理论分析的基本步骤主要有哪些？答：有限元分析的主要步骤主要有：（1）结构的离散化，即单元的划分；（2）单元分析，包括选择位移模式、根据几何方程建立应变与位移的关系、根据虚功原理建立节点力与节点位移的关系，最后得到单元刚度方程；（3）等效节点载荷计算；（4）整体分析，建立整体刚度方程；（5）引入约束，求解整体平衡方程。
2、有限元网格划分的基本原则是什么？指出图示网格划分中不合理的地方。
题2图答：一般选用三角形或四边形单元，在满足一定精度情况，尽可能少一些单元。
有限元划分网格的基本原则:1.拓扑正确性原则。
即单元间是靠单元顶点、或单元边、或单元面连接2.几何保持原则。
即网络划分后，单元的集合为原结构近似3.特性一致原则。
即材料相同，厚度相同4.单元形状优良原则。
单元边、角相差尽可能小5.密度可控原则。
即在保证一定精度的前提下，网格尽可能的稀疏一些。
（a）(b)中节点没有有效的连接，且（b）中单元边差相差很大。
（c）中没有考虑对称性，单元边差很大。
3、分别指出图示平面结构划分为什么单元？有多少个节点？多少个自由度？题3图答：（a ）划分为杆单元， 8个节点，12个自由度。
（b ）划分为平面梁单元，8个节点，15个自由度。
（c ）平面四节点四边形单元，8个节点，13个自由度。
（d ）平面三角形单元，29个节点，38个自由度。
4、什么是等参数单元？。
答：如果坐标变换和位移插值采用相同的节点，并且单元的形状变换函数与位移插值的形函数一样，则称这种变换为等参变换，这样的单元称为等参单元。
5、在平面三节点三角形单元中，能否选取如下的位移模式，为什么？(1).⎪⎩⎪⎨⎧++=++=26543221),(),(y x y x v yx y x u αααααα (2). ⎪⎩⎪⎨⎧++=++=2652423221),(),(yxy x y x v yxy x y x u αααααα 答：（1）不能，因为位移函数要满足几何各向同性，即单元的位移分布不应与人为选取的 坐标方位有关，即位移函数中的坐标x,y 应该是能够互换的。
西工大作业机考《有限元及程序设计》标准
试卷总分:100 得分:96一、单选题 (共 11 道试题,共 22 分)1.下列关于高精度单元描述正确的是（）。
A.等参元的位移模式和坐标变换采用不同的形函数B.矩形单元形状规则，因而使用范围较广C.6结点三角形单元、10结点三角形单元、8结点矩形单元和12结点矩形单元的单元刚度矩阵的建立过程是不一样的D.6结点三角形单元较容易模拟物体的边界形状正确答案:2.φ=cxy能解决矩形板（）问题。
A.左右均布拉压B.上下均布拉压C.纯剪切D.纯弯曲正确答案:3.下列关于等参元的叙述不正确的是（）。
A.精度较高B.能较好的模拟边界条件C.输入的信息量较少D.输入的信息量较多正确答案:4.薄板的边界不包括（）。
A.简支边界B.固定边界C.自由边界和荷载边界D.非固定边界正确答案:5.下列属于平面应力问题的是（）。
A.平板坝的平板支墩B.挡土墙C.重力水坝D.受内水压力作用的圆管正确答案:6.在应力函数上任意增减一个（），对应力分量无影响。
B.二次项C.三次项D.常数项正确答案:7.下列不属于提高单元精度的方法是（）。
A.增加单元结点数目B.在单元内增设结点C.减少单元结点数目D.设等参元正确答案:8.空间问题的基本平衡微分方程有（）个。
A.2B.3C.4D.5正确答案:9.φ=by2能解决矩形板（）问题。
A.左右均布拉压B.上下均布拉压C.纯剪切D.纯弯曲正确答案:10.下列属于不规则单元的有（）。
A.正四面体单元B.正三棱体单元C.任意四面体单元D.正六面体单元正确答案:11.空间问题的基本未知位移分量有（）个。
A.2B.3C.4D.5二、多选题 (共 16 道试题,共 32 分)1.薄板小挠度弯曲理论的基本假定是（）。
A.直法线假定B.法向位移假定C.中面位移假定D.板内无挤压假定2.弹性力学平面问题按应力求解具体可分为（）两种。
A.逆解法B.顺解法C.半逆解法D.半顺解法3.弹性力学的边界条件有（）。
有限元大作业
一、试题一1、问题描述图示无限长刚性地基上的三角形大坝，受齐顶的水压力，试用三节点常应变单元和六节点三角形单元对坝体进行有限元分析，并对以下集中方案进行比较：1）分别采用相同单元数目的三节点常应变单元和六节点三角形单元计算；2）分别采用不同数量的三节点常应变单元计算；3）当选常应变三角形单元时，分别采用不同划分方案计算。
图 1.12、数学建模及有限元建模2.1数学建模将无限长的三维问题转化为二维问题，取出大坝的一个截面，建立模型如图2.1.1所示，图 2.1.12.2单元选择分别采用相同单元数目的三节点常应变单元和六节点三角形单元计算，分别采用不同数量的三节点常应变单元计算，当选常应变三角形单元时，分别采用不同划分方案计算。
由于是平面应力问题，故Element Behavior K3选择Plane stress。
E=210Gpa,u=0.3。
2.3网格划分按照题意，采用不同数量的三节点常应变单元计算，划分具体的方案见第三部分计算结果。
2.4载荷及边界条件处理边界条件：底端将X，Y方向全约束。
载荷：由于水的压强是随着深度越来越大的P=密度*重力加速度*深度，故对大坝的梯梯形面施加载荷Result = 9800*(10-{Y})3、计算结果及结果分析3.1三节点常应变单元网格划分及结果图 1.3.1.1 图 1.3.1.2图 1.3.1.3 图1.3.1.4 由图可知：DMX=0.143E-4MSMN=49520PASMX=150581PA3.2六节点常应变单元网格划分及结果图1.3.2.1 图 1.3.2.2图 1.3.2.3 图 1.3.2.4由图可知：DMX=0.264E-4MSMN=1893PASMX=306986PA3.3三节点常应变单元网格划分及结果（24个单元）图 1.3.3.1 图 1.3.3.2图 1.3.3.3 图 1.3.3.4由图可知：DMX=0.203E-4MSMN=24617PASMX=199290PA3.4三节点常应变单元网格划分及结果（48个单元）图 1.3.4.1 图 1.3.4.2图 1.3.4.3 图 1.3.4.4由图可知：DMX=0.228E-4MSMN=16409PASMX=243234PA3.5三节点划分方案A图 1.3.5.1 图 1.3.5.2图 1.3.5.3 图 1.3.5.4由图可知：DMX=0.106E-4MSMN=50110PASMX=18194PA3.6三节点划分方案B图 1.3.6.1 图 1.3.6.2图 1.3.6.3 图 1.3.6.44、多方案分析比较4.1相同数目三节点和六节点单元比较DMX(mm) SMN（pa）SMX(pa) 方案数值三节点三角形单元0.143E-4 49520 150581六节点三角形单元0.264E-4 1893 306986表 1.4.1结论：相同的单元划分方案和单元规模，采用不同的阶次的相同形状的单元，分析结果不同。
西安交通大学有限元分析word版第一章
第一章 引言§1-1概述１、有限元方法（The Finite Element Method, FEM ）是计算机问世以后迅速发展起来的一种分析方法。
众所周知，每一种自然现象的背后都有相应的物理规律，对物理规律的描述可以借助相关的定理或定律表现为各种形式的方程（代数、微分、或积分）。
这些方程通常称为控制方程（Governing equation ）。
针对实际的工程问题推导这些方程并不十分困难，然而，要获得问题的解析的数学解却很困难。
人们多采用数值方法给出近似的满足工程精度要求的解答。
有限元方法就是一种应用十分广泛的数值分析方法。
有限元方法是处理连续介质问题的一种普遍方法，离散化是有限元方法的基础。
然而，这种思想自古有之。
齐诺（Zeno 公元前5世纪前后古希腊埃利亚学派哲学家）曾说过：空间是有限的和无限可分的。
故，事物要存在必有大小。
亚里士多德(Aristotle 古希腊大哲学家,科学家)也讲过：连续体由可分的元素组成。
古代人们在计算圆的周长或面积时就采用了离散化的逼近方法：即采用内接多边形和外切多边形从两个不同的方向近似描述圆的周长或面积，当多边形的边数逐步增加时近似值将从这两个方向逼近真解。
图１－２可以用来表示这一过程。
工程中的问题 （力学、物理）各种方程及相应的定解条件（边界条件及初始条件）线性的、边界规则的问题 数值分析法 精确解 近似解 非线性的、边界不规则的问题 解析法 图1-1 工程问题的求解思路图1-2 离散逼近有限单元法 有限差分法图1-3 有限元法与有限差分法比较近代，这一方法首先在航空结构分析中取得了明显的效果：一种称为框架分析法（framework method ）被用来分析平面弹性体（将平面弹性体描述为杆和梁的组合体）（1941，Hrenikoff ）；在采用三角形单元及最小势能原理研究St.Venant 扭转问题时，分片连续函数被用来在子域中近似描述未知函数（1943, Courant ）。
有限元-计算结构力学-大作业
有限元-计算结构力学-大作业本页仅作为文档页封面，使用时可以删除This document is for reference only-rar21year.MarchSHANGHAI JIAO TONG UNIVERSITY 平面应力问题解的Matlab实现姓名: heiya168 学号: 帆哥班级：指导老师：目录1绪论 (4)2平面问题的四节点四边形单元 (4)2.1单元的构造 (4)2.2等参变换 (5)2.3边界条件的处理——置“1”法 (7)3有限元分析流程 (8)3.1程序原理和流程 (8)3.2使用的函数 (9)3.3文件管理 (9)3.4数据文件格式 (10)4算例——开方孔的矩形板拉伸分析 (11)4.1问题的具体参数与载荷 (11)4.2Matlab程序计算 (11)4.3ANSYS建模计算 (13)4.4误差分析 (15)5总结 (15)参考文献 (16)附录 (17)1绪论有限元方法(finite element method)，是求取复杂微分方程近似解的一种非常有效的工具，是现代数字化科技的一种重要基础性原理。
将它用于在科学研究中，可成为探究物质客观规律的先进手段。
将它应用于工程技术中，可成为工程设计和分析的可靠工具。
弹性体在载荷作用下，其基本方程可写成以下的三类方程和两种边界条件。
平衡方程——应力与外载荷的关系；几何方程——应变位移关系；物理方程——应力应变关系；力的边界条件；几何边界条件。
应用最小位能原理，并利用上述关系，最终建立由刚度方程，节点位移和等效节点载荷所构成的求解方程。
带入边界条件求解方程，就可以得出弹性力学问题的一般性解答。
本次大作业基于有限元方法的基本原理，使用Matlab这一平台，针对平面应力问题，采用四节点四边形单元编写了求解单元节点位移的程序。
主要内容包括：1）介绍有限元的基本原理；2）编程基本思路及流程介绍；3）程序原理及说明； 4）具体算例这四个部分。
有限元分析题及大作业题答案
有限元分析及应用作业报告试题10一、问题描述确定图示扳手中的应力, E=210Gpa，μ=0.3, 假设厚度为10mm；并讨论采用何种处理可降低最大应力或改善应力分布。
图1为扳手的基本形状和基本尺寸图二、数学建模与分析由图1及问题描述可知，板手的长宽尺寸远远大于厚度，研究结构为一很薄的等厚度薄板，满足平面应力的几何条件；作用于薄板上的载荷平行于板平面且作用在沿厚度方向均匀分布在办手柄的左边缘线，而在两板面上无外力作用，满足平面应力的载荷条件。
故该问题属于平面应力问题，薄板所受的载荷为面载荷，分布情况及方向如图1所示，建立几何模型，并进行求解。
薄板的材料为钢，则其材料参数：弹性模量E=2.1e11，泊松比σ=0.3三、有限元建模1、单元选择：选取三节点常应变单元来计算分析薄板扳手的位移和应力。
由于此问题为平面应力问题，：三节点常应变单元选择的类型是PLANE42（Quad 4node42），该单元属于是四节点单元类型，在网格划分时可以对节点数目控制使其蜕化为三节点单元。
2、定义材料参数：ANSYS Main Menu: Preprocessor →Material Props →Material Models →Structural →Linear →Elastic →Isotropic →input EX:2.1e11, PRXY:0.3 →OK3、生成几何模型：a.创建关键点点：ANSYS Main Menu: Preprocessor →Modeling →Create →Keypoints →In Active CS →依次输入16个点的坐标→OKb、将这16个关键点有直线依次连起来，成为线性模型4、生成实体模型：ANSYS Main Menu: Preprocessor →Modeling →Create →Areas →Arbitrary →Through KPS →连接特征点→生成两个area→Operate→Subtract→拾取整个扳手区域→OK→生成扳手模型5、结点布置及规模6、网格划分方案ANSYS Main Menu: Preprocessor →Meshing →Mesh Tool →Mesh: Areas, Shape: Tri,Free →Mesh →Pick All (in Picking Menu) →Close( the Mesh Tool window)7、载荷及边界条件处理8、求解控制A、模型施加约束给模型施加x方向约束ANSYS Main Menu: Solution →Define Loads →Apply →Structural →Displacement →On Lines →拾取模型左部的竖直边→OKB、给模型施加载荷ANSYS Main Menu: Solution →Define Loads →Apply →Structural →force→on keypoints→拾取上面左端关键点→700N/mm→okC、分析计算：ANSYS Main Menu: Solution →Solve →Current LS →OK(to close the solve Current Load Step window) →OK6)结果显示：ANSYS Main Menu: General Postproc →Plot Results →Deformed Shape… →select Def + Undeformed →OK (back to Plot Results window) →Contour Plot →Nodal Solu →select: DOF solution →displacement vector sum，von mises stress→OK四、计算结果及结果分析1、三节点常应变单元1）三节点单元的网格划分图2 常应变三节点单元的网格划分平面图图3 常应变三节点单元的网格划分立体图2）三节点单元的约束受载情况图4 常应变三节点单元的约束受载图3)三节点单元的位移分析图5 常应变三节点单元的位移分布图4)三节点单元的应力分析图6 常应变三节点单元的应力分布图2、六节点三角形单元1)六节点三角形单元网格划分图7 六节点三角形单元网格划分图2)六节点三角形单元约束和受载情况分析图8 六节点三角形单元约束受载图3）六节点三角形单元位移分析图9 六节点三角形单元的变形分布图4) 六节点三角形单元的应力分析图9 六节点三角形单元的应力分布图图10 六节点三角形单元的局部应力分布图根据以上位移和应力图，可以得出常应变三节点单元和六节点三角形单元的最小最大位移应力如表1-1所示。
有限元分析与应用大作业
有限元分析及应用大作业课程名称: 有限元分析及应用班级:姓名:试题2：图示薄板左边固定，右边受均布压力P=100Kn/m作用，板厚度为0.3cm；试采用如下方案，对其进行有限元分析，并对结果进行比较。
1)三节点常应变单元；（2个和200个单元）2)四节点矩形单元；（1个和50个单元）3)八节点等参单元。
（1个和20个单元）图2-1 薄板结构及受力图一、建模由图2-1可知，此薄板长和宽分别为2m和1.5m，厚度仅为0.3cm，本题所研究问题为平面应力问题。
经计算，平板右边受均匀载荷P=33.33MPa，而左边被固定，所以要完全约束个方向的自由度，如图2-2所示。
取弹性模量E=2.1×11Pa，泊松比μ=0.3。
P=33.33MPa图2-2 数学模型二、第一问三节点常应变单元（2个和200个单元）三节点单元类型为PLANE42，设置好单元类型后，实常数设置板厚为0.3M。
采用2个单元的网格划分后的结果如图2-3,200个单元的网格划分图如图2-6所示。
约束的施加方式和载荷分布如图2-2中所示。
约束右边线上节点全部自由度。
计算得到的位移云图分别如图2-4、7所示，应力云图如图2-5、8所示。
图2-3 2个三角形单元的网格划分图图2-4 2个三角形单元的位移云图图2-5 2个三角形单元的应力云图图2-6 200个三角形单元的网格划分图图2-7 200个三角形单元的位移云图图2-8 200个三角形单元的应力云图三、第二问四节点矩形单元的计算四节点单元类型为PLANE42，设置好单元类型后，实常数设置板厚为0.3M。
采用1个单元的网格划分后的结果如图2-9,50个单元的网格划分图如图2-12所示。
约束的施加方式和载荷分布如图2-2中所示。
约束右边线上节点全部自由度。
计算得到的位移云图分别如图2-10、11所示，应力云图如图2-13、14所示。
图2-9 1个四边形单元的网格划分图图2-10 1个四边形单元的位移云图图2-11 1个四边形单元的应力云图图2-12 50个四边形单元的网格划分图图2-13 50个四边形单元的位移云图图2-14 50个四边形单元的应力云图四、第三问八节点等参单元的计算四节点单元类型为PLANE82，设置好单元类型后，实常数设置板厚为0.3M。
有限元大作业
有限元大作业第一篇：有限元大作业有限元应力分析报告大作业机械与运载工程学院车辆四班龙恒 20110402415 2014年8月30日一、问题描述桦木板凳材料参数如图形状参数：长40mm，宽30mm，高45mm（其他详细参数见零件图）通过施加垂直于板凳上表面的均匀载荷600N，分析板凳的应变和应力？二、使用inventor进行建模及应力分析1、通过inventor建立板凳3D模型利用草图拉伸等方法建立与零件图中尺寸一致的三维立体板凳模型2、点选环境下的应力分析开始对板凳进行应力分析3、根据所给条件设置材料等参数、将安全系数设为屈服强度，因为板凳主要受压变形点开“木材（桦木）”根据前面所给参数对其进行参数设置4、固定约束如图板凳的4个脚底面设置为固定约束，使得板凳受载后，脚底面不会沿垂直方向位移，模拟真实情况5、施加载荷在板凳上表面施加大小为600N的垂直均布载荷（这里是模拟一个成人坐上去的重力）6、划分网格通过设置网格的尺寸参数来划分出5种不同网格数量，从而得出5种不同网格数划分得出的应力应变分布图，最后分析划分不同网格数对结果的影响。
（1）网格最大（2）网格较大（3）网格一般大小（4）网格较小（5）网格最小7、求解得出结果得出5组不同网格数所得数据（应力云图，应变云图，所有结果数据）（1）网格数1437根据应力云图可知，红色地方所受的应力最大，最大应力为：15.48Mpa 根据应变云图可知，红色地方的应变最大，最大应变为：0.001434μl（2）网格数8651根据应力云图可知，红色地方所受的应力最大，最大应力为：18.88Mpa 根据应变云图可知，红色地方的应变最大，最大应变为：0.001755μl（3）网格数20484根据应力云图可知，红色地方所受的应力最大，最大应力为：22.62Mpa 根据应变云图可知，红色地方的应变最大，最大应变为：0.002103μl（4）网格数41578根据应力云图可知，红色地方所受的应力最大，最大应力为：23.76Mpa 根据应变云图可知，红色地方的应变最大，最大应变为：0.002206μl（5）网格数68788根据应力云图可知，红色地方所受的应力最大，最大应力为：25.97Mpa 根据应变云图可知，红色地方的应变最大，最大应变为：0.002454μl综合上述5种请况可知随着网格的细分，所得的应变以及应力的结果是收敛的。
西南交大《有限元方法》课件—赵华主讲习题全
7. 8. 9.
u a1 a2 x a3 y a4 xy a5 x 2 v a6 a7 x a8 y a9 xy a10 y 2
10. 采用面积坐标导出六节点三角形单元的应变矩阵。
Y 4 6 3 b 1 5 a 2 X
第3题
第4题
Y
c
两节点梁单元 过度单元 四节点单元
Y 1 (1) (3)
2
(2)
3
习题5附图
X
习题6附图
第2章习题
1. 单元形函数有何物理意义、特性？在单元分析中的作用是什么？合理 的形函数须满足什么条件？为什么？ 2. 在平面三节点三角形单元中，能否取如下的位移模式，为什么？ u a1 x 2 a2 xy a3 y 2 u a1 a2 x 2 a3 y 2 2 2 v a x a xy a y v a a x a y 或 ， 4 5 6 4 5 6 3. 三角形单元ABC（尺寸如图），已知材料的弹性模量 E 2 *105 MPa， 5 泊松比 0.3，如A点沿X的位移 u A 2 *10 mm ，B点沿Y的位移 5 vB 10 mm ，而C点位移为零。计算单元内应力以及节点力。 4. 对图示六节点矩形单元，应取什么样的形状函数来表示位移模式？ 试写出并检验是否满足收敛性条件。 5. 图为局部系 , 表示的8节点矩形单元。写出形函数表示的位移模 式并证明满足收敛性条件。
第四章习题
1、论述Ritz法与有限元法的相同点和主要差异。 2、利用最小势能原理推导空间四面体单元的刚度方程。
3、导出相应于下列泛函驻值条件的Euler方程，
2 2 2 1 2C d 2 x y z
有限元分析报告与应用大作业
有限元分析及应用大作业课程名称: 有限元分析及应用班级:姓名:试题2：图示薄板左边固定，右边受均布压力P=100Kn/m作用，板厚度为0.3cm；试采用如下方案，对其进行有限元分析，并对结果进行比较。
1)三节点常应变单元；（2个和200个单元）2)四节点矩形单元；（1个和50个单元）3)八节点等参单元。
（1个和20个单元）图2-1 薄板结构及受力图一、建模由图2-1可知，此薄板长和宽分别为2m和1.5m，厚度仅为0.3cm，本题所研究问题为平面应力问题。
经计算，平板右边受均匀载荷P=33.33MPa，而左边被固定，所以要完全约束个方向的自由度，如图2-2所示。
取弹性模量E=2.1×11Pa，泊松比μ=0.3。
P=33.33MPa图2-2 数学模型二、第一问三节点常应变单元（2个和200个单元）三节点单元类型为PLANE42，设置好单元类型后，实常数设置板厚为0.3M。
采用2个单元的网格划分后的结果如图2-3,200个单元的网格划分图如图2-6所示。
约束的施加方式和载荷分布如图2-2中所示。
约束右边线上节点全部自由度。
计算得到的位移云图分别如图2-4、7所示，应力云图如图2-5、8所示。
图2-3 2个三角形单元的网格划分图图2-4 2个三角形单元的位移云图图2-5 2个三角形单元的应力云图图2-6 200个三角形单元的网格划分图图2-7 200个三角形单元的位移云图图2-8 200个三角形单元的应力云图三、第二问四节点矩形单元的计算四节点单元类型为PLANE42，设置好单元类型后，实常数设置板厚为0.3M。
采用1个单元的网格划分后的结果如图2-9,50个单元的网格划分图如图2-12所示。
约束的施加方式和载荷分布如图2-2中所示。
约束右边线上节点全部自由度。
计算得到的位移云图分别如图2-10、11所示，应力云图如图2-13、14所示。
图2-9 1个四边形单元的网格划分图图2-10 1个四边形单元的位移云图图2-11 1个四边形单元的应力云图图2-12 50个四边形单元的网格划分图图2-13 50个四边形单元的位移云图图2-14 50个四边形单元的应力云图四、第三问八节点等参单元的计算四节点单元类型为PLANE82，设置好单元类型后，实常数设置板厚为0.3M。
【作业参考】[西安交通大学]《有限元及工程应用》在线作业-000015
![【作业参考】[西安交通大学]《有限元及工程应用》在线作业-000015](https://img.taocdn.com/s3/m/c11af31b6d175f0e7cd184254b35eefdc8d315f0.png)
西交《有限元及工程应用》在线作业-00001
等参数(Isoparametric)变换能够实行,必须使变换行列式(雅克比(Jacobi)行列式)的值在整个单元上( )。
A:等于0
B:为负数
C:为正数
D:不等于0
答案：D
下述那个不是单元刚度矩阵的特点:( )。
A:连续性
B:对称性
C:奇异性
D:按节点分块
答案：A
等参数单元划分时单元形状应该尽力接近( )。
A:长方形
B:菱形
C:正方形
D:梯形
答案：C
若某函数的假设解表示为,R表示余量,加权余量法的表达式为,其中表示求解域,W表示权函数,那么选取W为以下哪种权函数时的加权余量法被称为Galerkin法( )。
A:
B:
C:
D:
答案：B
平面应力问题与平面应变问题都具有( )个独立的应力分量。
A:1
B:2
C:3
D:4
答案：C
下面关于杆单元公共节点处位移的变形协调条件描述正确的说法是( )。
A:相邻单元公共节点处,位移叠加,力相同
B:相邻单元公共节点处,位移叠加,力叠加
C:相邻单元公共节点处,位移相同,力叠加
D:相邻单元公共节点处,位移相同,力相同
答案：C
对于常见的梯子那样的框架结构进行模型化处理应使用( )。
A:梁单元
B:板单元
C:壳单元
D:其它单元
答案：A
1。
有限元原理与应用小作品
有限元原理与应用小作品CRH2动车组拖车转向架车轮CAD/CAE设计班级：车辆工程073班（教学1班）姓名：杨文卓学号： 200703837摘要:在我国，铁路是国家的重要基础设施、国民经济的大动脉和大众化交通工具，在交通运输行业中起着主导性作用。
随着时速200km/h高速动车组的引进，如何评价引进动车组的安全性、可靠性等成为十分紧迫和重要的问题。
对于高速动车组而言，转向架的重要性更是不言而喻。
转向架是铁路车辆的基本零件之一，其性能直接关系到铁路车辆的安全性、稳定性和可靠性。
论文以CRH2动车组车轮为研究对象，并建立其有限元模型。
按照UIC标准，应用Solidworks对其进行CAD/CAE设计分析。
表1 CRH2动车组拖车转向架基本参数项目规格设计最高速度/km/h 250营业最高速度/km/h 200额定轴重/kN 137.2（14t）满员时最大轴重(200%定员/kN) 156.8（16t）编组能通过的最小曲线半径/m 180转向架转角/度 4.0轴距/mm 2500车轮直径/mm 860关键词:solidworks 轮对 CAD/CAE车轮的CAD/CAE分析一、问题描述车轮是轮对的重要组成部分，其疲劳强度直接关系到动车组运行的安全性、可靠性、稳定性等，故动车组车轮需进行车轮静强度、动强度和轮轴过盈配合强度等3个方面的分析。
本处只做静强度分析：UIC 510—5/2003《整体车轮技术条件》；EN 13979—l/2003《铁路应用轮对和转向架车轮技术验收程序》第一部分：锻制和轧制车轮。
二、CAD设计（CAD零件工程图）车轮的建模过程（零件图如图1）：拉伸轮胚→切左腹板→切右腹板→旋踏面→打轮毂孔→倒角→完成图1 CRH2拖车车轮零件工程图车轮的造型图如图2：图2 CRH2 拖车车轮造型图三、CAE设计1.受力分析根据UIC 510-5：2003（整体车轮技术）标准进行车轮设计，对于安装到动轴上的车轮，考虑车轮通过直线、曲线和道岔时的载荷。
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
有限元原理及工程应用——大作业学院：机械工程学院班级：硕4002班小组成员：李追3114001089陈草31140010802015.5.19作业题目：利用有限元方法对简支梁问题进行求解，梁的横截面为矩形，其约束情况如图1所示。
已知梁的几何尺寸和物理参数如下：（1）几何尺寸：长度40cm L =，截面尺寸2cm 02cm .b t ⨯=⨯；（2）物理参数：弹性模量70E =GPa ，泊松比0.3ν=，密度3=2700kg/m ρ。
图1．梁及其横截面示意图要求：(1) 至少划分五个节点（四个单元）； (2) 给出单元节点信息；(3) 给出单元刚度矩阵和质量矩阵； (4) 给出总刚度矩阵和总质量矩阵； (5) 求出梁各界固有频率及振型（五阶）；(6) 将所得结果与理论值进行对比，验证方法的可行性。
解：由有限元知识，根据Rayleigh-Ritz 法，解有限元分为四步：建立离散化、单元分析、形成总体方程、解方程，具体步骤如下： （1）建立离散化这里我们将矩形截面简支梁等分四等分，即分为六节点的五个杆单元，如图2所示：每个单元尺寸40cm=8cm 55L l ==，这里只考虑杆在竖直平面的弯曲，每个节点只有y 方向位移和绕z 轴的旋转自由度。
（2）单元分析构造一组Lagrange 插值基函数，在本节点值为1，其他节点值为0。
从Rayleigh-Ritz 法可以看到，插值函数要p 次可微，最高阶导数出现在应变能表达式中；同样，我们可以这一原则适用于基函数的选择以及形状函数，否则我们将无法正确计算应变能当我们使用有限元逼近方法。
梁的弯曲问题，应变能计算公式：()2220,12Lz v x t U EI dx x ⎛⎫∂= ⎪∂⎝⎭⎰ （1-1） 其中，E 为弹性模量，I z 为截面惯性矩。
从公式可知，位移函数必须连续，并且二阶导数平方可积。
如图3，是一维杆单元模型，每个节点两个自由度，该单元含有四个自由度，即（,,,i zi j zj v v θθ）。
本题中我们采用三次多项式插值函数：()231234u x x x x αααα=+++ （1-2）因此，我们必须给出四个形函数（位移模式）。
图3 一维杆单元模型 1) 构造Hermite 插值函数。
选择局部坐标系x ξ=(()0,0i n ,(),0j n l )，其中l 是单元长度，转角z θ是挠度值v 的一阶导数，定义边界条件：()()''i i i iv x v v x v== （1-3）因此，我们给出变形的Hermite 的多项式插值函数：()()()22(0)(1)'11ii i i i i v Hv H v ξξξ===+∑∑ （1-4）其中，()(0)i H ξ和()(1)i H ξ分别满足如下条件，对应的图形如图4所示：()()()()(0)(0)(1)(1)d 0d d 0d jji ijijii j ijH HHH ξξξξδξξξδξ==== （1-5）图4 一维Hermite 插值多项式2) 基于Langrage 和Hermite 插值多项式，写出单元形函数1.00.5()()()()()()()()()()(0)23(1)231121(0)3(1)323242==132==2==32==N H N H l N HN Hlξξξξξξξξξξξξξξξξξ-+-+-- （1-6）节点位移值也可以得出()()()()()1234i zi j zj v N v N N v N ξξξθξξθ=+++ （1-7）同时，表达式（1-7）用矩阵表示为()(){}e v ξξ=⎡⎤⎣⎦N v （1-8）其中，{}{}T,,,e i zi j zj v v θθ=v ，()()()()()1234N N N N ξξξξξ=⎡⎤⎡⎤⎣⎦⎣⎦Ν3) 用能量表达式替代表达式中的v 和ϕ。
动能表达式：(){}[]{}1201122Te e T Av d ρξαξ==⎰e v M v （1-9）将（1-8）带入（1-9），得到{}()(){}{}1012T T e eT A d ρξξξ=⎡⎤⎡⎤⎣⎦⎣⎦⎰v N N v （1-10）从而获得质量矩阵：[]()()10Te A d ρξξξ=⎡⎤⎡⎤⎣⎦⎣⎦⎰M N N （1-11）带入x ξ=，[]2222156225413224133541315622420133224e ll l ll l Al l l l ll l ρ-⎡⎤⎢⎥-⎢⎥=⎢⎥-⎢⎥---⎣⎦M （1-12） 应变能表达式：(){}[]{}{}()(){}221201''''3011221=2e v K v v N N v ⎛⎫∂== ⎪∂⎝⎭⎧⎫⎡⎤⎡⎤⎨⎬⎣⎦⎣⎦⎩⎭⎰⎰Tz e e T T z e e v U EI d EI d l ξξξξξξ （1-13）刚度矩阵表达式：[]()()1''''30TzEId lξξξ⎡⎤⎡⎤=⎣⎦⎣⎦⎰e K N N （1-14）带入x ξ=，可以得到[]2232212612664621261266264z l l l l l l EI l l l l ll l -⎡⎤⎢⎥-⎢⎥=⎢⎥---⎢⎥-⎣⎦e K （1-15） （3） 形成总体方程将每个杆单元的能量方程组装。
完整梁上的的总动能和能量的和所做的总功梁上的外力作用，所有的自由度的位移矢量可以给出：{}{}112233445566Tz z z z z z v v v v v v θθθθθθ=v （1-16）将位移矢量转换为全局坐标系下的位移矢量，变换矩阵为：{}[]{}e e v v =a （1-17）[]1100000000000010000000000001000000000000100000000⎡⎤⎢⎥⎢⎥=⎢⎥⎢⎥⎣⎦a []2001000000000000100000000000010000000000001000000⎡⎤⎢⎥⎢⎥=⎢⎥⎢⎥⎣⎦a []3000010000000000001000000000000100000000000010000⎡⎤⎢⎥⎢⎥=⎢⎥⎢⎥⎣⎦a （1-18） []4000000100000000000010000000000001000000000000100⎡⎤⎢⎥⎢⎥=⎢⎥⎢⎥⎣⎦a []5000000001000000000000100000000000010000000000001⎡⎤⎢⎥⎢⎥=⎢⎥⎢⎥⎣⎦a 分别以矩阵形式给出动能和质量矩阵：{}[][][]{}{}[]{}511122T T Ttot e e e e T ===∑v a M a v v M v （1-19）[]222222222221562254130000022413300000541331205413000133081330000054133120541300013308133000054133120540000133081300000054133120000001330000000005420000M --------------=l l l l l l l l l l l l l l l l l l l l l A l l l l l l l l ρ222220000000000000000001300300054138133413156220133224⎡⎤⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎣⎦--------l l l l l l ll l l ll [][][][]51Te e e e ==∑M a M a （1-20）因此，总体质量矩阵为总体刚度矩阵：求固有频率。
总应变能{}[][][]{}{}[]{}411122T T Ttote e e e U ===∑v a K a v v K v (1-21) [][][][]51Te e e e ==∑K a K a (1-22)[]3222222222221261260000000064620000000012624012600000062086200000000126240126000000620862000000001262401260000006208600000012600000062000000000000K ---------------=-z EI l l l l l ll l l l l l l l l l l l l l l ll l l l l ll l 2222220024012608621261266264⎡⎤⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢-----⎣⎦-⎥l l l l l l l l l l l利用Lagrange 方程，推导简支梁自由振动方程：[]{}[]{}0+=M u K u (1-23)这里，我们假设简支梁做简谐振动，则(){}j tt e ω=u A (1-24)因此，特征方程为：[][](){}0λ-=K M A (1-25)其中，2λω=，ω为固有频率。
（4） 解方程，有限元分析结果基于上述理论,我们获得了采用MATLAB 的有限元分析程序代码。
提交边界条件、材料特性和几何参数到上面的方程, 使用MATLAB 代码我们得到以下结果。
（1）离散简支梁为5单元6节点,那么我们得到的单元质量矩阵和刚度矩阵,如下所示：（2）总质量矩阵和总刚度矩阵如下：[]0.00320.00000.0011-0.0000000000000.00000.00000.0000-0.00000000M =00000.00110.00000.006400.0011-0.0000000000-0.0000-0.000000.00000.0000-0.0000 000000000.00110.00000.006400.0011-0.0000000000-0.0000-0.000000.00000.0000-0.0000000000000.00110.00000.006400.0011-0.000000000-0.0000-0.000000.00000.0000-0.0000000000000.00110.00000.0064 00.0011-0.0000000000-0.0000-0.000000.00000.0000-0.000000000 0000.00110.00000.0032-0.000000000000-0.0000-0.0000-0.00000.0000⎡⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎣⎤⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎢⎥⎦[] 0.0032 0.0000 0.0011 -0.0000 0.0000 0.0000 0.0000 -0.0000e 0.0011 0.0000 0.0032 -0.0000-0.0000 -0.0000 -0.0000 0.0000M ⎡⎤⎢⎥⎢⎥=⎢⎥⎢⎥⎣⎦[]4 2.1875 0.0875 -2.1875 0.0875 0.0875 0.0047 -0.0875 0.00231.010-2.1875 -0.0875 2.1875 -0.0875 0.0875 0.0023 -0.0875 0.0047Ke ⎡⎤⎢⎥⎢⎥=⨯⎢⎥⎢⎥⎣⎦[]4 2.1875 0.0875 -2.1875 0.0875 0 0 0 0 0 0 0 0 0.0875 0.0047 -0.08751.010 K=⨯0.0023 0 0 0 0 0 0 0 0-2.1875 -0.0875 4.3750 0 -2.1875 0.0875 0 0 0 0 0 00.0875 0.0023 0 0.0093 -0.0875 0.0023 0 0 0 0 0 00 0 -2.1875 -0.0875 4.3750 0 -2.1875 0.0875 0 0 0 00 0 0.0875 0.0023 0 0.0093 -0.0875 0.0023 0 0 0 00 0 0 0 -2.1875 -0.0875 4.3750 0 -2.1875 0.0875 0 00 0 0 0 0.0875 0.0023 0 0.0093 -0.0875 0.0023 0 00 0 0 0 0 0 -2.1875 -0.0875 4.3750 0 -2.1875 0.0875 0 0 0 0 0 0 0.0875 0.0023 0 0.0093 -0.0875 0.0023 0 0 0 0 0 0 0 0 -2.1875 -0.0875 2.1875 -0.0875 0 0 0 0 0 0 0 0 0.0875 0.0023 -0.0875 0.0047⎡⎤⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎣⎦(3)简支梁的振动分析表1列出了简支梁的五阶固有振动频率。