有限元与有限差分法基础ppt课件

1、第二讲第二讲 有限元与有限差分法根底有限元与有限差分法根底 CAE的工具：的工具： 有限元法有限元法FEM、有限差分法、有限差分法FDM、边境元法、边境元法BEM、有限体积法、有限体积法FVM、无网格法等等、无网格法等等 在资料成形的在资料成形的CAE中主要运用的是有限元法和有限差分法中主要运用的是有限元法和有限差分法“ 有限元法有限元法 的根本思想早在的根本思想早在20世纪世纪40年代初期就年代初期就有人提出，但真正用于工程中那么是电子计算机出现以后。有人提出，但真正用于工程中那么是电子计算机出现以后。 “ 有限元法有限元法 这一称号是这一称号是1960年美国的克拉夫年美国的克拉夫Cloug

2、h，R.W.在一篇题为在一篇题为 “平面应力分析的有限元法平面应力分析的有限元法 论文中首先运用。以后，有限元法的运用得到蓬勃开论文中首先运用。以后，有限元法的运用得到蓬勃开展。展。 到到20世纪世纪80年代初期国际上较大型的构造分析有限元年代初期国际上较大型的构造分析有限元通用程序多达几百种，从而为工程运用提供了方便条件。通用程序多达几百种，从而为工程运用提供了方便条件。由于有限元通用程序运用方便，计算精度高，其计算结果由于有限元通用程序运用方便，计算精度高，其计算结果已成为各类工业产品设计和性能分析的可靠根据。已成为各类工业产品设计和性能分析的可靠根据。 有限元法最初用于飞机构造的强度设计

3、，由于它在有限元法最初用于飞机构造的强度设计，由于它在实际上的通用性，因此它可用于处理工程中的许多问题。实际上的通用性，因此它可用于处理工程中的许多问题。目前，它可以处理几乎一切的延续介质和场的问题，目前，它可以处理几乎一切的延续介质和场的问题，包括热传导、电磁场、流体动力学、地质力学、原子工程包括热传导、电磁场、流体动力学、地质力学、原子工程和生物医学等方面的问题。和生物医学等方面的问题。 机械设计中，从齿轮、轴、轴承等通用零部件到机机械设计中，从齿轮、轴、轴承等通用零部件到机床、汽车、飞机等复杂构造的应力和变形分析包括热应床、汽车、飞机等复杂构造的应力和变形分析包括热应力和热变形分析。力和

4、热变形分析。有限元法不仅可以处理工程中的线性问题、非线性问有限元法不仅可以处理工程中的线性问题、非线性问题，而且对于各种不同性质的固体资料，如各向同性和各题，而且对于各种不同性质的固体资料，如各向同性和各向异性资料，粘弹性和粘塑性资料以及流体均能求解；向异性资料，粘弹性和粘塑性资料以及流体均能求解；对于工程中最有普遍意义的非稳态问题也能求解。对于工程中最有普遍意义的非稳态问题也能求解。2.1 有限元法根底有限元法根底根本思想：根本思想：将一个延续求解域对象离散剖分成有限个外形简将一个延续求解域对象离散剖分成有限个外形简单的子域单元单的子域单元利用有限个节点将各子域衔接起来利用有限个节点将各子域

5、衔接起来在给定的初始条件和边境条件下进展综合计算求解，从而在给定的初始条件和边境条件下进展综合计算求解，从而获得对复杂工程问题的近似数值解获得对复杂工程问题的近似数值解 物理系统举例 几何体几何体 载荷载荷 物理系统物理系统构造构造热热电磁电磁有限元模型真实系统真实系统有限元模型有限元模型 有限元模型有限元模型 是真实系统理想化的数学笼统。是真实系统理想化的数学笼统。定义定义自在度DOFs自在度自在度(DOFs) 用于描画一个物理场的呼应特性。用于描画一个物理场的呼应特性。构造构造 DOFs 构造构造 位移位移 热热 温度温度 电电 电位电位 流体流体 压力压力 磁磁 磁位磁位 方向方向 自在

6、度自在度ROTZUYROTYUXROTXUZ节点(node)和单元(element) 网格grid节点节点: 空间中的坐标位置，具有一定自在度空间中的坐标位置，具有一定自在度和存在相互物理作用。和存在相互物理作用。单元单元: 一组节点自在度间相互作用的数值、一组节点自在度间相互作用的数值、矩阵描画称为刚度或系数矩阵矩阵描画称为刚度或系数矩阵)。单元有线。单元有线、面或实体以及二维或三维的单元等种类。、面或实体以及二维或三维的单元等种类。有限元模型由一些简单外形的单元组成，单元之间经有限元模型由一些简单外形的单元组成，单元之间经过节点衔接，并接受一定载荷。过节点衔接，并接受一定载荷。载荷载荷载荷

7、载荷节点和单元节点和单元信息是经过单元之间的公共节点传送的。信息是经过单元之间的公共节点传送的。分别但节点重叠的单元分别但节点重叠的单元A和和B之间没有信息传送之间没有信息传送需进展节点合并处置需进展节点合并处置具有公共节点的单元具有公共节点的单元之间存在信息传送之间存在信息传送 .AB.AB.1 node2 nodes每个单元的特性是经过一些线性方程式来描画的。每个单元的特性是经过一些线性方程式来描画的。作为一个整体，单元构成了整体构造的数学模型。作为一个整体，单元构成了整体构造的数学模型。节点和单元节点和单元节点自在度是随衔接该节点节点自在度是随衔接该节点 单元类型单元类型 变化的。变化的

8、。JIIJJKLILKIPOMNKJIL三维杆单元三维杆单元 (铰接铰接)UX, UY, UZ三维梁单元三维梁单元二维或轴对称实体单元二维或轴对称实体单元UX, UY三维四边形壳单元三维四边形壳单元UX, UY, UZ,三维实体热单元三维实体热单元TEMPJPOMNKJIL三维实体构造单元三维实体构造单元ROTX, ROTY, ROTZROTX, ROTY, ROTZUX, UY, UZ,UX, UY, UZ为什么要离散？为什么要离散？ 1.无法得到复杂实践问题的解析解无法得到复杂实践问题的解析解 2.将域划分成一些微小而外形规那么的单元后，便于在一将域划分成一些微小而外形规那么的单元后，便于

9、在一个单元内得到近似解个单元内得到近似解 3.域中一切单元的解可视为该复杂问题的近似解域中一切单元的解可视为该复杂问题的近似解有限元分析的过程有限元分析的过程 1.延续体离散化延续体离散化 2.单元分析单元分析 3.整体分析整体分析 4.确定约束条件确定约束条件 5.方程求解方程求解 6.结果分析与讨论结果分析与讨论1.延续体离散化延续体离散化 延续体：是指所求解的对象如物体或构造。延续体：是指所求解的对象如物体或构造。离散化划分网格或网络化：是将所求解的对象划分离散化划分网格或网络化：是将所求解的对象划分为有限为有限个具有规那么外形的微小块体，把每个微小块体称为单元，个具有规那么外形的微小块

10、体，把每个微小块体称为单元，相邻两个相邻两个单元之间只经过假设干点相互衔接，每个衔接点称为节点。单元之间只经过假设干点相互衔接，每个衔接点称为节点。相邻单元只在节点处衔接，载荷也只经过节点在各单元相邻单元只在节点处衔接，载荷也只经过节点在各单元之间传之间传递，这些有限个单元的集合体，即原来的延续体。递，这些有限个单元的集合体，即原来的延续体。 *单元划分后，给每个单元及节点进展编号；单元划分后，给每个单元及节点进展编号； *选定坐标系，计算各个节点坐标；选定坐标系，计算各个节点坐标； *确定各个单元的形状和性态参数以及边境条件等。确定各个单元的形状和性态参数以及边境条件等。单元的划分根本上是恣

11、意的，一个构造体可以有多种单元的划分根本上是恣意的，一个构造体可以有多种划分结果。但应遵照以下划分原那么：划分结果。但应遵照以下划分原那么：(1) 分析清楚所讨论对象的性质，例如，是桁架构造分析清楚所讨论对象的性质，例如，是桁架构造还是构造物，是平面问题还是空间问题等等。还是构造物，是平面问题还是空间问题等等。(2) 单元的几何外形取决于构造特点和受力情况，单单元的几何外形取决于构造特点和受力情况，单元的几何尺寸元的几何尺寸(大小大小)要按照要求确定。普通来说，单元几要按照要求确定。普通来说，单元几何形体各边的长度比不能相差太大。何形体各边的长度比不能相差太大。 (3) 有限元模型的网格划分越

12、密，其计算结果越准确，有限元模型的网格划分越密，其计算结果越准确，但计算任务量就越大。因此，在保证计算精度的前提下，但计算任务量就越大。因此，在保证计算精度的前提下，单元网格数量应尽量少。单元网格数量应尽量少。(4) 在进展网格疏密规划时，应力集中或变形较大的在进展网格疏密规划时，应力集中或变形较大的部位，单元网格应取小一些，网格应划分得密一些，而其部位，单元网格应取小一些，网格应划分得密一些，而其他部分那么可疏一些。他部分那么可疏一些。(5) 在设计对象的厚度或者弹性系数有突变的情况下，在设计对象的厚度或者弹性系数有突变的情况下，应该取相应的突变线作为网格的边境限；应该取相应的突变线作为网格

13、的边境限；(6) 相邻单元的边境必需相容，不能从一单元的边或者相邻单元的边境必需相容，不能从一单元的边或者面的内部产生另一个单元的顶点。面的内部产生另一个单元的顶点。(7) 网格划分后，要将全部单元和节点按顺序编号，不网格划分后，要将全部单元和节点按顺序编号，不允许有错漏或者反复。允许有错漏或者反复。(8) 划分的单元集合成整体后，应准确逼近原设计对象。划分的单元集合成整体后，应准确逼近原设计对象。原设计对象的各个顶点都应该取成单元的顶点。原设计对象的各个顶点都应该取成单元的顶点。 一切网格的外表顶点都应该在原设计对象的外表上。一一切网格的外表顶点都应该在原设计对象的外表上。一切原设计对象的边

14、和面都应被单元的边和面所逼近。切原设计对象的边和面都应被单元的边和面所逼近。 有限元分析模型图例有限元分析模型图例 将悬臂梁划分为许多三角形单元将悬臂梁划分为许多三角形单元 三角形单元的三个顶点都是节点三角形单元的三个顶点都是节点 载荷直接施加在节点上载荷直接施加在节点上悬臂梁及其有限元模型悬臂梁及其有限元模型 2.单元分析单元分析 延续体离散化后，即可对单元体进展特性分析，简称为延续体离散化后，即可对单元体进展特性分析，简称为单元分析。单元分析。单元分析任务主要有两项：单元分析任务主要有两项：(1)选择单元位移方式选择单元位移方式(位移函数位移函数) 用节点位移来表示单元体内任一点的位移、应

15、变和用节点位移来表示单元体内任一点的位移、应变和应力，就需应力，就需搞清各单元中的位移分布。搞清各单元中的位移分布。 普通是假定单元位移是坐标的某种简单函数，用其模普通是假定单元位移是坐标的某种简单函数，用其模拟内位移的分布规律，这种函数就称为位移方式或位移函拟内位移的分布规律，这种函数就称为位移方式或位移函数。通常采用的函数方式多为多项式。数。通常采用的函数方式多为多项式。 根据所选定的位移方式，就可以导出用节点位移来根据所选定的位移方式，就可以导出用节点位移来表示单元体内任一点位移的关系式。表示单元体内任一点位移的关系式。2.单元分析单元分析(2) (2) 分析单元的特性，建立单元刚度矩阵

16、分析单元的特性，建立单元刚度矩阵 进展单元力学特性分析，将作用在单元上的一切力进展单元力学特性分析，将作用在单元上的一切力外表外表 力、体积力、集中力等效地移置为节点载荷；力、体积力、集中力等效地移置为节点载荷； 采用有关的力学原理建立单元的平衡方程，求得单元采用有关的力学原理建立单元的平衡方程，求得单元内节内节 点位移与节点力之间的关系矩阵点位移与节点力之间的关系矩阵单元刚度矩阵。单元刚度矩阵。 3. 整体分析整体分析 把各个单元的刚度矩阵集成为总体刚度矩阵，以及将把各个单元的刚度矩阵集成为总体刚度矩阵，以及将各单元的节点力向量集成总的力向量，求得整体平衡各单元的节点力向量集成总的力向量，求

17、得整体平衡方程。方程。 集成过程所根据的原理是节点变形协调条件和平衡条集成过程所根据的原理是节点变形协调条件和平衡条件。件。 4. 确定约束条件确定约束条件 由上述所构成的整体平衡方程是一组线性代数方程，由上述所构成的整体平衡方程是一组线性代数方程，在求解之前，必修根据详细情况分析，确定求解对象在求解之前，必修根据详细情况分析，确定求解对象问题的边境约束条件，并对这些方程进展适当修正。问题的边境约束条件，并对这些方程进展适当修正。5. 有限元方程求解有限元方程求解 运用有限元法求解机械构造应力类问题时，根据未知运用有限元法求解机械构造应力类问题时，根据未知量和分析量和分析 有三种根本解法：有三

18、种根本解法： 位移法位移法 力法力法 混合法混合法 (1)位移法位移法以节点位移作为根本未知量，经过选择适当的位移函以节点位移作为根本未知量，经过选择适当的位移函数，进展单元的力学特性分析。在节点处建立单元刚度方数，进展单元的力学特性分析。在节点处建立单元刚度方程，再组合成整体刚度矩阵，求解出节点位移后，进而由程，再组合成整体刚度矩阵，求解出节点位移后，进而由节点位移求解出应力。节点位移求解出应力。 位移法优点是比较简单，规律性强，易于编写计算机程位移法优点是比较简单，规律性强，易于编写计算机程序。所以得到广泛运用，其缺陷是精度稍低。序。所以得到广泛运用，其缺陷是精度稍低。 (2)力法力法以节

19、点力作为根本未知量，在节点处建立位移延续方以节点力作为根本未知量，在节点处建立位移延续方程，求解出节点力后，再求解节点位移和单元应力。程，求解出节点力后，再求解节点位移和单元应力。力法的特点是计算精度高。力法的特点是计算精度高。 (3)混合法混合法取一部分节点位移和一部分节点力作为根本未知量，取一部分节点位移和一部分节点力作为根本未知量，建立平衡方程进展求解。建立平衡方程进展求解。单元特性的推导方法单元特性的推导方法 单元刚度矩阵的推导是有限元分析的根本步骤之单元刚度矩阵的推导是有限元分析的根本步骤之一。目前，建立单元刚度矩阵的方法主要有以下四种：一。目前，建立单元刚度矩阵的方法主要有以下四种

20、： 直接刚度法直接刚度法 虚功原理法虚功原理法 能量变分法能量变分法 加权残数法加权残数法1. 直接刚度法直接刚度法 直接刚度法是直接运用物理概念来建立单元直接刚度法是直接运用物理概念来建立单元的有限元方程和分析单元特性的一种方法。这一的有限元方程和分析单元特性的一种方法。这一方法仅能适用于简单外形的单元，如梁单元。但方法仅能适用于简单外形的单元，如梁单元。但它可以协助了解有限元法的物理概念。它可以协助了解有限元法的物理概念。 图图1所示是所示是xoy平面中的一简支梁简图，现以它为例，平面中的一简支梁简图，现以它为例，来阐明用直接刚度法建立单元刚度矩阵的思想和过程。来阐明用直接刚度法建立单元刚

21、度矩阵的思想和过程。图图1平面简支梁元及其计算模型平面简支梁元及其计算模型 梁在横向外载荷可以是集中力或分布力或力矩等作用下产梁在横向外载荷可以是集中力或分布力或力矩等作用下产生弯曲变形，在程度载荷作用下产生线位移。生弯曲变形，在程度载荷作用下产生线位移。 对于该平面简支梁问题：对于该平面简支梁问题：梁上任一点受有三个力的作用：梁上任一点受有三个力的作用： 程度力程度力Fx， 剪切力剪切力Fy ， 弯矩弯矩Mz。相应的位移为：相应的位移为： 程度线位移程度线位移u， 挠度挠度v ， 转角转角 z 。 由上图可见：由上图可见： 程度线位移和程度力向右为正，程度线位移和程度力向右为正， 挠度和剪切

22、力向上为正，挠度和剪切力向上为正， 转角和弯矩逆时针方向为正。转角和弯矩逆时针方向为正。 通常规定：通常规定：为使问题简化，可把图示的梁看作是一个梁单元。为使问题简化，可把图示的梁看作是一个梁单元。如图如图1所示，当令左支承点为节点所示，当令左支承点为节点 i ，右支承点为节点，右支承点为节点 j 时，时，那么该单元的节点位移和节点力可以分别表示为：那么该单元的节点位移和节点力可以分别表示为：称为单元的节点位移列阵。称为单元的节点位移列阵。称为单元的节点力列阵；假设称为单元的节点力列阵；假设 F 为外载荷，那么称为载荷为外载荷，那么称为载荷列阵。列阵。 1-11-2写成矩阵方式为写成矩阵方式为

23、q(e)=ui ,vi , zi ,vj ,uj , zjTui ,vi , zi ,vj ,uj , zjF(e)=Fxi ,Fyi ,Mzi ,Fxj ,Fyj ,MzjTFxi ,Fyi ,Mzi ,Fxj ,Fyj ,Mzj显然，梁的节点力和节点位移是有联络的。在弹性显然，梁的节点力和节点位移是有联络的。在弹性小变形范围小变形范围内，这种关系是线性的，可用下式表示内，这种关系是线性的，可用下式表示 1112131415162122232425263132333435364142434445465152535455566162636 xiyizixjyjzjFkkkkkkFkkkkkkM

24、kkkkkkFkkkkkkFkkkkkkkkkkM46566 iizijjzjuvuvkk ( )( )( )eeeFKq或或1-3b1-3a上式上式(1-3b)称为单元有限元方程，或称为单元刚度方程，它代表称为单元有限元方程，或称为单元刚度方程，它代表了单元的载荷与位移之间或力与变形之间的联络；了单元的载荷与位移之间或力与变形之间的联络；式中，式中，K(e)称为单元刚度矩阵，它是单元的特性矩阵。称为单元刚度矩阵，它是单元的特性矩阵。 对于图对于图1所示的平面梁单元问题，利用资料力学中的杆所示的平面梁单元问题，利用资料力学中的杆件受力与变形间的关系及叠加原理，可以直接计算出单元件受力与变形间的

25、关系及叠加原理，可以直接计算出单元刚度矩阵刚度矩阵K(e)中的各系数中的各系数 kst( s, t = i, j ) 的数值的数值2. 虚功原理法虚功原理法下面以平面问题中的三角形单元为例，阐明利用虚功原理法来建下面以平面问题中的三角形单元为例，阐明利用虚功原理法来建立单元刚度矩阵的步骤。立单元刚度矩阵的步骤。如前所述，将一个延续的弹性体分割为一定外形和数量的单元，如前所述，将一个延续的弹性体分割为一定外形和数量的单元，从而使延续体转换为有限个单元组成的组合体。单元与单元之间仅经从而使延续体转换为有限个单元组成的组合体。单元与单元之间仅经过节点连结，除此之外再无其他连结。也就是说，一个单元上的

26、只能过节点连结，除此之外再无其他连结。也就是说，一个单元上的只能经过节点传送到相邻单元。经过节点传送到相邻单元。 从分析对象的组合体中任取从分析对象的组合体中任取一个三角形单元：一个三角形单元：设其编号为设其编号为 e ，三个节点的编号为三个节点的编号为i、j、m，在定义的坐标系在定义的坐标系 xoy 中，节点坐中，节点坐标分别为标分别为(x j , y j)、(xi , y i)、(xm, ym)，如图，如图2所示。所示。图图2三节点三角形单元三节点三角形单元由弹性力学平面问题的特点可知，单元每个节点有两个位移分由弹性力学平面问题的特点可知，单元每个节点有两个位移分量，即每个单元有量，即每个

27、单元有6个自在度，相应有个自在度，相应有6个节点载荷，写成矩阵方式，个节点载荷，写成矩阵方式，即即 单元节点载荷矩阵：单元节点载荷矩阵：F(e)=Fxi ,Fyi ,Fxj ,Fyj ,Fxm ,FymT单元节点位移矩阵：单元节点位移矩阵：q(e)=ui ,vi ,uj ,vj ,um ,vmT图图2三节点三角形单元三节点三角形单元(1)设定位移函数设定位移函数 按照有限元法的根本思想：首先需设定一种函数来近似表达单元按照有限元法的根本思想：首先需设定一种函数来近似表达单元内部的实践位移分布，称为位移函数，或位移方式。内部的实践位移分布，称为位移函数，或位移方式。 三节点三角形单元有三节点三角

28、形单元有6个自在度，可以确定个自在度，可以确定 6个待定系数，故三角形单元个待定系数，故三角形单元的位移函数为的位移函数为 1-4式式(1-4)为线性多项式，称为线性位移函数，相应的单元称为线性单元。为线性多项式，称为线性位移函数，相应的单元称为线性单元。u=u(x,y)= 1+ 2x+ 3yv=v(x,y)= 4+ 5x+ 6y 1234561 0 0 00 0 0 1 uxydsvx y 上式上式(5-5)也可用矩阵方式表示，即也可用矩阵方式表示，即 式中，式中，d为单元内恣意点的位移列阵。为单元内恣意点的位移列阵。 1-5由于节点由于节点 i、j、m 在单元上，它们的位移自然也就满足位移

29、在单元上，它们的位移自然也就满足位移函数式函数式(1-4)。设三个节点的位移值分别为设三个节点的位移值分别为( ui, vi)、( uj, vj )、( um, vm )，将节将节点位移和节点坐标代入式点位移和节点坐标代入式(1-4)，得，得 123123123iiijjjmmmuxyuxyuxy456456456iiijjjmmmvxyvxyvxy( )12345600000000010002000000eijmiijmiijnjijmjijmmijmmaaaubbbvcccuaaavbbbucccv1-6111112221iijjijjmmijimjimmmxyxyx yx yx yx y

30、x yx yxy 式中1-7由上可知，共有6个方程，可以求出6个待定系数。解方程，求得各待定系数和节点位移之间的表达式为 为三角形单元的面积。其中： , , , , , , ijmmjjmiimmijjiijmjmimijimjjimmjiax yx yax yx yax yx ybyybyybyycxxcxxcxx1-8将式(1-7)及式(1-8)、式(1-9)代入式(1-6)中，得到 1-91-10式中，矩阵N 称为单元的形函数矩阵； 为单元节点位移列阵。其中， 为单元的形函数，它们反映单元内位移的 分布形状，是x, y 坐标的延续函数，且有 ( ) eq, , ijmNNN222iiii

31、jjjjmmmmNab xc yNab xc yNab xc y1-11式(1-10)又可以写成, , ,iijjmmiii i j miijjmmiii i j muN uN uN uN uvN vN vN vN v1-12上式清楚地表示了单元内恣意点位移可由节点位移插值求出。 (2) 利用几何方程由位移函数求应变根据弹性力学的几何方程 ，线应变 剪切应变 那么应变列阵可以写成 /, / ,xyuxuy / ,xyuyux ( )( ) 0001 0002iixijmeejyijmjxyiijjmmmmuuvxbbbuucccB qvycbcbcbuuuyxv式中，B称为单元应变矩阵，它是仅

32、与单元几何尺寸有关的常量矩阵，即 1-13 00010002ijmijmiijjmmbbbBccccbcbcb1-14 上述方程(1-13)称为单元应变方程，它的意义在于： 单元内恣意点的应变分量亦可用根本未知量即节点位移分量来表示。 (3)利用广义虎克定律求出单元应力方程根据广义虎克定律，对于平面应力问题1()1()12(1)xxyyyxxyxyxyEEGE上式1-15也可写成 ( )( )eeD1-151-16式中， 为应力列阵； D 称为弹性力学平面问题的弹性矩阵，并有 ,Txyxy 2101011002ED ( )( )( )( )exeeeyxyDDBq那么有如下单元应力方程由式(1

33、-18)可求单元内恣意点的应力分量，它也可用根本未知量即节点位移分量来表示。1-171-18(4)由虚功原理求单元刚度矩阵 根据虚功原理，当弹性构造遭到外载荷作用途于平衡形状时，在恣意给出的微小的虚位移上，外力在虚位移上所做的虚功 AF等于构造内应力在虚应变上所存储的虚变情势能 A ，即FAA设处于平衡形状的弹性构造内任一单元发生一个微小的虚位移，那么单元各节点的虚位移 为 ( )eq ( ), , , , , eTiijjmmquvuvuv1-201-211-19那么单元内部必定产生相应的虚应变，故单元内任一点的虚应变 为 ( )e ( ), , eTxyxy显然，虚应变和虚位移之间关系为

34、( )( )eeBq设节点力为 ( ), , , , , TexiyixjyjxmymFFFFFFF那么外力虚功为 ( )( )TeeFAqF1-241-221-23单元内的虚变情势能为 ( ) TeVAdv根据虚功原理 ( )( )( )TTeeeVqFdv ( )( )eeDDBq ( )( )( )TTTeeeTBqBq由于1-261-25代人式(1-25)，那么有 ( )( )( )( )TTeeTeeVqFBqDBqdv式中， ，均与坐标 x, y 无关，故可以从积分符号中提出，可得： ( )Teq ( )eq ( )( )( )( )TeeeeVFBDB dvqKq其中，单元刚度矩

35、阵 ( )eTVKBDB dv1-27式(1-27)称为单元有限元方程，或称单元刚度方程，其中 是单元刚度矩阵。 ( )eK1-28由于三角形单元是常应变单元，其应变矩阵B 、弹性矩阵D均为常量，而 ，所以式(1-28)可以写成 Vdvtdxdyt ( )eTKtBDB 1-29式中，t 为三角形单元的厚度； 为三角形单元的面积。 ( )( ) eiiijimejijjjmmimjmmKKKKKKKKKK 对于图2所示的三角形单元，将D 及B代入式(1-28)，可以得到单元刚度为 1-30式中: K为66阶矩阵，其中每个子矩阵为22阶矩阵，由下式给出 211 22114(1) c22 ( ,

36、,)rsrsrsrsrsrsrsrsrsb bc cb cc bEtKc bb ccb br si j m 1-31按照力学的普通说法，任何一个实践形状的弹性体的总位能是这个系统从实践形状运动到某一参考形状(通常取弹性体外载荷为零时形状为参考形状)时它的一切作用力所做的功。弹性体的总位能 是一个函数的函数，即泛函，位移是泛函的允许函数。从能量原理思索，变形弹性体受外力作用途于平衡形状时，在很多能够的变外形状中，使总位能最小的就是弹性体的真正变形，这就是最小位能原理。用变分法求能量泛函的极值方法就是能量变分原理。能量变分原理除了可解机械构造位移场问题以外，还扩展到求解热传导、电磁场、流膂力学等延

37、续性问题。3. 能量变分原理法该方法是将假设的场变量的函数(称为试函数)引入问题的控制方程式及边境条件，利用最小二乘法等方法使残差最小，便得到近似的场变量函数方式。该方法的优点是不需求建立要处理问题的泛函式，所以，即使没有泛函表达式也能解题。 4. 加权残数法 有限元解的收敛性有限元解的收敛性 有限元解是近似解有限元解是近似解 近似解能否收敛于真实解、近似解收敛速度、近似解的稳近似解能否收敛于真实解、近似解收敛速度、近似解的稳定性定性 近似解的收敛条件近似解的收敛条件: 1.完备性准那么必要条件完备性准那么必要条件 试探函数插值函数的次数试探函数插值函数的次数m不小于场函数的最高不小于场函数的

38、最高可导阶数可导阶数 2.协调性准那么充分条件协调性准那么充分条件 试探函数在试探函数在m-1次延续可导。次延续可导。有限元分析的误差有限元分析误差建模误差计算误差离散误差物理离散误差几何离散误差边境条件误差单元外形误差舍入误差截断误差插值函数与真实函数之间的差别1.减小单元特征尺寸，称为h法2.提高插值函数的阶次，称为p法单元组合体与求解对象几何外形的差别1.网格部分加密2.选用边或面上带有节点的单元边境条件的复杂性1.准确测定，完善模型2.细分边境网格单元严重畸变而退化细分部分网格或者控制调整关键区域的网格数据储存计算方法、解题性质、解题规模留意网格的划分留意网格的划分 选择适宜的解算方法

39、选择适宜的解算方法 控制解题的规模控制解题的规模减少运算次数，降低解题规模减少运算次数，降低解题规模选择适宜的解算方法，控制解题规模选择适宜的解算方法，控制解题规模资料成形中的非线性问题 1.资料非线性 资料本构方程非线性 弹塑性、刚弹性、刚黏塑性、黏弹塑性 2.几何非线性 3.边境非线性2.2 有限差分法根底 一种直接将微分问题转变成代数问题的近似数值解法。 根本思想 数值微分法 是把延续的定解区域划分成差分网格，用有限个节点替代原延续求解域。 把原方程和定解条件中的微商用差商来近似 把原微分方程和定解条件近似地用代数方程组替代，即有限差分方程差分网格通常为矩形在边境不规那么或者外形复杂时精

40、度降低有限元网格有限差分网格差分概念自变量x的解析函数 y =f (x)，那么有：dx,dy自变量和函数的微分 函数对自变量的一阶导数 函数对自变量的一阶差商00()( )limlimxxdyyf xxf xdxxx yxdydx差商差商差分方向 向前差分 向后差分 中心差分()( )yf xxf x ( )()yf xf xx 11()()22yf xxf xx 差商的截断误差差商的截断误差 将函数将函数f (x+x)按按Taylor级数展开级数展开 向前向前 向后向后 中心中心234()()()( )( )( )( ) 0() )1!2!3!xxxf xxf xf xfxfxx 23()(

41、 )()( )( )( ) 0() )2!3!f xxf xxxf xfxfxxx 23( )()()( )( )( )0() )2!3!f xf xxxxfxfxfxxx23()()()22( )( )0() )2! 3!xxf xf xxfxfxxx 二阶中心差商 通常采用向前差商的向后差商 截断误差与(x)2 同一数量级 一阶向前差商 一阶向后差商 一阶计算精度 一阶中心差商 二阶中心差商 二阶计算精度2222232()2 ( )()()2()( )( )0() )4!yf xxf xf xxxxyxfxfxxx 我们在弹性体上，用相隔等间距h而平行于坐标轴的两组平行线织成正方形网格，x

42、=y=h，如图。 设f=f(x,y)为弹性体内的某一个延续函数。该函数在平行于x轴的一根网线上，如在-上,它只随x坐标的改动而变化。在临近结点处，函数f可展为泰勒级数如下：.)(! 31)(! 21)(3003320022000 xxxfxxxfxxxfff 我们将只思索分开结点充分近的那些结点，即x-x0充分小。于是可不计x-x0的三次及更高次幂的各项，那么上式简写为：在结点，x=x0-h, 在结点1, x=x0+h，代入(b) 得：)(.)(! 21)(20022000bxxxfxxxfff)(.20222003cxfhxfhff)(.20222001dxfhxfhff联立(c),(d),

43、解得差分公式： )11(2310hffxf)21 (22031022hfffxf同理，在网线-上可得到差分公式)41 (22)31 (2042022420hfffyfhffyf从而可导出其它的差分公式如下：)()(461)()(241)()(461)()(41121042040448765432104022411931040447586022fffffhyffffffffffhyxffffffhxfffffhyxf 相隔2h的两结点处的函数值来表示中间结点处的一阶导数值，可称为中点导数公式。 以相邻三结点处的函数值来表示一个端点处的一阶导数值，可称为端点导数公式。 中点导数公式与端点导数公式相

44、比，精度较高。由于前者反映了结点两边的函数变化，而后者却只反映了结点一边的函数变化。边境元法简介边境元法简介 边境元法边境元法boundary element method 一种结合有限元法和边境积分法开展起来的一种新数值方一种结合有限元法和边境积分法开展起来的一种新数值方法法 只在定义域的边境上划分单元，用满足控制方程的函数去只在定义域的边境上划分单元，用满足控制方程的函数去逼近边境条件。逼近边境条件。 适用于应力薄板、流膂力学、声场、电磁场等问题适用于应力薄板、流膂力学、声场、电磁场等问题边境元法根本思想边境元法根本思想 以微分控制方程的根本解为权函数，利用加权余量法将区以微分控制方程的根

45、本解为权函数，利用加权余量法将区域积分转化为边境积分，并结合求解域边境的离散，构建域积分转化为边境积分，并结合求解域边境的离散，构建基于边境单元的代数方程组，然后进展计算求解基于边境单元的代数方程组，然后进展计算求解 以定义在边境上的边境积分方程为控制方程，经过对边境以定义在边境上的边境积分方程为控制方程，经过对边境元插值离散，化为代数方程组求解元插值离散，化为代数方程组求解 降低了问题的维数，从而显著降低了自在度数降低了问题的维数，从而显著降低了自在度数 边境的离散比区域的离散方便得多，可用较简单的单元准边境的离散比区域的离散方便得多，可用较简单的单元准确地模拟边境外形，最终得到阶数较低的线

46、性代数方程组确地模拟边境外形，最终得到阶数较低的线性代数方程组加权余量法简介加权余量法简介一种直接从所需求解的微分方程及边境条件出发，寻求边值问题近似一种直接从所需求解的微分方程及边境条件出发，寻求边值问题近似解的数学方法。解的数学方法。从静力开展到了动力、稳定、资料非线性和几何非线性等各方面。从静力开展到了动力、稳定、资料非线性和几何非线性等各方面。在求解域上建立一个试函数在求解域上建立一个试函数 试函数由完备函数集的子集构成。已被采用过的试函数有幂级数、试函数由完备函数集的子集构成。已被采用过的试函数有幂级数、三角级数、样条函数、贝赛尔函数、切比雪夫和勒让德多项式等等。三角级数、样条函数、

47、贝赛尔函数、切比雪夫和勒让德多项式等等。试函数与真实解之间的偏向，即余量内部和边境试函数与真实解之间的偏向，即余量内部和边境引入权函数，定义消除余量的条件引入权函数，定义消除余量的条件加权余量法就是一种定义近似解与真解之间余量，并设法使其最小的加权余量法就是一种定义近似解与真解之间余量，并设法使其最小的方法。方法。设问题的控制微分方程为:在V域内 在S边境上 式中 : L、B分别为微分方程和边境条件中的微分算子； f、g 为与未知函数u无关的知函数域值； u为问题待求的未知函数。( )0B ug( )0L uf当利用加权余量法求近似解时，首先在求解域上建立一个试函数 ，普通具有如下方式：u 1

48、niiiuC NNC式中： 待定系数，也可称为广义坐标；iC取自完备函数集的线性无关的基函数。iN由于 一 般只是待求函数u的近似解，因此代入后将得不到满足，假设记：u ( )( )IBRL ufRB ug在V域内在S边境上显然 反映了试函数与真实解之间的偏向，它们分别称做内部和边境余量。BIRR、 假设在域V内引入内部权函数 ，在边境S上引入边境权函数那么可建立n个消除余量的条件，普通可表示为：IWBW0(1,2, )IiIBiBVSW R dVWR dSin不同的权函数 和 反映了不同的消除余量的准那么。从上式可以得到求解待定系数矩阵C的代数方程组。一经解得待定系数，由式(5.1.3)即可

49、得所需求解边值问题的近似解。BiWIiWIiW 由于试函数 的不同，余量 和 可有如下三种情况，依此加权余量法可分为：1内部法试函数满足边境条件，也即 此时消除余量的条件成为:2边境法试函数满足控制方程，也即 此时消除余量的条件为:( )0BRB ug( )0IRL uf0(1,2, )IiIVW R dVin0(1,2, )BiBSW R dSinu IRBR3混合法 试函数不满足控制方程和边境条件，此时消除余量的条件为: 显然，混合法对于试函数的选取最方便，但在一样精度条件下，任务量最大。对内部法和边境法必需使基函数事先满足一定条件，这对复杂构造分析往往有一定困难，但试函数一经建立，其任务量较小。0(1,2, )IiIBiBVSW R dVW R dSin边境元法特点边境元法特点 1 前处置任务量小前处置任务量小 2 解算规模小解算规模小 3 求解奇特性问题时计算精度高求解奇特性问题时计算精度高 4 在载荷集中和半无限域等问题上有优势在载荷集中和半无限域等问题上有优势 相对于有限元法，边境元法开展较慢，相对于有限元法，边境元法开展较慢，有限元法处理应