下载文档原格式
弹簧单元与梁单元实例计算1.绪论有限元法也叫有限单元法（finite element method, FEM），是随着电子计算机的发展而迅速发展起来的一种弹性力学问题的数值求解方法。
五十年代初，它首先应用于连续体力学领域—飞机结构静、动态特性分析中，用以求得结构的变形、应力、固有频率以及振型。
由于这种方法的有效性，有限单元法的应用已从线性问题扩展到非线性问题，分析的对象从弹性材料扩展到塑性、粘弹性、粘塑性和复合材料，从连续体扩展到非连续体。
关键词：有限元方法，数值求解，动态分析2.有限元方法2.1有限元法概述有限元法是把要分析的连续体假想地分割成有限个单元所组成的组合体，简称离散化。
这些单元仅在顶角处相互联接，称这些联接点为结点。
离散化的组合体与真实弹性体的区别在于：组合体中单元与单元之间的联接除了结点之外再无任何关联。
但是这种联接要满足变形协调条件，即不能出现裂缝，也不允许发生重叠。
显然，单元之间只能通过结点来传递内力。
通过结点来传递的内力称为结点力，作用在结点上的荷载称为结点荷载。
当连续体受到外力作用发生变形时，组成它的各个单元也将发生变形，因而各个结点要产生不同程度的位移，这种位移称为结点位移。
在有限元中，常以结点位移作为基本未知量。
并对每个单元根据分块近似的思想，假设一个简单的函数近似地表示单元内位移的分布规律，再利用力学理论中的变分原理或其他方法，建立结点力与位移之间的力学特性关系，得到一组以结点位移为未知量的代数方程，从而求解结点的位移分量。
然后利用插值函数确定单元集合体上的场函数。
显然，如果单元满足问题的收敛性要求，那么随着缩小单元的尺寸，增加求解区域内单元的数目，解的近似程度将不断改进，近似解最终将收敛于精确解。
2.2有限元法的优点1、物理概念浅显清晰，易于掌握。
有限元法不仅可以通过非常直观的物理解释来被掌握，而且可以通过数学理论严谨的分析掌握方法的本质。
2、描述简单，利于推广。
“有限元法基础及应用”补充讲义一、弹簧单元与弹簧系统1、 弹簧单元分析 1）单元描述弹性系统受力平衡时，从中隔离出一个典型弹簧单元进行分析。
单元节点编号：j i ,节点位移（基本未知量）：j i u u ,单元节点力（单元在节点处受到的作用力）：j i f f ,已知弹簧的物理特性：∆⋅=k F其中：为弹簧力（拉伸为正）为弹簧伸长量，为弹簧刚度F u u k i j -=∆, 2）建立弹簧单元的有限元特性方程考虑弹簧元在系统中变形平衡时的条件：力平衡条件和弹簧物理特性，得到下列方程：ji i j j j i i j i ku ku u u k F f ku ku u u k F f +-=-==-=--=-=)()(写成矩阵形式：⎭⎬⎫⎩⎨⎧⎥⎦⎤⎢⎣⎡--=⎭⎬⎫⎩⎨⎧j i j i u u k k k k f f 上式的矩阵符号形式为：kd f =方程（1-2）或（1-3）称为弹簧单元的刚度方程，反映了单元特性，即节点力~节点位移之间的关系。
式（1-3）中：称为单元节点力列阵称为单元节点位移列阵称为单元刚度矩阵,,,⎭⎬⎫⎩⎨⎧=⎭⎬⎫⎩⎨⎧=⎥⎦⎤⎢⎣⎡--=j i j i f f u u k k k k f d k（1-1）（1-2）（1-3）3）弹簧单元刚度方程的讨论a. 有何特点？k对称、奇异、主对角元素恒正。
b . 么？中元素的物理意义是什k刚度矩阵元素的大小等于弹簧刚度。
从对方程（1-2）分析的分析可以看出，矩阵中某列的各元素代表列序号对应节点有单位位移，其它节点位移为零时，单元各节点上的节点力；某行的各元素分别是单元各节点的位移对行序号对应节点的节点力贡献系数。
因此，矩阵中任意一个元素ij k 的物理意义是：j 节点的位移对i 节点的节点力贡献系数，或者j 节点有单位位移，其他节点位移为零时，i 节点上的节点力。
c. 单元刚度方程可以求解吗？为什么？不可以。
单元刚度方程仅仅表征一个单元的力学特性，单元水平上无法确定单元节点位移。
“有限元法基础及应用”补充讲义一、弹簧单元与弹簧系统1、 弹簧单元分析 1）单元描述弹性系统受力平衡时，从中隔离出一个典型弹簧单元进行分析。
单元节点编号：j i ,节点位移（基本未知量）：j i u u ,单元节点力（单元在节点处受到的作用力）：j i f f ,已知弹簧的物理特性：∆⋅=k F其中：为弹簧力（拉伸为正）为弹簧伸长量，为弹簧刚度F u u k i j -=∆, 2）建立弹簧单元的有限元特性方程考虑弹簧元在系统中变形平衡时的条件：力平衡条件和弹簧物理特性，得到下列方程：ji i j j j i i j i ku ku u u k F f ku ku u u k F f +-=-==-=--=-=)()(写成矩阵形式：⎭⎬⎫⎩⎨⎧⎥⎦⎤⎢⎣⎡--=⎭⎬⎫⎩⎨⎧j i j i u u k k k k f f 上式的矩阵符号形式为：kd f =方程（1-2）或（1-3）称为弹簧单元的刚度方程，反映了单元特性，即节点力~节点位移之间的关系。
式（1-3）中：称为单元节点力列阵称为单元节点位移列阵称为单元刚度矩阵,,,⎭⎬⎫⎩⎨⎧=⎭⎬⎫⎩⎨⎧=⎥⎦⎤⎢⎣⎡--=j i j i f f u u k k k k f d k（1-1）（1-2）（1-3）3）弹簧单元刚度方程的讨论a. 有何特点？k对称、奇异、主对角元素恒正。
b . 么？中元素的物理意义是什k刚度矩阵元素的大小等于弹簧刚度。
从对方程（1-2）分析的分析可以看出，矩阵中某列的各元素代表列序号对应节点有单位位移，其它节点位移为零时，单元各节点上的节点力；某行的各元素分别是单元各节点的位移对行序号对应节点的节点力贡献系数。
因此，矩阵中任意一个元素ij k 的物理意义是：j 节点的位移对i 节点的节点力贡献系数，或者j 节点有单位位移，其他节点位移为零时，i 节点上的节点力。
c. 单元刚度方程可以求解吗？为什么？不可以。
单元刚度方程仅仅表征一个单元的力学特性，单元水平上无法确定单元节点位移。
《有限元分析基础教程》（曾攀）笔记⼆-梁单元有限元⽅程推导不得不说，Mathematica 真是个好东西，以前学习有限元的时候，对于书中的⽅程推导，看到了就看过去了，从没有想过要⾃⼰推导⼀遍，原因是⼿⼯推导太复杂。
有了MM ，原来很复杂的东西突然变得简单了。
1.单元⼏何描述上图是纯弯梁单元，长度l ，弹模E ，⾯积A ，惯性矩I 。
两个节点1和2的位移列阵为q e =[v 1,θ1,v 2,θ2]Tv 是挠度(defection)，或者叫位移；θ是转⾓(slope)。
需注意的是v 和θ的⽅向，⼀个是向上，⼀个是逆时针。
两个节点的节点⼒矩阵为P e =[P v 1,M 1,P v 2,M 2]T当然实际情况往往是在梁的长度⽅向上作⽤有荷载，⽽不是只在节点处有，这时就要进⾏荷载等效，后⾯会有说明。
注意这两个矩阵都是列矩阵。
需要注意的是，节点⼒矩阵表⽰的的是节点上的所有的⼒，不仅包括荷载引起的等效节点⼒，还包括节点的反⼒，反⼒矩等。
2.单元位移场表达由于有4个位移节点的已知条件，那么假设纯弯曲梁单元的位移挠度函数具有四个待定系数，如下形式v (x )=a 0+a 1x +a 2x 2+a 3x 3对于两端节点，位移和转⾓分别为v 1,θ1,v 2,θ2，注意挠曲线⽅程在⼀点出的导数值即为改点的转⾓，所以四个边界条件为v (0)=v 1v ′(0)=θ1v (L )=v 2v ′(L )=θ2使⽤MM 求解⽅程组将求得的待定系数带⼊原⽅程，可得将四个位移合并同类项，可以得到即最终的挠曲线⽅程vfea 为 vfea =θ1x 3L 2−2x 2L +x +θ2x 3L 2−x 2L +v12x 3L 3−3x 2L 2+1+v23x 2L 2−2x 3L 3如果令ζ=x L ，上式中位移前的系数组成的矩阵称之为形函数矩阵，也就是常说的形函数。
即v (x )=N (x )q e 3.单元应变场，应⼒场的表达应变的表达式为ε=−yv ″其中B(x)=-yN''(x)，B(x)叫做单元的⼏何矩阵，表⽰应变与位移的⼏何关系。
基于有限元分析软件的弹簧、质量、阻尼振动系统的瞬态动力分析本文对振动系统瞬态动力学分析方法进行了阐述。
以有限元分析软件ANSYS 10.0作为平台，对弹簧、质量、阻尼系统进行瞬态动力学求导与分析，详细论述了分析的过程，结果与理论分析吻合得很好。
本文的研究可以为制造业的信息化过程提供一定的参考。
0 振动力学简介振动是一种运动形态，是指物体在平衡位置附近作往复运动。
从广义上讲，如果表征一种运动的物理量作时而增大时而减小的反复变化，就可以称这种运动为振动。
如果变化的物理量是一些机械量或力学量，例如物体的位移、速度，加速度、应力及应变等等，这种振动便称为机械振动。
振动力学是指借助数学、物理、实验和计算技术，探讨各种振动现象，阐明振动的基本规律，以便克服振动的消极因素，利用其积极因素，为合理解决各种振动问题提供理论依据的一门科学。
振动是普遍存在的物理现象，是受外界激励而使系统包含的质量、弹性、阻尼等元件对外界激励的响应。
在所有科学领域和日常生活中都会遇到各种不同程度的振动，基于振动对工业生产的重要影响，国内外许多学者在此领域进行了大量的研究。
在机械结构的动力学特性研究上主要体现在以下几方面：(1) 建立振动模型；(2) 确定结构系统的动态特性；(3) 采用非比例阻尼方法准确估计系统的阻尼矩阵；(4) 基于实验数据结构的有限元模型修正等方面。
1 振动系统瞬态动力学分析方法图1 振动模型关系图一般振动问题是由振动系统、激励和响应三部分组成，三者间的关系可表示为如图1所示。
振动问题的研究对象即为振动系统，外界激振力等因素叫做激励(输入)，作用于系统使之产生振动响应(输出)。
振动问题就是从以上三者中，已知两个量来求解另一个参数。
瞬态动力学分析（也称时间历程分析）是用于确定承受任意的随时间变化载荷的结构的动力学响应的一种方法。
可以用瞬态动力学分析确定结构在静载荷，瞬态载荷和简谐载荷的任意组合下的随时间变化的位移、应变、应力及力。
有限元方法对工程问题的研究具有重要帮助，本文介绍了有限元方法及其基本原理以及优势，对其解决问题的一般步骤进行了总结分析。
并具体用 MATLAB 对弹簧单元和梁单元实例进一步分析计算，得出了弹簧单元各个节点位移、支反力以及单元内力；计算了梁单元所受位移、转角、力和力矩大小并绘制其剪力图和转矩图。
有限元方法；弹簧单元；梁单元;matlab1有限单元法最初作为结构力学位移法的拓展,它的基本思路是将复杂的结构或物体看成由有限数目的单元体彼此仅在结点处相联系而构成的集合体。
首先,对每一个单元分析其特性,建立相关物理量之间的相互联系;然后，依据单元之间在结点处的联系，再将各单元组装成整体,从而获得整体特性方程;最后,应用与所得方程相应的解法，即可完成整个问题的分析。
这种先“化整为零”，然后再“集零为整"和“化未知为已知”的研究方法，是有普遍意义的，是科学研究的基本方法之- -。
有限单元法作为一种近似的(除杆件体系结构静力分析外)数值分析方法,它借助于矩阵等数学工具,尽管计算工作量很大,但是分析流程是- -致的,有很强的规律性和统--的模式，因此特别适合于编制计算机程序来处理。
-一般来说,一定前提条件下分析的近似性,随着离散化单元网格的不断细化,计算精度将随之提高。
随着计算机软硬件技术的飞速发展,有限单元分析技术得到了越来越多的应用，五十多年来的发展几乎涉及了各类科学、工程领域中的问题。
从应用的深度和广度来看,有限单元法的研究和应用正继续不断地向前探索和推进。
2基础工程学科中的各种力学问题,最终归结为求解数学物理方程边值或初值问题。
而传统的研究成果只能对较简单、规则的问题进行解析求解,大量的实际科学、工程计算问题,由于数学上的困难无法得到解析的解答。
从有限单元法正式提出至今已经历了半个多世纪的发展,用有限单元法来解决问题,从理论上讲，无论是简单的十维杆件体系结构,还是承受复杂荷载和不规则边界情况的二维平面轴对称问题、三维空间块体等问题的静力、动力和稳定性分析,考虑材料具有非线性力学行为和有限变形的分析,温度场、电磁场,流体、液-固、结构与土壤相互作用等工程复杂问题的分析,都可以得到满意的解决。
作为教材,本书将只介绍应用最广泛的、以结点位移作为基本未知量的“位移型有限元”，简称“位移元”,讨论范围仅限于最基本的线弹性问题。
虽然有限单元法能解答各种各样的复杂工程问题,但其求解的思路和过程是基本相同的,大致分为以下4个步骤:(1)结构离散化及数据化。
用--些假想的线或面对连续的结构进行分割，使其成为由结点连接起来的有限数量的单元体(注意单元体和材料力学中的微元体是根本不同的,它的尺度是有限值而不是微量)集合体,用这个单元集合体代替原来要分析的物体。
再通过某种编码格式，将单元集合体的所有相关信息编码成为数据，以供计算机分析。
(2) 单元分析。
这一步工作将对结构中的任- -典型单元进行分析。
以平面静力问题为例,单元分析的目的是将单元内任一点的位移、应变、应力都表示为单元结点位移的函数,并最终建立结点力与结点位移之间的关系方程，即单元刚度方程。
(3)整体分析。
有了单元特性分析的结果，像结构力学中解超静定结构的位移法一样,可以建立起表示整个结构(实质上，更确切地说是单元集合体)结点平衡的方程组,即整体刚度方程。
本步骤计算的细节取决于所求解的问题和所编制计算机程序的处理方法,对于一些问题将存在坐标(局部与整体)转换问题、位移边界条件的引人等。
(4)解方程组和输出计算结果。
对线弹性问题,整体刚度方程一.般是- -组高阶的线性代数方程组。
由于整体刚度矩阵具有带状、稀疏和对称等特性，在有限单元法发展过程中,人们通过研究，建立了许多不同的存储方式和计算方法，目的是节省计算机的存储空间和提高计算效率。
利用相应的计算方法即可求出全部未知的结点位移。
求出结构全部结点位移后,利用分析过程中已建立的一-些关系,即可进--步计算单元中的应力或内力,并以数表或图形的方式输出计算结果。
2.31.可以很容易地模拟不规则形状的结构。
2.可以毫无困难地处理一般的荷载条件。
3.因为单元方程是单个地建立的，因此可以模拟由几种不同材料构成的物体。
4.可以处理数量不受限制的和各种类型的边界条件。
5.单元的尺寸大小可以变化，必要时可使用小单元。
6.改变有限元模型比较容易，花费不大。
3该步骤包括将域分解成单元和节点。
对于像桁架和刚架这类离散系统，已经离散化，这一步就不需要了。
此处获得的结果应该已经是精确的。
然而，对于连续系统，如板壳，这一步就变得至关重要，因为它只能得到近似的结果。
因此解决方案的精确度取决于所使用的离散化方法。
离散化之后，应当选择最适当的单元类型来最接近地模拟实际的物理性能。
所用的单元总数和给定物体内单元大小和类型的变化是需要工程判断的主要问题。
单元必须小到可以给出有用的结果, 又必须足够大以节省计算费用。
选择每个单元内的位移函数。
该函数是用单元的节点值在单元内部定义的。
线性、二次和三次多项式是常常使用的位移函数，因为用它们建立有限元公式比较简单。
每个单元可重复使用同一个通用的位移函数。
因此有限元方法是这样一种方法:一个连续量，如整个物体内的位移，用一个离散的模型来近似，而此离散模型是由每个有限域或有限单元内定义的分片连续函数组成的。
为了推导被划分的有限个单元的方程，那么需要应变位移和应力应变关系，常用的方式有胡克定律。
单元刚度矩阵和单元方程的建立是根据刚度影响系数的概念,需要有结构分析的背景知识，常用的方法有之间平衡法、功和能量法以及加权残余法。
可以使用叠加法(称为直接刚度法，其理论基础为节点力平衡)将第 3.4 得出的单个单元方程加在一起得出整个结构的总体方程。
在直接刚度法中所隐含的是连续和协调概念，要求结构保持完整，在结构任何一处不发生撕开。
3.6方程在修改考虑了边界条件之后，形成一组联立代数方程组，可以写为扩展的矩阵阶数是未知自由度的结果总数；可用消元或迭代求解该方程。
3.73.844.1求出如图 1 所示，各个节点的位移，单元内力和支反力。
(弹簧弹性系数和节点受力均已在图中标识)图 1 弹簧单元实例图4.1.14.1.2>> k1=SpringElementStiffness(2500) >> k2=SpringElementStiffness(2500) >> k3=SpringElementStiffness(5000)4.1.3>> K=zeros(4,4)定义4×4 的空白矩阵。
>> K=SpringAssemble(K,k1,1,3)>> K=SpringAssemble(K,k2,2,3)>> K=SpringAssemble(K,k3,3,4)经过整合，得到总体刚度矩阵如下：2500 0 -2500 00 2500 -2500 0-2500 -2500 10000 -50000 0 -5000 5000 4.1.4对节点 1、2、4 的位移进行分析，U1 = 0,U2= 0,U4= 0 。
节点 3 所受合外力 F=10KN。
4.1.5>> k=K(3，3)>> f= [2]>> u=k\f>> U= [0;0;u(3) ;0]U=5.0E-03*5.0000如上矩阵所示，四个节点的位移均已用矩阵形式表示出来了 (向右为正) ，接着通过通过计算求出受力大小。
>> F=K*UF=-2.500-2.50010.000-5.000如上矩阵所示，四个节点的受力均已用矩阵形式表示出来了 (向右为正) ，接着通过通过计算求解各个单元内力。
4.1.6>> u1= [0;U(3)]>> u2= [0;U(3)]>> u3= [U(3) ;U(4)]>> f1=SpringElementForces(k1,u1)>> f2=SpringElementForces(k2,u2)>> f3=SpringElementForces(k3,u3)F1=-2.5002.500F2=-2.5002.500F3=5-5经由 MATLAB，解出各结果如下：以上即为利用MATLAB 与有限元法来对解弹簧单元的力学分析。
当然，有限元法不是只能用来分析弹簧单元，下面就让我们来看看，如何利用有限元法和 MATLAB 软件来求解弹簧元与平面桁架组成的结构。
如图 2 所示梁单元结构，假定弹性模量 E=200GPa 、惯性矩 I=70 × 10-5m4 、M=15kN m求：(1) 系统的总体刚度矩阵；(2) 各节点的位移、转角；(3) 各节点的支反力；(5) 每个单元的内力 (剪力、弯矩) ；(6) 每个单元的剪力图；(7) 每个单元的弯矩图。
图 2 梁单元实例图4.2.14.2.2>> E=200e6>> I=70e-5>> L1=3.5>> L2=2>> k1=BeamElementStiffness(E,I,L1)>> k2=BeamElementStiffness(E,I,L2)4.2.3建立6×6 矩阵：>> K=zeros(6,6)>> K=BeamAssemble(K,k 1, 1,2)>> K=BeamAssemble(K,k2,2,3)4.2.4位移、U1=0、U2=0、U3=0；弯矩：M1=0 、M2= 15 、M3=0。
4.2.5>> k=K([2,4,6], [2,4,6])>> f= [0;15;0]>> u=k\f>> U= [0;u(1) ;0;u(2) ;0;u(3)] >> F=K*U利用MATLAB 所得结果如下：U=1.0e-04*-0.4752-0.9504-0.4752F=3.1276-0.00006.255232.0000-7.42580.00004.2.6>> u1= [U(1) ;U(2) ;U(3) ;U(4)] >> u2= [U(3) ;U(4) ;U(5) ;U(6)] >> f1=BeamElementForces(k1,u1)>> f2=BeamElementForces(k2,u2)所得计算结果如下：U1=1.0e-04-0.47520.9504U2=1.0E-04*0.9504-0.4752>>BeamElementShearDiagram(f1, L)>>BeamElementShearDiagram(f2, L)>>BeamElementMomentDiagram(f1, L)>>BeamElementMomentDiagram(f2, L)用上述程序运行得出各单元的剪力图和弯矩图5在特殊条件下，通过十几周的学习，我已经对有限元的基础有了一个大致掌握，关于用有限元进行具体分析也掌握了一些最基本的方法。
