4.弹性力学轴对称问题的有限元法 本章包括以下内容: 4.1用虚功方程建立有限元方程 42三结点单元位移函数 43三结点单元刚度矩阵 4.4载荷移置 45轴对称分析举例 4.1用虚功方程建立有限元方程 物体的几何形状、约束情况及所受的外力都对称于空间的某一根轴,因此在物体中通过 该轴的任何平面都是对称面,所有应力、应变和位移也对称于该轴,这类问题称为轴对称问 题。研究轴对称问题时通常采用圆柱坐标系(r,θ,z),以z轴为对称轴 图41受均布内压作用的长圆筒 如图41所示的受均布内压作用的长圆筒,通过Z轴的一个纵截面就是对称面。由于 对称性,轴对问题共有4个应力分量 (4-1) 其中σ,表示沿半径方向的正应力,称为径向应力:O表示沿0方向的正应力,称为 环向应力或切向应力;:表示沿z方向的正应力,称为轴向应力;rx表示在圆柱面上沿 z方向作用的剪应力 同样,轴对称问题共有4个应变分量:
4. 弹性力学轴对称问题的有限元法 本章包括以下内容： 4.1 用虚功方程建立有限元方程 4.2 三结点单元位移函数 4.3 三结点单元刚度矩阵 4.4 载荷移置 4.5 轴对称分析举例 4.1 用虚功方程建立有限元方程 物体的几何形状、约束情况及所受的外力都对称于空间的某一根轴，因此在物体中通过 该轴的任何平面都是对称面，所有应力、应变和位移也对称于该轴，这类问题称为轴对称问 题。研究轴对称问题时通常采用圆柱坐标系（r，θ，z），以 z 轴为对称轴。 图 4.1 受均布内压作用的长圆筒 如图 4.1 所示的受均布内压作用的长圆筒，通过 Z 轴的一个纵截面就是对称面。由于 对称性，轴对问题共有 4 个应力分量： = zr z r { } （4-1） 其中 r 表示沿半径方向的正应力，称为径向应力； 表示沿θ方向的正应力，称为 环向应力或切向应力； z 表示沿 z 方向的正应力，称为轴向应力； zr 表示在圆柱面上沿 z 方向作用的剪应力。 同样，轴对称问题共有 4 个应变分量：
E (4-2) 其中E,表示沿半径方向的正应变,称为径向正应变:Ea表示沿方向的正应变,称为 环向正应变或切向正应变;E=表示沿z方向的正应变,称为轴向正应变;yx表示沿r和z 方向的剪应变。 在轴对称问题中,弹性体内任意一点上,不存在切向位移,只存在径向位移u和轴向位 移w,两个位移分量表示为, f (4-3) 在讨论弹性力学平面问题的有限元法时,我们先由将弹性体划分为有限个单元的组合 体,由虚功方程得到单元刚度矩阵,集成后得到整体刚度矩阵。在这里,我们用虚功方程直 接得到轴对称问题的有限元列式 由虚功方程可得,外力虚功等于内力虚功或虚应变能 (4-4) 其中{F}为体力,{p}为面力 将弹性体离散后,作用在弹性体上的外载荷移置到节点上,在每个节点上外力只有径向 分量U1,U2…,Un,轴向分量W,W2…,Wn, {F}={W2 (4-5) 每个节点的虚位移也只有径向分量a,12…,n,轴向位移分量w1,w2…,n
= zr z r { } （4-2） 其中 r 表示沿半径方向的正应变，称为径向正应变； 表示沿θ方向的正应变，称为 环向正应变或切向正应变； z 表示沿 z 方向的正应变，称为轴向正应变； zr 表示沿 r 和 z 方向的剪应变。 在轴对称问题中，弹性体内任意一点上，不存在切向位移，只存在径向位移 u 和轴向位 移 w，两个位移分量表示为， = w u { f } （4-3） 在讨论弹性力学平面问题的有限元法时，我们先由将弹性体划分为有限个单元的组合 体，由虚功方程得到单元刚度矩阵，集成后得到整体刚度矩阵。在这里，我们用虚功方程直 接得到轴对称问题的有限元列式。 由虚功方程可得，外力虚功等于内力虚功或虚应变能， dxdydz f F dxdydz f p ds T s T T { } { } { } { } { } { } * * * = + （4-4） 其中{F}为体力，{p}为面力。 将弹性体离散后，作用在弹性体上的外载荷移置到节点上，在每个节点上外力只有径向 分量 U U Un , ,..., 1 2 ，轴向分量 W W Wn , ,..., 1 2 ， = n n W U W U W U F ... { } 2 2 1 1 （4-5） 每个节点的虚位移也只有径向分量 * * 2 * 1 , ,..., u u un ，轴向位移分量 * * 2 * 1 , ,..., w w wn
6}=1w2 (4-6) 在单元中由虚位移引起的虚应变为, {E"}=[B]{G6"}° (4-7) 单元中的实际应力为 {G}°=[DB]{G}2 (4-8) 离散后的单元组合体的虚功方程为, }{F}=∑(6}°) B DIBS dxdydz (4-9) (F)=2(6))'J. [B][DILB]dxdyd-48 (4-10 K=.BDJB就是单元刚度矩阵 对于轴对称问题 [K]° [B] [DIB]rdrdEd0=27 [B] [DIB]rdrd= (4-11) 将(4-11)代入(4-10)可得 6}{F}=('}∑(GKrG)o (4-12) K]=∑(Kr(G)为整体刚度矩阵,得到方程组, [K]{G}={F (4-13)
= * * * 2 * 2 * 1 * 1 * ... { } n n w u w u w u （4-6） 在单元中由虚位移引起的虚应变为， e e { } [B]{ } * * = （4-7） 单元中的实际应力为， e e {} = [D][B]{} （4-8） 离散后的单元组合体的虚功方程为， = = n i e T T e e T F B D B dxdydz 1 * * { } { } ({ } ) [ ] [ ][ ]{} （4-9） e n i T e T e T { } {F} ({ } ) [B] [D][B]dxdydz{ } 1 * * = = （4-10） K B D B dxdydz T e e [ ] [ ] [ ][ ] = 就是单元刚度矩阵。 对于轴对称问题， K B D B rdrdzd B D B rdrdz e T T [ ] [ ] [ ][ ] 2 [ ] [ ][ ] 2 0 = = （4-11） 将（4-11）代入（4-10）可得 { } { } { } ([ ] [ ] [ ]){ } * * = e T T T e F G K G （4-12） = e T e [K] ([G] [K] ([G]) 为整体刚度矩阵，得到方程组， [K]{} = {F} （4-13）
42三结点单元位移函数 轴对称问题分析中所使用的三结点单元,在对称面上是三角形,在整个弹性体中是三棱 圆环,各单元中圆环形铰相联接。参照平面问题的三角形单元位移函数,轴对称问题的三结 点三角形单元位移函数取为 u=a1+ar+a3 (4-14) M (m, zm) 2,21) 图4-2三结点单元 按照平面问题三角形单元的分析过程,将结点坐标和结点位移代入(4-14)得到, b, b, bm ru (4-15) (4-16) C 其中 (4-17) b 定义形态函数为, N=立(a1+b+c12)(下标im轮换)(4-18) 用矩阵表示的单元位移为
4.2 三结点单元位移函数 轴对称问题分析中所使用的三结点单元，在对称面上是三角形，在整个弹性体中是三棱 圆环，各单元中圆环形铰相联接。参照平面问题的三角形单元位移函数，轴对称问题的三结 点三角形单元位移函数取为， = + + = + + a a r a z u a a r a z 4 5 6 1 2 3 w （4-14） 图 4-2 三结点单元 按照平面问题三角形单元的分析过程，将结点坐标和结点位移代入（4-14）得到， = m j i i j m i j m i j m u u u c c c b b b a a a A a a a 2 1 3 2 1 （4-15） = m j i i j m i j m i j m w w w c c c b b b a a a A a a a 2 1 6 5 4 （4-16） 其中， m m j j i i r z r z r z A 1 1 1 2 1 = （4-17） i j m m j a = r z − z r ，bi = zj − zm ， i m j c = r − r 定义形态函数为， ( ) 2 1 a b r c z A Ni = i + i + i (下标 i,j,m 轮换) (4-18) 用矩阵表示的单元位移为
uli 1「N0N,0N Ni Ni0 N 43三结点单元刚度矩阵 轴对称问题的几何方程 E E (4-20) 02 au a 由(4-19)式得 24ov+ (4-2la) (4-2lb) r 2A 其中,∫=二+b2+二(下标轮换) (c1;+c (4-21c) 2-24c1l+c1+cm4m) (4-21d) ar 2A (b; wi +b (4-2le) 用几何矩阵表示单元的应变 {}=[B]{o} (4-22) B=BBB. (423)
= = m m j j i i i j m i j m w u w u w u N N N N N N w u f 0 0 0 0 0 0 (4-19) 4.3 三结点单元刚度矩阵 轴对称问题的几何方程： + = r w z u z w r u r u zr z r （4-20） 由（4-19）式得， ( ) 2 1 biui b ju j bmum r A u = + + （4-21a） ( ) 2 1 f iui f ju j fmum r A u = + + （4-21b） 其中， r cz b r a f i i i i = + + （下标轮换） ( ) 2 1 ciwi c jwj cmwm z A w = + + （4-21c） ( ) 2 1 ciui c ju j cmum z A u = + + （4-21d） ( ) 2 1 biwi b jwj bmwm r A w = + + （4-21e） 用几何矩阵表示单元的应变， e {} = [B] {} （4-22） [ ] [ ] B = Bi Bj Bm （4-23）
6. 0 f10 [B;]= (下标轮换) (424) 2A0 Ci 由于在f是坐标r、z的函数,E分量在单元中不为常量,其它三个应变分量在单元 中仍为常量。 由轴对称问题的物理方程,得到弹性矩阵 1-1- E(1-4)1 (1+)(1-21) 0 0 0 2(1-) 令 2(1 )42,则弹性矩阵为, 1A1A10 [D]= E(-)4140 (4-26) (1+)(1-2)4141 0 000A 由弹性矩阵[D]和几何矩阵[B可以得到应力矩阵[S],并计算出单元内的应力分量, [S]=[D[B (4-27) S=Si Si smI b+f141 S;]=[DIB]= E(1-)4b+f (4-27) 2(1+1)(1-2u)4A1(b+f)c A2c A2b 下标轮换,可得到[S;1[S 由应力矩阵可知,除剪应力rx为常量,其它三个正应力分量都是r、z的函数
= i i i i i i c b c f b A B 0 0 0 2 1 [ ] （下标轮换） （4-24） 由于在 i f 是坐标 r、z 的函数， 分量在单元中不为常量，其它三个应变分量在单元 中仍为常量。 由轴对称问题的物理方程，得到弹性矩阵， − − − − − − − − + − − = 2(1 ) 1 2 0 0 0 1 0 1 1 0 1 1 1 0 1 1 1 (1 )(1 2 ) (1 ) [ ] E D （4-25） 令 1 1 = A − ， 2 2(1 ) 1 2 = A − − ，则弹性矩阵为， + − − = 2 1 1 1 1 1 1 0 0 0 1 0 1 0 1 0 (1 )(1 2 ) (1 ) [ ] A A A A A A A E D （4-26） 由弹性矩阵[D]和几何矩阵[B]可以得到应力矩阵[S]，并计算出单元内的应力分量， [S] = [D][B] （4-27） （4-28） [ ] [ m] S = Si S j S + + + + − − = = i i i i i i i i i i i i i A c A b A b f c A b f A c b f A c A E S D B 2 2 1 1 1 1 2(1 )(1 2 ) ( ) (1 ) [ ] [ ][ ] （4-27） 下标轮换，可得到 [ ], [ ] Si S j 。 由应力矩阵可知，除剪应力 zr 为常量，其它三个正应力分量都是 r、z 的函数
单元刚度矩阵为 kr=2』[B[ DIB]rrde (4-28) 单元刚度矩阵的分块矩阵为, KK]=2[ [B] [D][B],rdrd-s (4-29) 由于几何矩阵中的元素不是常量,单元刚度矩阵需要通过积分得到,为简化计算,可以 用三角形单元形心位置的坐标r,=c代替[B矩阵中的变量r、z r +r+r 应变矩阵变成 B]=园 b b;+c1;c0 B (4-30) b 单元刚度矩阵的近似表达式为 [K=2mB]ID B 单元刚度矩阵的分块矩阵近似表达式为, 区]=2mBD[B [K] E(1-)「b+,,+A(b,+fb,)+Ac 2(1+p)(1-21)4 A(c, b,+crfs)+A,b,c (4-33) A,(b, cs +f, c,)+A,c, b C,+a,bb 44载荷移置 单元上的体力为{p},与平面问题相同,由虚功方程可以得到结点载荷, (R =2T[INI(p)rdrd= (4-34) 作用在单元上的面力为何},结点载荷为
单元刚度矩阵为， K B D B rdrdz e T [ ] 2 [ ] [ ][ ] = （4-28） 单元刚度矩阵的分块矩阵为， K B D B rdrdz s T rs r 2 [ ] [ ] [ ] = （4-29） 由于几何矩阵中的元素不是常量，单元刚度矩阵需要通过积分得到，为简化计算，可以 用三角形单元形心位置的坐标 c c r ,z 代替[B]矩阵中的变量 r、z。 ( ) 3 1 c i j m r = r + r + r ， ( ) 3 1 zc = zi + z j + zm 应变矩阵变成， m [B] = Bi Bj B + + = i i i i i c c i i i c b c b c z r a b A B 0 0 0 2 1 [ ] （4-30） 单元刚度矩阵的近似表达式为： [K] 2 r [B] [D] [B] T c e = （4-31） 单元刚度矩阵的分块矩阵近似表达式为， s T K rs = 2rc [B] r [D] [B] （4-32） + + + + + + + + + + − − = r s r s r s r s r s r s r s r s c r s r s r s r s r s r s c c A b b A b c f c A c b A c b c f A b c b b f f A b f f b A c c A E r K 2 1 2 1 2 1 2 ( ) ( ) ( ) 2(1 )(1 2 ) (1 ) [ ] (4-33) 4.4 载荷移置 单元上的体力为{p}，与平面问题相同，由虚功方程可以得到结点载荷， R = N p rdrdz e T { } 2 [ ] { } (4-34) 作用在单元上的面力为 P ，结点载荷为
(R; =27. INT(P)rds (4-35) 轴对称问题分析中,如果直接定义结点载荷,载荷值是实际弹性体上绕对称轴一周的载 荷的累计结果。 45轴对称分析实例 图4- 座封头的结构 图44坯料形状 图4-5成形分析的轴对称有限元模型 封头作为压力容器中的重要受力部件,用户对其质量、强度、安全性等有很高的要求 带裙座封头的结构如图4-3所示,其优点是可以避免直接在封头壁上进行焊接,提高了封头 的可靠性,但也增加了成形过程的难度。成形的难点在于: 1)如何保证锻件的厚度
= s e T {R} 2 [N] {P}rds （4-35） 轴对称问题分析中，如果直接定义结点载荷，载荷值是实际弹性体上绕对称轴一周的载 荷的累计结果。 4.5 轴对称分析实例 图 4-3 带裙座封头的结构 图 4-4 坯料形状 图 4-5 成形分析的轴对称有限元模型 封头作为压力容器中的重要受力部件，用户对其质量、强度、安全性等有很高的要求。 带裙座封头的结构如图 4-3 所示，其优点是可以避免直接在封头壁上进行焊接，提高了封头 的可靠性，但也增加了成形过程的难度。成形的难点在于： 1） 如何保证锻件的厚度；
2)如何保证成形后的裙座位置。 厚壁封头在热冲压成形过程中还会出现明显的局部减薄或增厚现象,严重的会导致封头 撕裂、起皱、模具涨裂等问题。 制造带裙座封头关键之一是如何设计出一个特殊形状的坯料。普通的半球形封头采用圆 饼形坯料,制造带裙座封头要采用如图4-4所示的坯料。 分析整个成形过程可以发现,封头的底部明显变薄,会使封头的最小壁厚达不到设计要 求。在制作坯料时,要在坯料的中心部分加厚。封头边缘部分,在成形过程中明显增厚,壁 厚的增加量会超过10%,制作坯料时要在坯料的边缘部分减薄。在图4-5中,可以看出,我 们制作了一个心部增厚,边缘减薄的坯料 坯料上预制的凸台位置与成形后的裙座位置密切相关,由于成形过程中封头的底部变薄 导致凸台外移,合理的凸台位置要通过有限元分析来选择 TIME2.75 图4-6成形初期的等效应力分布 图47成形中间阶段的等效应力分布
2）如何保证成形后的裙座位置。 厚壁封头在热冲压成形过程中还会出现明显的局部减薄或增厚现象，严重的会导致封头 撕裂、起皱、模具涨裂等问题。 制造带裙座封头关键之一是如何设计出一个特殊形状的坯料。普通的半球形封头采用圆 饼形坯料，制造带裙座封头要采用如图 4-4 所示的坯料。 分析整个成形过程可以发现，封头的底部明显变薄，会使封头的最小壁厚达不到设计要 求。在制作坯料时，要在坯料的中心部分加厚。封头边缘部分，在成形过程中明显增厚，壁 厚的增加量会超过 10%，制作坯料时要在坯料的边缘部分减薄。在图 4-5 中，可以看出，我 们制作了一个心部增厚，边缘减薄的坯料。 坯料上预制的凸台位置与成形后的裙座位置密切相关，由于成形过程中封头的底部变薄 导致凸台外移，合理的凸台位置要通过有限元分析来选择。 图 4-6 成形初期的等效应力分布 图 4-7 成形中间阶段的等效应力分布
图48成形结束阶段的等效应力分布 图4-9等效应变分布与成形缺陷 通过有限元分析还发现,如果坯料上的凸台尺寸过大,会在封头的内壁上产生图49所 示的凹限,导致封头内表面尺寸超出设计要求。 采用 ANSYS软件,对坯料形状和尺寸、模具的尺寸、成形缺陷进行了综合分析得到了 优化的坯料设计和制造工艺
图 4-8 成形结束阶段的等效应力分布 图 4-9 等效应变分布与成形缺陷 通过有限元分析还发现，如果坯料上的凸台尺寸过大，会在封头的内壁上产生图 4-9 所 示的凹限，导致封头内表面尺寸超出设计要求。 采用 ANSYS 软件，对坯料形状和尺寸、模具的尺寸、成形缺陷进行了综合分析得到了 优化的坯料设计和制造工艺