有限元与有限差分法基础.ppt

1、有限元法基础及有限差分法基础,有限元法,有限差分法,有限元法基础,有限元发展过程 有限元应用 有限元发展方向,有限元法的基本思想,基本思想,1）将连续的求解系统离散为一组由节点相互联在一起的单元组合体,2）在每个单元内假设近似函数来分片表示系统的求解场函数,有限元法的基本思想,有限元法的基本思想,有限元法的基本思想,有限元法的基本思想,有限元法的基本思想,离散为单元网格的冲压件仍然要保证是一个连续体，单元与单元之间没有裂缝、不能重叠，所有单元通过单元节点相互关联着 板料无论产生多大的塑性变形，单元与单元之间依然不会产生裂缝、交叉和重叠，关联单元的节点也不能脱开,有限元法的基本思想,不合格单元,

2、单元裂缝,单元重叠,有限元法的基本思想,变形前后单元之间都是连续的,变形前的网格,变形后的网格,有限元法的基本思想,基本思想 通过在单元内假设不同的插值函数，建立不同的单元模型，适应各种各样的变形模式和受力模式,有限元法的基本思想,有限元法分类,1）位移法：基于最小势能原理或虚功原理,2）力法： 基于最小余能原理,3）杂交法：基于修正余能原理,4）混合法：基于Reissner变分原理,有限元法的基本思想,位移法基本过程,1）离散化过程,3）约束处理过程,2）单元平衡方程组装过程,5）应变、应力回代过程,4）方程组求解过程,离散化过程,最小势能原理,弹性体的势能,为弹性体变形后所具有的内能,为弹

3、性体所受的外力功,离散化过程, 为弹性体的应变, 为弹性体的应力,u为弹性体的可容位移,弹性体处于平衡状态时，其势能应为最小,0,离散化过程,单元插值关系,单元几何关系,单元本构关系,N为单元形函数矩阵,L为单元几何微分算子,为单元弹性矩阵,单元节点自由度向量,离散化过程,B 称为应变矩阵,单元平衡方程或单元刚度方程,k 称为单元刚度矩阵,f 称为单元载荷向量,单元刚度矩阵的特性,对称性,奇异性,主元恒正且对角占优,离散化过程,线弹性问题几何方程三维问题,三维问题,线弹性问题几何方程二维问题,二维问题,平面应力和平面应变状态,线弹性问题几何方程二维问题,二维问题,轴对称状态,线弹性问题几何方程

4、一维问题,一维问题,线弹性问题本构方程三维问题,三维问题,E为弹性模量；为泊松比,线弹性问题本构方程平面应力,二维问题,平面应力状态,线弹性问题本构方程平面应力,平面应力状态,线弹性问题本构方程平面应变,二维问题,平面应变状态,线弹性问题本构方程平面应变,平面应变状态,线弹性问题本构方程轴对称,二维问题,轴对称状态,线弹性问题本构方程轴对称,二维问题,轴对称状态,线弹性问题本构方程轴对称,轴对称状态,线弹性问题本构方程一维问题,一维问题,常用单元模型,单元模型,插值关系,一一对应,单元类型,一维单元、二维单元、三维单元,等参单元、超参单元、次参单元,常用单元模型,一维单元,2节点线单元,3节点

5、线单元,梁单元,常用单元模型,二维单元,3节点三角形线性单元,6节点三角形二次单元,常用单元模型,二维单元,10节点三角形三次单元,4节点四边形双线性单元,常用单元模型,二维单元,8节点四边形二次单元,12节点四边形三次单元,常用单元模型,三维单元,4节点四面体线性单元,10节点四面体二次单元,常用单元模型,三维单元,8节点六面体线性单元,20节点六面体二次单元,常用单元模型,准三维空间单元,桁架单元,一维2节点线单元+单元局部随体坐标系,为什么要建立单元局部随体坐标系 ？,简化分析问题的复杂程度。 在局部坐标系中，空间桁架的每根杆每变成了一维2节点线单元,常用单元模型,准三维空间单元,框架单

6、元,三维梁单元+一维2节点线单元+单元局部随体坐标系,两端都是刚性联结,可以要承受拉压、弯曲、扭转3种变形模式,框架单元的特点,常用单元模型,准三维空间单元,板单元,薄板单元,中厚板单元,弯曲和横向剪切2种变形模式抵抗板的变形,如果板很薄，忽略横向剪切抗力，认为抵抗载荷的主要因素是弯矩,常用单元模型,准三维空间单元,壳单元,抵抗拉压变形的二维单元+板单元+单元局部随体坐标系。适合于薄壳单元和中厚壳单元,从几何上分为薄壳单元和中厚壳单元, 组合单元,常用单元模型,准三维空间单元, 壳理论单元,由空间壳理论严格构造的壳单元。适合于薄壳单元和中厚壳单元, 退化单元,由三维实体单元退化成的壳单元。只适

7、合于中厚壳单元,单元模型构造,有限元法的基本思想,通过单元分片近似，在每个单元内假设近似函数来分片表示系统的场函数,选择近似函数,简单、实用的原则 在有限元法中，近似函数称为插值函数,单元模型构造,插值函数,一般都采用多项式函数，主要原因是:,采用多项式插值函数比较容易推导单元平衡方程，特别是易于进行微分和积分运算。,随着多项式函数阶次的增加，可以提高有限元法的计算精度。从理论上说，无限提高多项式的阶数，可以求得系统的精确解。,单元模型构造方法,整体坐标系法 局部坐标系法,Lagrange插值方法 Hermite插值方法,单元模型构造方法,2节点线单元,1,2,u1,u2,x1,x2,u,x,

8、1. 假设插值多项式,2. 利用节点值求 a0 和 a1,单元模型构造方法,3. 代入a0 和 a1，得插值多项式 u(x),4. 按u1 和 u2合并同类项，设 l = x2- x1,单元模型构造方法,关键,如何构造插值多项式 u ?,二维问题三维问题，如何构造插值多项式?,收敛性条件 在单元内，场函数必须是连续的； 完备性：插值多项式的阶次必须由低到高依次增加，不能出现跳跃现象； 协调性：各单元边界必须连续，单元边界不能出现开裂现象。,插值多项式收敛性条件,收敛：当单元逐渐缩小时，如果插值多项式满足收敛性条件，则数值解将收敛于精确解,插值多项式收敛性条件,协调单元 满足插值多项式收敛性条件

9、和的单元 完备单元 满足插值多项式收敛性条件的单元 cr 阶连续性 插值多项式的第r阶导数是连续的,插值多项式收敛性条件,非协调单元与部分协调单元 对于一般固体力学问题来说，协调性要求单元在变形时，相邻单元之间不应引起开裂、重叠或其它不连续现象。例如，梁、板、壳等单元，在单元边界不但要求位移是连续的，而且其一阶导数也必须是连续的。板、壳单元位移函数沿单元边界的法向导数（转角）的连续性一般比较难实现，因此出现了许多不完全满足协调性要求的“非协调单元”或“部分协调单元”，有时它们的精度也很好。,插值多项式选择条件,插值多项式应该尽可能满足其收敛性条件（收敛性） 由插值多项式所确定的场函数变化应该与

10、局部坐标系的选择无关（各向同性） 假设的插值多项式系数的数量应该等于单元的节点数（解的唯一性）,选择条件,插值多项式选择条件,深入分析,由收敛性条件可知，插值多项式中必须含有常数项（刚体位移项），高阶项的次数必须依次增加，不允许有跳跃,插值多项式选择条件,由选择条件可知，插值多项式函数在所有自由度方向上要满足各向同性性，这样就不会随局部坐标系变化而改变了,深入分析,插值多项式选择条件,深入分析,选择条件是为了能由单元节点值唯一确定插值多项式,4节点四边形的插值多项式应该是,插值多项式系数i (i = 0,1,2,3) 也是4个,单元模型构造整体坐标系法,基本思想,针对弹性体有限元网格建立一个统

11、一的坐标系，每个单元的插值多项式都在这个坐标系上建立,单元模型构造整体坐标系法,2节点线单元,1,2,u1,u2,x1,x2,u,x,1. 假设插值多项式,2. 利用节点值求 a0 和 a1,单元模型构造整体坐标系法,3. 代入a0 和 a1，得插值多项式 u(x),4. 按u1 和 u2合并同类项，设 l = x2- x1,单元模型构造整体坐标系法,N1 和 N2 称为单元的形函数； N 称为单元的形函数矩阵； ue 称为单元节点位移向量。,2节点线的单元形函数,单元模型构造整体坐标系法,二维3节点三角形单元,建立整体坐标系oxy,单元模型构造整体坐标系法,1. 假设插值多项式,2. 首先，

12、利用节点值求 0 、 1 和 2,二维3节点三角形单元,单元模型构造整体坐标系法,A为单元面积,单元模型构造整体坐标系法,3. 将 0 、 1 和 2 代入插值多项式，按u1、u2、u3合并同类项,单元模型构造整体坐标系法,4. 同理可得,单元模型构造整体坐标系法,5. 单元插值多项式为,单元模型构造整体坐标系法,6. 单元插值多项式写成矩阵形式（常用）,单元模型构造整体坐标系法,7. 单元插值多项式的另一种矩阵形式（不常用）,单元模型构造整体坐标系法,4节点四面体单元,单元模型构造整体坐标系法,1. 假设插值多项式,2. 插值多项式为,单元模型构造整体坐标系法,（i=1,2,3,4）,循环轮

13、换脚标1、2、3、4，相应可以得到a2，b2 ， c2 ， d2 、 a3 ， b3 ， c3 ， d3 、 a4 ， b4 ， c4 ， d4,单元模型构造整体坐标系法,3. 单元插值多项式写成矩阵形式（常用）,单元模型构造整体坐标系法,4. 单元插值多项式另一种矩阵形式（不常用）,单元模型构造整体坐标系法,从理论上讲，整体坐标系法可以求任意单元的形函数，但计算过程太复杂 只能求一维2节点线单元、二维3节点三角形单元和三维4节点四面体单元3种简单单元的形函数 复杂的或二次以上的单元必须采用局部坐标系法求 位移场 u 是形函数 Ni 的线性组合，因此形函数Ni同样具有插值多项式的特性,单元刚度

14、矩阵2节点线单元,一维2节点线单元,单元插值关系,单元几何关系,单元本构关系,N=N1 N2,De=E,单元刚度矩阵2节点线单元,单元刚度矩阵,A为单元截面积；l为单元长度,矩阵B,单元刚度矩阵三角形单元,二维3角形单元,单元插值关系,单元刚度矩阵三角形单元,单元几何关系,单元刚度矩阵三角形单元,单元本构关系,平面应力问题,单元刚度矩阵三角形单元,矩阵B,单元刚度矩阵三角形单元,单元刚度矩阵,h为单元厚度,k为对称的6*6常数矩阵,A为单元面积,作业,求4节点四面体单元的单元刚度矩阵,单元模型构造整体坐标系法,单元形函数的特性,正规性：单元形函数之和等于1。,正交性：形函数在本节点的值等于1，

15、在其它节点的值等于0。,例如: 2节点线单元形函数,单元模型等参单元,等参单元,单元内任意一点的位移u与单元节点位移ue之间的关系为,一般单元坐标的插值关系也采用与位移插值关系相同的变换关系即单元内任意一点的坐标x与单元节点坐标xe之间的关系为,单元模型等参单元,等参单元,凡是几何形状和位移场采用同阶同参数插值关系来描述的单元，称为等参单元,前面介绍的所有单元都属于等参单元,在描述单元的几何形状和位移场时，并不一定非采用同阶插值关系,单元模型等参单元,等参单元,3节点三角形等参单元,单元模型等参单元,超参单元,如果几何形状插值函数的阶数高于位移场插值函数的阶数，称为超参单元,次参单元,如果几何

16、形状插值函数的阶数低于位移场插值函数的阶数，称为次参单元,单元平衡方程组装过程,为什么要组装 ?,消除内力,组装的原则是什么 ?,单元自由度与结构自由度对应,单元平衡方程组装过程,2,F,1,3,U3,U4,U2,U1,U5,U6,结构自由度向量U,单元平衡方程组装过程,2,单元平衡方程组装过程,2,单元平衡方程组装过程,单元平衡方程组装过程,单元平衡方程组装过程,总体刚度方程,K 称为总体刚度矩阵,U 称为位移向量,F 称为载荷向量,总体刚度矩阵K的特性,对称性,奇异性,稀疏性,非零元素带状分布,约束处理过程,为什么要约束处理 ?,总体平衡方程组是奇异的,消除无限制的刚体运动,使总体平衡方程

17、组存在唯一一组解,约束处理过程边界条件,边界条件分类,力(载荷)边界条件,位移边界条件,集中载荷力,表面分布力,自重力,热交换引起的温度载荷,固定位移约束,强制位移约束,关联位移约束,约束处理过程模型简化,约束处理过程模型简化,约束处理过程约束方程,1,2,3,4,5,6,7,8,9,10,11,12,约束处理过程约束处理方法,位移约束处理方法,赋0赋1法,乘大数法,约束处理过程赋0赋1法,强制位移约束条件处理 U4=C,约束处理过程赋0赋1法,强制位移约束条件处理 U4=C,约束处理过程赋0赋1法,有6个方程，5个未知数，如果约束方程可以消除有限元平衡方程组的奇异性，则取任意5个方程联立求解

18、，都会得到方程组的唯一一组解。,系数矩阵由原来的对称的变成了非对称的，这对于大规模有限元方程组求解是十分不利的，采用相同的求解方法，在求解时间和矩阵存贮容量方面都增加了一倍。,约束处理过程赋0赋1法,为了保证系数矩阵的对称性，去掉方程组第4行,约束处理过程赋0赋1法,引入强制位移约束方程 U4=C，使方程组求解时直接将自由度U4求出,约束处理过程赋0赋1法,固定位移约束条件处理 U4=0,约束处理过程赋0赋1法,基本原理 利用初等变换对求解方程组进行相同的行列变换，既保证方程组解不会改变，又可以保持方程组系数矩阵的对称性。 在进行初等变换时，只要保证对方程组系数矩阵做相同的行列变换，就可以保持

19、方程组系数矩阵的对称性。,约束处理过程乘大数法,乘大数法,基本原理 利用矩阵的初等变换不改变方程组解的思想。,约束处理过程乘大数法,强制位移边界条件,约束处理过程乘大数法,强制约束方程,A是一个大数，是系数矩阵中对角线元素K44的1010倍量级以上,为什么要乘以大数A ？,放大位移约束方程的优势,约束处理过程乘大数法,强制位移边界条件,约束处理过程乘大数法,固定位移边界条件,C = 0,约束后的方程组简化为,约束处理过程两种方法比较,赋0赋1法在约束处理过程中是严格精确的，而乘大数法是一种近似约束处理方法，它的精度取决于所乘大数A值,两种方法都可以消除有限元平衡方程的奇异性，得到符合实际边界条

20、件的唯一一组解。但两种方法还是有很大的区别,约束处理过程两种方法比较,采用乘大数法约束处理后的有限元平衡方程在求解时可能造成解的失真，大数A值越大可能解的偏差会越大，而赋0赋1法就不会出现类似的问题，它在约束过程和求解过程都是精确的 乘大数法相对于赋0赋1法在约束处理过程上简单一些,约束处理过程两种方法比较,赋0赋1法实际上是将关联位移约束方程代入到有限元平衡方程中的，是代入法。而乘大数是将占绝对优势的关联位移约束方程合并到有限元平衡方程中的，是罚方法，计算误差来自于合并过程，计算精度取决于关联位移约束方程的优势大小 商业软件中，位移边界条件的约束处理都采用赋0赋1法，乘大数很少被采用主要原因

21、是它是一种近似方法，而且大数的大小也不好确定，有时还会造成求解失败,约束处理过程弹簧单元,假设柔性弹簧 k,f = kU4,弹簧约束方程 f = kU4,约束处理过程弹簧单元,方程组求解过程特点,方程组求解有限元计算过程中很重要的一部分，在有限元法的发展过程中，有限元方程的求解效率一直是其应用的最大瓶颈之一 有限元方程组的特点： 有限元方程组的系数矩阵具有对称、稀疏、带状分布以及正定、主元占优。有效地利用这些特点，以减少系数矩阵的存贮量，提高方程组求解效率,方程组求解过程分类比较,线性方程组的解法主要分两大类: 直接解法：以高斯消去法基础，以等带宽或变带宽方式存贮系数矩阵内元素，对于求解规模比

22、较大的问题，要存贮的元素非常巨大。 迭代解法：只需要存贮系数矩阵中非零元素，存贮量很小，一般是变带宽存贮量的20%或更少，有些算法的求解效率也非常高，适合求解大规模线性方程组。但是这种解法对接近病态的方程组很难保证收敛性。,方程组求解过程带宽定义,有限元方程组系数矩阵是稀疏的、非零元素呈带状分布，带宽就是它的宽度，带宽的大小是由系统有限元网格的节点号排序决定的，具体求法是,带宽=(单元最大节点号之差+1)*节点自由度数,带宽是网格节点标注方法直接决定的，不同标注方法带宽可能相关很大,方程组求解过程带宽,带宽是网格节点标注方法直接决定的，不同标注方法带宽可能相关很大,方程组求解过程带宽,所示四边

23、形网格的三种节点号标注方法，每个节点是2个自由度 结构的带宽分别是12，18，56，相差很大，其中12和56之间相差近5倍，这就意味着系数矩阵的存贮量也是相差5倍，因此，对于大规模复杂系统的节点号优化是十分必要的,方程组求解过程系数矩阵存贮,系数矩阵存贮 如果节点号排序优化的比较好，系数矩阵的存贮量就会减少很多。根据系数矩阵的对称性，一般都是按半带宽存贮。 系数矩阵存贮的方法 二维等带宽存贮 一维变带宽存贮,方程组求解过程二维等带宽存贮,二维等带宽存贮,方程组求解过程二维等带宽存贮,二维等带宽存贮消除了最大带宽以外的全部零元素，节省了系数矩阵元素的存贮量。但是由于取最大带宽为存贮范围，因此不能

24、排除在带宽内的大量零元素。当系数矩阵的各行带宽变化不大时，适合采用二维等带宽存贮，方程组求解过程中系数矩阵元素的寻址也比较方便，求解效率较高。 当出现局部带宽特别大的情况时，采用二维等带宽存贮时，将由于局部带宽过大而使整体系数矩阵的存贮大大增加。,方程组求解过程一维变带宽存贮,一维变带宽存贮 一维变带宽存贮方法就是把变化的带宽内的元素按一定的顺序存贮在一个一维数组中。由于它不按最大带宽存贮，因此比二维等带宽存贮更节省内存。按照解法可分为按行一维变带宽存贮和按列一维变带宽存贮。,按行一维变带宽存贮,方程组求解过程一维变带宽存贮,辅助的寻址数组M,一维变带宽存贮是最节省内存的一种方法，但是由于要借

25、助于寻址数组寻找系数矩阵元素的位置，相对二维等带宽存贮方法来说要复杂一些，而且在程序实现时也要复杂得多，方程组求解过程中也要消耗一些数组寻址时间。因此，在选用存贮方法时要权衡二者的利弊，统盘考虑。一般当带宽变化不大，计算机内存允许时，采用二维等带宽存贮方法是比较合适的。,方程组求解过程一维变带宽存贮,方程组求解过程求解方法,方程组求解方法 高斯消去法 三角分解法 雅可比（Jacobi）迭代法 高斯-赛德尔（Gauss-Seidel）迭代法,应变、应力回代过程,单元应变和应力回代求解 通过求解有限元平衡方程得到有限元节点位移后，就可以进行系统的刚度校核。如果所分析问题要进行强度校核，就要回代求解

26、单元的应变和应力。 由插值关系和几何关系可得单元应变，再通过本构关系得到单元应力,有限差分法,从弹性力学的基本方程建立以来，这些方程在各种问题的边界条件下如何求解，一直是很多数学工作者和力学工作者研究的内容。即弹性力学的经典解法存在一定的局限性，当弹性体的边界条件和受载情况复杂一点，往往无法求得偏微分方程的边值问题的解析解，许多工程重要问题，不能够得出函数式的解答。 因此，弹性力学问题的各种数值解法便具有重要的实际意义。,工程中常用的数值解法有有限单元法和差分法。,有限单元法 是以有限个单元的集合体来代替连续体，属于物理上的近似。,差分法 是把弹性力学的基本方程和边界条件（一般均为微分方程）近

27、似地改用差分方程（代数方程）来表示，把求解微分方程的问题改换成为求解代数方程的问题，属于数学上的近似。,第一节 差分方程,第二节 应力函数的差分解,第三节 深梁应力函数的差分解,第一节 差分方程,差分法是沿用已久的一种数值解法。随着计算机的普及和相应的软件发展，此法成为解弹性力学问题的一种有效的方法。,我们在弹性体上，用相隔等间距h而平行于坐标轴的两组平行线织成正方形网格，x=y=h，如图。,设f=f(x,y)为弹性体内的某一个连续函数。该函数在平行于x轴的一根网线上，如在-上,它只随x坐标的改变而变化。在邻近结点处，函数f可展为泰勒级数如下：,我们将只考虑离开结点充分近的那些结点，即（x-x

28、0）充分小。于是可不计（x-x0）的三次及更高次幂的各项，则上式简写为：,在结点，x=x0-h, 在结点1, x=x0+h，代入(b) 得：,联立(c),(d),解得差分公式：,同理，在网线-上 可得到差分公式,差分公式(-)及(-)是以相隔2h的两结点处的函数值来表示中间结点处的一阶导数值，可称为中点导数公式。,以相邻三结点处的函数值来表示一个端点处的一阶导数值，可称为端点导数公式。,应当指出：中点导数公式与端点导数公式相比，精度较高。因为前者反映了结点两边的函数变化，而后者却只反映了结点一边的函数变化。因此，我们总是尽可能应用前者，而只有在无法应用前者时才不得不应用后者。,以上(-)(-)

29、是基本差分公式,从而可导出其它的差分公式如下：,第二节 应力函数的差分解,当不计体力时，我们已把弹性力学平面问题归结为在给定边界条件下求解双调和方程的问题。用差分法解平面问题，就应先将双调和方程变换为差分方程，而后求解之。,一旦求得弹性体全部节点的值后，就可按应力分量差分公式（对节点0）算得弹性体各节点的应力。,可见，用差分法解平面问题，共有两大任务： 一、建立差分方程 将（1-68）代入双调和方程,对于弹性体边界以内的每一结点，都可以建立这样一个差分方程。,整理即得,二、联立求解这些线性代数方程，就能求得各内结点处的值。,为了求得边界上各结点处的值，须要应用应力边界条件，即：,一般建立和求解

30、差分方程，在数学上不会遇到很大困难。但是，当对于边界内一行的（距边界为h的）结点，建立的差分方程还将涉及边界上各结点处的值，并包含边界外一行的虚结点处的值。,代入上式，即得：,l1=cos(N,x)=cos=dy/ds, l2=cos(N,y)=sin=-dx/ds,于是，式（a）可改写为：,由右图可见，,关于边界上任一点处,由此得：,的值，可将（b）式从,A点到B点对s积分得到：,将此式亦从A点到B点沿s进行积分，就得到边界上任一点B处的值。为此利用分部积分法，得：,由高等数学可知，,将式(b),(c)代入，整理得：,由前知，把应力函数加上一个线性函数，并不影响应力。因此，可设想把应力函数加上a+bx+cy，然后调整a，b，c三个数值，使得,由式(d)及式(c)可见，设,即可根据面力分量及导数求得,为已知，,从图易看出，式（23）右边的积分式表示A与B之间的x方向的面力之和；式（24）右边的积分式表示A与B之间的y方向的面力之和；式（25）右边的积分式表示A与B之间的面力对于B点的矩。,于