3 弹性力学平面问题的有限元法
本章包括以下的内容:
3.1弹性力学平面问题的基本方程
3.2单元位移函数
3.3单元载荷移置
3.4单元刚度矩阵
3.5单元刚度矩阵的性质与物理意义
3.6整体分析
3.7约束条件的处理
3.8整体刚度矩阵的特点与存储方法
3.9方程组解法
3.1弹性力学平面问题的基本方程
弹性力学是研究弹性体在约束和外载荷作用下应力和变形分布规律的一门学科。在弹性力学中针对微小的单元体建立基本方程,把复杂形状弹性体的受力和变形分析问题归结为偏微分方程组的边值问题。弹性力学的基本方程包括平衡方程、几何方程、物理方程。
弹性力学的基本假定如下:
1)完全弹性,2)连续,3)均匀,4)各向同性,5)小变形。
3.1.1基本变量
弹性力学中的基本变量为体力、面力、应力、位移、应变,各自的定义如下。
体力
体力是分布在物体体积内的力,例如重力和惯性力。
面力
面力是分布在物体表面上的力,例如接触压力、流体压力。
应力
物体受到约束和外力作用,其内部将产生内力。物体内某一点的内力就是应力。
图3.1
如图3.1假想用通过物体内任意一点p 的一个截面mn 将物理分为Ⅰ、Ⅱ两部分。将部分Ⅱ撇开,根据力的平衡原则,部分Ⅱ将在截面mn 上作用一定的内力。在mn 截面上取包含p 点的微小面积A ?,作用于A ?面积上的内力为Q ?。
令A ?无限减小而趋于p 点时,Q ?的极限S 就是物体在p 点的应力。
S A Q
A =??→?0lim
应力S 在其作用截面上的法向分量称为正应力,用σ表示;在作用截面上的切向分量称为剪应力,用τ表示。
显然,点p 在不同截面上的应力是不同的。为分析点p 的应力状态,即通过p 点的各个截面上的应力的大小和方向,在p 点取出的一个平行六面体,六面体的各楞边平行于坐标轴。
图3.2
将每个上的应力分解为一个正应力和两个剪应力,分别与三个坐标轴平行。用六面体表面的应力分量来表示p 点的应力状态。应力分量的下标约定如下:
第一个下标表示应力的作用面,第二个下标表示应力的作用方向。
xy τ,第一个下标x 表示剪应力作用在垂直于X 轴的面上,第二个下标y 表示剪应力指
向Y 轴方向。
正应力由于作用表面与作用方向垂直,用一个下标。x σ表示正应力作用于垂直于X 轴的面上,指向X 轴方向。
应力分量的方向定义如下:
如果某截面上的外法线是沿坐标轴的正方向,这个截面上的应力分量以沿坐标轴正方向为正;
如果某截面上的外法线是沿坐标轴的负方向,这个截面上的应力分量以沿坐标轴负方向为正。
剪应力互等:xz zx zy yz yx xy ττττττ===,,
物体内任意一点的应力状态可以用六个独立的应力分量x σ、y σ、z σ、xy τ、yz τ、zx
τ
来表示。 位移
位移就是位置的移动。物体内任意一点的位移,用位移在x ,y ,z 坐标轴上的投影u 、v 、w 表示。
应变
物体的形状改变可以归结为长度和角度的改变。 各线段的单位长度的伸缩,称为正应变,用ε表示。
两个垂直线段之间的直角的改变,用弧度表示,称为剪应变,用γ表示。
物体内任意一点的变形,可以用zx yz xy z y x γγγεεε、、、、、六个应变分量表示。
3.1.2平面应力和平面应变问题
弹性体在满足一定条件时,其变形和应力的分布规律可以用在某一平面内的变形和应力的分布规律来代替,这类问题称为平面问题。平面问题分为平面应力问题和平面应变问题。 1)平面应力问题
设有很薄的等厚薄板,只在板边上受到平行于板面并且不沿厚度变化的面力,体力也平行于板面且不沿厚度变化。
图3.3
设板的厚度为t ,在板面上: ()02
=±=t
z z σ, ()02
=±=t z zx τ, ()02
=±=t z zy τ
由于平板很薄,外力不沿厚度变化,因此在整块板上有,
0=z σ, 0=zx τ, 0=zy τ
剩下平行于XY 平面的三个应力分量xy y x τσσ、、未知。
2)平面应变问题
设有很长的柱形体,支承情况不沿长度变化,在柱面上受到平行于横截面而且不沿长度变化的面力,体力也如此分布。
图3.4
以柱体的任一横截面为XY 平面,任一纵线为Z 轴。假定该柱体为无限长,则任一截面都可以看作对称面。由对称性, 0=zx τ,0=zy τ,0=w
由于没有Z 方向的位移,Z 方向的应变0=z ε。
未知量为平行于XY 平面的三个应力分量xy y x τσσ、、,物体在Z 方向处于自平衡状
态。
3.1.3平衡方程
弹性力学中,在物体中取出一个微小单元体建立平衡方程。平衡方程代表了力的平衡关系,建立了应力分量和体力分量之间的关系。对于平面问题,在物体内的任意一点有,
00=+??+??=+??+??Y x
y X y
x xy
y yx
x τστσ
(3-1)
3.1.4几何方程
由几何方程可以得到位移和变形之间的关系。对于平面问题,在物体内的任意一点有,
x
v y u y v x u
xy y x ??+
??=??=??=
γεε
(3-2)
刚体位移
由位移u=0,v=0可以得到应变分量为零,反过来,应变分量为零则位移分量不为零。应变分量为零时的位移称为刚体位移。刚体位移代表了物体在平面内的移动和转动。
由
000=??+??=??=??x
v y u y v
x u 可以得到刚体位移为以下形式,
x
v v y
u u ωω+=-=00
由
0,0=??=??y
v
x u 可得,
)(),
(21x f v y f u ==
将21,f f 代入
0=??+??x
v
y u 可得, ω==-
dx
x df dy y df )
()(21 积分后得到,
x
v x f y u y f ωω+=-=0201)()(
得到位移分量,
x
v v y
u u ωω+=-=00
当0,0,000==≠ωv u 时,物体内任意一点都沿x 方向移动相同的距离,可见0u 代
表物体在x 方向上的刚体平移。
当0,0,000=≠=ωv u 时,物体内任意一点都沿y 方向移动相同的距离,可见0v 代表
物体在y 方向上的刚体平移。
当0,0,000≠==ωv u 时,可以假定0>ω,此时的物体内任意一点P (x ,y )的位
移分量为, x v y u ωω=-=,
P 点位移与y 轴的夹角为α,
θωωαtg x
y
x y tg ===
P 点合成位移为,
r y x x y v u ωωωω=+=+-=+222222)()(
r 为P 点到原点的距离,可见ω代表物体绕z 轴的刚体转动。
3.1.5物理方程
弹性力学平面问题的物理方程由广义虎克定律得到。 1)平面应力问题的物理方程
()y x x E μσσε-=
1
()x y y E
μσσε-=1
(3-3)
xy xy E
τμγ)1(2+=
平面应力问题有,
0=z σ
()y x
z E
σσ
μ
ε+-
=
2)平面应变问题的物理方程
???? ??---=
y x x E σμμσμε112
???
? ??---=
x y y E
σμμσμε112
(3-4)
xy xy E
τμγ)1(2+=
平面应变问题有,
0=z ε
()y x z σσμσ+=
在平面应力问题的物理方程中,将E 替换为
2
1μ-E
、μ替换为μμ-1,可以得到平面
应变问题的物理方程;在平面应变问题的物理方程中,将E 替换为
2
)
1()
21(μμ++E 、μ替换为μ
μ
+1,可以得到平面应力问题的物理方程。
图3.5
求解弹性力学平面问题,可以归结为在任意形状的平面区域Ω内已知控制方程、在位移边界u S 上约束已知、在应力边界σS 上受力条件已知的边值问题。然后以应力分量为基本未知量求解,或以位移作为基本未知量求解。
如果以位移作为未知量求解,求出位移后,由几何方程可以计算出应变分量,得到物体的变形情况;再由物理方程计算出应力分量,得到物体的内力分布,就完成了对弹性力学平面问题的分析。
3.2单元位移函数
根据有限元法的基本思路,将弹性体离散成有限个单元体的组合,以结点的位移作为未知量。弹性体内实际的位移分布可以用单元内的位移分布函数来分块近似地表示。在单元内的位移变化可以假定一个函数来表示,这个函数称为单元位移函数、或单元位移模式。
对于弹性力学平面问题,单元位移函数可以用多项式表示,
...26524321++++++=y a xy a x a y a x a a u ...26524321++++++=y b xy b x b y b x b b v
(3-5)
多项式中包含的项数越多,就越接近实际的位移分布,越精确。具体取多项,由单元形式来确定。即以结点位移来确定位移函数中的待定系数。
图3.6
如图3.6所示的3结点三角形单元,结点I 、J 、M 的坐标分别为),(i i y x 、),(j j y x 、
),(m m y x ,结点位移分别为i u 、i v 、j u 、j v 、m u 、m v 。六个节点位移只能确定六个多
项式的系数,所以3结点三角形单元的位移函数如下,
?
??
++=++=y a x a a y a x a a u 654321v
(3-6)
将3个结点上的坐标和位移分量代入公式(3-6)就可以将六个待定系数用结点坐标和位移分量表示出来。
将水平位移分量和结点坐标代入(3-6)中的第一式,
m
m m j j j i i i y a x a a u y a x a a u y a x a a u 321321321++=++=++= 写成矩阵形式,
??
???
?????????????
??=??????????321111a a a y x
y x y x u u u m m j j i i m j i
(3-7)
令
[]T 111=????
?
?????m m j j i i y x y x
y x , 则有
[]??
???
?????=??????????-m j i u u u a a a 1321T (3-8)
T
]T []T [*
1
=-
A 2T =,A 为三角形单元的面积。
[T]的伴随矩阵为,
[]T
*T ?????
????
?---------=i j j
i i j j i m i i m m
i i m j m m j j
m m j x x y y y x y x x x y y y x y x x x y y y x y x (3-9)
令 ????
????
?
?=??
???
?????=m j
i m j i m j i m m
m
j j j
i i i
c c c b b b a a a c b a c b a c b a T
*
]T [
(3-10)
则 ??
?????????????????
?=?????
?????m j i m j
i
m j i
m j i u u u c c c b b b a a a A a a a 21321
(3-11)
同样,将垂直位移分量与结点坐标代入公式(3-6)中的第二式,可得,
??
???
????????????
?
??=?????
?????m j i m j
i
m j i m j i v v v c c c b b b a a a A a a a 21654
(3-12)
将(3-11)、(3-12)代回(3-6)整理后可得,
])()()[(21
m m m m j j j j i i i i u y c x b a u y c x b a u y c x b a A u ++++++++=
])()()[(21m m m m j j j j i i i i v y c x b a v y c x b a v y c x b a A
v ++++++++=
令)(21
y c x b a A
N i i i i ++=
(下标i ,j ,m 轮换)
可得
????
?
???????????????????????=??????m m j j i i m j
i
m j i v u v u v u N N N N N N v u 0
000 (3-13)
单元内的位移记为
{}?
??
??
?=v u f 单元的结点位移记为
{}???
??
?????????????????=??????????=m m j j i
i m j i e
v u v u v u δδδδ 单元内的位移函数可以简写成,
{}[]{}e N f δ=
(3-14)
把[N]称为形态矩阵,N i 称为形态函数。
选择单元位移函数应满足以下条件: 1)反映单元的刚体位移与常量应变。
2)相邻单元在公共边界上的位移连续,即单元之间不能重叠,也不能脱离。
由(3-6)可以将单元位移表示成以下的形式,
y a a y a a x a a u 223
53521++--
+= x a a x a a y a a v 2
23
53564++-+
+= 反映了刚体位移和常应变。
单元位移函数是线性插值函数,因此单元边界上各点的位移可以由两个结点的位移完全确定。两个单元的边界共用两个结点,所以边界上的位移连续。
形态函数N i 具有以下性质:
1)在单元结点上形态函数的值为1或为0。
2)在单元中的任意一点上,三个形态函数之和等于1。
用T 来计算三角形面积时,要注意单元结点的排列顺序,当三个结点i ,j ,m 取逆时
针顺序时,0T 21>=
A ;当三个结点i ,j ,m 取顺时针顺序时,0T 2
1
<=A 。
例题:如图3.7所示等腰三角形单元,求其形态矩阵[N]。
解: 由 j m m j i y x y x a -= m j i y y b -=
j m i x x c -=
在公式中轮换下标可以计算得
0000=?-?=-=a y x y x a j m m j i ,a a y y b m j i =-=-=0, 00-=-=j m i x x c
0000=?-?=-=a y x y x a m i i m j ,000=-=-=i m j y y b , a a x x c m i j =-=-=0
200a a a y x y x a i j j i m =?-?=-=,a a y y b j i m -=-=-=0 a a x x c i j m -=-=-=0
三角形积为2
2a A =
形态函数为
a x ax a y c x
b a A N i i i i =++=++=
)00(1)(212 a y ay a
y c x b a A N j j j j =++=++=)00(1)(212
a y a x ay ax a a
y c x b a A N m m m m --=--=++=
1)(1)(2122 形态矩阵为
????
?????
?----=a y a x a
y a
x a y a x a y a
x N 10
00
100][
三角形面积的计算公式可得,
221001010
121112
1
a a a y x y x y x A m m j j i
i
===
如果把三个结点按顺时针方向排列,即i (a ,0),j (0,0),m (0,a )
221010010
121112
1
a a
a y x y x y x A m m j j i
i
-===
3.3单元载荷移置
有限元法的求解对象是单元的组合体,因此作用在弹性体上的外力,需要移置到相应的结点上成为结点载荷。载荷移置要满足静力等效原则。静力等效是指原载荷与结点载荷在任意虚位移上做的虚功相等。 单元的虚位移可以用结点的虚位移{}e
*δ表示为,
{}{}e
N f **
][δ=
(3-15)
令结点载荷为{}
?????
?????????????????=m m j j i i e
Y X Y X Y X R 1)集中力的移置
如图3.7所示,在单元内任意一点作用集中力{}?
??
???=y x P P P
图3.8
由虚功相等可得,
()()}{][}{}{}{**P N R T
T
e e
T
e δδ=
由于虚位移是任意的,则 }{][}{P N R T
e =
(3-16)
例题1:在均质、等厚的三角形单元ijm 的任意一点p (x p ,y p )上作用有集中载荷。
????????????
?
???
??????????=??????????????????????y x y x m m j j i i m m j j i i P P N N N N N N Y X Y X Y X p p )
,(0
00000
2)体力的移置
令单元所受的均匀分布体力为?
?????=y x p }{ 由虚功相等可得,
()()
??=tdxdy p N R T T
e
T
e }{][}{}
{}{**δδ
??=tdxdy p N R T e }{][}{
(3-17)
3)分布面力的移置
设在单元的边上分布有面力{}
T Y X P ],[=,同样可以得到结点载荷,
?=s T e tds P N R }{][}{
(3-18)
例题2:设有均质、等厚的三角形单元ijm ,受到沿y 方向的重力载荷q y 的作用。求均布体力移置到各结点的载荷。
?????
??????????
??
?
?????
?????=??
?
?
??
????????????????tdxdy q N N N N N N Y X Y X Y X y m m j j i i m m j j i i 0000000 0,0,0===m j i X X X
????==dxdy N t q tdxdy q N Y i y y i i
A y c x b a A A Ay c Ax b A a A dxdy y c x b a A
dxdy N c i c i i c i c i i i i i i 3
1)(21][21)(21
=++=++=++=??
??
At q Y y i 3
1
=
同理,At q Y At q Y y m y j 31
,31==
例题3:在均质、等厚的三角形单元ijm 的ij 边上作用有沿x 方向按三角形分布的载荷,求移置后的结点载荷。
????
??????????
??
?
??????????=??
?
???
????????????????s x m m j j i i m m j j i i tdxdy q N N N N N N Y X Y X Y X 0000
000 取局部坐标s ,在i 点s=0,在j 点s=l ,L 为ij 边的长度。在ij 边上,以局部坐标表示
的插值函数为, L s N i -
=1,L
s
N j =,0=m N
载荷为L s
q q x =
qtL L s L s qt tds L s q L s X L L
i 6
1
)
32()1(0
23
20=
-=-=? qtL L s qt tds L s q L s X L L
j 3
130
2
3
===?
3.4单元刚度矩阵
根据单元的位移函数,
????
?
???????????????????????=??????m m j j i i m j
i
m j i v u v u v u N N N N N N v u 0
000 由几何方程可以得到单元的应变表达式,
{}?????
??????
??????????????????
?
??=??
??
??
?
??
???????????+??????=m m j j i i m m
j
j
i
i m j i m j i v u v u v u b c b c b c
c c c b b b A x v y u y v x u 00000021ε
(3-19)
记为
e B }]{[}{δε=,[B]矩阵称为几何矩阵。
[B]矩阵可以表示为分块矩阵的形式[][
]
m j i
B B B B =
[]????
?
?????=i i i i
i b c c b A B 0
021
(3-20)
由物理方程,可以得到单元的应力表达式,
{}[]{}[][]{}e B D D δεσ==
(3-21)
[D]称为弹性矩阵,对于平面应力问题,
[]?????
??????
?--=210
00101
)1(2μμ
μ
μE D 定义[][][]B D S =为应力矩阵。
将应力矩阵分块表示为,
[][]
m j i
S S S S =
[][][]??
?
??
?
?
???????---==i i i i i i
i i b c c b c b A E B D S 212
1)1(22μμ
μμμ
(3-22)
应用虚功原理可以建立单元结点位移与结点力的关系矩阵,单元刚度矩阵。
虚功原理:在外力作用下处于平衡状态的弹性体,如果发生了虚位移,则所有外力在虚位移上做的虚功等于内应力在虚应变上做的虚功。
单元的结点力记为{}[
]
T
m m j j i i
e
V U V U V U F =
单元的虚应变为{}[]{}e
B **
δε=
单元的外力虚功为,
{}(){}e
T e
F *δ
单元的内力虚功为,
{
}{}??tdxdy T
σε* 由虚功原理可得,
{}(){}{}{}??=tdxdy F T
e
T e
σεδ** (3-23)
{}[]{}{}()T
T e
T
e T
B B ][)(***δδε==
{}{}{}e e B D S δδσ]][[][==
{}(){}{}(){}??=e
T
T e
e
T e
tdxdy B D B F δδδ]][[][**
{}{}??=e T e tdxdy B D B F δ]][[][
(3-24)
定义[]??=tdxdy B D B K T e
]][[][为单元刚度矩阵。
在3结点等厚三角形单元中[B]和[D]的分量均为常量,则单元刚度矩阵可以表示为,
[]tA B D B K T e ]][[][=
(3-25)
单元刚度矩阵表示为分块矩阵:
[]????
??????=][][][][][][][][][mm mj mi jm jj ji im ij ii e
K K K K K K K K K K
]][[][][s T r rs B D B K =
(3-26)
??
?
?
??=sy ry sx
ry sy rx sx
rx rs K K K K K ,,,,][
(3-27)
3.5单元刚度矩阵的性质与物理意义
(一)单元刚度矩阵的物理意义
假设单元的结点位移如下:{}T e
]00000
1[=δ
由e e e K F }{][}{δ=,得到结点力如下:
??
??
?
?
????????????????=??????????????????????ix my ix mx ix jy ix jx ix iy ix ix m m j j i i K K K K K K V U V U V U ,,,,,,
(3-28)
ix ix K ,表示i 结点在水平方向产生单位位移时,在结点i 的水平方向上需要施加的结点力。
ix iy K ,表示i 结点在水平方向产生单位位移时,在结点i 的垂直方向上需要施加的结点力。
选择不同的单元结点位移,可以得到单元刚度矩阵中每个元素的物理含义:
sx rx K ,表示s 结点在水平方向产生单位位移时,在结点r 的水平方向上需要施加的结点力。
sx ry K ,表示s 结点在水平方向产生单位位移时,在结点r 的垂直方向上需要施加的结点力。 sy rx K ,表示s 结点在垂直方向产生单位位移时,在结点r 的水平方向上需要施加的结点力。 sy ry K ,表示s 结点在垂直方向产生单位位移时,在结点r 的垂直方向上需要施加的结点力。
因此单元刚度矩阵中每个元素都可以理解为刚度系数,即在结点产生单位位移时需要施
加的力。
(二)单元刚度矩阵的性质 1)对称性
利用分块矩阵的性质证明如下:
]][[][][s T r rs B D B K = ]][[][][r T s sr B D B K =
][]][[][][][][])][[]([][rs s T r s T T r T r T s T sr K B D B B D B B D B K ====
即T
e e K K )]([][= 2)奇异性
即单元刚度矩阵的行列式为零,0=e
K
。
将定单元产生了x 方向的刚体移动,T e ]010101[}{=δ,此时对应的单元结点力为零。
?????
?????????????????=??????????????????????010101][000000e K 可以得到,在单元刚度矩阵中1,3,5列中对应行的系数相加为零,由行列式的性质可知,0=e
K
。
同样如果假定单元产生了y 方向上的刚体位移T e ]101010[}{=δ,可以得到,在单元刚度矩阵中2,4,6列中对应行的系数相加为零。
3.6整体分析
得到了单元刚度矩阵后,要将单元组成一个整体结构,根据结点载荷平衡的原则进行分析,即整体分析。
整体分析包括以下4个步骤:
1)建立整体刚度矩阵,
2)根据支承条件修改整体刚度矩阵,
3)解方程组,求出结点的位移,
4)根据结点位移,求出单元的应变和应力。
在这里把结点位移作为基本未知量求解。
如何得到整体刚度矩阵?基本方法是刚度集成法,即整体刚度矩阵是单元刚度矩阵的集成。
图3.9
如图3.9所示,一个划分为6个结点、4个单元的结构。得到了每个单元的单元刚度矩阵后,要集成为整体刚度矩阵。
3.6.1刚度集成法的物理意义
由单元刚度矩阵的物理意义可知,单元刚度矩阵的系数是由单元结点产生单位位移时引起的单元结点力。
在如图3.9所示的结构中,使结点3产生单位位移时,在单元(1)中的结点2上引起结点力。由于结点2、3同时属于单元(1)、(3),在单元(2)中的结点2上同样也引起结点力,因此,在整体结构中当结点3产生位移时,结点2上的结点力应该是单元(1)、(2)在结点2上的结点力的叠加。
刚体集成法即结构中的结点力是相关单元结点力的叠加,整体刚度矩阵的系数是相关单元的单元刚度矩阵系数的集成。结点3在整体刚度矩阵的对应系数,应该是单元(1)、(3)、(4)中对应系数的集成。
3.6.2刚度矩阵集成的规则
1)将单元刚度矩阵中的每个分块放到在整体刚度矩阵中的对应位置上,得到单元的扩大刚度矩阵。
第2章 弹性力学平面问题有限单元法 2.1 三角形单元(triangular Element) 三角形单元是有限元分析中的常见单元形式之一,它的优点是: ①对边界形状的适应性较好,②单刚形式及其推导比较简单,故首先介绍之。 一、结点位移和结点力列阵 设右图为从某一结构中取出的一典型三角形单元。 在平面应力问题中,单元的每个结点上有沿x 、y 两个方向的力和位移,单元的结点位移列阵规定为: 相应结点力列阵为: (式2-1-1) 二、单元位移函数和形状函数 前已述及,有限单元法是一种近似方法,在单元分析中,首先要求假定(构 造)一组在单元内有定义的位移函数作为近似计算的基础。即以结点位移为已知量,假定一个能表示单元内部(包括边界)任意点位移变化规律的函数。 构造位移函数的方法是:以结点(i,j,m)为定点。以位移(u i ,v i ,…u m v m )为定点上的函数值,利用普通的函数插值法构造出一个单元位移函数。 在平面应力问题中,有u,v 两个方向的位移,若假定单元位移函数是线性的,则可表示成: (,)123 u u x y x y ααα==++ 546(,)v v x y x y ααα==++ (2-1-2)a 式中的6个待定常数α1 ,…, α6 可由已知的6个结点位移分量(3个结点的坐标) {}??? ?? ?????=????? ???? ?????????????=m j i m e d d d d m j j i v u v u v u i {} i i j j m X Y X (2-1-1)Y X Y i e j m m F F F F ?? ?? ???? ???? ??==??????????????????
第三章 平面问题的有限元法作业 1. 图示一个等腰三角形单元及其节点编码情况,设μ=0,单元厚度为t 。求 1)形函数矩阵[]N ;2)应变矩阵[]B ;3)应力矩阵[]S 。 4 第1题图 第2题图 2. 如题图所示,结构为边长等于a 的正方形,已知其节点位移分别为:11(,)u v 、 22(,)u v 、33(,)u v 、44(,)u v 。试求A 、B 、C 三点的位移。其中A 为正方形形心,B 为三角形形心。 3.直角边边长为l 的三角形单元,如题图所示。试计算单元等效节点载荷列阵(单元厚度为t ,不计自重)。 第3题图 第4题图 4. 如题图所示,各单元均为直角边边长等于l 的直角三角形。试计算(1)单元等效节点载荷列阵;(2)整体等效节点载荷列阵。已知单元厚度为t ,不计自重。
5.下列3个有限元模型网格,哪种节点编号更合理?为什么? 9 34 6 7912 11 34 6 12142 (a) (b) (c) 第5题图 6.将图示结构画出有限元模型;标出单元号和节点号;给出位移边界条件;并计算半带宽(结构厚度为t )。 2a (a) (b) 无限长圆筒 (c) 第6题图 7. 结构如图所示,已知结构材料常数E 和 ,单元厚度为t 。利用结构的对称性,采用一个单元,分别计算节点位移和单元应力。 第7题图
答案: 1. 1)形函数 i x N a = , j y N a = , 1m x y N a a =-- 2)应变矩阵 []1000101 000101011011B a -????=-??--???? 3)应力矩阵 []100010100 01 0111 110022 2 2S a ? ???-? ?=-????- -? ?? ? 2. A 点的位移为 ()2312A u u u = + , ()231 2A v v v =+ B 点的位移为 ()24313B u u u u = ++ , ()2431 3B v v v v =++ C 点的位移为 ()1223C a u u u = + , ()C 1223 a v v v =+ 3. 单元等效节点载荷列阵为 {}11 11 00003 663 T e i j i j R q q q q ?? =++?? ?? 4. (2)整体等效节点载荷向量为 {}111100006 322T R qlt P qlt P P qlt qlt ?? =-???? 7. (1) 减缩后的整体刚度方程 22 12 2 1222 22221110222021102(1)2 2102x x b b ab R b ab b P v Et ab a b ab ab R v b a μμμ μμμμμμ---??- - ??????????--?????? -??? ?=????---+ +? ???? ?????????-????+?? ? ? 节点位移
3 弹性力学平面问题的有限元法 本章包括以下的内容: 3.1弹性力学平面问题的基本方程 3.2单元位移函数 3.3单元载荷移置 3.4单元刚度矩阵 3.5单元刚度矩阵的性质与物理意义 3.6整体分析 3.7约束条件的处理 3.8整体刚度矩阵的特点与存储方法 3.9方程组解法 3.1弹性力学平面问题的基本方程 弹性力学是研究弹性体在约束和外载荷作用下应力和变形分布规律的一门学科。在弹性力学中针对微小的单元体建立基本方程,把复杂形状弹性体的受力和变形分析问题归结为偏微分方程组的边值问题。弹性力学的基本方程包括平衡方程、几何方程、物理方程。 弹性力学的基本假定如下: 1)完全弹性,2)连续,3)均匀,4)各向同性,5)小变形。 3.1.1基本变量 弹性力学中的基本变量为体力、面力、应力、位移、应变,各自的定义如下。 体力 体力是分布在物体体积内的力,例如重力和惯性力。 面力 面力是分布在物体表面上的力,例如接触压力、流体压力。 应力 物体受到约束和外力作用,其内部将产生内力。物体内某一点的内力就是应力。 图3.1
如图3.1假想用通过物体内任意一点p 的一个截面mn 将物理分为Ⅰ、Ⅱ两部分。将部分Ⅱ撇开,根据力的平衡原则,部分Ⅱ将在截面mn 上作用一定的内力。在mn 截面上取包含p 点的微小面积A ?,作用于A ?面积上的内力为Q ?。 令A ?无限减小而趋于p 点时,Q ?的极限S 就是物体在p 点的应力。 S A Q A =??→?0lim 应力S 在其作用截面上的法向分量称为正应力,用σ表示;在作用截面上的切向分量称为剪应力,用τ表示。 显然,点p 在不同截面上的应力是不同的。为分析点p 的应力状态,即通过p 点的各个截面上的应力的大小和方向,在p 点取出的一个平行六面体,六面体的各楞边平行于坐标轴。 图3.2 将每个上的应力分解为一个正应力和两个剪应力,分别与三个坐标轴平行。用六面体表面的应力分量来表示p 点的应力状态。应力分量的下标约定如下: 第一个下标表示应力的作用面,第二个下标表示应力的作用方向。 xy τ,第一个下标x 表示剪应力作用在垂直于X 轴的面上,第二个下标y 表示剪应力指 向Y 轴方向。 正应力由于作用表面与作用方向垂直,用一个下标。x σ表示正应力作用于垂直于X 轴的面上,指向X 轴方向。 应力分量的方向定义如下: 如果某截面上的外法线是沿坐标轴的正方向,这个截面上的应力分量以沿坐标轴正方向为正; 如果某截面上的外法线是沿坐标轴的负方向,这个截面上的应力分量以沿坐标轴负方向为正。 剪应力互等:xz zx zy yz yx xy ττττττ===,, 物体内任意一点的应力状态可以用六个独立的应力分量x σ、y σ、z σ、xy τ、yz τ、zx τ
Mmm 3 弹性力学平面问题的有限元法 本章包括以下的内容: 3.1弹性力学平面问题的基本方程 3.2单元位移函数 3.3单元载荷移置 3.4单元刚度矩阵 3.5单元刚度矩阵的性质与物理意义 3.6整体分析 3.7约束条件的处理 3.8整体刚度矩阵的特点与存储方法 3.9方程组解法 3.1弹性力学平面问题的基本方程 弹性力学是研究弹性体在约束和外载荷作用下应力和变形分布规律的一门学科。在弹性力学中针对微小的单元体建立基本方程,把复杂形状弹性体的受力和变形分析问题归结为偏微分方程组的边值问题。弹性力学的基本方程包括平衡方程、几何方程、物理方程。 弹性力学的基本假定如下: 1)完全弹性,2)连续,3)均匀,4)各向同性,5)小变形。 3.1.1基本变量 弹性力学中的基本变量为体力、面力、应力、位移、应变,各自的定义如下。 体力 体力是分布在物体体积内的力,例如重力和惯性力。 面力 面力是分布在物体表面上的力,例如接触压力、流体压力。 应力 物体受到约束和外力作用,其内部将产生内力。物体内某一点的内力就是应力。
图3.1 如图3.1假想用通过物体内任意一点p 的一个截面mn 将物理分为Ⅰ、Ⅱ两部分。将部分Ⅱ撇开,根据力的平衡原则,部分Ⅱ将在截面mn 上作用一定的内力。在mn 截面上取包含p 点的微小面积A ?,作用于A ?面积上的内力为Q ?。 令A ?无限减小而趋于p 点时,Q ?的极限S 就是物体在p 点的应力。 S A Q A =??→?0lim 应力S 在其作用截面上的法向分量称为正应力,用σ表示;在作用截面上的切向分量称为剪应力,用τ表示。 显然,点p 在不同截面上的应力是不同的。为分析点p 的应力状态,即通过p 点的各个截面上的应力的大小和方向,在p 点取出的一个平行六面体,六面体的各楞边平行于坐标轴。 图3.2 将每个上的应力分解为一个正应力和两个剪应力,分别与三个坐标轴平行。用六面体表面的应力分量来表示p 点的应力状态。应力分量的下标约定如下: 第一个下标表示应力的作用面,第二个下标表示应力的作用方向。 xy τ,第一个下标x 表示剪应力作用在垂直于X 轴的面上,第二个下标y 表示剪应力指 向Y 轴方向。
[] ()1 ,2 1 1112 1 =∴---++== i i i j k i j k i i k k j j i k k j j i i y x N y x y x y x y x y x y x y x y x y x A 即:()()()1,,,===k k k j j j i i i y x N y x N y x N (由i,j,k 轮换性知) 同理可证:()()0,,==k k i j j i y x N y x N (作业:证明:()()k j i j i y x N j j i ,,0,=≠=) 因此 ()()()()()()()()()()?? ??? ??===???≠===?======1 ,,0,,0,,1,00,,1,,0,0 ,,0,,1,k k k j j k i i k j i i k k j j j j i i j k k i j j i i i i y x N y x N y x N j i j i j N y x N y x N y x N y x N y x N y x N δ (2-12) 即形函数在自己节点上为1,在其余节点上为0。 2. 在单元上任意一点,三个形函数之和为1,即 ()()()1,,,=++y x N y x N y x N k j i 。 证明: ()()()()()()()[] y c c c x b b b a a a A y c x b a y c x b a y c x b a A y x N y x N y x N k j i k j i k j i k k k j j j i i i k j i ++++++++= ++++++++=++21 21 ,,,
3弹性力学平面问题的有限元法 本章包括以下的内容: 3.1弹性力学平面问题的基本方程 3.2单元位移函数 3.3单元载荷移置 3.4单元刚度矩阵 3.5单元刚度矩阵的性质与物理意义 3.6整体分析 3.7约束条件的处理 3.8整体刚度矩阵的特点与存储方法 3.9方程组解法 3.1弹性力学平面问题的基本方程 弹性力学是研究弹性体在约束和外载荷作用下应力和变形分布规律的一门学科。在弹性力学中针对微小的单元体建立基本方程,把复杂形状弹性体的受力和变形分析问题归结为偏 微分方程组的边值问题。弹性力学的基本方程包括平衡方程、几何方程、物理方程。 弹性力学的基本假定如下: 1)完全弹性,2)连续,3)均匀,4)各向同性,5)小变形。 3.1.1基本变量 弹性力学中的基本变量为体力、面力、应力、位移、应变,各自的定义如下。 体力 体力是分布在物体体积内的力,例如重力和惯性力。 面力 面力是分布在物体表面上的力,例如接触压力、流体压力。 应力 物体受到约束和外力作用,其内部将产生内力。物体内某一点的内力就是应力。
图3.1
如图3.1假想用通过物体内任意一点 p 的一个截面 mn 将物理分为I 、n 两部分。将部 分n 撇开,根据力的平衡原则, 部分n 将在截面 mn 上作用一定的内力。 在mn 截面上取包含 p 点的微小面积 A ,作用于:A 面积上的内力为:Q 。 令.\A 无限减小而趋于p 点时, Q 的极限S 就是物体在p 点的应力。 应力S 在其作用截面上的法向分量称为正应力,用 b 表示;在作用截面上的切向分量 称为剪应力,用T 表示。 显然,点p 在不同截面上的应力是不同的。 为分析点p 的应力状态,即通过p 点的各个 截面上的应力的大小和方向,在p 点取出的一个平行六面体,六面体的各楞边平行于坐标轴。 将每个上的应力分解为一个正应力和两个剪应力, 分别与三个坐标轴平行。用六面体表 面的应力分量来表示 p 点的应力状态。应力分量的下标约定如下: 第一个下标表示应力的作用面,第二个下标表示应力的作用方向。 xy ,第一个下标 x 表示剪应力作用在垂直于 X 轴的面上,第二个下标 y 表示剪应力指 向Y 轴方向。 正应力由于作用表面与作用方向垂直,用一个下标。 二x 表示正应力作用于垂直于 X 轴 的面上,指向X 轴方向。 应力分量的方向定义如下: 如果某截面上的外法线是沿坐标轴的正方向, 这个截面上的应力分量以沿坐标轴正方向 为正; 如果某截面上的外法线是沿坐标轴的负方向, 这个截面上的应力分量以沿坐标轴负方向 为正。 图3.2
一个受到集中力P 作用的结构,泊松比ν=61,m N P y 16=,t=1cm ,试按平面应力问题计算,采用三角形单元,求出节点位移。 解: 如图所示,划分三角形单元为四部分,并进行单元坐标编号,编程进行求解
单元①的刚度矩阵为: ???? ??????=333231232221 131211 1K K K K K K K K K K ()3,2,1===m j i 其中子矩阵表达式为: ???? ??????-+-+-+-+?-=s r s r s r s r s r s r s r s r rs b b v c c c b b c b c c b c c v b b Et K 21212121)1(42ννννν()m j i s r ,,,= E E E E E Et 2,,22210944.110944.121)611(401 .0)1(4--==?≈???? ??-??=? -ν 调用 Triangle2D3Node_Stiffness 函数,求出单元刚度矩阵: )3,2,1(0.2143 0 0 0.2143- 0.2143- 0.2143 0 0.5143 0.0857- 0 0.0857 0.5143- 0 0.0857- 0.5143 0 0.5143- 0.0857 0.2143- 0 0 0.2143 0.2143 0.2143- 0.2143- 0.0857 0.5143- 0.2143 0.7286 0.3000- 0.2143 0.5143- 0.0857 0.2143- 0.3000- 0.7286 '1===????????? ???????????=m j i E K
计算力学(有限元方法部分) 程序设计 程序说明书 程序1:平面问题的有限元分析 文件名:h01.m 算例文本:h01.txt 输出文本:result1.txt 使用前请先将h01.txt放入默认的文本读取路径(我的要求与m文件在同一文件夹内)! 文本输入顺序: 材料信息(编号、弹性模量、泊松比) 注意:材料编号须按1、2、3、4……的顺序排列 节点信息(编号、x坐标、y坐标) 注意:节点编号须按1、2、3、4……的顺序排列 约束信息(约束节点号、x方向有无约束、y方向有无约束、x方向位移、y 方向位移)有约束处填一正数,无约束处填0。无约束处请勿输入位移。 单元信息(厚度、材料编号、节点编号,若为3节点单元,则第四个编号填0) 注意:单元编号须按1、2、3、4……的顺序排列 集中力(作用节点号、x方向力、y方向力) 线荷载(作用边上的两个节点号、x方向分布力、y方向分布力) 面荷载(作用单元号、x方向分布力、y方向分布力) 程序可调部分: 4-6行中可以指定输出哪些图像(按顺序依次为节点、单元图像,x、y方向位移图像,xx、yy、xy方向应力图像),第7行中可以指定输入的.txt文本名称。 文本输出内容: 结点位移信息(节点号、x方向位移、y方向位移) 单元形心处的应变信息(单元号、x方向正应变、y方向正应变、xy方向工程切应变)
单元形心处的应力信息(单元号、x方向正应力、y方向正应力、xy方向切应力) 本程序附有三角形单元自动加密前处理部分h01auto.m,其算例文本: h01coarse.txt,输出文本:h01new.txt。它可以适用于本题的要求,在已给定本题所需全部信息的情况下将已有的单元加密为三角形单元。其输出文本可直接作为上面程序的输入文本。 h01.m h01.txt h01auto.m h01coarse.txt 欢迎交流与提问!附上邮箱:x67891@https://www.wendangku.net/doc/5812014986.html,。祝力学学习顺利!
第三章作业 3-1:试证明平面三角形单元内任一点的形函数之和恒等于1。 证明1:设单元发生X 方向的刚体位移0u ,则单元内到处应有位移0u ,有 0u u u u m j i === ()00u u N N N u N u N u N u m j i m m j j i i =++=++= 1=++m j i N N N 若位移函数不满足此要求,则不能反映单元的刚体位移,不能得到正确的结果。# 证明2:设P 是三角形内任一点,可用面积坐标表示为() m j i L L L P 。由面积坐标的定义和性质知1=++m j i L L L ,且三节点三角形的一点的面积坐标即为其形函数,故平面三角形单元内任一点的形函数之和恒等于1。# 3-2:试证明三角形单元的任一边上的一点的三个形函数与第三个顶点的坐标无关。 证明1:设k 是三角形ij 边上的任一点,点k 面积坐标得 0==m m L N # 证明2:三角形单元是协调单元,必须在单元边界上保持连续性,所以在单元边界上的点的位移只能由边上两个节点的形函数来贡献,否则就会撕裂和重叠,即(如在ij 边上的点) j j i i j j i i v N v N v u N u N u +=+= 故三角形的三边上的点的形函数只与边上节点的坐标有关,而与第三点无关。# 3-3:证明三角形单元是常应变单元。 证明: y x u 321ααα++=,y x v 654ααα++= 2αε=??= x u x 6αε=??= y v y
53ααγ+=??+??= x v y u xy # 即三角形单元是常应变单元。 3-4:已知单元刚度矩阵[][][][]tdxdy B D B k T A e ??= ,试说明[][]D B ,分别是什么矩阵,与单 元的那些特性有关?若厚度为t 的平面三角形常应变单元ijm 的单元刚度矩阵记为: [][][] [][][] []??? ? ? ?? ???=mm jm jj im ij ii k k k k k k k 说明子块[] ij k 的物理意义,并证明[]k 为对称矩阵。 解:[]B 是应变矩阵又称几何矩阵,与单元节点坐标有关;[]D 为弹性矩阵,与材料的弹性常数E 、μ有关。 []ij k 表示当节点j 处产生单位位移,其余节点完全被约束时,在节点i 处引起的节点力。 利用矩阵的运算关系 [] [][][][] [][][] []tA B D B tA B D B k T T T T T T T == 由于[]D 是对称矩阵,[][]D D T = 所以[][][][][]k tA B D B k T T ==,即[]k 为对称矩阵。# 3-5:图示平面等腰三角形单元,若3.0=μ,弹性模量为E ,厚度为t ,求形函数矩阵[]N 、应变矩阵[]B 及单元刚度矩阵[]K 。(补充题意:平面应力情况) 解:对平面等腰直角三角形建立图示坐标系。 x y