有限元法讲义

有限元分析与矩阵根底目录TOC\o"1-3"\u第一局部、矩阵与线性代数3第一章、矩阵的根本概念31.1 矩阵入门31.2 特殊矩阵41.3 矩阵相等、矩阵的加法和矩阵与数的乘法61.4 矩阵的乘法71.5 逆矩阵与伴随矩阵91.6 矩阵的分块101.7矩阵的迹和矩阵的行列式10第二章、矩阵和向量空间122.1 引言122.2 向量空间、子空间和矩阵的张成122.3 线性变换的矩阵表示172.4 基的改变202.5 向量模和矩阵模23第二局部、有限单元法25第三章、有限元法公式的建立253.1 引言253.2 利用虚位移原理建立有限元法的公式293.3 建立有限元步骤归纳36第四章、杆件结构的有限单元法384.1 平面刚架结构384.2 平面桁架结构41第五章、平面问题的有限单元法435.1 前言435.2 连续体的离散化与三角形常应变单元435.3 面积坐标495.4 采用高次位移函数的三角形单元515.5 矩形单元57第六章、等参单元616.1 等参单元616.2 等参杆单元616.3 平面线性等参单元636.4 高斯积分法666.5 8节点平面等参单元67第七章、空间问题有限单元法707.1 前言707.2 四面体单元717.3 正六面体单元777.4 等参单元80第八章、板弯曲问题有限单元法838.1 前言838.2 薄板弯曲838.3 矩形薄板单元878.4 三角形薄板单元918.5 中厚板弯曲928.6 等参中厚板单元93第九章、壳体有限单元法979.1 前言979.2 坐标系统989.3 单元的几何描述1009.4 位移模式1019.5 应力与应变1029.6 单元刚度矩阵与数值积分103第三局部、有限元分析的方程组解法104第一〇章、静力分析中线性方程组的解法10410.1 前言10410.2 基于高斯消去法的直接解法10510.3 应用正交矩阵的直接解法11210.4 高斯—塞德尔迭代法119第一一章、动力分析中平衡方程组的解法12111.1 前言12111.2 直接积分法12111.3 振型叠加法128矩阵与线性代数矩阵的根本概念矩阵入门通过研究线性联立方程组的解法，就容易体会到在实际计算中使用矩阵的有效性。例如〔1.1〕式中x是未知量，利用矩阵的表示法，方程组〔1.1〕可写成为〔1.2〕然而，用矩阵符号表示式(1.2)中的各个数组，方程组就可写成为〔1.3〕式中A是线性方程组的系数矩阵，x是未知量矩阵，b是量矩阵；即，，〔1.4〕现在显然可以对矩阵作出如下的正式定义。定义矩阵是由有序数组成的一个数组。一般的矩阵是由个数排成m行和n列的数组组成：〔1.5〕我们说这个矩阵为阶。当A矩阵只有一行(m＝1)或一列(n＝1)时，可称A为一个向量。在本书里矩阵用黑体字母表示，当矩阵不是向量时，那么用大写字母表示，而向量可以用大写或小写黑体字母表示。根据定义，以下的都是矩阵〔1.6〕矩阵A的第i行第j列的元素表示为,例如在式(1.4)中的第一个矩阵里，，。应当注意，式(1.5)中的元素的下角标变化范围是i从1到m，j从1到n。特殊矩阵当矩阵元素服从某种规律时，我们可以把它作为一个特殊矩阵来考虑。假设矩阵的元素都是实数，那么称为实矩阵；假设其中有复数的元素，那么称为复矩阵。下面我们只讨论实矩阵。此外，矩阵住往是对称的，它有如下所定义的性质。定义将矩阵A的行和列交换所得的矩阵称为A的转置矩阵，记作。假设，可知A的行数和列数相等，且.因为，故将A称为n阶方阵；又因,故称A为对称矩阵。应当指出，对称性即意昧着A是方阵，但反过来就不一定成立，即方阵不一定就是对称的。另一个特殊矩阵是单位矩阵，它是n阶方阵，除对角线上的元素等于1外，其余元素均等于零。例如三阶单位矩阵为〔1.7〕在实际计算中，单位矩阵的阶常常是非常明确的，并且不用写出它的下角标。类似于单位矩阵，我们也要用到n阶单位向量，记为，其中下角标i表示这个向量即是单位矩阵的第i列。我们将更多地和对称带状矩阵打交道。带状意味着在矩阵带宽以外的所有元素为零。因为A是对称的，所以可把带状矩阵的条件规定为〔1.8〕这里是A的带宽，例如下面的矩阵是一个五阶对称带状矩阵，半带宽为2。〔1.9〕假设矩阵的半带宽是零，那么只在对角线上有非零元素，称为对角矩阵。例如，单位矩阵就是一个对角矩阵。图1-1矩阵的一维存储在计算机上进行的计算时。需要用到一个将矩阵元素存人内存中的存贮方案。存储矩阵A的元素的一个明显的方案，是在FORTRAN程序中规定一个给出维数的数组，其中M=m，N=n,而每一个元素存储在的元素之中，但在许多计算中，按这样的存储方案，就会多余地存储很多根本不参加运算的零元素。如果A是对称的。那么应当利用对称性的特点只存储矩阵的上半局部元素(包括对角线上的元素)。通常在内存中可用的存储元素是有限的，为了能把尽可能大的矩阵存入内存，必须使用有效的存储方案。假设矩阵太大以致内存不能容纳时，那么求解过程需要从外存上读出和写入，这样做会显著地增加机器解题的费用。幸而，在有限元分析中，系统矩阵是对称带状的，因此，采用一种有效的存储方案，就可以把比拟高阶的距阵存人内存里。我们用A(I)表示一维存储的数组A中的第i个元素样，一个n阶对角矩阵可简单地按图1-1(a)存储。〔1.10〕假设有如图1-1(b)所示的带状矩阵，以后我们将会看到，在矩阵的带状边界线内的零元素在求解过程中会变为非零元素。例如可能是零元素，但在求解过程中会变为非零元素。因此，我们必须对带状边界线内的零元素分配存储元素，但不必存储边界线以外的零元素。图1-1示出了在有限元求解过程中将要使用的存储方案，我们将在以后的章节中作进一步说明。矩阵相等、矩阵的加法和矩阵与数的乘法现在我们对矩阵相等、矩阵的加法和矩阵与数的乘法作如下的定义。定义矩阵A和B是相等的。当且仅当(1)A和B有相同的行数和列数。(2)所有对应的元素相等，即对所有的i和j有定义当且仅当矩阵A和矩阵B具有相同的行数和列数，矩阵A和B才可以相加。矩阵的相加就是将所有对应的元素分别相加。即假设和表示A和B的一般元素，那么表示C的一般元素。此时。由此可知。C与A和B具有相同的行数和列数。矩阵减法可用类似的方法来定义。根据矩阵减法的定义，可知一个矩阵本身相减得到一个只有零元素的矩阵，这样的矩阵定义为零矩阵O。现在来定义矩阵与数的乘法。定义一个数与一个矩阵相乘。是把该数与矩阵的每一个元素相乘，即意味着.应该注意，到目前为止的所有定义，完全类似于计算普通数所用的定义，并且，两个阶的一般矩阵相加(或相减)需要次加法(或减法)运算。而一个数乘一个阶的一般矩阵需要次乘法运算。因此，当矩阵是特殊形式时(例如对称带状)，我们就可利用这些特点，只计算矩阵C在带状边界线以下的元素，因为其余的元素都是零。矩阵的乘法现考虑两个矩阵A和B，要求出矩阵的乘积C＝AB。定义当且仅当矩阵A的列数等于矩阵B的行数时，两个矩阵A和B才可以相乘得到C＝AB。假设A是阶。而B是阶、那么C中的每一个元素均按下式计算：(1。11)式中C是阶；即式(1.11)中i和j分别从1到p和从1到q变化。因此，计算C的第个元素，是把A的第i行元素和B的第j列元素对应相乘后相加。由于取A的每一行相B的每一列的乘积，可知C一定是阶的。例1.1计算矩阵的乘积，其中可得=(5)(1)十(3)(2)十(1)(3)＝14＝(4)(1)十(6)(2)十(2)(3)=22=(10)(1)十(3)(2)十(4)(3)=28等等因此得到容易证明，矩阵乘法需要进行次乘法运算。但是，在实际处理矩阵时，常常可以利用矩阵内的零元素而减少运算的次数。众所周知，普通数的乘法是可交换的，即。我们要研究上述结论对距阵乘法是否同样成立。考虑以下矩阵〔1.12〕可得〔1.13〕可见，乘积AB和BA是不相同的，因此，矩阵的乘法是不可交换的。实际上乘积与A和B的阶数有关。两个乘积矩阵AB和BA的阶可以是不同的，乘积AB可以有定义，但乘积BA却可能是不能计算的。为了区别矩阵相乘的次序。对于乘积AB，我们说是对B左乘以A或对A右乘以B。虽然通常，但对于特殊的A和B，可能有，这时，我们说A和B是可交换的。虽然在矩阵乘法中交换律不成立，但分配律和结合律都是成立的。分配律表示为(1.14)分配律可以用式〔1.11〕来证明〔1.15〕可得〔1.16〕结合律可以表示为〔1.17〕由于结合律成立，实际上对于一系列矩阵的相乘可以取任一种顺序来进行，如果顺序选择得好。往往可以节省许屡次运算。当计算一系列矩阵相乘时，必须记住的是括号可以去除或参加，幂次可以合并，促必须保持相乘的次序。除了矩阵乘法中的交换律一般不成立外，在矩阵方程中一般也不能像普通数相消那样来进行矩阵相消。特别是，如果AB＝CB，不一定能得出A＝C，这很容易由下面的特殊情况来证明：〔1.18〕但〔1.19〕然而，必须注意到，假设方程AB＝CB时所有可能的B均成立。那么有A＝C。此时我们简单地选取B为单位矩阵I，使得出A＝C还应看到，假设AB＝0，并不能由此得出A或B是零矩阵的结论。下面的特殊情况说明了这一点。〔1.20〕必须指出，在矩阵乘法中有关使用转置矩阵的一些特殊规那么。应当注意的是,两个矩阵A和B的乘积的转置等于按相反顺序的转置矩阵的乘积，即(1.21)逆矩阵与伴随矩阵定义矩阵A的逆矩阵用表尔。假设A的逆矩阵存在，那么的元素满足和。假设一个矩阵有逆矩阵，那么称该矩阵为非奇异的；假设没有逆矩阵，那么称该矩阵为奇异矩阵。计算乘积AB的逆矩阵，可按下述方法进行。设，式中A和B都是方阵。于是〔1.22〕两边均右乘以和，可得〔1.23〕(1.24)因而(1.25)读者当然会注意到,矩阵反置的同一规那么曾在计算矩阵乘积的转置时用过。矩阵的逆可以由，其中为矩阵A的伴随矩阵，是矩阵A的行列式值。定义记n阶方阵A的第i行第j列元素的代数余子式为，那么由元素组成的矩阵的转置，称为矩阵A的伴随矩阵，记为，即矩阵的分块把一个矩阵划分为假设干个子矩阵,可有效地简化矩阵的运算和利用矩阵的特殊形式。子矩阵只是从原矩阵中取出某些行和列的元素所组成的矩阵，这个概念可用下面的一个具体例子来说明。其中虚线是块的划分线。(1.26)应该注意，块的划分线必须穿过整个原矩阵。利用分块，矩阵A可写为(1.27)其中，等等。矩阵的分块对节省计算机存储量是有利的，即如果子矩阵是要重复使用的，就只需要将它存储一次。上述方法同样也可以用在算术运算中，利用子矩阵可以表示一个屡次重复的典型运算。因此，只需进行一次这样的运算，每当需要它时就可以直接调用其结果。用于分块矩阵计算的规那么，是根据矩阵加法、减法和乘法的定义得出的。假设原矩阵已分成可以进行各个子矩阵块相加、相减或相乘的形式，那么，只要把子矩阵看作普通的矩阵元素来进行加法、减法或乘法运算即可。只要我们记住矩阵的分块仅仅是简化矩阵的运算而一点也不改变计算结果的一种手段，那么上述规那么是容易证明并容易记忆的。矩阵的迹和矩阵的行列式仅当矩阵是方阵时，矩阵的迹和行列式才有定义。迹和行列式都是一个数，其值由矩阵的元素来计算，因而它们都是矩阵元素的函数。定义矩阵A的迹以来表示,它等于,其中n是A的阶。矩阵A的行列式可由A的子矩阵的行列式来定义。注意一阶矩阵的行列式就是矩阵的元素本身，即假设，那么。定义阶矩阵A的行列式以表示，它由下面的递推关系来定义：〔1.28〕式中是由删除A的第l行和第j列后的矩阵。许多定理与行列式的应用有关。有代表性的是，通过计算一系列行列式可以求得一个方程组的解。但从现代观点看，利用行列式所得到的许多结果，用矩阵理论可以更有效地求得。例如，用行列式解方程组不是很有效的。正如以后我们会看到的那样，应用行列式的主要作用，在于讨论像距阵的逆存在那样的一些问题时能用简便的表示法。特别是在求解特征值问题时，我们将要用到行列式。计算矩阵的行列式的较为有效的方法是，先把矩阵分解为假设干个矩阵的乘积，然后利用下面的结果：〔1.29〕关系式(1.29)说明，一系列矩阵的乘积的行列式等于每个矩阵行列式的乘积。这个结论的证明是比拟冗长的(利用(1.28)中行列式的定义来证明)，因此在这里我们就不加讨论。在特征值解法中，当需要求矩阵(例如矩阵A)的行列式时，我们将经常用(1.29)式的结果。而所用的特殊分解是，其中L是单位下三角矩阵，D是对角矩阵。此时，〔1.30〕因为，那么有〔1.31〕矩阵的行列式和迹都是矩阵元素的函数。但重要的是，非对角元素不影响矩阵的迹，而行列式是矩阵全部元素的函数。虽然我们能得出这样的结论，即行列式或迹的值大意味着某些矩阵元素的值也大。但并不能得出这样的结论，即行列式或迹的值小就意味着所有的矩阵元素的值也小。矩阵和向量空间引言在第一章中，我们对矩阵作了简单的介绍。矩阵作为由有次序排列的数组成的数组，服从加法、乘法等特定的规那么。在以后各章中会广泛地使用根本的代数规那么，完全熟识它们是非常重要的。然而，至此我们没有讨论过如何从一个实际的物理问题得出矩阵的元素和在什么条件下得到这些元素；如果需要的话，就要用到己介绍过的矩阵代数规那么。换句话说，我们不打算让读者对矩阵在计算中表示的是什么有较深入的了解。然而，如果要很好的理解以后所介绍的数值方法，那么对矩阵有这些了解是非常重要的(确实也是必要的)。本章的任务是介绍一些重要的线性代数概念，懂得这些概念对于完全理解在有限元分析中所用的数值方法是必要的。然而，除了扩大矩阵和线性代数的重要知识外，还希望通过本章介绍的概念，除去那些把矩阵单纯作为一堆数来使用而显得枯燥的一面。即我们只有增加对矩阵应用的了解，才能充分评价求解的方法。向量空间、子空间和矩阵的张成在第一章我们把一个n阶向量定义为用矩阵形式写出的n个数的数组。现在图2-1向量的几何表示我们要把向量的元素与其几何解释联系起来。讨论一个三阶向量的例子，如(2.1)从初等几何学可知，在三维空间所选定的坐标系下，x表示一个几何向量。图2-1表示在该坐标系下所采用的坐标轴及对应于(2.1)式的向量。应该指出,在(2.1)式中x的几何表示完全依赖于所选取的坐标系,即如果(2.1)式给出的是与此不同的坐标系下向量的分量，那么x的几何表示与图2.1所表示的应当不同。因此,仅仅是坐标(或向量的分量)不能定义实际的几何量，而必须同时结出量度它们的具体坐标系。三维几何向量的概念可以推广到有限的n阶向量。假设n＞3，我们再不可能得到向量图。然而，我们将会看到，在数学上所有关于向量的概念都与n无关。如上所述，当我们考虑n＝3的特殊情形时，n阶向量表示一个n维空间的具体坐标系中的一个几何量。假设要处理假设干个n阶向量，这些向量是在一个确定的坐标系中被定义的。在以后各章中将会极其有用的一些根本概念，可由下面的一些定义和论据加以概述。定义如果存在不全为零的数使〔2.2〕那么称向量组是线性相关的。如果向量组不是线性相关的，那么称它们为线性无关向量组。我们用下面的例子来说明该定义的合义。例2.1设n＝3，确定向量组是线性相关还是线性无关。根据线性相关的定义，我们需要检查是否存在满足下述方程的三个不全为零的常数：(a)而从方程组〔a〕得只有当时才满足方程组(a)，因此向量组是线性无关的。考虑这个问题的另一种方法可能更具有吸引力，即如果一个向量组中的任何一个向量都可以用其余的向量来表示时，那么称该向量组是线性相关的。这就是说，假设在(2.2)式中不全为零，例如，那么我们可写成〔2.3〕在几何上，当n＜3时，我们就能将向量组画出，并且如果它们是线性相关的，那么它们当中任一个向量可以用其余向量的组合描绘出来。例如，画出例2.1的向量组，我们立即就看到它们当中没有一个向量能用其余向量的组合来表示，因此，该向量组是线性无关的。假定有q个线性相关构n阶向量，但我们只研究它们当中的任意个向量时，这个向量可能仍是线性相关的。然而，继续减少所研究的向量组的向量个数，我们会得到p个线性无关向量，通常。其它个向量能用这p个向量来表示。于是就得出下面的定义定义假设有p个n阶线性无关向量，其中，这p个向量构成p维向量空间的一组基。因为在该空间中任何一个向量都能用p个基向量的线性组合来表示，所以我们只讨论一个p维向量空间。应该指出，对于所考虑的特定空间的基向量不是唯一的。基向量的一些线性组合可能给出同一个空间的另一组基。特别是，如果，那么所考虑的空间的一组基是，(称作自然基)。由此也得出p不能大于n的结论。定义假设q个向量中有p个向量是线性无关的，就说这q个向量张成一个p维向量空间。因此，我们认识到，基向量组是非常重要的，因为它们是张或所考虑空间的最少向量个数。所有q个向量都能用基向量组来表示，而且q可以是很大(甚至q会大于n)的。假设一个p维向量空间，记为，选为基向量，。然后，我们同样只考虑所有能用、表示的那些向量,而向量和也组成一个向量空间的基，称该向量空间为。假设p=2，那么与会重合。我们称是的子空间，其简明的含义如下面定义所述。定义一个向量空间的子空间也是一个向量空间。子空间上任一个向量也是原来空间的一个向量。如果是原来空间的基向量，这些基向量中的任何一个子集组成一个子空间的基；该子空间的维数等于所选择的基向量的个数。研究了向量空间的一些概念后，现在我们可以认为，任何矩阵A的列也张成一个向量空间，我们称这个空间为A的列空间。同样，矩阵的行张成一个称为A的行空间的向量空间。相反，我们可以将任意q个n阶向量集合成为一个阶矩阵,所用的线性无关向量的个数等于A的列空间的维数。假设矩阵A，计算A的列空间的维数，即要计算在A中有多少个列向量是线性无关的。取A的列向量的任意线性组合，既不会增加也不会减少在A中的线性无关列向量的个数。因此，为了确定A的列空间，我们试图用对A的列向量进行线性组合得到单位向量的方法来变换矩阵A。例2.2讨论向量组成的矩阵A的列空间的维数。把矩阵A的第二列和第三列相应地替换为其第—列和第二列，而用原第一列替换第三列，可得从第三列中威去第一列，再把两倍的第二列加上第三列，最后把第二列乘以(-1)，可得至此，巳将矩阵化为可使我们识别出这三列是线性无关的形式，就是说，因为这些向量的头三个元素是三阶单位矩阵的列，所以它们是线性无关的。而且，由于我们通过A的列交换及其线性组合的过程中得到，并在解算过程中没有因此而增加由矩阵的列所张成的空间的维数，故得到A的列空间的维数是3。在上面的讨论中，为了判断A的列是否线性无关，我们对A的列向量，，……，进行线性组合。另一方面，为了确定由向量组,,..,张成的空间的维数,我们可以利用在式(2.2)中向量线性无关的定义，并考虑以下齐次方程组(2.4)其矩阵形式为(2.5)式中a是一个向量。其元素是,,..,，而A的列是向量,,..,。对A的任意行进行线性组合或乘以因子都不会改变未知量,,..,的解。因此，我们可以通过对A的的各行乘以因子并进行组合，使A简化为只含有单位列向量的矩阵。这个简化了的矩阵称为A的行-梯阵形式。在A的行-梯阵中，单位列向量的个数等于A的列空间的维数。而且可以看出，它也等于A的行空间的维数。由此可见，A的列空间的维数等于A的行空间的维数，即在A中线性无关的列数等于线性无关的行数。这个结论概括在下面关于矩阵的秩的定义中。定义距阵A的秩等于A的列空间的维数。也等于A的行空间的维数。例2.3考虑下面的三个向量：将这些向量作为矩阵A的列，并将矩阵化为行－梯阵形式。有为了把第一列化为单位向量，从第二行起逐行减去乘以不同倍数的第一行，可得为了把第二列化为单位向量，用〔-5〕除第二行，然后从其余各行减去不同倍数的第二行，可得因此我们可以得出以下相应结论：〔1〕的解是，〔2〕向量,,是线性相关的，它们组成一个二维向量空间。向量,是线性无关的，它们组成,,所在的二维空间的一组基。〔3〕A的秩是2.〔4〕A的列空间的维数是2。〔5〕A的行空间的维数是2。线性变换的矩阵表示工程分析中，我们一般需要把一些量与另一些量，例如应力与应变，应变与位移，力和位移，温度与热流等等联系起来。每一种这样的关系都可以认为是从一种性质的量到另一种性质的量的变换，并可简要地定义和用矩阵表示为一种紧凑的形式。我们暂不列举所要讨论的具体的变换，而是用一个符号来标记它。这样在讨论中就包括了许多可能的重要变换。设x为集合X的一个元素，y是集合y的一个元素。那么把x变换成y可写成(2.6)式中是变换算子。如果对每一个元素x有唯一的元素y，反之亦然，那么我们说算子是非奇异的,否那么，算子就是奇异的。如果下面的两个关系成立，那么变换是线性的。(2.7)(2.8)我们只讨论线性变换。例2.4算子，其中表示在z中所有阶次的多项式，计算，其中。试讨论算子是否表示一个线性变换和它是否奇异。把算子用于上，因此因为求导数过程是线性的，所以算子是线性的。为了判断算子是否奇异，可以研究一个简单的元素，如。此时可以看出变换后。是否有唯一的元素x与y对应呢？然而x并不唯一，可知该算子奇异。用矩阵表示一个线性变换，第一步是选择两个基向量组，一组是用来表示x的元素，另一组是用来表示y的元素。设，，……是在x中的一组基，，，……是在y中的一组基。因此x的每一个元素x可表示为(2.9)式中,,..,是x的坐标。y的每一个元素y可表示为(2.10)式中,,..,是y的坐标。于是，线性变换的矩阵表示为〔2.11〕式中A是一个n阶的方阵，y和x分别是由y和x的坐标所排列成的向量，即〔2.12〕当算子作用于基向量时，矩阵A的元素就是由式〔2.6〕所求得的坐标。因此，为了计算A的第j列元素。我们进行变换，即计算，并用基向量(i＝1，2，……，n)来表示y。的系数就是A的第j列元素。例2.5试用矩阵来表例如2.4中的变换。选定下面两组基：〔1〕，，，，而；〔2〕，，，，而。例2.4中的算子是，其中表示在z中所有阶次的多项式。为了计算在算子k的表示矩阵A中的第j列元素，我们需要形成，并用基向量〔j=1，2，3，4〕来表示其结果。假设考虑(1)的一组基，有因此，在选定基下，算子的矩阵表示为为了变换元素x，其中，我们要求出在选定的基下x的坐标于是，利用，可得也可得到，这与例2.4的结果相同。〔第二种情况，同学们自己动手完成。〕基的改变假设A是在基，，……下算子k的阶矩阵表示。并假设只用一组基,即(i＝1，2，…,n)，已算出了矩阵A。这是最重要的情形。在(2.11)式中所给出的联系x和y的矩阵A中,向量x和y的第i个元素(即和)是基向量的乘子。假定新的基向量是(i＝1，…,n)，且与它相对应的坐标是落在向量和中。为了求得把与联系起来的矩阵，第一步是用旧的基向量(i＝1,2..,n)来表示全部新的基向量(i=1,2,..,n)，即计算，(i=1,2,……,n)(2.13)利用这个关系，那么有〔2.14〕式中P是一个阶的矩阵。把上式代入〔2.11〕式，可得〔2.15〕因而得到〔2.16〕在上面的推导中。假设了P的逆矩阵存在，这是规定了已经选取了一组适当的基(i＝1,2…,n)的情形，亦即全部(i＝1,2…,n)是线性无关的情形。应该指出，矩阵P的逆矩阵存在，保证了对于任何向量x和y存在着如式(2.14)所示的唯一的向量和，反之亦然。例2.6在2.5的例子中，用基，，，，得到线性算子的矩阵形式为现在，我们用，，，，而，求该算子的转换矩阵。根据式〔2.13〕的关系，我们可以写出：因此，有其逆矩阵为新的基下的矩阵表示为〔同学们，可以将该结果，与第二种基直接计算的结果比照。〕定义如果满足，那么矩阵P是一个正交矩阵，因此对正交矩阵有。这个定义说明，如果在基的改变时利用正交矩阵，那么有。对于实际应用来说，存在一些易于建立和使用的正交矩阵。下面介绍两种应用比拟广泛的正交矩阵。最常用的一种正交矩阵是旋转矩阵。〔2.17〕其中i和j可以是任意的，但。旋转矩阵的名称可通过讨论n＝2的情况来说明，在这种情况下，〔2.18〕另一种常用的正交矩阵是映射矩阵，〔2.19〕其中V可以是任意的。这个矩阵称为映射矩阵，因为向量Pw是向量w在与V正交的平面上的映射。例2.7确定对应于向量的映射矩阵和对向量的映射。在这种情况下，有向量w的映射Pw为向量模和矩阵模在使用向量和矩阵的迭代求解过程中，我们也需要一个关于收敛性的度量。认识到一个向量或矩阵的“大小”必然会依赖于数组中所有元素的大小，我们既能得出向量模和矩阵模的定义。模是一个数，它依赖于向量或矩阵的所有元素的大小。定义一个n阶向量v的模可写成，它是一个数。模是v的所有元素的函数，且满足以下的条件：(1)，当且仅当时才有(2.20)(2),对于任意常数C(2.21)(3)，对于向量V和W(2.22)关系式(2.22)是三角形不等式。通常使用的有以下三种向量模，它们分别叫做向量的无穷模、1模和2模：(2.23)(2.24)(2.25)现在，我们可以度量向量序列,……,趋向于向量x的收敛性。序列收敛于x的充分必要条件是对于任何一种向量模都满足(2.26)与向量模的定义相类似，我们也能定义矩阵模。定义一个阶矩阵A的模可写成，它是一个数。这个模是A的元素的函数，并且有以下的关系成立：(1)，当且仅当时才有(2.27)(2)，对于任何常数c(2.28)(3)，对于矩阵A和B(2.29)(4)，对于矩阵A和B(2.30)关系式(2.29)是一个与式(2.30)等效的三角形不等式，在向量模的定义中不要求附加的条件(2.30)式，但是为了能在有矩阵乘积时使用矩阵模，那么必须满足式(2.30)。以下是经常用到的一些矩阵模：(2.31)(2.32)；为的最大特征值(2.33)其中对于对称矩阵A，我们有。模称为矩阵A的谱模。像向量序列中的情形一样，我们可以度量矩阵序列A1，A2…，At趋向于矩阵A的收敛性。例如，收敛的充分必要条件是对任一种给出的矩阵模均满足(2.34)在遇到矩阵与向量的乘积时，我们也需要使用模。在这种情况下，为了得到使用模的有用的资料，我们需要把具体的向量模和具体的矩阵模一起使用。矩阵模和向量模可以一起使用，这是由对于任意矩阵A和向量v下面的关系式均成立的条件来决定的：(2.35)式中用矩阵模来计算，和用向量模来计算。有限单元法有限元法公式的建立引言作为一种分析工具的有限元方法的开展，实质上是从数字电子计算机的出现而开始的，一个连续介质问题的数值解法，根本上需要建立和求解一个代数方程组。在计算机上使用有限无法，就有可能用非常有效的方式建立和求解复杂系统的控制方程组。有限元法之所以具有广泛的吸引力，主要是因为它能分析一般的结构或连续体，能够比拟容易建立其控制方程，而且所建立的系统矩阵具有良好的数值性质。有限元法首先是在结构力学问题分析的物理根底上开展起来的，但不久人们就认识到这个方法同样可以用来解决许多其它类型的问题。当给出有限元法的变分公式并在分析中采用了以后，有限元法的普遍性就更为明显了。为了介绍和从物理上阐述这个方法，我们将只研究结构力学问题的解法。最初开展的情况注往是这样，要引证“创造”有限元法的准确日期是相当困难的。但这个方法的根源可追溯到应用数学家们、物理学家们、和工程师们所做的有关工作。虽然方法的原理早已发表，但有限元法通过工程师们所进行的独立开展工作才获得了真正的推动力。最早的重要著作有Turner等人的论文和Argyris与kelsey合写的论文。“有限元”这个名词是Clough在一篇介绍平面应力分析的有限元法的论文中提出来的。自从有限元法开始开展以来，由于人们很快就认识到它的潜力，因此对这个方法已进行了大量的研究工作。目前，有限元法的概念是一个非常广泛的概念。虽然我们只限于结构力学问题的分，但这种方法能用于各种不同的行业。当讨论到变分公式时，这一点就变得特别明显。在实际问题的求解中广泛使用的—种最重要的方法，是基于位移的有限元法。因为这种方法具有简单性、普遍性和良好的数值性质，所以所有较大的分析程序实际上都是根据它来编写的。基于位移的有限元法可以看作是位移分析法的一个推广，而位移分沂法在梁和桁架结构的分析中已经用了许多年。通过一个小的例子来回忆一下位移分析法并同时引入有限元法的根本概念是有好处的．图3-1简单的梁和桁架结构现考虑对图3-1中结构的分析。在位移分析法中我们是把该结构看成是两个梁单元、一个桁架单元和一个弹簧单元的分割体，第一步是计算对应于结构总体自由度的单元刚度矩阵。在这种情况下，对于梁单元、弹簧单元和桁架单元，我们分别有；；；；整个分割体的刚度矩阵可以由各个单元刚度矩阵通过直接刚度法有效地求得。在这个过程中，结构刚度矩阵K是通过各单元刚度矩阵直接相加而算得，即而系统的平衡方程为式中是系统的总体位移向量，R是作用在结构总体位移方向上的外力向量：，在求解结构的位移之前，我们需要利用边界条件和。这意味着我们可以只考虑含有五个未知位移的五个方程，即，式中是从K中删去第一和第七行以及第一和第七列后得到的，而，上面的讨论说明了有限元法的一些重要特点．根本的处理过程是先把整个结构看成为各个结构单元的分割体，计算对应于结构分割体总体自由度的各单元刚度矩阵，然后通过将各单元刚度矩阵叠加的方法形成结构刚度矩阵，求解单元分割体的平衡方程组就得到单元的位移，然后利用它们来计算单元的应力。虽然原来并不认为用位移法分析梁和衍架单元的分割体就是有限元分析，但以后我们将会看到，实际上我们可以把桁架和梁单元看作两种特殊的有限元。在这个定义上，上述的分析就是梁和衍架结构的有限元分析。而用另一种方法，我们可不通过求解平衡微分方程，而是用虚功原理计算其刚度系数。导出表示弹性体平衡的相应方程的一个等效方法是利用虚位移原理，这个原理表示物体处于平衡的要求是：对于强加在该物体上的任意相容的微小的虚位移，总的内虚功应等于总的外虚功，即〔3.1〕式中，，〔3.2〕是作用在弹性体上的体力、外表力和集中力。从未受荷载时的位置开始的弹性体的位移以表示，其中〔3.3〕相应于的应变为〔3.4〕相应的应力为〔3.5〕和表示虚应变和虚位移。介绍了一些根本概念之后，现在我们扼要地表达图3-2中的平面应力问题的有限元分析是有益的。在该问题的分析中首先使用了“有限元法”这个专门术语。在这个问题中，可假定位移和为未知量，由它们可以算出应变和应力值。分析的根本步骤和上述桁架和梁结构的分析步骤一样。图3-2有限元平面分析该问题的求解可按以下的步骤进行：(1)假设每一单元节点i的两个未知位移为和，而和用简单多项式函数来表示，其中的末知参数是单元节点位移。对于图3-3中的单元，未知位移是。(2)利用虚位移原理计算每个单元对应于节点自由度的刚度矩阵。(3)将各单元刚度矩阵叠加得到结构的刚度矩阵，利用边界条件解平衡方程组求出节点位移，然后算出单元的应力。图3-3局部坐标系统中典型的二维四节点有限元上述过程与本节开头所讨论的简单梁和桁架结构的分析在概念上是相同的，但是，它们之间有一个重要的差异。重要的是应认识到在连续介质的有限元分析中在一般荷载条件下我们将只能得到一个近似解，其原因是在每个单元上被假定为位移模式的多项式函数只能近似表示连续介质的精确位移。因此，利用虚功原理所导出的是近似的刚度系数，而所算出的是近似的节点位移和应力。这与本节开头所讨论的简单梁和衍架结构的分析完全不同，在梁和析架结构的分析中，单位位移的多项式函数可精确地表示单元范围内的结构位移，因而所算出的是结构的精确反响。实际上，就一般有限元分析来说，需要近似表示连续介质的精确反响，以提供各种各样的问题和要求，使该领域得到开展。利用虚位移原理建立有限元法的公式上节我们已经指出，最初是作为桁架和梁结构的位移分析法的推广来介绍有限元法的。在利用位移分析法时，一个框架或桁架结构被看作是在结构的节点处相互联结的离散杆元的分割体，首先算出对应于结构总自由度的杆元刚度矩阵，叠加成结构刚度矩阵，形成荷载向量，然后解节点的平衡方程组，求出结构的节点位移。在结构的一般有限元分析中，所用的概念根本上与上述相同。在这种分析中是将结构理想化为许多离散的有限元的分割体，并假定各单元是在有限个节点上相互联结，这些节点相当于在桁架和梁结构分析中的节点，然后假定在每个有限元中的位移是单元节点位移的函数。因此，可通过该分割体的节点位移来表示结构单元的分割体内任何一点的位移，这样就能够以像桁架和梁结构分析中一样的方式来使用虚位移原理，以便导出对应于节点位移的单元分割体的平衡方程组。〔一〕平面应力分析的位移和应变－位移的变换矩阵为了便于说明，再考虑一个平面内荷载作用下悬臂板分析的例子(图3-2)。该结构是处于平面应力的状态，因此可以将平桓理想化为二维平面应力有限元的分割体，如图3-2所示。所需要的根本矩阵是对于分割体的每个单元的位移变换矩阵和应变-位移变换矩阵。为了推导这些矩阵，我们着重研究如图3-3所示的典型单元，并假设局部单元位移和是以局部坐标变量和的多项式的形式给出：〔3.6〕〔3.7〕未知系数也称为广义坐标，它们可以用未知的单元节点位移和来表示。定义〔3.8〕我们可以把式〔3.6〕和〔3.7〕写成矩阵形式〔3.9〕其中〔3.10〕〔3.11〕〔3.12〕对于单元的各个节点，方程(3.9)都一定是成立的。因此，对于所有四个节点，利用式(3.9)，我们有〔3.13〕其中〔3.14〕因而广义坐标是〔3.15〕其中〔3.16〕将式〔3.9〕代入，可以得到〔3.17〕其中〔3.18〕我们用虚功原理来得到单元刚度矩阵。假设是平面应力状态，单元应变为〔3.19〕其中〔3.20〕利用式〔3.17〕，我们得〔3.21〕式中〔3.22〕〔3.23〕单元应力是〔3.24〕假定是各向同性的线弹性材料，对于平面应力状态，有〔3.25〕其中〔3.26〕而和是材料的弹性常数，为杨氏模量，为泊松比。要进行实际的有限元分析，就需要求出式(3.18)和式(3.22)中分别给出的分割体中每个单元的位移和应变-位移的变换矩阵，而算出每个单元的矩阵A，才能求出这些矩阵。至今，我们所考虑的是通过单元局部节点位移来决定的每个有限元的位移。然而，在推导整个单元分割体的平衡方程的时候，通过整个单元分割体节点位移来表示单元的位移是方便的，即对单元m可写出〔3.27〕而对于单元应变和应力类似地可写成〔3.28〕〔3.29〕式(3.27)至式(3.29)中的量与向量有关，而存储整个有限元分割体的总体坐标系下的全部节点位移。作为一个例子，考虑，即确定图3.2中悬臂板理想化后的单元2上位移矩阵。利用图3.2给出的整个单元分割体和图3.3给出的单元的节点位移的定义，我们有〔3.30〕式中是式〔3.17〕中考虑单元2时对应的单元。〔二〕一般公式的建立为了建立基于位移的有限元法的一般公式，我们来考虑一般的三维弹性体，以后再讨论如何将一般公式具体应用到具体问题中。假设三维弹性体上作用有体力、外表力和集中力，其中，并有初始应力，用顶上加一横扛表示虚量。虚位移原理给出下面的关系：〔3.1〕虚功方程(3.1)是弹性体平衡的一种表述。在有限元分析中，我们把弹性体近似表示为离散的有限元分割体，而这些单元在其边界的节点上互相联结。把每个单元内的局部坐标系度量的位移以及应变假定为n个有限元节点上的位移的函数。因此，对于单元m，我们有〔3.31〕〔3.32〕是全部节点的三个位移分量组成的节点位移向量。即是一个3n维的向量〔3.33〕虽然将所有节点位移都列入中，但应该认识到，对一个具体的单元而言，只是该单元的节点位移对单元内的位移和应变的分布有影响。在有限元离散化过程中，利用式(3.31)和式(3.32)，能使从单元矩阵集合为结构矩阵的过程自动地进行〔这种方法称为直接刚度法〕。单元中的应力与单元的应变和初应力的关系为〔3.34〕现在我们利用在每一个有限元内如式(3.31)所表示的位移的假定，可以推导对应于有限元分割体节点位移的平衡方程。首先，我们把式(3.1)以在全部有限元的体积和面积上积分之和的形式重新写出。即〔3.35〕利用式〔3.31〕和〔3.34〕代入，可得〔3.36〕式中F是作用在单元分割体节点上的外力向量，其中单元结合体的节点位移向量与单元无关，可以提到求和号的外面。要从式(3.36)得到需要用来求解节点位移的平衡方程，我们引用虚位移原理，并将单位虚位移依次强加在所有的位移分昼上。采用这一方法，我们有，而对应于节点位移的单元分割体的平衡方程组为〔3.37〕式中〔3.38〕矩阵K是单元分割体的刚度矩阵〔3.39〕应该注意，在式(3.39)中的单元体积积分的求和意味着将单元刚度短阵直接相加，以得到整个单元分割体的刚度矩阵。以同样的方式将单元的体积力向量直接相加可求得分割体的体积力向量，类似地可得到、和。因此，按上面所进行的方式来建立平衡方程组(3.37)，就包括从单元刚度矩阵求得结构刚度矩阵的组合过程；这个过程通常称为直接刚度法。应该指出，到目前为止，有限元系统平衡方程组是对没有支座的单元分割体建立的，所以，下一步是象梁和桁架结构分析那样的方式加上位移的边界条件。即如果我们把平衡方程组(3.37)写成下面的形式：〔3.40〕式中是未知的节点位移，是的边界位移，我们有〔3.41〕式中，未知的节点位移可以从式〔3.41〕算出，那么对应于的节点力为〔3.42〕方程(3.37)是单元分割体静力平衡的一种表述。但应该认识到，所作用的力可能随时间变化。在这种情况下位移也随时间变化，因而式(3.37)是对任一特定时刻平衡的一种表述。但是，当荷载随时间变化时，在大多数情况下必须考虑惯性力，即需要解一个真正的动力问题。利用达朗贝尔原理，我们可简单地将单元的惯性力作为体力的一局部，如果用与式(3.31)中单元位移相同的方式来表示单元加速度，那么总的体力对荷载向量R的奉献为〔3.43〕式中不再包括惯性力，表示节点加速度，而是单元m的质量密度。在这种情况下，平衡方程组是〔3.44〕式中和都是随时间变化的，矩阵是结构的质量矩阵〔3.45〕在实测的结构动力反响中可以发现，振动时能量会有损耗，在振动分析中通常是引入随速度变化的阻尼力来考虑这一因素的。引入作为对体力的另一类奉献的阻尼力，我们可得对应于式(3.43)的表达式〔3.46〕在这种情形下向量不再包括惯性力和随速度变化的阻尼力，是节点速度向量，而是单元m的阻尼特性参数。这时，平衡方程组为〔3.47〕式中C是结构的阻尼矩阵，即形式上写为〔3.48〕实际上，对于一般有限元分割体来说，要确定阻尼参数如果不是不可能也是相当困难的，特别是因为阻尼特性是随频率变化的。由于这个原因，矩阵C一般不是由单元阻尼矩阵组合而成，而是利用整个单元分割体的质旦矩阵和刚度矩阵以及阻尼量的实验结果构造出来的。〔三〕集中质量当考虑单元质量矩阵的推导时，我们应想到惯性力是作为体积力的一局部来考虑的。因此我们也可以通过把总的质量分成相等的几局部集中到节点上来建立近似的质量矩阵，各节点的质量实质上是对应于该节点周围的有奉献的单元体积的质量。利用这个集中质量的方法，我们大体上就可假定对一个节点奉献体积的加速度都是常数，并等于节点上的值。利用集中质量矩阵的一个重要优点是质量矩阵呈对角形，并且以后将会看到，求解动力平衡方程的数值运算可以减少。例3.1计算图3-4的杆单元组合体的集中质量矩阵，，。图3-4两个杆单元的组合体集中质量是建立有限元步骤归纳有限元法与结构力学中的位移法相似。首先将连续体转化为离散化结构，即将连续体代之以仅在节点互相连结的许多单元组成的结构。这种离散化结构类似于结构力学中的桁架或刚架，但其中的单元不一定是杆件，可以是平面块体或空间块体等。然后对此离散化结构按类似结构力学的位移法进行分析，步骤如下：取每个单元的节点位移作为根本未知数。在单元内建立位移模式，。根据几何关系建立应变矩阵，。根据物理方程建立应力矩阵，。其中称为应力矩阵。根据虚位移原理建立单元刚度矩阵，。建立单元节点荷载向量，。建立平衡方程，，其中，，是全部节点位移向量。代入约束条件求解平衡方程，。位移模式中的N称为形态函数。在有限单元法中，各种计算公式都依赖于位移模式，位移模式选择朗恰当与否，与有限单元法的计算精度和收敛性有关。为了使位移模式尽可能地反映物体中的真实位移形态，它应满足以下条件：(1)位移模式必须能反映单元的刚体位移。(2)位移模式必须能反映单元的常量应变。(3)位移模式应尽可能地反映位移的连续性。在单元之间，除了节点处有共同的节点位移值外，还应尽可能反映在单元之间边界上位移的连续性。杆件结构的有限单元法杆件结构的有限单元法又称结构矩阵分析法、其中以矩阵位移法应用最广。这里只讨论用该法分析平面刚架、桁架结构。平面刚架结构平面刚架结构的有限单元法分析，通常采用两端固接的平面固结单元(见图4-1)，现针对这种单元讨论有限元分析计算公式的建立。图4-1平面刚架单元〔a〕坐标系，单元节点位移与节点力向量杆件结构的有限元分析需要建立单元局部坐标系与结构整体坐标系。对平面固结单元(如图4-1所示)，结构矩阵位移法取结构的节点位移为根本未知量。单元节点位移与根本未知量之间满足变形协调条件，在局部坐标和整体坐标下单元节点位移、节点力分别记为〔4.1〕〔4.2〕〔4.3〕〔4.4〕图4-1中各量值所示方向均为正方向。设单元的弹性模量为E，惯性矩为I，截面积为A，杆长为L。〔b〕位移模式，应力与应变矩阵由材料力学知识可知，位移模式可以选择〔4.5〕〔4.6〕注意，可得单元的位移模式〔4.7〕其中，形函数单元的线应变分为拉压应变和弯曲应变两局部(略去剪切变形)，于是有〔4.8〕其中对应的应力为〔4.9〕〔c〕单元刚度矩阵根据虚位移原理可知，在局部坐标系中〔4.10〕具体表达式为〔d〕坐标转换矩阵，整体坐标下单元刚度矩阵从图4-1可知，局部坐标系相对于整体坐标系逆时针旋转了θ角。根据单元的节点位移或节点力在两个坐标系中的投影关系，可得坐标转换关系〔4.11〕于是，单元节点位移或节点力的坐标转换矩阵(从整体坐标转到局部坐标)为〔4.12〕根据正交矩阵的特性和矩阵的相关知识，可知在整体坐标中的单元刚度矩阵为〔4.13〕〔e〕节点荷载向量单元上的荷载按虚位移原理移置到节点上、得到等效节点荷载向量〔4.14〕其中是作用在杆上的均布荷载，是集中荷载。整体坐标中的节点荷载同局部坐标中的关系为〔4.15〕〔f〕结构的平衡方程根据虚位移原理，结构的平衡方程表示为〔4.16〕其中，，分别代表整体总刚度矩阵，总节点位移向量和总节点力向量。平面桁架结构平面桁架结构的有限元分析通常采用两端铰接的平面铰结单元，其坐标系如图4-2所示。〔a〕坐标系，单元节点位移与节点力向量杆件结构的有限元分析需要建立单元局部坐标系与结构整体坐标系。对平面固结单元(如图4-2所示)，在局部坐标和整体坐标下单元节点位移、节点力分别记为〔4.17〕〔4.18〕图4-2平面桁架单元〔4.18〕〔4.19〕〔b〕局部坐标中的单元刚度矩阵〔4.20〕〔c〕坐标转换矩阵〔4.21〕〔d〕整体坐标中的单元刚度矩阵〔4.22〕其中，。平面问题的有限单元法前言杆系在传统上即认为是离散体系，而且它的组成元件(即单元)一般具有明确的受力特性，在材料力学中已经仔细地研究过。从本章开始我们转入连续体系的有限单元法。我们主要是从工程观点出发，来提出有限单元法。也就是说，我们要把分析离散体系的矩阵位移法用之于分析连续弹性体。这就遇到两个新问题：1．由于连续体上没有自然的节点和单元分界面，所以必须人为地设置假设干分割线，将连续体划分成一个个的单元并选定假设干个点作为节点，将它们的广义位移做为根本未知数。由此随之又会产生许多问题，例如单元的形状、大小和节点数日如何确定等等。2．这些单元的受力特性如何?我们是不清楚的。没有现成的公式能把节点位移和单元内任一点处的位移和应力等联系起来。这就需要选择适当的形状函数来建立这种联系。由此又会产生准确度、可靠性等等问题。总之，与杆系比拟，连续体有限单元法的主要问题是单元特性的研究。至于整体分析的原理和步骤(即总刚度方程的组成规律和求解方法等等)，两者是相同的。在本章开始我们将重复表达全部单元分析和整体分析的原理，之后仅限于单元特性的研究。将连续体分割成许多单元并建立起描述其受力特性的各种公式，这项工作称为连续体的离散化。在有限单元法里，离散化的原理和途径有很多种类，我们将只限于介绍最常用的、假定位移函数的位移法这一类。连续体的离散化与三角形常应变单元图5-1表示一平板一端固定，另一端受分布载荷的作用。类似的悬臂板问题，我们在第三章曾经见过，不过当时没有仔细的讨论这样的事。这样的问题看起来似乎并不复杂，但在弹性理论中仍然是一个难题。为了计算它的变形和应力，最好使用有限单元法。〔a〕三角形单元的节点编码图5-1平面应力问题我们用一组网格线把平板划分成假设干三角形区域，每个三角形就算作是一个单元，网格线的交点就算是节点。把这些单元和节点按一定的顺序编号，并把各个节点沿坐标轴的位移取为根本未知量：对于节点1对于节点2对于节点i设有n个节点，那么整体的总节点位移向量为对每个单元来说，它与三个节点发生联系，设其整体序号是，依次标记为，那么此单元的节点位移向量可以表达为〔参看图5-2〕图5-2三角形常应变单元单元的节点位移向量是整体的总节点位移向量的一局部。〔b〕位移模式和桁架不同，这里没有现成的用节点位移向量来表示单元内任一点位移的现成公式，因而必须采用近似假设。我们知道，任一函数常常可以用有限项的幂级数来迫近。所以不妨假设单元上任何一点的位移为〔5.1〕也就是把高阶的项全都略去，只保存线性项，这自然会带来截尾误差。从几何意义上来解释：函数可认为是定义在单元上的某种复杂曲面，而现在用一平面来代表它，这样做当然是有误差的。可是我们从直觉上可以想像得出来，假设把单元划分得很小，那末在足够小的范围内平面和曲面是相差不大的。所以可以用式〔5.1〕来代表真正的位移函数而不致造成很大的误差。现在要用节点位移来表示式〔5.1〕中的系数，将节点坐标代入该式，位移应为节点位移，即〔5.2〕由上式可以得到用节点位移表示待定系数的式子，之后再代回到式〔5.1〕之第一式，就有〔5.3〕其中〔5.4〕，，〔5.5〕为单元的面积。其他的式子可用下述方法得到：将上边四个式子小的下标按“循环交换”的规那么依次更改。所谓循环交换就是的交换方式。同理可得〔5.6〕综合起来，就是〔5.7〕其中〔5.8〕通过以上处理方法，连续体上任一点处的位移皆可用假设干选定的节点位移来表示。一个寻求未知函数的问题已经转化成一个寻求假设干个未知数的问题。这样，我们就完成了将连续体离散化的关健一步。由此可见形状函数在有限单元法中占据的重要地位。式〔5.4〕和式〔5.6〕说明，这种单元的形状函数矩阵是由三个节点处的形状函数组合而成的。每个形状函数都是线性函数，在几何上代表一个平面。在所代表的节点处，函数的值为单位1，在另外两个节点处，函数的值为零。图5-3是这类函数的图形表示(以为例)。图5-3形状函数平面〔c〕应力应变矩阵利用前面的知识，可知应变〔5.9〕其中应变矩阵〔5.10〕；〔5.11〕可见由于采用了线性的形状函数，单元上的应变(以及应力)为常数。所以这种单元也可称之为常应变(或常应力)单元。这当然是近似的，这种应力实际上只是真实应力的某种平均值。对于平面应力问题，应力为〔5.12〕其中为材料的泊松比，为材料的弹性模量。对于平面应变问题，只要将弹性矩阵中的用、用代替即可。将应变代入，可得〔5.13〕其中称为应力矩阵〔针对平面应力情形〕〔5.14〕；〔5.15〕〔d〕单元刚度矩阵根据虚位移原理，由式〔3.39〕可知，单元的刚度矩阵为当平面问题的板厚度为时，上式改写为此处的标记均表示第个单元。将式〔5.9〕到式〔5.15〕代入上式，可得〔5.16〕式中每个子矩阵都是阶的，其下标与单元的节点编号相对应，具体的计算公式如下：；〔e〕节点荷载向量对于连续体来说，在单元的边界上有着分布作用的边界力，它们是其他单元或是结构的外部环境对该单元的作用力。而在节点上，一般却没有集中的节点力，这和杆系结构不一样。不过出于习慨在物理意义上我们常常可以设想把这些边界力向节点处简化，形成一组等价的虚拟的节点力。显然当单元尺寸很小时用这种等价节点力代表真实的边界力不致引起很大的误差。由第三章的知识，我们知道三节点三角形单元的节点荷载向量表示为下面给出几种情况下的节点荷载向量的表达式：设在三角形单元的点受到集中荷载作用，那么〔5.17a〕设单元受方向的体力作用，例如受重力作用。其合力，为容重。〔5.17b〕设单元的边上受有沿方向的集中力，其作用点距离节点的距离分别为，，，那么〔5.17c〕设单元的边上受有沿方向的三角形分布荷载，荷载在1点的集度为，在2点的集度为零，边长为，那么〔5.17d〕设单元的边上受有均匀分布的侧压力，那么〔5.17e〕〔f〕温度应力假定单元有温度增量，那么弹性应变为其中为总的应变，为弹性应变，为温度变化引起的应变。对于平面问题温度变化不引起剪切变形，故而上式中的剪切应变为零。代入物理方程对应第三章式〔3.34〕上式中的相当于初始应力。因此，温度变化对应的平面应力问题的节点力向量为〔5.17f〕对于平面应变问题将弹性常数作变换即可。面积坐标在采用较精密的三角形单元时，利用面积坐标，可以大大简化荷载向量、应力矩阵和刚度矩阵矩阵的推导工作。在如图5-4所示的三角形单元中，任意一点P的位置，可以用如下的三个比值来确定图5-4面积坐标，，〔5.18〕其中A为三角形的面积，、、分别为三角形、、的面积。这三个比值就称为P点的面积坐标。注意，三个面积坐标并不是互相独立的，即〔5.19〕这里所引用的面积坐标，只限于用在一个三角形单元之内，在该三角形之外并没有定义，因而是一种所谓局部坐标。与此相反，以前所用的直角坐标和，那么是一种所谓整体坐标，它通用于所有单元的总体也就是通用于整个结构物。根据面积坐标的定义，在图5.4中不难看出，在平行于边的一根直线上的所有各点，都具有相同的坐标而且这个坐标就等于“该直线至边的距离”与“节点i至边的距离”的比值。图中示出的一些等值线。根据面积的计算公式，可知，〔5.20〕同式〔5.4〕相比，可知三角形单元的形状函数，就是对应的。用矩阵表示为〔5.21〕将上三式分别乘以然后相加，可得同样，用坐标，可得上二式说明了直角坐标与面积坐标之间的关系，写成矩阵形式〔5.22〕对面积坐标函数求导数时，可以利用下面的公式〔5.23〕求面积坐标的幂函数在三角形单元上的积分值时，可应用积分公式〔5.24〕求面积坐标的幂函数在三角形莱一边上的积分值时，可应用积分公式〔5.25〕其中为边的长度。采用高次位移函数的三角形单元在上一节中，介绍了最简的三角形单元。即在每个单元范围内，设位移是线性变化的，应变和应力作为位移的一阶导数那么是不变的常量。这样的单元通常称为低阶(低次)元或是称为常应变(常应力)元。它的公式简单，计算省时，但精度较差。当弹性体内的应力有急剧变化时，除非把网格划分得很密，否那么难以逼近实际情况。可是划分很密时，相应单元的数日大大增加，总计算工作量相应增大，以至于可能不如采用较复杂的形状函数更为合算。从一些实际计算经验看来，采用二次多项式函数去逼近未知的位移场较为理想。它不太复杂，但又有较好的精度，能够以粗细适当的网格去逼近那种变化陡峭的应力场，到达工程上令人满意的精度。〔a〕位移模式假定位移模式按照二阶多项式变化，即〔5.26〕图5-5六节点三角形单元上式中包含着十二个待定系数，为了把它们用节点位移表示出来，每个单元应当有六个节点。通常的作法是在三个顶点之外再取三个边的个点，如图5-5所示。对应的形状函数矩阵是〔5.27〕其中，并且〔5.28〕式中六个系数与三角形的几何尺寸有关式中的下标仍然按照循环。〔b〕应变矩阵对应的应变矩阵为〔5.29〕式中，可见应变是线性变化的，并可以表示为。〔c〕单元刚度矩阵同理，单元刚度矩阵表达为〔5.30〕式中任意一个子矩阵可以写成其中此外，计算时应注意到，分别为三角形单元的体积、静矩、惯性矩和惯性积等等。对于节点荷载向量可以用类似的方法求得，此处不再给出，请同学们自行推导。图5-6Pascal三角形从理论上说，我们还可以采用更高次的多项式函数作为假设的位移函数。为此，要相应地增加节点的数量，除了在边上加点之外，在单元内部也可以设置节点。图5-6是一示意图，叫做帕斯卡(Pascal)三角形。它形象地指明了多项式的次数与节点安排之间的关系。例如，当采取三次多项式时，标明为的三个点就是三角形的三个顶点，在每个边有两个节点，在三角形中间还有一个节点，共计十个节点，正好与双变量的三次多项式的所包含的项数相等。从本节所导出的公式可见，高次单元的计算公式较低次单元复杂，故计算每个单元的特性时，花费的工作量亦较大。不过在同等精度要求下高次单元所对应的网格却比拟粗，即总的单元个数较少，所以总的计算工作量可能比采用低次单元财来得节省。〔d〕用面积坐标表示的高阶单元使用面积坐标，对于高阶单元是非常方便的，以二次单元为例，设〔5.31〕利用式〔5.19〕，上式可以表示为〔5.32〕这说明以面积坐标表示的多项式可以化成三个坐标的齐次型式。我们可以很方便地用节点位移来表示上式中的待定系数。例如，对于节点1。有，而，所以立即得到。最终可得形函数为〔5.33a〕，是临近的角点号〔5.33b〕这种表达式较上节的形式简洁。利用式〔5.23〕，单元的应变矩阵为；〔5.34a〕；〔5.34b〕采用面积处标的另一个突出优点，是便于利用插值函数来直接构成形状函数。方法如下：设我们要求采用n次多项式的位移函数,总节点数应为〔5.35〕在三角形的每一条边上，都有个节点。所有节点沿着任何一条坐标轴的方向(参看图5-4)均可分成层。现在我们约定：的那一层算作是第0层，而后是第1层、第2层……第n层。于是，任何一个节点，例如节点，其位置完全可以由它所在处的三个层序号来定义。设点沿上方向的层序号依次为，那么该点的坐标显然是；；〔5.36〕相应的，将点的形函数记为。现在，我们讨论一下这个函数的性质。根据以往的经历，我们不难知道应当满足两个条件：〔1〕为n次多项式；〔2〕在节点处取值1，其它节点取值为0。具有这些性质的函数可以用下述方法找到，把写成别离变量的形式〔5.37〕并令右端的三个函数具有如下的特征：分别是次的多项式，因而满足上述第一条。在的坐标线上分别取如下的特定值，即在层上取0，层上取1，其它层取任意值。在的坐标线上分别取如下的特定值，即在层上取0，层上取1，其它层取任意值。类似。这三个函数是很容易直接构成的，因为它们显然就是一维的拉格朗日插值函数。；〔5.38〕利用式〔5.36〕可以简化为；〔5.39〕另外，当层序号为0时，为满足前述要求，可直接取就可以了。最后，将式〔5.38〕代入式〔5.37〕就可以球得所要的形函数。例如构成三次单元的形函数〔图5-7〕：对于角点：；〔5.40a〕对于边点：；〔5.40b〕上式中的是两旁的节点，临近。对于中点：〔5.40c〕例如：，等等。最后，我们来指出形状函数的另外一些性质：节点对称性和沿边界的变化规律。节点对称性从上述的形状函数公式中可见，当以层序号计算的下标交换顺序时，函数表达式中自变量的下标也作对应的交换。例如，将的第一下标与第二下标对换，同时在表达式中的换成，即得的式子，等等。由式(5-36)知，以层序号命名的下标说明节点的位置。下标交换次序意昧着节点移动到新的位置。新旧位置相对于角点是“对称的”(即：或者是从某一角点移动到另一角点；或者是从距某一顶点为边长处移动到距另一顶点为边长处；等等)。另一方面函数表达式中自变量下标的交换，意味着函数的图形随着节点的移动一起“旋转”，变成另一个形状类似的曲面。形状函数的上述特征，称作“节点对称性”。值得注意的是，形状函数的节点对称性与位移模式对于坐标的对称性是一致的。沿边界的变化规律如研究以前介绍过的假设干种单元的形函数图形，将会发现：〔1〕内部节点的形状函数沿边界上恒为零。所以，这种函数称做“气泡状函数”或称“零函数”。显然它不影响单元之间的位移协调关系。并且可以使用缩聚的方法将对应的节点自由度消去。〔2〕角点的形状函数在它的两个邻边(即构成该角点的两边上的变化和一维的拉格朗日插值函数相同，在它的对边上恒为零。〔3〕边点的形状函数在所在边上是一个一维的拉格朗日插值函数，在其它两边上恒为零。〔4〕由此可以推知，在单元的任一条边上的位移，仅仅取决于这条边上那些节点的位移值及相应的形状函数。形状函数具备上述性质，对于相邻单元之间的位移协凋有重要意义。矩形单元图5-8矩形单元三角形单元的优点之一是它的“适应性”，任何复杂边界的弹性体，总是可以划分成三角形。可是在规那么边界的情况下，显然划分成矩形更加方便。另外，许多计算实例已经证明，矩形单元的计算精度也比三角形单元好。这是因为在每个小矩形范围内，矩形单元有连续变化的应力场，而对应的两个三角形单元却只有阶梯变化的应力场。所以，在复杂边界的情况下，同时使用三角形和矩形单元将是可取的计算方案。图5-8矩形单元与三角形单元不问，不可能采用完全的多项式作为位移函数。因为一次的完全多项式只有3个待定系数而矩形单元却有4个角点。为使二者的数目一致应在二次的多项式中选取补充项。设我们补充的项，那么在矩形的上、下边界上位移曲线是二次的，它不可能由两角点处的位移来唯一地确定，因而是不相容的。同样的道理，也不可以取的项。所以，应取〔5.41〕这个函数虽然含有二次项，但在单元的任何一条边界上都是线性的，所以只要两点就可以唯一确实定位移，因而满足相容条件。根据上式可以写出：；；；假设改用自然坐标表示〔见图5-8〕，，可得；；或者概括成；〔5.42〕单元的应变为式中为单元节点位移向量。；〔5.43〕对于平面应力问题，单元的刚度矩阵为式中；〔5.44〕其中，h为单元厚度。对于平面应变问题按照前述同样处理。从式〔5.43〕可以看出，在单元内部，应力〔应变〕不是常数。它虽然不是完全的一次多项式表达的，但却是一种线性变化的规律。为逼近更加复杂的应力场，可以采用更复杂的位移模式。例如采用9节点的单元，位移模式可取相应的节点安排和形函数图形见图5-9。图5-9二次的9节点矩形单元这种形函数可以用拉格朗日插值直接构造。设需要在一个矩形域上构造某点p的拉格朗日插值函数，并要求它在x方向是m次的，在y方向是n次的。为此，我们有个节点，它们在x方向分成排，在y方向分成排，设p点在x方向的第i排y方向的第j排，那么〔5.45a〕其中R代表一维的拉格朗日插值函数，具体表达式为；〔5.45b〕将式〔5.45〕用于图5-9的9节点单元，那么有；；；；；；〔5.46〕其中，，，对于的，具有相同的公式。如果将位移模式中的项除掉，那么为8节点的矩形单元。具体的建立方法在等参单元中再详细介绍。虽然对于矩形单元采取不完全的多项式作为位移函数，但是我们仍然可以把图5-8的四节点单元称为一次的或线性的(因为是二维情况，所以又叫双线性的)，把图5-9所示的九〔或八〕节点单元称为二次的。这是因为就每一个坐标方向来说，特别是就每一条边界来说位移的变化是二次的。其他更高次的单元也是如此，并可用图