1、第六章第六章 用有限元法解平面问题用有限元法解平面问题胡 衡武汉大学土木建筑工程学院弹 性 力 学 及 有 限 元二零零八二零零八年六月年六月有限元法的概念有限单元法的目的是求解连续体的力学问题，它将连续体变换为由有限个有限大小的单元组成的离散化结构：这样的离散化结构与桁架相比，其区别在于桁架的单元是杆件，而这里的单元是三角形块体。块体与块体之间只在结点处相连接，弹性体受的面力和体力都按静力等效移植到结点上。有限元法的基本量 有限元法是二十世纪五十年代以来随着计算机的广泛应用发展起来的一种数值解法。 在有限元法中，为了简洁清晰地表明各基本量之间的关系，也为了便于程序编写，广泛地采用矩阵表示和矩
2、阵运算，因此，我们先来介绍一些基本量和基本方程的矩阵表示:体力列阵：面力列阵：有限元法的基本量应力列阵：应变列阵：位移列阵：运用几何方程的应变列阵：平面应力状态下的物理方程的矩阵表示：有限元法的量对于平面应变问题，在D中做如下变换：换为换为有限元法的基本方程有限元法的基本方程在有限单元法中，作用于弹性体的各种外力常以作用于某些点的等效集中力来代替。这些点上的集中力以及它们相应的虚位移可用列阵表示为：于是，外力在虚位移上的虚功为：则集中力下的虚功方程为：有限元法的求解步骤(1)（1）取有限单元的结点位移为基本未知量，如对三角形三结点有限元，它们是：ijm（2）应用插值公式，由单元的结点位移求出单
3、元的位移函数，即求出关系式：这种插值公式表示单元中的位移分布形式，在有限元中称为位移模式，其中N为形函数矩阵。有限元法的求解步骤(2)（3）应用几何方程，由单元的位移模式求出单元的应变，即求出关系式：（4）应用物理方程，由单元的应变求出单元的应力：（5）应用虚功方程，由单元的应力求出单元的结点力。ijm结点对单元的作用力称为结点力，在三角形单元中，结点力列阵为：有限元法的求解步骤(3)（5）应用虚功方程，由单元的应力求出单元的结点力。ijm集中力下的虚功方程为：单元劲度矩阵有限元法的求解步骤(4)（6）将单元中的各种外力载荷向结点转移，化为结点载荷：ijm（7）列出各结点的平衡方程，组成整个结
4、构的平衡方程组。例如，由于结点i受有环绕它的单元转移过来的结点载荷和结点力，因此结点i的平衡方程为： 表示对那些环绕i结点的单元求和。i=1,2,3nn表示应列平衡方程的结点数。有限元法的求解步骤(5)（8）将代入下式：得到：这里 称为整体劲度矩阵， 是整体结点载荷列阵， 是整体结点位移列阵。位移模式单元的位移模式单元的应变列阵和应力列阵（1）单元的应变列阵和应力列阵（2）其中：i, j, m可交替轮换单元的应变列阵和应力列阵（3）单元的应变列阵和应力列阵（4）在三结点三角形单元中，当位移函数取为线性位移模式时，实际上就是将位移函数在单元中用泰勒级数展开并忽略 ， 二阶以上的项得到的结果。再经
5、过求导运算得出的应变和应力都是常量。由此可见，线性位移的误差是 , 的二阶小量，而应力和应变则是一阶小量。故应力的精度低于位移的精度。为了提高计算精度，一般可以采取两种方法，一是将单元的尺寸减小，而是采用包含更高此项的位移模式。单元的结点力列阵与劲度矩阵（1）：单元劲度矩阵单元的结点力列阵与劲度矩阵（2）单元的结点力列阵与劲度矩阵（3）单元刚度矩阵的特点：1. 单元刚度矩阵取决于单元的形状和弹性常数，与单元位置无关，与单元的大小无关。2. 单元刚度矩阵是对称矩阵。单元的结点力列阵与劲度矩阵（4）- 举例xyijmaai, j, m可交替轮换单元的结点力列阵与劲度矩阵（5）xyijmaa载荷向结
6、点移置（1）xyijmaa对于变形体，包括弹性体在内，所谓静力等效，是指原载荷与结点载荷在任何虚位移上的虚功相等。原载荷与结点载荷向任意点简化时，它们具有相同的主矢量和主矩。假设各结点的虚位移为：假设将集中力移置各结点上后，各结点上的结点载荷为：载荷向结点移置（2）xyijmaa则作用于P点的虚位移为：按照静力等效原则，结点载荷在结点虚位移上的虚功，应当等于集中载荷在其作用点上的虚功：载荷向结点移置（3）此为三结点单元上的集中载荷向结点移置的计算公式，其中的各项N为形函数在集中载荷作用处的数值。xyijmaa载荷向结点移置（4）设单元受分布体力 ，可将微分体积 上的体力 当作集中力，利用集中载
7、荷的计算公式，将体力也移至结点：例如，设单元 ijm 的密度为 ，试求自重的等效结点载荷。载荷向结点移置（5）移置到各结点的载荷均为1/3自重。设单元体受分布面力 ，可将微分面积 上的面力 当作集中力，利用集中载荷的计算公式，将面力也移置结点：载荷向结点移置（6）例如，设单元 ij 边上受有x方向上的均布面力q，试求等效结点载荷载荷向结点移置（7）结构整体分析（1）对于每个单元，我们已经知道了如何计算单元的劲度矩阵以及载荷列阵：结构整体分析（2）根据虚功原理，我们也推导了结点力与结点位移的关系：对于 i 点, 一个单元上的结点力为:i 点的力平衡要求围绕 i 点的各单元产生的结点力与各单元分配
8、到 i 点的结点载荷相等。结构整体分析（3）综上所述，i 点的平衡方程为。将所有结点的平衡方程写在一起，就得到了整个结构的结点平衡方程：这里 称为整体劲度矩阵， 是整体结点载荷列阵， 是整体结点位移列阵。结构整体分析（4）其中：n为结点数, 下标1, 2, 3是结点在整体结构中的编号。注：i, j, m 是结点在单元中的局部编号。整体刚度矩阵的一个元素Krs 就是由按整体编码的，同下标的单元劲度矩阵元素的叠加得到的：有限元计算简例(1)xy2m2m2m2m2N/m2N/m123456有限元计算简例(2) 单元刚度矩阵的意义i-3Ij-1m-212356IIIIIIIV4kij 表示由于j结点沿
9、x或y方向上的位移而引起的在i结点的沿x或y方向结点力。有限元计算简例(3)iIIIIVjmIIIijm12356IIIIIIIV4IIIIIIIV结构整体分析（4）整体刚度矩阵中的一个元素Krs表示，由于s结点的位移而引起的在r结点的结点力，下面举例说明：结构整体分析（5）iIIIIVjmIIIijm12356IIIIIIIV4整体刚度矩阵中的K23, 表示3号结点的位移由于结构劲度在2号结点上引起的结点力，因此，它将由两部分组成，一部分来自I单元，一部分来自III单元。观察以下表格，可知 。IIII结构整体分析（6）整体劲度矩阵的获得过程：首先将整体刚度全部充零，然后逐个单元的建立单元劲度
10、矩阵，最后根据结点的局部编码与整体编码的关系，将所有单元劲度矩阵的每一个子矩阵都叠加到整体劲度矩阵中去。结构整体分析（7）- 计算单元刚度矩阵E结构整体分析（8）- 得到整体刚度矩阵E结构整体分析（9）- 在整体刚度矩阵中引入边界条件12356IIIIIIIV4需求解的结点还剩：因此关于这六个零分量的六个平衡方程不用建立，须将整体刚度矩阵的第1，3，7，8，10，12以及同序列的各列去掉。最后得到：结构整体分析（10）- 结点载荷12356IIIIIIIV4iIIIIVjmIIIijm不考虑位移边界条件的情况下：1N/m结构整体分析（11）12356IIIIIIIV4iIIIIVjmIIIij
11、m1N/m考虑边界条件，去掉零位移对应的结点载荷：至此，结构的整体结构刚度矩阵以及整体载荷列阵都已经得到，下面根据平衡方程 求解结点位移。结构整体分析（12）- 求解结构整体分析（13）- 整理结果求应力求应力所需关系式如下：结构整体分析（14）- 整理结果求应力12356IIIIIIIV4iIIIIVjmIIIijmI, II，IVII结构整体分析（14）- 整理结果求应力12356IIIIIIIV4对于单元III：IIIijm结构整体分析（15）1. 有限元法的求解步骤：2. 划分有限元，3. 利用已知的结点坐标以及结构的物理特性写出单元劲度矩阵，4. 利用整体编码与局部编码的关系写出整体
12、刚度矩阵以及力列阵，5. 在整体刚度矩阵以及力列阵中将对应于零位移的行与列划去，得到引入边界条件后的平衡方程组。6. 求解平衡方程组，得到结点位移，并由此分析应力分布。有限单元法的单元划分（1）有限单元的划分需要注意以下几点： 对于结构的不同部位，可以根据需要，采用不同大小的单元。 应使三角形单元的内角大致相当，因为应力与位移的误差与最小内角的正玄成反比。 应该充分利用结构特性，减少计算工作量。 在结构的厚度或弹性常数有突变处，除了应该采用较小的单元外，还应该把突变线作为单元的界限。 如果结构受有集度突变的载荷，也应该在突变部位采用较小的有限元，并在突变处设置结点，使应力的突变有所反应。有限单
13、元法的单元划分（2）l当结构具有凹槽或孔洞时，为了正确地描述应力集中效应，必须把该处的网格画得很密。l当计算容量不允许时，可以分两次计算。第一次计算时，将需要细化网格的目标区域的网格画得稀疏一点，甚至和其他区域的网格大致相同，第二次计算时，将需要细化的部分区域（区域边界上的结点位移是第一次计算后的已知值）取出，利用第一次计算的计算结果，就可以计算分析网格很密的目标区域了。计算结果的整理（1）为了得到结构内某一点的应力，必须通过平均计算，通常可采用绕结点平均法或二单元平均法。所谓绕结点平均法就是把环绕某一结点的各单元的常量应力加以平均，用来表征该结点处的应力。12356IIIIIIIV4a所谓二单元平均法就是把两个相邻单元的常量应力加以平均，用来表征公共边中点的应力。计算结果的整理（2）1234在结构有应力集中的区域，应力往往数值较大且变化剧烈，在推导最大应力时，必须充分注意如何达到最高精度。设已知上图凹槽边界处的结点应力已知，有两种方法可以推导边界上的最大应力：第一种：手绘应力曲线，如果曲线变化不是很强烈，则直接可以观察出最大应力所在处；第二种，利用插值公式，求得应力的函数表达式，并利用该表达式求得最大应力及其所在坐标。
