有限元_4-薄板弯曲问题

1、有限元讲义第4章弹性薄板弯曲问题的有限元法板和壳是指厚度比其他尺寸要小得多的平面或曲面构件。由于它的这种几何特点，三维单元并不适合用来分析它们的变形。因为三维单元在三个方向的尺寸应尽量接近，否则求解精度由于“剪切自锁”(shear locking)或系统矩阵病态而大大降低，甚至得到错误的结果。所以必须采用很细密的网格来适应板和壳的几何特征，但是这将导致有限元模型的自由度疯狂地增长。仿照根据梁理论建立梁单元的思路，自然想到根据板理论建立板单元。这里讨论两种板理论，一是薄板理论，也被称为Kirchhoff 板理论，它忽略了板的横向剪切变形；另一种是Mindlin板理论，它考虑了板的横向剪切变形的影

2、响，适合于板的厚跨比较大的情形。后者也常被称为Reissner 板理论8或中厚板理论。根据这两种理论可以建立不同的板单元。薄板弯曲问题在理论上和应用上都具有重要意义,并有专门著作加以论述（如杨耀乾平板理论）。象其它弹性力学问题一样，用微分方程、差分法等经典方法所能求解的薄板问题很有限，一般只能解决等厚、小孔口、支承情况较简单的单跨板。故工程设计中以往多采用简化、近似、图表等方法来解决板的设计问题。在板的分析中，常取板的中面为xoy平面(如图)。平板结构按其厚度t与短边a的比值大小而分为:厚板（Thick plate）和薄板（Thin plate)两种。当ta5为绝对柔性板。)4.1wt基本理论

3、一、基本假定、略去垂直于中面的法向应力。（s z = 0），即以中面上沿Z方向的挠度W代表板的挠度)、变形前垂直中面的任意直线，变形后仍保持为垂直中面的直线。（法向假定t zx= 0 ,t zy = 0）、板弯曲时，中面不产生应力。(中面中性层假定)上述假定常称为薄板小挠度问题假定（or 柯克霍夫假定）。符合上述假定的平板即为刚性板。二、基本方法1将 U . 代 x2 e x e e y = - z 2 e xy w x2y2xy有限元讲义以上述假定为基础，板分析中常用挠度w 作为基本未知量,下面介绍以w 为基本未知量所导出的有关方程。、几何方程（应变挠度关系）弹性曲面沿x, y 方向的倾角从

4、中面取出一微小矩形ABCD，如图所示，设其边长为dx,dy，变形后弯曲成曲面ABCD设A点挠度w ,则沿x方向倾角(绕y轴)wxq y =wx(B点绕度 w +dx )沿y方向倾角(绕x轴)wyqx =wy(D点绕度 w +dy ) 沿x, y 方向位移作平行于 xoz 平面,设中面上点A到A1的距离为Z，变形后,A点有挠度W, 同时发生弯曲，曲面沿x方向的倾角为wx,根据法线假定,则A1点沿x方向的位移:xu = -zw(负号为方向与x相反)同理取 yoz 平面得：yv = -zw(4-1-1) Z平面的应变分量和曲、扭率基本假定，由于s z=t zx =t zy = 0, 故板内任意点的应

5、变与平面问题相同:vx+uxvyuye x =e y =e xy =V 入 2w - z 2w y2-2z xy(4-1-2)2w 2w 2w, ,为曲面在X,Y方向的曲、此为Z平面的应变挠度度几何方程。上式中的扭率，记为：22wcx -c= c y = - wy 2 c xy -22w 1m1- m有限元讲义 2x2 xy(4-1-3)所以,e= zc、物理方程（应力挠度关系）由于忽略z 对变形的影响, 因此z平面的应力应变关系具有与平面问题相同的形式: EEE 2(1-m)将(4-1-2)代入得:2 2s x 2 2 xy 或简写为：s= zD0 x(4-1-4)式中弹性矩阵:0 0 2

6、m10E1- m 2 0D0 =、内力方程（内力挠度关系）从板内取微元体tdxdy , 由其上正应力s x,s y 和剪应力t xy ,可在截面上合成合力矩:M x ( y0z 面上由s x产生的绕Y轴弯矩)3F= t D0cs=12 3 zFw(x, y)= Amn sinq(x,y)=qmn sin 2 + 2 绕x轴转角 q x(上节为沿 y方向倾角有限元讲义M y (x0z 面上由s y产生的绕X轴弯矩)扭矩:假定M xy（由剪应力产生，如图）M x,M y,M xy 分别表示单位宽度上的内力矩。如是，内力矩阵：M x x2 t t 3212 y2 2 简写成312(4-1-5)比较(

7、4-1-4)和(4-1-5)可得用内力矩表示的平板应力:t由此可见，平板上、下表面处的应力最大： sz= t = 6F22以上是薄板弯曲问题中的基本公式,从中可见其挠度W是弯曲问题中的基本未知函数。且由于忽略了z方向的变化,因此它只是x，y 的函数: w=w(x, y)。若w已知，则位移，内力、应力均可按上述相应公式求出。在经典解析法中，W(x, y)常设为三角级数形式。例如，四边简支矩形板的W(x, y)设为: (纳维尔解)sin m=1 n=1npybmpxa式中Amn 为待定系数。假定荷载sin m=1 n=1npybmpxa则可得位移函数:qmnnpybmpxasinsin1Dpw(x

8、, y)=2m2 n2 a b 44.2有限元分析方法一、矩形单元的典型形式将图示矩形薄板沿x,y方向划分成若干小矩形(常取等分)从中取出一小矩形(单元)，共有四个结点，每个结点有三个位移分量：挠度w,绕x、y轴转角) 挠度 w绕y轴转角 q y(上节为沿 x方向倾角即结点i的位移4wi wdi= qxi = q yi -wde=d 1 L d4 =w qx1 qy1 L w qx4 qy4L F 4 =f1 Mx1 My1 L f4 Mx4 My4在上式中,前10项取到了三次项的全部,最后两项则是从五个四次项(x第二项 q y =qx =wq y = -a2表示刚体转动 有限元讲义wy xi

9、 (i =1,K4)同理,相应的结点力 fi 竖向力 yi xT符号重新定义是为了有限元表示的方便，由此得单元结位移向量T节点力T TFe =F 1二、 位移模式（函数）、位移模式的选取插值多项式取为：w(x,y)=a1+a2x+a3y+a4x2+a5xy+a6y2+a7x3+a8x2y+a9xy2+a10y3 +a11x3y +a12xy3(4-2-1)4x 3y x 2y 2 xy 3 y 4)中选用了两个。没2选x y2 是因为它没有多一项与其配对，没选x4,y4它们在边界上结出的挠度函数是四次的，比x3y 和xy3 要高一次,较之更难满足边界的协调和条件。、位移模式的检验（三个基本要求

10、： 刚体位移，常应变，尽可能的边界协调） 前三项含单元的刚体位移状态：第一项a1与坐标x, y无关, 表示z方向的挠度是常量, 刚体移动第三项qx =a3 - wxy 二次项代表均匀变形状态:曲率2wx2= -2a4 ，2wy2= -2a6 ，= -2a52wxy 能保证相邻单元在公共边界上挠度的连续性。5w=NdN= N1 /yN y4 /yNxi = - bhi(1+xix)(1+hih)1-h 2N yi = axi(1+xix)(1+hih)(1-x 2)6有限元讲义 不能保证相邻单元在公共边界上法线转角的连续性。y = -b =常量，代入位移模式4-2-1，可知边界上的挠度W是x的以

11、单元边界为例，在此边界上三次函数，合并整理后可得：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代入后，此时qx = d1 +d2x+d3x2 +d4x3但对x来讲，1, 2结点只能提供2个已知条件，不能完全确定上式，故边界的法线转角不能保证连续性。因此，这种单元是非协调元，但可以验证这种非协调

12、远是能通过分片试验的。（即当单元划分不断缩小时，计算结果仍能收敛于精确解。）三、形函数和形函数矩阵。wy分别将单元结点1, 2, 3, 4的坐标值代入(4-2-1)，并事先求出qx =wx，q y = -，便可得到各结点的位移值。一共可得12个关于ai 的方程组，联立求解可得：e(4-2-2)形函数矩阵:N y4 N y4 /x N1N1/xNx4LLNx1LLN y1 L N4L L LL L L式中形函数:18(1+ xix)(1+hih)(2 + xix +hih -x 2 -h 2)Ni =)1818(4-2-3)(i=1 2 3 4)在上面的推导中,我们仍然选用了局部坐标(无因次坐标

13、)。局部坐标与整体坐标的关系为:1a16(y - y0)h = 2Nx2B= - N1 y22 N12N y4 N y42N y4xxy7有限元讲义V = z四、单元的几何矩阵B和内力矩阵S1几何矩阵B由前可知e= zc, 将（423）代入(4-2-4)得到几何矩阵：xyx 2222 y2 2x212N y1LL2Nx1LLLLL 2（4-2-5）或以子块形式表示：B=B1 B2 B3 B4。式中：2.内力矩阵S由基本方程（4-2-5）可得到：F=DBde =Sde（4-2-6） tB D BdxdyT有限元讲义S称为内力矩阵，把单元的四个结点坐标分别代入4-2-4，求得B后，即可获得S，各节

14、点内力矩阵S的显式：五、单元刚度矩阵a b-a -b由一般公式得:K=。将几何矩阵B和弹性矩阵D的表达式代入，积分可得薄板弯曲问题矩形单元的单元刚度矩阵的显示：8 ( ) - - = a b T dxdy y x q y x N R , ) , ( a N b 3 = q R = qab dxdy 1 N3 b x N 3 N 3 - y N 4 a 3有限元讲义9六、荷载等效变换由荷载等效变换的一般公式可得a b1法向均布荷载q代入上述公式得：- -a -b 3- 3N 3 1 b N y3 a 4 1 Nx4 ba b N y2 x2 a 1 N2 N1 3 x1 - N y1 3有限元讲

15、义102单元中心点受法向集中力P代入上述公式可得：七、位移边界条件对称、固定边和简支边上支点的已知位移条件如下：对称轴: 法线转角=0固定边: 挠度=0 (或已知值)边线转角=0 (或已知值)法线转角=0 (或已知值)简支边: 挠度=0 (或已知值)边线转角=0 (或已知值)自由边上节点的挠度、边线和法线转角均为特定参数，同内部节点一样。与板铰接的固定立柱，其节点挠度 = 0，也可以是已知值。八、计算例题例题1:计算图示四边固定方板方板的边长为l，厚度为t，弹性模型量为E，波松比=0.3，全板承受均布法向荷载，求薄板中的挠度和内力。单元划分：为了说明解题方法，采用最简单的网络22，即把方板分成

16、四个矩形单元。由于对称性，只需计算一个单元，例如，计算图中有阴影的单元，单元的节点编号为，。此时，单元的a, b是l4a =b =计算节点荷载:由前面的均布荷载计算公式得：ql219212 -l -l T12 l -l 12 l l 12 -l lR=边界条件：边界23和34为固定边，因此节点2, 3, 4的挠度、边线和法线转角均为零。边界12和14为对称轴，因此x1 =0、y1 =0。于是，在4个节点和12个位移分量中,只有一个待求的未知量 w1 。结构的代数方程组：这是一个单元的计算题目，单元刚度矩阵在此处即为总刚度矩阵。引入支承条件后，在总刚度矩阵中只取第一行、列元素，在方程组右端项中只

17、保留第一个元素。于是结构的代数方程为：ql 2168D015l28D015l 2k1w1 =(81-6m)w1 =P =2 b- a 2 b a 2 - b a2 - b - a P 8 同此解出 w1 = 0.00148ql 。其中有限元讲义4D0= 0.09158Et 3D0 =Et 312(1-m 2)内力:利用式(4-2-6)可求得方板中点力矩为：由表看出，网格越密，计算结果越接近于精确答案。还可看出，位移的精度一般比内力的精度高，这是因为在位移法中，位移是由基本方程直接求出的，而内力则是根据位移间接求出的。11有限元讲义4.3薄板有限元程序设计一、总框图根据弯曲板有限元分析方法的解题

18、过程，可写出其总框图如下： 输入原始数据 or CAI算等效结点力 形成荷载列阵 形成单元 定位向量 形 成 总 刚 单刚 解方程输出位移 几何矩阵B 弹性矩阵D 计算单元内力等 结束 下面结合程序对框图中的内容加以说明。二、子框图、单元坐标结点编号及单刚形式。为了取挠度向下为正，又能与前述坐标系统统一，特将前述坐标前翻180(如图)为了能适用板的弹型性分析,程序采用了应力元和弯曲元的组合形式，即每个结点考虑5个位移分量: U,V, W, x , y , 前2个为平面应力问题的未知量,后3个为弯曲板的结点未知量。当只作弹性分析时，平面应力元和弯曲元是非藕连的，即单刚的两个副块垣为0, 单刚的形式为：12有限元讲义w1 x1 y4 w4 x4 y4u1 v1 u4 v4 平面应力元 0Ke =（88）弯曲元0（1212）程序中单刚数组为 DK(20, 20), 子程序:Subroutine DG(A, B, E, T, U)为其形成单刚的子程序。、自动形成单元编号信息（ 单元信息数组：IB）。、结点定位向量。、形成荷载列向量。（ a. 结点力；b. 非结点力(只考虑均布力)）、总刚，Subroutine ZG(M, N, LD, A, B, E,