有限元分析的数学求解原理.ppt

第三章第三章 有限元分析的数学求解原理有限元分析的数学求解原理 前一章针对任意形状变形体，基于物体内的微小体元 dxdydz定义了描述弹性变形体的所有基本力学信息（ui，ij ，ij）、基本方程（平衡、几何、物理）及边界条件。接 下来的任务就是对这些方程在具体的条件下进行求解，也 就是说在已知边界条件下，由基本方程求出相应的位移场 、应力场和应变场。 一般来说，求解方程的途径有两大类：（1）直接针对 原始方程进行求解，方法有：解析法(analytical method)、半 逆解法(semi-inverse method)、有限差分法(finite difference method)等；（2）间接针对原始方程进行求解，方法有： 加权残值法、虚功原理、最小势能原理、变分方法等 主要内容 o 3.1 简单问题的解析求解 o 3.2 虚功原理 o 3.3 应用举例 o 3.4 基本步骤 3.1 简单问题的解析求解 o 1D拉杆问题 有一个左端固定的拉杆，其右端承受一外力P。该拉杆 的长度为l，横街面积为A，弹性模量为E，如图所示： (1)基本变量 由于该问题视为沿x方向的一维问题，因此只有沿x方向 的变量，而其它变量为零。即： 位移：u(x) 应变：x(x) 应力：x(x) 3.1 简单问题的解析求解 (2)基本方程 对原三维问题的所有基本方程进行简化，只保留沿x方 向的方程，得到该问题的三大类基本方程和边界条件 o 平衡方程（无体力） o 几何方程 o 物理方程 3.1 简单问题的解析求解 3.1 简单问题的解析求解 边界条件（BC） 上述方程中，力的边界条件为一种近似，因为在x=l的端面 ，x(x)不应是均匀分布的。由圣维南原理（Saint-Venant principle），在远离x=l的截面，力的边界条件才较好的满 足。 (3)求解 对上述的方程直接求解，可以得到以下的结果： 3.1 简单问题的解析求解 其中c及c1为待定系数，由边界条件可以 求出上式中的常数为c1=0，c=P/A，因此 有最后的结果： (4) 讨论1 上述问题若用经验方法求解（如材料力学的方法），则 需要先作平面假设，即假设x为均匀分布，这样可以得到 再由Hooke定律算出： 再计算右端的伸长量为： 3.1 简单问题的解析求解 通过比较可以看出，经验方法求解的结果和弹性力学的 解析结果完全一致。 比较以上解析方法和经验可以看出： o 解析方法的求解过程严谨，可以得到物体内各点力学变 量的表达，是场变量。 o 经验方法的求解过程比较简单，但需要事先进行假定， 往往只能得到一些特定位置的力学变量表达，而且只能应 用于一些简单情形。 3.1 简单问题的解析求解 (5) 讨论2 根据计算能量的方法，得到： 应变能 外力功 势能 3.1 简单问题的解析求解 o 平面梁的弯曲问题 假设有一个受分布载荷作用的简支梁如图所示，由 于简支梁的厚度较小，外载沿厚度方向无变化，那么该问 题可以认为是一个oxy平面内的问题 3.1 简单问题的解析求解 1.基本方程 有两种方法来建立基本方程。 方法一：采用一般建模及分析方法，即从对象取出 dxdy微元体进行分析，建立最一般的基于（ui，ij，ij）描 述的方程，类似于2D问题的基本变量及方程，这样，所用 的变量较多，方程复杂。 方法二：采用特征建模（characterized modeling）的 简化方法来推导的三大方程，其基本思想是采用工程宏观 特征量来进行问题的描述。 3.1 简单问题的解析求解 应此简支梁问题的特征为：梁为细长梁(long beam)， 因此可只用x坐标来刻画；主要形变为垂直于x的挠度， 可只用挠度（deflection）来描述位移场。 3.1 简单问题的解析求解 补充概念：挠度（ deflection ），弯曲变形时横截面形 心沿与轴线垂直方向的线位移称之为挠度。简言之， 就是指梁、桁架等受弯构件在载荷作用下的最大变形 ，通常指竖向，就是构件的竖向变形 针对这两个特征，可以做出以下假定： l 直法线假定 l 小变形与平面假定 该问题的三类基本变量： 位移： （中层性挠度） 应力：（采用 x，其他应力分量很小，不考虑），该变 量对应于梁截面上的弯矩M 应变：（采用x，满足直线假设） 3.1 简单问题的解析求解 下面取具有全高度梁的dx“微段”来推导三大类方程 3.1 简单问题的解析求解 平衡方程 首先是x方向的合力平衡 然后是y方向的合力平衡 y为距梁中性 层的坐标 最后是弯矩平衡 o 几何方程 由变形后的几何关系，可得到 其中，y为距中性层的坐标，k为梁挠度的曲率，即： 3.1 简单问题的解析求解 物理方程 由Hooke定律有： 对上述方程整理，就得到了平面简支梁弯曲问题的基本 方程： 3.1 简单问题的解析求解 式中， 为梁截面的惯性 矩（moment of inertia）。可 以看出：将原始基本变量定位 中性层的挠度v(x),则可以求出 其他参数 o 边界条件 该简支梁的边界为梁的两端，作用在梁上的q(x)已在 平衡方程中考虑，因此不作为力的边界条件。 两端位移： 两端力（弯矩）： 将弯矩以挠度的二阶导数来表示，即： 3.1 简单问题的解析求解 2.求解 若用基于dxdy微体所建立的原始方程（即原平面应 力问题中的三大类方程）进行直接求解，比较麻烦，并且 很困难，若用基于简化的“特征建模”方法所得到的基本方 程进行直接求解则比较简单，对简支梁问题求解，其方程 为： 3.1 简单问题的解析求解 对上述的常微分方程，其解的形式有： 其中c0c3为待定系数，可由四个边界条件BC求出，有 结果： 3.1 简单问题的解析求解 3.2 虚功原理 o 虚位移与虚功原理如图所示的平衡力系，由于该系统处于 平衡状态，则有： 假想在该平衡力系上作用有微小的扰动 （不影响原平衡条件），且外力所作用 的位置产生了微小的位移变化，即A ，B。该假想的位移如果不影响原平衡 条件，应满足以下几何关系： 进一步分析。当杠杆处于平衡状态时，A和B这两个位移是不存在的 ，但是如果某种原因，例如人为地振一下让它倾斜，一定满足这种关系。 将这个客观存在的关系抽象成一个普遍的原理，去指导分析和计算结 构。 对于在力的作用下处于平衡状态的任何物体，不用考虑它是否真正发 生了位移，而假想它发生了位移，(由于是假想，故称为虚位移)，那么， 物体上所有的力在这个虚位移上的总功必定等于零。这就叫做虚位移原理 ，也称虚功原理。 在图a中的PA和PB所作的功就不是发生在它本身(状态a)的位移上，(因 为它本身是平衡的，不存在位移)，而是在状态(b)的位移上作的功。可见 ，这个位移对于状态(a)来说就是虚位移，亦即是状态(a)假象的位移。 3.2 虚功原理 必须指出，虚功原理的应用范围是有条件的，它所涉及到的两个方 面，力和位移并不是随意的。对于力来讲，它必须是在位移过程中处于 平衡的力系；对于位移来讲，虽然是虚位移，但并不是可以任意发生的 。它必须是和约束条件相符合的微小的刚体位移。 还要注意，当位移是在某个约束条件下发生时，则在该约束力方向 的位移应为零，因而该约束力所作的虚功也应为零。这时该约束力叫做 被动力。(如图中的反力Rc由于支点C没有位移，故Rc所作的虚功对于零 )。反之，如图中的PA和PB是在位移过程中作功的力，称为主动力。因 此，在平衡力系中应当分清楚哪些是主动力，哪些是被动力，而在写虚 功方程时，只有主动力作虚功，而被动力是不作虚功的。 3.2 虚功原理 虚功原理表述如下： 在力的作用下处于平衡状态的体系，当发生与约束条件相 符合的任意微小的刚体位移时，体系上所有的主动力在位移 上所作的总功(各力所作的功的代数和)恒对于零。 虚功原理用公式表示为： 这就是虚功方程，其中P和相应的代表力和虚位移。 3.2 虚功原理 虚功方程是按刚体的情况得出的，即假设图示杠杆是绝对刚性，没有任何 的变形，因而在方程中没有内功项出现，而只有外功项。 将虚功原理用于弹性变形时，总功W要包括外力功(T)和内力功(U)两部分 ，即： W = T - U ；内力功(-U)前面有一负号，是由于弹性体在变形过程 中，内力是克服变形而产生的，所有内力的方向总是与变形的方向相反， 所以内力功取负值。 根据虚功原理，总功等于零得： T - U = 0 外力虚功 T = 内力虚功 U 弹性力学中的虚功原理可表达为：在外力作用下处于平衡状态的弹性体 ，如果发生了虚位移，那么所有的外力在虚位移上的虚功(外力功)等于整 个弹性体内应力在虚应变上的虚功(内力功)。 3.2 虚功原理 虚应变能 虚应变分量 外力虚功 内力虚功即应力在虚应变上做的的虚功，也称虚应变能 外力虚功即作用于弹性体上的外力在虚位移上做的功 由于虚位移是微小的，可认为在虚位移发生过程中外力保持 为常量，则上式的变分符号可提到积分号外。 3.2 虚功原理 最小势能原理： 表明在满足位移边界条件的所有可能位移中，实 际发生的位移使弹性体的势能最小。即对于稳定平衡状态，实际发生的 位移使弹性体总势能取极小值。显然，最小势能原理与虚功原理完全等 价。 最小势能原理 在有限元的理论中，最小势能原理是在所有满足给定边界条件的 位移时，满足平衡微分方程的位移使得势能取得最小值。 最小势能原理就是说当一个体系的势能最小时，系统会处于稳定平衡 状态。或者说在所有几何可能位移中，真实位移使得总势能取最小值 3.2 虚功原理 虚功原理的矩阵表示 i点外力分量 j点外力分量 外力分量用 表示； 引起的应力分量用 表示 3.2 虚功原理 假设发生了虚位移 虚位移分量为 用 表示；引起的虚应变分量 用 表示 虚功原理的矩阵表示 3.2 虚功原理 虚功原理的矩阵表示 在虚位移发生时，外力在虚位移上的虚功是： 式中 是 的转置矩阵。 同样，在虚位移发生时，在弹性体单位体积内，应力在虚应变上的虚 功是： 因此，在整个弹性体内，应力在虚应变上的虚功是： 根据虚功原理得到： 这就是弹性变形体的虚功方程，它通过虚位移和虚应变表明外力与应力 之间的关系。这是以后推导有限元方程的基础。 3.3 应用实例 A1, l1 A2, l2 R3 1. 离散化 1 2 3 2. 位移函数 A, l i j F i ui F j uj x x 2. 位移函数 A, l i j F i ui F j uj （单元内位移线性分布） 形函数矩阵 3.3 应用实例 34 3. 单元刚度矩阵方程 A 虚功原理 外力虚功 虚应变能 应变 应力 应变矩阵 弹性矩阵 3.3 应用实例 单元刚度矩阵 单元的刚度方程 单元刚度矩阵 Element Element 3.3 应用实例 36 B 最小势能定理 外力虚功 虚应变能 由势能变分原理（势能最小原理）得 势能变分，整理得平衡方程 3.3 应用实例 4 整体分析 整体分析就是建立整个离散结构所有节点位移与外力之间 的关系，实现未知节点位移的求解 整体平衡方程 整体刚度方程可基于势能变分原理建立，也可根据节点 的静力平衡来实现（即每个节点静力平衡）。 节点i的平衡为 三个节点三个自由度，即 3.3 应用实例 （A）扩充单元刚度方程法 位移协调性 载荷的叠加性 A1, l1 A2, l2 R3 1 2 3 3.3 应用实例 整体刚度方程 3.3 应用实例 A1, l1 A2, l2 P 1 2 3 （B）“对号入座”法 （方便编程） i j 1 2 2 3 Total 1 2 3 1 2 3 3.3 应用实例 5. 引入边界条件求解 A1, l1 A2, l2 R3 1 2 3 边界条件 支反力 3.3 应用实例 o 结构离散 o 单元分析 o 整体分析 3.4 基本步骤为三大步骤 1、结构离散:就是用假想的线或面将连续物体分割成有 限个单元组成的集合体且单元之间仅在节点处连接，单元 之间的作用仅由节点传递。（基本要求） 注意的问题 o 单元：满足一定几何特性和物理特性的最小结构域 o 节点：单元与单元间的连接点 o 节点力：单元与单元间通过节点的相互作用力 o 节点载荷：作用于节点上的外载 3.4 基本步骤为三大步骤 2单元分析:1），选择插值（位移）函数；2），构造位移 函数。 插值函数：用以表示单元内物理量变化（如位移或位移场 ）的近似函数。由于该近似函数常由单元节点物理量值插值 构成，故称为插值函数，如单元内物理量为位移，则该函数 称为位移函数。 选择位移函数的一般原则 o 位移函数在单元节点的值应等于节点位移（即单元内部是 连续的）； o 所选位移函数必须保证有限元的解收敛于真实解。 位移函数一般采用多项式形式，在单元内选适当阶次的多 项式可得到与真实解接近的近似解 3.4 基本步骤为三大步骤 构造位移函数：如平面问题位移函数的一般形式为 1.多项式项数越多，则逼近真实位移的精度越高，项数的多 少由单元的自由度数决定。 2多项式选取应由低阶到高阶，尽量选择完全多项式以提高单 元精度。 3.选取多项式时，还应使所选取的多项式具有坐标的对称性 ，即按Pascal（帕斯卡）三角形来选择 3.4 基本步骤为三大步骤 位移函数构造方法： 1.广义坐标法： 2.插值函数法：即将位移函数表示为各个节点位移与已知插 值基函数积的和 3.4 基本步骤为三大步骤 3.单元特性分析：单元特性分析的基本任务就是建立单元的 平衡方程，也称为刚度方程。在选择了单元类型和相应的位移 函数后，即可按弹性力学的几何方程、物理方程导出单元应变 与应力的表达式，最后利用虚位移原理或最小势能原理或直接 法或加权残值法建立单元的平衡方程，即单元节点力与节点位 移间的关系。 3.4 基本步骤为三大步骤 3.4 基本步骤为三大步骤 整体分析 整体分析的基本任务包括建立整体平衡方程，引入边 界条件，完成整体方程求解。 整体平衡方程的建立有多种方法，可基于能量原理（ 势能变分或虚位移原理）推导，也可基于节点力平衡得 到。 在引入边界条件之前，整体平衡方程是奇异的，这意 味着整体方程是不可解的。 方程求解包括边界条件引入和数值计算，一旦利用适 当的数值方法求出未知的节点位移，则可按前述的应力 应变公式计算出各个单元的应变、应力等物理量。 3.4 基本步骤为三大步骤 o 刚度由使其产生单位变形所需的外力值来量度，刚度是 指零件在载荷作用下抵抗弹性变形的能力。 o 单元的刚度矩阵：单元刚度矩阵反应的是单元节点力与 单元节点位移的关系； o 总刚度矩阵反应的是整体的节点力与节点位移的关系； 刚度矩阵将总体坐标下的节点位移与整个结构的总体力联 系在一起。 补充实例（刚度矩阵的理解） o 单元刚度矩阵行数等于位移向量的分量个数，列数等于 为位移的列向量的分量个数，由于两者相等所以单刚是个 方阵。 o 结构的总体刚度矩阵即结构的原始刚度矩阵，每1个元 素的物理意义就是当其所在列对应的节点位移分量等于单 位位移（其余结点位移分量为0）时，其所在行对应的节点 力的数值。 表示由于第j个自由度的单位位移dj在第i个自由 度需要的力 补充实例（刚度矩阵的理解） 单元的刚度矩阵的性质 a. 单元刚度矩阵仅与单元的几何特征和材料性质 有关。仅与单元的横截面积A、惯性矩I、单元长度l 、单元的弹性模量E有关。 b. 单元刚度矩阵是一个对称阵。在单元刚度矩阵 对角线两侧对称位置上的两个元素数值相等，即，根 据是反力互等定理。 c. 单元刚度矩阵是一个奇异阵。 d. 单元刚度矩阵可以分块矩阵的形式表示。具有 确定的物理意义。 整体刚度矩阵的性质 整体刚度矩阵 中位于主对角线上的子块 ，称 为主子块，其余 为副子块。 a. 中主子块 由结点i的各相关单元的主子块扩 展之后叠加求得，即 b. 当结点i、 j为单元e的相关结点时， 中副子块 为该单元e相应的副子块，即 。 c. 当结点i、 j为非相关结点时， 中副子块 为 零子块，即 。 d. 仅与各单元的几何特性、材料特性，即A、I 、l、E等因素有关。 e. 为对称方阵， f. 为奇异矩阵，其逆矩阵不存在，因为建立整 体刚度矩阵时没有考虑结构的边界约束条件。 g. 为稀疏矩阵，整体刚度矩阵中的非零元素分布区域的宽度 与结点编号有关，非零元素分布在以对角线为中心的带状区域内 ，称为带状分布规律，见图a。在包括对角线元素在内的区域中， 每行所具有的元素个数叫做把半带宽