弹性力学与有限元

1、郭吉坦郭吉坦Email：dl-QQ： 1628295008 手机：第一篇 弹性力学第一章 弹性力学基本方程1.1 绪论1.2 弹性力学的基本假定1.3 几个基本概念1.4 弹性力学基本方程第二章 弹性力学平面问题2.1 平面应力问题2.2 平面应变问题2.3 平面问题的基本方程第三章 弹性力学问题求解方法简述 第一章 弹性力学基本方程1.1 绪论1.2 弹性力学的基本假定1.3 几个基本概念1.4 弹性力学基本方程应力 应变 位移弹性体外界作用 弹性力学基本内容外力温度变化 弹性力学弹性力学，又称弹性理论。 是研究弹性体由于外力载荷或者温度改变，物体内部所产生的位移、

2、变形和应力分布等。为解决工程结构的强度，刚度和稳定性问题作准备 。 弹性力学的研究对象：是完全弹性体，包括构件、板和三维弹性体，比材料力学和结构力学的研究范围更为广泛 。 研究的内容： 外力作用下 应力、应变、位移1.1 弹性力学绪论 物体变形物体变形弹性变形、塑性变形弹性变形、塑性变形 弹性变形：弹性变形： 当外力撤去以后恢复到原始状态，没有变形残留，材当外力撤去以后恢复到原始状态，没有变形残留，材料的应力和应变之间具有一一对应的关系。与时间无料的应力和应变之间具有一一对应的关系。与时间无关，也与变形历史无关关，也与变形历史无关。 塑性变形：塑性变形： 当外力撤去以后尚残留部分变形量，不能恢

3、复到原始当外力撤去以后尚残留部分变形量，不能恢复到原始状态，状态，即存在永久变形。应力和应变之间的关系即存在永久变形。应力和应变之间的关系不再一一对应，与时间、与加载历程有关不再一一对应，与时间、与加载历程有关。 弹性：弹性：假定“完全弹性”关系，是抽象出来的理想模型。 完全弹性是指在一定温度条件下，材料的应力和应变之间具有一一对应的关系。 应力应变关系称为本构关系。本构关系。 材料模型包括：线性弹性体非线性弹性体1.2 弹性力学的基本假定1. 连续性假设连续性假设 根据这一假设，物体的所有物理量，例如位移、应变和应力等均成为物体所占空间的连续函数。2. 均匀性假设均匀性假设 假设弹性物体是由

4、同一类型的均匀材料组成的，物体各个部分的物理性质都是相同的，不随坐标位置的变化而改变。在处理问题时，可以取出物体的任意一个小部分讨论。 3. 各向同性假设各向同性假设 假定物体在各个不同的方向上具有相同的物理性质，物体的弹性常数不随坐标方向变化。 像木材、竹子以及纤维增强材料等，属于各向异性材料，它们是复合材料力学研究的对象。 4. 完全弹性假设完全弹性假设 应力和应变之间存在一一对应关系，与时间及变形历史无关。满足胡克定理。5. 小变形假设小变形假设 在弹性体的平衡等问题讨论时，不考虑因变形所引起的几何尺寸变化，使用物体变形前的几何尺寸来替代变形后的尺寸。采用这一假设，在基本方程中，略去位移

5、、应变和应力分量的高阶小量，使基本方程成为线性的偏微分方程组。1.3 几个基本概念1. 外力2. 一点的应力状态3. 一点的形变4. 位移分量 作用于物体的外力可以分为3种类型： 体力、面力、集中力。 体力体力就是分布在物体整个体积内部各个质点上的力，又称为质量力。例如物体的重力，惯性力，电磁力等等。 面力面力是分布在物体表面上的力，例如风力，静水压力，物体之间的接触力等。 集中力集中力作用物体一点上的力。（在弹性力学中一般不用，而在有限元中经常出现）1 外力 体力 物体任意一点P 所受体力的大小和方向，在P点区域取一微小体积元素V， 设V 的体力合力为F，则V 的平均体力为当V 趋近于0，则

6、为P点的体力 体力是矢量：一般情况下，物体每个点体力的大小和方向不同。 体力分量：将体力沿三个坐标轴xyz 分解，用X、Y、Z表示，称为体力分量。 符号规定：与坐标轴方向一致为正，反之为负。 应该注意的是：在弹性力学中，体力是指单位体积的力 。 体力的因次：力/长度3 表示：F=X Y Z 面力 与体力相似，在物体表面上任意一点P 所受面力的大小和方向，在P点区域取微小面积元素S ，当S 趋近于0，则为P点的面力面力分量符号规定：与坐标轴方向一致为正，反之为负。面力的因次：力/长度2 集中力 体力与面力都是分布力，集中力则只是作用在一个点上，作用区域V或S很小，但数值很大，这种形式的力可以认为

7、是集中力。 集中力分量：集中力直接将其沿三个坐标轴分解，用X0、Y0、Z0表示，即集中力力分量。 符号规定：与坐标轴方向一致为正，反之为负。 体力的因次：力2 一点的应力状态 应力表示方法 材料力学中接触过斜截面上的应力，斜截面上应力可以分成正应力、剪应力； 复杂物体任意截面上的应力可分为 1个与平面垂直的正应力、 2个平面内剪应力。X面Y面Z面正应力分量 3个：xyz、xyxzyxyzzxzy、剪应力分量 6个：正面 负面X面Y面Z面应力符号意义xyz、xy剪应力：正应力： 由法线方向确定作用面 作用方向 符号规定： 正面上与坐标轴正向一致，为正； 负面上与坐标轴负向一致，为正。剪应力互等定

8、理xyzxyyzzx、剪应力不再区分哪个是作用面或作用方向 。应力分量：xyyxyzzyxzzx相等 xyzxyyzzx3 一点应变分量 微分单元体的变形： 微分单元体棱边的伸长和缩短；正应变正应变 棱边之间夹角的变化；剪应变剪应变 xyzxyyzzx正应变分量 3个：剪应变分量 3个：xyz、xyyzzx、 应变的定义（自学）设平行六面体单元，3个轴棱边： 变形前为MA，MB，MC； 变形后变为MA，MB，MC。xyz、正应变（小变形） （自学）符号规定： 正应变以伸长为正。剪应变（自学）符号规定： 正应变以伸长为正；剪应变以角度变小为正。4 位移分量 位移：由于载荷作用或者温度变化等外界因

9、素等影响，物体内各点在空间的位置将发生变化，位置移动即产生位移。 位移刚体位移刚体位移、变形变形 刚体位移刚体位移物体内部各个点仍然保持初始状态的相对位置不变，由于物体整体在空间做刚体运动引起的位置改变。 变形变形物体整体位置不变，弹性体在外力作用下发生形状的变化，而改变了物体内部各个点的相对位置，引起位移。 后者与弹性体的应力有着直接的关系弹性力学研究的主要变形，通常叫位移。u=x（x，y，z）-x=u（x，y，z） v=y（x，y，z）-y=v（x，y，z） w=z（x，y，z）-z=w（x，y，z）根据连续性假设，弹性体在变形前和变形后仍保持为连续体。 弹性体中某点在变形过程中由M（x，

10、y，z）移动至M(x，y，z)，这一过程也是连续的，为 x、y、z的单值连续函数 ufvw 形变和位移之间的关系： 位移确定位移确定 形变完全确定：形变完全确定： 从物理概念看，各点的位置确定，则微分线段上的形变确定 。 从数学推导看，位移函数确定，则其导数（形变）确定 。 形变确定，位移不完全确定形变确定，位移不完全确定 ： 从物理概念看，、确定，物体还可作刚体位移。 从数学推导看，、确定，求位移是积分运算，出现待定函数。应力 应变 位移弹性力学各个量之间的关系平衡方程物理方程几何方程外力 弹性力学分析过程中： 通过静力平衡、几何变形和本构关系建立起外力、应力、应变、位移之间相互关联。 再必

11、须根据已知物理量，（一般外力、结构几何形状和约束条件等），推导和确定基本未知量（应力、应变、位移。1.4 弹性力学基本方程1. 平衡方程（应力外力之间的关系）2. 物理方程（应变应力之间的关系）3. 几何方程（柯西方程 ）（应变位移之间的关系）4、变形协调方程 5、边界条件如果物体表面的面力已知，则称为应力边界条件： 第一类边界条件 如果物体表面的位移已知，则称为位移边界条件： 第二类边界条件混合边界条件 = 第一类+第二类5、边界条件应力边界条件：位移边界条件：cos,cos,cos,NxlNymNzn外法线的方向余弦方程数量：平衡方程3个物理方程6个几何方程6个合计 15xyzxyyzzx

12、、未知量： 应力分量6个应变分量6个位移分量3个 u、v、w 合计 15xyzxyyzzx、 、空间问题第二章 弹性力学平面问题2.1 平面应力问题2.2 平面应变问题2.3 平面问题的基本方程2.1 平面应力问题1、平面应力问题的概念 平面应力问题讨论的弹性体为薄板。薄壁厚度远小于结构另外两个方向的尺度。薄板的中面为平面，其所受外力，包括体力均平行于中面O-xy面内，并沿厚度方向z不变。而且薄板的两个表面不受外力作用。 平面应力问题 几何特征 薄壁厚度为h远小于结构另外两个方向的尺寸 等厚度 中心层平直 受力特征 外力平行于中心层 外力沿厚度不变化 根据薄板的表面面力边界条件，即表面不受外力

13、作用，则 由于板很薄，外力沿厚度均匀分布，因此应力分量也沿厚度均匀分布，应力分量不随z改变。 2、 平面应力问题的应力 应力分量 应变分量xyxy、 xyxyz 0 xyxy、 xyxy0yzzx=3、平面应力问题应力、应变1 平面应变问题的概念 弹性体是具有很长的纵向轴的柱形物体，横截面大小和形状沿轴线长度不变；作用外力与纵向轴垂直，并且沿长度不变；柱体的两端受固定约束。 可以认为柱体是无限长的。如果从中任取一个横截面，则柱形物体的形状和所受载荷将对此横截面是对称的。因此物体变形时，横截面上的各点只能在其自身平面内移动。2.2 平面应变问题 几何特征 一个尺寸远大于结构另外两个方向的尺寸 中

14、心轴平直 沿中心轴截面不变化 受力特征 外力垂直于中心轴 外力沿中心轴长度方向不变化 平面应变问题2、平面应变问题的位移 沿纵向轴的位移恒等于零； 由于无限长，所以任一个横截面都是一样的，与z轴无关。 只要是x、y坐标函数 应力分量xyxy、 xyxyxyxy、 xyxy0zyzzx=应变分量0zyzzx0 3、平面应变问题的应力、应变2.3 平面问题的基本方程1.平衡方程（应力外力之间的关系）2. 几何方程（应变位移之间的关系）3. 物理方程（应变应力之间的关系）平面应力与平面应变问题的： 平衡方程、几何方程相同。 但物理方程不同。 从空间问题推得。 平面应力的物理关系 平面应力的物理关系2

15、10101 1 002ED D 平面应变的物理关系0zyzzx= 平面应变的物理关系 D1010(1)(1-2 )1 2002ED平面应变问题平面应力问题z向应力分量 z =n ( x + y ) z0z向位移分量w0w0正应变分量二者主要不同在于z向应变，位移和正应力的计算公式 两种平面问题的区别 两种平面问题的内在关系平面应力平面应变.1 ,12EE.1 ,)1 ()21 (2EE平面应力平面应变平面应变平面应力 两种平面问题的内在关系1010(1)(1-2 )1 2002ED平面应力平面应变.1 ,12EE平面应力平面应变210101 1 002ED D 4 变形协调方程平面应力平面应变

16、由6个简化为1个 调和方程方程数量：平衡方程2个物理方程3个几何方程3个合计 8xyxy、未知量： 应力分量3个应变分量3个位移分量2个 u、v 合计 8xyzxy、平面问题第三章弹性力学问题求解方法简述应力 应变 位移弹性力学各个量之间的关系平衡方程物理方程几何方程外力3.1 概述 根据几何方程和本构方程可见： 位移、应力和应变分量之间不是相互独立的。 假如已知位移分量，通过几何方程可以得到应变分量，然后通过物理方程可以得到应力分量。 如果已知应力分量，通过物理方程得到应变分量，再由几何方程的积分求出位移分量，不过这时的应变分量必须满足一组补充方程，即变形协调方程。应力 应变 位移 位移解法

17、：若以位移函数作为基本未知量求解，根据物理方程和几何方程，应力分量及平衡方程均由位移分量表达； 应力解法： 若以应力函数作为基本未知量，称为应力解法，对于应力解法，应力分量必须满足平衡微分方程和变形协调方程 ； 混合解法 ： 若以位移分量和应力分量作为基本未知量，通过物理方程中消去应变分量，表述基本方程，称为混合解法。基本方程的求解方法弹性力学是对整个研究对象建立平衡方程、几何方程、物理方程，再根据外力作用下求整体的应力、应变、位移。 解答的途径有两大类：1. 精确解（解析解、理论解法） 逆法、半逆法、复变函数法、 级数法、 特殊函数法等2. 近似解法（数值解法） 1 位移解法 当位移分量作为

18、基本未知函数求解时，变形协调方程是自然满足的。根据物理方程和几何方程，可以得到：以位移表示的平衡微分方程，称为拉梅（拉梅（Lam）方程。）方程。拉普拉斯运算符号， 3.2 解析解法 2 应力法 主要介绍应力函数法，应力函数法， 称为艾里（Airy）应力函数。 设( , )x y应力表示的变形协调方程 双调和方程双调和方程 应力函数（1） 一次多项式一次多项式应力函数对应无应力的应力状态。 这个结论说明在应力函数中增加或减少一个x，y 的线性函数，将不影响应力分量的值。( , )x yaxbyc（2) 二次多项式如仅a，b，c0，分别表示单向拉伸或者纯剪切应力状态。22( , )x yaxbxy

19、cy(3) 三次多项式如果仅考虑d不为零的情况，即a=b=c=0，其对应于矩形梁的纯弯曲应力状态。 3223( , )x yaxbx ycxydy 解析解的难点： 弹性力学研究对象是弹性体，形体复杂，是偏微分方程的边值问题。在数学上求解困难重重，除了少数特殊边界问题，一般弹性体问题很难得到解答。 要得到解析解： 1、简化形体，譬如材料力学的研究对象是杆件，常微分方程，可以求解；平面问题，忽略次要因素，简化应力状态。 2、简化边界约束条件，放松某些限制等。 结果： 寻求求解偏微分方程在特定条件下的数学解法，而造成所得到的结果并非实际问题的真实状态。结果误差很大，甚至是错误的结论。 近似解法（数值

20、解法） 差分法 加权余量法 变分法 有限元法（FEM） 边界元法（BEM）3.3 数值解法 有限元法与边界元法的比较 离散化，FEM在区域上，BEM在边界上； 维数， BEM降维，3D 2D；2D 1D； 通用性， FEM格式统一， BEM特定问题； 对使用者数学要求， FEM低， BEM高； 目前应用状况，FEM一统天下。1 有限元基本思想2 离散化（建立计算模型）3 位移插值函数4 单元分析5 等效结点载荷6 整体分析7 有限元方程求解方法8 应力结果 9 举例第二篇 有限元法基础应力 应变 位移弹性力学各个量之间的关系平衡方程物理方程几何方程外力1 有限元基本思想应力 应变 位移放弃物理

21、方程几何方程外力有限元的基本思路能量原理只要位移场确定，就可得到应变、应力。 有限元的基本思想： 在弹性体内选取足够多、有限个点，假定这些点的位移已知，再用这些假定的位移量描述其它位置点的位移，就得到了用特定点位移表示的弹性体的位移场。 这些选定的有代表性的点结点，（node） 结点：代表性尖点、拐角、截面改变处等 集中载荷作用、位移约束位置等。 位移场：某个点（非结点）位移不是由所有结点位移来表述的，而是划分成小区域/小块上的结点来表示的，这些小区域/小块单元。 有限元处理问题的方法连续体剖分小块（单元），即离散体。 有限元法特点：1. 概念浅显，容易掌握，可以在不同程度上理解与应用2. 通

22、用性强，应用广泛，几乎所有领域；3. 计算格式统一，便于编程计算；4. 大型通用程序成熟商业化，无需专门知识编程5. 先进的前处理，网格自动划分， 完善的后处理，可视或动态显示，直观形象。误差难估计2 离散化（计算模型）单元的形式是多样的实体单元模型单元类型单元类型维数维数主要应用主要应用杆单元杆单元2-D承受轴向力作用，承受轴向力作用，平面桁架结构3-D空间桁架、网架等梁单元梁单元2-D主要承受横向载荷作用，即承受弯矩；主要承受横向载荷作用，即承受弯矩；也可也可横向横向+轴向轴向3-D固体固体2-D平面应力、平面应变问题：平面应力、平面应变问题：三角形、四边形3-D一般三维（空间）问题：一般

23、三维（空间）问题：四面体、八面体板壳板壳2-D主要承受横向载荷作用主要承受横向载荷作用 三角形、四边形3-D2.1单元类型与作用单元类型与作用杆单元杆单元梁单元梁单元二维单元二维单元线性单元线性单元二次单元二次单元三维单元三维单元线性单元线性单元二次单元二次单元板壳单元板壳单元2.2 离散化应注意的问题：首要的问题是根据结构的几何特点、受力特征选择合理的单元形式。对称性的利用，在划分单元之前，有必要先研究一下计算对象的对称或反对称的情况，以便确定是取整个物体，还是部分物体作为计算模型。取四分之一作为计算模型（以平面三角形单元为例）1.共边：覆盖求解区域，单元间既不允许相互重叠，也不允许相互脱开

24、；2.共点：任意三角形的顶点必须是相邻单元的顶点，不能为相邻单元的内点。3.边长接近：单元的边长尽可能接近，采用锐角三角形4.数目与精度兼顾：单元划分细，计算精度越高，但结点数增加，计算时间加长。单元大小过渡，应力梯度大的区域单元尺寸小，应力变化小的区域，单元可以划分大些。或在初步计算的基础上对于高应力区，在进一步细化网格，进行二次分析。5.适当简化。不可以不可以可以可以较差较差较好较好 节点编号顺序节点编号顺序 在进行节点编号时，应该注意要尽量使同一单元的相邻节点的号码差尽可能地小，以便最大限度地缩小刚度矩阵的带宽，节省存储、提高计算效率。 平面问题的半带宽为 B =2 (d+1)1 2 3

25、 4 5 6 71 3 5 7 9 11 138 9 10 11 12 13 142 4 6 8 10 12 143 位移模式 iiijejjmmmuvuvuv要求：要求：i、j、m按逆时针排序按逆时针排序单元的结点位移向量单元的结点位移向量用来描述单元内各点位移变化规律的函数，称为位移模式 三角形单元的位移模式假定为 yxvyxu654321123123123:iiijjjmmmPiuxyPjuxyPmuxymmjjiimmjjiivNvNvNvuNuNuNu位移模式：000000iiijmjijmjmmuvNNNuuNNNvvuv 位移模式矩阵表达：位移模式通式 efN f单元内任一点的位

26、移；单元内任一点的位移； e单元的结点位移向量；单元的结点位移向量； N单元的形函数矩阵。单元的形函数矩阵。 形函数的性质 1()21()21()2iiiijjjjmmmmNab xc yANab xc yANab xc yAijjAPjmAPmiAPij jimijmAAALLLAAA面积坐标iijjmmNLNLNL形函数与面积坐标的关系1211iijjmmxyAxyxy三角形的面积 efN位移模式 反映了单元内任意一点的位移与结点位移之间的关系。是有限元计算精度的关键。4 单元分析 4.1单元上任意一点的应变 4.2单元上任意一点的应力 4.3单元的能量 4.4 单元刚度矩阵的性质 4.1

27、单元上任意一点的应变 xyxyuxvyuyvx mmjjiimmjjiivNvNvNvuNuNuNu几何方程000000ijimixjjimyjxymjjiimmmuNNNvxxxNuNNvyyyuNNNNNNvyxyxyx 12000000bbbccccbcbcbijmijmiijjmme1()2iiiiNab xc yA Be或写成通式B 矩阵叫做单元几何矩阵，矩阵叫做单元几何矩阵，反映了单元内任意一点的应变分量与结点位移之间的关系 00010002ijmijmiijjmmbbbBcccAcbcbcb 几何矩阵几何矩阵B 中的每个元素，中的每个元素，均为常数，它们由结点坐标确定。单元内任意

28、一点P的应变分量与坐标（x，y）无关，说明单元中应变是常量。 4.2单元上任意一点的应力物理方程物理方程 D D D弹性矩阵弹性矩阵210101 1 002ED1010(1)(1-2 )12002ED平面应力问题平面应力问题平面应变问题平面应变问题 D D Be SD B Se S S 叫做应力矩阵叫做应力矩阵4.3单元的能量 1、单元的应变能121()2xxyyxyxy 一维问题应变能密度为平面问题应变能密度为 12TeUdxdy 1122TTTTeeeeeUSBdxdyBDB dxdy TeKBDB dxdy 12TeeeeUK eK 称为单元刚度矩阵，简称单刚，它反映了单元应变能与单元结

29、点向量之间的关系。 2、外力势能 (1)、体力势能 1()TeVXuYv dxdyfp dxdy (2)、面力势能 2()TeSSVXuYv dsfp ds (3)、集中力势能 3TeoeeVP 3 、单元的总势能 123()eeeeeeeUVUVVV 01()2()TeeeeTTTeeeeeeKPPP 101101020200103121 1213014002020101101eEk4.4 单元刚度矩阵的性质 某单元的刚度矩阵，仔细看看，会发现该矩阵有哪些特点？4.4 单元刚度矩阵的性质1、对称性 单元刚度矩阵是对称方阵，其元素都对称于主对角线。2、奇异性 单元刚度矩阵中任意一行或列元素之和

30、为零。其物理意义 是在没有给单元施加任何约束时，单元可有刚体运动，位移不能唯一的确定。3、主对角线元素恒为正值 主对角线元素是正值说明结点位移方向与施加结点荷载的方向是一致的。 4、单元刚度矩阵与单元位置无关 单元刚度矩阵与单元位置无关，也就是单元在平移时，Ke不变；单元结点排列顺序不同时，Ke中元素大小不变，而排列顺序相应改变。 弹性体所受外力包括体积力、表面力、集中力弹性体所受外力包括体积力、表面力、集中力。分别作用在弹性体内部、物体表面上、物体的一。分别作用在弹性体内部、物体表面上、物体的一个点上。个点上。 载荷列阵载荷列阵RR，是由弹性体的全部单元的等效，是由弹性体的全部单元的等效节点

31、力集合而成，是将全部载荷转移到单元的节点节点力集合而成，是将全部载荷转移到单元的节点上，它们的作用位置发生了变化上，它们的作用位置发生了变化载荷移置。载荷移置。 它们的作用效果是等效的，故称等效节点力向它们的作用效果是等效的，故称等效节点力向量量RRe e 。 各种载荷各种载荷 分别移置到节点上，再逐点加以合分别移置到节点上，再逐点加以合成求得单元的等效结点载荷。成求得单元的等效结点载荷。 1、体力等效结点载荷 PNp tdxdyeT自重情况下：自重情况下： TP000333eWWWy0 xijmlq2ql2ql TeqltQ0010102y0 xijml0q20lq3l TeltqQ0320

32、3100202、面力等效结点载荷 QNq tdseT6.1 结构的结点位移向量1 2 3 4 5 6 71 3 5 7 9 11 132 4 6 8 10 12 14 假设弹性体被划分为假设弹性体被划分为N N个单元和个单元和n n个节点，整个弹性个节点，整个弹性体的节点位移向量体的节点位移向量 2n2n1 1 2112nTTnTT RRRRnTTnTT2112整个弹性体的载荷列阵整个弹性体的载荷列阵 R R 2 2n n1 16 整体分析 矢量有方向性，外力、应力，不能直接相加 标量没有方向，只有大小，可以相加。 弹性体的能量是标量，可以直接相加。6.2 结构的总势能 eeeUV 01()(

33、)2TTeeeeeeeeeKPPP 单刚的扩充为了实现上述运算 6*6eK 2n*2neK（） nmjinmjikkkkkkkkkkmmmjmijmjjjiimijiinn1 122扩展 eeKK eePP eePP 00eePP结构的总刚度矩阵 结构总的体力列阵 结构总的面力列阵 结构总的集中力列阵 0FPPP 12TTKF 结构的总势能 12TTKF 6.3 6.3 整体刚度矩阵形成方法整体刚度矩阵形成方法123421q图图5组装总刚组装总刚k的一般规则：的一般规则：1.当当krs中中r=s时，该点被哪几个单元所共有，则总刚子矩时，该点被哪几个单元所共有，则总刚子矩阵阵krs就是这几个单元

34、的刚度矩阵子矩阵就是这几个单元的刚度矩阵子矩阵krse的相加。的相加。2.当当krs中中r s时，若时，若rs边是组合体的内边，则总体刚度矩阵边是组合体的内边，则总体刚度矩阵krs就是共用该边的两相邻单元单刚子矩阵就是共用该边的两相邻单元单刚子矩阵krse的相加。的相加。3.当当krs中中r和和s不同属于任何单元时，则总体刚度矩阵不同属于任何单元时，则总体刚度矩阵krs=0。下面，我们考查一个组装总刚的实例：下面，我们考查一个组装总刚的实例：1.整体刚度矩阵及载荷列阵的组集整体刚度矩阵及载荷列阵的组集根据叠加原理，整体结构的各根据叠加原理，整体结构的各个刚度矩阵的元素显然是由有关单个刚度矩阵的

35、元素显然是由有关单元的单元刚度矩阵的元素组集而成元的单元刚度矩阵的元素组集而成的，为了便于理解，现结合图的，为了便于理解，现结合图5说说明组集过程。明组集过程。 子块33333231232221131211KKKKKKKKKKe图中有两种编码：一是节点总码：图中有两种编码：一是节点总码：1、2、3、4；二是节；二是节点局部码，是每个单元的三个节点按逆时针方向的顺序各自点局部码，是每个单元的三个节点按逆时针方向的顺序各自编码为编码为1，2，3。图中两个单元的局部码与总码的对应关系为：图中两个单元的局部码与总码的对应关系为：单元单元1：1，2，31，2，3单元单元2：1，2，33，4，1或：或：单

36、元单元1：1，2，31，2，3单元单元2：1，2，31，3，4单元单元e的刚度矩阵分块形式为：的刚度矩阵分块形式为： 子块4444434241343332312423222114131211KKKKKKKKKKKKKKKKK整体刚度矩阵分块形式为：整体刚度矩阵分块形式为：其中每个子块是按照节点总码排列的。其中每个子块是按照节点总码排列的。通常，采用刚度集成法或直接刚度法来组集整体结构刚通常，采用刚度集成法或直接刚度法来组集整体结构刚度矩阵。刚度集成法分两步进行。度矩阵。刚度集成法分两步进行。 eK eK第一步，把单元刚度矩阵第一步，把单元刚度矩阵扩大成单元的贡献矩阵扩大成单元的贡献矩阵，使单元

37、刚度矩阵的四个子块按总体编号排列，空白处作零子使单元刚度矩阵的四个子块按总体编号排列，空白处作零子块填充。块填充。 2K第二步，以单元第二步，以单元2为例，局部码为例，局部码1，2，3对应于总码对应于总码3，4，1，按照这个对应关系扩充后，可得出单元，按照这个对应关系扩充后，可得出单元2的贡献矩阵的贡献矩阵。 213444342413433323124232221141312112KKKKKKKKKKKKKKKKK总码总码1234234312局部码局部码 1K用同样的方法可得单元用同样的方法可得单元1的贡献矩阵的贡献矩阵。 K第三步，把各单元的贡献矩阵对应行和列的子块相叠加第三步，把各单元的贡

38、献矩阵对应行和列的子块相叠加，即可得出整体结构的刚度矩阵，即可得出整体结构的刚度矩阵，如（，如（42）式。）式。 K22nn22 在这里应该指出，整体刚度矩阵在这里应该指出，整体刚度矩阵中每个子块为中每个子块为阶阶矩阵，所以若整体结构分为矩阵，所以若整体结构分为n个节点，则整体刚度矩阵的阶个节点，则整体刚度矩阵的阶数是数是。 TeynxnyxyxeFFFFFFFF 212211 总码总码1234123（42）123局部码局部码 32132100432123323223122322213313222113112312212121321211311221111121KKKKKKKKKKKKKKKK

39、KKKKee F eF至于整体结构的节点载荷列阵至于整体结构的节点载荷列阵的组集，只需将各单元的组集，只需将各单元的等效节点力列阵的等效节点力列阵扩大成扩大成2n行的列阵，然后按各单元的节行的列阵，然后按各单元的节点位移分量的编号，对应相叠加即可点位移分量的编号，对应相叠加即可6.4整体刚度矩阵的性质整体刚度矩阵的性质刚度矩阵刚度矩阵K中每一列元素的物理意义为：欲使弹中每一列元素的物理意义为：欲使弹性体的某一节点在坐标轴方向发生单位位移，而其它节性体的某一节点在坐标轴方向发生单位位移，而其它节点都保持为零的变形状态，在各节点上所需要施加的节点都保持为零的变形状态，在各节点上所需要施加的节点力。

40、点力。正定性，刚度矩阵正定性，刚度矩阵K中主对角元素总是正的。中主对角元素总是正的。刚度矩阵刚度矩阵K是一个对称矩阵，即是一个对称矩阵，即Krs =Ksr T。刚度矩阵刚度矩阵K是一个稀疏矩阵。如果遵守一定的节点是一个稀疏矩阵。如果遵守一定的节点编号规则，就可使矩阵的非零元素都集中在主对角线编号规则，就可使矩阵的非零元素都集中在主对角线附附近呈带状。近呈带状。5.奇异性。刚度矩阵奇异性。刚度矩阵K是一个奇异矩阵，在排除刚体是一个奇异矩阵，在排除刚体位移后，它是正定阵。位移后，它是正定阵。 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 1 2 3 4 5 6 7 8 9

41、 10 11 12 13 14 15 半带存储半带存储 半带宽B=（相邻节点号的最大差值D+1）*27 有限元方程及求解方法 7.1 有限元方程 结构的总势能 12TTKF 最小势能原理，对于线弹性体，某一变形可能位移状态为真实位移状态的必要和充分条件是，此位移状态的变形体势能取最小值 。 0 结构总势能泛函对结点位移的变分为0. 结构有限元方程 KF 它是一个2n阶的线性代数方程组。因为该方程中K是结构的总刚度矩阵，F是外荷载列阵，都通过计算求得，因此可以根据有限元方程可以确定结点位移。 7.2 位移边界条件的处理 由于总体刚度矩阵是奇异的，物理意义是结构中存在刚体位移，不能直接求解。必须引

42、入限制结构刚体位移的位移边界条件，即位移约束条件，消除总体刚度矩阵的奇异性，才能求解结构有限方程。 位移边界条件是指结构的某些区域位移已知，对于离散体来说，位移约束条件是某些结点的位移分量受到限制，包括位置限制和方向限制两个方面。具体哪些结点受到限制，受限制结点哪个方向位移分量受到限制，要根据结构受力后变形特征来确定。 处理的方法，主要有三种： 降阶法（紧缩法）降阶法（紧缩法） 置大数法置大数法 改改1法法 1. 降阶法降阶法 降阶法也称紧缩法或直接代入法，降阶法也称紧缩法或直接代入法，该法是将结构该法是将结构有限元方程中已知结点位移的自由度全部消去，得到有限元方程中已知结点位移的自由度全部消

43、去，得到一组降阶的修正方程，用以求解其它未知的结点位移。一组降阶的修正方程，用以求解其它未知的结点位移。 如果给定的位移均为零位移 ，则 只需将总刚K、荷载列阵F中与该位移所对应的行和列全部划去即可。 如果给定的位移不为零位移， 也只保留了待定的结点位移作为未知量，但需对右端荷载列阵进行相应的修正。 2. 置大数法置大数法 将结构总刚度矩阵中与被约束的位移分量相对应的主对角线元素赋予一个大数A，如取A=10e30或更大。 再将右端荷载列阵对应的荷载值换成已知的位移值与该大数的乘积。 设结点位移分量r为已知，则有限元方程变为：11112111221222221212rhrhrrrrhhhhhrh

44、hhFkkkkFkkkkkkAkuAkkkkF 经过修改后第r个方程的为1 122rrrrhhkkAkuA(1,2,., ,)riAKihir但ru 方程两边同时除以A，除第r项外，其余各项均为微小量可略去。3. 对角元素改对角元素改1法法 当给定的位移值为零时，将总刚中与之相对应主对角线元素改为1，相对应的行和列中其余所有元素改为0，荷载列阵对应的元素也改为0即可。 111121122122221200000100hhrhhhhhhkkkFkkkFFkkk应力 应变 位移物理方程几何方程外力有限元方程计算模型中：位移场已经确定，就可得到应变、应力。计算模型中：位移场已经确定，就可得到应变、应

45、力。8 应力结果 网格化模型8.1 单元应力计算步骤 有限元方程求解之后，得到了所有结点的位移， 单元应力计算对每个单元循环；对于任一单元 e根据结点i、j、m的实际编号，从结构结点位移向量中选出单元结点位移向量 计算单元的应变分量， Be计算单元的应力分量： D Se8.2 应力分析 以上分析得到了所有单元的应力分量， 为了强度分析，进一步计算主应力或等效应力。xyxy、xx1223()22yyxy主应力取“”号为最大应力，取“”号为最小应力 min180arctan()xyy最大应力与x轴的夹角)()()(21213232221eMISES应力222222()3()exyzxyyzzxxy

46、yzzx 由应力分量表示的三维MISES应力由主应力表示的三维MISES应力由应力分量表示的二维MISES应力2223exyxyxy 8.3 应力显示x应力mises应力 确定根据工程实际情况确定根据工程实际情况确定问题的力学模型确定问题的力学模型，并按一定，并按一定比例绘制结构图、注明尺寸、载荷和约束情况等。比例绘制结构图、注明尺寸、载荷和约束情况等。 将计算对象进行将计算对象进行离散化离散化，即弹性体划分为许多三角形单，即弹性体划分为许多三角形单元，并对节点进行编号。确定全部节点的坐标值，对单元，并对节点进行编号。确定全部节点的坐标值，对单元进行编号，并列出各单元三个节点的节点号。元进行编

47、号，并列出各单元三个节点的节点号。计算载荷的计算载荷的等效节点力。等效节点力。单元分析，由各单元的相关参数，计算单元分析，由各单元的相关参数，计算单元的几何矩阵、单元的几何矩阵、刚度矩阵。刚度矩阵。 组集组集整体刚度矩阵整体刚度矩阵，即形成总刚的非零子矩阵。，即形成总刚的非零子矩阵。 组装各单元的等效结点载荷，形成总的外载荷向量。组装各单元的等效结点载荷，形成总的外载荷向量。 处理约束处理约束，消除刚体位移，消除刚体位移，求解线性方程组求解线性方程组，得到节，得到节点位移。点位移。 计算计算应力矩阵应力矩阵，求得，求得单元应力单元应力，并根据需要计算主应力，并根据需要计算主应力和主方向。和主方

48、向。 整理计算结果（后处理部分）。整理计算结果（后处理部分）。图图1所示为一厚度所示为一厚度t=1cm的均质正方形薄板，上的均质正方形薄板，上下受均匀拉力下受均匀拉力q=106N/m，材料弹性模量为，材料弹性模量为E，泊松比，泊松比，不记自重，试用有限元法求其应力分量。，不记自重，试用有限元法求其应力分量。3/1123421x图 2y2q2q2myxq=106N/m图 1解：解：.力学模型的确定力学模型的确定.结构离散结构离散由于此结构长、宽远大于厚度，而载荷作用于板平面由于此结构长、宽远大于厚度，而载荷作用于板平面内，且沿板厚均匀分布，故可按内，且沿板厚均匀分布，故可按平面应力问题处理平面应

49、力问题处理。考虑到结构和载荷的对称性，可取结构的考虑到结构和载荷的对称性，可取结构的1/4来研究。来研究。该该1/4结构被离散为两个三角结构被离散为两个三角形单元，节点编号，单元划分及形单元，节点编号，单元划分及取坐标如图取坐标如图2所示，其各节点的坐所示，其各节点的坐标值见表标值见表1。节点坐标xy表1123231312110byybyybyy .求单元的刚度矩阵求单元的刚度矩阵1) 计算单元的节点坐标差及单元面积计算单元的节点坐标差及单元面积单元（单元（i、j、m 1，2，3）123231312011cxxcxxcxx 12 33 21111 1 01222b cb c 单元面积单元面积2

50、）组装单元的几何矩阵）组装单元的几何矩阵 00010002101000000101011110ijmijmiijjmmbbbBcccAcbcbcb3）计算单元的应力矩阵）计算单元的应力矩阵2 100001 100001 21 002ijmijmiijjmmSD BbbbEScccAcbcbcb弹性矩阵弹性矩阵210101 1 002ED应力矩阵应力矩阵 srsrsrsrsrsrsrsrsTrrsbbcccbbcbccbccbbEttBDBk21212121142应力矩阵也可应用公式计算应力矩阵也可应用公式计算先计算用到的常数先计算用到的常数89)1 (2169)1 (4312122EEEEt单

51、元的刚度矩阵中各个子矩阵单元的刚度矩阵中各个子矩阵 310011691131000131103110310131003111169111EEK34323234169;031310169;3131311169122113112EKEKEK10031169;1313131169133123EKEK单元单元1的刚度矩阵为的刚度矩阵为: 1031131343131323403131313131031101169133132131123122121113112111166称对EKKKKKKKKKK123123（i、j、m=1，2，3）单元单元2 2：若按：若按i i、j j、m=3m=3、4 4、1 1

52、顺序，对应单元顺序，对应单元1 1的的123123排码排码时时, ,则这两个单元刚度矩阵内容完全一样，故有则这两个单元刚度矩阵内容完全一样，故有: : 1031131343131323403131313131031101169266称对EK3413414. 组集整体刚度矩阵组集整体刚度矩阵 由于Krs=KsrT，又单元1和单元2的节点号按123对应341，则可得： 24424324121331322131122121211188KKKKKKKKKK011016311131631003163113213131112243121233111EKKKEKKKEKKTT按刚度集成法可得整体刚度矩阵为：

53、按刚度集成法可得整体刚度矩阵为：31111630110163300116331111634224163214123241213113231211133123214132244122EKKKEKKKEKKEKKKEKKTTT所以组集的整体刚度矩阵为：所以组集的整体刚度矩阵为： 42110031413001140310240120421141340416388称对EK5.引入约束条件，修改刚度方程并求解引入约束条件，修改刚度方程并求解根据约束条件：u1 =v1=0；v2=0；u4=0和等效节点力列阵： ， 并 代 入 刚 度 方 程 ： ，划去K中与0位移相对应的1，2，4，7的行和列，则刚度方程

54、变为： TqqF2/02/00000 FK2/2/0041104014141634332qqvvuuETTEqEqEqEqvvuu/3/3/4332 TEq1013/103/100/求解上面方程组可得出节点位移为：求解上面方程组可得出节点位移为：所以所以先求出各单元的应力矩阵先求出各单元的应力矩阵S1、S2，然后再求得各单元的，然后再求得各单元的应力分量：应力分量：6.计算各单元应力矩阵，求出各单元应力计算各单元应力矩阵，求出各单元应力 01003/8083/3/03/0001111030310110130383111qqEqEqEqESxyyx 0100)3/(808300/0/3/0111

55、1030310110130383222qEqEEqEqEqESxyyx单元应力可看作是单元形心处的应力值。单元应力可看作是单元形心处的应力值。图中所示为一平面应力问题离散化以后的结构图，图中所示为一平面应力问题离散化以后的结构图，其中图（其中图（a）为离散化后的总体结构，图（）为离散化后的总体结构，图（b）为单元）为单元1，2，3，4的结构，图（的结构，图（c）为单元）为单元3的结构。用有限元法计的结构。用有限元法计算节点位移、单元应变及单元应力（为简便起见，取泊算节点位移、单元应变及单元应力（为简便起见，取泊松比松比，单元厚度，单元厚度t=1）。）。0 xyp1234651234a3ijmi

56、vaiXiujXjuiYjvjYmvmYmXmuaa1,2,4ijmiXiviujXjuiYjvjYmvmYmXmu 图 计算实例2的结构图 各单元节点号与总体节点号对应表各单元节点号与总体节点号对应表 单元号 1 2 3 4 节点号 节 点 总 编 号 I 1 2 2 3 j 2 4 5 5 m 3 5 3 6解：解：1、求确定各单元刚度所需的系数、求确定各单元刚度所需的系数及及面积面积A，对于单元，对于单元1，2，4有：有：mjimjicccbbb,2/0, 02aAcacacababbmjimji对于单元对于单元3有有:2/, 0, 0,2aAacaccabbabmjimji2、求出各单

57、元的单元刚度矩阵。对于、求出各单元的单元刚度矩阵。对于1，2，4单元，单元，其单元刚度矩阵为：其单元刚度矩阵为：10110102020010312112130100202010110144, 2, 1Ekijmijm各单元的节点编号与总体结构的总编号之间的对应关各单元的节点编号与总体结构的总编号之间的对应关系见表系见表2。 31211013011220200011011011011002000243Ek对于单元对于单元3，其单元刚度矩阵为：，其单元刚度矩阵为：ijmijm3、总刚、总刚将各单元刚度矩阵按节点总数及相应的节将各单元刚度矩阵按节点总数及相应的节点号关系扩充成点号关系扩充成12*12矩阵矩阵,分别如下分别如下:654321000000000000000000000000000000000000000000000000000000000000000000000000000000101101000000020200000000103121000000121301000000002020000000101101465432111212Ek65432100000000000000