平面问题有限元解法公式推导讲解课件

平面问题的有限单元解法2/7/2023有限元单元法基本思想有限单元法的思想是将物体（连续的求解域）离散成有限个且按一定方式相互联结在一起的单元组合，来模拟或逼近原来的物体，从而将一个连续的无限自由度问题简化为离散的有限自由度问题求解的一种数值分析法。物体被离散后，通过对其中各个单元进行单元分析，最终得到对整个物体的分析。有限单元法的分析步骤如下：物体离散化单元特性分析单元组集，整体分析求解未知节点的位移由节点的位移求解各单元的位移和应力2/7/2023南京农业大学工学院机械工程系基本变量uεσ

（位移）（应变）（应力）基本方程力的平衡方程几何方程物理方程求解方法经典解析半解析传统数值解法现代数值解法（计算机硬件、规范化、标准化、规模化）物体变形及受力情况的描述三大方面三大方程即：σ=EεE弹性模量2/7/2023南京农业大学工学院机械工程系平面问题有限单元法基本概念有限单元法(FEM)是20世纪50年代以来随着计算机的广泛应用而发展起来的一种数值解法。简单地说，就是用结构力学方法求解弹性力学问题。平面问题的有限单元法求解将连续体变换成为离散化结构。即将连续体划分为有限多个有限大小的单元，这些单元仅在一些结点连接起来，构成一个所谓离散化结构。（对于平面问题，常用的单元是三角形单元）用结构力学方法进行求解2/7/2023南京农业大学工学院机械工程系有限元单元法分析步骤（一）结构离散化

将结构分成有限个小的单元体，单元与单元、单元与边界之间通过节点连接。结构的离散化是有限元法分析地第一步，关系到计算精度和效率，包括以下三个方面：单元类型的选择。选定单元类型，确定单元形状、单元节点数、节点自由度数等。单元划分。网格划分越细，节点越多，计算结果越精确，但计算量越大。网格加密到一定程度后计算精度提高就不明显，对应应力变化平缓区域不必要细分网格。节点编码。

注意：有限元分析的结构已不是原有的物体或结构物，而是由同样材料、众多单元以一定方式连接成的离散物体。所以，用有限元分析计算所获得的结果是近似的（满足工程要求即可）。2/7/2023南京农业大学工学院机械工程系有限元单元法分析步骤（二）单元特性分析

选择未知量模式选择节点位移作为基本未知量时，称为位移法；选节点力作为基本未知量时，称为力法；取一部分节点位移和一部分节点力作为未知量，称为混合法。分析单元力学性质根据单元材料性质、形状、尺寸、节点数目、位置等，找出单元节点力和节点位移关系式，应用几何方程和物理方程建立力和位移的方程式，从而导出单元刚度矩阵。计算等效节点力作用在单元边界上的表面力、体积力或集中力都需要等效地移到节点上去，即用等效力来替代所有作用在单元上的力。2/7/2023南京农业大学工学院机械工程系弹性力学中的几个基本概念作用于物体的外力可以分为体积力和表面力。体力：分布在物体体积内的力，如重力、惯性力。为了表明物体在某一点P所受体力的大小和方向，在这一点取物体的一小部分，它包含P点，而它的体积为△V，作用于其上的体力为△F，则体力的平均集度为△F/△V。当△V不断减小，假定体力为连续分布，则△F/△V将趋于一定的极限f，即:这个极限矢量f就是该物体在P点所受体力在集度。f的方向就是△F的方向，矢量f在坐标轴x,y,z上的投影fx,fy,fz称为该物体在P点的体力分量，以沿坐标轴正方向为正，沿坐标轴负方向为负。2/7/2023南京农业大学工学院机械工程系弹性力学中的几个基本概念面力：分布在物体表面上的力，如流体压力和接触力。为了表明物体在某一点P所受面力的大小和方向，在这一点取物体表面的一小部分，它包含P点，而它的面积为△S，作用于其上的面力为△F，则面力的平均集度为△F/△S。当△S不断减小，假定体力为连续分布，则△F/△S将趋于一定的极限

，即:这个极限矢量

就是该物体在P点所受面力在集度。

的方向就是△F的方向，矢量

在坐标轴x,y,z上的投影

称为该物体在P点的面力分量，以沿坐标轴正方向为正，沿坐标轴负方向为负。2/7/2023南京农业大学工学院机械工程系弹性力学中应力的方向规定每一个面上的应力可以分解为一个正应力和两个切应力。正应力用σ表示，加上一个下标字母，表示作用面和作用方向。切应力用τ表示，并加上两个下标字母，表示作用面和作用方向。前一个字母表示作用面垂直于哪一个坐标轴，后一个字母表示作用方向沿着哪一个坐标轴。2/7/2023南京农业大学工学院机械工程系平面问题的基本理论任何一个实际的弹性力学问题都是空间问题，但是如果所考察的弹性体具有某种特殊的形状，并且承受的是某些特殊的外力和约束，就可以把空间问题简化为近似的平面问题。两种典型的平面问题平面应力问题平面应变问题2/7/2023南京农业大学工学院机械工程系由于板很薄，外力不沿厚度变化，应力沿板的厚度又是连续分布的，所以可以认为在整个薄板的所有各点：只剩下平行于xy面的三个平面应力分量，即：这种问题成为平面应力问题。平面应力问题设有很薄的等厚度薄板，只在板边上受有平行于板面并不沿厚度变化的面力或约束。同时，体力也平行于板面不沿厚度变化。设薄板的厚度为δ。以薄板的中面为xy面，以垂直于中面的任何一直线为z轴。所以有：2/7/2023南京农业大学工学院机械工程系只剩下平行于xy面的三个形变分量，即：这种问题成为平面应变问题。由于z方向的位移处处为0，所以：，由于z方向的伸缩被阻止，一般平面应变问题设有很长的柱形体，它的横截面不沿长度变化，在柱面上受有平行于横截面而且不沿长度变化的面力或约束。同时，体力也平行于横截面不沿长度变化。假想该柱体为无限长，以任一横截面为xy面，以任一纵线为z轴，则所有一切应力分量、形变分量和位移分量都不沿z方向变化，而只是xy的函数，所有各点的位移矢量都平行于xy面，这种问题称为平面位移问题。由对称条件可知：由胡克定律，相应的切应变：2/7/2023南京农业大学工学院机械工程系平衡微分方程从弹性体中取出一个微分体，根据平衡条件导出应力分量与体力分量之间的关系式，也就是平面问题的平衡微分方程。从弹性体中取出一个微小的正平行六面体，它在x和y方向的尺寸分别为dx和dy，在z方向的尺寸为一个单位长度。以x为投影轴，列出投影的平衡方程：约简以后，两边除以dxdy，得：同理，以y为投影轴，列出投影的平衡方程，化简得：2/7/2023南京农业大学工学院机械工程系假定已知任一点P处坐标面上的应力分量σx，σy，τxy=τyx。求经过该点的，平行于z轴而倾斜于x轴和

y轴的任何倾斜面上应力。从在P点附近取一个平面AB，它平行于上述斜面，并经过P点划出一个微小的三棱柱PAB。当AB无限小而趋于P点时，平面AB上的应力就成为斜面上的应力。平面问题中一点的应力状态设斜面AB的长度为ds，则PB面及A面的长度分别为lds及mds，而PAB的面积为

ldsmds/2，棱柱的厚度设为1。由x轴平衡条件，得：其中，fx为体力分量。将上式除以ds，并令ds趋于0（斜面AB趋于P点），即得：由y轴平衡条件，得：用n表示斜面AB的外法线方向，其方向余弦为：2/7/2023南京农业大学工学院机械工程系几何方程经过弹性体内的任意一点P，沿x轴和y轴的正方向取两个微小长度的线段PA＝dx和PB＝dy。假定弹性体受力后，P,A,B三点分别移动到P’,A’,B’.线段PA的线应变是：注:由于位移微小，y方向的位移v引起的PA的伸缩，是高一阶微量，略去不计。线段PB的线应变是：线段PA与

PB之间的直角的改变，即切应变线段PA的转角α是：线段PB的转角β是：2/7/2023南京农业大学工学院机械工程系边界条件若在su部分边界上给定了约束位移分量和，则对于此边界上的每一点，位移函数u和v应满足条件：其中（u）s和（v）s是位移的边界值，和在边界上是坐标的已知函数。边界条件表示在边界上位移与约束，或应力与面力之间的关系式。它可以分为位移边界条件、应力边界条件和混合边界条件。位移边界条件：应力边界条件：若在su部分边界上给定了面力和，则由平衡条件得出平面应力问题的应力(或面力）边界条件为：其中，l,m是边界面外法线的方向余弦。2/7/2023南京农业大学工学院机械工程系圣维南原理在求解弹性力学问题时，应力分量、形变分量和位移分量必须满足区域内的三套基本方程，还必须满足边界上的边界条件。但是，要使边界条件得到完全满足，往往遇到很大的困难。圣维南原理可为简化局部边界上的应力边界条件提供很大方便。圣维南原理表明，如果把物体的一小部分边界上的面力，变换为分布不同但静力等效的面力（主矢相同，对同一点的主矩也相同），那么，近处的应力分布将有显著的改变，但是远处所受的影响可以不计。2/7/2023南京农业大学工学院机械工程系圣维南原理的应用例，设有柱形构件，在两端截面的形心受到大小相等而方向相反的拉力F（a）。如果把一端或两端的拉力变换为静力等效的力，则只有虚线划出的部分的应力分布有显著的改变，而其余部分所受影响是可以不计的。由于（d）图中，面力连续分布，边界条件简单，应力容易求得。其它三种情况，应力难以求得。把d情况下的应力解答应用到其它三个情况，虽不能满足两端的应力边界条件，但仍然可以表明离杆端较远处的应力状态，而没有显著的误差。图e，构件右端有位移边界条件，，d情况的解答，不能满足位移边界条件，但e图右端的面力，一定是合成为经过截面形心的力F。所以把图d情况的解答应用于图e时，仍然只是在靠近两端处有显著的误差，而在离两端较远之处，误差可以不计。2/7/2023南京农业大学工学院机械工程系圣维南原理的应用(续)应力的主矢量和主矩的绝对值分别等于面力的主矢量和主矩的绝对值；面力的主矢量和主矩的方向就是应力的主矢量和主矩的方向。2/7/2023南京农业大学工学院机械工程系有限单元法中基本量的矩阵表示有限单元法(FEM)中，为了简洁清晰地表示各个基本量以及它们之间的关系，也为了便于编制程序利用计算机进行计算，广泛采用矩阵表示和矩阵运算。平面问题中，物体受体力，可用体力列阵表示：(1)物体受面力，可用面力列阵表示：(2)3个应力分量的应力列阵表示：(3)3个形变分量的应变列阵表示：(4)2个位移分量的位移列阵表示：(5)2/7/2023南京农业大学工学院机械工程系虚位移原理(续)代入虚功方程，得：上式为集中力作用下的虚功方程。集中力列阵（13）虚位移列阵（14）外力在虚位移上所做的功为：2/7/2023南京农业大学工学院机械工程系（1）取三角形单元的结点位移为基本位置量：

(a)其中，

称为单元的结点位移列阵；（2）应用插值公式，由单元结点位移求出单元的位移函数： (b)其中，N称为形函数矩阵；（3）应用几何方程，由单元的结点位移求出单元的应变： (c)其中，B是表示与之间关系的矩阵；三角形单元离散化结构分析步骤2/7/2023南京农业大学工学院机械工程系

(f)其中，Fe

是单元的结点力，k称为单元劲度列阵；

对三角形板单元，节点力为：

(e)

（5）应用虚功方程，由单元的结点应力求出单元的结点力。假设把单元和节点切开，对右图中的i节点：节点对单元的作用力为节点力，作用于单元上。三角形单元离散化结构分析步骤(续）（4）应用物理方程，由单元的结点位移求出单元的应力：

(d)其中，S称为应力转换矩阵；

Fe是作用于单元的外力，此外，单元内部还作用有应力。根据虚功方程，可以将单元的节点力Fe用应力来表示，从而得到节点力的公式：2/7/2023南京农业大学工学院机械工程系（7）列出各结点的平衡方程，组成整个结构的平衡方程组。由于节点i受有环绕节点的单元移置而来的节点载荷和节点力因而i节点的平衡方程为： （i=1,2,…,n）(h)三角形单元离散化结构分析步骤(续）（6）应用虚功方程，将单元中的外力载荷向结点移置，化为结点载荷（即求出单元的节点载荷）： (g)将（f）代入（h），整理得：（j)其中，K称为整体刚度矩阵，FL是整体结点载荷列阵，δ是整体结点位移列阵。在上述求解步骤中，（2）至（6）是针对每个单元进行的，称为单元分析；（7）是针对整个结构进行的称为整体分析。2/7/2023南京农业大学工学院机械工程系对三角形三个结点i,j,m结点，位移函数应当等于该节点的位移值，即：

三角形单元的位移模式对每个单元，只要求得单元中的位移函数，就可以应用几何方程求得应变，再应用物理方程求得应力。有限单元法中常取结点位移为基本未知量，由单元的结点位移求出单元中的位移函数是首先必须解决的问题。可以假定一个位移模式，来表示单元中的位移函数（即在单元中做出位移插值函数）。三角形单元中，可以假定位移分量只是坐标的线性函数，即假定：6个方程解出α1-6，代入u，v式整理得：其中：2/7/2023南京农业大学工学院机械工程系三角形单元的位移模式Ni也可以该写成为：其中系数ai,bi,ci是：其中A就等于三角形ijm的面积：按照解析几何学，在图示的坐标系中，为了得出的面积A不致成为负值，节点i，j，m的次序必须是逆时针转向的。

Ni，Nj，Nm这三个函数，表明了单元ijm的位移次形态（也就是位移在单元内的变化规律），因而称为形态函数，简称形函数。2/7/2023南京农业大学工学院机械工程系三角形单元的位移模式位移模式的表示式可用矩阵表示为：简写为：其中是单元的节点位移列阵。是形态函数矩阵或形函数矩阵。有限单元法中，应力转换矩阵和劲度矩阵的建立以及载荷的移置等，都依赖于位移模式。2/7/2023南京农业大学工学院机械工程系简写为：其中矩阵B可写成分块形式：

其子矩阵为：单元的应变列阵和应力列阵利用几何方程和物理方程，求出单元中的应变和应力，用结点位移表示：将位移函数(16)和(18)代入几何方程(6)，得出用结点位移表示单元应变。2/7/2023南京农业大学工学院机械工程系将D表达式(9)和B表达式(27)代入上式，并写成分块形式，即得到平面应力问中的应力转换矩阵：单元的应力列阵（续）再将单元的应变式(26)代入物理方程(8)，得出用结点位移表示单元中应力的表达式。

其中子矩阵为：简写为：其中，2/7/2023南京农业大学工学院机械工程系由式(26)引起的虚应变为：由于结点力在虚位移上的虚功应当等于应力在虚应变上的虚功，即：

单元的结点列阵与劲度矩阵对于任一单元，均假设所受的外力载荷已经被移置到结点上，并且单元已经切开，如右图所示：单元只受到结点对单元的作用力，即结点力：假想在结点i,j,m处发生了虚位移，即：对单元而言，这些结点力是外力，使单元内部产生应力。2/7/2023南京农业大学工学院机械工程系从而建立了单元结点力和结点位移之间的关系。对于三角形单元，B中的元素为常量。并且，因此，k可简写为：

k称为单元的劲度矩阵。单元的结点列阵与劲度矩阵(续）由于中的元素是常量，并且虚位移的值可以是任意的：则将B和D表达式代入上式，得：令则式可以简写为2/7/2023南京农业大学工学院机械工程系载荷向节点移置，单元的载荷列阵设单元ijm在坐标为（x，y）的任意一点M，在单位厚度上受有集中载荷fP，其坐标方向的分量为fPx和fPy，用矩阵表示为fP=（fPx

fPy）T，将此集中力移置到单元的节点处，转换为节点载荷，并且单元节点载荷列阵表示为：假想单元的各点发生了虚位移：由位移模式，相应于集中力fP的作用点（x,y）的虚位移为：集中载荷的移置2/7/2023南京农业大学工学院机械工程系载荷向节点移置，单元的载荷列阵（续）由于虚位移可以是任意的，所以：把N的表达式（25）代入上式，上式改写为：其中，Ni,

Nj,

Nm，为它们在M点的函数值：根据静力等效原则，节点载荷在节点虚位移上的虚功，等于原载荷集中力在其作用点的虚位移上的虚功，即：2/7/2023南京农业大学工学院机械工程系载荷向节点移置，单元的载荷列阵（续）例，设单元ijm的密度为ρ，试求自重的等效节点载荷。分析：因为fx=0，

fy=-ρg，故由式（43）得：由设上述单元受有分布的体力f=（fx

fy）T，可将微分体积tdxdy上的体力ftdxdy当作集中力，利用（40）式积分，得到：体力的移置注意单元的自重为-ρgtA，可见移置到每个节点的载荷均为1/3自重。2/7/2023南京农业大学工学院机械工程系载荷向节点移置，单元的载荷列阵（续）由设上述单元的某一边上受有分布的面力，可将微分面积tds上的面力当作集中载荷，利用（40）式积分，得到：面力的移置例，设在ij边上受有沿x方向的均布面力q，试求等效节点载荷。分析：因为

，故由式（45）得：2/7/2023南京农业大学工学院机械工程系注意：式(46)和(48）中的编码i,j,m仅是每个单元的局部编码，对于整个结构，则将结点的平衡方程按整体结点编码1，2，…，n排列起来，就组成整个结构的结点平衡方程组：

整体的结构分析节点平衡方程组因此，结点i的平衡方程是：以上几节的分析都是针对单元进行的，即一将单元上的外力载荷都向节点移置而成为节点载荷；另一方面求出节点载荷与单元之间的相互作用力，如左图所示。节点对单元的作用力是节点力，相反，单元对节点的作用力是节点力的负值。于是，作用于结点i上的力，有结点载荷FLi

，和结点力的负值，即：其中，是对环绕结点i的单元求和，写成标量形式：

2/7/2023南京农业大学工学院机械工程系由整体平衡方程组，解出结点位移δ，便可由式(23)和(30)求出每个单元的位移函数和应力。整体的结构分析节点平衡方程组其中，整体结点位移列阵：整体结点载荷列阵：K是整体刚度矩阵，其元素是： 整个结构的结点平衡方程组即整体劲度矩阵的元素,Krs就是按整体节点编码的、同下标rs的单元劲度矩阵元素叠加而得到的。2/7/2023南京农业大学工学院机械工程系平面有限元解法（例）设有对角受压的正方形薄板（如上图所示），载荷沿厚度均匀分布，为2N/m。试对该结构进行整体分析，建立整体刚度矩阵和整体结点载荷列阵，建立整体结点方程组，通过编程求解出结点的位移，并从而求出各单元的应力。（为简单起见，取板的厚度t=1,弹性常数E=1,泊松比μ＝0）2/7/2023南京农业大学工学院机械工程系平面有限元解法——划分单元由于平面薄板沿xz面和yz面均对称，所以只取1/4之一部分作为分析和计算对象。将对象划分成4个单元，共有6个结点，单元和结点上均编上号码，其中结点的整体编码1至6，以及个单元的结点局部编码i,j,m,均示于上图中。单元号ⅠⅡⅢⅣ局部编码整体编码i3526j1253m24352/7/2023南京农业大学工学院机械工程系平面有限元解法——整体劲度矩阵每个单元，结点的局部编码和整体编码对应关系已经确定，每个单元劲度矩阵中任一子矩阵在整体劲度矩阵中的位置及其力学意义也就明确了。如单元Ⅰ的kii，即k33，它的四个元素表示当结构的结点3沿x或y方向有单位位移时，在结点3的x方向或y方向引起的结点力。暂时不考虑位移边界条件，把所分析结构的整体结点平衡方程组列出：整体劲度矩阵写成6×6的矩阵，它的每个子块是2×2的矩阵，实际它是一个12×12的矩阵。如K23，它的四个元素表示当结构的结点3沿x或y方向有单位位移时，在结点2的x方向或y方向引起的结点力。2/7/2023南京农业大学工学院机械工程系平面有限元解法——整体劲度矩阵续由于于结点3和结点2在结构中是通过Ⅰ和Ⅲ这两个单元相联系，因而K23应是单元Ⅰ的k23和单元Ⅲ的k23之和。同理，可以找到各单元劲度矩阵中所有子矩阵在整体劲度矩阵K中的位置，得到整体劲度矩阵。式中k的上标1,2,3,4表示是哪一个单元的劲度矩阵中的子矩阵，空白处是2×2的零矩阵。2/7/2023南京农业大学工学院机械工程系平面有限元解法——整体劲度矩阵续对于单元Ⅰ、Ⅱ、Ⅳ，根据公式，可求得A=0.5m2，将上式中各子块的具体数值代入整体刚度矩阵K表达式中，得出整体刚度矩阵。对于单元Ⅲ，根据公式，可求得A=0.5m2，把μ＝0，t＝1m，代入单元的劲度矩阵，得两种单元的劲度矩阵k都是：(37’)2/7/2023南京农业大学工学院机械工程系平面有限元解法——整体劲度矩阵续整体刚度矩阵K(38’)2/7/2023南京农业大学工学院机械工程系平面有限元解法——