结构分析的有限元法-第七章

1. syms x e xi ei
2. syms xj ej
3. syms b a
4. syms E mu h
5. x0=xi*x;
6. e0=ei*e;
7. Ni=(1+x0)*(1+e0)*(2+x0+e0-x^2-e^2)/8;
8. Nxi=-b*ei*(1+x0)*(1+e0)*(1-e^2)/8;
薄板理论的应力分量
xxyy 1E21 0

1 0
0 0
xy D xy
(1)2 xy xy

x y




zD

w,xx w, yy
x w y b w b 1 (3 5 2 6 8 2 2 9 3 1 0 2 1 1 3 3 1 22 )
y w x a w a 1 ( 2 2 4 5 3 7 2 2 8 9 2 2 1 1 2 1 2 3 )
ww(x, y) u zw,x v zw,y
z

w z

0
x z u ,z w ,x w ,x w ,x 0
y z v ,z w ,y w ,y w ,y 0
x

u x
zw,xx
y

v y
zw,yy
xyv,xu,y2zw ,xy
p3183ipp44)
a(15i p1 i p2 3 p3 5ii p4 )
(i 1,2,3,4)
(7.39) (7.40) (7.41) (7.42)
1. syms x e xi ei 2. syms b a 3. syms p1 p2 p3 p4 4. x0=xi*x; 5. e0=ei*e; 6. Ni=(1+x0)*(1+e0)*(2+x0+e0-x^2-e^2)/8; 7. Nxi=-b*ei*(1+x0)*(1+e0)*(1-e^2)/8; 8. Nyi=a*xi*(1+x0)*(1+e0)*(1-x^2)/8; 9. N=[Ni Nxi Nyi]; 10. N1=(1-x)*(1-e)/4 ; 11. N2=(1+x)*(1-e)/4 ; 12. N3=(1+x)*(1+e)/4 ; 13. N4=(1-x)*(1+e)/4 ; 14. p = N1*p1+N2*p2+N3*p3+N4*p4 ; 15. f=int( int( p*transpose(N), x, -1, 1 ), e, -1, 1 )*a*b; 16. f=expand(f); 17. f=subs(subs(f,xi^2,1),ei^2,1) ; 18. f=simple(f)
18. kij=int(int(f,x,-1,1),e,-1,1)*a*b*h^3/12;
19. kij=subs(kij,xi^2,1);
20. kij=subs(kij,xj^2,1);
21. kij=subs(kij,ei^2,1);
22. kij=subs(kij,ej^2,1);
等效结点力
Ni (1 i )(1i) 4 (i 1,2,3,4)
那么式(7.39)的被积函数是多项式，可以积出显式如下
f
e pi


f M
zi xi
M yi


ab 180
4b5(p115i7p1iipi2p2
18i p3 5ii
第七章 板壳单元
板和壳是指厚度比其他尺寸要小得多的平面或曲面构 件，在工程中应用广泛。由于它的这种几何特点，前面所 述的三维单元并不十分适合用来分析它们的力学特性。因 为三维单元在三个方向的尺寸应尽量接近，否则求解精度 由于系统矩阵病态而大大降低。所以必须采用很细密的网 格来适应板和壳的几何特征，但是这将导致有限元模型的 自由度疯狂地增长，花费大量的计算和前后处理时间。因 此开发适合于板壳结构的专用单元是十分必要的。事实上， 60~70年代大量的关于有限元的研究其中很大一部分是在 板壳方面的工作。
9. Nyi=a*xi*(1+x0)*(1+e0)*(1-x^2)/8;
10. N1=[Ni Nxi Nyi];
11. N2=subs(subs(N1,xi,xj),ei,ej);
12. Bi=-1/a/b*[diff(N1,x,2)*b/a; diff(N1,e,2)*a/b; diff(diff(N1,x,1),e,1)*2];
薄板理论(Kirchhoff板理论)
y
z
中面
x
基本假设
（1）原先垂直于中面的直线，变形后仍然垂直于弯曲 的中面（直法线假设） （2）在横向荷载的作用下，中面既不伸长也不缩短
薄板理论的位移假设
ww(x, y)
u z w x
v z w y
z, w u
w,x
z
w
x, u
(b)
薄板理论的应变分量
d
(i, j 1,2,3,4)
(7.27)
把应变矩阵式(7.23)和式(7.3)中的弹性矩阵 D 代入上式，并完成积分运算得到矩
形薄板单元显式的单元刚度子矩阵
a11 a12 a13
K ij


a21
a22
a23

a31 a32 a33
(7.28)
单元刚度矩阵K的推导程序
x w y b w b 1 (3 5 2 6 8 2 2 9 3 1 0 2 1 1 3 3 1 22 )
y w x a w a 1 ( 2 2 4 5 3 7 2 2 8 9 2 2 1 1 2 1 2 3 )

b a
N i ,

1 4

3b a

i
(1


i
)
0
bi (1 3i )(1 i)
(7.24)

b a
N i ,

1 4

3a b

i
(1


i
)
ai (1 i )(1 3i)
0
(7.25)

2Ni ,

1 4
13. Bj=-1/a/b*[diff(N2,x,2)*b/a; diff(N2,e,2)*a/b; diff(diff(N2,x,1),e,1)*2];
14. D=E/(1-mu^2)*[1 mu 0
15.
mu 1 0
16.
0 0 (1-mu)/2] ;
17. f=transpose(Bi)*D*Bj ;
位移插值
Nxi, bi2 (1i)(12) 8bi(1i)(1i) 4
Nyi, ai2(1i)(12) 8ai(1i)(1i) 4
4
w (Niwi Nxixi Nyiyi)
Nyi, aii (1i)(12) 8
1
2
w (Niwi Nyiyi) i1
y1 ai 21(Ni,wi Nyi,yi)
单元应变矩阵B
将式(7.17)代入几何方程(7.4)，可以将单元应变用结点位移列阵表示为
ε Bδe zB1 B2 B3 B4 δe
(7.22)
式中
Bi




i1
xw ,yb 1i 4 1(N i,w iN xi,xiN yi,yi)
4

3
y w ,x 1 ai 4 1(N i,w i N x i, x i N y i, y i)

形函数
1
2
Ni (1i)1(i)2 (ii22)8 Nxibi (1i)1(i)1(2)8 Nyiai (1i)1(i)1(2)8
y
x
x
δi

wxii



wi w , yi


y
i

w , xi
fi


M
f
zi xi


M
yi
矩形单元的位移模式
w 1 2 3 4 2 5 6 2 7 3 8 2 9 2 1 0 3 1 1 3 1 2 3
Ni ,xx Ni ,yy






Ni , Ni ,
2Ni
,
xy

2Ni ,
a2
b2

(ab)


1 ab
bN aN
i i
, ,

2Ni ,
a b

(i 1, 2,3, 4) (7.23)
记号 N i , xx ， N i , 等分别表示 2 N i x 2 、 2 N i 2 等。按照式(7.21)容易得到
将矩形单元的 4 个结点坐标 (i ,i ) 和结点位移 (wi ,xi ,yi ) 分别代入上式，就可以得到 关于这 12 个参数的联立方程组。从中解出1 至12 ，经整理后，我们得到
4
4
w (N iw iN xi xiN yi yi) N iδi
i 1
i 1
ii (3
2
3 2

4)
bi (3 2 2i 1)
ai (3 2 2i 1)
(7.26)
单元刚度矩阵K
单元刚度矩阵可以写成分块形式，其子矩阵的公式是
Kij z2BTi DB jdxdydz
h3 12
1 1
1 1
BTi
DB
j
abd
D

Eh3
12(1
2)
薄板理论的控制方程 4xw 4 2x24wy2 4yw 4 D q
D Eh3
12(1 2 )
矩形单元
y
xi (Mxi)
z
ij
yi (Myi)
w (fzi)
o
x
图 7-4 平板划分成矩形单元
w,x ( x ) y
z
w, y ( y )
如果平板单元受分布横向面载荷 p 的作用，那么等效结点力是
f
e pi


f zi M xi



M

yi

1 1
1 1
pN Ti
abd
d
如果 p 在单元表面上线性分布，则可以表示成
(i 1,2,3,4)
4
p Ni pi i 1
式中 pi (i 1, 2,3, 4) 是四个结点上的载荷值， Ni 是四结点矩形平面单元的形函数，即
刚体位移
常应变
矩形的板单元是完备的
矩形单元的协调性 Ni, i(1i)(2ii2 2) 8(1i)(1i)(i 2) 8 Ni, i (1i)(2ii2 2) 8(1i)(1i)(i 2) 8
Nxi, bii(1i)(12) 8


xy

பைடு நூலகம்

2
w
,
xy

薄板理论的内力矩分量
h2
h2
h2
M x x zdz ， M y y zdz ， M xy xy zdz
h 2
h 2
h 2
M M MxxyyD10

1 0
0 0
(1)
22w ww,,,xyxxyy
Ni (1i)1(i)2 (ii22)8 Nxibi (1i)1(i)1(2)8 Nyiai (1i)1(i)1(2)8
矩形单元的完备性
w 1 2 3 4 2 5 6 2 7 3 8 2 9 2 1 0 3 1 1 3 1 2 3
1
1
3


6
2

2
10
3
2
2
3
15 4
3
22
3
4
矩形单元的形函数
w 1 2 3 4 2 5 6 2 7 3 8 2 9 2 1 0 3 1 1 3 1 2 3