Finite Element Method
3.面力的移置
设三角形单元某边界s 上受面力q 作用,分量为x q ,y q ,则
{}{}T y x q q q =
取ds 则
{}{}s R d q d +=
由一般公式:
{}[]{}s T
s
e d q N p +=? 积分在边界s 上
以上三种载荷的等效节点荷载由公式e 导出 通常我们称:
{}{}{}{}[]{}{}[]{}???
?
?????===???s s
T e e A y x T e e e
T m e e td q N p d td p N p R N p 为荷载移量的一般公式: 几点说明: 1. 虚功等效?静力等效。 []?N 唯一性 2. 一般{}[]{}[]{}[]{}s
T
s
A
y
x
T
e
T
m e
e td q N d td P N R N
p ???++=
3. 更多节点的单元公式形式不变,但
[]N 不同
4. 虽然公式e 导出但对于面力和体力的计算都是很麻烦和困难
的
N 为x,y 的函数,若p, q 再为 x, y 的函数则更难,且?单移分限不好定。
因此,我们将来还要进一步把这个问题解决好。 四. 三角形单元的面积坐标(自然坐标,局部坐标) 1. 面积坐标的定义:
图示三角形单元 I ,j ,k 中任意一点m ,其位置可由xoy 坐标系中两个坐标来确定,即m(x,y)
若我们连接m i ,m j ,m k ,则形成了3 个小三角形?ijm, ?ikm, ?jkm. 则有:若m(x,y)确定??ijm, ?ikm, ?jkm.面积确定。 反之,?ijm, ?ikm, ?jkm.面积确定?m(x,y)确定 (用同底等高的概念解释!!) 因此,三角形单元内任一点可以??
?用三个三角形面积描述
用直角坐标描述
我们如何用三角形面积来描述m 点的位置呢? 定义:节点I 对边为底的三角形面积为i
A ;
节点j 对边为底的三角形面积为j
A ;
节点k 对边为底的三角形面积为k
A ;
设三角形单元的面积为A 令
???
?
?
????
==
=A A L A A L A A L k
K j J i i (2-37)
则三个比值i
L ,j L ,k
L 称为三角形单元中m 点的面积坐标.
2.三角形面积坐标的性质:
1》 面积坐标为三角形单元的局部坐标,与三角形的形状及位置无
关。其定义域为10,,≤≤k j i L L L ; 2》
三个面积坐标之和:i
L +j L +k
L =1.即只有两个面积坐标是独立
的。(2-38)
证明:i
L +j L +k
L =A A i +A A j +A A k
=A
1(i A +j A +k
A )=1 (亦可几何解
释)。 3》
三角形单元内与jk 边平行的直线上各点i
L 相同(轮换)。(同
底等高三角形i
L =A
A i
)
4》 形心处的面积坐标为: i
L =j L =k
L =1/3 (2-39)
5》
三角形单元节点的面积坐标为:
??
?
??=========1
,0,0:0
,1,0:0
,0,1:K j i k j i k j i
L L L k L L L j L L L
i 节点节点节点 (2-40)
证:节点I: i
A =A. j A =k
A =0.
3.三角形面积坐标与直角坐标及形函数的关系 下面我们来推导面积坐标与直角坐标的关系: 设m 点的坐标为m(x,y),m 为任一点 则:i
A =2
1 k
k
j j y x y x y x 1
1
1=2
1(j k j k k k j j y x y x xy y x y x xy ---++)
=2
1[(j k k j y x y x -)+(k j y y -)x +(j k x x -)
y ]
显然: =
i
a
j k k j y x y x -,=i
b
k j y y -,=
i
c j k x x -
∴ i
A =2
1(y c x b a i i i ++)
∴
????
???
??++=++=++==)(21)(21)(21
y c x b a A
L
y c x b a A L y c x b a A A A L k k k K
j j j J i i i i i (2-41)
i
L 与i N 表达式比较可知:三节点三角形单元的面积坐标就
是其形函数。(对于一般的情况:面积坐标永远是线性坐标而形函数可以是非线性的,以后我们可以把形函数用面积坐标表示) 即i N =i
L ,j j
L N
=,k
k L N = (2-42)
∴i
L 具有i N 的全部性质
式(2-41)还可写成矩阵的形式:
????
?
???????????????=??????????y x c b a c b a c b a A L L L k k
k
j j
j
i i i k j i 121 直?面 (2-44) 这就是直角坐标与面积坐标的转换关系。 下面的结果留给大家自己证明:
??????
??????????????=??????????k j i k j
i
k j i L L L y y y x x x y x 111
1 面?直 (2-45) 4.
面积坐标函数的运算
我们可以不加证明得地给出面积坐标函数的微积运算结果。(证明复杂麻烦用Г函数等)
1.偏导 设z=f(i
L ,j L ,k L ) i
L =g(x,y) (I= I ,j ,k)
则:?????????=??+??+??=????=??+??+??=??∑∑i
i k k j j i
i k j i i
i k k j j i
i L z c A L z c L z c L z c A y z L z
b A L z b L z b L z b A x z
21
)(2121)(21,, (2-46)
2. 面积分
A d d L L L A
y x k j i 2)!
2(!
!!?+++=
??
γβαγβαγβα (2-47)
其中α,β,γ为正整数; 0!?1, A: 三角形面积 ex: 3
2)!2001(!0!0!1A
A d d L A y x i =+++=
?? (I= I ,j ,k)
6
2)!2002(!0!0!22
A
A d d L A
y x i =+++=
??
122)!2011(!0!1!1A
A d d L L A y x j i =+++=??
60
2)!2111(!1!1!1A
A d d L L L A
y x k j i =+++=
??
3.线积分:
s d L L s
s j i )!
1(!
!++=
?βαβαβ
α (s 为直线长) (2-48)
以上公式要会用 注意s d 表示的边 五. 三角形单元的荷载移置
有了面积坐标与形函数的关系,我们即可对荷载移置进行计算了。 1.
集中力的移置
设m 点)(m m y x 作用有集中力{}{}T y x R R R = m 点的形函数为:
)(21
m i m i i m i y c x b a A
N ++=
(I= I ,j ,k) 等效节点荷载为:
?????
??????
???????????=?????????
??????
??
???????????=??????????????????????y m k x m k y m j
x m
j
y m i
x m i y x m k m k m j m j
m i m i yk xk yj xj yi xi R N R N R N R N R N R N R R N N N N N N p p p p p p 00000
0 这就是三角形单元内m 点作用有{}{}T y x R R R =的等效节点荷载。只要计算出m i N (I= I ,j ,k)即可。作为特例,考虑三角形单元形心处重力的移置。 形心坐标:
?????
++=++=)
(31
)
(31k j i c k j i c y y y y x x x x
)
(21
c i c i i c i y c x b a A N ++=)])((3
1)(31)()[(21k j i j k k j i k j j k k j y y y x x x x x y y y x y x A ++-+++-+-= )](31
)(3
1
)[(21k j j j i j k k j k i k k k j k k j j j k i j i j k k j y x y x y x y x y x y x y x y x y x y x y x y x y x y x A ---+++-+-+-+-=
)](31
)(32)[(21i j i k k i j i k j j k j k k j y x y x y x y x y x y x y x y x A -+-+-+-=
)](31[21i j i k k i j i j k k j y x y x y x y x y x y x A -+-+-= 3
1= {}{}T
R R -=0
∴
xi p =xj p =xk p =0 yi p =yj p =yk p =-3
1
R
故:重力作用于形心时各节点均担。 2. 体积力的移置
设单元作用有体力{}{}T
y x p p p =
则等效节点荷载为:
{}[]{}y x T
A d td P N p ??=={}T
A y k x
k y j x j y i x
i P N P N P N P N P N P N ??y x d td
i i L N =
∴
y x x A
i y x A
x i d td P L d td P N ????
= 若x p 为x, y 的函数,则把x p 用面积坐
标表示(转换) 在常体力的作用下有:
y A
x x i x i d td P L P N ??==A t
p x 2)!2001(!0!0!1?+++=A t p x 2!3!
1?=3At P x
{}{}T
y x
y x y x T
yk xk
yj xj yi xi
p p p p p p At
p p p p p p
3
=
即:常体力作用下,总体力均分三节点。 2.
面力的移置。
设三角形单元I ,j 边上作用有梯形分布的面力q
由面力移置公式得:(可分别由节点合力表示及用节点分力表示)
{}e p e
k j i p p p ??
?
??
?????=[][]s s T td q N ?=[
][]s
s T
k
j i
td
q N N N ? (q 为合力,非分力)
则 e i P =s s i td q N ?=s s i td q L ?
q 为x, y 的函数,把q 表示面积坐标的函数有q=j j i i q L q L + i L ,
j L 在门边上是线性坐标,可利用两点式方程写出。
则
:
e
i p =s s j j i i i td q L q L L ?+)(=s j j i s s i i td q L L td q L ??+2
=
t
q t q j i ?+++?++)!
111(!
1!1)!102(!0!2=t q t q j i ?+?!3!1!3!
2=)2(6
1j i q q st +
同理:
=e
j P s s j j i i j td q L q L L ?+)(=)2(6
1j i q q st + e
k P =s s j j i i k td q L q L L ?+)(=s s i i k td q L L ?+s s j j k td q L L ?
注意到在s 上k L =0
∴e
k
P =0
故:{}e p ={}T k j i P P P ={}T j i j i q q q q st 0226
1++ 或:{}e p =
{}T
yj
yi xj
xi yj
yi xj xi o o q q q q q q q q st
22226
++++
此法:1. 避免复杂的分离。 2. 便于编程计算。 特例:若分布荷载为三角形分布。令0=i q (或0=j q ) 则有:{}e p =
{}T
yj
xj
yj
xj o o q q q q st 226
(近端为2 ,远端为1)
说明:用以上积分的方法求等效节点荷载适用于任意节点的三角形单元,形函数也未必是线性的。 六. 三角形单元节点荷载的形成
经过荷载处理后,我们已把非节点的荷载转化为常点荷载。 实际计算的荷载为:
计算荷载=原节点荷载+等效节点荷载 即:{}=p {}+0p {}e p (2-49) 等效节点荷载要注意: 1. 同时贡献的问题 2.
用哪个单元计算的问题。
七.计算结果的整理:
有限元计算提供的结果一般为:1。节点位移 2。单元应力
1. 节点位移的处理:
一般把节点位移按比例标出,提供出结构常点位移分布规律 1. 连成折线(线性位移函数) 2. 连成光滑曲线(实际变形)
2.
单元应力的处理:
输出的单元应力一般为x σ,y σ,xy τ(形心处)(三节点Ⅱ单元为常应力元,无所谓) 1. 变换为单元的主应力。x σy σ,xy τ?1σ,2σ (材力) 2.
变换为节点应力的主应力。 (max σ) (平均法)
Ex. 节点5的应力为:
???
?
?????+++=+++=+++=)(41)
(41)(41432154321543215xy xy xy xy xy y y y y y x x x x x τττττσσσσσσσσσσ
即:σ=
n
n
n ∑σ1 (x, y, xy ) (x σ,y σ,xy τ)?(1σ,2σ)?max σ 然后标出应力变化曲线。
计算结果的工作量随结构的单元,节点划分增加面增大。 要关注的是:1)位移的变化规律 2)应力的最大值及发生地点。
小概念:位移最大的地方,应力未必最大。
八:有限元计算小结: 1. 基本原理:
连续法??→?离散有限个节点连接,有限大小的单元的组合法。 建立的节点位移{}?为未知数,总刚[]k 为系的n阶线性代数方程组。 2. 研究方法
(确定)节点位移???→?位移函数单元位移???→?几何方程单元应变???→?物理方程单元应力
{}? i i U N u ∑= {}[]{}?=B ε
{}[]{}?=S σ
???→?虚功原理单刚???→?节点平衡总刚???→?虚功等效计算荷载(等效节点荷载)
{}[]{}e e e k p ?= {}[]{}?=k p {}{}{}e p p p +=
????→?引入边界条件约束处理→求解方程→整理结果
3.解答特点:
1.假定单元内的位移分布规律,近似离散的数值解。
2.误差主要来源于:结构离散(连续?离散),假定位移分布。 3.收敛性:单元缩小?划分细密?收敛于精确解。
Chap3. 平面问题较精密单元的分析(矩形,高阶单元,等参单元) 3-1. 问题的提出
在三角形单元中,我们假定位移函数是线性的。即:单元内的位移按线性规律变化。这是最简单,最基本的一个有限单元。而实际结构中在外载荷的作用下位移分布常非按线性变化。
设单元位移曲线为图示的f(x).显然,用线性插值解的精度较差。提示解的精度的方法:
1.增加单元,节点数(工作量大,费用高)
2.提高插值阶数
因此,提出了用高阶插值,高阶单元的问题
我们大家都知道,位移是一个连续函数(连续体),而任意的连续函数都是可以展成幂级数,用幂级数来表示的。
因此,一个单元内的位移分布为f(x)时,我们就可以取级数的前几项来表示它。用二次三次函数来插值,以改善计算结果。
至于等参数单元(等参单元)是一种为清除曲边误差而出的一种单元。
如果实际结构为曲线边界,则无论怎样提高位移函数(插值)的阶数也不能使解得到多大的改善。有限个直边代替曲边,终究是代替,而不会是相等。为了处理曲边问题,人们提出了等参单元的概念。有平面等参单元,空间等参单元等。我们只向大家介绍平面等参单元,以供了解。
3-2四节点矩形单元的有限元分析。
矩形单元常用于规则边界的有限元分析,它也是常用的一种有限单元。
一.单元的位移函数。
设单元e为矩形单元,边长为2a,2b;其节点为I,j,k,m;为研究方便我们取局部坐标系x-y如图(原点在形心)。
1.单元的自由度及位移函数:
4个节点,每个节点2 个自由度(位移)u,v,则单元的总自
由度为8个。为保证单元的收敛性准则,位移函数必须保证有常数项,线性项。
设位移函数为:(对称性与坐标选择无关)
??
?+++=+++=xy
y x u xy y x u 87654321αααααααα (U 关于坐标对称) (3-1)
其中xy 项是根据pascal 三角形及xy 的对称性选取(选二次项尚有协调问题,选2x ,2y 项不行)
这样,已知i i v u (i=1,2,3,4)这八个节点位移??→?确定i α(i=1,2…8)
2.形函数的推导: i
i
i
i
i
y x y x u 43
21
αααα+++=
∴同理:??
???
?
?-+-=+++=--+=+--=4
321432143214
321ααααααααααααααααab b a u ab b a u ab b a u ab b a u m k j i ?????
??-+-=+++=--+=+--=8765876587658765ααααααααααααααααab b a v ab b a v ab b a v ab b a v m k j i 我们不难从前4个方程中解出1α,2α,3α,4α,具体做法:
①+②:??
?+=+-=+3
1312222ααααb u u b u u m k j i ? ①+②:)(4
11m k j i u u u u +++=α?v ③+④: ②-①:
)(41
3m k j i u u u u b
++--=
α?v
①-
②: ??
?+=-+-=-4
2422222ααααab a u u ab a u u m k j i ? )(412m k j i u u u u a
-++-=α?v ③-④ )(414m k j i u u u u ab
-+-=α?v
∴ xy y x u 4321αααα+++=
=)]()(1)(1)[(4
1m k j i m k j i m k j i m k j i u u u u ab
xy y u u u u b
x u u u u a
u u u u -+-+++--+-++-++++
=])1()1()1()1[(4
1m k j i u ab
xy b y a x u ab
xy b y a x u ab
xy b y a x u ab
xy b y a x -+-+++++--+++--
?????
??????+-=-+-++=+++-+=--+--=---=+
--)1)(1(1)1)(1(1)
1)(1(1)1)(1()1()1()1(b
y a x ab xy b y a x b y a x ab xy b y a x b y
a x a
b xy b y a x b y
a x a x
b y a x ab xy b y a x ∴u=])1)(1()1)(1()1)(1()1)(1[(4
1m k j i u b
y a
x u b
y a
x u b
y a
x u b
y a
x +-++++-++--
令:???
??
??????
+-=++=-+=--=)1)(1(41)
1)(1(41)1)(1(41)1)(1(41b y a x N b y a x N b
y
a x N
b y a x N m k j i (3-2) 则:
??
??
?=+++==++=∑∑m k j i i
i m m k k j j i i i
m k j i i k k j j i i v N v N v N v N v N v u N u N u N u N u ,,,,,, (3-3) 式(3-2)称为节点矩形单元的形函数;式(3-2)尚可写为:
)1)(1(41i
i i y y
x x N ++=
(I =I, j, k, m) (3-4) 式(3-3)为节点位移表示的单元位移函数。 式(3-3)还可以写成矩阵的形式
[]{}[]
{}
e
m k
j i
e
e
N I N I N I N
I N v u ?=?=?
?????][][][][
=
??
???
?m k
j
i
m k j i
N N N N N N N N 0
00000{}
e
? ??
?
???=1001][I (3-5)
式中:{}e ?={}e m m k k j j i i v u v u v u v u (3-6) 二. 形函数的性质:
1.形函数i
N (I =I, j, k, m )在节点I 上i
N =1,在其余节点上i
N =0
(轮换)
即在 节点I :i N =1 0=j
N
0=k N 0=m N
节点j: i
N =0 1=j
N
0=k N 0=m N ?i N (j)=
j i δ=??
?≠=i
j i j 0
1
节点k: i
N =0 0=j
N 1=k N 0=m N 节点m: i
N =0 0=j
N
0=k N 0=m N
证明:)1)(1(41i
i i
y y
x x N
++=
(i=I, j, k, m ) 在节点i
x=i x ,y=i y 时:
)1)(1(41i
i i y y
x x N ++=
=1 j, k, m 各节点至少有一个节点坐标x=-i x 或y=-i y 故i N =0 (i=j, k, m )
同理可得到全部结果。 2.4个形函数之和:i
N +j N +k N +1=m
N
证明:
写出形函数(a
x ,b
y 的符号与该点的i
x ,i
y 相同)
)1)(1(41b y a x N i --= )1)(1(41b y
a x N j -+= )1)(1(41
b y a x N k ++=
)1)(1(41b
y a x N m +-= ∴i
N +j
N +k
N +m
N =)]1)(1()1)(1()1)(1()1)(1[(41b
y a x b y a x b y a x b y a
x +-++++-++--
=)]1()1()1()1[(41ab
xy b y a x ab xy b y a x ab xy b y a x ab xy b y a x -+-+++++--+++--=1
因此:4个形函数只有3 个是独立的。
在I j 边上,i
N ,j N ≠0(节点除外),k N =m
N =0 (轮换)(一条边上
的四个i
N )
证明: )1)(1(41i
i i
y y x x N
++=
在I j 边上,y=i y , ∴)1(21)11)(1(4
1i
i i x x
x x N +=++= (x ≠-i x ) )1)(1(41j
j j
y y
x x N ++=
在I, j 边上,y=j y
∴)1(21)11)(1(41j
j j
x x
x x N
+=++=
(x ≠-j x ) )1)(1(41k
k k y y
x x N ++=
∴I, j 边上 y=i
y =j
y =-k
y
∴)1)(1(4
1k
k k k y y
x x N -++
==0 同理:I, j 边上:m N =0。 证毕 4.i N 在4 条边界上的性质(节点除外)
在包含节点I 的边界上,i N ≠0,否则i N =0。(轮换)(四条边上的一个i N )
证明:性质3 的另一种表述
)1)(1(4
1i
i i y y x x N ++
= 在包含节点 I 的边界上: X=i x , 或 y=i y 显然有:)1(21i
i
y y N
+=
1?i y y
(y ≠
i y )或)
1(2
1i
i x x N +
= (节点j, m 除外)
在不包含节点I 的边界上: X=- i x 或y=-i y 显然:i N =0 得证
由以上i N 的性质,我们可以描述i N 的几何图形。
三. 位移函数的性质
1.
位移函数是双线性的
位移函数:xy y x u 4321αααα+++= xy y x v 8765αααα+++=
显然, u, v 包含坐标x y 的二次项x y ,但当x=const 或y=const 时
U,V 都是一个线性函数。 即:
在单元内任一点,无论沿x 方向变化(此时y=const )或沿y 方向变化(x=const ) u, v 都是线性的。 2.
位移函数解满足收敛准则:
① 解反映单元的刚体位移。(位移=刚+弹) 位移函数 xy y x u 4321αααα+++= xy y x v 8765αααα+++= 写成如下的形式: xy y y x u 43
663212
2ααααααα+--++
+= xy x x v 87366352
2ααααααα++-+++= 显然:第一项1α,5α与x, y 无关。反映了单元内各点沿x, y 方向的刚体位移。
若令: w=2
3
6αα-,则w 反映单元内各点绕z 轴的刚体转动。
② 解反映单元的常应变 由几何方程:
y x
u
x 42ααξ+=??=
x y
v
y 87ααξ+=??=
)()(8463y x x
v
y u r xy αααα+++=??+??= 显然,2α,7α分别反映了沿x, y 方向的常应变;63αα+反映了(2121,,γγγγγ+=xy )常量的剪切应变。
面4α和8α分别反映了线性变化的x ξ,y ξ,xy γ。(4α,8α≠0)
y x 42ααξ+=
x y 87ααξ+=
)()(8463y x r xy αααα+++=
1α,5α反映单元的刚体移动。 6α-3α——反映单元的刚体转动。
一般:位移函数只要包含1α+x 2α+y 3α 选择变量的一次式,则必须保证收敛性。
单元内各应变都不是常量,你如何能解释其收敛? 解释:当单元尺寸逐步缩小时,单元内各点x, y 的变化必然很小。 以x ξ为例:
单元尺寸逐步缩小?单元内min max y y -逐步缩小
可以保证单元内x ξ以2α为基准,收敛于2α附近。 (2α≠0)
否则:
若2α=0,则x ξ=y 4α?在y=0处x ξ=0?y=0处,永远仅发生刚体位移。
同理可知:7α,63αα+的意义。
由于有4α,8α≠0的存在,四边形单元称为非常应变单元。 ③ 位移函数在单元内连续,在边界上与相邻单元协调。 证明:u, v 是连续的——明显
由于u, v 是双线性函数,而单元的每条边界都满足x=const 或y=const.
因此:在每条边界上: u, v 都是一个线性函数。
设1e ,2e 为相邻单元; 1e 的节点为:I, j, k, m ,2e 的节点为:I, p, q, j. 则I, j 为公共边界。
第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/4b11840245.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