有限元与有限差分法基础

1、1第二讲第二讲 有限元与有限差分法基础有限元与有限差分法基础 CAE的工具：的工具： 有限元法（有限元法（FEM）、有限差分法（有限差分法（FDM）、边界元法边界元法（BEM）、有限体积法（、有限体积法（FVM）、无网格法等等）、无网格法等等 在材料成形的在材料成形的CAE中主要使用的是中主要使用的是有限元法和有限差分法有限元法和有限差分法2“ ” 的的早在早在20世纪世纪40年代初期就年代初期就有人提出，但真正用于工程中则是有人提出，但真正用于工程中则是电子计算机电子计算机出现以后。出现以后。 “ ” 这一名称是这一名称是1960年美国的年美国的克拉夫克拉夫（Clough，R.W.）在一篇题

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

3、，由于它在理论上的通用性，因而它可用于解决工程中的许多问题。理论上的通用性，因而它可用于解决工程中的许多问题。目前，它可以解决几乎所有的目前，它可以解决几乎所有的连续介质连续介质和和场场的问题，的问题，包括热传导、电磁场、流体动力学、地质力学、原子工程包括热传导、电磁场、流体动力学、地质力学、原子工程和生物医学等方面的问题。和生物医学等方面的问题。 中，从齿轮、轴、轴承等通用零部件到机中，从齿轮、轴、轴承等通用零部件到机床、汽车、飞机等床、汽车、飞机等复杂结构复杂结构的应力和变形分析（包括热应的应力和变形分析（包括热应力和热变形分析）。力和热变形分析）。不仅可以解决工程中的线性问题、非线性问不

4、仅可以解决工程中的线性问题、非线性问题，而且对于各种不同性质的固体材料，如各向同性和各题，而且对于各种不同性质的固体材料，如各向同性和各向异性材料，粘弹性和粘塑性材料以及流体均能求解；向异性材料，粘弹性和粘塑性材料以及流体均能求解；对于工程中最有普遍意义的对于工程中最有普遍意义的非稳态问题非稳态问题也能求解。也能求解。42.1 有限元法基础有限元法基础基本思想：基本思想： 将一个连续求解域（对象）将一个连续求解域（对象）离散离散（剖分）成（剖分）成有限有限个形状个形状简单的简单的子域（单元）子域（单元） 利用有限个利用有限个节点节点将各将各子域子域连接起来连接起来 在给定的在给定的初始条件初始

5、条件和和边界条件边界条件下进行综合计算求解，从下进行综合计算求解，从而获得对复杂工程问题的而获得对复杂工程问题的近似数值解近似数值解 5物理系统举例 几何体几何体 载荷载荷 物理系统物理系统结构结构热热电磁电磁6有限元模型真实系统真实系统有限元模型有限元模型 有限元模型有限元模型 是真实系统理想化的数学抽象是真实系统理想化的数学抽象。定义定义7自由度（DOFs）自由度自由度(DOFs) 用于描述一个物理场的响应特性用于描述一个物理场的响应特性。结构结构 DOFs 结构结构 位移位移 热热 温度温度 电电 电位电位 流体流体 压力压力 磁磁 磁位磁位 方向方向 自由度自由度ROTZUYROTYU

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

7、单元信息是通过单元之间的公共节点传递的。信息是通过单元之间的公共节点传递的。分离但节点重叠的单元分离但节点重叠的单元A和和B之间没有信息传递之间没有信息传递（需进行节点合并处理）（需进行节点合并处理）具有公共节点的单元具有公共节点的单元之间存在信息传递之间存在信息传递 .AB.AB.1 node2 nodes每个单元的特性是通过一些线性方程式来描述的。每个单元的特性是通过一些线性方程式来描述的。作为一个整体，单元形成了整体结构的数学模型。作为一个整体，单元形成了整体结构的数学模型。10节点和单元节点和单元节点自由度是随连接该节点节点自由度是随连接该节点 单元类型单元类型 变化的。变化的。JII

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

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

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

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

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

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

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

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

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

17、据的原理是节点变形协调条件节点变形协调条件和和平衡条平衡条件件。 204. 确定约束条件确定约束条件 由上述所形成的由上述所形成的是是一组线性代数方程一组线性代数方程，在求解之前，必修根据具体情况分析，确定在求解之前，必修根据具体情况分析，确定的的边界约束条件边界约束条件，并对，并对这些方程这些方程进行适当修正。进行适当修正。215. 有限元方程求解有限元方程求解 应用应用求解机械结构应力类问题时，根据求解机械结构应力类问题时，根据未知未知量量和和分析分析 有有三种基本解法三种基本解法： 位移法位移法 力法力法 混合法混合法22 (1)位移法位移法以以节点位移节点位移作为作为基本未知量基本未知

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

19、力后，再求解节点位移和单元应力。后，再求解节点位移和单元应力。力法的特点力法的特点是计算精度高。是计算精度高。 (3)混合法混合法取一部分取一部分节点位移节点位移和一部分和一部分节点力节点力作为作为基本未知量基本未知量，建立建立平衡方程平衡方程进行求解。进行求解。23单元特性的推导方法单元特性的推导方法 的推导是的推导是的的基本步骤之基本步骤之一一。目前，建立。目前，建立单元刚度矩阵单元刚度矩阵的方法主要有以下的方法主要有以下四种四种： 直接刚度法直接刚度法 虚功原理法虚功原理法 能量变分法能量变分法 加权残数法加权残数法241. 直接刚度法直接刚度法 是直接应用物理概念来建立是直接应用物理概

20、念来建立单元的有限元单元的有限元方程方程和分析单元特性的一种方法。和分析单元特性的一种方法。这一方法这一方法仅能适用于仅能适用于简单形状的单元简单形状的单元，如，如梁单元梁单元。但它可以帮助理解。但它可以帮助理解有限元有限元法法的物理概念。的物理概念。 图图1所示是所示是xoy平面中的一平面中的一，现以它为例，现以它为例，来说明用来说明用建立建立单元刚度矩阵单元刚度矩阵的思想和过程。的思想和过程。图图1平面简支梁元及其计算模型平面简支梁元及其计算模型25 梁在梁在横向外载荷横向外载荷（可以是集中力或分布力或力矩等）作用下产（可以是集中力或分布力或力矩等）作用下产生生弯曲变形弯曲变形，在，在水平

21、载荷水平载荷作用下产生作用下产生线位移线位移。 对于对于该平面简支梁该平面简支梁问题：问题：梁上任一点梁上任一点受有受有三个力三个力的作用：的作用： 水平力水平力Fx， 剪切力剪切力Fy ， 弯矩弯矩Mz。相应的位移相应的位移为：为： 水平线位移水平线位移u， 挠度挠度v ， 转角转角 z 。 由由上图上图可见：可见： 水平线位移水平线位移和和水平力水平力向右为正，向右为正， 挠度挠度和和剪剪切力切力向上为正，向上为正， 转角转角和和弯矩弯矩逆时针方向为正。逆时针方向为正。 通常规定：通常规定：26为使为使问题简化问题简化，可把，可把图示的梁图示的梁看作是一个看作是一个梁单元梁单元。如如图图1

22、所示，当令所示，当令左支承点左支承点为为节点节点 i ，右支承点右支承点为为节点节点 j 时，时，则则该单元该单元的的节点位移节点位移和和节点力节点力可以分别表示为：可以分别表示为：称为称为单元的节点位移列阵单元的节点位移列阵。称为称为；若；若 F 为为，则称为，则称为载荷列阵载荷列阵。 （1-1）（1-2）写成写成矩阵形式矩阵形式为为q(e)=ui ,vi , zi ,vj ,uj , zjTui ,vi , zi ,vj ,uj , zjF(e)=Fxi ,Fyi ,Mzi ,Fxj ,Fyj ,MzjTFxi ,Fyi ,Mzi ,Fxj ,Fyj ,Mzj27显然，显然，梁的节点力梁的

23、节点力和和节点位移节点位移是有联系的。在是有联系的。在弹性弹性小变形范围小变形范围内，内，这种关系这种关系是线性的，可用是线性的，可用下式下式表示表示 1112131415162122232425263132333435364142434445465152535455566162636 xiyizixjyjzjFkkkkkkFkkkkkkMkkkkkkFkkkkkkFkkkkkkkkkkM46566 iizijjzjuvuvkk ( )( )( )eeeFKq或或（1-3b）（1-3a）28上上式式(1-3b)称为称为，或称为，或称为单元刚度方单元刚度方程程，它代表了，它代表了单元的载荷单元的

24、载荷与与位移位移之间（或力与变形之间）之间（或力与变形之间）的联系；的联系；式中，式中，K(e)称为称为单元刚度矩阵单元刚度矩阵，它是单元的特性矩阵。，它是单元的特性矩阵。 对于对于图图1所示的所示的平面梁单元问题平面梁单元问题，利用材料力学中的杆，利用材料力学中的杆件受力与变形间的关系及叠加原理，可以直接计算出件受力与变形间的关系及叠加原理，可以直接计算出单元单元刚度矩阵刚度矩阵K(e)中的中的各系数各系数 kst( s, t = i, j ) 的数值的数值292. 虚功原理法虚功原理法下面以下面以平面问题平面问题中的中的三角形单元三角形单元为例，说明利用为例，说明利用虚功原理法虚功原理法来

25、建来建立立单元刚度矩阵单元刚度矩阵的步骤。的步骤。如前所述，将一个如前所述，将一个连续的弹性体连续的弹性体分割为分割为一定形状一定形状和和数量数量的的单元单元，从而使从而使连续体连续体转换为转换为有限个单元有限个单元组成的组成的组合体组合体。单元单元与与单元单元之间仅通之间仅通过过节点节点连结，除此之外再无其他连结。也就是说，一个连结，除此之外再无其他连结。也就是说，一个单元单元上的只能上的只能通过通过节点节点传递到传递到相邻单元相邻单元。 从分析对象的从分析对象的组合体组合体中任取中任取一个一个三角形单元三角形单元：设设其编号其编号为为 e ，三个节点三个节点的的编号编号为为i、j、m，在定

26、义的在定义的坐标系坐标系 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

27、图图2三节点三角形单元三节点三角形单元31(1)设定位移函数设定位移函数 按照按照有限元法有限元法的的基本思想基本思想：首先需设定：首先需设定一种函数一种函数来近似表达单元来近似表达单元内部的内部的实际位移分布实际位移分布，称为，称为位移函数位移函数，或，或位移模式位移模式。 三节点三角形单元三节点三角形单元有有6个自由度个自由度，可以确定，可以确定 6个待定系数个待定系数，故，故三角三角形单元形单元的的位移函数位移函数为为 （1-4）式式(1-4)为为线性多项式线性多项式，称为，称为线性位移函数线性位移函数，相应的单元相应的单元称为称为线性单元线性单元。u=u(x,y)= 1+ 2x+ 3y

28、v=v(x,y)= 4+ 5x+ 6y32 1234561 0 0 00 0 0 1 uxydsvx y 上上式式(5-5)也可用也可用矩阵形式矩阵形式表示，即表示，即 式中，式中，d为为单元内任意点单元内任意点的的位移列阵位移列阵。 （1-5）33由于由于节点节点 i、j、m 在在单元单元上，上，它们的位移它们的位移自然也就满足自然也就满足位移位移函数式函数式(1-4)。设设三个节点三个节点的的位移值位移值分别为分别为( ui, vi)、( uj, vj )、( um, vm )，将，将节节点位移点位移和和节点坐标节点坐标代入代入式式(1-4)，得，得 123123123iiijjjmmmu

29、xyuxyuxy456456456iiijjjmmmvxyvxyvxy34( )12345600000000010002000000eijmiijmiijnjijmjijmmijmmaaaubbbvcccuaaavbbbucccv（1-6）111112221iijjijjmmijimjimmmxyxyx yx yx yx yx yx yxy 式中（1-7）由上可知，共有6个方程，可以求出6个待定系数。解方程，求得各待定系数和节点位移之间的表达式为 为三角形单元的面积。其中： 35, , , , , , ijmmjjmiimmijjiijmjmimijimjjimmjiax yx yax yx

30、yax yx ybyybyybyycxxcxxcxx（1-8）将式(1-7)及式(1-8)、式(1-9)代入式(1-6)中，得到 （1-9）（1-10）36式中，矩阵N 称为单元的形函数矩阵； 为单元节点位移列阵。其中， 为单元的形函数，它们反映单元内位移的 分布形态，是x, y 坐标的连续函数，且有 ( ) eq, , ijmNNN222iiiijjjjmmmmNab xc yNab xc yNab xc y（1-11）式(1-10)又可以写成, , ,iijjmmiii i j miijjmmiii i j muN uN uN uN uvN vN vN vN v（1-12）上式清楚地表示了

31、单元内任意点位移可由节点位移插值求出。 37(2) 利用几何方程由位移函数求应变根据弹性力学的几何方程 ，线应变 剪切应变 则应变列阵可以写成 /, / ,xyuxuy / ,xyuyux ( )( ) 0001 0002iixijmeejyijmjxyiijjmmmmuuvxbbbuucccB qvycbcbcbuuuyxv式中，B称为单元应变矩阵，它是仅与单元几何尺寸有关的常量矩阵，即 （1-13）38 00010002ijmijmiijjmmbbbBccccbcbcb（1-14） 上述方程(1-13)称为单元应变方程，它的意义在于： 单元内任意点的应变分量亦可用基本未知量即节点位移分量来

32、表示。 39(3)利用广义虎克定律求出单元应力方程根据广义虎克定律，对于平面应力问题1()1()12(1)xxyyyxxyxyxyEEGE上式（1-15）也可写成 ( )( )eeD（1-15）（1-16）式中， 为应力列阵； D 称为弹性力学平面问题的弹性矩阵，并有 ,Txyxy 40 2101011002ED ( )( )( )( )exeeeyxyDDBq则有如下单元应力方程由式(1-18)可求单元内任意点的应力分量，它也可用基本未知量即节点位移分量来表示。（1-17）（1-18）41(4)由虚功原理求单元刚度矩阵 根据虚功原理，当弹性结构受到外载荷作用处于平衡状态时，在任意给出的微小的

33、虚位移上，外力在虚位移上所做的虚功 AF等于结构内应力在虚应变上所存储的虚变形势能 A ，即FAA设处于平衡状态的弹性结构内任一单元发生一个微小的虚位移，则单元各节点的虚位移 为 ( )eq ( ), , , , , eTiijjmmquvuvuv（1-20）（1-21）（1-19）则单元内部必定产生相应的虚应变，故单元内任一点的虚应变 为 ( )e ( ), , eTxyxy42显然，虚应变和虚位移之间关系为 ( )( )eeBq设节点力为 ( ), , , , , TexiyixjyjxmymFFFFFFF则外力虚功为 ( )( )TeeFAqF（1-24）（1-22）（1-23）单元内的

34、虚变形势能为 ( ) TeVAdv43根据虚功原理 ( )( )( )TTeeeVqFdv ( )( )eeDDBq ( )( )( )TTTeeeTBqBq因为（1-26）（1-25）代人式(1-25)，则有 ( )( )( )( )TTeeTeeVqFBqDBqdv式中， ，均与坐标 x, y 无关，故可以从积分符号中提出，可得： ( )Teq ( )eq44 ( )( )( )( )TeeeeVFBDB dvqKq其中，单元刚度矩阵 ( )eTVKBDB dv（1-27）式(1-27)称为单元有限元方程，或称单元刚度方程，其中 是单元刚度矩阵。 ( )eK（1-28）因为三角形单元是常应

35、变单元，其应变矩阵B 、弹性矩阵D均为常量，而 ，所以式(1-28)可以写成 Vdvtdxdyt ( )eTKtBDB （1-29）45式中，t 为三角形单元的厚度； 为三角形单元的面积。 ( )( ) eiiijimejijjjmmimjmmKKKKKKKKKK 对于图2所示的三角形单元，将D 及B代入式(1-28)，可以得到单元刚度为 （1-30）式中: K为66阶矩阵，其中每个子矩阵为22阶矩阵，由下式给出 211 22114(1) c22 ( , ,)rsrsrsrsrsrsrsrsrsb bc cb cc bEtKc bb ccb br si j m （1-31）46按照力学的一般说

36、法，任何一个实际状态的是这个系统从实际状态运动到某一参考状态(通常取弹性体外载荷为零时状态为参考状态)时。弹性体的总位能 是一个函数的函数，即泛函，位移是泛函的容许函数。从考虑，变形弹性体受外力作用处于平衡状态时，在很多可能的变形状态中，使总位能最小的就是弹性体的真正变形，这就是最小位能原理。用求能量泛函的极值方法就是能量变分原理。能量变分原理除了可解机械结构位移场问题以外，还扩展到求解热传导、电磁场、流体力学等连续性问题。3. 能量变分原理法47该方法是将假设的场变量的函数(称为试函数)引入问题的控制方程式及边界条件，利用最小二乘法等方法使残差最小，便得到近似的场变量函数形式。该方法的优点是

37、不需要建立要解决问题的泛函式，所以，即使没有泛函表达式也能解题。 4. 加权残数法48 有限元解的收敛性有限元解的收敛性有限元解是近似解近似解是否收敛于真实解、近似解收敛速度、近似解的稳定性 近似解的收敛条件:1.完备性准则（必要条件）试探函数（插值函数）的次数（m）不小于场函数的最高可导阶数2.协调性准则（充分条件）试探函数在m-1次连续可导。49有限元分析的误差有限元分析误差建模误差计算误差离散误差物理离散误差几何离散误差边界条件误差单元形状误差舍入误差截断误差插值函数与真实函数之间的差异1.减小单元特征尺寸，称为减小单元特征尺寸，称为h法法2.提高插值函数的阶次，称为提高插值函数的阶次，

38、称为p法法单元组合体与求解对象几何形状的差异1.网格局部加密网格局部加密2.选用边或面上带有节点的单元选用边或面上带有节点的单元边界条件的复杂性1.准确测定，完善模型准确测定，完善模型2.细分边界网格细分边界网格单元严重畸变而退化细分局部网格或者控制调整关键区域的网格细分局部网格或者控制调整关键区域的网格数据储存计算方法、解题性质、解题规模注意网格的划分注意网格的划分 选择合适的解算方法选择合适的解算方法 控制解题的规模控制解题的规模减少运算次数，降低解题规模减少运算次数，降低解题规模选择合适的解算方法，控制解题规模选择合适的解算方法，控制解题规模50材料成形中的非线性问题 1.材料非线性材料

39、本构方程非线性 弹塑性、刚弹性、刚黏塑性、黏弹塑性 2.几何非线性 3.边界非线性512.2 有限差分法基础 一种直接将微分问题转变成代数问题的近似数值解法。 基本思想基本思想 数值微分法数值微分法 是把连续的定解区域划分成差分网格，用有限个节点代替原连续求解域。 把原方程和定解条件中的微商用差商来近似 把原微分方程和定解条件近似地用代数方程组代替，即有限差分方程52差分网格通常为矩形在边界不规则或者形状复杂时精度降低有限元网格有限差分网格53差分概念自变量x的解析函数 y =f (x)，则有：dx,dy自变量和函数的微分 函数对自变量的一阶导数 函数对自变量的一阶差商00()( )limli

40、mxxdyyf xxf xdxxx yxdydx差商差商54差分方向 向前差分 向后差分 中心差分()( )yf xxf x ( )()yf xf xx 11()()22yf xxf xx 55差商的截断误差差商的截断误差 将函数将函数f (x+x)按按Taylor级数展开级数展开 向前向前 向后向后 中心中心234()()()( )( )( )( ) 0() )1!2!3!xxxf xxf xf xfxfxx23()( )()( )( )( ) 0() )2!3!f xxf xxxf xfxfxxx 23( )()()( )( )( )0() )2!3!f xf xxxxfxfxfxxx23

41、()()()22( )( )0() )2! 3!xxf xf xxfxfxxx56 二阶中心差商 通常采用向前差商的向后差商 截断误差与(x)2 同一数量级 一阶向前差商 一阶向后差商 一阶计算精度 一阶中心差商 二阶中心差商 二阶计算精度2222232()2 ( )()()2()( )( )0() )4!yf xxf xf xxxxyxfxfxxx57 我们在弹性体上，用相隔等间距h而平行于坐标轴的两组平行线织成正方形网格，x=y=h，如图。 设f=f(x,y)为弹性体内的某一个连续函数。该函数在平行于x轴的一根网线上，如在-上,它只随x坐标的改变而变化。在邻近结点处，函数f可展为泰勒级数如

42、下：.)(! 31)(! 21)(3003320022000 xxxfxxxfxxxfff58 我们将只考虑离开结点充分近的那些结点，即（x-x0）充分小。于是可不计（x-x0）的三次及更高次幂的各项，则上式简写为：在结点，x=x0-h, 在结点1, x=x0+h，代入(b) 得：)(.)(! 21)(20022000bxxxfxxxfff)(.20222003cxfhxfhff)(.20222001dxfhxfhff59联立(c),(d),解得差分公式： )11(2310hffxf)21 (22031022hfffxf同理，在网线-上可得到差分公式)41 (22)31 (2042022420

43、hfffyfhffyf60从而可导出其它的差分公式如下：)()(461)()(241)()(461)()(41121042040448765432104022411931040447586022fffffhyffffffffffhyxffffffhxfffffhyxf61 相隔2h的两结点处的函数值来表示中间结点处的一阶导数值，可称为中点导数公式。 以相邻三结点处的函数值来表示一个端点处的一阶导数值，可称为端点导数公式。 中点导数公式与端点导数公式相比，精度较高。因为前者反映了结点两边的函数变化，而后者却只反映了结点一边的函数变化。62边界元法简介边界元法简介 边界元法（边界元法（bounda

44、ry element method） 一种结合有限元法和边界积分法发展起来的一种一种结合有限元法和边界积分法发展起来的一种新数值方新数值方法法 只在定义域的只在定义域的边界边界上上划分单元划分单元，用满足控制方程的函数去，用满足控制方程的函数去逼近逼近边界条件边界条件。 适用于应力（薄板）、适用于应力（薄板）、流体力学流体力学、声场、电磁场等问题、声场、电磁场等问题63边界元法基本思想边界元法基本思想 以微分控制方程的基本解为权函数，利用以微分控制方程的基本解为权函数，利用加权余量法加权余量法将区将区域积分转化为边界积分，并结合求解域边界的离散，构建域积分转化为边界积分，并结合求解域边界的离散

45、，构建基于边界单元的代数方程组，然后进行计算求解基于边界单元的代数方程组，然后进行计算求解以定义在边界上的边界积分方程为控制方程，通过对边界元插值离散，以定义在边界上的边界积分方程为控制方程，通过对边界元插值离散，化为代数方程组求解化为代数方程组求解 降低了问题的降低了问题的维数维数，从而显著降低了，从而显著降低了自由度数自由度数边界的离散比区域的离散方便得多，可用较简单的单元准确地模拟边边界的离散比区域的离散方便得多，可用较简单的单元准确地模拟边界形状，最终得到阶数较低的线性代数方程组界形状，最终得到阶数较低的线性代数方程组64加权余量法简介加权余量法简介一种直接从所需求解的一种直接从所需求

46、解的微分方程微分方程及及边界条件边界条件出发，寻求出发，寻求边值问题近似边值问题近似解解的数学方法。的数学方法。从从静力静力发展到了发展到了动力、稳定、材料非线性动力、稳定、材料非线性和和几何非线性几何非线性等各方面。等各方面。在求解域上建立一个在求解域上建立一个试函数试函数 试函数由完备函数集的子集构成。已被采用过的试函数有幂级数、三角级数、试函数由完备函数集的子集构成。已被采用过的试函数有幂级数、三角级数、样条函数、贝赛尔函数、切比雪夫和勒让德多项式等等。样条函数、贝赛尔函数、切比雪夫和勒让德多项式等等。试函数与真实解之间的偏差，即试函数与真实解之间的偏差，即余量余量（内部和边界）（内部和

47、边界）引入权函数，定义消除余量的条件引入权函数，定义消除余量的条件加权余量法就是一种定义近似解与真解之间余量，并设法使其最小的加权余量法就是一种定义近似解与真解之间余量，并设法使其最小的方法。方法。65设问题的控制微分方程为:在V域内 在S边界上 式中 : L、B分别为微分方程和边界条件中的微分算子； f、g 为与未知函数u无关的已知函数域值； u为问题待求的未知函数。( )0B ug( )0L uf66当利用加权余量法求近似解时，首先在求解域上建立一个试函数 ，一般具有如下形式：u 1niiiuC NNC式中： 待定系数，也可称为广义坐标；iC取自完备函数集的线性无关的基函数。iN由于 一

48、般只是待求函数u的近似解，因此代入后将得不到满足，若记：u ( )( )IBRL ufRB ug在V域内在S边界上显然 反映了试函数与真实解之间的偏差，它们分别称做内部和边界余量余量。BIRR、67 若在域V内引入内部权函数 ，在边界S上引入边界权函数则可建立n个消除余量的条件，一般可表示为：IWBW0(1,2, )IiIBiBVSW R dVW R dSin不同的权函数 和 反映了不同的消除余量的准则。从上式可以得到求解待定系数矩阵C的代数方程组。一经解得待定系数，由式(5.1.3)即可得所需求解边值问题的近似解。BiWIiW68IiW 由于试函数 的不同，余量 和 可有如下三种情况，依此加

49、权余量法可分为：1内部法试函数满足边界条件，也即 此时消除余量的条件成为:2边界法试函数满足控制方程，也即 此时消除余量的条件为:( )0BRB ug( )0IRL uf0(1,2, )IiIVW R dVin0(1,2, )BiBSW R dSinu IRBR693混合法 试函数不满足控制方程和边界条件，此时消除余量的条件为: 显然，混合法对于试函数的选取最方便，但在相同精度条件下，工作量最大。对内部法和边界法必须使基函数事先满足一定条件，这对复杂结构分析往往有一定困难，但试函数一经建立，其工作量较小。0(1,2, )IiIBiBVSW R dVW R dSin70边界元法特点边界元法特点 1 前处理工作量小前处理工作量小 2 解算规模小解算规模小 3 求解奇异性问题时计算精度高求解奇异性问题时计算精度高 4 在在载荷集中载荷集中和和半无限域半无限域等问题上有优势等问题上有优势 相对于有限元法，边界元法发展较慢，相对于有限元法，边界元法发展较慢，71有限元法解决应力集中问题 在应力分析中对于应力集中区域必须划分很多的单元，从而增加了求解方程的阶数，计算费用也就随之增加 用位移型有限元法求解出的应力的精度低于位移的精度，对于一个比较复杂的问题必须划分很多单元，相应的数据