有限元薄板弯曲问题分析

《有限元》讲义第4章弹性薄板曲折问题的有限元法薄板曲折问题在理论上和应用上都拥有重要意义,并有特地著作加以论述（如杨耀乾《平板理论》）。象其余弹性力学问题相同，用微分方程、差分法等经典方法所能求解的薄板问题很有限，一般只好解决等厚、小孔口、支承状况较简单的单跨板。故工程设计中过去多采纳简化、近似、图表等方法来解决板的设计问题。在板的分析中，常取板的中面为xoy平面(如图)。平板结构按其厚度t与短边a的比值大小而分为:厚板（Thickplate）和薄板（Thinplate)两种。t当a1时称为薄板平板上所承受的荷载平常有两种:面内拉压荷载。由面内拉压刚度担当,属平面应力问题。2.垂直于板的法向荷载,弯扭变形为主,拥有梁的受力特色,即常说的曲折问题。平板在垂直于板面的荷载作用下产生挠度W。当最大挠度w远小于t时,称为小挠度问题(or刚性板)(stiffnessplate)当最大挠度w与t相差不大时,称为大挠度问题(or柔性板)(flexureplate)w11w(工程定义:为刚性板；55为柔性板；t5tw)5为绝对柔性板。t4.1基本理论一、基本假定１、略去垂直于中面的法向应力。（z0），即以中面上沿Z方向的挠度W代表板的挠度)２、变形前垂直中面的任意直线，变形后仍保持为垂直中面的直线。（─法向假定zx0,zy0）３、板曲折时，中面不产生应力。(─中面中性层假定)上述假定常称为薄板小挠度问题假定（or柯克霍夫假定）。符合上述假定的平板即为刚性板。二、基本方法1《有限元》讲义以上述假定为基础，板分析中常用挠度w作为基本未知量,下边介绍以w为基本未知量所导出的有关方程。１、几何方程（应变─挠度关系）①弹性曲面沿x,y方向的倾角从中面拿出一细小矩形ABCD，以以下图，设其边长为dx,dy，变形后曲折成曲面A'B'C'D'设A点挠度w,则沿x方向倾角(绕y轴)yw(B’点绕度wwdx)xx沿y方向倾角(绕x轴)xw(D’点绕度wwdy)yy②沿x,y方向位移作平行于xoz平面,设中面上点A到A1的距离为Z，变形后,A点有挠度W,同时发生曲折，曲面沿x方向的倾角为w,则A点沿x,依据法线假定1x方向的位移:uzw(负号为方向与x相反)x同理取yoz平面得：vzw(4-1-1)y③Z平面的应变重量和曲、扭率xyxy

基本假定，因为zzxzy0,故板内任意点的应变与平面问题相同:u2wxz2vxx将U.V代入2w＝yyz2(4-1-2)yuvxy2z2wyxxy2《有限元》讲义2w2w2w此为Z平面的应变─挠度度几何方程。上式中的x2,y2,xy为曲面在X,Y方向的曲、扭率，记为：2w2x2wx(4-1-3)yy22wxyxy所以,

z２、物理方程（应力─挠度关系）因为忽视σz对变形的影响,所以z平面的应力─应变关系拥有与平面问题相同的形式:xExy12Ey12yxExyxy21将(4-1-2)代入得:Ez2w2wx12x2y2Ez2w2wD0y12y2x22xyEzw1xy或简写为：zD0x(4-1-4)式中弹性矩阵:10D0E10120012３、内力方程（内力─挠度关系）从板内取微元体tdxdy,由其上正应力x,y和剪应力xy,3《有限元》讲义可在截面上合成合力矩:Mx(y0z面上由x产生的绕Y轴弯矩)My(x0z面上由y产生的绕X轴弯矩)扭矩:Mxy（由剪应力产生，如图）假定Mx,My,Mxy分别表示单位宽度上的内力矩。如是，内力矩阵：Mx2w2ttt3xFMy2tzdz2tz2D0dzD02w2Mxy22122y2wxy简写成t3(4-1-5)FD012比较(4-1-4)和(4-1-5)可得用内力矩表示的平板应力:12zFt3因而可知，平板上、下表面处的应力最大：z6ttF2以上是薄板曲折问题中的基本公式,从中可见其挠度W是曲折问题中的基本未知函数。且因为忽视了z方向的变化,所以它不过x，y的函数:w=w(x,y)。若w已知，则位移，内力、应力均可按上述相应公式求出。在经典分析法中，W(x,y)常设为三角级数形式。例如，四边简支矩形板的W(x,y)设为:(纳维尔解)wx,yAmnsinmxnyasinm1n1b式中Amn为待定系数。假定荷载qx,yqsinmxsinnym1n1mnab4《有限元》讲义则可得位移函数:wx,y1m2qmn2sinmxsinnyD4n2aba2b24.2有限元分析方法一、矩形单元的典型形式将图示矩形薄板沿x,y方向划分成若干小矩形(常取均分)从中拿出一小矩形(单元)，共有四个结点，此时不可以象在平面问题中相同，将结点视为“铰”，而是“刚性的”，即每个结点有三个位移重量：挠度w,绕x、y轴转角挠度w绕x轴转角x上节为沿y方向倾角绕y轴转角y上节为沿x方向倾角即结点i的位移wiwdixiwi1,4yyiwxi同理,相应的结点力fi竖向力FiMxi绕x轴力偶（上节中的My）Myi绕y轴力偶（上节中的Mx）符号重新定义是为了有限元表示的方便，由此得单元结位移向量deddTwx1y1wx4T1414y4节点力FeF1F4Tf1Mx1My1f4Mx4My4T二、位移模式（函数）１、位移模式的采纳插值多项式取为：5《有限元》讲义wx,y12x3y4x25xy6y27x38x2y9xy210y311x3y12xy3(4-2-1)在上式中,前10项取到了三次项的所有432234,最后两项则是从五个四次项xxyxyxyy中选用了两个。没选x2y2是因为它没有多一项与其配对，没选x4,y4它们在界限上结出的挠度函数是四次的，比x3y和xy3要高一次,较之更难满足界限的协调停条件。２、位移模式的检验（三个基本要求：刚体位移，常应变，尽可能的界限协调）①前三项含单元的刚体位移状态：第一项1与坐标x,y没关,表示z方向的挠度是─常量,刚体挪动第二项ywy2x表示刚体转动w第三项xyx3②二次项代表均匀变形状态:2w24，2w22w曲率2y26，25xxy③能保证相邻单元在公共界限上挠度的连续性。④不可以保证相邻单元在公共界限上法线转角的连续性。以单元１～２界限为例，在此界限上yb=常量，代入位移模式4-2-1，可知界限上的挠度W是x的三次函数，合并整理后可得：w12c1c2xc3x2c4x3两个端点共有4个界限条件,(结点1,2的挠度W1,W2,和转角y1,y2。利用他们可唯一确立四个常数C1～C4。因为相邻单元在结点1,2的W,θy对应相同，则两个单元依照四个条件获取的C1～C4亦相同，即两单元在界限拥有同一挠度函数W。⑤法线转角仍以1-2界限为例，将y=-b代入后，此时xdd2xd3x2d4x31但对θx来讲，1,2结点只好供给2个已知条件，不可以完整确立上式，故界限的法线转6《有限元》讲义角不可以保证连续性。所以，这类单元是非协调元，但可以考据这类非协调远是能经过分片试验的。（即当单元划分不停减小时，计算结果还能收敛于精确解。）三、形函数和形函数矩阵。分别将单元结点1,2,3,4的坐标值代入(4-2-1)，并早先求出w，yxy即可获取各结点的位移值。一共可得12个关于i的方程组，联立求解可得：{w}Nde(4-2-2)

w，xN1Nx1Ny1N4Nx4Ny4形函数矩阵:NN1/yNy4/yN1/xNy4/x式中形函数:Ni11i1i2ii228Nxi1bi1i1i128Nyi1i1112(4-2-3)(i=1234)aii8在上边的推导中,我们依旧采纳了局部坐标(无因次坐标)。局部坐标与整体坐标的关系为:1a

xx01yy06z四、单元的几何矩阵[B]和内力矩阵[S]1．几何矩阵[B]7《有限元》讲义由前可知z,将（4－2－3）代入(4-2-4)获取几何矩阵：2N1x22Nx1x22Ny1x22Ny4x22N1y22Ny4（4-2-5）By22N122Ny42xyxy或以子块形式表示：[B]=[BBBB]。1234式中：内力矩阵[S]由基本方程（4-2-5）可获取：FDBdeSde（4-2-6）8《有限元》讲义S称为内力矩阵，把单元的四个结点坐标分别代入4-2-4，求得B后，即可获取S，各节点内力矩阵S的显式：五、单元刚度矩阵abTDBdxdy。将几何矩阵[B]和弹性矩阵[D]的表达式由一般公式得:KtBab代入，积分可得薄板曲折问题矩形单元的单元刚度矩阵的显示：9《有限元》讲义10《有限元》讲义六、荷载等效变换由荷载等效变换的一般公式可得Rab[N(x,y)]Tqx,ydxdyab1．法向均布荷载q代入上述公式得：1N1b3Nx1aNy131N2bNx23aNy2abqab3{R}qdxdyabN31bNx33Ny3aN43Nx41bNy43a32．单元中心点受法向集中力P代入上述公式可得：

2ba2bPaP28ba2ba七、位移界限条件对称、固定边和简支边上支点的已知位移条件以下：对称轴:法线转角=0固定边:挠度=0(或已知值)边线转角=0(或已知值)法线转角=0(或已知值)简支边:挠度=0(或已知值)边线转角=0(或已知值)自由边上节点的挠度、边线和法线转角均为特定参数，同内部节点相同。与11《有限元》讲义板铰接的固定立柱，其节点挠度Ｗ=0，也可以是已知值。八、计算例题例题1:计算图示四边固定方板方板的边长为l，厚度为t，弹性模型量为E，波松比μ=0.3，全板承受均布法向荷载ｑ，求薄板中的挠度和内力。单元划分：为了说明解题方法，采纳最简单的网络2×2，即把方板分成四个矩形单元。因为对称性，只需计算一个单元，比方，计算图中有暗影的单元，单元的节点编号为１，２，３，４。此时，单元的a,b是abl4计算节点荷载:由前面的均布荷载计算公式得：{R}ql2ll12ll12ll12lT[12l]192界限条件：界限23和34为固定边，所以节点2,3,4的挠度、边线和法线转角均为零。边界12和14为对称轴，所以θx1=0、θy1=0。于是，在4个节点和12个位移重量中,只有一个待求的未知量w1。结构的代数方程组：这是一个单元的计算题目，单元刚度矩阵在此处即为总刚度矩阵。引入支承条件后，在总刚度矩阵中只取第一行、列元素，在方程组右端项中只保留第一个元素。于是结构的代数方程为：8D0k1w18D0(816)w1ql215l215l216同此解出w10.00148ql4。此中D0Et32)0.0915Et83D012(1内力:利用式(4-2-6)可求得方板中点力矩为：12《有限元》讲义由表看出，网格越密，计算结果越凑近于精确答案。还可看出，位移的精度一般比内力的精度高，这是因为在位移法中，位移是由基本方程直接求出的，而内力则是依据位移间接求出的。13《有限元》讲义4.3薄板有限元程序设计一、总框图依据曲折板有限元分析方法的解题过程，可写出其总框图以下：┌───────┐│输入原始数据││orCAI│└───┳───┘┌──────┐↓┌──┤算等效结点力│┌───┻───┐│└──────┘│形成荷载列阵├←┘┌─────┐│├←───┤形成单元│└───┳───┘┌─┤定位向量├─┐↓│└─────┘│┌───┻───┐│││形成总刚├←─┘┌─────┐││├←───┤单刚││└───┳───┘└─────┘│↓│┌───┻────┐││解方程输出位移││└───┳────┘│↓┌──────┐│┌───┻───┐│几何矩阵[B]│││├←───┤弹性矩阵[D]│││计算单元内力等│└──────┘││├←───────────┘└───┳───┘↓┌──┻──┐│结束│└─────┘下边结合程序对框图中的内容加以说明。二、子框图１、单元坐标结点编号及单刚形式。14《有限元》讲义为了取挠度向下为正，又能与前述坐标系通通一，特将前述坐标前翻180°(如图)为了能适用板的弹型性分析,程序采纳了应力元和曲折元的组合形式，即每个结点考虑5个位移重量:U,V,W,θx,θy,前2个为平面应力问题的未知量,后3个为曲折板的结点未知量。当只作弹性分析时，平面应力元和曲折元是非藕连的，即单刚的两个副块垣为0,单刚的形式为：u1v1u4v4w1θx1θy4w4θx4θy4┌┐│平面应力元│0│[K]e=│（8×8）││├──────┼──────────────┤││曲折元││0│（12×12）│└┘程序中单刚数组为DK(20,20),子程序:SubroutineDG(A,B,E,T,U)为其形成单刚的子程序。２、自动形成单元编号信息（单元信息数组：[IB]）。３、结点定位向量。４、形成荷载列向量。（a.结点力；b.非结点力(只考虑均布力)）５、总刚，SubroutineZG(M,N,LD,A,B,E,T,U)６、解方程。FJZG( ),HUD( )７、算单元力。SubroutineDYL( )８、算等效结点力。９、弹性矩阵[D]。１０、几何矩阵[B]。三、输入数听闻明１、总信息。共11个(见程序)２、结点拘束信息数组[JB]JB(I,1)──i结点的结点号JB(I,2)──i结点的拘束重量号(1～5)结点拘束信息应依据支承条件或对称条件决定，如算例中所给出的四边简支方板，承受满布均布力，此时可只取板的1/4作为分析对象，以以下图只取右上角1/4板，采纳6×6网络，则每个单元的边长为1米(A=0.5,B=0.5)。设结点编号如图示:在y=0的界限上(1-7结点):挠度w=0(第3个重量)绕y轴转角θy=0(第5个重量)同理,在平形于y轴的x=6m界限