有限单元法的基本原理.ppt

2.有限单元法的基本原理2.1 虚位移原理所谓虚位移可以是无限小的位移，它在结构内部必须是连续的，在结构的边界上必须满足运动学边界条件，例如对悬臂梁来说，在固定端处，虚位移及其斜率必须等于零外力在虚位移上所做的虚功单位体积内的虚应变能整个物体的的虚应变能虚位移原理：如果在虚位移发生之前，物体处于平衡状态，那么在虚位移发生时，外力所做的虚功等于物体的虚应变能2.2 变分原理泛函 如果对于某一类函数y(x)中的每一个函数y(x)，变量 有一个值和它对应，则变量 称为依赖于函数y(x)的泛函记为变分法就是研究泛函的极大值和极小值的方法 如图在xy平面内连接A、B两点的任一曲线的长度为因此，长度L就是函数y(x)的泛函只要积分的上下限保持不变，变分的运算与定积分的运算可以交换次序一般泛函定义泛函的变分泛函的极值问题——变分问题如果泛函 在 的邻近任意一根曲线上的值都不大于或都不小于 即则称泛函 在曲线 达到极大值或极小值，而必要的极值条件为例2.3 弹性力学平面问题 连续介质的离散连续介质的离散 对于二维连续介质，以图所示的建筑在岩石基础上的支墩坝为例，用有限单元法进行分析的步骤如下： (1)用虚拟的直线把原介质分割成有限个三角形单元，这些直线是单元的边界，几条直线的交点称为结点。

(2)假定各单元在结点上互相铰接，结点位移是基本的未知量 (3)选择位移函数 (4)通过位移函数，用结点位移唯一地表示单元内任一点的应变；再利用广义虎克定律，用结点位移可唯一地表示单元内任一点的应力 (5)利用能量原理，找到与单元内部应力状态等效的结点力，再利用单元应力与结点位移的关系，建立等效结点力与结点位移的关系 (6)将每一单元所承受的荷载，按静力等效原则移置到结点上 (7)在每一结点建立用结点位移表示的静力平衡方程，得到一个线性方程组：解出这个方程组，求出结点位移，然后可求得每个单元的应力 连续介质的有限单元分析包含三个基本方面：介质的离散化、单元特性计算以及单元组合体的结构分析位移函数位移函数对三角形单元，假定单元内的位移分量是坐标的线性函数形函数位移模式位移模式需满足以下三个条件：1位移模式必须反映单元的刚体位移2位移模式必须反映单元的常量应变3位移模式应尽可能反映位移的连续性单元应变单元应变 （几何方程）应变分量是常量单元应力单元应力 （物理方程）单元刚度矩阵单元刚度矩阵结点力和结点位移的关系平面应力问题平面应变问题等效结点力等效结点力静力等效原则：指原荷载与结点荷载在任何虚位移上的虚功都相等。

在一定的位移模式下这样的移置结果是唯一的，而且总能符合通常理解的对刚体而言的静力等效原则分布边界力的等效结点荷载ij边上均布力pxij边上三角形荷载px分布体积力的等效结点荷载结点平衡方程与整体刚度矩阵对单元e，所受结点力为结点i受单元e的力为Ui，Vi，环绕i结点其他单元一起所施加的力结点i从周围各单元移置的结点荷载为以 代入整体刚度矩阵的集成 结点平衡法 按整体编码表示为：直接刚度法 把每个单元的单刚阶数扩大为整体刚度矩阵阶数把单刚中按局部编码的子块搬到整体刚度矩阵中整体编码的位置中去，余下的部分用零子块填充 把各单元刚度矩阵的贡献矩阵叠加在一起 整体刚度矩阵的特点整体刚度矩阵的特点 整体刚度矩阵的任一元素的物理意义是：结构第个结点位移为单位个结点位移方向上施加的结点力的大小 值而其它结点位移为零时，需在第①整体刚度矩阵具有对称性②整体刚度矩阵具有稀疏性 ③整体刚度矩阵具有带状性④整体刚度矩阵是奇异矩阵 边界条件的处理 边界的约束情况 (1)基础支承结构 (2)具有对称轴的结构（3）具有给定位移边界的结构 边界条件的处理方法 (1)直接代入法 按结点位移已知和待定重新组合方程 对角元素改1法 只能用于给定零位移。

对角元素乘大数法 输入离散模型数据按选择的单元计算单元刚度矩阵按总刚存储模式集成总刚按单元循环计算单元等效结点荷载集成结点荷载列阵引入位移边界条件解方程组其他辅助计算结果输出、结束形成K形成P2.4 弹性力学空间问题 四面体单元位移模式位移模式常应变四面体单元2.4 弹性力学空间问题 右手坐标系中，当按照i→j→m的方向转动时，右手螺旋应向p的方向前进常应变四面体单元单元应变单元应变单元应力单元应力单元刚度矩阵单元刚度矩阵结点荷载结点荷载体积力面力 设单元e是靠近边界的，它的某一边界表面ijm，承受线性分布面力在结点i、j、m上的集度分布为2.4 弹性力学轴对称问题对于轴对称问题，采用圆柱坐标(r，θ，z)较为方便如果以弹性体的对称轴作为z轴，所有应力、应变和位移都与θ无关，只是r和z的函数任一点只有2个位移分量，即沿r方向的径向位移u和沿z方向的轴向位移w由于对称，θ方向的环向位移等于零 在轴对称问题中，采用的单元是一些圆环这些圆环和 rz平面正交的截面通常取为三角形，如图所示的ijm各单元之间用圆环形的铰链互相连接，每一个铰与rz平面的交点称为结点。

但是在轴对称问题中，每个单元的体积都是一个圆环的体积，这点与平面问题是不同的 由于对称，只须取出一个截面进行分析，但在计算中应注意到所采用的单元是圆环，所有结点力和结点荷载都是施加在圆环形的铰上 如果弹性体的几何形状是轴对称的，但荷载不是轴对称的，我们可以把荷载在θ方向展成富氏级数，然后分解为轴对称及反轴对称问题求解，即把一个三维问题分解为一组二维问题求解如果弹性体的几何形状、约束条件及荷载都对称于某一轴，例如z轴，则所有的位移、应变及应力也对称于此轴这种问题称为轴对称应力问题在竖井、压力容器及机械制造中，经常遇到轴对称应力问题位移模式单元应变单元应力单元刚度矩阵各向同性体结点荷载对于轴对称问题，结点荷载是作用在整圈圆环形铰上的如设结点的半径为r，单位长度的铰上作用的荷载为 （径向）和 （轴向），则计算中采用的荷载应为 （径向）和 （轴向）结点力移置的一般公式体积力表面力（ij边r方向）惯性力2.5 单元和形函数2.5.1 形函数定义2.5.2 一维形函数形函数是用自然坐标在母单元中定义的一维形函数 -1≤ξ≤1二次单元（3结点）2.5.3 二维形函数二维母单元是（ξ，η）平面中的2×2的正方形，其中 -1≤ξ≤1，-1≤η≤1线性单元（4结点）令二次单元（8结点）在结点i ,ξ0＝η0＝1，因此Ni＝1，而在其他点，Ni＝0，满足条件（1）。

在单元的4条边上，形函数是二次函数，而每边有3个结点，足以保证用形函数定义的未知量在相邻单元的连续性，故满足条件（2）把形函数展开，Ni中包括了线性项ξ和η，这些形函数的线性组合可以充分反映用形函数定义的未知量的任意线性变化，从而满足条件（3）形函数验证 (8结点二次单元）2.5.4 三维形函数三维母单元是（ξ,η, ζ）平面中的2×2×2的正六面体，其中 -1≤ξ≤1，-1≤η≤1， -1≤ ζ ≤1线性单元（8结点）二次单元（20结点）角点典型边中点2.5.5 坐标变换通过进行坐标变换，使(ξ，η，ζ)坐标系中形状简单的母单元，在(x，y，z)坐标系中变换为具有曲线(面)边界的形状复杂的单元，变换后的单元称为子单元子单元在几何上可以适应各种实际结构的复杂外形经过这样处理，单元具有双重特性：一方面，子单元的几何特征、荷载等等，都来自实际结构，充分反映了实际情况，另一方面，大量计算工作是在母单元内进行的，由于它的形状简单而且规则，计算比较方便，并便于循环，特别有利于在电子计算机上进行计算因此兼有两方面的优点平面坐标变换平面坐标变换二维线性单元坐标变换公式为直线24的方程形心坐标子单元的4条边都是二次曲线，局部坐标系（ξ，η）是曲线坐标空间坐标变换空间坐标变换经过空间坐标变换后，原来的直线将变成空间曲线，原来的平面将变成空间曲面。

母单元正六面体，将变为具有曲棱、曲面的六面体子单元例 相邻单元公共边连续性验证2.5.6 位移函数单元位移用位移函数表示如下：如果单元坐标变换和位移函数所采用形函数的阶次相等，那么用以规定单元形状的结点数应等于用以规定单元位移的结点数，这种单元称为等参单元如果坐标变换所用形函数的阶次高于位移函数中的形函数，坐标变换的结点数应超过用以规定单元位移的结点数，这种单元称为超参数单元反之，如果坐标变换所用形函数的阶次低于位移函数中形函数的阶次，则称为逊参数单元写成矩阵形式等参单元的位移函数满足刚体位移和常应变条件，满足连续性条件，也满足收敛条件空间单元的位移函数为2.5.7 单元应变 空间问题的应变可表示为：代入位移函数 形函数 是用局部坐标表示的，根据偏微分法则，有雅可比矩阵计算示例单元a单元b单元c2.5.8 刚度矩阵单元刚度矩阵的一般公式令将整体坐标系中的积分转换成局部坐标系中积分（注意是矢量乘积）2.5.9 结点荷载体积力 产生的结点荷载为作用于单元边界 上的分布力 所产生的结点荷载按下式计算式中：[N]为表面 Ω 的形函数矩阵，是3×3s阶矩阵，s为单元表面Ω 的结点数设 所作用的表面Ω是ζ＝1的面，在坐标变换公式中令ζ＝1可得到Ω的方程表面Ω上任一点的压强为设表面Ω任一点的法线的方向余弦为l，m，n，则表面力在各方向的分力为2.5.10 等参数单元的退化4结点平面等参数单元的退化结点平面等参数单元的退化令8结点空间等参单元的退化结点空间等参单元的退化高次单元的退化高次单元的退化高次单元退化后必须修改形函数高次单元退化后必须修改形函数2.5.11 数值积分在求解刚度矩阵和结点荷载时，需计算如 的积分。

但一般是很复杂的，通常难以用显式表示其积分，一般都用数值积分方法计算积分值，即在单元内选出某些点，称为积分点，求出被积函数 在这些点的值，然后根据这些数值求出积分值数值积分有两类方法，一类方法积分点是等间距的，如辛普生方法；另一类方法积分点是不等间距的，如高斯方法一维高斯积分公式一维高斯积分公式±ξiHi0.5773502692n＝21.00000000000.77459666920.0000000000n＝30.55555555560.88888888890.86113631160.3399810436n＝40.34785484510.6521451549 和 是根据计算精度最高而选定的，积分点 应是勒让德多项式 的根加权系数 按下式计算二维及三维高斯积分公式二维及三维高斯积分公式先令 保持常数，计算沿 方向的积分再沿 方向积分对三重积分有一般采用2×2×2高斯积分2.6 非线性有限元分析方法线性弹性力学采用两个基本假定：1材料的应力应变关系是线性的，即假定材料符合胡克定律2应变位移关系是线性的，即小位移假定。

例如当钢材的应力超过其比例极限后，应力应变关系便是非线性的又如土壤和岩石的应力应变关系也是非线性的这些称为材料非线性又如梁、板及薄壳等结构失稳后，由于产生了大位移，其应变位移关系是非线性的，这些称为几何非线性在热传导问题中，某些情况下材料的导温系数及内部热源与温度有关在流体力学中，粘滞系数与流速有关，或者由于出现紊动，达西定律不再适用这些问题都是非线性的当材料的应力应变关系是非线性的时，刚度矩阵不是常数，而与应变和变位值有关可记为 这时结构的整体平衡方程是如下的非线性方程组：（1）2.6.1 增量法 采用增量法分析非线性问题时，把荷载划分为许多荷载增量，这些增量可以相等，也可以不等每次施加一个荷载增量在每一步计算中，假定方程是线性的，刚度矩阵是常数，在不同的荷载增量中，刚度矩阵可以具有不同的数值每步施加一个荷载增量{ΔP}，得到一个位移增量{Δδ } ，累积后即得到位移{ δ} 增量法是用一系列线性问题去近似非线性问题，实质上是用分段线性的折线去代替非线性曲线一、始点刚度法刚度矩阵[Ki-1]是根据应力应变关系在第i步的开始计算的二、中点刚度法求出 →改进方法或2.6.2 迭代法用迭代法求解非线性问题时，一次施加全部荷载，然后逐步调整位移，使基本方程 得到满足一、直接迭代法 先给出一个近似解 如 ，由应力应变关系求出 第一次近似解为从第n次近似解求出第n+1次近似解的公式收敛准则直接迭代法每步采用的都是割线刚度矩阵二、牛顿（Newton-Raphson）法由此得第n+1次近似解为设 是第n次近似解，一般地，有在 附近将 式作泰勒展开，并只保留线性项，有三、修正牛顿法对大型问题来说形成刚度矩阵并求逆是很费计算时间的。

牛顿法在每次迭代中都要重新建立刚度矩阵并求逆，一次计算时间较长如果只在第1次迭代时计算刚度矩阵 并求出逆阵 在以后的迭代中都用这个逆阵进行计算那么第n步的迭代公式为该法每次迭代节省计算时间较多，虽然迭代过程中的收敛速度有所降低，但在大多数情况下，总的计算时间还是比牛顿法省为提高收敛速度，可以在每经过k次迭代以后重新计算一个 这样，在第1步计算中对 三角分解并存储，在以后个步迭代中只需按上式进行简单回代就行了这种方法称为修正牛顿法 四、 的计算初应力法初应力法材料的应力应变关系为初应力引起的单元结点力为将结点周围有关单元的结点力加以集合，得到初应力引起的结点失衡力为初应变法初应变法在某些问题中，难以用应变明显地表示应力，如徐变问题相反，可以用应力明显地表示应变线弹性应力应变关系为在应力应变关系中引进初应变 ，使得初应变 由下式计算初应变引起的结点失衡力2.6.3 混合法混合法是同时采用增量法和迭代法把荷载划分成较少的几个增量，对每一荷载增量进行迭代计算增量法的优点:适用范围广泛，即其通用性强，收敛性好；另一个优点是它可提供荷载——位移过程线。

除个别情况外，适用于各种类型和各种程度的非线性问题增量法的缺点：一是它比迭代法通常要消耗更多的计算时间；二是不知道近似解与真解相差多少 迭代法的优点是：计算量比增量法小一些，对计算精度也能加以控制，比较适合与加载无关的材料非线性问题和一般的几何非线性问题迭代法不能给出荷载——位移过程线，适用范围也小一些，例如当材料变形特性与加荷过程有关时（加荷与卸荷异性)，以及动力问题等，迭代法均不能使用2.7 材料非线性问题的有限元法2.7.1 材料非线性本构关系一、非线性弹性介质本构关系线弹性应力应变关系非线性弹性应力应变关系全量式应力应变关系 和弹性矩阵 有相同的形式，但它的元素不再是常数，而是应变或应力的函数其中的 ，如果从单向应力应变关系，就分别为割线弹性模量和割线泊松比，因此 也称为割线弹性矩阵增量式应力应变关系称为切线弹性矩阵二、弹塑性介质本构关系 由于弹塑性材料受外部作用的反应和加载路径有关，因此，本构关系应写成增量形式，又因弹塑性状态下加载和卸载有不同的规律，所以其本构关系的表述要比非线性弹性情况复杂 在复杂应力状态下，物体内某一点产生塑性变形时，应力必须满足一定的条件，即屈服条件（屈服准则，如Tresca，Mises，Mohr-Coulomb，Drucker-Prager等等）。

一般来说，它是应力分量的函数式中： 为与材料有关的常数； 为屈服函数弹性状态 塑性状态岩石、土体，混凝土等，其屈服条件受静水压力影响，屈服函数一般为金属等，其屈服条件不受静水压力影响，屈服函数一般为2.7.2 弹塑性增量理论的有限元解法一、增量理论的弹塑性矩阵 根据流动理论，在一个无限小应力增量间隔中，应变增量可看作是由弹性和塑性两部分组成的，即塑性位势理论广义胡克定律根据屈服准则硬化材料 理想塑性体 在Drucker公设成立的条件下，，这时的塑性本构关系为相关连才是对称矩阵在一般情况下，是不对称矩阵的流动法则，只有在这种情况下，二、弹塑性增量理论有限元的数值解法设屈服面为在本增量步开始时，已知应力 ，有效应变 硬化参数 和变位弹塑性本构方程1.计算变位增量2.计算试探应力按弹性关系计算用屈服函数验证3.比例因子设其中r 是一个比例因子因为C在屈服面上，它应该满足屈服条件4.塑性应力增量应变增量分成两部分：弹性应变增量 和塑性应变增量 塑性应力增量弹性应力增量塑性应力增量近似计算如果增量步开始时，应力状态就在屈服面上，只需在上述计算中令 r＝05.应力拉回屈服面第n+1步增量末尾的应力如下：假设应力的修正是沿着屈服面的法线方向进行设 暂时为常量，把屈服函数做一阶泰勒展开，有 应力修正量为一、混凝土的应力应变关系2.7.3 混凝土结构的徐变应力分析单向应力混凝土弹性模量徐变泊松比近似取弹性应变泊松比徐变变形可逆徐变 不可逆徐变 应力松弛复杂应力状态下混凝土的应力应变关系 1。

单向应力作用下应变增量的计算瞬时弹性应变徐变应变二、混凝土结构的徐变应力分析 弹性应变增量中点龄期徐变应变增量不同时刻的徐变变形为徐变应变增量2复杂应力作用下的应变增量弹性应变增量徐变增量应变增量应力增量应力增量和应变增量的关系四、平衡方程组结构刚度矩阵徐变引起的荷载增量温度变形引起的荷载增量一、单向应力作用下的应力应变关系2.7.4 粘弹性问题的有限单元法常泊松比粘弹性体的应力分析 弹性应变增量列阵 徐变应变增量列阵 二、一般粘弹性体的应力分析剪切变形的徐变柔量 体积变形的徐变柔量 应变偏张量 体积变形 弹性应变偏量增量为徐变应变偏量增量为弹性体积变形增量为徐变体积应变增量为二、一般粘弹性体的应力分析虚位移原理 如何用虚位移原理推导单元刚度矩阵形函数的基本要求（单元）刚度矩阵的特点非线性有限元的求解方法弹塑性增量理论有限元求解的主要步骤2.2 变分原理  根据流动理论，在一个无限小应力增量间隔中，应变增量可看作是由弹性和塑性两部分组成的，即而弹性应变增量 与应力增量 间仍为线性关系根据一致性条件和流动法则，可得正则屈服面的增量本构方程为其中切线弹塑性矩阵切线塑性矩阵。