第四章轴对称问题有限元法

第四章 轴对称问题有限元法

在工程问题中经常会遇到一些实际结构，它们的几何形状、约束条件和外载荷均对称某一固定轴，我们把该固定轴称为对称轴。则在载荷作用下产生的应力、应变和位移也都对称此轴。这种问题就称为轴对称问题。在离心机械、压力容器、矿山机械、飞行器中经常遇到轴对称问题。

第一节 轴对称问题弹性力学基本方程

对于轴对称问题，宜采用圆柱坐标系（,,r z θ）。如果将

y

弹性体的对称轴作为Z 轴，则所有应力、应变和位移分量都只是r 和Z 轴的函数，而与θ无关，即不随θ变化。弹性体内任意一点只有两个位移：即沿r 方向的径向位移u 和沿Z 方向的轴向位移w 。由于轴对称，沿θ方向的环向（周向）位移v 等于零。因此轴对称问题是二维问题。

在轴对称弹性体内用相距dr 的两个圆柱面和过轴线互

成d θ角的两个铅垂面切割出一个高为dz 的微元体，如图2所示。

(a)

σ(b)

沿r 方向作用的正应力r σ称为径向应力 沿θ方向作用的正应力θσ称为环向应力 沿z 方向作用的正应力z σ称为轴向应力 rz 面内的剪应力 zr τ=rz τ

故轴对称弹性体内任意一点的应力分量

{}[]T

r z rz θσσσστ=

对应的轴对称弹性体内任意一点的应变分量

{}[]

T

r z rz θεεεεγ=

其中

r ε ------ 沿r 方向径向线应变 θε ------ 沿θ方向环向线应变 z ε ------ 沿z 方向轴向线应变

rz γ------ rz 面内的剪应变

与平面问题相比，轴对称问题多了一个环向应变θε。弹性体受载时，点（,,r z θ）产生径向位移u ，使过点（,,r z θ）的周长增加了2()2r u r ππ+-，因而产生相对伸长，即环向应变:

2()22r u r u r r

θππεπ+-==

轴对称问题的几何方程（应变与位移之间的关系）为

,,,r z zr u u w w u

r r z r z

θεεεγ????====+????

写成矩阵形式

{}r z rz u r u r

w u w z r θεεεεγ????

?????

???????????????????????

????

==?????+??

根据虎克定律，应力与应变的关系为

1

()r r z E

θεσμσσ????=-+ 1

()z r E

θθεσμσσ????=-+ 1

()z z r E

θεσμσσ????=-+ 12(1)rz rz rz r G E

μττ+==

由上式得

[]10111

011(1)(1)(12)10

111200

2(1)r z zr r z rz E θθσσσστμ

μ

μμ

εμ

μ

εμμμμ

μ

εμμμμγμμ??????=??

??????

?

?

??

--??

????????---??????+-????????--????-????-?

?

= （4-2） 这里弹性矩阵[D]为

[D]＝

10111

011(1)(1)(12)10111200

2(1)E μ

μ

μμ

μ

μ

μ

μμμ

μ

μμμμμμ?

?

??

??????

????

????

?????

?

-----+-----

第二节 三角形截面环单元

一、 结构离散化

离散化轴对称体时，采用的单元是一些圆环。这些圆环单元与rz 平面（子午面）正交的截面可以有不同的形状：3节点三角形、6节点三角形、4节点四边形和8节点四边形等等。单元的节点是圆周状的铰链，各单元rz 平面（子午面）内形成网格。在我们这里研究的是3节点三角形轴对称单元，这些圆环单元与rz 平面（子午面）正交的截面是三角形，如图3所示。

r

(u)

图4-3 轴对称结构

注：（1）对轴对称问题进行计算时，只需取出一个截面进行网格划分和分析，即在rz 平面（子午面）截面进行网格划分和分析。但是应注意到单元是圆环状的，所有节点载荷都

应为作用在单元节点所在的圆周上；同样，位移边界条件也是如此。

（2）轴对称体受非轴对称载荷时，成为三维问题。此时，采用将载荷沿θ方向展成富氏级数的半解析方法，把三维问题化为一组二维问题。 轴对称问题离散化例1

如图所示是一承受内压和外压的无限长厚壁圆筒，可取单位长度圆筒进行分析，有限元模型：

r

轴对称问题离散化例2

取四分之一模型研究，有限元模型（网格未划，只给出位移边界条件）：

r

二、位移模式

采用三节点三角形单元，单元节点位移列阵为：

{}T

i i j j m m e

u w u w u w δ????= （4-3）

仿照平面三角形单元，取线性位移模式

123456u r z

w r z

αααααα=++=++ （4-4）

类似平面三角形单元的推导，将节点坐标

,,,,,i i j j m m r z r z r z 和节点位移,,,,,i i j j m m u w u w u w 代入位移

模式（4-3）中，可解得126,,...,ααα。再将这些系数代回式（4-4）中，得

i i j j m m i i j j m m

u N u N u N u w N w N w N w =++=++ （4-5）

其中形函数

1

()2i i i i N a b r c z =++?

（i,j,m 轮换） (4-6)

而

11

12

1i i j

j m m

r z r z r z ?=

m m i j j

m i j m i j

a r z r z

b z z

c r r ??

???

=-=-=- （i,j,m 轮换） (4-7) 三、单元应变

为了求单元应变，将式（4-5）代入轴对称问题的几何方 程，得

{}

00000010002i i j m i r i j m j e

i j m j z rz i i j

j m m m m u u r b b b w u

f f f u r

w c c c w z c b c b c b u u w w z r θεεεεγ????

??

??

????????

?

???????

????????

????????

??????

????????????????

??????

??

??===?????+?? （4-9）

式中

i i i i a c z

f b r r

=++

（i,j,m 轮换）

式（4-9）可简写为

{}[]{}e

e

B εδ= （4-10） 其中，[]B 为变矩阵，且可写为

[]00102i i

i i i i b f B c c b ????

??=

???????

（i,j,m ） 由式（4-9）可以看出，单元应变只有环向应变θε是坐标的函数，其它应变都是常量。

四、单元应力

为求单元应力，把式（4-5）代入轴对称问题的物理方程 式，得

{}

[]{}[][]{}{}

r e

e z rz e

e

i j m D D B S S S θσσσεστδδ??

??

??

==????????

??==?? （4-11）

式中，应力矩阵的子矩阵为：

[]111131222()i i i i i

i i i i i i i b A f A c Ab f A c A S A b f c A c

A b +??

??

+??=

+??????

?

（i,j,m ） 其中

12312(1),,12(1)4(1)(12)

E

A A A μ

μμμμμμ--===

--+- 从式（4-11）可知，只有剪应力rz τ在单元中是常量，而其它应力分量在单元中都不是常量，与坐标r 和z 有关。

为了简化计算和消除由r=0引起的奇异，通常取单元形心点的r 、z 坐标值作为其近似值，即

i i i i a c z

f f b r r

≈=++ (i,j, m) 其中

11

(),()33

i j m i j m r r r r z z z z =++=++

单元网格确定后，各单元的,r z 就是定值。这样按轴对称三

节点三角形单元所求的应变和应力，是单元形心处的应变和应力近似值，它们都是常数。当单元较小时，误差很小，特别是单元离z 轴较远时,误差就更小。

第三节 单元刚度矩阵

运用虚功原理来推导轴对称问题的单元刚度矩阵。单元在节点力作用下处于平衡状态，节点力列阵为

{}

T

e

i r i z j r j z

m r

m z

R R

R R R R

R ????

= ir

假设单元e 的三个节点的虚位移为

{}

******

*

e

T i i j j m m u w u w u w δ??=??

单元上任一点的虚位移为

{}[]{}*

*e

f N δ= （4-12）

单元的虚应变为 {}[]{}*

*e

B ε

δ= （4-13）

根据虚功原理，单元体所吸收的虚应变能等于单元节点力所

做的虚功，即 {}{}{}{}*

*T e T

e

R rdrdzd δεσθ??

=???

?

??? （4-14）

将式（4-11）和式（4-13）代入式（4-14），则得

{}{}

{}

[][][]{}{}

[][][]{}*

*

*

2T T e e e

e T T T

e e R B D B rdrdzd B D B rdrdz δ

δδθ

δπδ??

??=?????

?

??

??=???

?

?????

由于虚位移列阵

{}*e

δ的任意性，所以有

{}

[][][]{}[]{}2T

e

e e

e R B D B rdrdz K πδδ==??

（4-16） 式中，[]e

K 就是单元刚度矩阵

[]

[]

[][]

2T

e

K B D B r d r d z π=?? （4-17）

写成分块形式为

[]

66

ii

ij

im e

ji jj jm mi mj

mm k k k K k k k k k k ???

??

=??????

（4-18） 其中每个子矩阵为

[]

[][]2T

s t s t

K B D B r d r d z

π??=????（s,t=i,j,m ） (4-19) 把上式中的坐标用单元形心点的坐标代替，有

[]112123122()()()2()s t t s t t s t

t s s s t st

s t t s t

s t s t

b b A f f f A b A

c c A c b f A c b rA K A c b f A b c c c A b b π++++++=

?

+++?????

?

（s,t=i,j,m ） (4-20)

求得单元刚度矩阵后，就可向平面问题一样“子块搬家，对号入座”组集整体刚度矩阵。

注：实践表明，只要网格剖分不太粗，这样近似计算引起的误差是很小的。

第四节 等效节点载荷计算

与平面问题类似，当结构外载荷不作用在节点上时，也需要将这些作用在环形单元上的集中力、表面力、和体积力分别移置到节点上。移置的原则也是要求这些外力和等效节点载荷在任意位移上所作的虚功相等，即

{}()

{}{}{}{}{}{}{}

****2T e

T T T

e c

f p rd f q rd f G R drdz ds r δ

θθπ=++?????

其中：

{}

*******e

T

u u u m m i i j j

w w w δ=???

?

且单元内任意一点的位移与节点位移之间有

{}[]{}e u N v f δ=??

????

=

故

{}

[]{}

****

e u N w

f δ=??????????

=

所以有：

{}(){}{}()[]{}[]{}[]{}*

*

222T

T

T

c l T e e

T e drdz R

N p r N q rds r N G πππδδ?=

++?

?

??

?

??? （4-27） 注意

{}*

e

δ的任意性，所以有

{}

[]{}[]{}[]{}

222T T

e

T

c

l

R drdz r N p r N q rds N G π

ππ?

=++??? （4-28） 1. 体积力 （1）自重

假设单元体力分量为

{}[]0

T

p ρ=-

式中ρ为材料重量密度。由式（4-28），节点i 的等效节点载

荷为

{}02ir i i i z R R N rdrdz R π

ρ?

????

??????

??

==-?? 和平面问题一样，用面积坐标建立如下关系

,

i i i i j j m m

L N r L r L r L r ==++ 利用面积坐标积分，则得

()(3)

6121212i j i m i m m i i j j N rdrdz Li drdz r r r r r L r L r L r ??

?? ???

=?

=++?=+++????

{}0(3)6i r i i z i r r R R R πρ??????

??????????

?-+== （i ，j ，m 轮换）

{}

0303036i j m e

r r r r r r R πρ????+????????+??????+????

?=-

（2）离心力 设

{}20r p ρω??

????

=

式中ρ为材料质量密度，ω为结构绕对称轴旋转的角速度。 节点i 的等效节点载荷为 {}220e

i i g r R N rdrdz ρωπ?

??

????=??

积分

222

22()(6)30

i j m i j m m m i i i j j N r drdz Li r L r L r L drdz r r r rr r r ?

?

=++?=++++????

所以节点i 的等效节点载荷为

{}

22

(92)150i j m e

i r r r r R πρω???

???????

?+-= （i ，j ，m 轮换）

{}

222222292092150920m m i j e

m j

i i j

r r r r r r r r R r r rr πρω????

??

????????

??

??

????

+-+-?=+-

2. 表面力

设轴对称问题三角形截面环形单元的ij 边上作用有线性分布的r 方向面力。面力在节点i 处的集度为i q ，在节点j 处的集度为j q ，i j 边长为l 。

则在i j 边上，表面力为{}0i i j j q N q N q ??

????

+=

节点i 的等效节点载荷为

{}

02i i j

j

i l

e

i q N q N N R r d s

π???

???

+=? 因在i j 边上，面积坐标0m

L =，于是

22(3)12()i j i

i

i i j j l l l

r r N rds L r L r L ds +=+=?? 22()12()i j j

j

i i j j l l

l

r r N rds L r L r L ds +=+=?? 则

{}

{

}

6

(3)()i i j j i j e

i l

q r r q r r R π+++=

同理

{}{}6()(3)i

i

j

j

i

j

e

j

q r r q r r l

R π+++=

{}0m R =

单元等效节点载荷为

{}

(3)()0()(3)0600i i j j i j e

i i j j i j q r r q r r q r r q r r l R π????????????

????

??

++++++=

当各种力作用下的节点载荷列阵计算后，与平面问题类