有限元原理结构矩阵分析(平面桁架平面应力)变分.doc

第二章 结构矩阵分析由于有限元方法起源于力学中的结构分析，本章的作用是通过三个典型问题说明有限元方法应用于结构分析时的一般步骤，并借此了解有限元方法的一些基本概念。2-1平面桁架（直接法，结构矩阵分析中常用的力法，处理静定问题，位移法，可处理静定&静不定）本节讨论的对象是图2-1所示的平面桁架。组成桁架的各杆为等截面直杆，外载荷p直接作用于杆的铰接点（结点）上。为简单起见不妨设各杆的截面积均为，材料的弹性模量均为。我们可按下述步骤求得桁架的变形和内力。pp图y x图、 结构的离散化对结点及单元编号取组成桁架的每根杆为一个单元(该问题本身为一离散结构的力学问题)，以, , 加以编号；取杆的铰接点为结点，以、加以编号（总体结点序号）。如图所示，即：我们所讨论的桁架包括三个单元、三个结点。各单元（杆）仅在结点处连接。、 建立总体坐标系并确定结点坐标和自由度为了描述结构的平衡需要建立一个坐标系，称为总体坐标系，以区别于以后出现的“局部坐标系”。总体坐标系的选择原则上不受限制，但希望使用方便。本节所选的总体坐标系示于图，坐标原点与结点重合。以u, v 分别表示沿 x, y方向的位移分量， p, q分别表示力沿 x, y 轴的力分量（投影）。在总体坐标系中各结点的坐标为：x u py v qijy v qx u pss图（x1, y1）=(0, 0 )、（x2, y2）=(a, a )、（x3, y3）=(a, 0 )它们将作为程序的输入数据（几何参数）。每个结点有两个自由度，对结点、分别为u1， v1T 、 u2 ，v2T 、 u3， v3T若暂时不考虑支承约束条件，整个结构的结点自由度为u1 v1 u2 v2u3 v3T、 单元分析（建立结点力与结点位移之间的关系）取一个一般性的单元，设它的两个结点在结构中的编号为i, j （单元内部的结点序号）。由材料力学可知，杆的轴向刚度为EA/L。其中为杆的长度：（）单元局部坐标系现选取一典型单元对其进行单元分析，对所分析的单元按如下方式建立一个坐标系：原点：与结点i重合,x轴：沿 i ,j方向,y轴：与x轴垂直。如图2-3所示。这个坐标系只属于一个单元，故称为单元局部坐标系，不同单元的单元局部坐标系一般是不相同的。在单元局部坐标系中可以规定：结点自由度 ui vi T， uj vj T ；单元结点自由度 u ui vi uj vj T。（）局部坐标系中的单元刚度矩阵在外载荷作用下，结构发生变形，单元必受到来自结点的作用力。桁架中的杆只承受轴向力，大小与杆的轴向伸长成正比在局部坐标系中这种特性可以得到清楚的表述（这一点也是引入局部坐标系的理由之一）。若以pi , qi , pj , qj 分别表示结点i, j作用于单元的力在 x, y 轴上的投影，由号单元的静力平衡有（图2-3）有(2-1-1)用矩阵的形式可以写成若引入单元广义力矢量： 则上式可缩写为(2-1-2)其中(2-1-3)称为局部坐标系中的单元刚度矩阵,它只与杆的几个参数E、A、L有关，与杆的方位无关。（）坐标变换局部坐标系中的单元刚度矩阵公式简捷。但不同单元的局部坐标系一般不同，为了研究结构整体的平衡，必须将结点给单元的力以及相应的单元刚度矩阵转换到统一的坐标系总体坐标系。在总体坐标系中单元结点自由度 u ui vi uj vj T结点给单元的力r pi qi p j qj T在图2-3中，x 轴与x轴的夹角为 结点的位移分量的坐标变换为单元的位移分量的坐标变换为(2-1-)或缩写为类似，r 与 r 之间的转换关系为(2-1-)由于(2-1-6)是正交矩阵，因此(2-1-7)也是正交矩阵。所以有将（2-1-4）、（2-1-5）代入（2-1-2）有从上式可得到(2-1-8)其中(2-1-9)称为单元在总体坐标系中的单元刚度矩阵。以后将会看到，（2-1-9）是一个具有普遍意义的公式。它表明，当单元的自由度由一种形式换成另一种形式时，单元刚度矩阵只需进行一次相似变换。对于平面桁架单元，将（2-1-3）、（2-1-6）、（2-1-7）代入（2-1-9）可得到更便于应用的单元刚度矩阵公式(2-1-10)（）具体结果由（2-1-10）可求得各单元的刚度矩阵的具体形式如下：单元：单元自由度u1 v1 u2 v2 T ，单元刚度矩阵为 (2-1-11)单元：单元自由度u v u v T ，单元刚度矩阵为(2-1-12)单元：单元自由度u1 v1 u v T ，单元刚度矩阵为(2-1-13)请注意，单元刚度矩阵与单元自由度中位移分量的排列次序有关。如果改动这种排列次序，例如对号单元，将单元自由度次序由u1 v1 u2 v2 T改为 u2 v2 u1 v1 T，必然导致刚度矩阵（2-1-11）元素位置的变动。（）单元刚度矩阵的物理意义和特点设平面桁架单元在总体坐标系中刚度矩阵的一般形式为由（2-1-8），当单元结点位移为1 0 0 0 T时，在单元各结点上施加的力刚好为单元刚度矩阵中的第一列：k11 k21 k31 k41 T。对k的其他各列也可做出类似的解释。即单元刚度矩阵的每一列相当于一组特定位移下的结点力，如表2-1所示。由图2-4可以获得更为直观的理解。表2-1平面桁架单元刚度矩阵的物理意义单元结点位移作用于单元的结点力1 0 0 0 Tk11 k21 k31 k41 T0 1 0 0 Tk12 k22 k32 k42 T0 0 1 0 Tk13 k23 k33 k43 T0 0 0 1 Tk14 k24 k34 k44 Tk11k41k21k311jik12k42k22k321jik44jk24ik34k141k33k43k23jik131图对图2-4中的各种情况，据平面力系的平衡条件应有这三个关系说明，k的四个行向量中只有一个线性独立（四个元素有三个约束方程）。从以上分析可以看出，一般的单元刚度矩阵均具备以下两个特征。（对平面桁架单元而言，从（2-1-10）也可以得出这些结论）（i）单元刚度矩阵是对称矩阵，这是线性系统互易定理的具体体现。由于对称性，对行向量或列向量两者之一得到的结论，对另一个也适用。(ii)单元刚度矩阵是奇异矩阵。它的行向量（或列向量）线性相关，具有零特征值，detk=0。对平面桁架的单元刚度矩阵而言，它的四个行向量（或列向量）中只有一个线性独立，而k有三个零特征值。这三个零特征值对应的特征向量相当于三种独立的刚体位移模式：两个平移，一个旋转。这是我们在单元分析中不考虑位移约束条件的自然结果。、总体刚度矩阵的组装总体平衡方程将图2-1所示的桁架中的支承约束以约束反力代替，如图2-5所示。下面来建立平衡问题的有限元方程。（） 结点平衡条件作用于图2-5每个结点上的外载荷、支座反力以及来自单元的力应处于平衡。R3YpR1X图R1Y以pi(m)、qi(m)示结点i作用于单元m的力在x, y轴上的投影，则单元m给结点 i的力在x, y轴上的投影应为pi(m)、qi(m)。对结点：对结点：对结点；可以合并成(2-1-1)式（2-1-14）的右边为外载荷和支反力。左边则为单元给结点的力，它们是未知的，但可以借助单元刚度矩阵以结点位移来表示。（） 单元刚度矩阵的扩充为了表示（2-1-14）左边的各个列向量，设想将每个单元的自由度扩充到与结构总体自由度相同（本例为），并在单元刚度矩阵中补充零元素，由（2-1-11）、（2-1-12）、（2-1-13）和（2-1-8）可以用结点位移表示（2-1-14）左边的各列向量。由单元(2-1-15)由单元(2-1-16)由单元(2-1-17)（） 组装总体刚度矩阵将（2-1-15）、（2-1-16）、（2-1-17）代入（2-1-14）得到或(2-1-18)(2-1-19)其中上式称为没有考虑位移约束条件情况下的总体刚度矩阵（求和对所有单元进行）。对本节所分析的平面桁架有(2-1-20)要形成总体平衡方程（2-1-18），只需组装出它的右端项和总体刚度矩阵k就足够了，这只是同一件事的两种不同的提法而已。（2-1-19）用来“书写“组装总刚度矩阵的过程是简单而明了的。但事实上不能这样做，将所有km补充零元素扩充为km极大地浪费了宝贵的存贮空间，这些零元素仅起到使km的元素在总刚度矩阵中就位的作用。实际上采用的是另一种方法。如果已选定各结点位移在结构总体自由度中的排列次序为u1 v1 u2 v2 u3 v3 T，对每个单元在形成单元刚度矩阵的同时还形成了个定位数组LM，它将指出单元自由度的各分量在总体自由度中的序号，如下表：单元号单元自由度LM(1)LM(2)LM(3)LM(4)u1 v1 u2 v2T1234u2 v2 u3 v3 T3456u1 v1 u3 v3T1256单元刚度矩阵中第s行第t列的元素kst加到总刚度矩阵的第LM(s)行LM(t)列即可。这一组装总体刚度矩阵的方法被形象的称为“对号入座”。从（2-1-20）可以看出总刚阵k是奇异阵，它的六个行向量（或列向量）中只有三个线性独立。这是尚未考虑位移约束条件，结构的刚体位移未受到限制的必然结果。（） 引入位移约束条件(2-1-21)由图2-1，平面桁架的位移约束条件为 u1=v1=v3=0 代入方程（2-1-18），得到它显然可以分成两个方程组(2-1-22)和(2-1-23)(2-1-24)（2-1-22）就是平常所说的平衡问题的有限元方程，一般把它写成K为考虑了约束条件（2-1-21）后的总体刚度矩阵。U为非约束自由度位移向量，F为作用于这些自由度上的载荷向量。矩阵K相当从矩阵k中划去了第、行和第、列。这说明（2-1-21）这种零位移约束条件，可以通过对总刚度矩阵划行划列的方法来实现。至于其他形式的位移约束条件的实现方法，后面将有专门的一章加以讨论。实际上，在组装总刚度矩阵的过程中，只要不组装应划去的行和列即可直接得到K和F，即方程（2-1-24）。至于方程（2-1-23），尽管求得结点位移后可以由它求得支承反力，但是有限元分析中一般不形成这个方程。如确有必要求支座反力，将另找其他途径。、解有限元方程对于平衡问题，归结为解一个代数方程。本例的自由度较少，从（2-1-22）可以方便地求得非约束自由度的位移被约束自由度根据（2-1-21）可直接赋零、 求单元内力桁架单元的内力只有一个轴力。由（2-1-1）S的大小和正负与 pj相同。由（2-1-2）和（2-1-4）取出它的第三行即得到(2-1-24)对单元、求得内力于下表单元号单元结点位移内力（以拉为正） 0 0 T 0 0T 0 0 0 0 T0由于本节所讨论的桁架是一个十分简单的静定桁架，用理论力学的知识即可得到各杆的内力，结果相同。但是，若桁架为静不定桁架且杆的数目有上千个，那么本节所讨论的方法原则上不会遇到任何困难，我们要解决的课题将转为如何管理有关的大量数据和如何解一个数千阶的代数方程组这样一些技术问题。2-3 平面应力问题常应变三角形qqt图图2-8为一边长为a、厚度为t的正方形薄板。其中AB边固定，BC、CD边自由，AD边作用均布压力q。对这一问题，有限元分析的步骤是：、将ABCD划分（离散）为个三角形（单元），编号。各单元仅在顶点（结点）铰接，结点编号19。建立坐标系后，不难定出各结点的坐标（xi, yi）。、单元分析任取一个一般性的单元，如图所示。三个结点的编号为i, j, k。结点位移为( ui, vi )、( uj, vj ) 、( uk, vk )单元结点位移为（） 假定单元内位移场u, v 是x, y的一次函数y, v, qijk图（）x, u, p(2-3-1)为待定常数，在结点处应有可解出：其中当 i, j, k的位置为逆时针排列时，2恒正，且等于三角形单元面积的两倍。将这些结果代入（2-3-1）有类似可得到可以合并成(2-3-2)（）单元的应变、应力利用（2-3-1）不难求得(2-3-3)其中(2-3-4)在假定单元内位移场u、v是 x, y的一次函数的前提下，单元内的应变和应力将是常数，故这种单元又称为常应变三角元。图（10）xyijk0(a)xyijkF1xF3xF2xF3yF1yF2yxxy0y(b)（）为了在单元内构成均匀应力场，必须在单元的各边施加均布载荷，它们的合力一定作用在各边的中点，如图2-10(a)所示。再将各边上的合力平分到这边的两点结点，由图2-10 (b)不难得出类似可求得 qi、pj、qj、pk、qk 并可合并写成根据单元刚度矩阵k的直观意义，常应变三角元的单元刚度矩阵即为(2-3-5)（）单元和的边上作用着均布的外载荷，可以把它们的合力平分到两结点，如图2-11所示。4414f1f2q14878557qf2f1图（）（）为了组装总体刚度矩阵，每个单元还应形成一个数组LM。元素LM(1)LM(6)分别为 ui、vi、uj、vj、uk、vk 在总体自由度中的序号。由于单元的结点位移和结点力都是用总体坐标描述的，单元刚度矩阵不必再进行坐标变换。、组装总体刚度矩阵和载荷向量。为了实现位移约束条件这些自由度对应的行和列可以不必组装。总体平衡方程为(2-3-6)其中,非约束自由度位移为相应自由度的载荷向量为解方程（2-3-6）可得到各非约束自由度的位移，再由（2-3-3）可求得各单元的应力。2-4结束语、本章所讨论的三具体问题尽管力学背景不同，但分析步骤大体相同。区别仅在于结点参数不同，单元分析公式不同。而组装总体矩阵的方法、处理约束条件的方法以及代数方程组的求解方法完全通用。、从建立有限元方程的方法看，本章属于直接方法一类。这种方法只在简单情况下有效，但是物理意义清晰。尽管今天有限元方法已发展成为一种抽象的数学分析方法，本章提供的直观解释仍然常常被用来说明一些问题。、对于第一、二两个问题（平面桁架和平面框架）本章求的将是精确解（真实解）。而对第三个问题（平面应力问题），本章求得的只是近似解（位移场是假设的），这样就自然引出以下两个问题：（） 所求的近似解与真实解的接近程度如何？或者说，怎样才能接近得更好一些？（） 在进行单元分析时，曾作了两条关键性的假设：（i）单元内位移u、v是坐标x, y的一次函数； (ii)各单元仅在结点处铰接。这样的假定有没有依据？是否必需？这些问题将在后续的章节中讨论。变分原理简介一、引言变分涉及求以函数为变量的函数泛函的极值的问题。我们首先以三个非常简单的极值问题为例，引出泛函和变分的概念。这三个问题是：最短路径问题；最速下落问题；最小旋转曲面问题。（a）最短路径问题 求连接给定两点的最短平面曲线。坐标系取为直角坐标系，其中的一个点位于原点，另一个点的坐标为。如果是一条通过（0,0）和（x1,y1）的曲线，I是两点间曲线的长度，显然或 （1）所求问题为，求函数y，使上述积分值为最小。 （b）最速下落问题 求铅直面内一曲线，当粒子在重力作用下从曲线的一点由静止状态滑落到曲线上另一点时所用时间最少。建立如下图所示的坐标，起点位于原点，终点坐标为（x1,y1）。任一时刻粒子运动的速度为。粒子从原点下落到（x1,y1）所需时间为所求问题为，求函数y，使上述积分值为最小。（c）最小旋转曲面问题 给定两点和一