有限元法简介.ppt

1.2有限元法简介,1求解偏微分方程的一种数值分析方法上世纪50年代起源于飞机结构分析迅速扩展到固体力学各领域及非结构分析问题2国外成熟的商业软件及开发公司MSC公司ADINA公司等NASTRANANSYSADINA,3在我国的发展上世纪70年代SAP5ADINA4在汽车及发动机领域的应用结构分析传热及热负荷分析冷却润滑及磨损分析,第二章有限元法基础,现代内燃不断地追求更高地经济指标，常常采用强化和燃用重油等措施来降低燃油耗量。这导致了内燃机零部件的机械负荷与热负荷的增加。如何保证内燃机结构的可靠性与寿命便成为一个问题。内燃机的大部分零部件不但结构复杂，而且所受的载荷类型也各不相同。零部件的结构分析和设计，长期以来主要采用实验分析和经验设计的方法。由于测量方法的局限性和某些测试技术的复杂性，加上测量方法的同期较长、费用较高，要想通过大量的不同方案的测试研究来寻找一个最佳的设计，往往是困难的和不合算的。,而且对这些复杂的结构件，如活塞、曲轴、缸盖等，也不可能得到应力分析的解析解。在这种情况下，只有寻求结构的数值分析方法，才能为内燃机零部件的结构分析和设计提供有力的手段。,电子计算机的出现，使内燃机结构分析和设计的多种新理论和方法得以推广应用。有限元法就是其中一种强有力的数值计算方法，在结构形状和所受载荷类型相当复杂的情况下，它都可以进行求解，并能得到较准确的结果,有限元方法进行结构分析的思想是：将一个连续的弹性体进行离散化，分割成彼此用节点相连接的有限个单元，然后对单元进行分析，用节点位移来表示结构的变形，再建立整个结构总位能关于结构位移的表达式；根据变分原理，可以得到以节点位移为示知数的一大型方程，用消元法或迭代法，即可求出各节点处的位移近似值，进一步可求出各节点的应力值。这种先分后合，以有限的单元代替连续的弹性体的方法就是有限元法的基本思想。,、结构的离散化,有用限元法分析变截面直杆（图2）时，先将该杆进行离散化处理，即将其分割成若干个互不重叠的等断面直杆。这里，为讨论问题简便起见，将它会为两段直杆，如图2中的段和段，每一段称为一个单元。单元的端点称为节点。,因此，结构离散化后，单元之间只依靠节点连接。在节点上可以有载荷作用，图2中的p就是作用在节点上的载荷，称此载荷为节点载荷。最后将节点和单元分别编上号码。这样，结构被离散后，分割成由三个节点和两个单元组成的等价系统。问题就变成，只要求出这些节点处的位移值，就可以近似地表示这个变截面杆的变形情况，进行可以求出它的应力分布。因此，节点位移就成为有限元分析中待求的示知参数。,、选择位移插值函数,用节点上的位移值，来建立整个求解区间的位移函数。这种函数称为位移插值函数。它的作用就是将节点之间的位移用适当的方式补充起来。为了在整个区间上构造位移插值函数，可以在每个单元上分别进行。通常，假设位移函数为多项式的形式，其中最简单的位移函数是线性函数。在本例中，由于杆中应变处于常应变状态，故可采用线性位移插值函数。,在本例中，由于杆中应变处于常应变状态，故可采用线性位移插值函数。例如，对于第i个单元，若其两个节点i、i+1上的位移值分别为、（图）那么，位移插值函数的形式为,式中，。显然，这样构造的位移插值函数在相邻两个单元的公共节点上保持连续。,、用变分原理推导出有限元法计算格式,为了得出有限元法计算格式，常采用固体力学的变分原理，如最小总位能原理。最小总位能原理：受到外载荷的弹性体的总位能，为弹性体的变形能和外载荷的位能之和，如果弹性体处于稳定平衡状态，那么总位能为极小值。在结构的有限元法分析中，结构的总位能为各单元的位能之和。对于本例，单元的变形能按下式计算：,（22）,式中，为单元e的截面积；为材料的弹性模量；为应力；为应变。将式（2）所表示的位移插值函数代入式（2），得到,（2）,对于第一个单元，有,对于第二个单元，有,整个变截面直杆的变形能等于各个单元变形能的总和，即,（2）,节点载荷的位能,（2）,式中，为作用在第个节点上的载荷。于是，变截面直杆的总位能,（2）,根据最小位能原理,得到,（2）,把,(2),将方程()写成矩阵的形式：,（2）,再将上式写成如下的形式,式中为结构的总刚度矩阵，,不难看出，是奇异的，方程（2）没有唯一的解。必须进行一些处理。,式（2）中的是节点处的反力，因为节点是被约束住的，因此位移为零。通过划去最后一行和最后一列，得到修正的方程为：,（2-10）,4.求节点未知位移有限元法的最后结果，是要求解形式如（2-9）的线性方程组.求解线性方程（2-10），就可以得到未知的节点位移量,5.计算单元的应变和应力在结构各节点的位移求出后，可进一步求出各单元的应变和应力.两个单元的应变分别为,因此，两个单元的轴向应力分别为,式中的负号说明应力都是压应力.显然，计算结果和解析解完全是一样的,从上面这个简单的例子可以看出，用有限元进行结构分析时，是直接对结构进行离散化处理，在单元分析的基础上，利用变分原理，建立计算格式。由于单元和节点的划分比较任意，用有限元法求解结构复杂的零部件，可以得到较好的结果.,2-2弹性力学的平面问题,内燃机中的一部分结构件，如连杆，主轴承盖，机架的横隔板等，可以近似看作弹性力学的平面问题，采用平面问题的有限元法其解.由于平面问题的有限元法，无论在程序结构的复杂性，计算前后数据处理的工作量方面，还是在计算机时间和费用上，都要比三维问题的有限元法优越的多.本章只着重讨论三节点平面三角形单元的情形，这种单元的分割比较灵活，适宜处理边界比较复杂的零件，对内燃机结构分析较为适用.此外，还介绍一种八节点的四边形单元.,一.弹性力学平面问题的基本方程,弹性力学的平面问题有两种类型，即平面应力问题与平面应变问题.因为两者平面问题都具有相似的基本方程，在以后的叙述中，将主要讨论平面应力问题.,设有一很薄的均匀板，如果板的厚度比其他两个方向（长和宽）的尺寸要小得多，且只在板的厚度断面上受到平行于版面的均匀分布的外力作用，这类问题称为平面应力问题.例如内燃机的主轴承盖,平面应力问题有三个应力分量，即,写成应力矩阵的形式为,（2-11）,1.应变位移方程,弹性体在受力作用后，将产生位移和应变.弹性体内任一点的位移，可用它在坐标轴x,y上的分量u和v来表示，其应变用表示，描述应变与位移之间的几何关系的微分方程，可以用一个矩阵表示：,（1-22）,2.应力应变方程对于连续，均匀，各向同性的弹性体，应力分量与应变分量之间的关系由下式表示：,（2-13）,式中，为材料的泊松系数；G为剪切弹性模数.式（213）还可以写成另一种形式,（2-14）,引入,则式（2-14）可写成,D称为弹性矩阵，它是一个只与弹性体的弹性常数有关的常数矩阵.式（1-15）是描述应力与应变关系的物理过程.,（2-15）,3.初应变,在结构的应力分析中，如果还要考虑温度变化所引起的热应力时，那么，要在上述的总应变中考虑温度变化所引起的初应变.对于各向同性物体，若线性膨胀系数为，温度的增加为T，那么初应变为:,(2-16),式（2-5）所表示的物理方程为,（2-17）,4.平面应变问题,平面应变问题的基本方程，只要在平面应力问题的基本方程中，用代替E,用代替，用（1+）代替，就可以得到.,2.弹性体的离散化,对于平面问题，结构离散化处理时，常采用三角形单元（图1-4）.分割过程的目的主要是决定单元的数目，形状，大小以及排列，使得原来的结构尽可能精确地得到模拟.划分单元时，,图2-4,要注意使任一个三角形单元的顶点同时是其相邻三角形单元的顶点；要尽量使每个三角形单元三边边长不要相差太大；要使节点和分割线安置在集合形状和载荷发生突变处.在内燃机中，大部分零件都有应力变化较大的部位.因此，在可能会出现应力集中或应力梯度比较大的区域，应将单元划分得小些，以便在这些地方能得到较精确的计算结果.,一般来说，结构分割的单元越多，用有限元法所计算的结果就越精确.但是单元和节点的数目要受到计算机存储容量和计算时间及费用的限制.即节点和单元数目越多，所占用的计算机存储量越大，计算时间越长，费用越高.所以在进行结构的有限元法计算时，要根据电子计算机存储量的大小和工程上的计算精度要求，决定划分单元和节点的总数目.,对于已分割好的结构，要将所有的单元和节点按一定顺序编号.单元的编号除根据人为需要约定外，原则上是可以任意的；而节点编号，则要根据一些算法的特殊要求进行，有的需要进行优化，以节省电子计算机存储量.这种编号的原则是，每一节点编号与其相邻的所有节点编号的差值应尽可能地小（图1-5）.,对于大型复杂的结构分析问题，上述人工划分网格的工作量是很大的.这时，可采用自动网格划分程序，在电子计算机上实现单元的划分和节点，单元的编号以及节点坐标的自动生成.通常把这一过程叫做计算机的前处理.,三位移插值函数先考察一编号为e的单元（图2-6所示），设其三个节点按逆时针次序排列的编号为，j,m三个节点的坐标分别为由于是平面问题，每个节点有两个位移分量.三个节点的位移分量分别为,图2-6,X,Y,0,j,m,将三角形单元e的节点位移写成,（2-18）,下面，根据节点位移值来构造位移插值函数.这里采用最简单,的线性插值的方法.设单元e上的位移是坐标x,y的线性函数,（2-19）,式中，为待定系数.将各节点的坐标代入式（1-19），就有,联立求解这两个方程组，可得到,（2-20）,（2-21）,式中，,为三角形单元的面积，,将式（2-20）和（2-21）代入式（2-19）中，得到位移插值函数的表达式：,（2-22）,（2-23）,得到位移函数的简单形式,（2-24）,将式（2-24）写成矩阵的形式,式中，是二阶单位矩阵,（2-25）,此处从第九页开始,式中的是坐标的函数，它反映了单元的位移形态，一般称为形状函数。,根据式（2）所表示的位移函数可以知道，在单元的边界上位移是线性变化的。由于两个相邻的单元在其公共节点上具有相同的节点位移，因此在它们的公共边上，也具有相同的位移。这就说明了上述位移插值函数保证了任意相邻单元之间位移的连续性，即满足协调条件。,四、单元的应变与应力利用上面构造的单元的位移插值函数，可以计算单元的应变和应力。,1.单元的应变将式（2）代入式（2），可得到单元上的应变为,式中，,上式后面的表示这个公式代表三个公式，其中角码轮换代替就可以得到这三个公式。由于单元的面积和系数等都是常数，单元的应变分量也是常数。,（2）,2.单元的初应变单元内的初应变一般来说是不均匀的，但是当单元划分得较小时，可以用一个平均值来表示单元的初应变。若单元的线膨胀系数为a，单元上温度的平均值,（2）,式中，分别为节点的温度，则单元的初应变,（2）,3.单元的应力把式（）代入（），可得到单元应力的表达式为,（2）,令式中,，,称为应力矩阵，则式（2）可写成,当考虑由于温度变化引起的初应变时，单元的应力为：,（2）,=,（）,五、单元刚度矩阵在用固体力学的变分原理（如最小总位能原理）导出有限元计算格式时，需要求出结构总位能的表达式。结构的总位能为各单元总位能之和。因此，首先必须求结构的每单元总位能。,1.单元的变形能设单元具有单位厚度，其变形能为,（）,式中，,（2）,C是与节点位移无关的项。（二）单元的刚度矩阵式(2-34)表示的称为单元的刚度矩阵，它是一个对称矩阵。单元刚度矩阵由下式计算：,(1-35),式中，,为阶子矩阵，它按下式计算,(2-36),(2-37),在上式中,已考虑了单元的厚度t。,（三）单元上外载荷的位能通常，作用在单元上的外载荷有集中力、表面力和体积力。设单元上的三个节点所受的集中力为,在节点发生位移时，集中力的位能,结构上有集中力作用时，一般取集力的作用点为节点，如果集中力未能取在节点上，则需要按静力等效的原则将其移置到节点上。,(2-38),结构常受到边界上的分布载荷的作用，如连杆大、小端的轴承负荷，这种力称表面力。设边界上单位长度的表面力,则在单元发生位移时，表面力的位能,(2-40),式中，为单元的边界；s为边界的弧长。结构上有时还受到重力、惯性力等体积力的作用，设单位体积所受的体积力,在单元发生位移时，体积力的位能,于是，最后可得到单元e的总位能的表达式为,(2-41),六、整体刚度矩阵完成了结构的单元分析后，接着要进行的是整体分析，完成结构的有限元计算格式。,（一）结构的总位能现在，设所研究的结构被分为NE个三角形单元，NP个节点。节点和单元分别统一编号。结构的总位能由下式得到：,(2-44),(2-42),令,(2-43),(2-45),(1-46),(2-47),(2-48),式中，为结构的整体刚度矩阵；为作用在结构上的热负荷向量；为作用在结构上的体积力向量；为作用在结构上的表面力向量；为作用在结构上的集中力向量；F为作用在结构上的总的负荷向量。因此，结构总位能表达式可简写为,(2-49),式中，为结构的节点位移向量，而单元e的三个节点的位移在结构位移向量中的位置为,(2-50),式(1-49)所表示的总位能是节点位移的二次函数。其中的整体刚度矩阵是一个奇异矩阵。总刚度矩阵的奇异性是由于允许有刚体位移而引起的。不难证明单元刚度矩阵具有奇异性，由单元刚度矩阵叠加的总刚度矩阵也具有奇异性。在计算中，必须排除这种刚体运动。通常根据结构的约束条件，可选择结构上的某一点的位移分量为零，再选取一点的某一位移分量(或)亦为零，对刚度矩阵k、位移向量和负荷向量F也进行相应的处理，即所谓约束处理。这样就得到进行约束处理后的和。经处理后的整体刚度矩阵是一个对称正定的矩阵。,根据最小总位能原理，节点位移必须满足：,这里，,表示对位移向量的各个分量,的偏导数,于是可以得到：,这是一个以对称正定矩阵,为系数，以节点位移为未知数的线性代数方程组。求解此方程，就可以得到节点位移向量，进一步就可以求出结构的应力,（二）整体刚度矩阵的形成,将结构的各单元刚度矩阵中的子矩阵,、,、,按式(2-37)计算后，分别叠加到整体刚度矩阵K,中的相应位置上去，就形成了结构的整体刚度矩阵,。下面，举一个简单的例子来说明整体刚度矩阵的形成过程。,有一对角受压的正方形薄板，载荷沿厚度方向均匀分布。由于两对角线是对称轴，所以只取四分之一的部分作为计算对象。将计算部分划分为四个单元，六个节点（图1-7）。,计算出每一个单元的刚度矩阵后，将其扩大为的矩阵。对于第一个单元，单元刚度矩阵的子矩阵分别为、，单元刚度矩阵,对于第二个单元，。单元刚度矩阵的子矩阵分别为,单元刚度矩阵,对于第三个单元，,单元刚度矩阵的子矩为,单元刚度矩阵,对于第四个单元，单元刚度矩阵的子矩阵分别为,单元刚度矩阵,于是，结构的整体刚度矩阵,（三）约束处理如前所述，对于已形成的整体刚度矩阵，必须进行相应的约束处理，使线性代数方程组（2-51）具有唯一的解。在选择位移约束条件时有两种情况：其一是结构本身已具有位移约束条件可以消除刚体位移，如上例；其二是结构原有的位移约束条件不够或没有位移约束条件。对于后者，必须选择足够的位移约束条件来消除刚体位移。处理位移约束常用的方法有两种，下面分别讨论。当已知结构的某些节点位移分量为零时，在线性代数方程组（包括整体刚度矩阵）中，划去和节点零位移相对应的行和列，就可以消除刚体运动。现仍以图2-7为例。从图中可以看出，这里有六个位移分量为零，即,划去K中第1、3、7、8、10、12行及第1、3、7、8、10、12列的元素，并将和F中的第1、3、7、8、10、12行的元素划去，这样就得到一个6阶的线性方程组。,式中，为刚度矩阵中的非零元素。,刚度矩阵式一对称正定矩阵，方程组具有唯一的解。从上述过程可以看出，这种约束处理方法的优点是，最后要求解的方程组的阶数下降了，这样相应的可以节省计算时间。另一种处理已知零位移的方法是，将主题刚度矩阵中，相应于零位移分量的那些行的主对角线元素改为1，其余的元素，连同右端项中的相应元素都改为零。在上例中，将第1、3、7、8、10、12行主对角线元素改为1。其他元素连同右端项的相应元素都改为零，即,这种处理约束方法的优点是，编制程序要比上述划行划列的处理约束方法简单,七、负合向量的计算（一）热负荷向量的计算,由于单元的初应变是常数向量，和都是元素为常数的矩阵，所以单元e上的热负荷向量,（2-52）,将矩阵,（2-53）,按上式将结构的全部单元的热负荷向量叠加起来，就可得到结构的热负荷向量,（二）体积力向量的计算作用在单元e（厚度为t）上的体积力向量,在计算中，单元的体积力可按下式平分到其三个节点上,（2-54）,按上式将结构各单元的体积力向量叠加起来，就可以得到结构的体积力向量。,（三）表面力向量单元的表面力向量,式中，为作用在边界上的均布载荷数密度；为边界单元在边界上的一边的长度，其长度,（2-55）,单元上的表面力向量可按下式平分到边界的两个节点（j,m）上,（2-56）,将全部边界单元上的表面力向量按上式叠加起来，就可以得到结构的表面力向量。,（四）集中力向量如果作用在第p个节点上的集中力分量为,（2-57）,则集中力向量可直接得到,2-3平面问题的等参数单元,前面介绍的常应变三角形单元是最低阶单元。在内燃机工程中，许多构件的应力场是随着坐标的不同而发生急剧的变化的，采用常应变三角形单元时，须将结构分割为大量的单元和节点，才能得到较好的计算精度。本节介绍一种四边形等参数单元，单元位移函数采用二次以上的插值函数，能较好地反映单元内应力变化的情况。采用较少的单元也能获得较精确的计算结果。,一、矩形单元在求解平面问题时，可以采用一种矩形的单元。把一个平面结构的求解域分割成若干矩形单元后，再按照与三角形单元分析相似的过程，可得到矩形单元的计算格式。,图2-8为一矩形单元e其边长分别用2a和2b表示，以平行于两临边的两个对称轴,图2-8,为x轴和y轴。按逆时针顺序排列的四个顶点1、2、3、4的坐标分别为。设矩形单元e的节点位移向量,将矩形单元的位移插值函数取为,（2-58）,（2-59）,在矩形单元的四个节点上，有,求解以上两方程组可得到再将这些系数代入式（2-59），就能得到矩形单元e的位移插值函数：,(2-60),或简写成,(2-61),式中，,从以上可以看出，矩形单元的插值函数是双线性函数。由于其中有二次项，所以单元的应变和应力不再是常量，而是线性变化的。设，则式(2-61)可以写成统一的形式,(2-64),，,(2-63),插值函数可以写成如下形式：,由于插值函数是双线性的，在矩形的每一边上，插值函数分别是x或y的线性函数，其值只取决于这边上两节点的函数值。这样构造的插值函数，在相邻的两矩形单元的公共边上，只要在公共节点上有同样的值，就能保证连续性的要求。这样的双线性插值，对于任何线性函数都是成立的，于是有等式,(2-65),式(2-64)是用节点坐标来表示位置坐标变量的表达式，式中的系数是形状函数。这样，在单元e上，位移插值函数公式与位置坐标的表达式具有相同的形式，即它们具有完全相同的形状函数作系数。这种单元称为等参数单元。矩形单元虽然比三角形单元能更好地反映单元内应力的变化情况，但由于它不能处理具有曲线边界或非正交直线边界的结构，并且不能随意改变大小，所以在实际工程中应用并不多。,二、四边形等参数单元（一）四节点四边形等参数单元,在计算边界不规则的平面问题时，可采用任意四边形单元。如图2-9所示，一任意四边形单元，其四个顶点，的坐标分别为如果仍采用矩形单元的位移插值公式，在相邻两单元的公共边上由于位移不是线性变化的，因此不能保证位移插值函数的连续性。这时，可以采用坐标变换，将任意四边形单元变换到矩形单元上去。,图2-9,在图2-9所示的任意四边形单元上，以分别通过四边形两对边中点的两直线的交点为原点，两直线分别为轴和轴，令四个边上的和的值分别为和，就得到一个、坐标系，称为单元的局部坐标系。而原坐标系为整体坐标系。为了实现坐标变换，引入一个如图2-10所示的四节点正方形单元，并使正方形单元和任意四边形单元之间的节点一一对应。,图2-10,在局部坐标下，矩形单元位移的插值函数为,（2-66）,形状函数,（i=1,2,3,4）,(2-67),上式是单元在局部坐标（）下位移插值函数。在实际计算中，需要用整体坐标下的位移插值函数。所以必须给出整体坐标（x,y）与局部坐标（）之间的转换公式。仿照位移插值函数式（2-65），坐标变换式可以写成,利用上市，就可以将图2-10所示的正方形单元变换为图2-9所示的任意四边形单元。比较式（2-65）和式（2-67），可以看出位移插值函数式和坐标变换公式具有完全相同的形式，即它们用相同的节点值作为参考，并具有完全相同的形状函数。这种任意四边形单元是一种等参数单元，称为四节点四边形等参数单元。,（2-68）,（二）八节点四边形等参数单元,为了改进四节点四边形等参数单元的精度，可增加节点数目，提高位移插值多项式的次数在工程中，常采用八节点四边形等参数单元。图1-11为在局部坐标下的八节点正方形单元，边长为2.除原来四个节点外，又取各边的中点为节点。为了满足协调性条件，位移函数取为,由节点处的函数值，可以决定系数。若节点处的位移值分别为和，可以得到位移插值函数的公式为,根据在节点i为1，其余节点处为零，并考虑到八个节点的局部坐标值.,-,形状函数,有如下表达式,（1-71）,式（1-71）是单元在局部坐标下的位移插值函数公式。需要用下列坐标变换式，将其转换到整体坐标（x,y）上去，,（1,(1-72）,局部坐标下的正方形单元，经坐标变换后，变成整体坐标下的八节点曲边四边形单元。,如图1-12所示，这说明对边界为曲线的结果进行有限元分析时，应用八节点四边形等参数单元可以在边界上得到很好的逼近，从而提高了计算精度。,三）位移-应变方程在求有限元计算格式时，必须求出单元刚度矩阵。这就要计算单元的应变，对于四边形单元，根据位移-应变方程有,（1-73）,上式中，对于四节点四边形单元p=4，对于八节点四边形单元p=8。这时，矩阵,(1-74),而子矩阵,（1-75）,但是，四边形参数单元的形状函数式局部坐标下的函数，它们不能直接对x及y求导。这里，需要采用复合函数求导法则，,即,（1-76）,将上式写成矩阵形式，有,（1-77）,令,（1-78）,则式（1-77）变成,（1-79）,式中，为雅可比矩阵，它按下式计算：,（1-80）,由上式可求出的逆阵，就可以得到,（1-81）,利用式（1-81），就可以将计算转换到局部坐标下进行。这里需要指出，计算式（1-81）的必要条件是能求逆。即单元的每一个内角应小于180度，在某一内角大于180度时，是奇异的。因此任何内角趋近180度，单元应力的计算精度都会受到影响。在划分单元时，应该特别注意这一点。,（四）单元刚度矩阵平面问题的单元刚度矩阵,计算式，用,代替,dx,dy,（1-82）,式中，,是雅可比矩阵,的行列式，由于在局部坐标下，单元是正方形，,积分区间为,，故式（1-82）可写成,上式中的被积函数往往比较复杂，需采用数值积分的方法，通常采用高斯积分方法。,（1-83）,这种方法将在1-5中叙述。根据式（1-74）和（1-83），可以将单元刚度矩阵写成如下形式：,（1-84）,式中，子矩阵,（1-85）,式中，,为高斯求积公式中的积分点数,；,（五）单元的应力矩阵在节点的位移求出后，单元的应力分量按下式计算：,其余的分析和计算与平面三角形单元的类似，此处不再复述。,（1-87）,26平面稳定温度场,一个在内部没有热源的平面物体，其温度分布只是坐标x和y的函数，与坐标z无关，即温度沿z向是不变的。这种只随坐标x和y而改变，但不随时间变化的温度场，称为平面稳定温度场。,一.导热基本方程与边界条件结构处于平面稳定温度场的情况下，其温度分布函数T(x,y)在平面区域内满足调和方程,要惟一地确定结构上的温度分布，必须知道结构边界上的边界条件，边界条件是指结构的表面与周围介质相互作用的规律，常以下列两种方式给出：,2-168,（一）温度边界条件在计算时，如果结构边界上得温度值为已知，这种边界条件称为温度边界条件，即,式中，f为边界上的已知温度函数值。一般通过实测的方法得到边界上若干点的温度值，再进行插值而得到。,（2169）,（二）热交换边界条件,在计算时，如果已知结构表面与周围介质进行热交换的情况，这种边界条件称为热交换边界条件。设结构表面周围介质的温度为，介质与结构表面的热交换系数为，结构的热传导系数为，则在结构表面上得热交换条件为,（2170）,式中，为外法向导数。若令，则上式可写成,当结构的某部分边界条件与周围介质不进行热交换时（即=0），这部分的边界条件上式称为绝热边界条件。如果在结构的一部分边界上给出了温度边界条件，在另一部分边界上给出了热交换边界条件，这种情况称为混合边界条件。,（2171）,二.变分原理求解结构的稳定温度场问题，可以化为一个求泛函极值的变分问题。这种变分原理在数学上已给出了证明。下面，给出结论。,1.满足热传导方程式（2168）和边界条件式（2171）的温度函数T(x,y)，应该使泛函,（2172）,取极小值的函数。如果温度函数T(x,y)使泛函式（2172）取最小值，则T(x,y)在区域内满足热传导方程式（2168），并在边界上满足条件式（2171）。,2.满足热传导方程式（2168）和边界条件式（2169）的温度函数T(x,y)，应该是使泛函,取极小值的函数。如果温度函数T(x,y)使泛函式（2173）取极小值，则T(x,y)在区域内满足热传导方程式（2168），并在边界上满足条件式（2169）。,（1173）,根据上述变分原理，利用求泛函式（2172）或式（2173）极小值的方法，来求解结构的温度函数T(x,y)。具体方法是，先将结构离散化，建立单元的泛函关于温度的表达式。由各单元的泛函叠加而得到整个结构的泛函关于温度的表达式。再由求泛函极值的方法，得到以结构节点温度为未知数的线性代数方程组，解之可以求得结构节点的温度值。这就是用有限元法计算结构稳定温度场的基本思想。下面，以热交换边界条件为例，讨论平面稳定温度场有限元计算方法。,三.结构的离散化,在推导平面稳定温度的有限元计算格式时，先将结构的求解区域离散化，把它分割成若干个单元（这仍旧令其为三角形单元）。划分单元的规律和平面应力分析中叙述的地方相类似。为了统一计算格式，在划分单元时，规定每个边界单元只有一条边落在边界上，并规定与这条边相对的顶点编号为i.对于单元编号，则规定内部单元编号在前，边界单元编号在后。如果结构还要计算热应力和机械应力，则应使结构的有限元分割最好与应力分析时的分割相一致。,四.温度插值函数,设e为结构上任一三角形单元（见图26），其三个节点i,j,m的坐标分为，。三个节点处的温度分别为，用矩阵表示：,与构造位移插值函数一样，可以构造单元e上的温度插值函数。这里，温度插值函数仍取为线性函数，即,由节点的温度可求出上式中的系数。由此而得到温度插值函数为,式中：和与平面应力分析中计算式完全一样。如果令,（2174）,则式（2174）可以写成,（2175）,五。结构的泛函根据单元上的温度插值函数，可以求得单元的泛函。由单元的泛函叠加就可以得到结构的泛函。由式（2174）有,对于边界单元，设为边界段的长度，其值按式（255）计算，并设s为从j点起到上某点的距离，又令，则在边界上段上的温度T也可以用和的线性插值函数来构造（图216），即,（2176），（2177）,如果在边界上，用平均的热交换系数和平均的温度函数来代替和，将式（2176），（2177）和式（255）代入式（2172），可以得到单元e上的泛函,如果单元e是内部单元，即，则上式没有第二项。由式（2178）可以看出，泛函是节点温度的二次函数。将结构上的所有单元的泛函加起来，就可以得到结构的泛函，即,（2178）,（2179）,式中，为二次项，是由式（2178）中叠加而成的，而是由一次项叠加而成。由式（2178）可以看出，只有当温度T(x,y)恒为零时,泛函的二次项才为零，否则必大于零，所以其二次型是正定二次型。由变分原理（i=1,2,NP),（2180）根据式（2179）可以得到,（i=1,2,NP),（2181),或写成,上式是以对称正定矩阵K为系数，以节点温度T为未知数的线性代数方程组。求解这个方程组，就可以得到节点的温度值。其中，K称为温度刚度矩阵。,KT=P。,（2182）,六.单元刚度矩阵,用有限元法计算结构的稳定温度，主要的过程是形成温度刚度矩阵K和右端项P。它们分别是由各单元的温度刚度矩阵和右端项叠加而成的。由式（2178）有,（2183）,由上式易看出，单元温度刚度矩阵为,而右端项为,（2184）,（2185）,对于内部单元，=0.按式（2184）和（2185）将每个单元的温度刚度矩阵和右端项叠加起来，就可以形成结构的温度刚度矩阵和右端项。于是，求解式（2182）而得到结构的温度分布。,七.温度边界条件时的有限元计算在已知结构温度边界条件时，将温度插值函数式（2172）代入式（2170），可得到结构在求解区域上的泛函,而其右端项P中的所有分量全为零。,（2186）,按照上述同样的过程，根据变分原理，可以得到单元的温度刚度矩阵,（2187）,把结构的全部单元的温度刚度矩阵叠加形成总温度刚度矩阵K后，仿照平面应力分析中处理已知位移约束的方法，即将线性代数方程中对应于已知温度节点处的行，列划去；或者将温度刚度矩阵中相应于已知温度节点编号的行和列除主对角元素改为1以外，其余元素均改为零，右端项则令其相应于已知温度节点编号处为已知温度值，其余为，这里，为编号在前的内部单元的个数。求解经过处理后的线性代数方程组，就可以得到节点的温度。,于是,单元e的热负荷向量,将结构的各单元热负荷向量叠加起来,就可以得到结构的热负荷向量L,(二)体积力向量的计算,在轴对称问题中,作用在节点上的载荷是指作用在节点的一个圆周上的载荷,设单位体积力,则作用在单元上的体积力向量,当体积力是常量是时,被积函数中的r和z可取单元形心处的值,这时,单元的体积力向量,这相当于把单元e的体积力平均分配到三个节点上去.把各单元的体积力叠加起来,就可以得到机构的体积力向量,(三)表面力向量的计算同平面问题相类似,规定边界单元只有一边(如边)在边界上,边界上受的表面力,上式相当于把单元边界上的表面力平均分到其两个节点j和m上,将所有的边界单元的表面力向量叠加起来,就可以得到机构的表面力向量Q,并取,则单元的表面力向量,(四)集中力向量若节点p所受的集中力,这里,集中力是指作用在节点所在的圆周上的力的合力,则集中力向量,七.单元应力分量的计算在形成机构的整体刚度矩阵K和负荷向量F后,可求解以节点位移为未知数的线性代数方程组,从而得到各节点的位移值,进一步可以计算各单元的应力分量.将等的表达式代入式(1-104)可以得到单元的各个应力分量:,则单元的应力分量可写成,在计算单元应力分量时,为了便于计算,常用单元形心的坐标值,来代替r和z式(1-125)可以用来计算单元的热应力,机械应力和综合应力,当时,按式(1-116)形成线性方程右端项,求出节点位移后,再按式(1-125)计算,就得到单元的热应力;当不考虑热负荷向量时,令,按式(1-125)计算出的应力就是单元的机械应力,所以,结构的这三种应力计算很容易用一个程序实现.,26平面稳定温度场,一个在内部没有热源的平面物体，其温度分布只是坐标x和y的函数，与坐标z无关，即温度沿z向是不变的。这种只随坐标x和y而改变，但不随时间变化的温度场，称为平面稳定温度场。,一.导热基本方程与边界条件结构处于平面稳定温度场的情况下，其温度分布函数T(x,y)在平面区域内满足调和方程,要惟一地确定结构上的温度分布，必须知道结构边界上的边界条件，边界条件是指结构的表面与周围介质相互作用的规律，常以下列两种方式给出：,2-168,（一）温度边界条件在计算时，如果结构边界上得温度值为已知，这种边界条件称为温度边界条件，即,式中，f为边界上的已知温度函数值。一般通过实测的方法得到边界上若干点的温度值，再进行插值而得到。,（2169）,（二）热交换边界条件,在计算时，如果已知结构表面与周围介质进行热交换的情况，这种边界条件称为热交换边界条件。设结构表面周围介质的温度为，介质与结构表面的热交换系数为，结构的热传导系数为，则在结构表面上得热交换条件为,（2170）,式中，为外法向导数。若令，则上式可写成,当结构的某部分边界条件与周围介质不进行热交换时（即=0），这部分的边界条件上式称为绝热边界条件。如果在结构的一部分边界上给出了温度边界条件，在另一部分边界上给出了热交换边界条件，这种情况称为混合边界条件。,（2171）,二.变分原理求解结构的稳定温度场问题，可以化为一个求泛函极值的变分问题。这种变分原理在数学上已给出了证明。下面，给出结论。,1.满足热传导方程式（2168）和边界条件式（2171）的温度函数T(x,y)，应该使泛函,（2172）,取极小值的函数。如果温度函数T(x,y)使泛函式（2172）取最小值，则T(x,y)在区域内满足热传导方程式（2168），并在边界上满足条件式（2171）。,2.满足热传导方程式（2168）和边界条件式（2169）的温度函数T(x,y)，应该是使泛函,取极小值的函数。如果温度函数T(x,y)使泛函式（2173）取极小值，则T(x,y)在区域内满足热传导方程式（2168），并在边界上满足条件式（2169）。,（1173）,根据上述变分原理，利用求泛函式（2172）或式（2173）极小值的方法，来求解结构的温度函数T(x,y)。具体方法是，先将结构离散化，建立单元的泛函关于温度的表达式。由各单元的泛函叠加而得到整个结构的泛函关于温度的表达式。再由求泛函极值的方法，得到以结构节点温度为未知数的线性代数方程组，解之可以求得结构节点的温度值。这就是用有限元法计算结构稳定温度场的基本思想。下面，以热交换边界条件为例，讨论平面稳定温度场有限元计算方法。,三.结构的离散化,在推导平面稳定温度的有限元计算格式时，先将结构的求解区域离散化，把它分割成若干个单元（这仍旧令其为三角形单元）。划分单元的规律和平面应力分析中叙述的地方相类似。为了统一计算格式，在划分单元时，规定每个边界单元只有一条边落在边界上，并规定与这条边相对的顶点编号为i.对于单元编号，则规定内部单元编号在前，边界单元编号在后。如果结构还要计算热应力和机械应力，则应使结构的有限元分割最好与应力分析时的分割相一致。,四.温度插值函数,设e为结构上任一三角形单元（见图26），其三个节点i,j,m的坐标分为，。三个节点处的温度分别为，用矩阵表示：,与构造位移插值函数一样，可以构造单元e上的温度插值函数。这里，温度插值函数仍取为线性函数，即,由节点的温度可求出上式中的系数。由此而得到温度插值函数为,式中：和与平面应力分析中计算式完全一样。如果令,（2174）,则式（2174）可以写成,（2175）,五。结构的泛函根据单元上的温度插值函数，可以求得单元的泛函。由单元的泛函叠加就可以得到结构的泛函。由式（2174）有,对于边界单元，设为边界段的长度，其值按式（255）计算，并设s为从j点起到上某点的距离，又令，则在边界上段上的温度T也可以用和的线性插值函数来构造（图216），即,（2176），（2177）,如果在边界上，用平均的热交换系数和平均的温度函数来代替和，将式（2176），（2177）和式（255）代入式（2172），可以得到单元e上的泛函,如果单元e是内部单元，即，则上式没有第二项。由式（2178）可以看出，泛函是节点温度的二次函数。将结构上的所有单元的泛函加起来，就可以得到结构的泛函，即,（2178）,（2179）,式中，为二次项，是由式（2178）中叠加而成的，而是由一次项叠加而成。由式（2178）可以看出，只有当温度T(x,y)恒为零时,泛函的二次项才为零，否则必大于零，所以其二次型是正定二次型。由变分原理（i=1,2,NP),（2180）根据式（2179）可以得到,（i=1,2,NP),（2181),或写成,上式是以对称正定矩阵K为系数，以节点温度T为未知数的线性代数方程组。求解这个方程组，就可以得到节点的温度值。其中，K称为温度刚度矩阵。,KT=P。,（2182）,六.单元刚度矩阵,用有限元法计算结构的稳定温度，主要的过程是形成温度刚度矩阵K和右端项P。它们分别是由各单元的温度刚度矩阵和右端项叠加而成的。由式（2178）有,（2183）,由上式易看出，单元温度刚度矩阵为,而右端项为,（2184）,（2185）,对于内部单元，=0.按式（2184）和（2185）将每个单元的温度刚度矩阵和右端项叠加起来，就可以形成结构的温度刚度矩阵和右端项。于是，求解式（2182）而得到结构的温度分布。,七.温度边界条件时的有限元计算在已知结构温度边界条件时，将温度插值函数式（2172）代入式（2170），可得到结构在求解区域上的泛函,而其右端项P中的所有分量全为零。,（2186）,按照上述同样的过程，根据变分原理，可以得到单元的温度刚度矩阵,（2187）,把结构的全部单元的温度刚度矩阵叠加形成总温度刚度矩阵K后，仿照平面应力分析中处理已知位移约束的方法，即将线性代数方程中对应于已知温度节点处的行，列划去；或者将温度刚度矩阵中相应于已知温度节点编号的行和列除主对角元素改为1以外，其余元素均改为零，右端项则令其相应于已知温度节点编号处为已知温度值，其余为，这里，为编号在前的内部单元的个数。求解经过处理后的线性代数方程组，就可以得到节点的温度。,1-7轴对称稳定温度场,对于一个轴对称结构，如果作用在它上面的热负荷也是轴对称的，那么这个结构上的温度分布在圆柱坐标系中仅与坐标z、r有关，而与角度无关。这种结构的稳定温度称为轴对称稳定温度场。例如内燃机组合式活塞的活塞头、气阀和缸套的稳定温度场就是轴对称稳定温度场。,一、导热基本方程与边界条件,在计算结构的轴对称稳定温度场时，取结构的子午面作为求解区域。设z为对称轴，在结构内部无热源的情况下，温度函数在求解区域内满足导热微分方程,（2-188）,要确定结构体内的温度值，还必须给出结构表面的边界条件。与平面稳定温度场的情况相类似，通常也有两种形式的边界条件：,1、温度边界条件,以上各式的符号意义均与平面稳定温度场时的相同。,（2-189）,2、热交换边界条件,（2-190）,因,因,，故式（1-190）可写成,（2-192）,二、变分原理这里，仍以热交换边界条件为例，讨论轴对称稳定温度场的有限元法计算。由变分原理可以知道，满足导热方程式（2-188）和边界条件式（2-192）的温度函数，应该是使泛函,取极小值的函数。,（2-193）,三、温度插值函数用有限元法求解结构的轴对称稳定温度时，首先要将结构进行离散化，将它分割成若干个单元，与轴对称应力分析时一样，这里仍采用三角形单元，划分单元时的规定仍与平面稳定温度场时的相同。设单元的温度插值函数仍为线性函数,根据节点的温度，所得到的温度插值函数的算式为,（2-194）,（2-195）,式中，,（2-196）,四、单元刚度矩阵,如果所考虑的是边界单元，则温度插值函数T和r值在边界上是线性变化的，,由结构的泛函求极值而得到有限元计算格式时，首先仍计算单元的泛函。由式（2-195）有,(2-197),(2-198),将式（2-197）式（2-200）代入式（2-193），得到单元e上的泛函值为,（2-199）,（2-200）,（2-201）,式中，,（2-202）,从而有,（2-203）,由上式可以看出，单元的刚度矩阵为,（2-204）,（2-205）,而右端项为,如果是内部单元，计算时，令，将每个单元的温度刚度矩阵和右端项叠加起来，就可以形成结构的温度刚度矩阵和右端项，再求解线性代数方程，就可以求得结构的节点温度。在式（2-205）、（2-206）中，如果令，则可以得到平面稳定温度场的有限元法计算式。,（2-206）,2-8平面不稳定温度场,内燃机在变工况运行过程（如起动过程）中，常使其一些零件产生较大的温度梯度，而导致热应力过大或超过允许值，结果使这些零件受到破坏。因此，在这种情况下，需要计算结构的不稳定温度场。,一、热传导基本方程在结构内部无热源的情况下，对于各向同性的材料，平面热传导方程可写成如下的形式：,式中，为热传导系数；为材料的比热；为材料的比重；为时间。要唯一地确定式（2-207）的解，必须给出初始条件和边界条件。,（2-207）,1、初始条件,初始条件是指在初始瞬时，结构内部的温度分布规律，即当时,（2-208）,2、边界条件边界条件有两种形式。第一种是：结构表面的温度是时间的已知函数，即在边界s上有,。,（2-209）,第二种是：结构表面与周围介质接触时，通过结构表面的热流密度与结构表面温度和介质温度之差成正比，即,因此，对于不稳定温度场问题，温度场必须满足热传导方程式（2-207）和初始条件式（2-208）以及边界条件式（2-209）（2-210）。根据变分原理，这个问题可以化为求如下泛函的极小值问题：,（2-210）,