有限元_4-薄板弯曲问题

《有限元》讲义
第4章
弹性薄板弯曲问题的有限元法
板和壳是指厚度比其他尺寸要小得多的平面或曲面构件。由于它的这种几何特点，三维单 元并不适合用来分析它们的变形。因为三维单元在三个方向的尺寸应尽量接近，否则求解精度 由于“剪切自锁”(shear locking)或系统矩阵病态而大大降低，甚至得到错误的结果。所以 必须采用很细密的网格来适应板和壳的几何特征，但是这将导致有限元模型的自由度疯狂地增 长。 仿照根据梁理论建立梁单元的思路，自然想到根据板理论建立板单元。这里讨论两种板理 论，一是薄板理论，也被称为Kirchhoff 板理论，它忽略了板的横向剪切变形；另一种是Mindlin 板理论，它考虑了板的横向剪切变形的影响，适合于板的厚跨比较大的情形。后者也常被称为 Reissner 板理论[8]或中厚板理论。根据这两种理论可以建立不同的板单元。 薄板弯曲问题在理论上和应用上都具有重要意义,并有专门著作加以论述（如杨耀乾《平 板理论》）。象其它弹性力学问题一样，用微分方程、差分法等经典方法所能求解的薄板问题 很有限，一般只能解决等厚、小孔口、支承情况较简单的单跨板。故工程设计中以往多采用简 化、近似、图表等方法来解决板的设计问题。 在板的分析中，常取板的中面为xoy平面(如图)。平板结构按其厚度t与短边a的比值大小 而分为: 厚板（Thick plate）和 薄板（Thin plate)两种。
阵
[S]的显式：
五、单元刚度矩阵
由一般公式得:[K]=
t[B] [D] [B]dxdy。将几何矩阵[B]和弹性矩阵[D]的表达式代入，积分可
a b T -a -b
得薄板弯曲问题矩形单元的单元刚度矩阵的显示：
8
《有限元》讲义
六、荷载等效变换
由荷载等效变换的一般公式可得
{}
1．法向均布荷载q
代入上述公式得：
当
t << 1 时称为薄板 a
平板上所承受的荷载通常有两种: 1. 面内拉压荷载。 由面内拉压刚度承担, 属平面应力问题。 2. 垂直于板的法向荷载, 弯扭变形为主,具有梁的受力特征, 即常说的弯曲问题。平板在 垂直于板面的荷载作用下产生挠度W。 当最大挠度w远小于t时, 称为小挠度问题(or刚性板)(stiffness plate) 当最大挠度w与t相差不大时,称为大挠度问题(or柔性板)(flexure plate) (工程定义:
② 二次项代表均匀变形状态: 曲率
2w 2w 2w = -2a4 ， = -2a6 ， = -2a5 x2 y2 xy
③ 能保证相邻单元在公共边界上挠度的连续性。
5
《有限元》讲义
④ 不能保证相邻单元在公共边界上法线转角的连续性。 以单元１～２边界为例，在此边界上 三次函数，合并整理后可得：
比较(4-1-4)和(4-1-5)可得用内力矩表示的平板应力:
[s]=12 3 z{F}
t
由此可见，平板上、下表面处的应力最大： {s}z=± t = ± 6{ F} 2
2
以上是薄板弯曲问题中的基本公式,从中可见其挠度W是弯曲问题中的基本未知函数。且由于忽略了 z方向的变化,因此它只是x，y 的函数: w=w(x, y)。若w已知，则位移，内力、应力均可按上述相应公 式求出。在经典解析法中，W(x, y)常设为三角级数形式。例如，四边简支矩形板的W(x, y)设为: (纳 维尔解)
３、板弯曲时，中面不产生应力。(─中面中性层假定) 上述假定常称为薄板小挠度问题假定（or 柯克霍夫假定）。符合上述假定的平板即为刚性板。
二、基本方法
1
《有限元》讲义
以上述假定为基础，板分析中常用挠度w 作为基本未知量,下面介绍以w 为基本未知量所导出的 有关方程。
１、几何方程（应变─挠度关系）
①弹性曲面沿x, y 方向的倾角 从中面取出一微小矩形ABCD，如图所示，设其边长为dx, dy，变形后弯曲成曲面A'B'C'D' 设A点挠度
e x = x
u
此为Z平面的应变─挠度度几何方程。上式中的 扭率，记为：
2w 2w 2w 为曲面在X,Y方向的曲、 , , x2 y2 xy
2
《有限元》讲义
c x
2w x 2
2
{c}= c y =
- wy 2 c xy -22w x y
w w ，便可得到 ，q y = y x
{w}=[N]{d}
e
(4-2-2)
形函数矩阵:
N1 [ N]= N1 /y N1/x
Nx1 L L
N y1 L L
L N4 N x4 L L L L L L
N y4 /y N y4 /x N y4
式中形函数:
2
[B]= -
（4-2-5）
L 2
或以子块形式表示： 式中：
[B]=[B1 B2 B3 B4]。
2.内力矩阵[S]
由基本方程（4-2-5）可得到：
{F}=[D][B]{d}e =[S]{d}e
7
（4-2-6）
《有限元》讲义
[S]称为内力矩阵，把单元的四个结点坐标分别代入4-2-4，求得[B]后，即可获得[S]，各节点内力矩
a10y3 +a11x3y +a12xy3
(4-2-1)
4
在上式中 x 3y x,前 2y10 2 项取到了三次项的全部 xy 3 y 4 中选用了两个。没 ,最后两项则 选x
2
)
y2 是因为它没有多一项与其配对，没选x4,y4它们在边界上结出的挠度函数是四次的，比x3y 和
xy3 要高一次,较之更难满足边界的协调和条件。
(4-1-3)
所以,
{e}= z{c}
由于忽略σz 对变形的影响, 因此z平面的应力─应变关系具有与平面问题相同的形
２、物理方程（应力─挠度关系）
式:
E E 2(1-m) E
将(4-1-2)代入得:
2
s x
xy


2
2 2


或简写为： 式中弹性矩阵:
{s}= z[D0]{x}
{d}e=[d 1
T L d4] = w q x1 qy1 L w qx4 qy4
[ ]
节点力
{F}e =[FL 1
二、 位移模式（函数） １、位移模式的选取
插值多项式取为：
F 4] =[f1 Mx1 My1 L f4 Mx4 My4]
T
T
w(x,y)=a1+a2x+a3y+a4x2+a5xy+a6y2+a7x3+a8x2y+a9xy2+
w 1 1 w 5为柔性板； 为刚性板； t 5 5 t w >5 t
为绝对柔性板。)
4.1
基本理论
一、基本假定
１、略去垂直于中面的法向应力。（s z = 0），即以中面上沿Z方向的挠度W代表板的挠度) ２、变形前垂直中面的任意直线，变形后仍保持为垂直中面的直线。（─法向假定t zx
= 0 ,t zy = 0）
1 (y - y0) 6
6
《有限元》讲义
V=z
四、单元的几何矩阵[B]和内力矩阵[S] 1．几何矩阵[B]
由前可知
{e}= z{c}, 将（4－2－3）代入(4-2-4)得到几何矩阵：
2N 1 2 x2 2 2 N1 xy N1 y2 2Nx1 L L x2 2N y1 L L x 2 L L 2N y4 2 N y4 2N y4 x y2 xy
w(x, y)= Amn sin
m=1 n=1


mpx npy sin a b
式中
Amn 为待定系数。
假定荷载
q(x,y)=qmn sin
m=1 n=1


mpx npy sin a b
则可得位移函数:
w(x, y)=
mpx npy 1 qmn sin sin 4 2 a b Dp m2 n2 2 + 2÷÷ a b
) )
4
《有限元》讲义
wi w {di}=qxi = w y (i =1,K4) - qyi w x i
同理,相应的结点力
fi
yi
竖向力

T
x
符号重新定义是为了有限元表示的方便，由此得单元结位移向量
２、位移模式的检验
（三个基本要求： 刚体位移，常应变，尽可能的边界协调） ① 前三项含单元的刚体位移状态： 第一项a1与坐标x, y无关, 表示z方向的挠度是─常量, 刚体移动
第二项 q y = 第三项
- w x q y = -a2 表示刚体转动 w qx =a3 qx = y
w,
则沿x方向倾角(绕y轴)
q y = x (B’点绕度 w + dx ) x
沿y方向倾角(绕x轴)
w
w
qx =
w w (D’点绕度 w +于 xoz 平面,设中面上点A到A1的距离为Z，变形 后,A点有挠度W, 同时发生弯曲， 曲面沿x方向的倾角为 的位移:
y = -b =常量，代入位移模式 4-2-1，可知边界上的挠度W是x的
w1-2 = c1 + c2x + c3x2 + c4x3
两个端点共有4个边界条件,(结点1,2的挠度W1 , W2 ,和转角
qy1, qy2 。利用他们可唯一确定四个常
数C1 ～C4。因为相邻单元在结点1, 2的W, θy对应相同，则两个单元依据四个条件得到的C1 ～C4 亦 相同，即两单元在边界具有同一挠度函数W。 ⑤ 法线转角 仍以1-2边界为例，将y=-b代入后，此时
(4-1-4)
m 0 E 1 m 1 0 D0 = 1- m 2 1 m 0 0 2
３、内力方程（内力─挠度关系）
从板内取微元体tdxdy , 由其上正应力s x,s y 和剪应力t xy ,
可在截面上合成合力矩:
M x ( y0z 面上由s x产生的绕Y轴弯矩)
Ni = 1 (1+ xix)(1+hih)(2 + xix +hih -x 2 -h 2) 8
1 8
Nxi = - bhi(1+xix)(1+hih)1-h 2
(4-2-3)(i=1 2 3 4)
)
N yi =
1 axi(1+xix)(1+hih)(1-x 2) 8
在上面的推导中,我们仍然选用了局部坐标(无 因次坐标)。局部坐标与整体坐标的关系为: 1 a h=
qx = d1 +d2x+d3x2 +d4x3
但对θx来讲，1, 2结点只能提供2个已知条件，不能完全确定上式，故边界的法线转角不能保证连 续性。 因此，这种单元是非协调元，但可以验证这种非协调远是能通过分片试验的。（即当单元划分不断 缩小时，计算结果仍能收敛于精确解。）
三、形函数和形函数矩阵。
分别将单元结点1, 2, 3, 4的坐标值代入(4-2-1)，并事先求出qx = 各结点的位移值。一共可得12个关于ai 的方程组，联立求解可得：
w , x
根据法线假定,则A1点沿x方向
u = -zw
同理取 yoz 平面得：
x
(负号为方向与x相反) (4-1-1)
v = -zw
y
③ Z平面的应变分量和曲、扭率 基本假定，由于s z
=t zx =t zy = 0, 故板内任意点的应变与平面问题相同:
2w -z x2 e x v 2w V 入 ey = { e } ＝ e y = - z 2 将 U . 代 y y (4-1-2) e 2 u v xy w e xy = + -2z xy y x
a
b
(
) 1 b 3 - 3 1 b 3 a 3 qab dxdy b - 3 a 1 3 b a3 - 3
4.2
有限元分析方法
一、矩形单元的典型形式
将图示矩形薄板沿x,y方向划分成若干小矩形(常取等分) 从中取出一小矩形(单元)，共有四个结点，每个结点有三个位移 分量： 挠度w,绕x、y轴转角
挠度 w 绕x轴转角 q x(上节为沿 y方向倾角 绕y轴转角 q y(上节为沿 x方向倾角
即结点i的位移
3
《有限元》讲义
M y (x0z 面上由s y
扭矩: 假定
产生的绕X轴弯矩)
M xy（由剪应力产生，如图） M x,M y,M xy 分别表示单位宽度上的内力矩。如是，内力矩阵：
M x
t 2
t 2
3
2

y
x2
12
简写成
{F}= t [D0]{c}
12
3
(4-1-5)