《《有限元基础教程》》由会员分享，可在线阅读，更多相关《《有限元基础教程》（3页珍藏版）》请在金锄头文库上搜索。
1、基于4节点四边形单元的矩形薄板分析(Quad2D4Node)如 图 4-21 所 示的 一个 薄 平板 ，在 右 端部 受 集中 力 F 作 用 ，其 中 的参 数为 ：E二lx 107Pa,卩=1/3, t=O.lm, F=1 x 105N 。基于MATLAB平台，按平面应力问题计算各个节点位移、支座反力以及单元的应力。F(a) 问题描述(b) 有限元分析模型图 4-21 右端部受集中力作用的薄平板解答：对该问题进行有限元分析的过程如下。(1) 结构的离散化与编号将结构离散为二个4节点矩形单元，单元编号及节点编号如图4-21(b)所示，连接关系 见表4-6，节点的几何坐标见表4-7，载荷F按
2、静力等效原则向节点1，2移置。表 4-6 结构的单元连接关系单元号节点号13564213 4 2表4-7节点的坐标节点节点坐标/mxy121220311410501600节点位移列阵q = u1vu12vuvuvuvuv r(4-186)233445566节点外载列阵F =0 -F 0F00000000T(4-187)22 R=0000000总的节点载荷列阵P = F + R = R =0F0F0 022约束的支反力列阵0 R R R RT(4-188)x5y5x6y 6T0 0 R R R R(4-189)x5y 5x 6y 6其中，(R ,R )和(R ,R )分别为节点5和节点6的两个方
3、向的支反力。x5 y5x 6y 6(2) 计算各单元的刚度矩阵(以国际标准单位)首先在MATLAB环境下，输入弹性模量E、泊松比NU、薄板厚度h和平面应力问题 性质指示参数ID,然后针对单元1和单元2,分别两次调用函数Quad2D4Node_Stiffness, 就可以得到单元的刚度矩阵k1(8X 8)和k2(8X 8)。 E=1e7; NU=1/3; t=0.1; ID=1; k1= Quad2D4Node_Stiffness(E,NU,t, 1,1, 0,1,0,0,1,0, ID); k2=Quad2D4Node_Stiffness(E,NU,t, 2,1,1,1,1,0,2,0, ID
4、);(3) 建立整体刚度方程由于该结构共有6个节点，则总共的自由度数为12，因此，结构总的刚度矩阵为KK(12 X12),先对KK清零，然后两次调用函数Quad2D4Node_Assembly进行刚度矩阵的组装。 KK = zeros(12,12);KK= Quad2D4Node_Assembly(KK,k1, 3,5,6,4);KK= Quad2D4Node_Assembly(KK,k2, 1,3,4,2);(4) 边界条件的处理及刚度方程求解由图4-21(b)可以看出，节点5和节点6的两个方向的位移将为零，即u = 0, v = 0, u = 0, v = 0，因此，将针对节点1, 2，3
5、和4的位移进行求解。节点1, 2, 55663和4的位移将对应KK矩阵中的前8行和前8列，则需从KK(12X12)中提出，置给k,然后 生成对应的载荷列阵P，再采用高斯消去法进行求解。注意：MATLAB中的反斜线符号“” 就是采用高斯消去法。k=KK(1:8,1:8);P=0;-50000;0;-50000;0;0;0;0;u=kPu = 0.8000 -2.5333 -0.8000 -2.53330.6000-0.8667 -0.6000 -0.8667 将列排成行 由此可以看出，所求得的结果为：u = 0.800 0, v = 2.533 3, u = -0.800 0, v =2.533
6、 3,1 1 2 2u = 0.600 0, v = 0.866 7, u = 0.600 0, v = 0.866 7 。3344(5) 支反力的计算在得到整个结构的节点位移后，由原整体刚度方程就可以计算出对应的支反力。先将 上面得到的位移结果与位移边界条件的节点位移进行组合(注意位置关系)，可以得到整体的位移列阵U(12X1),再代回原整体刚度方程，计算出所有的节点力P(12X1),按式(4-189) 的对应关系就可以找到对应的支反力。U=u;0;0;0;0;P=KK*UP = 1.0e+005 *-0.0000-0.50000.0000-0.50000.00000.0000将列排成行0.
7、0000 -0.0000-2.00000.50002.00000.5000将列排成行由 式 (4-189) 的对应关系,可以得到 对 应 的 支 反 力 为R = -200 000, R = 50 000, R = 200 000, R = 50 000。x5y5x6y 6(6) 各单元的应力计算先从整体位移列阵U(12X 1)中提取出单元的位移列阵, 数Quad2D4Node_Stress,就可以得到各个单元的应力分量。 u1= U(5);U(6) ;U(9);U(10);U(11);U(12);U(7);U(8)u1 = 0.6000 -0.8667 0 0 0 0 -0.6000 -0.
8、8667 stress1=Quad2D4Node_Stress(E,NU, 1,1,0,1,0,0,1,0, u1,ID) stress1 = 1.0e+005 *0.0000 0.0000 -10.0000 u2= U(1);U(2) ;U(5);U(6); U(7);U(8);U(3);U(4)u2 = 0.8000-2.53330.6000-0.8667然后，调用计算单元应力的函将列排成行将列排成行将列排成行将列排成行-0.6000-0.8667-0.8000-2.5333 stress2=Quad2D4Node_Stress(E,NU, 2,1,1,1,1,0,2,0, u2,ID) stress2 = 1.0e+005 *0.0000 0.0000 -10.0000可以看出：计算得到的单元1的应力分量为b二0Pa, b x将列排成行二 0Pa, t =-1 000 000Pa ；xy单元 2 的应力分量为b = 0Pa, b = 0Pa, t = 1 000 000Pa。x y xy
