第五章 形函数与等参单元ppt课件

1 第五章形函数与等参单元 2 内容回顾 通过前面的学习认识到了形函数在有限单元法中的重要性 形函数作为单元的内插函数 在单元的位移模式中有重要作用 对单元刚度的推导 外力的节点等效起着关键作用 推导原理和过程明确 但是推导繁琐 只能适应简单的少节点单元 常应变三角形单元等 思变 形函数仅是单元内部的位移插值函数 利用节点位移达到内部位移 可考虑利用数学中的 插值函数 方法 直接给出形函数 从而避开繁琐的推导 3 一 面积坐标 对于三角形单元 用面积坐标代替一般的直角坐标 不仅可以简化应力矩阵 刚度矩阵和载荷矩阵等的运算 而且它不随三角形单元形状和方位改变 对于计算机的应用也十分有利 如图所示的三角形单元ijm 任意一点P x y 的位置 可以用如下三个比值来确定 式中为三角形ijm的面积 分别是三角形Pjm Pmi Pij的面积 这三个比值称为P点的面积坐标 节点的对向三角形面积 5 1 5 1面积坐标与自然坐标 1 面积坐标的定义 4 i 三个面积坐标并不完全是独立的 只有两个是独立的 所以 iv 根据面积坐标的定义 不难看出 在平行jm边的直线上的所有各点 都有相同的Li坐标 并且这个坐标就等于 该直线至jm边的距离 与 结点i至jm边的距离 的比值 图中示出了Li的一些等值线 5 2 5 3 2 面积坐标的性质 ii 三角形的三个顶点处 iii 三角形的三条边上 5 3面积坐标和直角坐标之间的关系 于是 5 5a 5 4 类似 5 5b 三角形Pjm的面积是 6 它们的矩阵形式可写为 与 5 4 式对比 可见前述三角形常应变单元中的形函数Ni Nj Nm就等于面积坐标Li Lj Lm 5 6 则三节点三角形单元的位移函数可以写成面积坐标的表达式 7 将式 5 6 求逆 可得到 同理有 以及 以上三式就是面积坐标与直角坐标之间的变换公式 设Li Lj为独立变量 则Lm 1 Li Lj 变换式可以把平面上的任意三角形ijm变换为LiLj平面上的三角形i1 j1 m1 8 4面积坐标的求导和积分 当面积坐标的函数对直角坐标求导时 可以应用下列公式 可写为 9 求面积坐标的幂函数在三角形单元上的积分时 可以应用积分公式 式中为整函数 求面积坐标的幂函数在三角形某一边的积分值时 可以应用积分公式 式中l为该边的长度 可得到 5 7 5 8 10 二 四边形自然坐标 归一化处理 母体单元 实际单元 以两对边中点建立局部坐标系 并且规定坐标值的有效范围0 1 正方形内部的任意一点位置可由其局部坐标来定义 即 四节点任意四边形 按照上述方法 两对边中点 定义局部坐标系 显然在这种局部坐标系下 二者是等同的 二者之间有映射关系 这种基于四边形自然形状而定义的坐标系成为 自然坐标系 二者之间可通过一定的坐标转换公式进行转换 11 优点 该种坐标系下 无论四边形的大小和形状如何 其坐标特征是相同的 因此可用统一的表达式描述 可以推导统一的有限元公式 以局部坐标系导出的公式 有利于数值积分运算 可克服高精度单元的单刚矩阵 等效节点力矩阵等因无法导出显式而必须进行积分所遇到的困难 对高精度单元 单元边界可以使直线 曲线 能更好的逼近实际物体的边界 12 六节点三角形单元具有12个节点自由度 位移函数采用如下的完全二次多项式 5 9 i 1 0 0 j 0 1 0 m 0 0 1 1 0 1 2 1 2 2 1 2 0 1 2 2 1 2 1 2 0 x y 1 位移函数 三六节点三角形单元 13 待定常数a1 a6可由六个节点水平位移分量确定 a7 a12可由六个节点垂直位移分量确定 a1 a2 a3和a7 a8 a9反应了刚体位移和常应变 由于位移函数采用了二次曲线 所以在任意两个单元的相邻边界上 位移分量是按抛物线变化的 抛物线的形状由三个系数确定 而边界上三个节点的水平位移u与垂直位移v完全可以分别确定三个系数值 从而保证了边界上位移的连续性 因此上述位移函数可以满足解答的收敛性条件 为了计算简便起见 把边界上的节点取在三角形三条边的中点 分别以1 2 3表示 如图所示 位移函数可表示为 5 10 14 式中的Ni 和N1是用面积坐标表示的形函数角点 边中点 5 12 5 11 15 利用上述形函数可将位移函数用矩阵形式表示为 式中 5 13 16 利用式 5 11 与式 5 12 并将六个节点的面积坐标依次代入式 5 10 将得出u分别等于ui uj um u1 u2 u3 v分别等于vi vj vm v1 v2 v3 同时由式 5 6 可见 面积坐标与直角坐标是线性相关的 因而式 5 10 中的形状函数即是面积坐标的二次式 也是直角坐标的二次式 都能在六个节点处给出六个节点位移 因此式 5 10 与式 5 9 是等同的 同样都能保证解答的收敛性 17 根据应变公式 2 单元的应变与应力 由于形函数是面积坐标的函数 而面积坐标又是坐标x y的函数 因此上述求导需按照复合函数求导法则进行 18 可导出用节点位移表示单元应变的表达式 式中 19 从而得 5 14 同理可得 5 15 20 5 16 写成矩阵形式为 21 式中 5 17 22 而 23 由上述各式可知 应变分量是面积坐标的一次式 因而也是直角坐标的一次式 故应变是按线性变化的 单元的应力为 式中 5 18 24 而 由于应力矩阵 S 的元素都是面积坐标的一次式 也是直线坐标的一次式 因此单元中的应力是按线性变化的 精度比常应变单元精度高 25 对于六节点三角形单元 节点力与节点位移之间的关系仍然可写成 式中 k 为单元刚度矩阵 它是12 12的方阵 其通式仍为 5 19 3 单元刚度矩阵 26 将式 5 57 转置与式 5 58 相乘 根据式 5 47 应用下列面积坐标的幂函数在单元上的积分公式 i j m 5 20 i j m 5 21 i j m 5 22 27 对其中各元素进行积分 再利用关系式bi bj bm 0和ci cj cm 0加以化简 最后得 5 23 28 式中 i j m i j m 对于平面应变问题 须在应力矩阵及刚度矩阵的各个公式中 将E改换为 将 改换为 29 由于位移函数是非线性的 非节点荷载的移置只能应用普遍公式求载荷列阵 1 分布体积力设单元重力为W 其体积力分量为 移置到各个结点的体积力荷载列阵为 5 24 4 荷载移置 30 利用积分公式 5 8 5 11 5 12 可求得 31 将以上两式代入式 5 64 得载荷列阵为 即Xi Xj Xm X1 X2 X3 0 Yi Yj Ym 0 Y1 Y2 Y3 也就是说 只须向节点1 2 3 而不是i j m 分别移置三分之一的重力 5 25 32 y x i j m 1 2 3 l p 2 分布面力 设单元在边界上作用有沿x方向线性变化的面力 分布的面力在i点的密度为p 在j点为零 如图所示 面积坐标Li在i点为1 Xi Xj Xm X1 X2 X3 0并沿边界按线性变化 面力分量可表示为 33 移置到各个节点的面力荷载列阵为 5 26 由于分布面力仅作用在边界上 所以 N1 N2 Nm 0 而 34 代入式 5 26 并应用面积坐标的幂函数在单元边界的积分公式 5 8 对其中各元素进行积分得 将积分结果代入式 5 26 得 5 27 35 这就是说 只须把总面力的移置到节点i 移置到节点3 根据这一原则 可以应用叠加原理求得边界上受任意线性分布面力的荷载列阵 整个弹性体的平衡方程仍然可取如下形式 式中 K 总体刚度矩阵 结点位移列阵 R 节点荷载列阵 采用六节点三角形单元进行计算 精度比简单三角形单元高 整理计算成果也较简单 但由于节点较多 总体刚度矩阵的带宽较大 所以计算时间稍长 36 5 2等参数单元 一等参单元的基本概念 在进行有限元分析时 单元离散化会带来计算误差 主要采用两种方法来降低单元离散化产生的误差 1 提高单元划分的密度 被称为h方法 h method 2 提高单元位移函数多项式的阶次 被称为p方法 p method 图5 1整体坐标系下四结点矩形单元 同第四章的方法类似 将单元的位移模式选为 5 28 1整体坐标系下 在整体坐标系中设四结点的矩形单元 如图所示 该矩形单元在x及y方向的边长分别为2a和2b 厚度为t 单元编号如图 一 四节点平面矩形单元 37 可得到 5 29 形态函数为 5 30 上述单元位移模式满足位移模式选择的基本要求 反映了单元的刚体位移和常应变 单元在公共边界上位移连续 在矩形单元的边界上 坐标x和y的其中一个取常量 因此在边界上位移是线性分布的 由两个结点上的位移确定 38 现选择如下的局部坐标系 单元局部坐标系的原点取在中心 并以 x y 坐标作为位移插值函数 并将单元坐标函数归一化 采用无量纲坐标 则有局部坐标与整体坐标的关系如下 因此 无论原单元的变长多大 其四个节点的分别为 2单元局部坐标系 a 注意 5 30 和 a 式的关系 39 与三结点三角形单元相比 四结点矩形单元的位移模式是坐标的二次函数 能够提高计算精度 但也有显著的缺点 两种单元的比较如下 如果任意形状的四边形四结点单元仍采用矩形单元的位移模式 则在公共边界上不满足位移连续性条件 为了既能得到较高的计算精度 又能适应复杂的边界形状 可以采用坐标变换 即通过坐标变换将局部坐标中与四边形方向无关的边长为2的正方形 下图b 变换为整体坐标中的四边形 下图a 3平面四节点等参单元 40 图5 2a整体坐标下的任意四结点四边形单元 图5 2b局部坐标下的四结点正方形单元 如上图所示的任意四边形单元上 用等分四条边的两族直线分割四边形 以两族直线的中心为原点 建立局部坐标系 沿和 增大的方向作为 和轴 由于可将局部坐标系下的正方形单元变换为整体坐标系的四边形单元 对于局部坐标下的正方形单元 可取如下的位移函数 a平面四节点等参单元的形函数与坐标变换 41 参照矩形单元 四结点正方形单元的位移模式为 其中 式中是矩形单元四个节点的局部坐标 其坐标值分别为 其中 5 31 5 32 5 34 上式可简记为 形态函数 5 33 母单元 42 由式 5 31 的位移表达式以及式 5 34 可知 在单元的四条边上 分别有 位移线性变化 固边界上的位移呈线性变化 保证了单元公共边界上位移的连续性 因此给出任意四边形单元的结点位移就能得到整个单元上的位移 5 31 的位移模式就是所要找的正确的位移模式 矩形单元的形函数具有与三角形单元形函数相类似的特点 既具有如下关系 在单元的形心上 43 对于母单元的计算 他本身并没有太大的使用价值 但可以利用他得到实际的计算单元 整体坐标下的单元 利用形函数 5 33 作如下的变换 5 35 使得图5 2b中平面上的几个角点分别映射为图5 2a中平面上的四个角点 如果对实际的计算单元的位移函数仍然采用母单元的形函数 即 可以证明他满足完备性和协调性的要求 由于描述单元变形的函数和描述单元几何形状的函数相同 故称计算单元为等参单元 采用等参单元 使我们可以在局部坐标系中的规则单元上进行单元分析 然后在映射笛卡尔坐标系的实际单元上 等参单元同时具有计算精度高和适用性好的特点 是有限元程序中主要采用的单元形式 实际上就是等参变换 44 将位移函数 5 35 代入几何方程 可以得到位移场 b平面四节点等参单元的几何矩阵 式中 45 根据复合函数求导法则 记为 J称为Jacobi矩阵 可得到 46 可以得到应力场 单元刚度矩阵 其子矩阵的计算公式为 47 二四边形八节点等参单元 为了更好地反映物体内的应力变化 适应曲线边界 在弹性力学平面问题的分析中经常使用四边形八节点等参单元 如图所示 由于每条边上增加了一个结点 单元的边是一条二次曲线 可以更好地适应曲线边界 图四边形八结点单元 图八结点基本单元 48 对于等参单元 先在图所示的八结点基本单元上进行分析 八结点单元一共有16个已知的结点位移分量 基本单元中取如下的位移模式 5 36 该位移模式实际上是一个双二次函数 待定系数由结点位移分量确定 在单元的每条边上 局部坐标 位移是局部坐标 的二次函数 完全由边上的三个结点的位移值确定 所以这个位移模式满足位移连续性条件 实际单元内的位移用形函数表示为 5 37 49 其中的形函数为 将形函数归纳为 5 38 50 形函数 在单元的i结点上的值为1 在其它结点上的值均为0 坐标变换式采用如下相似的公式 5 39 代入公式 5 39 可以得到单元345边在整体坐标下的参数方程 5 40 可见在整体坐标系中 单元的边是一条抛物线或退化为一条直线 51 如图所示 ANSYS提供的PLANE82单元是一个四边形八结点等参单元 局部坐标定义为s和t 如图所示 PLANE82单元可以退化为三角形六结点单元 ANSYS提供的Plane82单元 Plane82的基本单元 ANSYS理论手册中给出的PLANE82单元的位移模式如图5 8所示 位移模式与公式 5 18 展开后是一样的 52 三等参单元的单元分析 以平面问题的四边形八结点等参单元为例 介绍构造等参单元的单元刚度矩阵的基本过程 弹性力学平面问题的单元刚度矩阵为 单元的应变为 单元的结点位移 将形函数代入后 可以得到应变的矩阵表达式 5 41 53 可得应变矩阵的分块矩阵 5 42 由于等参单元的形函数是局部坐标 的函数 因此应变矩阵 B 也是局部坐标 的函数 形成等参单元的单元刚度矩阵需要在整体坐标系中对局部坐标的函数进行积分 包括以下三个基本步骤 1 计算用局部坐标表示的形函数 对整体坐标x y的偏导数 2 将整体坐标系中的面积积分转换为在局部坐标系中的面积积分 3 用数值积分计算出单元刚度矩阵中的元素 由于局部坐标与整体坐标之间存在坐标转换关系 因此形函数Ni是局部坐标的函数 同时也可以看作是整体坐标的函数 由复合函数求导法则可得 一 计算形函数对整体坐标x y的偏导数 54 5 43 5 44 定义 或 则有 5 45 矩阵 J 称为雅可比矩阵 JacobianMatrix 单元的整体坐标可以形函数来表示 因此用坐标变换公式可以计算雅可比矩阵 55 雅可比矩阵 5 46 5 47 由 5 45 可得 将 5 48 代入 5 47 得