有限元与有限差分法基础

1,第二讲 有限元与有限差分法基础,CAE的工具： 有限元法（FEM）、有限差分法（FDM）、边界元法（BEM）、有限体积法（FVM）、无网格法等等 在材料成形的CAE中主要使用的是有限元法和有限差分法,2,“ 有限元法 ” 的基本思想早在20世纪40年代初期就有人提出，但真正用于工程中则是电子计算机出现以后。 “ 有限元法 ” 这一名称是1960年美国的克拉夫（Clough，R.W.）在一篇题为 “平面应力分析的有限元法” 论文中首先使用。此后，有限元法的应用得到蓬勃发展。 到20世纪80年代初期国际上较大型的结构分析有限元通用程序多达几百种，从而为工程应用提供了方便条件。由于有限元通用程序使用方便，计算精度高，其计算结果已成为各类工业产品设计和性能分析的可靠依据。,3,有限元法最初用于飞机结构的强度设计，由于它在理论上的通用性，因而它可用于解决工程中的许多问题。 目前，它可以解决几乎所有的连续介质和场的问题，包括热传导、电磁场、流体动力学、地质力学、原子工程和生物医学等方面的问题。 机械设计中，从齿轮、轴、轴承等通用零部件到机床、汽车、飞机等复杂结构的应力和变形分析（包括热应力和热变形分析）。 有限元法不仅可以解决工程中的线性问题、非线性问题，而且对于各种不同性质的固体材料，如各向同性和各向异性材料，粘弹性和粘塑性材料以及流体均能求解； 对于工程中最有普遍意义的非稳态问题也能求解。,4,2.1 有限元法基础,基本思想： 将一个连续求解域（对象）离散（剖分）成有限个形状简单的子域（单元） 利用有限个节点将各子域连接起来 在给定的初始条件和边界条件下进行综合计算求解，从而获得对复杂工程问题的近似数值解,5,物理系统举例,几何体 载荷 物理系统,6,有限元模型,真实系统,有限元模型,有限元模型 是真实系统理想化的数学抽象。,定义,7,自由度（DOFs）,自由度(DOFs) 用于描述一个物理场的响应特性。,结构 DOFs,ROTZ,UY,ROTY,UX,ROTX,UZ,8,节点(node)和单元(element) 网格（grid）,节点: 空间中的坐标位置，具有一定自由度和存在相互物理作用。,单元: 一组节点自由度间相互作用的数值、矩阵描述（称为刚度或系数矩阵)。单元有线、面或实体以及二维或三维的单元等种类。,有限元模型由一些简单形状的单元组成，单元之间通过节点连接，并承受一定载荷。,载荷,9,节点和单元,信息是通过单元之间的公共节点传递的。,.,.,.,A,B,.,.,.,.,.,.,.,.,A,B,.,.,.,2 nodes,每个单元的特性是通过一些线性方程式来描述的。 作为一个整体，单元形成了整体结构的数学模型。,10,节点和单元,节点自由度是随连接该节点 单元类型 变化的。,J,I,I,J,J,K,L,I,L,K,I,P,O,M,N,K,J,I,L,三维杆单元 (铰接),UX, UY, UZ,三维梁单元,二维或轴对称实体单元,UX, UY,三维四边形壳单元,UX, UY, UZ,三维实体热单元,TEMP,J,P,O,M,N,K,J,I,L,三维实体结构单元,ROTX, ROTY, ROTZ,ROTX, ROTY, ROTZ,UX, UY, UZ,UX, UY, UZ,11,为什么要离散？,1.无法得到复杂实际问题的解析解 2.将域划分成一些微小而形状规则的单元后，便于在一个单元内得到近似解 3.域中所有单元的解可视为该复杂问题的近似解,12,有限元分析的过程,1.连续体离散化 2.单元分析 3.整体分析 4.确定约束条件 5.方程求解 6.结果分析与讨论,13,1.连续体离散化,连续体：是指所求解的对象（如物体或结构）。 离散化（划分网格或网络化）：是将所求解的对象划分为有限 个具有规则形状的微小块体，把每个微小块体称为单元，相邻两个 单元之间只通过若干点互相连接，每个连接点称为节点。 相邻单元只在节点处连接，载荷也只通过节点在各单元之间传 递，这些有限个单元的集合体，即原来的连续体。 *单元划分后，给每个单元及节点进行编号； *选定坐标系，计算各个节点坐标； *确定各个单元的形态和性态参数以及边界条件等。,14,单元的划分基本上是任意的，一个结构体可以有多种划分结果。但应遵循以下划分原则： (1) 分析清楚所讨论对象的性质，例如，是桁架结构还是结构物，是平面问题还是空间问题等等。 (2) 单元的几何形状取决于结构特点和受力情况，单元的几何尺寸(大小)要按照要求确定。一般来说，单元几何形体各边的长度比不能相差太大。 (3) 有限元模型的网格划分越密，其计算结果越精确，但计算工作量就越大。因此，在保证计算精度的前提下，单元网格数量应尽量少。 (4) 在进行网格疏密布局时，应力集中或变形较大的部位，单元网格应取小一些，网格应划分得密一些，而其他部分则可疏一些。,15,(5) 在设计对象的厚度或者弹性系数有突变的情况下，应该取相应的突变线作为网格的边界线； (6) 相邻单元的边界必须相容，不能从一单元的边或者面的内部产生另一个单元的顶点。 (7) 网格划分后，要将全部单元和节点按顺序编号，不允许有错漏或者重复。 (8) 划分的单元集合成整体后，应精确逼近原设计对象。原设计对象的各个顶点都应该取成单元的顶点。 所有网格的表面顶点都应该在原设计对象的表面上。所有原设计对象的边和面都应被单元的边和面所逼近。,16,有限元分析模型图例,将悬臂梁划分为许多三角形单元 三角形单元的三个顶点都是节点 载荷直接施加在节点上,悬臂梁及其有限元模型,17,2.单元分析,连续体离散化后，即可对单元体进行特性分析，简称为单元分析。 单元分析工作主要有两项： (1)选择单元位移模式(位移函数) 用节点位移来表示单元体内任一点的位移、应变和应力，就需 搞清各单元中的位移分布。 一般是假定单元位移是坐标的某种简单函数，用其模拟内位移的分布规律，这种函数就称为位移模式或位移函数。通常采用的函数形式多为多项式。 根据所选定的位移模式，就可以导出用节点位移来表示单元体内任一点位移的关系式。,18,2.单元分析(2),(2) 分析单元的特性，建立单元刚度矩阵 进行单元力学特性分析，将作用在单元上的所有力（表面 力、体积力、集中力）等效地移置为节点载荷； 采用有关的力学原理建立单元的平衡方程，求得单元内节 点位移与节点力之间的关系矩阵单元刚度矩阵。,19,3. 整体分析,把各个单元的刚度矩阵集成为总体刚度矩阵，以及将各单元的节点力向量集成总的力向量，求得整体平衡方程。 集成过程所依据的原理是节点变形协调条件和平衡条件。,20,4. 确定约束条件,由上述所形成的整体平衡方程是一组线性代数方程，在求解之前，必修根据具体情况分析，确定求解对象问题的边界约束条件，并对这些方程进行适当修正。,21,5. 有限元方程求解,应用有限元法求解机械结构应力类问题时，根据未知量和分析 有三种基本解法：,位移法 力法 混合法,22,(1)位移法 以节点位移作为基本未知量，通过选择适当的位移函数，进行单元的力学特性分析。在节点处建立单元刚度方程，再组合成整体刚度矩阵，求解出节点位移后，进而由节点位移求解出应力。 位移法优点是比较简单，规律性强，易于编写计算机程序。所以得到广泛应用，其缺点是精度稍低。 (2)力法 以节点力作为基本未知量，在节点处建立位移连续方 程，求解出节点力后，再求解节点位移和单元应力。 力法的特点是计算精度高。 (3)混合法 取一部分节点位移和一部分节点力作为基本未知量，建立平衡方程进行求解。,23,单元特性的推导方法,单元刚度矩阵的推导是有限元分析的基本步骤之一。目前，建立单元刚度矩阵的方法主要有以下四种：, 直接刚度法 虚功原理法 能量变分法 加权残数法,24,1. 直接刚度法 直接刚度法是直接应用物理概念来建立单元的有限元方程和分析单元特性的一种方法。这一方法仅能适用于简单形状的单元，如梁单元。但它可以帮助理解有限元法的物理概念。,图1所示是xoy平面中的一简支梁简图，现以它为例，来说明用直接刚度法建立单元刚度矩阵的思想和过程。,图1平面简支梁元及其计算模型,25, 梁在横向外载荷（可以是集中力或分布力或力矩等）作用下产生弯曲变形，在水平载荷作用下产生线位移。 对于该平面简支梁问题： 梁上任一点受有三个力的作用： 水平力Fx， 剪切力Fy ， 弯矩Mz。 相应的位移为： 水平线位移u， 挠度v ， 转角 z 。,由上图可见：,水平线位移和水平力向右为正， 挠度和剪切力向上为正， 转角和弯矩逆时针方向为正。, 通常规定：,26,为使问题简化，可把图示的梁看作是一个梁单元。 如图1所示，当令左支承点为节点 i ，右支承点为节点 j 时， 则该单元的节点位移和节点力可以分别表示为：,称为单元的节点位移列阵。,称为单元的节点力列阵；若 F 为外载荷，则称为载荷列阵。,（1-1）,（1-2）,写成矩阵形式为,q(e)=ui ,vi , zi ,vj ,uj , zjT,ui ,vi , zi ,vj ,uj , zj,F(e)=Fxi ,Fyi ,Mzi ,Fxj ,Fyj ,MzjT,Fxi ,Fyi ,Mzi ,Fxj ,Fyj ,Mzj,27,显然，梁的节点力和节点位移是有联系的。在弹性小变形范围 内，这种关系是线性的，可用下式表示,或,（1-3b）,（1-3a）,28,上式(1-3b)称为单元有限元方程，或称为单元刚度方程，它代表了单元的载荷与位移之间（或力与变形之间）的联系； 式中，K(e)称为单元刚度矩阵，它是单元的特性矩阵。,对于图1所示的平面梁单元问题，利用材料力学中的杆件受力与变形间的关系及叠加原理，可以直接计算出单元刚度矩阵K(e)中的各系数 kst( s, t = i, j ) 的数值,29,2. 虚功原理法,下面以平面问题中的三角形单元为例，说明利用虚功原理法来建立单元刚度矩阵的步骤。 如前所述，将一个连续的弹性体分割为一定形状和数量的单元，从而使连续体转换为有限个单元组成的组合体。单元与单元之间仅通过节点连结，除此之外再无其他连结。也就是说，一个单元上的只能通过节点传递到相邻单元。,从分析对象的组合体中任取一个三角形单元： 设其编号为 e ， 三个节点的编号为i、j、m， 在定义的坐标系 xoy 中，节点坐 标分别为(x j , y j)、(xi , y i)、(xm, ym)，如图2所示。,图2三节点三角形单元,30,由弹性力学平面问题的特点可知，单元每个节点有两个位移分 量，即每个单元有6个自由度，相应有6个节点载荷，写成矩阵形式， 即,单元节点载荷矩阵： F(e)=Fxi ,Fyi ,Fxj ,Fyj ,Fxm ,FymT,单元节点位移矩阵： q(e)=ui ,vi ,uj ,vj ,um ,vmT,图2三节点三角形单元,31,(1)设定位移函数,按照有限元法的基本思想：首先需设定一种函数来近似表达单元内部的实际位移分布，称为位移函数，或位移模式。,三节点三角形单元有6个自由度，可以确定 6个待定系数，故三角形单元的位移函数为,（1-4）,式(1-4)为线性多项式，称为线性位移函数，相应的单元称为线性单元。,u=u(x,y)= 1+ 2x+ 3y v=v(x,y)= 4+ 5x+ 6y,32,上式(5-5)也可用矩阵形式表示，即,式中，d为单元内任意点的位移列阵。,（1-5）,33,由于节点 i、j、m 在单元上，它们的位移自然也就满足位移函数式(1-4)。 设三个节点的位移值分别为( ui, vi)、( uj, vj )、( um, vm )，将节 点位移和节点坐标代入式(1-4)，得,34,（1-6）,式中,（1-7）,由上可知，共有6个方程，可以求出6个待定系数。解方程，求得 各待定系数和节点位移之间的表达式为,为三角形单元的面积。其中：,35,（1-8）,将式(1-7)及式(1-8)、式(1-9)代入式(1-6)中，得到,（1-9）,（1-10）,36,式中，矩阵N 称为单元的形函数矩阵； 为单元节点位移列阵。 其中， 为单元的形函数，它们反映单元内位移的 分布形态，是x, y 坐标的连续函数，且有,（1-11）,式(1-10)又可以写成,（1-12）,上式清楚地表示了单元内任意点位移可由节点位移插值求出。,37,(2) 利用几何方程由位移函数求应变,根据弹性力学的几何方程 ， 线应变 剪切应变 则应变列阵可以写成,式中，B称为单元应变矩阵，它是仅与单元几何尺寸有关的常量矩阵，即,（1-13）,38,（1-14）,上述方程(1-13)称为单元应变方程，它的意义在于： 单元内任意点的应变分量亦可用基本未知量即节点位移分量来表示。,39,(3)利用广义虎克定律求出单元应力方程,根据广义虎克定律，对于平面应力问题,上式（1-15）也可写成,（1-15）,（1-16）,式中， 为应力列阵； D 称为弹性力学平面问题的弹性矩阵，并有,40,则有如下单元应力方程,由式(1-18)可求单元内任意点的应力分量，它也可用基本未知量即 节点位移分量来表示。,（1-17）,（1-18）,41,(4)由虚功原理求单元刚度矩阵,根据虚功原理，当弹性结构受到外载荷作用处于平衡状态时， 在任意给出的微小的虚位移上，外力在虚位移上所做的虚功 AF等 于结构内应力在虚应变上所存储的虚变形势能 A ，即,设处于平衡状态的弹性结构内任一单元发生一个微小的虚位移， 则单元各节点的虚位移 为,（1-20）,（1-21）,（1-19）,则单元内部必定产生相应的虚应变，故单元内任一点的虚应变 为,42,显然，虚应变和虚位移之间关系为,设节点力为,则外力虚功为,（1-24）,（1-22）,（1-23）,单元内的虚变形势能为,43,根据虚功原理,因为,（1-26）,（1-25）,代人式(1-25)，则有,式中， ，均与坐标 x, y 无关，故可以从积分符号中提出，可得：,44,其中，单元刚度矩阵,（1-27）,式(1-27)称为单元有限元方程，或称单元刚度方程，其中 是单元刚度矩阵。,（1-28）,因为三角形单元是常应变单元，其应变矩阵B 、弹性矩阵D 均为常量，而 ，所以式(1-28)可以写成,（1-29）,45,式中，t 为三角形单元的厚度； 为三角形单元的面积。,对于图2所示的三角形单元，将D 及B代入式(1-28)，可以得到单元刚度为,（1-30）,式中: K为66阶矩阵，其中每个子矩阵为22阶矩阵，由下式给出,（1-31）,46,按照力学的一般说法，任何一个实际状态的弹性体的总位能是 这个系统从实际状态运动到某一参考状态(通常取弹性体外载荷为 零时状态为参考状态)时它的所有作用力所做的功。弹性体的总位 能 是一个函数的函数，即泛函，位移是泛函的容许函数。,从能量原理考虑，变形弹性体受外力作用处于平衡状态时，在 很多可能的变形状态中，使总位能最小的就是弹性体的真正变形， 这就是最小位能原理。用变分法求能量泛函的极值方法就是能量变 分原理。能量变分原理除了可解机械结构位移场问题以外，还扩展 到求解热传导、电磁场、流体力学等连续性问题。,3. 能量变分原理法,47,该方法是将假设的场变量的函数(称为试函数)引入问题的控制方程式及边界条件，利用最小二乘法等方法使残差最小，便得到近似的场变量函数形式。 该方法的优点是不需要建立要解决问题的泛函式，所以，即使没有泛函表达式也能解题。,4. 加权残数法,48,有限元解的收敛性 有限元解是近似解 近似解是否收敛于真实解、近似解收敛速度、近似解的稳定性 近似解的收敛条件: 1.完备性准则（必要条件） 试探函数（插值函数）的次数（m）不小于场函数的最高可导阶数 2.协调性准则（充分条件） 试探函数在m-1次连续可导。,49,有限元分析的误差,有限元 分析误差,建模误差,计算误差,离散误差,物理离散误差,几何离散误差,边界条件误差,单元形状误差,舍入误差,截断误差,插值函数与真实函数之间的差异 1.减小单元特征尺寸，称为h法 2.提高插值函数的阶次，称为p法,单元组合体与求解对象几何形状的差异 1.网格局部加密 2.选用边或面上带有节点的单元,边界条件的复杂性 1.准确测定，完善模型 2.细分边界网格,单元严重畸变而退化 细分局部网格或者控制调整关键区域的网格,数据储存,计算方法、解题性质、解题规模,注意网格的划分 选择合适的解算方法 控制解题的规模,减少运算次数，降低解题规模,选择合适的解算方法，控制解题规模,50,材料成形中的非线性问题,1.材料非线性 材料本构方程非线性 弹塑性、刚弹性、刚黏塑性、黏弹塑性 2.几何非线性 3.边界非线性,51,2.2 有限差分法基础,一种直接将微分问题转变成代数问题的近似数值解法。 基本思想 数值微分法 是把连续的定解区域划分成差分网格，用有限个节点代替原连续求解域。 把原方程和定解条件中的微商用差商来近似 把原微分方程和定解条件近似地用代数方程组代替，即有限差分方程,52,差分网格通常为矩形在边界不规则或者形状复杂时精度降低,有限元网格,有限差分网格,53,差分概念,自变量x的解析函数 y =f (x)，则有： dx,dy自变量和函数的微分 函数对自变量的一阶导数 函数对自变量的一阶差商,差商,54,差分方向,向前差分 向后差分 中心差分,55,差商的截断误差,将函数f (x+x)按Taylor级数展开 向前 向后 中心,56,二阶中心差商 通常采用向前差商的向后差商 截断误差与(x)2 同一数量级 一阶向前差商 一阶向后差商 一阶计算精度 一阶中心差商 二阶中心差商 二阶计算精度,57,我们在弹性体上，用相隔等间距h而平行于坐标轴的两组平行线织成正方形网格，x=y=h，如图。,设f=f(x,y)为弹性体内的某一个连续函数。该函数在平行于x轴的一根网线上，如在-上,它只随x坐标的改变而变化。在邻近结点处，函数f可展为泰勒级数如下：,58,我们将只考虑离开结点充分近的那些结点，即（x-x0）充分小。于是可不计（x-x0）的三次及更高次幂的各项，则上式简写为：,在结点，x=x0-h, 在结点1, x=x0+h，代入(b) 得：,59,联立(c),(d),解得差分公式：,同理，在网线-上可得到差分公式,60,从而可导出其它的差分公式如下：,61,相隔2h的两结点处的函数值来表示中间结点处的一阶导数值，可称为中点导数公式。,以相邻三结点处的函数值来表示一个端点处的一阶导数值，可称为端点导数公式。,中点导数公式与端点导数公式相比，精度较高。因为前者反映了结点两边的函数变化，而后者却只反映了结点一边的函数变化。,62,边界元法简介,边界元法（boundary element method） 一种结合有限元法和边界积分法发展起来的一种新数值方法 只在定义域的边界上划分单元，用满足控制方程的函数去逼近边界条件。 适用于应力（薄板）、流体力学、声场、电磁场等问题,63,边界元法基本思想,以微分控制方程的基本解为权函数，利用加权余量法将区域积分转化为边界积分，并结合求解域边界的离散，构建基于边界单元的代数方程组，然后进行计算求解 以定义在边界上的边界积分方程为控制方程，通过对边界元插值离散，化为代数方程组求解 降低了问题的维数，从而显著降低了自由度数 边界的离散比区域的离散方便得多，可用较简单的单元准确地模拟边界形状，最终得到阶数较低的线性代数方程组,64,加权余量法简介,一种直接从所需求解的微分方程及边界条件出发，寻求边值问题近似解的数学方法。 从静力发展到了动力、稳定、材料非线性和几何非线性等各方面。 在求解域上建立一个试函数 试函数由完备函数集的子集构成。已被采用过的试函数有幂级数、三角级数、样条函数、贝赛尔函数、切比雪夫和勒让德多项式等等。 试函数与真实解之间的偏差，即余量（内部和边界） 引入权函数，定义消除余量的条件 加权余量法就是一种定义近似解与真解之间余量，并设法使其最小的方法。,65,设问题的控制微分方程为: 在V域内 在S边界上 式中 : L、B分别为微分方程和边界条件中的微分算子； f、g 为与未知函数u无关的已知函数域值； u为问题待求的未知函数。,66,当利用加权余量法求近似解时，首先在求解域上建立一个试函数 ， 一般具有如下形式：,式中： 待定系数，也可称为广义坐标；,取自完备函数集的线性无关的基函数。,由于 一 般只是待求函数u的近似解，因此代入后将得不到满足，若记：,在V域内,在S边界上,显然 反映了试函数与真实解之间的偏差，它们分别称 做内部和边界余量。,67,若在域V内引入内部权函数 ，在边界S上引入边界权函数 则可建立n个消除余量的条件，一般可表示为：,不同的权函数 和 反映了不同的消除余量的准则。从上 式可以得到求解待定系数矩阵C的代数方程组。一经解得待定 系数，由式(5.1.3)即可得所需求解边值问题的近似解。,68,由于试函数 的不同，余量 和 可有如下三种情况， 依此加权余量法可分为： 1内部法 试函数满足边界条件，也即 此时消除余量的条件成为: 2边界法 试函数满足控制方程，也即 此时消除余量的条件为:,69,3混合法 试函数不满足控制方程和边界条件，此时消除余量的条件为: 显然，混合法对于试函数的选取最方便，但在相同精度条件下，工作量最大。对内部法和边界法必须使基函数事先满足一定条件，这对复杂结构分析往往有一定困难，但试函数一经建立，其工作量较小。,70,边界元法特点,1 前处理工作量小 2 解算规模小 3 求解奇异性问题时计算精度高 4 在载荷集中和半无限域等问题上有优势 相对于有限元法，边界元法发展较慢，,71,有限元法解决应力集中问题,在应力分析中对于应力集中区域必须划分很多的单元，从而增加了求解方程的阶数，计算费用也就随之增加 用位移型有限元法求解出的应力的精度低于位移的精度，对于一个比较复杂的问题必须划分很多单元，相应的数据输人量就很大，同时，在输出的大量信息中，又有许多并不是人们所需要的。,72,1.使问题的维数降低一维 原为三维空间的可降为二维空间，原为二维空间的问题可降为一维。 2.只需将边界离散而需将区域离散化 所划分的单元数目远小于有限元，这样它减少了方程组的方程个数和求解问题所需的数据，不但减少了准备工作，而且节约了计算时间。,边界元法与有限元法比较,73,3.由于直接建立在问题控制微分方程和边界条件上的，不需要事先寻找任何泛函。有限元法是以变分问题为基础，如果泛函不存在就难于使用。 可以求解经典区域法无法求解的无限域类问题。 4.由于边界元法引入基本解，具有解析与离散相结合的特点，因而具有较高的精度。,74,边界元法的缺点,对于非均值和非线性问题求解比较困难 用边界元法解非线性问题时，遇到同非线性项相对应的区域积分，这种积分在奇异点附近有强烈的奇异性，使求解遇到困难。,75,用迭代法可逐步精确方程 根的近似值，但必须要找到 的等价方程 ,如果 选得不合适,不仅影响收敛速度,而且有可能造成迭代格式发散。能否找到一种迭代方法, 结构简单,收敛速度快。这里通过牛顿迭代法的简单介绍，使大家了解数值求解过程。,牛顿迭代法,计算方法,76,计算方法,取x0作为初始近似值，将f(x)在x0做Taylor展开:,重复上述过程 ,作为第一次近似值,一、牛顿迭代法,基本思想：将非线性方程f(x)=0 线性化,Newton 迭代公式,77,计算方法,二、牛顿法的几何意义,x 1,x 2,牛顿法也称为切线法,78,计算方法,例 1 用Newton迭代法求方程xex-1=0在0.5附近的根,精度要求=10-5.,解 Newton迭代格式为,79,计算方法,例2,解：设,取,则由Newton迭代公式,用 Newton 迭代法求,选取初始值,