有限元4薄板弯曲问题

-第4章 弹性薄板弯曲问题的有限元法 薄板弯曲问题在理论上和应用上都具有重要意义,并有专门著作加以论述如耀乾平板理论。象其它弹性力学问题一样，用微分方程、差分法等经典方法所能求解的薄板问题很有限，一般只能解决等厚、小孔口、支承情况较简单的单跨板。故工程设计中以往多采用简化、近似、图表等方法来解决板的设计问题。在板的分析中，常取板的中面为*oy平面(如图)。平板构造按其厚度t与短边a的比值大小而分为:厚板Thick plate和薄板Thin plate)两种。当时称为薄板 平板上所承受的荷载通常有两种: 1. 面拉压荷载。 由面拉压刚度承担, 属平面应力问题。 2. 垂直于板的法向荷载, 弯扭变形为主,具有梁的受力特征, 即常说的弯曲问题。平板在垂直于板面的荷载作用下产生挠度W。当最大挠度w远小于t时, 称为小挠度问题(or刚性板)(stiffness plate)当最大挠度w与t相差不大时,称为大挠度问题(or柔性板)(fle*ure plate)(工程定义: 为刚性板；为柔性板； 为绝对柔性板。)4.1 根本理论一、根本假定、略去垂直于中面的法向应力。，即以中面上沿Z方向的挠度W代表板的挠度)、变形前垂直中面的任意直线，变形后仍保持为垂直中面的直线。法向假定,、板弯曲时，中面不产生应力。(中面中性层假定) 上述假定常称为薄板小挠度问题假定or 柯克霍夫假定。符合上述假定的平板即为刚性板。二、根本方法 以上述假定为根底，板分析中常用挠度作为根本未知量,下面介绍以为根本未知量所导出的有关方程。 、几何方程应变挠度关系弹性曲面沿*, y 方向的倾角从中面取出一微小矩形ABCD，如以下图，设其边长为d*, dy，变形后弯曲成曲面ABCD设A点挠度, 则沿*方向倾角(绕y轴) (B点绕度 )沿y方向倾角(绕*轴) (D点绕度 ) 沿*, y 方向位移作平行于平面,设中面上点A到A1的距离为Z，变形后,A点有挠度W, 同时发生弯曲，曲面沿*方向的倾角为, 根据法线假定,则A1点沿*方向的位移: (负号为方向与*相反) 同理取平面得： (4-1-1) Z平面的应变分量和曲、扭率 根本假定，由于, 故板任意点的应变与平面问题一样: (4-1-2) 此为Z平面的应变挠度度几何方程。上式中的,为曲面在*,Y方向的曲、扭率，记为： (4-1-3)所以, 、物理方程应力挠度关系 由于忽略z对变形的影响, 因此z平面的应力应变关系具有与平面问题一样的形式:将(4-1-2)代入得:或简写为： (4-1-4)式中弹性矩阵: 、力方程力挠度关系从板取微元体, 由其上正应力,和剪应力,可在截面上合成合力矩:(面上由产生的绕Y轴弯矩)(面上由产生的绕*轴弯矩)扭矩: 由剪应力产生，如图假定 分别表示单位宽度上的力矩。如是，力矩阵：简写成 (4-1-5)比较(4-1-4)和(4-1-5)可得用力矩表示的平板应力: 由此可见，平板上、下外表处的应力最大： 以上是薄板弯曲问题中的根本公式,从中可见其挠度W是弯曲问题中的根本未知函数。且由于忽略了z方向的变化,因此它只是*，y 的函数: w=w(*, y)。假设w，则位移，力、应力均可按上述相应公式求出。在经典解析法中，W(*, y)常设为三角级数形式。例如，四边简支矩形板的W(*, y)设为: (纳维尔解)式中为待定系数。假定荷载 则可得位移函数: 4.2 有限元分析方法一、矩形单元的典型形式 将图示矩形薄板沿*,y方向划分成假设干小矩形(常取等分)从中取出一小矩形(单元)，共有四个结点，此时不能象在平面问题中一样，将结点视为铰，而是刚性的，即每个结点有三个位移分量：挠度,绕*、y轴转角 即结点i的位移同理,相应的结点力符号重新定义是为了有限元表示的方便，由此得单元结位移向量节点力二、 位移模式函数、位移模式的选取插值多项式取为： (4-2-1)在上式中,前10项取到了三次项的全部,最后两项则是从五个四次项中选用了两个。没选是因为它没有多一项与其配对，没选它们在边界上结出的挠度函数是四次的，比和要高一次,较之更难满足边界的协调和条件。、位移模式的检验三个根本要求： 刚体位移，常应变，尽可能的边界协调 前三项含单元的刚体位移状态：第一项与坐标*, y无关, 表示z方向的挠度是常量, 刚体移动 二次项代表均匀变形状态:曲率 ， ， 能保证相邻单元在公共边界上挠度的连续性。 不能保证相邻单元在公共边界上法线转角的连续性。以单元边界为例，在此边界上=常量，代入位移模式 4-2-1，可知边界上的挠度W是*的三次函数，合并整理后可得：两个端点共有4个边界条件,(结点1,2的挠度W1 , W2 ,和转角。利用他们可唯一确定四个常数C1C4。因为相邻单元在结点1, 2的W, y对应一样，则两个单元依据四个条件得到的C1C4 亦一样，即两单元在边界具有同一挠度函数W。 法线转角仍以1-2边界为例，将y=-b代入后，此时 但对*来讲，1, 2结点只能提供2个条件，不能完全确定上式，故边界的法线转角不能保证连续性。 因此，这种单元是非协调元，但可以验证这种非协调远是能通过分片试验的。即当单元划分不断缩小时，计算结果仍能收敛于准确解。三、形函数和形函数矩阵。 分别将单元结点1, 2, 3, 4的坐标值代入(4-2-1)，并事先求出，便可得到各结点的位移值。一共可得12个关于的方程组，联立求解可得：(4-2-2) 形函数矩阵: 式中形函数:(4-2-3)(i=1 2 3 4) 在上面的推导中,我们仍然选用了局部坐标(无因次坐标)。局部坐标与整体坐标的关系为:四、单元的几何矩阵B和力矩阵S1几何矩阵B由前可知 , 将423代入(4-2-4)得到几何矩阵： 4-2-5或以子块形式表示： B=B1 B2 B3 B4。式中：2.力矩阵S由根本方程4-2-5可得到： 4-2-6称为力矩阵，把单元的四个结点坐标分别代入4-2-4，求得后，即可获得，各节点力矩阵的显式：五、单元刚度矩阵由一般公式得: 。将几何矩阵B和弹性矩阵D的表达式代入，积分可得薄板弯曲问题矩形单元的单元刚度矩阵的显示：六、荷载等效变换由荷载等效变换的一般公式可得1法向均布荷载q代入上述公式得：2单元中心点受法向集中力P代入上述公式可得：七、位移边界条件 对称、固定边和简支边上支点的位移条件如下：对称轴: 法线转角=0固定边: 挠度=0 (或值) 边线转角=0 (或值) 法线转角=0 (或值)简支边: 挠度=0 (或值) 边线转角=0 (或值)自由边上节点的挠度、边线和法线转角均为特定参数，同部节点一样。与板铰接的固定立柱，其节点挠度 = 0，也可以是值。八、计算例题例题1: 计算图示四边固定方板 方板的边长为l，厚度为t，弹性模型量为E，波松比=0.3，全板承受均布法向荷载，求薄板中的挠度和力。 单元划分： 为了说明解题方法，采用最简单的网络22，即把方板分成四个矩形单元。由于对称性，只需计算一个单元，例如，计算图中有阴影的单元，单元的节点编号为，。 此时，单元的a, b是计算节点荷载: 由前面的均布荷载计算公式得：边界条件：边界23和34为固定边，因此节点2, 3, 4的挠度、边线和法线转角均为零。边界12和14为对称轴，因此*1=0、y1=0。于是，在4个节点和12个位移分量中,只有一个待求的未知量。构造的代数方程组：这是一个单元的计算题目，单元刚度矩阵在此处即为总刚度矩阵。引入支承条件后，在总刚度矩阵中只取第一行、列元素，在方程组右端项中只保存第一个元素。于是构造的代数方程为：同此解出 。其中 力: 利用式(4-2-6)可求得方板中点力矩为： 由表看出，网格越密，计算结果越接近于准确答案。还可看出，位移的精度一般比力的精度高，这是因为在位移法中，位移是由根本方程直接求出的，而力则是根据位移间接求出的。4.3 薄板有限元程序设计一、总框图 根据弯曲板有限元分析方法的解题过程，可写出其总框图如下： 输入原始数据 or CAI 算等效结点力 形成荷载列阵 形成单元 定位向量 形 成 总 刚 单 刚 解方程输出位移 几何矩阵B 弹性矩阵D 计算单元力等 结 束 下面结合程序对框图中的容加以说明。二、子框图 、单元坐标结点编号及单刚形式。 为了取挠度向下为正，又能与前述坐标系统统一，特将前述坐标前翻180(如图)为了能适用板的弹型性分析,程序采用了应力元和弯曲元的组合形式，即每个结点考虑5个位移分量: U, V, W, * , y , 前2个为平面应力问题的未知量,后3个为弯曲板的结点未知量。当只作弹性分析时，平面应力元和弯曲元是非藕连的，即单刚的两个副块垣为0, 单刚的形式为： u1 v1 u4 v4 w1 *1 y4 w4 *4 y4 平面应力元 0 Ke = 88 弯曲元 0 1212 程序中单刚数组为 DK(20, 20), 子程序:Subroutine DG(A, B, E, T, U)为其形成单刚的子程序。 、自动形成单元编号信息 单元信息数组：IB。 、结点定位向量。 、形成荷载列向量。 a. 结点力； b. 非结点力(只考虑均布力) 、总刚，Subroutine ZG(M, N, LD, A, B, E, T,U) 、解方程。 FJZG( ), HUD( ) 、算单元力。 Subroutine DYL( ) 、算等效结点力。 、弹性矩阵D 。 、几何矩阵B。 三、输入数据说明 、总信息。共11个(见程序)、结点约束信息数组 (I,1) i结点的结点号 (I,2) i结点的约束分量号(15)结点约束信息应根据支承条件或对称条件决定，如算例中所给出的四边简支方板，承受满布均布力，此时可只取板的1/4作为分析对象，如以下图只取右上角1/4板，采用66网络，则每个单元的边长为1米(A=0.5, B=0.5)。设结点编号如图示: 在y=0的边界上(1-7结点): 挠度 w=0 (第3个分量) 绕y轴转角y =0 (第5个分量) 同理,在平形于y轴的*=6m边界上: w=0, *=0 (3,4分量) 在对称轴*=0 边界上 u=0, y =0 在对称轴y=6 边界上 v=0, * =0 中点(43点)除W外, 其余均=0 同时，在简支边上,也可设u和v均为 0 这样便共有77个约束。另外，也可通过改变约束信息来改变计算简图，如同样网格数的1/2、1/4板，固支边界板等。、结点荷载信息 pp(I, 1) 荷载值 pp(I, 2) 结点号位移分量号(如竖向力的位移分量为3) 、非结点荷载信息 程序中均为考虑了满布均布力 PF(I, 1) 荷载值PF(I, 2) 荷载作用单元号四、输出信息1.结点位移；2.单元力王勖成P144,5章2节应力计算结果的处理与改善。. z.