有限元方法课件 第四章 平面三角形单元.ppt

第四章平面三角形单元 4 1有限元法的基本思想 4 2三角形常应变单元 4 3形函数的性质 4 4刚度矩阵 4 5等效节点力载荷列阵 4 6有限元分析的实施步骤 4 7计算实例 第四章平面三角形单元 一 有限元法的基本思想 假想的把一连续体分割成数目有限的小体 单元 彼此间只在数目有限的指定点 结点 出相互连结 组成一个单元的集合体以代替原来的连续体 再在结点上引进等效力以代替实际作用于单元上的外力 选择一个简单的函数来近似地表示位移分量的分布规律 建立位移和节点力之间的关系 有限元法的实质是 把有无限个自由度的连续体 理想化为只有有限个自由度的单元集合体 使问题简化为适合于数值解法的结构型问题 4 1有限元法的基本思想 二 经典解与有限元解的区别 为平面应力问题 由于结构的对称性可取结构的1 4来研究 故所取的力学模型 三 有限元法算题的基本步骤 1 力学模型的选取 平面问题 平面应变问题 平面应力问题 轴对称问题 空间问题 板 梁 杆或组合体等 对称或反对称等 例如 根据题目的要求 可选择适当的单元把结构离散化 对于平面问题可用三角元 四边元等 2 单元的选取 结构的离散化 例如 结构离散化后 要用单元内结点的位移通过插值来获得单元内各点的位移 在有限元法中 通常都是假定单元的位移模式是多项式 一般来说 单元位移多项式的项数应与单元的自由度数相等 它的阶数至少包含常数项和一次项 至于高次项要选取多少项 则应视单元的类型而定 3 选择单元的位移模式 4 1 单元内任一点的位移列阵 单元的结点位移列阵 单元的形函数矩阵 它的元素是任一点位置坐标的函数 4 单元的力学特性分析 把 4 1 式代入几何方程可推导出用单元结点位移表示的单元应变表达式 4 2 式中 单元内任一点应变列阵 单元的应变矩阵 它的元素仍为位置坐标的函数 再把 4 式代入物理方程 可导出用单元结点位移列阵表示的单元应力表达式 4 3 最后利用弹性体的虚功方程建立单元结点力阵与结点位移列阵之间的关系 即形成单元的刚度方程式 式中 单元内任一点的应力列阵 单元的弹性矩阵 它与材料的特性有关 式中 单元刚度矩阵 4 4 4 5 考虑整体结构的约束情况 修改整体刚度方程之后 4 6 式就变成以结点位移为未知数的代数方程组 解此方程组可求出结点位移 用直接刚度法将单刚组集成总纲 并将组集成总载荷列阵 形成总体结构的刚度方程 4 6 解出整体结构的结点位移列阵后 再根据单元结点的编号找出对应于单元的位移列阵 将代入 4 3 式就可求出各单元的应力分量值 5 建立整体结构的刚度方程 6 求解修改后的整体结构刚度方程 7 由单元的结点位移列阵计算单元应力 求解出整体结构的位移和应力后 可有选择地整理输出某些关键点的位移值和应力值 特别要输出结构的变形图 应力图 应变图 结构仿真变形过程动画图及整体结构的弯矩 剪力图等等 8 计算结果输出 一 离散化 在运用有限单元法分析弹性力学平面问题时 第一步就是要对弹性体进行离散化 把一个连续的弹性体变换为一个离散的结构物 对于平面问题 三角形单元是最简单 也是最常用的单元 在平面应力问题中 单元为三角形板 而在平面应变问题中 则是三棱柱 假设采用三角形单元 把弹性体划分为有限个互不重叠的三角形 这些三角形在其顶点 即节点 处互相连接 组成一个单元集合体 以替代原来的弹性体 同时 将所有作用在单元上的载荷 包括集中载荷 表面载荷和体积载荷 都按虚功等效的原则移置到节点上 成为等效节点载荷 由此便得到了平面问题的有限元计算模型 如图4 1所示 4 2三角形常应变单元 图4 1弹性体和有限元计算模型 图4 2平面三角形单元 二 位移 首先 我们来分析一下三角形单元的力学特性 即建立以单元节点位移表示单元内各点位移的关系式 设单元e的节点编号为i j m 如图4 2所示 由弹性力学平面问题可知 每个节点在其单元平面内的位移可以有两个分量 所以整个三角形单元将有六个节点位移分量 即六个自由度 用列阵可表示为 其中的子矩阵 i j m轮换 a 式中ui vi是节点i在x轴和y轴方向的位移 4 7 从弹性力学平面问题的解析解法中可知 如果弹性体内的位移分量函数已知 则应变分量和应力分量也就确定了 但是 如果只知道弹性体中某几个点的位移分量的值 那么就不能直接求得应变分量和应力分量 因此 在进行有限元分析时 必须先假定一个位移模式 由于在弹性体内 各点的位移变化情况非常复杂 很难在整个弹性体内选取一个恰当的位移函数来表示位移的复杂变化 但是如果将整个区域分割成许多小单元 那么在每个单元的局部范围内就可以采用比较简单的函数来近似地表示单元的真实位移 将各单元的位移式连接 在有限单元法中 虽然是用离散化模型来代替原来的连续体 但每一个单元体仍是一个弹性体 所以在其内部依然是符合弹性力学基本假设的 弹性力学的基本方程在每个单元内部同样适用 起来 便可近似地表示整个区域的真实位移函数 这种化繁为简 联合局部逼近整体的思想 正是有限单元法的绝妙之处 基于上述思想 我们可以选择一个单元位移模式 单元内各点的位移可按此位移模式由单元节点位移通过插值而获得 线性函数是一种最简单的单元位移模式 故设 b 式中 1 2 6是待定常数 因三角形单元共有六个自由度 且位移函数u v在三个节点处的数值应该等于这些点处的位移分量的数值 假设节点i j m的坐标分别为 xi yi xj yj xm ym 代入 b 式 得 c 由 c 式左边的三个方程可以求得 d 其中 4 8 从解析几何可知 式中的 就是三角形i j m的面积 为保证求得的面积为正值 节点i j m的编排次序必须是逆时针方向 如图4 2所示 图4 2平面三角形单元 将 d 式代入 b 式的第一式 经整理后得到 e 其中 同理可得 若令 这样 位移模式 e 和 f 就可以写为 i j m轮换 4 10 i j m轮换 4 9 f 式中I是二阶单位矩阵 Ni Nj Nm是坐标的函数 它们反映了单元的位移状态 所以一般称之为形状函数 简称形函数 矩阵 N 叫做形函数矩阵 三节点三角形单元的形函数是坐标的线性函数 单元中任一条直线发生位移后仍为一条直线 即只要两单元在公共节点处保持位移相等 则公共边线变形后仍为密合 4 11 也可写成矩阵形式 4 12 三 应变 有了单元的位移模式 就可以利用平面问题的几何方程 求得应变分量 将 e f 两式代入上式 即得 g e f 可简写成 其中 B 矩阵叫做单元应变矩阵 可写成分块形式 而子矩阵 由于 和bi bj bm ci cj cm等都是常量 所以矩阵 B 中的诸元素都是常量 因而单元中各点的应变分量也都是常量 通常称这种单元为常应变单元 i j m轮换 4 15 4 14 4 13 四 应力 求得应变之后 再将 4 13 式代入物理方程 便可推导出以节点位移表示的应力 即 4 16 h 4 17 令 则 其中 S 叫做应力矩阵 若写成分块形式 有 对于平面应力问题 弹性矩阵 D 为 4 18 i 所以 S 的子矩阵可记为 i j m轮换 4 19 对于平面应变问题 只要将 i 式中的E换成E 1 2 换成 1 即得到其弹性矩阵 j i j m轮换 4 20 注意到 4 7 式 则有 4 21 由 4 19 4 20 式不难看出 S 中的诸元素都是常量 所以每个单元中的应力分量也是常量 可见 对于常应变单元 由于所选取的位移模式是线性的 因而其相邻单元将具有不同的应力和应变 即在单元的公共边界上应力和应变的值将会有突变 但位移却是连续的 在上节中 提出了形函数的概念 即 其中 i j m轮换 现在我们来讨论一下形函数所具有的一些性质 根据行列式的性质 行列式的任一行 或列 的元素与其相应的代数余子式的乘积之和等于行列式的值 而任一行 或列 的元素与其他行 或列 对应元素的代数余子式乘积之和为零 并注意到 4 9 式中的常数ai bi ci aj bj 4 3形函数的性质 cj和am bm cm分别是行列式的第一行 第二行和第三行各元素的代数余子式 我们有 形函数在各单元节点上的值 具有 本点是1 它点为零 的性质 即 在节点i上 在节点j m上 a b c 类似地有 d 在单元的任一节点上 三个形函数之和等于1 即 e 简记为 4 22 这说明 三个形函数中只有二个是独立的 三角形单元任意一条边上的形函数 仅与该边的两端节点坐标有关 而与其它节点坐标无关 例如 在ij边上 有 4 23 例如 对图4 3所示的单元ijm和ijn 具有公共边ij 这样 不论按哪个单元来计算 根据 4 11 式 公共边ij上的位移均由下式表示 图4 3 由 4 23 式可知 在ij边上 式中Ni Nj的表达形式如 4 23 式所示 i 利用形函数的这一性质可以证明 相邻单元的位移分别进行线性插值之后 在其公共边上将是连续的 由此可见 在公共边上的位移u v将完全由公共边上的两个节点i j的位移所确定 因而相邻单元的位移是保持连续的 为了在以后讨论问题中能够比较方便地确定单元中任意一点处的形函数数值 这里引入面积坐标的概念 在图4 4所示的三角形单元ijm中 任意一点P x y 的位置可以用以下三个比值来确定 图4 4 式中 为三角形单元ijm的面积 i j m分别是三角形Pjm Pmi Pij的面积 这三个比值就叫做P点的面积坐标 4 24 显然这三个面积坐标并不是完全独立的 由于 所以有 而三角形pjm的面积为 故有 类似地有 4 25 4 26 由此可见 前述的三角形常应变单元中的形函数Ni Nj Nm就是面积坐标Li Lj Lm 根据面积坐标的定义 我们不难发现 在平行jm边的直线上的所有各点 都有相同的坐标Li 并且该坐标就等于 该直线至jm边的距离 与 节点i至jm边的距离 之比 图4 4中给出了Li的一些等值线 容易看出 单元三个节点的面积坐标分别为 节点i Li 1Lj 0Lm 0 节点j Li 0Lj 1Lm 0 节点m Li 0Lj 0Lm 1 不难验证 面积坐标与直角坐标之间存在以下变换关系 4 27 当面积坐标的函数对直角坐标求导时 可利用下列公式 4 28 一 单元刚度矩阵 为了推导单元的节点力和节点位移之间的关系 可应用虚位移原理对图4 2中的单元e进行分析 单元e是在等效节点力的作用下处于平衡的 而这种节点力可采用列阵表示为 a 假设在单元e中发生有虚位移 则相应的三个节点i j m的虚位移为 且假设单元内各点的虚位移为 f 并具有与真实位移相同的位移模式 4 4刚度矩阵 故有 c 参照 4 13 式 单元内的虚应变 为 于是 作用在单元体上的外力在虚位移上所做的功可写为 d f 而单元内的应力在虚应变上所做的功为 g 这里我们假定单元的厚度t为常量 把 d 式及 4 16 式代入上式 并将虚位移提到积分号的前面 则有 根据虚位移原理 由 f 和 h 式可得到单元的虚功方程 即 注意到虚位移是任意的 所以等式两边与相乘的项应该相等 即得 记 4 32 则有 4 33 上式就是表征单元的节点力和节点位移之间关系的刚度方程 k e就是单元刚度矩阵 如果单元的材料是均质的 那么矩阵 D 中的元素就是常量 并且对于三角形常应变单元 B 矩阵中的元素也是常量 当单元的厚度也是常量时 所以 4 32 式可以简化为 k e B T D B t 4 34 二整体刚度矩阵 讨论了单元的力学特性之后 就可转入结构的整体分析 假设弹性体被划分为N个单元和n个节点 对每个单元按前述方法进行分析计算 便可得到N组形如 4 33 式的方程 将这些方程集合起来 就可得到表征整个弹性体的平衡关系式 为此 我们先引入整个弹性体的节点位移列阵 2n 1 它是由各节点位移按节点号码以从小到大的顺序排列组成 即 其中子矩阵 j i 1 2 n k 是节点i的位移分量 继而再引入整个弹性体的载荷列阵 R 2n 1 它是移置到节点上的等效节点载荷依节点号码从小到大的顺序排列组成 即 l 其中子矩阵 i 1 2 n m 是节点i上的等效节点载荷 q 同样 将六阶方阵 k 加以扩充 使之成为2n阶的方阵 组装总刚 k 的一般规则 1 当 krs 中r s时 该点被哪几个单元所共有 则总刚子矩阵 krs 就是这几个单元的刚度矩阵子矩阵 krs e的相加 2 当 krs 中rs时 若rs边是组合体的内边 则总体刚度矩阵 krs 就是共用该边的两相邻单元单刚子矩阵 krs e的相加 3 当 krs 中r和s不同属于任何单元时 则总体刚度矩阵 krs 0 下面 我们考查一个组装总刚的实例 1 整体刚度矩阵及载荷列阵的组集 根据叠加原理 整体结构的各个刚度矩阵的元素显然是由有关单元的单元刚度矩阵的元素组集而成的 为了便于理解 现结合图4 5说明组集过程 图中有两种编码 一是节点总码 1 2 3 4 二是节点局部码 是每个单元的三个节点按逆时针方向的顺序各自编码为1 2 3 图中两个单元的局部码与总码的对应关系为 单元e的刚度矩阵分块形式为 47 可编辑 整体刚度矩阵分块形式为 其中每个子块是按照节点总码排列的 通常 采用刚度集成法或直接刚度法来组集整体结构刚度矩阵 刚度集成法分两步进行 第一步 把单元刚度矩阵扩大成单元的贡献矩阵 使单元刚度矩阵的四个子块按总体编号排列 空白处作零子块填充 第二步 把各单元的贡献矩阵对应行和列的子块相叠加 即可得出整体结构的刚度矩阵 如 4 42 式 在这里应该指出 整体刚度矩阵中每个子块为阶矩阵 所以若整体结构分为n个节点 则整体刚度矩阵的阶数是 至于整体结构的节点载荷列阵的组集 只需将各单元的等效节点力列阵扩大成2n行的列阵 然后按各单元的节点位移分量的编号 对应相叠加即可 三整体刚度矩阵的性质 由总刚度方程可知 欲使弹性体的某一节点在坐标轴方向发生单位位移 而其它节点都保持为零的变形状态 在各节点上所需要施加的节点力 刚度矩阵 K 中每一列元素的物理意义为 令节点1在坐标轴x方向的位移u1 1 而其余的节点位移v1 u2 v2 u3 v3 u2n v2n 0 这样就可得到节点载荷列阵等于 K 的第一列元素组成的列阵 即 即表示 是在j节点有单位位移时 而在I节点所需施加的力 s 刚度矩阵 K 中主对角元素总是正的 例如 刚度矩阵 K 中的元素k33是表示节点2在x方向产生单位位移 而其它位移均为零时 在节点2的x方向上必须施加的力 很显然 力的方向应该与位移方向一致 故应为正号 刚度矩阵 K 是一个对称矩阵 即 Krs Ksr T 由 4 32 4 36 式得 所以 可以只存储上三角或下三角矩阵 t 刚度矩阵 K 是一个稀疏矩阵 如果遵守一定的节点编号规则 就可使矩阵的非零元素都集中在主对角线附近呈带状 前面在讨论总刚子矩阵的计算时曾指出 总刚中第r双行的子矩阵 Krs 有很多位置上的元素都等于零 只有当第二个下标s等于r或者s与r同属于一个单元的节点号码时才不为零 这就说明 在第r双行中非零子矩阵的块数 应该等于节点r周围直接相邻的节点数目加一 可见 K 的元素一般都不是填满的 而是呈稀疏状 带状 以图4 6a所示的单元网格为例 其整体刚度矩阵中的非零子块 每个子块为2行2列 的分布情况如图4 6b所示 图4 6a 图4 6b 半带宽B 相邻节点号的最大差值D 1 2 若第r双行的第一个非零元素子矩阵是 Krl 则从 Krl 到 Krr 共有 r l 1 个子矩阵 于是 K 的第2r行从第一个非零元素到对角元共有2 r l 1 个元素 显然 带状刚度矩阵的带宽取决于单元网格中相邻节点号码的最大差值D 我们把半个斜带形区域中各行所具有的非零元素的最大个数叫做刚度矩阵的半带宽 包括主对角元 用B表示 即B 2 D 1 通常的有限元程序 一般都利用刚度矩阵的对称和稀疏带状的特点 在计算求解中 只存储上半带的元素 即所谓的半带存储 因此 在划分完有限元网格进行节点编号时 要采用合理的编码方式 使同一单元中相邻两节点的号码差尽可能小 以便节省存储空间 提高计算效率 刚度矩阵 K 是一个奇异矩阵 在排除刚体位移后 它是正定阵 弹性体在 R 的作用下处于平衡 R 的分量应该满足三个静力平衡方程 这反映在整体刚度矩阵 K 中就意味着存在三个线性相关的行或列 所以 K 是个奇异阵 不存在逆矩阵 在上节讨论整体刚度矩阵时已经指出载荷列阵 R 是由弹性体的全部单元的等效节点力集合而成 而其中单元的等效节点力 R e则是由作用在单元上的集中力 表面力和体积力分别移置到节点上 再逐点加以合成求得 根据虚位移原理 等效节点力的大小 应按其所做的功与作用在单元上的三种力在任何虚位移上所做的功相等这一原则来确定 即 上式中等号的左边表示单元的等效节点力 R e所做的虚功 等号右边的第一项是集中力 G 所做的虚功 第二项的积分是沿着单元的边界进行 表示面力 q 所做的虚功 第三项的积分则是遍及整个单元 表示体积力 p 所做的虚功 t为单元的厚度 假定为常量 a 4 5等效节点力载荷矩阵 一 体力的移置 如果任意三角形单元ijk的重心c上受有自重则按刚体静力等效原理可把W直接移置到i j m三个节点上而组成 如果单元的重心c受有惯性力Pu作用 且 则Pu移置到i j m节点上的等效结点力为 式中 旋转角速度 是单元重心处x坐标 二 面力的移置 已知在ij边受有面力q 则移置到i j结点上的等效节点力为 当某一边上有三角形分布的面力时 可由刚体静力等效直接写出 三 集中力的移置 如集中力G做用于其一边界上如图 先将G分解为 后 分别按线段的比例把和分别移置到i j两点上 即 即按静力平衡方法分配 根据前面的讨论 现以三角形常应变单元为例来说明应用有限元法求解弹性力学平面问题的具体步骤 确定力学模型 根据工程实际情况确定问题的力学模型 并按一定比例绘制结构图 注明尺寸 载荷和约束情况等 将计算对象进行离散化 即弹性体划分为许多三角形单元 并对节点进行编号 确定全部节点的坐标值 对单元进行编号 并列出各单元三个节点的节点号 计算载荷的等效节点力 要求的输入信息 由各单元的常数bi ci bj cj bm cm及行列式2 计算单元刚度矩阵 组集整体刚度矩阵 即形成总刚的非零子矩阵 处理约束 消除刚体位移 4 6有限元分析的实施步骤 求解线性方程组 得到节点位移 计算应力矩阵 求得单元应力 并根据需要计算主应力和主方向 整理计算结果 后处理部分 为了提高有限元分析计算的效率 达到一定的精度 应该注意以下几个方面 一 对称性的利用 在划分单元之前 有必要先研究一下计算对象的对称或反对称的情况 以便确定是取整个物体 还是部分物体作为计算模型 例如 图4 11 a 所示受纯弯曲的梁 其结构对于x y轴都是几何对称的 而所受的载荷则是对于y轴对称 对于x轴反对称 可知 梁的应力和变形也将具有同样的对称特性 所以只需取1 4梁进行计算即可 取分离体如图4 11 b 所示 对于其它部分结构对此分离体的影响 可以作相应的处理 即对处于y轴对称面内各节点的x方向位移都设置为零 而对于在x轴反对称面上的各节点的x方向位移也都设置为零 这些条件就等价于在图4 11 b 中相应节点位置处施加约束 图中o点y方向施加的约束是为了消除刚体位移 图4 11 节点的布置是与单元的划分互相联系的 通常 集中载荷的作用点 分布载荷强度的突变点 分布载荷与自由边界的分界点 支承点等都应该取为节点 并且 当物体是由不同的材料组成时 厚度不同或材料不同的部分 也应该划分为不同的单元 节点的多少及其分布的疏密程度 即单元的大小 一般要根据所要求的计算精度等方面来综合考虑 从计算结果的精度上讲 当然是单元越小越好 但计算所需要的时间也要大大增加 另外 在微机上进行有限元分析时 还要考虑计算机的容量 因此 在保证计算精度的前提下 应力求采用较少的单元 为了减少单元 在划分单元时 对于应力变化梯度较大的部位单元可小一些 而在应力变化比较平缓的区域可以划分得粗一些 二 节点的选择及单元的划分 a b 图4 12 还有一点值得注意的是 单元各边的长度不要相差太大 以免出现过大的计算误差或出现病态矩阵 例如 图4 12所示的 a b 两种单元划分 虽然都是同样的四个节点 但 a 的划分方式显然要比 b 的方式好 在进行节点编号时 应该注意要尽量使同一单元的相邻节点的号码差尽可能地小 以便最大限度地缩小刚度矩阵的带宽 节省存储 提高计算效率 如前所述 平面问题的半带宽为B 2 d 1 三 节点的编号 若采取带宽压缩存储 则整体刚度矩阵的存储量N最多为N 2nB 4n d 1 其中d为相邻节点的最大差值 n为节点总数 例如在图3 13中 a 与 b 的单元划分相同 且节点总数都等于14 但两者的节点编号方式却完全不同 a 是按长边进行编号 d 7 N 488 而 b 是按短边进行编号 d 2 N 168 显然 b 的编号方式可比 a 的编号方式节省280个存储单元 a b 图4 13 四 单元节点i j m的次序 在前面章节中 我们曾指出 为了在计算中保证单元的面积 不会出现负值 节点i j m的编号次序必须是逆时针方向 事实上 节点i j m的编号次序是可以任意安排的 只要在计算刚度矩阵的各元素时 对 取绝对值 即可得到正确的计算结果 在实际计算时 应该注意所选有限元分析软件的使用要求 五 边界条件的处理及整体刚度矩阵的修正 在前面讨论整体刚度矩阵时 已经提到 整体刚度矩阵的奇异性可以采用提高考虑边界约束条件的方法来排除弹性体的刚体位移 以达到求解的目的 一般情况下 求解的问题 其边界往往已有一点的