第五章弹性力学平面问题的有限单元法

1、 严格地说，实际的弹性结构都是空间结构，并处于空间受力状 态，属于空间问题，然而，对于某些特定问题，根据其结构和外力 特点可以简化为平面问题来处理。这种近似，可大大减少计算工作 工作量，为有限元分析提供方便。弹性力学平面问题可分为两类： 如图柱形管道和长柱形坝体，具有如下特点：a纵向尺寸远大 于横向尺寸，且各横截面尺寸都相同；b 载荷和约束沿纵向不变， 因此可以认为，沿纵向的位移分量 等于零。 P x y (1) 平面应变问题： x y 等厚或不等厚平板，具有如下特点：a 长宽尺寸远大于厚度，b 载荷只沿板面，且沿厚度均匀分布，因此可以认为沿厚度方向的应 力分量等于零。 上述两类问题有许多共同

2、特点，合成为弹性力学平面问题。 (2) 平面应力问题： 0,0 Y yx X yx yxyyx x x v y u y v x u xyyx , xyxyxy yxyyxx EE EE 2 1 1)1 (2 , )( 1 , )( 1 2 22 1.平衡方程 2.几何方程 对于平面应力问题 3.物理方程 T xyyx T xyyx , D 2 1 00 01 01 1 2 E D D E 2 1 E 1 若令 则 而 称为弹性矩阵，它是一个对称矩阵，它的元素只与弹性常数与 有关。 换成，把换成 对于平面应变问题，须把 。 vvuu, u S YlmXml xyyyxx , S 位移边界条件 在

3、边界 上 在边界 上 应力边界条件 4.边界条件 YX , S YX , xyyx , * ,vu u S * u * v x v y u y v x u xyyx * * * * * * , U UW 设变形体处于平衡受力状态：体积力为，在自由边界 上的表面力为应力为 设变形体产生虚位移，在固定边界上的位移及 为零，相应的虚应变为 则体积力和表面力在虚位移上作的外力虚功W 恒等于应力在虚 应变上作的虚变形功 即 ， 。 A xyxyyyxx SA dydxtU dstvYuXdydxtYvXuW )( ,)()( * * W SU dydxtdU xyxyyyxx )( * 其中 上面三个积

4、分的意义为： 中的第一个积分表示全部体积力作的虚功；第二个积分表示 中的积分为 它表示单面体四个侧面上的应力在虚应变上作的虚功。 上的表面力作的虚功。自由边界 用有限元法研究实际工程结构的强度与刚度问题，首先要从工 程实际问题中抽象出力学模型，即要对实际问题的边界条件，约束 条件和外载荷进行简化，这种简化应尽可能反映实际情况，使简化 后的弹性力学问题的解答与实际相近，但也不要带来运算上的过分 复杂。 在力学模型简化过程中，必须明确以下几点 判断实际结构的问题类型，是 问题还是 问题；对于 平面问题，是 问题还是 问题。 结构是否 。如果是对称的，要充分利用对称条件，以简 化计算。 简化的力学模

5、型必是 的或 的。 1 力学模型的简化 二维三维 平面应变平面应力 对称 静定超静定 E 在力学模型简化中，还要指出结构的 ：材料的弹性 常数 ，外 大小和作用位置；以及结构的几何 和几何 。 ，泊松比 将已经简化的结构力学模型（或称连续结构）划分成只在一些 点连续的有限个单元，把每个单元看成是一个连续的、均质的、完 全弹性的、各向同性的单元体，把这些连续点称作结点，每个小单 元体称为一个单元；把外载荷按静力等效原理移置到有关受载的结 点上，构成结点载荷，把连续结构进行这样的分割，称为结构的离 散化。对于离散化的每个单元，都认为是符合弹性力学的基本假设 的，因此，弹性力学的基本方程对每一个单元

6、都同样适用。 2 结构的离散化 力学参数 载荷外形 尺寸 一悬臂梁的力学模型简化和单元划分如图： 在确立了力学模型的基础上，再把原来连续的弹性体离散化， 分为有限个单元，这些单元可以是三结点三角形、四结点任意四边 形、八结点曲边四边形等等。单元之间只在结点处相联结。平面问 题的结点为铰结点。完成单元划分以后，需要对所有单元按次序编 号，就得到了有限元的计算模型。 mji, 从现在开始，我们来分析一个典型三角形单元的力学特性。首 先建立以单元结点位移表示单元内各点位移的关系式。设单元的结 点号码为，见图。每个结点在其单元内的位移可以有两个 分量，整个单元将有六个结点位移分量。可用列阵表示为 1

7、单元的位移模式 yx , yaxaavyaxaau 654321 621 ,aaa 由于单元体也是一个二维的弹性体，单元内各点的位移分量是 坐标 按此位移模式，单元内各点的位移可以由单元结点位移通过插值来 获得。这是我们选择的最简单的线性函数，即 式中是待定常数，它可以确定于下： 的函数，在进行有限元分析时，需要假定一个位移模式： (b) T mmjjii T m T j T i e vuvuvu ),mjivu T iii （ ii vu ,ix y 其中子矩阵 式中是结点在轴和 (5-1) (a) 轴方向的位移。 mji,),( ii yx),( jj yx),( mm yx mmmmmm

8、 jjjjjj iiiiii yaxaavyaxaau yaxaavyaxaau yaxaavyaxaau 654321 654321 654321 u mm jj ii mm jj ii mmm jjj iii ux ux ux a yu yu yu a yxu yxu yxu a 1 1 1 2 1 , 1 1 1 2 1 , 2 1 321 设结点 的坐标分别为 、 、 将它们代入(b)式得 联立解(c)式关于的三个方程，可以求得 (c) (d) ， mm jj ii yx yx yx 1 1 1 2 其中 (5-2) mji, mji, 从解析几何知，(2-2)式中的等于三角形 为使求

9、得面积的值不致成负值，结点 转向，如图所示。 的面积， 的次序必须是逆时针 mmmm jjjjiiii uycxba uycxbauycxbau )( )()( 2 1 )( 1 1 , 1 1 , mj m j i mj m j ijmmj mm jj i xx x x c yy y y byxyx yx yx a ),(mji mmmm jjjjiiii vycxba vycxbavycxbav )( )()( 2 1 将(d)式代入(b)式中的第一式，并稍加整理得 其中 (5-3) 同理得到 (e) (f) ),(,)( 2 1 mjiycxbaN iiii mmjjii mmjjii

10、vNvNvNv uNuNuNu ee mji NINININ v u f I mji NNN, N 如令 位移模式(e)、(f)就可以写成 上两式可合并写成矩阵形式如下 式中是二阶单位阵； 形函数矩阵 (5-4) (5-5) (5-6) 位移状态，因而称为 ，矩阵 则称为 。 是坐标的函数，它们反映单元 形函数 x v y u y v x u xy y x 有了单元的位移模式，就可以利用平面问题的几何方程 (g) 2 单元的应变 e mmjjii mji mji bcbcbc ccc bbb 000 000 2 1 e B 求得应变分量，将(e)、(d)两式代入上式即得 或简写成 (5-7)

11、(g) B mji BBBB ),(,0 0 2 1 mji bc c b B ii i i i B mjimji cccbbb, B, yx 其中可写成分块形式 而子矩阵 公式(2-7)是用结点位移表示单元应变的矩阵方程，矩阵 是单元 。由于和 所以中的元素都是常量，因而单元中各点的应变分量 也都是常量，故通常称这种单元为 。 (5-9) 就 等都是常量， xy 常应变单元 应变矩阵 (5-8) D e BD BDS e S S 在得到应变之后，再利用物理方程 便可导出以结点位移表示应力的关系式。把(2-7)式代入上式， 得到 令 则(2-10)式写成 这就是应力与结点位移的关系式，其中 称

12、为 。 (5-10) (h) (5-11) 3 单元的应力 应力矩阵 S ii ii ii ii bc cb cb E BDS 2 1 2 1 )1 (2 2 对于平面应力问题，的子矩阵可写成 (5-13) 矩阵 S可写成分块形式 mjimji SSSBBBDS(5-12) ),(mji ii ii ii i bc cb cb E S )1 (2 21 )1 (2 21 1 1 )21)(1 (2 )1 ( 对于平面应变问题 (5-14) ),(mji mmjjii SSS S 如果注意到(5-1)式，则(5-11)式可写成 从(5-13)、(5-14)式可以看出， 常量。因而，相邻单元将具有

13、不同的应力 和应变。这样，越过公共边界，从一个单元到另一个与它相邻的单 元，应力和应变的值都将有突变，但是位移是连续的（参阅下节）， 常应变单元的这些性质实际上都是由于选取线性的位移模式所造成 的。 (5-15) 个单元中的应力分量也是 中的元素都是常量，所以每 ),(,)( 2 1 mjiycxbaN iiii mm jj ii yx yx yx 1 1 1 2 m j i m j i x x c y y b 1 1 , 1 1 我们在上节讨论常应变三角形单元时，曾提出形函数 其中 mm jj i yx yx a , 坐标轮换 i j m jjjiii cbacba, mmm cba, 2

14、i Ni 1)( 2 1 ),( iiiiiiii ycxbayxN mj , 0)( 2 1 ),( 0)( 2 1 ),( mimiimmi jijiijji ycxbayxN ycxbayxN 由(5-3)式可知，常数 和 依次式行列式 的第一行、第二行和第三行各元素的代数余子式， 根据行列式的性质：行列式的任一行（或列）的元素与其它行（或 列）的元素的代数余子式乘积之和则等于零，从而可以推出形函数 的许多性质如下： 1 形函数 在结点 上的值 而在其余两结点上的值 ,0),( iij yxN 0),(,1),( mmjjjj yxNyxN ,0),( iim yxN 1),(,0),(

15、 mmmjjm yxNyxN 类似地有 ycccxbbbaaa ycxbaycxbaycxba yxNyxNyxN mjimjimji mmmjjjiii mji )()()( 2 1 )( 2 1 ),(),(),( 1),(),(),(yxNyxNyxN mji 2 在单元任一点上三个形函数之和等于1 证明： 根据前述行列式的性质，第一圆括号等于2 ，而第二、第三圆 由此可见，三个形函数中只有二个是独立的。 括号都等于零。故有 ijmij 0),(,),( 1),( yxN xx xx yxN xx xx yxN m ij i j ij i i ij ij ii m m yxx c b y

16、)( 3 在三角形单元的一边上，例如边上有 也就是说，在边上的形函数与第三个顶点的坐标无关。事实上， 边的方程 到(2-4)式，即可得到上面结果（学生自己证明） 代 ijmijn ijij 0),(),(yxNyxN nm 利用这一性质，很容易证明相邻单元的位移，分别进行线性插 值之后 ，在公共边是连续的，例如图中单元 和 具有公共边 ，由(f)式在边上 ji NN ,vu , ji , 式中如(f)式所示，可见在公共边上的位移 公共边的两个结点的位移所确定，所示相邻单元的位移是连续 完全由 连续的。 jjiijjii vNvNvuNuNu, 不论按照哪个单元来计算，根据(5-5)式公共边上的

17、位移均由下 式表示 yyNyNyN xxNxNxN NNN mmjjii mmjjii mji 1 1试证：在三结点三角形单元内的任意一点都有 2试证：在三结点三角形单元mji的一边上，例如ji边上有 0),( ),( 1),( yxN xx xx yxN xx xx yxN m ij i j ij i i 3求所示三角形的二次插值位移模式。该单元有三个主结点， 两个副结点。 ee T mmjjii T T m T j T i e VUVUVU RRRR e mj , T mmjjii e vuvuvu * * f e Nf * 现在来推导单元结点力和结点位移之间的关系。为此，对图5-2 中的

18、典型单元 应用虚位移原理，单元 是在等效结点力的作用下 处于平衡，这种结点力可用列阵表示为 假定弹性体的所有结点都产生一虚位移，单元 的三个结点 的虚位移可表示为 并且，设单元内各点的虚位移 ，具有与5-5(b)式相同的 位移模式，因而有 , i 1 单元的刚度矩阵 * e B * e T e R)( * dydxt T * 而单元内的虚应变 参照(5-7)式，有 于是，作用在弹性体上的外力在虚位移上所作的虚功为 单元内的应力在虚应变上所作的虚功为 (d) (e) e B * e BD t T e )( * dydxtBDB eT T e )( * e T e R)( * dydxtBDB e

19、T T e )( * T e )( * e R eT dydxtBDB 这里假定单元的厚度 为常量，将(d)式和(5-10)式代入，并把 根据虚位移原理：在外力作用下处于平衡状态的弹性体，当发 生约束所允许的任意微小虚位移时，外力在虚位移上所作的功等于 弹性体内应力在虚应变上所作的功，于是，由(e)和(g)式得到单元的 虚功方程为 由于虚位移是任意的，等式两边与 相乘的矩阵应当相 等，故得 提到积分号的前面，则可化成 (g) (5-27) dydxtBDBk T ee kR k D B dydx tBDBk T 记 则上式将有前面曾谈到过的形式，即 这就是表征单元的结点力与位移之间关系的刚度方

20、程。 就是 单元刚度矩阵。如果单元的材料是均质的，矩阵 中的元素是常 量，而且在三角形常应变单元情况下，矩阵 中的元素也是常量， 当单元的厚度也是常量时，再注意到 ，于是(5-28)式 可以简化为 (5-30) (5-28) k mmmjmi jmjjji imijii kkk kkk kkk k 物理意义：单元刚度矩阵 中的任一列的元素分别等于该单元 的某个结点沿坐标方向发生单位位移时，在各结点上所引起的结点 力，它决定于该单元的形状、大小、方位和弹性常数，而与单元的 位置无关。即不随单元或坐标轴的平行移动而改变。 将表达式(5-8)代入(5-30)式，即该平面应力问题中三角形单元 的刚度矩

21、阵，写成分块形式如下 (5-31) srsrsrsr srsrsrsr s T rrs bbcccbbc bccbccbb Et tBDBk 2 1 2 1 2 1 2 1 )1 (4 2 E 2 1E 1 其中 对于平面应变问题，上式中的 应换成 ， 换成 至此单元的力学特性分析已告完成，下面就可以转入结构的整 体分析。 (5-32) 。 e nn e n 12 n T T n TT n 2112 ), 2 , 1(nivu T iii i 假说弹性体被划分为 个单元和 个结点，并对每一个单元 都进行上述运算，则得到 组形如 式的方程，把 这些方程集合起来，便可得到表征整个弹性体平衡的表达式

22、。为此 目的，我们首先引进整个弹性体的结点位移列阵 ，它是由 各结点位移按给定的号码从小到大顺序排列组成的，即 其中子矩阵 是结点 的位移分量。 2 整体的刚度矩阵 ee kR 12 n R T T n TT n RRRR 2112 ), 2 , 1( 11 niVUYXR T n e e i n e e i T iii ee i 再引进整个弹性体的载荷列阵 ，它是移置到结点上的等 效结点载荷，按照点的号码从小到大顺序排列组成的，即 其中子矩阵 是结点 上的等效结点载荷 e R 16 12 n T Te m Te j Te i e n RRRR nmji )()()( 1 12 ),(mjiV

23、UR T e i e i e i 将各单元的结点力列阵 加以扩大，使之成为 阶列 阵 其中子矩阵 是单元结点上的等效结点力 (l) T mmjjii T T m T j T i e VUVUVURRRR mji, i Rmji, T T n TT n e e RRRRR e 21 1 (l)式中的圆点元素均为零，矩阵号上面的 表示在分块 矩阵意义下 所在 的列的位置。这里以假定的次序恰和结 点号码的次序从小到大的排列是一致的。各单元的结点力列阵经过 这样扩大以后便可相加，将全部单元的结点力列阵叠加在一起，就 得到(j)式表示的弹性体的载荷列阵，即 这是因为相邻单元公共边内力引起的等效结点力，在

24、叠加过程 中必然互相抵消，只剩下载荷所引起的等效结点力。 (n) kn2 n m j i kkk kkk kkk k nmji mmmjmi jmjjji imijii nn 1 1 22 将(5-31)式确定的六阶方阵 加以扩大，使成为 阶的方阵 (o) mmmjmi jmjjji imijii kkk kkk kkk k ij ki j 可见，(5-32)式中的22阶子矩阵 被放到(o)式中的第 双 行、第 双列中。 这样，(5-29)式可改写为 e n e nnn Rk 121222 k mji, e n 12 12 n e n ee n e n n e Rk 1 12 1 考虑到 扩大

25、以后，除了对应 双行和双列上的九个子矩阵 外，其余都为零，故上式左边的单元位移列阵 已可用整体的 位移列阵 替代，把上式对 个单元作和，则得 (p) e n e k 1 K ee n e T n e dydxtBDBkK 11 上式左边 是弹性体所有单元刚度矩阵的总和，称为弹性 体的整体刚度矩阵（或称总刚度矩阵），通常都记作 。注意到 (5-28)式，则有 (5-34) 如写成分块矩阵的形式 nnnmnjnin mnmmmjmim jnjmjjjij inimijiii nmji KKKKK KKKKK KKKKK KKKKK KKKKK K 1 1 1 1 111111 (5-35) ),2

26、,1,(, 1 22 nsrkK e n e rsrs nn22 rs Ksr rs K RK n2 n n2 显然，其中子矩阵 它是单元的刚度矩阵扩大到 阶之后，在同一位置上的 子矩阵之和。由于(o)式中很多位置上的子矩阵都为零，实际上， (5-36)式不必对全部单元求和，只有当 的下标 或者属 于同一个单元的结点号码时， 才不等于零，否则都等于零。 (5-37) (5-36) 将(5-34)和(n)式代入(p)式，便可写成 上式包含有关于结点位移的 个线性方程，实际它就是 个结点 上列出的全部 个平衡方程。 0 mn K K 1 1 u0 223221 nn vuuvuv K T nn T

27、 nynxyxyx KKKKKK RRRRRR 1 ,21 , 1241312111 2211 整体刚度矩阵具有很多重要性质 这可以从 中看出，令结点1在坐标轴x方向的位移 而其余的结点位移 ， 便得到结点载荷列阵等于 的第一列元素的列阵， 3 整体刚度矩阵的性质 RK 刚度矩阵 中每一列元素的物理意义为：要迫使弹性体的 某一结点在坐标轴方向发生单位位移，而其它结点位移都保持为零 的变形状态，在所有各结点上需要施加的结点力。 K K 33 K x x K T srrs kK rs n e rss n e T r T r n e T s n e T sr T sr KktBDB tBDBkK e

28、e ee 11 11 )( 刚度矩阵 的主元素总是正的。 例如，由性质1可知，刚度矩阵 中的元素 表示结点2在 方向产生单位位移而其它位移均为零时，在结点2的 方向上必须 施加的力，它自然应顺着位移方向，因而为正号。 刚度矩阵 是对称矩阵 为此只要证明 ，事实上根据(5-32)、(5-36)式有 所以刚度矩阵是对称矩阵。因而在实际计算时，只需计算在对角线 上以及在其一边的元素。 K rs K K rl K rl K rr K K ) 1(2lrB )(lrd K 在讨论(5-36)式的求和号时已经指出，整体刚度矩阵 中第r双 行的子矩阵 ，有些都等于零，只有第二个下标s等于r或者s与r 同属一

29、个单元的结点号码时，才不是零，即在r双行中非零子矩阵的 块数，等于结点r周围直接相邻的结点数目加一。 的元素不是都 填满的，它是一个稀疏阵。 若第r双行的第一个非零子矩阵是 ，则从 到 共有 (r-l+1)个子矩阵，于是 的第2r行从第一个非零元素到对角元共有 个元素，这里B称为整体刚度矩阵的带宽。为了节约 计算机的存储单元，实际计算时应使B尽可能地小。这就要求在结 点编号时，尽量使直接相邻的两结点（属于同一个单元）的号码差 尽可能地地小，以达到 具有较小的带宽。 举例说明： 它是一个稀疏阵，如果遵守一定的结点编号准则，可使非零 元素集中在主对角线附近呈带状。 141414131412 131

30、4131313121311 12141213121212111210 1113111211111110119 7710111010109108 911910999897 81089888786 7978777675 6867666564 5756555453 4645444342 3534333231 24232221 131211 KKK KKKK KKKKK KKKKK KKKKK KKKKK KKKKK KKKKK KKKKK KKKKK KKKKK KKKKK KKKK KKK R R K 因为弹性体在 的作用下处于平衡， 的分量应满足三个静力 平衡方程。这反映在 中就存在三个线性相关

31、的行或列，因而它 是奇异的，不存在逆阵。 K e n e kK 1 我们指出，在排除刚体位移后， 将是正定阵。因为 将(5-30)式代入，得 e n e T tBDBK 1 整体刚度矩阵是一个奇异阵，在排除刚体位移后，它是正定 阵。 T B 12231223 nn e nn BB e n e TT tDK 1 D 0 0 1 e n e T tD 0 K T K 上式右乘 ，并注意到(5-7)式，在集合过程中 将扩大到 32n阶后，有 得 由于弹性矩阵 是正定的，而且t和 都是正的，所以只有在每 个单元中都有 时；才有 ，否则它大于 零。也就是说，当弹性体排除了刚体位移之后，若 ，则二次 型

32、恒大于零，于是 必定为正定阵，有关排除整体刚 度矩阵奇异性的方法将在以后讨论。 用有限元解题，需把弹性体离散化后的每个单元所受的体力、 面力、集中力都移置到有关的受载结点上，形成单元等效结点载荷 列阵 ，这种移置时按静力等效原理进行的。 e R 根据圣维南原理，这种载荷移置所引起的应力误差是局部的， 不至于影响到整个结构，而且随着单元的细分，这种影响将逐渐减 少。但载荷的移置，确实是有限元法的误差来源之一。 所谓静力等效原理，对弹性体来说，是指弹性体上的原载荷与 移置后的结点载荷，在弹性体的任意微小虚位移上所做的虚功相等； 对于刚体，是指刚体上的原载荷与移置后的结点载荷向任一点 简化时，具有相

33、同的主矢和主矩，即对任意坐标轴的载荷投影之和 相等，对任意轴的力矩和也相等。 P j i ),(yxM T yx PPP P R M T vuf * mji, * i e u T mmjji vuvuv * e M Nf * M NM 1 集中力 的移置： 设单元 的边界上 点受有集中力 ，则 利用虚功原理将 移置到三个结点上，并组成单元结点载荷列阵 。 设单元的 点上发生任意微小虚位移 ，则单 元三结点 上的虚位移组成的列阵 由(2-6)式可写出 式中 为形函数矩阵在 点处的取值。 PNPfR T M T eT e T e )()( * PNR T M e ),(mji PNPNR yiMx

34、iM e i 根据弹性体的虚功原理，可写出单元的虚功 故有 则 (5-38) Q mji, ),(yxc T yx QQQ c dydxtt dydxtQ P dydxtQNR Te 2 体积力 的移置。 (5-39) 三角形单元面积。 如果单元 上受有分布状态的体力（单位体积力），且 其上任意一点 的体力为 ，则可将 点的微 小体积 （ 为单元厚度）上的体力载荷 当作一个 小集中力载荷列阵来代替(5-38)式中的 ，然后对三角形单元面 积积分，即得单元的体积移置到结点上的结点载荷列阵。 q),(yxq j i),(yxq T yx yxqyxqq),(),( dlt dltq P j il

35、l Te dltqNR 3 面力 的移置 如果单元 的边上受有任意分布的面力 ，且 则可将微分面积 上的面力载荷 当作一个微小集中载荷 列阵代替式(2-38)的 ，然后对三 角形单元 边长 积分，得 (2-40) T T n TT n e e RRRRR e 21 1 N N 有了单元载荷列向量后，就可以按(n)式迭加到载荷列向量中去， 即 上述载荷移置列阵具有普遍意义，对于其它类型的平面单元也 是适用的，只需把形函数矩阵 换成该种类型单元的 即可。但是，这三个公式在计算 上是相当麻烦的。对线性位移模式的单元，其载荷移置方法也可以 利用刚体静力等效原理进行之。 3 1 关于刚体静力等效原理这里

36、不在赘述；现将四种特殊载荷简化 结果说明如下： 1. 均质等厚度的三角形单元所受的重力，只需把 的重量移置 到每个结点上。 lj i q l tq 2 1 i j 2. 对于作用在长度 的 边上强度为 的均布表面力，只需 把 移置到结点 及 上。 mqj l tq 2 1 3 1 j 3 2 m 3. 线性分布载荷，在结点 处强度为 ，在结点 处强度为0， 则合力大小为 只需将合力的 移置到结点 ， 移置到 结点。 j iPji , l l P 2 l l P 1 4. 在 边上作用在集中力 ，则 点的等效结点载荷分别 为 和 。 有点杠杆原理的味道。 对于一个数值方法，我们总是希望随着网格的

37、逐步细分，得到 的解答收敛于问题的精确解。从上面对于有限单元法的分析中可以 看出，在单元形状确定以后，位移模式的选择是关键。载荷的移置 、应力矩阵和刚度矩阵的建立等等，都依赖于位移模式。很难想象 ，这样一个与真正位移分布有很大差别的位移模式而能得到良好的 数值解。已经证明，对于一个给定的位移模式，其刚度系数的数值 比精确的要大。这样一来，在给定的载荷之下，计算模型的变形比 实际结构的要小。因此，当单元网格分割得越来越细时，位移的近 似解将由下方收敛于精确解，即得到真实解的下界。为了保证解答 的收敛性，要求位移模式必须满足以下三个条件。 一、收敛准则 每个单元的应变一般 总是包含着两个部分：一部

38、分是与该单元中各点的位置坐标有关的 （即所谓各点的变应变）；另一部分是与位置坐标无关的，即所谓 常应变。从物理意义上看，当单元尺寸无限缩小时，每个单元中的 应变应趋于常量。除非我们的位移模式包含着这些常应变，否则就 没有可能收敛于正确解。不难看出，位移模式(b)中与 6532 ,aaaa 有关的线性项是提供单元中的常应变的。 即当结点位移是由某 个刚体位移所引起时，弹性体内不会有应变。这样，位移模式就不 但要具有描述单元本身形变的能力。而且还要具有描述由于其他单 元形变而通过结点位移引起单元刚体位移的能力。例如，在位移模 式(b)中，常数项 就是提供刚体位移的。 41,a a 1位移模式必须包

39、含单元的刚体位移。 2位移模式必须能包含单元的常应变。 。当选择多项式来构成位移模式时，单元内的连续性要求总是满 足的，单元间的协调性要求单元之间不开裂也不重叠，对于在以后 要讨论的梁、板和壳单元，还要求单元之间有斜率的连续性。通常 ，当单元交界面上的位移取决于该交界面上结点的位移时，可以保 证位移的协调性。 在有限单元法中，满足条件1和2的单元，称为完备单元；满足 条件3的单元上称为协调单元或保续单元。我们已经讨论过的三角 形单元，显然同时满足三个条件，因此是完备的协调单元。在某些 梁、板以及壳体分析中，要使单元满足条件3比较困难，实践中也 出现只满足条件1，2的单元，其收敛性也是令人满意的

40、。尤其是放 松条件3的单元，即完备而不协调的单元，已经获得了很多成功的 应用。不协调单元的主要缺点是，不能事先肯定其刚度与真实刚度 的大小关系。但是另一方面，不协调单元一般没有协调单元那样刚 硬（换句话说，一般较柔软），因此可能比协调单元收敛得快。 3位移模式在单元内要连续、并使相邻单元间的位移必须协 调 选择多项式位移模式的阶次时，要考虑到解的收敛性，即要考 虑到完备性和协调性的要求。实践证明，这两项是所要考虑的重要 因素，但并不是唯一的因素。另一个因素是，位移模式应该与局部 坐标系的方位无关，即为几何各向同性。对于线性多项式，各向同 性的要求通常就等阶于必须包含常应变状态。对高次模式，就是

41、不 应有一个偏惠的坐标方向，也就是位移模式不应随局部坐标的更换 而改变。经验证明，实现几何各向同性的一种方法是，根据以下巴 斯卡三角形来选择二维多项式的各项： 二、多项式位移模式阶次的选择 1 常数项 x y 一次项 x2 xy y2 二次项 x3 x2y xy2 y3 三次项 x4 x3y x2y2 xy3 y4 四次项 x5 x4y x3y2 x2y3 xy4 y5 五次项 在二维多项式中，若包含三角形对称轴一边的任一项，则必须 同时包含另一边的对称项。 选择多项式模式要考虑的最后一个因素，就是多项式中的项数 必须等于或稍大于单元边界上的外结点的自由度数。通常是取得与 单元的外自由度数相等

42、，取过多的项是不恰当的。 ),(,mjicb ii 2 解题步骤可归纳为9步： 1划分单元，对结点编号，输入结点坐标； 2对单元编号，按逆时针列出单元三个结点的号码； 3计算载荷的等效结点力； 4计算单元常数 和 ； 5计算各单元的刚度矩阵； 6形成整体刚度矩阵中的非零子矩阵； 7处理约束，消除刚体位移； 8解总刚方程，求得结点位移； 9计算单元应力。 通常，49步由计算机完成，13步用手工完成或计算机完成。 在实现以上各步骤时，为了达到一定的计算精度，节约计算机 存储量，缩短计算机运行时间等目的，还需要注意下列事项。 2 1 4 1 对于具有结构对称、荷载对称和反对称情况，可取对称部分作 为

43、计算模型。如 或 等。对称边的边界条件为垂直于边界的位移 为零；反对称边的边界条件为平行于边界的位移为零。 网格最好也要按对称性划分（对整体而言）。 一对称性的利用 对称 反对称 通常集中荷载的作用点、分布荷载强度的突变点、分布荷载与 自由边界的交界点、支承点等都应取作为结点。如果物体的厚度有 突变点或者物体由不同材料组成时，不要把厚度不同或材料不同的 区域划分在同一单元里。 至于结点多少，要根据计算精度荷机算机容量等综合考虑。从 结果精度来看，当然划分得越细越好，但是，这样做要增加准备工 作和电子计算机的运算时间，甚至超出计算机的容量。因此，在保 证精度的前提下，力求采用较少的单元。故在划分

44、单元时对应力 变化急剧的区域可分得细一些，应力变化平缓得区域可以分得粗一 些。此外，单元三条边的长度不要悬殊太大，以免在计算中出现过 大的误差。 二、结点的选择和单元的划分 不好 好 应尽量使同一单元的相邻结点的号码差尽可能小，以便缩小刚 度矩阵的带宽B=2(d+1),节约计算机存储。若采用带宽压缩存放，整 体刚度矩阵的存储量N最多为 N = 2 n B = 4 n ( d + 1 ) 其中d为相邻结点的最大差值，n为结点总数。 d = 7，N = 448 d = 2，N = 168 三、结点的编号 mji, mji, 四、单元结点 的次序 为了 不出现负值，结点 的次序应为逆时针转向。 若对

45、 取绝对值，次序随便。 对于位移等于零的情况，这是常见的约束形式。 当约束条件为 结点的水平位移 时，则在整体刚度矩 阵 中，应把与位移 对应的 行码和 列码的主对 角线元素改为1，其它元素改为0；右边载荷向量 中，第 个元素改为0。例如： k0 k u K k u12 k12 k R12 k 五、边界条件处理和整体刚度矩阵的修正 0 0 2 0 2 0 0 0 0 0 0 0 10000000 01000000 0000 0000 0000 0000 00000010 00000001 3 3 2 2 4 4 3 3 2 2 1 1 66656463 56555453 46454443 36

46、353433 P P F F F F v u v u v u v u KKKK KKKK KKKK KKKK y x y x 修改后的刚度矩阵为 jj uu j u 12 j 12, 12jj K R j X jjj uAK 12, 12 n jjj n j nnnn njjjjj n n Y uAK Y X v u v u KKK KAKKK KKK KKK 12, 12 1 1 1 1 2,22,21 ,2 2, 1212, 122, 121 , 12 2, 22221 2, 11211 当有结点为给定值 （ 为已知位移值，通常不为 零），则对应刚度方程的第 个方程做如下修改： 主对角线刚

47、度系数 乘以一个大数A（A可取10 左右量 级），右端载荷向量 中 用 取代，其它各项不变 ，即 10 12 j jjj jjjjj uAK uAKvKuK 12, 12 12, 1212, 1211 , 12 12, 12jj AK kj K , 12 nk2, 112 jk jjj uAK 12, 12 jjjjjj uAKuAK 12, 1212, 12 jj uu 修改后的第 个方程变为 由于 （ ，但 ），方 程左端的 项较其它项要大得多，因此近似得到 则有 对于多个给定位移时，则按序将每个位移都做上述修正，得到 全部进行修正后得 和 ，然后解方程即可得到包括给定位移在 内的全部结点

48、位移值。 KR 这个方法使用简单，对任何给定位移（零或非零 ）都适用。 这两种方法引入强制边界条件时，方程阶数不变，结点位移顺 序不变，编制程序十分方便。因此，在有限单元法中经常采用，尤 其是后一种方法。 三角形常应变单元也是常应力单元。算出的应力，通常都作为 单元形心处的应力。 两种方法：绕结点平均法和两单元平均法 1. 所谓绕结点平均法，就是把环绕某一结点的各单元常应力加 以平均，用以表示该结点的应力。 六、应力计算结果的整理 2 cxbxa cba,o o 注意环绕该结点的各个单元面积不应相差太大。边界结点处的 应力应由内结点的应力推算出来。使用抛物线公式 由1、2、3点坐标和应力确定

49、系数。再将 点坐标代入 即可求得 点应力。 2. 两单元平均法 就是把两个相邻单元中 的常应力加以平均，用来表 示公共边界中点处的应力。 同样要求两个相邻单元的面 积不应相差太大。 mKNP/200 E 0 4 1 已知： （单位厚度）， ， 解：1. 力学模型的确定，平面应 力问题。对称性，取 。 123456 x0.00.01.00.01.02.0 y2.01.01.00.00.00.0 2. 划分单元，选取坐标，并对结 点和单元编号。 。求板应力 结点坐标值 2 1 1231232113322 1 323312231 213132321 )( 0,1,1 1,1,0 yxyxyxyxyx

50、yx xxcxxcxxc yybyybyyb 100 020 002 2 00 01 01 1 2 1 2 EE D 3. 求单元刚度矩阵。单元的各几何量如下 101101 001010 010100 000 000 2 mmjjii mji mji bcbcbc ccc bbb B 101101 002020 020200 2 E BDS 101101 020200 103121 121301 002020 101101 4 )1( Et tSBtBDBk TT e k 四个单元，几何形状和尺寸相同，结点编码的方式也相同，因 此，它们的 彼此相同。 1301101 020200 103121

51、 121301 002020 101101 4 4365109 )3( Et k 4 3 6 5 10 9 4. 求整体刚度矩阵 按单元刚度矩阵的总体编号向总刚组集。例如单元单元刚度 矩阵按总体编号为 按上面总体编号的单刚用此编号送入总刚，最后得 1 02 106 1216 103 1213 00416 101216 0120216 10111416 00202 101101 4 Et K 对号入座 同号相加 0 654421 vvvuuu K 0 0 0 0 0 10 220000 261210 016120 021610 012162 000022 4 6 5 3 3 2 1 u u v

52、u v v E 根据约束条件，应有 ， 可划去 中与0位移相对应的1，3，7，8，10，12的行和列，则刚度方 程变为 5. 引入约束条件 求解上面的方程组，即得结点位移 176. 0 176. 0 372. 0 088. 0 252. 1 252. 3 1 6 5 3 3 2 1 E u u v u v v Pa v u v v S xy y x 4 3 3 2 1 )1( 10 44. 0 00. 2 08. 0 0 0 Pa u v S xy y x 4 5 2 )2( 10 0 252. 1 176. 0 0 0 0 0 6. 求各单元应力 Pa v u u v S xy y x 4

53、2 2 3 5 )3( 10 308. 0 372. 0 08. 0 0 0 Pa u u v u S xy y x 4 6 5 3 3 )4( 10 132. 0 372. 0 0 0 0 该程序由一个主程序、七个主要子程序和两个辅助子程序组成 。 一. INPUT：输入数据子程序 1控制数据 NPROB：在一次作业中欲求解的题目数； TITLE：题目的标题，限制为48个字母数字串； NPOIN：结点总数； NELEM：单元总数； NTYPE：问题的类型码（0：平面应力问题，1：平面应变问 题）； NVFIX：受到约束的自由度总数，至少为3； NPFOR：结点荷载数。 2几何数据 NUMEL

54、：单元号； LNODS（NUMEL，INODE），INODE=1，3：单元整体结 点编号（按逆时针） IPOIN：结点号； COORD（IPOIN，INODE），INODE=1，2：结点总体x,y坐 标。 3边界条件 NUFIX（IVFIX）：受约束自由度的总体编号； PRESC（INFIX）：受约束自由度的限定位移值。 4材料性质 PROPS（IPOIN），IPROP=1，4 1弹性模量E，2泊松比 ，3容重 ，4单元厚度 。 p t 5载荷卡 POINT（IPFOR）：结点荷载的数值； NUFOR（IPFOR）：结点荷载的自由度号。 二. 辅助子程序：BCARE（IELEM，USAGE）

55、分离某一单元的结点号，计算该单元的面积。 三. 辅助子程序：STIFP（IELEM，ASK） 形成任一单元的应力矩阵SMATX（3，6）和单刚ESTIF（6，6 ）。 四. PIVOT：总刚的存储 一维变带宽存储。用NPIVO（IGVAB）表示第IGVAB行主对角 元在变带宽一维数组GSTIF（ISTIF）中的地址。 五. ASSEM子程序：总刚的组集 KR KR 六. LOAD子程序 形成等效结点荷载矩阵GLOAD（NGVAB），并能处理集中力 和边界面力。 七. TREAT子程序 用行列修正法，通过引入给定的位移边界条件修正总刚，消除 原总刚的奇异性。 八. SOLVE子程序 采用修正的平

56、方根法求解 ，即LU分解。解出的结 点位移 存放在数组GLOAD中。 九. STREPS子程序 根据结点位移求各单元的应力。 基本要求： 1. 读懂各程序段； 2. 会填数据卡（按格式）； 3. 会用Fortran编辑软件编译、连接、执行程序； 必须做到： 1. 写出程序总框图； 2. 算一个例子，并给出计算结果； 3. 编写计算应力子程序； 4. 将原程序求解总刚方程的1-0法转换成大数法。 23 qj x aa a qm i j m y 第二次作业：第二次作业： 4对平面应力状态下的三结点等边三角形单元，其边长为a， E = 1，= 1/4 （a）写出B，并用矩阵乘法求S和ke。 （b）求

57、出ke中每行每列元素之和，并说明其原因。 （c）设单元发生结点位移 （1）ui = uj = um = 1, vi = vj = vm = 0 （2）ui = uj = vi = 0, vj = 1, um = - , vm = 1/2 试分别求出单元中的应力并说明其原因 （d）设该单元在ij边上受有线性分布 压力，已知qj和qm。试求等效结点荷载。设t = 1 5对于图示的离散结构。试求结点1, 2的位移及铰支座3, 4的 内力（按平面应力计算，= 1/6，E、t已知）。 1 2 3 4 5 6 7 8 1m1m 1m 1m P 简单的三结点三角形单元能较好适应各种形状的边界，但它是 一种常

58、应变单元，单元内应力都是常数，这样便不能较好地反映实 际应力的变化。为了弥补这一缺点，可采用三角形高次单元，即在 三角形单元的三条边上各增加一个节点，如图所示称为六节点三角 形单元。它具有12个节点位移分量，采用完全二次多项式的位移函 数。这种单元的应变和应力都 不再是常量，而是按线性变化 的，能更好地反映实际的应力 状态。这种单元不仅适应性好， 而且具有计算精度高的优点。 i j m 1 2 3 m m j j i i LLL mji , 一面积坐标 对于六节点三角形单元，用面积坐标代替一般的直角坐标，不 仅可以简化应力矩阵、刚度矩阵和载荷矩阵等的运算，而且它不随 三角形单元形状和方位改变，

59、对于计算机的应用也十分有利。 如图所示的三角形单元ijm，任意一点P ( x , y )的位置，可以用 如下三个比值来确定 式中 为三角形ijm的面积， 分别是 三角形Pjm、Pmi、Pij的面积。这三个比值称 为P点的面积坐标。 （5-41） mji 1 mji LLL 0,0,1 mji LLL0,1,0 mji LLL 1,0,0 mji LLL 显然，三个面积坐标并不完全是独立的，由于 所以 根据面积坐标的定义，不难看出，在平行jm边的直线上的所有 各点，都有相同的Li坐标，并且这个坐标就等于“该直线至jm边的 距离”与“结点i至jm边的距离”的比值。图中示出了Li的一些等值 线。容易

60、看出，三个结点的面积坐标分别为 结点i：；结点j： 结点m： 。 （5-42） （5-43） ； )( 2 1 1 1 1 2 1 ycxba yx yx yx iii mm jji )( 2 1 ycxbaL iii i i 面积坐标和直角坐标之间的关系： 三角形Pjm的面积是 于是 （5-45a） （5-44） 类似 )( 2 1 , )( 2 1 ycxbaLycxbaL mmmmjjjj （5-45b） 它们的矩阵形式可写为 y x cba cba cba L L L mmm jjj iii m j i 1 2 1 与(5-4)式对比，可见前述三角形常应变单元中的形函数Ni、Nj、 N