结构有限元分析

习题一习题一 （第一章 杆系结构有限元分析的基本原理） 1.1 试用材料力学方法建立式（1-56）所示的单元平衡方程。 提示：利用式（1-50）和材料力学公式 yEIM ，yEIQ ，并注意 按照材力规定， 单元两端截面上的剪力和弯矩的正方向如题图 1-1 所 示。 , 题图 1-1 解：先求应变能,记 i，j 端的位移列向量 = ,应变能 = 1 2 0 (d 2y dx2) 2 dx 根据结点位移构造单元的位移函数 y，设 y = H1U+ H2+ H3V+ H4= He H = H1 H2 H3 H4 H为 Hermite 插值函数，满足以下条件： H1(0) = 1 H1 (0) = 0 H1(l) = 0 H1(l) = 0 H2(0) = 0 H2 (0) = 1 H2(l) = 0 H2(l) = 0 H3(0) = 0 H3 (0) = 0 H3(l) = 1 H3(l) = 0 H4(0) = 0 H4 (0) = 0 H4(l) = 0 H4(l) = 1 对于每个 Hermite 插值函数H都可以假设为三次多项式 H= a0+ a1x + a2x2+ a3x3 i=1,2,3,4 H满足插值条件，例如： i=1 时，H1= a10+ a11x + a12x2+ a13x3 H1 = a11+ 2a12x + 3a13x2 代入插值条件得：a10= 1； a11= 0 ； a10+ a11 + a122+ a133=0 a11+ 2a12 + 3a132=0 解得 a10= 1； a11= 0；a12= 3 2；a13 = 2 3； 故 H1= 1 3 (3 3x2+ 2x3)； 同理得 H2= 1 3 (3x2 2x3)；H3= 1 2 (2x 2x2+ x3)； H4= 1 2 ( 1 2 x3 x2) 令 = x ，则 H1= 1 32+ 23；H2= 32 23； H1= ( 22+ 3)；H4= (3 2) 外力势能 W = (QV+ M+ QV+ M) = = Q M Q M 故 = U + W = 1 2 0 (d 2y dx2) 2 dx = 1 2 I 1 0 d 2H d2 2 d 1 3 = 1 2 I 3 1 0 (d 2H d2 ) 2 d 根据势能驻值原理得 eT = I l3 1 0 d2H d2 d2H d2 ed Fe = 0 因此若令 = I l3 1 0 d2H d2 d2H d2 d 则上式变为 = Fe 单元刚度矩阵 = I l3 1 0 d2H d2 d2H d2 d 如求k11，此时H1=1 32+ 23 d2H d2 = 6 + 12 故k11 = I l3 1 0 (6 + 12)(6 + 12)d = 12EI 3 再如求k23，H2= 32 23； d2H d2 = 6 6 H1= ( 22+ 3)；d 2H d2 = (4 + 6) k23 = I l3 1 0 (6 6)(4 + 6)d = 6EI 2 同理用公式k = I l3 1 0 d2H d2 d2H d2 d可计算其余单元刚度矩阵元素 求得单元刚度矩阵为： ke = I 3 126126 642622 126126 622642 平衡方程为 = Fe. 1.2 在式（1-55）所示平面梁单元的弯曲刚度矩阵 e K中存在以下关 系 3 = 1 1 = 2 + 4 ( = 1,2,3,4) 说明以上关系式的物理意义。 解：式（1-55）为对称矩阵 = = 3 126126 642622 126126 622642 3 = 1 的物理意义：在单元的近端（远端）发生单元线位移 时引起的近端（远端）沿该线位移方向的力与在远端(近端)引起的同 方向的力是一对平衡力， 由于在单元坐标系下结点力都规定沿坐标系 正方向为正，故方向相反需加负号。 1 = 2 + 4 的物理意义：第 j 个位移分量发生单位位移时， 引起的近端力为 = 1 2 ，引起的远端力 = 3 4 ，根据单元力矩的平衡，对远端取矩得 1 + 2 + 4 = 0，故1 = 2 + 4 ，上式关系表明发生单 元位移分量时，引起的单元的力和力矩是平衡的。 1.3 某个杆件的单元局部坐标系如题图 1-2 所示， 平面与格栅所 在平面重合，轴正方向与整体坐标系轴相同。写出局部坐标 系下的格栅单元刚度矩阵。 , , 0 ,Q , ,Q 题图 1-2 解：由于该杆件只有沿方向的线位移、绕轴角位移和绕轴的角位 移，所以根据空间梁单元的单元刚度矩阵，划去第 1,2,6,7,8,12 行 和第 1,2,6,7,8,12 列对应的元素，即得该栅格单元的刚度矩阵： 66= 12061206 0 00 0 60(4 + 2)260(2 2)2 12061206 0 00 0 60(2 2)260(4 + 2)2 其中 = 622 23 , = 22 2(1+2) . 1.4 用本章1-10 所述杆系结构有限元法的求解步骤求解题图 1-3 所示各题。为便于与本书给出的答案相比较，建议每个梁单元的局部 坐标系均取为与整体坐标系相同， 且一律取单元左端结点为端。 要求给出： 整体坐标系下的结点总位移矢量 各单元在局部坐标系下的单元结点力 如有跨间荷载或变温荷载，应给出各单元局部坐标系下的等 效结点荷载 ， 绘制变形示意图 绘制梁的弯矩分布（弯矩正方向以使梁的下缘受拉，上缘受 压为正）全梁划分为一个单元。 () 12 , 1 2 () 12 温度变化=0/0.5 ) 12 O O /2 /2 =2o/ 题图 1-3 解：(a)图：单元刚度矩阵 = = = 3 126126 642622 126126 622642 结点位移列矢量： = 0 0 2 2T； 结点荷载列矢量： = = F1 M1 0T； 由于 = = 3 126 642 ， = 2 2T， = 0 T 3 (122 62) = 3 (62+ 422) = 0 2 = 33 2= 22 故 = 0 0 33 22 T 所以由 = 得结点 1 处的结点力为： F1= 3 (122+ 62) = F2= 3 (62+ 422) = 因此单元结点力为： = 0T 变形图和弯矩图： (b)图：单元刚度矩阵同(a) = = = 3 126126 642622 126126 622642 结点位移列矢量： = 0 0 2 2T； 直接结点荷载列向量： = = P 1 M1 0 0T 等效结点荷载列向量： = 1 2 1 12 2 1 2 1 12 2 T = + = F1 1 2 1 12 2 T 根据平衡方程得： 3 (122 62) = 1 2 3 (62+ 422) = 1 12 2 2= 48 2= 36 代入平衡方程得： F1= 3 (122+ 62) = 1 2 = 3 (62+ 222) = 5 12 2 因此结点位移矢量为： = 0 0 48 36 T； 单元结点力为： = = 1 122 0 0 T 变形图和弯矩图： (c)图：单元刚度矩阵 = 00 00 0 12 3 6 2 0 12 3 6 2 0 6 2 4 0 6 2 2 00 00 0 12 3 6 2 0 12 3 6 2 0 6 2 2 0 6 2 4 结点位移列向量 = 0 0 0 2 2 2T = 0 1 0.5 d 2 2 = 20 d 2 2 = 0 = 0 1 0.5 2d 2 2 = 20 2 1 3 3 8 = 1 6 02 等效结点荷载列向量 = 0 0 1 6 02 0 0 1 6 02 T 由于单元无结点荷载和跨间荷载，所以2 = 2 由单元的平衡方程得： 2= 0 12 3 2 6 2 2= 0 6 2 + 4 2= 1 6 02 2= 0 2= 022 12 2= 02 6 故 = 0 0 0 0 02212 02 6T 所以 1 号结点杆端内力为： N1= F1= 2= 0 Q1= F2= 12 3 2+ 6 2 2= 0 M1= F3= 6 2 + 2 2= 1 6 02 因此 = 0 0 1 6 02 0 0 1 6 02 T 0 0 1 6 02 0 0 1 6 02 T = 0 0 0 0 0 0T 即仅在温度变化情况下只发生变形，不引起内力。 变形图和弯矩图： (d)图：单元刚度矩阵同(c) = 00 00 0 12 3 6 2 0 12 3 6 2 0 6 2 4 0 6 2 2 00 00 0 12 3 6 2 0 12 3 6 2 0 6 2 2 0 6 2 4 等效结点荷载列向量 = 0 0 1 6 02 0 0 1 6 02 T 结点位移列向量 = 0 0 0 2 0 2T 由单元的平衡方程得： 2= 0 4 2= 1 6 02 2= 0 2= 1 24 02 故 = 0 0 0 0 0 02 24T 所以 1 号结点杆端内力为： N1= F1= 2= 0 Q1= F2= 6 2 2= 6 2 ( 1 24 02) = 02 4 M1= F3= 2 2= 2 ( 1 24 02) = 02 12 2 号结点杆端内力为： Q2= F4= 6 2 2= 6 2 ( 1 24 02) = 02 4 因此杆端内力为： = 0 02 4 02 12 0 02 4 02 6 T 0 0 02 6 0 0 02 6 T = 0 02 4 02 4 0 02 4 0 T 变形图和弯矩图： 1.5 题图 1.4 为长度为L，弯曲刚度为的梁，左端固定，右端为活 动铰支座，跨中受有集中力 P 。给出 1.4 题的 5 点要求。 图 解： （1）由单元刚度矩阵组集结构刚度矩阵： 设 = 2 K = K = 3 126126 642622 126126 622642 写成分块形式有： K = K11 K12 K21K22 = 3 126126 642622 126126 622642 K = K22 K23 K32K33 = 3 126126 642622 126126 622642 刚度矩阵膨胀后 K = 3 1261260 6426220 1261260 6226420 00000 K = 3 0 0000 0126126 0642622 0126126 0622642 则结构的刚度矩阵： K = 3 12612600 64262200 126240126 622082622 00126126 00622642 结点位移列向量 = 0 0 2 2 0 3T 结点荷载列向量 = Q1 M1 0 Q3 0T 由于跨间没有分布荷载和温度变化，则 = = 0 则由结构的平衡方程得 3 (242+ 63) = P 3 (822+ 223) = 0 3 (62+ 422+ 423) = 0 2 = 73 96 2= 2 32 3= 2 8 故 = 0 0 7396 232 0 28T 求单元结点力： F = Q1 M1 = K 所以 Q1= 3 (122+ 62) = 11 16 M1= 3 (62+ 222) = 3 8 = 6 16 同理得： Q3= 3 (122 62 63) = 5 16 所以 = 11 16 6 16 0 5 16 0 T 1.6 题图 1-5 所示为一由个杆件构成的平面铰结杆系， 假定各杆均具 有相同的拉压刚度EA，并已给定所有结点的坐标(，)，求结点 在整体坐标系下的位移矢量 = UA，VA T。 12 , n 21 题图 1-5 解：设第根杆的长度为，与水平线夹角为，则该杆的单元刚度 矩阵为： = 1010 0000 1010 0000 坐标变换矩阵为： = cossin00 sincos00 00cossin 00sincos = R 0 0R R = cossin sin cos 则杆在整体坐标系下的刚度矩阵为： = T = 22 22 22 22 其中 C = cos；S = sin 在刚度矩阵的组集和叠加中只有 A 点的各杆的刚度进行相加， 并在求 解方程中予以保留，而其余的不参与求解，所以 = 2 2 其中表示对所有杆与水平夹角的相应值求和； 故 = 2 2 1 = 1 22()2 2 2 1.7（a）在式（1-181）中，若给定强迫位移= 0，能否同时给定 相应的荷载分量？应如何确定。 （b）题图 1-6 所示两跨连续梁的中间支座 B 发生下沉，已知沉陷 量为。假设两跨梁具有相同的跨度 L 和弯曲刚度，按本章所述强 迫位移的处理方法， 求整体坐标系下所有其它未知的结点位移分量和 各支座反力。 z 题图 1-6 解： （a）可以同时给定相应的荷载分量,相当与在该支座处设置了 一个弹簧在力作用下， 发生了位移0， 设弹簧刚度为k， 并充分大， 并假设该外力加在第个结点上，则刚度方程变为 K11+ + K111+ (K+ K)+ K= 其中= 0 而其余的刚度方程不作修改 （b）单元刚度矩阵为： K = K = 3 126126 642622 126126 622642 则结构的刚度矩阵为： K = 3 12612600 64262200 126240126 622082622 00126126 00622642 荷载列向量 = RA 0 RB MB RC 0 T 位移列向量 = 0 0 0 T 引入强迫位移 0= ，并对刚度方程作修改，可以得到 2 K230= K211+ K222+ 03+ K244+ K255+ K266 故 60 3 = 3 42 = 30 2 = 3 2 6 K630= K611+ K622+ 03+ K644+ K655+ K666 故 6= 42 = 3 2 所以结点的位移分量为 = 0 3 2 0 0 3 2 T RA= 3 (6+ 12) = 3 3 RB= 3 (6+ 24 + 6) = 24 3 MB= 3 (22+ 22) = 0 RC= 3 (12 6) = 3 3 所以各支座反力为：RA= 3 3 ；RB= 24 3 ；RC= 3 3 习题习题二二 （第二章 结构刚度矩阵的贮存和组集） 2.1 对题图 2-1 所示平面桁架建立其结构总刚度矩阵66和结点荷载 总列阵61，假定各杆的 EA 相同。 3 4 1 2 5 题图 2-1 解：各单元的刚度矩阵 K = K = K = 1010 00 00 1010 0000 对于、号单元，与其 X 轴夹角为90，则 R = 0100 10 00 0001 0010 故 K = K = T 1010 00 00 1010 0000 = 0000 01 01 0000 0101 对于、号单元 R = cossin00 sincos00 00cossin 00sincos 其中cos = 2+2；sin = 2+2 故 K = K = T 1010 00 00 1010 0000 = 22 22 22 22 其中 C = cos ;S = sin; = 2+ 2 初始ID表时仅考虑线位移 ID = 1 11 001 011 001 101 则最终ID表：ID = 0 00 120 300 450 060 单元定位数组： 单元 LM =0 00120 单元 LM =1 20450 单元 LM =0 00450 单元 LM =0 00300 单元 LM =3 00450 单元 LM =4 50060 单元 LM =3 00060 各单元刚度矩阵膨胀后的位置，只写出非零元素 单元 II=1 JJ=2 单元 II=2 JJ=4 K = 00 00 0101 0000 0101 K = 1010 00 00 1010 0000 单元 II=1 JJ=4 单元 II=1 JJ=3 K = c2cs c2cs css2 c2s2 c2csc2cs css2css2 K = 1010 00 00 1010 0000 单元 II=3 JJ=4 单元 II=4 JJ=5 K = 00 00 0101 0000 0101 K = 1010 00 00 1010 0000 单元 II=3 JJ=5 K = c2cs c2cs css2 css2 c2csc2cs s2s2css2 由定位数将各单元对应的元素分别送入总刚中得 66= 1 00 1 00 0 1 + 2 00 0 1 + 2 0 1 2 1 00 2 + 2 0 00 1 2 + 1 0 0 2 00 2 其中 C = cos ;S = sin; = 2+ 2 结点荷载列向量 61= 0 1 2 0 0 1 2 T 2.2 题图 2-2 所示为一简单平面刚架，柱子和横梁的拉压刚度与弯曲 刚度分别为1，1和2，2，建立其总刚度矩阵66和结点荷 载总列阵61。 3 4 1 2 2 1 2 ,2 1 ,11 4 题图 2-2 解：单元刚度矩阵 K = K = = 010000 100000 001000 000010 000100 000001 1 1 00 1 1 00 0 121 13 61 12 0 121 13 61 12 0 61 12 41 1 0 61 12 21 1 1 1 00 1 1 00 0 121 13 61 12 0 121 13 61 12 0 61 12 21 1 0 61 12 41 1 010000 100000 001000 000010 000100 000001 = 121 13 0 61 12 121 13 0 61 12 0 1 1 00 1 1 0 61 12 0 41 1 61 12 0 21 1 121 13 0 61 12 00 61 12 0 1 1 00 1 1 0 61 12 0 21 1 61 12 0 41 1 K = 2 2 00 2 2 00 0 122 23 62 22 0 122 23 62 22 0 62 22 42 2 0 62 22 22 2 2 2 00 2 2 00 0 122 23 62 22 0 122 23 62 22 0 62 22 22 2 0 62 22 42 2 初始ID表: 最终ID表: ID = 111111 001110 111111 001110 ID = 000000 120003 000000 450006 定位数组： = 0 0 0 0 0 0 1 2 0 0 0 3 = 1 2 0 0 0 3 4 5 0 0 0 6 = 0 0 0 0 0 0 4 5 0 0 0 6 根据单元定位数组形成总刚度矩阵: K66= 2 2 0 61 12 2 2 00 0 1 1 + 122 23 62 22 0 122 23 62 22 61 12 62 22 41 1 + 42 2 0 62 22 22 2 2 2 00 2 2 0 61 12 0 122 23 62 22 0 1 1 + 122 23 62 22 0 62 22 22 2 61 12 62 22 41 1 + 42 2 结点荷载列向量 单元的等效荷载： = = 0 1 2 1 12 2 0 1 2 1 12 2 T 所以结点荷载列向量为： P61= 1 1 2 1 12 2 1 1 2 1 12 2 T 2.3 对题图 2-1 所示平面桁架和结点编号，采用等带宽方式存贮其总 刚度矩阵K66时，式（2-16）给出的等于多少？它与K的实际 最大半带宽 是否相等？原因何在？这对计算结果有无影响？ 3 4 1 2 5 题图 2-1 解：按式（2-10）计算得： = (1 + ) = (1 + 3) 2 = 8 实际最大半带宽： = 5 是在结构无约束下的刚度矩阵的最大半带宽，而 是 引入约束下的刚度矩阵的最大半带宽。 2.4 对题图 2-1 所示桁架和结点编号，采用变带宽一维数组A存贮其 总刚度矩阵，规定各杆单元均取小号端为端。按本章方法计算的 各列列高 j h，并参照图 2-15 的形式表示，和。计算所得的 j h 与的实际列高 j h是否相等？原因何在？这对计算结果有无影响？ 3 4 1 2 5 题图 2-1 解：将各单元的定位数重写： = 0 0 0 0 0 0 1 2 0 0 0 0 = 1 2 0 0 0 0 4 5 0 0 0 0 = 0 0 0 0 0 0 4 5 0 0 0 0 = 0 0 0 0 0 0 3 0 0 0 0 0 = 3 0 0 0 0 0 4 5 0 0 0 0 = 4 5 0 0 0 0 0 6 0 0 0 0 = 3 0 0 0 0 0 0 6 0 0 0 0 列高计算： 1= 1 2= 2 2+ 1 = 2 1 + 1 = 2 3= 3 3+ 1 = 3 3 + 1 = 1 4= 4 4+ 1 = 4 1 + 1 = 4 5= 3 5+ 1 = 5 1 + 1 = 5 6= 3 6+ 1 = 6 3 + 1 = 4 K = 1 00 1 00 0 1 + 2 00 00 1 + 2 0 1 2 000 2 + 2 0 0000 2 + 1 0 00000 2 A = (1) 00(8)00 0(2) (4)(7)0(16) 00(3)(6)(11)(15) 000(5)(10)(14) 0000(9)(13) 00000(12) 2.5 写出题图 2-3 所示刚架结构总刚按变带宽存贮时： (1) 各列列高 j h (2)主元在一维排列中的地址MAXA( ) (3)变带宽存贮一维数组总容量 S 1 4 2 3 1 题图 2-3 解：初始 ID 表： 最终 ID 表： ID = 111110 001110 001110 011110 ID = 000001 230004 560007 800009 单位定位数组： = 0 0 0 0 0 1 2 3 0 0 0 4 = 2 3 0 0 0 4 5 6 0 0 0 7 = 5 6 0 0 0 7 8 0 0 0 0 9 （1） 1= 1 2= 2 2+ 1 = 2 1 + 1 = 2 3= 3 3+ 1 = 3 3 + 1 = 1 4= 4 4+ 1 = 4 1 + 1 = 4 5= 5 5+ 1 = 5 2 + 1 = 4 6= 6 6+ 1 = 6 2 + 1 = 5 7= 7 7+ 1 = 7 2 + 1 = 6 8= 8 8+ 1 = 8 5 + 1 = 4 9= 9 9+ 1 = 9 5 + 1 = 5 （2） MAXA(1) = 1 MAXA(2) = 2 MAXA(3) = MAXA(2) + 2= 2 + 2 = 4 MAXA(4) = 4 + 3= 7 MAXA(5) = 7 + 4= 11 MAXA(6) = 11 + 5= 15 MAXA(7) = 15 + 6= 20 MAXA(8) = 20 + 7= 26 MAXA(9) = 26 + 8= 30 （3）S = MAXA(10) MAXA(11) = 30 + 5 1 = 34 习题习题三三 （第三章 线性方程组的若干直接解法） 3.1 证明若A为阶对称正定矩阵， 则A的所有对角线元素必为正数。 解：假设矩阵A为如下形式： A = 11111 21222 12nn 为一个n n的对称矩阵，即 = ， ， = 1,2,3,n 现证明如果该矩阵A为阶对称正定矩阵，则由正定矩阵的定义 出发证明其所有对角线元素为正数； 定义指出对于一个阶对称矩阵A，若存在任意一个非零列向量 Q，使得QTAQ恒大于 0，则矩阵A为一个对称正定矩阵； 假定存在一个列向量Q， 所有其余分量为 1， 所有其余分量为 0， 则计算QTAQ的值； QTAQ = 0,0,1,0T 11111 21222 12nn 0 1 0 = 1 21 0 1 0 = 0 由于对任一列向量Q，QTAQ恒大于 0，则证明A的所有对角线 元素必为正数。 3.2 证明下面的矩阵 A = 210 141 012 是一个对称正定矩阵。 证明：利用A为对称正定矩阵的充分必要条件: （A） 0 （j = 1,2,3) （A1） = | 210 141 012 | = 12 0 （A2） = | 210 141 012 | = 7 0 （A3） = 2